{"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":"# CITEseq LGBM + Optuna Baseline\nThis notebook is based on https://www.kaggle.com/code/swimmy/lgbm-baseline-msci-citeseq.\n\n# Please vote it if useful!\n\n# Optuna\nOptuna is a great hyperparameter optimization framework. \nI use Optuna to find params below:\n```\nparams = {'metric': 'mae', 'random_state': 42, 'n_estimators': 2000, 'reg_alpha': 0.03645857751758206, 'reg_lambda': 0.0025972855120393492, 'colsample_bytree': 1.0, 'subsample': 0.6, 'learning_rate': 0.013262872399411381, 'max_depth': 10, 'num_leaves': 186, 'min_child_samples': 263, 'min_data_per_groups': 46}\n```\nAnd this Notebook score 0.830 in leaderboard.","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\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\nfrom sklearn.linear_model import Lasso\nfrom sklearn.multioutput import MultiOutputRegressor\nimport lightgbm as lgb\nfrom sklearn.compose import TransformedTargetRegressor\nfrom sklearn.metrics import mean_squared_error\nimport optuna\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-26T13:07:44.107562Z","iopub.execute_input":"2022-08-26T13:07:44.108214Z","iopub.status.idle":"2022-08-26T13:07:47.679736Z","shell.execute_reply.started":"2022-08-26T13:07:44.108086Z","shell.execute_reply":"2022-08-26T13:07:47.678727Z"},"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-26T13:07:47.681534Z","iopub.execute_input":"2022-08-26T13:07:47.682124Z","iopub.status.idle":"2022-08-26T13:08:00.562716Z","shell.execute_reply.started":"2022-08-26T13:07:47.682087Z","shell.execute_reply":"2022-08-26T13:08:00.560823Z"},"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-26T13:08:00.565359Z","iopub.execute_input":"2022-08-26T13:08:00.565848Z","iopub.status.idle":"2022-08-26T13:08:01.047158Z","shell.execute_reply.started":"2022-08-26T13:08:00.565799Z","shell.execute_reply":"2022-08-26T13:08:01.045255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cross-validation\n\nThe note I referred to had the following description, but I confirmed that 13000 rows can be rotated in memory, so I changed columns_to_use = 13000. In addition, the search is performed by changing the starting point of the line to be acquired.\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- Of the remaining columns, we keep only the last 12000.\n- We do a PCA and keep only the 240 most important components.\n- We use PCA(copy=False), which overwrites its input in fit_transform().\n- We fit a ridge regression model with 70988\\*240 inputs and 70988\\*140 outputs. ","metadata":{}},{"cell_type":"code","source":"%%time\n# Preprocessing\ncol_start = 10000\n\nclass PreprocessCiteseq(BaseEstimator, TransformerMixin):\n    columns_to_use = 13000\n    \n    @staticmethod\n    def take_column_subset(X):\n        return X[:,-(PreprocessCiteseq.columns_to_use+col_start):-col_start]\n    \n    def transform(self, X):\n        print(X.shape)\n        X = X[:,~self.all_zero_columns]\n        print(X.shape)\n        X = PreprocessCiteseq.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        gc.collect()\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 = PreprocessCiteseq.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=240, copy=False, random_state=1)\n        X = self.pca.fit_transform(X)\n        print(X.shape)\n        return X\n\npreprocessor = PreprocessCiteseq()\n\ncite_train_x = None\ncite_train_x = preprocessor.fit_transform(pd.read_hdf(FP_CITE_TRAIN_INPUTS).values)\n\ncite_train_y = pd.read_hdf(FP_CITE_TRAIN_TARGETS).values\nprint(cite_train_y.shape)","metadata":{"execution":{"iopub.status.busy":"2022-08-26T13:08:01.053433Z","iopub.execute_input":"2022-08-26T13:08:01.054661Z","iopub.status.idle":"2022-08-26T13:10:15.525740Z","shell.execute_reply.started":"2022-08-26T13:08:01.054598Z","shell.execute_reply":"2022-08-26T13:10:15.524327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Optuna","metadata":{"execution":{"iopub.status.busy":"2022-08-27T00:02:09.267311Z","iopub.execute_input":"2022-08-27T00:02:09.267864Z","iopub.status.idle":"2022-08-27T00:02:09.293108Z","shell.execute_reply.started":"2022-08-27T00:02:09.267747Z","shell.execute_reply":"2022-08-27T00:02:09.291532Z"}}},{"cell_type":"code","source":"small_train_x, small_train_y =  cite_train_x[:1000,:],cite_train_y[:1000,:]\ndef objective(trial):\n    params = {\n        'metric': 'mae', \n        'random_state': 42,\n        'n_estimators': 2000,\n        'reg_alpha': trial.suggest_loguniform('reg_alpha', 1e-3, 10.0),\n        'reg_lambda': trial.suggest_loguniform('reg_lambda', 1e-3, 10.0),\n        'colsample_bytree': trial.suggest_categorical('colsample_bytree', [0.3,0.4,0.5,0.6,0.7,0.8,0.9, 1.0]),\n        'subsample': trial.suggest_categorical('subsample', [0.4,0.5,0.6,0.7,0.8,1.0]),\n        'learning_rate': trial.suggest_loguniform('learning_rate', 0.01, 0.1),\n        'max_depth': trial.suggest_categorical('max_depth', [10,20,100]),\n        'num_leaves' : trial.suggest_int('num_leaves', 1, 1000),\n        'min_child_samples': trial.suggest_int('min_child_samples', 1, 300),\n        'cat_smooth' : trial.suggest_int('min_data_per_groups', 1, 100)\n    }\n    model = MultiOutputRegressor(lgb.LGBMRegressor(**params))\n\n    model.fit(small_train_x, small_train_y)\n\n    y_va_pred = model.predict(cite_train_x)\n    mse = mean_squared_error(cite_train_y, y_va_pred)\n    \n    return mse","metadata":{"execution":{"iopub.status.busy":"2022-08-26T13:10:15.536896Z","iopub.execute_input":"2022-08-26T13:10:15.537370Z","iopub.status.idle":"2022-08-26T13:10:15.549391Z","shell.execute_reply.started":"2022-08-26T13:10:15.537331Z","shell.execute_reply":"2022-08-26T13:10:15.548160Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"find_params = False\nif find_params:\n    study = optuna.create_study(\n        direction='minimize', \n        pruner=optuna.pruners.MedianPruner(n_warmup_steps=20),\n        study_name='small')\n    study.optimize(objective, n_trials=20)","metadata":{"execution":{"iopub.status.busy":"2022-08-26T12:45:00.181175Z","iopub.execute_input":"2022-08-26T12:45:00.181630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling&Prediction\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":"params = {'metric': 'mae', 'random_state': 42, 'n_estimators': 2000, 'reg_alpha': 0.03645857751758206, 'reg_lambda': 0.0025972855120393492, 'colsample_bytree': 1.0, 'subsample': 0.6, 'learning_rate': 0.013262872399411381, 'max_depth': 10, 'num_leaves': 186, 'min_child_samples': 263, 'min_data_per_groups': 46}\nmodel = MultiOutputRegressor(lgb.LGBMRegressor(**params))\n\nmodel.fit(cite_train_x, cite_train_y)\n\ny_va_pred = model.predict(cite_train_x)\nmse = mean_squared_error(cite_train_y, y_va_pred)\n\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2022-08-26T13:12:57.905055Z","iopub.execute_input":"2022-08-26T13:12:57.905912Z","iopub.status.idle":"2022-08-26T20:53:24.531678Z","shell.execute_reply.started":"2022-08-26T13:12:57.905868Z","shell.execute_reply":"2022-08-26T20:53:24.526953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del cite_train_x, cite_train_y\ngc.collect()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_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","metadata":{"execution":{"iopub.status.busy":"2022-08-24T01:13:01.665885Z","iopub.execute_input":"2022-08-24T01:13:01.666729Z","iopub.status.idle":"2022-08-24T01:14:01.684733Z","shell.execute_reply.started":"2022-08-24T01:13:01.666685Z","shell.execute_reply":"2022-08-24T01:14:01.683388Z"},"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-24T01:14:01.686208Z","iopub.execute_input":"2022-08-24T01:14:01.686618Z","iopub.status.idle":"2022-08-24T01:14:01.794494Z","shell.execute_reply.started":"2022-08-24T01:14:01.686585Z","shell.execute_reply":"2022-08-24T01:14:01.793329Z"},"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","metadata":{"execution":{"iopub.status.busy":"2022-08-24T01:14:01.79634Z","iopub.execute_input":"2022-08-24T01:14:01.796828Z","iopub.status.idle":"2022-08-24T01:16:06.586099Z","shell.execute_reply.started":"2022-08-24T01:14:01.796778Z","shell.execute_reply":"2022-08-24T01:16:06.58487Z"},"trusted":true},"execution_count":null,"outputs":[]}]}