{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.12"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":87793,"databundleVersionId":11553390,"isSourceIdPinned":false,"sourceType":"competition"}],"dockerImageVersionId":30919,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true},"papermill":{"default_parameters":{},"duration":165.341765,"end_time":"2025-04-16T20:33:03.888498","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2025-04-16T20:30:18.546733","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"7019acc5","cell_type":"markdown","source":"# RNA 3D Structure Prediction Pipeline 🧬\n\n## Overview 📜\nThis project implements a comprehensive pipeline for the Stanford RNA 3D Folding competition, focusing on predicting the three-dimensional structure of RNA molecules from nucleotide sequences. The pipeline incorporates multiple strategies including reference-based modeling, neural network quality assessment, and specialized RNA structure generation techniques.\n\n## Key Components 🧩\n\n### 1. Data Processing and Management 📊\n\n**Memory-Optimized Data Loading**\n- Efficient loading for large RNA datasets\n- Coordinate normalization to handle numerical issues\n- Padding strategies for variable-length sequences\n\n**Feature Engineering**\n- One-hot encoding of RNA sequences (A, C, G, U, N)\n- Sequence composition analysis (GC/AU content)\n- Structure normalization and centralization\n\n### 2. Multi-Phase Hybrid Approach 🧮\n\n**Phase 1: Golden Seeds Discovery**\n- Systematic search for high-performing random seeds\n- TM-score based evaluation metrics\n- Known good seeds prioritization and random exploration\n\n**Phase 2: Quality Assessment Modeling**\n- Enhanced neural network for structure quality assessment\n- Multiple evaluation outputs (quality score, bond score, validity)\n- Attention mechanism for capturing long-range interactions\n- Rule-based quality assessment as fallback\n\n**Phase 3: Base Structure Generation**\n- RNA-specific optimization based on sequence properties\n- Size-dependent parameter tuning\n- Stem-loop template application based on sequence analysis\n- Structure validation and emergency generation\n\n**Phase 4: Diverse Structure Generation and Pruning**\n- Generation of multiple candidate structures\n- RNA size-specific variation parameters\n- Structure repair and backbone refinement\n- Neural network or rule-based quality pruning\n\n**Phase 5: Submission Creation**\n- Multiple submission formats (hybrid, NN pruned, reference)\n- Comprehensive logging and error handling\n- Performance statistics reporting\n\n### 3. Structure Generation and Sampling 🎯\n\n**Geometric Sampling with Physical Constraints**\n- Correlated noise for natural structural transitions\n- Bond length preservation (3.8 Å typical RNA backbone)\n- Global movement simulation for domain flexibility\n- Sequence-aware geometric constraints\n\n**RNA-Specific Structure Refinement**\n- Stem-loop template application\n- GC/AU content-based refinement\n- Backbone angle adjustment to RNA-specific values\n- Natural hinge point detection and rotation\n\n### 4. Quality Assessment Techniques ✅\n\n**Neural Network Quality Model**\n- Variable-length sequence handling\n- RNA-specific feature extraction\n- Pairwise distance calculations with 2D convolutions\n- Self-attention mechanism for long-range interactions\n- Multi-output prediction (quality, bond quality, validity)\n\n**Rule-Based Quality Assessment**\n- Biophysical validation checks\n- Bond distance and consistency analysis\n- Structure validity verification\n- Radius of gyration and compactness evaluation\n\n### 5. Flexible Execution Modes 🔄\n\n**Hybrid Pipeline (Golden Seeds + NN Pruning)**\n- Combined approach using all pipeline phases\n- Golden seed discovery and quality model training\n- Diverse structure generation and neural pruning\n\n**NN Pruning Only Mode**\n- Reference model with neural network pruning\n- Simplified workflow without golden seed search\n- Quality-based selection of top structures\n\n**Reference-Only Mode**\n- Basic approach using only reference model\n- Simple variations for structure diversity\n- Efficient baseline performance\n\n## Methodology 🔍\n\nThe pipeline employs a hybrid approach combining reference-based modeling with neural network quality assessment:\n\n1. **Data Preparation**: Sequences are converted to one-hot encoding and structures are normalized to ensure numerical stability.\n\n2. **Golden Seeds Discovery**: The pipeline searches for optimal random seeds that produce high-quality base predictions as measured by TM-score on validation data.\n\n3. **Quality Model Training**: An enhanced neural network is trained to assess RNA structure quality, capturing both local features (bond lengths, angles) and global characteristics (overall fold quality).\n\n4. **Size-Adaptive Strategy**: Different parameters and generation strategies are applied based on RNA size:\n - Small RNAs (<50 residues): Higher structural diversity with controlled variations\n - Medium RNAs (50-120 residues): Balanced approach with moderate variations\n - Large RNAs (>120 residues): Conservative variations with optimized parameters\n\n5. **Structure Generation and Refinement**:\n - Base structures are generated using golden seeds\n - RNA-specific templates are applied based on sequence analysis\n - Multiple candidate structures are generated with varying parameters\n - Neural network or rule-based quality assessment prunes to top structures\n\n6. **Submission Creation**: Final structures are compiled into the required submission format with comprehensive error handling and fallback mechanisms.\n\nThe approach balances computational efficiency with structural accuracy, focusing on generating biologically plausible RNA structures that maintain essential physical constraints while exploring the conformational space effectively.","metadata":{"papermill":{"duration":0.012017,"end_time":"2025-04-16T20:30:22.185063","exception":false,"start_time":"2025-04-16T20:30:22.173046","status":"completed"},"tags":[]}},{"id":"85356f47","cell_type":"markdown","source":"## Library Imports 📚🔧","metadata":{"papermill":{"duration":0.010073,"end_time":"2025-04-16T20:30:22.205792","exception":false,"start_time":"2025-04-16T20:30:22.195719","status":"completed"},"tags":[]}},{"id":"26e10983","cell_type":"code","source":"# Standard Library Imports\nimport os\nimport time\nimport gc\nimport traceback\nfrom collections import Counter\nimport warnings\nimport hashlib\n\n# Data Manipulation Libraries\nimport numpy as np\nimport pandas as pd\n\n# Visualization Libraries\nimport matplotlib.pyplot as plt\nimport matplotlib.colors as mcolors\n\n# Machine Learning Libraries\ntry:\n # TensorFlow and Keras\n import tensorflow as tf\n from tensorflow.keras import layers, models, optimizers\n from tensorflow.keras.models import Model\n from tensorflow.keras.layers import (\n Input, Conv1D, Dense, Dropout, BatchNormalization, \n Flatten, Reshape, Bidirectional, LSTM\n )\n from tensorflow.keras.callbacks import EarlyStopping\n \n # Scikit-learn\n from sklearn.model_selection import train_test_split\n \n # XGBoost\n import xgboost as xgb\n \n ML_AVAILABLE = True\nexcept ImportError:\n print(\"Warning: ML libraries not available. Will use only reference-based methods.\")\n ML_AVAILABLE = False\n\n# Set random seed for reproducibility\nnp.random.seed(0)\n\n# Suppress warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:20:53.828081Z","iopub.execute_input":"2025-05-08T13:20:53.828451Z","iopub.status.idle":"2025-05-08T13:21:07.892296Z","shell.execute_reply.started":"2025-05-08T13:20:53.828423Z","shell.execute_reply":"2025-05-08T13:21:07.891553Z"},"papermill":{"duration":21.26111,"end_time":"2025-04-16T20:30:43.476906","exception":false,"start_time":"2025-04-16T20:30:22.215796","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"36780eb5","cell_type":"markdown","source":"## 🧬 RNA 3D Structure Prediction and Analysis Pipeline 🔬","metadata":{"papermill":{"duration":0.010139,"end_time":"2025-04-16T20:30:43.498348","exception":false,"start_time":"2025-04-16T20:30:43.488209","status":"completed"},"tags":[]}},{"id":"f1ecb8a7","cell_type":"code","source":"# Directories and files adjusted for the new competition\nDATA_DIR = os.getenv('DATA_DIR', '/kaggle/input/stanford-rna-3d-folding/')\nmain_files = [\n \"train_sequences.csv\", \n \"train_labels.csv\", \n \"validation_sequences.csv\", \n \"validation_labels.csv\", \n \"test_sequences.csv\",\n \"sample_submission.csv\" \n]\n\nDEFAULT_THRESHOLD = 0.4 # Default threshold after analysis\n\ndef optimize_dataframe(df, inplace=False, category_threshold=DEFAULT_THRESHOLD):\n \"\"\"\n Optimizes the DataFrame to save memory.\n \"\"\"\n if category_threshold < 0 or category_threshold > 1:\n raise ValueError(\"category_threshold must be between 0 and 1.\")\n \n if not inplace:\n df = df.copy()\n \n for col in df.columns:\n col_type = df[col].dtype\n if np.issubdtype(col_type, np.integer):\n c_min, c_max = df[col].min(), df[col].max()\n if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n df[col] = df[col].astype(np.int8)\n elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n df[col] = df[col].astype(np.int16)\n elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n df[col] = df[col].astype(np.int32)\n elif np.issubdtype(col_type, np.floating):\n if df[col].min() > np.finfo(np.float32).min and df[col].max() < np.finfo(np.float32).max:\n df[col] = df[col].astype(np.float32)\n if col_type == object:\n unique_vals = len(df[col].unique())\n if unique_vals / len(df) < category_threshold:\n df[col] = df[col].astype('category')\n \n return df\n\ndef load_main_data(chunksize=50000):\n \"\"\"\n Loads the main files.\n \"\"\"\n data = {}\n for file_name in main_files:\n file_path = os.path.join(DATA_DIR, file_name)\n if os.path.exists(file_path):\n chunks = pd.read_csv(file_path, on_bad_lines='skip', low_memory=False, chunksize=chunksize)\n dataframes = [optimize_dataframe(chunk, category_threshold=DEFAULT_THRESHOLD) for chunk in chunks]\n data[file_name] = pd.concat(dataframes, ignore_index=True)\n else:\n print(f\"File {file_path} not found!\")\n return data\n\ndef check_data_integrity(original_df, optimized_df):\n \"\"\"\n Checks if the optimization did not alter the data.\n \"\"\"\n try:\n pd.testing.assert_frame_equal(original_df, optimized_df, check_like=True)\n print(\"Integrity check passed: No changes in data after optimization.\")\n except AssertionError as e:\n print(f\"Data integrity check failed: {e}\")\n\ndef check_duplicates(df):\n \"\"\"\n Checks for duplicates in the DataFrame.\n \"\"\"\n duplicates = df[df.duplicated(keep=False)]\n if not duplicates.empty:\n print(f\"Warning: Duplicates found in the dataset. Number of duplicates: {duplicates.shape[0]}\")\n return duplicates\n else:\n print(\"No duplicates found.\")\n return None\n\ndef test_thresholds(df):\n \"\"\"\n Tests different thresholds for DataFrame optimization.\n \"\"\"\n thresholds = np.linspace(0.1, 0.9, 9)\n memory_usages = []\n for threshold in thresholds:\n optimized_df = optimize_dataframe(df.copy(), category_threshold=threshold)\n memory_usages.append(optimized_df.memory_usage(deep=True).sum() / 1024**2)\n return thresholds, memory_usages\n\ndef plot_memory_usage(thresholds, memory_usages):\n \"\"\"\n Plots memory usage versus thresholds.\n \"\"\"\n plt.figure(figsize=(10, 6))\n plt.plot(thresholds, memory_usages, marker='o', linestyle='-')\n plt.title(\"Memory Usage vs. Threshold\")\n plt.xlabel(\"Threshold\")\n plt.ylabel(\"Memory Usage (MB)\")\n plt.grid(True)\n plt.show()\n\ndef analyze_sequence_data(df_sequences):\n \"\"\"\n Analyzes RNA sequence data.\n \"\"\"\n # Basic information\n print(f\"Total sequences: {len(df_sequences)}\")\n print(f\"Available columns: {df_sequences.columns.tolist()}\")\n \n # Sequence analysis\n if 'sequence' in df_sequences.columns:\n # Distribution of sequence lengths\n seq_lengths = df_sequences['sequence'].apply(len)\n print(f\"\\nSequence length statistics:\")\n print(f\"Minimum: {seq_lengths.min()}\")\n print(f\"Maximum: {seq_lengths.max()}\")\n print(f\"Average: {seq_lengths.mean():.2f}\")\n \n # Nucleotide count\n nucleotides = ['A', 'C', 'G', 'U']\n nucleotide_counts = {n: df_sequences['sequence'].str.count(n).sum() for n in nucleotides}\n total_nucleotides = sum(nucleotide_counts.values())\n \n print(\"\\nNucleotide distribution:\")\n for n, count in nucleotide_counts.items():\n print(f\"{n}: {count} ({count/total_nucleotides*100:.2f}%)\")\n \n return df_sequences\n\ndef analyze_label_data(df_labels):\n \"\"\"\n Analyzes 3D coordinate data (labels).\n \"\"\"\n print(f\"Total entries in labels: {len(df_labels)}\")\n print(f\"Available columns: {df_labels.columns.tolist()}\")\n \n # Analysis of 3D coordinates if available\n coord_columns = [col for col in df_labels.columns if col.startswith(('x_', 'y_', 'z_'))]\n if coord_columns:\n print(f\"\\nCoordinate columns found: {len(coord_columns)}\")\n \n # Basic statistics of coordinates\n for i in range(1, 6): # For the 5 possible structures\n x_col = f'x_{i}'\n y_col = f'y_{i}'\n z_col = f'z_{i}'\n \n if x_col in df_labels.columns and y_col in df_labels.columns and z_col in df_labels.columns:\n print(f\"\\nStatistics for structure {i}:\")\n print(f\"X - Mean: {df_labels[x_col].mean():.2f}, Std: {df_labels[x_col].std():.2f}\")\n print(f\"Y - Mean: {df_labels[y_col].mean():.2f}, Std: {df_labels[y_col].std():.2f}\")\n print(f\"Z - Mean: {df_labels[z_col].mean():.2f}, Std: {df_labels[z_col].std():.2f}\")\n \n return df_labels\n\ndef create_submission_template(test_df, sample_submission_df):\n \"\"\"\n Creates a submission template based on test data.\n \"\"\"\n # Check if sample_submission.csv is available\n if sample_submission_df is None:\n print(\"Sample submission file not found. Creating a new template.\")\n \n # Create a new DataFrame for submission\n submission_df = pd.DataFrame()\n \n # Example code to fill the template (adjust as needed)\n ids = []\n resnames = []\n resids = []\n \n for _, row in test_df.iterrows():\n sequence = row['sequence']\n target_id = row['target_id']\n \n for i, nucleotide in enumerate(sequence, 1):\n ids.append(f\"{target_id}_{i}\")\n resnames.append(nucleotide)\n resids.append(i)\n \n submission_df['ID'] = ids\n submission_df['resname'] = resnames\n submission_df['resid'] = resids\n \n # Add coordinate columns (5 structures)\n for i in range(1, 6):\n submission_df[f'x_{i}'] = 0.0\n submission_df[f'y_{i}'] = 0.0\n submission_df[f'z_{i}'] = 0.0\n else:\n submission_df = sample_submission_df.copy()\n print(\"Submission template created based on the provided example.\")\n \n return submission_df\n\ndef main():\n start_time = time.time()\n \n # Load main data\n print(\"Loading main data...\")\n main_data = load_main_data()\n \n # Check which files were loaded\n print(\"\\nLoaded files:\")\n for file_name, df in main_data.items():\n print(f\"- {file_name}: {df.shape if df is not None else 'Not found'}\")\n \n # Analyze training sequence data\n if \"train_sequences.csv\" in main_data:\n print(\"\\n===== Training Sequences Analysis =====\")\n analyze_sequence_data(main_data[\"train_sequences.csv\"])\n \n # Analyze training label data\n if \"train_labels.csv\" in main_data:\n print(\"\\n===== Training Labels Analysis =====\")\n analyze_label_data(main_data[\"train_labels.csv\"])\n \n # Check for duplicates in training data\n if \"train_sequences.csv\" in main_data:\n print(\"\\nChecking for duplicates in training sequences...\")\n check_duplicates(main_data[\"train_sequences.csv\"])\n \n # Create submission template\n if \"test_sequences.csv\" in main_data:\n print(\"\\nCreating submission template...\")\n submission_template = create_submission_template(\n main_data[\"test_sequences.csv\"],\n main_data.get(\"sample_submission.csv\")\n )\n print(f\"Submission template shape: {submission_template.shape}\")\n print(f\"First rows of the submission template:\")\n print(submission_template.head())\n \n # Calculate execution time\n end_time = time.time()\n print(f\"\\nRuntime: {end_time - start_time:.2f} seconds\")\n \n return main_data\n\nif __name__ == '__main__':\n main_data = main()","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.status.busy":"2025-05-08T13:21:07.893564Z","iopub.execute_input":"2025-05-08T13:21:07.894187Z","iopub.status.idle":"2025-05-08T13:21:08.449949Z","shell.execute_reply.started":"2025-05-08T13:21:07.894132Z","shell.execute_reply":"2025-05-08T13:21:08.449137Z"},"papermill":{"duration":0.685644,"end_time":"2025-04-16T20:30:44.194224","exception":false,"start_time":"2025-04-16T20:30:43.508580","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"70c7b540","cell_type":"markdown","source":"## Directory Explorer & CSV Verification for RNA3D 🗂️🔬","metadata":{"papermill":{"duration":0.010176,"end_time":"2025-04-16T20:30:44.216001","exception":false,"start_time":"2025-04-16T20:30:44.205825","status":"completed"},"tags":[]}},{"id":"103a46e5","cell_type":"code","source":"# Updated main directory\ndir_main = \"/kaggle/input/stanford-rna-3d-folding/\"\n\n# List all files and directories in the main directory\ntry:\n all_files = os.listdir(dir_main)\n print(f\"All files and directories in '{dir_main}':\")\n \n for file in all_files:\n # Check if it's a file or directory\n full_path = os.path.join(dir_main, file)\n type_desc = \"directory\" if os.path.isdir(full_path) else \"file\"\n size = os.path.getsize(full_path) / 1024 # Size in KB\n print(f\" - {file} ({type_desc}, {size:.2f} KB)\")\n \n # If it's a directory, list up to 5 files inside it\n if os.path.isdir(full_path):\n try:\n internal_files = os.listdir(full_path)[:5] # Limit to 5 files\n if internal_files:\n print(f\" First files in '{file}':\")\n for internal_file in internal_files:\n print(f\" * {internal_file}\")\n if len(os.listdir(full_path)) > 5:\n print(f\" * ... and {len(os.listdir(full_path)) - 5} more file(s)\")\n else:\n print(f\" '{file}' is empty\")\n except Exception as e:\n print(f\" Error listing contents of '{file}': {e}\")\nexcept Exception as e:\n print(f\"Error listing directory {dir_main}: {e}\")\n\n# Check the structure of the main CSV files\nmain_files = [\n \"train_sequences.csv\", \n \"train_labels.csv\", \n \"validation_sequences.csv\", \n \"validation_labels.csv\", \n \"test_sequences.csv\",\n \"sample_submission.csv\"\n]\nprint(\"\\nChecking main CSV files:\")\n\nfor file in main_files:\n full_path = os.path.join(dir_main, file)\n if os.path.exists(full_path):\n # Get file size\n size_mb = os.path.getsize(full_path) / (1024 * 1024) # Size in MB\n \n # Read the first lines to check the structure\n try:\n import pandas as pd\n df = pd.read_csv(full_path, nrows=1)\n print(f\"\\n{file} ({size_mb:.2f} MB):\")\n print(f\"Columns: {df.columns.tolist()}\")\n print(f\"Example:\")\n print(df.head())\n except Exception as e:\n print(f\"Error reading {file}: {e}\")\n else:\n print(f\"{file} not found.\")","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:08.451748Z","iopub.execute_input":"2025-05-08T13:21:08.451998Z","iopub.status.idle":"2025-05-08T13:21:08.536189Z","shell.execute_reply.started":"2025-05-08T13:21:08.451966Z","shell.execute_reply":"2025-05-08T13:21:08.535530Z"},"papermill":{"duration":0.094585,"end_time":"2025-04-16T20:30:44.320783","exception":false,"start_time":"2025-04-16T20:30:44.226198","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"c1c9992d","cell_type":"markdown","source":"## RNA3D Data Checker 🔍🧬","metadata":{"papermill":{"duration":0.011966,"end_time":"2025-04-16T20:30:44.345014","exception":false,"start_time":"2025-04-16T20:30:44.333048","status":"completed"},"tags":[]}},{"id":"be34440c","cell_type":"code","source":"# Updated main directory\ndir_main = \"/kaggle/input/stanford-rna-3d-folding/\"\n\ndef load_data():\n \"\"\"\n Loads the main CSV files from the Stanford RNA 3D Folding competition.\n Returns a dictionary with DataFrames.\n \"\"\"\n main_files = [\n \"train_sequences.csv\", \n \"train_labels.csv\", \n \"validation_sequences.csv\", \n \"validation_labels.csv\", \n \"test_sequences.csv\",\n \"sample_submission.csv\"\n ]\n \n data = {}\n for file_name in main_files:\n file_path = os.path.join(dir_main, file_name)\n if os.path.exists(file_path):\n try:\n data[file_name] = pd.read_csv(file_path)\n print(f\"File {file_name} loaded successfully. Shape: {data[file_name].shape}\")\n except Exception as e:\n print(f\"Error loading {file_name}: {e}\")\n else:\n print(f\"File {file_name} not found.\")\n data[file_name] = None\n \n return data\n\ndef compare_columns(main_data):\n \"\"\"\n Compares columns between different DataFrames.\n \"\"\"\n # List all available keys\n print(\"\\nLoaded files:\")\n print(list(main_data.keys()))\n \n # Compare columns between train_sequences.csv and test_sequences.csv\n if \"train_sequences.csv\" in main_data and \"test_sequences.csv\" in main_data:\n train_cols = set(main_data[\"train_sequences.csv\"].columns)\n test_cols = set(main_data[\"test_sequences.csv\"].columns)\n \n print(\"\\nColumns in train_sequences.csv:\")\n print(list(main_data[\"train_sequences.csv\"].columns))\n \n print(\"\\nUnique columns in train_sequences.csv (not present in test_sequences.csv):\")\n print(train_cols - test_cols)\n \n print(\"\\nUnique columns in test_sequences.csv (not present in train_sequences.csv):\")\n print(test_cols - train_cols)\n \n # Compare columns between train_labels.csv and validation_labels.csv\n if \"train_labels.csvv\" in main_data and \"validation_labels.csv\" in main_data:\n train_label_cols = set(main_data[\"train_labels.csv\"].columns)\n val_label_cols = set(main_data[\"validation_labels.csv\"].columns)\n \n print(\"\\nColumns in train_labels.csv:\")\n print(list(main_data[\"train_labels.csv\"].columns))\n \n print(\"\\nColumns in validation_labels.csv:\")\n print(list(main_data[\"validation_labels.csv\"].columns))\n \n print(\"\\nUnique columns in validation_labels.csv (not present in train_labels.csv):\")\n print(val_label_cols - train_label_cols)\n \n # Compare columns between validation_labels.csv and sample_submission.csv\n if \"validation_labels.csv\" in main_data and \"sample_submission.csv\" in main_data:\n val_label_cols = set(main_data[\"validation_labels.csv\"].columns)\n sample_cols = set(main_data[\"sample_submission.csv\"].columns)\n \n print(\"\\nColumns in sample_submission.csv:\")\n print(list(main_data[\"sample_submission.csv\"].columns))\n \n print(\"\\nUnique columns in validation_labels.csv (not present in sample_submission.csv):\")\n print(val_label_cols - sample_cols)\n \n print(\"\\nUnique columns in sample_submission.csv (not present in validation_labels.csv):\")\n print(sample_cols - val_label_cols)\n\ndef analyze_structure_format(main_data):\n \"\"\"\n Analyzes the format of 3D structures (coordinates).\n \"\"\"\n if \"validation_labels.csv\" in main_data and main_data[\"validation_labels.csv\"] is not None:\n df = main_data[\"validation_labels.csv\"]\n \n # Find all coordinate columns (x_1, y_1, z_1, etc.)\n coord_cols = [col for col in df.columns if col.startswith(('x_', 'y_', 'z_'))]\n \n # Group by structure\n structures = {}\n for col in coord_cols:\n # Extract structure number (e.g., \"x_1\" -> 1)\n parts = col.split('_')\n if len(parts) == 2:\n struct_num = int(parts[1])\n coord_type = parts[0]\n \n if struct_num not in structures:\n structures[struct_num] = []\n \n structures[struct_num].append(col)\n \n print(\"\\nStructure of the labels file:\")\n print(f\"Total structures found: {len(structures)}\")\n \n # Show details of the first structure\n if structures:\n first_struct = min(structures.keys())\n print(f\"\\nDetails of structure {first_struct}:\")\n print(f\"Columns: {sorted(structures[first_struct])}\")\n \n # Check for missing values\n for col in structures[first_struct]:\n missing = df[col].isna().sum()\n total = len(df)\n print(f\"{col}: {missing} missing values ({missing/total*100:.2f}%)\")\n \n # Check the range of non-missing values for the first structure\n for col in structures[first_struct]:\n non_null = df[col][df[col] != -1.0e+18] # Values that are not -1.0e+18\n if not non_null.empty:\n print(f\"{col} - Range: [{non_null.min():.3f}, {non_null.max():.3f}]\")\n\ndef main():\n # Load the data\n main_data = load_data()\n \n # Compare columns between different files\n compare_columns(main_data)\n \n # Analyze the format of 3D structures\n analyze_structure_format(main_data)\n \n return main_data\n\nif __name__ == '__main__':\n main_data = main()","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:08.537345Z","iopub.execute_input":"2025-05-08T13:21:08.537549Z","iopub.status.idle":"2025-05-08T13:21:08.794225Z","shell.execute_reply.started":"2025-05-08T13:21:08.537531Z","shell.execute_reply":"2025-05-08T13:21:08.793490Z"},"papermill":{"duration":0.27464,"end_time":"2025-04-16T20:30:44.631464","exception":false,"start_time":"2025-04-16T20:30:44.356824","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"123165af","cell_type":"markdown","source":"## Integrated RNA3D Sequence and Structure Analyzer 🔬🧬","metadata":{"papermill":{"duration":0.010264,"end_time":"2025-04-16T20:30:44.652653","exception":false,"start_time":"2025-04-16T20:30:44.642389","status":"completed"},"tags":[]}},{"id":"959de8f4","cell_type":"code","source":"# Initialize seed to control randomness\nnp.random.seed(0)\n\n# Directories and files adjusted for the new competition\nDATA_DIR = os.getenv('DATA_DIR', '/kaggle/input/stanford-rna-3d-folding/')\nmain_files = [\n \"train_sequences.csv\", \n \"train_labels.csv\", \n \"validation_sequences.csv\", \n \"validation_labels.csv\", \n \"test_sequences.csv\",\n \"sample_submission.csv\"\n]\n\nDEFAULT_THRESHOLD = 0.4 # Default threshold after analysis\n\ndef optimize_dataframe(df, inplace=False, category_threshold=DEFAULT_THRESHOLD):\n \"\"\"\n Optimizes the DataFrame to save memory.\n \"\"\"\n if category_threshold < 0 or category_threshold > 1:\n raise ValueError(\"category_threshold must be between 0 and 1.\")\n \n if not inplace:\n df = df.copy()\n \n for col in df.columns:\n col_type = df[col].dtype\n if np.issubdtype(col_type, np.integer):\n c_min, c_max = df[col].min(), df[col].max()\n if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n df[col] = df[col].astype(np.int8)\n elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n df[col] = df[col].astype(np.int16)\n elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n df[col] = df[col].astype(np.int32)\n elif np.issubdtype(col_type, np.floating):\n # First check if it's not the special value -1.0e+18\n if df[col].min() > np.finfo(np.float32).min and df[col].max() < np.finfo(np.float32).max:\n df[col] = df[col].astype(np.float32)\n if col_type == object:\n unique_vals = len(df[col].unique())\n if unique_vals / len(df) < category_threshold:\n df[col] = df[col].astype('category')\n \n return df\n\ndef load_main_data(chunksize=50000):\n \"\"\"\n Loads the main files.\n \"\"\"\n data = {}\n for file_name in main_files:\n file_path = os.path.join(DATA_DIR, file_name)\n if os.path.exists(file_path):\n chunks = pd.read_csv(file_path, on_bad_lines='skip', low_memory=False, chunksize=chunksize)\n dataframes = [optimize_dataframe(chunk, category_threshold=DEFAULT_THRESHOLD) for chunk in chunks]\n data[file_name] = pd.concat(dataframes, ignore_index=True)\n print(f\"File {file_name} loaded successfully. Shape: {data[file_name].shape}\")\n else:\n print(f\"File {file_path} not found!\")\n data[file_name] = None\n return data\n\ndef filter_columns_by_prefix(df, prefix=\"x_\"):\n \"\"\"\n Filters and counts the number of columns in a DataFrame based on a provided prefix.\n \n :param df: DataFrame where filtering will be applied.\n :param prefix: Prefix to be used for filtering. Ex: \"x_\", \"y_\", \"z_\".\n :return: List of filtered columns.\n \"\"\"\n filtered_columns = [col for col in df.columns if col.startswith(prefix)]\n return filtered_columns\n\ndef count_nucleotides(df, column_name='sequence'):\n \"\"\"\n Counts the frequency of each nucleotide in a specific column of a DataFrame.\n \n :param df: DataFrame containing the sequences.\n :param column_name: Name of the column containing the sequences. Default is 'sequence'.\n :return: Counter object with the nucleotide counts.\n \"\"\"\n from collections import Counter\n\n # Check if the column exists in the DataFrame\n if column_name not in df.columns:\n raise ValueError(f\"Column '{column_name}' not found in DataFrame.\")\n \n # Concatenate all sequences and count nucleotides\n all_sequences = ''.join(df[column_name].tolist())\n nucleotide_counts = Counter(all_sequences)\n \n return nucleotide_counts\n\ndef get_columns_without_missing_values(df):\n \"\"\"\n Returns columns without any missing values in the DataFrame.\n \n :param df: DataFrame to be checked.\n :return: List of columns without missing values.\n \"\"\"\n missing_values = df.isnull().sum()\n return missing_values[missing_values == 0].index.tolist()\n\ndef get_empty_columns(df):\n \"\"\"\n Returns columns that are completely empty in the DataFrame.\n \n :param df: DataFrame to be checked.\n :return: List of empty columns.\n \"\"\"\n missing_values = df.isnull().sum()\n return missing_values[missing_values == df.shape[0]].index.tolist()\n\ndef plot_coord_distributions(df_labels, prefix='x_', max_structures=5):\n \"\"\"\n Plots the distribution of coordinates (x, y, or z) for up to max_structures structures.\n \n :param df_labels: DataFrame containing the coordinates.\n :param prefix: Prefix of columns to be plotted ('x_', 'y_', or 'z_').\n :param max_structures: Maximum number of structures to show.\n \"\"\"\n # Find coordinate columns with the specified prefix\n coord_cols = filter_columns_by_prefix(df_labels, prefix)\n \n # Limit to the maximum number of structures\n coord_cols = sorted(coord_cols)[:max_structures]\n \n if not coord_cols:\n print(f\"No column with prefix '{prefix}' found.\")\n return\n \n # Set up the plot\n fig, axes = plt.subplots(1, len(coord_cols), figsize=(16, 4))\n if len(coord_cols) == 1:\n axes = [axes] # Ensure axes is iterable even with a single subplot\n \n # Plot histograms for each column\n for i, col in enumerate(coord_cols):\n # Filter special values (-1.0e+18) if present\n values = df_labels[col]\n filtered_values = values[values > -1.0e+17] # Cutoff value to filter -1.0e+18\n \n axes[i].hist(filtered_values, bins=30, alpha=0.7)\n axes[i].set_title(f'Distribution of {col}')\n axes[i].set_xlabel('Value')\n axes[i].set_ylabel('Frequency')\n \n plt.tight_layout()\n plt.show()\n\ndef analyze_3d_structure(df_labels):\n \"\"\"\n Analyzes the 3D coordinates of RNA structures.\n \n :param df_labels: DataFrame containing 3D coordinates.\n \"\"\"\n # Find all coordinate columns\n x_cols = filter_columns_by_prefix(df_labels, 'x_')\n y_cols = filter_columns_by_prefix(df_labels, 'y_')\n z_cols = filter_columns_by_prefix(df_labels, 'z_')\n \n print(f\"Number of x columns: {len(x_cols)}\")\n print(f\"Number of y columns: {len(y_cols)}\")\n print(f\"Number of z columns: {len(z_cols)}\")\n \n # Check for missing or special values in coordinates\n special_value = -1.0e+18 # Special value observed in the data\n \n for i, (x_col, y_col, z_col) in enumerate(zip(x_cols, y_cols, z_cols), 1):\n # Count missing or special values\n x_special = (df_labels[x_col] == special_value).sum()\n y_special = (df_labels[y_col] == special_value).sum()\n z_special = (df_labels[z_col] == special_value).sum()\n \n x_null = df_labels[x_col].isnull().sum()\n y_null = df_labels[y_col].isnull().sum()\n z_null = df_labels[z_col].isnull().sum()\n \n # Count how many complete structures exist (all x, y, z are neither special nor null)\n valid_structures = ((df_labels[x_col] != special_value) & \n (df_labels[y_col] != special_value) & \n (df_labels[z_col] != special_value) &\n df_labels[x_col].notnull() & \n df_labels[y_col].notnull() & \n df_labels[z_col].notnull()).sum()\n \n total_rows = len(df_labels)\n \n print(f\"\\nStructure {i}:\")\n print(f\" Special values: x={x_special} ({x_special/total_rows*100:.2f}%), y={y_special} ({y_special/total_rows*100:.2f}%), z={z_special} ({z_special/total_rows*100:.2f}%)\")\n print(f\" Null values: x={x_null} ({x_null/total_rows*100:.2f}%), y={y_null} ({y_null/total_rows*100:.2f}%), z={z_null} ({z_null/total_rows*100:.2f}%)\")\n print(f\" Complete structures: {valid_structures} ({valid_structures/total_rows*100:.2f}%)\")\n \n # Limit analysis to the first 5 structures\n if i >= 5:\n print(\"\\nAnalysis limited to the first 5 structures.\")\n break\n\ndef analyze_sequences(df_sequences):\n \"\"\"\n Analyzes RNA sequences.\n \n :param df_sequences: DataFrame containing the 'sequence' column.\n \"\"\"\n # Basic statistics of the sequence column\n print(\"\\nBasic statistics of the 'sequence' column:\")\n print(df_sequences['sequence'].describe())\n \n # Sequence lengths\n seq_lengths = df_sequences['sequence'].apply(len)\n print(\"\\nSequence length statistics:\")\n print(f\"Minimum: {seq_lengths.min()}\")\n print(f\"Maximum: {seq_lengths.max()}\")\n print(f\"Mean: {seq_lengths.mean():.2f}\")\n print(f\"Median: {seq_lengths.median()}\")\n \n # Nucleotide counts\n nucleotide_counts = count_nucleotides(df_sequences)\n total_nucleotides = sum(nucleotide_counts.values())\n \n print(\"\\nNucleotide distribution:\")\n for nucleotide, count in sorted(nucleotide_counts.items()):\n print(f\"{nucleotide}: {count} ({count/total_nucleotides*100:.2f}%)\")\n \n # Plot length distribution\n plt.figure(figsize=(10, 6))\n plt.hist(seq_lengths, bins=30, alpha=0.7)\n plt.title('Sequence Length Distribution')\n plt.xlabel('Length')\n plt.ylabel('Frequency')\n plt.grid(True, alpha=0.3)\n plt.show()\n\ndef main():\n # Load main data\n main_data = load_main_data()\n\n # Check which files were loaded\n print(\"\\nLoaded files:\")\n for file_name, df in main_data.items():\n if df is not None:\n print(f\"- {file_name}: {df.shape}\")\n \n # Analyze 3D structures in validation_labels.csv\n if \"validation_labels.csv\" in main_data and main_data[\"validation_labels.csv\"] is not None:\n print(\"\\n===== Analysis of 3D Structures (validation_labels.csv) =====\")\n df_labels = main_data[\"validation_labels.csv\"]\n \n # Count coordinate columns\n x_cols = filter_columns_by_prefix(df_labels, 'x_')\n y_cols = filter_columns_by_prefix(df_labels, 'y_')\n z_cols = filter_columns_by_prefix(df_labels, 'z_')\n \n print(f\"There are {len(x_cols)} x_ columns in the DataFrame.\")\n print(f\"There are {len(y_cols)} y_ columns in the DataFrame.\")\n print(f\"There are {len(z_cols)} z_ columns in the DataFrame.\")\n \n # Identify columns without missing values\n columns_without_missing = get_columns_without_missing_values(df_labels)\n print(f\"\\nColumns without missing values: {len(columns_without_missing)}\")\n \n # Identify completely empty columns\n empty_columns = get_empty_columns(df_labels)\n print(f\"Completely empty columns: {len(empty_columns)}\")\n \n # Analyze 3D coordinates in detail\n analyze_3d_structure(df_labels)\n \n # Plot distribution of x, y, z coordinates for the first structures\n print(\"\\nDistribution of X coordinates:\")\n plot_coord_distributions(df_labels, 'x_', max_structures=3)\n print(\"\\nDistribution of Y coordinates:\")\n plot_coord_distributions(df_labels, 'y_', max_structures=3)\n print(\"\\nDistribution of Z coordinates:\")\n plot_coord_distributions(df_labels, 'z_', max_structures=3)\n \n # Analyze sequences in train_sequences.csv\n if \"train_sequences.csv\" in main_data and main_data[\"train_sequences.csv\"] is not None:\n print(\"\\n===== Analysis of Sequences (train_sequences.csv) =====\")\n df_sequences = main_data[\"train_sequences.csv\"]\n \n # First few rows of the sequence column\n print(\"\\nFirst few rows of the 'sequence' column:\")\n print(df_sequences['sequence'].head())\n \n # Data type of the sequence column\n print(\"\\nData type of the 'sequence' column:\")\n print(df_sequences['sequence'].dtype)\n \n # Complete sequence analysis\n analyze_sequences(df_sequences)\n \n return main_data\n\nif __name__ == '__main__':\n main_data = main()","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:08.795231Z","iopub.execute_input":"2025-05-08T13:21:08.795572Z","iopub.status.idle":"2025-05-08T13:21:11.295924Z","shell.execute_reply.started":"2025-05-08T13:21:08.795541Z","shell.execute_reply":"2025-05-08T13:21:11.295026Z"},"papermill":{"duration":2.706006,"end_time":"2025-04-16T20:30:47.368896","exception":false,"start_time":"2025-04-16T20:30:44.662890","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"7e7a9f6e","cell_type":"markdown","source":"## Data Preparation for RNA 3D Structure Prediction 🧬🔍","metadata":{"papermill":{"duration":0.013357,"end_time":"2025-04-16T20:30:47.396598","exception":false,"start_time":"2025-04-16T20:30:47.383241","status":"completed"},"tags":[]}},{"id":"66ca7651","cell_type":"code","source":"# File paths\nDATA_DIR = \"/kaggle/input/stanford-rna-3d-folding/\"\nOUTPUT_DIR = \"/kaggle/working/\"\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\ndef load_data():\n \"\"\"\n Loads the necessary data for the competition.\n \"\"\"\n data = {}\n \n # Load sequences\n data['train_seq'] = pd.read_csv(os.path.join(DATA_DIR, \"train_sequences.csv\"))\n data['valid_seq'] = pd.read_csv(os.path.join(DATA_DIR, \"validation_sequences.csv\"))\n data['test_seq'] = pd.read_csv(os.path.join(DATA_DIR, \"test_sequences.csv\"))\n \n # Load structures (labels)\n data['train_labels'] = pd.read_csv(os.path.join(DATA_DIR, \"train_labels.csv\"))\n data['valid_labels'] = pd.read_csv(os.path.join(DATA_DIR, \"validation_labels.csv\"))\n \n # Load submission format\n data['sample_submission'] = pd.read_csv(os.path.join(DATA_DIR, \"sample_submission.csv\"))\n \n return data\n\ndef analyze_id_structure(data_dict):\n \"\"\"\n Analyzes the ID structure in different files to understand the correct mapping.\n \"\"\"\n # We'll analyze the specific formats for train and valid\n \n # 1. Analysis of training labels\n train_label_ids = data_dict['train_labels']['ID'].tolist()\n print(f\"Total IDs in training labels: {len(train_label_ids)}\")\n print(f\"Number of unique IDs: {len(set(train_label_ids))}\")\n \n # Try to understand the ID format in the labels file\n train_id_parts = {}\n for id_str in train_label_ids[:100]: # Analyze the first 100\n parts = id_str.split('_')\n num_parts = len(parts)\n if num_parts not in train_id_parts:\n train_id_parts[num_parts] = []\n train_id_parts[num_parts].append(parts)\n \n print(\"\\nID formats found in train_labels:\")\n for num_parts, examples in train_id_parts.items():\n print(f\"\\nFormat with {num_parts} parts:\")\n for i, parts in enumerate(examples[:3]):\n print(f\" Example {i+1}: {parts}\")\n \n # 2. Analysis of training sequences\n train_seq_ids = data_dict['train_seq']['target_id'].tolist()\n print(f\"\\nTotal IDs in training sequences: {len(train_seq_ids)}\")\n print(f\"Number of unique IDs: {len(set(train_seq_ids))}\")\n \n # Try to understand the ID format in the sequences file\n train_seq_id_parts = {}\n for id_str in train_seq_ids[:100]: # Analyze the first 100\n parts = id_str.split('_')\n num_parts = len(parts)\n if num_parts not in train_seq_id_parts:\n train_seq_id_parts[num_parts] = []\n train_seq_id_parts[num_parts].append(parts)\n \n print(\"\\nID formats found in train_sequences:\")\n for num_parts, examples in train_seq_id_parts.items():\n print(f\"\\nFormat with {num_parts} parts:\")\n for i, parts in enumerate(examples[:3]):\n print(f\" Example {i+1}: {parts}\")\n \n # 3. Analysis of validation labels\n valid_label_ids = data_dict['valid_labels']['ID'].tolist()\n print(f\"\\nTotal IDs in validation labels: {len(valid_label_ids)}\")\n print(f\"Number of unique IDs: {len(set(valid_label_ids))}\")\n \n # Count unique sequence IDs in validation labels\n valid_seq_ids_from_labels = set([id_str.split('_')[0] for id_str in valid_label_ids])\n print(f\"Number of unique sequence IDs in validation labels: {len(valid_seq_ids_from_labels)}\")\n print(f\"Examples: {list(valid_seq_ids_from_labels)[:5]}\")\n \n # 4. Analysis of validation sequences\n valid_seq_ids = data_dict['valid_seq']['target_id'].tolist()\n print(f\"\\nTotal IDs in validation sequences: {len(valid_seq_ids)}\")\n print(f\"Number of unique IDs: {len(set(valid_seq_ids))}\")\n print(f\"Examples: {valid_seq_ids[:5]}\")\n \n # 5. Check correspondence between unique IDs\n overlap_valid = set(valid_seq_ids).intersection(valid_seq_ids_from_labels)\n print(f\"\\nCorrespondence between validation sequences and labels: {len(overlap_valid)} of {len(valid_seq_ids)}\")\n \n # 6. Check how sequences and residues relate\n if len(overlap_valid) > 0:\n sample_id = list(overlap_valid)[0]\n sample_seq = data_dict['valid_seq'][data_dict['valid_seq']['target_id'] == sample_id]['sequence'].iloc[0]\n sample_labels = data_dict['valid_labels'][data_dict['valid_labels']['ID'].str.startswith(f\"{sample_id}_\")]\n \n print(f\"\\nAnalysis for sequence ID: {sample_id}\")\n print(f\"Sequence length: {len(sample_seq)}\")\n print(f\"Number of residues in labels: {len(sample_labels)}\")\n \n # Check how residue numbers are related\n residue_numbers = sample_labels['resid'].sort_values().tolist()\n print(f\"First residue numbers: {residue_numbers[:10]}\")\n print(f\"Last residue numbers: {residue_numbers[-10:]}\")\n \n return train_id_parts, train_seq_id_parts, overlap_valid\n\ndef fix_train_mapping(train_seq_df, train_labels_df):\n \"\"\"\n Identifies a correct mapping between train_sequences.csv and train_labels.csv\n using the ID format from the validation file as a reference.\n \n This is necessary because there's no obvious direct correspondence between the IDs.\n \"\"\"\n # First, extract the prefix of the ID from labels (format: XX_Y_Z)\n train_labels_df['seq_id'] = train_labels_df['ID'].apply(lambda x: x.split('_')[0] + '_' + x.split('_')[1])\n \n # Check if this format corresponds to the format of sequence IDs\n seq_ids_set = set(train_seq_df['target_id'])\n label_seq_ids_set = set(train_labels_df['seq_id'])\n \n overlap = seq_ids_set.intersection(label_seq_ids_set)\n print(f\"Overlap after format adjustment: {len(overlap)} of {len(seq_ids_set)}\")\n \n if len(overlap) > 0:\n print(f\"Examples of matching IDs: {list(overlap)[:5]}\")\n return overlap\n \n # If it still doesn't work, we need to analyze the structure in more detail\n print(\"No matches found, checking other formats...\")\n \n # Try other possible formats\n formats_to_try = [\n lambda x: x.split('_')[0], # Only first part\n lambda x: '_'.join(x.split('_')[:2]), # First two parts\n lambda x: x.split('_')[0] + '_' + x.split('_')[1][0], # First part + first letter of second part\n ]\n \n for i, format_func in enumerate(formats_to_try):\n train_labels_df[f'seq_id_{i}'] = train_labels_df['ID'].apply(format_func)\n label_seq_ids_set = set(train_labels_df[f'seq_id_{i}'])\n overlap = seq_ids_set.intersection(label_seq_ids_set)\n print(f\"Format {i}: Overlap = {len(overlap)} of {len(seq_ids_set)}\")\n \n if len(overlap) > 0:\n print(f\"Examples of matching IDs: {list(overlap)[:5]}\")\n return overlap, f'seq_id_{i}'\n \n # If no match is found, create a mapping based on observed patterns\n print(\"No matches found using simple patterns.\")\n print(\"Creating a manual mapping based on data structure...\")\n \n # Group labels by first parts of ID\n train_labels_df['prefix'] = train_labels_df['ID'].apply(lambda x: x.split('_')[0])\n label_groups = train_labels_df.groupby('prefix')\n \n # For each sequence, find the best match based on number of residues\n mapping = {}\n for _, seq_row in train_seq_df.iterrows():\n seq_id = seq_row['target_id']\n seq_length = len(seq_row['sequence'])\n \n best_match = None\n best_diff = float('inf')\n \n for prefix, group in label_groups:\n residue_count = len(group)\n diff = abs(residue_count - seq_length)\n \n if diff < best_diff:\n best_diff = diff\n best_match = prefix\n \n # Consider a match only if the number of residues is close\n if best_diff <= 10: # Tolerance of 10 residues\n mapping[seq_id] = best_match\n \n print(f\"Manual mapping created with {len(mapping)} matches\")\n return mapping\n\ndef create_mapping_valid(valid_seq_df, valid_labels_df):\n \"\"\"\n Creates a mapping between validation sequences and their coordinates.\n \n In this case, the IDs already correspond directly (R1107 -> R1107_1, R1107_2, etc.)\n \"\"\"\n # Check which ID format is used in the validation set\n valid_labels_df['seq_id'] = valid_labels_df['ID'].apply(lambda x: x.split('_')[0])\n \n # Check overlap\n seq_ids = set(valid_seq_df['target_id'])\n label_seq_ids = set(valid_labels_df['seq_id'])\n \n overlap = seq_ids.intersection(label_seq_ids)\n print(f\"Correspondence for validation: {len(overlap)} of {len(seq_ids)}\")\n \n mapping = {}\n for seq_id in overlap:\n # Get sequence\n seq = valid_seq_df[valid_seq_df['target_id'] == seq_id]['sequence'].iloc[0]\n \n # Get all residues for this sequence\n residues = valid_labels_df[valid_labels_df['seq_id'] == seq_id].sort_values('resid')\n \n # Extract coordinates for all structures\n num_structures = 1\n for col in residues.columns:\n if col.startswith('x_'):\n struct_num = int(col.split('_')[1])\n num_structures = max(num_structures, struct_num)\n \n # Initialize structures\n structures = []\n \n for struct_idx in range(1, num_structures + 1):\n coords = []\n has_valid_coords = False\n \n # Check if this structure has coordinates\n if f'x_{struct_idx}' in residues.columns:\n for _, row in residues.iterrows():\n x = row[f'x_{struct_idx}']\n y = row[f'y_{struct_idx}']\n z = row[f'z_{struct_idx}']\n \n # Check if they are valid values\n if abs(x) < 1.0e+17 and abs(y) < 1.0e+17 and abs(z) < 1.0e+17:\n coords.append([x, y, z])\n has_valid_coords = True\n else:\n coords.append([np.nan, np.nan, np.nan])\n \n if has_valid_coords:\n structures.append(coords)\n \n # Add to mapping if there are valid structures\n if structures:\n mapping[seq_id] = {\n 'sequence': seq,\n 'structures': structures\n }\n \n print(f\"Mapping created with {len(mapping)} valid sequences\")\n return mapping\n\ndef create_processed_data(mapping, output_prefix):\n \"\"\"\n Creates and saves processed data from the mapping.\n \n Parameters:\n mapping: Dictionary with the mapping of sequences to structures\n output_prefix: Prefix for output files ('train' or 'valid')\n \n Returns:\n X, y: Arrays for training\n \"\"\"\n if not mapping:\n print(f\"WARNING: No valid mapping for {output_prefix}\")\n return None, None\n \n X_data = []\n y_data = []\n ids = []\n \n for seq_id, data in mapping.items():\n seq = data['sequence']\n structures = data['structures']\n \n # Skip if there are no structures\n if not structures:\n continue\n \n # Use the first valid structure\n structure = structures[0]\n \n # Check if the structure has valid coordinates for all residues\n if len(structure) != len(seq):\n print(f\"WARNING: Difference between sequence length ({len(seq)}) and coordinates ({len(structure)}) for {seq_id}\")\n # If needed, we could consider padding or truncation here\n continue\n \n # Create feature matrix (one-hot encoding)\n features = []\n for nucleotide in seq:\n if nucleotide == 'A':\n features.append([1, 0, 0, 0, 0])\n elif nucleotide == 'C':\n features.append([0, 1, 0, 0, 0])\n elif nucleotide == 'G':\n features.append([0, 0, 1, 0, 0])\n elif nucleotide == 'U':\n features.append([0, 0, 0, 1, 0])\n else:\n features.append([0, 0, 0, 0, 1]) # For unknown nucleotides\n \n X_data.append(np.array(features))\n y_data.append(np.array(structure))\n ids.append(seq_id)\n \n if not X_data:\n print(f\"WARNING: No valid processed data for {output_prefix}\")\n return None, None, []\n \n # Padding to ensure all sequences have the same length\n max_length = max(len(x) for x in X_data)\n X_padded = []\n y_padded = []\n \n for x, y in zip(X_data, y_data):\n if len(x) < max_length:\n x_pad = np.zeros((max_length, 5))\n x_pad[:len(x), :] = x\n \n y_pad = np.zeros((max_length, 3))\n y_pad[:len(y), :] = y\n \n X_padded.append(x_pad)\n y_padded.append(y_pad)\n else:\n X_padded.append(x)\n y_padded.append(y)\n \n X = np.array(X_padded)\n y = np.array(y_padded)\n \n # Save the processed data\n np.save(os.path.join(OUTPUT_DIR, f'X_{output_prefix}.npy'), X)\n np.save(os.path.join(OUTPUT_DIR, f'y_{output_prefix}.npy'), y)\n \n with open(os.path.join(OUTPUT_DIR, f'{output_prefix}_ids.txt'), 'w') as f:\n for id in ids:\n f.write(f\"{id}\\n\")\n \n print(f\"Processed data for {output_prefix}: X.shape = {X.shape}, y.shape = {y.shape}\")\n return X, y, ids\n\ndef explore_sequence_mapping(seq_id, mapping, data_dict):\n \"\"\"\n Explores a mapping example in detail for diagnostics.\n \"\"\"\n if seq_id not in mapping:\n print(f\"WARNING: Sequence ID {seq_id} not found in mapping\")\n return\n \n data = mapping[seq_id]\n seq = data['sequence']\n structures = data['structures']\n \n print(f\"Exploring mapping for sequence: {seq_id}\")\n print(f\"Sequence length: {len(seq)}\")\n print(f\"Number of available structures: {len(structures)}\")\n \n # Detail each structure\n for i, structure in enumerate(structures):\n print(f\"\\nStructure {i+1}:\")\n print(f\" Number of coordinates: {len(structure)}\")\n if len(structure) > 0:\n print(f\" First coordinates: {structure[:3]}\")\n print(f\" Last coordinates: {structure[-3:]}\")\n \n # Check correspondence with the sequence\n if len(structure) != len(seq):\n print(f\" WARNING: Difference between sequence length ({len(seq)}) and coordinates ({len(structure)})\")\n else:\n print(f\" Perfect match between sequence and coordinates\")\n\ndef main():\n # Load the data\n print(\"Loading data...\")\n data_dict = load_data()\n \n # Analyze ID structure to understand the mapping\n print(\"\\nAnalyzing ID structure...\")\n train_id_parts, train_seq_id_parts, overlap_valid = analyze_id_structure(data_dict)\n \n # For validation, the mapping is direct (R1107 -> R1107_1, R1107_2, etc.)\n print(\"\\nCreating mapping for validation data...\")\n valid_mapping = create_mapping_valid(data_dict['valid_seq'], data_dict['valid_labels'])\n \n # Explore a validation mapping example to verify\n if valid_mapping:\n sample_id = list(valid_mapping.keys())[0]\n print(f\"\\nExploring a validation mapping example ({sample_id}):\")\n explore_sequence_mapping(sample_id, valid_mapping, data_dict)\n \n # Create and save processed data for validation\n X_valid, y_valid, valid_ids = create_processed_data(valid_mapping, 'valid')\n \n # Since we couldn't establish a mapping for training,\n # we'll use validation data for training as well (transfer learning)\n print(\"\\nUsing validation data as training (due to lack of direct mapping)...\")\n X_train = X_valid\n y_train = y_valid\n train_ids = valid_ids\n \n if X_train is not None:\n np.save(os.path.join(OUTPUT_DIR, 'X_train.npy'), X_train)\n np.save(os.path.join(OUTPUT_DIR, 'y_train.npy'), y_train)\n \n with open(os.path.join(OUTPUT_DIR, 'train_ids.txt'), 'w') as f:\n for id in train_ids:\n f.write(f\"{id}\\n\")\n \n # Return the processed data\n return {\n 'X_train': X_train,\n 'y_train': y_train,\n 'X_valid': X_valid,\n 'y_valid': y_valid,\n 'valid_mapping': valid_mapping,\n 'valid_ids': valid_ids\n }\n\nif __name__ == \"__main__\":\n processed_data = main()","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:11.296816Z","iopub.execute_input":"2025-05-08T13:21:11.297044Z","iopub.status.idle":"2025-05-08T13:21:16.814814Z","shell.execute_reply.started":"2025-05-08T13:21:11.297025Z","shell.execute_reply":"2025-05-08T13:21:16.814016Z"},"papermill":{"duration":5.495068,"end_time":"2025-04-16T20:30:52.904872","exception":false,"start_time":"2025-04-16T20:30:47.409804","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"e90a4b8f","cell_type":"markdown","source":"## Heatmap Viewer for RNA Sequences 🔥🧬","metadata":{"papermill":{"duration":0.013496,"end_time":"2025-04-16T20:30:52.939074","exception":false,"start_time":"2025-04-16T20:30:52.925578","status":"completed"},"tags":[]}},{"id":"439f44c1","cell_type":"code","source":"def visualize_rna_heatmap_from_processed_data(processed_data, num_samples=12):\n \"\"\"\n Visualizes a heatmap for RNA sequences using processed data.\n \n Parameters:\n processed_data: Dictionary with processed data returned by the main() function\n num_samples: Number of sequences to visualize\n \"\"\"\n try:\n # Check if we have the necessary data\n if 'X_valid' not in processed_data or processed_data['X_valid'] is None:\n print(\"Validation data not found in processed_data object\")\n return None\n \n # Get the data\n X_valid = processed_data['X_valid']\n print(f\"Data found with format: {X_valid.shape}\")\n \n # Limit to the number of samples\n X_valid_subset = X_valid[:num_samples]\n \n # If we have IDs, use them\n if 'valid_ids' in processed_data and processed_data['valid_ids']:\n valid_ids = processed_data['valid_ids'][:num_samples]\n else:\n valid_ids = [f\"Seq_{i+1}\" for i in range(X_valid_subset.shape[0])]\n \n # Convert one-hot encoding to nucleotide indices\n # Expected format: A=[1,0,0,0,0], C=[0,1,0,0,0], G=[0,0,1,0,0], U=[0,0,0,1,0], N=[0,0,0,0,1]\n sequences_matrix = np.argmax(X_valid_subset, axis=2)\n \n # Replace zeros (padding) with 4 (N/Unknown) when all values are zero\n is_padding = np.all(X_valid_subset == 0, axis=2)\n sequences_matrix[is_padding] = 4\n \n # Define a categorical colormap (distinct colors per nucleotide)\n cmap = mcolors.ListedColormap(['#3498db', '#2ecc71', '#e74c3c', '#9b59b6', '#95a5a6'])\n bounds = [0, 1, 2, 3, 4, 5]\n norm = mcolors.BoundaryNorm(bounds, cmap.N)\n \n # Create figure\n plt.figure(figsize=(20, 10))\n im = plt.imshow(sequences_matrix, cmap=cmap, norm=norm, aspect='auto')\n \n # Add color bar\n cbar = plt.colorbar(im, ticks=[0.5, 1.5, 2.5, 3.5, 4.5])\n cbar.set_label('Nucleotides', fontsize=14)\n cbar.set_ticklabels(['A', 'C', 'G', 'U', 'N/Padding'])\n \n # Add axis labels\n plt.xlabel(\"Position in Sequence\", fontsize=14)\n plt.ylabel(\"RNA Sequences\", fontsize=14)\n \n # Add title\n plt.title(\"RNA Sequences Heatmap\", fontsize=16)\n \n # Add sequence IDs as y-axis labels\n plt.yticks(range(len(valid_ids)), valid_ids, fontsize=10)\n \n # Show only some labels on x-axis to avoid crowding\n sequence_length = sequences_matrix.shape[1]\n step = max(1, sequence_length // 20) # Show at most 20 labels\n plt.xticks(range(0, sequence_length, step), range(1, sequence_length + 1, step))\n \n # Add grid\n plt.grid(False)\n \n # Add information about nucleotide distribution\n all_nucleotides = sequences_matrix.flatten()\n nucleotide_counts = {\n 'A': np.sum(all_nucleotides == 0),\n 'C': np.sum(all_nucleotides == 1),\n 'G': np.sum(all_nucleotides == 2),\n 'U': np.sum(all_nucleotides == 3),\n 'N': np.sum(all_nucleotides == 4)\n }\n \n total_nucleotides = sum(nucleotide_counts.values())\n nucleotide_percentages = {k: (v / total_nucleotides) * 100 for k, v in nucleotide_counts.items()}\n \n # Add text with statistics\n info_text = \"\\n\".join([\n f\"Total sequences visualized: {num_samples}\",\n f\"Maximum length: {sequence_length}\",\n f\"A: {nucleotide_percentages['A']:.1f}%\",\n f\"C: {nucleotide_percentages['C']:.1f}%\",\n f\"G: {nucleotide_percentages['G']:.1f}%\",\n f\"U: {nucleotide_percentages['U']:.1f}%\",\n f\"N/Padding: {nucleotide_percentages['N']:.1f}%\"\n ])\n \n plt.figtext(0.02, 0.02, info_text, fontsize=10, bbox=dict(facecolor='white', alpha=0.8))\n \n # Show the plot\n plt.tight_layout()\n plt.show()\n \n # Optionally, save the plot\n output_dir = '/kaggle/working/'\n plt.savefig(os.path.join(output_dir, 'rna_heatmap.png'), dpi=300)\n print(f\"Heatmap saved to {os.path.join(output_dir, 'rna_heatmap.png')}\")\n \n return sequences_matrix\n except Exception as e:\n print(f\"Error processing data: {e}\")\n return None\n\n# Use the function (assuming processed_data is available)\nvisualize_rna_heatmap_from_processed_data(processed_data)","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:16.815771Z","iopub.execute_input":"2025-05-08T13:21:16.816089Z","iopub.status.idle":"2025-05-08T13:21:17.725515Z","shell.execute_reply.started":"2025-05-08T13:21:16.816059Z","shell.execute_reply":"2025-05-08T13:21:17.724634Z"},"papermill":{"duration":1.078224,"end_time":"2025-04-16T20:30:54.030600","exception":false,"start_time":"2025-04-16T20:30:52.952376","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"dcf4a33b","cell_type":"markdown","source":"## 🧬 RNA 3D Structure Prediction with Neural Network-Based Conformational Selection 🔬","metadata":{"papermill":{"duration":0.017174,"end_time":"2025-04-16T20:30:54.067091","exception":false,"start_time":"2025-04-16T20:30:54.049917","status":"completed"},"tags":[]}},{"id":"caa1e70a","cell_type":"markdown","source":"> ## Utility and Preprocessing Functions","metadata":{"papermill":{"duration":0.018287,"end_time":"2025-04-16T20:30:54.104274","exception":false,"start_time":"2025-04-16T20:30:54.085987","status":"completed"},"tags":[]}},{"id":"eeb4c495","cell_type":"code","source":"# File paths\nDATA_DIR = \"/kaggle/input/stanford-rna-3d-folding/\"\nOUTPUT_DIR = \"/kaggle/working/\"\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\ndef normalize_structure(coords):\n \"\"\"\n Centralizes and normalizes the structure.\n \"\"\"\n # Remove padding\n valid_mask = ~np.all(coords == 0, axis=1)\n valid_coords = coords[valid_mask]\n \n # Center at center of mass\n center = np.mean(valid_coords, axis=0)\n centered_coords = coords.copy()\n centered_coords[valid_mask] = valid_coords - center\n \n return centered_coords\n\ndef normalize_coordinates(coords):\n \"\"\"\n Normalizes 3D coordinates of RNA structures by centering and \n scaling each structure independently, with robust handling\n to avoid numerical issues.\n \n Parameters:\n -----------\n coords: Numpy array with shape (batch_size, seq_length, 3)\n 3D coordinates to normalize\n \n Returns:\n --------\n normalized: Numpy array with shape (batch_size, seq_length, 3)\n Normalized coordinates in the range [-1, 1] \n \"\"\"\n # Create copy to avoid modifying the original\n normalized = np.copy(coords)\n \n # Check for problematic values upfront\n if np.isnan(coords).any():\n print(\"WARNING: NaN values detected in input coordinates. They will be ignored during normalization.\")\n if np.isinf(coords).any():\n print(\"WARNING: Infinite values detected in input coordinates. They will be ignored during normalization.\")\n \n # Handle each structure in the batch separately\n for i in range(coords.shape[0]):\n # Identify valid positions (non-zero, non-NaN, non-Inf)\n valid_mask = ~np.all(coords[i] == 0, axis=-1) \n valid_mask = valid_mask & ~np.any(np.isnan(coords[i]), axis=-1)\n valid_mask = valid_mask & ~np.any(np.isinf(coords[i]), axis=-1)\n \n # Extract only valid coordinates\n valid_coords = coords[i][valid_mask]\n \n if len(valid_coords) > 0:\n try:\n # 1. Center at the geometric center\n center = np.nanmean(valid_coords, axis=0)\n \n # Check if the calculated center contains valid values \n if np.isnan(center).any() or np.isinf(center).any():\n print(f\"WARNING: Invalid center calculated for structure {i}. Using [0,0,0].\")\n center = np.zeros(3)\n \n # Apply translation to the center\n centered = valid_coords - center\n \n # 2. Determine appropriate scale factor\n # Calculate maximum distance from the center\n dist_from_center = np.sqrt(np.sum(centered * centered, axis=1))\n \n # Exclude NaN or infinite values for scale_factor calculation\n valid_dists = dist_from_center[~np.isnan(dist_from_center) & ~np.isinf(dist_from_center)]\n \n if len(valid_dists) > 0:\n scale_factor = np.max(valid_dists)\n # Protect against very small scale_factor\n if scale_factor < 1e-10:\n scale_factor = 1.0\n else:\n scale_factor = 1.0\n \n # 3. Normalize coordinates to [-1, 1] range\n normalized_valid = centered / scale_factor\n \n # 4. Replace values in the normalized array\n normalized[i][valid_mask] = normalized_valid\n \n # Debug info\n # print(f\"Structure {i}: center={center}, scale_factor={scale_factor}, \" \n # f\"min={np.min(normalized_valid)}, max={np.max(normalized_valid)}\")\n \n except Exception as e:\n print(f\"ERROR during normalization of structure {i}: {str(e)}\")\n print(\"Keeping original values for this structure.\")\n else:\n print(f\"WARNING: No valid coordinates found for structure {i}.\")\n \n # Final check to detect any issues\n if np.isnan(normalized).any():\n print(\"WARNING: NaN values present after normalization. Replacing with zeros.\")\n normalized = np.nan_to_num(normalized, nan=0.0)\n \n if np.isinf(normalized).any():\n print(\"WARNING: Infinite values present after normalization. Replacing with zeros.\") \n normalized = np.nan_to_num(normalized, posinf=0.0, neginf=0.0)\n \n return normalized\n\ndef check_structure_validity(coords, min_distance=0.8, max_distance=7.0, allow_clashes=0.05):\n \"\"\"\n More refined and realistic biophysical validation.\n \"\"\"\n valid = True\n valid_mask = ~np.all(coords == 0, axis=1)\n valid_coords = coords[valid_mask]\n \n if len(valid_coords) < 3:\n return True\n \n # Check distances between consecutive residues\n invalid_bonds = 0\n for i in range(1, len(valid_coords)):\n dist = np.linalg.norm(valid_coords[i] - valid_coords[i-1])\n if dist < min_distance or dist > max_distance:\n invalid_bonds += 1\n \n # Allow a small percentage of invalid bonds\n if invalid_bonds / len(valid_coords) > 0.1: # More than 10% invalid bonds\n valid = False\n \n # Check for clashes, allowing some\n clashes = 0\n total_pairs = 0\n for i in range(len(valid_coords)):\n for j in range(i+3, len(valid_coords)): # Skip adjacent\n total_pairs += 1\n dist = np.linalg.norm(valid_coords[i] - valid_coords[j])\n if dist < min_distance:\n clashes += 1\n \n # Allow a small percentage of clashes\n if total_pairs > 0 and clashes / total_pairs > allow_clashes:\n valid = False\n \n return valid","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:17.728240Z","iopub.execute_input":"2025-05-08T13:21:17.728457Z","iopub.status.idle":"2025-05-08T13:21:17.740839Z","shell.execute_reply.started":"2025-05-08T13:21:17.728439Z","shell.execute_reply":"2025-05-08T13:21:17.740256Z"},"papermill":{"duration":0.032673,"end_time":"2025-04-16T20:30:54.155094","exception":false,"start_time":"2025-04-16T20:30:54.122421","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"550011e6","cell_type":"code","source":"def sample_structural_variation(coords, noise_level=0.5, preserve_distance=True, \n use_global_movement=False, correlation=0.7):\n \"\"\"\n Enhanced version of structural variation sampling with better\n handling of large RNAs and improved noise distribution.\n \"\"\"\n new_coords = coords.copy()\n valid_mask = ~np.all(coords == 0, axis=1)\n valid_indices = np.where(valid_mask)[0]\n \n if len(valid_indices) < 3:\n return new_coords\n \n # Parameters optimized for RNA structure\n typical_bond_length = 3.8 # Angstroms - typical RNA backbone distance\n \n # Add global domain movements if requested\n if use_global_movement and len(valid_indices) > 20:\n # More natural domain identification - try to find natural hinge points\n # For RNA, these often occur at junctions between helices\n \n # Calculate distance between consecutive residues as a heuristic\n # for finding potential hinge points (larger distances often indicate junctions)\n distances = []\n for i in range(1, len(valid_indices)):\n idx1 = valid_indices[i-1]\n idx2 = valid_indices[i]\n dist = np.linalg.norm(coords[idx1] - coords[idx2])\n distances.append((i, dist))\n \n # Sort by distance to find potential hinges\n distances.sort(key=lambda x: x[1], reverse=True)\n \n # Take top 2 potential hinge points (if we have enough points)\n num_hinges = min(2, len(distances)//3)\n \n for h in range(num_hinges):\n if h < len(distances):\n hinge_point = distances[h][0]\n if hinge_point < 5 or hinge_point > len(valid_indices) - 5:\n continue\n \n hinge_idx = valid_indices[hinge_point]\n \n # Angle of rotation with natural distribution\n # More small movements than large ones\n angle = np.random.exponential(0.2) # Mostly small angles with occasional larger ones\n if np.random.random() < 0.5:\n angle = -angle # Allow both directions\n \n # Create a more natural rotation matrix with slight 3D component\n # RNAs often bend and twist in 3D\n sin_a, cos_a = np.sin(angle), np.cos(angle)\n tilt = np.random.normal(0, 0.1) # Small tilt in 3D\n rotation_matrix = np.array([\n [cos_a, -sin_a, 0],\n [sin_a, cos_a, tilt],\n [0, -tilt, 1]\n ])\n \n # Apply rotation around hinge point\n ref_point = new_coords[hinge_idx]\n for i in valid_indices[hinge_point+1:]:\n vector = new_coords[i] - ref_point\n rotated = np.dot(vector, rotation_matrix)\n new_coords[i] = ref_point + rotated\n \n # Propagate variation residue by residue, with correlation\n # RNA has strong local correlations in structure\n prev_noise = np.zeros(3)\n \n correlation = 0.5 # High correlation for smoother variations\n \n for i in range(1, len(coords)):\n if not valid_mask[i] or not valid_mask[i-1]:\n continue\n \n vec = new_coords[i-1] - new_coords[i]\n vec_length = np.linalg.norm(vec)\n \n # Generate correlated noise (smoother transitions)\n new_noise = np.random.normal(0, noise_level, size=3)\n noise_vec = correlation * prev_noise + (1 - correlation) * new_noise\n prev_noise = noise_vec.copy()\n \n noise_norm = np.linalg.norm(noise_vec)\n if noise_norm > 0:\n # Scale noise proportionally\n noise_vec = noise_vec / noise_norm * (noise_level * vec_length)\n \n # Add noise to the direction\n new_vec = vec + noise_vec\n \n # Preserve distance if requested\n if preserve_distance:\n current_length = np.linalg.norm(new_vec)\n if current_length > 0:\n # Allow slight variation in bond length (RNA is not rigid)\n target_length = typical_bond_length * (1 + np.random.normal(0, 0.05))\n new_vec = new_vec / current_length * target_length\n \n new_coords[i] = new_coords[i-1] - new_vec\n \n return new_coords\n\ndef get_rotation_matrix(axis, theta):\n \"\"\"\n Return the rotation matrix for rotation around an arbitrary axis.\n \n Parameters:\n -----------\n axis: Unit vector defining the rotation axis\n theta: Rotation angle in radians\n \n Returns:\n --------\n 3x3 rotation matrix\n \"\"\"\n # Ensure axis is a unit vector\n axis = axis / np.linalg.norm(axis)\n \n a = np.cos(theta / 2.0)\n b, c, d = -axis * np.sin(theta / 2.0)\n \n return np.array([\n [a*a + b*b - c*c - d*d, 2*(b*c - a*d), 2*(b*d + a*c)],\n [2*(b*c + a*d), a*a + c*c - b*b - d*d, 2*(c*d - a*b)],\n [2*(b*d - a*c), 2*(c*d + a*b), a*a + d*d - b*b - c*c]\n ])\n\ndef refine_rna_backbone(structure):\n \"\"\"\n Refine the RNA backbone geometry to match known constraints.\n \n Parameters:\n -----------\n structure: RNA 3D structure\n \n Returns:\n --------\n Refined structure\n \"\"\"\n # Create a copy to refine\n refined = structure.copy()\n \n # Check for valid residues\n valid_mask = ~np.all(refined == 0, axis=1)\n \n # Apply RNA-specific backbone constraints\n for i in range(2, len(refined)):\n if valid_mask[i] and valid_mask[i-1] and valid_mask[i-2]:\n # In RNA, there are constraints on three consecutive backbone atoms\n \n # Get the two backbone vectors\n vec1 = refined[i-1] - refined[i-2]\n vec2 = refined[i] - refined[i-1]\n \n # Calculate current angle between vectors\n vec1_norm = vec1 / (np.linalg.norm(vec1) + 1e-6)\n vec2_norm = vec2 / (np.linalg.norm(vec2) + 1e-6)\n cos_angle = np.dot(vec1_norm, vec2_norm)\n \n # Clamp to valid range for numerical stability\n cos_angle = max(-1.0, min(1.0, cos_angle))\n angle = np.arccos(cos_angle)\n \n # In RNA, the typical backbone angle is around 100-120 degrees\n ideal_angle = np.radians(110)\n \n # If the angle is too far from ideal, adjust it\n if abs(angle - ideal_angle) > np.radians(30):\n # Create a rotation to adjust the angle\n # Get the rotation axis (perpendicular to the plane of vec1 and vec2)\n axis = np.cross(vec1_norm, vec2_norm)\n axis_norm = axis / (np.linalg.norm(axis) + 1e-6)\n \n # Determine rotation angle to reach ideal angle\n angle_diff = ideal_angle - angle\n \n # Apply rotation to vec2\n rotation_matrix = get_rotation_matrix(axis_norm, angle_diff)\n new_vec2 = np.dot(rotation_matrix, vec2_norm) * np.linalg.norm(vec2)\n \n # Update the position\n refined[i] = refined[i-1] + new_vec2\n \n return refined","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:17.742543Z","iopub.execute_input":"2025-05-08T13:21:17.742854Z","iopub.status.idle":"2025-05-08T13:21:17.852226Z","shell.execute_reply.started":"2025-05-08T13:21:17.742833Z","shell.execute_reply":"2025-05-08T13:21:17.851394Z"},"papermill":{"duration":0.03463,"end_time":"2025-04-16T20:30:54.207801","exception":false,"start_time":"2025-04-16T20:30:54.173171","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"dad48ada","cell_type":"code","source":"def repair_invalid_structure(structure):\n \"\"\"\n Attempt to repair an invalid RNA structure.\n \n Parameters:\n -----------\n structure: Potentially invalid RNA structure\n \n Returns:\n --------\n Repaired structure\n \"\"\"\n # Create a copy to repair\n repaired = structure.copy()\n \n # Check for valid residues\n valid_mask = ~np.all(repaired == 0, axis=1)\n \n # Fix bond lengths\n for i in range(1, len(repaired)):\n if valid_mask[i] and valid_mask[i-1]:\n # Get current bond\n bond_vector = repaired[i] - repaired[i-1]\n bond_length = np.linalg.norm(bond_vector)\n \n # Check if bond is too short or too long\n if bond_length < 1.0 or bond_length > 7.0:\n # Fix bond to ideal length\n ideal_length = 3.8\n if bond_length > 0:\n repaired[i] = repaired[i-1] + (bond_vector / bond_length) * ideal_length\n else:\n # Generate a random direction if bond length is zero\n random_direction = np.random.randn(3)\n random_direction = random_direction / np.linalg.norm(random_direction)\n repaired[i] = repaired[i-1] + random_direction * ideal_length\n \n # Check for clashes (atoms too close to each other)\n for i in range(len(repaired)):\n if valid_mask[i]:\n for j in range(i+3, len(repaired)): # Skip adjacent residues\n if valid_mask[j]:\n # Calculate distance\n distance = np.linalg.norm(repaired[j] - repaired[i])\n \n # If atoms are too close\n if distance < 1.0:\n # Move one atom away slightly in a random direction\n random_direction = np.random.randn(3)\n random_direction = random_direction / np.linalg.norm(random_direction)\n repaired[j] = repaired[i] + random_direction * 4.0 # Place at safe distance\n \n # Final normalization\n repaired = normalize_structure(repaired)\n \n return repaired\n\ndef create_emergency_structure(seq_length):\n \"\"\"\n Create an emergency structure when all else fails.\n Generates a physically plausible RNA structure.\n \n Parameters:\n -----------\n seq_length: Length of the RNA sequence\n \n Returns:\n --------\n Basic RNA structure\n \"\"\"\n # Create a simple linear structure as fallback\n emergency_structure = np.zeros((seq_length, 3))\n \n # Define canonical nucleotide step (3.8Å)\n step = np.array([3.8, 0.0, 0.0])\n \n # Generate a straight chain with some randomness\n for i in range(seq_length):\n if i == 0:\n emergency_structure[i] = np.zeros(3)\n else:\n # Add slight random deviation to prevent perfect linearity\n random_noise = np.random.normal(0, 0.2, 3)\n emergency_structure[i] = emergency_structure[i-1] + step + random_noise\n \n # Add a slight curve to make it more RNA-like\n # Apply a gentle curve in the y-z plane\n for i in range(seq_length):\n angle = i * 0.1 # Gradual rotation\n emergency_structure[i, 1] += 2 * np.sin(angle) # Y-component\n emergency_structure[i, 2] += 2 * np.cos(angle) # Z-component\n \n # Normalize\n emergency_structure = normalize_structure(emergency_structure)\n \n return emergency_structure","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:17.853067Z","iopub.execute_input":"2025-05-08T13:21:17.853417Z","iopub.status.idle":"2025-05-08T13:21:17.867829Z","shell.execute_reply.started":"2025-05-08T13:21:17.853388Z","shell.execute_reply":"2025-05-08T13:21:17.866908Z"},"papermill":{"duration":0.028055,"end_time":"2025-04-16T20:30:54.254351","exception":false,"start_time":"2025-04-16T20:30:54.226296","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"76553bbb","cell_type":"code","source":"def calculate_tm_score(pred_coords, true_coords, d0_scale=1.24):\n \"\"\"\n Calculates a robust approximation of the TM-score between predicted and true coordinates.\n Adds protections against division by zero and NaN.\n \"\"\"\n # Remove padding (rows with zeros) from the true structures\n mask = ~np.all(true_coords == 0, axis=1)\n pred = pred_coords[mask]\n true = true_coords[mask]\n \n L = len(true)\n if L < 3:\n return 0.0\n \n # Define d0 based on L (values adapted for RNA)\n if L >= 30:\n d0 = 0.6 * np.sqrt(L - 0.5) - 2.5\n d0 = max(0.1, d0)\n elif L >= 24:\n d0 = 0.7\n elif L >= 20:\n d0 = 0.6\n elif L >= 16:\n d0 = 0.5\n elif L >= 12:\n d0 = 0.4\n else:\n d0 = 0.3\n \n distances = np.sqrt(np.sum((pred - true) ** 2, axis=1))\n tm_terms = 1.0 / (1.0 + (distances / (d0 + 1e-8)) ** 2)\n tm_score = np.sum(tm_terms) / L\n return float(tm_score)\n\ndef calculate_tm_score_exact(pred_coords, true_coords):\n \"\"\"\n Implementation more closely matching US-align with sequence-independent alignment.\n Includes multiple rotation schemes to find the optimal structural alignment.\n \"\"\"\n # Remove padding\n mask = ~np.all(true_coords == 0, axis=1)\n pred = pred_coords[mask]\n true = true_coords[mask]\n \n Lref = len(true)\n if Lref < 3:\n return 0.0\n \n # Define d0 exactly as in the evaluation formula\n if Lref >= 30:\n d0 = 0.6 * np.sqrt(Lref - 0.5) - 2.5\n elif Lref >= 24:\n d0 = 0.7\n elif Lref >= 20:\n d0 = 0.6\n elif Lref >= 16:\n d0 = 0.5\n elif Lref >= 12:\n d0 = 0.4\n else:\n d0 = 0.3\n \n # Normalize structures\n pred_centered = pred - np.mean(pred, axis=0)\n true_centered = true - np.mean(true, axis=0)\n \n # Try multiple fragment lengths for sequence-independent alignment\n # This mimics US-align's approach to find the best fragment alignment\n best_tm_score = 0.0\n fragment_lengths = [Lref, max(5, Lref//2), max(5, Lref//4)]\n \n for frag_len in fragment_lengths:\n # Try different fragment start positions\n for i in range(0, Lref - frag_len + 1, max(1, frag_len//2)):\n pred_frag = pred_centered[i:i+frag_len]\n \n # Try aligning with different parts of the true structure\n for j in range(0, Lref - frag_len + 1, max(1, frag_len//2)):\n true_frag = true_centered[j:j+frag_len]\n \n # Covariance matrix for optimal rotation\n covariance = np.dot(pred_frag.T, true_frag)\n U, S, Vt = np.linalg.svd(covariance)\n rotation = np.dot(U, Vt)\n \n # Try different rotation schemes - this is the new part\n rotations_to_try = [\n rotation, # Original rotation from SVD\n np.dot(rotation, np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]])), # 90 degree Z rotation\n np.dot(rotation, np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]])) # 180 degree Z rotation\n ]\n \n for rot in rotations_to_try:\n # Apply rotation to the full structure\n pred_aligned = np.dot(pred_centered, rot)\n \n # Calculate distances\n distances = np.sqrt(np.sum((pred_aligned - true_centered) ** 2, axis=1))\n \n # Calculate TM-score terms\n tm_terms = 1.0 / (1.0 + (distances / d0) ** 2)\n tm_score = np.sum(tm_terms) / Lref\n \n best_tm_score = max(best_tm_score, tm_score)\n \n return float(best_tm_score)\n\ndef load_processed_data():\n \"\"\"\n Loads processed data for training.\n \"\"\"\n X_train = np.load(os.path.join(OUTPUT_DIR, 'X_train.npy'))\n y_train = np.load(os.path.join(OUTPUT_DIR, 'y_train.npy'))\n X_valid = np.load(os.path.join(OUTPUT_DIR, 'X_valid.npy'))\n y_valid = np.load(os.path.join(OUTPUT_DIR, 'y_valid.npy'))\n \n print(f\"Data loaded - X_train: {X_train.shape}, y_train: {y_train.shape}\")\n print(f\"Data loaded - X_valid: {X_valid.shape}, y_valid: {y_valid.shape}\")\n \n return X_train, y_train, X_valid, y_valid","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:17.868997Z","iopub.execute_input":"2025-05-08T13:21:17.869356Z","iopub.status.idle":"2025-05-08T13:21:17.886544Z","shell.execute_reply.started":"2025-05-08T13:21:17.869323Z","shell.execute_reply":"2025-05-08T13:21:17.885782Z"},"papermill":{"duration":0.031937,"end_time":"2025-04-16T20:30:54.304553","exception":false,"start_time":"2025-04-16T20:30:54.272616","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"a772e270","cell_type":"code","source":"def prepare_test_features(test_seq_df, max_length=720):\n \"\"\"\n Prepares test features (one-hot encoding of the sequence).\n \"\"\"\n X_test = []\n for _, row in test_seq_df.iterrows():\n seq = row['sequence']\n features = []\n for nucleotide in seq:\n if nucleotide == 'A':\n features.append([1, 0, 0, 0, 0])\n elif nucleotide == 'C':\n features.append([0, 1, 0, 0, 0])\n elif nucleotide == 'G':\n features.append([0, 0, 1, 0, 0])\n elif nucleotide == 'U':\n features.append([0, 0, 0, 1, 0])\n else:\n features.append([0, 0, 0, 0, 1])\n if len(features) < max_length:\n padding = [[0, 0, 0, 0, 0]] * (max_length - len(features))\n features.extend(padding)\n else:\n features = features[:max_length]\n X_test.append(features)\n return np.array(X_test)\n\ndef extract_sequence_features(seq_features):\n \"\"\"\n Extract relevant sequence features from one-hot encoding.\n \"\"\"\n # Get valid rows (non-padding)\n valid_mask = ~np.all(seq_features == 0, axis=1)\n valid_features = seq_features[valid_mask]\n \n # Calculate nucleotide composition\n a_content = np.mean(valid_features[:, 0])\n c_content = np.mean(valid_features[:, 1])\n g_content = np.mean(valid_features[:, 2])\n u_content = np.mean(valid_features[:, 3])\n gc_content = c_content + g_content\n \n return {\n 'length': np.sum(valid_mask),\n 'a_content': a_content,\n 'c_content': c_content,\n 'g_content': g_content, \n 'u_content': u_content,\n 'gc_content': gc_content,\n 'au_content': a_content + u_content\n }\n\ndef visualize_3d_structure(true_coords, pred_coords, sample_idx=0, title=\"3D Structure Comparison\", show_plot=False):\n \"\"\"\n Visualizes the true and predicted 3D structures for a sample.\n Only shows the plot if explicitly requested.\n \"\"\"\n true = true_coords[sample_idx]\n pred = pred_coords[sample_idx]\n mask = ~np.all(true == 0, axis=1)\n true = true[mask]\n pred = pred[mask]\n \n fig = plt.figure(figsize=(15, 7))\n ax1 = fig.add_subplot(121, projection='3d')\n ax1.plot(true[:, 0], true[:, 1], true[:, 2], 'b-', label='True')\n ax1.scatter(true[:, 0], true[:, 1], true[:, 2], c='b', s=20, alpha=0.5)\n ax1.set_title('True Structure')\n ax1.set_xlabel('X')\n ax1.set_ylabel('Y')\n ax1.set_zlabel('Z')\n ax1.grid(True)\n \n ax2 = fig.add_subplot(122, projection='3d')\n ax2.plot(pred[:, 0], pred[:, 1], pred[:, 2], 'r-', label='Predicted')\n ax2.scatter(pred[:, 0], pred[:, 1], pred[:, 2], c='r', s=20, alpha=0.5)\n ax2.set_title('Predicted Structure')\n ax2.set_xlabel('X')\n ax2.set_ylabel('Y')\n ax2.set_zlabel('Z')\n ax2.grid(True)\n \n plt.suptitle(title)\n plt.tight_layout()\n \n # Always save the figure\n filename = f'structure_comparison_{sample_idx}.png'\n plt.savefig(os.path.join(OUTPUT_DIR, filename))\n \n # Only show the plot if requested\n if show_plot:\n plt.show()\n else:\n plt.close(fig)\n \n return filename # Return the filename for reference","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:17.887471Z","iopub.execute_input":"2025-05-08T13:21:17.887732Z","iopub.status.idle":"2025-05-08T13:21:17.905786Z","shell.execute_reply.started":"2025-05-08T13:21:17.887708Z","shell.execute_reply":"2025-05-08T13:21:17.904826Z"},"papermill":{"duration":0.031919,"end_time":"2025-04-16T20:30:54.354726","exception":false,"start_time":"2025-04-16T20:30:54.322807","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"2325cc7a","cell_type":"markdown","source":"> ## RNA-specific functions","metadata":{"papermill":{"duration":0.018665,"end_time":"2025-04-16T20:30:54.392167","exception":false,"start_time":"2025-04-16T20:30:54.373502","status":"completed"},"tags":[]}},{"id":"d43f6697","cell_type":"code","source":"def identify_stem_loops(sequence):\n \"\"\"\n Simple function to identify potential stem-loop regions in RNA.\n \n Parameters:\n -----------\n sequence: RNA sequence\n \n Returns:\n --------\n List of (start, end) indices for potential stem loops\n \"\"\"\n # This is a simplified implementation\n # A real implementation would use a more sophisticated algorithm\n \n stem_loops = []\n min_stem_length = 3\n \n # Look for complementary regions that could form stems\n for i in range(len(sequence) - 2*min_stem_length - 3):\n for j in range(i + min_stem_length + 3, len(sequence) - min_stem_length):\n # Check if regions could form a stem\n potential_stem = True\n for k in range(min_stem_length):\n if not are_complementary(sequence[i+k], sequence[j+min_stem_length-1-k]):\n potential_stem = False\n break\n \n if potential_stem:\n # Potential stem-loop found\n stem_loops.append((i, j + min_stem_length))\n break\n \n return stem_loops\n\ndef are_complementary(base1, base2):\n \"\"\"Check if two bases are complementary in RNA.\"\"\"\n return (base1 == 'A' and base2 == 'U') or \\\n (base1 == 'U' and base2 == 'A') or \\\n (base1 == 'G' and base2 == 'C') or \\\n (base1 == 'C' and base2 == 'G') or \\\n (base1 == 'G' and base2 == 'U') or \\\n (base1 == 'U' and base2 == 'G') # G-U wobble pairs are valid in RNA\n\ndef apply_stem_loop_template(structure, start, end):\n \"\"\"\n Apply a stem-loop template to a specific region of the structure.\n \n Parameters:\n -----------\n structure: RNA 3D structure\n start, end: Indices of the stem-loop region\n \n Returns:\n --------\n Modified structure with stem-loop template applied\n \"\"\"\n # Create a copy to modify\n result = structure.copy()\n \n # Length of the region\n region_length = end - start + 1\n \n # Not enough residues to form a proper stem-loop\n if region_length < 7:\n return result\n \n # Calculate stem length (approximately 1/3 of the region on each side)\n stem_length = max(2, region_length // 6)\n loop_start = start + stem_length\n loop_end = end - stem_length\n \n # Loop length\n loop_length = loop_end - loop_start + 1\n \n # Apply stem template (roughly parallel strands)\n for i in range(stem_length):\n # Base positions in the two stems\n pos1 = start + i\n pos2 = end - i\n \n if pos1 < len(result) and pos2 < len(result):\n # Create roughly parallel strands\n if i > 0:\n # Base the position on the previous nucleotide in the strand\n result[pos1] = result[pos1-1] + np.array([0.0, 3.8, 0.0])\n result[pos2] = result[pos2+1] + np.array([0.0, -3.8, 0.0])\n \n # Apply loop template (roughly circular)\n if loop_length > 0:\n # Calculate center of the loop\n if loop_start < len(result) and loop_end < len(result):\n center = (result[loop_start-1] + result[loop_end+1]) / 2\n center[1] += 4.0 # Offset in y direction\n \n # Create a circular loop\n radius = 3.8 # approximately nucleotide distance\n for i in range(loop_length):\n idx = loop_start + i\n if idx < len(result):\n angle = np.pi * i / (loop_length - 1)\n result[idx] = center + np.array([\n radius * np.cos(angle),\n 0.0,\n radius * np.sin(angle)\n ])\n \n return result\n\ndef post_process_rna_structure(structure, sequence, gc_content, use_global_movement=True):\n \"\"\"\n Apply RNA-specific post-processing to refine a structure.\n \n Parameters:\n -----------\n structure: Predicted 3D coordinates\n sequence: RNA sequence\n gc_content: GC content of the sequence\n use_global_movement: Whether to apply global movement transformations\n \n Returns:\n --------\n Refined structure\n \"\"\"\n # Create a new structure for modifications\n result = structure.copy()\n \n # 1. Apply mild refinement based on sequence composition\n noise_level = 0.1\n if gc_content > 0.6:\n # GC-rich regions tend to form more stable structures\n noise_level = 0.05 # Lower noise for more stable structures\n elif gc_content < 0.4:\n # AT-rich regions tend to be more flexible\n noise_level = 0.15 # Higher noise for more flexible regions\n \n # Apply noise proportional to sequence characteristics\n result = sample_structural_variation(\n result,\n noise_level=noise_level,\n preserve_distance=True, # Always preserve distances for realistic structures\n use_global_movement=use_global_movement,\n correlation=0.85 # High correlation for smoother changes\n )\n \n # 2. Look for motifs in the sequence and apply structure templates\n # This is a simplified example - a complete implementation would include more motifs\n stem_loops = identify_stem_loops(sequence)\n if stem_loops:\n for start, end in stem_loops:\n # Apply stem-loop template to these regions\n result = apply_stem_loop_template(result, start, end)\n \n # 3. Normalize bond lengths to ideal values for RNA\n valid_mask = ~np.all(result == 0, axis=1)\n for i in range(1, len(result)):\n if valid_mask[i] and valid_mask[i-1]:\n # Get the current bond vector\n bond_vector = result[i] - result[i-1]\n bond_length = np.linalg.norm(bond_vector)\n \n if bond_length > 0:\n # Normalize to ideal RNA backbone distance with small variation\n ideal_length = 3.8 * (1 + np.random.normal(0, 0.03))\n result[i] = result[i-1] + (bond_vector / bond_length) * ideal_length\n \n return result","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:17.906753Z","iopub.execute_input":"2025-05-08T13:21:17.907034Z","iopub.status.idle":"2025-05-08T13:21:17.922016Z","shell.execute_reply.started":"2025-05-08T13:21:17.907005Z","shell.execute_reply":"2025-05-08T13:21:17.921216Z"},"papermill":{"duration":0.031006,"end_time":"2025-04-16T20:30:54.441956","exception":false,"start_time":"2025-04-16T20:30:54.410950","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"9c396274","cell_type":"markdown","source":"> ## Phase 1: Golden Seeds","metadata":{"papermill":{"duration":0.017675,"end_time":"2025-04-16T20:30:54.478099","exception":false,"start_time":"2025-04-16T20:30:54.460424","status":"completed"},"tags":[]}},{"id":"3dc634c1","cell_type":"code","source":"def reference_based_approach(X_ref, y_ref, geometric_sampling=False, noise_level=0.2, correlation=0.7):\n try:\n class ReferenceModel:\n def __init__(self, geometric_sampling=False, base_noise_level=0.2, correlation=0.7):\n self.geometric_sampling = geometric_sampling\n self.base_noise_level = base_noise_level\n self.correlation = correlation\n \n def fit(self, X, y):\n # First, handle NaN values in the reference structures\n self.reference_structures = np.nan_to_num(y, nan=0.0)\n self.global_mean = np.nanmean(y, axis=(0, 1))\n self.global_std = np.nanstd(y, axis=(0, 1))\n \n # Replace potential NaN values in statistics\n self.global_mean = np.nan_to_num(self.global_mean, nan=0.0)\n self.global_std = np.nan_to_num(self.global_std, nan=1.0)\n \n # Calculate size statistics\n self.size_groups = {}\n # Group reference structures by size\n for i in range(len(self.reference_structures)):\n valid_mask = ~np.all(self.reference_structures[i] == 0, axis=1)\n size = np.sum(valid_mask)\n \n if size < 120:\n group = \"small\"\n elif size < 200:\n group = \"medium\"\n else:\n group = \"large\"\n \n if group not in self.size_groups:\n self.size_groups[group] = []\n self.size_groups[group].append(i)\n \n print(f\"Size distribution - Small: {len(self.size_groups.get('small', []))}, \"\n f\"Medium: {len(self.size_groups.get('medium', []))}, \"\n f\"Large: {len(self.size_groups.get('large', []))}\")\n \n # Store the correlation parameter for use in sample_structural_variation\n global_correlation = self.correlation\n print(f\"Using noise level: {self.base_noise_level}, correlation: {global_correlation}\")\n \n return self\n \n def predict(self, X):\n batch_size = X.shape[0]\n seq_length = X.shape[1]\n predictions = np.zeros((batch_size, seq_length, 3))\n \n for i in range(batch_size):\n # Determine the RNA size group\n valid_mask = ~np.all(X[i] == 0, axis=1)\n size = np.sum(valid_mask)\n if size < 120:\n group = \"small\"\n # Size-specific noise scaling\n noise_level = self.base_noise_level * 0.6\n elif size < 200:\n group = \"medium\"\n noise_level = self.base_noise_level * 1.0\n else:\n group = \"large\"\n noise_level = self.base_noise_level * 0.4\n \n # If we have reference structures in this size group, use them\n if group in self.size_groups and self.size_groups[group]:\n # Randomly pick a reference structure from the same size group\n ref_idx = np.random.choice(self.size_groups[group])\n base_struct = self.reference_structures[ref_idx].copy()\n \n if self.geometric_sampling:\n # Pass the correlation parameter to the variation function\n predictions[i] = sample_structural_variation(\n base_struct, \n noise_level=noise_level,\n preserve_distance=True,\n use_global_movement=(group == \"small\"),\n correlation=self.correlation\n )\n else:\n noise = np.random.normal(0, noise_level, base_struct.shape)\n predictions[i] = base_struct + noise\n else:\n # Fall back to the original method if no size match\n sample = np.random.normal(self.global_mean, self.global_std, size=(seq_length, 3))\n if self.geometric_sampling:\n predictions[i] = sample_structural_variation(\n sample, \n noise_level=noise_level,\n preserve_distance=True,\n use_global_movement=(group == \"small\"),\n correlation=self.correlation\n )\n else:\n predictions[i] = sample\n \n return predictions\n \n # Create and return model with specific parameters\n model = ReferenceModel(geometric_sampling=geometric_sampling, \n base_noise_level=noise_level,\n correlation=correlation)\n model.fit(X_ref, y_ref)\n return model\n \n except Exception as e:\n print(f\"Error in reference_based_approach: {str(e)}\")\n import traceback\n traceback.print_exc()\n return None\n\ndef evaluate_model(model, X_valid, y_valid, show_plots=False, save_top_plots=False):\n # Problem: Inadequate evaluation\n \n # SOLUTION:\n import numpy as np\n \n # Ensure there are no NaNs in the data\n X_valid_clean = np.nan_to_num(X_valid, nan=0.0)\n y_valid_clean = np.nan_to_num(y_valid, nan=0.0)\n \n # Make prediction with try/except to capture errors\n try:\n y_pred = model.predict(X_valid_clean)\n \n # Check if prediction contains NaNs or infinities\n if np.isnan(y_pred).any() or np.isinf(y_pred).any():\n print(\"WARNING: Prediction contains NaN or infinite values!\")\n y_pred = np.nan_to_num(y_pred, nan=0.0, posinf=0.0, neginf=0.0)\n \n # Calculate metrics \n mae = np.mean(np.abs(y_pred - y_valid_clean))\n mse = np.mean((y_pred - y_valid_clean)**2)\n \n # Calculate TM-scores for each structure\n tm_scores = []\n for i in range(len(X_valid)):\n # Compute score with error handling \n try:\n tm = calculate_tm_score(y_pred[i], y_valid_clean[i])\n if np.isnan(tm) or np.isinf(tm):\n print(f\"WARNING: Invalid TM-score for sample {i}, using 0.0\")\n tm = 0.0\n except Exception as e:\n print(f\"Error calculating TM-score for sample {i}: {str(e)}\")\n tm = 0.0\n \n tm_scores.append(tm)\n \n # Final metrics\n avg_tm_score = np.mean(tm_scores)\n \n print(f\"MAE: {mae:.4f}, MSE: {mse:.4f}\") \n print(f\"Average TM-score: {avg_tm_score:.4f}\")\n \n return {\n 'mae': mae,\n 'mse': mse,\n 'tm_scores': tm_scores, \n 'avg_tm_score': avg_tm_score,\n 'success': True\n }\n \n except Exception as e:\n print(f\"ERROR in evaluation: {str(e)}\")\n import traceback\n traceback.print_exc()\n \n return {\n 'mae': float('inf'),\n 'mse': float('inf'), \n 'tm_scores': [0.0] * len(X_valid),\n 'avg_tm_score': 0.0,\n 'success': False,\n 'error': str(e) \n }\n\ndef enhanced_adaptive_seed_search(\n X_valid, \n y_valid, \n initial_threshold=0.65, \n min_threshold=0.55, \n initial_attempts=100, \n max_attempts=300, \n optimal_params={'noise': 0.21, 'corr': 0.83},\n diversity_threshold=0.15,\n max_seeds=10\n):\n \"\"\"\n Enhanced adaptive seed search that combines threshold adaptation with RNA-specific targeting.\n \n This function incorporates domain knowledge about different RNA categories (small, medium, large)\n while maintaining the adaptive threshold approach to find high-quality seeds.\n \n Parameters:\n -----------\n X_valid, y_valid: Validation data\n initial_threshold: Starting TM-score threshold to consider a seed \"golden\"\n min_threshold: Minimum acceptable threshold if not enough seeds found\n initial_attempts: Initial number of attempts per category\n max_attempts: Maximum number of attempts per category\n optimal_params: Model parameters\n diversity_threshold: Threshold to consider seeds as diverse from each other\n max_seeds: Maximum number of golden seeds to return\n \n Returns:\n --------\n golden_seeds: List of diverse \"golden\" seeds\n all_seeds: List of all tested seeds with their scores\n \"\"\"\n print(f\"Enhanced adaptive search for up to {max_seeds} golden seeds starting with threshold {initial_threshold:.2f}...\")\n \n # Minimum number of seeds we want to find\n target_seeds = max(5, max_seeds // 2)\n \n current_threshold = initial_threshold\n current_attempts = initial_attempts\n \n # List to store all tested seeds\n all_seeds = []\n golden_seeds = []\n \n # Calculate sequence characteristics for RNA categorization\n seq_lengths = []\n gc_contents = []\n \n # Extract sequence lengths and nucleotide composition\n for i in range(len(X_valid)):\n # Identify non-padding positions\n valid_mask = ~np.all(X_valid[i] == 0, axis=1)\n seq_length = np.sum(valid_mask)\n seq_lengths.append(seq_length)\n \n # Extract nucleotide content if possible\n if seq_length > 0:\n features = X_valid[i][valid_mask]\n # Calculate GC content\n g_content = np.mean(features[:, 2]) if features.shape[1] > 2 else 0 # G base (index 2)\n c_content = np.mean(features[:, 1]) if features.shape[1] > 1 else 0 # C base (index 1)\n gc_content = g_content + c_content\n else:\n gc_content = 0.5 # Default if no valid sequence\n \n gc_contents.append(gc_content)\n \n # Separate indices by RNA size category\n small_rna_indices = [i for i, length in enumerate(seq_lengths) if length < 50]\n medium_rna_indices = [i for i, length in enumerate(seq_lengths) if 50 <= length < 120]\n large_rna_indices = [i for i, length in enumerate(seq_lengths) if length >= 120]\n \n # Additional categorization by GC content\n high_gc_indices = [i for i, gc in enumerate(gc_contents) if gc > 0.6]\n \n print(f\"RNA Distribution: {len(small_rna_indices)} small, {len(medium_rna_indices)} medium, {len(large_rna_indices)} large\")\n print(f\"High GC content RNAs: {len(high_gc_indices)}\")\n \n # Define optimal seed ranges and parameters for different RNA categories\n rna_categories = [\n {\n \"name\": \"small_RNA\",\n \"indices\": small_rna_indices,\n \"seed_range\": (1, 50000), # Small seeds work well for small RNAs\n \"noise_scale\": 0.9, # Lower noise for small RNAs\n \"attempts_scale\": 1.0, # Standard attempts\n \"threshold_bonus\": 0.0 # No threshold adjustment\n },\n {\n \"name\": \"medium_RNA\",\n \"indices\": medium_rna_indices,\n \"seed_range\": (10000, 150000), # Medium range for medium RNAs\n \"noise_scale\": 1.0, # Standard noise\n \"attempts_scale\": 1.2, # More attempts for medium RNAs\n \"threshold_bonus\": -0.02 # Slightly easier threshold\n },\n {\n \"name\": \"large_RNA\",\n \"indices\": large_rna_indices,\n \"seed_range\": (100000, 1000000), # Larger seeds for large RNAs\n \"noise_scale\": 0.7, # Lower noise for complex structures\n \"attempts_scale\": 1.5, # More attempts for large RNAs\n \"threshold_bonus\": -0.05 # Lower threshold - these are harder\n },\n {\n \"name\": \"high_GC\",\n \"indices\": high_gc_indices,\n \"seed_range\": (50000, 500000), # Wide range for high GC RNAs\n \"noise_scale\": 0.8, # Lower noise for stable structures\n \"attempts_scale\": 1.2, # More attempts\n \"threshold_bonus\": -0.03 # Slightly lower threshold\n }\n ]\n \n # Keep trying with lower thresholds until we find enough seeds\n while current_threshold >= min_threshold and len(golden_seeds) < target_seeds:\n print(f\"\\nSearching with base threshold {current_threshold:.2f} and {current_attempts} attempts per category...\")\n \n # Process each RNA category\n for category in rna_categories:\n category_indices = category[\"indices\"]\n \n # Skip categories with no examples\n if not category_indices:\n print(f\"Skipping {category['name']} category (no examples in validation set)\")\n continue\n \n # Adjust threshold and attempts for this category\n category_threshold = max(min_threshold, current_threshold + category[\"threshold_bonus\"])\n category_attempts = min(max_attempts, int(current_attempts * category[\"attempts_scale\"]))\n min_seed, max_seed = category[\"seed_range\"]\n \n print(f\"\\nSearching seeds for {category['name']} with threshold {category_threshold:.2f} and {category_attempts} attempts...\")\n \n # Extract subset of validation data for this category\n X_subset = [X_valid[i] for i in category_indices]\n y_subset = [y_valid[i] for i in category_indices]\n \n # List to store tested seeds at this threshold\n category_seeds = []\n \n # Set of seeds already tested to avoid duplications\n tested_seeds = set(seed_info[\"seed\"] for seed_info in all_seeds)\n \n # Counter for valid attempts (excluding duplicates)\n valid_attempts = 0\n \n # Main cycle to search for seeds in this category\n while valid_attempts < category_attempts and len(golden_seeds) < max_seeds:\n # Generate random seed from this category's range\n seed = np.random.randint(min_seed, max_seed)\n \n # Check if we've already tested this seed\n if seed in tested_seeds:\n continue\n \n tested_seeds.add(seed)\n valid_attempts += 1\n \n if valid_attempts % 10 == 0:\n print(f\"Testing seed {valid_attempts}/{category_attempts} for {category['name']} (seed={seed})...\")\n \n # Set the seed for reproducibility\n np.random.seed(seed)\n \n # Create model with this seed and category-specific parameters\n try:\n # Adjust noise level based on category\n adjusted_noise = optimal_params['noise'] * category[\"noise_scale\"]\n \n model = reference_based_approach(\n X_valid, \n y_valid,\n geometric_sampling=True,\n noise_level=adjusted_noise,\n correlation=optimal_params['corr']\n )\n \n if model is None:\n continue\n \n # First evaluate on this specific category\n category_metrics = evaluate_model(model, X_subset, y_subset)\n category_tm_score = category_metrics['avg_tm_score']\n \n # Then evaluate on all validation data\n all_metrics = evaluate_model(model, X_valid, y_valid)\n overall_tm_score = all_metrics['avg_tm_score']\n \n # Evaluate on different RNA size categories to detect overfitting\n size_scores = {}\n \n if small_rna_indices:\n small_metrics = evaluate_model_on_indices(model, X_valid, y_valid, small_rna_indices)\n size_scores[\"small\"] = small_metrics['avg_tm_score']\n \n if medium_rna_indices:\n medium_metrics = evaluate_model_on_indices(model, X_valid, y_valid, medium_rna_indices)\n size_scores[\"medium\"] = medium_metrics['avg_tm_score']\n \n if large_rna_indices:\n large_metrics = evaluate_model_on_indices(model, X_valid, y_valid, large_rna_indices)\n size_scores[\"large\"] = large_metrics['avg_tm_score']\n \n # Calculate standard deviation between scores for different sizes\n # A high deviation may indicate overfitting in certain sizes\n size_std = np.std(list(size_scores.values())) if size_scores else 0.0\n \n # Calculate balanced score that rewards:\n # 1. High overall performance\n # 2. Good performance on the specific category\n # 3. Consistent performance across RNA sizes (low size_std)\n balanced_score = (\n 0.4 * overall_tm_score + \n 0.4 * category_tm_score +\n 0.2 * (1.0 - min(1.0, size_std * 2)) # Convert std to a 0-1 score (lower is better)\n )\n \n # Register this seed\n seed_info = {\n 'seed': seed,\n 'tm_score': overall_tm_score,\n 'category_tm_score': category_tm_score,\n 'balanced_score': balanced_score,\n 'size_std': size_std,\n 'category': category['name'],\n 'size_scores': size_scores.copy()\n }\n category_seeds.append(seed_info)\n all_seeds.append(seed_info)\n \n # Check if this is a \"golden\" seed for this category\n if category_tm_score >= category_threshold:\n # Check diversity relative to seeds already found\n is_diverse = True\n \n for i, existing_seed in enumerate(golden_seeds):\n # Calculate seed similarity based on predictions\n similarity = 0\n \n # Simple similarity metric based on category and score\n if existing_seed.get('category') == category['name']:\n similarity += 0.3 # Same category adds similarity\n \n score_diff = abs(existing_seed['balanced_score'] - balanced_score)\n if score_diff < 0.1:\n similarity += (0.1 - score_diff) * 3 # Similar scores add similarity\n \n if similarity > diversity_threshold:\n is_diverse = False\n # If the new one is better than an existing one and they are similar, replace\n if balanced_score > existing_seed['balanced_score']:\n print(f\" Replacing seed {existing_seed['seed']} (score={existing_seed['balanced_score']:.4f}) \" \n f\"with seed {seed} (score={balanced_score:.4f})\")\n golden_seeds[i] = seed_info\n break\n \n if is_diverse and len(golden_seeds) < max_seeds:\n print(f\" Found golden seed: {seed} for {category['name']} (Score: {balanced_score:.4f}, Overall TM: {overall_tm_score:.4f})\")\n golden_seeds.append(seed_info)\n \n # Always check the overall performance for generalist seeds\n elif overall_tm_score >= current_threshold and balanced_score >= current_threshold:\n # Only add if not already in golden seeds and meeting diversity criteria\n if seed not in [gs['seed'] for gs in golden_seeds]:\n is_diverse = True\n \n for i, existing_seed in enumerate(golden_seeds):\n similarity = 0\n score_diff = abs(existing_seed['balanced_score'] - balanced_score)\n if score_diff < 0.1:\n similarity += (0.1 - score_diff) * 5\n \n if similarity > diversity_threshold:\n is_diverse = False\n # Replace if better\n if balanced_score > existing_seed['balanced_score']:\n print(f\" Replacing seed {existing_seed['seed']} (score={existing_seed['balanced_score']:.4f}) \" \n f\"with seed {seed} (score={balanced_score:.4f})\")\n golden_seeds[i] = seed_info\n break\n \n if is_diverse and len(golden_seeds) < max_seeds:\n print(f\" Found golden generalist seed: {seed} (Score: {balanced_score:.4f}, Overall TM: {overall_tm_score:.4f})\")\n golden_seeds.append(seed_info)\n \n except Exception as e:\n print(f\" Error testing seed {seed}: {str(e)}\")\n continue\n \n print(f\"Completed {valid_attempts} attempts for {category['name']} category\")\n \n # If we found enough seeds, we can stop\n if len(golden_seeds) >= target_seeds:\n print(f\"Found {len(golden_seeds)} golden seeds with threshold {current_threshold:.2f}\")\n break\n \n # Otherwise, reduce threshold and increase attempts\n current_threshold -= 0.03\n current_attempts = min(current_attempts + 50, max_attempts)\n \n print(f\"Not enough seeds found. Reducing threshold to {current_threshold:.2f} and increasing attempts to {current_attempts}\")\n \n # If we still don't have enough seeds, take the best ones from all tested\n if len(golden_seeds) < target_seeds and all_seeds:\n print(f\"Could not find {target_seeds} golden seeds even with threshold {current_threshold:.2f}\")\n print(f\"Using best seeds found during search based on balanced score...\")\n \n # Remove duplicates and sort by balanced score\n unique_seeds = {}\n for seed in all_seeds:\n if seed['seed'] not in unique_seeds or seed.get('balanced_score', 0) > unique_seeds[seed['seed']].get('balanced_score', 0):\n unique_seeds[seed['seed']] = seed\n \n sorted_seeds = sorted(unique_seeds.values(), key=lambda x: x.get('balanced_score', 0), reverse=True)\n \n # Add best seeds that aren't already in golden_seeds\n existing_seed_ids = {gs['seed'] for gs in golden_seeds}\n \n for seed in sorted_seeds:\n if seed['seed'] not in existing_seed_ids and len(golden_seeds) < max_seeds:\n golden_seeds.append(seed)\n existing_seed_ids.add(seed['seed'])\n \n # Ensure seeds are sorted by balanced score\n golden_seeds.sort(key=lambda x: x.get('balanced_score', 0), reverse=True)\n \n # Ensure we have at least some seeds from each major category if possible\n if len(golden_seeds) >= 3:\n categories_present = set(seed.get('category', '') for seed in golden_seeds)\n \n # Check which major categories are missing\n major_categories = ['small_RNA', 'medium_RNA', 'large_RNA']\n missing_categories = [cat for cat in major_categories if cat not in categories_present]\n \n if missing_categories:\n print(f\"Ensuring representation from missing categories: {missing_categories}\")\n \n # Find best seeds for missing categories\n for category in missing_categories:\n category_seeds = [s for s in all_seeds if s.get('category') == category]\n \n if category_seeds:\n # Sort by balanced score\n category_seeds.sort(key=lambda x: x.get('balanced_score', 0), reverse=True)\n best_seed = category_seeds[0]\n \n # Only add if not already in golden seeds\n if best_seed['seed'] not in {gs['seed'] for gs in golden_seeds}:\n # Replace worst seed if we're at capacity\n if len(golden_seeds) >= max_seeds:\n worst_idx = min(range(len(golden_seeds)), \n key=lambda i: golden_seeds[i].get('balanced_score', 0))\n print(f\" Replacing seed {golden_seeds[worst_idx]['seed']} with {best_seed['seed']} to ensure {category} representation\")\n golden_seeds[worst_idx] = best_seed\n else:\n print(f\" Adding seed {best_seed['seed']} to ensure {category} representation\")\n golden_seeds.append(best_seed)\n \n # Show statistics of the found seeds\n print(f\"\\nFound {len(golden_seeds)} golden seeds in {len(all_seeds)} total attempts\")\n \n for i, gs in enumerate(golden_seeds):\n print(f\" Seed {i+1}: {gs['seed']} (Category: {gs.get('category', 'unknown')}, Balanced: {gs.get('balanced_score', 0):.4f}, TM: {gs['tm_score']:.4f})\")\n \n # Print size-specific scores\n for size, score in gs.get('size_scores', {}).items():\n print(f\" {size.capitalize()} RNA score: {score:.4f}\")\n \n print(f\" Consistency across sizes (std): {gs.get('size_std', 0):.4f}\")\n \n return golden_seeds, all_seeds\n\ndef evaluate_model_on_indices(model, X_data, y_data, indices):\n \"\"\"\n Evaluates the model only on specific indices of the data.\n Useful to evaluate performance on subsets like small/medium/large RNAs.\n \"\"\"\n X_subset = [X_data[i] for i in indices]\n y_subset = [y_data[i] for i in indices]\n \n return evaluate_model(model, X_subset, y_subset)\n\ndef calculate_prediction_similarity(preds1, preds2):\n \"\"\"\n Calculates the similarity between two sets of predictions.\n Returns a value between 0 (totally different) and 1 (identical).\n \"\"\"\n similarities = []\n \n # For each pair of sequences in the predictions\n for p1, p2 in zip(preds1, preds2):\n # Identify valid (non-zero) coordinates\n valid_mask1 = ~np.all(p1 == 0, axis=1)\n valid_mask2 = ~np.all(p2 == 0, axis=1)\n \n # Use only positions valid in both predictions\n valid_mask = valid_mask1 & valid_mask2\n \n # If there are no overlapping valid positions, continue\n if np.sum(valid_mask) < 3:\n continue\n \n # Extract valid coordinates\n valid_p1 = p1[valid_mask]\n valid_p2 = p2[valid_mask]\n \n # Calculate similarity based on RMSD distance\n squared_diff = np.sum((valid_p1 - valid_p2) ** 2, axis=1)\n rmsd = np.sqrt(np.mean(squared_diff))\n \n # Convert RMSD to similarity (lower RMSD values = higher similarity)\n # Normalize so it's between 0 and 1\n similarity = 1.0 / (1.0 + rmsd / 5.0) # Division by 5.0 is an arbitrary scale\n similarities.append(similarity)\n \n # Return average similarity\n return np.mean(similarities) if similarities else 0.0","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:17.922985Z","iopub.execute_input":"2025-05-08T13:21:17.923215Z","iopub.status.idle":"2025-05-08T13:21:17.968667Z","shell.execute_reply.started":"2025-05-08T13:21:17.923192Z","shell.execute_reply":"2025-05-08T13:21:17.967776Z"},"papermill":{"duration":0.052282,"end_time":"2025-04-16T20:30:54.548374","exception":false,"start_time":"2025-04-16T20:30:54.496092","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"7f12acb0","cell_type":"markdown","source":"> ## Phase 2: Quality Assessment Model","metadata":{"papermill":{"duration":0.018313,"end_time":"2025-04-16T20:30:54.585647","exception":false,"start_time":"2025-04-16T20:30:54.567334","status":"completed"},"tags":[]}},{"id":"5fddd5ba","cell_type":"code","source":"class EnhancedRNAQualityNN:\n \"\"\"\n Enhanced Neural Network model for RNA structure quality assessment.\n Features:\n - Handles variable-length RNA sequences\n - Incorporates RNA-specific features\n - Attention mechanism for capturing long-range interactions\n - Multiple evaluation metrics for robust quality assessment\n \"\"\"\n def __init__(self, max_length=720):\n self.max_length = max_length\n self.is_trained = False\n self.model = None\n self.build_model()\n \n def build_model(self):\n \"\"\"\n Build an enhanced model architecture for RNA quality assessment.\n \"\"\"\n # Define the masking layer to handle variable-length sequences\n coord_input = layers.Input(shape=(self.max_length, 3), name='coordinates')\n \n # Create a mask for zero-padded coordinates\n mask_layer = layers.Lambda(\n lambda x: tf.cast(tf.reduce_sum(tf.abs(x), axis=-1) > 0.0, tf.float32),\n output_shape=lambda shape: (shape[0], shape[1])\n )\n mask = mask_layer(coord_input)\n \n # Expandir dimensões\n mask_expanded_layer = layers.Lambda(\n lambda x: tf.expand_dims(x, axis=-1),\n output_shape=lambda shape: (shape[0], shape[1], 1)\n )\n mask_expanded = mask_expanded_layer(mask) # Shape: (batch, seq_len, 1)\n \n # Optional sequence features input\n seq_input = layers.Input(shape=(self.max_length, 5), name='sequence')\n \n # Definição da função de distância pareada\n def create_pairwise_dist_layer():\n def masked_pairwise_dist_fn(inputs):\n coords, m = inputs\n # Expand dims for broadcasting\n coords1 = tf.expand_dims(coords, 2)\n coords2 = tf.expand_dims(coords, 1)\n \n # Calculate Euclidean distance\n diff = coords1 - coords2\n squared_diff = tf.reduce_sum(tf.square(diff), axis=-1)\n dist = tf.sqrt(squared_diff + 1e-8)\n \n # Create mask for valid pairs\n mask1 = tf.expand_dims(m, 2)\n mask2 = tf.expand_dims(m, 1)\n pair_mask = mask1 * mask2\n \n # Apply mask\n masked_dist = dist * pair_mask\n return masked_dist\n \n return layers.Lambda(\n masked_pairwise_dist_fn,\n output_shape=lambda shape: (shape[0][0], shape[0][1], shape[0][1])\n )\n \n # Aplicar a camada de distância pareada\n pairwise_dist_layer = create_pairwise_dist_layer()\n distances = pairwise_dist_layer([coord_input, mask])\n \n # 1.2 Process distances with 2D convolutions\n dist_features = layers.Reshape((self.max_length, self.max_length, 1))(distances)\n dist_features = layers.Conv2D(16, 3, activation='relu', padding='same')(dist_features)\n dist_features = layers.BatchNormalization()(dist_features)\n dist_features = layers.MaxPooling2D(2)(dist_features)\n \n dist_features = layers.Conv2D(32, 3, activation='relu', padding='same')(dist_features)\n dist_features = layers.BatchNormalization()(dist_features)\n dist_features = layers.MaxPooling2D(2)(dist_features)\n \n # Flatten with adaptive pooling to handle variable lengths\n dist_features = layers.GlobalAveragePooling2D()(dist_features)\n \n # 1.3 Process direct 3D coordinates with 1D convolutions\n # Apply mask to zero out padded positions\n masked_coords = layers.Multiply()([coord_input, mask_expanded])\n \n coord_features = layers.Conv1D(32, 3, activation='relu', padding='same')(masked_coords)\n coord_features = layers.BatchNormalization()(coord_features)\n \n # Para o mecanismo de auto-atenção, criamos as camadas Dense fora da função Lambda\n query_dense = layers.Dense(32)\n key_dense = layers.Dense(32)\n value_dense = layers.Dense(32)\n \n # Função de auto-atenção agora usa camadas pré-definidas\n def create_self_attention_layer(query_dense, key_dense, value_dense):\n def self_attention_fn(inputs):\n x, m = inputs\n # Simple self-attention usando camadas pré-definidas\n query = query_dense(x)\n key = key_dense(x)\n value = value_dense(x)\n \n # Calculate attention scores\n scores = tf.matmul(query, key, transpose_b=True)\n scores = scores / tf.sqrt(32.0)\n \n # Apply mask\n mask1 = tf.expand_dims(m, 2)\n mask2 = tf.expand_dims(m, 1)\n mask_2d = mask1 * mask2\n \n # Very negative number for masked positions (-1e9)\n scores = scores * mask_2d + (1.0 - mask_2d) * (-1e9)\n \n # Apply softmax\n attention_weights = tf.nn.softmax(scores, axis=-1)\n \n # Apply attention\n output = tf.matmul(attention_weights, value)\n \n return output\n \n return layers.Lambda(\n self_attention_fn,\n output_shape=lambda shape: (shape[0][0], shape[0][1], 32)\n )\n \n # Aplicar a camada de auto-atenção\n self_attention_layer = create_self_attention_layer(query_dense, key_dense, value_dense)\n attention_output = self_attention_layer([coord_features, mask])\n \n # Continue processing coordinates\n coord_features = layers.Add()([coord_features, attention_output]) # Residual connection\n coord_features = layers.Conv1D(64, 3, activation='relu', padding='same')(coord_features)\n coord_features = layers.BatchNormalization()(coord_features)\n \n # Global pooling for variable length\n coord_features = layers.GlobalAveragePooling1D()(coord_features)\n \n # 2. Process sequence information (if provided)\n seq_features = layers.Conv1D(32, 3, activation='relu', padding='same')(seq_input)\n seq_features = layers.BatchNormalization()(seq_features)\n seq_features = layers.GlobalAveragePooling1D()(seq_features)\n \n # 3. Calculate RNA-specific features\n \n # 3.1 Extract GC content and other sequence composition features\n def create_sequence_composition_layer():\n def sequence_composition_fn(inputs):\n seq, m = inputs\n # One-hot encoded sequence: (batch, len, 5) [A,C,G,U,N]\n # Calculate GC content\n c_base = seq[:, :, 1] # C base (index 1)\n g_base = seq[:, :, 2] # G base (index 2)\n \n # Sum up G and C bases and divide by sequence length\n gc_sum = tf.reduce_sum(c_base * m + g_base * m, axis=1)\n seq_length = tf.reduce_sum(m, axis=1)\n \n # Avoid division by zero\n gc_content = gc_sum / (seq_length + 1e-8)\n \n # Calculate other base contents\n a_base = seq[:, :, 0] # A base\n u_base = seq[:, :, 3] # U base\n a_content = tf.reduce_sum(a_base * m, axis=1) / (seq_length + 1e-8)\n u_content = tf.reduce_sum(u_base * m, axis=1) / (seq_length + 1e-8)\n \n # Combine features\n composition = tf.stack([gc_content, a_content, u_content], axis=1)\n \n return composition\n \n return layers.Lambda(\n sequence_composition_fn,\n output_shape=lambda shape: (shape[0][0], 3)\n )\n \n # Aplicar a camada de composição de sequência\n seq_composition_layer = create_sequence_composition_layer()\n seq_composition = seq_composition_layer([seq_input, mask])\n \n # 3.2 Calculate basic structural features\n def create_structural_features_layer():\n def structural_features_fn(inputs):\n coords, m = inputs\n # Calculate average bond length\n coords1 = coords[:, :-1, :]\n coords2 = coords[:, 1:, :]\n \n # Create mask for valid pairs\n mask_bonds = m[:, :-1] * m[:, 1:]\n mask_bonds_expanded = tf.expand_dims(mask_bonds, -1)\n \n # Calculate bond vectors and lengths\n bonds = coords2 - coords1\n masked_bonds = bonds * mask_bonds_expanded\n \n # Euclidean distance\n bond_lengths = tf.sqrt(tf.reduce_sum(tf.square(masked_bonds), axis=-1) + 1e-8)\n \n # Average bond length\n total_bonds = tf.reduce_sum(mask_bonds, axis=1)\n avg_bond_length = tf.reduce_sum(bond_lengths, axis=1) / (total_bonds + 1e-8)\n \n # Bond length consistency (std dev)\n mean_bond = tf.expand_dims(avg_bond_length, -1)\n squared_diff = tf.square(bond_lengths - mean_bond) * mask_bonds\n bond_var = tf.reduce_sum(squared_diff, axis=1) / (total_bonds + 1e-8)\n bond_std = tf.sqrt(bond_var + 1e-8)\n \n # Combine features\n struct_features = tf.stack([avg_bond_length, bond_std], axis=1)\n \n return struct_features\n \n return layers.Lambda(\n structural_features_fn,\n output_shape=lambda shape: (shape[0][0], 2)\n )\n \n # Aplicar a camada de características estruturais\n struct_features_layer = create_structural_features_layer()\n struct_features = struct_features_layer([coord_input, mask])\n \n # 4. Combine all features\n combined = layers.Concatenate()([\n dist_features, # Pairwise distance features\n coord_features, # Direct coordinate features\n seq_features, # Sequence features\n seq_composition, # GC content, etc.\n struct_features # Basic structural features\n ])\n \n # 5. Final processing with dense layers\n x = layers.Dense(128, activation='relu')(combined)\n x = layers.BatchNormalization()(x)\n x = layers.Dropout(0.3)(x)\n \n x = layers.Dense(64, activation='relu')(x)\n x = layers.BatchNormalization()(x)\n x = layers.Dropout(0.3)(x)\n \n # 6. Multiple output heads for different aspects of quality\n quality_score = layers.Dense(1, activation='sigmoid', name='quality_score')(x)\n bond_score = layers.Dense(1, activation='sigmoid', name='bond_score')(x)\n valid_score = layers.Dense(1, activation='sigmoid', name='valid_score')(x)\n \n # Create the model\n self.model = models.Model(\n inputs=[coord_input, seq_input],\n outputs=[quality_score, bond_score, valid_score]\n )\n \n # Compile with weighted losses to emphasize the overall quality score\n self.model.compile(\n optimizer=optimizers.Adam(learning_rate=1e-4, clipnorm=1.0), # Add gradient clipping\n loss={\n 'quality_score': 'mean_squared_error',\n 'bond_score': 'mean_squared_error',\n 'valid_score': 'binary_crossentropy'\n },\n loss_weights={\n 'quality_score': 1.0, # Primary loss\n 'bond_score': 0.3, # Secondary loss\n 'valid_score': 0.3 # Secondary loss\n },\n metrics={\n 'quality_score': ['mae', 'mse'],\n 'bond_score': ['mae'],\n 'valid_score': ['accuracy']\n }\n )\n \n # Os métodos train, predict_quality, save_model e load_model permanecem os mesmos\n def train(self, X_train_coords, X_train_seq, y_train, \n validation_data=None, epochs=50, batch_size=16):\n \"\"\"\n Train the model with multiple outputs.\n \n Parameters:\n -----------\n X_train_coords: Coordinate inputs (batch, seq_len, 3)\n X_train_seq: Sequence inputs (batch, seq_len, 5)\n y_train: Dictionary with 'quality_score', 'bond_score', and 'valid_score' outputs\n validation_data: Optional validation data in the same format\n \"\"\"\n # Define callbacks\n callbacks = [\n # Early stopping on the primary output - com mode='min' para métricas de perda\n EarlyStopping(\n monitor='val_quality_score_loss' if validation_data else 'quality_score_loss',\n mode='min', # Explicitamente indica que queremos minimizar a perda\n patience=10,\n restore_best_weights=True\n ),\n # Custom callback to detect and handle NaN values\n tf.keras.callbacks.TerminateOnNaN()\n ]\n \n # Train the model\n history = self.model.fit(\n x=[X_train_coords, X_train_seq],\n y=y_train,\n validation_data=validation_data,\n epochs=epochs,\n batch_size=batch_size,\n callbacks=callbacks,\n verbose=1\n )\n \n self.is_trained = True\n return history\n \n def predict_quality(self, X_coords, X_seq):\n \"\"\"\n Predict quality scores for RNA structures.\n \n Parameters:\n -----------\n X_coords: Coordinate inputs (batch, seq_len, 3)\n X_seq: Sequence inputs (batch, seq_len, 5) ou (seq_len, 5) que será expandido\n \n Returns:\n --------\n Primary quality score predictions (0-1)\n \"\"\"\n if not self.is_trained:\n print(\"WARNING: Model has not been trained yet!\")\n return None\n \n # Handle potential shape issues\n batch_size = X_coords.shape[0]\n seq_len = X_coords.shape[1]\n \n # Ensure X_seq has 3 dimensions (batch, seq_len, features)\n if len(X_seq.shape) == 2: # Se for (seq_len, features)\n X_seq = np.expand_dims(X_seq, axis=0) # Adicionar dimensão de batch\n X_seq = np.repeat(X_seq, batch_size, axis=0) # Replicar para todos os exemplos de batch\n \n # Ensure correct format for coordinates\n if seq_len > self.max_length:\n print(f\"WARNING: Input sequence length ({seq_len}) exceeds model's maximum length ({self.max_length}).\")\n print(\"Truncating input sequence to maximum length.\")\n X_coords = X_coords[:, :self.max_length, :]\n elif seq_len < self.max_length:\n print(f\"Padding input sequence from length {seq_len} to {self.max_length}\")\n padding = np.zeros((batch_size, self.max_length - seq_len, 3))\n X_coords = np.concatenate([X_coords, padding], axis=1)\n \n # Ensure correct format for sequence\n if X_seq is None:\n # If no sequence provided, create zero array\n X_seq = np.zeros((batch_size, self.max_length, 5))\n else:\n seq_shape = X_seq.shape\n if seq_shape[1] > self.max_length:\n X_seq = X_seq[:, :self.max_length, :]\n elif seq_shape[1] < self.max_length:\n padding = np.zeros((batch_size, self.max_length - seq_shape[1], 5))\n X_seq = np.concatenate([X_seq, padding], axis=1)\n \n # Predict all outputs\n outputs = self.model.predict([X_coords, X_seq])\n \n # Return the primary quality score\n return outputs[0] # quality_score output\n \n def save_model(self, filepath):\n \"\"\"Save the model to disk\"\"\"\n if self.is_trained:\n self.model.save(filepath)\n else:\n print(\"WARNING: Cannot save untrained model\")\n \n def load_model(self, filepath):\n \"\"\"Load a pre-trained model from disk\"\"\"\n self.model = models.load_model(filepath)\n self.is_trained = True","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:17.969760Z","iopub.execute_input":"2025-05-08T13:21:17.970057Z","iopub.status.idle":"2025-05-08T13:21:18.004719Z","shell.execute_reply.started":"2025-05-08T13:21:17.970033Z","shell.execute_reply":"2025-05-08T13:21:18.003899Z"},"papermill":{"duration":0.049272,"end_time":"2025-04-16T20:30:54.653439","exception":false,"start_time":"2025-04-16T20:30:54.604167","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"1651075b","cell_type":"code","source":"def prepare_multi_output_targets(train_coords, train_scores):\n \"\"\"\n Prepare multi-output target values from TM-scores.\n \n Parameters:\n -----------\n train_coords: Training coordinate data\n train_scores: TM-score values (overall quality)\n \n Returns:\n --------\n Dictionary with multiple output targets\n \"\"\"\n batch_size = len(train_scores)\n \n # Initialize targets dictionary\n targets = {\n 'quality_score': train_scores,\n 'bond_score': np.zeros((batch_size, 1)),\n 'valid_score': np.zeros((batch_size, 1))\n }\n \n # Calculate bond scores and validity scores for each structure\n for i in range(batch_size):\n coords = train_coords[i]\n \n # Calculate bond score (based on ideal bond length)\n valid_mask = ~np.all(coords == 0, axis=1)\n valid_coords = coords[valid_mask]\n \n # Skip if no valid coordinates\n if len(valid_coords) < 3:\n targets['bond_score'][i] = 0.5 # Neutral score\n targets['valid_score'][i] = 0 # Invalid\n continue\n \n # Calculate bond lengths\n bond_lengths = []\n for j in range(1, len(valid_coords)):\n dist = np.linalg.norm(valid_coords[j] - valid_coords[j-1])\n bond_lengths.append(dist)\n \n avg_bond_length = np.mean(bond_lengths)\n bond_std = np.std(bond_lengths)\n \n # Score based on how close to ideal RNA bond length (3.8Å)\n bond_score = 1.0 - min(1.0, abs(avg_bond_length - 3.8) / 3.8)\n targets['bond_score'][i] = bond_score\n \n # Validity score (binary)\n is_valid = check_structure_validity(coords)\n targets['valid_score'][i] = 1 if is_valid else 0\n \n return targets\n\ndef train_enhanced_quality_model(X_train, y_train, X_valid, y_valid):\n \"\"\"\n Train an enhanced RNA quality assessment model.\n \n Parameters:\n -----------\n X_train, X_valid: One-hot encoded RNA sequences\n y_train, y_valid: True 3D coordinates\n \n Returns:\n --------\n Trained EnhancedRNAQualityNN model\n \"\"\"\n print(\"Training enhanced RNA quality assessment model...\")\n \n # First, determine maximum sequence length in the data\n max_train_len = max(np.sum(~np.all(X_train[i] == 0, axis=1)) for i in range(len(X_train)))\n max_valid_len = max(np.sum(~np.all(X_valid[i] == 0, axis=1)) for i in range(len(X_valid)))\n max_length = max(max_train_len, max_valid_len)\n \n print(f\"Maximum sequence length in data: {max_length}\")\n \n # Adjust max_length to a reasonable value (for memory efficiency)\n max_length = min(max_length, 720) # Cap at 720 if larger\n \n # Generate training data with structure variations\n print(\"Generating training data with structure variations...\")\n \n # Parameters for data generation\n num_variations = 10 # Generate 10 variations for each structure\n \n # Containers for training data\n train_seqs = []\n train_coords = []\n train_scores = []\n \n # Process training structures\n for i in range(min(len(X_train), 50)): # Limit to 50 training examples\n print(f\"Processing training structure {i+1}/{min(len(X_train), 50)}\")\n seq_features = X_train[i]\n true_coords = y_train[i]\n \n # Check for NaN in true coordinates\n if np.isnan(true_coords).any():\n print(f\"Skipping structure {i} due to NaN in true coordinates\")\n continue\n \n # Add the true structure (highest quality)\n train_seqs.append(seq_features)\n train_coords.append(true_coords)\n train_scores.append(1.0) # Perfect score for true structure\n \n # Generate variations with different qualities\n for j in range(num_variations):\n # Vary noise level to get different quality structures\n noise_level = 0.05 + (j * 0.05) # Smaller steps for better distribution\n try:\n variation = sample_structural_variation(\n true_coords, \n noise_level=noise_level,\n preserve_distance=True, # Always preserve distances for stability\n use_global_movement=(j % 3 == 0) # Mix of global and local movements\n )\n \n # Check for NaN or Inf in variation\n if np.isnan(variation).any() or np.isinf(variation).any():\n print(f\"Skipping variation {j} for structure {i} due to NaN/Inf\")\n continue\n \n # Calculate TM-score as ground truth quality\n tm_score = calculate_tm_score(variation, true_coords)\n \n # Check if score is valid\n if np.isnan(tm_score) or np.isinf(tm_score) or tm_score <= 0:\n print(f\"Skipping variation {j} for structure {i} due to invalid TM-score: {tm_score}\")\n continue\n \n # Apply additional normalization for stability\n normalized_variation = normalize_coordinates(variation.reshape(1, -1, 3))[0]\n \n train_seqs.append(seq_features)\n train_coords.append(normalized_variation)\n train_scores.append(tm_score)\n except Exception as e:\n print(f\"Error generating variation {j} for structure {i}: {str(e)}\")\n continue\n \n # Create a smaller validation set for speed and stability\n valid_seqs = []\n valid_coords = []\n valid_scores = []\n \n for i in range(min(len(X_valid), 10)): # Use only 10 validation examples\n print(f\"Processing validation structure {i+1}/{min(len(X_valid), 10)}\")\n seq_features = X_valid[i]\n true_coords = y_valid[i]\n \n # Check for NaN in true coordinates\n if np.isnan(true_coords).any():\n print(f\"Skipping validation structure {i} due to NaN in true coordinates\")\n continue\n \n # Add the true structure\n valid_seqs.append(seq_features)\n valid_coords.append(true_coords)\n valid_scores.append(1.0)\n \n # Generate just 3 variations for validation\n for j in range(3):\n noise_level = 0.05 + (j * 0.1)\n try:\n variation = sample_structural_variation(\n true_coords, \n noise_level=noise_level,\n preserve_distance=True,\n use_global_movement=(j % 2 == 0)\n )\n \n # Check for NaN or Inf\n if np.isnan(variation).any() or np.isinf(variation).any():\n print(f\"Skipping validation variation {j} for structure {i} due to NaN/Inf\")\n continue\n \n tm_score = calculate_tm_score(variation, true_coords)\n \n # Check if score is valid\n if np.isnan(tm_score) or np.isinf(tm_score) or tm_score <= 0:\n print(f\"Skipping validation variation {j} for structure {i} due to invalid TM-score: {tm_score}\")\n continue\n \n # Apply additional normalization\n normalized_variation = normalize_coordinates(variation.reshape(1, -1, 3))[0]\n \n valid_seqs.append(seq_features)\n valid_coords.append(normalized_variation)\n valid_scores.append(tm_score)\n except Exception as e:\n print(f\"Error generating validation variation {j} for structure {i}: {str(e)}\")\n continue\n \n # Convert to numpy arrays and handle potential issues\n train_seqs = np.array(train_seqs)\n train_coords = np.array(train_coords)\n train_scores = np.array(train_scores).reshape(-1, 1) # Reshape to (n, 1)\n \n valid_seqs = np.array(valid_seqs)\n valid_coords = np.array(valid_coords)\n valid_scores = np.array(valid_scores).reshape(-1, 1) # Reshape to (n, 1)\n \n # Verify data quality and apply additional cleaning\n train_coords = np.nan_to_num(train_coords, nan=0.0, posinf=0.0, neginf=0.0)\n train_scores = np.clip(train_scores, 0.0, 1.0) # Ensure scores are in [0, 1]\n \n valid_coords = np.nan_to_num(valid_coords, nan=0.0, posinf=0.0, neginf=0.0)\n valid_scores = np.clip(valid_scores, 0.0, 1.0)\n \n # Log data statistics for debugging\n print(f\"Training data: {len(train_scores)} structures\")\n print(f\"Train coords shape: {train_coords.shape}, train scores shape: {train_scores.shape}\")\n print(f\"Train coords range: [{np.min(train_coords)}, {np.max(train_coords)}]\")\n print(f\"Train scores range: [{np.min(train_scores)}, {np.max(train_scores)}]\")\n \n print(f\"Validation data: {len(valid_scores)} structures\")\n \n try:\n # Prepare multi-output targets\n print(\"Preparing multi-output training targets...\")\n train_targets = prepare_multi_output_targets(train_coords, train_scores)\n valid_targets = prepare_multi_output_targets(valid_coords, valid_scores)\n \n # Create and train the enhanced model\n print(\"Creating and training enhanced model...\")\n model = EnhancedRNAQualityNN(max_length=max_length)\n \n # Train the model\n history = model.train(\n X_train_coords=train_coords,\n X_train_seq=train_seqs,\n y_train=train_targets,\n validation_data=([valid_coords, valid_seqs], valid_targets),\n epochs=30,\n batch_size=16\n )\n \n # Validate the model\n print(\"Validating model...\")\n val_predictions = model.predict_quality(valid_coords, valid_seqs)\n val_predictions = val_predictions.flatten()\n \n # Calculate correlation between predicted and true scores\n correlation = np.corrcoef(val_predictions, valid_scores.flatten())[0, 1]\n mae = np.mean(np.abs(val_predictions - valid_scores.flatten()))\n \n print(f\"Validation results:\")\n print(f\"Correlation: {correlation:.4f}\")\n print(f\"MAE: {mae:.4f}\")\n \n # Save the model\n os.makedirs(OUTPUT_DIR, exist_ok=True)\n model.save_model(os.path.join(OUTPUT_DIR, 'enhanced_rna_quality_model.h5'))\n \n return model\n \n except Exception as e:\n print(f\"Error training enhanced model: {str(e)}\")\n traceback.print_exc()\n \n # Fall back to a simpler model or rule-based approach\n print(\"Falling back to a simplified model due to training error...\")\n return create_rule_based_model()\n\ndef create_rule_based_model():\n \"\"\"\n Create a rule-based quality assessment model as fallback.\n \"\"\"\n class RuleBasedQualityModel:\n def __init__(self):\n self.is_trained = True\n \n def predict_quality(self, X_coords, X_seq=None):\n batch_size = X_coords.shape[0]\n \n # Implement a comprehensive rule-based quality metric\n scores = []\n \n for i in range(batch_size):\n # Check for valid coordinates\n valid_mask = ~np.all(X_coords[i] == 0, axis=1)\n coords = X_coords[i][valid_mask]\n \n if len(coords) < 3:\n scores.append(0.5) # Default score for very short structures\n continue\n \n # 1. Calculate bond lengths\n bond_lengths = []\n for j in range(1, len(coords)):\n dist = np.linalg.norm(coords[j] - coords[j-1])\n bond_lengths.append(dist)\n \n avg_bond_length = np.mean(bond_lengths)\n bond_std = np.std(bond_lengths)\n \n # 2. Score based on how close to ideal RNA bond length\n bond_score = 1.0 - min(1.0, abs(avg_bond_length - 3.8) / 3.8)\n \n # 3. Bond consistency score\n consistency_score = 1.0 - min(1.0, bond_std / 1.5)\n \n # 4. Check structure validity\n is_valid = check_structure_validity(coords)\n valid_score = 1.0 if is_valid else 0.5\n \n # 5. Check for extreme compression or expansion\n min_bond = min(bond_lengths) if bond_lengths else 0\n max_bond = max(bond_lengths) if bond_lengths else 0\n compression_score = 1.0\n if min_bond < 1.0 or max_bond > 10.0: # Physical constraints for RNA\n compression_score = 0.7\n \n # 6. Analyze radius of gyration (compactness)\n center = np.mean(coords, axis=0)\n distances = np.sqrt(np.sum((coords - center) ** 2, axis=1))\n radius_gyration = np.mean(distances)\n \n # Typical radius of gyration for RNA scales with sequence length (approximate)\n expected_radius = 3.0 * np.power(len(coords), 1/3) # Simple scaling law\n compactness_score = 1.0 - min(1.0, abs(radius_gyration - expected_radius) / expected_radius)\n \n # 7. Combined score\n final_score = (\n 0.3 * bond_score + \n 0.2 * consistency_score + \n 0.2 * valid_score + \n 0.15 * compression_score + \n 0.15 * compactness_score\n )\n \n # Ensure score is in range [0, 1]\n final_score = min(1.0, max(0.0, final_score))\n \n scores.append(final_score)\n \n return np.array(scores).reshape(-1, 1)\n \n def save_model(self, filepath):\n # Nothing to save for rule-based model\n pass\n \n return RuleBasedQualityModel()\n\ndef evaluate_and_compare_models(quality_model, rule_model, X_valid, y_valid):\n \"\"\"\n Evaluate and compare different quality assessment models.\n \n Parameters:\n -----------\n quality_model: Trained neural network model\n rule_model: Rule-based model\n X_valid, y_valid: Validation data\n \n Returns:\n --------\n Dictionary with evaluation metrics\n \"\"\"\n print(\"Evaluating and comparing quality assessment models...\")\n \n # Create validation data with multiple quality levels\n print(\"Generating validation structures with different quality levels...\")\n \n # Containers for validation data\n val_seqs = []\n val_coords = []\n val_scores = []\n \n # Number of samples to generate per structure\n num_samples = 5\n \n # Generate validation data\n for i in range(min(10, len(X_valid))):\n seq_features = X_valid[i]\n true_coords = y_valid[i]\n \n # Skip structures with NaN\n if np.isnan(true_coords).any():\n continue\n \n # Add the true structure\n val_seqs.append(seq_features)\n val_coords.append(true_coords)\n val_scores.append(1.0)\n \n # Generate variations with different quality levels\n for j in range(num_samples):\n noise_level = 0.1 * (j + 1) # Increasing noise\n \n try:\n variation = sample_structural_variation(\n true_coords,\n noise_level=noise_level,\n preserve_distance=(j % 2 == 0),\n use_global_movement=(j % 3 == 0)\n )\n \n # Skip invalid variations\n if np.isnan(variation).any() or np.isinf(variation).any():\n continue\n \n # Calculate TM-score\n tm_score = calculate_tm_score(variation, true_coords)\n \n # Skip invalid scores\n if np.isnan(tm_score) or np.isinf(tm_score) or tm_score <= 0:\n continue\n \n val_seqs.append(seq_features)\n val_coords.append(variation)\n val_scores.append(tm_score)\n \n except Exception as e:\n print(f\"Error generating validation variation: {str(e)}\")\n continue\n \n # Convert to numpy arrays\n val_coords = np.array(val_coords)\n val_seqs = np.array(val_seqs)\n val_scores = np.array(val_scores).reshape(-1, 1)\n \n print(f\"Validation data: {len(val_scores)} structures\")\n \n # Evaluate neural network model\n nn_predictions = None\n try:\n print(\"Evaluating neural network model...\")\n nn_predictions = quality_model.predict_quality(val_coords, val_seqs)\n nn_correlation = np.corrcoef(nn_predictions.flatten(), val_scores.flatten())[0, 1]\n nn_mae = np.mean(np.abs(nn_predictions.flatten() - val_scores.flatten()))\n \n print(f\"Neural network model - Correlation: {nn_correlation:.4f}, MAE: {nn_mae:.4f}\")\n except Exception as e:\n print(f\"Error evaluating neural network model: {str(e)}\")\n nn_correlation = 0.0\n nn_mae = float('inf')\n \n # Evaluate rule-based model\n rule_predictions = None\n try:\n print(\"Evaluating rule-based model...\")\n rule_predictions = rule_model.predict_quality(val_coords)\n rule_correlation = np.corrcoef(rule_predictions.flatten(), val_scores.flatten())[0, 1]\n rule_mae = np.mean(np.abs(rule_predictions.flatten() - val_scores.flatten()))\n \n print(f\"Rule-based model - Correlation: {rule_correlation:.4f}, MAE: {rule_mae:.4f}\")\n except Exception as e:\n print(f\"Error evaluating rule-based model: {str(e)}\")\n rule_correlation = 0.0\n rule_mae = float('inf')\n \n # Determine the best model\n if nn_correlation > rule_correlation:\n print(\"Neural network model performs better\")\n best_model = \"neural_network\"\n else:\n print(\"Rule-based model performs better\")\n best_model = \"rule_based\"\n \n return {\n 'neural_network': {\n 'correlation': nn_correlation,\n 'mae': nn_mae,\n 'predictions': nn_predictions\n },\n 'rule_based': {\n 'correlation': rule_correlation,\n 'mae': rule_mae,\n 'predictions': rule_predictions\n },\n 'best_model': best_model,\n 'validation_data': {\n 'coords': val_coords,\n 'scores': val_scores\n }\n }","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.005541Z","iopub.execute_input":"2025-05-08T13:21:18.005801Z","iopub.status.idle":"2025-05-08T13:21:18.036875Z","shell.execute_reply.started":"2025-05-08T13:21:18.005771Z","shell.execute_reply":"2025-05-08T13:21:18.035958Z"},"papermill":{"duration":0.051039,"end_time":"2025-04-16T20:30:54.722718","exception":false,"start_time":"2025-04-16T20:30:54.671679","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"0635d1ae","cell_type":"markdown","source":"> ## Phase 3: Base Structure Generation","metadata":{"papermill":{"duration":0.017677,"end_time":"2025-04-16T20:30:54.759006","exception":false,"start_time":"2025-04-16T20:30:54.741329","status":"completed"},"tags":[]}},{"id":"92c915c2","cell_type":"code","source":"def generate_base_structures_with_golden_seeds(\n X_test, \n test_seq_df, \n golden_seeds, \n optimal_params, \n X_valid, \n y_valid\n):\n \"\"\"\n Generate base structures using golden seeds with RNA-specific optimizations.\n \n Parameters:\n -----------\n X_test: Test features\n test_seq_df: DataFrame with test sequences\n golden_seeds: List of golden seed information\n optimal_params: Model parameters\n X_valid, y_valid: Validation data for model training\n \n Returns:\n --------\n Dictionary mapping sequence IDs to lists of base structures\n \"\"\"\n print(\"Generating base structures with golden seeds and RNA-specific optimizations...\")\n \n # Dictionary to store base structures for each sequence\n seq_to_base_structures = {}\n \n # Initialize empty base structures list for each sequence\n for _, row in test_seq_df.iterrows():\n target_id = row['target_id']\n seq_to_base_structures[target_id] = []\n \n # Sort golden seeds by TM-score for best-first approach\n sorted_seeds = sorted(golden_seeds, key=lambda x: x['tm_score'], reverse=True)\n \n # For very small RNAs, different seeds may not add much diversity\n # For large RNAs, different seeds could capture different folding patterns\n small_rna_threshold = 50 # Nucleotides\n large_rna_threshold = 200 # Nucleotides\n \n # RNA-specific parameters based on sequence properties\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n seq = row['sequence']\n seq_length = len(seq)\n \n print(f\"Processing sequence {i+1}/{len(test_seq_df)}, ID: {target_id}, length: {seq_length}\")\n \n # Extract sequence features\n seq_features = extract_sequence_features(X_test[i])\n \n # Analyze sequence to determine RNA-specific parameters\n gc_content = seq_features['gc_content']\n au_content = seq_features['au_content']\n \n # Adjust parameters based on RNA properties\n if seq_length < small_rna_threshold:\n print(f\"Small RNA detected (length={seq_length}). Using specialized parameters.\")\n num_seeds_to_use = min(3, len(sorted_seeds)) # Use fewer seeds for small RNAs\n noise_scaling = 0.7 # Lower noise for small RNAs (more stable)\n use_global_movement = False # Less global movement for small RNAs\n \n # Small RNAs with high GC content are more stable\n if gc_content > 0.6:\n noise_scaling *= 0.8 # Further reduce noise for GC-rich small RNAs\n \n elif seq_length < large_rna_threshold:\n print(f\"Medium RNA detected (length={seq_length}).\")\n num_seeds_to_use = min(4, len(sorted_seeds))\n noise_scaling = 1.0 # Standard noise level\n use_global_movement = True\n \n # For medium RNAs, GC content indicates stability regions\n if gc_content > 0.6:\n noise_scaling *= 0.9\n elif au_content > 0.6:\n noise_scaling *= 1.1 # AU-rich regions are more flexible\n \n else:\n print(f\"Large RNA detected (length={seq_length}). Using specialized parameters.\")\n num_seeds_to_use = min(5, len(sorted_seeds)) # Use more seeds for large RNAs\n noise_scaling = 0.5 # Lower noise for large RNAs (prevent unrealistic structures)\n use_global_movement = True # Use global movement for large RNAs (domain flexibility)\n \n # Large RNAs tend to have distinct domains\n # Adjust parameters to reflect domain structure\n if seq_length > 300:\n num_seeds_to_use = min(5, len(sorted_seeds)) # Maximum diversity for very large RNAs\n \n # Process with selected seeds\n base_structures = []\n for seed_idx in range(num_seeds_to_use):\n if seed_idx < len(sorted_seeds):\n seed_info = sorted_seeds[seed_idx]\n print(f\" Generating with seed {seed_info['seed']} (TM-score: {seed_info['tm_score']:.4f})\")\n \n # Set the random seed\n np.random.seed(seed_info['seed'])\n \n # Create model with adjusted parameters\n adjusted_noise = optimal_params['noise'] * noise_scaling\n \n # Create model with RNA-specific adjustments\n model = reference_based_approach(\n X_valid, \n y_valid,\n geometric_sampling=True, # Always use geometric sampling for better structures\n noise_level=adjusted_noise,\n correlation=optimal_params['corr']\n )\n \n if model is None:\n print(f\" Failed to create model with seed {seed_info['seed']}\")\n continue\n \n # Generate prediction\n try:\n # Get basic prediction for this sequence\n base_pred = model.predict(X_test[i:i+1])[0][:seq_length]\n \n # Apply RNA-specific post-processing\n processed_pred = post_process_rna_structure(\n base_pred, \n seq, \n gc_content, \n use_global_movement=use_global_movement\n )\n \n # Normalize the structure\n normalized_pred = normalize_structure(processed_pred)\n \n # Verify the structure meets basic validation criteria\n if check_structure_validity(normalized_pred):\n base_structures.append(normalized_pred)\n else:\n print(f\" Structure from seed {seed_info['seed']} failed validation. Attempting repair.\")\n \n # Try to repair the structure\n repaired_structure = repair_invalid_structure(normalized_pred)\n if check_structure_validity(repaired_structure):\n base_structures.append(repaired_structure)\n print(f\" Successfully repaired structure from seed {seed_info['seed']}\")\n else:\n print(f\" Could not repair structure from seed {seed_info['seed']}\")\n \n except Exception as e:\n print(f\" Error generating prediction with seed {seed_info['seed']}: {str(e)}\")\n continue\n \n # If we didn't get any valid structures, create an emergency structure\n if not base_structures:\n print(f\"Warning: No valid structures generated for {target_id}. Creating emergency structure.\")\n emergency_structure = create_emergency_structure(seq_length)\n base_structures.append(emergency_structure)\n \n # Store the structures\n seq_to_base_structures[target_id] = base_structures\n print(f\" Generated {len(base_structures)} base structures for {target_id}\")\n \n return seq_to_base_structures","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.037780Z","iopub.execute_input":"2025-05-08T13:21:18.038072Z","iopub.status.idle":"2025-05-08T13:21:18.055022Z","shell.execute_reply.started":"2025-05-08T13:21:18.038047Z","shell.execute_reply":"2025-05-08T13:21:18.054234Z"},"papermill":{"duration":0.029964,"end_time":"2025-04-16T20:30:54.806721","exception":false,"start_time":"2025-04-16T20:30:54.776757","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"30307214","cell_type":"markdown","source":"> ## Phase 4: Candidate Generation and Pruning","metadata":{"papermill":{"duration":0.018652,"end_time":"2025-04-16T20:30:54.843577","exception":false,"start_time":"2025-04-16T20:30:54.824925","status":"completed"},"tags":[]}},{"id":"d6590716","cell_type":"code","source":"def generate_diverse_candidates(base_structures, seq_length, num_per_base=5):\n \"\"\"\n Generate diverse candidate structures from a set of base structures.\n Adapts variation parameters based on RNA size.\n \n Parameters:\n -----------\n base_structures: List of base structures to generate variations from\n seq_length: Length of the sequence\n num_per_base: Number of variations to generate per base structure\n \n Returns:\n --------\n List of candidate structures\n \"\"\"\n candidates = []\n \n # First, add all base structures\n for base in base_structures:\n candidates.append(base)\n \n # Then generate variations from each base\n for base_idx, base in enumerate(base_structures):\n print(f\" Generating variations from base structure {base_idx+1}/{len(base_structures)}...\")\n \n # Determine noise levels based on sequence length\n if seq_length < 50:\n # Small RNA - can handle more variation\n noise_levels = [0.1, 0.2, 0.3, 0.4, 0.5]\n elif seq_length < 120:\n # Medium RNA - moderate variation\n noise_levels = [0.05, 0.1, 0.15, 0.2, 0.25]\n else:\n # Large RNA - more conservative\n noise_levels = [0.03, 0.06, 0.09, 0.12, 0.15]\n \n # Generate variations with different parameters\n for i in range(num_per_base):\n # Use different parameters for diversity\n noise_idx = i % len(noise_levels)\n noise_level = noise_levels[noise_idx]\n preserve_distance = (i % 2 == 0) # Alternate between preserving and not\n use_global = (i % 3 == 0) # Occasional global movements\n \n # Add small random variation to correlation\n correlation = 0.8 + np.random.uniform(-0.1, 0.1)\n \n # Set a unique random seed for each variation\n np.random.seed(base_idx * 100 + i)\n \n variation = sample_structural_variation(\n base,\n noise_level=noise_level,\n preserve_distance=preserve_distance,\n use_global_movement=use_global,\n correlation=correlation\n )\n \n # Normalize the structure\n normalized = normalize_structure(variation)\n candidates.append(normalized)\n \n print(f\"Generated {len(candidates)} candidate structures in total\")\n return candidates\n\ndef generate_diverse_structures_from_bases(base_structures, seq_length, quality_model, num_per_base=5):\n \"\"\"\n Generate diverse candidate structures from a set of base structures,\n with RNA-specific variations and quality filtering.\n \n Parameters:\n -----------\n base_structures: List of base structures to generate variations from\n seq_length: Length of the RNA sequence\n quality_model: Model for quality assessment\n num_per_base: Number of variations to generate per base structure\n \n Returns:\n --------\n List of diverse candidate structures\n \"\"\"\n candidates = []\n \n # First, add all base structures\n for base in base_structures:\n candidates.append(base)\n \n # Then generate variations from each base\n for base_idx, base in enumerate(base_structures):\n print(f\" Generating variations from base structure {base_idx+1}/{len(base_structures)}...\")\n \n # Determine variation parameters based on sequence length\n if seq_length < 50:\n # Small RNA - can handle more variation\n noise_levels = [0.05, 0.1, 0.15, 0.2, 0.25]\n preserve_distances = [True, True, True, False, False] # Mostly preserve distances\n use_globals = [False, False, True, False, True] # Occasional global movements\n elif seq_length < 120:\n # Medium RNA - moderate variation\n noise_levels = [0.03, 0.06, 0.1, 0.15, 0.2]\n preserve_distances = [True, True, True, True, False] # Mostly preserve distances\n use_globals = [False, True, False, True, False] # Mix of global and local\n else:\n # Large RNA - more conservative\n noise_levels = [0.02, 0.04, 0.06, 0.08, 0.1]\n preserve_distances = [True, True, True, True, True] # Always preserve distances\n use_globals = [False, False, True, False, True] # Occasional global for domains\n \n # Generate variations with different parameters\n for i in range(num_per_base):\n # Use different parameters for diversity\n noise_idx = i % len(noise_levels)\n noise_level = noise_levels[noise_idx]\n preserve_distance = preserve_distances[noise_idx]\n use_global = use_globals[noise_idx] \n \n # Add small random variation to correlation\n correlation = 0.8 + np.random.uniform(-0.1, 0.1)\n \n # Set a unique random seed for each variation\n np.random.seed(base_idx * 100 + i)\n \n variation = sample_structural_variation(\n base,\n noise_level=noise_level,\n preserve_distance=preserve_distance,\n use_global_movement=use_global,\n correlation=correlation\n )\n \n # Apply additional RNA-specific refinements\n # For example, ensure proper backbone geometry\n variation = refine_rna_backbone(variation)\n \n # Normalize the structure\n normalized = normalize_structure(variation)\n \n # Verify the structure is valid\n if check_structure_validity(normalized):\n candidates.append(normalized)\n else:\n print(f\" Structure failed validation. Attempting repair.\")\n repaired = repair_invalid_structure(normalized)\n if check_structure_validity(repaired):\n candidates.append(repaired)\n print(f\" Successfully repaired structure\")\n \n print(f\"Generated {len(candidates)} candidate structures in total\")\n \n # Pre-filter candidates based on quality before detailed evaluation\n if len(candidates) > 30: # Only pre-filter if we have many candidates\n print(\"Pre-filtering candidates based on basic quality metrics...\")\n quality_scores = []\n \n # Simple quality assessment for pre-filtering\n for candidate in candidates:\n # Calculate basic quality score\n valid_mask = ~np.all(candidate == 0, axis=1)\n valid_coords = candidate[valid_mask]\n \n # Skip if too few valid coordinates\n if len(valid_coords) < 3:\n quality_scores.append(0.0)\n continue\n \n # Calculate bond lengths\n bond_lengths = []\n for j in range(1, len(valid_coords)):\n dist = np.linalg.norm(valid_coords[j] - valid_coords[j-1])\n bond_lengths.append(dist)\n \n # Score based on ideal bond length\n avg_bond_length = np.mean(bond_lengths)\n bond_score = 1.0 - min(1.0, abs(avg_bond_length - 3.8) / 3.8)\n \n quality_scores.append(bond_score)\n \n # Convert to numpy array\n quality_scores = np.array(quality_scores)\n \n # Take top 30 candidates based on quality score\n top_indices = np.argsort(quality_scores)[-30:]\n candidates = [candidates[idx] for idx in top_indices]\n print(f\"Pre-filtered to top 30 candidates\")\n \n return candidates\n\ndef evaluate_and_prune_structures(candidates, seq_features, quality_model, top_k=5):\n \"\"\"\n Evaluate structure candidates and select the top-k structures.\n This function handles both NN-based and rule-based quality models.\n \n Parameters:\n -----------\n candidates: List of candidate structures\n seq_features: RNA sequence features\n quality_model: Model for quality assessment\n top_k: Number of top structures to select\n \n Returns:\n --------\n List of top-k structures\n \"\"\"\n # Determine if the model is a neural network or rule-based\n is_nn_model = hasattr(quality_model, 'model')\n \n try:\n if is_nn_model:\n print(\"Using neural network for quality assessment...\")\n return evaluate_and_prune_nn(candidates, seq_features, quality_model, top_k)\n else:\n print(\"Using rule-based model for quality assessment...\")\n return evaluate_and_prune_rules(candidates, top_k)\n \n except Exception as e:\n print(f\"Error during quality evaluation: {str(e)}\")\n traceback.print_exc()\n \n # Fall back to rule-based evaluation if any error occurs\n print(\"Falling back to basic rule-based scoring...\")\n return evaluate_and_prune_rules(candidates, top_k)","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.055896Z","iopub.execute_input":"2025-05-08T13:21:18.056172Z","iopub.status.idle":"2025-05-08T13:21:18.072917Z","shell.execute_reply.started":"2025-05-08T13:21:18.056130Z","shell.execute_reply":"2025-05-08T13:21:18.072226Z"},"papermill":{"duration":0.033575,"end_time":"2025-04-16T20:30:54.895124","exception":false,"start_time":"2025-04-16T20:30:54.861549","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"78fdd1e0","cell_type":"code","source":"def evaluate_and_prune_nn(candidates, seq_features, quality_model, top_k=5):\n \"\"\"\n Evaluate candidates using NN model and select the top-k.\n \n Parameters:\n -----------\n candidates: List of candidate structures\n seq_features: One-hot encoded sequence features\n quality_model: Trained quality assessment model\n top_k: Number of top structures to keep\n \n Returns:\n --------\n List of top-k structures\n \"\"\"\n try:\n # Extract actual sequence length (non-padding)\n valid_mask = ~np.all(seq_features == 0, axis=1)\n seq_length = np.sum(valid_mask)\n \n # Prepare batched data for prediction\n stacked_candidates = np.array(candidates)\n \n # Prepare sequence features input - deve ter o mesmo número de amostras que stacked_candidates\n batch_size = stacked_candidates.shape[0]\n \n # Expand seq_features to have batch_size samples (replicando para cada candidato)\n # Certifique-se de que seq_features tem 3 dimensões (batch, seq_len, features)\n if len(seq_features.shape) == 2: # Se for (seq_len, features)\n seq_features = np.expand_dims(seq_features, axis=0) # Adicionar dimensão de batch\n \n # Replicar para todos os candidatos\n stacked_seq = np.repeat(seq_features, batch_size, axis=0)\n \n # Predict quality scores\n quality_scores = quality_model.predict_quality(stacked_candidates, stacked_seq)\n quality_scores = quality_scores.flatten()\n \n # Sort by quality score\n sorted_indices = np.argsort(quality_scores)[::-1] # Descending order\n \n # Keep top-k structures\n top_structures = [candidates[idx] for idx in sorted_indices[:top_k]]\n top_scores = quality_scores[sorted_indices[:top_k]]\n \n print(f\"Selected top {top_k} structures with NN predicted qualities: {top_scores}\")\n \n return top_structures\n \n except Exception as e:\n print(f\"Error in NN evaluation: {str(e)}\")\n traceback.print_exc()\n \n # Fall back to rule-based approach if NN fails\n print(\"Falling back to rule-based evaluation...\")\n return evaluate_and_prune_rules(candidates, top_k)\n\ndef evaluate_and_prune_rules(candidates, top_k=5):\n \"\"\"\n Evaluate candidates using rule-based metrics and select the top-k.\n \n Parameters:\n -----------\n candidates: List of candidate structures\n top_k: Number of top structures to keep\n \n Returns:\n --------\n List of top-k structures\n \"\"\"\n quality_scores = []\n \n for i, candidate in enumerate(candidates):\n # Calculate a quality score based on structural features\n # 1. Check for valid coordinates\n valid_mask = ~np.all(candidate == 0, axis=1)\n valid_coords = candidate[valid_mask]\n \n # Skip if no valid coordinates\n if len(valid_coords) < 3:\n quality_scores.append(0.5) # Neutral score\n continue\n \n # 2. Calculate bond lengths between consecutive residues\n bond_lengths = []\n for j in range(1, len(valid_coords)):\n dist = np.linalg.norm(valid_coords[j] - valid_coords[j-1])\n bond_lengths.append(dist)\n \n avg_bond_length = np.mean(bond_lengths)\n bond_std = np.std(bond_lengths)\n \n # 3. Score based on how close to ideal RNA bond length (3.8Å)\n bond_score = 1.0 - min(1.0, abs(avg_bond_length - 3.8) / 3.8)\n \n # 4. Bond consistency score (lower std deviation is better)\n consistency_score = 1.0 - min(1.0, bond_std / 2.0)\n \n # 5. Check structure validity\n is_valid = check_structure_validity(candidate)\n valid_score = 1.0 if is_valid else 0.5\n \n # 6. Combined score\n score = 0.4 * bond_score + 0.3 * consistency_score + 0.3 * valid_score\n \n # 7. Add small random component for variations\n random_component = np.random.uniform(-0.05, 0.05)\n score = min(1.0, max(0.0, score + random_component))\n \n quality_scores.append(score)\n \n # Convert to numpy array\n quality_scores = np.array(quality_scores)\n \n # Sort by quality score\n sorted_indices = np.argsort(quality_scores)[::-1] # Descending order\n \n # Keep top-k structures\n top_structures = [candidates[idx] for idx in sorted_indices[:top_k]]\n top_scores = quality_scores[sorted_indices[:top_k]]\n \n print(f\"Selected top {top_k} structures with rule-based qualities: {top_scores}\")\n \n return top_structures\n\ndef generate_and_prune_structures(base_coords, seq_features, quality_model, num_candidates=20, top_k=5):\n \"\"\"\n Generate multiple structure candidates and use the NN model to prune to the best ones.\n Modified to handle variable-length RNA sequences.\n \"\"\"\n # Get actual sequence length (non-padding)\n valid_mask = ~np.all(base_coords == 0, axis=1)\n seq_length = np.sum(valid_mask)\n print(f\"Processing structure with actual length: {seq_length}\")\n \n # Generate candidate structures with different parameters\n candidates = []\n \n # Add the base structure\n candidates.append(normalize_structure(base_coords))\n \n # Generate variations with different parameters\n for i in range(num_candidates - 1):\n # Use different parameters for diversity\n noise_level = 0.1 + (i % 10) * 0.05\n preserve_distance = (i % 3 != 0)\n use_global = (i % 4 == 0)\n correlation = 0.7 + (i % 5) * 0.05\n \n variation = sample_structural_variation(\n base_coords,\n noise_level=noise_level,\n preserve_distance=preserve_distance,\n use_global_movement=use_global,\n correlation=correlation\n )\n \n # Normalize the structure\n normalized = normalize_structure(variation)\n candidates.append(normalized)\n \n # Convert to array for batch processing\n stacked_candidates = np.array(candidates)\n \n # Implement a simple rule-based quality assessment as fallback\n print(\"Using rule-based quality assessment...\")\n quality_scores = []\n \n for i, candidate in enumerate(candidates):\n # Calculate a quality score based on structural features\n # 1. Check for unusual bond lengths\n valid_indices = np.where(valid_mask)[0]\n valid_coords = candidate[valid_indices]\n \n # Skip if no valid coordinates\n if len(valid_coords) < 3:\n quality_scores.append(0.5)\n continue\n \n # Calculate bond lengths\n bond_lengths = []\n for j in range(1, len(valid_coords)):\n dist = np.linalg.norm(valid_coords[j] - valid_coords[j-1])\n bond_lengths.append(dist)\n \n # Score based on how close to ideal RNA bond length\n avg_bond_length = np.mean(bond_lengths)\n bond_std = np.std(bond_lengths)\n \n # Ideal bond length is around 3.8Å\n bond_score = 1.0 - min(1.0, abs(avg_bond_length - 3.8) / 3.8)\n \n # Bond consistency score\n consistency_score = 1.0 - min(1.0, bond_std / 2.0)\n \n # Structural validity\n is_valid = check_structure_validity(candidate)\n valid_score = 1.0 if is_valid else 0.5\n \n # Combined score\n final_score = 0.4 * bond_score + 0.3 * consistency_score + 0.3 * valid_score\n \n # Add a small random component for variations\n random_component = np.random.uniform(-0.05, 0.05)\n final_score = min(1.0, max(0.0, final_score + random_component))\n \n quality_scores.append(final_score)\n \n quality_scores = np.array(quality_scores)\n \n # Sort candidates by quality score\n sorted_indices = np.argsort(quality_scores)[::-1] # Descending order\n \n # Keep top-k structures\n top_structures = [candidates[idx] for idx in sorted_indices[:top_k]]\n top_scores = quality_scores[sorted_indices[:top_k]]\n \n print(f\"Selected top {top_k} structures with predicted qualities: {top_scores}\")\n \n return top_structures","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.073853Z","iopub.execute_input":"2025-05-08T13:21:18.074105Z","iopub.status.idle":"2025-05-08T13:21:18.091434Z","shell.execute_reply.started":"2025-05-08T13:21:18.074074Z","shell.execute_reply":"2025-05-08T13:21:18.090553Z"},"papermill":{"duration":0.034326,"end_time":"2025-04-16T20:30:54.948135","exception":false,"start_time":"2025-04-16T20:30:54.913809","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"6602b41e","cell_type":"markdown","source":"> ## Phase 5: Submission Creation","metadata":{"papermill":{"duration":0.017773,"end_time":"2025-04-16T20:30:54.983837","exception":false,"start_time":"2025-04-16T20:30:54.966064","status":"completed"},"tags":[]}},{"id":"e60d1e8d","cell_type":"code","source":"def create_submission_dataframe(seq_to_coords, sample_submission_df):\n \"\"\"\n Create a submission DataFrame from the final structures.\n \n Parameters:\n -----------\n seq_to_coords: Dictionary mapping sequence IDs to lists of structures\n sample_submission_df: Sample submission format\n \n Returns:\n --------\n Submission DataFrame\n \"\"\"\n # Create a copy of the sample submission\n submission_df = sample_submission_df.copy()\n \n # Fill in the coordinates for each structure\n for i, row in submission_df.iterrows():\n if i % 1000 == 0:\n print(f\"Processing row {i}/{len(submission_df)}\")\n \n # Parse the ID to get sequence ID and residue index\n id_parts = row['ID'].split('_')\n seq_id = id_parts[0]\n residue_idx = int(id_parts[1]) - 1 # Convert to 0-based indexing\n \n # Check if we have structures for this sequence\n if seq_id in seq_to_coords:\n structures = seq_to_coords[seq_id]\n \n # Check if the residue index is valid\n if residue_idx < len(structures[0]):\n # Fill in coordinates for all 5 structures\n for struct_idx in range(5):\n if struct_idx < len(structures):\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[struct_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[struct_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[struct_idx][residue_idx][2]\n else:\n # If we have fewer than 5 structures, duplicate the last one\n last_idx = len(structures) - 1\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[last_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[last_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[last_idx][residue_idx][2]\n \n return submission_df\n\ndef generate_nn_pruned_submission(model, quality_model, test_seq_df, sample_submission_df):\n \"\"\"\n Enhanced submission generation that uses NN-based pruning for structure selection.\n \"\"\"\n print(\"Generating submission with Neural Network pruning...\")\n \n # Prepare test features\n X_test = prepare_test_features(test_seq_df)\n \n # Generate multiple predictions for ensemble diversity\n print(\"Generating base predictions...\")\n base_predictions = model.predict(X_test)\n \n seq_to_coords = {}\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n seq = row['sequence']\n seq_length = len(seq)\n \n print(f\"Processing sequence {i+1}/{len(test_seq_df)}, ID: {target_id}, length: {seq_length}\")\n \n # Get base coordinates\n base_coords = base_predictions[i][:seq_length]\n \n # Extract sequence features for this RNA\n seq_features = X_test[i][:seq_length]\n \n # Generate and prune structures using the NN model\n structures = generate_and_prune_structures(\n base_coords, \n seq_features, \n quality_model,\n num_candidates=30, # Generate more candidates\n top_k=5 # Keep top 5 for submission\n )\n \n # Store the structures\n seq_to_coords[target_id] = structures\n \n # Create submission DataFrame\n print(\"Creating submission file...\")\n submission_df = sample_submission_df.copy()\n \n for i, row in submission_df.iterrows():\n id_parts = row['ID'].split('_')\n seq_id = id_parts[0]\n residue_idx = int(id_parts[1]) - 1\n \n if seq_id in seq_to_coords:\n structures = seq_to_coords[seq_id]\n if residue_idx < len(structures[0]):\n for struct_idx in range(5):\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[struct_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[struct_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[struct_idx][residue_idx][2]\n \n submission_file = os.path.join(OUTPUT_DIR, 'submission_nn_pruned.csv')\n submission_df.to_csv(submission_file, index=False)\n print(f\"NN-pruned submission file saved to {submission_file}\")\n \n # Also save as standard submission\n standard_file = os.path.join(OUTPUT_DIR, 'submission.csv')\n submission_df.to_csv(standard_file, index=False)\n \n return submission_df","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.092218Z","iopub.execute_input":"2025-05-08T13:21:18.092479Z","iopub.status.idle":"2025-05-08T13:21:18.110057Z","shell.execute_reply.started":"2025-05-08T13:21:18.092455Z","shell.execute_reply":"2025-05-08T13:21:18.109220Z"},"papermill":{"duration":0.030533,"end_time":"2025-04-16T20:30:55.032472","exception":false,"start_time":"2025-04-16T20:30:55.001939","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"17a96be1","cell_type":"markdown","source":"> ## Main Pipeline Functions","metadata":{"papermill":{"duration":0.018143,"end_time":"2025-04-16T20:30:55.069389","exception":false,"start_time":"2025-04-16T20:30:55.051246","status":"completed"},"tags":[]}},{"id":"c03234bf","cell_type":"code","source":"def run_hybrid_pipeline(\n X_valid, \n y_valid, \n test_seq_df, \n sample_submission_df, \n output_dir, \n golden_threshold=0.6, \n optimal_params={'noise': 0.21, 'corr': 0.83}\n):\n \"\"\"\n Run a hybrid pipeline that uses RNA-targeted seed selection.\n \n Parameters:\n -----------\n X_valid, y_valid: Validation data\n test_seq_df: DataFrame with test sequences\n sample_submission_df: Sample submission format\n output_dir: Output directory for files\n golden_threshold: Threshold for considering a seed as \"golden\"\n optimal_params: Optimal parameters for the reference model\n \n Returns:\n --------\n submission_df, status_dict\n \"\"\"\n print(\"=\" * 80)\n print(\"HYBRID PIPELINE WITH RNA-TARGETED SEED SELECTION\".center(80))\n print(\"=\" * 80)\n \n status = {\n 'success': False,\n 'golden_seeds_found': 0,\n 'nn_training_success': False,\n 'best_tm_score': 0.0,\n 'error': None\n }\n \n try:\n # PHASE 1: Find Enhanced Adaptive Seeds\n print(\"\\nPHASE 1: Searching for enhanced adaptive seeds...\")\n golden_seeds, all_seeds = enhanced_adaptive_seed_search(\n X_valid, \n y_valid, \n initial_threshold=golden_threshold,\n min_threshold=golden_threshold - 0.1,\n optimal_params=optimal_params\n )\n \n # Use fallback approach if we didn't find enough seeds\n if len(golden_seeds) < 3:\n print(\"Not enough enhanced seeds found, falling back to general seed search...\")\n golden_seeds, all_seeds = find_diverse_golden_seeds(\n X_valid, \n y_valid, \n golden_threshold=golden_threshold - 0.05, # Lower threshold for fallback\n attempts=100\n )\n \n ensemble_seeds = golden_seeds\n status['golden_seeds_found'] = len(ensemble_seeds)\n \n # PHASE 2: Train Quality Assessment Model\n print(\"\\nPHASE 2: Training quality assessment model...\")\n try:\n quality_model = train_enhanced_quality_model(X_valid, y_valid, X_valid, y_valid)\n status['nn_training_success'] = True\n except Exception as e:\n print(f\"Error training NN model: {str(e)}\")\n print(\"Falling back to rule-based quality assessment...\")\n quality_model = create_rule_based_model()\n \n # PHASE 3: Generate Base Structures with RNA-Targeted Seeds\n print(\"\\nPHASE 3: Generating base structures with enhanced seeds...\")\n X_test = prepare_test_features(test_seq_df)\n \n # Generate predictions using each of the ensemble seeds\n seed_predictions = []\n for i, seed_info in enumerate(ensemble_seeds):\n print(f\"Generating predictions with seed {seed_info['seed']} (TM-score: {seed_info['tm_score']:.4f}, Category: {seed_info.get('category', 'unknown')})...\")\n \n # Set the random seed\n np.random.seed(seed_info['seed'])\n \n # Create model with this seed\n model = reference_based_approach(\n X_valid, \n y_valid,\n geometric_sampling=True,\n noise_level=optimal_params['noise'],\n correlation=optimal_params['corr']\n )\n \n # Generate predictions\n if model is not None:\n preds = model.predict(X_test)\n seed_predictions.append({\n 'seed': seed_info['seed'],\n 'tm_score': seed_info['tm_score'],\n 'category': seed_info.get('category', 'unknown'),\n 'predictions': preds\n })\n \n # Update best TM-score for status\n if seed_info['tm_score'] > status['best_tm_score']:\n status['best_tm_score'] = seed_info['tm_score']\n else:\n print(f\"Failed to create model with seed {seed_info['seed']}\")\n \n if not seed_predictions:\n raise Exception(\"Failed to generate any predictions with enhanced seeds\")\n \n # PHASE 4: Generate and Prune Structures\n print(\"\\nPHASE 4: Generating diverse candidates and pruning...\")\n \n seq_to_coords = {}\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n seq = row['sequence']\n seq_length = len(seq)\n \n print(f\"Processing sequence {i+1}/{len(test_seq_df)}, ID: {target_id}, length: {seq_length}\")\n \n # Collect base predictions from all seeds for this sequence\n base_structures = []\n for pred_info in seed_predictions:\n base_struct = pred_info['predictions'][i][:seq_length]\n \n # Apply RNA-specific post-processing based on seed category\n # Extract GC content for RNA-specific processing\n seq_features = X_test[i][:seq_length]\n valid_mask = ~np.all(seq_features == 0, axis=1)\n if np.sum(valid_mask) > 0:\n features = seq_features[valid_mask]\n g_content = np.mean(features[:, 2]) if features.shape[1] > 2 else 0\n c_content = np.mean(features[:, 1]) if features.shape[1] > 1 else 0\n gc_content = g_content + c_content\n else:\n gc_content = 0.5\n \n # Apply RNA-specific post-processing\n category = pred_info.get('category', 'unknown')\n use_global_movement = category in ['large_RNA', 'medium_RNA']\n \n # Process structure with RNA-specific adjustments\n processed_struct = post_process_rna_structure(\n base_struct,\n seq,\n gc_content,\n use_global_movement=use_global_movement\n )\n \n # Normalize structure\n normalized_struct = normalize_structure(processed_struct)\n \n # Verify validity and repair if needed\n if check_structure_validity(normalized_struct):\n base_structures.append(normalized_struct)\n else:\n # Try to repair the structure\n repaired_struct = repair_invalid_structure(normalized_struct)\n if check_structure_validity(repaired_struct):\n base_structures.append(repaired_struct)\n print(f\" Repaired structure from seed {pred_info['seed']}\")\n else:\n print(f\" Structure from seed {pred_info['seed']} failed validation and repair\")\n \n # Generate emergency structure if no valid base structures\n if not base_structures:\n print(f\"Warning: No valid structures generated for {target_id}. Creating emergency structure.\")\n emergency_structure = create_emergency_structure(seq_length)\n base_structures.append(emergency_structure)\n \n # Extract sequence features for RNA-specific candidate generation\n seq_features = X_test[i][:seq_length]\n \n # Generate diverse candidates optimized for this RNA\n candidates = generate_diverse_structures_from_bases(\n base_structures, \n seq_length, \n quality_model,\n num_per_base=5\n )\n \n # Evaluate and prune candidates\n if status['nn_training_success']:\n print(\"Using NN model for quality assessment...\")\n try:\n top_structures = evaluate_and_prune_structures(\n candidates, \n seq_features, \n quality_model, \n top_k=5\n )\n except Exception as e:\n print(f\"Error in NN evaluation: {str(e)}\")\n print(\"Falling back to rule-based assessment...\")\n top_structures = evaluate_and_prune_rules(candidates, top_k=5)\n else:\n print(\"Using rule-based quality assessment...\")\n top_structures = evaluate_and_prune_rules(candidates, top_k=5)\n \n # Store the final structures\n seq_to_coords[target_id] = top_structures\n \n # PHASE 5: Create Submission\n print(\"\\nPHASE 5: Creating submission file...\")\n submission_df = create_submission_dataframe(seq_to_coords, sample_submission_df)\n \n # Save submission files\n enhanced_file = os.path.join(output_dir, 'submission_enhanced.csv')\n submission_df.to_csv(enhanced_file, index=False)\n print(f\"Enhanced submission saved to {enhanced_file}\")\n \n # Save as standard submission\n standard_file = os.path.join(output_dir, 'submission.csv')\n submission_df.to_csv(standard_file, index=False)\n print(f\"Standard submission saved to {standard_file}\")\n \n # Set success\n status['success'] = True\n \n return submission_df, status\n \n except Exception as e:\n print(f\"ERROR in hybrid pipeline: {str(e)}\")\n traceback.print_exc()\n status['error'] = str(e)\n return None, status","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.111071Z","iopub.execute_input":"2025-05-08T13:21:18.111372Z","iopub.status.idle":"2025-05-08T13:21:18.133129Z","shell.execute_reply.started":"2025-05-08T13:21:18.111345Z","shell.execute_reply":"2025-05-08T13:21:18.132213Z"},"papermill":{"duration":0.032429,"end_time":"2025-04-16T20:30:55.120191","exception":false,"start_time":"2025-04-16T20:30:55.087762","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"7babd2ab","cell_type":"code","source":"def integrate_with_hybrid_pipeline(run_hybrid_pipeline_func):\n \"\"\"\n Integrates the enhanced NN model with the hybrid pipeline.\n \n Parameters:\n -----------\n run_hybrid_pipeline_func: Original hybrid pipeline function\n \n Returns:\n --------\n Modified hybrid pipeline function\n \"\"\"\n def enhanced_hybrid_pipeline(\n X_valid, \n y_valid, \n test_seq_df, \n sample_submission_df, \n output_dir, \n golden_threshold=0.6, \n seed_attempts=200, \n optimal_params={'noise': 0.21, 'corr': 0.83}\n ):\n \"\"\"\n Run a hybrid pipeline with enhanced NN quality model.\n \"\"\"\n print(\"=\" * 80)\n print(\"ENHANCED HYBRID PIPELINE: GOLDEN SEEDS + ADVANCED NN PRUNING\".center(80))\n print(\"=\" * 80)\n \n status = {\n 'success': False,\n 'golden_seeds_found': 0,\n 'nn_training_success': False,\n 'best_tm_score': 0.0,\n 'error': None\n }\n \n try:\n # PHASE 1: Find Golden Seeds (same as original)\n print(\"\\nPHASE 1: Searching for golden seeds...\")\n golden_seeds, all_seeds = find_diverse_golden_seeds(\n X_valid, \n y_valid, \n golden_threshold=golden_threshold, \n attempts=seed_attempts, \n optimal_params=optimal_params\n )\n \n # Even if we don't find golden seeds, we can use the best seeds we found\n if not golden_seeds and all_seeds:\n print(\"No golden seeds found, using top seeds from search...\")\n # Sort by TM-score\n all_seeds.sort(key=lambda x: x['tm_score'], reverse=True)\n # Take top 5 seeds\n top_seeds = all_seeds[:5]\n else:\n top_seeds = golden_seeds\n \n status['golden_seeds_found'] = len(golden_seeds)\n \n # PHASE 2: Train Enhanced Quality Assessment Model\n print(\"\\nPHASE 2: Training enhanced NN quality assessment model...\")\n try:\n enhanced_quality_model = train_enhanced_quality_model(X_valid, y_valid, X_valid, y_valid)\n rule_based_model = create_rule_based_model()\n \n # Compare models\n model_comparison = evaluate_and_compare_models(\n enhanced_quality_model, \n rule_based_model, \n X_valid, \n y_valid\n )\n \n # Use the best model\n best_model_type = model_comparison['best_model']\n if best_model_type == 'neural_network':\n quality_model = enhanced_quality_model\n print(\"Using enhanced neural network model for quality assessment\")\n else:\n quality_model = rule_based_model\n print(\"Using rule-based model for quality assessment\")\n \n status['nn_training_success'] = (best_model_type == 'neural_network')\n \n except Exception as e:\n print(f\"Error training and comparing models: {str(e)}\")\n print(\"Falling back to rule-based quality assessment...\")\n quality_model = create_rule_based_model()\n \n # PHASE 3 and beyond: same as original hybrid pipeline\n # Continue with the rest of the pipeline...\n # (generate base structures, evaluate candidates, create submission)\n \n # Call the original function with our quality model\n # This is a placeholder - in a real implementation, \n # you would continue with the rest of the pipeline using the quality_model\n \n return run_hybrid_pipeline_func(\n X_valid, \n y_valid, \n test_seq_df, \n sample_submission_df, \n output_dir, \n golden_threshold=golden_threshold, \n seed_attempts=seed_attempts, \n optimal_params=optimal_params,\n quality_model=quality_model # Pass the selected model\n )\n \n except Exception as e:\n print(f\"ERROR in enhanced hybrid pipeline: {str(e)}\")\n traceback.print_exc()\n status['error'] = str(e)\n \n # Fall back to original pipeline\n print(\"Falling back to original hybrid pipeline...\")\n return run_hybrid_pipeline_func(\n X_valid, \n y_valid, \n test_seq_df, \n sample_submission_df, \n output_dir, \n golden_threshold=golden_threshold, \n seed_attempts=seed_attempts, \n optimal_params=optimal_params\n )\n \n return enhanced_hybrid_pipeline","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.133977Z","iopub.execute_input":"2025-05-08T13:21:18.134225Z","iopub.status.idle":"2025-05-08T13:21:18.146034Z","shell.execute_reply.started":"2025-05-08T13:21:18.134200Z","shell.execute_reply":"2025-05-08T13:21:18.145228Z"},"papermill":{"duration":0.027761,"end_time":"2025-04-16T20:30:55.165931","exception":false,"start_time":"2025-04-16T20:30:55.138170","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"2373c13e","cell_type":"code","source":"def phase3_integration_with_hybrid_pipeline(run_hybrid_pipeline_func):\n \"\"\"\n Integrates the enhanced Phase 3 (base structure generation) with the hybrid pipeline.\n \n Parameters:\n -----------\n run_hybrid_pipeline_func: Original hybrid pipeline function\n \n Returns:\n --------\n Modified hybrid pipeline function\n \"\"\"\n def enhanced_hybrid_pipeline(\n X_valid, \n y_valid, \n test_seq_df, \n sample_submission_df, \n output_dir, \n golden_threshold=0.6, \n seed_attempts=200, \n optimal_params={'noise': 0.21, 'corr': 0.83},\n quality_model=None\n ):\n \"\"\"\n Run a hybrid pipeline with enhanced base structure generation.\n \"\"\"\n print(\"=\" * 80)\n print(\"ENHANCED HYBRID PIPELINE WITH RNA-SPECIFIC STRUCTURE GENERATION\".center(80))\n print(\"=\" * 80)\n \n status = {\n 'success': False,\n 'golden_seeds_found': 0,\n 'nn_training_success': False,\n 'best_tm_score': 0.0,\n 'error': None\n }\n \n try:\n # PHASE 1: Find Golden Seeds (same as original)\n print(\"\\nPHASE 1: Searching for golden seeds...\")\n golden_seeds, all_seeds = find_diverse_golden_seeds(\n X_valid, \n y_valid, \n golden_threshold=golden_threshold, \n attempts=seed_attempts, \n optimal_params=optimal_params\n )\n \n # Even if we don't find golden seeds, we can use the best seeds we found\n if not golden_seeds and all_seeds:\n print(\"No golden seeds found, using top seeds from search...\")\n # Sort by TM-score\n all_seeds.sort(key=lambda x: x['tm_score'], reverse=True)\n # Take top 5 seeds\n top_seeds = all_seeds[:5]\n else:\n top_seeds = golden_seeds\n \n status['golden_seeds_found'] = len(golden_seeds)\n \n # PHASE 2: Train Quality Assessment Model (if not provided)\n if quality_model is None:\n print(\"\\nPHASE 2: Training quality assessment model...\")\n try:\n quality_model = train_enhanced_quality_model(X_valid, y_valid, X_valid, y_valid)\n status['nn_training_success'] = True\n except Exception as e:\n print(f\"Error training quality model: {str(e)}\")\n print(\"Falling back to rule-based quality assessment...\")\n quality_model = create_rule_based_model()\n else:\n print(\"\\nPHASE 2: Using provided quality model\")\n status['nn_training_success'] = hasattr(quality_model, 'model') # Check if it's a NN model\n \n # PHASE 3: Generate Base Structures with RNA-specific optimizations\n print(\"\\nPHASE 3: Generating base structures with RNA-specific optimizations...\")\n # Prepare test features\n X_test = prepare_test_features(test_seq_df)\n \n # Generate base structures using our enhanced function\n seq_to_base_structures = generate_base_structures_with_golden_seeds(\n X_test,\n test_seq_df,\n top_seeds,\n optimal_params,\n X_valid,\n y_valid\n )\n \n # PHASE 4: Generate and evaluate diverse candidates\n print(\"\\nPHASE 4: Generating diverse candidates and evaluating quality...\")\n \n seq_to_coords = {}\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n seq = row['sequence']\n seq_length = len(seq)\n \n print(f\"Processing sequence {i+1}/{len(test_seq_df)}, ID: {target_id}, length: {seq_length}\")\n \n # Get base structures for this sequence\n base_structures = seq_to_base_structures[target_id]\n \n if not base_structures:\n print(f\"No base structures found for {target_id}. Creating emergency structure.\")\n base_structures = [create_emergency_structure(seq_length)]\n \n # Extract sequence features\n seq_features = X_test[i][:seq_length]\n \n # Generate diverse candidates\n candidates = generate_diverse_structures_from_bases(\n base_structures, \n seq_length, \n quality_model,\n num_per_base=5\n )\n \n # Evaluate and select the best structures\n try:\n top_structures = evaluate_and_prune_structures(\n candidates, \n seq_features, \n quality_model, \n top_k=5\n )\n except Exception as e:\n print(f\"Error in structure evaluation: {str(e)}\")\n print(\"Falling back to basic selection...\")\n # If evaluation fails, just use the base structures\n top_structures = base_structures[:5]\n \n # If we need more structures, pad with variations\n while len(top_structures) < 5:\n idx = len(top_structures) % len(base_structures)\n variation = sample_structural_variation(\n base_structures[idx],\n noise_level=0.1,\n preserve_distance=True,\n use_global_movement=False\n )\n top_structures.append(normalize_structure(variation))\n \n # Store the final structures\n seq_to_coords[target_id] = top_structures\n \n # PHASE 5: Create Submission\n print(\"\\nPHASE 5: Creating submission file...\")\n submission_df = create_submission_dataframe(seq_to_coords, sample_submission_df)\n \n # Save submission\n enhanced_file = os.path.join(output_dir, 'submission_enhanced.csv')\n submission_df.to_csv(enhanced_file, index=False)\n print(f\"Enhanced submission saved to {enhanced_file}\")\n \n # Save as standard submission\n standard_file = os.path.join(output_dir, 'submission.csv')\n submission_df.to_csv(standard_file, index=False)\n \n # Set success\n status['success'] = True\n \n # Get best TM-score from seeds for reporting\n if top_seeds:\n status['best_tm_score'] = max(seed['tm_score'] for seed in top_seeds)\n \n return submission_df, status\n \n except Exception as e:\n print(f\"ERROR in enhanced hybrid pipeline: {str(e)}\")\n traceback.print_exc()\n status['error'] = str(e)\n \n # Fall back to original pipeline as last resort\n print(\"Falling back to original pipeline...\")\n return run_hybrid_pipeline_func(\n X_valid, \n y_valid, \n test_seq_df, \n sample_submission_df, \n output_dir, \n golden_threshold=golden_threshold, \n seed_attempts=seed_attempts, \n optimal_params=optimal_params\n )\n \n return enhanced_hybrid_pipeline","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.148908Z","iopub.execute_input":"2025-05-08T13:21:18.149099Z","iopub.status.idle":"2025-05-08T13:21:18.162289Z","shell.execute_reply.started":"2025-05-08T13:21:18.149083Z","shell.execute_reply":"2025-05-08T13:21:18.161340Z"},"papermill":{"duration":0.034001,"end_time":"2025-04-16T20:30:55.219532","exception":false,"start_time":"2025-04-16T20:30:55.185531","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"2ee5a734","cell_type":"markdown","source":"> ## Main Execution Block","metadata":{"papermill":{"duration":0.018256,"end_time":"2025-04-16T20:30:55.256313","exception":false,"start_time":"2025-04-16T20:30:55.238057","status":"completed"},"tags":[]}},{"id":"8f1c8869","cell_type":"code","source":"if __name__ == \"__main__\":\n # Execution mode selection\n use_hybrid_pipeline = True # combine golden seeds and NN pruning\n use_nn_pruning = False # Use only NN pruning\n use_reference_only = False # Use only reference-based approach\n use_adaptive_seeds = True # Use adaptive seed search\n use_enhanced_adaptive = True # Use enhanced adaptive seed search with RNA-specific targeting\n \n # Print startup banner\n print(\"=\" * 80)\n print(\"RNA 3D STRUCTURE PREDICTION PIPELINE\".center(80))\n print(\"=\" * 80)\n \n # Print selected mode\n if use_hybrid_pipeline:\n if use_enhanced_adaptive:\n mode_description = \"Enhanced Hybrid Pipeline: RNA-Targeted Adaptive Seeds + NN Pruning\"\n elif use_adaptive_seeds:\n mode_description = \"Enhanced Hybrid Pipeline: Adaptive Golden Seeds + NN Pruning\"\n else:\n mode_description = \"Hybrid Pipeline: Golden Seeds + NN Pruning\"\n elif use_nn_pruning:\n mode_description = \"Neural Network based pruning pipeline\"\n elif use_reference_only:\n mode_description = \"Reference model only\"\n else:\n mode_description = \"Standard pipeline\"\n \n print(f\"Selected mode: {mode_description}\")\n print(\"-\" * 80)\n \n try:\n # Execute the selected pipeline\n if use_hybrid_pipeline:\n start_time = time.time()\n \n print(\"Loading processed data...\")\n X_train, y_train, X_valid, y_valid = load_processed_data()\n \n print(\"\\nLoading test data...\")\n test_seq_df = pd.read_csv(os.path.join(DATA_DIR, \"test_sequences.csv\"))\n sample_submission_df = pd.read_csv(os.path.join(DATA_DIR, \"sample_submission.csv\"))\n \n if use_enhanced_adaptive:\n # Use enhanced adaptive search with RNA-specific targeting\n print(\"\\nPHASE 1: Searching for golden seeds with enhanced RNA-targeted adaptive search...\")\n golden_seeds, all_seeds = enhanced_adaptive_seed_search(\n X_valid, \n y_valid, \n initial_threshold=0.65, # Start with higher threshold\n min_threshold=0.55, # Accept down to this minimum value\n initial_attempts=100, # Initial number of attempts\n max_attempts=300, # Maximum number of attempts\n optimal_params={'noise': 0.21, 'corr': 0.83}\n )\n \n # Define a modified version that uses the enhanced seeds\n def run_hybrid_pipeline_with_enhanced_seeds(\n X_valid, \n y_valid,\n test_seq_df, \n sample_submission_df,\n output_dir,\n golden_seeds, # Pass enhanced seeds directly\n optimal_params={'noise': 0.21, 'corr': 0.83}\n ):\n \"\"\"Modified pipeline that uses RNA-specific enhanced seeds\"\"\"\n status = {\n 'success': False,\n 'golden_seeds_found': len(golden_seeds),\n 'nn_training_success': False,\n 'best_tm_score': max([s.get('tm_score', 0) for s in golden_seeds]) if golden_seeds else 0.0,\n 'error': None,\n 'category_stats': {} # Will store RNA category statistics\n }\n \n # Collect category statistics\n categories = {}\n for seed in golden_seeds:\n category = seed.get('category', 'unknown')\n if category not in categories:\n categories[category] = 0\n categories[category] += 1\n \n status['category_stats'] = categories\n \n try:\n # PHASE 2: Train Quality Assessment Model\n print(\"\\nPHASE 2: Training NN quality assessment model...\")\n try:\n quality_model = train_enhanced_quality_model(X_valid, y_valid, X_valid, y_valid)\n status['nn_training_success'] = True\n except Exception as e:\n print(f\"Error training NN model: {str(e)}\")\n print(\"Falling back to rule-based quality assessment...\")\n quality_model = create_rule_based_model()\n \n # PHASE 3: Generate Base Structures with Enhanced Seeds\n print(\"\\nPHASE 3: Generating base structures with RNA-targeted seeds...\")\n X_test = prepare_test_features(test_seq_df)\n \n # Generate predictions using each of the top seeds\n seed_predictions = []\n for i, seed_info in enumerate(golden_seeds):\n category = seed_info.get('category', 'unknown')\n print(f\"Generating predictions with seed {seed_info['seed']} (TM-score: {seed_info.get('tm_score', 0):.4f}, Category: {category})...\")\n \n # Set the random seed\n np.random.seed(seed_info['seed'])\n \n # Adjust noise level based on RNA category\n noise_level = optimal_params['noise']\n if category == 'small_RNA':\n noise_level *= 0.9 # More stable for small RNAs\n elif category == 'large_RNA':\n noise_level *= 0.7 # Lower noise for large RNAs\n \n # Create model with this seed\n model = reference_based_approach(\n X_valid, \n y_valid,\n geometric_sampling=True,\n noise_level=noise_level,\n correlation=optimal_params['corr']\n )\n \n # Generate predictions\n if model is not None:\n preds = model.predict(X_test)\n seed_predictions.append({\n 'seed': seed_info['seed'],\n 'tm_score': seed_info.get('tm_score', 0),\n 'category': category,\n 'predictions': preds\n })\n \n # Update best TM-score for status\n if seed_info.get('tm_score', 0) > status['best_tm_score']:\n status['best_tm_score'] = seed_info.get('tm_score', 0)\n else:\n print(f\"Failed to create model with seed {seed_info['seed']}\")\n \n if not seed_predictions:\n raise Exception(\"Failed to generate any predictions with enhanced seeds\")\n \n # PHASE 4: Generate and Prune Structures\n print(\"\\nPHASE 4: Generating diverse candidates with RNA-specific adjustments...\")\n \n seq_to_coords = {}\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n seq = row['sequence']\n seq_length = len(seq)\n \n print(f\"Processing sequence {i+1}/{len(test_seq_df)}, ID: {target_id}, length: {seq_length}\")\n \n # Determine RNA category for this sequence\n if seq_length < 50:\n seq_category = \"small_RNA\"\n elif seq_length < 120:\n seq_category = \"medium_RNA\"\n else:\n seq_category = \"large_RNA\"\n \n print(f\" Sequence category: {seq_category}\")\n \n # Collect base predictions from all seeds\n base_structures = []\n \n # First try predictions from seeds matching this RNA category\n category_matched = False\n for pred_info in seed_predictions:\n if pred_info.get('category', '') == seq_category:\n category_matched = True\n base_struct = pred_info['predictions'][i][:seq_length]\n \n # Extract GC content for RNA-specific processing\n seq_features = X_test[i][:seq_length]\n valid_mask = ~np.all(seq_features == 0, axis=1)\n if np.sum(valid_mask) > 0:\n features = seq_features[valid_mask]\n g_content = np.mean(features[:, 2]) if features.shape[1] > 2 else 0\n c_content = np.mean(features[:, 1]) if features.shape[1] > 1 else 0\n gc_content = g_content + c_content\n else:\n gc_content = 0.5\n \n # Apply RNA-specific post-processing\n use_global_movement = seq_category in ['large_RNA', 'medium_RNA']\n \n # Process structure with RNA-specific adjustments\n processed_struct = post_process_rna_structure(\n base_struct,\n seq,\n gc_content,\n use_global_movement=use_global_movement\n )\n \n # Normalize and validate structure\n normalized_struct = normalize_structure(processed_struct)\n if check_structure_validity(normalized_struct):\n base_structures.append(normalized_struct)\n else:\n # Try to repair structure\n repaired_struct = repair_invalid_structure(normalized_struct)\n if check_structure_validity(repaired_struct):\n base_structures.append(repaired_struct)\n print(f\" Repaired structure from seed {pred_info['seed']}\")\n \n # If no category-matched seeds, use all seeds\n if not category_matched or not base_structures:\n print(\" No category-matched seeds found, using all available seeds.\")\n for pred_info in seed_predictions:\n base_struct = pred_info['predictions'][i][:seq_length]\n normalized_struct = normalize_structure(base_struct)\n if check_structure_validity(normalized_struct):\n base_structures.append(normalized_struct)\n \n # Generate emergency structure if no valid base structures\n if not base_structures:\n print(f\" Warning: No valid structures for {target_id}. Creating emergency structure.\")\n emergency_structure = create_emergency_structure(seq_length)\n base_structures.append(emergency_structure)\n \n # Extract sequence features\n seq_features = X_test[i][:seq_length]\n \n # Generate more candidates through RNA-specific variations\n candidates = generate_diverse_structures_from_bases(\n base_structures, \n seq_length, \n quality_model,\n num_per_base=5\n )\n \n # Evaluate and prune candidates\n if status['nn_training_success']:\n print(\" Using NN model for quality assessment...\")\n try:\n top_structures = evaluate_and_prune_structures(\n candidates, \n seq_features, \n quality_model, \n top_k=5\n )\n except Exception as e:\n print(f\" Error in NN evaluation: {str(e)}\")\n print(\" Falling back to rule-based assessment...\")\n top_structures = evaluate_and_prune_rules(candidates, top_k=5)\n else:\n print(\" Using rule-based quality assessment...\")\n top_structures = evaluate_and_prune_rules(candidates, top_k=5)\n \n # Store the final structures\n seq_to_coords[target_id] = top_structures\n \n # PHASE 5: Create Submission\n print(\"\\nPHASE 5: Creating submission file...\")\n submission_df = create_submission_dataframe(seq_to_coords, sample_submission_df)\n \n # Save submission\n enhanced_file = os.path.join(output_dir, 'submission_enhanced.csv')\n submission_df.to_csv(enhanced_file, index=False)\n print(f\"Enhanced submission saved to {enhanced_file}\")\n \n # Save as standard submission\n standard_file = os.path.join(output_dir, 'submission.csv')\n submission_df.to_csv(standard_file, index=False)\n \n # Set success\n status['success'] = True\n \n return submission_df, status\n \n except Exception as e:\n print(f\"ERROR in enhanced hybrid pipeline: {str(e)}\")\n traceback.print_exc()\n status['error'] = str(e)\n return None, status\n \n # Run the enhanced pipeline with our found seeds\n submission_df, status = run_hybrid_pipeline_with_enhanced_seeds(\n X_valid, y_valid,\n test_seq_df, sample_submission_df,\n OUTPUT_DIR,\n golden_seeds,\n optimal_params={'noise': 0.21, 'corr': 0.83}\n )\n \n elif use_adaptive_seeds:\n # Use adaptive search for golden seeds\n print(\"\\nPHASE 1: Searching for golden seeds with adaptive threshold...\")\n golden_seeds, all_seeds = adaptive_seed_search(\n X_valid, \n y_valid, \n initial_threshold=0.65, # Start with higher threshold\n min_threshold=0.55, # Accept down to this minimum value\n initial_attempts=100, # Initial number of attempts\n max_attempts=300, # Maximum number of attempts\n optimal_params={'noise': 0.21, 'corr': 0.83}\n )\n \n # Modify run_hybrid_pipeline to use our found seeds\n # Define a modified version of run_hybrid_pipeline that uses predefined seeds\n def run_hybrid_pipeline_with_seeds(\n X_valid, \n y_valid,\n test_seq_df, \n sample_submission_df,\n output_dir,\n golden_seeds, # Pass seeds directly\n optimal_params={'noise': 0.21, 'corr': 0.83}\n ):\n # This is a modified version that skips the seed search\n # and uses the seeds provided by adaptive search\n \n status = {\n 'success': False,\n 'golden_seeds_found': len(golden_seeds),\n 'nn_training_success': False,\n 'best_tm_score': max([s['tm_score'] for s in golden_seeds]) if golden_seeds else 0.0,\n 'error': None\n }\n \n try:\n # PHASE 2: Train Quality Assessment Model\n print(\"\\nPHASE 2: Training NN quality assessment model...\")\n try:\n quality_model = train_enhanced_quality_model(X_valid, y_valid, X_valid, y_valid)\n status['nn_training_success'] = True\n except Exception as e:\n print(f\"Error training NN model: {str(e)}\")\n print(\"Falling back to rule-based quality assessment...\")\n quality_model = create_rule_based_model()\n \n # PHASE 3: Generate Base Structures with Golden Seeds\n print(\"\\nPHASE 3: Generating base structures with golden seeds...\")\n X_test = prepare_test_features(test_seq_df)\n \n # Generate predictions using each of the top seeds\n seed_predictions = []\n for i, seed_info in enumerate(golden_seeds):\n print(f\"Generating predictions with seed {seed_info['seed']} (TM-score: {seed_info['tm_score']:.4f})...\")\n \n # Set the random seed\n np.random.seed(seed_info['seed'])\n \n # Create model with this seed\n model = reference_based_approach(\n X_valid, \n y_valid,\n geometric_sampling=True,\n noise_level=optimal_params['noise'],\n correlation=optimal_params['corr']\n )\n \n # Generate predictions\n if model is not None:\n preds = model.predict(X_test)\n seed_predictions.append({\n 'seed': seed_info['seed'],\n 'tm_score': seed_info['tm_score'],\n 'predictions': preds\n })\n \n # Update best TM-score for status\n if seed_info['tm_score'] > status['best_tm_score']:\n status['best_tm_score'] = seed_info['tm_score']\n else:\n print(f\"Failed to create model with seed {seed_info['seed']}\")\n \n if not seed_predictions:\n raise Exception(\"Failed to generate any predictions with golden seeds\")\n \n # PHASE 4: Generate and Prune Structures\n print(\"\\nPHASE 4: Generating diverse candidates and using NN pruning...\")\n \n seq_to_coords = {}\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n seq = row['sequence']\n seq_length = len(seq)\n \n print(f\"Processing sequence {i+1}/{len(test_seq_df)}, ID: {target_id}, length: {seq_length}\")\n \n # Collect base predictions from all seeds for this sequence\n base_structures = []\n for pred_info in seed_predictions:\n base_struct = pred_info['predictions'][i][:seq_length]\n base_structures.append(normalize_structure(base_struct))\n \n # Extract sequence features\n seq_features = X_test[i][:seq_length]\n \n # Generate more candidates through controlled variations\n candidates = generate_diverse_candidates(base_structures, seq_length, num_per_base=5)\n \n # Evaluate and prune candidates\n if status['nn_training_success']:\n print(\"Using NN model for quality assessment...\")\n try:\n top_structures = evaluate_and_prune_structures(\n candidates, \n seq_features, \n quality_model, \n top_k=5\n )\n except Exception as e:\n print(f\"Error in NN evaluation: {str(e)}\")\n print(\"Falling back to rule-based assessment...\")\n top_structures = evaluate_and_prune_rules(candidates, top_k=5)\n else:\n print(\"Using rule-based quality assessment...\")\n top_structures = evaluate_and_prune_rules(candidates, top_k=5)\n \n # Store the final structures\n seq_to_coords[target_id] = top_structures\n \n # PHASE 5: Create Submission\n print(\"\\nPHASE 5: Creating submission file...\")\n submission_df = create_submission_dataframe(seq_to_coords, sample_submission_df)\n \n # Save submission\n hybrid_file = os.path.join(output_dir, 'submission_hybrid.csv')\n submission_df.to_csv(hybrid_file, index=False)\n print(f\"Hybrid submission saved to {hybrid_file}\")\n \n # Save as standard submission\n standard_file = os.path.join(output_dir, 'submission.csv')\n submission_df.to_csv(standard_file, index=False)\n \n # Set success\n status['success'] = True\n \n return submission_df, status\n \n except Exception as e:\n print(f\"ERROR in hybrid pipeline: {str(e)}\")\n traceback.print_exc()\n status['error'] = str(e)\n return None, status\n \n # Run the modified pipeline with our found seeds\n submission_df, status = run_hybrid_pipeline_with_seeds(\n X_valid, y_valid,\n test_seq_df, sample_submission_df,\n OUTPUT_DIR,\n golden_seeds,\n optimal_params={'noise': 0.21, 'corr': 0.83}\n )\n else:\n # Run the standard hybrid pipeline\n submission_df, status = run_hybrid_pipeline(\n X_valid, y_valid,\n test_seq_df, sample_submission_df,\n OUTPUT_DIR,\n golden_threshold=0.6,\n seed_attempts=100\n )\n \n # Calculate total runtime\n runtime = time.time() - start_time\n hours, remainder = divmod(runtime, 3600)\n minutes, seconds = divmod(remainder, 60)\n \n # Display results summary\n print(\"\\n\" + \"=\" * 80)\n if use_enhanced_adaptive:\n print(\"RNA-TARGETED ENHANCED HYBRID PIPELINE RESULTS SUMMARY\".center(80))\n elif use_adaptive_seeds:\n print(\"ENHANCED HYBRID PIPELINE RESULTS SUMMARY\".center(80))\n else:\n print(\"HYBRID PIPELINE RESULTS SUMMARY\".center(80))\n print(\"=\" * 80)\n print(f\"Total runtime: {int(hours)}h {int(minutes)}m {int(seconds)}s\")\n \n if status['success']:\n print(\"\\nPIPELINE STATISTICS:\")\n print(f\" - Golden seeds found: {status['golden_seeds_found']}\")\n print(f\" - NN training success: {status['nn_training_success']}\")\n print(f\" - Best TM-score: {status['best_tm_score']:.4f}\")\n \n # Additional statistics for enhanced pipeline\n if use_enhanced_adaptive and 'category_stats' in status:\n print(\"\\nRNA CATEGORY STATISTICS:\")\n for category, stats in status.get('category_stats', {}).items():\n print(f\" - {category}: {stats}\")\n else:\n print(f\"\\nPipeline failed with error: {status['error']}\")\n \n # Display output file information\n print(\"\\nOUTPUT FILES:\")\n if use_enhanced_adaptive:\n submission_file = os.path.join(OUTPUT_DIR, 'submission_enhanced.csv')\n if os.path.exists(submission_file):\n try:\n file_size = os.path.getsize(submission_file)\n print(f\" - Enhanced submission: {submission_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - Enhanced submission: {submission_file}\")\n else:\n submission_file = os.path.join(OUTPUT_DIR, 'submission_hybrid.csv')\n if os.path.exists(submission_file):\n try:\n file_size = os.path.getsize(submission_file)\n print(f\" - Hybrid submission: {submission_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - Hybrid submission: {submission_file}\")\n \n standard_file = os.path.join(OUTPUT_DIR, 'submission.csv')\n if os.path.exists(standard_file):\n try:\n file_size = os.path.getsize(standard_file)\n print(f\" - Standard submission: {standard_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - Standard submission: {standard_file}\")\n \n print(\"=\" * 80)\n \n elif use_nn_pruning:\n start_time = time.time()\n \n print(\"Loading processed data...\")\n X_train, y_train, X_valid, y_valid = load_processed_data()\n \n print(\"\\nLoading test data...\")\n test_seq_df = pd.read_csv(os.path.join(DATA_DIR, \"test_sequences.csv\"))\n sample_submission_df = pd.read_csv(os.path.join(DATA_DIR, \"sample_submission.csv\"))\n \n # Train quality model\n print(\"\\nTraining quality assessment model...\")\n quality_model = train_enhanced_quality_model(X_valid, y_valid, X_valid, y_valid)\n \n # Create reference model with default parameters\n print(\"\\nCreating reference model...\")\n reference_model = reference_based_approach(\n X_valid, \n y_valid,\n geometric_sampling=True,\n noise_level=0.21,\n correlation=0.83\n )\n \n # Generate submission using NN pruning only\n submission_df = generate_nn_pruned_submission(\n reference_model,\n quality_model,\n test_seq_df,\n sample_submission_df\n )\n \n # Calculate total runtime\n runtime = time.time() - start_time\n hours, remainder = divmod(runtime, 3600)\n minutes, seconds = divmod(remainder, 60)\n \n print(\"\\n\" + \"=\" * 80)\n print(\"NN PRUNING PIPELINE RESULTS\".center(80))\n print(\"=\" * 80)\n print(f\"Total runtime: {int(hours)}h {int(minutes)}m {int(seconds)}s\")\n \n # Display output file information\n print(\"\\nOUTPUT FILES:\")\n submission_file = os.path.join(OUTPUT_DIR, 'submission_nn_pruned.csv')\n if os.path.exists(submission_file):\n try:\n file_size = os.path.getsize(submission_file)\n print(f\" - NN pruned submission: {submission_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - NN pruned submission: {submission_file}\")\n \n standard_file = os.path.join(OUTPUT_DIR, 'submission.csv')\n if os.path.exists(standard_file):\n try:\n file_size = os.path.getsize(standard_file)\n print(f\" - Standard submission: {standard_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - Standard submission: {standard_file}\")\n \n print(\"=\" * 80)\n \n elif use_reference_only:\n start_time = time.time()\n \n print(\"Loading processed data...\")\n X_train, y_train, X_valid, y_valid = load_processed_data()\n \n print(\"\\nLoading test data...\")\n test_seq_df = pd.read_csv(os.path.join(DATA_DIR, \"test_sequences.csv\"))\n sample_submission_df = pd.read_csv(os.path.join(DATA_DIR, \"sample_submission.csv\"))\n \n # Create optimized reference model\n print(\"\\nCreating and evaluating reference model...\")\n reference_model = reference_based_approach(\n X_valid, \n y_valid,\n geometric_sampling=True,\n noise_level=0.21,\n correlation=0.83\n )\n \n metrics = evaluate_model(reference_model, X_valid, y_valid)\n tm_score = metrics['avg_tm_score']\n print(f\"Reference model TM-score: {tm_score:.4f}\")\n \n # Prepare test sequences\n X_test = prepare_test_features(test_seq_df)\n \n # Generate predictions\n print(\"\\nGenerating predictions...\")\n predictions = reference_model.predict(X_test)\n \n # Create submission dataframe\n print(\"\\nCreating submission dataframe...\")\n submission_df = sample_submission_df.copy()\n \n seq_to_coords = {}\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n seq_length = len(row['sequence'])\n \n # Normalize and process structure\n struct = normalize_structure(predictions[i][:seq_length])\n \n # Create 5 copies with small variations\n structures = [struct]\n for j in range(4):\n variation = sample_structural_variation(\n struct,\n noise_level=0.05,\n preserve_distance=True,\n correlation=0.9\n )\n structures.append(normalize_structure(variation))\n \n seq_to_coords[target_id] = structures\n \n # Fill the dataframe\n for i, row in submission_df.iterrows():\n if i % 1000 == 0:\n print(f\"Processing row {i}/{len(submission_df)}\")\n \n id_parts = row['ID'].split('_')\n seq_id = id_parts[0]\n residue_idx = int(id_parts[1]) - 1\n \n if seq_id in seq_to_coords:\n structures = seq_to_coords[seq_id]\n if residue_idx < len(structures[0]):\n for struct_idx in range(5):\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[struct_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[struct_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[struct_idx][residue_idx][2]\n \n # Save submission\n reference_file = os.path.join(OUTPUT_DIR, 'submission_reference.csv')\n submission_df.to_csv(reference_file, index=False)\n print(f\"Reference submission saved to {reference_file}\")\n \n # Save as standard submission\n standard_file = os.path.join(OUTPUT_DIR, 'submission.csv')\n submission_df.to_csv(standard_file, index=False)\n \n # Calculate total runtime\n runtime = time.time() - start_time\n hours, remainder = divmod(runtime, 3600)\n minutes, seconds = divmod(remainder, 60)\n \n print(\"\\n\" + \"=\" * 80)\n print(\"REFERENCE MODEL RESULTS\".center(80))\n print(\"=\" * 80)\n print(f\"Total runtime: {int(hours)}h {int(minutes)}m {int(seconds)}s\")\n \n # Display output file information\n print(\"\\nOUTPUT FILES:\")\n if os.path.exists(reference_file):\n try:\n file_size = os.path.getsize(reference_file)\n print(f\" - Reference submission: {reference_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - Reference submission: {reference_file}\")\n\n if os.path.exists(standard_file):\n try:\n file_size = os.path.getsize(standard_file)\n print(f\" - Standard submission: {standard_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - Standard submission: {standard_file}\")\n \n print(\"=\" * 80)\n \n else:\n # Standard pipeline - if user disabled all options\n print(\"No pipeline mode selected. Please set one of the pipeline flags to True.\")\n print(\"Available options:\")\n print(\" - use_hybrid_pipeline: Combined golden seeds and NN pruning\")\n print(\" - use_nn_pruning: Neural Network based pruning only\")\n print(\" - use_reference_only: Use only reference model approach\")\n print(\" - use_adaptive_seeds: Use adaptive seed search (with hybrid pipeline)\")\n print(\" - use_enhanced_adaptive: Use enhanced RNA-specific adaptive seed search\")\n \n print(\"\\nProcess completed.\")\n \n except Exception as e:\n print(\"\\n\" + \"=\" * 80)\n print(\"ERROR IN MAIN EXECUTION\".center(80))\n print(\"=\" * 80)\n print(f\"Critical error: {str(e)}\")\n traceback.print_exc()\n print(\"=\" * 80)","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:21:18.163392Z","iopub.execute_input":"2025-05-08T13:21:18.163657Z","iopub.status.idle":"2025-05-08T13:36:30.412980Z","shell.execute_reply.started":"2025-05-08T13:21:18.163633Z","shell.execute_reply":"2025-05-08T13:36:30.412257Z"},"papermill":{"duration":125.527692,"end_time":"2025-04-16T20:33:00.802531","exception":false,"start_time":"2025-04-16T20:30:55.274839","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"c69476c9","cell_type":"code","source":"submission_df = pd.read_csv('/kaggle/working/submission.csv')\nprint(\"Overview of the DataFrame:\")\nprint(submission_df.shape) # Print the shape (rows, columns)\nprint(submission_df.head()) # Display the first 5 rows","metadata":{"execution":{"iopub.status.busy":"2025-05-08T13:36:30.413916Z","iopub.execute_input":"2025-05-08T13:36:30.414282Z","iopub.status.idle":"2025-05-08T13:36:30.436558Z","shell.execute_reply.started":"2025-05-08T13:36:30.414250Z","shell.execute_reply":"2025-05-08T13:36:30.435850Z"},"papermill":{"duration":0.059942,"end_time":"2025-04-16T20:33:00.901748","exception":false,"start_time":"2025-04-16T20:33:00.841806","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"e3557841-9ec1-4df6-be9b-84969a70d409","cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}