{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?\n\nLearn pytorch using CAFA5 as an example. Educational notebook. \n\nHere we implement simple neural network (multi-layer perceptron)  in pytorch, train it on CAFA5 data and submit predictions to the competition server.\n\nWe also add some educational comments.\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\n#### Educational task:\n\nChange parameters for the neural network and try to improve the score. \n\n\n#### Background:\n\n\nSimple example of NN on pytorch: https://www.kaggle.com/code/alexandervc/moa-nn-04/notebook\n\nThe most simple exercises to implement Linear and Logistic regression in pytorch: \nhttps://www.kaggle.com/code/alexandervc/pytorch-03-basics-implement-linear-regressions\nhttps://www.kaggle.com/code/alexandervc/pytorch-04-basics-implement-logregression\n\n\n#### Further steps: \n\n\nLater on try understand/modify/apply to CAFA5: \n\nhttps://www.kaggle.com/code/romanvinogradov/cafa-mlp-inference\n\nIdea by \"tmp\" - use cnn for such types of data, but make a first dense layer which arranges data to be suitable for cnn:  https://www.kaggle.com/competitions/lish-moa/discussion/202256\n\nRelated ideas: Senkin13: https://www.kaggle.com/code/senkin13/2nd-place-gru-cite , modifications by Dmitry Ershov: \nhttps://www.kaggle.com/code/bejeweled/rotate-2nd-place-cite-2d-cnn and other 3 notebooks by him (top-scored on Kaggle Open Problems (late submissions)): https://www.kaggle.com/competitions/open-problems-multimodal/code?competitionId=38128&sortBy=scoreDescending\n\n\n#### Related Notebooks and Results of Experiments :\n\nLB Score 0.45155 Version 9 here. Pay attention - if you rerun the same notebook and submit (change NOTHING) if you are lucky you get scores up to 0.47+ at least - because LB contains about 100 proteins and quite unstable - see e.g. 0.47509 https://www.kaggle.com/code/lidiashishina/pytorch-01-basics?scriptVersionId=138304155\n\nNote: LB falls to 0.39970 Version 16 if we increase number of epochs. In that case local BCE and AUCROC scores also decreases - so LB score fall down is not a complete surprise, however local F1|0.25 grows with epochs - that means F1|0.25 is not that much good as LB score estimator. Unfortunately. See version 16 here:  https://www.kaggle.com/code/alexandervc/pytorch-01-basics?scriptVersionId=138385096\n\nChange ReLU to GELU - LB: 0.47884: https://www.kaggle.com/code/lidiashishina/pytorch-01-basics?scriptVersionId=138309928\n\nChanging the network and LR scheduler - LB: 0.49087: \nhttps://www.kaggle.com/code/lidiashishina/pytorch-01-basics?scriptVersionId=138484544\n\nVarious losses - like focal, Tversky are explored here: https://www.kaggle.com/code/ilyasokolov/pytorch-01-basics\n\n\n","metadata":{}},{"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-01T08:50:14.036663Z","iopub.execute_input":"2023-08-01T08:50:14.037108Z","iopub.status.idle":"2023-08-01T08:50:15.750858Z","shell.execute_reply.started":"2023-08-01T08:50:14.037074Z","shell.execute_reply":"2023-08-01T08:50:15.749715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture \n!pip install torchmetrics","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:50:15.756137Z","iopub.execute_input":"2023-08-01T08:50:15.75674Z","iopub.status.idle":"2023-08-01T08:50:31.575348Z","shell.execute_reply.started":"2023-08-01T08:50:15.756691Z","shell.execute_reply":"2023-08-01T08:50:31.573289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture\n!pip install torchsummary","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:50:31.579901Z","iopub.execute_input":"2023-08-01T08:50:31.58032Z","iopub.status.idle":"2023-08-01T08:50:45.444872Z","shell.execute_reply.started":"2023-08-01T08:50:31.580282Z","shell.execute_reply":"2023-08-01T08:50:45.442278Z"},"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\nfrom torch.utils.data import Dataset, DataLoader\nfrom torchmetrics import AUROC,F1Score\nfrom torchmetrics.classification import BinaryF1Score\nfrom torchsummary import summary as torchsummary\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-08-01T08:50:45.450981Z","iopub.execute_input":"2023-08-01T08:50:45.451525Z","iopub.status.idle":"2023-08-01T08:51:05.450297Z","shell.execute_reply.started":"2023-08-01T08:50:45.451482Z","shell.execute_reply":"2023-08-01T08:51:05.448127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Metrics \n\nSome exercises and comments on metrics used in pytorch and torchmetrics. \n\n\"torchmetrics\" - is additional package to pytorch - contains many metrics which basic pytorch does not \n\nPay attention:  input arguments should be torch.tensor not lists or np.arrays \n\nPay attention:  auc function - order of arguments is opposite to sklearn rocauc function, and on some other details of intialization - different from the sklearn syntax. The syntax is auroc - not auCroc - as one may expect\n\nFunny note: same losses e.g. BCE can be imported from different modelues with differnt names, but they are the same - e.g. torch.nn.functional under the name binary_cross_entropy and import torch.nn as nn - as  nn.BCELoss()","metadata":{}},{"cell_type":"code","source":"\nimport torch.nn.functional as F\n\n# F.binary_cross_entropy(out, labels)   # Calculate loss\ns =  F.binary_cross_entropy( torch.tensor([0.,1,0.5]), torch.tensor([0,1.,1] ) )  # Important to have float - on [0,1] it crashes - \"Long not supported\"\nprint( s )\nprint( s.grad_fn )\n\nimport torch.nn as nn\npreds = torch.tensor([0.,1,0.5])\ntarget = torch.tensor([0,1.,1])\nloss_func = nn.BCELoss()\ns = loss_func(preds, target )\nprint( s )\nprint( s.grad_fn )\n\nimport torch.nn as nn\npreds = torch.tensor([0.,1,0.5])\ntarget = torch.tensor([0,1.,1])\nloss_func  = torch.nn.MSELoss() \ns = loss_func(preds, target )\nprint( s )\nprint( s.grad_fn )\n\n\nfrom torchmetrics import AUROC,F1Score\n\nauroc = AUROC(task = 'binary')\n# order of arguments is different from the sklearn : sklearn - y_true, y_score, while in pytorch - y_score, y_true\n# Also we first need to initialize  aucroc = AUROC(task = 'binary'), and only then apply  \ns = auroc(torch.tensor([0,.6,0]), torch.tensor([0,1,1]) ) \nprint('\\n AUC')\nprint( s )\nprint( s.grad_fn )\ns = auroc(torch.tensor([0,.01,0.02,0]), torch.tensor([0,1,1,1]) ) \nprint( s )\ns = auroc(torch.tensor([0,.1,0.2,0]), torch.tensor([0,1,1,1]) ) \nprint( s )\nprint()\n\ntarget = torch.tensor([0, 1, 1, 0, 1, 1])\npreds = torch.tensor([0, 1, 1, 0, 0, 1])\nf1 = F1Score(task=\"multiclass\", num_classes=2)\ns = f1(preds, target)\nprint( s )\nprint( s.grad_fn )\n\nfrom torchmetrics.classification import BinaryF1Score\nf1 = BinaryF1Score(threshold=0.5, multidim_average= 'global'  )\ns = f1(preds, target)\nprint( s )\nprint( s.grad_fn )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:05.45166Z","iopub.execute_input":"2023-08-01T08:51:05.452055Z","iopub.status.idle":"2023-08-01T08:51:05.701292Z","shell.execute_reply.started":"2023-08-01T08:51:05.451999Z","shell.execute_reply":"2023-08-01T08:51:05.699702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## How multi-target cases are computed \n\nMany  metrics support cases of many tartgets, i.e. Y_pred and Y_true are both matrices : n_samples x n_targets \n\nThere is NO such support in e.g. sklearn.metrics - where we typically work with single target tasks. \n\nThere are several ways how standard metrics can be extended to multi-target cases,\nand for some metrics there is a parameter \"multidim_average\", taking values 'global'  / 'samplewise' which control what way to choose.\n\n\nThe other detail - for the metrics like F1 predictions expected to be in [0,1] interval , if this fail, then torch and sklearn metrics may disagree. In standard case the agree as shown below. \n","metadata":{}},{"cell_type":"code","source":"target =  torch.tensor( (np.random.randn(100,10)>0).astype(float) )\npreds = torch.tensor( np.clip( np.random.randn(100,10) + 0.5, 0 ,1) )\n\n#  multidim_average= 'global'\nprint( \" multidim_average= 'global'   \", 'Scalar will be returned')\nf1 = BinaryF1Score(threshold=0.5, multidim_average= 'global'  )\ns = f1(preds, target)\nprint(s.shape)\nprint( s )\n\nprint()\nprint(\" multidim_average= 'samplewise'  \", 'vector of values for each sample will be returned')\nf1 = BinaryF1Score(threshold=0.5, multidim_average= 'samplewise'  )\ns = f1(preds, target)\nprint(s.shape)\nprint(s.mean() )\nprint( s )\n\n\nprint()\nprint()\nprint('Implement the same with sklearn')\nfrom sklearn.metrics import f1_score\nprint('Global case ')\ns = f1_score( target.numpy().ravel() ,   preds.numpy().ravel() >= 0.5  )\nprint(np.round(s,4))\n\nprint('Samplewise case:')\nl = [ f1_score( target.numpy()[i,:], preds.numpy()[i,:] >= 0.5 ) for i in range(preds.shape[0]) ]\ns = np.mean(l)\nprint( np.round(s,4))","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:05.703583Z","iopub.execute_input":"2023-08-01T08:51:05.704035Z","iopub.status.idle":"2023-08-01T08:51:05.938224Z","shell.execute_reply.started":"2023-08-01T08:51:05.703986Z","shell.execute_reply":"2023-08-01T08:51:05.936344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ns = auroc(preds, target ) \nprint( s )\n\n\nfrom sklearn.metrics import roc_auc_score\ns = roc_auc_score(target.numpy().ravel() , preds.numpy().ravel()) \nprint( np.round(s,4) )","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:05.94018Z","iopub.execute_input":"2023-08-01T08:51:05.940932Z","iopub.status.idle":"2023-08-01T08:51:05.964531Z","shell.execute_reply.started":"2023-08-01T08:51:05.940886Z","shell.execute_reply":"2023-08-01T08:51:05.963096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## do NOT do like that - canNOT apply metrics  to list/np.arrays etc  - we can  apply metrics to pytorch tensors only","metadata":{}},{"cell_type":"code","source":"try: \n    F.binary_cross_entropy( np.array([0.,1.,0.5]), np.array([0.,1.,1.] ) )\n    F.binary_cross_entropy( [0.,1.,0.5] ,  [0.,1.,1.]  )\n    aucroc(([0,.6,0]), ([0,1,1]) )\n    aucroc(np.array([0,.6,0]), np.array([0,1,1]) )\nexcept Exception as inst:\n    print(type(inst))    # the exception type\n    print(inst.args)     # arguments stored in .args\n    print(inst)      ","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:05.966579Z","iopub.execute_input":"2023-08-01T08:51:05.967969Z","iopub.status.idle":"2023-08-01T08:51:05.977706Z","shell.execute_reply.started":"2023-08-01T08:51:05.967916Z","shell.execute_reply":"2023-08-01T08:51:05.976352Z"},"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################ load  features #########################\n\nfn = '/kaggle/input/4637427/train_embeds_esm2_t36_3B_UR50D.npy'\nX = np.load(fn).astype(np.float32)\nprint(X.shape)\nprint(X[:2,:3])\nprot_ids  = np.load('/kaggle/input/4637427/train_ids_esm2_t36_3B_UR50D.npy')\nprint(prot_ids.shape)\nprint(prot_ids[:15])\n\n############################ load targets and their ids  ######################################\nfn = '/kaggle/input/cafa5-features-etc/Y_1499_labels.npy'\nY_labels = np.load(fn)\nprint(Y_labels.shape)\nprint(Y_labels[:20])\n# %%time \nfn = '/kaggle/input/cafa5-features-etc/Y_1499.npy'\nY = np.load(fn)\nprint(Y.shape)\nprint(Y[:20])\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# 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    ","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:05.979963Z","iopub.execute_input":"2023-08-01T08:51:05.980403Z","iopub.status.idle":"2023-08-01T08:51:47.962975Z","shell.execute_reply.started":"2023-08-01T08:51:05.980369Z","shell.execute_reply":"2023-08-01T08:51:47.961373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Define neural network \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-01T08:51:47.969055Z","iopub.execute_input":"2023-08-01T08:51:47.970057Z","iopub.status.idle":"2023-08-01T08:51:47.975952Z","shell.execute_reply.started":"2023-08-01T08:51:47.969988Z","shell.execute_reply":"2023-08-01T08:51:47.974097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/code/alexandervc/moa-nn-04/notebook\n\nclass MyModel(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, 800),\n            nn.ReLU(),\n\n            nn.BatchNorm1d(800),            # nn.Dropout(0.1),\n            nn.Linear(800, 600),\n            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(400, out_features),\n            nn.Sigmoid()\n        )\n\n    def forward(self, x): # \n        return self.model(x)","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:47.97817Z","iopub.execute_input":"2023-08-01T08:51:47.978658Z","iopub.status.idle":"2023-08-01T08:51:47.995982Z","shell.execute_reply.started":"2023-08-01T08:51:47.978615Z","shell.execute_reply":"2023-08-01T08:51:47.993898Z"},"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-01T08:51:47.998382Z","iopub.execute_input":"2023-08-01T08:51:47.998868Z","iopub.status.idle":"2023-08-01T08:51:48.032036Z","shell.execute_reply.started":"2023-08-01T08:51:47.998825Z","shell.execute_reply":"2023-08-01T08:51:48.029989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model to device","metadata":{}},{"cell_type":"code","source":"model = MyModel(X.shape[1],Y.shape[1])\nmodel.to(device)","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:48.034337Z","iopub.execute_input":"2023-08-01T08:51:48.034954Z","iopub.status.idle":"2023-08-01T08:51:48.097532Z","shell.execute_reply.started":"2023-08-01T08:51:48.034905Z","shell.execute_reply":"2023-08-01T08:51:48.095821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# That does not work for us - why ? \ntry:\n    torchsummary(model, XX.size(), batch_size=-1, device= 'cuda')\nexcept    Exception as inst:\n    print(type(inst))    # the exception type\n    print(inst.args)     # arguments stored in .args\n    print(inst)          # __str__ allows args to be printed directly,\n                         # but may be overridden in exception subclasses\n","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:48.100031Z","iopub.execute_input":"2023-08-01T08:51:48.101036Z","iopub.status.idle":"2023-08-01T08:51:48.110693Z","shell.execute_reply.started":"2023-08-01T08:51:48.100954Z","shell.execute_reply":"2023-08-01T08:51:48.108774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train test split ","metadata":{}},{"cell_type":"code","source":"%%time \nfrom sklearn.model_selection import train_test_split\n\nIX_train,IX_val = train_test_split( np.arange(len(X)), train_size=0.7, random_state=42)\nprint(IX_train.shape,IX_val.shape) \n","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:48.1126Z","iopub.execute_input":"2023-08-01T08:51:48.113077Z","iopub.status.idle":"2023-08-01T08:51:48.134659Z","shell.execute_reply.started":"2023-08-01T08:51:48.113034Z","shell.execute_reply":"2023-08-01T08:51:48.133571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create X/Y_train, X/Y_val as pytorch tensors and send to device","metadata":{}},{"cell_type":"code","source":"%%time\nX_train = torch.tensor(X[IX_train,:], dtype=torch.float32).to(device)\nY_train = torch.tensor(Y[IX_train,:], dtype=torch.float32).to(device)\n\nX_val = torch.tensor(X[IX_val,:], dtype=torch.float32).to(device)\nY_val = torch.tensor(Y[IX_val,:], dtype=torch.float32).to(device)\n\nprint(X_train.shape, Y_train.shape,X_val.shape, Y_val.shape )\n","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:48.136441Z","iopub.execute_input":"2023-08-01T08:51:48.136852Z","iopub.status.idle":"2023-08-01T08:51:51.12131Z","shell.execute_reply.started":"2023-08-01T08:51:48.13682Z","shell.execute_reply":"2023-08-01T08:51:51.119461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Wrap X,Y into \"TensorDataset\" and create \"DataLoader\"\n\n\"TensorDataset\" allows to unify several arrays to work with them as a whole. For example for pair X_train and Y_train we can create  \"train_dataset\"  which combines both. \n\nAs a next step we create \"DataLoader\" object from \"train_dataset\" which automatically take care about batches - see the main training loop of the neural network","metadata":{}},{"cell_type":"code","source":"# BATCH_SIZE = 5120\n# EPOCHS = 50\n# LEARNING_RATE = 0.001\n# MOMENTUM = 0.9\n# OPT_FUNC = torch.optim.Adam","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:51.123596Z","iopub.execute_input":"2023-08-01T08:51:51.124118Z","iopub.status.idle":"2023-08-01T08:51:51.130412Z","shell.execute_reply.started":"2023-08-01T08:51:51.124074Z","shell.execute_reply":"2023-08-01T08:51:51.129099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nBATCH_SIZE = 256\n\nfrom torch.utils.data import DataLoader, TensorDataset\ntrain_dataset = TensorDataset(X_train, Y_train)\ntrain_dataloader = DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True)\n\n# val_dataset = TensorDataset(X_val, Y_val)\n# val_dataloader = DataLoader(val_dataset, batch_size=BATCH_SIZE, shuffle=True)","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:51.132529Z","iopub.execute_input":"2023-08-01T08:51:51.133451Z","iopub.status.idle":"2023-08-01T08:51:51.150714Z","shell.execute_reply.started":"2023-08-01T08:51:51.133407Z","shell.execute_reply":"2023-08-01T08:51:51.14902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train neural network \n\nTraining neural networks in pytorch is NOT so simple as just write one command \"fit\" as in many other ML-libraries.\n\n\nBut also relatively simple. Training implements the gradient descent method. Each step updates the weights of the NN to minimize loss function. \nPytorch essentially is a blackbox which automatically calculates gradients and user should not care about that part. We just only need to do several steps:  \n\n    1) We need to write loop over epochs, and the inner loop over batches. \n\n    2) Technical step - optimizer.zero_grad() - set gradients to zero for each step \n    \n    3) Calcualte model predictions on each batch:   preds = model(x_batch) \n\n    4) Calculate the for that batch:  loss = criterion(preds, y_batch) \n    \n    5) loss.backward() - calulation of the gradient (i.e. loss function with respect to weights of the neural network)\n    \n    6)  optimizer.step() - update the weights along the gradients\n    \n    \nThat is basically all.\n\nMore comments: \n\n    model.train() - turns model into the standard mode when gradients will be computed - it is opposite to \"model.eval()\"\n    \n    model.eval() - turns model into the mode when gradients are NOT computed - that is done when we for example want NOT to train model, but just look what we get and evaluate metrics. So when we turn off graidents computate we do not waster resources - speed and RAM on their computation.\n    with torch.no_grad(): - similar \n    \n\n\n\n\n    ","metadata":{}},{"cell_type":"code","source":"import datetime\ncurrent_datetime = datetime.datetime.now()\n\nprint(\"Current datetime:\", current_datetime)","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:51.153156Z","iopub.execute_input":"2023-08-01T08:51:51.153795Z","iopub.status.idle":"2023-08-01T08:51:51.17049Z","shell.execute_reply.started":"2023-08-01T08:51:51.15374Z","shell.execute_reply":"2023-08-01T08:51:51.169116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nmax_epoch = 15\n\n\nimport time \nt0 = time.time()\n\nfrom torchmetrics.classification import BinaryF1Score\nf1 = BinaryF1Score(threshold=0.25, multidim_average= 'samplewise'  )\nfrom torchmetrics import AUROC\nauroc = AUROC(task = 'binary')\n\n\nprint('X.shape, Y.shape', X.shape, Y.shape )\nprint()\n\n# do not fortget to reinitialize model when you debug the notebook, since otherwise - each time you rerun the current cell \n# - you will update existing weights  - it would be easy to see that starting loss is unexpectdely small \nmodel = MyModel(X.shape[1],Y.shape[1])\nmodel.to(device)\n\n\n\noptimizer = torch.optim.Adam(model.parameters(), lr=1e-3) #SGD(model.parameters(), lr=1e-3)\n# scheduler = torch.optim.lr_scheduler.StepLR(optimizer, 0.05, 10)\ncriterion = nn.BCELoss()\n\n\ndf_stat = pd.DataFrame(); IX_stat = -1;\nprint()\nprint( str(model))\nprint()\ncurrent_datetime = datetime.datetime.now()\nprint('Start training NN',current_datetime)\nprint()\nfor epoch in range(max_epoch):\n    t0_epoch = time.time()\n    model.train() # switch model into train mode i.e. allow gradients computation - opposite of model.eval() - no grad computation\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        preds = model(x_batch)# Compute predictions only for batch samples \n\n        optimizer.zero_grad() # technical - set gradients to zero, otherwise they will be accumulated         \n\n        loss = criterion(preds, y_batch) # Compute loss function \n        loss.backward() # Compute gradients\n        optimizer.step() # Update NN weights using gradients\n\n        if i_batch % 100 == 0:\n            print(f'Epoch: {epoch}, batch: {i_batch},  train loss on batch: {loss.item():12.5f} , time: {time.time() - t0:.1f} ' )\n\n    t0_epoch_train = time.time() - t0_epoch            \n    model.eval()\n    with torch.no_grad():\n        train, y_train = train_dataset.tensors\n        # train, y_train = train.to(device), y_train.to(device) # They already should on device, why we need that ?\n        preds = model(train)\n        preds_train = preds\n        loss = criterion(preds, y_train)\n        loss1 = loss.item()\n        print()\n        print(f'Epoch: {epoch} finished,  train loss: {loss.item():12.5f} , time: {time.time() - t0:.1f} seconds ' )\n        \n        preds = model(X_val)\n        preds_val = preds\n        loss = criterion(preds, Y_val)\n        loss2 = loss.item()\n        print(f'Epoch: {epoch} finished,  VALidation loss: {loss.item():12.5f} , time: {time.time() - t0:.1f} seconds ' )\n        print()\n        \n        IX_stat += 1\n        df_stat.loc[IX_stat, 'epoch'] = epoch\n        df_stat.loc[IX_stat, 'BCE Val'] = loss2\n        df_stat.loc[IX_stat, 'BCE Train'] = loss1\n        \n        #from torchmetrics.classification import BinaryF1Score\n        f1 = BinaryF1Score(threshold=0.25, multidim_average= 'samplewise'  )\n        s = f1(preds_val, Y_val) # vector of values for each sample\n        df_stat.loc[IX_stat, 'F1|0.25 Val'] = s.mean().item()#  average over samples\n        s = f1(preds_train, y_train) \n        df_stat.loc[IX_stat, 'F1|0.25 Train'] =  s.mean().item()\n        \n        #from torchmetrics import AUROC\n        auroc = AUROC(task = 'binary')\n        s = auroc(preds_val, Y_val) # auroc between flattened arguments \n        df_stat.loc[IX_stat, 'AUC Val'] = s.item()#  average over samples\n        s = auroc(preds_train, y_train) # auroc between flattened arguments \n        df_stat.loc[IX_stat, 'AUC Train'] =  s.item()\n        \n        \n        for threshold_for_f1 in [0.2,0.3]:\n            #from torchmetrics.classification import BinaryF1Score\n            f1 = BinaryF1Score(threshold=threshold_for_f1, multidim_average= 'samplewise'  )\n            s = f1(preds_val, Y_val) # vector of values for each sample\n            df_stat.loc[IX_stat, 'F1|'+str(threshold_for_f1)+' Val'] = s.mean().item()#  average over samples\n            s = f1(preds_train, y_train) \n            df_stat.loc[IX_stat, 'F1|'+str(threshold_for_f1)+' Train'] =  s.mean().item()\n        \n        \n        df_stat.loc[IX_stat, 'Time epoch train'] = np.round(t0_epoch_train , 1 )\n        df_stat.loc[IX_stat, 'Time epoch full'] = np.round(time.time() - t0_epoch , 1 )\n        display(df_stat.tail(2))\n        \n\n        \nprint('Training finished', '%.1f seconds passed'%(time.time() - t0 ))\ndf_stat = df_stat.round(4)\ndisplay(df_stat)","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:51.173516Z","iopub.execute_input":"2023-08-01T08:51:51.174201Z","iopub.status.idle":"2023-08-01T08:51:51.632776Z","shell.execute_reply.started":"2023-08-01T08:51:51.174147Z","shell.execute_reply":"2023-08-01T08:51:51.631173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_metrics_keywords = ['BCE', 'AUC']\n\nplt.figure(figsize = (20,8) )\nfor i,kw in enumerate(list_metrics_keywords):\n    l = [col for col in df_stat if kw in col]\n    plt.subplot(1,len(list_metrics_keywords ),i+1)\n    plt.plot( df_stat[l], label = l )\n    plt.grid()\n    plt.legend(fontsize = 12)\n    plt.title(kw, fontsize = 20)\n    plt.xlabel('epoch',fontsize = 15)\n    \nplt.show() \n\nlist_metrics_keywords = ['F1|0.25','F1|0.2', 'F1|0.3' ]\n\nplt.figure(figsize = (20,8) )\nfor i,kw in enumerate(list_metrics_keywords):\n    l = [col for col in df_stat if kw in col]\n    plt.subplot(1,len(list_metrics_keywords ),i+1)\n    plt.plot( df_stat[l], label = l )\n    plt.grid()\n    plt.legend(fontsize = 12)\n    plt.title(kw, fontsize = 20)\n    plt.xlabel('epoch',fontsize = 15)\n    \nplt.show() \n","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:51.637414Z","iopub.execute_input":"2023-08-01T08:51:51.638146Z","iopub.status.idle":"2023-08-01T08:51:52.931939Z","shell.execute_reply.started":"2023-08-01T08:51:51.638086Z","shell.execute_reply":"2023-08-01T08:51:52.929981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Comments\n\nThe actual metric in the CAFA5 competition is similar to F1 above, but more sophisticated and computation time is signinificantly large. One may use F1 with threshold around 0.2-0.3 as a kind of \"fast-proxi\" for it.\n\nThe next step to get close to actual CAFA5 metric is to calculate F1 with different thresholds for different subontologies and choose not fixed threhold, but optimal one. Still it would not be exactly CAFA5 F1. Extra steps are: weighting and propagation. We will not go into these details here.","metadata":{}},{"cell_type":"markdown","source":"# Clear memory ","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()","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:52.934655Z","iopub.execute_input":"2023-08-01T08:51:52.936058Z","iopub.status.idle":"2023-08-01T08:51:53.471745Z","shell.execute_reply.started":"2023-08-01T08:51:52.935937Z","shell.execute_reply":"2023-08-01T08:51:53.469495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prepare submission","metadata":{}},{"cell_type":"markdown","source":"## Load features for submission","metadata":{"execution":{"iopub.status.busy":"2023-05-30T15:38:23.584425Z","iopub.execute_input":"2023-05-30T15:38:23.584794Z","iopub.status.idle":"2023-05-30T15:59:43.181211Z","shell.execute_reply.started":"2023-05-30T15:38:23.584742Z","shell.execute_reply":"2023-05-30T15:59:43.179867Z"}}},{"cell_type":"code","source":"%%time \n################ load  features #########################\n\n# fn = '/kaggle/input/4637427/train_embeds_esm2_t36_3B_UR50D.npy'\nfn = '/kaggle/input/4637427/test_embeds_esm2_t36_3B_UR50D.npy'\nprint(fn)\nX = np.load(fn).astype(np.float32)\nprint(X.shape)\nprint(X[:2,:3])","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:51:53.474653Z","iopub.execute_input":"2023-08-01T08:51:53.475449Z","iopub.status.idle":"2023-08-01T08:52:21.06254Z","shell.execute_reply.started":"2023-08-01T08:51:53.475318Z","shell.execute_reply":"2023-08-01T08:52:21.060728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nX = torch.tensor(X, dtype=torch.float32).to(device)","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:52:21.065064Z","iopub.execute_input":"2023-08-01T08:52:21.066069Z","iopub.status.idle":"2023-08-01T08:52:21.830078Z","shell.execute_reply.started":"2023-08-01T08:52:21.065973Z","shell.execute_reply":"2023-08-01T08:52:21.828348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Compute predictions","metadata":{}},{"cell_type":"code","source":"%%time\nmodel.eval()\nwith torch.no_grad():\n    preds = model(X)","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:52:21.831642Z","iopub.execute_input":"2023-08-01T08:52:21.832085Z","iopub.status.idle":"2023-08-01T08:52:41.277064Z","shell.execute_reply.started":"2023-08-01T08:52:21.832044Z","shell.execute_reply":"2023-08-01T08:52:41.275572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(preds.shape)\nprint(preds[:4,:3])","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:52:41.279662Z","iopub.execute_input":"2023-08-01T08:52:41.280381Z","iopub.status.idle":"2023-08-01T08:52:41.291664Z","shell.execute_reply.started":"2023-08-01T08:52:41.280339Z","shell.execute_reply":"2023-08-01T08:52:41.28961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ntest_protein_ids = np.load('/kaggle/input/4637427/test_ids_esm2_t36_3B_UR50D.npy')\nprint(test_protein_ids.shape, test_protein_ids[:10])","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:52:41.294501Z","iopub.execute_input":"2023-08-01T08:52:41.29531Z","iopub.status.idle":"2023-08-01T08:52:41.366908Z","shell.execute_reply.started":"2023-08-01T08:52:41.295253Z","shell.execute_reply":"2023-08-01T08:52:41.365917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare submission tsv file","metadata":{}},{"cell_type":"code","source":"%%time\ndf_submission = pd.DataFrame(columns = ['Protein Id', 'GO Term Id','Prediction'])\n\nn_targets_predicted = preds.shape[1]\nn_samples_predicted = preds.shape[0]\nprint('n_samples_predicted, n_targets_predicted',  n_samples_predicted, n_targets_predicted )\n\n\nprotein_list = []\nfor k in list(test_protein_ids):\n    protein_list += [k] * n_targets_predicted\ndf_submission['Protein Id'] = protein_list\n\ndf_submission['GO Term Id'] = list(Y_labels) * n_samples_predicted\ndf_submission['Prediction'] = preds.ravel()\n\ndf_submission = df_submission.round(3)\ndf_submission = df_submission[ df_submission['Prediction'] > 0.01  ]\n\nmemory_usage_per_column = df_submission.memory_usage(deep=True)\ntotal_memory_usage = memory_usage_per_column.sum()\nprint(\"\\nTotal memory usage:\", total_memory_usage/1e6, \"Megabytes\")\n\nprint(df_submission.shape)\ndisplay(df_submission)","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:52:41.376046Z","iopub.execute_input":"2023-08-01T08:52:41.378062Z","iopub.status.idle":"2023-08-01T08:56:55.124602Z","shell.execute_reply.started":"2023-08-01T08:52:41.377994Z","shell.execute_reply":"2023-08-01T08:56:55.122752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport gc\nif 0:\n    del preds \n    \ngc.collect()\n\ndf_submission.to_csv(\"submission.tsv\",header=False, index=False,sep='\\t')\n","metadata":{"execution":{"iopub.status.busy":"2023-08-01T08:56:55.12643Z","iopub.execute_input":"2023-08-01T08:56:55.126975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(df_submission.shape)\nplt.figure(figsize = (15,4))\nplt.hist(df_submission['Prediction'].values, bins = 1000 )\nplt.show()\nprint(df_submission.shape)\ndisplay(df_submission.describe())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"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    \nprint()    \ntry:\n    print( Y.sum(),  Y.sum()/ (Y.shape[0] * Y.shape[1]) )\nexcept:\n    pass    \n\nprint('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')","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":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}