cvr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "interstate-fitting",
   "metadata": {},
   "source": [
    "## Theory\n",
    "Define a simple linear model for the MRI signal $s$:\n",
    "* define a signal $s_0$ at $t=0$; assume that $\\newcommand\\etco{\\mathit{ETCO2}} \\newcommand\\cvr{\\mathit{CVR}}\\etco=\\etco_0$ for long enough that a steady-state signal has been reached.\n",
    "* $\\cvr$ is the change in $s$ (as a fraction of $s_0$) per unit change in $ETCO2$\n",
    "* there is a drift $d$ that causes a change in $s$ (as a fraction of $s_0$) per unit time\n",
    "* there is a delay $\\delta$ such that the signal depends on $\\etco(t-\\delta)$\n",
    "\n",
    "Therefore, $ s = s_0 + \\cvr\\cdot\\left(\\etco(t-\\delta)-\\etco_0\\right)\\cdot s_0 + d\\cdot t\\cdot s_0$.  \n",
    "Rearranging, $$ s = s_0\\cdot(1-\\cvr\\cdot\\etco_0) + \\cvr\\cdot s_0\\cdot\\etco(t-\\delta) + d\\cdot s_0\\cdot t$$\n",
    "Using multiple linear regression we obtain an intercept and a coefficient for each regressor:\n",
    "* $\\beta_0 = s_0 (1-\\cvr\\,\\etco_0)$\n",
    "* $\\beta_\\etco = \\cvr\\,s_0$\n",
    "* $\\beta_t = d\\,s_0$\n",
    "\n",
    "So, $\\cvr=\\frac{\\beta_\\etco}{s_0}$. We've defined $s_0$ as the signal at $t=0$, however for noisy data it may be more reliable to take the average signal over multiple early time points."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "liquid-valentine",
   "metadata": {},
   "source": [
    "## Define Python functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "lightweight-mouse",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipywidgets as widgets\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "def delay_etco2(etco2, delay):\n",
    "    # apply a delay to etco2 time course\n",
    "    # etco2 = ETCO2 array, delay = delay (points)\n",
    "    # returns: shifted ETCO2 array\n",
    "\n",
    "    # shift etco2\n",
    "    etco2_shifted = np.roll(etco2, delay)\n",
    "    \n",
    "    # above step is circular shift - correct start/end values (simplistic)\n",
    "    if delay > 0:\n",
    "        etco2_shifted[:delay] = etco2[0]\n",
    "    elif delay < 0:\n",
    "        etco2_shifted[delay::] = etco2[-1]\n",
    "    \n",
    "    return etco2_shifted\n",
    "\n",
    "\n",
    "def gen_cvr_signal(etco2, s0, cvr, delay = 0, sigma = 0, drift = 0, n = 1):\n",
    "    # simulate multiple CVR signals\n",
    "    # etco2 = ETCO2 array, s0 = signal at time zero, cvr = CVR (%/mmHg), delay (points), sigma = std of Gaussian noise, drift = signal drift / s0 per unit time, n = number of signals to generate\n",
    "    # returns a 2D array; one signal per column\n",
    "\n",
    "    # shift etco2\n",
    "    etco2_shifted = delay_etco2(etco2, delay)\n",
    "    \n",
    "    # generate the pure signals\n",
    "    s = np.tile( s0 * (1 + (cvr/100)*(etco2_shifted - etco2[0])), [n, 1] )\n",
    "    \n",
    "    # add drift\n",
    "    t = np.arange(etco2.size)\n",
    "    \n",
    "    s += (drift/100.) * t * s0\n",
    "    \n",
    "    # add noise\n",
    "    s += np.random.normal(loc = 0., scale = sigma, size = s.shape)\n",
    "    \n",
    "    s = s.transpose()\n",
    "    \n",
    "    return s # each signal is a column\n",
    "\n",
    "\n",
    "def fit_cvr(etco2, s, delay = 0, baseline = [0]):\n",
    "    # fit CVR signal(s) with a specified delay\n",
    "    # etco2 = ETCO2 array, s = array of signals (one column per signal), delay = delay (points), baseline = list of indices to use to estimate s0\n",
    "    # returns dict of results {cvr, drift, residuals}; each entry is an array contain 1 value per signal analysed\n",
    "   \n",
    "    # shift etco2\n",
    "    etco2_shifted = delay_etco2(etco2, delay)\n",
    "    \n",
    "    # organise regressors into matrix\n",
    "    a1 = np.ones(etco2_shifted.shape).reshape(-1,1)\n",
    "    a2 = etco2_shifted.reshape(-1,1)\n",
    "    a3 = np.arange(etco2.size).reshape(-1,1)\n",
    "    regs = np.concatenate((a1,a2,a3),1) # matrix with n rows and 2 columns (ones and shifted etco2)\n",
    "    \n",
    "    # do multiple linear regression\n",
    "    x, residuals, rank, svs = np.linalg.lstsq(regs, s, rcond=None)\n",
    "    s0 = np.mean(s[baseline], axis=0)\n",
    "    cvr = (100. * x[1]) / s0\n",
    "    drift = (100. * x[2]) / s0\n",
    "    \n",
    "    results = {'cvr': cvr, 'drift': drift, 'res': residuals}\n",
    "    \n",
    "    return results\n",
    "\n",
    "\n",
    "def fit_cvr_delay(etco2, s, min_delay=0, max_delay=0, baseline = [0]):\n",
    "    # fit CVR signal(s) with optimised delay\n",
    "    # etco2 = ETCO2 array, s = array of signals (one column per signal), min/max_delay = min/max allowed delay (points), baseline = list of indices to use to estimate s0\n",
    "    # returns dict of results {cvr, delay, drift, residuals}; each entry is an array contain 1 value per signal analysed\n",
    "    \n",
    "    cvr, delay, drift, res = [], [], [], []\n",
    "    for signal in s.transpose():\n",
    "        delay_range = range(min_delay, max_delay + 1)\n",
    "        \n",
    "        results_all_shifts = [ fit_cvr(etco2, signal, this_delay, baseline = baseline) for this_delay in delay_range ]\n",
    "\n",
    "        residuals = [ result['res'][0] for result in results_all_shifts ]\n",
    "        res_best = min(residuals)\n",
    "        idx_best = residuals.index(res_best)\n",
    "        \n",
    "        cvr.append(results_all_shifts[idx_best]['cvr'])        \n",
    "        delay.append(delay_range[idx_best])\n",
    "        drift.append(results_all_shifts[idx_best]['drift'])\n",
    "        res.append(res_best)\n",
    "        \n",
    "    results = { 'cvr': cvr, 'delay': delay, 'drift': drift, 'res': res }\n",
    "    \n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "digital-strip",
   "metadata": {},
   "source": [
    "## Do some simulations\n",
    "### Set up a paradigm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "brief-guatemala",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAW2ElEQVR4nO3df6xfd13H8eert91WYNCN3WGlLXdkA0PmaOddGaI4y6Zlkk1N1BlJpmIaCSH8iMLmEpKZaAANYowxWQRtAmIG8mNpHNAMRzJlnbesnZt0jOmAuUG7YUPKj8Lat3+cc28v9faec27P+Zzv59PXI7np/Z5+b+/7nNPv63zO53zO5ygiMDOz/KwauwAzM1sZB7iZWaYc4GZmmXKAm5llygFuZpap1Sl/2QUXXBAzMzMpf6WZWfb27t37VERMn7w8aYDPzMwwNzeX8leamWVP0leXWu4uFDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTLUOcElTku6XtKt+fb6k3ZIeqf88b7gyzczsZF1a4G8BvrTo9U3AXRFxCXBX/drMzBJpNQ5c0gbgl4A/Ad5eL74euKr+fidwN/DOfssrz73/9TT/9pWnxi4jG6tWiV+f3ciPr1s7dimn7QuPPs0XHvW+b2vVKvEbV2xk/fPy3/dDaXsjz/uBdwDnLlr2goh4EiAinpR04VI/KGkHsANg06ZNp1FqGd7z6QPc/7XDSGNXMvnmp6pfM7WKN/38xeMW04P3fPoA+77ufd9Gaft+KI0BLul1wMGI2Cvpqq6/ICJuA24DmJ2dPeOfHnH8eHDVS6f5+9/ZOnYpE+/48eDFf/TPHDtexn+b4xFs+4kL+eBvXzF2KRPvmWPHufiWOzleyL4fSpsW+KuA6yRdC5wDPFfSh4BvSlpft77XAweHLLQU/u/YXSkPjSplPVLyJlte40XMiLg5IjZExAxwA/C5iHg9cAdwY/22G4FPDVZlQSLAZ9DtzHc1RCEf4yC871tSvfN90Fve6YwDfzdwjaRHgGvq12Zmlkin2Qgj4m6q0SZExNPAa/ovqWxBLLQubHmltcIi8AXMluY3UylnX0PxnZhmZplygCfmPvDuSmmDVWcS3vttLFz/KGXnD8QBbmaWKQd4Yu4H7UaimGZY4H3f1sL1j5HrmHQOcDOzTDnAE6taFG6GtSXKaYVFeBx4Z4WcfQ3FAT4Cn0a35yGXZy7v+mYO8MTCLYrOStpkDqVuCtr1g3CAj8Cf4faqLpQyPsbVEFLv/bYKun49GAe4mVmmHOCJeRhhN1I5rbBqGoWxq8iHpGLOvobiADczy5QDPLFqSlE3w9oSKqYN5rOvbtwH3swBbmaWKQd4Ym6FdVRUH7hHoXQheRhhEwe4mVmmHOCJeUKjbsoaB+65hLsQKubsaygO8MSq+TD8KW5LBU2GUshqpKNyDt5DcYCbJeRDt/XJAZ6YH8rSTUnDCAlPztWFwKctDRzgZmaZcoCn5utYnVS30pfRDKuGEVpbHkbYzAFuZpYpB3hi1TBCt8PaKul26ghPZtVFNYywkJ0/EAe4mVmmHOCJ+bmI3VRTipbBfeDdlDSV8FAc4CPwaXR73lRnLu/7Zg7wxNyg6K6UVlh4HHhnhez6wTjAE/N0GB0VdDt1KeuRiuS5UJo4wM0S8sHb+uQAT6x6LqI/xm2VNYwQJ3gHJc1EOZTGAJd0jqT7JO2X9JCkW+vlmyXdK2mfpDlJW4cv18zM5q1u8Z6jwLaIOCJpDXCPpDuBPwZujYg7JV0LvBe4arhSy+A+8G5KOlup9n056zM4DyNs1BjgUd0KdaR+uab+ivrrufXy5wFPDFGgmZktrU0LHElTwF7gYuCvI2KPpLcCn5H051RdMT99ip/dAewA2LRpUy9F58z9oN2UNJkV+B6ALrypmrW6iBkRxyJiM7AB2CrpUuCNwNsiYiPwNuADp/jZ2yJiNiJmp6en+6rbzOyM12kUSkQcBu4GtgM3Ah+v/+qjgC9ituR+0PYKeqKap1HoqBoHXsreH0abUSjTktbV368FrgYOUPV5/1z9tm3AI0MVWRqfRrdX0kVM68a7vlmbPvD1wM66H3wVcHtE7JJ0GPhLSauB71P3c9vy3KLorpRNVk0lPHYVeSlk1w+mzSiUB4AtSyy/B/ipIYoqma9hdlPSzRylHIhSKekmrqH4TkyzhHz9w/rkAE+smpFu7CryUdKc0KWcSaRSzQXvbbYcB7hZQj54W58c4IkF4dPoTgp6Io/PvjpxH3gzB7iZWaYc4Im5FdZNWX3g4DFI7UkeRtjEAW5mlikHeGK+maObalOV0Q7z2VdXfqRaEwd4YtV/SH+K2yqpC6WUA1Eq1cHO22w5DnCzhHzotj45wJMLn0Z3oIJOo0tZj1Q8jLCZA9wsIR+8rU8O8MT8TMxuqqFkZTTDyliLdMq6/jEMB7hZQr4L1/rkAE/Mwwi7KakfNMLXP7oQnsyqiQPczCxTDvDEquciuhnWVjWlaBn8MI9u3AfezAE+Ap9GmzXzx6SZAzwxNyi6K6UVVt1K71jqopBdPxgHeGIeRthNUcMISzkSJSKVcxPXUBzgZmaZcoAnVg0lcxu8LVWPpS9CIauRVClnX0NxgJsl5GO39ckBnpjbE92ooGdilrMiaZR09jUUB7hZQr4HwPrkAE/NT2XppLqZo4xmmKdR6MbPxGzmADczy5QDPLHqdmo3w9oqqRu0mkbB2qoe5lHK3h+GA3wEPo1uz0Muz1ze9c0c4Im5RdFdKZvMfeDdFbLrB9MY4JLOkXSfpP2SHpJ066K/e7Okh+vl7x221DJ4RrpuyupCGbuCvJQ0F/xQVrd4z1FgW0QckbQGuEfSncBa4Hrgsog4KunCIQs1K4G7hKxPjQEe1Tn/kfrlmvorgDcC746Io/X7Dg5VZEnCwwi7KWoYYRnrkUpJc8EPpVUfuKQpSfuAg8DuiNgDvAT4WUl7JH1e0hWn+NkdkuYkzR06dKi/ys0y5GO39alVgEfEsYjYDGwAtkq6lKr1fh5wJfCHwO1a4vwwIm6LiNmImJ2enu6x9DwFnsyqC/eBn7mqPnBvtOV0GoUSEYeBu4HtwOPAx6NyH3AcuKD3Cs1K4mO39ajNKJRpSevq79cCVwMHgE8C2+rlLwHOAp4artQy+IEO3aig+6l9E1dH5ez6wbQZhbIe2Clpiirwb4+IXZLOAj4o6UHgB8CN4fMdM7Nk2oxCeQDYssTyHwCvH6KokgW4Cd5B1QdeSLvAI5A6EbgJ3sB3YqYWPo3uopqNcOwq+lHMgSiRahiht9lyHOBmCfnQbX1ygCdWDSMcu4p8VDPSjV1FP0pZj1R8K30zB7hZQj54W58c4Il5GGE31SjCMpphZaxFOiVd/xiKA9wsIV/Atj45wBPznNDdldIK820S3QiPQmniADdLyAdv65MDPLHquYj+FLdV0pSipaxHKu4Db+YAH4FbYe2VtqlKWx8blwM8MTcouiulFRa+ANJZIbt+MA7wxDyMsBt5QowzllTOTVxDcYCbJeSDt/XJAT4Gn0a3VsqFLA8h7K76lHi7LccBbpaQj93WJwd4QvOtMH+G26tu5sifG+DdlXL2NSQHuFlCvgfA+uQAT2i+NeHT6PaqVlj+zbD81yC9gh6HOhgHuFlCPnhbnxzgCc23Jnwa3V71TMz8lXAWkVr1MA9vt+U4wBNauIjp/G6vkJs5Thy8rS13oTRzgJsl5IO39ckBnpBbYd2V04UydgX58TMxmznAzRKSm+DWIwd4Qh5G2F05wwjzX4fkCpoLfigOcDOzTDnAE5pvhfk0ur1StlQBJxHJVX3g3nDLcYCbJeRjt/XJAZ6QGxPdeVL/M5cPds0c4GYJ+S5c65MDfARuWbRXjQPPvwnus4juPA68WWOASzpH0n2S9kt6SNKtJ/39H0gKSRcMV2YZFoYRuhXWWilzQp+4gD1yIRmRVMTBe0irW7znKLAtIo5IWgPcI+nOiLhX0kbgGuBrg1ZpVgjnt/WpsQUelSP1yzX11/xh8S+Ad1DG3c6Dcyusu2pGurGrOH0lrENq7kJp1qoPXNKUpH3AQWB3ROyRdB3wPxGxv+Fnd0iakzR36NChHko2y5cP3tanVgEeEcciYjOwAdgq6TLgFuBdLX72toiYjYjZ6enp06s2cyf6wK01FXIRc+wCMlTK9Y8hdRqFEhGHgbuB64GLgP2SHqMK9i9K+rG+CzQriS9gW5/ajEKZlrSu/n4tcDVwf0RcGBEzETEDPA5cHhHfGLTazC1MJ+vPcGul9IP6lvDuhEehNGkzCmU9sFPSFFXg3x4Ru4Yty6xMPnhbnxoDPCIeALY0vGemr4JKtvBINZ9GtybB8QIaYQWsQnruA2/kOzFH4FZYez7Ynbm855s5wBNyY2KFCthwbkmujDfb8hzgCflD3J0KGUY4vwqeC749lfJA1AE5wM0ScnxbnxzgKbkV1lkpN3MUcRaRmIcRNnOAmyXkY7f1yQGe0MJkViPXkZOqFZa/Es4iUivl7GtIDnCzhHzwtj45wBNamMzKn+LWqlZY/s2w/NcgvWoEki3HAW6WkC9gW58c4AktTGY1ahX5KaEVVsJZRGrVwzy83ZbjAE9oYS4Ut8Jakwp5Ik/9p3d9e+5CaeYANzPLlAM8IbfCuivlbmo/jWllSjj7GpID3MwsUw7whNwK606FPJIn8BjSrqQybuIakgPcLCHHt/XJAZ6QW2HdldIHXsZKpCUo4uxrSA5ws4R87LY+OcBTch94Z6WNA7f2PA68mQPcLCE/49P65ABPyOPAu6v6wPNvh5VwFpFaIQOQBuUAT+jEMEIneFulzAm9MBe8d31r1TDCAnb+gBzgZmaZcoAn5FbYShRyEdMXsDtzF0ozB7iZWaYc4Am5FdZdKUPJfAG7u1KufwzJAW6WkC9gW58c4Am5FdZd1Q+afzOshHVIz5NZNXGAm6Xkg7f1qDHAJZ0j6T5J+yU9JOnWevmfSTog6QFJn5C0bvhy87bwSDV/ilsr5WzFDfDuqj5wb7jltGmBHwW2RcTLgc3AdklXAruBSyPiMuDLwM3DlVmYQkIphdIOdmWtzbC8rZo1BnhUjtQv19RfERGfjYhn6uX3AhsGqrEYbkysjLeb2dJa9YFLmpK0DzgI7I6IPSe95XeBO0/xszskzUmaO3To0OlVWwi3LNqrhhHmn+ALQ0hL6RNKwMMIm7UK8Ig4FhGbqVrZWyVdOv93km4BngE+fIqfvS0iZiNidnp6uo+azcyMjqNQIuIwcDewHUDSjcDrgN8KX21o5FZYd6W0whamURi5jpwIT2bVpM0olOn5ESaS1gJXAwckbQfeCVwXEd8dtkwzMzvZ6hbvWQ/slDRFFfi3R8QuSV8BzgZ21y3KeyPi94crNX9uhXWnQm7mOHH2NW4dOSnl7GtIjQEeEQ8AW5ZYfvEgFZkVzAFuffKdmAm5FbYChdzMkf8apFfKRGZDcoCbJVTajUk2Lgd4Qp7MqrvqmZj5K+EsIjUhb7cGDvCEPBdKdyrkPNoH7xUoY9cPygFuZpYpB3hCboV1V04XytgV5EdQxs4fkAPczCxTDvCE3Arrrpw5oevrHz79ak0q4yauITnAzcwy5QBPyq2wrkrrA/eeb6+U56EOyQFulpCP3dYnB3hCboV1J6mIawcFrEJyhdwCMCgHuFlCvonL+uQAT8jjwLur+sDzb4eVcBaRWtUHPnYVk80BntCJLhQneGuFzAm9MBe8d31r1TDCAnb+gBzgZmaZcoAn5FZYd9WMdGNXcfp8Abs7d6E0c4CbmWXKAZ6QW2HdlXK24qcxrUAh1z+G5AA3M8uUAzwht8K6K+V26hOjKbzz2/JorWYOcDOzTDnAE3IrrLtSbqf22Vd35UwlPBwH+Aj8IW6vtNPostZmWN5WzRzgCbkxsTLebmcu7/rlOcBH4JZFe1UXSv4f4xNdKN77bcnDCBs5wM3MMuUAT8itsO5KaYUtTKMwch05EZ7MqokD3MwsUw7whNwKW4kynkzuYYTdlXL2NSQHuJlZphoDXNI5ku6TtF/SQ5JurZefL2m3pEfqP88bvty8uRXWXSmtMD+NqbtSbuIaUpsW+FFgW0S8HNgMbJd0JXATcFdEXALcVb82M7NEVje9Iap7WY/UL9fUXwFcD1xVL98J3A28s/cKgb+66xHu2P/EEP90Ut9/5hjgVlgXAr71naNc877Pj13KafneD+t97ysgHYhvfecH2e/7eX/6qz/JFTPn9/pvNgY4gKQpYC9wMfDXEbFH0gsi4kmAiHhS0oWn+NkdwA6ATZs2rajI6XPP5pIXPGdFPztpts48n80b3dvU1i9veSGHv/vDIoaTveKi57N547qxy8jGr2x5Id/+Xhn7HmDtmqne/011mSxG0jrgE8CbgXsiYt2iv/vfiFg2mWZnZ2Nubm6ltZqZnZEk7Y2I2ZOXdxqFEhGHqbpKtgPflLS+/sfXAwd7qNPMzFpqMwplum55I2ktcDVwALgDuLF+243Ap4Yq0szM/r82feDrgZ11P/gq4PaI2CXpC8Dtkt4AfA34tQHrNDOzk7QZhfIAsGWJ5U8DrxmiKDMza+Y7Mc3MMuUANzPLlAPczCxTDnAzs0x1upHntH+ZdAj46gp//ALgqR7LGUoudUI+tbrO/uVSq+usvCgipk9emDTAT4ekuaXuRJo0udQJ+dTqOvuXS62uc3nuQjEzy5QD3MwsUzkF+G1jF9BSLnVCPrW6zv7lUqvrXEY2feBmZvajcmqBm5nZIg5wM7NMZRHgkrZLeljSVySN+uxNSR+UdFDSg4uWnfIBz5Jurut+WNIvJqxzo6R/kfSl+mHUb5nEWlfy0Oyxtmn9u6ck3S9p14TX+Zik/5C0T9LcpNYqaZ2kj0k6UP9ffeWE1vnSelvOf31b0ltHrzUiJvoLmAIeBV4MnAXsB142Yj2vBi4HHly07L3ATfX3NwHvqb9/WV3v2cBF9XpMJapzPXB5/f25wJfreiaqVqrHXj6n/n4NsAe4ctLqXFTv24F/AHZN6r6vf/9jwAUnLZu4Wqmep/t79fdnAesmsc6Tap4CvgG8aOxak674CjfWK4HPLHp9M3DzyDXN8KMB/jCwvv5+PfDwUrUCnwFeOVLNnwKumeRagWcBXwReMYl1AhuAu4BtiwJ84uqsf99SAT5RtQLPBf6bejDFpNa5RN2/APzrJNSaQxfKC4GvL3r9eL1skvzIA56B+Qc8T0Ttkmao5nTfwwTWWndL7KN6LN/uiJjIOoH3A+8Aji9aNol1AgTwWUl76weLw+TV+mLgEPB3dbfU30p69gTWebIbgI/U349aaw4BriWW5TL2cfTaJT0H+CfgrRHx7eXeusSyJLVGxLGI2EzVwt0q6dJl3j5KnZJeBxyMiL1tf2SJZSn3/asi4nLgtcCbJL16mfeOVetqqu7Iv4mILcB3qLohTmXsbYqks4DrgI82vXWJZb3XmkOAPw5sXPR6A/DESLWcyqke8Dxq7ZLWUIX3hyPi45NcK7R+aPZYdb4KuE7SY8A/AtskfWgC6wQgIp6o/zwIfALYOoG1Pg48Xp9xAXyMKtAnrc7FXgt8MSK+Wb8etdYcAvzfgUskXVQf/W6geqDyJDnVA57vAG6QdLaki4BLgPtSFCRJwAeAL0XE+ya1VnV/aPYodUbEzRGxISJmqP4Pfi4iXj9pdQJIerakc+e/p+qzfXDSao2IbwBfl/TSetFrgP+ctDpP8puc6D6Zr2m8WlNfAFjhRYNrqUZRPArcMnItHwGeBH5IdZR9A/B8qotbj9R/nr/o/bfUdT8MvDZhnT9Ddcr2ALCv/rp20moFLgPur+t8EHhXvXyi6jyp5qs4cRFz4uqk6lveX389NP+ZmdBaNwNz9f7/JHDeJNZZ/+5nAU8Dz1u0bNRafSu9mVmmcuhCMTOzJTjAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8vU/wG2Bfdg/KZHFAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "etco2_mss3 = np.concatenate((30.*np.ones(120), 40.*np.ones(180), 30.*np.ones(120), 40.*np.ones(180), 30.*np.ones(120)))\n",
    "p = plt.plot(etco2_mss3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sharing-probability",
   "metadata": {},
   "source": [
    "### Generate and fit MRI signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "boring-switch",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean CVR: 0.04001365530413544, std CVR: 0.001675312412216947\n",
      "mean delay: 19.95, std delay: 0.2958039891549808\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydZ5gUVdaA31MdZoAhg2RBzAqCLGDCnNOa17SucVlcP+Oua845rau76pox57gGxBVZQEWC5CQ5Z5hhYJjpUPf7UV3dVdVV3dUTCE29zzPPdFfdunW7wrnnnnvuOaKUIiAgICCgeNG2dgMCAgICAhqWQNAHBAQEFDmBoA8ICAgocgJBHxAQEFDkBII+ICAgoMgJBH1AQEBAkRMI+oB6QUQ2ikj3rd2OfIjIESKyZGu3o74QkVtF5KUtcJ5LRGRUQ58noGEIb+0GBGxfiMgCoB2QtGzeQylVVsv6jgDeVEp1rnvrdjyUUg9u7TYEbPsEGn1AbThVKVVm+VuWq7CIhLZUwwICArIJBH1AvSAiSkR2S30eLCLPichXIrIJOFJEThKR6SJSKSJLReSvItIE+BromDL9bBSRji51nywiE0Rkg4gsFpG7Lfu6pc59sYgsEpE1InKbZX+jVHvWi8h0oJ+P3/FnEZmdaut9IrKriPyUOv/7IhK1lD9FRCaKSLmI/Cgi+1n23Swic1P1TBeRMyz7LhGRUSLyeKpt80XkxBztuil13SpFZJaIHJ3afreIvGkp9wcRWSgia0XkDhFZICLHWMq+LyKvp+qZJiJ9/bQ3YDtHKRX8BX++/4AFwDEu2xWwW+rzYKACOARDmSgFlgOHpva3BPqkPh8BLMlzziOAnqm69gNWAqen9nVLnftFoBHQC6gB9k7tfxgYCbQCugBTc50vVdfnQDNg31Rd3wHdgebAdODiVNk+wCrgACAEXJy6PiWp/ecAHVPtPhfYBHRI7bsEiAN/TB17JbAMEJc27QksBjpafvOuqc93Y5i+APYBNgIDgCjweOocx1jKVgMnpc75EDDacp587R21tZ+/4K92f4FGH1AbPk1psOUi8qlHmc+UUj8opXSlVDWGwNlHRJoppdYrpX7xezKl1HCl1JRUXZOBd4DDHcXuUUptVkpNAiZhCHyA3wEPKKXWKaUWA0/7OOUjSqkNSqlpGB3DUKXUPKVUBcYIZP9UuT8CzyulflZKJZVSr2F0DAem2v2BUmpZqt3vAbOB/pbzLFRKvaiUSgKvAR0w5j+cJIESjOsXUUotUErNdSl3NvAfpdQopVQMuBOj47IySin1Veqcb1iuk5/2BmynBII+oDacrpRqkfo73aPMYsf3szA0yYUi8j8ROcjvyUTkABH5XkRWi0gFMAho4yi2wvK5CjAnhzs62rLQxylXWj5vdvlu1t0V+Iul0yvHGDV0TLX7DxazTjnQw9HudJuVUlWpj1mT2kqpOcB1GBr5KhF5183EheO3pupc6yjjvE6lIhL22d6A7ZRA0Ac0FDZNUik1Vil1GrAT8Cnwvls5D97GMKd0UUo1B/4NiM92LMcQviY7+zzOD4sxRgstLH+NlVLviEhXDHPS/wGtlVItMEYHftttQyn1tlJqAEbnooBHXIotB9LeSyLSCGjtp/76bm/AtkUg6AMaHBGJisiFItJcKRUHNpBxz1wJtBaR5jmqaAqsU0pVi0h/4IICTv8+cIuItBSRzsDVtfkNHrwIDEqNOEREmqQmjpsCTTAE8moAEbkUQ0MuGBHZU0SOEpESDBv7ZuzurSYfAqeKyMGpCeN78C+o6629AdsegaAP2FJcBCwQkQ0YppffAyilZmLY3OelTAZuJok/A/eKSCWG3fl9lzJe3INhrpkPDMWwS9cLSqlxGHb6fwHrgTkYk5YopaYDTwA/YXRmPYEfanmqEoxJ5TUYppedgFtd2jMNoyN7F0O7r8SYLK7x8Vvqs70B2xiiVJB4JCCgGBGRMqAc2F0pNX9rtydg6xFo9AEBRYSInCoijcVYo/A4MAXD5TNgByYQ9AEBxcVpGP74y4DdgfNUMGzf4QlMNwEBAQFFTqDRBwQEBBQ522T0yjZt2qhu3bpt7WYEBAQEbDeMHz9+jVKqrdu+bVLQd+vWjXHjxm3tZgQEBARsN4iI56rvwHQTEBAQUOQEgj4gICCgyAkEfUBAQECREwj6gICAgCInEPQBAQEBRU4g6AMCAgKKnEDQBwQEBBQ5gaAP2GoMnbaCVZXVW7sZAQFFTyDoA7YK8aTOwDfGc/4Lo7d2UwICip5A0AdsFfRUML3F6zZv5ZYEBBQ/gaAP2KooXyljAwIC6oIvQS8iC0RkSipD/LjUtlYi8q2IzE79b+lxbAsR+VBEZorIDBE5qD5/QEBAQEBAbgrR6I9USvVWSvVNfb8Z+E4ptTvwXeq7G08BQ5RSewG9gBm1bm1A0RCkQQgI2HLUxXRzGvBa6vNrwOnOAiLSDDgMeBlAKRVTSpXX4ZwBAQEBAQXiV9ArYKiIjBeRgalt7ZRSywFS/3dyOa47sBp4VUQmiMhLqVyWATs4gUYfELDl8CvoD1FK9QFOBK4SkcN8HhcG+gDPKaX2BzbhYeIRkYEiMk5Exq1evdpn9QHbK8EkbEDAlsOXoFdKLUv9XwV8AvQHVopIB4DU/1Uuhy4Bliilfk59/xBD8Lud4wWlVF+lVN+2bV2TpAQUEXpKzguydRsSELADkFfQi0gTEWlqfgaOA6YCnwMXp4pdDHzmPFYptQJYLCJ7pjYdDUyvh3YHbOcESel3PCYtLqe8KsZLI+fR7eYvqYoltnaTdhj8pBJsB3wiImb5t5VSQ0RkLPC+iFwOLALOARCRjsBLSqmTUsdfDbwlIlFgHnBpPf+GgO0Qlf4fCPwdhdOe+YG92jelvCoOQHlVnMbRumczfWTITJ4bPpcFD59c57qKlbxXWSk1D8Mt0rl9LYaG7ty+DDjJ8n0i0NdZLmDHJlDod0xmrqikXbMSAKSerHbPDZ9bPxUVMcHK2ICtQyDodyh0PXPDVQPNzyT14KHyIhD0AVuFwGSzY6FbhnDmp/rS6E1iCb1+KywiAkEfsFUITDc7FknlptHXjkteHcM170zI2l6TSNayxuInEPQBW4VAzu9Y6C7K9osj59WqruGzVvP5pGXp75GQ0WXUBBq9J4GgD9gqBO6VOxa67X4bn18cOb9e6o6EDDFWEw8EvReBoA/YKgRifsfCzXRTX5iCPpYMTDdeBII+YKsQKPQ7Fjavm3quOxo2xFh1oNF7Egj6gK2C6XUTCPwdg6TNvbJ+b3rUNN0ENnpPAkEfsHUIBPwOhdXFvb5vfWYyNjDdeBEI+oCtQiDndyz0LWGjDzR6TwJBvwMxZOpybvl48tZuBhCYbHY0GtJ0EwlMN3kJBP0OxKA3f+GdMYu3djOAYGXsjobegD17JBwI+nwEgj5gqxBo9DsW1gVTyra97g9C1LTRxwMbvReBoA/YKgRyfsfC6kdvvfnJeujxUyHUG3TUsL0TCPqArYKXnVbXFdWBZtagKKUYMnU58eSWM3W4BTWDwiNOJlzabMbM2YI/Z7sjEPQBWwXzvXe+5rd8PIW97hiyxduzIzFs5ioGvfkL/xw2Z4udU/eYjPWrhV/zzgQufXUMu932ddY+MwpmoNF7U/f0LgEB9ch747aNyeJiZu2mGABL12/eYuf0MtEkfGr01iBmTsy49oGg9ybQ6AO2CvneySDoWcNRSHjg0fPWUrE5XudzJj1CINTHZKyp0QeJR7wJBH3AViGfe2Ug5xuefPdgcyzJeS+M5o+vjXPd/+PcNcxYvsHfuZT757oK5/u/mM6Pc9cC9tW3AXYC003AViGfINeVQqvnVHMBBuIztVMi5RM53UOYX/DizwC2pNxDpi4HhBN6tLeVtWv0mc918bpZtynGS6MyoY7rY3RQrAQa/Q7ItmAWcbagOp7klo+neO7PR3lVbIdZAq+UYurSinqoKPdus0Mo5HkZ9OYvDHpzfNZ2rzDFbglJ/NLnvm89zxFgJxD0OyDbwvvgFB4fjl/CO2MWpb/7mVjTdcUdn05lzqqN9L73W1cBU4x8MG4Jp/xzFP+dvrJWx5v6fL4rbN6j+nhclJd7ZT0+jMFkrDeBoN8B2RZeCGcLnG3y08Q5qzfyxuiFXJkS8MNmrqqn1m3bzFpZCcCCtZs8y1TFElnbhs1cyTfTVvhOym1aQqpiSW78YFLB7bRi83G3afTGlx/nrGH3276ivCpW63MEphtvAkG/A7ItvA9pP3pTa1Tu+3Nh2n01v5KrSKmsjttC9I5dsI597vyG//262lbussHj+NMbmVFPPpOMdf8H45fUqY2eNvrU9qe+m008qTznA/ydo/btK3YCQb8Dsm0EFMsjZHy0MS3otR1b0Pe8eyjn/Pun9PexC9YBhleMG2//vMh1u5P6VAi8OhXTdFOdml8pjYRqfY5tYaS6rRII+h2QbeF9MIWI14SfHyHjxzWvYnOc3/37Jxavqyq4jdsqbvdv8pIKfp63lunLNvBMnhWv4xau99y3oTrOqNlGB1GfgtN7Mtb4YgYkKw0Hgr4hCAT9DoiuFM8On8O81Ru3WhuyTDVZ+31o9D7KfDN1BWMWrOPp72YX0Lrtk3NfGM2p/xrFppi/WEFuV+/wR7/n9y//THlVrE6C09kJey2YMu/h5pSg1+ogkYIFU94Egn4HpGJznEeHzOK8F0ZvtTY4TTNOmeLnnTW1QTfDzYbqOBVV8bSw2lHs+HUVduurjFWwNQndt7uNUoql5fZwCs7AdPYMU9k2erN8XdwtAznvTSDod0DMl6uyOtszY0tRLxq97u3+t9/dQ+l179D0y18XTXFbo776rFyXOKkr34LzuxmrOOThYXw3I+PuWRVLcuazP6QzmnnFo88IeqNAXUYRgenGmyJ6/AP8Yi4s2poLTLK9bAp3r/TTfrOM39Wg2wN+b5vkWVmcqxpD0Ps70fqUS6Q1e1nF5hi/LCpPb/Oy0Sd0xUNfz0jH06lPc1FAhkDQFxmTl5TnXTVpplzbmi9GPq8afwumfJwnbbrx1ayAFPkEvdVnvWlpBIAl6zMT3pMWV3iWt7Kioprn/zcvU85SbPisVVTHk77zEwQavTdBrJsi47f/+gGwxx9xEtsWBH2eU/tpmR+N3hQwoSLS6OuLXOaxpFKElPc1s4YXdhOwzgVbXo+aM/nJvf+ZxpPn9qayOsElr47lggN29u8OGmj0nvjS6EVkgYhMEZGJIjIuta2ViHwrIrNT/1vmOD4kIhNE5Iv6anhA7YltQytLvBZM+dHOkimVPpcId7pxFgPW0VAu4VboT7ZqzrqLRj/w9XGsqKgG7EqC7nIPY0mH143H/dzs8BD6ZVE5f3pjfDpm/txV/j3DzHOMnreWbjd/WVQutXWlENPNkUqp3kqpvqnvNwPfKaV2B75LfffiWmBGLdsYUM8UGvyrvCrGwhzL7WtD9mRsntlZF8z+KlfRYva6EZGcPvH5sF63VRuqbZm9Ei6TsUOnr+SPr49L7c88Q5lJ8cwBTk3dq0NycwWtrE5kPKoKuG3mKT5MreL9ad5a/wcXOXWx0Z8GvJb6/BpwulshEekMnAy8VIdzBdQjhWr0x/9jBIc/Nrxe21Af7pVpAeM42OazbXrdFJ+cZ+Licn73/E/5C/pgkUP79bLRr0tp2m4avZW4RZmYtqzCc4RWVZPt+bW0fDOTlpQD+SeUrZidg3mvA1NOBr+CXgFDRWS8iAxMbWunlFoOkPq/k8ex/wD+BuSULiIyUETGici41atX5yoaUEcKjXq4ckNNvZ176tIK3vhpQX73ygJCIFjNMp9OWGqzD5vD+WIMk7CsvP5SATqvtq6Uqw3fvOY2G705srKZbjKv+8lPj/KcD/Ja3DUytTq3EI3ePEdIM1ML+j+22PEr6A9RSvUBTgSuEpHD/BwkIqcAq5RSeePHKqVeUEr1VUr1bdu2rc9mBdSGt3xObjUEp/xzFHd8Ns3Fb97+feaKSq56+5csE4CVhIvbzXXvTbTZfU1NspgsN+a1qrOXieVwZ1VJXTF8VrbCZZ7TKrjNztRahXPUOGN5pWsT3KJsAjQtLdxPxDkfE8Snz+BL0CullqX+rwI+AfoDK0WkA0Dqv1uM2EOA34rIAuBd4CgRebMe2h3QgJz57A/c9OHkBj1Hlt+8Q/Rf8/YEvpy8nF9XVvLG6IV87whBPOLX1Wk3USdWLfHRIbOA4vS6qatpQqGoqIpz2eCxWaODpK64/8vsaTVT0Ns1ejfTjX3bKz/Md23Dphp3jb5x1Ih5U5iN3u5htXhd1Q6TjCYfeQW9iDQRkabmZ+A4YCrwOXBxqtjFwGfOY5VStyilOiulugHnAcOUUr+vp7YHNBC/LCrnvXGL8xd0obI6zkUv/5wVItdJPhFVmbLdlpWEuePTqVw6eGx637CZK/nDK2N49nv34F3WkL0mxTgZW1eNVSkYNmslw2au4oGv7ELdy9SiKxg6bQWXW+6HWdTaeecaiVnZ5GKjh4w5qBAbfdJho39hxDxu/qhhFZbtBT8afTtglIhMAsYAXyqlhgAPA8eKyGzg2NR3RKSjiHzVUA0O2LZZtK6KkbPXcPErY9LbKjbH6Xbzl7w3NmMy8rsS1u1FX1FhzBksWOvuPuemxRWhib5OcWFM2pSVALC60j4P49WJxBM6A98Yz8wVlVllbaYbn5p0ZY2xIvbwPezm2qqUq2dtNHrrfMyQaSv8V1DE5DWEKaXmAb1ctq8FjnbZvgw4yWX7cGB4bRoZsG2g6yrvpGYimS0gzBWTr/6wIL0tnrQLh3zxyq3km6h10ybz+dFvrEkwa8UGftO1Vc5y2xL5bPR+ZGTCQ3P30ugrXTRwN9ONX8+u8lQQtcsH7GIbBbp54+TDzZV2s89VtcVOEAIhwDem0B05ezXjF65zLeMUshe/Mobr35sI2IWtW4fgxpeTl2XVnc807Wa7f+q72dz3xXTPY65++xfOeu4nKlKCZ3ugriublfK+D4XUne5wPLxucmEKetMmb1LlM9Ryd1lGBwx/eV0HZnxB84ThsRMhwd9C78Dmcl91FTOBoC9Syqti/OGVMayqrK63Ok2hcNHLYzjrObv/9uZYkn8Nm52OQgiwZmMN//t1Nb+uNFY3WjXMuMPu4KWcPj701/TndLTNLLOP/XvcQ3i9PMp9QhBgylIjhV1NcvvRAOvDqyThIZAL8ehxiyIa92m6MXPENsoS9Alb3e4ohpX8lW9LbjS+6jF470IumjGIUmqYXfoHrgz/h43f3MvGWowQiolA0BcJFVVx29D3vbGLGfHral4cMS/HUYXh5s5o8tR3s3l86K+8a7HD973/v951WYSxUoonvv3Vs6xJrrDEVord08IrbISTfPbt9VUx4h6C1O+IC6whEDLHDLWs1ThQm86J2s9Zx7VlPZtiSaJhjUhIS20rp7OsSptcvDyrAHYVY7RXJoYyE00Yq7db1iylm2TO/834WfzuuR99/55iJAhqViRcOngMvyzKDFHNl9yvYvbh+CU0joY4qWcHzzK5Xn5TA/Mb495q4nEmrfAisxLWf91+SQvF7cj1Op/WPWFRObd9MoU/Htqdrq0bZ+3/ef46+nR1D1FViEZv9hVeyve70fsB6Fb9dnpbF1nJyJLr+UfiTEpCGiUVHSmlhrGlf2ajKuW2+LVMkp2oiTclSpzfaL8yRt+LJBnNv6vYF/61jWcSmP8r8nT681mhUSxc3R7wtfynKAkEfZEwYbHdDunXLS2pK0Ka8NcPJgH2qJdKKbtdXVeeQ32zVC5t2qphWoWxX3vsJa+O4YWL+uYMeZCvDZ5tS/3fjuR8Xjv6j3PX8uPctbz18yKO2buda5kZyze4bveapM3Vjvx2fQUITdhMZzHs6NeFPzZ2vf0hJ2h/BgwN/Sn1CJTAnPJdmBNpwwmhsfxf7Gp+1Pflu7ZPcP/ao9lLy4weJ5QMpOWKTAC03bTM3A7AiaHsEcWORCDoiwQvBUxhxKg37eRO4kmdkJbRktZszLjZ1SR0SiOZfQld99TYzQ4hlzZtSzhhGR34FfQzV1Ty4FczOKC73TMmS9DXITrn9rSYspAJ0//OyA57ESVOk+oVgN0+foQ2EW1Tx6zygs5p2o98oR9IIiU6SoixaOE8W3v6ykzKpJoEISbr3dPHt6SSw7TJPBV9lpcSJ2bV/4/os1nbdtPns1vImFv5V/SfxsZKeCJqzwHcUuzP9zpVRivLtmlqF/bGUF5eHjWfM/t0plWTaNb5ipVA0G/nvDNmEd+6xK4xtecVG6rTMerdcGpuVkFeFUvaBX1SURHP7ZWSS9BbzQFWe7/Xohk3NC1bGDsndmul0W9Hfvbmz/fSutuxjg+i93BR/BYWqva2fYLOo+EXeCd5FBeGv+OslSP5L4OpwRB6YRIMjj5KxQ8fs5dcwTrVlFUY5p3jtHH8I/os/RMzeT15HM3YxPXhjzhowXTe423QE1wf/oBrw5+kz/dBImMumVA6KP357NAIAG6IDeLv0X/n/L0/JffhoJC3x5Qbt8Uv57noU+nvMWU8x9OWbeD+L2cwcvYaXrusf0F1bs8Egn4755aPp+Tcn8/7wWmKsX53Cu2Erqisdhf0pqB0xiG3Yu1E0n70Cs863c8jWeYV59xBbWz0Jn6CqW0reJnRTgmNZmdtNZeEvuGexMVo6OgIIDRnE+eER3BiaAxhjJHUIdpUemtz+EfibFphLIRqvnEOQ0qMyOM/63txbuxOmomxHuKC8DAuCA+znfO68Idcl/g4S6KcEx7h2sYWYkycfq4fzKSaXRna6glCm9wXN50fv53xbR6l9fqJrvvn6B1tppozau5hgtrdVkbDHvNo5Yb680bbHgi8booU05SSb3TvNHNc917mZXJqjEldz6st5xKy1klXa7kNBSQp10TyulPWzkbfsBEPlVL1nggjnlQcIDNYUHoBnWU1zdlIGVVspBEAZRjXe1bJxbwaeZTOsopXo48Z+6QalfrNr0Qf55rwp9wafoveWnZYiQO0mYRJ0B73tRNgsbUXQLlqQoIwc1UnKg/6W3r7GtUMgKQSXk6ZeCb1f5xvkn3TZX7VOzEwdj3XxP6P82J3cGbN3el9K5UxAnkjcUx6WygVPLckbGj2ubx5ipFA0Bc5+bwnnNrwtGWZyblklqasPM0FpqD0q03f85/MUNyvp45xHiMOuxWn26fVl993vQ0cw/yN0Qs59NHvmbyk7ot3yuJreCbyD0oSlfw+/C0AJ2ujmVQ6kDejD5JImSmaSDXnh74jIkmODE3i4fCL7G8R5I0kZqv3ivDXvBB90vWcgyOP8JfIh2xUpeg5UgxauTv+h/TnE2oe5i+xQbb9E/Td0p+re1zARbGbeTB+Pv9KGKkt+tU8x32JiwCINe3Cn+I3pMsfF3uMoXo/PtcPZg3N+UXtwXytK3EVYh1NAbgjcRndqt9msd4WTdJrsI3z7WArZgNBX6S8OXohkH/CLpfLpFOAJnVlq++hr2YwY/kGhs9alTHd1EJTqqiK5S+UYmNNgi8mL7e3y/EbarPsPe1100Aa/Zj5hjY8f03dM3WdtOQpTg6N4RhGc2poNAC9tLkA9NbmcV/kVaNcaAwPRV5OHzcgNK3W5zSPfSt5NIuUkXpik9Ysvd8qtE1G6PulP89X7flIP4y1qml622vJ49OftZAwUt+PF5KnMjh5PHtWD2YdmfrNyBsTLZO7Tv7U+An2r3meakps23WEsJieQca2QgT9qNlrOOKx77frziEQ9EWKKVDyafTOiUwrzk4intRtGv3zI+Zx4lMjueTVsSxIna829vHVG/0LereOxPkb6vJC1jnG+xagbfUCAI7WfklvOymUCSLXWOovUYyVj5MDeCxxLjPVzgBs1gzf/BvjAzkjdk+63DuJIwFYqtpQfvILzG12QHqyd0DNU/Sufp5Da57kR+mTPsYeXVTS5Z37z43dydJBv9LXxf+/KhliI9nrBZJoaUFvKi+FjPru+c80FqytYqFHED2T6niSR4fMLPj5GzJ1RYN3IoGgL3LyafRHP/E/z1CuTjPN37/9lWEubnoA5ZuNCdXa2D6tLp21wTkqcSac9kNmTqNhBL1Za30kKY9ppQAcF8rO5/NF8kDf9dwYH5i17dH47xit723bdkXsL3Srfpsb4n8mQZhhem+m6N0Y3OpaVqvmfJ3sjzXAxe2Jy+hd/Tw1RIntdRrv753xftlMKeU0ZbGy+/XnyxdgZo2qIYqUNmPwZf05qHtrWxmvZ08hhCUVJylVxBz1JZI6Q6et8AyqB/49sl4cMY9nh8/l9Z8W+DsAmLKkgkFvjufOz6b6PqY2BIK+yPHja/3uWPfY804BOnL2Gl77aaFrWVOTL1SjFylM0LtlJHJ2SLWNpQ8NqNGnqq0PL86aUBPX7X+N/4k5KuP//kNyXwBuj1/KJlVCv+pnqVKGWePB+Pl8mhyQVcf7ySP5Kml3O1yi2maVOTX2IL9E+tCv5rm0Fn16zb3cGB9IkhDlKTu5JkI4JaQbRUKc1LM9D53Z06hI4Ppj9uDVS/rlzRdg3a2JUFYSZr/OzW1lvNx0k2iEUsc74wM98/1cBr4xnmEz3fIm2cnnkVWdyoPgFWvJDbPDmbe67ia9XATulUWOl+BqWhrOOwmaK7aNE9OkUshDblKIoM/X5m6ynE2qEatpUVAb0pOx277lJmsiYYrejZ7aAr5N/oZdmtRwXdLwgLkkfhPN4ptYS3PeTB4LQCKl232ePJi4y+tfThPKVZltW+O2XWFV9oVxatAT1W5MTNpt9ZpIWltvFA3x7IW/Yd2mmOEWrODaYww3yA15XGytHYGWUk9LwnY91WvhnY6Gpox9ScczvXCdIWDNpOdu+F1lbt6WQgZtkVBhTgy1JRD0RY6XRt8kml/QF7Ly0nTTrE3o3I0FeN0429xNltOCTcxTHZhQMpCQKKpVhL1qXgPg9cv68wdLEpR8NJzpxu7HXRe0lKvgffELeTl5Es3YRG9tLhWUUV3WHiqMcnHCrMWu9W6mhGZspoYIAIfWPEkLNtFdlnF2aAQJwozUe7JaNWNQ7HoWqva0aVEGZOd89TPxrgmEUpLZ/Olul8Cv6de7g3wAACAASURBVMao0/h81m868/Qw9yxjJtGQho6GrpuC3r7f9LIK+chMk+/RSJvnChi3mcHcGtrdMzDdFDleCnaTkpD7DgsJXeW0XVqpcZnccmpcbiiVecgHaFNYUHoBrXCPvwLZmtfwkr/wacmd7CWLCKUm3EolzlGhCXx1VT8O696cv4TfpznuISBMMu6VeZtcK9LaXj0Yb0TpjNb35uXkyYCwgTJG6L3o160lN520d85jL4rdwr8Tp7I+ZVpZrNoxRXXnM30AF8VvBWA9zehX82/Gqz1ZQ3PbvMKZfTqlPx+6exvXczx5biZPkYgQcjwGae1cXLZ5/WaXsl1bN+G/N+QOVPbQmT3REVo1NnTaLE+y1H3JJej9BgisjUYfTmn0dQnb4YdA0Bc5XoK6SUn+wVxSV75NGW5Dz4jzDXehHetoETfso1enls5bg1WBsSy/eyokrflC7MR67gy/ni7TRew5al+JPMY+X54Jr57I1eFPuTX8Nn7YHrxuBN3Vl/3qo3anLM99/VV14eHE+RQyW2CVgafsl4lu+pfj9qRj81KX8lbtO6PRpzVel8dCy/OoWOu0av/hPAd2a9OEkmgkvWDK2ZG7ZaXyw+ZYMstTxnzXCqnJfNwa2nQTCPrtiKSueG/sIs+l717HuNEkml/QJ3Tl2yfd7UEdoE2mGZsoIcaQ6E0crE2lI2sYoGXCNvxc+n/8WHI1J2hj6JSKaGguVz9Im8bfwu/yQPgVhpX8lbZkFhvdFnmLy8JD0t+fcIuXsmIKLB0HwO7akuz9FjIrYxvIdFMLbc8LUYqky6sbDgmawF3xi3kxkZXNs/bns7Q5GsqMBEOa0Kllo6zydkGf0ehzCcJ8gta639pRRPKMGkOaoBDMwBnOyVjTdONH0FsnY/e5awj9HvivY3+qfQXcY/P9jCcaVsEIbPTbEe+OXcRtn0ylYnOcgYft6usYT0HvQ6NPJHXW+fRxt9prG0dDdI/P5t/qfl4PH8vLyRPZS1vM8dpYToiMpZ2Us0v1myiLsPp39B/pz80xJsgGRx6hRDI2+V7aXEbo+7GPLOS0UGGJJIywuPld6LaHyVghiSLb9BYJaYiIbSFSbbjqyF155vu56e9WIWhOHpq4adRWM0hIkyyN3k2o5hO01sGhrT15bOshEXRCnpOxybSN3rsO03Rl7SOMGE32+aJMR+Zf0pt1BqabgDRmfs31LnlNnS+giZeG6rTRi8BvtR84VJvM8Oj1dJGVJHTFmk12j5gwCS4PfUkzh807lkyyqyxFUuGzTkwt4AmTpCnGQpNe2lzaiaGVt6aSMO6TsM9En+bK0Oc2IQ9wW/hN7g4P5tOSO12Ps3Jz/I9w8yLY7zw2qEbsJOXpGOi52B5MN5rS3TV6TQo2QbjhNLlZaww79oVdnjur7BUh7dqYa1STby5UHOYgr7Zm1auBEg2zm8majE016tEhs/LOR+V7NPQCR23Lyjen36+GzooWCPoiIerxwHtp9I0dppuwBk9Hn+GN6MN001ZySWgoSV3RaeifuDs8OC2Ue8sc7oi8xVvRBwmR5DRtFM3YyFFqLN+V3MiD4ZfZQxZzmGYswrogPIwnIoZZpbeWSWvYW5vDuJIrPX/PTZF3bd//kzyQXbSVXBD+3rZ9lcp2o3wyfhbv6UdAaXM483kuid0EwB7i7V+fCYHQwF43ecrFkzr3fTGdtTlcTgUd3eXVjYS0vALTD1mCXuwaeq6yYNfyNRFCqTJOjde2FtaH6cZ8xq0as1tHY8Uw3VjdKx2mm9TXeWs2sdbDxdI8Q74cvZmOLHebut38JY8MmcnBDw/j0lfHAoGgD3BBgEVrq/j3/zLD65KwhpD9sFif6zZlmRggZRaNvilVvBB6zHbc5eGvabJ2Ku2WfMMl4aHclZr4NBM89NQW8E70fp6KPsuzkac4LWTEvD8//D2fyF/pqS1I17Wni338pegT6VC1fhit72P7/vvYLXSrfpsBNU9llX0qeSZWMWIuIjJzjOai4aJXGv/zaXtfT13By6Pm89DXMz3LiHIX9IaNvu6SPuwQ5tavzo7EzVslFLJq35LlOpn2YimgTSERHjijB80bRYha7PL5NPqQCLpNo3dGZM1897pyZnvT9nTHsOD7Wav4ed5aXvlhflY9m2oSXPHa2HTkVnN+7bnhc211NLTpJrDRb6cc9phdsx2pX8SoSA9bhL9GVFOShJHRa1lHU55tejPfbDQWw1wQ+4DrSp5hjurEKL0HR2oTss5x+PCz05/PCI3i0cR5tLBk7emvzQJqHyzrncSRHKhNZxctE1bh3vhF3Bl5I6usGXp2qWrNgJqn0vb9GBHOi93OUdoEBoa/TJW2v7IbKKNCNaazwzPHiqmF1WYdQH1irilwClswAtW1KYvSS5nx5e2ENQ29HvxDncLbbjYRDuremjNSbpZuJsOwzec9870uV1YEzunbhXP6drFtz2+6ERBBS12XrMnYgnLjGmWdtnlTK7e21WTI1BX8d8YqmpbO4slze1Ox2X/uhfokEPTbEbnMCk2o5vjQOIgbdvT/lVxPJ1nLS/Hz6KKtpgureX79QP4eOptm7bqyy6S/g0AvmUcvi0nFizKpZkrpFenvc/UO7Kot573EEZwW+oFS8f8AJyXCMdUPs0C147rwx1yrfcyzid/yROIckoT4KHkok0qNOCz9qp8hSoKdNcMFc2Syp20SFwxtf7S+D2P1PdMdkfNSxYgQwduDyHw3GzrWTT7jjRniwWlaA7j9UyMeyujm3jb6hGNzz07NmbK0oqC2Zgl662eBdwZm4umYZppoWEubH6zHi4ghbCF9EaIhja6tG3PDsXv4bpPXSMV6LrfV3uZkrCjTvdJbo8+HWXZDHmFtbanpt292dm7za1uCQNBvx7RjHeeGhvNB8vD0tl4yhxgROslaAM5PfAbAdL0r+2gLuSHyITnyR7gyqcnB9Npk93I5OvZ46pPwWOJcHoq8xMHaVJrkiJz4h9hNHKxNo8Vhg5g/zFgUtUEZ7nkT9V1JpjxJKihjUOw6yiljdSqN3VK9DbfGL+fT5CGe9X+r9/Xcl0RLryh1xWNRzOyVlezUrJTmjSLex7pQURVnVWU1U5ZWcMP7kzhsDyNejCmvVlVW06w0YkvVCJll/JGw0O3mL3n4zJ6c139ne1NzmW6SdoF40UFdad0kyuWvGW6mVo+aRpGQq/usU6i6LVZKnzMlwDo2L2VBKrqj0xPHqdFrmvC/G4/MOm8u8vnZn9qrI/+blR2vJqQJSjJmzV9X2lf4Wjt2662/9NUxdGjRiAfP6JkpmxL0+bRyzdL5mHGYzLmE8lRIbpEtm584EPTbMJtjScIhSQ9PrUPoa0IfG0Ib0v8BPnN4pDRJZRk6O3YXZ+8u3Lv40uzzqGhWEoqrYtdwQ7cFPLFsX5p0Opxev57gOCrTljU054/xv7CrLOW+8KscHJoOjVszbI87GDf2R2qIsFy1ZoTeixF6L+5u0hkwEo+8ljyeuaojw/XettqH6M58nsLbyaPdL5SFM/bvxCcTlmZtT6KlF83UJJJENM32Qpo4NfpjnxzB7juV8e0Nh2eVzcU5z//Irys30r6ZsaBoncN7qf8D3wFww7F7cM3RmbR3m1IavTnheN8X07MFPd6mG03sgjsSEo7eOxMp8rh92qcF/SNn78c172Sb7JyXxWm6sZ0zJcAaWUYgzhGB23UulFwhEqbdczwlYY2+Dr9289xKMp38sw7buNXSVVmdoDQSoqwkzPezDDPfg2f0zHK9zReXx2zp97NWMXWpodCYnZ+p0ZeG3TvZhiIQ9Nswe985hP7dWvH+oIOy9lmFez7WqqZUUUqrrruzvsXltJySSUZxUexmRur78bvQ96hQCYviLWkhG/lG70e/fS9j+JJZnBJtzFLVmg2qMXNVJ3rIfNfzzFWduDB+K6MH/Eq7fmeweKbGs8mdsspZF7kkCDNc39/3b8mHU0M20ZWGlgpVu+ftQzitd0eeOi9z3ozpJvvY2avcwycc8vAwRGDUTUdl7ft1pXFMMoenCRihn62C3gyx3CRq/I5NLoG6vBZMaVq2x4dTMNt83FP7nNqlUzDnmow13S2ttnrn/IJ5nnweTYMO35UDurfKsnkbbfQW9OaaEDfzTkgMrxtRumv2MKvN/sjHhyMC8x862fU8ZtkNm1OdcZ7ZW+vvCIeMFJjDZhrzUZGQsCXN9YGg38YZsyBjZ3F7UUxbuZPJ+i7M0Ltybng4q1MuiJoILU97CPY4GH3TGg77rIQlqWxB7yePpHVplLWxWHoMm9CN1IEl4RBH1DyJwhDMuVBobOzzJ9q1KiMaXuRaJpJvHO5Cs9Iwz1zYh4tezh2grDTi4WZq0egBPpu4zC7ozXj0BdhsrTlwvahJaW3p8C4e0mHtxhpaNYlaTDfe10gjmdboO7VoRFUswfqqOCHJ9rrJ9d2UxyEREpZnK8t0Y+menO03hbpzkZSVkM/J2JtP3Mv2vTSi0SgSYn1V3JfbqFsZ049e0F0Dhzlt9M5XrCqWSGvlTtONV7gJt6ZGQhovjJjHO2MWu56nodmh3St/nrd2u0wPZn3XZqou7F79OhtPtocA+GfiDIbo/QDS3iYCEC6BnmejHTgoLeRNGkXt2rCZOrAkrBEnnFfIm5j+zmf16Uy31tkZf/xECnTyzIV9OHT3tnnLmcmfWzS229STaIRz2OgbajLW1Mhz/eKFazfxm/v/y8uj5qcnY3O1Q5SOrjKToObKU00kr/ujLQqki5A26vFuq3Ofea+tAeycvu1pQV/gpU0kVbrT8fPMuHWixmhCQ1O667ue737/+38ZR4X0ZGzKdNPUS9C7NDWsic2kWOkRO7+h2GEF/bzVGzn3hdENntmloalQZZSWltLoN+exRhk5Nr9IHsBIvSc/6D0AeD95BJDfh7uxQ9AnUoI+6iMKpRXzpY+GNe45rUfWfrdFLi9c9JucdZrD/1N7dcxZLhoS7ju9Bx9febBte97J2BRWBa8uicJNU4ZTY3S7BXNXG2aeUXPWEEvFPHHmwbXXkZmMFYG3/3gAVx25Ky0aR/Jq9FZvRPOaOl0UncfoObR901RmNZk5TTeZydjCrmdCV+ln1s/6ALcSmcnY7LhN1fFk3vwJNYnMMeurYtz9+TTWVBrzLSWeZsLsbeGQVisFp77YYU035vBr1src4Wu3JdyUj400Ip7UCWnC233f59MfpzLPkmVon+pXqE7l3zx4N/ewsiaNHA+uqQGZWrJfrB1DmUs4ZLf4KL27tKA0onnm8jS1z3+evz//meS98ElEuOjArlnbnaab7OOM/1ahlm8lZC5KwyHiyYzWlqsu0y0xGtLS53dmzbJihEDIaMl7tGvKjccbZg+nPMzW8DXLZ3Et4xSqKsc+01RmnSwNOe6vVkuN3iAzj5AP1xg6lslYp0a/1x1Dssrn4vGhs1i5oSbTcXn8ILftEU2ywkdsSXydWUQWiMgUEZkoIuNS21qJyLciMjv1Pytbr4h0EZHvRWSGiEwTkWvr+wcUKzlzWFp0l9Uqk1iiOtLSJuQBqihFR+Praw+lz87ZCZWtODV306YZCRemiVh9wMtKst0S3TT6cEjLmXwhnzZ0Vp/OgLdA0FMaff54JhZBXweNvsQxV5BrHZP5u60aYq5zCyqt0Tt/T67JV7ALZFMAOwVQVj9sOYXz+pqavLUVXhp9bUjPI/iowxljx0SJRiKR4Kgn/ufrnLZravloptY0O2GvW+R278IhLW8AtivfHM9Pc9f6amOhFNLFHKmU6q2UMp2Vbwa+U0rtDnyX+u4kAfxFKbU3cCBwlYjs41KuKLnxg0l0u/nL/AVdsA4pq2IJPhi32PYyTdG7AfB68jg+HGSYKXJpPX5eNqfmbmpAhU6e2jT60uxBo+tqypDk1PjaWsI3uNEoapzTy9xiavReL6fZeVpXotdJ0DuuZS5bcI1FozfJrdEn04I+XxOzPGhcTDfZNvrM9zP370TLJpnO2lmfm5eTsz6zTKFZlI7Ys21BIX+tNvoDdzESh4c1Ack9mnNiffdyXV6v58NtcyQkeTurr6euYOAb4/w0sWDqMpY4DXgt9fk14HRnAaXUcqXUL6nPlcAMoJOzXLHywfjcMdBzYY2ncd8X07nxw8mMnpfp7WuIMjLZgz06taVHJ0OrzxUe1W2iKmtizanRx7NXOhaKm2eCm+kmV2cy6a7j6NYmkxB7hMtiG9PsFPOwuSYJESbpKXDdTDdOYfv2z4u46OWfPdtpxen9Y9Y7YfH6rATnGY1eS3d2OSdjybhXOss5hU8uDd/8mB3bxvjetmkJfz+3t6unjomboHfWt1PT3J20F//+/W/S5/Zj9rE+Qs/9vg9fXD2AxtEwSkJo4r/T9ht3RinF9e9NdN3uJKSJr0Q8fvJE1Aa/gl4BQ0VkvIgMTG1rp5RaDoZAB7Idpi2ISDdgf8D1TRGRgSIyTkTGrV7tHZNkW6Y+41gkLAJr5QZj8sfqUx1CRyF0bJHJ8FOItwRkC3Bn6r/3xhmuYE4NvF0z/y9uk2h+QQC5oxA6V6Xu7OLJYwp6r6Qs5mRsPi8L637n6ODWT6YwcvYaut38JXM8fOtNnALQrOqZ7+dy7bt24WC10ZtnTOSajFVJFO4C0Blwy7nQyM3rxqszSEcusPrYe9jorcLN+Vy1yTMa86I0EirIG8qq6JSVhNMKUFVM9zURbzJqdiactfV3OaNb6gr3xXkuKr0m4mt00thHis/a4FfQH6KU6gOciGF+yZ2o0YGIlAEfAdcppVwTgiqlXlBK9VVK9W3bNr8b3bbGhEXr6XXPUL6eku3T/uzw3AmM3cinVQhGvBOrppvrSfKT7MFLc3dOrpkaeZdW2dmFnLhNQLmtlKyLHRcy9m2vlGyG6Ub5zvsJuc0n4xfmjiPhHB1ZO43JS8pt+0zPDqtd35ogw9l5WTV6p/bYrpk9tZ9zoGSfNHW/5lb/erDfZ+cRponKbqO3n9TptlsIbkk/vLj1pIwfvvXZrtYzWcv8MOjN8enPuc7r1fm4PTZJPf+z14YKmkYaxjPHl6BXSi1L/V8FfAL0B1aKSAeA1P/sQBPGvgiGkH9LKfVxfTR6W2TqMqP/GjUnO7nFo0NmFVyfNYmxWwq2iBgTcn8+Yrf0tlyPiB9B77WYJ2tyLaV9n7F/5xxnzPDKJfYYNG5tyRfDOx/5TDfmytjvZ7o+pmms2lguLdJrBa6J89dY63Ka2Kannp2SkJa+19ZOxtnpa5YQCM4WNikJs+Dhk9ltJyNKaW7TjbtHi6S3G///evyeWfsydZBFKE+M+Hx0aF6aTjzeM6WV++ksTuiRyWdrVSZ0VZiN3iTikRgHjGvmLeizt+tKedr0d5WlNKKacaVXck31swW30w95DUIi0gTQlFKVqc/HAfcCnwMXAw+n/n/mcqwALwMzlFJ/r8+G+2HWikp236msXmJteKGUYnVlTc4l9LUhXw7JiKZo36yxzayRS1a67XNqdO08bKleC2ByxR+xctRe7WzfC5Hpbq6SbpTmNd0IEXSufOsX1/3plbE5bPRWnGaufFircv7+LyYbo8BoOKN3WttRE9d5/SdLaj/LZKxXX2SeIlfcGXOftT3H79sufV9NxbwkHDIiY+rKJQ5O9rndRmefXXUITV0m5t346ZZMTKO/n9uLPy7fpdbmHzDuvXho9K9FHqaxVHNO7G7b9t+HvuX+yKs8Ef8qve1YbRwj9P3oKfPoEVnC4I3HABDCMKWZ96Rs02LUJ4/TkYOIETGWGsY3pV1sO8tq1qmmlLGZxyLPc3hoMh8lBwBwdFVhLp9+8XPl2wGfpF6EMPC2UmqIiIwF3heRy4FFwDkAItIReEkpdRJwCHARMEVETMPkrUqpr5wnqW+mL9vASU+P5Ppj9uDaY3b3LFdRFeOTCUt8a6dO3h27mFs+nsKFB5iBp9wfqM8mLqVft1Z0bJHf3AF2Lc5Zo4ixOlKJIx1gDp3eNUGEY1v75qWcsG97hkxbkbNcWhDUsv/0e9yCh91jjrhhmkq8hLMxGVvtebzZJKvgzLVgymuxjBd+PHgiocxkrNVGv7R8Mw9bEpEY7pWmRp97ctkr2iTYO+p5D55kPFci6Xgsbh25V3gE63Vze9Z6dcnOBOaHxtEwv+naqlbHmuhKbBr9MXvvxH9nGCO7w0OTbWX7yUzOCI3igvAwAFpULwWa0Fdm8mLUrqt+ygGUU8aEkj+xVLXmk+QA5qqOXDzuCQD+W/Ixjc1oriPgEuAg/smokmtZplrRUTLmv7NCozIVK1WYNuSDvIJeKTUP6OWyfS2QFU4wZeY5KfV5FPkCcDcQKzcYL/Uvi9bnLLdgbRXXvzeJ3l1asovV3u0T0+/VXN3o5S997bsT6dKqESP/lh0Ay41EngQSGjqIM+VbjvJuy8MdL2RSV3Rumd0ROW2uWlrjy67zi6sHeDci3c76fyRMQe9moy8Ja3lXxroumMol6C1zD7qu8o4a/aYoNEtZF1htciyX19AtXjfu9ZgCONdiKOtttbZfctxf53NkBhSzTtD7HeltKZz3/sDurdOC3skHJffavpfF19BPFmdtBxhecgPj9T1oJlU0kyr21t6x7W/sErL7p9KrAWxC/r7477kj8iYAT3R7kb/4/F2FULQhEMyH1G/sEq9JvHw4Nadc51tZ4R6rffzC9Vz08s+25dZW041blYLKEvTW97J72yb8rm9n132ZbQ5Br5T7UNxhujHLuGlupqeDk5F/OzKtTdZHujsn5mIUt/vYpCTsuTI2ltC58YNJLK8wFAOr4LSODnIJ/biPrE5+THq6ynQICds6CuuKToWWY8GUidc9crXRO3SxzIjNZS7FITH6dWvJY2fvx12n7pupVxO+vvZQXrvMGWp665BEbJOxIU04TRvFgtIL0tsiJDhTG5F1bIeqWa5CHqCFbOLoUHaY50J5OXkSTyXO5NrYn1naeI961+ahiAW9+bD6XfSSS/gsWlvF5zmW3YNVI/TXPiuXvzaWkbPXsLw8Y1qwCg9TeFndNzWlp7LbW9uQ+Q1XHbEbj57dy3Vfug7HpmRSuZZzyhKvpfO56NKqMXu2b1rwcX4xfZTd3BIbRULplbFORs1ZzQfjl7Au5Tpn1aStnfaut9qtjdb7PGGR3YvGDT8Kh65Uupz1ud1k8bs3BZYZ1MzTRu8hrEMi9OhkxERys9Fbj/GjHIgI5/TtktbsTfbu0IzD99g2vOeSyrz3iucjf+f33x2Y1qBNrg9/yN+jRmDA1ao5A2qeIkGIw1YMrtO5+1c/w1Wxa1z3/TtxKkfUGGaeJxNn85k+oMFSWRatoDcf0vq4cKc9M8o1QQNk26UKDdwEUJ5KRmA9Mm5ZRWhqlvPXbEqf0zDdeNuJnV4hfvzoveKxOK9hqI6aea65hNpijjrc3FIbR0N5Y92Y2PymN8Y8y1mv1XkvjM47IvQj6Gcs38APcwxToPWaV9VkNHqzs/JaMGXidYU1TfjoyoOZcvdxnvfPfCzcPbVy/oRtEqUUTaSGodG/cXxoHBG9mjZi9/L+c/jz9OcnEuewRLVlRmhPZ1U5+b/EdRxd85ht22qa86V+oGv5L5MHsEB1sG3L5QBQF4pX0BdousllEzezwrh5dIhj+FuXKLdWYWFdhu02KZjPRn/MPvb1a35t9G44r42WY2ifi4fP3I8Bu7Vhj/ZlBR3nB1OjdxO4uQS9s9Mxr/Wo2Ws4/8XRnudzPlfO++4cGfl5f03vG7C/8NaVtOZvaJ+a1Peq1jy92/NYEg7RtDTi7UefoyNvCLNbfeK2Eru9bjgX7KHZFzctV61Y3e23WeXXK+P5/ClshPmeqXfJKrNRmjBa39u2bVGz/ZmrOnFMzaPpbWZ+48NqnsyqY5HKXmOaK2ppXSja6JWhtKD3Vz6fOyMY2qJzAZD52LtN5hVKzKLFW003bpq2oLJWw5hC60+Hdc+KtVIXQe8M5VpIoCkrPTs3580rDshbrv8urfjb8YVpUzu3MlbLusWsb1SARm9eArcVj1bGLbAvmBo1Z3WWG6kV67X1c9Wsz5F1RbTpJqhpxqvr5RmUDh2QY4Tpdfsyk+3Z+7ZxOc+wvxyenm8x6aDcJ14n6LvR9oinGPbSBiar7lwV/pSOso5yZZgYP9qwNwNLDM37xNhDCNBR1jCq5DpAOC92B83YyC3hd3gkcR5NGzcDNjNHdeb++IXpRPUAi1Q7bo9fSqtdenPDkmuZpXemgmznj4bS6ItW0Jv4Nd34mVCLJxSpiL8ZHA9+oXLe+qJab7LVdJP1Movk1OjdmuCcRAOPyViXY5OOa1MbG30hDL60ny0Cph+6tGrMmNuOpk2TbH/rXdo0IbnYXdA7O2bz+8K1m3Kez8y7anLZ4HE2d1DnJGmhCoD1WdjsCH0BEAplr0i1kkujz7QxVdaxPZfppiHMbvXJTs1K2cmxOrijWuXau14dv5pPwho3JYyoLl1kFYPCX1CeEsCz1M4cW/Mo81SH1BS4kf8BYL1muHxuoIxbEn8EoJFFIXopme0a/GbyWE5u1AHuruB4j2CH+bztakvRmm5MAe/3Bbv+vYkMfD0TOU4plfWy1iRd8nem42UXZioysWaasZodbB4fWWYCZbzw4m6Hdw2q5HMy1o2s+OI53O/qA7egZ37YqWmpa5vuOnVfmpaWpHPGWnFGVDTvX20Cud312VTGL3R357XeEj/updbO1RpHPS3ozdADeR632uiHbhO55oKlXE2/9JBuBcVB2lJcrd+Q/vxWIuMRniRkG/k+mjiPM2vuZpbKJGOfrTqTJFOmksbcFP8jD7a+P+s8+ZKYQH5BHkzGFoib90IuFq6tYuj0lbw4wkgd9siQWexyy1e24494bDi3fDzF9Xjz+c/54rm8JBssnjRWLd4q9J1m56RurMPzWjDl1gY/ppvf9u7o2sbTettj3Nd1MjYfdY1746Q0EqJrF2CDXAAAIABJREFU22auGr3VpRUyppvaxGd57aeFnPXcj64dbaEv8Jj5GdNQdSLbdBMKpUw3eSZj/fjvOzse5dLZfTDoIB49e7+cERjvOnVffr71mLzn29KMTOzD3+NnA1BOE57p8xWHVD8F2DNu6Wj8ovYA3IPxmbyXPJKKSLaZzo82ns+LO1cwu7pQvII+dUG9XjCvy/nAVzNYs7GGV36YDxgC13wPqmJJ3hljT3jtlElWm6hzXyyRneXG6jJp9Rix2uudppNkSqMXD83X1XTj6nWTOf6M/TuxX+fs1Yt/PW4Pz1Rz9bEw5v7Te/D0+fvbthUyUnjgjB5cdsgu+QtqIUJkj8icGa3M56W0wKxaVtw0O7/Zqn4f+paussJWh7WNd0deAyAcyq3RZ0aYxvcPBh3EbSfZJw+9WmQeY70NXVo15nd9sycltwfO7NOZ95JHMEPvwtuJo/n9sf1ZSm7XT6e7qBO3EV/cR7z9fCP+QKMvkGQejT7X9Q5JJsmyrlROzdXclZ6Mtdxrt+OudrhpWlc9WnvzXIt1dD1lMXRq9DkmhN01epff4VDpc0WfrA/TzSn7deC3efLA5uLCA7py56n5c9mIFnLX6B0dr9nB1qUPq05kdyhemnUH1nJ9+EM0dEqp4f7Iq7wdfcBeX6qNYRKcEfoBAC2l0XtNtt556j7s1b4p+3QwfOb7dWvFHw/rbitjdtTOGDTm89YQK5i3Bg+f2ZOVtOLE2CMspS3NG0XSnV7bslLXY9y8d6y4vU9uHfydp9ifzXyC3M9cYW0o2slYcwIzqRSv/biAo/baiS6tMnHMc6bqk8yNTOqKkAjJ1AvVonF2ajyjvtR5rVEKXd6Tb6evTH++7ZMptnC2dvdKa/RKex1J3cu90vvFzLf4xUs7dxuqmwtlNTEmOndp04RheaJCerGlJvc8Bb1DCzNHVXXxnqqOZQt66/u9tHxz+vNNkXc4PfQjv+i7M0s3VjK3ZCNR4sQI8170PipX7M0XnEMLMhPE4Tw2+j47t2TIdbmjiXdp1YhbT9orK+F62nRTJILeTVm54tBduOLQXTzfmXwavdW8uHOrxixaV+UqpFuX2b03Ao2+njEvWEVVnLs+n5blE53rcuoqY+NM6vawAHu3b2Yrawoq8wbqCtZvirG8YnNe1863fl7Eqz8sSH//24eTWbyuCrBrB85hv64MjV40p40+9dt82uhtgt5DOzeTjrRuknlgrX703//1CF6+uK/rsX6oS6zygtBCrjHJnaYbU6O3XsMmbKY52YlG+ssM19W2bgnOdaWIEqcjayghxrWhj9hXFlClDI2yh8ynuRiCvLHU8GvpxVwR+ooDtJkcs+ETAFpKZbq+cCiS1c5CEREGHrYrHZrb4xuZz1st58S3C0Qkp2LUOM9zaX1fzBGRn/cunw0+sNEXiPmwmllhlqzfbNufKzKh1VyT1O2mG+fDn20uURz40Hcc9NCwnL2zcxIQDA+cv3wwCXBOxtrrSeqm142He6Wb372rjd76u3Jr9GNuOyYt9DPxyrF9rw3OBB0NhWhhX5OxpqC3XvL/lVzPpNKBtnL9ZCbvl9zHa5GHs+p0mm4O1SZzpfYZt4ff5MfSa7g49A3XRz7iy5Jb6SxGNrUbI+/zftQeU+X2yFsAxAnTWVbTCougj+Q23dSFjI2+ODT62mCabjp5RJu1vjuFPMP55moayr2y6E03XuTW6DMqfdIRg9vZ44qlnHGsvyTIG6vdkxoIsLEmwaKUZg/ZvyWZ0ujx0ujd6nWz0Vs7MMdcg4kp6EOaqQFlxyTfHjAEvV0Ab6xJUJPQac9aVtAKEGIJndd/WpDOINVdlqWXy2voREjQgo3clhLCh4amsiB0AXtWD6YmtchicyxpuwdvRO2dwa2RTJTDw0IZL67mUoWTH5P7cHBoOqNKrk2bdj5KDiDW4WiYsLje8h9YMZWWHVnQm6Ybr1GN1QSTyxMpa52Gxw0bd/sx3PjBJBauy34G6oOiFfT5es5ctjKlLDZ6pRzZapzHZTT/fPVaqfQS9AKnP/ODLSep87ckkkYOTPGw0fsdzlsfYvP5c77aVjezp8/rzbPD56Yf7K05WXdu3y4FmX3cbPT3fPAT3bWVjC69mtvjl/Jm8lhiSZ07P5uWLvNw5MX0559LrqKtVDBZ34X9tPm2umaVXgJATIWo/qIfzZP/x1HaRH7U96UuDE4ez8Gh6QDsqRnJ5u+N/4HbGrUEFnNuv/r3hNmznbEy9Kzf1C5HQzFgCnq3eYrf9upoSzyTyx04a37N4+VsU1ZC80aRBrPRF6+gz3HBlpVv5p7Pp+c81rx3iaTddONcouyMWulXyG6s8dLoJSvxtNP3NpE03Ss9Fkz5HM5bh59emZk6WIauJ/TowAk9OnBVKkvT1tT3Hjl7v4LKa6EQYcuCqZZs4LE5mTC1h2uTGaBNZeaGCzBy7UBXWUFMZV6RtlIBkCXkrUQlSXTlaAYzGqLwz8TpBbXTyVC9L2eVPM9HNX8C4Myau6mgjGhYY8a9JxSc6coPHVs0Kijpy/bEsxf28VXONFNqmvD0+fuzU9MSznvBmOd7+vz9+duHk9JlTeXuiD3bMnzW6vTx8aTKehdzWRpCmtZgNvqiFfROzdra61737kRmrax0HmI71mqjt3bqXjeqUI1+Q3Xcdbubkrxmoz2OfTypo0m26QYXjf7CA3bmrZ/tvv8mtg7M4wHr2Dzb/cx8eK1tNVPNbavoEcPjaifWU0kj/ltyo23/sSEjIfQJS8dyXSlUqwil4n6PCuHq8Kc593evfpP9ZB6fltzpUUJYmGxDUgkhUfyijGxpkZC25Sayi4gT9m3vq5x1rYib+6/1UU/qigl3HEvjkhB73m6kAgxrGvFkMisRUa53xHiHAht9QTgV1EaWsL1rN7knADExMnkZN3r0vLW2MLFZGn3qvyngrUK2TVlJlpA28bTR+1CTk8nUsY4ho5uN/oEzevLAGT1d67Fq9KZbofP8brk6MzFSMoUn3nVcgw0764NEEyMc7JjSq1iqWtNavDt6wLeQHxS7jnmqA+tUM8aVXgnAgJqn2KAa8XjkeY5LdSBOPkgcxgOJC9HRmKgyCd7/L3Y1/4r+E4Bf9U6A4Ud/ROzvNKMK8y47k8EE+MPv2g+3xC2H7dGWjSkFzZmJrGUTuxulqVg634hc70hZaTivW2dtKVpB79S8rTfYy2ySPtaSaelmS8iDSEiybpRzMtY6VMtlu1tV6d4B+PErT6aFsvuCKb9p66z2R9PL54z9O/HM93N5+vz9qaiKub4YaUFv2ZVvgcnWRm+aifvdSdb6Pm7K0a/z4NezeCf6ANP1rvwrcRrPRp8G4ISah5lpiYtislo1p4YoP+g9PAX9VLUL5TTN2v6FfhCxWJh2sp6PkoYffHU8yWJlX3IfzTEBGFB3MqkYMw/569aMWQ6N3kkoZHe7Nsk14r/jlH2445T8i/9qw7b9dtYBt0BgJl7atImushcYndjDGPI57efinIy1jCScuT6tTFtW4brdj0av60a92X70hU7GZptudtupaV77rNmZ1cX75sojdqV7LXL01pZEWfbwO0GIMEluSw7kgdAL6e0Xxm7h3vBgdtWWs6nTAH7Sw3Srfju9//Z4JfNV+ywh/9ua+xigTUl737yWPI5pelfaSTnPpDqHY2oepTWVjFHeYZiH6v3s7XQRJIFG37DkC8VtvSNuE6xmYD6n0uW8lw+d2dNTFtQnxSvodaegz3ze5LJy0Ypb2INwSEMplb6pn01cSueWjbP86K09dmUOQT95Se1vbjLlp50l6NMavb96zu3bJb1SN16A2SVTf+2FzU0n7FXrY2uDataFIcl+7CZLma86cLT2C0uj3egan8vKcEcOrX6SStWYstIwS/RSTo09wEE7N+JKlxf9zeSxWdsiIWFyclcmJ3e1bBXGqb04TIyJu8GJ45ijOjPHpX0XxG6lCdUue9zJ5dIXUHfyhfmwm24y20WM98MczTtfK6el4aSeHTi/f/aosL4pWkFfl3jgup4dlz2sCboyOpCkrrj23YkAXHxQ11T95nntx5mz71aq40l+zTEZnI+0jd7hXtm5peEhs9tO/jI4HbNPO9664gAufOlnXwGZTMxfsz25WYcjEQbFrwdA0Amh8330LgAS4SZp04iZtreKUipDLX3/RnPyzY0R+n7cHL+CT5OHeB7/o97D5y8xiAQafYOSttF7XGarvLYKb02EpFLpEZdT7ji1/9qEw64NRasWZAUCK2CeULdMxpqENTFi3uiKuauzl8PrLjZ6cLelzlu9iXhS2cIKmPjxTdd1d43+0N3b8uGgg7h8wC556zBJJ9UuYLbfK2HFtox1hKbQSBBmXNmRAFRE3TNDGe+gv1+ZW/AK7yaPopr6i9UeaPQNS9rrxodGb313zHkv8/445U6/bq0c56lzU31RVE9Ln/u+5dEhMwFwegtab4zz4jq/WydjTcIhIaQZgn7RWmP1WqcWjdKCeV4qcbfzxkZc/JzNJfLOgEdgFyv9ZQb/ijyNOBb66ElT0GfX3bdbq4KiSmaSahfiMWO6V24/ot6tqd+0vJArOn5MvLRV9k6MOQ+/P3FLhXIwqW1ylgB/5MqyBdiM9NZ33ix+3TG7c8SebW25HIZefxgPOjzgttTq46J6WuIJnc2pkK5OW5gp5zfHklnCeKemdl9xNxv9D3PWpgW96aESCWULgmz/ffslbst69vrkeA7VJrumy7PWNzj6KKeERtMce1q7ZNJdo68NkVT7vBZMudFQGv2YW49m3O0Nk7jCra2hkEalamxzvbWiaf5/45YagptEw9tPJ7utY947q5ec+f57TXpb3/MenZpnHdeheSMGX9qfZqWZaLd7tGtKqeNZC0w3tSAa1tJBqdwiPgKc+NSIrOOc8biVi9fNonVVhDRBVyrtc57L9TDdJseDMrb0Khqvn8VDkZcojWRf/lkrMrb7CIYtvkwyAdkaUc3ilWsAI35LXYmkBEa8EEGf+l/fyshOzUpd/fbrA7e2ihgmPi9t3NDoffpdb2FDVqDR1x+/3G5MrluFuvnJS+O+8fg96d+tFW9dcQCPnpVZpW0Kbr+WtS2l0RfVZKxN0GfZ6BXdPBLyNmsUySrrvAH/OLc3vyxaT0JX6XOERLJecKUU+8tsWkolw/T9+Z3+Na/Qjw2U0YSMwN6oGtnyVZosr6gmTIIWbCIihubeko1UUMajkec5MTQ2XVarh5c9nNboC/G6yV4Zu+2T3diQJsR1RWOP6yhixNvfmvzfkbvRu0sLflm0nmeHZxKSu5kEA2qH6fO+d4dmTFhUDmTMkl4ad/e2Zbw/6KCs7ZnggP5ejsBGXwsiIS2tbZs3zCSXGHPLsOO8Tyf0aG8x3Ri1lRJD0GlHJr9nWK/hk5K7eCX6OH1kNtfFX+SxiOGj3VvLONZtojQdp+Ta0Ef8LvQ9YERIvC/8anqVJcB/Sm5nSukVNiEPULpxSY5f5Y/2qRAHFx/czfcxaY1+O5qO7dq6cdY2Y3Jd91zYponQvFGEBQ+fzO4+PJmuOXp3njinV53a+ZwjFssJPdpzzD7tsqOKbo8hRLdRykrCvHF5f169JLN+wVT0Cr3K5ijfb9KWLTXPVXQafTypo+uK/85YadvnNKk0opo2UsEK1Zp99NkMJ5Mv1c3rpiSsEQJGqsuYM3cgEXrznw0XUDWxGXeUbuDkmgeZoXbmo3Vnpo9plVpmf3xoHL9JzOJfkX+m94VJUhZK8H70HvprswA4XfuB/tpMW/CtXMSbdvJVLhdlJeGCA1jpGUm/3RAJadx20t488NWM9DZNExJJ5SnoC3kHReCGY/dgecXm/IVz0LTUPrr0akPgdVO/HLq7PYdsoZp55rjcI4GtRVE9LdGQYbpxy9np5PnIk4wsuZ4nI8/yt8V/pg2ZBUy7fnsZ91bdzz8i/+K96L2UUoNsWkPb+BJayEb6/n97Zx4lR13t8c+t7lmyTdZJSDJJJjE7CWQjgWwSEhIImCBrWDSCGmWRAE8wEQ8PfEcFXAg+fSpH5Lg8BYHEBVFAwKc+PWJAwr4pgYQgCZt5BDJLz+/9UVXd1dXV63RPV9fczzlzurq6qvtOV9e3bt3f/d373NcYK/aFpG/CrlX+q4bPcnv9NWmfMUleSS7f2XANg+UdPtPxcV5pPZnp8hJfe35VUuQBFsaeyiryq9u+mPb8h50reGvq2Xn/z0qQDN1U5dNLZ5SviURM7EJs2QbcComfutUQXUeiu+33/ALh3jVl9vKttW+/tnA986I9+nzZOlUiWkIft2jr7OJdZ+brUmtHckAT7MHNb9fdwCx5Idnw4cSYXXp0oKRy4wftfoCFnQ9xUuxPLLCe4ZnGc+ErE/nEY6cnt7nfV/0QYK71fNrzT8R/mbHNHYml/PH1fsn4u59dXbZncUvnKm6QDwHwmhnEU6aVlW3XcX3H6Tzb1cKWzlOIxat7Q1ZL6ZUAq2emVy50B2NjlsWN62ax7YKFRb+nf55Edz25eEzSQwhZzlD16CtLqntacfsl22yG7PBELnTT3tnFe+0JDpWd/KD+On7cuYw7E0sRDPXSyXGxv3KcL9YN8Jv6TWzq+Dij5fVu2/H7xEyWxh5nkBzgD43LWHLwQd4z9Sxtu4EEMf785gDOqAeDIL7Rg+8kTuSj3M33Eyt5p3Ecb71ncX/C9hqfM2N4LjGG/0rYNc6rdXtYixOmwL4wHdYyMFl+ItFl6Ozqos4S1s7KDIMVUhzOfwy6e0wsEZZNHZ6y2fXo/TF6FfqKkvLISwzdhMwJKkjoRWQn8H9AAug0xswTkSHAbUArsBM43RjzVsC+xwE3AjHgu8aYzCabZaIhbtnt4Q68zWGWnaFwVvxBzorbA51f6Tgt6751kuCr9d8u6HNOaPsC36rbwrauJWyMb0177ZSG7/D4vxr5pVzJFGs3tw86j8t2nsq7NHIAO3TwD2NXUnyg9RKG/2MbM62dyf1vSyxL1lIZYQk/SKzKakfVhD5Z1CxcP+ZC8N6FJLqM3cSlgMJV2d8v/bGYyWpB+FNuk7Fi33ZhiwGHnctXTSmq7Ei21pqF7tfd30G5KcYtWGaMmWWMmec83wTcb4yZBNzvPE9D7Dq63wSOB6YDZ4pIZepwYt9Gz3/390y8eRpfqrs54/Xp1ksZ637ceUzW93vZCaM83ZXeru1JM56l7TdyQ+epaesvbr+Ql7uG0U4dq9qvZ9rB77G/4RD2MTgp8gCPmfexrO2rPNlyFqe0X8O0g9/jL11T+U3iCDo81958XkG1hHbMYDuDpalP7d0Qes+/hDE5Y/SFlEfKKH5Xwgl+0bJUPXr/5K0avJaGkguXTeTGdbML2nZ+6xBPmeLiPidZDC1kB647Z+pa4Ghn+fvA74DP+LaZD7xgjPkHgIjc6uyXvY9fN6iLWWx+J/sNw+rYQ2nPr+84gzsTSzgr/kDg9lu7lnCJtZVticVMs+yGzrd3Lk3b5tS2q2iR13mgazb76UdTRyr2/h6NWSe2vGhG0lgfo506oI4z2jM7DOWLgZciKuXg6jWHsmzqcA5rGZR/45DhPQEPtHWS6DJZj5FX54MOxS0fOSJjfSknuPcw+rtGSanBYqUkdvz7ShrrLH74Z9spLDaFOBYwiHvCzJFVP3yFevQGuFdEHhaRDc66EcaYVwGcx+EB+40Gdnme73bWZSAiG0Rku4hs37dvX4FmpdPamd7L8/Gu1pzb7zQj2E9mfjXA1ztPYkvnqSw6eCM3JT4AG3dw/oT7uLzzk2nbbTdT+VnXYvZjT6zZ76t1n6vYVdCEKS/5bs+rdfveWBdjVYEt2cLGCTNTDUh++/Re3jzQnj1043Hpg7z7ZVOHZwh7scdk2ZTmtDkMjXF/j4H0R6WyDOxTR0M8lsqjL3Ew1lsi4Ztnz+EbZxXWq7ZSFCr0i4wxc7BDMBeKyNJ8OzgEfU2BN8TGmJuMMfOMMfOam5uDNsnLJbsuSnv+p8FrgVRLNoCXZqS2edUM5T2nouB7h6/ns/WbuLxjA5s7PpoMy7yCY8vgVg7k7lcSiHd6/eKJw9JeCyqBAHDuolYgv8cetjhgLXDuolae+nz6uIf3e964fBKD+9b5d8uKXwiKGYS7++Il3HLufIZ6yj74PfqwhQB6C6XG6P39KcJCQUJvjNnjPO4FtmGHZF4TkZEAzuPegF13A94AdwuwpzsG56KxKzVZ5aqO9bwwYD5/TBzKpsQnuSuxAIDOgePhs3u4oP1i/mYmAcKuC16kz9otrPvw+dyeOJqfJJZjAr6ag3kalrg0D2jgfc22h+/tATln3OC07bJ59G5bvnxCXq3QTS0jIhnF5OKeDJZLj53MxuV2A27vuZrthPcLcTEX36CIUYOvtIHqfHVI5dEXGbpJevRlN6lb5BV6EeknIgPcZWAl8ATwC2C9s9l64OcBu/8VmCQi40WkHljn7Fdxbk0cQ1u/0ZzTcSVPWZN4pGsyAJ1DJ0N9P+7uOjK5rVXfBywrq/e01cmv7ttQWLXIDxw2Kvle3l6q/oqa/pPaxY0Zh3UwNmr4L5hu+MXfWyCIUo6A20IxSET84zKlhhCU7pEaGyluvxXT7d4Gg/tmliCvJoV49COAP4rIDuAh4FfGmN8A1wLHisjzwLHOc0RklIjcDWCM6QQuAu4BngZ+aox5svz/hs3BmH0CXdB+Me3UJW/B45bFLYlVrGi7nq6RszL2cwU1m3C6HZu+fOrhXP2B/ElDfeqt5InpFfoOX3OPhiyhGzcLJJ93qCl25cFfvdI9yQvKuilhBqX7eUFNpf2k9EaPdU+SLa01H59eOYWHrlxO84DKVGEtlbxZN07GTEalJmPMG8DygPV7gNWe53cDd3fPzMKwjOHmzuOT3rp7Sx6zBIPFC6YlcHDUFcxss9ncC0HzgAY+smg83/qfv/Pa/rasdtTHYskT0xu6SST8Hn3wHYJrY74fWdgmZdQqfeuDa4QXll5pPxZzp97gpFC2FVCqo9ZmH0eNYp2pmCUZ/S3CQHSm1xnDgxM+za+cWDykTkKvuAedOLE8ua/+g/2HK45hxuimrKbUxcXj0adExN8BPttgbLKDfNZPcOzSeidlwR+zdy+gBYVuShBiN2TXVkCf3lIHBZXusdPpGDdhWGH9l8NOdIRehOdHreERMzm5KqiSXJCYp0I3wW/t36c+bhHLUczCW//E69H7+7Lm8+jzORPq0ZeHfg3Bk5RyefR3nm/XIi/lCKyeYaemjvYVWQsiVdRM6Un6N9hh38WThuXZsjaovamNOfBmT1hC8uzwTogJOmHcl7N5Z0G3b7nqoHgHWb0xen9zD3e7xjqLgx2eBsOOQUP61XPRsol848FUHXuRlACFrXBSreKfjZoajA1m8cRhzB1n95ktpW75+oWtfHBOCwM9DW9+vXEJBzsyQzk6DFMdPnn0BJZMHsacsYPzb1wCY4b0Ydeb3StpXQyRkgo3e2JY/wb+vHl58kzN69HnCd0EnWy58mTr41byouENC/hDN+6Fyd9H0jsl/6JjJqa95vXi1aMvD967LvD8DrzplR4p937tpWQ+idPQxMu0kU3MDhIV5+29v5HDWwZmbqeUlYZ4rGIiD/Dby97P058/rmLv7ydSHr1b0W9Ccz9GNKUGRLzpc0HnZb6Kc0GefleO8GrMspKy4L3ITD1kQNp2rl32bMiOQHv9KZiWJckkXc26KQ8Zk5QC0iuzxesDrgmBDOtfz+vvtBdtm/vb/PDCcbzT1smqQw9hbEC3LKW2yDcrvtxEy6P3DU6agPVBWp6sT1GEbuby6C2xb80gPf573qLx3PWpxcnnrofvT7OMpV2Y/HnVZH1NKZyPeMoOZGTdFJFeWcghuPbkmWy7YFEx5qXe33lsiMe49NjJTB/VlBYOVJRCiJTQ17lBa+cEdePo3oHTnIOxRXjIuYQ+ZglfPu1wblw3iykjUl68ZQkzRqduuwc5t++nzW1J298V8MAMIRX3snD1mkOTs5f9E6aCUibTQzeZocBcR2Xd/LGMGVKaF66T4pRyECmhdz139zbb1eJ4nhh9qrRo4Z+Va65LzBKaGutYO2t0Tq97cL96nvr8Ki5cNpGffiLVUT6XHZYljBwYvjzdWmTFNHsW48A+6bMYk6GbAlz6Sgux6rxSDiJ1DxjP0nXHK5w5RbSIsyqnR1/E+7iDtfPHD0muy5VSF7OEX35qMS+9caDgz1CCueK4qaxf2Joxi9H9HXiP8OlHjOE/7sqsrl3sMIn3Dq8QNDynlINIefT12Zo8pwW2s+9fVIze59J/8YMzgz+vBPyTZIZ7hMgSYVj/hmR6n1I6MUsyGobb6+1H77X8vEWt3LrBnnHd7Kk2ma1OfVD/2e2fW8G2C4vrS6s6r5SDSAl9ckapG6N31lsBMdUgivHE/aGbsxaMZbnT67O7cXS/F/fQlSuSyxqzrTwS4NGLCAvGD+G6U2ZyzdpDM7b1smzq8MBUyWH9GzJm4eZDj7dSDqIl9Nk8+gLznrsbukmYwtIe547LnZ+ba3ftCV15khdq3zEWEc44Ymxa1kulhVgzaJVyEKkYfZ1PBZMzSD0nY67zppiTNihEn5qxmvt9bttwZPKiUKwd6uFVnnwzY73kOtQPf25FQRUqs71vl9GqlUp5iJR/GPedoG72TaGhGyni2zh+RmYrvS+dPJMz549l4fuG5rYzZuWcMJGrFLYKfeUppNZNctscQjy0fwPDm0rLkHITC/RwK+UgUkLvD5m4J2raBKQc/3ExIrp59TTuvTS9o+KoQX340skzM+4siiWXGTobtvIUM8aS6t1d3uOy1CmmpUKvlINIhW6yYeXJo0+9Vvh7xixhaL/KdJFxL1BB4qE6X3mCSiBkI+X9l7d33H+eOYdX3n63x6fKK9EkUh69i/+kKzC7suiwSLxC5SOTQh/w2sThxeVhK8Xj/g5y1TNyqVSee5/6mB6IYv3NAAAJXklEQVRrpWxESuiznXOxAmP0/qnwq2cewmNXr8y6fTUaf9xwRkazL6XMFBMeK7cnryiVIFJCnw2v15XLAYvHLHZeewIznXo0jXUxmhrrsm5fqbozuaRjQA57lPJQSntAncGqhJlIxuiTWTfJvPbUa4Wcj/nq0/u3y8cXPjiDsSUUtVLtqA7JEgjqrSsRIVJCf1jLIE6aNYqLjpmUtr7Q9EqXXW++C8BbB3LXDy9U6M9eMK6g7ZRwUEoKq14UlDATqdBNXcxiy7rZTBxuN/QttgSCy7xWe+bqo7vezrmdZsBEE3eMPVfhOkWpJSIl9H6CZqq6S18/czYAl66YjJ+vnj4LgDfyePTljsv+euMStl6w0OMdpr//+qP0zqAnKGU2qsbolTATqdBNNtK7MtmPaw4fxZrDRwVu378hzozRTZw8uyXw9UoxbWQTAPc8+U8gPUa/89oTetSW3kwxM2PrnQGgEU0NebZUlOoRaaF3J7zE0rJuCvO87vrUkorYpIQf9xdSSOBmeFMjW86YxWJnJquihJFIC71Lrd1Wa2i4uhQ72/Wk2aMraI2idJ/eEaOvLZ1PUqNmRwA3vbLKZihKmYi00LtoITBFUXozkRb6ZHqlCr1SBKXMjFWUMBNtoXfOVH8Nm/BjG15jQwuRQXRmrBIxIi30rmDWWugmVb2ytuyOCsVk3ShKLRBxobepVPGxSjHcycmeNKJ/lS3pnfRzesJOd+Y1KEqtU3B6pYjEgO3AK8aYE0XkcODbQH9gJ3C2MWZ/wH6XAh/DdpAeB841xhwsg+0FU2se/dxxQ7htw5F5m4grlaF5QAN3nn9UcgKbotQ6xXj0G4GnPc+/C2wyxswEtgGX+3cQkdHAxcA8Y8wMIAasK93c4misi6U91hILJgxN9g1Vep6544bQt75XTDNRegEFKYmItAAnYIu7yxTg987yfcApWXaPA31EJA70BfaUZmrx/NvKKVy8fBKnzOnZUgaKoihholCXcQtwBeBtrvYEsMZZPg0Y49/JGPMK8BXgZeBV4F/GmHuDPkBENojIdhHZvm/fvgLNyk3/hjiXHTuZhjr1jBVF6b3kVUARORHYa4x52PfSecCFIvIwMADIKPUoIoOBtcB4YBTQT0TOCfocY8xNxph5xph5zc3NRf4buSmlvriiKEpUKCQIuQhYIyKrgUagSUR+ZIw5B1gJICKTsUM7flYALxpj9jnbbQUWAj8qh/GFUmuDsYqiKOUkr0dvjNlsjGkxxrRiD6Q+YIw5R0SGA4iIBXwOOwPHz8vAkSLSV+xZKMtJH9DtEWotvVJRFKWcdCd4faaIPAc8gz3AeguAiIwSkbsBjDF/Ae4AHsFOrbSAm7plcQlYFQzRN9ZZNOoYgKIoIUbCOM173rx5Zvv27WV7v3fbO5l+1T1A+Rt4tHUmAGiI114Kp6Io0UFEHjbGzAt6rVckCldyMFYFXlGUsNMrYg61V9RMURSlfPQKodesG0VRejO9QuhrrZWgoihKOekVQq8oitKbUaFXFEWJOCr0iqIoEUeFXlEUJeKo0CuKokQcFXpFUZSIo0KvKIoScVToFUVRIo4KvaIoSsRRoVcURYk4KvSKoigRR4VeURQl4qjQK4qiRBwVekVRlIijQq8oihJxVOgVRVEijgq9oihKxFGhVxRFiTjxahvQU3znQ3OJaUtBRVF6Ib1G6Fcdeki1TVAURakKGrpRFEWJOCr0iqIoEUeFXlEUJeKo0CuKokQcFXpFUZSIo0KvKIoScVToFUVRIo4KvaIoSsQRY0y1bchARPYBL5W4+zDg9TKaUynUzvJTK7aqneWnVmytpJ3jjDHNQS+EUui7g4hsN8bMq7Yd+VA7y0+t2Kp2lp9asbVadmroRlEUJeKo0CuKokScKAr9TdU2oEDUzvJTK7aqneWnVmytip2Ri9EriqIo6UTRo1cURVE8qNAriqJEnMgIvYgcJyLPisgLIrIpBPZ8T0T2isgTnnVDROQ+EXneeRzseW2zY/uzIrKqB+0cIyIPisjTIvKkiGwMo60i0igiD4nIDsfOa8Jop+ezYyLyNxG5K+R27hSRx0XkURHZHlZbRWSQiNwhIs84v9WjwmaniExxvkf3b7+IXBIKO40xNf8HxIC/AxOAemAHML3KNi0F5gBPeNZdD2xyljcB1znL0x2bG4Dxzv8S6yE7RwJznOUBwHOOPaGyFRCgv7NcB/wFODJsdnrsvQz4MXBXWI+98/k7gWG+daGzFfg+8DFnuR4YFEY7PfbGgH8C48JgZ4/94xX+Uo8C7vE83wxsDoFdraQL/bPASGd5JPBskL3APcBRVbL558CxYbYV6As8AiwIo51AC3A/cIxH6ENnp/N5QUIfKluBJuBFnOSRsNrps20l8L9hsTMqoZvRwC7P893OurAxwhjzKoDzONxZHwr7RaQVmI3tLYfOVicc8iiwF7jPGBNKO4EtwBVAl2ddGO0EMMC9IvKwiGxw1oXN1gnAPuAWJxz2XRHpF0I7vawDfuIsV93OqAi9BKyrpbzRqtsvIv2BO4FLjDH7c20asK5HbDXGJIwxs7A95vkiMiPH5lWxU0ROBPYaYx4udJeAdT157BcZY+YAxwMXisjSHNtWy9Y4dhj0W8aY2cAB7BBINqr6nYpIPbAGuD3fpgHrKmJnVIR+NzDG87wF2FMlW3LxmoiMBHAe9zrrq2q/iNRhi/x/G2O2htlWAGPM28DvgOMIn52LgDUishO4FThGRH4UQjsBMMbscR73AtuA+SG0dTew27mDA7gDW/jDZqfL8cAjxpjXnOdVtzMqQv9XYJKIjHeupuuAX1TZpiB+Aax3ltdjx8Pd9etEpEFExgOTgId6wiAREeBm4GljzNfCaquINIvIIGe5D7ACeCZsdhpjNhtjWowxrdi/wweMMeeEzU4AEeknIgPcZey48hNhs9UY809gl4hMcVYtB54Km50eziQVtnHtqa6dPTlAUeHBj9XYGSN/B64MgT0/AV4FOrCv3B8FhmIP0j3vPA7xbH+lY/uzwPE9aOdi7NvFx4BHnb/VYbMVOAz4m2PnE8BVzvpQ2emz+WhSg7GhsxM79r3D+XvSPW9CaussYLtz/H8GDA6pnX2BN4CBnnVVt1NLICiKokScqIRuFEVRlCyo0CuKokQcFXpFUZSIo0KvKIoScVToFUVRIo4KvaIoSsRRoVcURYk4/w/AxA0sd5p2ngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYJklEQVR4nO3df5BdZ33f8ffHKxsbByNcrxmhFZFThFvFAVe9lUUhA1FwKsmOlB8F7ODKOBk0IqhjIMHIkBLchI5LmsbxoFoYV8WqAeHOQLpRRBVPPS5Dg2xd1T+CsA1b1UQbCbQELHBsbGR/+sd9NFyur7RntWvtrp7Pa+bM3vM833Puc2ztfu4599z7yDYREVGf06Z7ABERMT0SABERlUoARERUKgEQEVGpBEBERKUSABERlUoARLwAJG2W9G+mexwRx5MAiBlN0m9Iakt6QtJBSV+U9AZJV0p6TJJ66udIOiTpcklvkvRc2fYHkh6VdM0LMMZ3SPpyd5vt9bb/YKqfqzyfJb3qhdh31CUBEDOWpPcBNwH/Dng58ErgPwFrgC8Ac4E39my2AjDwP8r6Ads/BZwDvBf4pKQLX/jRR8x8CYCYkSS9FPi3wLttf97239v+ke0/t/1+2z8E7gTW9my6Fvi07SPdje7YAXwXeM1xnneZpL+S9LikByW9qavvHZL2lbOJ/yfp7ZL+MbAZeF0503i81H5K0h+Wx2+SNCrpunJ2clDSr0haJenrkr4r6YNdz7NU0lfKGA5K+rikM0rfl0rZg+X53lbaL5f0QNnmryS9pmt/H5D0t11nQb84kf8XcQqznSXLjFvovJI/Asw5Ts3rge8DZ5X1lwJPAReX9TcBo+XxacBq4Dngnxxjf/OBvwNWlfpLy/ogcHZ5rgtL7TzgZ8vjdwBf7tnXp4A/7BrHEeDDwOnAO4Ex4DPAS4CfBX4I/Eyp/6fAMmAOsBB4GHhP174NvKprfQlwCLgEGACuBh4DXgRcCOwHXlFqFwL/cLr//2aZGUvOAGKm+gfAd9zzSr6b7f8NfBv41dL0VuDrth/oKntFeVX+FJ3LRu+zff8xdnkVsMP2DtvP2b4LaNMJBOiEx0WSzrJ90PbeCRzPj4CP2v4RsA04D/hT2z8o+9lLOTOxvcf2LttHbD8GfILnX+rq9k7gE7bvtf2s7duBp+mEyLN0gmCxpNNtP2b7/05g3HEKSwDETPV3wHmS5oxTt5UfXwb6V8DtPf0HbM+l8x7AzcDy4+zrp4G3lMsoj5fgeAMwz/bfA28D1gMHJf2FpH80keOx/Wx5/FT5+e2u/qeAnwKQ9GpJ2yV9S9L36bwHct444/6dnnEvoPOqfwR4D/AR4JCkbZJeMYFxxyksARAz1VfoXBb5lXHqtgK/KOl1dF7xfqZfke2ngQ8APyfpWPvcD/xX23O7lrNt31j2sdP2pXQu/zwCfPLo7idyYA3cUva/yPY5wAcBHad+P52zi+5xv9j2Z8u4P2P7DXSCwsC/n+LxxiyVAIgZyfZhOtfMN5U3TF8s6XRJKyV9rKvum8CXgc8Cd9n+1nH2+Qzwx2W//dwB/LKkfyFpQNKZ5Q3cIUkvl7Ra0tl0Lq88QefyCnReyQ8dfaN2CryEzvsNT5SzjHf19H8b+Jmu9U8C6yVdoo6zJV0m6SWSLpS0XNKL6ATqU13jjsolAGLGsv0fgfcBv0fnTdP9wAbgz3pKb6fz6nZrg91uAV4p6Zf7PN9+OreYfrDr+d5P5/fkNOB3gAN07iR6I/DbZdO76VzD/5ak7zQ/wmP6XeA3gB/Q+eP+uZ7+jwC3l8s9b7XdpvM+wMeB7wEjdN6Yhs71/xuB7wDfAs4vxxeB7EwIExFRo5wBRERUKgEQEVGpBEBERKUSABERlRrvQzYzynnnneeFCxdO9zAiImaVPXv2fMf2YG/7rAqAhQsX0m63p3sYERGziqRv9mvPJaCIiEolACIiKpUAiIioVKMAkLSiTCQxImljn35Jurn0PyRpSU//gKT7JW3vs+3vlinujvdthxERMcXGDQBJA8AmYCWwGLhS0uKespXAorKso/Ntht2upTOpRe++F9CZdONvJjzyiIiYlCZnAEuBEdv7yrcpbqPzhVnd1gBb3bELmCtpHoCkIeAy4LY++/4T4Dqm/ut0IyJiHE0CYD6db0U8arS0Na25ic4f+ee6N5C0Gvhb2w8e78klrZPUltQeGxtrMNyIiGiiSQD0m4ii9xV73xpJlwOHbO/5iWLpxcCHOPb3sv94J/attlu2W4ODz/scQ0REnKAmHwQbpTO93FFDdL4TvUnNvwRWS1oFnAmcI+kOOjMSXQA8KOlo/f+RtPR4E3pEnCzl3+ULLl/HHtOpyRnAbmCRpAvKjEdXAMM9NcPA2nI30DLgcJk0+3rbQ7YXlu3utn2V7b+2fb7thaVvFFiSP/4xU9ie0HIi2+SPf0y3cc8AbB+RtAHYCQwAW2zvlbS+9G8GdgCr6MxE9CRwzQs35IiImAqzakawVqvlfBdQzESS8oo+ZixJe2y3etvzSeCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVAIgIqJSjQJA0gpJj0oakbSxT78k3Vz6H5K0pKd/QNL9krZ3tf1BqX1A0l9KesXkDyciIpoaNwAkDQCbgJXAYuBKSYt7ylYCi8qyDrilp/9a4OGetj+y/RrbFwPbgQ9PfPgREXGimpwBLAVGbO+z/QywDVjTU7MG2OqOXcBcSfMAJA0BlwG3dW9g+/tdq2cDmU8vIuIkahIA84H9Xeujpa1pzU3AdcBzvTuW9FFJ+4G3c4wzAEnrJLUltcfGxhoMNyIimmgSAOrT1vtqvW+NpMuBQ7b39Nux7Q/ZXgB8GthwjJpbbbdstwYHBxsMNyIimmgSAKPAgq71IeBAw5rXA6slPUbn0tFySXf0eY7PAL/ecMwRETEFmgTAbmCRpAsknQFcAQz31AwDa8vdQMuAw7YP2r7e9pDthWW7u21fBSBpUdf2q4FHJnswERHR3JzxCmwfkbQB2AkMAFts75W0vvRvBnYAq4AR4EngmgbPfaOkC+m8N/BNYP2JHUJERJwI2bPn5ptWq+V2uz3dw4h4HknMpt+lqIukPbZbve35JHBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpRgEgaYWkRyWNSNrYp1+Sbi79D0la0tM/IOl+Sdu72v5I0iOl/guS5k7+cCIioqlxA0DSALAJWAksBq6UtLinbCWwqCzrgFt6+q8FHu5puwu4yPZrgK8D10949BERccKanAEsBUZs77P9DLANWNNTswbY6o5dwFxJ8wAkDQGXAbd1b2D7L20fKau7gKFJHEdERExQkwCYD+zvWh8tbU1rbgKuozP5+7H8JvDFfh2S1klqS2qPjY01GG5ERDTRJADUp6139uu+NZIuBw7Z3nPMnUsfAo4An+7Xb/tW2y3brcHBwQbDjYiIJpoEwCiwoGt9CDjQsOb1wGpJj9G5dLRc0h1HiyRdDVwOvN12b6hERMQLqEkA7AYWSbpA0hnAFcBwT80wsLbcDbQMOGz7oO3rbQ/ZXli2u9v2VdC5swj4ALDa9pNTdUAREdHMnPEKbB+RtAHYCQwAW2zvlbS+9G8GdgCrgBHgSeCaBs/9ceBFwF2SAHbZXn9CRxEREROm2XTlpdVqud1uT/cwIp5HErPpdynqImmP7VZvez4JHBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVCoBEBFRqUYBIGmFpEcljUja2Kdfkm4u/Q9JWtLTPyDpfknbu9reImmvpOckPW+mmoipcu655yLpBV2AF/w5zj333Gn+LxmnmnHnBJY0AGwCLgVGgd2Shm1/ratsJbCoLJcAt5SfR10LPAyc09X2VeDXgE9M5gAixvO9733vlJiu8WjQREyVJmcAS4ER2/tsPwNsA9b01KwBtrpjFzBX0jwASUPAZcBt3RvYftj2o5M+goiIOCFNAmA+sL9rfbS0Na25CbgOeO5EBihpnaS2pPbY2NiJ7CIiIvpoEgD9zjt7z6f71ki6HDhke8+ER3Z0J/attlu2W4ODgye6m4iI6NEkAEaBBV3rQ8CBhjWvB1ZLeozOpaPlku444dFGRMSUaRIAu4FFki6QdAZwBTDcUzMMrC13Ay0DDts+aPt620O2F5bt7rZ91VQeQEREnJhxA8D2EWADsJPOnTx32t4rab2k9aVsB7APGAE+Cfz2ePuV9KuSRoHXAX8haecJHkNERJwAzabb41qtltvt9nQPI2YZSafMbaCnwnHEySdpj+3nfd4qnwSOiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIq1SgAJK2Q9KikEUkb+/RL0s2l/yFJS3r6ByTdL2l7V9u5ku6S9I3y82WTP5yIiGhqzngFkgaATcCldCZ/3y1p2PbXuspWAovKcglwS/l51LV0ppM8p6ttI/A/bd9YQmUj8IFJHEtEX/79c+AjL53uYUyaf/+c8YsiJmDcAACWAiO29wFI2gasAboDYA2w1Z356nZJmitpnu2DkoaAy4CPAu/r2eZN5fHtwD0kAOIFoBu+f0pMpSgJf2S6RxGnkiaXgOYD+7vWR0tb05qbgOuA53q2ebntgwDl5/n9nlzSOkltSe2xsbEGw42IiCaaBID6tPW+nOpbI+ly4JDtPRMe2dGd2LfabtluDQ4OnuhuIiKiR5MAGAUWdK0PAQca1rweWC3pMWAbsFzSHaXm25LmAZSfhyY8+oiIOGFNAmA3sEjSBZLOAK4AhntqhoG15W6gZcBh2wdtX297yPbCst3dtq/q2ubq8vhq4L9P9mAiIqK5cd8Etn1E0gZgJzAAbLG9V9L60r8Z2AGsAkaAJ4FrGjz3jcCdkn4L+BvgLSd2CBERcSI0m+6OaLVabrfb0z2MmGUknTp3AZ0CxxEnn6Q9tlu97fkkcEREpRIAERGVSgBERFQqARARUakEQEREpRIAERGVSgBERFQqARARUakEQEREpRIAERGVSgBERFQqARARUakEQEREpRIAERGVSgBERFQqARARUalGASBphaRHJY1I2tinX5JuLv0PSVpS2s+UdJ+kByXtlXRD1zavlfQVSX8t6c8lnTN1hxUREeMZNwAkDQCbgJXAYuBKSYt7ylYCi8qyDriltD8NLLf9WuBiYEWZMxjgNmCj7Z8DvgC8f5LHEhERE9DkDGApMGJ7n+1ngG3Amp6aNcBWd+wC5kqaV9afKDWnl+XonHYXAl8qj+8Cfn0yBxIRERPTJADmA/u71kdLW6MaSQOSHgAOAXfZvrfUfBVYXR6/BVjQ78klrZPUltQeGxtrMNyIiGiiSQCoT1vvzNTHrLH9rO2LgSFgqaSLSv9vAu+WtAd4CfBMvye3favtlu3W4OBgg+FGREQTcxrUjPKTr86HgAMTrbH9uKR7gBXAV20/AvwSgKRXA5dNaOQRETEpTc4AdgOLJF0g6QzgCmC4p2YYWFvuBloGHLZ9UNKgpLkAks4C3gw8UtbPLz9PA34P2DwlRxQREY2MewZg+4ikDcBOYADYYnuvpPWlfzOwA1gFjABPAteUzecBt5c7iU4D7rS9vfRdKend5fHngf8yRccUERENyO69nD9ztVott9vt6R5GzDKSmE3/zo/lVDmOOPkk7bHd6m3PJ4EjIirV5E3giFlP6nej2uzyspe9bLqHEKeYBECc8k7GZZNcnonZKJeAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUowCQtELSo5JGJG3s0y9JN5f+hyQtKe1nSrpP0oOS9kq6oWubiyXtkvSApLakpVN3WBERMZ5xA6DM57sJWAkspjOX7+KespXAorKsA24p7U8Dy22/FrgYWFEmjQf4GHCD7YuBD5f1iIg4SZqcASwFRmzvs/0MsA1Y01OzBtjqjl3AXEnzyvoTpeb0shydNcPAOeXxS4EDkzmQiIiYmCYzgs0H9netjwKXNKiZDxwsZxB7gFcBm2zfW2reA+yU9B/oBNE/7/fkktbROavgla98ZYPhRkREE03OAPpNpto7990xa2w/Wy7zDAFLJV1U+t8FvNf2AuC9wH/u9+S2b7Xdst0aHBxsMNyIiGiiSQCMAgu61od4/uWacWtsPw7cA6woTVcDny+P/xudS00REXGSNAmA3cAiSRdIOgO4AhjuqRkG1pa7gZYBh20flDQoaS6ApLOANwOPlG0OAG8sj5cD35jksURExASM+x6A7SOSNgA7gQFgi+29ktaX/s3ADmAVMAI8CVxTNp8H3F7eBzgNuNP29tL3TuBPJc0Bfki5zh8RESeH7N7L+TNXq9Vyu92e7mFEPI8kZtPvUtRF0h7brd72fBI4IqJSCYCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVKMAkLRC0qOSRiRt7NMvSTeX/ockLSntZ0q6T9KDkvZKuqFrm89JeqAsj0l6YOoOKyIixjPulJBlOsdNwKV0Jn/fLWnY9te6ylYCi8pyCXBL+fk0sNz2E5JOB74s6Yu2d9l+W9dz/DFweKoOKiIixtfkDGApMGJ7n+1ngG3Amp6aNcBWd+wC5kqaV9afKDWnl+Un5s2TJOCtwGcncyARETExTQJgPrC/a320tDWqkTRQLu8cAu6yfW/Ptj8PfNv2N/o9uaR1ktqS2mNjYw2GGxERTTQJAPVp6539+pg1tp+1fTEwBCyVdFFP3ZUc59W/7Vttt2y3BgcHGww3IiKaaBIAo8CCrvUh4MBEa2w/DtwDrDjaJmkO8GvA5xqPOCIipkSTANgNLJJ0gaQzgCuA4Z6aYWBtuRtoGXDY9kFJg5LmAkg6C3gz8EjXdm8GHrE9OukjiYiICRn3LiDbRyRtAHYCA8AW23slrS/9m4EdwCpgBHgSuKZsPg+4vdxJdBpwp+3tXbu/grz5GxExLWT3Xs6fuVqtltvt9nQPI+J5JDGbfpeiLpL22G71tueTwBERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUaBYCkFZIelTQiaWOffkm6ufQ/JGlJaT9T0n2SHpS0V9INPdv967LfvZI+NjWHFBERTYw7JWSZznETcCmdyd93Sxq2/bWuspXAorJcAtxSfj4NLLf9hKTTgS9L+qLtXZJ+AVgDvMb205LOn9Iji5gESSdlm8wiFtNp3AAAlgIjtvcBSNpG5w93dwCsAba68695l6S5kubZPgg8UWpOL8vRf/HvAm60/TSA7UOTPpqIKZI/zFGDJpeA5gP7u9ZHS1ujGkkDkh4ADgF32b631Lwa+HlJ90r6X5L+Wb8nl7ROUltSe2xsrMFwIyKiiSYB0O+8tvfl0TFrbD9r+2JgCFgq6aLSPwd4GbAMeD9wp/qcQ9u+1XbLdmtwcLDBcCMiookmATAKLOhaHwIOTLTG9uPAPcCKrm0+7477gOeA8xqPPCIiJqVJAOwGFkm6QNIZwBXAcE/NMLC23A20DDhs+6CkQUlzASSdBbwZeKRs82fA8tL3auAM4DuTPqKIiGhk3DeBbR+RtAHYCQwAW2zvlbS+9G8GdgCrgBHgSeCasvk84PZyJ9FpwJ22t5e+LcAWSV8FngGudt55i4g4aTSb/ua2Wi232+3pHkZExKwiaY/tVm97PgkcEVGpBEBERKVm1SUgSWPAN6d7HBF9nEduYoiZ66dtP+8++lkVABEzlaR2v2usETNZLgFFRFQqARARUakEQMTUuHW6BxAxUXkPICKiUjkDiIioVAIgIqJSCYCISZC0RdKh8p1WEbNKAiBicj7Fj7/iPGJWSQBETILtLwHfne5xRJyIBEBERKUSABERlUoARERUKgEQEVGpBEDEJEj6LPAV4EJJo5J+a7rHFNFUvgoiIqJSOQOIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISv1/ywNvob2GrkoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "s0, cvr, delay = 50, 0.04, 20\n",
    "n, sigma, drift = 100, 0.1, 0.0005\n",
    "\n",
    "# generate signals\n",
    "s = gen_cvr_signal(etco2_mss3, s0, cvr, delay = delay, sigma = sigma, drift = drift, n = n)\n",
    "\n",
    "# plot signals\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(s[:,0]);\n",
    "ax.plot(np.mean(s, axis=1));\n",
    "ax.set_title('First and mean signal');\n",
    "\n",
    "# fit sigals\n",
    "results = fit_cvr_delay(etco2_mss3, s, min_delay=-100, max_delay=100, baseline=list(range(10)))\n",
    "\n",
    "# plot results\n",
    "fig, ax = plt.subplots()\n",
    "bp = ax.boxplot(results['cvr']);\n",
    "ax.set_title('CVR estimates');\n",
    "print(\"mean CVR: {mean}, std CVR: {std}\".format(mean=np.mean(results['cvr']), std=np.std(results['cvr'])))\n",
    "print(\"mean delay: {mean}, std delay: {std}\".format(mean=np.mean(results['delay']), std=np.std(results['delay'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interim-charge",
   "metadata": {},
   "source": [
    "### Interactive simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "noted-creature",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "881d5a1eb90c43959fad116b95b3e9ec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=50, description='s0'), FloatSlider(value=0.04, description='cvr', max=0.…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# interactive plot for fitting results\n",
    "%matplotlib widget\n",
    "\n",
    "etco2, n = etco2_mss3, 100\n",
    "\n",
    "@widgets.interact(s0 = (0, 100, 1), cvr = (0, 0.1, 0.005), delay = (-30, 30, 1), sigma = (0, 0.5, 0.025), drift = (0, 0.001, 0.0001))\n",
    "def fit_signals(s0 = 50, cvr = 0.04 , delay = 20, sigma = 0.1, drift = 0.0005):\n",
    "    # generate signals\n",
    "    s = gen_cvr_signal(etco2, s0, cvr, delay = delay, sigma = sigma, drift = drift, n = n)\n",
    "    \n",
    "    # plot signals\n",
    "    if plt.fignum_exists(10):\n",
    "        plt.close(10)    \n",
    "    fig, ax = plt.subplots(1,2)\n",
    "    ax[0].plot(s[:,0]);\n",
    "    ax[0].plot(np.mean(s, axis=1));\n",
    "    ax[0].set_title('First and mean signal');\n",
    "    \n",
    "    # fit signals\n",
    "    results = fit_cvr_delay(etco2_mss3, s, min_delay=-40, max_delay=40, baseline=list(range(10)))\n",
    "    #print(\"mean CVR: {mean}, std CVR: {std}\".format(mean=np.mean(results['cvr']), std=np.std(results['cvr'])))\n",
    "    #print(\"mean delay: {mean}, std delay: {std}\".format(mean=np.mean(results['delay']), std=np.std(results['delay'])))\n",
    "    ax[1].set_title('CVR estimates');\n",
    "    ax[1].set_ylim([0,0.1])\n",
    "    ax[1].boxplot(results['cvr']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "average-sphere",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close('all')"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}