{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Dataset Loading","metadata":{}},{"cell_type":"code","source":"# RDKit is not part of the current Kaggle base image, so install it if it is missing.\n# (Needs \"Internet\" switched on in the notebook settings.)\nimport importlib.util, subprocess, sys\nif importlib.util.find_spec('rdkit') is None:\n    print('installing rdkit ...')\n    subprocess.run([sys.executable, '-m', 'pip', 'install', '-q', 'rdkit'], check=True)\n\nimport numpy as np\nimport pandas as pd\nimport duckdb\nimport time\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.metrics import average_precision_score, roc_auc_score\n\nfrom rdkit import Chem, RDLogger\nfrom rdkit.Chem import Descriptors, rdFingerprintGenerator\n\nRDLogger.DisableLog('rdApp.*')   # silence RDKit SMILES parsing chatter\npd.set_option('display.max_columns', 50)\nsns.set_theme(style='whitegrid')\n\nprint('Imports OK')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:09:26.080785Z","iopub.execute_input":"2026-08-13T06:09:26.081181Z","iopub.status.idle":"2026-08-13T06:09:36.110508Z","shell.execute_reply.started":"2026-08-13T06:09:26.081147Z","shell.execute_reply":"2026-08-13T06:09:36.109684Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_path = '/kaggle/input/competitions/leash-BELKA/train.parquet'\ntest_path = '/kaggle/input/competitions/leash-BELKA/test.parquet'\n\n# Load a class-balanced sample: 90,000 non-binders + 90,000 binders\ncon = duckdb.connect()\nquery = f\"\"\"\n(SELECT * FROM parquet_scan('{train_path}') WHERE binds = 0 ORDER BY random() LIMIT 90000)\nUNION ALL\n(SELECT * FROM parquet_scan('{train_path}') WHERE binds = 1 ORDER BY random() LIMIT 90000)\n\"\"\"\ndf = con.query(query).df()\ncon.close()\n\nprint('Loaded balanced sample:', df.shape)\ndf['binds'].value_counts()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:09:36.112548Z","iopub.execute_input":"2026-08-13T06:09:36.113078Z","iopub.status.idle":"2026-08-13T06:10:13.539792Z","shell.execute_reply.started":"2026-08-13T06:09:36.113048Z","shell.execute_reply":"2026-08-13T06:10:13.539077Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Basic EDA","metadata":{}},{"cell_type":"code","source":"print('Shape:', df.shape)\ndf.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:10:13.540826Z","iopub.execute_input":"2026-08-13T06:10:13.541274Z","iopub.status.idle":"2026-08-13T06:10:13.568265Z","shell.execute_reply.started":"2026-08-13T06:10:13.541248Z","shell.execute_reply":"2026-08-13T06:10:13.567484Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- EDA 1: structure & data quality of the sample ----------\nprint('Rows / columns :', df.shape)\n\nprint('\\nColumn dtypes:')\nprint(df.dtypes)\n\nprint('\\nMissing values per column:')\nprint(df.isna().sum())\n\nprint('\\nUnique values per column:')\nprint(df.nunique())\n\nprint('\\nDuplicated molecule_smiles      :', df['molecule_smiles'].duplicated().sum())\nprint('Duplicated (molecule, protein) :', df.duplicated(subset=['molecule_smiles', 'protein_name']).sum())\n\nprint('\\nProtein targets present:')\nprint(df['protein_name'].value_counts())\n\ndf.describe(include='all').T","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:10:13.569278Z","iopub.execute_input":"2026-08-13T06:10:13.569542Z","iopub.status.idle":"2026-08-13T06:10:14.255012Z","shell.execute_reply.started":"2026-08-13T06:10:13.569519Z","shell.execute_reply":"2026-08-13T06:10:14.254012Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- EDA 2: target balance, overall and per protein ----------\nprint('binds distribution (this sample was drawn 50/50 on purpose):')\nprint(df['binds'].value_counts(normalize=True).rename('share'))\n\nper_protein = (df.groupby('protein_name')['binds']\n                 .agg(n='size', binders='sum', bind_rate='mean')\n                 .sort_values('bind_rate', ascending=False))\nprint('\\nBinding rate per protein target:')\nprint(per_protein)\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 4))\nsns.countplot(x='binds', data=df, ax=axes[0])\naxes[0].set_title('Class balance of the loaded sample')\nsns.barplot(x=per_protein.index, y=per_protein['bind_rate'], ax=axes[1])\naxes[1].set_title('Binding rate by protein')\naxes[1].set_ylabel('P(binds = 1)')\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:10:14.257444Z","iopub.execute_input":"2026-08-13T06:10:14.257881Z","iopub.status.idle":"2026-08-13T06:10:14.899701Z","shell.execute_reply.started":"2026-08-13T06:10:14.257849Z","shell.execute_reply":"2026-08-13T06:10:14.898473Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- EDA 3: physico-chemical profile, binders vs non-binders ----------\n# RDKit descriptors are expensive, so profile a random subsample.\nsub = df.sample(6000, random_state=42).copy()\n\ndef describe_molecule(smiles):\n    mol = Chem.MolFromSmiles(smiles)\n    if mol is None:\n        return pd.Series([np.nan] * 6)\n    return pd.Series([\n        Descriptors.MolWt(mol),\n        Descriptors.MolLogP(mol),\n        Descriptors.TPSA(mol),\n        Descriptors.NumHDonors(mol),\n        Descriptors.NumHAcceptors(mol),\n        Descriptors.NumRotatableBonds(mol),\n    ])\n\ndesc_cols = ['mol_wt', 'logp', 'tpsa', 'h_donors', 'h_acceptors', 'rot_bonds']\nsub[desc_cols] = sub['molecule_smiles'].apply(describe_molecule)\nsub['smiles_len'] = sub['molecule_smiles'].str.len()\n\nprint('Unparsable SMILES in subsample:', int(sub['mol_wt'].isna().sum()))\nprint('\\nMean descriptor value by class:')\nprint(sub.groupby('binds')[desc_cols + ['smiles_len']].mean().T)\n\nfig, axes = plt.subplots(2, 3, figsize=(15, 7))\nfor ax, col in zip(axes.ravel(), desc_cols):\n    sns.kdeplot(data=sub, x=col, hue='binds', common_norm=False, fill=True, ax=ax)\n    ax.set_title(col)\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:10:14.900791Z","iopub.execute_input":"2026-08-13T06:10:14.901207Z","iopub.status.idle":"2026-08-13T06:10:25.646531Z","shell.execute_reply.started":"2026-08-13T06:10:14.901177Z","shell.execute_reply":"2026-08-13T06:10:25.645444Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- EDA 4: building-block analysis ----------\nbb_cols = ['buildingblock1_smiles', 'buildingblock2_smiles', 'buildingblock3_smiles']\nprint('Distinct building blocks per position:')\nprint(df[bb_cols].nunique())\n\nfig, axes = plt.subplots(1, 3, figsize=(16, 4))\nfor ax, col in zip(axes, bb_cols):\n    rate = (df.groupby(col)['binds']\n              .agg(n='size', bind_rate='mean')\n              .query('n >= 50')\n              .sort_values('bind_rate', ascending=False))\n    sns.histplot(rate['bind_rate'], bins=30, ax=ax)\n    ax.set_title(col.replace('_smiles', '') + '\\nbind-rate spread (blocks with n>=50)')\n    ax.set_xlabel('P(binds = 1)')\n    print('\\nTop 5 blocks by bind rate for', col)\n    print(rate.head(5))\nplt.tight_layout()\nplt.show()\n\n# How many building blocks are shared between binders and non-binders?\nfor col in bb_cols:\n    pos = set(df.loc[df['binds'] == 1, col].unique())\n    neg = set(df.loc[df['binds'] == 0, col].unique())\n    print(f'{col}: binder-only={len(pos - neg)}, non-binder-only={len(neg - pos)}, shared={len(pos & neg)}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:10:25.647910Z","iopub.execute_input":"2026-08-13T06:10:25.648315Z","iopub.status.idle":"2026-08-13T06:10:26.449621Z","shell.execute_reply.started":"2026-08-13T06:10:25.648276Z","shell.execute_reply":"2026-08-13T06:10:26.448614Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- EDA 4b: RDKit descriptors per building block (bb1 / bb2 / bb3) ----------\nfrom rdkit.Chem import rdMolDescriptors\n\nbb_cols = ['buildingblock1_smiles', 'buildingblock2_smiles', 'buildingblock3_smiles']\n\nBB_DESCRIPTORS = {\n    'mol_wt':        Descriptors.MolWt,\n    'logp':          Descriptors.MolLogP,\n    'tpsa':          Descriptors.TPSA,\n    'h_donors':      Descriptors.NumHDonors,\n    'h_acceptors':   Descriptors.NumHAcceptors,\n    'rot_bonds':     Descriptors.NumRotatableBonds,\n    'heavy_atoms':   Descriptors.HeavyAtomCount,\n    'rings':         rdMolDescriptors.CalcNumRings,\n    'arom_rings':    rdMolDescriptors.CalcNumAromaticRings,\n    'frac_csp3':     Descriptors.FractionCSP3,\n    'hetero_atoms':  rdMolDescriptors.CalcNumHeteroatoms,\n    'n_count':       lambda m: sum(a.GetSymbol() == 'N' for a in m.GetAtoms()),\n    'o_count':       lambda m: sum(a.GetSymbol() == 'O' for a in m.GetAtoms()),\n    'halogens':      lambda m: sum(a.GetSymbol() in ('F', 'Cl', 'Br', 'I') for a in m.GetAtoms()),\n    'stereocenters': lambda m: len(Chem.FindMolChiralCenters(m, includeUnassigned=True)),\n}\n\n# There are only ~1,800 distinct blocks in total, so describe each one once.\ntables = {}\nfor col in bb_cols:\n    uniq = pd.Series(df[col].unique(), name=col)\n    rows = []\n    for smi in uniq:\n        mol = Chem.MolFromSmiles(smi)\n        rows.append({k: (np.nan if mol is None else fn(mol)) for k, fn in BB_DESCRIPTORS.items()})\n    tables[col] = pd.DataFrame(rows, index=uniq)\n    print(col, '->', len(uniq), 'unique blocks described')\n\nprint()\nprint('=== Mean RDKit descriptor value per building-block position ===')\nsummary = pd.DataFrame({c.replace('_smiles', ''): t.mean() for c, t in tables.items()})\nprint(summary.round(2))\n\nprint()\nprint('=== Correlation of each block descriptor with that block bind rate ===')\nfor col in bb_cols:\n    rate = df.groupby(col)['binds'].agg(n='size', bind_rate='mean')\n    joined = tables[col].join(rate)\n    joined = joined[joined['n'] >= 50]\n    corr = (joined.drop(columns=['n']).corr(numeric_only=True)['bind_rate']\n                  .drop('bind_rate').sort_values(key=abs, ascending=False))\n    print()\n    print(col, '(' + str(len(joined)) + ' blocks with n>=50)')\n    print(corr.round(3))\n\nfig, axes = plt.subplots(1, 3, figsize=(16, 4))\nfor ax, col in zip(axes, bb_cols):\n    rate = df.groupby(col)['binds'].agg(n='size', bind_rate='mean')\n    joined = tables[col].join(rate)\n    joined = joined[joined['n'] >= 50]\n    sns.scatterplot(data=joined, x='mol_wt', y='bind_rate', hue='arom_rings', palette='viridis', ax=ax)\n    ax.set_title(col.replace('_smiles', '') + ': block size vs bind rate')\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:10:26.450734Z","iopub.execute_input":"2026-08-13T06:10:26.451420Z","iopub.status.idle":"2026-08-13T06:10:29.251431Z","shell.execute_reply.started":"2026-08-13T06:10:26.451372Z","shell.execute_reply":"2026-08-13T06:10:29.250380Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- EDA 5 / feature engineering: Morgan (ECFP) fingerprints ----------\n# This is the step that was missing before: the preprocessing cell below expects `X`\n# and a `smiles_to_fp` helper (also reused for the test set).\nFP_BITS, FP_RADIUS = 2048, 2\nfp_gen = rdFingerprintGenerator.GetMorganGenerator(radius=FP_RADIUS, fpSize=FP_BITS)\n\ndef smiles_to_fp(smiles):\n    \"\"\"SMILES string -> uint8 numpy vector of Morgan fingerprint bits.\"\"\"\n    mol = Chem.MolFromSmiles(smiles)\n    if mol is None:\n        return np.zeros(FP_BITS, dtype=np.uint8)\n    return np.asarray(fp_gen.GetFingerprintAsNumPy(mol), dtype=np.uint8)\n\nstart = time.time()\nX = np.stack(df['molecule_smiles'].apply(smiles_to_fp).values)\nprint(f'Fingerprint matrix {X.shape} ({X.nbytes / 1e6:.0f} MB) built in {time.time() - start:.1f}s')\n\nbits_per_mol = X.sum(axis=1)\nbit_freq = X.mean(axis=0)\nprint(f'Bits set per molecule: mean {bits_per_mol.mean():.1f}, min {bits_per_mol.min()}, max {bits_per_mol.max()}')\nprint('Never-active bits:', int((bit_freq == 0).sum()), 'of', FP_BITS)\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 4))\nsns.histplot(bits_per_mol, bins=40, ax=axes[0])\naxes[0].set_title('Active fingerprint bits per molecule')\nsns.histplot(bit_freq[bit_freq > 0], bins=40, ax=axes[1])\naxes[1].set_title('How often each bit is switched on')\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:10:29.252783Z","iopub.execute_input":"2026-08-13T06:10:29.253177Z","iopub.status.idle":"2026-08-13T06:11:56.316514Z","shell.execute_reply.started":"2026-08-13T06:10:29.253138Z","shell.execute_reply":"2026-08-13T06:11:56.315410Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Data Preprocessing","metadata":{}},{"cell_type":"code","source":"# Encode protein_name as a numeric feature and combine with the fingerprint bits.\nprotein_encoder = LabelEncoder()\nprotein_encoded = protein_encoder.fit_transform(df['protein_name']).astype(np.uint8).reshape(-1, 1)\nprint('Protein classes:', dict(zip(protein_encoder.classes_, range(len(protein_encoder.classes_)))))\n\n# Keep everything uint8: hstack with an int64 column would blow the matrix up ~8x in RAM.\nX_full = np.hstack([X, protein_encoded])\ny_full = df['binds'].values\nprint(f'X_full {X_full.shape} dtype={X_full.dtype} ({X_full.nbytes / 1e6:.0f} MB)')\n\n# Split into train/validation. We stratify on binds to keep the class ratio similar\n# in both splits (this sample is already artificially balanced 50/50).\nX_train, X_val, y_train, y_val, protein_train, protein_val = train_test_split(\n    X_full, y_full, df['protein_name'].values,\n    test_size=0.2, random_state=42, stratify=y_full\n)\n\nprint('Train shape:', X_train.shape, 'Val shape:', X_val.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:11:56.317783Z","iopub.execute_input":"2026-08-13T06:11:56.318205Z","iopub.status.idle":"2026-08-13T06:11:56.734964Z","shell.execute_reply.started":"2026-08-13T06:11:56.318169Z","shell.execute_reply":"2026-08-13T06:11:56.734177Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Model Training","metadata":{}},{"cell_type":"code","source":"# A Random Forest baseline on Morgan fingerprint features + encoded protein target.\n# n_jobs=-1 uses all available CPU cores to speed up training.\nmodel = RandomForestClassifier(\n    n_estimators=200,\n    max_depth=None,\n    min_samples_leaf=2,\n    n_jobs=-1,\n    random_state=42,\n    class_weight='balanced'\n)\n\nstart = time.time()\nmodel.fit(X_train, y_train)\nprint(f'Trained in {time.time() - start:.1f} seconds')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:11:56.736181Z","iopub.execute_input":"2026-08-13T06:11:56.736514Z","iopub.status.idle":"2026-08-13T06:14:32.283732Z","shell.execute_reply.started":"2026-08-13T06:11:56.736479Z","shell.execute_reply":"2026-08-13T06:14:32.282729Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Evaluation","metadata":{}},{"cell_type":"code","source":"val_probs = model.predict_proba(X_val)[:, 1]\n\n# The competition metric is average precision computed PER PROTEIN, then averaged.\nresults = []\nfor protein in np.unique(protein_val):\n    mask = protein_val == protein\n    ap = average_precision_score(y_val[mask], val_probs[mask])\n    results.append({'protein_name': protein, 'average_precision': ap})\n\nresults_df = pd.DataFrame(results)\nprint(results_df)\nprint('\\nMean average precision across proteins:', results_df['average_precision'].mean())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:14:32.285019Z","iopub.execute_input":"2026-08-13T06:14:32.285409Z","iopub.status.idle":"2026-08-13T06:14:34.290074Z","shell.execute_reply.started":"2026-08-13T06:14:32.285368Z","shell.execute_reply":"2026-08-13T06:14:34.289201Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- Extra diagnostics on the validation split ----------\nprint('Micro-averaged average precision:', average_precision_score(y_val, val_probs))\nprint('ROC AUC                         :', roc_auc_score(y_val, val_probs))\n\nfig, axes = plt.subplots(1, 3, figsize=(16, 4))\n\nsns.barplot(x='protein_name', y='average_precision', data=results_df, ax=axes[0])\naxes[0].set_ylim(0, 1)\naxes[0].set_title('Average precision per protein')\n\nsns.histplot(x=val_probs, hue=y_val, bins=50, stat='density',\n             common_norm=False, element='step', ax=axes[1])\naxes[1].set_title('Predicted P(binds) by true class')\naxes[1].set_xlabel('predicted probability')\n\nimp = pd.Series(model.feature_importances_)\ntop = imp.sort_values(ascending=False).head(20)\nlabels = ['protein_name' if i == X_full.shape[1] - 1 else f'bit {i}' for i in top.index]\nsns.barplot(x=top.values, y=labels, ax=axes[2])\naxes[2].set_title('Top 20 most important features')\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:14:34.291263Z","iopub.execute_input":"2026-08-13T06:14:34.291723Z","iopub.status.idle":"2026-08-13T06:14:34.975408Z","shell.execute_reply.started":"2026-08-13T06:14:34.291660Z","shell.execute_reply":"2026-08-13T06:14:34.974310Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Testing on Test Data / Submission","metadata":{}},{"cell_type":"code","source":"# Load the official test set directly from parquet (no full copy needed beyond this query)\ncon = duckdb.connect()\ntest_df = con.query(\n    f\"SELECT id, molecule_smiles, protein_name FROM parquet_scan('{test_path}')\"\n).df()\ncon.close()\nprint('Test rows:', len(test_df))\n\n# protein_encoder was fit on train; the same 3 proteins appear in train and test,\n# so transform() is safe here.\ntest_protein_encoded = protein_encoder.transform(test_df['protein_name']).astype(np.uint8).reshape(-1, 1)\n\n# Featurize and predict in chunks: building one dense matrix for the whole test set\n# would need tens of GB of RAM, which is what used to make this step blow up.\nCHUNK = 100_000\ntest_probs = np.empty(len(test_df), dtype=np.float32)\nstart = time.time()\nfor lo in range(0, len(test_df), CHUNK):\n    hi = min(lo + CHUNK, len(test_df))\n    fps = np.stack(test_df['molecule_smiles'].iloc[lo:hi].apply(smiles_to_fp).values)\n    chunk = np.hstack([fps, test_protein_encoded[lo:hi]])\n    test_probs[lo:hi] = model.predict_proba(chunk)[:, 1]\n    print(f'  {hi:,}/{len(test_df):,} rows scored ({time.time() - start:.0f}s elapsed)')\n\nprint(f'Predicted {len(test_probs):,} probabilities in {time.time() - start:.0f}s')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:14:34.978358Z","iopub.execute_input":"2026-08-13T06:14:34.978698Z","iopub.status.idle":"2026-08-13T06:28:23.967073Z","shell.execute_reply.started":"2026-08-13T06:14:34.978669Z","shell.execute_reply":"2026-08-13T06:28:23.965633Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- Build the submission file ----------\nsubmission = pd.DataFrame({'id': test_df['id'].values, 'binds': test_probs})\nsubmission.to_csv('submission.csv', index=False)\n\nprint('submission.csv written:', submission.shape)\nprint(submission['binds'].describe())\nsubmission.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:28:23.968805Z","iopub.execute_input":"2026-08-13T06:28:23.969329Z","iopub.status.idle":"2026-08-13T06:28:26.948088Z","shell.execute_reply.started":"2026-08-13T06:28:23.969271Z","shell.execute_reply":"2026-08-13T06:28:26.947018Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}