{"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 + toda: 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 - realization of CAFA5 F1 (weighted) metric Kaggle\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\n\ntodo: add more \n\n#### Versions: \n \n 25 - option to concatenate embeddings/other-features is added \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","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' # 'quick_run' # \nif mode_loc == 'quick_run': # quick-run mode - for debug and testing \n # Crude downsampling - just for speed run/check - about 2-3 minutes for run \n n_samples_to_consider = 10000# 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 n_folds_to_process = 1 # we can only use 1 fold - for fast checks\nelse:\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 = 2200 # Up to 31466 but more than 3000-5000 may crash RAM \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# cfg1 = {'model': {'id':'Ridge' } }\n# list_main_config_model_feature_etc.append( cfg1 )\n# cfg2 = {'model': {'id':'pMLP_Andrey' , 'n_selfblend':10 } } # Pytorch MLP model\n# list_main_config_model_feature_etc.append( cfg2 )\n# cfg3 = {'model': {'n_selfblend':5 , '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':10, 'id':'KMLP_simple' , 'epochs':15 , 'batch_size':128 , 'verbose':0} } # Simple Keras MLP\n# list_main_config_model_feature_etc.append( cfg4 )\ncfg5 = {'model': {'n_selfblend':5, 'id':'KMLP', 'Layers': [1350,450],'Dropouts': [0.1 ], 'BatchNormalizations': [] , 'epochs':15 , 'batch_size':128 , 'verbose':0 } } \n# Keras MLP - with adjustable params \nlist_main_config_model_feature_etc.append( cfg5 )\n\n\n\n######################################### Feature sets ########################################################\n\n# Features: \n# default_features = [ {'id':'t5'}, {'id':'esm2S2560','preproc':'PCA100'} ] # 'esm2S2560' # 'esm2S1280' # 'esm2S320' # 't5'# \nlist_features_id = ['t5', 'esm2S1280', 'esm2S320' ] # 'esm2S640' # 'esm2S480' # 't5'# 'esm2S2560' # 'esm2S1280' # 'esm2S320' # \n# Concatenate several embeddings (other features)\n# !!! Pay attention proteins ids should be the same as in all the files !!!!!!!!!!!!!!!!!!\n\n######################################### Furthter params ########################################################\n\nflag_compute_oof_predictions = True # Necessary to compute local CV scores \n\nmode_submit = True # Compute prediction for submission part and prepare submission file in required CAFA5 format. Set to False if you only interested in local score\ncutoff_threshold_low = 0.01 # prediction < cutoff_threshold_low will be set to zero (i.e. no need to save to submission file)\n\nflag_compute_each_blend_stat = True \nflag_compute_cafa_f1_for_each_blend = False \n\nflag_compute_final_model_stat = True\n\nRANDOM_SEED = None # Fix or Not random seed \n\n\nstr_id = str(list_features_id)+'_Y'+str(n_labels_to_consider)\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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"logs_file_path = 'logs.txt'","metadata":{"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":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture\n!pip install torchsummary\nfrom torchsummary import summary as torchsummary\n","metadata":{"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","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":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load X,Y, etc \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 return fn_X, fn_protein_ids, fn_X_submit, fn_submit_protein_ids \n \n################ load features #########################\nprint(); print('!!! Pay attention proteins ids should be the same as in all the files !!!!!!!!!!!!!!!!!! '); print();\nfor 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 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 print(X.shape)\n print(X[:2,:3])\n prot_ids = np.load(fn_protein_ids)[:n_samples_to_consider ]\n vec_train_protein_ids = prot_ids\n print('prot_ids.shape:', prot_ids.shape)\n print('prot_ids[:15]:', prot_ids[:15])\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 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 print(X_submit.shape)\n print(X_submit[:2,:3])\n\n fn = fn_submit_protein_ids\n submit_protein_ids = np.load(fn)\n print(submit_protein_ids.shape, submit_protein_ids[:10])\n\n\n############################ load targets and their ids ######################################\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,:n_labels_to_consider].toarray()\nprint('Y', Y.shape, 'truncated')\nn_labels_to_consider = Y.shape[1] # in case n_labels_to_consider is greater that Y.shape we decrease it\n\nfn = '/kaggle/input/cafa5-features-etc/Y_31466_labels.npy'\nY_labels = np.load(fn, allow_pickle=True )[:n_labels_to_consider]\nlabels_to_consider = Y_labels\nprint(Y_labels.shape)\nprint(Y_labels[:20])\n\n# %%time\nif 1:\n v = Y.sum(axis = 0)\n plt.figure(figsize = (20,6))\n plt.plot(v[:500], '*-')\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","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"log_available_ram('After data load')","metadata":{"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":{"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":{"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":{"trusted":true},"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":{"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 \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 return model, str_model_id\n\nmodel_config_tmp = {'id':'Ridge'}\nmodel_config_tmp = {'id':'Ridge', 'alpha':10, }\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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fit model function ","metadata":{}},{"cell_type":"code","source":"def 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 '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', 10 )\n BATCH_SIZE = model_config.get('batch_size', 128 )\n LEARNING_RATE = model_config.get( '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'Right After Model, X_train initialization. Training loop fold: {ix_fold}') \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} ' ) ","metadata":{"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":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Auxilliary functions compute/save scores, etc ","metadata":{}},{"cell_type":"code","source":"%%time\nimport time \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 dict_scores_etc = get_F1_etc_scores_official_CAFA_evaluation( Y_pred[IX_loc,:], IX_loc, 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 \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 df_stat.to_csv('df_stat.csv') \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":{"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":{"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######################### 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 \ncnt_blend_submit = 0 ; \n\nif flag_compute_oof_predictions:\n Y_pred_oof_blend = np.zeros( ( Y.shape ) , dtype = np.float16 )\ncnt_blend_oof = -1;\n\n########################## Preparations ###########################################\nlog_available_ram('Right before modeling')\n\n\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 model_config = main_config_model_feature_etc['model']\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 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 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 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\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 \n ##################### Compute predictions for submission and blend with the previous one ###################################################\n if mode_submit : \n t0 = time.time()\n Y_pred_submit_current = model_predict(model , X_submit, model_config , str_model_id , verbose = 0 )\n Y_submit = (Y_submit * cnt_blend_submit + Y_pred_submit_current )/ (cnt_blend_submit + 1); # Average predictions from different folds/models\n cnt_blend_submit += 1 \n time_pred_submit = time.time() - t0\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 del X_train, Y_train, model , X_val, Y_pred_val\n torch.cuda.empty_cache()\n gc.collect()\n\n if flag_compute_each_blend_stat:\n if 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 = flag_compute_cafa_f1_for_each_blend , str_model_id = str_model_id, dict_optional_info = {'Time': time_one_model }, verbose = 0)\n log_available_ram('After Stat Calculation')\n \n if flag_compute_final_model_stat:\n if 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', dict_optional_info = { }, verbose = 0)\n log_available_ram('After Final Stat Calculation')\n\n \ndisplay(df_stat) ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_stat.shape)\ndisplay( df_stat )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(df_stat['AUC'].values,'*-')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ncol = 'AUC' \ncol2 = 'Fold'\nd2 = df_stat.groupby(col2).mean()\ndisplay(d2)\nd2.to_csv('df_stat_folds_mean.csv')\nplt.plot(d2[col].values,'*-')\nplt.xlabel(col2,fontsize = 20 )\nplt.title(col,fontsize = 20)\nplt.show()\n\ncol2 = 'Model'\nd3 = df_stat.groupby(col2).mean()\ndisplay(d3)\nd3.to_csv('df_stat_models_mean.csv')\nplt.plot(d3[col].values,'*-')\nplt.title(col,fontsize = 20)\nplt.xlabel(col2,fontsize = 20 )\nplt.show()\ndisplay()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = df_stat['Model'] == 'Final'\ndisplay(df_stat[m] )\ndisplay(df_stat[m].mean() )\nprint( df_stat[m].mean().loc['F1w'] )\ndf_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:\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":[]}]}