{"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,"sourceType":"competition"}],"dockerImageVersionId":30919,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true},"papermill":{"default_parameters":{},"duration":144.373007,"end_time":"2025-03-26T03:44:56.636459","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2025-03-26T03:42:32.263452","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# ","metadata":{"papermill":{"duration":0.012176,"end_time":"2025-03-26T03:42:35.898136","exception":false,"start_time":"2025-03-26T03:42:35.88596","status":"completed"},"tags":[]}},{"cell_type":"code","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.\n\n## Library Imports 📚🔧","metadata":{"papermill":{"duration":0.012176,"end_time":"2025-03-26T03:42:35.898136","exception":false,"start_time":"2025-03-26T03:42:35.88596","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Standard Library Imports\nimport os\nimport time\nimport gc\nimport traceback\nfrom collections import Counter\nimport warnings\nimport hashlib\nimport random\nimport datetime\nimport json\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-04-05T16:07:41.557542Z","iopub.execute_input":"2025-04-05T16:07:41.55809Z","iopub.status.idle":"2025-04-05T16:07:59.982736Z","shell.execute_reply.started":"2025-04-05T16:07:41.55802Z","shell.execute_reply":"2025-04-05T16:07:59.981991Z"},"papermill":{"duration":21.560049,"end_time":"2025-03-26T03:42:57.488491","exception":false,"start_time":"2025-03-26T03:42:35.928442","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🧬 RNA 3D Structure Prediction and Analysis Pipeline 🔬","metadata":{"papermill":{"duration":0.009708,"end_time":"2025-03-26T03:42:57.508811","exception":false,"start_time":"2025-03-26T03:42:57.499103","status":"completed"},"tags":[]}},{"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-04-05T16:07:59.983731Z","iopub.execute_input":"2025-04-05T16:07:59.984199Z","iopub.status.idle":"2025-04-05T16:08:00.541001Z","shell.execute_reply.started":"2025-04-05T16:07:59.984177Z","shell.execute_reply":"2025-04-05T16:08:00.540022Z"},"papermill":{"duration":0.867643,"end_time":"2025-03-26T03:42:58.386472","exception":false,"start_time":"2025-03-26T03:42:57.518829","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Directory Explorer & CSV Verification for RNA3D 🗂️🔬","metadata":{"papermill":{"duration":0.010576,"end_time":"2025-03-26T03:42:58.408104","exception":false,"start_time":"2025-03-26T03:42:58.397528","status":"completed"},"tags":[]}},{"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-04-05T16:08:00.542727Z","iopub.execute_input":"2025-04-05T16:08:00.543129Z","iopub.status.idle":"2025-04-05T16:08:00.615756Z","shell.execute_reply.started":"2025-04-05T16:08:00.543102Z","shell.execute_reply":"2025-04-05T16:08:00.614943Z"},"papermill":{"duration":0.094882,"end_time":"2025-03-26T03:42:58.513093","exception":false,"start_time":"2025-03-26T03:42:58.418211","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## RNA3D Data Checker 🔍🧬","metadata":{"papermill":{"duration":0.010053,"end_time":"2025-03-26T03:42:58.533702","exception":false,"start_time":"2025-03-26T03:42:58.523649","status":"completed"},"tags":[]}},{"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.csv\" 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-04-05T16:08:00.616969Z","iopub.execute_input":"2025-04-05T16:08:00.617178Z","iopub.status.idle":"2025-04-05T16:08:00.890912Z","shell.execute_reply.started":"2025-04-05T16:08:00.617159Z","shell.execute_reply":"2025-04-05T16:08:00.890025Z"},"papermill":{"duration":0.270583,"end_time":"2025-03-26T03:42:58.814635","exception":false,"start_time":"2025-03-26T03:42:58.544052","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Integrated RNA3D Sequence and Structure Analyzer 🔬🧬","metadata":{"papermill":{"duration":0.010382,"end_time":"2025-03-26T03:42:58.836058","exception":false,"start_time":"2025-03-26T03:42:58.825676","status":"completed"},"tags":[]}},{"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-04-05T16:08:00.892095Z","iopub.execute_input":"2025-04-05T16:08:00.892468Z","iopub.status.idle":"2025-04-05T16:08:03.673011Z","shell.execute_reply.started":"2025-04-05T16:08:00.892431Z","shell.execute_reply":"2025-04-05T16:08:03.671964Z"},"papermill":{"duration":2.579715,"end_time":"2025-03-26T03:43:01.426154","exception":false,"start_time":"2025-03-26T03:42:58.846439","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Data Preparation for RNA 3D Structure Prediction 🧬🔍","metadata":{"papermill":{"duration":0.013504,"end_time":"2025-03-26T03:43:01.454026","exception":false,"start_time":"2025-03-26T03:43:01.440522","status":"completed"},"tags":[]}},{"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-04-05T16:08:03.673871Z","iopub.execute_input":"2025-04-05T16:08:03.674199Z","iopub.status.idle":"2025-04-05T16:08:09.685632Z","shell.execute_reply.started":"2025-04-05T16:08:03.674168Z","shell.execute_reply":"2025-04-05T16:08:09.684623Z"},"papermill":{"duration":5.431349,"end_time":"2025-03-26T03:43:06.89879","exception":false,"start_time":"2025-03-26T03:43:01.467441","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Heatmap Viewer for RNA Sequences 🔥🧬","metadata":{"papermill":{"duration":0.013025,"end_time":"2025-03-26T03:43:06.925988","exception":false,"start_time":"2025-03-26T03:43:06.912963","status":"completed"},"tags":[]}},{"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-04-05T16:08:09.686616Z","iopub.execute_input":"2025-04-05T16:08:09.686961Z","iopub.status.idle":"2025-04-05T16:08:10.613457Z","shell.execute_reply.started":"2025-04-05T16:08:09.686928Z","shell.execute_reply":"2025-04-05T16:08:10.612652Z"},"papermill":{"duration":0.854076,"end_time":"2025-03-26T03:43:07.793328","exception":false,"start_time":"2025-03-26T03:43:06.939252","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🧬 RNA 3D Structure Prediction with Neural Network-Based Conformational Selection 🔬","metadata":{"papermill":{"duration":0.018177,"end_time":"2025-03-26T03:43:07.831159","exception":false,"start_time":"2025-03-26T03:43:07.812982","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"> ## Utility and Preprocessing Functions","metadata":{"papermill":{"duration":0.018692,"end_time":"2025-03-26T03:43:07.869111","exception":false,"start_time":"2025-03-26T03:43:07.850419","status":"completed"},"tags":[]}},{"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\nclass StructureWrapper:\n \"\"\"\n Wrapper for RNA structure arrays that provides both quality attributes\n and compatibility with NumPy array operations.\n \"\"\"\n def __init__(self, structure, quality_score=0.5):\n self.structure = structure\n self.quality = {'quality_score': quality_score}\n # Store shape from the underlying structure for numpy compatibility\n self.shape = structure.shape if hasattr(structure, 'shape') else None\n \n def __getitem__(self, idx):\n return self.structure[idx]\n \n def __len__(self):\n return len(self.structure)\n \n # Implement arithmetic operators for numpy compatibility\n def __sub__(self, other):\n \"\"\"Implement subtraction between structures\"\"\"\n if isinstance(other, StructureWrapper):\n # Subtract the underlying structures\n return StructureWrapper(self.structure - other.structure)\n else:\n # Subtract a scalar or numpy array directly\n return StructureWrapper(self.structure - other)\n \n def __add__(self, other):\n \"\"\"Implement addition between structures\"\"\"\n if isinstance(other, StructureWrapper):\n return StructureWrapper(self.structure + other.structure)\n else:\n return StructureWrapper(self.structure + other)\n \n def __mul__(self, other):\n \"\"\"Implement multiplication between structures\"\"\"\n if isinstance(other, StructureWrapper):\n return StructureWrapper(self.structure * other.structure)\n else:\n return StructureWrapper(self.structure * other)\n \n def __truediv__(self, other):\n \"\"\"Implement division between structures\"\"\"\n if isinstance(other, StructureWrapper):\n return StructureWrapper(self.structure / other.structure)\n else:\n return StructureWrapper(self.structure / other)\n \n def __neg__(self):\n \"\"\"Implement negation\"\"\"\n return StructureWrapper(-self.structure)\n \n def __abs__(self):\n \"\"\"Implement absolute value\"\"\"\n return StructureWrapper(abs(self.structure))\n \n # Implement reverse operations (for scalar op structure)\n def __radd__(self, other):\n return StructureWrapper(other + self.structure)\n \n def __rsub__(self, other):\n return StructureWrapper(other - self.structure)\n \n def __rmul__(self, other):\n return StructureWrapper(other * self.structure)\n \n def __rtruediv__(self, other):\n return StructureWrapper(other / self.structure)\n \n # Implement comparison operators\n def __eq__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure == other.structure\n else:\n return self.structure == other\n \n def __lt__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure < other.structure\n else:\n return self.structure < other\n \n def __gt__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure > other.structure\n else:\n return self.structure > other\n \n def __le__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure <= other.structure\n else:\n return self.structure <= other\n \n def __ge__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure >= other.structure\n else:\n return self.structure >= other\n \n # Implement numpy compatibility methods\n def __array__(self):\n \"\"\"Allow numpy to automatically convert to array when needed\"\"\"\n import numpy as np\n return np.array(self.structure)\n \n def sum(self, *args, **kwargs):\n \"\"\"Implement sum method for numpy compatibility\"\"\"\n return self.structure.sum(*args, **kwargs)\n \n def mean(self, *args, **kwargs):\n \"\"\"Implement mean method for numpy compatibility\"\"\"\n return self.structure.mean(*args, **kwargs)\n \n def max(self, *args, **kwargs):\n \"\"\"Implement max method for numpy compatibility\"\"\"\n return self.structure.max(*args, **kwargs)\n \n def min(self, *args, **kwargs):\n \"\"\"Implement min method for numpy compatibility\"\"\"\n return self.structure.min(*args, **kwargs)\n \n def reshape(self, *args, **kwargs):\n \"\"\"Implement reshape method for numpy compatibility\"\"\"\n reshaped = self.structure.reshape(*args, **kwargs)\n return StructureWrapper(reshaped, self.quality.get('quality_score', 0.5))\n \n def transpose(self, *args, **kwargs):\n \"\"\"Implement transpose method for numpy compatibility\"\"\"\n transposed = self.structure.transpose(*args, **kwargs)\n return StructureWrapper(transposed, self.quality.get('quality_score', 0.5))\n \n # String representation\n def __repr__(self):\n return f\"StructureWrapper(shape={self.shape}, quality_score={self.quality.get('quality_score', 0.5)})\"\n\nclass ParameterOptimizer:\n \"\"\"\n Meta-learning system for continuous parameter optimization\n based on historical results.\n \"\"\"\n \n def __init__(self, history_file=None):\n \"\"\"\n Initializes the optimizer, optionally loading previous history.\n \n Parameters:\n -----------\n history_file: str, optional\n Path to a file containing the parameter and result history\n \"\"\"\n self.history = []\n if history_file and os.path.exists(history_file):\n self.load_history(history_file)\n \n # Parameter bounds\n self.param_bounds = {\n 'divisor_mean': (3.0, 4.5),\n 'divisor_std': (0.5, 1.5),\n 'noise_base': (0.01, 0.3),\n 'correlation': (0.7, 0.95)\n }\n \n def load_history(self, filename):\n \"\"\"Loads previous parameter and result history\"\"\"\n try:\n with open(filename, 'r') as f:\n self.history = json.load(f)\n except Exception as e:\n print(f\"Error loading history: {str(e)}\")\n \n def save_history(self, filename):\n \"\"\"Saves the current history to a file\"\"\"\n with open(filename, 'w') as f:\n json.dump(self.history, f)\n \n def record_result(self, params, size_category, mode, quality_score):\n \"\"\"\n Records a new result into the history\n \n Parameters:\n -----------\n params: dict\n Parameters used\n size_category: str\n Size category ('small', 'medium', 'large')\n mode: str\n Mode used ('adaptive' or 'fixed')\n quality_score: float\n Quality score obtained\n \"\"\"\n self.history.append({\n 'params': params,\n 'size_category': size_category,\n 'mode': mode,\n 'quality_score': quality_score,\n 'timestamp': datetime.datetime.now().isoformat()\n })\n \n def suggest_parameters(self, size_category, mode):\n \"\"\"\n Suggests optimized parameters based on the history\n for a given size category and mode\n \n Parameters:\n -----------\n size_category: str\n Size category ('small', 'medium', 'large')\n mode: str\n Operation mode ('adaptive' or 'fixed')\n \n Returns:\n --------\n dict: Suggested parameters\n \"\"\"\n # Filter history for the specified category and mode\n relevant_history = [\n entry for entry in self.history \n if entry['size_category'] == size_category and entry['mode'] == mode\n ]\n \n if len(relevant_history) < 5:\n # Not enough history, use default values\n return self._get_default_params(size_category, mode)\n \n # Sort by quality score, from best to worst\n relevant_history.sort(key=lambda x: x['quality_score'], reverse=True)\n \n # Extract parameters from the top N results\n top_n = min(5, len(relevant_history))\n top_params = [entry['params'] for entry in relevant_history[:top_n]]\n \n # Compute weighted average of parameters\n weights = [0.4, 0.25, 0.15, 0.1, 0.1][:top_n] # Weights for top N results\n \n suggested_params = {}\n for param in self.param_bounds.keys():\n if all(param in p for p in top_params):\n weighted_sum = sum(w * p[param] for w, p in zip(weights, top_params))\n suggested_params[param] = weighted_sum / sum(weights[:top_n])\n \n # Ensure suggested parameters are within bounds\n for param, (min_val, max_val) in self.param_bounds.items():\n if param in suggested_params:\n suggested_params[param] = max(min_val, min(suggested_params[param], max_val))\n \n return suggested_params\n \n def _get_default_params(self, size_category, mode):\n \"\"\"Returns default parameters when historical data is insufficient\"\"\"\n # Default values for different size categories and modes\n defaults = {\n 'small': {\n 'adaptive': {'divisor_mean': 3.6, 'divisor_std': 0.9, 'noise_base': 0.15},\n 'fixed': {'noise_base': 0.12, 'correlation': 0.85}\n },\n 'medium': {\n 'adaptive': {'divisor_mean': 3.8, 'divisor_std': 1.0, 'noise_base': 0.1},\n 'fixed': {'noise_base': 0.08, 'correlation': 0.85}\n },\n 'large': {\n 'adaptive': {'divisor_mean': 4.0, 'divisor_std': 1.1, 'noise_base': 0.05},\n 'fixed': {'noise_base': 0.04, 'correlation': 0.9}\n }\n }\n \n return defaults.get(size_category, {}).get(mode, {})\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-04-05T16:08:10.615927Z","iopub.execute_input":"2025-04-05T16:08:10.616148Z","iopub.status.idle":"2025-04-05T16:08:10.646961Z","shell.execute_reply.started":"2025-04-05T16:08:10.616129Z","shell.execute_reply":"2025-04-05T16:08:10.646056Z"},"papermill":{"duration":0.033629,"end_time":"2025-03-26T03:43:07.921623","exception":false,"start_time":"2025-03-26T03:43:07.887994","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"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\n# Auxiliary function to calculate the dihedral angle (in degrees)\ndef calculate_dihedral(p0, p1, p2, p3):\n \"\"\"\n Calculates the dihedral angle (in degrees) defined by the points p0, p1, p2, and p3.\n \"\"\"\n b0 = p1 - p0\n b1 = p2 - p1\n b2 = p3 - p2\n\n # Normalize b1 so its length does not influence the calculation\n b1 /= np.linalg.norm(b1) + 1e-8\n\n # Normal vectors to the planes formed by (p0,p1,p2) and (p1,p2,p3)\n v = b0 - np.dot(b0, b1) * b1\n w = b2 - np.dot(b2, b1) * b1\n\n x = np.dot(v, w)\n y = np.dot(np.cross(b1, v), w)\n angle = np.degrees(np.arctan2(y, x))\n return angle\n\n# Auxiliary function to generate a rotation matrix\ndef get_rotation_matrix(axis, theta):\n \"\"\"\n Returns the 3x3 rotation matrix for a rotation of theta radians around the given 'axis'.\n \"\"\"\n a = np.cos(theta / 2.0)\n b, c, d = -axis * np.sin(theta / 2.0)\n aa, bb, cc, dd = a*a, b*b, c*c, d*d\n bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d\n return np.array([\n [aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],\n [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],\n [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]\n ])\n\n# New function to refine the RNA backbone with dihedral angle adjustment\ndef refine_rna_backbone_with_dihedrals(structure, ideal_dihedral=180.0):\n \"\"\"\n Refines the geometry of the RNA backbone by adjusting the dihedral angles to an ideal value.\n\n Parameters:\n structure: np.array of shape (seq_length, 3) containing the coordinates.\n ideal_dihedral: Ideal angle (in degrees) for the backbone segments (e.g., 180°).\n\n Returns:\n np.array with the refined structure.\n \"\"\"\n refined = structure.copy()\n n = len(refined)\n if n < 4:\n return refined # There are no dihedral angles to correct\n\n for i in range(n - 3):\n p0, p1, p2, p3 = refined[i], refined[i+1], refined[i+2], refined[i+3]\n current_angle = calculate_dihedral(p0, p1, p2, p3)\n # Compute the difference (in radians) between the ideal angle and the current one\n angle_diff = np.radians(ideal_dihedral - current_angle)\n \n # Define the rotation axis as the direction of the segment (p3 - p2)\n axis = p3 - p2\n norm_axis = np.linalg.norm(axis)\n if norm_axis < 1e-6:\n continue\n axis /= norm_axis\n \n # Get the rotation matrix for the correction angle\n R = get_rotation_matrix(axis, angle_diff)\n \n # Apply the rotation to all points starting from p3\n for j in range(i+3, n):\n vec = refined[j] - p2\n refined[j] = p2 + np.dot(vec, R.T)\n \n return refined","metadata":{"execution":{"iopub.status.busy":"2025-04-05T16:08:10.648965Z","iopub.execute_input":"2025-04-05T16:08:10.649221Z","iopub.status.idle":"2025-04-05T16:08:10.675591Z","shell.execute_reply.started":"2025-04-05T16:08:10.649201Z","shell.execute_reply":"2025-04-05T16:08:10.674943Z"},"papermill":{"duration":0.039629,"end_time":"2025-03-26T03:43:07.980725","exception":false,"start_time":"2025-03-26T03:43:07.941096","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"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-04-05T16:08:10.676487Z","iopub.execute_input":"2025-04-05T16:08:10.676778Z","iopub.status.idle":"2025-04-05T16:08:10.697462Z","shell.execute_reply.started":"2025-04-05T16:08:10.67675Z","shell.execute_reply":"2025-04-05T16:08:10.69688Z"},"papermill":{"duration":0.031523,"end_time":"2025-03-26T03:43:08.03412","exception":false,"start_time":"2025-03-26T03:43:08.002597","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"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-04-05T16:08:10.698324Z","iopub.execute_input":"2025-04-05T16:08:10.698612Z","iopub.status.idle":"2025-04-05T16:08:10.719982Z","shell.execute_reply.started":"2025-04-05T16:08:10.698584Z","shell.execute_reply":"2025-04-05T16:08:10.719374Z"},"papermill":{"duration":0.031127,"end_time":"2025-03-26T03:43:08.084534","exception":false,"start_time":"2025-03-26T03:43:08.053407","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"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-04-05T16:08:10.720731Z","iopub.execute_input":"2025-04-05T16:08:10.720958Z","iopub.status.idle":"2025-04-05T16:08:10.743724Z","shell.execute_reply.started":"2025-04-05T16:08:10.720938Z","shell.execute_reply":"2025-04-05T16:08:10.743093Z"},"papermill":{"duration":0.029709,"end_time":"2025-03-26T03:43:08.132009","exception":false,"start_time":"2025-03-26T03:43:08.1023","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## RNA-specific functions","metadata":{"papermill":{"duration":0.017717,"end_time":"2025-03-26T03:43:08.167465","exception":false,"start_time":"2025-03-26T03:43:08.149748","status":"completed"},"tags":[]}},{"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-04-05T16:08:10.744473Z","iopub.execute_input":"2025-04-05T16:08:10.744658Z","iopub.status.idle":"2025-04-05T16:08:10.765855Z","shell.execute_reply.started":"2025-04-05T16:08:10.744642Z","shell.execute_reply":"2025-04-05T16:08:10.765219Z"},"papermill":{"duration":0.031324,"end_time":"2025-03-26T03:43:08.217678","exception":false,"start_time":"2025-03-26T03:43:08.186354","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Phase 1: Golden Seeds","metadata":{"papermill":{"duration":0.017715,"end_time":"2025-03-26T03:43:08.25309","exception":false,"start_time":"2025-03-26T03:43:08.235375","status":"completed"},"tags":[]}},{"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 find_diverse_golden_seeds(\n X_valid, \n y_valid, \n golden_threshold=0.6, \n attempts=200, \n optimal_params={'noise': 0.21, 'corr': 0.83},\n diversity_threshold=0.15,\n max_seeds=10\n):\n \"\"\"\n Searches for \"golden\" seeds that produce good results, ensuring diversity\n and controlling overfitting.\n \n Parameters:\n -----------\n X_valid: Validation data for features\n y_valid: Validation data for target structures\n golden_threshold: TM-score threshold to consider a seed as \"golden\"\n attempts: Number of attempts to find good seeds\n optimal_params: Optimal parameters for the reference model\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\"Searching for up to {max_seeds} diverse golden seeds with TM-score threshold of {golden_threshold}...\")\n \n # List to store all tested seeds\n all_seeds = []\n \n # List to store the \"golden\" seeds\n golden_seeds = []\n \n # List to store the predicted structures for each golden seed\n golden_predictions = []\n \n # Set of seeds already tested to avoid duplications\n tested_seeds = set()\n \n # Counter for valid attempts (excluding duplicates)\n valid_attempts = 0\n \n # Define parameter search ranges for different seed ranges\n seed_ranges = [\n (1, 1000), # Initial range\n (1001, 10000), # Medium seeds\n (10001, 100000), # Larger seeds\n (100001, 1000000) # Very large seeds\n ]\n \n # Alternating between different ranges to promote diversity\n range_index = 0\n \n # Keep track of the best seed for each RNA size range\n best_small_rna_seed = {'seed': None, 'tm_score': 0.0} # <50 residues\n best_medium_rna_seed = {'seed': None, 'tm_score': 0.0} # 50-120 residues\n best_large_rna_seed = {'seed': None, 'tm_score': 0.0} # >120 residues\n \n # Calculate sequence length statistics\n seq_lengths = []\n for coords in y_valid:\n valid_mask = ~np.all(coords == 0, axis=1)\n seq_length = np.sum(valid_mask)\n seq_lengths.append(seq_length)\n \n # Separate indices by size\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 print(f\"RNA Distribution: {len(small_rna_indices)} small, {len(medium_rna_indices)} medium, {len(large_rna_indices)} large\")\n \n # Main cycle to search for seeds\n while valid_attempts < attempts and len(golden_seeds) < max_seeds:\n # Select seed range\n min_seed, max_seed = seed_ranges[range_index]\n range_index = (range_index + 1) % len(seed_ranges)\n \n # Generate random seed from this 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}/{attempts} (seed={seed})...\")\n \n # Set the seed for reproducibility\n np.random.seed(seed)\n \n # Create model with this seed\n try:\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 if model is None:\n print(f\" Failed to create model with seed {seed}\")\n continue\n \n # Evaluate the model on different validation subsets\n # Calculate overall TM-score\n metrics = evaluate_model(model, X_valid, y_valid)\n tm_score = metrics['avg_tm_score']\n \n # Check for overfitting using TM-score on different subsets\n if len(small_rna_indices) > 0:\n small_metrics = evaluate_model_on_indices(model, X_valid, y_valid, small_rna_indices)\n small_tm_score = small_metrics['avg_tm_score']\n else:\n small_tm_score = 0.0\n \n if len(medium_rna_indices) > 0:\n medium_metrics = evaluate_model_on_indices(model, X_valid, y_valid, medium_rna_indices)\n medium_tm_score = medium_metrics['avg_tm_score']\n else:\n medium_tm_score = 0.0\n \n if len(large_rna_indices) > 0:\n large_metrics = evaluate_model_on_indices(model, X_valid, y_valid, large_rna_indices)\n large_tm_score = large_metrics['avg_tm_score']\n else:\n large_tm_score = 0.0\n \n # Calculate standard deviation between scores for different sizes\n # A high deviation may indicate overfitting in certain sizes\n size_scores = [s for s in [small_tm_score, medium_tm_score, large_tm_score] if s > 0]\n size_std = np.std(size_scores) if len(size_scores) > 1 else 0.0\n \n # Penalize the score for high variability between sizes (possible overfitting)\n adjusted_tm_score = tm_score - size_std\n \n # Register this seed\n seed_info = {\n 'seed': seed,\n 'tm_score': tm_score,\n 'adjusted_tm_score': adjusted_tm_score,\n 'small_tm_score': small_tm_score,\n 'medium_tm_score': medium_tm_score,\n 'large_tm_score': large_tm_score,\n 'size_std': size_std\n }\n all_seeds.append(seed_info)\n \n # Update the best seeds by size\n if small_tm_score > best_small_rna_seed['tm_score'] and small_tm_score > golden_threshold:\n best_small_rna_seed = {'seed': seed, 'tm_score': small_tm_score}\n \n if medium_tm_score > best_medium_rna_seed['tm_score'] and medium_tm_score > golden_threshold:\n best_medium_rna_seed = {'seed': seed, 'tm_score': medium_tm_score}\n \n if large_tm_score > best_large_rna_seed['tm_score'] and large_tm_score > golden_threshold:\n best_large_rna_seed = {'seed': seed, 'tm_score': large_tm_score}\n \n # Check if this is a \"golden\" seed\n if adjusted_tm_score >= golden_threshold:\n # Generate predictions for diversity comparison\n preds = model.predict(X_valid)\n \n # Check diversity relative to seeds already found\n is_diverse = True\n for i, existing_preds in enumerate(golden_predictions):\n similarity = calculate_prediction_similarity(preds, existing_preds)\n if similarity > (1.0 - diversity_threshold):\n is_diverse = False\n # If the new one is better than an existing one and they are similar, we replace\n if adjusted_tm_score > golden_seeds[i]['adjusted_tm_score']:\n print(f\" Replacing seed {golden_seeds[i]['seed']} (score={golden_seeds[i]['adjusted_tm_score']:.4f}) \" \n f\"with seed {seed} (score={adjusted_tm_score:.4f})\")\n golden_seeds[i] = seed_info\n golden_predictions[i] = preds\n break\n \n if is_diverse and len(golden_seeds) < max_seeds:\n print(f\" Found golden seed: {seed} (TM-score: {tm_score:.4f}, Adjusted: {adjusted_tm_score:.4f})\")\n golden_seeds.append(seed_info)\n golden_predictions.append(preds)\n \n if len(golden_seeds) >= max_seeds:\n print(f\" Reached maximum number of {max_seeds} golden seeds.\")\n break\n \n except Exception as e:\n print(f\" Error testing seed {seed}: {str(e)}\")\n continue\n \n # If we didn't find enough golden seeds, include the best by size\n if len(golden_seeds) < max_seeds:\n # Add the best seeds from each size category, if not already included\n special_seeds = [\n best_small_rna_seed,\n best_medium_rna_seed,\n best_large_rna_seed\n ]\n \n for special in special_seeds:\n if special['seed'] is not None:\n # Check if this seed is already in the golden ones\n if not any(gs['seed'] == special['seed'] for gs in golden_seeds):\n # Find the complete details of this seed in all_seeds\n for seed_detail in all_seeds:\n if seed_detail['seed'] == special['seed']:\n golden_seeds.append(seed_detail)\n break\n \n if len(golden_seeds) >= max_seeds:\n break\n \n # Sort golden seeds by adjusted TM-score (for better diversity and less overfitting)\n golden_seeds.sort(key=lambda x: x['adjusted_tm_score'], reverse=True)\n \n # Show statistics of the found seeds\n print(f\"Found {len(golden_seeds)} golden seeds in {valid_attempts} attempts\")\n for i, gs in enumerate(golden_seeds):\n print(f\" Seed {i+1}: {gs['seed']} (TM-score: {gs['tm_score']:.4f}, Adjusted: {gs['adjusted_tm_score']:.4f})\")\n print(f\" TM-scores by size - Small: {gs['small_tm_score']:.4f}, Medium: {gs['medium_tm_score']:.4f}, Large: {gs['large_tm_score']:.4f}\")\n print(f\" Standard deviation between sizes: {gs['size_std']:.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-04-05T16:08:10.766664Z","iopub.execute_input":"2025-04-05T16:08:10.766876Z","iopub.status.idle":"2025-04-05T16:08:10.799929Z","shell.execute_reply.started":"2025-04-05T16:08:10.766859Z","shell.execute_reply":"2025-04-05T16:08:10.798939Z"},"papermill":{"duration":0.051213,"end_time":"2025-03-26T03:43:08.322437","exception":false,"start_time":"2025-03-26T03:43:08.271224","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Phase 2: Quality Assessment Model","metadata":{"papermill":{"duration":0.017635,"end_time":"2025-03-26T03:43:08.358315","exception":false,"start_time":"2025-03-26T03:43:08.34068","status":"completed"},"tags":[]}},{"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-04-05T16:08:10.800961Z","iopub.execute_input":"2025-04-05T16:08:10.801262Z","iopub.status.idle":"2025-04-05T16:08:10.831973Z","shell.execute_reply.started":"2025-04-05T16:08:10.801233Z","shell.execute_reply":"2025-04-05T16:08:10.831104Z"},"papermill":{"duration":0.046737,"end_time":"2025-03-26T03:43:08.423264","exception":false,"start_time":"2025-03-26T03:43:08.376527","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"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-04-05T16:08:10.832878Z","iopub.execute_input":"2025-04-05T16:08:10.833164Z","iopub.status.idle":"2025-04-05T16:08:10.865935Z","shell.execute_reply.started":"2025-04-05T16:08:10.833143Z","shell.execute_reply":"2025-04-05T16:08:10.865187Z"},"papermill":{"duration":0.049052,"end_time":"2025-03-26T03:43:08.490394","exception":false,"start_time":"2025-03-26T03:43:08.441342","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Phase 3: Base Structure Generation","metadata":{"papermill":{"duration":0.017469,"end_time":"2025-03-26T03:43:08.525807","exception":false,"start_time":"2025-03-26T03:43:08.508338","status":"completed"},"tags":[]}},{"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-04-05T16:08:10.866883Z","iopub.execute_input":"2025-04-05T16:08:10.867157Z","iopub.status.idle":"2025-04-05T16:08:10.890865Z","shell.execute_reply.started":"2025-04-05T16:08:10.867138Z","shell.execute_reply":"2025-04-05T16:08:10.890266Z"},"papermill":{"duration":0.029881,"end_time":"2025-03-26T03:43:08.573358","exception":false,"start_time":"2025-03-26T03:43:08.543477","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Phase 4: Candidate Generation and Pruning","metadata":{"papermill":{"duration":0.01798,"end_time":"2025-03-26T03:43:08.609405","exception":false,"start_time":"2025-03-26T03:43:08.591425","status":"completed"},"tags":[]}},{"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 adaptive_noise_adjustment(structure, base_noise, ideal_length=3.8, lower_bound=0.05, divisor_mean=3.8, divisor_std=1.0):\n \"\"\"\n Adaptively adjusts the noise level based on the average and standard deviation of bond lengths.\n\n Parameters:\n -----------\n structure: Array containing the structure coordinates (shape: (seq_length, 3))\n base_noise: The previously selected base noise level\n ideal_length: Ideal bond length (default: 3.8 Å)\n lower_bound: Minimum noise level to maintain\n divisor_mean: Scaling factor for the mean deviation from the ideal bond length\n divisor_std: Scaling factor for the standard deviation of bond lengths\n\n Returns:\n --------\n The adjusted noise level (float)\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n if len(valid_coords) < 2:\n return lower_bound\n bond_lengths = [np.linalg.norm(valid_coords[i] - valid_coords[i - 1]) for i in range(1, len(valid_coords))]\n avg_length = np.mean(bond_lengths)\n std_length = np.std(bond_lengths)\n # Compute the adjustment factor based on the mean deviation and standard deviation\n adjustment_factor = 1 + abs(avg_length - ideal_length) / divisor_mean + std_length / divisor_std\n # Apply the adjustment, ensuring the noise level does not fall below lower_bound\n adjusted_noise = max(base_noise * adjustment_factor, lower_bound)\n return adjusted_noise\n\ndef calculate_bond_lengths(coords):\n \"\"\"\n Calculates the lengths of consecutive bonds based on the valid residues of the structure.\n \"\"\"\n valid_mask = ~np.all(coords == 0, axis=1)\n valid_indices = np.where(valid_mask)[0]\n if len(valid_indices) < 2:\n return []\n bond_lengths = []\n for i in range(1, len(valid_indices)):\n idx1 = valid_indices[i - 1]\n idx2 = valid_indices[i]\n bond_lengths.append(np.linalg.norm(coords[idx2] - coords[idx1]))\n return bond_lengths\n\ndef compute_scaling_factor(bond_lengths, group_size, approach=\"adaptive\", lower_bound=0.05):\n \"\"\"\n Computes a scaling factor to adjust the noise level based on the mean and standard deviation\n of bond lengths and the RNA size group.\n\n Parameters:\n - bond_lengths: list or array with the bond lengths.\n - group_size: string indicating \"small\", \"medium\" or \"large\".\n - approach: 'fixed' to use fixed values or 'adaptive' to adjust according to the group.\n - lower_bound: minimum value for the scaling factor (prevents too low scales).\n\n Returns:\n - scaling_factor: value that can be multiplied by the base noise level.\n \"\"\"\n mean_length = np.mean(bond_lengths)\n std_length = np.std(bond_lengths)\n \n # Fixed (ideal) values for comparison\n fixed_divisor_mean = 3.8\n fixed_divisor_std = 1.0\n\n if approach == \"fixed\":\n divisor_mean = fixed_divisor_mean\n divisor_std = fixed_divisor_std\n elif approach == \"adaptive\":\n if group_size == \"small\":\n divisor_mean = 3.5 # Can be adjusted for small RNAs\n divisor_std = 0.8\n elif group_size == \"medium\":\n divisor_mean = 3.8 # For medium RNAs\n divisor_std = 1.0\n elif group_size == \"large\":\n divisor_mean = 4.2 # For large RNAs\n divisor_std = 1.2\n else:\n divisor_mean = fixed_divisor_mean\n divisor_std = fixed_divisor_std\n else:\n raise ValueError(\"Approach must be 'fixed' or 'adaptive'\")\n\n scaling_factor = lower_bound + (mean_length / divisor_mean) + (std_length / divisor_std)\n return scaling_factor\n\ndef determine_optimal_noise(seq, seq_length, gc_content=None):\n \"\"\"\n Determine the optimal noise level based on sequence characteristics.\n \n Parameters:\n -----------\n seq: str\n RNA sequence\n seq_length: int\n Length of the sequence\n gc_content: float, optional\n GC content (calculated if not provided)\n \n Returns:\n --------\n float: Recommended noise level\n \"\"\"\n # Calculate GC% if not provided\n if gc_content is None:\n gc_content = sum(1 for n in seq if n in 'GC') / seq_length\n \n # Base noise adjustments based on sequence length\n if seq_length < 50:\n base_noise = 0.15\n elif seq_length < 200:\n # Linear reduction as sequence gets longer\n base_noise = 0.15 - 0.05 * (seq_length - 50) / 150\n else:\n base_noise = 0.10 - 0.05 * min(1.0, (seq_length - 200) / 500)\n \n # Adjust based on GC content (higher GC% typically means more structure)\n gc_factor = 1.0\n if gc_content > 0.6: # High GC content\n gc_factor = 0.8 # Reduce noise for highly structured RNAs\n elif gc_content < 0.4: # Low GC content\n gc_factor = 1.2 # Increase noise for less structured RNAs\n \n # Final noise level with constraints\n noise = base_noise * gc_factor\n return max(0.02, min(0.25, noise)) # Keep within reasonable bounds\n\ndef select_optimal_mode_and_params(seq_length, sequence_features, parameter_optimizer=None):\n \"\"\"\n Selects the optimal modeling mode and parameters based on sequence length\n and specific sequence features.\n \n Parameters\n ----------\n seq_length : int\n Length of the sequence\n sequence_features : dict\n Sequence characteristics (e.g., GC content, complexity, etc.)\n parameter_optimizer : ParameterOptimizer, optional\n Parameter optimizer with performance history\n \n Returns\n -------\n mode : str\n Either 'adaptive', 'fixed', or 'hybrid'\n params : dict\n Optimized parameters\n \"\"\"\n # Determine size category\n if seq_length < 50:\n size_category = 'small'\n elif seq_length < 200:\n size_category = 'medium'\n else:\n size_category = 'large'\n \n # If a parameter optimizer with historical data is available, use it\n if parameter_optimizer is not None:\n best_params = parameter_optimizer.suggest_parameters(size_category)\n if best_params:\n # Use historical parameters if available\n return best_params.get('mode', 'hybrid'), best_params\n \n # Otherwise, select strategy based on size and features\n gc_content = sequence_features.get('gc_content', 0.5)\n \n if size_category == 'small':\n # For small RNAs, use hybrid mode with parameters based on GC content\n mode = 'hybrid'\n noise_base = determine_optimal_noise(None, seq_length, gc_content)\n params = {\n 'noise_base': noise_base,\n 'divisor_mean': 3.4 if gc_content > 0.5 else 3.6,\n 'divisor_std': 0.7,\n 'correlation': 0.85,\n 'merge_strategy': 'segment_quality'\n }\n elif size_category == 'medium':\n # Medium RNAs benefit most from adaptive mode\n mode = 'adaptive'\n params = {\n 'noise_base': determine_optimal_noise(None, seq_length, gc_content),\n 'divisor_mean': 3.6,\n 'divisor_std': 0.85, # Increased for flexibility\n 'correlation': 0.82\n }\n else: # large\n # Large RNAs perform better with fixed mode plus adaptive regions\n mode = 'hybrid'\n params = {\n 'noise_base': determine_optimal_noise(None, seq_length, gc_content),\n 'divisor_mean': 3.3,\n 'divisor_std': 0.6, # Reduced for stability\n 'correlation': 0.9,\n 'merge_strategy': 'regional' # Specialized strategy for large RNAs\n }\n \n return mode, params\n\ndef sensitivity_analysis(parameter_ranges, evaluate_func, n_samples=100):\n \"\"\"\n Performs sensitivity analysis to identify the relative importance\n of parameters for the structural quality.\n\n Parameters:\n -----------\n parameter_ranges: dict\n Dictionary with lower and upper bounds for each parameter.\n Example: {'divisor_mean': (3.0, 4.5), 'noise_level': (0.01, 0.3)}\n evaluate_func: function\n Function that evaluates quality given a set of parameters.\n n_samples: int\n Number of samples to generate for the analysis.\n\n Returns:\n --------\n dict: Results of the sensitivity analysis\n \"\"\"\n # Generate samples within the specified ranges\n param_names = list(parameter_ranges.keys())\n samples = []\n\n for _ in range(n_samples):\n sample = {}\n for param, (min_val, max_val) in parameter_ranges.items():\n sample[param] = min_val + random.random() * (max_val - min_val)\n samples.append(sample)\n\n # Evaluate each sample\n results = []\n for sample in samples:\n quality_score = evaluate_func(sample)\n results.append((sample, quality_score))\n\n # Sort results by quality\n results.sort(key=lambda x: x[1], reverse=True)\n\n # Compute relative importance\n importances = {}\n correlations = {}\n quartiles = {}\n\n # Collect parameter values and quality scores\n param_values = {param: [] for param in param_names}\n quality_scores = []\n\n for sample, score in results:\n for param in param_names:\n param_values[param].append(sample[param])\n quality_scores.append(score)\n\n # Compute correlation\n for param in param_names:\n correlation = np.corrcoef(param_values[param], quality_scores)[0, 1]\n correlations[param] = correlation\n\n # Split into quartiles based on quality\n quartile_size = n_samples // 4\n quartiles = {\n 'Q1': results[:quartile_size], # Top 25%\n 'Q2': results[quartile_size:2*quartile_size],\n 'Q3': results[2*quartile_size:3*quartile_size],\n 'Q4': results[3*quartile_size:] # Bottom 25%\n }\n\n # Analyze parameter values by quartile\n quartile_stats = {}\n for q_name, q_results in quartiles.items():\n quartile_stats[q_name] = {}\n for param in param_names:\n param_vals = [r[0][param] for r in q_results]\n quartile_stats[q_name][param] = {\n 'mean': sum(param_vals) / len(param_vals),\n 'std': (sum((x - sum(param_vals)/len(param_vals))**2 \n for x in param_vals) / len(param_vals))**0.5,\n 'min': min(param_vals),\n 'max': max(param_vals)\n }\n\n # Estimate importance based on variation between quartiles\n for param in param_names:\n q1_mean = quartile_stats['Q1'][param]['mean']\n q4_mean = quartile_stats['Q4'][param]['mean']\n q1_std = quartile_stats['Q1'][param]['std']\n q4_std = quartile_stats['Q4'][param]['std']\n\n # The greater the difference between means of extreme quartiles and\n # the smaller the within-quartile variance, the more important the parameter\n mean_diff = abs(q1_mean - q4_mean)\n variance_ratio = (q1_std + q4_std) / 2\n if variance_ratio == 0:\n variance_ratio = 0.001 # Avoid division by zero\n\n importances[param] = mean_diff / variance_ratio\n\n # Normalize importance scores\n total_importance = sum(importances.values())\n if total_importance > 0:\n for param in importances:\n importances[param] /= total_importance\n\n return {\n 'importances': importances,\n 'correlations': correlations,\n 'quartile_stats': quartile_stats,\n 'top_samples': results[:10] # Top 10 best combinations\n }\n\ndef generate_diverse_structures_from_bases(base_structures, seq_length, quality_model, num_per_base=5, mode='adaptive'):\n \"\"\"\n Generates diverse candidate structures from a list of base structures,\n incorporating RNA-specific variations and quality filtering.\n \n Parameters:\n -----------\n base_structures: List of base structures (each as a coordinate array, shape (seq_length, 3))\n seq_length: Length of the RNA sequence (number of valid residues)\n quality_model: Model used for quality assessment (retained for compatibility)\n num_per_base: Number of structural variations to generate for each base\n mode: 'adaptive' to dynamically adjust noise, or 'fixed' to use base noise values as-is\n \n Returns:\n --------\n A list of candidate structures (coordinate arrays) after variation, refinement, and normalization.\n \"\"\"\n candidates = []\n\n # Include the original base structures\n for base in base_structures:\n candidates.append(base)\n\n # Generate structural 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 # Define variation parameters based on RNA size\n if seq_length < 50:\n noise_levels = [0.05, 0.1, 0.15, 0.2, 0.25]\n preserve_distances = [True, True, True, False, False]\n use_globals = [False, False, True, False, True]\n elif seq_length < 120:\n noise_levels = [0.03, 0.06, 0.1, 0.15, 0.2]\n preserve_distances = [True, True, True, True, False]\n use_globals = [False, True, False, True, False]\n else:\n noise_levels = [0.02, 0.04, 0.06, 0.08, 0.1]\n preserve_distances = [True, True, True, True, True]\n use_globals = [False, False, True, False, True]\n\n # Generate the desired number of variations\n for i in range(num_per_base):\n noise_idx = i % len(noise_levels)\n base_noise = noise_levels[noise_idx]\n preserve_distance = preserve_distances[noise_idx]\n use_global = use_globals[noise_idx]\n\n # Set a unique seed for reproducibility\n np.random.seed(base_idx * 100 + i)\n\n # Adjust noise dynamically if in 'adaptive' mode\n if mode == 'adaptive':\n adjusted_noise = adaptive_noise_adjustment(\n base,\n base_noise,\n ideal_length=3.8,\n lower_bound=0.05,\n divisor_mean=3.8,\n divisor_std=1.0\n )\n else:\n adjusted_noise = base_noise # Fixed mode: use base noise as-is\n\n # Apply structural variation\n variation = sample_structural_variation(\n base,\n noise_level=adjusted_noise,\n preserve_distance=preserve_distance,\n use_global_movement=use_global,\n correlation=0.8 + np.random.uniform(-0.1, 0.1)\n )\n\n # Refine the backbone using dihedral angle correction\n variation = refine_rna_backbone_with_dihedrals(variation, ideal_dihedral=180.0)\n # Normalize the structure (center, remove padding, etc.)\n normalized = normalize_structure(variation)\n\n # Check for validity and attempt to repair if needed\n if check_structure_validity(normalized):\n candidates.append(normalized)\n else:\n print(\" Invalid structure detected. Attempting repair.\")\n repaired = repair_invalid_structure(normalized)\n if check_structure_validity(repaired):\n candidates.append(repaired)\n print(\" Structure successfully repaired.\")\n\n print(f\"Generated {len(candidates)} candidate structures in total.\")\n\n # Pre-filter based on a simple bond-length quality metric\n if len(candidates) > 30:\n print(\"Pre-filtering candidates based on basic quality metrics...\")\n quality_scores = []\n for candidate in candidates:\n valid_mask = ~np.all(candidate == 0, axis=1)\n valid_coords = candidate[valid_mask]\n if len(valid_coords) < 3:\n quality_scores.append(0.0)\n continue\n # Calculate bond lengths\n bond_lengths = [np.linalg.norm(valid_coords[j] - valid_coords[j - 1]) for j in range(1, len(valid_coords))]\n avg_bond_length = np.mean(bond_lengths)\n # Simple score: the closer to 3.8 Å, the better\n bond_score = 1.0 - min(1.0, abs(avg_bond_length - 3.8) / 3.8)\n quality_scores.append(bond_score)\n quality_scores = np.array(quality_scores)\n top_indices = np.argsort(quality_scores)[-30:]\n candidates = [candidates[idx] for idx in top_indices]\n print(\"Pre-filtering completed: top 30 candidates retained.\")\n\n return candidates\n\ndef group_label(seq_length):\n if seq_length < 50:\n return 'small'\n elif seq_length < 120:\n return 'medium'\n else:\n return 'large'\n\ndef compute_quality_metrics(candidates):\n scores = []\n for candidate in candidates:\n valid_mask = ~np.all(candidate == 0, axis=1)\n valid_coords = candidate[valid_mask]\n if len(valid_coords) < 3:\n scores.append(0.0)\n continue\n bond_lengths = [np.linalg.norm(valid_coords[j] - valid_coords[j - 1]) for j in range(1, len(valid_coords))]\n avg_bond_length = np.mean(bond_lengths)\n score = 1.0 - min(1.0, abs(avg_bond_length - 3.8) / 3.8)\n scores.append(score)\n return np.mean(scores), np.std(scores)\n\ndef evaluate_candidates_by_group(candidates):\n \"\"\"\n Evaluate candidates by grouping them into size categories and calculating metrics.\n \n Parameters:\n -----------\n candidates: List of candidate structures\n \n Returns:\n --------\n Dictionary with metrics for each group\n \"\"\"\n # Group candidates by size\n small = [c for c in candidates if len(c) < 50]\n medium = [c for c in candidates if 50 <= len(c) < 100]\n large = [c for c in candidates if len(c) >= 100]\n \n # Calculate metrics for each group\n results = {}\n \n if small:\n results['small'] = calculate_group_metrics(small)\n \n if medium:\n results['medium'] = calculate_group_metrics(medium)\n \n if large:\n results['large'] = calculate_group_metrics(large)\n \n return results\n\n\ndef calculate_group_metrics(structures):\n \"\"\"\n Calculate comprehensive quality metrics for a group of structures,\n with improved handling of structures of all sizes.\n \n Parameters:\n -----------\n structures: List of structures in the group\n \n Returns:\n --------\n Dictionary of metrics\n \"\"\"\n if not structures:\n return {\n 'count': 0,\n 'avg_energy': 0,\n 'avg_compactness': 0,\n 'structural_diversity': 0,\n 'quality_score': 0,\n 'size_category': 'unknown'\n }\n \n # Determine size category\n first_valid_mask = ~np.all(structures[0] == 0, axis=1)\n seq_length = np.sum(first_valid_mask)\n \n if seq_length < 50:\n size_category = 'small'\n # Small structures: prioritize compactness more\n weight_energy = 0.35\n weight_compactness = 0.40\n weight_diversity = 0.25\n elif seq_length < 200:\n size_category = 'medium'\n # Medium structures: balanced weights\n weight_energy = 0.40\n weight_compactness = 0.30\n weight_diversity = 0.30\n else:\n size_category = 'large'\n # Large structures: prioritize energy and diversity\n weight_energy = 0.45\n weight_compactness = 0.25\n weight_diversity = 0.30\n \n # Calculate energy (lower is better, so we use 1-energy for the final score)\n energies = [calculate_energy(s) for s in structures]\n avg_energy = sum(energies) / len(structures)\n energy_score = 1 - avg_energy # Invert so higher is better\n \n # Calculate compactness (higher is better)\n compactness_values = [calculate_compactness(s) for s in structures]\n avg_compactness = sum(compactness_values) / len(structures)\n \n # Calculate structural diversity (higher is better)\n structural_diversity = calculate_structural_diversity(structures)\n \n # Calculate distributions and variability\n energy_std = np.std(energies) if len(energies) > 1 else 0\n compactness_std = np.std(compactness_values) if len(compactness_values) > 1 else 0\n \n # Combine metrics into overall quality score with size-adaptive weights\n quality_score = (\n weight_energy * energy_score + \n weight_compactness * avg_compactness + \n weight_diversity * structural_diversity\n )\n \n # Add bonus for consistency in energy (if std is low)\n if energy_std < 0.1:\n quality_score += 0.02\n \n # Add bonus for consistency in compactness (if std is low)\n if compactness_std < 0.1:\n quality_score += 0.01\n \n # Detailed metrics for analysis\n return {\n 'count': len(structures),\n 'avg_energy': avg_energy,\n 'energy_std': energy_std,\n 'avg_compactness': avg_compactness,\n 'compactness_std': compactness_std,\n 'structural_diversity': structural_diversity,\n 'size_category': size_category,\n 'quality_score': quality_score\n }\n\ndef calculate_energy(structure):\n \"\"\"\n Calculate the \"energy\" of the structure based on bond lengths.\n \n Parameters:\n structure: numpy array (n, 3) with the structure coordinates.\n \n Returns:\n Normalized energy score (float) between 0-1, where lower is better.\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n if len(valid_coords) < 2:\n return 1.0 # Worst energy score for invalid structures\n \n # Calculate bond lengths\n bond_lengths = [np.linalg.norm(valid_coords[i] - valid_coords[i-1]) for i in range(1, len(valid_coords))]\n \n # Calculate energy based on deviation from ideal bond length (3.8 Å)\n deviations = [abs(length - 3.8) for length in bond_lengths]\n avg_deviation = np.mean(deviations)\n \n # Normalize to 0-1 range (0 = best, 1 = worst)\n energy_score = min(1.0, avg_deviation / 3.8)\n return energy_score\n\ndef calculate_compactness(structure):\n \"\"\"\n Calculate the compactness of the structure using radius of gyration with scaling\n based on structure size to handle structures of all sizes properly.\n \n Parameters:\n structure: numpy array (n, 3) with the structure coordinates.\n \n Returns:\n Normalized compactness score (float) between 0-1, where higher is better.\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n \n if len(valid_coords) < 3:\n return 0.0\n \n # Calculate the radius of gyration\n centroid = np.mean(valid_coords, axis=0)\n sq_dists = np.sum(np.square(valid_coords - centroid), axis=1)\n rg = np.sqrt(np.mean(sq_dists))\n \n # Calculate sequence length (number of valid coordinates)\n seq_length = len(valid_coords)\n \n # Scale the expected radius of gyration based on sequence length\n # For RNA, compactness follows approximately a power law\n # These coefficients are empirically derived for RNA structures\n expected_rg = 3.5 * (seq_length ** 0.33)\n \n # For very large structures, apply an additional correction\n if seq_length > 200:\n expected_rg *= 1.1\n \n # Calculate score: how close is the actual Rg to the expected Rg?\n # Lower variance for larger structures\n scaling_factor = min(1.0, 50.0 / seq_length) # Reduces tolerance for larger structures\n tolerance = 0.3 * scaling_factor * expected_rg\n \n # Exponential scoring function - peaks at expected_rg and falls off with distance\n # This gives better discrimination than a simple linear function\n score = np.exp(-0.5 * ((rg - expected_rg) / tolerance) ** 2)\n \n # Add secondary assessment for local compactness\n # For large structures, also assess compactness of local regions\n local_compactness = 1.0\n if seq_length > 100:\n # Sample local regions and calculate their compactness\n window_size = min(50, seq_length // 4)\n n_samples = min(5, seq_length // window_size)\n \n local_scores = []\n for i in range(n_samples):\n start_idx = random.randint(0, seq_length - window_size)\n window = valid_coords[start_idx:start_idx + window_size]\n window_centroid = np.mean(window, axis=0)\n window_rg = np.sqrt(np.mean(np.sum(np.square(window - window_centroid), axis=1)))\n window_expected_rg = 3.5 * (window_size ** 0.33)\n window_score = np.exp(-0.5 * ((window_rg - window_expected_rg) / (0.3 * window_expected_rg)) ** 2)\n local_scores.append(window_score)\n \n if local_scores:\n local_compactness = np.mean(local_scores)\n \n # Weight global and local compactness\n if seq_length > 100:\n final_score = 0.6 * score + 0.4 * local_compactness\n else:\n final_score = score\n \n return final_score\n\ndef rmsd(struct1, struct2):\n \"\"\"\n Calculate the RMSD (Root-Mean-Square Deviation) between two structures.\n Optimized for handling large structures efficiently.\n \n Parameters:\n struct1, struct2: numpy arrays (n, 3).\n \n Returns:\n RMSD (float).\n \"\"\"\n valid_mask1 = ~np.all(struct1 == 0, axis=1)\n valid_mask2 = ~np.all(struct2 == 0, axis=1)\n \n # Use only points that are valid in both structures\n common_mask = valid_mask1 & valid_mask2\n \n if np.sum(common_mask) < 3:\n return float('inf') # Not enough valid points for comparison\n \n points1 = struct1[common_mask]\n points2 = struct2[common_mask]\n \n # For very large structures, subsample points to improve efficiency\n n_points = len(points1)\n if n_points > 200:\n # Select regularly spaced points for efficiency\n step = n_points // 200\n indices = np.arange(0, n_points, step)\n points1 = points1[indices]\n points2 = points2[indices]\n \n # Center the structures to account for translational differences\n centroid1 = np.mean(points1, axis=0)\n centroid2 = np.mean(points2, axis=0)\n \n centered1 = points1 - centroid1\n centered2 = points2 - centroid2\n \n # Optionally, implement Kabsch algorithm for optimal alignment\n # For simplicity, we'll use a direct RMSD calculation\n diff = centered1 - centered2\n squared_dist = np.sum(diff * diff, axis=1)\n rmsd_val = np.sqrt(np.mean(squared_dist))\n \n return rmsd_val\n\ndef calculate_structural_diversity(structures):\n \"\"\"\n Calculate the average structural diversity between candidates with optimizations\n for large structures and more meaningful scaling.\n \n Parameters:\n structures: list of numpy arrays, each representing a structure (n, 3).\n \n Returns:\n Average diversity (float), calculated as the mean RMSD between all pairs.\n Normalized to 0-1 range where higher values indicate more diversity.\n \"\"\"\n n = len(structures)\n if n < 2:\n return 0.0\n \n # For large structures, limit the number of comparisons to avoid excessive computation\n max_pairs = 100\n \n # Determine valid residues for each structure\n valid_masks = [~np.all(struct == 0, axis=1) for struct in structures]\n \n # Calculate sequence length from the first structure\n seq_length = np.sum(valid_masks[0])\n \n # Adaptive scaling based on sequence length\n # Scale the expected RMSD based on sequence length\n if seq_length < 50:\n max_expected_rmsd = 10.0\n elif seq_length < 200:\n max_expected_rmsd = 15.0\n else:\n max_expected_rmsd = 20.0 + (seq_length - 200) * 0.05 # Additional scaling for very large structures\n \n # Sampling strategy based on structure size\n if n * (n-1) // 2 <= max_pairs:\n # For small enough sets, compare all pairs\n pair_indices = [(i, j) for i in range(n) for j in range(i+1, n)]\n else:\n # For larger sets, sample random pairs\n pair_indices = set()\n while len(pair_indices) < max_pairs:\n i, j = random.sample(range(n), 2)\n if i > j:\n i, j = j, i\n if (i, j) not in pair_indices:\n pair_indices.add((i, j))\n pair_indices = list(pair_indices)\n \n # For large structures, use landmark-based RMSD to improve efficiency\n use_landmarks = seq_length > 200\n n_landmarks = min(50, seq_length // 4) if use_landmarks else seq_length\n \n rmsd_values = []\n \n for i, j in pair_indices:\n struct1, struct2 = structures[i], structures[j]\n mask1, mask2 = valid_masks[i], valid_masks[j]\n \n # Use only points that are valid in both structures\n common_mask = mask1 & mask2\n \n if np.sum(common_mask) < 3:\n continue # Not enough valid points for comparison\n \n if use_landmarks:\n # Select landmark points (regularly spaced)\n landmark_indices = np.where(common_mask)[0]\n if len(landmark_indices) > n_landmarks:\n # Take regularly spaced points\n step = len(landmark_indices) // n_landmarks\n landmark_indices = landmark_indices[::step][:n_landmarks]\n \n landmark_mask = np.zeros_like(common_mask)\n landmark_mask[landmark_indices] = True\n \n points1 = struct1[landmark_mask]\n points2 = struct2[landmark_mask]\n else:\n points1 = struct1[common_mask]\n points2 = struct2[common_mask]\n \n # Calculate RMSD\n diff = points1 - points2\n rmsd_val = np.sqrt(np.mean(np.sum(diff**2, axis=1)))\n rmsd_values.append(rmsd_val)\n \n if not rmsd_values:\n return 0.0\n \n # Calculate average RMSD\n avg_rmsd = np.mean(rmsd_values)\n \n # Scale to 0-1 range using sigmoid function for better distribution\n # This gives more meaningful differentiation between moderate and high diversity\n normalized_diversity = 2.0 / (1.0 + np.exp(-avg_rmsd / (max_expected_rmsd/4.0))) - 1.0\n \n return normalized_diversity\n\ndef identify_flexible_regions(structure, ideal_angle=np.deg2rad(120), angle_threshold=np.deg2rad(20)):\n \"\"\"\n Identify flexible regions (such as loops and junctions) in a structure based on bond angles.\n \n Parameters:\n structure: np.ndarray of shape (n, 3) representing coordinates.\n ideal_angle: ideal bond angle in radians.\n angle_threshold: allowed deviation from the ideal angle in radians.\n \n Returns:\n A list of indices (or segments) where the angle deviates more than the threshold.\n \"\"\"\n n = structure.shape[0]\n flexible_indices = []\n if n < 3:\n return flexible_indices\n\n # Compute bond vectors\n vectors = structure[1:] - structure[:-1]\n \n # Compute angles between consecutive bond vectors\n for i in range(1, len(vectors)):\n v1 = vectors[i-1]\n v2 = vectors[i]\n # Avoid division by zero\n norm1 = np.linalg.norm(v1)\n norm2 = np.linalg.norm(v2)\n if norm1 == 0 or norm2 == 0:\n continue\n # Compute cosine of angle and then the angle\n cos_angle = np.clip(np.dot(v1, v2) / (norm1 * norm2), -1.0, 1.0)\n angle = np.arccos(cos_angle)\n # If deviation from the ideal angle is larger than threshold, mark the middle residue as flexible\n if np.abs(angle - ideal_angle) > angle_threshold:\n flexible_indices.append(i)\n \n return flexible_indices\n\ndef merge_regions(base_struct, donor_struct, regions):\n \"\"\"\n Merge two structures by replacing specified regions of the base structure with the donor structure.\n \n Parameters:\n base_struct: np.ndarray, the base structure coordinates.\n donor_struct: np.ndarray, the donor structure coordinates.\n regions: list of tuples (start, end) specifying the indices to replace.\n \n Returns:\n A new structure with regions merged.\n \"\"\"\n merged = base_struct.copy()\n for (start, end) in regions:\n merged[start:end+1] = donor_struct[start:end+1]\n return merged\n\ndef weighted_merge(struct_a, struct_b, weight_a, weight_b):\n \"\"\"\n Merge two structures using a weighted average of coordinates.\n \n Parameters:\n struct_a: np.ndarray, coordinates of structure A.\n struct_b: np.ndarray, coordinates of structure B.\n weight_a: float, weight for structure A.\n weight_b: float, weight for structure B.\n \n Returns:\n A new structure that is the weighted average of A and B.\n \"\"\"\n total_weight = weight_a + weight_b\n return (weight_a * struct_a + weight_b * struct_b) / total_weight\n\ndef assess_structure_quality(structure, ideal_bond=3.8):\n \"\"\"\n Assess the quality of a structure by comparing its bond lengths to an ideal value.\n \n Parameters:\n structure: np.ndarray of shape (n, 3).\n ideal_bond: float, ideal bond length.\n \n Returns:\n A quality score (higher is better). Here, we use an inverse error metric.\n \"\"\"\n # Compute bond lengths\n bonds = np.linalg.norm(np.diff(structure, axis=0), axis=1)\n error = np.mean(np.abs(bonds - ideal_bond))\n # A simple quality metric: lower error yields higher score\n quality = 1 / (1 + error)\n return quality\n\ndef get_best_structures(structures, top_k=5):\n \"\"\"\n Return the top_k best structures based on their quality scores.\n \n Parameters:\n structures: list of np.ndarray structures.\n top_k: int, number of best structures to return.\n \n Returns:\n A list of the top_k structures sorted by quality (highest first).\n \"\"\"\n scored_structs = []\n for struct in structures:\n score = assess_structure_quality(struct)\n scored_structs.append((score, struct))\n # Sort descending by score\n scored_structs.sort(key=lambda x: x[0], reverse=True)\n best_structures = [s for score, s in scored_structs[:top_k]]\n return best_structures\n\ndef remove_duplicate_structures(structures, threshold=0.5):\n \"\"\"\n Remove structures that are too similar based on RMSD.\n \n Parameters:\n structures: list of np.ndarray structures.\n threshold: float, RMSD threshold below which structures are considered duplicates.\n \n Returns:\n A list of unique structures.\n \"\"\"\n unique = []\n for s in structures:\n duplicate = False\n for u in unique:\n # Compute RMSD between s and u\n rmsd = np.sqrt(np.mean((s - u) ** 2))\n if rmsd < threshold:\n duplicate = True\n break\n if not duplicate:\n unique.append(s)\n return unique\n\ndef calculate_comprehensive_quality(structure, seq_length, ideal_bond=3.8):\n \"\"\"\n Calculate comprehensive quality metrics for a given structure.\n \n Parameters:\n structure: np.ndarray of shape (n, 3).\n seq_length: int, length of the RNA sequence.\n ideal_bond: float, ideal bond length.\n \n Returns:\n A dictionary with keys:\n - energy: average absolute deviation of bond lengths from ideal.\n - avg_compactness: average distance from the centroid.\n - structural_diversity: standard deviation of bond lengths.\n - quality_score: combined quality score.\n \"\"\"\n # Compute bond lengths\n bonds = np.linalg.norm(np.diff(structure, axis=0), axis=1)\n energy = np.mean(np.abs(bonds - ideal_bond))\n \n # Compute compactness: mean distance from centroid\n centroid = np.mean(structure, axis=0)\n distances = np.linalg.norm(structure - centroid, axis=1)\n avg_compactness = np.mean(distances)\n \n # Structural diversity: standard deviation of bond lengths\n diversity = np.std(bonds)\n \n # Combined quality score (a simple weighted inverse function)\n quality_score = 1 / (1 + energy + 0.1 * avg_compactness)\n \n return {\n \"energy\": energy,\n \"avg_compactness\": avg_compactness,\n \"structural_diversity\": diversity,\n \"quality_score\": quality_score\n }\n\ndef evaluate_with_comprehensive_metrics(candidates, seq_length):\n \"\"\"\n Evaluate candidate structures using comprehensive quality metrics.\n \n Parameters:\n candidates: list of np.ndarray candidate structures.\n seq_length: int, length of the RNA sequence.\n \n Returns:\n A dictionary summarizing the average metrics over all candidates.\n \"\"\"\n metrics = [calculate_comprehensive_quality(candidate, seq_length) for candidate in candidates]\n \n avg_energy = np.mean([m[\"energy\"] for m in metrics])\n avg_compactness = np.mean([m[\"avg_compactness\"] for m in metrics])\n avg_diversity = np.mean([m[\"structural_diversity\"] for m in metrics])\n avg_quality = np.mean([m[\"quality_score\"] for m in metrics])\n \n return {\n \"count\": len(candidates),\n \"avg_energy\": avg_energy,\n \"avg_compactness\": avg_compactness,\n \"structural_diversity\": avg_diversity,\n \"quality_score\": avg_quality\n }\n\ndef predict_secondary_structure(structure):\n \"\"\"\n Predicts RNA secondary structure from 3D coordinates.\n\n Parameters:\n -----------\n structure: numpy.ndarray\n 3D coordinates of the structure (N x 3)\n\n Returns:\n --------\n list: Predicted secondary structure (e.g., 'H' for helix, 'S' for strand, 'L' for loop)\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 4:\n return ['L'] * len(valid_coords)\n\n secondary_structure = []\n dist_matrix = np.zeros((len(valid_coords), len(valid_coords)))\n\n for i in range(len(valid_coords)):\n for j in range(len(valid_coords)):\n dist_matrix[i, j] = np.linalg.norm(valid_coords[i] - valid_coords[j])\n\n for i in range(len(valid_coords)):\n if i+3 < len(valid_coords) and dist_matrix[i, i+3] < 10.0:\n secondary_structure.append('H')\n elif i+2 < len(valid_coords) and dist_matrix[i, i+2] > 8.0:\n secondary_structure.append('S')\n else:\n secondary_structure.append('L')\n\n smoothed = secondary_structure.copy()\n for i in range(1, len(secondary_structure)-1):\n if secondary_structure[i-1] == secondary_structure[i+1] and secondary_structure[i] != secondary_structure[i-1]:\n smoothed[i] = secondary_structure[i-1]\n\n return smoothed\n\ndef evaluate_secondary_structure(secondary_structure):\n \"\"\"\n Evaluates the quality of a secondary structure based on regularity and patterns.\n\n Parameters:\n -----------\n secondary_structure: list\n List of secondary structure elements ('H', 'S', 'L')\n\n Returns:\n --------\n float: Quality score between 0 and 1\n \"\"\"\n if not secondary_structure:\n return 0.0\n\n segment_lengths = []\n current_type = secondary_structure[0]\n current_length = 1\n\n for ss_type in secondary_structure[1:]:\n if ss_type == current_type:\n current_length += 1\n else:\n segment_lengths.append(current_length)\n current_type = ss_type\n current_length = 1\n segment_lengths.append(current_length)\n\n avg_segment_length = np.mean(segment_lengths)\n max_segment_length = max(segment_lengths)\n\n type_counts = {\n 'H': secondary_structure.count('H'),\n 'S': secondary_structure.count('S'),\n 'L': secondary_structure.count('L')\n }\n\n total = len(secondary_structure)\n type_fractions = {k: v/total for k, v in type_counts.items()}\n\n helix_bonus = type_fractions.get('H', 0) * 0.5\n loop_penalty = max(0, type_fractions.get('L', 0) - 0.4) * 0.3\n length_bonus = min(1.0, avg_segment_length / 5.0) * 0.2\n\n score = 0.5 + helix_bonus + length_bonus - loop_penalty\n return min(1.0, max(0.0, score))\n\ndef calculate_packing_density(structure):\n \"\"\"\n Calculates the packing density of a 3D structure.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure (shape: N x 3)\n\n Returns\n -------\n float\n Normalized packing density score (range: 0 to 1)\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 4:\n return 0.5 # Default value for very small structures\n\n try:\n # Compute the structure volume using convex hull\n from scipy.spatial import ConvexHull\n hull = ConvexHull(valid_coords)\n volume = hull.volume\n\n # Normalize by the number of residues\n normalized_volume = volume / len(valid_coords)\n\n # Normalize to range [0, 1] (lower volume implies better packing)\n # Typically, good packing volume per residue is around 30–50 ų\n if normalized_volume < 30:\n return 1.0\n elif normalized_volume > 100:\n return 0.0\n else:\n return max(0.0, min(1.0, 1.0 - (normalized_volume - 30) / 70))\n except:\n # Fallback if ConvexHull fails — use radius of gyration as a proxy\n centroid = np.mean(valid_coords, axis=0)\n radii = np.linalg.norm(valid_coords - centroid, axis=1)\n rg = np.mean(radii)\n\n # Normalize radius of gyration to a packing score\n # Smaller radius implies tighter packing\n expected_rg = 3.5 * (len(valid_coords) ** 0.33)\n return max(0.0, min(1.0, 1.0 - abs(rg - expected_rg) / expected_rg))\n\ndef calculate_contact_order(structure):\n \"\"\"\n Calculates the contact order of a structure, which measures\n the sequence separation of spatially close residues.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure (shape: N x 3)\n\n Returns\n -------\n float\n Normalized contact order score (range: 0 to 1)\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 4:\n return 0.5 # Default value for very small structures\n\n contact_threshold = 8.0 # Ångströms\n n_residues = len(valid_coords)\n contacts = []\n\n for i in range(n_residues):\n for j in range(i + 3, n_residues): # Ignore sequential contacts\n if np.linalg.norm(valid_coords[i] - valid_coords[j]) < contact_threshold:\n contacts.append(abs(i - j))\n\n if not contacts:\n return 0.3 # Few contacts usually indicates poor folding\n\n # Contact order (CO) = average sequence separation of contacting pairs\n co = np.mean(contacts) / n_residues\n\n # Normalize CO to a score (typical CO values are 0.1–0.3)\n # Higher CO generally indicates more complex topology\n if co < 0.1:\n return 0.3 # Too low — overly simple structure\n elif co > 0.4:\n return 0.7 # High CO — potentially interesting topology\n else:\n return 0.3 + co # Intermediate values are scaled\n\ndef identify_domains(structure):\n \"\"\"\n Identifies structural domains based on coordinate clustering.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure (shape: N x 3)\n\n Returns\n -------\n list\n Domain assignment for each residue\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 30:\n return [0] * len(valid_coords) # Small structures considered a single domain\n\n try:\n from sklearn.cluster import KMeans, DBSCAN\n\n # Estimate number of domains based on structure size\n n_domains = max(1, len(valid_coords) // 100)\n\n # Try DBSCAN clustering\n clustering = DBSCAN(eps=15.0, min_samples=5).fit(valid_coords)\n domains = clustering.labels_\n\n # Fallback to KMeans if DBSCAN fails\n if len(set(domains)) <= 1 or -1 in domains:\n kmeans = KMeans(n_clusters=n_domains, random_state=42).fit(valid_coords)\n domains = kmeans.labels_\n except:\n # Simple fallback using coordinate slicing\n domains = []\n for i in range(len(valid_coords)):\n domains.append(min(i * n_domains // len(valid_coords), n_domains - 1))\n\n return domains\n\ndef evaluate_domain_separation(domains):\n \"\"\"\n Evaluates the separation of structural domains. A well-formed structure \n should have clearly defined and well-separated domains.\n\n Parameters\n ----------\n domains : list or numpy.ndarray\n List of domain assignments for each residue\n\n Returns\n -------\n float\n Domain separation quality score (range: 0 to 1)\n \"\"\"\n # Check if input is None or empty\n if domains is None or len(domains) == 0:\n return 0.5\n\n # Convert to regular Python list if it's a NumPy array\n if isinstance(domains, np.ndarray):\n domains = domains.tolist()\n\n # If only a single domain is present\n unique_domains = set(domains)\n if len(unique_domains) <= 1:\n return 0.5\n\n # Identify domain transition boundaries\n boundaries = []\n prev_domain = domains[0]\n for i, domain in enumerate(domains[1:], 1):\n if domain != prev_domain:\n boundaries.append(i)\n prev_domain = domain\n\n # Assess the domain size distribution\n domain_sizes = {}\n for d in domains:\n domain_sizes[d] = domain_sizes.get(d, 0) + 1\n\n sizes = list(domain_sizes.values())\n mean_size = np.mean(sizes)\n size_variance = 1.0 if mean_size == 0 else np.std(sizes) / mean_size\n\n # Penalty for high domain size variance\n size_penalty = min(1.0, size_variance)\n\n # Penalty for excessive domain transitions (fragmentation)\n fragment_penalty = min(1.0, len(boundaries) / (len(domains) / 20))\n\n # Bonus for ideal number of domains (typically 2–4)\n n_domains = len(domain_sizes)\n if n_domains <= 1:\n domain_bonus = 0.0\n elif n_domains <= 4:\n domain_bonus = 0.2 * n_domains\n else:\n domain_bonus = 0.8 - (n_domains - 4) * 0.1\n\n score = 0.5 + domain_bonus - 0.25 * size_penalty - 0.25 * fragment_penalty\n return max(0.0, min(1.0, score))\n\ndef evaluate_tertiary_packing(structure):\n \"\"\"\n Evaluates the quality of tertiary structure packing.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure (shape: N x 3)\n\n Returns\n -------\n float\n Tertiary packing quality score (range: 0 to 1)\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 10:\n return 0.5 # Default for very small structures\n\n long_range_threshold = 8.0 # Ångströms\n seq_distance_threshold = 7 # Residues\n\n n_residues = len(valid_coords)\n n_possible_contacts = 0\n n_actual_contacts = 0\n\n for i in range(n_residues):\n for j in range(i + seq_distance_threshold, n_residues):\n n_possible_contacts += 1\n if np.linalg.norm(valid_coords[i] - valid_coords[j]) < long_range_threshold:\n n_actual_contacts += 1\n\n if n_possible_contacts == 0:\n return 0.5\n\n contact_density = n_actual_contacts / n_possible_contacts\n\n # Normalize to score\n if contact_density < 0.05:\n return 0.3 # Too few contacts\n elif contact_density > 0.3:\n return 0.7 # Very good packing\n else:\n return 0.3 + contact_density * 2.0\n\ndef evaluate_long_range_contacts(structure):\n \"\"\"\n Specifically evaluates long-range contacts in the structure.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure (shape: N x 3)\n\n Returns\n -------\n float\n Score for quality of long-range contacts (range: 0 to 1)\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 20:\n return 0.5 # Default for small structures\n\n contact_threshold = 10.0 # Ångströms\n short_range_threshold = 10 # Residues\n medium_range_threshold = 30 # Residues\n\n n_residues = len(valid_coords)\n short_range_contacts = 0\n medium_range_contacts = 0\n long_range_contacts = 0\n\n for i in range(n_residues):\n for j in range(i + short_range_threshold, n_residues):\n dist = np.linalg.norm(valid_coords[i] - valid_coords[j])\n if dist < contact_threshold:\n seq_dist = j - i\n if seq_dist < medium_range_threshold:\n short_range_contacts += 1\n elif seq_dist < n_residues // 2:\n medium_range_contacts += 1\n else:\n long_range_contacts += 1\n\n short_range_density = short_range_contacts / max(1, n_residues - short_range_threshold)\n medium_range_density = medium_range_contacts / max(1, n_residues - medium_range_threshold)\n long_range_density = long_range_contacts / max(1, n_residues // 2)\n\n # Weighted combination (long-range contacts are more valuable)\n combined_score = (\n 0.2 * short_range_density +\n 0.3 * medium_range_density +\n 0.5 * long_range_density\n )\n\n # Normalize to final score\n return min(1.0, combined_score * 5.0)\n\ndef calculate_bond_angles(structure):\n \"\"\"\n Calculates bond angles in the structure.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure (shape: N x 3)\n\n Returns\n -------\n list\n Bond angles in degrees\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 3:\n return []\n\n angles = []\n for i in range(1, len(valid_coords) - 1):\n # Bond vectors\n v1 = valid_coords[i - 1] - valid_coords[i]\n v2 = valid_coords[i + 1] - valid_coords[i]\n\n # Normalize vectors\n v1 = v1 / np.linalg.norm(v1)\n v2 = v2 / np.linalg.norm(v2)\n\n # Compute angle\n cos_angle = np.clip(np.dot(v1, v2), -1.0, 1.0)\n angle = np.arccos(cos_angle) * 180 / np.pi\n angles.append(angle)\n\n return angles\n\ndef calculate_dihedral_angles(structure):\n \"\"\"\n Calculates dihedral (torsion) angles in the structure.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure (shape: N x 3)\n\n Returns\n -------\n list\n Dihedral angles in degrees\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 4:\n return []\n\n dihedrals = []\n for i in range(len(valid_coords) - 3):\n p1, p2, p3, p4 = valid_coords[i:i + 4]\n\n # Bond vectors\n b1 = p2 - p1\n b2 = p3 - p2\n b3 = p4 - p3\n\n # Normal vectors to planes\n n1 = np.cross(b1, b2)\n n2 = np.cross(b2, b3)\n\n # Normalize\n n1 = n1 / np.linalg.norm(n1)\n n2 = n2 / np.linalg.norm(n2)\n\n # Rotation direction\n m1 = np.cross(n1, b2 / np.linalg.norm(b2))\n\n # Calculate angle\n x = np.dot(n1, n2)\n y = np.dot(m1, n2)\n\n angle = np.arctan2(y, x) * 180 / np.pi\n dihedrals.append(angle)\n\n return dihedrals\n\n\ndef bond_angle_regularity(angles):\n \"\"\"\n Evaluates the regularity of bond angles.\n\n Parameters\n ----------\n angles : list\n List of bond angles in degrees\n\n Returns\n -------\n float\n Regularity score (range: 0 to 1)\n \"\"\"\n if not angles:\n return 0.5\n\n std_dev = np.std(angles)\n score = max(0.0, min(1.0, 1.0 - std_dev / 20.0))\n return score\n\ndef dihedral_distribution_score(dihedrals):\n \"\"\"\n Evaluates the distribution of dihedral angles.\n\n Parameters\n ----------\n dihedrals : list\n List of dihedral angles in degrees\n\n Returns\n -------\n float\n Distribution quality score (range: 0 to 1)\n \"\"\"\n if not dihedrals:\n return 0.5\n\n # Bin the angles into 12 bins (30 degrees each)\n bins = np.linspace(-180, 180, 13)\n hist, _ = np.histogram(dihedrals, bins=bins)\n\n distribution = hist / np.sum(hist)\n\n # Calculate entropy of the distribution\n entropy = -np.sum([p * np.log(p + 1e-10) for p in distribution])\n max_entropy = np.log(len(bins) - 1)\n normalized_entropy = entropy / max_entropy\n\n # Score based on entropy (moderate entropy preferred)\n if normalized_entropy < 0.4:\n score = normalized_entropy / 0.4\n elif normalized_entropy > 0.8:\n score = 1.0 - (normalized_entropy - 0.8) / 0.2\n else:\n score = 1.0\n\n return score\n\ndef detect_local_clashes(structure):\n \"\"\"\n Detects and scores local steric clashes in the structure.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure (shape: N x 3)\n\n Returns\n -------\n float\n Clash score based on absence of steric overlaps (range: 0 to 1)\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n valid_coords = structure[valid_mask]\n\n if len(valid_coords) < 4:\n return 0.5\n\n # Approximate van der Waals radius for RNA nucleotides\n vdw_radius = 3.5 # Ångströms\n clash_threshold = vdw_radius * 0.7\n\n n_clashes = 0\n n_pairs = 0\n\n for i in range(len(valid_coords)):\n for j in range(i + 3, len(valid_coords)): # Skip close sequential residues\n n_pairs += 1\n dist = np.linalg.norm(valid_coords[i] - valid_coords[j])\n if dist < clash_threshold:\n n_clashes += 1\n\n if n_pairs == 0:\n return 1.0 # No pairs to check\n\n clash_ratio = n_clashes / n_pairs\n score = max(0.0, 1.0 - clash_ratio * 20.0)\n return score\n\ndef calculate_comprehensive_quality(structure, seq_length=None, reference=None):\n \"\"\"\n Calculates a comprehensive quality score for an RNA structure,\n with metrics adapted to the structure's size.\n\n Parameters:\n -----------\n structure: numpy.ndarray\n Structure coordinates (shape: N x 3)\n seq_length: int, optional\n Sequence length (automatically determined if None)\n reference: numpy.ndarray, optional\n Reference structure for comparative metrics\n\n Returns:\n --------\n dict: Detailed quality metrics\n \"\"\"\n # Determine size if not provided\n if seq_length is None:\n valid_mask = ~np.all(structure == 0, axis=1)\n seq_length = np.sum(valid_mask)\n\n # Determine size category\n if seq_length < 50:\n size_category = 'small'\n elif seq_length < 200:\n size_category = 'medium'\n else:\n size_category = 'large'\n\n # Compute basic metrics\n energy = calculate_energy(structure)\n compactness = calculate_compactness(structure)\n\n # Compute size-specific metrics\n size_specific_metrics = {}\n\n if size_category == 'small':\n # For small structures, detailed local metrics are important\n bond_angles = calculate_bond_angles(structure)\n dihedral_angles = calculate_dihedral_angles(structure)\n\n size_specific_metrics = {\n 'bond_angle_regularity': bond_angle_regularity(bond_angles),\n 'dihedral_distribution': dihedral_distribution_score(dihedral_angles),\n 'local_clash_score': detect_local_clashes(structure)\n }\n\n elif size_category == 'medium':\n # For medium structures, a balance between local and global metrics\n secondary_structure = predict_secondary_structure(structure)\n packing_density = calculate_packing_density(structure)\n\n size_specific_metrics = {\n 'secondary_structure_score': evaluate_secondary_structure(secondary_structure),\n 'packing_density': packing_density,\n 'contact_order': calculate_contact_order(structure)\n }\n\n else: # large\n # For large structures, domain and global metrics are crucial\n domains = identify_domains(structure)\n\n size_specific_metrics = {\n 'domain_separation': evaluate_domain_separation(domains),\n 'tertiary_packing': evaluate_tertiary_packing(structure),\n 'long_range_contact_score': evaluate_long_range_contacts(structure)\n }\n\n # Reference-based metrics, if available\n comparative_metrics = {}\n if reference is not None:\n rmsd_val = rmsd(structure, reference)\n gdt_ts = calculate_gdt_ts(structure, reference)\n lddt = calculate_lddt(structure, reference)\n\n comparative_metrics = {\n 'rmsd': rmsd_val,\n 'gdt_ts': gdt_ts,\n 'lddt': lddt\n }\n\n # Combine all metrics\n all_metrics = {\n 'size_category': size_category,\n 'energy': energy,\n 'compactness': compactness,\n **size_specific_metrics,\n **comparative_metrics\n }\n\n # Compute final weighted score based on size category\n weights = {\n 'small': {\n 'energy': 0.3,\n 'compactness': 0.25,\n 'bond_angle_regularity': 0.2,\n 'dihedral_distribution': 0.15,\n 'local_clash_score': 0.1\n },\n 'medium': {\n 'energy': 0.35,\n 'compactness': 0.2,\n 'secondary_structure_score': 0.2,\n 'packing_density': 0.15,\n 'contact_order': 0.1\n },\n 'large': {\n 'energy': 0.4,\n 'compactness': 0.15,\n 'domain_separation': 0.2,\n 'tertiary_packing': 0.15,\n 'long_range_contact_score': 0.1\n }\n }\n\n # Apply category-specific weights\n category_weights = weights.get(size_category, {})\n weighted_score = 0.0\n weight_sum = 0.0\n\n for metric, weight in category_weights.items():\n if metric in all_metrics:\n weighted_score += all_metrics[metric] * weight\n weight_sum += weight\n\n # Normalize if needed\n if weight_sum > 0:\n final_score = weighted_score / weight_sum\n else:\n # Fallback to simple average of available numeric metrics\n available_metrics = [v for k, v in all_metrics.items() \n if k != 'size_category' and isinstance(v, (int, float))]\n final_score = sum(available_metrics) / max(1, len(available_metrics))\n\n all_metrics['final_quality_score'] = final_score\n\n return all_metrics\n\ndef create_variation(base_structure, scale=0.2):\n \"\"\"\n Generate a variation of the base structure using either adaptive or fixed parameters.\n \n Parameters:\n base_structure: numpy array (n, 3) representing the base structure.\n scale: noise scale to use (equivalent to base_noise in the original function)\n \n Returns:\n A new candidate structure (numpy array, n, 3).\n \"\"\"\n # Default to adaptive mode and use scale as base_noise\n base_noise = scale\n preserve_distance = True\n use_global = False\n correlation = 0.8 + np.random.uniform(-0.1, 0.1)\n \n # Call the adaptive function\n adjusted_noise = adaptive_noise_adjustment(\n base_structure, \n base_noise=base_noise, \n ideal_length=3.8, \n lower_bound=0.05, \n divisor_mean=3.8, \n divisor_std=1.0\n )\n \n variation = sample_structural_variation(\n base_structure,\n noise_level=adjusted_noise,\n preserve_distance=preserve_distance,\n use_global_movement=use_global,\n correlation=correlation\n )\n \n # Apply additional refinement\n variation = refine_rna_backbone_with_dihedrals(variation, ideal_dihedral=180.0)\n \n # Normalize the structure\n variation = normalize_structure(variation)\n \n return variation\n\ndef stratified_cross_validation(structures, labels, size_groups, n_folds=5):\n \"\"\"\n Performs stratified cross-validation by RNA size groups.\n \n Parameters:\n -----------\n structures: list\n List of RNA structures\n labels: list\n True quality scores for each structure\n size_groups: list\n Size category for each structure ('small', 'medium', 'large')\n n_folds: int\n Number of folds for cross-validation\n \n Returns:\n --------\n dict: Cross-validation results by size group\n \"\"\"\n # Group indices by size category\n grouped_indices = {\n 'small': [],\n 'medium': [],\n 'large': []\n }\n \n for i, group in enumerate(size_groups):\n if group in grouped_indices:\n grouped_indices[group].append(i)\n \n # Results per group\n results = {\n 'small': {'adaptive': [], 'fixed': []},\n 'medium': {'adaptive': [], 'fixed': []},\n 'large': {'adaptive': [], 'fixed': []}\n }\n \n # Perform separate cross-validation for each group\n for group_name, indices in grouped_indices.items():\n if not indices:\n continue\n \n # Shuffle the indices\n random.shuffle(indices)\n \n # Split into folds\n fold_size = max(1, len(indices) // n_folds)\n folds = [indices[i:i + fold_size] for i in range(0, len(indices), fold_size)]\n \n # Ensure exactly n_folds by merging the smallest folds\n while len(folds) > n_folds:\n folds.sort(key=len)\n folds[0].extend(folds[1])\n folds.pop(1)\n \n # Run cross-validation\n for test_fold_idx in range(len(folds)):\n test_indices = folds[test_fold_idx]\n train_indices = [idx for fold_idx, fold in enumerate(folds) \n if fold_idx != test_fold_idx for idx in fold]\n \n # Extract training and testing data\n X_train = [structures[i] for i in train_indices]\n y_train = [labels[i] for i in train_indices]\n X_test = [structures[i] for i in test_indices]\n y_test = [labels[i] for i in test_indices]\n \n # Evaluate adaptive mode\n adaptive_predictions = evaluate_structures_adaptive(X_test)\n adaptive_score = calculate_metrics(adaptive_predictions, y_test)\n results[group_name]['adaptive'].append(adaptive_score)\n \n # Evaluate fixed mode\n fixed_predictions = evaluate_structures_fixed(X_test)\n fixed_score = calculate_metrics(fixed_predictions, y_test)\n results[group_name]['fixed'].append(fixed_score)\n \n # Compute averages by group and mode\n summary = {}\n for group in results:\n summary[group] = {}\n for mode in results[group]:\n if results[group][mode]:\n avg_score = sum(results[group][mode]) / len(results[group][mode])\n std_score = (sum((x - avg_score) ** 2 for x in results[group][mode]) \n / len(results[group][mode])) ** 0.5\n summary[group][mode] = {\n 'mean_score': avg_score,\n 'std_score': std_score,\n 'individual_scores': results[group][mode]\n }\n \n return summary\n\ndef evaluate_segment_quality(segment):\n \"\"\"\n Evaluates the quality of a specific structural segment.\n\n Parameters\n ----------\n segment : numpy.ndarray\n 3D coordinates of the structural segment to be evaluated\n\n Returns\n -------\n float\n Quality score of the segment (range: 0 to 1)\n \"\"\"\n if len(segment) < 3:\n return 0.5 # Neutral score for very short segments\n\n # Compute local metrics for the segment\n\n # 1. Local energy based on bond lengths\n bond_lengths = []\n for i in range(1, len(segment)):\n if not (np.all(segment[i] == 0) or np.all(segment[i - 1] == 0)):\n bond_lengths.append(np.linalg.norm(segment[i] - segment[i - 1]))\n\n if not bond_lengths:\n return 0.5\n\n avg_bond_length = np.mean(bond_lengths)\n bond_std = np.std(bond_lengths)\n\n # Score based on closeness to ideal bond length (3.8 Å)\n bond_score = 1.0 - min(1.0, abs(avg_bond_length - 3.8) / 3.8)\n\n # 2. Bond angle regularity (if at least 3 valid points)\n angle_score = 0.5\n if len(segment) >= 3:\n angles = []\n for i in range(1, len(segment) - 1):\n if not (np.all(segment[i] == 0) or np.all(segment[i - 1] == 0) or np.all(segment[i + 1] == 0)):\n v1 = segment[i - 1] - segment[i]\n v2 = segment[i + 1] - segment[i]\n\n norm_v1 = np.linalg.norm(v1)\n norm_v2 = np.linalg.norm(v2)\n\n if norm_v1 > 0 and norm_v2 > 0:\n v1 = v1 / norm_v1\n v2 = v2 / norm_v2\n\n cos_angle = np.clip(np.dot(v1, v2), -1.0, 1.0)\n angle = np.arccos(cos_angle) * 180 / np.pi\n angles.append(angle)\n\n if angles:\n angle_std = np.std(angles)\n angle_score = max(0.0, min(1.0, 1.0 - angle_std / 20.0))\n\n # 3. Local clash detection\n clash_score = 1.0\n if len(segment) >= 4:\n clash_count = 0\n valid_mask = ~np.all(segment == 0, axis=1)\n valid_coords = segment[valid_mask]\n\n if len(valid_coords) >= 4:\n vdw_radius = 3.5 # Ångströms\n clash_threshold = vdw_radius * 0.7\n\n for i in range(len(valid_coords)):\n for j in range(i + 3, len(valid_coords)):\n dist = np.linalg.norm(valid_coords[i] - valid_coords[j])\n if dist < clash_threshold:\n clash_count += 1\n\n clash_ratio = clash_count / max(1, (len(valid_coords) * (len(valid_coords) - 3) // 2))\n clash_score = max(0.0, 1.0 - clash_ratio * 10.0)\n\n # 4. Local compactness\n compactness_score = 0.5\n if len(segment) >= 4:\n valid_mask = ~np.all(segment == 0, axis=1)\n valid_coords = segment[valid_mask]\n\n if len(valid_coords) >= 4:\n centroid = np.mean(valid_coords, axis=0)\n distances = np.linalg.norm(valid_coords - centroid, axis=1)\n rg = np.mean(distances)\n\n # Expected radius of gyration for a well-folded RNA segment\n expected_rg = 3.5 * (len(valid_coords) ** 0.33)\n compactness_score = max(0.0, min(1.0, 1.0 - abs(rg - expected_rg) / expected_rg))\n\n # Combine scores with appropriate weights\n final_score = (\n 0.4 * bond_score +\n 0.3 * angle_score +\n 0.2 * clash_score +\n 0.1 * compactness_score\n )\n\n return final_score\n\ndef calculate_residue_confidence(structure, all_structures):\n \"\"\"\n Calculates per-residue confidence based on consistency across multiple structures.\n\n Parameters\n ----------\n structure : numpy.ndarray\n Reference structure (3D coordinates)\n\n all_structures : list of tuples\n List of (structure, score) pairs, ordered by quality\n\n Returns\n -------\n list\n Confidence score for each residue (range: 0 to 1)\n \"\"\"\n if not all_structures or len(all_structures) < 2:\n return [0.5] * len(structure) # Not enough structures for comparison\n\n # Extract only the structures from the (structure, score) pairs\n structures = [s[0] for s in all_structures]\n\n n_residues = len(structure)\n confidence = np.zeros(n_residues)\n\n # Calculate positional variability of each residue across structures\n for i in range(n_residues):\n if np.all(structure[i] == 0):\n confidence[i] = 0.0\n continue\n\n # Collect valid coordinates for this residue from all structures\n coords = []\n for s in structures:\n if i < len(s) and not np.all(s[i] == 0):\n coords.append(s[i])\n\n if len(coords) < 2:\n confidence[i] = 0.3 # Low confidence due to insufficient data\n continue\n\n coords = np.array(coords)\n mean_coord = np.mean(coords, axis=0)\n rmsd = np.sqrt(np.mean(np.sum((coords - mean_coord) ** 2, axis=1)))\n\n # Convert RMSD to confidence score (lower RMSD = higher confidence)\n # Typically, RMSDs < 2 Å indicate high confidence\n if rmsd < 1.0:\n confidence[i] = 1.0\n elif rmsd > 5.0:\n confidence[i] = 0.0\n else:\n confidence[i] = max(0.0, 1.0 - (rmsd - 1.0) / 4.0)\n\n return confidence\n\ndef evaluate_residue_quality(structure):\n \"\"\"\n Evaluates the individual quality of each residue in the structure.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure\n\n Returns\n -------\n list\n Quality score for each residue (range: 0 to 1)\n \"\"\"\n n_residues = len(structure)\n quality = np.zeros(n_residues)\n\n # Assign 0.0 to invalid residues\n valid_mask = ~np.all(structure == 0, axis=1)\n quality[~valid_mask] = 0.0\n\n for i in range(n_residues):\n if not valid_mask[i]:\n continue\n\n # 1. Bond length (for non-terminal residues)\n bond_score = 0.5\n if i > 0 and valid_mask[i - 1]:\n dist = np.linalg.norm(structure[i] - structure[i - 1])\n bond_score = 1.0 - min(1.0, abs(dist - 3.8) / 3.8)\n\n # 2. Bond angle (for residues with neighbors)\n angle_score = 0.5\n if i > 0 and i < n_residues - 1 and valid_mask[i - 1] and valid_mask[i + 1]:\n v1 = structure[i - 1] - structure[i]\n v2 = structure[i + 1] - structure[i]\n\n norm_v1 = np.linalg.norm(v1)\n norm_v2 = np.linalg.norm(v2)\n\n if norm_v1 > 0 and norm_v2 > 0:\n v1 = v1 / norm_v1\n v2 = v2 / norm_v2\n\n cos_angle = np.clip(np.dot(v1, v2), -1.0, 1.0)\n angle = np.arccos(cos_angle) * 180 / np.pi\n\n # Ideal RNA bond angles are close to 110°\n angle_score = 1.0 - min(1.0, abs(angle - 110) / 50)\n\n # 3. Local packing environment\n packing_score = 0.5\n contact_count = 0\n\n for j in range(max(0, i - 10), min(n_residues, i + 11)):\n if abs(i - j) > 2 and valid_mask[j]:\n dist = np.linalg.norm(structure[i] - structure[j])\n if dist < 8.0:\n contact_count += 1\n\n contact_density = contact_count / min(20, n_residues)\n if contact_density < 0.1:\n packing_score = 0.3\n elif contact_density > 0.5:\n packing_score = 0.9\n else:\n packing_score = 0.3 + contact_density * 1.2\n\n # Combine the metrics with appropriate weights\n quality[i] = 0.3 * bond_score + 0.3 * angle_score + 0.4 * packing_score\n\n return quality\n\ndef identify_low_quality_regions(residue_quality, threshold=0.4, min_region_size=3):\n \"\"\"\n Identifies continuous low-quality regions in the structure.\n\n Parameters\n ----------\n residue_quality : list\n Per-residue quality scores\n\n threshold : float\n Threshold below which quality is considered low\n\n min_region_size : int\n Minimum size of a region to be considered\n\n Returns\n -------\n list of tuples\n (start, end) pairs representing low-quality regions\n \"\"\"\n n_residues = len(residue_quality)\n regions = []\n\n i = 0\n while i < n_residues:\n if residue_quality[i] < threshold:\n start = i\n while i < n_residues and residue_quality[i] < threshold:\n i += 1\n end = i\n if end - start >= min_region_size:\n regions.append((start, end))\n else:\n i += 1\n\n return regions\n\ndef refine_structure_continuity(structure):\n \"\"\"\n Refines the structure to ensure continuity of the main chain.\n\n Parameters\n ----------\n structure : numpy.ndarray\n 3D coordinates of the structure\n\n Returns\n -------\n numpy.ndarray\n Refined structure with adjusted backbone continuity\n \"\"\"\n valid_mask = ~np.all(structure == 0, axis=1)\n\n # If structure is too small or has no valid residues, return as is\n if np.sum(valid_mask) < 2:\n return structure\n\n refined = structure.copy()\n\n # Detect discontinuities (large jumps between consecutive residues)\n for i in range(1, len(structure)):\n if valid_mask[i] and valid_mask[i - 1]:\n dist = np.linalg.norm(structure[i] - structure[i - 1])\n\n # If the distance is much larger than expected (typically 3.8 Å)\n if dist > 5.0:\n # Interpolate to smooth the transition\n direction = structure[i] - structure[i - 1]\n direction = direction / dist\n\n # Adjust the position to be at a more typical bond distance\n refined[i] = structure[i - 1] + direction * 3.8\n\n return refined\n\ndef enhanced_merge_structures(adaptive_structures, fixed_structures, sequence_features=None, merge_strategy='auto'):\n \"\"\"\n Performs intelligent merging of structures generated by the adaptive and fixed modes,\n using the best regions from each.\n\n Parameters:\n -----------\n adaptive_structures: list\n List of structures generated by the adaptive mode\n fixed_structures: list\n List of structures generated by the fixed mode\n sequence_features: dict, optional\n Additional sequence features\n merge_strategy: str\n Strategy for merging ('segment_quality', 'regional', 'weighted', or 'auto')\n\n Returns:\n --------\n list: Optimized merged structures\n \"\"\"\n if not adaptive_structures or not fixed_structures:\n return adaptive_structures or fixed_structures\n \n # Get sequence length from first structure\n seq_length = len(adaptive_structures[0])\n \n # Determine best strategy if auto is selected\n if merge_strategy == 'auto':\n if seq_length < 50:\n merge_strategy = 'segment_quality'\n elif seq_length < 200:\n merge_strategy = 'weighted'\n else:\n merge_strategy = 'regional'\n \n # Get best structures from each approach\n best_adaptive = get_best_structures(adaptive_structures, top_k=3)\n best_fixed = get_best_structures(fixed_structures, top_k=3)\n \n merged_structures = []\n \n # Strategy 1: Segment-based quality merge\n if merge_strategy in ['segment_quality', 'auto']:\n segment_size = min(20, seq_length // 5)\n for a_struct in best_adaptive[:1]:\n for f_struct in best_fixed[:1]:\n merged = np.zeros_like(a_struct)\n for start in range(0, seq_length, segment_size):\n end = min(start + segment_size, seq_length)\n \n # Calculate segment quality using our new function\n a_quality = assess_structure_quality(a_struct[start:end])\n f_quality = assess_structure_quality(f_struct[start:end])\n \n # Choose the better segment\n if a_quality > f_quality:\n merged[start:end] = a_struct[start:end]\n else:\n merged[start:end] = f_struct[start:end]\n \n merged_structures.append(merged)\n \n # Strategy 2: Regional merge based on flexibility\n if merge_strategy in ['regional', 'auto'] and seq_length >= 100:\n for f_struct in best_fixed[:2]:\n for a_struct in best_adaptive[:1]:\n # Identify flexible regions using our new function\n flexible_indices = identify_flexible_regions(f_struct)\n \n # Convert indices to regions (start, end)\n flexible_regions = []\n if flexible_indices:\n start = flexible_indices[0]\n prev = start\n for idx in flexible_indices[1:]:\n if idx > prev + 3: # Gap of more than 3 residues\n flexible_regions.append((start, prev))\n start = idx\n prev = idx\n flexible_regions.append((start, prev)) # Add last region\n \n # If no flexible regions found, use regular segmentation\n if not flexible_regions:\n segment_size = min(20, seq_length // 5)\n flexible_regions = [(start, min(start + segment_size - 1, seq_length - 1)) \n for start in range(0, seq_length, segment_size)]\n \n # Merge regions\n merged = merge_regions(f_struct, a_struct, flexible_regions)\n merged_structures.append(merged)\n \n # Strategy 3: Weighted average merge\n if merge_strategy in ['weighted', 'auto']:\n for a_struct in best_adaptive[:2]:\n for f_struct in best_fixed[:2]:\n # Calculate quality scores\n a_quality = assess_structure_quality(a_struct)\n f_quality = assess_structure_quality(f_struct)\n total_quality = a_quality + f_quality\n \n if total_quality > 0:\n a_weight = a_quality / total_quality\n f_weight = f_quality / total_quality\n merged = weighted_merge(a_struct, f_struct, a_weight, f_weight)\n merged_structures.append(merged)\n \n # Ensure we have enough structures\n if not merged_structures:\n merged_structures = best_adaptive[:2] + best_fixed[:2]\n \n # Add top original structures for diversity\n merged_structures.extend(best_adaptive[:2])\n merged_structures.extend(best_fixed[:2])\n \n # Remove duplicates and refine\n unique_structures = remove_duplicate_structures(merged_structures)\n \n # Refine structures for continuity (ensure bonds are reasonable)\n refined_structures = []\n for structure in unique_structures:\n # Simple refinement: ensure bond lengths are reasonable\n refined = structure.copy()\n bonds = np.linalg.norm(np.diff(refined, axis=0), axis=1)\n \n # Identify and fix unreasonable bonds\n for i in range(len(bonds)):\n if bonds[i] > 5.0 or bonds[i] < 2.0: # Unreasonable bond length\n if i > 0 and i < len(refined) - 2:\n # Average with neighboring residues\n refined[i+1] = (refined[i] + refined[i+2]) / 2.0\n \n refined_structures.append(refined)\n \n # Get the best structures using our quality assessment\n # We're using calculate_comprehensive_quality instead of get_best_structures \n # to take advantage of the more detailed quality metrics\n scored_structures = []\n for structure in refined_structures:\n quality = calculate_comprehensive_quality(structure, seq_length)\n scored_structures.append((quality[\"quality_score\"], structure))\n \n # Sort by quality score (descending)\n scored_structures.sort(key=lambda x: x[0], reverse=True)\n \n # Return the top 5 structures\n return [s for _, s in scored_structures[:5]]","metadata":{"execution":{"iopub.status.busy":"2025-04-05T16:08:10.892524Z","iopub.execute_input":"2025-04-05T16:08:10.892759Z","iopub.status.idle":"2025-04-05T16:08:11.027839Z","shell.execute_reply.started":"2025-04-05T16:08:10.892741Z","shell.execute_reply":"2025-04-05T16:08:11.027047Z"},"papermill":{"duration":0.033009,"end_time":"2025-03-26T03:43:08.66035","exception":false,"start_time":"2025-03-26T03:43:08.627341","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def 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)\n \ndef 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-04-05T16:08:11.02864Z","iopub.execute_input":"2025-04-05T16:08:11.028938Z","iopub.status.idle":"2025-04-05T16:08:11.045513Z","shell.execute_reply.started":"2025-04-05T16:08:11.02891Z","shell.execute_reply":"2025-04-05T16:08:11.044656Z"},"papermill":{"duration":0.033522,"end_time":"2025-03-26T03:43:08.712388","exception":false,"start_time":"2025-03-26T03:43:08.678866","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Phase 5: Submission Creation","metadata":{"papermill":{"duration":0.017635,"end_time":"2025-03-26T03:43:08.747861","exception":false,"start_time":"2025-03-26T03:43:08.730226","status":"completed"},"tags":[]}},{"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-04-05T16:08:11.04633Z","iopub.execute_input":"2025-04-05T16:08:11.046617Z","iopub.status.idle":"2025-04-05T16:08:11.067866Z","shell.execute_reply.started":"2025-04-05T16:08:11.046585Z","shell.execute_reply":"2025-04-05T16:08:11.067315Z"},"papermill":{"duration":0.029933,"end_time":"2025-03-26T03:43:08.79562","exception":false,"start_time":"2025-03-26T03:43:08.765687","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Main Pipeline Functions","metadata":{"papermill":{"duration":0.017442,"end_time":"2025-03-26T03:43:08.831104","exception":false,"start_time":"2025-03-26T03:43:08.813662","status":"completed"},"tags":[]}},{"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 seed_attempts=200, \n optimal_params={'noise': 0.21, 'corr': 0.83},\n ab_test=True\n):\n \"\"\"\n Run a hybrid pipeline that combines golden seeds approach with NN pruning.\n \n Parameters:\n -----------\n X_valid, y_valid: Validation data for training models\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 seed_attempts: Number of seeds to try\n optimal_params: Optimal parameters for the reference model\n ab_test: Whether to run A/B testing between 'adaptive' and 'fixed' modes\n \n Returns:\n --------\n submission_df, status_dict\n \"\"\"\n print(\"=\" * 80)\n print(\"HYBRID PIPELINE: GOLDEN SEEDS + NN PRUNING\".center(80))\n print(\"=\" * 80)\n \n # Initialize parameter optimizer for continuous learning\n parameter_optimizer = ParameterOptimizer(history_file=\"parameter_history.json\")\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 'ab_test_results': {}\n }\n \n try:\n # PHASE 1: Find Golden Seeds\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\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(top_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 # Calculate GC content for mode selection\n gc_content = sum(1 for n in seq if n in 'GC') / len(seq)\n sequence_features = {'gc_content': gc_content, 'length': seq_length}\n \n # Determine optimal noise level for this sequence\n optimal_noise = determine_optimal_noise(seq, seq_length, gc_content)\n \n if ab_test:\n # A/B TEST: Comparison between 'adaptive' and 'fixed' modes\n print(\"\\n=== A/B Test: ADAPTIVE Mode ===\")\n # Get optimized parameters from history for adaptive mode\n adaptive_params = parameter_optimizer.suggest_parameters(\n get_size_category(seq_length), 'adaptive')\n \n # Generate candidates using 'adaptive' mode only\n candidates_adaptive = generate_diverse_structures_from_bases(\n base_structures,\n seq_length,\n quality_model,\n num_per_base=5,\n mode='adaptive' # Only passing the mode parameter\n )\n print(f\"Adaptive mode: {len(candidates_adaptive)} candidates generated.\")\n \n # Evaluate candidates by size group (small, medium, large)\n print(\"Candidate evaluation (adaptive) by size group:\")\n metrics_adaptive = evaluate_candidates_by_group(candidates_adaptive)\n for group, metrics in metrics_adaptive.items():\n print(f\"Group {group}: {metrics}\")\n \n print(\"\\n=== A/B Test: FIXED Mode ===\")\n # Get optimized parameters from history for fixed mode\n fixed_params = parameter_optimizer.suggest_parameters(\n get_size_category(seq_length), 'fixed')\n \n # Generate candidates using 'fixed' mode only\n candidates_fixed = generate_diverse_structures_from_bases(\n base_structures,\n seq_length,\n quality_model,\n num_per_base=5,\n mode='fixed' # Only passing the mode parameter\n )\n print(f\"Fixed mode: {len(candidates_fixed)} candidates generated.\")\n \n # Evaluate candidates by size group\n print(\"Candidate evaluation (fixed) by size group:\")\n metrics_fixed = evaluate_candidates_by_group(candidates_fixed)\n for group, metrics in metrics_fixed.items():\n print(f\"Group {group}: {metrics}\")\n \n # Store A/B test results in status\n status['ab_test_results'] = {\n 'adaptive': {\n 'num_candidates': len(candidates_adaptive),\n 'metrics_by_group': metrics_adaptive\n },\n 'fixed': {\n 'num_candidates': len(candidates_fixed),\n 'metrics_by_group': metrics_fixed\n }\n }\n \n # Choose the best mode based on results\n adaptive_score = sum([m.get('quality_score', 0) for m in metrics_adaptive.values()])\n fixed_score = sum([m.get('quality_score', 0) for m in metrics_fixed.values()])\n \n # Determine optimal merge strategy based on sequence size\n if seq_length < 50:\n merge_strategy = 'segment_quality'\n elif seq_length < 200:\n merge_strategy = 'weighted'\n else:\n merge_strategy = 'regional'\n \n print(f\"Merging structures from both modes using {merge_strategy} strategy...\")\n \n # Implement multi-layered approach with improved StructureWrapper\n combined_candidates = None\n \n # APPROACH 1: Try the enhanced_merge_structures with improved StructureWrapper\n try:\n print(\"Attempting enhanced structure merging...\")\n # Convert candidates to format compatible with enhanced_merge_structures\n fixed_scored_wrappers = []\n adaptive_scored_wrappers = []\n \n # Get quality scores for each candidate if possible\n adaptive_quality_scores = []\n fixed_quality_scores = []\n \n # Try to get quality scores for adaptive candidates\n try:\n for candidate in candidates_adaptive:\n if 'calculate_comprehensive_quality' in globals():\n quality_metrics = calculate_comprehensive_quality(candidate, seq_length)\n quality_score = quality_metrics.get('quality_score', 0.5)\n else:\n # Use neural network if available\n if status['nn_training_success'] and quality_model is not None:\n # This depends on your actual implementation of quality prediction\n quality_score = 0.5 # Default if prediction not possible\n else:\n quality_score = 0.5\n adaptive_quality_scores.append(quality_score)\n except Exception as e:\n print(f\"Error getting quality scores for adaptive candidates: {e}\")\n # Use default scores\n adaptive_quality_scores = [0.5] * len(candidates_adaptive)\n \n # Try to get quality scores for fixed candidates\n try:\n for candidate in candidates_fixed:\n if 'calculate_comprehensive_quality' in globals():\n quality_metrics = calculate_comprehensive_quality(candidate, seq_length)\n quality_score = quality_metrics.get('quality_score', 0.5)\n else:\n # Use neural network if available\n if status['nn_training_success'] and quality_model is not None:\n # This depends on your actual implementation of quality prediction\n quality_score = 0.5 # Default if prediction not possible\n else:\n quality_score = 0.5\n fixed_quality_scores.append(quality_score)\n except Exception as e:\n print(f\"Error getting quality scores for fixed candidates: {e}\")\n # Use default scores\n fixed_quality_scores = [0.5] * len(candidates_fixed)\n \n # Create wrapper objects with quality scores\n for i, candidate in enumerate(candidates_adaptive):\n adaptive_scored_wrappers.append(\n StructureWrapper(candidate, quality_score=adaptive_quality_scores[i])\n )\n \n for i, candidate in enumerate(candidates_fixed):\n fixed_scored_wrappers.append(\n StructureWrapper(candidate, quality_score=fixed_quality_scores[i])\n )\n \n # Try using the enhanced_merge_structures function with our wrappers\n merged_candidates = enhanced_merge_structures(\n adaptive_scored_wrappers, \n fixed_scored_wrappers, \n sequence_features,\n merge_strategy=merge_strategy\n )\n \n # Extract the original structures from the wrappers\n combined_candidates = [c.structure for c in merged_candidates]\n print(\"Enhanced structure merging successful!\")\n except Exception as e:\n print(f\"Enhanced structure merging failed: {e}\")\n combined_candidates = None\n \n # APPROACH 2: If enhanced merge failed, try custom merging strategy\n if combined_candidates is None:\n try:\n print(\"Falling back to custom merging strategy...\")\n # Use the predefined quality scores from the metrics\n adaptive_quality = sum([m.get('quality_score', 0) for m in metrics_adaptive.values()])\n fixed_quality = sum([m.get('quality_score', 0) for m in metrics_fixed.values()])\n \n # Weight the number of candidates from each approach based on quality\n total_quality = adaptive_quality + fixed_quality\n if total_quality > 0:\n adaptive_ratio = adaptive_quality / total_quality\n else:\n adaptive_ratio = 0.5 # Equal if no quality data\n \n # Calculate how many to take from each\n target_total = 10 # We want 10 total candidates\n num_adaptive = max(1, min(len(candidates_adaptive), int(target_total * adaptive_ratio)))\n num_fixed = max(1, min(len(candidates_fixed), target_total - num_adaptive))\n \n # Get top candidates from each approach based on evaluation\n try:\n if status['nn_training_success'] and quality_model is not None:\n adaptive_top = evaluate_and_prune_structures(\n candidates_adaptive, seq_features, quality_model, top_k=num_adaptive\n )\n fixed_top = evaluate_and_prune_structures(\n candidates_fixed, seq_features, quality_model, top_k=num_fixed\n )\n else:\n # Use rule-based evaluation\n adaptive_top = evaluate_and_prune_rules(candidates_adaptive, top_k=num_adaptive)\n fixed_top = evaluate_and_prune_rules(candidates_fixed, top_k=num_fixed)\n except Exception as inner_e:\n print(f\"Error in evaluation for custom merge: {inner_e}\")\n # Just take the first n structures\n adaptive_top = candidates_adaptive[:num_adaptive]\n fixed_top = candidates_fixed[:num_fixed]\n \n # Combine the top selections\n combined_candidates = list(adaptive_top) + list(fixed_top)\n print(f\"Custom merging successful! Selected {len(combined_candidates)} candidates.\")\n except Exception as e:\n print(f\"Custom merging strategy failed: {e}\")\n combined_candidates = None\n \n # APPROACH 3: If both custom and enhanced merging failed, just combine the candidates\n if combined_candidates is None:\n print(\"Falling back to simple combination of candidates...\")\n combined_candidates = list(candidates_adaptive) + list(candidates_fixed)\n # Limit to 10 candidates if there are too many\n if len(combined_candidates) > 10:\n combined_candidates = combined_candidates[:10]\n \n # Final candidates\n candidates = combined_candidates\n \n # Register results for continuous learning\n for group, metrics in metrics_adaptive.items():\n parameter_optimizer.record_result(\n {'mode': 'adaptive', **adaptive_params},\n group,\n 'adaptive',\n metrics.get('quality_score', 0)\n )\n \n for group, metrics in metrics_fixed.items():\n parameter_optimizer.record_result(\n {'mode': 'fixed', **fixed_params},\n group,\n 'fixed',\n metrics.get('quality_score', 0)\n )\n \n if adaptive_score > fixed_score:\n print(\"A/B Test Result: ADAPTIVE mode performed better\")\n else:\n print(\"A/B Test Result: FIXED mode performed better\")\n else:\n # Use optimal mode and parameters based on sequence characteristics\n mode, params = select_optimal_mode_and_params(\n seq_length, \n sequence_features,\n parameter_optimizer\n )\n print(f\"Using {mode} mode for this sequence based on its characteristics\")\n \n # Extract parameters for use elsewhere if needed\n merge_strategy = params.get('merge_strategy', 'auto') if isinstance(params, dict) else 'auto'\n \n # Generate candidates using only the optimal mode, without additional parameters\n candidates = generate_diverse_structures_from_bases(\n base_structures,\n seq_length,\n quality_model,\n num_per_base=5,\n mode=mode # Only passing the mode parameter\n )\n \n # Evaluate and prune candidates - implement all three approaches with fallback\n print(\"Evaluating candidates with multi-layer fallback strategy...\")\n \n # Try all three evaluation approaches in sequence, fallback if one fails\n top_structures = None\n \n # APPROACH 1: Try comprehensive quality metrics first (most sophisticated)\n if top_structures is None:\n try:\n print(\"Attempting comprehensive quality assessment...\")\n quality_scores = []\n for candidate in candidates:\n # Safely calculate comprehensive quality\n try:\n quality_metrics = calculate_comprehensive_quality(candidate, seq_length)\n score = quality_metrics.get('quality_score', 0.0)\n # Handle NaN or inf values\n if np.isnan(score) or np.isinf(score):\n score = 0.0\n except Exception as inner_e:\n print(f\"Error calculating quality for a candidate: {inner_e}\")\n score = 0.0\n quality_scores.append(score)\n \n # Only proceed if we have valid scores\n if quality_scores and any(s > 0 for s in quality_scores):\n # Convert quality_scores to numpy array for safe indexing\n quality_scores = np.array(quality_scores)\n # Select top structures (handle ties by using stable sort)\n sort_indices = np.argsort(-quality_scores, kind='stable') # Negative for descending\n top_indices = sort_indices[:5] # Take top 5\n top_structures = [candidates[i] for i in top_indices]\n print(f\"Selected top 5 structures with comprehensive quality scores\")\n else:\n print(\"No valid quality scores found, falling back...\")\n top_structures = None\n except Exception as e:\n print(f\"Comprehensive quality assessment failed: {e}\")\n top_structures = None\n \n # APPROACH 2: Try neural network assessment next if comprehensive failed\n if top_structures is None and status['nn_training_success']:\n try:\n print(\"Falling back to neural network quality assessment...\")\n top_structures = evaluate_and_prune_structures(\n candidates, \n seq_features, \n quality_model, \n top_k=5\n )\n print(\"Neural network assessment successful!\")\n except Exception as e:\n print(f\"Neural network assessment failed: {e}\")\n top_structures = None\n \n # APPROACH 3: If both failed, use rule-based assessment as final backup\n if top_structures is None:\n try:\n print(\"Falling back to rule-based quality assessment...\")\n top_structures = evaluate_and_prune_rules(candidates, top_k=5)\n print(\"Rule-based assessment successful!\")\n except Exception as e:\n print(f\"Rule-based assessment failed: {e}\")\n # Last resort: just take first 5 candidates\n print(\"Using first 5 candidates as last resort...\")\n top_structures = candidates[:min(5, len(candidates))]\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 # Save parameter history for future runs\n parameter_optimizer.save_history(\"parameter_history.json\")\n \n # Set success\n status['success'] = True\n \n # Run sensitivity analysis on parameters if requested\n if ab_test:\n print(\"\\nPHASE 6: Parameter Sensitivity Analysis...\")\n try:\n parameter_ranges = {\n 'divisor_mean': (3.0, 4.5),\n 'divisor_std': (0.5, 1.5),\n 'noise_base': (0.01, 0.3),\n 'correlation': (0.7, 0.95)\n }\n \n def evaluate_param_set(params):\n # Simple evaluation function for sensitivity analysis\n # In practice, this would use more complex metrics\n return np.mean([m.get('quality_score', 0) for m in metrics_adaptive.values()])\n \n sensitivity_results = sensitivity_analysis(parameter_ranges, evaluate_param_set, n_samples=50)\n print(\"Parameter importance ranking:\")\n for param, importance in sorted(sensitivity_results['importances'].items(), \n key=lambda x: x[1], reverse=True):\n print(f\" {param}: {importance:.4f}\")\n \n print(\"\\nBest parameter combinations found:\")\n for i, (params, score) in enumerate(sensitivity_results['top_samples'][:3]):\n print(f\" #{i+1}: {params} (score: {score:.4f})\")\n \n # Store sensitivity results in status\n status['sensitivity_analysis'] = sensitivity_results\n except Exception as e:\n print(f\"Error in sensitivity analysis: {str(e)}\")\n print(\"Continuing without sensitivity analysis.\")\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# Custom merge function that works with structure score tuples instead of expecting quality attributes\ndef custom_merge_structures(adaptive_scored_structures, fixed_scored_structures, sequence_features, merge_strategy='auto'):\n \"\"\"\n Custom function to merge structures from adaptive and fixed modes.\n \n Parameters:\n -----------\n adaptive_scored_structures: List of (score, structure) tuples from adaptive mode\n fixed_scored_structures: List of (score, structure) tuples from fixed mode\n sequence_features: Dict with sequence features like gc_content and length\n merge_strategy: Strategy for merging ('segment_quality', 'weighted', 'regional', or 'auto')\n \n Returns:\n --------\n List of merged structures\n \"\"\"\n print(f\"Using custom merge with strategy: {merge_strategy}\")\n \n # If no structures from one mode, return the other\n if not adaptive_scored_structures:\n return [s[1] for s in fixed_scored_structures]\n if not fixed_scored_structures:\n return [s[1] for s in adaptive_scored_structures]\n \n # Extract scores and structures\n adaptive_scores = [s[0] for s in adaptive_scored_structures]\n adaptive_structures = [s[1] for s in adaptive_scored_structures]\n fixed_scores = [s[0] for s in fixed_scored_structures]\n fixed_structures = [s[1] for s in fixed_scored_structures]\n \n # Auto-select merge strategy based on sequence length if 'auto'\n if merge_strategy == 'auto':\n seq_length = sequence_features.get('length', 0)\n if seq_length < 50:\n merge_strategy = 'segment_quality'\n elif seq_length < 200:\n merge_strategy = 'weighted'\n else:\n merge_strategy = 'regional'\n print(f\"Auto-selected merge strategy: {merge_strategy}\")\n \n # Different merge strategies\n if merge_strategy == 'segment_quality':\n # Take top structures from each approach based on quality\n # For small sequences, quality is more important than diversity\n top_adaptive = [s for _, s in sorted(zip(adaptive_scores, adaptive_structures), reverse=True)[:3]]\n top_fixed = [s for _, s in sorted(zip(fixed_scores, fixed_structures), reverse=True)[:3]]\n return top_adaptive + top_fixed\n \n elif merge_strategy == 'weighted':\n # For medium sequences, blend approaches with more weight to higher quality\n # Take more from the approach with higher average quality\n avg_adaptive = sum(adaptive_scores) / len(adaptive_scores) if adaptive_scores else 0\n avg_fixed = sum(fixed_scores) / len(fixed_scores) if fixed_scores else 0\n \n if avg_adaptive > avg_fixed:\n # Adaptive is better, take more from it\n top_adaptive = [s for _, s in sorted(zip(adaptive_scores, adaptive_structures), reverse=True)[:5]]\n top_fixed = [s for _, s in sorted(zip(fixed_scores, fixed_structures), reverse=True)[:3]]\n else:\n # Fixed is better, take more from it\n top_adaptive = [s for _, s in sorted(zip(adaptive_scores, adaptive_structures), reverse=True)[:3]]\n top_fixed = [s for _, s in sorted(zip(fixed_scores, fixed_structures), reverse=True)[:5]]\n \n return top_adaptive + top_fixed\n \n elif merge_strategy == 'regional':\n # For large sequences, regional analysis - take best from each mode\n # and potentially combine structures regionally\n # This is a simplified implementation\n top_adaptive = [s for _, s in sorted(zip(adaptive_scores, adaptive_structures), reverse=True)[:4]]\n top_fixed = [s for _, s in sorted(zip(fixed_scores, fixed_structures), reverse=True)[:4]]\n return top_adaptive + top_fixed\n \n else:\n # Default: just combine the top half of each approach\n half_adaptive = len(adaptive_scored_structures) // 2\n half_fixed = len(fixed_scored_structures) // 2\n \n top_adaptive = [s for _, s in sorted(zip(adaptive_scores, adaptive_structures), reverse=True)[:half_adaptive]]\n top_fixed = [s for _, s in sorted(zip(fixed_scores, fixed_structures), reverse=True)[:half_fixed]]\n \n return top_adaptive + top_fixed\n\n# Helper function to determine size category based on sequence length\ndef get_size_category(seq_length):\n \"\"\"Determine the size category of a sequence based on its length\"\"\"\n if seq_length < 50:\n return 'small'\n elif seq_length < 200:\n return 'medium'\n else:\n return 'large'","metadata":{"execution":{"iopub.status.busy":"2025-04-05T16:08:11.068677Z","iopub.execute_input":"2025-04-05T16:08:11.068969Z","iopub.status.idle":"2025-04-05T16:08:11.114044Z","shell.execute_reply.started":"2025-04-05T16:08:11.068937Z","shell.execute_reply":"2025-04-05T16:08:11.113143Z"},"papermill":{"duration":0.031083,"end_time":"2025-03-26T03:43:08.88051","exception":false,"start_time":"2025-03-26T03:43:08.849427","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"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-04-05T16:08:11.115024Z","iopub.execute_input":"2025-04-05T16:08:11.115302Z","iopub.status.idle":"2025-04-05T16:08:11.13502Z","shell.execute_reply.started":"2025-04-05T16:08:11.115274Z","shell.execute_reply":"2025-04-05T16:08:11.134172Z"},"papermill":{"duration":0.02683,"end_time":"2025-03-26T03:43:08.925387","exception":false,"start_time":"2025-03-26T03:43:08.898557","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"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-04-05T16:08:11.138209Z","iopub.execute_input":"2025-04-05T16:08:11.138484Z","iopub.status.idle":"2025-04-05T16:08:11.157119Z","shell.execute_reply.started":"2025-04-05T16:08:11.138449Z","shell.execute_reply":"2025-04-05T16:08:11.156443Z"},"papermill":{"duration":0.030956,"end_time":"2025-03-26T03:43:08.97437","exception":false,"start_time":"2025-03-26T03:43:08.943414","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Main Execution Block","metadata":{"papermill":{"duration":0.053971,"end_time":"2025-03-26T03:43:09.046403","exception":false,"start_time":"2025-03-26T03:43:08.992432","status":"completed"},"tags":[]}},{"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 \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 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 # Run the 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 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(\"\\nHYBRID PIPELINE 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 else:\n print(f\"\\nPipeline failed with error: {status['error']}\")\n \n # Display output file information\n print(\"\\nOUTPUT FILES:\")\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 \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-04-05T16:08:11.158216Z","iopub.execute_input":"2025-04-05T16:08:11.158513Z","iopub.status.idle":"2025-04-05T16:18:00.560817Z","shell.execute_reply.started":"2025-04-05T16:08:11.158487Z","shell.execute_reply":"2025-04-05T16:18:00.559995Z"},"papermill":{"duration":103.685128,"end_time":"2025-03-26T03:44:52.749557","exception":false,"start_time":"2025-03-26T03:43:09.064429","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"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-04-05T16:18:00.561766Z","iopub.execute_input":"2025-04-05T16:18:00.562046Z","iopub.status.idle":"2025-04-05T16:18:00.586783Z","shell.execute_reply.started":"2025-04-05T16:18:00.562022Z","shell.execute_reply":"2025-04-05T16:18:00.586104Z"},"papermill":{"duration":0.068636,"end_time":"2025-03-26T03:44:52.857355","exception":false,"start_time":"2025-03-26T03:44:52.788719","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"papermill":{"duration":0.037766,"end_time":"2025-03-26T03:44:52.935922","exception":false,"start_time":"2025-03-26T03:44:52.898156","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null}]}