{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Figure S4: Non-independent cost of infection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "use czrecursion, cztogrowthrate\n",
      "use cstepmarkov\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.append('../lib/')\n",
    "from cycler import cycler\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import palettable\n",
    "import plotting\n",
    "import projgrad\n",
    "import scipy.optimize\n",
    "import evolimmune\n",
    "import analysis\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "plt.style.use(['paper'])\n",
    "plt.rc('axes', prop_cycle=cycler('color', palettable.colorbrewer.qualitative.Dark2_4.mpl_colors))\n",
    "black = matplotlib.rcParams['text.color']\n",
    "eps = 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define growth rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Lambda(p00, p01, p10, p11, pienv1, pienv2, lambda_, mu1, mu2, nu):\n",
    "    return pienv1*pienv2*np.log(np.exp(-2*lambda_-nu)*p00+np.exp(-lambda_-mu2)*(p01+p10)+np.exp(-2*mu2)*p11) \\\n",
    "    +pienv2*(1-pienv1)*np.log(np.exp(-lambda_)*p00+np.exp(-mu2)*p01+np.exp(-lambda_-mu1)*p10+np.exp(-mu1-mu2)*p11)\\\n",
    "    +(1-pienv2)*pienv1*np.log(np.exp(-lambda_)*p00+np.exp(-lambda_-mu1)*p01+np.exp(-mu2)*p10+np.exp(-mu1-mu2)*p11)\\\n",
    "    +(1-pienv2)*(1-pienv1)*np.log(p00+np.exp(-mu1)*(p01+p10)+p11*np.exp(-2*mu1))\n",
    "def Lambda_ni(x, *args):\n",
    "    p00, p01, p10, p11 = x\n",
    "    return -Lambda(p00, p01, p10, p11, *args)\n",
    "def Lambda_i(x, *args):\n",
    "    pi1, pi2 = x\n",
    "    p00, p01, p10, p11 = (1-pi1)*(1-pi2), pi2*(1-pi1), pi1*(1-pi2), pi1*pi2\n",
    "    return -Lambda(p00, p01, p10, p11, *args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimize non-factorizing case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pienv1, pienv2, lambda_, mu1, mu2 = 0.4, 0.4, 2.0, 1.0, 1.0\n",
    "nus = np.linspace(0, 2, 20)\n",
    "ps = np.zeros((len(nus), 4))\n",
    "fopts = np.zeros(len(nus))\n",
    "for i, nu in enumerate(nus):\n",
    "    res = projgrad.minimize(Lambda_ni, 0.25*np.ones(4), args=(pienv1, pienv2, lambda_, mu1, mu2, nu),\n",
    "                            jac=False, method='fast', disp=False, reltol=1e-6, nboundupdate=200)\n",
    "    ps[i] = res.x\n",
    "    fopts[i] = -res.fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimize independent solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ps_ind = np.zeros((len(nus), 4))\n",
    "fopts_ind = np.zeros(len(nus))\n",
    "for i, nu in enumerate(nus):\n",
    "    res = scipy.optimize.minimize(Lambda_i, 0.5*np.ones(2), args=(pienv1, pienv2, lambda_, mu1, mu2, nu),\n",
    "                                  bounds = [(0, 1), (0, 1)],\n",
    "                                  method='L-BFGS-B')\n",
    "    pi1, pi2 = res.x\n",
    "    ps_ind[i] = [(1-pi1)*(1-pi2), pi2*(1-pi1), pi1*(1-pi2), pi1*pi2]\n",
    "    fopts_ind[i] = -res.fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_pcov(ps):\n",
    "    fig, axes = plt.subplots(figsize=(4, 4), nrows=2, sharex=True)\n",
    "    \n",
    "    E1 = ps[:, 2] + ps[:, 3]\n",
    "    E2 = ps[:, 1] + ps[:, 3]\n",
    "    \n",
    "    ax = axes[0]\n",
    "    ax.plot(nus, E1)\n",
    "    ax.set_xlim(min(nus), max(nus))\n",
    "    ax.set_ylabel('fraction protected')\n",
    "    \n",
    "    ax = axes[1]\n",
    "    corr = (ps[:, 3]-E1*E2)/((E1*(1-E1))**.5 * (E2*(1-E2))**.5)\n",
    "    ax.plot(nus, corr)\n",
    "    ax.set_ylabel('protection\\ncorrelation coefficient')\n",
    "    ax.set_xlabel(r'non-additivity of costs $\\nu$')\n",
    "    ax.set_ylim(-1.02, 1.02)\n",
    "    for ax in axes:\n",
    "        ax.locator_params(nbins=5)\n",
    "        ax.grid()\n",
    "        plotting.despine(ax)\n",
    "    fig.tight_layout()\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAGICAYAAAC0gRwYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt8U+dhP/7PsXzBIMnmlkCQTEJIwMZKEyAXyW0uTcDG\n+24t7jAkbZJCXFi7fSHfAfuuaSEBs3272mzxtl87wM7Cmq6xXFi2bLJlkoaksUXS5tLIl1yAgCTC\nPZZ0BL7K5/eHrIOFLVs3S7L9eb9eflnn/vgg9NE5z3meR5AkSQIREdEwUhJdACIiSl4MCSIiCooh\nQUREQTEkiIgoKIYEEREFxZAgIqKgGBJERBQUQ4KIiIJKTXQBAMBsNsPlckEQBGg0Guj1+mHXMxqN\n0Gq1sNvtAIDS0lIAgCiKqK2tRU5ODpxOpzyfiIiiJEXgu9/9rrRo0SKppKQkks0DuN1uafPmzfL0\n+vXrg667YsUKeZu7775bnr9582ZJFEVJkqSYlImIiHzCvt0kiiIsFgsEQUBbWxva29ujCqnm5mZk\nZ2cHzAu2z8OHDwMAmpqaUFxcDACw2+3weDxQKpUAgEOHDkVVHiIiuibskHj55ZchCALKysogSRJM\nJlNUBbDb7VCr1fK0SqWC0+kcdl2lUgmz2YyGhgZs27YNANDW1gaVSgWLxQKj0Qiz2RxVeYiI6Jqw\nQ6K+vh5qtRqbNm0CgDH5UBZFMeiywsJClJeXY/Xq1QAAl8sFu90OvV6P0tJSVFZWwuPxBGzT2dmJ\n1tZWdHZ2xrysREQTWVgV13a7HW1tbVi3bh2USiUMBgMsFgva29uRm5sbUQEGV0QDvoDQarUjbqNS\nqQD4bktptdqA9bVaLaxWa0Dl98mTJ1FSUoK/+Zu/wS233BJROYmIxptly5ZFvY+wQqKhoQGCIKCw\nsBCA71t9c3MzTCZTxCFhMBhQX18vTwuCIO/LbrfLAWA2m9HU1ITdu3cDANxuN1QqFXJzc1FdXS1v\n73K5oNPphj3WwoULgy6j4Xm9XlitVuh0OigUikQXZ9zh+Yscz114JElCR/dVODwdcFxx4tJVD6KP\niDBDwv9hvmHDBkiSBEEQIAgCzGYztm7dGlEBVCoViouLYTQaIYoiysrK5GXPPvsstm/fjtzcXBgM\nBgiCAIvFgqamJpSXl0Oj0QAAysrKUF1dDUEQsGnTJrkS+3oKhYJvtgjx3EWH5y9yPHc+kiThYqcH\njisdcIi+IHB4OuQfu6cDnX29Aduszy+I+rghh4T/VtOqVasCAmHnzp1R33JauXLlsPNfeOEF+bVK\npZLXu74dhV6vD9q2gohoPJAkCc7uqzglfonT4mXYxcAAOHPFiW5vX9DtUwQBN03LglY5HfOU06Gd\nlh103XCEHBL+x143bdokf4MHgO9973uwWCxobm6OOCSIiCaDfqkf56+KOCVexmn3ZZwWv8Qp92Wc\nFn0/rp6uoNsqhBRoldOhUU4fCILsa4GgnI6507KQlnLtisvr9cakzCGHRGlp6bAtmfV6fdRtJYiI\nJoq+fi/sno7AAHBf9gWD+OWIVwM3TlXjZtVM3KyeAa1yhhwIGuV03DhVhdSU+N92S4puOYiIxhtn\n91WccF3CCdcFnHBdwnHXBZx0XcIp8TJ6+4f/Fq8QUjBfNRM3q2fiZtUMzFfNxHzVDNysnoUc1XRk\npqbH+a8YHUOCiCgIb38/HFc6cNx5ESdcvp/jA78vdXmG3SZDkYpF2TcOBMFMzFf7g2Am5k3LTsjV\nQDQYEkQ06fV4+3DSfQmfdJzHJx3n5CA4JV4OenvohkwVFmTNwsKsG7AwazYWZM3GwqzZmDctG4qU\nidPBNkOCiCaNfqkfDo8Tn3Scw8cd5/Gx8xw+/vIcTrovDXuLKC1Fgduzb8AC9WwszJ4th8Gt6tnI\nyshMwF8QfwwJIpqQLnV68EnHObR3nMMnzvP4uOMcPu04jyt9PUPWzVCkIn/GTVg8fQ5un36jHAw5\nqunj7vZQrDEkiGhc6+v34rjrIloun0HL5S98odBxftg6gxRBwAL1LCyePgeLpt+IxdPnYPH0ObhZ\nNXNC3SKKJYYEEY0bPd4+fOo8D+vlL9By+Qysl79A25dn0eXtHbLunKnqgTCYg9zpN2JR9hwszL4B\nmalpCSj5+MWQIKKk1NXXi0+uXob1k3fQ2nEWLZe/wMcd54bUHaQKKVgyYy7yZ85D/oybkDdjLhZN\nvxHZGVMTVPKJhSFBRAnX1deLlstf4A+XHWi5fAYfXTqD466L8Er9wPFr62UoUvGVWRroZs5D/syb\noJs5D4uyb8QUXh2MGYYEEcWVJEmwezrw3gUb3r/o+2n78uyQK4SpqemYn67CvTm3445Z85A/cx5u\ny74hoOsJGnsMCSIaU1d7e/DhJTvev2jH+xdO4/2L9iGVyhmKVNxz4824c5YGupka6GbehJxp02H9\n6CPceeed7AU2gRgSRBQzkiThc/elgKuE9o5z6JekgPXmq2Zi6Wwtls7OwbIbcpA7Y+6QK4RYdVBH\n0WFIEFHEerx9+OjSGVjOncS750/hg0t2OLuvBqwzNTUdd87SYNkN87F0thZ3zc7BrMzhx3yh5MOQ\nIKKQdXv78IdLDljOnsCx85/j9xdODxnoZmHWbCydneP7uSEHi7JvZBuEcYwhQURBdfX14sNLdljO\nncSxc75QGNyXUaqQgqWzc6CfswD3zrkFd83WYjofPZ1QGBJEJOvq68X7F21yKLx/0RYQCmkpCiy/\nYT70cxbgvjm3YPkN8zEtLSOBJaaxxpAgmsT6+r14/6Idb33xGSxnT+KDizb0DHoUNS1FgXtuvFkO\nhWWz52NqWvKNeUBjhyFBNMmc8Thx9MynOHrmEzSdPQH3oCEz01MUuPfGW6CfuwD6G2/B0htyknIg\nHIofhgTRBNfV14t3zn/uCwbHp/jMdUFeJkDAnbO0eGDebSiYeyvump3Dvo0oAEOCaIKRJAkn3Zfw\nhuMTHD3zKY6d+zygA7zZmUo8OO92PDhvEb5200LMmDItgaWlZMeQIJoAxJ4uNJ09Id9Gcnic8rJU\nIQX6OQsGguF25M6YgxSBj6RSaBgSROOUw9OBhtOtMNva8Lvzp9An9cvLcpQz8KDGFwqGubdCySeQ\nKEIMCaJxQpIkfOI8j4bTrWg43YqWL7+Ql01RpOH+ubfJt5FuUc+EIAgJLC1NFAwJoiTWL/Xjg4t2\n1A8Ewynxsrxs5pRpWJmTh6KcJSiYeyu7y6YxMeYh0djYiJUrV471YYgmjB5vHyznTqLhdCsabW04\n3ynKy7TK6SiavwRFOUuw/Ib57O6CxlxUIWE0GlFbWxt0uSRJOHPmDEOCaBRXe3vwxplP0HC6Fa87\nPg5ou7B4+hwUzV+CVTlLkDdjLm8jUVxFFRJutxvbtm2DRqMBAJjNZhgMBqhUKgCAw+FAa2tr9KUk\nmoCu9HajwdYG0ykrjp75VO7+QoCA5TfMR1HOEhTOz8Mt6lkJLilNZlGFRFlZWcB0Tk4O8vLy5Gmt\nVhvN7okmnL5+L9764jgOn/gAZlur3INqqpCCB266DUXzl2BlTh5unKpOcEmJfGJaJ2G324fMa2tr\ng16vj+VhiMYVSZLw4SU7/vPzj/Bfn38kj8qmEFLwdc0ifGPBnXhEsxhZGZkJLinRUDENiby8PJSU\nlKCgoACALyDWrl076nZmsxkulwuCIECj0QQNFaPRCK1WK4dRaWnpkOXZ2dmsA6GkcMp9GYePv4+X\nP3kHX1ivDdd512wtShbchT++5Q4OvkNJL6YhodfrUVVVBbPZDAB47rnnRr3lJIoiTCYTqqqqAAAb\nNmwIGhLV1dVobGxEfn4+Hn744YCQEEURDQ0NWLduXYz+GqLwfdl1Ba9+/hEOn/gA7120yfNvVs3E\n6lvvxOoFd2FBFusYaPyI+SOwgiCgqKgIGo0GDodj1PWbm5uRnZ0dMK+9vR25ublD1j18+DAAoKmp\nCcXFxQHL6uvr5SsYonjq7OvBEVs7Dp/8AEcdn8otn2dkTMMf36zDEu9UlOq/jtRUNkui8Sem71qj\n0Yimpibk5ORg69atcLvdo7aTsNvtUKuvVdKpVCo4nc5h11UqlTCbzWhoaMCePXvk+RaLBcXFxXj5\n5ZdHLJ/X6+Xg6mHyny+et0CSJOHdC6dQd/x91Nta4entBuBr+Vyszcc3F3wF9990G1IkwGq1or+/\nn+cwTHzvRcfr9UKhUES9n5iGhCiKAbeb8vLyQrqaGG4/wRQWFsJgMGD16tU4cuSIvK5SOfq93ePH\nj/MNFyGr1ZroIiQFj7cHRztOo+HySdi73QCAFAB3Km/EA9NzcJ96HqYq0oBLXWi9dO2c8fxFjucu\ncsuWLYt6HzENCf8VweDGPlardcQricEV0YAvIEarx/C3w2hvb4fNZoPb7YbRaERzczPsdjvy8vLk\nthuDLVy4EDqdLqy/abLzer2wWq3Q6XQx+VYyXv3hkgMvffou/uvzj+Rut29WzcSjt92N1Qu+EvSR\nVZ6/yPHcRSdWX4hjGhKSJGHLli0QBAF2ux0mkwmbNm0acRuDwYD6+np5WhAEuT7CbrfLgWE2m9HU\n1ITdu3cD8DXkU6lUKCwslLe12Wy44447hg0IAFAoFHyzRWgynrurvT145fMP8dLH7+Cjy2cA+B5b\nLZ6fj8cX34uCubeG3OX2ZDx/scJzl1gxDYnS0lLk5+fDZDKho6MD5eXlAY3rhqNSqVBcXAyj0QhR\nFAMa6D377LPYvn07cnNzYTAYIAgCLBYLmpqaUF5eHhAGFosFzc3NaG9vD3olQRSKjzvO4Rcfv4PD\nJ96HOFDXMHdqFr696B6su/1uzGFDN5pEYv64RV5eHvLy8mC320PuYybY7agXXnhBfq1SqeT1hntE\nVq/Xy08/EYWrq68XptMt+MXHx/C7C6cB+LrHeGjeIjyx+F48pFmE1BR+m6XJJ6YhUVdXhzVr1gC4\nVtdQU1ODp556KpaHIYqZk65L+OWn78L42e/R0X0VADBrihLrbl+Ox26/BzmqGQkuIVFiRR0SoijK\nTzDZbDa0t7fLy5xOJ06fPh3tIYhiSpIkvHHmUxxo/S1++8Vxeb5hzgI8vvg+FObkIV3BNg1EQIyu\nJGw2G/bt2wePxxPQ66tarR614pooXnr7vfivzz/Cz61v4uOOcwCArPQpWLNwGb6z6F4szL4hwSUk\nSj5Rh4T/CSODwYDm5uaAp42IksHV3h78+6fv4kDr2zhzxddQ82bVTPxZ/v341sKlyOSIbkRBxeya\n2h8WHo9HbtjmcDj4lBElzOUuD/613YIX2y1wDtQ3fGWWBj/QPYCinCUc1Y0oBDG98VpdXY36+noc\nOnQIgO/e7+DKbKJ4OC1exv6Wt1H72e/lhm8PzrsdP9A9AP2cBRzZjSgMMQ0JrVYrB4R/+tixY7E8\nBFFQLZfP4OfWt/DqqY/QL0lQCClYveBOfF/3APJmzE108YjGpZiGREtLy5A6CZvNFmRtouhJkoS3\nzx7Hz61v4a0vPgMAZKam4dHb7sbG/K9Bo5ye4BISjW8xDQmDwYCSkhK5f6SWlhY+3URjoq/fC9Op\nFvy85S1YB7rMmJExDevz9PjuYj2mT5mW4BISTQxjOuhQWVkZx7mmmOqX+vHKyT/g7z94DafEywCA\nHOUMbMz/GtbetgyZqekJLiHRxJLwQYeIQiFJEl53fIyfvGeW2zgsmTEXP9A9iD+6OZ9dZhCNkYQP\nOkQ0mmPnTuIn75nx+4E+lW7LugHbl67EqvlL+KQS0RhLykGHiADf00o/ec+Mo2c+BQDMm5aNrXc9\ngm/dupRtHIjiJOGDDhFd76TrIireP4JXT30EAJg5ZRo2f+Xr+M6ie5HBPpWI4irhgw4R+X1xxYV/\n+PA1GD97D16pH6q0DGzKvx9lS74KZVpGootHNCklfNAhoi+7ruCfP3oDBz8+hm5vHzIUqSjL/Sr+\nQvcAH2UlSrCkGHSIJidPbzf2t/wW+1t/C09vNxRCCr59+z14+s6HMXdaVqKLR0QAYlr7V1dXJ7/W\narWQJAk1NTWxPARNAF19vahufRuGup/i7z98DZ7ebnzjlq/gjdV/ib8rKGFAECURDjpEcfWavR07\nj70Km+dLAMBD8xbh/y5bifyZ8xJcMiIaDgcdoriwiV/i2XdexRG770vE0tk5+NHyVbh3zi0JLhkR\njYSDDtGY6urrxb+0vIV/+ugNdHv7MHPKNPx4eTG+tfAupAhs60CU7GI+6FBlZaVcN1FUVIRdu3bF\n6hA0zrzh+AQ7jv0XTomXkSIIeHKxHtuXrkB2xtREF42IQhTTp5tqampQUFAg32JqaWnB3r17sXXr\n1lgehpLcGY8Tz737KupP+2493jVbi7+975vQzWK9A9F4E9OQ0Gg00Ov18rRer4fdbo/lISiJ9Xj7\nsL/1bVT94XV09vViesZU/HB5Edbdtpy3lojGqZiGhNvtHjJPFMVYHoKS1G+/+Aw/PvZfOOG6CAEC\nvn37PfjrZYVsDEc0zsW876ann35aHnTIarWiuLg4loegJPPFFRfK3/0fuZ+lO2bOw9/ov4m7ZnMc\nEaKJIKYhUVhYCK1WC5PJBLfbjY0bN7Jbjgmqt9+LmtYm/P2Hr+FqXw+y0jPx18sK8djt97CHVqIJ\nJKYhUVdXB71ej23btsVyt5Rkms+ewI+P/Sc+dV4AAKy9bTmeWV6EmVOUCS4ZEcVaTEPCarVi1apV\nAfMcDgc0Gk0sD0MJIvZ0Yde7/42XP/s9AN/IcH+r/yaW3TA/wSUjorES05AoKChAfX09NBoNsrOz\nAQD79u3D888/H8vDUAK8e/4U/rLp17B7OqBMy8D/XVqIJxbfx1tLRBNcTENix44dyM/Ph0qlkucd\nO3YsloegOOv29uHFsx/hlY8+hQQJ+jkL8A9fWwONcnqii0ZEcRDTkCgvLx/SLYfFYhl1O7PZDJfL\nBUEQhrS1GMxoNEKr1cptL0pLSwO2b21thcFgYNcgMdL25VlsfvNlfOw8j4yUVPz18kI8lVfANg9E\nk0jMn26yWCxoaGgAADz66KNBP/D9RFGEyWRCVVUVAGDDhg1Bt6murkZjYyPy8/Px8MMPo7S0FHa7\nHS6XSw6Me+65BwUFBVAqWYkaKW9/P/6l5S1UfnAEvf1e3DIlG/tXPoHcmTclumhEFGcx/UpoNBpR\nUVEBrVYLrVaLZ555Bo2NjSNu09zcLNdf+A3ubnyww4cPAwCamprk9hcOhyOg59msrCw4nc5o/oxJ\nzSZ+iTUN+/H/3muAV+rHn+seQMXCh3F79o2JLhoRJUBMryRsNpv8QQ4AZWVlqKysxMqVK4NuY7fb\noVar5WmVShX0Q16pVMJsNqOhoQF79uwB4Ov6w994z+12w+12B32ayuv1wuv1hv13TQaSJKH2+HvY\n/bv/wZW+HsxXzcDfF/wp7pqpgdVq5XmLkP+88fyFj+cuOl6vFwqFIur9xDQkcnJyhszzf4ADgMfj\nCek20Ehdefi7JV+9ejWOHDkCAPI+d+zYgYMHDwbd9vjx43zDDcPZ24V/dvwevxPPAgAKZyzA+rlf\nQeoXHbB+0QHA93gzRY7nL3I8d5FbtmxZ1PuIaUjY7XY0NjbKTzeJooiWlhb5SqG2tnbI47CDK6L9\n22i1I3fp4N9/e3s7cnNzAfjqKx599FEsXrw46HYLFy4MCC0CGmyt+KHlf/Bl91XMzlTip/oSfF2z\nSF7u9XphtVqh0+li8q1ksuH5ixzPXXRi9YU4piHR0NAAm802ZL5/3nCPwxoMBtTX18vTgiDIH/x2\nu10ODLPZjKamJuzevRuA79aSPywaGhqwZMkS3HfffbBYLNBqtcPeclIoFHyzDRB7uvDsO6/CePw9\nAEDx/Hz8xLAaM4J0yMdzFx2ev8jx3CVWTENi9+7dIz7NNNzjsCqVCsXFxTAajRBFEWVlZfKyZ599\nFtu3b0dubi4MBgMEQYDFYkFTUxPKy8uh0Whgt9vx9NNPQxAESJIEQRCCVnyTj+XcSfyf3xrh8Dih\nSstA+X3fwLduvQuCICS6aESUZGIaEqM97hpsebCK7RdeeEF+rVKp5PUG70er1eLjjz8Ot6iTUo+3\nD3/3nhn7W9+GBAmGOQvwD18rxTxl9ugbE9GkFNOQoOR17qobm37zEt67aEOGIhU/XFaEDXkGNowj\nohExJCaB350/hU1v/BIXOkXcmjUb+x/6DhZNZ7sHIhodQ2ICkyQJ//bxMTz7zqvok/pRmJOH579W\nClX6lEQXjYjGiTG/1zBai2saG119vdjW9Gv86Nh/witJ2H7XChz4+ncYEEQUlpheSZjNZuzbtw8e\njweA75usw+Hg00Zx9oXHie+98RL+cMkBdfoU/NP96/CwNnj7ESKiYGI+6FBVVZXcF5MkSdi/f38s\nD0GjaD57At8/+u+43HUFi7JvRPXDj+MW9axEF4uIxqmYhoROpxvSWnrdunWxPAQFIUkSatqaUP47\nE7xSP/7XzTrs/eqfYlpaRqKLRkTjWExDwu12Y+/evQFdX5hMJo5MN8Y6+3qwvekwXjn5IVIEAT9a\nvgp/ln8/G8cRUdRiGhL79+9HXl4eXC6XPK+trS2Wh6Dr2MQv8b3f/AKtX55FdsZU/OyBR3H/vNsS\nXSwimiDGvFsOhsTYeevMZ/jBm7+Cs/sqlsyYiwNffxw5qhmJLhYRTSAx75ajsrISRqMRgiCgqKgI\nu3btiuUhCL76h5+3vIWfvNeAfknC6gV34qcFJchMTU900YhogolpSNTU1KCgoACbNm0CALS0tGDv\n3r3YunVrLA8zqV3p7cbWt3+N/z5lhUJIwXP3/BGeyitg/QMRjYmYhoRGowm43aTX6wPGiqDofO6+\nhLLXf4FPnOcxc8o0/PzBx2CYe2uii0VEE1jMn2663kijzFHofnf+FL772otw9XThK7M0OPDQd3AT\ne28lojEW05BQq9V4+umn5UdgrVYriouLY3mISel1+8fY9MYv0eXtxbduvQt/ZyjBlNS0RBeLiCaB\nmIZEYWEhtFotTCYT3G43Nm7ciLy8vFgeYtI5dOID/OVv6+CV+vFn+ffjR8tXsf6BiOIm5r3A5uXl\nMRhipLr1bTz37n8DAH60fBW+r3sgwSUioskm6pBob2+Xx6QebnjS2tpatrgOkyRJ+On7jfinj96A\nQkjBTwtKsPa25YkuFhFNQlGHxObNm+VGdBUVFdDpdJAkSV7OxnTh8fb34xnLK/jlp+8iQ5GKnz/4\nGFbm8MqMiBIj6pA4cuSI/HrPnj1DbjUxJELX7e3D/37zZZhOt0CVloF/feRJ3DdnQaKLRUSTWEzr\nJAYHhN1uhyAIrJ8Ikae3G0+9/m9oOnsCszOVeGnFBiyZeVOii0VEk1xMR6arq6uTX2u1Wl/31TU1\nsTzEhHSp04PS+v1oOnsCOcoZ+I/i7zMgiCgpRH0lIYoiHA4HAMBmswWMQud0OnH69OloDzGhOTwd\neMxcg5PuS8idPgcvrdyAG6eqE10sIiIAMbrdZLPZ5GFLW1tb5flqtVrux4mG+qTjPB5rrMH5q27c\nc+PN+NeHn0RWRmaii0VEJIs6JFQqFQoLC2EwGNDc3IzCwsJYlGvCe+/CaTxx5EW4ejqxQpuLnz34\nGDLZipqIkkzM6iT8YeHxeOR5/ttQFOg3jk+wtqEarp5OrFm4FAe+/h0GBBElpZhWXFdXV+PJJ5+U\npyVJCqjMJuDwiQ+w4bWD6PL2YtOSr2HvV/8UqSmKRBeLiGhYMQ0JrVaLQ4cOBUzTNTVtTdj8Vi36\npH48s3wVdtzzR0gRYvpPQEQUUzH9hGppaRkyz2azxfIQ49Y/fPgann3nVaQIAioLvoUfsB8mIhoH\nYtqYzmAwoKSkRO4qvKWlhU83AfgX61vY+8FrSE9R4GcPPoai+UsSXSQiopDEfIzrqqoqmM1mAMBd\nd90Fg8EQy0OMOy998g72/N6EVCEF+x76NlawHyYiGkdi3lW4VqtFWVmZ/JRTRUUFdu3aNeI2ZrMZ\nLpcLgiAMGQJ1MKPRCK1WKw+JWlpaGtb28fbKyQ/xw+ZXIEDA8/eXMiCIaNyJaZ2ExWLBihUrcO+9\n92L16tW4++67kZOTM+I2oijCZDKhtLQUa9aswYEDB4KuW11dDb1ej1WrVqGysjLs7ePpNXs7nn7L\nCAkS/p/hm/jmgjsTXSQiorDFNCSam5tx5MgR7N69G0eOHEF7eztUKtWo22RnB47VPLhrj8EOHz4M\nAGhqapKHRQ1n+3hpOnsCm974Jfqkfvxo+Sp8Z9G9CS0PEVGkYl5xDQAulwsejwdKpXLIB/j17HY7\n1OprfRWpVCo4nc5h11UqlTCbzWhoaMCePXvC3t7r9cLr9Yb1N4Xrg4t2bHjtILq9ffgL3YPYmPfV\nMT/mWPKXfTz/DYnE8xc5nrvoeL1eKBTRt8GKaUjY7XZUVlbixRdfxBNPPAGdTgeXy4WVK1eGtR9R\nFIMu83cBsnr16oCxLELZ/vjx42P6hjvV6cKPTr6BK95e/NHMhXgEs/Dhhx+O2fHiyWq1JroI4xrP\nX+R47iK3bNmyqPcR05AoLS2VK5OrqqpgsViwatWqEbcZXBEN+D7gR2uE57+F1d7eHtb2CxculB/P\njbVT7svY01APj7cXf3rrXagwlEyIhnJerxdWqxU6nS4m30omG56/yPHcRSdWX4hjGhJ1dXXQ6/XQ\naDTQarUhtbg2GAyor6+XpwVBkMfMttvt8j7MZjOampqwe/duAIDb7YZKpRpx++spFIoxebN94XHi\n26+9gItdHqyavwSVE7CrjbE6d5MFz1/keO4SK6YhYbVah1w5OBwOaDSaoNuoVCoUFxfDaDRCFEWU\nlZXJy5599lls374dubm5MBgMEAQBFosFTU1NKC8vl/cbbPt4uNTpwaONNXB4nHjgptvwzw88OuEC\ngogmr5iGREFBAerr66HRaOQK63379uH5558fcbtgdRYvvPCC/FqlUsnrXd8OItw6j1hxdXfi2401\nOOG6iLtvmI8DX38cGYqYNz0hIkqYmH6i7dixA/n5+QGPvR47diyWh0gaV3t78ORrL6L1y7PIn3ET\nXnzku5hmRa/wAAAgAElEQVSalp7oYhERxVTUIWGxWJCdnY3c3FyUl5cPGXTIYrFEe4ik0+3tQ9lv\nfoHfXziNhVmz8cvCDRxRjogmpKgfv2lubpZfC4IwZPlo7STGm75+L/786K/w1hefQaPMxr8XlmHm\nFGWii0VENCaivpKQJAl2ux1OpxNNTU1DWljX1taOWicxXvRL/dj69q/RYGvFDZkq/KqwDDdNy0p0\nsYiIxkzUIbFp0ybs27cPbrcbLS0tkCQpYHlbW1u0h0gKkiRhx7FXcejEB8jOmIp/L3wKt6hnJbpY\nRERjKuqQUKlU2LZtGwBfIOTlBfZ0OlFC4qfvN+LgxxZMS03HSyvWY/H0OYkuEhHRmItpk+DrAyLY\nvPGmpq0J//TRG8hQpOJfH3kSd87msKxENDmM/34jxthvv/gMu9/9H6QIAvY99G0Y5t6a6CIREcUN\nQ2IEp8XL+P7RX8E70OX3I9rhu/sgIpqoGBJBXOntxobX/g3O7qsoufUubFzytUQXiYgo7hgSw5Ak\nCf/nt3X4xHkeX5mlwd8ZSoZtA0JENNExJIbxj3/4DUynWzBrihIHHvoOMlPTEl0kIqKEYEhcp9HW\nhooPjiAtRYH9X/8OblJOrBbjREThYEgM8pnzAja/VQsAKL/vT3DPjTcntkBERAnGkBjg6u7Ehtf/\nDZ7ebjy+6F58Z9G9iS4SEVHCMSQAePv78edv/gqfuy/hnhtvxq57/zjRRSIiSgoMCQA/ec+Mo2c+\nxdypWdj30LeRzoGDiIgAMCTwyskP8fOWN5GhSEXNw49jdqZq9I2IiCaJSR0SLZfPYNvbhwAAlQV/\nijtmBR+Lm4hoMpq0IXGp04OnXv8Fury92JR/P1bfemeii0RElHQmZUj09nux6Y1f4swVJx646TY8\ns6wo0UUiIkpKkzIknnvnVbxz/nPcrJqJ/+/BR6FImZSngYhoVJPu0/HfP30XBz8+hmmp6ah5+Alk\nZ0xNdJGIiJLWpHrWs+3Ls/jRqTcBAP94/1osmn5jgktERJTcJtWVxN/+vh69/V785Z2PoHD+kkQX\nh4go6U2qkOjo6URRzhI8fefXE10UIqJxYVKFRI5qOp6/vxQpwqT6s4mIIjapPi13Li+GMi0j0cUg\nIho3JlVIzJ3GsSGIiMIxqUKCiIjCkxSPwJrNZrhcLgiCAI1GA71eP+J6ra2tMBgMKCwsBABYLBaI\noghJkgBAnk9ERNFJeEiIogiTyYSqqioAwIYNG4YNCbvdDpfLhdLSUgDAPffcg4KCAiiVSrS2tqKs\nrAwAsHPnToYEEVGMJPx2U3NzM7KzA+sK2tvbh6zncDjQ2toqT2dlZcHpdAIA9u/fj8bGRgCAIAhj\nWFoioskl4VcSdrsdarVanlapVPKH/2B6vR46nQ4A4Ha74Xa7odH4uvauqqrC+vXrkZWVhddffz0+\nBScimgQSHhLDEUVx2PlKpRIAsGPHDhw8eFCeX19fj/Lychw4cABPPvkkDh06FLBdd3c3AOCzzz4b\noxJPXF6vF59//jkUCgUUCkWiizPu8PxFjucuOl6vF2lpaViwYAEyMzMj3k/CQ0Kr1cJut8vToihC\nq9UGXb+6uhqPPvooFi9eDMBXma3T6bBmzRqsWbMGTz31FNrb25Gbmytv43A4AAA//OEPx+ivICJK\nTocPH8aSJZF3Q5TwkDAYDKivr5enBUGQP+DtdntAYDQ0NGDJkiW47777YLFYoNVqIQgCsrKy5HWK\nioqgUgUOQfrVr34VFRUV0Gg0yMhgYzoimjwWLFgQ1faC5H9uNIEaGxvhdDohiiLy8vLkp5s2bNiA\n7du3Izc3F3a7HStWrIAgCJAkCYIgyBXc1dXVEAQBarUaWVlZWLlyZSL/HCKiCSMpQoKIiJJTwh+B\nJSKi5MWQICKioBgSREQUFEOCiIiCYkgQEVFQDAkiIgqKIUFEREExJIiIKCiGBBERBcWQICKioBgS\nREQUFEOCiIiCYkgQEVFQDAkiIgoq5JAoKSnB4sWLkZubi8WLF2Px4sXYsmVL0KFG46mtrQ3V1dWJ\nLgYR0YQT8sh0DocDWVlZeP311+F0OtHW1obNmzfD4/GgpqZmLMs4IqPRiKamJtxxxx0JKwMR0UQV\n9vClSqUSSqUSGo0GWq0WLS0tY1GukJWWlgLAiFc0nZ2dOHnyZNQDghMRTTYR10k0NzfDbrdj1apV\nsSzPmDh58iRKSkpw/PjxRBdl3Onv78dHH32E/v7+RBdlXOL5ixzPXXRidd7CupJwu92499574XK5\nIAgCioqKsG3btpgUJB68Xi+8Xm+iizGueL1e9Pb2oq+vDwqFItHFGXd4/iLHcxcdr9eL9PT0qPcT\nVkio1Wq88847AHx1FOvXr8fDDz8sz0t2x48fZ0hEyGq1JroI4xrPX+R47iK3bNmyqPcRdp2En0aj\nwdq1a7F3715YLBbo9fqoCxMNSZJGXWfhwoXQ6XRxKM3E4fV6YbVaodPp+G0uAjx/keO5i06svhBH\nHBJ2ux0mkwkAEvrBazab0dDQAJfLhaysLKxZsybougqFgm+2CPHcRYfnL3I8d4kVckhoNBq0t7cj\nNzdXnqfVarF7924olcoxKVwoCgsLUVhYmLDjExFNZCGHxOHDh8eyHERElITYLQcREQXFkCAioqAY\nEkREFBRDgoiIgmJIEBFRUAwJIiIKiiFBRERBMSSIiCgohgQREQXFkCAioqAYEkREFBRDgoiIgmJI\nEBFRUAwJIiIKiiFBRERBJSQk7HY7HA5HIg5NRERhiFtI1NXVya+1Wi0kSUJNTU28Dk9ERBGIeIzr\nUIiiKF8x2Gw2tLe3y8ucTidsNttYHp6IiKI0piEB+MJh37598Hg8aGtrk+erVCrodLqxPjwREUVh\nTG83qVQqFBYW4uDBg3jooYfQ0dEBm80Gm82G1tZWVFZWjuXhiYgoSmN+JQH4wiItLQ1VVVXIzs4G\nAEiShP3798fj8EREFKG4hAQA6HQ6aLXagHnr1q2L1+GJiCgCcQsJt9uNvXv3BtRDmEwmPP/88/Eq\nAhERhSluIbF//37k5eXB5XLJ8wZXZBMRUfKJW0js3r0ber0+YB5DgogoucWtMZ1er4fD4ZDbTTgc\nDuTl5cXr8EREFIG4hYTRaERFRQVqa2sB+OooGhsb43V4IiKKQNxCQhRFVFVVIT8/HwB4FUFENA7E\nLSTUajUAQBAEeZ7Vao3X4YmIKAJxq7iWJAlbtmyBIAiw2+0wmUzYtGlTvA5PREQRiFtIlJaWIj8/\nHyaTCR0dHSgvL+ctJyKiJBe3kAB89RAMBiKi8WNMQ6K9vR25ubkAAIvFMmR5bW0tW1wTESWxMQ2J\nzZs3y43oKioqoNPpIEmSvJyN6YiIktuYhsSRI0fk13v27Blyq4khQUSU3OL2CGxeXh48Ho88zRbX\nRETJL24hUV1djSeffFKeliQpYNxrIiJKPnELCa1Wi0OHDgVMExFRcotbSLS0tAyZZ7PZ4nV4IiKK\nQNzaSRgMBpSUlMiDDrW0tLDFNRFRkotbSOj1elRVVcFsNgMAysrKeMuJiCjJxbXFtVarRVlZmTzt\n8XigVCrjWQQiIgrDmIeEw+GARqNhi2sionFoTEOisrISgiBg69at2LlzJwwGQ0CL69bW1rE8PBER\nRWlMQ0Kn06GwsBAA8L3vfQ+lpaUBy9nimogoucXkEdhgw5D6x7MGAgcbGm45EREln7CvJMxmM/bt\n2yd3sSFJEhwOB9rb24es29HRgZUrV0Kr1cJms+Hll1+Wl0mShDNnzmDlypVRFJ+IiMZS2CFhtVpR\nVVWF7OxsAL4P+/379w+7bk5ODnbt2gWNRgOz2SzfevI7cOBABEUmIqJ4CTskdDrdkPYN69atG3bd\n7Oxs6PV6AEBWVtaQ7f7qr/4q3MMTEVEchR0Sbrcbe/fulVtOA4DJZBr2UVabzQaLxSLfbrr+ltS+\nffv4CCwRURILOyT279+PvLw8uFwueV6wp5TWrl2LiooKuFwutLW1Demr6dixY+EenoiI4ijskPCP\nNDdYsJBQqVTYvXs3AN/wpddvN1wDOyIiSh5hPwKr1+tRWVmJe+65B/feey+effbZkAYP0uv1cDgc\n8mOvDodjSGgQEVFyCTskampqUFBQgNdffx2vvfYaioqKsHfv3lG3MxqNqKioQG1tLQBf3Uaw9hVE\nRJQcwg4JjUYDvV4PlUoFlUoFvV4fUm+uoiiiqqoK+fn5AMChS4mIxoGwQ8Ltdg+ZJ4riqNup1WoA\ngS2vrVZruIcnIqI4CrviWq1W4+mnn5YfgbVarSguLh51O0mSsGXLFgiCALvdDpPJxEGHiIiSXNgh\nUVhYCK1WC5PJBLfbjY0bN4Z066i0tBT5+fkwmUxwOp3Ys2cPcnNzIyo0ERHFR0S9wObl5UVUp2Ay\nmWA0GgH4blvt2rUrksMTEVGchBQS7e3t8rf+SAcP8j8V5b/F1NLSgr1792Lr1q3hlpmIiOIkpJDY\nvHmz3IiuoqICOp0uYPCgUMaF8D8V5afX62G32yMoMhERxUtIIXHkyBH59Z49e4bcagolJCJ9KoqI\niBIn7DqJwQFht9shCEJI9RORPhVFRESJE3Y7ibq6Ovm1VquFJEmoqakZdbvCwkJs3LgRHR0dOH36\nNDZu3MgBh4iIklxIVxKiKMp9Ll3f5bfT6cTp06dDOlhOTg62bdsGgEOXEhGNByHfbrLZbPKwpa2t\nrfJ8tVodUqO46upq1NfX49ChQwB8jevq6uqwZs2aCIpNRETxEFJIqFQqFBYWwmAwoLm5ecgwpKHI\nycmRAwLw3arieBJERMktrDoJf1h4PB55Xqi3jYbrp+n6QYiIiCi5hP10U6S3jQwGA0pKSuSnm1pa\nWmLWd5PZbIbL5YIgCEPaYxARUeTCfrpJq9UOuW0UCr1ej6qqKmi1Wmi1Wjz//PMxebpJFEWYTCaU\nlpZizZo1OHDgQNT7JCIin7CvJFpaWobUSYR620ir1aKsrCzcQ46oubkZ2dnZAfMGdyNCRESRCzsk\nxvK2USTsdrs8VgXgqzdxOp3Druvq7kRfvxepKYp4FY+IaFwLOyT8t43MZjMAoKysLORbTvESrLuP\nx468gH7rYajTpiArIxPZGZmYnj4V2RlTfdPpmcjOmIrpGVORnZGJrIHp7PRMZGVkIm0ShovX6w34\nTeHh+Yscz110vF4vFIroP7Mi6ipcEAQUFRVBo9EkvFGcVqsN6ChQFMWgoTU9dQo6U1Lh7u2Cu7cL\ndk9HWMeampIKpSIdqtQMqBTpA6/ToVL4fpTy6wwoFelQp6ZjmiINCiHsqp+kw1EEo8PzFzmeu8gt\nW7Ys6n2EHRJGoxFNTU3IycnB1q1b4Xa70djYGFEldKTbDWYwGFBfXy9PC4IQtD7i4IrvQqfTobff\nC3dPF5zdV9HRfRWunk44uzvh7L4KZ8/A74Bp32t3bxeu9vfhQu/VsMqoTpviuyIZdLWSnXHtKiU7\nia9cvF4vrFYrdDpdTL6VTDY8f5HjuYtOrK7Awg4JURQDbjfl5eWFdDVhNpvlFtuA79FZh8MR0MVH\nJFQqFYqLi2E0GiGK4ogV4wqFQv6ZkpaOG6apg647nL6BcOnovioHSUf3lYHfg+ddlQPIHy7u3i7Y\nPKMfI+BvS8sYCBB/wAyEyaBA8f3OlJdnZWRiamp6wFjiseA/bxQZnr/I8dwlVkRjXAMI+BCyWq2j\nXhFYrVZUVVXJTyJJkoT9+/eHe/hhxaujwNQUBWZMmYYZU6aFtV1vvxeugfBwdXfC2TM4RIZeuXR0\nXYWz5yrcPV0Qe7vDvi2WlqKQwyNr4KpkcKBkpV+bn5WeCbV/nfRMZKamxTxgiGj8CjskJEnCli1b\nIAgC7HY7TCZTSE836XS6IXUF69atC/fw41JaigKzMpWYlakMa7u+gNtigbfDXAPT/ltlrp7A6Utd\nHlzqCvPSZaCs6vQp18IjfQr6r3RhftdpZE+ZGjBfnZ4JVfqUgNcZioiquYgoSYX9P7q0tBT5+fkw\nmUzo6OhAeXl5SONJuN1u7N27V350FvCNeT3asKeTWaRXLpIkobOv91qgBARLJ9w9A6HS0wl3dxdc\nA9Punk64ujtxuesKLnddCdjnb12hjSKYoUhF1kBgqAfCQ5026PXAfP9yZVqGbzptCpTpU6BKy+Aj\nykRJJKKvfXl5ecjLy5MHHQrF/v37kZeXB5fLJc8LZUQ7Cp8gCJialo6paem4aVpWWNtKkoQuby9c\nPV2+W2SdHnz4SRtmzpsLd1+3HCRibxdc3V1w93RC7B2Y39MFsacLFzpFXOiMfNTBqanpUKVlQDUQ\nJqo0/++MQdMZvmBJy4AyfQqUqRlQpmf4pgd+GDZE0Qs7JAb30+R//LSmpgZPPfXUiNv5x8gejCGR\nfARBQGZqOjJT0zFnqhpe9SyknXXizlvvDKnyUJIkXO3rkQPD3dMJd0/XwE9guHh6u+Hu6YKnt2ug\n/qULnt5uiD3duNrXg/NRBA0ATFGkQZWegWmpvnCZlpZ+LVjSMjBtUKBMS8vAtNR03++09EHzfMtZ\nV0OTVdwGHdLr9aisrITRaJTbWezatSvCYlOyEgRh4IM2AwjzKsavX+rHld4eiL3dEAeHyEBFvjgo\nUDy9XRB7unGltxviwLRvvu/nYmcvLiL8upkhfxcETEtLl8NEmZaBqf4wSfWFyLS0DExNTce0tHRM\nTfX9ZCrScM59Fl3nsqDKmOKbn5aBqQPrJ8NjzkQjidugQzU1NSgoKJDXbWlpwd69e7F169YIik0T\nWYqQIt9qijRogGt1M2Jv17UQ6Rm4WhkIkSv+n74eeHq7cbV34Heff3kPrvR1BwRPRE69Pezs9BTF\nwJVb2rVgGXgdMD/t2rLM69a7Nj8NUxTpmDLwOlORhgxFKq+AKCpxG3To+i689Xp9QEtpolgbXDcT\nC97+flzt68GVvh45XDy93ejs68WV3m55WWdfjy9cBtY5c/ECMlRTcbWvF1f7enB1IHiu9vXiam+3\n/ODAWBAgYEpqKjIVviCZMhAe/rCZokj1/U5NwxRFGqYoUq+9Tk0dmJc2MG/QMkUaMoZZnpaiYChN\nMGHVSfjDorKyEnV1dQAQ8m0jt9s9ZF6wPpaIkpEiZdAVToi8Xi8+/PBD3Hln8Dqd3n4vOgcCpHPg\n52pvLzq9vkC52nft9eD1rg563dXXi05vLzr7egde98ivr/b55iPCi6BwpAgCMhSpyBi4ipmiSB34\nnRYw3/8zJfXafP+66Sm+32lCCs52nIHjlMK3Xmoa0lMUA/tQIH3Q+lMUqUgf+EkVUhhUMRR2xXWk\nt43UajWefvpp+RFYq9WK4uLiCIpMNLGkpSiQlu5rnzIWJElCT7/XFybevoEgGggU70Co9PUMvO7z\n/fYHjrcPXf71BtYNmOefHnjd09+Hrr4+XyjFiv13Ya0uQEC6QiEHSFqKP1AUvumB3xmKVKQrFNeW\np/hCxj+dnuJfNvBb3laBtIDXw62fitSUFPn4aSkpSBsIsHSFAinjqD+3sEMi0ttGhYWF0Gq1MJlM\ncLvd2LhxY0jtK4goOoL87T4+DR0lSUKf1I/ugTDp9vah2+sLE//r7oHX15b7gqbb24ee/mvLzpw/\nB/X0bPT0e9HT70W3txc9Xu/APrzyuj2D9t3T7x1Ypysuf28kUgRhIDz8AaRA6qDptIFgShVS5Pmp\nKSkBv9MGbSMvE3yBlJqiwBRFKn5wx4NRlzXsd000t4387SuIaOISBAFpgu/DS5mWEfF+QrlVN+L2\n/f3o6e9Dz0CY+IKjD739XvQMCiP/dLd/vrcPfQPb9g4Ejv91r9cXVtf24R2YP+h1v+84ffJ0P3r7\n+3y/vV70Sb59+tcbSwkJiXBuGw0eIc5isQxZXltbyxbXRDQmFCkpyExJR2YS9xTTL/Wjt78ffQNX\nSn2DgqhvIGD6+r3olbzo6+9HrzzfO/B68Lx+eVlffz8kSYpJGcM+feHcNtq8ebPciK6iogI6nS6g\n4GxMR0STWYqQggxFCjIUqQiv853RJayr8Lq6Ouj1emzbtm3UdY8cOSK/3rNnz5AwYUgQESW3sKvY\nrVar3N23XyjjSQwOCLvdDofDwfoJIqIkF3ZIFBQUoL6+HhaLBe3t7Whvb0dlZeWo2/nbVQC+Pp8k\nSUJNTU24hyciojgK+3bTjh07kJ+fD5VKJc87duzYsOvGos8nIiJKnLBDory8fEi3HMM9ueQ3Up9P\nGzduDPfwREQURxE93WSxWNDQ0AAAePTRR4d0Ae43Up9PFosFOTk5ERabiIjiIew6CaPRiIqKCmi1\nWmi1WjzzzDNobGwccRt/WPh5PB7odDpUVFSEX2IiIoqbsK8kbDYbDh8+LE+XlZWhsrISK1euHHE7\ni8WCnTt3wu12Q61Ww+FwhPQYLRERJU7YITHcLaLB41Z7PB4olcoh6zQ3N+PIkSMwm83yVYXRaAz3\n8EREFEdhh4TdbkdjY6P8dJMoimhpaYFarQYQvKsNg8EAAHC5XHKQXN/egoiIkkvYIdHQ0ACbzTZk\nvn9esMdh7XY7Kisr8eKLL+KJJ56ATqeDy+Ua9TYVERElTtgh4e+LKZhgj8OWlpaitLQUAFBVVQWL\nxYJVq1aFe3giIoqjsENipIAIZTkA+ckojyf6AeqJiGjsjGknuiM1sgPYVTgRUbIb05DYsmULVq1a\nFbRfc/YCS0SU3MY0JKqqqka8/cSQICJKbmM6Gvf1AeFwOOQO/9hVOBFR8hvTkBjM351HbW0tAN9Y\n2aN150FERIkVt5AQRRFVVVXIz88HAF5FEBGNA3ELCX+LbEEQ5HlWqzVehyciogiMacX1YJIkYcuW\nLRAEAXa7HSaTCZs2bYrX4YmIKAJxC4nS0lLk5+fDZDKho6MD5eXlvOVERJTk4hYSdXV10Ov17B6c\niGgciVudhNVqHdLrq/9xWCIiSk5xC4mCggLU19fDYrGgvb0d7e3tqKysjNfhiYgoAnG73bRjxw7k\n5+fL41AAwbsVJyKi5BC3kCgvLw8Y5xoYvQNAIiJKrLjdbnK73UPqIELpVpyIiBKHFddERBQUK66J\niCgoVlwTEVFQrLgmIqKg4hYShYWFqKyshNFohCAIKCoqwq5du+J1eCIiikDc6iRqampQUFCA119/\nHa+99hqKioqwd+/eeB2eiIgiELcrCY1GE/DIq16vh91uj9fhiYgoAnFtJ3E9URTjdXgiIopA3K4k\n1Go1nn76aeh0OgC+dhPFxcXxOjwREUUgrhXXWq0WJpMJbrcbGzdu5HgSRERJLm4hAQA5OTnyeBJs\nbU1ElPziVidRXV2NJ598Up6WJAl1dXXxOjwREUUgbiGRk5ODQ4cOydNarTZehyYiogjFtYO/69ls\ntngdnoiIIhC3OgmDwYCSkhL56aaWlhZs2rQpXocnIqIIxC0k9Ho9qqqqYDabAQBlZWW85URElOTi\n+nSTVqtFWVlZPA9JRERRiFudBBERjT8MCSIiCoohQUREQTEkiIgoKIYEEREFxZAgIqKgGBJERBQU\nQ4KIiIJiSBARUVAMCSIiCoohQUREQTEkiIgoKIYEEREFNe5Doq2tDdXV1YkuBhHRhDSuQ8JoNGLf\nvn0QBCHRRSEimpDiOp5ErJWWlgIARFEMaX2v1wuv1zuWRZpw/OeL5y0yPH+R47mLjtfrhUKhiHo/\n4zokQtXd3Q0AePPNN3H8+PEEl2Z8+vzzzxNdhHGN5y9yPHeRy83NxYIFC5CZmRnxPiZFSDgcDgDA\nz372swSXhIgovg4fPowlS5ZEvH3ShoTRaITNZguob5AkCdOnT8dTTz0V1r6++tWvoqKiAhqNBhkZ\nGbEuKhFR0lqwYEFU2ydtSPjrG0IhSdKIy2fMmIE/+ZM/ibZIRESTTtKGRCjMZjMaGhrgcrmQlZWF\nNWvWJLpIREQTiiCN9jWciIgmrXF9JXE9s9kMl8sFQRCg0Wig1+ujWm+yCfW8VFZWYtOmTXA6nTh2\n7Biv4OBr1Nnc3IyysrKg6/B9N7xQzh3fc8H531etra0wGAwoLCwccb2w33/SBOF2u6XNmzfL0+vX\nr49qvckmnPOyfv16afHixdKGDRskURTjUbykVltbK23evFmqrq4Oug7fd8ML5dxJEt9zwdhsNqm2\ntlaevvvuu4c9P9G8/8Z1i+vBmpubkZ2dHTCvvb094vUmm3DOy6pVq9De3o6amhoolcp4FC+plZaW\noqCgYMR1+L4bXijnDuB7LhiHw4HW1lZ5OisrC06nc8h60bz/JsztJrvdDrVaLU+rVKphT1ao6002\n4ZwXm82GxsZGOJ1OZGVlBb28pWv4vosO33PD0+v10Ol0AAC32w232w2NRjNkvWjefxMmJIYTancd\noa432QQ7L1u3bpVfr1ixgv9hI8T3Xej4ngvOf2W1Y8cOHDx4MOTtQn3/TZjbTVqtFm63W54WRRFa\nrTbi9SabUM+L3W6H2WyWp1UqFW+bhIDvu8jxPTe66upqPProo1i8ePGwy6N5/02YkDAYDHC5XPK0\nIAjIzc0F4HuThbLeZBbq+RNFMeDN5vF4eP4GSNc9Tc73XehGOnd8z42soaEBS5YswX333QeLxSJ3\nQxSr99+Eaifhv2cpiiLy8vLkR7w2bNiA7du3yycl2HqTXajnr66uDpIkwW63o7i4eNL/hzWbzait\nrYXL5cK6devkxzP5vhtdqOeO77nh2e12rFixAoIgQJIkCIIgX2XF6v03oUKCiIhia8LcbiIiothj\nSBARUVAMCSIiCoohQRNWdXU1jEZj2MuI6BqGBE1YS5YsQW1tbdBlL7/8sjy9YsWKgGfvRVEc0tjo\n+nVGEsr+oiWKIiorK7F3715YLJaY7vv647Dh3+TFkKAJa7juCQYvGzzq4e7duwMeq2xubpafNw+2\nzkhC2V+09u3bhzvuuAPFxcVD+uWJpbEoO40fDAkiIOCZcVEUsW/fvhHXicX+ouV2u6FSqZCbmztm\n7VeRdFgAAATVSURBVAbGquw0fiiee+655xJdCJpYLBYLNmzYgMzMTDgcDrzyyisQRRELFy4E4GtA\n9cEHH8DhcODNN9/E0qVLR93mekajEW63G0ajEYIgyF0MWCwWNDY24tKlS3j//fdx/PhxrF27dsRl\nbW1t2LRpE3JycqDVanH06FG8+eabkCQJZ86cwZIlS4asYzab8dhjj8HtdsNgMMjlv/vuu3Hx4sUR\n9+dwOAK2NZvN+P73v4+7774bs2fPHvK3Dne+2tra8PLLL0MURfT29uLWW28N+m9hNBrh8Xjwxhtv\nYOnSpUH3abfb8etf/xrd3d0wm83o7u7G8ePHcfTo0YBzMdx6w3XxsHPnTrz00kv4xje+AcA3bsTB\ngwdhMBjCfEdRQoXfgznR6DZv3ixVVlZKkuTr83716tWSJElSa2urtGXLFnm95uZmqaKiYsRthlNS\nUiK/fuSRR4bdprW1VV5vpGWSJEk7d+6Umpub5en169dLbW1tAce8fp3KykrJaDTK+x+8bLT9jbTt\nYCOdrx07dgTdbri/uaSkRBJFMeg+q6urpdbWVnlbf3mvL/uBAweGXe96DQ0N8r+nv7xmszloeSk5\n8XYTjYmsrCy5C+Ps7Gx4PB4AgMlkkucDQH5+vvyUUbBthvPiiy/CaDTCbDbL/fqYzeaAsQlUKpX8\neqRlwNC+g4Zz/TqlpaX41a9+BcD3jX3wLabR9jfStoONdL5GYzabA7Y9dOgQlEpl0H0WFhbiu9/9\nLp566ilYLJagt7BCXa+pqUm+anC73WhtbcXKlStDKjslD4YExdXgjtr8BlcgD2fDhg1YsWIFVq5c\niZqaGtjtdjz55JMoLi5GYWEh1Gr1iIESrcEdpQ2m1WohCALa2tpG/RsGczgc8rZ2u33EbSM5X6MJ\nts/s7Gz85je/wdq1a2E2m1FXVzdkPbvdHtJ6gC/8/GG0Y8cO/OM//mNU5abEYEjQmJMkSf5mvW7d\nOlitVnlZS0sLioqKRtzmhRdewJEjR9DY2Ch/e9XpdHI/+v4nb2bNmoW2tjZ5Hw6HQ95HYWFh0GXD\nUavVcDqdo36Ir127Fj/+8Y9HrdQevL/B227ZsmXEbUM9X8MpLCxES0uLPG232+FwOILuc9++fZAk\nCStXrsTzzz8Pm802bNmDrTeYKIoQBAFKpRI7d+7E9u3bMW/evJDKTcmFFdcUc21tbaiurobT6UR+\nfj4OHjyIN998E7fffjvuu+8+pKSk4OjRozh+/DhaW1vxzDPPjLjN9ZWyOTk5qK+vR0pKCt577z0s\nX74c9fX1ePzxx5GVlYWjR4+iu7sbra2tOHr0KObPn4+lS5ciMzNz2GXd3d04cOAALly4IH8Az58/\nH//xH/+By5cv48EHH5TLN3gdwNfewmKxoLS0dMjfP9L+/H/Hm2++iccffzzouZw9e/aI5+vChQuY\nP3/+sBXeWVlZmD59Ol555RW5wlyv1wfdZ3t7uxwkra2tKCoqglqtHlL2YOsN9v777+PEiRPo6enB\nX/zFXwxZTuMHe4ElSqDGxkbep6ekxttNRHFWXV0tj7bGgKBkx5AgirOsrCwcO3aMt2BoXODtJiIi\nCopXEkREFBRDgoiIgmJIEBFRUAwJIiIKiiFBRERBMSSIiCgohgQREQXFkCAioqD+f43WRzyup04/\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8784122cd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plot_pcov(ps)\n",
    "plotting.label_axes(fig, xy=(-0.15, 0.97))\n",
    "fig.savefig('SIniinfection.pdf')\n",
    "fig.savefig('SIniinfection.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Optimal protection strategy against two equally frequent pathogens $\\pi_{\\rm env, 1} = \\pi_{\\rm env, 2} = 0.4$ as a function of the degree of non-additivity of the cost of infection $\\nu$.**\n",
    "**(A)** Fraction of population protected against a particular pathogen. **(B)** Pearson correlation coefficient between the protection states against the two pathogens. As costs are non-additive, the problem no longer factorizes and the optimal strategy no longer chooses protections against different pathogens independently. However, here the optimal strategy treats each pathogen almost indendently, as measured by the low correlation coefficient. With an increasing cost of co-infection, more protection is needed, in agreement with our intuition that co-infection leads to higher effective costs. Parameters: $c_{\\rm infection} = 2$, $c_{\\rm defense} = c_{\\rm constitutive} = 1$, optimization of the distribution over protection states respecting the probability simplex constraints using an accelerated projected gradient algorithm as described in [Mayer et.al. 2015]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:evolimmune]",
   "language": "python",
   "name": "conda-env-evolimmune-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "41px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}