{"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":"# Note: \n\nThis NB is taken from AmbrosM QuickStart CITEseg model: https://www.kaggle.com/code/ambrosm/msci-citeseq-quickstart\n\nIt reads the whole training CITEseq file and runs PCA with 512 components.\n\nIt implements the LGBM from: https://www.kaggle.com/code/vuonglam/lgbm-baseline-optuna-drop-constant-cite-task/notebook\n\nThe predictions for the Multiome come from: https://www.kaggle.com/ambrosm/msci-multiome-quickstart\n\nA cool and useful EDA from AmbrosM here: https://www.kaggle.com/ambrosm/msci-eda-which-makes-sense","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import os, gc, pickle\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\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\nfrom sklearn.dummy import DummyRegressor\nfrom sklearn.pipeline import make_pipeline, Pipeline\n#from sklearn.linear_model import Ridge\nfrom sklearn.compose import TransformedTargetRegressor\nfrom sklearn.metrics import mean_squared_error\n\nfrom sklearn.multioutput import MultiOutputRegressor\nimport lightgbm as lgb\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-08-30T03:09:04.792449Z","iopub.execute_input":"2022-08-30T03:09:04.792810Z","iopub.status.idle":"2022-08-30T03:09:06.296559Z","shell.execute_reply.started":"2022-08-30T03:09:04.792738Z","shell.execute_reply":"2022-08-30T03:09:06.295270Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# If you see a warning \"Failed to establish a new connection\" running this cell,\n# go to \"Settings\" on the right hand side, and turn on internet. \n# We need this library to read HDF files.\n!pip install --quiet tables\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-30T03:09:13.156638Z","iopub.execute_input":"2022-08-30T03:09:13.156919Z","iopub.status.idle":"2022-08-30T03:09:23.624743Z","shell.execute_reply.started":"2022-08-30T03:09:13.156896Z","shell.execute_reply":"2022-08-30T03:09:23.623888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loading the common metadata table\n\nThe current version of the model is so primitive that it doesn't use the metadata, but we load it anyway.","metadata":{}},{"cell_type":"code","source":"#df_cell = pd.read_csv(FP_CELL_METADATA)\n#df_cell_cite = df_cell[df_cell.technology==\"citeseq\"]\n#df_cell_multi = df_cell[df_cell.technology==\"multiome\"]\n#df_cell_cite.shape, df_cell_multi.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-27T04:44:47.017709Z","iopub.execute_input":"2022-08-27T04:44:47.018155Z","iopub.status.idle":"2022-08-27T04:44:47.479599Z","shell.execute_reply.started":"2022-08-27T04:44:47.018118Z","shell.execute_reply":"2022-08-27T04:44:47.478502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The scoring function\n\nThis competition has a special metric: For every row, it computes the Pearson correlation between y_true and y_pred, and then all these correlation coefficients are averaged.","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    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-08-30T03:09:27.259820Z","iopub.execute_input":"2022-08-30T03:09:27.260132Z","iopub.status.idle":"2022-08-30T03:09:27.265888Z","shell.execute_reply.started":"2022-08-30T03:09:27.260107Z","shell.execute_reply":"2022-08-30T03:09:27.265037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cross-validation\n\nData size:\n- The training input has shape 70988\\*22050 (10.6 GByte).\n- The training labels have shape 70988\\*140.\n- The test input has shape 48663\\*22050 (4.3 GByte).\n\nTo get a result with only 16 GByte RAM, we simplify the problem as follows:\n- We ignore the complete metadata (donors, days, cell types).\n- We drop all feature columns which are constant.\n- We do a PCA and keep only the 512 most important components.\n- We use PCA(copy=False), which overwrites its input in fit_transform().\n- We fit a LGBM regression model with 70988\\*240 inputs and 70988\\*140 outputs. ","metadata":{}},{"cell_type":"code","source":"%%time\nclass PreprocessCiteseq(BaseEstimator, TransformerMixin):\n    def transform(self, X):\n        print(X.shape)\n        gc.collect()\n        X = self.pca.transform(X)\n        print(X.shape)\n        return X\n\n    def fit_transform(self, X):\n        print(X.shape)\n        gc.collect()\n        self.pca = PCA(n_components=512, copy=False, random_state=42)\n        X = self.pca.fit_transform(X)\n        print(X.shape)\n        return X","metadata":{"execution":{"iopub.status.busy":"2022-08-30T03:10:05.739031Z","iopub.execute_input":"2022-08-30T03:10:05.739556Z","iopub.status.idle":"2022-08-30T03:10:05.745088Z","shell.execute_reply.started":"2022-08-30T03:10:05.739532Z","shell.execute_reply":"2022-08-30T03:10:05.744344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preprocessor = PreprocessCiteseq()\ncite_train_x = None\ncite_train_x = preprocessor.fit_transform(pd.read_hdf(FP_CITE_TRAIN_INPUTS).values)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T03:10:09.683733Z","iopub.execute_input":"2022-08-30T03:10:09.684053Z","iopub.status.idle":"2022-08-30T03:13:44.941735Z","shell.execute_reply.started":"2022-08-30T03:10:09.684029Z","shell.execute_reply":"2022-08-30T03:13:44.940233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_train_y = pd.read_hdf(FP_CITE_TRAIN_TARGETS).values\nprint(cite_train_y.shape)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T03:16:13.836564Z","iopub.execute_input":"2022-08-30T03:16:13.836916Z","iopub.status.idle":"2022-08-30T03:16:14.478091Z","shell.execute_reply.started":"2022-08-30T03:16:13.836890Z","shell.execute_reply":"2022-08-30T03:16:14.477043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"params = {\n     'learning_rate': 0.1,\n     'objective' : 'regression',\n     'metric': 'rmse',#mae', \n     'random_state': 4223,\n     'reg_alpha': 0.03, \n     'reg_lambda': 0.002, \n     'colsample_bytree': 0.8, \n     'subsample': 0.6, \n     'max_depth': 10, \n     'num_leaves': 186, \n     'min_child_samples': 263\n    }","metadata":{"execution":{"iopub.status.busy":"2022-08-30T03:16:18.897837Z","iopub.execute_input":"2022-08-30T03:16:18.898257Z","iopub.status.idle":"2022-08-30T03:16:18.903640Z","shell.execute_reply.started":"2022-08-30T03:16:18.898148Z","shell.execute_reply":"2022-08-30T03:16:18.902715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Cross-validation\n\nhaz_CV = False\n\nif haz_CV:\n    small_cite_train_x, small_cite_train_y =  cite_train_x[:1000,:],cite_train_y[:1000,:]\n    kf = KFold(n_splits = 2, shuffle=True, random_state=1)\n    score_list = []\n    for fold, (idx_tr, idx_va) in enumerate(kf.split(small_cite_train_x)):\n        model = None\n        gc.collect()\n        X_tr = small_cite_train_x[idx_tr] \n        y_tr = small_cite_train_y[idx_tr]\n\n        model = MultiOutputRegressor(lgb.LGBMRegressor(**params))\n        model.fit(X_tr, y_tr)\n        del X_tr, y_tr\n        gc.collect()\n\n        # We validate the model\n        X_va = small_cite_train_x[idx_va]\n        y_va = small_cite_train_y[idx_va]\n    \n        y_va_pred = model.predict(X_va)\n        mse = mean_squared_error(y_va, y_va_pred)\n        corrscore = correlation_score(y_va, y_va_pred)\n        del X_va, y_va\n\n        print(f\"Fold {fold}: mse = {mse:.5f}, corr =  {corrscore:.3f}\")\n        score_list.append((mse, corrscore))\n\n    # Show overall score\n    result_df = pd.DataFrame(score_list, columns=['mse', 'corrscore'])\n    print(f\"{Fore.GREEN}{Style.BRIGHT}Average  mse = {result_df.mse.mean():.5f}; corr = {result_df.corrscore.mean():.3f}{Style.RESET_ALL}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-30T03:33:17.717451Z","iopub.execute_input":"2022-08-30T03:33:17.717813Z","iopub.status.idle":"2022-08-30T03:56:56.166605Z","shell.execute_reply.started":"2022-08-30T03:33:17.717787Z","shell.execute_reply":"2022-08-30T03:56:56.165551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This ridge regression is far better than DummyRegressor, which has the score `mse = 4.87618; corr = 0.805`.","metadata":{}},{"cell_type":"markdown","source":"# Retraining\n\nWe retrain the model on all training rows, delete the training data, load the test data and compute the predictions.","metadata":{}},{"cell_type":"code","source":"model = None # free the RAM occupied by the old model\nmodel = MultiOutputRegressor(lgb.LGBMRegressor(**params, n_estimators = 300))\n\nmodel.fit(cite_train_x, cite_train_y)\ndel cite_train_x, cite_train_y\ngc.collect()\n\ncite_test_x = preprocessor.transform(pd.read_hdf(FP_CITE_TEST_INPUTS).values)\ntest_pred = model.predict(cite_test_x)\ndel cite_test_x\ntest_pred.shape\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-08-27T03:25:07.817061Z","iopub.execute_input":"2022-08-27T03:25:07.81816Z","iopub.status.idle":"2022-08-27T03:26:00.975027Z","shell.execute_reply.started":"2022-08-27T03:25:07.818114Z","shell.execute_reply":"2022-08-27T03:26:00.973425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission\n\nWe save the CITEseq predictions so that they can be merged with the Multiome predictions in the [Multiome quickstart notebook](https://www.kaggle.com/ambrosm/msci-multiome-quickstart).\n\nThe CITEseq test predictions produced by the ridge regressor have 48663 rows (i.e., cells) and 140 columns (i.e. proteins). 48663 * 140 = 6812820.\n","metadata":{}},{"cell_type":"code","source":"with open('citeseq_pred.pickle', 'wb') as f: pickle.dump(test_pred, f) # float32 array of shape (48663, 140)","metadata":{"execution":{"iopub.status.busy":"2022-08-21T05:22:55.629053Z","iopub.execute_input":"2022-08-21T05:22:55.629877Z","iopub.status.idle":"2022-08-21T05:22:55.693746Z","shell.execute_reply.started":"2022-08-21T05:22:55.629823Z","shell.execute_reply":"2022-08-21T05:22:55.692043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The final submission will have 65744180 rows, of which the first 6812820 are for the CITEseq predictions and the remaining 58931360 for the Multiome predictions. \n\nWe now read the Multiome predictions and merge the CITEseq predictions into them:","metadata":{}},{"cell_type":"code","source":"with open(\"../input/msci-multiome-quickstart/partial_submission_multi.pickle\", 'rb') as f: submission = pickle.load(f)\nsubmission.iloc[:len(test_pred.ravel())] = test_pred.ravel()\nassert not submission.isna().any()\nsubmission = submission.round(6) # reduce the size of the csv\nsubmission.to_csv('submission.csv')\nsubmission\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T05:22:55.695798Z","iopub.execute_input":"2022-08-21T05:22:55.69667Z","iopub.status.idle":"2022-08-21T05:25:00.380518Z","shell.execute_reply.started":"2022-08-21T05:22:55.696623Z","shell.execute_reply":"2022-08-21T05:25:00.379105Z"},"trusted":true},"execution_count":null,"outputs":[]}]}