{"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":"# What is about ?\n\n### Briefly - downsample and make quick experiments\n\nHere we will downsample the data. I.e. select 10% of data as a kind of \"Playground\" - for quick experiments. \nAnd will train/tune different models on that playground first.\n\n\n### About CV scheme\n\nPublic and private test sets are quite different. (Private - contains new DAY (and donor), while public ONLY new donor).\nSo one should be careful with the validation schemes.\nSome proposals are described here: \n\nNotebook: https://www.kaggle.com/code/alexandervc/mmscel-crossvalidation-schemes\nTopic: https://www.kaggle.com/competitions/open-problems-multimodal/discussion/358860\n\nWe will be based on them. \n\n### Start with selection of \"Playground\" - 10% part of data - to quickly test models, ideas\n\nData is quite big, and so training models, tuning params might take long time. That is not always affordable.\nIn the present script we first choose some 10% part of data - to make quick experiments.\nThat part is chosen to have SAME proporitions of key characteristics: donors, days, cell types as the initial data\n\nSo we can create CV scheme for that \"playground\" part.\nBut we can also simplify even further - for start - use  not 6-fold scheme but just splite by days. And have only 1 train subset for quick experiments. \n\nTo simplify even further we can start play with just one target only. \n\n\n\n### Versions\n\n\n#### 1 Test several models - Ridge,LGB, SVR, RF etc...\n\n    Target chosen - only CD31\n    No params tuning - just fist look \n","metadata":{}},{"cell_type":"markdown","source":"# Install/import modules, load technical data\n","metadata":{}},{"cell_type":"code","source":"import time\nt0start = time.time()\n\nimport pandas as pd\nimport numpy as np\nimport os\nimport sys\n\nimport matplotlib.pyplot as plt\n#plt.style.use('dark_background')\nimport seaborn as sns\n\n#If you see a urllib warning running this cell, go to \"Settings\" on the right hand side, \n#and turn on internet. Note, you need to be phone verified.\n!pip install --quiet tables\n\n\nimport h5py\n!pip install hdf5plugin~=2.0 # https://forum.hdfgroup.org/t/cant-open-directory-usr-local-hdf5-lib-plugin/9738/4\nimport hdf5plugin\n\n# !pip install scanpy\n# import scanpy as sc\n# import anndata\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\")\n\ndf_cell = pd.read_csv(FP_CELL_METADATA)\ndf_cell","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:52:14.62687Z","iopub.execute_input":"2022-10-30T15:52:14.627222Z","iopub.status.idle":"2022-10-30T15:52:33.996076Z","shell.execute_reply.started":"2022-10-30T15:52:14.627193Z","shell.execute_reply":"2022-10-30T15:52:33.994832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import time\nt0start = time.time()\n\nimport pandas as pd\nimport numpy as np\nimport os\nimport sys\n\nimport matplotlib.pyplot as plt\n#plt.style.use('dark_background')\nimport seaborn as sns\n\n#If you see a urllib warning running this cell, go to \"Settings\" on the right hand side, \n#and turn on internet. Note, you need to be phone verified.\n!pip install --quiet tables\n\n\nimport h5py\n!pip install hdf5plugin~=2.0 # https://forum.hdfgroup.org/t/cant-open-directory-usr-local-hdf5-lib-plugin/9738/4\nimport hdf5plugin\n\n# !pip install scanpy\n# import scanpy as sc\n# import anndata\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\")\n\ndf_cell = pd.read_csv(FP_CELL_METADATA)\ndf_cell","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:52:33.998138Z","iopub.execute_input":"2022-10-30T15:52:33.99855Z","iopub.status.idle":"2022-10-30T15:52:52.858159Z","shell.execute_reply.started":"2022-10-30T15:52:33.998507Z","shell.execute_reply":"2022-10-30T15:52:52.856908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load prepared Features for CITE-seq part of task","metadata":{}},{"cell_type":"code","source":"%%time\n\nprint('Load prepared features for CITE-seq')\n# These files contain both train and test parts .\n# For CITEseq part - first 70988 elements - train, and later 48663 - test. Overall 119651 samples.\nfn = '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_train_and_test_TruncatedSVD200_niter7_rs42.csv'\nfn = '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_train_and_test_PCA500.csv'\ndf_cite = pd.read_csv(fn,index_col = 0)\ndisplay(df_cite)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:52:52.861491Z","iopub.execute_input":"2022-10-30T15:52:52.86189Z","iopub.status.idle":"2022-10-30T15:53:10.576538Z","shell.execute_reply.started":"2022-10-30T15:52:52.861852Z","shell.execute_reply":"2022-10-30T15:53:10.575492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_cite.mean(axis = 0)","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:53:10.577964Z","iopub.execute_input":"2022-10-30T15:53:10.578689Z","iopub.status.idle":"2022-10-30T15:53:10.738748Z","shell.execute_reply.started":"2022-10-30T15:53:10.57865Z","shell.execute_reply":"2022-10-30T15:53:10.737727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_cite.std(axis=0)","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:53:10.740473Z","iopub.execute_input":"2022-10-30T15:53:10.740928Z","iopub.status.idle":"2022-10-30T15:53:11.362272Z","shell.execute_reply.started":"2022-10-30T15:53:10.740893Z","shell.execute_reply":"2022-10-30T15:53:11.361214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Targets for CITE-seq","metadata":{}},{"cell_type":"code","source":"%%time\n#if 1:\nprint('Load CITE-seq targets and ')\ndf_cite_train_y = pd.read_hdf(FP_CITE_TRAIN_TARGETS)\ndisplay(df_cite_train_y)","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:53:11.364249Z","iopub.execute_input":"2022-10-30T15:53:11.36473Z","iopub.status.idle":"2022-10-30T15:53:12.100588Z","shell.execute_reply.started":"2022-10-30T15:53:11.364685Z","shell.execute_reply":"2022-10-30T15:53:12.099524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load and start prepare some metadata (cut only CITE-seq train part)","metadata":{}},{"cell_type":"code","source":"%%time\nfn2 = '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/_citeseq_meta_all_text_also.csv'\ndf_meta_full = pd.read_csv(fn2,index_col = 0)\ndisplay(df_meta_full)\n#if 1:\ndf_meta = pd.DataFrame(index = df_cite_train_y.index) \ndf_meta = df_meta.join(df_cell.set_index('cell_id') )\ndisplay(df_meta)","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:53:12.102152Z","iopub.execute_input":"2022-10-30T15:53:12.102786Z","iopub.status.idle":"2022-10-30T15:53:12.412446Z","shell.execute_reply.started":"2022-10-30T15:53:12.102748Z","shell.execute_reply":"2022-10-30T15:53:12.411487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create \"Playground\"  (i.e. downsample)\n\nSmall 10% of data of data where we can make prelimanary experiments. \n\nit will be labeled by special column \"Playground\" in df_meta (meta data for CITE-seq train only)","metadata":{}},{"cell_type":"code","source":"# Prepare for creation of additional holdout folds with 10% of samples \n# We will use stratified Kfold to achieve that days, cell_types and donors are equally distributed \nimport numpy as np\nfrom sklearn.model_selection import StratifiedKFold\nscol = 'donor&day&CT'\ndf_meta[scol] =df_meta['donor'].apply(lambda x:str(x)+'_') + df_meta['day'].apply(lambda x:str(x)+'_') + df_meta['cell_type']\n\n\nskf = StratifiedKFold(n_splits=10,  shuffle=True, random_state=40)\nskf.get_n_splits(df_meta, df_meta[scol] )\n\n\n\ny = df_meta[scol] \nfor train_index, test_index in skf.split(df_meta, df_meta[scol]):\n    print(\"TRAIN:\", len(train_index), \"TEST:\", len(test_index) ); \n    break\nprint(test_index)\nprint(df_meta[scol].value_counts().head(5)   )\nprint(df_meta.iloc[test_index,:][scol].value_counts().head(5)    )\n\n\nflagged_column_name = 'Playground'\ndf_meta[flagged_column_name] = 0 \ndf_meta.loc[df_meta.index[test_index],flagged_column_name]  = 1\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:53:12.414218Z","iopub.execute_input":"2022-10-30T15:53:12.415024Z","iopub.status.idle":"2022-10-30T15:53:12.653409Z","shell.execute_reply.started":"2022-10-30T15:53:12.414984Z","shell.execute_reply":"2022-10-30T15:53:12.652308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create X,y, X_train, y_train, etc - INSIDE \"Playground\"\n","metadata":{}},{"cell_type":"code","source":"selected_target = 'CD31'\n\nX= df_cite.iloc[:70988,:][df_meta['Playground']==1]\ny= df_cite_train_y[df_meta['Playground']==1][selected_target]\nprint('X.shape, y.shape', X.shape, y.shape )\n\n# Create simplfied validation scheme - like real test data - with two test-sets private-like, public-like:\n# Private like test - new DAY, and donor, \n# While public like - only new donor (days are the same as in train):\n# Step 1: \nmask_train = (df_meta['Playground']==1)&(df_meta['day']!=4)&(df_meta['donor']!=31800) \nX_train = df_cite.iloc[:70988,:][mask_train]\ny_train = df_cite_train_y[mask_train][ selected_target ]\n# Step 2:\nmask_test_private_like = (df_meta['Playground']==1)&(df_meta['day']==4)\nX_test_private_like = df_cite.iloc[:70988,:][ mask_test_private_like  ]\ny_test_private_like = df_cite_train_y[mask_test_private_like][ selected_target ]\nX_test = X_test_private_like\ny_test = y_test_private_like\n# Step 3: \nmask_test_public_like = (df_meta['Playground']==1)&(df_meta['day']!=4)  &(df_meta['donor']==31800) \nX_test_public_like = df_cite.iloc[:70988,:][mask_test_public_like ]\ny_test_public_like = df_cite_train_y[mask_test_public_like][ selected_target ]\n\nX.shape,y.shape, X_train.shape, X_test_private_like.shape, X_test_public_like.shape, y_train.shape, y_test_private_like.shape, y_test_public_like.shape\n\n","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:53:12.656463Z","iopub.execute_input":"2022-10-30T15:53:12.656842Z","iopub.status.idle":"2022-10-30T15:53:12.780408Z","shell.execute_reply.started":"2022-10-30T15:53:12.656815Z","shell.execute_reply":"2022-10-30T15:53:12.779386Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CB model","metadata":{}},{"cell_type":"code","source":"%%time\nfrom hyperopt import hp, tpe, Trials\nfrom hyperopt.fmin import fmin\nimport hyperopt\nfrom sklearn.model_selection import KFold\nfrom catboost import CatBoost, CatBoostRegressor, Pool\nfrom sklearn import metrics\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import r2_score","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:53:55.64623Z","iopub.execute_input":"2022-10-30T15:53:55.647185Z","iopub.status.idle":"2022-10-30T15:53:55.693842Z","shell.execute_reply.started":"2022-10-30T15:53:55.647135Z","shell.execute_reply":"2022-10-30T15:53:55.692846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nmodel = CatBoostRegressor(task_type=\"CPU\", iterations=30, verbose=True, random_state=42)#**params)\nmodel.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2022-10-30T16:12:57.323308Z","iopub.execute_input":"2022-10-30T16:12:57.324189Z","iopub.status.idle":"2022-10-30T16:13:01.424036Z","shell.execute_reply.started":"2022-10-30T16:12:57.324143Z","shell.execute_reply":"2022-10-30T16:13:01.423018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_pred = model.predict(X_test)\nr2_score(y_test, y_pred), mean_squared_error(y_test, y_pred)","metadata":{"execution":{"iopub.status.busy":"2022-10-30T16:13:01.429447Z","iopub.execute_input":"2022-10-30T16:13:01.431846Z","iopub.status.idle":"2022-10-30T16:13:01.475211Z","shell.execute_reply.started":"2022-10-30T16:13:01.431809Z","shell.execute_reply":"2022-10-30T16:13:01.474373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#check if given parameter can be interpreted as a numerical value\ndef is_number(s):\n    if s is None:\n        return False\n    try:\n        float(s)\n        return True\n    except ValueError:\n        return False\n\n#convert given set of paramaters to integer values\n#this at least cuts the excess float decimals if they are there\ndef convert_int_params(names, params):\n    for int_type in names:\n        #sometimes the parameters can be choices between options or numerical values. like \"log2\" vs \"1-10\"\n        raw_val = params[int_type]\n        if is_number(raw_val):\n            params[int_type] = int(raw_val)\n    return params\n\n#convert float parameters to 3 digit precision strings\n#just for simpler diplay and all\ndef convert_float_params(names, params):\n    for float_type in names:\n        raw_val = params[float_type]\n        if is_number(raw_val):\n            params[float_type] = '{:.3f}'.format(raw_val)\n    return params\n\n\n# how many CV folds to do on the data\nn_folds = 5\n# max number of rows to use for X and y. to reduce time and compare options faster\nmax_n = None\n# max number of trials hyperopt runs\nn_trials = 200\n#verbosity in LGBM is how often progress is printed. with 100=print progress every 100 rounds. 0 is quite?\nverbosity = False\nprint_summary = False\n\nall_scores = []\nall_params = []\n\ndef fit_cb(params):\n    train_dataset = Pool(X_train, y_train)\n    val_dataset = Pool(X_test_public_like, y_test_public_like)\n\n    model = CatBoostRegressor(**params)\n\n    model.fit(train_dataset, eval_set = val_dataset, early_stopping_rounds = 15, \n              use_best_model = True, verbose = 0)\n    oof_preds = model.predict(X_test_private_like)\n    score = mean_squared_error(y_test_private_like, oof_preds)\n    features = X.columns\n\n    all_scores.append(score)\n    all_params.append(params)\n    if print_summary:\n        print(f\"score: {score}\")\n    return score\n\n\n# this is the objective function the hyperopt aims to minimize\n# i call it objective_sklearn because the lgbm functions called use sklearn API\ndef objective_sklearn(params):\n    int_types = [\"depth\", \"max_bin\"]\n    params = convert_int_params(int_types, params)\n\n    # Extract the boosting type\n    params['boosting_type'] = params['boosting_type']['boosting_type']\n    #    print(\"running with params:\"+str(params))\n    \n    score = fit_cb(params)\n    if verbosity == 0:\n        if print_summary:\n            print(\"Score {:.3f}\".format(score))\n    else:\n        print(\"Score {:.3f} params {}\".format(score, params))\n    result = {\"loss\": score, \"score\": score, \"params\": params, 'status': hyperopt.STATUS_OK}\n    return result\n\ndef optimize_cb(max_n_search=None, boosting_type=None):\n    space = {\n        'boosting_type': hp.choice('boosting_type',\n                                   [{'boosting_type': 'Plain',\n                                     }]),\n#         'num_leaves': hp.quniform('num_leaves', 4, 127, 12),\n        'depth': hp.quniform('max_depth', 5, 12, 1),\n        'learning_rate': hp.uniform('learning_rate', 0.1, 0.2),\n        'l2_leaf_reg': hp.uniform('l2_leaf_reg', 0, 2), \n        'loss_function': 'RMSE', \n        'eval_metric': 'RMSE', \n        'task_type': 'CPU', \n        'iterations': 100,\n        'od_type': 'Iter', \n        'bootstrap_type': 'Bayesian', \n        'allow_const_label': False, \n        'bagging_temperature': hp.uniform('bagging_temperature', 1, 2), \n        'random_state': 42,\n        'od_wait': 20,\n        'max_bin': hp.quniform('max_bin', 500, 1501, 100),\n    }\n    \n    if boosting_type:\n        space['boosting_type'] = boosting_type\n\n    global max_n\n    max_n = max_n_search\n    trials = Trials()\n    best = fmin(fn=objective_sklearn,\n                space=space,\n                algo=tpe.suggest,\n                max_evals=n_trials,\n                trials=trials,\n                verbose= 1)\n\n    # find the trial with lowest loss value. this is what we consider the best one\n    idx = np.argmin(trials.losses())\n    print(idx)\n\n    print(trials.trials[idx])\n\n    # these should be the training parameters to use to achieve the best score in best trial\n    params = trials.trials[idx][\"result\"][\"params\"]\n    max_n = None\n\n    print('==============================')\n    print('= PARAMS')\n    print('==============================')\n    print(params)\n    return params, trials\n\n# run a search\ndef run_lgb(max_n=120000):\n    # the param is the number of rows to use for training\n    params, trials = optimize_lgbm(max_n)\n    print(params)\n\n    return params, trials","metadata":{"execution":{"iopub.status.busy":"2022-10-30T16:17:25.693706Z","iopub.execute_input":"2022-10-30T16:17:25.694061Z","iopub.status.idle":"2022-10-30T16:17:25.713228Z","shell.execute_reply.started":"2022-10-30T16:17:25.69403Z","shell.execute_reply":"2022-10-30T16:17:25.71225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"params, trials = optimize_cb()\n# boosting_type={\"boosting_type\": \"rf\"}","metadata":{"execution":{"iopub.status.busy":"2022-10-30T16:17:27.012625Z","iopub.execute_input":"2022-10-30T16:17:27.013389Z","iopub.status.idle":"2022-10-30T16:20:27.877097Z","shell.execute_reply.started":"2022-10-30T16:17:27.013326Z","shell.execute_reply":"2022-10-30T16:20:27.87448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nmodel = CatBoostRegressor(**params)#**params)\nmodel.fit(train_dataset)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_pred = model.predict(X_test)\nr2_score(y_test, y_pred), mean_squared_error(y_test, y_pred)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"execution":{"iopub.status.busy":"2022-10-28T11:40:31.076516Z","iopub.execute_input":"2022-10-28T11:40:31.077004Z","iopub.status.idle":"2022-10-28T11:40:31.0841Z","shell.execute_reply.started":"2022-10-28T11:40:31.076961Z","shell.execute_reply":"2022-10-28T11:40:31.082408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}