{"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":"## RESOURCES:\n\nhttps://towardsdatascience.com/5-reasons-why-you-should-use-cross-validation-in-your-data-science-project-8163311a1e79\n\nhttps://www.kaggle.com/code/alexandervc/mmscel-crossvalidation-schemes\n\n\nTWO test sets - 1) similar to public LB 2) similar to private LB","metadata":{}},{"cell_type":"code","source":"## RESOURCES:\n\n# https://towardsdatascience.com/5-reasons-why-you-should-use-cross-validation-in-your-data-science-project-8163311a1e79\n\n# https://www.kaggle.com/code/alexandervc/mmscel-crossvalidation-schemes\n\n\n# TWO test sets - 1) similar to public LB 2) similar to private LB","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rescale_Y_to_mean0_std1 = True","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:26:52.82114Z","iopub.execute_input":"2022-11-02T13:26:52.821824Z","iopub.status.idle":"2022-11-02T13:26:52.832018Z","shell.execute_reply.started":"2022-11-02T13:26:52.821754Z","shell.execute_reply":"2022-11-02T13:26:52.830832Z"},"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-11-02T13:26:52.839306Z","iopub.execute_input":"2022-11-02T13:26:52.839782Z","iopub.status.idle":"2022-11-02T13:26:52.853482Z","shell.execute_reply.started":"2022-11-02T13:26:52.839731Z","shell.execute_reply":"2022-11-02T13:26:52.852176Z"},"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-11-02T13:26:52.855325Z","iopub.execute_input":"2022-11-02T13:26:52.855796Z","iopub.status.idle":"2022-11-02T13:27:16.441344Z","shell.execute_reply.started":"2022-11-02T13:26:52.85575Z","shell.execute_reply":"2022-11-02T13:27:16.439951Z"},"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-11-02T13:27:16.446931Z","iopub.execute_input":"2022-11-02T13:27:16.447692Z","iopub.status.idle":"2022-11-02T13:27:16.457575Z","shell.execute_reply.started":"2022-11-02T13:27:16.447293Z","shell.execute_reply":"2022-11-02T13:27:16.456314Z"},"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-11-02T13:27:16.459145Z","iopub.execute_input":"2022-11-02T13:27:16.45961Z","iopub.status.idle":"2022-11-02T13:27:30.695362Z","shell.execute_reply.started":"2022-11-02T13:27:16.45956Z","shell.execute_reply":"2022-11-02T13:27:30.694104Z"},"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-11-02T13:27:30.697249Z","iopub.execute_input":"2022-11-02T13:27:30.698068Z","iopub.status.idle":"2022-11-02T13:27:31.47479Z","shell.execute_reply.started":"2022-11-02T13:27:30.698013Z","shell.execute_reply":"2022-11-02T13:27:31.47369Z"},"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-11-02T13:27:31.476267Z","iopub.execute_input":"2022-11-02T13:27:31.476635Z","iopub.status.idle":"2022-11-02T13:27:31.788821Z","shell.execute_reply.started":"2022-11-02T13:27:31.4766Z","shell.execute_reply":"2022-11-02T13:27:31.787739Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta['donor'].unique()","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:27:31.790045Z","iopub.execute_input":"2022-11-02T13:27:31.790377Z","iopub.status.idle":"2022-11-02T13:27:31.799039Z","shell.execute_reply.started":"2022-11-02T13:27:31.790346Z","shell.execute_reply":"2022-11-02T13:27:31.797623Z"},"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-11-02T13:27:31.800356Z","iopub.execute_input":"2022-11-02T13:27:31.800696Z","iopub.status.idle":"2022-11-02T13:27:31.810726Z","shell.execute_reply.started":"2022-11-02T13:27:31.800664Z","shell.execute_reply":"2022-11-02T13:27:31.809519Z"},"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-11-02T13:27:31.812246Z","iopub.execute_input":"2022-11-02T13:27:31.812689Z","iopub.status.idle":"2022-11-02T13:27:31.827346Z","shell.execute_reply.started":"2022-11-02T13:27:31.812658Z","shell.execute_reply":"2022-11-02T13:27:31.826131Z"},"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-11-02T13:27:31.828329Z","iopub.execute_input":"2022-11-02T13:27:31.828659Z","iopub.status.idle":"2022-11-02T13:27:31.841266Z","shell.execute_reply.started":"2022-11-02T13:27:31.828629Z","shell.execute_reply":"2022-11-02T13:27:31.840159Z"},"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-11-02T13:27:31.842915Z","iopub.execute_input":"2022-11-02T13:27:31.84326Z","iopub.status.idle":"2022-11-02T13:27:31.851547Z","shell.execute_reply.started":"2022-11-02T13:27:31.843229Z","shell.execute_reply":"2022-11-02T13:27:31.850306Z"},"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-11-02T13:27:31.857351Z","iopub.execute_input":"2022-11-02T13:27:31.857794Z","iopub.status.idle":"2022-11-02T13:27:31.879603Z","shell.execute_reply.started":"2022-11-02T13:27:31.85776Z","shell.execute_reply":"2022-11-02T13:27:31.878308Z"},"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-11-02T13:27:31.881511Z","iopub.execute_input":"2022-11-02T13:27:31.881994Z","iopub.status.idle":"2022-11-02T13:27:31.889014Z","shell.execute_reply.started":"2022-11-02T13:27:31.881946Z","shell.execute_reply":"2022-11-02T13:27:31.887924Z"},"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-11-02T13:27:31.890719Z","iopub.execute_input":"2022-11-02T13:27:31.891191Z","iopub.status.idle":"2022-11-02T13:27:32.181649Z","shell.execute_reply.started":"2022-11-02T13:27:31.891145Z","shell.execute_reply":"2022-11-02T13:27:32.180482Z"},"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-11-02T13:27:32.183885Z","iopub.execute_input":"2022-11-02T13:27:32.184352Z","iopub.status.idle":"2022-11-02T13:27:32.224151Z","shell.execute_reply.started":"2022-11-02T13:27:32.184305Z","shell.execute_reply":"2022-11-02T13:27:32.222893Z"},"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-11-02T13:27:32.225679Z","iopub.execute_input":"2022-11-02T13:27:32.227007Z","iopub.status.idle":"2022-11-02T13:27:32.271486Z","shell.execute_reply.started":"2022-11-02T13:27:32.226953Z","shell.execute_reply":"2022-11-02T13:27:32.270216Z"},"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-11-02T13:27:32.273352Z","iopub.execute_input":"2022-11-02T13:27:32.273775Z","iopub.status.idle":"2022-11-02T13:29:01.262822Z","shell.execute_reply.started":"2022-11-02T13:27:32.273734Z","shell.execute_reply":"2022-11-02T13:29:01.261409Z"},"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-11-02T13:29:01.264395Z","iopub.execute_input":"2022-11-02T13:29:01.264745Z","iopub.status.idle":"2022-11-02T13:29:01.275523Z","shell.execute_reply.started":"2022-11-02T13:29:01.264704Z","shell.execute_reply":"2022-11-02T13:29:01.274281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Concatenating additional features ","metadata":{}},{"cell_type":"code","source":"X = r\n","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:29:01.276536Z","iopub.execute_input":"2022-11-02T13:29:01.276889Z","iopub.status.idle":"2022-11-02T13:29:01.895565Z","shell.execute_reply.started":"2022-11-02T13:29:01.276842Z","shell.execute_reply":"2022-11-02T13:29:01.894284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# codes for bold text in print to highlight hyperparameter value\nstart = \"\\033[1m\"; end = \"\\033[0;0m\"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.metrics import r2_score\nfrom sklearn.linear_model import Ridge\n\ndf_fold_score_stat_total = []\n\nfor alpha2test in np.logspace(-5, 2, 8):\n    print(f'\\n{start}Ridge_alpha = {alpha2test}{end}\\n')\n    cc = 0\n    for list_folds_indices in  [list_folds_indices_by_days_and_donors,\n                                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#         X = r\n#         print('X.shape:', X.shape,'Y.shape:', Y.shape)\n\n        model = Ridge(alpha=alpha2test, random_state=42)\n        #, (train_index, test_index )\n\n        df_fold_score_stat = pd.DataFrame()\n        df_fold_score_stat.index.name = 'Fold'\n        t0 = time.time()\n    \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,\n                  'Shapes of train:', train_index.shape,\n                  'Tests: ',[t.shape for t in indices_tuple[1::2] ] )\n            t1 = time.time()\n\n            # Train model: \n            model.fit(X[train_index], Y[train_index])\n\n            # Calculate metrics on test folds \n            list_scores = []\n            list_scores_r2 = []\n            \n            for i_loc in range(1, len(indices_tuple) , 2 ):\n                # score only on test like private LB\n                indices_loc = indices_tuple[i_loc]\n                y_pred = model.predict( X[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(0, len(list_scores)):\n                 df_fold_score_stat.loc[fold, 'Score Test'+str(i_loc+1)] = list_scores[i_loc]\n            df_fold_score_stat.loc[fold, 'Time'] = np.round( (time.time() - t1), 2) \n            for i_loc in range(0, len(list_scores )):\n                 df_fold_score_stat.loc[fold, 'R2 Score Test'+str(i_loc+1)] = list_scores_r2[i_loc]\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        df_fold_score_stat = df_fold_score_stat.describe(percentiles=[] ).iloc[1:,:]\n        df_fold_score_stat['alpha'] = alpha2test\n        df_fold_score_stat_total.append(df_fold_score_stat)\n\ndf_fold_score_stat_total = pd.concat(df_fold_score_stat_total, axis=0)\ndf_fold_score_stat_total.reset_index(inplace=True); df_fold_score_stat_total.set_index(['alpha', 'index'], inplace=True)\n","metadata":{},"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-11-02T13:30:37.309278Z","iopub.execute_input":"2022-11-02T13:30:37.310097Z","iopub.status.idle":"2022-11-02T13:30:37.318186Z","shell.execute_reply.started":"2022-11-02T13:30:37.310042Z","shell.execute_reply":"2022-11-02T13:30:37.316884Z"},"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-11-02T13:30:37.32003Z","iopub.execute_input":"2022-11-02T13:30:37.320515Z","iopub.status.idle":"2022-11-02T13:30:37.337887Z","shell.execute_reply.started":"2022-11-02T13:30:37.320467Z","shell.execute_reply":"2022-11-02T13:30:37.336438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nmodel.fit(X[:70988,:], Y)","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:30:37.341205Z","iopub.execute_input":"2022-11-02T13:30:37.342097Z","iopub.status.idle":"2022-11-02T13:30:37.70231Z","shell.execute_reply.started":"2022-11-02T13:30:37.342043Z","shell.execute_reply":"2022-11-02T13:30:37.701029Z"},"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-11-02T13:30:37.703947Z","iopub.execute_input":"2022-11-02T13:30:37.70474Z","iopub.status.idle":"2022-11-02T13:30:37.849727Z","shell.execute_reply.started":"2022-11-02T13:30:37.704677Z","shell.execute_reply":"2022-11-02T13:30:37.848467Z"},"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-11-02T13:30:37.851558Z","iopub.execute_input":"2022-11-02T13:30:37.852302Z","iopub.status.idle":"2022-11-02T13:31:34.536486Z","shell.execute_reply.started":"2022-11-02T13:30:37.852253Z","shell.execute_reply":"2022-11-02T13:31:34.535307Z"},"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-11-02T13:31:34.538414Z","iopub.execute_input":"2022-11-02T13:31:34.539114Z","iopub.status.idle":"2022-11-02T13:31:34.5647Z","shell.execute_reply.started":"2022-11-02T13:31:34.539074Z","shell.execute_reply":"2022-11-02T13:31:34.563621Z"},"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')","metadata":{"execution":{"iopub.status.busy":"2022-11-02T13:31:34.566814Z","iopub.execute_input":"2022-11-02T13:31:34.567706Z","iopub.status.idle":"2022-11-02T13:33:57.643893Z","shell.execute_reply.started":"2022-11-02T13:31:34.567657Z","shell.execute_reply":"2022-11-02T13:33:57.642631Z"},"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-11-02T13:33:57.64528Z","iopub.execute_input":"2022-11-02T13:33:57.645617Z","iopub.status.idle":"2022-11-02T13:33:57.650241Z","shell.execute_reply.started":"2022-11-02T13:33:57.645584Z","shell.execute_reply":"2022-11-02T13:33:57.649063Z"},"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-11-02T13:33:57.652648Z","iopub.execute_input":"2022-11-02T13:33:57.653026Z","iopub.status.idle":"2022-11-02T13:33:57.665591Z","shell.execute_reply.started":"2022-11-02T13:33:57.652989Z","shell.execute_reply":"2022-11-02T13:33:57.664588Z"},"trusted":true},"execution_count":null,"outputs":[]}]}