{"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":"Learning pytorch using CAFA5 as an example.\n\nHere we implement simple neural network (multi-layer perceptron) in pytorch, train it on CAFA5 data to submit predictions\n\nUsing already computed emebdding for the proteins sequences created by Andrey Shevtsov esm2-model: ( https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/406168 - Might use multiple embedding models in future work.\n\nThis model was prepared especially for use in this Cafa5 task and may be more relevant than the previously used T5 model\n\n**Ongoing task:** Change parameters for the neural network and try to improve the score.\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: https://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","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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-08-21T20:10:57.464332Z","iopub.execute_input":"2023-08-21T20:10:57.465041Z","iopub.status.idle":"2023-08-21T20:10:58.804705Z","shell.execute_reply.started":"2023-08-21T20:10:57.464992Z","shell.execute_reply":"2023-08-21T20:10:58.803365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture \n!pip install torchmetrics","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:10:58.807483Z","iopub.execute_input":"2023-08-21T20:10:58.808369Z","iopub.status.idle":"2023-08-21T20:11:11.082896Z","shell.execute_reply.started":"2023-08-21T20:10:58.808322Z","shell.execute_reply":"2023-08-21T20:11:11.081698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture\n!pip install torchsummary","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:11:11.084276Z","iopub.execute_input":"2023-08-21T20:11:11.084601Z","iopub.status.idle":"2023-08-21T20:11:22.218783Z","shell.execute_reply.started":"2023-08-21T20:11:11.084572Z","shell.execute_reply":"2023-08-21T20:11:22.217205Z"},"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","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:11:22.221715Z","iopub.execute_input":"2023-08-21T20:11:22.222063Z","iopub.status.idle":"2023-08-21T20:11:36.056988Z","shell.execute_reply.started":"2023-08-21T20:11:22.222031Z","shell.execute_reply":"2023-08-21T20:11:36.05591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Metrics**\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":"import 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 )","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:11:36.058591Z","iopub.execute_input":"2023-08-21T20:11:36.059241Z","iopub.status.idle":"2023-08-21T20:11:36.208967Z","shell.execute_reply.started":"2023-08-21T20:11:36.05921Z","shell.execute_reply":"2023-08-21T20:11:36.207846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**How multi-target cases are computed**\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, and for some metrics there is a parameter \"multidim_average\", taking values 'global' / 'samplewise' which control what way to choose.\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.","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-21T20:11:36.210212Z","iopub.execute_input":"2023-08-21T20:11:36.210555Z","iopub.status.idle":"2023-08-21T20:11:36.369277Z","shell.execute_reply.started":"2023-08-21T20:11:36.210526Z","shell.execute_reply":"2023-08-21T20:11:36.368266Z"},"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-21T20:11:36.370443Z","iopub.execute_input":"2023-08-21T20:11:36.370743Z","iopub.status.idle":"2023-08-21T20:11:36.387355Z","shell.execute_reply.started":"2023-08-21T20:11:36.370715Z","shell.execute_reply":"2023-08-21T20:11:36.386324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"^doesn't work- 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-21T20:11:36.388711Z","iopub.execute_input":"2023-08-21T20:11:36.389035Z","iopub.status.idle":"2023-08-21T20:11:36.396504Z","shell.execute_reply.started":"2023-08-21T20:11:36.388999Z","shell.execute_reply":"2023-08-21T20:11:36.395407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Loading precalulcated features - embeddings and targets transformed in 0,1 multi-target task.\n\nWe use already computed emebdding for the proteins sequences 'esm2-model' - thanks to Andrey Shevtsov \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. Note: original input file text describing targets have been trasnformed to np.array Y. See the notebook https://www.kaggle.com/code/alexandervc/baseline-multilabel-to-multitarget-binary for details.","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'","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:11:36.397993Z","iopub.execute_input":"2023-08-21T20:11:36.398847Z","iopub.status.idle":"2023-08-21T20:12:21.191991Z","shell.execute_reply.started":"2023-08-21T20:11:36.398817Z","shell.execute_reply":"2023-08-21T20:12:21.19087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Define neural network**\nThis is a simple pytourch NN\nIt will need to be optimised for Cafa5 in future versions","metadata":{}},{"cell_type":"code","source":"import torch.nn as nn","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:12:21.195485Z","iopub.execute_input":"2023-08-21T20:12:21.195822Z","iopub.status.idle":"2023-08-21T20:12:21.200802Z","shell.execute_reply.started":"2023-08-21T20:12:21.195792Z","shell.execute_reply":"2023-08-21T20:12:21.199784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class 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-21T20:12:21.201895Z","iopub.execute_input":"2023-08-21T20:12:21.202181Z","iopub.status.idle":"2023-08-21T20:12:21.216825Z","shell.execute_reply.started":"2023-08-21T20:12:21.202155Z","shell.execute_reply":"2023-08-21T20:12:21.21561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Choose device, CPU or GPU","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-21T20:12:21.218244Z","iopub.execute_input":"2023-08-21T20:12:21.218645Z","iopub.status.idle":"2023-08-21T20:12:21.232458Z","shell.execute_reply.started":"2023-08-21T20:12:21.218614Z","shell.execute_reply":"2023-08-21T20:12:21.231197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Model to deice","metadata":{}},{"cell_type":"code","source":"model = MyModel(X.shape[1],Y.shape[1])\nmodel.to(device)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:12:21.233779Z","iopub.execute_input":"2023-08-21T20:12:21.234092Z","iopub.status.idle":"2023-08-21T20:12:21.286163Z","shell.execute_reply.started":"2023-08-21T20:12:21.234067Z","shell.execute_reply":"2023-08-21T20:12:21.28511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"try:\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 subclasse","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:12:21.287353Z","iopub.execute_input":"2023-08-21T20:12:21.287684Z","iopub.status.idle":"2023-08-21T20:12:21.293411Z","shell.execute_reply.started":"2023-08-21T20:12:21.287656Z","shell.execute_reply":"2023-08-21T20:12:21.292439Z"},"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) ","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:12:21.294833Z","iopub.execute_input":"2023-08-21T20:12:21.295235Z","iopub.status.idle":"2023-08-21T20:12:21.310957Z","shell.execute_reply.started":"2023-08-21T20:12:21.295198Z","shell.execute_reply":"2023-08-21T20:12:21.309933Z"},"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 )","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:12:21.312282Z","iopub.execute_input":"2023-08-21T20:12:21.312638Z","iopub.status.idle":"2023-08-21T20:12:23.540953Z","shell.execute_reply.started":"2023-08-21T20:12:21.312607Z","shell.execute_reply":"2023-08-21T20:12:23.53981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Wrap X,Y into \"TensorDataset\" and create \"DataLoader\"","metadata":{}},{"cell_type":"markdown","source":"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-21T20:12:23.542546Z","iopub.execute_input":"2023-08-21T20:12:23.543201Z","iopub.status.idle":"2023-08-21T20:12:23.546713Z","shell.execute_reply.started":"2023-08-21T20:12:23.543164Z","shell.execute_reply":"2023-08-21T20:12:23.546004Z"},"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-21T20:12:23.547741Z","iopub.execute_input":"2023-08-21T20:12:23.5485Z","iopub.status.idle":"2023-08-21T20:12:23.565382Z","shell.execute_reply.started":"2023-08-21T20:12:23.548466Z","shell.execute_reply":"2023-08-21T20:12:23.563936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Training**","metadata":{}},{"cell_type":"markdown","source":"Training Pytorch NNs is relatively simple. Training implements the gradient descent method. Each step updates the weights of the NN to minimize loss function. Pytorch essentially is a blackbox which automatically calculates gradients and user should not care about that part. We just only need to do several steps:","metadata":{}},{"cell_type":"markdown","source":"model.train() - turns model into the standard mode when gradients will be computed - it is opposite to \"model.eval()\"\n\nmodel.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.\nwith torch.no_grad(): - similar ","metadata":{}},{"cell_type":"markdown","source":"1) We need to write loop over epochs, and the inner loop over batches. \n\n2) Technical step - optimizer.zero_grad() - set gradients to zero for each step \n\n3) Calcualte model predictions on each batch:   preds = model(x_batch) \n\n4) Calculate the for that batch:  loss = criterion(preds, y_batch) \n\n5) loss.backward() - calulation of the gradient (i.e. loss function with respect to weights of the neural network)\n\n6)  optimizer.step() - update the weights along the gradients","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-21T20:12:23.566788Z","iopub.execute_input":"2023-08-21T20:12:23.567193Z","iopub.status.idle":"2023-08-21T20:12:23.57999Z","shell.execute_reply.started":"2023-08-21T20:12:23.567154Z","shell.execute_reply":"2023-08-21T20:12:23.579105Z"},"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-21T20:12:23.581306Z","iopub.execute_input":"2023-08-21T20:12:23.581832Z","iopub.status.idle":"2023-08-21T20:19:51.985502Z","shell.execute_reply.started":"2023-08-21T20:12:23.581799Z","shell.execute_reply":"2023-08-21T20:19:51.984202Z"},"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() ","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:19:51.987621Z","iopub.execute_input":"2023-08-21T20:19:51.988115Z","iopub.status.idle":"2023-08-21T20:19:53.468708Z","shell.execute_reply.started":"2023-08-21T20:19:51.988073Z","shell.execute_reply":"2023-08-21T20:19:53.467472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The 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 yet.","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-21T20:19:53.470282Z","iopub.execute_input":"2023-08-21T20:19:53.471201Z","iopub.status.idle":"2023-08-21T20:19:53.816802Z","shell.execute_reply.started":"2023-08-21T20:19:53.471167Z","shell.execute_reply":"2023-08-21T20:19:53.815609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load features for submission**","metadata":{}},{"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])\n","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:19:53.818559Z","iopub.execute_input":"2023-08-21T20:19:53.818958Z","iopub.status.idle":"2023-08-21T20:20:19.645712Z","shell.execute_reply.started":"2023-08-21T20:19:53.818925Z","shell.execute_reply":"2023-08-21T20:20:19.644609Z"},"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-21T20:20:19.647652Z","iopub.execute_input":"2023-08-21T20:20:19.648062Z","iopub.status.idle":"2023-08-21T20:20:20.241622Z","shell.execute_reply.started":"2023-08-21T20:20:19.648023Z","shell.execute_reply":"2023-08-21T20:20:20.240472Z"},"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-21T20:20:20.243193Z","iopub.execute_input":"2023-08-21T20:20:20.244042Z","iopub.status.idle":"2023-08-21T20:20:25.59657Z","shell.execute_reply.started":"2023-08-21T20:20:20.244005Z","shell.execute_reply":"2023-08-21T20:20:25.59562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(preds.shape)\nprint(preds[:4,:3])","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:20:25.597999Z","iopub.execute_input":"2023-08-21T20:20:25.598344Z","iopub.status.idle":"2023-08-21T20:20:25.605332Z","shell.execute_reply.started":"2023-08-21T20:20:25.598313Z","shell.execute_reply":"2023-08-21T20:20:25.604304Z"},"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-21T20:20:25.606599Z","iopub.execute_input":"2023-08-21T20:20:25.606906Z","iopub.status.idle":"2023-08-21T20:20:25.679338Z","shell.execute_reply.started":"2023-08-21T20:20:25.606879Z","shell.execute_reply":"2023-08-21T20:20:25.678423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Prepare submission 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-21T20:20:25.686308Z","iopub.execute_input":"2023-08-21T20:20:25.686697Z","iopub.status.idle":"2023-08-21T20:21:44.124906Z","shell.execute_reply.started":"2023-08-21T20:20:25.686663Z","shell.execute_reply":"2023-08-21T20:21:44.12379Z"},"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')","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:21:44.126422Z","iopub.execute_input":"2023-08-21T20:21:44.126753Z","iopub.status.idle":"2023-08-21T20:27:09.08113Z","shell.execute_reply.started":"2023-08-21T20:21:44.126722Z","shell.execute_reply":"2023-08-21T20:27:09.08008Z"},"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":{"execution":{"iopub.status.busy":"2023-08-21T20:27:09.082597Z","iopub.execute_input":"2023-08-21T20:27:09.082902Z","iopub.status.idle":"2023-08-21T20:32:45.384776Z","shell.execute_reply.started":"2023-08-21T20:27:09.082876Z","shell.execute_reply":"2023-08-21T20:32:45.383659Z"},"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\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')","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:32:45.386285Z","iopub.execute_input":"2023-08-21T20:32:45.386768Z","iopub.status.idle":"2023-08-21T20:34:23.71233Z","shell.execute_reply.started":"2023-08-21T20:32:45.386736Z","shell.execute_reply":"2023-08-21T20:34:23.711465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"execution":{"iopub.status.busy":"2023-08-21T20:34:23.71458Z","iopub.execute_input":"2023-08-21T20:34:23.715202Z","iopub.status.idle":"2023-08-21T20:34:23.72021Z","shell.execute_reply.started":"2023-08-21T20:34:23.715171Z","shell.execute_reply":"2023-08-21T20:34:23.719264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}