{"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":"none","dataSources":[{"sourceType":"competition","sourceId":118765,"databundleVersionId":15231210},{"sourceType":"datasetVersion","sourceId":15129160,"datasetId":9688076,"databundleVersionId":16017634},{"sourceType":"datasetVersion","sourceId":10855324,"datasetId":6742586,"databundleVersionId":11219268},{"sourceType":"kernelVersion","sourceId":290004465}],"dockerImageVersionId":31286,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport re\nfrom collections import deque\nfrom sklearn.ensemble import RandomForestRegressor\nfrom numba import njit \nimport warnings\nwarnings.filterwarnings('ignore')\n\n# --- PATHS ---\nTRAIN_PATH = \"/kaggle/input/datasets/pjleek/sovereign-gold-star/sovereign_gold_star.csv\"\nFALLBACK_PATH = \"/kaggle/input/sovereign-gold-star/sovereign_gold_star.csv\"\nTEST_PATH = \"/kaggle/input/competitions/stanford-rna-3d-folding-2/test_sequences.csv\"\nOUTPUT_PATH = \"submission.csv\"\n\nnp.random.seed(42)\nVOCAB = {'A': 0, 'C': 1, 'G': 2, 'U': 3, 'N': 4}\nNUM_CHARS = len(VOCAB)\n\n# ========================================================\n# 1. METADATA & ROBUST STOICHIOMETRY PARSING\n# ========================================================\ndef parse_stoich_and_seq(seq, stoich):\n    parts = str(stoich).split(';')\n    total_copies = sum([int(p.split(':')[1]) if ':' in p else 1 for p in parts])\n    unique_chains = set([p.split(':')[0] for p in parts if ':' in p])\n            \n    if total_copies > 1:\n        chunk = len(seq) // total_copies\n        first = seq[:chunk]\n        if len(unique_chains) == 1 and len(parts) == 1:\n            return total_copies, chunk, np.zeros(chunk, dtype=bool), first\n            \n        breaks = np.zeros(len(seq), dtype=bool)\n        for i in range(1, total_copies): \n            if i * chunk < len(seq):\n                breaks[i * chunk] = True\n        return 1, len(seq), breaks, seq\n            \n    return 1, len(seq), np.zeros(len(seq), dtype=bool), seq\n\n# ========================================================\n# 2. UNIVERSAL ASSEMBLY\n# ========================================================\ndef get_universal_assembly(monomer_coords, total_copies, avg_hoog_pot=0.0):\n    if total_copies <= 1: \n        return monomer_coords - np.mean(monomer_coords, axis=0)\n\n    if avg_hoog_pot <= 0.6:\n        monomer_coords = monomer_coords - np.mean(monomer_coords, axis=0) \n\n    rg = np.sqrt(np.mean(np.sum(monomer_coords**2, axis=1)))\n    rg_modifier = np.clip(1.2 - avg_hoog_pot, 0.5, 1.0)\n    offset = max(1.5, (total_copies * rg * rg_modifier) / 6.28) \n    \n    assembled = []\n    for i in range(total_copies):\n        theta = (2 * np.pi * i) / total_copies\n        cos_t, sin_t = np.cos(theta), np.sin(theta)\n        rot = np.array([[cos_t, -sin_t, 0], [sin_t, cos_t, 0], [0, 0, 1]])\n        \n        c_i_base = np.copy(monomer_coords)\n        \n        if avg_hoog_pot > 0.6:\n            z_shift = 3.4 * (i % 2)\n            if i % 2 == 1:\n                c_i_base = c_i_base * np.array([1.0, -1.0, -1.0])\n            c_i = (c_i_base @ rot.T) + np.array([offset * cos_t, offset * sin_t, z_shift])\n        else:\n            c_i = (c_i_base @ rot.T) + np.array([offset * cos_t, offset * sin_t, 0])\n            \n        assembled.append(c_i)\n        \n    complex_coords = np.vstack(assembled)\n    complex_coords -= np.mean(complex_coords, axis=0)\n    return complex_coords\n\n# ========================================================\n# 3. FAST DP & FEATURE EXTRACTION \n# ========================================================\n@njit(fastmath=True)\ndef _numba_dp(seq_int):\n    N = len(seq_int)\n    dp, match = np.zeros((N, N), dtype=np.int32), np.zeros((N, N), dtype=np.bool_)\n    for i in range(N):\n        for j in range(i+4, N):\n            p1, p2 = seq_int[i], seq_int[j]\n            if (p1==0 and p2==3) or (p1==3 and p2==0) or (p1==1 and p2==2) or \\\n               (p1==2 and p2==1) or (p1==2 and p2==3) or (p1==3 and p2==2): match[i, j] = True\n    for k in range(4, N):\n        for i in range(N - k):\n            j = i + k\n            opt1, opt2 = dp[i, j-1], 0\n            for t in range(i, j-3):\n                if match[t, j]:\n                    val1 = dp[i, t-1] if t-1>=i else 0\n                    val2 = dp[t+1, j-1] if j-1>=t+1 else 0\n                    if val1 + val2 + 1 > opt2: opt2 = val1 + val2 + 1\n            dp[i, j] = max(opt1, opt2)\n            \n    pairs, stack_i, stack_j = np.full(N, -1, dtype=np.int32), [0], [N-1]\n    while len(stack_i) > 0:\n        i, j = stack_i.pop(), stack_j.pop()\n        if i >= j - 4: continue\n        if dp[i, j] == dp[i, j-1]: stack_i.append(i); stack_j.append(j-1)\n        else:\n            for t in range(i, j-3):\n                if match[t, j]:\n                    val1 = dp[i, t-1] if t-1>=i else 0\n                    val2 = dp[t+1, j-1] if j-1>=t+1 else 0\n                    if val1 + val2 + 1 == dp[i, j]:\n                        pairs[t] = j; pairs[j] = t\n                        stack_i.extend([i, t+1]); stack_j.extend([t-1, j-1])\n                        break\n    return pairs\n\ndef compute_pairing_features(seq):\n    seq_map = {'A':0, 'C':1, 'G':2, 'U':3}\n    seq_int = np.array([seq_map.get(c, 4) for c in seq], dtype=np.int32)\n    pairs = _numba_dp(seq_int)\n    is_paired = np.array([1.0 if p != -1 else 0.0 for p in pairs])\n    is_loop = 1.0 - is_paired\n    is_g_tract = np.zeros(len(seq))\n    for m in re.finditer(r'G{2,}', seq): is_g_tract[m.start():m.end()] = 1.0\n    return seq_int, is_paired, is_loop, is_g_tract, pairs\n\ndef detect_3d_motifs(seq):\n    L = len(seq)\n    motif_flags = np.zeros((L, 2)) \n    for match in re.finditer(r'(?=G[ACGU][AG]A)', seq): motif_flags[match.start():match.start()+4, 0] = 1.0\n    for match in re.finditer(r'(?=GA|AG)', seq): motif_flags[match.start():match.start()+2, 1] = 1.0\n    return motif_flags\n\ndef get_encoded(seq, win):\n    mapped = np.array([VOCAB.get(c, 4) for c in seq])\n    pad_len = win // 2\n    padded = np.pad(mapped, (pad_len, pad_len), constant_values=4)\n    strides = (padded.strides[0], padded.strides[0])\n    windows = np.lib.stride_tricks.as_strided(padded, shape=(len(seq), win), strides=strides)\n    return np.eye(NUM_CHARS)[windows].reshape(len(seq), -1)\n\ndef compute_chemical_features(seq):\n    L = len(seq)\n    if L == 0: return np.zeros((0, 7))\n        \n    n_weights = {'A': 0.6, 'G': 1.0, 'C': 0.2, 'U': 0.2}\n    n_exp = np.array([n_weights.get(c, 0.0) for c in seq])\n    w_n = max(1, min(101, L if L % 2 == 1 else L - 1))\n    n_score = np.convolve(n_exp, np.ones(w_n), mode='same')\n    if np.max(n_score) > 0: n_score /= np.max(n_score)\n        \n    hoogsteen = np.array([1.0 if c in['A', 'G'] else 0.0 for c in seq])\n    is_motif = np.array([1.0 if seq[i:i+2] in ['GG', 'GA'] else 0.0 for i in range(L-1)] + [0.0])\n    w_m = max(1, min(21, L if L % 2 == 1 else L - 1))\n    motif_dens = np.convolve(is_motif, np.ones(w_m), mode='same')\n    if np.max(motif_dens) > 0: motif_dens /= np.max(motif_dens)\n        \n    idx = np.arange(L)\n    weights = np.exp(-(np.abs(idx[:, None] - idx[None, :]) * 5.45)**2 / 16.0)\n    nitrogen_cluster_score = weights @ n_exp \n    rg_expected = 3.0 * (L ** 0.35)\n    \n    return np.column_stack([\n        np.full(L, L), (np.arange(L) + 1) / L, nitrogen_cluster_score, \n        np.full(L, rg_expected), n_score, hoogsteen, motif_dens\n    ])\n\ndef compute_v62_lw_features(seq, breaks=None):\n    L = len(seq)\n    if L == 0: return np.zeros((0, 12))\n    \n    max_chain = L\n    if breaks is not None and np.any(breaks):\n        b_idx = np.where(breaks)[0]\n        max_chain = np.max(np.diff(np.concatenate(([0], b_idx, [L]))))\n        \n    win_size = 3 if max_chain < 50 else 7\n    half_win = win_size // 2\n    \n    lw_logic = {\n        'A': [1.0, 1.0, 1.0, 0.8, 1.0, 1.0], \n        'G': [1.0, 1.0, 1.0, 1.0, 2.0, 1.0], \n        'C': [0.0, 1.0, 0.2, 0.5, 1.0, 2.0], \n        'U': [0.0, 1.0, 0.2, 0.5, 1.0, 2.0], \n        'N': [0.0, 0.5, 0.2, 0.2, 1.0, 1.0]  \n    }\n    \n    base_lw = np.array([lw_logic.get(c, lw_logic['N']) for c in seq], dtype=np.float32)\n    v62_feats = np.zeros((L, 12), dtype=np.float32)\n    \n    for i in range(L):\n        start = max(0, i - half_win)\n        end = min(L, i + half_win + 1)\n        window = base_lw[start:end]\n        \n        v62_feats[i, 0:6] = np.mean(window, axis=0)\n        v62_feats[i, 6] = i / max(1, L - 1)\n        v62_feats[i, 7] = 1.0 if i < 6 else 0.0\n        v62_feats[i, 8] = 1.0 if (L - i - 1) < 6 else 0.0\n        v62_feats[i, 9] = np.sum(window[:, 0]) / max(1, end - start)\n        v62_feats[i, 10] = np.mean(np.var(window, axis=0)) \n        v62_feats[i, 11] = 1.0 if max_chain < 50 else 0.0 \n        \n    return v62_feats\n\n# ========================================================\n# 4. DATA-DRIVEN TRAINING\n# ========================================================\nprint(\"[*] Loading Sovereign Gold Star dataset...\")\ntry:\n    df_robust = pd.read_csv(TRAIN_PATH)\nexcept FileNotFoundError:\n    try:\n        df_robust = pd.read_csv(FALLBACK_PATH)\n    except FileNotFoundError:\n        df_robust = pd.DataFrame({'pdb_id':['x'], 'chain':['A'], 'res_num':[1], 'res_name':['A'], 'kappa':[0], 'tau':[0], 't_x':[0], 't_y':[0], 't_z':[0], 'pcv_x':[0], 'pcv_y':[0], 'pcv_z':[0], 'topo_density':[0], 'topo_avg_dist':[0]})\n\ndf_test_raw = pd.read_csv(TEST_PATH)\n\nmetadata_rows = []\nfor _, row in df_test_raw.iterrows():\n    target_id, seq_raw = row['target_id'], row['sequence']\n    stoich = str(row.get('stoichiometry', 'A:1'))\n    ligands_raw = str(row.get('ligand_ids', '')).upper()\n    has_ligand = pd.notna(row.get('ligand_ids')) and ligands_raw not in ['NAN', 'NONE', '']\n    assembly_copies, L_mono, breaks, seq_mono = parse_stoich_and_seq(seq_raw, stoich)\n    \n    v62_feat = compute_v62_lw_features(seq_mono, breaks)\n    avg_hoog_pot = np.mean(v62_feat[:, 2]) if len(v62_feat) > 0 else 0.0\n    \n    metadata_rows.append({\n        'target_id': target_id, 'seq_raw': seq_raw, 'seq_mono': seq_mono, \n        'L_mono': L_mono, 'assembly_copies': assembly_copies, 'breaks': breaks,\n        'has_ligand': has_ligand, 'ligands_raw': ligands_raw, 'avg_hoog_pot': avg_hoog_pot\n    })\ndf_metadata = pd.DataFrame(metadata_rows)\n\nprint(\"[*] Extracting Structural Environment Features...\")\nX_stack, Y_stack, W_stack, X_anchor, Y_anchor = [], [], [], [], []\n\nfor _, group in df_robust.groupby(['pdb_id', 'chain']):\n    group = group.sort_values('res_num').reset_index(drop=True)\n    seq = ''.join(group['res_name'].astype(str).str[0].values)\n    \n    if len(seq) > 600: continue \n        \n    steps = np.full(len(group), 5.45)\n    if 'c1_x' in group.columns:\n        coords = group[['c1_x', 'c1_y', 'c1_z']].values\n        raw_steps = np.linalg.norm(coords[1:] - coords[:-1], axis=1)\n        steps[1:] = np.clip(raw_steps, 0.0, 8.5) \n\n    motifs = detect_3d_motifs(seq)\n    seq_int, is_paired, is_loop, is_g_tract, _ = compute_pairing_features(seq)\n    chem_feat = compute_chemical_features(seq)\n    v62_feat = compute_v62_lw_features(seq)\n    \n    w = np.ones(len(seq))\n    w[is_loop == 1.0] = 3.0\n    w[is_g_tract == 1.0] = 4.0\n    if 'bond_count' in group.columns:\n        bond_counts = group['bond_count'].fillna(0).values.astype(float)\n        w *= (1.0 + np.log1p(bond_counts)) \n    \n    features = np.hstack([get_encoded(seq, 5), motifs, is_paired[:, None], is_loop[:, None], is_g_tract[:, None], v62_feat])\n    X_stack.append(features)\n    Y_stack.append(group[['kappa', 'tau', 't_x', 't_y', 't_z']].assign(step_size=steps).values.astype(float))\n    W_stack.extend(w)\n    \n    X_anchor.append(np.hstack([get_encoded(seq, 7), chem_feat, motifs, is_loop[:, None], is_g_tract[:, None], v62_feat]))\n    Y_anchor.append(group[['pcv_x', 'pcv_y', 'pcv_z', 'topo_density', 'topo_avg_dist']].values.astype(float))\n\nX_stack, Y_stack = np.vstack(X_stack), np.vstack(Y_stack)\nX_anchor, Y_anchor = np.vstack(X_anchor), np.vstack(Y_anchor)\nW_stack = np.array(W_stack)\n\nprint(\"[*] Training Dual Engine Pathfinders...\")\nm_stack = RandomForestRegressor(n_estimators=100, max_depth=20, random_state=42, n_jobs=-1).fit(X_stack, Y_stack, sample_weight=W_stack)\nm_anchor = RandomForestRegressor(n_estimators=80, max_depth=15, random_state=42, n_jobs=-1).fit(X_anchor, Y_anchor, sample_weight=W_stack)\n\n# ========================================================\n# 5. THE GLOBAL PATHFINDER ENGINE (The Selective Anarchist)\n# ========================================================\ndef rodrigues_rot(v, k, theta):\n    cos_t, sin_t = np.cos(theta), np.sin(theta)\n    return v * cos_t + np.cross(k, v) * sin_t + k * np.dot(k, v) * (1 - cos_t)\n\ndef apply_clutch(base_arr, chain_breaks, lw_variance, win):\n    \"\"\"Dynamic multi-tier smoothing array generator.\"\"\"\n    if win <= 1: return np.copy(base_arr)\n    smoothed = np.copy(base_arr)\n    N = len(base_arr)\n    for i in range(N):\n        if lw_variance[i] > 0.05: continue\n        start = max(0, i - win//2)\n        end = min(N, i + win//2 + 1)\n        for k in range(i, start-1, -1):\n            if chain_breaks[k]: start = k; break\n        for k in range(i+1, end):\n            if chain_breaks[k]: end = k; break\n        if end > start:\n            smoothed[i] = np.mean(base_arr[start:end])\n    return smoothed\n\ndef predict_ensemble(seq, L, chain_breaks, m_stack, m_anchor, row_meta):\n    if L <= 1: return np.zeros((L, 5, 3))\n    \n    seq_int, is_paired, is_loop, is_g_tract, pairs = compute_pairing_features(seq)\n    motifs = detect_3d_motifs(seq)\n    chem_feat = compute_chemical_features(seq)\n    v62_feat = compute_v62_lw_features(seq, chain_breaks)\n    \n    max_chain = L\n    if np.any(chain_breaks):\n        b_idx = np.where(chain_breaks)[0]\n        max_chain = np.max(np.diff(np.concatenate(([0], b_idx, [L]))))\n    \n    lw_variance = v62_feat[:, 10]\n    \n    features_stack = np.hstack([get_encoded(seq, 5), motifs, is_paired[:, None], is_loop[:, None], is_g_tract[:, None], v62_feat])\n    features_anchor = np.hstack([get_encoded(seq, 7), chem_feat, motifs, is_loop[:, None], is_g_tract[:, None], v62_feat])\n    \n    p_stack = m_stack.predict(features_stack)\n    p_anchor = m_anchor.predict(features_anchor)\n    base_pcvs, base_dens, base_dist = p_anchor[:, 0:3], p_anchor[:, 3], p_anchor[:, 4]\n    \n    # Target Coordinates\n    target_coords = np.zeros((L, 3))\n    for i in range(L):\n        v = base_pcvs[i]\n        n = np.linalg.norm(v)\n        if n > 1e-8: target_coords[i] = (v / n) * base_dist[i]\n        else: target_coords[i] = np.array([0., 0., 0.])\n            \n    target_coords_smooth = np.copy(target_coords)\n    win_t = max(3, min(9, L // 20))\n    for i in range(L):\n        start = max(0, i - win_t//2)\n        end = min(L, i + win_t//2 + 1)\n        for k in range(i, start-1, -1):\n            if chain_breaks[k]: start = k; break\n        for k in range(i+1, end):\n            if chain_breaks[k]: end = k; break\n        target_coords_smooth[i] = np.mean(target_coords[start:end], axis=0)\n\n    target_coords_smooth -= target_coords_smooth[0]\n\n    # Base Physics Values\n    step_base, tau_base, kappa_base = np.zeros(L), np.zeros(L), np.zeros(L)\n    for i in range(L):\n        tau_rf, kappa_rf = p_stack[i, 1], p_stack[i, 0]\n        if is_g_tract[i]:\n            step_base[i] = np.clip(p_stack[i, 5], 4.0, 6.8)\n            tau_base[i], kappa_base[i] = tau_rf, kappa_rf\n        elif is_paired[i]:\n            step_base[i] = 5.45\n            tau_base[i] = (tau_rf * 0.5) + (0.55 * 0.5)\n            kappa_base[i] = (kappa_rf * 0.5) + (0.30 * 0.5)\n        else:\n            step_base[i] = np.clip(p_stack[i, 5], 3.5, 8.0)\n            tau_base[i], kappa_base[i] = tau_rf, kappa_rf\n            \n    # 🌟 V69: The 3-Tier Smoothing Engine\n    win_mid = 5 if max_chain >= 50 else 3\n    win_max = 9 if max_chain >= 50 else 5\n    \n    step_mid = apply_clutch(step_base, chain_breaks, lw_variance, win_mid)\n    tau_mid = apply_clutch(tau_base, chain_breaks, lw_variance, win_mid)\n    kappa_mid = apply_clutch(kappa_base, chain_breaks, lw_variance, win_mid)\n    \n    step_max = apply_clutch(step_base, chain_breaks, lw_variance, win_max)\n    tau_max = apply_clutch(tau_base, chain_breaks, lw_variance, win_max)\n    kappa_max = apply_clutch(kappa_base, chain_breaks, lw_variance, win_max)\n    \n    step_raw, tau_raw, kappa_raw = np.copy(step_base), np.copy(tau_base), np.copy(kappa_base)\n    \n    mean_dens = np.mean(base_dens) if L > 0 else 0.0\n    \n    ensemble_pos = np.zeros((5, 3))\n    ensemble_fwds = np.array([[1.0, 0.0, 0.0]] * 5)\n    ensemble_ups = np.array([[0.0, 1.0, 0.0]] * 5)\n    \n    history_pulls = [deque(maxlen=3) for _ in range(5)]\n    all_coords = np.zeros((L, 5, 3))\n    all_ups = np.zeros((L, 5, 3)) \n    \n    for i in range(L):\n        all_coords[i, :, :] = ensemble_pos\n        if i == L - 1: break\n                \n        if i > 0 and chain_breaks[i]:\n            for j in range(5):\n                history_pulls[j].clear()\n                \n                # Fetch profile specific global weight for chain break logic\n                w_g, _, _ = (0.5, 0.0, 0.0)\n                if j == 1: w_g = 0.8\n                elif j == 2: w_g = 0.5 if L > 50 else 0.7\n                elif j == 3: w_g = 0.6 if mean_dens > 4.5 else (0.3 if mean_dens < 3.0 else 0.5)\n                elif j == 4: w_g = 0.2\n                \n                # The Anarchist Chain Scatter (only if long enough to warrant it)\n                if j == 4 and L > 50:\n                    core_center = np.mean(all_coords[:i, j, :], axis=0) if i > 0 else np.zeros(3)\n                    ensemble_pos[j] = core_center + (np.random.randn(3) * 15.0)\n                elif w_g >= 0.5:\n                    ensemble_pos[j] = target_coords_smooth[i]\n                else:\n                    core_center = np.mean(all_coords[:i, j, :], axis=0) if i > 0 else np.zeros(3)\n                    ensemble_pos[j] = core_center + np.array([12.0 * np.cos(i), 12.0 * np.sin(i), 4.0])\n        \n        tx, ty, tz = p_stack[i][2:5]\n        partner_idx = pairs[i]\n        t_vec = np.array([tx, ty, tz])\n        t_vec_norm = (t_vec / np.linalg.norm(t_vec)) if np.linalg.norm(t_vec) > 1e-8 else None\n        \n        for j in range(5):\n            curr_kappa, curr_tau, curr_step = kappa_mid[i], tau_mid[i], step_mid[i]\n            w_glob, k_mult, t_mult, w_repel, c_lim = 0.5, 1.0, 1.0, 0.0, 0.0\n            \n            # 🌟 Profile 1: The Conservative (v67 Pure Anchor)\n            if j == 0:\n                pass \n                \n            # 🌟 Profile 2: The Refiner (High Target Pull, Mild Repulsion)\n            elif j == 1:\n                w_glob = 0.8\n                w_repel = 0.1\n                c_lim = 2.5\n                \n            # 🌟 Profile 3: The Shifter (Phase/Chiral modifier for L > 50)\n            elif j == 2:\n                if L > 50:\n                    curr_kappa = kappa_mid[min(i+1, L-1)]\n                    curr_tau = tau_mid[max(0, i-1)]\n                    if not is_paired[i] and not is_g_tract[i]:\n                        curr_tau = -curr_tau # Mirror chirality in loops\n                else:\n                    w_glob = 0.7 # Safe fallback for short sequences\n                    \n            # 🌟 Profile 4: The Specialist (Density-Adaptive Smoothing)\n            elif j == 3:\n                if mean_dens > 4.5:\n                    curr_kappa, curr_tau, curr_step = kappa_max[i], tau_max[i], step_max[i]\n                    w_glob = 0.6\n                elif mean_dens < 3.0:\n                    curr_kappa, curr_tau, curr_step = kappa_raw[i], tau_raw[i], step_raw[i]\n                    w_glob = 0.3\n                    w_repel, c_lim = 0.1, 2.5\n                    \n            # 🌟 Profile 5: The Anarchist (High Jitter, Loop Teleporting if L > 50)\n            elif j == 4:\n                curr_kappa, curr_tau, curr_step = kappa_raw[i], tau_raw[i], step_raw[i]\n                k_mult = 1.0 + (np.random.rand() - 0.5) * 0.2\n                t_mult = 1.0 + (np.random.rand() - 0.5) * 0.2\n                w_glob, w_repel, c_lim = 0.2, 0.3, 3.5\n                \n                if L > 50 and i > 0 and not is_paired[i] and is_loop[i] and (i % 8 == 0):\n                    ensemble_pos[j] += (np.random.rand(3) - 0.5) * 3.5\n            \n            fwd_j, up_j = ensemble_fwds[j], ensemble_ups[j]\n            t_vec_local = t_vec_norm if t_vec_norm is not None else fwd_j\n            right = np.cross(up_j, fwd_j)\n            right = (right / np.linalg.norm(right)) if np.linalg.norm(right) > 1e-8 else np.array([0., 0., 1.])\n            \n            f_stack_vec = rodrigues_rot(fwd_j, right, curr_kappa * k_mult)\n            u_next = rodrigues_rot(up_j, f_stack_vec, curr_tau * t_mult)\n            if np.dot(up_j, u_next) < 0: u_next = -u_next \n                \n            if is_paired[i] and not is_g_tract[i] and partner_idx < i:\n                partner_up = all_ups[partner_idx, j, :]\n                u_next = (u_next * 0.6) - (partner_up * 0.4)\n                u_next /= (np.linalg.norm(u_next) + 1e-8)\n                \n            repulsion_pull = np.zeros(3)\n            if w_repel > 0.0 and i > 4:\n                prev_coords = all_coords[:i-2, j, :]\n                vecs_to_prev = ensemble_pos[j] - prev_coords\n                dists_to_prev = np.linalg.norm(vecs_to_prev, axis=1)\n                clashes = dists_to_prev < c_lim \n                if np.any(clashes):\n                    repulse_vecs = vecs_to_prev[clashes] / (dists_to_prev[clashes, None]**3 + 1e-8)\n                    repulsion_pull = np.sum(repulse_vecs, axis=0)\n                    r_norm = np.linalg.norm(repulsion_pull)\n                    if r_norm > 1.0: repulsion_pull /= r_norm\n\n            f_local = ((0.5 * f_stack_vec) + (0.5 * t_vec_local) + (w_repel * repulsion_pull))\n            f_local /= (np.linalg.norm(f_local) + 1e-8)\n            \n            if w_glob < 0.5:\n                if is_paired[i] and partner_idx < i:\n                    v_part = all_coords[partner_idx, j, :] - ensemble_pos[j]\n                    d_part = np.linalg.norm(v_part)\n                    if d_part > 11.5: f_local += (v_part / d_part) * 0.5\n                    elif d_part < 9.5: f_local -= (v_part / d_part) * 0.5\n                f_local /= (np.linalg.norm(f_local) + 1e-8)\n\n            f_global = target_coords_smooth[i+1] - ensemble_pos[j]\n            n_glob = np.linalg.norm(f_global)\n            f_global = (f_global / n_glob) if n_glob > 1e-8 else fwd_j\n\n            fwd_final = (f_local * (1.0 - w_glob)) + (f_global * w_glob)\n            fwd_final /= (np.linalg.norm(fwd_final) + 1e-8)\n            \n            next_pos = ensemble_pos[j] + (fwd_final * curr_step)\n\n            ensemble_pos[j] = next_pos\n            ensemble_fwds[j] = fwd_final \n            ensemble_ups[j] = u_next - np.dot(u_next, fwd_final) * fwd_final\n            ensemble_ups[j] = (ensemble_ups[j] / np.linalg.norm(ensemble_ups[j])) if np.linalg.norm(ensemble_ups[j]) > 1e-8 else np.cross(right, fwd_final)\n            \n            all_ups[i, j, :] = ensemble_ups[j]\n            \n    return all_coords\n\n# ========================================================\n# 6. ASSEMBLY & OUTPUT\n# ========================================================\nprint(\"[*] Generating Selective Anarchist Ensembles...\")\nfinal_coords = []\n\nfor _, row in df_metadata.iterrows():\n    target_id = row['target_id']\n    seq_raw = row['seq_raw']\n    \n    p_opt_ensemble = predict_ensemble(\n        row['seq_mono'], row['L_mono'], row['breaks'], m_stack, m_anchor, row\n    )\n    \n    p_final_full = []\n    for j in range(5):\n        p_final_full.append(get_universal_assembly(p_opt_ensemble[:, j, :], row['assembly_copies'], row['avg_hoog_pot']))\n    \n    p_final = np.stack(p_final_full, axis=1) \n\n    for i in range(len(seq_raw)):\n        idx = min(i, p_final.shape[0] - 1)\n        row_id = str(target_id) + \"_\" + str(i + 1)\n        row_dict = {'ID': row_id, 'target_id': target_id, 'resname': seq_raw[i], 'resid': i + 1}\n        \n        for j in range(5):\n            row_dict['x_' + str(j+1)] = p_final[idx, j, 0]\n            row_dict['y_' + str(j+1)] = p_final[idx, j, 1]\n            row_dict['z_' + str(j+1)] = p_final[idx, j, 2]\n            \n        final_coords.append(row_dict)\n\nexpected_columns = [\n    'ID', 'target_id', 'resname', 'resid', \n    'x_1', 'y_1', 'z_1', 'x_2', 'y_2', 'z_2', \n    'x_3', 'y_3', 'z_3', 'x_4', 'y_4', 'z_4', 'x_5', 'y_5', 'z_5'\n]\n\npd.DataFrame(final_coords)[expected_columns].to_csv(OUTPUT_PATH, index=False)\nprint(\"\\n[+] V110 (The Selective Anarchist) saved to \" + str(OUTPUT_PATH))\nprint(\"    -> Profile 1: The Conservative (Restored v67 Base)\")\nprint(\"    -> Profile 2: The Refiner (Smooth + PCV Pull)\")\nprint(\"    -> Profile 3: The Shifter (Phase/Chiral Mod, gated to L > 50)\")\nprint(\"    -> Profile 4: The Specialist (Adaptive Topo Density Smoothing)\")\nprint(\"    -> Profile 5: The Anarchist (Jitter + Loop Teleport, gated to L > 50)\")\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-12T21:25:55.070467Z","iopub.execute_input":"2026-03-12T21:25:55.071493Z","iopub.status.idle":"2026-03-12T21:26:08.414881Z","shell.execute_reply.started":"2026-03-12T21:25:55.071459Z","shell.execute_reply":"2026-03-12T21:26:08.413935Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport subprocess\n\n# 1. Define the correct path to the binary folder\nusalign_dir = \"/kaggle/input/datasets/metric/usalign\"\nusalign_bin = os.path.join(usalign_dir, \"USalign\")\n\n# 2. Add the directory to the system PATH so metric.py finds it\nos.environ['PATH'] += f\":{usalign_dir}\"\n\n# 3. Ensure the binary is executable (Kaggle inputs can be finicky with permissions)\n# We copy it to /tmp to ensure we have 'chmod' control\nif not os.path.exists(\"/tmp/USalign\"):\n    !cp {usalign_bin} /tmp/USalign\n    !chmod +x /tmp/USalign\n    # Add /tmp to PATH as well to be safe\n    os.environ['PATH'] += \":/tmp\"\n\n# 4. Verify the system can now \"see\" USalign\nresult = subprocess.run([\"which\", \"USalign\"], capture_output=True, text=True)\nprint(f\"USalign found at: {result.stdout.strip()}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-12T21:26:08.416302Z","iopub.execute_input":"2026-03-12T21:26:08.41672Z","iopub.status.idle":"2026-03-12T21:26:08.436275Z","shell.execute_reply.started":"2026-03-12T21:26:08.416692Z","shell.execute_reply":"2026-03-12T21:26:08.434981Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport runpy\n\n# 1. Load the metric\nmodule_globals = runpy.run_path(\"/kaggle/usr/lib/notebooks/rhijudas/tm-score-permutechains/metric.py\")\nscore = module_globals['score']\n\n# Read in the validation solution and the output submission.csv from your notebook\nimport pandas as pd\nsol = pd.read_csv('/kaggle/input/competitions/stanford-rna-3d-folding-2/validation_labels.csv')\nsub = pd.read_csv('/kaggle/working/submission.csv')\n\n# Run the eval code on each target and get the score, if this is a notebook run using the public test set (whose size matches the solution file, `validation_labels.csv`]\nsol['target_id'] = sol['ID'].apply(lambda x: '_'.join(str(x).split('_')[:-1]))\nsub['target_id'] = sub['ID'].apply(lambda x: '_'.join(str(x).split('_')[:-1]))\n\nif len(sol)==len(sub): # This tests if we're looking at public val\n    results = []\n    for target_id, group_native in sol.groupby('target_id'):\n        group_predicted = sub[sub['target_id'] == target_id]\n        result = score(group_native,group_predicted,'ID')\n        print(target_id,result)\n        results.append( result )\n    print( 'Mean score:',  \n          float(sum(results) / len(results)) if len(results)>0 else 0.0, \n          f'(n={len(results)})' )\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-12T21:26:08.437909Z","iopub.execute_input":"2026-03-12T21:26:08.438281Z","iopub.status.idle":"2026-03-12T21:35:58.348579Z","shell.execute_reply.started":"2026-03-12T21:26:08.438244Z","shell.execute_reply":"2026-03-12T21:35:58.34765Z"}},"outputs":[],"execution_count":null}]}