{"metadata":{"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"}],"dockerImageVersionId":30684,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Baseline model - Simple NN","metadata":{}},{"cell_type":"markdown","source":"This notebook presents a baseline Neural Network model for predicting the binding affinity of small molecules to specific protein targets, as part of the BELKA competition on Kaggle.","metadata":{}},{"cell_type":"code","source":"%pip install --upgrade pip\n\nimport subprocess\nimport sys\n\ndef import_or_install(package):\n    global_vars = globals()\n    try:\n        global_vars[package] = __import__(package)\n        print(f'Package {package} has been imported')\n    except ImportError:\n        subprocess.check_call([sys.executable, '-m', 'pip', 'install', package])\n        global_vars[package] = __import__(package)\n        print(f'Package {package} has been installed and imported')","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:05:54.512348Z","iopub.execute_input":"2024-05-09T09:05:54.512675Z","iopub.status.idle":"2024-05-09T09:06:08.866860Z","shell.execute_reply.started":"2024-05-09T09:05:54.512643Z","shell.execute_reply":"2024-05-09T09:06:08.865484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pathlib\nimport seaborn as sns\nfrom sklearn.preprocessing import OneHotEncoder \nfrom sklearn.feature_extraction.text import HashingVectorizer\nfrom tensorflow import keras\nimport tensorflow as tf\n\n# import_or_install(\"rdkit\")\nimport_or_install(\"duckdb\")","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:06:08.872914Z","iopub.execute_input":"2024-05-09T09:06:08.873271Z","iopub.status.idle":"2024-05-09T09:06:29.921795Z","shell.execute_reply.started":"2024-05-09T09:06:08.873229Z","shell.execute_reply":"2024-05-09T09:06:29.920683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The data are too large to handle directly. Efficient OLAP database is required","metadata":{}},{"cell_type":"code","source":"%%time\ninput_path = pathlib.Path(\"/kaggle/input/leash-BELKA\")\ntrain_file = str(input_path / \"train.parquet\")\ntest_file = str(input_path / \"test.parquet\")","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:06:29.923170Z","iopub.execute_input":"2024-05-09T09:06:29.923776Z","iopub.status.idle":"2024-05-09T09:06:29.930548Z","shell.execute_reply.started":"2024-05-09T09:06:29.923745Z","shell.execute_reply":"2024-05-09T09:06:29.929450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%time\nconn = duckdb.connect(database=':memory:')\nprint(f\"Connected to {conn}\")\nconn.execute(\n    f\"CREATE OR REPLACE VIEW filtered_train AS SELECT * \"\n    f\"FROM read_parquet('{train_file}')\"\n    \"WHERE binds IS NOT NULL AND molecule_smiles IS NOT NULL\"\n)\nprint(\"Data has been loaded and filtered in DuckDB\")\n\ncolumns_df = conn.execute(\n    \"SELECT column_name FROM information_schema.columns \"\n    \"WHERE table_name = 'filtered_train'\"\n).fetchdf()\nprint(columns_df)","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:06:29.934515Z","iopub.execute_input":"2024-05-09T09:06:29.935082Z","iopub.status.idle":"2024-05-09T09:06:30.033471Z","shell.execute_reply.started":"2024-05-09T09:06:29.935041Z","shell.execute_reply":"2024-05-09T09:06:30.032307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_types = conn.execute(\n    \"DESCRIBE filtered_train\"\n).fetchdf()\nprint(data_types)\n","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:06:30.034802Z","iopub.execute_input":"2024-05-09T09:06:30.035164Z","iopub.status.idle":"2024-05-09T09:06:30.046266Z","shell.execute_reply.started":"2024-05-09T09:06:30.035134Z","shell.execute_reply":"2024-05-09T09:06:30.045050Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get the list of column names\ncolumns_df = conn.execute(\n    \"SELECT column_name FROM information_schema.columns \"\n    \"WHERE table_name = 'filtered_train'\"\n).fetchdf()\ncolumn_names = columns_df['column_name'].tolist()\n\n# Prepare SQL to count non-nulls for each column dynamically\nsql_query = \"\"\"\nSELECT COUNT(*) AS total_count,\n\"\"\"\n# Adding each column dynamically to the SQL query\nsql_query += \",\\n\".join([f\"COUNT({col}) AS count_{col}\" for col in column_names])\n\nsql_query += f\"\\nFROM filtered_train\"\n\n# Execute the query\ncount_df = conn.execute(sql_query).fetchdf()\n\n# Calculate missing values for each column dynamically\nfor col in column_names:\n    count_df[f'missing_{col}'] = count_df['total_count'] - count_df[f'count_{col}']\n\n# Display the results\nmissing_columns = [f'missing_{col}' for col in column_names]\nprint(count_df[['total_count'] + missing_columns])\n","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:06:30.048099Z","iopub.execute_input":"2024-05-09T09:06:30.048425Z","iopub.status.idle":"2024-05-09T09:07:02.463105Z","shell.execute_reply.started":"2024-05-09T09:06:30.048396Z","shell.execute_reply":"2024-05-09T09:07:02.461783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"summary_stats = conn.execute(\n    \"SELECT AVG(binds) AS mean, STDDEV_SAMP(binds) AS stddev, \"\n    \"MIN(binds) AS min, MAX(binds) AS max \"\n    \"FROM filtered_train\"\n).fetchdf()\nprint(summary_stats)\n\nunique_values = conn.execute(\n    \"SELECT protein_name, COUNT(*) AS frequency \"\n    \"FROM filtered_train \"\n    \"GROUP BY protein_name \"\n    \"ORDER BY frequency DESC\"\n).fetchdf()\nprint(unique_values)","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:07:02.464886Z","iopub.execute_input":"2024-05-09T09:07:02.465408Z","iopub.status.idle":"2024-05-09T09:07:48.611849Z","shell.execute_reply.started":"2024-05-09T09:07:02.465364Z","shell.execute_reply":"2024-05-09T09:07:48.610672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Continuous Data Plotting histograms or density plots for continuous data like binds:\nbucketed_data = conn.execute(\"\"\"\n  SELECT\n    floor(binds / 100.0) * 100 AS binds_bucket,\n    COUNT(*) AS count\n  FROM filtered_train\n  GROUP BY binds_bucket\n  ORDER BY binds_bucket\n\"\"\").fetchdf()\n\n# Plotting\nplt.figure(figsize=(10, 6))\nsns.barplot(x='binds_bucket', y='count', data=bucketed_data)\nplt.title('Distribution of Binds')\nplt.xlabel('Binds Bucket')\nplt.ylabel('Frequency')\nplt.xticks(rotation=45)  # Rotate x labels for better readability\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:07:48.613990Z","iopub.execute_input":"2024-05-09T09:07:48.614476Z","iopub.status.idle":"2024-05-09T09:08:13.055729Z","shell.execute_reply.started":"2024-05-09T09:07:48.614432Z","shell.execute_reply":"2024-05-09T09:08:13.054706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# #  Approximate the percentiles with a sample of the data\n# percentiles = conn.execute(\n#     \"\"\"\n#     WITH percentile_values AS (\n#         SELECT \n#             percentile_disc(0.01) WITHIN GROUP (ORDER BY binds) AS p1,\n#             percentile_disc(0.25) WITHIN GROUP (ORDER BY binds) AS p25,\n#             percentile_disc(0.50) WITHIN GROUP (ORDER BY binds) AS p50,\n#             percentile_disc(0.75) WITHIN GROUP (ORDER BY binds) AS p75,\n#             percentile_disc(0.99) WITHIN GROUP (ORDER BY binds) AS p99\n#         FROM filtered_train\n#     ),\n#     percentile_counts AS (\n#         SELECT\n#             (CASE \n#                 WHEN binds <= p1 THEN '1st Percentile'\n#                 WHEN binds <= p25 THEN '1-25th Percentile'\n#                 WHEN binds <= p50 THEN '25-50th Percentile'\n#                 WHEN binds <= p75 THEN '50-75th Percentile'\n#                 WHEN binds <= p99 THEN '75-99th Percentile'\n#                 ELSE 'Above 99th Percentile'\n#              END) AS percentile_bucket,\n#              COUNT(*) AS count\n#         FROM filtered_train, percentile_values\n#         GROUP BY percentile_bucket\n#     )\n#     SELECT * FROM percentile_counts ORDER BY percentile_bucket\n#     \"\"\"\n# ).fetchdf()\n\n# # Plotting with adjusted buckets\n# plt.figure(figsize=(12, 8))\n# sns.barplot(x='percentile_bucket', y='count', data=percentiles)\n# plt.title('Distribution of Binds Across Percentiles')\n# plt.xlabel('Percentile Bucket')\n# plt.ylabel('Frequency')\n# plt.xticks(rotation=45)\n# plt.show()\n\n","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:08:13.056774Z","iopub.execute_input":"2024-05-09T09:08:13.057130Z","iopub.status.idle":"2024-05-09T09:08:13.064203Z","shell.execute_reply.started":"2024-05-09T09:08:13.057100Z","shell.execute_reply":"2024-05-09T09:08:13.063142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Assuming `unique_values` contains the protein names and their frequencies\n# plt.figure(figsize=(12, 8))\n# sns.barplot(x='frequency', y='protein_name', data=unique_values)\n# plt.title('Frequency of Protein Names')\n# plt.xlabel('Frequency')\n# plt.ylabel('Protein Names')\n# plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:08:13.065431Z","iopub.execute_input":"2024-05-09T09:08:13.065747Z","iopub.status.idle":"2024-05-09T09:08:13.076098Z","shell.execute_reply.started":"2024-05-09T09:08:13.065719Z","shell.execute_reply":"2024-05-09T09:08:13.075047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Box Plots to Identify Outliers For numerical columns to identify outliers:\n# plt.figure(figsize=(10, 6))\n# sns.boxplot(data=binds_data['binds'])\n# plt.title('Box Plot of Binds')\n# plt.xlabel('Binds')\n# plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:08:13.077592Z","iopub.execute_input":"2024-05-09T09:08:13.078253Z","iopub.status.idle":"2024-05-09T09:08:13.083678Z","shell.execute_reply.started":"2024-05-09T09:08:13.078214Z","shell.execute_reply":"2024-05-09T09:08:13.082582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%time\nn_features = 100000\nvectorizer = HashingVectorizer(n_features=n_features, norm=None, alternate_sign=False)","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:08:13.084929Z","iopub.execute_input":"2024-05-09T09:08:13.085267Z","iopub.status.idle":"2024-05-09T09:08:13.093705Z","shell.execute_reply.started":"2024-05-09T09:08:13.085242Z","shell.execute_reply":"2024-05-09T09:08:13.092690Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train simple model in chunks\n\nFor the reason mentioned, the model can be trained in chunks","metadata":{}},{"cell_type":"code","source":"%%time\nmodel = keras.models.Sequential([\n    keras.layers.Dense(\n        128,\n        activation='relu',\n        input_shape=(n_features,)\n    ),\n    keras.layers.Dropout(0.5),\n    keras.layers.Dense(64, activation='relu'),\n    keras.layers.Dense(1, activation='sigmoid')\n])\n\nmodel.compile(\n    optimizer=keras.optimizers.Adam(),\n    loss='binary_crossentropy',\n    metrics=['accuracy']\n)\n","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:08:13.099245Z","iopub.execute_input":"2024-05-09T09:08:13.100035Z","iopub.status.idle":"2024-05-09T09:08:13.976135Z","shell.execute_reply.started":"2024-05-09T09:08:13.099974Z","shell.execute_reply":"2024-05-09T09:08:13.975058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Training the model with data fetched in chunks","metadata":{}},{"cell_type":"code","source":"%%time\nchunksize = 100000\nsteps_per_epoch = 10  # This should be adjusted based on the total number of rows / chunksize\n\ndef fetch_dataset(query, offset=0, chunksize=100000):\n    while True:\n        chunk_query = f\"{query} LIMIT {chunksize} OFFSET {offset}\"\n        chunk_df = conn.execute(chunk_query).fetchdf()\n        if chunk_df.empty:\n            break\n        X = vectorizer.transform(chunk_df['molecule_smiles'].values)\n        X_dense = X.todense()  # Convert sparse matrix to dense\n        y = chunk_df['binds'].values\n        yield (X_dense, y)\n        offset += chunksize\n\ndef tf_data_generator(query, chunksize):\n    return tf.data.Dataset.from_generator(\n        lambda: fetch_dataset(query, chunksize=chunksize),\n        output_signature=(\n            tf.TensorSpec(shape=(None, n_features), dtype=tf.float32),\n            tf.TensorSpec(shape=(None,), dtype=tf.float32)\n        )\n    ).prefetch(tf.data.AUTOTUNE)\n\n# Training loop\nquery = \"SELECT binds, molecule_smiles FROM filtered_train\"\ntotal_data_processed = 0\nall_histories = []\n\nwhile True:\n    train_dataset = tf_data_generator(query, chunksize)\n    history = model.fit(train_dataset, epochs=5, steps_per_epoch=steps_per_epoch)\n    all_histories.append(history)\n    \n    # Check if we need to stop\n    processed_in_this_epoch = steps_per_epoch * chunksize\n    total_data_processed += processed_in_this_epoch\n    if processed_in_this_epoch < chunksize:\n        print(\"All data processed.\")\n        break","metadata":{"execution":{"iopub.status.busy":"2024-05-09T09:08:13.977606Z","iopub.execute_input":"2024-05-09T09:08:13.978055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# def query_data(query, chunksize=100000):\n#     result = conn.execute(query)\n#     while True:\n#         chunk_df = result.fetch_df_chunk(chunksize)\n#         if chunk_df is None:\n#             break\n#         yield chunk_df\n\n# import gc\n# chunksize = 30000\n# all_histories = [] \n\n# query = \"SELECT binds, molecule_smiles FROM filtered_train\"\n# print(\"Starting model training\")\n# for i, chunk_df in enumerate(query_data(query, chunksize)):\n#     X_train = vectorizer.transform(chunk_df['molecule_smiles'].values)\n#     y_train = chunk_df['binds'].to_numpy()\n\n#     history = model.fit(X_train, y_train, epochs=5, batch_size=32, verbose=1)\n#     all_histories.append(history)\n#     print(f\"Completed chunk {i+1}\")\n\n#     # Clear memory\n#     del X_train, y_train, chunk_df\n#     gc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Initialize lists to store combined loss and accuracy\ncombined_loss = []\ncombined_accuracy = []\n\n# Loop through all collected histories and append their values\nfor history in all_histories:\n    combined_loss.extend(history.history['loss'])\n    combined_accuracy.extend(history.history['accuracy'])\n\n# Plotting the results\nplt.figure(figsize=(14, 5))\n\n# Plot training loss\nplt.subplot(1, 2, 1)\nplt.plot(combined_loss, label='Training Loss')\nplt.title('Training Loss Over Chunks')\nplt.xlabel('Training Iterations')\nplt.ylabel('Loss')\nplt.legend()\n\n# Plot training accuracy\nplt.subplot(1, 2, 2)\nplt.plot(combined_accuracy, label='Training Accuracy')\nplt.title('Training Accuracy Over Chunks')\nplt.xlabel('Training Iterations')\nplt.ylabel('Accuracy')\nplt.legend()\n\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Build prediction","metadata":{}},{"cell_type":"code","source":"%%time\nconn.register('test', test_file)\npredictions = []\nids = []\nfor chunk_df in query_data(\"SELECT id, molecule_smiles FROM test\", 100000):\n    X_test = encoder.transform(chunk_df[['molecule_smiles']])\n    chunk_predictions = model.predict(X_test).flatten()\n    predictions.extend(chunk_predictions)\n    ids.extend(chunk_df['id'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsubmission_df = pd.DataFrame({\n    'id': ids,\n    'binds': predictions\n})\nsubmission_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_df.to_csv('sample_submission.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}