{"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":"# Important note\n\n1. This notebook is built upon this quickstarter by https://www.kaggle.com/code/ambrosm/msci-multiome-quickstart\n\n2. It reads the Multime traning data in sparse mode,from this dataset: https://www.kaggle.com/datasets/sbunzini/open-problems-msci-multiome-sparse-matrices\n\n3. During CV this NB reads the whole training dataset,in four chunks. It takes 33 minutes to complete.\n\n4. It uses truncatedSVD for reducing the number of features, to 512.\n\n5. The model is Ridge. I tried LinearRegression and KNN  (very fast) and DecisionTrees, LinearSVD (dual=False),LGBM and CatBoost (excruciatingly slow).\n\n6. The CV is about 0.658, which is not an improvement. \n\nI hope this helps.","metadata":{}},{"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":{}},{"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\nimport scipy.sparse as sps\nfrom scipy.sparse.linalg import lsqr\n\nfrom sklearn.base import BaseEstimator, TransformerMixin\nfrom sklearn.model_selection import KFold\nfrom sklearn.decomposition import TruncatedSVD\n#from sklearn.preprocessing import StandardScaler, scale\n#from sklearn.decomposition import PCA\n#from sklearn.dummy import DummyRegressor\nfrom sklearn.pipeline import make_pipeline, Pipeline\nfrom sklearn.linear_model import Ridge, LinearRegression, Lasso, HuberRegressor\nfrom sklearn.metrics import mean_squared_error\n\n#import lightgbm as lgb\n#import catboost as cb\n#from catboost import CatBoost,CatBoostRegressor, Pool\n#from sklearn.multioutput import MultiOutputRegressor\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-09-02T15:51:47.086909Z","iopub.execute_input":"2022-09-02T15:51:47.087282Z","iopub.status.idle":"2022-09-02T15:51:47.097462Z","shell.execute_reply.started":"2022-09-02T15:51:47.087252Z","shell.execute_reply":"2022-09-02T15:51:47.096345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The scoring function\n\nIt is a slight modification of the original scoring function. No averages.","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   \n    return corrsum, len(y_true)\n    #return corrsum / len(y_true)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-02T15:51:51.326523Z","iopub.execute_input":"2022-09-02T15:51:51.326947Z","iopub.status.idle":"2022-09-02T15:51:51.335221Z","shell.execute_reply.started":"2022-09-02T15:51:51.326915Z","shell.execute_reply":"2022-09-02T15:51:51.334052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessing and cross-validation","metadata":{}},{"cell_type":"code","source":"%%time\n# Preprocessing\n\ndef fit_transform(X):\n        \n    print(X.shape)\n\n    tsvd = TruncatedSVD(n_components = 32,random_state = 4223)\n    X = tsvd.fit_transform(X)\n    plt.plot(tsvd.explained_variance_ratio_.cumsum())\n    plt.title(\"Cumulative explained variance ratio\")\n    plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))\n    plt.xlabel('tSVD component')\n    plt.ylabel('Cumulative explained variance ratio')\n    plt.show()\n    print(X.shape)\n    \n    return X","metadata":{"execution":{"iopub.status.busy":"2022-09-02T15:51:55.290123Z","iopub.execute_input":"2022-09-02T15:51:55.291120Z","iopub.status.idle":"2022-09-02T15:51:55.298685Z","shell.execute_reply.started":"2022-09-02T15:51:55.291078Z","shell.execute_reply":"2022-09-02T15:51:55.297693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)","metadata":{"execution":{"iopub.status.busy":"2022-09-02T15:51:58.500341Z","iopub.execute_input":"2022-09-02T15:51:58.501194Z","iopub.status.idle":"2022-09-02T15:51:58.505622Z","shell.execute_reply.started":"2022-09-02T15:51:58.501161Z","shell.execute_reply":"2022-09-02T15:51:58.504595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Apply the tSVD to the whole training set\nmulti_train_x = sps.load_npz('../input/open-problems-msci-multiome-sparse-matrices/train_multiome_input_sparse.npz')\nmulti_train_x = fit_transform(multi_train_x)","metadata":{"execution":{"iopub.status.busy":"2022-09-02T15:52:01.652317Z","iopub.execute_input":"2022-09-02T15:52:01.652719Z","iopub.status.idle":"2022-09-02T15:59:30.583967Z","shell.execute_reply.started":"2022-09-02T15:52:01.652686Z","shell.execute_reply":"2022-09-02T15:59:30.582341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Cross-validation\n\nstart = 0\nchunksize = 27000\ntotal_rows = 0\n\nfor i in range(2):\n        \n    multi_train_y = None\n    multi_train_y = sps.load_npz('../input/open-problems-msci-multiome-sparse-matrices/train_multi_targets_sparse.npz')\n    multi_train_y = multi_train_y[total_rows:chunksize + total_rows, :]\n    multi_train_y = multi_train_y.todense()\n    \n    kf = KFold(n_splits = 3, shuffle = True, random_state=4223)\n    score_list = []\n    for fold, (idx_tr, idx_va) in enumerate(kf.split(multi_train_x[total_rows:chunksize + total_rows, :])):\n        model = None\n        gc.collect()\n        X_tr = multi_train_x[idx_tr] \n        y_tr = multi_train_y[idx_tr]\n        del idx_tr\n    \n        model = Ridge(copy_X = True)\n        model.fit(X_tr, y_tr)\n   \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    \n        y_va_pred = model.predict(X_va)\n        mse = mean_squared_error(y_va, y_va_pred)\n        corrSum, nrows = correlation_score(y_va, y_va_pred)\n        del X_va, y_va\n\n        print(f\"Fold {fold}: mse = {mse:.5f}, corrSum =  {corrSum:.3f}, nrows = {nrows:.1f}, corrscore = {(corrSum/nrows):.3f}\")\n        score_list.append((mse, corrSum, nrows))\n\n    # Show overall score\n    if i == 0:\n        result_df = pd.DataFrame(score_list, columns=['mse', 'corrSum', 'nrows'])\n    else: \n        result_df = pd.concat((result_df, pd.DataFrame(score_list, columns=['mse', 'corrSum', 'nrows'])), axis = 0)\n    \n    if len(multi_train_x[total_rows:chunksize + total_rows, :]) < chunksize: break # this is the last chunk\n    total_rows += len(multi_train_x[total_rows:chunksize + total_rows, :])\n    print(total_rows)\n    start += chunksize\n    \nprint(f\"{Fore.GREEN}{Style.BRIGHT}{multi_train_x.shape} Average  mse = {result_df.mse.mean():.5f}; corr = {result_df.corrSum.sum()/result_df.nrows.sum():.3f}{Style.RESET_ALL}\")","metadata":{"execution":{"iopub.status.busy":"2022-09-02T15:59:47.335466Z","iopub.execute_input":"2022-09-02T15:59:47.335989Z","iopub.status.idle":"2022-09-02T16:01:17.625712Z","shell.execute_reply.started":"2022-09-02T15:59:47.335933Z","shell.execute_reply":"2022-09-02T16:01:17.624487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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()\n\nmulti_train_y = sps.load_npz('../input/open-problems-msci-multiome-sparse-matrices/train_multi_targets_sparse.npz')\nmulti_train_y = multi_train_y.todense()\n\nmodel = Ridge(copy_X = False) \nmodel.fit(multi_train_x[:27000,:], multi_train_y[:27000,:])","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:01:35.153278Z","iopub.execute_input":"2022-09-02T16:01:35.153753Z","iopub.status.idle":"2022-09-02T16:02:09.799978Z","shell.execute_reply.started":"2022-09-02T16:01:35.153715Z","shell.execute_reply":"2022-09-02T16:02:09.798953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del multi_train_x, multi_train_y \n_ = gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:02:40.080940Z","iopub.execute_input":"2022-09-02T16:02:40.081409Z","iopub.status.idle":"2022-09-02T16:02:40.641782Z","shell.execute_reply.started":"2022-09-02T16:02:40.081371Z","shell.execute_reply":"2022-09-02T16:02:40.640726Z"},"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":"markdown","source":"# Predicting","metadata":{}},{"cell_type":"code","source":"%%time\nmulti_test_x = sps.load_npz(\"../input/open-problems-msci-multiome-sparse-matrices/test_multi_inputs_sparse.npz\")\nmulti_test_x = fit_transform(multi_test_x)\ntest_pred = model.predict(multi_test_x)\n\ndel multi_test_x\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:06:03.406106Z","iopub.execute_input":"2022-09-02T16:06:03.406961Z","iopub.status.idle":"2022-09-02T16:10:27.545394Z","shell.execute_reply.started":"2022-09-02T16:06:03.406914Z","shell.execute_reply":"2022-09-02T16:10:27.544208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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#eval_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))\n\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-02T16:15:35.339888Z","iopub.execute_input":"2022-09-02T16:15:35.340808Z","iopub.status.idle":"2022-09-02T16:17:38.340246Z","shell.execute_reply.started":"2022-09-02T16:15:35.340770Z","shell.execute_reply":"2022-09-02T16:17:38.339053Z"},"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":"submission = pd.Series(name='target',\n                       index=pd.MultiIndex.from_frame(eval_ids), \n                       dtype=np.float32)\nsubmission","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:19:37.983741Z","iopub.execute_input":"2022-09-02T16:19:37.984209Z","iopub.status.idle":"2022-09-02T16:19:38.119544Z","shell.execute_reply.started":"2022-09-02T16:19:37.984164Z","shell.execute_reply":"2022-09-02T16:19:38.118277Z"},"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":"%%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-02T16:19:56.375882Z","iopub.execute_input":"2022-09-02T16:19:56.376546Z","iopub.status.idle":"2022-09-02T16:19:56.472432Z","shell.execute_reply.started":"2022-09-02T16:19:56.376501Z","shell.execute_reply":"2022-09-02T16:19:56.471055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-02T16:20:05.272948Z","iopub.execute_input":"2022-09-02T16:20:05.273415Z","iopub.status.idle":"2022-09-02T16:20:05.302433Z","shell.execute_reply.started":"2022-09-02T16:20:05.273378Z","shell.execute_reply":"2022-09-02T16:20:05.301183Z"},"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-02T16:20:09.427055Z","iopub.execute_input":"2022-09-02T16:20:09.428262Z","iopub.status.idle":"2022-09-02T16:20:11.712878Z","shell.execute_reply.started":"2022-09-02T16:20:09.428214Z","shell.execute_reply":"2022-09-02T16:20:11.711853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.iloc[valid_multi_rows] = test_pred[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-02T16:20:20.611430Z","iopub.execute_input":"2022-09-02T16:20:20.612191Z","iopub.status.idle":"2022-09-02T16:20:22.948574Z","shell.execute_reply.started":"2022-09-02T16:20:20.612152Z","shell.execute_reply":"2022-09-02T16:20:22.947286Z"},"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-09-02T16:20:27.072700Z","iopub.execute_input":"2022-09-02T16:20:27.073386Z","iopub.status.idle":"2022-09-02T16:20:27.213421Z","shell.execute_reply.started":"2022-09-02T16:20:27.073309Z","shell.execute_reply":"2022-09-02T16:20:27.212155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:20:30.212063Z","iopub.execute_input":"2022-09-02T16:20:30.212887Z","iopub.status.idle":"2022-09-02T16:20:30.223604Z","shell.execute_reply.started":"2022-09-02T16:20:30.212845Z","shell.execute_reply":"2022-09-02T16:20:30.222548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Merging with CITEseq predictions","metadata":{}},{"cell_type":"code","source":"submission.reset_index(drop=True, inplace=True)\nsubmission.index.name = 'row_id'","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:20:41.410948Z","iopub.execute_input":"2022-09-02T16:20:41.411613Z","iopub.status.idle":"2022-09-02T16:20:41.421312Z","shell.execute_reply.started":"2022-09-02T16:20:41.411576Z","shell.execute_reply":"2022-09-02T16:20:41.420356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_submission = pd.read_csv(\"../input/tune-lgbm-only-final-cite-task/submission.csv\")\ncite_submission = cite_submission.set_index(\"row_id\")\ncite_submission = cite_submission[\"target\"]","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:21:00.465445Z","iopub.execute_input":"2022-09-02T16:21:00.465854Z","iopub.status.idle":"2022-09-02T16:21:21.310112Z","shell.execute_reply.started":"2022-09-02T16:21:00.465819Z","shell.execute_reply":"2022-09-02T16:21:21.309002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission[submission.isnull()] = cite_submission[submission.isnull()]\nsubmission","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:21:35.048106Z","iopub.execute_input":"2022-09-02T16:21:35.048556Z","iopub.status.idle":"2022-09-02T16:21:37.350334Z","shell.execute_reply.started":"2022-09-02T16:21:35.048519Z","shell.execute_reply":"2022-09-02T16:21:37.349085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.isnull().any()","metadata":{"execution":{"iopub.status.busy":"2022-09-02T16:22:25.088506Z","iopub.execute_input":"2022-09-02T16:22:25.088934Z","iopub.status.idle":"2022-09-02T16:22:25.142866Z","shell.execute_reply.started":"2022-09-02T16:22:25.088901Z","shell.execute_reply":"2022-09-02T16:22:25.141750Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv(\"submission.csv\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!head submission.csv","metadata":{},"execution_count":null,"outputs":[]}]}