{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Note\n\nThis is Fabien Crom's NB except that I increase the tSVD from 16 to 32 components. CV = 0.664 (from 0.662) and overall LB = 0.848 (fom 0.847).","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"# Multiome Quickstart With Sparse Matrices\n\nThis notebook is mostly for demonstrating the utility of sparse matrices in this competition. (Especially for the Multiome dataset).\n\nAs the Multiome dataset is  very sparse (about 98% of cells are zeros), it benefits greatly from being encoded as sparse matrices. \n\nThis notebook is largely based on [this notebook](https://www.kaggle.com/code/ambrosm/msci-multiome-quickstart) by AmbrosM. It is a nice first attempt at handling Multiome data, and I thought it would informative for kagglers to be able to contrast directly the performances of sparse vs dense representations. \n\nMostly, the differences with AmbrosM's notebooks are:\n- We use a representation of the data in sparse CSR format, which let us load all of the training data in memory (using less than 8GB memory instead of the >90GB it would take to represent the data in a dense format)\n- We perform PCA (actually, TruncatedSVD) on the totality of the training data (while AmbrosM's notebook had to work with a subset of 6000 rows and 4000 columns). \n- We keep 16 components (vs 4 in AmbrosM's notebook)\n- We apply Ridge regression on 50000 rows (vs 6000 in AmbrosM's notebook)\n- Despite using much more data, this notebook should run in a bit more than 10 minutes (vs >1h for AmbrosM's notebook)\n\nThe competition data is pre-encoded as sparse matrices in [this dataset](https://www.kaggle.com/datasets/fabiencrom/multimodal-single-cell-as-sparse-matrix) generated by [this notebook](https://www.kaggle.com/code/fabiencrom/multimodal-single-cell-creating-sparse-data/).\n\nSince we will only generate the multiome predictions in this notebook, I am taking the CITEseq predictions from [this notebook](https://www.kaggle.com/code/vuonglam/lgbm-baseline-optuna-drop-constant-cite-task) by VuongLam, which is the public notebook with the best score at the time I am publishing.","metadata":{}},{"cell_type":"code","source":"import os, gc, pickle\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom colorama import Fore, Back, Style\nfrom matplotlib.ticker import MaxNLocator\n\nfrom sklearn.base import BaseEstimator, TransformerMixin\nfrom sklearn.model_selection import KFold\nfrom sklearn.preprocessing import StandardScaler, scale\nfrom sklearn.decomposition import PCA, TruncatedSVD\nfrom sklearn.dummy import DummyRegressor\nfrom sklearn.pipeline import make_pipeline, Pipeline\nfrom sklearn.linear_model import Ridge, LinearRegression, Lasso\nfrom sklearn.metrics import mean_squared_error\n\nimport scipy\nimport scipy.sparse","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-09T15:45:24.797608Z","iopub.execute_input":"2022-09-09T15:45:24.798532Z","iopub.status.idle":"2022-09-09T15:45:26.240397Z","shell.execute_reply.started":"2022-09-09T15:45:24.798393Z","shell.execute_reply":"2022-09-09T15:45:26.239396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessing and cross-validation\n\nWe first load all of the training input data for Multiome. It should take less than a minute.","metadata":{}},{"cell_type":"code","source":"%%time\ntrain_inputs = scipy.sparse.load_npz(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_inputs_values.sparse.npz\")","metadata":{"execution":{"iopub.status.busy":"2022-09-09T15:45:26.246446Z","iopub.execute_input":"2022-09-09T15:45:26.249016Z","iopub.status.idle":"2022-09-09T15:46:30.588370Z","shell.execute_reply.started":"2022-09-09T15:45:26.248950Z","shell.execute_reply":"2022-09-09T15:46:30.586349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## PCA / TruncatedSVD\nIt is not possible to directly apply PCA to a sparse matrix, because PCA has to first \"center\" the data, which destroys the sparsity. This is why we apply `TruncatedSVD` instead (which is pretty much \"PCA without centering\"). It might be better to normalize the data a bit more here, but we will keep it simple.","metadata":{}},{"cell_type":"code","source":"%%time\npca = TruncatedSVD(n_components = 128, random_state = 42)\ntrain_inputs = pca.fit_transform(train_inputs)\n\nwith open(\"multiome_train_x.pickle\", \"wb\") as f:\n    pickle.dump(train_inputs , f)\ndel train_inputs\ngc.collect()\n\nmulti_test_x = scipy.sparse.load_npz(\"../input/multimodal-single-cell-as-sparse-matrix/test_multi_inputs_values.sparse.npz\")\nmulti_test_x = pca.transform(multi_test_x)\n\nwith open(\"multiome_test_x.pickle\", \"wb\") as f:\n    pickle.dump(multi_test_x, f, protocol=pickle.HIGHEST_PROTOCOL)\ndel multi_test_x\ngc.collect()\n\nwith open(\"pca.pickle\", \"wb\") as f:\n    pickle.dump(pca, f, protocol=pickle.HIGHEST_PROTOCOL)\n\ndel pca\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-09T15:46:30.590452Z","iopub.execute_input":"2022-09-09T15:46:30.590944Z","iopub.status.idle":"2022-09-09T16:11:51.339644Z","shell.execute_reply.started":"2022-09-09T15:46:30.590899Z","shell.execute_reply":"2022-09-09T16:11:51.338418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Random row selection and conversion of the target data to a dense matrix\n\nUnfortunately, although sklearn's `Ridge` regressor do accept sparse matrices as input, it does not accept sparse matrices as target values. This means we will have to convert the targets to a dense format. Although we could fit in memory both the dense target data and the sparse input data, the Ridge regression process would then lack memory. Therefore, from now on, we will work with a subset of 50 000 rows from the training data.","metadata":{}},{"cell_type":"code","source":"%%time\ntrain_target = scipy.sparse.load_npz(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_targets_values.sparse.npz\")","metadata":{"execution":{"iopub.status.busy":"2022-09-09T16:11:51.342713Z","iopub.execute_input":"2022-09-09T16:11:51.343210Z","iopub.status.idle":"2022-09-09T16:12:16.687531Z","shell.execute_reply.started":"2022-09-09T16:11:51.343154Z","shell.execute_reply":"2022-09-09T16:12:16.686162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_target = train_target.todense()\n# gc.collect()\noutput_pca = TruncatedSVD(n_components = 128, random_state = 42)\ntrain_target = output_pca.fit_transform(train_target)\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-09T16:12:16.690575Z","iopub.execute_input":"2022-09-09T16:12:16.691250Z","iopub.status.idle":"2022-09-09T16:18:08.635395Z","shell.execute_reply.started":"2022-09-09T16:12:16.691143Z","shell.execute_reply":"2022-09-09T16:18:08.634254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def correlation_score(y_true, y_pred):\n    \"\"\"Scores the predictions according to the competition rules. \n    \n    It is assumed that the predictions are not constant.\n    \n    Returns the average of each sample's Pearson correlation coefficient\"\"\"\n    if type(y_true) == pd.DataFrame: y_true = y_true.values\n    if type(y_pred) == pd.DataFrame: y_pred = y_pred.values\n    if y_true.shape != y_pred.shape: raise ValueError(\"Shapes are different.\")\n    corrsum = 0\n    for i in range(len(y_true)):\n        corrsum += np.corrcoef(y_true[i], y_pred[i])[1, 0]\n    return corrsum / len(y_true)\n\noriginal_train_target = scipy.sparse.load_npz(\n    \"../input/multimodal-single-cell-as-sparse-matrix/train_multi_targets_values.sparse.npz\"\n)[:5000]\ncorrelation_score(output_pca.inverse_transform(train_target[:5000]), original_train_target.todense())","metadata":{"execution":{"iopub.status.busy":"2022-09-09T16:18:08.636995Z","iopub.execute_input":"2022-09-09T16:18:08.637718Z","iopub.status.idle":"2022-09-09T16:18:28.751727Z","shell.execute_reply.started":"2022-09-09T16:18:08.637667Z","shell.execute_reply":"2022-09-09T16:18:28.750429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nwith open(\"multiome_train_y.pickle\", \"wb\") as f:\n    pickle.dump(train_target, f, protocol=pickle.HIGHEST_PROTOCOL)\ndel train_target\ngc.collect()\n\nwith open(\"output_pca.pickle\", \"wb\") as f:\n    pickle.dump(output_pca, f, protocol=pickle.HIGHEST_PROTOCOL)\ndel output_pca\ngc.collect()\n","metadata":{"execution":{"iopub.status.busy":"2022-09-09T16:18:28.753561Z","iopub.execute_input":"2022-09-09T16:18:28.753967Z","iopub.status.idle":"2022-09-09T16:18:29.150995Z","shell.execute_reply.started":"2022-09-09T16:18:28.753929Z","shell.execute_reply":"2022-09-09T16:18:29.149858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}