{"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":"code","source":"import numpy as np \nimport pandas as pd \n\nimport matplotlib.pyplot as plt\nimport os\n\n!pip install lightgbm\n\nimport os\nimport gc\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom collections import defaultdict\nfrom IPython.display import FileLink\n\n\nimport catboost as cb\nimport lightgbm as lgbm\nfrom sklearn.ensemble import GradientBoostingRegressor","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-01-07T10:17:35.929285Z","iopub.execute_input":"2023-01-07T10:17:35.929926Z","iopub.status.idle":"2023-01-07T10:17:52.536194Z","shell.execute_reply.started":"2023-01-07T10:17:35.929796Z","shell.execute_reply":"2023-01-07T10:17:52.534849Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Define model","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport xgboost as xgb\nimport operator\nimport warnings\n\n\n########################################################################################\n#\n# Main Class and Methods\n#\n########################################################################################\n\n\nclass BoostARoota(object):\n\n    def __init__(self, metric=None, clf=None, cutoff=4, iters=10, max_rounds=100, delta=0.1, silent=False):\n        self.metric = metric\n        self.clf = clf\n        self.cutoff = cutoff\n        self.iters = iters\n        self.max_rounds = max_rounds\n        self.delta = delta\n        self.silent = silent\n        self.keep_vars_ = None\n        \n        \n        #Throw errors if the inputted parameters don't meet the necessary criteria\n        if (metric is None) and (clf is None):\n            raise ValueError('you must enter one of metric or clf as arguments')\n        if cutoff <= 0:\n            raise ValueError('cutoff should be greater than 0. You entered' + str(cutoff))\n        if iters <= 0:\n            raise ValueError('iters should be greater than 0. You entered' + str(iters))\n        if (delta <= 0) | (delta > 1):\n            raise ValueError('delta should be between 0 and 1, was ' + str(delta))\n\n        #Issue warnings for parameters to still let it run\n        if (metric is not None) and (clf is not None):\n            warnings.warn('You entered values for metric and clf, defaulting to clf and ignoring metric')\n        if delta < 0.02:\n            warnings.warn(\"WARNING: Setting a delta below 0.02 may not converge on a solution.\")\n        if max_rounds < 1:\n            warnings.warn(\"WARNING: Setting max_rounds below 1 will automatically be set to 1.\")\n\n    def fit(self, x, y):\n        self.keep_vars_ = _BoostARoota(x, y,\n                                       metric=self.metric,\n                                       clf = self.clf,\n                                       cutoff=self.cutoff,\n                                       iters=self.iters,\n                                       max_rounds=self.max_rounds,\n                                       delta=self.delta,\n                                       silent=self.silent)\n        return self\n    \n    \n    def transform(self, x):\n        if self.keep_vars_ is None:\n            raise ValueError(\"You need to fit the model first\")\n        return x[self.keep_vars_]\n\n    def fit_transform(self, x, y):\n        self.fit(x, y)\n        return self.transform(x)\n\n########################################################################################\n#\n# Helper Functions to do the Heavy Lifting\n#\n########################################################################################\n\n\ndef _create_shadow(x_train):\n    \"\"\"\n    Take all X variables, creating copies and randomly shuffling them\n    :param x_train: the dataframe to create shadow features on\n    :return: dataframe 2x width and the names of the shadows for removing later\n    \"\"\"\n    x_shadow = x_train.copy()\n    for c in x_shadow.columns:\n        np.random.shuffle(x_shadow[c].values)\n    # rename the shadow\n    shadow_names = [\"ShadowVar\" + str(i + 1) for i in range(x_train.shape[1])]\n    x_shadow.columns = shadow_names\n    # Combine to make one new dataframe\n    new_x = pd.concat([x_train, x_shadow], axis=1)\n    return new_x, shadow_names\n\n########################################################################################\n#\n# BoostARoota\n#\n########################################################################################\n\n\ndef _reduce_vars_xgb(x, y, metric, this_round, cutoff, n_iterations, delta, silent):\n    \"\"\"\n    Function to run through each\n    :param x: Input dataframe - X\n    :param y: Target variable\n    :param metric: Metric to optimize in XGBoost\n    :param this_round: Round so it can be printed to screen\n    :return: tuple - stopping criteria and the variables to keep\n    \"\"\"\n    #Set up the parameters for running the model in XGBoost - split is on multi log loss\n    if metric == 'mlogloss':\n        param = {'objective': 'multi:softmax',\n                 'eval_metric': 'mlogloss',\n                 'num_class': len(np.unique(y)),\n                 'silent': 1}\n    else:\n        param = {'eval_metric': metric,\n                 'silent': 1}\n    for i in range(1, n_iterations+1):\n        # Create the shadow variables and run the model to obtain importances\n        new_x, shadow_names = _create_shadow(x)\n        dtrain = xgb.DMatrix(new_x, label=y)\n        bst = xgb.train(param, dtrain, verbose_eval=False)\n        if i == 1:\n            df = pd.DataFrame({'feature': new_x.columns})\n            pass\n\n        importance = bst.get_score(importance_type='weight')\n        importance = sorted(importance.items(), key=operator.itemgetter(1), reverse=True)\n        df2 = pd.DataFrame(importance, columns=['feature', 'fscore'+str(i)])\n        df2['fscore'+str(i)] = df2['fscore'+str(i)] / df2['fscore'+str(i)].sum()\n        df = pd.merge(df, df2, on='feature', how='outer')\n        if not silent:\n            print(\"Round: \", this_round, \" iteration: \", i)\n\n    df['Mean'] = df.mean(axis=1)\n    #Split them back out\n    real_vars = df[~df['feature'].isin(shadow_names)]\n    shadow_vars = df[df['feature'].isin(shadow_names)]\n\n    # Get mean value from the shadows\n    mean_shadow = shadow_vars['Mean'].mean() / cutoff\n    real_vars = real_vars[(real_vars.Mean > mean_shadow)]\n    real_vars.sort_values(by='Mean', inplace=True, ascending=False, ignore_index=True)\n\n    \n\n    #Check for the stopping criteria\n    #Basically looking to make sure we are removing at least 10% of the variables, or we should stop\n    if (len(real_vars['feature']) / len(x.columns)) > (1-delta):\n        criteria = True\n    else:\n        criteria = False\n\n    return criteria, real_vars['feature']\n\n\ndef _reduce_vars_sklearn(x, y, clf, this_round, cutoff, n_iterations, delta, silent):\n    \"\"\"\n    Function to run through each\n    :param x: Input dataframe - X\n    :param y: Target variable\n    :param clf: the fully specified classifier passed in by user\n    :param this_round: Round so it can be printed to screen\n    :return: tuple - stopping criteria and the variables to keep\n    \"\"\"\n    #Set up the parameters for running the model in XGBoost - split is on multi log loss\n\n    for i in range(1, n_iterations+1):\n        # Create the shadow variables and run the model to obtain importances\n        new_x, shadow_names = _create_shadow(x)\n        clf = clf.fit(new_x, np.ravel(y))\n\n        if i == 1:\n            df = pd.DataFrame({'feature': new_x.columns})\n            df2 = df.copy()\n            pass\n\n        try:\n            importance = clf.feature_importances_\n            df2['fscore' + str(i)] = importance\n        except ValueError:\n            print(\"this clf doesn't have the feature_importances_ method.  Only Sklearn tree based methods allowed\")\n\n        # importance = sorted(importance.items(), key=operator.itemgetter(1))\n\n        # df2 = pd.DataFrame(importance, columns=['feature', 'fscore'+str(i)])\n        df2['fscore'+str(i)] = df2['fscore'+str(i)] / df2['fscore'+str(i)].sum()\n        df = pd.merge(df, df2, on='feature', how='outer')\n        if not silent:\n            print(\"Round: \", this_round, \" iteration: \", i)\n\n    df['Mean'] = df.mean(axis=1)\n    #Split them back out\n    real_vars = df[~df['feature'].isin(shadow_names)]\n    shadow_vars = df[df['feature'].isin(shadow_names)]\n\n    # Get mean value from the shadows\n    mean_shadow = shadow_vars['Mean'].mean() / cutoff\n    real_vars = real_vars[(real_vars.Mean > mean_shadow)]\n    real_vars.sort_values(by='Mean', inplace=True, ascending=False, ignore_index=True)\n\n    #Check for the stopping criteria\n    #Basically looking to make sure we are removing at least 10% of the variables, or we should stop\n    if (len(real_vars['feature']) / len(x.columns)) > (1-delta):\n        criteria = True\n    else:\n        criteria = False\n\n    return criteria, real_vars['feature']\n\n#Main function exposed to run the algorithm\ndef _BoostARoota(x, y, metric, clf, cutoff, iters, max_rounds, delta, silent):\n    \"\"\"\n    Function loops through, waiting for the stopping criteria to change\n    :param x: X dataframe One Hot Encoded\n    :param y: Labels for the target variable\n    :param metric: The metric to optimize in XGBoost\n    :return: names of the variables to keep\n    \"\"\"\n\n    new_x = x.copy()\n    #Run through loop until \"crit\" changes\n    i = 0\n    weight = []\n    while True:\n        #Inside this loop we reduce the dataset on each iteration exiting with keep_vars\n        i += 1\n        \n        if clf is None:\n            crit, keep_vars = _reduce_vars_xgb(new_x,\n                                               y,\n                                               metric=metric,\n                                               this_round=i,\n                                               cutoff=cutoff,\n                                               n_iterations=iters,\n                                               delta=delta,\n                                               silent=silent)\n        else:\n            crit, keep_vars = _reduce_vars_sklearn(new_x,\n                                                   y,\n                                                   clf=clf,\n                                                   this_round=i,\n                                                   cutoff=cutoff,\n                                                   n_iterations=iters,\n                                                   delta=delta,\n                                                   silent=silent)\n\n        if crit | (i >= max_rounds):\n            break  # exit and use keep_vars as final variables\n        else:\n            new_x = new_x[keep_vars].copy()\n\n    if not silent:\n        print(\"BoostARoota ran successfully! Algorithm went through \", i, \" rounds.\")\n    return keep_vars","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:18:01.024137Z","iopub.execute_input":"2023-01-07T10:18:01.024808Z","iopub.status.idle":"2023-01-07T10:18:01.089931Z","shell.execute_reply.started":"2023-01-07T10:18:01.024754Z","shell.execute_reply":"2023-01-07T10:18:01.088260Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"code","source":"train_cite_input_df = pd.read_hdf(\"/kaggle/input/open-problems-multimodal/train_cite_inputs.h5\")\ntrain_cite_target_df = pd.read_hdf(\"/kaggle/input/open-problems-multimodal/train_cite_targets.h5\")\nmetadata_df = pd.read_csv(\"/kaggle/input/open-problems-multimodal/metadata.csv\")\n\nday_4_cell_ids = set(metadata_df[metadata_df[\"day\"] == 4][\"cell_id\"])\n\nday_4_test_index = train_cite_input_df.index[train_cite_input_df.index.isin(day_4_cell_ids)]\ntrain_index = train_cite_input_df.index[~train_cite_input_df.index.isin(day_4_cell_ids)]\n#X = pd.read_csv('/kaggle/input/preproccessed-data-for-boostaroota/X_forBoostARoota_day_4_EryP_CD36.csv')\ndf = pd.read_csv('/kaggle/input/research-project-01-around-multimodal-singlecell/HVG_cite_train_test_donor_3000.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:18:06.429658Z","iopub.execute_input":"2023-01-07T10:18:06.430178Z","iopub.status.idle":"2023-01-07T10:19:07.791042Z","shell.execute_reply.started":"2023-01-07T10:18:06.430115Z","shell.execute_reply":"2023-01-07T10:19:07.789578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"new_features = df['seurat_v3']\nnew_feat = train_cite_input_df[train_cite_input_df.columns[train_cite_input_df.columns.isin(new_features)]]","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:19:12.922028Z","iopub.execute_input":"2023-01-07T10:19:12.922537Z","iopub.status.idle":"2023-01-07T10:19:13.418292Z","shell.execute_reply.started":"2023-01-07T10:19:12.922500Z","shell.execute_reply":"2023-01-07T10:19:13.416979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#clf = 'xgb'\n#clf = cb.CatBoostRegressor()\n#clf= GradientBoostingRegressor()\nclf= lgbm.LGBMRegressor()","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:54:50.331607Z","iopub.execute_input":"2023-01-07T10:54:50.332047Z","iopub.status.idle":"2023-01-07T10:54:50.339035Z","shell.execute_reply.started":"2023-01-07T10:54:50.332013Z","shell.execute_reply":"2023-01-07T10:54:50.337609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cell_types = [\n#    \"HSC\",\n    \"EryP\",\n#    \"BP\",\n#    \"MasP\",\n#    \"MkP\",\n#    \"Mop\",\n#    \"NeuP\",\n#    \"hidden\"\n]","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:19:17.550342Z","iopub.execute_input":"2023-01-07T10:19:17.552438Z","iopub.status.idle":"2023-01-07T10:19:17.560673Z","shell.execute_reply.started":"2023-01-07T10:19:17.551515Z","shell.execute_reply":"2023-01-07T10:19:17.559003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"target_cds = [\n    \"CD36\",\n#    \"CD41\",\n#    \"CD32\",\n#    \"CD45\",\n#     \"CD88\",\n#     \"CD48\",\n#     \"CD62L\",\n#     \"CD49b\",\n#     \"CD45RA\",\n#     \"CD82\",\n#     \"CD38\",\n#     \"CD71\",\n#     \"CD115\",\n#     \"CD11a\",\n#     \"CD244\",\n]","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:19:18.899764Z","iopub.execute_input":"2023-01-07T10:19:18.900326Z","iopub.status.idle":"2023-01-07T10:19:18.907456Z","shell.execute_reply.started":"2023-01-07T10:19:18.900289Z","shell.execute_reply":"2023-01-07T10:19:18.905715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for cd_it in target_cds:\n    for cell_type_it in cell_types:\n        set_cell_type_it = set(metadata_df[metadata_df[\"cell_type\"] == cell_type_it][\"cell_id\"].tolist())\n        cell_type_test_index = day_4_test_index[day_4_test_index.isin(set_cell_type_it)]\n        cell_type_train_index = train_index[train_index.isin(set_cell_type_it)]    \n\n        del set_cell_type_it","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:19:20.755653Z","iopub.execute_input":"2023-01-07T10:19:20.756343Z","iopub.status.idle":"2023-01-07T10:19:20.891741Z","shell.execute_reply.started":"2023-01-07T10:19:20.756287Z","shell.execute_reply":"2023-01-07T10:19:20.890128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = new_feat.loc[cell_type_train_index]\ny = train_cite_target_df.loc[cell_type_train_index][cd_it]","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:19:23.013731Z","iopub.execute_input":"2023-01-07T10:19:23.014198Z","iopub.status.idle":"2023-01-07T10:19:23.187494Z","shell.execute_reply.started":"2023-01-07T10:19:23.014130Z","shell.execute_reply":"2023-01-07T10:19:23.186075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model","metadata":{}},{"cell_type":"code","source":"for i in range(34):\n    br = BoostARoota(metric=None,\n                     clf=clf,\n                     cutoff=4,\n                     iters=10,\n                     max_rounds=100,\n                     delta=0.1,\n                     silent=False\n                    )\n\n    br.fit(X,y)\n    br.keep_vars_.to_csv('br_lgb_'+str(i+66)+'_.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-07T10:37:08.900513Z","iopub.execute_input":"2023-01-07T10:37:08.900947Z","iopub.status.idle":"2023-01-07T10:54:50.329016Z","shell.execute_reply.started":"2023-01-07T10:37:08.900912Z","shell.execute_reply":"2023-01-07T10:54:50.327512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TO DO: \n# предсказать рнк(CD36) без мэйн белка(CD36) и наоборот\n# сд 36 и рнк из других белков\n# белки => рнк\n# сд 36 белок из всех рнк (без сд36)\n# проверка фиче импотанса cross-validation + holdout + проверку на то, что улучшение больше, чем std (НАМ НЕ ПОДХОДИТ)","metadata":{},"execution_count":null,"outputs":[]}]}