{"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":"# Multiome LGBM Baseline\n\n* This notebook will be implemented in the LGBM model using the data processed in the quick start.\n* LGBM models usually cannot output multiple target variables, but this method can output\n* The reference notes for data processing are below. https://www.kaggle.com/code/ambrosm/msci-multiome-quickstart","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 numpy as np\nfrom colorama import Fore, Back, Style\nfrom matplotlib.ticker import MaxNLocator\nimport warnings\nwarnings.simplefilter('ignore')\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\nfrom sklearn.linear_model import Ridge, LinearRegression\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.multioutput import MultiOutputRegressor\nimport lightgbm as lgb\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-25T01:41:02.85346Z","iopub.execute_input":"2022-08-25T01:41:02.854131Z","iopub.status.idle":"2022-08-25T01:41:04.591197Z","shell.execute_reply.started":"2022-08-25T01:41:02.853986Z","shell.execute_reply":"2022-08-25T01:41:04.590218Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install --quiet tables","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-25T01:41:04.593056Z","iopub.execute_input":"2022-08-25T01:41:04.593636Z","iopub.status.idle":"2022-08-25T01:41:18.830442Z","shell.execute_reply.started":"2022-08-25T01:41:04.593598Z","shell.execute_reply":"2022-08-25T01:41:18.828967Z"},"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)\ndf_cell_cite = df_cell[df_cell.technology==\"citeseq\"]\ndf_cell_multi = df_cell[df_cell.technology==\"multiome\"]\ndf_cell_cite.shape, df_cell_multi.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-25T01:41:18.833261Z","iopub.execute_input":"2022-08-25T01:41:18.833791Z","iopub.status.idle":"2022-08-25T01:41:19.353376Z","shell.execute_reply.started":"2022-08-25T01:41:18.833733Z","shell.execute_reply":"2022-08-25T01:41:19.352108Z"},"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":"markdown","source":"# Preprocessing and cross-validation\n\nThe Multiome dataset is way too large to fit into 16 GByte RAM:\n- train inputs:  105942 * 228942 float32 values (97 GByte)\n- train targets: 105942 *  23418 float32 values (10 GByte)\n- test inputs:    55935 * 228942 float32 values (13 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 read only 6000 rows of the training data.\n- We drop all feature columns which are constant.","metadata":{}},{"cell_type":"code","source":"#%%time\n# Preprocessing\n\nclass PreprocessMultiome(BaseEstimator, TransformerMixin):\n    columns_to_use = slice(10000, 14000)\n    \n    @staticmethod\n    def take_column_subset(X):\n        return X[:,PreprocessMultiome.columns_to_use]\n    \n    def transform(self, X):\n        print(X.shape)\n        X = X[:,~self.all_zero_columns]\n        print(X.shape)\n        X = PreprocessMultiome.take_column_subset(X) # use only a part of the columns\n        print(X.shape)\n        gc.collect()\n\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        self.all_zero_columns = (X == 0).all(axis=0)\n        X = X[:,~self.all_zero_columns]\n        print(X.shape)\n        X = PreprocessMultiome.take_column_subset(X) # use only a part of the columns\n        print(X.shape)\n        gc.collect()\n\n        self.pca = PCA(n_components=4, copy=False, random_state=1)\n        X = self.pca.fit_transform(X)\n#         plt.plot(self.pca.explained_variance_ratio_.cumsum())\n#         plt.title(\"Cumulative explained variance ratio\")\n#         plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))\n#         plt.xlabel('PCA component')\n#         plt.ylabel('Cumulative explained variance ratio')\n#         plt.show()\n        print(X.shape)\n        return X\n\npreprocessor = PreprocessMultiome()\n\nmulti_train_x = None\nstart, stop = 0, 6000\nmulti_train_x = preprocessor.fit_transform(pd.read_hdf(FP_MULTIOME_TRAIN_INPUTS, start=start, stop=stop).values)\n\nmulti_train_y = pd.read_hdf(FP_MULTIOME_TRAIN_TARGETS, start=start, stop=stop)\ny_columns = multi_train_y.columns\nmulti_train_y = multi_train_y.values\nprint(multi_train_y.shape)","metadata":{"execution":{"iopub.status.busy":"2022-08-25T01:41:19.366198Z","iopub.execute_input":"2022-08-25T01:41:19.366622Z","iopub.status.idle":"2022-08-25T01:41:53.141235Z","shell.execute_reply.started":"2022-08-25T01:41:19.366585Z","shell.execute_reply":"2022-08-25T01:41:53.139672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling&Prediction\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()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T01:41:53.151314Z","iopub.execute_input":"2022-08-25T01:41:53.151862Z","iopub.status.idle":"2022-08-25T01:41:53.382831Z","shell.execute_reply.started":"2022-08-25T01:41:53.151808Z","shell.execute_reply":"2022-08-25T01:41:53.38127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If there are many objective variables and many n_estimators, the model will become large, so it will not run to the end due to memory constraints.Confirm that it turns around 10","metadata":{}},{"cell_type":"code","source":"params={'learning_rate': 0.1,\n        'objective':'mae', \n        'metric':'mae',\n        'num_leaves': 256,\n        'verbose': -1 ,\n        \"seed\": 42,\n#         'bagging_fraction': 0.7,\n#         'feature_fraction': 0.7\n       }\nmodel = MultiOutputRegressor(lgb.LGBMRegressor(**params, n_estimators=12))\n\nmodel.fit(multi_train_x, multi_train_y)\n\ny_va_pred = model.predict(multi_train_x)\nmse = mean_squared_error(multi_train_y, y_va_pred)\nprint(mse)\n","metadata":{"execution":{"iopub.status.busy":"2022-08-25T01:41:53.384526Z","iopub.execute_input":"2022-08-25T01:41:53.384931Z","iopub.status.idle":"2022-08-25T02:26:23.866447Z","shell.execute_reply.started":"2022-08-25T01:41:53.384894Z","shell.execute_reply":"2022-08-25T02:26:23.86534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del multi_train_x, multi_train_y\ngc.collect()","metadata":{},"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":"code","source":"%%time\n# Read the table of rows and columns required for submission\neval_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))\neval_ids.cell_id = eval_ids.cell_id.astype(pd.CategoricalDtype())\neval_ids.gene_id = eval_ids.gene_id.astype(pd.CategoricalDtype())\ndisplay(eval_ids)\n\n# Create the set of needed cell_ids\ncell_id_set = set(eval_ids.cell_id)\n\n# Convert the string gene_ids to a more efficient categorical dtype\ny_columns = pd.CategoricalIndex(y_columns, dtype=eval_ids.gene_id.dtype, name='gene_id')\n","metadata":{"execution":{"iopub.status.busy":"2022-08-25T02:26:24.217581Z","iopub.execute_input":"2022-08-25T02:26:24.218342Z","iopub.status.idle":"2022-08-25T02:28:35.39384Z","shell.execute_reply.started":"2022-08-25T02:26:24.2183Z","shell.execute_reply":"2022-08-25T02:28:35.39158Z"},"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)\n# submission","metadata":{"execution":{"iopub.status.busy":"2022-08-25T02:28:35.396577Z","iopub.execute_input":"2022-08-25T02:28:35.397029Z","iopub.status.idle":"2022-08-25T02:28:35.508365Z","shell.execute_reply.started":"2022-08-25T02:28:35.396985Z","shell.execute_reply":"2022-08-25T02:28:35.507209Z"},"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":"%%time\n# Process the test data in chunks of 5000 rows\n\nstart = 0\nchunksize = 5000\ntotal_rows = 0\nwhile True:\n    multi_test_x = None # Free the memory if necessary\n    gc.collect()\n    # Read the 5000 rows and select the 30 % subset which is needed for the submission\n    multi_test_x = pd.read_hdf(FP_MULTIOME_TEST_INPUTS, start=start, stop=start+chunksize)\n    rows_read = len(multi_test_x)\n    needed_row_mask = multi_test_x.index.isin(cell_id_set)\n    multi_test_x = multi_test_x.loc[needed_row_mask]\n    \n    # Keep the index (the cell_ids) for later\n    multi_test_index = multi_test_x.index\n    \n    # Predict\n    multi_test_x = multi_test_x.values\n    multi_test_x = preprocessor.transform(multi_test_x)\n    test_pred = model.predict(multi_test_x)\n    \n    # Convert the predictions to a dataframe so that they can be matched with eval_ids\n    test_pred = pd.DataFrame(test_pred,\n                             index=pd.CategoricalIndex(multi_test_index,\n                                                       dtype=eval_ids.cell_id.dtype,\n                                                       name='cell_id'),\n                             columns=y_columns)\n    gc.collect()\n    \n    # Fill the predictions into the submission series row by row\n    for i, (index, row) in enumerate(test_pred.iterrows()):\n        row = row.reindex(eval_ids.gene_id[eval_ids.cell_id == index])\n        submission.loc[index] = row.values\n    print('na:', submission.isna().sum())\n\n    #test_pred_list.append(test_pred)\n    total_rows += len(multi_test_x)\n    print(total_rows)\n    if rows_read < chunksize: break # this was the last chunk\n    start += chunksize\n    \ndel multi_test_x, multi_test_index, needed_row_mask\n","metadata":{"execution":{"iopub.status.busy":"2022-08-25T02:28:35.510048Z","iopub.execute_input":"2022-08-25T02:28:35.510404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"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":"submission.reset_index(drop=True, inplace=True)\nsubmission.index.name = 'row_id'\nwith open(\"partial_submission_multi.pickle\", 'wb') as f: pickle.dump(submission, f)\nsubmission","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}