{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":null,"isInternetEnabled":false,"language":"python","sourceType":"notebook"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"761041c5","cell_type":"markdown","source":"# 68th Place Solution for the CAFA 6 Protein Function Prediction Competition\n\n**Silver Medal — 68th of 2,259 teams**\n\nThis notebook presents the prediction-stage ensemble used for my final CAFA 6 submission. It combines two public prediction pipelines and focuses on the part I developed: **score harmonization, agreement-aware stacking, Gene Ontology processing, Information Accretion-aware ranking, memory-efficient inference, and submission validation**.\n\n### Context\n\n- [Competition overview](https://www.kaggle.com/competitions/cafa-6-protein-function-prediction/overview)\n- [Competition data](https://www.kaggle.com/competitions/cafa-6-protein-function-prediction/data)\n\n### Public inputs\n\n| Source | Public notebook | Contribution used here |\n|---|---|---|\n| **Model A** | [CAFA ensemble 0.378](https://www.kaggle.com/code/antonoof/cafa-ensemble-0-378) by **antonoof** | Fold-level ensemble with GO ancestor propagation |\n| **Model B** | [CAFA6 GOA + ProtT5 IC-Aware](https://www.kaggle.com/code/maryamkhaled67/cafa6-goa-prott5-ic-aware) by **maryamkhaled67** | GOA transfer, ProtT5/InterPro evidence, GO propagation, and IA-aware filtering |\n\nThis is a **stacking and post-processing solution**, not a claim that I trained the upstream models.\n","metadata":{}},{"id":"eeb35f91","cell_type":"markdown","source":"## Overview of the approach\n\nThe two public submissions contain complementary evidence but cannot be averaged blindly:\n\n1. **Score harmonization.** Model A may export values outside the probability range, whereas Model B is probability-like. I inspect both files first and apply a temperature-scaled sigmoid only when a source is unbounded.\n2. **Memory-aware loading.** Each file contains tens of millions of rows, so predictions are read in chunks and pruned per protein.\n3. **Ontology consistency.** Alternative GO IDs are canonicalized and each prediction is propagated to all ancestors by hard maximum.\n4. **Agreement-aware stacking.** A third signal rewards terms supported by both sources. The three signals are combined in logit space.\n5. **IA-aware ranking.** Low-information generic terms are pruned unless their confidence is high; remaining terms are ranked using probability and Information Accretion.\n6. **Strict output checks.** Scores, identifiers, duplicate pairs, formatting, and the 1,500-term limit are validated before publication.\n\n```text\nModel A ─ score-space check ─ top-k ─┐\n                                      ├─ agreement-aware logit stack\nModel B ─ score-space check ─ top-k ─┘\n                    │\n                    ▼\n       GO canonicalization and propagation\n                    │\n                    ▼\n          IA-aware pruning and ranking\n                    │\n                    ▼\n         submission.tsv + validation\n```\n","metadata":{}},{"id":"278fb180","cell_type":"markdown","source":"## Validation caveat\n\nCAFA 6 is prospective: the final labels are annotations that appeared after the submission deadline. Current training annotations can therefore be used only as a **proxy consistency check**, not as a reliable estimate of the final leaderboard score.\n\nThe optional proxy evaluation below:\n\n- uses target proteins that also occur in `train_terms.tsv`;\n- adds GO ancestor closure to the proxy truth;\n- separates calibration and holdout proteins deterministically;\n- reports IA-weighted F-max by ontology.\n\nBecause one public source contains GOA-derived information, this proxy can favor that source. The final competition parameters are therefore recorded explicitly rather than presented as the result of a trustworthy hidden-label validation scheme.\n","metadata":{}},{"id":"7506e97e","cell_type":"code","source":"# Imports and reproducibility\nimport gc\nimport hashlib\nimport math\nimport os\nimport random\nfrom collections import defaultdict\nfrom functools import lru_cache\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm.auto import tqdm\n\nSEED = 7\nrandom.seed(SEED)\nRNG = np.random.default_rng(SEED)\n\nprint(\"NumPy:\", np.__version__)\nprint(\"pandas:\", pd.__version__)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:26:23.491701Z","iopub.execute_input":"2026-07-28T15:26:23.492083Z","iopub.status.idle":"2026-07-28T15:26:23.832318Z","shell.execute_reply.started":"2026-07-28T15:26:23.492053Z","shell.execute_reply":"2026-07-28T15:26:23.831507Z"}},"outputs":[],"execution_count":null},{"id":"a80c0832","cell_type":"code","source":"# Input resolution and configuration\ndef first_existing(label, candidates):\n    for candidate in candidates:\n        if Path(candidate).exists():\n            print(f\"{label}: {candidate}\")\n            return candidate\n    raise FileNotFoundError(\n        f\"Could not find {label}. Add the required Kaggle input and check one of:\\n\"\n        + \"\\n\".join(f\"  - {path}\" for path in candidates)\n    )\n\n\nOBO_PATH = first_existing(\n    \"GO ontology\",\n    [\n        \"/kaggle/input/cafa-6-protein-function-prediction/Train/go-basic.obo\",\n        \"/kaggle/input/competitions/cafa-6-protein-function-prediction/Train/go-basic.obo\",\n    ],\n)\nIA_PATH = first_existing(\n    \"Information Accretion\",\n    [\n        \"/kaggle/input/cafa-6-protein-function-prediction/IA.tsv\",\n        \"/kaggle/input/competitions/cafa-6-protein-function-prediction/IA.tsv\",\n    ],\n)\nTRAIN_TERMS_PATH = first_existing(\n    \"training annotations\",\n    [\n        \"/kaggle/input/cafa-6-protein-function-prediction/Train/train_terms.tsv\",\n        \"/kaggle/input/competitions/cafa-6-protein-function-prediction/Train/train_terms.tsv\",\n    ],\n)\nTEST_FASTA_PATH = first_existing(\n    \"test superset\",\n    [\n        \"/kaggle/input/cafa-6-protein-function-prediction/Test/testsuperset.fasta\",\n        \"/kaggle/input/competitions/cafa-6-protein-function-prediction/Test/testsuperset.fasta\",\n    ],\n)\n\nMODEL_A_PATH = first_existing(\n    \"Model A — CAFA ensemble 0.378\",\n    [\n        \"/kaggle/input/cafa-ensemble-0-378/submission.tsv\",\n        \"/kaggle/input/notebooks/antonoof/cafa-ensemble-0-378/submission.tsv\",\n    ],\n)\nMODEL_B_PATH = first_existing(\n    \"Model B — GOA + ProtT5\",\n    [\n        \"/kaggle/input/cafa6-goa-prott5-ic-aware/submission.tsv\",\n        \"/kaggle/input/notebooks/maryamkhaled67/cafa6-goa-prott5-ic-aware/submission.tsv\",\n    ],\n)\n\nOUTPUT_PATH = \"/kaggle/working/submission.tsv\"\n\n# Fail early if the same file was accidentally assigned to both sources.\ndef partial_signature(path, n_bytes=8 * 1024 * 1024):\n    digest = hashlib.sha256()\n    with open(path, \"rb\") as handle:\n        digest.update(handle.read(n_bytes))\n    return os.path.getsize(path), digest.hexdigest()\n\nif Path(MODEL_A_PATH).resolve() == Path(MODEL_B_PATH).resolve():\n    raise ValueError(\"Model A and Model B point to the same file.\")\n\nif partial_signature(MODEL_A_PATH) == partial_signature(MODEL_B_PATH):\n    raise ValueError(\n        \"Model A and Model B have the same file size and initial checksum. \"\n        \"Check the attached notebook outputs before continuing.\"\n    )\n\n# Memory controls\nCHUNK_SIZE = 2_000_000\nMIN_READ_SCORE = 1e-4\nTOPK_LOAD_A = 220\nTOPK_LOAD_B = 500\n\n# Score harmonization\nTEMP_A = 8.0\nTEMP_B = 1.0\nUNBOUNDED_THRESHOLD = 1.2\n\n# Conditional tie-breaking for unusually flat Model B profiles\nUSE_B_FLAT_IA_BOOST = True\nB_FLAT_STD_THRESHOLD = 0.01\nB_IA_BOOST_ALPHA = 0.12\n\n# Final prediction rules\nMIN_BLEND_SCORE = 0.005\nMIN_OUTPUT_SCORE = 0.001\nTOPK_FINAL = 650\nHARD_CAP = 1500\nMIN_IA = 1.0\nKEEP_LOW_IA_IF_SCORE_GE = 0.60\nIA_EPS = 0.10\nRANK_DECAY_TAU = 90.0\nSIGNIFICANT_DIGITS = 3\n\n# Proxy evaluation\nRUN_TUNING = False\nRUN_PROXY_ABLATION = True\nTUNING_TRIALS = 40\nPROXY_SAMPLE_N = 2500\nPROXY_HOLDOUT_FRACTION = 0.20\nPROXY_THRESHOLDS = np.linspace(0.05, 0.80, 25, dtype=np.float32)\n\n# Parameters recorded from the final competition notebook\nRECORDED_PARAMETERS = {\n    \"A\": 0.08729884028434753,\n    \"B\": 0.7224695682525635,\n    \"agreement\": 0.190231591463089,\n    \"temperature\": 0.8547518464920547,\n}\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:26:23.833562Z","iopub.execute_input":"2026-07-28T15:26:23.834050Z","iopub.status.idle":"2026-07-28T15:26:24.096565Z","shell.execute_reply.started":"2026-07-28T15:26:23.833980Z","shell.execute_reply":"2026-07-28T15:26:24.096023Z"}},"outputs":[],"execution_count":null},{"id":"03e57d90","cell_type":"markdown","source":"## Parse the Gene Ontology and official metadata\n\nThe parser keeps primary GO terms, namespaces, `is_a` and `part_of` parents, alternative identifiers, and valid replacements for obsolete terms. Ancestors are cached because the same closure is reused many times.\n","metadata":{}},{"id":"0ae1e0f7","cell_type":"code","source":"def parse_fasta_ids(path):\n    ids = []\n    with open(path, \"r\") as handle:\n        for line in handle:\n            if line.startswith(\">\"):\n                ids.append(line[1:].strip().split()[0])\n    return ids\n\n\ndef parse_go_obo(path):\n    raw_terms = {}\n    current = None\n\n    def new_record():\n        return {\n            \"id\": None,\n            \"namespace\": None,\n            \"parents\": [],\n            \"alt_ids\": [],\n            \"is_obsolete\": False,\n            \"replaced_by\": None,\n        }\n\n    with open(path, \"r\") as handle:\n        for raw_line in handle:\n            line = raw_line.strip()\n\n            if line.startswith(\"[\") and line.endswith(\"]\"):\n                if current and current[\"id\"]:\n                    raw_terms[current[\"id\"]] = current\n                current = new_record() if line == \"[Term]\" else None\n                continue\n\n            if current is None or not line:\n                continue\n\n            if line.startswith(\"id: \"):\n                current[\"id\"] = line.split(\"id: \", 1)[1]\n            elif line.startswith(\"namespace: \"):\n                current[\"namespace\"] = line.split(\"namespace: \", 1)[1]\n            elif line.startswith(\"alt_id: \"):\n                current[\"alt_ids\"].append(line.split(\"alt_id: \", 1)[1])\n            elif line.startswith(\"is_a: \"):\n                current[\"parents\"].append(line.split()[1])\n            elif line.startswith(\"relationship: part_of \"):\n                current[\"parents\"].append(line.split()[2])\n            elif line == \"is_obsolete: true\":\n                current[\"is_obsolete\"] = True\n            elif line.startswith(\"replaced_by: \"):\n                current[\"replaced_by\"] = line.split(\"replaced_by: \", 1)[1]\n\n    if current and current[\"id\"]:\n        raw_terms[current[\"id\"]] = current\n\n    valid_terms = {\n        term for term, record in raw_terms.items()\n        if not record[\"is_obsolete\"]\n    }\n\n    namespaces = {\n        term: record[\"namespace\"]\n        for term, record in raw_terms.items()\n        if term in valid_terms and record[\"namespace\"] is not None\n    }\n\n    parents = defaultdict(list)\n    for term, record in raw_terms.items():\n        if term in valid_terms:\n            parents[term] = [\n                parent for parent in record[\"parents\"]\n                if parent in valid_terms\n            ]\n\n    aliases = {}\n    for term, record in raw_terms.items():\n        if term in valid_terms:\n            for alias in record[\"alt_ids\"]:\n                aliases[alias] = term\n        elif record[\"replaced_by\"] in valid_terms:\n            aliases[term] = record[\"replaced_by\"]\n            for alias in record[\"alt_ids\"]:\n                aliases[alias] = record[\"replaced_by\"]\n\n    return parents, namespaces, aliases, valid_terms\n\n\ndef namespace_code(namespace):\n    return {\n        \"biological_process\": \"BPO\",\n        \"cellular_component\": \"CCO\",\n        \"molecular_function\": \"MFO\",\n    }.get(namespace)\n\n\ntarget_ids = parse_fasta_ids(TEST_FASTA_PATH)\ntarget_set = set(target_ids)\nparents, namespace_map, alias_to_primary, valid_terms = parse_go_obo(OBO_PATH)\n\n\ndef canonicalize(term):\n    return alias_to_primary.get(term, term)\n\n\n@lru_cache(maxsize=None)\ndef ancestors(term):\n    term = canonicalize(term)\n    seen = set()\n    stack = list(parents.get(term, ()))\n\n    while stack:\n        parent = canonicalize(stack.pop())\n        if parent in seen or parent not in valid_terms:\n            continue\n        seen.add(parent)\n        stack.extend(parents.get(parent, ()))\n\n    return tuple(seen)\n\n\nia_df = pd.read_csv(IA_PATH, sep=\"\\t\", header=None, names=[\"term\", \"ia\"])\nia_df[\"term\"] = ia_df[\"term\"].astype(str).map(canonicalize)\nia_df[\"ia\"] = pd.to_numeric(ia_df[\"ia\"], errors=\"coerce\").fillna(0.0).astype(np.float32)\nterm_ia = dict(zip(ia_df[\"term\"], ia_df[\"ia\"]))\n\nprint(f\"Targets: {len(target_ids):,}\")\nprint(f\"Valid GO terms: {len(valid_terms):,}\")\nprint(f\"GO edges: {sum(len(values) for values in parents.values()):,}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:26:24.097359Z","iopub.execute_input":"2026-07-28T15:26:24.097549Z","iopub.status.idle":"2026-07-28T15:26:27.151256Z","shell.execute_reply.started":"2026-07-28T15:26:24.097530Z","shell.execute_reply":"2026-07-28T15:26:27.150436Z"}},"outputs":[],"execution_count":null},{"id":"a8ef7cec","cell_type":"markdown","source":"## Inspect the source score spaces\n\nThis diagnostic is part of the method, not only an EDA display. A source with values above 1 cannot be treated as a probability file. The transform is selected from the observed range and displayed before loading the full submissions.\n","metadata":{}},{"id":"9b716d7b","cell_type":"code","source":"def inspect_score_file(path, nrows=1_000_000):\n    sample = pd.read_csv(\n        path,\n        sep=\"\\t\",\n        header=None,\n        names=[\"protein\", \"term\", \"score\"],\n        nrows=nrows,\n        usecols=[0, 1, 2],\n    )\n    scores = pd.to_numeric(sample[\"score\"], errors=\"coerce\").dropna().to_numpy(np.float64)\n    quantiles = np.quantile(scores, [0.00, 0.01, 0.50, 0.99, 1.00])\n    return pd.Series(\n        {\n            \"rows_sampled\": len(scores),\n            \"min\": quantiles[0],\n            \"p01\": quantiles[1],\n            \"median\": quantiles[2],\n            \"p99\": quantiles[3],\n            \"max\": quantiles[4],\n            \"file_size_MB\": os.path.getsize(path) / 1e6,\n        }\n    )\n\n\ndiagnostics = pd.DataFrame(\n    {\n        \"Model A\": inspect_score_file(MODEL_A_PATH),\n        \"Model B\": inspect_score_file(MODEL_B_PATH),\n    }\n).T\n\ndiagnostics[\"score_space\"] = np.where(\n    diagnostics[\"max\"] > UNBOUNDED_THRESHOLD,\n    \"unbounded → sigmoid\",\n    \"probability\",\n)\ndisplay(diagnostics.round(4))\n\nMODEL_A_SCORE_SPACE = (\n    \"unbounded_sigmoid\"\n    if diagnostics.loc[\"Model A\", \"max\"] > UNBOUNDED_THRESHOLD\n    else \"probability\"\n)\nMODEL_B_SCORE_SPACE = (\n    \"unbounded_sigmoid\"\n    if diagnostics.loc[\"Model B\", \"max\"] > UNBOUNDED_THRESHOLD\n    else \"probability\"\n)\n\nif MODEL_A_SCORE_SPACE == MODEL_B_SCORE_SPACE:\n    print(\n        \"Note: both files appear to use the same numeric range. \"\n        \"The distinct-file checks above have passed, so processing continues.\"\n    )\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:26:27.152916Z","iopub.execute_input":"2026-07-28T15:26:27.153236Z","iopub.status.idle":"2026-07-28T15:26:28.183059Z","shell.execute_reply.started":"2026-07-28T15:26:27.153212Z","shell.execute_reply":"2026-07-28T15:26:28.182351Z"}},"outputs":[],"execution_count":null},{"id":"e983639d","cell_type":"markdown","source":"## Load both submissions in chunks\n\nOnly valid test targets and GO terms are retained. Duplicate `(protein, term)` pairs are merged by maximum score. Periodic pruning controls peak memory, followed by an exact final top-k pass.\n","metadata":{}},{"id":"6d5bdcdb","cell_type":"code","source":"def sigmoid(values):\n    values = np.clip(values, -60.0, 60.0)\n    return 1.0 / (1.0 + np.exp(-values))\n\n\ndef safe_logit(probability):\n    probability = np.clip(probability, 1e-6, 1.0 - 1e-6)\n    return np.log(probability / (1.0 - probability))\n\n\ndef transform_scores(raw_scores, score_space, temperature):\n    raw_scores = np.asarray(raw_scores, dtype=np.float32)\n\n    if score_space == \"unbounded_sigmoid\":\n        probabilities = sigmoid(raw_scores / float(temperature))\n    elif score_space == \"probability\":\n        probabilities = raw_scores\n    else:\n        raise ValueError(f\"Unknown score space: {score_space}\")\n\n    return np.clip(probabilities, MIN_OUTPUT_SCORE, 1.0).astype(np.float32)\n\n\ndef prune_prediction_dict(values, topk, ia_rank=False):\n    if len(values) <= topk:\n        return values\n\n    if ia_rank:\n        ranked = sorted(\n            values.items(),\n            key=lambda item: item[1] * (float(term_ia.get(item[0], 0.0)) + IA_EPS),\n            reverse=True,\n        )[:topk]\n    else:\n        ranked = sorted(\n            values.items(),\n            key=lambda item: item[1],\n            reverse=True,\n        )[:topk]\n\n    return dict(ranked)\n\n\ndef load_submission(path, *, name, topk, score_space, temperature, ia_rank=False):\n    predictions = defaultdict(dict)\n    rows_read = 0\n    canonical_rows = 0\n\n    reader = pd.read_csv(\n        path,\n        sep=\"\\t\",\n        header=None,\n        names=[\"protein\", \"term\", \"score\"],\n        dtype={\"protein\": str, \"term\": str},\n        chunksize=CHUNK_SIZE,\n        on_bad_lines=\"skip\",\n    )\n\n    for chunk in tqdm(reader, desc=f\"Load {name}\"):\n        rows_read += len(chunk)\n\n        chunk = chunk[\n            chunk[\"protein\"].isin(target_set)\n            & chunk[\"term\"].str.startswith(\"GO:\", na=False)\n        ].copy()\n        if chunk.empty:\n            continue\n\n        chunk[\"score\"] = pd.to_numeric(chunk[\"score\"], errors=\"coerce\")\n        chunk = chunk.dropna(subset=[\"score\"])\n        if chunk.empty:\n            continue\n\n        chunk[\"prob\"] = transform_scores(\n            chunk[\"score\"].to_numpy(),\n            score_space=score_space,\n            temperature=temperature,\n        )\n        chunk = chunk[chunk[\"prob\"] >= MIN_READ_SCORE]\n        chunk[\"term\"] = chunk[\"term\"].map(canonicalize)\n        chunk = chunk[chunk[\"term\"].isin(valid_terms)]\n        if chunk.empty:\n            continue\n\n        chunk = chunk.groupby([\"protein\", \"term\"], as_index=False)[\"prob\"].max()\n        canonical_rows += len(chunk)\n        touched = set()\n\n        for protein, term, probability in chunk.itertuples(index=False, name=None):\n            previous = predictions[protein].get(term)\n            if previous is None or probability > previous:\n                predictions[protein][term] = float(probability)\n            touched.add(protein)\n\n        for protein in touched:\n            if len(predictions[protein]) > int(1.5 * topk):\n                predictions[protein] = prune_prediction_dict(\n                    predictions[protein],\n                    topk=topk,\n                    ia_rank=ia_rank,\n                )\n\n        del chunk\n        gc.collect()\n\n    for protein in tqdm(list(predictions), desc=f\"Final top-k {name}\"):\n        predictions[protein] = prune_prediction_dict(\n            predictions[protein],\n            topk=topk,\n            ia_rank=ia_rank,\n        )\n\n    for protein in target_ids:\n        predictions.setdefault(protein, {})\n\n    print(\n        f\"{name}: rows read={rows_read:,} | canonical rows={canonical_rows:,} \"\n        f\"| proteins={len(predictions):,}\"\n    )\n    return predictions\n\n\npred_a = load_submission(\n    MODEL_A_PATH,\n    name=\"Model A\",\n    topk=TOPK_LOAD_A,\n    score_space=MODEL_A_SCORE_SPACE,\n    temperature=TEMP_A,\n    ia_rank=False,\n)\n\npred_b = load_submission(\n    MODEL_B_PATH,\n    name=\"Model B\",\n    topk=TOPK_LOAD_B,\n    score_space=MODEL_B_SCORE_SPACE,\n    temperature=TEMP_B,\n    ia_rank=True,\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:26:28.183921Z","iopub.execute_input":"2026-07-28T15:26:28.184218Z","iopub.status.idle":"2026-07-28T15:29:54.213572Z","shell.execute_reply.started":"2026-07-28T15:26:28.184189Z","shell.execute_reply":"2026-07-28T15:29:54.212803Z"}},"outputs":[],"execution_count":null},{"id":"bf1036ea","cell_type":"code","source":"# Optional tie-breaking for unusually flat Model B profiles\ndef boost_flat_profiles_inplace(predictions):\n    changed = 0\n\n    for terms in tqdm(predictions.values(), desc=\"Flat-profile adjustment\"):\n        if len(terms) < 20:\n            continue\n\n        values = np.fromiter(terms.values(), dtype=np.float32)\n        if float(values.std()) >= B_FLAT_STD_THRESHOLD:\n            continue\n\n        for term in list(terms):\n            ia = float(term_ia.get(term, 0.0))\n            terms[term] = float(\n                np.clip(\n                    terms[term] * (1.0 + B_IA_BOOST_ALPHA * np.log1p(ia)),\n                    MIN_OUTPUT_SCORE,\n                    1.0,\n                )\n            )\n        changed += 1\n\n    print(f\"Adjusted profiles: {changed:,}\")\n\n\nif USE_B_FLAT_IA_BOOST:\n    boost_flat_profiles_inplace(pred_b)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:29:54.215224Z","iopub.execute_input":"2026-07-28T15:29:54.215943Z","iopub.status.idle":"2026-07-28T15:30:01.290237Z","shell.execute_reply.started":"2026-07-28T15:29:54.215916Z","shell.execute_reply":"2026-07-28T15:30:01.289172Z"}},"outputs":[],"execution_count":null},{"id":"7386fba0","cell_type":"markdown","source":"## Source complementarity\n\nThe stack is useful only when the two sources are not identical. The table below summarizes term-set overlap on a deterministic sample of proteins. A very high overlap is not automatically an error, but exact equality should prompt a check of the attached inputs.\n","metadata":{}},{"id":"baf2517c","cell_type":"code","source":"sample_size = min(5_000, len(target_ids))\nsample_proteins = RNG.choice(target_ids, size=sample_size, replace=False)\n\noverlap_rows = []\nfor protein in sample_proteins:\n    terms_a = set(pred_a[protein])\n    terms_b = set(pred_b[protein])\n    union = terms_a | terms_b\n    intersection = terms_a & terms_b\n\n    overlap_rows.append(\n        {\n            \"terms_A\": len(terms_a),\n            \"terms_B\": len(terms_b),\n            \"intersection\": len(intersection),\n            \"union\": len(union),\n            \"jaccard\": len(intersection) / len(union) if union else 0.0,\n        }\n    )\n\nsource_overlap = pd.DataFrame(overlap_rows)\ndisplay(source_overlap.agg([\"mean\", \"median\", \"max\"]).round(3))\n\nexact_match_fraction = np.mean(\n    [\n        set(pred_a[protein]) == set(pred_b[protein])\n        for protein in sample_proteins\n    ]\n)\nprint(f\"Exact term-set match fraction: {exact_match_fraction:.3f}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:30:01.291362Z","iopub.execute_input":"2026-07-28T15:30:01.291715Z","iopub.status.idle":"2026-07-28T15:30:01.683749Z","shell.execute_reply.started":"2026-07-28T15:30:01.291690Z","shell.execute_reply":"2026-07-28T15:30:01.682893Z"}},"outputs":[],"execution_count":null},{"id":"0847b433","cell_type":"markdown","source":"## Agreement-aware logit stacking\n\nFor each candidate term:\n\n- `pA` is the harmonized Model A score;\n- `pB` is the Model B probability;\n- `pAgreement` is a product-of-experts-style signal defined only when both sources predict the term.\n\nMissing evidence receives a very small probability. The blend therefore rewards consensus while still allowing a strong source-specific prediction to survive. After blending, scores are propagated to all GO ancestors, filtered using IA, and limited to the highest-ranked terms.\n","metadata":{}},{"id":"3331ff4d","cell_type":"code","source":"ASPECTS = (\"BPO\", \"CCO\", \"MFO\")\nMISSING_PROBABILITY = 1e-6\nAGREEMENT_TEMPERATURE = 1.0\n\n\ndef agreement_probability(prob_a, prob_b):\n    value = (\n        safe_logit(prob_a) + safe_logit(prob_b)\n    ) / AGREEMENT_TEMPERATURE\n    return float(sigmoid(value))\n\n\ndef blended_probability(prob_a, prob_b, prob_agreement, weights, temperature):\n    logit_value = (\n        weights[\"A\"] * safe_logit(prob_a)\n        + weights[\"B\"] * safe_logit(prob_b)\n        + weights[\"agreement\"] * safe_logit(prob_agreement)\n    ) / max(float(temperature), 1e-6)\n\n    return float(sigmoid(logit_value))\n\n\ndef predict_protein(protein, weights, blend_temperature):\n    source_a = pred_a[protein]\n    source_b = pred_b[protein]\n    candidate_terms = set(source_a).union(source_b)\n\n    if not candidate_terms:\n        return []\n\n    scores = {}\n    for term in candidate_terms:\n        prob_a = source_a.get(term, MISSING_PROBABILITY)\n        prob_b = source_b.get(term, MISSING_PROBABILITY)\n\n        if term in source_a and term in source_b:\n            prob_agreement = agreement_probability(prob_a, prob_b)\n        else:\n            prob_agreement = MISSING_PROBABILITY\n\n        score = blended_probability(\n            prob_a,\n            prob_b,\n            prob_agreement,\n            weights=weights,\n            temperature=blend_temperature,\n        )\n\n        if score >= MIN_BLEND_SCORE:\n            scores[term] = score\n\n    # Hard-max propagation to all ancestors\n    for term, score in sorted(scores.items(), key=lambda item: item[1], reverse=True):\n        for ancestor in ancestors(term):\n            if score > scores.get(ancestor, 0.0):\n                scores[ancestor] = score\n\n    ranked = []\n    for term, score in scores.items():\n        aspect = namespace_code(namespace_map.get(term))\n        if aspect is None:\n            continue\n\n        ia = float(term_ia.get(term, 0.0))\n        if ia < MIN_IA and score < KEEP_LOW_IA_IF_SCORE_GE:\n            continue\n\n        rank_score = score * (ia + IA_EPS)\n        ranked.append((term, score, rank_score, aspect))\n\n    ranked.sort(key=lambda item: item[2], reverse=True)\n    ranked = ranked[: min(TOPK_FINAL, HARD_CAP)]\n\n    output = []\n    for rank, (term, score, _, aspect) in enumerate(ranked):\n        shaped_score = score * math.exp(-rank / RANK_DECAY_TAU)\n        shaped_score = float(np.clip(shaped_score, MIN_OUTPUT_SCORE, 1.0))\n        output.append((term, shaped_score, aspect))\n\n    return output\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:30:01.684873Z","iopub.execute_input":"2026-07-28T15:30:01.685275Z","iopub.status.idle":"2026-07-28T15:30:01.696537Z","shell.execute_reply.started":"2026-07-28T15:30:01.685248Z","shell.execute_reply":"2026-07-28T15:30:01.695911Z"}},"outputs":[],"execution_count":null},{"id":"1654281d","cell_type":"markdown","source":"## Proxy evaluation and recorded parameters\n\nThe code below implements an IA-weighted proxy F-max over BP, CC, and MF. It is included to make the calibration logic inspectable, not to reconstruct the prospective leaderboard.\n\n`RUN_TUNING=False` uses the parameters from the final competition notebook. Setting it to `True` performs a deterministic search on the calibration subset.\n","metadata":{}},{"id":"543259ca","cell_type":"code","source":"def normalize_aspect(value):\n    return {\n        \"P\": \"BPO\", \"BP\": \"BPO\", \"BPO\": \"BPO\",\n        \"C\": \"CCO\", \"CC\": \"CCO\", \"CCO\": \"CCO\",\n        \"F\": \"MFO\", \"MF\": \"MFO\", \"MFO\": \"MFO\",\n    }.get(str(value))\n\n\ndef build_proxy_truth(path):\n    frame = pd.read_csv(path, sep=\"\\t\")\n\n    rename = {}\n    for column in frame.columns:\n        lower = column.lower()\n        if lower in {\"entryid\", \"entry_id\", \"protein_id\", \"id\"}:\n            rename[column] = \"protein\"\n        elif lower in {\"term\", \"go_term\", \"goterm\"}:\n            rename[column] = \"term\"\n        elif lower in {\"aspect\", \"ontology\"}:\n            rename[column] = \"aspect\"\n\n    frame = frame.rename(columns=rename)\n    required = {\"protein\", \"term\", \"aspect\"}\n    missing = required.difference(frame.columns)\n    if missing:\n        raise ValueError(f\"Missing columns in train_terms.tsv: {missing}\")\n\n    frame[\"protein\"] = frame[\"protein\"].astype(str)\n    frame[\"term\"] = frame[\"term\"].astype(str).map(canonicalize)\n    frame[\"aspect\"] = frame[\"aspect\"].map(normalize_aspect)\n\n    frame = frame[\n        frame[\"protein\"].isin(target_set)\n        & frame[\"term\"].isin(valid_terms)\n        & frame[\"aspect\"].isin(ASPECTS)\n    ]\n\n    proteins = sorted(frame[\"protein\"].unique())\n    if len(proteins) > PROXY_SAMPLE_N:\n        proteins = sorted(\n            RNG.choice(proteins, size=PROXY_SAMPLE_N, replace=False).tolist()\n        )\n\n    truth = {\n        protein: {aspect: set() for aspect in ASPECTS}\n        for protein in proteins\n    }\n\n    subset = frame[\n        frame[\"protein\"].isin(proteins)\n    ][[\"protein\", \"term\", \"aspect\"]]\n\n    for protein, term, aspect in subset.itertuples(index=False, name=None):\n        truth[protein][aspect].add(term)\n        truth[protein][aspect].update(ancestors(term))\n\n    return proteins, truth\n\n\ndef proxy_wfmax(proteins, truth, weights, blend_temperature):\n    prediction_cache = {\n        protein: predict_protein(protein, weights, blend_temperature)\n        for protein in proteins\n    }\n\n    best_by_aspect = {aspect: 0.0 for aspect in ASPECTS}\n    best_threshold = {aspect: None for aspect in ASPECTS}\n\n    for aspect in ASPECTS:\n        aspect_targets = [\n            protein\n            for protein in proteins\n            if sum(\n                float(term_ia.get(term, 0.0))\n                for term in truth[protein][aspect]\n            ) > 0\n        ]\n\n        for threshold in PROXY_THRESHOLDS:\n            precision_sum = 0.0\n            predicted_proteins = 0\n            recall_sum = 0.0\n\n            for protein in aspect_targets:\n                true_terms = truth[protein][aspect]\n                true_weight = sum(\n                    float(term_ia.get(term, 0.0))\n                    for term in true_terms\n                )\n\n                predicted_terms = [\n                    term\n                    for term, score, term_aspect in prediction_cache[protein]\n                    if term_aspect == aspect and score >= threshold\n                ]\n\n                true_positive_weight = sum(\n                    float(term_ia.get(term, 0.0))\n                    for term in predicted_terms\n                    if term in true_terms\n                )\n\n                predicted_weight = sum(\n                    float(term_ia.get(term, 0.0))\n                    for term in predicted_terms\n                )\n\n                if predicted_weight > 0:\n                    precision_sum += true_positive_weight / predicted_weight\n                    predicted_proteins += 1\n\n                recall_sum += true_positive_weight / true_weight\n\n            weighted_precision = (\n                precision_sum / predicted_proteins\n                if predicted_proteins\n                else 0.0\n            )\n            weighted_recall = (\n                recall_sum / len(aspect_targets)\n                if aspect_targets\n                else 0.0\n            )\n\n            if weighted_precision + weighted_recall > 0:\n                weighted_f = (\n                    2.0 * weighted_precision * weighted_recall\n                    / (weighted_precision + weighted_recall)\n                )\n            else:\n                weighted_f = 0.0\n\n            if weighted_f > best_by_aspect[aspect]:\n                best_by_aspect[aspect] = weighted_f\n                best_threshold[aspect] = float(threshold)\n\n    return (\n        float(np.mean(list(best_by_aspect.values()))),\n        best_by_aspect,\n        best_threshold,\n    )\n\n\nproxy_proteins, proxy_truth = build_proxy_truth(TRAIN_TERMS_PATH)\nshuffled = np.array(proxy_proteins, dtype=object)\nRNG.shuffle(shuffled)\n\nholdout_n = max(1, int(len(shuffled) * PROXY_HOLDOUT_FRACTION))\nproxy_holdout = sorted(shuffled[:holdout_n].tolist())\nproxy_calibration = sorted(shuffled[holdout_n:].tolist())\n\nprint(f\"Proxy proteins: {len(proxy_proteins):,}\")\nprint(f\"Calibration: {len(proxy_calibration):,}\")\nprint(f\"Holdout: {len(proxy_holdout):,}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:30:01.697552Z","iopub.execute_input":"2026-07-28T15:30:01.697983Z","iopub.status.idle":"2026-07-28T15:30:02.806250Z","shell.execute_reply.started":"2026-07-28T15:30:01.697961Z","shell.execute_reply":"2026-07-28T15:30:02.805518Z"}},"outputs":[],"execution_count":null},{"id":"f8b95771","cell_type":"code","source":"def random_parameter_candidate():\n    values = RNG.dirichlet(np.ones(3)).astype(np.float64)\n    return {\n        \"A\": float(values[0]),\n        \"B\": float(values[1]),\n        \"agreement\": float(values[2]),\n        \"temperature\": float(RNG.uniform(0.85, 1.35)),\n    }\n\n\nif RUN_TUNING:\n    selected = None\n    best_calibration_score = -np.inf\n\n    for _ in tqdm(range(TUNING_TRIALS), desc=\"Proxy tuning\"):\n        candidate = random_parameter_candidate()\n        score, _, _ = proxy_wfmax(\n            proxy_calibration,\n            proxy_truth,\n            weights={\n                \"A\": candidate[\"A\"],\n                \"B\": candidate[\"B\"],\n                \"agreement\": candidate[\"agreement\"],\n            },\n            blend_temperature=candidate[\"temperature\"],\n        )\n\n        if score > best_calibration_score:\n            best_calibration_score = score\n            selected = candidate\n\n    print(\"Best calibration proxy F-max:\", best_calibration_score)\nelse:\n    selected = RECORDED_PARAMETERS.copy()\n\nselected_weights = {\n    \"A\": selected[\"A\"],\n    \"B\": selected[\"B\"],\n    \"agreement\": selected[\"agreement\"],\n}\nselected_temperature = selected[\"temperature\"]\n\nprint(\"Selected parameters:\")\nprint(selected)\n\nholdout_score, holdout_aspects, holdout_thresholds = proxy_wfmax(\n    proxy_holdout,\n    proxy_truth,\n    weights=selected_weights,\n    blend_temperature=selected_temperature,\n)\n\nprint(\"Proxy holdout mean weighted F-max:\", holdout_score)\nprint(\"By ontology:\", holdout_aspects)\nprint(\"Best thresholds:\", holdout_thresholds)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:30:02.808259Z","iopub.execute_input":"2026-07-28T15:30:02.808465Z","iopub.status.idle":"2026-07-28T15:30:07.891146Z","shell.execute_reply.started":"2026-07-28T15:30:02.808445Z","shell.execute_reply":"2026-07-28T15:30:07.890420Z"}},"outputs":[],"execution_count":null},{"id":"cb31735d","cell_type":"code","source":"if RUN_PROXY_ABLATION:\n    variants = {\n        \"Model A only\": (\n            {\"A\": 1.0, \"B\": 0.0, \"agreement\": 0.0},\n            1.0,\n        ),\n        \"Model B only\": (\n            {\"A\": 0.0, \"B\": 1.0, \"agreement\": 0.0},\n            1.0,\n        ),\n        \"Equal A/B\": (\n            {\"A\": 0.5, \"B\": 0.5, \"agreement\": 0.0},\n            1.0,\n        ),\n        \"Final stack\": (\n            selected_weights,\n            selected_temperature,\n        ),\n    }\n\n    rows = []\n    for name, (weights, temperature) in variants.items():\n        mean_score, by_aspect, _ = proxy_wfmax(\n            proxy_holdout,\n            proxy_truth,\n            weights=weights,\n            blend_temperature=temperature,\n        )\n        rows.append(\n            {\n                \"variant\": name,\n                \"proxy_mean_wfmax\": mean_score,\n                \"BPO\": by_aspect[\"BPO\"],\n                \"CCO\": by_aspect[\"CCO\"],\n                \"MFO\": by_aspect[\"MFO\"],\n            }\n        )\n\n    proxy_ablation = pd.DataFrame(rows).sort_values(\n        \"proxy_mean_wfmax\",\n        ascending=False,\n    )\n    display(proxy_ablation.round(4))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:30:07.892289Z","iopub.execute_input":"2026-07-28T15:30:07.892619Z","iopub.status.idle":"2026-07-28T15:30:28.081897Z","shell.execute_reply.started":"2026-07-28T15:30:07.892586Z","shell.execute_reply":"2026-07-28T15:30:28.081190Z"}},"outputs":[],"execution_count":null},{"id":"e4ec6f89","cell_type":"markdown","source":"## Generate and validate the submission\n\nThe final file is tab-separated, has no header, uses scores in `(0, 1]` with up to three significant digits, and contains no more than 1,500 GO terms per target.\n","metadata":{}},{"id":"9c38626c","cell_type":"code","source":"line_count = 0\nproteins_with_predictions = 0\nmaximum_terms = 0\n\nwith open(OUTPUT_PATH, \"w\") as handle:\n    for protein in tqdm(target_ids, desc=\"Write submission\"):\n        predictions = predict_protein(\n            protein,\n            weights=selected_weights,\n            blend_temperature=selected_temperature,\n        )\n\n        if not predictions:\n            continue\n\n        predictions = predictions[:HARD_CAP]\n        proteins_with_predictions += 1\n        maximum_terms = max(maximum_terms, len(predictions))\n\n        for term, score, _ in predictions:\n            score_text = format(\n                float(np.clip(score, MIN_OUTPUT_SCORE, 1.0)),\n                f\".{SIGNIFICANT_DIGITS}g\",\n            )\n            handle.write(f\"{protein}\\t{term}\\t{score_text}\\n\")\n            line_count += 1\n\nprint(f\"Output: {OUTPUT_PATH}\")\nprint(f\"Lines: {line_count:,}\")\nprint(f\"Proteins with predictions: {proteins_with_predictions:,}\")\nprint(f\"Maximum terms for one protein: {maximum_terms:,}\")\nprint(f\"File size: {os.path.getsize(OUTPUT_PATH) / 1e6:.1f} MB\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T15:30:28.082778Z","iopub.execute_input":"2026-07-28T15:30:28.083157Z","iopub.status.idle":"2026-07-28T16:01:23.877540Z","shell.execute_reply.started":"2026-07-28T15:30:28.083117Z","shell.execute_reply":"2026-07-28T16:01:23.876531Z"}},"outputs":[],"execution_count":null},{"id":"c6b25e47","cell_type":"code","source":"def validate_submission(path):\n    total_lines = 0\n    current_protein = None\n    current_terms = set()\n    maximum_terms = 0\n\n    with open(path, \"r\") as handle:\n        for line_number, raw_line in enumerate(handle, start=1):\n            fields = raw_line.rstrip(\"\\n\").split(\"\\t\")\n            if len(fields) != 3:\n                raise ValueError(\n                    f\"Line {line_number}: expected 3 columns, found {len(fields)}\"\n                )\n\n            protein, term, score_text = fields\n\n            if protein not in target_set:\n                raise ValueError(f\"Line {line_number}: unknown target {protein}\")\n            if canonicalize(term) not in valid_terms:\n                raise ValueError(f\"Line {line_number}: invalid GO term {term}\")\n\n            score = float(score_text)\n            if not np.isfinite(score) or not (0.0 < score <= 1.0):\n                raise ValueError(f\"Line {line_number}: invalid score {score_text}\")\n\n            if protein != current_protein:\n                if current_protein is not None:\n                    maximum_terms = max(maximum_terms, len(current_terms))\n                    if len(current_terms) > HARD_CAP:\n                        raise ValueError(\n                            f\"{current_protein}: more than {HARD_CAP} terms\"\n                        )\n\n                current_protein = protein\n                current_terms = set()\n\n            if term in current_terms:\n                raise ValueError(\n                    f\"Line {line_number}: duplicate pair ({protein}, {term})\"\n                )\n\n            current_terms.add(term)\n            total_lines += 1\n\n    if current_protein is not None:\n        maximum_terms = max(maximum_terms, len(current_terms))\n        if len(current_terms) > HARD_CAP:\n            raise ValueError(f\"{current_protein}: more than {HARD_CAP} terms\")\n\n    print(\"Validation passed\")\n    print(f\"Total lines: {total_lines:,}\")\n    print(f\"Maximum terms per protein: {maximum_terms:,}\")\n\n\nvalidate_submission(OUTPUT_PATH)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-28T16:01:23.878570Z","iopub.execute_input":"2026-07-28T16:01:23.878886Z","iopub.status.idle":"2026-07-28T16:02:04.443286Z","shell.execute_reply.started":"2026-07-28T16:01:23.878850Z","shell.execute_reply":"2026-07-28T16:02:04.442483Z"}},"outputs":[],"execution_count":null},{"id":"fdd14302","cell_type":"markdown","source":"## Details of the submission\n\n### What mattered most\n\n**1. Diagnosing score spaces before blending**  \nThe largest practical mistake would have been to average unbounded values with probabilities. The source diagnostics determine whether a temperature-scaled sigmoid is required.\n\n**2. Rewarding agreement rather than taking a large union**  \nThe agreement signal favors GO terms supported by both systems. Strong unique predictions can remain, but weak source-specific terms are suppressed.\n\n**3. Using the ontology and evaluation weights during post-processing**  \nAlternative-ID handling, hard-max ancestor propagation, IA-aware filtering, and the output cap are part of the predictive strategy rather than formatting details.\n\n**4. Keeping the pipeline reproducible and memory-aware**  \nChunked loading makes the workflow practical on Kaggle. The public version fixes the random-number generation and records the final parameters.\n\n### What did not work well\n\n- A direct arithmetic average was not appropriate because the source score ranges differed.\n- A pure union of both files retained too many weak predictions.\n- The competition-time random search was not fully reproducible because NumPy had not been seeded.\n- The proxy based on current annotations is optimistic for GOA-derived predictions and should not be treated as a leaderboard estimator.\n- One global parameter set was used for BP, CC, and MF; ontology-specific calibration is a natural next step.\n\n### Limitations\n\nThis notebook reconstructs the ensemble and post-processing stage. It depends on two public upstream outputs and does not train a new protein language model. A stronger future version would use time-aware validation, ontology-specific calibration, taxonomy-aware reliability, and direct comparison with the official CAFA evaluator.\n","metadata":{}},{"id":"aa3b9d5a","cell_type":"markdown","source":"## Sources and acknowledgments\n\n- [CAFA 6 Protein Function Prediction](https://www.kaggle.com/competitions/cafa-6-protein-function-prediction)\n- [CAFA ensemble 0.378](https://www.kaggle.com/code/antonoof/cafa-ensemble-0-378) — antonoof\n- [CAFA6 GOA + ProtT5 IC-Aware](https://www.kaggle.com/code/maryamkhaled67/cafa6-goa-prott5-ic-aware) — maryamkhaled67\n- [My original competition notebook](https://www.kaggle.com/code/ismailcanoguz/mybestscore-372?scriptVersionId=295439875)\n- [CAFA-evaluator](https://github.com/BioComputingUP/CAFA-evaluator)\n- Gene Ontology and Information Accretion files supplied by the competition organizers\n\nThank you to the public notebook authors and the CAFA community for making their work available.\n","metadata":{}}]}