{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"rescale_Y_to_mean0_std1 = True\nadditional_feature_analysis_part = True #  False # Optional part at the end of the notebook \n    # which explores effects from the additional feautures\n    # It can be time consuming part takeas about 1 h 15 minutes for RIDGE and about 1500 features ","metadata":{"execution":{"iopub.status.busy":"2022-10-25T13:35:33.780777Z","iopub.execute_input":"2022-10-25T13:35:33.781188Z","iopub.status.idle":"2022-10-25T13:35:33.786477Z","shell.execute_reply.started":"2022-10-25T13:35:33.781138Z","shell.execute_reply":"2022-10-25T13:35:33.785112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Install/import","metadata":{}},{"cell_type":"code","source":"import 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-25T13:35:33.797283Z","iopub.execute_input":"2022-10-25T13:35:33.798061Z","iopub.status.idle":"2022-10-25T13:36:01.002048Z","shell.execute_reply.started":"2022-10-25T13:35:33.798027Z","shell.execute_reply":"2022-10-25T13:36:01.000376Z"},"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-25T13:36:01.005205Z","iopub.execute_input":"2022-10-25T13:36:01.005605Z","iopub.status.idle":"2022-10-25T13:36:01.014809Z","shell.execute_reply.started":"2022-10-25T13:36:01.005569Z","shell.execute_reply":"2022-10-25T13:36:01.013081Z"},"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-25T13:36:01.016611Z","iopub.execute_input":"2022-10-25T13:36:01.017105Z","iopub.status.idle":"2022-10-25T13:36:23.015254Z","shell.execute_reply.started":"2022-10-25T13:36:01.017061Z","shell.execute_reply":"2022-10-25T13:36:23.013745Z"},"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-25T13:36:23.018702Z","iopub.execute_input":"2022-10-25T13:36:23.019263Z","iopub.status.idle":"2022-10-25T13:36:23.866278Z","shell.execute_reply.started":"2022-10-25T13:36:23.019212Z","shell.execute_reply":"2022-10-25T13:36:23.864608Z"},"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-25T13:36:23.867985Z","iopub.execute_input":"2022-10-25T13:36:23.86837Z","iopub.status.idle":"2022-10-25T13:36:24.327331Z","shell.execute_reply.started":"2022-10-25T13:36:23.868336Z","shell.execute_reply":"2022-10-25T13:36:24.325841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta['donor'].unique()","metadata":{"execution":{"iopub.status.busy":"2022-10-25T13:36:24.329071Z","iopub.execute_input":"2022-10-25T13:36:24.329374Z","iopub.status.idle":"2022-10-25T13:36:24.341027Z","shell.execute_reply.started":"2022-10-25T13:36:24.329345Z","shell.execute_reply":"2022-10-25T13:36:24.339582Z"},"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-25T13:36:24.343074Z","iopub.execute_input":"2022-10-25T13:36:24.343379Z","iopub.status.idle":"2022-10-25T13:36:24.355574Z","shell.execute_reply.started":"2022-10-25T13:36:24.34335Z","shell.execute_reply":"2022-10-25T13:36:24.354133Z"},"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-25T13:36:24.357334Z","iopub.execute_input":"2022-10-25T13:36:24.357739Z","iopub.status.idle":"2022-10-25T13:36:24.371977Z","shell.execute_reply.started":"2022-10-25T13:36:24.357696Z","shell.execute_reply":"2022-10-25T13:36:24.37064Z"},"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-25T13:36:24.373129Z","iopub.execute_input":"2022-10-25T13:36:24.374219Z","iopub.status.idle":"2022-10-25T13:36:24.386121Z","shell.execute_reply.started":"2022-10-25T13:36:24.374184Z","shell.execute_reply":"2022-10-25T13:36:24.384826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check:\nprint( train_index.shape, test_index.shape, np.intersect1d(test_index, train_index ) )\n","metadata":{"execution":{"iopub.status.busy":"2022-10-25T13:36:24.392664Z","iopub.execute_input":"2022-10-25T13:36:24.393046Z","iopub.status.idle":"2022-10-25T13:36:24.401919Z","shell.execute_reply.started":"2022-10-25T13:36:24.393015Z","shell.execute_reply":"2022-10-25T13:36:24.400533Z"},"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-25T13:36:24.403706Z","iopub.execute_input":"2022-10-25T13:36:24.404127Z","iopub.status.idle":"2022-10-25T13:36:24.42564Z","shell.execute_reply.started":"2022-10-25T13:36:24.404095Z","shell.execute_reply":"2022-10-25T13:36:24.42456Z"},"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-25T13:36:24.427163Z","iopub.execute_input":"2022-10-25T13:36:24.428114Z","iopub.status.idle":"2022-10-25T13:36:24.435783Z","shell.execute_reply.started":"2022-10-25T13:36:24.428082Z","shell.execute_reply":"2022-10-25T13:36:24.434321Z"},"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-25T13:36:24.437809Z","iopub.execute_input":"2022-10-25T13:36:24.438332Z","iopub.status.idle":"2022-10-25T13:36:24.68339Z","shell.execute_reply.started":"2022-10-25T13:36:24.43829Z","shell.execute_reply":"2022-10-25T13:36:24.682143Z"},"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-25T13:36:24.685268Z","iopub.execute_input":"2022-10-25T13:36:24.685793Z","iopub.status.idle":"2022-10-25T13:36:24.723332Z","shell.execute_reply.started":"2022-10-25T13:36:24.685724Z","shell.execute_reply":"2022-10-25T13:36:24.722037Z"},"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-25T13:36:24.724976Z","iopub.execute_input":"2022-10-25T13:36:24.725884Z","iopub.status.idle":"2022-10-25T13:36:24.763701Z","shell.execute_reply.started":"2022-10-25T13:36:24.725848Z","shell.execute_reply":"2022-10-25T13:36:24.762236Z"},"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\nfrom sklearn.neighbors import KNeighborsRegressor\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    model = KNeighborsRegressor(n_neighbors=600)\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-25T13:52:00.906792Z","iopub.execute_input":"2022-10-25T13:52:00.90798Z","iopub.status.idle":"2022-10-25T13:55:19.592363Z","shell.execute_reply.started":"2022-10-25T13:52:00.907934Z","shell.execute_reply":"2022-10-25T13:55:19.59139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Features by TruncatedSVD, modeling - KNNRegressor. \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-25T13:56:12.736237Z","iopub.execute_input":"2022-10-25T13:56:12.736735Z","iopub.status.idle":"2022-10-25T13:56:12.748587Z","shell.execute_reply.started":"2022-10-25T13:56:12.736701Z","shell.execute_reply":"2022-10-25T13:56:12.747256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Concatenating additional features ","metadata":{}},{"cell_type":"code","source":"list_f = [ '/kaggle/input/top-100-closest-to-target-genes-with-features/split_by_target/anti-human_CD63_semantic_similarity.csv',\n'/kaggle/input/top-100-closest-to-target-genes-with-features/split_by_target/anti-human_CD49b_semantic_similarity.csv',\n'/kaggle/input/top-100-closest-to-target-genes-with-features/split_by_target/anti-human_CD62L_semantic_similarity.csv',\n'/kaggle/input/top-100-closest-to-target-genes-with-features/split_by_target/anti-human_HLA-DR_semantic_similarity.csv'    ,\n]\n\n# [ '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_TRIMAPfromTSVD200_n_components10.csv',\n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_UMAPfromPCA500_n_components50_n_neighbors5_min_dist0d1_rs42.csv',\n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_NCVISfromPCA500.csv',\n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_UMAPfromTSVD200_n_components10_n_neighbors5_min_dist0d9_rs42.csv', \n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_UMAPfromPCA500_n_components10_n_neighbors5_min_dist0d9_rs42.csv', \n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_UMAPfromTSVD200_n_components5_n_neighbors5_min_dist0d9_rs42.csv',\n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_UMAPfromPCA500_n_components10_n_neighbors15_min_dist0d9_rs42.csv', \n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/_citeseq_meta_Gender_and_CellType1Hot.csv', \n# '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/citeseq_UMAPfromPCA500_n_components50_n_neighbors15_min_dist0d9_rs42.csv',          \n          #    '/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    print(d.shape)\n    r = np.concatenate( (r,d.values), axis = 1)\nprint(r.shape)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-25T13:56:19.668467Z","iopub.execute_input":"2022-10-25T13:56:19.668894Z","iopub.status.idle":"2022-10-25T13:56:24.857844Z","shell.execute_reply.started":"2022-10-25T13:56:19.668862Z","shell.execute_reply":"2022-10-25T13:56:24.856058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.metrics import r2_score\n\nverbose = 100\n\nmodel = KNeighborsRegressor(n_neighbors=600)\n#, (train_index, test_index )\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\n    X = r # df_cite_train_x.iloc[:,:N_features].values\n    if verbose>= 100:\n        print('X.shape:', X.shape,'Y.shape:', Y.shape)\n\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\n        if verbose >= 100:\n            print('Correlation scores:', np.round(list_scores,4),  'time:', '%.2f'%(time.time() - t1))\n        \n\n    if verbose >= 1:\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-25T13:56:33.047772Z","iopub.execute_input":"2022-10-25T13:56:33.048247Z","iopub.status.idle":"2022-10-25T14:15:01.735387Z","shell.execute_reply.started":"2022-10-25T13:56:33.04821Z","shell.execute_reply":"2022-10-25T14:15:01.73415Z"},"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 = KNeighborsRegressor(n_neighbors=600)\n\nmodel","metadata":{"execution":{"iopub.status.busy":"2022-10-25T14:15:01.737457Z","iopub.execute_input":"2022-10-25T14:15:01.737825Z","iopub.status.idle":"2022-10-25T14:15:01.750391Z","shell.execute_reply.started":"2022-10-25T14:15:01.737791Z","shell.execute_reply":"2022-10-25T14:15:01.749094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nmodel.fit(X[:70988,:], Y)","metadata":{"execution":{"iopub.status.busy":"2022-10-25T14:15:01.751843Z","iopub.execute_input":"2022-10-25T14:15:01.752506Z","iopub.status.idle":"2022-10-25T14:15:01.786214Z","shell.execute_reply.started":"2022-10-25T14:15:01.752473Z","shell.execute_reply":"2022-10-25T14:15:01.785426Z"},"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-25T14:15:01.788334Z","iopub.execute_input":"2022-10-25T14:15:01.789488Z","iopub.status.idle":"2022-10-25T14:17:21.573705Z","shell.execute_reply.started":"2022-10-25T14:15:01.789455Z","shell.execute_reply":"2022-10-25T14:17:21.572527Z"},"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-25T14:17:21.575299Z","iopub.execute_input":"2022-10-25T14:17:21.577114Z","iopub.status.idle":"2022-10-25T14:18:34.197653Z","shell.execute_reply.started":"2022-10-25T14:17:21.577065Z","shell.execute_reply":"2022-10-25T14:18:34.196003Z"},"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-25T14:18:34.200084Z","iopub.execute_input":"2022-10-25T14:18:34.20063Z","iopub.status.idle":"2022-10-25T14:18:34.227884Z","shell.execute_reply.started":"2022-10-25T14:18:34.200581Z","shell.execute_reply":"2022-10-25T14:18:34.22639Z"},"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 = 'CVschemesV31NewFeaturesAndAlfa_multiome_sskknts_KerDrop'\ndf_submission_full.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2022-10-25T14:18:34.23003Z","iopub.execute_input":"2022-10-25T14:18:34.230623Z","iopub.status.idle":"2022-10-25T14:20:50.121405Z","shell.execute_reply.started":"2022-10-25T14:18:34.23057Z","shell.execute_reply":"2022-10-25T14:20:50.120485Z"},"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-25T14:20:50.122815Z","iopub.execute_input":"2022-10-25T14:20:50.124277Z","iopub.status.idle":"2022-10-25T14:20:50.132893Z","shell.execute_reply.started":"2022-10-25T14:20:50.124229Z","shell.execute_reply":"2022-10-25T14:20:50.130873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploration of additional features ","metadata":{}},{"cell_type":"code","source":"%%time\nif additional_feature_analysis_part: \n\n    # Notebook versions 6 and above use prepared features, so we do not need to make PCA/SVD\n    N_features = 100\n    r = df_cite.values[:,:N_features]\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    print(r[:3,:5])\n    print(r.shape)","metadata":{"execution":{"iopub.status.busy":"2022-10-25T14:20:50.136032Z","iopub.execute_input":"2022-10-25T14:20:50.137405Z","iopub.status.idle":"2022-10-25T14:20:50.1501Z","shell.execute_reply.started":"2022-10-25T14:20:50.137321Z","shell.execute_reply":"2022-10-25T14:20:50.148235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if additional_feature_analysis_part: \n    \n    list_f = [ '/kaggle/input/top-100-closest-to-target-genes-with-features/split_by_target/anti-human_CD63_semantic_similarity.csv',\n    '/kaggle/input/top-100-closest-to-target-genes-with-features/split_by_target/anti-human_CD49b_semantic_similarity.csv',\n    '/kaggle/input/top-100-closest-to-target-genes-with-features/split_by_target/anti-human_CD62L_semantic_similarity.csv',\n    '/kaggle/input/top-100-closest-to-target-genes-with-features/split_by_target/anti-human_HLA-DR_semantic_similarity.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)","metadata":{"execution":{"iopub.status.busy":"2022-10-25T14:20:50.154599Z","iopub.execute_input":"2022-10-25T14:20:50.155816Z","iopub.status.idle":"2022-10-25T14:20:53.643654Z","shell.execute_reply.started":"2022-10-25T14:20:50.15574Z","shell.execute_reply":"2022-10-25T14:20:53.642116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nif additional_feature_analysis_part: \n\n    dn = '/kaggle/input/feature-shop-for-multimodal-singlecell-competition/'\n\n    list_f = os.listdir(dn)\n    list_f.remove('_citeseq_cell_ids_only.csv'); list_f = ['_citeseq_cell_ids_only.csv'] +  list_f # Put empty file at the first position\n    # Skip technical and basic  feature files, keep only additional feature candidates \n    for f in ['_citeseq_meta_all_text_also.csv', # '_citeseq_cell_ids_only.csv',\n                'citeseq_train_and_test_TruncatedSVD200_niter7_rs42.csv',\n                'citeseq_train_and_test_TruncatedSVD500.csv',\n                'citeseq_train_and_test_PCA200.csv',\n             'citeseq_train_and_test_PCA500.csv', \n             'citeseq_TRIMAPfromTSVD200_n_components10.csv',\n             'citeseq_UMAPfromPCA500_n_components50_n_neighbors5_min_dist0d1_rs42.csv']: \n            list_f.remove(f)\n    print(len(list_f ))\n\n    list_folds_indices = list_folds_indices_by_days_and_donors_with2holdouts\n    verbose = 10\n\n\n    df_stat_new_feat = pd.DataFrame()\n    cc = -1\n    for f in list_f[:1000]:\n        if f in ['citeseq_UMAPfromPCA500_n_components50_n_neighbors100_min_dist0d9_rs42.csv']: continue # broken file\n        ff = os.path.join(dn, f)\n\n        print(f)\n        try:\n            d = pd.read_csv(ff, index_col = 0)\n            X = np.concatenate( (r,d.values), axis = 1)\n        except:\n            print('Problem with file', f)\n            continue\n\n        cc += 1\n        df_stat_new_feat.loc[cc,'Features'] = f\n\n        df_stat_new_feat.loc[cc,'N_features'] = X.shape[1]\n        if verbose >= 100:\n            print('X.shape:', X.shape,'Y.shape:', Y.shape)\n\n        model = KNeighborsRegressor(n_neighbors=600)\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            if verbose >= 100:\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            if verbose >= 100:\n                print('Correlation scores:', np.round(list_scores,4),  'time:', '%.2f'%(time.time() - t1))\n\n        df_descr = df_fold_score_stat.describe(percentiles=[] )\n        for ii in range(df_descr.shape[1]):\n            col = df_descr.columns[ii]\n            df_stat_new_feat.loc[cc,col] = df_descr.iloc[1,ii]\n        df_stat_new_feat.loc[cc, 'Time'] = np.round( (time.time() - t0), 2)      \n        if verbose >= 100:\n            display(df_fold_score_stat)\n            display(df_descr.iloc[1:,:])\n        if verbose >= 10:\n            display(df_stat_new_feat.tail(1))\n\n\n    df_stat_new_feat.sort_values( df_stat_new_feat.columns[2], ascending = False ) ","metadata":{"execution":{"iopub.status.busy":"2022-10-25T14:20:53.64567Z","iopub.execute_input":"2022-10-25T14:20:53.64662Z","iopub.status.idle":"2022-10-25T14:26:15.877389Z","shell.execute_reply.started":"2022-10-25T14:20:53.646558Z","shell.execute_reply":"2022-10-25T14:26:15.875849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if additional_feature_analysis_part: \n    df_sorted_stat =  df_stat_new_feat.sort_values( df_stat_new_feat.columns[2], ascending = False )\n\n    print( list( df_sorted_stat['Features'])[:20] )\n\n    print( list(df_sorted_stat['Score Test1'].values - df_stat_new_feat['Score Test1'].iat[0])[:20] )\n\n    display(df_stat_new_feat.describe())\n    display(df_stat_new_feat.head(1) )\n    display( df_stat_new_feat.sort_values( df_stat_new_feat.columns[2], ascending = False ).head(60) )\n    display( df_stat_new_feat.sort_values( df_stat_new_feat.columns[2], ascending = False ).tail(60) )\n\n\n    plt.figure(figsize = (20,6))\n    df_stat_new_feat['Score Test1'].plot()\n    plt.plot(df_stat_new_feat['Score Test1'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    df_stat_new_feat['Score Test2'].plot()\n    plt.plot(df_stat_new_feat['Score Test2'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    plt.legend(); plt.show()\n\n    plt.figure(figsize = (20,6))\n    df_stat_new_feat['Score Test3'].plot()\n    plt.plot(df_stat_new_feat['Score Test3'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    df_stat_new_feat['Score Test4'].plot()\n    plt.plot(df_stat_new_feat['Score Test4'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    plt.legend(); plt.show()\n\n    plt.figure(figsize = (20,6))\n    df_stat_new_feat['Score Train'].plot()\n    plt.plot(df_stat_new_feat['Score Train'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    plt.legend(); plt.show()\n\n\n    plt.figure(figsize = (20,6))\n    df_stat_new_feat['R2 Score Test1'].plot()\n    plt.plot(df_stat_new_feat['R2 Score Test1'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    df_stat_new_feat['R2 Score Test2'].plot()\n    plt.plot(df_stat_new_feat['R2 Score Test2'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    plt.legend(); plt.show()\n\n    plt.figure(figsize = (20,6))\n    df_stat_new_feat['R2 Score Test3'].plot()\n    plt.plot(df_stat_new_feat['R2 Score Test3'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    df_stat_new_feat['R2 Score Test4'].plot()\n    plt.plot(df_stat_new_feat['R2 Score Test4'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    plt.legend(); plt.show()\n\n    plt.figure(figsize = (20,6))\n    df_stat_new_feat['R2 Score Train'].plot()\n    plt.plot(df_stat_new_feat['R2 Score Train'].iat[0]*np.ones(df_stat_new_feat.shape[0] ) )\n    plt.legend(); plt.show()\n\n    display( df_stat_new_feat.corr() )\n    df_stat_new_feat.to_csv('df_stat_new_feat.csv')","metadata":{"execution":{"iopub.status.busy":"2022-10-25T14:26:15.888158Z","iopub.execute_input":"2022-10-25T14:26:15.888929Z","iopub.status.idle":"2022-10-25T14:26:16.004177Z","shell.execute_reply.started":"2022-10-25T14:26:15.888897Z","shell.execute_reply":"2022-10-25T14:26:16.00213Z"},"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-25T14:26:15.878965Z","iopub.execute_input":"2022-10-25T14:26:15.879408Z","iopub.status.idle":"2022-10-25T14:26:15.886162Z","shell.execute_reply.started":"2022-10-25T14:26:15.879359Z","shell.execute_reply":"2022-10-25T14:26:15.88491Z"},"trusted":true},"execution_count":null,"outputs":[]}]}