{"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 Valuation Decomposition for CITE prediction","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"## 1. Introduction","metadata":{}},{"cell_type":"markdown","source":"In this notebook, we will develop a SVD-based approach to predict CITE problem. We have applied this approach on Multiomer prediction previously in [this notebook](https://www.kaggle.com/code/stautxie/svd-decomposition-for-multiomer-prediction) and achieved very good results.\n\nThis structure of this notebook is as follows.\n\n1. We presented Singular Value Decomposition in Section 2.\n2. We present a 5-fold cross-validation (CV) study in Section 3 and compare the solution quality of Ridge and Lasso regression.\n3. We apply the model to the testing set to get prediction.\n4. We merge the prediction on CITE with our previous prediction on Multiomer from [this notebook](https://www.kaggle.com/code/stautxie/svd-decomposition-for-multiomer-prediction) to submit the full prediction.","metadata":{}},{"cell_type":"markdown","source":"## 2. Singular Value 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","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:05:24.269691Z","iopub.execute_input":"2022-09-13T02:05:24.270100Z","iopub.status.idle":"2022-09-13T02:05:24.277570Z","shell.execute_reply.started":"2022-09-13T02:05:24.270065Z","shell.execute_reply":"2022-09-13T02:05:24.276351Z"},"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-13T02:05:27.098441Z","iopub.execute_input":"2022-09-13T02:05:27.099369Z","iopub.status.idle":"2022-09-13T02:05:40.663303Z","shell.execute_reply.started":"2022-09-13T02:05:27.099330Z","shell.execute_reply":"2022-09-13T02:05:40.661955Z"},"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 = '.'\nMULT_DIR = \"../input/output-from-svd-for-multiomer-prediction\"\nseed = 42","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:43:42.869341Z","iopub.execute_input":"2022-09-13T02:43:42.869837Z","iopub.status.idle":"2022-09-13T02:43:42.876248Z","shell.execute_reply.started":"2022-09-13T02:43:42.869797Z","shell.execute_reply":"2022-09-13T02:43:42.875021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_ci_inputs_sp = load_npz(os.path.join(OUT_DIR, \"tr_ci_inputs_sp.npz\"))\ntr_ci_targets_sp = load_npz(os.path.join(OUT_DIR, \"tr_ci_targets_sp.npz\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:06:46.272452Z","iopub.execute_input":"2022-09-13T02:06:46.272860Z","iopub.status.idle":"2022-09-13T02:07:07.159670Z","shell.execute_reply.started":"2022-09-13T02:06:46.272828Z","shell.execute_reply":"2022-09-13T02:07:07.158580Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here, we first apply SVD on the input matrix.","metadata":{}},{"cell_type":"code","source":"\nsvd_comp = 128\n\nt1 = time.time()\nif os.path.exists(os.path.join(SVD_DIR, f\"tr_ci_inputs_svd_{str(svd_comp)}.pkl\")):\n    t_svd = pickle.load(open(os.path.join(SVD_DIR, f\"tr_ci_inputs_svd_{str(svd_comp)}.pkl\"), \"rb\"))\nelse:\n    t_svd = TruncatedSVD(n_components = svd_comp, random_state = seed)\n    t_svd.fit(tr_ci_inputs_sp)\n    pickle.dump(t_svd, open(os.path.join(SVD_DIR, f\"tr_ci_inputs_svd_{str(svd_comp)}.pkl\"), \"wb\"))\n    \nprint(f'CPU time used : {time.time() - t1}')\ntr_ci_inputs_svd = t_svd.transform(tr_ci_inputs_sp)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:21:06.263433Z","iopub.execute_input":"2022-09-13T02:21:06.263873Z","iopub.status.idle":"2022-09-13T02:25:30.235811Z","shell.execute_reply.started":"2022-09-13T02:21:06.263838Z","shell.execute_reply":"2022-09-13T02:25:30.234793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We also apply SVD on the target matrix.","metadata":{}},{"cell_type":"code","source":"\nsvd_comp_targets = 128\n\nt1 = time.time()\nif os.path.exists(os.path.join(SVD_DIR, f\"tr_ci_targets_svd_{str(svd_comp_targets)}.pkl\")):\n    t_svd_targets = pickle.load(open(os.path.join(SVD_DIR, \n                                                  f\"tr_ci_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_ci_targets_sp)\n    pickle.dump(t_svd_targets, open(os.path.join(SVD_DIR, \n                                                  f\"tr_ci_targets_svd_{str(svd_comp_targets)}.pkl\"), \"wb\"))\n    \nprint(f'CPU time used : {time.time() - t1}')    \ntr_ci_targets_svd = t_svd_targets.transform(tr_ci_targets_sp)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:27:36.227702Z","iopub.execute_input":"2022-09-13T02:27:36.228520Z","iopub.status.idle":"2022-09-13T02:27:44.953902Z","shell.execute_reply.started":"2022-09-13T02:27:36.228472Z","shell.execute_reply":"2022-09-13T02:27:44.952710Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Cross-validation","metadata":{}},{"cell_type":"markdown","source":"In this section, we will conduct cross-validation. We will also compare Ridge and Lasso regressions.","metadata":{}},{"cell_type":"markdown","source":"### 3.1 Ridge Regression","metadata":{}},{"cell_type":"code","source":"\ndef 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))","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:35:47.492111Z","iopub.execute_input":"2022-09-13T02:35:47.493470Z","iopub.status.idle":"2022-09-13T02:35:47.499276Z","shell.execute_reply.started":"2022-09-13T02:35:47.493407Z","shell.execute_reply":"2022-09-13T02:35:47.498408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nnum_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_ci_inputs_svd)):\n    tr_X = tr_ci_inputs_svd[tr_index, :]\n    va_X = tr_ci_inputs_svd[va_index, :]\n    tr_Y = tr_ci_targets_svd[tr_index, :]\n    va_Y = tr_ci_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_ci_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    cv_list.append((va_pred, mse_val, mean_corr_val))\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-13T02:35:50.287517Z","iopub.execute_input":"2022-09-13T02:35:50.288169Z","iopub.status.idle":"2022-09-13T02:36:13.640053Z","shell.execute_reply.started":"2022-09-13T02:35:50.288130Z","shell.execute_reply":"2022-09-13T02:36:13.638700Z"},"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_ci_inputs_svd)):\n    tr_X = tr_ci_inputs_svd[tr_index, :]\n    va_X = tr_ci_inputs_svd[va_index, :]\n    tr_Y = tr_ci_targets_svd[tr_index, :]\n    va_Y = tr_ci_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_ci_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    cv_list.append((va_pred, mse_val, mean_corr_val))\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-13T02:41:39.573339Z","iopub.execute_input":"2022-09-13T02:41:39.574965Z","iopub.status.idle":"2022-09-13T02:42:17.989575Z","shell.execute_reply.started":"2022-09-13T02:41:39.574886Z","shell.execute_reply":"2022-09-13T02:42:17.987872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.3 Refit Model","metadata":{}},{"cell_type":"markdown","source":"As we can see, Ridge regression performs better than Lasso. Therefore, 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_ci_inputs_svd, tr_ci_targets_svd)\n\ndel tr_ci_inputs_sp, tr_ci_targets_sp, tr_ci_inputs_svd, tr_ci_targets_svd\ngc.collect()\n","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:43:14.938325Z","iopub.execute_input":"2022-09-13T02:43:14.938821Z","iopub.status.idle":"2022-09-13T02:43:19.480583Z","shell.execute_reply.started":"2022-09-13T02:43:14.938784Z","shell.execute_reply":"2022-09-13T02:43:19.479425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Predicting on the Test Data","metadata":{}},{"cell_type":"markdown","source":"In this section, we will first use the fitted model to predict the test CITE data. Once the prediction is available, we will combine it with our prediction on Multiomer to obtain the final submission file.","metadata":{}},{"cell_type":"markdown","source":"### 4.1 Predicting on CITE testing data","metadata":{}},{"cell_type":"code","source":"te_ci_inputs_sp = load_npz(os.path.join(OUT_DIR, \"te_ci_inputs_sp.npz\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:43:51.705693Z","iopub.execute_input":"2022-09-13T02:43:51.706153Z","iopub.status.idle":"2022-09-13T02:44:06.421886Z","shell.execute_reply.started":"2022-09-13T02:43:51.706116Z","shell.execute_reply":"2022-09-13T02:44:06.420613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nte_ci_inputs_svd = t_svd.transform(te_ci_inputs_sp)\ndel te_ci_inputs_sp; gc.collect()\n\nte_Y_pred_svd = model.predict(te_ci_inputs_svd)\nte_Y_pred = t_svd_targets.inverse_transform(te_Y_pred_svd)\ndel te_ci_inputs_svd, te_Y_pred_svd; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:44:35.443528Z","iopub.execute_input":"2022-09-13T02:44:35.444394Z","iopub.status.idle":"2022-09-13T02:44:47.170665Z","shell.execute_reply.started":"2022-09-13T02:44:35.444348Z","shell.execute_reply":"2022-09-13T02:44:47.169481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will reformat the prediction matrix into the submission file's required format.","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\"))\n\nmetadata_df = pd.read_csv(os.path.join(DAT_DIR, \"metadata.csv\"))   ","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:45:44.459256Z","iopub.execute_input":"2022-09-13T02:45:44.459697Z","iopub.status.idle":"2022-09-13T02:46:54.821208Z","shell.execute_reply.started":"2022-09-13T02:45:44.459662Z","shell.execute_reply":"2022-09-13T02:46:54.819828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_cell_ids = metadata_df.loc[metadata_df.technology == 'citeseq', \"cell_id\"].unique()   \ndel metadata_df; gc.collect() \n \ncite_eval_df = eval_df[eval_df.cell_id.isin(cite_cell_ids)]","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:47:21.861246Z","iopub.execute_input":"2022-09-13T02:47:21.861720Z","iopub.status.idle":"2022-09-13T02:47:24.821194Z","shell.execute_reply.started":"2022-09-13T02:47:21.861681Z","shell.execute_reply":"2022-09-13T02:47:24.819901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submit = pd.read_csv(os.path.join(DAT_DIR, \"sample_submission.csv\"))   \n\ntest_submit_df = pd.merge(sample_submit, cite_eval_df,\n                          on=\"row_id\",\n                          how=\"inner\")","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:47:28.822290Z","iopub.execute_input":"2022-09-13T02:47:28.822722Z","iopub.status.idle":"2022-09-13T02:48:12.837670Z","shell.execute_reply.started":"2022-09-13T02:47:28.822688Z","shell.execute_reply":"2022-09-13T02:48:12.836519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del eval_df, cite_cell_ids, cite_eval_df, sample_submit\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:48:24.996669Z","iopub.execute_input":"2022-09-13T02:48:24.997095Z","iopub.status.idle":"2022-09-13T02:48:25.926790Z","shell.execute_reply.started":"2022-09-13T02:48:24.997061Z","shell.execute_reply":"2022-09-13T02:48:25.925877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To join the prediction in the te_Y_matrix to the submission file, we will need to build lookup function from cell_id and gene_id to respective row or column position in the te_Y_matrix. Detailed discussion can be found in [this notebook](https://www.kaggle.com/code/stautxie/svd-decomposition-for-multiomer-prediction).","metadata":{}},{"cell_type":"code","source":"te_ci_inputs_h5 = tb.open_file(os.path.join(DAT_DIR, \"test_cite_inputs.h5\"), \"r\")\n\nte_ci_inputs_axis0 = te_ci_inputs_h5.root.test_cite_inputs.axis0\nprint(f'te_ci_inputs_axis0.shape = {te_ci_inputs_axis0.shape}')\nprint(te_ci_inputs_axis0[:5])\n\nte_ci_inputs_axis1 = te_ci_inputs_h5.root.test_cite_inputs.axis1\nprint(f'te_ci_inputs_axis1.shape = {te_ci_inputs_axis1.shape}')\nprint(te_ci_inputs_axis1[:5])\n\n\ntr_ci_targets_h5 = tb.open_file(os.path.join(DAT_DIR, \"train_cite_targets.h5\"), \"r\")\n\ntr_ci_targets_axis0 = tr_ci_targets_h5.root.train_cite_targets.axis0\nprint(f'tr_ci_targets_axis0.shape = {tr_ci_targets_axis0.shape}')\nprint(tr_ci_targets_axis0[:5])\n\ntr_ci_targets_axis1 = tr_ci_targets_h5.root.train_cite_targets.axis1\nprint(f'tr_ci_targets_axis1.shape = {tr_ci_targets_axis1.shape}')\nprint(tr_ci_targets_axis1[:5])\n\nte_ci_targets_rows = np.array([x.decode('utf-8') for x in te_ci_inputs_axis1])\nte_ci_targets_cols = np.array([x.decode('utf-8') for x in tr_ci_targets_axis0])","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:48:32.759265Z","iopub.execute_input":"2022-09-13T02:48:32.759735Z","iopub.status.idle":"2022-09-13T02:48:32.984978Z","shell.execute_reply.started":"2022-09-13T02:48:32.759698Z","shell.execute_reply":"2022-09-13T02:48:32.983599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nte_ci_cell_dict = dict(zip(te_ci_targets_rows,\n                           range(te_ci_targets_rows.shape[0])))\n\nte_ci_gene_dict = dict(zip(te_ci_targets_cols,\n                           range(te_ci_targets_cols.shape[0])))\n\n\ntest_submit_df['cell_pos'] = test_submit_df['cell_id'].map(te_ci_cell_dict)\ntest_submit_df['gene_pos'] = test_submit_df['gene_id'].map(te_ci_gene_dict)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:48:40.492218Z","iopub.execute_input":"2022-09-13T02:48:40.492896Z","iopub.status.idle":"2022-09-13T02:48:41.778883Z","shell.execute_reply.started":"2022-09-13T02:48:40.492841Z","shell.execute_reply":"2022-09-13T02:48:41.777843Z"},"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])\n\ntest_submit_df['target'] = vec_test_pred_lookup(test_submit_df['cell_pos'].values,\n                                                test_submit_df['gene_pos'].values)\n\ndel te_Y_pred; gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:48:46.585406Z","iopub.execute_input":"2022-09-13T02:48:46.585883Z","iopub.status.idle":"2022-09-13T02:48:50.238105Z","shell.execute_reply.started":"2022-09-13T02:48:46.585847Z","shell.execute_reply":"2022-09-13T02:48:50.236880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we load the saved prediction from [our previous work on Multiomer prediction](https://www.kaggle.com/code/stautxie/svd-decomposition-for-multiomer-prediction) to merge.","metadata":{}},{"cell_type":"code","source":"#~ load prediction of multiomer\n\nmu_test_submit_df = pd.read_parquet(os.path.join(MULT_DIR, \"te_Y_pred_mu.parquet\"))\n\nfinal_submit_df = pd.concat([test_submit_df[['row_id', 'target']],\n                             mu_test_submit_df],\n                             axis=0)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:49:30.896777Z","iopub.execute_input":"2022-09-13T02:49:30.897576Z","iopub.status.idle":"2022-09-13T02:49:39.888996Z","shell.execute_reply.started":"2022-09-13T02:49:30.897522Z","shell.execute_reply":"2022-09-13T02:49:39.887734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_submit_df.to_csv(\"submission.csv\", index=None)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T02:50:40.839365Z","iopub.execute_input":"2022-09-13T02:50:40.839900Z"},"trusted":true},"execution_count":null,"outputs":[]}]}