scripts/contact_map_filter.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "91f60c89-690c-47b5-8af0-e1ba34cae9c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import List, Tuple, Optional, Dict, NamedTuple, Union, Callable\n",
    "import numpy as np\n",
    "import torch\n",
    "from scipy.spatial.distance import squareform, pdist, cdist\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "import biotite.structure as bs\n",
    "from biotite.structure.io.pdbx import PDBxFile, get_structure\n",
    "from biotite.database import rcsb\n",
    "from tqdm import tqdm\n",
    "import pandas as pd\n",
    "import torch.nn.functional as F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4156b2ee-ff33-478f-b7da-f52a356ce860",
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.optim import AdamW\n",
    "from transformers import get_scheduler\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "from transformers import EsmTokenizer, EsmForSequenceClassification, EsmModel, EsmConfig,EsmForMaskedLM\n",
    "from sklearn.metrics import accuracy_score\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "32496cc7-faac-4a9e-a832-a70a092510d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run it in cuda or cpu\n",
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "34a4296f-9ffd-4eaa-9490-c705707a0714",
   "metadata": {},
   "outputs": [],
   "source": [
    "def symmetrize(x):\n",
    "    \"Make layer symmetric in final two dimensions, used for contact prediction.\"\n",
    "    return x + x.transpose(-1, -2)\n",
    "\n",
    "def apc(x):\n",
    "    \"Perform average product correct, used for contact prediction.\"\n",
    "    a1 = x.sum(-1, keepdims=True)\n",
    "    a2 = x.sum(-2, keepdims=True)\n",
    "    a12 = x.sum((-1, -2), keepdims=True)\n",
    "\n",
    "    avg = a1 * a2\n",
    "    avg.div_(a12)  # in-place to reduce memory\n",
    "    normalized = x - avg\n",
    "    return normalized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "40a9a20d-d8fa-4e23-b209-1639dd31c265",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Adapted from https://github.com/facebookresearch/esm/blob/main/esm/modules.py\n",
    "class AttentionLogisticRegression(nn.Module):\n",
    "    \"\"\"Performs symmetrization, apc, and computes a logistic regression on the output features\"\"\"\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        in_features:int,\n",
    "        prepend_bos: bool,\n",
    "        append_eos: bool,\n",
    "        bias=True,\n",
    "        eos_idx: Optional[int] = None,\n",
    "    ):\n",
    "        super().__init__()\n",
    "        self.in_features=in_features\n",
    "        self.prepend_bos = prepend_bos\n",
    "        self.append_eos = append_eos\n",
    "        if append_eos and eos_idx is None:\n",
    "            raise ValueError(\"Using an alphabet with eos token, but no eos token was passed in.\")\n",
    "        self.eos_idx = eos_idx\n",
    "        self.regression = nn.Linear(in_features, 1, bias)\n",
    "        self.activation = nn.Sigmoid()\n",
    "    \n",
    "    def forward(self, tokens, attentions):\n",
    "        # remove eos token attentions\n",
    "        if self.append_eos:\n",
    "            eos_mask = tokens.ne(self.eos_idx).to(attentions)\n",
    "            eos_mask = eos_mask.unsqueeze(1) * eos_mask.unsqueeze(2)\n",
    "            attentions = attentions * eos_mask[:, None, None, :, :]\n",
    "            attentions = attentions[..., :-1, :-1]\n",
    "        # remove cls token attentions\n",
    "        if self.prepend_bos:\n",
    "            attentions = attentions[..., 1:, 1:]\n",
    "        batch_size, layers, heads, seqlen, _ = attentions.size()\n",
    "        attentions = attentions.view(batch_size, layers * heads, seqlen, seqlen)\n",
    "        attentions = attentions.to(self.regression.weight.device)  # attentions always float32, may need to convert to float16\n",
    "        attentions= apc(symmetrize(attentions))\n",
    "        attentions = attentions.permute(0, 2, 3, 1)\n",
    "        \n",
    "        return self.activation(self.regression(attentions).squeeze(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b86953e8-b30c-4f76-aa3e-9ad2ac43c5f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_fasta_dict(fasta_file):\n",
    "    fasta_dict = {}\n",
    "    with open(fasta_file, 'r') as infile:\n",
    "        for line in infile:\n",
    "            if line.startswith(\">\"):\n",
    "                head = line.replace(\"\\n\", \"\").replace(\">\", \"\")\n",
    "                fasta_dict[head] = ''\n",
    "            else:\n",
    "                fasta_dict[head] += line.replace(\"\\n\", \"\")\n",
    "    return fasta_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ac38cfbf-4941-4a87-8a47-363eebc5f5dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pred_contact(attention):\n",
    "    attention=torch.where(attention < 0.9, torch.tensor(0), torch.tensor(1))\n",
    "    \n",
    "    return attention"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "eca8ccbc-7dc4-4462-99eb-5d2cf2edc88d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AttentionLogisticRegression(\n",
       "  (regression): Linear(in_features=660, out_features=1, bias=True)\n",
       "  (activation): Sigmoid()\n",
       ")"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mymodel=torch.load('../models/contact-based model.pt')\n",
    "mymodel.eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5de583bd-f0d5-48ed-b93d-ad5cb55f786e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "model_name = \"../models/esm2_650M\"\n",
    "num_classes=2\n",
    "max_length=24\n",
    "tokenizer = EsmTokenizer.from_pretrained(model_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a2df473a-baf7-499c-bbb9-3f127d9b5826",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Some weights of the model checkpoint at ../models/esm2_650M were not used when initializing EsmForSequenceClassification: ['lm_head.dense.weight', 'lm_head.layer_norm.weight', 'lm_head.dense.bias', 'lm_head.layer_norm.bias', 'lm_head.decoder.weight', 'lm_head.bias']\n",
      "- This IS expected if you are initializing EsmForSequenceClassification from the checkpoint of a model trained on another task or with another architecture (e.g. initializing a BertForSequenceClassification model from a BertForPreTraining model).\n",
      "- This IS NOT expected if you are initializing EsmForSequenceClassification from the checkpoint of a model that you expect to be exactly identical (initializing a BertForSequenceClassification model from a BertForSequenceClassification model).\n",
      "Some weights of EsmForSequenceClassification were not initialized from the model checkpoint at ../models/esm2_650M and are newly initialized: ['classifier.out_proj.weight', 'classifier.dense.weight', 'classifier.dense.bias', 'classifier.out_proj.bias']\n",
      "You should probably TRAIN this model on a down-stream task to be able to use it for predictions and inference.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "PeftModelForSequenceClassification(\n",
       "  (base_model): LoraModel(\n",
       "    (model): EsmForSequenceClassification(\n",
       "      (esm): EsmModel(\n",
       "        (embeddings): EsmEmbeddings(\n",
       "          (word_embeddings): Embedding(33, 1280, padding_idx=1)\n",
       "          (dropout): Dropout(p=0.0, inplace=False)\n",
       "          (position_embeddings): Embedding(1026, 1280, padding_idx=1)\n",
       "        )\n",
       "        (encoder): EsmEncoder(\n",
       "          (layer): ModuleList(\n",
       "            (0-32): 33 x EsmLayer(\n",
       "              (attention): EsmAttention(\n",
       "                (self): EsmSelfAttention(\n",
       "                  (query): Linear(\n",
       "                    in_features=1280, out_features=1280, bias=True\n",
       "                    (lora_dropout): ModuleDict(\n",
       "                      (default): Dropout(p=0.6, inplace=False)\n",
       "                    )\n",
       "                    (lora_A): ModuleDict(\n",
       "                      (default): Linear(in_features=1280, out_features=48, bias=False)\n",
       "                    )\n",
       "                    (lora_B): ModuleDict(\n",
       "                      (default): Linear(in_features=48, out_features=1280, bias=False)\n",
       "                    )\n",
       "                    (lora_embedding_A): ParameterDict()\n",
       "                    (lora_embedding_B): ParameterDict()\n",
       "                  )\n",
       "                  (key): Linear(\n",
       "                    in_features=1280, out_features=1280, bias=True\n",
       "                    (lora_dropout): ModuleDict(\n",
       "                      (default): Dropout(p=0.6, inplace=False)\n",
       "                    )\n",
       "                    (lora_A): ModuleDict(\n",
       "                      (default): Linear(in_features=1280, out_features=48, bias=False)\n",
       "                    )\n",
       "                    (lora_B): ModuleDict(\n",
       "                      (default): Linear(in_features=48, out_features=1280, bias=False)\n",
       "                    )\n",
       "                    (lora_embedding_A): ParameterDict()\n",
       "                    (lora_embedding_B): ParameterDict()\n",
       "                  )\n",
       "                  (value): Linear(\n",
       "                    in_features=1280, out_features=1280, bias=True\n",
       "                    (lora_dropout): ModuleDict(\n",
       "                      (default): Dropout(p=0.6, inplace=False)\n",
       "                    )\n",
       "                    (lora_A): ModuleDict(\n",
       "                      (default): Linear(in_features=1280, out_features=48, bias=False)\n",
       "                    )\n",
       "                    (lora_B): ModuleDict(\n",
       "                      (default): Linear(in_features=48, out_features=1280, bias=False)\n",
       "                    )\n",
       "                    (lora_embedding_A): ParameterDict()\n",
       "                    (lora_embedding_B): ParameterDict()\n",
       "                  )\n",
       "                  (dropout): Dropout(p=0.0, inplace=False)\n",
       "                  (rotary_embeddings): RotaryEmbedding()\n",
       "                )\n",
       "                (output): EsmSelfOutput(\n",
       "                  (dense): Linear(in_features=1280, out_features=1280, bias=True)\n",
       "                  (dropout): Dropout(p=0.0, inplace=False)\n",
       "                )\n",
       "                (LayerNorm): LayerNorm((1280,), eps=1e-05, elementwise_affine=True)\n",
       "              )\n",
       "              (intermediate): EsmIntermediate(\n",
       "                (dense): Linear(in_features=1280, out_features=5120, bias=True)\n",
       "              )\n",
       "              (output): EsmOutput(\n",
       "                (dense): Linear(in_features=5120, out_features=1280, bias=True)\n",
       "                (dropout): Dropout(p=0.0, inplace=False)\n",
       "              )\n",
       "              (LayerNorm): LayerNorm((1280,), eps=1e-05, elementwise_affine=True)\n",
       "            )\n",
       "          )\n",
       "          (emb_layer_norm_after): LayerNorm((1280,), eps=1e-05, elementwise_affine=True)\n",
       "        )\n",
       "        (contact_head): EsmContactPredictionHead(\n",
       "          (regression): Linear(in_features=660, out_features=1, bias=True)\n",
       "          (activation): Sigmoid()\n",
       "        )\n",
       "      )\n",
       "      (classifier): ModulesToSaveWrapper(\n",
       "        (original_module): EsmClassificationHead(\n",
       "          (dense): Linear(in_features=1280, out_features=1280, bias=True)\n",
       "          (dropout): Dropout(p=0.0, inplace=False)\n",
       "          (out_proj): Linear(in_features=1280, out_features=2, bias=True)\n",
       "        )\n",
       "        (modules_to_save): ModuleDict(\n",
       "          (default): EsmClassificationHead(\n",
       "            (dense): Linear(in_features=1280, out_features=1280, bias=True)\n",
       "            (dropout): Dropout(p=0.0, inplace=False)\n",
       "            (out_proj): Linear(in_features=1280, out_features=2, bias=True)\n",
       "          )\n",
       "        )\n",
       "      )\n",
       "    )\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from peft import PeftModel, PeftConfig\n",
    "config = PeftConfig.from_pretrained('../models/esm2_650M_LORA_SEQ_CLS_0.99')\n",
    "model = EsmForSequenceClassification.from_pretrained('../models/esm2_650M', num_labels=2)\n",
    " \n",
    "model = PeftModel.from_pretrained(model, '../models/esm2_650M_LORA_SEQ_CLS_0.99')\n",
    "model.eval()\n",
    "model.to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "37d2d1d8-4821-49c5-b64c-20b76945e386",
   "metadata": {},
   "outputs": [],
   "source": [
    "lbdin_dict=get_fasta_dict('../database/LBD.fasta')\n",
    "seqs=[]\n",
    "for header,seq in lbdin_dict.items():\n",
    "    seqs.append(seq)\n",
    "\n",
    "tokens = tokenizer(seqs, return_tensors='pt', padding=\"max_length\", truncation=True, max_length=max_length).to(\"cuda\")\n",
    "with torch.no_grad():\n",
    "    outputs = model(**tokens)\n",
    "    logits = outputs.logits\n",
    "    predictions = torch.argmax(logits, dim=1)\n",
    "    predictions=predictions.cpu()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c7c8ccf6-fc73-44a3-a8d3-eaccc304214b",
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs = tokenizer(seqs, return_tensors='pt', padding=\"max_length\", truncation=True, max_length=max_length).to(device)\n",
    "outputs=model.esm(**inputs,output_attentions=True,output_hidden_states=True)\n",
    "lbdin_lora_attention=torch.stack(outputs.attentions,1)\n",
    "lbdin_outputs=mymodel(inputs['input_ids'],lbdin_lora_attention).cpu()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "66f5647e-db2d-4a00-b43e-9d0987b305b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "index=torch.where(predictions==0)\n",
    "lbdin_filter=lbdin_outputs[index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3bec976e-b5d5-4b64-8201-b54d1170f1c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "lbdsy_dict=get_fasta_dict('../database/LBD_test.fasta')\n",
    "seqs=[]\n",
    "for header,seq in lbdsy_dict.items():\n",
    "    seqs.append(seq)\n",
    "\n",
    "inputs = tokenizer(seqs, return_tensors='pt', padding=\"max_length\", truncation=True, max_length=max_length).to(device)\n",
    "outputs=model.esm(**inputs,output_attentions=True,output_hidden_states=True)\n",
    "lbdsy_lora_attention=torch.stack(outputs.attentions,1)\n",
    "lbdsy_outputs=mymodel(inputs['input_ids'],lbdsy_lora_attention).cpu()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aafd5a27-398b-4c3c-8e09-4ca98594fa81",
   "metadata": {},
   "source": [
    "## Plan A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ec3c9347-c135-4688-ac2d-acadc3e2308b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def binary_matrix_intersection(matrices):\n",
    "    # 将第一个矩阵作为基准\n",
    "    intersection_matrix = np.array(matrices[0])\n",
    "\n",
    "    # 逐个矩阵进行逻辑与运算\n",
    "    for matrix in matrices[1:]:\n",
    "        intersection_matrix = np.logical_and(intersection_matrix, matrix)\n",
    "\n",
    "    # 将布尔型矩阵转换为0和1组成的矩阵\n",
    "    intersection_matrix = intersection_matrix.astype(int)\n",
    "\n",
    "    return intersection_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7f758166-eb2c-4b97-b0c3-13f172eda38e",
   "metadata": {},
   "outputs": [],
   "source": [
    "lbdin_pred_contact=get_pred_contact(lbdin_filter)\n",
    "lbdsy_pred_contact=get_pred_contact(lbdsy_outputs)\n",
    "lbd_p=torch.cat([lbdin_pred_contact[:8],lbdsy_pred_contact[:6]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "8419cd57-f37a-4680-ba98-f58cffa658da",
   "metadata": {},
   "outputs": [],
   "source": [
    "ac_map=binary_matrix_intersection(np.array(lbd_p))\n",
    "ne_map=binary_matrix_intersection(np.array(lbdin_pred_contact[8:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "cbc7761d-2676-4b80-89d6-ad07b4ccc955",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 8,  9, 13, 15, 16, 16, 18, 19, 20, 21]),\n",
       " array([ 8,  9, 16, 15, 13, 18, 16, 19, 21, 20]))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.where(ac_map != ne_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e914aec9-47b3-4184-8bb8-f5a23ec76a17",
   "metadata": {},
   "source": [
    "## Plan B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d4c937e8-7bf7-4463-8549-d2d542cb29c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "class ContactMapRegression(nn.Module):\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        in_features:int,\n",
    "        bias=True,\n",
    "    ):\n",
    "        super().__init__()\n",
    "        self.in_features=in_features\n",
    "        self.regression = nn.Linear(in_features, 2, bias)\n",
    "        self.activation = nn.Sigmoid()\n",
    "    \n",
    "    def forward(self,contact_map):\n",
    "        contact_map=contact_map.reshape((contact_map.shape[0],contact_map.shape[1]*contact_map.shape[1]))\n",
    "        outputs=self.activation(self.regression(contact_map))\n",
    "        return outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "bdb7d10e-4931-47a3-8672-1fe8f8f3b054",
   "metadata": {},
   "outputs": [],
   "source": [
    "contactmodel=ContactMapRegression(22*22)\n",
    "contactmodel.to(device)\n",
    "\n",
    "# Define loss function and optimizer\n",
    "criterion = nn.CrossEntropyLoss()\n",
    "optimizer = optim.Adam(contactmodel.parameters(), lr=1e-2)\n",
    "\n",
    "# Training loop\n",
    "num_epochs = 500  # Adjust the number of epochs as needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "64c537c7-c468-49bd-829b-0bc712681a09",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_158084/312588753.py:4: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  cmin=torch.tensor(cmin,dtype=torch.float32).to(device)\n"
     ]
    }
   ],
   "source": [
    "lbd_p=torch.cat([lbdin_pred_contact[:8],lbdsy_pred_contact[:6]])\n",
    "lbd_n=lbdin_pred_contact[8:]\n",
    "cmin=torch.cat([lbd_p,lbd_n,lbdsy_pred_contact[6:]])\n",
    "cmin=torch.tensor(cmin,dtype=torch.float32).to(device)\n",
    "cmin=cmin.detach()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "c258e12f-b7d5-437d-a564-6f3545120fcb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_158084/4231155913.py:5: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  y=torch.tensor(y,dtype=torch.long).to(device)\n"
     ]
    }
   ],
   "source": [
    "yp=torch.ones(len(lbd_p))\n",
    "yn=torch.zeros(len(lbd_n))\n",
    "y_=torch.zeros(len(lbdsy_pred_contact[6:]))\n",
    "y=torch.cat([yp,yn,y_])\n",
    "y=torch.tensor(y,dtype=torch.long).to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "08766497-5de1-4f53-a386-5a504cb61fa6",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 500/500 [00:00<00:00, 1638.47it/s]\n"
     ]
    }
   ],
   "source": [
    "loss_list=[]\n",
    "accurate_list=[]\n",
    "for epoch in tqdm(range(num_epochs)):\n",
    "    contactmodel.train()\n",
    "    outputs=contactmodel(cmin)\n",
    "    \n",
    "    loss = criterion(outputs, y)\n",
    "    _, prediction = outputs.max(dim=1)\n",
    "    accurate = (prediction==y).sum().item()\n",
    "        \n",
    "    loss.backward()# calulate loss\n",
    "    optimizer.step()# update gradient\n",
    "    optimizer.zero_grad()# reset gradient\n",
    "        \n",
    "    train_accuracy = accurate / len(outputs)\n",
    "    average_loss = loss / len(outputs)\n",
    "        \n",
    "    loss_list.append(average_loss)\n",
    "    accurate_list.append(train_accuracy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "8bf8e639-df47-4427-9afa-5b6b52ce257f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc492d2fb20>]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAT8klEQVR4nO3de5Cdd13H8fc3m3t6SZPdltCkpEAciNxad0oVHItcJu04rQ4q7cCoDJo/tA6OjFjUqVpHx8so4lgvVRFBpRYVjBCmIBQvDIVspYWmtbCEYBLBbDZpa5Imezlf/zjPbk42u9mT7dmcPL/zfs2c2ef5Pb9zzveXbj/55blGZiJJqr8l3S5AktQZBrokFcJAl6RCGOiSVAgDXZIKsbRbX9zf35+bN2/u1tdLUi099NBDhzJzYLZtXQv0zZs3MzQ01K2vl6RaiohvzLXNXS6SVAgDXZIKYaBLUiEMdEkqhIEuSYWYN9Aj4r0RcTAiHp1je0TEH0TEcER8KSKu7XyZkqT5tDNDfx+w7SzbbwS2VK/twB8/+7IkSedq3vPQM/PfImLzWbrcArw/m/fhfTAi1kbEhsz8ZqeKrJOJyQZ/+dm9/N+J8W6XIukC9doXX8HLN63t+Od24sKiK4F9Lev7q7YzAj0ittOcxXPVVVd14KsvPI/sf4pf3/k4ABFdLkbSBenyS1ZesIHetsy8B7gHYHBwsMgnaxw6ehKAf7791bx046VdrkZSL+nEWS4HgE0t6xurtp50+NgYAOsuWt7lSiT1mk4E+g7gR6qzXa4HnurV/edwKtDXrzHQJZ1f8+5yiYgPAjcA/RGxH/hlYBlAZv4JsBO4CRgGjgNvXaxi62D06Birl/excllft0uR1GPaOcvltnm2J/BTHauo5g4fO8k6Z+eSuqBrt8+tg6eeGefk+OQ5vedbT59wd4ukrjDQ5zB88ChvePe/0ljAuTive/EVnS9IkuZhoM/h64eO0Ui4/TUvZMPalef03u96Qf8iVSVJczPQ53D4WPN88luv28TGy1Z3uRpJmp93W5zD6PTphyu6XIkktcdAn8Po0TFWLetj1XJPP5RUDwb6HA4fG/P0Q0m1YqDPYfTYGOu9fF9SjfTMQdHhg//H737iK4xPtnce4iP7nuSaq9YublGS1EE9E+ifeOx/+fij3+LFGy6hnbvaXrl2FTe9ZMOi1yVJndIzgT51j5WPv/27u12KJC2KntmHfvjYGJetdp+4pHL1TKB7kFNS6Xom0L0LoqTS9U6gH/W8ckllK/qg6InxSXbtPcxkIzl0bMzb2koqWtGB/v7P7eU3dv7X9Pqmdd5kS1K5ig70/3nyBGuW9/GBH38lS5cEWzdc0u2SJGnRFB3oo8fG6L94BddedVm3S5GkRVf0QVHPbJHUS4oO9NGjHgiV1DuKDnRvgSuplxQb6JnJkeNjrPOJQ5J6RLGB/vSJCcYnk34v95fUI4oN9MPVM0Hd5SKpVxQc6CcBA11S7yg20A8dbc7Q17sPXVKPKDbQp3e5uA9dUo8oPtA9D11Sryg20KceObdyWV+3S5Gk86K4QN/55W9y5NgYf/3gNzwgKqmntHVzrojYBrwH6AP+PDN/c8b2q4C/AtZWfe7IzJ2dLXV+D33jCD/5N/9J/0UrGJts8JxLVp7vEiSpa+adoUdEH3A3cCOwFbgtIrbO6PZLwH2ZeQ1wK/BHnS60HU8/Mw7AoaPNUxb/6M3XdqMMSeqKdna5XAcMZ+aezBwD7gVumdEngambjV8K/E/nSmzfM+OTp61f5i4XST2knUC/EtjXsr6/amv1K8BbImI/sBP46dk+KCK2R8RQRAyNjIwsoNyzG63ObAG4dNUylvUVd4hAkubUqcS7DXhfZm4EbgI+EBFnfHZm3pOZg5k5ODAw0KGvPuXw0VOB7umKknpNO4F+ANjUsr6xamv1NuA+gMz8HLAS6O9Egedi6nJ/8JJ/Sb2nnUDfBWyJiKsjYjnNg547ZvT5b+C1ABHxYpqB3vl9KvNo3eWy3itEJfWYeQM9MyeA24H7gcdpns2yOyLuioibq27vAH4iIh4BPgj8WGbmYhXd6uuHjrH5jo/xH189NH11KMDAxd7DRVJvaes89Oqc8p0z2u5sWX4MeFVnS2vPrr2HAfjIwwc4fGyMl29ay/e/4rm8fusV3ShHkrqmrUC/oFX/Dshs7nJ57Ysu562vurq7NUlSF9T+vL7jYxMAJMkRnyEqqYfVPtCn9ps/eXyciUYa6JJ6Vu0DferMlj0jRwHPbpHUu2of6FMz9L2jxwFY5xOKJPWo2h4U/dzXRvnQ0D4e3DN6WrtXiErqVbUN9A8N7eMfv3jqgtXn969h7eplXN2/potVSVL31DbQJxqnrlt68yuv4td/4KVdrEaSuq+2+9AnWy5EdTeLJNU50CdPBbqnKkpSnQO9ZYa+7iLPbJGk+gZ6w10uktSqiEC/dNWyLlYiSReGIgK9310uklTvQP+O513Gv7/zNTzn0pXdLkeSuq7Wgb68bwmb1q3udimSdEGobaBPNBos7YtulyFJF4zaBvpkwpIw0CVpSn0DvdFg6RIDXZKm1DjQYYmBLknTahzoztAlqVWNAz2doUtSi1oHujN0STqlvoGeSZ9nuUjStPoG+mTS5wxdkqbVN9DTQJekVvUN9IaBLkmtDHRJKkRtA33CQJek09Q20BsNz3KRpFa1DfSJRtLn3RYlaVpbgR4R2yLiiYgYjog75ujzwxHxWETsjoi/7WyZZ2p4HroknWbpfB0iog+4G3g9sB/YFRE7MvOxlj5bgHcBr8rMIxFx+WIVPGXCK0Ul6TTtzNCvA4Yzc09mjgH3ArfM6PMTwN2ZeQQgMw92tszTNRpJpndblKRW7QT6lcC+lvX9VVurbwO+LSI+GxEPRsS22T4oIrZHxFBEDI2MjCysYpoXFQHO0CWpRacOii4FtgA3ALcBfxYRa2d2ysx7MnMwMwcHBgYW/GWTjWagO0OXpFPaCfQDwKaW9Y1VW6v9wI7MHM/MrwNfoRnwi2Iq0J2hS9Ip7QT6LmBLRFwdEcuBW4EdM/p8hObsnIjop7kLZk/nyjzdxNQM3bNcJGnavIGemRPA7cD9wOPAfZm5OyLuioibq273A6MR8RjwAPBzmTm6WEU3nKFL0hnmPW0RIDN3AjtntN3ZspzAz1avRTc1Q/fSf0k6pZZXijZyKtBrWb4kLYpaJuKpGXqXC5GkC0gtI7HRcIYuSTPVMhGdoUvSmWoZiZPO0CXpDLVMRC8skqQz1TrQvbBIkk6pdaA7Q5ekU2oZ6BONBuCFRZLUqpaBPnp0DDDQJalVLQP9Fz78ZQDWrGjrzgWS1BNqGegnxif59udewjWb1na7FEm6YNQu0McnGzx9YoLXb73CB1xIUovaBfqR48395+vXLO9yJZJ0YaldoE8dEF23ZkWXK5GkC0vtAv3wsalAd4YuSa1qF+ijVaCvv8hAl6RWtQv0w0dPAs7QJWmm2gX65Zes5Lu39HPZagNdklrV7sqcm166gZteuqHbZUjSBad2M3RJ0uwMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFaCvQI2JbRDwREcMRccdZ+r0xIjIiBjtXoiSpHfMGekT0AXcDNwJbgdsiYuss/S4G3g58vtNFSpLm184M/TpgODP3ZOYYcC9wyyz9fg34LeBEB+uTJLWpnUC/EtjXsr6/apsWEdcCmzLzY2f7oIjYHhFDETE0MjJyzsVKkub2rA+KRsQS4PeAd8zXNzPvyczBzBwcGBh4tl8tSWrRTqAfADa1rG+s2qZcDLwE+ExE7AWuB3Z4YFSSzq92An0XsCUiro6I5cCtwI6pjZn5VGb2Z+bmzNwMPAjcnJlDi1KxJGlW8wZ6Zk4AtwP3A48D92Xm7oi4KyJuXuwCJUntaeuZopm5E9g5o+3OOfre8OzLkiSdK68UlaRCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEG0FekRsi4gnImI4Iu6YZfvPRsRjEfGliPhURDyv86VKks5m3kCPiD7gbuBGYCtwW0RsndHti8BgZr4M+HvgtztdqCTp7NqZoV8HDGfmnswcA+4FbmntkJkPZObxavVBYGNny5QkzaedQL8S2Neyvr9qm8vbgI/PtiEitkfEUEQMjYyMtF+lJGleHT0oGhFvAQaB35lte2bek5mDmTk4MDDQya+WpJ63tI0+B4BNLesbq7bTRMTrgF8EviczT3amPElSu9qZoe8CtkTE1RGxHLgV2NHaISKuAf4UuDkzD3a+TEnSfOYN9MycAG4H7gceB+7LzN0RcVdE3Fx1+x3gIuBDEfFwROyY4+MkSYuknV0uZOZOYOeMtjtbll/X4bokSefIK0UlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUiKXdLuCcfeUTsPvDsP4FsP6Fzde658Py1d2uTJK6qn6B/vQB2PMAPPK3p7dfcmUz5C/dBKsug9XrYNW605dXXAzL1zRfS1dCRHfGIEmLoH6BPvjW5uvkUTi8B0aHYfRr1c9h2PMZOH4YJp45++fEEli2pjmzX76mWq7W+1ZA3zLoW958LV1+ann6VW1f2tI3+mDJUliypGW5r1rua1mer31J8ydR/aUTzbbp5ZlttNlvljb/UpOK0VagR8Q24D1AH/DnmfmbM7avAN4PfAcwCrwpM/d2ttQZVlwEG17WfM1m/Bl45kgz3J853Pw5dqz5Gj92anm67TiMHYUTT8PkWMtrvPlz4uSp5cb4og6tO+YI/unNcXrfWdvPtq2d9llqWvBndbKuhXwWC7SANy74L+UFvu98fl8txraA99zw8/CSNy7gu85u3kCPiD7gbuD1wH5gV0TsyMzHWrq9DTiSmS+MiFuB3wLe1PFqz8WyVc3XJc/t/Gdnzh74OQmN6pWT0Jiolhsty1PtjZb+E7O8dxLI5ndN/cxGSxuztM3Wb2YbbfbL1gGfPvbZ2s/YNlf7Qj6rjfect8+aY+Vs72nXXH9+Z3/Tefyu8/19NRjbQr9r5dqFvW8e7czQrwOGM3MPQETcC9wCtAb6LcCvVMt/D/xhRETmgv9kL2wRzV0tS1d0uxJJmtbOaYtXAvta1vdXbbP2ycwJ4Clg/cwPiojtETEUEUMjIyMLq1iSNKvzeh56Zt6TmYOZOTgwMHA+v1qSitdOoB8ANrWsb6zaZu0TEUuBS2keHJUknSftBPouYEtEXB0Ry4FbgR0z+uwAfrRa/kHg08XuP5ekC9S8B0UzcyIibgfup3na4nszc3dE3AUMZeYO4C+AD0TEMHCYZuhLks6jts5Dz8ydwM4ZbXe2LJ8AfqizpUmSzoU355KkQhjoklSI6Naxy4gYAb6xwLf3A4c6WE4dOObe4Jh7w7MZ8/Myc9bzvrsW6M9GRAxl5mC36zifHHNvcMy9YbHG7C4XSSqEgS5JhahroN/T7QK6wDH3BsfcGxZlzLXchy5JOlNdZ+iSpBkMdEkqRO0CPSK2RcQTETEcEXd0u55OiYj3RsTBiHi0pW1dRHwyIr5a/bysao+I+IPqz+BLEXFt9ypfuIjYFBEPRMRjEbE7It5etRc77ohYGRFfiIhHqjH/atV+dUR8vhrb31U3wiMiVlTrw9X2zV0dwAJFRF9EfDEiPlqtFz1egIjYGxFfjoiHI2KoalvU3+1aBXrL4/BuBLYCt0XE1u5W1THvA7bNaLsD+FRmbgE+Va1Dc/xbqtd24I/PU42dNgG8IzO3AtcDP1X99yx53CeB783MlwOvALZFxPU0H9v47sx8IXCE5mMdoeXxjsC7q3519Hbg8Zb10sc75TWZ+YqWc84X93c7M2vzAr4TuL9l/V3Au7pdVwfHtxl4tGX9CWBDtbwBeKJa/lPgttn61fkF/BPNZ9f2xLiB1cB/Aq+kedXg0qp9+vec5l1Ov7NaXlr1i27Xfo7j3FiF1/cCH6X5VOVix9sy7r1A/4y2Rf3drtUMnfYeh1eSKzLzm9Xyt4ArquXi/hyqf1pfA3yewsdd7X54GDgIfBL4GvBkNh/fCKePq63HO17gfh94J9Co1tdT9ninJPCJiHgoIrZXbYv6u93W7XPVfZmZEVHkOaYRcRHwD8DPZObTETG9rcRxZ+Yk8IqIWAt8GHhRdytaPBHxfcDBzHwoIm7ocjnn26sz80BEXA58MiL+q3XjYvxu122G3s7j8EryvxGxAaD6ebBqL+bPISKW0Qzzv8nMf6yaix83QGY+CTxAc5fD2urxjXD6uOr+eMdXATdHxF7gXpq7Xd5DueOdlpkHqp8Haf7FfR2L/Ltdt0Bv53F4JWl9tN+P0tzHPNX+I9WR8euBp1r+GVcb0ZyK/wXweGb+XsumYscdEQPVzJyIWEXzmMHjNIP9B6tuM8dc28c7Zua7MnNjZm6m+f/rpzPzzRQ63ikRsSYiLp5aBt4APMpi/253+8DBAg403AR8heZ+x1/sdj0dHNcHgW8C4zT3n72N5r7DTwFfBf4FWFf1DZpn+3wN+DIw2O36FzjmV9Pcz/gl4OHqdVPJ4wZeBnyxGvOjwJ1V+/OBLwDDwIeAFVX7ymp9uNr+/G6P4VmM/Qbgo70w3mp8j1Sv3VNZtdi/2176L0mFqNsuF0nSHAx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVIj/B/MGeid5+s9lAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(torch.tensor(accurate_list))\n",
    "plt.plot(torch.tensor(loss_list))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7067b0e8-f962-4870-b5a6-62a1836f5dbb",
   "metadata": {},
   "source": [
    "### Evaluate the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "745832f9-0a4f-4982-93e8-7206d203f5ba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ContactMapRegression(\n",
       "  (regression): Linear(in_features=484, out_features=2, bias=True)\n",
       "  (activation): Sigmoid()\n",
       ")"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contactmodel=torch.load('../models/contactmap_filter_planb.pt')\n",
    "contactmodel.eval()\n",
    "contactmodel.to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "f70b3fd5-bf42-4e83-9644-8413c06d17bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_158084/2454322881.py:3: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  lbdall_pred_contact=torch.tensor(lbdall_pred_contact,dtype=torch.float32).to(device)\n"
     ]
    }
   ],
   "source": [
    "index=torch.where(predictions!=0)\n",
    "lbdall_pred_contact=get_pred_contact(lbdin_outputs[index])\n",
    "lbdall_pred_contact=torch.tensor(lbdall_pred_contact,dtype=torch.float32).to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "46d1c357-b497-49bb-8e35-d5e7d4f0326d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_158084/2803984175.py:5: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  labels=torch.tensor(labels,dtype=torch.long).numpy()\n"
     ]
    }
   ],
   "source": [
    "predictions=contactmodel(lbdall_pred_contact)\n",
    "_,prediction=predictions.max(dim=1)\n",
    "prediction=prediction.cpu()\n",
    "labels=torch.cat([torch.ones(2),torch.zeros(9)])\n",
    "labels=torch.tensor(labels,dtype=torch.long).numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "3ea9d395-9362-4e33-a79b-51e2a44ce91b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Confusion matrix: \n",
      "[[7 2]\n",
      " [0 2]]\n",
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.78      1.00      0.88         7\n",
      "           1       1.00      0.50      0.67         4\n",
      "\n",
      "    accuracy                           0.82        11\n",
      "   macro avg       0.89      0.75      0.77        11\n",
      "weighted avg       0.86      0.82      0.80        11\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.metrics import classification_report\n",
    "confusion_mat = confusion_matrix(labels, prediction.cpu())\n",
    "report = classification_report(prediction.cpu(), labels)\n",
    "print(f\"Confusion matrix: \\n{confusion_mat}\")\n",
    "print(report)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "187ea7c6-2710-493a-bdca-aed159b508d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Save the model\n",
    "# torch.save(contactmodel, './models/contactmap_filter_planb.pt')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "AMP (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}