{"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":"# Singular Value Decomposition for Multiomer Prediction","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"## 1. Introduction","metadata":{}},{"cell_type":"markdown","source":"In this work, we developed a Singular Value Decomposition (SVD) approach to predict the gene expression level of [the Kaggle Open Problem Multimodal competition](https://www.kaggle.com/competitions/open-problems-multimodal).\n\nDue to large-scale nature of this problem, SVD decomposition has been very popular for this competition, particularly for the Multiomer part, where it becomes essentially impossible to load all the data into the computer memory. Further, huge size of the data also prevents statistical models from completing in time. Therefore, a few works, [this one](https://www.kaggle.com/code/xiafire/lb-t15-msci-multiome-catboostregressor), [this one](https://www.kaggle.com/code/fabiencrom/msci-multiome-torch-quickstart-submission), [this one](https://www.kaggle.com/code/jsmithperera/multiome-w-sparse-m-tsvd-32) and [this one](https://www.kaggle.com/code/ambrosm/msci-multiome-quickstart) all either use SVD to decompose the input matrix or both the input and target matrix.\n\nIn this work, we will develop a novel approach to incorporate SVD and statistical algorithms. The contributions of this work are \n\n1. We systematically explore the tradeoff between more SVD components and model accuracy through cross-validation.\n\n2. We develop a highly efficient (both in CPU and memory footprint) and intuitive approach to generate the testing submission final. This is not a trivial task due to the size of the problem and the complexity of sampled entries for submission. Current approach is through a highly complex and proprietory algorithm to format the data. On the contrary, our approach makes clear uses of python dictionary and numpy vectorized functions to complete the work within a pandas table.\n\nThe rest of the notebook is organized as follows:\n\n1. We demonstrate SVD on both input and target matrics in Section 2.\n\n2. We show how to build a cross-validation in Section 3.\n\n3. We show how to prepare testing data set and generate prediction in Section 4.\n\n4. We present our work to format the prediction result in Section 5.\n\n5. Finally, we conclude this work and present some closing thoughts in Section 6.","metadata":{}},{"cell_type":"markdown","source":"## 2. SVD Decomposition","metadata":{}},{"cell_type":"markdown","source":"Singular Value Decomposition is a mathematical technique widely used in various Science and Engineering applications. It aims at decomposing a rectangular matrix $A$ (of size $m$ by $n$) into a product of three matrices: $A = UDV^T$ where $V^T$ is the transpose of $V$. Here, the matrix $U$ is of size $m$ by $m$, matrix $D$ is of size $m$ by $n$, and matrix $V$ is of size $n$ by $n$. Highly efficient algorithm and implementation exists to achieve it.\n\nA closely related algorithm, Truncated Singular Value Decomposition, uses only the first few SVD components and is widely used for applications involves large-scale matrices (e.g., [this notebook](https://www.kaggle.com/code/stautxie/multiomer-prediction-with-svd-regression) and [this notebook](https://www.kaggle.com/code/stautxie/cite-seq-baseline-prediction)). \n\nInterested readers can visit these introductory pages [here](https://en.wikipedia.org/wiki/Singular_value_decomposition) and [here](https://langvillea.people.cofc.edu/DISSECTION-LAB/Emmie'sLSI-SVDModule/p5module.html) for more details.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport scipy\nfrom scipy.sparse import csr_matrix, csc_matrix, vstack, save_npz, load_npz\nfrom scipy.stats import pearsonr\nfrom sklearn.multioutput import MultiOutputRegressor\nfrom sklearn.linear_model import Ridge, Lasso\nfrom sklearn.model_selection import KFold\nfrom sklearn.decomposition import TruncatedSVD\nfrom sklearn.metrics import mean_absolute_percentage_error, mean_squared_error\nfrom sklearn.preprocessing import LabelEncoder\nimport pickle\nimport os\nimport random\nimport gc\nimport time\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:48:34.183867Z","iopub.execute_input":"2022-09-12T16:48:34.184778Z","iopub.status.idle":"2022-09-12T16:48:35.574679Z","shell.execute_reply.started":"2022-09-12T16:48:34.184327Z","shell.execute_reply":"2022-09-12T16:48:35.573478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"try:\n    import tables as tb\nexcept:\n    !pip install tables\n    import tables as tb","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:48:37.331926Z","iopub.execute_input":"2022-09-12T16:48:37.332366Z","iopub.status.idle":"2022-09-12T16:48:37.539545Z","shell.execute_reply.started":"2022-09-12T16:48:37.332332Z","shell.execute_reply":"2022-09-12T16:48:37.538437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DAT_DIR = '../input/open-problems-multimodal'\nOUT_DIR = \"../input/python-sparse-matrix-open-problem-multimodal\"\nSVD_DIR = '.'\nseed = 42","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:48:39.708467Z","iopub.execute_input":"2022-09-12T16:48:39.709303Z","iopub.status.idle":"2022-09-12T16:48:39.715647Z","shell.execute_reply.started":"2022-09-12T16:48:39.709255Z","shell.execute_reply":"2022-09-12T16:48:39.714087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will load the input data in sparse matrix form to minimize memory footprint. Interested readers can visit [the notebook](https://www.kaggle.com/code/stautxie/lower-memory-footprint-via-python-sparse-matrix/notebook) and the [dataset](https://www.kaggle.com/datasets/stautxie/python-sparse-matrix-open-problem-multimodal).","metadata":{}},{"cell_type":"code","source":"# tr_mu_inputs_sp = pickle.load(open(os.path.join(OUT_DIR, \"tr_mu_inputs_sp.pkl\"), \"rb\"))\ntr_mu_inputs_sp = load_npz(os.path.join(OUT_DIR, \"tr_mu_inputs_sp.npz\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T15:41:21.985658Z","iopub.execute_input":"2022-09-12T15:41:21.986107Z","iopub.status.idle":"2022-09-12T15:42:23.911848Z","shell.execute_reply.started":"2022-09-12T15:41:21.986072Z","shell.execute_reply":"2022-09-12T15:42:23.909443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will first apply SVD on the training input matrix.","metadata":{}},{"cell_type":"code","source":"svd_comp = 32\n\nt1 = time.time()\nif os.path.exists(os.path.join(SVD_DIR, f\"tr_mu_inputs_svd_{str(svd_comp)}.pkl\")):\n    t_svd = pickle.load(open(os.path.join(SVD_DIR, f\"tr_mu_inputs_svd_{str(svd_comp)}.pkl\"), \"rb\"))\nelse:\n    t_svd = TruncatedSVD(n_components = svd_comp, random_state = seed)\n    t_svd.fit(tr_mu_inputs_sp)\n    pickle.dump(t_svd, open(os.path.join(SVD_DIR, f\"tr_mu_inputs_svd_{str(svd_comp)}.pkl\"), \"wb\"))\n    \nprint(f'CPU time used : {time.time() - t1}')","metadata":{"execution":{"iopub.status.busy":"2022-09-12T15:43:34.340652Z","iopub.execute_input":"2022-09-12T15:43:34.341055Z","iopub.status.idle":"2022-09-12T15:51:01.211529Z","shell.execute_reply.started":"2022-09-12T15:43:34.341023Z","shell.execute_reply":"2022-09-12T15:51:01.208875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_mu_inputs_svd = t_svd.transform(tr_mu_inputs_sp)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:04:51.247701Z","iopub.execute_input":"2022-09-12T16:04:51.248235Z","iopub.status.idle":"2022-09-12T16:05:11.226378Z","shell.execute_reply.started":"2022-09-12T16:04:51.248193Z","shell.execute_reply":"2022-09-12T16:05:11.225085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del tr_mu_inputs_sp; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:05:20.705611Z","iopub.execute_input":"2022-09-12T16:05:20.706043Z","iopub.status.idle":"2022-09-12T16:05:20.933141Z","shell.execute_reply.started":"2022-09-12T16:05:20.706005Z","shell.execute_reply":"2022-09-12T16:05:20.931802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To save time, we track the cpu time usage of computing SVD components for 16, 32, 64, and 128 and the results are shown below. The computation time increases approximately linearly with increase of SVD components.","metadata":{}},{"cell_type":"code","source":"svd_input_cpu_df = pd.DataFrame({'svd_comp':[16, 32, 64, 128],\n                                'cpu_time':[459.94, 603.46, 939.77, 1749.85]})\n\ng = sns.relplot(data=svd_input_cpu_df,\n               x = 'svd_comp',\n               y = 'cpu_time',\n               kind = 'line').set(title=\"CPU time vs. SVD Components for Training Input Matrix\")","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:05:24.728441Z","iopub.execute_input":"2022-09-12T16:05:24.728856Z","iopub.status.idle":"2022-09-12T16:05:25.204647Z","shell.execute_reply.started":"2022-09-12T16:05:24.728823Z","shell.execute_reply":"2022-09-12T16:05:25.203221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can also apply SVD on the training target matrix. Doing so will save memory footprint as well as computation time as fitting more output takes more time.","metadata":{}},{"cell_type":"code","source":"# tr_mu_targets_sp = pickle.load(open(os.path.join(OUT_DIR, \"tr_mu_targets_sp.pkl\"), \"rb\"))\ntr_mu_targets_sp = load_npz(os.path.join(OUT_DIR, \"tr_mu_targets_sp.npz\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:05:30.089628Z","iopub.execute_input":"2022-09-12T16:05:30.090079Z","iopub.status.idle":"2022-09-12T16:05:56.285986Z","shell.execute_reply.started":"2022-09-12T16:05:30.090041Z","shell.execute_reply":"2022-09-12T16:05:56.284517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"svd_comp_targets = 32\n\nt1 = time.time()\nif os.path.exists(os.path.join(SVD_DIR, f\"tr_mu_targets_svd_{str(svd_comp_targets)}.pkl\")):\n    t_svd_targets = pickle.load(open(os.path.join(SVD_DIR, \n                                                  f\"tr_mu_targets_svd_{str(svd_comp_targets)}.pkl\"), \"rb\"))\nelse:\n    t_svd_targets = TruncatedSVD(n_components = svd_comp_targets, random_state=seed)\n    t_svd_targets.fit(tr_mu_targets_sp)\n    pickle.dump(t_svd_targets, open(os.path.join(SVD_DIR, \n                                                  f\"tr_mu_targets_svd_{str(svd_comp_targets)}.pkl\"), \"wb\"))\n    \nprint(f'CPU time used : {time.time() - t1}')    ","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:07:47.887249Z","iopub.execute_input":"2022-09-12T16:07:47.888441Z","iopub.status.idle":"2022-09-12T16:09:39.017657Z","shell.execute_reply.started":"2022-09-12T16:07:47.888376Z","shell.execute_reply":"2022-09-12T16:09:39.016181Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The computation time also increases approximately linearly with increase of SVD components.","metadata":{}},{"cell_type":"code","source":"svd_target_cpu_df = pd.DataFrame({'svd_comp':[16, 32, 64, 128],\n                                'cpu_time':[80.57, 121.43, 208.91, 380.32]})\n\ng = sns.relplot(data=svd_target_cpu_df,\n               x = 'svd_comp',\n               y = 'cpu_time',\n               kind = 'line').set(title=\"CPU time vs. SVD Components for Training Target Matrix\")","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:12:31.252084Z","iopub.execute_input":"2022-09-12T16:12:31.252618Z","iopub.status.idle":"2022-09-12T16:12:31.497956Z","shell.execute_reply.started":"2022-09-12T16:12:31.252582Z","shell.execute_reply":"2022-09-12T16:12:31.496764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_mu_targets_svd = t_svd_targets.transform(tr_mu_targets_sp)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:12:38.354251Z","iopub.execute_input":"2022-09-12T16:12:38.354754Z","iopub.status.idle":"2022-09-12T16:12:44.191999Z","shell.execute_reply.started":"2022-09-12T16:12:38.354712Z","shell.execute_reply":"2022-09-12T16:12:44.190944Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Cross-Validation Study","metadata":{}},{"cell_type":"markdown","source":"In this section, we will employ 5-fold CV to study the quality of fitting. More specifically, we will compare Ridge and Lasso regression as well as the impact of number of SVD components. We implement the competition metric, mean pearson correlation coefficient, in the function mean_corr() below","metadata":{}},{"cell_type":"code","source":"def mean_corr(Y_true, Y_pred):\n    n_row = Y_true.shape[0]\n    row_corr_list = map(lambda i: pearsonr(Y_true[i, :], Y_pred[i, :])[0],\n                        range(n_row))\n    row_corr_list = np.array(list(row_corr_list))\n    \n    return np.mean(list(row_corr_list))\n","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:13:01.818847Z","iopub.execute_input":"2022-09-12T16:13:01.819267Z","iopub.status.idle":"2022-09-12T16:13:01.826850Z","shell.execute_reply.started":"2022-09-12T16:13:01.819233Z","shell.execute_reply":"2022-09-12T16:13:01.825311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We also present comparison of Ridge vs. Lasso Regression below.","metadata":{}},{"cell_type":"markdown","source":"### 3.1 Ridge Regression","metadata":{}},{"cell_type":"code","source":"num_folds = 5\nkfolds = KFold(n_splits = num_folds, shuffle=True, random_state = seed)\n\ncv_list = []\nfor i, (tr_index, va_index) in enumerate(kfolds.split(tr_mu_inputs_svd)):\n    tr_X = tr_mu_inputs_svd[tr_index, :]\n    va_X = tr_mu_inputs_svd[va_index, :]\n    tr_Y = tr_mu_targets_svd[tr_index, :]\n    va_Y = tr_mu_targets_svd[va_index, :]\n    \n    model = MultiOutputRegressor(Ridge(alpha=1e1), n_jobs=5)\n    model.fit(tr_X, tr_Y)\n    va_pred_svd = model.predict(va_X)\n    va_pred = t_svd_targets.inverse_transform(va_pred_svd)\n    va_Y_true = tr_mu_targets_sp[va_index, :].toarray()\n    mse_val = np.sqrt(mean_squared_error(va_Y_true, va_pred))\n    mean_corr_val = mean_corr(va_Y_true, va_pred)\n    \n    print(f'Iteration = {i}, MSE = {mse_val:.3f}, mean_corr = {mean_corr_val:.6f}')\n    \n    del model, tr_X, va_X, tr_Y, va_Y, va_Y_true, va_pred; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:13:06.151808Z","iopub.execute_input":"2022-09-12T16:13:06.152233Z","iopub.status.idle":"2022-09-12T16:14:50.500057Z","shell.execute_reply.started":"2022-09-12T16:13:06.152201Z","shell.execute_reply":"2022-09-12T16:14:50.498930Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.2 Lasso Regression","metadata":{}},{"cell_type":"code","source":"num_folds = 5\nkfolds = KFold(n_splits = num_folds, shuffle=True, random_state = seed)\n\ncv_list = []\nfor i, (tr_index, va_index) in enumerate(kfolds.split(tr_mu_inputs_svd)):\n    tr_X = tr_mu_inputs_svd[tr_index, :]\n    va_X = tr_mu_inputs_svd[va_index, :]\n    tr_Y = tr_mu_targets_svd[tr_index, :]\n    va_Y = tr_mu_targets_svd[va_index, :]\n    \n    model = MultiOutputRegressor(Lasso(alpha=1e1), n_jobs=5)\n    model.fit(tr_X, tr_Y)\n    va_pred_svd = model.predict(va_X)\n    va_pred = t_svd_targets.inverse_transform(va_pred_svd)\n    va_Y_true = tr_mu_targets_sp[va_index, :].toarray()\n    mse_val = np.sqrt(mean_squared_error(va_Y_true, va_pred))\n    mean_corr_val = mean_corr(va_Y_true, va_pred)\n    \n    print(f'Iteration = {i}, MSE = {mse_val:.3f}, mean_corr = {mean_corr_val:.6f}')\n    \n    del model, tr_X, va_X, tr_Y, va_Y, va_Y_true, va_pred; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:19:27.878917Z","iopub.execute_input":"2022-09-12T16:19:27.879430Z","iopub.status.idle":"2022-09-12T16:21:14.401150Z","shell.execute_reply.started":"2022-09-12T16:19:27.879379Z","shell.execute_reply":"2022-09-12T16:21:14.400015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see the Ridge regression actually performs better than the Lasso and is different than [the conclusion from our previous experiment with an R implementation](https://www.kaggle.com/code/stautxie/multiomer-prediction-with-svd-regression). This is likely due to the underlying implementation differences in the SVD algorithm and regression algorithm.","metadata":{}},{"cell_type":"markdown","source":"### 3.3 Model Score vs. SVD Components","metadata":{}},{"cell_type":"markdown","source":"In this section, we will analyze the impact of SVD components on model quality.\n\nIn the first experiment, we will fix the SVD component of target matrix to 16 while changing the SVD component of input matrix from 16 to 128. While the increase of fitting quality with regard to increase of SVD components gradually moderates, we still see definite increase. Therefore, it may make sense to experiment with even higher SVD components in input matrix if cpu and memory resources are still available.","metadata":{}},{"cell_type":"code","source":"cv_input_df = pd.DataFrame({'svd_comp':[16, 32, 64, 128],\n                             'mean_corr':[0.661501, 0.663247, 0.663602, 0.663845]})\n\ng = sns.relplot(data=cv_input_df,\n               x = 'svd_comp',\n               y = 'mean_corr',\n               kind = 'line').set(title=\"Mean Correlation vs. SVD Components for Training Input Matrix\")","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:27:55.680441Z","iopub.execute_input":"2022-09-12T16:27:55.680924Z","iopub.status.idle":"2022-09-12T16:27:56.011365Z","shell.execute_reply.started":"2022-09-12T16:27:55.680886Z","shell.execute_reply":"2022-09-12T16:27:56.010487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the second experiment, we will fix the SVD component of the input matrix to 16 while changing the SVD component of the target matrix from 16 to 128. While CV score increases with increases of SVD components on the target matrix, the increment is marginal after 64 components. As computation effort continually increases, it may make sense to stop at 64 components for the target matrix.","metadata":{}},{"cell_type":"code","source":"cv_target_df = pd.DataFrame({'svd_comp':[16, 32, 64, 128],\n                             'mean_corr':[0.661501, 0.662344, 0.662393, 0.662400]})\n\ng = sns.relplot(data=cv_target_df,\n               x = 'svd_comp',\n               y = 'mean_corr',\n               kind = 'line').set(title=\"Mean Correlation vs. SVD Components for Training Target Matrix\")","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:27:13.020704Z","iopub.execute_input":"2022-09-12T16:27:13.021159Z","iopub.status.idle":"2022-09-12T16:27:13.352877Z","shell.execute_reply.started":"2022-09-12T16:27:13.021124Z","shell.execute_reply":"2022-09-12T16:27:13.351398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.4 Refit Model","metadata":{}},{"cell_type":"markdown","source":"As we conclude the CV study, we will refit the model with all available data.","metadata":{}},{"cell_type":"code","source":"print(f'Refit with all data')\nmodel = MultiOutputRegressor(Ridge(alpha=1), n_jobs=5)\nmodel.fit(tr_mu_inputs_svd, tr_mu_targets_svd)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:30:09.959931Z","iopub.execute_input":"2022-09-12T16:30:09.960456Z","iopub.status.idle":"2022-09-12T16:30:12.376923Z","shell.execute_reply.started":"2022-09-12T16:30:09.960391Z","shell.execute_reply":"2022-09-12T16:30:12.375356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del tr_mu_targets_sp, tr_mu_inputs_svd, tr_mu_targets_svd\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:30:19.105200Z","iopub.execute_input":"2022-09-12T16:30:19.105689Z","iopub.status.idle":"2022-09-12T16:30:19.383126Z","shell.execute_reply.started":"2022-09-12T16:30:19.105646Z","shell.execute_reply":"2022-09-12T16:30:19.381835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Test Data Preparation & Prediction","metadata":{}},{"cell_type":"code","source":"# te_mu_inputs_sp = pickle.load(open(os.path.join(OUT_DIR, \"te_mu_inputs_sp.pkl\"), \"rb\"))\nte_mu_inputs_sp = load_npz(os.path.join(OUT_DIR, \"te_mu_inputs_sp.npz\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:30:24.888138Z","iopub.execute_input":"2022-09-12T16:30:24.888601Z","iopub.status.idle":"2022-09-12T16:31:01.112393Z","shell.execute_reply.started":"2022-09-12T16:30:24.888564Z","shell.execute_reply":"2022-09-12T16:31:01.110952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"te_mu_inputs_svd = t_svd.transform(te_mu_inputs_sp)\ndel te_mu_inputs_sp; gc.collect()\n\nte_Y_pred_svd = model.predict(te_mu_inputs_svd)\nte_Y_pred = t_svd_targets.inverse_transform(te_Y_pred_svd)\ndel te_mu_inputs_svd, te_Y_pred_svd; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:31:57.854160Z","iopub.execute_input":"2022-09-12T16:31:57.854697Z","iopub.status.idle":"2022-09-12T16:32:22.451335Z","shell.execute_reply.started":"2022-09-12T16:31:57.854654Z","shell.execute_reply":"2022-09-12T16:32:22.448866Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.save(\"te_Y_pred.npy\", te_Y_pred)\n#pickle.dump(te_Y_pred, open(\"te_Y_pred.pkl\", \"wb\"))\ndel te_Y_pred; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:35:02.569120Z","iopub.execute_input":"2022-09-12T16:35:02.569745Z","iopub.status.idle":"2022-09-12T16:35:18.016225Z","shell.execute_reply.started":"2022-09-12T16:35:02.569707Z","shell.execute_reply":"2022-09-12T16:35:18.014338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5. Prepare Submission File","metadata":{}},{"cell_type":"markdown","source":"The official submission instructions suggest that:\n\n\"Your task is to predict the labels corresponding to the inputs in the test set. To facilitate submission scoring, we only require predictions on a subset of the Multiome data. This subset was created by sampling 30% of the Multiome rows, and for each row, 15% of the columns. The sample of columns varies from row-to-row. All of the CITEseq labels are scored.\n\nevaluation_ids.csv - Identifies the labels from the test set to be evaluated. It provides a join key from the cell_id / gene_id identifiers of the label matrix to the row_id needed for the submission file.\nsample_submission.csv - A sample submission file in the correct format. See the Evaluation page for more information.\"","metadata":{}},{"cell_type":"markdown","source":"To convert matrix te_Y_pred from the last section to the required format, we observe the following:\n\n1. We can drop all the cell_ids that do not exist in the evaluation_id file.\n\n2. We only need to include all the (cell_id, gene_id) combinations that exist in the evaluation_ids.csv file.","metadata":{}},{"cell_type":"markdown","source":"Implementing #1 is fairly straightforward while implementing #2 taking quite some effort. In fact, existing notebooks, e.g., [this one](), [this one]() and [this one](), all rely on complex numpy array operations, which is quite challenging to follow. We will show how to do it intuitively within a Pandas dataframe.","metadata":{}},{"cell_type":"markdown","source":"The block of code below essentially implement idea #1 through inner join with the evaluation_ids data set.","metadata":{}},{"cell_type":"code","source":"if os.path.exists(os.path.join(SVD_DIR, \"evaluation_ids.parquet\")):\n    eval_df = pd.read_parquet(os.path.join(SVD_DIR, \"evaluation_ids.parquet\"))\nelse:\n    eval_df = pd.read_csv(os.path.join(DAT_DIR, \"evaluation_ids.csv\"))\n    eval_df.to_parquet(os.path.join(SVD_DIR, \"evaluation_ids.parquet\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:48:58.969473Z","iopub.execute_input":"2022-09-12T16:48:58.969930Z","iopub.status.idle":"2022-09-12T16:49:12.462081Z","shell.execute_reply.started":"2022-09-12T16:48:58.969893Z","shell.execute_reply":"2022-09-12T16:49:12.461100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"metadata_df = pd.read_csv(os.path.join(DAT_DIR, \"metadata.csv\"))    \n\n#~ Find only multiome related cells\n\nmulti_cell_ids = metadata_df.loc[metadata_df.technology == 'multiome', \"cell_id\"].unique()    \nmulti_eval_df = eval_df[eval_df.cell_id.isin(multi_cell_ids)]\n\ndel metadata_df; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:50:12.374756Z","iopub.execute_input":"2022-09-12T16:50:12.375239Z","iopub.status.idle":"2022-09-12T16:50:19.325899Z","shell.execute_reply.started":"2022-09-12T16:50:12.375200Z","shell.execute_reply":"2022-09-12T16:50:19.324814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#~ further reduce the size of test matrix by including only required in the test.\n\nsample_submit = pd.read_csv(os.path.join(DAT_DIR, \"sample_submission.csv\"))   \n\ntest_submit_df = pd.merge(sample_submit, multi_eval_df,\n                          on=\"row_id\",\n                          how=\"inner\")","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:50:27.138225Z","iopub.execute_input":"2022-09-12T16:50:27.138668Z","iopub.status.idle":"2022-09-12T16:51:34.660936Z","shell.execute_reply.started":"2022-09-12T16:50:27.138633Z","shell.execute_reply":"2022-09-12T16:51:34.659722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del eval_df, multi_cell_ids, multi_eval_df, sample_submit\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:51:39.335772Z","iopub.execute_input":"2022-09-12T16:51:39.337149Z","iopub.status.idle":"2022-09-12T16:51:40.672795Z","shell.execute_reply.started":"2022-09-12T16:51:39.337087Z","shell.execute_reply":"2022-09-12T16:51:40.671482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To implement #2, we make use of the fact that cell_id and gene_id columns essentially serve the roles of (row_index, column_index) in a sparse matrix format. Therefore, if we can find the mapping from cell_id to row index of te_Y_pred and gene_id to column index to te_Y_pred then we can use them to look up the predicted value within te_Y_pred matrix. More specifically, the blocks of code below implement this idea.","metadata":{}},{"cell_type":"markdown","source":"Firstly, we will retrieve row and column names of testing target matrix. Note that, while there is no file for the testing target, we can use the column names from the training target as the column names and the row nams from the testing input file as the row names.  However, these strings are stored as byte-strings. Therefore, we need to convert them into unicode-encoded string so that we can join them to the cell_id and gene_id columns from the metadata file.","metadata":{}},{"cell_type":"code","source":"\nte_mu_inputs_h5 = tb.open_file(os.path.join(DAT_DIR, \"test_multi_inputs.h5\"), \"r\")\n\nte_mu_inputs_axis0 = te_mu_inputs_h5.root.test_multi_inputs.axis0\nprint(f'te_mu_inputs_axis0.shape = {te_mu_inputs_axis0.shape}')\nprint(te_mu_inputs_axis0[:5])\n\nte_mu_inputs_axis1 = te_mu_inputs_h5.root.test_multi_inputs.axis1\nprint(f'te_mu_inputs_axis1.shape = {te_mu_inputs_axis1.shape}')\nprint(te_mu_inputs_axis1[:5])\n\n\ntr_mu_targets_h5 = tb.open_file(os.path.join(DAT_DIR, \"train_multi_targets.h5\"), \"r\")\n\ntr_mu_targets_axis0 = tr_mu_targets_h5.root.train_multi_targets.axis0\nprint(f'tr_mu_targets_axis0.shape = {tr_mu_targets_axis0.shape}')\nprint(tr_mu_targets_axis0[:5])\n\ntr_mu_targets_axis1 = tr_mu_targets_h5.root.train_multi_targets.axis1\nprint(f'tr_mu_targets_axis1.shape = {tr_mu_targets_axis1.shape}')\nprint(tr_mu_targets_axis1[:5])\n","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:51:50.415742Z","iopub.execute_input":"2022-09-12T16:51:50.416224Z","iopub.status.idle":"2022-09-12T16:51:50.479272Z","shell.execute_reply.started":"2022-09-12T16:51:50.416182Z","shell.execute_reply":"2022-09-12T16:51:50.477825Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"te_mu_targets_rows = np.array([x.decode('utf-8') for x in te_mu_inputs_axis1])\nte_mu_targets_cols = np.array([x.decode('utf-8') for x in tr_mu_targets_axis0])","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:51:56.454380Z","iopub.execute_input":"2022-09-12T16:51:56.454817Z","iopub.status.idle":"2022-09-12T16:51:56.731305Z","shell.execute_reply.started":"2022-09-12T16:51:56.454783Z","shell.execute_reply":"2022-09-12T16:51:56.729921Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will build dictionaries to map the row names to row index (from 0 to row_number - 1) and column names to column index (from 0 to column_number - 1).","metadata":{}},{"cell_type":"code","source":"te_mu_cell_dict = dict(zip(te_mu_targets_rows,\n                           range(te_mu_targets_rows.shape[0])))\n\nte_mu_gene_dict = dict(zip(te_mu_targets_cols,\n                           range(te_mu_targets_cols.shape[0])))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:52:01.354537Z","iopub.execute_input":"2022-09-12T16:52:01.355868Z","iopub.status.idle":"2022-09-12T16:52:01.421183Z","shell.execute_reply.started":"2022-09-12T16:52:01.355814Z","shell.execute_reply":"2022-09-12T16:52:01.419596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will also pre-compute the row and column indices and store them into columns so as to save subsequent search time.","metadata":{}},{"cell_type":"code","source":"test_submit_df['cell_pos'] = test_submit_df['cell_id'].map(te_mu_cell_dict)\ntest_submit_df.drop(columns=[\"cell_id\"], inplace=True)\n\ntest_submit_df['gene_pos'] = test_submit_df['gene_id'].map(te_mu_gene_dict)\ntest_submit_df.drop(columns=[\"gene_id\"], inplace=True)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:53:18.024663Z","iopub.execute_input":"2022-09-12T16:53:18.025078Z","iopub.status.idle":"2022-09-12T16:53:39.427483Z","shell.execute_reply.started":"2022-09-12T16:53:18.025045Z","shell.execute_reply":"2022-09-12T16:53:39.426305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To expedite computaton, we make use of the fact that calculation involving numpy arrays typically runs much faster than pandas. This is because \n\n1. numpy code is often run in compiled language,\n\n2. numpy code allows multiprocessing and therefore takes advantage of multiple CPUs.\n\nTherefore, we vectorize the lookup function test_pred_lookup() to make it work with the numpy arrays underlying the pandas data set test_submit_df.","metadata":{}},{"cell_type":"code","source":"# te_Y_pred = pickle.load(open(\"te_Y_pred.pkl\", \"rb\"))\nte_Y_pred = np.load(\"te_Y_pred.npy\")","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:53:49.737995Z","iopub.execute_input":"2022-09-12T16:53:49.738503Z","iopub.status.idle":"2022-09-12T16:53:56.706645Z","shell.execute_reply.started":"2022-09-12T16:53:49.738463Z","shell.execute_reply":"2022-09-12T16:53:56.705400Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def test_pred_lookup(cell_pos, gene_pos):\n    return te_Y_pred[cell_pos, gene_pos]\n\nvec_test_pred_lookup = np.vectorize(test_pred_lookup, otypes=[float])","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:54:21.195157Z","iopub.execute_input":"2022-09-12T16:54:21.195637Z","iopub.status.idle":"2022-09-12T16:54:21.201703Z","shell.execute_reply.started":"2022-09-12T16:54:21.195595Z","shell.execute_reply":"2022-09-12T16:54:21.200154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t1 = time.time()\ntest_submit_df['target'] = vec_test_pred_lookup(test_submit_df['cell_pos'].values,\n                                                test_submit_df['gene_pos'].values)\nprint(f'Processing time (sec): {time.time() - t1}')","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:54:24.291344Z","iopub.execute_input":"2022-09-12T16:54:24.292537Z","iopub.status.idle":"2022-09-12T16:55:00.591959Z","shell.execute_reply.started":"2022-09-12T16:54:24.292486Z","shell.execute_reply":"2022-09-12T16:55:00.591031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_submit_df.drop(columns=[\"cell_pos\", \"gene_pos\"], inplace=True)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:56:52.675028Z","iopub.execute_input":"2022-09-12T16:56:52.675534Z","iopub.status.idle":"2022-09-12T16:56:53.460768Z","shell.execute_reply.started":"2022-09-12T16:56:52.675492Z","shell.execute_reply":"2022-09-12T16:56:53.459332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we will save the prediction data set to be merged with the CITE prediction.","metadata":{}},{"cell_type":"code","source":"test_submit_df.to_parquet(\"te_Y_pred_mu.parquet\")","metadata":{"execution":{"iopub.status.busy":"2022-09-12T16:57:07.145743Z","iopub.execute_input":"2022-09-12T16:57:07.147000Z","iopub.status.idle":"2022-09-12T16:57:11.447090Z","shell.execute_reply.started":"2022-09-12T16:57:07.146937Z","shell.execute_reply":"2022-09-12T16:57:11.444845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 6. Conclusios & Discussions","metadata":{}},{"cell_type":"markdown","source":"In this notebook, we have present a systematic study on the prediction model for Multiomer prediction as well as \n\n1. Presenting an integration solution with SVD and regression. \n2. Analyzing the impact of SVD components on model quality.\n3. Analyzing \n3. Developing a novel test data formatting approach.","metadata":{}},{"cell_type":"markdown","source":"We can see a number of interesting research questions:\n\n1. Does it help to use more sophisticated algorithms instead of Ridge regression?\n\n2. Any model transformation techniques that could help improve the model performance?\n\n3. How to improve the computational efficiency to improve productivity?","metadata":{}}]}