timvieira/Exponential-jumps.ipynb

## Exponential-jumps.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiments with the exponential jumps algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jump(stream):\n",
    "    \"Weighted-reservoir sampling by jumping\"\n",
    "    R = None\n",
    "    T = np.inf\n",
    "    J = 0.0\n",
    "    for i, w in enumerate(stream):\n",
    "        J -= w\n",
    "        if J <= 0:\n",
    "            # Sample the key for item i, given that it is smaller than the current threshold\n",
    "            T = Exponential.sample_truncated(w, 0, T)\n",
    "            # i enters the reservoir\n",
    "            R = i\n",
    "            # sample the waiting time (size of the jump)\n",
    "            J = Exponential.sample(T)\n",
    "    return R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def walk(stream):\n",
    "    \"Weighted-reservoir sampling by walking\"\n",
    "    R = None\n",
    "    T = np.inf\n",
    "    J = 0.0\n",
    "    for i, w in enumerate(stream):\n",
    "        X = Exponential.sample(w)\n",
    "        if X < T:\n",
    "            R = i   # i enters the reservoir\n",
    "            T = X   # threshold to enter the reservoir\n",
    "    return R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, pylab as pl\n",
    "\n",
    "def base(w):\n",
    "    \"Baseline algorithm for sampling from weights w\"\n",
    "    c = np.cumsum(w)\n",
    "    u = np.random.uniform()\n",
    "    return c.searchsorted(u * c[-1])\n",
    "\n",
    "def empirical(gen, s, n, reps):\n",
    "    \"Compute empirical categorical distribution\"\n",
    "    q = np.zeros(n)\n",
    "    for _ in range(reps):\n",
    "        q[gen(s)] += 1\n",
    "    q /= reps\n",
    "    return q\n",
    "\n",
    "def error(p, q, reps):\n",
    "    \"Compute total variation between two (possibly unnormalized) distributions `p` and `q`.\"\n",
    "    assert p.shape == q.shape\n",
    "    p = p / p.sum()\n",
    "    q = q / q.sum()\n",
    "    return 0.5 * np.abs(p-q).sum()  # total variation\n",
    "\n",
    "class Exponential:\n",
    "    \n",
    "    @staticmethod\n",
    "    def pdf(w, x):\n",
    "        return w * np.exp(-x * w)\n",
    "    \n",
    "    @staticmethod\n",
    "    def cdf(w, x):\n",
    "        return 1 - np.exp(-x * w)\n",
    "        \n",
    "    @staticmethod\n",
    "    def ppf(w, u):\n",
    "        return -np.log1p(-u) / w\n",
    "\n",
    "    @classmethod\n",
    "    def sample(cls, w, u=None):\n",
    "        \"Generate a random variate.\"\n",
    "        if u is None: u = np.random.uniform()\n",
    "        return cls.ppf(w, u)\n",
    "\n",
    "    @classmethod\n",
    "    def sample_truncated(cls, w, a, b, u=None):\n",
    "        \"Generate a random variate such that `a <= X <= b`\"\n",
    "        assert a <= b\n",
    "        if u is None: u = np.random.uniform()\n",
    "        return cls.ppf(w, cls.cdf(w, a) + u * (cls.cdf(w, b) - cls.cdf(w, a)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Does it work?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate some test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 20\n",
    "w = np.random.uniform(0, 10, size=n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "base error: 0.005302418883592938\n",
      "walk error: 0.005337774053666295\n",
      "jump error: 0.005126327265828854\n"
     ]
    }
   ],
   "source": [
    "reps = 100000\n",
    "print('base error:', error(w, empirical(base, w, n, reps), reps))\n",
    "print('walk error:', error(w, empirical(walk, w, n, reps), reps))\n",
    "print('jump error:', error(w, empirical(jump, w, n, reps), reps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Is it faster?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from arsenal.timer import timers\n",
    "jobs = []\n",
    "methods = {'walk': walk, 'jump': jump} #, 'base': base}\n",
    "for i in range(20):\n",
    "    n = 2 ** i\n",
    "    for rep in range(3):\n",
    "        w = np.random.uniform(0, 10, size=n)\n",
    "        for name, method in methods.items():\n",
    "            jobs.append([name, n, rep, method, w])\n",
    "np.random.shuffle(jobs)  # shuffle jobs to avoid accidental correlations\n",
    "T = timers()\n",
    "for [name, n, rep, method, w] in jobs:\n",
    "    with T[name](n=n):\n",
    "        method(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "jump is 4.2050x faster than walk \u001b[0;33m(median: walk: 0.00163066 jump: 0.000387788)\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "T.compare()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Where is the speed-up coming from?** Both algorithms are asymptotically linear time to sample.  However, the constant factors associated with the `jump` algorithm are smaller because it samples substantially fewer random variates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5jkV3ng+++p3FWdcw4TNTlqpFEcSZYlQIEoJJCxEEYXbNl+2Ot9jO19Hu7d9a5YX3z3gQs2K4PA7Bowy+ViSSiABoVRnpE0mpxD51g5p3P/ONXTrWFCT3dXV3X3+3mefmbq11W/OpPqnXPe875Haa0RQgghZsJS6AEIIYSY/ySYCCGEmDEJJkIIIWZMgokQQogZk2AihBBixiSYCCGEmDFboQeQD7W1tbqzs7PQwxBCiHnlnXfeGdVa103ntQsqmCil7gbuXrZsGXv27Cn0cIQQYl5RSp2d7msX1DKX1voprfUjFRUVhR6KEEIsKgsqmAghhCgMCSZCCCFmbEHlTC4llUrR29tLPB4v9FDmhMvlorW1FbvdXuihCCEWgUUTTHp7eykrK6OzsxOlVKGHk1daa8bGxujt7aWrq6vQwxFCLAKLZpkrHo9TU1Oz4AMJgFKKmpqaRTMLE2Kx01rjiyTp9UXxRZIUoht80c9MlFIe4B+AJPCS1vpfZnCvWRtXsVtMv1YhFjOtNceGwvT7Y3R7o6xoKKPa42BFQ+mcfg4UZGailHpCKTWslDpw3vU7lVJHlVInlFJfzV3+OPBzrfUXgXvmfLAFtmPHjnM1M6WlpQUejRCi2PijKY4OBnnh8BDvnPVyejSMN5LEH03N6TgKtcz1Q+DOyReUUlbgO8CHgNXAA0qp1UAr0JN7WmYOxyiEEEVNa83rJ0d56dgIQ8E4sWSWcpfZdBNJpud0LAUJJlrrVwDveZe3ASe01qe01kngp8C9QC8moMA8zvH83d/9Hd/61rcA+MpXvsKtt94KwM6dO3nwwQf58pe/zNatW1mzZg1f+9rXLnmv0dFRtm/fzq9+9au8j1sIUZxiyQy/3NvHu2d9DAZiuB1WbltVT7XHAYDHMbdZjGLKmbQwMQMBE0SuAb4FfFsp9RHgqYu9WCn1CPAIQHt7+yXf6L/95thMx3pBX7l9xUW/d9NNN/H3f//3/Nmf/Rl79uwhkUiQSqV49dVXufHGG/nUpz5FdXU1mUyG2267jX379rF+/frfuc/Q0BD33HMPf/u3f8vtt9+el1+HEKK49ftjPLN/gKFgnNFwko1tVWxqq6LUZT7Sqz0OKt1zWxZQTMHkQpkirbWOAJ+/3Iu11o8rpQaAux0Ox5ZZH90MbdmyhXfeeYdQKITT6WTz5s3s2bOHXbt28a1vfYuf/exnPP7446TTaQYGBjh06NDvBJNUKsVtt93Gd77zHW6++eYC/UqEEIWitebdbj+7jo0wFIqjs3DX+iZuWlFHNJkhkkzjcdiodNvnfBNOMQWTXqBt0uNWoP9KbqC1fgp4auvWrV+81PMuNYPIF7vdTmdnJz/4wQ+47rrrWL9+PS+++CInT56kpKSEb3zjG+zevZuqqioeeuihC27rtdlsbNmyheeff16CiRCLTDyV4deHhjg6EKTHF6W1ys0nNreypM6DUgqn3UpVbomrEIopB7EbWK6U6lJKOYD7gSev5AZKqbuVUo8HAoG8DHCmbrrpJr7xjW9w0003ceONN/Ld736XjRs3EgwG8Xg8VFRUMDQ0xLPPPnvB1yuleOKJJzhy5Ahf//rX53j0QohCGQrG+fFb3ezr8dPrj7G1o4ov71jK0vq53f57KYXaGvwT4A1gpVKqVyn1Ba11GngUeB44DPxMa33wSu5b7F2Db7zxRgYGBti+fTsNDQ24XC5uvPFGNmzYwKZNm1izZg0PP/ww119//UXvYbVa+elPf8qLL77IP/zDP8zh6IUQc01rzfs9fn7ydjdHB4OEEinu3djMZ67poMxVXK2SVCEqJfNl0nkmXzx+/PgHvnf48GFWrVpVmIEVyGL8NQuxUCTSGXYeHuZAX+DcstYD21ppr85fvZlS6h2t9dbpvLaYlrlmrNhnJkIIMRUjoQQ/eaubPWe89PrMstaf3rosr4FkpoopAS+EEIuS1hp/NEUkmabXG+Pt02P0+GJYFHxqSyvXL6vFYimO3MjFLKhgMvnYXiGEmA/Ge2sNh+Ls7fZzYjhMIJbiqqYy/nB7J02VJYUe4pTIMpcQQhSQP5ri+FCInYeHOdgfYCSUoLPGzUPXzZ9AAgtsZiKEEPNJKpPlxaPDvHJslJFQHLvNwq1X1bG6uYJ0dn5tjlpQwUSWuYQQ80W/P8ZvDg1xZCDIYDBOW5WL21c3ntvyO9e9tWZKlrnm2HXXXVfoIQghCiiVyfLKsRF+/FY373b7yGrNXeub+Oim1nOBpBC9tWZqfoW+BeD1118v9BCEEAXS74/x/MFBuscijISTrGgo48Fr26lyO87t5ipUb62ZWlAzk2JvpwLmgKuXXnqJu+6669y1Rx99lB/+8IcAdHZ28td//dds376drVu38u6773LHHXewdOlSvvvd7wLw0ksvcdNNN/Gxj32M1atX86UvfYlsNluIX44QYgpSmSwv52Yje3v8RJNZPrWllUdvWUa1x4lSiiqPg9YqN1Uex7wLJLDAZiZTbfTIi4/lZwC3/NWs3KatrY033niDr3zlKzz00EO89tprxONx1qxZw5e+9CUA3n77bQ4dOkRHRwd33nknv/jFL/jkJz85K+8vhJg9ff4Yv87NRobDSVbmZiPVHmehhzarFlQwWSjuucecTrxu3TrC4TBlZWWUlZXhcrnw+/0AbNu2jSVLlgDwwAMP8Oqrr0owEaKIpDJZXj85xu7TXnp8UZw2K/dtaeXG5XVFX4A4HYszmMzSDGK6bDbbB5alzm8373Sa/7FYLJZzPx9/nE6bozjPnwbPx2mxEAvF5Ap2j8NGJJnmN4eGzs1Grmoo47MLcDYy2eIMJgXW0dHBoUOHSCQSxONxdu7cyQ033HBF93j77bc5ffo0HR0d/Ou//iuPPPJInkYrhLiU8Qp2byRJOpvlYF+Qs94IkUSaEruN+7a0csPyOqwLcDYy2YIKJvOhzkQpRVtbG/fddx/r169n+fLlbNq06Yrvs337dr761a+yf//+c8l4IcTc80dTeCNJxsIJdp/xMhiME4il2NxWyRduWkLNAp6NTLagWtCP27p1q96zZ88HrhVDO/axsTE2b97M2bNnZ3Sfl156iW984xs8/fTTl3xeMfyahVjour0RXjhkWsUPBePYrRa2L6nhznWNdNR4Cj28KzKTFvQLamZSzPr7+9mxYwd/8Rd/UeihCCFmSSie4sUjI7x71sdoOEFrVQm3r26kvMROeZEdXpVvEkzmSHNzM8eOHZuVe+3YsYMdO3bMyr2EENNzYjjMbw4N0u2NEktl2NJZxfYltVgtal5WsM+UBBMhhLgC6UyWXcdH2XPWHFxV7bbz1TtXUlvmmtcV7DO1qIKJ1nrR/AEvxFyYEIU2Fk7wzIFBzo5FGAklWNtcwWevacftNB+lVR5HgUdYOIsmmLhcLsbGxqipqVnwAUVrzdjYGC6Xq9BDEWJB0FpzoC/IS0eH6fVFAfjoxmZuXlFfHAWIWkPMB8kwOEqhpArm+HNuQQWTS20Nbm1tpbe3l5GRkbkfWAG4XC5aW1sLPQwh5r14KsMLh4c43B+kzx+jpaqEB6/toKmiSA6u0hqGD0N4GCIjUNYI7hqoXzWnAWVBBZNL9eay2+10dXUVYFRCiPmqzx/j2f0D9Ptj+GMptnRW8emt7ThsRdQjN+aD0ePQ8xZEx2D5HRPX3dVzNowFFUyEEGI2ZLOat894eePkGH2+KG6njc9e087m9qriWibPpOHEC3DyRRM8dAbSZhmOZFiCiRBCzKXJvbW0hjdPjXF6NEK/P8aKhlLu39ZBdbEl1/09cPQZ8J6CYD+UNUHXTVDZZr7vKJ3T4UgwEUIsapN7a/X5Y7xzxks0mcFpt3Dj8lo+uqm1uPpqpRNw6mXo3QOhAZMXWXIzNG0Ea662xV1jkvBzSIKJEGJR80dTDAfj7OsLcHI4zEgoQYXbzqe2tHJ1V02hh/dBYyfh2HMQGjRftVfBxk9DeYvs5hJCiEI6PRrmhcNDjEWSBGMpuuo87FhRR1NlkezWAkhG4eROGNgHgV5weGD1R2HlnWCxmue4q+c0R3K+og8mSqklwN8AFVprOf1JCDErMlnNW6fGePnYCGe9UZxWC9u6qs8l2T2OIvh41BpGjsDxX0NoyOzWalgHG+8HT22hR/cBef3dUko9AdwFDGut1066fifwTcAKfE9r/fWL3UNrfQr4glLq5/kcqxBi8RgNJ3j+4CC93ijDoThdtR6u7qymym2S7EXRWysRgmPPm2Di7zZ5kI2fha4b53wJayryHXp/CHwb+NH4BaWUFfgOcDvQC+xWSj2JCSznH87+sNZ6OM9jFEIsElpr3u328/qJUQaDMTSwY0U9H1nfRCieLlxvrckV7HYPBLrNdt/wEMRD0LIVNjwArrK5G9MVymsw0Vq/opTqPO/yNuBEbsaBUuqnwL1a68cwsxghhJh1gViKXx8c5MxohP5AjPYaN5/c3EZbtRswfbUK0ltrvII9OmZmIz1vQ8wL6bjZ7rvt09Cyee7HdYUKsSjYAvRMetwLXHOxJyulaoD/DGxSSv1VLuhc6HmPAI8AtLe3z95ohRDzmtaaQwNBXjo6wlAwTiKdYfuSGj62uRW7tQgq2WM+iIzC6DEY2AtRL2RSsPQW2PJ5cBTRRoBLKEQwudDc8aItbrXWY8CXLndTrfXjSqkB4G6Hw7FlBuMTQiwQ0WSanYeHOTpo+mo1Vrh48JoOVjQW0XJR1AtnXzPFh5ERKG2Arpuh64Z5E0igMMGkF2ib9LgV6J+NG1+qN5cQYnE5ORLmhUNDDAbjhOIptnZU8fHNrZQUwy6tceFhOPykSbInwtCwFpb9ntnuO8cV7DNViN/V3cBypVQX0AfcD3xmNm58qa7BQoiFa3I7FLvFwvu9fvb3Bej1RanxOHnwmg42ts9tRfhlDewzBYj+HrDYYPnvmWACBalgn6l8bw3+CbADqFVK9QJf01p/Xyn1KPA8ZgfXE1rrg7PxfjIzEWLxmdwOZSSUYPcZL5FkGjSsa6ng01e34XEW0RG6mZSpG+l7FwI9UNVlciM2R0Er2Gcq37u5HrjI9WeAZ2b7/WRmIsTi44+mGA0lONAf4OhQiNFQghK7lY9tbuGONY3F1eU36oWDvwDvGYiOQts1sPlzk3pqFa6CfaaKaPFw5mRmIsTi0+uPsvPIMCOhOL5oitaqEm67qp5VzeXFFUiGj8CRX5kCRIsF1t0Hy24t9KhmzYIKJkKIxUNrzfu9AX59YIAzY2EUiqs7q9jaWY2lWNqhAGQzpgCx+00TSCpaYOsXoHJhnYRaJL/bs0OWuYRYHCKJNL85NMTxoRD9gThN5SVcs6SaujIXUCTtUADiATj4Sxg7YarZWzbDlofA7ir0yGad0vqiJR7z1tatW/WePXsKPQwhRB6Mb/kdCMSJJtNsaa/iY5tbiCQyhWuHciFjJ+HQk+A/a2YnV30EVtxR1Il1pdQ7Wuut03ntgpqZCCEWrmQ6y67jI+zt8dPri1Jb6uSh6zpZ3VwBgMNmLUw7FDivt5Ybhg7A6VdNICltgKu/ANVdhRnbHFlQwUSWuYRYmIaCcZ7dP0CfP8ZYJMn6lgo+fXU7bmcRfIRN7q2VisHZ1yEybE5EbFoPV38RHO5CjzLviuBPYvbIbi4hFpZsVrPnrI83To7S74/htFv52MZmblheV/hlrHExnwkk4SE485o5SjeThLUfgw2fNTu3FoEFFUyEEAtHIJbi+YODnB6NMBiIsayulPuubjuXZC8aMR/0vQNDh0xAcZbBqrvNuSOLJJDAAgsmsswlxPw0uR2Kx2FjIBDjxaMjDAbipDIZblxWx72bWrBaimQ2Mi7QC4f+Dfreg7jf5EVWfsQsa82z3loztaCCiSxzCTH/TG6Hkkxnea/Hx1AwTjyVoamihE9u6WRJXZF9MGfScGaX6fbrO2t6oS+7HZo3mt1a87C31kwtqGAihJh//NFUrq9WnN1nfIyFE0SSaa5fWsvnruvAZS+yj6nQIBx+ymz9DQ9D3SrY9jCg5nVvrZkqsj8lIcRiE4gn2dfrP9dXy+O0sWNFHbevaSyuQJLNmJ1ap18B3xmwOWH9/abb73huZB731pqpIvqTEkIsNiOhBM8fGGJ/X4BANEVbjZtbVtRR6rIXTzsUgPAIHHkKRk9AsB/qVsDVfwSl9YUeWdEooj+tmZMEvBDzQzarea/Hx6vHR+n1RdEarllSzeb2KpRSxdMOJZuFnrfg1EvgOw3KCus+CSs/ZA6wEucsqGAiCXghil8gluLXBwc5MRxmKBinq9bDA9vasFutxdUOJeo1uZHRYxDog+olZjZS0VzYcRWpBRVMhBDFS2vN4YEQvz0yRI83SlbDjpV13LOhBUtuy29xtEPxgPcknNhp2qFobepGVt8LVvnIvBj5nRFC5F0smWHnkSEO9AboD8RorCjhD67toK26CNqMTG6HkghDz5um+DAVh8p22PbIgmsXnw8STIQQeXVmNMJzBwc5MxomnspyzZJq7tvSht1WJDmHmA8io6ZNfN8es/U3mzazkc1/aI7TFZc1pWCilKoHrgeagRhwANijtc7mcWxCiHkslTFdft8+5aXPH6PKbefh65eworGs0EP7oGAfnNxpig8jI+CpNVXsHddKILkClwwmSqlbgK8C1cB7wDDgAj4KLFVK/Rz4e611MN8DnQrZzSVEcRgMxHlm/wDHh0MEY2k2tlXw4LUdlBTTdt9s1vTUOvI0DOwDnYbmTbD0VnMm+yJrhzJTlzwcSyn1fwH/j9a6+wLfswF3AVat9f+bvyFeOTkcS4i5Nd5bK5RIcXQwxLtnffT4YngcNh7Y1sqGtqrC786aLDJqzmMfOWJ2ajnc0LXDHKkLph1K/apFV8Wet8OxtNb//hLfSwO/nM6bCiEWjvHeWmfGIrx9ysvpsTCJVJaru6p4+PouykuKaKkomzUJ9pO/NctaOgtXfRhWfxRS0UXdDmWmppoz+XPgB0AI+B6wCfiq1vrXeRybEGIeGAkl2HV8hH29AYaDcaxWxc0r6vjEltbiCiShodxs5CiE+k2H361fgKoO8327a1G3Q5mpqS5gPqy1/qZS6g6gDvg8JrhIMBFiETs9GuHn7/RwsC9IOJGmpbKED61rpKLEQSyVKfTwjGzGdPc99ZKZjaDMTGTVXaa/lpgVUw0m4/O9DwM/0Fq/r4pqAVQIMZeC8RQvHx1mb7efs2NRslpzx5pGVjeVnytALIreWsF+OPy0qWIPDULtctj6MFS2FXpkC85U/7TfUUr9GugC/kopVQbItmAhFplMVvNet49dx0fp8UaJJjNc3VXF5vZq0tmJzTwF762VSZnzRk6/Ar5u00dr/adhxR2y3TdPphpMvgBsBE5praNKqRrMUpcQYpHo9UXZeXiYUyNhBoNxGstdfOGGLlY1laMUHzgpcc57a01uhxIPwOldpggxPAx1K01uRHpq5dXl6kw2n3dpSSFWt5RSHwU+AtQD35HEvxBzJ5pMs+v4KO91++j3xwC4fVUDH1rXhMc58RFS5XEUprfWeDuU0CAMvA+D+80Rup5a2PSgqRuR2UjeXW5m8ve5H13AFmAfJn+yHngLuOFyb6CUegJTjzKstV476fqdwDcBK/A9rfXXL3YPrfUvgV8qpaqAbyCJfyHyLpvVHOgP8MqxEQYCccbCCZbXl/LxLa101niKp24k5jM7tM68as5kjwegqgtu+DNTKyLmxOXqTG4BUEr9FHhEa70/93gt8BdTfI8fAt8GfjR+QSllBb4D3A70AruVUk9iAstj573+Ya31cO7n/yH3OiHELBovOhxfpkqkM7x4dIRTI2H6/TFKnTY+taWN65bVUuIokp5aYOpGTr8CR5+B0ABY7LDiQ9C8HhyeQo9uUZlqzuSq8UACoLU+oJTaOJUXaq1fUUp1nnd5G3BCa30KzgWre7XWj2FmMR+Q2zn2deBZrfW7UxyzEGIKxosOvZEkyXSWg/0Ben0x0tkskUSa9W2V3L2uiabKkuKZjQDEg6YVSs9uCPRAeStcdRd4asz3pR3KnJpqMDmslPoe8D8BDTwIHJ7B+7YAPZMe9wLXXOL5fwr8HlChlFqmtf7u+U9QSj0CPALQ3t4+g6EJsbj4oym8kSR9vhjvdvvwRRJ4IymW1nv43PZONrRV4rIX0WwEzLLW4afBe8pUri/7fWjeaHpqgWmHUlJV2DEuMlMNJp8Hvgz8ee7xK8A/zuB9L/Tfm4s2CdNafwv41qVuqLV+XCk1ANztcDi2zGBsQiwqwXiK93v8HBkM4o0k0RrWtZbzkfVNbGovsorwTMocWtX9pjm4yl0D1/4J1F9lku7SDqVgphRMtNZx4L/lvmZDLzC5aqgV6J/pTeXYXiGuTCSR5sUjwxzsDzASStBQ7mRrRzUdtR46a4psmSg0BId+CWMnzeFVLVtgy+fBlWtp766WdigFNNXeXNcD/wfQMfk1Wusl03zf3cBypVQX0AfcD3xmmvc6R1rQCzF1ff4Yv9rXT78vRjCeYmVjGds6q6kudRa+6HAyraF3t5mReE+Zpawtn4eum8BiKfToRM5Ul7m+D3wFeAe4ooY7SqmfADuAWqVUL/A1rfX3lVKPAs9jdnA9obU+eCX3vRCZmQhxeVpr3uvx8/LRYQYDcbIaPr21jQ1tlWgoTNHhxSTCpjnj0EGTZK9eCtu+CGUNhR6ZOM9Ug0lAa/3sdN5Aa/3ARa4/AzwznXtejMxMhLi0ZDrLC4eHONAXoMcbpaHcxUc3NbOqqQKrpQiCx2RjJ02S3XfG5EKW3gYbPj2RZBdF5ZKHY517klJfx8wgfgEkxq8X6zZdORxLiN/ljSR5el8/3WNRBgNx1rVW8NGNzTRUlBR6aB9sh2JzmUr27jfB322S6Zv/EBrXFHqUC17eDseaZHzb7uQ30cCt03lTIcTcOjEc4rkDgwwE4kSTaW5bVc8daxtxF0Nn3/F2KNExiPnh7OsmsOgUNG40XX6dUoBY7Ka6m+uWfA9kNsgylxAflM1qXjs5ypunxujzxShz2fnc9k7Wt1YWz7JWzGcaMo4eg4G9ZteWxQobHoA1H5UtvvPEVHdzVQBfA27KXXoZ+I9a60C+BjYdkoAXYkIkkeaZ/QOcHA7T64+xoqGM+7a00FjpLvTQjGzWFByeehmOPWdmJukElDfDig9D6yYJJPPIVOe4TwAHgPtyj/8Ac9Lix/MxKCHEzPT7Yzyd2/bri6W4cVktd29sLsyyVjZjgkYqZnZnxbxmOct3CgYPmG6/kRGTG+m4AdqvMTMTaYcyr0z1b9ZSrfUnJj3+P5VSe/MxoJmQZS6xWI03agwnUpwZjbLnrJc+XwybVXH/1W1s66qZm2WtZBQSQUiETNCI+iEZmvi+skA6OXEOezoOzjJo3w6Na03yHaQdyjw01WASU0rdoLV+Fc4VMcbyN6zpkWUusRiNN2ocDsZ5t9vPqZEw/liStc0VfOaadlqq5mhZKx4w56xnUiZoWJ3mjHVPvVmuSsfNTGTsuJml6Cx0bIe1nzTBY3w3l7RDmZemGky+DPxzLncC4AMeysuIhBBXxB9NcXggyHvdfoZDcfyxFCvqS/nsNR00V83Rtt9kxJwnYrFDyXktTbIZMxMZOmiel0lAwxrY8Bmo6ph4nrRDmdemuptrL7BBKVWeexzM66imSZa5xGITSaT51f5+dp/2MRpOYLdauGlZLetbK8levHfq7ErFTSDRWXBVTlzX2jRjHHjfzFrSMShvg433Q+N6mXksMFPdzfVfgL/TWvtzj6uA/11r/R/yObgrJctcYrHIZjX7+wLsOj7C2bEIQ8E47TUebllRR2Xu6FzPXCTbMynofsMsW1kdZouvvQRSEejfC5FRk3wvqYJNfwCd15vkulhwpvq37UNa678ef6C19imlPow5+VAIMYeGg3FeODzEqdEIA/4YVW4Hn9jSQlOF+1ySfU4aNWYz0LvHbOlNxSERMMtYYyfMwVVWh8mZrPk4rPwQ2F35HY8oqKkGE6tSyqm1TgAopUoAZ/6GJYQ4XyKd4c1TXt46NcZgIE4yk2VbZzX3bmqhosT+gWN3896oUWuzfBXoBbvbHJk7dtK0P0nnOi6tvBM2PSi7shaJqQaT/wnsVEr9ANNG5WHgn/M2KiHEOVprTo6E2Xl4mB5vlJFwgo5qN5/Y0sby+lIsudlIlcdBVW6JK++GDsHocShrMo0Ye9+BuM8EkprlsPRWaN0igWQRmWoC/u+UUvswR+cq4D9prZ/P68imQRLwYqEJxFK8eGSI/b0BBgJxHDYLd21o4taVDXicBeqrNXYShg9AaYOZmZx9DaKjYLHB6o+bUw+VkqLDReZK/jYeBtJa6xeUUm6lVJnWOnTZV80hScCLhSKT1bzb7ePlYyP0+2OE42nWtZRz78YW2qrdhTtrxN9j8iTuGuh/F0aOmQR7RTusuGNia68UHS46U93N9UXgEaAaWAq0AN8Fbsvf0IRYPMYr2CPJNMFYijdPeTkxHGYkHKe21MnHt3ewtbMal72AO6HCI6YtvM0JJ39rZiOpGCy7Ddbfb4KKFB0uWlOdmfwJsA14C0BrfVwpVZ+3UQmxiIxXsA8EYuzvDXCoP0A4mcHjsHL9slo+tKax8GeOxHxwZpfZpTW417RNsdhg25eg8zoTOOwuKTpcxKYaTBJa6+T41FopZYO5qogSYmEbr2B//eQoA/44kWSapgoXn766lWuX1OGwFfic80TIzESGj5ijc1MxqGiFa/8EKpoLOzZRNKYaTF5WSv01UKKUuh34Y+Cp/A1LiMUhncmy88gQLxwawhtJYrda2L6khvVtlSypKyt8IEnFzBns3W+an6fj0HWjaYXiKJJW9qIoTDWYfBX4ArAf+N8wZ7d/L1+Dmi7ZzSXmk7Fwgiff72dvt5/hUIKWyhKuW1pDW7UHq0bmHQkAACAASURBVEXNTQX7paSTsPfHcPYN89hihav/CDquB0uBg5woOlM6A/4DL1CqGmjVWu/Lz5BmTs6AF8VMa82+3gDPHhigxxslk9WsairnqqZyyl2mar3a42BFQ2nhdm0lo/DGt02XX2Uxy1nbvgRV7YUZj5gTeT8DXin1EnBP7vl7gRGl1Mta6383nTcVYrGKJtM8d2CQPWe8DAUTNFW4+NjmFtY2VxBNZuaugv18Wk+0gE9G4d0fQaAblBU6rjNH6LrK5m48Yt6Z6jy6QmsdVEr9EfADrfXXckWMQogpOjMa4cn3+zk5HCKSzLB9aQ33bGimvtz0rHLarXNXwT6Z1jB82DRlHDsJJ35jDrgqb4H1n4aum8Ba4CU3UfSm+jfEppRqwhzb+zd5HI8QC046k2XX8VFePDrMQCBOudPG57a3sq2rprB1I+NiPvCeNkWIQwfNEbrlrbDtEdMSRYgpmGow+Y/A88CrWuvdSqklwPH8DUuIhWE0nODJvX0c7A/ii5rTDz+5pbWwVeyTJaNw9Bk4sdMEkXQCmjeZQsTSukKPTswjU+3N9b+A/zXp8SngExd/hRCLm9aavT1+ntlvkuxWi4VPbm5lx1X1uAu9SwvMOSS9u+HUy6ZhY6AHSpvMTq26FWC1S28tcUUu+bdaKfUfgH/QWnsv8v1bAbfW+ul8DE6I+WJyOxQFvH5ylD1n/IyE4nTWenhgWztL6yY6/BZMNgtDB+Dki6bbb3jIVK2v/RRUdU7UjkhvLXGFLvdfpP3AU0qpOPAuMAK4gOXARuAF4L/kc4BKqVXAnwO1wE6t9T/m8/2EuFLj7VC8kSQDgRivHh9lNJzAabNw++oG7trQTJkrzwdVTYX3FJz4rTm8KtgPNges/LA5uMpdPbGbS3priWm4ZDDRWv8b8G9KqeXA9UATEMScb/KI1jp2qdcrpZ4A7gKGtdZrJ12/E/gmYAW+p7X++iXGcBj4klLKAvzTlH5VQswhfzTFUDDO+z1+3u/x440mqXDZ+dx1HVy3tO7c6YcFEx427VCGDkGwz5yQ2HY1rP4YVLRMBA13tfTWEtM21ZzJcaaXcP8h8G3gR+MXlFJW4DvA7UAvsFsp9SQmsDx23usf1loPK6XuwVThf3saYxAir44MBnl6Xz+93hjxdIaVDWXcsrKertrSwgaSeBBOvwL975nlrHgA6lebY3Rrlsp2XzGr8vq3SWv9ilKq87zL24ATuSQ+SqmfAvdqrR/DzGIudJ8ngSeVUr8Cfpy/EQsxdbFkhhePDvPikWFOjYQpcVi5c20jq5sq5r4dyuSiQ4sdRo5C79tmVhIegsoO2PwQNK0zLeSFmGWF+K9JC9Az6XEvcM3FnqyU2gF8HHPm/DOXeN4jmDNXaG+Xlg8if7TWHB8O8+yBAU6PRAjG06xtqeCaJTVUuU3RYbXHQaV7jvIk40WH4WGTDxncD5mEaczoqYWNn4XOG8Apu7NE/hQimFxo3n/RBmFa65eAly53U63140qpAeBuh8MhlVYiL8KJNDsPDfHOWS8DwQTlLhufv66Tq7uqiCWzhWmHEvOZIHL2DYgMm3oRqwOW3Q7rPyV5EDEnptqbawXwj0CD1nqtUmo9cI/W+m+n8Z69QNukx61A/zTu8zvk2F6RL1prDvYH+c2hIbq9UULxFBvbKvn4phaaq8x2WpeduW+Hks3AqRfh+AsQ80MmDrUrTb1IyyYJJGLOTHVm8k/Avwf+O4DWep9S6sfAdILJbmC5UqoL6APuBz4zjfv8DmlBL/IhEE3xm0OD7O8LmHYoJXY+e20713bVUuIoYDuUyCgcfgpGjkCg11Ssd95h6kUsVnBKY0Yxd6YaTNxa67fPm7anL/cipdRPgB1ArVKqF/ia1vr7SqlHMe1ZrMATWuuDVzbsC5OZiZhN2axmb6+fl44M0+uPEY6n2dBawV3rm2iv8RSuHUo2m6tefwlCg5AKQ/u10LJVig5FwUw1mIwqpZaSy20opT4JDFzuRVrrBy5y/RkukUyfLpmZiNkyGk7wm0NDHBsM0R+IUVni4L6tbVy7tIZSZwG31Ea95uTDsRPg74HyJtj8h2aXVjwgRYeiYKZ0OFausePjwHWADzgNPKi1PpPX0U2THI4lrtR4O5RgPMWxoRDv9/gZCMSJJjOsaynn91c30lVXwLoRrU1X3xM7IThgWsS3X2cS7K7ywoxJLDh5PxwrVxPye0opD2DRWoem82ZCFKPxdignhsPsOeNlIBAjnEjTXu3m91c3cO2SWirmapvvhcQDcPRZGD5kZiOl9bDpUWjeIMfniqIx1d1clcDngE7M2SYAaK3/LG8jmwZZ5hLT4Ysk2X3Gy7vdPsbCCRLpLF21Hu7e0MzVndXYrAX6wNba1Iwce960QUkETW5kw/3gqijMmIS4iKku/j4DvIlp/JjN33BmRhLw4kol01mePTDIm6fGGA0n8DisbOyqYVVTOW3V7sIFkkQYjj0HA/tMe3h3NVz7J+asEZmNiCI01WDikvPexULjjSR5+v1+DvYHGA7Gaa4sYXNHFZ01HuxWS+HaoQT74exr4Os2S1xt22DjZ6Ckcu7GI8QVmuq/lv+hlPoi8DSQGL94sXNOCkWWucRUHRsK8ez+AXp9MaLJNNd0VbO8oYzaUidKqcK0Qwn0QN+7MHrMBJHKdjMbad1s6kaEKGJT3c31J8B/BvxMtD7RWusleRzbtMluLnExmaxm1/ER3jg5Rq8vSqnLxh1rGrm6s5p0Rs9tO5RsFkID0L8XzuwyNSPJCKQiUHsVXP+nUN2V3zEIMUned3MB/w5YprUenc6bCFEMQvEUv9o3wOGBIAOBOMsbSrlnQzPL6svObfnNezuUZAS8p02dyMhRUzcS6jf9tLJZU2y46l6oWymzETGvTDWYHASi+RzIbJBlLnEx3WNRntrXz5nRCOFEmptX1PGhdU1Uz1bwmJzzmFw0OD778J6EkePgOw2JkNmZlU2B3QPuWnDXQf0q0wrF7jL3lDPYxTwy1WCSAfYqpV7kgzmTotoaLLu5xPm01rx92suLR4bp9kVxO6w8eG07Wzurcdpm6X/+4zmP6Jh5nI6b9u86a/pmRb0mB5KKgrJAaQPUXQ2uSjP7qOo019PxiXtKOxQxz0w1mPwy9yXEvBFPZXjuwADvdfvp98dY3lDG/Ve30Vbtnt18SMxnAon3tCksDA2awGGxmSDhKDWHU1V3mUCSipm8SHkTNKw1QeNiMxsh5ompVsD/c74HIsRsGgrG+eXePo4Phggl0ty6qoG7NzTnp69WPAA9b8PA3onZSUkV1K+BpvXgLDeBIRmB6Ch46kzxoadm4h5KyRnsYl675L8spdTPtNb3KaX2c4EDrLTW6/M2MiGuwHhvrXAiRbc3xqvHRzjrjVJit/LIjUvY2F6Vn75aiTAcfhrOvm5mIw1robLVzEqqukwb+FTMnDXiroKum007FJl1iAXmcv9N+/Pcjxc8m73YSAJ+cRrvrTUcivPuWR8H+4OE4ik2tFby+Rs6aSgvyc8bB3rhnX82dSFWO7TcMFFY6KwAq80seTlKofN6KGuS6nWxYF0ymGitx9vM/7HW+i8nf08p9V+Bv/zdVxWOJOAXJ380xdHBIG+d9nJ2LEoinWFjeyWfv76LhnJXft60Zze89z/MslZlGyy5xZx6mIqZpHo2BamEqV6vaJVtvmLBm+oC8u38buD40AWuCTGnYskMv9rfzyvHRvFFkpQ4rNyzsZmltaWksnloI5fNwL6fmeaLOgvNm6FlSy6A5IKJUuZ6ZYeZnQixCFwuZ/Jl4I+BJUqpfZO+VQa8ls+BCXEp42ey//rgIEcGg4yFk1zVVM4tK+vx5JLss95bKzoGb/6j2QZs90DHdVCz1HwvETIJ9vrVULscbHN8FrwQBXa5f20/Bp4FHgO+Oul6qNj6conFwxtJ8tyBAd7vCTAaTtBQ7uTu9S1Ulzqw5BLbs95ba+B9eOu/m4R7aQN03gieWsikIDYGrmpYvl1qQ8SidbmcSQAIABc8fleIuZTOZHnrlJffHhmizx9DA7etauCeDU14nDb80dTs99ZKJ+HgL+HIU2B1QFWHCST2ElMXkklB40aoWSJ5EbGoLagFXdnNtXB1j0X51f5+jg6G8MdSLK3zcN+WNpY3lGGZ1FdrVntrhYZg9z+ZHlp2t+mX1bwZsmlzdG55ozlfxFk2e+8pxDw1pa7B8410DV44osk0Ow8P8fqJMUYjSZw2xe2rG7ltVUN+ChDBzDb63oH3fmwq1W1OaN0G1Usgllvdbd5kWsRLvYhYQOaia7AQc0przYG+AM/sH6DbGyWRzrK6qZx7N7bQVes5NxuZxTc0PbQCPaaavXe3aYXiqjDLWg6PadhY3QWN68wylxDiHAkmoiiMV7BHkmkSqSy7ToywrydAKJGi2u3gng0NXL+8Nj+zkVQczrxqzlv3ngR/t9mtVbsMOq43s5NMEpbcDGWNs//+QiwAEkxEwU2uYD/YH+CdMz6iqQyVJXa2tFfx+2sa6ayZ5dnIeGPFsVPQ/QaMHIJk1LSGt1ihtBYa1pjHdVeZL9nuK8RFSTARBeePpjjYH2DXsREGAnHSWU19mZM71zSy46p6ylyzuMU3GTMJ9e7XwXfG9MxKRc3sA8zSVs1y0wIlm4aVH5bmi0JMgQQTUVCBWJJ/eessu46PorXG7bDRUeNmfWslG9orpx5ILtXCPR6Agf3QtweGDppzQ6xOU51uyeVFXOWmu6/dBcpqciPLb5dAIsQUSTARBRGMp3j7lJdnDwwQiKawWxS1pU6uaixneWMZboeNMucVBJLzD6dKJ0GnYfCgyYNkEqBsZmdWSZUJEs4KcFWB0w3BftMexVUB5S1mp5anLn+/AUIsMPMimCilPMArwNe01k8XejxierTW+KIpDg8E2HVslB5fDKfNQm2Zg9pSJ0vrSqkvd2FRauoV7Fqb5aq+dyAyCuFhCA+a9ibKamYfNpfp2OuuMQHCVWF+7q41QcXhMV+JkBxOJcQ05TWYKKWewLSvH9Zar510/U7gm4AV+J7W+uuXudVfAj/L20BFXmWzmtFwgqODIU6MhDkyECKVyWK3WmiscHHd0ho2tlWSyuiLV7BnM2ZHVSJsWr8HeiDYZ34eHIC4HzJpk+fQGRMcqpeYXlk1S00LFFeFuW53X7gVvBxOJcS05Xtm8kPg28CPxi8opazAdzCdiHuB3UqpJzGB5bHzXv8wsB44BOSpl7jIl3Qmy2AwzpHBIMFomh5flG5vDK01dquF9S3l3LSyntaqEhRAdIyq0Emz3GS1mWR5ZAgC/RAZNsnyRIjfOactmzKJ8/ImKG8zwcNVZqrVK5oL8CsXYvHJazDRWr+ilOo87/I24ITW+hSAUuqnwL1a68e4wCFcSqlbAA+wGogppZ7RWueht7iYLYl0hj5fjGNDIRLpLOmMZn9fgFA8RSaraaxwsWNlPRtaKylxWM1SVc9b5rRCf7dp457NmC26OmsCxfgsxVlqEuW2EpM0L22A2hUmT5JNmRMOwSxjlTcV7jdBiEWmEDmTFqBn0uNe4JqLPVlr/TcASqmHgNGLBRKl1CPAIwDt7e2zNVZxBbJZzYmRMMcGQ2TRlDrsnB2LcmwojNYam1WxpaOSW1Y20FjhMstYmRSceQ32/sTMPianKVxVJnfhqQVHGdidpg+Wo8w0XCxrMicbWqyX3s0lhMi7QgSTC/0Lv2yDMK31Dy/z/ceVUgPA3Q6HY8s0xyamKZPV7Ov1c9YbpcbjIBhLs+v4CMF4inRW01Th4vdWNbC2pQKXPffh7zsL+39utuyOHyrlrDBnpTvLoaLN/FxnzePKTijL5T7ODxRKSc5DiAIqRDDpBdomPW4F+mfjxnJsb2GkMlneOevj1EgYl93Cu2d9dHujaA1Wi+LaJVXcdlUD9eNH6Mb8cOAXcPY1QJklK5vLJMxtLrONV2uz86pxvalGl868QhS1QgST3cBypVQX0AfcD3xmNm4sLejnXjyVYfcZL4f6g0STaY4OhogmM1gsipUNZdyxpoG1LZU4bBZT+3H0WTj6K0gnTBDx1EDDWoj6ITpqigmrusz5IM2bZKlKiHki31uDfwLsAGqVUr2YOpHvK6UeBZ7H7OB6Qmt9cDbeT2YmcyuWzPDWqTEGg3F6fVG6vVGyWY1SiuUNpdy3tY2l9aUmmX7yJTj4/5nOvPYSUxjYuM5s080koWmtyYFY7WYWIjkPIeaVfO/muuAJjVrrZ4BnZvv9FvXMJJ00CWzfWbNUVFJtdj453GYGcKG6ihkIJ9K8eXIUfyzF3h4//f4Y6UyWmlIn1y6pob3GjdOqoWc3HPqlGZfVARWt0LjWBJFs2uRA6lbKcbdCzHNyONZ8ls2YFiK+M6aIT2vzIa2zZhlpfF+DUhM7o9zVZreTvcR8TeN//4FoildPDNPni3NiJEI0kcIfS7GhFq5ptlPiKcMRHmB54FXc/uOmiLCsERrWmQCns2Ypq3aZSaYLIYqCHI61mIxvgQ30gvc0ZJNgdUFJzcXPIB8PLoEe06dqnLJOBJiSSrNNd7w/lfvC9xsLJ/jtkWGODobwRZOkMpraMif3NgdpsAZwewcoPb4PT2qMEqcdSuuhaYMJImCS7DXLJh4LIRaEBRVMFvQyVyJk2oaMnYBkxBTnlVROFOldirJMzEQmy2ZMkBk7CSPHIO41ASWbMQ0R3dW5yU0WdJZQPMWxgQiWiGI1WSzZNJ3lWVY6MrgGvaTSmkw2iy0ZwuZyo9pvNOekKyvUrjSdeB3ufPzuCCEKTJa5ilkqBqEh8J6C2BhgNVXfNucFnqxN36pUzAQNZymgzEwmkzRndiRz53Yko7nHEfNjzJ+7/3lclSbPgSYQS3N4NEU4ZQE0JdYsyyuhxm3Dms7dS2HyNVVLTBK9vMkceVvVaVq7CyGKmixzLSSZlOl+6z1lzhzXmA/m0ksdF6th7DT4z5iTAdMJQJmgkoqZNiOXfM+ECRo2p9maa7GaGY+7Du2p4ag3zc5BTdrhIGm3sqLaypauMqrKPKaHVioB/tNmE0AmaV5b1QlLb5FjboVYJBZUMJnXy1yJsNnxNHrMLDPZS8BTf+kEeSZlzuEYOw7Dh8zuqMlcVeaoWYvdLC/ZPbkf3RPdc+1ukyAP9PzO7VPlnfzytIXX+5LYHVZKbWnu7LKzsaUUh9ITsxutTRBylE/kX0obzJcQYlFYUMFk3tWZaG3qLsaOmw/z8YT4pfIg6bgJIP6e3MwlYz7Qs+lcfUa1CUQ2l2lHUtVhgskld21piAchETh3ZShbwb/sg55gEqtSrCxPc/cSRZ0zhopFzVjLGs056SVVZodY3C+9sYRYpBZUMJk3MmkTCEaOmJ1ZthLwNFz8wzcZNcEm0GtqSSb3uvTUQd1VJpicn2Ava8zlPC5HoSvbiQR9xGIRDoY87OyzEo7FqbAlub05ydZWN86yOlMn4q42O77O3+0lvbGEWLQWVDAp+mWuZNTMKEaOmCUqZ5mp+gZMAj00kUBHm+Dh7zFtRsYpi3lNRRtUtJhlKjT4uj8ws8BZMeXttxpNjy/GgD/LG31WTvnCpDKatQ1O7t60grb2NlRJlalOF0KIC5DdXHMh5oOxU+A7bYKBq/K8D+ZcAj3YCzGvOXo2k4aSXEGfspqdURVtpg3JlezmuhitzZJZMsJYOM5zZzK873cxknKSsrjY0FHPA9ctp7VatvIKsVjIbq5ilM2YoDByxOzOsjrMkpSymA/yRMhcj46ZIOI788HlK4tt4uCn8Z5Vl6TMTOdi3XV11gSaVPTc+2Sclbwda+XJ0ylO+zPEMorGChe/v7Ke9hr3JWOREEJMJsFktown02NjEA+ZJHk6YXZNuashMgZDB03wiI7mtu/mjO+IcpSaXISn3pzjUdZsKsinI5uZCB5oQJmq9uolZF2VvDuseebQKP3+GJGElRKXjRs7qljdnDtvBPA45K+HEGJq5NPiSowXAKbjJhik42ZpKR4wgSI8nPsQj5jdVI4SE2Digd+9l80F7lrzAW93mfqQ83dxnZ9Qv5RMygSOdBzQ5l6eulwTxUpwlpPCwoG+AP/2Rj+DgRhZDR6nldtXN9BR7SGVnVjyrPY4qHRLjkQIMTULKpjMSgI+nasWTycgHTPBIhE0y1LJyMRSVDadqyCPmQryQLd5TTJituvCRJ2Hspgtu55c8HDXmhnLud1bV5hA19q8VyqSK0hUJo9S2mhmMq4KczJhrlNwJJHmUHeA5w4O0ueLYbEo7FYLmzuq+Mi6JurLXWit8UdTRJJpPA4blW67OVZXCCGmQBLw48YbKB5/wQQTpyf3oZ2cmIVMDi6p6MRrkxFTXzHOXmKCQXUX1CzP1Y5cpAnjxAAunkDPZsx7p2ITwcxVYXIpntwphHb3B7YWa63xRpIcGQzy6vFRur0xrBZIZzVXNZRxz8ZmOmo8EjCEEOdIAn6mtIahQ9D7NvTuzi0VWcBmN9+7EGWZSHhbbCbA2HNV5rZcbUdV1xUcNzspgZ5JQSxgltTQJhC5aycCk7Ns4j3Ok8pkGQrGOToQ4kB/gDNjUdCajM7SWuHh3o3NrGwsx2qRICKEmD0STMDMSGJeGNxvajtQZolo/AwQV5lZNnKW5z7wy01bEjV+4NQ06zyymYn8y+T+WbYSsxV4fMnKUXbZw60ymSyHB4Ls7fExEk4yHEqQSGWJpzLUeBzcd3UbWzqqzfG5QggxyySYwMQSVe0K8wHv8JjuvFVdUN48hRsoqGq/+DJVJmWaKaYTpn7Ekuvma7Gb5Hh5i/nR4blwq/jLGAnFeXb/AMeHwoxGEoTjaZRS1JU5uWdDEztW1lMiO7OEEHkknzBgtuSC6TOVipoZx/hZ5BeidS53kftx/Esps0srnfhgTsXmNMGissPMNOwlJnBcsPhw6kLxFAf7g+w6NsyxoTCJdJZUxuRU1rVU8Nlr2+mo8czoPYQQYioWVDCZ9m6ukiqzyyo6BiiT/3BW5nZ1xXIHROVqNciaH5XVBByL1cwwrE6wWsHiMAHDVT7Rlfci+Y3pSqQzHOgN8uLRIfr8cWLJNJFEGg0sqSvluqU1VHuckhcRQswZ2c01bnw3V2jAVKu7a8wMxWI1gePczy25r7n/oM5kNe/3+HJbfOM4bRYyWY1GU+V2cHVnNdUex7kdWlc1llHlmd1AJoRYuGQ312xQqmi73mYyWfac9fHcgUH6AzFK7FYcNgtlLhsrGsu4aVkN4WQWf3QiiS9Fh0KIuSTBpIjFUxnePu3lhcNDDAXNTKSixE5tqZOuWjdbO2toqyrBZrVI0aEQoqAkmBSJycEglc5ybCjEayfGGA7FsVoU9WUuGsqdNFWWsKapnK46D07bRCGkUooqj0OWtYQQBSHBpAhorTk2FObYUIjD/UFOjoSJJtN4XDY6azw0VriodNtZUlfKsvpS3LLNVwhRZORTqcC01hwfCvHzPb0cHwkTT2bQQFtVCRtaKygtsdNa5WZlYxnlLsmBCCGKkwSTOaa1JpxIE4ilGAzE2XPGy+4zXkLxNDaLhSV1HhrKnSTTmlKXnZtX1FMtS1dCiCJX9MFEKbUD+E/AQeCnWuuXCjqgK5TNasLJNIFoiqFQnOFggnQ2SyCa4uhgiEgyjdtho9Rpp7PGTSarsVksLGv2sH1JjQQSIcS8kNdgopR6ArgLGNZar510/U7gm4AV+J7W+uuXuI0GwoAL6M3jcGdFNqsJJdIEY0mGggmGQgkyuap0p92K02bhzGCEkyMRQONx2FndXIY/miSdhSX1bmo9TmrLnJJMF0LMG/memfwQ+Dbwo/ELSikr8B3gdkxw2K2UehITWB477/UPA7u01i8rpRqA/xv4bJ7HPGXZrCaSTBNNZvBFkoyGE/hjKbLjXeLtFipcdqwWhdaabm+U/X0B4qkMSim6aktpLHficdq4fmkNbqeNRDorW3uFEPNOXoOJ1voVpVTneZe3ASe01qcAlFI/Be7VWj+GmcVcjA+4aDMrpdQjwCMA7e3tMxj1hWXGA0cigzeaYCycxB9NodGgwWqx4LJbqCpxYDmvjUkwluK9Hj8joTgAVR4HXTUeKt12VjWV01btxm6Vbr5CiPmrEDmTFqBn0uNe4JqLPVkp9XHgDqASM8u5IK3148DjYNqpzGSA6UyWSDJDNJnGGzYzjkAsTa5JF1aLhRK7lWqPA8slZg/pbJYjAyGODYXIao3daqGj2k1bVQkrm8rpqPFIS3ghxIJQiGByoU/fi374a61/AfxiSjeehWN73+/xc9YbAa0Ajc1qAkdtqeMyy06acDxNLJ2hxGYlGE+xtydANJlGa6grc7KioYw1zeV01npw2S938qIQQswfhQgmvUDbpMetQP9s3Fhr/RTw1NatW7843XuMRhJUuBxXOGPQ9HhjhOIp4qkMJ4bDhBJpyl02HDYrq5rK2NpZzZLaUkocEkSEEAtPIYLJbmC5UqoL6APuBz4zGzeejZnJdPgiSc6MhvFFUwwE4mSyWTIaltWVcue6RpbXl+FxFv0ubCGEmLZ8bw3+CbADqFVK9QJf01p/Xyn1KPA8ZgfXE1rrg7PxfrMxM5mKdDbLWO5o3OFQnD5/jEg8jdaaVEZT6XawuaOSm1fUsbKxPJ9DEUKIopDv3VwPXOT6M8Azs/1+M5mZjDdaHAslKHXacNgcjKd3MlmNN5JgJJRgOJTAG0mSzZ0Do7XZImyzKipcDpbUl7KqsRyP00Z9mWsWf3VCCFG8FtTay3RnJuONFr2RJCPhJN5wgl5fDLtVMRxKMhZJkhkvHsHsFnDbrZSV2KkrdbCqqRwN56rXQc4TEUIsLgsqmEyXP5piNJTgxEiYQ/1BwvEUGqgosZ9LxJe77JSX2Cl1WqkuddJaVUJHtYdqj0nWy3kiQojFbEEFk+kuc0WSaZSC40NhfNEkFsDjtNFY4aK1yo3HYcFhs1JX46b+ewAABYNJREFU6qK9poTaMucHzhLJvbecJyKEWLQWVDCZ7jKXx2FDKcXq5nKsVih1WLFYLDRXlNBWXUJnrYfaUqfUhgghxEUsqGAyXZVu+7nuvNH/v727C5WijsM4/n3sRetUx7Iw0sgko5ebkrAICoOoKHvBoMIuMsQyqOuUriW66CI1EKOoi16Ig5RCQVdi0E0aBYYUVkSHgrJUKMHUfl3slMth9+zO+Y/z3519PrAw8585u7/H2d3fzou7x05wIoJrLj2fZYsu4lxf0mtm1lOj3ilnephLElfPP4/DR49z2dw5XDDnLJ/zMDMroVFfDBUROyPiqfHx8dJ/+985jyvmjXHhWK+vTjEzs3aNaiZmZpZHo5qJpPslbTty5EjuUszMRkqjmknKYS4zM5u5RjUTMzPLw83EzMySuZmYmVmyRjUTn4A3M8tDEUk/lz6QJP0G/FjMjgPt3aV9vtP0xcDBhIef+nhl1+m0rNfYIGSabr1RydRpvEwmGI7nXz8Zm5ipfbqpmcYi4pIedXUWEY2+Adu6zXeaBvZU+Xhl1+m0rNfYIGSabr1RydRpvEymunKlPv/6ydjETFPyOdOUW6MOc3Wxc5r5btNVPl7ZdTot6zU2CJmmW29UMnUarzNTv/eV+vzrJ2MTM/VbSz8al6mRh7lSSNoTETflrqNKzjQ8mpjLmYZDaqZR2DMpa1vuAk4DZxoeTczlTMMhKZP3TMzMLJn3TMzMLJmbiZmZJXMzMTOzZG4mPUgak/SWpNckPZ67nipIWizpdUkTuWupiqSHim30oaS7ctdTBUnXStoqaULSM7nrqUrxmtoraUXuWqoiabmkT4vttTx3PVWQNEvSRkmbJT3Ra/2RbCaS3pD0q6R9U8bvkfSNpAOS1hfDK4GJiFgLPFB7sX0qkykivo+INXkq7V/JTB8U22g18GiGcvtSMtP+iFgHPAIM7GWoJV9PAM8D79dbZXklcwXwJzAHmKy71n6VzPQgsAA4Tj+ZUv7H47DegNuBpcC+trEzgO+AxcDZwFfAdcAG4IZinXdy115FprblE7nrPg2ZXgaW5q69qky0PsB8BqzKXXsVmYA7gcdoNf0VuWuvMNesYvl84O3ctVeUaT3wdLFOz/eKkdwziYjdwB9ThpcBB6L1qf1v4D1anXkSWFisM7D/XiUzDYUymdTyEvBxRHxRd639KrudImJHRNwKDOwh1pKZ7gBuAVYBayU14jUVEf8Uyw8Bs2sss5QZvPcdKtY52eu+z6yy0CG3APipbX4SuBnYBGyRdB/VfkVEHTpmkjQP2AjcKGlDRLyYpbqZ6badnqP1qXdc0lURsTVHcTPUbTstp3WYdTbwUYa6UnTMFBHPAkhaDRxsexMeFt221UrgbmAusCVHYQm6vaZeATZLug3Y3etO3ExOUYexiIi/gCfrLqYi3TL9Dqyru5iKdMu0iVbjH0bdMu0CdtVbSmU6Zvp/IuLN+kqpVLdttR3YXncxFemW6SjQ97nVgd3FzGASuLxtfiHwc6ZaquJMw8GZhkcTc1WSyc3klM+BJZKulHQ2rZOEOzLXlMqZhoMzDY8m5qomU+6rCzJd0fAu8AunLnlbU4zfC3xL68qGF3LX6UzONAy3JmZqaq7Tmclf9GhmZsl8mMvMzJK5mZiZWTI3EzMzS+ZmYmZmydxMzMwsmZuJmZklczMxM7NkbiZmZpbMzcSsBpIWSdpf/Brk15I+kXRO7rrMquJmYlafJcCrEXE9cBh4OHM9ZpVxMzGrzw8R8WUxvRdYlLEWs0q5mZjV51jb9En8e0LWIG4mZmaWzM3EzMyS+SvozcwsmfdMzMwsmZuJmZklczMxM7NkbiZmZpbMzcTMzJK5mZiZWTI3EzMzS+ZmYmZmyf4FaYbfO23mdRkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = T.plot_feature('n')\n",
    "ax.set_xscale('log'); ax.set_yscale('log')"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}