{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":87793,"databundleVersionId":11403143,"sourceType":"competition"}],"dockerImageVersionId":30919,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Ribonanza: Enhanced RNA Structure Prediction\n\n**Key Improvements:**\n1. Hybrid architecture with coordinate-specific heads\n2. Enhanced feature engineering\n3. Temporal trend incorporation\n4. Rigorous submission validation","metadata":{"_uuid":"1d34afad-4891-4b4a-935d-3bf9be92469c","_cell_guid":"e6e18754-e2ff-408b-9813-0ae3a1ba8096","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false}}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom tensorflow.keras.layers import (\n    Input, Dense, LSTM, Embedding, Bidirectional, \n    LayerNormalization, MultiHeadAttention, GlobalAveragePooling1D,\n    Dropout, Add, concatenate, Layer\n)\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau, TerminateOnNaN, ModelCheckpoint,  TensorBoard\nfrom tensorflow.keras.preprocessing.sequence import pad_sequences\nfrom sklearn.base import BaseEstimator, TransformerMixin\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder\nimport warnings\nfrom datetime import datetime\nwarnings.filterwarnings('ignore')","metadata":{"_uuid":"9fecbd1d-19b1-4c41-b078-2029519725a5","_cell_guid":"4c5778b3-29ca-4445-ba4f-201708569d78","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-03-20T08:24:29.431221Z","iopub.execute_input":"2025-03-20T08:24:29.431589Z","iopub.status.idle":"2025-03-20T08:24:29.437317Z","shell.execute_reply.started":"2025-03-20T08:24:29.431561Z","shell.execute_reply":"2025-03-20T08:24:29.436591Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class SequenceFeatureEngineer(BaseEstimator, TransformerMixin):\n    \"\"\"Creates and encodes sequence position features\"\"\"\n    def __init__(self, max_seq_len=206):\n        self.max_seq_len = max_seq_len\n        self.residues = ['A','C','G','U']\n        self.combinations = [a+b for a in self.residues for b in self.residues]\n        self.vocab = {'A':0, 'C':1, 'G':2, 'U':3}\n        self.encoder = LabelEncoder()\n        self.expected_features = None\n\n    def fit(self, X, y=None):\n        # Create temporary features for encoding\n        X_temp = X.copy()\n        X_temp['start_res'] = X_temp['sequence'].str[0]\n        X_temp['end_res'] = X_temp['sequence'].str[-1]\n        \n        # Keep all original columns except specific ones\n        base_features = [col for col in X.columns \n                        if col not in ['temporal_cutoff', 'description', 'all_sequences']]\n        \n        new_features = [\n            'start_res', 'end_res', 'sequence',\n            *[f'{r}_count' for r in self.residues],\n            *[f'{combo}_count' for combo in self.combinations]\n        ]\n        \n        self.expected_features = list(set(base_features + new_features))\n        self.encoder.fit(pd.concat([X_temp['start_res'], X_temp['end_res']]))\n        return self\n    \n    def transform(self, X):\n        df = X.copy()\n\n        # Preserve original sequence column\n        if 'sequence' not in df.columns:\n            raise ValueError(\"Missing required 'sequence' column\")\n        \n        # Create position features\n        df['start_res'] = df['sequence'].str[0]\n        df['end_res'] = df['sequence'].str[-1]\n        \n        # Encode start/end residues\n        df['start_res'] = self.encoder.transform(df['start_res'])\n        df['end_res'] = self.encoder.transform(df['end_res'])\n        \n        # Add residue counts\n        for r in self.residues:\n            df[f'{r}_count'] = df['sequence'].str.count(r).fillna(0)\n        \n        # Add combination counts\n        for combo in self.combinations:\n            df[f'{combo}_count'] = df['sequence'].str.count(combo).fillna(0)\n            \n        # Add encoded sequence\n        df['encoded'] = df['sequence'].apply(\n            lambda s: [self.vocab.get(b.upper(), 0) for b in s[:self.max_seq_len]] +\n                     [0]*(self.max_seq_len - len(s)))\n\n        # Ensure consistent feature set\n        for col in self.expected_features:\n            if col not in df.columns:\n                df[col] = 0\n                \n        return df[self.expected_features + ['encoded']]","metadata":{"_uuid":"e2631d81-64e1-421f-ba3e-8c31f98dbbf5","_cell_guid":"6a6db586-69c5-4c3f-89ee-a970ac4d0da4","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-03-20T08:24:29.438488Z","iopub.execute_input":"2025-03-20T08:24:29.438810Z","iopub.status.idle":"2025-03-20T08:24:29.454947Z","shell.execute_reply.started":"2025-03-20T08:24:29.438771Z","shell.execute_reply":"2025-03-20T08:24:29.454248Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Enhanced Feature Engineering\nclass EnhancedFeatureEngineer(SequenceFeatureEngineer):\n    def transform(self, X):\n        df = super().transform(X)\n        \n        # Add features while preserving identifiers\n        df['gc_content'] = (df['G_count'] + df['C_count']) / df['sequence'].str.len()\n        \n        # Calculate dinucleotide frequencies\n        for combo in self.combinations:\n            df[f'{combo}_freq'] = df[f'{combo}_count'] / (df['sequence'].str.len() - 1)\n            \n        # Proper structural feature calculation\n        df['stem_potential'] = df['sequence'].apply(lambda s: s.count('G') + s.count('C'))\n        df['loop_potential'] = df['sequence'].apply(lambda s: s.count('A') + s.count('U'))\n        \n        return df\n\n    def add_structural_features(self, df, sequence):\n        # Placeholder for actual structural calculations\n        df['stem_potential'] = sequence.count('G') + sequence.count('C')\n        df['loop_potential'] = sequence.count('A') + sequence.count('T')\n\n# Enhanced Data Augmentation\ndef augment_sequences(features_df, labels_df):\n    \"\"\"Proper data augmentation with coordinate list handling\"\"\"\n    # Store original feature columns\n    feature_cols = features_df.columns.tolist()\n    \n    # Ensure that the labels DataFrame has an 'ID' column\n    if 'ID' not in labels_df.columns:\n        raise ValueError(\"Labels DataFrame does not have an 'ID' column.\")\n    \n    # Extract join key from labels: extract target_id from ID (format: targetid_chainid_resid)\n    labels_df = labels_df.copy()  # Avoid modifying original DataFrame\n    labels_df['sequence_id'] = labels_df['ID'].str.split('_').apply(lambda x: '_'.join(x[:2]))\n    \n    # Debug: Print unique join keys from both DataFrames\n    # print(\"Unique target_id in features:\", features_df['target_id'].unique())\n    # print(\"Unique sequence_id in labels:\", labels_df['sequence_id'].unique())\n    \n    # Group labels by the extracted sequence_id\n    labels_processed = (\n        labels_df.groupby('sequence_id')\n        .agg({\n            'x_1': list,\n            'y_1': list,\n            'z_1': list,\n            'ID': 'first',   \n            'resname': 'first', \n            'resid': 'first' \n        })\n        .reset_index()\n    )\n    \n    # Merge using features_df.target_id and labels_processed.sequence_id\n    merged = pd.merge(\n        features_df,\n        labels_processed,\n        left_on='target_id',\n        right_on='sequence_id',\n        how='inner'\n    )\n    \n    if merged.empty:\n        raise ValueError(\"No matching records found between features and labels. \"\n                         \"Please check that your join keys (target_id and sequence_id) match.\")\n    \n    augmented = []\n    for _, row in merged.iterrows():\n        # Original entry\n        augmented.append(row.to_dict())\n        \n        # Reverse complement with coordinate reversal\n        rev_entry = row.to_dict()\n        rev_entry['sequence'] = rev_entry['sequence'][::-1].translate(str.maketrans('ACGU', 'UGCA'))\n        rev_entry['x_1'] = rev_entry['x_1'][::-1]\n        rev_entry['y_1'] = rev_entry['y_1'][::-1]\n        rev_entry['z_1'] = rev_entry['z_1'][::-1]\n        # Update both target_id and ID so that grouping later treats this as a separate sample\n        rev_entry['target_id'] = f\"{rev_entry['target_id']}_rev\"\n        rev_entry['ID'] = f\"{rev_entry['ID']}_rev\"\n        augmented.append(rev_entry)\n\n    \n    augmented_df = pd.DataFrame(augmented)\n    \n    # Ensure all original feature columns are still present\n    missing = [col for col in feature_cols if col not in augmented_df.columns]\n    if missing:\n        raise ValueError(f\"Missing columns after augmentation: {missing}\")\n    \n    return augmented_df[feature_cols], augmented_df[['x_1', 'y_1', 'z_1']]\n\n\n# Enhanced Training Configuration\ndef configure_training(model):\n    return [\n        EarlyStopping(patience=15, restore_best_weights=True),\n        ReduceLROnPlateau(factor=0.2, patience=5, min_lr=1e-6),\n        TerminateOnNaN(),\n        ModelCheckpoint('best_model.h5', save_best_only=True)\n    ]","metadata":{"_uuid":"7e3483a6-3ced-4522-855d-125cea2f04ae","_cell_guid":"07ebcb7b-3697-4165-8023-d296c81e1543","trusted":true,"collapsed":false,"execution":{"iopub.status.busy":"2025-03-20T08:24:29.456109Z","iopub.execute_input":"2025-03-20T08:24:29.456378Z","iopub.status.idle":"2025-03-20T08:24:29.476523Z","shell.execute_reply.started":"2025-03-20T08:24:29.456348Z","shell.execute_reply":"2025-03-20T08:24:29.475729Z"},"jupyter":{"outputs_hidden":false}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Define the GPU device\ndevice_name = tf.test.gpu_device_name()","metadata":{"_uuid":"b75d52a0-de5c-4860-a977-c40f0fceee5a","_cell_guid":"a9aaf56c-c09f-4324-a785-970ad4b6db35","trusted":true,"collapsed":false,"execution":{"iopub.status.busy":"2025-03-20T08:24:29.556519Z","iopub.execute_input":"2025-03-20T08:24:29.556728Z","iopub.status.idle":"2025-03-20T08:24:29.562493Z","shell.execute_reply.started":"2025-03-20T08:24:29.556711Z","shell.execute_reply":"2025-03-20T08:24:29.561703Z"},"jupyter":{"outputs_hidden":false}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class RNAHybridModel:\n    \"\"\"Combines sequence embedding with engineered features\"\"\"\n    def __init__(self, max_seq_len, feature_dim):\n        self.max_seq_len = max_seq_len\n        self.feature_dim = feature_dim\n        self.model = self.build_model()\n\n    # Use the GPU device for model creation and training\n    def build_model(self):\n        # Use the GPU device for model creation\n        with tf.device(device_name):\n            # Enhanced sequence processing\n            seq_input = Input(shape=(self.max_seq_len,), name='seq_input')\n            x = Embedding(4, 256, mask_zero=False)(seq_input)\n            \n            # Stacked BiLSTMs with dropout\n            x = Bidirectional(LSTM(512, return_sequences=True, dropout=0.3))(x)\n            x = Bidirectional(LSTM(256, return_sequences=True, dropout=0.2))(x)\n            \n            # Transformer block with positional encoding\n            pos_enc = self.positional_encoding(self.max_seq_len, 512)\n            x = Add()([x, pos_enc])\n            x = MultiHeadAttention(num_heads=12, key_dim=64)(x, x)\n            x = GlobalAveragePooling1D()(x)\n            \n            # Feature fusion\n            feat_input = Input(shape=(self.feature_dim,), name='feat_input')\n            merged = concatenate([x, feat_input])\n            \n            # Enhanced dense processing\n            merged = Dense(1024, activation='swish', kernel_regularizer='l2')(merged)\n            merged = Dropout(0.5)(merged)\n            \n            # Create named output heads\n            x_head = self.create_coord_head(merged, 'x_out')\n            y_head = self.create_coord_head(merged, 'y_out')\n            z_head = self.create_coord_head(merged, 'z_out')\n            \n            return Model(inputs=[seq_input, feat_input], \n                       outputs=[x_head, y_head, z_head])\n\n    def create_coord_head(self, x, name):\n        x = Dense(768, activation='swish')(x)\n        x = Dense(384, activation='swish')(x)\n        return Dense(self.max_seq_len, activation='linear', name=name)(x)\n\n    def positional_encoding(self, length, depth):\n        positions = np.arange(length)[:, np.newaxis]\n        depths = np.arange(depth)[np.newaxis, :]/depth\n        angle_rates = 1 / (10000**depths)\n        angle_rads = positions * angle_rates\n        pos_encoding = np.concatenate(\n            [np.sin(angle_rads[:, ::2]), np.cos(angle_rads[:, ::2])],\n            axis=-1\n        )\n        return tf.cast(pos_encoding[np.newaxis, ...], tf.float32)\n\n    def compile_model(self):\n        self.model.compile(\n            optimizer=Adam(learning_rate=1e-4),\n            loss={'x_out': 'huber', 'y_out': 'huber', 'z_out': 'huber'},\n            loss_weights=[0.35, 0.35, 0.3],\n            metrics={'x_out': ['mae'], 'y_out': ['mae'], 'z_out': ['mae']} \n        )","metadata":{"_uuid":"316414ea-0807-441e-b793-805ee6a8bd2a","_cell_guid":"a3251d74-8159-439b-bc98-d0c9e4d8961e","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-03-20T08:24:29.571397Z","iopub.execute_input":"2025-03-20T08:24:29.571597Z","iopub.status.idle":"2025-03-20T08:24:29.586146Z","shell.execute_reply.started":"2025-03-20T08:24:29.571579Z","shell.execute_reply":"2025-03-20T08:24:29.585305Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class DataProcessor:\n    \"\"\"Handles temporal features and data validation\"\"\"\n    def __init__(self, max_seq_len=206):\n        self.max_seq_len = max_seq_len\n        self.feature_engineer = SequenceFeatureEngineer(max_seq_len)\n        \n    def add_temporal_features(self, df):\n        processed_df = df.copy()\n        # Ensure target_id is preserved\n        if 'target_id' not in processed_df.columns:\n            raise ValueError(\"target_id column missing in input data\")\n            \n        # Existing temporal feature calculations\n        processed_df['cutoff_year'] = pd.to_datetime(processed_df['temporal_cutoff']).dt.year\n        processed_df['seq_length'] = processed_df['sequence'].str.len()\n        \n        # Yearly averages\n        yearly_avg = processed_df.groupby('cutoff_year')['seq_length'].mean().to_dict()\n        processed_df['yearly_avg_length'] = processed_df['cutoff_year'].map(yearly_avg)\n        \n        return processed_df\n    \n    def preprocess_labels(self, label_df):\n        \"\"\"Process labels at the sequence level\"\"\"\n        coord_arrays = []\n        \n        for _, row in label_df.iterrows():\n            # Convert coordinate lists to numpy arrays\n            x = np.array(row['x_1'], dtype='float32')\n            y = np.array(row['y_1'], dtype='float32')\n            z = np.array(row['z_1'], dtype='float32')\n            \n            # Pad each coordinate sequence individually\n            x_pad = pad_sequences([x], maxlen=self.max_seq_len, padding='post')[0]\n            y_pad = pad_sequences([y], maxlen=self.max_seq_len, padding='post')[0]\n            z_pad = pad_sequences([z], maxlen=self.max_seq_len, padding='post')[0]\n            \n            # Combine into single array (x1, x2..., y1, y2..., z1, z2...)\n            coord_arrays.append(np.concatenate([x_pad, y_pad, z_pad]))\n            \n        return np.array(coord_arrays)","metadata":{"_uuid":"6d9a6120-7bd4-4674-8287-53eaea4f8175","_cell_guid":"6eafe82c-f4ac-4a18-bb1c-f952c7492a3f","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-03-20T08:24:29.587198Z","iopub.execute_input":"2025-03-20T08:24:29.587514Z","iopub.status.idle":"2025-03-20T08:24:29.607548Z","shell.execute_reply.started":"2025-03-20T08:24:29.587493Z","shell.execute_reply":"2025-03-20T08:24:29.606789Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def predict_structures(model, inputs, num_samples=5):\n    \"\"\"Returns predictions in (num_samples, batch_size, 3, max_seq_len) shape\"\"\"\n    mc_model = Model(inputs=model.inputs, outputs=model.outputs)\n    ensemble_preds = []\n    \n    for _ in range(num_samples):\n        preds = mc_model(inputs, training=True)\n        # Stack predictions as (batch_size, max_seq_len, 3)\n        stacked = np.stack(preds, axis=-1)\n        ensemble_preds.append(stacked)\n    \n    # Combine to (num_samples, batch_size, max_seq_len, 3)\n    combined = np.stack(ensemble_preds, axis=0)\n    # Reshape to (num_samples, batch_size, 3, max_seq_len)\n    return combined.transpose(0, 1, 3, 2)\n\ndef create_submission(model, test_df, sample_path, feature_engineer, max_seq_len):\n    \"\"\"Create submission using pre-fitted feature engineer\"\"\"\n    # Preserve original IDs and sequences\n    original_ids = test_df['target_id'].copy()\n    original_sequences = test_df['sequence'].copy()\n    \n    # Process test data with existing feature engineer\n    test_processed = feature_engineer.transform(test_df)\n    \n    # Prepare features\n    X_seq = np.array(test_processed['encoded'].tolist())\n    feature_columns = [col for col in test_processed.columns \n                      if col not in ['target_id', 'sequence', 'encoded', \n                                    'temporal_cutoff', 'description', 'all_sequences']]\n    X_feat = test_processed[feature_columns].values.astype('float32')\n    \n    # Generate predictions\n    ensemble_preds = predict_structures(model, [X_seq, X_feat])\n    \n    # Format predictions using original IDs and sequences\n    submission_rows = []\n    for idx in range(len(test_df)):\n        target_id = original_ids.iloc[idx]\n        sequence = original_sequences.iloc[idx]\n        seq_len = len(sequence)\n        \n        # Get predictions (num_samples, 3, max_seq_len)\n        mc_samples = ensemble_preds[:, idx, :, :]\n        \n        for pos in range(seq_len):\n            pos_idx = min(pos, max_seq_len-1)\n            entry = {\n                'ID': f\"{target_id}_{pos+1}\",\n                'resname': sequence[pos].upper(),  # Ensure uppercase\n                'resid': pos+1\n            }\n            \n            # Add all 5 samples\n            for sample_num in range(5):\n                entry.update({\n                    f'x_{sample_num+1}': mc_samples[sample_num, 0, pos_idx],\n                    f'y_{sample_num+1}': mc_samples[sample_num, 1, pos_idx],\n                    f'z_{sample_num+1}': mc_samples[sample_num, 2, pos_idx]\n                })\n            \n            submission_rows.append(entry)\n    \n    submission_df = pd.DataFrame(submission_rows)\n\n    print(submission_df.columns)\n    print(submission_df.shape)\n    print(submission_df.head())\n    \n    return submission_df\n\ndef compute_calibration_factors(true_x, true_y, true_z, pred_x, pred_y, pred_z):\n    \"\"\"\n    Compute calibration factors for each coordinate axis.\n    For each axis, the calibration factor is defined as:\n    \n        factor = mean(true_coordinate) / mean(predicted_coordinate)\n    \n    If the predicted mean is zero (to avoid division by zero), a factor of 1.0 is used.\n    \n    Parameters:\n        true_x, true_y, true_z: Arrays of true coordinates (flattened or over all residues)\n        pred_x, pred_y, pred_z: Arrays of predicted coordinates (flattened or over all residues)\n    \n    Returns:\n        A dictionary with keys 'x', 'y', 'z' and corresponding calibration factors.\n    \"\"\"\n    factors = {}\n    factors['x'] = np.mean(true_x) / np.mean(pred_x) if np.mean(pred_x) != 0 else 1.0\n    factors['y'] = np.mean(true_y) / np.mean(pred_y) if np.mean(pred_y) != 0 else 1.0\n    factors['z'] = np.mean(true_z) / np.mean(pred_z) if np.mean(pred_z) != 0 else 1.0\n    return factors\n\ndef apply_calibration(submission_df, calibration_factors):\n    \"\"\"\n    Adjust the submission DataFrame coordinates by applying calibration factors.\n    \n    For every coordinate column (x, y, and z) and for each sample (x_1, x_2, …, x_5),\n    multiply the predicted value by the corresponding calibration factor.\n    \n    Parameters:\n        submission_df (DataFrame): The submission DataFrame with columns like 'x_1', 'y_1', 'z_1', etc.\n        calibration_factors (dict): A dictionary containing calibration factors for 'x', 'y', and 'z'.\n    \n    Returns:\n        The adjusted submission DataFrame.\n    \"\"\"\n    for axis in ['x', 'y', 'z']:\n        for sample_num in range(1, 6):  # For samples 1 through 5\n            col = f'{axis}_{sample_num}'\n            submission_df[col] = submission_df[col] * calibration_factors.get(axis, 1.0)\n    return submission_df\n","metadata":{"_uuid":"9a73d0df-1366-4f1f-bf2e-73237418331e","_cell_guid":"078fc1dd-8306-4f85-9a14-64d3892ef357","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-03-20T08:24:29.609140Z","iopub.execute_input":"2025-03-20T08:24:29.609367Z","iopub.status.idle":"2025-03-20T08:24:29.630179Z","shell.execute_reply.started":"2025-03-20T08:24:29.609325Z","shell.execute_reply":"2025-03-20T08:24:29.629410Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Define a TensorBoard callback (log directory can be adjusted)\nlog_dir = \"./logs/fit/\" + datetime.now().strftime(\"%Y%m%d-%H%M%S\")\ntensorboard_cb = TensorBoard(log_dir=log_dir, histogram_freq=1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-20T08:24:29.631185Z","iopub.execute_input":"2025-03-20T08:24:29.631462Z","iopub.status.idle":"2025-03-20T08:24:29.648084Z","shell.execute_reply.started":"2025-03-20T08:24:29.631434Z","shell.execute_reply":"2025-03-20T08:24:29.647411Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"callbacks = [\n    EarlyStopping(patience=10, restore_best_weights=True),\n    ReduceLROnPlateau(factor=0.5, patience=3, min_lr=1e-6),\n    TerminateOnNaN(),\n    tensorboard_cb\n]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-20T08:24:29.648797Z","iopub.execute_input":"2025-03-20T08:24:29.648980Z","iopub.status.idle":"2025-03-20T08:24:29.669200Z","shell.execute_reply.started":"2025-03-20T08:24:29.648964Z","shell.execute_reply":"2025-03-20T08:24:29.668608Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def main():\n    # Configuration\n    # MAX_SEQ_LEN = 15\n    \n    # Load raw data with original columns\n    train_seq = pd.read_csv('/kaggle/input/stanford-rna-3d-folding/train_sequences.csv')\n    train_labels = pd.read_csv('/kaggle/input/stanford-rna-3d-folding/train_labels.csv')\n\n    # Handle Missing Values\n    train_labels['x_1'] = train_labels['x_1'].fillna(train_labels['x_1'].median())\n    train_labels['y_1'] = train_labels['y_1'].fillna(train_labels['y_1'].median())\n    train_labels['z_1'] = train_labels['z_1'].fillna(train_labels['z_1'].median())\n\n    # Compute lengths of sequences\n    seq_lengths = train_seq['sequence'].str.len().values\n\n    MAX_SEQ_LEN = int(np.percentile(seq_lengths, 95))  # Covers 95% of sequences\n    \n\n    # After loading data\n    # print(\"Original features:\", train_seq.columns.tolist())\n    \n    # 1. Add temporal features (requires original columns)\n    processor = DataProcessor(MAX_SEQ_LEN)\n    train_seq = processor.add_temporal_features(train_seq)\n    # print(\"After temporal processing:\", train_seq.columns.tolist())\n    \n    # 2. Add sequence features while preserving temporal features\n    feature_engineer = EnhancedFeatureEngineer(MAX_SEQ_LEN)\n    feature_engineer.fit(train_seq) \n    train_seq = feature_engineer.transform(train_seq)\n    # print(\"After feature engineering:\", train_seq.columns.tolist())\n\n    # Data augmentation (2x dataset size)\n    # print(\"Features going into augmentation:\", train_seq.columns.tolist())\n    # print(\"Applying data augmentation...\")\n    augmented_features, augmented_labels = augment_sequences(train_seq, train_labels)\n    # print(\"Feature columns after augmentation:\", train_seq.columns.tolist())\n    \n    # After augmentation\n    X_seq = np.stack(augmented_features['encoded'].values)\n    feature_columns = [col for col in augmented_features.columns \n                      if col not in ['encoded', 'target_id', 'ID', 'sequence_id', 'sequence']]  # Added 'sequence'\n    X_feat_array = augmented_features[feature_columns].values.astype('float32')\n    \n    # Process labels and split coordinates\n    padded_labels = processor.preprocess_labels(augmented_labels)\n    y = [\n        padded_labels[:, :MAX_SEQ_LEN],         # x coordinates\n        padded_labels[:, MAX_SEQ_LEN:2*MAX_SEQ_LEN],  # y coordinates\n        padded_labels[:, 2*MAX_SEQ_LEN:]        # z coordinates\n    ]\n    \n    # Add validation before training\n    # print(f\"Feature shapes: {X_seq.shape}, {X_feat_array.shape}\")\n    # print(f\"Label shapes: {[arr.shape for arr in y]}\")\n    \n    # Build and train model (assuming EnhancedRNAHybridModel is defined)\n    model_wrapper = RNAHybridModel(MAX_SEQ_LEN, X_feat_array.shape[1])\n    model_wrapper.compile_model()\n\n    # After model initialization\n    # print(model_wrapper.model.output_names)  # Should show ['x_out', 'y_out', 'z_out']\n    \n    history = model_wrapper.model.fit(\n        [X_seq, X_feat_array], y,\n        epochs=50,\n        batch_size=128,\n        validation_split=0.2,\n        callbacks=callbacks\n    )\n\n    # Predict on your validation set (or a subset of training data)\n    preds = model_wrapper.model.predict([X_seq, X_feat_array])\n    # preds is a list: [pred_x, pred_y, pred_z] each of shape (num_samples, MAX_SEQ_LEN)\n    \n    # Flatten predictions (ignoring padded zeros if needed)\n    pred_x = preds[0].flatten()\n    pred_y = preds[1].flatten()\n    pred_z = preds[2].flatten()\n    \n    # Similarly, flatten the true coordinates from padded_labels\n    true_x = padded_labels[:, :MAX_SEQ_LEN].flatten()\n    true_y = padded_labels[:, MAX_SEQ_LEN:2*MAX_SEQ_LEN].flatten()\n    true_z = padded_labels[:, 2*MAX_SEQ_LEN:].flatten()\n    \n    calibration_factors = compute_calibration_factors(true_x, true_y, true_z, pred_x, pred_y, pred_z)\n    print(\"Calibration factors:\", calibration_factors)\n\n    \n    # Create submission with the fitted feature engineer\n    test_df = pd.read_csv('/kaggle/input/stanford-rna-3d-folding/test_sequences.csv')\n    submission = create_submission(\n        model_wrapper.model,\n        test_df,\n        '/kaggle/input/stanford-rna-3d-folding/sample_submission.csv',\n        feature_engineer,\n        MAX_SEQ_LEN\n    )\n    # Post-prediction calibration if applicable\n    submission = apply_calibration(submission, calibration_factors)\n    submission.to_csv('submission.csv', index=False)\n\nif __name__ == \"__main__\":\n    main()","metadata":{"_uuid":"44907575-19d3-41fc-8325-305934a77ced","_cell_guid":"25471c88-30cf-4a90-b78f-78e1d92111f0","trusted":true,"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2025-03-20T08:24:29.670123Z","iopub.execute_input":"2025-03-20T08:24:29.670414Z","iopub.status.idle":"2025-03-20T08:32:58.011206Z","shell.execute_reply.started":"2025-03-20T08:24:29.670386Z","shell.execute_reply":"2025-03-20T08:32:58.010417Z"}},"outputs":[],"execution_count":null}]}