{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.0"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":118765,"databundleVersionId":15231210,"sourceType":"competition"}],"dockerImageVersionId":31234,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!touch submission.csv\nprint('submission.csv created')\n!pip install ViennaRNA -q 2>/dev/null || echo \"ViennaRNA not available, using simple secondary structure prediction\"\n!sleep 5","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom scipy.spatial.transform import Rotation as R\nfrom scipy.optimize import minimize\nimport random\nfrom Bio import pairwise2\nfrom Bio.Seq import Seq\nimport time\nfrom collections import defaultdict\nfrom sklearn.preprocessing import normalize\nfrom scipy.spatial import distance_matrix\nfrom scipy.spatial.distance import cdist\n\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nimport seaborn as sns\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Try to import ViennaRNA for secondary structure prediction\ntry:\n    import RNA\n    HAS_VIENNA = True\n    print(\"ViennaRNA available for secondary structure prediction\")\nexcept ImportError:\n    HAS_VIENNA = False\n    print(\"ViennaRNA not available, using Nussinov algorithm\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_seqs = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/train_sequences.csv')\nvalid_seqs = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/validation_sequences.csv')\ntest_seqs = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/test_sequences.csv')\ntrain_labels = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/train_labels.csv')\nvalid_labels = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/validation_labels.csv')\n\nprint(f\"Loaded {len(train_seqs)} training sequences, {len(valid_seqs)} validation sequences, and {len(test_seqs)} test sequences\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Process training labels to create a dictionary mapping target_id to its 3D coordinates\ndef process_labels(labels_df):\n    \"\"\"\n    Process labels dataframe to create a dictionary mapping target_id to coordinates\n    \"\"\"\n    coords_dict = {}\n    \n    # Group by target ID\n    for id_prefix, group in labels_df.groupby(lambda x: labels_df['ID'][x].rsplit('_', 1)[0]):\n        # Extract just the coordinates columns for the first structure (x_1, y_1, z_1)\n        coords = []\n        for _, row in group.sort_values('resid').iterrows():\n            coords.append([row['x_1'], row['y_1'], row['z_1']])\n        \n        coords_dict[id_prefix] = np.array(coords)\n    \n    return coords_dict\n\n# Process training labels\nprint(\"Processing training labels...\")\ntrain_coords_dict = process_labels(train_labels)\nvalid_coords_dict = process_labels(valid_labels)\nprint(f\"Processed coordinates for {len(train_coords_dict)} training structures and {len(valid_coords_dict)} validation structures\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class KmerIndex:\n    \"\"\"K-mer based index for fast sequence similarity search.\"\"\"\n    \n    def __init__(self, k=6):\n        self.k = k\n        self.index = defaultdict(set)\n        self.sequences = {}\n        \n    def add_sequence(self, target_id, sequence):\n        \"\"\"Add a sequence to the index.\"\"\"\n        self.sequences[target_id] = sequence\n        kmers = self._get_kmers(sequence)\n        for kmer in kmers:\n            self.index[kmer].add(target_id)\n    \n    def _get_kmers(self, sequence):\n        \"\"\"Extract all k-mers from a sequence.\"\"\"\n        return set(sequence[i:i+self.k] for i in range(len(sequence) - self.k + 1))\n    \n    def find_candidates(self, query_seq, min_shared_kmers=3, max_candidates=50):\n        \"\"\"Find candidate sequences with shared k-mers.\"\"\"\n        query_kmers = self._get_kmers(query_seq)\n        \n        # Count shared k-mers for each sequence\n        candidate_scores = defaultdict(int)\n        for kmer in query_kmers:\n            for target_id in self.index.get(kmer, []):\n                candidate_scores[target_id] += 1\n        \n        # Filter and sort by shared k-mer count\n        candidates = [(tid, score) for tid, score in candidate_scores.items() \n                     if score >= min_shared_kmers]\n        candidates.sort(key=lambda x: x[1], reverse=True)\n        \n        return [tid for tid, _ in candidates[:max_candidates]]\n\n# Build k-mer index\nprint(\"Building k-mer index...\")\nkmer_index = KmerIndex(k=6)\nfor _, row in train_seqs.iterrows():\n    if row['target_id'] in train_coords_dict:\n        kmer_index.add_sequence(row['target_id'], row['sequence'])\nprint(f\"Indexed {len(kmer_index.sequences)} sequences\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def nussinov_fold(sequence, min_loop_size=3):\n    \"\"\"Simple Nussinov algorithm for RNA secondary structure prediction.\"\"\"\n    n = len(sequence)\n    \n    # Base pair scores\n    def can_pair(b1, b2):\n        pairs = {('A', 'U'), ('U', 'A'), ('G', 'C'), ('C', 'G'), ('G', 'U'), ('U', 'G')}\n        return (b1, b2) in pairs\n    \n    # DP table\n    dp = np.zeros((n, n), dtype=int)\n    \n    # Fill DP table\n    for length in range(min_loop_size + 2, n + 1):\n        for i in range(n - length + 1):\n            j = i + length - 1\n            \n            # Case 1: i unpaired\n            dp[i, j] = dp[i + 1, j] if i + 1 <= j else 0\n            \n            # Case 2: j unpaired\n            if i <= j - 1:\n                dp[i, j] = max(dp[i, j], dp[i, j - 1])\n            \n            # Case 3: i pairs with j\n            if can_pair(sequence[i], sequence[j]):\n                inner = dp[i + 1, j - 1] if i + 1 <= j - 1 else 0\n                dp[i, j] = max(dp[i, j], inner + 1)\n            \n            # Case 4: bifurcation\n            for k in range(i + 1, j):\n                dp[i, j] = max(dp[i, j], dp[i, k] + dp[k + 1, j])\n    \n    # Traceback to get structure\n    structure = ['.'] * n\n    \n    def traceback(i, j):\n        if i >= j:\n            return\n        \n        if dp[i, j] == dp[i + 1, j]:\n            traceback(i + 1, j)\n        elif dp[i, j] == dp[i, j - 1]:\n            traceback(i, j - 1)\n        elif can_pair(sequence[i], sequence[j]):\n            inner = dp[i + 1, j - 1] if i + 1 <= j - 1 else 0\n            if dp[i, j] == inner + 1:\n                structure[i] = '('\n                structure[j] = ')'\n                traceback(i + 1, j - 1)\n                return\n        \n        # Bifurcation\n        for k in range(i + 1, j):\n            if dp[i, j] == dp[i, k] + dp[k + 1, j]:\n                traceback(i, k)\n                traceback(k + 1, j)\n                return\n    \n    traceback(0, n - 1)\n    return ''.join(structure)\n\n\ndef clean_sequence(sequence):\n    \"\"\"Clean RNA sequence, replacing invalid characters with valid ones.\"\"\"\n    valid_bases = {'A', 'C', 'G', 'U'}\n    cleaned = []\n    for base in sequence.upper():\n        if base == 'T':\n            cleaned.append('U')  # Convert DNA T to RNA U\n        elif base in valid_bases:\n            cleaned.append(base)\n        else:\n            # Replace unknown bases with A (most common)\n            cleaned.append('A')\n    return ''.join(cleaned)\n\n\ndef predict_secondary_structure(sequence):\n    \"\"\"Predict RNA secondary structure using ViennaRNA or Nussinov.\"\"\"\n    # Clean sequence first\n    clean_seq = clean_sequence(sequence)\n    \n    if HAS_VIENNA:\n        try:\n            fc = RNA.fold_compound(clean_seq)\n            structure, mfe = fc.mfe()\n            return structure\n        except Exception as e:\n            # Fall back to Nussinov if ViennaRNA fails\n            return nussinov_fold(clean_seq)\n    else:\n        return nussinov_fold(clean_seq)\n\n\ndef get_base_pairs(structure):\n    \"\"\"Extract base pairs from dot-bracket notation.\"\"\"\n    pairs = []\n    stack = []\n    \n    for i, char in enumerate(structure):\n        if char == '(':\n            stack.append(i)\n        elif char == ')':\n            if stack:\n                j = stack.pop()\n                pairs.append((j, i))\n    \n    return pairs\n\n\ndef structure_similarity(struct1, struct2):\n    \"\"\"Calculate similarity between two secondary structures.\"\"\"\n    pairs1 = set(get_base_pairs(struct1))\n    pairs2 = set(get_base_pairs(struct2))\n    \n    if not pairs1 or not pairs2:\n        return 0.0\n    \n    intersection = len(pairs1 & pairs2)\n    union = len(pairs1 | pairs2)\n    \n    return intersection / union if union > 0 else 0.0\n\n\n# Cache secondary structures for training sequences\nprint(\"Computing secondary structures for training sequences...\")\ntrain_structures = {}\nfor idx, row in train_seqs.iterrows():\n    if row['target_id'] in train_coords_dict:\n        train_structures[row['target_id']] = predict_secondary_structure(row['sequence'])\n    if (idx + 1) % 500 == 0:\n        print(f\"  Processed {idx + 1}/{len(train_seqs)} sequences...\")\nprint(f\"Computed {len(train_structures)} secondary structures\")","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def find_similar_sequences_fast(query_seq, train_seqs_df, train_coords_dict, \n                                kmer_index, train_structures, temporal_cutoff=None, top_n=5):\n    \"\"\"\n    Find similar sequences using k-mer pre-filtering and combined sequence+structure scoring.\n    \"\"\"\n    # Get candidates using k-mer index\n    candidates = kmer_index.find_candidates(query_seq, min_shared_kmers=2, max_candidates=100)\n    \n    if not candidates:\n        # Fall back to all sequences if no k-mer matches\n        candidates = list(train_coords_dict.keys())\n    \n    # Predict secondary structure for query\n    query_structure = predict_secondary_structure(query_seq)\n    \n    # Filter by temporal cutoff and score candidates\n    similar_seqs = []\n    query_seq_obj = Seq(query_seq)\n    \n    # Create lookup for temporal cutoff filtering\n    if temporal_cutoff:\n        valid_targets = set(train_seqs_df[train_seqs_df['temporal_cutoff'] < temporal_cutoff]['target_id'])\n        candidates = [c for c in candidates if c in valid_targets]\n    \n    for target_id in candidates:\n        if target_id not in train_coords_dict:\n            continue\n            \n        train_seq = kmer_index.sequences.get(target_id)\n        if not train_seq:\n            continue\n        \n        # Skip if sequence is too different in length\n        len_ratio = abs(len(train_seq) - len(query_seq)) / max(len(train_seq), len(query_seq))\n        if len_ratio > 0.5:\n            continue\n        \n        # Sequence alignment score\n        alignments = pairwise2.align.globalms(query_seq_obj, train_seq, 2, -1, -10, -0.5, one_alignment_only=True)\n        \n        if not alignments:\n            continue\n            \n        alignment = alignments[0]\n        seq_score = alignment.score / (2 * min(len(query_seq), len(train_seq)))\n        \n        # Secondary structure similarity score\n        struct_score = 0.0\n        if target_id in train_structures:\n            struct_score = structure_similarity(query_structure, train_structures[target_id])\n        \n        # Combined score (weighted)\n        combined_score = 0.7 * seq_score + 0.3 * struct_score\n        \n        similar_seqs.append((\n            target_id, \n            train_seq, \n            combined_score,\n            seq_score,\n            struct_score,\n            train_coords_dict[target_id]\n        ))\n    \n    # Sort by combined score\n    similar_seqs.sort(key=lambda x: x[2], reverse=True)\n    return similar_seqs[:top_n]","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def superimpose_structures(ref_coords, mobile_coords):\n    \"\"\"Superimpose mobile structure onto reference using Kabsch algorithm.\"\"\"\n    # Center both structures\n    ref_center = np.mean(ref_coords, axis=0)\n    mobile_center = np.mean(mobile_coords, axis=0)\n    \n    ref_centered = ref_coords - ref_center\n    mobile_centered = mobile_coords - mobile_center\n    \n    # Compute optimal rotation using SVD\n    H = mobile_centered.T @ ref_centered\n    U, S, Vt = np.linalg.svd(H)\n    \n    # Handle reflection case\n    d = np.sign(np.linalg.det(Vt.T @ U.T))\n    D = np.diag([1, 1, d])\n    \n    R_opt = Vt.T @ D @ U.T\n    \n    # Apply transformation\n    mobile_aligned = (mobile_centered @ R_opt) + ref_center\n    \n    return mobile_aligned\n\n\ndef ensemble_average_templates(adapted_structures, weights=None):\n    \"\"\"\n    Create ensemble average from multiple template-adapted structures.\n    \"\"\"\n    if len(adapted_structures) == 0:\n        return None\n    \n    if len(adapted_structures) == 1:\n        return adapted_structures[0]\n    \n    # Use first structure as reference\n    reference = adapted_structures[0]\n    n_residues = len(reference)\n    \n    if weights is None:\n        weights = np.ones(len(adapted_structures)) / len(adapted_structures)\n    else:\n        weights = np.array(weights)\n        weights = weights / weights.sum()\n    \n    # Superimpose all structures onto reference\n    aligned_structures = [reference]\n    for i, struct in enumerate(adapted_structures[1:], 1):\n        if len(struct) == n_residues:\n            aligned = superimpose_structures(reference, struct)\n            aligned_structures.append(aligned)\n        else:\n            # Skip if length mismatch\n            weights[i] = 0\n    \n    # Renormalize weights\n    weights = weights[:len(aligned_structures)]\n    weights = weights / weights.sum()\n    \n    # Compute weighted average\n    ensemble_coords = np.zeros((n_residues, 3))\n    for struct, weight in zip(aligned_structures, weights):\n        ensemble_coords += weight * struct\n    \n    return ensemble_coords","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compute_energy(coords, sequence, base_pairs, \n                   seq_dist_weight=1.0, clash_weight=5.0, bp_weight=0.5):\n    \"\"\"\n    Compute pseudo-energy for RNA structure.\n    \"\"\"\n    energy = 0.0\n    n = len(coords)\n    \n    # Sequential distance penalty\n    target_seq_dist = 5.9  # Target C1'-C1' distance\n    for i in range(n - 1):\n        dist = np.linalg.norm(coords[i+1] - coords[i])\n        energy += seq_dist_weight * (dist - target_seq_dist) ** 2\n    \n    # Steric clash penalty\n    min_dist = 3.8\n    for i in range(n):\n        for j in range(i + 2, n):\n            dist = np.linalg.norm(coords[j] - coords[i])\n            if dist < min_dist:\n                energy += clash_weight * (min_dist - dist) ** 2\n    \n    # Base pair distance penalty\n    target_bp_dist = 10.5  # Target base pair C1'-C1' distance\n    for i, j in base_pairs:\n        if i < n and j < n:\n            dist = np.linalg.norm(coords[j] - coords[i])\n            energy += bp_weight * (dist - target_bp_dist) ** 2\n    \n    return energy\n\n\ndef refine_coordinates_iterative(coords, sequence, base_pairs, n_iterations=50, learning_rate=0.1):\n    \"\"\"\n    Refine coordinates using gradient-based optimization.\n    \"\"\"\n    refined = coords.copy()\n    n = len(coords)\n    \n    target_seq_dist = 5.9\n    min_dist = 3.8\n    target_bp_dist = 10.5\n    \n    for iteration in range(n_iterations):\n        gradients = np.zeros_like(refined)\n        \n        # Sequential distance gradients\n        for i in range(n - 1):\n            diff = refined[i+1] - refined[i]\n            dist = np.linalg.norm(diff)\n            if dist > 1e-6:\n                grad = 2 * (dist - target_seq_dist) * (diff / dist)\n                gradients[i] -= grad * 0.5\n                gradients[i+1] += grad * 0.5\n        \n        # Steric clash gradients (only for severe clashes)\n        for i in range(n):\n            for j in range(i + 2, n):\n                diff = refined[j] - refined[i]\n                dist = np.linalg.norm(diff)\n                if dist < min_dist and dist > 1e-6:\n                    grad = 2 * 5.0 * (min_dist - dist) * (-diff / dist)\n                    gradients[i] += grad * 0.5\n                    gradients[j] -= grad * 0.5\n        \n        # Base pair gradients\n        for i, j in base_pairs:\n            if i < n and j < n:\n                diff = refined[j] - refined[i]\n                dist = np.linalg.norm(diff)\n                if dist > 1e-6:\n                    grad = 2 * 0.5 * (dist - target_bp_dist) * (diff / dist)\n                    gradients[i] -= grad * 0.3\n                    gradients[j] += grad * 0.3\n        \n        # Update coordinates with decreasing learning rate\n        lr = learning_rate * (0.95 ** iteration)\n        refined -= lr * gradients\n    \n    return refined","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def adapt_template_to_query(query_seq, template_seq, template_coords, alignment=None):\n    \"\"\"\n    Adapt template coordinates to fit the query sequence based on sequence alignment.\n    \"\"\"\n    if alignment is None:\n        query_seq_obj = Seq(query_seq)\n        template_seq_obj = Seq(template_seq)\n        alignments = pairwise2.align.globalms(query_seq_obj, template_seq_obj, 2, -1, -10, -0.5, one_alignment_only=True)\n        \n        if not alignments:\n            return generate_basic_structure(query_seq)\n            \n        alignment = alignments[0]\n    \n    aligned_query = alignment.seqA\n    aligned_template = alignment.seqB\n    \n    query_coords = np.zeros((len(query_seq), 3))\n    query_coords.fill(np.nan)\n    \n    query_idx = 0\n    template_idx = 0\n    \n    for i in range(len(aligned_query)):\n        query_char = aligned_query[i]\n        template_char = aligned_template[i]\n        \n        if query_char != '-' and template_char != '-':\n            if template_idx < len(template_coords):\n                query_coords[query_idx] = template_coords[template_idx]\n            template_idx += 1\n            query_idx += 1\n        elif query_char != '-' and template_char == '-':\n            query_idx += 1\n        elif query_char == '-' and template_char != '-':\n            template_idx += 1\n    \n    # Interpolate NaN positions\n    query_coords = interpolate_missing_coords(query_coords)\n    \n    return query_coords\n\n\ndef interpolate_missing_coords(coords):\n    \"\"\"Interpolate missing (NaN) coordinates.\"\"\"\n    n = len(coords)\n    typical_step = 5.9\n    \n    # First pass: interpolate between valid points\n    for i in range(n):\n        if np.isnan(coords[i, 0]):\n            prev_valid = next((j for j in range(i-1, -1, -1) if not np.isnan(coords[j, 0])), -1)\n            next_valid = next((j for j in range(i+1, n) if not np.isnan(coords[j, 0])), -1)\n            \n            if prev_valid >= 0 and next_valid >= 0:\n                weight = (i - prev_valid) / (next_valid - prev_valid)\n                coords[i] = (1 - weight) * coords[prev_valid] + weight * coords[next_valid]\n    \n    # Second pass: handle endpoints\n    for i in range(n):\n        if np.isnan(coords[i, 0]):\n            if i == 0:\n                first_valid = next((j for j in range(1, n) if not np.isnan(coords[j, 0])), -1)\n                if first_valid >= 0:\n                    direction = np.random.normal(0, 1, 3)\n                    direction = direction / (np.linalg.norm(direction) + 1e-10) * typical_step\n                    for j in range(first_valid-1, -1, -1):\n                        coords[j] = coords[j+1] - direction\n                else:\n                    for j in range(n):\n                        angle = j * 0.6\n                        coords[j] = [10.0 * np.cos(angle), 10.0 * np.sin(angle), j * 2.5]\n                    break\n            else:\n                prev_valid = next((j for j in range(i-1, -1, -1) if not np.isnan(coords[j, 0])), -1)\n                if prev_valid >= 0:\n                    if prev_valid > 0:\n                        direction = coords[prev_valid] - coords[prev_valid-1]\n                    else:\n                        direction = np.random.normal(0, 1, 3)\n                    direction = direction / (np.linalg.norm(direction) + 1e-10) * typical_step\n                    coords[i] = coords[prev_valid] + direction\n    \n    return np.nan_to_num(coords)\n\n\ndef generate_basic_structure(sequence):\n    \"\"\"Generate a simple helical structure.\"\"\"\n    n = len(sequence)\n    coords = np.zeros((n, 3))\n    \n    radius = 10.0\n    rise = 2.5\n    angle_step = 0.6\n    \n    for i in range(n):\n        angle = i * angle_step\n        coords[i] = [radius * np.cos(angle), radius * np.sin(angle), i * rise]\n    \n    return coords","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def predict_rna_structures_improved(sequence, target_id, train_seqs_df, train_coords_dict,\n                                    kmer_index, train_structures, n_predictions=5, \n                                    temporal_cutoff=None):\n    \"\"\"\n    Improved RNA structure prediction with ensemble averaging and refinement.\n    \"\"\"\n    predictions = []\n    \n    # Get secondary structure and base pairs\n    structure = predict_secondary_structure(sequence)\n    base_pairs = get_base_pairs(structure)\n    \n    # Find similar sequences\n    similar_seqs = find_similar_sequences_fast(\n        sequence, train_seqs_df, train_coords_dict,\n        kmer_index, train_structures, temporal_cutoff=temporal_cutoff, top_n=10\n    )\n    \n    # Adapt templates\n    adapted_templates = []\n    template_weights = []\n    \n    for template_id, template_seq, combined_score, seq_score, struct_score, template_coords in similar_seqs:\n        adapted = adapt_template_to_query(sequence, template_seq, template_coords)\n        if adapted is not None and not np.isnan(adapted).any():\n            adapted_templates.append(adapted)\n            template_weights.append(combined_score)\n    \n    if adapted_templates:\n        # Create ensemble average (first prediction)\n        ensemble_coords = ensemble_average_templates(adapted_templates, template_weights)\n        if ensemble_coords is not None:\n            # Refine ensemble\n            refined_ensemble = refine_coordinates_iterative(ensemble_coords, sequence, base_pairs)\n            predictions.append(refined_ensemble)\n        \n        # Add top individual templates with variations\n        for i, (adapted, weight) in enumerate(zip(adapted_templates[:4], template_weights[:4])):\n            # Refine each template\n            refined = refine_coordinates_iterative(adapted, sequence, base_pairs, n_iterations=30)\n            \n            # Add small random perturbation based on template quality\n            noise_scale = max(0.05, 0.5 - weight)\n            perturbed = refined + np.random.normal(0, noise_scale, refined.shape)\n            \n            predictions.append(perturbed)\n            \n            if len(predictions) >= n_predictions:\n                break\n    \n    # Fill remaining with de novo predictions\n    while len(predictions) < n_predictions:\n        seed = hash(target_id) % 10000 + len(predictions) * 1000\n        np.random.seed(seed)\n        random.seed(seed)\n        \n        de_novo = generate_rna_structure_improved(sequence, base_pairs)\n        refined_de_novo = refine_coordinates_iterative(de_novo, sequence, base_pairs)\n        predictions.append(refined_de_novo)\n    \n    return predictions[:n_predictions]\n\n\ndef generate_rna_structure_improved(sequence, base_pairs, seed=None):\n    \"\"\"\n    Generate RNA structure with secondary structure guidance.\n    \"\"\"\n    if seed is not None:\n        np.random.seed(seed)\n        random.seed(seed)\n    \n    n = len(sequence)\n    coords = np.zeros((n, 3))\n    \n    # Initialize with helix\n    for i in range(min(3, n)):\n        angle = i * 0.6\n        coords[i] = [10.0 * np.cos(angle), 10.0 * np.sin(angle), i * 2.5]\n    \n    # Create structure aware of base pairs\n    bp_dict = {i: j for i, j in base_pairs}\n    bp_dict.update({j: i for i, j in base_pairs})\n    \n    current_direction = np.array([0.0, 0.0, 1.0])\n    \n    for i in range(3, n):\n        # Check if this residue has a base pair\n        if i in bp_dict and bp_dict[i] < i:\n            # Position near base-paired residue\n            pair_idx = bp_dict[i]\n            pair_pos = coords[pair_idx]\n            \n            center = np.mean(coords[:i], axis=0)\n            direction = center - pair_pos\n            direction = direction / (np.linalg.norm(direction) + 1e-10)\n            \n            bp_dist = 10.5 + random.uniform(-1.0, 1.0)\n            offset = np.random.normal(0, 1, 3) * 1.5\n            coords[i] = pair_pos + direction * bp_dist + offset\n        else:\n            # Continue backbone\n            if random.random() < 0.3:\n                angle = random.uniform(0.2, 0.6)\n                axis = np.random.normal(0, 1, 3)\n                axis = axis / (np.linalg.norm(axis) + 1e-10)\n                rotation = R.from_rotvec(angle * axis)\n                current_direction = rotation.apply(current_direction)\n            else:\n                current_direction += np.random.normal(0, 0.15, 3)\n                current_direction = current_direction / (np.linalg.norm(current_direction) + 1e-10)\n            \n            step = random.uniform(5.5, 6.5)\n            coords[i] = coords[i-1] + step * current_direction\n    \n    return coords","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def calculate_rmsd(coords1, coords2):\n    \"\"\"Calculate RMSD between two coordinate sets after optimal superposition.\"\"\"\n    # Superimpose\n    aligned = superimpose_structures(coords1, coords2)\n    \n    # Calculate RMSD\n    diff = aligned - coords1\n    rmsd = np.sqrt(np.mean(np.sum(diff ** 2, axis=1)))\n    \n    return rmsd\n\n\ndef evaluate_on_validation(valid_seqs_df, valid_coords_dict, train_seqs_df, train_coords_dict,\n                          kmer_index, train_structures, n_samples=20):\n    \"\"\"\n    Evaluate prediction quality on validation set.\n    \"\"\"\n    rmsds = []\n    \n    # Sample validation sequences\n    sample_ids = random.sample(list(valid_coords_dict.keys()), min(n_samples, len(valid_coords_dict)))\n    \n    for target_id in sample_ids:\n        row = valid_seqs_df[valid_seqs_df['target_id'] == target_id]\n        if len(row) == 0:\n            continue\n            \n        sequence = row['sequence'].values[0]\n        true_coords = valid_coords_dict[target_id]\n        temporal_cutoff = row['temporal_cutoff'].values[0] if 'temporal_cutoff' in row else None\n        \n        # Generate predictions\n        predictions = predict_rna_structures_improved(\n            sequence, target_id, train_seqs_df, train_coords_dict,\n            kmer_index, train_structures, n_predictions=5, temporal_cutoff=temporal_cutoff\n        )\n        \n        # Calculate RMSD for best prediction\n        min_rmsd = float('inf')\n        for pred in predictions:\n            if len(pred) == len(true_coords):\n                rmsd = calculate_rmsd(true_coords, pred)\n                min_rmsd = min(min_rmsd, rmsd)\n        \n        if min_rmsd < float('inf'):\n            rmsds.append(min_rmsd)\n    \n    if rmsds:\n        print(f\"Validation Results (n={len(rmsds)}):\")\n        print(f\"  Mean RMSD: {np.mean(rmsds):.2f} Å\")\n        print(f\"  Median RMSD: {np.median(rmsds):.2f} Å\")\n        print(f\"  Min RMSD: {np.min(rmsds):.2f} Å\")\n        print(f\"  Max RMSD: {np.max(rmsds):.2f} Å\")\n        \n        plt.figure(figsize=(10, 5))\n        plt.hist(rmsds, bins=20, edgecolor='black', alpha=0.7)\n        plt.axvline(np.mean(rmsds), color='red', linestyle='--', label=f'Mean: {np.mean(rmsds):.2f} Å')\n        plt.xlabel('RMSD (Å)')\n        plt.ylabel('Count')\n        plt.title('RMSD Distribution on Validation Set')\n        plt.legend()\n        plt.grid(True, alpha=0.3)\n        plt.show()\n    \n    return rmsds\n\n# Run validation\nprint(\"Evaluating on validation set...\")\nvalidation_rmsds = evaluate_on_validation(\n    valid_seqs, valid_coords_dict, train_seqs, train_coords_dict,\n    kmer_index, train_structures, n_samples=30\n)","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Generate predictions for test set\nall_predictions = []\nstart_time = time.time()\ntotal_targets = len(test_seqs)\n\nprint(f\"Generating predictions for {total_targets} test sequences...\")\n\nfor idx, row in test_seqs.iterrows():\n    target_id = row['target_id']\n    sequence = row['sequence']\n    temporal_cutoff = row.get('temporal_cutoff', None)\n    \n    # Progress tracking\n    if idx % 10 == 0:\n        elapsed = time.time() - start_time\n        targets_processed = idx + 1\n        if targets_processed > 1:\n            avg_time = elapsed / targets_processed\n            remaining = avg_time * (total_targets - targets_processed)\n            print(f\"Processing {targets_processed}/{total_targets}: {target_id} ({len(sequence)} nt), \"\n                  f\"elapsed: {elapsed:.1f}s, remaining: {remaining:.1f}s\")\n    \n    # Generate predictions\n    predictions = predict_rna_structures_improved(\n        sequence, target_id, train_seqs, train_coords_dict,\n        kmer_index, train_structures, n_predictions=5, temporal_cutoff=temporal_cutoff\n    )\n    \n    # Store predictions\n    for j in range(len(sequence)):\n        pred_row = {\n            'ID': f\"{target_id}_{j+1}\",\n            'resname': sequence[j],\n            'resid': j + 1\n        }\n        \n        for i in range(5):\n            pred_row[f'x_{i+1}'] = predictions[i][j][0]\n            pred_row[f'y_{i+1}'] = predictions[i][j][1]\n            pred_row[f'z_{i+1}'] = predictions[i][j][2]\n        \n        all_predictions.append(pred_row)\n\n# Create submission DataFrame\nsubmission_df = pd.DataFrame(all_predictions)\n\n# Ensure correct column order\ncolumn_order = ['ID', 'resname', 'resid']\nfor i in range(1, 6):\n    for coord in ['x', 'y', 'z']:\n        column_order.append(f'{coord}_{i}')\nsubmission_df = submission_df[column_order]\n\n# Save submission\nsubmission_df.to_csv('submission.csv', index=False)\n\nprint(f\"\\nGenerated predictions for {len(test_seqs)} RNA sequences\")\nprint(f\"Total runtime: {time.time() - start_time:.1f} seconds\")\nprint(f\"Submission shape: {submission_df.shape}\")\nsubmission_df.head()","metadata":{},"outputs":[],"execution_count":null}]}