{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmarking predictors\n",
    "\n",
    "Compares three prediction tools fot MHC-I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys, math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "pd.set_option('display.width', 130)\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_context(\"notebook\", font_scale=1.4)\n",
    "import epitopepredict as ep\n",
    "from epitopepredict import sequtils, base, peptutils, mhclearn\n",
    "from IPython.display import display, HTML, Image\n",
    "from importlib import reload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate_predictor(P, allele):\n",
    "\n",
    "    data = mhclearn.get_evaluation_set1(allele, length=9)\n",
    "    print (len(data))\n",
    "    if len(data) < 200:\n",
    "        return None,None,None\n",
    "    P.predict_peptides(list(data.peptide), alleles=allele, cpus=14)\n",
    "    x = P.get_scores(allele)\n",
    "    #x = P.data\n",
    "    x = data.merge(x,on='peptide') \n",
    "    #print (x[:4])\n",
    "    #x.plot(x='ic50',y='score',kind='scatter',s=20,)\n",
    "    #auc = auc_score(x.log50k_x,x.log50k_y,cutoff=.426)\n",
    "    auc = round(ep.auc_score(x.ic50,x.score,cutoff=500),3)\n",
    "    import scipy\n",
    "    pr = scipy.stats.pearsonr(x.ic50, x.score)[0]\n",
    "    return auc, pr, data\n",
    "\n",
    "reload(base)\n",
    "reload(mhclearn)\n",
    "\n",
    "def run_tests():\n",
    "    preds = [base.get_predictor('basicmhc1'),\n",
    "             base.get_predictor('netmhcpan',scoring='affinity'),\n",
    "             ep.get_predictor('mhcflurry')]\n",
    "    comp=[]\n",
    "    test_alleles = mhclearn.get_allele_names()#[:20]\n",
    "    print (len(test_alleles))\n",
    "    for P in preds:\n",
    "        m=[]\n",
    "        for a in test_alleles:\n",
    "            print (a)\n",
    "            if not a.startswith('HLA'): continue\n",
    "            try:\n",
    "                auc,pr,df = evaluate_predictor(P, a)\n",
    "                if auc==None:\n",
    "                    continue\n",
    "                m.append((a,auc,pr,len(df)))            \n",
    "            except Exception as e:\n",
    "                print (a,e)\n",
    "                pass\n",
    "            print (P, auc, pr)\n",
    "        m=pd.DataFrame(m,columns=['allele','auc','pearson r','size'])\n",
    "        m['name'] = P.name\n",
    "        comp.append(m)\n",
    "    return comp\n",
    "\n",
    "comp = run_tests()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "123\n"
     ]
    }
   ],
   "source": [
    "c=pd.concat(comp)\n",
    "c.to_csv('benchmarks.csv')\n",
    "print (len(c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "40\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABEAAAAG4CAYAAABFBpKcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABK7UlEQVR4nO3de7xcVX3//9c7IYSL1gQCokZBAQWrlXqtrbYKBQ1q1YriPfFW2q+CtFK8VhSsF6rIL7EqVZDYqqDVFosEQQSk1gveqnKRQA0IghJIEARiIJ/fH3sfGIaTc51z5pzJ6/l4zGPOrL3W3p+dM2dYfGZdUlVIkiRJkiQNsjn9DkCSJEmSJGmqmQCRJEmSJEkDzwSIJEmSJEkaeCZAJEmSJEnSwDMBIkmSJEmSBp4JEEmSJEmSNPBMgEiSJEmSpIFnAkSSJEnSFifJq5NUkp9t5ngl+fhmjj2tPf7iYY49KsnJSa5MsiHJTUkuSHJYkm16fR+Sxm6rfgcgSZIkSX3wcmAN8PAkT6iqCyd7wiTLgE8ANwD/CvwM2BZ4CvBB4FHAX032OpImxgSIpBknyfZV9dt+xzGa2RKnJEm6pySLgT8DXkqTmHg5MKkESJIn0iQ/LgSWVNVNHYdXJNkTeNZkriFpcpwCI4kk72qHcT4yyWfboZrrkpyQ5D7D1D8gyflJbmkfZybZp6vOHyT5VJIrktyeZG2SU5I8pKvesvba+yZZnuRXwC3tsa2SvCPJZUluS3Jjku8k+cuuc/xZG89v29hPT/KozdzjI9phqevbup9Kst0Y/o1Obu9jtyRfTvIb4Ctj/keWJEkzyUuBW4EvA6cCByeZO8lzvgsI8NKu5AcAVbW6qo6f5DUkTYIjQCR1OgW4BngbsA/NEM0HAwcOVUjyUuDfgLOBtwLz23oXtMNHL22r7g88Avg08Etgd+CvgScmeVRV3dp17RXAOuAfgfu1ZUcBbwdOBL4LbA/8IfBE4EttPE8HzgJ+TtPx2AZ4PfDNNp7LhrnH/2tjfyzwWuDXwJvH8O8zp73Wd4G/B+4YQxtJkjTzvBw4rapuS/I54E00fZczJ3KyJNsCfw58o6rW9CxKST1lAkRSp2uAA6uqAJJcC/xDkj+vqq8l2R74CHByVb16qFGSE2nmuL6T5hsVgI9V1Yc6T57ky8A3gb+kSaJ0ugV4WlV1JhWeDZxRVa8bIeYPATcBT66qG9rrnAJcBLwXOKir/g+7Yt8ReA1jS4DMA06vqr8bQ11JkjQDJfkD4NE0X/hQVd9PspomKTKhBAiwJ00/4cc9CVLSlHAKjKROHxlKfrSWt8/Pbp/3BxYCn02yaOgBzAUuAJ4+1LBzhEeS+7SJhsuA9cDjhrn2J7qSH9AkNn4/ycOHCzbJA2hGhKwcSn60115NM6T1mcMMZ/1E1+sLgB2T/N5w1xjGR8dYT5IkzUwvB24EvtpR9jngee2XPRMx1I+4eTKBSZpaJkAkdVrd+aKq1tJMS9mtLRpKRJwNXN/1+Etg56G2SRa2a4jcQNMZWNvWW8DdU1w6XTFM2Tvbuj9LclGS45I8vuP4ru3zcNvXXUIzZWZRV/lVXa/Xtc8LhzlHt000q8VLkqRZKMkc4CXA+cCuSfZIsgd3T7V93jhPOfTF0W/a5/v2Ik5JU8MpMJLGYyhpuoxmusxIPg/8Cc0UlR/SJEGKZg2O4ZKvt3UXVNU3kuwOPAc4AHglcHiSt1TVsRO5AeDOzZRnDG03DjNKRZIkzR5PAxa3j+cPc/zlwGfanzfQbGE7nKEF1G9vny8HNtJMrZE0Q5kAkdRpT5ppKgC001sWcveoh6FRGtdX1dc2d5IkC2kWAntXVb27o3wbxjbS4i5VtY5mIdVPtwuMnQG8O8mHgCvbao8YpulewG9pRp5IkiRBk+BYC/zNMMeeASxLsnNV/ZqmnzFcHwOafgZtHarq1iTnAPsn2bWqrtxMO0l95BQYSZ3ekKRzJMRh7fPQdq9fpVnD421Jtu5unGSn9sehURbdoyr+lnF87rTrhtylqm4DLqXZ6WXbqroW+AHwyiQ7dLTbHfgLYFVVbW7EhyRJ2oK0X8S8gGaB9X/vftCMWt0KeHHb5AzgCUme2HWe7YFXA78AftJx6Oj2+d+GW1ssye5J3tjbu5I0Ho4AkdTpQcAZSU4HHgO8Djirqs4GqKrfJPlrmqGhP2y3jfsV8BDgmTQ7ryxr650HHNkmSq4EngL8GXADY3dJkm8AF9J8W/MYmm1rT6+qW9o6R9BsTfutJJ/g7m1wb6fZQleSJAmaL0d+j2ah9Hupqks7doNZDrwfeCFwXtvHuBi4f3t8d+CFnV+0VNW3khwCnECzftm/0qxTti3wxzTJl5On5tYkjYUJEEmdXkKTNHgvzXodnwTe1Fmhqk5N8kuarePeRJNw+CXN9rYndFR9KfD/AYfQbAv3DWBfYLNTZ4ZxPE1nZV+azsMvaDojH+iI59wk+9N863I0cAfNzi5vqarLuk8oSZK2WC8HfkfzxcnmnAYckeThVXVZO/rjKJr1Qv6GZnrtd4FDqurc7sZVdWKSC2n6SC+mSZjcTrM97t/R9K0k9UnuueOlpC1RknfR/Mf9AVV1XZ/DkSRJkqSecw0QSZIkSZI08EyASJIkSZKkgWcCRJIkSZIkDTzXAJEkSZIkSQNvi90FJsl84AnAtcCdo1SXJEnDmws8ALiwqjb0O5jZyn6JJEk9MWK/ZItNgNB0Mi7odxCSJA2IpwL/3e8gZjH7JZIk9c6w/ZItOQFyLcAFF1zA4sWL+x2LJEmz0tVXX81Tn/pUaP+7qgmzXyJJ0iSN1i/ZkhMgdwIsXryY3Xbbrc+hSJI06zltY3Lsl0iS1DvD9kvcBUaSJEmSJA08EyCSJEmSJGngmQCRJEmSJEkDzwSIJEmSJEkaeCZAJEmSJEnSwDMBIkmSJEmSBp4JEEmSJEmSNPBMgEiSJEmSpIFnAkSSJA20JHsk+XiSHyW5I8lPx9H2lUkuTXJ7kouSHDxMnXlJ3pfkl0luTXJ+kn16ehOSJGnSTIBIkqRB9/vAs4DLgYvH2ijJQcBK4D+AJcDXgM8lWdJV9cPA64GjgOcCvwPOSfLAyYcuSZJ6xQSIJEkadP9VVQ+uqoOAH4yj3THAF6rqrVV1blW9kSYJ8u6hCkkeBPw18Jaq+kRVnQ38JVDA4T27A0mSNGkmQCRJ0kCrqk3jbZPkocBewCldhz4LPCHJTu3rA4C5wKkd17sZOB04cEIBS5KkKbFVvwOQJEmagfZun7unzFzUPj8CuL6t96uqumGYei9NMqc7AZNkAbCgq/7iyQbcL8uXL+fyyy/vy7WvvvpqABYv7t8/3x577MFhhx3Wt+tLksbOBIjG5cgjj+S6665jl1124dhjj+13OJIkTZWF7fP6rvJ17fMOHfW66wzVmwfcB/hN17HDadYL0STddttt/Q5BkjSLmADRuFx33XVcc801/Q5DkqTZ7Hjg5K6yxcAF0x5JD/Rz9MPQtZcvX963GDS9tuQRR442kibPBIgkSdK9DY30WABc11E+NDLkxo56C4ZpvxDYCNzSfaCq1tM1aiTJROOUNE0ccSTNfiZAJEmS7u2S9nlv4NKO8ke2zz/rqLdzkh2q6sauepdNZAFWSZvniCNNJ0ccDd6II3eBkSRJ6lJVP6dJfBzcdeglwIVVdX37+ixgE/CioQpJ7gM8BzhjGkKVJA2g2267zVFHU8ARIJIkaaAl2Y67t6TdFfi9JAe1ry+sqiuTnAgsrarOvtE7gVOTXAGcDTyXZtvbZw1VqKprknwc+ECSO4ArgSOA0Kz1IUmapRxxNHhMgEiSpEG3M/CFrrKh16+iWZB0bvu4S1V9oU2evI0mqXEF8NKqWtV1rr+lWevjPcD9gAuBP6+qX/bwHiRJ0iSZAJEkSQOtqtbQjMgYqc4yYNkw5SuBlaO03Qi8pX1IkqQZyjVAJEmSJEnSwDMBIkmSJEmSBp4JEEmSJEmSNPBMgEiSJEmSpIFnAkSSJEmSJA08EyCSJEmSJGngmQCRJEmSJEkDzwSIJEmSJEkaeCZAJEmSJEnSwDMBIkmSJEmSBp4JEEmSJEmSNPBMgEiSJEmSpIFnAkSSJEmSJA08EyCSJEmSJGngmQCRJEmSJEkDr68JkCR7JjkzyS1Jrk+yIsl2Y2h3XpIa5vH46YhbkiRJkiTNLlv168JJFgDnAlcCBwE7A8cBOwEvHsMpvgkc0VV2SQ9DlCRJkiRJA6JvCRDgEGAhsE9VrQVIcgfwmSTHVNVFo7RfX1XfnuogJUmSJEnS7NfPKTAHAucMJT9aXwQ2AEv6E5IkSZIkSRpE/UyA7A1c3FlQVRuAK4C9xtD+z9q1Q25P8t9J9ttcxSQLkuzW+QAWTyZ4SZIkSZI0e/RzCsxCYP0w5euAHUZpez7wr8Bq4P7A4cBZSfavqq8PU/9w4KiJBipJkiRJkma3fiZAJqyq7pHMSPJl4H+BdwHDJUCOB07uKlsMXND76CRJkiRJ0kzTzwTIOmDBMOULgUvHc6Kq2pDkNOANmzm+nq7RJknGcwlJkiRJkjSL9XMNkEto1gG5S5L5wO6MMwEiSZIkSZI0kn4mQM4A9kuyY0fZ84H57bExaxMnzwMu7Fl0kiRJkiRpYPQzAXICzbSU05I8I8krgBXAqVV11+4wSU5MckfH66cm+XKSVyV5epKX0CyK+jCaNUAkSZIkSZLuoW9rgFTV+iT7AsuBLwG3AacAR3ZVnds+hlwLbA28F9gRuBX4NvC0qvrmVMctSZIkSZJmn77uAlNVlwHPHKXOMmBZx+vLR2sjSZIkSZLUqZ9TYCRJkiRJkqaFCRBJkiRJkjTwTIBIkqSBl2TPJGcmuSXJ9UlWJNluDO3mJfnHJL9IcnuS/03ygmHqrUlSwzwWTc0dSZKk8errGiCSJElTLckC4FzgSuAgYGfgOGAn4MWjND8BeBHwduBSYCnwhSTPqqpVXXX/HfhQV9n6ycQuSZJ6xwSIJEkadIcAC4F9qmotQJI7gM8kOaaqLhquUZJdaRZif2NVrWjLzgIeSbMbXXcC5FdV9e2puQVJkjRZToGRJEmD7kDgnKHkR+uLwAZgyQjtngAEOGuooKoKOBvYJ8mDpyBWSZI0RUyASJKkQbc3cHFnQVVtAK4A9hqh3ab2eUNX+e/a50d2lb+sXSfkt0m+muSxEw1YkiT1nlNgJEnSoFvI8GtxrAN2GKHdZe3zE4E1HeVPbJ87234Z+A5wFbAr8FbggiRPqKp7JF/aNUkWdF1r8QhxSJKkHjABIkmSNIyq+mmSC4APJLmauxdBfXpbZVNH3cM6ml6QZFVb/y3AK7tOfThw1FTFLUmShmcCZBIe9/ef7ncI0+6+a29mLnDV2pu3uPv//j91918lSbPEOu494gKakSGXjtJ2KfB54Jvt6zXAu4GjgWs316iqbkjydeBxwxw+Hji5q2wxcMEosUiSpEkwASJJkgbdJTTrgNwlyXxgd+BTIzWsqp8DT0iyG7AtzbSYv6NZF+SHEwmmqtbTNSUnyUROJUmSxsFFUCVJ0qA7A9gvyY4dZc8H5rfHRlVVa6rqEmBr4DXA56rq5s3VT7II2A+4cMJRS5KknnIEiCRJGnQnAIcCpyU5BtgZOA44tXOB0iQnAkuraquOsjcAv6FZ3HQ3mtEf82nW9hiq8xLg2cAq4Jq23pvbeu+fwvuSJEnjYAJEkiQNtKpan2RfYDnwJeA24BTgyK6qc9tHp/k0C5Yuppm2cjrw9qr6VUednwMPpEmqLARuAs4HDqqq0dYYkSRJ08QEiCRJGnhVdRnwzFHqLAOWdZV9CPjQKO2+zd07w0iSpBnKNUAkSZIkSdLAMwEiSZIkSZIGngkQSZIkSZI08FwDRNKMdeSRR3Ldddexyy67cOyxx/Y7HEmSJEmzmAkQSTPWddddxzXXXNPvMCRJkiQNAKfASJIkSZKkgWcCRJIkSZIkDTwTIJIkSZIkaeCZAJEkSZIkSQPPBIgkSZIkSRp4JkAkSZIkSdLAMwEiSZIkSZIGngkQSZIkSZI08LbqdwCSJEmanOXLl3P55Zf3O4xpt3r1agAOO+ywPkfSH3vssccWe++SNBEmQCRJkma5yy+/nB/+5GI2bbdDv0OZVvldAfD9K67rcyTTb86tN/Y7BEmadUyASJIkDYBN2+3A7Y98dr/D0DTZ5uLT+3ZtRxxtmaNuHHGkQWACRJIkSdKYXX755Vz20x/wkPvc2e9QptXWG5vlE29fc2GfI5l+V90yt98hSD1hAkSSJEnSuDzkPnfyjsff0u8wNE3e87379DsEqSfcBUaSJEmSJA08EyCSJEmSJGngmQCRJEmSJEkDzwSIJEmSJEkaeCZAJEmSJEnSwDMBIkmSJEmSBp4JEEmSJEmSNPBMgEiSJEmSpIFnAkSSJEmSJA28rfodgCRJM8GRRx7Jddddxy677MKxxx7b73AkSZLUYyZAJEkCrrvuOq655pp+hyFJkqQp4hQYSZIkSZI08EyASJIkSZKkgWcCRJIkDbwkeyY5M8ktSa5PsiLJdmNoNy/JPyb5RZLbk/xvkhdspu4RSX6e5LYk30uyX+/vRJIkTZQJEEmSNNCSLADOBe4LHAS8CXgJcNIYmp8AvBH4IPBc4CLgC0mWdF3jCOC9wD8DzwJWA19J8pje3IUkSZosF0GVJEmD7hBgIbBPVa0FSHIH8Jkkx1TVRcM1SrIrsAx4Y1WtaMvOAh5Jk+xY1ZbNB94BHF9VH2zLzgd+ArwdeNHU3ZokSRorR4BIkqRBdyBwzlDyo/VFYAOwZPgmADwBCHDWUEFVFXA2sE+SB7fFfwzcDzilo96dwOeBJUnSi5uQJEmT09cEyETn43ad4/lJKslPpypOSZI0q+0NXNxZUFUbgCuAvUZot6l93tBV/rv2+ZEd5we4pKveRcB9gAd1FiZZkGS3zgeweLSbkCRJk9O3KTAd83GvpJmPuzNwHLAT8OIxnmM74HjgV1MSpCRJGgQLgfXDlK8Ddhih3WXt8xOBNR3lT2yfh9ouBDZU1W3DnH+o3tUd5YcDR40UsCRJ6r1+rgEyofm4Xf4B+D+aJMrjpyxSSZK0xamqnya5APhAkquBS4GlwNPbKps223hkxwMnd5UtBi6Y4PkkaWAtX76cyy+/vN9hTLvVq1cDcNhhh/U5kv7YY489puTe+5kA2dx83JNo5uOOmABJshdwGPAk4IipClKSJM1664AFw5QvpElqjGQpzVoe32xfrwHeDRwNXNtx/vlJtqmq27vOD3Bj5wmraj1dI1JcJkSShnf55Zfzw4t+OPyn+CBrU+w/vOaH/Y2jH9ZP3an7mQDZm67t56pqQ5LR5uMO+Wfgk+23MyNWbKfbLOgqdq6tZpWrjn50v0OYdnfcuAOwFXfceOUWd/8PeedP+h2CNEgu4e51OoC7dm7ZHfjUSA2r6ufAE9p1OralmRbzdzTrggz1SofW/ti7owyaNUJuBq6ZXPiStIVbAJueNtFBd5pt5pw3dUuV9jMBMtH5uCR5MfBo4AVjvNbhONe2JzZtvf09niVJmgXOAP4hyY5VdUNb9nxgfntsVFW1BiDJtsBrgM9V1c3t4f8BbgIOpk2AJJlLs/3tme3OMdLAuPrqq/ntzXN5z/fu0+9QNE2uvHku21999egVpRmunwmQCUlyX+BDwNvaIaRjcTzOte2J3+55QL9DkCRpvE4ADgVOS3IMdy+8fmpV3bU7TJITgaVVtVVH2RuA3wBXAbvRjP6YD7xlqE47gvU9wHuTXA/8AHgtzQiTl07trUmSpLHqZwJkovNx304zl/ZL7dQWgK2BOe3r29qt7e7iXFtJkrZcVbU+yb7AcuBLwG3AKcCRXVXnto9O82lGkS6m6UucDry9qu6xA11VfbDtWxwG3J9mLbNnVdX/9vRmpBlg8eLF3H7Htbzj8bf0OxRNk/d87z5ss9gVBDT79TMBMtH5uHsBjwJuGObYOuBvaUZ8SJIkAVBVlwHPHKXOMmBZV9mHaEaejuUaHwQ+OLEIJ+fqq69mzq03sc3Fp/fj8uqDObfewNVX39HvMCRpVulnAmSi83Hfwb0THG8BHgG8Ctjy9kiSJEmSJEkj6mcCZELzcavqp90nSrIMWFxV501D3JIkSTPK4sWL+dWGrbj9kc/udyiaJttcfDqLF+/S7zAkaVaZuv1lRtGuy7EvcAvNfNwPA6cCr+6qOtx8XEmSJEmSpDHr6y4wE52Pu5k6kiRJkiRJw+rbCBBJkiRJkqTpYgJEkiRJkiQNPBMgkiRJkiRp4JkAkSRJkiRJA6+vi6BKkiRJkrQ5V199NdwEc87zu/stxnq4uq6eklP7LpIkSZIkSQPPESCSJEmSpBlp8eLFXJ/r2fS0Tf0ORdNkznlzWPygxVNz7ik5qyRJkiRJ0gxiAkSSJEmSJA08EyCSJEmSJGngmQCRJEmSJEkDzwSIJEmSJEkaeCZAJEmSJEnSwHMbXEnSvfzJij/pdwjTbuv1WzOHOfxi/S+2uPv/5qHf7HcIkiRJU84RIJIkSZIkaeCZAJEkSZIkSQPPBIgkSZIkSRp4JkAkSZIkSdLAMwEiSZIkSZIGngkQSZIkSZI08EyASJIkSZKkgWcCRJIkSZIkDTwTIJIkSZIkaeCZAJEkSZIkSQPPBIgkSZIkSRp4JkAkSZIkSdLAG1MCJMl2Sb6e5FVTHZAkSVKvJdkzyZlJbklyfZIVSbYbQ7vtk7w/yRVJbk2yOsk7k2zdUWe3JDXM46dTe1eSJGk8thpLpaq6NcnjgM9NcTySJEk9lWQBcC5wJXAQsDNwHLAT8OJRmn8MeB7wduCnwBOBY4CFwN921X1be50ht04uckmS1EtjSoC0zgeeCnxiimKRpHtYtM0m4I72WZIm7BCahMU+VbUWIMkdwGeSHFNVFw3XKMlWwAuBY6tqRVt8bpJdgZdy7wTI6qr69pTcgSRJmrTxJEAOBc5K8k8034asqSr/r0TSlDniD9b3OwRJg+FA4Jyh5Efri8BJwBJg2AQIEJq+0k1d5evbY5IkaRYZzyKolwK7An8HrAY2tHNhOx+/nZIoJUmSJm5v4OLOgqraAFwB7LW5RlW1Efg0cGiSJyW5T5KnA68DPjJMk48muSPJ2iSfSrJz725BkiRN1nhGgJwK1FQFIkmSNEUW0oza6LYO2GGUtocAHwc6p7Z8uKqO7ni9gWZ07FntOR9Hs2bIHyV5bFXd1nnCdk2SBV3XWTxKHJIkaZLGnACpqmVTGIckSdJM9D7gWTSjPi4D/gg4Ksl1VXUsQFVdC/y/jjbnJ/k+cB7wEpqpNp0OB46a2rAlSVK38YwAkSRJmo3Wce8RF9CMDLl0c42SPAo4AnhuVX25Lf5GknnA0Uk+VlU3D9e2qs5P8mua0SDdCZDjgZO7yhYDF4x8G5IkaTJMgEiSpEF3Cc06IHdJMh/YHfjUCO0e2T7/qKv8h8B8mqTFJeMNpqrW0zUlJ3FNVUmSptp4FkGVJEmajc4A9kuyY0fZ82mSGGeM0O7K9vlxXeWPo1kX7Uo2o10sdWfgwnFHK0mSpoQjQCRJ0qA7ATgUOC3JMTSJieOAU6vqrt1hkpwILK2qof7R94DvAh9vd3RZDTwJeCtwUlXd2rb7ELCJZqHUdcDj2zo/BU6Z+tuTJEljYQJEkiQNtKpan2RfYDnwJeA2msTEkV1V57aPoXZ3JnkOcAxNQuP+wC+AD9IsjjrkYppFUF8HbA/8Evg34J1VdftU3JMkSRo/EyCSJGngVdVlwDNHqbMMWNZV9muarXBHancicOLkIpy8ObfeyDYXn97vMKZVbv8NALXN7/U5kuk359YbgV36HYYkzSoTSoAkuQ/Nyun3WrGrqq6abFCSJEkauz322KPfIfTF6tXNJjx77r4lJgJ22WJ/75I0UWNOgCTZhmbP+tcAO45Qde4IxyRJktRjhx12WL9D6Iuh+16+fHmfI5EkzQbjGQHyUWAp8J80+9Svm4qAJEmSJEmSem08CZC/BD5ZVSPOg5UkSZIkSZpp5oyjbgE/mKpAJEmSJEmSpsp4EiCnAX8+VYFIkiRJkiRNlfEkQN4L7J7kE0melOQBSXbufkxVoJIkSZIkSRM1njVALm2f9wFePUI9d4GRJM06tV2xiU3UdtXvUCRJkjQFxpMAOZpmHRBJkgbOxj/Z2O8QJEmSNIXGnACpqndNYRySJEmSJElTZjwjQO6SJMCi9uXaqnJkiCRJkiRJmrHGswgqSfZI8nngJuC69nFTklOS7DHeiyfZM8mZSW5Jcn2SFUm2G0O7f05ySZKbk/wmyXeTvHi815ckSZIkSVuGMY8ASfL7wDeBbYEvA5e0h/YGngcckOSpVXXRGM+3ADgXuBI4CNgZOA7YCRgtmbEd8DHgZ0CAFwKfSzKnqj471nuSJEmSJElbhvFMgXk/cCvw+Kq6vPNAkt2BC4D3AX8xxvMdAiwE9qmqte157gA+k+SYkRIpVfWqrqIzk+wNLANMgEiSJEmSpHsYzxSYpwL/3J38AKiqK4CPAn86jvMdCJwzlPxofRHYACwZx3mGrAW2nkA7SZI0wyR5S5JvjnD8giRHTGdMkiRpdhvPCJCtgNtHOH7bOM+3N3BSZ0FVbUhyBbDXaI3bhVjnAvcFngMcALx8M3UXAAu6ihePI1ZJkjS9XgZ8dYTj3wJeAXxwesKR1OmqW+bynu/dp99hTKtf3dp8d3z/7Tb1OZLpd9Utc3l4v4OQemA8CYvvA69LcmJVre880CYYXgd8bxznWwisH6Z8HbDDGNo/F/iP9uc7gDdU1b9vpu7hwFHjiE2SJPXXw2jW+tqc1cDfTFMskjrssce49z4YCL9bvRqAbXbbs8+RTL+Hs+X+3jVYxpMAeSdwNnBZkpOBy9ryRwCvpBlh8Ve9DG4U5wFPaK+7BPhIkjuq6sRh6h4PnNxVtphm3RJJkjTz3A48YITjDwTunKZYJHU47LDD+h1CXwzd9/Lly/sciaSJGnMCpKrOT/IM4ENA95zbHwAHV9U3xnHtddx7Wgo0I0MuHUM867l7xMnXkmwNHJfk5Kq6c5i66zvLmhk0kiRphvof4NVJ/r+quqnzQJKFwKvbOpIkSWMynhEgVNW5wGOT7ALs2hZfWVXXTeDal9CsA3KXJPOB3YFPTeB83wfeQLON7kTikSRJM8e7aUZq/jjJ8cDQ7nCPAt4I3B94UX9CkyRNq/Uw57zx7N8xAG5pn7espXYa64EHTc2px5UAGdImPK6DZjHSJNtV1a3jPM0ZwD8k2bGqbmjLng/Mb4+N11OA39DsBiNJkmaxqvpekmcD/0Iz+rTaQwH+D3h2VX2nX/FJkqbHlrr2yOp2zZk9H7TlrTnDg6bu9z7mBEiS5wFPrKq3dZQdQfMNzTZJ/gt46TgSIScAhwKnJTkG2Bk4Dji1qi7uuMaJwNKq2qp9/VSaKTj/AVwJ/B7NLjCvAd5SVXeM9Z4kSdLMVVXnJNkDeCzNCFGAK4AfVFVtvqUkaVC45oxrzvTSeEaAvIVm2goASR4HfAA4n2aV9tcARwLvGsvJqmp9kn2B5cCXaLbRPaU9R6e57WPIL4DfAUNJk3VtXM+rqtPGcT+SJGmGaxMd328fkiRJEzaeBMiewGc7Xr8UuAFYUlUbkmwEXswYEyAAVXUZ8MxR6iwDlnW8XgO8cKzXkCRJs0+SPx1LvXEuwC5JkrZg40mAbAv8tuP1M4Azq2pD+/pHNKNAJEmSJus87l73YyRzR68iSZI0vgTIL4AnACcm2RN4JPD+juOLaKaxSJIkTdbThymbC+wG/BUwh2Z6riRJ0piMJwHyr8C7kzwQ+H3gRuC/Oo4/Abish7FJkqQtVFWdv7ljSU6m2SL3acDXpykkSZI0y41nM+X3Ae8FFgNXAc+vqpsAkuwA/Cnw5Z5HKEmS1KGqNtEsnP7afsciSZJmjzGPAKmqO4F/aB/dx24E7t/DuCRJkkayA7Cg30FIkqTZY0wJkCTbATcD/1BV753akCRJ0pYuyUM2c2gBzajTv6eZBiNJkjQmY0qAVNWtSa4HfjPF8UiSJAGsYfO7wAT4NnDItEUjSZJmvfEsgvp54EVJPtrOvZUkSZoqr+beCZAC1gFXVNXF4zlZu4PdCuApNLvWnQK8uapuHaXd9jTTf18IPAC4hmZh+PdX1e866t0X+CfgIGAb4Fzg0KpaM544JUnS1BlPAuQ/gH2B/07yCeD/GGbb26r6bo9ikyRJW6iqOrlX50qygCYhcSVNgmJn4DhgJ+DFozT/GPA84O3AT4EnAscAC4G/7aj3OeCxwKE0I2aPBs5J8ujRkiySJGl6jCcBck7Hz3/Evb+VSVs2d7JBSZIk9dAhNAmLfapqLUCSO4DPJDmmqi4arlGSrWhGfhxbVSva4nOT7Aq8lDYBkuRJwLOAZ1XVGW3ZT4ArgGXAR6fqxiRJ0tiNJwHyqimLQpIkqUuS+wOvAR4H3A+Y01Wlqmq/MZzqQOCcoeRH64vAScASYNgECM2XO1sBN3WVr2+PdZ7/JuDMjsCuSvLN9pgJEEmSZoDxbIO7cioDkSRJGpLkUcB5wPbAz4BHAxfTjOR4IM3oil+M8XR70yQ77lJVG5JcAey1uUZVtTHJp4FD22TGRcATgNfRrCfSef5Lh1kj7SLgGcPc2wLuvYXv4jHdiSRJmrDxjACRJEmaLu8DbgceD9wC/Bp4Y1V9PclLaBIQo63fMWQhzaiNbuuAHUZpewjwcZpdZ4Z8uKqOnsT5DweOGuW6kiSpx8aVAEmyDfCXjDwU9TU9ik2SJG25ngIcV1VrkgwlEeYAVNXnkjyFZteVfac4jvfRrO/xOuAymnXQjkpyXVUdO8FzHg+c3FW2GLhggueTJEljMOYESJIH06yg/jCabznuB9xI863HHGAtzTc0kiRJk7U18Mv256Fd5xZ0HP8R8Moxnmsd955yAk0f5tLNNWqn4RwBPLeqvtwWfyPJPODoJB+rqpvb8z9kM+e/sbuwqtbTNWIkSXc1SZLUY90jOEZyLLAj8MfAw2kW/zqYZm7u24FbgT/vdYCSJGmLdCVtUqGqbgOuBZ7ccfxRjP2Ll0to1um4S5L5wO6MkAABHtk+/6ir/IfAfO5et+MS4BG5dxbjkaOcX5IkTaPxJED+HPhYVX0bGFrkK1W1oareB3wD+HCvA5QkSVukc4Hndbz+DHBYkk8mOQn4f8BpYzzXGcB+SXbsKHs+TRLjjBHaXdk+P66r/HFAdRw/g2aEyV0LnrYjZ58yyvklSdI0Gs8aIPehWXEdYEP7fN+O4/9NM0pEkiRpsj4AfD3J/KraAPwDzZSSg4A7gX+lmZ4yFicAhwKnJTkG2Bk4Dji1qi4eqpTkRGBpVQ31j74HfBf4eJKdgdXAk4C3AidV1a0AVfWdJF8BTkzyJuA3wNHAVdx7rQ9JktQn40mAXEOz7RxV9dskNwL7AP/ZHt8V2NjL4CRJ0papqq6iSSAMvd5AsxDp6yZwrvVJ9gWWA1+iWVPkFODIrqpz28dQuzuTPAc4hibpcX+arXc/SLM4aqeXtOUfpRlZci7wwqEkiSRJ6r/xJEC+QTO085j29ReBI5LcQTOV5o3A6b0NT5IkafKq6jLgmaPUWQYs6yr7Nc1WuKOd/+a23qh1JUlSf4wnAfJhYP8k21TV7cCbaXaEObo9fh5NEkSSJEmSJGlGGXMCpKp+Avyk4/V6moTIAuDO9psPSZIkSZKkGWc8I0CG1SZCJEmSJEmSZqzxbINLkocn+bck1yT5XbugGEkWJTkpyZOmJkxJkiRJkqSJG3MCJMljgAuB/YFvcc9V0tcCjwL+ptcBSpIkSZIkTdZ4RoC8H7gWeDjw10C6jp8J/EmP4pIkSZIkSeqZ8SRAngL8S1XdBNQwx68CHtiTqCRJkiRJknpoXGuAABtGOHZ/4PZJxCJJkiRJkjQlxpMA+T7w7OEOJJkHvAT4di+CkiRJkiRJ6qXxJEDeCxyQ5BPAY9qyByZ5JnAOzdog7+1xfJIkSZIkSZO21VgrVtVZSV4BrABe3RavpFkMdT3w8qr6Zs8jlCRJkiRJmqQxJ0AAquqzSf4TOADYk2YEyRXAV6vq5t6HJ0mSJEmSNHnjSoAAVNWtwH/2PhRJkiRJkqSpMe4ESJJ9aRZD3a0tWgN8parO6V1YkiRJkiRJvTPmBEiS7YFTgSU0636saw89D3hjkq8CL6qqW3odpCRJkiRJ0mSMZxeYDwEHAu8BdqqqHatqR2An4B+BZwIf7H2IkiRJkiRJkzOeBMiLgE9U1VFVdcNQYVXdUFXvBD7Z1pEkSZIkSZpRxpMAmQP8aITjP6KZGiNJkiRJkjSjjCcBcgbN4qeb8+y2jiRJkiRJ0owynl1gjgFOSXI68BHg8rZ8T+ANwAOBNyXZubNRVf26F4FKkiRJkiRN1HgSIBe1z4+m2Qmm09DUl58O027ueIOSJEmSJEnqpfEkQI4GaqoCkSRJkiRJmipjToBU1bumMA5JkiTNQsuXL+fyyy8fveIUWL16NQCHHXZYX64PsMcee/T1+pKksRvPCBBJkiRpxth22237HYIkaRYxASJJkqQJc/SDJGm2GM82uJIkSZIkSbOSCRBJkiRJkjTwTIBIkiRJkqSBZwJEkiRJkiQNPBMgkiRJkiRp4PU1AZJkzyRnJrklyfVJViTZbpQ2v5fkXUm+k2R92+7MJI+drrglSdLsMsE+x25JaoTHAzrqrtlMnUVTf3eSJGks+rYNbpIFwLnAlcBBwM7AccBOwItHaPoQ4BDgJOCdwDzgjcD/JPnjqvrBFIYtSZJmmUn0Oa4FnjxM+SnAjVV1bVf5vwMf6ipbP/6IJUnSVOhbAoQmibEQ2Keq1gIkuQP4TJJjquqizbT7ObB7Vd06VJDka8D/AYcCr5rasCVJ0iwzoT5HVW0Avt1ZlmRvYFdg+TBNflVV3x6mXJIkzQD9nAJzIHDOUEek9UVgA7Bkc42q6redyY+27HbgEuCBUxGoJEma1SbU59iMlwN3Ap/rUWySJGma9DMBsjdwcWdB+03LFcBe4zlRku2BP6RJggx3fEE7j/euB7B4QlFLkqTZpid9jiQBXgp8fZjpLwAvS3J7kt8m+arrk0mSNLP0cwrMQoafF7sO2GGc53oPsB3wkc0cPxw4apznlCRJg6FXfY6nALsxfJ/iy8B3gKtopsi8FbggyROq6h7Jl3ZNkgVd7f1iRpKkKdbPBEhPJHkpTYLj9VV1+WaqHQ+c3FW2GLhgygKTJEmD5mXArcCXug9U1WEdLy9Isgq4FHgL8Mqu6ofjFzOSJE27fiZA1nHvbz+g+Zbm0rGcIMn+wKeAf6qqj26uXlWtp+ubn2YUqyRJ2gL0os+xNfBC4LSqumW0+lV1Q5KvA48b5vDx+MWMJEnTrp8JkEto5uTeJcl8YHeapMaIkjyR5huYzwNvnooAJUnSQJhUn6N1IM10mX+bbDB+MSNJUn/0cxHUM4D9kuzYUfZ8YH57bLPaLejOAL4JvLqqasqilCRJs92E+xwdXgb8GjhrLJWTLAL2Ay4cR5ySJGkK9TMBcgLNtx+nJXlGklcAK4BTOxcLS3Jikjs6Xu8MfBX4HfBPwOOS/FH7+MNpvQNJkjQbTKjP0VF+P+DZbf3hjr8kyWeSvDzJ05O8CvhvmgTL+6fmliRJ0nj1bQpMVa1Psi+wnGYqy23AKcCRXVXnto8hjwQe3P78ta66V9Kszi5JkgRMqs8x5AXANmx++svPgQcCx9GsK3ITcD5wUFWNaY0RSZI09fq6C0xVXQY8c5Q6y4BlHa/PA5woK0mSxmwifY6O8pOAk0Zo923g6ZOLUJIkTbV+ToGRJEmSJEmaFiZAJEmSJEnSwDMBIkmSJEmSBp4JEEmSJEmSNPBMgEiSJEmSpIFnAkSSJEmSJA08EyCSJEmSJGngmQCRJEmSJEkDzwSIJEmSJEkaeCZAJEmSJEnSwDMBIkmSJEmSBp4JEEmSJEmSNPBMgEiSJEmSpIG3Vb8DkCRJkqSxWL58OZdffnlfrr169WoADjvssL5cf4899ujbtaVBYQJEkiRJkkax7bbb9jsESZNkAkSSJEnSrOAICEmTYQJEkiRJkqQuTrkavISjCRBJkiRJkmYQp1xNDRMgkiRJkiR1GcQREFs6t8GVJEmSJEkDzwSIJEmSJEkaeCZAJEmSJEnSwDMBIkmSJEmSBp4JEEmSJEmSNPBMgEiSJEmSpIFnAkSSJEmSJA08EyCSJEmSJGngmQCRJEmSJEkDzwSIJEmSJEkaeCZAJEmSJEnSwDMBIkmSJEmSBp4JEEmSNPCS7JnkzCS3JLk+yYok243SZrckNcLjAV31j0jy8yS3Jflekv2m9q4kSdJ4bNXvACRJkqZSkgXAucCVwEHAzsBxwE7Ai0doei3w5GHKTwFurKprO65xBPBe4G3AD4DXAV9J8qSq+t8e3IYkSZokEyCSJGnQHQIsBPapqrUASe4APpPkmKq6aLhGVbUB+HZnWZK9gV2B5R1l84F3AMdX1QfbsvOBnwBvB17U8zuSJEnj5hQYSZI06A4EzhlKfrS+CGwAlozzXC8H7gQ+11H2x8D9aEaGAFBVdwKfB5YkyUSCliRJvWUCRJIkDbq9gYs7C9rRHVcAe431JG0i46XA1zunv7TnB7ikq8lFwH2AB403YEmS1HtOgZEkSYNuIbB+mPJ1wA7jOM9TgN2Ao4Y5/4aqum2Y89Ne4+qhwnZNkgVddRePIw5JkjQBJkAkSZLG5mXArcCXJnmew7l3EkWSJE0xp8BIkqRBt457j7iAZuTGjWM5QZKtgRcCp1XVLcOcf36SbYY5P8Nc43jgoV2Pp44lDkmSNHGOAJEkSYPuEu5epwO4a+eW3YFPjfEcB9JMZfm3zZyf9ho/7Ch/JHAzcE1n5apaT9eUHNdJlSRp6jkCRJIkDbozgP2S7NhR9nxgfntsLF4G/Bo4a5hj/wPcBBw8VJBkLs32t2dWVU0kaEmS1FsmQCRJ0qA7gWbExWlJnpHkFcAK4NSqumt3mCQnJrmju3GS+wHPbuvf63i7o8x7gL9L8qYkTwc+TTPC5B+n4obUWLt2LYceeig33HBDv0ORJM0CJkAkSdJAa6ec7AvcQrOA6YeBU4FXd1Wd2z66vQDYhuGnvwxd44PA24DDgFU02+s+q6r+d5LhawQrV67kxz/+MStXrux3KJKkWcAEiCRJGnhVdVlVPbOqtq+qRVX1hqq6tavOsqq612IcVXVSVaWqvjvKNT5YVbtW1TZV9biqOqfX96G7rV27llWrVlFVrFq1ylEgkqRRmQCRJEnSrLNy5UqGllfZtGmTo0AkSaMyASJJkqRZ5+yzz2bjxo0AbNy4kbPOGm59WkmS7mYCRJIkSbPO/vvvz7x58wCYN28eBxxwQJ8jkiTNdCZAJEmSNOssXbqUpFmyZc6cOSxdurTPEUmSZjoTIJIkSZp1Fi1axJIlS0jCkiVL2HHHHfsdkiRphjMBIkmSpFlp6dKl/MEf/IGjPyQNnLVr13LooYe6w1WP9TUBkmTPJGcmuSXJ9UlWJNluDO0OTvLFJFcnqSRHTEe8kiRJmjkWLVrEihUrHP0haeCsXLmSH//4x+5w1WN9S4AkWQCcC9wXOAh4E/AS4KQxND8IeBhw+lTFJ0mSJEnSdFu7di2rVq2iqli1apWjQHqonyNADgEWAs+tqjOr6tPAYcDBSX5/lLYHV9UfVtVfT3mUkiRJkiRNk5UrV1JVAGzatMlRID3UzwTIgcA5VbW2o+yLwAZgyUgNq2rTVAYmSZIkSVI/nH322WzcuBGAjRs3ctZZZ/U5osHRzwTI3sDFnQVVtQG4AtirlxdKsiDJbp0PYHEvryFJkiRJ0mTtv//+zJs3D4B58+ZxwAEH9DmiwdHPBMhCYP0w5euAHXp8rcOBn3c9LujxNSRJkiRJmpSlS5eSBIA5c+a401UPbSnb4B4PPLTr8dR+BiRJkiRJUrdFixaxZMkSkrBkyRJ3uuqhrfp47XXAgmHKFwKX9vJCVbWertEmQxk1SZIkSZJmkqVLl7JmzRpHf/RYP0eAXEKzDshdkswHdqfHCRBJkiRJkmaLRYsWsWLFCkd/9Fg/EyBnAPsl6fyNPh+Y3x6TJEmSJEnqiX4mQE6gmZZyWpJnJHkFsAI4taru2h0myYlJ7uhsmOSRSQ5KclBb9Oiu15IkSZIkSXfp2xogVbU+yb7AcuBLwG3AKcCRXVXnto9OLwKO6nj9yvYB4OIekiRJkiTpHvq6C0xVXVZVz6yq7atqUVW9oapu7aqzrKrSVfauqspwj+m9A0mSJEmSNBtsKdvgSpIkSZKkLZgJEEmSJEmSNPBMgEiSJEmSpIFnAkSSJEmSpBlk7dq1HHroodxwww39DmWgmACRJEmSJGkGWblyJT/+8Y9ZuXJlv0MZKCZAJEmSJEmaIdauXcuqVauoKlatWuUokB4yASJJkiRJ0gyxcuVKqgqATZs2OQqkh0yASJIkSZI0Q5x99tls3LgRgI0bN3LWWWf1OaLBYQJEkiRJkqQZYv/992fevHkAzJs3jwMOOKDPEQ0OEyCSJEmSJM0QS5cuJQkAc+bMYenSpX2OaHCYAJEkSZIkaYZYtGgRS5YsIQlLlixhxx137HdIA2OrfgcgSZIkSZLutnTpUtasWePojx5zBIgkSRp4SfZMcmaSW5Jcn2RFku3G2PZ+SY5PcnWSDUnWJDm6q04N87hlau5GkjToFi1axIoVKxz90WOOAJEkSQMtyQLgXOBK4CBgZ+A4YCfgxaO03R44HyjgSOCXwMOABw9TfQXw2Y7Xd04ydEmS1EMmQCRJ0qA7BFgI7FNVawGS3AF8JskxVXXRCG3fAiwAHlVVQyM6zttM3auq6tu9CVmSJPWaU2AkSdKgOxA4Zyj50foisAFYMkrb1wKf7Eh+SJKkWcoEiCRJGnR7Axd3FlTVBuAKYK/NNUqyG7ALsDbJl5PcnmR9kk8nWThMk7ck2djW+WKS3Tdz3gVJdut8AIsneG+SJGmMnAIjSZIG3UJg/TDl64AdRmi3S/v8T8CXgWcDuwLvp1lH5JkddT8NnA5cR5NweQfwzSSPqapfdZ33cOCocd2BJEmaNBMgkiRJwxsaKXs58PKqKoAkNwFfSPKEqroQoKo69ym8IMl5wE+A1wPv7Drv8cDJXWWLgQt6GbwkSbonEyCSJGnQraNZyLTbQuDSUdpBs35IdZSf0z4/CrhwuIZVdVmSHwGPG+bYerpGpCQZIQxJktQLrgEiSZIG3SU001LukmQ+sDsjJ0CuoFkodXO2mXxokiRpupgAkSRJg+4MYL8kO3aUPR+Y3x4bVlX9DjgL+PPcc4jG/u3z9zfXNskjgH3YzAgRSZI0/UyASJKkQXcCzZST05I8I8krgBXAqVV11+4wSU5MckdX23fTjB75XNv2r4CPAV+tqu+27Y5I8rEkByd5epL/B3wdWAt8dMrvTpIkjYlrgEiSpIFWVeuT7AssB74E3AacAhzZVXVu++hs+/0kz6TZ+eU04Ddt2zd3VPsZ8ALgRcDv0SQ+zgLeUVW/7vkNSZKkCTEBIkmSBl5VXcY9t60drs4yYNkw5ecCTxqh3X8B/zW5CCVJ0lRzCowkSZIkSRp4JkAkSZIkSdLAMwEiSZIkSZIGngkQSZIkSZI08EyASJIkSZKkgWcCRJIkSZIkDTwTIJIkSZIkaeCZAJEkSZIkSQPPBIgkSZIkSRp4JkAkSZIkSdLAMwEiSZIkSaNYu3Ythx56KDfccEO/Q5E0QSZAJEmSJGkUK1eu5Mc//jErV67sdyiSJsgEiCRJkiSNYO3ataxatYqqYtWqVY4CkWYpEyCSJEmSNIKVK1dSVQBs2rTJUSDSLGUCRJIkSZJGcPbZZ7Nx40YANm7cyFlnndXniCRNhAkQSZIkSRrB/vvvz7x58wCYN28eBxxwQJ8jkjQRJkAkSZIkaQRLly4lCQBz5sxh6dKlfY5I0kSYAJEkSZKkESxatIglS5aQhCVLlrDjjjv2OyRJE7BVvwOQJEmSpJlu6dKlrFmzxtEf0ixmAkSSJEmSRrFo0SJWrFjR7zAkTYJTYCRJkiRJ0sAzASJJkiRJkgaeCRBJkiRJkjTwTIBIkiRJkqSBZwJEkiRJkiQNvL4mQJLsmeTMJLckuT7JiiTbjbHtK5NcmuT2JBclOXiq45UkSZIkSbNT3xIgSRYA5wL3BQ4C3gS8BDhpDG0PAlYC/wEsAb4GfC7JkqmKV5IkSZIkzV5b9fHahwALgX2qai1AkjuAzyQ5pqouGqHtMcAXquqt7etzk+wNvBtYNZVBS5IkSZKk2aefU2AOBM4ZSn60vghsoBnVMawkDwX2Ak7pOvRZ4AlJdup1oJIkSZIkaXbr5wiQvema7lJVG5JcQZPgGKkdwMVd5UMjRh4BXN95oJ1us6Cr/q4AV1999ZgD7rbhputHr6SBsWbNmr5e/5p1v+vr9TW9NvX5/Xb7Dbf39fqaXpP5fOv47+jcXsSyBZsLk+uXSJK0pRutX9LPBMhCYP0w5euAHUZpxzBt17XPw7U9HDhquJM99alPHeFS0t0e+okj+h2CtiTHP7TfEWgL8tB39+T99gDgil6caAv1ALBfIklSjwzbL+lnAmQ6HQ+c3FW2NfAwYDVw5zTHM5stBi4Angr4NZWmmu83TSffbxMzl6aTcWG/A5nlLqR5712L/ZLx8O9W08n3m6aT77eJGbFf0s8EyDruPS0FmhEel47SjrbtdV3tAG7sblBV6xl+tMllI4eobkmGfry6qtb0MRRtAXy/aTr5fpsUR35MUlVtAP6733HMNv7dajr5ftN08v02KZvtl/RzEdRLuHs9DwCSzAd2Z+QEyCXt895d5Y9sn3/Wk+gkSZIkSdLA6GcC5AxgvyQ7dpQ9H5jfHhtWVf2cJkFycNehlwAXVpUrk0qSJEmSpHvo5xSYE4BDgdOSHAPsDBwHnFpVd+3wkuREYGlVdcb6TuDUdseYs4HnAgcAz5qu4CVJkiRJ0uzRtwRIVa1Psi+wHPgScBtwCnBkV9W5dG1hU1VfSLId8DbgCJo5Pi+tqlVTHrjWA+9m+DVVpF5bj+83TZ/1+H6TZpv1+Her6bMe32+aPuvx/dZzqap+xyBJkiRJkjSl+rkGiCRJkiRJ0rQwASJJkiRJkgaeCZAZKMnJSX46DddZlqSSLJrqa7XX68l9Jdk6ybFJvpHkt9N5D7q3JPskeVe7Lk+vzjmt703NPr38nEyyf5IfJ7k9SbVl5yU5vRfnl2Y7+yWjnsd+yQxiv0T9YL9k9jABsmX7CvBkZt/COtsBrwNuBy7ocyyCfYCjaH4v0mz0aeBqmt3EntznWKQtmf0S9cI+2C/R7Ga/ZAr1cxtc9VlVXQ9c3+84xqvdQWiHqqoky4Bn9DsmSbNTkgXALsC/V9U3puD8c4G5VfW7rvIA86vq9l5fU5qt7JdI2tLZL5l6jgCZwZI8I8lP2uFP30/y5I5jL09yQZIbkqxvf35KV/sHJTklya/ac6xJ8vGO4/cazpdkfpL3JPm/JBuSXJ3k5I7jJyf5aZL9kvwoyW1JvpVkryT3TbIyyU1JrkzyN5u5rz9N8oMktyb5YZI/HabOK9tjtydZm+SMJLsOHS+3L5q0jt/liL+P9r32o/Z3cV2S45LMb48tAz7VVr2+fT+tGTrWvn5ckjPbYcFXJPmLNN6a5Jr2PfzRJFsPE+biJKd3tH39MPfx5CRnJflNkpuTfCfJ/u2x3doYlib5RPu3cmOSDyeZ13GO+yc5sX3f35bk8iQfTLJt17Uqyd8neWeSa9tznZJk4SR+FVu0yX6mTObzpH3/rmurndj+fk8eKc5hyivJER2vz2vfsy9LcimwAXhSx30ekOQHbfnBSX6d5B+HOe+/JPnZ2P4VpekR+yX2S6ZQ7JfYL5kBJvuZMpnPk9gvmRYmQGauBwAfBz4IHAxsBL6aZOf2+G7AvwEvAl4MXAacm+QxHef4NM0wwMNovo14xxiu+0Xg74CTgGcBfw9s31VnF+B44P3AS4D7A6cCn6EZrnUQcBbw0SSPHabtPwMfBl4A/A74zyT3HaqQ5O+BlcAP2jqvAVYDO40hfo3PiL+PJIcBJwNfB/4CeBewDFjetv8K8J7252fSDNN7ftc1/g04sy1fDXwe+BDwWJrf7T8CfwW8YZj4PgucDzwX+AbwkSR/MnSw/fk8YBvgte09nAY8pOs8/wjMo/lb+lB7rWM6ju9IM+T679r7eB/wQpq/oW5vAH4feDXwVuDZwD8NU09jN9HPlMl+nnyF5vcNzfv4ydzzfTFRjwPeDrwbWAL8vC1/IM3n+kdoPpO/08a2NM03MkMxb0/zuX5iD2KResV+if2S6WC/pGG/pL/slwxyv6SqfMywB80HewH7dpQtBG4B3jdM/Tk005m+DyzvKL8FOHSE6yxrr7Oofb1/+/olo8S2CXj0MOf5aEfZPODGzng303aftu3z2tf3A34LnDDGf6t73IOPcb/PNvv7AO4L/Ab4QFe7FwN3ALuN9DvoKH99R9lubdlPgTkd5acD3xqm7Ru63lPXA8d3lH0TuIhmKN9w9zh0vW90lb+3fZ8t3Ey7rWg6N5uAHTvKC7iwq+7xwPp+/z5n62Oinym9+jwBFrVtlnWVnwec3hXnT4dpX8ARXe1+B+w6zH0W8Mdd5Xu25Qd2lL2K5n8u79/v348PH1X2S0b7HBnpHnyM+31mv2T4dvZL+vs+HPUzpVefJ9gvmfKHI0Bmrpuq6utDL6pqHU22+0kASfZO8qUk1wF30rwpHws8vOMcPwCOSPL/kuw5hmvuB9wKnDJKvWur6icdry9rn8/qiHcjsAZ48ChtL26fF7fPT6ZZtGowMowz30i/jyfTdDZOTbLV0AM4B5hL834bi873xRqaD+GvVdWmjjqXce/3CsBXO9pupMmQLwZIs7r7HwErq+rOUWL4j67X/07zPnt0e64kOTzJxUluo/l7+k8gNP8hGPZ+WhcD90tyn1Fi0OZN9DNlpn6e/Liqrhym/Iaq+p/OgqpaTdM5eU1H8WtoOjm/mroQpXGzX6LpYL8E+yUzgP2SAe6XmACZuYZbBOxXwAPaYVRnAQ8DjgD+FHgC8C2aIXdDDga+RjN06rIkq5O8eIRr7kjzh1ujxLau6/XQIjrrhynfpqvsHm3r7gV4hurt2D7/cpQY1Bsj/T6GhvZ+n+Y/vEOPX7fl3cM5x3SN9hzru8qGe68M17az3kKaz7CxvFd+3fV66AP8Ae3z4cBxwH/RfMv0JJpsN8PEtbn3/3Dxa2wm+pkyUz9PNtdB2Fz5vwDPSbJTkr2AP8H/2dLMY79E08F+SeNw7Jf0k/2SAe6XuAvMzDXcvNL7A9fSZA8XA8+pqh8NHWw7IDcMva6qa4HXJHktTVb8zcBnkvy4qi7m3m6g6chkDJ2NqTIU/wNp5tipf25sn18AXDXM8V9MYyzDWU8z1PCBY6i7c9fr+7fP17bPLwS+XFVvHqqQZLhvfjS79PLz5HbgHgviJdlhM3U39/m5ufIvASuAV9LMH/4lsGoCMUpTyX6J/ZJ+s1+i2c5+yQzgCJCZ635J9h16kWY1531pFqcZWgF6Q8fxx9AsgHQv1fg+8Baa3/lem7nm12iGZb1o0tFP3Ldohru+arSKmnL/QzNP8cFV9b1hHkNZ475801BVv6V5v7yyc6GmzeheAO0gmvfZ0DDFben4e2q9bNJBqt96+XnyC+DBnQuZ0aOtLqtqA82iY6+l6WycPIbh09J0s1+ifrNfotnOfskM4AiQmetGmu2P3kUznOptNFm649vjtwAfS/I+mqzx0XRkEpPcj2Y46r8CP6OZG/k3wM00nZV7qaqvJTkDOCnJ7m29HYCDqurgHt/fsKrqpiTvBj7Q/sfjP2k6R08HPldV3wNIsoRmFfjHt02fk+RmYM1QHU1O+7t4B83vYjFwLk2n4qE0K4y/vqquBi5pmxya5EvArV3zH6fSW2jmoJ+T5J9p/lYeC6ytqpM66j0syado5pE/nmYXgQ+3c9gBzgYOb1eXv5Tmm5d9pucWNFXG+nkyRl+k+Zz9VJITaNY1eF0Pw/0XmtX+i2a3C2mmsV9iv6Sv7JfYL5nt7JfMDCZAZq5raT4MPwjsQbOi9DOHsttJDmqP/SdwBfC3NAvUDC14dDvwv8DrgV3b198HnlFV14xw3RcARwGH0Gwt9ivuvbjSlKqqY5NcT3NPS2k6R9/invMlP0ZzX0OG/jBX0qzUrB6oquOTXA28iea9dAfNgk+raOc5VtUP2w7xa2nmfv+CZpXz6Yjvv5M8jWarsJNpFt67iHtvrfh24M9otrrbBHy0LRtyNM28zKNo/kP0ZZot8M6esuA1Lcb4eTKW81ya5OU0n4unARfSrGdwaY/i/FmSi4FfV9UVvTin1GP2S+yX9J39Evsls539kv5L/6ZUStLUSrIbzV7nL6yqf+9zONJmte/VK4BXVNVn+xyOJGkK2C/RbDHI/RJHgEiS1CdJdqQZtvpO4BqarRAlSZKm3ZbQL3ERVEmS+uc5wDeB3YGXd2yZJ0mSNN0Gvl/iFBhJkiRJkjTwHAEiSZIkSZIGngkQSZIkSZI08EyASJIkSZKkgWcCRJIkSZIkDTwTIJIkSZIkaeCZAJEkSZIkSQPPBIgkSZIkSRp4JkAk9USSdyWpJI9IcnKS9UluSvKpJNt11FuW5GtJrkuyIcnqJG9NMqfrfOcluTTJo5Ocn+TWJP+X5OD2+FOSfDvJbUl+luQZw8T0gCSf7LjWJUn+Zur/NSRJUj/ZL5E0nK36HYCkgXMK8H/AW4HHAq8Ffg28uT3+euAS4AzgdmA/4L3A/YC3dJ3rfsBXgM8DXwD+GvhMkgDHAx8HPgccAXwhyYOr6iaAJDsD3wbmAh9tY9gP+GiSHavqPb2+cUmSNOPYL5F0l1RVv2OQNACSvAs4CvhUVb26o/xLwJ9W1aL29XZVdWtX238BXgrsWFUb2rLzgD8DXllV/9qWPQK4FCjgqVX1zbb8AOCrwOuq6pMd5/wL4NFVdX3HtT4BvAx4YFWt7/E/gyRJmgHsl0gajlNgJPXaJ7peXwDsmOT3AIY6GUnmJlmYZBFwPrA9sFdX29uAzwy9qKqfAeuBy4Y6Ga3vtM8Pa88d4CCab2kqyaKhB3AWsC3wpMneqCRJmvHsl0i6i1NgJPXaVV2v17XPC4HfJHkKzdDSJwFbd9W9X9fra6pqU1fZTcAvOguq6qamb8HCtmin9udXt4/h7DzCPUiSpMFgv0TSXUyASOq1OzdTniQPA74GXAb8LU2n5HaaObkf4N6j0jZ3rs1eo30eOs/ngJM2U/eizZRLkqTBYb9E0l1MgEiaTn8BzAeeU1VXDhUmeWiPr3M9cDOwVVV9rcfnliRJg8F+ibSFcQ0QSdNp6BuSoW9ESDIfeEMvL1JVdwL/DjwvyWO6jyfZqZfXkyRJs5L9EmkL4wgQSdPpq8DvgNOTnEDzrcsrgO75tL3wFuBpwLfaFdYvopl/uw/wfGCbKbimJEmaPeyXSFsYR4BImjZVdRnwPGAjcCzwRuB04MgpuNavaRY0+2R7zY/QzO/dBXhTr68nSZJmF/sl0pYnVdXvGCRJkiRJkqaUI0AkSZIkSdLAMwEiSZIkSZIGngkQSZIkSZI08EyASJIkSZKkgWcCRJIkSZIkDTwTIJIkSZIkaeCZAJEkSZIkSQPPBIgkSZIkSRp4JkAkSZIkSdLAMwEiSZIkSZIG3v8PPmY1IRrj3dQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c = pd.read_csv('benchmarks.csv')\n",
    "a=pd.pivot_table(c,index=['allele','size'],columns='name',values='auc')\n",
    "r=pd.pivot_table(c,index=['allele','size'],columns='name',values='pearson r')\n",
    "print (len(a))\n",
    "def highlight_max(s):\n",
    "    is_max = s == s.max()\n",
    "    return ['background-color: yellow' if v else '' for v in is_max]\n",
    "\n",
    "fig = plt.figure(constrained_layout=True,figsize=(15,6))\n",
    "gs = fig.add_gridspec(1, 2, hspace=1)\n",
    "ax = fig.add_subplot(gs[0])\n",
    "#x.plot(x='basicmhc1',y='netmhcpan',kind='scatter',s=50,c='orange',ax=ax)\n",
    "#x.plot(x='basicmhc1',y='mhcflurry',kind='scatter',s=50,c='green',ax=ax)\n",
    "#ax.plot((0,1), (0,1), ls=\"--\", lw=2, c=\".2\")\n",
    "#ax.set_xlim(.76,.94);ax.set_ylim(.76,.94)\n",
    "sns.barplot(data=c,y='pearson r',x='name',ax=ax)\n",
    "ax.set_title('pearson r')\n",
    "\n",
    "ax = fig.add_subplot(gs[1])\n",
    "sns.boxplot(data=c,y='auc',x='name',ax=ax)\n",
    "#ax = fig.add_subplot(gs[2:])\n",
    "#g=sns.barplot(data=c,y='score',x='allele',hue='name', ax=ax)\n",
    "#plt.legend(bbox_to_anchor=(1.1, 1.05),fontsize=16)\n",
    "#plt.setp(ax.get_xticklabels(), rotation=90)\n",
    "#plt.tight_layout()\n",
    "ax.set_title('AUC')\n",
    "plt.savefig('basicmhc1_benchmarks.jpg', dpi=150)\n"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}