{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.12.12"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":118765,"databundleVersionId":15231210},{"sourceType":"datasetVersion","sourceId":11025231,"datasetId":6865700,"databundleVersionId":11407570},{"sourceType":"datasetVersion","sourceId":14834274,"datasetId":9487358,"databundleVersionId":15692465}],"dockerImageVersionId":31259,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":735.869346,"end_time":"2026-02-23T12:45:46.73553","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2026-02-23T12:33:30.866184","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# RNA 3D Folding Part 2 — Persistent Homology + XGBoost (Train: no PKL)\n\nThis notebook **builds features on-the-fly without saving intermediate pkl files**, and trains directly in memory. To prevent PH computation from exploding, **landmark subsampling** and **Rips max_edge percentile clipping** are applied.\n\n- Objective: Regress per-residue 3D coordinates $(x,y,z)$ from RNA sequences\n- Training: **NO PKL** (builds `train_features` in memory)\n- Inference: Uses the same feature pipeline to generate `submission.csv`\n\n*Designed for Kaggle offline environments (`gudhi` installed from a local wheel).*","metadata":{"papermill":{"duration":0.006099,"end_time":"2026-02-23T12:33:34.227694","exception":false,"start_time":"2026-02-23T12:33:34.221595","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"# RNA 3D Folding — Persistent Homology + XGBoost\n\n## Overview\n\nGiven an RNA sequence $s = (s_1, s_2, \\ldots, s_n),\\ s_i \\in \\{A, U, G, C\\}$, the task is to regress the 3D coordinates $\\hat{y}_i = (\\hat{x}_i, \\hat{y}_i, \\hat{z}_i)$ of each residue $i$.\n\n**Pipeline:**\n\n$$\n\\text{RNA sequence} \\xrightarrow{\\text{Nussinov}} \\text{base pairs} \\xrightarrow{\\text{Dijkstra}} \\text{graph distance matrix} \\xrightarrow{\\text{PH}} \\text{topological features} \\xrightarrow{\\text{XGBoost}} \\hat{x}, \\hat{y}, \\hat{z}\n$$\n\n---\n\n## Step 1: Secondary Structure Prediction — Nussinov DP\n\nThe set of base pairs $P \\subseteq \\{(i,j) \\mid i < j\\}$ is computed via dynamic programming.\n\n$$\ndp[i][j] = \\max \\begin{cases}\ndp[i][j-1] & \\text{($j$ is unpaired)} \\\\\n\\max_{k: (s_k, s_j) \\in \\mathcal{P},\\, k < j - L} \\left( dp[i][k-1] + 1 + dp[k+1][j-1] \\right)\n\\end{cases}\n$$\n\n- $\\mathcal{P} = \\{(A,U),\\,(U,A),\\,(G,C),\\,(C,G),\\,(G,U),\\,(U,G)\\}$: Watson–Crick + Wobble base pairs\n- $L = 3$: minimum loop length (stereochemical constraint)\n- Time complexity: $O(n^3)$\n\n---\n\n## Step 2: Weighted Graph Distance Matrix — Dijkstra\n\nA weighted graph $G = (V, E, w)$ is constructed with residues as nodes, and edges between adjacent residues and base-paired residues.\n\n$$\nw(i, j) = \\begin{cases}\n1.0 & \\text{if } |i - j| = 1 \\quad \\text{(backbone bond)} \\\\\n\\alpha & \\text{if } (i,j) \\in P \\quad \\text{(base pair, } \\alpha = 0.7 \\text{)}\n\\end{cases}\n$$\n\nThe all-pairs shortest path matrix $D \\in \\mathbb{R}^{n \\times n}$ is computed with Dijkstra's algorithm:\n\n$$\nD[s][t] = \\min_{\\text{path } s \\to t} \\sum_{(u,v) \\in \\text{path}} w(u,v)\n$$\n\nSetting $\\alpha < 1$ makes base-paired residues appear closer to each other.\n\n---\n\n## Step 3: Landmark Subsampling\n\nFor long sequences, computing the $n \\times n$ distance matrix requires $O(n^2)$ memory. We select $N_s = 120$ landmark nodes to subsample.\n\n$$\nI = I_{\\text{paired}} \\cup I_{\\text{uniform}}, \\quad |I| \\leq N_s\n$$\n\n- $I_{\\text{paired}}$: nodes involved in base pairs (structurally important)\n- $I_{\\text{uniform}}$: uniformly spaced nodes to fill the remainder\n- Subsampled distance matrix: $D_s = D[I, I] \\in \\mathbb{R}^{N_s \\times N_s}$\n\n---\n\n## Step 4: Topological Data Analysis — Persistent Homology\n\nA Vietoris–Rips complex is built from $D_s$, and the persistence diagram of 1-dimensional homology $H_1$ (loops/holes) is computed.\n\n$$\n\\text{PD}_1 = \\{(b_k, d_k)\\}_{k=1}^{K}\n$$\n\n- $b_k$: scale at which the $k$-th loop is born (birth)\n- $d_k$: scale at which the $k$-th loop dies (death)\n- $\\text{life}_k = d_k - b_k$: persistence (robustness of the loop)\n\nThe following features are extracted from each loop (top $K=12$):\n\n$$\n\\text{features}_k = \\left( b_k,\\; d_k,\\; \\underbrace{d_k - b_k}_{\\text{life}},\\; \\underbrace{d_k / (b_k + \\varepsilon)}_{\\text{ratio}} \\right)\n$$\n\nAdditional global statistics:\n\n$$\n\\text{life}_{\\max} = \\max_k \\text{life}_k, \\quad\n\\text{life}_{\\text{sum}} = \\sum_k \\text{life}_k, \\quad\n\\bar{b} = \\frac{1}{K}\\sum_k b_k, \\quad\n\\bar{d} = \\frac{1}{K}\\sum_k d_k\n$$\n\nThe `max_edge_length` is clipped at the $p = 75$th percentile of the upper-triangular entries of the distance matrix to prevent computational explosion.\n\n---\n\n## Step 5: Local Features\n\nLocal features are appended for each residue $i$.\n\n$$\n\\mathbf{f}_i^{\\text{local}} = \\left[\n\\underbrace{\\mathbf{1}[s_i = A],\\, \\mathbf{1}[s_i = C],\\, \\mathbf{1}[s_i = G],\\, \\mathbf{1}[s_i = U]}_{\\text{one-hot}},\\;\n\\underbrace{\\frac{i-1}{n-1}}_{\\text{position}},\\;\n\\underbrace{\\deg_P(i)}_{\\text{base-pair degree}},\\;\n\\underbrace{\\mathbf{1}[\\deg_P(i) > 0]}_{\\text{paired flag}}\n\\right]\n$$\n\nFinal feature vector:\n\n$$\n\\mathbf{f}_i = \\left[ \\mathbf{f}_i^{\\text{local}} \\;\\|\\; \\mathbf{f}^{\\text{PH}} \\right] \\in \\mathbb{R}^{59}\n$$\n\nThe global topological features $\\mathbf{f}^{\\text{PH}}$ are shared across all residues within the same RNA.\n\n---\n\n## Step 6: XGBoost Regression\n\nAn independent regression model is trained for each coordinate axis:\n\n$$\n\\hat{x}_i = f_x(\\mathbf{f}_i), \\quad \\hat{y}_i = f_y(\\mathbf{f}_i), \\quad \\hat{z}_i = f_z(\\mathbf{f}_i)\n$$\n\nLoss function (RMSE):\n\n$$\n\\mathcal{L} = \\sqrt{\\frac{1}{N} \\sum_{i=1}^{N} (\\hat{c}_i - c_i)^2}, \\quad c \\in \\{x, y, z\\}\n$$\n\n","metadata":{"papermill":{"duration":0.004553,"end_time":"2026-02-23T12:33:34.237044","exception":false,"start_time":"2026-02-23T12:33:34.232491","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# =============\n# 1) Install\n# =============\n# Assumes Kaggle offline wheels (adjust the path to match your wheel location)\n!pip install --no-index --find-links=/kaggle/input/datasets/takahiro3110/gungun-cp312/gudhi_wheels_cp312/wheels gudhi","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:33:34.248043Z","iopub.status.busy":"2026-02-23T12:33:34.247647Z","iopub.status.idle":"2026-02-23T12:33:40.027917Z","shell.execute_reply":"2026-02-23T12:33:40.026549Z"},"papermill":{"duration":5.789193,"end_time":"2026-02-23T12:33:40.03057","exception":false,"start_time":"2026-02-23T12:33:34.241377","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 2) Imports\n# =============\nimport os, gc, heapq, warnings\nimport numpy as np\nimport pandas as pd\nimport gudhi\nimport xgboost as xgb\nfrom tqdm import tqdm\n\nwarnings.filterwarnings(\"ignore\")\nnp.set_printoptions(suppress=True)","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:33:40.043873Z","iopub.status.busy":"2026-02-23T12:33:40.04284Z","iopub.status.idle":"2026-02-23T12:33:43.290161Z","shell.execute_reply":"2026-02-23T12:33:43.289009Z"},"papermill":{"duration":3.256533,"end_time":"2026-02-23T12:33:43.29246","exception":false,"start_time":"2026-02-23T12:33:40.035927","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 3) Config\n# =============\nSEED = 42\nrng = np.random.default_rng(SEED)\n\n# ---- Start small for quick testing (increase gradually) ----\nSAMPLE_N_TARGETS = 300  # None = use all. Recommend 100-300 to start.\n\n# ---- PH stabilization ----\nSUBSAMPLE_N  = 120        # Max nodes for TDA (long sequences use landmarks)\nEDGE_PCTILE  = 75         # Clip Rips max_edge_length at this percentile\nTOP_K        = 12         # Number of top H1 barcodes to use\nBP_ALPHA     = 0.7        # Weight of base-pair edges (smaller = closer)\nMIN_LOOP     = 3          # Nussinov minimum loop length\n\n# ---- Training settings ----\nNUM_ROUND_MAX = 3000\nEARLY_STOP    = 100\n\nprint(\"Config loaded.\")","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:33:43.304217Z","iopub.status.busy":"2026-02-23T12:33:43.303746Z","iopub.status.idle":"2026-02-23T12:33:43.310646Z","shell.execute_reply":"2026-02-23T12:33:43.309462Z"},"papermill":{"duration":0.015226,"end_time":"2026-02-23T12:33:43.312749","exception":false,"start_time":"2026-02-23T12:33:43.297523","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 4) Load data (Kaggle Stanford RNA 3D Folding 2)\n# =============\ntrain_seq_file   = '/kaggle/input/stanford-rna-3d-folding-2/train_sequences.csv'\ntrain_labels_file= '/kaggle/input/stanford-rna-3d-folding-2/train_labels.csv'\nval_seq_file     = '/kaggle/input/stanford-rna-3d-folding-2/validation_sequences.csv'\nval_labels_file  = '/kaggle/input/stanford-rna-3d-folding-2/validation_labels.csv'\ntest_seq_file    = '/kaggle/input/stanford-rna-3d-folding-2/test_sequences.csv'\nsample_sub_file  = '/kaggle/input/stanford-rna-3d-folding-2/sample_submission.csv'\n\ntrain_sequences = pd.read_csv(train_seq_file).dropna(subset=[\"sequence\"]).reset_index(drop=True)\ntrain_labels_raw= pd.read_csv(train_labels_file)\n\nval_sequences   = pd.read_csv(val_seq_file).dropna(subset=[\"sequence\"]).reset_index(drop=True)\nval_labels      = pd.read_csv(val_labels_file)\n\ntest_sequences  = pd.read_csv(test_seq_file).dropna(subset=[\"sequence\"]).reset_index(drop=True)\nsample_sub      = pd.read_csv(sample_sub_file)\n\nprint(\"train_sequences:\", train_sequences.shape)\nprint(\"train_labels_raw:\", train_labels_raw.shape)\nprint(\"val_sequences:\", val_sequences.shape)\nprint(\"val_labels:\", val_labels.shape)\nprint(\"test_sequences:\", test_sequences.shape)\nprint(\"sample_sub:\", sample_sub.shape)","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:33:43.324251Z","iopub.status.busy":"2026-02-23T12:33:43.323919Z","iopub.status.idle":"2026-02-23T12:33:54.826635Z","shell.execute_reply":"2026-02-23T12:33:54.825484Z"},"papermill":{"duration":11.5112,"end_time":"2026-02-23T12:33:54.829016","exception":false,"start_time":"2026-02-23T12:33:43.317816","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 5) Normalize train_labels_raw\n#    ID -> (target_id, resid), coord columns -> x_1,y_1,z_1\n# =============\ncand_cols = [c for c in [\"ID\",\"x_1\",\"y_1\",\"z_1\",\"x\",\"y\",\"z\"] if c in train_labels_raw.columns]\ntrain_labels_raw = train_labels_raw.dropna(subset=cand_cols).reset_index(drop=True)\n\nspl = train_labels_raw[\"ID\"].astype(str).str.split(\"_\", n=1, expand=True)\ntrain_labels_raw[\"target_id\"] = spl[0]\nif \"resid\" not in train_labels_raw.columns:\n    train_labels_raw[\"resid\"] = pd.to_numeric(spl[1], errors=\"coerce\")\n\nrename_map = {}\nif \"x_1\" not in train_labels_raw.columns and \"x\" in train_labels_raw.columns: rename_map[\"x\"]=\"x_1\"\nif \"y_1\" not in train_labels_raw.columns and \"y\" in train_labels_raw.columns: rename_map[\"y\"]=\"y_1\"\nif \"z_1\" not in train_labels_raw.columns and \"z\" in train_labels_raw.columns: rename_map[\"z\"]=\"z_1\"\nif rename_map:\n    train_labels_raw = train_labels_raw.rename(columns=rename_map)\n\ntrain_labels_raw[\"resid\"] = pd.to_numeric(train_labels_raw[\"resid\"], errors=\"coerce\")\ntrain_labels_raw = train_labels_raw.dropna(subset=[\"resid\"]).copy()\ntrain_labels_raw[\"resid\"] = train_labels_raw[\"resid\"].astype(np.int32)\n\nneed = {\"target_id\",\"resid\",\"x_1\",\"y_1\",\"z_1\"}\nmissing = need - set(train_labels_raw.columns)\nassert not missing, f\"train_labels_raw missing columns: {missing}\"\n\nprint(\"Normalized train_labels_raw:\", train_labels_raw.shape)\ntrain_labels_raw.head()","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:33:54.841794Z","iopub.status.busy":"2026-02-23T12:33:54.840839Z","iopub.status.idle":"2026-02-23T12:34:10.208651Z","shell.execute_reply":"2026-02-23T12:34:10.207675Z"},"papermill":{"duration":15.376273,"end_time":"2026-02-23T12:34:10.210829","exception":false,"start_time":"2026-02-23T12:33:54.834556","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 6) Optional: sample targets (for quick debug)\n# =============\nif SAMPLE_N_TARGETS is not None:\n    keep_ids = train_sequences[\"target_id\"].unique()[:SAMPLE_N_TARGETS]\n    train_sequences = train_sequences[train_sequences[\"target_id\"].isin(keep_ids)].reset_index(drop=True)\n    train_labels_raw= train_labels_raw[train_labels_raw[\"target_id\"].isin(keep_ids)].reset_index(drop=True)\n\nprint(\"train_sequences:\", train_sequences.shape, \"targets:\", train_sequences[\"target_id\"].nunique())\nprint(\"train_labels_raw:\", train_labels_raw.shape, \"targets:\", train_labels_raw[\"target_id\"].nunique())","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:34:10.224315Z","iopub.status.busy":"2026-02-23T12:34:10.223936Z","iopub.status.idle":"2026-02-23T12:34:10.707958Z","shell.execute_reply":"2026-02-23T12:34:10.706366Z"},"papermill":{"duration":0.493465,"end_time":"2026-02-23T12:34:10.710157","exception":false,"start_time":"2026-02-23T12:34:10.216692","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 7) Base-pair prediction: Nussinov DP\n# =============\n_PAIRS = frozenset({\n    ('A','U'), ('U','A'),\n    ('G','C'), ('C','G'),\n    ('G','U'), ('U','G'),\n})\n\ndef _nussinov_core(seq: str, min_loop: int):\n    \"\"\"Original Nussinov DP (for short sequences only)\"\"\"\n    n = len(seq)\n    if n < min_loop + 2:\n        return []\n    dp = np.zeros((n,n), dtype=np.int16)\n\n    for span in range(min_loop+1, n):\n        for i in range(n-span):\n            j = i+span\n            best = int(dp[i, j-1])\n            bj = seq[j]\n            for k in range(i, j-min_loop):\n                if (seq[k], bj) in _PAIRS:\n                    left  = int(dp[i, k-1]) if k>i else 0\n                    inner = int(dp[k+1, j-1]) if k+1<=j-1 else 0\n                    val = left + 1 + inner\n                    if val > best:\n                        best = val\n            dp[i,j] = best\n\n    pairs = []\n    stack = [(0, n-1)]\n    while stack:\n        i,j = stack.pop()\n        if i>=j:\n            continue\n        if dp[i,j] == dp[i, j-1]:\n            stack.append((i, j-1))\n            continue\n        bj = seq[j]\n        for k in range(i, j-min_loop):\n            if (seq[k], bj) in _PAIRS:\n                left  = int(dp[i, k-1]) if k>i else 0\n                inner = int(dp[k+1, j-1]) if k+1<=j-1 else 0\n                if left + 1 + inner == dp[i,j]:\n                    pairs.append((k,j))\n                    if k>i:\n                        stack.append((i, k-1))\n                    stack.append((k+1, j-1))\n                    break\n    return pairs\n\n\ndef predict_base_pairs(seq: str, min_loop: int = 3, max_n_nuss: int = 600):\n    seq = seq.upper().replace(\"T\",\"U\")\n    n = len(seq)\n\n    if n <= max_n_nuss:\n        # Short sequences: run Nussinov directly\n        return _nussinov_core(seq, min_loop)\n\n    # Long sequences: uniform subsample -> Nussinov -> map back to original indices\n    step = max(1, n // max_n_nuss)\n    indices = list(range(0, n, step))[:max_n_nuss]\n    sub_seq = \"\".join(seq[i] for i in indices)\n    raw_pairs = _nussinov_core(sub_seq, min_loop)\n    return [(indices[i], indices[j]) for i, j in raw_pairs]\n\n# quick test\nprint(predict_base_pairs(\"GGGAAAUCC\", min_loop=MIN_LOOP)[:10])","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:34:10.722864Z","iopub.status.busy":"2026-02-23T12:34:10.722499Z","iopub.status.idle":"2026-02-23T12:34:10.73817Z","shell.execute_reply":"2026-02-23T12:34:10.737172Z"},"papermill":{"duration":0.024943,"end_time":"2026-02-23T12:34:10.74035","exception":false,"start_time":"2026-02-23T12:34:10.715407","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 8) Weighted graph distance matrix (Dijkstra)\n# =============\ndef graph_distance_matrix_weighted(n: int, pairs, alpha: float=0.7) -> np.ndarray:\n    adj = [[] for _ in range(n)]\n    for i in range(n-1):\n        adj[i].append((i+1, 1.0))\n        adj[i+1].append((i, 1.0))\n    w_bp = float(alpha)\n    for i,j in pairs:\n        if 0 <= i < n and 0 <= j < n and i != j:\n            adj[i].append((j, w_bp))\n            adj[j].append((i, w_bp))\n\n    D = np.full((n,n), np.inf, dtype=np.float32)\n    for s in range(n):\n        D[s,s] = 0.0\n        pq = [(0.0, s)]\n        while pq:\n            d,v = heapq.heappop(pq)\n            if d != D[s,v]:\n                continue\n            for u,w in adj[v]:\n                nd = d + w\n                if nd < D[s,u]:\n                    D[s,u] = nd\n                    heapq.heappush(pq, (nd,u))\n\n    fin = np.isfinite(D)\n    if not fin.all():\n        max_f = float(D[fin].max()) if fin.any() else float(n)\n        D[~fin] = max_f + w_bp\n    return D\n\n# quick test\nseq = \"GGGAAAUCC\"\npairs = predict_base_pairs(seq, MIN_LOOP)\nD = graph_distance_matrix_weighted(len(seq), pairs, BP_ALPHA)\nprint(\"D shape:\", D.shape, \"finite:\", np.isfinite(D).all())","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:34:10.752604Z","iopub.status.busy":"2026-02-23T12:34:10.752248Z","iopub.status.idle":"2026-02-23T12:34:10.764337Z","shell.execute_reply":"2026-02-23T12:34:10.763032Z"},"papermill":{"duration":0.020695,"end_time":"2026-02-23T12:34:10.766438","exception":false,"start_time":"2026-02-23T12:34:10.745743","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 9) Landmark subsampling (prevents PH explosion)\n# =============\ndef subsample_D(D_full: np.ndarray, pairs, max_n: int) -> np.ndarray:\n    n = D_full.shape[0]\n    if n <= max_n:\n        return D_full\n    paired_nodes = sorted({i for i,j in pairs} | {j for i,j in pairs})\n    keep = list(dict.fromkeys(paired_nodes))\n    remaining = max_n - len(keep)\n    if remaining > 0:\n        step = max(1, n // remaining)\n        extra = [i for i in range(0, n, step) if i not in set(keep)]\n        keep += extra\n    keep = sorted(set(keep))[:max_n]\n    idx = np.array(keep, dtype=np.int32)\n    return D_full[np.ix_(idx, idx)]\n\nD2 = subsample_D(D, pairs, max_n=5)\nprint(\"subsample:\", D2.shape)","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:34:10.779005Z","iopub.status.busy":"2026-02-23T12:34:10.778636Z","iopub.status.idle":"2026-02-23T12:34:10.787444Z","shell.execute_reply":"2026-02-23T12:34:10.786446Z"},"papermill":{"duration":0.017798,"end_time":"2026-02-23T12:34:10.789751","exception":false,"start_time":"2026-02-23T12:34:10.771953","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 10) PH(H1) features via Vietoris-Rips\n# =============\ndef h1_features(D_full: np.ndarray, pairs, top_k:int=12, edge_pctile:float=75) -> dict:\n    D = subsample_D(D_full, pairs, max_n=SUBSAMPLE_N)\n    upper = D[np.triu_indices_from(D, k=1)]\n    upper = upper[np.isfinite(upper)]\n    if len(upper) == 0:\n        max_edge = 1.0\n    else:\n        max_edge = float(np.percentile(upper, edge_pctile))\n    del upper\n\n    rips = gudhi.RipsComplex(distance_matrix=D, max_edge_length=max_edge)\n    st = rips.create_simplex_tree(max_dimension=2)\n    st.persistence()\n\n    h1 = st.persistence_intervals_in_dimension(1)\n    if len(h1) == 0:\n        feats = {}\n        for k in range(top_k):\n            feats[f\"h1_birth_{k}\"] = 0.0\n            feats[f\"h1_death_{k}\"] = 0.0\n            feats[f\"h1_life_{k}\"]  = 0.0\n            feats[f\"h1_ratio_{k}\"] = 0.0\n        feats[\"h1_max_life\"] = 0.0\n        feats[\"h1_sum_persistence\"] = 0.0\n        feats[\"h1_mean_birth\"] = 0.0\n        feats[\"h1_mean_death\"] = 0.0\n        return feats\n\n    births = h1[:,0].astype(np.float32)\n    deaths = h1[:,1].astype(np.float32)\n    life   = (deaths - births).astype(np.float32)\n\n    idx = np.argsort(life)[::-1]\n    births, deaths, life = births[idx], deaths[idx], life[idx]\n\n    eps = 1e-6\n    feats = {}\n    K = min(top_k, len(life))\n    for k in range(top_k):\n        if k < K:\n            b = float(births[k]); d = float(deaths[k]); l = float(life[k])\n            feats[f\"h1_birth_{k}\"] = b\n            feats[f\"h1_death_{k}\"] = d\n            feats[f\"h1_life_{k}\"]  = l\n            feats[f\"h1_ratio_{k}\"] = float(d/(b+eps))\n        else:\n            feats[f\"h1_birth_{k}\"] = 0.0\n            feats[f\"h1_death_{k}\"] = 0.0\n            feats[f\"h1_life_{k}\"]  = 0.0\n            feats[f\"h1_ratio_{k}\"] = 0.0\n\n    feats[\"h1_max_life\"] = float(life.max())\n    feats[\"h1_sum_persistence\"] = float(life.sum())\n    feats[\"h1_mean_birth\"] = float(births.mean())\n    feats[\"h1_mean_death\"] = float(deaths.mean())\n    return feats\n\nph = h1_features(D, pairs, top_k=TOP_K, edge_pctile=EDGE_PCTILE)\n{k: ph[k] for k in list(ph.keys())[:8]}","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:34:10.803389Z","iopub.status.busy":"2026-02-23T12:34:10.802119Z","iopub.status.idle":"2026-02-23T12:34:10.823515Z","shell.execute_reply":"2026-02-23T12:34:10.822357Z"},"papermill":{"duration":0.030405,"end_time":"2026-02-23T12:34:10.825708","exception":false,"start_time":"2026-02-23T12:34:10.795303","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 10b) Persistence diagram visualization\n# =============\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\n\ndef plot_persistence_diagram(seq: str, target_id: str = \"\"):\n    seq = seq.upper().replace(\"T\", \"U\")\n    n = len(seq)\n    pairs = predict_base_pairs(seq, min_loop=MIN_LOOP)\n    D = graph_distance_matrix_weighted(n, pairs, alpha=BP_ALPHA)\n    D_sub = subsample_D(D, pairs, max_n=SUBSAMPLE_N)\n\n    upper = D_sub[np.triu_indices_from(D_sub, k=1)]\n    upper = upper[np.isfinite(upper)]\n    max_edge = float(np.percentile(upper, EDGE_PCTILE)) if len(upper) > 0 else 1.0\n\n    rips = gudhi.RipsComplex(distance_matrix=D_sub, max_edge_length=max_edge)\n    st = rips.create_simplex_tree(max_dimension=2)\n    st.persistence()\n\n    h0 = st.persistence_intervals_in_dimension(0)\n    h1 = st.persistence_intervals_in_dimension(1)\n\n    fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n    title = f\"Persistence Diagram — {target_id} (n={n}, pairs={len(pairs)})\"\n    fig.suptitle(title, fontsize=13)\n\n    for ax, intervals, label, color in zip(\n        axes,\n        [h0, h1],\n        [\"H₀ (connected components)\", \"H₁ (loops)\"],\n        [\"steelblue\", \"tomato\"],\n    ):\n        # diagonal line\n        all_finite = intervals[np.isfinite(intervals[:, 1])] if len(intervals) > 0 else np.array([]).reshape(0,2)\n        lim_max = float(all_finite.max()) * 1.1 if len(all_finite) > 0 else 2.0\n        ax.plot([0, lim_max], [0, lim_max], \"k--\", lw=0.8, alpha=0.5)\n\n        if len(intervals) > 0:\n            finite_mask = np.isfinite(intervals[:, 1])\n            # finite points\n            ax.scatter(intervals[finite_mask, 0], intervals[finite_mask, 1],\n                       c=color, s=40, alpha=0.7, zorder=3, label=\"finite\")\n            # infinite points (death=inf) shown at top\n            inf_births = intervals[~finite_mask, 0]\n            if len(inf_births) > 0:\n                ax.scatter(inf_births, [lim_max] * len(inf_births),\n                           c=color, s=60, marker=\"^\", alpha=0.9, label=\"∞\")\n\n        ax.set_xlabel(\"Birth\", fontsize=11)\n        ax.set_ylabel(\"Death\", fontsize=11)\n        ax.set_title(label, fontsize=11)\n        ax.set_xlim(-0.05, lim_max)\n        ax.set_ylim(-0.05, lim_max * 1.05)\n        ax.legend(fontsize=9)\n        ax.grid(True, alpha=0.3)\n\n    plt.tight_layout()\n    plt.show()\n    print(f\"H0: {len(h0)} bars  |  H1: {len(h1)} bars  |  top-H1 life: \"\n          + \", \".join(f\"{v:.3f}\" for v in sorted((h1[:,1]-h1[:,0])[np.isfinite(h1[:,1])], reverse=True)[:5]))\n\n# --- Plot: specify target_id here ---\nTARGET_ID = \"1ASZ\"  # <- change this to any target_id\n\nrow = train_sequences[train_sequences[\"target_id\"] == TARGET_ID]\nif len(row) == 0:\n    print(f\"{TARGET_ID} not found. Available target_ids:\")\n    print(train_sequences[\"target_id\"].tolist())\nelse:\n    plot_persistence_diagram(row[\"sequence\"].iloc[0], target_id=TARGET_ID)\n\n# Uncomment to compare multiple RNAs\nfor i in range(30):\n     plot_persistence_diagram(\n         train_sequences[\"sequence\"].iloc[i],\n         target_id=train_sequences[\"target_id\"].iloc[i]\n     )","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:34:10.838858Z","iopub.status.busy":"2026-02-23T12:34:10.838501Z","iopub.status.idle":"2026-02-23T12:34:24.300151Z","shell.execute_reply":"2026-02-23T12:34:24.299181Z"},"papermill":{"duration":13.471955,"end_time":"2026-02-23T12:34:24.303281","exception":false,"start_time":"2026-02-23T12:34:10.831326","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 11) Local features (one-hot + position + pairing flags)\n# =============\ndef one_hot_base(b: str):\n    return {\n        \"oh_A\": 1.0 if b==\"A\" else 0.0,\n        \"oh_C\": 1.0 if b==\"C\" else 0.0,\n        \"oh_G\": 1.0 if b==\"G\" else 0.0,\n        \"oh_U\": 1.0 if b==\"U\" else 0.0,\n    }\n\ndef build_rows_for_target(target_id: str, seq: str):\n    seq = seq.upper().replace(\"T\",\"U\")\n    n = len(seq)\n    pairs = predict_base_pairs(seq, min_loop=MIN_LOOP)\n\n    D = graph_distance_matrix_weighted(n, pairs, alpha=BP_ALPHA)\n    ph = h1_features(D, pairs, top_k=TOP_K, edge_pctile=EDGE_PCTILE)\n\n    deg = np.zeros(n, dtype=np.int16)\n    for i,j in pairs:\n        deg[i]+=1; deg[j]+=1\n\n    rows = []\n    for resid in range(1, n+1):\n        i = resid - 1\n        b = seq[i]\n        r = {\n            \"ID\": f\"{target_id}_{resid}\",\n            \"target_id\": target_id,\n            \"resid\": int(resid),\n            \"pos_norm\": float(i / max(n-1,1)),\n            \"deg_pair\": float(deg[i]),\n            \"is_paired\": 1.0 if deg[i] > 0 else 0.0,\n        }\n        r.update(one_hot_base(b))\n        r.update(ph)  # global PH features are shared across all residues in the same RNA\n        rows.append(r)\n\n    return rows\n\n# quick test\ntid = train_sequences[\"target_id\"].iloc[0]\nseq = train_sequences[\"sequence\"].iloc[0]\nrows = build_rows_for_target(tid, seq)\nprint(\"rows:\", len(rows), \"keys:\", len(rows[0]))","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:34:24.365469Z","iopub.status.busy":"2026-02-23T12:34:24.365101Z","iopub.status.idle":"2026-02-23T12:34:24.427783Z","shell.execute_reply":"2026-02-23T12:34:24.42667Z"},"papermill":{"duration":0.095756,"end_time":"2026-02-23T12:34:24.430158","exception":false,"start_time":"2026-02-23T12:34:24.334402","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 12) Build train features (NO PKL)\n# =============\ndef build_train_features_no_pkl(train_sequences: pd.DataFrame, train_labels_raw: pd.DataFrame):\n    feature_rows = []\n    yx, yy, yz = [], [], []\n\n    grp = {tid: df for tid, df in train_labels_raw.groupby(\"target_id\", sort=False)}\n    tids = train_sequences[\"target_id\"].unique()\n\n    for tid in tqdm(tids, desc=\"Build train features\"):\n        seq = train_sequences.loc[train_sequences[\"target_id\"]==tid, \"sequence\"].iloc[0]\n        lab = grp.get(tid)\n        if lab is None or len(lab)==0:\n            continue\n\n        # resid -> (x,y,z)\n        m = {int(r): (float(x),float(y),float(z))\n             for r,x,y,z in zip(lab[\"resid\"].values, lab[\"x_1\"].values, lab[\"y_1\"].values, lab[\"z_1\"].values)}\n\n        rows = build_rows_for_target(tid, seq)\n\n        for r in rows:\n            resid = int(r[\"resid\"])\n            if resid not in m:\n                continue\n            x,y,z = m[resid]\n            feature_rows.append(r)\n            yx.append(x); yy.append(y); yz.append(z)\n\n        del rows, m\n        gc.collect()\n\n    train_features = pd.DataFrame(feature_rows)\n    yx = np.asarray(yx, dtype=np.float32)\n    yy = np.asarray(yy, dtype=np.float32)\n    yz = np.asarray(yz, dtype=np.float32)\n    return train_features, yx, yy, yz\n\ntrain_features, y_x, y_y, y_z = build_train_features_no_pkl(train_sequences, train_labels_raw)\nprint(\"train_features:\", train_features.shape, \"y:\", y_x.shape)\ntrain_features.head()","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:34:24.49443Z","iopub.status.busy":"2026-02-23T12:34:24.494048Z","iopub.status.idle":"2026-02-23T12:43:57.767971Z","shell.execute_reply":"2026-02-23T12:43:57.766956Z"},"papermill":{"duration":573.305617,"end_time":"2026-02-23T12:43:57.770096","exception":false,"start_time":"2026-02-23T12:34:24.464479","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 13) Make X / feature_columns & impute NaN/inf\n# =============\ndef impute_X(X: np.ndarray) -> np.ndarray:\n    X = np.where(np.isinf(X), np.nan, X.copy())\n    col_medians = np.nanmedian(X, axis=0)\n    col_medians = np.where(np.isnan(col_medians), 0.0, col_medians)\n    nan_mask = np.isnan(X)\n    X[nan_mask] = np.take(col_medians, np.where(nan_mask)[1])\n    return X.astype(np.float32)\n\ndrop_cols = [\"ID\",\"target_id\",\"resid\"]\nfeature_columns = [c for c in train_features.columns if c not in drop_cols]\n\nX = train_features[feature_columns].to_numpy(np.float32)\nX = impute_X(X)\n\ny_x = np.asarray(y_x, dtype=np.float32)\ny_y = np.asarray(y_y, dtype=np.float32)\ny_z = np.asarray(y_z, dtype=np.float32)\n\nprint(\"X:\", X.shape, \"nan:\", np.isnan(X).sum(), \"inf:\", np.isinf(X).sum())\nprint(\"n_features:\", len(feature_columns))","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:43:57.844665Z","iopub.status.busy":"2026-02-23T12:43:57.844308Z","iopub.status.idle":"2026-02-23T12:43:57.944911Z","shell.execute_reply":"2026-02-23T12:43:57.943681Z"},"papermill":{"duration":0.140612,"end_time":"2026-02-23T12:43:57.946928","exception":false,"start_time":"2026-02-23T12:43:57.806316","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 14) Train/Valid split (simple: random residue-level split)\n#     NOTE: use target_id-level split to avoid data leakage\n# =============\nfrom sklearn.model_selection import train_test_split\n\nidx = np.arange(X.shape[0])\ntr_idx, va_idx = train_test_split(idx, test_size=0.15, random_state=SEED)\n\nX_tr, X_va = X[tr_idx], X[va_idx]\nyx_tr, yx_va = y_x[tr_idx], y_x[va_idx]\nyy_tr, yy_va = y_y[tr_idx], y_y[va_idx]\nyz_tr, yz_va = y_z[tr_idx], y_z[va_idx]\n\ndtr_x = xgb.DMatrix(X_tr, label=yx_tr)\ndva_x = xgb.DMatrix(X_va, label=yx_va)\n\ndtr_y = xgb.DMatrix(X_tr, label=yy_tr)\ndva_y = xgb.DMatrix(X_va, label=yy_va)\n\ndtr_z = xgb.DMatrix(X_tr, label=yz_tr)\ndva_z = xgb.DMatrix(X_va, label=yz_va)\n\nprint(\"split:\", X_tr.shape, X_va.shape)","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:43:58.019962Z","iopub.status.busy":"2026-02-23T12:43:58.019611Z","iopub.status.idle":"2026-02-23T12:43:58.171461Z","shell.execute_reply":"2026-02-23T12:43:58.169641Z"},"papermill":{"duration":0.191114,"end_time":"2026-02-23T12:43:58.173655","exception":false,"start_time":"2026-02-23T12:43:57.982541","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 15) Train boosters (xgb.train + early stopping)\n# =============\nparams = {\n    \"tree_method\": \"hist\",\n    \"objective\": \"reg:squarederror\",\n    \"eval_metric\": \"rmse\",\n    \"eta\": 0.03,\n    \"max_depth\": 10,\n    \"subsample\": 0.75,\n    \"colsample_bytree\": 0.65,\n    \"min_child_weight\": 5,\n    \"gamma\": 0.1,\n    \"lambda\": 1.0,\n    \"alpha\": 0.1,\n    \"seed\": SEED,\n}\n\ndef train_booster(dtr, dva, name=\"x\"):\n    booster = xgb.train(\n        params=params,\n        dtrain=dtr,\n        num_boost_round=NUM_ROUND_MAX,\n        evals=[(dtr,\"train\"), (dva,\"valid\")],\n        early_stopping_rounds=EARLY_STOP,\n        verbose_eval=200\n    )\n    print(f\"[{name}] best_iteration =\", booster.best_iteration, \"best_score =\", booster.best_score)\n    return booster\n\nbx = train_booster(dtr_x, dva_x, \"x\")\nby = train_booster(dtr_y, dva_y, \"y\")\nbz = train_booster(dtr_z, dva_z, \"z\")","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:43:58.249084Z","iopub.status.busy":"2026-02-23T12:43:58.248668Z","iopub.status.idle":"2026-02-23T12:44:17.385849Z","shell.execute_reply":"2026-02-23T12:44:17.384141Z"},"papermill":{"duration":19.178269,"end_time":"2026-02-23T12:44:17.388136","exception":false,"start_time":"2026-02-23T12:43:58.209867","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 16) Build test features (same pipeline)\n# =============\ndef build_features_for_sequences(seqs_df: pd.DataFrame):\n    rows_all = []\n    for tid, seq in tqdm(zip(seqs_df[\"target_id\"].values, seqs_df[\"sequence\"].values),\n                         total=len(seqs_df), desc=\"Build test features\"):\n        rows = build_rows_for_target(str(tid), str(seq))\n        rows_all.extend(rows)\n        del rows\n        gc.collect()\n    return pd.DataFrame(rows_all)\n\ntest_features = build_features_for_sequences(test_sequences)\nprint(\"test_features:\", test_features.shape)\ntest_features.head()","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:44:17.464244Z","iopub.status.busy":"2026-02-23T12:44:17.463907Z","iopub.status.idle":"2026-02-23T12:45:44.578338Z","shell.execute_reply":"2026-02-23T12:45:44.57705Z"},"papermill":{"duration":87.154793,"end_time":"2026-02-23T12:45:44.580513","exception":false,"start_time":"2026-02-23T12:44:17.42572","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 17) Predict (x_1,y_1,z_1) and create decoys (x_2..z_5)\n# =============\nX_test = (\n    test_features\n    .reindex(columns=([\"ID\",\"target_id\",\"resid\"] + feature_columns), fill_value=0.0)\n)\n\nX_test_model = X_test[feature_columns].to_numpy(np.float32)\nX_test_model = impute_X(X_test_model)\ndtest = xgb.DMatrix(X_test_model)\n\npred_x = bx.predict(dtest)\npred_y = by.predict(dtest)\npred_z = bz.predict(dtest)\n\ntest_features = test_features.copy()\ntest_features[\"x_1\"] = pred_x\ntest_features[\"y_1\"] = pred_y\ntest_features[\"z_1\"] = pred_z\n\nfor k in range(2, 6):\n    noise_scale = 0.5\n    test_features[f\"x_{k}\"] = pred_x + rng.normal(0, noise_scale, size=len(pred_x))\n    test_features[f\"y_{k}\"] = pred_y + rng.normal(0, noise_scale, size=len(pred_y))\n    test_features[f\"z_{k}\"] = pred_z + rng.normal(0, noise_scale, size=len(pred_z))\n\nprint(\"Pred done:\", len(test_features))","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:45:44.666113Z","iopub.status.busy":"2026-02-23T12:45:44.665731Z","iopub.status.idle":"2026-02-23T12:45:45.102014Z","shell.execute_reply":"2026-02-23T12:45:45.101018Z"},"papermill":{"duration":0.485379,"end_time":"2026-02-23T12:45:45.10452","exception":false,"start_time":"2026-02-23T12:45:44.619141","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =============\n# 18) Build submission.csv (robust, fixed)\n# =============\nsub = sample_sub.copy()\n\n# prediction columns (excluding ID)\npred_cols = (\n    [f\"x_{k}\" for k in range(1,6)] +\n    [f\"y_{k}\" for k in range(1,6)] +\n    [f\"z_{k}\" for k in range(1,6)]\n)\n\n# ensure required columns exist in test_features (fill missing with 0)\nneed_cols = [\"ID\"] + pred_cols\nfor c in need_cols:\n    if c not in test_features.columns:\n        test_features[c] = 0.0\n\n# drop prediction columns from sample_sub if they already exist (keep ID!)\nsub = sub.drop(columns=[c for c in pred_cols if c in sub.columns], errors=\"ignore\")\n\n# merge\nsub = sub.merge(test_features[need_cols], on=\"ID\", how=\"left\")\n\n# reorder columns to match sample_submission\nsub = sub.reindex(columns=sample_sub.columns)\n\nsub.to_csv(\"submission.csv\", index=False)\nprint(\"saved: submission.csv\", sub.shape)\nsub.head()\n","metadata":{"execution":{"iopub.execute_input":"2026-02-23T12:45:45.183797Z","iopub.status.busy":"2026-02-23T12:45:45.183418Z","iopub.status.idle":"2026-02-23T12:45:45.471031Z","shell.execute_reply":"2026-02-23T12:45:45.470196Z"},"papermill":{"duration":0.330063,"end_time":"2026-02-23T12:45:45.473065","exception":false,"start_time":"2026-02-23T12:45:45.143002","status":"completed"},"tags":[]},"outputs":[],"execution_count":null}]}