{"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":"# Still under construction \n\n# What is about ?\n\nLearn ML using CAFA5 as an example. Educational notebook 03. \n\nSee the previous notebooks: https://www.kaggle.com/code/alexandervc/pytorch-01-basics for basics of pytorch - metrics, setting neural networks, training neural networks in the most simple setup.  See also https://www.kaggle.com/code/alexandervc/pytorch-2-cv-focal-sophia-etc for more advanced topics - optimizers, custom losses, etc .\n\n\nHere we  learn how to do some additional things: explore models from differnt frameworks - Pytorch, Keras, Sklearn, blend them , ....\n\n\nWe also extend options realted to the CAFA5 challenge:\n\n    + CAFA F1 metric computation (local)\n    + concatenate several embeddings available: t5, esm2 with different sizes  \n    + todo: more models\n    + todo: preprocessings - add features obtained by pca, ica, umap, etc ... - on-fly\n    \n\n#### Educational task:\n\nChange models to improve the score . \n\n#### Links and Acknowledgements:\n\nPlease upvote useful works: \n\nhttps://www.kaggle.com/code/andreylalaley/pytorch-cafa-5-prediction?scriptVersionId=138595845\nLB 0.52875 Andrey Shevtsov pytorch model with skip connections and Sophia non-standard optimizer\n\nhttps://www.kaggle.com/code/sergeifironov/validate-ridge  - Sergei Fironov - CAFA5 F1 (weighted) metric - ported to Kaggle environment\n\nhttps://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/420241 - Anton Vakhrushev - correcting error in the initial code of the metric computation\n\nhttps://www.kaggle.com/code/simonveitner/simple-mlp - SIMON VEITNER -  Keras-\"Simple MLP\"\n\nhttps://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/406168 - Andrey Shevtsov - datasets with esm2 embeddings\n\nhttps://www.kaggle.com/datasets/sergeifironov/t5embeds - Sergei Fironov - dataset with t5  embeddings\n\nhttps://www.kaggle.com/code/szabo7zoltan/combineembeddings - Zoltan - demonstrated effeciency of combining embeddings, and several other ideas \n\nhttps://www.kaggle.com/datasets/henriupton/protbert-embeddings-for-cafa5   - HENRI UPTON -      ProtBERT Vector Embeddings for Protein Sequences\n\nhttps://www.kaggle.com/code/geraseva/protbert-embedding - Liza Geraseva - protbert embeddings \n\nhttps://www.kaggle.com/datasets/viktorfairuschin/cafa-5-ems-2-embeddings-numpy - VIKTOR FAIRUSCHIN - esm2-size 1280  (Pay attention order of samples is different from the standard one - reoder before using). \n\nhttps://www.kaggle.com/datasets/daehunbae/cafa5-esm-2-3b-various-poolings - DAEHUN BAE - esm2-size2560 with 3 different poolings\n\nhttps://www.kaggle.com/datasets/visualcomments/uniprot-description-only - Anton Kostin - embeddings based on text descriptions of proteins , using Sentence Bert  \n\n\ntodo: add more \n\n\n\n\n#### Main Versions and updates: \n\n    244 Added:  protbert-embedding - by Liza Geraseva https://www.kaggle.com/code/geraseva/protbert-embedding/data \n    239 Added: dict_labels_to_consider_by_subontology e.g. {'BPO':1100, 'MFO':450, 'CCO':300} - how many top targets to take from each subontology\n    217 Option - compute cafa_f1 on downsamle  - mode_downsample_validation_for_cafa_f1 = e.g. 'top_300' # 'random_1000'\n    201 namepostfix option - to add to str_model_id\n    199 mode_downsample_train_default - option to downsample only train part - while OOF will be on full data, \n        supported:  '43k','random_subsample_percent_XXX'  \n    177 Option to calculate stat for individual models: flag_compute_stat_for_each_model flag_compute_cafa_f1_for_each_model, \n        Plots for analysis of the individual models and blend \n    169 Protbert embeddings (thxs: Henry Upton ) added to options \n    167-168 special modeling code for Keras models - trying to understand memory leaking with Keras \n    150 LB 0.52414 - 12 times self-blend for Andrey's pytorch NN with Sophia optimizer , with 20 epochs(!) each. \n    77 \"quick_mode\" - downsampling for parameter tuning \"very_quick_mode\" - crude downsampling just for checks that code does not have bugs\n    71-74 Added: option to have different features in each configuration \n    55 Added:  flag_save_numpy_Y_pred_..,flag_save_numpy_Y_submit - options to save predictions in numpy format for possible further blend outside the notebook  \n    39 Added config: 'custom_model_fit_function_name': e.g.: 'model_fit_pytorch_Sophia_Andrey1' - option to customize model train \n    27 LB 0.52189 CV: 0.5482, 0.96668 5*KMLP t5+esm2s1260 8h50min\n    26 - option to concatenate embeddings/other-features is added \n    25 - CV: 0.52366, 0.959 5*KMLP 'esm2S640' 7р52min\n    18 - local CAFA-F1 calculation added\n    3 option for oof computation and aucroc score for oof added \n    2 Simple version with many model training and blending submission predictions. No OOF-predictions, no local scoring . \n\n#### Others versions with some results and remarks: \n    \n    231 LB 0.51641 6.5h Cpu trained on 43k, selfblend=10, model = kMLP 'Layers': [800,1500],'Dropouts': [0.2]\n        Keras crashes by RAM, only 1 out many runs can produce results - run 5 may be 1 will be successful. Reasons still not clear. \n    215 LB 0.52931 1h GPU trained on 43k, selfblend=12, model = pMLP-Andrey(Sophia,ResCon), 20 epochs, no scoring saved.  \n    209 Full-mode 43min (mainly F1-score - 32min, 8min - save submission ): \n        Ridge10 t5+esm2-1280\t 0.929654\t0.518779\t0.573918\t0.380896\t0.601523\t0.216\t0.134\t0.174\t326.64\n    208 Full-mode 43k 1h33min Ridge10\t\t0.920402\t0.503298\t0.555992\t0.359805\t0.594098\n    150 LB 0.52414 - 12 times self-blend for Andrey's pytorch NN with Sophia optimizer , with 20 epochs(!) each. \n    78-... Modeling in quick_mode: Self-blend gives uplift even upto 40 (or more) times  - not clear is not only due to down sampling or in general\n    78-... Modeling in quick_mode: Epochs 20 is better than 10,30,40 for pMLP-Andrey . See version 150 for submit - LB 0.52414 \n    78-... Modeling in quick_mode: Seems 1000,2000Dr0.2 is top for Keras 2-Layer MLP \n    Later change threshold to 0.1 \n    Measurments: about 30 min for CAFA-F1 on each fold - so 2.5h for folds -- for threshold 0.01\n    Seems results are quite stable - two runs with same params - does not change much the scores : \n    47 3*skMLP 6h24min F1w 0.52884 AUC 0.94253 F1 CCO  0.59000 F1 BPO 0.4 F1 MFO 0.596 thres  CCO 0.206 thres BPO 0.138 thres MFO 0.17\n    46 3*skMLP 7h19min F1w 0.52864 AUC 0.94284 F1 CCO  0.58987 F1 BPO 0.4 F1 MFO 0.5957 thres CCO 0.212 thres BPO 0.142 thres MFO 0.172","metadata":{}},{"cell_type":"markdown","source":"# Key params \n\nNotebook allows to train several different models in the same framework and blend them.\n\nlist_main_config_model_feature_etc - list containing specification of each model -  key parameter for the notebook\n","metadata":{}},{"cell_type":"code","source":"mode_loc =   'full_run' #  'very_quick_run' #     'full'# 'quick_run' #      'full' #         'quick_run' #              \n\nif mode_loc == 'very_quick_run': # CRUDE downsampling. For debug and testing \n    \n    # Crude downsampling - just for speed run/check - about 2-3 minutes for run \n    n_samples_to_consider =  1000#    downsampling - might be useful for debug:  small number - fast run \n    #n_labels_to_consider = 10 # Up to 31466 but more than 3000-5000 may crash RAM \n    dict_labels_to_consider_by_subontology = { 'BPO':5, 'MFO':5, 'CCO':5} # how many targets to take from each subontology\n    n_folds_to_process = 1 # we can only use 1 fold - for fast checks\n    \nelif mode_loc == 'quick_run': # Downsampling. You can tune params for the models using that mode. Later compute the best model in the \"full\" mode.\n    \n    n_samples_to_consider =   20_000 #  1000 #   50_000#   142246 #    downsampling - might be useful for debug:  small number - fast run \n    #n_labels_to_consider = 1000 # Up to 31466 but more than 3000-5000 may crash RAM \n    dict_labels_to_consider_by_subontology = {'BPO':1100, 'MFO':450, 'CCO':300} # how many targets to take from each subontology\n    n_folds_to_process = 2 # can reduce number of folds to speed-up - set 1,2,3 .. , , here 100 folds does NOT mean 100 folds - it just will be effectively clipped down the same number as in loaded folds file\n\nelse: # Full mode for final computations\n    n_samples_to_consider =   142246 #  1000 #   50_000#   142246 #    downsampling - might be useful for debug:  small number - fast run \n    #n_labels_to_consider = 2000 # Up to 31466 but more than 3000-5000 may crash RAM \n    dict_labels_to_consider_by_subontology = {'BPO':1100, 'MFO':450, 'CCO':300} # how many targets to take from each subontology\n    n_folds_to_process = 5 # can reduce number of folds to speed-up - set 1,2,3 .. , , here 100 folds does NOT mean 100 folds - it just will be effectively clipped down the same number as in loaded folds file\n\n    \n######################################### Set the models, their params, etc  ########################################################\n######################################### Set the models, their params, etc  ########################################################\nlist_main_config_model_feature_etc = [] #     \n# for list_features_id in (['t5'],['esm2S2560'],['esm2S1280'] , ['esm2S640'], ['esm2S480'], ['esm2S320'] ):\n#     cfg = {'model': {'id':'Ridge' }  , 'list_features_id': list_features_id }\n#     list_main_config_model_feature_etc.append( cfg )\n# for list_features_id in (['t5','esm2S2560'], ['t5','esm2S1280'] , ['t5', 'esm2S640', 'esm2S480']  ):\n#     cfg = {'model': {'n_selfblend':2, 'id':'KMLP_simple' , 'epochs':15 ,   'batch_size':128 , 'verbose':0} , 'list_features_id': list_features_id }\n#     list_main_config_model_feature_etc.append( cfg )\n\n# cfg = {'model': {'id':'Ridge' }  }\n# list_main_config_model_feature_etc.append( cfg )\n\n# cfg = {'model': {'id':'pMLP_simple1', 'Layers':[500,1000] }  }\n# list_main_config_model_feature_etc.append( cfg )\n\n# for dr in [None,  0.1, 0.15, 0.2,  0.25,0.3,  0.35, 0.4,0.45 ]:\n#     #Default:  'epochs': 10, 'batch_size': 128, 'LR':0.001 \n#     cfg = {'model': {'n_selfblend':1, 'id':'pMLP_Configurable1', 'Layers':[2000,1500], 'BatchNormalizations': [True ],'Dropouts':[0.35 ] ,\n#                     'epochs': 15, 'batch_size': 128, 'LR':0.001 }  }\n#     list_main_config_model_feature_etc.append( cfg )\n\n##########################\n\"\"\"( ['t5', 'esm2S1280' ], ['t5', 'esm2S1280','esm2S320' ],   ['t5', 'esm2S1280','protbert' ] ,  ['t5' ] ,  ['esm2S640' ] ,\n     ['esm2S1280' ] ,  ['esm2S640' ]  , ['esm2S480' ]  , ['esm2S320' ]  ):\"\"\"\n\n##########################\nfor list_features_id_loc in  ( ['t5', 'esm2S1280' ], ['t5', 'esm2S1280','esm2S320' ],   ['t5', 'esm2S1280','protbert' ] ,  ['t5' ], ['esm2S1280' ] ):\n    cfg = {'model': {'n_selfblend':1, 'id':'pMLP_Configurable1', 'Layers':[2000,1500], 'BatchNormalizations': [True],'Dropouts':[0.35 ] ,\n                    'epochs': 20, 'batch_size': 128, 'LR':0.001,  'namepostfix': str( list_features_id_loc ) } ,  'list_features_id':list_features_id_loc }\n    list_main_config_model_feature_etc.append( cfg )\n\n# for L1 in [500,1000,1500,2000]:\n#     for L2 in [500,1000,1500,2000]:\n#         cfg = {'model': {'id':'pMLP_simple1', 'Layers':[L1,L2] }  }\n#         list_main_config_model_feature_etc.append( cfg )\n\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' t5 esm2-1280,320' }, 'list_features_id':  ['t5', 'esm2S1280','esm2S320' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' t5 esm2-1280 protbert' }, 'list_features_id':  ['t5', 'esm2S1280','protbert' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' t5' }, 'list_features_id':  ['t5' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' esm2S640' }, 'list_features_id':  ['esm2S640' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' esm2S1280' }, 'list_features_id':  ['esm2S1280' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' esm2S640' }, 'list_features_id':  ['esm2S640' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' esm2S480' }, 'list_features_id':  ['esm2S480' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' esm2S320' }, 'list_features_id':  ['esm2S320' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n\n# cfg = {'model': {'id':'Ridge', 'namepostfix':'protbert' } , 'list_features_id':  ['protbert' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':'protbert_2' } , 'list_features_id':  ['protbert_2' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n\n# cfg = {'model': {'id':'Ridge' ,  'n_selfblend':1, 'mode_downsample_train':'43k', 'namepostfix':' 43k t5+esm2S1280'}  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge' ,  'n_selfblend':1, 'mode_downsample_train':'random_subsample_percent_30' , 'namepostfix':' RS30 t5+esm2S1280'}  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' 43k t5' } , 'list_features_id':  ['t5' ] }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge' ,  'n_selfblend':1, 'mode_downsample_train':'43k', 'namepostfix':' 43k t5'}, 'list_features_id':  ['t5' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge' ,  'n_selfblend':1, 'mode_downsample_train':'random_subsample_percent_30' , 'namepostfix':' RS30 t5'} , 'list_features_id':  ['t5' ] }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge', 'namepostfix':' 43k esm2S1280' } , 'list_features_id':  ['esm2S1280' ] }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge' ,  'n_selfblend':1, 'mode_downsample_train':'43k', 'namepostfix':' 43k esm2S1280'}, 'list_features_id':  ['esm2S1280' ]  }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge' ,  'n_selfblend':1, 'mode_downsample_train':'random_subsample_percent_30' , 'namepostfix':' RS30 esm2S1280'} , \n#        'list_features_id':  ['esm2S1280' ] }\n# list_main_config_model_feature_etc.append( cfg )\n\n# cfg = {'model': {'id':'Ridge' }  , 'list_features_id': ['esm2S320'] }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge' }  , 'list_features_id': ['protbert'] }\n# list_main_config_model_feature_etc.append( cfg )\n# cfg = {'model': {'id':'Ridge' }  , 'list_features_id': ['t5'] }\n# list_main_config_model_feature_etc.append( cfg )\n\n# cfg2 = {'model': {'id':'pMLP_Andrey' , 'n_selfblend':2, 'epochs':20, 'custom_model_fit_function_name':'model_fit_pytorch_Sophia_Andrey1' } } # Pytorch MLP model\n# list_main_config_model_feature_etc.append( cfg2 )\n# cfg2b = {'model': {'id':'pMLP_Andrey' , 'n_selfblend':12, 'epochs':15, 'custom_model_fit_function_name':'model_fit_pytorch_Sophia_Andrey1' } } # Pytorch MLP model\n# list_main_config_model_feature_etc.append( cfg2b )\n\n#cfg3 = {'model': {'n_selfblend':10 , 'id':'skMLP', 'Layers': [ 500,800] ,'LR':  0.001 , 'alpha':1e-4, 'max_iter':500 } } # Sklearn MLP model\n#list_main_config_model_feature_etc.append( cfg3 )\n# cfg4 = {'model': {'n_selfblend':1, 'id':'KMLP_simple' , 'epochs':15 ,   'batch_size':128 , 'verbose':0} } # Simple Keras MLP\n# list_main_config_model_feature_etc.append( cfg4 )\n# cfg5 = {'model': {'n_selfblend':1, 'id':'KMLP', 'Layers': [400,600],'Dropouts': [0.1 ],  'BatchNormalizations': [] , 'epochs':15 ,   'batch_size':128 , 'verbose':0 } } \n# #  Keras MLP - with adjustable params \n# list_main_config_model_feature_etc.append( cfg5 )\n# cfg6 = {'model': {'Keras': True, 'n_selfblend':10, 'id':'KMLP', 'Layers': [1000,2000],'Dropouts': [0.2],  'BatchNormalizations': [False] , 'epochs':15 ,   'batch_size':128 , 'verbose':0 } } \n# # Keras MLP - with adjustable params \n# list_main_config_model_feature_etc.append( cfg6 )\n\n# list_layers_sizes = [500,1000,1500]\n# for L1 in list_layers_sizes:\n#     for L2 in list_layers_sizes:\n# #         for L3 in list_layers_sizes:\n#         cfg6b = {'model': {'Keras': True, 'n_selfblend':1, 'id':'KMLP', 'Layers': [L1,L2],'Dropouts': [0.3],  'BatchNormalizations': [False] , 'epochs':15 ,   'batch_size':128 , 'verbose':0 } } \n#         # Keras MLP - with adjustable params \n#         list_main_config_model_feature_etc.append( cfg6b )\n\n# list_layers_sizes = [500,1000,1500]\n# for L1 in list_layers_sizes:\n#     for L2 in list_layers_sizes:\n# #         for L3 in list_layers_sizes:\n#         cfg6b = {'model': {'Keras': True, 'n_selfblend':1, 'id':'KMLP', 'Layers': [L1,L2],'Dropouts': [0.3],  'BatchNormalizations': [False] , 'epochs':15 ,   'batch_size':128 , 'verbose':0 } } \n#         # Keras MLP - with adjustable params \n#         list_main_config_model_feature_etc.append( cfg6b )\n\n\n\n#########################################  Feature sets  ########################################################\n\n# Features: \n# default_features =  [ {'id':'t5'}, {'id':'esm2S2560','preproc':'PCA100'} ] # 'esm2S2560' #  'esm2S1280' #  'esm2S320' #   't5'#  \n#list_features_id =  ['t5', 'esm2S1280', 'esm2S320', 'esm2S640']   # 'esm2S640'  #  'esm2S480' # 't5'#   'esm2S2560' #  'esm2S1280' #  'esm2S320' #   \n# list_features_id =  ['protbert'  ] #  ['protbert']#  ] \nlist_features_id =  ['t5', 'esm2S1280' ] #  ['protbert']#  ] \n# Concatenate several embeddings (other features)\n# !!! Pay attention proteins ids should be the same as in all the files !!!!!!!!!!!!!!!!!!\n\n#########################################  Compute/not/when the scores  ########################################################\n# Params to control scoring computations.\n# Computation of CAFA-F1 is quite time and RAM consuming. So we need to have options to restrict its computation only to desired cases. \n\n\nflag_compute_oof_predictions = True # Compute local OOF predictions.  Necessary to compute local CV scores/other statistics        \n\n# Compute statistics on OOF for each model \nflag_compute_stat_for_each_model = ( False ) and (flag_compute_oof_predictions)\nflag_compute_cafa_f1_for_each_model = ( True ) and ( flag_compute_stat_for_each_model )  and ( flag_compute_oof_predictions )\n\n# Compute statistics on blend of OOF each blend  \nflag_compute_each_blend_stat = ( False) and (flag_compute_oof_predictions)\nflag_compute_cafa_f1_for_each_blend = ( True ) and (flag_compute_each_blend_stat) and (flag_compute_oof_predictions)\n\n# Compute stat on final blend   \nflag_compute_final_model_stat = ( True ) and (flag_compute_oof_predictions)\n\n\n#########################################  Furthter params   ########################################################\n\nmode_downsample_train_default = '43k' ###'43k'  #  None #    'random_subsample_percent_90' #  None # \n\nmode_downsample_validation_for_cafa_f1 = 'top_5000' #  'random_1000'#  \n\ncutoff_threshold_low = 0.1 # prediction < cutoff_threshold_low will be set to zero (i.e. no need to save to submission file)\n\nmode_submit = True # True # Compute prediction for submission part and prepare submission file in required CAFA5 format. Set to False if you only interested in local score\n\n###  Save/not/what predictions   #####\nflag_save_final_submit_file = (True ) and ( mode_submit ) # Prepare and save final submission txt file. \n#  mode_submit - controls two things - 1) computation Y_submit  (and it save as numpy array) 2) final submission file, and \n\n# Set these flags to False if do not plan to blend current prediction with other models outside the notebook:\nflag_save_numpy_Y_pred_oof_blend = (True)  and (flag_compute_oof_predictions)  # Save  Y_pred_oof_blend in numpy format \"npy\" for possible further blend outside current notebook \nflag_save_numpy_Y_submit = (True) and (mode_submit) # Save Y_submit matrix in numpy  format \"npy\" for possible further blend outside current notebook \n\n\nRANDOM_SEED = None # Fix or Not random seed \n\n\n#########################################  Change some params  ########################################################\n\nif mode_loc == 'quick_run': # Downsampling. You can tune params for the models using that mode. Later \n    flag_save_numpy_Y_pred_oof_blend = False # Save  Y_pred_oof_blend in numpy format \"npy\" for possible further blend outside current notebook \n    flag_save_numpy_Y_submit = False # Save Y_submit matrix in numpy  format \"npy\" for possible further blend outside current notebook \n    mode_submit = False # Compute prediction for submission part and prepare submission file in required CAFA5 format. Set to False if you only interested in local score\n\n\n#########################################  Information string  ########################################################\n\nstr_id = str(list_features_id)+'_Y'+str(dict_labels_to_consider_by_subontology)\nstr_id += '_S'+str(n_samples_to_consider)\nstr_id += '_CUT'+str(cutoff_threshold_low)\nprint(str_id, len(str_id))\nprint(str_id[:48])","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:44:34.971839Z","iopub.execute_input":"2023-08-12T19:44:34.972391Z","iopub.status.idle":"2023-08-12T19:44:35.032263Z","shell.execute_reply.started":"2023-08-12T19:44:34.972343Z","shell.execute_reply":"2023-08-12T19:44:35.030818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"logs_file_path = 'logs.txt'","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:44:35.035083Z","iopub.execute_input":"2023-08-12T19:44:35.035548Z","iopub.status.idle":"2023-08-12T19:44:35.043123Z","shell.execute_reply.started":"2023-08-12T19:44:35.035506Z","shell.execute_reply":"2023-08-12T19:44:35.041808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for cfg in list_main_config_model_feature_etc:\n    print(cfg)","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:44:35.045192Z","iopub.execute_input":"2023-08-12T19:44:35.046052Z","iopub.status.idle":"2023-08-12T19:44:35.055832Z","shell.execute_reply.started":"2023-08-12T19:44:35.046007Z","shell.execute_reply":"2023-08-12T19:44:35.054385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparations \n\nTechnical installs/imports","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport time\nt0start = time.time()\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:44:35.057674Z","iopub.execute_input":"2023-08-12T19:44:35.058085Z","iopub.status.idle":"2023-08-12T19:44:35.105929Z","shell.execute_reply.started":"2023-08-12T19:44:35.058055Z","shell.execute_reply":"2023-08-12T19:44:35.104702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import psutil\nimport datetime\n\ndef get_available_ram():\n    virtual_memory = psutil.virtual_memory()\n    available_ram = virtual_memory.available\n    return available_ram\ndef log_available_ram( str_for_logging_optional = None):\n    try:\n        virtual_memory = psutil.virtual_memory()\n        available_ram_bytes = virtual_memory.available\n        # Convert bytes to other units if needed (e.g., megabytes, gigabytes)\n        available_ram_megabytes = available_ram_bytes / (1024 ** 2)\n        available_ram_gigabytes = available_ram_bytes / (1024 ** 3)\n\n    #     print(f\"Available RAM: {available_ram_bytes} bytes\")\n    #     print(f\"Available RAM: {available_ram_megabytes:.2f} MB\")\n        if str_for_logging_optional is not None:\n            print(str_for_logging_optional)\n        current_datetime = datetime.datetime.now()\n        str1 = f\"Available RAM: {available_ram_gigabytes:.2f} G  Current datetime: {current_datetime}\"\n        print( str1 ) \n        \n        with open(logs_file_path, 'a') as file:\n            if str_for_logging_optional is not None:\n                file.write(str_for_logging_optional + '\\n')\n            file.write(str1 + '\\n')\n        # print(\"Data appended successfully.\")\n    except Exception as e:\n        print(f\"Error while appending data: {e}\")        \n    \n#     return available_ram\nlog_available_ram('On start')\n\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:44:35.109077Z","iopub.execute_input":"2023-08-12T19:44:35.109414Z","iopub.status.idle":"2023-08-12T19:44:35.123447Z","shell.execute_reply.started":"2023-08-12T19:44:35.109386Z","shell.execute_reply":"2023-08-12T19:44:35.121986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture \n!pip install torchmetrics\n\nfrom torchmetrics import AUROC as torch_AUCROC\nfrom torchmetrics import F1Score as torch_F1Score\n\n# from torchmetrics import AUROC,F1Score\n# from torchmetrics.classification import BinaryF1Score\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:44:35.124984Z","iopub.execute_input":"2023-08-12T19:44:35.125324Z","iopub.status.idle":"2023-08-12T19:44:48.971643Z","shell.execute_reply.started":"2023-08-12T19:44:35.125297Z","shell.execute_reply":"2023-08-12T19:44:48.969934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture\n!pip install torchsummary\nfrom torchsummary import summary as torchsummary\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:44:48.973939Z","iopub.execute_input":"2023-08-12T19:44:48.974461Z","iopub.status.idle":"2023-08-12T19:45:02.443742Z","shell.execute_reply.started":"2023-08-12T19:44:48.974408Z","shell.execute_reply":"2023-08-12T19:45:02.441892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## https://www.kaggle.com/code/alexandervc/baseline-multilabel-to-multitarget-binary#Load-train-features---precalculated-embeddings-for-the-proteins\nimport os\nimport gc\nfrom sklearn.model_selection import train_test_split\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\ntqdm.pandas()\nimport torch\nimport warnings\nwarnings.filterwarnings('ignore')\nfrom sklearn.model_selection import train_test_split\nimport torch.nn.functional as F\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, TensorDataset , DataLoader\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-08-12T19:45:02.445623Z","iopub.execute_input":"2023-08-12T19:45:02.446057Z","iopub.status.idle":"2023-08-12T19:45:02.460809Z","shell.execute_reply.started":"2023-08-12T19:45:02.446016Z","shell.execute_reply":"2023-08-12T19:45:02.459543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport random \ndef seed_all(RANDOM_SEED):\n    if RANDOM_SEED is not None: \n        try:\n            SEED = RANDOM_SEED\n            random.seed(SEED)\n            np.random.seed(SEED)\n            torch.manual_seed(SEED)\n            torch.cuda.manual_seed_all(SEED)\n            torch.backends.cudnn.deterministic = True\n            torch.backends.cudnn.benchmark = False\n\n        except Exception as e:\n            print(f\"Exception: {e}\")\n            \nseed_all(RANDOM_SEED)            ","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:02.462854Z","iopub.execute_input":"2023-08-12T19:45:02.463307Z","iopub.status.idle":"2023-08-12T19:45:02.478097Z","shell.execute_reply.started":"2023-08-12T19:45:02.463273Z","shell.execute_reply":"2023-08-12T19:45:02.476819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Optimizer \"Sophia\" sometimes better than Adam","metadata":{}},{"cell_type":"code","source":"%%capture\n!git clone https://github.com/kyegomez/Sophia.git\n! python Sophia/setup.py install\n!rm Sophia/Sophia/__init__.py\nfrom Sophia.Sophia.Sophia import SophiaG","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:02.480296Z","iopub.execute_input":"2023-08-12T19:45:02.48094Z","iopub.status.idle":"2023-08-12T19:45:09.47122Z","shell.execute_reply.started":"2023-08-12T19:45:02.480897Z","shell.execute_reply":"2023-08-12T19:45:09.469403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Choose device - GPU or CPU and assign model to it","metadata":{}},{"cell_type":"code","source":"device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nprint(device)","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:09.473305Z","iopub.execute_input":"2023-08-12T19:45:09.473698Z","iopub.status.idle":"2023-08-12T19:45:09.483037Z","shell.execute_reply.started":"2023-08-12T19:45:09.473658Z","shell.execute_reply":"2023-08-12T19:45:09.481215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load/prepare features: X \n\nLoad precalulcated features - embeddings and targets transformed in 0,1 multi-target task.\n\nWe use already computed emebdding for the proteins sequences - thanks to Andrey Shevtsov for sharing the embeddings by esm2-model: \n( https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/406168 - please upvote his work ).\n\nWe load targets matrix \"Y\" which contains NOT all the CAFA5 targets but only top 1499. That is quite enough to get not so bad score. \nNote: original input file text describing targets have been trasnformed to np.array Y.\nSee the notebook https://www.kaggle.com/code/alexandervc/baseline-multilabel-to-multitarget-binary for details. \n","metadata":{}},{"cell_type":"code","source":"%%time \n\ndef get_paths_to_features(features_id):\n    # Thanks to Andrey Shevtsov https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/406168\n    # Sergei Fironov: https://www.kaggle.com/datasets/sergeifironov/t5embeds\n    # Please upvote ! \n    \n    if features_id == 'esm2S2560':\n        fn_X = '/kaggle/input/4637427/train_embeds_esm2_t36_3B_UR50D.npy'\n        fn_protein_ids = '/kaggle/input/4637427/train_ids_esm2_t36_3B_UR50D.npy'\n        fn_X_submit = '/kaggle/input/4637427/test_embeds_esm2_t36_3B_UR50D.npy'\n        fn_submit_protein_ids = '/kaggle/input/4637427/test_ids_esm2_t36_3B_UR50D.npy'\n        \n    elif features_id == 'esm2S1280':\n        fn_X = '/kaggle/input/23468234/train_embeds_esm2_t33_650M_UR50D.npy'\n        fn_protein_ids = '/kaggle/input/23468234/train_ids_esm2_t33_650M_UR50D.npy'\n        fn_X_submit = '/kaggle/input/23468234/test_embeds_esm2_t33_650M_UR50D.npy'\n        fn_submit_protein_ids = '/kaggle/input/23468234/test_ids_esm2_t33_650M_UR50D.npy'\n        \n    elif features_id == 't5':\n        fn_X = '/kaggle/input/cafa5-features-etc/T5_train_embeds_float32.npy'\n        fn_protein_ids = '/kaggle/input/cafa5-features-etc/train_ids.npy'\n        fn_X_submit = '/kaggle/input/cafa5-features-etc/T5_test_embeds_float32.npy'\n        fn_submit_protein_ids = '/kaggle/input/cafa5-features-etc/test_ids.npy'\n        \n    elif features_id == 'esm2S320':\n        fn_X = '/kaggle/input/315701375/train_embeds_esm2_t6_8M_UR50D.npy'\n        fn_protein_ids = '/kaggle/input/315701375/train_ids_esm2_t6_8M_UR50D.npy'\n        fn_X_submit = '/kaggle/input/315701375/test_embeds_esm2_t6_8M_UR50D.npy'\n        fn_submit_protein_ids = '/kaggle/input/315701375/test_ids_esm2_t6_8M_UR50D.npy'\n        \n    elif features_id == 'esm2S640':\n        dn = '/kaggle/input/8947923/'\n        fn_X = os.path.join(  dn , 'train_embeds_esm2_t30_150M_UR50D.npy' )\n        fn_protein_ids = os.path.join(  dn , 'train_ids_esm2_t30_150M_UR50D.npy' )\n        fn_X_submit = os.path.join(  dn , 'test_embeds_esm2_t30_150M_UR50D.npy' )\n        fn_submit_protein_ids = os.path.join(  dn ,  'test_ids_esm2_t30_150M_UR50D.npy' )\n        \n    elif features_id == 'esm2S480':\n        dn = '/kaggle/input/3023750/'\n        fn_X = os.path.join(  dn , 'train_embeds_esm2_t12_35M_UR50D.npy' )\n        fn_protein_ids = os.path.join(  dn , 'train_ids_esm2_t12_35M_UR50D.npy' )\n        fn_X_submit = os.path.join(  dn , 'test_embeds_esm2_t12_35M_UR50D.npy' )\n        fn_submit_protein_ids = os.path.join(  dn ,  'test_ids_esm2_t12_35M_UR50D.npy' )\n        \n    elif features_id == 'protbert':\n        # HENRI UPTON https://www.kaggle.com/datasets/henriupton/protbert-embeddings-for-cafa5        \n        dn = '/kaggle/input/protbert-embeddings-for-cafa5/'\n        fn_X = os.path.join(  dn , 'train_embeddings.npy' )\n        fn_protein_ids = os.path.join(  dn , 'train_ids.npy' )\n        fn_X_submit = os.path.join(  dn , 'test_embeddings.npy' )\n        fn_submit_protein_ids = os.path.join(  dn ,  'test_ids.npy' )\n\n    elif features_id == 'protbert_2':\n        # Liza Geraseva https://www.kaggle.com/code/geraseva/protbert-embedding/data       \n        dn = '/kaggle/input/protbert-embedding/'\n        fn_X = os.path.join(  dn , 'train_embeddings.npy' )\n        fn_protein_ids = os.path.join(  dn , 'train_ids.npy' )\n        fn_X_submit = os.path.join(  dn , 'test_embeddings.npy' )\n        fn_submit_protein_ids = os.path.join(  dn ,  'test_ids.npy' )\n        \n    return fn_X, fn_protein_ids, fn_X_submit, fn_submit_protein_ids \n    \n################ load  features #########################\n\nprint(); print('!!! Pay attention proteins ids should be the same as in all the files !!!!!!!!!!!!!!!!!! '); print();\n\ndef get_features(list_features_id, verbose = 0):\n    \n    if verbose >= 100: print(list_features_id ); \n        \n    X_submit,  submit_protein_ids = None, None\n    for i0,features_id in enumerate(list_features_id):\n        # Pay attention proteins ids should be the same as in all the files !!!!!!!!!!!!!!!!!!\n\n        fn_X, fn_protein_ids, fn_X_submit, fn_submit_protein_ids  = get_paths_to_features(features_id)\n        fn = fn_X #  '/kaggle/input/4637427/train_embeds_esm2_t36_3B_UR50D.npy'\n        if verbose >= 100:  print(fn)\n        if i0 == 0:\n            X = np.load(fn).astype(np.float32)[:n_samples_to_consider, :]\n        else:\n            X = np.concatenate( [X, np.load(fn).astype(np.float32)[:n_samples_to_consider, :] ] , axis = 1 )\n        if verbose >= 100: print(X.shape)\n        if verbose >= 100: print(X[:2,:3])\n        protein_ids  = np.load(fn_protein_ids)[:n_samples_to_consider ]\n        vec_train_protein_ids = protein_ids\n        if verbose >= 100: print('protein_ids.shape:', protein_ids.shape)\n        if verbose >= 100: print('protein_ids[:15]:', protein_ids[:15])\n\n            \n        ################ load  features for submit #########################\n        if mode_submit:\n            # fn = '/kaggle/input/4637427/train_embeds_esm2_t36_3B_UR50D.npy'\n            # fn = '/kaggle/input/4637427/test_embeds_esm2_t36_3B_UR50D.npy'\n            fn = fn_X_submit\n            if verbose >= 100: print(fn)\n            # X_submit = np.load(fn).astype(np.float32)\n            if i0 == 0:\n                X_submit = np.load(fn).astype(np.float32)\n            else:\n                X_submit = np.concatenate( [X_submit, np.load(fn).astype(np.float32) ] , axis = 1 )\n            if verbose >= 100: print(X_submit.shape)\n            if verbose >= 100: print(X_submit[:2,:3])\n\n            fn = fn_submit_protein_ids\n            submit_protein_ids = np.load(fn)\n            if verbose >= 100: print(submit_protein_ids.shape, submit_protein_ids[:10])\n    return X,vec_train_protein_ids,X_submit,  submit_protein_ids \n\nX,vec_train_protein_ids,X_submit,  submit_protein_ids = get_features(list_features_id, verbose = 100)\n\n\n\nimport gc\ngc.collect()\nprint('X mbytes:', X.nbytes/1024/1024)\ntry :\n    print('X_submit mbytes:', X_submit.nbytes/1024/1024)\nexcept:\n    pass\n\nlog_available_ram('After features load')","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:09.486562Z","iopub.execute_input":"2023-08-12T19:45:09.487146Z","iopub.status.idle":"2023-08-12T19:45:13.906017Z","shell.execute_reply.started":"2023-08-12T19:45:09.486987Z","shell.execute_reply":"2023-08-12T19:45:13.904767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prepare Y \n\nWill take top frequent labels by subontology. \n\nMain param: dict_labels_to_consider_by_subontology{'BPO':1100, 'MFO':450, 'CCO':300}\n\nWe load  full \"Y\" in sparse matric format where targets already sorted by frequency. We just need to take top frequent according to param above.\n\nIdea to sort by subontology is taken from publinc notebook by Zoltan : https://www.kaggle.com/code/szabo7zoltan/combineembeddings \n","metadata":{}},{"cell_type":"markdown","source":"## Load auxiliary information on labels - subontology markers - BPO, MFO, CCO","metadata":{}},{"cell_type":"code","source":"%%time \ndf_labels_info = pd.read_csv('/kaggle/input/cafa5-features-etc/df_terms_counts_weights_names_etc.csv',index_col = 0)\nprint(df_labels_info.shape )\ndisplay( df_labels_info.head(5) )\nprint( df_labels_info['namespace'].value_counts()   )\n\ndict_labels_ontologies = {}\nfor ont in ['BPO',  'MFO', 'CCO']:\n    m = df_labels_info['namespace'] == ont\n    dict_labels_ontologies[ont] = set( df_labels_info[m].index )\n    print(len(dict_labels_ontologies[ont]), list(dict_labels_ontologies[ont])[:10]   )\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:13.90774Z","iopub.execute_input":"2023-08-12T19:45:13.908099Z","iopub.status.idle":"2023-08-12T19:45:14.315098Z","shell.execute_reply.started":"2023-08-12T19:45:13.90807Z","shell.execute_reply":"2023-08-12T19:45:14.313707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare Y itself ","metadata":{}},{"cell_type":"code","source":"%%time\n\n#    dict_labels_to_consider_by_subontology{'BPO':1100, 'MFO':450, 'CCO':300}\n\nprint(dict_labels_to_consider_by_subontology)\n\n############################ Load Y_labels  ######################################\n\nfn = '/kaggle/input/cafa5-features-etc/Y_31466_labels.npy'\nY_labels = np.load(fn, allow_pickle=True )#[:n_labels_to_consider]\nprint(Y_labels.shape)\nprint(Y_labels[:20])\n\n############################ prepare list_reindex_labels  ######################################\n\nlist_reindex_labels = []\nseries_labels = pd.Series(Y_labels)\nfor ont in ['BPO',  'MFO', 'CCO']:\n    m = series_labels.isin(  dict_labels_ontologies[ont]  )\n    labels_limit = np.min( [dict_labels_to_consider_by_subontology[ont], m.sum() ] )\n    dict_labels_ontologies[ont]  = labels_limit # Just in case parameter is greater than maximum possible \n    list_reindex_labels += list( series_labels[m].index[:labels_limit] )\n    \nprint(len(list_reindex_labels), list_reindex_labels[:10])\n\n############################ take appropriate part of Y_labels  ######################################\n\nY_labels = Y_labels[ np.array(list_reindex_labels ) ]\nlabels_to_consider = Y_labels\nprint(Y_labels.shape)\nprint(Y_labels[:20])\n\n############################ load targets and select appropriate part  ######################################\n\n\n# Load targets Y\nfrom scipy import sparse\nfn = '/kaggle/input/cafa5-features-etc/Y_31466_sparse_float32.npz'\nY = sparse.load_npz(fn )\nprint('Y', Y.shape, 'loaded')\nY = Y[:n_samples_to_consider,np.array(list_reindex_labels )].toarray()\nprint('Y', Y.shape, 'truncated')\nn_labels_to_consider = Y.shape[1] \n\n# %%time\nif 1:\n    v = Y.sum(axis = 0)\n    plt.figure(figsize = (20,6))\n    plt.plot(v, '*-')\n    plt.grid()\n    plt.title(' Number of 1 in targets',fontsize = 20 )\n    plt.xlabel('target index', fontsize = 20 )\n    plt.show()\n\n    \n\n\n# fn_train_terms = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv'\n# fn_train_taxonomy ='/kaggle/input/cafa-5-protein-function-prediction/Train/train_taxonomy.tsv'\n\nimport gc\ngc.collect()\n\nlog_available_ram('After targets load')","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:14.322132Z","iopub.execute_input":"2023-08-12T19:45:14.32254Z","iopub.status.idle":"2023-08-12T19:45:15.791603Z","shell.execute_reply.started":"2023-08-12T19:45:14.322509Z","shell.execute_reply":"2023-08-12T19:45:15.789824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_borders_tmp = pd.Series(dict_labels_ontologies).to_frame().T\ndisplay(df_borders_tmp)\ndf_borders_tmp.to_csv('borders.csv', header = None, index = None)","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:15.793136Z","iopub.execute_input":"2023-08-12T19:45:15.793488Z","iopub.status.idle":"2023-08-12T19:45:15.811503Z","shell.execute_reply.started":"2023-08-12T19:45:15.793458Z","shell.execute_reply":"2023-08-12T19:45:15.810303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('X mbytes:', X.nbytes/1024/1024)\nprint('Y mbytes:', Y.nbytes/1024/1024)\ntry :\n    print('X_submit mbytes:', X_submit.nbytes/1024/1024)\nexcept:\n    pass\nlog_available_ram('After data load')","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:15.813083Z","iopub.execute_input":"2023-08-12T19:45:15.813431Z","iopub.status.idle":"2023-08-12T19:45:15.821382Z","shell.execute_reply.started":"2023-08-12T19:45:15.813402Z","shell.execute_reply":"2023-08-12T19:45:15.820069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare possible downsampling of the train ","metadata":{}},{"cell_type":"code","source":"%%time\ndict_set_allowed_train_indexes = {}\n\nfn = '/kaggle/input/cafa5-features-etc/train_ids_cut43k.npy'\nallowed_train_ids = np.load(fn)\nprint(allowed_train_ids.shape, allowed_train_ids[:10])\nvec_allowed_train_indexes_43k =  [ ix for ix in range(len(vec_train_protein_ids)) if  vec_train_protein_ids[ix] in ( allowed_train_ids ) ] \nset_allowed_train_indexes_43k = set( vec_allowed_train_indexes_43k  )\ndict_set_allowed_train_indexes['43k'] = set_allowed_train_indexes_43k\nprint(len(dict_set_allowed_train_indexes['43k']), list(dict_set_allowed_train_indexes['43k'])[:10] )\n\ndef get_downsampled_IX_train(IX_train, mode_downsample_train ):\n    if mode_downsample_train in dict_set_allowed_train_indexes.keys():\n        set_allowed_train_indexes = dict_set_allowed_train_indexes[ mode_downsample_train ]\n        IX_train = [t for t in IX_train if t in set_allowed_train_indexes ]\n    elif 'random_subsample_percent' in str(mode_downsample_train):\n        random_subsample_percent = float( str(mode_downsample_train).split('_')[-1] ) # 'random_subsample_percent_90'\n        IX_train = np.random.permutation(IX_train)[ :int(len(IX_train)*random_subsample_percent / 100  )]\n    \n    return IX_train\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:15.823203Z","iopub.execute_input":"2023-08-12T19:45:15.823627Z","iopub.status.idle":"2023-08-12T19:45:30.074401Z","shell.execute_reply.started":"2023-08-12T19:45:15.823595Z","shell.execute_reply":"2023-08-12T19:45:30.073478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load partition to folds for cross-validation\n\nfolds_gkf - 5-folds random split","metadata":{}},{"cell_type":"code","source":"%%time \nfn = '/kaggle/input/cafa5-features-etc/random_folds/folds_gkf.npy'\nfolds = np.load(fn)[:n_samples_to_consider]\nprint(folds.shape, len(set(folds)))\nfor k in set(folds):\n    m = folds == k\n    print(k, m.sum() )\nfolds\n\n\n# from sklearn.model_selection import train_test_split\n# IX_train,IX_val = train_test_split( np.arange(len(X)), train_size=0.7, random_state=42)\n# print(IX_train.shape,IX_val.shape) \n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.076142Z","iopub.execute_input":"2023-08-12T19:45:30.076899Z","iopub.status.idle":"2023-08-12T19:45:30.104789Z","shell.execute_reply.started":"2023-08-12T19:45:30.076854Z","shell.execute_reply":"2023-08-12T19:45:30.103477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Define neural networks in Pytorch\n\nExample of simple neural network on pytorch. \n\n\nJust an example randomly taken from previous Kaggle comeptitions, nothing have been optmized for CAFA5 - you are welcome to do it. \n\n\n","metadata":{}},{"cell_type":"code","source":"import torch.nn as nn","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.106264Z","iopub.execute_input":"2023-08-12T19:45:30.107141Z","iopub.status.idle":"2023-08-12T19:45:30.113161Z","shell.execute_reply.started":"2023-08-12T19:45:30.107107Z","shell.execute_reply":"2023-08-12T19:45:30.111355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## https://www.kaggle.com/code/andreylalaley/pytorch-cafa-5-prediction?scriptVersionId=138595845\n## 0.52875 Andrey \nclass Model5(nn.Module):\n    def __init__(self,input_features,output_features):\n        super().__init__()\n        \n        self.activation = nn.PReLU()\n        \n        self.bn1 = nn.BatchNorm1d(input_features)\n        self.fc1 = nn.Linear(input_features, 800)\n        self.ln1 = nn.LayerNorm(800, elementwise_affine=True)\n        \n        self.bn2 = nn.BatchNorm1d(800)\n        self.fc2 = nn.Linear(800, 600)\n        self.ln2 = nn.LayerNorm(600, elementwise_affine=True)\n        \n        self.bn3 = nn.BatchNorm1d(600)\n        self.fc3 = nn.Linear(600, 400)\n        self.ln3 = nn.LayerNorm(400, elementwise_affine=True)\n        \n        self.bn4 = nn.BatchNorm1d(1200)\n        self.fc4 = nn.Linear(1200, output_features)\n        self.ln4 = nn.LayerNorm(output_features, elementwise_affine=True)\n        \n        self.sigm = nn.Sigmoid()\n    def forward(self,inputs):\n#         print(inputs.shape)\n\n        fc1_out = self.bn1(inputs)\n        fc1_out = self.ln1(self.fc1(inputs))\n        fc1_out = self.activation(fc1_out)\n        \n        x = self.bn2(fc1_out)\n        \n        x = self.ln2(self.fc2(x))\n        x = self.activation(x)\n        \n        x = self.bn3(x)\n        \n        x = self.ln3(self.fc3(x))\n        x = self.activation(x)\n        \n        x = torch.cat([x, fc1_out], axis = -1)\n        \n        x = self.bn4(x)\n        \n        x = self.ln4(self.fc4(x))\n        out = self.sigm(x)\n        return out\n    \n    \n## https://www.kaggle.com/code/mrgobus/pytorch-01-basics-49cdbc?scriptVersionId=138589166&cellId=27\n## 0.49651 Ivan\nclass Model4(nn.Module):\n\n    def __init__(self, in_features, out_features, neurons_per_layer = 1000):\n\n        super().__init__()\n\n        self.in_features = in_features\n        self.out_features = out_features\n\n        self.model = nn.Sequential(\n            nn.BatchNorm1d(in_features),\n            nn.Dropout(0.2),\n            nn.Linear(in_features, neurons_per_layer),\n            nn.LayerNorm(neurons_per_layer),\n            nn.PReLU(init = 0.5),\n\n            nn.BatchNorm1d(neurons_per_layer),\n            nn.Dropout(0.2),\n            nn.Linear(neurons_per_layer, neurons_per_layer),\n            nn.LayerNorm(neurons_per_layer),\n            nn.PReLU(init = 0.5),\n\n            nn.BatchNorm1d(neurons_per_layer),\n            nn.Linear(neurons_per_layer, out_features),\n            nn.LayerNorm(out_features),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x):\n        return self.model(x)\n    \n    \n## https://www.kaggle.com/code/sebastian157/pytorch-01-basics-3a6104?scriptVersionId=138569077&cellId=18\n## LB 0.49545  Boris\nclass Model3(nn.Module):\n    def __init__(self, in_features, out_features):\n        super().__init__()\n        self.in_features = in_features\n        self.out_features = out_features\n\n        self.model = nn.Sequential(\n            nn.BatchNorm1d(in_features),\n         \n            nn.Linear(in_features, 800),\n            nn.LayerNorm(800, elementwise_affine=True),\n            nn.PReLU(),\n            \n            nn.BatchNorm1d(800),  \n            \n            nn.Linear(800, 600),\n            nn.LayerNorm(600, elementwise_affine=True),\n            nn.PReLU(),\n            \n            nn.BatchNorm1d(600),\n          \n            nn.Linear(600, 400),\n            nn.LayerNorm(400, elementwise_affine=True),\n            nn.PReLU(),\n            \n            nn.BatchNorm1d(400),\n       \n            nn.Linear(400, out_features),\n            nn.LayerNorm(out_features, elementwise_affine=True),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x):\n        return self.model(x)\n\n\n## https://www.kaggle.com/code/lidiashishina/pytorch-01-basics?scriptVersionId=138484544&cellId=18\n# LB 0.4908 Lidia\nclass Model2(nn.Module):\n    def __init__(self, in_features, out_features):\n        super().__init__()\n        self.in_features = in_features\n        self.out_features = out_features\n\n        self.model = nn.Sequential(\n            nn.BatchNorm1d(in_features),\n         \n            nn.Linear(in_features, 800),\n            nn.LayerNorm(800, elementwise_affine=True),\n            nn.GELU(),\n      \n            nn.BatchNorm1d(800),  \n            \n            nn.Linear(800, 600),\n            nn.LayerNorm(600, elementwise_affine=True),\n            nn.GELU(),\n\n            nn.BatchNorm1d(600),\n          \n            nn.Linear(600, 400),\n            nn.LayerNorm(400, elementwise_affine=True),\n            nn.GELU(),\n\n            nn.BatchNorm1d(400),\n       \n            nn.Linear(400, out_features),\n            nn.LayerNorm(out_features, elementwise_affine=True),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x):\n        return self.model(x)\n    \n# https://www.kaggle.com/code/alexandervc/moa-nn-04/notebook\n\nclass Model1(nn.Module): # The parent class for the models is nn.Module \n    \n    def __init__(self, in_features, out_features): # constructor \n        \n        super().__init__() # the constructor of the upper class is first called\n\n        self.in_features = in_features\n        self.out_features = out_features\n\n        self.model = nn.Sequential( #  Sequential addition of layers -  multi-layer perceptron \n            nn.BatchNorm1d(in_features),\n            nn.Linear(in_features, 1000),\n            nn.GELU(), # nn.ReLU(),\n\n            nn.BatchNorm1d(1000),            # \n            nn.Dropout(0.35),\n            nn.Linear(1000, 2000),\n            nn.GELU(),# nn.ReLU(),\n\n#             nn.BatchNorm1d(600),            # nn.Dropout(0.1),\n#             nn.Linear(600, 400),\n#             nn.ReLU(),\n\n#             nn.BatchNorm1d(400),\n            nn.Linear(2000, out_features),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x): # \n        return self.model(x)","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.115457Z","iopub.execute_input":"2023-08-12T19:45:30.116247Z","iopub.status.idle":"2023-08-12T19:45:30.151468Z","shell.execute_reply.started":"2023-08-12T19:45:30.116186Z","shell.execute_reply":"2023-08-12T19:45:30.150114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ModelPytorchSimple1(nn.Module):\n    '''\n    MLP with two hidden layers \n    '''\n    def __init__(self, in_features, out_features, neurons_per_hidden_layer1 = 1000,neurons_per_hidden_layer2 = 1000):\n\n        super().__init__()\n\n        self.in_features = in_features\n        self.out_features = out_features\n\n        self.model = nn.Sequential(\n            #nn.BatchNorm1d(in_features),\n            #nn.Dropout(0.2),\n            nn.Linear(in_features, neurons_per_hidden_layer1),\n            #nn.LayerNorm(neurons_per_layer),\n            #nn.PReLU(init = 0.5),\n            nn.ReLU(),\n            \n            #nn.BatchNorm1d(neurons_per_layer),\n            #nn.Dropout(0.2),\n            nn.Linear(neurons_per_hidden_layer1, neurons_per_hidden_layer2),\n            #nn.LayerNorm(neurons_per_layer),\n            #nn.PReLU(init = 0.5),\n            nn.ReLU(),\n            \n            #nn.BatchNorm1d(neurons_per_layer),\n            nn.Linear(neurons_per_hidden_layer2, out_features),\n            #nn.LayerNorm(out_features),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x):\n        return self.model(x)    \n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.153126Z","iopub.execute_input":"2023-08-12T19:45:30.153607Z","iopub.status.idle":"2023-08-12T19:45:30.170507Z","shell.execute_reply.started":"2023-08-12T19:45:30.153563Z","shell.execute_reply":"2023-08-12T19:45:30.169393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"    \nmodel_config_test = {'id':'pMLP_Configurable1',  'Layers':[100,200] }\n\nclass ModelPytorchSequentialConfigurable1(nn.Module):\n    '''\n    MLP with congfigurable params - layers, dropouts, BatchNormalizations, LayerNorm\n    '''\n    def __init__(self, in_features, out_features, model_config):\n\n        super().__init__()\n\n        self.in_features = in_features\n        self.out_features = out_features\n        self.str_model_id = 'update at end'\n        str_model_id = 'pMLP_Configurable1'\n\n        layers_sizes = model_config.get( 'Layers', [500,1000] )\n        list_droupouts = model_config.get('Dropouts' ,  []  )\n        list_batchnormalization = model_config.get( 'BatchNormalizations', [] )\n        list_LayerNorm = model_config.get( 'LayerNorms', [] )\n        \n        self.model = nn.Sequential()\n        \n        for current_layer in range(0, len(layers_sizes) + 1):\n            if current_layer == 0:\n                size_in = in_features; size_out = layers_sizes[current_layer]\n                str_model_id += '_L'+str(size_out)\n            elif current_layer ==  len(layers_sizes):\n                size_in = layers_sizes[current_layer-1]; size_out = out_features\n            else:\n                size_in = layers_sizes[current_layer-1]; size_out = layers_sizes[current_layer]\n                str_model_id += '_L'+str(size_out)\n\n            if (len(list_batchnormalization) > current_layer) and list_batchnormalization[current_layer]:\n                self.model.add_module('BatchNorm1d ' + str( current_layer ), nn.BatchNorm1d(size_in) )\n                str_model_id += '_BN'\n            if (len(list_droupouts) > current_layer) and (list_droupouts[current_layer] is not None):\n                self.model.add_module('Dropout ' + str( current_layer ), nn.Dropout(list_droupouts[current_layer]) )\n                str_model_id += '_DR'+str(list_droupouts[current_layer])\n            self.model.add_module('Dense ' + str( current_layer ), nn.Linear(size_in, size_out) )\n            if (len(list_LayerNorm) > current_layer) and (list_LayerNorm[current_layer] is not None):\n                self.model.add_module('LayerNorm ' + str( current_layer ), nn.LayerNorm(size_out) )\n                str_model_id += '_LN'\n            if current_layer < len(layers_sizes):\n                self.model.add_module('ReLU ' + str( current_layer ), nn.ReLU() )\n            else: \n                self.model.add_module('Sigmoid ' + str( current_layer ),  nn.Sigmoid()  )\n\n        self.epochs = model_config.get('epochs', 10)\n        if self.epochs != 10: str_model_id += '_E'+str( self.epochs )\n        self.batch_size = model_config.get('batch_size', 128)\n        if self.batch_size != 128: str_model_id += '_BS'+str( self.batch_size )\n        self.LR = model_config.get('LR', 0.001)\n        if self.LR != 0.001: str_model_id += '_LR'+str( self.LR )\n                \n        self.str_model_id = str_model_id\n        \n    def forward(self, x):\n        return self.model(x)    \nmodel_config_test = {'id':'pMLP_Configurable1',  'Layers':[100,200], 'BatchNormalizations': [True, False ],'Dropouts':[0.5, 0.2]   }\nmodel = ModelPytorchSequentialConfigurable1(1280,1850, model_config_test)\nprint(model.str_model_id)\nmodel\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T20:01:14.879096Z","iopub.execute_input":"2023-08-12T20:01:14.879903Z","iopub.status.idle":"2023-08-12T20:01:14.91811Z","shell.execute_reply.started":"2023-08-12T20:01:14.879837Z","shell.execute_reply":"2023-08-12T20:01:14.916948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport gc\nimport time \nimport numpy as np\ntry: \n    from sklearn.metrics import roc_auc_score\n    from sklearn.metrics import f1_score\n    from sklearn.linear_model import Ridge\n    from sklearn.neural_network import MLPRegressor\n    import lightgbm as lgbm\n    from sklearn.multioutput import MultiOutputRegressor\n    from sklearn.multioutput import MultiOutputClassifier  \n    \n    from catboost import CatBoostRegressor\n    from catboost import CatBoostClassifier    \nexcept Exception as e:\n    print(f'Exception importing models {e} ')\n    pass\n\n\nimport keras\nfrom keras.models import Sequential, Model\nfrom keras.layers import Dense, Input,  Concatenate, Dropout, BatchNormalization, Activation\nfrom keras.layers import LayerNormalization\n\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.215343Z","iopub.execute_input":"2023-08-12T19:45:30.216006Z","iopub.status.idle":"2023-08-12T19:45:30.226694Z","shell.execute_reply.started":"2023-08-12T19:45:30.215973Z","shell.execute_reply":"2023-08-12T19:45:30.225373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Get model function: Pytorch, Keras, Sklearn, etc models ","metadata":{}},{"cell_type":"code","source":"def get_model(model_config):\n    \n    str_model_id = str(model_config['id'])\n    if model_config['id'] == 'pMLP_Andrey':\n        model = Model5(X.shape[1],Y.shape[1])\n        model.to(device)\n    elif model_config['id'] == 'Ridge':\n        alpha = model_config.get('alpha',10)\n        str_model_id = 'Ridge'+str(alpha)\n        model = Ridge(alpha=alpha)\n    elif model_config['id'] == 'pMLP_simple1':\n        \n        layers_sizes = model_config.get(  'Layers',[500,1000] )\n        model = ModelPytorchSimple1(X.shape[1],Y.shape[1],neurons_per_hidden_layer1 = layers_sizes[0],neurons_per_hidden_layer2 = layers_sizes[1])\n        model.to(device)\n        str_model_id += '_L1_'+str(layers_sizes[0])+'_L2_'+str(layers_sizes[1])\n        \n    if model_config['id'] == 'pMLP_Configurable1':\n        model = ModelPytorchSequentialConfigurable1(X.shape[1],Y.shape[1], model_config)\n        str_model_id = model.str_model_id\n        model.to(device)\n        \n    elif model_config['id'] == 'skMLP':\n        # model_config = {'id':'skMLP', 'Layers': [ 500,800] ,'LR':  0.001 , 'alpha':1e-4, 'max_iter':500 }\n        str_model_id = 'skMLP'\n        max_iter = model_config.get('max_iter',500); str_model_id += '_MI'+str(max_iter) \n        random_state = model_config.get('random_state',np.random.randint(0,100)); str_model_id += '_RS'+str(random_state)\n        hidden_layer_sizes = model_config.get('Layers', (100,) ); str_model_id += '_HL'+str(hidden_layer_sizes)\n        alpha = model_config.get('alpha', 1e-4 ); str_model_id += '_alpha'+str(alpha)\n        learning_rate_init = model_config.get('LR', 0.001 ); str_model_id += '_LR'+str(learning_rate_init)\n        model = MLPRegressor( max_iter=max_iter , random_state=random_state, hidden_layer_sizes = hidden_layer_sizes, alpha= alpha, learning_rate_init = learning_rate_init  )\n        \n    elif model_config['id'] == 'KMLP_simple':\n        model = Sequential(); str_model_id = 'KMLP_simple'\n        \n        nfeats = X.shape[1]#  model_config['input_dim']\n        nlabels = Y.shape[1]#  model_config['output_dim']\n        \n        layer_dim = 500\n        model.add(Dense(layer_dim, activation='relu', input_dim=nfeats   ) )  ; str_model_id += '_L1_'+str(layer_dim) \n        model.add(BatchNormalization()) ; str_model_id += '_BN'\n        droupout_rate = 0.1\n        model.add(Dropout(droupout_rate)) ; str_model_id += '_DR'+str( np.round(droupout_rate ,2) ) \n        \n        layer_dim = 800\n        model.add(Dense(layer_dim, activation='relu'   ) )  ; str_model_id += '_L2_'+str(layer_dim) \n        model.add(BatchNormalization()) ; str_model_id += '_BN'\n        droupout_rate = 0.1\n        model.add(Dropout(droupout_rate)) ; str_model_id += '_DR'+str( np.round(droupout_rate ,2) ) \n        \n        ######### Last layer ########################################################################################\n        model.add(Dense(nlabels, activation='sigmoid'   ))\n        \n        model.compile(loss='binary_crossentropy', optimizer='adam',  metrics=[keras.metrics.AUC()] )\n        \n        \n    elif model_config['id'] == 'KMLP':   \n        # {'id':'KMLP', 'Layers': [500,800],'Dropouts': [0.1 ],  'BatchNormalizations': [] }\n        model = Sequential(); str_model_id = 'KMLP'\n        \n        layers_sizes = model_config['Layers']\n        list_droupouts = model_config.get('Dropouts' ,  []  )\n        list_batchnormalization = model_config.get( 'BatchNormalizations', [] )\n        \n        nfeats = X.shape[1]#  model_config['input_dim']\n        nlabels = Y.shape[1]#  model_config['output_dim']\n        \n        ######### First Layer #################################################################################\n        i_layer = 0 \n        layer_dim = layers_sizes[i_layer]\n        model.add(Dense(layer_dim, activation='relu', input_dim=nfeats   ) )  ; str_model_id += '_L1_'+str(layer_dim) \n        #model.add(Dense(layer_dim, activation='relu', input_dim=nfeats, kernel_regularizer = keras.regularizers.l2( 0  )  ) )  ; str_model_id += '_L1_'+str(layer_dim) \n        if len( list_batchnormalization ) > i_layer:\n            if list_batchnormalization[i_layer]:\n                model.add(BatchNormalization())\n                str_model_id += '_BN'\n        if (len( list_droupouts ) > i_layer) and ( list_droupouts[i_layer] is not None ):\n            droupout_rate  = list_droupouts[i_layer]\n            model.add(Dropout(droupout_rate))\n            str_model_id += '_DR'+str( np.round(droupout_rate ,2) ) \n            \n        ######### Middle layers ##########################################################################################\n        for i_layer in range(1, len( layers_sizes ) ) : \n            layer_dim = layers_sizes[i_layer]\n            if layer_dim is None: break  \n            # model.add(Dense(layer_dim, activation='relu' ,  kernel_regularizer = keras.regularizers.l2( 0 ) ));                 str_model_id += '_L'+str(i_layer+1)+'_'+str(layer_dim)    \n            model.add(Dense(layer_dim, activation='relu' ));                 str_model_id += '_L'+str(i_layer+1)+'_'+str(layer_dim)    \n#             model.add(LayerNormalization())\n            if len( list_batchnormalization ) > i_layer:\n                if list_batchnormalization[i_layer]:\n                    model.add(BatchNormalization())\n                    str_model_id += '_BN'\n            if (len( list_droupouts ) > i_layer) and ( list_droupouts[i_layer] is not None ):\n                droupout_rate  = list_droupouts[i_layer]\n                model.add(Dropout(droupout_rate))\n                str_model_id += '_DR'+str( np.round(droupout_rate ,2) )         \n        \n        ######### Last layer ########################################################################################\n        model.add(Dense(nlabels, activation='sigmoid'   ))\n#         model.add(Dense(nlabels, activation='sigmoid' ,  kernel_regularizer = keras.regularizers.l2( 0  ) ))\n#         model.add(Dense(nlabels, activation='tanh'))\n        #model.add(Dense(nlabels, activation='relu'))\n        \n        \n        model.compile(loss='binary_crossentropy',\n                        optimizer='adam',\n                        metrics=[keras.metrics.AUC()])           \n        \n    elif model_config['id'] == 'gpuLogReg':\n        import torch\n        if torch.cuda.is_available():\n            from cuml.linear_model import LogisticRegression as CULogisticRegression        \n            model = MultiOutputClassifier(estimator= CULogisticRegression( ) )# **params ) )\n            str_model_id = 'gpuLogReg'\n            \n    elif model_config['id'] == 'gpuCatBClasdefault':\n        model = MultiOutputClassifier(estimator= CatBoostClassifier(task_type = 'GPU', verbose = 0 ))# **params ) )\n        str_model_id = model_cfg[0] \n\n    elif model_config['id'] == 'gpuCatBdefault': \n        model = MultiOutputRegressor(estimator= CatBoostRegressor(task_type = 'GPU', verbose = 0 ))# **params ) )\n        str_model_id = model_cfg[0] \n\n    elif model_config['id'] == 'CatBdefault': \n        model = MultiOutputRegressor(estimator= CatBoostRegressor(verbose = 0 ))# **params ) )\n        str_model_id = model_cfg[0] \n    \n    elif model_config['id'] == 'LGBdefault': \n        model = MultiOutputRegressor(estimator=lgbm.LGBMRegressor())# **params ) )\n        str_model_id = 'LGBdefault'\n            \n    namepostfix = model_config.get('namepostfix',\"\")\n    str_model_id += namepostfix\n    return model, str_model_id\n\nmodel_config_tmp = {'id':'Ridge'}\nmodel_config_tmp = {'id':'Ridge', 'alpha':10,'namepostfix':'_test' }\n# model_config_tmp = {'id':'skMLP', 'alpha':1e-4, 'Layers': [500,1000] }\n# model_config_tmp = {'id':'pMLP_Andrey' , 'n_selfblend':1 }  # Pytorch model\n# model_config_tmp = {'id':'skMLP', 'Layers': [ 500,800] ,'LR':  0.001 , 'alpha':1e-4, 'max_iter':500 }  # Sklearn MLP model\n# model_config_tmp = {'id':'KMLP_simple' , 'epochs':15 ,   'batch_size':128 , 'verbose':0}  # Simple Keras MLP\n# model_config_tmp = {'id':'KMLP', 'Layers': [400,600],'Dropouts': [0.1 ],  'BatchNormalizations': [] , 'epochs':15 ,   'batch_size':128 , 'verbose':0 }  \n\nmodel, str_model_id = get_model(model_config_tmp)\nprint(model, str_model_id )","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.22908Z","iopub.execute_input":"2023-08-12T19:45:30.229903Z","iopub.status.idle":"2023-08-12T19:45:30.277488Z","shell.execute_reply.started":"2023-08-12T19:45:30.229866Z","shell.execute_reply":"2023-08-12T19:45:30.275923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fit model function ","metadata":{}},{"cell_type":"code","source":"%%time\n\n# Call function by name: \n# https://stackoverflow.com/questions/3061/calling-a-function-of-a-module-by-using-its-name-a-string\n# def tt2(a,b):\n#     print(a+b)\n#     return a+b\n# locals()[\"tt\"](1,2)\n\ndef model_fit(model , X_train, Y_train, model_config , str_model_id = '', verbose = 1000 ):\n    \n    if verbose >= 100:\n        print('model_fit', str_model_id,  model_config )\n\n    if 'custom_model_fit_function_name' in model_config.keys():\n        str_func_name = model_config['custom_model_fit_function_name']\n        if verbose >= 1000: print('str_func_name', str_func_name)\n        #locals()[str_func_name](model , X_train, Y_train, model_config , str_model_id = '', verbose = 1000 ) # Call a function with name \"str_func_name\"\n        globals()[str_func_name](model , X_train, Y_train, model_config , str_model_id = '', verbose = 1000 ) # Call a function with name \"str_func_name\"\n        # https://stackoverflow.com/questions/3061/calling-a-function-of-a-module-by-using-its-name-a-string\n\n    elif  'pMLP' in str(model_config['id']):\n        model_fit_pytorch( model , X_train, Y_train, model_config , str_model_id, verbose )\n        \n    else:\n        \n        dict_prm_for_fit = {t : model_config[t] for t in ['epochs' ,   'batch_size' , 'verbose']   if t in model_config.keys()  }\n        if len( dict_prm_for_fit ) == 0:\n            # For models like sklearn we just write : \n            model.fit( X_train, Y_train  )\n        else:\n            # For Keras models we can specify epochs, batch_size, etc \n            model.fit( X_train, Y_train , **dict_prm_for_fit  )\n            \ndef model_fit_pytorch( model , X_train, Y_train, model_config , str_model_id = '', verbose = 1000 ):\n    \n    criterion = model_config.get( 'criterion', nn.BCELoss() )\n    max_epoch = model_config.get('epochs', getattr(model,'epochs',10) )\n    BATCH_SIZE = model_config.get('batch_size',  getattr(model,'batch_size' ,128) )\n    LEARNING_RATE = model_config.get( 'LR' , getattr(model, 'LR' , 0.001) )\n    \n    optimizer = model_config.get( 'optimizer' ,  torch.optim.Adam(model.parameters(), lr=LEARNING_RATE)  )\n#         optimizer = torch.optim.Adam(model.parameters(), lr=LEARNING_RATE) \n#         optimizer = SophiaG(model.parameters(),lr=LEARNING_RATE, betas=(0.965, 0.99), rho = 0.01, weight_decay=1e-1)\n#         optimizer = torch.optim.RMSprop(model.parameters(), lr=LEARNING_RATE) \n#         optimizer = torch.optim.SGD(model.parameters(), lr=LEARNING_RATE)\n    \n    lr_sched = model_config.get( 'lr_scheduler' , None ) \n#     lr_sched = torch.optim.lr_scheduler.CosineAnnealingWarmRestarts(optimizer, T_0=10, T_mult=2, eta_min=0.001, last_epoch=-1) # https://www.kaggle.com/code/lidiashishina/pytorch-01-basics?scriptVersionId=138484544&cellId=33\n#     lr_sched = torch.optim.lr_scheduler.StepLR(optimizer, 0.5, 5)    \n\n    X_train = torch.tensor(X_train, dtype=torch.float32).to(device)\n    Y_train = torch.tensor(Y_train, dtype=torch.float32).to(device)\n    train_dataset = TensorDataset(X_train, Y_train)\n    train_dataloader = DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True) \n\n    if verbose >= 100:\n        print(str_model_id, 'Start model training. X_train.shape, Y_train.shape', X_train.shape, Y_train.shape )\n    log_available_ram(f'After X_train initialization. Right before model train {str_model_id}')        \n    \n    ##################### Model Training ###################################################\n    for epoch in range(max_epoch):\n        t0_epoch = time.time()\n        model.train() # switch model into train mode. This helps inform layers such as Dropout and BatchNorm, which are designed to behave differently during training and evaluation. For instance, in training mode, BatchNorm updates a moving average on each new batch; whereas, for evaluation mode, these updates are frozen.\n        for i_batch, (x_batch, y_batch) in enumerate(train_dataloader): # Loop ove batches\n            # x_batch, y_batch = x_batch.to(device), y_batch.to(device) # do we need it ? may be already on device \n            optimizer.zero_grad() # technical - set gradients to zero, otherwise they will be accumulated         \n\n            preds = model(x_batch)# Compute predictions only for batch samples \n\n            loss = criterion(preds, y_batch) # Compute loss function for batch predictions\n\n            loss.backward() # Compute gradients\n            optimizer.step() # Update NN weights using gradients\n\n        if lr_sched is not None:\n            lr_sched.step() # Step LR scheduler\n\n        if (verbose >= 10 ) and (i_batch % 100 == 0)  :\n            print(str_model_id,  f'Epoch: {epoch}, batch: {i_batch},  train loss on batch: {loss.item():12.5f} , time: {time.time() - t0:.1f} ' )         \n            \n            \ndef model_fit_pytorch_Sophia_Andrey1(  model , X_train, Y_train, model_config , str_model_id = '', verbose = 1000 ):\n    '''\n    Fit pytorch model with Sophia optimizer, with certain params borrowed from: https://www.kaggle.com/code/andreylalaley/pytorch-cafa-5-prediction?scriptVersionId=138595845&cellId=32\n    '''\n    criterion = model_config.get( 'criterion', nn.BCELoss() )\n    max_epoch = model_config.get('epochs', 39 ) #  EPOCHS = 39\n    BATCH_SIZE = model_config.get('batch_size',  128 )\n    LEARNING_RATE = model_config.get( 'LR' , 0.001 )\n    \n    criterion = model_config.get( 'criterion', nn.BCELoss() )\n\n    lr = LEARNING_RATE # learning rate \n    optimizer = SophiaG(model.parameters(),lr , betas=(0.965, 0.99), rho = 0.01, weight_decay=1e-1)\n\n    \n    X_train = torch.tensor(X_train, dtype=torch.float32).to(device)\n    Y_train = torch.tensor(Y_train, dtype=torch.float32).to(device)\n    train_dataset = TensorDataset(X_train, Y_train)\n    train_dataloader = DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True) \n\n    if verbose >= 100:\n        print(str_model_id, 'Start model training. X_train.shape, Y_train.shape', X_train.shape, Y_train.shape )\n    log_available_ram(f'After X_train initialization. Right before model train {str_model_id}')        \n        \n    ##################### Model Training ###################################################\n    for epoch in range(max_epoch):\n        t0_epoch = time.time()\n        model.train() # switch model into train mode. This helps inform layers such as Dropout and BatchNorm, which are designed to behave differently during training and evaluation. For instance, in training mode, BatchNorm updates a moving average on each new batch; whereas, for evaluation mode, these updates are frozen.\n        for i_batch, (x_batch, y_batch) in enumerate(train_dataloader): # Loop ove batches\n            # x_batch, y_batch = x_batch.to(device), y_batch.to(device) # do we need it ? may be already on device \n            optimizer.zero_grad() # technical - set gradients to zero, otherwise they will be accumulated         \n\n            preds = model(x_batch)# Compute predictions only for batch samples \n\n            loss = criterion(preds, y_batch) # Compute loss function for batch predictions\n\n            loss.backward() # Compute gradients\n            optimizer.step() # Update NN weights using gradients\n\n        # lr_sched is not used in original code     \n        #if lr_sched is not None:\n        #    lr_sched.step() # Step LR scheduler\n\n        if (verbose >= 10 ) and (i_batch % 100 == 0)  :\n            print(str_model_id,  f'Epoch: {epoch}, batch: {i_batch},  train loss on batch: {loss.item():12.5f} , time: {time.time() - t0:.1f} ' )         \n    ","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.279821Z","iopub.execute_input":"2023-08-12T19:45:30.280191Z","iopub.status.idle":"2023-08-12T19:45:30.317823Z","shell.execute_reply.started":"2023-08-12T19:45:30.280159Z","shell.execute_reply":"2023-08-12T19:45:30.316434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Models predict function","metadata":{}},{"cell_type":"code","source":"def model_predict(model , XX,  model_config , str_model_id = '', verbose = 1000 ):\n    \n    if verbose >= 100:\n        print('model_predict',  str_model_id,  model_config )\n\n    if  'pMLP' in str(model_config['id']):\n        Y_pred = model_predict_pytorch( model ,XX,  model_config , str_model_id, verbose )\n    else:\n        Y_pred = model.predict( XX ) \n    \n    return Y_pred\n\n\ndef model_predict_pytorch( model ,XX,  model_config , str_model_id = '', verbose = 1000 ):\n    \n    t0_submit = time.time()\n    XX = torch.tensor(XX, dtype=torch.float32).to(device)\n    \n    model.eval()\n    with torch.no_grad():\n        Y_pred = model(XX).cpu().numpy()\n        \n    if verbose >= 100: print(str_model_id,  f'Y_pred.shape: {Y_pred.shape}, type(Y_pred): {type(Y_pred)}, predict on submit time: {time.time() - t0_submit :.1f} ' )  \n\n    return Y_pred","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.319513Z","iopub.execute_input":"2023-08-12T19:45:30.320002Z","iopub.status.idle":"2023-08-12T19:45:30.334882Z","shell.execute_reply.started":"2023-08-12T19:45:30.31994Z","shell.execute_reply":"2023-08-12T19:45:30.333777Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CAFA5 metric (speed-upped version by sparse matrices)\n\nEvaluation for CAFA \nhttps://github.com/BioComputingUP/CAFA-evaluator\n\nVersion updated to work on Kaggle and speed-upped by sparse matric usage \n\nPay attention - that metric computation is very slow and RAM consuming - be careful ! \n\nIt is quite technical - no need to go into details - just use as a blacbox.\n\nhttps://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/420241 - Anton Vakhrushev - correcting error in the initial code of the metric computation - please upvote \n","metadata":{}},{"cell_type":"code","source":"%%time\n\n# Evaluation for CAFA \n# https://github.com/BioComputingUP/CAFA-evaluator\n\nflag_correct_metric_computation_bug_found_by_Anton = True\n# https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/420241 - Anton Vakhrushev - correcting error in the initial code of the metric computation - please upvote \n\n\nimport numpy as np\nimport pandas as pd\nimport multiprocessing as mp\nimport copy\nimport logging\n\nclass Graph:\n    \"\"\"\n    Ontology class. One ontology == one namespace\n    DAG is the adjacence matrix (sparse) which represent a Directed Acyclic Graph where\n    DAG(i,j) == 1 means that the go term i is_a (or is part_of) j\n    Parents that are in a different namespace are discarded\n    \"\"\"\n    def __init__(self, namespace, terms_dict, ia_dict=None, orphans=False):\n        \"\"\"\n        terms_dict = {term: {name: , namespace: , def: , alt_id: , rel:}}\n        \"\"\"\n        self.namespace = namespace\n        self.dag = []  # [[], ...] terms (rows, axis 0) x parents (columns, axis 1)\n        self.terms_dict = {}  # {term: {index: , name: , namespace: , def: }  used to assign term indexes in the gt\n        self.terms_list = []  # [{id: term, name:, namespace: , def:, adg: [], children: []}, ...]\n        self.idxs = None  # Number of terms\n        self.order = None\n        self.toi = None\n        self.ia = None\n\n        rel_list = []\n        for self.idxs, (term_id, term) in enumerate(terms_dict.items()):\n            rel_list.extend([[term_id, rel, term['namespace']] for rel in term['rel']])\n            self.terms_list.append({'id': term_id, 'name': term['name'], 'namespace': namespace, 'def': term['def'],\n                                 'adj': [], 'children': []})\n            self.terms_dict[term_id] = {'index': self.idxs, 'name': term['name'], 'namespace': namespace, 'def': term['def']}\n            for a_id in term['alt_id']:\n                self.terms_dict[a_id] = copy.copy(self.terms_dict[term_id])\n        self.idxs += 1\n\n        self.dag = np.zeros((self.idxs, self.idxs), dtype='bool')\n\n        # id1 term (row, axis 0), id2 parent (column, axis 1)\n        for id1, id2, ns in rel_list:\n            if self.terms_dict.get(id2):\n                i = self.terms_dict[id1]['index']\n                j = self.terms_dict[id2]['index']\n                self.dag[i, j] = 1\n                self.terms_list[i]['adj'].append(j)\n                self.terms_list[j]['children'].append(i)\n                logging.debug(\"i,j {},{} {},{}\".format(i, j, id1, id2))\n            else:\n                logging.debug('Skipping branch to external namespace: {}'.format(id2))\n        logging.debug(\"dag {}\".format(self.dag))\n        # Topological sorting\n        self.top_sort()\n        logging.debug(\"order sorted {}\".format(self.order))\n\n        if orphans:\n            self.toi = np.arange(self.dag.shape[0])  # All terms, also those without parents\n        else:\n            self.toi = np.nonzero(self.dag.sum(axis=1) > 0)[0]  # Only terms with parents\n        logging.debug(\"toi {}\".format(self.toi))\n\n        if ia_dict is not None:\n            self.set_ia(ia_dict)\n\n        return\n\n    def top_sort(self):\n        \"\"\"\n        Takes a sparse matrix representing a DAG and returns an array with nodes indexes in topological order\n        https://en.wikipedia.org/wiki/Topological_sorting\n        \"\"\"\n        indexes = []\n        visited = 0\n        (rows, cols) = self.dag.shape\n\n        # create a vector containing the in-degree of each node\n        in_degree = self.dag.sum(axis=0)\n        # logging.debug(\"degree {}\".format(in_degree))\n\n        # find the nodes with in-degree 0 (leaves) and add them to the queue\n        queue = np.nonzero(in_degree == 0)[0].tolist()\n        # logging.debug(\"queue {}\".format(queue))\n\n        # for each element of the queue increment visits, add them to the list of ordered nodes\n        # and decrease the in-degree of the neighbor nodes\n        # and add them to the queue if they reach in-degree == 0\n        while queue:\n            visited += 1\n            idx = queue.pop(0)\n            indexes.append(idx)\n            in_degree[idx] -= 1\n            l = self.terms_list[idx]['adj']\n            if len(l) > 0:\n                for j in l:\n                    in_degree[j] -= 1\n                    if in_degree[j] == 0:\n                        queue.append(j)\n\n        # if visited is equal to the number of nodes in the graph then the sorting is complete\n        # otherwise the graph can't be sorted with topological order\n        if visited == rows:\n            self.order = indexes\n        else:\n            raise Exception(\"The sparse matrix doesn't represent an acyclic graph\")\n\n    def set_ia(self, ia_dict):\n        self.ia = np.zeros(self.idxs, dtype='float')\n        for term_id in self.terms_dict:\n            if ia_dict.get(term_id):\n                self.ia[self.terms_dict[term_id]['index']] = ia_dict.get(term_id)\n            else:\n                logging.debug('Missing IA for term: {}'.format(term_id))\n        # Convert inf to zero\n        np.nan_to_num(self.ia, copy=False, nan=0, posinf=0, neginf=0)\n        self.toi = np.nonzero(self.ia > 0)[0]\n\n\nclass Prediction:\n    \"\"\"\n    The score matrix contains the scores given by the predictor for every node of the ontology\n    \"\"\"\n    def __init__(self, ids, matrix, idx, namespace=None):\n        self.ids = ids\n        self.matrix = matrix  # scores\n        self.next_idx = idx\n        # self.n_pred_seq = idx + 1\n        self.namespace = namespace\n\n    def __str__(self):\n        return \"\\n\".join([\"{}\\t{}\\t{}\".format(index, self.matrix[index], self.namespace) for index, _id in enumerate(self.ids)])\n\n\nclass GroundTruth:\n    def __init__(self, ids, matrix, namespace=None):\n        self.ids = ids\n        self.matrix = matrix\n        self.namespace = namespace\n\n\ndef propagate(matrix, ont, order, mode='max'):\n    \"\"\"\n    Update inplace the score matrix (proteins x terms) up to the root taking the max between children and parents\n    \"\"\"\n    if matrix.shape[0] == 0:\n        raise Exception(\"Empty matrix\")\n\n    deepest = np.where(np.sum(matrix[:, order], axis=0) > 0)[0][0]\n    if deepest.size == 0:\n        raise Exception(\"The matrix is empty\")\n\n    # Remove leaves\n    order_ = np.delete(order, [range(0, deepest)])\n\n    for i in order_:\n        # Get direct children\n        children = np.where(ont.dag[:, i] != 0)[0]\n        if children.size > 0:\n            cols = np.concatenate((children, [i]))\n            if mode == 'max':\n                matrix[:, i] = matrix[:, cols].max(axis=1)\n            elif mode == 'fill':\n                rows = np.where(matrix[:, i] == 0)[0]\n                if rows.size:\n                    idx = np.ix_(rows, cols)\n                    if flag_correct_metric_computation_bug_found_by_Anton:\n                        # Corrected way (see https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/420241 )\n                        matrix[rows, i] = matrix[idx].max(axis=1) #  matrix[idx].max(axis=1)[0] # Correction: https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/420241\n                    else:\n                        # Old way - not corrected\n                        matrix[rows, i] = matrix[idx].max(axis=1)[0] # Correction: https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/420241\n                        \n    return\n\n\ndef obo_parser(obo_file, valid_rel=(\"is_a\", \"part_of\")):\n    \"\"\"\n    Parse a OBO file and returns a list of ontologies, one for each namespace.\n    Obsolete terms are excluded as well as external namespaces.\n    \"\"\"\n    term_dict = {}\n    term_id = None\n    namespace = None\n    name = None\n    term_def = None\n    alt_id = []\n    rel = []\n    obsolete = True\n    with open(obo_file) as f:\n        for line in f:\n            line = line.strip().split(\": \")\n            if line and len(line) > 1:\n                k = line[0]\n                v = \": \".join(line[1:])\n                if k == \"id\":\n                    # Populate the dictionary with the previous entry\n                    if term_id is not None and obsolete is False and namespace is not None:\n                        term_dict.setdefault(namespace, {})[term_id] = {'name': name,\n                                                                       'namespace': namespace,\n                                                                       'def': term_def,\n                                                                       'alt_id': alt_id,\n                                                                       'rel': rel}\n                    # Assign current term ID\n                    term_id = v\n\n                    # Reset optional fields\n                    alt_id = []\n                    rel = []\n                    obsolete = False\n                    namespace = None\n\n                elif k == \"alt_id\":\n                    alt_id.append(v)\n                elif k == \"name\":\n                    name = v\n                elif k == \"namespace\" and v != 'external':\n                    namespace = v\n                elif k == \"def\":\n                    term_def = v\n                elif k == 'is_obsolete':\n                    obsolete = True\n                elif k == \"is_a\" and k in valid_rel:\n                    s = v.split('!')[0].strip()\n                    rel.append(s)\n                elif k == \"relationship\" and v.startswith(\"part_of\") and \"part_of\" in valid_rel:\n                    s = v.split()[1].strip()\n                    rel.append(s)\n\n        # Last record\n        if obsolete is False and namespace is not None:\n            term_dict.setdefault(namespace, {})[term_id] = {'name': name,\n                                                          'namespace': namespace,\n                                                          'def': term_def,\n                                                          'alt_id': alt_id,\n                                                          'rel': rel}\n    return term_dict\n\n\ndef gt_parser(gt_file, ontologies):\n    \"\"\"\n    Parse ground truth file. Discard terms not included in the ontology.\n    \"\"\"\n    gt_dict = {}\n    with open(gt_file) as f:\n        for line in f:\n            line = line.strip().split()\n            if line:\n                p_id, term_id = line[:2]\n                for ont in ontologies:\n                    if term_id in ont.terms_dict:\n                        gt_dict.setdefault(ont.namespace, {}).setdefault(p_id, []).append(term_id)\n                        break\n\n    gts = {}\n    for ont in ontologies:\n        if gt_dict.get(ont.namespace):\n            matrix = np.zeros((len(gt_dict[ont.namespace]), ont.idxs), dtype='bool')\n            ids = {}\n            for i, p_id in enumerate(gt_dict[ont.namespace]):\n                ids[p_id] = i\n                for term_id in gt_dict[ont.namespace][p_id]:\n                    matrix[i, ont.terms_dict[term_id]['index']] = 1\n            propagate(matrix, ont, ont.order, mode='max')\n            gts[ont.namespace] = GroundTruth(ids, matrix, ont.namespace)\n\n    return gts\n\n\ndef pred_parser(f, ontologies, gts, prop_mode, max_terms=None):\n    \"\"\"\n    Parse a prediction file and returns a list of prediction objects, one for each namespace.\n    If a predicted is predicted multiple times for the same target, it stores the max.\n    This is the slow step if the input file is huge, ca. 1 minute for 5GB input on SSD disk.\n    \"\"\"\n    ids = {}\n    matrix = {}\n    ns_dict = {}  # {namespace: term}\n    onts = {ont.namespace: ont for ont in ontologies}\n    for ns in gts:\n        matrix[ns] = np.zeros(gts[ns].matrix.shape, dtype='float')\n        ids[ns] = {}\n        for term in onts[ns].terms_dict:\n            ns_dict[term] = ns\n\n    for line in f:\n        p_id, term_id, prob = line\n        ns = ns_dict.get(term_id)\n        if ns in gts and p_id in gts[ns].ids:\n            i = gts[ns].ids[p_id]\n            if max_terms is None or np.count_nonzero(matrix[ns][i]) <= max_terms:\n                j = onts[ns].terms_dict.get(term_id)['index']\n                ids[ns][p_id] = i\n                matrix[ns][i, j] = max(matrix[ns][i, j], float(prob))\n\n    predictions = []\n    for ns in ids:\n        if ids[ns]:\n            propagate(matrix[ns], onts[ns], onts[ns].order, mode=prop_mode)\n            predictions.append(Prediction(ids[ns], matrix[ns], len(ids[ns]), ns))\n\n    if not predictions:\n        raise Exception(\"Empty prediction, check format\")\n\n    return predictions\n\n\ndef ia_parser(file):\n    ia_dict = {}\n    with open(file) as f:\n        for line in f:\n            if line:\n                term, ia = line.strip().split()\n                ia_dict[term] = float(ia)\n    return ia_dict\n\n# Computes the root terms in the dag\ndef get_roots_idx(dag):\n    return np.where(dag.sum(axis=1) == 0)[0]\n\n\n# Computes the leaf terms in the dag\ndef get_leafs_idx(dag):\n    return np.where(dag.sum(axis=0) == 0)[0]\n\n\n# Return a mask for all the predictions (matrix) >= tau\ndef solidify_prediction(pred, tau):\n    return pred >= tau\n\n\n# computes the f metric for each precision and recall in the input arrays\ndef compute_f(pr, rc):\n    n = 2 * pr * rc\n    d = pr + rc\n    return np.divide(n, d, out=np.zeros_like(n, dtype=float), where=d != 0)\n\n\ndef compute_s(ru, mi):\n    return np.sqrt(ru**2 + mi**2)\n    # return np.where(np.isnan(ru), mi, np.sqrt(ru + np.nan_to_num(mi)))\n\nimport time\nfrom scipy.sparse import csr_matrix\n\ndef compute_metrics_(tau_arr, g, pred, toi, n_gt, wn_gt=None, ic_arr=None):\n\n    verbose = 0;         \n\n    if verbose >= 10:\n        t0 = time.time()\n    \n    metrics = np.zeros((len(tau_arr), 7), dtype='float')  # cov, pr, rc, wpr, wrc, ru, mi\n\n    if verbose >= 10:\n        print('type(toi), toi', type(toi), toi )\n    tmp = pred.matrix[:, toi]\n    if verbose >= 10:\n        print('type(tmp), tmp.shape', type(tmp), tmp.shape )\n    p_s = csr_matrix(tmp )\n    ic_arr_toi = ic_arr[toi]\n    if verbose >= 10:\n        print('type(ic_arr_toi), ic_arr_toi.shape', type(ic_arr_toi), ic_arr_toi.shape )\n\n    \n    g_s = csr_matrix( g )\n    \n    if verbose >= 10:\n        print( 'csr_matrix done %.1f'%(time.time( ) - t0 ), 'p_s.shape, g_s.shape:', p_s.shape, g_s.shape )    \n\n\n    \n    for i, tau in enumerate(tau_arr):\n\n#         p = solidify_prediction(pred_matrix_toi, tau)\n\n#         # number of proteins with at least one term predicted with score >= tau\n#         metrics[i, 0] = (p.sum(axis=1) > 0).sum()\n\n#         # Terms subsets\n#         intersection = np.logical_and(p, g)  # TP                                # SLOW PART !!!\n\n#         # Subsets size\n#         n_pred = p.sum(axis=1)\n#         n_intersection = intersection.sum(axis=1)\n\n#         # Precision, recall\n#         metrics[i, 1] = np.divide(n_intersection, n_pred, out=np.zeros_like(n_intersection, dtype='float'),\n#                                   where=n_pred > 0).sum()\n#         metrics[i, 2] = np.divide(n_intersection, n_gt, out=np.zeros_like(n_gt, dtype='float'), where=n_gt > 0).sum()\n\n\n        if verbose >= 100:\n            t0 = time.time()\n            print()\n            print(i, tau, 'Start %.1f'%(time.time( ) - t0 ) )\n        p = p_s > tau # solidify_prediction(p, tau)\n        if verbose >= 100:\n            print(i, tau, 'solidify done %.1f'%(time.time( ) - t0 ), 'p.shape:', p.shape,  )\n#         p_s = csr_matrix(p)\n#         print(i, tau, 'csr_matrix done %.1f'%(time.time( ) - t0 ), 'p.shape:', p.shape,  )\n        \n#         print(p.shape)\n        # number of proteins with at least one term predicted with score >= tau\n        metrics[i, 0] = (p.sum(axis=1) > 0).sum()\n\n        # Terms subsets\n#         intersection = np.logical_and(p, g)  # TP\n        intersection = p.multiply( g_s)  # TP\n\n\n        if ic_arr is not None:\n            \n            # Weighted precision, recall\n            # wn_pred = (p * ic_arr_toi).sum(axis=1) # \n#             wn_pred = np.dot(p , ic_arr_toi)\n            wn_pred = p.dot( ic_arr_toi)\n            # wn_intersection = (intersection *ic_arr_toi).sum(axis=1)\n#             wn_intersection = np.dot( intersection , ic_arr_toi )\n            wn_intersection =  intersection.dot( ic_arr_toi )\n            \n            if verbose >= 100:\n                print(i, tau, 'After w_pred wn_intersection  %.1f'%(time.time( ) - t0 ) )\n            \n            metrics[i, 3] = np.divide(wn_intersection, wn_pred, out=np.zeros( wn_intersection.shape, dtype='float'),\n                                      where=wn_pred > 0).sum()\n            metrics[i, 4] = np.divide(wn_intersection, wn_gt, out=np.zeros(wn_intersection.shape, dtype='float'),\n                                      where=n_gt > 0).sum()\n            if verbose >= 100:\n                print(i, tau, 'After metrics 3,4   %.1f'%(time.time( ) - t0 ) )\n\n#             # Terms subsets\n#             remaining = np.logical_and(np.logical_not(p), g)  # FN --> not predicted but in the ground truth\n#             mis = np.logical_and(p, np.logical_not(g))  # FP --> predicted but not in the ground truth\n\n#             print(i, tau, 'After remining and miss  %.1f'%(time.time( ) - t0 ) )\n            \n#             # Misinformation, remaining uncertainty\n#             metrics[i, 5] = (remaining * ic_arr_toi).sum(axis=1).sum()\n#             metrics[i, 6] = (mis * ic_arr_toi).sum(axis=1).sum()\n\n    return metrics\n\n\ndef compute_metrics(pred, gt, toi, tau_arr, ic_arr=None, n_cpu=0):\n    \"\"\"\n    Takes the prediction and the ground truth and for each threshold in tau_arr\n    calculates the confusion matrix and returns the coverage,\n    precision, recall, remaining uncertainty and misinformation.\n    Toi is the list of terms (indexes) to be considered\n    \"\"\"\n    g = gt.matrix[:, toi]\n    n_gt = g.sum(axis=1)\n    wn_gt = None\n    if ic_arr is not None:\n        wn_gt = (g * ic_arr[toi]).sum(axis=1)\n\n    # Parallelization\n    if n_cpu == 0:\n        n_cpu = mp.cpu_count()\n\n    arg_lists = [[tau_arr, g, pred, toi, n_gt, wn_gt, ic_arr] for tau_arr in np.array_split(tau_arr, n_cpu)]\n    if 0:\n        # Original parallel way (# It does not work on Kaggle)\n        arg_lists = [[tau_arr, g, pred, toi, n_gt, wn_gt, ic_arr] for tau_arr in np.array_split(tau_arr, n_cpu)]\n        with mp.Pool(processes=n_cpu) as pool:\n            metrics = np.concatenate(pool.starmap(compute_metrics_, arg_lists), axis=0)\n    else: \n        # no-parallel: \n        metrics = compute_metrics_(tau_arr, g, pred, toi, n_gt, wn_gt, ic_arr )\n\n    return pd.DataFrame(metrics, columns=[\"cov\", \"pr\", \"rc\", \"wpr\", \"wrc\", \"ru\", \"mi\"])\n\n\ndef evaluate_prediction(prediction, gt, ontologies, tau_arr, normalization='cafa', n_cpu=0):\n    dfs = []\n    for p in prediction:\n        ns = p.namespace\n        ne = np.full(len(tau_arr), gt[ns].matrix.shape[0])\n\n        ont = [o for o in ontologies if o.namespace == ns][0]\n\n        # cov, pr, rc, wpr, wrc, ru, mi\n        metrics = compute_metrics(p, gt[ns], ont.toi, tau_arr, ont.ia, n_cpu)\n\n        for column in [\"pr\", \"rc\", \"wpr\", \"wrc\", \"ru\", \"mi\"]:\n            if normalization == 'gt' or (column in [\"rc\", \"wrc\"] and normalization == 'cafa'):\n                metrics[column] = np.divide(metrics[column], ne, out=np.zeros_like(metrics[column], dtype='float'), where=ne > 0)\n            else:\n                metrics[column] = np.divide(metrics[column], metrics[\"cov\"], out=np.zeros_like(metrics[column], dtype='float'), where=metrics[\"cov\"] > 0)\n\n        metrics['ns'] = [ns] * len(tau_arr)\n        metrics['tau'] = tau_arr\n        metrics['cov'] = np.divide(metrics['cov'], ne, out=np.zeros_like(metrics['cov'], dtype='float'), where=ne > 0)\n        metrics['f'] = compute_f(metrics['pr'], metrics['rc'])\n        metrics['wf'] = compute_f(metrics['wpr'], metrics['wrc'])\n        metrics['s'] = compute_s(metrics['ru'], metrics['mi'])\n\n        dfs.append(metrics)\n\n    return pd.concat(dfs)\n\n# Tau array, used to compute metrics at different score thresholds\nth_step = 0.01\ntau_arr = np.arange(0.01, 1, th_step)\n#Consider terms without parents, e.g. the root(s), in the evaluation\nno_orphans = False\n# Parse and set information accretion (optional)\nia_dict = ia_parser('/kaggle/input/cafa-5-protein-function-prediction/IA.txt')\n\n# Parse the OBO file and creates a different graph for each namespace\nontologies = []\nobo_file = '/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo'\nfor ns, terms_dict in obo_parser(obo_file).items():\n    ontologies.append(Graph(ns, terms_dict, ia_dict, not no_orphans))","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:30.337125Z","iopub.execute_input":"2023-08-12T19:45:30.338105Z","iopub.status.idle":"2023-08-12T19:45:34.927441Z","shell.execute_reply.started":"2023-08-12T19:45:30.338062Z","shell.execute_reply":"2023-08-12T19:45:34.92611Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Wrapper to call CAFA5 metric computation\n\nFunction to call CAFA5 metric computation from current notebook environment.\n\nPay attention - that metric computation is very slow and RAM consuming - be careful ! \n\nIt is quite technical - no need to go into details - just use as a blacbox.\n\nThanks to Sergei Fironov https://www.kaggle.com/code/sergeifironov/validate-ridge  - please upvote his work.\nWe are based on his code. \n\nhttps://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/420241 - Anton Vakhrushev - correcting error in the initial code of the metric computation - please upvote \n","metadata":{}},{"cell_type":"code","source":"%%time\nimport os.path\n\n######################################################################################3\n###############  Load trainTerms\n######################################################################################3\n\nprint()\nfn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv'\nprint(fn)\ntrainTerms = pd.read_csv(fn, sep=\"\\t\")\nprint(trainTerms.shape)\nprint('trainTerms memory_usage Mb:', trainTerms.memory_usage().sum()/1e6  )\ndisplay(trainTerms.head(3))\n\ndef get_F1_etc_scores_official_CAFA_evaluation( Y_pred, IX, cutoff_threshold_low = 0.01,    make_plots = True ,  verbose = 0 ): \n    '''\n    Computation of F1-weighted scores are called here. \n    Here we prepare Y_pred, Y in format required by functions provided by organizers - see github: https://github.com/BioComputingUP/CAFA-evaluator\n    Y_pred  -  predictions\n    IX  - indices selecting part which correspond to Y_pred in full Y \n    Params:\n    cutoff_threshold_low - predictions lower (strictly) will be dropped (effectively set to zero)\n        (!) higher cutoff_threshold_low will improve  both RAM/speed. For most models 0.1 is Okay, and even 0.18. Consider using 0.1-0.15.   \n    make_plots = True  - create plots of  F1,precision,recall (weighted) depending on threshold  \n    verbose = 0\n    \n    Function uses external variables: \n    trainTerms - training labels provided by orgs: /kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv\n    vec_train_protein_ids - ids of the proteins in the current train - should correspond to \"X\" - features \n    ontologies - data from: /kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo\n    tau_arr - array of thresholds \n    \n    Note: pay attention - Y_true is NOT an input argument. (Unusual for metric computation). \n    We will create Y_true here from \"trainTerms\" cutting only those proteins which correspond to Y_pred indexes \"IX\".\n    Thus it is highly important pass here the correct \"IX\" i.e. corresponding to Y_pred, otherwise results will not be correct\n    '''\n\n    t00 = time.time()\n    if verbose >= 100:\n        print('Scoring starts. n_samples:', len(IX) )\n\n    ##########################################################################################\n    # Prepare \"ground truth\" - \"gt\" terms(labels) in required format  \n    ##########################################################################################\n\n    # First save to file, because function \"gt_parser\" works with files as input \n    # Only part corresponding to providex indices IX will be generated \n    t0 = time.time()\n    trainTerms[ trainTerms.EntryID.isin(vec_train_protein_ids[IX]) ].to_csv('valid.tsv', sep='\\t', index=False) # Wall time: 4.11 s  for 28k samples\n    if verbose >= 1000:\n        print('save valid.csv %.1f'%(time.time() - t0 )) \n\n    # Prepare \"gt\" labels \n    t0 = time.time()\n    gt = gt_parser('valid.tsv', ontologies) # Wall time: 1min 22s  for 28k samples\n    if verbose >= 100:\n        print('gt_parser %.1f'%(time.time() - t0 ))\n\n    ##########################################################################################\n    # prepare predicitons as list of triples - (protein, term(label), prediction) \n    ##########################################################################################\n\n    t0 = time.time()\n    vec_train_protein_ids_loc = vec_train_protein_ids[IX]\n    preds = []\n    for i in range(len(vec_train_protein_ids_loc)):\n        for j in range(len(labels_to_consider)):\n            if Y_pred[i,j] >= cutoff_threshold_low:\n                preds.append((vec_train_protein_ids_loc[i], \n                              labels_to_consider[j],\n                              Y_pred[i,j]                        ))\n    if verbose >= 1000:            \n        print('create preds %.1f'%(time.time() - t0 ))       \n\n    ##########################################################################################\n    # Parse predictions - propagation happens here  \n    ##########################################################################################\n    t0 = time.time()\n    preds = pred_parser(preds, ontologies, gt, prop_mode='fill', max_terms=500) # \n    if verbose >= 1000:            \n        print('pred_parser %.1f'%(time.time() - t0 ), 'len(preds)', len(preds) )            \n\n    gc.collect()\n\n    ##########################################################################################\n    # Main scores calculations happends here: \n    ##########################################################################################\n    # %%time\n    t0 = time.time()\n    df_metrics = evaluate_prediction(preds, gt, ontologies, tau_arr, n_cpu=1) # Wall time: 37.7 s for 28k samples\n    if verbose >= 1000:            \n        print('evaluate_prediction %.1f'%(time.time() - t0 ), 'got df_metrics with shape:', df_metrics.shape )            \n    if verbose >= 10000:            \n        display( df_metrics.head(2) )\n\n        \n    ##########################################################################################\n    # Comptutations finished. Below are optional plots, output preparartions etc.  \n    ##########################################################################################\n    \n    ##########################################################################################\n    ##########################################################################################\n    ##########################################################################################\n    ##########################################################################################\n    ##########################################################################################\n    \n    \n    if verbose >= 100:\n        _t = df_metrics.groupby('ns').agg({'wf':'max'})\n        display( _t )\n        print( _t.mean() ) \n\n    if verbose >= 100:\n        print('F1-scoring finished. %.1f secs passed'%(time.time() - t00 ))\n\n    # %%time\n    if make_plots:\n        try:\n            list_uv = list(df_metrics['ns'].unique() )\n            #print(list_uv)\n            fig = plt.figure(figsize = (20,4))\n            i0 = 0;\n            for  col in  ['wf', 'wpr', 'wrc' ] :\n                i0+=1\n    #             print(i0,col)\n                fig.add_subplot(1,3,i0)\n\n                for uv in list_uv:\n                    mask = df_metrics['ns'] == uv\n                    v = df_metrics[mask][col]\n                    plt.plot(v.values, label = uv)\n                plt.title(col, fontsize  = 20)\n                plt.legend()\n                plt.grid()\n            plt.show()        \n        except:\n            print('Exception in plot')\n    \n    ########################################################################################\n    # Prepare output of scores : \n    ########################################################################################\n    _t = {'cellular_component':'CCO', 'biological_process':'BPO','molecular_function':'MFO'}\n    dict_scores_etc = {}\n    df_s = df_metrics.groupby('ns').agg({'wf':'max'})\n    dict_scores_etc['F1w'] = np.round( df_s.mean().iloc[0], 6) \n    for k in _t:\n        k2 = _t[k]\n        # print(k,dict_scores_etc )\n        if k in  df_s.index:\n            dict_scores_etc['F1 '+ k2 ] = np.round( df_s.loc[k].iat[0], 6) \n        else:\n            dict_scores_etc['F1 '+ k2 ] = 0\n\n    ########################################################################################\n    # Prepare output of thresholds : \n    ########################################################################################\n    for k in _t:\n        k2 = _t[k]\n        m = df_metrics['ns'] == k\n        if m.sum()>0:\n            IX = df_metrics[m]['wf'].argmax()\n            thres_optimal = df_metrics[m]['tau'].iat[IX]\n            dict_scores_etc['thres '+ k2 ] = thres_optimal\n        else:\n            dict_scores_etc['thres '+ k2 ] = 0\n            \n    dict_scores_etc['F-Scores Time'] = np.round( time.time() - t00   ,1)         \n    if verbose >= 100:\n        print('Scores: ', dict_scores_etc)        \n\n    if os.path.isfile('valid.tsv') :\n        os.remove('valid.tsv')\n        \n    return  dict_scores_etc   ","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:34.929073Z","iopub.execute_input":"2023-08-12T19:45:34.929542Z","iopub.status.idle":"2023-08-12T19:45:38.48513Z","shell.execute_reply.started":"2023-08-12T19:45:34.929498Z","shell.execute_reply":"2023-08-12T19:45:38.483827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Auxilliary functions compute/save scores, etc ","metadata":{}},{"cell_type":"code","source":"%%time\nimport time \nimport gc\n\ndef update_modeling_stat( df_stat, Y_pred,  Y,  flag_compute_cafa_f1 = False ,  str_model_id = '',  dict_optional_info = {}, verbose = 0):\n    '''\n    Compute/store/save scores/metrics/statistics on modelling.\n    '''\n    if verbose >= 100:\n        print('Scoring starts')          \n\n    list_folds_ix =  np.sort(list ( set(folds))  )                \n    for ix_fold  in  list_folds_ix[:n_folds_to_process]:\n        t0 = time.time()\n        IX_df_stat = len(df_stat)+1\n        mask_fold = folds == ix_fold\n        IX_loc = np.where(mask_fold >  0)[0]; \n        df_stat.loc[IX_df_stat,'Model'] = str_model_id\n        df_stat.loc[IX_df_stat,'Fold'] = ix_fold\n\n        #from torchmetrics import AUROC as torch_AUCROC\n        #torch_auroc = torch_AUCROC(task = 'binary')\n#             s = torch_auroc(Y_pred,  Y) # auroc between flattened arguments \n#             df_stat.loc[IX_df_stat, 'AUC'] = np.round( s.item(),5)\n        from sklearn.metrics import roc_auc_score\n        s = roc_auc_score(Y[IX_loc,:].ravel(), Y_pred[IX_loc,:].ravel() )\n        df_stat.loc[IX_df_stat, 'AUC'] = np.round( s,5)\n\n        \n        ####################################################################\n        #### Call CAFA-F1 computation - slow and RAM consuming - be careful \n        ####################################################################\n        if flag_compute_cafa_f1:\n            IX_for_cafa_f1 = IX_loc\n            if str(mode_downsample_validation_for_cafa_f1).startswith('top_'):\n                nn = int(mode_downsample_validation_for_cafa_f1.split('_')[1]) # 'top_300'\n                IX_for_cafa_f1 = IX_loc[:nn]\n            elif str(mode_downsample_validation_for_cafa_f1).startswith('random_'):\n                nn = int(mode_downsample_validation_for_cafa_f1.split('_')[1]) # 'top_300'\n                IX_for_cafa_f1 = np.random.permutation(IX_loc)[:nn]\n            dict_scores_etc = get_F1_etc_scores_official_CAFA_evaluation( Y_pred[IX_for_cafa_f1,:], IX_for_cafa_f1,  cutoff_threshold_low = cutoff_threshold_low ,\n                                                                         make_plots = True ,  verbose = 10000 )\n            for k in dict_scores_etc:\n                df_stat.loc[IX_df_stat,k] = dict_scores_etc[k]\n\n            for t in [0.2, 0.3,0.4,0.5]:\n                _c = (Y_pred >= t).ravel().sum()\n                df_stat.loc[IX_df_stat,'GE%.1f per prot'%(t)] = np.round(_c/Y.shape[0])\n        \n            log_available_ram(f'CAFA-F1 scoring fold {ix_fold} finished. Model {str_model_id}')\n        \n        df_stat.loc[IX_df_stat, 'n_targets'] = Y.shape[1]\n        df_stat.loc[IX_df_stat, 'n_samples Val'] = Y.shape[0]       \n        df_stat.loc[IX_df_stat, 'Time scoring'] = np.round(time.time() - t0, 1 )    \n        for k in dict_optional_info:\n            val = dict_optional_info[k]\n            df_stat.loc[IX_df_stat, k] = val   \n        \n        try:\n            df_stat.loc[IX_df_stat, 'Features'] = str( list_features_id )    \n            df_stat.loc[IX_df_stat, 'Model Features'] = str_model_id.split(' ')[-1] + ' ' + str( list_features_id )\n            \n        except:\n            pass\n        \n        \n        torch.cuda.empty_cache()\n        gc.collect()\n        df_stat.to_csv('df_stat.csv')     \n        \n        \n    if verbose > 0:\n        display(df_stat.tail(n_folds_to_process))\n        print('Scoring finished. Seconds passed:  %.1f'%(time.time() - t0)  )\n\n    return df_stat\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:38.487182Z","iopub.execute_input":"2023-08-12T19:45:38.487956Z","iopub.status.idle":"2023-08-12T19:45:38.510765Z","shell.execute_reply.started":"2023-08-12T19:45:38.487911Z","shell.execute_reply":"2023-08-12T19:45:38.509402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:38.512514Z","iopub.execute_input":"2023-08-12T19:45:38.512961Z","iopub.status.idle":"2023-08-12T19:45:38.522157Z","shell.execute_reply.started":"2023-08-12T19:45:38.512928Z","shell.execute_reply":"2023-08-12T19:45:38.520646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train the model","metadata":{}},{"cell_type":"code","source":"%%time\nimport gc \nimport time \nimport datetime\n# current_datetime = datetime.datetime.now();  print(\"Current datetime:\", current_datetime)\nlog_available_ram('Before modeling')","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:38.52604Z","iopub.execute_input":"2023-08-12T19:45:38.526491Z","iopub.status.idle":"2023-08-12T19:45:38.535952Z","shell.execute_reply.started":"2023-08-12T19:45:38.526455Z","shell.execute_reply":"2023-08-12T19:45:38.534773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n######################### Params ##################################################3\n\nif 0: # Code updated - now most params are specified for each model or in the top section \"Key params\" \n    mode_submit = True \nverbose = 0\n\n\n######################### Output ##################################################\ndf_stat = pd.DataFrame()\n\nif (mode_submit is not None) and ( mode_submit != False ):\n    Y_submit = np.zeros( (141865, Y.shape[1] )  , dtype = np.float16 )  # Predictions for submission will be stored here \n    # Results from all models and all folds will be blended \n    print('Y_submit mbytes:', Y_submit.nbytes/1024/1024)\ncnt_blend_submit = 0 ;  \n\nif  flag_compute_oof_predictions:\n    Y_pred_oof_blend  = np.zeros( ( Y.shape )  , dtype = np.float16 )\n    print('Y_pred_oof_blend mbytes:', Y_pred_oof_blend.nbytes/1024/1024)\ncnt_blend_oof = -1;\n\n\n########################## Preparations ###########################################\nlog_available_ram('Right before modeling')\n\nif flag_compute_stat_for_each_model:  # Predictions OOF for each particular model - will be rewritten for each modelling \n    Y_pred_oof = np.zeros( ( Y.shape )  , dtype = np.float16 )\n    print('Y_pred_oof mbytes:', Y_pred_oof.nbytes/1024/1024)\n\ni_model = -1 # \ni_config = -1 # conter for configurations \ni_config_processed = -1\nt0modeling = time.time()\nlist_folds_ix =  np.sort(list ( set(folds))  )\nprint(); print('Start training models',datetime.datetime.now()) ; print()\n########################## Main modelling  ###########################################\nfor main_config_model_feature_etc  in list_main_config_model_feature_etc:\n    i_config += 1 \n    model_config = main_config_model_feature_etc['model']\n    if ('Keras' in model_config.keys() ) and ( model_config['Keras'] ): continue   # Keras models will be processed in the next cell - RAM leak problem\n    i_config_processed += 1\n    if 'list_features_id' in main_config_model_feature_etc.keys():\n        print()\n        X,vec_train_protein_ids,X_submit,  submit_protein_ids = get_features(main_config_model_feature_etc['list_features_id'], verbose = 100)\n        gc.collect()\n        log_available_ram(f\"New features loaded:  {str(main_config_model_feature_etc['list_features_id'])}\" )  \n        print()\n        \n    mode_downsample_train = model_config.get('mode_downsample_train', mode_downsample_train_default)\n    \n    n_selfblend = model_config.get( 'n_selfblend' , 1)\n    if verbose >= 100:\n        print(); print('Starting model_config:', model_config, f'time from start: {(time.time() - t0modeling ):.1f}' )\n    for i_selfblend in range( n_selfblend ): # train-predict same model several times and blend predictions - especially useful for NN, but do not fix random seed (!)\n        i_model += 1 # Models count\n        t0one_model_all_folds = time.time()\n        for ix_fold  in  list_folds_ix[:n_folds_to_process]:\n        \n            model, str_model_id = get_model(model_config)\n            str_model_id_pure_save =  str_model_id\n            str_model_id = str( i_model) + ' ' + str_model_id\n            if n_selfblend > 1:  str_model_id += ' ' + str( i_selfblend )\n\n            ##################### Prepare train data ###################################################\n            mask_fold = folds == ix_fold\n            IX_train = np.where(mask_fold ==  0)[0]; \n            # IX_train = [ix for ix in IX_train if ix in  set_allowed_train_indexes]\n            IX_train = get_downsampled_IX_train(IX_train, mode_downsample_train )\n            X_train = X[IX_train,:]; Y_train = Y[IX_train,:]\n\n            if verbose >= 10:\n                print(f'fold {ix_fold}, model: {str_model_id},  X_train.shape: {X_train.shape}, Y_train.shape: {Y_train.shape}, time: { (time.time() - t0modeling):12.1f} ')\n                print('X_train Mbytes:', X_train.nbytes/1024/1024, 'Y_train Mbytes:', Y_train.nbytes/1024/1024,  )\n            ##################### Call train model ###################################################\n            t0 = time.time()\n            model_fit(model , X_train, Y_train, model_config , str_model_id, verbose = 0 )   \n            time_fit = time.time() - t0\n            if verbose >= 1000:\n                print(f'time_fit {time_fit:.1f}' ) \n            del X_train, Y_train\n            torch.cuda.empty_cache()\n            gc.collect()\n            log_available_ram(f'After Model fit. ix_fold {ix_fold}, i_selfblend, {i_selfblend}, {str_model_id}')\n            \n            ##################### Compute predictions for submission and blend with the previous one ###################################################\n            if  mode_submit : \n                t0 = time.time()\n                Y_submit = (Y_submit * cnt_blend_submit  + model_predict(model , X_submit,  model_config , str_model_id , verbose = 0 ) )/ (cnt_blend_submit + 1);  # Average predictions from different folds/models\n                cnt_blend_submit += 1 \n                time_pred_submit = time.time() - t0\n                torch.cuda.empty_cache()\n                gc.collect()\n\n            #flag_compute_oof_predictions = True                 \n            if  flag_compute_oof_predictions:\n                t0 = time.time()\n                IX_val = np.where(mask_fold > 0 )[0]; \n                X_val = X[IX_val,:];#  Y_val = Y[IX_val,:]\n                Y_pred_val = model_predict(model , X_val,  model_config , str_model_id , verbose = 0 )\n                time_pred_val = time.time() - t0\n                if verbose >= 10000:\n                    print('Y_pred_val.shape', Y_pred_val.shape, f'time_pred_val {time_pred_val:.1f}')\n                    \n                if ix_fold == 0: cnt_blend_oof += 1 \n                Y_pred_oof_blend[IX_val,:] = (Y_pred_oof_blend[IX_val,:] * cnt_blend_oof  + Y_pred_val )/ (cnt_blend_oof + 1); \n            \n                if  flag_compute_stat_for_each_model:\n                    Y_pred_oof[IX_val,:] = (Y_pred_val ) \n                    \n                del X_val, Y_pred_val                    \n                torch.cuda.empty_cache()\n                gc.collect()\n                \n            del model \n            torch.cuda.empty_cache()\n            gc.collect()\n            log_available_ram(f'At fold end. Model {str_model_id}, i_selfblend {i_selfblend}, fold {ix_fold} ')\n\n        time_one_model = np.round( time.time() - t0one_model_all_folds )\n        if flag_compute_stat_for_each_model and flag_compute_oof_predictions: \n            update_modeling_stat(df_stat, Y_pred_oof,  Y, flag_compute_cafa_f1 = flag_compute_cafa_f1_for_each_model , \n                                 str_model_id = str_model_id, dict_optional_info = {'Time': time_one_model , 'i_selfblend':i_selfblend,\n                                'ModelID Pure':str_model_id_pure_save, 'i_config':i_config}, verbose = 0)\n            torch.cuda.empty_cache()\n            gc.collect()\n            log_available_ram(f'After OOF-Stat Calculation. Model {str_model_id}, i_selfblend {i_selfblend}' )\n            \n        if flag_compute_each_blend_stat and flag_compute_oof_predictions:\n            update_modeling_stat(df_stat, Y_pred_oof_blend,  Y, flag_compute_cafa_f1 = flag_compute_cafa_f1_for_each_blend , \n                                 str_model_id = str(cnt_blend_oof )+ ' Blend'+ ' ' +str_model_id, dict_optional_info = {'Time': time_one_model, \n                                        'Blend': cnt_blend_oof , 'i_selfblend':i_selfblend}, verbose = 0)\n            torch.cuda.empty_cache()\n            gc.collect()\n            log_available_ram(f'After Blend-Stat Calculation. Model {str_model_id}, i_selfblend {i_selfblend}' )\n                \nif (i_config_processed >= 0) and flag_save_numpy_Y_pred_oof_blend and flag_compute_oof_predictions:\n    t0 = time.time()\n    fn = 'Y_pred_oof_blend.npy'\n    np.save(fn,Y_pred_oof_blend)\n    print(f'File {fn} saved. Y_pred_oof_blend.shape: {Y_pred_oof_blend.shape}. Time: {(time.time()-t0):.1f}')\n    t0 = time.time()\n    fn = 'Y_labels.npy'\n    np.save(fn,Y_labels)\n    print(f'File {fn} saved. Time: {(time.time()-t0):.1f}')\n\nif (i_config_processed >= 0) and flag_save_numpy_Y_submit and mode_submit:\n    t0 = time.time()\n    fn = 'Y_submit.npy'\n    np.save(fn,Y_submit)\n    print(f'File {fn} saved. Y_submit.shape: {Y_submit.shape}. Time: {(time.time()-t0):.1f}')\n    t0 = time.time()\n    fn = 'Y_labels.npy'\n    np.save(fn,Y_labels)\n    print(f'File {fn} saved. Time: {(time.time()-t0):.1f}')\n\n\ndisplay(df_stat)            ","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:45:38.537329Z","iopub.execute_input":"2023-08-12T19:45:38.537669Z","iopub.status.idle":"2023-08-12T19:58:58.073107Z","shell.execute_reply.started":"2023-08-12T19:45:38.537623Z","shell.execute_reply":"2023-08-12T19:58:58.071812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\nlog_available_ram('After Modelling 1 Finished')\n\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:58:58.075118Z","iopub.execute_input":"2023-08-12T19:58:58.075906Z","iopub.status.idle":"2023-08-12T19:58:58.083515Z","shell.execute_reply.started":"2023-08-12T19:58:58.075861Z","shell.execute_reply":"2023-08-12T19:58:58.08253Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Modeling specific to Keras models\n\nKeras models leak RAM in the code above - so here we try to create a specific  code \nfor modelling with Keras \n\nThe cell below is almost same as the one above, but we simplify things specifically to Keras models setup\n\nThat is not intended technical modifaction - if we will find the ways to avoid memory leak in Keras models will return to uniform code as in the previous cell. \n","metadata":{}},{"cell_type":"code","source":"%%time\n######################### Params ##################################################3\n\nverbose = 1000\n\n######################### Output ##################################################\n\n########################## Preparations ###########################################\nlog_available_ram('Right before Keras modeling')\n\n\ni_config_keras = -1\nt0modeling = time.time()\nlist_folds_ix =  np.sort(list ( set(folds))  )\nprint(); print('Start training Keras models',datetime.datetime.now()) ; print()\n########################## Main modelling  ###########################################\nfor main_config_model_feature_etc  in list_main_config_model_feature_etc:\n    model_config = main_config_model_feature_etc['model']\n    if 'Keras' not in model_config.keys(): continue  \n    if  model_config['Keras'] == False : continue \n    i_config_processed += 1\n    i_config_keras += 1\n    if 'list_features_id' in main_config_model_feature_etc.keys():\n        print()\n        X,vec_train_protein_ids,X_submit,  submit_protein_ids = get_features(main_config_model_feature_etc['list_features_id'], verbose = 100)\n        gc.collect()\n        log_available_ram(f\"New features loaded:  {str(main_config_model_feature_etc['list_features_id'])}\" )  \n        print()\n    \n    mode_downsample_train = model_config.get('mode_downsample_train', mode_downsample_train_default)\n        \n    n_selfblend = model_config.get( 'n_selfblend' , 1)\n    if verbose >= 100:\n        print(); print('Starting model_config:', model_config, f'time from start: {(time.time() - t0modeling ):.1f}' )\n    for i_selfblend in range( n_selfblend ): # train-predict same model several times and blend predictions - especially useful for NN, but do not fix random seed (!)\n        i_model += 1 # Models count\n        t0one_model_all_folds = time.time()\n        for ix_fold  in  list_folds_ix[:n_folds_to_process]:\n            if verbose >= 10: \n                print('---------------------------------------------- Fold', ix_fold, '----------------------------------------------')\n                \n            model, str_model_id = get_model(model_config)\n            str_model_id_pure_save = str_model_id\n            log_available_ram(f'After Keras model init. ix_fold {ix_fold}, n_selfblend, {n_selfblend}, {str_model_id}')\n            str_model_id = str( i_model) + ' ' + str_model_id\n            if n_selfblend > 1:  str_model_id += ' ' + str( i_selfblend )\n                \n            ##################### Prepare train data ###################################################\n            mask_fold = folds == ix_fold\n            IX_train = np.where(mask_fold ==  0)[0]; \n            # IX_train = [ix for ix in IX_train if ix in  set_allowed_train_indexes]\n            IX_train = get_downsampled_IX_train(IX_train, mode_downsample_train )\n\n            #X_train = X[IX_train,:]; Y_train = Y[IX_train,:]\n            log_available_ram(f'After setting X_train. ix_fold {ix_fold}, n_selfblend, {n_selfblend}, {str_model_id}')\n\n            if verbose >= 10:\n                print(f'fold {ix_fold}, model: {str_model_id},  X_train.shape: {X[IX_train,:].shape}, Y_train.shape: {Y[IX_train,:].shape}, time: { (time.time() - t0modeling):12.1f} ')\n\n            ##################### Call train model ###################################################\n            t0 = time.time()\n            #model_fit(model , X_train, Y_train, model_config , str_model_id, verbose = 0 )   \n            epochs = model_config.get( 'epochs' , 15)\n            batch_size = model_config.get( 'batch_size' , 128 )\n            model.fit( X[IX_train,:], Y[IX_train,:] ,epochs = epochs,   batch_size = batch_size, verbose = 0  )\n            time_fit = time.time() - t0\n            if verbose >= 1000:\n                print(f'time_fit {time_fit:.1f}' ) \n            #del X_train, Y_train\n            torch.cuda.empty_cache()\n            gc.collect()\n            log_available_ram(f'After Model fit. ix_fold {ix_fold}, i_selfblend, {i_selfblend}, {str_model_id}')\n            \n            ##################### Compute predictions for submission and blend with the previous one ###################################################\n            if  mode_submit : \n                t0 = time.time()\n                Y_submit = (Y_submit * cnt_blend_submit  + model_predict(model , X_submit,  model_config , str_model_id , verbose = 0 ) )/ (cnt_blend_submit + 1);  # Average predictions from different folds/models\n                cnt_blend_submit += 1 \n                time_pred_submit = time.time() - t0\n                torch.cuda.empty_cache()\n                gc.collect()\n                log_available_ram(f'After Predict on submit. ix_fold {ix_fold}, i_selfblend, {i_selfblend}, {str_model_id}')\n                \n            #flag_compute_oof_predictions = True                 \n            if  flag_compute_oof_predictions:\n                t0 = time.time()\n                IX_val = np.where(mask_fold > 0 )[0]; \n                #X_val = X[IX_val,:];#  Y_val = Y[IX_val,:]\n                Y_pred_val = model_predict(model , X[IX_val,:],  model_config , str_model_id , verbose = 0 )\n                time_pred_val = time.time() - t0\n                if verbose >= 10000:\n                    print('Y_pred_val.shape', Y_pred_val.shape, f'time_pred_val {time_pred_val:.1f}')\n                if ix_fold == 0: cnt_blend_oof += 1 \n                Y_pred_oof_blend[IX_val,:] = (Y_pred_oof_blend[IX_val,:] * cnt_blend_oof  + Y_pred_val )/ (cnt_blend_oof + 1); \n            \n                if  flag_compute_stat_for_each_model:\n                    Y_pred_oof[IX_val,:] = (Y_pred_val ) \n                    \n                del  Y_pred_val                    \n                torch.cuda.empty_cache()\n                gc.collect()\n                log_available_ram(f'After Predict on OOF. ix_fold {ix_fold}, i_selfblend, {i_selfblend}, {str_model_id}')\n            \n            keras.backend.clear_session()\n            del model\n            gc.collect()\n            keras.backend.clear_session()\n            torch.cuda.empty_cache()\n            gc.collect()\n            log_available_ram(f'At fold end. Model {str_model_id}, i_selfblend {i_selfblend}, fold {ix_fold} ')\n\n        time_one_model = np.round( time.time() - t0one_model_all_folds )\n        if flag_compute_stat_for_each_model and flag_compute_oof_predictions: \n            update_modeling_stat(df_stat, Y_pred_oof,  Y, flag_compute_cafa_f1 = flag_compute_cafa_f1_for_each_model , \n                 str_model_id = str_model_id, dict_optional_info = {'Time': time_one_model, \n                'i_selfblend':i_selfblend, 'ModelID Pure':str_model_id_pure_save, 'i_config':i_config }, verbose = 0)\n            torch.cuda.empty_cache()\n            gc.collect()\n            log_available_ram(f'After OOF-Stat Calculation. Model {str_model_id}, i_selfblend {i_selfblend}' )\n            \n        if flag_compute_each_blend_stat and flag_compute_oof_predictions:\n            update_modeling_stat(df_stat, Y_pred_oof_blend,  Y, flag_compute_cafa_f1 = flag_compute_cafa_f1_for_each_blend , \n                                 str_model_id = str(cnt_blend_oof )+ ' Blend'+ ' ' +str_model_id, dict_optional_info = {'Time': time_one_model, \n                                'Blend': cnt_blend_oof, 'i_selfblend':i_selfblend}, verbose = 0)\n            torch.cuda.empty_cache()\n            gc.collect()\n            log_available_ram(f'After Blend-Stat Calculation. Model {str_model_id}, i_selfblend {i_selfblend}' )\n                \nif  (i_config_keras >= 0) and (i_config_processed >= 0) and  flag_save_numpy_Y_pred_oof_blend and flag_compute_oof_predictions:\n    t0 = time.time()\n    fn = 'Y_pred_oof_blend.npy'\n    np.save(fn,Y_pred_oof_blend)\n    print(f'File {fn} saved. Y_pred_oof_blend.shape: {Y_pred_oof_blend.shape}. Time: {(time.time()-t0):.1f}')\n    t0 = time.time()\n    fn = 'Y_labels.npy'\n    np.save(fn,Y_labels)\n    print(f'File {fn} saved. Time: {(time.time()-t0):.1f}')\n\nif  (i_config_keras >= 0) and (i_config_processed >= 0) and  flag_save_numpy_Y_submit and mode_submit:\n    t0 = time.time()\n    fn = 'Y_submit.npy'\n    np.save(fn,Y_submit)\n    print(f'File {fn} saved. Y_submit.shape: {Y_submit.shape}. Time: {(time.time()-t0):.1f}')\n    t0 = time.time()\n    fn = 'Y_labels.npy'\n    np.save(fn,Y_labels)\n    print(f'File {fn} saved. Time: {(time.time()-t0):.1f}')\n    log_available_ram(f'After Save Y_submit' )\n\n\ndisplay(df_stat)            ","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:58:58.085388Z","iopub.execute_input":"2023-08-12T19:58:58.086235Z","iopub.status.idle":"2023-08-12T19:58:58.162365Z","shell.execute_reply.started":"2023-08-12T19:58:58.086186Z","shell.execute_reply":"2023-08-12T19:58:58.161029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if flag_compute_stat_for_each_model: \n    del Y_pred_oof\n    \ntorch.cuda.empty_cache()    \ngc.collect()\n\nlog_available_ram('After Modelling Keras Finished')\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:58:58.164099Z","iopub.execute_input":"2023-08-12T19:58:58.164572Z","iopub.status.idle":"2023-08-12T19:58:58.702108Z","shell.execute_reply.started":"2023-08-12T19:58:58.164527Z","shell.execute_reply":"2023-08-12T19:58:58.700768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compute metrics for the final blend","metadata":{}},{"cell_type":"code","source":"%%time\nif  (i_config_processed >= 0) and  flag_compute_final_model_stat and  flag_compute_oof_predictions:  \n    #time_one_model = np.round( time.time() - t0one_model_all_folds )\n    update_modeling_stat(df_stat, Y_pred_oof_blend,  Y, flag_compute_cafa_f1 = True , str_model_id= 'Final Blend', dict_optional_info = { }, verbose = 0)\n    gc.collect()\n    log_available_ram('After Final Stat Calculation')\n            \ndisplay(df_stat)            ","metadata":{"execution":{"iopub.status.busy":"2023-08-12T19:58:58.703786Z","iopub.execute_input":"2023-08-12T19:58:58.705891Z","iopub.status.idle":"2023-08-12T19:59:11.537692Z","shell.execute_reply.started":"2023-08-12T19:58:58.705792Z","shell.execute_reply":"2023-08-12T19:59:11.536406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ntorch.cuda.empty_cache()    \ngc.collect()\n\nlog_available_ram('After Final stat computation')\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Show/plot modelling results ","metadata":{}},{"cell_type":"code","source":"%%time\nif len(df_stat) > 1:\n    if ('AUC' in df_stat.columns) and ( 'F1w' in df_stat.columns ):\n        print(); print('Pearson correlations:')\n        display( df_stat[['AUC','F1w']].corr() )\n\n        print(); print('Spearman correlations:')\n        display( df_stat[['AUC','F1w']].corr(method = 'spearman') )\n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_stat.shape)\ndisplay( df_stat )\n# if 'F1w' in df_stat.columns:\n#     print('Top Sorted by F1w')\n#     display( df_stat.sort_values('F1w', ascending = False).head(10) )    \n# if 'AUC' in df_stat.columns:\n#     print('Top Sorted by AUC')\n#     display( df_stat.sort_values('AUC', ascending = False).head(10) )\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plots in particular Fold-averaged ","metadata":{}},{"cell_type":"code","source":"%%time\nif len(df_stat) > 0:\n    print('Fold averaged')\n    df_stat = df_stat.copy()\n    df_stat['RankTmp'] = range(len(df_stat))\n    if 'Model' in df_stat.columns:\n        d = df_stat.groupby('Model').mean()\n        d = d.sort_values('RankTmp') # Make sure ordering is the same as in the intial df_stat\n        for col2 in ['F1w', 'AUC']:\n            if  (col2 not in df_stat.columns ): continue \n            print(col2, 'top  data:', )\n            display( d.sort_values(col2, ascending = False).head(10) )\n            plt.figure(figsize = (20,5) )\n            plt.suptitle(col2 + ' Fold averaged', fontsize = 20 )\n\n            plt.subplot(1,2,1)\n            m = ['Blend' not in t for t in d.index ]\n            plt.plot(d[col2][m],'*-')\n            plt.title('Models (no blend)', fontsize = 20)\n            plt.grid()# b='on')\n            plt.xticks(rotation=90)\n\n            plt.subplot(1,2,2)\n            m = ['Blend' in t for t in d.index ]\n            plt.plot(d[col2][m],'*-')\n            plt.title('Blend', fontsize = 20)\n            plt.grid()# b='on')\n            plt.xticks(rotation=90)\n\n            plt.show()    \n\n            d.sort_values(col2, ascending = False).to_csv(col2+'_sorted_fold_averaged.csv' )\n\n        print()\n        print('Output top and tail sorted data')    \n        for col2 in ['F1w', 'AUC']:\n            if  (col2 not in df_stat.columns ): continue \n            print(col2, 'top 50  data:', )\n            display( d.sort_values(col2, ascending = False).head(50) )\n            print(col2, 'Tail 50  data:', )\n            display( d.sort_values(col2, ascending = False).tail(50) )\n            print(); print(); print(); print(); print(); print(); \n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Grouping by model type and averageing (i.e. self-blend will be averaged) \n\n","metadata":{}},{"cell_type":"code","source":"%%time\nif len(df_stat) > 0:\n    print('Group for each model type')\n    print('Beware of possible grouping models on different input features - which may not be desirable - if so look on i_config grouping - below')\n    print('If all models on the same data - not at problem at all')\n    df_stat['RankTmp'] = range(len(df_stat))\n    if 'ModelID Pure' in df_stat.columns:\n        d = df_stat.groupby('ModelID Pure').mean()\n        d = d.sort_values('RankTmp') # Make sure ordering is the same as in the intial df_stat\n        for col2 in ['F1w', 'AUC']:\n            if  (col2 not in df_stat.columns ): continue \n            print(col2, 'top  data:', )\n            display( d.sort_values(col2, ascending = False).head(10) )\n            plt.figure(figsize = (20,5) )\n            plt.suptitle(col2 + ' Fold averaged', fontsize = 20 )\n\n            plt.subplot(1,2,1)\n            m = ['Blend' not in t for t in d.index ]\n            plt.plot(d[col2][m],'*-')\n            plt.title('Models (no blend)', fontsize = 20)\n            plt.grid()# b='on')\n            plt.xticks(rotation=90)\n\n    #         plt.subplot(1,2,2)\n    #         m = ['Blend' in t for t in d.index ]\n    #         plt.plot(d[col2][m],'*-')\n    #         plt.title('Blend', fontsize = 20)\n    #         plt.grid()# b='on')\n    #         plt.xticks(rotation=90)\n\n            plt.show()    \n\n            d.sort_values(col2, ascending = False).to_csv(col2+'_sorted_groupped_by_model_type.csv' )\n\n        print()\n        print('Output top and tail sorted data')    \n        for col2 in ['F1w', 'AUC']:\n            if  (col2 not in df_stat.columns ): continue \n            print(col2, 'top 50  data:', )\n            display( d.sort_values(col2, ascending = False).head(50) )\n            print(col2, 'Tail 50  data:', )\n            display( d.sort_values(col2, ascending = False).tail(50) )\n            print(); print(); print(); print(); print(); print(); \n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Grouping by i_cofing (all self-blends will be averaged) , different features will be preserved\n\nby the naming will disapped in groupped dataset - may be will update it later ","metadata":{}},{"cell_type":"code","source":"%%time\nprint('Group for each i_config')\ndf_stat = df_stat.copy()\ndf_stat['RankTmp'] = range(len(df_stat))\nif 'i_config' in df_stat.columns:\n    d = df_stat.groupby('i_config').mean()\n    d = d.sort_values('RankTmp') # Make sure ordering is the same as in the intial df_stat\n    for col2 in ['F1w', 'AUC']:\n        if  (col2 not in df_stat.columns ): continue \n        print(col2, 'top  data:', )\n        display( d.sort_values(col2, ascending = False).head(10) )\n        plt.figure(figsize = (20,5) )\n        plt.suptitle(col2 + ' Fold averaged', fontsize = 20 )\n\n        plt.subplot(1,2,1)\n        m = np.ones( len(d) ).astype(bool) # ['Blend' not in t for t in d.index ]\n        plt.plot(d[col2][m],'*-')\n        plt.title('Models (no blend)', fontsize = 20)\n        plt.grid()# b='on')\n        plt.xticks(rotation=90)\n\n#         plt.subplot(1,2,2)\n#         m = ['Blend' in t for t in d.index ]\n#         plt.plot(d[col2][m],'*-')\n#         plt.title('Blend', fontsize = 20)\n#         plt.grid()# b='on')\n#         plt.xticks(rotation=90)\n\n        plt.show()    \n\n        d.sort_values(col2, ascending = False).to_csv(col2+'_sorted_groupped_by_i_config.csv' )\n\n    print()\n    print('Output top and tail sorted data')    \n    for col2 in ['F1w', 'AUC']:\n        if  (col2 not in df_stat.columns ): continue \n        print(col2, 'top 50  data:', )\n        display( d.sort_values(col2, ascending = False).head(50) )\n        print(col2, 'Tail 50  data:', )\n        display( d.sort_values(col2, ascending = False).tail(50) )\n        print(); print(); print(); print(); print(); print(); \n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plots without averaging by folds","metadata":{}},{"cell_type":"code","source":"print('Plots without averaging by folds')\ncol = 'Model'\nfor col2 in ['F1w', 'AUC']:\n    if ( col not in df_stat.columns) or (col2 not in df_stat.columns ): continue \n    print(col2, 'top  data:', )\n    display( df_stat.sort_values(col2, ascending = False).head(5) )\n    plt.figure(figsize = (20,5) )\n    plt.suptitle(col2, fontsize = 20 )\n    \n    plt.subplot(1,2,1)\n    m = ['Blend' not in t for t in df_stat[col]]\n    plt.plot(df_stat[col2][m],'*-')\n    plt.title('Models (no blend)', fontsize = 20)\n    plt.grid()# b='on')\n    plt.xticks(rotation=90)\n\n    \n    plt.subplot(1,2,2)\n    m = ['Blend' in t for t in df_stat[col]]\n    plt.plot(df_stat[col2][m],'*-')\n    plt.title('Blend', fontsize = 20)\n    plt.grid()# b='on')\n    plt.xticks(rotation=90)\n    \n    plt.show()    \n    \n    \nprint()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if 'AUC' in df_stat.columns:\n    plt.plot(df_stat['AUC'].values,'*-')\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ncol = 'AUC' \ncol2 = 'Fold'\nif (col in df_stat.columns) and ( col2 in df_stat.columns ):\n    d2 = df_stat.groupby(col2).mean()\n    display(d2)\n    d2.to_csv('df_stat_folds_mean.csv')\n    plt.plot(d2[col].values,'*-')\n    plt.xlabel(col2,fontsize = 20 )\n    plt.title(col,fontsize = 20)\n    plt.show()\n\ncol2 = 'Model'\nif (col in df_stat.columns) and ( col2 in df_stat.columns ):\n    d3 = df_stat.groupby(col2).mean()\n    display(d3)\n    d3.to_csv('df_stat_models_mean.csv')\n    plt.plot(d3[col].values,'*-')\n    plt.title(col,fontsize = 20)\n    plt.xlabel(col2,fontsize = 20 )\n    plt.show()\n    display()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif len(df_stat) > 0:\n    if 'Model' in df_stat.columns:\n        m = df_stat['Model'] == 'Final Blend'\n        if m.sum() == 0:\n            m = ( df_stat['Model'] == 'Final' )\n        if m.sum() == 0:\n            m =  ( df_stat['Model'] == 'FinalKeras Blend' )\n        if m.sum() == 0:\n            m = ( df_stat['Model'] == 'Final1 Blend' )\n\n\n        display(df_stat[m] )\n        display(df_stat[m].mean() )\n        if 'F1w' in df_stat.columns:\n            print( df_stat[m].mean().loc['F1w']  )\n        df_stat[m].mean().to_csv('final_means.csv')\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prepare submission","metadata":{}},{"cell_type":"code","source":"%%time\nimport gc\nif 0:\n    del X_train,Y_train, X_val, Y_val, train_dataset, train_dataloader, X, Y\ngc.collect()\ntorch.cuda.empty_cache()\n\nlog_available_ram()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare submission tsv file","metadata":{}},{"cell_type":"code","source":"%%time\n\nmode_submit_prepare = 'slow_less_RAM_consuming'\n# 'slow_less_RAM_consuming' - works slower but consumes less RAM\n\nimport time \nt0 = time.time()\n\nprint( mode_submit_prepare , mode_submit)\nif (mode_submit ) and ( flag_save_final_submit_file ):\n    print(Y_submit.shape)\n    if mode_submit_prepare == 'slow_less_RAM_consuming':\n\n\n        file_path = \"submission.tsv\"\n        cc = 0\n        cc2 = 0\n        with open(file_path, 'w') as file:\n            for i in range(Y_submit.shape[0]):\n                for j in range(Y_submit.shape[1]):\n                    val = Y_submit[i,j]\n                    if val >= cutoff_threshold_low:\n                        str_go_term = str(Y_labels[j])\n                        str_protein_id = str( submit_protein_ids[i] )\n                        str_save = str_protein_id+'\\t'+str_go_term + '\\t' + '%.3f'%val + '\\n'\n                        file.write(str_save)   \n                        cc2 +=1\n                        if cc2 <= 10:\n                            if cc2 == 1: print('First 10 examples of the saved data:')\n                            print(str_save)\n                    cc += 1\n                    if cc % 30_000_000  == 0: \n                        sz = Y_submit.shape[0]*Y_submit.shape[1]\n                        print(cc, 'out of',sz, 'percent %.2f'%(cc/sz*100), 'saved:'  ,cc2, 'time %.1f'%(time.time() - t0 ))\n\n        print(cc2,'results saved to submission file', 'time  %.1f'%(time.time() - t0 )  )                \n\n    else:\n\n        # That is widely used way to preparase submission , but it might crash by RAM \n\n        df_submission = pd.DataFrame(columns = ['Protein Id', 'GO Term Id','Prediction'])\n\n        n_targets_predicted = Y_submit.shape[1]\n        n_samples_predicted = Y_submit.shape[0]\n        print('n_samples_predicted, n_targets_predicted',  n_samples_predicted, n_targets_predicted )\n\n\n        protein_list = []\n        for k in list(submit_protein_ids):\n            protein_list += [k] * n_targets_predicted\n        df_submission['Protein Id'] = protein_list\n\n        df_submission['GO Term Id'] = list(Y_labels) * n_samples_predicted\n        df_submission['Prediction'] = Y_submit.ravel()\n\n        df_submission = df_submission.round(3)\n        df_submission = df_submission[ df_submission['Prediction'] >= cutoff_threshold_low  ]\n\n        memory_usage_per_column = df_submission.memory_usage(deep=True)\n        total_memory_usage = memory_usage_per_column.sum()\n        print(\"\\nTotal memory usage:\", total_memory_usage/1e6, \"Megabytes\")\n\n        print(df_submission.shape)\n        display(df_submission)\n\n        import gc\n        if 0:\n            del preds \n\n        gc.collect()\n\n        df_submission.to_csv(\"submission.tsv\",header=False, index=False,sep='\\t')\n\n\n    log_available_ram('After saving submission')    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plot histograms on submission predictions","metadata":{}},{"cell_type":"code","source":"%%time\ntry:\n    print(df_submission.shape)\n    plt.figure(figsize = (15,4))\n    plt.hist(df_submission['Prediction'].values, bins = 1000 )\n    plt.show()\n    print(df_submission.shape)\n    display(df_submission.describe())\n\n    for t in [0.1,0.2, 0.25,0.28, 0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]:\n        m = df_submission['Prediction'] > t\n        print(t, m.sum(), m.sum()/ (n_samples_predicted * n_targets_predicted ) )\n\n    print()    \n    try:\n        print( Y.sum(),  Y.sum()/ (Y.shape[0] * Y.shape[1]) )\n    except:\n        pass    \n\n    print('Here is fast rationale why we should think of threshold for F1 is around 0.28 - number of 1 in that case corresponds to train data')\nexcept:\n    pass\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exercises. Take home messages. Check yourself questions","metadata":{}},{"cell_type":"markdown","source":"## Exercises","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}