{"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\nCross-validation schemes and their analysis. \n\nCorrect cross-validation seems to be quite importand for current competition, because\npublic LB and private LB are different and that would cause shake-up.\n\nThe main difficulty is that private LB - made on additional DAY(!) (and donor), while public LB ONLY on additional donor. \n\n#### CV with TWO tests. \n\nSo we propose to make a cross-validiation which will take that into account.\n\nThus we propose to have TWO test sets - 1) similar to public LB 2) similar to private LB\n\nWe are lucky that there seems to be natural way to organize it.\nFor CITE-seq it goes as follows:\n\n#### The scheme contains 6-\"folds\" (but warning - it is not completely usual CV scheme ).  \n\n    Fold 0: Train: excludes Day 2 and Donor 32606 Sizes: train: 32536 Test Like Priv 21942 Test Like Publ 16510\n    Fold 1: Train: excludes Day 2 and Donor 31800 Sizes: train: 32638 Test Like Priv 21942 Test Like Publ 16408\n    Fold 2: Train: excludes Day 3 and Donor 32606 Sizes: train: 33100 Test Like Priv 20901 Test Like Publ 16987\n    Fold 3: Train: excludes Day 3 and Donor 31800 Sizes: train: 31543 Test Like Priv 20901 Test Like Publ 18544\n    Fold 4: Train: excludes Day 4 and Donor 32606 Sizes: train: 28368 Test Like Priv 28145 Test Like Publ 14475\n    Fold 5: Train: excludes Day 4 and Donor 31800 Sizes: train: 28189 Test Like Priv 28145 Test Like Publ 14654\n\nEach \"Fold\" is parametrized by excluded from train day and donor. And folds are like that: \n\n    train_index = np.where( (df_meta['day']  != day2exclude) & ( df_meta['donor']  != donor2exclude  ) )  [0]\n    test_index1_like_private_lb = np.where( df_meta['day']  == day2exclude)[0]\n    test_index2_like_public_lb = np.where( (df_meta['day']  != day2exclude) &  (df_meta['donor']  == donor2exclude ) ) [0]\n\n#### Example:  assume day2exclude = 2, donor2exlude = 32606,\n     We train on days 3,4 and without donor 32606,\n     as a private test we consider day 2 , INCLUDING donor 32606 - thus it is quite similar to real private LB - where we have unseen day and unseen donor\n     as a public test we consider days 3,4 donor is ONLY 32606 - thus it is similar to real public LB - the days are the same, but donor is unseen\n     \n#### Male/Female are taken into account: \n    We train always on donors where we have male and female. \n    And predict only for MALE.\n    That is exactly as on real LB. \n    Thus we have 6 folds - we only exclude from train male donors 32606,  31800 - thus MALE is always in test - like on real LB.\n    \n### Update. From version 4 we add advanced option - additional holdout folds for advanced control (10%).\n\n    Such scheme have same 6 folds, but each has 4 tests - see details below. \n    It is an option - one may not use it  - both schemes are presereved in the code. \n  \n\n### Versions:\n\n#### 14 try PCA for main features instead of tSVD (100 also)\n    \n    citeseq_train_and_test_PCA500.csv\n    \n#### 13 Started to play with adding features - adding meta + Kmeans - improve both CV&LB\n\n#### 12 Started to play with adding features - meta data: cell type, gender - some small improve on CV and LB\n\n    Improved LB\n    \n    \n#### 11 Played with number of dimensions in tSVD - seems 100 is optimal, both 75 and 150 are a bit worse on CV difference is 0.000055 - no improvement on LB ! 0.805 but not the best. \n\n     ! no improvement on LB !  0.805 but not the best. \n    \n\n#### 10 Played with Alpha for Rigde seems 1e4 is optimal on CV, but improve is very small 0.0001 ,  LB: 0.805 \n    \n    Indeed LB is improved 0.805 V10\n    But we do not know how much exactly - as expected\n   \n#### 9 Option to rescale Y to mean 0, std 1  LB: 0.805 \n\n    if rescale_Y_to_mean0_std1:\n        Y -= Y.mean(axis=1).reshape(-1, 1)\n        Y /= Y.std(axis=1).reshape(-1, 1)\n        print('Rescaling to mean 0 and std 1 has been done')\n        \n    It improves Ridge's results\n    \n    CV corresponds to LB.\n    It seems \"public\" like test and \"private\" like test are quite similar. So may be no shake-up. \n    \n\n#### 8 Added multiome part predictions from sskknt msci-citeseq-keras-quickstart-dropout LB: 0.803\n\n    https://www.kaggle.com/code/sskknt/msci-citeseq-keras-quickstart-dropout\n\n#### 7 Added code to prepare submissions\n\n    Multome part was submitted as ZEROs all. \n\n#### 6 Slight modifications - now use engineered features from the dataset \"Feature shop...\"\n    \n    Only small updates in code. \n    \n\n#### 5 Added new - correlation metric scoring - the same as used on LB \n\n    previously only R2 has been used\n\n#### 4 Added new CV-scheme - additionally created 2 holdouts - like private and like public LB\n    \n    Can be used to additional testing - conrol that we are not overfitting on our CV,\n    and for futher possible blend of the models. \n    Holdouts are about 10% of original tests. \n    \n    Note: that does NOT decrease size of train sets. Only tests are now splitted. \n\n#### 3 simple modeling (Ridge) example with TruncatedSVD is added\n\n    Strange conclusion - the model on private-like-test is better, than on  public-like-test\n\n#### 2 - 6-fold validation scheme with TWO TESTs is described\n\n    \n    CITE-seq part only is considered.\n    \n#### 1  - Draft\n    \n\n\n\n","metadata":{}},{"cell_type":"code","source":"rescale_Y_to_mean0_std1 = False","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:08.104179Z","iopub.execute_input":"2022-10-31T19:48:08.104623Z","iopub.status.idle":"2022-10-31T19:48:08.110184Z","shell.execute_reply.started":"2022-10-31T19:48:08.104537Z","shell.execute_reply":"2022-10-31T19:48:08.108849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Install/import","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 numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\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":"2022-10-31T19:48:08.135460Z","iopub.execute_input":"2022-10-31T19:48:08.136124Z","iopub.status.idle":"2022-10-31T19:48:08.187047Z","shell.execute_reply.started":"2022-10-31T19:48:08.136081Z","shell.execute_reply":"2022-10-31T19:48:08.186310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nimport 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-31T19:48:08.188719Z","iopub.execute_input":"2022-10-31T19:48:08.189240Z","iopub.status.idle":"2022-10-31T19:48:30.170553Z","shell.execute_reply.started":"2022-10-31T19:48:08.189205Z","shell.execute_reply":"2022-10-31T19:48:30.169544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def correlation_score(y_true, y_pred):\n    \"\"\"Scores the predictions according to the competition rules. \n    \n    It is assumed that the predictions are not constant.\n    \n    Returns the average of each sample's Pearson correlation coefficient\"\"\"\n    if type(y_true) == pd.DataFrame: y_true = y_true.values\n    if type(y_pred) == pd.DataFrame: y_pred = y_pred.values\n    corrsum = 0\n    for i in range(len(y_true)):\n        corrsum += np.corrcoef(y_true[i], y_pred[i])[1, 0]\n    return corrsum / len(y_true)","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:30.172321Z","iopub.execute_input":"2022-10-31T19:48:30.172663Z","iopub.status.idle":"2022-10-31T19:48:30.179473Z","shell.execute_reply.started":"2022-10-31T19:48:30.172632Z","shell.execute_reply":"2022-10-31T19:48:30.177810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data\n","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\n\nif 0: # Notebook version from V6 will use genereated features, not the original data \n    # Here df_cite contains ONLY train part - otherwise it will crash RAM \n    print('Load original data')\n    \n    #df_cite_train_x = pd.read_hdf(FP_CITE_TRAIN_INPUTS)\n    df_cite = pd.read_hdf(FP_CITE_TRAIN_INPUTS)\n    print('8.2G RAM consumed')\n    display(df_cite.head())    \n    \n    filename = '/kaggle/input/open-problems-multimodal/test_cite_inputs.h5'\n\n    f2 = h5py.File(filename,'r')#, mode)\n    print(f2.keys())\n    first_key = 'test_cite_inputs'\n    print(f2[first_key].keys())\n    list_cells_barcodes = [t.decode() for t in f2[first_key]['axis1'] ]\n    print(len(list_cells_barcodes), list_cells_barcodes[:5], list_cells_barcodes[-5:])\n\n","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:30.180662Z","iopub.execute_input":"2022-10-31T19:48:30.180977Z","iopub.status.idle":"2022-10-31T19:48:45.843306Z","shell.execute_reply.started":"2022-10-31T19:48:30.180942Z","shell.execute_reply":"2022-10-31T19:48:45.841588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif 1:\n    print('Load CITE-seq targets and ')\n    df_cite_train_y = pd.read_hdf(FP_CITE_TRAIN_TARGETS)\n    display(df_cite_train_y)","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:45.847555Z","iopub.execute_input":"2022-10-31T19:48:45.848689Z","iopub.status.idle":"2022-10-31T19:48:46.602319Z","shell.execute_reply.started":"2022-10-31T19:48:45.848640Z","shell.execute_reply":"2022-10-31T19:48:46.600834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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)\nif 1:\n    df_meta = pd.DataFrame(index = df_cite_train_y.index) \n    df_meta = df_meta.join(df_cell.set_index('cell_id') )\n    display(df_meta)","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:46.603841Z","iopub.execute_input":"2022-10-31T19:48:46.604604Z","iopub.status.idle":"2022-10-31T19:48:46.913466Z","shell.execute_reply.started":"2022-10-31T19:48:46.604567Z","shell.execute_reply":"2022-10-31T19:48:46.912052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta['donor'].unique()","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:46.915090Z","iopub.execute_input":"2022-10-31T19:48:46.915465Z","iopub.status.idle":"2022-10-31T19:48:46.925787Z","shell.execute_reply.started":"2022-10-31T19:48:46.915434Z","shell.execute_reply":"2022-10-31T19:48:46.924179Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create Folds for Cross Validation \n\nWe will create several example of cross validation schemes.\nThe information will be stored like that:\n\n    Example with 2 folds; and train and test:\n    list_folds_indices = [ (train_index1, test_index1 ), (train_index2, test_index2  )  ]\n\n    Example with 3 folds;  and train and  TWO tests:\n    list_folds_indices = [ (train_index1, test_index1a,  test_index1b ), (train_index2, test_index2a,  test_index2b ), (train_index3, test_index3a,  test_index3b )  ]\n","metadata":{}},{"cell_type":"markdown","source":"## Simple illustrative examples","metadata":{}},{"cell_type":"code","source":"# Create indices \"train_index\", \"test_index\" \nmask = ( (df_meta['day'] == 2) |  (df_meta['day'] == 3) )\ntrain_index = np.where(mask > 0 )[0]\ntest_index = np.where((~mask) > 0 )[0]\n\n# Store  as 1-fold scheme \n# Storage convention - first in tuple - is always a \nlist_folds_indices = [ (train_index, test_index )  ]\n\n# Store \"2-fold\" scheme\ntrain_index2 , test_index2  = test_index , train_index\nlist_folds_indices2 = [ (train_index, test_index ), (train_index2, test_index2  )  ]\n","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:46.927634Z","iopub.execute_input":"2022-10-31T19:48:46.928054Z","iopub.status.idle":"2022-10-31T19:48:46.939577Z","shell.execute_reply.started":"2022-10-31T19:48:46.928017Z","shell.execute_reply":"2022-10-31T19:48:46.938009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check:\nprint( train_index.shape, test_index.shape, np.intersect1d(test_index, train_index ) )","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:46.942473Z","iopub.execute_input":"2022-10-31T19:48:46.942933Z","iopub.status.idle":"2022-10-31T19:48:46.953957Z","shell.execute_reply.started":"2022-10-31T19:48:46.942899Z","shell.execute_reply":"2022-10-31T19:48:46.952609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Folds corresponding to days - relevant to private LB (but not perfect). (BAD for public LB)  ","metadata":{}},{"cell_type":"code","source":"list_folds_indices_by_days = []\nfor day2exclude in [2,3,4]:\n    train_index = np.where( df_meta['day']  != day2exclude)[0]\n    test_index = np.where( df_meta['day']  == day2exclude)[0]\n    list_folds_indices_by_days.append( (train_index,  test_index ))\n\nlen(list_folds_indices_by_days)","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:46.957047Z","iopub.execute_input":"2022-10-31T19:48:46.957673Z","iopub.status.idle":"2022-10-31T19:48:46.972129Z","shell.execute_reply.started":"2022-10-31T19:48:46.957627Z","shell.execute_reply":"2022-10-31T19:48:46.970712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check:\nprint( train_index.shape, test_index.shape, np.intersect1d(test_index, train_index ) )","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:46.977968Z","iopub.execute_input":"2022-10-31T19:48:46.978756Z","iopub.status.idle":"2022-10-31T19:48:46.988551Z","shell.execute_reply.started":"2022-10-31T19:48:46.978718Z","shell.execute_reply":"2022-10-31T19:48:46.985711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Should be good: 6-fold scheme with two tests - for private and public LB","metadata":{}},{"cell_type":"code","source":"list_folds_indices_by_days_and_donors = []\nc = 0\nfor day2exclude in [2,3,4]:\n    for donor2exclude in [32606,  31800]: # We will need to predict always MALE (not female) - like on LB. (# donor 13176 - female)\n        train_index = np.where( (df_meta['day']  != day2exclude) & ( df_meta['donor']  != donor2exclude  ) )  [0]\n        test_index1_like_private_lb = np.where( df_meta['day']  == day2exclude)[0]\n        test_index2_like_public_lb = np.where( (df_meta['day']  != day2exclude) &  (df_meta['donor']  == donor2exclude ) ) [0]\n        list_folds_indices_by_days_and_donors.append( (train_index,  test_index1_like_private_lb , test_index2_like_public_lb) )\n    \n        str_fold_inf = 'Fold ' +str(c) + ': Train: excludes Day '+str(day2exclude) + ' and Donor ' + str( donor2exclude )\n        print(str_fold_inf, 'Sizes: train:',len(train_index), 'Test Like Priv'  ,len(test_index1_like_private_lb), 'Test Like Publ',   len(test_index2_like_public_lb),  ); c+=1\n    \nprint(len(list_folds_indices_by_days_and_donors))","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:46.990110Z","iopub.execute_input":"2022-10-31T19:48:46.990910Z","iopub.status.idle":"2022-10-31T19:48:47.011791Z","shell.execute_reply.started":"2022-10-31T19:48:46.990853Z","shell.execute_reply":"2022-10-31T19:48:47.010360Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Advanced CV - 6-folds like above with some additional hold-out \n\nTo be done \n\nSmall additional hold out - might be used for 1) check we not overting for current CV 2) to blend models \n\nSo we can just holdout small random part from test parts of the previous folds - something \n\nThus scheme have THREE tests - for public LB, for private LB, and hold out for additional internal validation / blending\n","metadata":{}},{"cell_type":"code","source":"df_meta.columns","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.013250Z","iopub.execute_input":"2022-10-31T19:48:47.014040Z","iopub.status.idle":"2022-10-31T19:48:47.022669Z","shell.execute_reply.started":"2022-10-31T19:48:47.014006Z","shell.execute_reply":"2022-10-31T19:48:47.021353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Preparations for folds creation","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 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']\ndf_meta\n\ny = df_meta[scol] \nskf = StratifiedKFold(n_splits=10,  shuffle=True, random_state=40)\nskf.get_n_splits(df_meta, 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) ); break\nprint(test_index)\nprint(df_meta[scol].value_counts()   )\nprint(df_meta.iloc[test_index,:][scol].value_counts()    )\n\ndf_meta['HoldOut'] = 0 \ndf_meta.loc[df_meta.index[test_index],'HoldOut']  = 1\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.024284Z","iopub.execute_input":"2022-10-31T19:48:47.024724Z","iopub.status.idle":"2022-10-31T19:48:47.211935Z","shell.execute_reply.started":"2022-10-31T19:48:47.024688Z","shell.execute_reply":"2022-10-31T19:48:47.210734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_folds_indices_by_days_and_donors_with2holdouts = []\nc = 0\nfor day2exclude in [2,3,4]:\n    for donor2exclude in [32606,  31800]: # We will need to predict always MALE (not female) - like on LB. (# donor 13176 - female)\n        train_index = np.where( (df_meta['day']  != day2exclude) & ( df_meta['donor']  != donor2exclude  ) )  [0]\n        mask_holdout = (df_meta['HoldOut']==1)\n        test_index1_like_private_lb = np.where( (df_meta['day']  == day2exclude) & (~mask_holdout) ) [0]\n        test_index1_like_private_lb_holdout = np.where( (df_meta['day']  == day2exclude) & (mask_holdout) ) [0]\n        test_index2_like_public_lb = np.where( (df_meta['day']  != day2exclude) &  (df_meta['donor']  == donor2exclude ) & (~mask_holdout) ) [0]\n        test_index2_like_public_lb_holdout = np.where( (df_meta['day']  != day2exclude) &  (df_meta['donor']  == donor2exclude ) & mask_holdout ) [0]\n        \n        list_folds_indices_by_days_and_donors_with2holdouts.append( (train_index,  test_index1_like_private_lb , test_index2_like_public_lb, \n                                                      test_index1_like_private_lb_holdout, test_index2_like_public_lb_holdout ) )\n    \n        str_fold_inf = 'Fold ' +str(c) + ': Train: excludes Day '+str(day2exclude) + ' and Donor ' + str( donor2exclude )\n        print(str_fold_inf, 'Sizes: train:',len(train_index), 'Test Like Priv'  ,len(test_index1_like_private_lb), \n              'Test Like Publ',   len(test_index2_like_public_lb),  \n              'Test Like Priv HoldOut',   len(test_index1_like_private_lb_holdout),  \n              'Test Like Publ HoldOut',   len(test_index2_like_public_lb_holdout),  \n              \n             ); c+=1\n    \nprint(len(list_folds_indices_by_days_and_donors_with2holdouts))","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.213737Z","iopub.execute_input":"2022-10-31T19:48:47.214870Z","iopub.status.idle":"2022-10-31T19:48:47.242806Z","shell.execute_reply.started":"2022-10-31T19:48:47.214816Z","shell.execute_reply":"2022-10-31T19:48:47.241192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Optional Rescaling of Targets","metadata":{}},{"cell_type":"code","source":"# Since metric is - correlation coefficient - any rescaling aY+b will not be change it , so one can do like that:\n# it is not clear is it optimal or not (see https://www.kaggle.com/competitions/open-problems-multimodal/discussion/360253 )\n\nY = df_cite_train_y.values\n\nif rescale_Y_to_mean0_std1:\n    Y -= Y.mean(axis=1).reshape(-1, 1)\n    Y /= Y.std(axis=1).reshape(-1, 1)\n    print('Rescaling to mean 0 and std 1 has been done')\n    ","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.244588Z","iopub.execute_input":"2022-10-31T19:48:47.244949Z","iopub.status.idle":"2022-10-31T19:48:47.250065Z","shell.execute_reply.started":"2022-10-31T19:48:47.244920Z","shell.execute_reply":"2022-10-31T19:48:47.249133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling","metadata":{}},{"cell_type":"markdown","source":"# Simple modeling example :","metadata":{}},{"cell_type":"code","source":"# %%time\n# from sklearn.metrics import r2_score\n\n# N_features = 1 # 200 # For the fast test - just take some small number of features \n# # When number of features is very low like 1 - the results are more or less the same as just MEAN over the targets - kind of baseline solution (no modeling in fact)\n\n# cc = 0\n# for list_folds_indices in  [list_folds_indices_by_days_and_donors, list_folds_indices_by_days_and_donors_with2holdouts]:\n#     if cc == 0: print('CV withOUT holdout\\n' ); \n#     else : print('\\n\\nCV with holdout\\n' )\n#     cc += 1\n    \n    \n#     alpha4Ridge  = 1\n\n#     X = df_cite.iloc[:,:N_features].values\n#     print('X.shape:', X.shape,'Y.shape:', Y.shape)\n\n#     from sklearn.linear_model import Ridge\n#     model = Ridge(alpha=alpha4Ridge)\n#     #, (train_index, test_index )\n\n#     df_fold_score_stat = pd.DataFrame();df_fold_score_stat.index.name = 'Fold'\n#     t0 = time.time()\n#     for fold, indices_tuple  in enumerate( list_folds_indices ):\n#         train_index = indices_tuple[0]\n#         main_test_index = indices_tuple[1]\n#         print('Fold:', fold, 'Shapes of train:', train_index.shape, 'Tests: ',[t.shape for t in indices_tuple[1:] ] )\n#         t1 = time.time()\n\n#         # Train model: \n#         model.fit(X[train_index], Y[train_index])\n\n#         # Calculate metrics on test and train folds \n#         list_scores = []; list_scores_r2 = []\n#         for i_loc in range(0,len( indices_tuple)  ):\n#             indices_loc = indices_tuple[i_loc]\n#             y_pred = model.predict( X[indices_loc ])\n#             #s = model.score( X[indices_loc ], Y[indices_loc]  ) \n#             s = correlation_score( y_pred , Y[indices_loc]  )\n#             list_scores.append(s)\n#             s = r2_score( Y[indices_loc] , y_pred  )\n#             list_scores_r2.append(s)\n\n#         # Just save statistic for output\n#         for i_loc in range(1, len(list_scores  )):\n#              df_fold_score_stat.loc[fold, 'Score Test'+str(i_loc)] = list_scores[i_loc]\n#         df_fold_score_stat.loc[fold, 'Score Train'] = list_scores[0]\n#         df_fold_score_stat.loc[fold, 'Time'] = np.round( (time.time() - t1), 2) \n#         for i_loc in range(1, len(list_scores  )):\n#              df_fold_score_stat.loc[fold, 'R2 Score Test'+str(i_loc)] = list_scores_r2[i_loc]\n#         df_fold_score_stat.loc[fold, 'R2 Score Train'] = list_scores_r2[0]\n\n#         print('Correlation scores:', np.round(list_scores,4),  'time:', '%.2f'%(time.time() - t1))\n        \n#     display(df_fold_score_stat)\n#     display(df_fold_score_stat.describe(percentiles=[] ).iloc[1:,:])\n#     #display(df_fold_score_stat.mean())","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.251318Z","iopub.execute_input":"2022-10-31T19:48:47.252035Z","iopub.status.idle":"2022-10-31T19:48:47.268897Z","shell.execute_reply.started":"2022-10-31T19:48:47.252007Z","shell.execute_reply":"2022-10-31T19:48:47.267783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nStratificationSplitter class\n\"\"\"\nimport math\n\nclass StratificationSplitter:\n    def __init__(self, data, stratification_column_names, max_categories=10,\n                 test_size=None, val_size=None, n_train_cv_splits=3, random_state=None):\n    \n        print(\n            \"Initializing StratificationSplitter with test_size {}, val_size {} and {} train data cross-validation \"\n            \"splits. Stratification by {} columns with max_categories {}\".format(\n                test_size if test_size else 0, val_size if val_size else 0, n_train_cv_splits,\n                stratification_column_names, max_categories\n        ))\n        assert data.shape[0] > n_train_cv_splits\n        self._data = data[~data.index.duplicated(keep=\"last\")]\n        self.stratification_column_names = stratification_column_names\n        assert max_categories > 2\n        self.max_categories = max_categories\n        self.test_size = test_size\n        self.val_size = val_size\n        if test_size is not None:\n            self.test_size = math.ceil(self._data.shape[0] * test_size) if test_size < 1 else int(test_size)\n        if val_size is not None:\n            self.val_size = math.ceil(self._data.shape[0] * val_size) if val_size < 1 else int(val_size)\n        assert self.test_size is None or 0 < self.test_size <= self._data.shape[0] - n_train_cv_splits\n        assert self.val_size is None or 0 < self.val_size <= self._data.shape[0] - n_train_cv_splits\n        if test_size is not None and val_size is not None:\n            if self.test_size + self.val_size > data.shape[0] - n_train_cv_splits:\n                raise ValueError(\n                    \"Sum of test and validation sizes must be less than dataset size minus `n_train_cv_splits`\"\n                )\n        assert n_train_cv_splits >= 2\n        self.n_train_cv_splits = n_train_cv_splits\n        self.random_state = random_state\n        self._stratification = pd.DataFrame([], index=self._data.index.values)\n        self._train_ids = None\n        self._test_ids = None\n        self._val_ids = None\n        self._cv_ids = None\n        for col in stratification_column_names:\n            if col not in self._data.columns:\n                raise ValueError(\"Column `{}` wasn't found in dataset\".format(col))\n            self._add_stratification_column(self._data[col])\n        np.random.seed(self.random_state)\n        if self._stratification.empty:\n            random_col = np.zeros_like(self._stratification.index.values, dtype=np.int64)\n            random_col[:len(random_col) // 2] = 1\n            np.random.shuffle(random_col)\n            self._stratification[\"no_stratification (random)\"] = random_col\n\n    def _add_stratification_column(self, col):\n        unique, unique_counts = np.unique(col, return_counts=True)\n        if len(unique) <= self.max_categories and np.all(unique_counts >= self.n_train_cv_splits):\n            self._stratification[col.name] = col\n        elif len(unique) > self.max_categories and col.dtype.kind in {\"f\", \"u\", \"i\"}:\n            self._stratification[col.name] = pd.qcut(col, self.max_categories, duplicates=\"drop\")\n            buckets = list(zip(*np.unique(self._stratification[col.name], return_counts=True)))\n            print(\n                \"Column `{}` is numeric with more than {} unique values, so it was quantile-discretized \"\n                \"into {} buckets: {}\".format(col.name, self.max_categories, len(buckets), buckets)\n            )\n        elif ~np.all(unique_counts >= self.n_train_cv_splits):\n            print(\n                \"Column `{}` was removed from stratification because it had categories with less than {} members. The \"\n                \"minimum number of members in any category cannot be less than number of cross-validation \"\n                \"splits\".format(col.name, self.n_train_cv_splits)\n            )\n        else:\n            print(\n                \"Column `{}` was removed from stratification because it had more than {} categories and couldn't be \"\n                \"discretized\".format(col.name, self.max_categories)\n            )\n\n    def _split(self):\n        # train/test/validation split\n#         print(\"Creating train/test/validation set splits...\")\n        train_ids = self._stratification.index.values\n        all_train_test_val_ids = [self._stratification.index.values]\n        all_train_test_val_ids_names = [\"full dataset\"]\n        if self.test_size is not None:\n            train_ids, test_ids = train_test_split(\n                train_ids, stratify=self._stratification.loc[train_ids],\n                test_size=self.test_size, random_state=self.random_state\n            )\n            self._test_ids = test_ids\n            all_train_test_val_ids.append(test_ids)\n            all_train_test_val_ids_names.append(\"test\")\n        if self.val_size is not None:\n            train_ids, val_ids = train_test_split(\n                train_ids, stratify=self._stratification.loc[train_ids],\n                test_size=self.val_size, random_state=self.random_state\n            )\n            self._val_ids = val_ids\n            all_train_test_val_ids.append(val_ids)\n            all_train_test_val_ids_names.append(\"validation\")\n        self._train_ids = train_ids\n        all_train_test_val_ids.append(train_ids)\n        all_train_test_val_ids_names.append(\"train\")\n        # log stratification percentages\n#         for feature in self._stratification.columns:\n#             messages = list()\n#             messages.append(\"Stratification on feature: %s\" % feature)\n#             for i, ids in enumerate(all_train_test_val_ids):\n#                 messages.append(\n#                     \"%s sample:  Shape: %s, statistics: %s\" % (\n#                         all_train_test_val_ids_names[i].upper(), len(ids),\n#                         (self._stratification.loc[ids, feature].value_counts().sort_index() / len(ids)).to_dict()\n#                     )\n#                 )\n#             print(\"\\n\".join(messages))\n\n        # train cross validation split\n        print(\"Creating KFold cross-validation train data splits...\")\n        cv = []\n        np.random.seed(self.random_state)\n        groups = self._stratification.loc[train_ids].groupby(list(self._stratification.columns)).groups.items()\n        reminder = 0\n        for key, idx in sorted(groups, key=lambda x: x[0]):\n            group_arr = idx.values.copy()\n            min_examples = int(len(group_arr) / self.n_train_cv_splits)\n            if min_examples < 1:\n                print(\n                    \"Combined category `{}` of columns {} cannot be split into cross-validation folds equally because \"\n                    \"the minimum number of members in any category cannot be less than number of cross-validation \"\n                    \"splits: {} < {}\".format(\n                        key, list(self._stratification.columns), len(group_arr), self.n_train_cv_splits\n                    )\n                )\n            np.random.shuffle(group_arr)\n            split = np.array_split(group_arr, self.n_train_cv_splits)\n            cv.append(split[self.n_train_cv_splits - reminder:] + split[0:self.n_train_cv_splits - reminder])\n            reminder = (reminder + len(group_arr)) % self.n_train_cv_splits\n            # log stratification percentages\n#             message = (\n#                 \"KFold cross-validation stratification on combined category `{}` of columns {} split \"\n#                 \"sizes: {}\".format(key, list(self._stratification.columns), list(map(len, split)))\n#             )\n#             print(message)\n        cv = list(map(np.concatenate, zip(*cv)))\n        fold_sizes = list(map(len, cv))\n        assert sum(fold_sizes) == len(train_ids)\n        message = \"KFold cross-validation split sizes: {}\".format(fold_sizes)\n        print(message)\n        self._cv_ids = []\n        for idx in range(len(cv)):\n            val = cv[idx]\n            train = np.concatenate(cv[0:idx] + cv[idx+1:])\n            assert len(set(train).intersection(set(val))) == 0, \"No such indices that are in both train and val sets\"\n            self._cv_ids.append((train, val))\n\n    @property\n    def train_ids(self):\n        if self._train_ids is None:\n            self._split()\n        return self._train_ids\n\n    @property\n    def test_ids(self):\n        if self._test_ids is None and self.test_size is not None:\n            self._split()\n        return self._test_ids\n\n    @property\n    def val_ids(self):\n        if self._val_ids is None and self.val_size is not None:\n            self._split()\n        return self._val_ids\n\n    @property\n    def cv_ids(self):\n        if self._cv_ids is None:\n            self._split()\n        return self._cv_ids\n    \n\n\ndef get_data_by_ids(features, ids, target_col):\n    samples = features.loc[ids, :]\n    y = samples[target_col].values\n    x = samples.drop(columns=[target_col])\n    return x, y\n\ndef get_indices_by_ids(features, ids):\n    \"\"\"\n    Retrieves indices of features for specified ids\n    :param features:                    DataFrame of all features with id index\n    :param ids:                         list of ids to retrieve indices for\n    :return:                            array of indices for requested ids\n    \"\"\"\n    features[\"indexing_column\"] = np.arange(features.shape[0], dtype=np.int64)\n    samples = features.loc[ids, :]\n    indices = samples[\"indexing_column\"].values\n    features.drop(columns=[\"indexing_column\"], inplace=True)\n    return indices","metadata":{"execution":{"iopub.status.busy":"2022-10-31T21:43:18.345433Z","iopub.execute_input":"2022-10-31T21:43:18.345750Z","iopub.status.idle":"2022-10-31T21:43:18.381606Z","shell.execute_reply.started":"2022-10-31T21:43:18.345724Z","shell.execute_reply":"2022-10-31T21:43:18.379813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nFunctions for hyper-parameter optimization\n\"\"\"\nfrom scipy import stats\nfrom sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score, make_scorer, mean_gamma_deviance\nfrom sklearn.model_selection import train_test_split, GridSearchCV, RandomizedSearchCV\nfrom sklearn import clone\nfrom functools import reduce\nfrom operator import mul\nfrom lightgbm import early_stopping\nimport eli5\nfrom eli5.permutation_importance import get_score_importances\nfrom eli5.sklearn import PermutationImportance\n\n\ndef format_time(seconds):\n    \"\"\"\n    Format time in seconds to time string, including minutes and hours when appropriate\n    :param seconds:                     float, seconds\n    :return:                            formatted time string\n    \"\"\"\n    if seconds < 60:\n        return \"0:{:0>2}\".format(int(seconds))\n    elif seconds < 3600:\n        minutes = int(seconds / 60)\n        seconds = int(seconds % 60)\n        return \"{}:{:0>2}\".format(minutes, seconds)\n    else:\n        hours = int(seconds / 3600)\n        minutes = int((seconds % 3600) / 60)\n        seconds = int((seconds % 3600) % 60)\n        return \"{}:{:0>2}:{:0>2}\".format(hours, minutes, seconds)\n\ndef tune_hyperparameters_cv(estimator, x, y, search_spaces, optimizers, n_iter, scoring=None,\n                            cv=3, n_jobs=1, \n                            verbose_output=None, \n                            random_state=None, sample_weight=None, **kwargs):\n    allowed_optimizers = {\"grid\", \"random\"}\n    invalid_optimizers = [optimizer for optimizer in optimizers if optimizer not in allowed_optimizers]\n    if len(invalid_optimizers) > 0:\n        raise ValueError(\n            \"Values in `optimizers` must be one of {}, found: {}\".format(allowed_optimizers, invalid_optimizers)\n        )\n    if len(search_spaces) != len(optimizers):\n        raise ValueError(\n            \"`search_spaces` and `optimizers` must have the same length, found: \"\n            \"{} != {}\".format(len(search_spaces), len(optimizers))\n        )\n    n_steps = len(search_spaces)\n    cv_splits = cv if isinstance(cv, int) else len(cv)\n    processed_best_params = {}\n    processed_best_score = -float(\"inf\")\n    for step in range(n_steps):\n        step_start_time = time.time()\n        print(\n            \"Step {} of {} of hyper-parameter optimization \"\n            \"(`{}` optimizer)\".format(step + 1, n_steps, optimizers[step])\n        )\n        if optimizers[step] == \"grid\":\n            n_fits = reduce(mul, [len(par_values) for par_values in search_spaces[step].values()]) * cv_splits\n            optimizer_kwargs = {}\n            optimizer = GridSearchCV\n        else:\n            n_fits = n_iter * cv_splits\n            optimizer_kwargs = {\"random_state\": random_state, \"n_iter\": n_iter}\n            optimizer = RandomizedSearchCV\n        print(\"Models fitting time start: {}\".format(time.strftime('%H:%M:%S', time.gmtime(time.time()))))\n        print(\"Optimizer is fitting {} models...\".format(n_fits))\n        estimator_clone = clone(estimator)\n        estimator_clone.set_params(**processed_best_params)\n        model = optimizer(\n            estimator_clone, search_spaces[step], scoring=scoring, cv=cv,\n            n_jobs=n_jobs, refit=False, **optimizer_kwargs\n        )\n        model.fit(x, y, **kwargs, callbacks=[early_stopping(stopping_rounds=100, first_metric_only=False, verbose=False)])\n        iter_params = model.cv_results_[\"params\"]\n        iter_scores = model.cv_results_[\"mean_test_score\"]\n        iter_score_stds = model.cv_results_[\"std_test_score\"]\n#         for iter_idx in range(len(iter_params)):\n#             print(\n#                 \"Parameters {}: mean-score={:.12f}, \"\n#                 \"std-score={:.12f}\".format(iter_params[iter_idx], iter_scores[iter_idx], iter_score_stds[iter_idx])\n#             )\n        print(\n            \"Step {} optimization time: {}\".format(step + 1, format_time(time.time() - step_start_time))\n        )\n        print(\n            \"Step {} completed. Best parameter tuning score is {:.12f} with \"\n            \"parameters: {}\".format(step + 1, model.best_score_, model.best_params_)\n        )\n        if model.best_score_ >= processed_best_score:\n            processed_best_score = model.best_score_\n            processed_best_params.update(model.best_params_)\n        else:\n            print(\n                \"Step {} best tuning score {:.12f} is worse than the previous step score {:.12f}, so parameters \"\n                \"will be discarded and the correspondent default parameter values from the previous step will \"\n                \"be used instead\".format(step + 1, model.best_score_, processed_best_score)\n            )\n    print(\n        \"Best parameter tuning score overall is {:.12f} with \"\n        \"parameters: {}\".format(processed_best_score, processed_best_params)\n    )\n    print(\"Refitting model on the whole train dataset with best parameters..\")\n    estimator.set_params(**processed_best_params)\n    estimator.fit(x, y, **kwargs, callbacks=[early_stopping(stopping_rounds=100, first_metric_only=False, verbose=False)])\n    print(\"Hyper-parameter tuning completed\")\n    print(processed_best_params)\n    \n    perm = PermutationImportance(estimator,scoring=None, n_iter=1, random_state=42, cv=None, refit=False).fit(x,y)\n    tmp = eli5.show_weights(perm)\n    display(eli5.show_weights(perm, top = len(list(x.columns)), feature_names = list(x.columns)))\n\n    return estimator, processed_best_params","metadata":{"execution":{"iopub.status.busy":"2022-10-31T20:30:14.832747Z","iopub.execute_input":"2022-10-31T20:30:14.833148Z","iopub.status.idle":"2022-10-31T20:30:14.858482Z","shell.execute_reply.started":"2022-10-31T20:30:14.833113Z","shell.execute_reply":"2022-10-31T20:30:14.856742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-31T19:48:47.335729Z","iopub.execute_input":"2022-10-31T19:48:47.336213Z","iopub.status.idle":"2022-10-31T19:48:47.517515Z","shell.execute_reply.started":"2022-10-31T19:48:47.336126Z","shell.execute_reply":"2022-10-31T19:48:47.516301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import lightgbm as lgbm\nfrom sklearn.multioutput import MultiOutputClassifier, MultiOutputRegressor\n\n\nLGB_PARAMETERS = {\n#     \"random_state\": 0, #+\n    \"learning_rate\": 0.1 ,  #+ hypersensitive - change to 0.0671 - worsens about 2% from 0.43590181064083766\n    \"num_leaves\": 12, #+ hypersensitive - change to 11, worsens about 2%\n    \"max_depth\": 9, #+ hypersensitive - change to 10 - worsens about 1%\n    \"n_estimators\": 300, #+ senstitive - change by 10 - worsents about 0.5%\n    \"min_child_samples\": 19,#+ hypersensitive - change to 21 - worsens about 3% \n    \"min_split_gain\": 17.254716573280355,#+ sometimes hypersenstive - change by 0.1 worsens by 0.9%, But for 100 features changes 11.8-12.2 - not change AT ALLL !!!  # Uplifts to  0.43590.. !!! \n    \"reg_lambda\": 0.22423081744595857, #+ seems only makes worse\n    \"reg_alpha\": 0.5069721277162882, #+ seems only makes worse\n    \"subsample_for_bin\": 500, #+ seems  less 10000 - worsens, but after 10 000 does not influence  \n#     \"colsample_bytree\": 1, # only worsens\n#     \"subsample\": 1, # Does not seem to influence at all\n#     \"other_rate\": 1,# no influence ?  \n    \"min_child_weight\":  0.3707947176588889, #+ no influence ?\n    \"subsample_freq\": 9,#+ no influence ? # Integer; alias: bagging_freq; k means perform bagging at every k iteration\n     \"metric\": 'rmse',\n#     \"boosting\": 'dart',\n}\n# {'max_depth': 4, 'min_child_samples': 19, 'n_estimators': 300, 'num_leaves': 13, 'subsample_for_bin': 500, \n# 'subsample_freq': 9, 'min_child_weight': 0.3707947176588889, 'min_split_gain': 17.254716573280355, \n# 'reg_alpha': 0.5069721277162882, 'reg_lambda': 0.22423081744595857}\n\n# {'max_depth': 8, 'min_child_samples': 20, 'n_estimators': 150, 'num_leaves': 9, 'subsample_for_bin': 500, \n# 'subsample_freq': 9, 'min_child_weight': 0.3707947176588889, 'min_split_gain': 17.254716573280355}\n\nHYPER_PARAMETER_TUNE_RANDOM_N_ITER = 100\nLGB_TUNABLE_PARAMETERS_STEP_1 = {\n                                \"max_depth\": list(range(8,10)), \n                                 \"num_leaves\": list(range(8, 11)),\n                                 \"min_child_samples\": list(range(19, 21)),\n                                 \"subsample_for_bin\":[450,500,550],\n                                 \"subsample_freq\": [9,10,11], \n                                \"n_estimators\":[50,150,300],\n#                                 \"random_state\":[0,32]\n                                }\nLGB_TUNABLE_PARAMETERS_STEP_2 = {\n        \"min_split_gain\": stats.loguniform(1, 20),\n        \"min_child_weight\": stats.uniform(0.00, 0.99),\n    }\nXGB_TUNABLE_PARAMETERS_STEP_3 = {\n    \"reg_lambda\": stats.uniform(0.00, 0.99), \n    \"reg_alpha\": stats.uniform(0.00, 0.99)}#LogUniform(1, 1000)\nXGB_TUNABLE_PARAMETERS_STEP_4 = {\n    \"learning_rate\": [1, 0.09, 0.08, 0.07, 0.06]\n}\nLGB_TUNABLE_PARAMETERS = [\n        LGB_TUNABLE_PARAMETERS_STEP_1,\n        LGB_TUNABLE_PARAMETERS_STEP_2,\n        XGB_TUNABLE_PARAMETERS_STEP_3,\n        XGB_TUNABLE_PARAMETERS_STEP_4\n    ]\nLGB_TUNING_OPTIMIZERS = [\"grid\",\"random\",\"random\",\"grid\"] #\"grid\",\"random\"\n\nLGBestimator = (lgbm.LGBMRegressor(random_state = 0, **LGB_PARAMETERS))#gbdt default, goss, rf\n\n# LGBestimator = MultiOutputRegressor(xgb.sklearn.XGBRegressor(**LGB_PARAMETERS))","metadata":{"execution":{"iopub.status.busy":"2022-10-31T23:42:25.000771Z","iopub.execute_input":"2022-10-31T23:42:25.001350Z","iopub.status.idle":"2022-10-31T23:42:25.015212Z","shell.execute_reply.started":"2022-10-31T23:42:25.001315Z","shell.execute_reply":"2022-10-31T23:42:25.013319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selected_target = 'CD31'\nn_features = 40\nNOT_FEATURES = [\n#     \"14\",\"26\",\"30\",\"19\",\"29\",\"35\",\"32\",\"24\",\"25\",\"37\",\"27\",\"15\"\n               ]\n\n#\"22\",\"34\",\"7\",\"40\"\n\n\nFEATURES = [col for col in df_cite.iloc[:1,:n_features].columns if col not in NOT_FEATURES]\n# print(df_cite.iloc[:70988,:n_features].loc[:,FEATURES])\n\nX= df_cite.iloc[:70988,:n_features].loc[:,FEATURES][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,:n_features].loc[:,FEATURES][mask_train]\ny_train = df_cite_train_y[mask_train][ selected_target ]\n#random feature\n# X_train['randNumCol'] = np.random.rand(len(X_train),1)\ndf = pd.concat([X_train, y_train], axis=1)\n# print(X_train.iloc[:,:41])\n# Step 2:\nmask_test_private_like = (df_meta['Playground']==1)&(df_meta['day']==4)\nX_test_private_like = df_cite.iloc[:70988,:n_features].loc[:,FEATURES][ mask_test_private_like  ]\ny_test_private_like = df_cite_train_y[mask_test_private_like][ selected_target ]\nX_test = X_test_private_like\n#random feature\n# X_test['randNumCol'] = np.random.rand(len(X_test),1)\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,:n_features][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","metadata":{"execution":{"iopub.status.busy":"2022-10-31T23:42:28.231335Z","iopub.execute_input":"2022-10-31T23:42:28.231806Z","iopub.status.idle":"2022-10-31T23:42:28.326148Z","shell.execute_reply.started":"2022-10-31T23:42:28.231768Z","shell.execute_reply":"2022-10-31T23:42:28.324973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# splitter = StratificationSplitter(\n#         df, [selected_target], val_size=0.01, n_train_cv_splits=5, max_categories=29,\n#         random_state=42\n#     )\n\n# x_train, y_train = get_data_by_ids(df, splitter.train_ids, selected_target)\n# x_val, y_val = get_data_by_ids(df, splitter.val_ids, selected_target)\n# cv_splits = []\n# for train_ids, val_ids in splitter.cv_ids:\n#     train_idxs = get_indices_by_ids(x_train, train_ids)\n#     val_idxs = get_indices_by_ids(x_train, val_ids)\n#     cv_splits.append((train_idxs, val_idxs))","metadata":{"execution":{"iopub.status.busy":"2022-10-31T23:42:29.395610Z","iopub.execute_input":"2022-10-31T23:42:29.395972Z","iopub.status.idle":"2022-10-31T23:42:29.400395Z","shell.execute_reply.started":"2022-10-31T23:42:29.395941Z","shell.execute_reply":"2022-10-31T23:42:29.399234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"current_best_r2 = 0.4684255866834217\n\nmodel = LGBestimator\nmodel.fit(X_train,y_train)\n# model, best_params = tune_hyperparameters_cv(\n#             LGBestimator, x_train, y_train, LGB_TUNABLE_PARAMETERS, LGB_TUNING_OPTIMIZERS,\n#             HYPER_PARAMETER_TUNE_RANDOM_N_ITER, scoring=make_scorer(mean_squared_error, greater_is_better=False),\n#     #mean_squared_log_error, r2_score mean_squared_error mean_gamma_deviance\n#             n_jobs=1, cv=10, verbose_output=0,\n#             # model key-word params:\n#             eval_set=[(X_test, y_test)],\n#             random_state=42\n#         )\n    \ny_pred = model.predict(X_test)\n\nprint(\"CD31 LGB tuned. R2: {},\\n MSE: {},\\n current_best_r2 improvement: {}\".format(r2_score(y_test, y_pred), \n                mean_squared_error(y_test, y_pred), \n                (r2_score(y_test, y_pred) - current_best_r2 ) / current_best_r2 * 100))\n","metadata":{"execution":{"iopub.status.busy":"2022-10-31T23:42:30.152068Z","iopub.execute_input":"2022-10-31T23:42:30.152683Z","iopub.status.idle":"2022-10-31T23:42:30.339705Z","shell.execute_reply.started":"2022-10-31T23:42:30.152648Z","shell.execute_reply":"2022-10-31T23:42:30.338788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"current_best_r2 = 0.4684255866834217\n\nmodel = LGBestimator\n# model.fit(X_train,y_train)\nmodel, best_params = tune_hyperparameters_cv(\n            LGBestimator, X_train, y_train, LGB_TUNABLE_PARAMETERS, LGB_TUNING_OPTIMIZERS,\n            HYPER_PARAMETER_TUNE_RANDOM_N_ITER, scoring=make_scorer(mean_squared_error, greater_is_better=False),\n    #mean_squared_log_error, r2_score mean_squared_error mean_gamma_deviance\n            n_jobs=1, cv=5, verbose_output=0,\n            # model key-word params:\n            eval_set=[(X_test, y_test)],\n            random_state=42\n        )\n    \ny_pred = model.predict(X_test)\n\nprint(\"CD31 LGB tuned. R2: {},\\n MSE: {},\\n current_best_r2 improvement: {}\".format(r2_score(y_test, y_pred), \n                mean_squared_error(y_test, y_pred), \n                (r2_score(y_test, y_pred) - current_best_r2 ) / current_best_r2 * 100))\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-31T23:03:37.617567Z","iopub.execute_input":"2022-10-31T23:03:37.617921Z","iopub.status.idle":"2022-10-31T23:08:02.591324Z","shell.execute_reply.started":"2022-10-31T23:03:37.617891Z","shell.execute_reply":"2022-10-31T23:08:02.590543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# !pip install xgboost\n# import xgboost as xgb\n# from scipy import stats\n# from sklearn.multioutput import MultiOutputClassifier, MultiOutputRegressor\n\n\n# XGB_PARAMETERS = {\n#             \"booster\": \"gbtree\", \"verbosity\": 1,\n# #             \"objective\": \"reg:squarederror\",# \"eval_metric\": [\"mape\"],\n# #             \"objective\": \"multi:softprob\", \"eval_metric\": \"mlogloss\",#\"merror\",#\n#             \"importance_type\": \"gain\",\n#             \"learning_rate\": 0.05, \"n_estimators\": 5,\n#             \"max_depth\": 5, \"min_child_weight\": 6, \"gamma\": 0,\n#             \"subsample\": 0.8, \"colsample_bytree\": 0.2,\n#             \"reg_lambda\": 10, \"reg_alpha\": 10\n#         }\n  \n# HYPER_PARAMETER_TUNE_RANDOM_N_ITER = 10\n# XGB_TUNABLE_PARAMETERS_STEP_1 = {\"max_depth\": [4], \"min_child_weight\": [6]}#list(range(1, 8))\n# XGB_TUNABLE_PARAMETERS_STEP_2 = {\n#         \"gamma\": stats.uniform(0, 10), \"colsample_bytree\": stats.uniform(0.01, 0.99), \"subsample\": stats.uniform(0.7, 0.3)\n#     }\n# # XGB_TUNABLE_PARAMETERS_STEP_3 = {\"reg_lambda\": LogUniform(1, 1000), \"reg_alpha\": LogUniform(1, 1000)}\n# XGB_TUNABLE_PARAMETERS_STEP_4 = {\"learning_rate\": [0.4, 0.2, 0.1, 0.06, 0.03, 0.01, 0.005]}\n# XGB_TUNABLE_PARAMETERS = [\n#         XGB_TUNABLE_PARAMETERS_STEP_1,\n#         XGB_TUNABLE_PARAMETERS_STEP_2,\n# #         XGB_TUNABLE_PARAMETERS_STEP_3,\n# #         XGB_TUNABLE_PARAMETERS_STEP_4\n#     ]\n# XGB_TUNING_OPTIMIZERS = [\"grid\",\"random\"]\n\n# XGBestimator = xgb.sklearn.XGBRegressor(**XGB_PARAMETERS)\n# XGBestimator = xgb.sklearn.XGBRegressor(**XGB_PARAMETERS)\n\n# model, best_params = tune_hyperparameters_cv(\n    \n#     estimator, x_train, y_train, XGB_TUNABLE_PARAMETERS, XGB_TUNING_OPTIMIZERS,\n#     HYPER_PARAMETER_TUNE_RANDOM_N_ITER, scoring=make_scorer(mean_squared_error, greater_is_better=False),\n#     n_jobs=1, cv=cv_splits, #verbose_output=args.verbose,\n#     # model key-word params:\n#     eval_set=[(x_val, y_val)], early_stopping_rounds=50, verbose=False,\n#     # random_state=42\n#         )\n","metadata":{"execution":{"iopub.status.busy":"2022-10-31T22:00:28.823387Z","iopub.execute_input":"2022-10-31T22:00:28.824340Z","iopub.status.idle":"2022-10-31T22:00:28.830841Z","shell.execute_reply.started":"2022-10-31T22:00:28.824291Z","shell.execute_reply":"2022-10-31T22:00:28.830083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.metrics import r2_score\n\nN_features = 1 # 200 # For the fast test - just take some small number of features \n# When number of features is very low like 1 - the results are more or less the same as just MEAN over the targets - kind of baseline solution (no modeling in fact)\n\ncc = 0\nfor list_folds_indices in  [list_folds_indices_by_days_and_donors, list_folds_indices_by_days_and_donors_with2holdouts]:\n    if cc == 0: print('CV withOUT holdout\\n' ); \n    else : print('\\n\\nCV with holdout\\n' )\n    cc += 1\n    \n    \n    alpha4Ridge  = 1\n\n    X = df_cite.iloc[:,:N_features].values\n    print('X.shape:', X.shape,'Y.shape:', Y.shape)\n\n    from sklearn.linear_model import Ridge\n#     model = Ridge(alpha=alpha4Ridge)\n    #, (train_index, test_index )\n    model = MultiOutputRegressor(LGBestimator)\n\n    df_fold_score_stat = pd.DataFrame();df_fold_score_stat.index.name = 'Fold'\n    t0 = time.time()\n    for fold, indices_tuple  in enumerate( list_folds_indices ):\n        train_index = indices_tuple[0]\n        main_test_index = indices_tuple[1]\n        print('Fold:', fold, 'Shapes of train:', train_index.shape, 'Tests: ',[t.shape for t in indices_tuple[1:] ] )\n        t1 = time.time()\n\n        # Train model: \n        model.fit(X[train_index], Y[train_index]), #eval_set=[(X[train_index], Y[train_index])],\n#         model = tune_hyperparameters_cv(\n#             LGBestimator, X_train,y_train, LGB_TUNABLE_PARAMETERS, LGB_TUNING_OPTIMIZERS,\n#             HYPER_PARAMETER_TUNE_RANDOM_N_ITER, scoring=make_scorer(mean_squared_error, greater_is_better=False),\n#             n_jobs=1, cv=5, #verbose_output=args.verbose,\n#             # model key-word params:\n#             eval_set=[(X_train,y_train)], verbose=False,\n#             random_state=42\n#         )\n        \n\n        # Calculate metrics on test and train folds \n        list_scores = []; list_scores_r2 = []\n        for i_loc in range(0,len( indices_tuple)  ):\n            indices_loc = indices_tuple[i_loc]\n            y_pred = model.predict( X[indices_loc ])\n            #s = model.score( X[indices_loc ], Y[indices_loc]  ) \n            s = correlation_score( y_pred , Y[indices_loc]  )\n            list_scores.append(s)\n            s = r2_score( Y[indices_loc] , y_pred  )\n            list_scores_r2.append(s)\n\n        # Just save statistic for output\n        for i_loc in range(1, len(list_scores  )):\n             df_fold_score_stat.loc[fold, 'Score Test'+str(i_loc)] = list_scores[i_loc]\n        df_fold_score_stat.loc[fold, 'Score Train'] = list_scores[0]\n        df_fold_score_stat.loc[fold, 'Time'] = np.round( (time.time() - t1), 2) \n        for i_loc in range(1, len(list_scores  )):\n             df_fold_score_stat.loc[fold, 'R2 Score Test'+str(i_loc)] = list_scores_r2[i_loc]\n        df_fold_score_stat.loc[fold, 'R2 Score Train'] = list_scores_r2[0]\n\n        print('Correlation scores:', np.round(list_scores,4),  'time:', '%.2f'%(time.time() - t1))\n        \n    display(df_fold_score_stat)\n    display(df_fold_score_stat.describe(percentiles=[] ).iloc[1:,:])\n    display(df_fold_score_stat.mean())","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.590897Z","iopub.status.idle":"2022-10-31T19:48:47.591333Z","shell.execute_reply.started":"2022-10-31T19:48:47.591113Z","shell.execute_reply":"2022-10-31T19:48:47.591134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Features by TruncatedSVD, modeling - Ridge. \n\n\nConclusion - model peforms better on \"priviate like\" test (unseen day&donor), rather than on just unseen donor - quite strange. \nBut it is only for \"r2\"-metric, not for correlation metric \n","metadata":{}},{"cell_type":"code","source":"# %%time\n\n# # Notebook versions 6 and above use prepared features, so we do not need to make PCA/SVD\n\n# N_features = 100\n# r = df_cite.values[:,:N_features]\n\n\n# if 0: # Notebook versions 6 and above use prepared features, so we do not need to make PCA/SVD\n#     from sklearn.decomposition import TruncatedSVD\n#     from scipy.sparse import csr_matrix\n#     import numpy as np\n\n#     # Remark n_components=75 - slightly better results \n#     # For 75 and 100 not monotone difference - on one test is one better, on another - another \n\n#     reducer = TruncatedSVD(n_components=50, n_iter=7, random_state=42)\n#     r = reducer.fit_transform(df_cite_train_x)\n    \n# r.shape","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.592193Z","iopub.status.idle":"2022-10-31T19:48:47.592623Z","shell.execute_reply.started":"2022-10-31T19:48:47.592392Z","shell.execute_reply":"2022-10-31T19:48:47.592411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Concatenating additional features ","metadata":{}},{"cell_type":"code","source":"# list_f = ['/kaggle/input/feature-shop-for-multimodal-singlecell-competition/_citeseq_meta_Gender_and_CellType1Hot.csv',\n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_KmeansFromTruncatedSVD_rs0.csv' \n#     #'/kaggle/input/feature-shop-for-multimodal-singlecell-competition/_citeseq_meta_day_only.csv',\n#          ]\n\n# for f in list_f:\n#     d = pd.read_csv(f, index_col = 0)\n#     r = np.concatenate( (r,d.values), axis = 1)\n# print(r.shape)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.594748Z","iopub.status.idle":"2022-10-31T19:48:47.595265Z","shell.execute_reply.started":"2022-10-31T19:48:47.595003Z","shell.execute_reply":"2022-10-31T19:48:47.595026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %%time\n# from sklearn.metrics import r2_score\n# from sklearn.neural_network import MLPRegressor\n\n# cc = 0\n# for list_folds_indices in  [list_folds_indices_by_days_and_donors, list_folds_indices_by_days_and_donors_with2holdouts]:\n#     if cc == 0: print('CV withOUT holdout\\n' ); \n#     else : print('\\n\\nCV with holdout\\n' )\n#     cc += 1\n    \n\n#     X = r #df_cite_train_x.iloc[:,:N_features].values\n#     print('X.shape:', X.shape,'Y.shape:', Y.shape)\n\n#     model = MLPRegressor(max_iter=500, activation='logistic', early_stopping=True,\n#                          solver='adam', alpha=1e-5, random_state=42, \n#                          hidden_layer_sizes=(300, 200))\n#     #, (train_index, test_index )\n\n#     df_fold_score_stat = pd.DataFrame();df_fold_score_stat.index.name = 'Fold'\n#     t0 = time.time()\n#     for fold, indices_tuple  in enumerate( list_folds_indices ):\n#         train_index = indices_tuple[0]\n#         main_test_index = indices_tuple[1]\n#         print('Fold:', fold, 'Shapes of train:', train_index.shape, 'Tests: ',[t.shape for t in indices_tuple[1:] ] )\n#         t1 = time.time()\n\n#         # Train model: \n#         model.fit(X[train_index], Y[train_index])\n\n#         # Calculate metrics on test and train folds \n#         list_scores = []; list_scores_r2 = []\n#         for i_loc in range(0,len( indices_tuple)  ):\n#             indices_loc = indices_tuple[i_loc]\n#             y_pred = model.predict( X[indices_loc ])\n#             #s = model.score( X[indices_loc ], Y[indices_loc]  ) \n#             s = correlation_score( y_pred , Y[indices_loc]  )\n#             list_scores.append(s)\n#             s = r2_score( Y[indices_loc] , y_pred  )\n#             list_scores_r2.append(s)\n\n#         # Just save statistic for output\n#         for i_loc in range(1, len(list_scores  )):\n#              df_fold_score_stat.loc[fold, 'Score Test'+str(i_loc)] = list_scores[i_loc]\n#         df_fold_score_stat.loc[fold, 'Score Train'] = list_scores[0]\n#         df_fold_score_stat.loc[fold, 'Time'] = np.round( (time.time() - t1), 2) \n#         for i_loc in range(1, len(list_scores  )):\n#              df_fold_score_stat.loc[fold, 'R2 Score Test'+str(i_loc)] = list_scores_r2[i_loc]\n#         df_fold_score_stat.loc[fold, 'R2 Score Train'] = list_scores_r2[0]\n\n#         print('Correlation scores:', np.round(list_scores,4),  'time:', '%.2f'%(time.time() - t1))\n        \n#     display(df_fold_score_stat)\n#     display(df_fold_score_stat.describe(percentiles=[] ).iloc[1:,:])\n#     #display(df_fold_score_stat.mean())","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.596239Z","iopub.status.idle":"2022-10-31T19:48:47.596695Z","shell.execute_reply.started":"2022-10-31T19:48:47.596449Z","shell.execute_reply":"2022-10-31T19:48:47.596470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# cell type and gender added - some improve\n# 0.884071\t0.885336\t0.895725\t7.721667\t0.138083\t0.134252\t0.183134\n# 0.884075\t0.885241\t0.884042\t0.886190\t0.895725\t7.900000\t0.138056\t0.134239\t0.137856\t0.133646\t0.183134\n\n# cell type and gender AND Kmeans - some improve\n# 0.884117\t0.885439\t0.896065\t7.831667\t0.137676\t0.133914\t0.184766\n# 0.884125\t0.885347\t0.884047\t0.886261\t0.896065\t7.990000\t0.137668\t0.133907\t0.137274\t0.133242\t0.184766\n\n# Alpha Ridge = 1e4 ; N_features = 100 - better than oringal \n# 0.883975\t0.885260\t0.895664\t7.445000\t0.137777\t0.133981\t0.182901\n# 0.883978\t0.885165\t0.883948\t0.886112\t0.895664\t7.67000\t0.137748\t0.133969\t0.137561\t0.133367\t0.182901\n\n# Alpha Ridge = 1e4 ; N_features = 150 - worse than 100  \n# 0.883920\t0.885207\t0.895964\t7.670000\t0.136816\t0.133002\t0.184593\n# 0.883924\t0.885112\t0.883889\t0.886064\t0.895964\t7.810000\t0.136792\t0.132988\t0.136555\t0.132391\t0.184593\n\n# Alpha Ridge = 1e4 ; CV withOUT holdout; CV with holdout  \n# 0.883832\t0.885090\t0.896214\t8.015000\t0.135730\t0.131835\t0.186134\n# 0.883833\t0.884991\t0.883825\t0.885982\t0.896214\t8.13000\t0.135705\t0.131815\t0.135470\t0.131264\t0.186134\n\n# Alpha Ridge = 1e4 ; N_features = 50 - worse than original 200 \n# 0.882561\t0.883669\t0.893772\t7.408333\t0.135783\t0.131606\t0.177000\n# 0.882562\t0.883575\t0.882557\t0.884518\t0.893772\t7.691667\t0.135774\t0.131611\t0.135402\t0.130865\t0.177000\n\n\n# Alpha Ridge = 1 ; CV withOUT holdout; CV with holdout  \n# 0.883707\t0.884942\t0.896238\t8.005000\t0.134507\t0.130284\t0.186191\n# 0.883708\t0.884844\t0.883702\t0.885828\t0.896238\t8.141667\t0.134480\t0.130263\t0.134265\t0.129713\t0.18619\n\n# Alpha Ridge = 1e3 ; CV withOUT holdout; CV with holdout  \n# 0.883723\t0.884961\t0.896237\t8.036667\t0.134646\t0.130464\t0.186190\n# 0.883723\t0.884862\t0.883717\t0.885848\t0.896237\t8.216667\t0.134619\t0.130442\t0.134403\t0.129892\t0.186190\n\n# Alpha Ridge = 1e4 ; CV withOUT holdout; CV with holdout  \n# 0.883832\t0.885090\t0.896214\t8.015000\t0.135730\t0.131835\t0.186134\n# 0.883833\t0.884991\t0.883825\t0.885982\t0.896214\t8.13000\t0.135705\t0.131815\t0.135470\t0.131264\t0.186134\n\n# Alpha Ridge = 1e5 ; CV withOUT holdout; CV with holdout  \n# 0.883647\t0.884958\t0.895453\t8.150000\t0.138976\t0.135943\t0.183658\n# 0.883650\t0.884856\t0.883622\t0.885872\t0.895453\t8.113333\t0.138972\t0.135939\t0.138562\t0.135293\t0.18365","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.599339Z","iopub.status.idle":"2022-10-31T19:48:47.600432Z","shell.execute_reply.started":"2022-10-31T19:48:47.600097Z","shell.execute_reply":"2022-10-31T19:48:47.600137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calculation predictions for  submission","metadata":{}},{"cell_type":"markdown","source":"## Way 1 - train selected model on full train","metadata":{}},{"cell_type":"code","source":"# paths = ['../output/kaggle/working/submission.csv']\n# dfs = [pd.read_csv(x) for x in paths]\n# print(dfs)","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.601591Z","iopub.status.idle":"2022-10-31T19:48:47.601956Z","shell.execute_reply.started":"2022-10-31T19:48:47.601775Z","shell.execute_reply":"2022-10-31T19:48:47.601792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.602921Z","iopub.status.idle":"2022-10-31T19:48:47.603244Z","shell.execute_reply.started":"2022-10-31T19:48:47.603088Z","shell.execute_reply":"2022-10-31T19:48:47.603103Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nmodel.fit(X[:70988,:], Y)","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.604472Z","iopub.status.idle":"2022-10-31T19:48:47.604823Z","shell.execute_reply.started":"2022-10-31T19:48:47.604669Z","shell.execute_reply":"2022-10-31T19:48:47.604685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nY_pred4submit = model.predict(X[70988:,:])\nprint(Y_pred4submit.shape)\nprint(Y_pred4submit[:3,:3])","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.606230Z","iopub.status.idle":"2022-10-31T19:48:47.606589Z","shell.execute_reply.started":"2022-10-31T19:48:47.606396Z","shell.execute_reply":"2022-10-31T19:48:47.606412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparation of the submission file","metadata":{}},{"cell_type":"code","source":"%%time\n\nmode_subm = 'USE: sskknts MSCI CITEseq Keras Quickstart + Dropout' # 'put_zeros_to_multiome_part'\n\nif mode_subm == 'USE: sskknts MSCI CITEseq Keras Quickstart + Dropout' :\n    df_submission_full = pd.read_csv('/kaggle/input/msci-citeseq-keras-quickstart-dropout/submission.csv',\n                             index_col='row_id', squeeze=True)\n    df_submission_full = df_submission_full.to_frame()\nelif mode_subm == 'put_zeros_to_multiome_part':\n    df_submission_full = pd.DataFrame(index = range(65_744_180), columns = ['target'], data = np.zeros(65_744_180) )\n    df_submission_full.index.name = 'row_id'\n    \ndisplay(df_submission_full.info() )\nprint()\ndisplay(df_submission_full)","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.607451Z","iopub.status.idle":"2022-10-31T19:48:47.607857Z","shell.execute_reply.started":"2022-10-31T19:48:47.607703Z","shell.execute_reply":"2022-10-31T19:48:47.607719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# We can prepare CITE-seq part as: Y_pred4submit.values.ravel() -- see notebook: https://www.kaggle.com/code/alexandervc/mmscel-analysis-of-evaluation-and-submission\ndf_submission_full['target'].iloc[:6_812_820] = Y_pred4submit.ravel()\ndisplay(df_submission_full)","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.609132Z","iopub.status.idle":"2022-10-31T19:48:47.609428Z","shell.execute_reply.started":"2022-10-31T19:48:47.609284Z","shell.execute_reply":"2022-10-31T19:48:47.609299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save submission file - about 1.46 minutes ","metadata":{}},{"cell_type":"code","source":"%%time\n# Wall time: 1min 46s\nsubmit_filename_postfix = 'CVschemesV14tryPCA_addMetaAndKmeans_RescaleY_RidgeAlpha1e4Ndim100_multiome_sskknts_KerDrop'\n\ndf_submission_full.to_csv('submission.csv')\n# df_submission_full.to_csv('submission_cite_seq_'+submit_filename_postfix +'.csv')","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.610468Z","iopub.status.idle":"2022-10-31T19:48:47.610830Z","shell.execute_reply.started":"2022-10-31T19:48:47.610681Z","shell.execute_reply":"2022-10-31T19:48:47.610697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import os\n# os.listdir()\n# #os.remove( 'submission_cite_seq_CVschemesV09RescaleY_Ridge_multiome_sskknts_KerDrop.csv')","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.612434Z","iopub.status.idle":"2022-10-31T19:48:47.612755Z","shell.execute_reply.started":"2022-10-31T19:48:47.612608Z","shell.execute_reply":"2022-10-31T19:48:47.612622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\n","metadata":{"execution":{"iopub.status.busy":"2022-10-31T19:48:47.613824Z","iopub.status.idle":"2022-10-31T19:48:47.614121Z","shell.execute_reply.started":"2022-10-31T19:48:47.613975Z","shell.execute_reply":"2022-10-31T19:48:47.613990Z"},"trusted":true},"execution_count":null,"outputs":[]}]}