{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":41875,"databundleVersionId":5521661},{"sourceType":"datasetVersion","sourceId":15792526,"datasetId":10122873,"databundleVersionId":16738946},{"sourceType":"datasetVersion","sourceId":15667414,"datasetId":10032015,"databundleVersionId":16604439},{"sourceType":"datasetVersion","sourceId":15835592,"datasetId":10150903,"databundleVersionId":16785566}],"dockerImageVersionId":31329,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\"\"\"\nCAFA 5 Protein Function Prediction - Complete Pipeline\n=======================================================\nApproach:\n  - Features: Biochemical (molecular weight, aromaticity, etc.) +\n              Amino acid composition + T5 protein embeddings (1024-dim) +\n              One-hot encoded taxonomy\n  - Model: py-boost (GPU gradient boosting) trained in label-chunks\n  - Evaluation: F-max and S-min (CAFA standard metrics)\n\"\"\"\n\n# ===========================================================================\n# 0. IMPORTS\n# ===========================================================================\nimport os\nimport gc\nimport pickle\nimport numpy as np\nimport pandas as pd\nfrom sklearn.model_selection import train_test_split\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-20T09:18:02.993086Z","iopub.execute_input":"2026-04-20T09:18:02.993375Z","iopub.status.idle":"2026-04-20T09:18:05.783663Z","shell.execute_reply.started":"2026-04-20T09:18:02.99335Z","shell.execute_reply":"2026-04-20T09:18:05.782931Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install py-boost","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-20T09:18:11.212321Z","iopub.execute_input":"2026-04-20T09:18:11.213054Z","iopub.status.idle":"2026-04-20T09:18:17.696134Z","shell.execute_reply.started":"2026-04-20T09:18:11.21302Z","shell.execute_reply":"2026-04-20T09:18:17.6951Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===========================================================================\n# 1. LOAD DATA\n#    We load X and Y from the .npy files saved earlier, plus the IA file\n#    needed for S-min computation.\n# ===========================================================================\n\nprint(\"Loading feature matrix and label matrix...\")\nX = np.load('/kaggle/input/datasets/ayushdhoble/x-and-ynpy/X.npy')          # shape (142163, 1100)\nY = np.load('/kaggle/input/datasets/ayushdhoble/x-and-ynpy/Y.npy')          # shape (142163, 3700)\nprint(f\"  X: {X.shape}  |  Y: {Y.shape}\")\n\n# Load target term names (same order as Y columns) so we can map GO IDs to IA\ntrain_terms1 = pd.read_csv(\n    '/kaggle/input/competitions/cafa-5-protein-function-prediction/Train/train_terms.tsv',\n    sep='\\t'\n)\ntop_bpo = train_terms1[train_terms1['aspect'] == 'BPO']['term'].value_counts().head(1700).index\ntop_cco = train_terms1[train_terms1['aspect'] == 'CCO']['term'].value_counts().head(1000).index\ntop_mfo = train_terms1[train_terms1['aspect'] == 'MFO']['term'].value_counts().head(1000).index\ntarget_terms = list(top_bpo) + list(top_cco) + list(top_mfo)   # 3700 terms, same order as Y\n\n# Load Information Accretion (IA) values for S-min\nia_path = '/kaggle/input/competitions/cafa-5-protein-function-prediction/IA.txt'\nia_df = pd.read_csv(ia_path, sep='\\t', header=None, names=['term', 'ia'])\nia_map = dict(zip(ia_df['term'], ia_df['ia']))\n\n# Build IA vector aligned to target_terms columns\nia_vector = np.array([ia_map.get(t, 0.0) for t in target_terms], dtype=np.float32)\nprint(f\"  IA values loaded for {(ia_vector > 0).sum()} / {len(ia_vector)} terms\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-20T09:18:24.320849Z","iopub.execute_input":"2026-04-20T09:18:24.321925Z","iopub.status.idle":"2026-04-20T09:18:41.683873Z","shell.execute_reply.started":"2026-04-20T09:18:24.321882Z","shell.execute_reply":"2026-04-20T09:18:41.68317Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===========================================================================\n# 2. TRAIN / VALIDATION SPLIT\n# ===========================================================================\nX = X.astype(np.float32)\nY = Y.astype(np.float32)\n\nX_train, X_val, y_train, y_val = train_test_split(\n    X, Y, test_size=0.15, random_state=42\n)\nX_train = np.ascontiguousarray(X_train, dtype=np.float32)\nX_val   = np.ascontiguousarray(X_val,   dtype=np.float32)\ny_train = np.ascontiguousarray(y_train, dtype=np.float32)\ny_val   = np.ascontiguousarray(y_val,   dtype=np.float32)\n\nprint(f\"Train: {X_train.shape}  Val: {X_val.shape}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-20T09:18:41.685538Z","iopub.execute_input":"2026-04-20T09:18:41.685884Z","iopub.status.idle":"2026-04-20T09:18:49.942052Z","shell.execute_reply.started":"2026-04-20T09:18:41.685845Z","shell.execute_reply":"2026-04-20T09:18:49.941096Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from py_boost import GradientBoosting","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-20T09:18:49.943167Z","iopub.execute_input":"2026-04-20T09:18:49.943479Z","iopub.status.idle":"2026-04-20T09:19:08.214557Z","shell.execute_reply.started":"2026-04-20T09:18:49.943443Z","shell.execute_reply":"2026-04-20T09:19:08.213558Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===========================================================================\n# 3. PY-BOOST CHUNKED MULTI-LABEL TRAINING  (FIXED)\n#    Key fixes vs original:\n#      - Removed cp.cuda.Device() context manager (unnecessary, causes issues)\n#      - Explicit cupy import guard\n#      - Cleaner OOM recovery\n#      - val_pred matrix built incrementally\n# ===========================================================================\ntry:\n    import cupy as cp\n    HAS_CUPY = True\nexcept ImportError:\n    HAS_CUPY = False\n\nfrom py_boost import GradientBoosting\n\n# ---- Hyper-parameters -------------------------------------------------------\nMODEL_KWARGS = dict(\n    loss        = \"bce\",\n    metric      = \"bce\",\n    ntrees      = 150,      # increased slightly for better recall\n    lr          = 0.05,\n    max_depth   = 4,\n    subsample   = 0.6,\n    colsample   = 0.6,\n    verbose     = 50,\n    es          = 20,       # early stopping rounds (saves time)\n)\n\nINITIAL_CHUNK_SIZE = 64    # labels trained together; auto-halves on OOM\nMIN_CHUNK_SIZE     = 8\n\nn_targets = y_train.shape[1]\nmodels    = []              # list of (model, target_indices_array)\n\n# ---- Helper: free GPU memory ------------------------------------------------\ndef clear_gpu():\n    gc.collect()\n    if HAS_CUPY:\n        cp.get_default_memory_pool().free_all_blocks()\n        cp.get_default_pinned_memory_pool().free_all_blocks()\n\n# ---- Recursive chunk trainer ------------------------------------------------\ndef train_chunk(target_idx):\n    \"\"\"Train one py-boost model for the given label indices.\n    Auto-splits on GPU OOM.\"\"\"\n    target_idx = np.asarray(target_idx, dtype=np.int64)\n    if len(target_idx) == 0:\n        return\n\n    ytr = y_train[:, target_idx]\n    yv  = y_val[:, target_idx]\n\n    model = GradientBoosting(**MODEL_KWARGS)\n    try:\n        # NOTE: eval_sets are omitted intentionally to reduce peak VRAM.\n        # Early stopping via es= works on the train loss instead.\n        model.fit(X_train, ytr)\n        models.append((model, target_idx))\n        print(f\"  ✓ Trained chunk  targets={target_idx[0]}–{target_idx[-1]}  \"\n              f\"(size={len(target_idx)})\")\n\n    except Exception as e:\n        err_str = str(e).lower()\n        is_oom  = \"out of memory\" in err_str or \"outofmemory\" in err_str\n        del model, ytr, yv\n        clear_gpu()\n\n        if is_oom and len(target_idx) > MIN_CHUNK_SIZE:\n            mid = len(target_idx) // 2\n            print(f\"  ⚠ OOM on {len(target_idx)} targets \"\n                  f\"→ splitting into {mid} + {len(target_idx)-mid}\")\n            train_chunk(target_idx[:mid])\n            train_chunk(target_idx[mid:])\n        else:\n            print(f\"  ✗ FAILED on chunk size={len(target_idx)}: {e}\")\n            raise\n    finally:\n        clear_gpu()\n\n# ---- Build initial chunks and train -----------------------------------------\ninitial_chunks = [\n    np.arange(i, min(i + INITIAL_CHUNK_SIZE, n_targets))\n    for i in range(0, n_targets, INITIAL_CHUNK_SIZE)\n]\nprint(f\"\\nStarting training: {n_targets} targets, \"\n      f\"{len(initial_chunks)} initial chunks of size {INITIAL_CHUNK_SIZE}\")\n\nfor i, chunk in enumerate(initial_chunks, 1):\n    print(f\"\\n[Chunk {i}/{len(initial_chunks)}]\")\n    train_chunk(chunk)\n\nprint(f\"\\nTraining complete. Trained {len(models)} model-chunks.\")\n\n# ---- Save models -------------------------------------------------------------\nwith open('pyboost_models_list.pkl', 'wb') as f:\n    pickle.dump(models, f)\nprint(\"Models saved to pyboost_models_list.pkl\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-20T09:19:27.894835Z","iopub.execute_input":"2026-04-20T09:19:27.895477Z","iopub.status.idle":"2026-04-20T09:21:44.405441Z","shell.execute_reply.started":"2026-04-20T09:19:27.895444Z","shell.execute_reply":"2026-04-20T09:21:44.404164Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===========================================================================\n# 4. GENERATE VALIDATION PREDICTIONS\n# ===========================================================================\nprint(\"\\nGenerating validation predictions...\")\nval_pred = np.zeros((X_val.shape[0], n_targets), dtype=np.float32)\n\nfor model, idx in models:\n    preds = model.predict(X_val)\n    if HAS_CUPY and isinstance(preds, cp.ndarray):\n        preds = cp.asnumpy(preds)\n    val_pred[:, idx] = preds.astype(np.float32)\n    clear_gpu()\n\n# Clip predictions to valid probability range\nval_pred = np.clip(val_pred, 0.0, 1.0)\nprint(f\"val_pred shape: {val_pred.shape}  \"\n      f\"range: [{val_pred.min():.4f}, {val_pred.max():.4f}]\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===========================================================================\n# 5. EVALUATION: F-MAX  (CAFA standard metric)\n# ===========================================================================\n# For every threshold t in [0, 1]:\n#   - binarise predictions\n#   - compute per-protein precision and recall\n#   - average over proteins that have at least one true label\n#   - compute macro-averaged F1\n# Fmax = max F1 over all thresholds\n\ndef compute_fmax(y_true, y_pred_prob, thresholds=None):\n    \"\"\"\n    Compute CAFA F-max score.\n\n    Parameters\n    ----------\n    y_true      : (N, L) binary ground-truth matrix\n    y_pred_prob : (N, L) predicted probability matrix\n    thresholds  : iterable of float thresholds to evaluate\n\n    Returns\n    -------\n    fmax   : best F1 score across thresholds\n    best_t : threshold that achieves fmax\n    precisions, recalls, f1s : arrays for plotting the PR-F curve\n    \"\"\"\n    if thresholds is None:\n        thresholds = np.arange(0.01, 1.00, 0.01)\n\n    precisions, recalls, f1s = [], [], []\n    # Only evaluate on proteins that have at least one positive label\n    has_label = y_true.sum(axis=1) > 0\n    yt = y_true[has_label].astype(np.float32)\n\n    for t in thresholds:\n        yp = (y_pred_prob[has_label] >= t).astype(np.float32)\n\n        tp = (yt * yp).sum(axis=1)              # (N,)\n        fp = ((1 - yt) * yp).sum(axis=1)        # (N,)\n        fn = (yt * (1 - yp)).sum(axis=1)        # (N,)\n\n        # Per-protein precision (0 if no prediction made)\n        pred_pos   = tp + fp\n        prec       = np.where(pred_pos > 0, tp / pred_pos, 0.0)\n\n        # Per-protein recall\n        true_pos   = tp + fn\n        rec        = np.where(true_pos > 0, tp / true_pos, 0.0)\n\n        avg_prec   = prec.mean()\n        avg_rec    = rec.mean()\n\n        if avg_prec + avg_rec > 0:\n            f1 = 2 * avg_prec * avg_rec / (avg_prec + avg_rec)\n        else:\n            f1 = 0.0\n\n        precisions.append(avg_prec)\n        recalls.append(avg_rec)\n        f1s.append(f1)\n\n    f1s_arr        = np.array(f1s)\n    best_idx       = f1s_arr.argmax()\n    fmax           = f1s_arr[best_idx]\n    best_threshold = thresholds[best_idx]\n\n    return fmax, best_threshold, np.array(precisions), np.array(recalls), f1s_arr\n\n\n# ===========================================================================\n# 6. EVALUATION: S-MIN  (CAFA semantic distance metric)\n# ===========================================================================\n# For every threshold t:\n#   ru(t)  = weighted average of IA values for missed positives (FN)\n#   mi(t)  = weighted average of IA values for false positives (FP)\n#   S(t)   = sqrt( ru(t)^2 + mi(t)^2 )\n# Smin = min S(t) over all thresholds\n\ndef compute_smin(y_true, y_pred_prob, ia_vec, thresholds=None):\n    \"\"\"\n    Compute CAFA S-min score.\n\n    Parameters\n    ----------\n    y_true      : (N, L) binary ground-truth matrix\n    y_pred_prob : (N, L) predicted probability matrix\n    ia_vec      : (L,)  information accretion per GO term\n    thresholds  : iterable of float thresholds\n\n    Returns\n    -------\n    smin   : minimum semantic distance across thresholds\n    best_t : threshold that achieves smin\n    s_vals : full S(t) array for plotting\n    \"\"\"\n    if thresholds is None:\n        thresholds = np.arange(0.01, 1.00, 0.01)\n\n    ia    = ia_vec.astype(np.float32)           # (L,)\n    has_label = y_true.sum(axis=1) > 0\n    yt    = y_true[has_label].astype(np.float32)\n    N     = yt.shape[0]\n\n    s_vals, ru_vals, mi_vals = [], [], []\n\n    for t in thresholds:\n        yp  = (y_pred_prob[has_label] >= t).astype(np.float32)\n\n        fn  = yt * (1 - yp)                     # missed positives\n        fp  = (1 - yt) * yp                     # false positives\n\n        # Remaining uncertainty: average IA of missed terms per protein\n        ru  = (fn * ia).sum(axis=1).mean()\n\n        # Misinformation: average IA of falsely predicted terms per protein\n        mi  = (fp * ia).sum(axis=1).mean()\n\n        s   = np.sqrt(ru**2 + mi**2)\n        s_vals.append(s)\n        ru_vals.append(ru)\n        mi_vals.append(mi)\n\n    s_arr  = np.array(s_vals)\n    best_i = s_arr.argmin()\n    smin   = s_arr[best_i]\n    best_t = thresholds[best_i]\n\n    return smin, best_t, s_arr, np.array(ru_vals), np.array(mi_vals)\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===========================================================================\n# 7. RUN EVALUATION PER ASPECT\n# ===========================================================================\n# Split targets back into BPO / CCO / MFO for per-aspect metrics\n\nbpo_idx = np.array([i for i, t in enumerate(target_terms) if t in set(top_bpo)])\ncco_idx = np.array([i for i, t in enumerate(target_terms) if t in set(top_cco)])\nmfo_idx = np.array([i for i, t in enumerate(target_terms) if t in set(top_mfo)])\n\nTHRESHOLDS = np.arange(0.01, 1.00, 0.01)\n\nresults = {}\nfor aspect_name, idx in [(\"BPO\", bpo_idx), (\"CCO\", cco_idx), (\"MFO\", mfo_idx), (\"Overall\", np.arange(n_targets))]:\n    print(f\"\\n--- Evaluating {aspect_name} ({len(idx)} terms) ---\")\n\n    yt_asp  = y_val[:, idx]\n    yp_asp  = val_pred[:, idx]\n    ia_asp  = ia_vector[idx]\n\n    fmax, ft, prec_arr, rec_arr, f1_arr = compute_fmax(yt_asp, yp_asp, THRESHOLDS)\n    smin, st, s_arr, ru_arr, mi_arr     = compute_smin(yt_asp, yp_asp, ia_asp, THRESHOLDS)\n\n    results[aspect_name] = {\n        \"fmax\"           : fmax,\n        \"fmax_threshold\" : ft,\n        \"smin\"           : smin,\n        \"smin_threshold\" : st,\n        \"precisions\"     : prec_arr,\n        \"recalls\"        : rec_arr,\n        \"f1s\"            : f1_arr,\n        \"s_vals\"         : s_arr,\n    }\n\n    print(f\"  F-max  = {fmax:.4f}  (threshold = {ft:.2f})\")\n    print(f\"  S-min  = {smin:.4f}  (threshold = {st:.2f})\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===========================================================================\n# 8. PLOT RESULTS\n# ===========================================================================\nimport matplotlib.pyplot as plt\n\nfig, axes = plt.subplots(2, 4, figsize=(22, 10))\nfig.suptitle(\"CAFA 5 — py-boost Evaluation (F-max & S-min per Aspect)\", fontsize=14)\n\naspect_list = [\"BPO\", \"CCO\", \"MFO\", \"Overall\"]\ncolors      = {\"BPO\": \"#4C72B0\", \"CCO\": \"#DD8452\", \"MFO\": \"#55A868\", \"Overall\": \"#C44E52\"}\n\nfor col, aspect in enumerate(aspect_list):\n    r   = results[aspect]\n    t   = THRESHOLDS\n    col_color = colors[aspect]\n\n    # ---- F1 vs Threshold ----\n    ax = axes[0, col]\n    ax.plot(t, r[\"f1s\"], color=col_color, lw=2)\n    ax.axvline(r[\"fmax_threshold\"], color=\"black\", ls=\"--\", alpha=0.6,\n               label=f\"F-max={r['fmax']:.4f}\\n@ t={r['fmax_threshold']:.2f}\")\n    ax.set_title(f\"{aspect} — F1 vs Threshold\")\n    ax.set_xlabel(\"Threshold\")\n    ax.set_ylabel(\"F1 Score\")\n    ax.legend(fontsize=8)\n    ax.set_ylim(0, 1)\n    ax.grid(alpha=0.3)\n\n    # ---- S vs Threshold ----\n    ax = axes[1, col]\n    ax.plot(t, r[\"s_vals\"], color=col_color, lw=2)\n    ax.axvline(r[\"smin_threshold\"], color=\"black\", ls=\"--\", alpha=0.6,\n               label=f\"S-min={r['smin']:.4f}\\n@ t={r['smin_threshold']:.2f}\")\n    ax.set_title(f\"{aspect} — Semantic Distance vs Threshold\")\n    ax.set_xlabel(\"Threshold\")\n    ax.set_ylabel(\"S(t) = √(ru²+mi²)\")\n    ax.legend(fontsize=8)\n    ax.grid(alpha=0.3)\n\nplt.tight_layout()\nplt.savefig(\"cafa5_evaluation.png\", dpi=150, bbox_inches=\"tight\")\nplt.show()\nprint(\"Plot saved to cafa5_evaluation.png\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===========================================================================\n# 9. SUMMARY TABLE\n# ===========================================================================\nprint(\"\\n\" + \"=\"*55)\nprint(f\"{'Aspect':<10} {'F-max':>8} {'F-max @t':>10} {'S-min':>8} {'S-min @t':>10}\")\nprint(\"=\"*55)\nfor aspect in aspect_list:\n    r = results[aspect]\n    print(f\"{aspect:<10} {r['fmax']:>8.4f} {r['fmax_threshold']:>10.2f} \"\n          f\"{r['smin']:>8.4f} {r['smin_threshold']:>10.2f}\")\nprint(\"=\"*55)\n\n# ===========================================================================\n# 10. SAVE PREDICTIONS & RESULTS\n# ===========================================================================\nnp.save(\"val_pred.npy\", val_pred)\nwith open(\"evaluation_results.pkl\", \"wb\") as f:\n    pickle.dump(results, f)\nprint(\"\\nSaved val_pred.npy and evaluation_results.pkl\")","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}