{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59094,"databundleVersionId":7010844,"sourceType":"competition"}],"dockerImageVersionId":30587,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about \n\nHere is a componet of U900 team (#13) solution for the \"Open Problems – Single-Cell Perturbations\" challenge.\n\nSee main writeup here: https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/460858\n\nThe present notebook implements the CatBoost part of the solution.\nIt achieves .584 (0.776 private) scores - version 64 of the present notebook.\n\nModeling based on the TSVD-30 with Target Encoding by Qunatile 80 \nModel Id: tsvd30_modelCATB_NI250_MD6_LR0.03_SS1_CS0.5_encQuantileEncoder_quantile0.8\n\nTraining: without CD8 and removing 3 random samples each time. \n","metadata":{}},{"cell_type":"markdown","source":"## Preliminaries","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport time\nt0start = time.time() \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport pandas as pd\nimport tensorflow as tf\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-11-29T05:49:05.714438Z","iopub.execute_input":"2023-11-29T05:49:05.714790Z","iopub.status.idle":"2023-11-29T05:49:19.523772Z","shell.execute_reply.started":"2023-11-29T05:49:05.714762Z","shell.execute_reply":"2023-11-29T05:49:19.522230Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\ndf_de_train = pd.read_parquet(fn)# , index_col = 0)\nprint(df_de_train.shape)\ndisplay(df_de_train )\n\nplt.figure(figsize = (20,4) )\nv = df_de_train.iloc[:,5:].max(axis = 0 ).sort_values(ascending = False, key = abs )\nplt.plot(v.values,'*-')\nplt.title('Max-abs DE for genes',fontsize = 20 )\nplt.grid()\nplt.show()\ndisplay(v.head(15))\n\n# %%time\nfn = '/kaggle/input/open-problems-single-cell-perturbations/id_map.csv'\ndf_id_map = pd.read_csv(fn,index_col = 0)\nprint(df_id_map.shape)\ndisplay(df_id_map)\nfn = '/kaggle/input/open-problems-single-cell-perturbations/sample_submission.csv'\ndf_sample_submit = pd.read_csv(fn, index_col = 0)\nprint(df_sample_submit.shape)\ndisplay( df_sample_submit )\n\n\nfeatures_columns = ['cell_type','sm_name']\n# selected_features_columns = ['cell_type','sm_name']  # [ 'sm_name'] # What features to consider - others not used \n#     # For OP2 task we have two features - cell-type and compound - both categorical\n#     # Some simple models (like Ridge) better work when one uses only one feature - 'sm_name' (=\"compound\") - so you can try that option. \n#     # But seems pyboost is powerful enough to proper work with the both features.\n#     # Table with results for Ridge can be found here: https://www.kaggle.com/code/alexandervc/op2-target-encoders?scriptVersionId=148576642&cellId=1\n#     # Or here: https://docs.google.com/spreadsheets/d/1APN63PMaWZygVjYimK9Ivt0RvifdAU5JRYkxiDn4szw/edit?usp=sharing (sheet \"Encoders\") and discussed here:\n    \nX_submit_categorical = df_id_map[features_columns] # selected_features_columns] # pd.DataFrame(df_id_map, columns= features_columns )\nX_full_categorical = df_de_train[features_columns]# selected_features_columns ]\n\nY_full = df_de_train.iloc[:,5:].values\nprint('X_submit_categorical.shape, X_full_categorical.shape ,  Y_full.shape', X_submit_categorical.shape, X_full_categorical.shape ,  Y_full.shape)","metadata":{"execution":{"iopub.status.busy":"2023-11-29T05:49:19.526479Z","iopub.execute_input":"2023-11-29T05:49:19.527102Z","iopub.status.idle":"2023-11-29T05:49:24.823917Z","shell.execute_reply.started":"2023-11-29T05:49:19.527046Z","shell.execute_reply":"2023-11-29T05:49:24.822170Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_de_train['cell_type'].unique() )\ndf_de_train['cell_type'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-11-29T05:49:24.826250Z","iopub.execute_input":"2023-11-29T05:49:24.826733Z","iopub.status.idle":"2023-11-29T05:49:24.839403Z","shell.execute_reply.started":"2023-11-29T05:49:24.826690Z","shell.execute_reply":"2023-11-29T05:49:24.838669Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Custom CV-schemes\n\nSee the notebook: https://www.kaggle.com/code/alexandervc/op2-class-for-custom-cv-schemes#Ridge-CV-scoring.-Simple-example\n\nMore analysis of the CV schemes can be found in slides: https://docs.google.com/presentation/d/1wiz0Wmt4D54pqMMsIOyJHuQYMZ3hTBZQQnjbLzwoGYY/edit?usp=sharing and sheet: https://docs.google.com/spreadsheets/d/1APN63PMaWZygVjYimK9Ivt0RvifdAU5JRYkxiDn4szw/edit?usp=sharing\n\n\nCV schemes (see discussion https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/444494 ):\n\nAmbrosM: https://www.kaggle.com/code/ambrosm/scp-quickstart\n\nMT: https://www.kaggle.com/code/masato114/scp-quickstart-another-cv-strategy\n\nKishan's notebook: https://www.kaggle.com/code/liudacheldieva/neural-network-regression/notebook#Predicting-on-test-data ","metadata":{}},{"cell_type":"code","source":"%%time\nimport random \nfrom sklearn.model_selection import KFold\n\n# Prepare indexes for AmbrosM scheme:\n# For AmbrosM scheme: \nfolds_index_data_AmbrosM = [ ]\ntrain_sm_names = ['Idelalisib', 'Crizotinib', 'Linagliptin', 'Palbociclib', 'Dabrafenib', 'Alvocidib', 'LDN 193189', 'R428', 'Porcn Inhibitor III', \n  'Belinostat', 'Foretinib', 'MLN 2238', 'Penfluridol', 'Dactolisib', 'O-Demethylated Adapalene', 'Oprozomib (ONX 0912)', 'CHIR-99021']\nlist_fold_ids =  ['NK cells', 'T cells CD4+', 'T cells CD8+', 'T regulatory cells'] \nfor fold_id in list_fold_ids:\n    mask_va = (df_de_train.cell_type == fold_id) & ~df_de_train.sm_name.isin(train_sm_names)\n    mask_tr = ~mask_va # 485 or 487 training rows\n    IX_train = np.where( mask_tr > 0 )[0]\n    IX_test = np.where( mask_va > 0 )[0]\n    #print(fold_id,  len(IX_test), type(IX_test), IX_test[:3], len(IX_train), type(IX_train), IX_train[:3]  )\n    folds_index_data_AmbrosM.append( [IX_train, IX_test ])\n\n# For MT CV scheme:\nfolds_index_data_MT = []\nfold_to_compounds = {0: ['Alvocidib', 'Belinostat', 'Foretinib', 'LDN 193189',  'Linagliptin', 'O-Demethylated Adapalene'],\n 1: ['Dabrafenib', 'Dactolisib', 'Idelalisib', 'MLN 2238', 'Palbociclib', 'Porcn Inhibitor III'],\n 2: ['CHIR-99021', 'Crizotinib', 'Oprozomib (ONX 0912)', 'Penfluridol',  'R428']}\nfor fold_id in [0,1,2]:\n    mask_va = df_de_train['cell_type'].isin(['Myeloid cells', 'B cells']) & df_de_train['sm_name'].isin(fold_to_compounds[fold_id])\n    mask_tr = ~mask_va    \n    IX_train = np.where( mask_tr > 0 )[0]\n    IX_test = np.where( mask_va > 0 )[0]\n    #print(fold_id,  len(IX_test), type(IX_test), IX_test[:3], len(IX_train), type(IX_train), IX_train[:3]  )\n    folds_index_data_MT.append( [IX_train, IX_test ])    \n    \n    \ndict_folds_Antonina = {0: [4, 15, 18, 41, 47, 53, 56, 57, 59, 61, 62, 66, 78, 83, 86, 93, 100, 110, 115, 116, 117, 120, 123, 131, 148, 152, 161, 185, 189, 193, 197, 198, 205, 207, 209, 210, 213, 218, 219, 225, 228, 230, 250, 252, 257, 263, 286, 292, 294, 304, 306, 310, 325, 327, 336, 338, 339, 350, 351, 356, 357, 360, 362, 369, 379, 385, 395, 419, 424, 430, 432, 434, 438, 441, 443, 445, 448, 452, 456, 465, 467, 471, 472, 474, 498, 508, 523, 534, 542, 544, 574, 581, 584, 589, 601, 607], 1: [1, 3, 5, 14, 24, 38, 40, 49, 50, 54, 55, 60, 79, 82, 87, 102, 112, 114, 118, 119, 127, 132, 146, 149, 166, 179, 183, 184, 186, 190, 192, 201, 203, 204, 216, 227, 242, 253, 261, 264, 271, 289, 295, 296, 299, 300, 309, 323, 324, 329, 333, 337, 354, 387, 390, 393, 406, 408, 415, 416, 420, 421, 422, 436, 442, 449, 453, 458, 462, 463, 469, 484, 487, 489, 490, 492, 497, 501, 505, 506, 509, 513, 515, 516, 546, 548, 566, 569, 585, 587, 588, 592, 593, 595, 600, 605], 2: [0, 16, 19, 44, 48, 51, 52, 58, 63, 65, 70, 81, 85, 88, 89, 92, 111, 113, 122, 125, 135, 137, 139, 150, 164, 175, 191, 202, 206, 211, 217, 220, 221, 229, 239, 243, 247, 248, 251, 258, 259, 260, 267, 270, 274, 281, 288, 290, 291, 303, 322, 326, 328, 332, 347, 359, 364, 370, 383, 384, 389, 392, 401, 423, 428, 454, 455, 461, 464, 494, 496, 504, 510, 511, 524, 525, 526, 530, 540, 541, 551, 565, 568, 570, 572, 573, 575, 577, 579, 580, 597, 598, 603, 606, 610, 611, 613], 3: [21, 22, 26, 36, 42, 45, 46, 64, 67, 69, 71, 80, 84, 90, 101, 103, 134, 136, 140, 147, 151, 160, 163, 165, 167, 174, 176, 177, 180, 181, 182, 187, 188, 195, 196, 199, 200, 212, 215, 226, 231, 240, 241, 249, 255, 256, 262, 265, 269, 285, 287, 297, 298, 307, 312, 319, 321, 330, 340, 342, 344, 345, 348, 358, 366, 368, 382, 391, 400, 403, 404, 405, 425, 427, 431, 433, 435, 444, 450, 451, 457, 459, 460, 473, 485, 499, 500, 502, 503, 507, 512, 514, 527, 528, 529, 531, 543, 547, 553, 554, 564, 571, 576, 583, 586, 591, 596, 602, 604, 608, 609], 4: [2, 6, 7, 17, 20, 23, 25, 27, 28, 29, 37, 39, 43, 68, 91, 121, 124, 126, 128, 129, 130, 133, 138, 153, 162, 178, 194, 208, 214, 222, 223, 224, 232, 244, 245, 246, 254, 266, 268, 272, 273, 282, 283, 284, 293, 301, 302, 305, 308, 311, 320, 331, 334, 335, 341, 343, 346, 349, 352, 353, 355, 361, 363, 365, 367, 377, 378, 380, 381, 386, 388, 394, 396, 397, 398, 399, 402, 407, 417, 418, 426, 429, 437, 439, 440, 446, 447, 466, 468, 470, 481, 482, 483, 486, 488, 491, 493, 495, 532, 533, 545, 549, 550, 552, 555, 562, 563, 567, 578, 582, 590, 594, 599, 612], 5: [8, 9, 10, 11, 12, 13, 30, 31, 32, 33, 34, 35, 72, 73, 74, 75, 76, 77, 94, 95, 96, 97, 98, 99, 104, 105, 106, 107, 108, 109, 141, 142, 143, 144, 145, 154, 155, 156, 157, 158, 159, 168, 169, 170, 171, 172, 173, 233, 234, 235, 236, 237, 238, 275, 276, 277, 278, 279, 280, 313, 314, 315, 316, 317, 318, 371, 372, 373, 374, 375, 376, 409, 410, 411, 412, 413, 414, 475, 476, 477, 478, 479, 480, 517, 518, 519, 520, 521, 522, 535, 536, 537, 538, 539, 556, 557, 558, 559, 560, 561]}\n\nfolds_index_data_Antonina = []    \nfor i in range(5):  \n    l_valid = np.array( dict_folds_Antonina[i] )\n    l_train = np.array( [k for k in range(614) if k not in l_valid ] )\n    folds_index_data_Antonina.append( [l_train, l_valid ])    \nprint( len(folds_index_data_Antonina ) )\n\nclass KFold_custom:\n    '''\n    Class with similar to sklearn \"KFold\" class, interface.\n    Support of CV schemes proposed by AmbrosM, MT, Kishan and standard sklearn KFold\n    Supports options to return IX_train,IX_valid,IX_test - triplet - if valid_size > 0\n    Examples:\n    kf = KFold_custom('AmbrosM'):     \n    for i_fold,(IX_train,IX_test)   in enumerate( kf.split()) :\n        pass\n    kf = KFold_custom('AmbrosM', valid_size = 0.1 )     \n    for i_fold,(IX_train ,IX_valid,IX_test)   in enumerate( kf.split()) :\n        pass\n        \n    kf = KFold_custom(CV_scheme = 'Random') # wrapper for sklearn Kfolds with shuffle = True\n    kf = KFold_custom(CV_scheme = CV_scheme, random_state = 42 ) # fixing random_state ensures reprodicibility of folds \n    '''\n    def __init__(self, CV_scheme,  valid_size = 0, n_splits=5, random_state = 42, verbose = 0 ):\n        self.CV_scheme = CV_scheme\n        self.CV_scheme_inf =  CV_scheme\n        self.valid_size = np.clip(valid_size,0,1)\n        self.random_state = random_state\n\n        self.folds_index_data = [ [np.arange(0), np.arange(0), np.arange(0)] ] # Example data\n        self.n_splits = 1 # Example data\n        if CV_scheme == 'Kishan1':\n            # One fold scheme which contains train, valid, test parts \n            IX_train_Kishan1 = np.arange(429)\n            IX_val_Kishan1 = np.arange(429,521)\n            IX_test_Kishan1 = np.arange(521,614)\n            self.n_splits = 1\n            self.folds_index_data =  [   [IX_train_Kishan1, IX_val_Kishan1, IX_test_Kishan1]  ]\n        elif CV_scheme == 'Tonya':\n            self.n_splits = 5\n            self.folds_index_data = folds_index_data_Antonina\n        elif CV_scheme == 'AmbrosM':\n            self.n_splits = 4\n            self.folds_index_data = folds_index_data_AmbrosM\n        elif CV_scheme == 'MT':\n            self.n_splits = 3\n            self.folds_index_data = folds_index_data_MT\n        elif 'Random'.lower() in CV_scheme.lower():\n            self.n_splits = n_splits\n            self.CV_scheme_inf = 'Random_'+str(n_splits) +'_'+str(random_state )\n            kf = KFold(n_splits=n_splits, random_state = random_state, shuffle=True )#,\n            self.folds_index_data = list( kf.split( np.arange(614) ) )\n        elif 'Full'.lower() in CV_scheme.lower():\n            # Just return the full set as both train and test - it useful to re-train the model on the entire data\n            IX_train_full = np.arange(614)\n            IX_test_full = np.arange(614)\n            self.n_splits = 1\n            self.folds_index_data = [   [IX_train_full, IX_test_full]  ]\n        else:\n            s = 'Uncrecognized CV_scheme ' + str(CV_scheme) \n            raise ValueError(s)\n\n        self.verbose = verbose \n        if verbose >= 10:\n            print(self.CV_scheme, self.n_splits, len(self.folds_index_data) )\n            #print(self.folds_index_data)\n            \n    def split(self, X=None):\n        '''\n        X - NOT used, just for compatibility with sklearn \n        '''\n        for item in self.folds_index_data:\n            if self.CV_scheme in ['Kishan1']:\n                if self.valid_size > 0:\n                    yield item[0],item[1],item[2]\n                else :\n                    yield np.array( list(set(item[0])|set(item[1]))  ) , item[2] # Return train, test only \n            elif self.CV_scheme in ['Tonya', 'AmbrosM','MT', 'Random', 'Full']:\n                if self.valid_size == 0:\n                    yield item[0],item[1]\n                else:\n                    # Split \"full-train\"->( real-train, valid ) \n                    index_for_real_train_part = int(  (1-self.valid_size) * len(item[0]) )\n                    if self.random_state is None:\n                        IX_train =  item[0][:index_for_real_train_part]\n                        IX_valid =  item[0][index_for_real_train_part:]\n                    elif self.random_state == -1:\n                        p = np.random.permutation(len(item[0])) \n                        IX_train =  item[0][p][:index_for_real_train_part]\n                        IX_valid =  item[0][p][index_for_real_train_part:]\n                    else:\n                        np.random.seed(self.random_state) # Set temporary random seed\n                        p = np.random.permutation(len(item[0])) \n                        np.random.seed(random.randint(0,30000)) # Randomize seed again - use Python random, not numpy       \n                        IX_train =  item[0][p][:index_for_real_train_part]\n                        IX_valid =  item[0][p][index_for_real_train_part:]\n                    yield IX_train, IX_valid, item[1]\n                    \n        \n    def get_list_main_CV_schemes(self):\n        return ['AmbrosM','MT','Kishan1','Random']\n    def get_n_splits(self, X=np.array([])):\n        return self.n_splits","metadata":{"execution":{"iopub.status.busy":"2023-11-29T05:49:24.841453Z","iopub.execute_input":"2023-11-29T05:49:24.841884Z","iopub.status.idle":"2023-11-29T05:49:25.150351Z","shell.execute_reply.started":"2023-11-29T05:49:24.841851Z","shell.execute_reply":"2023-11-29T05:49:25.149459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kf = KFold_custom('Random',n_splits = 20,valid_size = 0 )     \nfor i_fold,(IX_train,IX_test)   in enumerate( kf.split()) :\n    print(i_fold, len(IX_train),len(IX_test))","metadata":{"execution":{"iopub.status.busy":"2023-11-29T05:55:20.323729Z","iopub.execute_input":"2023-11-29T05:55:20.324126Z","iopub.status.idle":"2023-11-29T05:55:20.332397Z","shell.execute_reply.started":"2023-11-29T05:55:20.324096Z","shell.execute_reply":"2023-11-29T05:55:20.331268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Aux: Get model, encoder, reducer - wrapper functions","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.metrics import r2_score\nfrom sklearn.decomposition import TruncatedSVD\nfrom sklearn.decomposition import FastICA\nfrom sklearn.decomposition import TruncatedSVD\nfrom sklearn.linear_model import Ridge\nfrom sklearn.svm import LinearSVR\nfrom sklearn.svm import SVR\nfrom sklearn.kernel_ridge import KernelRidge\n\nimport lightgbm as lgb\nimport catboost\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.ensemble import ExtraTreesRegressor\n\nfrom sklearn.multioutput import MultiOutputRegressor\n\n\nimport category_encoders as ce\n\nimport warnings\n\n\n\n\n\ndef get_encoder( main_config_model_feature_etc ):\n    '''\n    Returns category encoder, intialized with appropriate params from the input config\n    \n    enc, str_encoder_id   =  get_encoder( main_config_model_feature_etc )\n    '''\n    \n    features_mode = main_config_model_feature_etc['encoder']\n    str_encoder_id =  main_config_model_feature_etc['encoder']\n    if features_mode == 'TargetEncoder':\n        # Target encoding for categorical features.\n        # For the case of continuous target: \n        # features are replaced with a blend of the expected value of the target given particular categorical value \n        # and the expected value of the target over all the training data.\n        # The blend weight is controlled by \"smoothing\"\n        smoothing = main_config_model_feature_etc.get('smoothing', 10)\n        min_samples_leaf = main_config_model_feature_etc.get('min_samples_leaf', 20)\n        enc = ce.TargetEncoder(smoothing = smoothing, min_samples_leaf = min_samples_leaf  )\n    elif features_mode == 'QuantileEncoder':\n        # Quantile Encoding for categorical features.\n        # This a statistically modified version of target MEstimate encoder \n        # where selected features are replaced by the statistical quantile instead of the mean. \n        # Replacing with the median is a particular case where self.quantile = 0.5. \n        # In comparison to MEstimateEncoder it has two tunable parameter m and quantile\n\n        quantile = main_config_model_feature_etc.get('quantile', .5)\n        m = main_config_model_feature_etc.get('smoothing', 1) \n        m = main_config_model_feature_etc.get('m', 1) \n        # this is the “m” in the m-probability estimate. Higher value of m results into stronger shrinking. M is non-negative. 0 for no smoothing.\n        enc = ce.QuantileEncoder(quantile = quantile, m = m  )\n        str_encoder_id += str( int(100*quantile) )\n    elif features_mode == 'CatBoostEncoder':\n        # CatBoost Encoding for categorical features. \n        \n        a = main_config_model_feature_etc.get('smoothing', 1) \n        a = main_config_model_feature_etc.get('a', 1) \n        # additive smoothing (it is the same variable as “m” in m-probability estimate). By default set to 1.\n        random_state = main_config_model_feature_etc.get('random_state', None) \n        sigma = main_config_model_feature_etc.get('sigma', None)\n        # adds normal (Gaussian) distribution noise into training data in order to decrease overfitting (testing data are untouched). \n        # sigma gives the standard deviation (spread or “width”) of the normal distribution.\n        # See example here: https://www.kaggle.com/code/alexandervc/op2-target-encoders/notebook#Example:-add-Gaussian-noise-into-training-data-(%22sigma%22-parameter-of-encoder)\n        enc = ce.CatBoostEncoder(a = a, sigma = sigma  , random_state = random_state)\n    elif features_mode == 'JamesSteinEncoder':\n        #  James-Stein estimator - it is a kind of target encoder but uses theoretical  estimate for smoothing parameter.\n        sigma = main_config_model_feature_etc.get('sigma', None)\n        if (sigma is not None) and (sigma > 0):             randomized = True\n        else: randomized = False \n        # adds normal (Gaussian) distribution noise into training data in order to decrease overfitting (testing data are untouched). \n        # sigma gives the standard deviation (spread or “width”) of the normal distribution.\n        # It works only if randomized = True \n        # See example here: https://www.kaggle.com/code/alexandervc/op2-target-encoders/notebook#Example:-add-Gaussian-noise-into-training-data-(%22sigma%22-parameter-of-encoder)\n        random_state = main_config_model_feature_etc.get('random_state', None) \n        enc = ce.JamesSteinEncoder( sigma = sigma  , randomized = randomized, random_state = random_state)\n        \n    elif features_mode == 'LeaveOneOutEncoder':\n        # This is very similar to target encoding but excludes the current row’s target when calculating the mean target for a level to reduce the effect of outliers.\n        \n        sigma = main_config_model_feature_etc.get('sigma', None)\n        # adds normal (Gaussian) distribution noise into training data in order to decrease overfitting (testing data are untouched). \n        # sigma gives the standard deviation (spread or “width”) of the normal distribution.\n        # It works only if randomized = True \n        # See example here: https://www.kaggle.com/code/alexandervc/op2-target-encoders/notebook#Example:-add-Gaussian-noise-into-training-data-(%22sigma%22-parameter-of-encoder)\n        random_state = main_config_model_feature_etc.get('random_state', None) \n        enc = ce.LeaveOneOutEncoder( sigma = sigma  , random_state = random_state)\n        \n    elif features_mode in  ['OneHotEncoder', 'HelmertEncoder', 'BackwardDifferenceEncoder', 'CountEncoder', 'OrdinalEncoder'  ]:\n        enc = getattr(ce, features_mode )()\n        \n    return enc, str_encoder_id\n\ndef get_reducer( main_config_model_feature_etc , verbose = 0): \n    '''\n    Returns the dimensional reduction method initialized with appropriate params from the config.\n    \n    reducer, str_reducer_id  = get_reducer( main_config_model_feature_etc , verbose = 0)\n    '''\n    \n    str_reducer_id =  str( main_config_model_feature_etc['reducer'] )\n    if main_config_model_feature_etc['reducer']=='tsvd':\n        n_components = main_config_model_feature_etc.get('n_components',30)\n        reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n        str_reducer_id = 'tsvd'+str( n_components )\n    elif main_config_model_feature_etc['reducer']=='ica':\n        n_components = main_config_model_feature_etc.get('n_components',30)\n        reducer = FastICA(n_components=n_components, random_state=0, whiten='unit-variance')\n        str_reducer_id = 'ica'+str( n_components )\n    elif main_config_model_feature_etc['reducer']=='pca':\n        n_components = main_config_model_feature_etc.get('n_components',30)\n        reducer =  PCA(n_components=n_components)\n        str_reducer_id = 'pca'+str( n_components )\n    \n    if verbose >= 100:\n        print( str_reducer_id )\n        print( reducer )\n        print( main_config_model_feature_etc )\n\n    return reducer, str_reducer_id    \n\n\ndef get_model(main_config_model_feature_etc , verbose = 0):\n    '''\n    Returns model specified by the config\n    \n    model, str_model_id = get_model(main_config_model_feature_etc , verbose = 0)\n    '''\n    \n    str_model_id =  str( main_config_model_feature_etc['model'] )\n    if main_config_model_feature_etc['model'] == 'Ridge':\n        alpha = main_config_model_feature_etc.get('alpha',1)\n        fit_intercept = main_config_model_feature_etc.get('fit_intercept',True)\n        positive=main_config_model_feature_etc.get('positive',False)\n        model = Ridge( alpha ,  fit_intercept = fit_intercept,  positive = positive,)\n        str_model_id  = 'Ridge'+str(alpha)\n    elif main_config_model_feature_etc['model'] == 'KRRrbf':\n        # main_config_model_feature_etc = {'model': 'KRR', 'alpha':1, 'kernel': 'linear', 'gamma': None,  'degree': 3,  'coef0': 1, 'reducer': 'tsvd', 'n_components': n_components, 'features_mode': 'target_enc_i_th_target', 'list_features_in': ['cell_type', 'sm_name']}\n        alpha = main_config_model_feature_etc.get('alpha',1)\n        kernel = 'rbf'\n        gamma = main_config_model_feature_etc.get('gamma',None ) # Gamma parameter for the RBF, laplacian, polynomial, exponential chi2 and sigmoid kernels. Interpretation of the default value is left to the kernel; see the documentation for sklearn.metrics.pairwise. Ignored by other kernels.\n        coef0 = main_config_model_feature_etc.get('coef0',1.0 )  #  Independent term in kernel function. It is only significant in ‘poly’ and ‘sigmoid’.\n        degree = main_config_model_feature_etc.get('degree',3 ) ##Degree of the polynomial kernel function (‘poly’). Must be non-negative. Ignored by all other kernels.\n        model = KernelRidge(alpha=alpha, kernel = kernel,  gamma = gamma, coef0=coef0, degree = degree)\n    elif main_config_model_feature_etc['model'] in [ 'KRR', 'KRRrbf', 'KRRlin' ] :\n        # main_config_model_feature_etc = {'model': 'KRR', 'alpha':1, 'kernel': 'linear', 'gamma': None,  'degree': 3,  'coef0': 1, 'reducer': 'tsvd', 'n_components': n_components, 'features_mode': 'target_enc_i_th_target', 'list_features_in': ['cell_type', 'sm_name']}\n        alpha = main_config_model_feature_etc.get('alpha',1)\n        kernel = main_config_model_feature_etc.get('kernel','linear' ) # {‘linear’, ‘poly’, ‘rbf’, ‘sigmoid’, ‘precomputed’} or callable, default=’rbf’\n        if main_config_model_feature_etc['model'] == 'KRRrbf': kernel = 'rbf'\n        elif main_config_model_feature_etc['model'] == 'KRRlin': kernel = 'linear'\n        gamma = main_config_model_feature_etc.get('gamma',None ) # Gamma parameter for the RBF, laplacian, polynomial, exponential chi2 and sigmoid kernels. Interpretation of the default value is left to the kernel; see the documentation for sklearn.metrics.pairwise. Ignored by other kernels.\n        coef0 = main_config_model_feature_etc.get('coef0',1.0 )  #  Independent term in kernel function. It is only significant in ‘poly’ and ‘sigmoid’.\n        degree = main_config_model_feature_etc.get('degree',3 ) ##Degree of the polynomial kernel function (‘poly’). Must be non-negative. Ignored by all other kernels.\n        model = KernelRidge(alpha=alpha, kernel = kernel,  gamma = gamma, coef0=coef0, degree = degree)\n        \n    elif main_config_model_feature_etc['model'] == 'LSVR':\n        C = main_config_model_feature_etc.get('C',1)#  default=1.0 Regularization parameter. The strength of the regularization is inversely proportional to C. Must be strictly positive.\n        max_iter = main_config_model_feature_etc.get('max_iter',1000) # sklearn default default=1000 \n        epsilon = main_config_model_feature_etc.get('epsilon',0.1) # Epsilon parameter in the epsilon-insensitive loss function. Note that the value of this parameter depends on the scale of the target variable y. If unsure, set epsilon=0. \n        # 0.1 is used: https://www.kaggle.com/code/mehrankazeminia/1-op2-eda-linearsvr-regressorchain?scriptVersionId=145592242&cellId=88\n        model = MultiOutputRegressor( LinearSVR(max_iter= max_iter, epsilon= epsilon, C=C) )\n        \n    elif main_config_model_feature_etc['model'] in ['SVR','SVRrbf','SVRlin' ]:\n        # main_config_model_feature_etc = {'model':'SVR', 'kernel':'rbf','C':1 }\n        kernel = main_config_model_feature_etc.get('kernel','rbf' ) # {‘linear’, ‘poly’, ‘rbf’, ‘sigmoid’, ‘precomputed’} or callable, default=’rbf’\n        if 'rbf' in main_config_model_feature_etc['model']: kernel = 'rbf'\n        elif 'lin' in main_config_model_feature_etc['model']: kernel = 'linear'\n        C = main_config_model_feature_etc.get('C',1)#  default=1.0 Regularization parameter. The strength of the regularization is inversely proportional to C. Must be strictly positive.\n        epsilon = main_config_model_feature_etc.get('epsilon',0.1) # Epsilon parameter in the epsilon-insensitive loss function. Note that the value of this parameter depends on the scale of the target variable y. If unsure, set epsilon=0. \n        gamma = main_config_model_feature_etc.get('gamma','scale' ) # Kernel coefficient for ‘rbf’, ‘poly’ and ‘sigmoid’. \n        coef0 = main_config_model_feature_etc.get('coef0',0.0 )  #  Independent term in kernel function. It is only significant in ‘poly’ and ‘sigmoid’.\n        tol = main_config_model_feature_etc.get('tol',0.001 ) #\n        degree = main_config_model_feature_etc.get('degree',3 ) ##Degree of the polynomial kernel function (‘poly’). Must be non-negative. Ignored by all other kernels.\n        max_iter = main_config_model_feature_etc.get('max_iter',-1) # sklearn default default=-1 - it may cause too long\n        shrinking = main_config_model_feature_etc.get('shrinking', True) # Whether to use the shrinking heuristic. See the User Guide.        \n        model = MultiOutputRegressor( SVR(kernel = kernel, max_iter= max_iter, epsilon= epsilon, C=C, gamma = gamma, coef0=coef0, degree = degree, tol = tol, shrinking=shrinking) )\n        \n\n    elif main_config_model_feature_etc['model'] == 'CATB':\n        params_loc = main_config_model_feature_etc.get('CATB_params',{} )\n        params_loc['loss_function'] = main_config_model_feature_etc.get('loss_function', 'RMSE' )\n        categorical_features = main_config_model_feature_etc.get('categorical_features',[])\n        for prm_name in ['iterations', 'depth','learning_rate','subsample', 'colsample_bylevel', 'min_data_in_leaf', 'random_strength','random_seed','random_state', 'l2_leaf_reg' ]: #  \n            if prm_name in main_config_model_feature_etc.keys(): params_loc[prm_name] = main_config_model_feature_etc[prm_name]\n        model = MultiOutputRegressor( CatBoostRegressor(cat_features=categorical_features, verbose = 0, **params_loc ) )  \n        \n        # https://forecastegy.com/posts/catboost-hyperparameter-tuning-guide-with-optuna/#subsample-subsample\n        # def objective(trial):\n        #     params = {\n        #         \"iterations\": 1000,\n        #         \"learning_rate\": trial.suggest_float(\"learning_rate\", 1e-3, 0.1, log=True),\n        #         \"depth\": trial.suggest_int(\"depth\", 1, 10),\n        #         \"subsample\": trial.suggest_float(\"subsample\", 0.05, 1.0),\n        #         \"colsample_bylevel\": trial.suggest_float(\"colsample_bylevel\", 0.05, 1.0),\n        #         \"min_data_in_leaf\": trial.suggest_int(\"min_data_in_leaf\", 1, 100),\n        #     }\n\n        #     model = cb.CatBoostRegressor(**params, silent=True)\n        #     model.fit(X_train, y_train)\n        #     predictions = model.predict(X_val)\n        #     rmse = mean_squared_error(y_val, predictions, squared=False)\n        #     return rmse\n\n            \n    elif main_config_model_feature_etc['model'] == 'RFR':\n        #n_estimators=100,*, criterion='squared_error', max_depth=None, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features=1.0, max_leaf_nodes=None, min_impurity_decrease=0.0, bootstrap=True, oob_score=False, n_jobs=None, random_state=None, verbose=0, warm_start=False, ccp_alpha=0.0, max_samples=None\n        params_loc = main_config_model_feature_etc.get('RFR_params',{} )\n        for prm_name in ['criterion', 'n_estimators', 'max_depth', 'min_samples_split', 'min_samples_split', 'min_samples_leaf', 'min_weight_fraction_leaf',  'max_features', 'max_leaf_nodes','min_impurity_decrease','random_state','ccp_alpha','max_samples' ]:\n            if prm_name in main_config_model_feature_etc.keys(): params_loc[prm_name] = main_config_model_feature_etc[prm_name]\n        model = MultiOutputRegressor( RandomForestRegressor( **params_loc )  )\n        \n    elif main_config_model_feature_etc['model'] == 'ETR':\n        # n_estimators=100, *, criterion='squared_error', max_depth=None, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features=1.0, max_leaf_nodes=None, min_impurity_decrease=0.0, bootstrap=False, oob_score=False, n_jobs=None, random_state=None, verbose=0, warm_start=False, ccp_alpha=0.0, max_samples=None\n        params_loc = main_config_model_feature_etc.get('ETR_params',{} )\n        for prm_name in ['criterion', 'n_estimators', 'max_depth', 'min_samples_split', 'min_samples_split', 'min_samples_leaf', 'min_weight_fraction_leaf',  'max_features', 'max_leaf_nodes','min_impurity_decrease','random_state','ccp_alpha','max_samples' ]:\n            if prm_name in main_config_model_feature_etc.keys(): params_loc[prm_name] = main_config_model_feature_etc[prm_name]\n        model = MultiOutputRegressor(  ExtraTreesRegressor( **params_loc )  )\n        \n\n    elif main_config_model_feature_etc['model'] == 'LGB':\n        params_loc = main_config_model_feature_etc.get('LGB_params',{} )\n        for prm_name in ['n_estimators', 'max_depth','learning_rate', 'colsample_bytree', 'subsample', 'random_state', 'reg_alpha',  'reg_lambda', 'num_leaves','min_child_samples' ]:\n            if prm_name in main_config_model_feature_etc.keys(): params_loc[prm_name] = main_config_model_feature_etc[prm_name]\n        model = MultiOutputRegressor(  lgb.LGBMRegressor( **params_loc ) )\n    elif main_config_model_feature_etc['model'] == 'LGBcv1_036':\n        # Params for LGB just on two categorical features as it is  found by optimization in https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning#Optuna+LightGBM\n        # But LB score is terrible - \n        params_best1_cv1_036 = {'random_state': 0, 'n_estimators': 20, 'reg_alpha': 6.764079452929363, 'reg_lambda': 0.41900776876588564, \n        'colsample_bytree': 0.3, 'subsample': 0.7, 'max_depth': 1, 'learning_rate': 0.08456104070184789,  'num_leaves': 682, 'min_child_samples': 102}\n        model = MultiOutputRegressor( lgb.LGBMRegressor( **params_best1_cv1_036 )  )\n        \n    if verbose >= 100:\n        print( str_model_id )\n        print( model )\n        print( main_config_model_feature_etc )\n        \n    return model, str_model_id\n\ndef get_brief_string_info_on_config(main_config_model_feature_etc):\n\n    str_inf_cfg = ''\n    if 'reducer' in main_config_model_feature_etc.keys():\n        str_inf_cfg += main_config_model_feature_etc['reducer']\n        n_components = main_config_model_feature_etc.get('n_components', 25)\n        str_inf_cfg += str(n_components)\n        \n    dict_abbreviations = { 'iterations': 'NI','depth':'MD',  'learning_rate':'LR', 'subsample': 'SS',  'colsample_bylevel':'CS',  'encoder':'enc'    }    \n    for key_loc in ['model','alpha', 'iterations', 'depth',  'learning_rate', 'subsample', 'colsample_bylevel',  'encoder', 'quantile'  ]: # 'iterations': 500, 'depth': 3,  'learning_rate': 0.03, 'subsample': 1, 'colsample_bylevel': 0.5, \n        if key_loc in main_config_model_feature_etc.keys():\n            abbrevate_loc = str( dict_abbreviations.get( key_loc,  key_loc ) )\n            str_inf_cfg += '_'+abbrevate_loc+str(main_config_model_feature_etc[key_loc])\n            \n#     if 'sm_name' in list_features_to_encode: str_inf_cfg += 'Compound'\n#     if 'cell_type' in list_features_to_encode: str_inf_cfg += 'CellType'\n    \n    return str_inf_cfg\n\n\ncfg1 = {'model':'Ridge', 'reducer': 'tsvd',  'n_components': 2, 'encoder':'TargetEncoder'   }\nprint('cfg1', cfg1)\nstr_inf_cfg = get_brief_string_info_on_config(cfg1)\nprint(str_inf_cfg)\nmain_config_model_feature_etc_589 = {'reducer': 'tsvd', 'n_components': 30, 'model': 'CATB', 'encoder': 'QuantileEncoder', 'iterations': 250, 'depth': 6, 'learning_rate': 0.03, 'subsample': 1, 'colsample_bylevel': 0.5, 'min_data_in_leaf': 5, 'loss_function': 'RMSE', 'random_strength': None, 'random_seed': 42, 'l2_leaf_reg': 1, 'm': 1, 'quantile': 0.8}\nstr_inf_cfg = get_brief_string_info_on_config(main_config_model_feature_etc_589)\nprint(str_inf_cfg)\n\n\nmodel, str_model_id = get_model({'model':'Ridge'}, verbose = 100)\nprint(model, str_model_id  ); print()\nmodel, str_model_id = get_model({'model':'CATB'}, verbose = 100)\nprint(model, str_model_id  ); print()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T08:01:14.900289Z","iopub.execute_input":"2023-11-27T08:01:14.901139Z","iopub.status.idle":"2023-11-27T08:01:17.632440Z","shell.execute_reply.started":"2023-11-27T08:01:14.901104Z","shell.execute_reply":"2023-11-27T08:01:17.631164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Specify Model and params","metadata":{}},{"cell_type":"code","source":"%%time\nif 1:\n    # Catboost 0.589 Final config:\n    # https://www.kaggle.com/code/alexandervc/op2-gentle-param-tuner?scriptVersionId=150416991&cellId=36\n    main_config_model_feature_etc_589 = {'reducer': 'tsvd', 'n_components': 30, 'model': 'CATB', 'encoder': 'QuantileEncoder', 'iterations': 250, 'depth': 6, 'learning_rate': 0.03, 'subsample': 1, 'colsample_bylevel': 0.5, 'min_data_in_leaf': 5, 'loss_function': 'RMSE', 'random_strength': None, 'random_seed': 42, 'l2_leaf_reg': 1, 'm': 1, 'quantile': 0.8}\n    main_config_model_feature_etc = main_config_model_feature_etc_589\n\n#     main_config_model_feature_etc = {'reducer': 'tsvd', 'n_components': 30, 'encoder':'LeaveOneOutEncoder',# 'QuantileEncoder',  'quantile': 0.8 ,  'm': 1,  \n#                                      'model': 'CATB', 'iterations': 250, 'depth': 6, \n#                                      'learning_rate': 0.03, 'subsample': 1, 'colsample_bylevel': 0.5, 'min_data_in_leaf': 5, 'loss_function': 'RMSE', 'random_strength': None, 'random_seed': 42, \n#                                      'l2_leaf_reg': 1,}\n\n    print( main_config_model_feature_etc )\n    str_inf_cfg = get_brief_string_info_on_config(main_config_model_feature_etc)\n    print( str_inf_cfg )\n    \n    verbose = 10\n\n\n\n    model, str_model_id = get_model(main_config_model_feature_etc,  verbose = 100)\n    print(model, str_model_id  ); print()\n    #str_inf_cfg += str_model_id\n\n    reducer, str_reducer_id  = get_reducer( main_config_model_feature_etc , verbose = 0)\n    #str_inf_cfg += str_reducer_id\n    \n    enc, str_encoder_id   =  get_encoder( main_config_model_feature_etc )\n    #str_inf_cfg += '_'+ str_encoder_id\n\n    list_features_to_encode = ['cell_type', 'sm_name']# [ 'sm_name']  # [ 'sm_name'] # What features to consider - others not used \n    #     # For OP2 task we have two features - cell-type and compound - both categorical\n    #     # Some simple models (like Ridge) better work when one uses only one feature - 'sm_name' (=\"compound\") - so you can try that option. \n    #     # But seems pyboost is powerful enough to proper work with the both features.\n    #     # Table with results for Ridge can be found here: https://www.kaggle.com/code/alexandervc/op2-target-encoders?scriptVersionId=148576642&cellId=1\n    #     # Or here: https://docs.google.com/spreadsheets/d/1APN63PMaWZygVjYimK9Ivt0RvifdAU5JRYkxiDn4szw/edit?usp=sharing (sheet \"Encoders\") and discussed here:\n\n    if 'sm_name' in list_features_to_encode: str_inf_cfg += '_Dr'\n    if 'cell_type' in list_features_to_encode: str_inf_cfg += '_CT'\n\n\n\n    n_selfblend = 1\n    str_inf_cfg += '_SB'+str(n_selfblend)\n\n\n    CV_scheme = 'Tonya' #  'MT'\n    valid_size = 0 \n    kf = KFold_custom(CV_scheme, valid_size = valid_size)   \n    print('CV_scheme:', CV_scheme)\n    str_inf_cfg += '_'+CV_scheme\n\n\n    print( str_inf_cfg )","metadata":{"execution":{"iopub.status.busy":"2023-11-27T08:01:17.633983Z","iopub.execute_input":"2023-11-27T08:01:17.634848Z","iopub.status.idle":"2023-11-27T08:01:17.650886Z","shell.execute_reply.started":"2023-11-27T08:01:17.634798Z","shell.execute_reply":"2023-11-27T08:01:17.648880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif 0:\n    from sklearn.linear_model import Ridge\n    from sklearn.decomposition import TruncatedSVD\n    from sklearn.metrics import r2_score\n    from sklearn.metrics import mean_absolute_error\n    import category_encoders as ce\n\n    verbose = 10\n\n\n    str_inf_cfg = ''\n\n\n    alpha = 2e5 # 2e5 seems optimum for  LOO  CT&Drug     (seems for various tsvd n_comps - same is true)\n    # alpha = 5e4 # 5e4 seems optimum for  LOO  drug only   (seems for various tsvd n_comps - same is true)\n    model = Ridge(alpha=alpha)\n    str_model_id = 'Ridge'+str(alpha)\n    str_inf_cfg += str_model_id\n\n    n_components = 70# 450\n    reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n    str_inf_cfg += '_tsvd'+str(n_components)\n\n    # n_components = 300 # 300 seems to be optimum for encoding both CT and Drug\n    # n_components = 400 # 400 seems to be optimum for encoding both drug only\n    # V33 Ridge50000.0_tsvd400_Dr_LOO_SB1_Tonya_SB0      0.960772 0.67219  0.252819 0.343247 \n    # V31 Ridge200000.0_tsvd250_Dr_CT_LOO_SB1_Tonya_SB0  0.956498 0.668837 0.262345 0.363502\n\n\n\n    list_features_to_encode = ['cell_type', 'sm_name']# [ 'sm_name']  # [ 'sm_name'] # What features to consider - others not used \n    #     # For OP2 task we have two features - cell-type and compound - both categorical\n    #     # Some simple models (like Ridge) better work when one uses only one feature - 'sm_name' (=\"compound\") - so you can try that option. \n    #     # But seems pyboost is powerful enough to proper work with the both features.\n    #     # Table with results for Ridge can be found here: https://www.kaggle.com/code/alexandervc/op2-target-encoders?scriptVersionId=148576642&cellId=1\n    #     # Or here: https://docs.google.com/spreadsheets/d/1APN63PMaWZygVjYimK9Ivt0RvifdAU5JRYkxiDn4szw/edit?usp=sharing (sheet \"Encoders\") and discussed here:\n\n    if 'sm_name' in list_features_to_encode: str_inf_cfg += '_Dr'\n    if 'cell_type' in list_features_to_encode: str_inf_cfg += '_CT'\n\n\n    # smoothing = 10; enc = ce.TargetEncoder(smoothing = smoothing )\n    enc = ce.LeaveOneOutEncoder()\n    str_inf_cfg += '_LOO'\n\n    n_selfblend = 1\n    str_inf_cfg += '_SB'+str(n_selfblend)\n\n\n    CV_scheme = 'Tonya' #  'MT'\n    valid_size = 0 \n    kf = KFold_custom(CV_scheme, valid_size = valid_size)   \n    print('CV_scheme:', CV_scheme)\n    str_inf_cfg += '_'+CV_scheme\n\n    print(str_inf_cfg)\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T08:01:17.654057Z","iopub.execute_input":"2023-11-27T08:01:17.654565Z","iopub.status.idle":"2023-11-27T08:01:17.676503Z","shell.execute_reply.started":"2023-11-27T08:01:17.654519Z","shell.execute_reply":"2023-11-27T08:01:17.675024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.linear_model import Ridge\nfrom sklearn.decomposition import TruncatedSVD\nfrom sklearn.metrics import r2_score\nfrom sklearn.metrics import mean_absolute_error\nimport category_encoders as ce\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T08:01:17.678351Z","iopub.execute_input":"2023-11-27T08:01:17.678744Z","iopub.status.idle":"2023-11-27T08:01:17.689976Z","shell.execute_reply.started":"2023-11-27T08:01:17.678712Z","shell.execute_reply":"2023-11-27T08:01:17.689045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling","metadata":{}},{"cell_type":"code","source":"%%time\nprint('Modeling starts ')\nprint(reducer)\nprint(model)\nprint(enc)\nprint(str_inf_cfg)\n\n\ndf_stat = pd.DataFrame()\n\nfor i_cfg in [2]:#[0,3,4,5]: # [0,1,2,3,4]: # range(1):\n    print('i_cfg', i_cfg)\n    Y_pred_oof_blend = np.zeros( (614, 18211) ); \n    Y_submit = np.zeros( (255, 18211) ); i_blend_submit = 0\n    \n    for i_selfblend in range(n_selfblend ):\n        t0_start_folds_loop = time.time()\n        list_mrrmse = [];list_r2 = []; list_mae = [];  list_corr_row = []; list_corr_col = [];  list_n_train_sample = []\n        for i_fold,(IX_train,IX_test)   in enumerate( kf.split()) :\n            t0_fold = time.time()\n            # \n            if i_cfg == 1:\n                IX_train = np.array( [i for i in IX_train if df_de_train['cell_type'].iat[i] != 'T cells CD8+'  ] )\n            elif i_cfg == 2:\n                np.random.shuffle(IX_train) #  np.array( [i for i in IX_train if df_de_train['cell_type'].iat[i] != 'T cells CD8+'  ] )\n                IX_train = IX_train[:-10]\n            elif i_cfg == 3:\n                IX_train = np.array( list(IX_train) *2 )\n            elif i_cfg == 4:\n                IX_train = np.array( list(IX_train) *4 )\n            elif i_cfg == 5:\n                IX_train = np.array( list(IX_train) *6 )\n            elif i_cfg == 6:\n                IX_train_upd = []\n                dict_ct_mult = {'NK cells':1, 'T cells CD4+':1, 'T cells CD8+':1, 'T regulatory cells':1, 'B cells':3, 'Myeloid cells':3}\n                for ct in dict_ct_mult.keys():            \n                    IX_tmp = [i for i in IX_train if df_de_train['cell_type'].iat[i] == ct ]\n                    IX_train_upd += (IX_tmp * dict_ct_mult[ct] )\n                IX_train = np.array(IX_train_upd )\n            elif i_cfg == 7:\n                IX_train = np.arange(614 )\n                np.random.shuffle(IX_train)\n                IX_train = IX_train[:-3]\n            elif i_cfg == 8:\n                IX_train = np.arange(614 )\n                IX_train = np.array( [i for i in IX_train if df_de_train['cell_type'].iat[i] != 'T cells CD8+'  ] )\n                np.random.shuffle(IX_train)\n                IX_train = IX_train[:-3]\n                \n                \n            # Get train, test:\n            if hasattr(X_full_categorical ,'values'):\n                X_train_categorical,X_test_categorical = X_full_categorical[list_features_to_encode].values[IX_train], X_full_categorical[list_features_to_encode].values[IX_test]\n            else:\n                X_train_categorical,X_test_categorical = X_full_categorical[list_features_to_encode][IX_train], X_full_categorical[list_features_to_encode][IX_test]\n            if hasattr(X_submit_categorical,'values'):\n                X_submit_categorical_loc = X_submit_categorical[list_features_to_encode].values\n            else:\n                X_submit_categorical_loc = X_submit_categorical[list_features_to_encode]\n\n            # Reduce the targets \n            Y_train_red = reducer.fit_transform(Y_full[IX_train] ) \n\n            if verbose >= 100:\n                print('fold:', i_fold, 'X_train_categorical.shape, X_test_categorical.shape, Y_train_red.shape', X_train_categorical.shape, X_test_categorical.shape, Y_train_red.shape)\n\n            ### -------------- Target Encoding -------------------------------------------------------------------\n            for i_target in range(Y_train_red.shape[1]):\n                if i_target == 0:\n                    X_train_encoded = enc.fit_transform(X_train_categorical, Y_train_red[:,i_target])\n                    X_test_encoded = enc.transform(X_test_categorical)\n                    X_submit_encoded = enc.transform(X_submit_categorical_loc)\n                else:\n                    X_encoded_tmp = enc.fit_transform(X_train_categorical, Y_train_red[:,i_target])\n                    X_train_encoded = np.concatenate( [X_train_encoded, X_encoded_tmp], axis = 1)\n                    X_encoded_tmp = enc.transform(X_test_categorical)\n                    X_test_encoded = np.concatenate( [X_test_encoded, X_encoded_tmp], axis = 1)\n                    X_encoded_tmp = enc.transform(X_submit_categorical_loc)\n                    X_submit_encoded = np.concatenate( [X_submit_encoded, X_encoded_tmp], axis = 1)\n\n            list_n_train_sample.append( X_train_encoded.shape[0] ) # Save to inf \n            # Train model\n            model.fit(X_train_encoded, Y_train_red)\n            Y_test_pred_red = model.predict( X_test_encoded )\n            Y_test_pred = reducer.inverse_transform( Y_test_pred_red )\n            # Preparing OOF: \n            #Y_pred_oof_blend[ IX_test ] =  Y_test_pred\n            Y_pred_oof_blend[ IX_test ] = ( Y_pred_oof_blend[ IX_test ] *  i_selfblend + Y_test_pred ) / ( i_selfblend + 1)\n            # Prapare submit\n            Y_submit_pred_red = model.predict( X_submit_encoded )\n            Y_submit_pred = reducer.inverse_transform( Y_submit_pred_red )\n            Y_submit = ( Y_submit *  i_blend_submit  +  Y_submit_pred ) / ( i_blend_submit  + 1  )\n\n            # Scoring  (wiil score blended OOF)\n            Y_test = Y_full[IX_test]\n            Y_test_pred = Y_pred_oof_blend[ IX_test ]\n            mrrmse_one_fold = np.sqrt(np.square(Y_test - Y_test_pred).mean(axis=1)).mean();  r2_one_fold = r2_score( Y_test , Y_test_pred ) \n            mae_one_fold = mean_absolute_error( Y_test , Y_test_pred  )\n            list_mrrmse.append(mrrmse_one_fold); list_r2.append(r2_one_fold)\n            list_mae.append( mae_one_fold )\n            list_tmp =[ np.corrcoef( Y_test[k,:] , Y_test_pred[k,:]  )[0,1] for k in range(Y_test.shape[0]) ]\n            list_corr_row.append( np.mean(list_tmp ))\n            list_tmp_col =[ np.corrcoef( Y_test[:,k] , Y_test_pred[:,k] )[0,1] for k in range(Y_test.shape[1]) ]\n            list_corr_col.append( np.mean(list_tmp_col ))\n            if verbose >= 10:\n                print('fold:', i_fold, 'mrrmse:', np.round(mrrmse_one_fold,4), 'r2:', np.round(r2_one_fold,4), 'mae:', np.round(mae_one_fold,4), \n                      'corr row:', np.round(list_corr_row[-1],4), 'corr col:', np.round(list_corr_col[-1],4), \n                      'Time', np.round(time.time()- t0_fold,1) )\n\n\n\n        IX = len(df_stat)+1\n        df_stat.loc[IX,'model'] = str_inf_cfg +'_SB'+str(i_selfblend) #  str_model_id # 'NLPr_SB'+str(i_selfblend+1)+ '_MUL'+str(n_multiplex_train) + '_Epo'+str(epochs) + '_Noi'+str(prm_GaussianNoise) +\\\n        #     '_RBW' +str(int(restore_best_weights))\n        df_stat.loc[IX,'mrrmse'] = np.mean(list_mrrmse)\n        df_stat.loc[IX,'mae'] = np.mean(list_mae)\n        df_stat.loc[IX,'corr row'] = np.mean(list_corr_row)\n        df_stat.loc[IX,'corr col'] = np.mean(list_corr_col)\n        df_stat.loc[IX,'corr row med'] = np.median(list_corr_row)\n        df_stat.loc[IX,'corr col med'] = np.median(list_corr_col)\n        df_stat.loc[IX,'r2'] = np.mean(list_r2)\n        df_stat.loc[IX,'time'] = np.round(  time.time()- t0_start_folds_loop  ,1 )\n        for kk in range( len( list_mrrmse ) ):\n            df_stat.loc[IX,'mrrmse ' + str(kk)] = list_mrrmse[kk]\n        for kk in range( len( list_mae ) ):\n            df_stat.loc[IX,'mae ' + str(kk)] = list_mae[kk]\n        for kk in range( len( list_corr_row ) ):\n            df_stat.loc[IX,'corr row ' + str(kk)] = list_corr_row[kk]\n        for kk in range( len( list_corr_col ) ):\n            df_stat.loc[IX,'corr col ' + str(kk)] = list_corr_col[kk]    \n        df_stat.loc[IX,'n_feat'] = X_train_encoded.shape[1]\n        df_stat.loc[IX,'n_train'] = np.median( list_n_train_sample )\n\n\n        if verbose >= 1:\n            print('Average mrrmse:', np.round(np.mean(list_mrrmse ),4 ),  'mae:', np.round(np.mean(list_mae ),4 ), \n                  'corr row:',  np.round(np.mean(list_corr_row),4),  'corr col:',  np.round(np.mean(list_corr_col),4) )\n            print('Average r2:', np.round(np.mean(list_r2 ),4 ), 'Time:', np.round(time.time() - t0_start_folds_loop, 1))   \n            print()\n\npd.set_option('display.max_columns', None)\npd.set_option('display.max_rows', None)\n        \ndf_stat.round(6).to_csv('df_stat.csv')\ndisplay(df_stat.round(3) )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Show df_stat","metadata":{}},{"cell_type":"code","source":"%%time\npd.set_option('display.max_columns', None)\npd.set_option('display.max_rows', None)\ndisplay(df_stat  )\n# display(df_stat.sort_values(df_stat.columns[0]  ) )\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T08:01:18.738278Z","iopub.execute_input":"2023-11-27T08:01:18.738940Z","iopub.status.idle":"2023-11-27T08:01:18.751712Z","shell.execute_reply.started":"2023-11-27T08:01:18.738898Z","shell.execute_reply":"2023-11-27T08:01:18.750211Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Save OOF and submit ","metadata":{}},{"cell_type":"code","source":"%%time\nstr_inf4file = str_inf_cfg\nfn4save = 'Y_oof_Index_from1_'+str_inf4file \nfn4save\nnp.save( fn4save+'.npy', Y_pred_oof_blend )\ndf_oof_for_save =  pd.DataFrame( Y_pred_oof_blend )\ndf_oof_for_save.index = range(1,615)\ndf_oof_for_save.index.name = 'id'\ndf_oof_for_save.to_csv(  fn4save + '.csv')\nprint( df_oof_for_save.shape )\n\n\ndf_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\ndf_submit.index.name = 'id'\nprint( df_submit.shape )\ndisplay(df_submit)\nfn4save = 'Y_submit_'+str_inf4file \nprint( fn4save )\nnp.save( fn4save+'.npy', Y_submit )\ndf_submit.to_csv(  fn4save + '.csv')\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T08:01:18.753986Z","iopub.execute_input":"2023-11-27T08:01:18.754452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Final timing","metadata":{}},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\nprint('%.1f minutes passed total '%( (time.time()-t0start)/60)  )\nprint('%.2f hours passed total '%( (time.time()-t0start)/3600)  )","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}