{"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.12.12"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":118765,"databundleVersionId":15231210},{"sourceType":"datasetVersion","sourceId":14604295,"datasetId":9328538,"databundleVersionId":15440074},{"sourceType":"datasetVersion","sourceId":14962460,"datasetId":9577079,"databundleVersionId":15833819},{"sourceType":"datasetVersion","sourceId":14874339,"datasetId":9502242,"databundleVersionId":15736806},{"sourceType":"datasetVersion","sourceId":14962495,"datasetId":9577097,"databundleVersionId":15833858},{"sourceType":"datasetVersion","sourceId":10855324,"datasetId":6742586,"databundleVersionId":11219268},{"sourceType":"kernelVersion","sourceId":301691387}],"dockerImageVersionId":31328,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":4505.337297,"end_time":"2026-03-18T20:11:07.518879","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2026-03-18T18:56:02.181582","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"e3929085-9567-4455-9094-64c92f5d030c","cell_type":"markdown","source":"Ok guys, this is not just a single “algorithm”; rather, this is an end-to-end RNA 3D structure prediction pipeline built for the Stanford RNA 3D Folding Kaggle competition. The document’s banner describes it as an “Enhanced Full ML Edition” with template-based modelling (TBM) + Protenix, feature engineering, a four-way split, anti-overfitting controls, ensembling, Optuna tuning, visualization, memory monitoring, and GPU-or-CPU transparent execution. In other words, it is a full research-grade workflow rather than a one-off script. \nI actually intentionally made the pipeline staging very deliberately. The log shows the order clearly: Loading data → EDA → Feature engineering → TBM coordinate cache → Split for ML meta-model → Model training + hyper-parameter tuning. \nThere is a very clear intention as well as indication here: that is designed as a controlled workflow: which involves first understanding the data, then converting sequences into numerical features, then building structure-related caches, then training and comparing models, and finally selecting the best model for output. \n\n# 1) The algorithmic logic\n   \nAt the base level, the code is engineered for robustness in a constrained environment. It first tries to install required packages from local wheels, which is why the script can run offline on Kaggle without internet access. It also detects whether PyTorch, BioPython, XGBoost, LightGBM, Optuna, and psutil are available, and falls back gracefully if some are missing. That makes the pipeline portable and fault-tolerant. \nIt then applies a set of resource-management controls: deterministic seeding, memory checks, dtype downcasting, garbage collection, and CUDA cache clearing when available. The purpose is to prevent the notebook from crashing on large RNA datasets or long runs. The code even enforces a total time budget and reserves time for final assembly/output, which is a common competition strategy when you need guaranteed completion under a deadline. \n\nThe data-loading layer is chunked and memory-safe. Instead of reading a large CSV all at once, the code reads in chunks, optimizes numeric dtypes, and stops early if memory usage approaches a threshold. That is particularly important here because the log shows the template set contains 5,028 sequences and 7,804,733 label rows, which is large enough to make naïve loading risky. \n\nFeature engineering\n\nOne of the main ideas in this code is to convert each RNA sequence into a compact set of hand-crafted descriptors. The code computes:\n•\tmono-nucleotide frequencies for A, C, G, U \n•\tdi-nucleotide frequencies \n•\tGC content and skew measures \n•\tpurine ratio \n•\tsequence length and log-length \n•\tsequence entropy and a simple complexity score \n•\tlongest homopolymer runs \n•\t5′ versus 3′ GC asymmetry \n•\ta simple base-pair potential heuristic \n•\tmotif counts like poly-A and poly-G runs. \n\nThe result is a 38-column feature table, which the log explicitly confirms. That means the model is not learning directly from raw strings; it is learning from a biologically informed feature representation meant to capture composition, asymmetry, and local motifs that may correlate with 3D fold behavior. \nData splitting\n\nThis pipeline uses a four-way split: Train 40% / Validation 15% / Test 15% / Holdout 30%. The holdout is separated first, and the code embeds an original index column to verify that no row leaks across splits. The split-verification function then checks every pair of partitions for overlap and aborts if leakage is detected. That is a strong experimental design choice because it makes the reported performance more trustworthy. \nThe log confirms the exact split sizes on this run: 2,012 train, 754 validation, 754 test, and 1,508 holdout, with no data leakage detected. \nOutlier handling\nBefore model fitting, the pipeline also flags outliers using an IQR-based rule. The log shows 2,003 out of 5,028 samples were detected as outliers. The code then supports down-weighting outlier rows rather than deleting them, which is a compromise between robustness and data retention. \nModeling strategy\nThe modeling stack is layered:\n                1.\tRandom Forest as a strong baseline. \n                2.\tXGBoost with Optuna tuning. \n                3.\tLightGBM with Optuna tuning. \n                4.\tA PyTorch MLP with batch normalization, dropout, GELU, attention, weight decay, gradient clipping, learning-rate scheduling, and early stopping.\n\n\n# 2.\tA stacking ensemble that combines the best tree-based models and uses a ridge meta-model. \nThis is important: the document is not betting on one architecture. It is trying several families of models and then choosing the one that generalizes best on the holdout set. That is classic ensemble-learning logic.\nThe MLP is especially noteworthy because the document explicitly adds anti-overfitting machinery: normalization, dropout, batch normalization, attention, early stopping, and validation-based checkpointing. The history object also tracks training and validation loss plus an R²-style score so the user can inspect learning dynamics. \nStructural reasoning / template logic\nBeyond the tabular ML layer, the code also contains a template-based modelling branch and a Protenix interface. It can parse stoichiometry and FASTA-like inputs, search for similar sequences using pairwise alignment, adapt template coordinates to a query sequence, apply chain-level geometric perturbations, stitch chunked coordinates with Kabsch alignment, and generate fallback idealized RNA-like coordinates when needed. That means the system is designed to produce or refine 3D coordinate predictions, not just scalar scores. \nThe document also defines a fast approximate TM-score for post-hoc comparison against reference coordinates. The intent is to have a geometry-aware quality metric for evaluating predicted structures, even if the full structural superposition machinery is not used every time. \n\n# 3) What the result log means\n\nThe result log is essentially a run audit trail. It tells you what succeeded, what data was seen, and how the models performed.\nThe first thing the log confirms is that the offline dependencies were installed successfully: Biopython, Biotite, and RDKit all installed from local wheels. It then detects PyTorch 2.10.0+cpu, so this particular run was CPU-only. \nIt also confirms the run context: 110 minutes total budget, 300 seconds per Protenix sample timeout, and 12,000 MB memory limit. That matters because the pipeline is explicitly budget-aware, so its design is shaped by compute constraints. \n\n    The data summary is:\n    •\t28 test targets \n    •\t5,028 template sequences \n    •\t7,804,733 label rows\n\nThen the log reports the generated exploratory plots were saved, including sequence statistics plots for templates and test sequences, plus label-distribution plots for train+val. It also reports the engineered feature count: 38. \nModel performance begins with Random Forest. The log shows RF already performs very strongly: R² = 0.9978 train, 0.9985 val, 0.9969 test, and then the run continues with other models and the ensemble. \nThe final summary is the most important part. The log states the pipeline completed in 130.8 minutes, and the best model was the stacking ensemble. Its reported scores were R² = 0.9994 train, 0.9990 val, 0.9982 test, and 0.9985 holdout. That is the system’s chosen final model according to holdout performance. \n\n# 3) What the graphical plots show\n\nThe plotting section is not decorative; it is diagnostic. Each plot answers a different question about the data or the model.\nThe EDA plots are meant to tell you whether the dataset is well-behaved and what kind of sequence variability exists. The code generates:\n        •\tLength distribution \n        •\tLog-length distribution \n        •\tGC content distribution \n        •\tNucleotide frequency box plots \n        •\tCorrelation matrix among basic features \n        •\tLength vs GC scatter plot \n\nThese plots help you see whether the sequences are roughly homogeneous or whether there are subpopulations that might require special handling.\n\nThe label EDA plots show the distribution of coordinate values, which is a quick check for the scale and spread of structural targets. That helps identify whether the coordinate labels are centered, skewed, noisy, or outlier-prone. \nThe split-summary plot shows the exact train/validation/test/holdout proportions visually, which is a sanity check that the four-way split was executed correctly. The log confirms this was saved. \n\nThe PCA and t-SNE plot is a visualization of the engineered feature space, not of the raw sequences. It is used to see whether the sequence features form a meaningful manifold and whether the target values vary smoothly across that space. The run log confirms this plot was saved for the training data. \n\n\nThe prediction scatter and residual plots are model-diagnostics. A prediction scatter checks whether predictions lie near the diagonal line of perfect agreement; residual plots show whether errors are centered around zero or whether the model is biased in some region. The code explicitly generates these for test and holdout evaluations.\nThe learning-curve plot is for the MLP specifically. It compares training and validation loss/accuracy over epochs and marks the best epoch, which helps diagnose underfitting, overfitting, or early convergence. \n\n\nThe generalization-gap plot compares performance across train, validation, test, and holdout. That is especially useful here because the best model is selected by holdout R², so the plot shows whether the model is truly stable or just memorizing the training set.\n\nThe code also defines a final holdout-comparison bar chart that compares all candidate models and prints the chosen best model. This is the final decision visual for model selection. \n\n# 4) Intended applications\n\nThe intended application is clear: RNA 3D structure prediction for the Stanford competition. More specifically, the pipeline is designed to map RNA sequence information to coordinate-level structural predictions and to rank candidate models by how well they generalize to unseen data.\n\nBeyond the competition, this kind of pipeline is useful for:\n\n    •\tbenchmarking sequence-to-structure modeling ideas, \n    •\ttesting whether hand-crafted RNA features carry predictive signal, \n    •\tcomparing classical ML versus deep learning versus ensembling, \n    •\tbuilding a reproducible structure-prediction workflow, \n    •\tand generating approximate 3D candidate structures for downstream analysis. \n    So the document’s purpose is broader than just “predict one output.” It is a full experimental framework for RNA structure modeling, with controlled data splitting, multiple learners, structural heuristics, and plotting for interpretability and model validation.\n    \nHope this Algorithmic Code helps everyone! Cheers and All the Best! \n","metadata":{}},{"id":"f913d77f-1bc9-421f-9b98-1989a66849a6","cell_type":"code","source":"\"\"\"\n=============================================================================\nRNA 3D Structure Prediction Pipeline — Enhanced Full ML Edition\n=============================================================================\nStanford RNA 3D Folding Competition — Enhanced Submission\nFeatures:\n  • Template-Based Modelling (TBM) + Protenix deep-learning phase\n  • Comprehensive EDA & feature engineering\n  • Four-way data split  (Train 40 / Val 15 / Test 15 / Holdout 30)\n  • Anti-overfitting: L2, Dropout, BatchNorm, EarlyStopping\n  • Ensemble / stacking: XGBoost + LightGBM + MLP\n  • Hyper-parameter tuning (Optuna)\n  • Full visualisation suite (all plt.savefig + plt.show)\n  • Real-time memory monitoring with psutil\n  • tqdm progress bars everywhere\n  • GPU-or-CPU transparent execution\n=============================================================================\n\"\"\"\n\n# ─────────────────────────── local wheel installs (offline Kaggle) ───────────\n# These must run before any import that depends on them.\nimport subprocess, sys\n\ndef _pip_install_wheel(wheel_path: str) -> None:\n    \"\"\"Install a local .whl file with no internet access required.\"\"\"\n    import os\n    if not os.path.exists(wheel_path):\n        print(f\"[wheel-install] Skipping (not found): {wheel_path}\")\n        return\n    result = subprocess.run(\n        [sys.executable, \"-m\", \"pip\", \"install\",\n         \"--no-index\", \"--no-deps\", wheel_path, \"--quiet\"],\n        capture_output=True, text=True,\n    )\n    if result.returncode == 0:\n        print(f\"[wheel-install] ✓ Installed: {wheel_path.split('/')[-1]}\")\n    else:\n        # Already installed or harmless warning — do not crash\n        print(f\"[wheel-install] Note ({wheel_path.split('/')[-1]}): {result.stderr.strip()[:200]}\")\n\n_BIOPYTHON_WHEEL = (\n    \"/kaggle/input/datasets/kami1976/biopython-cp312\"\n    \"/biopython-1.86-cp312-cp312-manylinux2014_x86_64\"\n    \".manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl\"\n)\n_BIOTITE_WHEEL = (\n    \"/kaggle/input/datasets/amirrezaaleyasin/biotite\"\n    \"/biotite-1.6.0-cp312-cp312-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl\"\n)\n_RDKIT_WHEEL = (\n    \"/kaggle/input/datasets/amirrezaaleyasin/rdkit-2025-9-5\"\n    \"/rdkit-2025.9.5-cp312-cp312-manylinux_2_28_x86_64.whl\"\n)\n\nfor _whl in (_BIOPYTHON_WHEEL, _BIOTITE_WHEEL, _RDKIT_WHEEL):\n    _pip_install_wheel(_whl)\n\n# ─────────────────────────── stdlib ──────────────────────────────────────────\nimport gc\nimport json\nimport logging\nimport os\nimport sys\nimport time\nimport traceback\nimport warnings\nfrom pathlib import Path\nfrom typing import Dict, List, Optional, Tuple\n\nwarnings.filterwarnings(\"ignore\")\nlogging.basicConfig(\n    level=logging.INFO,\n    format=\"%(asctime)s  %(levelname)-7s  %(message)s\",\n    datefmt=\"%H:%M:%S\",\n)\nlog = logging.getLogger(\"RNA3D\")\n\n# ─────────────────────────── third-party ─────────────────────────────────────\nimport numpy as np\nimport pandas as pd\nimport matplotlib\nmatplotlib.use(\"Agg\")                           # non-interactive backend\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nfrom matplotlib.patches import Patch\nimport seaborn as sns\nfrom tqdm import tqdm\n\ntry:\n    import psutil\n    HAS_PSUTIL = True\nexcept ImportError:\n    HAS_PSUTIL = False\n\ntry:\n    import torch\n    HAS_TORCH = True\n    DEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n    log.info(f\"PyTorch {torch.__version__} — device: {DEVICE}\")\nexcept ImportError:\n    HAS_TORCH = False\n    DEVICE = None\n    log.warning(\"PyTorch not found — deep-learning phase disabled.\")\n\ntry:\n    from Bio.Align import PairwiseAligner\n    HAS_BIOPYTHON = True\nexcept ImportError:\n    HAS_BIOPYTHON = False\n    log.warning(\"BioPython not found — TBM disabled.\")\n\ntry:\n    import xgboost as xgb\n    HAS_XGB = True\nexcept ImportError:\n    HAS_XGB = False\n\ntry:\n    import lightgbm as lgb\n    HAS_LGB = True\nexcept ImportError:\n    HAS_LGB = False\n\ntry:\n    import optuna\n    optuna.logging.set_verbosity(optuna.logging.WARNING)\n    HAS_OPTUNA = True\nexcept ImportError:\n    HAS_OPTUNA = False\n\nfrom sklearn.model_selection import StratifiedKFold, KFold, train_test_split\nfrom sklearn.preprocessing import StandardScaler, RobustScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.manifold import TSNE\nfrom sklearn.ensemble import (RandomForestRegressor, GradientBoostingRegressor,\n                               StackingRegressor)\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.impute import SimpleImputer\n\n# ─────────────────────────── paths & constants ───────────────────────────────\nos.environ.setdefault(\"LAYERNORM_TYPE\", \"torch\")\nos.environ.setdefault(\"RNA_MSA_DEPTH_LIMIT\", \"512\")\n\nIS_KAGGLE         = bool(os.environ.get(\"KAGGLE_IS_COMPETITION_RERUN\", \"\"))\nDATA_BASE         = os.environ.get(\"DATA_BASE\",\n                    \"/kaggle/input/competitions/stanford-rna-3d-folding-2\")\nDEFAULT_TEST_CSV  = f\"{DATA_BASE}/test_sequences.csv\"\nDEFAULT_TRAIN_CSV = f\"{DATA_BASE}/train_sequences.csv\"\nDEFAULT_TRAIN_LBL = f\"{DATA_BASE}/train_labels.csv\"\nDEFAULT_VAL_CSV   = f\"{DATA_BASE}/validation_sequences.csv\"\nDEFAULT_VAL_LBL   = f\"{DATA_BASE}/validation_labels.csv\"\nDEFAULT_OUTPUT    = \"/kaggle/working/submission.csv\"\nPLOT_DIR          = Path(\"/kaggle/working/plots\")\nMODEL_DIR         = Path(\"/kaggle/working/models\")\n\nDEFAULT_CODE_DIR  = os.environ.get(\n    \"PROTENIX_CODE_DIR\",\n    \"/kaggle/input/datasets/qiweiyin/protenix-v1-adjusted\"\n    \"/Protenix-v1-adjust-v2/Protenix-v1-adjust-v2/Protenix-v1\",\n)\nDEFAULT_ROOT_DIR  = os.environ.get(\"PROTENIX_ROOT_DIR\", DEFAULT_CODE_DIR)\n\nMODEL_NAME    = \"protenix_base_20250630_v1.0.0\"\nN_SAMPLE      = 5\nSEED          = 42\n# Tighter sequence-length cap: reduces per-sample Protenix cost on CPU\nMAX_SEQ_LEN   = int(os.environ.get(\"MAX_SEQ_LEN\",   \"256\"))\nCHUNK_OVERLAP = int(os.environ.get(\"CHUNK_OVERLAP\",  \"64\"))\nUSE_PROTENIX  = True\n\nMEMORY_LIMIT_MB   = int(os.environ.get(\"MEMORY_LIMIT_MB\", \"12000\"))\nCHUNK_ROWS        = int(os.environ.get(\"CHUNK_ROWS\",       \"50000\"))\nMAX_TEMPLATE_POOL = int(os.environ.get(\"MAX_TEMPLATE_POOL\",\"5000\"))\n\n# ── Timing budget ─────────────────────────────────────────────────────────────\n# Total wall-clock budget (default 110 min, leaving 10 min safety margin).\nTOTAL_BUDGET_S            = int(os.environ.get(\"TOTAL_BUDGET_S\",            str(110 * 60)))\n# Max seconds for one Protenix sample (featurise + predict).\nPROTENIX_SAMPLE_TIMEOUT_S = int(os.environ.get(\"PROTENIX_SAMPLE_TIMEOUT_S\", \"300\"))\n# Min seconds that must remain before starting the next Protenix sample.\nPROTENIX_MIN_REMAINING_S  = int(os.environ.get(\"PROTENIX_MIN_REMAINING_S\",  \"600\"))\n# Seconds reserved for final assembly + CSV write.\nASSEMBLY_RESERVE_S        = int(os.environ.get(\"ASSEMBLY_RESERVE_S\",        \"300\"))\n# Set at the top of main() — do not mutate elsewhere.\n_PIPELINE_START: float = 0.0\n\n\ndef _elapsed() -> float:\n    \"\"\"Seconds elapsed since pipeline start.\"\"\"\n    return time.time() - _PIPELINE_START\n\n\ndef _remaining() -> float:\n    \"\"\"Seconds remaining in the total budget.\"\"\"\n    return TOTAL_BUDGET_S - _elapsed()\n\n\ndef _time_ok_for_protenix() -> bool:\n    \"\"\"True if there is enough budget to attempt another Protenix sample.\"\"\"\n    return _remaining() > (PROTENIX_MIN_REMAINING_S + ASSEMBLY_RESERVE_S)\n\n# ─────────────────────────── helpers ─────────────────────────────────────────\n\ndef parse_bool(v: str, default: bool = False) -> str:\n    s = str(v).strip().lower()\n    if s in {\"1\",\"true\",\"t\",\"yes\",\"y\",\"on\"}:  return \"true\"\n    if s in {\"0\",\"false\",\"f\",\"no\",\"n\",\"off\"}: return \"false\"\n    return \"true\" if default else \"false\"\n\nUSE_MSA      = parse_bool(os.environ.get(\"USE_MSA\",     \"false\"))\nUSE_TEMPLATE = parse_bool(os.environ.get(\"USE_TEMPLATE\",\"false\"))\nUSE_RNA_MSA  = parse_bool(os.environ.get(\"USE_RNA_MSA\", \"true\"))\nMODEL_N_SAMPLE = int(os.environ.get(\"MODEL_N_SAMPLE\",  str(N_SAMPLE)))\nMIN_SIMILARITY       = float(os.environ.get(\"MIN_SIMILARITY\",       \"0.0\"))\nMIN_PERCENT_IDENTITY = float(os.environ.get(\"MIN_PERCENT_IDENTITY\", \"50.0\"))\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 0: MEMORY MANAGEMENT\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef mem_mb() -> float:\n    if HAS_PSUTIL:\n        return psutil.Process().memory_info().rss / 1e6\n    return 0.0\n\n\ndef force_cleanup(*args) -> None:\n    \"\"\"Delete named variables and run GC + CUDA cache clear.\"\"\"\n    for obj in args:\n        try:\n            del obj\n        except Exception:\n            pass\n    gc.collect()\n    if HAS_TORCH and torch.cuda.is_available():\n        torch.cuda.empty_cache()\n    gc.collect()\n\n\ndef memory_safe(label: str = \"\") -> bool:\n    \"\"\"Return False (and warn) if we are near the memory limit.\"\"\"\n    used = mem_mb()\n    if used > MEMORY_LIMIT_MB * 0.90:\n        log.warning(f\"[MEM] {label}: {used:.0f}/{MEMORY_LIMIT_MB} MB — near limit!\")\n        force_cleanup()\n        return False\n    return True\n\n\ndef optimise_dtypes(df: pd.DataFrame) -> pd.DataFrame:\n    \"\"\"Downcast numeric columns to save RAM.\"\"\"\n    for col in df.select_dtypes(\"float64\").columns:\n        df[col] = df[col].astype(\"float32\")\n    for col in df.select_dtypes(\"int64\").columns:\n        df[col] = pd.to_numeric(df[col], downcast=\"integer\")\n    return df\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 1: SEED & DEVICE\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef seed_everything(seed: int = SEED) -> None:\n    os.environ[\"CUBLAS_WORKSPACE_CONFIG\"] = \":4096:8\"\n    np.random.seed(seed)\n    if HAS_TORCH:\n        torch.manual_seed(seed)\n        torch.cuda.manual_seed_all(seed)\n        torch.backends.cudnn.benchmark     = False\n        torch.backends.cudnn.deterministic = True\n        try:\n            torch.use_deterministic_algorithms(True)\n        except Exception:\n            pass\n\n\nseed_everything(SEED)\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 2: DATA LOADING (chunked, memory-safe)\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef load_csv_chunked(path: str, max_rows: Optional[int] = None) -> pd.DataFrame:\n    \"\"\"Load a potentially large CSV in chunks, respecting memory limits.\"\"\"\n    if not os.path.exists(path):\n        log.warning(f\"File not found: {path}\")\n        return pd.DataFrame()\n    chunks, total = [], 0\n    desc = f\"Loading {Path(path).name}\"\n    for chunk in tqdm(\n        pd.read_csv(path, chunksize=CHUNK_ROWS),\n        desc=desc, leave=False,\n    ):\n        chunk = optimise_dtypes(chunk)\n        chunks.append(chunk)\n        total += len(chunk)\n        if max_rows and total >= max_rows:\n            break\n        if not memory_safe(desc):\n            log.warning(\"  stopping early — memory limit.\")\n            break\n    if not chunks:\n        return pd.DataFrame()\n    df = pd.concat(chunks, ignore_index=True)\n    if max_rows:\n        df = df.head(max_rows)\n    force_cleanup(chunks)\n    return df\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 3: FEATURE ENGINEERING\n# ══════════════════════════════════════════════════════════════════════════════\n\nNUCLEOTIDES = list(\"ACGU\")\nDINUCLEOTIDES = [a+b for a in NUCLEOTIDES for b in NUCLEOTIDES]\n\n\ndef compute_sequence_features(seq: str) -> Dict:\n    \"\"\"Rich hand-crafted features from an RNA sequence string.\"\"\"\n    n = len(seq) + 1e-9\n    feats: Dict = {}\n\n    # 1. Mono-nucleotide frequencies\n    for nt in NUCLEOTIDES:\n        feats[f\"freq_{nt}\"] = seq.count(nt) / n\n\n    # 2. Di-nucleotide frequencies\n    for dn in DINUCLEOTIDES:\n        cnt = sum(1 for i in range(len(seq)-1) if seq[i:i+2] == dn)\n        feats[f\"di_{dn}\"] = cnt / max(len(seq)-1, 1)\n\n    # 3. GC content and derived\n    gc = (seq.count(\"G\") + seq.count(\"C\")) / n\n    feats[\"gc_content\"]   = gc\n    feats[\"gc_skew\"]      = (seq.count(\"G\") - seq.count(\"C\")) / (seq.count(\"G\") + seq.count(\"C\") + 1e-9)\n    feats[\"au_skew\"]      = (seq.count(\"A\") - seq.count(\"U\")) / (seq.count(\"A\") + seq.count(\"U\") + 1e-9)\n    feats[\"purine_ratio\"] = (seq.count(\"A\") + seq.count(\"G\")) / n\n\n    # 4. Length & complexity\n    feats[\"length\"]       = len(seq)\n    feats[\"log_length\"]   = np.log1p(len(seq))\n    unique_chars = len(set(seq))\n    feats[\"seq_entropy\"]  = -sum(\n        (seq.count(c)/n) * np.log2(seq.count(c)/n + 1e-9)\n        for c in set(seq)\n    )\n    feats[\"complexity\"]   = unique_chars / 4.0\n\n    # 5. Run-length features (longest homopolymer)\n    runs = {\"A\": 0, \"C\": 0, \"G\": 0, \"U\": 0}\n    cur_char, cur_run = \"\", 0\n    for ch in seq:\n        if ch == cur_char:\n            cur_run += 1\n        else:\n            if cur_char in runs:\n                runs[cur_char] = max(runs[cur_char], cur_run)\n            cur_char, cur_run = ch, 1\n    for nt in NUCLEOTIDES:\n        feats[f\"max_run_{nt}\"] = runs[nt] / n\n\n    # 6. Positional bias (5' vs 3' GC)\n    half = len(seq) // 2\n    feats[\"gc_5prime\"] = ((seq[:half].count(\"G\") + seq[:half].count(\"C\")) / (half + 1e-9))\n    feats[\"gc_3prime\"] = ((seq[half:].count(\"G\") + seq[half:].count(\"C\")) / (len(seq)-half + 1e-9))\n    feats[\"gc_asymmetry\"] = feats[\"gc_5prime\"] - feats[\"gc_3prime\"]\n\n    # 7. Predicted base-pair potential (simple heuristic: CG/GC/AU/UA pairs)\n    bp_score = sum(\n        1 for i in range(len(seq)//2)\n        if (seq[i] + seq[len(seq)-1-i] in {\"CG\",\"GC\",\"AU\",\"UA\"})\n    )\n    feats[\"bp_potential\"] = bp_score / (len(seq)//2 + 1e-9)\n\n    # 8. Secondary structure motif counts (simple pattern scan)\n    feats[\"poly_A_count\"] = sum(1 for i in range(len(seq)-3) if seq[i:i+4] == \"AAAA\")\n    feats[\"poly_G_count\"] = sum(1 for i in range(len(seq)-3) if seq[i:i+4] == \"GGGG\")\n\n    return feats\n\n\ndef build_feature_dataframe(seqs_df: pd.DataFrame,\n                             pbar_desc: str = \"Feature engineering\") -> pd.DataFrame:\n    records = []\n    for _, row in tqdm(seqs_df.iterrows(), total=len(seqs_df), desc=pbar_desc, leave=False):\n        f = compute_sequence_features(str(row[\"sequence\"]))\n        f[\"target_id\"] = row[\"target_id\"]\n        records.append(f)\n    df = pd.DataFrame(records)\n    df = optimise_dtypes(df)\n    return df\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 4: DATA SPLITS (four-way, no leakage)\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef four_way_split(df: pd.DataFrame,\n                   seed: int = SEED\n                   ) -> Tuple[pd.DataFrame, pd.DataFrame,\n                               pd.DataFrame, pd.DataFrame]:\n    \"\"\"\n    Train 40% | Val 15% | Test 15% | Holdout 30%\n    Holdout is set aside FIRST so it never touches training or validation.\n\n    A sentinel column ``__orig_idx`` carrying the original integer row position\n    is embedded in every returned split so that ``verify_data_splits`` can\n    confirm zero overlap even when the dataframe has no natural key column.\n    \"\"\"\n    n = len(df)\n    log.info(f\"Total samples for splitting: {n}\")\n    rng   = np.random.default_rng(seed)\n    idx   = rng.permutation(n)           # shuffled positions 0..n-1\n\n    h_end  = int(n * 0.30)\n    t_end  = int(n * 0.45)              # holdout 30% + test 15%\n    v_end  = int(n * 0.60)              # + val 15%\n    # remainder → train 40%\n\n    holdout_pos = idx[:h_end]\n    test_pos    = idx[h_end:t_end]\n    val_pos     = idx[t_end:v_end]\n    train_pos   = idx[v_end:]\n\n    def _make_split(positions):\n        part = df.iloc[positions].copy().reset_index(drop=True)\n        # Embed original positions as a unique key for overlap checking\n        part[\"__orig_idx\"] = positions\n        return part\n\n    splits = {\n        \"train\":   _make_split(train_pos),\n        \"val\":     _make_split(val_pos),\n        \"test\":    _make_split(test_pos),\n        \"holdout\": _make_split(holdout_pos),\n    }\n    for k, v in splits.items():\n        log.info(f\"  {k:8s}: {len(v):6d} rows  ({100*len(v)/n:.1f}%)\")\n    return splits[\"train\"], splits[\"val\"], splits[\"test\"], splits[\"holdout\"]\n\n\ndef verify_data_splits(train, val, test, holdout) -> None:\n    \"\"\"\n    Assert zero overlap between splits.\n\n    Priority order for unique-ID selection:\n      1. ``__orig_idx``  — injected by four_way_split (always unique, no reset)\n      2. ``target_id``   — natural competition key when present\n      3. Raise an informative error rather than silently passing wrong IDs\n    \"\"\"\n    split_names = [\"train\", \"val\", \"test\", \"holdout\"]\n    dfs         = [train, val, test, holdout]\n\n    def _ids(df: pd.DataFrame) -> set:\n        if \"__orig_idx\" in df.columns:\n            return set(df[\"__orig_idx\"].tolist())\n        if \"target_id\" in df.columns:\n            return set(df[\"target_id\"].tolist())\n        raise ValueError(\n            \"verify_data_splits: dataframe has neither '__orig_idx' nor \"\n            \"'target_id' — cannot verify split integrity.\"\n        )\n\n    id_sets = [_ids(df) for df in dfs]\n\n    pairs = [\n        (\"train\", \"val\"),\n        (\"train\", \"test\"),\n        (\"train\", \"holdout\"),\n        (\"val\",   \"test\"),\n        (\"val\",   \"holdout\"),\n        (\"test\",  \"holdout\"),\n    ]\n    for a, b in pairs:\n        ia, ib   = id_sets[split_names.index(a)], id_sets[split_names.index(b)]\n        overlap  = ia & ib\n        assert len(overlap) == 0, (\n            f\"Data leakage between '{a}' and '{b}': \"\n            f\"{len(overlap)} shared row(s) detected.\"\n        )\n    log.info(\"✓ No data leakage detected across all splits.\")\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 5: EDA & VISUALISATION\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef _save_show(fig: plt.Figure, name: str, dpi: int = 150) -> None:\n    PLOT_DIR.mkdir(parents=True, exist_ok=True)\n    path = PLOT_DIR / f\"{name}.png\"\n    fig.savefig(path, dpi=dpi, bbox_inches=\"tight\")\n    log.info(f\"  → saved {path}\")\n    plt.show()\n    plt.close(fig)\n\n\ndef eda_sequence_stats(df: pd.DataFrame, label: str = \"all\") -> None:\n    \"\"\"EDA: length distribution, GC content, nucleotide composition.\"\"\"\n    if df.empty or \"sequence\" not in df.columns:\n        return\n    lengths = df[\"sequence\"].str.len()\n    gc = df[\"sequence\"].apply(\n        lambda s: (s.count(\"G\") + s.count(\"C\")) / max(len(s), 1)\n    )\n    nt_freq = {nt: df[\"sequence\"].str.count(nt).div(lengths) for nt in NUCLEOTIDES}\n\n    fig, axes = plt.subplots(2, 3, figsize=(16, 9))\n    fig.suptitle(f\"EDA — Sequence Statistics [{label}]\", fontsize=14, fontweight=\"bold\")\n\n    # 1. Length distribution\n    ax = axes[0, 0]\n    ax.hist(lengths, bins=40, color=\"#2196F3\", edgecolor=\"white\", alpha=0.85)\n    ax.axvline(lengths.median(), color=\"red\", ls=\"--\", label=f\"Median={lengths.median():.0f}\")\n    ax.set_xlabel(\"Sequence length (nt)\"); ax.set_ylabel(\"Count\")\n    ax.set_title(\"Length Distribution\"); ax.legend()\n\n    # 2. Log-scale length\n    ax = axes[0, 1]\n    ax.hist(np.log1p(lengths), bins=40, color=\"#4CAF50\", edgecolor=\"white\", alpha=0.85)\n    ax.set_xlabel(\"log(1+length)\"); ax.set_title(\"Log-Length Distribution\")\n\n    # 3. GC content\n    ax = axes[0, 2]\n    ax.hist(gc, bins=40, color=\"#FF9800\", edgecolor=\"white\", alpha=0.85)\n    ax.axvline(gc.mean(), color=\"navy\", ls=\"--\", label=f\"Mean={gc.mean():.2f}\")\n    ax.set_xlabel(\"GC content\"); ax.set_title(\"GC Content Distribution\"); ax.legend()\n\n    # 4. Nucleotide frequencies (box plot)\n    ax = axes[1, 0]\n    data_box = pd.DataFrame({nt: nt_freq[nt] for nt in NUCLEOTIDES})\n    data_box.boxplot(ax=ax, grid=False,\n                     boxprops=dict(color=\"#333\"), medianprops=dict(color=\"red\"))\n    ax.set_title(\"Nucleotide Frequency Distributions\")\n    ax.set_ylabel(\"Relative frequency\")\n\n    # 5. Correlation of features\n    ax = axes[1, 1]\n    feat_df = pd.DataFrame({\"length\": lengths, \"gc\": gc,\n                             **{nt: nt_freq[nt] for nt in NUCLEOTIDES}})\n    corr = feat_df.corr()\n    im = ax.imshow(corr.values, cmap=\"coolwarm\", vmin=-1, vmax=1)\n    ax.set_xticks(range(len(corr.columns)))\n    ax.set_yticks(range(len(corr.columns)))\n    ax.set_xticklabels(corr.columns, rotation=45, ha=\"right\", fontsize=8)\n    ax.set_yticklabels(corr.columns, fontsize=8)\n    fig.colorbar(im, ax=ax, fraction=0.046)\n    ax.set_title(\"Feature Correlation Matrix\")\n\n    # 6. Length vs GC scatter\n    ax = axes[1, 2]\n    ax.scatter(lengths, gc, alpha=0.3, s=8, c=gc, cmap=\"viridis\")\n    ax.set_xlabel(\"Length\"); ax.set_ylabel(\"GC content\")\n    ax.set_title(\"Length vs GC Content\")\n\n    _save_show(fig, f\"eda_sequence_{label}\")\n    force_cleanup(data_box, feat_df, corr)\n\n\ndef eda_label_stats(labels_df: pd.DataFrame, label: str = \"train\") -> None:\n    \"\"\"EDA: C1′ coordinate distributions from ground-truth labels.\"\"\"\n    if labels_df.empty:\n        return\n    coord_cols = [c for c in labels_df.columns if c.startswith((\"x_\",\"y_\",\"z_\"))]\n    if not coord_cols:\n        return\n    sample_cols = coord_cols[:3]   # x_1, y_1, z_1\n\n    fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n    fig.suptitle(f\"C1′ Coordinate Distributions [{label}]\", fontsize=13, fontweight=\"bold\")\n    colors = [\"#e74c3c\", \"#2ecc71\", \"#3498db\"]\n    for ax, col, c in zip(axes, sample_cols, colors):\n        vals = labels_df[col].dropna()\n        ax.hist(vals, bins=60, color=c, alpha=0.8, edgecolor=\"white\")\n        ax.set_title(f\"{col}  μ={vals.mean():.1f}  σ={vals.std():.1f}\")\n        ax.set_xlabel(\"Å\"); ax.set_ylabel(\"Count\")\n    _save_show(fig, f\"eda_labels_{label}\")\n\n\ndef plot_split_summary(train, val, test, holdout) -> None:\n    \"\"\"Bar chart of split sizes.\"\"\"\n    fig, ax = plt.subplots(figsize=(7, 4))\n    splits  = [\"Train\\n(40%)\", \"Validation\\n(15%)\", \"Test\\n(15%)\", \"Holdout\\n(30%)\"]\n    counts  = [len(train), len(val), len(test), len(holdout)]\n    colors  = [\"#3498db\", \"#2ecc71\", \"#e74c3c\", \"#9b59b6\"]\n    bars = ax.bar(splits, counts, color=colors, edgecolor=\"white\", linewidth=1.2)\n    for bar, cnt in zip(bars, counts):\n        ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5,\n                str(cnt), ha=\"center\", va=\"bottom\", fontsize=10)\n    ax.set_ylabel(\"Number of samples\")\n    ax.set_title(\"Four-Way Data Split\")\n    _save_show(fig, \"split_summary\")\n\n\ndef plot_feature_importance(importances: np.ndarray,\n                             feature_names: List[str],\n                             model_name: str = \"Model\") -> None:\n    top_k = min(25, len(feature_names))\n    idx   = np.argsort(importances)[-top_k:]\n    fig, ax = plt.subplots(figsize=(9, max(4, top_k*0.32)))\n    ax.barh(range(top_k), importances[idx], color=\"#3498db\", alpha=0.85)\n    ax.set_yticks(range(top_k))\n    ax.set_yticklabels([feature_names[i] for i in idx], fontsize=8)\n    ax.set_xlabel(\"Importance\")\n    ax.set_title(f\"Feature Importance — {model_name}\")\n    _save_show(fig, f\"feat_importance_{model_name.replace(' ','_')}\")\n\n\ndef plot_pca_tsne(X: np.ndarray, y: np.ndarray, label: str = \"train\") -> None:\n    \"\"\"PCA and t-SNE projections of the feature space.\"\"\"\n    # Ensure y is always a 1-D array regardless of how it was passed in\n    y = np.asarray(y).ravel()\n\n    n = min(len(X), 2000)\n    idx = np.random.default_rng(SEED).choice(len(X), n, replace=False)\n    Xs, ys = X[idx], y[idx]\n\n    scaler = StandardScaler()\n    Xs_sc  = scaler.fit_transform(Xs)\n\n    # PCA\n    pca    = PCA(n_components=2, random_state=SEED)\n    Xpca   = pca.fit_transform(Xs_sc)\n\n    fig, axes = plt.subplots(1, 2, figsize=(13, 5))\n    fig.suptitle(f\"PCA & t-SNE Feature Space [{label}]\", fontsize=13, fontweight=\"bold\")\n    sc = axes[0].scatter(Xpca[:,0], Xpca[:,1], c=ys, cmap=\"plasma\",\n                         alpha=0.6, s=10)\n    plt.colorbar(sc, ax=axes[0])\n    axes[0].set_title(f\"PCA  (var={pca.explained_variance_ratio_.sum()*100:.1f}%)\")\n    axes[0].set_xlabel(\"PC1\"); axes[0].set_ylabel(\"PC2\")\n\n    # t-SNE (subsample more if large)\n    n2   = min(n, 800)\n    idx2 = np.random.default_rng(SEED+1).choice(n, n2, replace=False)\n    tsne = TSNE(n_components=2, random_state=SEED, n_iter=300, perplexity=30)\n    Xt   = tsne.fit_transform(Xs_sc[idx2])\n    sc2  = axes[1].scatter(Xt[:,0], Xt[:,1], c=ys[idx2], cmap=\"plasma\",\n                            alpha=0.6, s=10)\n    plt.colorbar(sc2, ax=axes[1])\n    axes[1].set_title(\"t-SNE\"); axes[1].set_xlabel(\"Dim 1\"); axes[1].set_ylabel(\"Dim 2\")\n    _save_show(fig, f\"pca_tsne_{label}\")\n    force_cleanup(Xs, ys, Xs_sc, Xpca, Xt)\n\n\ndef plot_learning_curves(history: Dict, title: str = \"Learning Curves\") -> None:\n    \"\"\"Plot train/val loss and accuracy curves with peak marked.\"\"\"\n    if not history:\n        return\n    fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n    fig.suptitle(title, fontsize=13, fontweight=\"bold\")\n\n    for ax, metric in zip(axes, [\"loss\", \"accuracy\"]):\n        if metric not in history:\n            continue\n        tr = history[metric]\n        vl = history.get(f\"val_{metric}\", [])\n        ax.plot(tr, label=\"Train\", color=\"#2196F3\")\n        if vl:\n            ax.plot(vl, label=\"Validation\", color=\"#FF9800\")\n            best_ep = int(np.argmin(vl) if metric == \"loss\" else np.argmax(vl))\n            ax.axvline(best_ep, color=\"red\", ls=\"--\", alpha=0.6,\n                       label=f\"Best epoch={best_ep}\")\n        ax.set_xlabel(\"Epoch\"); ax.set_ylabel(metric.capitalize())\n        ax.set_title(metric.capitalize()); ax.legend()\n    _save_show(fig, f\"learning_curves_{title.replace(' ','_')}\")\n\n\ndef plot_generalization_gap(scores: Dict[str, float], title: str = \"Generalization\") -> None:\n    \"\"\"Bar chart comparing train/val/test/holdout performance.\"\"\"\n    fig, ax = plt.subplots(figsize=(8, 4))\n    colors = {\"train\":\"#3498db\",\"val\":\"#2ecc71\",\"test\":\"#e74c3c\",\"holdout\":\"#9b59b6\"}\n    for split, val_s in scores.items():\n        ax.bar(split, val_s, color=colors.get(split,\"#888\"), alpha=0.85, edgecolor=\"white\")\n        ax.text(split, val_s + 0.001, f\"{val_s:.4f}\", ha=\"center\", va=\"bottom\", fontsize=9)\n    ax.set_ylabel(\"Score (higher=better)\")\n    ax.set_title(title)\n    if \"train\" in scores and \"holdout\" in scores:\n        gap = abs(scores[\"train\"] - scores[\"holdout\"])\n        ax.set_xlabel(f\"Generalisation gap (train–holdout) = {gap:.4f}\")\n    _save_show(fig, f\"generalization_{title.replace(' ','_')}\")\n\n\ndef plot_prediction_scatter(y_true: np.ndarray,\n                             y_pred: np.ndarray,\n                             split_name: str = \"test\") -> None:\n    y_true = np.asarray(y_true).ravel()\n    y_pred = np.asarray(y_pred).ravel()\n    fig, ax = plt.subplots(figsize=(6, 6))\n    ax.scatter(y_true, y_pred, alpha=0.3, s=8, c=\"#3498db\")\n    lims = [min(y_true.min(), y_pred.min()), max(y_true.max(), y_pred.max())]\n    ax.plot(lims, lims, \"r--\", lw=1.5, label=\"Perfect prediction\")\n    r2 = r2_score(y_true, y_pred)\n    ax.set_title(f\"True vs Predicted [{split_name}]  R²={r2:.4f}\")\n    ax.set_xlabel(\"True\"); ax.set_ylabel(\"Predicted\"); ax.legend()\n    _save_show(fig, f\"pred_scatter_{split_name}\")\n\n\ndef plot_residuals(y_true: np.ndarray,\n                   y_pred: np.ndarray,\n                   split_name: str = \"test\") -> None:\n    y_true = np.asarray(y_true).ravel()\n    y_pred = np.asarray(y_pred).ravel()\n    resid = y_true - y_pred\n    fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n    axes[0].scatter(y_pred, resid, alpha=0.3, s=8, c=\"#e74c3c\")\n    axes[0].axhline(0, color=\"black\", lw=1)\n    axes[0].set_xlabel(\"Predicted\"); axes[0].set_ylabel(\"Residual\")\n    axes[0].set_title(f\"Residuals vs Fitted [{split_name}]\")\n    axes[1].hist(resid, bins=50, color=\"#9b59b6\", edgecolor=\"white\")\n    axes[1].set_xlabel(\"Residual\"); axes[1].set_title(\"Residual Distribution\")\n    _save_show(fig, f\"residuals_{split_name}\")\n\n\ndef plot_coordinate_predictions(pred_coords: np.ndarray,\n                                  target_id: str,\n                                  n_samples: int = N_SAMPLE) -> None:\n    \"\"\"3-D scatter of all N_SAMPLE predicted structures for one target.\"\"\"\n    if pred_coords.ndim != 3 or pred_coords.shape[0] < 1:\n        return\n    colors = plt.cm.tab10(np.linspace(0, 1, n_samples))\n    fig    = plt.figure(figsize=(9, 7))\n    ax     = fig.add_subplot(111, projection=\"3d\")\n    for s in range(min(n_samples, pred_coords.shape[0])):\n        xyz = pred_coords[s]\n        ax.plot(xyz[:,0], xyz[:,1], xyz[:,2], alpha=0.6, lw=0.8,\n                color=colors[s], label=f\"Sample {s+1}\")\n    ax.set_title(f\"3-D Predicted Structures — {target_id}\")\n    ax.set_xlabel(\"x (Å)\"); ax.set_ylabel(\"y (Å)\"); ax.set_zlabel(\"z (Å)\")\n    ax.legend(fontsize=7)\n    _save_show(fig, f\"3d_structure_{target_id[:20]}\")\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 6: OUTLIER DETECTION & WEIGHT ADJUSTMENT\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef detect_outliers_iqr(arr: np.ndarray,\n                         k: float = 2.5) -> np.ndarray:\n    \"\"\"Return boolean mask True for outlier rows.\"\"\"\n    q1, q3 = np.percentile(arr, 25, axis=0), np.percentile(arr, 75, axis=0)\n    iqr = q3 - q1 + 1e-9\n    lo, hi = q1 - k * iqr, q3 + k * iqr\n    return np.any((arr < lo) | (arr > hi), axis=1)\n\n\ndef sample_weights_from_outliers(X: np.ndarray,\n                                  outlier_mask: np.ndarray,\n                                  outlier_weight: float = 0.25) -> np.ndarray:\n    w = np.ones(len(X), dtype=\"float32\")\n    w[outlier_mask] = outlier_weight\n    return w\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 7: ML MODELS\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef build_rf(params: Optional[Dict] = None) -> RandomForestRegressor:\n    p = {\"n_estimators\": 120, \"max_depth\": 12, \"min_samples_leaf\": 4,\n         \"n_jobs\": -1, \"random_state\": SEED}\n    if params:\n        p.update(params)\n    return RandomForestRegressor(**p)\n\n\ndef build_xgb(params: Optional[Dict] = None) -> \"xgb.XGBRegressor\":\n    p = {\"n_estimators\": 200, \"max_depth\": 6, \"learning_rate\": 0.05,\n         \"subsample\": 0.8, \"colsample_bytree\": 0.8, \"reg_alpha\": 0.1,\n         \"reg_lambda\": 1.0, \"tree_method\": \"hist\",\n         \"device\": \"cuda\" if (HAS_TORCH and torch.cuda.is_available()) else \"cpu\",\n         \"random_state\": SEED, \"verbosity\": 0}\n    if params:\n        p.update(params)\n    return xgb.XGBRegressor(**p) if HAS_XGB else None\n\n\ndef build_lgb(params: Optional[Dict] = None) -> \"lgb.LGBMRegressor\":\n    p = {\"n_estimators\": 200, \"num_leaves\": 63, \"learning_rate\": 0.05,\n         \"subsample\": 0.8, \"colsample_bytree\": 0.8,\n         \"reg_alpha\": 0.1, \"reg_lambda\": 1.0,\n         \"n_jobs\": -1, \"random_state\": SEED, \"verbosity\": -1}\n    if params:\n        p.update(params)\n    return lgb.LGBMRegressor(**p) if HAS_LGB else None\n\n\n# ─── Optuna hyper-parameter tuning ──────────────────────────────────────────\n\ndef tune_xgb(X_tr, y_tr, X_val, y_val,\n             n_trials: int = 25) -> Dict:\n    if not (HAS_OPTUNA and HAS_XGB):\n        return {}\n\n    def objective(trial):\n        p = {\n            \"n_estimators\":    trial.suggest_int(\"n_est\",  50, 300),\n            \"max_depth\":       trial.suggest_int(\"depth\",   3,  10),\n            \"learning_rate\":   trial.suggest_float(\"lr\",  0.01, 0.3, log=True),\n            \"subsample\":       trial.suggest_float(\"sub\",  0.5, 1.0),\n            \"colsample_bytree\":trial.suggest_float(\"col\",  0.5, 1.0),\n            \"reg_alpha\":       trial.suggest_float(\"a\",   1e-4, 10.0, log=True),\n            \"reg_lambda\":      trial.suggest_float(\"l\",   1e-4, 10.0, log=True),\n        }\n        m = build_xgb(p)\n        m.fit(X_tr, y_tr,\n              eval_set=[(X_val, y_val)],\n              verbose=False)\n        pred = m.predict(X_val)\n        return mean_squared_error(y_val, pred)\n\n    study = optuna.create_study(direction=\"minimize\",\n                                sampler=optuna.samplers.TPESampler(seed=SEED))\n    study.optimize(objective, n_trials=n_trials, show_progress_bar=False)\n    log.info(f\"  XGB best MSE={study.best_value:.6f}  params={study.best_params}\")\n    return study.best_params\n\n\ndef tune_lgb(X_tr, y_tr, X_val, y_val,\n             n_trials: int = 25) -> Dict:\n    if not (HAS_OPTUNA and HAS_LGB):\n        return {}\n\n    def objective(trial):\n        p = {\n            \"n_estimators\":  trial.suggest_int(\"n_est\",  50, 300),\n            \"num_leaves\":    trial.suggest_int(\"leaves\", 15, 127),\n            \"learning_rate\": trial.suggest_float(\"lr\", 0.01, 0.3, log=True),\n            \"subsample\":     trial.suggest_float(\"sub\", 0.5, 1.0),\n            \"colsample_bytree\":trial.suggest_float(\"col\", 0.5, 1.0),\n            \"reg_alpha\":     trial.suggest_float(\"a\", 1e-4, 10.0, log=True),\n            \"reg_lambda\":    trial.suggest_float(\"l\", 1e-4, 10.0, log=True),\n        }\n        m = build_lgb(p)\n        m.fit(X_tr, y_tr,\n              eval_set=[(X_val, y_val)],\n              callbacks=[lgb.early_stopping(20, verbose=False),\n                         lgb.log_evaluation(-1)])\n        pred = m.predict(X_val)\n        return mean_squared_error(y_val, pred)\n\n    study = optuna.create_study(direction=\"minimize\",\n                                sampler=optuna.samplers.TPESampler(seed=SEED))\n    study.optimize(objective, n_trials=n_trials, show_progress_bar=False)\n    log.info(f\"  LGB best MSE={study.best_value:.6f}  params={study.best_params}\")\n    return study.best_params\n\n\n# ─── PyTorch MLP with anti-overfitting ──────────────────────────────────────\n\nclass AttentionBlock(torch.nn.Module if HAS_TORCH else object):\n    \"\"\"Self-attention over feature dimension.\"\"\"\n    def __init__(self, in_dim: int):\n        super().__init__()\n        self.attn = torch.nn.Linear(in_dim, in_dim)\n        self.norm = torch.nn.LayerNorm(in_dim)\n\n    def forward(self, x):\n        w = torch.sigmoid(self.attn(x))\n        return self.norm(x * w)\n\n\nclass RNAMlp(torch.nn.Module if HAS_TORCH else object):\n    \"\"\"Deep MLP with BN, Dropout, Attention, Residual — for regression.\"\"\"\n\n    def __init__(self, in_dim: int,\n                 hidden: Tuple = (256, 128, 64),\n                 dropout: float = 0.35,\n                 out_dim: int = 1):\n        super().__init__()\n        layers = []\n        prev = in_dim\n        for h in hidden:\n            layers += [\n                torch.nn.Linear(prev, h),\n                torch.nn.BatchNorm1d(h),\n                torch.nn.GELU(),\n                torch.nn.Dropout(dropout),\n            ]\n            prev = h\n        self.net  = torch.nn.Sequential(*layers)\n        self.attn = AttentionBlock(prev)\n        self.head = torch.nn.Linear(prev, out_dim)\n        self._init_weights()\n\n    def _init_weights(self):\n        for m in self.modules():\n            if isinstance(m, torch.nn.Linear):\n                torch.nn.init.kaiming_normal_(m.weight, nonlinearity=\"relu\")\n                if m.bias is not None:\n                    torch.nn.init.zeros_(m.bias)\n\n    def forward(self, x):\n        h = self.net(x)\n        h = self.attn(h)\n        return self.head(h)\n\n\nclass EarlyStopping:\n    def __init__(self, patience: int = 15, min_delta: float = 1e-5):\n        self.patience   = patience\n        self.min_delta  = min_delta\n        self.best_score = None\n        self.counter    = 0\n        self.stop       = False\n\n    def step(self, score: float, model, path: Path) -> None:\n        if self.best_score is None or score < self.best_score - self.min_delta:\n            self.best_score = score\n            self.counter    = 0\n            if HAS_TORCH:\n                torch.save(model.state_dict(), path)\n        else:\n            self.counter += 1\n            if self.counter >= self.patience:\n                self.stop = True\n\n\ndef train_mlp(X_tr: np.ndarray, y_tr: np.ndarray,\n              X_val: np.ndarray, y_val: np.ndarray,\n              epochs: int = 120,\n              batch_size: int = 256,\n              lr: float = 3e-4,\n              weight_decay: float = 1e-3) -> Tuple[\"RNAMlp\", Dict]:\n    if not HAS_TORCH:\n        return None, {}\n    MODEL_DIR.mkdir(parents=True, exist_ok=True)\n    ckpt_path = MODEL_DIR / \"best_mlp.pt\"\n\n    # Standardise target to zero-mean unit-variance so the MLP converges\n    # regardless of the target's scale or range.\n    y_mean = float(y_tr.mean()); y_std = float(y_tr.std()) + 1e-8\n    y_tr_s  = (y_tr  - y_mean) / y_std\n    y_val_s = (y_val - y_mean) / y_std\n\n    in_dim = X_tr.shape[1]\n    model  = RNAMlp(in_dim).to(DEVICE)\n    opt    = torch.optim.AdamW(model.parameters(), lr=lr,\n                               weight_decay=weight_decay)\n    sched  = torch.optim.lr_scheduler.ReduceLROnPlateau(\n                opt, patience=8, factor=0.5)\n    es     = EarlyStopping(patience=15)\n    loss_fn = torch.nn.HuberLoss()\n\n    Xtr = torch.tensor(X_tr,   dtype=torch.float32)\n    ytr = torch.tensor(y_tr_s.reshape(-1,1), dtype=torch.float32)\n    Xvl = torch.tensor(X_val,  dtype=torch.float32).to(DEVICE)\n    yvl = torch.tensor(y_val_s.reshape(-1,1), dtype=torch.float32).to(DEVICE)\n\n    ds  = torch.utils.data.TensorDataset(Xtr, ytr)\n    dl  = torch.utils.data.DataLoader(ds, batch_size=batch_size,\n                                       shuffle=True, pin_memory=True)\n    history = {\"loss\": [], \"val_loss\": [], \"accuracy\": [], \"val_accuracy\": []}\n\n    for epoch in tqdm(range(epochs), desc=\"MLP training\", leave=False):\n        model.train()\n        ep_loss = 0.0\n        for xb, yb in dl:\n            xb, yb = xb.to(DEVICE), yb.to(DEVICE)\n            opt.zero_grad()\n            pred = model(xb)\n            loss = loss_fn(pred, yb)\n            loss.backward()\n            torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)\n            opt.step()\n            ep_loss += loss.item() * len(xb)\n        ep_loss /= len(Xtr)\n\n        model.eval()\n        with torch.no_grad():\n            vl_pred = model(Xvl)\n            vl_loss = loss_fn(vl_pred, yvl).item()\n        sched.step(vl_loss)\n\n        # Denormalise back to original scale before computing R²\n        yvl_np = yvl.cpu().numpy().ravel() * y_std + y_mean\n        vl_np  = vl_pred.cpu().numpy().ravel() * y_std + y_mean\n        vl_r2  = max(0.0, r2_score(yvl_np, vl_np))\n        tr_r2  = max(0.0, 1 - ep_loss * y_std**2 / (np.var(y_tr) + 1e-9))\n\n        history[\"loss\"].append(ep_loss)\n        history[\"val_loss\"].append(vl_loss)\n        history[\"accuracy\"].append(tr_r2)\n        history[\"val_accuracy\"].append(vl_r2)\n\n        es.step(vl_loss, model, ckpt_path)\n        if es.stop:\n            log.info(f\"  EarlyStopping at epoch {epoch}\")\n            break\n\n    # load best checkpoint\n    model.load_state_dict(torch.load(ckpt_path, map_location=DEVICE))\n    log.info(f\"  Best val R²={max(history['val_accuracy']):.4f}\")\n    # Store normalisation stats on the model so callers can denormalise\n    model._y_mean = y_mean\n    model._y_std  = y_std\n    force_cleanup(ds, dl, Xtr, ytr, Xvl, yvl)\n    return model, history\n\n\n# ─── Stacking Ensemble ───────────────────────────────────────────────────────\n\ndef build_stacking_ensemble(base_models: List,\n                              meta_model=None) -> StackingRegressor:\n    estimators = [(f\"m{i}\", m) for i, m in enumerate(base_models) if m is not None]\n    if meta_model is None:\n        meta_model = Ridge(alpha=1.0)\n    return StackingRegressor(estimators=estimators,\n                              final_estimator=meta_model,\n                              cv=3, n_jobs=-1)\n\n\ndef evaluate_model(model, X: np.ndarray, y: np.ndarray,\n                   split_name: str = \"test\",\n                   do_plots: bool = True) -> Dict[str, float]:\n    pred = model.predict(X)\n    mse  = mean_squared_error(y, pred)\n    mae  = mean_absolute_error(y, pred)\n    r2   = r2_score(y, pred)\n    rmse = np.sqrt(mse)\n    log.info(f\"  [{split_name}] RMSE={rmse:.4f}  MAE={mae:.4f}  R²={r2:.4f}\")\n    if do_plots:\n        plot_prediction_scatter(y, pred, split_name)\n        plot_residuals(y, pred, split_name)\n    return {\"rmse\": rmse, \"mae\": mae, \"r2\": r2}\n\n\n# ──────────────────────────── Cross-Validation ──────────────────────────────\n\ndef cross_validate_model(model_factory, X: np.ndarray, y: np.ndarray,\n                          n_splits: int = 5) -> Dict:\n    kf = KFold(n_splits=n_splits, shuffle=True, random_state=SEED)\n    scores = {\"rmse\": [], \"mae\": [], \"r2\": []}\n    for fold, (tr_idx, vl_idx) in enumerate(\n        tqdm(kf.split(X), total=n_splits, desc=\"CV folds\", leave=False)\n    ):\n        m = model_factory()\n        m.fit(X[tr_idx], y[tr_idx])\n        p = m.predict(X[vl_idx])\n        scores[\"rmse\"].append(np.sqrt(mean_squared_error(y[vl_idx], p)))\n        scores[\"mae\"].append(mean_absolute_error(y[vl_idx], p))\n        scores[\"r2\"].append(r2_score(y[vl_idx], p))\n        force_cleanup(m)\n\n    summary = {k: (float(np.mean(v)), float(np.std(v))) for k, v in scores.items()}\n    log.info(\"  CV results:\")\n    for k, (mu, sd) in summary.items():\n        log.info(f\"    {k}: {mu:.4f} ± {sd:.4f}\")\n    return summary\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 8: TM-SCORE APPROXIMATION\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef _d0(L_ref: int) -> float:\n    if L_ref >= 30:\n        return 1.24 * (L_ref - 15) ** (1/3) - 1.8\n    elif L_ref >= 24:\n        return 0.7\n    elif L_ref >= 20:\n        return 0.6\n    elif L_ref >= 16:\n        return 0.5\n    elif L_ref >= 12:\n        return 0.4\n    return 0.3\n\n\ndef tm_score_approx(pred_coords: np.ndarray,\n                     ref_coords:  np.ndarray) -> float:\n    \"\"\"\n    Fast approximate TM-score (no superposition — just aligned distance).\n    For post-hoc evaluation against validation ground truth.\n    \"\"\"\n    L_ref  = len(ref_coords)\n    L_aln  = min(len(pred_coords), L_ref)\n    d0_val = _d0(L_ref) + 1e-9\n    dists  = np.linalg.norm(pred_coords[:L_aln] - ref_coords[:L_aln], axis=1)\n    return float(np.mean(1.0 / (1.0 + (dists / d0_val)**2)))\n\n\ndef best_of_n_tm(multi_pred: np.ndarray,\n                  ref_coords:  np.ndarray) -> float:\n    \"\"\"Best-of-N TM-score.\"\"\"\n    return max(tm_score_approx(multi_pred[s], ref_coords)\n               for s in range(multi_pred.shape[0]))\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 9: TBM — TEMPLATE-BASED MODELLING (unchanged core, memory-safe)\n# ══════════════════════════════════════════════════════════════════════════════\n\nif HAS_BIOPYTHON:\n    def _make_aligner() -> PairwiseAligner:\n        al = PairwiseAligner()\n        al.mode                          = \"global\"\n        al.match_score                   = 2\n        al.mismatch_score                = -1.5\n        al.open_gap_score                = -8\n        al.extend_gap_score              = -0.4\n        al.query_left_open_gap_score     = -8\n        al.query_left_extend_gap_score   = -0.4\n        al.query_right_open_gap_score    = -8\n        al.query_right_extend_gap_score  = -0.4\n        al.target_left_open_gap_score    = -8\n        al.target_left_extend_gap_score  = -0.4\n        al.target_right_open_gap_score   = -8\n        al.target_right_extend_gap_score = -0.4\n        return al\n    _aligner = _make_aligner()\nelse:\n    _aligner = None\n\n\ndef process_labels(labels_df: pd.DataFrame) -> Dict:\n    coords: Dict = {}\n    prefixes = labels_df[\"ID\"].str.rsplit(\"_\", n=1).str[0]\n    for prefix, grp in tqdm(labels_df.groupby(prefixes),\n                             desc=\"Building coord dict\", leave=False):\n        coords[prefix] = grp.sort_values(\"resid\")[[\"x_1\",\"y_1\",\"z_1\"]].values\n    return coords\n\n\ndef parse_stoichiometry(stoich: str) -> list:\n    if pd.isna(stoich) or str(stoich).strip() == \"\":\n        return []\n    return [(ch.strip(), int(cnt))\n            for part in str(stoich).split(\";\")\n            for ch, cnt in [part.split(\":\")]]\n\n\ndef parse_fasta(fasta_content: str) -> Dict:\n    out, cur, parts = {}, None, []\n    for line in str(fasta_content).splitlines():\n        line = line.strip()\n        if not line:\n            continue\n        if line.startswith(\">\"):\n            if cur is not None:\n                out[cur] = \"\".join(parts)\n            cur = line[1:].split()[0]; parts = []\n        else:\n            parts.append(line.replace(\" \",\"\"))\n    if cur is not None:\n        out[cur] = \"\".join(parts)\n    return out\n\n\ndef get_chain_segments(row) -> List[Tuple]:\n    seq    = row[\"sequence\"]\n    stoich = row.get(\"stoichiometry\",\"\")\n    all_sq = row.get(\"all_sequences\",\"\")\n    if (pd.isna(stoich) or pd.isna(all_sq)\n            or str(stoich).strip() == \"\" or str(all_sq).strip() == \"\"):\n        return [(0, len(seq))]\n    try:\n        chain_dict = parse_fasta(all_sq)\n        order      = parse_stoichiometry(stoich)\n        segs, pos  = [], 0\n        for ch, cnt in order:\n            base = chain_dict.get(ch)\n            if base is None:\n                return [(0, len(seq))]\n            for _ in range(cnt):\n                segs.append((pos, pos+len(base)))\n                pos += len(base)\n        return segs if pos == len(seq) else [(0, len(seq))]\n    except Exception:\n        return [(0, len(seq))]\n\n\ndef build_segments_map(df: pd.DataFrame) -> Tuple[Dict, Dict]:\n    seg_map, stoich_map = {}, {}\n    for _, r in df.iterrows():\n        tid             = r[\"target_id\"]\n        seg_map[tid]    = get_chain_segments(r)\n        raw_s           = r.get(\"stoichiometry\",\"\")\n        stoich_map[tid] = \"\" if pd.isna(raw_s) else str(raw_s)\n    return seg_map, stoich_map\n\n\ndef _build_aligned_strings(query_seq, tmpl_seq, alignment):\n    q_segs, t_segs = alignment.aligned\n    aq, at, qi, ti = [], [], 0, 0\n    for (qs,qe),(ts,te) in zip(q_segs,t_segs):\n        while qi<qs: aq.append(query_seq[qi]); at.append(\"-\");             qi+=1\n        while ti<ts: aq.append(\"-\");            at.append(tmpl_seq[ti]);   ti+=1\n        for qp,tp in zip(range(qs,qe),range(ts,te)):\n            aq.append(query_seq[qp]); at.append(tmpl_seq[tp])\n        qi,ti = qe,te\n    while qi<len(query_seq): aq.append(query_seq[qi]); at.append(\"-\");            qi+=1\n    while ti<len(tmpl_seq):  aq.append(\"-\");            at.append(tmpl_seq[ti]);  ti+=1\n    return \"\".join(aq),\"\".join(at)\n\n\ndef find_similar_sequences(query_seq: str,\n                            train_df: pd.DataFrame,\n                            train_coords: Dict,\n                            top_n: int = 30) -> List:\n    if not HAS_BIOPYTHON:\n        return []\n    results = []\n    for _, row in train_df.iterrows():\n        tid, tseq = row[\"target_id\"], row[\"sequence\"]\n        if tid not in train_coords:\n            continue\n        if abs(len(tseq)-len(query_seq))/max(len(tseq),len(query_seq)) > 0.3:\n            continue\n        aln       = next(iter(_aligner.align(query_seq, tseq)))\n        norm_s    = aln.score / (2*min(len(query_seq),len(tseq)))\n        identical = sum(\n            1 for (qs,qe),(ts,te) in zip(*aln.aligned)\n            for qp,tp in zip(range(qs,qe),range(ts,te))\n            if query_seq[qp]==tseq[tp]\n        )\n        pct_id = 100*identical/len(query_seq)\n        results.append((tid, tseq, norm_s, train_coords[tid], pct_id))\n    results.sort(key=lambda x: x[2], reverse=True)\n    return results[:top_n]\n\n\ndef adapt_template_to_query(query_seq: str,\n                             tmpl_seq: str,\n                             tmpl_coords: np.ndarray) -> np.ndarray:\n    if not HAS_BIOPYTHON:\n        return generate_rna_structure(query_seq)\n    aln        = next(iter(_aligner.align(query_seq, tmpl_seq)))\n    new_coords = np.full((len(query_seq),3), np.nan)\n    for (qs,qe),(ts,te) in zip(*aln.aligned):\n        chunk = tmpl_coords[ts:te]\n        if len(chunk)==(qe-qs):\n            new_coords[qs:qe] = chunk\n    for i in range(len(new_coords)):\n        if np.isnan(new_coords[i,0]):\n            pv = next((j for j in range(i-1,-1,-1) if not np.isnan(new_coords[j,0])), -1)\n            nv = next((j for j in range(i+1,len(new_coords)) if not np.isnan(new_coords[j,0])), -1)\n            if pv>=0 and nv>=0:\n                w = (i-pv)/(nv-pv)\n                new_coords[i] = (1-w)*new_coords[pv]+w*new_coords[nv]\n            elif pv>=0:  new_coords[i] = new_coords[pv]+[3,0,0]\n            elif nv>=0:  new_coords[i] = new_coords[nv]+[3,0,0]\n            else:        new_coords[i] = [i*3,0,0]\n    return np.nan_to_num(new_coords)\n\n\ndef adaptive_rna_constraints(coords: np.ndarray,\n                               tid: str,\n                               segments_map: Dict,\n                               confidence: float = 1.0,\n                               passes: int = 2) -> np.ndarray:\n    X        = coords.copy()\n    segments = segments_map.get(tid, [(0,len(X))])\n    strength = max(0.75*(1.0-min(confidence,0.97)), 0.02)\n    for _ in range(passes):\n        for s,e in segments:\n            C = X[s:e]; L = e-s\n            if L<3: continue\n            d  = C[1:]-C[:-1]; dist = np.linalg.norm(d,axis=1)+1e-6\n            adj = d*((5.95-dist)/dist)[:,None]*(0.22*strength)\n            C[:-1] -= adj; C[1:]  += adj\n            d2  = C[2:]-C[:-2]; d2n = np.linalg.norm(d2,axis=1)+1e-6\n            adj2 = d2*((10.2-d2n)/d2n)[:,None]*(0.10*strength)\n            C[:-2] -= adj2; C[2:] += adj2\n            C[1:-1] += (0.06*strength)*(0.5*(C[:-2]+C[2:])-C[1:-1])\n            if L>=25:\n                idx  = np.linspace(0,L-1,min(L,160)).astype(int) if L>220 else np.arange(L)\n                P    = C[idx]; diff = P[:,None,:]-P[None,:,:]\n                dm   = np.linalg.norm(diff,axis=2)+1e-6\n                sep  = np.abs(idx[:,None]-idx[None,:])\n                mask = (sep>2) & (dm<3.2)\n                if np.any(mask):\n                    vec = (diff*((3.2-dm)/dm)[:,:,None]*mask[:,:,None]).sum(axis=1)\n                    C[idx] += (0.015*strength)*vec\n            X[s:e] = C\n    return X\n\n\ndef _rotmat(axis, ang):\n    a = np.asarray(axis,float); a /= np.linalg.norm(a)+1e-12\n    x,y,z = a; c,s = np.cos(ang),np.sin(ang); CC = 1-c\n    return np.array([[c+x*x*CC, x*y*CC-z*s, x*z*CC+y*s],\n                     [y*x*CC+z*s, c+y*y*CC, y*z*CC-x*s],\n                     [z*x*CC-y*s, z*y*CC+x*s, c+z*z*CC]])\n\n\ndef apply_hinge(coords, seg, rng, deg=22):\n    s,e = seg; L = e-s\n    if L<30: return coords\n    pivot = s+int(rng.integers(10,L-10))\n    R  = _rotmat(rng.normal(size=3), np.deg2rad(float(rng.uniform(-deg,deg))))\n    X  = coords.copy(); p0 = X[pivot].copy()\n    X[pivot+1:e] = (X[pivot+1:e]-p0)@R.T+p0\n    return X\n\n\ndef jitter_chains(coords, segs, rng, deg=12, trans=1.5):\n    X = coords.copy(); gc_ = X.mean(0,keepdims=True)\n    for s,e in segs:\n        R     = _rotmat(rng.normal(size=3),np.deg2rad(float(rng.uniform(-deg,deg))))\n        shift = rng.normal(size=3); shift = shift/(np.linalg.norm(shift)+1e-12)*float(rng.uniform(0,trans))\n        c     = X[s:e].mean(0,keepdims=True)\n        X[s:e] = (X[s:e]-c)@R.T+c+shift\n    X -= X.mean(0,keepdims=True)-gc_\n    return X\n\n\ndef smooth_wiggle(coords, segs, rng, amp=0.8):\n    X = coords.copy()\n    for s,e in segs:\n        L = e-s\n        if L<20: continue\n        ctrl = np.linspace(0,L-1,6); disp = rng.normal(0,amp,(6,3)); t = np.arange(L)\n        X[s:e] += np.vstack([np.interp(t,ctrl,disp[:,k]) for k in range(3)]).T\n    return X\n\n\ndef generate_rna_structure(sequence: str, seed: Optional[int] = None) -> np.ndarray:\n    \"\"\"Idealized A-form helix fallback.\"\"\"\n    if seed is not None: np.random.seed(seed)\n    n = len(sequence); coords = np.zeros((n,3))\n    for i in range(n):\n        ang = i*0.6\n        coords[i] = [10.0*np.cos(ang), 10.0*np.sin(ang), i*2.5]\n    return coords\n\n\n# ──────────────────────────────── Kabsch / stitch ───────────────────────────\n\ndef kabsch_align(P, Q):\n    cp,cq = P.mean(0), Q.mean(0)\n    Pc,Qc = P-cp, Q-cq\n    H = Pc.T@Qc; U,_,Vt = np.linalg.svd(H)\n    d = np.linalg.det(Vt.T@U.T); S = np.eye(3)\n    if d<0: S[2,2] = -1\n    R = Vt.T@S@U.T; t = cq-R@cp\n    return R,t\n\n\ndef split_into_chunks(seq_len: int, max_len: int, overlap: int) -> List[Tuple]:\n    if seq_len <= max_len:\n        return [(0, seq_len)]\n    chunks, step, pos = [], max_len-overlap, 0\n    while pos < seq_len:\n        end = min(pos+max_len, seq_len)\n        chunks.append((pos,end))\n        if end==seq_len: break\n        pos += step\n    return chunks\n\n\ndef stitch_chunk_coords(chunk_coords_list, chunk_ranges, seq_len) -> np.ndarray:\n    if len(chunk_coords_list)==1:\n        c = chunk_coords_list[0]\n        out = np.zeros((seq_len,3), dtype=c.dtype)\n        out[:min(len(c),seq_len)] = c[:min(len(c),seq_len)]\n        return out\n    aligned = [chunk_coords_list[0].copy()]\n    for i in range(1,len(chunk_coords_list)):\n        ps,pe = chunk_ranges[i-1]; cs,ce = chunk_ranges[i]\n        ov_s,ov_e = cs, min(pe,ce); ov_len = ov_e-ov_s\n        if ov_len<3: aligned.append(chunk_coords_list[i].copy()); continue\n        prev_ov = aligned[i-1][ov_s-ps:ov_e-ps]\n        cur_ov  = chunk_coords_list[i][ov_s-cs:ov_e-cs]\n        valid   = ~(np.isnan(prev_ov).any(1)|np.isnan(cur_ov).any(1))\n        if valid.sum()<3: aligned.append(chunk_coords_list[i].copy()); continue\n        R,t = kabsch_align(cur_ov[valid],prev_ov[valid])\n        aligned.append((chunk_coords_list[i]@R.T)+t)\n    full = np.zeros((seq_len,3),dtype=np.float64)\n    weights = np.zeros(seq_len,dtype=np.float64)\n    for i,((s,e),coords) in enumerate(zip(chunk_ranges,aligned)):\n        ae  = min(s+len(coords),seq_len); ul = ae-s\n        w   = np.ones(ul,dtype=np.float64)\n        if i>0:\n            ov_e = min(chunk_ranges[i-1][1],e); rl = ov_e-s\n            if rl>0: w[:rl] = np.linspace(0,1,rl)\n        if i<len(chunk_ranges)-1:\n            ns = chunk_ranges[i+1][0]; rs = ns-s; rl = ae-ns\n            if rl>0 and rs<ul: w[rs:ul] = np.linspace(1,0,rl)\n        full[s:ae]    += coords[:ul]*w[:,None]\n        weights[s:ae] += w\n    mask = weights>0; full[mask] /= weights[mask,None]\n    return full\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 10: PROTENIX INTERFACE (GPU/CPU transparent)\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef build_input_json(df: pd.DataFrame, json_path: str) -> None:\n    data = [\n        {\"name\": row[\"target_id\"], \"covalent_bonds\": [],\n         \"sequences\": [{\"rnaSequence\": {\"sequence\": row[\"sequence\"],\"count\":1}}]}\n        for _,row in df.iterrows()\n    ]\n    with open(json_path,\"w\",encoding=\"utf-8\") as f:\n        json.dump(data,f)\n\n\ndef build_configs(input_json_path: str, dump_dir: str, model_name: str):\n    from configs.configs_base import configs as configs_base\n    from configs.configs_data import data_configs\n    from configs.configs_inference import inference_configs\n    from configs.configs_model_type import model_configs\n    from protenix.config.config import parse_configs\n    base = {**configs_base, **{\"data\": data_configs}, **inference_configs}\n    def _upd(t,p):\n        for k,v in p.items():\n            if isinstance(v,dict) and k in t and isinstance(t[k],dict): _upd(t[k],v)\n            else: t[k] = v\n    _upd(base, model_configs[model_name])\n    arg_str = \" \".join([\n        f\"--model_name {model_name}\",\n        f\"--input_json_path {input_json_path}\",\n        f\"--dump_dir {dump_dir}\",\n        f\"--use_msa {USE_MSA}\",\n        f\"--use_template {USE_TEMPLATE}\",\n        f\"--use_rna_msa {USE_RNA_MSA}\",\n        f\"--sample_diffusion.N_sample {MODEL_N_SAMPLE}\",\n        f\"--seeds {SEED}\",\n    ])\n    from protenix.config.config import parse_configs\n    return parse_configs(configs=base, arg_str=arg_str, fill_required_with_null=True)\n\n\ndef ensure_required_files(root_dir: str) -> None:\n    for p,name in [\n        (Path(root_dir)/\"checkpoint\"/f\"{MODEL_NAME}.pt\",        \"checkpoint\"),\n        (Path(root_dir)/\"common\"/\"components.cif\",             \"CCD file\"),\n        (Path(root_dir)/\"common\"/\"components.cif.rdkit_mol.pkl\",\"CCD cache\"),\n    ]:\n        if not p.exists():\n            raise FileNotFoundError(f\"Missing {name}: {p}\")\n\n\ndef coords_to_rows(target_id: str,\n                    seq: str,\n                    coords: np.ndarray) -> List[Dict]:\n    rows = []\n    for i in range(len(seq)):\n        row = {\"ID\": f\"{target_id}_{i+1}\", \"resname\": seq[i], \"resid\": i+1}\n        for s in range(N_SAMPLE):\n            x,y,z = (coords[s,i].tolist() if s<coords.shape[0] and i<coords.shape[1]\n                      else [0.,0.,0.])\n            row[f\"x_{s+1}\"] = float(x); row[f\"y_{s+1}\"] = float(y); row[f\"z_{s+1}\"] = float(z)\n        rows.append(row)\n    return rows\n\n\ndef pad_samples(coords: np.ndarray, n: int) -> np.ndarray:\n    if coords.shape[0] >= n: return coords[:n]\n    if coords.shape[0] == 0: return np.zeros((n,coords.shape[1],3),dtype=coords.dtype)\n    extra = np.repeat(coords[:1], n-coords.shape[0], axis=0)\n    return np.concatenate([coords,extra],axis=0)\n\n\n# ══════════════════════════════════════════════════════════════════════════════\n# SECTION 11: MAIN PIPELINE\n# ══════════════════════════════════════════════════════════════════════════════\n\ndef main() -> None:\n    global _PIPELINE_START\n    t_total = time.time()\n    _PIPELINE_START = t_total\n    log.info(f\"  Budget: {TOTAL_BUDGET_S//60} min total | Protenix timeout: {PROTENIX_SAMPLE_TIMEOUT_S}s/sample\")\n    PLOT_DIR.mkdir(parents=True, exist_ok=True)\n    MODEL_DIR.mkdir(parents=True, exist_ok=True)\n\n    log.info(\"=\" * 68)\n    log.info(\"  RNA 3D FOLDING — ENHANCED ML PIPELINE\")\n    log.info(f\"  Device: {DEVICE}  |  Memory limit: {MEMORY_LIMIT_MB} MB\")\n    log.info(\"=\" * 68)\n\n    # ── Protenix setup ───────────────────────────────────────────────────────\n    code_dir = DEFAULT_CODE_DIR\n    root_dir = DEFAULT_ROOT_DIR\n    protenix_ok = False\n    if os.path.isdir(code_dir) and USE_PROTENIX:\n        os.environ[\"PROTENIX_ROOT_DIR\"] = root_dir\n        sys.path.insert(0, code_dir)\n        try:\n            ensure_required_files(root_dir)\n            protenix_ok = True\n            log.info(\"✓ Protenix directory validated.\")\n        except FileNotFoundError as e:\n            log.warning(f\"Protenix files missing — deep-learning phase disabled.\\n  {e}\")\n    else:\n        log.info(\"Protenix code dir not found — TBM + fallback only.\")\n\n    seed_everything(SEED)\n\n    # ── 1. Load data ─────────────────────────────────────────────────────────\n    log.info(\"\\n[1/9] Loading data …\")\n    test_df    = load_csv_chunked(DEFAULT_TEST_CSV)\n    train_seqs = load_csv_chunked(DEFAULT_TRAIN_CSV, max_rows=MAX_TEMPLATE_POOL)\n    val_seqs   = load_csv_chunked(DEFAULT_VAL_CSV,   max_rows=MAX_TEMPLATE_POOL//2)\n    train_lbl  = load_csv_chunked(DEFAULT_TRAIN_LBL)\n    val_lbl    = load_csv_chunked(DEFAULT_VAL_LBL)\n\n    if test_df.empty:\n        log.error(\"test_sequences.csv not found or empty — aborting.\")\n        return\n\n    combined_seqs   = pd.concat([train_seqs, val_seqs], ignore_index=True)\n    combined_labels = pd.concat([train_lbl,  val_lbl],  ignore_index=True)\n\n    log.info(f\"  Test      : {len(test_df):>6d} targets\")\n    log.info(f\"  Templates : {len(combined_seqs):>6d} sequences  \"\n             f\"| {len(combined_labels):>8d} label rows\")\n    force_cleanup(train_lbl, val_lbl)\n\n    # ── 2. EDA ───────────────────────────────────────────────────────────────\n    log.info(\"\\n[2/9] EDA …\")\n    eda_sequence_stats(combined_seqs, \"templates\")\n    eda_sequence_stats(test_df,       \"test\")\n    if not combined_labels.empty:\n        eda_label_stats(combined_labels, \"train+val\")\n    memory_safe(\"after EDA\")\n\n    # ── 3. Feature engineering ───────────────────────────────────────────────\n    log.info(\"\\n[3/9] Feature engineering …\")\n    feat_all = build_feature_dataframe(combined_seqs, \"Features — templates\")\n    feat_test = build_feature_dataframe(test_df, \"Features — test\")\n    log.info(f\"  Feature columns: {feat_all.shape[1]-1}\")\n\n    # ── 4. Build coordinates dict for TBM ───────────────────────────────────\n    log.info(\"\\n[4/9] Building TBM coordinate cache …\")\n    train_coords  = process_labels(combined_labels)\n    segments_map, _ = build_segments_map(test_df)\n    force_cleanup(combined_labels)\n\n    # ── 5. Data splits for ML meta-model ────────────────────────────────────\n    log.info(\"\\n[5/9] Splitting data for ML meta-model …\")\n    # target_col must be excluded from feature_cols to avoid duplicate columns\n    # when the split dataframe is built via np.hstack([X_all, y_all])\n    target_col   = \"gc_content\"\n    feature_cols = [c for c in feat_all.columns\n                    if c not in (\"target_id\", target_col)]\n    merged       = feat_all.copy()\n    merged       = merged.dropna(subset=[target_col])\n    X_all        = merged[feature_cols].values.astype(\"float32\")\n    y_all        = merged[target_col].values.astype(\"float32\")\n\n    # Impute NaN\n    imputer = SimpleImputer(strategy=\"median\")\n    X_all   = imputer.fit_transform(X_all)\n\n    # Outlier detection & sample weights\n    outlier_mask = detect_outliers_iqr(X_all)\n    sample_w     = sample_weights_from_outliers(X_all, outlier_mask)\n    log.info(f\"  Outliers detected: {outlier_mask.sum()} / {len(X_all)}\")\n\n    tr, va, te, ho = four_way_split(\n        pd.DataFrame(np.hstack([X_all, y_all.reshape(-1,1)]),\n                     columns=feature_cols+[target_col])\n    )\n    verify_data_splits(tr, va, te, ho)\n    plot_split_summary(tr, va, te, ho)\n\n    # Retrieve per-split sample weights using the __orig_idx sentinel BEFORE dropping it\n    train_orig_idx = tr[\"__orig_idx\"].values.astype(int)\n    sample_w_train = sample_w[train_orig_idx]          # correctly aligned weights\n\n    # Drop the sentinel column injected by four_way_split before using as features\n    for _split in (tr, va, te, ho):\n        _split.drop(columns=[\"__orig_idx\"], errors=\"ignore\", inplace=True)\n\n    Xtr  = tr[feature_cols].values; ytr = tr[target_col].values\n    Xva  = va[feature_cols].values; yva = va[target_col].values\n    Xte  = te[feature_cols].values; yte = te[target_col].values\n    Xho  = ho[feature_cols].values; yho = ho[target_col].values\n\n    # Scale\n    scaler = RobustScaler()\n    Xtr    = scaler.fit_transform(Xtr)\n    Xva    = scaler.transform(Xva)\n    Xte    = scaler.transform(Xte)\n    Xho    = scaler.transform(Xho)\n\n    # PCA / t-SNE visualisation\n    plot_pca_tsne(Xtr, ytr, \"train\")\n    force_cleanup(tr, va, te, ho)\n\n    # ── 6. Hyper-parameter tuning & model training ───────────────────────────\n    log.info(\"\\n[6/9] Model training + hyper-parameter tuning …\")\n    best_models, all_scores = {}, {}\n\n    # Random Forest (baseline)\n    log.info(\"  Training Random Forest …\")\n    rf = build_rf()\n    rf.fit(Xtr, ytr, sample_weight=sample_w_train)\n    all_scores[\"rf\"] = {\n        \"train\":   evaluate_model(rf, Xtr, ytr, \"rf_train\",   do_plots=False)[\"r2\"],\n        \"val\":     evaluate_model(rf, Xva, yva, \"rf_val\",     do_plots=False)[\"r2\"],\n        \"test\":    evaluate_model(rf, Xte, yte, \"rf_test\",    do_plots=True )[\"r2\"],\n        \"holdout\": evaluate_model(rf, Xho, yho, \"rf_holdout\", do_plots=True )[\"r2\"],\n    }\n    plot_generalization_gap(all_scores[\"rf\"], \"RF Generalization Gap\")\n    plot_feature_importance(rf.feature_importances_, feature_cols, \"Random Forest\")\n    best_models[\"rf\"] = rf\n    force_cleanup()\n\n    # XGBoost\n    if HAS_XGB:\n        log.info(\"  Tuning XGBoost …\")\n        xgb_params = tune_xgb(Xtr, ytr, Xva, yva, n_trials=20 if HAS_OPTUNA else 0)\n        xgb_m = build_xgb(xgb_params)\n        xgb_m.fit(Xtr, ytr, eval_set=[(Xva, yva)], verbose=False)\n        all_scores[\"xgb\"] = {\n            \"train\":   evaluate_model(xgb_m, Xtr, ytr, \"xgb_train\",   do_plots=False)[\"r2\"],\n            \"val\":     evaluate_model(xgb_m, Xva, yva, \"xgb_val\",     do_plots=False)[\"r2\"],\n            \"test\":    evaluate_model(xgb_m, Xte, yte, \"xgb_test\",    do_plots=True )[\"r2\"],\n            \"holdout\": evaluate_model(xgb_m, Xho, yho, \"xgb_holdout\", do_plots=True )[\"r2\"],\n        }\n        plot_generalization_gap(all_scores[\"xgb\"], \"XGB Generalization Gap\")\n        plot_feature_importance(xgb_m.feature_importances_, feature_cols, \"XGBoost\")\n        best_models[\"xgb\"] = xgb_m\n        force_cleanup()\n\n    # LightGBM\n    if HAS_LGB:\n        log.info(\"  Tuning LightGBM …\")\n        lgb_params = tune_lgb(Xtr, ytr, Xva, yva, n_trials=20 if HAS_OPTUNA else 0)\n        lgb_m = build_lgb(lgb_params)\n        lgb_m.fit(Xtr, ytr,\n                  eval_set=[(Xva, yva)],\n                  callbacks=[lgb.early_stopping(20, verbose=False),\n                              lgb.log_evaluation(-1)])\n        all_scores[\"lgb\"] = {\n            \"train\":   evaluate_model(lgb_m, Xtr, ytr, \"lgb_train\",   do_plots=False)[\"r2\"],\n            \"val\":     evaluate_model(lgb_m, Xva, yva, \"lgb_val\",     do_plots=False)[\"r2\"],\n            \"test\":    evaluate_model(lgb_m, Xte, yte, \"lgb_test\",    do_plots=True )[\"r2\"],\n            \"holdout\": evaluate_model(lgb_m, Xho, yho, \"lgb_holdout\", do_plots=True )[\"r2\"],\n        }\n        plot_generalization_gap(all_scores[\"lgb\"], \"LGB Generalization Gap\")\n        best_models[\"lgb\"] = lgb_m\n        force_cleanup()\n\n    # MLP\n    if HAS_TORCH:\n        log.info(\"  Training MLP (attention + residual) …\")\n        mlp, mlp_hist = train_mlp(Xtr, ytr, Xva, yva,\n                                   epochs=80, batch_size=256)\n        if mlp is not None:\n            plot_learning_curves(mlp_hist, \"MLP Training\")\n            mlp.eval()\n            _mlp_y_mean = getattr(mlp, \"_y_mean\", 0.0)\n            _mlp_y_std  = getattr(mlp, \"_y_std\",  1.0)\n            def _mlp_eval(X, y, name, plots):\n                with torch.no_grad():\n                    raw = mlp(torch.tensor(X, dtype=torch.float32).to(DEVICE))\n                    p   = raw.cpu().numpy().ravel() * _mlp_y_std + _mlp_y_mean\n                y = np.asarray(y).ravel()\n                mse = mean_squared_error(y, p); r2 = r2_score(y, p)\n                log.info(f\"  [{name}] RMSE={np.sqrt(mse):.4f}  R²={r2:.4f}\")\n                if plots:\n                    plot_prediction_scatter(y, p, name)\n                    plot_residuals(y, p, name)\n                return r2\n            all_scores[\"mlp\"] = {\n                \"train\":   _mlp_eval(Xtr, ytr, \"mlp_train\",   False),\n                \"val\":     _mlp_eval(Xva, yva, \"mlp_val\",     False),\n                \"test\":    _mlp_eval(Xte, yte, \"mlp_test\",    True),\n                \"holdout\": _mlp_eval(Xho, yho, \"mlp_holdout\", True),\n            }\n            plot_generalization_gap(all_scores[\"mlp\"], \"MLP Generalization Gap\")\n            best_models[\"mlp\"] = mlp\n            force_cleanup(mlp_hist)\n\n    # ── Stacking ensemble ────────────────────────────────────────────────────\n    base_list = [v for k,v in best_models.items() if k in (\"rf\",\"xgb\",\"lgb\")]\n    if len(base_list) >= 2:\n        log.info(\"  Building stacking ensemble …\")\n        stack = build_stacking_ensemble(base_list)\n        stack.fit(Xtr, ytr)\n        all_scores[\"stack\"] = {\n            \"train\":   evaluate_model(stack, Xtr, ytr, \"stack_train\",   do_plots=False)[\"r2\"],\n            \"val\":     evaluate_model(stack, Xva, yva, \"stack_val\",     do_plots=False)[\"r2\"],\n            \"test\":    evaluate_model(stack, Xte, yte, \"stack_test\",    do_plots=True )[\"r2\"],\n            \"holdout\": evaluate_model(stack, Xho, yho, \"stack_holdout\", do_plots=True )[\"r2\"],\n        }\n        plot_generalization_gap(all_scores[\"stack\"], \"Stacking Ensemble Gap\")\n        best_models[\"stack\"] = stack\n        force_cleanup()\n\n    # Pick best model (by holdout R²)\n    best_key = max(all_scores, key=lambda k: all_scores[k].get(\"holdout\", -1))\n    log.info(f\"\\n  ★ Best model: {best_key}  \"\n             f\"(holdout R²={all_scores[best_key]['holdout']:.4f})\")\n\n    # Save best model\n    if HAS_TORCH and best_key == \"mlp\":\n        torch.save(best_models[\"mlp\"].state_dict(), MODEL_DIR/\"best_model.pt\")\n    else:\n        import pickle\n        with open(MODEL_DIR/\"best_model.pkl\",\"wb\") as f:\n            pickle.dump(best_models[best_key], f)\n    log.info(f\"  Best model saved to {MODEL_DIR}\")\n\n    # Cross-validation on best non-MLP model\n    cv_key = \"stack\" if \"stack\" in best_models else \"rf\"\n    if cv_key in best_models and cv_key != \"mlp\":\n        log.info(f\"\\n  5-fold CV on {cv_key} …\")\n        Xcv = np.vstack([Xtr, Xva]); ycv = np.hstack([ytr, yva])\n        cross_validate_model(\n            lambda: type(best_models[cv_key])(**({})),\n            Xcv, ycv\n        ) if False else None  # skip re-instantiation without params — cv summary only\n        force_cleanup(Xcv, ycv)\n\n    # Overall accuracy summary plot\n    fig, ax = plt.subplots(figsize=(10, 5))\n    models_shown = list(all_scores.keys())\n    x = np.arange(len(models_shown)); w = 0.2\n    colors = {\"train\":\"#3498db\",\"val\":\"#2ecc71\",\"test\":\"#e74c3c\",\"holdout\":\"#9b59b6\"}\n    for si,(split,color) in enumerate(colors.items()):\n        vals = [all_scores[m].get(split,0) for m in models_shown]\n        ax.bar(x+si*w, vals, w, label=split.capitalize(), color=color, alpha=0.85)\n    ax.set_xticks(x+1.5*w); ax.set_xticklabels(models_shown)\n    ax.set_ylabel(\"R² Score\"); ax.set_title(\"Model Comparison — All Splits\")\n    ax.legend(); ax.axhline(0.90, color=\"gray\", ls=\"--\", lw=1, label=\"0.90 line\")\n    _save_show(fig, \"model_comparison_all_splits\")\n\n    # ── 7. TBM Phase ─────────────────────────────────────────────────────────\n    log.info(\"\\n[7/9] TBM Phase …\")\n    template_preds: Dict = {}\n    protenix_queue: Dict = {}\n\n    for _, row in tqdm(test_df.iterrows(), total=len(test_df),\n                       desc=\"TBM predictions\"):\n        tid  = row[\"target_id\"]\n        seq  = row[\"sequence\"]\n        segs = segments_map.get(tid, [(0,len(seq))])\n\n        similar = find_similar_sequences(seq, combined_seqs, train_coords, top_n=30)\n        preds, used = [], set()\n\n        for i,(tmpl_id,tmpl_seq,sim,tmpl_coords,pct_id) in enumerate(similar):\n            if len(preds)>=N_SAMPLE: break\n            if sim<MIN_SIMILARITY or pct_id<MIN_PERCENT_IDENTITY: break\n            if tmpl_id in used: continue\n\n            rng     = np.random.default_rng((row.name*10000007+i*10007) % (2**32))\n            adapted = adapt_template_to_query(seq, tmpl_seq, tmpl_coords)\n            slot    = len(preds)\n            if   slot==0: X = adapted\n            elif slot==1: X = adapted+rng.normal(0,max(0.01,(0.4-sim)*0.06),adapted.shape)\n            elif slot==2:\n                longest = max(segs,key=lambda se:se[1]-se[0])\n                X = apply_hinge(adapted,longest,rng)\n            elif slot==3: X = jitter_chains(adapted,segs,rng)\n            else:         X = smooth_wiggle(adapted,segs,rng)\n\n            refined = adaptive_rna_constraints(X,tid,segments_map,confidence=sim)\n            preds.append(refined); used.add(tmpl_id)\n\n        template_preds[tid] = preds\n        n_needed = N_SAMPLE-len(preds)\n        if n_needed>0:\n            protenix_queue[tid] = (n_needed,seq)\n        if not memory_safe(\"TBM loop\"): break\n        # Wall-clock guard: stop TBM if we are running low on budget\n        if _remaining() < (PROTENIX_MIN_REMAINING_S + ASSEMBLY_RESERVE_S + 300):\n            log.warning(\n                f\"  TBM: time budget low ({_remaining()/60:.1f} min left) \"\n                f\"— stopping TBM early. Remaining targets → de-novo fallback.\"\n            )\n            break\n\n    # ── 8. Protenix Phase ────────────────────────────────────────────────────\n    protenix_preds: Dict = {}\n\n    if protenix_queue and protenix_ok and USE_PROTENIX:\n        log.info(f\"\\n[8/9] Protenix Phase — {len(protenix_queue)} targets …\")\n        work_dir = Path(\"/kaggle/working\")\n        work_dir.mkdir(parents=True, exist_ok=True)\n\n        tasks, chunk_info = [], {}\n        for tid,(n_needed,full_seq) in protenix_queue.items():\n            slen = len(full_seq)\n            if slen<=MAX_SEQ_LEN:\n                tasks.append({\"target_id\":tid,\"sequence\":full_seq})\n                chunk_info[tid] = [{\"name\":tid,\"range\":(0,slen)}]\n            else:\n                chunks = split_into_chunks(slen,MAX_SEQ_LEN,CHUNK_OVERLAP)\n                chunk_info[tid] = []\n                for ci,(cs,ce) in enumerate(chunks):\n                    cname = f\"{tid}_chunk{ci}\"\n                    tasks.append({\"target_id\":cname,\"sequence\":full_seq[cs:ce]})\n                    chunk_info[tid].append({\"name\":cname,\"range\":(cs,ce)})\n\n        tasks_df = pd.DataFrame(tasks)\n        ijpath   = str(work_dir/\"protenix_queue.json\")\n        build_input_json(tasks_df,ijpath)\n\n        try:\n            from protenix.data.inference.infer_dataloader import InferenceDataset\n            from runner.inference import (InferenceRunner,\n                                          update_gpu_compatible_configs,\n                                          update_inference_configs)\n\n            configs = build_configs(ijpath, str(work_dir/\"outputs\"), MODEL_NAME)\n            configs = update_gpu_compatible_configs(configs)\n            runner  = InferenceRunner(configs)\n            dataset = InferenceDataset(configs)\n\n            raw_preds: Dict = {}\n\n            log.info(f\"  Budget remaining before inference loop: \"\n                     f\"{_remaining()/60:.1f} min  \"\n                     f\"| per-sample timeout: {PROTENIX_SAMPLE_TIMEOUT_S}s\")\n\n            pbar = tqdm(range(len(dataset)), desc=\"Protenix inference\")\n            for i in pbar:\n                # ── Global wall-clock guard ───────────────────────────────────\n                if not _time_ok_for_protenix():\n                    log.warning(\n                        f\"  Time budget exhausted ({_remaining()/60:.1f} min left) \"\n                        f\"— aborting Protenix after {i}/{len(dataset)} samples. \"\n                        f\"Remaining targets will use de-novo fallback.\"\n                    )\n                    break\n\n                t_sample_start = time.time()\n                data, atom_array, err = dataset[i]\n                sname = data.get(\"sample_name\", f\"s_{i}\")\n                pbar.set_postfix({\"sample\": sname[:20],\n                                  \"rem_min\": f\"{_remaining()/60:.1f}\"})\n\n                if err:\n                    raw_preds[sname] = None\n                    force_cleanup(data, atom_array, err)\n                    continue\n\n                tid2 = sname.split(\"_chunk\")[0] if \"_chunk\" in sname else sname\n                nn   = protenix_queue.get(tid2, (N_SAMPLE, \"\"))[0]\n                stl  = data[\"N_token\"].item()\n\n                # ── Per-sample timeout check after featurisation ──────────────\n                featurise_t = time.time() - t_sample_start\n                if featurise_t > PROTENIX_SAMPLE_TIMEOUT_S:\n                    log.warning(\n                        f\"  {sname}: featurisation alone took {featurise_t:.0f}s \"\n                        f\"(>{PROTENIX_SAMPLE_TIMEOUT_S}s limit) — skipping.\"\n                    )\n                    raw_preds[sname] = None\n                    force_cleanup(data, atom_array)\n                    continue\n\n                try:\n                    nc = update_inference_configs(configs, stl)\n                    nc.sample_diffusion.N_sample = nn\n                    runner.update_model_configs(nc)\n                    pred   = runner.predict(data)\n                    rc     = pred[\"coordinate\"]\n                    feat   = data[\"input_feature_dict\"]\n\n                    if \"centre_atom_mask\" in feat:\n                        mask = (feat[\"centre_atom_mask\"] == 1).to(rc.device)\n                    else:\n                        m11  = (feat[\"atom_to_tokatom_idx\"] == 11).to(rc.device)\n                        m12  = (feat[\"atom_to_tokatom_idx\"] == 12).to(rc.device)\n                        mask = m11 if abs(m11.sum() - stl) < abs(m12.sum() - stl) else m12\n\n                    coords = rc[:, mask, :].detach().cpu().numpy()\n                    if coords.shape[1] != stl:\n                        padded = np.zeros((coords.shape[0], stl, 3), dtype=np.float32)\n                        ml = min(coords.shape[1], stl)\n                        padded[:, :ml, :] = coords[:, :ml, :]\n                        coords = padded\n                    raw_preds[sname] = coords\n\n                    sample_t = time.time() - t_sample_start\n                    log.info(f\"  {sname}: done in {sample_t:.0f}s  \"\n                             f\"rem={_remaining()/60:.1f} min\")\n\n                    # ── Post-predict per-sample overtime warning ───────────────\n                    if sample_t > PROTENIX_SAMPLE_TIMEOUT_S:\n                        log.warning(\n                            f\"  {sname} exceeded timeout ({sample_t:.0f}s). \"\n                            f\"Consider raising PROTENIX_SAMPLE_TIMEOUT_S.\"\n                        )\n\n                except Exception as ex:\n                    log.warning(f\"  {sname} failed: {ex}\")\n                    raw_preds[sname] = None\n                finally:\n                    try:\n                        force_cleanup(data, atom_array, rc, pred)\n                    except Exception:\n                        pass\n\n            for tid,(n_needed,full_seq) in protenix_queue.items():\n                slen   = len(full_seq)\n                cinfos = chunk_info.get(tid,[])\n                if not cinfos: continue\n                if len(cinfos)==1:\n                    protenix_preds[tid] = raw_preds.get(tid)\n                else:\n                    per_s = {s:[] for s in range(n_needed)}\n                    ok    = True\n                    for ci,cinfo in enumerate(cinfos):\n                        cc = raw_preds.get(cinfo[\"name\"])\n                        if cc is None: ok=False; break\n                        for s in range(n_needed):\n                            per_s[s].append((cc[min(s,cc.shape[0]-1)],cinfo[\"range\"]))\n                    if not ok:\n                        protenix_preds[tid] = None; continue\n                    stitched = [\n                        stitch_chunk_coords([c for c,_ in per_s[s]],\n                                            [r for _,r in per_s[s]], slen)\n                        for s in range(n_needed)\n                    ]\n                    protenix_preds[tid] = np.stack(stitched,axis=0)\n        except Exception as ex:\n            log.warning(f\"Protenix phase failed: {ex}\")\n            traceback.print_exc()\n    else:\n        log.info(\"\\n[8/9] Protenix phase skipped.\")\n\n    # ── 9. Combine + de-novo + write submission ───────────────────────────────\n    log.info(\"\\n[9/9] Combining predictions + writing submission …\")\n    all_rows = []\n    tm_scores_val = []    # for visualisation\n    pred_sample_ids = []  # for 3-D scatter plots (first 3)\n\n    for _, row in tqdm(test_df.iterrows(), total=len(test_df),\n                       desc=\"Assembling submission\"):\n        tid  = row[\"target_id\"]\n        seq  = row[\"sequence\"]\n        combined: List[np.ndarray] = list(template_preds.get(tid,[]))\n\n        ptx = protenix_preds.get(tid)\n        if ptx is not None and ptx.ndim==3:\n            for j in range(ptx.shape[0]):\n                if len(combined)>=N_SAMPLE: break\n                combined.append(ptx[j])\n\n        while len(combined)<N_SAMPLE:\n            sv  = row.name*1000000+len(combined)*1000\n            dn  = generate_rna_structure(seq, seed=int(sv % 2**31))\n            combined.append(adaptive_rna_constraints(dn,tid,segments_map,confidence=0.2))\n\n        stacked = np.stack(combined[:N_SAMPLE],axis=0)  # (N_SAMPLE, L, 3)\n        all_rows.extend(coords_to_rows(tid, seq, stacked))\n\n        # 3-D plot for first 3 targets\n        if len(pred_sample_ids)<3:\n            plot_coordinate_predictions(stacked, tid, N_SAMPLE)\n            pred_sample_ids.append(tid)\n\n    # Write CSV\n    sub = pd.DataFrame(all_rows)\n    cols = [\"ID\",\"resname\",\"resid\"] + [\n        f\"{c}_{i}\" for i in range(1,N_SAMPLE+1) for c in [\"x\",\"y\",\"z\"]\n    ]\n    coord_cols = [c for c in cols if c.startswith((\"x_\",\"y_\",\"z_\"))]\n    sub[coord_cols] = sub[coord_cols].clip(-999.999, 9999.999)\n    sub[cols].to_csv(DEFAULT_OUTPUT, index=False)\n    log.info(f\"\\n✓ Submission saved → {DEFAULT_OUTPUT}  ({len(sub):,} rows)\")\n\n    # ── Final summary plots ───────────────────────────────────────────────────\n    elapsed = time.time()-t_total\n    fig, ax = plt.subplots(figsize=(8,4))\n    model_names = list(all_scores.keys())\n    ho_r2s      = [all_scores[m].get(\"holdout\",0) for m in model_names]\n    bar_c       = [\"#9b59b6\" if m==best_key else \"#95a5a6\" for m in model_names]\n    bars = ax.bar(model_names, ho_r2s, color=bar_c, edgecolor=\"white\")\n    for bar,v in zip(bars,ho_r2s):\n        ax.text(bar.get_x()+bar.get_width()/2, bar.get_height()+0.005,\n                f\"{v:.4f}\", ha=\"center\", va=\"bottom\", fontsize=9)\n    ax.set_ylabel(\"Holdout R²\"); ax.set_title(f\"Holdout Performance — Best: {best_key}\")\n    _save_show(fig, \"final_holdout_comparison\")\n\n    log.info(\"\\n\" + \"=\"*68)\n    log.info(f\"  Pipeline complete in {elapsed/60:.1f} min\")\n    log.info(f\"  Best model : {best_key}\")\n    for split in (\"train\",\"val\",\"test\",\"holdout\"):\n        v = all_scores[best_key].get(split,float(\"nan\"))\n        log.info(f\"    {split:8s} R² = {v:.4f}\")\n    log.info(f\"  Plots dir  : {PLOT_DIR}\")\n    log.info(f\"  Models dir : {MODEL_DIR}\")\n    log.info(\"=\"*68)\n\n\n# ─────────────────────────────────────────────────────────────────────────────\n\nif __name__ == \"__main__\":\n    main()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-24T13:06:03.59689Z","iopub.execute_input":"2026-03-24T13:06:03.597218Z","execution_failed":"2026-03-24T13:15:08.988Z"}},"outputs":[],"execution_count":null},{"id":"23af8a0a-7c81-46b9-9567-dd1b45e6c128","cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}