{"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":"- DIFFERENCE: tried another cv strategy  \n  \n- VER1:alpha=0.1\n- VER2:alpha=5","metadata":{}},{"cell_type":"markdown","source":"# SCP Quickstart\n\nThis notebook shows how to cross-validate a model for the *Open Problems – Single-Cell Perturbations* competition.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nimport seaborn as sns\n\nfrom sklearn.preprocessing import OneHotEncoder, OrdinalEncoder\nfrom sklearn.decomposition import TruncatedSVD\nfrom sklearn.compose import ColumnTransformer\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.linear_model import Ridge\nfrom sklearn.model_selection import StratifiedKFold\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-04T01:51:32.041408Z","iopub.execute_input":"2023-10-04T01:51:32.041772Z","iopub.status.idle":"2023-10-04T01:51:33.184242Z","shell.execute_reply.started":"2023-10-04T01:51:32.041739Z","shell.execute_reply":"2023-10-04T01:51:33.183381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading the data","metadata":{}},{"cell_type":"code","source":"id_map = pd.read_csv('/kaggle/input/open-problems-single-cell-perturbations/id_map.csv',\n                     index_col='id')\n\nde_train = pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet')\ngenes = de_train.columns[5:] # 18211 genes\nde_train","metadata":{"execution":{"iopub.status.busy":"2023-10-04T01:51:33.185461Z","iopub.execute_input":"2023-10-04T01:51:33.185798Z","iopub.status.idle":"2023-10-04T01:51:35.223449Z","shell.execute_reply.started":"2023-10-04T01:51:33.185778Z","shell.execute_reply":"2023-10-04T01:51:35.222876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cross-validation\n\nWe can see this competition as a multi-output regression task with 18211 targets and only two features (cell_type and sm_name). Both features are categorical.\n\nThe diagram shows in red which combinations of the two categorical features (cell_type, sm_name) are in the test set.\n\n**Insight:**\n- We want to predict the 18211 differential gene expressions for unseen (cell_type, sm_name) combinations. Cross-validation must simulate this setting. A possible cross-validation strategy makes four folds, of which the diagram shows the first fold in blue.\n","metadata":{}},{"cell_type":"code","source":"B_M_df = de_train[de_train[\"cell_type\"].isin([\"B cells\", \"Myeloid cells\"])].reset_index(drop=True)\nmean_score = B_M_df[[\"sm_name\"]+genes.tolist()].groupby(\"sm_name\").agg(\"mean\").mean(axis=1)\nclasses = np.digitize(mean_score.values, bins=[0.1, 0.5, 1])\ncpds = mean_score.index.values\nfold_arr = np.full(len(cpds), -1)\nskf = StratifiedKFold(n_splits=3, random_state=42, shuffle=True)\nfor fold, (trn_ind, val_ind) in enumerate(skf.split(classes, classes)):\n    fold_arr[val_ind] = fold\nfold_map = {c: f for c, f in zip(cpds, fold_arr)}\nfold_to_cpds = {fold: cpds[fold_arr==fold] for fold in range(3)}","metadata":{"execution":{"iopub.status.busy":"2023-10-04T01:51:40.707664Z","iopub.execute_input":"2023-10-04T01:51:40.707997Z","iopub.status.idle":"2023-10-04T01:51:40.747360Z","shell.execute_reply.started":"2023-10-04T01:51:40.707973Z","shell.execute_reply":"2023-10-04T01:51:40.745875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_cv_diagram(fold, fold_to_cpds):\n    cv_diagram = pd.concat([de_train[['cell_type', 'sm_name']],\n                            id_map[['cell_type', 'sm_name']]], axis=0, keys=[2, 1])\n    cv_diagram = cv_diagram.reset_index().drop(columns=['level_1'])\n    cv_diagram = cv_diagram.pivot(index='cell_type', columns='sm_name')\n    cv_diagram.fillna(0, inplace=True)\n    cv_diagram = cv_diagram.droplevel(level=0, axis=1)\n    cv_diagram = cv_diagram.sort_values('Myeloid cells', axis=1)\n    cv_diagram = cv_diagram.sort_index(ascending=False)\n#     cv_diagram.loc[val_cell_type] *= 1 + (cv_diagram.loc[['Myeloid cells', 'B cells']] == 1).any().astype(float) * 0.5\n    cpds_fold = fold_to_cpds[fold].tolist()\n    cv_diagram.loc[['Myeloid cells', 'B cells'], cpds_fold] = 3\n    # 0=Missing\n    # 1=Test\n    # 2=Training\n    # 3=Validation\n    \n    _, (ax1, ax2) = plt.subplots(1, 2, width_ratios=(20, 1), figsize=(36, 3))\n    sns.heatmap(cv_diagram, cbar=False, linewidths=1, ax=ax1, cmap=['k', 'r', 'g', 'b'])\n    \n    ax2.legend(handles=[mpatches.Patch(color='g', label='Training'),\n                        mpatches.Patch(color='b', label='Validation (fold 0)'),\n                        mpatches.Patch(color='r', label='Test'),\n                        mpatches.Patch(color='k', label='Missing')])\n    ax2.axis('off')\n    plt.show()\n    \nplot_cv_diagram(0, fold_to_cpds)\n# plot_cv_diagram('T cells CD8+')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-04T01:51:41.051515Z","iopub.execute_input":"2023-10-04T01:51:41.051852Z","iopub.status.idle":"2023-10-04T01:51:42.025986Z","shell.execute_reply.started":"2023-10-04T01:51:41.051830Z","shell.execute_reply":"2023-10-04T01:51:42.025317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our model is simple:\n- We denoise the targets by applying a singular value decomposition. As a consequence of this transformation, we must inverse-tranform the predictions.\n- We use sm_name as our unique feature, one-hot encode it and then use ridge regression to predict the transformed targets.\n- The two hyperparameters `n_components` and `alpha` were tuned for the best cross-validation score.\n","metadata":{}},{"cell_type":"code","source":"temp = de_train.groupby(['cell_type']).size() > 20\nvalidation_cell_types = temp[temp].index # 4 cell types\n\ntrain_sm_names = de_train.query(\"cell_type == 'B cells'\").sm_name.values # 17 compounds including the two control compounds\n\nfeatures = ['cell_type', 'sm_name']\n\ndef cross_val_svd(model, label, n_components=5, print_each=False):\n    mrrmse_list = []\n    for fold in range(3):\n#         mask_va = (de_train.cell_type == val_cell_type) & ~de_train.sm_name.isin(train_sm_names)\n        mask_va = de_train['cell_type'].isin(['Myeloid cells', 'B cells']) & de_train['sm_name'].isin(fold_to_cpds[fold])\n        mask_tr = ~mask_va # 485 or 487 training rows\n        \n        train = de_train[mask_tr]\n        val = de_train[mask_va]\n        y_true = val[genes]\n        \n        svd = TruncatedSVD(n_components=n_components, random_state=1)\n        z_tr = svd.fit_transform(train[genes])\n        \n        model.fit(train[features], z_tr)\n        y_pred = svd.inverse_transform(model.predict(val[features]))\n        \n        mrrmse = np.sqrt(np.square(y_true - y_pred).mean(axis=1)).mean()\n        if print_each:\n#             print(f\"# Fold {fold}: {mrrmse:5.3f} val='{val_cell_type}'\")\n            print(f\"# Fold {fold}: {mrrmse:5.3f}\")\n        mrrmse_list.append(mrrmse)\n    mrrmse = np.array(mrrmse_list).mean()\n    print(f\"# Overall {mrrmse:5.3f} {label}\")\n    return mrrmse","metadata":{"execution":{"iopub.status.busy":"2023-10-04T01:58:31.489460Z","iopub.execute_input":"2023-10-04T01:58:31.489815Z","iopub.status.idle":"2023-10-04T01:58:31.707049Z","shell.execute_reply.started":"2023-10-04T01:58:31.489785Z","shell.execute_reply":"2023-10-04T01:58:31.706056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# alpha_list = [0.1, 1, 5, 10]\n# n_comps_list = np.arange(10, 160, 10)\n# results = {alpha: [] for alpha in alpha_list}\n# for n_comps in n_comps_list:\n#     for alpha in alpha_list:\n#         model = make_pipeline(ColumnTransformer([('ohe', OneHotEncoder(), ['sm_name'])]),\n#                       Ridge(alpha=alpha, fit_intercept=False))\n#         mrrmse = cross_val_svd(model, f'svd {n_comps=} ridge {alpha=}',\n#                                n_components=n_comps)\n#         results[alpha].append(mrrmse)\n\n# fig, ax = plt.subplots(figsize=(8,4))\n# for alpha in [0.1, 1, 5, 10]:\n#     ax.plot(n_comps_list, results[alpha], \"-o\", label=f\"Alpha {alpha}\")\n# ax.set_xlabel(\"n_components\")\n# ax.set_ylabel(\"MRRMSE\")\n# ax.legend()\n# fig.tight_layout()\n# plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-04T01:51:50.621242Z","iopub.execute_input":"2023-10-04T01:51:50.621517Z","iopub.status.idle":"2023-10-04T01:56:37.720465Z","shell.execute_reply.started":"2023-10-04T01:51:50.621489Z","shell.execute_reply":"2023-10-04T01:56:37.719483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_components = 30\nalpha = 5\nmodel = make_pipeline(ColumnTransformer([('ohe', OneHotEncoder(), ['sm_name'])]),\n                      Ridge(alpha=alpha, fit_intercept=False))\ncross_val_svd(model, f'svd {n_components=} ridge {alpha=}',\n              n_components=n_components, print_each=True)\n# Fold 0: 1.089 val='NK cells'\n# Fold 1: 0.971 val='T cells CD4+'\n# Fold 2: 0.791 val='T cells CD8+'\n# Fold 3: 1.008 val='T regulatory cells'\n# Overall 0.965 svd n_components=100 ridge alpha=5","metadata":{"execution":{"iopub.status.busy":"2023-10-04T01:58:34.115288Z","iopub.execute_input":"2023-10-04T01:58:34.115573Z","iopub.status.idle":"2023-10-04T01:58:36.843063Z","shell.execute_reply.started":"2023-10-04T01:58:34.115553Z","shell.execute_reply":"2023-10-04T01:58:36.842434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission\n\nWe retrain the model on the full training data and create a submission file.","metadata":{}},{"cell_type":"code","source":"svd = TruncatedSVD(n_components=n_components, random_state=1)\nz_tr = svd.fit_transform(de_train[genes])\nmodel.fit(de_train[features], z_tr)\ny_pred = svd.inverse_transform(model.predict(id_map[features]))\n","metadata":{"execution":{"iopub.status.busy":"2023-10-04T01:59:14.572114Z","iopub.execute_input":"2023-10-04T01:59:14.572401Z","iopub.status.idle":"2023-10-04T01:59:15.521364Z","shell.execute_reply.started":"2023-10-04T01:59:14.572379Z","shell.execute_reply":"2023-10-04T01:59:15.520643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.DataFrame(y_pred, columns=genes, index=id_map.index)\ndisplay(submission)\nsubmission.to_csv('submission.csv')\n","metadata":{"execution":{"iopub.status.busy":"2023-10-04T01:59:17.013843Z","iopub.execute_input":"2023-10-04T01:59:17.014187Z","iopub.status.idle":"2023-10-04T01:59:21.988628Z","shell.execute_reply.started":"2023-10-04T01:59:17.014164Z","shell.execute_reply":"2023-10-04T01:59:21.987638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}