{"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).\n\nThe CITE predictions from this NB: https://www.kaggle.com/code/jsmithperera/msci-citeseq-quickstart-v3","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,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The scoring function (from AmbrosM)\n\nThis competition has a special metric: For every row, it computes the Pearson correlation between y_true and y_pred, and then all these correlation coefficients are averaged.","metadata":{}},{"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","metadata":{"execution":{"iopub.status.busy":"2022-09-02T18:36:36.503856Z","iopub.execute_input":"2022-09-02T18:36:36.504273Z","iopub.status.idle":"2022-09-02T18:36:36.513885Z","shell.execute_reply.started":"2022-09-02T18:36:36.504241Z","shell.execute_reply":"2022-09-02T18:36:36.512028Z"},"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-02T18:36:44.463945Z","iopub.execute_input":"2022-09-02T18:36:44.464372Z","iopub.status.idle":"2022-09-02T18:37:47.462308Z","shell.execute_reply.started":"2022-09-02T18:36:44.464333Z","shell.execute_reply":"2022-09-02T18:37:47.460478Z"},"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 = 32, random_state = 1)\ntrain_inputs = pca.fit_transform(train_inputs)","metadata":{"execution":{"iopub.status.busy":"2022-09-02T18:39:02.827816Z","iopub.execute_input":"2022-09-02T18:39:02.828289Z","iopub.status.idle":"2022-09-02T18:42:10.172824Z","shell.execute_reply.started":"2022-09-02T18:39:02.828247Z","shell.execute_reply":"2022-09-02T18:42:10.171561Z"},"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":"np.random.seed(42)\nall_row_indices = np.arange(train_inputs.shape[0])\nnp.random.shuffle(all_row_indices)\nselected_rows_indices = all_row_indices[:50000]","metadata":{"execution":{"iopub.status.busy":"2022-09-02T13:45:56.014400Z","iopub.execute_input":"2022-09-02T13:45:56.014885Z","iopub.status.idle":"2022-09-02T13:45:56.026551Z","shell.execute_reply.started":"2022-09-02T13:45:56.014845Z","shell.execute_reply":"2022-09-02T13:45:56.025560Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_inputs = train_inputs[selected_rows_indices]","metadata":{"execution":{"iopub.status.busy":"2022-09-02T13:46:01.307114Z","iopub.execute_input":"2022-09-02T13:46:01.307749Z","iopub.status.idle":"2022-09-02T13:46:01.317838Z","shell.execute_reply.started":"2022-09-02T13:46:01.307713Z","shell.execute_reply":"2022-09-02T13:46:01.316509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-02T13:46:05.874345Z","iopub.execute_input":"2022-09-02T13:46:05.875166Z","iopub.status.idle":"2022-09-02T13:46:32.083705Z","shell.execute_reply.started":"2022-09-02T13:46:05.875129Z","shell.execute_reply":"2022-09-02T13:46:32.082016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_target = train_target[selected_rows_indices]\ntrain_target = train_target.todense()\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-02T13:46:41.298748Z","iopub.execute_input":"2022-09-02T13:46:41.299157Z","iopub.status.idle":"2022-09-02T13:46:49.225947Z","shell.execute_reply.started":"2022-09-02T13:46:41.299123Z","shell.execute_reply":"2022-09-02T13:46:49.224788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## KFold Ridge regression\n`sklearn` complains that we should use array instead of matrices. Unfortunately, the old `scipy` version available on kaggle do not provide sparse arrays; only sparse matrices. So we suppress the warnings.","metadata":{}},{"cell_type":"code","source":"import warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)","metadata":{"execution":{"iopub.status.busy":"2022-09-02T13:46:53.485367Z","iopub.execute_input":"2022-09-02T13:46:53.485859Z","iopub.status.idle":"2022-09-02T13:46:53.492065Z","shell.execute_reply.started":"2022-09-02T13:46:53.485821Z","shell.execute_reply":"2022-09-02T13:46:53.490723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This Kfold ridge regression code is mostly taken from AmbrosM's [notebook](https://www.kaggle.com/code/ambrosm/msci-multiome-quickstart). Note that `sklearn`'s `Ridge` handles sparse matrices transparently. I found [this blog post](https://dziganto.github.io/Sparse-Matrices-For-Efficient-Machine-Learning/) that list the other algorithms of `sklearn` that accept sparse matrices.","metadata":{}},{"cell_type":"code","source":"%%time\n# Cross-validation\n\nkf = KFold(n_splits=3, shuffle=True, random_state=1)\nscore_list = []\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(train_inputs)):\n    model = None\n    gc.collect()\n    X_tr = train_inputs[idx_tr] # creates a copy, https://numpy.org/doc/stable/user/basics.copies.html\n    y_tr = train_target[idx_tr]\n    del idx_tr\n\n    model = Ridge(copy_X=False)\n    model.fit(X_tr, y_tr)\n    del X_tr, y_tr\n    gc.collect()\n\n    # We validate the model\n    X_va = train_inputs[idx_va]\n    y_va = train_target[idx_va]\n    del idx_va\n    y_va_pred = model.predict(X_va)\n    mse = mean_squared_error(y_va, y_va_pred)\n    corrscore = correlation_score(y_va, y_va_pred)\n    del X_va, y_va\n\n    print(f\"Fold {fold}: mse = {mse:.5f}, corr =  {corrscore:.3f}\")\n    score_list.append((mse, corrscore))\n\n# Show overall score\nresult_df = pd.DataFrame(score_list, columns=['mse', 'corrscore'])\nprint(f\"{Fore.GREEN}{Style.BRIGHT}{train_inputs.shape} Average  mse = {result_df.mse.mean():.5f}; corr = {result_df.corrscore.mean():.3f}{Style.RESET_ALL}\")\n","metadata":{"execution":{"iopub.status.busy":"2022-09-02T13:47:00.371305Z","iopub.execute_input":"2022-09-02T13:47:00.371794Z","iopub.status.idle":"2022-09-02T13:47:53.518035Z","shell.execute_reply.started":"2022-09-02T13:47:00.371756Z","shell.execute_reply":"2022-09-02T13:47:53.516334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Retraining\n","metadata":{}},{"cell_type":"code","source":"# We retrain the model and then delete the training data, which is no longer needed\nmodel, score_list, result_df = None, None, None # free the RAM occupied by the old model\ngc.collect()\nmodel = Ridge(copy_X=False) # we overwrite the training data\nmodel.fit(train_inputs, train_target)\n","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:06:45.810946Z","iopub.execute_input":"2022-08-30T21:06:45.811807Z","iopub.status.idle":"2022-08-30T21:06:49.967731Z","shell.execute_reply.started":"2022-08-30T21:06:45.811767Z","shell.execute_reply":"2022-08-30T21:06:49.966426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train_inputs, train_target # free the RAM\n_ = gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:06:49.970122Z","iopub.execute_input":"2022-08-30T21:06:49.971258Z","iopub.status.idle":"2022-08-30T21:06:50.429061Z","shell.execute_reply.started":"2022-08-30T21:06:49.971205Z","shell.execute_reply":"2022-08-30T21:06:50.427586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predicting","metadata":{}},{"cell_type":"code","source":"%%time\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)\ntest_pred = model.predict(multi_test_x)\ndel multi_test_x\ngc.collect()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Creating submission\n\nWe load the cells that will have to appear in submission.","metadata":{}},{"cell_type":"code","source":"%%time\n# Read the table of rows and columns required for submission\neval_ids = pd.read_parquet(\"../input/multimodal-single-cell-as-sparse-matrix/evaluation.parquet\")\n\n# Convert the string columns to more efficient categorical types\n#eval_ids.cell_id = eval_ids.cell_id.apply(lambda s: int(s, base=16))\n\neval_ids.cell_id = eval_ids.cell_id.astype(pd.CategoricalDtype())\neval_ids.gene_id = eval_ids.gene_id.astype(pd.CategoricalDtype())\n","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:07:06.485029Z","iopub.execute_input":"2022-08-30T21:07:06.485517Z","iopub.status.idle":"2022-08-30T21:07:26.212989Z","shell.execute_reply.started":"2022-08-30T21:07:06.485471Z","shell.execute_reply":"2022-08-30T21:07:26.211556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prepare an empty series which will be filled with predictions\nsubmission = pd.Series(name='target',\n                       index=pd.MultiIndex.from_frame(eval_ids), \n                       dtype=np.float32)\nsubmission","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:07:26.222677Z","iopub.execute_input":"2022-08-30T21:07:26.223535Z","iopub.status.idle":"2022-08-30T21:07:49.40937Z","shell.execute_reply.started":"2022-08-30T21:07:26.223441Z","shell.execute_reply":"2022-08-30T21:07:49.408498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We load the `index`  and `columns` of the original dataframe, as we need them to make the submission.","metadata":{}},{"cell_type":"code","source":"%%time\ny_columns = np.load(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_targets_idxcol.npz\",\n                   allow_pickle=True)[\"columns\"]\n\ntest_index = np.load(\"../input/multimodal-single-cell-as-sparse-matrix/test_multi_inputs_idxcol.npz\",\n                    allow_pickle=True)[\"index\"]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We assign the predicted values to the correct row in the submission file.","metadata":{}},{"cell_type":"code","source":"cell_dict = dict((k,v) for v,k in enumerate(test_index)) \nassert len(cell_dict)  == len(test_index)\n\ngene_dict = dict((k,v) for v,k in enumerate(y_columns))\nassert len(gene_dict) == len(y_columns)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:08:41.609646Z","iopub.execute_input":"2022-08-30T21:08:41.610028Z","iopub.status.idle":"2022-08-30T21:08:41.634105Z","shell.execute_reply.started":"2022-08-30T21:08:41.609997Z","shell.execute_reply":"2022-08-30T21:08:41.632255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"eval_ids_cell_num = eval_ids.cell_id.apply(lambda x:cell_dict.get(x, -1))\neval_ids_gene_num = eval_ids.gene_id.apply(lambda x:gene_dict.get(x, -1))\n\nvalid_multi_rows = (eval_ids_gene_num !=-1) & (eval_ids_cell_num!=-1)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:08:41.679451Z","iopub.execute_input":"2022-08-30T21:08:41.679803Z","iopub.status.idle":"2022-08-30T21:08:42.721245Z","shell.execute_reply.started":"2022-08-30T21:08:41.679773Z","shell.execute_reply":"2022-08-30T21:08:42.719597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.iloc[valid_multi_rows] = test_pred[eval_ids_cell_num[valid_multi_rows].to_numpy(),\neval_ids_gene_num[valid_multi_rows].to_numpy()]","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:08:43.976624Z","iopub.execute_input":"2022-08-30T21:08:43.977165Z","iopub.status.idle":"2022-08-30T21:08:46.309975Z","shell.execute_reply.started":"2022-08-30T21:08:43.9771Z","shell.execute_reply":"2022-08-30T21:08:46.308794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del eval_ids_cell_num, eval_ids_gene_num, valid_multi_rows, eval_ids, test_index, y_columns\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:08:46.311728Z","iopub.execute_input":"2022-08-30T21:08:46.312201Z","iopub.status.idle":"2022-08-30T21:08:46.451644Z","shell.execute_reply.started":"2022-08-30T21:08:46.312134Z","shell.execute_reply":"2022-08-30T21:08:46.449967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:08:46.455478Z","iopub.execute_input":"2022-08-30T21:08:46.455901Z","iopub.status.idle":"2022-08-30T21:08:46.46827Z","shell.execute_reply.started":"2022-08-30T21:08:46.455866Z","shell.execute_reply":"2022-08-30T21:08:46.467247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Merging with CITEseq predictions\n\nWe use the CITEseq predictions from [this notebook](https://www.kaggle.com/code/vuonglam/lgbm-baseline-optuna-drop-constant-cite-task) by VuongLam.","metadata":{}},{"cell_type":"code","source":"submission.reset_index(drop = True, inplace = True)\nsubmission.index.name = 'row_id'\n\n# with open(\"partial_submission_multi.pickle\", 'wb') as f:\n#     pickle.dump(submission, f)\n# submission","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:08:46.469655Z","iopub.execute_input":"2022-08-30T21:08:46.470305Z","iopub.status.idle":"2022-08-30T21:08:46.684431Z","shell.execute_reply.started":"2022-08-30T21:08:46.470257Z","shell.execute_reply":"2022-08-30T21:08:46.683314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#cite_submission = pd.read_csv(\"../input/lgbm-baseline-optuna-drop-constant-cite-task/submission.csv\")\n\nwith open(\"../input/msci-citeseg-by-cell-type/citeseq_pred_by_cell.pickle\", 'rb') as f: cite_pred = pickle.load(f)\n    \n#cite_submission = cite_submission.set_index(\"row_id\")\n#cite_submission = cite_submission[\"target\"]","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:08:46.685902Z","iopub.execute_input":"2022-08-30T21:08:46.686539Z","iopub.status.idle":"2022-08-30T21:09:09.624762Z","shell.execute_reply.started":"2022-08-30T21:08:46.686505Z","shell.execute_reply":"2022-08-30T21:09:09.623247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#submission[submission.isnull()] = cite_submission[submission.isnull()]\nsubmission.iloc[:len(cite_pred.ravel())] = cite_pred.ravel()","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:09:09.626409Z","iopub.execute_input":"2022-08-30T21:09:09.626837Z","iopub.status.idle":"2022-08-30T21:09:11.581751Z","shell.execute_reply.started":"2022-08-30T21:09:09.626806Z","shell.execute_reply":"2022-08-30T21:09:11.580764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:09:11.583267Z","iopub.execute_input":"2022-08-30T21:09:11.584331Z","iopub.status.idle":"2022-08-30T21:09:11.592472Z","shell.execute_reply.started":"2022-08-30T21:09:11.584292Z","shell.execute_reply":"2022-08-30T21:09:11.59144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.isnull().any()","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:09:11.59419Z","iopub.execute_input":"2022-08-30T21:09:11.595001Z","iopub.status.idle":"2022-08-30T21:09:11.655Z","shell.execute_reply.started":"2022-08-30T21:09:11.594949Z","shell.execute_reply":"2022-08-30T21:09:11.65403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv(\"submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:09:11.656341Z","iopub.execute_input":"2022-08-30T21:09:11.657248Z","iopub.status.idle":"2022-08-30T21:11:24.259006Z","shell.execute_reply.started":"2022-08-30T21:09:11.657213Z","shell.execute_reply":"2022-08-30T21:11:24.257429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!head submission.csv","metadata":{"execution":{"iopub.status.busy":"2022-08-30T21:11:24.26136Z","iopub.execute_input":"2022-08-30T21:11:24.261899Z","iopub.status.idle":"2022-08-30T21:11:25.569433Z","shell.execute_reply.started":"2022-08-30T21:11:24.261849Z","shell.execute_reply":"2022-08-30T21:11:25.56775Z"},"trusted":true},"execution_count":null,"outputs":[]}]}