Figure6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 6\n",
    "This notebook recreates Figure 6 in Rein & Tamayo 2016. The figure illustrates the scaling of second order variational equations as a function of the number of varied parameters.\n",
    "\n",
    "We start by import the REBOUND, numpy and matplotlib packages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feb 15 2016 12:13:46\n"
     ]
    }
   ],
   "source": [
    "import rebound\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "from timeit import Timer\n",
    "print(rebound.__build__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function initialises a simulation and adds a given number of variational parameters. It then runs the simulation for 100 time units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def evaluate(order=0,N=0, integrator=\"ias15\"):\n",
    "    sim = rebound.Simulation()\n",
    "    sim.integrator = integrator\n",
    "    sim.add(m=1.)\n",
    "    sim.add(primary=sim.particles[0],m=1e-3, a=1)\n",
    "    sim.add(primary=sim.particles[0],m=1e-3, a=1.8,f=1.4)\n",
    "    \n",
    "    vs = [\"a\",\"e\",\"i\",\"omega\",\"Omega\",\"f\",\"m\"]\n",
    "    var = []\n",
    "    if order>=1:\n",
    "        for i in range(N):\n",
    "            pi = 1\n",
    "            shifti = 0\n",
    "            if i>=7:\n",
    "                pi = 2\n",
    "                shifti = -7\n",
    "            var_d = sim.add_variation()\n",
    "            var_d.vary(pi,vs[i+shifti])\n",
    "            var.append(var_d)\n",
    "        if order >=2:\n",
    "            for i in range(N):\n",
    "                pi = 1\n",
    "                shifti = 0\n",
    "                if i>=7:\n",
    "                    pi = 2\n",
    "                    shifti = -7\n",
    "                for j in range(i,N):\n",
    "                    pj = 1\n",
    "                    shiftj = 0\n",
    "                    if j>=7:\n",
    "                        pj = 2\n",
    "                        shiftj = -7\n",
    "                    var_dd = sim.add_variation(order=2,first_order=var[i+shifti],first_order_2=var[j+shiftj])\n",
    "                    if pi==pj:\n",
    "                        var_dd.vary(pi,vs[i+shifti],vs[j+shiftj])\n",
    "    \n",
    "    sim.move_to_com()\n",
    "    sim.integrate(100.)\n",
    "             \n",
    "    return "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wrapper function to make `evaluate()` callable from a timer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def evaluateWithN(order,N,integrator=\"ias15\"):\n",
    "    def _e():\n",
    "        evaluate(order,N,integrator)\n",
    "        pass\n",
    "    return _e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now run the timing test for different numbers of parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Nmax = 15\n",
    "repeat = 1\n",
    "t = Timer(evaluateWithN(0,0))\n",
    "var_0 = t.timeit(repeat*2)/2.\n",
    "var_1 = np.zeros(Nmax)\n",
    "var_2 = np.zeros(Nmax)\n",
    "for i in range(Nmax):\n",
    "    t = Timer(evaluateWithN(1,i))\n",
    "    var_1[i] = t.timeit(repeat)\n",
    "    t = Timer(evaluateWithN(2,i))\n",
    "    var_2[i] = t.timeit(repeat)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following plot show the runtime of the simulation with variational equations compared to one without."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEOCAYAAAB8aOvdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FFXWwOHfSQhgICEJ+xICsqiIgAsuOEoQFWFEGVQU\nwRkQPxhHBUcdBRUIIuOGOIOiwogZFRERccENFcQVF9wQUUGJgCyOQCIJEJJ0zvdHdZokdCfdne50\nOjnv89RDV1XXrZMm6VN17617RVUxxhhjvImJdADGGGNqLksSxhhjfLIkYYwxxidLEsYYY3yyJGGM\nMcanepEOIBgiYl2yjDEmCKoqgbw/au8kVDVql6lTp0Y8Bos/8nFY/NG3RHPsqsFdW0dtkjDGGBN+\nliSMMcb4FLVJIiMjg1WrVkU6jKCkp6dHOoQqsfgjy+KPnGiNfdWqVWRkZAR1rARbTxVJIqLRGLcx\nxkSSiKABNlxHZe8mXxo1asT+/fsjHYYxtVJ8fDz79u2LdBimmtWqOwl3loxARMbUfvb3Ff2CuZOI\n2jYJY4wx4WdJwhhjjE+WJIwxxvgUtUkimrvAGmNMdbIusIe2W8NaEO666y6ysrKYN29eUMd3796d\nhx9+mDPPPDPEkZmaxP6+ol8wDdeWJExARo8eTWpqKnfccUekQzHVzP6+op/1bqqCrKwsRo4cSb9+\n/Rg5ciRZWVmRDqlauFyuSIdgjKnJIj0qYZAjGao3vrZXZtOmTdqpUycFPEunTp1006ZNQZXnS4cO\nHXTmzJnao0cPTUpK0ssuu0wPHjzo2T9v3jzt3LmzNm3aVC+88ELdvn2713IGDhyoc+bMKbOtZ8+e\n+sILL6iq6oQJEzQ1NVUTExP1pJNO0vfff9/zvoyMDL344ot15MiR2qRJE50/f75mZGToyJEjPe+5\n5JJLtFWrVpqUlKR9+/bV9evXe+KLi4vTBg0aaEJCgl5wwQWen2vFihWqqnrw4EGdMGGCtmnTRtu2\nbavXX3+9FhQUqKrqqlWrtF27dnr//fdrixYttE2bNpqZmek576uvvqrdunXThIQEz/tMzRHs35ep\nOdz/h4F93wZ6QE1YQp0kRowYUSZBlCwjRowIqjxfOnTooKeccoru3LlTs7Oz9ZhjjtG5c+eqquqK\nFSu0WbNm+tVXX2lBQYFed911euaZZ3ot58knn9TTTz/ds/7tt99qcnKy58v46aef1uzsbHW5XDpr\n1ixt1aqVJxllZGRo/fr19eWXX1ZV1QMHDmhGRoZeccUVnvIyMzN13759WlBQoH//+9+1V69enn2j\nRo3SyZMnH/ZzlSSJyZMn62mnnaa7du3SXbt2aZ8+fXTKlCmq6iSJevXqaUZGhhYVFelrr72m8fHx\nmpOTo6qqrVu31g8//FBVVXNycvTLL78M8pM24WBJIvpZkvCy3duXf6gXf3Xo0EEXLlzoWb/55pv1\n6quvVlXVMWPG6C233OLZl5eXp3Fxcbp58+bDysnNzdXGjRvrli1bVFX1tttu0zFjxvg8b3Jysq5d\nu1ZVnSTRt2/fMvvLJ4nSsrOzVUR07969qlp5kujUqZO+8cYbnn3Lly/Xjh07qqqTJOLj49Xlcnn2\nt2jRQj/55BNVVU1LS9N58+Z5zmVqFksS0S+YJFHr2yT8+RBGjBjh9dgRI0b4m7T81rJlS8/r+Ph4\n8vLyANi+fTtpaWmefY0aNaJp06Zs27btsDIaN27MoEGDWLRoEQDPPPNMmZ9h5syZdOvWjeTkZJKT\nk9m7dy+7du3y7E9NTfUZX3FxMRMnTqRz584kJSXRsWNHRKTM8RXZvn077du396ynpaWxfft2z3rT\npk2JiTn0a1f6M3j++ed59dVXSUtLo1+/fnz88cd+ndMYEz61Pkn4Y/r06XTq1KnMtk6dOjF9+vRq\ni6FNmzZs3rzZs75v3z52795N27Ztvb5/+PDhLFy4kI8//piDBw/Sr18/AD744APuu+8+lixZQnZ2\nNtnZ2SQmJpZJZiK+Ozc8/fTTLFu2jJUrV5KTk8PPP/9cJhlWdKy3n2Pz5s20adOm8g8AOPHEE3nx\nxRf57bffuPDCCxk2bJhfxxljwseSBNCxY0feeustRowYQb9+/RgxYgRvvfUWHTt2rLYYhg8fTmZm\nJmvXruXgwYPceuutnHrqqWWuyksbNGgQmzdvZsqUKVx66aWe7bm5ucTFxdG0aVMKCgq44447yM3N\n9TuOvLw8GjRoQHJyMvv27WPSpEllEkPLli3ZtGlThT/HnXfeya5du9i1axfTp0/niiuuqPS8hYWF\nLFy4kL179xIbG0tCQgKxsbF+x22MCY+oTRKhfuK6Y8eOLFiwgJUrV7JgwYKwJIiKrsL79+/P9OnT\nGTp0KG3btiUrK8tTneRN/fr1GTp0KCtWrODyyy/3bB8wYAADBgyga9eudOzYkfj4+Aqrl8r785//\nTPv27Wnbti3du3enT58+ZfaPGTOGb7/9lpSUFIYOHXrYz3X77bdz0kkn0aNHD3r27MlJJ53Ebbfd\n5vN8pY996qmn6NixI0lJScybN4+FCxf6Hbcxxjd74vrQ9oDbCIwx/rG/r+hnD9MZY4wJKUsSxhhj\nfLIkYYwxxidLEsYYY3yyJGGMMcYnSxLGGGN8siRhjDHGJ0sSxhhjfLIkYQDo168fjz/+eKTD8Oqu\nu+5i7NixQR/fvXt33nvvvRBGZEzdYUmiGs2ZM4fevXvTsGFDrrzySr+P69ixIytXrgxjZDXbpEmT\n/J5/e/To0UyZMqXMtnXr1tn828YEqV4gbxaRGKCxqu4NUzwRcffYseRv2HDY9oZduzLRzy8nf7Rt\n25bJkyezfPlyDhw4ELJyA+FyuUI6cF6oy6vu8o0xlfBjroSFQCLQCFgP/AL8I9CJK0K5EOKZ6ab2\n7evMv1RumVpucp5Quf3223X06NFltu3atUvPP/98TUpK0pSUFM+sdFdccYXGxMRofHy8JiQk6H33\n3ee1zIqmPhURnTNnjnbp0kWPPPJIVVV988039eijj9akpCS99tprtW/fvjp//nzPMfPnz9djjjlG\nU1JS9Lzzzisz+ZG38kqz6VVrp2D/vkzNQThmpgO+cv87ArgfiAPWBnqiUC6ATp06Vd955x1vH0DA\nakKSmDRpkl599dXqcrm0qKhIP/jgA8++Dh066MqVK32WV9nUpyKi5557rubk5Gh+fr7u2rVLExIS\ndOnSpVpUVKQPPPCA1qtXz5MkXnzxRe3SpYv+8MMP6nK5dMaMGdqnTx+f5ZVn06vWTpYkotc777yj\nU6dODVuS+NadGJ4D+rq3fR3oiUK5BHQnMXXqoS/+qVO97p+aluY9SVRynNftfvCWJKZMmaJDhgzR\nH3/88bD3l/6S86ayqU9FRFetWuXZ/+STT+ppp51Wpox27dp5ksTAgQP18ccf9+xzuVwaHx/vmS61\nfHnl2fSqtZMliegXTJLwp+F6LvCzu7rpPRFJA6KnTSIj49DXvrfx1DMyoEMH78f27VvxcUGOz+7N\nzTffTKdOnTj33HPp3Lkz99xzj9/H+jP1abt27cq8v/wcE6XXN2/ezIQJE0hJSSElJYWmTZsiIj7L\nK8+mVzWm9qg0SajqbFVtq6qD3MloM9CvGmKrUxo1asTMmTP56aefePnll5k1axbvvPMOEPiUoSVT\nn5b+Ii9dRuvWrdmyZUuZMrZu3ep5nZqayty5c9mzZw979uwhOzubvLw8Tj31VK/leWPTqxpTO1Sa\nJESkgYhcLiK3isgUEZkC3FoNsVWbhl27ktG372FLw65dQ3oel8tFfn4+LpeLoqIiDh48iMvlAuDV\nV1/lp59+AiAhIYF69ep5evX4M2Wot6lPfV2R//GPf2T9+vW8+OKLuFwu/v3vf7Nz507P/r/+9a/8\n85//ZP369QD8/vvvLFmyJKCf1aZXNaaWqKw+CngDeBa4GbixZAm0XiuUCyHu3VRdMjIyVEQ0JibG\ns0ybNk1VVR944AHt0KGDNm7cWFNTU3XGjBme41566SVt3769Jicn++xRM3fuXO3UqZM2bdpUBw8e\nrNu2bfPsi4mJ0Z9++qnM+5cvX65du3bVpKQkve666zQ9Pb1M76YFCxbocccdp02aNNH27duXaVPw\nVp43Y8aM0ZiYGF2zZo1nm8vl0iuvvFITExO1TZs2et9992nHjh09bQbe2h9Kb8vLy9MLL7xQExIS\ntEOHDvrUU0+ViWfjxo3aq1cvTU5O1j/96U+qqmXKz8/P1wkTJmjr1q21TZs2ev3113sazVetWqWp\nqallzl1ybEFBgZ533nmakpKiTZo00ZNPPlk/+uijSj+D2qSm/32ZyhFEm0Sl05eKyDpV7R7GPBUw\nm77UmOpnf1/RL1zTl34kIscFGZMxxpgaINgE788T138ARolIFnAQEOd82iOoMxpjjKl2zy97Pqjj\n/EkSA4Mq2RhjTMTNfXwus+fPprBFYVDH+9MFdjOQBAx2L0nubcYYY2oQVSX34KHego+ueZSf2v1E\nxs0Z5BfkB1WmP11gJwBPAy3cywIRuS6osxljjAmZg0UH+TnnZ8/64m8X8/flf/esj+wxkhn9ZyAi\n5OTlBHUOfxquxwCnqOoUVZ0CnAr8X1BnM8YY45WqMnHaxAobmPcV7GP11tWe9S93fsnkdyZ71i85\n9hIeu+Axz3rj+o2Ji41jY9ZGMm/KDCouf5KEAK5S6y73NmOMMSHy/LLneXjlwyx9Zaln24HCAyxa\nt8iz/vvB35nz2RzP+qntTuWpPz3lWY8R71/pkyZM4qLBFwUVlz9JIhP4REQyRCQD+BiYH9TZjDHG\nlDE3cy7Hnn4st2beSm6/XMY9NI5up3djbuZc6sXUY8WmFRRrMQBtEtqwYOiCao3Pn4brWcBoYI97\nGa2q/wp3YCZyqmsmvM2bNxMTE0NxcXFYyn/iiSc444wzwlJ2VW3duvWwMasCUdUpXU3NkdUuixsn\n3Og0LAsUuYqYdOMkxo4aS1xsHP+54D8+7xCqg88zi0ii+98UnFFgF7iXze5tJkC1YfrSmTNnctxx\nx5GYmEinTp2YOXNmlcqrbEC+qgp3+cFKTU1l7969fsX37rvvHjYOVyBTupqa5ea3bmbtr2s96ye2\nOZEYiSEnL4dun3ejuLCY+Lj4GvO7W1F6Wuj+93NgTamlZL3W8afhqCpKpi8dM2ZMWMqvLk899RQ5\nOTm8/vrrPPTQQyxevDjSIYX8biRcdzeAZ1BHf6lqjfnCMP4p/R0y6e1JLPthmWd96DFDaZd4aITm\nS469hB3bd5B5UybrXlpH5k2ZbMzaWK3xVijQwZ5qwkKYBvh77qXnNOHMBF3y8pIqlVOZUE9fmp2d\nreeff742b95cU1JS9Pzzz9dffvnFsz89PV0nT56sp59+uiYkJOiAAQN09+7dnv1PPvmkpqWlabNm\nzXTGjBmVTnJU2vjx43X8+PGedRHRRx99VLt06aLJycl6zTXXePa5XC698cYbtVmzZtqpUyedM2eO\nxsTElJnkp7TvvvtO09PTNSkpSbt37+6ZrU7VmVzo6quv1kGDBmnjxo11xYoVunv3bh08eLAmJibq\nKaecopMnT9YzzjijTHnnnHOOpqSk6NFHH62LFy+usLzSnn32WT3ppJPKbJs1a5ZeeOGFqupMbXr8\n8cdrYmKitm/fXjMyMjzv+/nnn1VEdP78+dq+fXvt27evZ1vJz56ZmanHHHOMJiQkaKdOnXTu3Lmq\nqrpv3z494ogjNDY2Vhs3bqwJCQm6Y8eOw6Z0femll/TYY4/V5ORk7devn3733XeefR06dNCZM2dq\njx49NCkpSS+77DLPoIa+fu+8qerfV7QrLi7WWzJu0eLi4sO27yvY51mftmqaPrD6Ac96VnaW7s2v\nGRNVEaaZ6Vb4sy2UC3AhMA94BjjHy/6KPoCAPfr4o9qtTzftMqSLMhXtMqSLduvTTR99/NGgyqtM\nqKcv3b17ty5dulTz8/M1Ly9Phw0bpkOGDPHsT09P186dO+uPP/6o+fn5mp6erpMmTVJVZ2rRxo0b\n6wcffKAFBQV6ww03aFxcnN9J4vjjj/d8oak6SWLw4MG6d+9e3bJlizZv3lyXL1+uqqqPPPKIHnPM\nMbpt2zbNzs7Wfv36+UwShYWF2rlzZ7377ru1sLBQV65cqQkJCbphwwZVdb7Uk5KSdPXq1arqjO56\n6aWX6qWXXqoHDhzQdevWadu2bT1JYt++fZqamqpPPPGEFhcX61dffaXNmjXzfJmWL6/kS7TE/v37\nNTExsczMgb179/YkmnfffVfXrVunqqrffPONtmrVSl966SVVPZQk/vKXv+j+/fs1Pz9ff/755zI/\n+2uvvaZZWVmqqvree+9pfHy8Z3pUb6PTlh4Z94cfftBGjRrpihUrtKioSO+9917t3LmzFhYWqqrz\n+3PKKafozp07NTs7W4855hjP/1lFv3fl1fUkUXIRuXDpQv3l90MXYQ+sfkBvW3GbZz3nQI4Wugoj\nEWKlQpokgIZACvA1kOx+nQJ0AL4P9ETBLDhPev/Hy/aKPoAypr4zVae+M7XC9Skrp+jiFxdr6qBU\nJQNN7J+oz730nBYXF1d6fDBCPX1peV9++aWmpKR41tPT08sMPf7www/rwIEDVVX1jjvu0OHDh3v2\n7du3T+vXr+/X+aZMmaK9evXyzF2t6iSJ0kNoDxs2TO+55x5VVT3rrLPKJJQ333zTZ5J4//33tXXr\n1mW2DR8+3DO0+qhRo/Qvf/mLZ5/L5dK4uDhPElFVvfXWWz1J4tlnnz3sKnncuHF6xx13eC3Pmyuu\nuEKnT5+uqqobNmzQxMREPXDggNf3Xn/99XrDDTeoqnoSws8//+zZXz5JlDdkyBCdPXu2qlaeJKZP\nn66XXnqpZ19xcbG2bdtW3333XVV1fn8WLlzo2X/zzTfr1VdfraoV/96VV1eTxIP/eVA7n9LZcxHZ\nelBrTemZ4rmILHIVRThC/wWTJCpqkxiH0/5wtPvfkuUl4KFAqrREZL6I/Coia8ttP09EvheRDSJy\ni5dDbwfmeNnut4z0DDLSMypcn9ZvmueJxG6fd0MLnTpgEan0+FCpyvSlBw4cYNy4cXTo0IGkpCT6\n9u1LTk5OSUIFoFWrVp7XpafkLD+VaXx8PE2bNq30nA899BALFizgtddeIy4ursy+li1b+nWu0lOu\nludtitW0tLQyU6iW3v/bb7/hcrnKzMZXuvzNmzfz8ccfe6ZkTU5OZuHChfz6669ey/Nm+PDhPPPM\nMwAsXLiQIUOG0LBhQwA+/fRTzjrrLFq0aEFSUhJz5849bErViqZ8ff311znttNNo2rQpycnJvP76\n6wFNyVr6ZxURUlNTy3xWvv5P/vGPfwT9e1dbuYpd/LjnR8/6aYNOo9VZrTy9j+pRj7l3zGXsKKd3\nWWxM7Z58ymeSUNV/q2pH4CZVPVJVO7qXnqoaUJLAedZiQOkNIhKDk2wGAMcCw0Xk6FL77wZeU9Wv\nAjxXUEqeSIxUw1FVpi+9//772bhxI5999hk5OTm89957gH9DA7du3brM1KX79+9n9+7dFR7z+OOP\nc++997Jy5Upat25d6Tl8nav0NKLltWnTpsx7AbZs2ULbtm0966U/l+bNm1OvXr0yx5SeojU1NZX0\n9PQyU7Lu3buXhx56yGt53pxzzjn89ttvfP311yxatIjLL7/cs+/yyy9nyJAhbNu2jZycHMaNG3fY\n5++r/IKCAi6++GJuvvlmfvvtN7Kzsxk4cKDn+ECnZAWni21FSalE48aNff7e1RWqSlZ2lmc9Oz+b\nK164wvP5n9jmRCacOsFzEZmTl+O5iKwL/HlO4kER6S4iw0TkzyVLICdR1Q+A7HKbTwY2qupmVS0E\nFuG0ReAeG6o/cLGIVEtn8JInEkWEiwZfxMTxE0N+jnBNX5qbm8sRRxxBYmIie/bsISMjw++YLr74\nYl555RU++ugjCgsLmTJlSoXJ5emnn+a2227jrbfeqvBOwJthw4Yxe/Zstm3bRnZ2doVXraeccgrx\n8fHce++9FBUVsWrVKl555RWGDx/u9f0xMTEMHTqUjIwMDhw4wPr163niiSc8+88//3w2bNjAggUL\nKCoqorCwkDVr1vDDDz/4HX+9evW45JJL+Mc//kF2djbnnHOOZ19eXh7JycnExcXx6aefsnDhwjLH\nevtMS7YVFBRQUFBAs2bNiImJ4fXXX+fNN9/0vK9ly5bs3r2bvXv3eo1r2LBhvPrqq7zzzjsUFRUx\nc+ZMGjZsyGmnnVbpz+Tt9y4mJnx98lXD24PQXztyd1BUXARAsRZz4aILyStw7q6axTdj9ZjVZZJA\npC8iI8mfAf6mAg+6l37AvcAFITh3W6D0peIv7m2o6oOq2ltV/6aqXjuDZ2RkeJZVq1aFIJzwu/PO\nO4mPj+eee+7h6aefJj4+nhkzZgCwceNGzj77bBISEjj99NO55pprOPPMMwGnT/z06dNJSUlh1qxZ\nh5V7/fXXs3//fpo1a0afPn0YNGhQmf0VXfF069aNOXPmMHz4cNq0aUPTpk0rvAKdPHkye/bsoXfv\n3iQkJJCYmMjf/vY3n+cqvf5///d/DBgwgJ49e3LSSSdx0UW+hwmIi4tj2bJlvPbaazRr1oxrr72W\np556ii5duvj8mR588EFyc3Np3bo1V155ZZlnURo3bsybb77JokWLaNOmDW3atGHixIkcPHjQZwze\nDB8+nBUrVjBs2LAyX6YPP/wwkydPpkmTJtx5551l5vX2FW/JtsaNGzN79mwuueQSUlJSWLRoERde\neKHnfUcddRTDhw/nyCOPJCUlpcx85ABdu3ZlwYIFXHvttTRv3pxXX32VZcuWUa9ePZ/nLuHt965v\n374BfSaB8Db0RCj5SkK/5//O/sL9nvXhzw9nU7Zz4RUbE8vaq9fSuH5jn+VWx0VkOKxatarMd2VQ\nKmu0AL7BSSZfu9dbAm8F2vgBpAFrS61fBMwrtT4SmO1nWRU1yhhjwqAqf1+Pzn9UjzrtKE/jb+r5\nqZrWO83T+LtpzybdsOtQp4O8g3llupX6q6QH0jMvPFOm2+llSy7Tt396O+j4awtC3HBd4oCqFgNF\n7qew/wdU3MLnn21A+1Lr7dzbjDG1TGLvRNr3b+9p/D1YeJDLrrzM0/j7ybZPWPXzKs/75385n/s+\nvM+zPvuT2dz9wd2e9WfXPcsTXx2qTrzxvhtJ653mGf/o6oev5ug+RzM3cy4AC4cupP+R/cP8U9ZO\n/sxMt0ZEkoD/4PRuygNWV3yIV0LZ0WM/AzqLSBqwA7gM8F7h7EVGRgbp6emkp6cHEYoxJpx+2fsL\nGasyPMNWX3LsJdQ7qR5j3htDt8+7sfXAVnq37e2pCrus+2Vljh9/yvgy66N7jcalh55UP6H1CWXW\nhwwbQkJSAo+/+DgIJMQlMGvqLM/Ip3WlkdmXVatWBV0tLxpAA5KIdAASVXVtJW8tf9xCIB1oCvwK\nTFXVTBEZCPwLpzprvqre7buUMuWpt7hFJOINYsbUVhX9fakqi79dzMXdLiY2Jpai4iJWbFrBuZ3O\n9XxB3/Xvu+h6ZFeGnj+Upa8sZWPWxpDW7S95eQlX3n8lqU1S2fr7VjJvygx6eOzayv1/GFDG9Jkk\nRORoVf1eRE7wtl9VvwgixpCwJGFM9Sv/91XoKkRR6sfWB+C6165jSt8pNG/UPCLxhTsJ1QahThLz\nVHWsiHjrNK2qelYwQYaCJQljql/5v68RS0cwvPtwzu96fgSjMoEIaZJwFxgDnKaqH1Y1uFASEZ06\ndephbRKWJIwJHxFh3pp5/N+JzuzFBa4Cz12EqdlK2iSmTZsW2iQBICJfqurxVYowxHzdSTRq1Ij9\n+/d7OcIYU1Xx8fF8kvUJ3Vt0j3QoJkjB3En407tphYhcBCz1+s1cg+zbty/SIRgTleZmzmX2Y7Mp\naF7Aj71+pMvXXdAdStof0nh75tuRDs9EkD93ErlAI6AIyMfpxqqqmhj+8HzGVNPzlTFRRVVZ9OIi\nrnjgClz9XaR+msqscbM8Txmb2iGYOwl/xm5KUNUYVa2vqonu9YglCGNM6Ix+aTTrf1uPiBAXG0e8\nxNfJQeyMb/5UNyEiyUAXnDkmAFDV98IVlD/sYTpjArf0u6W0SWjDqe1OBeD6U66nQ1IH4NAgdqW7\nkJraIawP04nIVcAEnGEzvgJOBVbXxC6wxpiyftn7C7v376Znq54ALP9xOS0ateD41jWqL4qpJmGp\nbsJJEL2BzaraDzgeyAkiPmNMNcg+cGhU/m9+/YZ3N7/rWR/QeYAlCBMQf5JEvqrmA4hIA1X9Hjgq\nvGEZY7zRSuZj+Hrn1wxdPNSzPrDLwMPGQTImEP4kiV/cA/y9CLwlIi8BvqcUM6YOq+xLvKrKz8eQ\nX5TPSfNO4mCRMy9Gj5Y9WPHnFWE5t6mbAh3gry/QBHhDVQvCFlXlcXh94tqYSCsZZK4qg8vtL9xP\n/dj61Itx+pX8sOsH3lj6BvOemEdhi0I29tzIkV8eScNdDRl/1Xj6XdCPLildrCeS8SncT1zPBhap\n6kdVCTKUrOHa1DQlD6OVfIl3/KIjR+w+gvFXjSe7SzZDjh7C0c2cKdyvfe1axhw/xtM2MGTREG4+\n/Wb6pPYBYMCCAUztO9WzPualMVzT+xp++uwnbpx3I1tP3krbT9ryr7/+y55jMAEJ1xPXnwO3i8hR\nwAs4CWNNMAEaUxs9880z9DqnFxkpGdw470YQpxro3lvu5aLBF/HGj2+Q1DDJ8/6xJ44lrcmh+cHn\nXzCfhAYJnvXlI5eXKX/+hfMB2CSbyMnLceZj2LfVnmMw1cKfh+meUNVBOD2cfgDuERHrQG3qrHX/\nW8fqrYfm3Uo+IplG9RshIp4v8bx9eZ4v8YFdBtKqcSvP+3u07EGThk08603jm/o1UF7JcwzrXlpH\n5k2Z9hyDqRZ+t0mIyMnApcCFwHeqOjicgVUSi1U3mWqz9+Befs75mR4tewDw5k9vkpOfw7Bjh5V5\nn81nYGq6kA8V7i70XuBPwE/AIuBFVY3ocxKWJEy4/W/f/2jRqAUAn277lEXrFjFrwKwIR2VM1YTr\nYbqfcOYFvHoeAAAgAElEQVSUOE9V/xvpBFEiIyMj6MfMTd3lq4tqsRZ7Xu/I3cHZT57tec/JbU+2\nBGGi2qpVq8jIyAjq2IC6wNYUdidhguWti2qxFtPlwS58PvZzTwOzqlqjsKl1wlLdVBNZkjCBKt9F\nNeHDBJrvbc7Nf72ZcaPHkZOfU6YHkjG1Ubi6wBoT9caOGktKSoqni2qjeo2Yfst0hv9pOIAlCGN8\n8Heo8J7AGe7V91X16/CFZEzobc/dzoz3Zhx6zmD/VhrENbAqJWMqUWnDtYhMAJ4GWriXBSJyXbgD\nMyaU2iS0IT0l3Z4zMCZA/nSBXYvTu2mfe70RznwSPaohPl8xWZuE8cuu/btoFt8s0mEYUyOEqwus\nAK5S6y73NmNqtJz8HM564izyi/IjHYoxUcufNolM4BMRecG9PgSYH76Q/GPTl5rKJDVMYs3YNX4N\neWFMbRbW6UsBRORE4HT36vuq+mVQZwsRq24yvqgqr2x4hfO7nm+N0saUE67qJnDmtl6CM/HQbhFp\nH2hwxlSHg66DvPD9C+QW5EY6FGNqBX8arq8DpgK/cqg9Qq3h2hhjoku47iQmAEep6rGq2kNVj4tk\ngjCmPFexiwmvT2DX/l2RDsWYWsefJLEV+D3cgRgTrNiYWP7Q/g8kNkiMdCjG1Dr+VDfNB44CXgUO\nlmxX1YgNi2nVTQacO4jYmNhIh2FM1AhXddMW4C2gPpBQajEmYlSVc546hw27N0Q6FGNqNRsF1kSt\nbXu30TaxbaTDMCZqhLMLrDE1QlZ2lmeCIEsQxoRf1CYJm5mubrrl7Vv4YscXkQ7DmKhiM9OZOsNm\njDMmeGGpbhKRdiLygoj8JiL/E5HnRaRd8GEaE5iXvn+JX/b+AmAJwphq5k91UybwMtAaaAMsc28z\nplpsy91G9oHsSIdhTJ3kz3MSX6lqr8q2VSerbjLGmMCFa47r3SIyEnjGvT4c2B1ocMb4Q1WZdMck\nmpzVhJ6tejKoy6BIh2RMneZPddOVwDBgJ7ADuBgYFcaYTB32/LLneXjlwzTc3JCT254c6XCMqfP8\nqW46XVU/rGxbdbLqptpnbuZc7nn0HmJbx/Jjrx/p8nUX4v4Xx/irxjNu9LhIh2dMrRCuh+ke9HOb\nMUEbO2osRw88mryDeSCQX5DPtFumMXbU2EiHZkyd5rNNQkROA/oAzUXkhlK7EgEbVc1U2e79u/lq\n51f0P7I/IsKVJ1zJle9cSbfPu7E1bysiYl1ejYmwiu4k6gONcRJJ6YH99uK0SxhTJbsP7GZl1krP\n+sasjWTelMm6l9aReVMmG7M2RjA6Ywz41yaRpqqbqykev1ibRPR6eu3TDOwykJQjUiIdijF1Tlja\nJGpagjDR7X/7/sfu/daD2phoYWM3mbD6ac9PvL/lfUb1GhXpUIyp8+rUUOE2Cmx0OCLuCApdhZEO\nw5g6LayjwIpIV+ARoKWqdheRHsAFqnpnUGcMAbuTqNnufO9ORvUaRbtEGwfSmJokXHcS/wEmAYUA\nqroWuCzw8Exd0a15N+Ji4iIdhjEmBPxJEvGq+mm5bUXhCMZEp692fsUd797hWR96zFBaNm4ZwYiM\nMaHiT5LYJSKdAAUQkYtxxnAydZCqMnHaREpX96U1SeP01NMjGJUxJlz8aZM4EpiH8/R1NpAFjIhk\n11hrk4icJS8v4cr7r6T3mb2ZM2EORzc7OtIhGWP8FEybhD9JIlZVXSLSCIhR1dyqBBkKliSq39zM\nucx+bDaFLQrZ2HMjqZ+n0mhPI66/6nobgM+YKBGuhussEZkHnArkBRWZiXpn/PEM0s5JI78gHwQo\nhum3TLcB+Iyp5fxJEkcDbwPX4CSMh0TkD+ENy9Q0R6YcycntTiYnL4dun3cjJy/HBuAzpg6odGY6\nVd0PLAYWi0gy8G/gXWwk2Frv2XXP0jG5Iye3PZmG9RrSYF8DMm/KZOj5Q1n6ylIbgM+YKHD32LHk\nb9gQ9PF+DcshIn2BS4HzgDXAs6r6fNBnrSJrk6gey39cTpuENhzX8rhIh2KMCVJGejoZ774LODXF\nIZ/jWkR+Br7EuZv4h6ruCzxMEw0KXYUsWreIkT1GIiIM6Dwg0iEZYyLMnzaJHqr6J1V9xhJE7SYi\nfLnzS/IKrH+CMcZR0cx0N6vqvcAMETmsbkdVx4c1MlMt9h7cy5bft9C9RXfqxdRj1oBZkQ7JmBrF\nV51+w65dmThvXgQiql4VVTd95/53TXUEYiJj9dbVfLj1Q7q36B7pUIypkfI3bPDU6ZeWUf2h+Ofz\nz+Hee+HZZ0NSnM8koarL3C/3q+pzpfeJyCUhObuJiJ15O2ke35zYmFgGdB5gbQ/G1CbHHgv33edZ\nbdi166GE5iXZVcafNolJfm4zUeK616/j023lx2w0xpShCp9+6vzrj127YP/+8MZUoqAAXngBhgyB\nX38tu69hQ2jf3rM6cd48MlatIiPI+XcqapMYCAwC2orI7FK7ErFRYKNOflE+Des1BODZi58lRqJ2\nviljwu/FF+Gf/4ScHGjWzL9jnnoKJk+GXr2cfweE8Q794ovh999h1ChISAjfeai4TWI7TnvEBcDn\npbbnAn8PV0Ai0hG4DUhU1WHhOk9dsj13OwOfHsiX474kRmIsQRhTmawsuPVWGDwY+vf375i//x3G\njoXVqyE19fD933/vbG/UqMzmShvGVaH8yAaLFzt3DNWgojaJr4GvRWShqlbb/JOqmgVcJSKLq+uc\ntYmqMumOSdw15S7A6dbaJqEN749+35KDiZga3UPI5YLYcgNI/P3QdXCZOv1SGnbtevjGRo3g7LO9\nn+f++50k0rt3mc0VNoy//josWQLz55c7efUkCPDjYTqgg4jcBXQDPJGp6pH+nEBE5gPnA7+qao9S\n288D/oXTLjJfVe8JJHDj3fPLnufhlQ+TnZhNp5M7cfPpNwOQ2CAxwpGZuqzG9hD6/nvni/u993y+\nJWRJ7D//OXybKnzzje9jzjgDTjstNOcPkj9JIhOYCjwA9ANG41+Dd+njHwSeLNkgIjHAQ0B/nGqt\nz0TkJVX9vtRxNnJcAOZmzuVfj/0LVwsXuf1yeXvV27z73Ls0+b8mNpQ3NfxKtgaos5/PUUfB8xEb\nYcjRvj3s2eN9X+PG1RuLF/4kiSNUdYU4AyZtBjJE5HNgij8nUNUPRCSt3OaTgY0lExeJyCLgQuB7\nEUkBZgC9ROQWu8Pwz4jLRzB59WQabGsAAoVFhdwz8R4uGnxRpEOrEWrslWwNEfTnk5MDO3Y4PXva\nt4e0cn/qs2dDhw6+j9+xA375Bdq1CzDiALlc8PLLcNxx0Lnzoe0i0Lx5eM9dERFo0iRy5/eDP0ni\noPvKf6OIXAtsA6qa3toCW0ut/4KTOFDVPcDVlRWQkZHheZ2enk56enoVQ4o+j33xGOd2Opf2TdrT\nuEFjHjjvAa7+99V0+7wbW/O22lDeJnQWLHC+zAYPLrt91ix47jlo2hRuvPHwJHH22ZBYQVVnTo7T\nSFw+SaxcCd26QatWfoVX4Z3Q0KEwfjwkJ8NDD/lVXm2xatUqVgXZ9bWEP0liAhAPjAemA2cBf6nS\nWUOgdJKoK1SVA0UHiI+LB6BBbAMKXYf6FGzZusWG8jbBOXiw4v3HHgtHHHH49jvucBZfunWruNxj\njnHq3cv76CNo0eLwJLFsmdPls1s35w7AfRFU4Z1Q69bw2GPOeWrgRVNADeMBKn8BPW3atIDL8Gc+\nic/cL/Nw2iNCYRvQvtR6O/c2U4GHP3uY7bnbmdF/BgBX9LyizP5JEw4942jVTCYgBQUV7z/++CoV\nH/AX4e23e9/+ySfOU8PffuskkqP9mGO9Z09/w4yImt7mU9HDdMsAn48aquoFAZxHKNsQ/RnQ2d1W\nsQO4DBgeQHl1wo7cHSz+djETTp0AwFUnXEX92PoRjqoWKipy+sT/7W8V15/XZmF+ICtkX4R33un8\na/PJVJuK7iRmhuIEIrIQSAeaisgWYKqqZorIdcCbHOoC+10FxRwmIyOjVrZF5B7MJaGB8wdb8m+J\nBvUaRCKk6FVUBNOnw7XXVnwlW1jo9Ds/8UQYN8550rY2e/pp5wr8xBM9m8JZ5REWNbDaqCarStuE\nvzPTHQG0V9UfgjpLiEXzzHSlH3Yr36isqvR8tCevXP4K7Zu091GC8cv+/TB8OBw44HRx9OdK+bff\nYO1a/5+wjVavveY0MB97bKQjCZnSs6+V2d63b9BjFtVGIhKWmekG49xV1Ac6ikgv4I4Aq5uMW8nD\nbr1P6M1Fgy/i+fXP0zqhNX1S+yAirBm7xqqUQmH1aqc3y3PPQX0/P8/mzb0niOJiiKlFT6sPGhTp\nCEIu6u6EokildxLuZyLOAlap6vHubd+oasQmPo7GO4m5mXOZ/dhsClsUsrHnRrp83YW4/8VxzgXn\nMOovo+jVqlekQzTeuFzOMArLl0e2P30wVq92qs6eeAJSUiIdjakBwnInARSq6u/lqkYi/g0dbW0S\nY0eNJSUlheseuQ4E8gvy+ect/+SiwRfZsww1WWys0+0y2hLE9dfD0qWQkVHxcwqmTghrm4R77KUV\nwETgIpznJeJU9a9BnTEEovFOAmDJy0u48v4raZvYlm17t5F5U6Z1VQ2VtWuhe/fqqxbaudN5gCwu\nrnrOF6jvv3d6alXjQHCm5gvXncR1OEN3HwQWAsuBOwMPr277Pf93NmZttIfdwkEVZsxwHuo66qjq\nOefs2c5wzRkZMHw4d199deTGPjp4EBqU6/nmz/MDxvihwjsJEYkF7lHVm6ovpMpF251EoauQE+ed\nyMq/rKRZvJ8TmJia7513YOpUeOwxMsaODWvvGp/DTnTqxMT334cvvqgRg8GZmi3kdxKq6hKRP1Qt\nrPCIpjaJuNg4Pv2/Tz0zw5laol8/Z6kGFQ47sWaNJQhToaq0SfhTgfuliLwsIleIyNCSJaizhVBJ\nkqjJfs//HVexC8ASRCgVFMD69ZGOInAvv+xUA/Xo4TzhXd777zvVWOV99ZUzCJ4v1jBtKpGenh70\neHf+tEk0BHbjdIMtocDSoM5Yh0x7dxontD6BkT1GRjqU2iM315nft3Vr+O9/Ix1NYNLTnbmTCwq8\nP9zXvr33hvB69ewJYxMx/gzwF6pB/eqcu8++m7iYGtr7JVpddpnTa2fOnEhHErjExIqv+tPSDh9q\nG5xeWx06wObNYQvNGF/8uZMwAVBV9hzYQ9P4pvbkdDj85z/OXUQNu7K2J35NbWVJIsQ+/uVj7v3o\nXl649IVIh1I7tWkT6Qi8Cnc3V0tCJlL8GuCvphERnTp1ao3t3VToKiQu1qqZQqK2jZtkTASU9G6a\nNm1awF1g/XniuiXwT6CNqg4UkW7Aaao6P/iQq6YmPifxa96vtGzcMtJh1C5ffOE8rPbyy5GOxJha\nIVxPXP8XyMR56hpgA/AsELEkUdMUugoZ+PRA3rziTXtYLkheHxZTpWG7dkyMTEjGGPxLEs1UdbGI\nTAJQ1SIRcYU5rqhS8rBcvRhr4gmWz4fF+vaNQDTGmBL+VPbuE5GmuEd+FZFTgd/DGlWU+G3fbxS4\nnLmBLUFUUWVzLBtjIsKfJHED8DLQSUQ+BJ7EGfSvzrv3w3t57tvnIh1GdFOFyy+HTz+NdCTGGC/8\neZjuCxHpCxwFCPCDqhaGPbJK1ISxm+4++25ixHreVIkIjB4NW7bAhx9GOhpjaqVwzycRC/wR6ECp\npKKqs4I6YwhEsndTsRazI3cHbRPbRuT8Ue3tt52Z3gYMOGyXzVFsTPiFq3fTMiAf+AYoDiaw2uTT\nbZ8y86OZLBm2JNKhRJ/4eOe5By/sYTFjaiZ/7iTWqmqPaorHL5F+TqJYi62aqSLbtzsjml56aaQj\nMcaUEsydhD/fdK+LyLlBxlRr/Jzzs+e1JQgfCgth8GA49lhYtcpplDbGRDV/7iT+BCzASSiFOI3X\nqqoRG8S+uu8kCl2FnJF5Bq9c/kqdfVjO58xo5afnfO016NsXGjWqxuiMMf4IV5vELOA04JsaNxZG\nNYmLjeOjMR/V6TuICmdGK23QoOoIxxhTTfxJEluBdXUxQWzP3U5yw2SOiDuibiYIl8uZw8BHY7Mx\npvbzJ0lsAlaJyOvAwZKNkewCC+F7TkJVmXTHJO6achcPffoQPVr24LLul4X0HOHgd3WQNy4XxMaW\n3bZkCfz5z9CsGfz1ryGM1BhT3arynIQ/SSLLvdR3LzVCsPO1Vub5Zc/z8MqH6X1Cb2acPwOpYZPb\n+FJpdZAq/P47JCWVfcMbb8Ajj8BLL5XdPnAg/PbbobaFN98MdcjGmGpSckE9bdq0gI/154nrwEuN\nQnMz5zL7sdkUtigkt18ukx6fxJS7pzD+qvGMGz0u0uEFb/9+OOkk+OEH6NMHli8vu//ss+FcL53X\nrOHZGEMFSUJE/qWq14vIMtyD+5WmqheENbJqNnbUWFJSUhj/6HgQyC/I55+3/JOLBl8U6dD846vd\noEEDmDkTjjrq8LsIgHr+DUxoD7sZUzdV9A3xlPvfmdURSKSJCCLCvv376PZ5N7bmbfVsiwq+4oyN\nhVNOqXLx4Z6e0xhTM/lMEqr6uftlL1X9d+l9IjIBOLwCPMptzNpI5k2ZDD1/KEtfWcrGrI2RDsk7\nlwvy88tWCUVLMjPGRBV/6hr+Avy73LZRXrZFvfb92nNKh1MQkZpZzZSTA48/DnPmwDXXwA03eHZZ\ndZAxJhwqapMYDlwOdBSR0pMMJwB7wh1YJOzM20msxFb+xkhZvRrWrIEFC+DUU8vssuogY0w4+ByW\nQ0TSgI7AXVBmmuFcYK2qFoU/PO8iPcBftVC1KiRjTEiFdFgOVd0MbMYZkqPGqQmTDoWNqnOnsGQJ\npKZGOhpjTJQL96RDQ4F7gBY4g/vVugH+Xt3wKh9t/YgZ/WeErMwq27YN2trERsaY0AnmTsKfJPEj\nMFhVv6tKcKEU6iSRV5DH9tztdG0avkZer8NmqNKwQQMmjh0LF18ctnMbYwyEbxTYX2tSggiHxvUb\nhzVBQAXDZiQkwHXXhfXcxhgTLH+SxBoReRZ4kbID/C0NW1TVaEfuDlontI5cACec4EzUY4wxNZA/\n418nAvuBc4HB7uX8cAZVXVzFLgYsGMDu/bvDe6L9+50H4IwxJsr4M8Df6OoIJBJiY2L5+q9fh3/o\njUmTYNeu8J7DGGPCoNIkISKZeB/g78qwRFTNqmVsplmzoH9/+P778J/LGGNCyJ82iVdKvW4I/AnY\nHp5wqs8nv3wCwCntqj74XRmrVzvPOJROPrGxNmyGMSYqVdoF9rADRGKAD1S1T3hC8iuGKneBfW3j\nawAM6hKiOZk//xwmToSff4b33oPWEWwMN8YYL8LynISXkxwFvKqqnQM6MIRq5LAcM2ZA06YwZgzE\nxUU6GmOMOUy4HqbLpWybxE5gkqo+H3iIoVEjk4QxxtRwwSSJCrvAitOqe6yqJpZaukYyQVRVTn4O\ng54ehKs4yC6phYXOTG/WpdUYUwdUmCTcl+uvVlMs1aJx/cZkpGcQGxPkkOD16jkJYv/+0AZmjDE1\nkD+9m74Qkd6q+lnYowlAsKPA1oupx8ltTz5su9exlXB6H5WZq0EEbrkl0HCNMSZiwj0K7PdAZ5xh\nw/dxaBTYHkGdMQSCbZPYvX83KUekeH02IiM93fvYSqmpZGzZElScxhhTk4RrgL8BQcZT49y+8nb6\nH9mfi7sFMOJqUlL4AjLGmBrOn2E5NldHINXh4T8+jB7+8HjFUlLCE4wxxkQBf+4kag0RQbApQY0x\nxl/+jAIb9XLyc1j87eJIh2GMMVGnTtxJ/LbvNzZlb6rwPQ2Tksg488yyYy5hYysZY+q2gIflqAlC\n/sS1Kgwb5ozWmpoaunKNMaYGqZaxm2oCG5bDGGMCF/JhOWqDccvG8cOuHyIdhjHGRKVanySu6HkF\nHZI6RDoMY4yJSnW3umnHDtiwAfr2DU1QxhhTw1l1UykHCg9wsOig953FxTBqFLz9drXGZIwx0abW\nJolXNrzC+NfHe9/5yCPw++8wZUr1BmWMMVGmVlc3FboKiYv1Mkvc9u1QUAAdOoQ+OGOMqaGsC6wx\nxhifrE3C7blvn+NA4YFIh2GMMVGv1iWJQlchb296mxgp96PZnYcxxgSsblQ3HTjgdHVdsQISEsIX\nmDHG1GC1ok1CROKBh4GDwLuqutDLewJvk8jKgo4dQxKjMcZEo9rSJjEUeE5VxwEXBHLg/C/m89y3\nz3nfaQnCGGMCFvYkISLzReRXEVlbbvt5IvK9iGwQkVtK7WoHbHW/dvkqd+TIkWRlZZXZdkbaGfRq\n1SvoWLOyshg5ciT9+vXzWn5VWflWvpUfneVHc+ylyw+KqoZ1Af4A9ALWltoWA/wIpAFxwFfA0e59\nI4BB7tcLfZSpgHbq1Ek3bdqkXhUUqBYVed/nxaZNm7RTp05aUnal5QfIyrfyrfzoLD+aY/dWvgb6\nHR7oAcEs7mRQOkmcCrxean0icIv7dTzwODAHGO6jPM8PPGLECFVVzTuYV/aTufFG1YwMvz/IESNG\nlPlPKl9+VVn5Vr6VH53lR3Ps3srXAL+/q6XhWkTSgGWq2sO9fhEwQFXHutdHAierqo9xNA4rr2zQ\nDYExOM3dCv2B/wI9gT0h+hmMMaY20FrQcB2QESNGoAeU/H/lo8VO5nv7/vtp9+ab7A4gW44YMcJ3\n+SG4m7LyrXwrPzrLj+bYKyrfb6EIorIF79VNb5Ra91Q3+Vle1NcLWvlWvpUfHeVHc+zeytdAv78D\nPSCYBegAfFNqPZZDDdf1cRqujwmgPD3uuON04cKF+sHmD3T73u0h+zBHjBih/fr10xEjRoTsP6l8\n+b169Qpr+RZ/xeVb/BWXb/H7LjsaY1dVXbhwoR533HE1M0kAC4HtOA/HbQFGu7cPBH4ANgITAyzT\n88Pf/f7d+skvn4Tu06wGU6dOjXQIVWLxR5bFHznRHLuqBpUk6gVcPxUgVb3cx/bXgderWv4tf7jF\nGfZ77ly46iqIja1qkcYYY9yivuEagNxc2LkTYmrHj2OMMTVFjRu7yR+HdYE1xhjjF432Af6MMcbU\nHFY/Y4wxxidLEsYYY3yyJGGMMcanqEsSFQwxXuOJSDsRWSki34rINyLi11hVNYmIxIjIFyLycqRj\nCYaINBGR50TkO/f/wymRjslfIvJ3EVknImtF5GkRqR/pmCribZoAEUkWkTdF5AcRWS4iTSIZY0V8\nxH+v+3fnKxF5XkQSIxljRXxN0+Ded6OIFItISmXlRFWSEJEY4CFgAHAsMFxEjo5sVAEpAm5Q1WOB\n04Broix+gAnA+kgHUQX/Bl5T1WNwxoD8LsLx+EVE2gDXASeoM1BmPeCyyEZVqUycv9XSJgJvq+pR\nwEpgUrVH5T9v8b8JHKuqvXAeBI62+BGRdsA5wGZ/ComqJAGcDGxU1c2qWggsAi6McEx+U9WdqvqV\n+3UezhdU28hG5T/3L9cg4LFIxxIM91XfGaqaCaCqRaq6N8JhBSIWaCQi9XCG1N8e4XgqpKofANnl\nNl8IPOF+/QQwpFqDCoC3+FX1bVUtdq9+jDNJWo3k4/MHeAD4h7/lRFuSaMuhWesAfiGKvmRLE5EO\nOJMxfRLZSAJS8ssVrf2mOwK7RCTTXWU2T0SOiHRQ/lDV7cD9OEPbbANyVPXtyEYVlBaq+is4F01A\niwjHUxVXEoJRI6qTiFwAbFXVb/w9JtqSRK0gIo2BJcAE9x1FjScifwR+dd8JiXuJNvWAE4A5qnoC\nsB+n+qPGE5EknKvwNKAN0FhEvA55E2Wi8oJDRG4DClV1YaRj8Zf7guhWYGrpzZUdF21JYhvQvtR6\nO/e2qOGuKlgCPKWqL0U6ngCcDlwgIpuAZ4B+IvJkhGMK1C84V1Fr3OtLcJJGNDgb2KSqe1TVBSwF\n+kQ4pmD8KiItAUSkFfC/CMcTMBEZhVPtGm1JuhPOiNxfi0gWzvfn5yJS4d1ctCWJz4DOIpLm7tlx\nGRBtvWweB9ar6r8jHUggVPVWVW2vqkfifO4rVfXPkY4rEO5qjq0i0tW9qT/R0wi/BThVRBqKiODE\nHg2N7uXvOl8GRrlf/wWo6RdKZeIXkfNwqlwvUNWDEYvKf574VXWdqrZS1SNVtSPORdPxqlphoo6q\nJOG+groWp4fBt8AiVY2GPxQAROR0YARwloh86a4XPy/ScdUx44GnReQrnN5N/4xwPH5R1U9x7ny+\nBL7G+cOfF9GgKiEiC4GPgK4iskVERgN3A+eIyA84ie7uSMZYER/xPwg0Bt5y//0+HNEgK+Aj/tIU\nP6qbbOwmY4wxPkXVnYQxxpjqZUnCGGOMT5YkjDHG+GRJwhhjjE+WJIwxxvhkScIYY4xPliSMMcb4\nZEnCGGOMT5YkTFRwD8Xi98iVVTzXeBFZLyJPhaHsDyrZ30RErg7kmFASketFpNA9rhIi0kdEPheR\nkdUVg6lZ7IlrExVEJA1Y5p5wJ5jjRf38ZReR74D+7uG5Q8afGNxDyC9T1eNCeW5/iUg6ztA3H6jq\nv9zbLlHV5yIRj4k8u5MwIeW+4l/vnqthnYi8ISINyt8JuKdPnOLe/p17jocfRGSBiPQXkQ/c6yeV\nKj7OvX+9iCwWkYbuskaIyCfusXQeEUeaONPcPuE+72GTw4jIDeJMI7tW3FPJisgjwJHA6yIyodz7\n7xKRv5Van+r+OV4Qkc/cZV1V7rMoHUOqiOS693k9BrgLONL9s9zjfm/JMaXjnVDR5+3eFy8ir7jH\nCVsrIpf48V/YAWco6ZHuMhoD0TQxkwk1VbXFlpAtOPMdFADHudefxRlSOQ1YW+p9NwJT3NsLgW7u\n7WuAx9yvLwBeKFVuMXCqe30+cANwNM7IorHu7XNwvuDSABfQ20ecJ+AMlNcQaASsA3q6920Ckr0c\n0wtYVWr9W5xJr5Lc6w2Bb0qOdcdQVDoGYK/734qOWVvuvHt9xevr83a/HgrMLVVOgh//f392//sZ\ncMJ+c4oAAAL2SURBVBTQF2gW6d8rWyK32J2ECYcsPTTz1ec4V6cV2aSqJUN2fwuscL/+BudLsMQW\nVf3Y/XoBcAbOSKInAp+JyJfAWTh3AgA/q+pnPs75B5wElK+q+3DmZzjDvc/rpErqTLjUXERaiUgP\nYI+qbgOud48qWzKdZZdSh232EUNFx5QnOPN5+IrX1+f9Dc6Iq3eJyB9UNbeCc5QomZrzCZxk20JV\nd/lxnKml6kU6AFMrlR5n34Vz9VuEM0dziYY+3l9car2Ysr+j5evzS77Q/quqt5Xe4W7D2BdY2H55\nDrgEaAU8KyJ9cRLTKap6UETeoezPdlgM7mP6V3CML76Gdfb2eaOqG0XkBJwJcu4UkbdV9U6fhTuN\n1SXTAz8DfIhzx2LqMLuTMOHg7cvsV5yr8GR3nfn5lbzf2740ETnF/fpy4ANgJXCxiDQHcJff3sux\n5b0PDBFnEp9GwJ+A9yp4f4nFOJMuXYSTMJoA2e4v+6OBUyuIv2Q9EecuxNsxuUCCn/G+7+MczkaR\n1sABdabYvA/3LHwi8rZ7X3kn4VQzoaq7cSY1Kvlc/yYiY0XkchHpLSLnisiJIjLGvX+ciFwnIoO8\nxWKil91JmHA4rAePqhaJyHScL6FfKDurmvp4XX79e+AaEcnEqZZ6RFXzReR24E0RicGpn78GJyn5\n7Emkql+KyH/d8SgwT1XX+oq/1HHrRSQB+EVVfxWRN4C/isi3wA/A6griL1l/A7ja2zGqukdEPhSR\ntcDrqnqLs1m/8hLv1+47Jl/xHgfcJyLF7s/lryIiONNY7in9RhHpB2Tg3IUscW9+HNjtfv2t+zwf\n4zSu/wunPaakmux7IEZV3/ERi4lS1gXWmDpERI4FRqvqTQEe1xenk8BvQBbOncwSYDhO+9DpQLy6\nu82a2sOShDGmUiIyFnCp6vxIx2KqlyUJY4wxPlnDtTHGGJ8sSRhjjPHJkoQxxhifLEkYY4zxyZKE\nMcYYnyxJGGOM8cmShDHGGJ8sSRhjjPHp/wEnpKUiu1YZBgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106b94350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "matplotlib.rcParams['legend.numpoints'] = 1\n",
    "fig = plt.figure(figsize=(6,4))\n",
    "ax = plt.subplot(111)\n",
    "#ax.set_xlim(extent[0],extent[1])\n",
    "ax.set_xlabel(\"number of variations, $N_{\\\\rm{par}}$\")\n",
    "ax.set_xlim(0,Nmax-1)\n",
    "ax.set_ylim(0.8,550)\n",
    "#ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "\n",
    "\n",
    "ax.set_ylabel(\"runtime relative to no variations\")\n",
    "ax.plot(range(Nmax),np.ones(Nmax),\"o-\", color=\"black\", label=\"no variations\")\n",
    "ax.plot(range(Nmax),var_1/var_0,\"s-.\", color=\"red\", label=\"1st order variations\")\n",
    "ax.plot(range(Nmax),var_2/var_0,\"*:\", color=\"green\", label=\"1st and 2nd order variations\")\n",
    "\n",
    "legend = plt.legend(loc=2)\n",
    "ax = plt.gca().add_artist(legend)\n",
    "plt.savefig('paper_test5.pdf',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}