{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":118765,"databundleVersionId":15231210,"sourceType":"competition"},{"sourceId":14604295,"sourceType":"datasetVersion","datasetId":9328538}],"dockerImageVersionId":31260,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true},"papermill":{"default_parameters":{},"duration":2354.045049,"end_time":"2026-01-13T06:54:44.524542","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2026-01-13T06:15:30.479493","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install --no-index /kaggle/input/biopython-cp312/biopython-1.86-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-01T07:09:56.53566Z","iopub.execute_input":"2026-02-01T07:09:56.536216Z","iopub.status.idle":"2026-02-01T07:10:01.884273Z","shell.execute_reply.started":"2026-02-01T07:09:56.536188Z","shell.execute_reply":"2026-02-01T07:10:01.883411Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os, sys, time, warnings\nimport numpy as np\nimport pandas as pd\nwarnings.filterwarnings(\"ignore\")\n\nDATA_PATH = \"/kaggle/input/stanford-rna-3d-folding-2/\"\nTRAIN_SEQS_CSV = os.path.join(DATA_PATH, \"train_sequences.csv\")\nTEST_SEQS_CSV  = os.path.join(DATA_PATH, \"test_sequences.csv\")\nTRAIN_LABELS_CSV = os.path.join(DATA_PATH, \"train_labels.csv\")\n\n# ---------- toggles ----------\nDEBUG = True\nDEBUG_PRINT_TOPK = 8\nDEBUG_EVERY_N = 1\nSAVE_DEBUG_CSV = True\nSAVE_DEBUG_MODELS_CSV = True\n\n# Candidate settings\nTOPN_POOL = 30\nLEN_RATIO_MAX = 0.30\nK_DIVERSITY = 12\n\n# Constraints\nCONSTRAINT_PASSES = 2\n\n# Output clip\nCLIP_MIN, CLIP_MAX = -999.999, 9999.999\n\n# RNG seed base\nGLOBAL_SEED = 20260207\n\n# Special handling for very long RNAs\nLONG_SEQ_THRESHOLD = 800  # >800 disable heavy transforms\n\nBOND_TARGET = 5.95\n\n# ============================================================\n# Robust import for Kaggle extra/parse_fasta_py.py\n# ============================================================\nsys.path.append(os.path.join(DATA_PATH, \"extra\"))\ntry:\n    import typing as _typing\n    import builtins as _builtins\n    _builtins.Dict  = getattr(_typing, \"Dict\")\n    _builtins.Tuple = getattr(_typing, \"Tuple\")\n    _builtins.List  = getattr(_typing, \"List\")\n    from parse_fasta_py import parse_fasta as _parse_fasta_raw\n    def parse_fasta(fasta_content: str):\n        d = _parse_fasta_raw(fasta_content)\n        out = {}\n        for k, v in d.items():\n            out[k] = v[0] if isinstance(v, tuple) else v\n        return out\nexcept Exception:\n    def parse_fasta(fasta_content: str):\n        out = {}\n        cur = None\n        seq_parts = []\n        for line in str(fasta_content).splitlines():\n            line = line.strip()\n            if not line:\n                continue\n            if line.startswith(\">\"):\n                if cur is not None:\n                    out[cur] = \"\".join(seq_parts)\n                cur = line[1:].split()[0]\n                seq_parts = []\n            else:\n                seq_parts.append(line.replace(\" \", \"\"))\n        if cur is not None:\n            out[cur] = \"\".join(seq_parts)\n        return out\n\ndef parse_stoichiometry(stoich: str):\n    if pd.isna(stoich) or str(stoich).strip() == \"\":\n        return []\n    out = []\n    for part in str(stoich).split(\";\"):\n        ch, cnt = part.split(\":\")\n        out.append((ch.strip(), int(cnt)))\n    return out\n\ndef get_chain_segments(row):\n    seq = row[\"sequence\"]\n    stoich = row.get(\"stoichiometry\", \"\")\n    all_seq = row.get(\"all_sequences\", \"\")\n\n    if pd.isna(stoich) or pd.isna(all_seq) or str(stoich).strip()==\"\" or str(all_seq).strip()==\"\":\n        return [(0, len(seq))]\n\n    try:\n        chain_dict = parse_fasta(all_seq)\n        order = parse_stoichiometry(stoich)\n        segs, pos = [], 0\n        for ch, cnt in order:\n            base = chain_dict.get(ch)\n            if base is None:\n                return [(0, len(seq))]\n            for _ in range(cnt):\n                L = len(base)\n                segs.append((pos, pos+L))\n                pos += L\n        if pos != len(seq):\n            return [(0, len(seq))]\n        return segs\n    except Exception:\n        return [(0, len(seq))]\n\ndef build_segments_map(df):\n    return {r[\"target_id\"]: get_chain_segments(r) for _, r in df.iterrows()}\n\n# ============================================================\n# Load data\n# ============================================================\ntrain_seqs = pd.read_csv(TRAIN_SEQS_CSV)\ntest_seqs  = pd.read_csv(TEST_SEQS_CSV)\ntrain_labels = pd.read_csv(TRAIN_LABELS_CSV)\n\nprint(\"Train sequences:\", len(train_seqs))\nprint(\"Test sequences:\", len(test_seqs))\nprint(\"Train labels rows:\", len(train_labels))\n\ntrain_segs_map = build_segments_map(train_seqs)\ntest_segs_map  = build_segments_map(test_seqs)\n\ndef process_labels(labels_df):\n    coords_dict = {}\n    prefixes = labels_df[\"ID\"].str.rsplit(\"_\", n=1).str[0]\n    for id_prefix, group in labels_df.groupby(prefixes):\n        coords_dict[id_prefix] = group.sort_values(\"resid\")[[\"x_1\",\"y_1\",\"z_1\"]].values.astype(np.float32)\n    return coords_dict\n\ntrain_coords_dict = process_labels(train_labels)\nprint(\"Train templates with coords:\", len(train_coords_dict))\n\ntrain_ids = []\ntrain_seq_list = []\ntrain_len = []\ntrain_coords_list = []\nfor _, r in train_seqs.iterrows():\n    tid = r[\"target_id\"]\n    if tid not in train_coords_dict:\n        continue\n    s = r[\"sequence\"]\n    train_ids.append(tid)\n    train_seq_list.append(s)\n    train_len.append(len(s))\n    train_coords_list.append(train_coords_dict[tid])\n\ntrain_len = np.array(train_len, dtype=np.int32)\nprint(\"Usable train templates:\", len(train_ids))\n\n# ============================================================\n# PairwiseAligner setup (anti-sliding)\n# ============================================================\nfrom Bio.Align import PairwiseAligner\n\naligner = PairwiseAligner()\naligner.mode = \"global\"\naligner.match_score = 2\naligner.mismatch_score = -1.5\naligner.open_gap_score = -8\naligner.extend_gap_score = -0.4\n\naligner.query_left_open_gap_score     = -8\naligner.query_left_extend_gap_score   = -0.4\naligner.query_right_open_gap_score    = -8\naligner.query_right_extend_gap_score  = -0.4\naligner.target_left_open_gap_score    = -8\naligner.target_left_extend_gap_score  = -0.4\naligner.target_right_open_gap_score   = -8\naligner.target_right_extend_gap_score = -0.4\n\n# ============================================================\n# Alignment stats\n# ============================================================\ndef alignment_stats(aln, qlen, tlen):\n    q_blocks, t_blocks = aln.aligned\n    if len(q_blocks) == 0:\n        return {\"aligned_len\":0,\"coverage\":0.0,\"terminal_gap\":qlen+tlen,\"blocks\":0,\"longest_block\":0}\n    aligned_len = int(sum(qe-qs for qs, qe in q_blocks))\n    coverage = aligned_len / max(1, qlen)\n    q_start = int(q_blocks[0][0]); q_end = int(q_blocks[-1][1])\n    t_start = int(t_blocks[0][0]); t_end = int(t_blocks[-1][1])\n    terminal_gap = (q_start + (qlen - q_end) + t_start + (tlen - t_end))\n    longest_block = int(max(qe-qs for qs, qe in q_blocks))\n    return {\n        \"aligned_len\": aligned_len,\n        \"coverage\": float(coverage),\n        \"terminal_gap\": int(terminal_gap),\n        \"blocks\": int(len(q_blocks)),\n        \"longest_block\": int(longest_block),\n    }\n\n# ============================================================\n# Template search (fast scan + topN)\n# ============================================================\ndef find_similar_sequences(query_seq, top_n=TOPN_POOL):\n    qlen = len(query_seq)\n    ratios = np.abs(train_len - qlen) / np.maximum(train_len, qlen)\n    cand_idx = np.where(ratios <= LEN_RATIO_MAX)[0]\n    if len(cand_idx) == 0:\n        return []\n    scored = []\n    for idx in cand_idx:\n        t_seq = train_seq_list[idx]\n        raw = aligner.score(query_seq, t_seq)  # fast score\n        norm = raw / (2.0 * max(1, min(qlen, len(t_seq))))\n        scored.append((idx, float(norm)))\n    scored.sort(key=lambda x: x[1], reverse=True)\n    out = []\n    for idx, sim in scored[:top_n]:\n        out.append((train_ids[idx], train_seq_list[idx], sim, train_coords_list[idx], idx))\n    return out\n\n# ============================================================\n# NEW: rerank top templates precisely (prevents bad choices)\n# ============================================================\ndef rerank_top_templates(query_seq, cands, topM=8):\n    M = min(topM, len(cands))\n    if M <= 1:\n        return cands\n    scored2 = []\n    for i in range(M):\n        t_id, t_seq, sim, t_coords, idx = cands[i]\n        aln = next(iter(aligner.align(query_seq, t_seq)))\n        st = alignment_stats(aln, len(query_seq), len(t_seq))\n        cov = st[\"coverage\"]\n        est_nan = max(0.0, 1.0 - cov)\n        combo = (\n            1.00*float(sim)\n            + 0.35*float(cov)\n            - 0.004*float(st[\"terminal_gap\"])\n            - 0.20*float(est_nan)\n        )\n        scored2.append((combo, i))\n    scored2.sort(key=lambda x: x[0], reverse=True)\n    return [cands[j] for _, j in scored2] + cands[M:]\n\n# ============================================================\n# Geometry helpers\n# ============================================================\ndef unit(v):\n    n = float(np.linalg.norm(v))\n    if n < 1e-8:\n        return None\n    return v / n\n\ndef random_walk_backbone(L, rng, step=BOND_TARGET, bend=0.35, jitter=0.15):\n    X = np.zeros((L, 3), dtype=np.float32)\n    d = rng.normal(size=3).astype(np.float32)\n    d = d / (np.linalg.norm(d) + 1e-8)\n    for i in range(1, L):\n        d = d + (bend * rng.normal(size=3)).astype(np.float32)\n        d = d / (np.linalg.norm(d) + 1e-8)\n        X[i] = X[i-1] + d * float(step)\n    if jitter > 0:\n        X += rng.normal(0, jitter, size=X.shape).astype(np.float32)\n    return X\n\ndef fill_missing_like_rna(new_coords, rng, step=BOND_TARGET):\n    X = new_coords.copy().astype(np.float32)\n    L = X.shape[0]\n    ok = np.isfinite(X[:, 0])\n    nan_ratio = 1.0 - (ok.sum() / max(1, L))\n\n    if ok.sum() == 0:\n        return random_walk_backbone(L, rng, step=step, jitter=0.15), nan_ratio\n\n    first = int(np.where(ok)[0][0])\n\n    # Fill left side\n    for i in range(first-1, -1, -1):\n        if i+2 < L and np.isfinite(X[i+2,0]):\n            d = unit(X[i+1] - X[i+2])\n        else:\n            d = unit(rng.normal(size=3))\n        if d is None:\n            d = np.array([1.0,0,0], dtype=np.float32)\n        X[i] = X[i+1] + d.astype(np.float32) * float(step)\n\n    # Fill middle + right\n    for i in range(first+1, L):\n        if np.isfinite(X[i,0]):\n            continue\n        prev = i-1\n        while prev >= 0 and not np.isfinite(X[prev,0]):\n            prev -= 1\n        nxt = i+1\n        while nxt < L and not np.isfinite(X[nxt,0]):\n            nxt += 1\n\n        if prev >= 0 and nxt < L:\n            w = (i - prev) / (nxt - prev)\n            X[i] = (1-w)*X[prev] + w*X[nxt]\n        elif prev >= 0:\n            if prev-1 >= 0 and np.isfinite(X[prev-1,0]):\n                d = unit(X[prev] - X[prev-1])\n            else:\n                d = unit(rng.normal(size=3))\n            if d is None:\n                d = np.array([1.0,0,0], dtype=np.float32)\n            X[i] = X[prev] + d.astype(np.float32) * float(step)\n        elif nxt < L:\n            if nxt+1 < L and np.isfinite(X[nxt+1,0]):\n                d = unit(X[nxt] - X[nxt+1])\n            else:\n                d = unit(rng.normal(size=3))\n            if d is None:\n                d = np.array([1.0,0,0], dtype=np.float32)\n            X[i] = X[nxt] + d.astype(np.float32) * float(step)\n        else:\n            X[i] = X[i-1] + np.array([step,0,0], dtype=np.float32)\n\n    return np.nan_to_num(X).astype(np.float32), nan_ratio\n\ndef clip_xyz(arr):\n    return np.clip(arr, CLIP_MIN, CLIP_MAX).astype(np.float32)\n\ndef centerize(X):\n    return (X - X.mean(axis=0, keepdims=True)).astype(np.float32)\n\n# ============================================================\n# Coordinate transfer via alignment.aligned\n# ============================================================\ndef adapt_template_to_query(query_seq, template_seq, template_coords, rng):\n    aln = next(iter(aligner.align(query_seq, template_seq)))\n    new_coords = np.full((len(query_seq), 3), np.nan, dtype=np.float32)\n\n    for (qs, qe), (ts, te) in zip(*aln.aligned):\n        ts = int(ts); te = int(te)\n        qs = int(qs); qe = int(qe)\n        chunk = template_coords[ts:te]\n        if len(chunk) == (qe - qs):\n            new_coords[qs:qe] = chunk\n\n    filled, nan_ratio = fill_missing_like_rna(new_coords, rng, step=BOND_TARGET)\n    return filled, aln, nan_ratio\n\n# ============================================================\n# NEW: kill catastrophic jumps (stepMax explosions)\n# ============================================================\ndef clamp_and_equalize_bonds(coords, target_id, step=BOND_TARGET, max_factor=2.2, passes=2):\n    X = coords.astype(np.float32).copy()\n    segments = test_segs_map.get(target_id, [(0, len(X))])\n\n    for _ in range(passes):\n        for (s, e) in segments:\n            if e - s < 2:\n                continue\n            Y = X[s:e]\n\n            # clamp huge jumps\n            d = Y[1:] - Y[:-1]\n            dist = np.linalg.norm(d, axis=1) + 1e-8\n            bad = dist > (max_factor * step)\n            if np.any(bad):\n                for i in np.where(bad)[0]:\n                    u = d[i] / dist[i]\n                    Y[i+1] = Y[i] + u * (max_factor * step)\n\n            # equalize forward\n            for i in range(1, len(Y)):\n                v = Y[i] - Y[i-1]\n                dv = float(np.linalg.norm(v)) + 1e-8\n                u = v / dv\n                Y[i] = Y[i-1] + u * float(step)\n\n            # light smoothing\n            if len(Y) >= 4:\n                Z = Y.copy()\n                Z[1:-1] = 0.25*Y[:-2] + 0.5*Y[1:-1] + 0.25*Y[2:]\n                Y = Z\n\n            X[s:e] = Y\n\n    return X.astype(np.float32)\n\n# ============================================================\n# Constraints: segment-aware (FIXED strength schedule)\n# ============================================================\ndef adaptive_rna_constraints(coordinates, target_id, confidence=1.0, passes=CONSTRAINT_PASSES, nan_ratio=0.0):\n    coords = coordinates.astype(np.float32).copy()\n    segments = test_segs_map.get(target_id, [(0, len(coords))])\n\n    conf = float(confidence)\n    nr = float(nan_ratio)\n\n    # near 0 if conf high\n    strength = 0.85 * max(0.0, (0.80 - min(conf, 0.80)) / 0.80)  # conf>=0.8 => 0\n    # more if gaps filled\n    strength += 0.35 * min(1.0, nr / 0.35)\n    strength = float(np.clip(strength, 0.0, 0.60))\n\n    # very confident & little gaps: don't touch\n    if conf >= 0.85 and nr <= 0.08:\n        return coords\n\n    for _ in range(passes):\n        for (s, e) in segments:\n            X = coords[s:e]\n            L = e - s\n            if L < 3:\n                coords[s:e] = X\n                continue\n\n            # bond i,i+1\n            d = X[1:] - X[:-1]\n            dist = np.linalg.norm(d, axis=1) + 1e-6\n            scale = (BOND_TARGET - dist) / dist\n            adj = (d * scale[:, None]) * (0.22 * strength)\n            X[:-1] -= adj\n            X[1:]  += adj\n\n            # i,i+2\n            d2 = X[2:] - X[:-2]\n            dist2 = np.linalg.norm(d2, axis=1) + 1e-6\n            scale2 = (10.2 - dist2) / dist2\n            adj2 = (d2 * scale2[:, None]) * (0.10 * strength)\n            X[:-2] -= adj2\n            X[2:]  += adj2\n\n            # smoothing\n            lap = 0.5 * (X[:-2] + X[2:]) - X[1:-1]\n            X[1:-1] += (0.05 * strength) * lap\n\n            # self-avoidance only when low conf / high gaps\n            if (conf < 0.35 or nr > 0.20) and L >= 25:\n                k = min(L, 160) if L > 220 else L\n                idx = np.linspace(0, L - 1, k).astype(int) if k < L else np.arange(L)\n                P = X[idx]\n                diff = P[:, None, :] - P[None, :, :]\n                distm = np.linalg.norm(diff, axis=2) + 1e-6\n                sep = np.abs(idx[:, None] - idx[None, :])\n                mask = (sep > 2) & (distm < 3.2)\n                if np.any(mask):\n                    force = (3.2 - distm) / distm\n                    vec = (diff * force[:, :, None] * mask[:, :, None]).sum(axis=1)\n                    X[idx] += (0.012 * strength) * vec\n\n            coords[s:e] = X\n\n    return coords\n\n# ============================================================\n# Diversity transforms + rotation helpers\n# ============================================================\ndef _rotmat(axis, ang):\n    axis = np.asarray(axis, float)\n    axis = axis / (np.linalg.norm(axis) + 1e-12)\n    x, y, z = axis\n    c, s = np.cos(ang), np.sin(ang)\n    C = 1.0 - c\n    return np.array([\n        [c + x*x*C,     x*y*C - z*s, x*z*C + y*s],\n        [y*x*C + z*s,   c + y*y*C,   y*z*C - x*s],\n        [z*x*C - y*s,   z*y*C + x*s, c + z*z*C]\n    ], dtype=float)\n\ndef apply_hinge(coords, seg, rng, max_angle_deg=16):\n    s, e = seg\n    L = e - s\n    if L < 30:\n        return coords\n    pivot = s + int(rng.integers(10, L - 10))\n    axis = rng.normal(size=3)\n    ang = np.deg2rad(float(rng.uniform(-max_angle_deg, max_angle_deg)))\n    R = _rotmat(axis, ang)\n    X = coords.copy()\n    p0 = X[pivot].copy()\n    X[pivot+1:e] = (X[pivot+1:e] - p0) @ R.T + p0\n    return X\n\ndef jitter_chains(coords, segments, rng, max_angle_deg=7, max_trans=0.7):\n    X = coords.copy()\n    global_center = X.mean(axis=0, keepdims=True)\n    for (s, e) in segments:\n        axis = rng.normal(size=3)\n        ang = np.deg2rad(float(rng.uniform(-max_angle_deg, max_angle_deg)))\n        R = _rotmat(axis, ang)\n        shift = rng.normal(size=3)\n        shift = shift / (np.linalg.norm(shift) + 1e-12) * float(rng.uniform(0.0, max_trans))\n        c = X[s:e].mean(axis=0, keepdims=True)\n        X[s:e] = (X[s:e] - c) @ R.T + c + shift\n    X -= X.mean(axis=0, keepdims=True) - global_center\n    return X\n\ndef smooth_wiggle(coords, segments, rng, amp=0.55):\n    X = coords.copy()\n    for (s, e) in segments:\n        L = e - s\n        if L < 20:\n            continue\n        n_ctrl = 6\n        ctrl_x = np.linspace(0, L - 1, n_ctrl)\n        ctrl_disp = rng.normal(0, amp, size=(n_ctrl, 3))\n        t = np.arange(L)\n        disp = np.vstack([np.interp(t, ctrl_x, ctrl_disp[:, k]) for k in range(3)]).T\n        X[s:e] += disp\n    return X\n\n# ============================================================\n# Diagnostics\n# ============================================================\ndef bond_mae(X):\n    d = X[1:] - X[:-1]\n    dist = np.linalg.norm(d, axis=1)\n    return float(np.mean(np.abs(dist - BOND_TARGET))) if len(dist) else 0.0\n\ndef step_max(X):\n    d = X[1:] - X[:-1]\n    dist = np.linalg.norm(d, axis=1)\n    return float(np.max(dist)) if len(dist) else 0.0\n\ndef radius_median(X):\n    c = X.mean(axis=0, keepdims=True)\n    r = np.linalg.norm(X - c, axis=1)\n    return float(np.median(r)) if len(r) else 0.0\n\n# ============================================================\n# NEW: diversity policy (protect hi-conf templates)\n# ============================================================\ndef choose_diverse_model(adapted, segments, rng, sim, L, i, nan_ratio):\n    conf = float(sim)\n    nr = float(nan_ratio)\n\n    # high-confidence: DO NOT hinge/jitter/wiggle\n    if conf >= 0.85 and nr <= 0.10:\n        sigma = 0.004 + 0.0015*i\n        return adapted + rng.normal(0, sigma, adapted.shape).astype(np.float32), f\"hi_conf_tiny_noise={sigma:.4f}\"\n\n    # mid confidence: small noise only\n    if conf >= 0.55 and nr <= 0.20:\n        if i == 0:\n            return adapted, \"base\"\n        sigma = 0.010 + 0.004*i\n        return adapted + rng.normal(0, sigma, adapted.shape).astype(np.float32), f\"mid_conf_noise={sigma:.3f}\"\n\n    # long RNAs: noise only\n    if L > LONG_SEQ_THRESHOLD:\n        if i == 0:\n            return adapted, \"base\"\n        sigma = 0.012 + 0.003*i\n        return adapted + rng.normal(0, sigma, adapted.shape).astype(np.float32), f\"long_noise={sigma:.3f}\"\n\n    # many gaps: avoid hinge/jitter; noise + mild wiggle\n    if nr > 0.35:\n        if i == 0:\n            return adapted, \"base\"\n        if i <= 2:\n            sigma = 0.018 + 0.004*i\n            return adapted + rng.normal(0, sigma, adapted.shape).astype(np.float32), f\"gap_noise={sigma:.3f}\"\n        X = smooth_wiggle(adapted, segments, rng, amp=0.45).astype(np.float32)\n        return X, \"gap_mild_wiggle\"\n\n    # otherwise: mild version of old scheme\n    if i == 0:\n        return adapted, \"base\"\n    if i == 1:\n        sigma = max(0.012, (0.42 - conf) * 0.05)\n        return adapted + rng.normal(0, sigma, adapted.shape).astype(np.float32), f\"noise={sigma:.3f}\"\n    if i == 2:\n        longest = max(segments, key=lambda se: se[1] - se[0])\n        return apply_hinge(adapted, longest, rng, max_angle_deg=16).astype(np.float32), \"hinge16\"\n    if i == 3:\n        return jitter_chains(adapted, segments, rng, max_angle_deg=7, max_trans=0.7).astype(np.float32), \"jitter_mild\"\n    return smooth_wiggle(adapted, segments, rng, amp=0.55).astype(np.float32), \"wiggle0.55\"\n\n# ============================================================\n# Prediction per target\n# ============================================================\ndef predict_rna_structures(row, test_index, n_predictions=5):\n    tid = row[\"target_id\"]\n    seq = row[\"sequence\"]\n    segments = test_segs_map.get(tid, [(0, len(seq))])\n    L = len(seq)\n\n    assert set(seq).issubset(set(\"ACGU\")), f\"Non-ACGU in {tid}; do not remap.\"\n\n    cands = find_similar_sequences(seq, top_n=TOPN_POOL)\n    cands = rerank_top_templates(seq, cands, topM=8)\n\n    dbg = {\n        \"target_id\": tid,\n        \"len\": L,\n        \"num_cands\": len(cands),\n        \"best_sim\": cands[0][2] if cands else 0.0,\n        \"best_template\": cands[0][0] if cands else \"\",\n        \"used_templates\": 0,\n    }\n\n    if DEBUG and (test_index % DEBUG_EVERY_N == 0):\n        print(f\"\\n[{tid}] L={L} segments={len(segments)} candidates={len(cands)}\")\n        if cands:\n            print(\"  Top candidates (id, sim):\")\n            for k in range(min(DEBUG_PRINT_TOPK, len(cands))):\n                print(f\"    {k:02d} {cands[k][0]}  sim={cands[k][2]:.4f}  tlen={len(cands[k][1])}\")\n\n    predictions = []\n    model_debug_rows = []\n    used = set()\n\n    # fallback if no cands\n    if not cands:\n        for i in range(n_predictions):\n            seed = (abs(hash(tid)) + i * 10007 + GLOBAL_SEED) % (2**32)\n            rng = np.random.default_rng(seed)\n            X = random_walk_backbone(L, rng, step=BOND_TARGET, jitter=0.15 + 0.05*i)\n            X = centerize(X)\n            predictions.append(clip_xyz(X))\n            model_debug_rows.append({\n                \"target_id\": tid, \"model\": i+1,\n                \"template\": \"\", \"sim\": 0.0,\n                \"cov\": 0.0, \"blocks\": 0, \"terminal_gap\": 0,\n                \"nan_ratio\": 1.0,\n                \"bond_mae\": bond_mae(X),\n                \"step_max\": step_max(X),\n                \"radius_median\": radius_median(X),\n                \"note\": \"no_cands_randomwalk\"\n            })\n        return predictions, dbg, model_debug_rows\n\n    for i in range(n_predictions):\n        seed = (abs(hash(tid)) + i * 10007 + GLOBAL_SEED) % (2**32)\n        rng = np.random.default_rng(seed)\n\n        # choose template\n        if i == 0:\n            t_id, t_seq, sim, t_coords, _ = cands[0]\n        else:\n            K = min(K_DIVERSITY, len(cands))\n            sims = np.array([cands[k][2] for k in range(K)], float)\n            w = np.exp((sims - sims.max()) / 0.08)\n            for k in range(K):\n                if cands[k][0] in used:\n                    w[k] *= 0.10\n            w = w / (w.sum() + 1e-12)\n            k = int(rng.choice(np.arange(K), p=w))\n            t_id, t_seq, sim, t_coords, _ = cands[k]\n\n        used.add(t_id)\n\n        adapted, aln, nan_ratio = adapt_template_to_query(seq, t_seq, t_coords, rng)\n        stats = alignment_stats(aln, len(seq), len(t_seq))\n\n        if DEBUG and (test_index % DEBUG_EVERY_N == 0) and i == 0:\n            print(f\"  Best alignment stats: cov={stats['coverage']:.3f}, terminal_gap={stats['terminal_gap']}, \"\n                  f\"blocks={stats['blocks']}, longest_block={stats['longest_block']}\")\n\n        # diversity (fixed policy)\n        X, note = choose_diverse_model(adapted, segments, rng, sim, L, i, nan_ratio)\n\n        # constraints (fixed schedule)\n        refined = adaptive_rna_constraints(\n            X, tid, confidence=sim, passes=CONSTRAINT_PASSES, nan_ratio=nan_ratio\n        ).astype(np.float32)\n\n        # NEW: kill catastrophic jumps\n        refined = clamp_and_equalize_bonds(refined, tid, step=BOND_TARGET, max_factor=2.2, passes=2)\n\n        # centerize\n        refined = centerize(refined)\n        predictions.append(refined)\n\n        model_debug_rows.append({\n            \"target_id\": tid, \"model\": i+1,\n            \"template\": t_id, \"sim\": float(sim),\n            \"cov\": stats[\"coverage\"], \"blocks\": stats[\"blocks\"], \"terminal_gap\": stats[\"terminal_gap\"],\n            \"nan_ratio\": float(nan_ratio),\n            \"bond_mae\": bond_mae(refined),\n            \"step_max\": step_max(refined),\n            \"radius_median\": radius_median(refined),\n            \"note\": note\n        })\n\n    # align translations to model1 (translation only)\n    c0 = predictions[0].mean(axis=0, keepdims=True)\n    for i in range(1, len(predictions)):\n        ci = predictions[i].mean(axis=0, keepdims=True)\n        predictions[i] = predictions[i] - ci + c0\n\n    predictions = [clip_xyz(p) for p in predictions]\n    dbg[\"used_templates\"] = len(used)\n\n    if DEBUG and (test_index % DEBUG_EVERY_N == 0):\n        print(\"  Model diagnostics:\")\n        for r in model_debug_rows:\n            danger = (r[\"step_max\"] > 2.2*BOND_TARGET) or (r[\"bond_mae\"] > 2.5)\n            flag = \" !!!\" if danger else \"\"\n            print(f\"    m{r['model']} sim={r['sim']:.4f} cov={r['cov']:.3f} blocks={r['blocks']:3d} \"\n                  f\"nan={r['nan_ratio']:.2f} bondMAE={r['bond_mae']:.2f} stepMax={r['step_max']:.2f} \"\n                  f\"radMed={r['radius_median']:.1f} note={r['note']}{flag}\")\n\n    return predictions, dbg, model_debug_rows\n\n# ============================================================\n# Run + build submission\n# ============================================================\nall_predictions = []\ndebug_rows = []\ndebug_model_rows = []\n\nstart_time = time.time()\nfor test_index, row in enumerate(test_seqs.itertuples(index=False)):\n    row = row._asdict()\n    if test_index % 10 == 0:\n        print(f\"\\nProcessing {test_index}/{len(test_seqs)} | {time.time()-start_time:.1f}s\")\n\n    preds, dbg, mrows = predict_rna_structures(row, test_index, n_predictions=5)\n    debug_rows.append(dbg)\n    debug_model_rows.extend(mrows)\n\n    tid = row[\"target_id\"]\n    seq = row[\"sequence\"]\n\n    for j in range(len(seq)):\n        res = {\"ID\": f\"{tid}_{j+1}\", \"resname\": seq[j], \"resid\": j+1}\n        for i in range(5):\n            res[f\"x_{i+1}\"], res[f\"y_{i+1}\"], res[f\"z_{i+1}\"] = map(float, preds[i][j])\n        all_predictions.append(res)\n\nsub = pd.DataFrame(all_predictions)\ncols = [\"ID\", \"resname\", \"resid\"] + [f\"{c}_{i}\" for i in range(1, 6) for c in [\"x\", \"y\", \"z\"]]\n\ncoord_cols = [c for c in cols if c.startswith((\"x_\",\"y_\",\"z_\"))]\nsub[coord_cols] = sub[coord_cols].clip(CLIP_MIN, CLIP_MAX)\n\nsub[cols].to_csv(\"submission.csv\", index=False)\nprint(\"\\nsubmission.csv saved!\")\n\ndbg_df = pd.DataFrame(debug_rows)\nif SAVE_DEBUG_CSV:\n    dbg_df.to_csv(\"debug_targets.csv\", index=False)\n    print(\"debug_targets.csv saved!\")\n\nif SAVE_DEBUG_MODELS_CSV:\n    pd.DataFrame(debug_model_rows).to_csv(\"debug_models.csv\", index=False)\n    print(\"debug_models.csv saved!\")\n\nprint(\"\\n=== Debug summary ===\")\nprint(dbg_df.describe(include=\"all\"))\nprint(\"\\nHead of submission:\")\nprint(sub.head())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-01T07:10:01.886103Z","iopub.execute_input":"2026-02-01T07:10:01.886358Z","iopub.status.idle":"2026-02-01T07:12:28.114333Z","shell.execute_reply.started":"2026-02-01T07:10:01.88633Z","shell.execute_reply":"2026-02-01T07:12:28.113707Z"}},"outputs":[],"execution_count":null}]}