SpeCollate.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from pandas import ExcelWriter\n",
    "from pandas import ExcelFile\n",
    "import numpy as np\n",
    "import random as rand\n",
    "import math\n",
    "import queue\n",
    "import csv\n",
    "from collections import OrderedDict\n",
    "from IPython.display import clear_output\n",
    "import csv\n",
    "from heapq import merge\n",
    "from sklearn import preprocessing\n",
    "import gc\n",
    "import os\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "import shutil\n",
    "from configparser import ConfigParser\n",
    "import ast\n",
    "import sys\n",
    "import re\n",
    "from pathlib import Path\n",
    "import pickle\n",
    "\n",
    "from sklearn.datasets import make_circles, make_moons\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import torch.optim as optim\n",
    "import torchvision\n",
    "import torchvision.transforms as transforms\n",
    "import torchvision.datasets as datasets\n",
    "from IPython.core.debugger import set_trace\n",
    "from torch.utils import data\n",
    "from torch.utils.data.sampler import SubsetRandomSampler\n",
    "from torch.autograd import Variable\n",
    "\n",
    "rand.seed(37)\n",
    "\n",
    "#from src.snapconfig import config\n",
    "from src.snaputils import simulatespectra as sim\n",
    "from src.snaptrain import process\n",
    "from src.snaputils import reader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Temporary config func. Original one in the project.\n",
    "class config:\n",
    "    \"\"\"Define constants\"\"\"\n",
    "    AAMass = OrderedDict([('A', 71.037114), ('C', 103.009185), ('D', 115.026943), ('E', 129.042593),\n",
    "                          ('F', 147.068414), ('G', 57.021464), ('H', 137.058912), ('I', 113.084064),\n",
    "                          ('K', 128.094963), ('L', 113.084064), ('M', 131.040485), ('N', 114.042927),\n",
    "                          ('P', 97.052764), ('Q', 128.058578), ('R', 156.101111), ('S', 87.032028),\n",
    "                          ('T', 101.047679), ('V', 99.068414), ('W', 186.079313), ('Y', 163.0633),\n",
    "                          ('p', 79.97), ('o', 15.99), ('h', 0.98), ('c', 57.02), ('a', 42.01),\n",
    "                          ('r', -17.03), ('y', 43.01), ('d', -18.01), ('t', 26.02)])\n",
    "    \n",
    "    aas = ['_PAD'] + list(AAMass.keys())\n",
    "    aa2idx = {a:i for i, a in enumerate(aas)}\n",
    "\n",
    "    ModMass = {\"Oxidation\": 15.994915, \"CAM\": 57.02146, \"Carbamidomethyl\": 57.02146, \"ICAT_light\": 227.12,\n",
    "               \"ICAT_heavy\": 236.12, \"AB_old_ICATd0\": 442.20, \"AB_old_ICATd8\": 450.20, \"Acetyl\": 42.0106,\n",
    "               \"Deamidation\": 0.9840, \"Pyro-cmC\": -17.026549, \"Pyro-glu\": -17.026549, \"Pyro_glu\": -18.010565,\n",
    "               \"Amide\": -0.984016, \"Phospho\": 79.9663, \"Methyl\": 14.0157, \"Carbamyl\": 43.00581}\n",
    "\n",
    "    ModCHAR = OrderedDict([(\"15.99\", \"o\"), (\"0.98\", \"h\"), (\"42.01\", \"a\"), (\"-17.03\", \"r\"),\n",
    "                           (\"79.97\", \"p\"), (\"43.01\", \"y\"), (\"-18.01\", \"d\"), (\"26.02\", \"t\"),\n",
    "                           (\"57.02\", \"c\")\n",
    "                          ])\n",
    "    # ModCHAR = {\"15.99\": \"o\", \"0.98\": \"h\", \"57.02\": \"c\", \"42.01\": \"a\", \"-17.03\": \"r\", \"79.97\": \"p\"}\n",
    "    Ignore = [\"U\", \"X\"]\n",
    "    Mods = [{\"mod_char\": \"p\", \"aas\": [\"S\", \"T\", \"Y\"]}\n",
    "            # {\"mod_char\": \"o\", \"aas\": [\"nt\", \"M\"]}\n",
    "           ]\n",
    "    H2O = 18.015\n",
    "    NH3 = 17.031\n",
    "    PROTON = 1.00727647\n",
    "    DEFAULT_PARAM_PATH = os.path.join(os.getcwd(), 'config.ini')\n",
    "    PARAM_PATH = None\n",
    "    l_config = None\n",
    "\n",
    "\n",
    "    def get_config(section='input', key=None):\n",
    "        \"\"\"Read the configuration parameters and return a dictionary.\"\"\"\n",
    "\n",
    "        # If file path is given use it otherwise use default.\n",
    "        file_path = config.PARAM_PATH if config.PARAM_PATH else config.DEFAULT_PARAM_PATH\n",
    "\n",
    "        # Read config and convert each value to appropriate type.\n",
    "        # Only for the first time.\n",
    "        if not config.l_config:\n",
    "            config.l_config = dict()\n",
    "            config_ = ConfigParser()\n",
    "            assert isinstance(file_path, str)\n",
    "            config_.read(file_path)\n",
    "            for section_ in config_.sections():\n",
    "                config.l_config[section_] = dict()\n",
    "                for key_ in config_[section_]:\n",
    "                    try:\n",
    "                        config.l_config[section_][key_] = ast.literal_eval(config_[section_][key_])\n",
    "                    except (ValueError, SyntaxError):\n",
    "                        config.l_config[section_][key_] = config_[section_][key_]\n",
    "\n",
    "        if section and section in config.l_config:\n",
    "            if key and key in config.l_config[section]:\n",
    "                return config.l_config[section][key]\n",
    "            return config.l_config[section]\n",
    "        return config.l_config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocess_msps(msp_dir, out_dir):\n",
    "    in_path = Path(msp_dir)\n",
    "    assert in_path.exists() and in_path.is_dir()\n",
    "    \n",
    "    msp_files = [join(msp_dir, f) for f in listdir(msp_dir) if\n",
    "                 isfile(join(msp_dir, f)) and f.split('.')[-1] == 'msp']\n",
    "    assert len(msp_files) > 0\n",
    "    \n",
    "    out_path = Path(out_dir)\n",
    "    if out_path.exists() and out_path.is_dir():\n",
    "        shutil.rmtree(out_path)\n",
    "    out_path.mkdir()\n",
    "    Path(join(out_path, 'spectra')).mkdir()\n",
    "    Path(join(out_path, 'peptides')).mkdir()\n",
    "        \n",
    "    print('reading {} files'.format(len(msp_files)))\n",
    "    \n",
    "    count = 0\n",
    "    max_peaks = max_moz = 0\n",
    "    for species_id, msp_file in enumerate(msp_files):\n",
    "        print('Reading: {}'.format(msp_file))\n",
    "        \n",
    "        f = open(msp_file, \"r\")\n",
    "        lines = f.readlines()\n",
    "        f.close()\n",
    "\n",
    "        pep_list = []\n",
    "        dataset = []\n",
    "        label = []\n",
    "\n",
    "        # FIXME: config should use only one get_config call.\n",
    "        spec_size = config.get_config(section='input', key='spec_size')\n",
    "        charge = config.get_config(section='input', key='charge')\n",
    "        use_mods = config.get_config(section='input', key='use_mods')\n",
    "        num_species = config.get_config(section='input', key='num_species')\n",
    "        seq_len = config.get_config(section='ml', key='pep_seq_len')\n",
    "\n",
    "        print('len of file: ' + str(len(lines)))\n",
    "        limit = 200000\n",
    "        pep = []\n",
    "        spec = []\n",
    "        pep_set = set()\n",
    "        is_name = is_mw = is_num_peaks = False\n",
    "        prev = 0\n",
    "        i = 0\n",
    "        while i < len(lines) and limit > 0:\n",
    "            line = lines[i]\n",
    "            i += 1\n",
    "            if line.startswith('Name:'):\n",
    "                name_groups = re.search(r\"Name:\\s(?P<pep>[a-zA-Z]+)/(?P<charge>\\d+)\"\n",
    "                                        r\"(?:_(?P<num_mods>\\d+)(?P<mods>.*))?\", line)\n",
    "                if not name_groups:\n",
    "                    continue\n",
    "                    \n",
    "                pep = name_groups['pep']\n",
    "                if len(pep) + 1 > seq_len:\n",
    "                    continue\n",
    "                    \n",
    "                l_charge = int(name_groups['charge'])\n",
    "                num_mods = int(name_groups['num_mods'])\n",
    "\n",
    "                is_name = True\n",
    "\n",
    "            if is_name and line.startswith('MW:'):\n",
    "                mass = float(re.findall(r\"MW:\\s([-+]?[0-9]*\\.?[0-9]*)\", line)[0])\n",
    "                if round(mass) < spec_size:\n",
    "                    is_mw = True\n",
    "                    # limit = limit - 1\n",
    "                else:\n",
    "                    is_name = is_mw = is_num_peaks = False\n",
    "                    continue\n",
    "\n",
    "            if is_name and is_mw and line.startswith('Num peaks:'):\n",
    "                num_peaks = int(re.findall(r\"Num peaks:\\s([0-9]*\\.?[0-9]*)\", line)[0])\n",
    "                if num_peaks > max_peaks:\n",
    "                    max_peaks = num_peaks\n",
    "\n",
    "                spec = np.zeros(spec_size)\n",
    "                while lines[i] != '\\n':\n",
    "                    mz_line = lines[i]\n",
    "                    i += 1\n",
    "                    mz_splits = mz_line.split('\\t')\n",
    "                    moz, intensity = float(mz_splits[0]), float(mz_splits[1])\n",
    "                    if moz > max_moz:\n",
    "                        max_moz = moz\n",
    "                    spec[round(moz)] += round(intensity)\n",
    "\n",
    "                # for k in range(1, charge + 1):\n",
    "                #     spec[-k] = 0\n",
    "                # spec[-l_charge] = 1000.0\n",
    "                spec = np.clip(spec, None, 1000.0)\n",
    "                # spec = preprocessing.scale(spec)\n",
    "\n",
    "                is_num_peaks = True\n",
    "\n",
    "            if is_name and is_mw and is_num_peaks:\n",
    "                is_name = is_mw = is_num_peaks = False\n",
    "                \n",
    "                #pep = '{}{}{}'.format(charge, species_id, pep)\n",
    "\n",
    "                \"\"\"output the data to \"\"\"\n",
    "                spec_tensor = torch.tensor((np.asarray(spec) - 3.725) / 51.479, dtype=torch.float)\n",
    "                \n",
    "                torch.save(spec_tensor, \n",
    "                           join(out_dir, 'spectra', '{}-{}-{}-{}-{}.pt'\n",
    "                                .format(count, species_id, mass, l_charge, int(num_mods > 0))))\n",
    "                \n",
    "                pep_file_name = '{}-{}-{}-{}-{}.pep'.format(count, species_id, mass, l_charge, int(num_mods > 0))\n",
    "                    \n",
    "                with open(join(out_path, 'peptides', pep_file_name), 'w+') as f:\n",
    "                    f.write(pep)\n",
    "\n",
    "                count = count + 1\n",
    "                pep = 0\n",
    "                spec = []\n",
    "                new = int((i / len(lines)) * 100)\n",
    "                if new > prev:\n",
    "                    clear_output(wait=True)\n",
    "                    print(str(new) + '%')\n",
    "                    prev = new\n",
    "\n",
    "        print('max peaks: ' + str(max_peaks))\n",
    "        print('count: ' + str(count))\n",
    "        print('max moz: ' + str(max_moz))\n",
    "#         return pep_list, dataset, label\n",
    "#         tmp_pep_list, tmp_dataset, tmp_labels = read_msp(msp_file, species_id, decoy)\n",
    "#         pep_list.extend(tmp_dataset)\n",
    "#         dataset.extend(tmp_dataset)\n",
    "#         label.extend(tmp_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocess_mgfs_back(mgf_dir, out_dir):\n",
    "    in_path = Path(mgf_dir)\n",
    "    assert in_path.exists() and in_path.is_dir()\n",
    "    \n",
    "    mgf_files = [join(mgf_dir, f) for f in listdir(mgf_dir) if\n",
    "                 isfile(join(mgf_dir, f)) and f.split('.')[-1] == 'mgf']\n",
    "    assert len(mgf_files) > 0\n",
    "    \n",
    "    out_path = Path(out_dir)\n",
    "    if out_path.exists() and out_path.is_dir():\n",
    "        shutil.rmtree(out_path)\n",
    "        \n",
    "    out_path.mkdir()\n",
    "    Path(join(out_path, 'spectra')).mkdir()\n",
    "    Path(join(out_path, 'peptides')).mkdir()\n",
    "        \n",
    "    print('reading {} files'.format(len(mgf_files)))\n",
    "    \n",
    "    ch = np.zeros(20)\n",
    "    modified = 0\n",
    "    unmodified = 0\n",
    "    unique_pep_set = set()\n",
    "    \n",
    "    summ = 0\n",
    "    sq_sum = 0\n",
    "    N = 0\n",
    "    \n",
    "    tot_count = 0\n",
    "    max_peaks = max_moz = 0\n",
    "    for species_id, mgf_file in enumerate(mgf_files):\n",
    "        print('Reading: {}'.format(mgf_file))\n",
    "        \n",
    "        f = open(mgf_file, \"r\")\n",
    "        lines = f.readlines()\n",
    "        f.close()\n",
    "        \n",
    "        count = lcount = 0\n",
    "\n",
    "        pep_list = []\n",
    "        dataset = []\n",
    "        label = []\n",
    "        \n",
    "        mass_ign = 0\n",
    "        pep_len_ign = 0\n",
    "        dup_ign = 0\n",
    "\n",
    "        # FIXME: config should use only one get_config call.\n",
    "        spec_size = config.get_config(section='input', key='spec_size')\n",
    "        charge = config.get_config(section='input', key='charge')\n",
    "        use_mods = config.get_config(section='input', key='use_mods')\n",
    "        num_species = config.get_config(section='input', key='num_species')\n",
    "        seq_len = config.get_config(section='ml', key='pep_seq_len')\n",
    "\n",
    "        print('len of file: ' + str(len(lines)))\n",
    "        limit = 200000\n",
    "        pep = []\n",
    "        spec = []\n",
    "        pep_set = set()\n",
    "        is_name = is_mw = is_charge = False\n",
    "        prev = 0\n",
    "        i = 0\n",
    "        while i < len(lines) and limit > 0:\n",
    "            line = lines[i]\n",
    "            i += 1\n",
    "\n",
    "            if line.startswith('PEPMASS'):\n",
    "                count += 1\n",
    "                mass = float(re.findall(r\"PEPMASS=([-+]?[0-9]*\\.?[0-9]*)\", line)[0])\n",
    "                if round(mass) < spec_size:\n",
    "                    is_mw = True\n",
    "                    # limit = limit - 1\n",
    "                else:\n",
    "                    is_name = is_mw = is_charge = False\n",
    "                    mass_ign += 1\n",
    "                    continue\n",
    "            \n",
    "            if is_mw and line.startswith('CHARGE'):\n",
    "                l_charge = int(re.findall(r\"CHARGE=([-+]?[0-9]*\\.?[0-9]*)\", line)[0])\n",
    "                is_charge = True\n",
    "                \n",
    "            if is_mw and is_charge and line.startswith(\"SEQ\"):\n",
    "                line = line.strip()\n",
    "                pep_segs = re.findall(r\"([A-Z]+)\", line)[1:]\n",
    "                num_mods = len(pep_segs) - 1 + (1 if line[-1] == ')' else 0)\n",
    "                pep = ''.join(pep_segs)\n",
    "                mods = re.findall(r\"([A-Z]\\([^\\(\\)]*\\))\", line)\n",
    "                \n",
    "                if len(pep) + 2 > seq_len:\n",
    "                    pep_len_ign += 1\n",
    "                    is_name = is_mw = is_charge = False\n",
    "                    continue\n",
    "                \n",
    "                rtinsec = lines[i-2]\n",
    "                rt = round(float(re.findall(r\"RTINSECONDS=([-+]?[0-9]*\\.?[0-9]*)\", rtinsec)[0]))\n",
    "                check_pep = line[4:] + str(l_charge) + \"-\" + str(rt)\n",
    "                # check_pep = pep + str(l_charge)\n",
    "                if check_pep in pep_set:\n",
    "                    dup_ign += 1\n",
    "                    is_mw = is_charge = is_name = False\n",
    "                    continue\n",
    "                else:\n",
    "                    ch[l_charge] += 1\n",
    "                    if num_mods > 0:\n",
    "                        modified += 1\n",
    "                    else:\n",
    "                        unmodified += 1\n",
    "                    if pep not in unique_pep_set:\n",
    "                        unique_pep_set.add(pep)\n",
    "                    pep_set.add(check_pep)\n",
    "                    \n",
    "\n",
    "                ind = [] # setting the precision to one decimal point.\n",
    "                val = []\n",
    "                #for ch_val in range(l_charge):\n",
    "                #    ind.append(ch_val)\n",
    "                #    val.append(1)\n",
    "                    \n",
    "                while 'END IONS' not in lines[i].upper():\n",
    "                    if lines[i] == '\\n':\n",
    "                        continue\n",
    "                    mz_line = lines[i]\n",
    "                    i += 1\n",
    "                    mz_splits = mz_line.split(' ')\n",
    "                    moz, intensity = float(mz_splits[0]), float(mz_splits[1])\n",
    "                    if moz > max_moz:\n",
    "                        max_moz = moz\n",
    "                    if 0 < round(moz) < spec_size:\n",
    "                        # spec[round(moz*10)] += round(intensity)\n",
    "                        ind.append(round(moz*10))\n",
    "                        val.append(intensity)\n",
    "                        \n",
    "                        \n",
    "                # for k in range(1, charge + 1):\n",
    "                #     spec[-k] = 0\n",
    "                # spec[-l_charge] = 1000.0\n",
    "                \n",
    "                ind = np.array(ind)\n",
    "                val = np.clip(val, None, 10000.0)\n",
    "                assert len(ind) == len(val)\n",
    "                spec = np.array([ind, val])\n",
    "                \n",
    "                summ += sum(val)\n",
    "                sq_sum += sum(s_val*s_val for s_val in val)\n",
    "                N += 80000\n",
    "                # spec = preprocessing.scale(spec)\n",
    "\n",
    "                is_name = True\n",
    "\n",
    "            if is_name and is_mw and is_charge:\n",
    "                is_name = is_mw = is_charge = False\n",
    "                \n",
    "                #pep = '{}{}{}'.format(charge, species_id, pep)\n",
    "\n",
    "                \"\"\"output the data to \"\"\"\n",
    "                #spec_tensor = torch.tensor((np.asarray(spec) - 3.725) / 51.479, dtype=torch.float)\n",
    "                # ind = torch.LongTensor([[0]*len(ind), ind])\n",
    "                # val = torch.FloatTensor(val)\n",
    "                # spec_tensor = torch.sparse.FloatTensor(ind, val, torch.Size([1, 80000]))\n",
    "                \n",
    "                # spec_tensor = torch.tensor(np.asarray(spec), dtype=torch.float)\n",
    "                \n",
    "#                 torch.save(spec_tensor, \n",
    "#                            join(out_dir, 'spectra', '{}-{}-{}-{}-{}.pt'\n",
    "#                                 .format(tot_count, species_id, mass, l_charge, int(num_mods > 0))))\n",
    "                \n",
    "                np.save(join(out_dir, 'spectra', '{}-{}-{}-{}-{}.npy'\n",
    "                             .format(tot_count, species_id, mass, l_charge, int(num_mods > 0))),\n",
    "                       spec)\n",
    "                \n",
    "                pep_file_name = '{}-{}-{}-{}-{}.pep'.format(tot_count, species_id, mass, l_charge, int(num_mods > 0))\n",
    "                    \n",
    "                with open(join(out_dir, 'peptides', pep_file_name), 'w+') as f:\n",
    "                    f.write(pep)\n",
    "\n",
    "                lcount += 1\n",
    "                tot_count += 1\n",
    "                \n",
    "                pep = 0\n",
    "                spec = []\n",
    "                new = int((i / len(lines)) * 100)\n",
    "                if new >= prev + 10:\n",
    "                    #clear_output(wait=True)\n",
    "                    print('count: ' + str(lcount))\n",
    "                    print(str(new) + '%')\n",
    "                    prev = new\n",
    "\n",
    "        #print('max peaks: ' + str(max_peaks))\n",
    "        print('In current file, read {} out of {}'.format(lcount, count))\n",
    "        print(\"Ignored: large mass: {}, pep len: {}, dup: {}\".format(mass_ign, pep_len_ign, dup_ign))\n",
    "        print('overall running count: ' + str(tot_count))\n",
    "        print('max moz: ' + str(max_moz))\n",
    "#         return pep_list, dataset, label\n",
    "#         tmp_pep_list, tmp_dataset, tmp_labels = read_msp(msp_file, species_id, decoy)\n",
    "#         pep_list.extend(tmp_dataset)\n",
    "#         dataset.extend(tmp_dataset)\n",
    "#         label.extend(tmp_labels)\n",
    "    print(\"Statistics:\")\n",
    "    print(\"Charge distribution:\")\n",
    "    print(ch)\n",
    "    print(\"Modified:\\t{}\".format(modified))\n",
    "    print(\"Unmodified:\\t{}\".format(unmodified))\n",
    "    print(\"Unique Peptides:\\t{}\".format(len(unique_pep_set)))\n",
    "    print(\"Sum: {}\".format(summ))\n",
    "    print(\"Sum-Squared: {}\".format(sq_sum))\n",
    "    print(\"N: {}\".format(N))\n",
    "    mean = summ / N\n",
    "    print(\"mean: {}\".format(mean))\n",
    "    std = math.sqrt((sq_sum / N) - mean**2)\n",
    "    print(\"std: {}\".format(std))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_out_dir(dir_path, exist_ok=True):\n",
    "    out_path = Path(dir_path)\n",
    "    if out_path.exists() and out_path.is_dir():\n",
    "        if not exist_ok:\n",
    "            shutil.rmtree(out_path)\n",
    "            out_path.mkdir()\n",
    "    else:\n",
    "        out_path.mkdir()\n",
    "        \n",
    "    Path(join(out_path, 'spectra')).mkdir()\n",
    "    Path(join(out_path, 'peptides')).mkdir()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_in_dir(dir_path, ext, ignore_list):\n",
    "    in_path = Path(dir_path)\n",
    "    assert in_path.exists() and in_path.is_dir()\n",
    "    \n",
    "    files = [join(dir_path, f) for f in listdir(dir_path) if\n",
    "                 isfile(join(dir_path, f)) and not f.startswith('.') \n",
    "                 and f.split('.')[-1] == ext and f not in ignore_list]\n",
    "    assert len(files) > 0\n",
    "    return files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def isfloat(str_float):\n",
    "    try:\n",
    "        float(str_float)\n",
    "        return True\n",
    "    except ValueError: \n",
    "        return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mod_repl(match):\n",
    "    lookup = str(round(float(match.group(0)), 2))\n",
    "    return config.ModCHAR[lookup] if lookup in config.ModCHAR else \"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mod_filt(pep, mods, count):\n",
    "    is_valid = True\n",
    "    current_count = 0\n",
    "    for aa in pep:\n",
    "        if aa.islower() and aa != \"c\":\n",
    "            if aa not in mods:\n",
    "                is_valid = False\n",
    "                break\n",
    "            else:\n",
    "                current_count += 1\n",
    "                if current_count > count:\n",
    "                    is_valid = False\n",
    "                    break\n",
    "    \n",
    "    return is_valid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gray_code(num):\n",
    "    return num ^ (num >> 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decimal_to_binary_array(num, arr_len):\n",
    "    bin_arr = [float(i) for i in list('{0:0b}'.format(num))]\n",
    "    assert len(bin_arr) <= arr_len\n",
    "    res = [0.] * (arr_len - len(bin_arr)) + bin_arr\n",
    "    inds = [int(i) for i, _ in enumerate(res) if res[i] > 0.1] # greater than zero. 0.1 for the floating pointing errors.\n",
    "    vals = [1.0] * len(inds)\n",
    "    return inds, vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocess_mgfs(mgf_dir, out_dir):\n",
    "    \n",
    "    mgf_files = verify_in_dir(mgf_dir, \"mgf\", [])\n",
    "    create_out_dir(out_dir, exist_ok=False)\n",
    "        \n",
    "    print('reading {} files'.format(len(mgf_files)))\n",
    "    \n",
    "    spec_size = config.get_config(section='input', key='spec_size')\n",
    "    print(\"spec size: {}\".format(spec_size))\n",
    "    charge = config.get_config(section='input', key='charge')\n",
    "    use_mods = config.get_config(section='input', key='use_mods')\n",
    "    num_species = config.get_config(section='input', key='num_species')\n",
    "    seq_len = config.get_config(section='ml', key='pep_seq_len')\n",
    "    aa2idx = config.aa2idx\n",
    "    \n",
    "    non_mod_c = 0\n",
    "    \n",
    "    ch = np.zeros(20)\n",
    "    modified = 0\n",
    "    unmodified = 0\n",
    "    unique_pep_set = set()\n",
    "    \n",
    "    pep_dict = {}\n",
    "    idx_spec_map = []\n",
    "    pep_spec = []\n",
    "    pep_idx = 0\n",
    "    \n",
    "#     pep_lens = np.zeros(64)\n",
    "#     missed_cleavs_vec = np.zeros(17)\n",
    "#     aas = np.zeros(30)\n",
    "    aux_summ = np.zeros(5)\n",
    "    aux_sq_sum = np.zeros(5)\n",
    "    aux_N = 0\n",
    "\n",
    "#     pairs = np.zeros(900)\n",
    "    \n",
    "    summ = np.zeros(spec_size + 5)\n",
    "    sq_sum = np.zeros(spec_size + 5)\n",
    "    N = 0\n",
    "    \n",
    "    tot_count = 0\n",
    "    max_peaks = max_moz = max_missed_cleavs = 0\n",
    "    for species_id, mgf_file in enumerate(mgf_files):\n",
    "        print('Reading: {}'.format(mgf_file))\n",
    "        \n",
    "        f = open(mgf_file, \"r\")\n",
    "        lines = f.readlines()\n",
    "        f.close()\n",
    "        \n",
    "        count = lcount = 0\n",
    "        \n",
    "        pep_list = []\n",
    "        dataset = []\n",
    "        label = []\n",
    "        \n",
    "        mass_ign = 0\n",
    "        pep_len_ign = 0\n",
    "        dup_ign = 0\n",
    "\n",
    "        print('len of file: ' + str(len(lines)))\n",
    "        limit = 200000\n",
    "        pep = []\n",
    "        spec = []\n",
    "        pep_set = set()\n",
    "        is_name = is_mw = is_charge = False\n",
    "        prev = 0\n",
    "        i = 0\n",
    "        while i < len(lines) and limit > 0:\n",
    "            line = lines[i]\n",
    "            i += 1\n",
    "\n",
    "            if line.startswith('PEPMASS'):\n",
    "                count += 1\n",
    "                mass = float(re.findall(r\"PEPMASS=([-+]?[0-9]*\\.?[0-9]*)\", line)[0])\n",
    "                is_mw = True\n",
    "#                 if round(mass)*10 < spec_size:\n",
    "#                     is_mw = True\n",
    "#                     # limit = limit - 1\n",
    "#                 else:\n",
    "#                     is_name = is_mw = is_charge = False\n",
    "#                     mass_ign += 1\n",
    "#                     continue\n",
    "            \n",
    "            if is_mw and line.startswith('CHARGE'):\n",
    "                l_charge = int(re.findall(r\"CHARGE=([-+]?[0-9]*\\.?[0-9]*)\", line)[0])\n",
    "                mass = (mass - config.PROTON) * l_charge\n",
    "                is_charge = True\n",
    "                if l_charge > charge or round(mass*10) > spec_size:\n",
    "                    is_name = is_mw = is_charge = False\n",
    "                    continue\n",
    "                \n",
    "            if is_mw and is_charge and line.startswith(\"SEQ\"):\n",
    "                line = re.sub(r\"[()]\", \"\", line.strip()).split('=')[-1]\n",
    "                mod_repl_rex = r'([-+]?\\d*\\.\\d+|[-+]?\\d+)'\n",
    "                pep, num_mods = re.subn(mod_repl_rex, mod_repl, line)\n",
    "#                 missed_cleavs = (pep.count(\"K\") + pep.count(\"R\")) - (pep.count(\"KP\") + pep.count(\"RP\")) - 1\n",
    "#                 max_missed_cleavs = max(missed_cleavs, max_missed_cleavs)\n",
    "                    \n",
    "#                 if re.search(r\"([a-z]{2,})\", pep):\n",
    "#                     print(pep)\n",
    "                mods = [\"p\", \"o\"]\n",
    "                count = 3\n",
    "                #if len(pep) + 2 > seq_len or \"O\" in pep or \"U\" in pep or re.search(r\"([a-z]{2,})\", pep) or not mod_filt(pep, mods, count):\n",
    "                if len(pep) + 2 > seq_len or \"O\" in pep or \"U\" in pep or re.search(r\"([a-z]{2,})\", pep):\n",
    "                    pep_len_ign += 1\n",
    "                    is_name = is_mw = is_charge = False\n",
    "                    continue\n",
    "                \n",
    "                p_charge = [0] * 5\n",
    "                for p_ch in range(l_charge):\n",
    "                    p_charge[p_ch] = 1\n",
    "                pep = ''.join(map(str, p_charge.copy())) + pep\n",
    "                p_charge = [0] * 5\n",
    "                #pep, num_mods = re.subn(mod_repl_rex, \"\", line)\n",
    "                \n",
    "                # rtinsec = lines[i-2]\n",
    "                # rt = round(float(re.findall(r\"RTINSECONDS=([-+]?[0-9]*\\.?[0-9]*)\", rtinsec)[0]))\n",
    "                # check_pep = line[4:] + str(l_charge)#  + \"-\" + str(rt)\n",
    "                # check_pep = pep + str(l_charge)\n",
    "                sub_spec_count = 1\n",
    "                pep_file_name = ''\n",
    "                spec_file_name = ''\n",
    "                \n",
    "                ch[l_charge] += 1\n",
    "                if num_mods > 0:\n",
    "                    modified += 1\n",
    "                else:\n",
    "                    unmodified += 1\n",
    "                        \n",
    "                if pep not in unique_pep_set:\n",
    "                    unique_pep_set.add(pep)\n",
    "#                     aa_len = len(config.aas)\n",
    "#                     for a1, a2 in zip(pep[:-1], pep[1:]):\n",
    "#                         pair_id = (aa2idx[a1] * aa_len) + aa2idx[a2]\n",
    "#                         pairs[pair_id] += 1\n",
    "#                     pep_lens[len(pep)] = 1\n",
    "#                     missed_cleavs_vec[missed_cleavs] = 1\n",
    "#                     for a in pep:\n",
    "#                         aa_id = config.aa2idx[a]\n",
    "#                         aas[aa_id + 1] += 1\n",
    "                    \n",
    "                    aux_N += 1\n",
    "#                     concat_vec = np.concatenate((pep_lens, missed_cleavs_vec, aas))\n",
    "                    aux_summ += np.array(p_charge)\n",
    "                    aux_sq_sum += np.array(p_charge)**2\n",
    "#                     pairs = np.zeros(900)\n",
    "#                     pep_lens = np.zeros(64)\n",
    "#                     missed_cleavs_vec = np.zeros(17)\n",
    "#                     aas = np.zeros(30)\n",
    "                    # pep_set.add(check_pep)\n",
    "\n",
    "                while not isfloat(re.split(' |\\t|=', lines[i])[0]):\n",
    "                    i += 1\n",
    "                \n",
    "                spec_ind = []\n",
    "                spec_val = []\n",
    "                num_peaks = 0\n",
    "                while 'END IONS' not in lines[i].upper():\n",
    "                    if lines[i] == '\\n':\n",
    "                        i += 1\n",
    "                        continue\n",
    "                    mz_line = lines[i]\n",
    "                    i += 1\n",
    "                    num_peaks += 1\n",
    "                    \n",
    "                    mz_splits = re.split(' |\\t', mz_line)\n",
    "                    moz = round(float(mz_splits[0]) * 10) # 32 because charge is len 8 and mass is len 24\n",
    "                    intensity = math.sqrt(float(mz_splits[1]) + 1.0) # adding 1 to avoid sqrt of zero\n",
    "                    if moz > max_moz:\n",
    "                        max_moz = moz\n",
    "                    if 0 < moz < spec_size:\n",
    "                        # spec[round(moz*10)] += round(intensity)\n",
    "                        if spec_ind and spec_ind[-1] == moz:\n",
    "                            spec_val[-1] = max(intensity, spec_val[-1])\n",
    "                        else:\n",
    "                            spec_ind.append(moz)\n",
    "                            spec_val.append(intensity) # adding one to avoid sqrt of zero\n",
    "                if num_peaks < 15:\n",
    "                    is_name = is_mw = is_charge = False\n",
    "                    continue\n",
    "                    \n",
    "                if pep not in pep_dict:\n",
    "                    pep_dict[pep] = pep_idx\n",
    "                    idx_spec_map.append(sub_spec_count)\n",
    "                    assert len(idx_spec_map) == len(pep_dict)\n",
    "                    tot_count = pep_idx\n",
    "                    pep_file_name = '{}-{}.pep'.format(tot_count, sim.get_pep_mass(pep[5:]))\n",
    "                    spec_file_name = '{}.{}-{}-{}-{}-{}.npy'.format(\n",
    "                        tot_count, sub_spec_count, species_id, mass, l_charge, num_mods)\n",
    "                    pep_spec.append([pep_file_name, [spec_file_name]])\n",
    "                    pep_idx += 1\n",
    "                else:\n",
    "                    tot_count = pep_dict[pep]\n",
    "                    sub_spec_count = idx_spec_map[tot_count] + 1\n",
    "                    idx_spec_map[tot_count] = sub_spec_count\n",
    "                    spec_file_name = '{}.{}-{}-{}-{}-{}.npy'.format(\n",
    "                        tot_count, sub_spec_count, species_id, mass, l_charge, num_mods)\n",
    "                    pep_spec[tot_count][1].append(spec_file_name)\n",
    "                # for k in range(1, charge + 1):\n",
    "                #     spec[-k] = 0\n",
    "                # spec[-l_charge] = 1000.0\n",
    "                \n",
    "                # val = preprocessing.scale(val)\n",
    "                \n",
    "                # TODO: also try with clipping\n",
    "                spec_ind = np.array(spec_ind)\n",
    "                spec_val = np.array(spec_val)\n",
    "                spec_val = (spec_val - np.amin(spec_val)) / (np.amax(spec_val) - np.amin(spec_val))\n",
    "                \n",
    "                gray_mass = gray_code(round(mass * 10000))\n",
    "                mass_ind, mass_val = decimal_to_binary_array(gray_mass, 27) # mass is gray encoded into array of length 24\n",
    "#                 ch_ind,   ch_val   = [l_charge - 1], [1.0]\n",
    "                ch_ind = []\n",
    "                ch_val = []\n",
    "                for l_ch in range(l_charge):\n",
    "                    ch_ind.append(l_ch)\n",
    "                    ch_val.append(1)\n",
    "#                 for ch_val in range(l_charge):\n",
    "#                     ind.append(ch_val+24)\n",
    "#                     val.append(1)\n",
    "                # Join mass, charge, and spec indicies and values.\n",
    "#                 ind = np.concatenate((np.array(mass_ind), np.array(ch_ind) + 27, spec_ind + 35))\n",
    "#                 val = np.concatenate((np.array(mass_val), np.array(ch_val), spec_val))\n",
    "\n",
    "                # Join only charge inside the spectrum\n",
    "                ind = np.concatenate((np.array(ch_ind), spec_ind + 5))\n",
    "                val = np.concatenate((np.array(ch_val), spec_val))\n",
    "            \n",
    "#                 ind = spec_ind\n",
    "#                 val = spec_val\n",
    "        \n",
    "                    \n",
    "                assert len(ind) == len(val)\n",
    "                spec = np.array([ind, val])\n",
    "                \n",
    "                summ[ind] += val\n",
    "                sq_sum[ind] += val**2\n",
    "                N += 1\n",
    "\n",
    "                is_name = True\n",
    "\n",
    "            if is_name and is_mw and is_charge:\n",
    "                is_name = is_mw = is_charge = False\n",
    "                \n",
    "                #pep = '{}{}{}'.format(charge, species_id, pep)\n",
    "\n",
    "                \"\"\"output the data to \"\"\"\n",
    "                #spec_tensor = torch.tensor((np.asarray(spec) - 3.725) / 51.479, dtype=torch.float)\n",
    "                # ind = torch.LongTensor([[0]*len(ind), ind])\n",
    "                # val = torch.FloatTensor(val)\n",
    "                # spec_tensor = torch.sparse.FloatTensor(ind, val, torch.Size([1, 80000]))\n",
    "                \n",
    "                # spec_tensor = torch.tensor(np.asarray(spec), dtype=torch.float)\n",
    "                \n",
    "#                 torch.save(spec_tensor, \n",
    "#                            join(out_dir, 'spectra', '{}-{}-{}-{}-{}.pt'\n",
    "#                                 .format(tot_count, species_id, mass, l_charge, int(num_mods > 0))))\n",
    "                \n",
    "#                 np.save(join(out_dir, 'spectra', '{}.{}-{}-{}-{}-{}.npy'\n",
    "#                              .format(tot_count, sub_spec_count, species_id, mass, l_charge, num_mods)),\n",
    "#                        spec)\n",
    "                lcount += 1\n",
    "#                 spec_file_name = '{}-{}-{}.npy'.format(lcount, mass, l_charge)\n",
    "                np.save(join(out_dir, 'spectra', spec_file_name), spec)\n",
    "                # sub_spec_count > 1 means more than one spectra for the same peptide \n",
    "                # therefore, no need to save the peptide again.\n",
    "                if sub_spec_count == 1: # Write the peptide only the first time.\n",
    "                    #pep_file_name = '{}.pep'.format(tot_count)\n",
    "                    \n",
    "                    with open(join(out_dir, 'peptides', pep_file_name), 'w+') as f:\n",
    "                        f.write(pep)\n",
    "                \n",
    "                pep = 0\n",
    "                spec = []\n",
    "                new = int((i / len(lines)) * 100)\n",
    "                if new >= prev + 10:\n",
    "                    #clear_output(wait=True)\n",
    "                    print('count: ' + str(lcount))\n",
    "                    print(str(new) + '%')\n",
    "                    prev = new\n",
    "\n",
    "        #print('max peaks: ' + str(max_peaks))\n",
    "        print('In current file, read {} out of {}'.format(lcount, count))\n",
    "        print(\"Ignored: large mass: {}, pep len: {}, dup: {}\".format(mass_ign, pep_len_ign, dup_ign))\n",
    "        print('overall running count: ' + str(tot_count))\n",
    "        print('max moz: ' + str(max_moz))\n",
    "#         return pep_list, dataset, label\n",
    "#         tmp_pep_list, tmp_dataset, tmp_labels = read_msp(msp_file, species_id, decoy)\n",
    "#         pep_list.extend(tmp_dataset)\n",
    "#         dataset.extend(tmp_dataset)\n",
    "#         label.extend(tmp_labels)\n",
    "\n",
    "    # save the map. this will be used to generate masks for hard positive/negative mining during training.\n",
    "    # np.save(join(out_dir, \"idx_spec_map.npy\"), idx_spec_map)\n",
    "    with open(join(out_dir, 'pep_spec.pkl'), 'wb') as f:\n",
    "        pickle.dump(pep_spec, f)\n",
    "    \n",
    "    print(\"Statistics:\")\n",
    "    print(\"Max Missed Cleaveages: {}\".format(max_missed_cleavs))\n",
    "    print(\"Charge distribution:\")\n",
    "    print(ch)\n",
    "    print(\"Modified:\\t{}\".format(modified))\n",
    "    print(\"Unmodified:\\t{}\".format(unmodified))\n",
    "    print(\"Unique Peptides:\\t{}\".format(len(unique_pep_set)))\n",
    "    print(\"Sum: {}\".format(summ))\n",
    "    print(\"Sum-Squared: {}\".format(sq_sum))\n",
    "    print(\"N: {}\".format(N))\n",
    "    means = summ / N\n",
    "    print(\"mean: {}\".format(means))\n",
    "    stds = np.sqrt((sq_sum / N) - means**2)\n",
    "    stds[stds < 0.0000001] = float(\"inf\")\n",
    "    print(\"std: {}\".format(stds))\n",
    "    np.save(join(out_dir, 'means.npy'), means)\n",
    "    np.save(join(out_dir, 'stds.npy'), stds)\n",
    "    aux_means = aux_summ / aux_N\n",
    "    aux_stds = np.sqrt((aux_sq_sum / aux_N) - aux_means**2)\n",
    "    aux_stds[aux_stds < 0.00000001] = float(\"inf\")\n",
    "    np.save(join(out_dir, 'aux_means.npy'), aux_means)\n",
    "    np.save(join(out_dir, 'aux_stds.npy'), aux_stds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# msp_dir = \"/disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/msp-labeled/\"\n",
    "# in_tensor_dir = \"/disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/train_lstm/\"\n",
    "\n",
    "# in_dir = \"/disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/\"\n",
    "# in_tensor_dir = \"/disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/train-ready/train_lstm_mods_mass_hcd_all/\"\n",
    "\n",
    "# in_dir = \"/disk/raptor-2/mtari008/data/deepsnap/pred-full/fixed\"\n",
    "in_dir = \"/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd\"\n",
    "in_tensor_dir = \"/disk/raptor-2/mtari008/data/deepsnap/train-ready/nist_massiv_80k_ch_pep_ch-semi/\"\n",
    "#in_tensor_dir = \"/disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/train-ready/test/\"\n",
    "\n",
    "# in_dir = \"/disk/raptor-2/mtari008/data/deepsnap/pxd000612/human-hcd-phospho-mgf\"\n",
    "# in_tensor_dir = \"/disk/raptor-2/mtari008/data/deepsnap/preprocessed/human-hcd-phospho-ch-graymass\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading 11 files\n",
      "spec size: 80000\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/cptac2_mouse_hcd_itraq_phospho_selected.mgf\n",
      "len of file: 281657\n",
      "count: 115\n",
      "10%\n",
      "count: 217\n",
      "20%\n",
      "count: 337\n",
      "30%\n",
      "count: 433\n",
      "40%\n",
      "count: 541\n",
      "50%\n",
      "count: 664\n",
      "60%\n",
      "count: 783\n",
      "70%\n",
      "count: 902\n",
      "80%\n",
      "count: 1032\n",
      "90%\n",
      "In current file, read 1134 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 1108\n",
      "max moz: 42039\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_synthetic_hcd_selected.mgf\n",
      "len of file: 213341830\n",
      "count: 67736\n",
      "10%\n",
      "count: 139353\n",
      "20%\n",
      "count: 210363\n",
      "30%\n",
      "count: 281378\n",
      "40%\n",
      "count: 348647\n",
      "50%\n",
      "count: 417265\n",
      "60%\n",
      "count: 486688\n",
      "70%\n",
      "count: 556090\n",
      "80%\n",
      "count: 626310\n",
      "90%\n",
      "In current file, read 696173 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 253422\n",
      "max moz: 42039\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_hcd_tryp_good.mgf\n",
      "len of file: 35103062\n",
      "count: 19572\n",
      "10%\n",
      "count: 39483\n",
      "20%\n",
      "count: 59449\n",
      "30%\n",
      "count: 78972\n",
      "40%\n",
      "count: 98098\n",
      "50%\n",
      "count: 118067\n",
      "60%\n",
      "count: 137859\n",
      "70%\n",
      "count: 157580\n",
      "80%\n",
      "count: 177702\n",
      "90%\n",
      "In current file, read 197977 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 439965\n",
      "max moz: 48835\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/cptac2_mouse_hcd_selected.mgf\n",
      "len of file: 4170040\n",
      "count: 1654\n",
      "10%\n",
      "count: 3444\n",
      "20%\n",
      "count: 5285\n",
      "30%\n",
      "count: 7041\n",
      "40%\n",
      "count: 8896\n",
      "50%\n",
      "count: 10682\n",
      "60%\n",
      "count: 12449\n",
      "70%\n",
      "count: 14253\n",
      "80%\n",
      "count: 15898\n",
      "90%\n",
      "In current file, read 17824 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 452593\n",
      "max moz: 54618\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_hair_selected_with_GVPs_passed.mgf\n",
      "len of file: 1920729\n",
      "count: 632\n",
      "10%\n",
      "count: 1276\n",
      "20%\n",
      "count: 1895\n",
      "30%\n",
      "count: 2536\n",
      "40%\n",
      "count: 3178\n",
      "50%\n",
      "count: 3801\n",
      "60%\n",
      "count: 4442\n",
      "70%\n",
      "count: 5007\n",
      "80%\n",
      "count: 5641\n",
      "90%\n",
      "In current file, read 6273 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 456990\n",
      "max moz: 54618\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/chinese_hamster_hcd_selected.mgf\n",
      "len of file: 21770862\n",
      "count: 15867\n",
      "10%\n",
      "count: 32082\n",
      "20%\n",
      "count: 48092\n",
      "30%\n",
      "count: 63785\n",
      "40%\n",
      "count: 79489\n",
      "50%\n",
      "count: 95308\n",
      "60%\n",
      "count: 110847\n",
      "70%\n",
      "count: 126257\n",
      "80%\n",
      "count: 141948\n",
      "90%\n",
      "In current file, read 157726 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 253422\n",
      "max moz: 54777\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/cptac2_human_hcd_itraq_phospho_selected.mgf\n",
      "len of file: 3483183\n",
      "count: 1263\n",
      "10%\n",
      "count: 2569\n",
      "20%\n",
      "count: 3837\n",
      "30%\n",
      "count: 5116\n",
      "40%\n",
      "count: 6422\n",
      "50%\n",
      "count: 7759\n",
      "60%\n",
      "count: 9156\n",
      "70%\n",
      "count: 10579\n",
      "80%\n",
      "count: 12003\n",
      "90%\n",
      "In current file, read 13322 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 560235\n",
      "max moz: 54777\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_hcd_tryp_best.mgf\n",
      "len of file: 53157094\n",
      "count: 39665\n",
      "10%\n",
      "count: 79384\n",
      "20%\n",
      "count: 119536\n",
      "30%\n",
      "count: 158967\n",
      "40%\n",
      "count: 198145\n",
      "50%\n",
      "count: 238150\n",
      "60%\n",
      "count: 277319\n",
      "70%\n",
      "count: 316416\n",
      "80%\n",
      "count: 356384\n",
      "90%\n",
      "In current file, read 397196 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 846442\n",
      "max moz: 54777\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/massiv-all.mgf\n",
      "len of file: 250114903\n",
      "count: 211161\n",
      "10%\n",
      "count: 421844\n",
      "20%\n",
      "count: 632613\n",
      "30%\n",
      "count: 844441\n",
      "40%\n",
      "count: 1055390\n",
      "50%\n",
      "count: 1266515\n",
      "60%\n",
      "count: 1477978\n",
      "70%\n",
      "count: 1689185\n",
      "80%\n",
      "count: 1900541\n",
      "90%\n",
      "In current file, read 2111350 out of 3\n",
      "Ignored: large mass: 0, pep len: 556, dup: 0\n",
      "overall running count: 2360013\n",
      "max moz: 54777\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/rat_qtof_consensus_final_true_lib.mgf\n",
      "len of file: 1551941\n",
      "count: 2316\n",
      "10%\n",
      "count: 4632\n",
      "20%\n",
      "count: 6839\n",
      "30%\n",
      "count: 9099\n",
      "40%\n",
      "count: 11408\n",
      "50%\n",
      "count: 13638\n",
      "60%\n",
      "count: 15906\n",
      "70%\n",
      "count: 18301\n",
      "80%\n",
      "count: 20536\n",
      "90%\n",
      "In current file, read 22737 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 452593\n",
      "max moz: 54777\n",
      "Reading: /lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_hcd_semitryp.mgf\n",
      "len of file: 54407007\n",
      "count: 30160\n",
      "10%\n",
      "count: 60115\n",
      "20%\n",
      "count: 91360\n",
      "30%\n",
      "count: 123714\n",
      "40%\n",
      "count: 154331\n",
      "50%\n",
      "count: 186443\n",
      "60%\n",
      "count: 217362\n",
      "70%\n",
      "count: 247855\n",
      "80%\n",
      "count: 278746\n",
      "90%\n",
      "In current file, read 311604 out of 3\n",
      "Ignored: large mass: 0, pep len: 0, dup: 0\n",
      "overall running count: 2518072\n",
      "max moz: 54777\n",
      "Statistics:\n",
      "Max Missed Cleaveages: 0\n",
      "Charge distribution:\n",
      "[      0.   18813. 1977540. 1428987.  433514.   98244.       0.       0.\n",
      "       0.       0.       0.       0.       0.       0.       0.       0.\n",
      "       0.       0.       0.       0.]\n",
      "Modified:\t1459405\n",
      "Unmodified:\t2497693\n",
      "Unique Peptides:\t2538352\n",
      "Sum: [3933316. 3914722. 1951191. ...       0.       0.       0.]\n",
      "Sum-Squared: [3933316. 3914722. 1951191. ...       0.       0.       0.]\n",
      "N: 3933316\n",
      "mean: [1.         0.99527269 0.49606769 ... 0.         0.         0.        ]\n",
      "std: [       inf 0.06859272 0.49998454 ...        inf        inf        inf]\n"
     ]
    }
   ],
   "source": [
    "# transformer = transforms.Normalize(mean=[3.725], std=[51.479])\n",
    "preprocess_mgfs(in_dir, in_tensor_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocess_mgf_mods(mgf_dir):\n",
    "    \n",
    "    mgf_files = verify_in_dir(mgf_dir, \"mgf\")\n",
    "        \n",
    "    print('reading {} files'.format(len(mgf_files)))\n",
    "    \n",
    "    mod_counters = np.zeros(20)\n",
    "    mod_type_counter_dict = {}\n",
    "    n_mods = []\n",
    "    i_mods = []\n",
    "    c_mods = []\n",
    "    modified = unmodified = 0\n",
    "    \n",
    "    mass_ign = pep_len_ign = 0\n",
    "    \n",
    "    for species_id, mgf_file in enumerate(mgf_files):\n",
    "        print('Reading: {}'.format(mgf_file))\n",
    "        \n",
    "        f = open(mgf_file, \"r\")\n",
    "        lines = f.readlines()\n",
    "        f.close()\n",
    "        \n",
    "        count = lcount = 0\n",
    "\n",
    "        # FIXME: config should use only one get_config call.\n",
    "        seq_len = config.get_config(section='ml', key='pep_seq_len')\n",
    "        spec_size = config.get_config(section='input', key='spec_size')\n",
    "\n",
    "        print('len of file: ' + str(len(lines)))\n",
    "        prev = 0\n",
    "        i = 0\n",
    "        while i < len(lines):\n",
    "            line = lines[i]\n",
    "            i += 1\n",
    "\n",
    "            if line.startswith('PEPMASS'):\n",
    "                count += 1\n",
    "                mass = float(re.findall(r\"PEPMASS=([-+]?[0-9]*\\.?[0-9]*)\", line)[0])\n",
    "                if not round(mass)*10 < spec_size:\n",
    "                    mass_ign += 1\n",
    "                    continue\n",
    "                \n",
    "            if line.startswith(\"SEQ\"):\n",
    "                line = re.sub(r\"[()]\", \"\", line.strip())\n",
    "                pep_segs = re.findall(r\"([A-Z]+)\", line)[1:]\n",
    "                pep = ''.join(pep_segs)\n",
    "                #mods = re.findall(r\"([A-Z]\\([^\\(\\)]*\\))\", line)\n",
    "                \n",
    "                \n",
    "                rex = r\"([A-Z]?)([-+]?\\d*\\.\\d+|\\d+)([A-Z]?)\"\n",
    "                mods = re.findall(rex, line)\n",
    "                num_mods = len(mods)\n",
    "                \n",
    "                if len(pep) + num_mods > seq_len:\n",
    "                    pep_len_ign += 1\n",
    "                    continue\n",
    "                \n",
    "                if num_mods <= len(mod_counters):\n",
    "                    mod_counters[num_mods] += 1\n",
    "                    \n",
    "                if num_mods > 0:\n",
    "                    modified += 1\n",
    "                else:\n",
    "                    unmodified += 1\n",
    "                    \n",
    "                for a in mods:\n",
    "                    mass = round(float(a[1]), 2)\n",
    "                    pre = a[0]\n",
    "                    post = a[2]\n",
    "                    if mass in mod_type_counter_dict:\n",
    "                        mod_type_counter_dict[mass][\"count\"] += 1\n",
    "                        if pre in mod_type_counter_dict[mass]:\n",
    "                            mod_type_counter_dict[mass][pre] += 1\n",
    "                        else:\n",
    "                            mod_type_counter_dict[mass][pre] = 1\n",
    "                    else:\n",
    "                        mod_type_counter_dict[mass] = {\"count\":1, pre:1}\n",
    "                    \n",
    "                    if not a[0]:\n",
    "                        n_mods.append(a)\n",
    "                    elif not a[2]:\n",
    "                        c_mods.append(a)\n",
    "                    elif a[0] and a[1] and a[2]:\n",
    "                        i_mods.append(a)\n",
    "                    else:\n",
    "                        print(\"Invalid mod\")\n",
    "                \n",
    "                lcount += 1\n",
    "                \n",
    "                pep = 0\n",
    "                spec = []\n",
    "                new = int((i / len(lines)) * 100)\n",
    "                if new >= prev + 10:\n",
    "                    #clear_output(wait=True)\n",
    "                    print('count: ' + str(lcount))\n",
    "                    print(str(new) + '%')\n",
    "                    prev = new\n",
    "\n",
    "        #print('max peaks: ' + str(max_peaks))\n",
    "        print('In current file, read {} out of {}'.format(lcount, count))\n",
    "    \n",
    "    print(\"Statistics:\")\n",
    "    print(\"Modified:\\t{}\".format(modified))\n",
    "    print(\"Unmodified:\\t{}\".format(unmodified))\n",
    "    print(\"Num mod distribution:\")\n",
    "    print(mod_counters)\n",
    "    print(\"Mod types:\")\n",
    "    print(\"{} types of mods found\".format(len(mod_type_counter_dict)))\n",
    "    print(\"Distribution:\")\n",
    "    print(mod_type_counter_dict)\n",
    "    print(\"Number of N mods: {}\".format(len(n_mods)))\n",
    "    print(\"Number of C mods: {}\".format(len(c_mods)))\n",
    "    print(\"Number of I mods: {}\".format(len(i_mods)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading 10 files\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/cptac2_mouse_hcd_itraq_phospho_selected.mgf\n",
      "len of file: 281657\n",
      "count: 116\n",
      "10%\n",
      "count: 218\n",
      "20%\n",
      "count: 338\n",
      "30%\n",
      "count: 434\n",
      "40%\n",
      "count: 542\n",
      "50%\n",
      "count: 665\n",
      "60%\n",
      "count: 784\n",
      "70%\n",
      "count: 903\n",
      "80%\n",
      "count: 1033\n",
      "90%\n",
      "In current file, read 1134 out of 1134\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_synthetic_hcd_selected.mgf\n",
      "len of file: 213341830\n",
      "count: 67748\n",
      "10%\n",
      "count: 139377\n",
      "20%\n",
      "count: 210414\n",
      "30%\n",
      "count: 281468\n",
      "40%\n",
      "count: 348745\n",
      "50%\n",
      "count: 417377\n",
      "60%\n",
      "count: 486829\n",
      "70%\n",
      "count: 556235\n",
      "80%\n",
      "count: 626455\n",
      "90%\n",
      "In current file, read 696341 out of 696341\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/cptac2_mouse_hcd_selected.mgf\n",
      "len of file: 4170040\n",
      "count: 1656\n",
      "10%\n",
      "count: 3446\n",
      "20%\n",
      "count: 5287\n",
      "30%\n",
      "count: 7043\n",
      "40%\n",
      "count: 8900\n",
      "50%\n",
      "count: 10688\n",
      "60%\n",
      "count: 12455\n",
      "70%\n",
      "count: 14259\n",
      "80%\n",
      "count: 15905\n",
      "90%\n",
      "In current file, read 17830 out of 17830\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_hair_selected_with_GVPs_passed.mgf\n",
      "len of file: 1920729\n",
      "count: 633\n",
      "10%\n",
      "count: 1277\n",
      "20%\n",
      "count: 1896\n",
      "30%\n",
      "count: 2537\n",
      "40%\n",
      "count: 3179\n",
      "50%\n",
      "count: 3802\n",
      "60%\n",
      "count: 4443\n",
      "70%\n",
      "count: 5012\n",
      "80%\n",
      "count: 5646\n",
      "90%\n",
      "In current file, read 6277 out of 6277\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/chinese_hamster_hcd_selected.mgf\n",
      "len of file: 21770862\n",
      "count: 15881\n",
      "10%\n",
      "count: 32100\n",
      "20%\n",
      "count: 48111\n",
      "30%\n",
      "count: 63824\n",
      "40%\n",
      "count: 79540\n",
      "50%\n",
      "count: 95373\n",
      "60%\n",
      "count: 110923\n",
      "70%\n",
      "count: 126347\n",
      "80%\n",
      "count: 142054\n",
      "90%\n",
      "In current file, read 157838 out of 157838\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/cptac2_human_hcd_itraq_phospho_selected.mgf\n",
      "len of file: 3483183\n",
      "count: 1264\n",
      "10%\n",
      "count: 2570\n",
      "20%\n",
      "count: 3838\n",
      "30%\n",
      "count: 5117\n",
      "40%\n",
      "count: 6423\n",
      "50%\n",
      "count: 7760\n",
      "60%\n",
      "count: 9157\n",
      "70%\n",
      "count: 10580\n",
      "80%\n",
      "count: 12004\n",
      "90%\n",
      "In current file, read 13322 out of 13322\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_hcd_tryp_best.mgf\n",
      "len of file: 53223140\n",
      "count: 39778\n",
      "10%\n",
      "count: 79604\n",
      "20%\n",
      "count: 119863\n",
      "30%\n",
      "count: 159438\n",
      "40%\n",
      "count: 198762\n",
      "50%\n",
      "count: 238889\n",
      "60%\n",
      "count: 278171\n",
      "70%\n",
      "count: 317402\n",
      "80%\n",
      "count: 357467\n",
      "90%\n",
      "In current file, read 398373 out of 398373\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/human_hcd_labelfree_phospho_selected_passed.mgf\n",
      "len of file: 22597575\n",
      "count: 6613\n",
      "10%\n",
      "count: 13083\n",
      "20%\n",
      "count: 19699\n",
      "30%\n",
      "count: 26154\n",
      "40%\n",
      "count: 32627\n",
      "50%\n",
      "count: 39208\n",
      "60%\n",
      "count: 45495\n",
      "70%\n",
      "count: 53141\n",
      "80%\n",
      "count: 60262\n",
      "90%\n",
      "In current file, read 66922 out of 66922\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/massiv-all.mgf\n",
      "len of file: 250114903\n",
      "count: 215408\n",
      "10%\n",
      "count: 430302\n",
      "20%\n",
      "count: 645472\n",
      "30%\n",
      "count: 861629\n",
      "40%\n",
      "count: 1076846\n",
      "50%\n",
      "count: 1292171\n",
      "60%\n",
      "count: 1507968\n",
      "70%\n",
      "count: 1723442\n",
      "80%\n",
      "count: 1939125\n",
      "90%\n",
      "In current file, read 2154269 out of 2154269\n",
      "Reading: /disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/rat_qtof_consensus_final_true_lib.mgf\n",
      "len of file: 1551941\n",
      "count: 2383\n",
      "10%\n",
      "count: 4781\n",
      "20%\n",
      "count: 7045\n",
      "30%\n",
      "count: 9377\n",
      "40%\n",
      "count: 11757\n",
      "50%\n",
      "count: 14067\n",
      "60%\n",
      "count: 16395\n",
      "70%\n",
      "count: 18894\n",
      "80%\n",
      "count: 21219\n",
      "90%\n",
      "In current file, read 23492 out of 23492\n",
      "Statistics:\n",
      "Modified:\t1402486\n",
      "Unmodified:\t2133312\n",
      "Num mod distribution:\n",
      "[2.133312e+06 1.112510e+06 2.283820e+05 4.711300e+04 9.696000e+03\n",
      " 2.861000e+03 1.164000e+03 4.650000e+02 2.000000e+02 7.100000e+01\n",
      " 1.900000e+01 5.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n",
      " 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00]\n",
      "Mod types:\n",
      "11 types of mods found\n",
      "Distribution:\n",
      "{79.97: {'count': 99632, 'S': 71756, 'T': 15925, 'Y': 2386, '': 9565}, 15.99: {'count': 395876, 'M': 352966, '': 42910}, 57.02: {'count': 743391, 'C': 701480, '': 41911}, -18.01: {'count': 688, '': 688}, -17.03: {'count': 40313, '': 40313}, 42.01: {'count': 88495, '': 88495}, 26.02: {'count': 284, '': 284}, 0.01: {'count': 1070, '': 1070}, 0.02: {'count': 12, '': 12}, 43.01: {'count': 123561, '': 123561}, 0.98: {'count': 283103, '': 5531, 'Q': 131437, 'N': 146135}}\n",
      "Number of N mods: 354340\n",
      "Number of C mods: 3980\n",
      "Number of I mods: 1418105\n"
     ]
    }
   ],
   "source": [
    "in_dir = \"/disk/raptor/lclhome/mtari008/Proteogenomics/DeepSNAP/data/mgfs-hcd/\"\n",
    "preprocess_mgf_mods(in_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistics  \n",
    "#### Modified:\t1259406  \n",
    "#### Unmodified:\t2422940  \n",
    "### Num mod distribution:  \n",
    "[2.422940e+06 1.023623e+06 1.836600e+05 3.832100e+04 9.358000e+03  \n",
    " 3.060000e+03 8.790000e+02 3.250000e+02 1.290000e+02 4.000000e+01  \n",
    " 8.000000e+00 3.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00  \n",
    " 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00]  \n",
    "### Mod types:  \n",
    "6 types of mods found  \n",
    "### Distribution:  \n",
    "{15.99: 290045, 0.98: 396500, 57.02: 640177, 42.01: 82730, 43.01: 123561, -17.03: 34679}  \n",
    "{15.99: {'count': 290045, 'M': 279077, '': 10968}, 0.98: {'count': 396500, 'N': 234278, 'Q': 153930, '': 8292}, 57.02: {'count': 640177, 'C': 615747, '': 24430}, 42.01: {'count': 82730, '': 82730}, 43.01: {'count': 123561, '': 123561}, -17.03: {'count': 34679, '': 34679}}  \n",
    "{15.99: {'count': 344652, 'M': 326428, '': 18224}, 0.98: {'count': 396500, 'N': 234278, 'Q': 153930, '': 8292}, 57.02: {'count': 745790, 'C': 711468, '': 34322}, 42.01: {'count': 86484, '': 86484}, 0.01: {'count': 702, '': 702}, -18.01: {'count': 18, '': 18}, -17.03: {'count': 34756, '': 34756}, 43.01: {'count': 123561, '': 123561}}  \n",
    "#### Number of N mods: 284660  \n",
    "#### Number of C mods: 4734  \n",
    "#### Number of I mods: 1278298  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = 325447965576.6338 / 25019920000\n",
    "print(\"mean: {}\".format(mean))\n",
    "std = math.sqrt((2885406787428155.0 / 25019920000) - mean**2)\n",
    "print(std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LabeledSpectra(data.Dataset):\n",
    "    'Characterizes a dataset for PyTorch'\n",
    "    def __init__(self, dir_path, filt, test=False):\n",
    "        'Initialization'\n",
    "        \n",
    "        self.aas         = ['_PAD'] + list(config.AAMass.keys())\n",
    "        self.aa2idx      = {a:i for i, a in enumerate(self.aas)}\n",
    "        self.idx2aa      = {i:a for i, a in enumerate(self.aas)}\n",
    "        \n",
    "        self.spec_path   = join(dir_path, 'spectra')\n",
    "        self.pep_path    = join(dir_path, 'peptides')\n",
    "        self.charge      = filt['charge'] if 'charge' in filt else config.get_config(section='input', key='charge')\n",
    "        self.num_species = config.get_config(section='input', key='num_species')\n",
    "        # self.vocab_size  = len(self.aa2idx) + self.charge + self.num_species + 1\n",
    "        self.vocab_size  = round(max(config.AAMass.values())) + 1\n",
    "        self.seq_len     = config.get_config(section='ml', key='pep_seq_len')\n",
    "        self.modified    = filt['modified'] if 'modified' in filt else False\n",
    "        self.test_size   = config.get_config(section='ml', key='test_size')\n",
    "        self.test        = test\n",
    "        \n",
    "        self.file_names  = []\n",
    "        for file in listdir(self.spec_path):\n",
    "            if self.apply_filter(file):\n",
    "                self.file_names.append(file)\n",
    "        \n",
    "        print('dataset size: {}'.format(len(self.file_names)))        \n",
    "        \n",
    "        self.train_files, self.test_files = train_test_split(\n",
    "            self.file_names, test_size = self.test_size, random_state = rand.randint(0, 1000), shuffle=True)\n",
    "        \n",
    "        if self.test:\n",
    "            print('test size: {}'.format(len(self.test_files)))\n",
    "        else:\n",
    "            print('train size: {}'.format(len(self.train_files)))\n",
    "        \n",
    "    def __len__(self):\n",
    "        'Denotes the total number of samples'\n",
    "        if self.test:\n",
    "            return len(self.test_files)\n",
    "        else:\n",
    "            return len(self.train_files)\n",
    "\n",
    "    def __getitem__(self, index):\n",
    "        'Generates one sample of data'\n",
    "        file_name = ''\n",
    "        # Select sample\n",
    "        if self.test:\n",
    "            file_name = self.test_files[index]\n",
    "        else:\n",
    "            file_name = self.train_files[index]\n",
    "            \n",
    "        spec_file_name = join(self.spec_path, file_name)\n",
    "        pep_file_name  = join(self.pep_path, file_name.replace('.pt', '.pep'))\n",
    "        \n",
    "        # Load data and get label\n",
    "        spec_torch = torch.load(spec_file_name)\n",
    "        \n",
    "        # Load peptide and convert to idx array\n",
    "        f = open(pep_file_name, \"r\")\n",
    "        pep = f.readlines()[0].strip()\n",
    "        f.close()\n",
    "        \n",
    "        pepl = np.zeros(len(pep))\n",
    "        file_parts = re.search(r\"(\\d+)-(\\d+)-(\\d+.\\d+)-(\\d)-(0|1).pt\", file_name)\n",
    "        pepl[0] = int(file_parts[4]) + len(self.aas)  # coded value of charge\n",
    "        pepl[1] = int(file_parts[2]) + self.charge + 1 + len(self.aas) # coded value of specie id\n",
    "        \n",
    "        # for i in range(2, len(pep)):\n",
    "        #     pepl[i] = self.aa2idx[pep[i]]\n",
    "        for i, aa in enumerate(pep[2:]):\n",
    "            pepl[i + 2] = self.aa2idx[aa]\n",
    "            # pepl[i + 2] = round(config.AAMass[aa])\n",
    "        \n",
    "        pepl = self.pad_left(pepl, self.seq_len)\n",
    "        pep_torch = torch.tensor(pepl, dtype=torch.long)\n",
    "        \n",
    "        return [spec_torch, pep_torch]\n",
    "    \n",
    "    def apply_filter(self, file_name):\n",
    "        file_parts = re.search(r\"(\\d+)-(\\d+)-(\\d+.\\d+)-(\\d)-(0|1).pt\", file_name)\n",
    "        charge = int(file_parts[4])\n",
    "        modified = bool(int(file_parts[5]))\n",
    "        \n",
    "        if ((self.charge == 0 or charge <= self.charge)\n",
    "            and (self.modified or self.modified == modified)):\n",
    "            return True\n",
    "        \n",
    "        return False\n",
    "    \n",
    "    def pad_left(self, arr, size):\n",
    "        out = np.zeros(size)\n",
    "        out[-len(arr):] = arr\n",
    "        return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_mean = np.mean(dataset)\n",
    "# data_std = np.std(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# batch_size = config.get_config(section='ml', key='batch_size')\n",
    "batch_size = 256\n",
    "charge = config.get_config(section='input', key='charge')\n",
    "use_mods = config.get_config(section='input', key='use_mods')\n",
    "filt = {'charge':charge, 'modified':use_mods}\n",
    "\n",
    "\n",
    "train_dataset = LabeledSpectra(in_tensor_dir, filt, test=False)\n",
    "test_dataset = LabeledSpectra(in_tensor_dir, filt, test=True)\n",
    "\n",
    "train_loader = torch.utils.data.DataLoader(\n",
    "    dataset=train_dataset, batch_size=batch_size, shuffle=True, drop_last=True, num_workers=8)\n",
    "\n",
    "test_loader = torch.utils.data.DataLoader(\n",
    "    dataset=test_dataset, batch_size=batch_size, shuffle=False, drop_last=True, num_workers=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "do_learn = True\n",
    "save_frequency = 2\n",
    "lr = 0.0001\n",
    "num_epochs = 200\n",
    "weight_decay = 0.0001\n",
    "margin = 0.2\n",
    "vocab_size = train_dataset.vocab_size\n",
    "#torch.manual_seed(0)\n",
    "#torch.cuda.manual_seed(0)\n",
    "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
    "#device = torch.device('cuda')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#train_loader = torch.load('train_loader.pt')\n",
    "#test_loader = torch.load('test_loader.pt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Net(nn.Module):\n",
    "    def __init__(self, vocab_size, output_size=512, embedding_dim=512, hidden_lstm_dim=1024, lstm_layers=2, drop_prob=0.5):\n",
    "        super(Net, self).__init__()\n",
    "        \n",
    "        self.spec_size = config.get_config(section='input', key='spec_size')\n",
    "        self.output_size = output_size\n",
    "        self.lstm_layers = lstm_layers\n",
    "        self.hidden_lstm_dim = hidden_lstm_dim\n",
    "        self.embedding_dim = embedding_dim\n",
    "        self.searching = False\n",
    "        \n",
    "        self.embedding = nn.Embedding(vocab_size, embedding_dim)\n",
    "        self.lstm = nn.LSTM(embedding_dim, self.hidden_lstm_dim, self.lstm_layers, \n",
    "                            dropout=drop_prob, batch_first=True, bidirectional=True)\n",
    "        \n",
    "        self.linear1_1 = nn.Linear(self.spec_size, 1024)\n",
    "        self.linear1_2 = nn.Linear(1024, 512)\n",
    "        self.linear1_3 = nn.Linear(512, 256)\n",
    "        \n",
    "        self.linear2_1 = nn.Linear(2048, 1024)\n",
    "        self.linear2_2 = nn.Linear(1024, 512)\n",
    "        self.linear2_3 = nn.Linear(512, 256)\n",
    "        \n",
    "        self.dropout1 = nn.Dropout(0.5)\n",
    "        self.dropout2 = nn.Dropout(0.2)\n",
    "        #self.dropout3 = nn.Dropout(0.3)\n",
    "        \n",
    "    def forward(self, data, hidden):\n",
    "        specs = data[0]\n",
    "        peps = data[1]\n",
    "        # print(peps.type())\n",
    "        # print('Input to the model size: {}'.format(specs.size()))\n",
    "        # print('Input to the model size: {}'.format(peps.size()))\n",
    "        # peps = peps.unsqueeze(-1).float()\n",
    "        # print(peps.type())\n",
    "        \n",
    "        embeds = self.embedding(peps)\n",
    "        lstm_out, hidden = self.lstm(embeds, hidden)\n",
    "        # print(lstm_out.size())\n",
    "        lstm_out = lstm_out[:, -1, :]\n",
    "        # print(lstm_out.size())\n",
    "        #lstm_out = torch.mean(lstm_out, dim=1)\n",
    "        lstm_out = lstm_out.contiguous().view(-1, self.hidden_lstm_dim * 2)\n",
    "        out = self.dropout2(lstm_out)\n",
    "        \n",
    "        out = self.linear2_1(out)\n",
    "        out = F.relu(out)\n",
    "        out = self.dropout2(out)\n",
    "        \n",
    "        out = self.linear2_2(out)\n",
    "        out = F.relu(out)\n",
    "        \n",
    "        # out = out.view(batch_size, peps.size()[1], 512)\n",
    "        # out = out[:,-1,:]\n",
    "        out_pep = F.normalize(out)\n",
    "        \n",
    "        out = self.linear1_1(specs.view(-1, self.spec_size))\n",
    "        out = F.relu(out)\n",
    "        out = self.dropout2(out)\n",
    "        out = self.linear1_2(out)\n",
    "        out = F.relu(out)\n",
    "        out_spec = F.normalize(out)\n",
    "        \n",
    "        res = out_spec, out_pep, hidden\n",
    "        return res\n",
    "    \n",
    "    def init_hidden(self, batch_size):\n",
    "        weight = next(self.parameters()).data\n",
    "        hidden = (weight.new(self.lstm_layers * 2, batch_size, self.hidden_lstm_dim).zero_().to(device),\n",
    "                      weight.new(self.lstm_layers * 2, batch_size, self.hidden_lstm_dim).zero_().to(device))\n",
    "        return hidden\n",
    "    \n",
    "    def name(self):\n",
    "        return \"Net\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hinge = torch.nn.HingeEmbeddingLoss()\n",
    "triplet_loss = nn.TripletMarginLoss(margin=margin, p=2, reduction='sum')\n",
    "l2_squared = nn.MSELoss(reduction='none')\n",
    "pdist = nn.PairwiseDistance(p=2)\n",
    "zero_tensor = torch.tensor(0.).to(device)\n",
    "clip = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(model, device, train_loader, epoch, optimizer):\n",
    "    model.train()\n",
    "    h = model.init_hidden(batch_size)\n",
    "    \n",
    "    accurate_labels = 0\n",
    "    all_labels = 0\n",
    "    for (batch_idx, data) in enumerate(train_loader):\n",
    "        h = tuple([e.data for e in h])\n",
    "        data[0], data[1] = data[0].to(device), data[1].to(device)\n",
    "        \n",
    "        optimizer.zero_grad()\n",
    "        \n",
    "        Q, P, h = model(data, h)\n",
    "        nn.utils.clip_grad_norm_(model.parameters(), clip)\n",
    "        \"\"\"Mine the hardest triplets. Get rid of N.\"\"\" \n",
    "        QxQ = process.pairwise_distances(Q)    # calculate distance matrix for spectra\n",
    "        PxP = process.pairwise_distances(P)    # calculate distance matrix for peptides\n",
    "        QxP_ = process.pairwise_distances(Q, P) # calculate distance matrix for spectra-peptides\n",
    "        \n",
    "        # Set the diagonal of all distance matrices to inf so we don't get self as the closest negative.\n",
    "        QxQ.fill_diagonal_(float(\"inf\"))\n",
    "        PxP.fill_diagonal_(float(\"inf\"))\n",
    "        QxP = QxP_.clone()\n",
    "        QxP.fill_diagonal_(float(\"inf\"))\n",
    "        \n",
    "        #print(QP.argmin(1)[:100])\n",
    "        \n",
    "#         pos = torch.sum(l2_squared(Q, P), dim=1) + margin\n",
    "        \n",
    "#         QxQ_min = QxQ.gather(1, torch.randint(len(Q), (len(Q),), device=device).view(-1,1))             # farthest spectrum for each spectrum\n",
    "#         PxP_min = PxP.gather(1, torch.randint(len(Q), (len(Q),), device=device).view(-1,1))             # farthest peptide for each peptide\n",
    "#         QxP_min = QxP.gather(1, torch.randint(len(Q), (len(Q),), device=device).view(-1,1))             # farthest peptide for each spectrum\n",
    "#         PxQ_min = QxP.gather(0, torch.randint(len(Q), (len(Q),), device=device).view(1,-1))             # farthest spectrum for each peptide\n",
    "        \n",
    "        QxQ_min = QxQ.min(1).values              # nearest spectrum for each spectrum\n",
    "        PxP_min = PxP.min(1).values              # nearest peptide for each peptide\n",
    "        QxP_min = QxP.min(1).values              # nearest peptide for each spectrum\n",
    "        PxQ_min = QxP.min(0).values              # nearest spectrum for each peptide\n",
    "        \n",
    "        #neg = QxQ_min + PxP_min + QxP_min + PxQ_min\n",
    "        \n",
    "#         divider = torch.tensor(float(len(pos)))\n",
    "#         loss = torch.sum(torch.max(pos - QxQ_min, zero_tensor)) / divider\n",
    "#         loss += torch.sum(torch.max(pos - PxP_min, zero_tensor)) / divider\n",
    "#         loss += torch.sum(torch.max(pos - QxP_min, zero_tensor)) / divider\n",
    "#         loss += torch.sum(torch.max(pos - PxQ_min, zero_tensor)) / divider\n",
    "        \n",
    "        #divider = torch.sum(pos - neg > 0)\n",
    "        #loss = torch.sum(torch.max(pos - neg, zero_tensor)) / divider\n",
    "        \n",
    "        QxQ_min = Q[QxQ.min(1).indices]              # nearest spectrum for each spectrum\n",
    "        PxP_min = P[PxP.min(1).indices]              # nearest peptide for each peptide\n",
    "        QxP_min = P[QxP.min(1).indices]              # nearest peptide for each spectrum\n",
    "        PxQ_min = Q[QxP.min(0).indices]              # nearest spectrum for each peptide\n",
    "        loss = triplet_loss(Q, P, QxQ_min)           # spectrum-spectrum negatives\n",
    "        loss += triplet_loss(Q, P, QxP_min)          # spectrum-peptide negatives\n",
    "        loss += triplet_loss(P, Q, PxP_min)          # peptide-peptide negatives\n",
    "        loss += triplet_loss(P, Q, PxQ_min)          # peptide-spectrum negatives\n",
    "        \n",
    "        loss = loss / 4\n",
    "                \n",
    "        loss.backward()\n",
    "            \n",
    "        optimizer.step()\n",
    "        \n",
    "        seq = torch.arange(0, len(Q), step=1, device=device, requires_grad=False)\n",
    "        accurate_labels = accurate_labels + torch.sum(QxP_.argmin(1) == seq) # use QP_ since it doesn't have diag set to zero\n",
    "        \n",
    "        all_labels = all_labels + len(Q)  \n",
    "    \n",
    "    accuracy = 100. * float(accurate_labels) / all_labels\n",
    "    train_accuracy.append(accuracy)\n",
    "    train_loss.append(loss)\n",
    "    print('Epoch: ' + str(epoch))\n",
    "    print('Train accuracy: {}/{} ({:.3f}%)\\tLoss: {:.6f}'.format(accurate_labels, all_labels, accuracy, loss))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test(model, device, test_loader):\n",
    "    model.eval()\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        accurate_labels = 0\n",
    "        all_labels = 0\n",
    "        loss = 0\n",
    "        h = model.init_hidden(batch_size)\n",
    "        \n",
    "        for (batch_idx, data) in enumerate(test_loader):\n",
    "            h = tuple([e.data for e in h])\n",
    "            data[0], data[1] = data[0].to(device), data[1].to(device)\n",
    "        \n",
    "            optimizer.zero_grad()\n",
    "            \n",
    "            Q, P, h = model(data, h)\n",
    "            \n",
    "            \"\"\"Mine the hardest triplets. Get rid of N.\"\"\" \n",
    "            QxQ = process.pairwise_distances(Q)    # calculate distance matrix for spectra\n",
    "            PxP = process.pairwise_distances(P)    # calculate distance matrix for peptides\n",
    "            QxP_ = process.pairwise_distances(Q, P) # calculate distance matrix for spectra-peptides\n",
    "\n",
    "            # Set the diagonal of all distance matrices to inf so we don't get self as the closest negative.\n",
    "            QxQ.fill_diagonal_(float(\"inf\"))\n",
    "            PxP.fill_diagonal_(float(\"inf\"))\n",
    "            QxP = QxP_.clone()    # clone to measure accuracy. can be done in a better way.\n",
    "            QxP.fill_diagonal_(float(\"inf\"))\n",
    "\n",
    "            #print(QP.argmin(1)[:100])\n",
    "\n",
    "            pos = 4 * (torch.sum(l2_squared(Q, P), dim=1) + margin)\n",
    "\n",
    "            QxQ_min = QxQ.min(1).values              # farthest spectrum for each spectrum\n",
    "            PxP_min = PxP.min(1).values              # farthest peptide for each peptide\n",
    "            QxP_min = QxP.min(1).values              # farthest peptide for each spectrum\n",
    "            PxQ_min = QxP.min(0).values              # farthest spectrum for each peptide\n",
    "\n",
    "            #neg = QxQ_min + PxP_min + QxP_min + PxQ_min\n",
    "        \n",
    "#             divider = torch.tensor(float(len(pos)))\n",
    "#             loss = torch.sum(torch.max(pos - QxQ_min, zero_tensor)) / divider\n",
    "#             loss += torch.sum(torch.max(pos - PxP_min, zero_tensor)) / divider\n",
    "#             loss += torch.sum(torch.max(pos - QxP_min, zero_tensor)) / divider\n",
    "#             loss += torch.sum(torch.max(pos - PxQ_min, zero_tensor)) / divider\n",
    "\n",
    "            QxQ_min = Q[QxQ.min(1).indices]              # nearest spectrum for each spectrum\n",
    "            PxP_min = P[PxP.min(1).indices]              # nearest peptide for each peptide\n",
    "            QxP_min = P[QxP.min(1).indices]              # nearest peptide for each spectrum\n",
    "            PxQ_min = Q[QxP.min(0).indices]              # nearest spectrum for each peptide\n",
    "            loss = triplet_loss(Q, P, QxQ_min)     # spectrum-spectrum negatives\n",
    "            loss += triplet_loss(Q, P, QxP_min)    # spectrum-peptide negatives\n",
    "            loss += triplet_loss(P, Q, PxP_min)    # peptide-peptide negatives\n",
    "            loss += triplet_loss(P, Q, PxQ_min)    # peptide-spectrum negatives\n",
    "            \n",
    "            loss = loss / 4\n",
    "\n",
    "            #divider = torch.tensor(float(len(pos)))\n",
    "            #divider = torch.sum(pos - neg > 0)\n",
    "            #loss = torch.sum(torch.max(pos - neg, zero_tensor)) / divider\n",
    "            \n",
    "#             loss =  torch.sum(torch.max(pos - QxQ_min, zero_tensor)) / divider\n",
    "#             loss += torch.sum(torch.max(pos - PxP_min, zero_tensor)) / divider\n",
    "#             loss += torch.sum(torch.max(pos - QxP_min, zero_tensor)) / divider\n",
    "#             loss += torch.sum(torch.max(pos - PxQ_min, zero_tensor)) / divider\n",
    "            \n",
    "            seq = torch.arange(0, len(Q), step=1, device=device, requires_grad=False)\n",
    "            accurate_labels = accurate_labels + torch.sum(QxP_.argmin(1) == seq) # use QP_ since it doesn't have diag set to zero\n",
    "            \n",
    "            all_labels = all_labels + len(Q)\n",
    "                \n",
    "        accuracy = 100. * float(accurate_labels) / all_labels\n",
    "        test_accuracy.append(accuracy)\n",
    "        test_loss.append(loss)\n",
    "        print('Test accuracy: {}/{} ({:.3f}%)\\tLoss: {:.6f}'.format(accurate_labels, all_labels, accuracy, loss))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def oneshot(model, device, data):\n",
    "    model.eval()\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        for i in range(len(data)):\n",
    "            data[i] = data[i].to(device)\n",
    "            \n",
    "        output = model(data)\n",
    "        return torch.squeeze(torch.argmax(output, dim=1)).cpu().item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#with redirect_output(\"deepSNAP_redirect.txt\"):\n",
    "train_loss = []\n",
    "test_loss = []\n",
    "train_accuracy = []\n",
    "test_accuracy = []\n",
    "#drop_prob=0.5\n",
    "print(vocab_size)\n",
    "model = Net(vocab_size, output_size=512, embedding_dim=256, hidden_lstm_dim=1024, lstm_layers=1).to(device)\n",
    "# model.linear1_1.weight.requires_grad = False\n",
    "# model.linear1_1.bias.requires_grad = False\n",
    "# model.linear1_2.weight.requires_grad = False\n",
    "# model.linear1_2.bias.requires_grad = False\n",
    "\n",
    "if do_learn: # training mode\n",
    "\n",
    "    optimizer = optim.Adam(model.parameters(), lr=lr, weight_decay=weight_decay)\n",
    "    #optimizer = optim.SGD(model.parameters(), lr=lr)\n",
    "\n",
    "    for epoch in range(num_epochs):\n",
    "        train(model, device, train_loader, epoch, optimizer)\n",
    "        test(model, device, test_loader)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.save(model, 'models/lstm_97.7%_v3.0.pt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "remove modifications  \n",
    "use one charge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = torch.tensor([[1., 1.], [1., 2.], [3., 4.]])\n",
    "S = torch.tensor([[2., 3.], [3., 2.], [4., 3.], [5, 3]])\n",
    "\n",
    "mul = torch.mm(T, S.t())\n",
    "print('mul: ')\n",
    "print(mul)\n",
    "adder = torch.tensor([1., 2., 3., 4.])\n",
    "added = adder + mul\n",
    "print('added: ')\n",
    "print(added)\n",
    "\n",
    "norm = T.pow(2).sum(1)\n",
    "print(norm)\n",
    "exp_norm = norm.expand(4, -1).t()\n",
    "print(exp_norm)\n",
    "# pdist = nn.PairwiseDistance(p=2)\n",
    "# output = pdist(T, S)\n",
    "# output\n",
    "# print(T)\n",
    "# print(T.t())\n",
    "# test = torch.tensor([1, 2, 3])\n",
    "# expanded_test = test.expand(3, -1)\n",
    "# print(expanded_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pairwise_distance(A, B):\n",
    "    \"\"\"Compute the 2D matrix of distances between all the embeddings.\n",
    "\n",
    "    Args:\n",
    "        embeddings: tensor of shape (batch_size, embed_dim)\n",
    "        squared: Boolean. If true, output is the pairwise squared euclidean distance matrix.\n",
    "                 If false, output is the pairwise euclidean distance matrix.\n",
    "\n",
    "    Returns:\n",
    "        pairwise_distances: tensor of shape (batch_size, batch_size)\n",
    "    \"\"\"\n",
    "    # Get the dot product between all embeddings\n",
    "    # shape (batch_size, batch_size)\n",
    "    dot_product = torch.mm(A, B.t())\n",
    "\n",
    "    # Get squared L2 norm for each embedding. We can just take the diagonal of `dot_product`.\n",
    "    # This also provides more numerical stability (the diagonal of the result will be exactly 0).\n",
    "    # shape (batch_size,)\n",
    "    \n",
    "    A_L2_norm = A.pow(2).sum(1)\n",
    "    B_L2_norm = B.pow(2).sum(1)\n",
    "    \n",
    "    # Compute the pairwise distance matrix as we have:\n",
    "    # ||a - b||^2 = ||a||^2  - 2 <a, b> + ||b||^2\n",
    "    # shape (batch_size, batch_size)\n",
    "    distances = A_L2_norm[:, None] - (2.0 * dot_product) + B_L2_norm\n",
    "    \n",
    "    # Because of computation errors, some distances might be negative so we put everything >= 0.0\n",
    "    distances = torch.max(distances, torch.tensor(0.0))\n",
    "\n",
    "#     if not squared:\n",
    "#         # Because the gradient of sqrt is infinite when distances == 0.0 (ex: on the diagonal)\n",
    "#         # we need to add a small epsilon where distances == 0.0\n",
    "#         mask = tf.to_float(tf.equal(distances, 0.0))\n",
    "#         distances = distances + mask * 1e-16\n",
    "\n",
    "#         distances = tf.sqrt(distances)\n",
    "\n",
    "#         # Correct the epsilon added: set the distances on the mask to be exactly 0.0\n",
    "#         distances = distances * (1.0 - mask)\n",
    "\n",
    "    return distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dists = pairwise_distance(T, T)\n",
    "print(T)\n",
    "print(S)\n",
    "print(dists)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = torch.tensor([[1., 1.], [1., 2.], [3., 4.]])  \n",
    "S = torch.tensor([[2., 3.], [3., 2.], [4., 3.], [5, 3]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pairwise_distances(x, y=None):\n",
    "    '''\n",
    "    Input: x is a Nxd matrix\n",
    "           y is an optional Mxd matirx\n",
    "    Output: dist is a NxM matrix where dist[i,j] is the square norm between x[i,:] and y[j,:]\n",
    "            if y is not given then use 'y=x'.\n",
    "    i.e. dist[i,j] = ||x[i,:]-y[j,:]||^2\n",
    "    '''\n",
    "    x_norm = (x**2).sum(1).view(-1, 1)\n",
    "    if y is not None:\n",
    "        y_t = torch.transpose(y, 0, 1)\n",
    "        y_norm = (y**2).sum(1).view(1, -1)\n",
    "    else:\n",
    "        y_t = torch.transpose(x, 0, 1)\n",
    "        y_norm = x_norm.view(1, -1)\n",
    "    \n",
    "    dist = x_norm + y_norm - 2.0 * torch.mm(x, y_t)\n",
    "    # Ensure diagonal is zero if x=y\n",
    "    if y is None:\n",
    "        dist = dist - torch.diag(dist.diag())\n",
    "    dist[dist != dist] = 0 # set all nan values to zero\n",
    "    return torch.clamp(dist, 0.0, np.inf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = torch.tensor([[ 1.4335, -1.0990, -0.8586],\n",
    "        [ 2.1553,  2.7028, -0.8020],\n",
    "        [ 1.0524,  0.1599, -0.0374]])\n",
    "b = torch.tensor([[ 3., 7., 2.],\n",
    "                [ 6.,  8., 2.],\n",
    "                [ 9.,  5., 7.]])\n",
    "dists = process.pairwise_distances(x=a, y=None)\n",
    "print(dists)\n",
    "print(b.max(1).values)\n",
    "print(b[b.max(1).indices])\n",
    "#print(b.max(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = torch.tensor([[ 3., 1., 4.],\n",
    "                [ 0.,  2., 3.],\n",
    "                [ 5.,  3., 6.]])\n",
    "b = torch.tensor([[ 1., 2., 1.],\n",
    "                [ 2.,  2., 1.],\n",
    "                [ 1.,  1., 1.]])\n",
    "dist = (a - b) ** 2\n",
    "print(l2_squared(a, b))\n",
    "print(torch.sum(l2_squared(a, b), 1))\n",
    "print(a.min(1).values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a0 = torch.tensor([1., -1., 3.])\n",
    "max0 = torch.max(a0, torch.tensor(0.))\n",
    "print(max0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.asarray([1, 2, 3, 4])\n",
    "b = np.asarray([0, 1, 0])\n",
    "a.append(b)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_print(i):\n",
    "    print(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pep_path = join(in_tensor_dir, 'peptides')\n",
    "charge1 = 0\n",
    "charge2 = 0\n",
    "charge3 = 0\n",
    "for file in listdir(pep_path):\n",
    "    file_parts = re.search(r\"(\\d+)-(\\d+)-(\\d+.\\d+)-(\\d)-(0|1).pep\", file)\n",
    "    charge = int(file_parts[4])\n",
    "    if charge == 1:\n",
    "        charge1 += 1\n",
    "    elif charge == 2:\n",
    "        charge2 += 1\n",
    "    elif charge == 3:\n",
    "        charge3 += 1\n",
    "print('charge 1 count: {}'.format(charge1))\n",
    "print('charge 2 count: {}'.format(charge2))\n",
    "print('charge 3 count: {}'.format(charge3))\n",
    "\n",
    "unmod = 0\n",
    "mod = 0\n",
    "for file in listdir(pep_path):\n",
    "    file_parts = re.search(r\"(\\d+)-(\\d+)-(\\d+.\\d+)-(\\d)-(0|1).pep\", file)\n",
    "    m = int(file_parts[5])\n",
    "    if m == 0:\n",
    "        unmod += 1\n",
    "    elif m == 1:\n",
    "        mod += 1\n",
    "print('Unmodified count: {}'.format(unmod))\n",
    "print('Modified count: {}'.format(mod))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scratch Work Below This:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = torch.FloatTensor([1., 2., 3., 4., 5.]).unsqueeze(dim=0)\n",
    "s2 = torch.FloatTensor([6., 7., 8., 9., 10.]).unsqueeze(dim=0)\n",
    "print(s1)\n",
    "torch.cat((s1, s2), dim=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(indices=tensor([[0, 0, 0, 0, 0],\n",
      "                       [1, 3, 5, 7, 9]]),\n",
      "       values=tensor([1., 1., 1., 1., 1.]),\n",
      "       size=(1, 10), nnz=5, layout=torch.sparse_coo)\n",
      "tensor(indices=tensor([[0, 0, 0, 0, 0],\n",
      "                       [1, 2, 5, 8, 9]]),\n",
      "       values=tensor([2., 2., 2., 2., 2.]),\n",
      "       size=(1, 10), nnz=5, layout=torch.sparse_coo)\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "sparse tensors do not have strides",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-14-4899f5aaab63>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0mspec_lst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mspec1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mspec2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0mspecs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspec_lst\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspecs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspec1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_dense\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspec2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_dense\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: sparse tensors do not have strides"
     ]
    }
   ],
   "source": [
    "ind = torch.LongTensor([[0]*5, [1, 3, 5, 7, 9]])\n",
    "val = torch.FloatTensor([1., 1., 1., 1., 1.])\n",
    "spec1 = torch.sparse_coo_tensor(ind, val, torch.Size([1, 10]))\n",
    "print(spec1)\n",
    "ind = torch.LongTensor([[0]*5, [1, 2, 5, 8, 9]])\n",
    "val = torch.FloatTensor([2., 2., 2., 2., 2.])\n",
    "spec2 = torch.sparse_coo_tensor(ind, val, torch.Size([1, 10]))\n",
    "print(spec2)\n",
    "spec_lst = [spec1, spec2]\n",
    "specs = torch.cat(spec_lst)\n",
    "print(specs[1, :])\n",
    "print(spec1.to_dense())\n",
    "print(spec2.to_dense())\n",
    "print(specs.to_dense())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "length = 5\n",
    "width = 15\n",
    "count = np.array([5, 3, 4, 1, 2])\n",
    "PQ_mask = torch.zeros(length, width)\n",
    "\n",
    "rows = []\n",
    "[rows.extend([i]*x) for i,x in enumerate(count)]\n",
    "cols = range(width)\n",
    "PQ_mask[rows, cols] = 1.\n",
    "print(c_sums)\n",
    "print(PQ_mask)\n",
    "QQ_mask = torch.zeros(width, width)\n",
    "QQ_mask[cols, :] = PQ_mask[rows, :]\n",
    "print(QQ_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[1., 4., 9., 7.],\n",
      "        [2., 3., 2., 6.],\n",
      "        [7., 1., 3., 5.]])\n",
      "tensor([3.3333, 2.6667, 4.6667, 6.0000])\n",
      "tensor([[-2.3333,  1.3333,  4.3333,  1.0000],\n",
      "        [-1.3333,  0.3333, -2.6667,  0.0000],\n",
      "        [ 3.6667, -1.6667, -1.6667, -1.0000]])\n"
     ]
    }
   ],
   "source": [
    "specs = torch.FloatTensor([[1., 4., 9., 7.], [2., 3., 2., 6.], [7., 1., 3., 5.]])\n",
    "print(specs)\n",
    "means = specs.mean(dim=0)\n",
    "print(means)\n",
    "print(specs - means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9b3/8deHsCfDmpBAIGySsMlmXOrSqtjWfUOrtnpdi7XVVttbW9tabe+jv262t3uVqq293ewFRErdwLreKgqIAsq+BsiCbFkg23x+f5whQAQZQmZOMvN+Ph48yJw5yXnnkPnk8J3z+X7N3RERkfTRIewAIiKSXCr8IiJpRoVfRCTNqPCLiKQZFX4RkTTTMewA8cjOzvYhQ4aEHUNEpF1ZuHDhNnfPab69XRT+IUOGsGDBgrBjiIi0K2a24VDbNdQjIpJmVPhFRNKMCr+ISJpR4RcRSTMq/CIiaSZhhd/MHjWzcjNb2mz7HWa2wsyWmdmPEnV8ERE5tERe8f8BOPfADWZ2FnAJMM7dxwAPJPD4IiJyCAkr/O7+MrC92ebbgB+4e21sn/JEHV9EpF2r3wNP3Q3V21r9Syd7jL8QOMPM5pvZS2Z24uF2NLOpZrbAzBZUVFQkMaKISBvw3LfgjYeg9J1W/9LJLvwdgd7AKcBXgb+bmR1qR3ef5u7F7l6ck/OBjmMRkdS1/J/w5sPwkdth+Nmt/uWTXfhLgJkeeAOIAtlJziAi0nbt3gJPfgHyxsHkbyfkEMku/LOAswHMrBDoDLT+AJaISHsUbYSZU6GhFq54FDp2SchhEjZJm5n9FTgTyDazEuA+4FHg0dgtnnXA9a5Ff0VEAv/3c1j/Clz8K8gekbDDJKzwu/s1h3nq2kQdU0Sk3SpZCC98D8ZcBhMTWybVuSsiEra9u2HGTRDpDxf+DA59z0uraRfz8YuIpLSnvgo7N8KNT0O3Xgk/nK74RUTC9M7f4Z2/wce+BgWnJOWQKvwiImHZvg7mfBkKPgJn/GfSDqvCLyIShsZ6mHEzWAe4fBpkJG/kXWP8IiJhePH7sHkhXPkH6FWQ1EPril9EJNnWvQyv/BQmXhfcvplkKvwiIslUsx1m3gp9h8N5PwwlgoZ6RESSxR2evB2qK+CaedA5M5QYKvwiIsmy4BFY8U/4xPdgwITQYmioR0QkGcrfg2e/CcMnwymfDzWKCr+ISKLV74XpN0OXCFz2IHQIt/RqqEdEJNHm3gvly+Az0yGrX9hpdMUvIpJQK56GN6YFwzsjPh52GkCFX0QkcXZvhVmfh7zj4Zz7w07TRIVfRCQRolGY9Tmo3wNTEreaVktojF9EJBH+/QtY+yJc9HPIKQw7zUESdsVvZo+aWXlsmcXmz/2nmbmZaaF1EUk9mxfCv/4LRl0Mk64PO80HJHKo5w/Auc03mtkg4OPAxgQeW0QkHLWVwa2bWXlw8S8SvppWSySs8Lv7y8D2Qzz138DdgBZZF5HU89TdsHMDTPkddOsddppDSuqbu2Z2MbDZ3d9O5nFFRJJiyXR4+y/w0a/C4FPDTnNYSXtz18y6A98EPhHn/lOBqQAFBcmdq1pE5KjtWA9z7oJBJ8NH7w47zYdK5hX/cGAo8LaZrQcGAovMLO9QO7v7NHcvdvfinJycJMYUETlKjfUw4xbA4PLfJXU1rZZIWjp3XwI09SrHin+xu29LVgYRkYR46YdQ8iZc8Sj0Hhx2miNK5O2cfwVeA4rMrMTMbk7UsUREQrP+VXj5AZhwLYydEnaauCTsit/drznC80MSdWwRkaSo2Q4zp0KfYaGtptUSbXsgSkSkrXKH2XdAVTncMhe6ZIWdKG6aq0dEpCUW/gGWz4HJ34YBE8NOc1RU+EVEjlb5cnjmHhh2Fnzk9rDTHDUVfhGRo1G/F2bcDJ27t4nVtFpCY/wiIkdj3n1QthQ+/XeIHLINqc1rf7+qRETCsvJZmP8gnPw5KPxk2GlaTIVfRCQelaUw6zbIHQvnfCfsNMdEhV9E5EiiUXjic1BXA1MegU5dw050TDTGLyJyJK/9Cta+ABf+N/QbGXaaY6YrfhGRD7PlLXj+uzDqIjjhxrDTtAoVfhGRw6mtClbTysyBi9rmalotoaEeEZHDefprsH0tXP8P6N4n7DStRlf8IiKHsnQGLP4TnPEVGHpG2GlalQq/iEhzOzbAP+6EgSfCmV8PO02rU+EXETlQY0OwmpY7THkYMjqFnajVaYxfRORAL/8ISt6Ayx+G3kPCTpMQuuIXEdlnw7/h5R/D+Gtg3JVhp0kYFX4REYA9O2DGZ4Or/PN/HHaahErkmruPmlm5mS09YNuPzWy5mb1jZk+YWa9EHV9EJG7uMPuLUFUajOt3iYSdKKESecX/B+DcZtvmAmPdfRywErgngccXEYnPoj/Ce7Ph7Hsh/4Sw0yRcwgq/u78MbG+27Tl3b4g9fB0YmKjji4jEpWIlPPN1GPoxOPWLYadJijDH+G8Cnj7ck2Y21cwWmNmCioqKJMYSkbTRUAszboKOXeGyh9rlalotEcp3aWbfBBqAPx9uH3ef5u7F7l6ck5OTvHAikj7m3Q+lS+DS30CP/mGnSZqk38dvZtcDFwKT3d2TfXwREQBWzYXXfwMnTYWi88JOk1RJLfxmdi7wNeBj7l6TzGOLiDSpKg9W0+o3Gj7+3bDTJF0ib+f8K/AaUGRmJWZ2M/ArIALMNbPFZvZgoo4vInJI+1bTqq2EKx6FTt3CTpR0Cbvid/drDrH5kUQdT0QkLq//BtY8Dxf8BPqNCjtNKNLjLWwREYAti4M3dIsugOKbw04TGhV+EUkPddUw42bIzIaLf5kyq2m1hGbnFJH08PTX4P01cP1syOwbdppQ6YpfRFLfsifgrf+B0++CoR8NO03oVPhFJLXt3AizvxTMwXPWN8JO0yao8ItI6mpsCKZa9mjKrqbVEhrjF5HU9coDsOl1uGwa9BkWdpo2Q1f8IpKaNrwGL/0Qxl0F468KO02bosIvIqlnzw6Y+VnoVQDnPxB2mjZHQz0iklrc4R93QuVWuOlZ6Noj7ERtjq74RSS1vPUneHdWcAfPwOKw07RJKvwikjq2rYKn74YhZ8Bpd4adps1S4ReR1NBQC9Nvgo5d4PJp0CEj7ERtlsb4RSQ1PP9dKH0Hrv4L9BgQdpo2TVf8ItL+rZ4Hr/0qmHFz5AVhp2nzVPhFpH2rqoAnboOcUfDJ74Wdpl3QUI+ItF/RaLCE4t5d8B+z0nI1rZZI5NKLj5pZuZktPWBbHzOba2arYn/3TtTxRSQNzH8QVs8NrvRzx4Sdpt1I5FDPH4Bzm237OvC8u48Ano89FhE5elvfhnn3QeF5cOItYadpVxJW+N39ZWB7s82XAI/FPn4MuDRRxxeRFFZXDTNugW594JJfp/VqWi0RV+E3sxlmdoGZHesvilx33woQ+7vfhxxzqpktMLMFFRUVx3hYEUkpz9wTNGtd/lDar6bVEvEW8t8CnwZWmdkPzGxkAjMB4O7T3L3Y3YtzcnISfTgRaS+WzYJFj8FpX4JhZ4adpl2Kq/C7+zx3/wwwCVgPzDWzf5vZjWZ2NCsblJlZf4DY3+VHG1hE0tjOTfCPL8KAiXDWN8NO027FPXRjZn2BG4BbgLeAnxP8Iph7FMebDVwf+/h64Mmj+FwRSWfRRpg5Nfh7yiPQsXPYidqtuO7jN7OZwEjgf4CL9o3TA4+b2YLDfM5fgTOBbDMrAe4DfgD83cxuBjYCVx5bfBFJG6/8BDb+Gy59EPoODztNuxZvA9ev3P1fh3rC3Q8576m7X3OYrzU5zmOKiAQ2zocXfwDHXwnjrw47TbsX71DPKDPrte+BmfU2s88nKJOIyH57dga3bvbMhwt+ols3W0G8hf+z7r5z3wN33wF8NjGRRERi3GHOXbB7M0x5FLr2DDtRSoi38Hcw2/9r1swyAL2zIiKJtfgvsGwmnHUPDDox7DQpI94x/mcJ3pR9EHDgc8AzCUslIrJtNTz1VRh8Opz+5bDTpJR4C//XgFuB2wADngMeTlQoEUlzDXUw42bI6BR052o1rVYVV+F39yhB9+5vExtHRAT413dh62K46k/Qc2DYaVJOvPfxjwC+D4wGuu7b7u7DEpRLRNLVmn/Bv38JJ9wIoy4KO01KivfN3d8TXO03AGcBfyRo5hIRaT3V2+CJz0F2EXzy/4WdJmXFW/i7ufvzgLn7Bne/Hzg7cbFEJO24w6zPw54dcMUj0Ll72IlSVrxv7u6NTcm8ysxuBzbzIVMqi4gctfkPwapn4dwfQt7xYadJafFe8d8JdAe+CJwAXMv+ydZERI5N6VKYey+M+CScfGvYaVLeEa/4Y81an3L3rwJVwI0JTyUi6aOuBqbfBN16w6W/0ZQMSXDEwu/ujWZ2gpmZu3syQolIGnn2G7BtBVw7EzKzw06TFuId438LeNLM/heo3rfR3WcmJJWIpId3Z8PC38Opd8Bxmrg3WeIt/H2A9zn4Th4HVPhFpGV2lcDsO6D/BDj722GnSSvxdu5qXF9EWk+0EWbeCo31Wk0rBPF27v6e4Ar/IO5+U6snEpHU9+pPYcOrcMmvIfu4sNOknXiHeuYc8HFX4DJgS+vHEZGUt+kNeOH7MOZymPCZsNOkpXiHemYc+Di2nu68lh7UzO4iWLTdgSXAje6+t6VfT0Taib27glk3e+TDhf+tWzdDEm8DV3MjgIKWfKKZ5RM0ghW7+1ggA9AimiKpzh3mfBl2bYYpD0O3Xkf+HEmIeMf4Kzl4jL+UYI7+YzluNzOrJ+gI1rCRSCqrrQru1186Hc76JhScHHaitBbvUE+ktQ7o7pvN7AFgI7AHeM7dn2u+n5lNBaYCFBS06D8XItIWbJwPT9wKO9bDqV+EM74SdqK0F9dQj5ldZmY9D3jcy8wubckBzaw3cAkwFBgAZJrZtc33c/dp7l7s7sU5OTktOZSIhKmhDp7/Lvz+XPBGuOGf8In/0mpabUC8Y/z3ufuufQ/cfSdwXwuPeQ6wzt0r3L2eoAns1BZ+LRFpi8rfg4cnwys/gfGfhs/9Hww5LexUEhPv7ZyH+gUR7+c2txE4xcy6Ewz1TAYWtPBriUhbEo3C/N/CvO9Alwhc/RcYeUHYqaSZeIv3AjP7KfBrgjd57wAWtuSA7j7fzKYDiwhW9HoLmNaSryUibcjOTTDrNlj/ChSeBxf/ArK0bEdbFG/hvwO4F3g89vg54FstPai730fLh4pEpC1xh3ceh6e+Ch6Fi38JE6/TPfptWLx39VQDX09wFhFpb2q2w5w74d0nYdApcNmD0Gdo2KnkCOK9q2eumfU64HFvM3s2cbFEpM1bNRd+cwosfwrOuR9ufEpFv52Id6gnO3YnDwDuvsPMNHgnko7qquG5b8GCR6HfaPjMdOg/LuxUchTiLfxRMytw940AZjaEQ8zWKSIpbtOb8MRU2L4uWDzlrG9Bp65hp5KjFG/h/ybwqpm9FHv8UWJdtSKSBhrq4KUfBtMp9xgIN8yBIaeHnUpaKN43d58xs2KCYr8YeJLgHnwRSXXly4Or/K1vB9Mon/sD6Noj7FRyDOKdpO0W4EvAQILCfwrwGgcvxSgiqSQahfkPwrz7oUsWXPUnGHVR2KmkFcQ71PMl4ETgdXc/y8xGAt9JXCwRCdWukqAZa93LUHguXPQLiOSGnUpaSbyFf6+77zUzzKyLuy83s6KEJhOR5HOHd/4eNGNFG+Cin8Ok69WMlWLiLfwlsfv4ZwFzzWwHmkNfJLXUbIc5d8G7s2DQybFmrGFhp5IEiPfN3ctiH95vZi8APYFnEpZKRJJr1Tx48gtQ8z5M/jacdqemT05hRz3Dpru/dOS9RKRdqKuG5+6FBY9Azkj4zN+h//iwU0mCtXRqZRFp70oWwMypsH0tfOR2OPteNWOlCRV+kXTTWA8v/ShYJCXSH66fDUM/GnYqSSIVfpF0UrEiuMrfuhjGXwPn/RC69jzy50lKUeEXSQfRKLwxDebdB526w6f+CKMvCTuVhESFXyTV7doMT34e1r4IIz4RLJQSyQs7lYQolMIf6wl4GBhLMMvnTe7+WhhZRFLakunwzy9DYwNc+DM44QY1Y0loV/w/B55x9yvMrDPQPaQcIqmpZjv88yuwbCYMPCloxuo7POxU0kYkvfCbWQ+CaZ1vAHD3OqAu2TlEUtbqefDk7VBdEdyiedqdkKFRXdkvjJ+GYUAF8HszGw8sBL4UW9e3iZlNJTbnf0FBQdJDirQ7dTUw99vw5u+CZqxr/gYDJoSdStqguNbcbWUdgUnAb919InDIhdzdfZq7F7t7cU5OTrIzirQvJQvhoTOCon/KF2Dqiyr6clhhXPGXACXuPj/2eDqHKPwiEofGenj5AXj5x0Ez1n/MhmEfCzuVtHFJL/zuXmpmm8ysyN1XAJOBd5OdQ6Td27YKZn4WtrwF466C834E3XqFnUragbDe8bkD+HPsjp61wI0h5RBpf6JRePPhYDy/U1e48jEYc2nYqaQdCaXwu/tioDiMY4u0a7u3wKzPw9oX4LiPwyW/UjOWHDXd4yXSXiyZHtyb31gHF/wUim9SM5a0iAq/SFu3Z0dQ8JfOgPxiuHyamrHkmKjwi7Rla/4Fs74A1eVw1rfg9LvUjCXHTD9BIm1RXQ3Mux/eeAiyC+Gav8CAiWGnkhShwi/S1mxeCDNvhfdXwcm3wTn3QaduYaeSFKLCL9JWNNYHq2K99KPgTp3rZsHws8JOJSlIhV+kLdi2Gp6YGlztH/8pOP/HasZKUzV1Dawqq2JFaSUryiq5/iNDKOjbuhMYq/CLhMk9aMZ67l7o2AWu+D2MvTzsVJIE9Y1R1m2rZnlpJStjRX5FaSWbdtTgHuzTpWMHzhiRrcIvkjJ2b4UnvwBrnofhk+GSX0OP/mGnklYWjTqbd+4JCnysuK8orWTttirqG4MKn9HBGJqdyfH5PbnihIEU5kYoyotQ0Kc7GR1av1dDhV8kDEtnwJwvQ0MtnP8AnHiLmrHaOXdnW1Vd0xDNytJKlpdVsqqskpq6xqb98nt1oygvwtmj+lGUG6EwN8Lwfpl06ZiRtKwq/CLJtGcHPPVVWPK/kH8CXDYNso8LO5Ucpcq99bGr9ypWllWyvHQ3K8uq2F69f02pPpmdKcqN8KniQRTlBQW+MDeLSNdOISYPqPCLJMuaF4J5dqrK4MxvwBlfUTNWG1fb0Mjq8qqmIr8iVuA379zTtE/3zhkU5kb4xOjcpiGaorwI2VldQkz+4fRTJ5Jo9XuCZqz5D0LfEXDLPMifFHYqOUBj1Nm4vYYVpbuDAl+2mxWllax/v4bGaDAO3ynDGJ6TRfGQ3nw6t4CiWJHP79WNDgkYh08kFX6RRNq8CJ64FbathJNuhXPuh86te4eGxM/dKdtdGxuaqWwq8qvKqqhtiALBWy0FfbpTlBvh/OP7U5gbYWRehCHZmXTKCGPRwtanwi+SCI0N8OpP4aUfQmY/uO4JGH522KnSys6a4I3WlWX7b5VcUVrJ7r0NTfv0i3ShKC/CdacMpjAvKPDH9cuie+fULo2p/d2JhGHb6uAqf/MCGHsFXPAAdOsddqqUtaeukVXllU1Fft9tk2W7a5v2iXTtyMi8CBeNHxCMwcfupumd2TnE5OFR4RdpLe6w4JGgGSujE0x5BI6/IuxUKaO+Mcr6bdX7b5WMFfgN2w9ueBqRm8Vpx2U3jcEX5UXI69EV0+2yTUIr/GaWASwANrv7hWHlEGkVu7fC7Nth9bxgSOeSX0OPAWGnapfcnZIdez4wRLO2opq6xmAcvoPB0OxMRg/owaUT8xkZu11ycN/MhDQ8pZowr/i/BLwH9Agxg8ixW/YEzLkL6veqGesovV9Vu7/hKTZMs6qsiqra/ePw+b26UZibxceKcpoK/PCcLLp2Sl7DU6oJpfCb2UDgAuB7wJfDyCByzHZsgBe+B+88DgMmBStjZY8IO1WbVFXbwMqyg+ekWVlWybaq/Q1Pvbt3oigvwpRJ+RTl9aAoL4sRuRF6tIGGp1QT1hX/z4C7gcjhdjCzqcBUgIKCgiTFEjmCyrLgCn/pDCh5AywDzrwn1oylAlXXEGVNRdVBc9KsKKukZMfBDU8jciOcPbJfUOBzIxTmZZGT1UXj8EmS9MJvZhcC5e6+0MzOPNx+7j4NmAZQXFzsSYon8kF7dsC7s4Niv/4V8CjkHg+T74OxU6D34LATJl10X8PTAXPSrCytZN22ahpiDU8dOwQNTxMLenPNSQVBV2tuhIG921/DU6oJ44r/NOBiMzsf6Ar0MLM/ufu1IWQRObTaKljxdFDsV8+DaD30GQZn/GdQ7PuNDDthUrg75ZW1H7hVcmVZJXvro037FfTpTmFuhE+OyaMwdrvk0OxMOndMjYanVJP0wu/u9wD3AMSu+P9TRV/ahIZaWDU3KPYrn4H6GuiRDyffGtyW2X9CSr9pu6umnpXllQfND7+yrJKdNfVN++REulCUG+EzJw+ODdFEGNEvi8wuujO8PdG/lqS3xgZY9xIsnQnv/QNqd0H3vjD+mqDYDzoFOqTWVeve+mDisebzw5fu3tu0T6RLRwrzgikL9jU7FeVF6JOmDU+pJtTC7+4vAi+GmUHSUDQKm+YHV/bvzoLqCujSA0ZeCMdPgaFnpsSsmQ2NUda/X3PQ/PAryirZ8H41sWF4OnfswHE5WZw6vG/TEE1RXoT+PdXwlMra/0+3SDzcYevbQbFfOhN2l0DHrlB4bnBlf9zHoVPXsFO2iLuzZdfeppkl913Fr66ooq5hf8PTkL6ZFOVGuDg2bUFhboQhfbvTMUUmHpP4qfBLaqtYGSv20+H91dChY7DM4Tn3QdF50OWwdxS3Sdur64KZJUsrWVG2f374Axue+vfsSmFuhDNGZDcN0RzXTw1Psp8Kv6SenRuDq/ql06F0CWAw5HQ49Q4YdTF07xN2wiOq3tfwdMDUwStKq9hWtX/isZ7dgoanyyflNxX4wtwIPbupn0A+nAq/pIaqclg2Kyj2m+YH2/KL4dwfwJjLIJIXbr7DqGuIsnZb1f7pg2Pj8Ju272946tqpA4W5Ec4qymkq7iPzIuRE1PAkLaPCL+3Xnh3w3pyg2K97OWis6jcGJn871lg1JOyETaJRZ9OOmg/cD7+2Yn/DU0YHY1h2JuMH9uJTJwxqmh9+UO/uaniSVqXCL+1LXfX+xqpVc4PGqt5DgykTxk6BfqNCjefuVFTWHjSrZDBkU8We+sam/Qb27sbIvAjnjMptmjp4aHYmXTpqHF4ST4Vf2r6GWlj9fHBlv+LpoLEqMiBorBo7BQZMDKWxateeelY1mzp4ZVklOw5oeMrO6kxRXoSrTxrUdKvkiNwIWWp4khDpp0/apsaGYF6cpdODxqq9u6BbHxh/dVDsC05NWmPVvoanA8fgV5ZWsmXX/oanzM4ZFOZFOHds3kFvtGZndUlKRpGjocIvbUc0CiVvBsV+2RNBY1XnCIy6MCj2w85M6AyYDY1RNmyvOWh1pxVllazfdkDDU0YHhuVkctLQPgc1POX36qY3WqXdUOGXcLkHt1wunR7cgrlrU6yx6pNBsR/xCejUrZUP6WzdtXf/vPCxq/hV5fsbnizW8FSYm8WFx/dveqN1cN9MOqnhSdo5FX4Jx7bVsWI/A7atjDVWnQ1nfwuKzoeurbMw247quv1j8AdMW1C5d3/DU16PrhTmRTh1eN+m+eGP65dFt856o1VSkwq/JM/OTbBsJiyZDqXv0NRYdcptMOoSyOzb4i9dU9fAqrKqg1Z3Wl5aSUXl/oanHl07MjKvB5dMGBAboulBYW4Wvbpr4jFJLyr8klhV5fDuk0Gx3/R6sC3/BPjk92HMpUe9IHl9Y5R126oPWt1pZVklG7fX4LFx+C4dg4anj46IrdEaG4vP7aGGJxFQ4ZdE2LMTls8Jiv26l2KNVaPh7Hth7OXBgiZHEI06m3fuaSru+67i11RUUd+4v+FpaHYmYwf0ZMqkgU130xT06U6GGp5EDkuFX1pHXQ2sfBqWzIDVc6GxLuicPf0uGHsF5I4+5Ke5O9uq6vZ3s+57o7Wskuq6/Q1P+b26UZQX4ayR/Zrmhx/eTw1PIi2hwi8t11AHa54PruxXPA311ZCVByfeEhT7/EkHNVZV7q1nZdkH56XZXl3XtE+fzM4U5Ua4snhQ073whblZRLpq4jGR1qLCL0cn2hg0Vi2ZDu/NjjVW9YZxVwbFfvCp1EZhTXk1KxZvPmh++M0790881r1zBoW5ET4xOrdpiKYoTw1PIsmgwi9H5h40Vi3Z11hVDp2ziBadT/ngi3i78wSWl9ey8t+VLJ/5Cuvfr6Ex1vHUKcMYnpNF8ZDefDq34KCGJ008JhKOpBd+MxsE/BHIA6LANHf/ebJzyBG4Q9lSWDIdXzoD27WJaIfObOh7Oi/3PpMna8aw7K16at+MAkswg4I+3SnKDdZp3Td18JBsNTyJtDVhXPE3AF9x90VmFgEWmtlcd383hCzSTGXJcnYt+BuRVbPoWb2ORjrwGuOYWXcBz0WLqarpTr9IF4ryIlx3SqSpo/W4fll076z/QIq0B0l/pbr7VmBr7ONKM3sPyAdU+JNoT3Ulm1e+xfb1bxMtXUbmrpXk7l1HP7aT6cabXsSzHW5hfb9z6N9/IBPyIlwVu5umd6YankTas1Av0cxsCDARmH+I56YCUwEKCgqSmiuVNNTVsnntUratXUz9lqV02bGCnJq1DIiWcpwF4/B7vRMlHQvY0PNEVvUbR4fRFzFkWCH39uiqhieRFGS+r90x2Qc2ywJeAr7n7jM/bN/i4mJfsGBBcoK1Ux5tpHTTaspWv8Xeknfo+P4K+lStZmBjCZ0tmJem0Y2SjHze7z6M2j4j6Zo/luxhExgwdAwZHTVMI5JqzGyhuxc33x7Kq93MOgEzgD8fqejLB20vK2HLqkVUb3oHq3iPnrtXM7B+Pf1tL/1j+5SSQ1m3oSzsdTqd+o+l95Dx5I8Yz+BumQwONb2IhC2Mu3oMeAR4z91/muzjtydVu3dQsnIRu9e/jZe/S+auVQyoXUsfdtMnts8OImzpPJSlORdA7hh6Dh7HgMJJ5PXqS9tcXlxEwhbGFf9pwHXAEuyIaj4AAAiiSURBVDNbHNv2DXd/KoQsbULt3ho2r17CjnVvUb91GV13rCR37xr6ewUjY/vUeBdKOg1mda/TieaMInPQePoXTqRvv4H0TtJKVCKSGsK4q+dVIC3fMYw2NLBl/XIq1r7F3s1L6fz+crJr1pDfuJlhFiwAUu8ZlGQMZEvW8azvW0S3/HHkDJ9I/8GFFGZoXhoROXZ6Ry8BPBplW+lGSlctombTO3TYtpxeVasZWL+BgVbHwNh+my2X8m7D2dJ7Mp0GjKXv0AnkDz+eoV26MjTU70BEUpkK/zHatb2CLasWsXvD21D+HpHdq8ivW0sO1eTE9tlGL7Z2GcrbfS4jI28MPYeMJ3/EBPIjvcgPNb2IpCMV/jjtramiZNVidqxbTGPpu3TfuZK8vWvox3Z6xvap9G5s7jyU5X0mQ7/RZA0aR/8RE8jul092qOlFRPZT4W+mob6OzWvfZduat6jbupQu21eQU7OGAdGtTQ1Ptd6Jko6D2NjjBNZmj6LbwOPJPW4iuQOHM1JvtIpIG5e2hd+jUcpK1lC2eiE1JcvotO09elevZlDDJgZbA4MJGp62dOhPRffhlPQ5ny4DxtJ32ETyh41meCdNWyAi7VNaFP4dFVvZvHIBVRuXxBqeVjGwfj15tqfpXvcy+lLWdRiLck4lI28MvYeOZ+CICQzqnsWgUNOLiLSulC78rz16NyM2Pk42O+kd27aLTDZ3HsaynPOg32h6DB7HgBGTyO2TQ26oaUVEkiOlC39GzwGs6/URVmePJHPQOPJGTCI7r4CeGocXkTSW0oX/pCl3hh1BRKTN0aWviEiaUeEXEUkzKvwiImlGhV9EJM2o8IuIpBkVfhGRNKPCLyKSZlT4RUTSjLl72BmOyMwqgA0t/PRsYFsrxmktynV0lOvoKNfRaau54NiyDXb3nOYb20XhPxZmtsDdi8PO0ZxyHR3lOjrKdXTaai5ITDYN9YiIpBkVfhGRNJMOhX9a2AEOQ7mOjnIdHeU6Om01FyQgW8qP8YuIyMHS4YpfREQOoMIvIpJmUqbwm9m5ZrbCzFab2dcP8XwXM3s89vx8MxvSRnLdYGYVZrY49ueWJGR61MzKzWzpYZ43M/tFLPM7ZjYp0ZnizHWmme064Fx9O0m5BpnZC2b2npktM7MvHWKfpJ+zOHMl/ZyZWVcze8PM3o7l+s4h9kn66zHOXEl/PR5w7Awze8vM5hziudY9X+7e7v8AGcAaYBjQGXgbGN1sn88DD8Y+vhp4vI3kugH4VZLP10eBScDSwzx/PvA0YMApwPw2kutMYE4IP1/9gUmxjyPAykP8Oyb9nMWZK+nnLHYOsmIfdwLmA6c02yeM12M8uZL+ejzg2F8G/nKof6/WPl+pcsV/ErDa3de6ex3wN+CSZvtcAjwW+3g6MNnMrA3kSjp3fxnY/iG7XAL80QOvA73MrH8byBUKd9/q7otiH1cC7wH5zXZL+jmLM1fSxc5BVexhp9if5neRJP31GGeuUJjZQOAC4OHD7NKq5ytVCn8+sOmAxyV88AXQtI+7NwC7gL5tIBfAlNjwwHQzG5TgTPGIN3cYPhL7r/rTZjYm2QeP/Rd7IsHV4oFCPWcfkgtCOGexYYvFQDkw190Pe76S+HqMJxeE83r8GXA3ED3M8616vlKl8B/qN1/z3+Tx7NPa4jnmP4Ah7j4OmMf+3+phCuNcxWMRwdwj44FfArOSeXAzywJmAHe6++7mTx/iU5Jyzo6QK5Rz5u6N7j4BGAicZGZjm+0SyvmKI1fSX49mdiFQ7u4LP2y3Q2xr8flKlcJfAhz4m3kgsOVw+5hZR6AniR9WOGIud3/f3WtjD38HnJDgTPGI53wmnbvv3vdfdXd/CuhkZtnJOLaZdSIorn9295mH2CWUc3akXGGes9gxdwIvAuc2eyqM1+MRc4X0ejwNuNjM1hMMB59tZn9qtk+rnq9UKfxvAiPMbKiZdSZ482N2s31mA9fHPr4C+JfH3ikJM1ezceCLCcZpwzYb+I/YnSqnALvcfWvYocwsb9+4ppmdRPDz+34SjmvAI8B77v7Tw+yW9HMWT64wzpmZ5ZhZr9jH3YBzgOXNdkv66zGeXGG8Ht39Hncf6O5DCGrEv9z92ma7ter56tjST2xL3L3BzG4HniW4k+ZRd19mZt8FFrj7bIIXyP+Y2WqC35RXt5FcXzSzi4GGWK4bEp3LzP5KcLdHtpmVAPcRvNGFuz8IPEVwl8pqoAa4MdGZ4sx1BXCbmTUAe4Crk/DLG4IrsuuAJbHxYYBvAAUHZAvjnMWTK4xz1h94zMwyCH7R/N3d54T9eowzV9Jfj4eTyPOlKRtERNJMqgz1iIhInFT4RUTSjAq/iEiaUeEXEUkzKvwiImlGhV8kwSyYIfMDMy6KhEWFX0Qkzajwi8SY2bWx+doXm9lDsQm9qszsJ2a2yMyeN7Oc2L4TzOz12GReT5hZ79j248xsXmxStEVmNjz25bNik34tN7M/J2FmWJHDUuEXAcxsFHAVcFpsEq9G4DNAJrDI3ScBLxF0EwP8EfhabDKvJQds/zPw69ikaKcC+6ZtmAjcCYwmWJ/htIR/UyKHkRJTNoi0gskEE3K9GbsY70YwdW8UeDy2z5+AmWbWE+jl7i/Ftj8G/K+ZRYB8d38CwN33AsS+3hvuXhJ7vBgYArya+G9L5INU+EUCBjzm7vcctNHs3mb7fdgcJx82fFN7wMeN6LUnIdJQj0jgeeAKM+sHYGZ9zGwwwWvkitg+nwZedfddwA4zOyO2/Trgpdhc+CVmdmnsa3Qxs+5J/S5E4qCrDhHA3d81s28Bz5lZB6Ae+AJQDYwxs4UEqx5dFfuU64EHY4V9Lftn47wOeCg2s2I9cGUSvw2RuGh2TpEPYWZV7p4Vdg6R1qShHhGRNKMrfhGRNKMrfhGRNKPCLyKSZlT4RUTSjAq/iEiaUeEXEUkz/x/72LMam6zmWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU5dn/8c+VEJIQkrAEkrALggKuiLjgjuIOLriLaPXBLrb6tD/buvSx2j6tba2ttX1ssai4KyCKiKKigisKFBWJC+ICGiCAJED2zP3740xCVhggZ84s3/frxcuZkzM5F8fMlzv3nHPd5pxDRESSR0rQBYiISHQp+EVEkoyCX0QkySj4RUSSjIJfRCTJdAi6gEjk5eW5AQMGBF2GiEhcWbJkyQbnXI/m2+Mi+AcMGMDixYuDLkNEJK6Y2VetbddUj4hIklHwi4gkGQW/iEiSUfCLiCQZBb+ISJLxLfjNrK+ZvWpmRWb2kZldG97+azP7xsyWhf+c5lcNIiLSkp+Xc9YCP3POLTWzbGCJmb0U/tpfnHN3+HhsERFpg28jfudcsXNuafjxFqAI6O3X8UREEkr1Nnj+F1Cxud2/dVTm+M1sAHAwsCi86Roz+8DM7jOzrm28ZrKZLTazxSUlJdEoU0QkNlRtgYcnwLtTYPWine+/i3wPfjPrDMwErnPOlQH3AIOAg4Bi4M+tvc45N8U5N9I5N7JHjxZ3HIuIJKaKzfDQ2V7gnzsVhpzc7ofwNfjNLA0v9B9xzj0F4Jxb55yrc86FgHuBUX7WICISN8o3wUNnwbfL4PxpsN85vhzGz6t6DJgKFDnn7my0vbDRbmcDy/2qQUQkbmzbCA+Og3UfwQUPw9AzfTuUn1f1jAYmAh+a2bLwthuBi8zsIMABXwJX+1iDiEjs27oepo2D776Aix6Hvcf4ejjfgt859wZgrXxprl/HFBGJO2XF3ki/dA1c/CQMPNb3Q8ZFW2YRkYRUugamnemN+C+dCf2PjMphFfwiIkH47iuYdoZ3Fc/EWdA3ete5KPhFRKJt4+fenH71VrjsGeg9IqqHV/CLiERTyafenH5dNUx6FgoPiHoJCn4RkWhZtwIeHA84mDQH8ocFUobaMouIRMPaD705fUuBy+cGFvqg4BcR8d+3/4EHzoAOGXDFXOgxJNByNNUjIuKn1e/Bw+dCZq43p991QNAVacQvIuKbr972eu906uZN78RA6IOCX0TEH18shIfPgexCb3qnS9+gK2qg4BcRaW8r58Mj50GXfnD5c5DTK+iKmlDwi4i0p0/nwWMXQvfBXuhn5wddUQsKfhGR9lI0Bx6/BHoOg0mzISsv6IpapeAXEWkPH82C6ZOg10FeG4ZO3YKuqE0KfhGRPfXBkzDje9DnULj0KcjsEnRFO6TgFxHZE/95GJ6aDP1HwyUzICMn6Ip2SsEvIrK7Ft8Pz/wIBh7nLaKS3jnoiiKi4BcR2R2LpsCc62Dwyd5yiR07BV1RxBT8IiK76q274fnrYd8zvIXR0zKCrmiXqFePiMiuWHgHvPIbGHYWnPtvSE0LuqJdpuAXEYmEc/Da7bDgdtj/fDjrHkiNzwiNz6pFRKLJOZh/G7xxJxx0CYy7G1JSg65qtyn4RUR2xDl48WZ4++9wyBVw+p2QEt8fjyr4RUTaEgrBC7+Ad6fAqKvh1D+AWdBV7TEFv4hIa0Ih73LNpdPgiGtg7G8TIvRBwS8i0lKoDmb/GJY9Akf/DE74VcKEPij4RUSaqquFp78PH06H426EY3+eUKEPCn4Rke3qamDmVbDiaRhzCxz906Ar8oWCX0QEoLYKpl8BnzwHY/8Xjrwm6Ip8o+AXEamphCcnwmcvwql/gsMmB12RrxT8IpLcqsvh8Yth1Wtwxl9h5BVBV+Q7Bb+IJK+qrd76uF++AeP/AQdfEnRFUeHb7Wdm1tfMXjWzIjP7yMyuDW/vZmYvmdln4f929asGEZE2VZbBw+fCV2/BOfcmTeiDv22Za4GfOeeGAocDPzKzYcAvgfnOucHA/PBzEZHoqdgMD50N3yyGCVPhgPOCriiqfAt+51yxc25p+PEWoAjoDYwHpoV3mwac5VcNIiItlG+CB8dB8ftw/oMw/OygK4q6qHQaMrMBwMHAIiDfOVcM3j8OQM82XjPZzBab2eKSkpJolCkiiW7bBph2Jqz/GC58FPY9PeiKAuF78JtZZ2AmcJ1zrizS1znnpjjnRjrnRvbo0cO/AkUkOWxZBw+cDhs/h4sfhyFjg64oML4Gv5ml4YX+I865p8Kb15lZYfjrhcB6P2sQEaHsW3jgNNi8Gi6ZDoNOCLqiQPl5VY8BU4Ei59ydjb40G5gUfjwJeMavGkRE2Lwa7j/NG/FPfAr2OjroigLn53X8o4GJwIdmtiy87UbgduBJM7sS+BpIro/TRSR6Nn0B08ZBZSlc9jT0GRl0RTHBt+B3zr0BtNXSboxfxxURAby5/GlnQk05THoGeh0cdEUxQ3fuikjiKfnEG+mHamDSs1Cwf9AVxRQFv4gklnUrvOv0Mbj8Oeg5NOiKYk58rxgsItJY8fveJZspHeCKuQr9Nij4RSQxfLPEm9NP6+SN9PMGB11RzNJUj4jEv9Xveg3XMrt6c/pd+wddUUzTiF9E4tuXb3oN17J6eNM7Cv2dUvCLSPxa9Ro8MgFyenmhn9sn6IrigoJfROLTypfh0Qug6wBvTj+7IOiK4oaCX0TizyfPw2MXeR/gTpoDnVtt8ittUPCLSHxZMRueuBTyh8NlsyGre9AVxR0Fv4jEj+UzYfrl0GsEXPYMdOoWdEVxScEvIvHh/cdh5lXQ9zCvy2ZGbtAVxS0Fv4jEvqUPwazvw4Cj4NIZkJ4ddEVxTcEvIrHtvX/D7Gu8xVMufhI6ZgVdUdxT8ItI7HrnHnjuZzDkFG+N3LTMoCtKCAp+EYlNb94FL/wShp4J5z8EaRlBV5Qw1KtHRGLPgj/Bq7+F4efAOVMgNS3oihKKgl9EYodz8OrvYOEf4YALYfw/IFUx1d50RkUkNjgHL/8a3vwrHHwpnPk3SEkNuqqEpOAXkeA5B/NuhHf+D0ZeCafdASn6CNIvCn4RCVYoBM9f7122edgP4JTfg1nQVSU0Bb+IBCcUgjnXwtIH4cifwEm3KfSjQMEvIsEI1cEzP4L3H4Njrofjb1LoR4mCX0Sir64WZk32mq4dfxMc+/OgK0oqCn4Ria7aaph5JRTNhhNvhaOuC7qipKPgF5Hoqa2CJyfBp8/Dyb+HI34YdEVJScEvItFRXQ5PTvSWTDztDhj1X0FXlLQU/CLiv5JPYfokWF/k3Zh1yKSgK0pqCn4R8dcH0+HZa70ma5fOgL1PDLqipKfgFxF/1FTCvBtg8X3Q93CYcB/k9g66KkHBLyJ+2LTK+xB37Qcw+lo44VfqsBlDFPwi0r5WzPZuzDKDCx+DfU8LuiJpxrcuSGZ2n5mtN7Pljbb92sy+MbNl4T/6iRBJFLXV8MIN3pU73feGq19X6McoP0f8DwB/Bx5stv0vzrk7fDyuiETb5tUw/XL4ZjGMuhrG/gY6pAddlbTBt+B3zi00swF+fX8RiRGfzoNZV3ttGM57AIafHXRFshNBNLy+xsw+CE8FdQ3g+CLSHupqvYVTHj0fcvrA1QsU+nEi2sF/DzAIOAgoBv7c1o5mNtnMFpvZ4pKSkmjVJyKRKCuGB8fBG3+BQy6Hq16C7oOCrkoiFNWrepxz6+ofm9m9wJwd7DsFmAIwcuRI5391IhKRz1+FmVdBTTmcPQUOvCDoimQXRXXEb2aFjZ6eDSxva18RiTGhOnjtdnjobMjKg/96VaEfp3wb8ZvZY8BxQJ6ZrQFuAY4zs4MAB3wJXO3X8UWkHW0tgaeuglWvwQEXwhl3QsesoKuS3eTnVT0XtbJ5ql/HExGffPkmzPgeVG72GqyNuEwrZcU53bkrIq0LheCtu2D+b6DrAK/BWsH+QVcl7UDBLyItlW+CWd+Hz+bBsLNg3N2QkRN0VdJOFPwi0tSaxd5duFvWegumHHqVpnYSjIJfRDzOwaJ/wou/gpxCuHIe9D4k6KrEBwp+EYHKUq+jZtGzMORUOPseyNSN9Ykqouv4zexaM8sxz1QzW2pmY/0uTkSi4Ntl8K9j4OO5MPa3cNFjCv0EF+kNXN9zzpUBY4EewBXA7b5VJSL+cw7emwpTx3otla+YC0f+WPP5SSDSqZ76n4TTgPudc++b6adDJG5VbYU518GH02HQGDhninc3riSFSIN/iZm9COwF3GBm2UDIv7JExDfrVsD0SbBxJZxwMxz1M0gJolGvBCXS4L8Sr6PmKudcuZl1w5vuEZF4suxRmPNTSM+Gy56BvY4JuiIJQKTBfwSwzDm3zcwuBUYAd/lXloi0q+pyeP56+M/DMOBoOHcqZOcHXZUEJNLf7+4Bys3sQODnwFe0XFJRRGLRhs/g3yd6oX/M9TDxaYV+kot0xF/rnHNmNh64yzk31cwm+VmYiLSDD2fAs9dCake4ZCYMPjHoiiQGRBr8W8zsBmAicLSZpQJp/pUlInukphLm3QiLp0Lfw2DCfZDbJ+iqJEZEOtVzAVCFdz3/WqA38CffqhKR3bfpC7hvrBf6R/4YLn9OoS9NRDTid86tNbNHgEPN7AzgXeec5vhFYk3RHHj6h96dNxc+CvueHnRFEoMibdlwPvAucB5wPrDIzCb4WZiI7IK6Gph3EzxxCXQfCFcvVOhLmyKd478JONQ5tx7AzHoALwMz/CpMRCK0eTXMuALWvAejJnv9djqkB12VxLBIgz+lPvTDNhLlhdpFpBWfvgizJkNdLUy4H/Y7J+iKJA5EGvwvmNk84LHw8wuAuf6UJCI7VVcLr/4vvHEn5O8H502DvL2DrkriRKQf7l5vZucCo/E+NprinJvla2Ui0rqyYph5JXz1prfw+al/hLTMoKuSOBLxQizOuZnATB9rEZGdWfUazLwKqrfBWf+Egy4KuiKJQzsMfjPbArjWvgQ455xWXxaJhlAdLLwDXvs95A2BSc9Cz6FBVyVxaofB75zLjlYhItKGrSXw1H/Bqldh//PhjL9Aeuegq5I4pjV3RWLZV2/BjO9B+SY48y4YMUkrZMkeU/CLxKJQCN76G8y/Dbr2h6tehsIDgq5KEoSCXyTWlG+Cp38An74Aw8bDuLshIzfoqiSBKPhFYsmaJd6yiFvWepdpjpqsqR1pdwp+kVjgHCz6F7x4M2QXwvfmQZ9Dgq5KEpSCXyRolaXwzDVQNBuGnAJn3QOdugVdlSQwBb9IkIrfhycnweav4aTb4IgfQ4raYIm/FPwiQXAOljwAz//CG91f/hz0PyLoqiRJKPhFoq1qK8z5b/jwSRh0ApxzL2TlBV2VJBHffqc0s/vMbL2ZLW+0rZuZvWRmn4X/29Wv44vEpPVFcO/xsHwGHH8TXDJDoS9R5+dk4gPAKc22/RKY75wbDMwPPxdJDssegynHQ8VmmPg0HPtzSEkNuipJQr4Fv3NuIbCp2ebxwLTw42nAWX4dXyRm1FR4V+08/X3ofQh8/3UYeGzQVUkSi/Ycf75zrhjAOVdsZj3b2tHMJgOTAfr16xel8kTa2YaV3g1Z65bD0T+D426EVH20JsGK2Z9A59wUYArAyJEjW2sNLRLbls+E2T+B1DRvLn/wSUFXJAJEP/jXmVlheLRfCKzf6StE4k1tFcy7Cd67F/qMgvPuh9w+QVcl0iDad4rMBiaFH08Cnony8UX89d2XMHWsF/pHXANXzFXoS8zxbcRvZo8BxwF5ZrYGuAW4HXjSzK4EvgbO8+v4IlH38XMw6wfe4wsegaFnBFuPSBt8C37nXFuLgY7x65gigairgZd/DW//HQoPgvMegG57BV2VSJti9sNdkbhQugamXwFr3oVDr4KTfwcd0oOuSmSHFPwiuyNUBx884X2IW1cN506F/ScEXZVIRBT8IrvCOfjsJW9qZ/1H0GsEnDMF8gYHXZlIxBT8IpFaswRevgW+fB267gUT7ofhZ2uFLIk7Cn6Rndn4ubfo+YqnoVMenHYHjJgEHToGXZnIblHwi7Rl63pY8Aevb35qOhz7SzjyGkjPDroykT2i4BdprmoLvPV3eOtuqKuCQy6HY34O2flBVybSLhT8IvXqarzR/YI/wLYSGHYWjPkf6D4o6MpE2pWCX8Q5+GgWvPIb2LQK+h8FFz0BfQ4JujIRXyj4Jbl9sRBeugW+XQo9h8HF070umrpSRxKYgl+S09rl3rX4K1+CnD5w1j1wwAVaEUuSgoJfksvm1fDq/8L7j0NGLpz0Gxg1GdIygq5MJGoU/JIcyjfB63+Gd+/1no/+CRz135DZNdi6RAKg4JfEVlMBi/4Jr/8FqsrgoEvg+BvUI1+SmoJfElOoDpY9Cq/+DrZ8C0NOgTG3QP6woCsTCZyCXxKLc/DpC/DyrVBSBL1Hwrn3woCjgq5MJGYo+CVxrH4PXvof+Pot6DYIzn8Qho7TpZkizSj4Jf5t+Azm3wpFz0JWTzj9ThhxGaSmBV2ZSExS8Ev82rIWXrsdlj4IaZlw/E1w+A8hvXPQlYnENAW/xJ/KMnjrb/D2P7zVrw69Co65Hjr3CLoykbig4Jf4UVsNi++DhX+E8o0w/BwY8yvoNjDoykTiioJfYl8oBB895TVR++5L2OsYOPFW6D0i6MpE4pKCX2Lbqte8JmrFyyB/f7h0Jgwaoyt1RPaAgl9iU/EH3vq2n78Cuf3g7Cmw/3mQkhJ0ZSJxT8EvseW7r+CV38KHT3p9dE7+HYy8Uk3URNqRgl9iw7aN8Pod8N6/wVK8Bmqjr4PMLkFXJpJwFPwSrOpyWHQPvPFXqN4abqJ2I+T0CroykYSl4Jdg1NXCskfgtd/DlmLY5zSviVrPfYOuTCThKfglupyDT+Z6TdQ2fAJ9RsGE+6H/EUFXJpI0FPwSPV+/412aufod6D4YLngE9j1dl2aKRJmCX/xX8ok3wv/kOehcAGf8FQ6eCKn68RMJgt554p+yYm8O/z8PQVoWnHCz10StY1bQlYkktUCC38y+BLYAdUCtc25kEHWITypL4c274O3/g1AtjLraa6KW1T3oykSEYEf8xzvnNgR4fGlvtVXw3lRY+Ceo2OTdaXvCzdB1QNCViUgjmuqRPRcKwfIZXhO1zV/DwOPhpFuh8MCgKxOJeZU1dawvq6K4tIK1ZZWsLa2kuNT779qySm4dN5wD+7bvjYxBBb8DXjQzB/zLOTel+Q5mNhmYDNCvX78olycRWznf66mz9kMoOAAm3gWDTgi6KpGYsLWqlrWlFduDvLSS4rJK1tWHe1klm7ZVt3hd5/QOFORmUJibQci5dq8rqOAf7Zz71sx6Ai+Z2cfOuYWNdwj/YzAFYOTIke3/N5c98+1/4OVfe90zu/SDc/4N+52rJmqSFJxzfFdeEx6VNw32tWVeqK8rrWRLVW2L13bL6khBTgYFuRkc1K8LhTkZ5IdDvjA3g/ycDLIz/F02NJDgd859G/7vejObBYwCFu74VRITNn3hTeksnwmZ3eCU22Hk96BDetCVibSLupBjw9aqcJhXNIzS66dg1oWDvbo21OR1KQY9stMpyM1k7x6dOWrvvIZRe33Q5+dkkJGWGtDfbLuoB7+ZZQEpzrkt4cdjgduiXYfsom0bYMEfvRWwUjrA0f8PRv8EMnKDrkwkYlW19fPpleH59IomYb62tJL1W6qoCzWdZOiYmkJ+bjqFOZkc0KcLJw/3wrwwd/tovUfndDqkxsdvvEGM+POBWebdrdkBeNQ590IAdUgkqrd5a9u++TeoKYcRE+HYX0JOYdCViTSxraq2WYi3/LB0Yyvz6Z06pnqj8twMjhyU1/C4fpRemJtBt6yOWALdYR714HfOrQJ0uUesq6uBpQ/Cgj/A1nWw7xleE7UeQ4KuTJKMc47N5TVNQ7yV0fqWypbz6V07pZEfHpkf0KdLwyi9PtDzczPITu+QUKEeCV3OKU05B0WzYf5tsHEl9D0czn8I+h0WdGWSgOpCjo3h+fQdjdarms2nm0GPzukU5mawV14WRwzsTkFuZsOHo/XhHgvz6bFIwS+eqq3w2Tx45x5Y8x7k7QMXPgb7nKomarJbqmtDrCtrOn++PdwrGubTa5vNp6elWkN479c7l5OG5VOQm9lk6qVHdjppcTKfHosU/MmsYjN8+gKseMa7Hr+uCnJ6w7i74cCL1URN2lReXbv9uvSGqZf6xxWsLa1iw9aqFq/LTNs+n374oO6NrnjZPlrvntWRlBQNNvykd3ay2bYBPn7Om85ZtQBCNZDdC0ZeAUPHQb/DIUW/Hicr5xxlFbUUh69NX1fa9C5SL9wrKGtlPj03M60h1PfvndtoymX7aD0nI/nm02ORgj8ZlBXDx3O8kf1Xb4ILQZf+cPj3Yeh46H2IbrxKAqGQY8O2KtaVNm0P0HzUXlFT1+R1ZpAXnk/v370Thw3stv3D0ZwMCsPBntlRA4Z4oeBPVJu/hqJnYcVsWL0IcJA3BI76KQwb57VX0MgrYdTUhVi/pWqH7QHWb6mkpq7pfHqHFG8+vSA3g2G9chizb0/vUsbc7aP1nppPTzgK/kSy8XNvVF8022upAJC/Hxx3Awwbr/Vs41RFdd320Xkb7QE2bK2ieUuXjLSUhtH4YXt1awj0gvAoPT83nbysdM2nJyEFfzxzDko+9sJ+xWxY/5G3vdcIOPHX3px990FBVig74JyjrLK26SWMpVVNw72sks3lNS1em5PRwQv13AyGFeY0uYSxIDeDwpxMcjI1ny6tU/DHG+eg+H1vVL9iNmz8DDDvQ9mTfw9Dz4QufYOuMumFQo5N5dU7vOFobWkl5dV1LV5bP5/et1snDh3QrUW/l4LcDDp11FtXdp9+euJBKATfLIEVT3vz9pu/AkuBAUfBYVd7YZ9dEHSVSaO2fj69lcZd9eG+vqyK6rqmNx3Vz6fn56QztCCH4/fp2eKGo57ZGXTsoPl08ZeCP1aF6uDrt71RfdGzsOVbSEmDgcfBMf8P9jldSxn6oLKmrsmIvPEljPXPS7ZUEWplPr1+RF4/Sm98w1FBTgbdO6eTqvl0iQEK/lhSVwNfLPSmcT5+DraVQIcMGDQGht0CQ06BzPZdiSeZbKmsaf2Go9IK1pZ5V8R818p8enZGh4YrXPYpyG642Wj7B6UZ5GamaT5d4oaCP2g1lbDqVW9k/8lcqNwMaVkwZKz34ezgsZDeOegqY5pzjk3bqlvMnzcfrW9rdT69IwW5GfTuksEh/bt4V7s0/qA0J4OsdL1NJLHoJzoI1dtg5cte2H86D6q3QHou7HOKd9nloBMgLTPoKmNCbV2Ikq1VTW40atL7payCdWVVLRbFSE0xemanU5CbwT4F2Rw7pCcFuenbR+s5GfTMSSe9g246kuSj4I+WyjIv5Iuegc9ehtoKbwWr4Wd5Yb/XsdChY9BVRtXOFpleG77pqPl8escOKQ3hPaJf1/Dli+HWAOE59TzNp4u0ScHvp/JN8Mnz3pz9569AXTV0zoeDL/GmcfqPTthGaLu7yHR2eoeGFY0G98xrssJRQY43Wu/SSfPpInsiMVMnSFvXh/vizIYvX4dQLeT2hUOv8sK+72Fx3RenfpHp4tKKpvPpje4iXVtaydYdLDJd2GiRaW+EnklBbnpUFpkWEQV/+yj9xrvksmi2dwmmC0G3gXDENV5fnF4j4qIvTl3IUdJwfXpFi6tfdrTIdM9sL8TrF5luetVLJj1z0rUohkiMUPDvru++DF9jP9tbuASgx77eIuTDxkP+8JgK+8aLTLcYrTfMp+94kekD+3ThlOEZLdoDxNMi0yKi4N81JZ96H86umA1rP/C2FRwAJ9zstTcOaD3a+kWmt4d4y9YArS0yndUxtWGqpX6R6fzcDHo1ups00RaZFhEF/445B+s+CvfFecZriAbQ51A46Tdeq4Rue/l4+FYWmQ5fAdM46NtaZNpbACOdA/p0aXHDUUGu5tNFkpWCvznn4Nul26dxNq0CDPofCaf8wQv73N57fJgdLTLd+Hlri0z3zE6nICeDgT2yGL13XpOpl/reL5pPF5G2KPjBa4K25t1wL/tnoXQ1WCrsdQwc+WPY9wzo3DPib1e/yHTLFY4qGp6va2U+vfEi0/v36cJJw7bfcFS/XYtMi8ieSt7gr6v1liEsmg1Fc2DrWkjtCAOP9xYu2edU6NStxcvKq2ubrkXaYpHpSjZsbTmfnpmWSmEXb6qlYZHp8CIZ9aP1bp20yLSI+C+5gr+2Gr5Y4I3sP5kL5RuhQyYMPhE3dBylfcewtirNC/HlW1lbuqHhxqO14X4vrS0y3aVTWkMnxv1751KQk9mkPUB+jhaZFpHYkfjBX1NBaOV8qj+YRdrnL5JaXUZ1ahYru4xmUdejWFB3IF+thuLlFVTWvNPkpY0XmR7QPYvDB3Zv1GY3s+HDUi0yLSLxJKGD/+37f8GBXz1AJyqpdFnMCR3C83WH8kZof+oq0sPz5h0Z1iujYZHp+rtItci0iCSqhA7+1NxevN/lJFYXnEhV39Hkd8nmutxMbs/NoHuW5tNFJDkldPCPOudaAI4IuA4RkViieQwRkSSj4BcRSTIKfhGRJBNI8JvZKWb2iZmtNLNfBlGDiEiyinrwm1kq8A/gVGAYcJGZDYt2HSIiySqIEf8oYKVzbpVzrhp4HBgfQB0iIkkpiODvDaxu9HxNeFsTZjbZzBab2eKSkpKoFScikuiCCP7W7ppyLTY4N8U5N9I5N7JHjx5RKEtEJDkEcQPXGqBvo+d9gG939IIlS5ZsMLOvdvN4ecCG3Xytn1TXrlFdu0Z17ZpYrQv2rLb+rW0051oMtn1lZh2AT4ExwDfAe8DFzrmPfDreYufcSD++955QXbtGde0a1bVrYrUu8Ke2qI/4nXO1ZnYNMA9IBe7zK/RFRKSlQHr1OOfmAnODOLaISLJLhjt3pwRdQBtU165RXbtGde2aWK0LfKgt6nP8IiISrGQY8YuISCMKfhGRJJMwwb+zxm9mlm5mT4S/vsjMBsRIXZebWYmZLQv/uSoKNd1nZuvNbHkbXzcz+1u45g/MbITfNUVY13FmVtroXP1PlOrqa2avmlmRmX1kZtNbDVIAAAUpSURBVNe2sk/Uz1mEdUX9nJlZhpm9a2bvh+u6tZV9ov5+jLCuqL8fGx071cz+Y2ZzWvla+54v51zc/8G7LPRzYCDQEXgfGNZsnx8C/ww/vhB4Ikbquhz4e5TP1zHACGB5G18/DXge7y7rw4FFMVLXccCcAH6+CoER4cfZePehNP//GPVzFmFdUT9n4XPQOfw4DVgEHN5snyDej5HUFfX3Y6Nj/xR4tLX/X+19vhJlxB9J47fxwLTw4xnAGDPze9HdmGxI55xbCGzawS7jgQed5x2gi5kVxkBdgXDOFTvnloYfbwGKaNlfKurnLMK6oi58DraGn6aF/zS/iiTq78cI6wqEmfUBTgf+3cYu7Xq+EiX4I2n81rCPc64WKAW6x0BdAOeGpwdmmFnfVr4ebZHWHYQjwr+qP29mw6N98PCv2AfjjRYbC/Sc7aAuCOCchactlgHrgZecc22eryi+HyOpC4J5P/4V+DkQauPr7Xq+EiX4I2n8FlFzuHYWyTGfBQY45w4AXmb7v+pBCuJcRWIp0N85dyBwN/B0NA9uZp2BmcB1zrmy5l9u5SVROWc7qSuQc+acq3POHYTXi2uUme3XbJdAzlcEdUX9/WhmZwDrnXNLdrRbK9t2+3wlSvBH0vitYR/z+gXl4v+0wk7rcs5tdM5VhZ/eCxzic02R2OVGetHgnCur/1XdeXd/p5lZXjSObWZpeOH6iHPuqVZ2CeSc7ayuIM9Z+JibgdeAU5p9KYj3407rCuj9OBoYZ2Zf4k0Hn2BmDzfbp13PV6IE/3vAYDPby8w64n34MbvZPrOBSeHHE4BXXPiTkiDrajYPPA5vnjZos4HLwleqHA6UOueKgy7KzArq5zXNbBTez+/GKBzXgKlAkXPuzjZ2i/o5i6SuIM6ZmfUwsy7hx5nAicDHzXaL+vsxkrqCeD86525wzvVxzg3Ay4hXnHOXNtutXc9XIL162ptro/Gbmd0GLHbOzcZ7gzxkZivx/qW8MEbq+omZjQNqw3Vd7nddZvYY3tUeeWa2BrgF74MunHP/xOujdBqwEigHrvC7pgjrmgD8wMxqgQrgwij84w3eiGwi8GF4fhjgRqBfo9qCOGeR1BXEOSsEppm3zGoK8KRzbk7Q78cI64r6+7Etfp4vtWwQEUkyiTLVIyIiEVLwi4gkGQW/iEiSUfCLiCQZBb+ISJJR8Iv4zLwOmS06LooERcEvIpJkFPwiYWZ2abhf+zIz+1e4oddWM/uzmS01s/lm1iO870Fm9k64mdcsM+sa3r63mb0cboq21MwGhb9953DTr4/N7JEodIYVaZOCXwQws6HABcDocBOvOuASIAtY6pwbASzAu5sY4EHgF+FmXh822v4I8I9wU7Qjgfq2DQcD1wHD8NZnGO37X0qkDQnRskGkHYzBa8j1XngwnonXujcEPBHe52HgKTPLBbo45xaEt08DpptZNtDbOTcLwDlXCRD+fu8659aEny8DBgBv+P/XEmlJwS/iMWCac+6GJhvNftVsvx31ONnR9E1Vo8d16L0nAdJUj4hnPjDBzHoCmFk3M+uP9x6ZEN7nYuAN51wp8J2ZHR3ePhFYEO6Fv8bMzgp/j3Qz6xTVv4VIBDTqEAGccyvM7GbgRTNLAWqAHwHbgOFmtgRv1aMLwi+ZBPwzHOyr2N6NcyLwr3BnxRrgvCj+NUQiou6cIjtgZludc52DrkOkPWmqR0QkyWjELyKSZDTiFxFJMgp+EZEko+AXEUkyCn4RkSSj4BcRSTL/H66cbeFIwseEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "train_accuracy = [1, 2, 3, 4, 5]\n",
    "test_accuracy = [1, 2, 4, 8, 16]\n",
    "train_loss = [1, 2, 3, 4, 5]\n",
    "test_loss = [1, 4, 9, 16, 25]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(train_accuracy)\n",
    "plt.plot(test_accuracy)\n",
    "plt.xlabel(\"epoch\")\n",
    "plt.ylabel(\"accuracy\")\n",
    "plt.savefig(\"accuracy.png\", dpi=600)\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(train_loss)\n",
    "plt.plot(test_loss)\n",
    "plt.xlabel(\"epoch\")\n",
    "plt.ylabel(\"loss\")\n",
    "plt.savefig(\"loss.png\", dpi=600)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('', '-41.312', 'K')\n",
      "N terminal mod\n",
      "('C', '+57.02', 'H')\n",
      "Internal mod\n",
      "('C', '+57.02', 'E')\n",
      "Internal mod\n",
      "('K', '41', '')\n",
      "C terminal mod\n"
     ]
    }
   ],
   "source": [
    "pep = \"-41.312KYSDASDC+57.02HGEDSQAFC+57.02EK41\"\n",
    "rex = r\"([A-Z]?)([-+]?\\d*\\.\\d+|\\d+)([A-Z]?)\"\n",
    "A = re.findall(rex, pep)\n",
    "for a in A:\n",
    "    print(a)\n",
    "    if not a[0]:\n",
    "        print(\"N terminal mod\")\n",
    "    elif not a[2]:\n",
    "        print(\"C terminal mod\")\n",
    "    else:\n",
    "        print(\"Internal mod\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_decoy(pep):\n",
    "    first = \"\"\n",
    "    if pep[0].islower():\n",
    "        first = pep[0]\n",
    "        pep = pep[1:]\n",
    "\n",
    "    pep_parts = re.findall(r\"([A-Z][a-z]?)\", pep)\n",
    "    decoy_pep = pep_parts[0] + \"\".join(pep_parts[-2:0:-1]) + pep_parts[-1]\n",
    "    decoy_pep = first + decoy_pep\n",
    "    if decoy_pep == pep:\n",
    "        decoy_pep = []\n",
    "\n",
    "    return decoy_pep\n",
    "    \n",
    "def pad_left(arr, size):\n",
    "    out = np.zeros(size)\n",
    "    out[-len(arr):] = arr\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n"
     ]
    }
   ],
   "source": [
    "pep = \"aMoQMTPTDNTyNPGSTGScGSSGR\"\n",
    "count = 3\n",
    "mods = []\n",
    "print(mod_filt(pep, mods, count))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "aas            = ['_PAD'] + list(config.AAMass.keys())# + list(config.ModCHAR.values())\n",
    "aa2idx         = {a:i for i, a in enumerate(aas)}\n",
    "idx2aa         = {i:a for i, a in enumerate(aas)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([19, 21, 23], [1.0, 1.0, 1.0])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def decimal_to_binary_array(num, arr_len):\n",
    "    bin_arr = [float(i) for i in list('{0:0b}'.format(num))]\n",
    "    assert len(bin_arr) <= arr_len\n",
    "    res = [0.] * (arr_len - len(bin_arr)) + bin_arr\n",
    "    inds = [int(i) for i, _ in enumerate(res) if res[i] > 0.1] # greater than zero. 0.1 for the floating pointing errors.\n",
    "    vals = [1.0] * len(inds)\n",
    "    return inds, vals\n",
    "decimal_to_binary_array(21, 24)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "if re.search(r\"([a-z]{2,})\", \"MaoEPAMEPETLEAR\"):\n",
    "    print(True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjrElEQVR4nO3debhcVZnv8e+PkBAkCUETBUIG0BAFtAWDCKGZFFsiw0VphVYUZFCQqVtU5PYVtG1tbVt9BAFphgBCFJkEGhAUBIMyBAgJIaB0mEIChMEkDCLB9/6x1sGdoqrOTjh7V86p3+d56qk9r7eqzqm39lp7r6WIwMzMutcanQ7AzMw6y4nAzKzLORGYmXU5JwIzsy7nRGBm1uWcCMzMupwTgXUtSQdImtHpOMpojFXSc5I26aNjHy/pjDw9QVJIWrOPjj0uxzqoL45n1XAi6CKStpf0O0lLJD0j6WZJW1dc5kOSPlBlGd0oIoZFxPx220jaSdKCEsf6ZkQc3BdxNX7eEfFIjvWVvji+VaNPsr6t/iSNAK4EDgMuBIYAfw+81OG41oyI5Z2MoS/019fRX+O2vuUzgu6xKUBETI+IVyLixYi4NiJmw6tVDzdLOimfMdwn6f09O0taV9KZkhZJekzSN4qn+5IOkTRP0jJJ90raStJ5wDjgilw98KVC1cNBkh4Brm/2y7X4y1LSiZJ+Lukn+fhzJG0q6SuSnpT0qKQPtnrhksZKukTSYklPSzq5Yf13JT0r6UFJuxWWH1h4TfMlfbawbidJCyR9WdLjwNmS1pZ0Tj7WvPx6FxT22VDSxTmOByUd1SbmN0m6XNJSSbcBb21YH5Lelqen5vd8Wf5sjpW0DnA1sGF+75/L5Z8o6aL8Xi4FDsjLftIQwmckLcyf9xcK5U6T9I3G9yFPt/u81yy8B5fnM9IHJB1SONaJki6UdG5+LXMlTW71HlnfcSLoHn8AXslfVLtJWq/JNtsA84FRwAnAJZLemNedAywH3gZsCXwQOBhA0j8CJwKfAkYAewJPR8T+wCPAHrl64DuFsnYE3gH8Q8n49wDOA9YD7gJ+Sfr7HQN8Hfhxs51ysroSeBiYkLf/acNrvj+/5u8AZ0pSXvcksHt+TQcC35e0VWHf9YE3AuOBQ0nv2QRgE2BX4JOFONYArgDuzjG8HzhGUqvX/yPgz8AGwGfyo5Uzgc9GxHBgC+D6iHge2A1YmN/7YRGxMG+/F3ARMBI4v8UxdwYmkj7n41Sieq+Xz7vHdGABsCGwD/DN4g8O0t/OT3NslwMnNx7AKhAR/e4BnEX6J72n5PYfA+4F5gIXdDr+Dr5v7wCmkf4Rl5P+0d6S1x0ALARU2P42YH/gLaQqpLUL6/YDbsjTvwSOblHmQ8AHCvMTgAA2KSzbCVjQaj9SkrmusG4P4DlgUJ4fno85skn52wKLgTWbrDsAeKAw/4Z8nPVbvJbLel5njvkvwNDC+vnAPxTmD+55XaSE80jD8b4CnN2knEHAy8DbC8u+CcwozAfwtjz9CPBZYETDcZq9rycCNzVZ9pOGz6dY9neAM/P0NOAbrcpo83mvCYwFXgGGF9Z/C5hWiONXhXWbAS92+v+mGx799YxgGvChMhtKmkj6h5sSEZsDx1QX1uotIuZFxAERsRHpl+OGwA8KmzwW+T8wezhvMx4YDCyS9CdJfyL9An9z3m4s8L8rGc6jK7n9E4XpF4Gn4m8NkC/m52FN9hsLPByt68Ef75mIiBeKx8lnTrfkaow/AVNJZw49FkfEnwvzG7Li6ypOjydV0/yp8B4eT0qyjUaTvjiL+z/cIn6Aj+bYHpZ0o6Rt22zbGFeZbXr+Dl6vDYFnImJZw7HHFOYfL0y/AAxVH13BZK31y0QQETcBzxSXSXqrpGsk3SHpt5LenlcdAvwoIp7N+z5Zc7irpYi4j5RQtygsHlOoFoFU37uQ9KXwEjAqIkbmx4icWMnrV6jDLhZVYvnzpF/jwKvVOaPLvpZePAqMW9kvE0lrARcD3yWdNY0ErgKK70/ja1sEbFSYH9sQx4OF929kRAyPiKlNil9MOmMr7j+uVawRcXtE7EVKzJeRLgZoFl+ruJtpLLunWmmFz4pUPVb22AuBN0oa3nDsx0rEYxXql4mghdOBIyPiPcCxwCl5+abApkoNobdIKnUmMdBIerukL0jaKM+PJVXv3FLY7M3AUZIG53r/dwBXRcQi4FrgvySNkLRGTrw75v3OAI6V9B4lb5M0Pq97glRn3s4fSL/8PixpMPCvwFp98bpJ1VuLgP+QtI6koZKmlNhvSI5hMbBcqRG5ZYN0diHwFUnrSRoDHNEQx9LcuLy2pEGStlCTy3fzmc4lwImS3iBpM+DTzQqUNETSJyStGxEvA0tJ1S+Q3vs3SVq3xOtt9P9y2ZuT2kd+lpfPAqZKeqOk9XntGXbLzzsiHgV+B3wrfw7vAg6idTuF1WRAJAJJw4DtgJ9LmkWqttggr16T1Oi1E+mL7wxJI+uPsuOWkeqpb5X0PCkB3AN8obDNraT36ing34F9IuLpvO5TpC/He4FnSY2NGwBExM/z9hfkci4jNaJCqgP+11wdcmyzwCJiCXA4KaE8RvrV2ev172XkL9U9SI3cj+TjfrzEfsuAo0hf7s8C/0RqU2nn6/n4DwK/Ir1HLzXE8e68/inS6231JX0EqYrqcdKZ29ltyt0feChfBfQ5ciN1PuubDszP7//KVO/cCDwA/Br4bkRcm5efR2rwfoj04+BnDfv19nnvR2o3WAhcCpwQEdeVCUjS1ZKOX4nXYCVpxSrh/kPSBODKiNhC6Rr5+yNigybbnQbcEhHT8vyvgeMi4vY6413dSToAODgitu90LAOFpMOAfSNix143NuugAXFGEBFLgQdzdQa5euLv8urLSJfCIWkUqaqo7R2ZZqtC0gaSpuSqs0mks61LOx2XWW/6ZSKQNB34PTBJ6aaeg4BPAAdJupt0meheefNfAk9Luhe4AfhiobrDrC8NIVVLLgOuB37B39qqzFZb/bZqyMzM+ka/PCMwM7O+0+9u1Bg1alRMmDCh02GYmfUrd9xxx1MR0fT+nH6XCCZMmMDMmTM7HYaZWb8iqeXd6a4aMjPrck4EZmZdzonAzKzLORGYmXU5JwIzsy7nRGBm1uWcCMzMupwTgZlZl3MiMDPrcv3uzmJbNXucNKPS419xpIcxMOuvfEZgZtblnAjMzLrcKiUCSXf1dSBmZtYZq5QIImLLvg7EzMw6o9dEIOnbZZaZmVn/VOaMYNcmy3brbSdJYyXdIGmepLmSjm6yzU6SlkialR9fLRO0mZn1nZaXj0o6DDgc2ETS7MKq4cDNJY69HPhCRNwpaThwh6TrIuLehu1+GxG7r2zgZmbWN9rdR3ABcDXwLeC4wvJlEfFMbweOiEXAojy9TNI8YAzQmAjMzKyDWlYNRcSSiHgoIvYDxgK7RMTDwBqSNl6ZQiRNALYEbm2yeltJd0u6WtLmLfY/VNJMSTMXL168MkWbmVkvyjQWnwB8GfhKXjQE+EnZAiQNAy4GjomIpQ2r7wTGR8TfAScBlzU7RkScHhGTI2Ly6NFNx142M7NVVKaxeG9gT+B5gIhYSGon6JWkwaQkcH5EXNK4PiKWRsRzefoqYLCkUSVjNzOzPlAmEfwlIgIIAEnrlDmwJAFnAvMi4nsttlk/b4ek9+Z4ni5zfDMz6xtlOp27UNKPgZGSDgE+A/x3if2mAPsDcyTNysuOB8YBRMRpwD7AYZKWAy8C++akY2ZmNek1EUTEdyXtCiwFJgFfjYjrSuw3A1Av25wMnFwyVjMzq0CviSBXBV0fEddJmgRMkjQ4Il6uPjwzM6tamTaCm4C1JI0BfgUcCEyrMigzM6tPmUSgiHgB+AhwUkTsDWxWbVhmZlaXUolA0rbAJ4D/ycs8spmZ2QBRJhEcTbqZ7NKImCtpE+CGasMyM7O6lLlq6CZSO0HP/HzgqCqDMjOz+pS5amg08CVgc2Boz/KI2KXCuMzMrCZlqobOB+4DNga+BjwE3F5hTGZmVqMyieBNEXEm8HJE3BgRnwHeV3FcZmZWkzJX//TcOLZI0oeBhcBG1YVkZmZ1KpMIviFpXeALpK6iRwDHVBmUmZnVp0wieDYilgBLgJ0BJE2pNCozM6tNmTaCk0ouMzOzfqjd4PXbAtsBoyX9S2HVCGBQ1YGZmVk92lUNDQGG5W2KI5ItJY0jYGZmA0DLRBARNwI3SpqWB603M7MBqExj8QuS/hPfWWxmNiD5zmIzsy7nO4vNzLqc7yw2M+tyq3pn8T9XGpWZmdWmzHgEV+bJV+8sNjOzgaPdDWUnAdFqfUR4cBozswGgXWPxTOAO0iWjWwF/zI93A69UHpmZmdWi3Q1l5wBIOgDYOSJezvOnAdfWEp2ZmVWuzOWjG7JiFxPD8jIzMxsAylw19B/AXZJuyPM7AidWFpGZmdWqzFVDZ0u6GtgmLzouIh6vNiwzM6tLmTMC8hf/LyqOxczMOqBMG4GZmQ1gTgRmZl2u10Qg6Y1NHoNL7DdW0g2S5kmaK+noJttI0g8lPSBptqStVvWFmJnZqilzRnAnsBj4A+mGssXAg5LulPSeNvstB74QEe8g9Vb6eUmbNWyzGzAxPw4FTl3J+M3M7HUqkwiuAaZGxKiIeBPpy/tC4HDglFY7RcSiiLgzTy8D5gFjGjbbCzg3kluAkZI2WIXXYWZmq6hMIpgcEb/smYmIa4Ed8hf3WmUKkTQB2BK4tWHVGODRwvwCXpsskHSopJmSZi5evLhMkWZmVlKZRPCMpC9LGp8fXwKelTQI+GtvO0saBlwMHBMRSxtXN9nlNR3dRcTpETE5IiaPHj26RMhmZlZWmUTwT6SBaC4j3UswLi8bBHys3Y65Ufli4PyIuKTJJguAsYX5jUgD35iZWU3K3Fn8FHBki9UPtNpPkoAzgXkR8b0Wm10OHCHpp6Q7l5dExKLeYjIzs77TayKQtClwLDChuH1E7NLLrlOA/YE5kmblZceTziiIiNOAq4CppITyAnDgSkVvZmavW5kuJn4OnAacwUqMQxARM2jeBlDcJoDPlz2mmZn1vTKJYHlE+Pp+M7MBqkxj8RWSDpe0QfHu4sojMzOzWpQ5I/h0fv5iYVkAm/R9OGZmVrcyVw1tXEcgZmbWGS0TgaRdIuJ6SR9ptr7FfQFmZtbPtDsj2BG4HtijyboAnAjMzAaAlokgIk7Iz76238xsACtzQ9lawEd57Q1lX68uLDMzq0uZq4Z+ASwB7gBeqjYcMzOrW5lEsFFEfKjySMzMrCPK3FD2O0nvrDwSMzPriDJnBNsDB0h6kFQ1JFI3Qe+qNDIzM6tFmUSwW+VRmJlZx7S7oWxEHlFsWY3xmJlZzdqdEVwA7E66WihYsUtp9zVkZjZAtLuhbPf87L6GzMwGsDJtBEhaD5gIDO1ZFhE3VRWUmZnVp8ydxQcDR5MGlp8FvA/4PdDbUJVmZtYPlLmP4Ghga+DhiNgZ2BJYXGlUZmZWmzKJ4M8R8WdI/Q5FxH3ApGrDMjOzupRpI1ggaSRwGXCdpGeBhVUGZWZm9SkzQtneefJESTcA6wLXVBqVmZnVpm0ikLQGMDsitgCIiBtricrMzGrTNhFExF8l3S1pXEQ8UldQNrDscdKMSo9/xZHbV3p8s4GuTBvBBsBcSbcBz/csjIg9K4vKzMxqUyYRfK3yKMzMrGPKJIKpEfHl4gJJ3wbcXmBmNgCUuY9g1ybL3DW1mdkA0a4b6sOAw4FNJM0urBoO3Fx1YGZmVo/euqG+GvgWcFxh+bKIeKbSqMzMrDYtq4YiYklEPBQR+0XEw4VHqSQg6SxJT0q6p8X6nSQtkTQrP766qi/CzMxWXaluqFfRNOBk4Nw22/y2Z9wDMzPrjJZnBJLWej0HzuMVuArJzGw11+6qod8DSDqvwvK3zXcuXy1p81YbSTpU0kxJMxcvdg/YZmZ9qV3V0BBJnwa2k/SRxpURccnrLPtOYHxEPCdpKql304nNNoyI04HTASZPnhyvs1wzMytolwg+B3wCGAns0bAugNeVCCJiaWH6KkmnSBoVEU+9nuOamdnKaTd4/QxghqSZEXFmXxcsaX3giYgISe8lVVM93dflmJlZe2WuGjpP0lHADnn+RuC0iHi53U6SpgM7AaMkLQBOAAYDRMRpwD7AYZKWAy8C+0aEq33MzGpWJhGcQvoCPyXP7w+cChzcbqeI2K+X9SeTLi81M7MOKpMIto6IvyvMXy/p7qoCMjOzepXpdO4VSW/tmZG0CfBKdSGZmVmdypwRfBG4QdJ8QMB44MBKozIzs9qUGbz+15ImApNIieC+iHip8sjMzKwWpfoayl/8s3vd0MzM+p0ybQRmZjaAORGYmXW5UlVDksaQGolf3T73LmpmZv1cr4kgD1T/ceBe/nbZaABOBGZmA0CZM4L/A0zylUJmZgNTmTaC+eQ+gszMbOApc0bwAjBL0q+BV88KIuKoyqIyM7PalEkEl+eHmZkNQGXuLD5H0hBg07zo/t66oDYzs/6jzFVDOwHnAA+RupgYK+nTvnzUzGxgKFM19F/AByPifgBJmwLTgfdUGZiZmdWjzFVDg3uSAEBE/AFfRWRmNmCUOSOYKelM4Lw8/wngjupCMjOzOpVJBIcBnweOIrUR3MTfhq00M7N+rsxVQy8B38sPMzMbYFomAkkXRsTHJM0h9S20goh4V6WRmZlZLdqdERydn3evIxAzM+uMllcNRcSiPHl4RDxcfACH1xOemZlVrczlo7s2WbZbXwdiZmad0a6N4DDSL/+3SiqOVzwcuLnqwMzMrB7t2gguAK4GvgUcV1i+LCKeqTQqMzOrTctEEBFLgCWSvtywapikYRHxSLWhmZlZHcrcUPY/pMtHBQwFNgbuBzavMC4zM6tJmRvK3lmcl7QV8NnKIjIzs1qVuWpoBRFxJ7B1BbGYmVkHlBmP4F8Ks2sAWwGLK4vIzMxqVeaMYHjhsRapzWCv3naSdJakJyXd02K9JP1Q0gOSZucqJzMzq1mZNoKvAUgakWZjWcljTwNOBs5tsX43YGJ+bAOcmp/NzKxGvZ4RSJqcO56bDcyRdLekXkcny0NZtrvfYC/g3EhuAUZK2qBs4GZm1jfKVA2dRepvaEJETCCNTXB2H5Q9Bni0ML8gL3sNSYdKmilp5uLFbp4wM+tLZRLBsoj4bc9MRMwAylYPtaMmy17T3XUu8/SImBwRk0ePHt0HRZuZWY92fQ31NN7eJunHpAHrA/g48Js+KHsBMLYwvxGwsA+Oa2ZmK6FdY/F/NcyfUJhu+st9JV0OHCHpp6RG4iWFrq/NzKwm7foa2vn1HFjSdGAnYJSkBaREMjgf+zTgKmAq8ADwAnDg6ynPzMxWTbuqoU9GxE8abih7VUS0HcM4IvbrZX2QGp7NzKyD2lUNrZOfh9cRiJmZdUa7qqEfSxoELI2I79cYk5mZ1ajt5aMR8QqwZ02xmJlZB5QZj+B3kk4GfgY837Mw90JqZmb9XJlEsF1+/nphWQC79H04ZmZWtzKJ4KCImF9cIGmTiuIxM7Oaleli4qImy37e14GYmVlntLuP4O2kcYnXlfSRwqoRpLGLzcxsAGhXNTQJ2B0YCexRWL4MOKTCmMzMrEbt7iP4BfALSdtGxO9rjMnMzGpUpo1gb0kjJA2W9GtJT0n6ZOWRmZlZLcokgg9GxFJSNdECYFPgi5VGZWZmtSmTCAbn56nA9IhoN/ykmZn1M2XuI7hC0n3Ai8DhkkYDf642LDMzq0uvZwQRcRywLTA5Il4mdTOxV9WBmZlZPdrdR7BLRFxfvIdAWmGY4UuqDMzMzOrRrmpoR+B6VryHoEfgRGBmNiC0u4/ghPzsISTNzAawdlVDTYeo7NHbUJVmZtY/tKsa6hmichKwNXB5nt8DuKnKoMzMrD7tqoa+BiDpWmCriFiW50/EvY+amQ0YZW4oGwf8pTD/F2BCJdGYmVntytxQdh5wm6RLSVcL7Q2cU2lUZmZWm14TQUT8u6Srgb/Piw6MiLuqDcvMzOpS5oygZ6B6D1ZvZjYAlWkjMDOzAcyJwMysy5VKBJLGS/pAnl5b0vDe9jEzs/6h10Qg6RDgIuDHedFGwGUVxmRmZjUqc0bweWAKsBQgIv4IvLnKoMzMrD5lEsFLEfHqDWWS1iTdT9ArSR+SdL+kByQd12T9TpKWSJqVH18tH7qZmfWFMpeP3ijpeGBtSbsChwNX9LaTpEHAj4BdSWMd3y7p8oi4t2HT30bE7isZt5mZ9ZEyZwTHAYuBOcBngauAfy2x33uBByJifj6j+Cke2czMbLVT5s7ivwL/nR8rYwzwaGF+AbBNk+22lXQ3sBA4NiLmNm4g6VDgUIBx48atZBhmZtZOr4lA0hxe2yawBJgJfCMinm61a5Nljce5ExgfEc9Jmkq6Gmnia3aKOB04HWDy5Mml2ifMzKycMm0EVwOvABfk+X3z81JgGs2HsoR0BjC2ML8R6Vf/qyJiaWH6KkmnSBoVEU+ViMvMzPpAmUQwJSKmFObnSLo5IqZI+mSb/W4HJkraGHiMlED+qbiBpPWBJyIiJL2X1GbR6gzDzMwqUCYRDJO0TUTcCpC/sIfldctb7RQRyyUdAfwSGAScFRFzJX0urz8N2Ac4TNJy4EVg34hw1Y+ZWY3KJIKDgLMl9Xz5LwMOkrQO8K12O0bEVaSrjIrLTitMnwycvFIRm5lZn2qbCPK9AH8fEe+UtC6giPhTYZMLqwzOzMyq1zYRRMQrkvYCvh8RS2qKyazP7HHSjEqPf8WR21d6fLM6lKkaulnSycDPgOd7FubBaszMrJ8rkwi2y89fLywLYJe+D8fMzOpW5s7inesIxMzMOqPUmMWSPgxsDgztWRYRX2+9h5mZ9RdlBqY5Dfg4cCSp24h/BMZXHJeZmdWkTO+j20XEp4BnI+JrwLas2HWEmZn1Y2USwYv5+QVJGwIvAxtXF5KZmdWpTBvBlZJGAv9J6i00gDOqDMrMzOpT5qqhf8uTF0u6Ehjqm8vMzAaOslcNbQdM6NleEhFxboVxmZlZTcoMTHMe8FZgFmlcAkjVQ04EZmYDQJkzgsnAZu4e2sxsYCpz1dA9wPpVB2JmZp3R8oxA0hWkKqDhwL2SbgNe6lkfEXtWH56ZmVWtXdXQd2uLwszMOqZdIngMeEtE3FxcKGmHvM7MzAaAdm0EPyANS9nohbzOzMwGgHaJYEJEzG5cGBEzSfcUmJnZANAuEQxts27tvg7EzMw6o10iuF3SIY0LJR0E3FFdSGZmVqd2jcXHAJdK+gR/++KfDAwB9q44LjMzq0nLRBARTwDbSdoZ2CIv/p+IuL6WyMzMrBZleh+9AbihhljMzKwDynQxYWZmA5gTgZlZl3MiMDPrcqUGpjGzlbfHSTMqL+OKI7evvAwb+HxGYGbW5ZwIzMy6XKWJQNKHJN0v6QFJxzVZL0k/zOtnS9qqynjMzOy1KmsjkDQI+BGwK7CA1GXF5RFxb2Gz3YCJ+bENcGp+HnBcX2xmq6sqG4vfCzwQEfMBJP0U2AsoJoK9gHPzeMi3SBopaYOIWFRhXGYDnn942MqoMhGMAR4tzC/gtb/2m20zBlghEUg6FDg0zz4n6f6+DbWtUcBTNZa3ymXrKJfdTWX3cfn9quw+1E1lj2+1ospEoCbLYhW2ISJOB07vi6BWlqSZETHZZbtsl+2yB0rZjapsLF4AjC3MbwQsXIVtzMysQlUmgtuBiZI2ljQE2Be4vGGby4FP5auH3gcscfuAmVm9Kqsaiojlko4AfgkMAs6KiLmSPpfXnwZcBUwFHiCNhXxgVfG8Dh2pknLZLttlu+y6KF2wY2Zm3cp3FpuZdTknAjOzLudE0EJv3WNUXPZZkp6UdE+d5eayx0q6QdI8SXMlHV1j2UMl3Sbp7lz21+oqO5c/SNJdkq6ss9xc9kOS5kiaJWlmzWWPlHSRpPvy575tTeVOyq+357FU0jF1lJ3L/+f8d3aPpOmShtZY9tG53Ll1vuaWIsKPhgepcft/gU2AIcDdwGY1lr8DsBVwTwde+wbAVnl6OPCHul476b6SYXl6MHAr8L4aX/u/ABcAV3bgfX8IGFV3ubnsc4CD8/QQYGQHYhgEPA6Mr6m8McCDwNp5/kLggJrK3gK4B3gD6YKdXwETO/HZ9zx8RtDcq91jRMRfgJ7uMWoRETcBz9RVXkPZiyLizjy9DJhH+qepo+yIiOfy7OD8qOVqBkkbAR8GzqijvNWFpBGkHx5nAkTEXyLiTx0I5f3A/0bEwzWWuSawtqQ1SV/Kdd3D9A7gloh4ISKWAzcCe9dUdlNOBM216vqiq0iaAGxJ+mVeV5mDJM0CngSui4i6yv4B8CXgrzWV1yiAayXdkbtUqcsmwGLg7FwtdoakdWosv8e+wPS6CouIx4DvAo+QurRZEhHX1lT8PcAOkt4k6Q2kS+jH9rJPpZwImivV9cVAJmkYcDFwTEQsravciHglIt5Nusv8vZK2qLpMSbsDT0bEHVWX1caUiNiK1CPv5yXtUFO5a5KqIU+NiC2B54G628SGAHsCP6+xzPVIZ/kbAxsC60j6ZB1lR8Q84NvAdcA1pKrn5XWU3YoTQXNd3fWFpMGkJHB+RFzSiRhy9cRvgA/VUNwUYE9JD5GqAXeR9JMayn1VRCzMz08Cl5KqJ+uwAFhQOPO6iJQY6rQbcGdEPFFjmR8AHoyIxRHxMnAJsF1dhUfEmRGxVUTsQKoG/mNdZTfjRNBcme4xBiRJItUXz4uI79Vc9mhJI/P02qR/1vuqLjcivhIRG0XEBNJnfX1E1PLrEEDSOpKG90wDHyRVH1QuIh4HHpU0KS96Pyt2FV+H/aixWih7BHifpDfkv/n3k9rDaiHpzfl5HPAR6n/9K/Dg9U1Ei+4x6ipf0nRgJ2CUpAXACRFxZk3FTwH2B+bkunqA4yPiqhrK3gA4Jw9qtAZwYUTUfilnB7wFuDR9H7EmcEFEXFNj+UcC5+cfPfOpsauXXEe+K/DZusoEiIhbJV0E3EmqlrmLert8uFjSm4CXgc9HxLM1lv0a7mLCzKzLuWrIzKzLORGYmXU5JwIzsy7nRGBm1uWcCMzMupwTgfULkv5v7qlxdu6pcps+Pv7xfXm8VSj/AEkndzIG615OBLbay90i707qFfVdpBvNHm2/10prmgjyeNqr/f9JvvfCbJWs9n/gZqQbzZ6KiJcAIuKpni4Zcj/+387jGNwm6W15+WhJF0u6PT+m5OXDJJ2d+/6fLemjkv6D1AvlLEnnS5qQ++U/hXTD0VhJPb2iImkfSdPy9DRJpyqN4TBf0o5K40nM69mmkaStJf1OadyF23ruKgY2lHSNpD9K+k5h+1MlzVTDGA35tX9V0gzgHyVNVRpTYIakHyqPq5DvXD4rvw93SaqtJ13rJzrZB7YffpR5AMOAWaSxEU4Bdiysewj4v3n6U+SxBEjjCmyfp8eRusyA1NnXDwr7r5efnyssm0DqhfR9hWXF9fsA0/L0NFL/RCJ1YrYUeCfpR9YdwLsbXkvP3btb5/kRpLuJD8jL1wWGAg8DY/M2b8zPg0j9L72r8Nq/lKeHks6SNs7z0wvvxTeBT+bpkfl9XKfTn6sfq8/DZwS22os0RsF7gENJXSb/TNIBhU2mF557Rtf6AHBy7ibjcmBE/uX9AeBHhWO3urX/4Yi4pWSIV0REAHOAJyJiTkT8FZhLSipFk4BFEXF7Ln9ppD7pAX4dEUsi4s+k/n7G5+Ufk3QnqRuEzYHNCsf7WX5+OzA/Ih7M88W+az4IHJffi9+Qksa4kq/NuoD7GrJ+ISJeIX2J/UbSHODTpF/jsGIX4T3TawDbRsSLxePkDsbK9KvyfGMIhenGIQ1fys9/LUz3zDf+j7Urv7jvK8CakjYGjiWdQTybq5uK5ffE2azr9GKZH42I+9tsY13MZwS22lMa23ZiYdG7SVUnPT5eeP59nr4WOKJwjHe3WL5ennw5d7/dyhOS3pEbjl/PaFL3kdoCts7lD1caIauVEaQv+yWS3kLqsrnVcTdRGkwI/vaeQOo88cicBJG05euI3wYgJwLrD4aReiW9V9JsUtXIiYX1a0m6FTga+Oe87Chgcm4Qvhf4XF7+DWA9pYHD7wZ2zstPB2ZLOr9FDMcBVwLXk0a0WiWRhj79OHBSLv86XnuGUdz+blKV0FzgLODmFtu9CBwOXJMbj58AluTV/0Ya9nO2pHvyvNmr3Puo9WtKg8lMjoinOh1Lp0kaFhHP5V/+PwL+GBHf73RctvrzGYHZwHFIbhCeS7r66MedDcf6C58RmJl1OZ8RmJl1OScCM7Mu50RgZtblnAjMzLqcE4GZWZf7/7ZB6wPDUxdKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_pos = range(0, 10)\n",
    "charges = [0.000000e+00, 2.295100e+04, 2.621179e+06, 1.636579e+06, 4.118380e+05, 7.507700e+04, 1.315300e+04,\n",
    "           1.878000e+03, 2.800000e+02, 1.000000e+00]\n",
    "plt.bar(x_pos, charges, align='center', alpha=0.8)\n",
    "plt.xticks(x_pos, x_pos)\n",
    "plt.ylabel('Charge distribution of training dataset.')\n",
    "plt.xlabel('Spectrum charge')\n",
    "plt.title('Spectrum charge distribution.')\n",
    "\n",
    "# plt.savefig('oxygen_group.png', dpi=600)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}