{"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":"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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-08-13T06:45:53.755501Z","iopub.execute_input":"2026-08-13T06:45:53.755746Z","iopub.status.idle":"2026-08-13T06:45:53.762207Z","shell.execute_reply.started":"2026-08-13T06:45:53.755726Z","shell.execute_reply":"2026-08-13T06:45:53.761586Z"}},"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# New strategy: use ALL 1.59M binders + 4x random non-binders for training.\n# Validation = natural ~0.5% prevalence sample, split by hash(molecule_smiles)\n# so all 3 protein rows of a molecule stay on the same side (no molecule overlap).\nCOLS = 'id, molecule_smiles, protein_name, binds'\nTRAIN_SIDE = \"hash(molecule_smiles) % 20 <> 0\"\nVAL_SIDE = \"hash(molecule_smiles) % 20 = 0\"\n\ncon = duckdb.connect()\nstart = time.time()\ntrain_df = con.query(f\"\"\"\nSELECT {COLS} FROM parquet_scan('{train_path}') WHERE binds = 1 AND {TRAIN_SIDE}\nUNION ALL\nSELECT * FROM (SELECT {COLS} FROM parquet_scan('{train_path}') WHERE binds = 0 AND {TRAIN_SIDE}) USING SAMPLE reservoir(6000000 ROWS) REPEATABLE (42)\n\"\"\"\"\").df()\nval_df = con.query(f\"\"\"\nSELECT * FROM (SELECT {COLS} FROM parquet_scan('{train_path}') WHERE {VAL_SIDE}) USING SAMPLE reservoir(1000000 ROWS) REPEATABLE (43)\n\"\"\"\"\").df()\ncon.close()\n\n# Shuffle so positives and negatives are interleaved\ntrain_df = train_df.sample(frac=1.0, random_state=42).reset_index(drop=True)\nprint(f'Loaded in {time.time()-start:.0f}s')\nprint('Train:', train_df.shape, 'positives:', int(train_df.binds.sum()), f'({train_df.binds.mean():.3%})')\nprint('Val:  ', val_df.shape, 'positives:', int(val_df.binds.sum()), f'({val_df.binds.mean():.3%})')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T07:15:45.671629Z","iopub.execute_input":"2026-08-13T07:15:45.671901Z","iopub.status.idle":"2026-08-13T07:17:43.695055Z","shell.execute_reply.started":"2026-08-13T07:15:45.671882Z","shell.execute_reply":"2026-08-13T07:17:43.694320Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Profile the FULL BELKA train set with DuckDB - scans parquet directly, no RAM blowup\nimport duckdb, os\ncon = duckdb.connect()\nprint('train.parquet size:', round(os.path.getsize(train_path)/1e9, 2), 'GB')\nprint(con.query(f\"SELECT COUNT(*) AS total_rows FROM parquet_scan('{train_path}')\").df())\nprint(con.query(f\"SELECT binds, COUNT(*) AS n FROM parquet_scan('{train_path}') GROUP BY binds ORDER BY binds\").df())\nprint(con.query(f\"SELECT protein_name, binds, COUNT(*) AS n FROM parquet_scan('{train_path}') GROUP BY protein_name, binds ORDER BY protein_name, binds\").df())\nprint(con.query(f\"SELECT COUNT(*) AS test_rows FROM parquet_scan('{test_path}')\").df())\ncon.close()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T07:06:40.092726Z","iopub.execute_input":"2026-08-13T07:06:40.093577Z","iopub.status.idle":"2026-08-13T07:06:43.560996Z","shell.execute_reply.started":"2026-08-13T07:06:40.093552Z","shell.execute_reply":"2026-08-13T07:06:43.560338Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print('Train sample (all binders kept, 1:4 neg ratio):')\ntrain_df.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T07:03:46.948645Z","iopub.execute_input":"2026-08-13T07:03:46.949792Z","iopub.status.idle":"2026-08-13T07:03:46.980339Z","shell.execute_reply.started":"2026-08-13T07:03:46.949743Z","shell.execute_reply":"2026-08-13T07:03:46.979705Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Sparse Morgan fingerprints (2048 bits, radius 2) + 3 one-hot protein columns.\n# Dense uint8 would be ~2 KB/row -> 15+ GB for 7.5M rows; sparse (~60 set bits\n# per molecule) is ~50x smaller and LightGBM consumes CSR matrices directly.\nfrom multiprocessing import Pool\nfrom scipy import sparse\n\nFP_SIZE = 2048\nPROTEIN_MAP = {'BRD4': 0, 'HSA': 1, 'sEH': 2}\nfp_gen = None\n\ndef _init_worker():\n    global fp_gen\n    fp_gen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=FP_SIZE)\n\ndef _onbits_chunk(smiles_chunk):\n    global fp_gen\n    if fp_gen is None:\n        _init_worker()\n    bits = []\n    lens = np.empty(len(smiles_chunk), dtype=np.int32)\n    for i, s in enumerate(smiles_chunk):\n        mol = Chem.MolFromSmiles(s)\n        ob = list(fp_gen.GetFingerprint(mol).GetOnBits()) if mol is not None else []\n        bits.extend(ob)\n        lens[i] = len(ob)\n    return np.array(bits, dtype=np.int32), lens\n\ndef featurize(smiles, proteins, n_workers=4, chunk_size=25000, tag=''):\n    chunks = [smiles[i:i+chunk_size] for i in range(0, len(smiles), chunk_size)]\n    all_bits, all_lens, done, t0 = [], [], 0, time.time()\n    with Pool(n_workers, initializer=_init_worker) as pool:\n        for b, l in pool.imap(_onbits_chunk, chunks):\n            all_bits.append(b)\n            all_lens.append(l)\n            done += len(l)\n            if done % 1000000 < chunk_size:\n                print(f'  {tag} {done:,}/{len(smiles):,} ({time.time()-t0:.0f}s)')\n    indices = np.concatenate(all_bits)\n    lens = np.concatenate(all_lens)\n    indptr = np.zeros(len(smiles) + 1, dtype=np.int64)\n    np.cumsum(lens, out=indptr[1:])\n    fps = sparse.csr_matrix((np.ones(len(indices), dtype=np.uint8), indices, indptr), shape=(len(smiles), FP_SIZE))\n    prot = np.array([PROTEIN_MAP[p] for p in proteins], dtype=np.int64)\n    onehot = sparse.csr_matrix((np.ones(len(prot), dtype=np.uint8), prot, np.arange(len(prot) + 1, dtype=np.int64)), shape=(len(prot), 3))\n    return sparse.hstack([fps, onehot], format='csr')\n\nX_train = featurize(train_df['molecule_smiles'].tolist(), train_df['protein_name'].tolist(), tag='train')\ny_train = train_df['binds'].values\nX_val = featurize(val_df['molecule_smiles'].tolist(), val_df['protein_name'].tolist(), tag='val')\ny_val = val_df['binds'].values\nprint('X_train:', X_train.shape, f'({X_train.data.nbytes/1e9:.2f} GB data)', ' X_val:', X_val.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T07:20:37.364283Z","iopub.execute_input":"2026-08-13T07:20:37.364968Z","iopub.status.idle":"2026-08-13T07:38:03.873323Z","shell.execute_reply.started":"2026-08-13T07:20:37.364942Z","shell.execute_reply":"2026-08-13T07:38:03.872236Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# LightGBM instead of RandomForest: histogram-based, handles 7.5M x 2051 sparse\n# comfortably, and early-stops on the competition metric (average precision).\nimport lightgbm as lgb\n\n# LightGBM requires float32/float64 sparse data (uint8 raises TypeError)\nif X_train.dtype != np.float32:\n    X_train = X_train.astype(np.float32)\nif X_val.dtype != np.float32:\n    X_val = X_val.astype(np.float32)\n\nparams = dict(\n    objective='binary',\n    metric='average_precision',\n    learning_rate=0.1,\n    num_leaves=127,\n    min_data_in_leaf=200,\n    feature_fraction=0.7,\n    bagging_fraction=0.8,\n    bagging_freq=1,\n    num_threads=4,\n    seed=42,\n    verbosity=-1,\n)\n\ndtrain = lgb.Dataset(X_train, label=y_train)\ndval = lgb.Dataset(X_val, label=y_val, reference=dtrain)\n\nstart = time.time()\nmodel = lgb.train(\n    params, dtrain, num_boost_round=1500,\n    valid_sets=[dval], valid_names=['val'],\n    callbacks=[lgb.early_stopping(50), lgb.log_evaluation(25)],\n)\nprint(f'Best iteration: {model.best_iteration}, trained in {(time.time()-start)/60:.1f} min')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T07:43:15.879814Z","iopub.execute_input":"2026-08-13T07:43:15.880820Z","iopub.status.idle":"2026-08-13T08:02:37.356554Z","shell.execute_reply.started":"2026-08-13T07:43:15.880782Z","shell.execute_reply":"2026-08-13T08:02:37.355626Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"val_probs = model.predict(X_val, num_iteration=model.best_iteration)\n\nprint(f'Val prevalence (baseline AP for a random model): {y_val.mean():.4f}')\nprint(f'Overall val AP:  {average_precision_score(y_val, val_probs):.4f}')\nprint(f'Overall val AUC: {roc_auc_score(y_val, val_probs):.4f}')\n\n# The competition metric is average precision computed PER PROTEIN, then averaged.\nresults = []\nfor protein in ['BRD4', 'HSA', 'sEH']:\n    mask = (val_df['protein_name'] == protein).values\n    ap = average_precision_score(y_val[mask], val_probs[mask])\n    results.append({'protein_name': protein, 'n': int(mask.sum()), 'positives': int(y_val[mask].sum()), 'average_precision': ap})\n\nresults_df = pd.DataFrame(results)\nprint(results_df)\nprint('\\nMean average precision across proteins:', round(results_df['average_precision'].mean(), 4))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T08:06:11.386058Z","iopub.execute_input":"2026-08-13T08:06:11.386485Z","iopub.status.idle":"2026-08-13T08:07:17.056827Z","shell.execute_reply.started":"2026-08-13T08:06:11.386457Z","shell.execute_reply":"2026-08-13T08:07:17.055196Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ---------- Extra diagnostics on the validation split ----------\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_yscale('log')\naxes[1].set_title('Predicted P(binds) by true class')\naxes[1].set_xlabel('predicted probability')\n\nimp = pd.Series(model.feature_importance(importance_type='gain'))\ntop = imp.sort_values(ascending=False).head(20)\nlabels = [['BRD4', 'HSA', 'sEH'][i - FP_SIZE] if i >= FP_SIZE 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 (gain)')\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T08:08:12.683430Z","iopub.execute_input":"2026-08-13T08:08:12.683922Z","iopub.status.idle":"2026-08-13T08:08:14.028788Z","shell.execute_reply.started":"2026-08-13T08:08:12.683869Z","shell.execute_reply":"2026-08-13T08:08:14.027799Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Score the official test set in chunks and write submission.csv\ncon = duckdb.connect()\ntest_df = con.query(f\"SELECT id, molecule_smiles, protein_name FROM parquet_scan('{test_path}')\").df()\ncon.close()\nprint('Test rows:', len(test_df))\n\nCHUNK = 200_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    Xt = featurize(test_df['molecule_smiles'].iloc[lo:hi].tolist(), test_df['protein_name'].iloc[lo:hi].tolist()).astype(np.float32)\n    test_probs[lo:hi] = model.predict(Xt, num_iteration=model.best_iteration)\n    print(f'  {hi:,}/{len(test_df):,} rows scored ({time.time()-start:.0f}s elapsed)')\n\nsub = pd.DataFrame({'id': test_df['id'], 'binds': test_probs})\nsub.to_csv('submission.csv', index=False)\nprint('Wrote submission.csv')\nsub.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T08:08:43.918475Z","iopub.execute_input":"2026-08-13T08:08:43.918925Z","iopub.status.idle":"2026-08-13T08:14:36.043408Z","shell.execute_reply.started":"2026-08-13T08:08:43.918893Z","shell.execute_reply":"2026-08-13T08:14:36.042091Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T07:11:15.238859Z","iopub.execute_input":"2026-08-13T07:11:15.239282Z","iopub.status.idle":"2026-08-13T07:11:15.248788Z","shell.execute_reply.started":"2026-08-13T07:11:15.239260Z","shell.execute_reply":"2026-08-13T07:11:15.248085Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Apples-to-apples check: score the SAME LightGBM on a 50/50 balanced\n# subset of val (positives + equal number of random negatives), which is\n# how the old RandomForest was evaluated.\nrng = np.random.default_rng(42)\npos_idx = np.where(y_val == 1)[0]\nneg_idx = rng.choice(np.where(y_val == 0)[0], size=len(pos_idx), replace=False)\nbal = np.concatenate([pos_idx, neg_idx])\n\nprint(f'Balanced 50/50 val subset: {len(bal):,} rows')\nprint(f'Micro-averaged AP on balanced val: {average_precision_score(y_val[bal], val_probs[bal]):.4f}')\nprint(f'ROC AUC on balanced val          : {roc_auc_score(y_val[bal], val_probs[bal]):.4f}')\nprot_bal = val_df['protein_name'].values[bal]\naps = []\nfor protein in ['BRD4', 'HSA', 'sEH']:\n    m = prot_bal == protein\n    ap = average_precision_score(y_val[bal][m], val_probs[bal][m])\n    aps.append(ap)\n    print(f'{protein}: AP = {ap:.4f}')\nprint(f'Mean AP across proteins (balanced val): {np.mean(aps):.4f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-13T08:22:25.147802Z","iopub.execute_input":"2026-08-13T08:22:25.148175Z","iopub.status.idle":"2026-08-13T08:22:25.291728Z","shell.execute_reply.started":"2026-08-13T08:22:25.148145Z","shell.execute_reply":"2026-08-13T08:22:25.290660Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}