{"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 = True","metadata":{"execution":{"iopub.status.busy":"2022-10-26T08:42:49.591118Z","iopub.execute_input":"2022-10-26T08:42:49.591511Z","iopub.status.idle":"2022-10-26T08:42:49.597719Z","shell.execute_reply.started":"2022-10-26T08:42:49.591471Z","shell.execute_reply":"2022-10-26T08:42:49.596283Z"},"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-26T08:42:49.599518Z","iopub.execute_input":"2022-10-26T08:42:49.600024Z","iopub.status.idle":"2022-10-26T08:42:49.656808Z","shell.execute_reply.started":"2022-10-26T08:42:49.599958Z","shell.execute_reply":"2022-10-26T08:42:49.655807Z"},"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-26T08:42:49.658909Z","iopub.execute_input":"2022-10-26T08:42:49.659606Z","iopub.status.idle":"2022-10-26T08:43:12.259915Z","shell.execute_reply.started":"2022-10-26T08:42:49.659568Z","shell.execute_reply":"2022-10-26T08:43:12.258501Z"},"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-26T08:43:17.847877Z","iopub.execute_input":"2022-10-26T08:43:17.848329Z","iopub.status.idle":"2022-10-26T08:43:17.856416Z","shell.execute_reply.started":"2022-10-26T08:43:17.848288Z","shell.execute_reply":"2022-10-26T08:43:17.855046Z"},"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-26T08:43:18.609090Z","iopub.execute_input":"2022-10-26T08:43:18.609915Z","iopub.status.idle":"2022-10-26T08:43:38.763556Z","shell.execute_reply.started":"2022-10-26T08:43:18.609871Z","shell.execute_reply":"2022-10-26T08:43:38.762381Z"},"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-26T08:43:38.765681Z","iopub.execute_input":"2022-10-26T08:43:38.766168Z","iopub.status.idle":"2022-10-26T08:43:39.649224Z","shell.execute_reply.started":"2022-10-26T08:43:38.766131Z","shell.execute_reply":"2022-10-26T08:43:39.647703Z"},"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-26T08:43:39.650872Z","iopub.execute_input":"2022-10-26T08:43:39.651284Z","iopub.status.idle":"2022-10-26T08:43:40.072199Z","shell.execute_reply.started":"2022-10-26T08:43:39.651239Z","shell.execute_reply":"2022-10-26T08:43:40.070923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta['donor'].unique()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T08:43:40.075029Z","iopub.execute_input":"2022-10-26T08:43:40.075930Z","iopub.status.idle":"2022-10-26T08:43:40.087662Z","shell.execute_reply.started":"2022-10-26T08:43:40.075879Z","shell.execute_reply":"2022-10-26T08:43:40.086266Z"},"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-26T08:43:40.089081Z","iopub.execute_input":"2022-10-26T08:43:40.089473Z","iopub.status.idle":"2022-10-26T08:43:40.105324Z","shell.execute_reply.started":"2022-10-26T08:43:40.089439Z","shell.execute_reply":"2022-10-26T08:43:40.104153Z"},"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-26T08:43:40.106719Z","iopub.execute_input":"2022-10-26T08:43:40.107077Z","iopub.status.idle":"2022-10-26T08:43:40.129146Z","shell.execute_reply.started":"2022-10-26T08:43:40.107045Z","shell.execute_reply":"2022-10-26T08:43:40.127879Z"},"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-26T08:49:11.367251Z","iopub.execute_input":"2022-10-26T08:49:11.368422Z","iopub.status.idle":"2022-10-26T08:49:11.383436Z","shell.execute_reply.started":"2022-10-26T08:49:11.368367Z","shell.execute_reply":"2022-10-26T08:49:11.381847Z"},"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-26T08:49:11.905811Z","iopub.execute_input":"2022-10-26T08:49:11.906577Z","iopub.status.idle":"2022-10-26T08:49:11.916113Z","shell.execute_reply.started":"2022-10-26T08:49:11.906493Z","shell.execute_reply":"2022-10-26T08:49:11.914864Z"},"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-26T08:49:14.183407Z","iopub.execute_input":"2022-10-26T08:49:14.183839Z","iopub.status.idle":"2022-10-26T08:49:14.204808Z","shell.execute_reply.started":"2022-10-26T08:49:14.183802Z","shell.execute_reply":"2022-10-26T08:49:14.203575Z"},"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-26T08:49:16.315232Z","iopub.execute_input":"2022-10-26T08:49:16.318150Z","iopub.status.idle":"2022-10-26T08:49:16.325185Z","shell.execute_reply.started":"2022-10-26T08:49:16.318104Z","shell.execute_reply":"2022-10-26T08:49:16.323922Z"},"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-26T08:49:18.057714Z","iopub.execute_input":"2022-10-26T08:49:18.058175Z","iopub.status.idle":"2022-10-26T08:49:18.442435Z","shell.execute_reply.started":"2022-10-26T08:49:18.058139Z","shell.execute_reply":"2022-10-26T08:49:18.441220Z"},"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-26T08:49:26.620101Z","iopub.execute_input":"2022-10-26T08:49:26.620501Z","iopub.status.idle":"2022-10-26T08:49:26.654334Z","shell.execute_reply.started":"2022-10-26T08:49:26.620470Z","shell.execute_reply":"2022-10-26T08:49:26.653253Z"},"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')","metadata":{"execution":{"iopub.status.busy":"2022-10-26T08:49:28.440640Z","iopub.execute_input":"2022-10-26T08:49:28.441423Z","iopub.status.idle":"2022-10-26T08:49:28.482003Z","shell.execute_reply.started":"2022-10-26T08:49:28.441373Z","shell.execute_reply":"2022-10-26T08:49:28.480755Z"},"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\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\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-26T08:49:41.056624Z","iopub.execute_input":"2022-10-26T08:49:41.057464Z","iopub.status.idle":"2022-10-26T08:51:11.001577Z","shell.execute_reply.started":"2022-10-26T08:49:41.057417Z","shell.execute_reply":"2022-10-26T08:51:11.000722Z"},"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\nN_features = 100\nr = df_cite.values[:,:N_features]\n\n\nif 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    \nr.shape","metadata":{"execution":{"iopub.status.busy":"2022-10-26T08:51:11.003229Z","iopub.execute_input":"2022-10-26T08:51:11.004129Z","iopub.status.idle":"2022-10-26T08:51:11.014028Z","shell.execute_reply.started":"2022-10-26T08:51:11.004094Z","shell.execute_reply":"2022-10-26T08:51:11.012827Z"},"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\nfor f in list_f:\n    d = pd.read_csv(f, index_col = 0)\n    r = np.concatenate( (r,d.values), axis = 1)\nprint(r.shape)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-26T08:51:52.864311Z","iopub.execute_input":"2022-10-26T08:51:52.864734Z","iopub.status.idle":"2022-10-26T08:51:53.885734Z","shell.execute_reply.started":"2022-10-26T08:51:52.864703Z","shell.execute_reply":"2022-10-26T08:51:53.884498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.metrics import r2_score\nfrom sklearn.neural_network import MLPRegressor\nfrom sklearn.linear_model import Ridge\nfrom sklearn.model_selection import GridSearchCV\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    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    \n    # new code\n    gs = GridSearchCV(\n            estimator=Ridge(),\n            param_grid = {'alpha': np.logspace(-3, 0, 15)})\n    gs.fit(X[train_index], Y[train_index])\n    model = gs.best_estimator_\n    \n    \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-26T09:06:38.463409Z","iopub.execute_input":"2022-10-26T09:06:38.463828Z","iopub.status.idle":"2022-10-26T09:08:42.756505Z","shell.execute_reply.started":"2022-10-26T09:06:38.463792Z","shell.execute_reply":"2022-10-26T09:08:42.755157Z"},"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-17T07:29:12.164318Z","iopub.execute_input":"2022-10-17T07:29:12.164942Z","iopub.status.idle":"2022-10-17T07:29:12.170818Z","shell.execute_reply.started":"2022-10-17T07:29:12.164908Z","shell.execute_reply":"2022-10-17T07:29:12.169758Z"},"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":"model","metadata":{"execution":{"iopub.status.busy":"2022-10-26T08:53:40.596582Z","iopub.execute_input":"2022-10-26T08:53:40.597548Z","iopub.status.idle":"2022-10-26T08:53:40.606340Z","shell.execute_reply.started":"2022-10-26T08:53:40.597506Z","shell.execute_reply":"2022-10-26T08:53:40.605120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nmodel.fit(X[:70988,:], Y)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T08:54:45.113459Z","iopub.execute_input":"2022-10-26T08:54:45.113945Z","iopub.status.idle":"2022-10-26T08:54:45.479463Z","shell.execute_reply.started":"2022-10-26T08:54:45.113909Z","shell.execute_reply":"2022-10-26T08:54:45.477442Z"},"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-26T08:54:46.908803Z","iopub.execute_input":"2022-10-26T08:54:46.909269Z","iopub.status.idle":"2022-10-26T08:54:47.047522Z","shell.execute_reply.started":"2022-10-26T08:54:46.909230Z","shell.execute_reply":"2022-10-26T08:54:47.045757Z"},"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-26T08:54:50.861763Z","iopub.execute_input":"2022-10-26T08:54:50.862190Z","iopub.status.idle":"2022-10-26T08:56:09.651125Z","shell.execute_reply.started":"2022-10-26T08:54:50.862151Z","shell.execute_reply":"2022-10-26T08:56:09.649611Z"},"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-26T08:56:09.653311Z","iopub.execute_input":"2022-10-26T08:56:09.653676Z","iopub.status.idle":"2022-10-26T08:56:09.678634Z","shell.execute_reply.started":"2022-10-26T08:56:09.653643Z","shell.execute_reply":"2022-10-26T08:56:09.677334Z"},"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'\ndf_submission_full.to_csv('submission_cite_seq_'+submit_filename_postfix +'.csv')\ndf_submission_full.to_csv('submission.csv')\n\n","metadata":{"execution":{"iopub.status.busy":"2022-10-26T08:57:31.807694Z","iopub.execute_input":"2022-10-26T08:57:31.808147Z","iopub.status.idle":"2022-10-26T08:59:49.134131Z","shell.execute_reply.started":"2022-10-26T08:57:31.808111Z","shell.execute_reply":"2022-10-26T08:59:49.132709Z"},"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-26T08:59:49.136171Z","iopub.execute_input":"2022-10-26T08:59:49.136544Z","iopub.status.idle":"2022-10-26T08:59:49.141809Z","shell.execute_reply.started":"2022-10-26T08:59:49.136509Z","shell.execute_reply":"2022-10-26T08:59:49.140608Z"},"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-17T07:32:55.706605Z","iopub.execute_input":"2022-10-17T07:32:55.706999Z","iopub.status.idle":"2022-10-17T07:32:55.720057Z","shell.execute_reply.started":"2022-10-17T07:32:55.706968Z","shell.execute_reply":"2022-10-17T07:32:55.718687Z"},"trusted":true},"execution_count":null,"outputs":[]}]}