{"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":"code","source":"Based on [notebook](www.kaggle.com/code/fabiencrom/msci-multiome-quickstart-w-sparse-matrices),\n\nI used TruncatedSVD(n_components=50). However, when n_components is 50 , its explained_varianceratio is only 0.008763663.\n\nTraining and test data were normalized before input.\n\nA training strategy of randomsampling. First, shuffled the train_norms and split it to n=5 parts. Trained 5 models on the 5 parts. Because of 5 fold cv, I get 25 models.","metadata":{},"execution_count":null,"outputs":[]},{"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\n\nimport gc\nimport pickle\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-07T15:09:14.95467Z","iopub.execute_input":"2022-09-07T15:09:14.955287Z","iopub.status.idle":"2022-09-07T15:09:16.366412Z","shell.execute_reply.started":"2022-09-07T15:09:14.955162Z","shell.execute_reply":"2022-09-07T15:09:16.365366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-07T15:09:16.368062Z","iopub.execute_input":"2022-09-07T15:09:16.369207Z","iopub.status.idle":"2022-09-07T15:09:16.378494Z","shell.execute_reply.started":"2022-09-07T15:09:16.369162Z","shell.execute_reply":"2022-09-07T15:09:16.376605Z"},"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-07T15:09:19.058296Z","iopub.execute_input":"2022-09-07T15:09:19.058754Z","iopub.status.idle":"2022-09-07T15:10:21.476562Z","shell.execute_reply.started":"2022-09-07T15:09:19.058702Z","shell.execute_reply":"2022-09-07T15:10:21.475259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_inputs = train_inputs.astype('float16', copy=False)","metadata":{"execution":{"iopub.status.busy":"2022-09-07T15:11:27.130595Z","iopub.execute_input":"2022-09-07T15:11:27.131095Z","iopub.status.idle":"2022-09-07T15:11:31.061605Z","shell.execute_reply.started":"2022-09-07T15:11:27.131056Z","shell.execute_reply":"2022-09-07T15:11:31.059828Z"},"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=50, random_state=42)\ntrain_inputs = pca.fit_transform(train_inputs)\nprint(pca.explained_variance_ratio_.sum())","metadata":{"execution":{"iopub.status.busy":"2022-09-07T15:12:16.97131Z","iopub.execute_input":"2022-09-07T15:12:16.971801Z","iopub.status.idle":"2022-09-07T15:24:44.044469Z","shell.execute_reply.started":"2022-09-07T15:12:16.971763Z","shell.execute_reply":"2022-09-07T15:24:44.042768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ntrain_targets = scipy.sparse.load_npz(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_targets_values.sparse.npz\")","metadata":{"execution":{"iopub.status.busy":"2022-09-07T14:40:08.562687Z","iopub.execute_input":"2022-09-07T14:40:08.56326Z","iopub.status.idle":"2022-09-07T14:40:37.484669Z","shell.execute_reply.started":"2022-09-07T14:40:08.563213Z","shell.execute_reply":"2022-09-07T14:40:37.483741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\npca2 = TruncatedSVD(n_components=50, random_state=42)\ntrain_target = pca2.fit_transform(train_targets)\nprint(pca2.explained_variance_ratio_.sum())","metadata":{"execution":{"iopub.status.busy":"2022-09-07T14:41:10.181433Z","iopub.execute_input":"2022-09-07T14:41:10.182054Z","iopub.status.idle":"2022-09-07T14:43:35.133871Z","shell.execute_reply.started":"2022-09-07T14:41:10.18199Z","shell.execute_reply":"2022-09-07T14:43:35.132538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.gaussian_process.kernels import RBF\nfrom sklearn.kernel_ridge import KernelRidge\nkernel = RBF(length_scale = 10)\nkrr = KernelRidge(alpha=0.2, kernel=kernel)","metadata":{"execution":{"iopub.status.busy":"2022-09-07T14:47:58.337708Z","iopub.execute_input":"2022-09-07T14:47:58.338493Z","iopub.status.idle":"2022-09-07T14:47:58.360671Z","shell.execute_reply.started":"2022-09-07T14:47:58.338446Z","shell.execute_reply":"2022-09-07T14:47:58.359635Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n = 5","metadata":{"execution":{"iopub.status.busy":"2022-09-07T14:56:43.381469Z","iopub.execute_input":"2022-09-07T14:56:43.381944Z","iopub.status.idle":"2022-09-07T14:56:43.388454Z","shell.execute_reply.started":"2022-09-07T14:56:43.381905Z","shell.execute_reply":"2022-09-07T14:56:43.387198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.random.seed(42)\nall_row_indices = np.arange(train_inputs.shape[0])\nnp.random.shuffle(all_row_indices)\n\nkf = KFold(n_splits=5, shuffle=True, random_state=42)\n\nindex = 0\nscore = []\n\nd = train_inputs.shape[0]//n\nfor i in range(0, n*d, d):\n    print(f'start [{i}:{i+d}]')\n    ind = all_row_indices[i:i+d]    \n    for idx_tr, idx_va in kf.split(ind):\n        X = train_inputs[ind]\n        Y = train_target[ind] #.todense()\n        Yva = train_targets[ind][idx_va]\n        Xtr, Xva = X[idx_tr], X[idx_va]\n        Ytr = Y[idx_tr]\n        del X, Y\n        gc.collect()\n        print('Train...')\n        model = krr #Ridge(copy_X=False)\n        model.fit(Xtr, Ytr)\n        del Xtr, Ytr\n        gc.collect()\n        s = correlation_score(Yva.todense(), model.predict(Xva)@pca2.components_)\n        score.append(s)\n        print(index, s)\n        del Xva, Yva\n        gc.collect()\n        pkl_filename = f\"model{index:02d}.pkl\"\n        index += 1\n        with open(pkl_filename, 'wb') as file:\n            pickle.dump(model, file)\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-07T14:56:52.297849Z","iopub.execute_input":"2022-09-07T14:56:52.299058Z","iopub.status.idle":"2022-09-07T14:57:48.818328Z","shell.execute_reply.started":"2022-09-07T14:56:52.299002Z","shell.execute_reply":"2022-09-07T14:57:48.817036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train_target, train_inputs, train_targets\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-07T15:02:21.569014Z","iopub.execute_input":"2022-09-07T15:02:21.569515Z","iopub.status.idle":"2022-09-07T15:02:21.986403Z","shell.execute_reply.started":"2022-09-07T15:02:21.569475Z","shell.execute_reply":"2022-09-07T15:02:21.985441Z"},"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)","metadata":{"execution":{"iopub.status.busy":"2022-09-07T15:03:15.463562Z","iopub.execute_input":"2022-09-07T15:03:15.464062Z","iopub.status.idle":"2022-09-07T15:04:17.758581Z","shell.execute_reply.started":"2022-09-07T15:03:15.46402Z","shell.execute_reply":"2022-09-07T15:04:17.757308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# test_sd = np.std(multi_test_x, axis=1).reshape(-1, 1)\n# test_sd[test_sd == 0] = 1\n# test_norm = (multi_test_x - np.mean(multi_test_x, axis=1).reshape(-1, 1)) / test_sd\n# test_norm = test_norm.astype(np.float16)\n# del multi_test_x\n# gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-06T02:04:05.12926Z","iopub.execute_input":"2022-09-06T02:04:05.131149Z","iopub.status.idle":"2022-09-06T02:04:05.293098Z","shell.execute_reply.started":"2022-09-06T02:04:05.131082Z","shell.execute_reply":"2022-09-06T02:04:05.291505Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_len = multi_test_x.shape[0]\nd = test_len//n\nx = []\nfor i in range(n):\n    x.append(multi_test_x[i*d:i*d+d])\ndel multi_test_x\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-07T15:04:17.761016Z","iopub.execute_input":"2022-09-07T15:04:17.761372Z","iopub.status.idle":"2022-09-07T15:04:18.141498Z","shell.execute_reply.started":"2022-09-07T15:04:17.761338Z","shell.execute_reply":"2022-09-07T15:04:18.139919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"index","metadata":{"execution":{"iopub.status.busy":"2022-09-07T15:05:37.684696Z","iopub.execute_input":"2022-09-07T15:05:37.685195Z","iopub.status.idle":"2022-09-07T15:05:37.693264Z","shell.execute_reply.started":"2022-09-07T15:05:37.685157Z","shell.execute_reply":"2022-09-07T15:05:37.692011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds = np.zeros((test_len, 23418), dtype='float16')\nfor i,xx in enumerate(x):\n    for ind in range(index):\n        print(ind, end=' ')\n        with open(f'model{ind:02}.pkl', 'rb') as file:\n            model = pickle.load(file)\n        preds[i*d:i*d+d,:] += (model.predict(xx)@pca2.components_)/index\n        gc.collect()\n    print('')\n    del xx\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-07T15:06:38.908087Z","iopub.execute_input":"2022-09-07T15:06:38.908535Z","iopub.status.idle":"2022-09-07T15:07:57.642497Z","shell.execute_reply.started":"2022-09-07T15:06:38.908497Z","shell.execute_reply":"2022-09-07T15:07:57.641605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del x\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-06T02:10:38.482932Z","iopub.execute_input":"2022-09-06T02:10:38.483316Z","iopub.status.idle":"2022-09-06T02:10:38.634478Z","shell.execute_reply.started":"2022-09-06T02:10:38.48329Z","shell.execute_reply":"2022-09-06T02:10:38.632399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.save('preds.npy', preds)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds = preds.astype('float16', copy=False)","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# 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())","metadata":{"execution":{"iopub.status.busy":"2022-09-06T03:21:01.932482Z","iopub.execute_input":"2022-09-06T03:21:01.932999Z","iopub.status.idle":"2022-09-06T03:21:35.622322Z","shell.execute_reply.started":"2022-09-06T03:21:01.932952Z","shell.execute_reply":"2022-09-06T03:21:35.621059Z"},"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-09-06T03:21:35.624004Z","iopub.execute_input":"2022-09-06T03:21:35.624421Z","iopub.status.idle":"2022-09-06T03:21:57.787288Z","shell.execute_reply.started":"2022-09-06T03:21:35.624391Z","shell.execute_reply":"2022-09-06T03:21:57.786245Z"},"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":{"iopub.status.busy":"2022-09-06T03:21:57.788566Z","iopub.execute_input":"2022-09-06T03:21:57.788879Z","iopub.status.idle":"2022-09-06T03:21:57.884422Z","shell.execute_reply.started":"2022-09-06T03:21:57.78885Z","shell.execute_reply":"2022-09-06T03:21:57.883033Z"},"trusted":true},"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-09-06T03:21:57.886943Z","iopub.execute_input":"2022-09-06T03:21:57.887407Z","iopub.status.idle":"2022-09-06T03:21:57.92178Z","shell.execute_reply.started":"2022-09-06T03:21:57.887372Z","shell.execute_reply":"2022-09-06T03:21:57.920446Z"},"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-09-06T03:21:57.924167Z","iopub.execute_input":"2022-09-06T03:21:57.925436Z","iopub.status.idle":"2022-09-06T03:22:00.076237Z","shell.execute_reply.started":"2022-09-06T03:21:57.925387Z","shell.execute_reply":"2022-09-06T03:22:00.075154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.iloc[valid_multi_rows] = preds[eval_ids_cell_num[valid_multi_rows].to_numpy(),\neval_ids_gene_num[valid_multi_rows].to_numpy()]","metadata":{"execution":{"iopub.status.busy":"2022-09-06T03:22:00.077516Z","iopub.execute_input":"2022-09-06T03:22:00.077855Z","iopub.status.idle":"2022-09-06T03:22:00.199966Z","shell.execute_reply.started":"2022-09-06T03:22:00.077825Z","shell.execute_reply":"2022-09-06T03:22:00.197716Z"},"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'","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/msci-citeseq-keras-quickstart/submission.csv\")\ncite_submission = cite_submission.set_index(\"row_id\")\ncite_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()]","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":[]}]}