{"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":"# Multiome Quickstart\n\nThis notebook shows how to cross-validate a baseline model and create a submission for the Multiome part of the *Multimodal Single-Cell Integration* competition without running out of memory.\n\nIt does not show the EDA - see the separate notebook [MSCI EDA which makes sense ⭐️⭐️⭐️⭐️⭐️](https://www.kaggle.com/ambrosm/msci-eda-which-makes-sense).\n\nThe baseline model for the other part of the competition (CITEseq) is [here](https://www.kaggle.com/ambrosm/msci-citeseq-quickstart).","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"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\nfrom sklearn.dummy import DummyRegressor\nfrom sklearn.pipeline import make_pipeline, Pipeline\nfrom sklearn.linear_model import Ridge, LinearRegression\nfrom sklearn.metrics import mean_squared_error\n\nDATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\nFP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\nFP_CITE_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_cite_inputs.h5\")\nFP_CITE_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_cite_targets.h5\")\nFP_CITE_TEST_INPUTS = os.path.join(DATA_DIR,\"test_cite_inputs.h5\")\n\nFP_MULTIOME_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_multi_inputs.h5\")\nFP_MULTIOME_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_multi_targets.h5\")\nFP_MULTIOME_TEST_INPUTS = os.path.join(DATA_DIR,\"test_multi_inputs.h5\")\n\nFP_SUBMISSION = os.path.join(DATA_DIR,\"sample_submission.csv\")\nFP_EVALUATION_IDS = os.path.join(DATA_DIR,\"evaluation_ids.csv\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-20T22:40:44.879348Z","iopub.execute_input":"2022-08-20T22:40:44.879868Z","iopub.status.idle":"2022-08-20T22:40:44.894334Z","shell.execute_reply.started":"2022-08-20T22:40:44.879828Z","shell.execute_reply":"2022-08-20T22:40:44.893012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# If you see a warning \"Failed to establish a new connection\" running this cell,\n# go to \"Settings\" on the right hand side, \n# and turn on internet. Note, you need to be phone verified.\n# We need this library to read HDF files.\n!pip install --quiet tables\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-20T22:40:44.907771Z","iopub.execute_input":"2022-08-20T22:40:44.909052Z","iopub.status.idle":"2022-08-20T22:40:56.401124Z","shell.execute_reply.started":"2022-08-20T22:40:44.908998Z","shell.execute_reply":"2022-08-20T22:40:56.399328Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loading the common metadata table\n\nThe current version of the model is so primitive that it doesn't use the metadata, but we load it anyway.","metadata":{}},{"cell_type":"code","source":"df_cell = pd.read_csv(FP_CELL_METADATA)\ndf_cell_cite = df_cell[df_cell.technology==\"citeseq\"]\ndf_cell_multi = df_cell[df_cell.technology==\"multiome\"]\ndf_cell_cite.shape, df_cell_multi.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-20T22:40:56.404482Z","iopub.execute_input":"2022-08-20T22:40:56.404925Z","iopub.status.idle":"2022-08-20T22:40:56.752398Z","shell.execute_reply.started":"2022-08-20T22:40:56.404858Z","shell.execute_reply":"2022-08-20T22:40:56.751124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The scoring function\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-08-20T22:40:56.754781Z","iopub.execute_input":"2022-08-20T22:40:56.755283Z","iopub.status.idle":"2022-08-20T22:40:56.763769Z","shell.execute_reply.started":"2022-08-20T22:40:56.755244Z","shell.execute_reply":"2022-08-20T22:40:56.762392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessing and cross-validation\n\nThe Multiome dataset is way too large to fit into 16 GByte RAM:\n- train inputs:  105942 * 228942 float32 values (97 GByte)\n- train targets: 105942 *  23418 float32 values (10 GByte)\n- test inputs:    55935 * 228942 float32 values (13 GByte)\n\nTo get a result with only 16 GByte RAM, we simplify the problem as follows:\n- We ignore the complete metadata (donors, days, cell types).\n- We read only 6000 rows of the training data.\n- We drop all feature columns which are constant.\n- Of the remaining columns, we keep only 4000.\n- We do a PCA and keep only the 4 most important components.\n- We fit a ridge regression model with 6000\\*4 inputs and 6000\\*23418 targets.","metadata":{}},{"cell_type":"code","source":"#%%time\n# Preprocessing\n\nclass PreprocessMultiome(BaseEstimator, TransformerMixin):\n    columns_to_use = slice(10000, 14000)\n    \n    @staticmethod\n    def take_column_subset(X):\n        return X[:,PreprocessMultiome.columns_to_use]\n    \n    def transform(self, X):\n        print(X.shape)\n        X = X[:,~self.all_zero_columns]\n        print(X.shape)\n        X = PreprocessMultiome.take_column_subset(X) # use only a part of the columns\n        print(X.shape)\n        gc.collect()\n\n        X = self.pca.transform(X)\n        print(X.shape)\n        return X\n\n    def fit_transform(self, X):\n        print(X.shape)\n        self.all_zero_columns = (X == 0).all(axis=0)\n        X = X[:,~self.all_zero_columns]\n        print(X.shape)\n        X = PreprocessMultiome.take_column_subset(X) # use only a part of the columns\n        print(X.shape)\n        gc.collect()\n\n        self.pca = PCA(n_components=4, copy=False, random_state=1)\n        X = self.pca.fit_transform(X)\n        plt.plot(self.pca.explained_variance_ratio_.cumsum())\n        plt.title(\"Cumulative explained variance ratio\")\n        plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))\n        plt.xlabel('PCA component')\n        plt.ylabel('Cumulative explained variance ratio')\n        plt.show()\n        print(X.shape)\n        return X\n\npreprocessor = PreprocessMultiome()\n\nmulti_train_x = None\nstart, stop = 0, 6000\nmulti_train_x = preprocessor.fit_transform(pd.read_hdf(FP_MULTIOME_TRAIN_INPUTS, start=start, stop=stop).values)\n\nmulti_train_y = pd.read_hdf(FP_MULTIOME_TRAIN_TARGETS, start=start, stop=stop)\ny_columns = multi_train_y.columns\nmulti_train_y = multi_train_y.values\nprint(multi_train_y.shape)","metadata":{"execution":{"iopub.status.busy":"2022-08-20T22:40:56.768404Z","iopub.execute_input":"2022-08-20T22:40:56.768973Z","iopub.status.idle":"2022-08-20T22:40:58.921768Z","shell.execute_reply.started":"2022-08-20T22:40:56.768922Z","shell.execute_reply":"2022-08-20T22:40:58.920440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Cross-validation\n\nkf = KFold(n_splits=5, shuffle=True, random_state=1)\nscore_list = []\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(multi_train_x)):\n    model = None\n    gc.collect()\n    X_tr = multi_train_x[idx_tr] # creates a copy, https://numpy.org/doc/stable/user/basics.copies.html\n    y_tr = multi_train_y[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 = multi_train_x[idx_va]\n    y_va = multi_train_y[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}{multi_train_x.shape} Average  mse = {result_df.mse.mean():.5f}; corr = {result_df.corrscore.mean():.3f}{Style.RESET_ALL}\")\n","metadata":{"execution":{"iopub.status.busy":"2022-08-20T22:40:58.923484Z","iopub.execute_input":"2022-08-20T22:40:58.924842Z","iopub.status.idle":"2022-08-20T22:41:04.399786Z","shell.execute_reply.started":"2022-08-20T22:40:58.924786Z","shell.execute_reply":"2022-08-20T22:41:04.398080Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"By the way, this ridge regression is not much better than DummyRegressor, which scores `mse = 2.01718; corr = 0.679`.","metadata":{}},{"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(multi_train_x, multi_train_y)\ndel multi_train_x, multi_train_y # free the RAM\n_ = gc.collect()\n","metadata":{"execution":{"iopub.status.busy":"2022-08-20T22:41:04.402167Z","iopub.execute_input":"2022-08-20T22:41:04.403562Z","iopub.status.idle":"2022-08-20T22:41:05.492237Z","shell.execute_reply.started":"2022-08-20T22:41:04.403490Z","shell.execute_reply":"2022-08-20T22:41:05.490742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The final submission will contain 65744180 predictions, of which the first 6812820 are CITEseq predictions and the remaining 58931360 are Multiome. \n\nThe Multiome test predictions have 55935 rows and 23418 columns. 55935 \\* 23418 = 1’309’885’830 predictions. We'll only submit 4.5 % of these predictions. According to the data description, this subset was created by sampling 30 % of the Multiome rows, and for each row, 15 % of the columns (i.e., 16780 rows and 3512 columns per row). Consequently, when reading the test data, we can immediately drop 70 % of the rows and keep only the remaining 16780.\n\nThe eval_ids table specifies which predictions are required for the submission file.","metadata":{}},{"cell_type":"code","source":"%%time\n# Read the table of rows and columns required for submission\neval_ids = pd.read_csv(FP_EVALUATION_IDS, index_col='row_id')\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))\neval_ids.cell_id = eval_ids.cell_id.astype(pd.CategoricalDtype())\neval_ids.gene_id = eval_ids.gene_id.astype(pd.CategoricalDtype())\ndisplay(eval_ids)\n\n# Create the set of needed cell_ids\ncell_id_set = set(eval_ids.cell_id)\n\n# Convert the string gene_ids to a more efficient categorical dtype\ny_columns = pd.CategoricalIndex(y_columns, dtype=eval_ids.gene_id.dtype, name='gene_id')\n","metadata":{"execution":{"iopub.status.busy":"2022-08-20T22:41:05.494037Z","iopub.execute_input":"2022-08-20T22:41:05.494613Z","iopub.status.idle":"2022-08-20T22:43:11.300770Z","shell.execute_reply.started":"2022-08-20T22:41:05.494561Z","shell.execute_reply":"2022-08-20T22:43:11.299312Z"},"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-20T22:43:11.303536Z","iopub.execute_input":"2022-08-20T22:43:11.303958Z","iopub.status.idle":"2022-08-20T22:43:11.501091Z","shell.execute_reply.started":"2022-08-20T22:43:11.303920Z","shell.execute_reply":"2022-08-20T22:43:11.499920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now compute the predictions in chunks of 5000 rows and match them with the eval_ids table row by row. The matching is very slow, but space-efficient.","metadata":{}},{"cell_type":"code","source":"%%time\n# Process the test data in chunks of 5000 rows\n\nstart = 0\nchunksize = 5000\ntotal_rows = 0\nwhile True:\n    multi_test_x = None # Free the memory if necessary\n    gc.collect()\n    # Read the 5000 rows and select the 30 % subset which is needed for the submission\n    multi_test_x = pd.read_hdf(FP_MULTIOME_TEST_INPUTS, start=start, stop=start+chunksize)\n    rows_read = len(multi_test_x)\n    needed_row_mask = multi_test_x.index.isin(cell_id_set)\n    multi_test_x = multi_test_x.loc[needed_row_mask]\n    \n    # Keep the index (the cell_ids) for later\n    multi_test_index = multi_test_x.index\n    \n    # Predict\n    multi_test_x = multi_test_x.values\n    multi_test_x = preprocessor.transform(multi_test_x)\n    test_pred = model.predict(multi_test_x)\n    \n    # Convert the predictions to a dataframe so that they can be matched with eval_ids\n    test_pred = pd.DataFrame(test_pred,\n                             index=pd.CategoricalIndex(multi_test_index,\n                                                       dtype=eval_ids.cell_id.dtype,\n                                                       name='cell_id'),\n                             columns=y_columns)\n    gc.collect()\n    \n    # Fill the predictions into the submission series row by row\n    for i, (index, row) in enumerate(test_pred.iterrows()):\n        row = row.reindex(eval_ids.gene_id[eval_ids.cell_id == index])\n        submission.loc[index] = row.values\n    print('na:', submission.isna().sum())\n\n    #test_pred_list.append(test_pred)\n    total_rows += len(multi_test_x)\n    print(total_rows)\n    if rows_read < chunksize: break # this was the last chunk\n    start += chunksize\n    \ndel multi_test_x, multi_test_index, needed_row_mask\n","metadata":{"execution":{"iopub.status.busy":"2022-08-20T22:49:34.321278Z","iopub.execute_input":"2022-08-20T22:49:34.321748Z","iopub.status.idle":"2022-08-20T23:51:50.284303Z","shell.execute_reply.started":"2022-08-20T22:49:34.321711Z","shell.execute_reply":"2022-08-20T23:51:50.281608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission\n\nAs we don't yet have the CITEseq predictions, we save the partial predictions so that they can be used in the [CITEseq notebook](https://www.kaggle.com/ambrosm/msci-citeseq-quickstart).","metadata":{}},{"cell_type":"code","source":"submission.reset_index(drop=True, inplace=True)\nsubmission.index.name = 'row_id'\nwith open(\"partial_submission_multi.pickle\", 'wb') as f: pickle.dump(submission, f)\nsubmission","metadata":{"execution":{"iopub.status.busy":"2022-08-20T23:53:51.263900Z","iopub.execute_input":"2022-08-20T23:53:51.265185Z","iopub.status.idle":"2022-08-20T23:53:52.528470Z","shell.execute_reply.started":"2022-08-20T23:53:51.265109Z","shell.execute_reply":"2022-08-20T23:53:52.527117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}