{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceType":"competition","sourceId":118765,"databundleVersionId":15231210,"isSourceIdPinned":false}],"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"### Reference\n[RNA 3D Folding Template-Based (0.321)](https://www.kaggle.com/code/jaejohn/rna-3d-folding-template-based-0-321)\n","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n!pip install biopython\nfrom Bio import BiopythonWarning\nfrom tqdm import tqdm\nimport warnings\nimport random\nwarnings.simplefilter('ignore', BiopythonWarning)\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nfrom collections import defaultdict\n\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicki ng run or pressing Shift+Enter) will list all files under the input directory\n  \nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:16:20.455982Z","iopub.execute_input":"2026-02-21T09:16:20.456276Z","iopub.status.idle":"2026-02-21T09:16:32.37636Z","shell.execute_reply.started":"2026-02-21T09:16:20.456251Z","shell.execute_reply":"2026-02-21T09:16:32.375324Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 📦 Load core datasets\n# ============================================================\n\nDATA_PATH = \"/kaggle/input/stanford-rna-3d-folding-2\"\n\ndef load_csv(name):\n    return pd.read_csv(f\"{DATA_PATH}/{name}.csv\", low_memory=False)\n\nprint(\"📥 Loading datasets...\")\n\ntrain_seqs   = load_csv(\"train_sequences\")\nvalid_seqs   = load_csv(\"validation_sequences\")\ntest_seqs    = load_csv(\"test_sequences\")\ntrain_labels = load_csv(\"train_labels\")\nvalid_labels = load_csv(\"validation_labels\")\n\nprint(\"✅ Dataset loading complete\")\nprint(f\"   Train sequences:      {len(train_seqs)}\")\nprint(f\"   Validation sequences: {len(valid_seqs)}\")\nprint(f\"   Test sequences:       {len(test_seqs)}\")\nprint(f\"   Train labels:         {len(train_labels)}\")\nprint(f\"   Validation labels:    {len(valid_labels)}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:16:32.378146Z","iopub.execute_input":"2026-02-21T09:16:32.378466Z","iopub.status.idle":"2026-02-21T09:16:44.475317Z","shell.execute_reply.started":"2026-02-21T09:16:32.378442Z","shell.execute_reply":"2026-02-21T09:16:44.474573Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧬 Preprocess coords and collect mean and std\n# ============================================================\ndef preprocess_coords(coords):\n    coords = np.asarray(coords, float)\n    mean = coords.mean(axis=0, keepdims=True)\n    std  = coords.std(axis=0, keepdims=True) + 1e-8\n    return (coords - mean) / std","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:16:44.476296Z","iopub.execute_input":"2026-02-21T09:16:44.476534Z","iopub.status.idle":"2026-02-21T09:16:44.481245Z","shell.execute_reply.started":"2026-02-21T09:16:44.476486Z","shell.execute_reply":"2026-02-21T09:16:44.480694Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧬 Prepare label structures\n# ============================================================\n\n# Training labels: one coordinate triplet\ntrain_coord_cols = [\"x_1\", \"y_1\", \"z_1\"]\n\n# Validation labels: many coordinate triplets (x_1 ... x_40)\nvalid_coord_cols = [\n    col for col in valid_labels.columns\n    if col.startswith((\"x_\", \"y_\", \"z_\"))\n]\n\n# Submission requires x_1 ... x_5 for each axis\nsubmission_coord_cols = [\n    f\"{axis}_{i}\"\n    for i in range(1, 6)\n    for axis in [\"x\", \"y\", \"z\"]\n]\n\n# Convert numeric columns safely\ntrain_labels[train_coord_cols] = train_labels[train_coord_cols].apply(\n    pd.to_numeric, errors=\"coerce\"\n)\nvalid_labels[valid_coord_cols] = valid_labels[valid_coord_cols].apply(\n    pd.to_numeric, errors=\"coerce\"\n)\n\n# Extract 4-character ID prefix\ntrain_seqs[\"short_id\"] = train_seqs[\"target_id\"].str[:4]\n\n# Keep only labels whose ID exists in train_seqs\nvalid_ids = set(train_seqs[\"short_id\"])\ntrain_labels = (\n    train_labels[train_labels[\"ID\"].str[:4].isin(valid_ids)]\n    .reset_index(drop=True)\n)\n\nprint(\"🔧 Label preprocessing complete\")\nprint(f\"   Training coord columns:    {train_coord_cols}\")\nprint(f\"   Validation coord columns:  {len(valid_coord_cols)}\")\nprint(f\"   Filtered training labels:  {len(train_labels)}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:16:44.483262Z","iopub.execute_input":"2026-02-21T09:16:44.483484Z","iopub.status.idle":"2026-02-21T09:16:47.467591Z","shell.execute_reply.started":"2026-02-21T09:16:44.483464Z","shell.execute_reply":"2026-02-21T09:16:47.466834Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧩 Convert label DataFrames into coordinate dictionaries\n# ============================================================\n\ndef process_labels(labels_df, coord_cols):\n    \"\"\"\n    Convert a label DataFrame into:\n        { target_id : ndarray of coordinates }\n    \"\"\"\n    df = labels_df.copy()\n    df[\"target_id\"] = df[\"ID\"].str[:4]\n    df = df.sort_values([\"target_id\", \"resid\"])\n\n    coords_dict = {}\n    for target_id, idx in df.groupby(\"target_id\").groups.items():\n        coords = df.loc[idx, coord_cols].to_numpy()\n        coords_dict[target_id] = coords\n\n    return coords_dict\n\ntrain_coords_dict = process_labels(train_labels, train_coord_cols)\nvalid_coords_dict = process_labels(valid_labels, valid_coord_cols)\n\nprint(\"📐 Coordinate dictionaries created\")\nprint(f\"   Training structures:   {len(train_coords_dict)}\")\nprint(f\"   Validation structures: {len(valid_coords_dict)}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:16:47.468591Z","iopub.execute_input":"2026-02-21T09:16:47.468897Z","iopub.status.idle":"2026-02-21T09:22:17.493484Z","shell.execute_reply.started":"2026-02-21T09:16:47.468865Z","shell.execute_reply":"2026-02-21T09:22:17.492657Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compute_geom_stats(coords, max_coord=100.0):\n    \"\"\"\n    Compute biologically meaningful RNA geometry statistics using a\n    LOCAL helical axis and robust smoothing.\n\n    Returns:\n        twist_mean, twist_std\n        rise_mean, rise_std\n        radius_mean, radius_std\n        noise_xy_std, noise_z_std\n    \"\"\"\n\n    coords = np.asarray(coords, float)\n\n    # ------------------------------------------------------------\n    # 1. Filter invalid coordinates\n    # ------------------------------------------------------------\n    valid = np.isfinite(coords).all(axis=1)\n    valid &= (np.abs(coords) < max_coord).all(axis=1)\n    coords = coords[valid]\n\n    if coords.shape[0] < 4:\n        return {\n            \"twist_mean\": 0.0, \"twist_std\": 0.01,\n            \"rise_mean\": 2.8,  \"rise_std\": 0.1,\n            \"radius_mean\": 11.0, \"radius_std\": 1.0,\n            \"noise_xy_std\": 0.1, \"noise_z_std\": 0.05,\n        }\n\n    # ------------------------------------------------------------\n    # 2. Center coordinates\n    # ------------------------------------------------------------\n    center = coords.mean(axis=0)\n    coords_c = coords - center\n\n    # ------------------------------------------------------------\n    # 3. Compute LOCAL helix axis using sliding-window PCA\n    # ------------------------------------------------------------\n    N = coords_c.shape[0]\n    window = 8  # small window for local axis\n    axes = np.zeros((N, 3))\n\n    for i in range(N):\n        lo = max(0, i - window)\n        hi = min(N, i + window)\n        segment = coords_c[lo:hi]\n\n        # PCA on local segment\n        U, S, Vt = np.linalg.svd(segment - segment.mean(axis=0), full_matrices=False)\n        axis = Vt[0]\n        axis /= np.linalg.norm(axis)\n        axes[i] = axis\n\n    # Smooth axis field\n    axes = np.apply_along_axis(lambda v: v / (np.linalg.norm(v) + 1e-8), 1,\n                               np.convolve(axes[:,0], np.ones(5)/5, mode='same')[:,None] *\n                               np.array([1,0,0]) +\n                               np.convolve(axes[:,1], np.ones(5)/5, mode='same')[:,None] *\n                               np.array([0,1,0]) +\n                               np.convolve(axes[:,2], np.ones(5)/5, mode='same')[:,None] *\n                               np.array([0,0,1]))\n\n    # ------------------------------------------------------------\n    # 4. Project onto local axis → rise\n    # ------------------------------------------------------------\n    proj = np.sum(coords_c * axes, axis=1)\n    rise = np.diff(proj)\n\n    # ------------------------------------------------------------\n    # 5. Compute twist using rotation around local axis\n    # ------------------------------------------------------------\n    # Remove axis component → radial vectors\n    radial = coords_c - axes * proj[:, None]\n\n    # Compute angles in radial plane\n    angles = np.unwrap(np.arctan2(radial[:, 1], radial[:, 0]))\n    twist = np.diff(angles)\n\n    # ------------------------------------------------------------\n    # 6. Compute radius\n    # ------------------------------------------------------------\n    radius = np.linalg.norm(radial, axis=1)\n\n    # ------------------------------------------------------------\n    # 7. Compute noise (deviation from ideal helical step)\n    # ------------------------------------------------------------\n    step_vecs = np.diff(coords_c, axis=0)\n    ideal_steps = axes[:-1] * rise[:, None]\n    noise_vecs = step_vecs - ideal_steps\n\n    noise_xy = np.linalg.norm(noise_vecs[:, :2], axis=1)\n    noise_z = np.abs(noise_vecs[:, 2])\n\n    # ------------------------------------------------------------\n    # 8. Return stats\n    # ------------------------------------------------------------\n    return {\n        \"twist_mean\": float(np.mean(twist)),\n        \"twist_std\":  float(np.std(twist) + 1e-8),\n\n        \"rise_mean\":  float(np.mean(rise)),\n        \"rise_std\":   float(np.std(rise) + 1e-8),\n\n        \"radius_mean\": float(np.mean(radius)),\n        \"radius_std\":  float(np.std(radius) + 1e-8),\n\n        \"noise_xy_std\": float(np.mean(noise_xy)),\n        \"noise_z_std\":  float(np.mean(noise_z)),\n    }","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:22:17.494618Z","iopub.execute_input":"2026-02-21T09:22:17.49486Z","iopub.status.idle":"2026-02-21T09:22:17.507206Z","shell.execute_reply.started":"2026-02-21T09:22:17.494836Z","shell.execute_reply":"2026-02-21T09:22:17.506616Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧬 Simple monomer → base mapper (fast, universal)\n# ============================================================\n\ndef map_monomer_to_base(mon):\n    mon = mon.upper().strip()\n\n    # Empty or invalid\n    if len(mon) == 0 or mon in (\".\", \"?\", \"-\", \"NULL\"):\n        return \"N\"\n\n    # Direct RNA\n    if mon in (\"A\", \"C\", \"G\", \"U\"):\n        return mon\n\n    # DNA-style (DA, DC, DG, DT)\n    if mon.startswith(\"D\"):\n        if len(mon) >= 2:\n            base = mon[1]\n            return \"U\" if base == \"T\" else base\n        else:\n            return \"N\"\n\n    # Modified nucleotides (m5C → C, 2MG → G, etc.)\n    if mon[0] in (\"A\", \"C\", \"G\", \"U\"):\n        return mon[0]\n\n    # Chemistry-style (*GP → G)\n    if mon.endswith(\"P\") and len(mon) >= 2:\n        return mon[-2]\n\n    # Unknown → ignore for GC%\n    return \"N\"\n\n\n# ============================================================\n# 🧬 Extract sequence from CIF using _entity_poly_seq.mon_id\n# ============================================================\n\ndef extract_sequence_from_cif(filepath):\n    \"\"\"\n    Extracts sequences from mmCIF using _entity_poly_seq.mon_id.\n    Returns {chain_id: \"ACGU...\"}.\n    \"\"\"\n    seqs = {}\n    with open(filepath, \"r\") as f:\n        lines = f.readlines()\n\n    in_loop = False\n    cols = []\n    data_start = None\n\n    # Find the mon_id loop\n    for i, line in enumerate(lines):\n        line = line.strip()\n\n        if line == \"loop_\":\n            in_loop = True\n            cols = []\n            continue\n\n        if in_loop:\n            if line.startswith(\"_entity_poly_seq.\"):\n                cols.append(line)\n                continue\n            else:\n                if cols:\n                    data_start = i\n                    break\n                else:\n                    in_loop = False\n\n    if not cols:\n        return {}\n\n    col_idx = {c: i for i, c in enumerate(cols)}\n\n    # Required field\n    if \"_entity_poly_seq.mon_id\" not in col_idx:\n        return {}\n\n    chain_col = col_idx.get(\"_entity_poly_seq.entity_id\", None)\n    mon_col   = col_idx[\"_entity_poly_seq.mon_id\"]\n\n    # Parse rows\n    for line in lines[data_start:]:\n        if not line or line.startswith(\"#\"):\n            continue\n\n        parts = line.split()\n        if len(parts) < len(cols):\n            continue\n\n        chain_id = parts[chain_col] if chain_col is not None else \"A\"\n        mon_id   = parts[mon_col]\n\n        base = map_monomer_to_base(mon_id)\n        if base == \"N\":\n            continue\n\n        seqs.setdefault(chain_id, \"\")\n        seqs[chain_id] += base\n\n    return seqs\n\n\n# ============================================================\n# 📊 Compute dataset geometry + GC distribution (streaming)\n# ============================================================\n\ndef compute_dataset_stats_streaming(cif_dir, compute_local_geometry):\n    \"\"\"\n    Streams over all CIFs:\n      - extracts backbone coords (your loader)\n      - extracts sequences\n      - computes geometry stats\n      - computes GC distribution\n    \"\"\"\n\n    import numpy as np\n    import os\n    from tqdm import tqdm\n\n    # Accumulators for geometry\n    acc = {\n        \"radius_mean\": [],\n        \"radius_std\": [],\n        \"rise_mean\": [],\n        \"rise_std\": [],\n        \"twist_mean\": [],\n        \"twist_std\": [],\n        \"noise_xy_std\": [],\n        \"noise_z_std\": [],\n    }\n\n    # Accumulators for GC distribution\n    gc_values = []\n    seq_lengths = []\n\n    # Collect CIF files\n    cif_files = []\n    for root, dirs, files in os.walk(cif_dir):\n        for fname in files:\n            if fname.endswith(\".cif\"):\n                cif_files.append(os.path.join(root, fname))\n\n    # Stream through CIFs\n    for path in tqdm(cif_files, desc=\"Processing CIF files\"):\n\n        # --- Extract sequences ---\n        seqs = extract_sequence_from_cif(path)\n        for chain_id, seq in seqs.items():\n            if len(seq) >= 3:\n                gc = (seq.count(\"G\") + seq.count(\"C\")) / len(seq)\n                gc_values.append(gc)\n                seq_lengths.append(len(seq))\n\n        # --- Extract backbone coords ---\n        chain_dict = load_backbone_from_cif_safe(path)\n\n        for chain_id, coords_12 in chain_dict.items():\n            P_coords = coords_12[:, 0, :]  # P atom positions\n            geom = compute_local_geometry(P_coords)\n            if geom is None:\n                continue\n\n            for key in acc:\n                acc[key].append(geom.get(key, np.nan))\n\n    # Final geometry stats\n    geom_stats = {\n        key: float(np.nanmean(values)) if len(values) > 0 else 0.0\n        for key, values in acc.items()\n    }\n\n    # GC distribution stats\n    if len(gc_values) > 0:\n        geom_stats[\"gc_mean\"] = float(np.mean(gc_values))\n        geom_stats[\"gc_std\"]  = float(np.std(gc_values))\n    else:\n        geom_stats[\"gc_mean\"] = 0.5\n        geom_stats[\"gc_std\"]  = 0.1\n\n    return geom_stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:22:17.508061Z","iopub.execute_input":"2026-02-21T09:22:17.508286Z","iopub.status.idle":"2026-02-21T09:22:17.524938Z","shell.execute_reply.started":"2026-02-21T09:22:17.508256Z","shell.execute_reply":"2026-02-21T09:22:17.524318Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# Safe mmCIF backbone loader (NaN-safe, keeps incomplete residues)\n# ============================================================\n\nBACKBONE_ORDER = [\n    \"P\", \"OP1\", \"OP2\",\n    \"O5'\", \"C5'\", \"C4'\", \"C3'\", \"O3'\",\n    \"C2'\", \"C1'\", \"O4'\", \"O2'\"\n]\n\nBACKBONE_ATOMS = set(BACKBONE_ORDER)\n\n\ndef load_backbone_from_cif_safe(filepath):\n    \"\"\"\n    Safe mmCIF parser that:\n      - reads atom_site labels dynamically\n      - accepts ATOM + HETATM\n      - keeps only backbone atoms\n      - keeps altLoc '.', '?', 'A'\n      - preserves insertion codes\n      - KEEPS incomplete residues (fills missing atoms with NaN)\n      - skips chains with < 3 residues\n      - returns {chain_id: (N,12,3) array}\n    \"\"\"\n\n    import numpy as np\n\n    with open(filepath, \"r\") as f:\n        lines = f.readlines()\n\n    # -------------------------------\n    # Step 1: Find the atom_site loop\n    # -------------------------------\n    in_loop = False\n    atom_labels = []\n    data_start = None\n\n    for i, line in enumerate(lines):\n        line = line.strip()\n\n        if line == \"loop_\":\n            in_loop = True\n            atom_labels = []\n            continue\n\n        if in_loop:\n            if line.startswith(\"_atom_site.\"):\n                atom_labels.append(line)\n                continue\n            else:\n                if atom_labels:\n                    data_start = i\n                    break\n                else:\n                    in_loop = False\n\n    if not atom_labels:\n        return {}\n\n    # -----------------------------------------\n    # Step 2: Build label → column index mapping\n    # -----------------------------------------\n    label_to_index = {label: idx for idx, label in enumerate(atom_labels)}\n\n    required = [\n        \"_atom_site.label_atom_id\",\n        \"_atom_site.label_alt_id\",\n        \"_atom_site.label_comp_id\",\n        \"_atom_site.label_asym_id\",\n        \"_atom_site.label_seq_id\",\n        \"_atom_site.Cartn_x\",\n        \"_atom_site.Cartn_y\",\n        \"_atom_site.Cartn_z\",\n    ]\n\n    if any(r not in label_to_index for r in required):\n        return {}\n\n    # -----------------------------------------\n    # Step 3: Parse atom rows safely\n    # -----------------------------------------\n    residues = {}  # (chain, resid) → {atom_name: (x,y,z)}\n\n    for line in lines[data_start:]:\n        if not line or line.startswith(\"#\"):\n            continue\n\n        parts = line.split()\n        if len(parts) < len(atom_labels):\n            continue\n\n        atom_name = parts[label_to_index[\"_atom_site.label_atom_id\"]]\n        altloc    = parts[label_to_index[\"_atom_site.label_alt_id\"]]\n        chain_id  = parts[label_to_index[\"_atom_site.label_asym_id\"]]\n        resid     = parts[label_to_index[\"_atom_site.label_seq_id\"]]\n\n        if altloc not in (\".\", \"?\", \"A\"):\n            continue\n\n        if atom_name not in BACKBONE_ATOMS:\n            continue\n\n        try:\n            x = float(parts[label_to_index[\"_atom_site.Cartn_x\"]])\n            y = float(parts[label_to_index[\"_atom_site.Cartn_y\"]])\n            z = float(parts[label_to_index[\"_atom_site.Cartn_z\"]])\n        except ValueError:\n            continue\n\n        key = (chain_id, resid)\n        residues.setdefault(key, {})[atom_name] = (x, y, z)\n\n    # -----------------------------------------\n    # Step 4: Group residues into chains\n    # -----------------------------------------\n    chains = {}\n    for (chain_id, resid), atom_dict in residues.items():\n        chains.setdefault(chain_id, []).append((resid, atom_dict))\n\n    # -----------------------------------------\n    # Step 5: Sort residues within each chain\n    # -----------------------------------------\n    for chain_id in chains:\n        chains[chain_id].sort(key=lambda x: x[0])\n\n    # -----------------------------------------\n    # Step 6: Convert to (N,12,3) arrays (NaN-safe)\n    # -----------------------------------------\n    final = {}\n\n    for chain_id, reslist in chains.items():\n        coords = []\n\n        for resid, atom_dict in reslist:\n            row = []\n            for atom in BACKBONE_ORDER:\n                if atom in atom_dict:\n                    row.append(atom_dict[atom])\n                else:\n                    row.append((np.nan, np.nan, np.nan))  # fill missing atoms\n\n            coords.append(row)\n\n        if len(coords) >= 3:\n            final[chain_id] = np.array(coords, dtype=float)\n\n    return {k: final[k] for k in sorted(final.keys())}","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:22:17.525829Z","iopub.execute_input":"2026-02-21T09:22:17.526069Z","iopub.status.idle":"2026-02-21T09:22:17.542763Z","shell.execute_reply.started":"2026-02-21T09:22:17.52604Z","shell.execute_reply":"2026-02-21T09:22:17.542045Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compute_geom_stats(coords, max_coord=100.0):\n    \"\"\"\n    Compute biologically meaningful RNA geometry statistics:\n      - helix axis via PCA\n      - radius (distance to axis)\n      - rise (projection onto axis)\n      - twist (rotation around axis)\n      - noise (residual motion)\n    \"\"\"\n\n    coords = np.asarray(coords, float)\n\n    # -------------------------------\n    # 1. Filter invalid coordinates\n    # -------------------------------\n    valid = np.isfinite(coords).all(axis=1)\n    valid &= (np.abs(coords) < max_coord).all(axis=1)\n    coords = coords[valid]\n\n    if coords.shape[0] < 4:\n        return {\n            \"twist_mean\": 0.0, \"twist_std\": 0.01,\n            \"rise_mean\": 2.8,  \"rise_std\": 0.1,\n            \"radius_mean\": 11.0, \"radius_std\": 1.0,\n            \"noise_xy_std\": 0.1, \"noise_z_std\": 0.05,\n        }\n\n    # -------------------------------\n    # 2. Center coordinates\n    # -------------------------------\n    center = coords.mean(axis=0)\n    coords_c = coords - center\n\n    # -------------------------------\n    # 3. Estimate helix axis via PCA\n    # -------------------------------\n    U, S, Vt = np.linalg.svd(coords_c, full_matrices=False)\n    axis = Vt[0]  # principal direction\n    axis /= np.linalg.norm(axis)\n\n    # -------------------------------\n    # 4. Project coordinates onto axis\n    # -------------------------------\n    proj = coords_c @ axis  # scalar projection\n    diffs = np.diff(proj)\n\n    # Rise = projection step\n    rise = diffs\n\n    # -------------------------------\n    # 5. Compute twist\n    # -------------------------------\n    # Remove axis component to get radial vectors\n    radial = coords_c - np.outer(proj, axis)\n\n    # Compute angles in radial plane\n    angles = np.unwrap(np.arctan2(radial[:,1], radial[:,0]))\n    twist = np.diff(angles)\n\n    # -------------------------------\n    # 6. Compute radius\n    # -------------------------------\n    radius = np.linalg.norm(radial, axis=1)\n\n    # -------------------------------\n    # 7. Noise estimation\n    # -------------------------------\n    # Noise = deviation from ideal helix step\n    step_vecs = np.diff(coords_c, axis=0)\n    ideal_step = np.outer(rise, axis)\n    noise_vecs = step_vecs - ideal_step\n\n    noise_xy = np.linalg.norm(noise_vecs[:, :2], axis=1)\n    noise_z = np.abs(noise_vecs[:, 2])\n\n    # -------------------------------\n    # 8. Return stats\n    # -------------------------------\n    return {\n        \"twist_mean\": float(np.mean(twist)),\n        \"twist_std\":  float(np.std(twist) + 1e-8),\n\n        \"rise_mean\":  float(np.mean(rise)),\n        \"rise_std\":   float(np.std(rise) + 1e-8),\n\n        \"radius_mean\": float(np.mean(radius)),\n        \"radius_std\":  float(np.std(radius) + 1e-8),\n\n        \"noise_xy_std\": float(np.mean(noise_xy)),\n        \"noise_z_std\":  float(np.mean(noise_z)),\n    }","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:22:17.543599Z","iopub.execute_input":"2026-02-21T09:22:17.543829Z","iopub.status.idle":"2026-02-21T09:22:17.554751Z","shell.execute_reply.started":"2026-02-21T09:22:17.543803Z","shell.execute_reply":"2026-02-21T09:22:17.554076Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧬 Training-ready biological backbone generator (Updated to be safer)\n# ============================================================\n\ndef create_bio_backbone_structure(sequence, sequence_id, length, gc_content, sample_stats=None):\n    \"\"\"\n    Training-ready biological backbone generator.\n    Creates a simple right-handed helix with constant rise and twist.\n    Now includes:\n      - stat clamping\n      - twist unit fix\n      - overflow protection\n    \"\"\"\n\n    # Default stable geometry\n    tamed_stats = {\n        \"radius_mean\": 5.3,\n        \"radius_std\": 1.0,\n        \"twist_mean\": 0.25,   # radians per residue\n        \"twist_std\": 0.3,\n        \"rise_mean\": 2.8,\n        \"rise_std\": 0.5,\n        \"noise_xy_std\": 0.5,\n        \"noise_z_std\": 0.25,\n    }\n\n    # Use provided stats but clamp them to safe ranges\n    raw = sample_stats if sample_stats is not None else tamed_stats\n\n    radius = np.clip(raw.get(\"radius_mean\", 5.3), 1.0, 20.0)\n    rise   = np.clip(raw.get(\"rise_mean\",   2.8), 0.5, 10.0)\n\n    # FIX: twist_mean is already in radians, so DO NOT treat it as degrees\n    twist  = np.clip(raw.get(\"twist_mean\", 0.25), -1.0, 1.0)\n\n    coords = []\n    for i in range(length):\n        angle = i * twist  # radians directly\n        x = radius * np.cos(angle)\n        y = radius * np.sin(angle)\n        z = i * rise\n        coords.append((float(x), float(y), float(z)))\n\n    return [{\"coords\": coords, \"angles\": None}]","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2026-02-21T09:22:17.556993Z","iopub.execute_input":"2026-02-21T09:22:17.557408Z","iopub.status.idle":"2026-02-21T09:22:17.573112Z","shell.execute_reply.started":"2026-02-21T09:22:17.557386Z","shell.execute_reply":"2026-02-21T09:22:17.572569Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ------------------------------------------------------------\n# 🌉 Warp biological backbone into sample geometry (SAFE)\n# ------------------------------------------------------------\ndef warp_bio_to_sample(coords, sample_stats, max_scale=5.0):\n    \"\"\"\n    Warp biological backbone into sample geometry space.\n    Now includes:\n      - realistic stat clamping\n      - safe scaling\n      - drift prevention\n    \"\"\"\n\n    coords = np.asarray(coords, float)\n\n    # ⭐ Clamp sample stats to realistic RNA ranges ⭐\n    rs_mean = np.clip(sample_stats[\"radius_mean\"], 8.0, 20.0)\n    rs_std  = np.clip(sample_stats[\"radius_std\"],  0.1, 3.0)\n\n    tw_mean = np.clip(sample_stats[\"twist_mean\"], -0.5, 0.5)\n    tw_std  = np.clip(sample_stats[\"twist_std\"],   0.1, 1.0)\n\n    ri_mean = np.clip(sample_stats[\"rise_mean\"], 2.8, 3.4)\n    ri_std  = np.clip(sample_stats[\"rise_std\"],  0.05, 0.3)\n\n    xs, ys, zs = coords[:, 0], coords[:, 1], coords[:, 2]\n\n    # Compute biological stats safely\n    bio_stats = compute_geom_stats(coords)\n    if bio_stats is None:\n        return coords\n\n    # Avoid division by zero\n    b_rstd = max(bio_stats[\"radius_std\"], 1e-6)\n    b_tstd = max(bio_stats[\"twist_std\"],  1e-6)\n    b_rizd = max(bio_stats[\"rise_std\"],   1e-6)\n\n    # --- radius scaling ---\n    r = np.sqrt(xs**2 + ys**2)\n    r_centered = r - bio_stats[\"radius_mean\"]\n\n    scale_r = np.clip(rs_std / b_rstd, 1/max_scale, max_scale)\n    r_scaled = r_centered * scale_r + rs_mean\n\n    angle_xy = np.arctan2(ys, xs)\n    xs = r_scaled * np.cos(angle_xy)\n    ys = r_scaled * np.sin(angle_xy)\n\n    # --- twist warping ---\n    angles = np.unwrap(np.arctan2(ys, xs))\n    twist = np.diff(angles)\n\n    scale_t = np.clip(tw_std / b_tstd, 1/max_scale, max_scale)\n    twist_scaled = (twist - bio_stats[\"twist_mean\"]) * scale_t + tw_mean\n\n    angles_warped = angles.copy()\n    angles_warped[1:] = angles_warped[0] + np.cumsum(twist_scaled)\n\n    # --- rise warping ---\n    rise = np.diff(zs)\n\n    scale_ri = np.clip(ri_std / b_rizd, 1/max_scale, max_scale)\n    rise_scaled = (rise - bio_stats[\"rise_mean\"]) * scale_ri + ri_mean\n\n    zs_warped = zs.copy()\n    zs_warped[1:] = zs_warped[0] + np.cumsum(rise_scaled)\n\n    xs = r_scaled * np.cos(angles_warped)\n    ys = r_scaled * np.sin(angles_warped)\n\n    warped = np.stack([xs, ys, zs_warped], axis=1)\n\n    # ⭐ Center to prevent global drift ⭐\n    warped -= warped.mean(axis=0, keepdims=True)\n\n    return warped\n\n# ------------------------------------------------------------\n# 🌀 Hybrid generator (SAFE & REALISTIC)\n# ------------------------------------------------------------\ndef create_hybrid_structure(sequence, sequence_id, length, gc_content, sample_stats):\n    \"\"\"\n    Hybrid generator:\n    1) Build biological backbone\n    2) Warp into sample geometry\n    3) Add noise (clamped to realistic RNA levels)\n    4) Repeat 5 times\n    \"\"\"\n\n    structures = []\n\n    # 1) Biological backbone\n    bio_structures = create_bio_backbone_structure(\n        sequence, sequence_id, length, gc_content, sample_stats\n    )\n    bio_coords = np.array(bio_structures[0][\"coords\"], float)\n\n    # ⭐ Clamp noise to realistic RNA levels ⭐\n    nxy = np.clip(sample_stats[\"noise_xy_std\"], 0.0, 2.0)\n    nz  = np.clip(sample_stats[\"noise_z_std\"],  0.0, 0.5)\n\n    # 2) Warp + noise\n    for _ in range(5):\n        warped_coords = warp_bio_to_sample(bio_coords, sample_stats)\n\n        xs = warped_coords[:, 0] + np.random.normal(0, nxy, len(warped_coords))\n        ys = warped_coords[:, 1] + np.random.normal(0, nxy, len(warped_coords))\n        zs = warped_coords[:, 2] + np.random.normal(0, nz,  len(warped_coords))\n\n        coords = np.stack([xs, ys, zs], axis=1)\n\n        # ⭐ Final safety clamp ⭐\n        coords = np.nan_to_num(coords, nan=0.0, posinf=0.0, neginf=0.0)\n        coords = np.clip(coords, -50, 50)\n\n        structures.append({\"coords\": coords.tolist(), \"angles\": None})\n\n    print(\"DEBUG sample_stats:\", sample_stats)\n    print(\"DEBUG bio first residue:\", bio_coords[0])\n    print(\"DEBUG warped first residue:\", warped_coords[0])\n\n    return structures\n\n# ============================================================\n# 🧬 Biological-style vortex backbone + hybrid warping (SAFE & REALISTIC)\n# ============================================================\n\n# 1. Load CIF geometry\ncif_dir = \"/kaggle/input/stanford-rna-3d-folding-2/PDB_RNA\"\n\nsample_stats = compute_dataset_stats_streaming(cif_dir, compute_geom_stats)\n\nprint(\"Dataset geometry means:\")\nprint(\" Radius Mean =\", sample_stats[\"radius_mean\"])\nprint(\" Rise Mean   =\", sample_stats[\"rise_mean\"])\nprint(\" Twist Mean  =\", sample_stats[\"twist_mean\"])\nprint(\" Noise XY    =\", sample_stats[\"noise_xy_std\"])\nprint(\" Noise Z     =\", sample_stats[\"noise_z_std\"])\n\n# ============================================================\n# 🌀 RNA VORTEX PARAMETERS (full version)\n# ============================================================\n\nRNA_VORTEX_PARAMS = {\n    \"A\": {\n        \"radius\": None,\n        \"twist\": None,\n        \"rise\": None,\n        \"noise_xy\": None,\n        \"noise_z\": None,\n\n        # Base-specific offsets (you can tune these)\n        \"radius_offset\": 0.00,\n        \"twist_offset\":  0.00,\n        \"rise_offset\":   0.00,\n\n        # Optional: backbone phase shift\n        \"phase_shift\": 0.0,\n    },\n\n    \"C\": {\n        \"radius\": None,\n        \"twist\": None,\n        \"rise\": None,\n        \"noise_xy\": None,\n        \"noise_z\": None,\n\n        \"radius_offset\": 0.00,\n        \"twist_offset\":  0.00,\n        \"rise_offset\":   0.00,\n\n        \"phase_shift\": 0.0,\n    },\n\n    \"G\": {\n        \"radius\": None,\n        \"twist\": None,\n        \"rise\": None,\n        \"noise_xy\": None,\n        \"noise_z\": None,\n\n        \"radius_offset\": 0.00,\n        \"twist_offset\":  0.00,\n        \"rise_offset\":   0.00,\n\n        \"phase_shift\": 0.0,\n    },\n\n    \"U\": {\n        \"radius\": None,\n        \"twist\": None,\n        \"rise\": None,\n        \"noise_xy\": None,\n        \"noise_z\": None,\n\n        \"radius_offset\": 0.00,\n        \"twist_offset\":  0.00,\n        \"rise_offset\":   0.00,\n\n        \"phase_shift\": 0.0,\n    }\n}\n\nfor base in RNA_VORTEX_PARAMS:\n    RNA_VORTEX_PARAMS[base][\"radius\"] = sample_stats[\"radius_mean\"]\n    RNA_VORTEX_PARAMS[base][\"twist\"]  = sample_stats[\"twist_mean\"]\n    RNA_VORTEX_PARAMS[base][\"rise\"]   = sample_stats[\"rise_mean\"]\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-21T09:22:17.574047Z","iopub.execute_input":"2026-02-21T09:22:17.574257Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def chunk_structure(coords, max_len=2048, stride=2048):\n    chunks = []\n    for i in range(0, len(coords), stride):\n        chunk = coords[i:i+max_len]\n        if len(chunk) > 4:\n            chunks.append(chunk)\n    return chunks\n\nreal_train_structs = []\nfor target_id, coords in train_coords_dict.items():\n    coords = np.asarray(coords, dtype=np.float32)\n    if coords.ndim == 2 and coords.shape[1] == 3 and len(coords) > 4:\n        real_train_structs.extend(chunk_structure(coords))\n\nprint(f\"Real train structures: {len(real_train_structs)}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧮 Build hybrid_stats from REAL TRAINING GEOMETRY (SAFE)\n# ============================================================\n\nif len(real_train_structs) > 0:\n\n    print(\"📦 Concatenating real training structures...\")\n    all_train_coords = np.concatenate(real_train_structs, axis=0)\n\n    # ⭐ Center to remove global drift across structures\n    print(\"📐 Centering coordinates...\")\n    all_train_coords = all_train_coords - all_train_coords.mean(axis=0, keepdims=True)\n\n    # ⭐ Compute raw stats with progress indicator\n    print(\"📊 Computing geometry statistics...\")\n    hybrid_stats = compute_geom_stats(all_train_coords)\n\n    # ========================================================\n    # ⭐ Clamp to biologically realistic RNA ranges\n    # ========================================================\n    hybrid_stats[\"radius_mean\"] = float(np.clip(hybrid_stats[\"radius_mean\"], 5.0, 20.0))\n    hybrid_stats[\"radius_std\"]  = float(np.clip(hybrid_stats[\"radius_std\"],  0.5, 5.0))\n\n    hybrid_stats[\"rise_mean\"]   = float(np.clip(hybrid_stats[\"rise_mean\"],  2.5, 3.5))\n    hybrid_stats[\"rise_std\"]    = float(np.clip(hybrid_stats[\"rise_std\"],   0.1, 1.0))\n\n    hybrid_stats[\"twist_mean\"]  = float(np.clip(hybrid_stats[\"twist_mean\"], -0.5, 0.5))\n    hybrid_stats[\"twist_std\"]   = float(np.clip(hybrid_stats[\"twist_std\"],   0.1, 1.0))\n\n    hybrid_stats[\"noise_xy_std\"] = float(np.clip(0.5 * hybrid_stats[\"radius_std\"], 0.1, 2.0))\n    hybrid_stats[\"noise_z_std\"]  = float(np.clip(0.5 * hybrid_stats[\"rise_std\"],   0.05, 0.5))\n\n    # ========================================================\n    # Print final stats\n    # ========================================================\n    print(\"\\n📊 hybrid_stats (clamped + centered):\")\n    for k, v in hybrid_stats.items():\n        print(f\"  {k}: {v:.4f}\")\n\nelse:\n    raise ValueError(\"No real_train_structs available to compute hybrid_stats\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport numpy as np\nfrom Bio.PDB import MMCIFParser\n\ncif_folder = \"/kaggle/input/stanford-rna-3d-folding-2/PDB_RNA\"\nparser = MMCIFParser(QUIET=True)\n\nc1p_coords = []\nother_coords = []\n\nmax_files = 80          # tune: 40–100 is a good range\nother_sample_rate = 30  # take 1 out of every 30 non-C1' atoms\n\nfile_count = 0\n\nfor filename in os.listdir(cif_folder):\n    if not filename.endswith(\".cif\"):\n        continue\n\n    path = os.path.join(cif_folder, filename)\n    structure = parser.get_structure(filename, path)\n    file_count += 1\n\n    # Loop to gather samples into chuncks\n    for pdb_model in structure:\n        for chain in pdb_model:\n            for residue in chain:\n\n                # C1' atoms (directly relevant to competition)\n                if \"C1'\" in residue:\n                    c1p_coords.append(residue[\"C1'\"].get_coord())\n\n                # Optional: lightly sample other atoms for richer geometry\n                for atom_idx, atom in enumerate(residue):\n                    if atom.get_name() == \"C1'\":\n                        continue\n                    if atom_idx % other_sample_rate == 0:\n                        other_coords.append(atom.get_coord())\n\n    if file_count >= max_files:\n        break\n\n# Convert to arrays\nc1p_coords = np.array(c1p_coords)\nother_coords = np.array(other_coords) if other_coords else np.empty((0, 3))\n\n# Combine\nsample_coords = (\n    np.concatenate([c1p_coords, other_coords], axis=0)\n    if other_coords.size > 0 else c1p_coords\n)\n\nprint(\"Used CIF files:\", file_count)\nprint(\"c1p_coords shape:\", c1p_coords.shape)\nprint(\"other_coords shape:\", other_coords.shape)\nprint(\"sample_coords shape:\", sample_coords.shape)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🚀 Generate hybrid predictions and build submission\n# ============================================================\n\nprint(\"\\n🚀 Creating hybrid biological + sample-matched predictions...\")\nsubmission_rows = []\nrow_count = 0\n\nfor _, row in test_seqs.iterrows():\n    sequence = row.get(\"sequence\", \"\")\n    target_id = row.get(\"target_id\", \"UNKNOWN\")\n    all_sequences = row.get(\"all_sequences\", \"\")\n\n    if not isinstance(sequence, str) or len(sequence) == 0:\n        print(f\"⚠️ Skipping invalid sequence: {target_id}\")\n        continue\n\n    id_from_all = None\n    if isinstance(all_sequences, str):\n        parts = all_sequences.split(\">\")\n        if len(parts) > 1:\n            id_from_all = parts[1].split(\"|\")[0]\n\n    if not id_from_all:\n        id_from_all = f\"{target_id}_1\"\n\n    length = len(sequence)\n    gc_content = (sequence.count(\"G\") + sequence.count(\"C\")) / length\n\n    vortex_structures = create_hybrid_structure(\n        sequence,\n        id_from_all,\n        length,\n        gc_content,\n        hybrid_stats,\n    )\n\n    if row_count == 0:\n        print(\"\\n🌀 Example hybrid predictions for first residue:\")\n        for pred_num in range(1, 6):\n            x, y, z = vortex_structures[pred_num - 1][\"coords\"][0]\n            print(f\"  Pred{pred_num}: ({x:.3f}, {y:.3f}, {z:.3f})\")\n\n    for i, residue in enumerate(sequence, 1):\n        row_data = {\n            \"ID\": f\"{id_from_all}_{i}\",\n            \"resname\": residue,\n            \"resid\": int(i),\n        }\n\n        for pred_num in range(1, 6):\n            x, y, z = vortex_structures[pred_num - 1][\"coords\"][i - 1]\n            row_data[f\"x_{pred_num}\"] = float(x)\n            row_data[f\"y_{pred_num}\"] = float(y)\n            row_data[f\"z_{pred_num}\"] = float(z)\n\n        submission_rows.append(row_data)\n        row_count += 1\n\nprint(f\"\\n✅ Created {len(submission_rows)} hybrid predictions\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 📦 Build structure-level real / fake / hybrid sets (Into Chunks)\n# ============================================================\n\ndef chunk_structure(coords, max_len=2048, stride=2048):\n    chunks = []\n    for i in range(0, len(coords), stride):\n        chunk = coords[i:i+max_len]\n        if len(chunk) > 4:\n            chunks.append(chunk)\n    return chunks\n\n\n# -----------------------------\n# Real structures from training labels\n# -----------------------------\nreal_train_structs = []\nprint(\"📦 Building REAL TRAIN structures...\")\nfor target_id, coords in tqdm(train_coords_dict.items()):\n    coords = np.asarray(coords, dtype=np.float32)\n    if coords.ndim == 2 and coords.shape[1] == 3 and len(coords) > 4:\n        real_train_structs.extend(chunk_structure(coords))\n\n\n# -----------------------------\n# Real structures from validation labels\n# -----------------------------\nreal_val_structs = []\nprint(\"📦 Building REAL VAL structures...\")\nfor target_id, coords in tqdm(valid_coords_dict.items()):\n    coords = np.asarray(coords, dtype=np.float32)\n    if coords.ndim == 2 and coords.shape[1] >= 3 and len(coords) > 4:\n        c1 = coords[:, :3]\n        real_val_structs.extend(chunk_structure(c1))\n\nprint(f\"Real train structures: {len(real_train_structs)}\")\nprint(f\"Real val structures:   {len(real_val_structs)}\")\n\n\n# -----------------------------\n# Fake structures (chunked)\n# -----------------------------\ndef make_fake_structs(struct_list, noise_scale=5.0):\n    fake_list = []\n    for coords in struct_list:\n        noise = np.random.normal(0, noise_scale, coords.shape).astype(np.float32)\n        fake_list.append(coords + noise)\n    return fake_list\n\nprint(\"📦 Building FAKE structures...\")\nfake_train_structs = make_fake_structs(real_train_structs)\nfake_val_structs   = make_fake_structs(real_val_structs)\n\nprint(f\"Fake train structures: {len(fake_train_structs)}\")\nprint(f\"Fake val structures:   {len(fake_val_structs)}\")\n\n\n# -----------------------------\n# Hybrid structures (chunked)\n# -----------------------------\nval_seq_map = {row[\"target_id\"][:4]: row[\"sequence\"] for _, row in valid_seqs.iterrows()}\n\nhybrid_val_structs = []\nprint(\"📦 Building HYBRID VAL structures...\")\nfor target_id, coords in tqdm(valid_coords_dict.items()):\n\n    coords = np.asarray(coords, dtype=np.float32)\n    if coords.ndim == 2 and coords.shape[1] >= 3:\n        coords = coords[:, :3]\n\n    short_id = target_id if target_id in val_seq_map else target_id[:4]\n    seq = val_seq_map.get(short_id, None)\n    if seq is None:\n        continue\n\n    seq = str(seq)\n    length = len(seq)\n    if length < 5:\n        continue\n\n    gc_content = (seq.count(\"G\") + seq.count(\"C\")) / length\n\n    hybrid_structs = create_hybrid_structure(\n        seq, target_id, length, gc_content, hybrid_stats\n    )\n\n    coords_h = np.array(hybrid_structs[0][\"coords\"], dtype=np.float32)\n    hybrid_val_structs.extend(chunk_structure(coords_h))\n\nprint(f\"Hybrid val structures: {len(hybrid_val_structs)}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 📊 Compare REAL vs HYBRID geometry distributions\n# ============================================================\n\ndef geom_from_struct_list(struct_list):\n    all_coords = np.concatenate(struct_list, axis=0)\n    all_coords = all_coords - all_coords.mean(axis=0, keepdims=True)\n    return compute_geom_stats(all_coords)\n\nreal_stats = geom_from_struct_list(real_val_structs)\nhyb_stats  = geom_from_struct_list(hybrid_val_structs)\n\nprint(\"\\n📊 REAL vs HYBRID geometry comparison:\")\nfor key in real_stats:\n    r = real_stats[key]\n    h = hyb_stats[key]\n    print(f\"{key:12s}  real={r:8.3f}   hybrid={h:8.3f}   diff={h-r:8.3f}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧪 Structure-level dataset + transformer config\n# ============================================================\n\nimport torch\nfrom torch.utils.data import Dataset, DataLoader\n\n# -----------------------------\n# Transformer size config\n# -----------------------------\nTRANSFORMER_SIZE = \"medium\"  # \"medium\" or \"large\"\n\nif TRANSFORMER_SIZE == \"medium\":\n    D_MODEL = 256\n    N_HEADS = 8\n    N_LAYERS = 6\n    D_FF = 512\nelif TRANSFORMER_SIZE == \"large\":\n    D_MODEL = 512\n    N_HEADS = 12\n    N_LAYERS = 12\n    D_FF = 1024\nelse:\n    raise ValueError(\"TRANSFORMER_SIZE must be 'medium' or 'large'\")\n\nMAX_LEN = 2048  # chunk/pad length for transformer\n\n\n# -----------------------------\n# Structure-level dataset\n# -----------------------------\nclass StructureDataset(Dataset):\n    def __init__(self, real_structs, fake_structs, hybrid_structs=None):\n        \"\"\"\n        Each item is (coords[N,3], label)\n        label: 1 = real, 0 = fake, 2 = hybrid (if provided)\n        \"\"\"\n        self.samples = []\n        for s in real_structs:\n            self.samples.append((s, 1))\n        for s in fake_structs:\n            self.samples.append((s, 0))\n        if hybrid_structs is not None:\n            for s in hybrid_structs:\n                self.samples.append((s, 2))\n\n    def __len__(self):\n        return len(self.samples)\n\n    def __getitem__(self, idx):\n        coords, label = self.samples[idx]\n        coords = torch.tensor(coords, dtype=torch.float32)\n        label = torch.tensor(label, dtype=torch.long)\n        return coords, label\n\n\ndef collate_structures(batch):\n    \"\"\"\n    Pads/truncates each structure to MAX_LEN.\n    Returns:\n        coords_padded: [B, L, 3]\n        mask:          [B, L] (True for valid positions)\n        labels:        [B]\n    \"\"\"\n    coords_list, labels_list = zip(*batch)\n    batch_size = len(coords_list)\n\n    coords_padded = torch.zeros(batch_size, MAX_LEN, 3, dtype=torch.float32)\n    mask = torch.zeros(batch_size, MAX_LEN, dtype=torch.bool)\n    labels = torch.stack(labels_list, dim=0)\n\n    for i, coords in enumerate(coords_list):\n        n = min(coords.shape[0], MAX_LEN)\n        coords_padded[i, :n] = coords[:n]\n        mask[i, :n] = True\n\n    return coords_padded, mask, labels\n\n\n# -----------------------------\n# Build train/val datasets\n# -----------------------------\ntrain_dataset = StructureDataset(\n    real_structs=real_train_structs,\n    fake_structs=fake_train_structs,\n    hybrid_structs=None,  # use hybrid as extra class in train\n)\n\nval_dataset = StructureDataset(\n    real_structs=real_val_structs,\n    fake_structs=fake_val_structs,\n    hybrid_structs=None,  # keep val as real vs fake\n)\n\ntrain_loader = DataLoader(\n    train_dataset,\n    batch_size=8,\n    shuffle=True,\n    drop_last=True,\n    collate_fn=collate_structures,\n)\n\nval_loader = DataLoader(\n    val_dataset,\n    batch_size=8,\n    shuffle=False,\n    drop_last=False,\n    collate_fn=collate_structures,\n)\n\nprint(\"Train structures:\", len(train_dataset))\nprint(\"Val structures:  \", len(val_dataset))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧠 Structure-level transformer encoder (3-class) + NaN checks\n# ============================================================\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\ndef check_nan(tensor, name):\n    if torch.isnan(tensor).any():\n        raise ValueError(f\"❌ NaN detected inside model at: {name}\")\n    if torch.isinf(tensor).any():\n        raise ValueError(f\"❌ Inf detected inside model at: {name}\")\n\n\nclass GeometryEmbedding(nn.Module):\n    def __init__(self, d_model=256):\n        super().__init__()\n        self.mlp = nn.Sequential(\n            nn.Linear(3, d_model),\n            nn.ReLU(),\n            nn.Linear(d_model, d_model),\n        )\n\n    def forward(self, coords):\n        x = self.mlp(coords)\n        check_nan(x, \"GeometryEmbedding output\")\n        return x\n\n\nclass StructureLevelTransformerEncoder(nn.Module):\n    def __init__(\n        self,\n        d_model=256,\n        nhead=8,\n        num_layers=6,\n        dim_feedforward=512,\n        dropout=0.1,\n        num_classes=3,   # 0=fake, 1=real, 2=hybrid\n        proj_dim=64,\n    ):\n        super().__init__()\n\n        self.embed = GeometryEmbedding(d_model)\n\n        encoder_layer = nn.TransformerEncoderLayer(\n            d_model=d_model,\n            nhead=nhead,\n            dim_feedforward=dim_feedforward,\n            dropout=dropout,\n            batch_first=True,\n        )\n        self.encoder = nn.TransformerEncoder(\n            encoder_layer,\n            num_layers=num_layers,\n        )\n\n        self.classifier = nn.Linear(d_model, num_classes)\n        self.projection = nn.Linear(d_model, proj_dim)\n\n    def forward(self, coords, mask):\n        # coords: [B, L, 3]\n        # mask:   [B, L] (True = valid)\n        x = self.embed(coords)                 # [B, L, d_model]\n        check_nan(x, \"After embedding\")\n\n        src_key_padding_mask = ~mask           # True = pad\n        x = self.encoder(x, src_key_padding_mask=src_key_padding_mask)\n        check_nan(x, \"After transformer encoder\")\n\n        # mean pool over valid positions\n        mask_f = mask.unsqueeze(-1).float()\n        pooled = (x * mask_f).sum(dim=1) / (mask_f.sum(dim=1) + 1e-8)\n        check_nan(pooled, \"Pooled representation\")\n\n        logits = self.classifier(pooled)\n        check_nan(logits, \"Logits\")\n\n        z = self.projection(pooled)\n        check_nan(z, \"Contrastive projection\")\n\n        return logits, z\n\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\nmodel = StructureLevelTransformerEncoder(\n    d_model=D_MODEL,\n    nhead=N_HEADS,\n    num_layers=N_LAYERS,\n    dim_feedforward=D_FF,\n    dropout=0.1,\n    num_classes=3,\n    proj_dim=64,\n).to(device)\n\nce_loss = nn.CrossEntropyLoss()\noptimizer = torch.optim.Adam(model.parameters(), lr=1e-4)\n\nprint(\"Model initialized on\", device)\nprint(\"Total parameters:\", sum(p.numel() for p in model.parameters()))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧲 Supervised contrastive loss\n# ============================================================\n\ndef supervised_contrastive_loss(z, labels, temperature=0.1):\n    \"\"\"\n    z:      [B, D]  (projection vectors)\n    labels: [B]\n    \"\"\"\n    z = F.normalize(z, dim=1)\n    sim = torch.matmul(z, z.T) / temperature  # [B, B]\n\n    labels = labels.view(-1, 1)\n    mask = torch.eq(labels, labels.T).float()  # [B, B]\n\n    # remove self-comparisons\n    logits_mask = torch.ones_like(mask) - torch.eye(mask.size(0), device=mask.device)\n    mask = mask * logits_mask\n\n    # log-softmax over rows\n    log_prob = F.log_softmax(sim, dim=1)\n\n    # mean over positive pairs\n    numerator = (mask * log_prob).sum(dim=1)\n    denom = mask.sum(dim=1).clamp(min=1.0)\n    loss = -(numerator / denom).mean()\n    return loss\n\nlambda_contrastive = 0.1\nprint(\"Supervised contrastive loss ready.\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🛑 Early Stopping Utility\n# ============================================================\n\nclass EarlyStopping:\n    def __init__(self, patience=5, min_delta=1e-4):\n        self.patience = patience\n        self.min_delta = min_delta\n        self.best_loss = float(\"inf\")\n        self.counter = 0\n        self.should_stop = False\n\n    def step(self, val_loss):\n        if val_loss < self.best_loss - self.min_delta:\n            self.best_loss = val_loss\n            self.counter = 0\n        else:\n            self.counter += 1\n            if self.counter >= self.patience:\n                self.should_stop = True","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🧠 Structure-level transformer training loop WITH NaN DEBUGGING\n# ============================================================\n\ndef dbg(msg):\n    print(msg, flush=True)\n\ndbg(\"🚀 Starting training loop with NaN detection...\")\nearly_stop = EarlyStopping(patience=5)\nEPOCHS = 3\n\nfor epoch in range(EPOCHS):\n    dbg(f\"\\n======================\\n🔥 EPOCH {epoch+1} START\\n======================\")\n\n    model.train()\n    total_loss = 0.0\n    total_ce = 0.0\n    total_con = 0.0\n\n    dbg(\"Entering training batches...\")\n\n    for batch_idx, (coords, mask, labels) in enumerate(train_loader):\n\n        coords = coords.to(device)   # [B, L, 3]\n        mask   = mask.to(device)     # [B, L]\n        labels = labels.to(device)   # [B]\n\n        # ⭐ SANITIZE INPUTS BEFORE MODEL FORWARD ⭐\n        coords = torch.nan_to_num(coords, nan=0.0, posinf=0.0, neginf=0.0)\n        coords = torch.clamp(coords, -1000, 1000)\n\n        # Forward\n        logits, z = model(coords, mask)\n\n        # NaN / Inf checks on logits\n        if torch.isnan(logits).any():\n            dbg(f\"❌ NaN detected in logits at batch {batch_idx}\")\n            dbg(f\"logits sample: {logits[0][:5]}\")\n            raise ValueError(\"NaN in logits\")\n\n        if torch.isinf(logits).any():\n            dbg(f\"❌ Inf detected in logits at batch {batch_idx}\")\n            raise ValueError(\"Inf in logits\")\n\n        # Losses\n        loss_ce = ce_loss(logits, labels)\n        loss_con = supervised_contrastive_loss(z, labels)\n        loss = loss_ce + lambda_contrastive * loss_con\n\n        # NaN checks on losses\n        if torch.isnan(loss_ce):\n            dbg(f\"❌ NaN in CE loss at batch {batch_idx}\")\n            dbg(f\"logits sample: {logits[0][:5]}\")\n            dbg(f\"labels: {labels}\")\n            raise ValueError(\"NaN in CE loss\")\n\n        if torch.isnan(loss_con):\n            dbg(f\"❌ NaN in contrastive loss at batch {batch_idx}\")\n            dbg(f\"z sample: {z[0][:5]}\")\n            dbg(f\"labels: {labels}\")\n            raise ValueError(\"NaN in contrastive loss\")\n\n        if torch.isnan(loss):\n            dbg(f\"❌ NaN in TOTAL loss at batch {batch_idx}\")\n            raise ValueError(\"NaN in total loss\")\n\n\n        # Backprop\n        optimizer.zero_grad()\n        loss.backward()\n\n        # Gradient NaN check\n        for name, param in model.named_parameters():\n            if param.grad is not None and torch.isnan(param.grad).any():\n                dbg(f\"❌ NaN in gradients of {name}\")\n                raise ValueError(\"NaN in gradients\")\n\n        optimizer.step()\n\n        total_loss += loss.item()\n        total_ce += loss_ce.item()\n        total_con += loss_con.item()\n\n        # ⭐ PRINT ONLY FINAL BATCH OF THE EPOCH ⭐\n        if batch_idx == len(train_loader) - 1:\n            dbg(\n                f\"Final Batch {batch_idx} | \"\n                f\"Loss={loss.item():.4f} | \"\n                f\"CE={loss_ce.item():.4f} | \"\n                f\"Con={loss_con.item():.4f}\"\n            )\n\n    # ---- End of training epoch ----\n    n_train = max(len(train_loader), 1)\n    train_log = (\n        f\"[Train] Epoch {epoch+1} | \"\n        f\"Total: {total_loss/n_train:.4f} | \"\n        f\"CE: {total_ce/n_train:.4f} | \"\n        f\"Con: {total_con/n_train:.4f}\"\n    )\n\n    # -----------------------------\n    # Validation\n    # -----------------------------\n    dbg(\"Starting validation...\")\n    model.eval()\n    val_ce = 0.0\n    correct = 0\n    total = 0\n\n    with torch.no_grad():\n        for v_idx, (coords, mask, labels) in enumerate(val_loader):\n\n            coords = coords.to(device)\n            mask   = mask.to(device)\n            labels = labels.to(device)\n\n            # ⭐ SANITIZE VALIDATION INPUTS TOO ⭐\n            coords = torch.nan_to_num(coords, nan=0.0, posinf=0.0, neginf=0.0)\n            coords = torch.clamp(coords, -1000, 1000)\n\n            logits, z = model(coords, mask)\n            loss_ce = ce_loss(logits, labels)\n            val_ce += loss_ce.item()\n\n            preds = logits.argmax(dim=1)\n            correct += (preds == labels).sum().item()\n            total += labels.size(0)\n\n    n_val = max(len(val_loader), 1)\n    val_loss = val_ce / n_val\n    val_log = (\n        f\"[Val] CE: {val_loss:.4f} | \"\n        f\"Acc: {correct/max(total,1):.4f}\"\n    )\n\n    dbg(train_log + \"  ||  \" + val_log)\n\n    # -----------------------------\n    # Early stopping check\n    # -----------------------------\n    early_stop.step(val_loss)\n\n    if early_stop.should_stop:\n        print(\"⛔ Early stopping triggered — training halted.\")\n        break\n\ndbg(\"✅ Training complete.\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 🔧 Refine hybrid_stats using validation geometry\n# ============================================================\n\n# Concatenate all validation real coords for stats\nif len(real_val_structs) > 0:\n    all_val_coords = np.concatenate(real_val_structs, axis=0)\n    val_geom_stats = compute_geom_stats(all_val_coords)\n    refined_hybrid_stats = hybrid_stats.copy()\n\n    # Blend means\n    for key in [\"radius_mean\", \"twist_mean\", \"rise_mean\"]:\n        refined_hybrid_stats[key] = 0.7 * hybrid_stats[key] + 0.3 * val_geom_stats[key]\n\n    # Blend std/noise gently\n    for key in [\"radius_std\", \"twist_std\", \"rise_std\", \"noise_xy_std\", \"noise_z_std\"]:\n        if key in val_geom_stats:\n            refined_hybrid_stats[key] = 0.7 * hybrid_stats[key] + 0.3 * val_geom_stats[key]\nelse:\n    refined_hybrid_stats = hybrid_stats.copy()\n\nprint(\"Refined hybrid_stats:\")\nfor k, v in refined_hybrid_stats.items():\n    print(f\"  {k}: {v:.4f}\")\n\n\n# ============================================================\n# 📈 Structure-level scoring utilities\n# ============================================================\n\n# Build real cluster center from real train + val structures\nwith torch.no_grad():\n    all_real_structs = real_train_structs + real_val_structs\n    embed_list = []\n\n    for coords in real_val_structs:\n        coords = np.asarray(coords, dtype=np.float32)\n    \n        # Convert to tensor\n        coords_t = torch.tensor(coords, dtype=torch.float32).unsqueeze(0).to(device)\n        n = min(coords_t.shape[1], MAX_LEN)\n    \n        # Pad\n        coords_pad = torch.zeros(1, MAX_LEN, 3, device=device)\n        mask_pad = torch.zeros(1, MAX_LEN, dtype=torch.bool, device=device)\n        coords_pad[0, :n] = coords_t[0, :n]\n        mask_pad[0, :n] = True\n    \n        # ⭐ SANITIZE BEFORE MODEL ⭐\n        coords_pad = torch.nan_to_num(coords_pad, nan=0.0, posinf=0.0, neginf=0.0)\n        coords_pad = torch.clamp(coords_pad, -1000, 1000)\n    \n        # Forward\n        _, z = model(coords_pad, mask_pad)\n        embed_list.append(z[0].cpu().numpy())\n    \n    # Compute center\n    real_center = torch.tensor(np.mean(embed_list, axis=0), dtype=torch.float32).to(device)\n\n\ndef score_structure_coords(coords):\n    \"\"\"\n    coords: np.ndarray [N,3]\n    Returns a single realness score for the entire structure.\n    \"\"\"\n    coords = np.asarray(coords, dtype=np.float32)\n    coords_t = torch.tensor(coords, dtype=torch.float32).unsqueeze(0).to(device)\n    n = min(coords_t.shape[1], MAX_LEN)\n    coords_pad = torch.zeros(1, MAX_LEN, 3, device=device)\n    mask_pad = torch.zeros(1, MAX_LEN, dtype=torch.bool, device=device)\n    coords_pad[0, :n] = coords_t[0, :n]\n    mask_pad[0, :n] = True\n\n    logits, z = model(coords_pad, mask_pad)\n    probs = torch.softmax(logits, dim=1)[0]\n    prob_real = probs[1].item()  # class 1 = real\n\n    z_mean = z[0]\n    dist = torch.norm(z_mean - real_center).item()\n    sim_score = 1.0 / (1.0 + dist)\n\n    final = 0.6 * prob_real + 0.4 * sim_score\n    return final\n\n\n# ============================================================\n# 🧬 Regenerate submission with per-target best-of-5 selection\n# ============================================================\n\nprint(\"\\n🚀 Regenerating hybrid predictions with refined stats and best-of-5 scoring...\")\nsubmission_rows = []\nrow_count = 0\n\nfor _, row in test_seqs.iterrows():\n    sequence = row.get(\"sequence\", \"\")\n    target_id = row.get(\"target_id\", \"UNKNOWN\")\n    all_sequences = row.get(\"all_sequences\", \"\")\n\n    if not isinstance(sequence, str) or len(sequence) == 0:\n        print(f\"⚠️ Skipping invalid sequence: {target_id}\")\n        continue\n\n    id_from_all = None\n    if isinstance(all_sequences, str):\n        parts = all_sequences.split(\">\")\n        if len(parts) > 1:\n            id_from_all = parts[1].split(\"|\")[0]\n    if not id_from_all:\n        id_from_all = f\"{target_id}_1\"\n\n    length = len(sequence)\n    gc_content = (sequence.count(\"G\") + sequence.count(\"C\")) / length\n\n    # Generate 5 candidate structures\n    vortex_structures = create_hybrid_structure(\n        sequence,\n        id_from_all,\n        length,\n        gc_content,\n        refined_hybrid_stats,\n    )\n\n    # Score all 5 and pick the best\n    scores = []\n    for s in vortex_structures:\n        coords = np.array(s[\"coords\"], dtype=np.float32)\n        scores.append(score_structure_coords(coords))\n    best_idx = int(np.argmax(scores))\n    best_struct = vortex_structures[best_idx]\n\n    if row_count == 0:\n        print(\"\\n🌀 Example best-of-5 scores for first target:\")\n        for i, s in enumerate(scores):\n            print(f\"  Struct {i}: score={s:.4f}\")\n        x, y, z = best_struct[\"coords\"][0]\n        print(f\"  Best struct first residue: ({x:.3f}, {y:.3f}, {z:.3f})\")\n\n    # Add best structure to submission (same coords for all 5 preds)\n    for i, residue in enumerate(sequence, 1):\n        row_data = {\n            \"ID\": f\"{id_from_all}_{i}\",\n            \"resname\": residue,\n            \"resid\": int(i),\n        }\n        x, y, z = best_struct[\"coords\"][i - 1]\n        for pred_num in range(1, 6):\n            row_data[f\"x_{pred_num}\"] = float(x)\n            row_data[f\"y_{pred_num}\"] = float(y)\n            row_data[f\"z_{pred_num}\"] = float(z)\n\n        submission_rows.append(row_data)\n        row_count += 1\n\nprint(f\"\\n✅ Created {len(submission_rows)} best-of-5 hybrid predictions\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 📊 Validation-based evaluation of hybrid generator (NaN-safe + centered + aligned)\n# ============================================================\n\ndef kabsch_align(P, Q):\n    \"\"\"\n    Align P to Q using the Kabsch algorithm.\n    P, Q: [N, 3]\n    Returns aligned P.\n    \"\"\"\n    P = np.asarray(P, dtype=np.float32)\n    Q = np.asarray(Q, dtype=np.float32)\n\n    # Center\n    Pc = P - P.mean(axis=0, keepdims=True)\n    Qc = Q - Q.mean(axis=0, keepdims=True)\n\n    C = Pc.T @ Qc\n    V, S, Wt = np.linalg.svd(C)\n    d = np.sign(np.linalg.det(V @ Wt))\n    D = np.diag([1, 1, d])\n    R = V @ D @ Wt\n    P_aligned = Pc @ R\n    return P_aligned\n\ndef rmsd(a, b):\n    a = np.asarray(a, dtype=np.float32)\n    b = np.asarray(b, dtype=np.float32)\n\n    a = np.nan_to_num(a, nan=0.0, posinf=0.0, neginf=0.0)\n    b = np.nan_to_num(b, nan=0.0, posinf=0.0, neginf=0.0)\n\n    a = np.clip(a, -1000, 1000)\n    b = np.clip(b, -1000, 1000)\n\n    # Center both\n    a = a - a.mean(axis=0, keepdims=True)\n    b = b - b.mean(axis=0, keepdims=True)\n\n    # Optional: align a to b (shape-only RMSD)\n    a = kabsch_align(a, b)\n\n    return np.sqrt(np.mean((a - b) ** 2))\n\nval_rmsds = []\nval_twist_err = []\nval_rise_err = []\nval_radius_err = []\n\nfor target_id, coords in valid_coords_dict.items():\n    coords = np.asarray(coords, dtype=np.float32)\n    if coords.ndim != 2 or coords.shape[1] < 3:\n        continue\n\n    true_c1 = coords[:, :3]\n\n    # Find matching sequence\n    mask = valid_seqs[\"target_id\"].str.startswith(target_id)\n    if not mask.any():\n        continue\n    seq_row = valid_seqs[mask].iloc[0]\n    sequence = str(seq_row[\"sequence\"])\n    length = len(sequence)\n    if length < 5:\n        continue\n\n    n = min(len(true_c1), length)\n    true_c1 = true_c1[:n]\n    sequence = sequence[:n]\n\n    gc_content = (sequence.count(\"G\") + sequence.count(\"C\")) / len(sequence)\n\n    # Generate one hybrid structure with refined stats\n    hybrid_structs = create_hybrid_structure(\n        sequence,\n        target_id,\n        len(sequence),\n        gc_content,\n        refined_hybrid_stats,\n    )\n    hybrid_coords = np.array(hybrid_structs[0][\"coords\"], dtype=np.float32)[:n]\n\n    # Sanitize\n    true_c1 = np.nan_to_num(true_c1, nan=0.0, posinf=0.0, neginf=0.0)\n    hybrid_coords = np.nan_to_num(hybrid_coords, nan=0.0, posinf=0.0, neginf=0.0)\n\n    true_c1 = np.clip(true_c1, -1000, 1000)\n    hybrid_coords = np.clip(hybrid_coords, -1000, 1000)\n\n    # RMSD (centered + aligned)\n    val_rmsds.append(rmsd(true_c1, hybrid_coords))\n\n    # Geometry stats (use raw sanitized coords, no alignment)\n    true_geom = compute_geom_stats(true_c1)\n    hybrid_geom = compute_geom_stats(hybrid_coords)\n\n    val_twist_err.append(abs(true_geom[\"twist_mean\"] - hybrid_geom[\"twist_mean\"]))\n    val_rise_err.append(abs(true_geom[\"rise_mean\"] - hybrid_geom[\"rise_mean\"]))\n    val_radius_err.append(abs(true_geom[\"radius_mean\"] - hybrid_geom[\"radius_mean\"]))\n\nif len(val_rmsds) > 0:\n    print(\"\\n📊 Validation geometry evaluation:\")\n    print(f\"  Mean RMSD:        {np.mean(val_rmsds):.4f}\")\n    print(f\"  Mean |twist_err|: {np.mean(val_twist_err)c:.4f}\")\n    print(f\"  Mean |rise_err|:  {np.mean(val_rise_err):.4f}\")\n    print(f\"  Mean |radius_err|:{np.mean(val_radius_err):.4f}\")\nelse:\n    print(\"⚠️ No validation structures could be evaluated.\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# 💾 Save submission file\n# ============================================================\n\nsubmission_df = pd.DataFrame(submission_rows)\nsubmission_df = submission_df[\n    [\"ID\", \"resname\", \"resid\"]\n    + [f\"{axis}_{i}\" for i in range(1, 6) for axis in [\"x\", \"y\", \"z\"]]\n]\n\nsubmission_df.to_csv(\"submission.csv\", index=False)\nprint(\"📁 Saved submission.csv\")","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}