Drift_Fourier_Analysis_Bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trend analysis with bootstrap resampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python / Numpy implementation of the trend analysis presented in [Gardiner et al., 2008][1]\n",
    "\n",
    "The following model is used to fit the annual trend (drift) + intra-annual variability (i.e., fourier series tuncated at $n$ degrees):\n",
    "\n",
    "$$ F(t, p, \\alpha, \\mathrm{{\\mathbf \\beta}}) = p + \\alpha t + \\left[\\beta_1 \\cos(\\frac{\\pi t}{L}) + \\beta_2 \\sin(\\frac{\\pi t}{L}) + \\beta_3 \\cos(\\frac{2 \\pi t}{L}) + \\beta_4 \\sin(\\frac{2 \\pi t}{L}) + \\cdots +  \\beta_{2n-1} \\cos(\\frac{n \\pi t}{L}) + \\beta_{2n} \\sin(\\frac{n \\pi t}{L}) \\right]$$\n",
    "\n",
    "where $p$ is the intercept, $\\alpha$ is the trend (slope), $\\mathrm{{\\mathbf \\beta}}$ are the truncated fourier series parameters, and $L$ is the half-period (= 0.5 if $t$ is in years).\n",
    "\n",
    "Model parameters are estimated using the weighted linear least squares method. Bootstrap resampling can be used to evaluate the confidence intervals on the parameter estimates (more robust in case of non-normally distributed random effects in the data).\n",
    "\n",
    "Author: B. Bovy | GIRPAS, ULg\n",
    "\n",
    "[1]: http://dx.doi.org/10.5194/acp-8-6719-2008"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pprint\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fourier_basis(x, degree, half_period):\n",
    "    \"\"\"Returns a 2-d array of fourier basis.\"\"\"\n",
    "    A = np.ones((x.size, 2 * degree + 1))\n",
    "    \n",
    "    for d in range(1, degree + 1):\n",
    "        A[:, 2*d-1] = np.cos(d * np.pi * x / half_period)\n",
    "        A[:, 2*d] = np.sin(d * np.pi * x / half_period)\n",
    "    \n",
    "    return A\n",
    "\n",
    "\n",
    "def fit_driftfourier(x, data, weights, degree, half_period=0.5):\n",
    "    \"\"\"\n",
    "    Fit y = f(x - x.min()) to data where f is given by\n",
    "    fourier series + drift.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x : 1-d array\n",
    "        x-coordinates\n",
    "    data : 1-d array\n",
    "        data values\n",
    "    weights : 1-d array\n",
    "        weights (>=0)\n",
    "    degree : int\n",
    "        degree of fourier series\n",
    "    half_period : float\n",
    "        half period\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    intercept : float\n",
    "        intercept at x.min()\n",
    "    slope : float\n",
    "        slope (drift) for the normalized data\n",
    "        (x - x.min())\n",
    "    pfourier : 1-d array\n",
    "        Fourier series parameters for the\n",
    "        normalized data\n",
    "    f_drift : callable\n",
    "        Can be used to calculate the drift\n",
    "        given any (non-normalized) x\n",
    "    f_fourier : callable\n",
    "        Can be used to calculate fourier series\n",
    "    f_driftfourier : callable\n",
    "        Can be used to calculate drift + fourier\n",
    "    residual_std : float\n",
    "        estimated standard deviation of residuals\n",
    "    A : 2-d array\n",
    "        matrix of \"coefficients\"\n",
    "    \n",
    "    \"\"\"\n",
    "    xmin = x.min()\n",
    "    xnorm = x - xmin\n",
    "    \n",
    "    # coefficient matrix\n",
    "    A = np.ones((x.size, 2 * degree + 2))\n",
    "    A[:, 1] = xnorm\n",
    "    A[:, 2:] = fourier_basis(xnorm, degree, half_period)[:, 1:]\n",
    "    \n",
    "    # linear weighted least squares\n",
    "    results = np.linalg.lstsq(A * weights[:, np.newaxis],\n",
    "                              data * weights)\n",
    "    \n",
    "    params = results[0]\n",
    "    intercept = params[0]\n",
    "    slope = params[1]\n",
    "    pfourier = params[2:]\n",
    "    \n",
    "    f_drift = lambda t: slope * (t - xmin) + intercept\n",
    "    f_fourier = lambda t: np.sum(fourier_basis(t - xmin, degree,\n",
    "                                               half_period)[:, 1:]\n",
    "                                 * pfourier[np.newaxis, :],\n",
    "                                 axis=1) + intercept\n",
    "    f_driftfourier = lambda t: f_drift(t) + f_fourier(t) - intercept\n",
    "    \n",
    "    residual_std = np.sqrt(results[1][0] / (x.size - 2 * degree + 2)) \n",
    "    \n",
    "    return (intercept, slope, pfourier,\n",
    "            f_drift, f_fourier, f_driftfourier,\n",
    "            residual_std, A)\n",
    "\n",
    "\n",
    "def cf_driftfourier(x, data, weights, degree,\n",
    "                    half_period=0.5, nboot=5000,\n",
    "                    percentiles=(2.5, 50., 97.5)):\n",
    "    \"\"\"\n",
    "    Calculate confidence intervals for the fitted\n",
    "    parameters from fourier series + drift modelling,\n",
    "    using bootstrap resampling.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    nboot : int\n",
    "        number of bootstrap replicates\n",
    "    percentiles : sequence of floats\n",
    "        percentiles of parameter estimate\n",
    "        distributions to return \n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    perc : dict\n",
    "        percentiles for of each parameter\n",
    "        distribution\n",
    "    intercept : 1-d array\n",
    "        intercept estimates from bootstraped\n",
    "        datasets.\n",
    "    slope : 1-d array\n",
    "        slope estimates\n",
    "    pfourier : 2-d array\n",
    "        fourier parameters estimates\n",
    "    \n",
    "    See Also\n",
    "    --------\n",
    "    :func:`fit_driftfourier`\n",
    "    \"\"\"\n",
    "    \n",
    "    # 1st fit without bootstraping\n",
    "    results = fit_driftfourier(x, data, weights,\n",
    "                               degree, half_period)\n",
    "    f_driftfourier = results[5]\n",
    "    A = results[7]\n",
    "    model = f_driftfourier(x)\n",
    "    residuals = data - model\n",
    "    \n",
    "    # generate bootstrap resamples of residuals\n",
    "    # and new datasets from these resamples\n",
    "    boot_dataset = np.empty((x.size, nboot))\n",
    "    for i in range(nboot):\n",
    "        resample_i = np.floor(np.random.rand(x.size) * x.size).astype(int)\n",
    "        resample_residuals = residuals[resample_i]\n",
    "        boot_dataset[:, i] = model + resample_residuals\n",
    "    \n",
    "    # fit all bootstrap datasets\n",
    "    results_boot = np.linalg.lstsq(A * weights[:, np.newaxis],\n",
    "                                   boot_dataset * weights[:, np.newaxis])\n",
    "    \n",
    "    params_boot = results_boot[0]\n",
    "    \n",
    "    # compute percentiles\n",
    "    perc_boot = np.column_stack(np.percentile(params_boot,\n",
    "                                              percentiles, axis=1))\n",
    "    \n",
    "    perc = {'intercept' : perc_boot[0],\n",
    "            'slope' : perc_boot[1],\n",
    "            'pfourier' : perc_boot[2:]}\n",
    "    \n",
    "    intercept = params_boot[0]\n",
    "    slope = params_boot[1]\n",
    "    pfourier = params_boot[2:]\n",
    "    \n",
    "    return perc, intercept, slope, pfourier\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A test with a dataset generated from given trend and intra-annual variability, with random noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fitted trend (slope): 0.10121480221453186\n",
      "fitted intercept at xmin: 1.9929242719095708\n",
      "fitted fourier parameters: [ 0.3725897  -0.13916781 -0.0016657   0.00325461]\n",
      "std of residuals: 0.04965929108657873\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7efc7cac72e8>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAFhCAYAAABK5GKRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4VGXWwH+TOmkkkIQWCJ0BBARFICqoWAApKwIyohQF\nRFYRlRVcWesKCK5YdhUV0QB+OoAQpChFUAISUESa4EAggEAgvWdS5/tjMiWTSYFkWnJ+z5Mnt7z3\n3vPOnbnnvuec9xyFXq9HEARBEATXwcPZAgiCIAiCUB5RzoIgCILgYohyFgRBEAQXQ5SzIAiCILgY\nopwFQRAEwcUQ5SwIgiAILoZXVTtVKpUnsAzoDOiBJ7Va7R8W+0cALwPFwOdarfYzO8oqCIIgCA2C\n6kbOw4FSrVZ7O/AvYL5xh0ql8gaWAPcCdwBPqFSqpvYSVBAEQRAaClUqZ61W+y0wvWy1LZBusbsr\nEK/VajO1Wm0RsBcYaA8hBUEQBKEhUaVZG0Cr1ZaoVKpoYBQwxmJXIyDTYj0bCK5T6QRBEAShAVKt\ncgbQarWTVSrVXOCASqXqqtVq8zEo5iCLZkGUH1lXQK/X6xUKxXULKwiCIAhuyDUrvuoCwiYArbRa\n7UIgHyjFEBgG8CfQSaVSNQZyMZi0365SOoWC5OTsa5XRbQgPD5L+uTHSP/elPvcNpH/uTnh4UPWN\nrKguIOwboJdKpdoNbAVmAaNUKtW0Mj/z88A2YB+wXKvVJl6zBIIgCIIglKPKkXOZ+XpcFfs3A5vr\nWihBEARBaMhIEhJBEARBcDFEOQuCIAiCiyHKWRAEQRBcDFHOgiAIguBiiHIWBEEQBBdDlLMgCIIg\nuBiinAVBqDUlJSXMnDmdGTMeZ82ar9i7NxaAdetWA1BYWMjmzRtqfL4lSxbx+++/2UXWmrJ8+Sds\n2LCu3LasrCx27Nhap9eZM+c5rlyRFBFCeWqUvlMQBPdhza54fv0zqcJ2T08FJSV6G0dUzy1dmvLQ\noI6V7k9OTiYvL4/ly1eV275y5eeMHj2O1NQUNm36luHDH6jR9Vwhza8tGeLjT7F3byz33jukrq9W\nx+cT3B1RzoIg1Jr//GcBFy9e4O23FxAaGkaTJqFkZWWSlZXFO+8sori4mHPnzhId/Rljx6pZuPAN\nsrKyAHj22X/Qvn1HNmz4ho0bYwgJaYJOl8+dd95d7hpnz8bzv/+9R0lJKZmZGfzjHy/SvXtP1OpR\n9O17C6dOxdO4cRPmz1/M1q1biIv7mYKCAi5fvsgjj0xi6NDhPP30E8yZM4/IyDZs2PANaWlpPP74\nE3z88f/Qak+SmZlJx46deOmlV232c+XKzzlzJp6NG2M4duyIqY+LF7/H//3fCo4ePUxpaSnjxo3n\nrrvu4emnn6BzZxVnz54hNzeXf/97Ec2bN2f58k/Yt28voaFhJCVdtfv9EdwPUc6CUM94aFBHm6Nc\ne+Yv/sc//smrr77ECy+8xOeff4pCoWDixMdZt241s2fP5cqVRM6ejWfy5Kl89NEH9OnTlwceGMNf\nf11g4cI3mD//bdas+ZqVK1fj4eHBzJnTK4xcExISePrpZ2nfviM7dmxly5ZNdO/ek8TEyzz77LN4\nePgzY8YUTp48gUKhIDc3lyVL/svFi38xd+5zDB063OqchuW8vFwaNWrEu+9+SGlpKRMnjiMlJdlm\nPydNmsK3365n5MhRHD9+lJtv7stDDz1MXNzPJCZe5qOPPqOgoIAnn3yMW27pj0KhoFu37jzzzGw+\n/fQjfvhhK7fc0p9Dhw6yfPkqCgoKmDix0iSMQgNGlLMgCLVGr6/aXG65/+zZeH7//SA7d+4AIDs7\ni0uX/qJNm3Z4eRkeST163FjhnGFh4URHL8fX15e8vFwCAgIBCA4OoVmzZiQnZ9O0aTMKCwsA6NSp\nMwDh4U0pLCysVCYfH1/S0tJ47bV5+Pn5k5eXR3FxcaX9sJQrMrKNqU9a7Z/MnDkdMPjgExMvA9C5\nswqApk2bkZaWyoUL51CpugDg6+tLly7dMNcTEgQDopwFQahTLBWYUY8pFB6UlpYC0KZNO7p06cq9\n9w4hOTmJHTu20qpVJAkJZyko0OHj48vJk3/Qv/+t5c77/vv/4dVX36RNm7YsX/6JKYiqMve0LZ+x\nj48vKSnJREa24dSpPwkPb8r+/ftITr7K668vJD09nT17fqz0ZcPT07PcPuM12rRpx0033cycOfMo\nLi5m1aoviIhoZWxV7hxt27Zn3bo1lJaWUlJSwunT2gptBEGUsyAIdYJRUSkUCtNy27bt+Pe/X2Hu\n3H9RXFzExx//j0mTHmfhwn+zcWMMubm5TJkynZCQECZNepwZM6bSqFEjPD0rPpoGDx7Kyy/PpWnT\nZnTp0o3U1JQayVO2BsCYMeNYsmQRTZs2Jzw8vMzsfAMrVnzGM888SZMmoXTr1t1k1rZW8BERrTh7\nNp41a74ut//22wfy+++/8dRT08jPz2PgwLvw9/e3KVOnTp25/faBTJs2icaNGxMcHFKDT1doaCiq\nM0fVMfr6XrNT+ue+SP/cl/rcN5D+uTvh4UHXbBqRec6CIAiC4GKIchYEQRAEF0OUsyAIgiC4GKKc\nBUEQBMHFEOUsCIIgCC6GKGdBEARBcDFEOQuCUGukKtX1cz1VqdLSUnnnnUVVtjF+9u7AoEGDKCoq\nqpNzvfrqS5VmeANspkv9/vvN7N0by6FDB3n11ZcAmDfvBQDOnInnyJHf60S2a0GSkAhCPWN9/GZ+\nTzpWYbunh4KS0uvLa9C7aQ8e7Di80v1SlarWV7um1k2ahDJ79twq2xg/+4bG668vuOZjhg41fLct\nXwjnz38bgJ9+2kloaBg33ti7bgSsIaKcBUGoNVKVyj5VqcaMGcFXX62rsD0x8TKvvTaPTz75gkmT\n1PTufTPx8adRKBS89dY7fPPNarKysliyZBFdu97A5s3fotfrmTJlOufOnSU29ify8/MJCQlhwYL/\nmHKaG1m3bnWFNtu3f1/pZ2rdR72+1CQfwBNPTOaNN97Cw0PBO++8RWFhIampKUybNoMBA+6s0L/4\n+NN88ME7fPDBxwDMmfMs06bN4OLFv4iJ+Ybi4mIUCgULFrzNmTPxLF36X3x8fBg5chSfffYxX321\njr/+Om/z+5KTk8NLL71AenoanTureO65OSxf/gmhoWG0adPWJMPIkYP5/PMv+f77zfj4+KBSdWHJ\nksUsW7YCgFde+ScPP/woXbveUOn9qw2inAWhnvFgx+E2R7lSlcp9qlItWPA6iYmXSUtL4/nnn8bP\nz5e33/6vTZny8vK4554hPPvsC7zxxsvs37+PSZOmsH79Gp5/fi7ffbeJRo2CWbjwP+j1eo4ePcx7\n732EQqHg+edncvLkH/TocaPpfHq9nqysrAptqvpMrft4zz2Dy8lo+Nz1XLhwAbX6UXr3vpnjx4+y\nfPknNpVzx46dKCws5MqVK3h5eZGZmUmnTir279/H22+/h6+vkrffXsCBA/sJDw+nqKjIpDQ/++xj\n9Hp9pd+XvLw85syZR0hICK+88k/27o21aSVRKAzFVu6/fwShoWF07XoDvr6+nDuXQJMmTUhMvGw3\nxQyinAVBqAOkKlXdVqUyjtzHjh3Ju+9+SMuWTap8sbK8hq2+tm4dCRiUpJeXF6+99hJ+fv4kJ1+l\nuLiYZcuWcvToYRQKBe+995HNNlV9ptZ9tPW5gYImTUJZufJzNm/+FoVCQUlJSaV9Gj58JFu3Gkat\nw4aNBCAkpDFvvvkafn5+XLhwnu7dewLm+2BEoVBU+n1p27YdISGGfObdu/fgwoXzlcpglN14z0eO\nHMV3322iWbPmDBlyf5XH1RZRzoIg1ClSlaouq1LVzBdtq6+WXfDwMMT+xsefZs+e3Xz6aTQ6nY6p\nUyeg1+uZNm2GqW1lbSq7Tnk5jS88PqSnp1FaWkpubm7Zi4qe5cs/ZsSIUfTvfytbtmzk++83V9qn\nu+8ezDPPPImnpyfvvvshOTk5fP75p6xfv4XS0lKef/7pSuXS6/WVfl8uXrxAVlYWQUFBHDnyO3/7\n22iOHz9aqRyW9/zOO+/mq69WERwcwptvVh2QV1tEOQuCUCdIVaq6r0q1du23NexfRQyf/cv06dPP\n1LZ169b4+fnx1FPTCA4OoXPnLhU+R1ttUlJSbFzT1vUN975Jk1BuuaUfU6dOJCKiFa1atQbgrrvu\n4cMP32PtWg033NCd7OysSuX38/OjU6fOlJaW4ufnBxgsKtOnP0bjxo1p3boNqakptGjR0qZclX1f\nQkIas3Dh66Snp9Or10307duf48ePlvv+Wp5HperChx9+QNu27ejd+2Z69bqJzMwMgoKCqvz8a4tU\npapDGkBlFemfG1Of+1ef+wbSP1fi3XcXc8cdg7jppj41PkaqUgmCIAiCnXj++afJzs6+JsV8vYhZ\nWxAEQRBqwJIl/3PYtWTkLAiCIAguRpUjZ5VK5Q18DrQBfIE3tVrtJov9o4CXMITofa7Vaj+2o6yC\nIAiC0CCobuT8CJCs1WoHAkMA6zH9EuBe4DZgtkqlCq57EQVBEAShYVGdz3kt8E3ZsgdgPTO/CAgB\nSjHEnTs09FsQBEEQ6iM1mkqlUqmCgG+BT7VarcZi+3RgAZALrNNqtc9VcypR3oLgIhQXF5Oenk5o\naKgpSUVdMWLECDZt2mRz3+rVqxk9enSFfM6CUI+55qlU1f46VCpVa2A98KGVYo4Ensbgj84DvlSp\nVGO0Wu03ts9kwF3msl0P7jRX73qQ/rk3lv376adj/Otff5GY2AqV6jQffRRF27Yt6+xaxcUllX6W\nS5d+zO2334OPj0+dXa8h3bv6SEPo37VSXUBYM2A78HetVvuj1W4lUAIUaLXaUpVKlYTBxC0Igouz\nYME5Tp0aD8DBg7exYMH/8emnI6/7fDqdjjfeeJnMzAwiIlpRWlrK4cOH+OKLZZSWlpKfn8+rr77J\nkSOHSE1N5bXX5vHmm4tYvHg+SUlJpKamlGXNmlH9xQShAVDdyPklIBh4RaVSvVK2bRkQoNVql6lU\nqhXAPpVKpQPigWi7SSoIQp2Rnl4+tWRmprJW59uw4RvatWvPtGkzuHDhHC+88Cznzp3l5Zf/TVhY\nGKtWfcGPP/7AxImPs2LF57z++gKSkq7SvXsPhg9/gIKCAkaPHibKWRDKqFI5a7XaWcCsKva/C7xb\n10IJgmBf+vRJ5/z5AsAXT8+rREXVLhzkwoXzREXdBkBkZFtCQhoTFhbOe++9jb+/P8nJSfTs2avc\nMUFBQZw8eYJDh37D3z+AwsKiWskgCPUJicgQhAbI+++PpGXL9Vy54sONN3oxbdp9tTpf27btOXbs\nCAMG3MmlSxfJyEhn0aL5rFnzLX5+fsyf/5qpKpVCoaC0tITvvttMYGAQL7zwEhcv/sWmTTF10TVB\nqBeIchaEBoiPjw8vvzy8zs73wAOjWbjwDWbMmEKLFi1p1CiYAQPu5KmnphIWFk5kZFtTVaAbb+zN\nCy88y3PPzeH11/+FVnuS5s1boFJ1JSUlhbCwsDqTS3ANdDodGk0sAGr1QJTK2rlRGgJSlaoOaQgR\nh9I/96U+968+9w3cu386nY5x42KIi3sMgKioL1i9elQ5Be3O/asJUpVKEARBcCk0mtgyxewNeBMX\nN9k0ihYqR5SzIAiCILgYopwFQRAEu6FWDyQq6gugECgkKioatXqgs8VyeSQgTBAEQbAbSqWS1atH\nodEY0rmq1aMkIKwGiHIWBEEQ7IpSqWTy5NpN12toiFlbEARBEFwMUc6CIAiC4GKIWVsQBEFwCsbk\nJEFBSoYN6yu+aAtEOQuCIAgOpybJSRoyYtYWBEEQHI4kJ6kaGTkLgiAIdY7k064dMnIWBEEQ6hSj\nyXrOnJHMmTOSceNi0Ol05dpIcpKqkZGzIAiCUKeUN1lDXNxknn56Af36qfD29jEp4eHDG9G8+SLu\nuqs7Dzwg/mZLRDkLglAvEDOqa7NxY282bkwGhrNu3UoUCk/++Ot+2t7hhe7MCUbqb3W2iC6FmLUF\nQXB7amJGFRyHtckavgQGAxOBWH491JyQhxIZ8s5CuozcSWDURTSrJRjMElHOgiC4PdHROyXy14Uw\n5tMeNWoR8D0wHjBYMrz8ihi96geCWiaZ2p/67g48KlFHOp2O6OjtREdvb1AvXGLWFgTBLRCztXuh\nVCp5//0nuXIlhrg4D6AQT59o/rbsaIW2ne/fzZ/egfzv8AUiAlvQKrAlEYEtCPYMYrx6o2kudExM\nw5kLrdDr9Y68nj45OduR13Mo4eFBSP/cF+mf61JdwoqgIG8GDVpJXNzksv3R9eoh7u73TqOJpaCo\ngKMd9pTbF1gQTCNdE5qpAjifcYk0XXq5/d7Fvmx47kXykpuXbSlk8eJNbldEIzw8SHGtx8jIWRAE\nl8dW9K9GY35IS1lC18VYkeqpXXPKbb8n8g5GdRwGmF8+8oryuJSTyMWcRM5mnuNQ0lFunqZhz4JZ\nwDXrN7dGlLMgCPUCKUvoupxOP1NuvV/zm3mgw/0V2vl7+9OpcQc6Ne7Ana1uI68wH3pqaTPwZ87H\n9i2bCz3KUWI7FQkIEwTB5akuYUVDDRpyB/R6Pe/9/olpvUOjthT/EsyKFTuqvFcKhYLxXcfg4+FD\n1PTVzP/PN/XKVVEdMnIWBMEtMCas6NOnIxMmmB/SOp2OsWPXsnv3BEDH0qXzmTq1FxMm3N1gHuSu\nzJd/ri23HrswmH17HgDMAV4QZPPYUL/G/K3DUNae/hb/W3QN6n7KyFkQBJfGGAw2b95oYmJeZPPm\n8oFRGk1smWIuAdaSkPBv5s0bLXOdXYBSfSn7Ew+a1rsl9mXfnilUN+XN0hLSN6w3zQOaEZd4kKu5\nSRXa1ldEOQuC4NJUVb1Ip9MRF3cC+A7YCkyw2U5wDh/8/qlpOSKwBQFFjao9xjqhzMPqbxna+m70\n6EnIumBPcV0KUc6CILglxod4TMyLwP3AAWeLJFhxOuOsaXl0xxE1KnZh62Xs+PYkXun/Arc06+1I\n8Z2KKGdBEFyayh7o1g9x+Cehof+u0E6wLzUJxvMrDCDSL8I05W3x4k0sXrzpmgK8mvmH4+nhWZei\nuzQSECYIgktT8znMSp5/vive3jLX2VFYJ4dZunQJ27aNJyQkBJ1OR/bBpgT1SWJP9AjGfbSeESOC\nTVWpKrs3avVAYmK+KJdQpqFMn7KkygxhKpXKG/gcaAP4Am9qtdpNFvtvAd7BMDv8EjBRq9UWVnE9\nyRDmxkj/3Jv61j+zYpgM1L+sYJa46r2Ljt7OnDkjMSaHgULatXuF3btn85XmR/aF/AIK2PLUAkoK\nwZBn+/4KGd6s+2edqhVw69St15MhrDqz9iNAslarHQgMAf5n3KFSqRTAp8BkrVY7ANgJtLtWAQTB\n1ZA5s+6BUqlkxYqhjBq1CLV6CStWDHW7h3Z9JCGhPxpNLKkBV/FtlEvCrtsoKfQp2+sFlBAX15RZ\nsz6u9PdlTChjTCrTECuOVWfWXgt8U7bsARRb7OsMpALPq1Sq7sAWrVarrXsRBcFxWJvpGlKifXdD\np9MxadL3xMW9CMD581+wYsVQNmz4haKiQkCBt7e3W4603AW1eiBLly4hIeG5si1fAmPQs53C1mno\ncyB+exSGOICVwBjgK2ACMTHDuHKl6nnOUH3q1vpKlcpZq9XmAqhUqiAMinqexe4w4FbgKeAMsFml\nUh3UarU/2klWQbA7DfVB4I5UvFdqBg9eSkLC34HVwCRAXrDsiVKpZNu28Qwe/AoJCf2BMdx6x5c0\nH9CCg+evcGP4DfSd+zNFRUVs3FjAgQM/YZ7uZv59vfDCaCf2wjWpNiBMpVK1BtYDH2q1Wo3FrlQg\n3jhaVqlUW4E+QJXKOTy88jek+oD0z70JCqr4AA8KUtabfteXfgAolR4Y5jd7AXcDu0hIeB7YjkEx\nmxXAli07ePLJYc4StU5w1XsXHh7EiROvER29E9jDrsbH2Hz+GABjbhxCl3s6AvD88zoee+w9NJry\n98H4m6usfzNnDmPLllXs3v0oAHfc8SUzZ46v9y9bVSpnlUrVDMM3/e82RsRngUCVStVBq9WeAQYA\nn1V3QVcMaqgrXDVoo65oCP0bNqwvUVHlI0WHDRtVL/rtzvfPVoDQV19dBabg4VVEUMtPadniHKeP\nDqI4v+Lx2dk6t+07uO69s74vei89u2I3A9A6sCVN9E3Lyb148RTOn6/4+4KqdcOqVSMsovVHkJ1d\nRHZ2kT26ZBeu58Wqumjt94GxgKUveRkQoNVql6lUqruAtzBEa/+s1Wqfs3EaSyRa241pKP2zfuDU\nlzd0d7x/Op2OVat+4LPPLpSNig21nIcPb8TLr41k9KrZFY4pyPbnh5fak5c8ray9+0dxu+K9s1Vj\ne+wCL+Ku/grAhK4P0b9FH5vHWf++XLF/dcn1RGtXqZztgChnN0b65964W/8yMjIYPPgrEhK6YMgA\nZp6uM2rUIvaeepA7X/5fpcd3vNSLxiVh9eIFyxXvna1pVGM1s0z737tzAd4eNUul4Yr9q0vsMZVK\nEATB4eh0OgYP/rBstFzxAd+nT0c6NtpL7IIZZJxvYdquwPwMjI84jF+/XLdXzO6Cd0CuafmeyDtq\nrJgF24hyFgTB5dBoYklIiCpbuw14FdiAp08Kt929jEceHcRzH0Ry14vLCGmTSPuQSF685Vn+c+vr\nXHy/l+k8sYlx7LkQ54wu1Hus06oOeXWhad+g1gOcJld9QV5tBEFwUQZhSFDoCbwOilIeXPksAEtP\n5HA28xxKb18e7DCM0TfeR2pqLtHR24mLexyvw8WM+sLgo9bEx9A/og/ent6VXkm4dqzTqv4aaTZL\nB/tWX31KqBoZOQuC4HIYRmVfA6EYp0Xd+ep/TfvPZp6jV3gPXu7/D+5sdRseHuUfZcX5fmz7x8um\n9Zf3LUSoe4yZvO4d08O07cmek50nUD1ClLMgCC6HcVQ2atRJACL6HSK8yxnT/id7TmZajwmE+AaX\nO87S1Jp1MYyzu0IByC7KQZt82mHy12dspbdd8Mu7pv2dAts7S7R6hShnQRBcEqVSyaJFk2nXbgm3\nPrfMtH1Bv3n0COtW6TEGpb4I+J7fPv2Xad8Hx5bh4Nkp9Q7j9CnLPNcZGRnl2hhzX0uO+tohylkQ\nBJfEmDu7zQRzBPZLPZ8lOCC4iqMMCjoqqhuG6VdK9r79pGlfbnFeufOL8rg2rGtox8VNZu7cL0z7\n9yz6O3Fxk1m16ocGWayiLhHlLAiCS2JUBOHdEkzbdmw+XqNjLc3bib91oSjNFwBdsUFB2BoBivK4\nPjz6ZJqWr/zeHYCDB+MrKHFj4hGhZohyFgTBZVGN3GFa/nbqghofZzRvL168icWLNzO53wgAfrgQ\ni06nY9asj4mLawaUIMqj5lhPn4qKisajQ55Fi3WEhs6gsLDQSRLWH2QqlSAILskDD/Tl10izQr65\nx1rU6lE1Pt4YSQxQUlrC1r92EXf5V1a/ls3PO18sa7UKGI+MU2qG9fQptXoUs/f9DkDaMT/gb6Sm\n/o0tWz4jKOgtsrPnAoYUqtdy7wRRzoJgwuiDhPqVU9sd0el0TH1hHc0eNqzHr/Rn9Yqh131PPD08\nuSfyDlafiiEzvAXmlJOPAt8RFZUkyqOGWL70/HBht2n7rrfmY/5cp5CdvYFRoxYRFdUNtdq9c5s7\nA3ldFAQMymDIkK8q+CAlaMg5aDSxnDxtrqF9eOubbNjwyzWdw/re9W/RB68SHzretwcvP3PpqlGj\njrl9YQxnkJGRQUz8FtO6vsR6rOdDVFQ3Jk++Tz7b60CUsyBgUAa7dxuLwHtLxKkL4O1vUKCnttyN\nvtTzmo61FfBVWlTCfe0G4u2vo8O9P2L0mb7//pOiPK6RjIwM+vRZUW5bnz6fYPRFw0r69btkKu8p\nXDuinIUGi+XIqqioYm1YiTh1Hmr1QG56YC0A53bfXOazrPmD3taUH40mlkFtB6D09OWmcdt4a/EG\nGTFfJ3PnRlOgmGFa3zLzFVq2TGL+/PWMGrWI+fMDWLv2Iflsa4H4nIUGSflatDratHmdjh3PEx9v\nKEkeFRVNnz6diIlxrpwNlWKPEgI75eBTEMi8pw/Vmc/Sz8uPARFR7LjwE9k3/oWvr28dSNswGfmJ\nOcFLXnIonp7eTJs2jGnTnChUPUJGzkKDxDyyKgHWcv78QuLj/067dq8wf/56Vq8exYQJd1eYNiJm\nOsdwMOkwJfoShna7E7V6IBpNbDm/f3WxALam/Bjv3V2tbwfgVHo8287vclSX6hVvvKEut96o0WIW\nLZrsHGHqKTJyFho432MorrAduJuEhDfw9t5kGqVZTxsRM51jOJD4Gx4KD3qGdLOwcEBMzBesWDGU\nSZO+L7dt166J5Y63NeXHeO8sKyZtOruNKz+VSHT+NaDT6XjiiR9pNcuwfnaLPz//PI6QkBDnClbP\nkJGz0CB54IG+tGnzTxSev9NjfBFjNZsZuewfDH7nTbJ800ztjNNGJOLUcSTmXuVC9kW6NenMlvWH\nbKSLjK6wLTp6Z4XzVHXv5vSaaVp+6dV7JdjvGtBoYvnz6t2m9d9WLeD77486UaL6iShnocFhzNnc\n/nFvxvxfMl1GGh7svkElNIpIJrB7CReyL3I1L9nJkjZMDiT+BkDf5jfb7Ro/bjxpWh6yZL4E+10j\n9yz4j8WaotJ2wvUjyllocGg0sXgOLCBMlWFz/2/JR1j06we8sf9t/kyTMoOOQKfTsWzZZqY/+S4/\nnd2H0lNJz7BuNn3HixZNrrBt8uS7q76ADU6sN8yj9mucVZddqfdYx11ILIZ9EJ+z0OAo8iikRe8T\nNWr738PL6HVxIN6lPuKXtBMZGRnce+8Kzp9vTs9HuqLy2kX64SaU9Cup1Hdsa1t2dsXpcJWhVg8k\nJma9af32IZ+gVqurOEKwRdfz/Xhn9TD5XdgBGTkL9Z6MjAymT3+P6dPfIyMjg8OtzObLlD/b8sf7\ngbRMaU/XsI42jz/cKlb8knZCp9MxePCHnD/fg4hbeqEaYYie/mXtRJOZ2ZbvuLaxAAYF/6Bpvf3U\nBFEwNUC09PfvAAAgAElEQVSn0/Hw5DWm9R0bpMCFvZCRs1CvMWYyysoyFDrwGveMaZ++UMEgn1ZM\nWv04SqWS8PAgzl9O4mSalvXxW8goMJfCG/XFXNaq30Oj2WzKKyzUHo0mloSEKAD6zzLXBU491Rao\nWXnI60WpVHJT054cSjpKfnF+9QcIaDSxNJ9itjoZfPWb5DdhB2TkLNRr5s6NJitrLuDNXW+8V27f\n+sdbsvzTw+VGw/7eftzcrBfzb5vHzRfuKte+92PrHCFyA2QQcAUPr1IAzu+9mf79v3CIH/Oh9n8z\nLadlp1XRUjDi4Wm4T5d/6+FkSeo3opyFBkFo5zOEdT5rWv/mkd6UFs/l3LkF3HfflzbN1ePVd3Hx\n/V6m9Y6D95gUhiGAaQvTp7/HsmWbxdx9nRgCvr4moJl55HVPCz1r1oy2u5lZp9Px6MObTesv//qW\n3MtqsHxhOvT5aAkGsyOinIV6zaJFkwkKeotBb5infnzzyH3oS8zzZM+dm21zGo3RL9nj8q2mbXuu\n7ken0/HQQ+uYN280MTEvMm9eLmPHrpGH+nVgDO4a/fZHpm3PTlM7xP9rzhJnZt48f4ktqIKrhebp\nha/P3Su5ye2IKGehXhMSEsLzs1ub1rfOvgl9Sc2/9kqlkicffcC0vuHMd2g0sezfPwWjcoeJHDgQ\nIfNkrxOlUkm+Tw4AnULaO+y6xmInu16dbd6o8JQ5z1Xw+R9fmZYlMY99EeUs1Bsqy7dc3MQQUXpu\ndz+yL00ALgIrMM6T7d//82pNc+/ducC0fCbUvoFKDY2kvBTT8rQeE6toWXfodDo2bkwFVpCqjTRt\n7/2Y7bnvgoGU/FQAWgY0d7Ik9R+J1hbqBeWrTBnyLRtNbt5d8iEN4rffBnjQrx8MHerDkSOL6NOn\nExMmVO/f9PbwolVgSy7mXCYt4ApRAz8lLtZYfmcl/foVoFY/ZN9O1lNe37/YtBzg7V/n59fpdKaR\nsHGuukYTy4EDT2AofPKdqW3H+35GeTIPtXpUnctRn3jshvHOFqHeI8pZqBeUr99rnuLRZ3hr/kj7\nk9YBLfnn9KPAUdTqmteZtXywzxgzmXm/GEbQ9/+rlOG/rufgwdP06dORCRMMUb/R0dsBJGGJi1DZ\nS5sZJfAAVw4n0LzXKQBWrBhaQZkLkJqfblpuGSgjZ3tTpXJWqVTewOdAG8AXeFOr1W6y0e5TIFWr\n1f7TLlIKwnVQoihh6VHD3Nk7Wt9GVL9brul4Ww/2G+c0JrUgnaNpJ/hgykKmTRtWaVsJlqmeUn2p\nafnFW56t8/NX9tJmyBD2BXFxkwEo3uMPZYH5U2evJ3bT3wG5j5Z8+edaZ4vQoKjO5/wIkKzVagcC\nQ4D/WTdQqVTTge6Avu7FE4SaYZ2DuX//z0kI+cO0/+ZmvSo9tjLKP9i9iYtTc3Wtj2n/87EvV9FW\ngopqwsm0U6bl1kEtHXZdY5T44sWbWLx4E6u/NmcLa/bIn8h9rMip9HgAlJ6+TpakYVCdcl4LvGLR\ntthyp0qluhXoC3yClCYRnIjlw/a111aTmBhPelASAKU6T0qLSmp5BR2whpiv5pm2FJcWV95cqBEf\nHfncrue3VTjDGPwn5UCvj+k9JzlbhAZBlcpZq9XmarXaHJVKFYRBUZueTCqVqgUGxf00opgFF0Cp\nVKJWD2TFissUhpkDVrbOtT2PuTrKP9i3ApMAb3b+a46pTULmeRttCyU5wzXSv3kfu5y3wgi5ChN1\n1ytmGRSe+UAh7dq9K/cRyqU3ba2MsDkrQqhbqg0IU6lUrYH1wIdarVZjsWsMEIYh1LE54K9SqU5q\ntdqVVZ0vPDyoFuK6Ppb9M0ztMdQKnjz57nrxdu7q9+/jj2NJSHiesQvNObRzrx4hKEhZI9nLtwli\n166JREfvYPfuP9FoDP7ltPh2phb/+e1D1oxbWq4twOTJE13yfrvS/SsuNVszfE6GEtTPu1afWeV9\nC+KFF0ZXe3wL3+YYqzzf8uQ7/PJhd557rjOtW4dft0x1iTPv3dQNb5iWJ07czO7dEwDYsmUVW7eO\nr5Pvuit9N12B6gLCmgHbgb9rtdofLfdptdr/Av8tazcJ6FKdYgZITs6+fmldnPDwIFP/rAOEVq50\n/8ASy/65KtnZOlCYg4xOfz+Qdu32M2zY7Gplr6x/o0cPZNiwvpw/bxFAlOGDV4hh/nRSUhYKhcLU\n1iBH0TWVMHQErnD/LKPfm99mNtzNnDGMr2rxG6mLvg0b1pdd+wzpPNsMuITiYGOys0N5++11To/a\ndva9yyrIMS0bFLN32fKj/Pe/tS984ez+2ZvrefGozuf8EhAMvKJSqX4s+xuvUqmm2WgrAWEWSICQ\nc1CrBzJikTlQK3PPebZte6rWD1Zr8+jiwWbT9ifHVtTq3A0F4wvrnDkjmTNnJJv+2m6x1/m/EaVS\niSrEXDZUofBk3rzRzJkzUlJ6ltE2tauzRWgwVDly1mq1s4BZ1Z1Eq9XK00lwCZRKJcpI8xv41Clt\n2bDhlzoZ+RgDiKw5lnLCRmvBGutpTUYOrxjjHIFs8ETPicyONcTAnkoZhPUUrIZYGvFixmXT8pR7\nHmBvlNmCZIitkIQt9kDSd9oJ6wChfv2WUVRUJEEUdkZXXGBaTo7pYNeRz5hOI+v0fA0JD+9C0/Lp\nrbfhKkF0Si/zC9zdb77rRElcA51Ox8JD5lKrj03eyooVQ2sUYCfUDskQZieMZlCNZhNFRYVs2mQw\nkYEkNrAnsy3mHv+0ehb2GPkY/abFiiJoDQFedZ9ysj5imfij18TVpu2LF20p2++KvwnDS0RDHSFq\nNLFgTj1OXNxkNmxomBYERyPK2Y4olUpuvL8FXx5az7GzsxETmftjHeg3YtEhaJNNmi6dJsrGTpbO\ntbF8Yf01cp9pu6v9Dp7sOZmPj0YDMGrUW2X5113xxcH+GCt3AZzdeXuFbYL9ELO2nfFQeKDzyeWm\nqauRmDnHMbPHVLvMO7YO9PtjxxAAjiT/UeVxggFrv/3TN051ojS26RDQ1rQc+8c4Nm+uv1HE1VGq\nMCfaObT8QWAl8hxzDKKc7UzPsG50CelE855aIm/fh6v41uojpfpSlJ5K/LyUdA7rWOPkE7Xh8sEe\noIfDycfq/Nz1lXSduSxjlyadnCiJbdas3mtaHvDPpU6PIncmZ1uay6PqS3cBD+Ht7VP5AUKdIcrZ\nzigUCsZ3HY2Phze3zljNqHELGT68kbPFqpecTjmDrkRHQEYTCgsK7ZKe0TrQr3fXGNo2iiQ+I0Gi\ntmvI+79/Ylo2zg93Na4cVQHg7VdQTcv6TY4y02LtXqKiNDKwcBCinB1AqF8ThkQOosSzCP2QdObN\nGy3zJusYnU7Hm59vBGDD0mF2+3xtpYPs3awHAB8fjS5XZUmwTXJ+KgCB3gFOlsQ2avVA9L95mtZv\nHbi8wSskn2KlRGc7GFHODuJyrCGIwtuvgH4zVzZoU5k90Ghi0Tf3oKTQm6Rj3e36+RpH5Gr1QDSa\nWOJ3m+vcHrx62C7XrC/o9WZ/5bM3PVmW4ta18jQrlUq+WqY2rbefGo9GE+sy8jkKy3v11t3/kgIh\nDkaitR1EcVEJScc70bT7aSJvO8iB/xpy01qmM3R2ikB3RueZR3BEIpd/60FJoQ/GKTB2u55V1PZY\njSH6eNXJNfRtfpNdr+1ulEvZOcA8Im3sGeyyNbAtZSjwyWPOnJEuJZ+90el0vPNdNIQY1n09xc/s\naGTk7AB0Oh0bN6ay+80gcpOaUFriwW33fMwDD/Qtl85QTN3XT/vbggG4/Fs3HBF0Zx21fWFfbwAx\na1tRIWXn+W2mfatW7XTpFLctMttarLmefPbCeM8uh5wtt01wLKKcHYBGE8uBA08AEzixvhUenqX0\nGp3Dhg2/uPTDyZ04mXkagCn3X3aKb+zkevP0oIyCzCpaNiysX2Is+ewz13YBNM9qa1r29LGvJcaV\nMN+z8tsExyLK2aEoufDzVHSZgaQ2SqREUVL9IUK1JOWlcCo9njZBrXli4giH+Maso7ZvaG0eEc77\neb5dr10f+P6510hIeIl27ZbgqjWwHx03yLSsGvmdy8nnKL6b9XL1jYQ6R3zODsAybWFpERRqA1D2\nzeF4m71ERWVIEvla8vr+xQComnSspmXdYZntCkCtfpDZ+8wjwejo7WXbG3YcgVo9kHXrPuXAgQha\n32YuoJCT2AwoZOrUSLy9jZ+ha/lzLWW5Ycw23rn1DZeSz16o1QP57rf3Tes922+W55ITEOXsAKwf\n5HfdP5bFx/9HYWkRn31xD1s2uubDyVWxDDAa/VCUafvNTW90qBzW2a6eu2kG7x5aCsCL84ZQWuTT\noIKIKkOh8ASG0X/mMxZbDSNlV0+LqVaNQqONAXBpOesSpVJJ+KgzpvWG/v11FqKcHYTxQW4MtmhV\nVohz3q4P+EA9T778NcQ6SvqHC28RbEj5S0RgC4dcv7Lo+vbBbUzLd732ATvn/bPB51HXaGLZv38K\n1v7mxYs3ucXL6K0t+pqUc1ZuFuvX7gcajkXEU+HptH6W6vWkZuq4lJzLxeQc0rJ0DOkXSdPGDaPQ\njChnB2MMtrilly9t79iPd5MCNJrYBvvwvlasawIH3242ldo725T1i4H1qNhDYQ7haNLhgl1lcWdu\nSOznNt93Tw/z1K+nP/6ErYsM/teGYhGZ1Xu6Q66TkV3AyXNpXEzO5VJKTtn/XAoKy8fldGvbRJSz\nULcYR1xxcSeAkRz89FHa3mF4Cy8sajiRoPbiYdWDdr+G9YuBrVHx7c37sffKgbK1wgYfR2CMtzhX\n0Me07fHRw5wo0fUT1DuZ+lpZztIi9ODY/qbtHULa1u11Cou5lJJrGg1fSs7lUnIOWXnlK115eiho\nHupPq/BAWoUHEBEWSOumgYQG1++XIUtEOTuA8iOu+wgKeovs7Lmc3RlF+7vj2H7oAhN1unr/Fl4X\nWAbX+YWaM3P1a9GniqMcx0NdHjAp53n/ica/oBEaTWyDMYNaY4y3mL3vlXLb3Ikp3R9l+fEvnS2G\n3bC2CP2Y9gaBvWp3zuKSUq6k5ZVTwheTc0jJrDhfOixYSb8bQglr5Eur8EAiwgNo3sQfL8+GPZlI\nlLMDsJ7rmZ2tAr7nrzg97e+GnJAIMW3XEMvguhPNfiG3bLu3h/2/ypYvBmA7ut7SDHqq1e+sVX8A\nNBwzqC3cvc+9w3uYlj19sykp8K1XFhFri1BabiSBpBDmF1rtsdZ+4Usphv9XUvMoKS1fWrKRvzdd\n2zQmIjzApIQjwgJQ+ngRHh5EcnLDLc1pC1HOTsEHuJ/kEx6UFv9Ch3t/pvDSAGcL5TYYg+ue2vUD\nAM0DmjnsuuWnTxmUrXWQWHnqpxn0evnHzU85W4RrxjKWYdpHS2if2t0tgtmuF32pYcQ6o+fkctuz\n8gq5lJRTrV/Y19uTNs2DDObo8EBahRn+NwqQFKDXgihnB2A94urX7xIKxXL2738cMPzwj0TsAUYA\nkm/7Wnn6xikOu5b19ClbQWIfr3yeN39bYmigKAV9wzbP5RblmZbbNop0oiTXj4fCg1J9KakBV3hj\nxPPOFqdOsX4+RfQ8iq+HEu3pYnalnKqRXzgizDwaDg1W4uGipUDdCVHODqDiiOshADSaTaSld+NM\n+DFT2+oighs6xheXdL8kCDdsa6wMcZo8toLEtsVsgjIddNvspfz8n2n1ygx6raw/vdm07Kr1m6vj\n9ai5vLxvodu+XNjC0i88/pmb6HznFnIpRBFcSF5aOCv2a01tw4KV9OoYXM4kLX5h+yLK2UFYj7gA\nJk++j1J9KTN/fBGAI1eO89GrPxAX9yJiDq2I5YvL3z6bi9FI5srZuFr2Oe42c3rtxf4rB50tQq1p\nomxMmF8oV/OSKNWXlps25+qU6vWkZerKm6OTc0is4Bf2ILB5NiVA+0Zt6D+0CxHhAbQMDcDPV1SF\no5FP3MlY/sg/PbGSmJghTpTGtbEcpWb+FUF413hOfdGUtdtGAs6xMlQWJNYs0ZPNCYZ828bkM678\nEuEIHu36kLNFqBWdQtoTl/grl3ISaR0U4WxxbGLyC6cYFPCl5Fwu2vALeyj0BHhBz27NiGwWbPIP\n77i8jV1/waib+9CpcUsn9UIAUc4uwZhOI/nm9MaytcHAKuBRQPJtV4ZfkwwKc7w4sv0lnGllqCxI\nbGCz/iblfDH5Ik9N3dMgXRUlpWal0LdZbydKUns6hrQjLvFXTmecdbpy1hUWczklj4vJOTWaL2z0\nCTcN9mHh63v4+aeJgIJzUebvok6n40D8ERQ+Cpr5hDmnY4IJUc4uwJ2tbjMp58AWmeQkjge+Y9So\nY7z//pMN4iFeE4yj1MOnRhLYLIX0E14uEWxlK0hs/MMbTSlaFx77gLi4xTREV8Xuiz+bli2nmbkj\nnULaAxCffpZBrR0zu6K4pJSraXmcvJjJybMpXEwymKaTM2zPFzb6hY2+YWu/cHT0dn7+aRLW30W1\neiDqR9bRYnoO6WcimfDBlgbzAumqiHJ2ASyDZIa++xpr1e8TFZUkitkK4yh17p5/UwxERbUjN6rq\necfOwGh+Hzvrd4utu4C/OUskp7EufnP1jdyEUL8mNPYNIT4zoc79zjX3C0OQjfnCtfULazSxnEoZ\nSITXEVK1bTkSp2TWrI/lGeRERDm7CF2bdOZk2inAfYoCOItiT4PprlVQS1avvruCSdnZFBUZ5Nv3\n7jRufW4ZAO3a/UxCwlDAdV4iHImjRpr2plPj9vxy5RBXcpNoGdj8us6RlVdYIX2lLb+wcb5wRFgA\nXdqFEuznRatazheuLEZCo4klrHMCACnai8AsYmKGceVKw3HBuBqinF2Ex1XjeSHuNQBuHdZBfgw2\nMERrrzeZi6PnZDJUgwuah/XACi4dmGDa8uhjKoKUrvUSUdcYp7kZXk70eHv7mKaU3d/uXqfKVluM\nfUsOyIdQOJ1xtlrlbOkXtsyglZVbPpe+tV/YOCK2nC9cVxm0KouRUKsHElvyHwBST02mIbpgXA1R\nzi6ATqfjYfW3JqXz/rFPuOHP/vj7BzbYyF5baDSxXORGWnEYgLh9j7nkg8Pb2wcYDvxg2pYccZGZ\nI9w7WrkqrOfnwwo8vAczetUPeCo88fNy3++wZd8Cmycx9L03OZH0Jwk/5AAwduztZOSVmk3S1fqF\nw6r0C9sbWzESX2t249cuC32OB7qMYIfJIlSOKGcXwOCjfLicj/K11xKBp4iJ+VrMShb0n7XCYu07\n4uJOuNwLjMF0+HWZ6dDgc00NuOJUmeyNdTIWmEjjdqsBCM1y7yk55r55UZIXQkFWAEfy40n4MYqg\nsBx2X45D4VE+uYrJLxwWQKumhgxaLcNcb76w8cUjp10TukSCIrCU/v2N2QsbpgvGVajym6JSqbyB\nz4E2gC/wplar3WSx/2FgFlAMHAP+rtVq9bbOJVTHdjZOf5ORn/wLAC+/GRTn7xSzkgVq9UB+3WcY\njSb8GAoMc0m/mKXpMDOpN6ea/l79QfWQMFUaAIEF7jcSs/QLn8r04Db1PgJDc/D2LcazKAjPRrm0\n6n2aoqxgMq6G0LV9OnfdqjKNht0lj7Q5ePEZ07YRI4J58MH67YJxB6p7jXsESNZqtRNUKlVj4DCw\nCUClUvkB/wa6a7VanUql+gqDLW9TpWcTbKJWD2TRohdJTR1p2jbqixdZq24D3O88wVwIU77xMh/m\n4RVzcWW/mKXp8KldBuWcV5SHv3f9LBRvHWgEKwlTGUaTjwwe7DS5qqOgsKSsvrAhQtq2X9iDkObp\nZKcFkZQQROPiw4TeCqeONke7eRBQxN8Wb+LeW1o7qxt1jiv9lhoq1SnntcA3ZcseGEbIRnRAlFar\nNTpWvID8uhWvYaBUKpk58w5ee22F1Z72YlYCMjIyGDz4Q7J9Vdz1mmFbsc793uZXnlzDk1aVfuoL\nltaCoqIi9PjzR9ujBPmG0DzYMVXDqsI4X9hY0tA4Kq6RXzgskJAAD75ZuxfI4NbBI3n/5KcEt00A\nitz6N6pWD2Tduk9N699OW8BfN6xmwgSpL+9sFHp99VZolUoVBHwLfKrVajU29s8Ehmi12mHVnEpM\n3pWg0+m4774VpDbL5oYxZwAYkHQf06YNbdA/Ep1OR8+e/+X06WcZ+ek/8W1kqOB8+OX2nD5tiKC7\n444v2bp1vMt+To9+8wyFJYbpVWvGLXWyNI7hYlYiz3//Bre36csz/R+r/oA6Qq/Xk5Sez/krWZxP\nzOJ8Yjbnr2RxMSmb4hKr+sIBPrRt0Yg2LRrRpnkj2rQIIrJZEP5K70rObvg+fhH9A3uCtwN6BmQO\n5rHJ97jsd68mzP3sPRKCDEUu1qqXAoUsXbqDJ5+s7nEuXAPXXPGl2ugElUrVGlgPfGitmFUqlQew\nGOgIjK7JBetzQe3aTnf4v/97kK81uzmIQTn731JMdnYR2dlF1RzpGJxRED06ejunTz8LeJsUc+pp\nHx55pCn+/ka/2Ig6+Zzs1b+Zvabxzm8fAc79/jvy/h28/AcAEb4Rdrtmdl6hKVlHak4h8X+lcyk5\nF52N+cKtm1rUFy77H2zDL5ybrSM3u+JoGspHbY/VbAGgXVSwQ36j9rx3RsVsSXa2zqHfVWc8WxxJ\neHjQNR9TXUBYM2A7hkCvH200+QSDeXuUBILVHqVSyWOTB3Nw104Atp3fxcgOUgjDgPmB+evSf5AX\nuYU1a+53ixFL++C2puXMgmyCfa/9h+punMk4B0CHkLa1Ppe1X9iYQcvmfOEm/mXmaLMSDquj+sIV\nI9Lh8z+/4uaWvWp9bteh0K3N9PWJ6kbOLwHBwCsqleqVsm3LgADgIPA4EAvsUqlUAO9rtdoNdpJV\naIAYA41Op5sfiNmXW7P/8uMuFwRWEzTa9UzvOcnZYtiVUn0pB678BsD2b46i4FiNprsVl5RyNT2/\nTAmb/cIpGboK/rCwYCU3dgg1TVPqoWqGr0LvsPnC8dsG0nFwrEOu5SjapHaR7IQuRJXKWavVzsIw\nVaoy3DuTvYvySr9/8MYBQ7ae4tJivDxca26kPTFFZWMurbh69Shm73ulQltjmkx34LaW/fj58gEC\n6mm0tk6nY9WqnRw8eJrWd/hAWfKsuXMM+cQtK3Hp9XpSs3Qmk7RBCedyJS23gl840M8bVWRIucxZ\ntuYLO8IsahmRfuTLkW6rnC1/YwNHdDZtnzP2cWeJJNig4Tz13YhmAU1Ny9r0M9wQqnKiNI7D7NN7\nGNjF0qXvsG3bUzbe4guBlWzcWOA2UaVD297Nz5cPkFNkyCpl6yXEXdHpdDz00Dr2758CwNhx5jmz\n3spSGoVnkZh/J298speAJiG19gs7C+vUl7+WbV+yehV//9tYt7iH1pncxkY+U80RgrMQ5ezifHRk\nOR8OWuxsMRyCMVOaYQbfBBIS7ufeexfRokUYrZ81tPn1Ez/gO+BRDhzwcBvTdmNlCGHKJsRnJJCX\nn8fD6m/rTX1njSaWXw9OIrhZDo3CskzbCxJbM/jvW03rV/I98EzMLu8XDgsgomnd+YXtjXH+uk6n\nMyXEORN+jHHjvNziHtrymwuuiShnF+Wmpj05lHTU2WI4gV3ABIwPj/Pnu3M55S5acwSAcz++BewA\nlBhG0O5Dp8YdiEv8lS/Wf1fuAemKSVQqo7xf2GCWPpnsyZCZ27HWrbl/diDrSjhZKUE0b/I7b70+\nkMjmjfH2cn4N7tqi0cQS/+ftdLxvL+Be99AWxlrVgusgytlFeajzAyblXFJa4vaF6muCWj2QpUvf\nISGhfFa0Yf971aplEe4YVdoppD1xib+SrUx3tijVcq1+4aK0XC6d7UyxQkfPvobtOz8ZSrt2C5g6\ntRcTJtSf+fo6nY64uBMc2XwDHU262D0mq1j6zX2DzVaOWb2nO1EqwRainF2UIJ9A0/JvSUfo2/wm\nJ0rjGJRKJdu2PcXgwUtISHgOgH79LuHtV2Bq07//54wYEYy3t/tFlXYsG5006uBJVFTFmrrOItuy\nvrBFBi1rv7CPtwetmwaZArOM/308Svjyy50cSN6A9/0ppvaLF+9ArZ7tVveoOsw+2xfLtuwAYNDE\n91CrXT+gytJvfqTlXpPtSeEGLoWGhihnF8Q6j/SKE5oGoZwBQkJC2L37cYt6sw8xe99xAILzQ3ln\nzWi3fdiH+jUmVNmYhKzzfK2Zy5rVji0uoCsoJiExi4tJ5ZVwZnXzhavwCxuU1QZzgNF4Q4BRY98Q\ntzXxVkVlPtvQ+8+h0cS6RXCf0W/+1K4fqm8sOA1Rzi6GZTRlj/F6uozc6WyRHI5l0QjL9LJv3v8C\nHgr39ld2DGnPgSu/kVacYTflZcsvfCk5l+TMfKyz9YY2MswXtoyQbt7Ev8Z+4cqU1dQej9ZRb1yb\nMzv60+He/QDMmTPS5YP7LGcJGF/+WwY0d55AQqWIcnYxLB92J74ZaVLORSVFeHs2vAjLr7XrTcvu\nrpjBEBR24MpvnM44S6ug2tU5NvqFy5mkkyr3C/foEEZ4sNKkhCPquL6wt3+eablto8g6O68rYV19\nK3PPCbi3bKfC06UDwyxf/JUhmYz42DByfu6mGU6WTLCFKGcXpqTQPMdzR8JPJO01PHDdwXRWW4xv\n+L9GHii3zd37bYyKPZ1xlrta317j467NLxxYzhzdKjyQRv7eNG3aqM4TdVgqq+Ef/dO0vT7cK1tY\nz3UuKurNYXYD0Krf71zc38OZ4lWJ5Yv/DWO3mbb7e/s5TyihUkQ5uxA6nY6ioiLatXuZhISXMEwX\nMrDlwg7WzvkAcP95sdVhmYxkrMbsFxs3LsZt+2182dCjp3HHYOIzzlKqL61gDSgoLOFyqln5Gk3T\n1n5hD4WC5qH+hlFwmDlAKyzEz6HzhY3KatWq9RxWmmV053tVHZZuF51Ox+F9BuUc9exnXHy/t1vM\nIGh/915niyBUgyhnF6F85p7RtGu3hKlTI2nb7n42JHxX1sr95sVeD+Y3/O2mbZtmLESX7u+W/bbO\nyshFtHQAACAASURBVHT/q4cI6JrHkUvnKcr2L+8XzsivMCmntn5he6NUKvH2NrtcDq8Yw+m4AW55\nr64VpVJJ/2Z92H/1IIBLv5BYm+SNREdvbxDWOHdDlLOLYB1Yk5DwHN7emxjUZoBJOXv55VOc33BM\nUM17XTUt69JDcLekI2DwC0f/Xyxnrgyjwy3nCArLorioA3CQj7b/RElSG1NbYx7pCItpSnXtF64r\nrNOP6i1eKU5vvQsodpJkjmeMaoRJOXv7uG5ciNHKsXLtN2UpfWDny8+x9nS7em+Nc0dc71cvlMMy\n+UgPdQy/fzHG6fNi7Y1aPZD165fT+tkjFltdP+lITn5RhWlKl1JyyC/wou+og6Z2JUUBABTodWh/\n6krrpvv48N17adokyC3mm1pbAmJivuD5dyM5eK6sgb7Y5e9VXeLnZX5hjr0Ud02xBI5GqVTS+95W\nHDlhWE87bSh8Ud+tce6IKGcXwdrkZOvh1nHwHh7u6ut2yTeuB72+fKDT/PnrmTDBNfpdUFTC5ZSa\n+4W7RirZu/M8f/w+kOyUQJo1eYtubYPxbpJMwqFOJNCW7we5z4PR2soTF6dmzbm5pv0NuezgN6c3\nurRyBjiVfgaAH176h5MlEapClLOLYB0Favlwm9L9UZYf/xLAbR7gtUGjieXAgSeIxJDQ4ue3p3LL\nuHMOf9iXlJZyNS2/XG3hSym5JKfb9gv37BBaLnOWpV946jAVGk0sRUVFHDjgx5UTrWkz4CBBEVfI\nvhTq0H7VLTpgjWktTxuC+jHxX7oypzPOoPT0pUvzH4g72xpwfpY6oSKinF0IyyhQS24Mu8G0nFGQ\nSYhvsCPFcgqWeX8v/9YDxp2z27X0ej1pWQWcS87lxJkUkzJOTK28vvC1+oWVSiVq9cAyc/Bs2g36\nL20GQHi3k9zQ+hRq9Wi79a+uKW/l2QpMAn4DYMfbc9D47WoQL5GWvB71Iq/GvQVAUWkx3i5agz1d\nl0FyfirdQ7swX/OgzcGA4Bq45jdIKIel3/nrP9cz48bHnCiN/VGrB/Kzcr5pvS7f6nPyi0xm6PJ+\n4UrmC4cFlqsx3CjA57r9wmZzcAnJJw1+56bd9qDfH1bbbjkUSytPXNwJfvjFnFq2MCfAiZI5jzC/\nJqblb+O/Y0znkU6UpnJOpGkBOJ76Z6WDAcE1EOXsJvh4+lBYUsjx1JPOFsXuKJVKfJrmm9avJ4rU\npl84JZfMHNt+4e7tAlC1bUKIvzet7D5feCc5iTPIT3uF8K55bHp/KhrNZrd6UBof7Gr1QGZuWGja\nLuZR+PHiXoZH3mf6zlpHtjtzhPrVn+ucdm3h2hDl7Cb8vedjvPf7J84Ww2H4eSnJL9Yxq/cTVT7M\nrs0v7EvPDqEmc7S1Xzg8PKjOM2hZYjYHNwUUJJ/sRORtBwlskWS3a9obyxcp72LfBjsdR6fTUVro\ngYdPKWBOwmJctoxsd4XP6OamNzr1+kL1iHJ2EzqEtDMt6/V6t5hyc73o9Xp8PX3x8vCiU0gH07a0\nrAIupZiLOVyszi8cFkhE0wBahQXSMiwAf6Vzv+7mbFo/8NlnS0g+0ZfI2w7Sf/gq1OppTpWtLnim\n3xSnKx1nsWrVTr5b8irDPzTUHv/tWAizZn1Mnz6drCLbnTdlKU1nriM+ssMQh19fuDZEObsJHgoP\nbgzvzpHk46Tq0gjzc+cI36q5mJlMRkEmLb06sGqb1qCMbfmFvTxMI+C68gvbG6VSybRpw5kwQccy\nzSb+BNr11TlbrOsmv9gse/vgts4TxInodDo+++ww+anmoL5ek7KJ+fBFDh9+GXCNYL+DVw+bluvz\n86O+IMrZjegc0oEjycc5lX6mXvy4jH5hkzm6zC+c43sOnw5w/qwXZ65cxkOhoFkTP7q3s4iQDg8g\n3MF5pOuaH2KKCJsURGJpAePGrWf16gfdbuS56exWZ4vgdDSaWBIS5gGrTNvaDDjILx9OISHhJdq1\nW0JCwnOA43zytvzc35753u7XFeoOUc5uROfGBhPvqfQz3Nqyr5OlqTlGv7ChpKE5g1ZlfuHwiHwy\ngRG9buLGFh1pERrgMnmk6wpD5PZ4+vd9l9ZRiRxLuBeNJtatgsIAdl/c52wRXAQlMJ7YhecZ+M+k\nctt79Mhl6tT1eHt7O2TKkq0Mbkb/N8DQtvfY9fpC3SDK2Y1oEdCMQO8AtOnxfBG9DQUKp0d/WqLX\n60nPLuBicvV+4QClF51bhxhGwVZ+4ad2fQvA0J498HLR+aK1paioEFhN8onbaR21lvCuaygqalPt\nca7K2M5/c7YITsNy3nfSsWeBlwAIjkwg88JONm58geTkrx0WCFYxg9tkVqxdAxGG/fdEDrS7DELt\nqZ9PvnqKQqGgQ6N2HEk9zr/fjSInsanToj8t5wsbzdGXknPJLyhf8MDSL2xpkg6uxC+cnJdqWq6v\nitmAAphE0okUAMK7hZVtcx9KSs0xAANa9neiJM7FOrvfr2Xb+z39Edvn/BtQOj13dWKj86ZlpZdr\nvMwLVVOfn371kozTxdAEmnY7S05ihN1/9JZ+4bTc85y+kM7F5Byb84UNfuEm5f3CwX54eNRc6aw4\noanrLrgkxhKL2Zeao8sMJLxbPAd/8C8zfbqONaQq9ieai3lYJsppiFgm9Lh04DiXc68QHJmDZU12\nR2ErT39y0EWHyyHUDlHObkaQrjEA4TdoObtzQJ2dt6S0lKT0/HLm6EvJOSRVN184zKCEW4T64+1V\n+wd0QpbhDd/fq36XxrR8gKb82ZZW/Y7z3b5niIkJd5m5sNXxlVYSWthiePv7+PTYSgAUnuvRl9xD\nVJTGYclZzFP21nPw4Glu7tPBVCJycreHHSKDUHtEObsZk0cP4bkdv9K0mxYoICpqxTX96C39wpcs\nUlheTs2juKS0XNtyfuHwALp3aoq/l8Ih84Wf6DHR7tdwJpam0P3JfwEQ3vUcuVdbOt0Eeq2EZ0cQ\nHb3dbUb89qaVV0vTcsQtEXhfXcqKFeMd/tls3pxFXNyLbN+XzP3v7wHg5maSfMRdEOX8/+3deXhU\n1fnA8W8SAhMgISxhURAiyAHEnR0B6wJihAKKRDQSFdRqFbWV+pOWWoWq4FJLrVUUA1EJoISyCLIa\nEAEFWV0OIEGw7EhICAyEJL8/7sydyWSSSSAzc2fyfp6H55ncO8s5OWTee8/ynhATExNDp1Yd2HR0\nK397I4P77yj7Dqty48J1fI4L+zuDVkFhgfn4MsfM9HDm7Ao99dFJvmMdCe13sOeLHsEuVoUUF7v6\nU9594gnOnY4JmTt+f7Lb7dx269t0cmQ07f5kGrOT32Tu3MBebLlPCrtupGvXsMiI8Fr1EM4kOIeg\n9o0uY9PRrbTv0wSbzcaZgkIOHMvnl8P5JTJo5ZQxLnx5YgMjaUejujRvXPlxYX9ZkL0k2EUIitQh\nt/L0im9I6LADOBsS+an18V3m43On44DgZr+yCmPNc3c6YZ31302uMDa7iC6sGeSSiMooNzgrpaKB\nqUBLoBYwXms93+38AOAvwDlgqtb6PT+WtdpzjgvbfzXGnRdv38TnCyMDPi7sL8v2ZgW7CAFnt9uZ\nNfNLGjSoz4mEo7z42ixGDLX+3ec7W9OCXQQLu5EVf93DjX/7EYC2V04iOfmhgJbANadhhHns0atT\nA1oGcWF83TnfAxzRWqcopeoDm4H5YAbu14FOwClgjVJqntY6dLP4W4RrXLjk5CzXuHAxtuuiyKu5\nl4jTV9C2RXyJ7uiLG9UNeh7pivDMYuQ08NLqkffXPVlEm/5xXDNiDu37NLZ8YAY4W2QMQdj3xAJG\nD00o3PH7mxEUZzi6lP8EwLBXYgLeps45De9/+jHfO461TWgT0DKIC+PrG3w28InjcSTGHbJTe2CX\n1voEgFLqS6C32/NFBZQYF3bb4rCsceGmDWws/2w3MV2MNaZ5xesZfccTIfGF7s5bFqPmo41zN1/S\nJ4glCxz3ccEj37UHYOm3X7F7hTGuHwoTrF668xF6RRqdaYHIfmV1rol+S831zj+e2BnwcjgvfHde\ntMX3k4UllRuctdb5AEqpWIxAPdbtdBxwwu3nPKBeVRcwXJwtKGTX/06wdvtB6taOPu9x4bS0Jayc\nM4yhd60BoO5VR0My7aNnFqM9ZzrRnE1A9Vwze2LfRZw9WZsD5w7zzpiBgHW2F/T0v5MHzMeN4xJC\n7v+evzkn+uV8u5udObsBKCouCthkLPcL36EZywBoHdcqIJ8tqo7Pvk+lVAtgDvCW1to9Q8QJINbt\n51jgOD4kJMT6ekpI86zfd7uPMXnWZvYfPUmxx8Bwo/gYOrVvQsumsbRsFkfLpnE0b1yXmtFlB6fY\nWOOL+stJD3L9M+8DYLNFBuz3WlWf46yHU88/uqYrBPP/SCA/+/HHk1i4MJ2srHuhGCKj7NRoUExc\n8yPk/mIsqVq4cCmPPJJUZZ9ZFfV78etJVfp+VcVKZQF4/uYnueeTJwA4VLSfK5u2L/O5drudtLTl\nAKSm3uT1gsxX/ZzvkZW1jbVr/4D713vTva1J6G+t348nq7VfsPmaENYEWAI8qrVe6XH6R+Ayx1h0\nPkaX9iR88OdSnGDzttRo3/4cTtkLuKx5PDVrRHLqzDmG3dimzHHhEzmnyv2MpKQudO36LuvX1zKP\nzZh1gFtv3cfcuV8D/usOrcqlVElJXeje3ZXFKCLSuHLp0FAF7f+Iv5eKeZOePsCV9jHGWGd+y8sT\n+PTetwDIy7NXWZmqqn4HTx4B4KI6Tdm370ip3Y+CIRhtVxnjs/7JWzdO9HrOc4hn+vTSPSa+6lfy\nPYwhsYu7bDLPn82NsPTvx+rtd6HO58LD153zcxhd1eOUUuMcx6YAdbTWU5RSTwOfY4xHv6+1PlDG\n+1Rb16nGXKcaV9n72Ww2Bg5syPr1dwAbATjTLoZ+/T4mO/tpwLrdoe7Kykd8T7s7g1eoIHBP+3hm\nywG2HvuOyBpFhMKSqpS2d3nd/cjK/+8CrffF3Vn1v7XlPsfbRhWVXZJW8j36AdPo8fRG87z7hEsR\nGnyNOY8GRpdzfgGwoKoLJcrnzMvs1Prmr5n93q1AIVZIsl9RzsB07PRxvllrjI3F16qe0xbsdjs1\nNidAC+PniRPnW3KC1XF7jvl4xX+/v+CgEu5+2/o2MzhPmb6A6KKaAehhsAEDcV68RxJhuf9HwjdJ\nFxOCkpN7G13C/xjudvQ24GPAHqRSnb/0H2YGuwhB5eyS/L9nhpjHOvZNsOQX6j83vxvsIoQUWw3X\n8NO8vTsYM6YGQ4fOwm53/Z06/56NJWnOHpPK3em63iMXmEV80zfMc9/8pQU5OTllvlZYkwTnEOTs\nEu7c7AvzmLG74r3AYuLiJjJoUJcgla7ynDNaa0VVzwxGnt2aAB/8OCN4BSrH4VNHzcdVEVSqhWIj\n+94lPQ4Ct7F+fS3S05eZp51/zxMnzmfixPnnNTRgs9mYNq0/iYlvA7W45R+uYLx75x/o1++tEhcE\nwvokOIcom81G9+4dzJ+vuPu/jkc/kpv7hDk5LJQ8dMUI308Kc7s+D4013nerIVUSVKqDwjXxrh8i\nooD72LBhV4nnOId4UlP7nvfvcO7crx3zTkqv9sjO7mZO3BOhQYJzCHO/S2mbtBz4EHiSYOwhe77O\nFbmSrbRrcFkQSxI87negWz4cYB4vKi4q+0VB4H7XfP3F3YCqCSrh7spG7czHl99pXER36uS//+vx\nrVyfN++h8RjfCzf67fOEf1g/x6MoV4vDbdnXeIfjpzuBSMvP8gVXBqP/1fup2qeuKWvm+u4TP9Mm\nPjF4BfPwt3XelwKJstntdhYvPkULR7zscMcS6vz8CykpyVX+Wc582k1GfW8eO5O7GrgzoPtJi6oh\nd84hyjmJ6PXRj5rHxk9YEBLdi86yjxkzkP31soNdHEtwvwN1jr2/8e3bQS6VuFAZGatYt+7BEsMV\nv/2tfyb72Ww2PpoxgOjaZwC4sn4HJk4sZOLEpZb/ThClyZ1ziPI2iWhni81MuvUvwStUBXkr+0U5\nlwavQBZit9tpdrg1exr+EOyilODexd7xQDfsdrt82VfCpg+SadPP2HWtoHbFJ2Y5e5hiY20kJXUp\n8Tv33DjGZrPx3Qltnr+3413UuaZ2FdVABJrcOYeRUzVDN8NOs9xWwS5C0Dl7FCY99oh57NeTPjPi\nBsSuo7vNx399ajjDhmXK7N8KcJ9PsCerMwDNOlfsosa9h+l3v+vLsGGZ5OTkkJa2hClTFnDXXZ8y\nZsxAxowZaLbHdLdliXWiJTCHMgnOIcr9j37/tx18Pt9KnGVvcqVrx5zhyTcEr0AW4epRcC0pm5T1\nTvAK5Gb2l0sB+PaDYUA0a9emkp6+PLiFsigjx/US0tKWAJgz2odcWpeakTXZdHQbxZ6J9r0o2cMU\nzdq1yfTr9zFjxgxk7NjarFv3oNu5VNJnLSv/DUVIkeAcotyXsQxs1Mw8XlBYEMRSVYyz7H2efbfE\nMeFy6mh9AHJjfg1ySQw5tmMUnYvi51VdzWPvvbdZ7p49uN/tOu9oAVJT+/LgiP5cldCRY/Zfyc79\n+TzefYVjqVQ03kYkdya4LnYn9Bxb6rwILRKcQ5hzEtH9qf3MYy99848glqjibDYbxZHGOGbb+rIJ\nPJTsDfnihd8HuziAEWymTJ/PqVq5nNxTg3OnozCSjnxIdvZzsnbWQ+m73dQSv6POTa8B4JuDm7y/\ngRvPJC+Jievczt4ETMM9AUx+LdcOvnNnrCctbYlcPIUwmRAWZg6dOhLsIlTafe3vCnYRLOP222Np\n2vQVrrqqJd85jm3ct5ltyw8Dgd35yXkXuL9mB7r8DnJ+LAYWYXxtDEeu7SuvXf02xEbX5dvDW7nz\nsoHl7l3uvsQuNtbGTTc9xn33vc+6dQ9g/O4LaNVqLKNGXUfH/m34cIcR8H+ZfDWz11h7T3Dhm/x1\nhYkrG10e7CJUivsGCvVt8eU8s3pwBsKxY+8kM/NJXn99r3lu6s6PS038CQTnXWDLnsYGCt+veJzE\nxB+BW3Ctp5d0ne58pTSNiozi2iZXcbIgny0Ht5tj02W1qbN37JFHkoiPj2fAgHoYF0hLgfvZs2cC\n0dHRfLhjtvmatWseoKw7dxE65M45TIy6IoXHVz4LwKmCU9S2+EzNP3/192AXwVJKdocuITf3T+gF\n81C3Oyf5BGfnp4ioQppc+SMARYUxjBx5CdHRRrIUK+6aFWyeCWW8/Y6uqt+BrF/W8OonmSwb/zJQ\n8Tvc6OiaGJvcOJchnqUg8ox5vt3ezlVVFRFkcuccJiIjXE350jdvBrEkoqpsmzHIfBwRda6cZ/pH\ncnJv+gx5y/z5yjYLSEm5WdJ1+lBeSlO73c6Ykd9wfHcc8R1OU6dJDpW5w/V2Z765+Wrz/JJPT+M+\nFt2t21Tp3QhREpzD0K92a6yNrYjEYx1k0gqeX7q9iIt7heLCQvP8FcmZBHrnJ5vNRmLyAfPnmTOG\nSEC+QEYPyQPoBa2JiCzmhnFv+H6RG8/NRqZ9fJt5rumJlqxf/xDG7nRLgc8YMKCetFmIkm7tMHJF\no/ZsO2qtzFJQOpPRiTO55rmJjz3C6g9l0opnd+igQSOYO3c+O081Iqf2UbrdtJ17r4oMeFfyyYJ8\nAC5v2K5at09V+2X9PcAmajc8TqN239O63koKChqSlrbE56Q/5505wKwdc83jzU+0dj4DSALOmkMQ\nIvTInXOIc094kNLGNev56GlrrI/1tu7z1TXuiTVk0oqTe3dofHw8qal9eb7/U0RHRmNLiAh4V/KR\nU8fMx/e0Gxqwzw1ngwZ1ITHx9RK9Ir95fjIREVGMHXtHpSb9FRQWkPXLV+bPdyffIPtrhxEJziHM\nM/ANv3ueeW7SN5N9zgQNBG/rPk/ZjDSj+UcaBK1cocBut/NR+kpq59fjQP6hgA5X2O12JmW5Nt6o\nVys2YJ8drux2OyNGLCI7+wHgbRY8+jfz3L6i66jsDOt1Bzeaj8d0elz21w4zEpxDmLfA53TyXL4Z\ntO+661NLjut++coo5ArfO9eFV19WzkgAYMM+34krLvQz3fM259cyhh/OHqxt5nQO9sVeKHP9vX4N\nPM3pXxuj598MQI+npxJdJ7/C71VYVMinO11d1i3jWgCyv3Y4keAcZi7Oae36ISIKiGbdugeClgfZ\nc3Zpr9v+Y5778xMb5Qq/DMYX+d3AbA5uNrqU/7tvsd8Co3svzNixtdmyY4h5buWrD5s5nQO91jrc\nfTd7gPl40Pt/JKpWXoUuVr89vJWCogJqRdXkxR7/5+9iiiCQ4BzCvC2reLJ/inm+fqIrkcWGDTsD\nX0BKzy694qHD5jm5wvdlBZBC/mFX7vSPZq4oMc+gqoJkyV6YGnR+9CPzXO4v37nldJY5AufL9ffa\nC+dyp8KzcDC9nfmcIdPGMP3jpHL/JoqKi1jy80oiIyIZ2+VpGtjq+73sIvBktnYI85XwoHm3DRzf\nfREwnU6dgpe/2n126WMrjKQaLepeFLTyhILk5N68/fZrZGcbS2VOHmpI3SbHOFbzIMOGZToCqb/S\nM95Ek44LAMg/Up/ExPVkZ99ehe9fPbn+XpdSUFAHmEN0dDTJycOZt3exObnrufXjeb3PeGpF1fT6\nPm9v/YD9+Qe5tF5LGsbIvI1wJXfOIc5zjMlut6P2XQdAg0s3AZ/RtesZUlJuDm5BocQ2efd3vCeI\nJbE+m83G558/RmLi68BZVr/8MADH6x4ud2OF85Wc3Jtu3d4HzlLvEtfa5usKW/P554/JLOAq4vx7\nHTXqdlJSbgKMXovbW5TM+PZ01p9Z+vMXnDzrGoc+V3SOzF0L+f6YBqDnRV0R4UvunMOIc9xw7dr7\nufmlAzRqt58X/57PiHvvskTXcXauq5t90SxjclMgN3IINfHx8WRlPWD2jByocxG/FB+gZmweZ/Ma\nAHZgMWvXfl8lv8fi4kJgEVff97V57Pd3DyUiIsJnSkpROe5/q2D0gDwzeQgzds0xnzP3p8+Y+9Nn\nXNnoclRCK2b/sBCAmpE1GdC6H92adQpK2UVgyJ1zGHEfNzy0rSORNcBe96Rlvkjf3jLVfCyTiyrG\nvWekS7NrIaKY3sPfBnKBD4EkMjOfveDfY0bGKtavf4iIqCQadzTWyEcWRREREVGqHFb5/xTKvK20\n2LUil7Fdni713K1HvzMDM0D2v9rRI6H8HNr+mJcgAkuCc5g6tNWYZHLCFrxkJJ5fEKfOnQbgwOb2\nyOSiyruuyVVEEME1AyMYPPhNYARV/Xtsdu1283HLQ+3kCz7ALqrblJevH1fm+R2f3chXqx4st629\nJf6R9gs9EpzDiPvs7VPH6gJwKG5v+S/yk1JfEMmu7rqNU5KDUqZQF1+rHm3rt2ZP3j6uvT6xSt/b\n+X9H3b7EPLZqZpR8wftJeVtLxtasy1s3TmRSr78RHWmMPG5JH8Ts5H+zZbrvTG3e7srlAjj0SHAO\nI+7LlsY884V5/MSZvICXxfMLYm/hNea5pnEyueh8dW5i/B4v7m6r0klaNpuN99JvoZHKNo+tWzcS\n+YL3j4pk84osjODqPTdw9U+9aZK7FyhA/maqD5kQFqZs51z7Oa/ct5pBbW4r59n+1+OpD8zH2dmP\nkpg4jpEjryYlRSYXVcbVjTuSsSOTTUe3kZHxGDNnVt0krS0535mPWxy/7ILeS/jmvsTQk91uZ+jQ\nWaxf3xyoQc+ex5gwwbn0qvy2Tk7uTWbmB2bGQCOYD/ZDDYQ/VejOWSnVVSm10svxwUqpb5RSXyul\nHqn64onKKNmV7PqjX7r3i4CXxbPbzmnX572BWLKzXyA6OloCcyXF1IihY8P2HDx1mKPnjlfZJK3i\n4mIyd7kmHT2ZlCLLp4IoPX0Z69fXAm4DbmPNGhtQXKG2lhzb4cHnnbNSagzGBqEnvZx+HbgGyAe+\nV0rN0FqfqNoiiopydSUXArNZM2kkPZ95D4DTp08TExPj9zK4bw85bVp/5s6dT2HEOb51nN8+a6Df\nyxDuOje9hs1HtpG2PJMWOZdVyTKqnTm7S/zsK8GN8K8NG3YBz2IMKwDcx4YNrzBqVMVeX95duQgN\nFblz3gUMASK8nCsA4oEYx/liL88RAbccSGH/xmvNI+9l+n9fV89JYCNGLCI5uTddk4zJS2ePxFCQ\nH43ciV2YNnVaUXQmiuxzOfzp2aQqmaz15ibXNp7juj0DyPKpYOrUqfSwgrdjInz5DM5a6znAuTJO\nvwZsBLYD87XWuVVYNlFJrq7kAscR1/XUj003+P3zS88STWb06P8w50tj040nbhgmXW1V4JNZX7H7\nix7ENDhB867bzN+zc8lTZde4Hso/XOLnJrUTzMeyXjY4hg7tSWzsyziHFeLiJjJ0aM9gF0sE0HlP\nCFNKXQL8HmgJnAI+VErdqbX+pLzXJSSE976wwa1fLCtW3MeUKYuYPPkf7Nw5usTZqihbee8RG+se\nbHOAV8jM7MYt3TZTr3YEnRI7cOMz1k45GAr/P2Njbez4rDutb/kSNWAJ+9ZuJjPzWTIzYf78qUAk\nq1enArBwYTqLFw83L4S81e+xFWPMx6/0fY6E+sZzjElJs8nKSvH6XlYTCm1XUZ9+uoq8vCeBpQDk\n5j7B8uWreeSRpOAWzI/Cqf2qwoXM1rZhDG6e0VoXKaUOY3Rxl+vIkcAv6wmUhIRYS9QvOflmBg26\nnoyM+Zw60I3vmq0DYN+Bo9hq1Drv9/VVv6SkLnTv/gFr1yYD7wAvUKteLvEtF3Nw62W8s3+JpcfB\nrNJ+ntzH8ZOTe5OU1IXp0zP55euraNFtMwkdbuPI98b65NWrGwK/xTlWmZV1L5Mnzyc5uTcLF35N\nXp69xBj19qM/lPisuufizd9BWtoSR2Au+V5WbEOrtt35ysuzY3zFOoPxWfLy7GFVR3fh1n6ezufC\nozLrnIsBlFJ3K6VGaa13YOx79pVSajVQD0irdAmEXzjHCx+9x7Uv77ivXvL7Z86cOdiRvcrYfiY0\nAwAAE0xJREFUYrBxx58AOLzN25QF4Yu3bE8AM2cOpm8rI2iq2zOAvo5/G0u9R0HBWYYNy+R3v+tr\nvkdOTg4fpH3O21tdS9weveqBANRIVITnaoc+fT6UORrVTIWCs9Z6j9a6h+PxDK31FMfjN7TWnbXW\nvbTW92utyxqbFhaQf+6U3z/DZrPRvXsH8+e2txnjzTV+3S9fLuehrGxPNpuN3w8fSp3T9Wh2bS71\nLjnkeM7/0arVq7gvgYKIUu/Rr99bfLSl5PzN1nValfi5vCxWwr88l0NZeThB+IdkCKsGRnZMMR8X\nFRdd8Pv5miTk+lI/Q4PWRvrQT6c+Ll8ufnBx3qUAXH3fbIzOLRupqReXmHgXHR1d6nX7DnRGDVhh\n/nxoe4NSs75lvWxwyWz56k0yhFUD1zS+wnyctn0GkRuN2bjnsz7W21Z3nl/azi/1v6543lwc36B+\n/QurRDXlK9vTw4MG8+hHu2jccQcXdd7A/m+2sGjRGWbP7m+2SXJyb+bMeZ9164xu61pxrzLw3X0l\nPmf137tRXNSPjIySY8qyXlaI4JA75zDmfofrtPHIlgvazKCiSfVtNhsnbTkXWAPh6+7VZrNxTZQC\noOcfplKv5W9Yv35UiTax2WwMGFAPWESz69JLBeaz+VEUF90SkPoIISpGgnOY8pxItH+9KzuYsdGN\nfzczKCwqNB+P6fS4Xz6juvDVvRnrtkii7yuvktBhZ6nnRNaM5LpRJ7j+mZLr3U/ujeS/D/4diKRr\n1ykUFBTImmYLkvXm1Y8E5zBV6g73Xy+Y5zrcubDM17nz9oVQ0UlCo794znzcMq7FBdRE+JKc3Jv9\n77iGLm4Y9y8a9zDyZQP8lLOHn9ts4dKb1pZ67fE5HZkwYSETJnxKcXERY8fGMGZMDYYOnSVBwCLs\ndju33vqxbN9ZzUhwriaKCmqaj9sP+hxfs2/L2rDds5t12rT+ZGSsKnVFX+yWyVWu9v3LZrPxwTsD\n+OljV2avhb8s44+LX2DmD5m88e3bHDtznN+07F7idatffox16x4kOtr4v/H11zE4N1pYv74W6enL\nAlgLUZaMjFVu681l+87qQoJzmPK8w42Lm8iRH1uZ51/4+6xyZ996G1tOSzOWRTm7WZOTezNixKJS\nATz7xF7zfRY9NVau9gNg7tyv+XbeX1jzqmutsr1WPqsOrDUvlFb+7Lpznp38Ngc3dzR/NjZaGIGz\nvY2NFnYFpvCiBOnCFiDBOWx53uFu2DCCW2tfap5vesOFb9dY1uSwVzf+y3zOyQPNkav9wNm/4TDL\nnvtDuc85tTMez2EJ2WjBGrz1WA0a1IU+fdKR9ebViyylCmOey2AeSL2VjSuMta1L937BoDa3lfla\nb0t4UlPvIy+voMzXABRFFJZ7XviHq70ac3x3S+bc9yYN1U/0GftPis5FsOjJcbz24kaGDOnFmR5n\nyKhVcivIlJSbmD/ftdyqW7eppKTcEcwqVUslL3hh7dpU5s41kpBMnizbd1YnEpyrmcviLzX37s07\ne5LYmnW9Pq+s/Xzdg7O3AH7qynw4bpw/OK09xtV+6fW5omo52ys9fRnvvfc62dlPcXhba2YndwaG\n4+wki4iI8Lp22WazMWvWHW7tfYcEAAuR9ebVT4RzRmeAFId7cnOr16+ouIjHVz4LQJ/mPbir7aAK\nv9Zb/dw3ZRg2rBd/XPtX89xrPV4osWGD1b/sQ6H9KsJut5Oevpz33ttMdvZzgI3u3dNYscJ3z0eo\nCqe2M5L8pALGRe3MmYNp0SIhLOpXlnBpv7IkJMRWenMBuXOuZiIjXNMMsn75qlLB2Rv3K/rvjv1o\nHh/Yql9IBeZwYrPZGDUqiZSUm8jIMLYc9NbzIaynrB4rUf1IcK6GJvZ6njGrnwfg6OljNIppeF7v\n47mV4b+3TDXP/fupw+Wm+BT+59kV6pwFDHLBZGXShS1AZmtXS3Wia3Nv+7sA2HBoM1D55Rues0qH\nPzjTPJd4rEOFUnyKwJFEFkKEFgnO1dTVCZdTI7IG3xzcxOnTp70mHCmP5zKqqG5nzHMN85v5t/Ci\n0iSRhRChRYJzNeF5ZxxTI4YrGrbn4KnDfDDnswu6042MPkvjy42EFY1jGnF3ch/ZB9hiCgo8x5rt\nfPnlNkl0IYRFyZhzNVDWNo+dml7DpiPbOFbnYKXf030Z1aU3ugL5qCvuk0ktllQMTAPuA+zAP5k3\nbyzz5smcACGsSIJzNeAtsUFGxnzuSfkNMTVsnG50nK7d3mH9uuYAdO36P5KT7yr3PZ0BeHr6J2xp\nvdo8flHdpuZ5mdRiHUb+7NuBpcA24E94/n+Q9hLCOqRbuxorLCgk9ngjTpzNpaDZASAJSCIiIqrC\n77H0W9ddd+43TaSL1KKSk3vTp89s4BagXbCLI4TwQYJzNeBtm8dBg7owbFgmU59NoegcNOkTSURU\nJBDNunUPMHr0f8jJyTHHqXNycvjPfxaWGKNMT19Ow9/uMT9n+b+fkklGFmWz2Vi8eDgTJ85nwoTT\ndOv2PjInQAjrkm7tasB9DNiYGBTLn/6Uxtq1zwJL+GnZJVx2615a981i16IbAcjMLGTZsjTy8oxs\nYuPHv0Ju7mjARmbmB0yb1p9ps7/lqmeNiUbHd7fg3OkY8zM910DLeGbwuQ81pKTYZU6AEBYm6Tur\nkNVT0JWcGPYZxt69C6kZ243+b4yHiJoseeZZTv86D6gPDMQ5LmncZS3F6Po+y+DBr1Bj2H7zvec+\n8BLNm7xDVpaxcYL7BLTu3UNjwpHV2+9ChXP9wrluIPULdeeTvlO6tauRkhPD+gHvAQc4m7eA7TMH\nULPOaW7/919JvOlSwEcgrX3OfLh9VlsKTn3ByJGXYLPZytxKUgghRMVIcK62bEACMBK4l5+WutbB\ndho1mzs+/Iy45v/BOS4ZFzcR6IUxRvkBNQYcNp//w5zH6N79MCkpNwe0BkIIEa6kW7sKWb1rxnPH\nm8TEcWRnv4iz6zquhabfpH+UeE39vCa0yG3DsAF9WL58K3l5di65wcanuxcAcM0vfahRFF1iXLms\nnXWkWzu4wrl+4Vw3kPqFOtmVSpTLMznIoEGPMWKEaz/my5tnMaRlEnN+Xmi+5njsIY7HHqL+Yegz\nuDMvrXqLb3YXEFUYzdguo2lyY2OfnyMTjoQQonLkzrkKheLVn7dZ1cv2ZpG5a2G5r1v+56dRCStD\n4o64okKx/SojnOsXznUDqV+okwlhotKcy2tSU/uaQfbmS/owvsdzZb5m2XN/4tddl8lELyGE8BPp\n1hZe1bfF03nvzYx7qQ9t+n1J9oqe5B9eg7GUKtrXy4UQQlwAuXMW5bIfr8f2jEHkH44FztKw4YtI\nZikhhPAvCc6iTK60n7nAh8Bgjh17hsTEcUyYMCesxpuFEMJKKhSclVJdlVIrvRzvrJRapZRarZTK\nUErVrPoiimBxzroePPhNYARGd3Ys2dkvEB0dLYFZCCH8xGdwVkqNAaYAtTyORwDvAqla617AciDR\nH4UUwWOz2ejevUOwiyGEENVKRe6cdwFDAM+p4G2BY8DTSqkvgHitta7a4gkr8LarlYw1CyGE/1Ro\nnbNSqhUwQ2vd3e1YT4ydEK4BfgIWAK9orUt1f7sJ6KJqUXXsdjtpacsBSE29Sbq0hRCi4gKaIewY\nsMt5t6yUWgx0AsoLzuG+0Dys6/fII0kcOZJHXl4BeXkFvl8QYsK9/cK5fuFcN5D6hbqEhNhKv+ZC\nZmvvBuoqpVo7fu4FbL+A9xNCCCEElbtzLgZQSt0N1NVaT1FKPQh87JgctkZrvcgfhRRCCCGqkwoF\nZ631HqCH4/EMt+Mrga5+KZkQQghRTUkSEiGEEMJiJDgLIYQQFiMbXwivvG0lKYQQIjAkOItS7HY7\nw4Zlsnbt/QBkZn7AzJmDgcovBxBCCFF50q0tSsnIWOUIzNFAtOzbLIQQASbBWQghhLAYCc6iFMml\nLYQQwSVjzqIU51aRGRnzAUhOln2bhRAikCQ4C69sNhupqX2DXQwhhKiWpFtbCCGEsBi5cxYVZmwb\nuQSQtc9CCOFPEpxFhdjtdoYOnU1WVgrgWvssAVoIIaqedGuLCsnIWOUIzLL2WQgh/E2CsxBCCGEx\nEpxFhSQn96ZPn3Rk7bMQQvifjDmLCrHZbCxePJzJk2XtsxBC+JsEZ1FhsvZZCCECQ7q1hRBCCIuR\n4CyEEEJYjARnIYQQwmIkOAshhBAWI8FZCCGEsBgJzkIIIYTFSHAWQgghLEaCsxBCCGExEpyFEEII\ni5HgLIQQQliMBGchhBDCYiQ4CyGEEBYjwVkIIYSwmAoFZ6VUV6XUynLOv6uUeqnqiiWEEEJUXz6D\ns1JqDDAFqFXG+YeBjkBx1RZNCCGEqJ4qcue8CxgCRHieUEr1ALoA73g7L4QQQojK8xmctdZzgHOe\nx5VSzYBxwO+RwCyEEEJUmRoX8No7gUbAZ0BToLZS6get9fRyXhORkBB7AR9pfVK/0Cb1C13hXDeQ\n+lU35x2ctdaTgckASqkRQDsfgVkIIYQQFVCZpVTFAEqpu5VSo8o6L4QQQogLE1FcLDFVCCGEsBJJ\nQiKEEEJYjARnIYQQwmIkOAshhBAWI8FZCCGEsJgLWedcJqVUJPBv4ErgDDBSa/2T2/mngAeBI45D\nD2utd/ijLP6ilOoKvKy1/o3H8QHAXzASt0zVWr8XjPJdqHLqF9Jtp5SKBqYCLTFS0o7XWs93Ox/S\n7VeB+oV6+0VhpBNui7FC5BGt9Xdu50O9/XzVL6TbD0Ap1RjYCNzkXvZQbzuncupXqbbzS3AGBgE1\ntdY9HF/yrzmOOV0LpGitN/np8/3KkW/8XuCkx/Fo4HWgE3AKWKOUmqe1Phz4Up6/surnENJtB9wD\nHNFapyil6gObgfkQNu1XZv0cQr39bgeKtNbXK6X6ABNwfLeESfuVWT+HkG4/Rxu9A+R7OR7qbVdm\n/Rwq1Xb+6tbuCSwG0Fqvx/iFu7sOeE4ptVop9ayfyuBPZeUbbw/s0lqf0FoXAF8CvQNduCpQZj51\nQr/tZmOknQXj/797atpwaL/y6gch3n5a6/8CDzt+bAUcdzsd8u3no34Q4u0HTALeBg54HA/5tnMo\nq35QybbzV3COA3Ldfi50dHU7zcD4D3gjcL1SKslP5fCLsvKNY9T7hNvPeUC9gBSqCpVTPwj9tsvX\nWp9USsViBLKxbqdDvv181A9CvP0AtNaFSqk04J/Ax26nQr79oNz6QQi3n1IqFaNXZ4njkPvFf8i3\nnY/6QSXbzl/BORdwT5QaqbUucvv5Ta31r44rpIXANX4qR6CdoGS9Yyl95RvqQr7tlFItgBXAdK11\nhtupsGi/cuoHYdB+AFrrVIxx2SlKqRjH4bBoPyizfhDa7Xc/cItSaiVwNTDNMT4L4dF25dUPKtl2\n/hpzXgMMAGYrpboBW50nlFL1gK1KqQ4YYws3Au/7qRyB9iNwmWOsLx+jW2ZScItUdcKh7ZRSTYAl\nwKNa65Uep0O+/cqrX5i0XwrQXGv9EnAaKMKVOjgc2q/M+oV6+2mt+zgfOwLYw25jyiHfduXV73za\nzl/BORPjCmKN4+f7lVJ3A3W11lMc/e0rMWZyL9NaL/ZTOfzNzDeOq25PA59j9Eq8r7X2NvYQKrzV\nL9Tb7jmM7rJxSinn2OwUoE6YtJ+v+oV6+30CpCmlsoBoYDQwWCkVLn9/vuoX6u3nLiKMvzuhdP0q\n1XaSW1sIIYSwGElCIoQQQliMBGchhBDCYiQ4CyGEEBYjwVkIIYSwGH/N1hZCCCGCwrF06UOM9dI1\ngae11us8njMKeAgj4dJ4rfVCx5ryD4EEjEQoI7TWR5VSPTHSUBdjzLT+Szmf3QF41/HjToy9JQor\nWwe5cxZCCBFungKWaq1vAFKBt9xPKqWaAo8DPYB+wEtKqZrA74AtWuvewHTgz46XvAakaq27Azco\npTqW89kTgGe11tc7fh5wPhWQO2chhBAhTyn1PHBAa/0O8AbGemIw1ouf9nh6F2CNI1tXgVJqF8Yu\nij2BVxzPWYyxSxaO1zd0BHAbjvTGSqmXgOuBKOB1rfUnwB1a6yLHc5sCOedTHwnOQgghQpZS6i6M\nO95WwFmlVDIwQWu9zHGHnI6RzMVdLN5zebvvC+Ge3/tVYAFwDNgCaKVUf6CV1rqXUsoGrFVKLdVa\nn1BKXQIswwjMWzkP0q0thBAiZGmtZzn2nU8DXtNa/8YRmK/ACJD/p7Ve7fEyz/0fYjECqfvxWCDH\nEXj/CbTXWrfB2LXvD0BH4DpHqs5FGDe7LR1l2qu1bouxfeTr51MvCc5CCCHCimNS1mzgbq31516e\n8jXQSylVyzF5rD2wHWNfiNscz+kPrMKIk9EYObEBDgLxGPnAVzouDG5xfN5updQ8pVQbx3NPApWe\nDAaSvlMIIUSYUUrNxRhD/tlxKEdrPVgp9RTGvtHzlVIjMWZrR2J0g2c6ZmtPA5phjFkP11ofduTI\nfgojQB/HmBx2Qin1GtAZqAvM0VqPV0p1x9i04yzGJh4jtdaHKlsHCc5CCCGExUi3thBCCGExEpyF\nEEIIi5HgLIQQQliMBGchhBDCYiQ4CyGEEBYjwVkIIYSwGAnOQgghhMX8Pxko77afsw3LAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efc7cadac18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set_context(\"notebook\")\n",
    "np.random.seed(seed=0)\n",
    "\n",
    "intercept = 2.\n",
    "trend = 0.1\n",
    "variability = 0.4\n",
    "noise = 0.05\n",
    "\n",
    "x = np.arange(2009, 2012, 0.01)\n",
    "x += np.random.random(size=x.size) * 0.1 \n",
    "y = variability * np.cos(x * np.pi / 0.5) + trend * (x - x.min()) + intercept\n",
    "y += np.random.normal(scale=noise, size=y.size)\n",
    "w = np.ones_like(x)\n",
    "\n",
    "res = fit_driftfourier(x, y, w, 2)\n",
    "\n",
    "print(\"fitted trend (slope): {}\".format(res[1]))\n",
    "print(\"fitted intercept at xmin: {}\".format(res[0]))\n",
    "print(\"fitted fourier parameters: {}\".format(res[2]))\n",
    "print(\"std of residuals: {}\".format(res[6]))\n",
    "\n",
    "f_drift, f_fourier, f_driftfourier = res[3:6]\n",
    "\n",
    "plt.scatter(x, y, label='data')\n",
    "plt.plot(x, f_drift(x), label='fitted annual trend')\n",
    "plt.plot(x, f_driftfourier(x),\n",
    "         label='fitted annual trend + intra-annual variability')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statistical significance using bootstrap resampling: 2.5%, 50% and 97.5% percentiles for each fitted parameter and histogram of the distribution obtained for the annual trend (slope) and intercept estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'intercept': array([ 1.98152376,  1.99281233,  2.00459845]),\n",
      " 'pfourier': array([[ 0.36480368,  0.37266487,  0.38059166],\n",
      "       [-0.14720898, -0.13925057, -0.13075194],\n",
      "       [-0.00928356, -0.00159154,  0.00640595],\n",
      "       [-0.00456552,  0.00319257,  0.01144754]]),\n",
      " 'slope': array([ 0.09453977,  0.10123101,  0.10780417])}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7efc7c9934a8>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAADKCAYAAACBt8VKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH2ZJREFUeJzt3X2QJHd52PHvjXRzuzs7tyDdSpRjCqeM9Sthl3gRLykZ\nI2HejN9E8Q/lYMqQshRilQyOCUkORyQpyVBxLJeFgSQ6sOLgxIkUhA0qBDJQpxfLErIwiQrzoAOM\nghHyHXe6nZ3dvT3dbP6Ymb3Zudmd2Znuedn9fqpUuvlNd+/TPT3PPN3961/vWVtbQ5IkSf0pjDoA\nSZKkSWYxJUmSNACLKUmSpAFYTEmSJA3AYkqSJGkAFlOSJEkDOH+rN1NK5wG3ApcAa8A7gSLwGeAb\njck+EhG3p5SuAa4FngFujIi7cotaknpgDpM0DHu2GmcqpXQ18AsR8asppSuB3wA+DcxFxM0t0z0H\n+DxwOTAN3A+8NCJW8wxekrZiDpM0DFuemYqIP00pfabx8keAp6knm9RIUo8D7wZeDjwQEaeB0yml\nI8BlwCN5BS5J3ZjDJA1D1z5TEXEmpXQb8PvAHwMPA++JiCuBbwHvB8rAyZbZKsBc5tFK0jaZwyTl\nbcszU00R8faU0sXAQ8AVEfG9xlt3Ah8C7qWejJrKwImtlrm2tra2Z8+e7UcsaZKN5EtvDpOUgU2/\n8N06oL8N+OGI+ACwDNSAT6aUro+ILwOvpX4a/GHgppTSPmAKuBR4bMuI9uzh6NHKttZiFObnyxMR\nJ0xOrMaZvUmJdX6+3H2iDJnD+jcp+1S/dvr6wc5fx2Gv31b5q9uZqTuA21JKh4G9wLuAJ4APp5RO\nA08C10bEYkrpFuA+6pcOD9pxU9IYMIdJyl23DujLwFs6vPXKDtMeAg5lFJckDcwcJmkYHLRTkiRp\nABZTkiRJA7CYkiRJGkBPQyNI46pWq7G4eO7dHLOzZQoFjxUkSfmzmNJEW1yscM9DR5ieKa23LS9V\ned0rns/+/Y65KEnKn8WUJt70TImZ0nDHL5IkqcnrIJIkSQOwmJIkSRqAxZQkSdIA7DMlSdr1vDNY\ng7CYkiTtet4ZrEFYTEmShHcGq3+eu5QkSRqAZ6Y0Fjr1V6jVagDn9FewD4MkaZxYTGksdOqvcPzY\nUxQK5/OsCy5cb7MPgyRp3FhMaSTaz0RVKgtMTc1s6K+wVF2kUDjPPgySxoZ3/amTLYuplNJ5wK3A\nJcAa8E7gFHAbUAMeA66LiLWU0jXAtcAzwI0RcVeOcWuCtCefYrHGk0/+HQ8+9hQzs7NA/SzUTGk/\npfL+UYWpHcgcpqx515866XZm6ueBWkS8MqV0JfDbjfaDEXFvSumjwNUppb8ErgcuB6aB+1NK90TE\nam6Ra2K0J5/Z0nGe+M4TzJT2r591WqoujjJE7VzmMGXOu/7UbstzkhHxp8A/bbz8EeAEcHlE3Nto\n+yzwWuBlwAMRcToiFoAjwGW5RKyJ1Ew+M6Uypdn9TE2Xus8kDcgcJmkYul7gjYgzKaXbgN8H/hjY\n0/J2BZgD9gMnO7RL0kiZwyTlracO6BHx9pTSxcDDwFTLW/uBp4EFoPWcZ5n6EeCW5ucn4zTppMQJ\n4xlrsVhjtnSc0uzZXadUKlIo7KXcaFuubny9WVuBVQ4cKDM3V9502e3TDGIct+dmJinWYdvtOaxf\nu2n9es0lnaajtsLevTWKxdqG5ZfLo++Uvps+w1Hq1gH9bcAPR8QHgGXgDPBISunKiDgMvBH4AvUE\ndVNKaR/1RHUp9Y6dWzp69Nw7IsbN/Hx5IuKE8Y11YaHCYvUUNVYAKM9OUa2uUiicYd90va399WZt\nS9VTHDtWYXW10HHZnabp17huz04mJdZhJz5zWP8mZZ/qV/v6dcoli5Vlvv3tv6Nc3njncaWyQo3i\netuxYye488+/P3bDuOy2z3AYf28z3c5M3QHcllI6DOwF3gV8Hbg1pVQEvgbc0bgT5hbgPuqXDg/a\ncVPSGDCHqW8ry0scfvTEhiJpszuPp6Y3Du1Sq9WoVBbOWaZDKOxMWxZTEbEMvKXDW1d1mPYQcCib\nsCRpcOYwDaq9SOr1zuNOhdg4nK1SPhy0U5KkHLQXYtq5PNcoSZI0AIspSZKkAXiZTxOlvVNnpbLA\nWm1thBFJknY7iylNlPZOnT7TT5I0ahZTmjitnTp9pp+kSeFwCTuXxZQkSUPgcAk7l8WUJGnXqNVq\nLC5WKBZrLCxsHNl8GP0vHS5hZ7KYkiTtGouLFe556Ajz8xewWD213m7/Sw3CYkqStKtMz5Qoze5v\ne6an/S/VP3u8SZIkDcAzU5KkHanZP6qVY9MpDxZTkqQdqdk/anqmtN7W7BslZcliSpK0Y03PlDbc\nPWffKOXBPlOSJEkDsJiSJEkawJaX+VJKe4GPA88D9gE3At8FPgN8ozHZRyLi9pTSNcC1wDPAjRFx\nV25Ra2x06uAJo308go9sUJM5TNIwdOsz9VbgaES8LaX0bOCrwL8Dfjcibm5OlFJ6DnA9cDkwDdyf\nUronIlZziltjolMHz1E/HsFHNqiFOUxS7roVU7cDdzT+XQBOU082KaV0NfA48G7g5cADEXEaOJ1S\nOgJcBjySS9QaK+0dPNvPDI3iVmQf2aAGc5ik3G1ZTEVEFSClVKaelN4HTAG3RsRXUkoHgfcDfw2c\nbJm1AngKYJdqPzPkYxo0KuYwScPQdWiElNJzgU8CH46IP0kpzUVEM+ncCXwIuBdoPQ1QBk50W/b8\n/GScOZiUOGH4sRaLNWZLxynNTq23LVeLlMslDszPA7CHVQqFvZRbpimVihvalqvFc6bppa3X+Qqs\ncuBAmbm57W0fP/vJZw7r3yStX/2M+Mb+m3v31ijNFJntkB+AbeebPNr6zU29mqTPsB/jsn7dOqBf\nDHwe+LWI+FKj+e6U0q9HxJeB11I/Df4wcFNKaR/1o75Lgce6/fGjR8/tuDxu5ufLExEnjCbWhYUK\ni9VTG55xVa2uUiicYd/0SsfX5dmprtP02tbrfEvVUxw7VmF1tfcO6H722Rt24jOH9W9S9qmmhYWT\nmw7QubZn33pbMz8cmIfK4vbyTR5t/eSmXk3aZ7hdw16/rfJXtzNTB6mf6r4hpXRDo+3dwO+llE4D\nTwLXRsRiSukW4D7q/RIO2nFT0hgwh+0iDtCpUenWZ+pdwLs6vPXKDtMeAg5lFJeUqU7DJThUws5n\nDpM0DD5ORrtCe6d4h0qQJGXFYkq7hsMlSJLyYDElSdKIbPbEBrArwiSxmJIkaUQ6PbEB7IowaSym\n1LNOz+EbxejmkrST2AVh8llMqWednsPn6OaSpN3OYkrb4jgukiRtZM82SZKkAVhMSZIkDcBiSpIk\naQAWU5IkSQOwmJIkSRqAxZQkSdIAHBpBu9Jmj3Dw8Q2SpO2ymNKu1OkRDj6+QZLUD4sp7Vo+wkGS\nlIUti6mU0l7g48DzgH3AjcDfALcBNeAx4LqIWEspXQNcCzwD3BgRd+UYtyR1ZQ6TNAzdOoe8FTga\nEa8Cfgb4MPC7wMFG2x7g6pTSc4DrgSuANwAfSCkV8wtbknpiDpOUu26X+W4H7mj8uwCcBl4SEfc2\n2j4LvB44AzwQEaeB0ymlI8BlwCPZhyxJPTOHScrdlsVURFQBUkpl6knpt4D/2DJJBZgD9gMnO7RL\n0siYwyQNQ9cO6Cml5wKfBD4cEf8jpfQfWt7eDzwNLACtPXnLwIluy56fn4zOv5MSJ+Qba7FYY7Z0\nnNLs1HrbcrVIobCX8hZtnaYplbpP0++y+52vwCoHDpSZmzu7Df3sJ585rH+TtH7bzU/AtvPGsNoA\nqK2wd2+NYrF2tqlW/3f78C3l8uZDukzSZ9iPcVm/bh3QLwY+D/xaRHyp0fyVlNKVEXEYeCPwBeBh\n4KaU0j5gCriUesfOLR09Whkk9qGYny9PRJyQf6wLCxUWq6eosbLeVq2uUiicYd/05m3tr8uzU12n\n6XfZg8y3VD3FsWMVVlfrScnPPnvDTnzmsP5Nyj7VtN38dGAeKovbyxvDagM4duwEd/759zcM33L8\n2FMUCuf3PKTLpH2G2zXs9dsqf3U7M3WQ+qnuG1JKNzTa3gXc0uic+TXgjsadMLcA91Hvl3AwIlYH\njlySBmMO24FqtRqLixt/RCuVBdZqayOKKB/tw7csVRcpFM5zSJcx1K3P1LuoJ552V3WY9hBwKJuw\nJGlw5rDJt1nh9OBjTzEzO7vedvzYU8yU9lMq7x92iJKDdkqSxtfiYoV7HjrC9Expva1ZOLWftZFG\nxWJKkjTWpmdKFk4aaz7RVZIkaQAWU5IkSQPwMp8kSROqVqtRqSyc0z476x1/w2QxJUnShFpZXuLw\noyc6jj118cUO4j8sFlOSJE2w9vGoNHwWU1JD++nyYrFGrXbuoxskSWplMSU1tJ8uL+z5f1zxE8/t\n+JgGSZKaLKa0qfaRh3fi4xratZ4uL+DTRCRJ3VlMaVPtIw/7uAZJks5lMaUttY487KjDkiSdy561\nkiRJA7CYkiRJGoCX+QSc29kcdkeHc0mSBmUxJeDczuZgh3NJknrRUzGVUnoF8MGIeHVK6cXAp4HH\nG29/JCJuTyldA1wLPAPcGBF35RKxctPa2RzscK6dwfwlKW9di6mU0nuBXwaav6yXAzdHxM0t0zwH\nuL7x3jRwf0rpnohwoB5JI2P+kjQMvZyZOgK8GfhvjdeXA5eklK6mfnT3buDlwAMRcRo4nVI6AlwG\nPJJ9yJLUM/OXpNx1vZsvIj5J/dR300PAeyLiSuBbwPuBMnCyZZoK4DM4NNGaz+pbWDi54b9arTbq\n0NQj85ekYeinA/qdEdFMPHcCHwLupZ6QmsrAiW4Lmp+fjKdcT0qc0H+sxWKN2dJxSrNT623L1SKF\nwl7Kjbb21722dZqmVMpv2VnN9/dP/YAvxwkuuODsb/HSUpVfvOoFzM2N3z4xSfvpCGWWv2Dnb/Nx\nWL9eclM/bcBA8+fZNuj8BVY5cKD+2Y3DZ5incVm/foqpu1NKvx4RXwZeS/1U+MPATSmlfcAUcCnw\nWLcFHT1a6TbJyM3PlyciThgs1oWFCovVU9RYWW+rVlcpFM6wb3ql4+te29pfl2enclt21vPVaudR\no3j29dopjh2rsLo6XkO0Tcp+OgaJL7P8BZORw/o1LvtUL7mpn7YD81BZzHaZWbUNOv9StZ6n5ubm\nxuIzzMuw99Gt8td2iqnmgEPvBD6cUjoNPAlcGxGLKaVbgPuoXzo8aOdNSWPE/CUpNz0VUxHxt8AV\njX9/FXhlh2kOAYeyDE6SBmX+kpS38bpWIUmSNGEcAV2SNBZ8rJUmlcWUJGks+FgrTSqLKWkbmmNP\ntZudLVMoeNVcGpSPtdIkspiStmFleYnDj57gWRdcuN62vFTlda94Pvv3O86jpNFrHvSdPHmShYWz\nl0096MuPxZS0TVPTMxuOnCVpnDQP+r559BkWq6cAD/ryZjElSdIOMzU9Q2l2/4bBTpUfi6ldyDtm\nJEnKjsXULuQdM5IkZcdiapfyjhlJkrJhMSVJGjq7G2gnsZiSJA2d3Q2GyzHy8mUxJUkaCbsbDI9j\n5OXLYkqSpF3AMfLyYzG1C7T3TbBfgiRJ2bGY2gXa+ybYL0GSpOz0VEyllF4BfDAiXp1Sej5wG1AD\nHgOui4i1lNI1wLXAM8CNEXFXTjGrD619E+yXoN3E/CUpb1278KeU3gvcCuxrNN0MHIyIVwF7gKtT\nSs8BrgeuAN4AfCClVMwnZEnqjflL0jD0cj/kEeDN1BMPwEsi4t7Gvz8LvBZ4GfBARJyOiIXGPJdl\nHawkbZP5S1LuuhZTEfFJ6qe+m/a0/LsCzAH7gZMd2iVpZMxfkoahnw7otZZ/7weeBhaA1vsty8CJ\nbguan5+MWzQnJU6ACy8sUalsHFV4794apZkis7NTACxXixQKeyk3XvfaluV8pVJ+y85qvuUqlEpT\nXZddYJUDB8rMzY12P5mk/XSEMstfsPO3eZ7rVyzWmC0dp7TN725WbcBQ/k4/bXms4zjnrkGMy3ew\nn2LqKymlKyPiMPBG4AvAw8BNKaV9wBRwKfXOnVs6erTSbZKRm58vT0ScUI/129/+3qajCq/tqXcb\nqVZXKRTOsG96ZX2aXtqymq88O5XbsrOfb6XrNEvVUxw7VmF1dXSjCE/KfjoGiS+z/AWTkcP6ldU+\n1emxMVAfoqVSWaHG2e5p/X6f+2k7MA+Vxfz/Tj9teazjuOauQQw7722Vv7ZTTDUHJvpN4NZGB82v\nAXc07oa5BbiP+qXDgxGx2me8GpCjCkvnMH+NSKfHxoBDtGhn6amYioi/pX6nCxHxOHBVh2kOAYcy\njE2aCJ2eeeXzrsaH+Wv02g/wwIM87SwO2ikNqP2ZVz7vSpJ2F4spKQM+80qSdi+vQ0iSJA3AYkqS\nJGkAXuaTMtapQzrYKV2SdiqLKSlj7R3SwU7p2h06jSlVqSywVlvbZA6Nkgd+2bGYknJgh3TtRp3G\nlHI8qfHlgV92LKYkSZlx0ODJ4oFfNjyPJ0mSNACLKUmSpAFYTEmSJA3AYkqSJGkAdkCXJEmb6jTk\nRZPDKNRZTEmSpE11GvICHEahlcXUhGs9YigWaw6QN6YcHE/SJGsf8kIbWUxNuNYjhtnScZ74zhMO\nkDeGHBxPO42jne9MnQ78/Fy767uYSik9CpxsvPwW8AHgNqAGPAZcFxFu/SFoHjGUZqeYmi51n0Ej\n4eB448P8NThHO9+ZOh34+bl211cxlVKaAoiIV7e0/RlwMCLuTSl9FLga+FQmUUpSRsxf2XG0852p\n/cDPz7W7fs9MvRCYSSl9rrGM9wEviYh7G+9/Fng9JiNJ48f8JSlT/RZTVeB3IuJjKaUfA+5ue38R\nsCOIpHFk/pIy4I01Z/VbTH0DOAIQEY+nlH4AvLjl/TLwdLeFzM9PRv+RcY6zWKwxWzpOaXYKgFKp\nSKGwl3LjNcBydWNb++te27Kcrz3OcYjp3GmgVJrKLaYCqxw4UGZuLpv9a5z30zGTSf6Cnb/Nt1q/\n9twDve/725k2rzZgZH97WNundR3zifEHfDme4oILnllvW1qq8otXvSCzvNbNuHwH+y2m3gFcBlyX\nUvoh6snn8ymlKyPiMPBG4AvdFnL0aOdBwMbJ/Hx5rONcWKiwWD1FjRXKs1NUq6sUCmfYN72yPk17\nWy/T5DlfpzhHHdPmbSu5LXupeopjxyqsrg5+BDfu+2nTmCS+TPIXTEYO61e3fao19zT1uu9vZ9q8\n2g7MQ2VxNH97WNundR3zi/E8ahTX22pr2eW1boad97bKX/0WUx8D/jCl1Oxj8A7gB8CtKaUi8DXg\njj6XLUl5Mn9JylRfxVREPAO8rcNbVw0UjSTlzPwlKWu7q4eYJElSxhwBfYI44rCkUTD3SFuzmJog\njjgsaRTMPdLWLKbG1GZHglNTjky7U3Qao2U3js+iyeBo59LmLKbGlEeCO1/7M7B88LEkTSaLqTHm\nkeDO58OPNUqdzoBD/QyppN5ZTI2J9qRm587dx0czaNg6nQGvLla44ieew9RUfWBOMB+pd7s1j1lM\njYn2pOYlvd2n/bIfeOlP+et0Bvzwo0/wzaPPsFg9BZiP1LvdmscspsZIa1Lzkt7u5GU/jYOp6RlK\ns/vXHxVjPtJ2tOex3XC2ymJKkiTlZjecrbKYkiRJudrpZ90tpqQJstndVxdeWOowtXSWo5hL+bGY\nGgGTmnrV3tegUlngwceeYmZ2dr1teanKLx0o46M2tRXHrtM42Wn9qCymRsCkpl619zVo7ic7+XS5\n8uPYdRoXO60flcVUznwsjAbV2tfA/UTSTrGT+lFZTGWs0+Cb7ZdlPAulLNVqNU6ePMnp04UNbcD6\n6fL2102Tekp9N9pqtHI/Q+0EnS79bZW7xkmmxVRKqQB8BLgMOAX8akR8M8u/MU42O+vUWjx1uizj\n2QVlaWV5ic89+E2K+zYW7IXC+RsuD7a+hrMjXZfbivrd+uM87vlrq9HKWz/DzX587Jepcdfp0t9W\nuat1lP6mUeWvrM9MvQkoRsQVKaVXAL/baJs4zUKpWKxt+LBaP6it+j55WUbDND1dYt/0xoK9UDhv\nw37Y+rrZdvjRJ3ZMn4UMjE3+OnXqFPc//H/Zu29qva1aXWTlFFw4v/Vn2OnHp9nuGXGNu/ZLf1vl\nrtZR+qFz/hrWGd2si6mfBO4GiIiHUkovzXj5XT311N+zsFhdf33mzBmePVdmevpsUurlEkjzDNNF\nF1+4/mG1HwXa90mTrlufhV12aWnk+aupVquxXJvivOKz19vWTp/P6tPHzpm2lx+fZru0U7SP0g+d\nLxNudgd01geNWRdT+4HWNTmTUipERC3jv7OpI9/5HieXz57Kri4ucPT4IhddfNF624njxygUzmPu\nWc/u+LrZViptPII7tbLE3Q/8zYb5SqX97CnsWZ9mZblKoXA+S9VKx9e9tvUzX4HV3Jad5Xyd4hx1\nTB2nWVpiZeXMeMW0yTTnnw9naoPvh8tL1XOGYvjSX32bqanps/OtLPOLV71wJ569Gnn+atqzZw+n\nV55midPrbUvVRVZXVra9X2132mZbdXGBpcaBZL/7dJbxZN3Wun7jEE8e22e3fYZPHz/K3U9+t+Nv\n+Qxni6k8ZF1MLQCth0NbJaI98/PZdyB708/9VObLnCyXjTqAHk1KnHrRi14w6hCGZTv5C3LKYXVl\nrv3lN+a0bElZy/o8/QPAzwKklP4R8H8yXr4k5cX8JakvWZ+ZuhN4XUrpgcbrd2S8fEnKi/lLUl/2\nrK15q6wkSVK/dtztOJIkScNkMSVJkjQAiylJkqQBZNYBvdujGFJKvwD8G+AZ4OMRcagxzyHgEqAG\nXBMRkVJ6MfBp4PHG7B+NiP81wjiLjTifD5wGfj0ivppSej5wWyP2x4DrIiKzTmgZxzqybdqYZga4\nB/gnjc+44zx5btOM48xte/Yba6NtnvpdaT8REasppWngE8A8UAF+JSLOHflx9HHuAb4LfKMx24MR\ncTCrOPvRGAX9gxHx6rb2XwL+BbAC3B4Rv5dSejvwK41JpoEXAhdHxMYRBMfINtevY64edszbsc31\n65g3hx1zr1JKe4GPA88D9gE3RsSnW97f7Pd2bB+X1K6fdWx5r+Nnn6csz0ytP4oB+FfUH8UArG+U\nm4HXAVcC16aULgJeD5Qi4pXAvwduasxyOXBzRLy68V9mP1J9xnkNsNSY5xrqHzCNaQ9GxKuAPcDV\nGcaZdawj2aaNWF8K3Av8Q2Ctyzx5btMs48xze/YVa0rpDcDngYtaJv1nwFcb2/OPgN8a0zh/FPir\nlu056kLqvcCt1JN4a/uFwG8DP019xPSrU0ovjojbmrEDjwDXj3khta31Y/NcPZb6WL/N8ua4eitw\ntPG9/hngD5pvbPHb8CZg32bf1THUzzpu+tnnLctiasOjGIDWRzFcChyJiJMRcRq4H3gVsAzMNY5K\n54DVxvSXAz+XUjqcUjqUUspy6NJ+4nxByzzfAP5BSmkOeElE3NuY97PAazOMM+tYR7VNAYrUv8jR\nwzx5btMs48xze/Yb6xngNcCJTstp/H+Y++h24ryc+r76xZTSXSmlSzKOc7uOAG+mXtC3+lHqxenT\njTOmf0n9ewesF48/3nqUPKa2u36b5epxtd3165Q3x/khhrcDNzT+XaB+dqZps9+Gn6SeUzf7ro6b\nftYRNv/sc5VlMdXxUQwt751sea9C/Qt5PzAFfB34L8AtjfcfAt4TEVcC3wLeP+I4/xr4eVgfzG8e\nKLHxw1psTJulLGMd1TYlIv4iIr7bwzznke82zTLOPLdnX7FGxJ9HxPEOy2nuJ819ZBzj/B7w2xHx\n09TPHHwi4zi3JSI+ycbk3fQ48OMppYsaly9fA8y0vH8Q+Lf5RziYPtavNVf/Z+BDw4q1H9tcvxKb\n582xFBHViFhMKZWpFx3va3l7s9+GLb+r46bPddzqs89Vlhtyq0cxnGx7rww8DfxL4IGISNT7GPxR\n49r1nRHxlca0nwJePKI491M/ev44sJBSuo+zR9nHqfcdaF+nLGUR6zeAHzC6bbqdec6Q7zbNMs5P\n5bg9N/u7/TwjboH6vgHD30e34xHgzwAi4gHghzKILXMRcQL4DeB/A/8deBQ4BpBSehZwSUQcHl2E\ng9lk/X7Axlz9IuC/NnL1RNlk/Y7SOW+2F/xjJaX0XOCLwB9FxJ+0vLXZ721W39Wh2eY6tp7pHros\ni6mtHsXwdeDHUkrPbnwBXwU8SL3yb1bKJ6h3iD8fuDul9LJG+2uoJ9pRxPlTjThfDnwxIn4KuAP4\nfkSsAF9JKV3ZmPeN1PuGZCmLWJ+MiFOMbptud548t2mWcX42x+251d/tezkMfx/djhuAdzeW80Lg\niUyiy1hK6XzgpY3v11uoHwR+ofH2q1r+PZG2WL/2XL0XOG8kQQ5gi/XbLG+OpZTSxdT7Hb43Im5r\ne7vT7+1fMGGPS+pjHR8ccogbZPk4mXMexdC4a2I2Im5NKf1z4HPUC7iPRcT3Ukq/A/xh42hgL/WO\nx0sppXcCH04pnQaeBK4dYZxPppROAf8zpXSQ+h0g1zTm/U3g1saH+TXqX8IsZRnryLZpr/M0/p/n\nNs0yzjy3Z7+xNrXe/fhR6mcR7qN+B88/HtM4Pwh8IqX0s9RP0b894zj71ew03/q9O5NS+ivqfb/+\nU0R8qzHtJcDY3h21iV7W75sdcvW/jojl0YXds54+v5TS03TOm+PqIPXLWjeklJr9im6lfpPAZr8N\nk/a4pG2vY9v8Q328i4+TkSRJGsDYdj6TJEmaBBZTkiRJA7CYkiRJGoDFlCRJ0gAspiRJkgZgMSVJ\nkjQAiylJkqQBWExJkiQN4P8DUWisf2GuXykAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efc96e6c630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cf, intercept, slope, pfourier = cf_driftfourier(x, y, w, 2,\n",
    "                                                 nboot=5000,\n",
    "                                                 percentiles=(2.5, 50., 97.5))\n",
    "\n",
    "pprint.pprint(cf)\n",
    "\n",
    "f, axes = plt.subplots(1, 2, figsize=(10, 3))\n",
    "sns.distributions.distplot(slope, kde=False, ax=axes[0])\n",
    "sns.distributions.distplot(intercept, kde=False, ax=axes[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Anaconda (py34)",
   "language": "python",
   "name": "anaconda_py34"
  },
  "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}