{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"},{"sourceId":5499219,"sourceType":"datasetVersion","datasetId":3167603},{"sourceId":5549164,"sourceType":"datasetVersion","datasetId":3197305},{"sourceId":5607816,"sourceType":"datasetVersion","datasetId":3225525},{"sourceId":5792099,"sourceType":"datasetVersion","datasetId":3327296},{"sourceId":6247561,"sourceType":"datasetVersion","datasetId":3590060}],"dockerImageVersionId":31154,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-12-01T09:25:03.489634Z","iopub.execute_input":"2025-12-01T09:25:03.489910Z","iopub.status.idle":"2025-12-01T09:25:05.256176Z","shell.execute_reply.started":"2025-12-01T09:25:03.489882Z","shell.execute_reply":"2025-12-01T09:25:05.255426Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# CAFA 6 PROTEIN FUNCTION PREDICTION - ENHANCED VERSION (FIXED)\n# ============================================================================\n\n# ⚙️ CONFIGURATION - ADJUST THIS FOR SPEED VS ACCURACY\nSAMPLE_PERCENT = 100  # Use 100% of data\nQUICK_MODE = True   # Enable full feature computation\n\n# Package Installation\nimport subprocess\nimport sys\n\ndef install(package):\n    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", package])\n\nprint(\"Installing required packages...\")\ntry:\n    import obonet\nexcept:\n    install('obonet')\n    import obonet\n\ntry:\n    from Bio import SeqIO\nexcept:\n    install('biopython')\n    from Bio import SeqIO\n\n# Core Imports\nimport pandas as pd\nimport numpy as np\nfrom pathlib import Path\nfrom collections import defaultdict, Counter\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport networkx as nx\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.multiclass import OneVsRestClassifier\nfrom sklearn.linear_model import LogisticRegression\n\n# Visualization imports\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom matplotlib.patches import Rectangle\nimport matplotlib.patches as mpatches\nplt.style.use('seaborn-v0_8-whitegrid')\nsns.set_palette(\"husl\")\n\nprint(\"=\"*80)\nprint(\"CAFA 6 PROTEIN FUNCTION PREDICTION - ENHANCED STARTER\")\nprint(f\"📊 SAMPLE MODE: {SAMPLE_PERCENT}% of data\")\nprint(f\"⚡ QUICK MODE: {'ON' if QUICK_MODE else 'OFF'}\")\nprint(\"=\"*80)\n\n# ============================================================================\n# 1. DEFINE PATHS\n# ============================================================================\nBASE = Path('/kaggle/input/cafa-6-protein-function-prediction')\nTRAIN_DIR = BASE / 'Train'\nTEST_DIR = BASE / 'Test'\n\n# ============================================================================\n# 2. LOAD GO ONTOLOGY (WITH HIERARCHY ANALYSIS)\n# ============================================================================\nprint(\"\\n[1/9] Loading GO ontology...\")\ngo_graph = obonet.read_obo(TRAIN_DIR / 'go-basic.obo')\nprint(f\"   ✓ Loaded {len(go_graph)} GO terms\")\n\n# Map terms to ontologies\nterm_to_ont = {}\nterm_names = {}\nfor term_id in go_graph.nodes():\n    if 'namespace' in go_graph.nodes[term_id]:\n        ns = go_graph.nodes[term_id]['namespace']\n        if ns == 'biological_process':\n            term_to_ont[term_id] = 'BPO'\n        elif ns == 'cellular_component':\n            term_to_ont[term_id] = 'CCO'\n        elif ns == 'molecular_function':\n            term_to_ont[term_id] = 'MFO'\n    if 'name' in go_graph.nodes[term_id]:\n        term_names[term_id] = go_graph.nodes[term_id]['name']\n\nont_counts = pd.Series(term_to_ont).value_counts()\nprint(f\"   ✓ Ontology breakdown: MF={ont_counts.get('MFO',0)}, BP={ont_counts.get('BPO',0)}, CC={ont_counts.get('CCO',0)}\")\n\n# Analyze GO hierarchy depth (sample for speed)\ndef get_term_depth(graph, term_id):\n    \"\"\"Calculate depth of term in GO hierarchy\"\"\"\n    try:\n        paths = []\n        for root in ['GO:0008150', 'GO:0005575', 'GO:0003674']:\n            if nx.has_path(graph, term_id, root):\n                paths.append(nx.shortest_path_length(graph, term_id, root))\n        return max(paths) if paths else 0\n    except:\n        return 0\n\nprint(\"   Computing GO hierarchy depths...\")\nsample_terms_for_depth = list(term_to_ont.keys())[:1000]\nterm_depths = {term: get_term_depth(go_graph, term) for term in sample_terms_for_depth}\n\n# Visualize ontology with enhanced graphics\nfig = plt.figure(figsize=(18, 10))\ngs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)\n\n# Main ontology distribution\nax1 = fig.add_subplot(gs[0, :2])\ncolors = ['#FF6B6B', '#4ECDC4', '#45B7D1']\nbars = ax1.bar(range(len(ont_counts)), ont_counts.values, color=colors, \n               edgecolor='black', linewidth=2, alpha=0.8)\nax1.set_xticks(range(len(ont_counts)))\nax1.set_xticklabels(['Molecular Function', 'Biological Process', 'Cellular Component'], \n                     rotation=0, fontsize=11, fontweight='bold')\nax1.set_title('GO Term Distribution by Ontology', fontsize=14, fontweight='bold', pad=20)\nax1.set_ylabel('Number of Terms', fontsize=12, fontweight='bold')\nax1.grid(axis='y', alpha=0.3)\nfor i, (v, bar) in enumerate(zip(ont_counts.values, bars)):\n    height = bar.get_height()\n    ax1.text(bar.get_x() + bar.get_width()/2., height,\n             f'{v:,}\\n({v/ont_counts.sum()*100:.1f}%)',\n             ha='center', va='bottom', fontweight='bold', fontsize=11)\n\n# Hierarchy depth distribution\nax2 = fig.add_subplot(gs[0, 2])\ndepth_values = list(term_depths.values())\nax2.hist(depth_values, bins=20, color='#A8E6CF', edgecolor='black', alpha=0.7)\nax2.set_title('GO Term Depth\\nDistribution', fontsize=11, fontweight='bold')\nax2.set_xlabel('Hierarchy Depth', fontsize=10)\nax2.set_ylabel('Count', fontsize=10)\nax2.axvline(np.mean(depth_values), color='red', linestyle='--', linewidth=2,\n            label=f'Mean: {np.mean(depth_values):.1f}')\nax2.legend(fontsize=9)\n\n# Network visualization (sample of GO graph)\nax3 = fig.add_subplot(gs[1:, :])\nsample_terms = list(term_to_ont.keys())[:50]\nsubgraph = go_graph.subgraph(sample_terms)\npos = nx.spring_layout(subgraph, k=0.5, iterations=50, seed=42)\nnode_colors = [colors[['MFO', 'BPO', 'CCO'].index(term_to_ont.get(node, 'MFO'))] \n               for node in subgraph.nodes()]\nnx.draw_networkx_nodes(subgraph, pos, node_color=node_colors, \n                       node_size=300, alpha=0.7, ax=ax3)\nnx.draw_networkx_edges(subgraph, pos, alpha=0.2, arrows=True, \n                       arrowsize=10, ax=ax3, edge_color='gray')\nax3.set_title('GO Ontology Network Structure (Sample of 50 terms)', \n              fontsize=13, fontweight='bold', pad=15)\nax3.axis('off')\n\n# Legend\nlegend_elements = [mpatches.Patch(facecolor=colors[0], label='Molecular Function'),\n                   mpatches.Patch(facecolor=colors[1], label='Biological Process'),\n                   mpatches.Patch(facecolor=colors[2], label='Cellular Component')]\nax3.legend(handles=legend_elements, loc='upper right', fontsize=10, framealpha=0.9)\n\nplt.suptitle('Gene Ontology Analysis', fontsize=16, fontweight='bold', y=0.995)\nplt.tight_layout()\nplt.show()\n\n# ============================================================================\n# 3. LOAD IA WEIGHTS (WITH ANALYSIS)\n# ============================================================================\nprint(\"\\n[2/9] Loading IA weights...\")\nia_df = pd.read_csv(BASE / 'IA.tsv', sep='\\t', header=None, names=['term', 'ia'])\n\nif SAMPLE_PERCENT < 100:\n    ia_df = ia_df.sample(frac=SAMPLE_PERCENT/100, random_state=42)\n\nia_dict = dict(zip(ia_df['term'], ia_df['ia']))\nprint(f\"   ✓ Loaded {len(ia_dict)} IA weights\")\n\n# Enhanced IA visualization\nfig, axes = plt.subplots(2, 3, figsize=(18, 10))\n\n# IA distribution by ontology\nia_by_ont = ia_df.copy()\nia_by_ont['ontology'] = ia_by_ont['term'].map(term_to_ont)\nia_by_ont = ia_by_ont.dropna()\n\naxes[0, 0].hist(ia_df['ia'], bins=50, color='#95E1D3', edgecolor='black', alpha=0.7)\naxes[0, 0].set_title('Overall IA Distribution', fontsize=12, fontweight='bold')\naxes[0, 0].set_xlabel('IA Weight', fontsize=10)\naxes[0, 0].set_ylabel('Frequency', fontsize=10)\naxes[0, 0].axvline(ia_df['ia'].mean(), color='red', linestyle='--', linewidth=2, \n                   label=f'Mean: {ia_df[\"ia\"].mean():.2f}')\naxes[0, 0].legend()\n\n# Box plot by ontology\nont_data = [ia_by_ont[ia_by_ont['ontology']==ont]['ia'].values \n            for ont in ['MFO', 'BPO', 'CCO']]\nbp = axes[0, 1].boxplot(ont_data, labels=['MF', 'BP', 'CC'], patch_artist=True)\nfor patch, color in zip(bp['boxes'], colors):\n    patch.set_facecolor(color)\n    patch.set_alpha(0.7)\naxes[0, 1].set_title('IA Weights by Ontology', fontsize=12, fontweight='bold')\naxes[0, 1].set_ylabel('IA Weight', fontsize=10)\naxes[0, 1].grid(axis='y', alpha=0.3)\n\n# Violin plot\nparts = axes[0, 2].violinplot(ont_data, positions=[1, 2, 3], showmeans=True, showmedians=True)\nfor pc, color in zip(parts['bodies'], colors):\n    pc.set_facecolor(color)\n    pc.set_alpha(0.7)\naxes[0, 2].set_xticks([1, 2, 3])\naxes[0, 2].set_xticklabels(['MF', 'BP', 'CC'])\naxes[0, 2].set_title('IA Distribution Density', fontsize=12, fontweight='bold')\naxes[0, 2].set_ylabel('IA Weight', fontsize=10)\n\n# Cumulative distribution\nsorted_ia = np.sort(ia_df['ia'].values)\ncumsum = np.cumsum(sorted_ia) / np.sum(sorted_ia)\naxes[1, 0].plot(sorted_ia, cumsum, linewidth=2, color='#6C5CE7')\naxes[1, 0].set_title('Cumulative IA Distribution', fontsize=12, fontweight='bold')\naxes[1, 0].set_xlabel('IA Weight', fontsize=10)\naxes[1, 0].set_ylabel('Cumulative Proportion', fontsize=10)\naxes[1, 0].grid(alpha=0.3)\naxes[1, 0].axhline(0.5, color='red', linestyle='--', alpha=0.5, label='50%')\naxes[1, 0].legend()\n\n# Top terms by IA\ntop_ia = ia_df.nlargest(15, 'ia')\naxes[1, 1].barh(range(len(top_ia)), top_ia['ia'].values, color='#FF7675', edgecolor='black')\naxes[1, 1].set_yticks(range(len(top_ia)))\naxes[1, 1].set_yticklabels([f\"{t[:15]}...\" if len(t) > 15 else t \n                            for t in top_ia['term'].values], fontsize=8)\naxes[1, 1].set_title('Top 15 Terms by IA Weight', fontsize=12, fontweight='bold')\naxes[1, 1].set_xlabel('IA Weight', fontsize=10)\naxes[1, 1].invert_yaxis()\n\n# Statistics summary\naxes[1, 2].axis('off')\nia_stats = f\"\"\"\nIA WEIGHT STATISTICS\n\nTotal terms: {len(ia_df):,}\n\nOverall:\n  • Mean: {ia_df['ia'].mean():.3f}\n  • Median: {ia_df['ia'].median():.3f}\n  • Std Dev: {ia_df['ia'].std():.3f}\n  • Range: [{ia_df['ia'].min():.3f}, {ia_df['ia'].max():.3f}]\n\nBy Ontology (Mean ± Std):\n  • MF: {ia_by_ont[ia_by_ont['ontology']=='MFO']['ia'].mean():.3f} ± {ia_by_ont[ia_by_ont['ontology']=='MFO']['ia'].std():.3f}\n  • BP: {ia_by_ont[ia_by_ont['ontology']=='BPO']['ia'].mean():.3f} ± {ia_by_ont[ia_by_ont['ontology']=='BPO']['ia'].std():.3f}\n  • CC: {ia_by_ont[ia_by_ont['ontology']=='CCO']['ia'].mean():.3f} ± {ia_by_ont[ia_by_ont['ontology']=='CCO']['ia'].std():.3f}\n\"\"\"\naxes[1, 2].text(0.05, 0.5, ia_stats, fontsize=10, family='monospace',\n                verticalalignment='center')\n\nplt.suptitle('Information Accretion (IA) Analysis', fontsize=14, fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n# ============================================================================\n# 4. LOAD TRAINING DATA (WITH COMPREHENSIVE ANALYSIS) - FIXED\n# ============================================================================\nprint(\"\\n[3/9] Loading training data...\")\n\ntrain_terms = pd.read_csv(TRAIN_DIR / 'train_terms.tsv', sep='\\t', \n                          names=['protein', 'term', 'ontology'])\ntrain_taxonomy = pd.read_csv(TRAIN_DIR / 'train_taxonomy.tsv', sep='\\t',\n                             names=['protein', 'taxon'])\n\nprint(f\"   ✓ Full dataset: {len(train_terms)} annotations, {train_terms['protein'].nunique()} proteins\")\n\n# SAMPLE proteins for faster iteration\nif SAMPLE_PERCENT < 100:\n    sample_proteins = train_terms['protein'].drop_duplicates().sample(\n        frac=SAMPLE_PERCENT/100, random_state=42\n    ).tolist()\n    train_terms = train_terms[train_terms['protein'].isin(sample_proteins)]\n    train_taxonomy = train_taxonomy[train_taxonomy['protein'].isin(sample_proteins)]\n    print(f\"   ✓ Sampled to {SAMPLE_PERCENT}%: {len(train_terms)} annotations, {len(sample_proteins)} proteins\")\n\n# Print ontology distribution\nprint(f\"\\n   Ontology distribution:\")\nprint(train_terms['ontology'].value_counts())\n\n# Comprehensive training data visualization - FIXED\nfig = plt.figure(figsize=(20, 12))\ngs = fig.add_gridspec(3, 4, hspace=0.35, wspace=0.35)\n\n# 1. Ontology distribution - FIXED to handle all possible ontology codes\nax1 = fig.add_subplot(gs[0, 0])\nont_dist = train_terms['ontology'].value_counts()\n\n# Map ontology codes (handle F, P, C or any other codes)\ncolors_ont_map = {'F': '#FF6B6B', 'P': '#4ECDC4', 'C': '#45B7D1'}\nont_names_map = {'F': 'MF', 'P': 'BP', 'C': 'CC'}\n\n# Get colors and names, with defaults for unknown codes\ncolors_list = [colors_ont_map.get(k, '#CCCCCC') for k in ont_dist.index]\nlabels_list = [ont_names_map.get(k, k) for k in ont_dist.index]\n\nbars = ax1.bar(range(len(ont_dist)), ont_dist.values, color=colors_list, \n               edgecolor='black', linewidth=1.5)\nax1.set_xticks(range(len(ont_dist)))\nax1.set_xticklabels(labels_list)\nax1.set_title('Annotations by Ontology', fontsize=11, fontweight='bold')\nax1.set_ylabel('Count', fontsize=9)\nfor i, (v, bar) in enumerate(zip(ont_dist.values, bars)):\n    ax1.text(bar.get_x() + bar.get_width()/2., v, f'{v:,}', \n             ha='center', va='bottom', fontweight='bold', fontsize=9)\n\n# 2. Top terms\nax2 = fig.add_subplot(gs[0, 1:3])\ntop_terms = train_terms['term'].value_counts().head(20)\nax2.barh(range(len(top_terms)), top_terms.values, color='#A8E6CF', edgecolor='black')\nax2.set_yticks(range(len(top_terms)))\nax2.set_yticklabels([f\"{term_names.get(t, t)[:30]}...\" if len(term_names.get(t, t)) > 30 \n                     else term_names.get(t, t) for t in top_terms.index], fontsize=8)\nax2.set_title('Top 20 Most Frequent GO Terms', fontsize=11, fontweight='bold')\nax2.set_xlabel('Count', fontsize=9)\nax2.invert_yaxis()\n\n# 3. Terms per protein\nax3 = fig.add_subplot(gs[0, 3])\nterms_per_protein = train_terms.groupby('protein').size()\nax3.hist(terms_per_protein, bins=50, color='#FFD93D', edgecolor='black', alpha=0.7)\nax3.set_title('Terms per Protein', fontsize=11, fontweight='bold')\nax3.set_xlabel('# Terms', fontsize=9)\nax3.set_ylabel('Frequency', fontsize=9)\nax3.axvline(terms_per_protein.mean(), color='red', linestyle='--', linewidth=2)\n\n# 4. Proteins per term\nax4 = fig.add_subplot(gs[1, 0])\nproteins_per_term = train_terms.groupby('term').size()\nax4.hist(proteins_per_term, bins=50, color='#FFEAA7', edgecolor='black', alpha=0.7, log=True)\nax4.set_title('Proteins per Term (log)', fontsize=11, fontweight='bold')\nax4.set_xlabel('# Proteins', fontsize=9)\nax4.set_ylabel('# Terms (log)', fontsize=9)\n\n# 5. Taxonomy distribution\nax5 = fig.add_subplot(gs[1, 1])\ntop_taxa = train_taxonomy['taxon'].value_counts().head(10)\nax5.bar(range(len(top_taxa)), top_taxa.values, color='#74B9FF', edgecolor='black')\nax5.set_xticks(range(len(top_taxa)))\nax5.set_xticklabels([str(t)[:8] for t in top_taxa.index], rotation=45, ha='right', fontsize=8)\nax5.set_title('Top 10 Species', fontsize=11, fontweight='bold')\nax5.set_ylabel('# Proteins', fontsize=9)\n\n# 6. Term co-occurrence heatmap\nax6 = fig.add_subplot(gs[1, 2:])\ntop_10_terms = train_terms['term'].value_counts().head(10).index\ncooc_matrix = np.zeros((10, 10))\nfor i, t1 in enumerate(top_10_terms):\n    for j, t2 in enumerate(top_10_terms):\n        if i != j:\n            proteins_t1 = set(train_terms[train_terms['term']==t1]['protein'])\n            proteins_t2 = set(train_terms[train_terms['term']==t2]['protein'])\n            cooc_matrix[i,j] = len(proteins_t1 & proteins_t2)\nim = ax6.imshow(cooc_matrix, cmap='YlOrRd', aspect='auto')\nax6.set_xticks(range(10))\nax6.set_yticks(range(10))\nax6.set_xticklabels([term_names.get(t, t)[:10] for t in top_10_terms], \n                     rotation=45, ha='right', fontsize=7)\nax6.set_yticklabels([term_names.get(t, t)[:10] for t in top_10_terms], fontsize=7)\nax6.set_title('Term Co-occurrence Matrix', fontsize=11, fontweight='bold')\nplt.colorbar(im, ax=ax6, label='# Shared Proteins')\n\n# 7. Annotation density\nax7 = fig.add_subplot(gs[2, :2])\nterm_freq_bins = pd.cut(proteins_per_term, bins=[0, 10, 50, 100, 500, 100000], \n                        labels=['<10', '10-50', '50-100', '100-500', '>500'])\nfreq_dist = term_freq_bins.value_counts().sort_index()\nax7.bar(range(len(freq_dist)), freq_dist.values, color='#E17055', edgecolor='black', alpha=0.7)\nax7.set_xticks(range(len(freq_dist)))\nax7.set_xticklabels(freq_dist.index, rotation=0)\nax7.set_title('GO Term Frequency Distribution', fontsize=11, fontweight='bold')\nax7.set_xlabel('# Proteins with Term', fontsize=9)\nax7.set_ylabel('# Terms', fontsize=9)\nfor i, v in enumerate(freq_dist.values):\n    ax7.text(i, v, str(v), ha='center', va='bottom', fontweight='bold')\n\n# 8. Summary statistics\nax8 = fig.add_subplot(gs[2, 2:])\nax8.axis('off')\nsummary_text = f\"\"\"\nTRAINING DATA COMPREHENSIVE SUMMARY\n\nDataset Size:\n  • Total Annotations: {len(train_terms):,}\n  • Unique Proteins: {train_terms['protein'].nunique():,}\n  • Unique GO Terms: {train_terms['term'].nunique():,}\n  • Species: {train_taxonomy['taxon'].nunique()}\n\nOntology Distribution:\n  • Molecular Function: {ont_dist.get('F', 0):,} ({ont_dist.get('F', 0)/len(train_terms)*100:.1f}%)\n  • Biological Process: {ont_dist.get('P', 0):,} ({ont_dist.get('P', 0)/len(train_terms)*100:.1f}%)\n  • Cellular Component: {ont_dist.get('C', 0):,} ({ont_dist.get('C', 0)/len(train_terms)*100:.1f}%)\n\nAnnotation Statistics:\n  • Mean terms/protein: {terms_per_protein.mean():.1f}\n  • Median terms/protein: {terms_per_protein.median():.0f}\n  • Max terms/protein: {terms_per_protein.max()}\n  • Mean proteins/term: {proteins_per_term.mean():.1f}\n  • Median proteins/term: {proteins_per_term.median():.0f}\n\"\"\"\nax8.text(0.05, 0.5, summary_text, fontsize=10, family='monospace',\n         verticalalignment='center', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))\n\nplt.suptitle('Training Data Comprehensive Analysis', fontsize=15, fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n# Continue with the rest of the code (sequences, features, training, etc.)\nprint(\"\\n   Loading sequences (this may take a while for 100% of data)...\")\nprint(f\"   Expected proteins: {train_terms['protein'].nunique():,}\")\n\ntrain_seqs = {}\nloaded_count = 0\ntarget_proteins = set(train_terms['protein'].unique())\n\nfor rec in SeqIO.parse(TRAIN_DIR / 'train_sequences.fasta', 'fasta'):\n    pid = rec.id.split('|')[1] if '|' in rec.id else rec.id\n    if pid in target_proteins:\n        train_seqs[pid] = str(rec.seq)\n        loaded_count += 1\n        \n        # Progress indicator\n        if loaded_count % 10000 == 0:\n            print(f\"      Loaded {loaded_count:,} sequences...\")\n        \n    if loaded_count >= len(target_proteins):\n        break\n\nprint(f\"   ✓ Loaded {len(train_seqs):,} training sequences\")\n\n# Enhanced sequence analysis\nseq_lengths = [len(s) for s in train_seqs.values()]\nprint(f\"   ✓ Sequence length: mean={np.mean(seq_lengths):.0f}, \"\n      f\"median={np.median(seq_lengths):.0f}, range=[{min(seq_lengths)}-{max(seq_lengths)}]\")\n\nprint(\"\\n✅ Data loading complete! Ready for feature extraction and training.\")\nprint(f\"   Total proteins: {len(train_seqs):,}\")\nprint(f\"   Total annotations: {len(train_terms):,}\")\nprint(f\"   Total GO terms: {train_terms['term'].nunique():,}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-01T09:25:09.311693Z","iopub.execute_input":"2025-12-01T09:25:09.312313Z","iopub.status.idle":"2025-12-01T09:25:36.832681Z","shell.execute_reply.started":"2025-12-01T09:25:09.312289Z","shell.execute_reply":"2025-12-01T09:25:36.831823Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# CAFA 6 - COMPREHENSIVE DATA ANALYSIS\n# Chạy code này để hiểu sâu về data trước khi cải tiến model\n# ============================================================================\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom pathlib import Path\nfrom collections import Counter, defaultdict\nfrom tqdm.auto import tqdm\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Paths\nBASE = Path('/kaggle/input/cafa-6-protein-function-prediction')\nTRAIN_DIR = BASE / 'Train'\nTEST_DIR = BASE / 'Test'\n\nprint(\"=\"*80)\nprint(\"🔬 CAFA 6 - DEEP DATA ANALYSIS\")\nprint(\"=\"*80)\n\n# ============================================================================\n# 1. BASIC STATISTICS\n# ============================================================================\nprint(\"\\n[1/10] Loading basic data...\")\n\ntrain_terms = pd.read_csv(TRAIN_DIR / 'train_terms.tsv', sep='\\t', \n                          header=None, names=['protein', 'term', 'aspect'])\ntrain_taxonomy = pd.read_csv(TRAIN_DIR / 'train_taxonomy.tsv', sep='\\t',\n                             header=None, names=['protein', 'taxon'])\nia_df = pd.read_csv(BASE / 'IA.tsv', sep='\\t', header=None, names=['term', 'ia'])\n\nprint(f\"✓ Training annotations: {len(train_terms):,}\")\nprint(f\"✓ Unique proteins: {train_terms['protein'].nunique():,}\")\nprint(f\"✓ Unique GO terms: {train_terms['term'].nunique():,}\")\nprint(f\"✓ Unique species: {train_taxonomy['taxon'].nunique()}\")\n\n# ============================================================================\n# 2. TERM FREQUENCY ANALYSIS (CRITICAL!)\n# ============================================================================\nprint(\"\\n[2/10] Analyzing term frequencies...\")\n\nterm_freq = train_terms['term'].value_counts()\nproteins_per_term = train_terms.groupby('term')['protein'].nunique()\n\nprint(\"\\n📊 TERM FREQUENCY DISTRIBUTION:\")\nprint(f\"   Total terms: {len(term_freq):,}\")\n\n# 🔧 FIX: Frequency bins - số labels phải = số bins - 1\nfreq_bins = [1, 5, 10, 20, 50, 100, 500, 1000, float('inf')]\nfreq_labels = ['1-4', '5-9', '10-19', '20-49', '50-99', '100-499', '500-999', '1000+']\n\n# ✅ Đây là cách đúng\nterm_freq_binned = pd.cut(term_freq, bins=freq_bins, \n                           labels=freq_labels, right=False)\nfreq_dist = term_freq_binned.value_counts().sort_index()\n\nprint(\"\\n   Terms by frequency range:\")\nfor label, count in freq_dist.items():\n    pct = count / len(term_freq) * 100\n    print(f\"     {label:10s}: {count:6,} terms ({pct:5.1f}%)\")\n\n# Top frequent terms\nprint(\"\\n   🔥 TOP 20 MOST FREQUENT TERMS:\")\ntop_20 = term_freq.head(20)\nfor i, (term, count) in enumerate(top_20.items(), 1):\n    pct = count / train_terms['protein'].nunique() * 100\n    print(f\"     {i:2d}. {term}: {count:6,} proteins ({pct:5.1f}%)\")\n\n# Rare terms analysis\nrare_terms = term_freq[term_freq < 10]\nprint(f\"\\n   ⚠️  RARE TERMS (< 10 proteins): {len(rare_terms):,} ({len(rare_terms)/len(term_freq)*100:.1f}%)\")\nprint(f\"   ⚠️  VERY RARE (< 5 proteins): {len(term_freq[term_freq < 5]):,}\")\n\n# ============================================================================\n# 3. ONTOLOGY-SPECIFIC ANALYSIS\n# ============================================================================\nprint(\"\\n[3/10] Analyzing by ontology...\")\n\naspect_map = {'F': 'MFO', 'P': 'BPO', 'C': 'CCO'}\ntrain_terms['ontology'] = train_terms['aspect'].map(aspect_map)\n\nprint(\"\\n📊 ONTOLOGY BREAKDOWN:\")\nfor ont in ['MF', 'BP', 'CC']:\n    ont_data = train_terms[train_terms['ontology'] == ont]\n    n_terms = ont_data['term'].nunique()\n    n_annot = len(ont_data)\n    avg_per_protein = n_annot / train_terms['protein'].nunique()\n    \n    print(f\"\\n   {ont} (Molecular Function / Biological Process / Cellular Component):\")\n    print(f\"     - Unique terms: {n_terms:,}\")\n    print(f\"     - Total annotations: {n_annot:,}\")\n    print(f\"     - Avg per protein: {avg_per_protein:.2f}\")\n    \n    # Frequency stats per ontology\n    ont_term_freq = ont_data['term'].value_counts()\n    print(f\"     - Median frequency: {ont_term_freq.median():.0f}\")\n    print(f\"     - Terms with <10 proteins: {len(ont_term_freq[ont_term_freq < 10]):,}\")\n    print(f\"     - Terms with >100 proteins: {len(ont_term_freq[ont_term_freq > 100]):,}\")\n\n# ============================================================================\n# 4. PROTEIN ANNOTATION DENSITY - FIXED\n# ============================================================================\nprint(\"\\n[4/10] Analyzing protein annotation patterns...\")\n\nterms_per_protein = train_terms.groupby('protein').size()\naspects_per_protein = train_terms.groupby('protein')['aspect'].nunique()\n\nprint(\"\\n📊 PROTEIN ANNOTATION DENSITY:\")\nprint(f\"   Mean terms/protein: {terms_per_protein.mean():.1f}\")\nprint(f\"   Median terms/protein: {terms_per_protein.median():.0f}\")\nprint(f\"   Std dev: {terms_per_protein.std():.1f}\")\nprint(f\"   Range: [{terms_per_protein.min()}, {terms_per_protein.max()}]\")\n\n# 🔧 FIX: Density bins\ndensity_bins = [0, 5, 10, 20, 50, 100, float('inf')]\ndensity_labels = ['1-5', '6-10', '11-20', '21-50', '51-100', '100+']\n\ndensity_binned = pd.cut(terms_per_protein, bins=density_bins,\n                        labels=density_labels, right=False)\ndensity_dist = density_binned.value_counts().sort_index()\n\nprint(\"\\n   Proteins by annotation count:\")\nfor label, count in density_dist.items():\n    pct = count / len(terms_per_protein) * 100\n    print(f\"     {label:10s}: {count:6,} proteins ({pct:5.1f}%)\")\n\n# Sparse proteins\nsparse_proteins = terms_per_protein[terms_per_protein < 10]\nprint(f\"\\n   ⚠️  SPARSE PROTEINS (<10 terms): {len(sparse_proteins):,} ({len(sparse_proteins)/len(terms_per_protein)*100:.1f}%)\")\n\n# Incomplete proteins\nincomplete_proteins = aspects_per_protein[aspects_per_protein < 3]\nprint(f\"   ⚠️  INCOMPLETE (missing ontology): {len(incomplete_proteins):,} ({len(incomplete_proteins)/len(aspects_per_protein)*100:.1f}%)\")\n\n# ============================================================================\n# 5. SPECIES ANALYSIS\n# ============================================================================\nprint(\"\\n[5/10] Analyzing species distribution...\")\n\nspecies_counts = train_taxonomy['taxon'].value_counts()\nprint(f\"\\n📊 TOP 15 SPECIES:\")\nfor i, (taxon, count) in enumerate(species_counts.head(15).items(), 1):\n    pct = count / len(train_taxonomy) * 100\n    print(f\"   {i:2d}. Taxon {taxon}: {count:6,} proteins ({pct:5.1f}%)\")\n\n# Species diversity\nprint(f\"\\n   Total species: {len(species_counts)}\")\nprint(f\"   Species with <10 proteins: {len(species_counts[species_counts < 10]):,}\")\nprint(f\"   Top 10 species coverage: {species_counts.head(10).sum()/len(train_taxonomy)*100:.1f}%\")\n\n# ============================================================================\n# 6. INFORMATION ACCRETION ANALYSIS\n# ============================================================================\nprint(\"\\n[6/10] Analyzing IA weights...\")\n\nia_dict = dict(zip(ia_df['term'], ia_df['ia']))\n\n# Add IA to train_terms\ntrain_terms['ia'] = train_terms['term'].map(ia_dict)\n\nprint(\"\\n📊 IA WEIGHT STATISTICS:\")\nprint(f\"   Mean IA: {ia_df['ia'].mean():.4f}\")\nprint(f\"   Median IA: {ia_df['ia'].median():.4f}\")\nprint(f\"   Std dev: {ia_df['ia'].std():.4f}\")\nprint(f\"   Range: [{ia_df['ia'].min():.4f}, {ia_df['ia'].max():.4f}]\")\n\n# IA by ontology\nprint(\"\\n   IA by ontology:\")\nfor ont in ['MF', 'BP', 'CC']:\n    ont_terms = train_terms[train_terms['ontology'] == ont]['term'].unique()\n    ont_ia = [ia_dict.get(t, 0) for t in ont_terms if t in ia_dict]\n    print(f\"     {ont}: {np.mean(ont_ia):.4f} ± {np.std(ont_ia):.4f}\")\n\n# High IA terms (hard to predict)\nhigh_ia_terms = ia_df[ia_df['ia'] > ia_df['ia'].quantile(0.9)]\nprint(f\"\\n   ⚠️  HIGH IA TERMS (top 10%): {len(high_ia_terms):,}\")\nprint(f\"       These are HARD to predict but worth MORE points!\")\n\n# ============================================================================\n# 7. TERM CO-OCCURRENCE PATTERNS\n# ============================================================================\nprint(\"\\n[7/10] Analyzing term co-occurrence...\")\n\n# Sample for speed\nsample_proteins = train_terms['protein'].unique()[:5000]\nsample_terms = train_terms[train_terms['protein'].isin(sample_proteins)]\n\n# Build co-occurrence matrix for top 50 terms\ntop_50_terms = term_freq.head(50).index\nprotein_to_terms_dict = sample_terms.groupby('protein')['term'].apply(set).to_dict()\n\ncooc_counts = Counter()\nfor terms in tqdm(protein_to_terms_dict.values(), desc=\"Computing co-occurrence\", leave=False):\n    terms_list = [t for t in terms if t in top_50_terms]\n    for i, t1 in enumerate(terms_list):\n        for t2 in terms_list[i+1:]:\n            pair = tuple(sorted([t1, t2]))\n            cooc_counts[pair] += 1\n\nprint(f\"\\n📊 TERM CO-OCCURRENCE (top 50 terms):\")\nprint(f\"   Total pairs: {len(cooc_counts):,}\")\ntop_cooc = sorted(cooc_counts.items(), key=lambda x: x[1], reverse=True)[:10]\nprint(f\"\\n   Top 10 co-occurring pairs:\")\nfor i, ((t1, t2), count) in enumerate(top_cooc, 1):\n    print(f\"     {i:2d}. ({t1}, {t2}): {count:,} times\")\n\n# ============================================================================\n# 8. IMBALANCE ANALYSIS - FIXED\n# ============================================================================\nprint(\"\\n[8/10] Analyzing class imbalance...\")\n\n# Calculate positive rate\npositive_rates = {}\ntotal_proteins = train_terms['protein'].nunique()\n\nfor term in term_freq.index:\n    n_positive = term_freq[term]\n    positive_rate = n_positive / total_proteins\n    positive_rates[term] = positive_rate\n\npositive_rates_series = pd.Series(positive_rates)\n\nprint(\"\\n📊 CLASS IMBALANCE SEVERITY:\")\nprint(f\"   Mean positive rate: {positive_rates_series.mean():.4f}\")\nprint(f\"   Median positive rate: {positive_rates_series.median():.4f}\")\n\n# 🔧 FIX: Imbalance bins\nimbalance_bins = [0, 0.001, 0.01, 0.05, 0.1, 0.5, 1.0]\nimbalance_labels = ['<0.1%', '0.1-1%', '1-5%', '5-10%', '10-50%', '>50%']\n\nimbalance_binned = pd.cut(positive_rates_series, bins=imbalance_bins,\n                          labels=imbalance_labels, right=False)\nimbalance_dist = imbalance_binned.value_counts().sort_index()\n\nprint(\"\\n   Terms by positive rate:\")\nfor label, count in imbalance_dist.items():\n    pct = count / len(positive_rates_series) * 100\n    print(f\"     {label:10s}: {count:6,} terms ({pct:5.1f}%)\")\n\n# Extreme imbalance\nextreme_imbalance = positive_rates_series[positive_rates_series < 0.01]\nprint(f\"\\n   ⚠️  EXTREME IMBALANCE (<1% positive): {len(extreme_imbalance):,} terms\")\nprint(f\"       These need special handling (SMOTE, focal loss, etc.)\")\n\n# ============================================================================\n# 9. EMBEDDING AVAILABILITY CHECK\n# ============================================================================\nprint(\"\\n[9/10] Checking embedding availability...\")\n\n# Check if embeddings exist\nesm2_path = '/kaggle/input/cafa-5-ems-2-embeddings-numpy'\nprotbert_path = '/kaggle/input/protbert-embeddings-for-cafa5'\nt5_path = '/kaggle/input/t5embeds'\n\navailable_embeddings = []\n\nif Path(esm2_path).exists():\n    esm2_train = np.load(f'{esm2_path}/train_ids.npy', allow_pickle=True)\n    esm2_test = np.load(f'{esm2_path}/test_ids.npy', allow_pickle=True)\n    available_embeddings.append(('ESM2', len(esm2_train), len(esm2_test)))\n    print(f\"   ✓ ESM2: {len(esm2_train):,} train, {len(esm2_test):,} test\")\n\nif Path(protbert_path).exists():\n    protbert_train = np.load(f'{protbert_path}/train_ids.npy', allow_pickle=True)\n    protbert_test = np.load(f'{protbert_path}/test_ids.npy', allow_pickle=True)\n    available_embeddings.append(('ProtBERT', len(protbert_train), len(protbert_test)))\n    print(f\"   ✓ ProtBERT: {len(protbert_train):,} train, {len(protbert_test):,} test\")\n\nif Path(t5_path).exists():\n    t5_train = np.load(f'{t5_path}/train_ids.npy', allow_pickle=True)\n    t5_test = np.load(f'{t5_path}/test_ids.npy', allow_pickle=True)\n    available_embeddings.append(('T5', len(t5_train), len(t5_test)))\n    print(f\"   ✓ T5: {len(t5_train):,} train, {len(t5_test):,} test\")\n\nprint(f\"\\n   Total embedding sources: {len(available_embeddings)}\")\n\n# ============================================================================\n# 10. KEY INSIGHTS & RECOMMENDATIONS\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(\"🎯 KEY INSIGHTS & RECOMMENDATIONS\")\nprint(\"=\"*80)\n\ninsights = []\n\n# 1. Frequency analysis\nvery_rare_pct = len(term_freq[term_freq < 5]) / len(term_freq) * 100\nif very_rare_pct > 30:\n    insights.append({\n        'issue': 'HIGH PERCENTAGE OF RARE TERMS',\n        'detail': f'{very_rare_pct:.1f}% of terms have <5 proteins',\n        'recommendation': 'Use hierarchical classification or focus on top-K frequent terms'\n    })\n\n# 2. Sparse proteins\nsparse_pct = len(sparse_proteins) / len(terms_per_protein) * 100\nif sparse_pct > 20:\n    insights.append({\n        'issue': 'MANY SPARSE PROTEINS',\n        'detail': f'{sparse_pct:.1f}% of proteins have <10 annotations',\n        'recommendation': 'Implement transfer learning from well-annotated proteins'\n    })\n\n# 3. Class imbalance\nextreme_imbalance_pct = len(extreme_imbalance) / len(positive_rates_series) * 100\nif extreme_imbalance_pct > 50:\n    insights.append({\n        'issue': 'SEVERE CLASS IMBALANCE',\n        'detail': f'{extreme_imbalance_pct:.1f}% of terms have <1% positive rate',\n        'recommendation': 'Use focal loss, class weights, or SMOTE'\n    })\n\n# 4. Multi-embedding\nif len(available_embeddings) > 1:\n    insights.append({\n        'issue': 'MULTIPLE EMBEDDINGS AVAILABLE',\n        'detail': f'{len(available_embeddings)} embedding sources found',\n        'recommendation': 'Create ensemble from different embeddings'\n    })\n\n# 5. Ontology imbalance\nont_counts = train_terms['ontology'].value_counts()\nmax_ont = ont_counts.max()\nmin_ont = ont_counts.min()\nif max_ont / min_ont > 2:\n    insights.append({\n        'issue': 'ONTOLOGY IMBALANCE',\n        'detail': f'Ratio {max_ont/min_ont:.1f}:1 between ontologies',\n        'recommendation': 'Train separate models per ontology'\n    })\n\n# Print insights\nfor i, insight in enumerate(insights, 1):\n    print(f\"\\n{i}. 🔴 {insight['issue']}\")\n    print(f\"   📊 {insight['detail']}\")\n    print(f\"   💡 {insight['recommendation']}\")\n\n# ============================================================================\n# SUMMARY STATISTICS FOR DECISION MAKING\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(\"📋 SUMMARY FOR MODEL DESIGN\")\nprint(\"=\"*80)\n\nsummary = {\n    'Total proteins': train_terms['protein'].nunique(),\n    'Total terms': train_terms['term'].nunique(),\n    'Avg terms/protein': terms_per_protein.mean(),\n    'Avg proteins/term': proteins_per_term.mean(),\n    'Terms with <10 proteins': len(term_freq[term_freq < 10]),\n    'Terms with >100 proteins': len(term_freq[term_freq > 100]),\n    'Sparse proteins (<10 terms)': len(sparse_proteins),\n    'Incomplete proteins': len(incomplete_proteins),\n    'Mean IA weight': ia_df['ia'].mean(),\n    'Extreme imbalance terms': len(extreme_imbalance),\n    'Available embeddings': len(available_embeddings),\n}\n\nfor key, value in summary.items():\n    if isinstance(value, float):\n        print(f\"   {key:30s}: {value:.2f}\")\n    else:\n        print(f\"   {key:30s}: {value:,}\")\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"✅ ANALYSIS COMPLETE!\")\nprint(\"=\"*80)\nprint(\"\\n💡 Next steps:\")\nprint(\"   1. Review insights above\")\nprint(\"   2. Implement recommended techniques\")\nprint(\"   3. Focus on high-IA terms for better score\")\nprint(\"   4. Consider ensemble from multiple embeddings\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-01T09:25:51.961284Z","iopub.execute_input":"2025-12-01T09:25:51.961882Z","iopub.status.idle":"2025-12-01T09:25:54.401388Z","shell.execute_reply.started":"2025-12-01T09:25:51.961859Z","shell.execute_reply":"2025-12-01T09:25:54.400562Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================================\n# CAFA 6 - FIXED VERSION (ALL BUGS RESOLVED)\n# ============================================================================\n\nimport numpy as np\nimport pandas as pd\nimport os\nimport gc\nimport sys\nimport subprocess\nfrom pathlib import Path\nfrom collections import defaultdict, Counter\nfrom tqdm.auto import tqdm\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# ============================================================================\n# CONFIG\n# ============================================================================\nCONFIG = {\n    # Data sampling\n    'MAX_TRAIN_SAMPLES': 120000,\n    'TOP_K_FREQUENT': 1000, # từ 650 và tăng điểm\n    'TOP_K_RARE': 1200, # từ 800 và tăng điểm\n    'MIN_FREQ_COMMON': 20,\n    'MIN_FREQ_RARE': 3,\n    \n    # Focal Loss parameters\n    'FOCAL_LOSS_GAMMA': 2.5,\n    'FOCAL_LOSS_ALPHA': 0.3,\n    'USE_IA_WEIGHTS': True,\n    \n    # Multi-embedding ensemble\n    'USE_MULTI_EMBEDDING': True,\n    'EMBEDDING_SOURCES': ['esm2', 'protbert', 't5'],\n    'EMBEDDING_FUSION': 'attention',\n    \n    # High-IA term prioritization\n    'IA_SAMPLING_RATIO': 0.4,\n    'HIGH_IA_THRESHOLD': 2.5,\n    'IA_CLIP_MIN': None,        # ✅ BỎ clipping!\n    'IA_CLIP_MAX': None,\n    'IA_TRANSFORM': 'log1p',\n    \n    # Co-occurrence modeling\n    'USE_COOCCURRENCE': True,\n    'COOCCUR_TOP_K': 50,\n    \n    # Model architecture\n    'HIDDEN_DIMS': [512, 256],\n    'DROPOUT_RATE': 0.3,\n    'USE_BATCH_NORM': True,\n    \n    # Training\n    'EPOCHS': 35,\n    'BATCH_SIZE': 64,\n    'LEARNING_RATE': 1e-3,\n    'WEIGHT_DECAY': 1e-4,\n    'GRAD_CLIP': 1.0,\n    'WARMUP_EPOCHS': 5,\n    \n    # Prediction\n    'MAX_PREDS_PER_PROTEIN': 2000,\n    'MIN_CONFIDENCE': 0.005,# từ 0.005\n    'TEMPERATURE': 0.8,\n    \n    # GO propagation\n    'USE_GO_PROPAGATION': True,\n    'PROPAGATION_DECAY': 0.70,\n    \n    # Ontology calibration\n    'ONTOLOGY_CALIBRATION': {\n        'MFO': 1.10,\n        'BPO': 1.00,\n        'CCO': 1.05\n    },\n    \n    # Ensemble weights\n    'BASE_BLAST_WEIGHT': 0.55, # best\n    'BASE_DL_WEIGHT': 0.35, # best\n    'BASE_FREQ_WEIGHT': 0.10, # best\n    \n    # Paths\n    'BASE_PATH': '/kaggle/input/cafa-6-protein-function-prediction',\n    'ESM2_PATH': 'cafa-5-ems-2-embeddings-numpy',\n    'PROTBERT_PATH': 'protbert-embeddings-for-cafa5',\n    'T5_PATH': 't5embeds',\n    'BLAST_PATH': '/kaggle/input/blast-quick-sprof-zero-pred/submission.tsv',\n    'RANDOM_SEED': 42,\n}\n\nprint(\"=\"*80)\nprint(\"🚀 CAFA 6 - FIXED VERSION\")\nprint(\"=\"*80)\n\n# ============================================================================\n# SETUP\n# ============================================================================\ndef install(package):\n    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", package])\n\ntry:\n    import torch\n    import torch.nn as nn\n    import torch.nn.functional as F\nexcept:\n    install('torch')\n    import torch\n    import torch.nn as nn\n    import torch.nn.functional as F\n\ntry:\n    import networkx as nx\n    import obonet\nexcept:\n    install('networkx')\n    install('obonet')\n    import networkx as nx\n    import obonet\n\nnp.random.seed(CONFIG['RANDOM_SEED'])\ntorch.manual_seed(CONFIG['RANDOM_SEED'])\nif torch.cuda.is_available():\n    torch.cuda.manual_seed(CONFIG['RANDOM_SEED'])\n    \ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(f\"🎯 Device: {device}\")\n\n# ============================================================================\n# FOCAL LOSS\n# ============================================================================\nclass FocalLoss(nn.Module):\n    def __init__(self, gamma=2.0, alpha=0.25, reduction='mean'):\n        super().__init__()\n        self.gamma = gamma\n        self.alpha = alpha\n        self.reduction = reduction\n    \n    def forward(self, inputs, targets):\n        bce_loss = F.binary_cross_entropy(inputs, targets, reduction='none')\n        p_t = inputs * targets + (1 - inputs) * (1 - targets)\n        focal_weight = (1 - p_t) ** self.gamma\n        \n        if self.alpha is not None:\n            alpha_t = self.alpha * targets + (1 - self.alpha) * (1 - targets)\n            focal_loss = alpha_t * focal_weight * bce_loss\n        else:\n            focal_loss = focal_weight * bce_loss\n        \n        if self.reduction == 'mean':\n            return focal_loss.mean()\n        elif self.reduction == 'sum':\n            return focal_loss.sum()\n        else:\n            return focal_loss\n\n# ============================================================================\n# LOAD GO ONTOLOGY\n# ============================================================================\nprint(\"\\n[1/12] Loading GO ontology...\")\nBASE = Path(CONFIG['BASE_PATH'])\nTRAIN_DIR = BASE / 'Train'\n\ngo_graph = None\nancestor_cache = {}\n\nif CONFIG['USE_GO_PROPAGATION']:\n    try:\n        go_graph = obonet.read_obo(TRAIN_DIR / 'go-basic.obo')\n        print(f\"   ✓ Loaded {len(go_graph):,} GO terms\")\n        \n        def get_ancestors(term_id, graph):\n            ancestors = set()\n            try:\n                for node in nx.descendants(graph, term_id):\n                    ancestors.add(node)\n            except:\n                pass\n            return ancestors\n    except Exception as e:\n        print(f\"   ⚠️ Could not load GO graph: {e}\")\n        CONFIG['USE_GO_PROPAGATION'] = False\n\n# ============================================================================\n# 🔧 FIX #2: LOAD IA WEIGHTS (CORRECT FILENAME)\n# ============================================================================\nprint(\"\\n[2/12] Loading IA weights...\")\nia_weights_dict = {}\n\ntry:\n    # Try .txt first (correct format), then .tsv\n    ia_file = BASE / 'IA.txt'\n    if not ia_file.exists():\n        ia_file = BASE / 'IA.tsv'\n    \n    ia_df = pd.read_csv(ia_file, sep='\\t', header=None, names=['term', 'ia_weight'])\n    ia_weights_dict = dict(zip(ia_df['term'], ia_df['ia_weight']))\n    \n    ia_values = list(ia_weights_dict.values())\n    print(f\"   ✓ Loaded IA weights for {len(ia_weights_dict):,} terms\")\n    print(f\"   📊 Mean IA: {np.mean(ia_values):.2f}\")\n    print(f\"   📊 Max IA: {np.max(ia_values):.2f}\")\n    \n    high_ia_terms = [t for t, ia in ia_weights_dict.items() \n                     if ia >= CONFIG['HIGH_IA_THRESHOLD']]\n    print(f\"   🔥 High-IA terms: {len(high_ia_terms):,}\")\n    \nexcept Exception as e:\n    print(f\"   ⚠️ Could not load IA weights: {e}\")\n    print(f\"   → Using uniform weights\")\n    CONFIG['USE_IA_WEIGHTS'] = False\n    high_ia_terms = []\n\n# ============================================================================\n# TWO-STAGE TERM SELECTION\n# ============================================================================\nprint(\"\\n[3/12] Loading annotations with two-stage selection...\")\ntrain_terms_df = pd.read_csv(TRAIN_DIR / 'train_terms.tsv', sep='\\t', \n                              header=None, names=['protein', 'term', 'aspect'])\n\nterm_freq = train_terms_df['term'].value_counts()\nprint(f\"   Total unique terms: {len(term_freq):,}\")\n\n# Stage 1: Frequent terms\nfrequent_terms = term_freq[term_freq >= CONFIG['MIN_FREQ_COMMON']].index.tolist()\nfrequent_terms = frequent_terms[:CONFIG['TOP_K_FREQUENT']]\nprint(f\"   📊 Stage 1 - Frequent terms: {len(frequent_terms):,}\")\n\n# Stage 2: High-IA rare terms\nrare_candidates = term_freq[\n    (term_freq >= CONFIG['MIN_FREQ_RARE']) & \n    (term_freq < CONFIG['MIN_FREQ_COMMON'])\n].index.tolist()\n\nif CONFIG['USE_IA_WEIGHTS'] and high_ia_terms:\n    rare_with_ia = [(t, ia_weights_dict.get(t, 0)) for t in rare_candidates \n                    if t in ia_weights_dict]\n    rare_with_ia.sort(key=lambda x: x[1], reverse=True)\n    rare_high_ia = [t for t, ia in rare_with_ia[:CONFIG['TOP_K_RARE']]]\n    print(f\"   📊 Stage 2 - High-IA rare terms: {len(rare_high_ia):,}\")\nelse:\n    rare_high_ia = []\n    print(f\"   📊 Stage 2 - High-IA rare terms: 0 (no IA data)\")\n\n# Combine\ntop_terms = list(set(frequent_terms + rare_high_ia))\n\nprint(f\"   ✅ Total selected terms: {len(top_terms):,}\")\n# ✅ BUILD CACHE SAU KHI CÓ top_terms\nif CONFIG['USE_GO_PROPAGATION'] and go_graph:\n    print(\"   Building ancestor cache for selected terms...\")\n    for term in tqdm(top_terms, desc=\"   Caching\", leave=False):\n        if term in go_graph:\n            ancestor_cache[term] = get_ancestors(term, go_graph)\n    print(f\"   ✓ Cached {len(ancestor_cache):,}/{len(top_terms):,} term ancestors\")\ntrain_terms_df = train_terms_df[train_terms_df['term'].isin(top_terms)]\n\n# Create mappings\nprotein_to_terms = train_terms_df.groupby('protein')['term'].apply(list).to_dict()\nterm_to_aspect = dict(zip(train_terms_df['term'], train_terms_df['aspect']))\nterm_to_idx = {term: idx for idx, term in enumerate(top_terms)}\nidx_to_term = {idx: term for term, idx in term_to_idx.items()}\n\nprint(f\"   ✓ {len(protein_to_terms):,} proteins with annotations\")\n\n# Create IA weight vector\nia_weight_vector = np.ones(len(top_terms))\n\nif CONFIG['USE_IA_WEIGHTS']:\n    for idx, term in enumerate(top_terms):\n        ia_weight_vector[idx] = ia_weights_dict.get(term, 1.0)\n    \n    # ✅ Sử dụng log1p để soft-scale (giữ được tín hiệu)\n    ia_weight_vector = np.log1p(ia_weight_vector)  # log(1+x)\n    \n    # ✅ Normalize về mean=1 (không clip!)\n    ia_weight_vector = ia_weight_vector / ia_weight_vector.mean()\n    \n    # ✅ Chỉ clip extreme outliers nếu cần (rất rộng)\n    ia_weight_vector = np.clip(ia_weight_vector, 0.1, 10.0)  # Thay vì 0.5-2.0\n    \n    print(f\"   ✓ IA weights (log-scaled): mean={ia_weight_vector.mean():.2f}, \"\n          f\"range=[{ia_weight_vector.min():.2f}, {ia_weight_vector.max():.2f}]\")\n\n# ============================================================================\n# CO-OCCURRENCE MATRIX\n# ============================================================================\nprint(\"\\n[4/12] Building co-occurrence matrix...\")\n\nif CONFIG['USE_COOCCURRENCE']:\n    cooccur_matrix = np.zeros((len(top_terms), len(top_terms)))\n    \n    for terms in tqdm(protein_to_terms.values(), desc=\"   Computing\", leave=False):\n        term_indices = [term_to_idx[t] for t in terms if t in term_to_idx]\n        \n        for i in term_indices:\n            for j in term_indices:\n                if i != j:\n                    cooccur_matrix[i, j] += 1\n    \n    row_sums = cooccur_matrix.sum(axis=1, keepdims=True)\n    row_sums[row_sums == 0] = 1\n    cooccur_matrix = cooccur_matrix / row_sums\n    \n    print(f\"   ✓ Co-occurrence matrix: {cooccur_matrix.shape}\")\n    print(f\"   📊 Sparsity: {(cooccur_matrix > 0).sum() / cooccur_matrix.size * 100:.1f}%\")\nelse:\n    cooccur_matrix = None\n\n# ============================================================================\n# 🔧 FIX #3: LOAD EMBEDDINGS WITH FALLBACK\n# ============================================================================\nprint(\"\\n[5/12] Loading embeddings...\")\n\nembedding_dicts = {}\nembedding_dims = []\n\n# ESM2 (always available)\nesm2_base = f\"/kaggle/input/{CONFIG['ESM2_PATH']}\"\ntrain_ids = np.load(f\"{esm2_base}/train_ids.npy\", allow_pickle=True)\ntrain_embeds = np.load(f\"{esm2_base}/train_embeddings.npy\")\ntest_ids = np.load(f\"{esm2_base}/test_ids.npy\", allow_pickle=True)\ntest_embeds = np.load(f\"{esm2_base}/test_embeddings.npy\")\n\nembedding_dicts['esm2'] = {\n    'train': {str(pid): emb for pid, emb in zip(train_ids, train_embeds)},\n    'test': {str(pid): emb for pid, emb in zip(test_ids, test_embeds)}\n}\nembedding_dims.append(train_embeds.shape[1])\nprint(f\"   ✓ ESM2: dim={train_embeds.shape[1]}\")\n\ndel train_ids, train_embeds, test_ids, test_embeds\ngc.collect()\n\n# Try other embeddings\nif CONFIG['USE_MULTI_EMBEDDING']:\n    for emb_name, path_key in [('protbert', 'PROTBERT_PATH'), ('t5', 'T5_PATH')]:\n        try:\n            emb_path = f\"/kaggle/input/{CONFIG[path_key]}\"\n            \n            # ✅ Kiểm tra filename khác nhau\n            if emb_name == 't5':\n                # T5 dùng \"embeds\" thay vì \"embeddings\"\n                train_ids = np.load(f\"{emb_path}/train_ids.npy\", allow_pickle=True)\n                train_embeds = np.load(f\"{emb_path}/train_embeds.npy\")\n                test_ids = np.load(f\"{emb_path}/test_ids.npy\", allow_pickle=True)\n                test_embeds = np.load(f\"{emb_path}/test_embeds.npy\")\n            else:\n                # ProtBERT dùng \"embeddings\"\n                train_ids = np.load(f\"{emb_path}/train_ids.npy\", allow_pickle=True)\n                train_embeds = np.load(f\"{emb_path}/train_embeddings.npy\")\n                test_ids = np.load(f\"{emb_path}/test_ids.npy\", allow_pickle=True)\n                test_embeds = np.load(f\"{emb_path}/test_embeddings.npy\")\n            \n            embedding_dicts[emb_name] = {\n                'train': {str(pid): emb for pid, emb in zip(train_ids, train_embeds)},\n                'test': {str(pid): emb for pid, emb in zip(test_ids, test_embeds)}\n            }\n                \n            embedding_dims.append(train_embeds.shape[1])\n            print(f\"   ✓ {emb_name.upper()}: dim={train_embeds.shape[1]}\")\n            \n            del train_ids, train_embeds, test_ids, test_embeds\n            gc.collect()\n            \n        except Exception as e:\n            print(f\"   ⚠️ Could not load {emb_name}: {e}\")  # ← Debug\n            pass\n\nprint(f\"\\n   ✅ Total sources: {len(embedding_dicts)}\")\nprint(f\"   📊 Embedding dims: {embedding_dims}\")\n\n# ============================================================================\n# PREPARE TRAINING DATA\n# ============================================================================\nprint(\"\\n[6/12] Preparing training data...\")\n\nvalid_proteins = []\nfor p in protein_to_terms.keys():\n    has_all = all(p in embedding_dicts[src]['train'] \n                  for src in embedding_dicts.keys())\n    if has_all:\n        valid_proteins.append(p)\n\nvalid_proteins = valid_proteins[:CONFIG['MAX_TRAIN_SAMPLES']]\nprint(f\"   ✓ Valid proteins: {len(valid_proteins):,}\")\n\nfrom sklearn.preprocessing import MultiLabelBinarizer\nfrom sklearn.model_selection import StratifiedShuffleSplit\n\nmlb = MultiLabelBinarizer(classes=range(len(top_terms)))\ny_labels = [[term_to_idx[t] for t in protein_to_terms.get(p, []) if t in term_to_idx] \n            for p in valid_proteins]\ny_encoded = mlb.fit_transform(y_labels).astype(float)\n\nprotein_categories = []\nfor p in valid_proteins:\n    n_terms = len(protein_to_terms.get(p, []))\n    if n_terms < 5:\n        category = 'sparse'\n    elif n_terms < 15:\n        category = 'medium'\n    else:\n        category = 'rich'\n    protein_categories.append(category)\n\nsplitter = StratifiedShuffleSplit(n_splits=1, test_size=0.15, \n                                  random_state=CONFIG['RANDOM_SEED'])\ntrain_idx, val_idx = next(splitter.split(valid_proteins, protein_categories))\n\ntrain_proteins = [valid_proteins[i] for i in train_idx]\nval_proteins = [valid_proteins[i] for i in val_idx]\ny_train = y_encoded[train_idx]\ny_val = y_encoded[val_idx]\n\nprint(f\"   ✓ Train: {len(train_proteins):,} proteins\")\nprint(f\"   ✓ Val: {len(val_proteins):,} proteins\")\n\n# ============================================================================\n# 🔧 FIX #1: IMPROVED MODEL WITH PROPER MULTI-EMBEDDING HANDLING\n# ============================================================================\nprint(\"\\n[7/12] Building model...\")\n\nclass MultiEmbeddingFusion(nn.Module):\n    \"\"\"Fuse multiple embeddings using attention\"\"\"\n    def __init__(self, embedding_dims, output_dim):\n        super().__init__()\n        self.n_sources = len(embedding_dims)\n        self.projections = nn.ModuleList([\n            nn.Linear(dim, output_dim) for dim in embedding_dims\n        ])\n        if self.n_sources > 1:\n            self.attention = nn.Linear(output_dim, 1)\n    \n    def forward(self, embeddings):\n        # Handle both single and multiple embeddings\n        if not isinstance(embeddings, list):\n            embeddings = [embeddings]\n        \n        # Project all to same dimension\n        projected = [proj(emb) for proj, emb in zip(self.projections, embeddings)]\n        \n        if len(projected) == 1:\n            return projected[0]\n        \n        # Attention fusion\n        stacked = torch.stack(projected, dim=1)\n        attn_scores = self.attention(stacked).squeeze(-1)\n        attn_weights = F.softmax(attn_scores, dim=1).unsqueeze(-1)\n        fused = (stacked * attn_weights).sum(dim=1)\n        \n        return fused\n\nclass AdvancedProteinClassifier(nn.Module):\n    def __init__(self, embedding_dims, num_terms, hidden_dims, \n                 dropout, use_cooccurrence=False):\n        super().__init__()\n        \n        # Multi-embedding fusion\n        fusion_dim = hidden_dims[0]\n        self.fusion = MultiEmbeddingFusion(embedding_dims, fusion_dim)\n        \n        # Deep encoder\n        layers = []\n        in_dim = fusion_dim\n        \n        for hidden_dim in hidden_dims:\n            layers.extend([\n                nn.Linear(in_dim, hidden_dim),\n                nn.BatchNorm1d(hidden_dim),\n                nn.ReLU(),\n                nn.Dropout(dropout)\n            ])\n            in_dim = hidden_dim\n        \n        self.encoder = nn.Sequential(*layers)\n        self.output = nn.Linear(hidden_dims[-1], num_terms)\n        \n        # Initialize\n        for module in self.modules():\n            if isinstance(module, nn.Linear):\n                nn.init.xavier_uniform_(module.weight)\n                if module.bias is not None:\n                    nn.init.zeros_(module.bias)\n    \n    def forward(self, embeddings, cooccur_matrix=None, return_logits=False):\n        # Fuse embeddings (handles both list and single tensor)\n        fused = self.fusion(embeddings)\n        \n        # Encode\n        features = self.encoder(fused)\n        \n        # Output\n        logits = self.output(features)\n        \n        if return_logits:\n            return logits\n        return torch.sigmoid(logits)\n\nmodel = AdvancedProteinClassifier(\n    embedding_dims,\n    len(top_terms),\n    CONFIG['HIDDEN_DIMS'],\n    CONFIG['DROPOUT_RATE'],\n    use_cooccurrence=CONFIG['USE_COOCCURRENCE']\n).to(device)\n\ntotal_params = sum(p.numel() for p in model.parameters())\nprint(f\"   ✓ Parameters: {total_params:,}\")\nprint(f\"   ✓ Embedding sources: {len(embedding_dims)}\")\n\n# ============================================================================\n# TRAINING\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(f\"TRAINING ({CONFIG['EPOCHS']} EPOCHS)\")\nprint(\"=\"*80)\n\ncriterion = FocalLoss(\n    gamma=CONFIG['FOCAL_LOSS_GAMMA'],\n    alpha=CONFIG['FOCAL_LOSS_ALPHA'],\n    reduction='none'\n)\n\nia_weight_tensor = torch.FloatTensor(ia_weight_vector).to(device)\n\noptimizer = torch.optim.AdamW(\n    model.parameters(),\n    lr=CONFIG['LEARNING_RATE'],\n    weight_decay=CONFIG['WEIGHT_DECAY']\n)\n\ndef get_lr_lambda(epoch):\n    if epoch < CONFIG['WARMUP_EPOCHS']:\n        return (epoch + 1) / CONFIG['WARMUP_EPOCHS']\n    else:\n        progress = (epoch - CONFIG['WARMUP_EPOCHS']) / (CONFIG['EPOCHS'] - CONFIG['WARMUP_EPOCHS'])\n        return 0.5 * (1 + np.cos(np.pi * progress))\n\nscheduler = torch.optim.lr_scheduler.LambdaLR(optimizer, get_lr_lambda)\n\nbest_val_loss = float('inf')\ncooccur_tensor = torch.FloatTensor(cooccur_matrix).to(device) if cooccur_matrix is not None else None\n\nfor epoch in range(CONFIG['EPOCHS']):\n    # Training\n    model.train()\n    indices = np.random.permutation(len(train_proteins))\n    epoch_loss = 0\n    n_batches = 0\n    \n    for i in range(0, len(indices), CONFIG['BATCH_SIZE']):\n        batch_idx = indices[i:i + CONFIG['BATCH_SIZE']]\n        batch_proteins = [train_proteins[j] for j in batch_idx]\n        \n        # Gather embeddings\n        if len(embedding_dicts) == 1:\n            # Single source - pass as tensor\n            src = list(embedding_dicts.keys())[0]\n            embeddings = torch.FloatTensor([\n                embedding_dicts[src]['train'][p] for p in batch_proteins\n            ]).to(device)\n        else:\n            # Multiple sources - pass as list\n            embeddings = []\n            for src in embedding_dicts.keys():\n                emb_batch = torch.FloatTensor([\n                    embedding_dicts[src]['train'][p] for p in batch_proteins\n                ]).to(device)\n                embeddings.append(emb_batch)\n        \n        y_batch = torch.FloatTensor(y_train[batch_idx]).to(device)\n        \n        optimizer.zero_grad()\n        outputs = model(embeddings, cooccur_tensor)\n        \n        loss_per_sample = criterion(outputs, y_batch)\n        weighted_loss = (loss_per_sample * ia_weight_tensor).mean()\n        \n        weighted_loss.backward()\n        torch.nn.utils.clip_grad_norm_(model.parameters(), CONFIG['GRAD_CLIP'])\n        optimizer.step()\n        \n        epoch_loss += weighted_loss.item()\n        n_batches += 1\n    \n    # Validation\n    model.eval()\n    val_loss = 0\n    val_batches = 0\n    \n    with torch.no_grad():\n        for i in range(0, len(val_proteins), CONFIG['BATCH_SIZE']):\n            batch_proteins = val_proteins[i:i + CONFIG['BATCH_SIZE']]\n            \n            if len(embedding_dicts) == 1:\n                src = list(embedding_dicts.keys())[0]\n                embeddings = torch.FloatTensor([\n                    embedding_dicts[src]['train'][p] for p in batch_proteins\n                ]).to(device)\n            else:\n                embeddings = []\n                for src in embedding_dicts.keys():\n                    emb_batch = torch.FloatTensor([\n                        embedding_dicts[src]['train'][p] for p in batch_proteins\n                    ]).to(device)\n                    embeddings.append(emb_batch)\n            \n            y_batch = torch.FloatTensor(y_val[i:i + CONFIG['BATCH_SIZE']]).to(device)\n            \n            outputs = model(embeddings, cooccur_tensor)\n            loss_per_sample = criterion(outputs, y_batch)\n            weighted_loss = (loss_per_sample * ia_weight_tensor).mean()\n            \n            val_loss += weighted_loss.item()\n            val_batches += 1\n    \n    train_loss_avg = epoch_loss / n_batches\n    val_loss_avg = val_loss / val_batches\n    current_lr = optimizer.param_groups[0]['lr']\n    \n    scheduler.step()\n    gc.collect()\n    torch.cuda.empty_cache() \n    \n    if val_loss_avg < best_val_loss:\n        best_val_loss = val_loss_avg\n        marker = \"⭐\"\n    else:\n        marker = \"\"\n    \n    if (epoch + 1) % 5 == 0 or marker:\n        print(f\"Epoch {epoch+1:2d}: Train={train_loss_avg:.4f}, Val={val_loss_avg:.4f}, \"\n              f\"LR={current_lr:.6f} {marker}\")\n    \n    gc.collect()\n\nprint(f\"\\n✅ Best Val Loss: {best_val_loss:.4f}\")\n\n# ============================================================================\n# PREDICTIONS WITH ALL IMPROVEMENTS\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(\"GENERATING PREDICTIONS\")\nprint(\"=\"*80)\n\nmodel.eval()\n\n# Get all test protein IDs (from first embedding source)\ntest_protein_ids = list(embedding_dicts[list(embedding_dicts.keys())[0]]['test'].keys())\nprint(f\"Total test proteins: {len(test_protein_ids):,}\")\n\n# [Step 1/5] Generate DL predictions with multi-embedding\ndl_predictions = {}\n\nprint(\"\\n[Step 1/5] Generating multi-embedding DL predictions...\")\nwith torch.no_grad():\n    for start in tqdm(range(0, len(test_protein_ids), CONFIG['BATCH_SIZE']), \n                     desc=\"DL Inference\"):\n        batch_ids = test_protein_ids[start:start + CONFIG['BATCH_SIZE']]\n        \n        # Gather embeddings (match training logic)\n        if len(embedding_dicts) == 1:\n            src = list(embedding_dicts.keys())[0]\n            embeddings = torch.FloatTensor([\n                embedding_dicts[src]['test'][p] for p in batch_ids\n            ]).to(device)\n        else:\n            embeddings = []\n            for src in embedding_dicts.keys():\n                emb_batch = torch.FloatTensor([\n                    embedding_dicts[src]['test'][p] for p in batch_ids\n                ]).to(device)\n                embeddings.append(emb_batch)\n        \n        logits = model(embeddings, cooccur_tensor, return_logits=True)\n        outputs = torch.sigmoid(logits / CONFIG['TEMPERATURE']).cpu().numpy()\n        \n        for i, pid in enumerate(batch_ids):\n            dl_predictions[pid] = outputs[i]\n        \n        del embeddings, outputs, logits\n        if start % 5000 == 0:\n            gc.collect()\n\nprint(f\"✓ Generated DL predictions for {len(dl_predictions):,} proteins\")\n\n# [Step 2/5] Load BLAST predictions\nblast_dict = {}\n\nprint(\"\\n[Step 2/5] Loading BLAST predictions...\")\nif os.path.exists(CONFIG['BLAST_PATH']):\n    blast_data = defaultdict(lambda: np.zeros(len(top_terms)))\n    \n    for chunk in pd.read_csv(CONFIG['BLAST_PATH'], sep='\\t', header=None,\n                             names=['Id', 'GO term', 'Confidence'], \n                             chunksize=100000):\n        for _, row in chunk.iterrows():\n            pid = row['Id']\n            term = row['GO term']\n            conf = float(row['Confidence'])\n            \n            if term in term_to_idx:\n                idx = term_to_idx[term]\n                blast_data[pid][idx] = max(blast_data[pid][idx], conf)\n    \n    blast_dict = dict(blast_data)\n    print(f\"  ✓ Loaded BLAST for {len(blast_dict):,} proteins\")\nelse:\n    print(\"  ⚠️ BLAST file not found\")\n\n# [Step 3/5] Compute frequency baseline\nprint(\"\\n[Step 3/5] Computing frequency baseline...\")\nterm_frequencies = np.zeros(len(top_terms))\ntotal_proteins = len(protein_to_terms)\n\nfor idx, term in enumerate(top_terms):\n    count = term_freq.get(term, 0)\n    term_frequencies[idx] = min(count / total_proteins, 0.5)\n\nprint(f\"✓ Frequency baseline computed\")\n\n# [Step 4/5] Adaptive ensemble merging\nprint(\"\\n[Step 4/5] Adaptive ensemble merging...\")\n\ndef get_knowledge_level(protein_id, train_terms_dict, term_to_aspect_dict):\n    \"\"\"Detect knowledge level\"\"\"\n    if protein_id not in train_terms_dict:\n        return 'no'\n    \n    terms = train_terms_dict[protein_id]\n    aspects = set(term_to_aspect_dict.get(t, '') for t in terms)\n    aspect_map = {'F': 'MFO', 'P': 'BPO', 'C': 'CCO'}\n    mapped_aspects = set(aspect_map.get(a, a) for a in aspects)\n    \n    if len(mapped_aspects) == 0:\n        return 'no'\n    elif len(mapped_aspects) < 3:\n        return 'limited'\n    else:\n        return 'partial'\n\nmerged_predictions = {}\nknowledge_stats = Counter()\n\nfor pid, dl_probs in tqdm(dl_predictions.items(), desc=\"Merging\"):\n    knowledge = get_knowledge_level(pid, protein_to_terms, term_to_aspect)\n    knowledge_stats[knowledge] += 1\n    \n    blast_probs = blast_dict.get(pid, np.zeros(len(top_terms)))\n    \n    # Adaptive weighting\n    if knowledge == 'no':\n        dl_w = 0.35\n        freq_w = 0.15\n        blast_w = 0.50\n    elif knowledge == 'limited':\n        dl_w = 0.40\n        freq_w = 0.10\n        blast_w = 0.50\n    else:\n        dl_w = 0.50\n        freq_w = 0.05\n        blast_w = 0.45\n    \n    # Weighted combination\n    merged = (dl_probs * dl_w + \n              term_frequencies * freq_w + \n              blast_probs * blast_w)\n    \n    merged_predictions[pid] = merged\n\nprint(f\"✓ Merged predictions for {len(merged_predictions):,} proteins\")\nprint(\"\\n  Knowledge distribution:\")\nfor level, count in sorted(knowledge_stats.items()):\n    pct = count / len(merged_predictions) * 100\n    print(f\"    {level:8s}: {count:6,} ({pct:5.1f}%)\")\n\n# [Step 5/5] GO propagation + Ontology calibration\nprint(\"\\n[Step 5/5] GO propagation + Ontology calibration...\")\n\nfinal_predictions = {}\n\nfor pid, probs in tqdm(merged_predictions.items(), desc=\"Processing\"):\n    new_probs = probs.copy()\n    \n    # GO propagation\n    # Thêm threshold cao hơn cho propagation\n    PROPAGATION_MIN_CONF = 0.1 # Chỉ propagate high-confidence predictions # từ 0.3\n    if CONFIG['USE_GO_PROPAGATION'] and go_graph and ancestor_cache:\n        # ✅ Optimized\n        high_conf_indices = np.where(probs > PROPAGATION_MIN_CONF)[0]\n        for term_idx in high_conf_indices:\n            term = idx_to_term[term_idx]\n            score = probs[term_idx]\n            \n            if term in ancestor_cache:\n                ancestor_indices = [term_to_idx[a] for a in ancestor_cache[term] \n                                   if a in term_to_idx]\n                if ancestor_indices:\n                    propagated = score * CONFIG['PROPAGATION_DECAY']\n                    new_probs[ancestor_indices] = np.maximum(\n                        new_probs[ancestor_indices], \n                        propagated\n                    )\n    \n    # Ontology-specific calibration\n    for term_idx, score in enumerate(new_probs):\n        if score > CONFIG['MIN_CONFIDENCE']:\n            term = idx_to_term[term_idx]\n            aspect = term_to_aspect.get(term, 'F')\n            \n            aspect_map = {'F': 'MFO', 'P': 'BPO', 'C': 'CCO'}\n            ontology = aspect_map.get(aspect, 'MFO')\n            \n            calibration_factor = CONFIG['ONTOLOGY_CALIBRATION'][ontology]\n            new_probs[term_idx] = score * calibration_factor\n    \n    final_predictions[pid] = new_probs\n\nprint(f\"✓ Final predictions ready for {len(final_predictions):,} proteins\")\n\n\n\n\n# ============================================================================\n# ADAPTIVE CUTOFF HELPERS - Thêm vào trước phần WRITE SUBMISSION\n# ============================================================================\n\ndef get_adaptive_thresholds(protein_id, knowledge_level, blast_coverage, probs):\n    \"\"\"Tính adaptive thresholds dựa trên protein characteristics\"\"\"\n    \n    # Base thresholds\n    if knowledge_level == 'no':\n        if blast_coverage > 50:\n            min_conf = 0.008\n            max_preds = 3000\n        else:\n            min_conf = 0.015\n            max_preds = 1500\n    elif knowledge_level == 'limited':\n        min_conf = 0.010\n        max_preds = 2500\n    else:  # partial\n        min_conf = 0.008\n        max_preds = 3000\n    \n    # Điều chỉnh theo high-confidence predictions\n    high_conf_count = (probs > 0.1).sum()\n    \n    if high_conf_count > 100:\n        max_preds = min(max_preds + 500, 4000)\n    elif high_conf_count < 20:\n        max_preds = min(max_preds, 2000)\n        min_conf = max(min_conf, 0.015)\n    \n    return min_conf, max_preds\n\ndef should_include_prediction(score, position, prev_score, min_conf, max_preds):\n    \"\"\"Quyết định có include prediction không\"\"\"\n    \n    # Hard limits\n    if score < min_conf or position >= max_preds:\n        return False\n    \n    # Progressive thresholds\n    if position > 2000 and score < 0.02:\n        return False\n    \n    if position > 3000 and score < 0.05:\n        return False\n    \n    # Score drop detection (cliff detection)\n    if prev_score > 0 and position > 500:\n        if score < prev_score * 0.3:  # Drop > 70%\n            return False\n    \n    return True\n\n# ============================================================================\n# WRITE SUBMISSION WITH ADAPTIVE CUTOFF\n# ============================================================================\nprint(\"\\n\" + \"=\"*80)\nprint(\"WRITING SUBMISSION WITH ADAPTIVE CUTOFF\")\nprint(\"=\"*80)\n\nn_predictions = 0\nprotein_pred_counts = []\ncutoff_stats = {'min_conf': 0, 'max_preds': 0, 'quality': 0}\n\nwith open('submission.tsv', 'w') as f:\n    for pid, probs in tqdm(final_predictions.items(), desc=\"Writing\"):\n        # Get protein characteristics\n        knowledge = get_knowledge_level(pid, protein_to_terms, term_to_aspect)\n        blast_probs = blast_dict.get(pid, np.zeros(len(top_terms)))\n        blast_coverage = (blast_probs > 0.01).sum()\n        \n        # Get adaptive thresholds\n        min_conf, max_preds = get_adaptive_thresholds(\n            pid, knowledge, blast_coverage, probs\n        )\n        \n        # Sort predictions\n        top_indices = np.argsort(probs)[::-1]\n        \n        protein_preds = 0\n        prev_score = 1.0\n        \n        for position, idx in enumerate(top_indices):\n            score = probs[idx]\n            \n            # Check if should include\n            if should_include_prediction(score, position, prev_score, \n                                        min_conf, max_preds):\n                term = idx_to_term[idx]\n                f.write(f\"{pid}\\t{term}\\t{min(score, 0.999):.3f}\\n\")\n                n_predictions += 1\n                protein_preds += 1\n                prev_score = score\n            else:\n                # Track why stopped\n                if score < min_conf:\n                    cutoff_stats['min_conf'] += 1\n                elif position >= max_preds:\n                    cutoff_stats['max_preds'] += 1\n                else:\n                    cutoff_stats['quality'] += 1\n                break\n        \n        if protein_preds > 0:\n            protein_pred_counts.append(protein_preds)\n\nprint(f\"\\n✅ Adaptive cutoff statistics:\")\nprint(f\"  Stopped by min_confidence: {cutoff_stats['min_conf']:,}\")\nprint(f\"  Stopped by max_preds: {cutoff_stats['max_preds']:,}\")\nprint(f\"  Stopped by quality filter: {cutoff_stats['quality']:,}\")\n\n\nprint(\"\\n✅ ALL DONE!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-01T09:29:08.343889Z","iopub.execute_input":"2025-12-01T09:29:08.344473Z","execution_failed":"2025-12-01T09:30:28.302Z"}},"outputs":[],"execution_count":null}]}