{"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":"gpu","dataSources":[{"databundleVersionId":15231210,"sourceId":118765,"sourceType":"competition"},{"databundleVersionId":15819611,"datasetId":9567922,"sourceId":14949595,"sourceType":"datasetVersion"},{"databundleVersionId":15440074,"datasetId":9328538,"sourceId":14604295,"sourceType":"datasetVersion"},{"databundleVersionId":15833819,"datasetId":9577079,"sourceId":14962460,"sourceType":"datasetVersion"},{"databundleVersionId":16259964,"datasetId":9582898,"sourceId":15350286,"sourceType":"datasetVersion"},{"databundleVersionId":16278719,"datasetId":9648466,"sourceId":15367418,"sourceType":"datasetVersion"},{"databundleVersionId":11442412,"datasetId":6888333,"sourceId":11056335,"sourceType":"datasetVersion"},{"databundleVersionId":11442459,"datasetId":6888367,"sourceId":11056379,"sourceType":"datasetVersion"},{"databundleVersionId":16283972,"datasetId":9827351,"sourceId":15372130,"sourceType":"datasetVersion"},{"databundleVersionId":15980951,"datasetId":9665334,"sourceId":15096066,"sourceType":"datasetVersion"},{"databundleVersionId":15846521,"datasetId":9584272,"sourceId":14974056,"sourceType":"datasetVersion"},{"databundleVersionId":15842735,"datasetId":9567965,"sourceId":14970606,"sourceType":"datasetVersion"},{"databundleVersionId":15736806,"datasetId":9502242,"sourceId":14874339,"sourceType":"datasetVersion"},{"databundleVersionId":15996885,"datasetId":9674779,"sourceId":15110328,"sourceType":"datasetVersion"},{"databundleVersionId":15833858,"datasetId":9577097,"sourceId":14962495,"sourceType":"datasetVersion"},{"databundleVersionId":12505876,"datasetId":7541592,"sourceId":11990166,"sourceType":"datasetVersion"},{"databundleVersionId":15948658,"datasetId":9568372,"sourceId":15066946,"sourceType":"datasetVersion"},{"databundleVersionId":15347344,"datasetId":9271415,"sourceId":14519720,"sourceType":"datasetVersion"},{"databundleVersionId":15980738,"datasetId":9660182,"sourceId":15095867,"sourceType":"datasetVersion"},{"databundleVersionId":16190954,"datasetId":9724076,"sourceId":15288472,"sourceType":"datasetVersion"}],"dockerImageVersionId":31286,"isGpuEnabled":true,"isInternetEnabled":false,"language":"python","sourceType":"notebook"},"papermill":{"default_parameters":{},"duration":12044.146287,"end_time":"2026-03-27T21:30:44.032898","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2026-03-27T18:09:59.886611","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# # This Python 3 environment comes with many helpful analytics libraries installed\n# # It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# # For example, here's several helpful packages to load\n\n# import numpy as np # linear algebra\n# import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# # Input data files are available in the read-only \"../input/\" directory\n# # For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n# for filename in filenames:\n# print(os.path.join(dirname, filename))\n\n# # You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# # You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.execute_input":"2026-03-27T18:10:03.188373Z","iopub.status.busy":"2026-03-27T18:10:03.187904Z","iopub.status.idle":"2026-03-27T18:10:03.1987Z","shell.execute_reply":"2026-03-27T18:10:03.197702Z"},"papermill":{"duration":0.039664,"end_time":"2026-03-27T18:10:03.200787","exception":false,"start_time":"2026-03-27T18:10:03.161123","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import os\n# import sys\n# import shutil\n\n# # 1. Paths (Ensure these match your actual dataset structure)\n# OFFLINE_DATASET_PATH = \"/kaggle/input/datasets/khushvaghela/gcosrna/python_packages/kaggle/working/python_packages\"\n# LIB_PATH = \"/kaggle/working/custom_lib\"\n\n# # 2. CRITICAL: Wipe the existing corrupted library\n# if os.path.exists(LIB_PATH):\n# print(\"๐Ÿงน Clearing corrupted custom library...\")\n# shutil.rmtree(LIB_PATH)\n# os.makedirs(LIB_PATH, exist_ok=True)\n\n# # 3. Fresh Install (Added --upgrade for extra safety)\n# print(\"โณ Installing clean NumPy 1.26 stack...\")\n# !pip install --no-index --find-links={OFFLINE_DATASET_PATH} \\\n# --target={LIB_PATH} \\\n# --upgrade \\\n# --no-deps \\\n# \"scipy==1.11.4\" \"scikit-learn==1.3.2\"\n\n# # 4. Re-install dependencies (if they are in your folder)\n# # This ensures things like 'six' and 'pytz' are available in the custom path\n# !pip install --no-index --find-links={OFFLINE_DATASET_PATH} \\\n# --target={LIB_PATH} \\\n# \"python-dateutil\" \"pytz\" \"tzdata\" \"joblib\" \"threadpoolctl\"\n\n# # 5. Hijack the Path\n# sys.path.insert(0, LIB_PATH)\n\n# # 6. The Nuclear Clear (Clear memory of any previous imports)\n# for module in list(sys.modules.keys()):\n# if any(m in module for m in ['numpy', 'pandas', 'scipy', 'sklearn']):\n# del sys.modules[module]\n\n# print(\"โœ… Environment Purified.\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:10:03.238413Z","iopub.status.busy":"2026-03-27T18:10:03.23811Z","iopub.status.idle":"2026-03-27T18:10:03.242907Z","shell.execute_reply":"2026-03-27T18:10:03.24207Z"},"papermill":{"duration":0.023642,"end_time":"2026-03-27T18:10:03.244691","exception":false,"start_time":"2026-03-27T18:10:03.221049","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## MxFold2","metadata":{"papermill":{"duration":0.017932,"end_time":"2026-03-27T18:10:03.278132","exception":false,"start_time":"2026-03-27T18:10:03.2602","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# 1. Define paths\nsource = \"/kaggle/input/datasets/khushvaghela/mxfold2/mxfold2-master/mxfold2-master\"\ntemp_copy = \"/kaggle/working/mxfold2-master\"\nfinal_name = \"/kaggle/working/mxfold2\"\n\n# 2. Copy the whole folder to working first\n!cp -r {source} {temp_copy}\n\n# 3. Rename the folder to exactly what you want\n!mv {temp_copy} {final_name}\n\n# Verify the structure (should show files, not another mxfold2-master folder)\n!ls {final_name}","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:10:03.339043Z","iopub.status.busy":"2026-03-27T18:10:03.337804Z","iopub.status.idle":"2026-03-27T18:10:03.951587Z","shell.execute_reply":"2026-03-27T18:10:03.950139Z"},"papermill":{"duration":0.648759,"end_time":"2026-03-27T18:10:03.953948","exception":false,"start_time":"2026-03-27T18:10:03.305189","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Install poetry-core and poetry from your dataset\n!pip install --no-index --find-links=/kaggle/input/datasets/khushvaghela/build-tools/poetry poetry-core poetry","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:10:03.985857Z","iopub.status.busy":"2026-03-27T18:10:03.985308Z","iopub.status.idle":"2026-03-27T18:10:12.181765Z","shell.execute_reply":"2026-03-27T18:10:12.180604Z"},"papermill":{"duration":8.215278,"end_time":"2026-03-27T18:10:12.183903","exception":false,"start_time":"2026-03-27T18:10:03.968625","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Install using the existing environment's poetry\n!pip install {final_name} --no-index --no-deps --no-build-isolation","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:10:12.217345Z","iopub.status.busy":"2026-03-27T18:10:12.216792Z","iopub.status.idle":"2026-03-27T18:11:16.780288Z","shell.execute_reply":"2026-03-27T18:11:16.779334Z"},"papermill":{"duration":64.582841,"end_time":"2026-03-27T18:11:16.782581","exception":false,"start_time":"2026-03-27T18:10:12.19974","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!mxfold2 --help","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:11:16.816626Z","iopub.status.busy":"2026-03-27T18:11:16.815773Z","iopub.status.idle":"2026-03-27T18:11:56.280767Z","shell.execute_reply":"2026-03-27T18:11:56.279451Z"},"papermill":{"duration":39.484214,"end_time":"2026-03-27T18:11:56.28295","exception":false,"start_time":"2026-03-27T18:11:16.798736","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!mxfold2 predict @/kaggle/input/datasets/khushvaghela/mxfold2/models-0.1.0/models/TrainSetA.conf /kaggle/input/datasets/khushvaghela/mxfold2/test.fasta 2>/dev/null","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:11:56.318521Z","iopub.status.busy":"2026-03-27T18:11:56.317503Z","iopub.status.idle":"2026-03-27T18:12:08.947693Z","shell.execute_reply":"2026-03-27T18:12:08.946849Z"},"papermill":{"duration":12.650549,"end_time":"2026-03-27T18:12:08.94989","exception":false,"start_time":"2026-03-27T18:11:56.299341","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Setting up IPKnot","metadata":{"papermill":{"duration":0.01624,"end_time":"2026-03-27T18:12:08.982797","exception":false,"start_time":"2026-03-27T18:12:08.966557","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"## Offline","metadata":{"papermill":{"duration":0.016708,"end_time":"2026-03-27T18:12:09.015646","exception":false,"start_time":"2026-03-27T18:12:08.998938","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import os\n\n# --- 1. INSTALL SYSTEM DEPENDENCIES (.deb files) ---\n# Install everything from your uploaded dataset\n!dpkg -i /kaggle/input/datasets/khushvaghela/ipknot/system_libs/system_libs/*.deb","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:09.051199Z","iopub.status.busy":"2026-03-27T18:12:09.050466Z","iopub.status.idle":"2026-03-27T18:12:18.690909Z","shell.execute_reply":"2026-03-27T18:12:18.689769Z"},"papermill":{"duration":9.660821,"end_time":"2026-03-27T18:12:18.693231","exception":false,"start_time":"2026-03-27T18:12:09.03241","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- 2. SETUP DIRECTORIES ---\n# Copy the repo from your dataset to the writable working directory\n!cp -r /kaggle/input/datasets/khushvaghela/ipknot/ipknot /kaggle/working/ipknot","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:18.734024Z","iopub.status.busy":"2026-03-27T18:12:18.733355Z","iopub.status.idle":"2026-03-27T18:12:19.775549Z","shell.execute_reply":"2026-03-27T18:12:19.77439Z"},"papermill":{"duration":1.064918,"end_time":"2026-03-27T18:12:19.777855","exception":false,"start_time":"2026-03-27T18:12:18.712937","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- 3. SETUP VIENNARNA LINK ---\n# Use a wildcard to ensure it finds the correct .so file for Python 3.12\n!ln -sf /kaggle/input/datasets/khushvaghela/ipknot/_RNA.cpython-312-x86_64-linux-gnu.so /kaggle/working/libRNA.so","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:19.82045Z","iopub.status.busy":"2026-03-27T18:12:19.819491Z","iopub.status.idle":"2026-03-27T18:12:19.942337Z","shell.execute_reply":"2026-03-27T18:12:19.941097Z"},"papermill":{"duration":0.146403,"end_time":"2026-03-27T18:12:19.944694","exception":false,"start_time":"2026-03-27T18:12:19.798291","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Ensure you have the .whl file in your 'wheels' dataset\n!pip install /kaggle/input/datasets/khushvaghela/ipknot/viennarna-2.7.2-cp312-cp312-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl --no-index --no-deps","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:19.987334Z","iopub.status.busy":"2026-03-27T18:12:19.986257Z","iopub.status.idle":"2026-03-27T18:12:22.192908Z","shell.execute_reply":"2026-03-27T18:12:22.191669Z"},"papermill":{"duration":2.231063,"end_time":"2026-03-27T18:12:22.195631","exception":false,"start_time":"2026-03-27T18:12:19.964568","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!ls /usr/local/lib/python3.12/dist-packages/RNA/","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:22.23627Z","iopub.status.busy":"2026-03-27T18:12:22.23593Z","iopub.status.idle":"2026-03-27T18:12:22.354485Z","shell.execute_reply":"2026-03-27T18:12:22.353389Z"},"papermill":{"duration":0.141891,"end_time":"2026-03-27T18:12:22.356712","exception":false,"start_time":"2026-03-27T18:12:22.214821","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\n\n# 1. Define the verified path\nPYTHON_SO = \"/usr/local/lib/python3.12/dist-packages/RNA/_RNA.cpython-312-x86_64-linux-gnu.so\"\nLIB_DEST = \"/kaggle/working/libRNA.so\"\n\n# 2. Create the symbolic link\n# This allows the compiler to find the library under a standard name\n!ln -sf {PYTHON_SO} {LIB_DEST}","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:22.399039Z","iopub.status.busy":"2026-03-27T18:12:22.398162Z","iopub.status.idle":"2026-03-27T18:12:22.524339Z","shell.execute_reply":"2026-03-27T18:12:22.523339Z"},"papermill":{"duration":0.149399,"end_time":"2026-03-27T18:12:22.526398","exception":false,"start_time":"2026-03-27T18:12:22.376999","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!dpkg-deb -x /kaggle/input/datasets/khushvaghela/ipknot/viennarna-dev_2.7.2-1_amd64.deb ./extracted_headers","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:22.570341Z","iopub.status.busy":"2026-03-27T18:12:22.569801Z","iopub.status.idle":"2026-03-27T18:12:23.675653Z","shell.execute_reply":"2026-03-27T18:12:23.674689Z"},"papermill":{"duration":1.128881,"end_time":"2026-03-27T18:12:23.677786","exception":false,"start_time":"2026-03-27T18:12:22.548905","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Run this after extraction but before zipping/uploading\n!mv ./extracted_headers/usr/include/ViennaRNA ./extracted_headers/usr/include/RNA","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:23.720405Z","iopub.status.busy":"2026-03-27T18:12:23.719873Z","iopub.status.idle":"2026-03-27T18:12:23.839144Z","shell.execute_reply":"2026-03-27T18:12:23.838118Z"},"papermill":{"duration":0.143802,"end_time":"2026-03-27T18:12:23.841406","exception":false,"start_time":"2026-03-27T18:12:23.697604","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 1. Define the verified path\nPYTHON_SO = \"/usr/local/lib/python3.12/dist-packages/RNA/_RNA.cpython-312-x86_64-linux-gnu.so\"\nLIB_DEST = \"/kaggle/working/libRNA.so\"\n\n# 2. Create the symbolic link\n# This allows the compiler to find the library under a standard name\n!ln -sf {PYTHON_SO} {LIB_DEST}","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:23.883623Z","iopub.status.busy":"2026-03-27T18:12:23.882572Z","iopub.status.idle":"2026-03-27T18:12:24.00167Z","shell.execute_reply":"2026-03-27T18:12:24.000427Z"},"papermill":{"duration":0.142594,"end_time":"2026-03-27T18:12:24.004155","exception":false,"start_time":"2026-03-27T18:12:23.861561","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Copy the repo from your dataset to the writable working directory\n!cp -r /kaggle/input/datasets/khushvaghela/ipknot/ipknot /kaggle/working/ipknot","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:24.046052Z","iopub.status.busy":"2026-03-27T18:12:24.045229Z","iopub.status.idle":"2026-03-27T18:12:24.57087Z","shell.execute_reply":"2026-03-27T18:12:24.569573Z"},"papermill":{"duration":0.549148,"end_time":"2026-03-27T18:12:24.573281","exception":false,"start_time":"2026-03-27T18:12:24.024133","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!ls /kaggle/working","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:24.614375Z","iopub.status.busy":"2026-03-27T18:12:24.614034Z","iopub.status.idle":"2026-03-27T18:12:24.733931Z","shell.execute_reply":"2026-03-27T18:12:24.732673Z"},"papermill":{"duration":0.142814,"end_time":"2026-03-27T18:12:24.736167","exception":false,"start_time":"2026-03-27T18:12:24.593353","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import sys\nimport os\n\n# Define all absolute paths\npy_ver = f\"{sys.version_info.major}.{sys.version_info.minor}\"\npy_lib = f\"/usr/lib/x86_64-linux-gnu/libpython{py_ver}.so\"\npy_inc = f\"/usr/include/python{py_ver}\"\nvienna_inc = \"/kaggle/working/extracted_headers/usr/include\"\nvienna_lib = \"/kaggle/working/libRNA.so\"\nglpk_inc = \"/usr/include\"\nglpk_lib = \"/usr/lib/x86_64-linux-gnu/libglpk.so\"\n\ncmake_path = \"/kaggle/working/ipknot/CMakeLists.txt\"\n\nfinal_cmake = f\"\"\"\ncmake_minimum_required(VERSION 3.10)\nproject(ipknot CXX C)\n\nset(CMAKE_CXX_STANDARD 17)\nset(CMAKE_CXX_STANDARD_REQUIRED ON)\nset(CMAKE_BUILD_TYPE Release)\n\n# --- THE SOLVER FLAGS (Enabling both common variants) ---\nadd_definitions(-DPACKAGE_VERSION=\"2.0.0-offline\")\nadd_definitions(-DENABLE_GLPK=1)\nadd_definitions(-DHAVE_GLPK=1)\nadd_definitions(-DHAVE_LIBGLPK=1)\n\n# --- Manual Dependency Setup ---\nfind_package(OpenMP REQUIRED)\ninclude_directories(\"{vienna_inc}\" \"{glpk_inc}\" \"{py_inc}\" \"src\" \"src/contrafold\" \"src/linearpartition\")\n\n# --- Linker Flags ---\nset(CMAKE_EXE_LINKER_FLAGS \"${{CMAKE_EXE_LINKER_FLAGS}} -Wl,--allow-multiple-definition -lgomp\")\n\n# --- Manual Linkage ---\nlink_libraries(\"{vienna_lib}\" \"{glpk_lib}\" \"{py_lib}\" OpenMP::OpenMP_CXX gsl gslcblas)\n\n# --- Source Selection ---\nset(SOURCES\n src/aln.cpp src/boltzmann_param.c src/bpseq.cpp src/fa.cpp src/fold.cpp\n src/ip.cpp src/ipknot.cpp src/nupack/nupack.cpp\n src/contrafold/InferenceEngine.cpp src/contrafold/ParameterManager.cpp\n src/contrafold/SStruct.cpp src/contrafold/Utilities.cpp\n src/linearpartition/LinearPartition.cpp\n)\n\nadd_executable(ipknot ${{SOURCES}})\n\"\"\"\n\nwith open(cmake_path, \"w\") as f:\n f.write(final_cmake)\n\nprint(\"CMakeLists.txt updated with exhaustive GLPK definitions.\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:24.77768Z","iopub.status.busy":"2026-03-27T18:12:24.777095Z","iopub.status.idle":"2026-03-27T18:12:24.785413Z","shell.execute_reply":"2026-03-27T18:12:24.784664Z"},"papermill":{"duration":0.031164,"end_time":"2026-03-27T18:12:24.78711","exception":false,"start_time":"2026-03-27T18:12:24.755946","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\n\nip_cpp_path = \"/kaggle/working/ipknot/src/ip.cpp\"\n\n# A fully corrected, consolidated version of your ip.cpp\nfinal_src = r\"\"\"\n#define WITH_GLPK 1\n#include \"ip.h\"\n#include \n#include \n#include \n#include \n#include \n#include \n\n#ifdef WITH_GLPK\nclass IPimpl\n{\npublic:\n IPimpl(IP::DirType dir, int n_th)\n : ip_(NULL), ia_(1), ja_(1), ar_(1)\n {\n ip_ = glp_create_prob();\n switch (dir)\n {\n case IP::MIN: glp_set_obj_dir(ip_, GLP_MIN); break;\n case IP::MAX: glp_set_obj_dir(ip_, GLP_MAX); break;\n }\n }\n\n ~IPimpl()\n {\n glp_delete_prob(ip_);\n }\n\n int make_variable(double coef)\n {\n int col = glp_add_cols(ip_, 1);\n glp_set_col_bnds(ip_, col, GLP_DB, 0, 1);\n glp_set_col_kind(ip_, col, GLP_BV);\n glp_set_obj_coef(ip_, col, coef);\n return col;\n }\n\n int make_variable(double coef, int lo, int hi)\n {\n int col = glp_add_cols(ip_, 1);\n glp_set_col_bnds(ip_, col, GLP_DB, lo, hi);\n glp_set_col_kind(ip_, col, GLP_IV);\n glp_set_obj_coef(ip_, col, coef);\n return col;\n }\n\n int make_constraint(IP::BoundType bnd, double l, double u)\n {\n int row = glp_add_rows(ip_, 1);\n switch (bnd)\n {\n case IP::FR: glp_set_row_bnds(ip_, row, GLP_FR, l, u); break;\n case IP::LO: glp_set_row_bnds(ip_, row, GLP_LO, l, u); break;\n case IP::UP: glp_set_row_bnds(ip_, row, GLP_UP, l, u); break;\n case IP::DB: glp_set_row_bnds(ip_, row, GLP_DB, l, u); break;\n case IP::FX: glp_set_row_bnds(ip_, row, GLP_FX, l, u); break;\n }\n return row;\n }\n\n void add_constraint(int row, int col, double val)\n {\n assert(row>=0);\n ia_.push_back(row);\n assert(col>=0);\n ja_.push_back(col);\n ar_.push_back(val);\n }\n\n void update() {}\n\n double solve()\n {\n glp_smcp smcp;\n glp_iocp iocp;\n glp_init_smcp(&smcp); smcp.msg_lev = GLP_MSG_ERR;\n glp_init_iocp(&iocp); iocp.msg_lev = GLP_MSG_ERR;\n glp_load_matrix(ip_, ia_.size()-1, &ia_[0], &ja_[0], &ar_[0]);\n glp_simplex(ip_, &smcp);\n glp_intopt(ip_, &iocp);\n return glp_mip_obj_val(ip_);\n }\n\n double get_value(int col) const\n {\n return glp_mip_col_val(ip_, col);\n }\n\nprivate:\n glp_prob *ip_;\n std::vector ia_;\n std::vector ja_;\n std::vector ar_;\n};\n\nIP::IP(DirType dir, int n_th) : impl_(new IPimpl(dir, n_th)) {}\nIP::~IP() { delete impl_; }\nint IP::make_variable(double coef) { return impl_->make_variable(coef); }\nint IP::make_variable(double coef, int lo, int hi) { return impl_->make_variable(coef, lo, hi); }\nint IP::make_constraint(BoundType bnd, double l, double u) { return impl_->make_constraint(bnd, l, u); }\nvoid IP::add_constraint(int row, int col, double val) { impl_->add_constraint(row, col, val); }\nvoid IP::update() { impl_->update(); }\ndouble IP::solve() { return impl_->solve(); }\ndouble IP::get_value(int col) const { return impl_->get_value(col); }\n\n#else\n#error \"WITH_GLPK must be defined for this offline build\"\n#endif\n\"\"\"\n\nwith open(ip_cpp_path, \"w\") as f:\n f.write(final_src)\n\nprint(\"ip.cpp rewritten with clean C++ structure and matched signatures.\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:24.830494Z","iopub.status.busy":"2026-03-27T18:12:24.830176Z","iopub.status.idle":"2026-03-27T18:12:24.838545Z","shell.execute_reply":"2026-03-27T18:12:24.837503Z"},"papermill":{"duration":0.032438,"end_time":"2026-03-27T18:12:24.84036","exception":false,"start_time":"2026-03-27T18:12:24.807922","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%cd /kaggle/working/ipknot\n!rm -rf build && mkdir build","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:24.884454Z","iopub.status.busy":"2026-03-27T18:12:24.883996Z","iopub.status.idle":"2026-03-27T18:12:25.009999Z","shell.execute_reply":"2026-03-27T18:12:25.008729Z"},"papermill":{"duration":0.151805,"end_time":"2026-03-27T18:12:25.012251","exception":false,"start_time":"2026-03-27T18:12:24.860446","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!cmake -S . -B build -DCMAKE_BUILD_TYPE=Release","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:25.05427Z","iopub.status.busy":"2026-03-27T18:12:25.053408Z","iopub.status.idle":"2026-03-27T18:12:26.627865Z","shell.execute_reply":"2026-03-27T18:12:26.626983Z"},"papermill":{"duration":1.597855,"end_time":"2026-03-27T18:12:26.630089","exception":false,"start_time":"2026-03-27T18:12:25.032234","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!cmake --build build --parallel $(nproc)","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:26.674048Z","iopub.status.busy":"2026-03-27T18:12:26.673025Z","iopub.status.idle":"2026-03-27T18:12:47.258541Z","shell.execute_reply":"2026-03-27T18:12:47.257351Z"},"papermill":{"duration":20.60982,"end_time":"2026-03-27T18:12:47.260675","exception":false,"start_time":"2026-03-27T18:12:26.650855","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\n\n# Update the environment variable for this session\n# Replace the path below with the result from the 'find' command above\nhidden_lib_path = \"/usr/local/lib/python3.12/dist-packages/viennarna.libs:\" \n\nos.environ['LD_LIBRARY_PATH'] = f\"{hidden_lib_path}:\" + os.environ.get('LD_LIBRARY_PATH', '')\n\n# Now test the binary\n!/kaggle/working/ipknot/build/ipknot --version","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:47.304352Z","iopub.status.busy":"2026-03-27T18:12:47.303852Z","iopub.status.idle":"2026-03-27T18:12:47.447226Z","shell.execute_reply":"2026-03-27T18:12:47.44601Z"},"papermill":{"duration":0.168136,"end_time":"2026-03-27T18:12:47.44969","exception":false,"start_time":"2026-03-27T18:12:47.281554","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Run the test again\nimport os\n# Ensure the loader finds the hidden ViennaRNA libs\n# Replace with your actual path if 'find' showed something else\nos.environ['LD_LIBRARY_PATH'] = \"/usr/local/lib/python3.12/dist-packages/viennarna.libs:\" + os.environ.get('LD_LIBRARY_PATH', '')\n\n!echo -e \">test\\nGACUGGCGAUUCGGAGGCCCUUAGGCCUUACGGGGCAAAUCUCGA\" > test.fa\n!/kaggle/working/ipknot/build/ipknot test.fa","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:47.499707Z","iopub.status.busy":"2026-03-27T18:12:47.4992Z","iopub.status.idle":"2026-03-27T18:12:47.75152Z","shell.execute_reply":"2026-03-27T18:12:47.750404Z"},"papermill":{"duration":0.280072,"end_time":"2026-03-27T18:12:47.753774","exception":false,"start_time":"2026-03-27T18:12:47.473702","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Test 1: Help command\n!/kaggle/working/ipknot/build/ipknot -h\n\n# Test 2: Dummy prediction\n!echo -e \">test\\nGACUAGGGACUCAGGAGCUUCAGGACGGUACGCGAUCGAGCUUCGAGGACUCUUUUCAGUCGGAU\" > test.fa\n!/kaggle/working/ipknot/build/ipknot test.fa","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:47.798794Z","iopub.status.busy":"2026-03-27T18:12:47.7983Z","iopub.status.idle":"2026-03-27T18:12:48.180766Z","shell.execute_reply":"2026-03-27T18:12:48.179941Z"},"papermill":{"duration":0.407406,"end_time":"2026-03-27T18:12:48.182886","exception":false,"start_time":"2026-03-27T18:12:47.77548","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import subprocess\nimport os\nimport tempfile\nimport re\n\n# --- Configuration Paths ---\nMXFOLD2_CONFIG = \"@/kaggle/input/datasets/khushvaghela/mxfold2/models-0.1.0/models/TrainSetA.conf\"\nIPKNOT_BIN = \"/kaggle/working/ipknot/build/ipknot\"\nVIENNA_LIB_PATH = \"/usr/local/lib/python3.12/dist-packages/viennarna.libs\"\n\ndef run_mxfold2(sequence):\n \"\"\"Runs MxFold2 and returns the dot-bracket string.\"\"\"\n with tempfile.NamedTemporaryFile(mode='w', suffix='.fa', delete=True) as tmp:\n tmp.write(f\">temp\\n{sequence}\\n\")\n tmp.flush()\n \n cmd = [\"mxfold2\", \"predict\", MXFOLD2_CONFIG, tmp.name]\n result = subprocess.run(cmd, capture_output=True, text=True)\n \n if result.returncode == 0:\n lines = result.stdout.strip().split('\\n')\n if len(lines) >= 3:\n # Line 3 format: \"((...)) (energy)\"\n # We split by space and take the first part\n structure_line = lines[2].strip()\n dot_bracket = structure_line.split(' ')[0]\n return dot_bracket\n print(f\"โŒ MxFold2 Error: {result.stderr}\")\n return \".\" * len(sequence)\n\ndef run_ipknot(sequence):\n \"\"\"Runs IPKnot with the required LD_LIBRARY_PATH environment variable.\"\"\"\n # Ensure the environment finds the hidden ViennaRNA libs\n env = os.environ.copy()\n env['LD_LIBRARY_PATH'] = f\"{VIENNA_LIB_PATH}:\" + env.get('LD_LIBRARY_PATH', '')\n \n with tempfile.NamedTemporaryFile(mode='w', suffix='.fa', delete=True) as tmp:\n tmp.write(f\">temp\\n{sequence}\\n\")\n tmp.flush()\n \n cmd = [IPKNOT_BIN, tmp.name]\n result = subprocess.run(cmd, env=env, capture_output=True, text=True)\n \n if result.returncode == 0:\n lines = result.stdout.strip().split('\\n')\n if len(lines) >= 3:\n # Line 3 is the structure\n return lines[2].strip()\n print(f\"โŒ IPKnot Error: {result.stderr}\")\n return \".\" * len(sequence)","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:48.227255Z","iopub.status.busy":"2026-03-27T18:12:48.226473Z","iopub.status.idle":"2026-03-27T18:12:48.23728Z","shell.execute_reply":"2026-03-27T18:12:48.236334Z"},"papermill":{"duration":0.034625,"end_time":"2026-03-27T18:12:48.239093","exception":false,"start_time":"2026-03-27T18:12:48.204468","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# FIX: force offline resolution (was missing --no-index/--no-deps,\n# which makes pip try to reach PyPI for dependency resolution and hang/fail air-gapped)\n!pip install /kaggle/input/datasets/kami1976/biopython-cp312/biopython-1.86-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl --no-index --no-deps\n","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:48.283414Z","iopub.status.busy":"2026-03-27T18:12:48.282264Z","iopub.status.idle":"2026-03-27T18:12:53.583932Z","shell.execute_reply":"2026-03-27T18:12:53.582958Z"},"papermill":{"duration":5.326578,"end_time":"2026-03-27T18:12:53.586696","exception":false,"start_time":"2026-03-27T18:12:48.260118","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport os\nfrom Bio import AlignIO\n\ndef generate_msa_covariance(msa_dir, target_id, original_target_id):\n \"\"\"\n Phase 1: Generates co-evolution signal from a directory of MSAs.\n Constructs path: {msa_dir}/{target_id}.MSA.fasta\n \"\"\"\n # 1. Construct and Verify Path\n file_path = os.path.join(msa_dir, f\"{original_target_id}.MSA.fasta\")\n \n if not os.path.exists(file_path):\n # FIX: report the ACTUAL path/target we look up (original_target_id, e.g. 8ZNQ),\n # not the per-chain label (8ZNQ_c0). The lookup itself already used the bare target.\n print(f\"โš ๏ธ MSA not found for {original_target_id} at {file_path}. Returning empty signal.\")\n # We can't determine the sequence length yet, so we return None \n # to let the consensus function handle it.\n return None\n\n # 2. Load MSA\n try:\n alignment = AlignIO.read(file_path, \"fasta\")\n except Exception as e:\n print(f\"โŒ Error reading MSA for {original_target_id} ({file_path}): {e}\")\n return None\n\n msa_array = np.array([list(rec.seq) for rec in alignment], dtype='U1')\n num_seqs, seq_len = msa_array.shape\n mi_matrix = np.zeros((seq_len, seq_len))\n\n # 3. Calculate Mutual Information (MI)\n for i in range(seq_len):\n for j in range(i + 1, seq_len):\n cols = msa_array[:, [i, j]]\n valid = (cols[:, 0] != '-') & (cols[:, 1] != '-')\n if np.sum(valid) < (num_seqs * 0.1): \n continue\n \n pair_data = cols[valid]\n (unique_pairs, counts) = np.unique(pair_data, axis=0, return_counts=True)\n p_xy = counts / len(pair_data)\n \n # Marginal probabilities\n _, counts_x = np.unique(pair_data[:, 0], return_counts=True)\n p_x = counts_x / len(pair_data)\n _, counts_y = np.unique(pair_data[:, 1], return_counts=True)\n p_y = counts_y / len(pair_data)\n \n mi = 0\n # Pre-calculate unique bases for index mapping\n u_x = np.unique(pair_data[:, 0])\n u_y = np.unique(pair_data[:, 1])\n\n for k, (base_x, base_y) in enumerate(unique_pairs):\n idx_x = np.where(u_x == base_x)[0][0]\n idx_y = np.where(u_y == base_y)[0][0]\n mi += p_xy[k] * np.log(p_xy[k] / (p_x[idx_x] * p_y[idx_y]))\n \n mi_matrix[i, j] = mi_matrix[j, i] = mi\n\n # 4. Apply APC (Average Product Correction)\n apc_matrix = np.zeros_like(mi_matrix)\n mean_row = np.mean(mi_matrix, axis=1)\n mean_all = np.mean(mi_matrix)\n \n # Avoid division by zero if MI matrix is all zeros\n if mean_all == 0: return mi_matrix\n\n for i in range(seq_len):\n for j in range(seq_len):\n if i != j:\n apc_matrix[i, j] = mi_matrix[i, j] - (mean_row[i] * mean_row[j] / mean_all)\n \n return apc_matrix","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:53.632631Z","iopub.status.busy":"2026-03-27T18:12:53.632108Z","iopub.status.idle":"2026-03-27T18:12:53.699111Z","shell.execute_reply":"2026-03-27T18:12:53.698368Z"},"papermill":{"duration":0.092577,"end_time":"2026-03-27T18:12:53.70123","exception":false,"start_time":"2026-03-27T18:12:53.608653","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def get_consensus_ss(seq, target_id, msa_dir, original_target_id):\n \"\"\"Phase 2: Consensus 2D Folding with Directory-based MSA lookup.\"\"\"\n \n # Get 2D predictions\n ss_mxfold = run_mxfold2(seq)\n ss_ipknot = run_ipknot(seq)\n \n # Phase 1: Signal Generation\n cov_signal = generate_msa_covariance(msa_dir, target_id, original_target_id)\n\n # FIX: the covariance matrix is indexed by MSA-alignment columns; those only\n # line up 1:1 with query residues when the matrix width == len(seq). If it does\n # not, cov_signal[i, partner] mis-indexes or raises IndexError (which would kill\n # the whole run). Disable the signal unless the shapes match.\n if cov_signal is not None and getattr(cov_signal, 'shape', (0,))[0] != len(seq):\n cov_signal = None\n \n consensus = list(\".\" * len(seq))\n \n for i in range(len(seq)):\n m_char = ss_mxfold[i]\n ip_char = ss_ipknot[i]\n \n if m_char == ip_char:\n consensus[i] = m_char\n else:\n # If no MSA signal was found, default to '.' for disagreements\n if cov_signal is None:\n consensus[i] = \".\"\n continue\n \n # Tie-breaker logic (find_partner helper is required here)\n m_partner = find_partner(ss_mxfold, i)\n ip_partner = find_partner(ss_ipknot, i)\n \n # Extract scores from the matrix\n score_m = cov_signal[i, m_partner] if m_partner != -1 else 0\n score_ip = cov_signal[i, ip_partner] if ip_partner != -1 else 0\n \n if score_m > score_ip and score_m > 0.05:\n consensus[i] = m_char\n elif score_ip > score_m and score_ip > 0.05:\n consensus[i] = ip_char\n else:\n consensus[i] = \".\"\n \n return \"\".join(consensus)\n\ndef find_partner(dot_bracket, index):\n \"\"\"\n Finds the index of the base-pairing partner in a dot-bracket string.\n Returns -1 if the base is unpaired ('.').\n \"\"\"\n if dot_bracket[index] == '.':\n return -1\n \n char = dot_bracket[index]\n # Define pairs and direction\n pairs = { '(': ')', '[': ']', '{': '}', '<': '>',\n ')': '(', ']': '[', '}': '{', '>': '<' }\n \n target = pairs[char]\n step = 1 if char in '([{<' else -1\n \n counter = 0\n for i in range(index, len(dot_bracket) if step == 1 else -1, step):\n if dot_bracket[i] == char:\n counter += 1\n elif dot_bracket[i] == target:\n counter -= 1\n \n if counter == 0:\n return i\n \n return -1","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:53.746773Z","iopub.status.busy":"2026-03-27T18:12:53.746265Z","iopub.status.idle":"2026-03-27T18:12:53.756676Z","shell.execute_reply":"2026-03-27T18:12:53.755774Z"},"papermill":{"duration":0.035565,"end_time":"2026-03-27T18:12:53.7585","exception":false,"start_time":"2026-03-27T18:12:53.722935","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"get_consensus_ss(\"ACCGUGACGGGCCUUUUGGCUAUACGCGGU\", \"8ZNQ\", \"/kaggle/input/stanford-rna-3d-folding-2/MSA\", \"8ZNQ\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:12:53.803324Z","iopub.status.busy":"2026-03-27T18:12:53.802754Z","iopub.status.idle":"2026-03-27T18:13:05.872343Z","shell.execute_reply":"2026-03-27T18:13:05.87135Z"},"papermill":{"duration":12.09461,"end_time":"2026-03-27T18:13:05.874588","exception":false,"start_time":"2026-03-27T18:12:53.779978","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## DRFold2","metadata":{"papermill":{"duration":0.021595,"end_time":"2026-03-27T18:13:05.919152","exception":false,"start_time":"2026-03-27T18:13:05.897557","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# The path will be /kaggle/input/[dataset-name]/wheels\n# Replace 'biopython-offline-files' with whatever you named your dataset\n!pip install biopython --no-index --find-links=/kaggle/input/datasets/khushvaghela/drfold2","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:05.963502Z","iopub.status.busy":"2026-03-27T18:13:05.96304Z","iopub.status.idle":"2026-03-27T18:13:09.861107Z","shell.execute_reply":"2026-03-27T18:13:09.860132Z"},"papermill":{"duration":3.923793,"end_time":"2026-03-27T18:13:09.863848","exception":false,"start_time":"2026-03-27T18:13:05.940055","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport shutil\nimport time\nimport subprocess\nimport torch\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nfrom Bio.Align import PairwiseAligner \n\n# --- PHASE 1: ENVIRONMENT SETUP ---\nWORKING_DIR = Path(\"/kaggle/working/DRfold2\")\n\ndef setup_environment():\n print(\"โณ Reconstructing DRfold2 environment...\")\n if WORKING_DIR.exists(): shutil.rmtree(WORKING_DIR)\n \n # 1. Copy Repo and Weights\n shutil.copytree(\"/kaggle/input/datasets/jaejohn/drfold2-repo/DRfold2\", WORKING_DIR)\n os.makedirs(WORKING_DIR / \"model_hub\", exist_ok=True)\n !cp -r /kaggle/input/datasets/jaejohn/drfold2/model_hub/* /kaggle/working/DRfold2/model_hub/\n \n # 2. Compile Arena Refiner with g++ (standard on Kaggle)\n %cd {WORKING_DIR}/Arena\n !g++ -O3 Arena.cpp -o Arena \n %cd /kaggle/working/\n print(\"โœ… Environment Ready.\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:09.908539Z","iopub.status.busy":"2026-03-27T18:13:09.908021Z","iopub.status.idle":"2026-03-27T18:13:12.262313Z","shell.execute_reply":"2026-03-27T18:13:12.261317Z"},"papermill":{"duration":2.379346,"end_time":"2026-03-27T18:13:12.264602","exception":false,"start_time":"2026-03-27T18:13:09.885256","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"setup_environment()","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:12.309795Z","iopub.status.busy":"2026-03-27T18:13:12.308729Z","iopub.status.idle":"2026-03-27T18:13:42.305091Z","shell.execute_reply":"2026-03-27T18:13:42.303841Z"},"papermill":{"duration":30.020668,"end_time":"2026-03-27T18:13:42.306932","exception":false,"start_time":"2026-03-27T18:13:12.286264","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile /kaggle/working/DRfold2/DRfold_infer.py\nimport os,sys\nimport torch\nimport numpy as np\nfrom subprocess import Popen, PIPE, STDOUT\n\n# Get the directory where the script is located\nexp_dir = os.path.dirname(os.path.abspath(__file__))\n\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n# dlexps = ['cfg_95','cfg_96','cfg_97','cfg_99']\ndlexps = ['cfg_97']\n\nprint(f\"[DRfold2] Starting prediction pipeline on {device} device\")\n\n# Get input FASTA file and output directory from command line arguments\nfastafile = os.path.realpath(sys.argv[1])\noutdir = os.path.realpath(sys.argv[2])\n\nprint(f\"[DRfold2] Input: {fastafile}\")\nprint(f\"[DRfold2] Output: {outdir}\")\n\n# Initialize clustering flag\npclu = False\n\n# If third argument is '1', enable clustering\nif len(sys.argv) == 4 and sys.argv[3] == '1': \n print('[DRfold2] Clustering enabled - will generate multiple models')\n pclu = True\nelse:\n print('[DRfold2] Clustering disabled - will generate single model')\n\n# Create output directory if it doesn't exist\nif not os.path.isdir(outdir):\n os.makedirs(outdir)\n print(f\"[DRfold2] Created output directory: {outdir}\")\n\n# Create subdirectories for different outputs\nret_dir = os.path.join(outdir,'rets_dir') # For return files\nif not os.path.isdir(ret_dir):\n os.makedirs(ret_dir)\n print(f\"[DRfold2] Created returns directory: {ret_dir}\")\n\nfolddir = os.path.join(outdir,'folds') # For folded structures\nif not os.path.isdir(folddir):\n os.makedirs(folddir)\n print(f\"[DRfold2] Created folds directory: {folddir}\")\n\nrefdir = os.path.join(outdir,'relax') # For relaxed structures\nif not os.path.isdir(refdir):\n os.makedirs(refdir)\n print(f\"[DRfold2] Created relaxation directory: {refdir}\")\n\n# Helper function to run commands and capture output\ndef run_cmd(cmd, description):\n print(f\"[DRfold2] {description}\")\n print(f\"[DRfold2] Command: {cmd}\")\n \n # Execute the command and capture output in real-time\n process = Popen(cmd, shell=True, stdout=PIPE, stderr=STDOUT, universal_newlines=True, bufsize=1)\n \n # Print output line by line as it becomes available\n for line in iter(process.stdout.readline, ''):\n line = line.strip()\n if line:\n print(f\"[DRfold2 subprocess] {line}\")\n \n # Get return code\n return_code = process.wait()\n if return_code == 0:\n print(f\"[DRfold2] {description} completed successfully\")\n else:\n print(f\"[DRfold2] {description} failed with return code {return_code}\")\n return return_code\n\n# Create paths for model directories and test scripts\ndlmains = [os.path.join(exp_dir, one_exp, 'test_modeldir.py') for one_exp in dlexps]\ndirs = [os.path.join(exp_dir, 'model_hub', one_exp) for one_exp in dlexps]\n\n# Check if processing has been done before\nif not os.path.isfile(ret_dir + '/done'): \n print(\"[DRfold2] Step 1/4: GENERATING INITIAL PREDICTIONS\")\n print(f\"[DRfold2] No previous predictions found, will generate e2e and geo files\")\n \n # Run each model configuration\n for idx, (dlmain, one_exp, mdir) in enumerate(zip(dlmains, dlexps, dirs)):\n # Construct command to run the model\n cmd = f'python {dlmain} {device} {fastafile} {ret_dir}/{one_exp}_ {mdir}'\n description = f\"Running model {idx+1}/{len(dlexps)}: {one_exp}\"\n run_cmd(cmd, description)\n\n # Mark processing as complete\n wfile = open(ret_dir+'/done','w')\n wfile.write('1')\n wfile.close()\n print(\"[DRfold2] Initial predictions generation completed\")\nelse:\n print(\"[DRfold2] Step 1/4: USING EXISTING PREDICTIONS\")\n print(f\"[DRfold2] Found previous predictions in {ret_dir}, using existing e2e and geo files\")\n\n# Helper function to get model PDB file\ndef get_model_pdb(tdir,opt):\n files = os.listdir(tdir)\n files = [afile for afile in files if afile.startswith(opt)][0]\n return files\n\n# Set up directory paths and configuration files\ncso_dir = folddir # Directory for coarse-grained structures\nclufile = os.path.join(folddir,'clu.txt') # Clustering results file\nconfig_sel = os.path.join(exp_dir,'cfg_for_selection.json') # Selection configuration\nfoldconfig = os.path.join(exp_dir,'cfg_for_folding.json') # Folding configuration\nselpython = os.path.join(exp_dir,'PotentialFold','Selection.py') # Selection script\noptpython = os.path.join(exp_dir,'PotentialFold','Optimization.py') # Optimization script\nclupy = os.path.join(exp_dir,'PotentialFold','Clust.py') # Clustering script\narena = os.path.join(exp_dir,'Arena','Arena') # Arena executable for structure refinement\n\n# Set up initial save prefixes for optimization and selection\noptsaveprefix = os.path.join(cso_dir, f'opt_0')\nsave_prefix = os.path.join(cso_dir, f'sel_0')\n\n# Get all .ret files from the return directory\nrets = os.listdir(ret_dir)\nrets = [afile for afile in rets if afile.endswith('.ret')]\nrets = [os.path.join(ret_dir,aret) for aret in rets ]\nret_str = ' '.join(rets)\n\nprint(\"[DRfold2] Step 2/4: SELECTION PROCESS\")\nprint(f\"[DRfold2] Found {len(rets)} return files for selection\")\nprint(f\"[DRfold2] Using selection config: {config_sel}\")\nprint(f\"[DRfold2] Output prefix: {save_prefix}\")\n\n# Run selection process\ncmd = f'python {selpython} {fastafile} {config_sel} {save_prefix} {ret_str}'\nrun_cmd(cmd, \"Running selection process\")\n\nprint(\"[DRfold2] Step 3/4: OPTIMIZATION PROCESS\")\nprint(f\"[DRfold2] Using fold config: {foldconfig}\")\nprint(f\"[DRfold2] Optimization output prefix: {optsaveprefix}\")\n\n# Run optimization process\ncmd = f'python {optpython} {fastafile} {optsaveprefix} {ret_dir} {save_prefix} {foldconfig}'\nrun_cmd(cmd, \"Running optimization process\")\n\n# Get the coarse-grained PDB and save refined structure\ncgpdb = os.path.join(folddir,get_model_pdb(folddir,'opt_0'))\nsavepdb = os.path.join(refdir,'model_1.pdb')\n\nprint(\"[DRfold2] Step 4/4: STRUCTURE REFINEMENT\")\nprint(f\"[DRfold2] Found optimized structure: {cgpdb}\")\nprint(f\"[DRfold2] Final output will be saved to: {savepdb}\")\n\ncmd = f'{arena} {cgpdb} {savepdb} 7'\nrun_cmd(cmd, \"Running structure refinement\")\n\n# If clustering is enabled (pclu=True)\nif pclu:\n print(\"[DRfold2] ADDITIONAL STEP: CLUSTERING\")\n print(f\"[DRfold2] Running clustering process, output: {clufile}\")\n \n # Run clustering process\n cmd = f'python {clupy} {ret_dir} {clufile}'\n run_cmd(cmd, \"Running clustering\")\n\n # Read clustering results\n lines = open(clufile).readlines()\n lines = [aline.strip() for aline in lines]\n lines = [aline for aline in lines if aline]\n \n cluster_count = len(lines) - 1\n print(f\"[DRfold2] Found {cluster_count} additional clusters to process\")\n\n # Process each cluster\n for i in range(1,len(lines)):\n print(f\"[DRfold2] PROCESSING CLUSTER {i}/{cluster_count}\")\n \n # Get return files for this cluster\n rets = lines[i].split()\n rets = [os.path.join(ret_dir,aret.replace('.pdb','.ret')) for aret in rets ]\n ret_str = ' '.join(rets)\n\n # Set up save prefixes for this cluster\n optsaveprefix = os.path.join(cso_dir,f'opt_{str(i+1)}')\n save_prefix = os.path.join(cso_dir,f'sel_{str(i+1)}')\n \n print(f\"[DRfold2] Cluster {i} Selection Process\")\n print(f\"[DRfold2] Found {len(rets)} return files for selection\")\n print(f\"[DRfold2] Selection output prefix: {save_prefix}\")\n\n # Run selection process for this cluster\n cmd = f'python {selpython} {fastafile} {config_sel} {save_prefix} {ret_str}'\n run_cmd(cmd, f\"Running selection for cluster {i}\")\n \n print(f\"[DRfold2] Cluster {i} Optimization Process\")\n print(f\"[DRfold2] Optimization output prefix: {optsaveprefix}\")\n\n # Run optimization process for this cluster\n cmd = f'python {optpython} {fastafile} {optsaveprefix} {ret_dir} {save_prefix} {foldconfig}'\n run_cmd(cmd, f\"Running optimization for cluster {i}\")\n\n # Get the coarse-grained PDB and save refined structure for this cluster\n cgpdb = os.path.join(folddir,get_model_pdb(folddir,f'opt_{str(i+1)}'))\n savepdb = os.path.join(refdir,f'model_{str(i+1)}.pdb')\n \n print(f\"[DRfold2] Cluster {i} Refinement Process\")\n print(f\"[DRfold2] Found optimized structure: {cgpdb}\")\n print(f\"[DRfold2] Final output will be saved to: {savepdb}\")\n\n cmd = f'{arena} {cgpdb} {savepdb} 7'\n run_cmd(cmd, f\"Running refinement for cluster {i}\")\n\nprint(\"[DRfold2] PREDICTION PIPELINE COMPLETED SUCCESSFULLY\")\n","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:42.353117Z","iopub.status.busy":"2026-03-27T18:13:42.352674Z","iopub.status.idle":"2026-03-27T18:13:42.363874Z","shell.execute_reply":"2026-03-27T18:13:42.362907Z"},"papermill":{"duration":0.036875,"end_time":"2026-03-27T18:13:42.365619","exception":false,"start_time":"2026-03-27T18:13:42.328744","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile /kaggle/working/DRfold2/PotentialFold/operations.py\n\"\"\"\noperations.py: Core Mathematical Operations for RNA Structure Analysis\n\nThis module provides essential mathematical operations for manipulating and analyzing\nRNA 3D structures, organized into four main categories:\n\n1. Basic Vector Operations:\n Functions for selecting coordinates and calculating distances between points,\n which form the foundation for all structural calculations.\n\n2. Angle Calculations:\n Functions for computing bond angles and dihedral (torsion) angles between atoms,\n with differentiable implementations suitable for gradient-based optimization.\n\n3. Rigid Body Transformations:\n Functions for determining optimal rotations and translations between sets of\n coordinates, enabling structure alignment and manipulation.\n\n4. Sequence Utilities:\n Functions for converting RNA sequence data into standard 3D coordinate templates,\n allowing sequence-structure mapping.\n\nThese operations support the core functionality of RNA structure prediction, analysis,\nand optimization throughout the codebase.\n\"\"\"\n\nimport os\nimport torch\nimport torch.nn as nn\nimport numpy as np \nimport math, sys, math\nfrom io import BytesIO\nimport torch.nn.functional as F\nfrom torch.autograd import Function\nfrom torch.nn.parameter import Parameter\nfrom subprocess import Popen, PIPE, STDOUT\n\n# Use consistent epsilon value across all functions\nEPS = 1e-8\n\n\n# === Basic Vector Operations ===\ndef coor_selection(coor,mask):\n #[L,n,3],[L,n],byte\n return torch.masked_select(coor,mask.bool()).view(-1,3)\n\ndef pair_distance(x1, x2, eps=1e-6, p=2):\n # Use torch.cdist for p=2 (Euclidean) which is highly optimized\n if p == 2:\n return torch.cdist(x1, x2, p=2)\n \n # For other p-norms, avoid memory expansion with broadcasting\n x1_ = x1.unsqueeze(1) # [n1, 1, dim]\n x2_ = x2.unsqueeze(0) # [1, n2, dim]\n diff = torch.abs(x1_ - x2_)\n out = torch.pow(diff + eps, p).sum(dim=2)\n return torch.pow(out, 1. / p)\n\n\n# === Angle Calculations ===\ndef angle(p0, p1, p2):\n # [b 3] \n b0 = p0-p1\n b1 = p2-p1\n\n b0 = b0 / (torch.norm(b0, dim =-1, keepdim=True) + EPS)\n b1 = b1 / (torch.norm(b1, dim =-1, keepdim=True) + EPS)\n \n recos = torch.sum(b0*b1, -1)\n recos = torch.clamp(recos, -0.9999, 0.9999)\n return torch.acos(recos)\n\nclass torsion(Function):\n #PyTorch class to calculate differentiable torsion angle\n #https://stackoverflow.com/questions/20305272/dihedral-torsion-angle-from-four-points-in-cartesian-coordinates-in-python\n #https://salilab.org/modeller/manual/node492.html\n @staticmethod\n def forward(ctx, p0, p1, p2, p3):\n # Save input points for backward pass\n ctx.save_for_backward(p0, p1, p2, p3)\n\n # Calculate bond vectors\n b0 = p0 - p1\n b1 = p2 - p1\n b2 = p3 - p2\n\n # Normalize the middle bond vector\n b1_norm = torch.norm(b1, dim=-1, keepdim=True) + 1e-8\n b1_unit = b1 / b1_norm\n\n # Project the other bonds onto the plane perpendicular to middle bond\n v = b0 - torch.sum(b0 * b1_unit, dim=-1, keepdim=True) * b1_unit\n w = b2 - torch.sum(b2 * b1_unit, dim=-1, keepdim=True) * b1_unit\n\n # Calculate torsion using the arctan2 formula (more stable than arccos)\n x = torch.sum(v * w, dim=-1) # cosine component\n y = torch.sum(torch.cross(b1_unit, v, dim=-1) * w, dim=-1) # sine component\n\n return torch.atan2(y, x)\n\n \n @staticmethod\n def backward(ctx, grad_output):\n # Retrieve saved tensors from forward pass\n p0, p1, p2, p3 = ctx.saved_tensors\n\n # Calculate bond vectors\n r01 = p0 - p1\n r12 = p2 - p1\n r23 = p3 - p2\n\n # Calculate bond lengths with numerical stability\n d01 = torch.norm(r01, dim=-1, keepdim=True) + 1e-8\n d12 = torch.norm(r12, dim=-1, keepdim=True) + 1e-8\n d23 = torch.norm(r23, dim=-1, keepdim=True) + 1e-8\n\n # Normalize bond vectors\n e01 = r01 / d01\n e12 = r12 / d12\n e23 = r23 / d23\n\n # Calculate normal vectors to the two planes\n n1 = torch.cross(e01, e12, dim=-1)\n n2 = torch.cross(e12, e23, dim=-1)\n\n # Normalize normal vectors\n n1_norm = torch.norm(n1, dim=-1, keepdim=True) + 1e-8\n n2_norm = torch.norm(n2, dim=-1, keepdim=True) + 1e-8\n n1 = n1 / n1_norm\n n2 = n2 / n2_norm\n\n # Calculate gradients for each atom\n # These are based on the analytical derivatives of dihedral angles\n g0 = torch.cross(e01, n1, dim=-1) / d01\n g1 = -g0 - torch.cross(e12, n1, dim=-1) / d12\n g2 = torch.cross(e12, n2, dim=-1) / d12 - torch.cross(e23, n2, dim=-1) / d23\n g3 = torch.cross(e23, n2, dim=-1) / d23\n\n # Apply chain rule with incoming gradient\n g0 = g0 * grad_output.unsqueeze(-1)\n g1 = g1 * grad_output.unsqueeze(-1)\n g2 = g2 * grad_output.unsqueeze(-1)\n g3 = g3 * grad_output.unsqueeze(-1)\n\n return g0, g1, g2, g3\n\n\ndef dihedral(input1, input2, input3, input4):\n return torsion.apply(input1, input2, input3, input4)\n\n\n\n# === Rigid Body Transformations ===\ndef rigidFrom3Points(x): \n x1, x2, x3 = x[:, 0], x[:, 1], x[:, 2]\n v1 = x3 - x2\n v2 = x1 - x2\n \n # Normalize v1 to get e1\n e1 = F.normalize(v1, p=2, dim=-1)\n \n # Project v2 onto e1 and subtract to get the component orthogonal to e1\n u2 = v2 - e1 * (torch.einsum('bn,bn->b', e1, v2)[:, None])\n \n # Normalize u2 to get e2\n e2 = F.normalize(u2, p=2, dim=-1)\n \n # Cross product to get e3\n e3 = torch.cross(e1, e2, dim=-1)\n \n return torch.stack([e1, e2, e3], dim=1)\n\n\n# return the direction from to_q to from_p\ndef Kabsch_rigid(bases,x1,x2,x3):\n # Early return for empty input\n if x1.shape[0] == 0:\n return torch.empty(0, 3, 3), torch.empty(0, 3)\n \n the_dim=1\n to_q = torch.stack([x1,x2,x3],dim=the_dim)\n biasq=torch.mean(to_q,dim=the_dim,keepdim=True)\n q=to_q-biasq\n m = torch.einsum('bnz,bny->bzy',bases,q)\n u, s, v = torch.svd(m)\n vt = torch.transpose(v, 1, 2)\n det = torch.det(torch.matmul(u, vt))\n det = det.view(-1, 1, 1)\n vt = torch.cat((vt[:, :2, :], vt[:, -1:, :] * det), 1)\n r = torch.matmul(u, vt)\n return r,biasq.squeeze()\n\n\n\n# === Sequence Utilities ===\ndef Get_base(seq,basenpy_standard):\n base_num = basenpy_standard.shape[1]\n basenpy = np.zeros([len(seq),base_num,3])\n seqnpy = np.array(list(seq))\n basenpy[seqnpy=='A']=basenpy_standard[0]\n basenpy[seqnpy=='a']=basenpy_standard[0]\n\n basenpy[seqnpy=='G']=basenpy_standard[1]\n basenpy[seqnpy=='g']=basenpy_standard[1]\n\n basenpy[seqnpy=='C']=basenpy_standard[2]\n basenpy[seqnpy=='c']=basenpy_standard[2]\n\n basenpy[seqnpy=='U']=basenpy_standard[3]\n basenpy[seqnpy=='u']=basenpy_standard[3]\n\n basenpy[seqnpy=='T']=basenpy_standard[3]\n basenpy[seqnpy=='t']=basenpy_standard[3]\n \n return torch.from_numpy(basenpy).double()","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:42.411175Z","iopub.status.busy":"2026-03-27T18:13:42.410748Z","iopub.status.idle":"2026-03-27T18:13:42.41883Z","shell.execute_reply":"2026-03-27T18:13:42.417904Z"},"papermill":{"duration":0.033123,"end_time":"2026-03-27T18:13:42.420731","exception":false,"start_time":"2026-03-27T18:13:42.387608","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile /kaggle/working/DRfold2/PotentialFold/Optimization.py\n#! /nfs/amino-home/liyangum/miniconda3/bin/python\nimport torch\nimport random\nimport numpy as np \nimport os, json, sys\n\nimport Cubic, Potential\nimport operations\nimport a2b, rigid\nimport torch.optim as opt\nfrom scipy.optimize import minimize\nimport pickle\n\ntorch.manual_seed(6)\nnp.random.seed(9)\nrandom.seed(9)\n\n\nScale_factor = 1.0\nUSEGEO = False\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\n\ndef readconfig(configfile=''):\n config=[]\n expdir=os.path.dirname(os.path.abspath(__file__))\n if configfile=='':\n configfile=os.path.join(expdir,'lib','ddf.json')\n config=json.load(open(configfile,'r'))\n return config \n\n \nclass Structure:\n def __init__(self,fastafile,geofiles,saveprefix,initial_ret,foldconfig):\n self.config=readconfig(foldconfig)\n self.seqfile=fastafile\n self.init_ret = initial_ret\n self.foldconfig = foldconfig\n self.geofiles = geofiles\n self.rets = [pickle.load(open(refile,'rb')) for refile in geofiles]\n self.txs=[]\n for ret in self.rets:\n self.txs.append( torch.from_numpy(ret['coor']).double().to(device))\n self.handle_geo()\n self.pair = []\n for ret in self.rets:\n self.pair.append(torch.from_numpy(ret['plddt']).double().to(device))\n self.saveprefix=saveprefix\n self.seq=open(fastafile).readlines()[1].strip()\n self.L=len(self.seq)\n basenpy = np.load(os.path.join(os.path.dirname(os.path.abspath(__file__)),'lib','base.npy'))\n self.basex = operations.Get_base(self.seq,basenpy).to(device)\n othernpy = np.load(os.path.join(os.path.dirname(os.path.abspath(__file__)),'lib','other2.npy'))\n self.otherx = operations.Get_base(self.seq,othernpy).to(device)\n sidenpy = np.load(os.path.join(os.path.dirname(os.path.abspath(__file__)),'lib','side.npy'))\n self.sidex = operations.Get_base(self.seq,sidenpy).to(device)\n \n self.init_mask()\n self.init_paras()\n self._init_fape()\n self.tx2ds = [td.to(device) for td in self.tx2ds]\n self.local_weight = torch.ones(self.L,self.L).to(device)\n \n for i in range(self.L):\n for j in range(i+1,min(self.L,i+2)):\n self.local_weight[i,j] = self.local_weight[j,i] = 4\n for j in range(i+2,min(self.L,i+3)):\n self.local_weight[i,j] = self.local_weight[j,i] = 3\n for j in range(i+3,min(self.L,i+4)):\n self.local_weight[i,j] = self.local_weight[j,i] = 2\n\n def _init_fape(self):\n self.tx2ds = []\n for tx in self.txs:\n true_rot,true_trans = operations.Kabsch_rigid(self.basex,tx[:,0],tx[:,1],tx[:,2])\n true_x2 = tx[:,None,:,:] - true_trans[None,:,None,:]\n true_x2 = torch.einsum('ijnd,jde->ijne',true_x2,true_rot.transpose(-1,-2))\n self.tx2ds.append(true_x2)\n \n def handle_geo(self):\n oldkeys=['dist_p','dist_c','dist_n']\n newkeys=['pp','cc','nn']\n self.geos=[]\n for ret in self.rets:\n geo = {}\n for nk,ok in zip(newkeys,oldkeys):\n geo[nk] = torch.from_numpy(ret[ok].astype(np.float64)).to(device) + 0\n self.geos.append(geo)\n\n\n def init_mask(self):\n halfmask=np.zeros([self.L,self.L])\n fullmask=np.zeros([self.L,self.L])\n for i in range(self.L):\n for j in range(i+1,self.L):\n halfmask[i,j]=1\n fullmask[i,j]=1\n fullmask[j,i]=1\n self.halfmask=(torch.DoubleTensor(halfmask) > 0.5).to(device)\n self.fullmask=(torch.DoubleTensor(fullmask) > 0.5).to(device)\n self.clash_mask = torch.zeros([self.L,self.L,22,22], device=device)\n for i in range(self.L):\n for j in range(i+1,self.L):\n self.clash_mask[i,j]=1\n\n for i in range(self.L):\n self.clash_mask[i,i,:6,7:]=1\n\n for i in range(self.L-1):\n self.clash_mask[i,i+1,:,0]=0\n self.clash_mask[i,i+1,0,:]=0\n self.clash_mask[i,i+1,:,5]=0\n self.clash_mask[i,i+1,5,:]=0\n\n self.side_mask = rigid.side_mask(self.seq).to(device)\n self.side_mask = (self.side_mask[:,None,:,None] * self.side_mask[None,:,None,:]).to(device)\n self.clash_mask = ((self.clash_mask > 0.5) * (self.side_mask > 0.5)).to(device)\n\n self.geo_confimask_cc = []\n self.geo_confimask_pp = []\n self.geo_confimask_nn = []\n for geo in self.geos:\n confimask_cc = geo['cc'][:,:,-1] < 0.5\n confimask_pp = geo['pp'][:,:,-1] < 0.5\n confimask_nn = geo['nn'][:,:,-1] < 0.5\n self.geo_confimask_cc.append(confimask_cc)\n self.geo_confimask_pp.append(confimask_pp)\n self.geo_confimask_nn.append(confimask_nn)\n\n\n def init_paras(self):\n self.geo_cc = []\n self.geo_pp = []\n self.geo_nn = []\n self.cs_coefs = {'cc': [], 'pp': [], 'nn': []}\n self.cs_knots = {'cc': [], 'pp': [], 'nn': []}\n for geo in self.geos:\n cc_cs,cc_decs=Cubic.dis_cubic(geo['cc'],2,40,36)\n pp_cs,pp_decs=Cubic.dis_cubic(geo['pp'],2,40,36)\n nn_cs,nn_decs=Cubic.dis_cubic(geo['nn'],2,40,36)\n self.geo_cc.append([cc_cs,cc_decs])\n self.geo_pp.append([pp_cs,pp_decs])\n self.geo_nn.append([nn_cs,nn_decs])\n \n L = self.L\n cc_coefs_np = np.stack([[cc_cs[i,j].c for j in range(L)] for i in range(L)], axis=0)\n cc_knots_np = np.stack([[cc_cs[i,j].x for j in range(L)] for i in range(L)], axis=0)\n self.cs_coefs['cc'].append(torch.from_numpy(cc_coefs_np).to(device))\n self.cs_knots['cc'].append(torch.from_numpy(cc_knots_np).to(device))\n \n pp_coefs_np = np.stack([[pp_cs[i,j].c for j in range(L)] for i in range(L)], axis=0)\n pp_knots_np = np.stack([[pp_cs[i,j].x for j in range(L)] for i in range(L)], axis=0)\n self.cs_coefs['pp'].append(torch.from_numpy(pp_coefs_np).to(device))\n self.cs_knots['pp'].append(torch.from_numpy(pp_knots_np).to(device))\n \n nn_coefs_np = np.stack([[nn_cs[i,j].c for j in range(L)] for i in range(L)], axis=0)\n nn_knots_np = np.stack([[nn_cs[i,j].x for j in range(L)] for i in range(L)], axis=0)\n self.cs_coefs['nn'].append(torch.from_numpy(nn_coefs_np).to(device))\n self.cs_knots['nn'].append(torch.from_numpy(nn_knots_np).to(device))\n\n\n def compute_bb_clash(self,coor,other_coor):\n com_coor = torch.cat([coor,other_coor],dim=1)\n com_dis = (com_coor[:,None,:,None,:] - com_coor[None,:,None,:,:]).norm(dim=-1)\n dynamicmask2_vdw= (com_dis <= 3.15) * (self.clash_mask)\n vdw_dynamic = Potential.LJpotential(com_dis[dynamicmask2_vdw],3.15)\n return vdw_dynamic.sum()*self.config['weight_vdw']\n\n def compute_full_clash(self,coor,other_coor,side_coor):\n com_coor = torch.cat([coor[:,:2],other_coor,side_coor],dim=1)\n com_dis = (com_coor[:,None,:,None,:] - com_coor[None,:,None,:,:]).norm(dim=-1)\n dynamicmask2_vdw= (com_dis <= 2.5) * (self.clash_mask)\n vdw_dynamic = Potential.LJpotential(com_dis[dynamicmask2_vdw],2.5)\n return vdw_dynamic.sum()*self.config['weight_vdw']\n\n\n def _cubic_pair_energy(self, atom_map, geo_cs, geo_confimask, weight_key):\n \"\"\"General cubic-spline energy for CC/PP/NN pairs.\"\"\"\n min_dis, max_dis, bin_num = 2, 40, 36\n dev = atom_map.device\n upper_th = max_dis - ((max_dis - min_dis) / bin_num) * 0.5\n lower_th = 2.5\n total = torch.zeros((), device=dev, dtype=torch.double)\n spline_key = weight_key.split('_')[1] # 'cc', 'pp', or 'nn'\n coeffs_list = self.cs_coefs[spline_key]\n knots_list = self.cs_knots[spline_key]\n for block_idx, mask_block in enumerate(geo_confimask):\n mask = (atom_map <= upper_th) & mask_block & self.fullmask & (atom_map >= lower_th)\n idx = mask.nonzero(as_tuple=True)\n if idx[0].numel() > 1:\n coef = coeffs_list[block_idx][idx]\n knots = knots_list[block_idx][idx]\n part1 = Potential.cubic_distance(atom_map[mask], coef, knots, min_dis, max_dis, bin_num).sum() * self.config[weight_key] * 0.5\n else:\n part1 = torch.zeros((), device=dev)\n part2 = ((atom_map <= lower_th) & mask_block & self.fullmask).sum() * self.config[weight_key]\n total = total + part1 + part2\n return total\n\n def compute_cc_energy(self, coor):\n atom_map = operations.pair_distance(coor[:,1], coor[:,1])\n return self._cubic_pair_energy(atom_map, self.geo_cc, self.geo_confimask_cc, 'weight_cc')\n \n def compute_pp_energy(self, coor):\n atom_map = operations.pair_distance(coor[:,0], coor[:,0])\n return self._cubic_pair_energy(atom_map, self.geo_pp, self.geo_confimask_pp, 'weight_pp')\n \n def compute_nn_energy(self, coor):\n atom_map = operations.pair_distance(coor[:,-1], coor[:,-1])\n return self._cubic_pair_energy(atom_map, self.geo_nn, self.geo_confimask_nn, 'weight_nn')\n\n def compute_pccp_energy(self,coor):\n p_atoms=coor[:,0]\n c_atoms=coor[:,1]\n pccpmap=operations.dihedral( p_atoms[self.pccpi], c_atoms[self.pccpi], c_atoms[self.pccpj] ,p_atoms[self.pccpj] )\n neg_log = Potential.cubic_torsion(pccpmap,self.pccp_coe,self.pccp_x,36)\n return neg_log.sum()*self.config['weight_pccp']\n\n def compute_cnnc_energy(self,coor):\n n_atoms=coor[:,-1]\n c_atoms=coor[:,1]\n pccpmap=operations.dihedral( c_atoms[self.cnnci], n_atoms[self.cnnci], n_atoms[self.cnncj] ,c_atoms[self.cnncj] )\n neg_log = Potential.cubic_torsion(pccpmap,self.cnnc_coe,self.cnnc_x,36)\n return neg_log.sum()*self.config['weight_cnnc']\n\n def compute_pnnp_energy(self,coor):\n n_atoms=coor[:,-1]\n p_atoms=coor[:,0]\n pccpmap=operations.dihedral( p_atoms[self.pnnpi], n_atoms[self.pnnpi], n_atoms[self.pnnpj] ,p_atoms[self.pnnpj] )\n neg_log = Potential.cubic_torsion(pccpmap,self.pnnp_coe,self.pnnp_x,36)\n return neg_log.sum()*self.config['weight_pnnp']\n\n def compute_pcc_energy(self,coor):\n p_atoms=coor[:,1]\n c_atoms=coor[:,2]\n pccmap=operations.angle( p_atoms[self.pcci], c_atoms[self.pcci], c_atoms[self.pccj] )\n neg_log = Potential.cubic_angle(pccmap,self.pcc_coe,self.pcc_x,12)\n return neg_log.sum()*self.config['weight_pcc']\n\n def compute_fape_energy(self,coor,ep=1e-3,epmax=20):\n energy= 0\n for tx in self.tx2ds:\n px_mean = coor[:,[1]]\n p_rot = operations.rigidFrom3Points(coor)\n p_tran = px_mean[:,0]\n pred_x2 = coor[:,None,:,:] - p_tran[None,:,None,:] # Lx Lrot N , 3\n pred_x2 = torch.einsum('ijnd,jde->ijne',pred_x2,p_rot.transpose(-1,-2)) # transpose should be equal to inverse\n errmap=torch.sqrt( ((pred_x2 - tx)**2).sum(dim=-1) + ep )\n energy = energy + torch.sum( torch.clamp(errmap,max=epmax) )\n return energy * self.config['weight_fape']\n\n def compute_bond_energy(self,coor,other_coor):\n # 3.87\n o3 = other_coor[:-1,-2]\n p = coor[1:,0]\n dis = (o3-p).norm(dim=-1)\n energy = ((dis-1.607)**2).sum()\n return energy * self.config['weight_bond']\n\n def tooth_func(self,errmap, ep = 0.05):\n return -1/(errmap/10+ep) + (1/ep)\n\n def reweight_func(self,ww):\n reweighting = torch.pow(ww,self.config['pair_weight_power'])\n reweighting[ww < self.config['pair_weight_min']] = 0\n return reweighting\n\n def compute_fape_energy_fromquat(self,x,coor,ep=1e-6,epmax=100):\n energy= 0\n p_rot,px_mean = a2b.Non2rot(x[:,:9],x.shape[0]),x[:,9:]\n pred_x2 = coor[:,None,:,:] - px_mean[None,:,None,:] # Lx Lrot N , 3\n pred_x2 = torch.einsum('ijnd,jde->ijne',pred_x2,p_rot.transpose(-1,-2)) # transpose should be equal to inverse\n for tx,weightplddt in zip(self.tx2ds,self.pair):\n\n tamplate_dist_map = torch.min( tx.norm(dim=-1), dim=2 )[0]\n errmap=torch.sqrt( ((pred_x2 - tx)**2).sum(dim=-1) + ep ) \n energy = energy + torch.sum( ( (torch.clamp(errmap,max=self.config['FAPE_max'])**self.config['pair_error_power']) * self.reweight_func(weightplddt[...,None]) * self.local_weight[...,None] )[tamplate_dist_map>self.config['pair_rest_min_dist']] )\n\n return energy * self.config['weight_fape']\n\n\n def energy(self,rama):\n coor=a2b.quat2b(self.basex,rama[:,9:])\n other_coor = a2b.quat2b(self.otherx,rama[:,9:])\n side_coor = a2b.quat2b(self.sidex,torch.cat([rama[:,:9],coor[:,-1]],dim=-1))\n \n if self.config['weight_cc']>0:\n E_cc= self.compute_cc_energy(coor) / len(self.rets)\n else:\n E_cc=0\n if self.config['weight_pp']>0:\n E_pp= self.compute_pp_energy(coor) / len(self.rets)\n else:\n E_pp=0\n if self.config['weight_nn']>0:\n E_nn= self.compute_nn_energy(coor) / len(self.rets)\n else:\n E_nn=0\n\n if self.config['weight_pccp']>0:\n E_pccp= self.compute_pccp_energy(coor) / len(self.rets)\n else:\n E_pccp=0\n\n if self.config['weight_cnnc']>0:\n E_cnnc= self.compute_cnnc_energy(coor) / len(self.rets)\n else:\n E_cnnc=0\n\n if self.config['weight_pnnp']>0:\n E_pnnp= self.compute_pnnp_energy(coor) / len(self.rets)\n else:\n E_pnnp=0\n\n if self.config['weight_vdw']>0:\n E_vdw= self.compute_full_clash(coor,other_coor,side_coor)\n else:\n E_vdw=0\n\n if self.config['weight_fape']>0:\n E_fape= self.compute_fape_energy_fromquat(rama[:,9:],coor) / len(self.rets)\n else:\n E_fape=0\n if self.config['weight_bond']>0:\n E_bond= self.compute_bond_energy(coor,other_coor)\n else:\n E_bond=0\n return E_vdw + E_fape + E_bond + E_pp + E_cc + E_nn + E_pccp + E_cnnc + E_pnnp\n\n\n def obj_func_grad_np(self,rama_):\n rama=torch.DoubleTensor(rama_)\n rama.requires_grad=True\n if rama.grad:\n rama.grad.zero_()\n f=self.energy(rama.view(self.L,21))*Scale_factor\n grad_value=autograd.grad(f,rama)[0]\n return grad_value.data.numpy().astype(np.float64)\n \n def obj_func_np(self,rama_):\n rama=torch.DoubleTensor(rama_)\n rama=rama.view(self.L,21)\n with torch.no_grad():\n f=self.energy(rama)*Scale_factor\n return f.item()\n\n\n def foldning(self):\n ilter = self.init_ret\n # 1) get initial quaternions (double precision)\n try:\n init_q = self.init_quat(ilter).double()\n except:\n init_q = self.init_quat_safe(ilter).double()\n\n # 2) move to target device (GPU if available), enable grad\n param = init_q.to(device).clone().detach().requires_grad_(True)\n\n # 3) set up PyTorch LBFGS optimizer over `param`\n optimizer = opt.LBFGS(\n [param],\n max_iter=self.config.get('max_iter', 300),\n tolerance_grad=1e-6,\n tolerance_change=1e-9,\n history_size=10,\n line_search_fn='strong_wolfe'\n )\n\n # 4) define the โ€œclosureโ€ that LBFGS will call to reevaluate loss + gradients\n def closure():\n optimizer.zero_grad() # clear old grads\n E = self.energy(param.view(self.L,21)) * Scale_factor # compute โˆ‚E/โˆ‚param\n E.backward()\n return E\n\n # 5) run LBFGS until convergence (it calls closure repeatedly)\n optimizer.step(closure)\n\n # 6) write out final PDB\n final_energy = self.energy(param.view(self.L,21)).item()\n self.outpdb(param, self.saveprefix + '.pdb', energystr=str(final_energy))\n\n\n def outpdb(self,rama,savefile,start=0,end=10000,energystr=''):\n # bring baseframes and quaternion data onto CPU to prevent device mismatch\n basex_cpu = self.basex.detach().cpu()\n otherx_cpu = self.otherx.detach().cpu()\n sidex_cpu = self.sidex.detach().cpu()\n shaped_rama = rama.view(self.L,21).detach().cpu()\n # compute backbone and other coords\n coor_np = a2b.quat2b(basex_cpu, shaped_rama[:,9:]).detach().cpu().numpy()\n other_np = a2b.quat2b(otherx_cpu, shaped_rama[:,9:]).detach().cpu().numpy()\n coor = torch.FloatTensor(coor_np)\n # compute side atom coords\n side_coor_NP = a2b.quat2b(sidex_cpu, torch.cat([shaped_rama[:,:9], coor[:,-1]], dim=-1)).detach().cpu().numpy()\n \n Atom_name=[' P ',\" C4'\",' N1 ']\n Other_Atom_name = [\" O5'\",\" C5'\",\" C3'\",\" O3'\",\" C1'\"]\n other_last_name = ['O',\"C\",\"C\",\"O\",\"C\"]\n\n side_atoms= [' N1 ',' C2 ',' O2 ',' N2 ',' N3 ',' N4 ',' C4 ',' O4 ',' C5 ',' C6 ',' O6 ',' N6 ',' N7 ',' N8 ',' N9 ']\n side_last_name = ['N', \"C\", \"O\", \"N\", \"N\", 'N', 'C', 'O', 'C', 'C', 'O', 'N', 'N', 'N','N']\n\n base_dict = rigid.base_table()\n last_name=['P','C','N']\n wstr=[f'REMARK {str(energystr)}']\n templet='%6s%5d %4s %3s %1s%4d %8.3f%8.3f%8.3f%6.2f%6.2f %2s%2s'\n count=1\n for i in range(self.L):\n if self.seq[i] in ['a','g','A','G']:\n Atom_name = [' P ',\" C4'\",' N9 ']\n #atoms = ['P','C4']\n\n elif self.seq[i] in ['c','u','C','U']:\n Atom_name = [' P ',\" C4'\",' N1 ']\n for j in range(coor_np.shape[1]):\n outs=('ATOM ',count,Atom_name[j],self.seq[i],'A',i+1,coor_np[i][j][0],coor_np[i][j][1],coor_np[i][j][2],0,0,last_name[j],'')\n if i>=start-1 and i < end:\n wstr.append(templet % outs)\n count+=1\n\n for j in range(other_np.shape[1]):\n outs=('ATOM ',count,Other_Atom_name[j],self.seq[i],'A',i+1,other_np[i][j][0],other_np[i][j][1],other_np[i][j][2],0,0,other_last_name[j],'')\n if i>=start-1 and i < end:\n wstr.append(templet % outs)\n count+=1\n \n wstr='\\n'.join(wstr)\n wfile=open(savefile,'w')\n wfile.write(wstr)\n wfile.close()\n \n def outpdb_coor(self,coor_np,savefile,start=0,end=1000,energystr=''):\n Atom_name=[' P ',\" C4'\",' N1 ']\n last_name=['P','C','N']\n wstr=[f'REMARK {str(energystr)}']\n templet='%6s%5d %4s %3s %1s%4d %8.3f%8.3f%8.3f%6.2f%6.2f %2s%2s'\n count=1\n for i in range(self.L):\n if self.seq[i] in ['a','g','A','G']:\n Atom_name = [' P ',\" C4'\",' N9 ']\n\n elif self.seq[i] in ['c','u','C','U']:\n Atom_name = [' P ',\" C4'\",' N1 ']\n for j in range(coor_np.shape[1]):\n outs=('ATOM ',count,Atom_name[j],self.seq[i],'A',i+1,coor_np[i][j][0],coor_np[i][j][1],coor_np[i][j][2],0,0,last_name[j],'')\n if i>=start-1 and i < end:\n wstr.append(templet % outs)\n count+=1\n \n wstr='\\n'.join(wstr)\n wfile=open(savefile,'w')\n wfile.write(wstr)\n wfile.close()\n\n\n def init_quat(self,ii):\n x = torch.rand([self.L,21])\n x[:,18:] = self.txs[ii].mean(dim=1)\n init_coor = self.txs[ii]\n biasq = torch.mean(init_coor,dim=1,keepdim=True)\n q = init_coor - biasq\n m = torch.einsum('bnz,bny->bzy',self.basex,q).reshape([self.L,-1])\n x[:,:9] = x[:,9:18] = m\n x.requires_grad_()\n return x\n\n def init_quat_safe(self,ii):\n x = torch.rand([self.L,21])\n x[:,18:] = self.txs[ii].mean(dim=1)\n init_coor = self.txs[ii]\n biasq = torch.mean(init_coor,dim=1,keepdim=True)\n q = init_coor - biasq + torch.rand([self.L,3,3])\n m = (torch.einsum('bnz,bny->bzy',self.basex,q) + torch.eye(3)[None,:,:]).reshape([self.L,-1])\n x[:,:9] = x[:,9:18] = m\n x.requires_grad_()\n return x\n\n\nif __name__ == '__main__': \n\n fastafile=sys.argv[1]\n saveprefix=sys.argv[2]\n retdirs =sys.argv[3]\n ret_score = sys.argv[4]\n foldconfig = sys.argv[5]\n\n savepare = os.path.dirname(saveprefix)\n if not os.path.isdir(savepare):\n os.makedirs(savepare)\n\n num_of_models = readconfig(foldconfig)['num_of_models']\n\n score_dict = readconfig(ret_score)\n sorted_items = sorted(score_dict.items(), key=lambda x: x[1])\n lowest_n_keys = [item[0] for item in sorted_items][:num_of_models]\n bestkey = lowest_n_keys[0] + ''\n print(\"Before sort:\", lowest_n_keys)\n lowest_n_keys.sort()\n print(\"After sort:\", lowest_n_keys)\n bestindex = lowest_n_keys.index(bestkey)\n\n current_ret = bestkey\n retfiles = [os.path.join(retdirs, afile) for afile in lowest_n_keys]\n stru = Structure(fastafile, retfiles, saveprefix + '_from_' + current_ret, bestindex, foldconfig)\n stru.foldning()","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:42.467781Z","iopub.status.busy":"2026-03-27T18:13:42.467474Z","iopub.status.idle":"2026-03-27T18:13:42.484166Z","shell.execute_reply":"2026-03-27T18:13:42.483298Z"},"papermill":{"duration":0.042904,"end_time":"2026-03-27T18:13:42.485894","exception":false,"start_time":"2026-03-27T18:13:42.44299","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile /kaggle/working/DRfold2/PotentialFold/Selection.py\n#! /nfs/amino-home/liyangum/miniconda3/bin/python\nimport numpy\nimport torch\nimport torch.autograd as autograd\nimport numpy as np \n\nimport random\nimport Cubic, Potential\nimport operations\nimport os, json, sys\n\nimport a2b, rigid\nimport torch.optim as opt\nfrom torch.nn.parameter import Parameter\nimport torch.nn as nn\nimport math\nfrom scipy.optimize import fmin_l_bfgs_b,fmin_cg,fmin_bfgs\nfrom scipy.optimize import minimize\nimport lbfgs_rosetta\nimport pickle\nimport shutil\n\ntorch.manual_seed(6)\ntorch.set_num_threads(4)\nnp.random.seed(9)\nrandom.seed(9)\n\nScale_factor = 1.0\nUSEGEO = False\n\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n\ndef readconfig(configfile=''):\n config=[]\n expdir=os.path.dirname(os.path.abspath(__file__))\n if configfile=='':\n configfile=os.path.join(expdir,'lib','ddf.json')\n config=json.load(open(configfile,'r'))\n return config \n\n \nclass Structure:\n def __init__(self, fastafile, geofiles, foldconfig, saveprefix):\n # Load Configuration and Inputs\n self.config = readconfig(foldconfig)\n self.seqfile = fastafile\n self.foldconfig = foldconfig\n self.geofiles = geofiles\n\n # Load Model Results\n self.rets = [pickle.load(open(refile, 'rb')) for refile in geofiles]\n \n # Extract Coordinates\n self.txs = []\n for ret in self.rets:\n self.txs.append(torch.from_numpy(ret['coor']).double().to(device))\n \n # Handle Geometrical Data\n self.handle_geo()\n\n # Extract pLDDT Scores\n self.pair = []\n for ret in self.rets:\n self.pair.append( torch.from_numpy(ret['plddt']).double().to(device))\n \n # Store Output and Sequence Info\n self.saveprefix = saveprefix\n self.seq = open(fastafile).readlines()[1].strip()\n self.L = len(self.seq)\n \n # Load Reference Arrays for Structure Construction\n basenpy = np.load(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'lib', 'base.npy'))\n self.basex = operations.Get_base(self.seq, basenpy).double().to(device)\n \n othernpy = np.load(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'lib', 'other2.npy'))\n self.otherx = operations.Get_base(self.seq, othernpy).double().to(device)\n \n sidenpy = np.load(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'lib', 'side.npy'))\n self.sidex = operations.Get_base(self.seq, sidenpy).double().to(device) \n \n # Initialize Masks, Parameters, and FAPE\n self.init_mask()\n self.init_paras()\n self._init_fape()\n \n\n def _init_fape(self):\n self.tx2ds = []\n for tx in self.txs:\n true_rot, true_trans = operations.Kabsch_rigid(self.basex, tx[:, 0], tx[:, 1], tx[:, 2])\n true_x2 = tx[:, None, :, :] - true_trans[None, :, None, :]\n true_x2 = torch.einsum('ijnd,jde->ijne', true_x2, true_rot.transpose(-1,-2))\n self.tx2ds.append(true_x2)\n \n\n def handle_geo(self):\n oldkeys = ['dist_p', 'dist_c', 'dist_n']\n newkeys = ['pp', 'cc', 'nn']\n self.geos = []\n geo = {'pp':0, 'cc':0, 'nn':0}\n \n for ret in self.rets: \n for nk, ok in zip(newkeys, oldkeys):\n geo[nk] = geo[nk] + (ret[ok].astype(np.float64) /(len(self.rets)))\n self.geos.append(geo)\n\n\n def init_mask(self):\n halfmask=np.zeros([self.L,self.L])\n fullmask=np.zeros([self.L,self.L])\n for i in range(self.L):\n for j in range(i+1,self.L):\n halfmask[i,j]=1\n fullmask[i,j]=1\n fullmask[j,i]=1\n self.halfmask=torch.DoubleTensor(halfmask) > 0.5\n self.fullmask=torch.DoubleTensor(fullmask) > 0.5\n self.clash_mask = torch.zeros([self.L,self.L,22,22])\n for i in range(self.L):\n for j in range(i+1,self.L):\n self.clash_mask[i,j]=1\n\n for i in range(self.L):\n self.clash_mask[i,i,:6,7:]=1\n\n for i in range(self.L-1):\n self.clash_mask[i,i+1,:,0]=0\n self.clash_mask[i,i+1,0,:]=0\n self.clash_mask[i,i+1,:,5]=0\n self.clash_mask[i,i+1,5,:]=0\n\n self.side_mask = rigid.side_mask(self.seq)\n self.side_mask = self.side_mask[:,None,:,None] * self.side_mask[None,:,None,:]\n self.clash_mask = (self.clash_mask > 0.5) * (self.side_mask > 0.5)\n\n self.geo_confimask_cc = []\n self.geo_confimask_pp = []\n self.geo_confimask_nn = []\n for geo in self.geos:\n confimask_cc = torch.DoubleTensor(geo['cc'][:,:,-1]) < 0.5\n confimask_pp = torch.DoubleTensor(geo['pp'][:,:,-1]) < 0.5\n confimask_nn = torch.DoubleTensor(geo['nn'][:,:,-1]) < 0.5\n self.geo_confimask_cc.append(confimask_cc)\n self.geo_confimask_pp.append(confimask_pp)\n self.geo_confimask_nn.append(confimask_nn)\n\n # Move masks and confimasks to the GPU/CPU device\n self.halfmask = self.halfmask.to(device)\n self.fullmask = self.fullmask.to(device)\n self.clash_mask = self.clash_mask.to(device)\n self.side_mask = self.side_mask.to(device)\n # geo_confimasks are lists\n self.geo_confimask_cc = [m.to(device) for m in self.geo_confimask_cc]\n self.geo_confimask_pp = [m.to(device) for m in self.geo_confimask_pp]\n self.geo_confimask_nn = [m.to(device) for m in self.geo_confimask_nn]\n\n\n def init_paras(self):\n self.geo_cc = []\n self.geo_pp = []\n self.geo_nn = []\n self.cs_coefs = {'cc': [], 'pp': [], 'nn': []}\n self.cs_knots = {'cc': [], 'pp': [], 'nn': []}\n for geo in self.geos:\n cc_cs, cc_decs = Cubic.dis_cubic(geo['cc'], 2, 40, 36)\n pp_cs, pp_decs = Cubic.dis_cubic(geo['pp'], 2, 40, 36)\n nn_cs, nn_decs = Cubic.dis_cubic(geo['nn'], 2, 40, 36)\n self.geo_cc.append([cc_cs, cc_decs])\n self.geo_pp.append([pp_cs, pp_decs])\n self.geo_nn.append([nn_cs, nn_decs])\n L = self.L\n cc_coefs_np = np.stack([[cc_cs[i,j].c for j in range(L)] for i in range(L)], axis=0)\n cc_knots_np = np.stack([[cc_cs[i,j].x for j in range(L)] for i in range(L)], axis=0)\n self.cs_coefs['cc'].append(torch.from_numpy(cc_coefs_np).to(device))\n self.cs_knots['cc'].append(torch.from_numpy(cc_knots_np).to(device))\n pp_coefs_np = np.stack([[pp_cs[i,j].c for j in range(L)] for i in range(L)], axis=0)\n pp_knots_np = np.stack([[pp_cs[i,j].x for j in range(L)] for i in range(L)], axis=0)\n self.cs_coefs['pp'].append(torch.from_numpy(pp_coefs_np).to(device))\n self.cs_knots['pp'].append(torch.from_numpy(pp_knots_np).to(device))\n nn_coefs_np = np.stack([[nn_cs[i,j].c for j in range(L)] for i in range(L)], axis=0)\n nn_knots_np = np.stack([[nn_cs[i,j].x for j in range(L)] for i in range(L)], axis=0)\n self.cs_coefs['nn'].append(torch.from_numpy(nn_coefs_np).to(device))\n self.cs_knots['nn'].append(torch.from_numpy(nn_knots_np).to(device))\n \n\n def _cubic_pair_energy(self, atom_map, geo_cs, geo_confimask, weight_key):\n \"\"\"General cubic-spline energy for CC/PP/NN pairs.\"\"\"\n min_dis, max_dis, bin_num = 2, 40, 36\n dev = atom_map.device\n upper_th = max_dis - ((max_dis - min_dis) / bin_num) * 0.5\n lower_th = 2.5\n total = torch.zeros((), device=dev, dtype=torch.double)\n spline_key = weight_key.split('_')[1]\n coeffs_list = self.cs_coefs[spline_key]\n knots_list = self.cs_knots[spline_key]\n for block_idx, mask_block in enumerate(geo_confimask):\n mask = (atom_map <= upper_th) & mask_block & self.fullmask & (atom_map >= lower_th)\n idx = mask.nonzero(as_tuple=True)\n if idx[0].numel() > 1:\n coef = coeffs_list[block_idx][idx]\n knots = knots_list[block_idx][idx]\n part1 = Potential.cubic_distance(atom_map[mask], coef, knots, min_dis, max_dis, bin_num).sum() * self.config[weight_key] * 0.5\n else:\n part1 = torch.zeros((), device=dev, dtype=torch.double)\n part2 = ((atom_map <= lower_th) & mask_block & self.fullmask).sum() * self.config[weight_key]\n total = total + part1 + part2\n return total\n\n # GPU-friendly torsion and angle energy helpers\n def _cubic_torsion_energy(self, atom_map, coef, x_vals, weight_key, num_bin):\n energy = Potential.cubic_torsion(atom_map, coef, x_vals, num_bin)\n return energy.sum() * self.config[weight_key]\n\n def _cubic_angle_energy(self, atom_map, coef, x_vals, weight_key, num_bin):\n energy = Potential.cubic_angle(atom_map, coef, x_vals, num_bin)\n return energy.sum() * self.config[weight_key]\n\n def compute_cc_energy(self, coor):\n atom_map = operations.pair_distance(coor[:,1], coor[:,1])\n return self._cubic_pair_energy(atom_map, self.geo_cc, self.geo_confimask_cc, 'weight_cc')\n\n def compute_pp_energy(self, coor):\n atom_map = operations.pair_distance(coor[:,0], coor[:,0])\n return self._cubic_pair_energy(atom_map, self.geo_pp, self.geo_confimask_pp, 'weight_pp')\n\n def compute_nn_energy(self, coor):\n atom_map = operations.pair_distance(coor[:,-1], coor[:,-1])\n return self._cubic_pair_energy(atom_map, self.geo_nn, self.geo_confimask_nn, 'weight_nn')\n\n def compute_pccp_energy(self, coor):\n # P-C-C-P dihedral energy on GPU\n p = coor[:, 0]\n c = coor[:, 1]\n dia = operations.dihedral(\n p[self.pccpi], c[self.pccpi], c[self.pccpj], p[self.pccpj]\n )\n return self._cubic_torsion_energy(dia, self.pccp_coe, self.pccp_x, 'weight_pccp', 36)\n\n def compute_cnnc_energy(self, coor):\n # C-N-N-C dihedral energy on GPU\n n = coor[:, -1]\n c = coor[:, 1]\n dia = operations.dihedral(\n c[self.cnnci], n[self.cnnci], n[self.cnncj], c[self.cnncj]\n )\n return self._cubic_torsion_energy(dia, self.cnnc_coe, self.cnnc_x, 'weight_cnnc', 36)\n\n def compute_pnnp_energy(self, coor):\n # P-N-N-P dihedral energy on GPU\n n = coor[:, -1]\n p = coor[:, 0]\n dia = operations.dihedral(\n p[self.pnnpi], n[self.pnnpi], n[self.pnnpj], p[self.pnnpj]\n )\n return self._cubic_torsion_energy(dia, self.pnnp_coe, self.pnnp_x, 'weight_pnnp', 36)\n\n def compute_pcc_energy(self, coor):\n # P-C-C angle energy on GPU\n p = coor[:, 1]\n c = coor[:, 2]\n ang = operations.angle(\n p[self.pcci], c[self.pcci], c[self.pccj]\n )\n return self._cubic_angle_energy(ang, self.pcc_coe, self.pcc_x, 'weight_pcc', 12)\n\n def compute_fape_energy(self,coor,ep=1e-3,epmax=20):\n energy= 0\n for tx in self.tx2ds:\n px_mean = coor[:,[1]]\n p_rot = operations.rigidFrom3Points(coor)\n p_tran = px_mean[:,0]\n pred_x2 = coor[:,None,:,:] - p_tran[None,:,None,:] # Lx Lrot N , 3\n pred_x2 = torch.einsum('ijnd,jde->ijne',pred_x2,p_rot.transpose(-1,-2)) # transpose should be equal to inverse\n errmap=torch.sqrt( ((pred_x2 - tx)**2).sum(dim=-1) + ep )\n energy = energy + torch.sum( torch.clamp(errmap,max=epmax) )\n return energy * self.config['weight_fape']\n\n def compute_bond_energy(self,coor,other_coor):\n # 3.87\n o3 = other_coor[:-1,-2]\n p = coor[1:,0]\n dis = (o3-p).norm(dim=-1)\n energy = ((dis-1.607)**2).sum()\n return energy * self.config['weight_bond']\n\n def tooth_func(self,errmap, ep = 0.05):\n return -1/(errmap/10+ep) + (1/ep)\n \n def reweight_func(self,ww):\n reweighting = torch.pow(ww,self.config['pair_weight_power'])\n reweighting[ww < self.config['pair_weight_min']] = 0\n return reweighting\n \n def compute_fape_energy_fromquat(self,x,coor,ep=1e-6,epmax=100):\n energy= 0\n p_rot,px_mean = a2b.Non2rot(x[:,:9],x.shape[0]),x[:,9:]\n pred_x2 = coor[:,None,:,:] - px_mean[None,:,None,:] # Lx Lrot N , 3\n pred_x2 = torch.einsum('ijnd,jde->ijne',pred_x2,p_rot.transpose(-1,-2)) # transpose should be equal to inverse\n\n for tx,weightplddt in zip(self.tx2ds,self.pair):\n tamplate_dist_map = torch.min( tx.norm(dim=-1), dim=2 )[0]\n errmap=torch.sqrt( ((pred_x2 - tx)**2).sum(dim=-1) + ep ) \n energy = energy + torch.sum( ( (torch.clamp(errmap,max=self.config['FAPE_max'])**self.config['pair_error_power']) * self.reweight_func(weightplddt[...,None]) )[tamplate_dist_map>self.config['pair_rest_min_dist']] )\n\n return energy * self.config['weight_fape']\n \n def compute_fape_energy_fromcoor(self,coor,ep=1e-6,epmax=100):\n energy= 0\n \n p_rot,px_mean = operations.Kabsch_rigid(self.basex,coor[:,0],coor[:,1],coor[:,2])\n pred_x2 = coor[:,None,:,:] - px_mean[None,:,None,:] # Lx Lrot N , 3\n pred_x2 = torch.einsum('ijnd,jde->ijne',pred_x2,p_rot.transpose(-1,-2)) # transpose should be equal to inverse\n \n for tx,weightplddt in zip(self.tx2ds,self.pair):\n tamplate_dist_map = torch.min( tx.norm(dim=-1), dim=2 )[0]\n errmap=torch.sqrt( ((pred_x2 - tx)**2).sum(dim=-1) + ep ) \n energy = energy + torch.sum( ( (torch.clamp(errmap,max=self.config['FAPE_max'])**self.config['pair_error_power']) * self.reweight_func(weightplddt[...,None]) )[tamplate_dist_map>self.config['pair_rest_min_dist']] )\n\n return energy * self.config['weight_fape']\n \n \n def energy(self, rama):\n coor = a2b.quat2b(self.basex, rama[:, 9:])\n other_coor = a2b.quat2b(self.otherx, rama[:, 9:])\n side_coor = a2b.quat2b(self.sidex, torch.cat([rama[:, :9], coor[:, -1]], dim=-1))\n\n E_cc = self.compute_cc_energy(coor) / len(self.geofiles) if self.config['weight_cc'] > 0 else 0\n E_pp = self.compute_pp_energy(coor) / len(self.geofiles) if self.config['weight_pp'] > 0 else 0\n E_nn = self.compute_nn_energy(coor) / len(self.geofiles) if self.config['weight_nn'] > 0 else 0\n E_pccp = self.compute_pccp_energy(coor) / len(self.geofiles) if self.config['weight_pccp'] > 0 else 0\n E_cnnc = self.compute_cnnc_energy(coor) / len(self.geofiles) if self.config['weight_cnnc'] > 0 else 0\n E_pnnp = self.compute_pnnp_energy(coor) / len(self.geofiles) if self.config['weight_pnnp'] > 0 else 0\n E_vdw = self.compute_full_clash(coor, other_coor, side_coor) if self.config['weight_vdw'] > 0 else 0\n E_fape = self.compute_fape_energy_fromquat(rama[:, 9:], coor) / len(self.geofiles) if self.config['weight_fape'] > 0 else 0\n E_bond = self.compute_bond_energy(coor, other_coor) if self.config['weight_bond'] > 0 else 0\n\n return E_vdw + E_fape + E_bond + E_pp + E_cc + E_nn + E_pccp + E_cnnc + E_pnnp\n\n\n def energy_from_coor(self, coor):\n E_cc = self.compute_cc_energy(coor) if self.config['weight_cc'] > 0 else 0\n E_pp = self.compute_pp_energy(coor) if self.config['weight_pp'] > 0 else 0\n E_nn = self.compute_nn_energy(coor) if self.config['weight_nn'] > 0 else 0\n E_fape = (self.compute_fape_energy_fromcoor(coor) / len(self.geofiles)) if self.config['weight_fape'] > 0 else 0\n print(E_fape, E_pp, E_cc, E_nn)\n return E_fape + E_pp + E_cc + E_nn \n\n def obj_func_grad_np(self,rama_):\n rama=torch.DoubleTensor(rama_)\n rama.requires_grad=True\n if rama.grad:\n rama.grad.zero_()\n f=self.energy(rama.view(self.L,21))*Scale_factor\n grad_value=autograd.grad(f,rama)[0]\n return grad_value.data.numpy().astype(np.float64)\n \n def obj_func_np(self,rama_):\n rama=torch.DoubleTensor(rama_)\n rama=rama.view(self.L,21)\n with torch.no_grad():\n f = self.energy(rama)*Scale_factor\n return f.item()\n\n def saveconfig(self,dict,confile):\n json_object = json.dumps(dict, indent = 4)\n wfile = open(confile,'w')\n wfile.write(json_object)\n wfile.close()\n \n def scoring(self):\n geoscale = self.config['geo_scale']\n self.config['weight_pp'] = geoscale * self.config['weight_pp']\n self.config['weight_cc'] = geoscale * self.config['weight_cc']\n self.config['weight_nn'] = geoscale * self.config['weight_nn']\n self.config['weight_pccp'] = geoscale * self.config['weight_pccp']\n self.config['weight_cnnc'] = geoscale * self.config['weight_cnnc']\n self.config['weight_pnnp'] = geoscale * self.config['weight_pnnp'] \n \n energy_dict = {}\n saveenergy_dict = {}\n \n with torch.no_grad():\n for retfile, tx in zip(self.geofiles, self.txs):\n one = self.energy_from_coor(tx)\n aaretfile = os.path.basename(retfile) \n energy_dict[aaretfile] = one.item()\n saveenergy_dict[retfile] = one.item()\n self.saveconfig(energy_dict, self.saveprefix)\n\n\n def foldning(self):\n minenergy=1e16\n count=0\n for tx in self.txs:\n count+=1\n \n minirama=None\n\n ilter = self.init_ret\n selected_ret = self.geofiles[ilter]\n try:\n rama=self.init_quat(ilter).data.numpy()\n self.config=readconfig(os.path.join(os.path.dirname(os.path.abspath(__file__)),'lib','vdw.json'))\n rama = fmin_l_bfgs_b(func=self.obj_func_np, x0=rama, fprime=self.obj_func_grad_np,iprint=10,maxfun=100)[0]\n rama = rama.flatten()\n except:\n rama=self.init_quat_safe(ilter).data.numpy()\n self.config=readconfig(os.path.join(os.path.dirname(os.path.abspath(__file__)),'lib','vdw.json'))\n rama = fmin_l_bfgs_b(func=self.obj_func_np, x0=rama, fprime=self.obj_func_grad_np,iprint=10,maxfun=100)[0]\n rama = rama.flatten()\n \n self.config=readconfig(self.foldconfig)\n geoscale = self.config['geo_scale']\n self.config['weight_pp'] =geoscale * self.config['weight_pp']\n self.config['weight_cc'] =geoscale * self.config['weight_cc']\n self.config['weight_nn'] =geoscale * self.config['weight_nn']\n self.config['weight_pccp'] =geoscale * self.config['weight_pccp']\n self.config['weight_cnnc'] =geoscale * self.config['weight_cnnc']\n self.config['weight_pnnp'] =geoscale * self.config['weight_pnnp']\n for i in range(3):\n line_min = lbfgs_rosetta.ArmijoLineMinimization(self.obj_func_np,self.obj_func_grad_np,True,len(rama),120)\n lbfgs_opt = lbfgs_rosetta.lbfgs(self.obj_func_np,self.obj_func_grad_np)\n rama=lbfgs_opt.run(rama,256,lbfgs_rosetta.absolute_converge_test,line_min,8000,self.obj_func_np,self.obj_func_grad_np,1e-9)\n newrama=rama+0.0\n newrama=torch.DoubleTensor(newrama) \n current_energy =self.obj_func_np(rama)\n\n if current_energy < minenergy:\n print(current_energy,minenergy)\n minenergy=current_energy\n self.outpdb(newrama,self.saveprefix+'.pdb',energystr=str(current_energy))\n\n\n def outpdb(self,rama,savefile,start=0,end=10000,energystr=''):\n coor_np=a2b.quat2b(self.basex,rama.view(self.L,21)[:,9:]).data.numpy()\n other_np=a2b.quat2b(self.otherx,rama.view(self.L,21)[:,9:]).data.numpy()\n shaped_rama=rama.view(self.L,21)\n coor = torch.FloatTensor(coor_np)\n side_coor_NP = a2b.quat2b(self.sidex,torch.cat([shaped_rama[:,:9],coor[:,-1]],dim=-1)).data.numpy()\n \n Atom_name=[' P ',\" C4'\",' N1 ']\n Other_Atom_name = [\" O5'\",\" C5'\",\" C3'\",\" O3'\",\" C1'\"]\n other_last_name = ['O',\"C\",\"C\",\"O\",\"C\"]\n\n side_atoms= [' N1 ',' C2 ',' O2 ',' N2 ',' N3 ',' N4 ',' C4 ',' O4 ',' C5 ',' C6 ',' O6 ',' N6 ',' N7 ',' N8 ',' N9 ']\n side_last_name = ['N', \"C\", \"O\", \"N\", \"N\", 'N', 'C', 'O', 'C', 'C', 'O', 'N', 'N', 'N','N']\n\n base_dict = rigid.base_table()\n \n last_name=['P','C','N']\n wstr=[f'REMARK {str(energystr)}']\n templet='%6s%5d %4s %3s %1s%4d %8.3f%8.3f%8.3f%6.2f%6.2f %2s%2s'\n count=1\n for i in range(self.L):\n if self.seq[i] in ['a','g','A','G']:\n Atom_name = [' P ',\" C4'\",' N9 ']\n\n elif self.seq[i] in ['c','u','C','U']:\n Atom_name = [' P ',\" C4'\",' N1 ']\n for j in range(coor_np.shape[1]):\n outs=('ATOM ',count,Atom_name[j],self.seq[i],'A',i+1,coor_np[i][j][0],coor_np[i][j][1],coor_np[i][j][2],0,0,last_name[j],'')\n if i>=start-1 and i < end:\n wstr.append(templet % outs)\n count+=1\n\n for j in range(other_np.shape[1]):\n outs=('ATOM ',count,Other_Atom_name[j],self.seq[i],'A',i+1,other_np[i][j][0],other_np[i][j][1],other_np[i][j][2],0,0,other_last_name[j],'')\n if i>=start-1 and i < end:\n wstr.append(templet % outs)\n count+=1\n \n wstr='\\n'.join(wstr)\n wfile=open(savefile,'w')\n wfile.write(wstr)\n wfile.close()\n \n \n def outpdb_coor(self,coor_np,savefile,start=0,end=1000,energystr=''):\n Atom_name=[' P ',\" C4'\",' N1 ']\n last_name=['P','C','N']\n wstr=[f'REMARK {str(energystr)}']\n templet='%6s%5d %4s %3s %1s%4d %8.3f%8.3f%8.3f%6.2f%6.2f %2s%2s'\n count=1\n for i in range(self.L):\n if self.seq[i] in ['a','g','A','G']:\n Atom_name = [' P ',\" C4'\",' N9 ']\n\n elif self.seq[i] in ['c','u','C','U']:\n Atom_name = [' P ',\" C4'\",' N1 ']\n \n for j in range(coor_np.shape[1]):\n outs=('ATOM ',count,Atom_name[j],self.seq[i],'A',i+1,coor_np[i][j][0],coor_np[i][j][1],coor_np[i][j][2],0,0,last_name[j],'')\n if i>=start-1 and i < end:\n wstr.append(templet % outs)\n count+=1\n \n wstr='\\n'.join(wstr)\n wfile=open(savefile,'w')\n wfile.write(wstr)\n wfile.close()\n\n\nif __name__ == '__main__': \n\n fastafile = sys.argv[1]\n foldconfig = sys.argv[2]\n save_prefix = sys.argv[3]\n retfiles = sys.argv[4:]\n\n save_parent_dir = os.path.dirname(save_prefix)\n if not os.path.isdir(save_parent_dir):\n os.makedirs(save_parent_dir)\n\n retfiles.sort()\n print(retfiles)\n\n stru = Structure(fastafile, retfiles, foldconfig, save_prefix) \n stru.scoring()","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:42.533007Z","iopub.status.busy":"2026-03-27T18:13:42.532327Z","iopub.status.idle":"2026-03-27T18:13:42.547858Z","shell.execute_reply":"2026-03-27T18:13:42.546924Z"},"papermill":{"duration":0.041616,"end_time":"2026-03-27T18:13:42.549693","exception":false,"start_time":"2026-03-27T18:13:42.508077","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile /kaggle/working/DRfold2/PotentialFold/Cubic.py\nimport numpy as np \nfrom scipy.interpolate import CubicSpline,UnivariateSpline\nimport os\nfrom torch.autograd import Function\nimport torch\nimport math\n\ndef fit_dis_cubic(dis_matrix,min_dis,max_dis,num_bin):\n # convert torch Tensor on GPU to numpy array for SciPy\n if isinstance(dis_matrix, torch.Tensor):\n dis_matrix = dis_matrix.detach().cpu().numpy()\n dis_region=np.zeros(num_bin)\n for i in range(num_bin):\n dis_region[i]=min_dis+(i+0.5)*(max_dis-min_dis)*1.0/num_bin\n L=dis_matrix.shape[0]\n csnp=[]\n decsnp=[]\n for i in range(L):\n css=[]\n decss=[]\n for j in range(L):\n y=-np.log( (dis_matrix[i,j,1:-1]+1e-8) / (dis_matrix[i,j,[-2]]+1e-8) )\n x=dis_region\n x[0]=-0.0001\n y[0]= max(10,y[1]+4)\n cs= CubicSpline(x,y)\n decs=cs.derivative()\n css.append(cs)\n decss.append(decs)\n csnp.append(css)\n decsnp.append(decss)\n return np.array(csnp),np.array(decsnp)\n\ndef dis_cubic(out,min_dis,max_dis,num_bin):\n print('fitting cubic distance')\n cs,decs=fit_dis_cubic(out,min_dis,max_dis,num_bin)\n return cs,decs\n\n\n\ndef cubic_matrix_torsion(dis_matrix,min_dis,max_dis,num_bin):\n dis_region=np.zeros(num_bin)\n bin_size=(max_dis-min_dis)/num_bin\n for i in range(num_bin):\n dis_region[i]=min_dis+(i+0.5)*(max_dis-min_dis)*1.0/num_bin\n L=dis_matrix.shape[0]\n csnp=[]\n decsnp=[]\n for i in range(L):\n css=[]\n decss=[]\n for j in range(L):\n y=-np.log( dis_matrix[i,j,:-1]+1e-8 )\n x=dis_region\n x=np.append(x,x[-1]+bin_size)\n y=np.append(y,y[0])\n cs= CubicSpline(x,y,bc_type='periodic')\n decs=cs.derivative()\n css.append(cs)\n decss.append(decs)\n csnp.append(css)\n decsnp.append(decss)\n return np.array(csnp),np.array(decsnp)\ndef torsion_cubic(out,min_dis,max_dis,num_bin):\n print('fitting cubic')\n cs,decs=cubic_matrix_torsion(out,min_dis,max_dis,num_bin)\n return cs,decs\n\ndef cubic_matrix_angle(dis_matrix,min_dis,max_dis,num_bin): # 0 - np.pi 12\n dis_region=np.zeros(num_bin)\n bin_size=(max_dis-min_dis)/num_bin\n for i in range(num_bin):\n dis_region[i]=min_dis+(i+0.5)*(max_dis-min_dis)*1.0/num_bin\n L=dis_matrix.shape[0]\n csnp=[]\n decsnp=[]\n for i in range(L):\n css=[]\n decss=[]\n for j in range(L):\n y=-np.log( dis_matrix[i,j,:-1]+1e-8 )\n x=dis_region\n\n x=np.concatenate([[x[0]-bin_size*3,x[0]-bin_size*2,x[0]-bin_size], x,[x[-1]+bin_size,x[-1]+bin_size*2,x[-1]+bin_size*3] ])\n y=np.concatenate([ [y[2],y[1],y[0]],y,[y[-1],y[-2],y[-3]] ])\n\n cs= CubicSpline(x,y)\n decs=cs.derivative()\n\n css.append(cs)\n decss.append(decs)\n csnp.append(css)\n decsnp.append(decss)\n\n return np.array(csnp),np.array(decsnp)\ndef angle_cubic(out,min_dis,max_dis,num_bin):\n\n print('fitting angle cubic')\n cs,decs=cubic_matrix_angle(out,min_dis,max_dis,num_bin)\n\n return cs,decs","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:42.595972Z","iopub.status.busy":"2026-03-27T18:13:42.595585Z","iopub.status.idle":"2026-03-27T18:13:42.60275Z","shell.execute_reply":"2026-03-27T18:13:42.601894Z"},"papermill":{"duration":0.032175,"end_time":"2026-03-27T18:13:42.604628","exception":false,"start_time":"2026-03-27T18:13:42.572453","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Replace 'your-dataset-name' with the actual folder name from Step 1\n!pip install --no-index --find-links=/kaggle/input/datasets/khushvaghela/gcosrna/gcosrna_wheels/gcosrna_wheels draw_rna","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:42.649404Z","iopub.status.busy":"2026-03-27T18:13:42.649066Z","iopub.status.idle":"2026-03-27T18:13:47.069014Z","shell.execute_reply":"2026-03-27T18:13:47.067683Z"},"papermill":{"duration":4.44498,"end_time":"2026-03-27T18:13:47.071284","exception":false,"start_time":"2026-03-27T18:13:42.626304","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import os\n# import sys\n\n# # 1. Define where your wheels are (UPDATE THIS PATH with your actual dataset path)\n# OFFLINE_DATASET_PATH = \"/kaggle/input/datasets/khushvaghela/gcosrna/python_packages/kaggle/working/python_packages\" # <--- Fix this!\n\n# # 2. Define a target folder to install them into so they don't conflict with system libs\n# LIB_PATH = \"/kaggle/working/custom_lib\"\n# os.makedirs(LIB_PATH, exist_ok=True)\n\n# # 3. Install specifically to that folder\n# # We use scipy==1.11.4 because 1.13+ is what's causing the 'numpy.strings' error\n# !pip install --no-index --find-links={OFFLINE_DATASET_PATH} \\\n# --target={LIB_PATH} \\\n# \"numpy==1.26.4\" \\\n# \"pandas==2.1.4\" \\\n# \"scikit-learn==1.3.2\" \\\n# \"scipy==1.11.4\"\n\n# # 4. FORCE Python to use these versions by putting this folder at the START of the path\n# sys.path.insert(0, LIB_PATH)\n\n# print(\"โœ… Custom libraries installed and prioritized.\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:47.119183Z","iopub.status.busy":"2026-03-27T18:13:47.118672Z","iopub.status.idle":"2026-03-27T18:13:47.123796Z","shell.execute_reply":"2026-03-27T18:13:47.122881Z"},"papermill":{"duration":0.031575,"end_time":"2026-03-27T18:13:47.12564","exception":false,"start_time":"2026-03-27T18:13:47.094065","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport shutil\nfrom pathlib import Path\n\n# Paths - Adjust 'simrna-standalone' to your actual dataset name\nSIMRNA_SRC = Path(\"/kaggle/input/datasets/khushvaghela/simrna/SimRNA_64bitIntel_Linux\") \nSIMRNA_PKG = Path(\"/kaggle/working/SimRNA_pkg\")\n\ndef setup_simrna():\n if SIMRNA_PKG.exists():\n shutil.rmtree(SIMRNA_PKG)\n \n # Copy the package to a writable directory\n shutil.copytree(SIMRNA_SRC, SIMRNA_PKG)\n \n # Ensure the 'data' folder is directly inside the working directory [cite: 13]\n # The manual states SimRNA will crash without this [cite: 18, 19]\n data_dir = SIMRNA_PKG / \"data\"\n if not data_dir.exists():\n print(\"โŒ Error: 'data' folder not found in the package!\")\n else:\n print(\"โœ… SimRNA data folder (with .hist files) is in place[cite: 12].\")\n\n # Set execution permissions for the binary\n simrna_bin = SIMRNA_PKG / \"SimRNA\"\n if simrna_bin.exists():\n os.chmod(simrna_bin, 0o755)\n print(\"โœ… SimRNA binary is executable.\")\n\nsetup_simrna()\n\nimport os\nimport shutil\nfrom pathlib import Path\n\ndef batch_repair_simrna_data(data_path):\n data_dir = Path(data_path)\n if not data_dir.exists():\n print(f\"โŒ Error: Directory {data_dir} not found.\")\n return\n\n hist_files = list(data_dir.glob(\"*.hist\"))\n repaired_count = 0\n \n print(f\"๐Ÿ” Scanning {len(hist_files)} .hist files for corruption...\")\n\n for target_file in hist_files:\n try:\n with open(target_file, 'r') as f:\n first_line = f.readline().strip()\n # Check if first line looks like a filename and line 2 is empty\n second_line = f.readline().strip()\n \n # If line 1 is another file in the same folder and line 2 is empty...\n source_file = data_dir / first_line\n if not second_line and source_file.exists() and source_file != target_file:\n print(f\"๐Ÿ›  Repairing: {target_file.name} (Placeholder for {first_line})\")\n \n # Backup and replace\n shutil.copy(source_file, f\"{target_file}.bak\")\n shutil.copy(source_file, target_file)\n repaired_count += 1\n \n except Exception as e:\n print(f\"โš ๏ธ Could not check {target_file.name}: {e}\")\n\n print(f\"\\nโœ… Finished. Repaired {repaired_count} files.\")\n\n# Execute repair on your specific folder\nbatch_repair_simrna_data(SIMRNA_PKG / \"data\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:47.172627Z","iopub.status.busy":"2026-03-27T18:13:47.17197Z","iopub.status.idle":"2026-03-27T18:13:49.203227Z","shell.execute_reply":"2026-03-27T18:13:49.202145Z"},"papermill":{"duration":2.057174,"end_time":"2026-03-27T18:13:49.205292","exception":false,"start_time":"2026-03-27T18:13:47.148118","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# gCoSRNA coaxial-stacking โ€” IN-PROCESS (returns residue indices, not just text)\n# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# Why in-process instead of the old `subprocess` call?\n# * The subprocess only printed lines like \"Multibranch Stem 4 and Stem 1 ...\".\n# Those \"Stem N\" are 1-based indices LOCAL to one multibranch junction's stem\n# list (a Python set -> order not even guaranteed) and are AMBIGUOUS when a\n# structure has more than one multibranch loop. You cannot reliably turn them\n# back into residues from stdout.\n# * Calling gCoSRNA's own functions here lets us read the actual stem objects,\n# so every predicted stack maps to exact 0-based residue indices + the pivot\n# loop between the two stems โ€” which is what the geometry step needs.\nimport os, sys, importlib.util, warnings\nimport numpy as np\n\n# (a) scikit-learn 1.3.2-pickle compatibility under the image's sklearn >=1.4 (1.6.1),\n# applied in THIS kernel (gCoSRNA_model.pkl is loaded here via joblib).\nimport sklearn.tree._classes as _t\nif not hasattr(_t.BaseDecisionTree, \"monotonic_cst\"):\n _t.BaseDecisionTree.monotonic_cst = None\nfor _cn in (\"DecisionTreeClassifier\", \"DecisionTreeRegressor\",\n \"ExtraTreeClassifier\", \"ExtraTreeRegressor\"):\n _c = getattr(_t, _cn, None)\n if _c is not None and \"monotonic_cst\" not in _c.__dict__:\n _c.monotonic_cst = None\ntry:\n from sklearn.exceptions import InconsistentVersionWarning\n warnings.filterwarnings(\"ignore\", category=InconsistentVersionWarning)\nexcept Exception:\n pass\n\nGCOSRNA_SRC = \"/kaggle/input/datasets/khushvaghela/gcosrna/gCoSRNA-main/gCoSRNA-main/gCoSRNA_predict.py\"\nGCOSRNA_LIB = \"/kaggle/working/gcosrna_lib.py\"\nGCOSRNA_MODEL = \"/kaggle/input/datasets/khushvaghela/gcosrna/gCoSRNA_model.pkl\"\n\n\ndef _build_gcosrna_lib(src=GCOSRNA_SRC, dst=GCOSRNA_LIB):\n \"\"\"Make an import-safe copy of gCoSRNA_predict.py:\n (1) stub the `from draw_rna...` import so no matplotlib/draw_rna is needed,\n (2) strip the module-level Example-4 demo block that otherwise executes\n (draw + a wall of prints) on every import.\"\"\"\n code = open(src, encoding=\"utf-8\", errors=\"replace\").read()\n code = code.replace(\"from draw_rna.ipynb_draw import draw_struct\",\n \"def draw_struct(*a, **k):\\n return None\")\n lines = code.split(\"\\n\")\n def_idx = next(i for i, l in enumerate(lines)\n if l.lstrip().startswith(\"def predict_coaxial_stacking\"))\n stubs = [i for i, l in enumerate(lines[:def_idx]) if l.strip() == \"def draw_struct(*a, **k):\"]\n demo_start = stubs[-1] if len(stubs) >= 2 else def_idx # 2nd stub == start of demo block\n kept = lines[:demo_start] + [\"# [demo block stripped for safe import]\"] + lines[def_idx:]\n open(dst, \"w\").write(\"\\n\".join(kept))\n return dst\n\n\n_GC = None\ndef _gc():\n \"\"\"Lazily load the sanitized gCoSRNA module.\"\"\"\n global _GC\n if _GC is None:\n path = _build_gcosrna_lib()\n spec = importlib.util.spec_from_file_location(\"gcosrna_lib\", path)\n mod = importlib.util.module_from_spec(spec)\n sys.modules[\"gcosrna_lib\"] = mod\n spec.loader.exec_module(mod)\n _GC = mod\n return _GC\n\n\ndef _branch_for_pair(a, b, ps):\n \"\"\"A returned cycle pair (a,b) <-> the one branch whose two flanking stems are\n exactly {a-1, b-1}. That branch carries the pivot (between-stems) loop.\"\"\"\n want = {a - 1, b - 1}\n for p in ps:\n if set(p[\"stems\"]) == want:\n return p\n return None\n\n\ndef get_coaxial_stem_pairs(sequence, dot_bracket, model_path=GCOSRNA_MODEL):\n \"\"\"Run gCoSRNA in-process and return a list of predicted coaxial stacks:\n {\n 'stemA', 'stemB' : lists of (i, j) base-pair tuples (0-based),\n 'stemA_res', 'stemB_res': sorted residue indices of each whole stem,\n 'pivot' : (start, end) of the junction loop between them,\n 'prob' : coaxial probability,\n 'label' : ('Stem a', 'Stem b') (gCoSRNA's own labels)\n }\n All indices are 0-based positions in `sequence` / `dot_bracket`, i.e. they\n line up 1:1 with the rows of the 3D `coords` array.\"\"\"\n G = _gc()\n out = []\n try:\n pairs = G.parse_dot_bracket(dot_bracket)\n stems = G.extract_stems(pairs)\n loop_info = G.identify_loops(sequence, dot_bracket, stems)\n mbs = G.extract_substructures(sequence, dot_bracket, stems, loop_info, \"multibranch\")\n except Exception as e:\n print(f\"โš ๏ธ gCoSRNA parse failed: {e}\")\n return out\n\n for sub in mbs:\n try:\n ps = G.extract_pseudo_two_way_substructures_v2(sub, sequence, dot_bracket)\n prob_df = G.predict_coaxial_probabilities(ps, model_path)\n probs = prob_df[\"coaxial_prob\"].tolist()\n _, cyc = G.assign_final_coaxial_stacking(ps, probs)\n except Exception as e:\n print(f\"โš ๏ธ gCoSRNA scoring failed on a junction: {e}\")\n continue\n if cyc == \"none\" or not cyc:\n continue\n for (a, b) in cyc:\n p = _branch_for_pair(a, b, ps)\n if p is None:\n continue\n stemA, stemB = sub[\"stems\"][a - 1], sub[\"stems\"][b - 1]\n resA = sorted({x for ij in stemA for x in ij})\n resB = sorted({x for ij in stemB for x in ij})\n piv = p[\"between_loops\"][0]\n try:\n prob = float(probs[ps.index(p)])\n except Exception:\n prob = None\n out.append({\n \"stemA\": stemA, \"stemB\": stemB,\n \"stemA_res\": resA, \"stemB_res\": resB,\n \"pivot\": (piv[0], piv[1]),\n \"prob\": prob,\n \"label\": (f\"Stem {a}\", f\"Stem {b}\"),\n })\n print(f\" ๐Ÿ”— Coaxial: Stem {a} โ†” Stem {b} \"\n f\"(res {resA[0]}-{resA[-1]} โ†” {resB[0]}-{resB[-1]}, pivot {piv[0]}-{piv[1]})\")\n return out\n\n\n# Backward-compatible thin wrapper: same return shape as the old subprocess version\n# (a list of ('Stem a','Stem b') tuples). `script_path` is accepted but ignored now.\ndef run_gcosrna_alignment(sequence, structure, model_path=GCOSRNA_MODEL, script_path=None):\n print(\"๐Ÿงฌ Running gCoSRNA Junction Alignment (in-process)...\")\n return [r[\"label\"] for r in get_coaxial_stem_pairs(sequence, structure, model_path)]\n","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:49.252554Z","iopub.status.busy":"2026-03-27T18:13:49.251687Z","iopub.status.idle":"2026-03-27T18:13:49.258399Z","shell.execute_reply":"2026-03-27T18:13:49.257687Z"},"papermill":{"duration":0.032019,"end_time":"2026-03-27T18:13:49.260124","exception":false,"start_time":"2026-03-27T18:13:49.228105","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\n\n# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# Coaxial-stacking GEOMETRY\n# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# Given a predicted coaxial stack between two stems, rotate the SMALLER (\"mobile\")\n# helix arm so its axis becomes collinear with the other (\"reference\") helix,\n# pivoting at the junction. We move the WHOLE arm the mobile stem subtends\n# (range [outer_i .. outer_j]), not just the stem's own paired residues โ€” moving\n# only the stem would tear the backbone apart.\n#\n# โš ๏ธ This is a HEURISTIC post-processing step. Rigidly re-orienting an arm of an\n# already-predicted 3D model can introduce clashes and may HURT the score as\n# easily as help it. It is therefore: (a) gated by ENABLE_COAXIAL_ALIGNMENT,\n# (b) guarded by a clash check that rejects a rotation bringing the moved arm\n# within COAXIAL_CLASH_THRESHOLD ร… of the rest of the structure, and\n# (c) fully wrapped in try/except by the caller. A/B test it against your\n# leaderboard score before trusting it.\nENABLE_COAXIAL_ALIGNMENT = True # master switch โ€” set False to keep raw 3D coords\nCOAXIAL_CLASH_THRESHOLD = 3.0 # ร…; reject rotations that clash this tightly\nCOAXIAL_ANTIPARALLEL = True # stacked helices extend opposite ways from junction\n\n\ndef _mid(c, i, j):\n return (c[i] + c[j]) / 2.0\n\ndef _stem_axis(c, stem):\n \"\"\"Helix axis: outer base pair (at the junction) -> innermost base pair.\"\"\"\n o_i, o_j = stem[0]\n n_i, n_j = stem[-1]\n return _mid(c, n_i, n_j) - _mid(c, o_i, o_j)\n\ndef _enclosure(stem):\n \"\"\"All residue indices the stem subtends (its whole arm).\"\"\"\n return set(range(stem[0][0], stem[0][1] + 1))\n\ndef _rodrigues(v_from, v_to):\n a = v_from / (np.linalg.norm(v_from) + 1e-9)\n b = v_to / (np.linalg.norm(v_to) + 1e-9)\n axis = np.cross(a, b)\n s = np.linalg.norm(axis)\n if s < 1e-8:\n return np.eye(3)\n axis /= s\n ang = np.arccos(np.clip(np.dot(a, b), -1.0, 1.0))\n K = np.array([[0, -axis[2], axis[1]],\n [axis[2], 0, -axis[0]],\n [-axis[1], axis[0], 0]])\n return np.eye(3) + np.sin(ang) * K + (1 - np.cos(ang)) * (K @ K)\n\n\ndef align_stems_coaxial(coords, stemA, stemB,\n antiparallel=COAXIAL_ANTIPARALLEL,\n clash=COAXIAL_CLASH_THRESHOLD):\n \"\"\"Make stemA & stemB collinear by rotating the smaller arm about the junction.\n Returns (new_coords, info).\"\"\"\n coords = np.asarray(coords, dtype=float)\n eA, eB = _enclosure(stemA), _enclosure(stemB)\n # mobile = smaller arm (the nested child); reference = the larger (parent/longer sibling)\n if len(eA) <= len(eB):\n mob, ref, emob = stemA, stemB, eA\n else:\n mob, ref, emob = stemB, stemA, eB\n\n a_ref = _stem_axis(coords, ref)\n a_mob = _stem_axis(coords, mob)\n target = -a_ref if antiparallel else a_ref\n R = _rodrigues(a_mob, target)\n pivot = _mid(coords, mob[0][0], mob[0][1]) # mobile arm's junction-proximal base\n\n new = coords.copy()\n mob_idx = sorted(emob)\n for idx in mob_idx:\n new[idx] = R @ (coords[idx] - pivot) + pivot\n\n # clash guard: nearest approach between moved atoms and everything that stayed put\n fixed = [i for i in range(len(coords)) if i not in emob]\n if fixed:\n M = new[mob_idx]\n Fx = new[fixed]\n dmin = float(np.sqrt(((M[:, None, :] - Fx[None, :, :]) ** 2).sum(-1)).min())\n if dmin < clash:\n return coords, {\"applied\": False, \"reason\": f\"min_clash {dmin:.2f}ร… < {clash}ร…\",\n \"mobile_arm\": mob_idx}\n else:\n dmin = None\n return new, {\"applied\": True, \"mobile_arm\": mob_idx, \"min_clash\": dmin}\n\n\ndef apply_coaxial_alignment(coords, sequence, coaxial_pairs,\n antiparallel=COAXIAL_ANTIPARALLEL,\n clash=COAXIAL_CLASH_THRESHOLD):\n \"\"\"Apply each predicted coaxial stack to `coords` (in order), with guards.\"\"\"\n if not ENABLE_COAXIAL_ALIGNMENT or not coaxial_pairs:\n return coords\n coords = np.asarray(coords, dtype=float).copy()\n if coords.shape[0] != len(sequence):\n print(f\" โš ๏ธ coords/sequence length mismatch \"\n f\"({coords.shape[0]} vs {len(sequence)}); skipping alignment.\")\n return coords\n for pr in coaxial_pairs:\n try:\n coords, info = align_stems_coaxial(coords, pr[\"stemA\"], pr[\"stemB\"],\n antiparallel=antiparallel, clash=clash)\n if info[\"applied\"]:\n print(f\" โœ… stacked {pr['label'][0]}โ†”{pr['label'][1]} \"\n f\"(moved {len(info['mobile_arm'])} res, min_clash=\"\n f\"{info['min_clash']:.2f}ร…)\" if info['min_clash'] is not None else \"\")\n else:\n print(f\" โญ๏ธ skipped {pr['label'][0]}โ†”{pr['label'][1]} ({info['reason']})\")\n except Exception as e:\n print(f\" โš ๏ธ alignment error on {pr.get('label')}: {e}\")\n return coords\n\n\n# โ”€โ”€ kept for backwards-compatibility (the ORIGINAL stem-only rotator) โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# NOTE: this one rotates only the stem's own residues (tears the backbone) and\n# took a different signature; it is no longer used by the pipeline. Prefer\n# apply_coaxial_alignment() above.\ndef get_unit_vector(coords, indices):\n \"\"\"Calculates the normalized direction vector of a stem.\"\"\"\n vec = coords[indices[-1]] - coords[indices[0]]\n norm = np.linalg.norm(vec)\n return vec / norm if norm > 1e-6 else vec\n","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:49.307291Z","iopub.status.busy":"2026-03-27T18:13:49.306382Z","iopub.status.idle":"2026-03-27T18:13:49.316036Z","shell.execute_reply":"2026-03-27T18:13:49.315097Z"},"papermill":{"duration":0.035257,"end_time":"2026-03-27T18:13:49.317866","exception":false,"start_time":"2026-03-27T18:13:49.282609","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import subprocess\nimport os\nimport glob\nfrom pathlib import Path\n\n# --- CONFIGURATION ---\nSIMRNA_PKG = Path(\"/kaggle/working/SimRNA_pkg\")\nSIMRNA_BIN = SIMRNA_PKG / \"SimRNA\"\n\ndef run_simrna_settling(target_id, sequence, coords, steps=1000):\n \"\"\"\n Phase 4.2: Relaxing the 3D draft using SimRNA.\n \"\"\"\n # 1. Use absolute paths to avoid 'cwd' confusion\n work_dir = Path(f\"/kaggle/working/simrna_{target_id}\").resolve()\n work_dir.mkdir(parents=True, exist_ok=True)\n \n seq_file = work_dir / f\"{target_id}.seq\"\n start_pdb = work_dir / \"start_input.pdb\"\n # We use a simple prefix for the -o flag\n output_prefix = \"settled\" \n \n with open(seq_file, \"w\") as f:\n f.write(sequence)\n \n # Save the DRFold2/gCoSRNA draft as the starting point\n save_coords_to_pdb(sequence, coords, target_id, str(start_pdb))\n \n # 2. Build the command\n # We run the command INSIDE the work_dir so all outputs land there\n cmd = [\n str(SIMRNA_BIN.resolve()),\n \"-s\", str(seq_file),\n \"-P\", str(start_pdb),\n \"-n\", str(steps),\n \"-o\", output_prefix \n ]\n \n print(f\"๐Ÿš€ SimRNA: Relaxing {target_id} in {work_dir}...\")\n os.chmod(SIMRNA_BIN, 0o755)\n \n # 3. Execute\n process = subprocess.Popen(\n cmd,\n cwd=str(work_dir), # Change CWD to the work_dir so files land here!\n stdout=subprocess.PIPE,\n stderr=subprocess.STDOUT,\n text=True,\n bufsize=1\n )\n\n # SimRNA needs to find its 'data' folder, so we link it if it's missing\n data_link = work_dir / \"data\"\n if not data_link.exists():\n os.symlink(SIMRNA_PKG / \"data\", data_link)\n\n for line in iter(process.stdout.readline, ''):\n if \"Energy\" in line:\n print(f\" [SimRNA]: {line.strip()}\", flush=True)\n \n process.wait()\n\n # 4. Find the Output PDB\n # Based on your logs, it will be named 'settled-000001.pdb'\n expected_pdb = work_dir / \"settled-000001.pdb\"\n \n if expected_pdb.exists():\n print(f\"โœ… Success! Found relaxed PDB: {expected_pdb.name}\")\n return extract_coords_from_pdb(str(expected_pdb), len(sequence))\n \n # 5. Ultimate Fallback: Search for ANY new PDB in the folder\n all_pdbs = list(work_dir.glob(\"*.pdb\"))\n for p in all_pdbs:\n if p.name != \"start_input.pdb\":\n print(f\"๐Ÿ” Found auto-named PDB: {p.name}\")\n return extract_coords_from_pdb(str(p), len(sequence))\n\n print(f\"โš ๏ธ No output PDB found in {work_dir}. Returning original coordinates.\")\n return coords","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:49.36568Z","iopub.status.busy":"2026-03-27T18:13:49.364762Z","iopub.status.idle":"2026-03-27T18:13:49.375798Z","shell.execute_reply":"2026-03-27T18:13:49.374869Z"},"papermill":{"duration":0.036591,"end_time":"2026-03-27T18:13:49.377593","exception":false,"start_time":"2026-03-27T18:13:49.341002","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\n\ndef kabsch_alignment(P, Q):\n \"\"\"\n Aligns P (template) to Q (target segment) using the Kabsch algorithm.\n P, Q: NumPy arrays of shape [N, 3].\n Returns the rotated and translated P.\n \"\"\"\n # 1. Center both sets of coordinates\n centroid_P = np.mean(P, axis=0)\n centroid_Q = np.mean(Q, axis=0)\n P_centered = P - centroid_P\n Q_centered = Q - centroid_Q\n\n # 2. Computation of the covariance matrix\n H = P_centered.T @ Q_centered\n\n # 3. Singular Value Decomposition (SVD)\n U, S, Vt = np.linalg.svd(H)\n V = Vt.T\n \n # 4. Correct for reflection (ensure a right-handed coordinate system)\n d = np.linalg.det(V @ U.T)\n step = np.eye(3)\n if d < 0:\n step[2, 2] = -1\n\n # 5. Optimal Rotation Matrix\n R = V @ step @ U.T\n\n # 6. Transform P to match Q's pose\n return (R @ P_centered.T).T + centroid_Q\n\ndef graft_functional_templates(coords, sequence, train_seqs, train_coords_dict, min_len=15):\n \"\"\"\n Phase 3: Finds exact structural motifs in training data and grafts \n them onto the predicted skeleton.\n \"\"\"\n new_coords = coords.copy()\n grafted_regions = [] # Track (start, end) to avoid overlapping grafts\n\n # Sort training sequences by length (longest first) \n # This ensures we graft the largest possible functional motifs first\n sorted_templates = train_seqs.copy()\n sorted_templates['len'] = sorted_templates['sequence'].str.len()\n sorted_templates = sorted_templates.sort_values('len', ascending=False)\n\n print(f\"๐Ÿงฉ Scanning for structural motifs (threshold: {min_len}nt)...\")\n\n for _, row in sorted_templates.iterrows():\n t_id, t_seq = row['target_id'], row['sequence']\n \n if len(t_seq) < min_len:\n continue\n \n # 1. Search for sequence match\n start_idx = sequence.find(t_seq)\n if start_idx != -1:\n end_idx = start_idx + len(t_seq)\n \n # 2. Safety Check: Avoid overlapping grafts\n if any(max(start_idx, g_start) < min(end_idx, g_end) for g_start, g_end in grafted_regions):\n continue\n \n # 3. Retrieve Template Coordinates\n t_coords = train_coords_dict.get(t_id)\n if t_coords is not None and len(t_coords) == len(t_seq):\n # 4. Perform Kabsch Alignment\n # We align the experimental template to the predicted 'rough' pose\n target_segment = new_coords[start_idx:end_idx]\n \n try:\n print(f\" โ†ณ Grafting {t_id} at residue {start_idx+1} (Length: {len(t_seq)})\")\n aligned_motif = kabsch_alignment(t_coords, target_segment)\n \n # 5. Overwrite coordinates\n new_coords[start_idx:end_idx] = aligned_motif\n grafted_regions.append((start_idx, end_idx))\n except Exception as e:\n print(f\" โš ๏ธ Alignment failed for {t_id}: {e}\")\n \n return new_coords","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:49.422985Z","iopub.status.busy":"2026-03-27T18:13:49.422589Z","iopub.status.idle":"2026-03-27T18:13:49.433726Z","shell.execute_reply":"2026-03-27T18:13:49.432767Z"},"papermill":{"duration":0.036065,"end_time":"2026-03-27T18:13:49.435538","exception":false,"start_time":"2026-03-27T18:13:49.399473","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport subprocess\nimport tempfile\nimport numpy as np\nimport Bio.PDB\nimport torch\nimport shutil\nfrom pathlib import Path\nfrom Bio.Align import PairwiseAligner\n\n# --- CONFIGURATION ---\nWORKING_DIR = Path(\"/kaggle/working\")\nDRFOLD_SCRIPT_DIR = Path(\"/kaggle/working/DRfold2\")\n\ndef run_drfold2_instance(sequence, consensus_ss, target_id, original_target_id, msa_dir, seed=42):\n msa_path = os.path.join(msa_dir, f\"{original_target_id}.MSA.fasta\")\n fasta_path = WORKING_DIR / f\"{target_id}.fasta\"\n output_dir = WORKING_DIR / f\"preds_{target_id}\"\n \n with open(fasta_path, \"w\") as f:\n f.write(f\">{target_id}\\n{sequence}\\n\")\n \n cmd = [\"python\", \"-u\", \"DRfold_infer.py\", str(fasta_path), str(output_dir), \"1\"]\n \n process = subprocess.Popen(\n cmd, \n cwd=str(DRFOLD_SCRIPT_DIR), \n stdout=subprocess.PIPE, \n stderr=subprocess.STDOUT, # Merges error messages into the same stream\n text=True,\n bufsize=1\n )\n\n for line in iter(process.stdout.readline, ''):\n # flush=True is critical for Jupyter Notebooks\n print(f\" [DRfold2 Log]: {line.strip()}\", flush=True)\n \n process.wait()\n\n if process.returncode == 0:\n relax_dir = output_dir / \"relax\"\n pdb_files = sorted(list(relax_dir.glob(\"model_*.pdb\")))\n if pdb_files:\n coords = extract_coords_from_pdb(pdb_files[0], len(sequence))\n shutil.rmtree(output_dir, ignore_errors=True)\n if fasta_path.exists(): os.remove(fasta_path)\n return coords\n return None\n\ndef predict_rna_structures_drfold2(sequence, target_id, original_target_id, msa_dir, train_seqs, train_coords_dict):\n \"\"\"Phase 3: Hybrid 3D Drafting with Logging, Grafting, and 2D Validation.\"\"\"\n \n # 1. 2D Consensus (From Phase 2 logic)\n consensus_ss = get_consensus_ss(sequence, target_id, msa_dir, original_target_id)\n hard_constraints = parse_pairs(consensus_ss)\n \n # 2. Generate Single Decoy\n coords = run_drfold2_instance(sequence, consensus_ss, target_id, original_target_id, msa_dir)\n \n if coords is None:\n return None # Trigger fallback\n\n # 3. Constraint Filtering\n\n results = check_physical_constraints(coords, hard_constraints)\n if not is_structure_acceptable(results):\n print(f\"โš ๏ธ {target_id} output violates 2D constraints. Falling back...\")\n return None\n\n # 4. Template Integration (Grafting)\n refined_coords = graft_functional_templates(coords, sequence, train_seqs, train_coords_dict)\n\n # 5. Competition Formatting (5 copies)\n return [refined_coords] * 5\n\ndef predict_rna_structures(sequence, target_id, train_seqs, train_coords_dict):\n \"\"\"Modernized template matching with connectivity gap-filling.\"\"\"\n predictions = []\n aligner = PairwiseAligner()\n aligner.match_score, aligner.mismatch_score = 2.0, -1.0\n aligner.open_gap_score, aligner.extend_gap_score = -0.5, -0.1\n aligner.mode = 'global'\n \n similarities = []\n for _, row in train_seqs.iterrows():\n t_id, t_seq = row['target_id'], row['sequence']\n # Filter by length to speed up alignment\n if abs(len(t_seq) - len(sequence)) > len(sequence) * 0.4: continue\n score = aligner.score(sequence, t_seq)\n similarities.append((t_id, t_seq, score))\n \n similarities.sort(key=lambda x: x[2], reverse=True)\n \n for t_id, t_seq, _ in similarities[:5]:\n t_coords = train_coords_dict.get(t_id)\n if t_coords is None: continue\n \n # Mapping coordinates query <-> template\n alignment = aligner.align(sequence, t_seq)[0]\n q_aln, t_aln = str(alignment[0]), str(alignment[1])\n \n new_coords, t_idx = [], 0\n for q_char, t_char in zip(q_aln, t_aln):\n if q_char != \"-\":\n if t_char != \"-\" and t_idx < len(t_coords):\n new_coords.append(t_coords[t_idx])\n else:\n # HEURISTIC: Maintain 3.8A backbone connectivity for gaps\n prev = new_coords[-1] if new_coords else np.array([0.0, 0.0, 0.0])\n new_coords.append(prev + np.array([3.8, 0, 0]))\n if t_char != \"-\": t_idx += 1\n \n predictions.append(new_coords[:len(sequence)])\n\n # Ensure we always return 5 models\n while len(predictions) < 5:\n predictions.append([(0.0, 0.0, 0.0)] * len(sequence))\n \n return predictions\n\n# --- ๐Ÿ› ๏ธ HELPER FUNCTIONS ---\n\ndef extract_coords_from_pdb(pdb_file, expected_len):\n parser = Bio.PDB.PDBParser(QUIET=True)\n structure = parser.get_structure(\"tmp\", pdb_file)\n coords = np.zeros((expected_len, 3))\n idx = 0\n for model in structure:\n for chain in model:\n for residue in chain:\n # Handle both prime and asterisk notations\n atom_key = \"C1'\" if \"C1'\" in residue else (\"C1*\" if \"C1*\" in residue else None)\n if atom_key:\n coords[idx] = residue[atom_key].get_coord()\n idx += 1\n if idx >= expected_len: break\n if idx >= expected_len: break\n return coords\n\ndef check_physical_constraints(coords, pairs, threshold=11.0):\n # Convert pairs to numpy array for indexing\n pairs = np.array(pairs)\n \n # Extract coordinates for all 'i' and 'j' in one go\n start_coords = coords[pairs[:, 0]]\n end_coords = coords[pairs[:, 1]]\n \n # Calculate all Euclidean distances simultaneously\n dists = np.linalg.norm(start_coords - end_coords, axis=1)\n \n # Calculate metrics\n num_violations = np.sum(dists > threshold)\n max_violation = np.max(dists) if len(dists) > 0 else 0\n accuracy_score = np.sum(dists <= threshold) / len(pairs) if len(pairs) > 0 else 1.0\n \n return {\n \"is_valid\": num_violations == 0,\n \"accuracy_score\": accuracy_score, # % of constraints met\n \"max_dist\": max_violation\n }\n\ndef is_structure_acceptable(results, min_accuracy=0.75, extreme_violation_limit=20.0):\n \"\"\"\n Determines if a DRFold2 structure is physically valid.\n \n Args:\n results (dict): Output from check_physical_constraints_pro\n min_accuracy (float): Minimum % of pairs that must be within threshold (default 85%)\n extreme_violation_limit (float): Absolute max distance allowed for ANY pair (default 20.0A)\n \n Returns:\n bool: True if structure is acceptable, False otherwise.\n \"\"\"\n # 1. Check if the global topology generally matches the 2D prediction\n has_high_accuracy = results['accuracy_score'] >= min_accuracy\n \n # 2. Check for \"Catastrophic Failures\" (extreme stretching)\n # Even if accuracy is high, one 50A jump means the chain is broken.\n is_physically_plausible = results['max_dist'] <= extreme_violation_limit\n \n # 3. Robust Decision\n return has_high_accuracy and is_physically_plausible\n\ndef parse_pairs(ss):\n stack, pairs = [], []\n for i, char in enumerate(ss):\n if char == '(': stack.append(i)\n elif char == ')':\n if stack: pairs.append((stack.pop(), i))\n return pairs","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:49.482017Z","iopub.status.busy":"2026-03-27T18:13:49.481609Z","iopub.status.idle":"2026-03-27T18:13:49.524425Z","shell.execute_reply":"2026-03-27T18:13:49.523744Z"},"papermill":{"duration":0.069165,"end_time":"2026-03-27T18:13:49.526424","exception":false,"start_time":"2026-03-27T18:13:49.457259","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import csv\nimport numpy as np\nimport pandas as pd\nfrom collections import defaultdict\n\n# Paths to your Kaggle datasets\nseq_path = '/kaggle/input/datasets/jaejohn/rna-all-data/merged_sequences_final.csv'\nlabels_path = '/kaggle/input/datasets/jaejohn/rna-all-data/merged_labels_final.csv'\n\n# ==========================================\n# 1. LOAD SEQUENCE METADATA (train_seqs)\n# ==========================================\nprint(\"โณ Loading sequence metadata...\")\ntrain_seqs_list = []\n\nwith open(seq_path, mode='r', encoding='utf-8') as f:\n reader = csv.DictReader(f)\n for row in reader:\n # We keep it as a list of dicts first to avoid the read_csv crash\n train_seqs_list.append(row)\n\n# Convert to DataFrame (Usually safer than read_csv in this environment)\ntrain_seqs = pd.DataFrame(train_seqs_list)\nprint(f\"โœ… Loaded {len(train_seqs)} sequences.\")\n\n# ==========================================\n# 2. LOAD COORDINATES (train_coords_dict)\n# ==========================================\ntemp_coords = defaultdict(list)\nprint(\"โณ Pre-indexing training coordinates (Handling missing values)...\")\n\nwith open(labels_path, mode='r', encoding='utf-8') as f:\n reader = csv.DictReader(f)\n for row in reader:\n # Validation: Skip rows where coordinates are empty strings\n if not row['x_1'] or not row['y_1'] or not row['z_1']:\n continue \n \n try:\n # Extract t_id from ID (e.g., 'target_1' -> 'target')\n # This matches your original .rsplit('_', 1)[0] logic\n t_id = row['ID'].rsplit('_', 1)[0]\n \n # Store [resid, x, y, z] for sorting\n temp_coords[t_id].append([\n int(row['resid']), \n float(row['x_1']), \n float(row['y_1']), \n float(row['z_1'])\n ])\n except (ValueError, KeyError):\n # Skip any malformed rows or headers\n continue\n\n# 3. Sort by resid and convert to final NumPy arrays\ntrain_coords_dict = {}\nfor t_id, data_list in temp_coords.items():\n # Sort by the residue number to ensure 3D order matches sequence order\n data_list.sort(key=lambda x: x[0])\n \n # Extract only the x, y, z values as a NumPy array\n train_coords_dict[t_id] = np.array([item[1:] for item in data_list])\n\nprint(f\"โœ… Successfully indexed {len(train_coords_dict)} structures into train_coords_dict.\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:13:49.574268Z","iopub.status.busy":"2026-03-27T18:13:49.573481Z","iopub.status.idle":"2026-03-27T18:14:46.832487Z","shell.execute_reply":"2026-03-27T18:14:46.831678Z"},"papermill":{"duration":57.308258,"end_time":"2026-03-27T18:14:46.857294","exception":false,"start_time":"2026-03-27T18:13:49.549036","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def save_coords_to_pdb(sequence, coords, target_id, output_path):\n \"\"\"Writes C1' coordinates to a PDB file format.\"\"\"\n with open(output_path, \"w\") as f:\n f.write(f\"HEADER RNA STRUCTURE PREDICTION - {target_id}\\n\")\n for i, (char, (x, y, z)) in enumerate(zip(sequence, coords)):\n # PDB ATOM format is very specific about column spacing\n # ATOM index name resname chain resID x y z\n f.write(f\"ATOM {i+1:>5} C1' {char:>3} A{i+1:>4} {x:>8.3f}{y:>8.3f}{z:>8.3f} 1.00 0.00 C\\n\")\n f.write(\"END\\n\")\n print(f\"๐Ÿ’พ Saved: {output_path}\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:46.902921Z","iopub.status.busy":"2026-03-27T18:14:46.90214Z","iopub.status.idle":"2026-03-27T18:14:46.908051Z","shell.execute_reply":"2026-03-27T18:14:46.907211Z"},"papermill":{"duration":0.030615,"end_time":"2026-03-27T18:14:46.90969","exception":false,"start_time":"2026-03-27T18:14:46.879075","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nfrom pathlib import Path\n\ndef process_single_target_to_pdb(target_id, original_target_id, sequence, msa_dir, train_seqs, train_coords_dict, output_folder):\n \"\"\"\n Predicts the structure for a single RNA target and saves the best model as a PDB.\n \n Args:\n target_id (str): The identifier for the sequence (e.g., 'R1107').\n sequence (str): The RNA sequence string.\n msa_dir (str/Path): Directory containing MSA files.\n train_seqs (pd.DataFrame): Reference sequences for template matching.\n train_coords_dict (dict): Indexed coordinates for reference sequences.\n output_folder (str/Path): Where to save the resulting .pdb file.\n \"\"\"\n output_folder = Path(output_folder) # โœ… FIX\n output_path = output_folder / f\"{target_id}_best.pdb\"\n preds = None\n\n # 1. Prediction Strategy Selection\n if len(sequence) < 400:\n print(f\"๐Ÿงฌ [{target_id}] Attempting DRFold2...\")\n try:\n preds = predict_rna_structures_drfold2(sequence, target_id, original_target_id, msa_dir, train_seqs, train_coords_dict)\n except Exception as e:\n print(f\"โš ๏ธ DRFold2 Error for {target_id}: {e}\")\n \n # If DRFold2 fails OR is rejected by constraint checks, trigger Fallback\n if preds is None:\n print(f\"๐Ÿ”„ [{target_id}] DRFold2 failed/rejected. Falling back to Template Matching...\")\n preds = predict_rna_structures(sequence, target_id, train_seqs, train_coords_dict)\n else:\n print(f\"๐Ÿ“ [{target_id}] Length {len(sequence)} > 400. Using Template Matching fallback.\")\n preds = predict_rna_structures(sequence, target_id, train_seqs, train_coords_dict)\n\n # 2. Save Resulting Coordinates\n if preds and len(preds) > 0:\n # We take the first model (index 0) as the 'best' candidate\n coords = preds[0]\n consensus_ss = get_consensus_ss(sequence, target_id, msa_dir, original_target_id)\n \n # --- PHASE 4.1: gCoSRNA Junction Alignment (in-process; moves atoms) ---\n try:\n model_pkl = \"/kaggle/input/datasets/khushvaghela/gcosrna/gCoSRNA_model.pkl\"\n\n # Predicted coaxial stacks WITH residue indices + pivots\n coaxial_pairs = get_coaxial_stem_pairs(sequence, consensus_ss, model_pkl)\n\n if coaxial_pairs:\n print(f\"๐Ÿ“ Adjusting 3D geometry based on {len(coaxial_pairs)} coaxial stack(s)...\")\n coords = apply_coaxial_alignment(coords, sequence, coaxial_pairs)\n\n except Exception as e:\n print(f\"โš ๏ธ Junction Alignment skipped due to error: {e}\")\n \n # --- PHASE 4.2: Monomer Settling (SimRNA) ---\n # print(f\"๐Ÿง˜ Phase 4.2: SimRNA Loop Relaxation for {target_id}\")\n # try:\n # # 1000 steps is a good \"settling\" start; \n # # increase for a deeper 15-minute soak if needed.\n # coords = run_simrna_settling(target_id, sequence, coords, steps=1000)\n # except Exception as e:\n # print(f\"โš ๏ธ SimRNA failed: {e}. Keeping aligned draft.\")\n \n # --- FINAL SAVE ---\n output_path = os.path.join(output_folder, f\"{target_id}_final_monomer.pdb\")\n save_coords_to_pdb(sequence, coords, target_id, output_path)\n else:\n print(f\"โŒ [{target_id}] Failed to generate any valid coordinates.\")\n\n # 3. Memory Cleanup (Critical for consecutive GPU calls)\n if torch.cuda.is_available():\n torch.cuda.empty_cache()","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:46.956767Z","iopub.status.busy":"2026-03-27T18:14:46.956411Z","iopub.status.idle":"2026-03-27T18:14:46.966882Z","shell.execute_reply":"2026-03-27T18:14:46.965943Z"},"papermill":{"duration":0.036628,"end_time":"2026-03-27T18:14:46.968682","exception":false,"start_time":"2026-03-27T18:14:46.932054","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import time\nimport pandas as pd\nimport torch\nimport numpy as np\n\ndef run_drfold2_inference(test_csv_path, msa_dir, train_seqs, train_coords_dict, output_csv=\"submission.csv\"):\n \"\"\"\n Orchestrates the full inference pipeline:\n 1. Loads test sequences.\n 2. Decides between DRFold2 and Template Matching based on length.\n 3. Handles time-outs and GPU memory.\n 4. Formats and saves the final submission.\n \"\"\"\n # Load the target sequences\n test_sequences = pd.read_csv(test_csv_path)\n all_results = []\n start_time = time.time()\n \n print(f\"๐Ÿ Starting Inference on {len(test_sequences)} targets...\")\n\n for idx, row in test_sequences.iterrows():\n # 1. 8-hour safety limit (Kaggle limit is 9h, leaving 1h for saving/cleanup)\n elapsed = time.time() - start_time\n if elapsed > (7 * 3600): \n print(\"๐Ÿ•’ Time limit approaching (8 hours reached). Wrapping up current results...\")\n break \n \n target_id = str(row['target_id'])\n sequence = str(row['sequence'])\n seq_len = len(sequence)\n \n print(f\"\\n[{idx+1}/{len(test_sequences)}] ๐Ÿš€ Processing {target_id} (Length: {seq_len})\")\n\n # 2. Hybrid Decision Logic\n preds = None\n \n if seq_len < 400:\n # Try DRFold2 first for smaller/medium RNAs\n try:\n # FIX: original call omitted `original_target_id`, silently shifting every\n # positional arg. (This run_drfold2_inference is legacy/unused โ€” run_main_loop\n # is the real entry point โ€” but corrected for safety.)\n preds = predict_rna_structures_drfold2(sequence, target_id, target_id, msa_dir, train_seqs, train_coords_dict)\n except Exception as e:\n print(f\"โš ๏ธ DRFold2 encountered an error for {target_id}: {e}\")\n \n # Fallback to Template Matching if DRFold2 fails or returns None\n if preds is None:\n print(f\"๐Ÿ”„ DRFold2 failed or was rejected. Falling back to Template Matching for {target_id}...\")\n preds = predict_rna_structures(sequence, target_id, train_seqs, train_coords_dict)\n else:\n # Direct Template Matching for long sequences to avoid OOM\n print(f\"๐Ÿ“ Length {seq_len} exceeds DRFold2 threshold. Using Template Matching.\")\n preds = predict_rna_structures(sequence, target_id, train_seqs, train_coords_dict)\n\n # 3. Final Result Verification & Formatting\n # If everything failed, provide zero-coordinates to ensure a valid submission row\n if preds is None or len(preds) < 5:\n print(f\"โŒ Critical Failure for {target_id}. Filling with zero-coordinates.\")\n preds = [[(0.0, 0.0, 0.0)] * seq_len for _ in range(5)]\n\n # 4. Row-wise Formatting (Competition Schema)\n for j in range(seq_len):\n # The 'ID' is usually target_id_residueIndex\n res_entry = {\n 'ID': f\"{target_id}_{j+1}\", \n 'resname': sequence[j], \n 'resid': j+1\n }\n # Unpack the 5 sets of x, y, z coordinates\n for i in range(5):\n x, y, z = preds[i][j]\n res_entry[f'x_{i+1}'] = x\n res_entry[f'y_{i+1}'] = y\n res_entry[f'z_{i+1}'] = z\n \n all_results.append(res_entry)\n \n # 5. Maintenance\n if torch.cuda.is_available():\n torch.cuda.empty_cache()\n\n # 6. Final Save\n if not all_results:\n print(\"Empty results. No submission file generated.\")\n return None\n\n submission_df = pd.DataFrame(all_results)\n submission_df.to_csv(output_csv, index=False)\n \n total_time = (time.time() - start_time) / 3600\n print(f\"\\nโœ… Final Submission Generated: {output_csv}\")\n print(f\"โฑ๏ธ Total runtime: {total_time:.2f} hours\")\n \n return submission_df\n\n# --- EXAMPLE CALL ---\n# submission = run_drfold2_inference(\n# test_csv_path=\"/kaggle/input/datasets/khushvaghela/drfold2/test_sequences.csv\",\n# msa_dir=\"/kaggle/input/competition-msas/MSA\",\n# train_seqs=train_seqs,\n# train_coords_dict=train_coords_dict\n# )","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.015248Z","iopub.status.busy":"2026-03-27T18:14:47.014792Z","iopub.status.idle":"2026-03-27T18:14:47.027567Z","shell.execute_reply":"2026-03-27T18:14:47.026448Z"},"papermill":{"duration":0.039297,"end_time":"2026-03-27T18:14:47.030233","exception":false,"start_time":"2026-03-27T18:14:46.990936","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# output_folder = \"/kaggle/working/output_pdbs\"\n# msa_dir = \"/kaggle/input/stanford-rna-3d-folding-2/MSA\"\n\n# process_single_target_to_pdb(\n# \"8ZNQ\",\n# \"8ZNQ1\",\n# \"ACCGUGACGGGCCUUUUGGCUAUACGCGGU\",\n# msa_dir=msa_dir,\n# train_seqs=train_seqs,\n# train_coords_dict=train_coords_dict,\n# output_folder=output_folder\n# )","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.079164Z","iopub.status.busy":"2026-03-27T18:14:47.078816Z","iopub.status.idle":"2026-03-27T18:14:47.082958Z","shell.execute_reply":"2026-03-27T18:14:47.082124Z"},"papermill":{"duration":0.030701,"end_time":"2026-03-27T18:14:47.084694","exception":false,"start_time":"2026-03-27T18:14:47.053993","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## SimRNA","metadata":{"papermill":{"duration":0.022219,"end_time":"2026-03-27T18:14:47.129389","exception":false,"start_time":"2026-03-27T18:14:47.10717","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# import os\n# import shutil\n# from pathlib import Path\n\n# # Paths - Adjust 'simrna-standalone' to your actual dataset name\n# SIMRNA_SRC = Path(\"/kaggle/input/datasets/khushvaghela/simrna/SimRNA_64bitIntel_Linux\") \n# SIMRNA_PKG = Path(\"/kaggle/working/SimRNA_pkg\")\n\n# def setup_simrna():\n# if SIMRNA_PKG.exists():\n# shutil.rmtree(SIMRNA_PKG)\n \n# # Copy the package to a writable directory\n# shutil.copytree(SIMRNA_SRC, SIMRNA_PKG)\n \n# # Ensure the 'data' folder is directly inside the working directory [cite: 13]\n# # The manual states SimRNA will crash without this [cite: 18, 19]\n# data_dir = SIMRNA_PKG / \"data\"\n# if not data_dir.exists():\n# print(\"โŒ Error: 'data' folder not found in the package!\")\n# else:\n# print(\"โœ… SimRNA data folder (with .hist files) is in place[cite: 12].\")\n\n# # Set execution permissions for the binary\n# simrna_bin = SIMRNA_PKG / \"SimRNA\"\n# if simrna_bin.exists():\n# os.chmod(simrna_bin, 0o755)\n# print(\"โœ… SimRNA binary is executable.\")\n\n# setup_simrna()\n\n# import os\n# import shutil\n# from pathlib import Path\n\n# def batch_repair_simrna_data(data_path):\n# data_dir = Path(data_path)\n# if not data_dir.exists():\n# print(f\"โŒ Error: Directory {data_dir} not found.\")\n# return\n\n# hist_files = list(data_dir.glob(\"*.hist\"))\n# repaired_count = 0\n \n# print(f\"๐Ÿ” Scanning {len(hist_files)} .hist files for corruption...\")\n\n# for target_file in hist_files:\n# try:\n# with open(target_file, 'r') as f:\n# first_line = f.readline().strip()\n# # Check if first line looks like a filename and line 2 is empty\n# second_line = f.readline().strip()\n \n# # If line 1 is another file in the same folder and line 2 is empty...\n# source_file = data_dir / first_line\n# if not second_line and source_file.exists() and source_file != target_file:\n# print(f\"๐Ÿ›  Repairing: {target_file.name} (Placeholder for {first_line})\")\n \n# # Backup and replace\n# shutil.copy(source_file, f\"{target_file}.bak\")\n# shutil.copy(source_file, target_file)\n# repaired_count += 1\n \n# except Exception as e:\n# print(f\"โš ๏ธ Could not check {target_file.name}: {e}\")\n\n# print(f\"\\nโœ… Finished. Repaired {repaired_count} files.\")\n\n# # Execute repair on your specific folder\n# batch_repair_simrna_data(SIMRNA_PKG / \"data\")\n\n","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.176158Z","iopub.status.busy":"2026-03-27T18:14:47.175642Z","iopub.status.idle":"2026-03-27T18:14:47.180951Z","shell.execute_reply":"2026-03-27T18:14:47.180212Z"},"papermill":{"duration":0.030932,"end_time":"2026-03-27T18:14:47.182645","exception":false,"start_time":"2026-03-27T18:14:47.151713","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import os\n# import shutil\n# from pathlib import Path\n\n# def batch_repair_simrna_data(data_path):\n# data_dir = Path(data_path)\n# if not data_dir.exists():\n# print(f\"โŒ Error: Directory {data_dir} not found.\")\n# return\n\n# hist_files = list(data_dir.glob(\"*.hist\"))\n# repaired_count = 0\n \n# print(f\"๐Ÿ” Scanning {len(hist_files)} .hist files for corruption...\")\n\n# for target_file in hist_files:\n# try:\n# with open(target_file, 'r') as f:\n# first_line = f.readline().strip()\n# # Check if first line looks like a filename and line 2 is empty\n# second_line = f.readline().strip()\n \n# # If line 1 is another file in the same folder and line 2 is empty...\n# source_file = data_dir / first_line\n# if not second_line and source_file.exists() and source_file != target_file:\n# print(f\"๐Ÿ›  Repairing: {target_file.name} (Placeholder for {first_line})\")\n \n# # Backup and replace\n# shutil.copy(source_file, f\"{target_file}.bak\")\n# shutil.copy(source_file, target_file)\n# repaired_count += 1\n \n# except Exception as e:\n# print(f\"โš ๏ธ Could not check {target_file.name}: {e}\")\n\n# print(f\"\\nโœ… Finished. Repaired {repaired_count} files.\")\n\n# # Execute repair on your specific folder\n# batch_repair_simrna_data(SIMRNA_PKG / \"data\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.23037Z","iopub.status.busy":"2026-03-27T18:14:47.229291Z","iopub.status.idle":"2026-03-27T18:14:47.234635Z","shell.execute_reply":"2026-03-27T18:14:47.233911Z"},"papermill":{"duration":0.031053,"end_time":"2026-03-27T18:14:47.236365","exception":false,"start_time":"2026-03-27T18:14:47.205312","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# # Define your path (assuming SimRNA_pkg is your working directory)\n# seq_file_path = \"/kaggle/working/SimRNA_pkg/test.seq\"\n# multi_seq_file_path = \"/kaggle/working/SimRNA_pkg/multi_test.seq\"\n\n# # 1. Generate a Single-Strand File\n# # Standard unbroken sequence\n# single_strand = \"GGAGGUCUUC\" \n\n# with open(seq_file_path, \"w\") as f:\n# f.write(single_strand)\n# print(f\"โœ… Created single-strand file: {seq_file_path}\")\n\n# # 2. Generate a Multi-Strand File\n# # Independent strands separated by a space\n# multi_strand = \"GGAGGUCUUC AACCGGUU\" \n\n# with open(multi_seq_file_path, \"w\") as f:\n# f.write(multi_strand)\n# print(f\"โœ… Created multi-strand file: {multi_seq_file_path}\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.282934Z","iopub.status.busy":"2026-03-27T18:14:47.282444Z","iopub.status.idle":"2026-03-27T18:14:47.286768Z","shell.execute_reply":"2026-03-27T18:14:47.285985Z"},"papermill":{"duration":0.02948,"end_time":"2026-03-27T18:14:47.288461","exception":false,"start_time":"2026-03-27T18:14:47.258981","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import subprocess\n# import os\n# from pathlib import Path\n\n# # --- UPDATE THIS LINE based on the 'find' command output ---\n# # Example: if it found /kaggle/working/SimRNA_pkg/bin/SimRNA\n# simrna_bin_abs = Path(\"/kaggle/working/SimRNA_pkg/SimRNA\").resolve() \n\n# # 1. Verification and Permission Fix\n# if not simrna_bin_abs.exists():\n# # If the direct path failed, try to auto-detect\n# potential_bins = list(Path(\"/kaggle/working/SimRNA_pkg\").glob(\"**/SimRNA*\"))\n# if potential_bins:\n# simrna_bin_abs = potential_bins[0].resolve()\n# print(f\"โœ… Auto-detected binary at: {simrna_bin_abs}\")\n# else:\n# raise FileNotFoundError(f\"โŒ Binary NOT found in SimRNA_pkg. Check your dataset.\")\n\n# os.chmod(simrna_bin_abs, 0o755)\n\n# # 2. Hardened Command Launch\n# # Ensure all paths in 'cmd' are strings and absolute\n# test_seq_abs = Path(\"/kaggle/working/SimRNA_pkg/test.seq\").resolve()\n# cmd = [str(simrna_bin_abs), \"-s\", str(test_seq_abs), \"-n\", \"100\", \"-o\", \"test_output\"]\n\n# print(f\"๐Ÿš€ Launching SimRNA (Raw Output Mode)...\")\n# process = subprocess.Popen(\n# cmd,\n# cwd=str(SIMRNA_PKG), # Crucial for finding the 'data' folder\n# stdout=subprocess.PIPE,\n# stderr=subprocess.STDOUT,\n# text=True\n# )\n\n# for line in process.stdout:\n# print(f\" [SimRNA Raw]: {line.strip()}\")\n\n# process.wait()\n# print(f\"\\nTerminated with Exit Code: {process.returncode}\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.33552Z","iopub.status.busy":"2026-03-27T18:14:47.335049Z","iopub.status.idle":"2026-03-27T18:14:47.339723Z","shell.execute_reply":"2026-03-27T18:14:47.338923Z"},"papermill":{"duration":0.030388,"end_time":"2026-03-27T18:14:47.341532","exception":false,"start_time":"2026-03-27T18:14:47.311144","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## HDock","metadata":{"papermill":{"duration":0.022326,"end_time":"2026-03-27T18:14:47.386771","exception":false,"start_time":"2026-03-27T18:14:47.364445","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import os\nimport shutil\nimport subprocess\nfrom pathlib import Path\n\n# Paths - Adjust 'hdock-lite' to your actual dataset name\nHDOCK_SRC = Path(\"/kaggle/input/datasets/khushvaghela/h-dock/HDOCKlite-v1.1\") \nHDOCK_PKG = Path(\"/kaggle/working/HDOCK_pkg\")\n\ndef setup_hdock():\n if HDOCK_PKG.exists():\n shutil.rmtree(HDOCK_PKG)\n \n # 1. Copy the package\n shutil.copytree(HDOCK_SRC, HDOCK_PKG)\n \n # 2. Grant permissions to the core binaries\n # Usually includes 'hdock' and 'createpl'\n binaries = [\"hdock\", \"createpl\"]\n for bin_name in binaries:\n bin_path = HDOCK_PKG / bin_name\n if bin_path.exists():\n os.chmod(bin_path, 0o755)\n print(f\"โœ… {bin_name} is ready.\")\n else:\n # Check for versions with extensions or in subfolders\n found = list(HDOCK_PKG.glob(f\"**/{bin_name}*\"))\n if found:\n os.chmod(found[0], 0o755)\n print(f\"โœ… Found and readied {bin_name} at: {found[0].name}\")\n\nsetup_hdock()","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.433901Z","iopub.status.busy":"2026-03-27T18:14:47.433385Z","iopub.status.idle":"2026-03-27T18:14:47.583126Z","shell.execute_reply":"2026-03-27T18:14:47.582105Z"},"papermill":{"duration":0.17567,"end_time":"2026-03-27T18:14:47.584979","exception":false,"start_time":"2026-03-27T18:14:47.409309","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# def run_hdock_test():\n# # Find the actual binaries\n# hdock_bin = next(HDOCK_PKG.glob(\"**/hdock*\"), None)\n# createpl_bin = next(HDOCK_PKG.glob(\"**/createpl*\"), None)\n \n# if not hdock_bin or not createpl_bin:\n# print(\"โŒ Missing binaries. Check setup.\")\n# return\n\n# # Use a sample PDB from your previous SimRNA test if available\n# # Or point to a DRfold2 output\n# sample_pdb = \"/kaggle/working/SimRNA_pkg/test_output-000001.pdb\" \n \n# if not Path(sample_pdb).exists():\n# print(\"โš ๏ธ No sample PDB found. Please run your SimRNA test first.\")\n# return\n\n# print(\"๐Ÿš€ Launching HDOCK ab initio docking...\")\n \n# # Step A: Perform Docking\n# # Command: hdock receptor.pdb ligand.pdb -out Hdock.out\n# cmd1 = [str(hdock_bin), sample_pdb, sample_pdb, \"-out\", \"Hdock.out\"]\n \n# try:\n# subprocess.run(cmd1, cwd=str(HDOCK_PKG), check=True, capture_output=True, text=True)\n# print(\"โœ… Docking calculation complete (Hdock.out generated).\")\n \n# # Step B: Create Models\n# # Command: createpl Hdock.out top10.pdb -nmax 10 -complex -models\n# cmd2 = [str(createpl_bin), \"Hdock.out\", \"top10.pdb\", \"-nmax\", \"10\", \"-complex\", \"-models\"]\n# subprocess.run(cmd2, cwd=str(HDOCK_PKG), check=True, capture_output=True, text=True)\n# print(\"โœ… Top 10 models generated.\")\n \n# except subprocess.CalledProcessError as e:\n# print(f\"โŒ HDOCK failed: {e.stderr}\")\n\n# run_hdock_test()","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.632843Z","iopub.status.busy":"2026-03-27T18:14:47.632054Z","iopub.status.idle":"2026-03-27T18:14:47.63729Z","shell.execute_reply":"2026-03-27T18:14:47.636462Z"},"papermill":{"duration":0.031793,"end_time":"2026-03-27T18:14:47.6393","exception":false,"start_time":"2026-03-27T18:14:47.607507","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import subprocess\nimport shutil\nimport os\nfrom pathlib import Path\n\n# --- CONFIGURATION ---\nHDOCK_PKG = Path(\"/kaggle/working/HDOCK_pkg\").resolve()\n\ndef run_hdock_assembly(receptor_pdb, ligand_pdb, output_folder, target_id):\n \"\"\"\n Phase 5: Docking with Absolute Path Enforcement.\n \"\"\"\n # Force everything to absolute paths immediately\n receptor_pdb = Path(receptor_pdb).resolve()\n ligand_pdb = Path(ligand_pdb).resolve()\n output_folder = Path(output_folder).resolve()\n \n # 1. Locate binaries\n hdock_bin = next(HDOCK_PKG.glob(\"**/hdock*\"), None)\n createpl_bin = next(HDOCK_PKG.glob(\"**/createpl*\"), None)\n \n if not hdock_bin or not createpl_bin:\n print(\"โŒ HDOCK binaries not found.\")\n return None\n\n # 2. Setup Intermediate Paths\n work_dir = Path(f\"/kaggle/working/dock_{target_id}\").resolve()\n work_dir.mkdir(parents=True, exist_ok=True)\n \n dock_out = work_dir / \"Hdock.out\"\n final_complex_path = output_folder / f\"{target_id}_complex_best.pdb\"\n\n print(f\"๐Ÿšข HDOCK: Docking {target_id}...\")\n\n try:\n # Step A: Run docking using absolute strings\n subprocess.run(\n [str(hdock_bin), str(receptor_pdb), str(ligand_pdb), \"-out\", str(dock_out)],\n cwd=str(work_dir),\n check=True,\n capture_output=True,\n text=True\n )\n\n # Step B: Generate top model\n subprocess.run(\n [str(createpl_bin), \"Hdock.out\", \"top1.pdb\", \"-nmax\", \"1\", \"-complex\", \"-models\"],\n cwd=str(work_dir),\n check=True,\n capture_output=True,\n text=True\n )\n\n model_1 = work_dir / \"model_1.pdb\"\n if model_1.exists():\n shutil.move(str(model_1), str(final_complex_path))\n shutil.rmtree(work_dir)\n return str(final_complex_path)\n \n except subprocess.CalledProcessError as e:\n print(f\"โŒ HDOCK failed: {e.stderr}\")\n return None","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.686417Z","iopub.status.busy":"2026-03-27T18:14:47.685868Z","iopub.status.idle":"2026-03-27T18:14:47.695385Z","shell.execute_reply":"2026-03-27T18:14:47.694494Z"},"papermill":{"duration":0.035508,"end_time":"2026-03-27T18:14:47.697289","exception":false,"start_time":"2026-03-27T18:14:47.661781","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# # Assuming you already have your relaxed RNA monomer from SimRNA\n# rna_monomer = \"/kaggle/working/SimRNA_pkg/test_output-000001.pdb\"\n# protein_receptor = \"/kaggle/working/SimRNA_pkg/test_output-000001.pdb\"\n\n# complex_pdb = run_hdock_assembly(\n# receptor_pdb=protein_receptor,\n# ligand_pdb=rna_monomer,\n# output_folder=\"/kaggle/working/final_results\",\n# target_id=\"R1107\"\n# )","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.743189Z","iopub.status.busy":"2026-03-27T18:14:47.742779Z","iopub.status.idle":"2026-03-27T18:14:47.747076Z","shell.execute_reply":"2026-03-27T18:14:47.746146Z"},"papermill":{"duration":0.029714,"end_time":"2026-03-27T18:14:47.748896","exception":false,"start_time":"2026-03-27T18:14:47.719182","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Tinker","metadata":{"papermill":{"duration":0.022548,"end_time":"2026-03-27T18:14:47.79387","exception":false,"start_time":"2026-03-27T18:14:47.771322","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import os\nimport shutil\nimport subprocess\nfrom pathlib import Path\n\n# --- 1. CONFIGURATION ---\n# Base path where you uploaded the Tinker source\nTINKER_SRC = Path(\"/kaggle/input/datasets/khushvaghela/tinker-src/tinker-main/tinker-main\")\nTINKER_CMAKE_DIR = TINKER_SRC / \"cmake\"\n\n# Path to your offline FFTW dataset\n# โš ๏ธ Double check this path in your Kaggle sidebar! \n# It might be /kaggle/input/fftw-offline-kit/fftw_offline_kit\nFFTW_KIT_ROOT = Path(\"/kaggle/input/datasets/khushvaghela/tinker-src/fftw_offline_kit/fftw_offline_kit\")\n\nBUILD_DIR = Path(\"/kaggle/working/tinker_build\")\nBIN_DEST = Path(\"/kaggle/working/tinker_bin\")\n\ndef build_tinker_offline_v2():\n # Path Validation\n if not FFTW_KIT_ROOT.exists():\n print(f\"โŒ ERROR: FFTW Kit not found at {FFTW_KIT_ROOT}\")\n print(\"Please verify the path in the Kaggle 'Data' sidebar.\")\n return False\n \n if BUILD_DIR.exists(): shutil.rmtree(BUILD_DIR)\n BUILD_DIR.mkdir(parents=True, exist_ok=True)\n BIN_DEST.mkdir(parents=True, exist_ok=True)\n\n # 2. Locate the specific .so files\n # We use the exact names the CMakeLists.txt logic expects\n fftw_lib = FFTW_KIT_ROOT / \"lib/libfftw3.so\"\n fftw_threaded = FFTW_KIT_ROOT / \"lib/libfftw3_omp.so\" # Or libfftw3_threads.so\n\n # Fallback: if omp version isn't there, use the main lib for both\n if not fftw_threaded.exists():\n fftw_threaded = fftw_lib\n\n print(\"๐Ÿ”จ Step 1: Configuring with CMake (Force-feeding FFTW paths)...\")\n \n # We use -D to set the exact variables used in your CMakeLists.txt lines 44-60\n cmake_cmd = [\n \"cmake\",\n str(TINKER_CMAKE_DIR),\n f\"-DFFTW_LIB={fftw_lib}\",\n f\"-DFFTW_THREADED_LIB={fftw_threaded}\",\n f\"-DCMAKE_INCLUDE_PATH={FFTW_KIT_ROOT}/include\",\n \"-DCMAKE_Fortran_FLAGS='-O3 -fno-range-check -fallow-argument-mismatch'\",\n \"-DCMAKE_INSTALL_RPATH_USE_LINK_PATH=TRUE\"\n ]\n \n # Execute configuration\n result_config = subprocess.run(\" \".join(cmake_cmd), cwd=BUILD_DIR, shell=True, capture_output=True, text=True)\n \n if result_config.returncode != 0:\n print(f\"โŒ Configuration Failed:\\n{result_config.stderr}\")\n return False\n\n print(\"๐Ÿ”จ Step 2: Compiling specific tools...\")\n # Building only what we need for the competition\n make_cmd = \"make pdbxyz minimize xyzpdb -j4\"\n result_make = subprocess.run(make_cmd, cwd=BUILD_DIR, shell=True, capture_output=True, text=True)\n\n if result_make.returncode == 0:\n for t in [\"pdbxyz\", \"minimize\", \"xyzpdb\"]:\n shutil.move(BUILD_DIR / t, BIN_DEST / t)\n os.chmod(BIN_DEST / t, 0o755)\n print(f\"โœ… Success! Binaries ready in {BIN_DEST}\")\n return True\n else:\n print(f\"โŒ Compilation Failed. Log:\\n{result_make.stderr}\")\n return False\n\n# Run it\nbuild_tinker_offline_v2()","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:14:47.842209Z","iopub.status.busy":"2026-03-27T18:14:47.841406Z","iopub.status.idle":"2026-03-27T18:16:35.524293Z","shell.execute_reply":"2026-03-27T18:16:35.523278Z"},"papermill":{"duration":107.732137,"end_time":"2026-03-27T18:16:35.548635","exception":false,"start_time":"2026-03-27T18:14:47.816498","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport shutil\nimport subprocess\nfrom pathlib import Path\n\n# --- Global Config ---\nTINKER_BIN_DIR = Path(\"/kaggle/working/tinker_bin\")\nTINKER_PARAMS = \"/kaggle/input/datasets/khushvaghela/tinker-src/tinker-main/tinker-main/params/amber99.prm\"\n\ndef tinker_refine_pdb(input_pdb, output_dir, gradient=1.0):\n \"\"\"\n Performs full Tinker refinement using target-specific temporary space.\n \"\"\"\n # 1. PATH CHANGE: Make work_dir a SUBFOLDER of the target folder\n # This prevents different targets from overwriting each other's temp files\n output_dir = Path(output_dir)\n work_dir = output_dir / \"tinker_work\" \n \n if work_dir.exists():\n shutil.rmtree(work_dir)\n work_dir.mkdir(parents=True, exist_ok=True)\n\n base_name = \"refine_task\"\n temp_pdb = work_dir / f\"{base_name}.pdb\"\n final_output_path = output_dir / \"final_refined.pdb\"\n\n # 2. PATH CHANGE: Ensure these are defined or passed in!\n # If these are globally defined at the top of your notebook, you are fine.\n # Otherwise, you need to define them here.\n # TINKER_BIN_DIR = Path(\"/opt/tinker/bin\") \n # TINKER_PARAMS = Path(\"/opt/tinker/params/amber99.prm\")\n\n # [Sanitization Logic stays the same...]\n rna_map = {'A': 'RA', 'C': 'RC', 'G': 'RG', 'U': 'RU'}\n print(f\"๐Ÿงน Sanitizing {Path(input_pdb).name} for Tinker...\")\n\n with open(input_pdb, 'r') as f_in, open(temp_pdb, 'w') as f_out:\n for line in f_in:\n if line.startswith(\"ATOM\"):\n res_name = line[17:20].strip()\n if res_name in rna_map:\n line = line[:17] + rna_map[res_name].ljust(3) + line[20:]\n \n atom_name = line[12:16]\n if \"'\" in atom_name:\n line = line[:12] + atom_name.replace(\"'\", \"*\") + line[16:]\n f_out.write(line)\n\n # 3. Create Key file inside work_dir\n with open(work_dir / f\"{base_name}.key\", \"w\") as k:\n k.write(f\"parameters {TINKER_PARAMS}\\n\")\n k.write(\"verbose\\n\")\n\n # 4. Phase A: pdbxyz (The Mapping)\n # LOGIC TWEAK: pdbxyz sometimes needs more/fewer \"1\" inputs \n # depending on the number of chains. \"1\\n\" is usually enough for RNA.\n print(\"โณ Phase A: Mapping atoms (pdbxyz)...\")\n subprocess.run(\n [str(TINKER_BIN_DIR / \"pdbxyz\"), base_name],\n cwd=work_dir,\n input=\"1\\n\", \n text=True,\n capture_output=True\n )\n\n if not (work_dir / f\"{base_name}.xyz\").exists():\n print(f\"โŒ Mapping failed for {output_dir.name}. Tinker couldn't recognize residues.\")\n return False\n\n # 5. Phase B: minimize\n print(f\"๐Ÿš€ Phase B: Minimizing (gradient={gradient})...\")\n os.environ[\"OMP_NUM_THREADS\"] = \"4\"\n subprocess.run(\n [str(TINKER_BIN_DIR / \"minimize\"), base_name, str(gradient)],\n cwd=work_dir,\n input=\"\\n\",\n text=True,\n capture_output=True\n )\n\n # 6. Phase C: xyzpdb (The Export)\n # Tinker versioning: The minimized file is usually base_name.xyz_2\n xyz_out = work_dir / f\"{base_name}.xyz_2\"\n if not xyz_out.exists():\n xyz_out = work_dir / f\"{base_name}.xyz\" # Fallback if no movement occurred\n\n if xyz_out.exists():\n print(\"โณ Phase C: Exporting PDB (xyzpdb)...\")\n subprocess.run(\n [str(TINKER_BIN_DIR / \"xyzpdb\"), xyz_out.name],\n cwd=work_dir,\n input=\"1\\n\", # Selection for the parameter set\n text=True,\n capture_output=True\n )\n\n # Tinker output can be .pdb_2 or just .pdb\n pdb_out = work_dir / f\"{base_name}.pdb_2\"\n if not pdb_out.exists():\n pdb_out = work_dir / f\"{base_name}.pdb\"\n\n if pdb_out.exists():\n shutil.move(str(pdb_out), str(final_output_path))\n print(f\"๐Ÿ† SUCCESS! Saved as: {final_output_path}\")\n return True\n\n print(f\"โŒ Refinement failed for {output_dir.name}.\")\n return False","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:16:35.595994Z","iopub.status.busy":"2026-03-27T18:16:35.595279Z","iopub.status.idle":"2026-03-27T18:16:35.608303Z","shell.execute_reply":"2026-03-27T18:16:35.607501Z"},"papermill":{"duration":0.038907,"end_time":"2026-03-27T18:16:35.610063","exception":false,"start_time":"2026-03-27T18:16:35.571156","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Protenix","metadata":{"papermill":{"duration":0.022289,"end_time":"2026-03-27T18:16:35.655016","exception":false,"start_time":"2026-03-27T18:16:35.632727","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import subprocess, glob, os\n\n# Install core bioinformatics libraries from local wheels\npackages = {\n 'biopython': '/kaggle/input/datasets/kami1976/biopython-cp312/*.whl',\n 'biotite': '/kaggle/input/datasets/amirrezaaleyasin/biotite/*.whl',\n 'rdkit': '/kaggle/input/datasets/amirrezaaleyasin/rdkit-2025-9-5/*.whl'\n}\n\nfor name, pattern in packages.items():\n whls = glob.glob(pattern)\n if whls:\n subprocess.run(['pip', 'install', '--no-index', '--no-deps', whls[0]], check=True)\n print(f'Successfully installed {name} from {whls[0]}')\n else:\n print(f'WARNING: Could not find wheels for {name}. Ensure the dataset is added to the notebook.')","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:16:35.701831Z","iopub.status.busy":"2026-03-27T18:16:35.701339Z","iopub.status.idle":"2026-03-27T18:16:44.669976Z","shell.execute_reply":"2026-03-27T18:16:44.669105Z"},"papermill":{"duration":8.994448,"end_time":"2026-03-27T18:16:44.672233","exception":false,"start_time":"2026-03-27T18:16:35.677785","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import sys\nfrom pathlib import Path\n\n# Paths to Protenix code and competition data\nPROTENIX_CODE_DIR = \"/kaggle/input/datasets/qiweiyin/protenix-v1-adjusted/Protenix-v1-adjust-v2/Protenix-v1-adjust-v2/Protenix-v1\"\nDATA_BASE = \"/kaggle/input/competitions/stanford-rna-3d-folding-2\"\n\n# Required for importing from the Protenix package\nif PROTENIX_CODE_DIR not in sys.path:\n sys.path.append(PROTENIX_CODE_DIR)\n\n# Set model-specific environment variables\nos.environ[\"LAYERNORM_TYPE\"] = \"torch\"\nos.environ[\"RNA_MSA_DEPTH_LIMIT\"] = \"512\"\nos.environ[\"PROTENIX_ROOT_DIR\"] = PROTENIX_CODE_DIR\n\nprint(f\"Environment set. Protenix directory: {PROTENIX_CODE_DIR}\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:16:44.721257Z","iopub.status.busy":"2026-03-27T18:16:44.720269Z","iopub.status.idle":"2026-03-27T18:16:44.726477Z","shell.execute_reply":"2026-03-27T18:16:44.72558Z"},"papermill":{"duration":0.03226,"end_time":"2026-03-27T18:16:44.728196","exception":false,"start_time":"2026-03-27T18:16:44.695936","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import json\nimport torch\nfrom protenix.config.config import parse_configs\n\ndef build_protenix_configs(input_json_path, dump_dir, model_name=\"protenix_base_20250630_v1.0.0\"):\n # Import config modules from the Protenix code\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\n # Merge base configurations\n base = {**configs_base, **{\"data\": data_configs}, **inference_configs}\n \n # Deep update with model-specific params\n model_spec = model_configs[model_name]\n for k, v in model_spec.items():\n if isinstance(v, dict) and k in base:\n base[k].update(v)\n else:\n base[k] = v\n\n # CLI-style arguments to override defaults\n arg_str = (\n f\"--model_name {model_name} \"\n f\"--input_json_path {input_json_path} \"\n f\"--dump_dir {dump_dir} \"\n f\"--use_msa false \"\n f\"--use_template false \"\n f\"--use_rna_msa true \"\n f\"--sample_diffusion.N_sample 5 \"\n f\"--seeds 42\"\n )\n return parse_configs(configs=base, arg_str=arg_str, fill_required_with_null=True)","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:16:44.776586Z","iopub.status.busy":"2026-03-27T18:16:44.776116Z","iopub.status.idle":"2026-03-27T18:16:45.165046Z","shell.execute_reply":"2026-03-27T18:16:45.164251Z"},"papermill":{"duration":0.415525,"end_time":"2026-03-27T18:16:45.167055","exception":false,"start_time":"2026-03-27T18:16:44.75153","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import biotite.structure.io.pdb as pdb\nfrom pathlib import Path\n\ndef save_protenix_output_to_pdb(atom_array, prediction_tensor, output_pdb_path):\n \"\"\"\n Saves Protenix inference samples to a specific PDB file path.\n \"\"\"\n out_file = Path(output_pdb_path)\n # Ensure the parent directory (the target_id folder) exists\n out_file.parent.mkdir(parents=True, exist_ok=True)\n \n # Convert coordinates [N_samples, N_atoms, 3] -> use the first sample\n coords_np = prediction_tensor.detach().cpu().numpy()\n \n # Create a copy of the metadata and update coordinates\n sample_struct = atom_array.copy()\n sample_struct.coord = coords_np[0]\n \n # Write to PDB\n pdb_file = pdb.PDBFile()\n pdb_file.set_structure(sample_struct)\n pdb_file.write(str(out_file))\n \n print(f\"Successfully saved: {out_file}\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:16:45.215052Z","iopub.status.busy":"2026-03-27T18:16:45.214126Z","iopub.status.idle":"2026-03-27T18:16:46.321495Z","shell.execute_reply":"2026-03-27T18:16:46.320553Z"},"papermill":{"duration":1.133132,"end_time":"2026-03-27T18:16:46.323431","exception":false,"start_time":"2026-03-27T18:16:45.190299","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import os\n# import sys\n\n# # 1. Define where your wheels are (UPDATE THIS PATH with your actual dataset path)\n# OFFLINE_DATASET_PATH = \"/kaggle/input/datasets/khushvaghela/gcosrna/python_packages/kaggle/working/python_packages\" # <--- Fix this!\n\n# # 2. Define a target folder to install them into so they don't conflict with system libs\n# LIB_PATH = \"/kaggle/working/custom_lib\"\n# os.makedirs(LIB_PATH, exist_ok=True)\n\n# # 3. Install specifically to that folder\n# # We use scipy==1.11.4 because 1.13+ is what's causing the 'numpy.strings' error\n# !pip install --no-index --find-links={OFFLINE_DATASET_PATH} \\\n# --target={LIB_PATH} \\\n# \"numpy==1.26.4\" \\\n# \"pandas==2.1.4\" \\\n# \"scikit-learn==1.3.2\" \\\n# \"scipy==1.11.4\"\n\n# # 4. FORCE Python to use these versions by putting this folder at the START of the path\n# sys.path.insert(0, LIB_PATH)\n\n# print(\"โœ… Custom libraries installed and prioritized.\")","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:16:46.372424Z","iopub.status.busy":"2026-03-27T18:16:46.371397Z","iopub.status.idle":"2026-03-27T18:16:46.376593Z","shell.execute_reply":"2026-03-27T18:16:46.375565Z"},"papermill":{"duration":0.031433,"end_time":"2026-03-27T18:16:46.378335","exception":false,"start_time":"2026-03-27T18:16:46.346902","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import json\nimport torch\nimport os\nimport shutil\nfrom pathlib import Path\nfrom runner.inference import InferenceRunner, update_gpu_compatible_configs\nfrom protenix.data.inference.infer_dataloader import InferenceDataset\n\n# --- 1. GLOBAL INITIALIZATION (Run this once at the start) ---\n# We initialize with dummy values just to get the model into VRAM\nDUMMY_OUT = \"/kaggle/working/protenix_init\"\nos.makedirs(DUMMY_OUT, exist_ok=True)\n\nprint(\"๐Ÿš€ Loading Protenix weights into GPU (Global Init)...\")\nglobal_cfg = build_protenix_configs(None, DUMMY_OUT) \nglobal_cfg = update_gpu_compatible_configs(global_cfg)\nshared_runner = InferenceRunner(global_cfg)\nprint(\"โœ… Protenix ready for high-speed inference.\")\n\n\n# --- 2. FAST INFERENCE FUNCTION ---\nimport json\nimport os\nimport tempfile\n\ndef run_protenix_fast(sequence, target_id, mol_type, output_pdb_path, runner):\n \"\"\"\n Predicts structure for RNA, DNA, or Protein using the correct JSON schema.\n \"\"\"\n # Map our internal ctype to Protenix JSON keys\n type_map = {\n \"PROTEIN\": \"proteinChain\",\n \"DNA\": \"dnaSequence\",\n \"RNA\": \"rnaSequence\"\n }\n json_key = type_map.get(mol_type, \"rnaSequence\")\n\n with tempfile.NamedTemporaryFile(mode='w', suffix='.json', delete=False) as tf:\n input_data = [{\n \"name\": target_id,\n \"sequences\": [{json_key: {\"sequence\": sequence, \"count\": 1}}]\n }]\n json.dump(input_data, tf)\n temp_json_path = tf.name\n\n try:\n target_cfg = build_protenix_configs(temp_json_path, \"/tmp\")\n dataset = InferenceDataset(target_cfg)\n \n if not dataset or len(dataset) == 0:\n return None\n \n data, atom_array, err = dataset[0]\n if err:\n print(f\"โŒ Featurization failed for {target_id}: {err}\")\n return None\n \n prediction = runner.predict(data)\n save_protenix_output_to_pdb(atom_array, prediction[\"coordinate\"], str(output_pdb_path))\n \n return str(output_pdb_path)\n \n finally:\n if os.path.exists(temp_json_path):\n os.remove(temp_json_path)","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:16:46.425984Z","iopub.status.busy":"2026-03-27T18:16:46.425572Z","iopub.status.idle":"2026-03-27T18:18:49.621663Z","shell.execute_reply":"2026-03-27T18:18:49.620537Z"},"papermill":{"duration":123.222376,"end_time":"2026-03-27T18:18:49.623644","exception":false,"start_time":"2026-03-27T18:16:46.401268","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport subprocess\nimport json\nimport torch\nimport shutil\nimport time\nimport os\nimport re\nfrom pathlib import Path\nfrom Bio.PDB import PDBParser, PDBIO, Select\nfrom Bio.Align import PairwiseAligner","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:18:49.674203Z","iopub.status.busy":"2026-03-27T18:18:49.67343Z","iopub.status.idle":"2026-03-27T18:18:49.678646Z","shell.execute_reply":"2026-03-27T18:18:49.6779Z"},"papermill":{"duration":0.032212,"end_time":"2026-03-27T18:18:49.680362","exception":false,"start_time":"2026-03-27T18:18:49.64815","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import re\nfrom Bio.PDB import PDBParser\n\ndef extract_all_atoms_from_pdb(pdb_file):\n parser = PDBParser(QUIET=True)\n structure = parser.get_structure(\"tmp\", pdb_file)\n data = []\n for model in structure:\n for chain in model:\n for residue in chain:\n resname = residue.get_resname().strip()\n resid = residue.get_id()[1]\n for atom in residue:\n coord = atom.get_coord()\n data.append({'resname': resname, 'resid': resid,\n 'atom_name': atom.get_name(),\n 'x': coord[0], 'y': coord[1], 'z': coord[2]})\n return data\n\n# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# FIXED chain parsing\n# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# The competition's `sequence` column is the AUTHORITATIVE, scored target: it is\n# the concatenation, in `all_sequences` entry order, of every polymer ENTITY\n# repeated once per author-chain copy. The number of copies of an entity equals\n# the number of `[auth X]` chain IDs listed in that entity's FASTA header โ€” NOT\n# the `stoichiometry` column (which the old identify_chains mis-parsed, producing\n# the wrong residue count for every multi-copy target and thus an invalid\n# submission). We rebuild the chains from the header and then VERIFY the RNA\n# concatenation equals `sequence`; the caller falls back to a single-chain model\n# if verification ever fails.\n\ndef parse_entities(all_sequences_str):\n \"\"\"Ordered list of unique entities: {'seq', 'n_chains', 'is_rna'}.\"\"\"\n entities = []\n for entry in str(all_sequences_str).strip().split('>'):\n if not entry.strip():\n continue\n parts = entry.split('|')\n n_chains = 0\n for p in parts:\n if 'Chain' in p:\n auth = re.findall(r'\\[auth\\s+([A-Za-z0-9]+)\\]', p)\n if auth:\n n_chains = len(auth)\n else:\n m = re.search(r'Chains?\\s+([A-Za-z0-9,\\s]+)', p)\n if m:\n n_chains = len([c for c in re.split(r'[,\\s]+', m.group(1)) if c])\n break\n if n_chains == 0:\n n_chains = 1\n seq = re.sub(r'[^A-Za-z]', '', parts[-1]).upper()\n if not seq:\n continue\n is_rna = all(c in 'ACGUN' for c in seq)\n entities.append({'seq': seq, 'n_chains': n_chains, 'is_rna': is_rna})\n return entities\n\ndef build_segments(all_sequences_str, sequence=None):\n \"\"\"Ordered list of (chain_seq, is_rna), one item per author-chain copy.\n If `sequence` is given and the RNA concatenation does not match it, we\n fall back to a single RNA chain equal to `sequence` (guarantees the\n submission always has exactly len(sequence) residues in the right order).\"\"\"\n ents = parse_entities(all_sequences_str)\n segs = []\n for e in ents:\n for _ in range(e['n_chains']):\n segs.append((e['seq'], e['is_rna']))\n if sequence is not None:\n rna_concat = ''.join(s for s, isrna in segs if isrna)\n if rna_concat != str(sequence):\n segs = [(str(sequence), True)]\n return segs\n\n# โ”€โ”€ legacy (kept for reference; no longer used by the main loop) โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\ndef identify_chains(all_sequences_str, stoichiometry_str):\n \"\"\"DEPRECATED: superseded by build_segments(). The old stoichiometry-based\n copy counting double-counted chains for every multi-copy target.\"\"\"\n segs = build_segments(all_sequences_str, None)\n return [{'chain': None, 'auth': None,\n 'type': 'RNA' if isrna else 'PROTEIN', 'seq': s}\n for s, isrna in segs]\n","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:18:49.73067Z","iopub.status.busy":"2026-03-27T18:18:49.729972Z","iopub.status.idle":"2026-03-27T18:18:49.74225Z","shell.execute_reply":"2026-03-27T18:18:49.741216Z"},"papermill":{"duration":0.039367,"end_time":"2026-03-27T18:18:49.744111","exception":false,"start_time":"2026-03-27T18:18:49.704744","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def format_submission_row(target_id, refined_atoms):\n formatted_rows = []\n residue_map = {}\n\n # Standard 3-letter codes for proteins to avoid truncation\n amino_acids = [\"ALA\", \"ARG\", \"ASN\", \"ASP\", \"CYS\", \"GLN\", \"GLU\", \"GLY\", \"HIS\", \n \"ILE\", \"LEU\", \"LYS\", \"MET\", \"PHE\", \"PRO\", \"SER\", \"THR\", \"TRP\", \"TYR\", \"VAL\"]\n\n for a in refined_atoms:\n resid = a['resid']\n if resid not in residue_map:\n residue_map[resid] = {\n 'resname': a['resname'],\n 'coord': None\n }\n\n # โœ… FIX: Look for C1' (Nucleic) OR CA (Protein)\n atom_name = a['atom_name'].strip()\n if atom_name in [\"C1'\", \"C1*\", \"C1\", \"CA\"]:\n residue_map[resid]['coord'] = (a['x'], a['y'], a['z'])\n\n for resid in sorted(residue_map.keys()):\n res = residue_map[resid]\n resname = res['resname'].strip()\n coord = res['coord']\n\n if coord is None:\n # Re-enable the warning to help debug missing chain segments\n print(f\"โš ๏ธ Warning: No representative atom (C1'/CA) found for {target_id} residue {resid}.\")\n continue \n\n # โœ… FIX: Normalization that respects Proteins\n if resname in [\"A\", \"C\", \"G\", \"U\"]:\n resname_sub = resname\n elif resname.upper() in amino_acids:\n resname_sub = resname # Keep full 3-letter code for proteins\n elif resname.startswith(\"R\"): # RA, RC...\n resname_sub = resname[1]\n elif resname.startswith(\"D\"): # DNA\n resname_sub = resname\n else:\n resname_sub = resname[-1] # Fallback for single-letter formats\n\n row_entry = {\n 'ID': f\"{target_id}_{resid}\",\n 'resname': resname_sub,\n 'resid': resid\n }\n\n # Repeat coordinates for all 5 sample columns\n for i in range(1, 6):\n row_entry[f'x_{i}'], row_entry[f'y_{i}'], row_entry[f'z_{i}'] = coord\n\n formatted_rows.append(row_entry)\n\n print(f\"๐Ÿ“Š Formatted {len(formatted_rows)} residues (RNA + Protein) for {target_id}\")\n return formatted_rows","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:18:49.794796Z","iopub.status.busy":"2026-03-27T18:18:49.794215Z","iopub.status.idle":"2026-03-27T18:18:49.803339Z","shell.execute_reply":"2026-03-27T18:18:49.802591Z"},"papermill":{"duration":0.036583,"end_time":"2026-03-27T18:18:49.804896","exception":false,"start_time":"2026-03-27T18:18:49.768313","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nfrom Bio.Align import PairwiseAligner\n\ndef run_msa_homology_prediction(sequence, target_id, train_seqs, train_coords_dict, output_path):\n \"\"\"\n Template-based prediction for large molecules with coordinate mapping.\n Ensures connectivity and sequence length matching.\n \"\"\"\n print(f\"๐Ÿ” Homology Mapping: Processing {target_id}...\")\n \n # 1. Setup Aligner\n aligner = PairwiseAligner()\n aligner.match_score, aligner.mismatch_score = 2.0, -1.0\n aligner.open_gap_score, aligner.extend_gap_score = -0.5, -0.1\n aligner.mode = 'global'\n \n similarities = []\n for _, row in train_seqs.iterrows():\n t_id, t_seq = row['target_id'], row['sequence']\n # Length filter to optimize performance\n if abs(len(t_seq) - len(sequence)) > len(sequence) * 0.5: \n continue\n \n # We check if coordinates exist before aligning to save time\n if t_id not in train_coords_dict:\n continue\n \n score = aligner.score(sequence, t_seq)\n similarities.append((t_id, t_seq, score))\n \n # Sort by alignment score\n similarities.sort(key=lambda x: x[2], reverse=True)\n \n if not similarities:\n print(f\"โš ๏ธ No suitable template found for {target_id}. Using dummy linear model.\")\n dummy_coords = [np.array([i * 3.8, 0, 0]) for i in range(len(sequence))]\n save_coords_to_pdb(sequence, dummy_coords, target_id, output_path)\n return output_path\n\n # 2. Map Coordinates from Best Template\n best_t_id, best_t_seq, _ = similarities[0]\n t_coords = train_coords_dict[best_t_id]\n \n alignment = aligner.align(sequence, best_t_seq)[0]\n q_aln, t_aln = str(alignment[0]), str(alignment[1])\n \n new_coords = []\n t_idx = 0\n \n for q_char, t_char in zip(q_aln, t_aln):\n if q_char != \"-\":\n if t_char != \"-\" and t_idx < len(t_coords):\n # Use coordinate from template\n new_coords.append(t_coords[t_idx])\n else:\n # HEURISTIC: Maintain 3.8A backbone connectivity for gaps or overhangs\n prev = new_coords[-1] if new_coords else np.array([0.0, 0.0, 0.0])\n new_coords.append(prev + np.array([3.8, 0.0, 0.0]))\n \n if t_char != \"-\":\n t_idx += 1\n\n # 3. Final Save\n save_coords_to_pdb(sequence, new_coords[:len(sequence)], target_id, output_path)\n print(f\"โœ… Homology model saved using template {best_t_id}\")\n return output_path","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:18:49.853903Z","iopub.status.busy":"2026-03-27T18:18:49.853502Z","iopub.status.idle":"2026-03-27T18:18:49.864276Z","shell.execute_reply":"2026-03-27T18:18:49.863511Z"},"papermill":{"duration":0.037309,"end_time":"2026-03-27T18:18:49.865985","exception":false,"start_time":"2026-03-27T18:18:49.828676","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# SUBMISSION-AUTHORITATIVE ORCHESTRATION\n# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n# The submission MUST contain exactly len(sequence) rows per target, IDs\n# {target_id}_{1..N}, resname = sequence[i]. So we build the output directly from\n# the `sequence` column: predict one 3D monomer per unique RNA chain, then lay the\n# chains into a full-length coordinate array in the exact order that reconstructs\n# `sequence`. Any failure degrades to a valid placeholder chain.\n#\n# HOW SCORING SHAPES THIS DESIGN (Stanford RNA 3D Folding 2):\n# โ€ข Scored with TM-score on C1' atoms only; the prediction is superimposed onto\n# the reference by US-align (a global rotation+translation, sequence-independent),\n# best-of-5, averaged over targets.\n# โ€ข Because US-align applies ONE global transform to the whole prediction, a rigid\n# move of the entire structure does not change the score. Consequently:\n# - For a SINGLE RNA chain, only its internal conformation matters.\n# - For MULTIPLE RNA chains (homo/hetero multimers: 9G4P, 9G4Q, 9WHV, 9I9W,\n# 9JGM, 9MME, 9ZCC), the RELATIVE placement of the chains matters, since one\n# global transform must fit them all โ€” this is where docking helps.\n# - A non-scored DNA/protein partner does NOT change the RNA's internal fold\n# (HDOCK is rigid-body; DRfold2 folds the RNA monomer independently), so\n# docking a partner is score-NEUTRAL for the RNA C1' residues. It can still\n# be used as a scaffold if you want (ENABLE_PARTNER_PROTENIX below).\n#\n# NOTE ON DRFOLD_MAX_LEN: this decision is now made ONCE PER TARGET, from the\n# length of the overall/full sequence (row.sequence), NOT per individual chain.\n# A target whose full sequence is > DRFOLD_MAX_LEN nt uses template matching for\n# every one of its chains; a target <= DRFOLD_MAX_LEN nt uses DRfold2 (falling\n# back to template matching per-chain only if DRfold2 errors/returns nothing).\n\nDRFOLD_MAX_LEN = 150 # overall target-sequence length above which we skip\n # DRfold2 for every chain and use template matching only\nINTER_CHAIN_OFFSET = 60.0 # ร… translation between assembled chains (avoid overlap)\n\n# --- assembly strategy for MULTI-CHAIN targets ---\n# False -> fast, always-valid: place each chain's monomer with a fixed offset\n# (correct internal folds; relative placement is only a guess).\n# True -> dock chains into one complex with HDOCK to get real relative placement,\n# then extract ONLY the RNA C1' atoms (in order) for the rows. Slower and\n# depends on the HDOCK binaries; falls back to offset assembly on any\n# failure or residue-count mismatch. (Requires setup_hdock() to have run.)\nENABLE_COMPLEX_DOCKING = True\n# When docking, also predict DNA/protein partner chains with Protenix and include\n# them in the docking scaffold. Score-neutral for the RNA rows (see note above);\n# only affects the docking geometry. Needs a loaded `shared_runner`.\nENABLE_PARTNER_PROTENIX = False\n\n_MONOMER_CACHE = {} # cache monomer coords by chain sequence\n\ndef _straight_chain(L):\n return np.array([[i * 3.8, 0.0, 0.0] for i in range(L)], dtype=float)\n\ndef _to_coords_array(coords, L):\n \"\"\"Coerce anything to a clean (L,3) float array, padding along +x if short.\"\"\"\n try:\n a = np.asarray(coords, dtype=float)\n except Exception:\n a = np.empty((0, 3))\n if a.ndim != 2 or a.shape[1] != 3:\n a = _straight_chain(L)\n a = a[:L]\n while a.shape[0] < L:\n last = a[-1] if a.shape[0] > 0 else np.array([0.0, 0.0, 0.0])\n a = np.vstack([a, last + np.array([3.8, 0.0, 0.0])])\n return a\n\ndef predict_monomer_coords(chain_seq, chain_id, original_target_id, msa_dir,\n train_seqs, train_coords_dict, use_drfold2):\n \"\"\"Predict (L,3) C1' coords for ONE RNA chain: DRfold2 -> template -> straight\n chain, then an optional gCoSRNA coaxial refinement. Cached per sequence.\n Note: MSA lookups use `original_target_id` (the bare target, e.g. 8ZNQ),\n never the per-chain label (8ZNQ_c0).\n\n `use_drfold2` is decided ONCE per target in run_main_loop (based on the\n length of the overall/full sequence, not this chain's own length), and\n passed down through assemble_full_coords / _assemble_offset / _assemble_docked\n so every chain belonging to the same target is treated consistently.\"\"\"\n if chain_seq in _MONOMER_CACHE:\n return _MONOMER_CACHE[chain_seq].copy()\n\n L = len(chain_seq)\n coords = None\n\n if use_drfold2:\n try:\n preds = predict_rna_structures_drfold2(\n chain_seq, chain_id, original_target_id, msa_dir,\n train_seqs, train_coords_dict)\n if preds:\n coords = np.asarray(preds[0], dtype=float)\n except Exception as e:\n print(f\" โš ๏ธ DRfold2 failed for {chain_id}: {e}\")\n\n if coords is None:\n try:\n preds = predict_rna_structures(chain_seq, chain_id, train_seqs, train_coords_dict)\n if preds:\n coords = np.asarray(preds[0], dtype=float)\n except Exception as e:\n print(f\" โš ๏ธ Template matching failed for {chain_id}: {e}\")\n\n if coords is None or coords.size == 0:\n print(f\" โ†ช๏ธ {chain_id}: using straight-chain placeholder ({L} nt).\")\n coords = _straight_chain(L)\n\n coords = _to_coords_array(coords, L)\n\n try:\n ss = get_consensus_ss(chain_seq, chain_id, msa_dir, original_target_id)\n cx = get_coaxial_stem_pairs(chain_seq, ss)\n if cx:\n print(f\" ๐Ÿ“ {chain_id}: applying {len(cx)} coaxial stack(s).\")\n coords = _to_coords_array(apply_coaxial_alignment(coords, chain_seq, cx), L)\n except Exception as e:\n print(f\" โš ๏ธ Coaxial refinement skipped for {chain_id}: {e}\")\n\n _MONOMER_CACHE[chain_seq] = coords.copy()\n return coords\n\n# โ”€โ”€ OFFSET assembly (default): fast + always valid โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\ndef _assemble_offset(rna_chains, target_id, msa_dir, train_seqs, train_coords_dict,\n use_drfold2):\n parts = []\n for idx, cseq in enumerate(rna_chains):\n mc = predict_monomer_coords(cseq, f\"{target_id}_c{idx}\", target_id,\n msa_dir, train_seqs, train_coords_dict, use_drfold2)\n parts.append(mc + np.array([idx * INTER_CHAIN_OFFSET, 0.0, 0.0]))\n return np.vstack(parts) if parts else None\n\n# โ”€โ”€ DOCKED assembly (optional): real relative placement via HDOCK โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€\n_NUC = {\"A\", \"C\", \"G\", \"U\", \"RA\", \"RC\", \"RG\", \"RU\", \"N\"}\ndef _is_rna_res(resname):\n r = resname.strip().upper()\n return r in _NUC or (len(r) == 1 and r in \"ACGU\") or (r.startswith(\"R\") and r[1:] in {\"A\",\"C\",\"G\",\"U\"})\n\ndef _extract_rna_c1_in_order(pdb_path):\n \"\"\"Return C1' coords of every RNA residue in the complex, in file order.\"\"\"\n atoms = extract_all_atoms_from_pdb(str(pdb_path))\n out, seen = [], set()\n for a in atoms:\n if not _is_rna_res(a['resname']):\n continue\n an = a['atom_name'].strip()\n if an in (\"C1'\", \"C1*\", \"C1\"):\n out.append((a['x'], a['y'], a['z']))\n return np.array(out, dtype=float) if out else np.empty((0, 3))\n\ndef _assemble_docked(segs, target_id, msa_dir, train_seqs, train_coords_dict,\n shared_runner, expected_len, use_drfold2):\n \"\"\"Predict every chain, dock them in sequence order into one complex, then\n pull back ONLY the RNA C1' atoms (in order). Falls back (returns None) on\n any failure or if the RNA residue count != expected_len.\"\"\"\n try:\n target_dir = Path(target_id); target_dir.mkdir(exist_ok=True)\n chain_pdbs = []\n for idx, (cseq, is_rna) in enumerate(segs):\n cid = f\"{target_id}_c{idx}\"\n pdb_path = target_dir / f\"{cid}.pdb\"\n if is_rna:\n mc = predict_monomer_coords(cseq, cid, target_id, msa_dir,\n train_seqs, train_coords_dict, use_drfold2)\n save_coords_to_pdb(cseq, mc, cid, str(pdb_path))\n chain_pdbs.append(str(pdb_path))\n elif ENABLE_PARTNER_PROTENIX and shared_runner is not None:\n mol = \"DNA\" if (\"T\" in cseq and \"U\" not in cseq) else \"PROTEIN\"\n run_protenix_fast(sequence=cseq, target_id=cid, mol_type=mol,\n output_pdb_path=str(pdb_path), runner=shared_runner)\n if pdb_path.exists():\n chain_pdbs.append(str(pdb_path))\n # else: partner chain skipped (score-neutral for RNA rows)\n\n if not chain_pdbs:\n return None\n\n current = chain_pdbs[0]\n for j in range(1, len(chain_pdbs)):\n docked = run_hdock_assembly(receptor_pdb=current, ligand_pdb=chain_pdbs[j],\n output_folder=target_dir,\n target_id=f\"{target_id}_dock_{j}\")\n if docked:\n current = docked\n else:\n print(f\" โš ๏ธ HDOCK step {j} failed; falling back to offset assembly.\")\n return None\n\n coords = _extract_rna_c1_in_order(current)\n if coords.shape[0] != expected_len:\n print(f\" โš ๏ธ docked RNA residue count {coords.shape[0]} != {expected_len}; \"\n f\"falling back to offset assembly.\")\n return None\n return coords\n except Exception as e:\n print(f\" โš ๏ธ docked assembly error ({e}); falling back to offset assembly.\")\n return None\n\ndef assemble_full_coords(segs, target_id, msa_dir, train_seqs, train_coords_dict,\n shared_runner=None, use_drfold2=True):\n \"\"\"Produce (len(sequence),3) coords. Uses docking if enabled (with fallback),\n otherwise the fast offset assembly. Returns None only if no RNA chains.\n\n `use_drfold2` reflects the per-target decision (based on OVERALL sequence\n length) made once in run_main_loop, and is forwarded to every chain.\"\"\"\n rna_chains = [s for s, is_rna in segs if is_rna]\n if not rna_chains:\n return None\n expected_len = sum(len(s) for s in rna_chains)\n\n if ENABLE_COMPLEX_DOCKING and len(rna_chains) > 1:\n docked = _assemble_docked(segs, target_id, msa_dir, train_seqs,\n train_coords_dict, shared_runner, expected_len,\n use_drfold2)\n if docked is not None:\n return docked\n return _assemble_offset(rna_chains, target_id, msa_dir, train_seqs,\n train_coords_dict, use_drfold2)\n\ndef rows_for_target(target_id, sequence, coords):\n \"\"\"Emit exactly len(sequence) rows: IDs 1..N, resname=sequence[i], coords in\n all 5 model columns.\"\"\"\n coords = _to_coords_array(coords, len(sequence))\n rows = []\n for i, ch in enumerate(sequence):\n x, y, z = float(coords[i][0]), float(coords[i][1]), float(coords[i][2])\n entry = {'ID': f\"{target_id}_{i+1}\", 'resname': ch, 'resid': i + 1}\n for k in range(1, 6):\n entry[f'x_{k}'], entry[f'y_{k}'], entry[f'z_{k}'] = x, y, z\n rows.append(entry)\n return rows\n\ndef _whole_structure_template_coords(target_id, sequence, msa_dir, train_seqs, train_coords_dict):\n \"\"\"Large targets (use_drfold2=False): run template/homology matching ONCE on\n the full, concatenated sequence โ€” not per chain โ€” so a matched homolog's\n relative chain placement is used as-is instead of being reconstructed via\n independent per-chain predictions + HDOCK docking.\"\"\"\n target_dir = Path(target_id); target_dir.mkdir(exist_ok=True)\n homology_pdb = target_dir / f\"{target_id}_homology.pdb\"\n try:\n run_msa_homology_prediction(\n sequence, target_id, train_seqs, train_coords_dict, str(homology_pdb))\n except Exception as e:\n print(f\" โš ๏ธ whole-structure template matching failed for {target_id}: {e}\")\n return None\n\n if not homology_pdb.exists():\n return None\n\n coords = _extract_rna_c1_in_order(str(homology_pdb))\n if coords.shape[0] != len(sequence):\n print(f\" โš ๏ธ homology RNA residue count {coords.shape[0]} != {len(sequence)}; \"\n f\"discarding.\")\n return None\n return coords\n\n\ndef run_main_loop(input_csv_path, msa_dir, train_seqs, train_coords_dict,\n shared_runner=None, output_csv=\"submission.csv\"):\n \"\"\"Full pipeline. Every target yields a schema-correct block of rows; a\n failing target degrades to a placeholder rather than aborting the run.\"\"\"\n df = pd.read_csv(input_csv_path)\n all_rows = []\n start = time.time()\n print(f\"๐Ÿ Starting inference on {len(df)} targets...\")\n\n for pos, row in enumerate(df.itertuples(), 1):\n target_id = str(row.target_id)\n sequence = str(row.sequence)\n L = len(sequence)\n\n # Decided ONCE per target, from the OVERALL sequence length:\n # L <= DRFOLD_MAX_LEN -> per-chain DRfold2 + optional HDOCK docking\n # L > DRFOLD_MAX_LEN -> template matching on the WHOLE structure,\n # no chain splitting, no docking\n use_drfold2 = L <= DRFOLD_MAX_LEN\n\n print(f\"\\n[{pos}/{len(df)}] ๐Ÿš€ {target_id} (len {L}) \"\n f\"[{'DRfold2+dock' if use_drfold2 else 'whole-structure template'}]\")\n\n if time.time() - start > 8 * 3600:\n print(\"๐Ÿ•’ Time budget reached; filling the rest with placeholders.\")\n all_rows.extend(rows_for_target(target_id, sequence, _straight_chain(L)))\n continue\n\n try:\n if use_drfold2:\n segs = build_segments(row.all_sequences, sequence) # verified vs sequence\n coords = assemble_full_coords(segs, target_id, msa_dir,\n train_seqs, train_coords_dict, shared_runner,\n use_drfold2=True)\n else:\n coords = _whole_structure_template_coords(\n target_id, sequence, msa_dir, train_seqs, train_coords_dict)\n\n if coords is None or len(coords) != L:\n print(f\" โš ๏ธ assembly length {None if coords is None else len(coords)} \"\n f\"!= {L}; using placeholder.\")\n coords = _straight_chain(L)\n except Exception as e:\n print(f\" โŒ {target_id} failed ({e}); using placeholder.\")\n coords = _straight_chain(L)\n\n all_rows.extend(rows_for_target(target_id, sequence, coords))\n print(f\" โœ… {target_id}: wrote {L} residues.\")\n if torch.cuda.is_available():\n torch.cuda.empty_cache()\n\n cols = ['ID', 'resname', 'resid',\n 'x_1', 'y_1', 'z_1', 'x_2', 'y_2', 'z_2',\n 'x_3', 'y_3', 'z_3', 'x_4', 'y_4', 'z_4',\n 'x_5', 'y_5', 'z_5']\n submission_df = pd.DataFrame(all_rows)[cols]\n submission_df.to_csv(output_csv, index=False)\n print(f\"\\n๐Ÿ“ฆ Wrote {len(submission_df)} rows for \"\n f\"{submission_df['ID'].str.rsplit('_', n=1).str[0].nunique()} targets -> {output_csv}\")\n print(f\"โฑ๏ธ Runtime: {(time.time()-start)/3600:.2f} h\")\n return submission_df","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:18:49.91729Z","iopub.status.busy":"2026-03-27T18:18:49.916321Z","iopub.status.idle":"2026-03-27T18:18:49.931206Z","shell.execute_reply":"2026-03-27T18:18:49.930321Z"},"papermill":{"duration":0.042053,"end_time":"2026-03-27T18:18:49.933073","exception":false,"start_time":"2026-03-27T18:18:49.89102","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"cell_type":"code","source":"input_csv_path = \"/kaggle/input/competitions/stanford-rna-3d-folding-2/test_sequences.csv\"\nmsa_dir = \"/kaggle/input/competitions/stanford-rna-3d-folding-2/MSA\"\n\n# shared_runner (Protenix) is optional now: the scored `sequence` is RNA-only, so\n# protein prediction is not needed to produce a valid submission. Pass it if present.\n_runner = shared_runner if \"shared_runner\" in globals() else None\n\nrun_main_loop(input_csv_path=input_csv_path,\n msa_dir=msa_dir,\n train_seqs=train_seqs,\n train_coords_dict=train_coords_dict,\n shared_runner=_runner)\n","metadata":{"execution":{"iopub.execute_input":"2026-03-27T18:18:49.981656Z","iopub.status.busy":"2026-03-27T18:18:49.981193Z","iopub.status.idle":"2026-03-27T21:30:39.153975Z","shell.execute_reply":"2026-03-27T21:30:39.152551Z"},"papermill":{"duration":11509.199477,"end_time":"2026-03-27T21:30:39.156034","exception":false,"start_time":"2026-03-27T18:18:49.956557","status":"completed"},"tags":[]},"outputs":[],"execution_count":null}]}