{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NOTEBOOK_HEADER-->\n",
    "*This notebook contains material from [cbe67701-uncertainty-quantification](https://ndcbe.github.io/cbe67701-uncertainty-quantification);\n",
    "content is available [on Github](https://github.com/ndcbe/cbe67701-uncertainty-quantification.git).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [4.0 Local Sensitivity Analysis Based on Derivative Approximations](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.00-Local-Sensitivity-Analysis-Based-on-Derivative-Approximations.html) | [Contents](toc.html) | [4.2 Advection-Diffusion-Reaction (ADR) Example](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html)<p><a href=\"https://colab.research.google.com/github/ndcbe/cbe67701-uncertainty-quantification/blob/master/docs/04.01-Local-sensitivity-analysis.ipynb\"> <img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.ipynb\"> <img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[4.1 Local sensitivity analysis and difference approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1-Local-sensitivity-analysis-and-difference-approximation)",
     "section": "4.1 Local sensitivity analysis and difference approximation"
    }
   },
   "source": [
    "# 4.1 Local sensitivity analysis and difference approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[4.1 Local sensitivity analysis and difference approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1-Local-sensitivity-analysis-and-difference-approximation)",
     "section": "4.1 Local sensitivity analysis and difference approximation"
    }
   },
   "source": [
    "Created by Jialu Wang (jwang44@nd.edu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[4.1 Local sensitivity analysis and difference approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1-Local-sensitivity-analysis-and-difference-approximation)",
     "section": "4.1 Local sensitivity analysis and difference approximation"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.sparse as sparse\n",
    "import scipy.sparse.linalg as linalg\n",
    "import scipy.integrate as integrate\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[4.1 Local sensitivity analysis and difference approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1-Local-sensitivity-analysis-and-difference-approximation)",
     "section": "4.1 Local sensitivity analysis and difference approximation"
    }
   },
   "source": [
    "This example was adapted from code from:\n",
    "\n",
    "\n",
    "McClarren, Ryan G (2018). *Uncertainty Quantification and Predictive Computational Science: A Foundation for Physical Scientists and Engineers*, *Chapter 4: Local Sensitivity Analysis Based on Derivative Approximations*, Springer, https://doi.org/10.1007/978-3-319-99525-0_4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.1 Simple advection-diffusion-reaction (ADR) example ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.1-Simple-advection-diffusion-reaction-(ADR)-example)",
     "section": "4.1.1 Simple advection-diffusion-reaction (ADR) example "
    }
   },
   "source": [
    "## 4.1.1 Simple advection-diffusion-reaction (ADR) example \n",
    "\n",
    "The steady advection-diffusion-reaction(ADR) equation in one-spatial dimension with a spatially constant, but uncertain, diffusion coefficient, a linear reaction term, and a prescribed uncertain source.\n",
    "\n",
    "$$\n",
    "\\nu \\frac{du}{dx} - \\omega \\frac{d^2u}{dx^2} + \\kappa(x)u = S(x)\n",
    "$$\n",
    "\n",
    "The QoI is total reaction rate: \n",
    "$$\n",
    "Q= \\int^{10}_{0} \\kappa(x_ u)u(x)dx\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.1 Simple advection-diffusion-reaction (ADR) example ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.1-Simple-advection-diffusion-reaction-(ADR)-example)",
     "section": "4.1.1 Simple advection-diffusion-reaction (ADR) example "
    }
   },
   "outputs": [],
   "source": [
    "def ADRSource(Lx, Nx, Source, omega, v, kappa):\n",
    "    '''\n",
    "    Solve the ADR example with algorithm 4.1 from the textbook\n",
    "    \n",
    "    Arguments:\n",
    "        Lx: the scale of x is [0, Lx]\n",
    "        Nx: the number of points in x scale\n",
    "        Source: S(x)\n",
    "        omega: the diffusion coefficient\n",
    "        v: the convection term constant\n",
    "        kappa: the linear reaction term coefficient\n",
    "        \n",
    "    Return:\n",
    "        Solution: u \n",
    "        Q: QoI, the total reaction rate \n",
    "    '''\n",
    "    #Solves the diffusion equation\n",
    "    A = sparse.dia_matrix((Nx,Nx),dtype=\"complex\")\n",
    "    dx = Lx/Nx\n",
    "    i2dx2 = 1.0/(dx*dx)\n",
    "    \n",
    "    #fill diagonal of A\n",
    "    A.setdiag(2*i2dx2*omega + np.sign(v)*v/dx + kappa)\n",
    "    \n",
    "    #fill off diagonals of A\n",
    "    A.setdiag(-i2dx2*omega[1:Nx] + \n",
    "              0.5*(1-np.sign(v[1:Nx]))*v[1:Nx]/dx,1)\n",
    "    \n",
    "    A.setdiag(-i2dx2*omega[0:(Nx-1)] - \n",
    "              0.5*(np.sign(v[0:(Nx-1)])+1)*v[0:(Nx-1)]/dx,-1)\n",
    "    \n",
    "    #solve A x = Source\n",
    "    Solution = linalg.spsolve(A,Source)\n",
    "    Q = integrate.trapz(Solution*kappa,dx=dx)\n",
    "    return Solution, Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.1 Simple advection-diffusion-reaction (ADR) example ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.1-Simple-advection-diffusion-reaction-(ADR)-example)",
     "section": "4.1.1 Simple advection-diffusion-reaction (ADR) example "
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(52.390252927707564+0j)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format\n",
      "  SparseEfficiencyWarning)\n"
     ]
    }
   ],
   "source": [
    "# q: mean, variance, function\n",
    "q_nom = 1\n",
    "q_var = 7.062353e-4\n",
    "Source_func = lambda x, q:  q*x*(10-x)\n",
    "\n",
    "# kappa: mean, variance, function\n",
    "kappal_nom = 0.1\n",
    "kappal_var = 8.511570e-6\n",
    "kappah_nom = 2\n",
    "kappah_var = 0.002778142\n",
    "kappa_func = lambda x, kappal, kappah: kappah + (kappal-kappah)*(x>5)*(x<7.5)\n",
    "\n",
    "# v: mean, variance, function\n",
    "v_nom = 10\n",
    "v_var = 0.0723493\n",
    "v_func = lambda x,v: v*np.ones(x.size)\n",
    "\n",
    "# omega: mean, variance, function\n",
    "omega_nom = 20\n",
    "omega_var = 0.3195214\n",
    "omega_func = lambda x,omega: omega*np.ones(x.size)\n",
    "\n",
    "\n",
    "# define the scale of x \n",
    "Lx = 10\n",
    "Nx = 2000\n",
    "dx = Lx/Nx\n",
    "xs = np.linspace(dx/2,Lx-dx/2,Nx)\n",
    "\n",
    "# Calculating Q with the mean value of parameters\n",
    "source = Source_func(xs, 1)\n",
    "kappa = kappa_func(xs, 0.1, 2)\n",
    "vs = v_func(xs, v_nom)\n",
    "omegas = omega_func(xs, omega_nom)\n",
    "\n",
    "sol,Q = ADRSource(Lx, Nx, source, omegas, vs, kappa)\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.2 First-order sensitivity approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.2-First-order-sensitivity-approximation)",
     "section": "4.1.2 First-order sensitivity approximation"
    }
   },
   "source": [
    "## 4.1.2 First-order sensitivity approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.2 First-order sensitivity approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.2-First-order-sensitivity-approximation)",
     "section": "4.1.2 First-order sensitivity approximation"
    }
   },
   "source": [
    "Forward difference derivative approximation:\n",
    "\n",
    "$\\frac{\\partial Q}{\\partial x_i} \\bigg|_{\\overline{x}} \\approx \\frac{Q(\\overline{x} + \\delta_i\\hat{e_i}) - Q(\\overline{x})}{\\delta_i}$\n",
    "\n",
    "Central difference derivative approximation:\n",
    "\n",
    "$\\frac{\\partial Q}{\\partial x_i} \\bigg|_{\\overline{x}} \\approx \\frac{Q(\\overline{x} + 0.5\\delta_i\\hat{e_i}) - Q(\\overline{x} - 0.5\\delta_i\\hat{e_i})}{\\delta_i}$\n",
    "\n",
    "The complex step formula: \n",
    "\n",
    "$\\frac{\\partial Q}{\\partial x_i} \\bigg|_{\\overline{x}} \\approx \\frac{\\mathscr{I} \\{Q(\\overline{x} + i\\delta_i \\hat{e_i})\\}}{\\delta_i} + O(\\delta_i^2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.2 First-order sensitivity approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.2-First-order-sensitivity-approximation)",
     "section": "4.1.2 First-order sensitivity approximation"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:27: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:39: ComplexWarning: Casting complex values to real discards the imaginary part\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEcCAYAAADtODJSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd8VFXawPHfk0IChJqE3lsCgRA6KE1QwAYWsLu2fRVdV11dFXV1XVdXUde1iwUWV11FEVxUBESaICIgvRNAOgRCJyHtef+4kzCEmWSSzGQS8nw/nzEz9557znNDzJNzz7nniqpijDHGFFdIsAMwxhhTvlkiMcYYUyKWSIwxxpSIJRJjjDElYonEGGNMiVgiMcYYUyKWSIwxxpSIJRJjjDElYonEGEBE1ohI/+LuL27ZkhCROBFZJiLHROS+QLdnjDeWSEyZJCK9ReQnETkiIqkiskBEugWqPVVNUNU5bu1vE5ELve33ta789fjZI8AcVa2mqq97KiAit4rIKhE5KSJ7ReRtEanhrcIAx2vOUZZITJkjItWBb4A3gNpAQ+BvwKlgxlUGNQXWeNspIg8Bo4GHgRpAT6AZMENEwksjQFNBqKq97FWmXkBX4HAB+xsAXwIpwFbgvnz7twF/BlYCR4AJQKRr36PALuAYsAEY6HbMha73HwE5QBpwHOcvf/f9o4CJ+dp8DXjdvS4v9TwMfJnv2DeAV72ca1tgDnAYJ2kMdW2fBWQD6a662+Q7rrpr+zX5tkcB+4FbPLR1VrwFxeAl3uHAz673dYCFwIPB/pmyV2BfQQ/AXvbK/3L9EjwIfAhcDNRy2xcCLAWeAioBLYAtwGC3MtuAX1wJpzawDhgJxAE7gAaucs2Alm7HXJivDo+fcXoCJ4Hqrs+hwB6gp4ey+eupD5wAaro+h7l+sXfx8H0IBzYDj7vOdQBOAoxz7Z8D/N7L93AIkAWEedj3IfCJl+Pyx1tgDB6OfwZ4F4gHVgNXBPvnyV6Bf9mlLVPmqOpRoDegwPtAiohMEZG6QDcgVlWfUdUMVd3iKnNdvmpeV9XdqpoKfA0k4fwFHwG0E5FwVd2mqsnFiO834FfgCtemAcBJVf3Zh2P3APOAEa5NQ4ADqrrUQ/GeOD2IF1znOgvnkt/1PoQZ46o3y8O+PUCsD3UUJ4aOQC3gJ+BOVf3Kx3ZMOWaJxJRJqrpOVW9V1UZAe5zexas4vYEGInI494Xz13LdfFXsdXt/EohS1c3AA8DTwH4R+UxEGhQzxP9y+pfpDa7PvvoQuMn1/iacS0qeNAB2qGqO27bfcMaMCnMAiBGRMA/76uNcFvRFUWPoCLTG6bXE+9iGKecskZgyT1XXA+NxEsoOYKuq1nR7VVPVS3ys67+q2hsnISnOYLTHooVU9QXQX0QaAVfiPZF4qucrIFFE2gOXAZ94OXY30FhE3P8/bYIzxlOYhTiTE65y3ygiVXEuF871MV6fY3DNBmsCDASeB57wksjMOcYSiSlzRCReRB5y/ZJGRBrj/PX/M87Yx1EReVREKotIqIi092VqsOu+iwEiEoEzSJ2Gc7nLk3044y8eqWoKzhjFv3ES2zpf61HVdGAiTvL5RVW3ezl2Ec54yiMiEu66N+Vy4DNvcbm1cQRnptsbIjLEdXwznAR4AO/JK3+8RYmhI873IhUYhzOmckthsZryzxKJKYuOAT2ARSJyAieBrAYeUtVsnF9kSTgztg4AH+BMby1MBPCC65i9OLOKHvdS9nngL67LZ3/2Uua/OLOzCrqs5a2eD4EOeL+shapmAENxehAHgLeB37l6aIVS1Rdxzu9lnO/pVqAKzmD6CV/iLWIMHXFmyuXG/gLWK6kQRNUetWtMaRORJsB6oJ5rckFptHk7Ti/l/AJ6QcYUmf2lYEwpc403PAh8VlpJBEBVx4lIJnAeYInE+I31SIwpRa7B7n04M5+GqOqOIIdkTIlZIjHGGFMiNthujDGmRCyRGGOMKZEKMdgeExOjzZo1C3YYxhhTbsTExDB9+vTpqjqksLIVIpE0a9aMJUuWBDsMY4wpV0QkxpdydmnLGGNMiVgiMcYYUyKWSIwxxpRIhRgjMcZ4lpmZyc6dO0lPTw92KCaIIiMjadSoEeHhxXsCsyUSL75atotx0xbyRNpLPFv5Ee4Y0pMrOvnyGAhjyo+dO3dSrVo1mjVrhogEOxwTBKrKwYMH2blzJ82bNy9WHXZpy4Ovlu3isUmrGHHiU7rJBkac+C+PTVrFV8t8eQyEMeVHeno60dHRlkQqMBEhOjq6RL1SSyQevDR9A1GZB7g2dA4hoowInUdU5kFemr4h2KEZ43eWRExJfwYskXiw+3Aa94VNIhzncdcRZPCv8LdIPXyIjfuOYeuTGeM/oaGhJCUl5b22bdsW7JCYM2cOl112WaHl+vfvn3eP2iWXXMLhw4cBeP3112nbti033ngjp06d4sILLyQpKYkJEyYENO5gsTESDzrUSGNE+jxyk3SIQO/QNSwJuZuv3+jFCxEXUaVFL3q1iqFXi2iax1S1v+pMhfDVsl28NH0Duw+n0aBmZR4eHFfiscPKlSuzfPnyIh+XlZVFWJh/foVlZ2cTGhpaojqmTp2a9/7tt9/mu+++o3nz5vz8889kZmYW6Rz9eW6lwXokHrxa/3sk36OrMzWUjKr1uTpiEeOyn+DPm25ix5R/cO0/v6LX87P404TlfL54BztSTwYpamMCK3fscNfhNBTYdTgtYGOH6enp3HbbbXTo0IFOnToxe/ZsAMaPH8+IESO4/PLLGTRoEPfccw9TpkwB4Morr+T2228HYOzYsfzlL38B4IorrqBLly4kJCTw3nvv5bURFRXFU089RY8ePVi4cCHTpk0jPj6e3r17M2nSJI9xpaWlcd1115GYmMi1115LWlpa3r5mzZpx4MABRo4cyZYtWxg6dCijR4/mpptuYvny5SQlJZGcnMzSpUvp168fXbp0YfDgwezZswdwejePP/44/fr147XXXiMlJYWrr76abt260a1bNxYsWADA008/ze23307//v1p0aIFr7/+el4M//nPf0hMTKRjx47cfPPNAF7r8afyk/JKUYv0NSBZZ2wLl2xqVa8Gty2AtV/RdNnHjNr+GY9U+oI14d35cENvnljWgUzCaFizMr1aRtOrRTS9WkbToGblIJ2JMb7729drWLvb+3O2lm0/TEZ2zhnb0jKzeWTiSj79xfNzsto1qM5fL08osN20tDSSkpIAaN68OZMnT+att94CYNWqVaxfv55BgwaxceNGABYuXMjKlSupXbs2n332GT/++CNDhw5l165deb+U58+fz3XXXQfAuHHjqF27NmlpaXTr1o2rr76a6OhoTpw4Qfv27XnmmWdIT0+ndevWzJo1i1atWnHttdd6jPWdd96hSpUqrFy5kpUrV9K5c+ezyowZM4Zp06Yxe/ZsYmJi6NGjBy+//DLffPMNmZmZ3Hzzzfzvf/8jNjaWCRMm8MQTTzBu3DgADh8+zNy5cwG44YYb+NOf/kTv3r3Zvn07gwcPZt26dQCsX7+e2bNnc+zYMeLi4rj77rvZuHEjzz33HAsWLCAmJobU1FQA7r//fq/1+IslEk9Gzi94f6ebkE43wYHNhCz/hA4rPuXlnIW8UDOaDXUv5svsC5i0LouJS3cC0Cy6Cr1aRtPTlVjqVIvMqyoQlwqMCYT8SaSw7b7ydGlr/vz5/PGPfwQgPj6epk2b5iWSiy66iNq1awPQp08fXn31VdauXUu7du04dOgQe/bsYeHChXl/qb/++utMnjwZgB07drBp0yaio6MJDQ3l6quvBpxfzM2bN6d169YA3HTTTWf0XnLNmzeP++67D4DExEQSExOLdK4bNmxg9erVXHTRRYBzSa1+/fp5+90T2MyZM1m7dm3e56NHj3Ls2DEALr30UiIiIoiIiKBOnTrs27ePWbNmMXz4cGJinOWxcr9H3uqpVq1akWIviCWSkohpBRf+FS54ApJnEbbsIxI2TCAh52OebNCZPS2uZmZIH+btyOCbFXv49BfnYXit6kTRq0U0oSHw2eIdpGc6/yPmXioALJmYUldYz+H8F2ax63DaWdsb1qzMhLt6+TWWgia0VK1a9XTbDRty6NAhpk2bRt++fUlNTeXzzz8nKiqKatWqMWfOHGbOnMnChQupUqUK/fv3z5vmGhkZeca4iK/jnCUZD1VVEhISWLhwYaHnlpOTw8KFC6lc+ewrGhEREXnvQ0NDycrKQlU9xlZQPf5iYyT+EBoGbQbBtR/BQxtgyAtI1ikazH+C3/00iA+i3mX5zRFM+UMvRl0cT8Oalfny152M/+m3vCSSKy0z26YZmzLp4cFxVA4/c0C6cngoDw+O83tbffv25ZNPPgFg48aNbN++nbg4z+306tWLV199lb59+9KnTx9efvll+vTpA8CRI0eoVasWVapUYf369fz8888e64iPj2fr1q0kJycD8OmnnxYa1+rVq1m5cmWRzisuLo6UlJS8RJKZmcmaNWs8lh00aBBvvvlm3ufCBusHDhzI559/zsGDBwHyLm0VtZ7isETib1WjoefdcPcCuHMOdLoJNs0g9ONhJE7sx0j9gg+vqseKvw7C2981uz381WdMsF3RqSHPX9WBhjUrIzg9keev6hCQ3vM999xDdnY2HTp04Nprr2X8+PFn/BXurk+fPmRlZdGqVSs6d+5MampqXiIZMmQIWVlZJCYm8uSTT9KzZ0+PdURGRvLee+9x6aWX0rt3b5o2beqx3N13383x48dJTEzkxRdfpHv37kU6r0qVKjFx4kQeffRROnbsSFJSEj/99JPHsq+//jpLliwhMTGRdu3aMWbMmALrTkhI4IknnqBfv3507NiRBx98sFj1FEeFeGZ7165dNajPI8lMh/XfwLKPYcscZ1uLfjy9oxOfHutIdU7wZqU3uDfjPlKoScOalVkwakDw4jUVxrp162jbtm2wwzBlgKefBRFZqqpdCzvWxkhKQ3gkdBjuvA5vh+WfwvKPeTpzDg9GVGGP1qK17OaPYZN4Vu8IyKUCY4wJFLu0VdpqNoH+j8J9K+CWrznW4HzahOwiRJTrQ2cTF3WSYUkNgh2lMcb4zBJJsISEQPO+NGzYBAmtBDj3qjxx8iXmbEwJcnDGGOO7gCYSERkiIhtEZLOIjPKwP0JEJrj2LxKRZq7t3UVkueu1QkSudDumpohMFJH1IrJORPw777A0HdsLyz+B7Iy8TT1D17Pj6+eDGJQxxhRNwBKJiIQCbwEXA+2A60WkXb5idwCHVLUV8C9gtGv7aqCrqiYBQ4B3RSR3POc1YJqqxgMdAf/eolma5r4Ieub0X0X43fFxbP3udS8HGWNM2RLIHkl3YLOqblHVDOAzYFi+MsOAD13vJwIDRURU9aSq5q5REgnOwlciUh3oC4wFUNUMVT0cwHMIrJ2/nNEbARCU41Sm+aInYfEHQQrMGGN8F8hE0hDY4fZ5p2ubxzKuxHEEiAYQkR4isgZYBYx07W8BpAD/FpFlIvKBiFTFAxG5U0SWiMiSlJQyOuYwcj48feSs18d9Z/F9dmf49iFLJuact3fvXq677jpatmxJu3btuOSSS/KWQymq8ePHs3v3br/Gl7sYI8B5552Xt/3hhx8mISGBhx9+mJSUFHr06EGnTp348ccf/dp+eRDIROLpfrv8N614LaOqi1Q1AegGPCYikTjTlTsD76hqJ+AEcNbYi+v491S1q6p2jY2NLe45BMWN57Xm0dCHWFGll5NMfnk/2CEZc9qxvfDvi+HYvhJXpapceeWV9O/fn+TkZNauXcs//vEP9u0rXt3FSSRZWVmFF3Jxv3nw3Xff5ddff+Wll17ihx9+ID4+nmXLluXdDFmY7OzsIsVZlgUykewEGrt9bgTk/xfOK+MaA6kBpLoXUNV1OAmjvav8TlVd5No9ESexnFOqRYZz43mtGZE6kuNNL4Kpf7ZkYsqOuS/C9p9h7ujCyxZi9uzZhIeHM3LkyLxtSUlJeb+MX3rpJbp160ZiYiJ//etfAdi2bRtt27bl//7v/0hISGDQoEGkpaUxceJElixZwo033khSUhJpaWk+L9nu7uDBgwwaNIhOnTpx1113nbHuV1RUFABDhw7lxIkT9OjRg9GjR/PII48wderUvHZnzJhBr1696Ny5MyNGjOD48eOA07t55pln6N27N1988QXJyckMGTKELl260KdPH9avXw/Arbfeyn333cd5551HixYtmDhxYl4ML774Ih06dKBjx46MGuX8He2tntISyBsSFwOtRaQ5sAu4DrghX5kpwC3AQmA4MEtV1XXMDlXNEpGmQBywTVUPiMgOEYlT1Q3AQGAt56Dbzm/OBz9u5enKj/JyXLiTTAC6/19wAzPnru9Gwd5VBZfJyoDdS5xJIkv/7ZR3TV/3qF4HuPgFr7tXr15Nly5dPO6bMWMGmzZt4pdffkFVGTp0KPPmzaNJkyZs2rSJTz/9lPfff59rrrmGL7/8kptuuok333yTl19+ma5du5KZmckf//hHn5Zsd/e3v/2N3r1789RTT/Htt996XAV4ypQpREVF5a1bVbduXZYsWcKbb77JgQMHePbZZ5k5cyZVq1Zl9OjRvPLKKzz11FOAsxzL/PnOCuMDBw5kzJgxtG7dmkWLFnHPPfcwa9YsAPbs2cP8+fNZv349Q4cOZfjw4Xz33Xd89dVXLFq0iCpVquStp3XnnXd6rac0BCyRuJLAvcB0IBQYp6prROQZYImqTsEZNP9IRDbj9ESucx3eGxglIplADnCPqh5w7fsj8ImIVAK2ALcF6hyCqXbVStzQownjf9rGfQ+8QxPutmRigu/Idsj9C13VWakhulVAmpoxYwYzZsygU6dOABw/fpxNmzbRpEkTmjdvnvcMky5dunh8PG9Rlmx3N2/evLwHW1166aXUqlWrSHH//PPPrF27lvPPPx+AjIwMevU6fZdCbrvHjx/np59+YsSIEXn7Tp06lff+iiuuICQkhHbt2uVd6ps5cya33XYbVapUAZyl4gurpzQEdIkUVZ0KTM237Sm39+nACA/HfQR85KXO5UCha7+cC+7s24KPFv7GmAU7+MeID+GLW51kogo97gx2eOZcU0DPAXDGRl7ryOmhToX0wzB8HFSrW6wmExISzrhs405Veeyxx7jrrrvO2L5t27azllF3f1Kh+/G+LtmeX0mXir/ooou8riCc225OTg41a9b0uhqv+znmXl7ztFR8YfWUBruzvQyrWz2S4V0bMXHJTvadzIER4yH+MvjuYVh0dnfbmIDycN8TmlOisZIBAwZw6tQp3n//9Bjg4sWLmTt3LoMHD2bcuHF54wu7du1i//79BdZXrVq1vIc/FWXJdnfuS8V/9913HDp0qEjn1LNnTxYsWMDmzZsBOHnypMdZaNWrV6d58+Z88cUXgJMkVqxYUWDdgwYNYty4cZw86TzSOzU1tVj1+JslkjJuZN+WZKvy/rwtEFYJhv/bLZm8G+zwTEXi4b4nsjOc7cUkIkyePJnvv/+eli1bkpCQwNNPP02DBg0YNGgQN9xwA7169aJDhw4MHz48L0l4c+uttzJy5EiSkpLIzs72ecl2d3/961+ZN28enTt3ZsaMGTRp0qRI5xQbG8v48eO5/vrrSUxMpGfPnl4Hvz/55BPGjh1Lx44dSUhI4H//+1+BdQ8ZMoShQ4fStWtXkpKSePnll4tVj7/ZMvLlwJ8mLGfa6r38NGoAtapWcgY8J97mLE0/ZDT0HFl4JcZ4YMvIm1wlWUbeeiTlwN39W5KWmc2/F2x1Nrj3TKY9Cj/7/0E1xhjjK0sk5UCbutUYnFCX8T9t41h6prMxrNLpMZNpj8LP7wQ1RmNMxWWJpJz4wwWtOJqexcc/bz+9MTTcSSZtL4dpoyyZGGOCwhJJOZHYqCZ9Wscwdv4W0jPdllYIDXcuc+Umk4VvBy9IY0yFZImkHLn3glYcOJ7BhMU7ztyRl0yGwvTHYOFbwQnQGFMhWSIpR7o3r03XprV4d24yGVn55vOHhjs3hrUbBtMft2RijCk1lkjKERHhDwNasftIOl8t33V2gdBwuHrs6WTy05ulH6QxReTPZeR9MWfOHC677LKA1H348GHefrviXV62RFLO9G8TS0KD6rwzJ5nsHA/3AOUlkytgxhPw0xulH6Q5p6WcTOHWabdyIO1A4YUL4e9l5IPNEokpF0SEP1zQiq0HTjB11R7PhULD4eoPXMnkL5ZMjF+NWTmGX/f9yjsrSj5LsKBl5FWVhx9+mPbt29OhQwcmTJgAOD2Kfv36cc0119CmTRtGjRrFJ598Qvfu3enQoQPJycnA6bvc+/TpQ5s2bfjmm2/Oav/EiRPcfvvtdOvWjU6dOuXdEf7KK69w++23A7Bq1Srat2+ftyxJrjVr1tC9e3eSkpJITExk06ZNjBo1iuTkZJKSknj44YcB70vhx8fHc8stt5CYmMjw4cPPqr88CeiijSYwhiTUo2VsVd6avZnLEut7XmAut2cCTjJRhfPvK91ATbky+pfRrE/1/hyLpfuWom7Ppvt8w+d8vuFzBKFLXc9LwcfXjufR7o96rbOgZeQnTZrE8uXLWbFiBQcOHKBbt2707dsXgBUrVrBu3Tpq165NixYt+P3vf88vv/zCa6+9xhtvvMGrr74KOL+w586dS3JyMhdccEHe+le5nnvuOQYMGMC4ceM4fPgw3bt358ILL+SBBx6gf//+TJ48meeee4533303b8XdXGPGjOH+++/nxhtvJCMjg+zsbF544QVWr16dt4BiQUvhb9iwgbFjx3L++edz++238/bbb/PnP//Z6/eqLLMeSTkUEiLc3b8V6/ceY9b6AhaxCw1zkknClfD9k7DgNb8+3c5ULB1iOlA7ojbierCpINSOrE1iTGJA2ps/fz7XX389oaGh1K1bl379+rF48WIAunXrRv369YmIiKBly5YMGjTIibFDhzOWlL/mmmsICQmhdevWtGjR4qw1r2bMmMELL7xAUlIS/fv3Jz09ne3btxMSEsL48eO5+eab6devX96S8O569erFP/7xD0aPHs1vv/1G5cqVzyrjvhR+586dWb9+PZs2bQKgcePGefXedNNNec8oKY+sR1JODUtqwL++38ibszczIL6O92WvQ8Pgqg8Age+fgnVfw66lzoqtl71SqjGbsq2gnkOuZxY+w8SNE6kUWonM7EwubHohT/Z8sthtFraMvDfuS6yHhITkfQ4JCTnj0bn5/7/I/1lV+fLLL4mLizurjU2bNhEVFeX10b033HADPXr04Ntvv2Xw4MF88MEHtGjR4qz6vS2FX1hs5Yn1SMqp8NAQRvZrwbLth1m45WDBhUPD4Kr3oc3FsHOxs/T38k+sV2KKLDU9lWviruG/l/yXa+Ku4WBaIT97hShoGfm+ffsyYcIEsrOzSUlJYd68eXTv3r1I9X/xxRfk5OSQnJzMli1bzkoYgwcP5o033shLWsuWLQPgyJEj3H///cybN4+DBw96THZbtmyhRYsW3HfffQwdOpSVK1eesYx9bv3elsLfvn173hL3n376Kb179y7SuZUl1iMpx0Z0bcxrP2zm7dnJnNcypuDCoWFQrd7pz7nPkbBeiSmCVy94Ne/9X3r+pcT15S4j/8ADD/DCCy8QGRlJs2bNePXVV+nbty8LFy6kY8eOiAgvvvgi9erVK9LzyOPi4ujXrx/79u1jzJgxREZGnrH/ySef5IEHHiAxMRFVpVmzZnzzzTf86U9/4p577qFNmzaMHTuWCy64gL59+1KnTp28YydMmMDHH39MeHg49erV46mnnqJ27dqcf/75tG/fnosvvpiXXnqJdevW5T0hMSoqio8//pjQ0FDatm3Lhx9+yF133UXr1q25++67S/z9DBZbRr6ce3duMs9/t56v/nA+SY1rei+Y+3S7rPTT28Ii4f6VxX66nSn/zuVl5G+99VYuu+wyhg8fHuxQzrJt2zYuu+wyVq9eHexQ8tgy8hXYjT2bUqNyOG/O2lxwwQA83c4YY8AubZV7URFh3HZ+M16duYn1e48SX6+654IBeLqdMWXZ+PHjgx2CV82aNStTvZGSsh7JOeDW85pRtVIob89O9l5o5Hx4+gjcs8j5fNX7zueR5XfKoTGmbLBEcg6oWaUSN/Vsyjcrd7PtwImCC9duASFhkOL7gKU5t1WEcVJTsJL+DFgiOUfc0bs5YaEhjJlbQK8EnCcr1m4JKRtKJzBTpkVGRnLw4EFLJhWYqnLw4MGzZrQVRUDHSERkCPAaEAp8oKov5NsfAfwH6AIcBK5V1W0i0h14L7cY8LSqTnY7LhRYAuxS1cAs41nO1KkeybVdG/PZ4u3cf2Fr6tc4+y7bPLFxsH9t6QVnyqxGjRqxc+dOUlJSgh2KCaLIyEgaNWpU7OMDlkhcv+zfAi4CdgKLRWSKqrr/BrsDOKSqrUTkOmA0cC2wGuiqqlkiUh9YISJfq2ruLav3A+sALyPLFdNd/Vrw6S/beW/eFv56eYL3grHxsP4byDoFYRHey5lzXnh4OM2bNw92GKacC+Slre7AZlXdoqoZwGfAsHxlhgEfut5PBAaKiKjqSbekEQmnV4oTkUbApcAHAYy9XGpUqwrDkhry6S/bOXD8lPeCsXHO1N+DhUwZNsYYHwQykTQE3J8Ju9O1zWMZV+I4AkQDiEgPEVkDrAJGuiWWV4FHgHw3RZxJRO4UkSUisqQiddvv7t+SU1k5/HvBVu+FYuOdrzbgbozxg0AmEk8rkOUf0fNaRlUXqWoC0A14TEQiReQyYL+qLi2scVV9T1W7qmrX2NjYosZebrWqE8XF7evxn59+40hapudC0a1AQmzA3RjjF4FMJDuBxm6fGwH5l9HMKyMiYUANINW9gKquA04A7YHzgaEisg3nUtkAEfk4EMGXZ/f0b8WxU1l8/PNvnguER0Kt5tYjMcb4RSATyWKgtYg0F5FKwHXAlHxlpgC3uN4PB2apqrqOCQMQkaZAHLBNVR9T1Uaq2sxV3yxVvSmA51AutW9Yg/5xsYydv5WTGVmeC8XGW4/EGOMXAUskrjGNe4HpODOsPlfVNSLyjIgMdRUbC0Q/Pr3FAAAgAElEQVSLyGbgQWCUa3tvnJlay4HJwD2qWvIHRFcg917QitQTGXz2yw7PBWLjnMH2bC+Xv4wxxkcBvY9EVacCU/Nte8rtfTowwsNxHwEfFVL3HGCOP+I8F3VtVpsezWvz3rwt3NizCRFhoWcWiI2HnCxI3eIkFWOMKSa7s/0c9ocLWrH3aDqTft119s7c5GHjJMaYErJEcg7r0zqGxEY1GDM3mazsfLOlY9oAYuMkxpgSs0RyDhMR7unfit8OnuTbVXvO3FmpCtRsYj0SY0yJWSI5xw1qV5c2daN4e3YyOTn5buOxmVvGGD+wRHKOCwlxeiUb9h1j5rp9Z+6MjYMDmyDbyxRhY4zxgSWSCuCyxPo0qV2Ft2ZvPnO58Nh4yD4Fh73cuGiMMT6wR+1WAGGhIYzs15LHJ6+i67MzST2RQYOalXmue236A+xfB9EtgxylMaa8sh5JBVEp1FnW7OCJDBTYdTiNh2anOzttwN0YUwKWSCqIf83cdNa2g5mV2EuMDbgbY0rEEkkFsftwmsftG7IbWI/EGFMilkgqiAY1PT96d3elZnBgI+Rkl25AxphzhiWSCuLhwXFUDj9zva3K4aHEd+gGWelweHuQIjPGlHeWSCqIKzo15PmrOlC/RiQAURFhPH9VBzp16ekUsHESY0wxWSKpQK7o1JCFjw2kT+sY6laP4IpODV1rbmHjJMaYYrNEUgENiK9DcsoJth04AZVrQrX61iMxxhSbJZIKaGB8XQB+WL/f2RAbZz0SY0yxWSKpgJpEV6F1nSh+yF17K3fxRtWCDzTGGA8skVRQA9vW5ZetqRxNz3R6JJkn4MjOYIdljCmHLJFUUAPb1iErR5m3McXpkYCNkxhjisUSSQXVuUktalYJZ9a6/W6JxMZJjDFFZ4mkggoNES6Iq8PsDfvJjqwFVWMtkRhjisUSSQU2sG0dDp3MZNn2Q/a0RGNMsQU0kYjIEBHZICKbRWSUh/0RIjLBtX+RiDRzbe8uIstdrxUicqVre2MRmS0i60RkjYjcH8j4z3V9WscSFiLMXLffNQXYZm4ZY4rOp0QiIl+KyKUi4nPiEZFQ4C3gYqAdcL2ItMtX7A7gkKq2Av4FjHZtXw10VdUkYAjwroiEAVnAQ6raFugJ/MFDncZHNSqH061ZbWat3+f0SE4dgWN7gx2WMaac8TUxvAPcAGwSkRdEJN6HY7oDm1V1i6pmAJ8Bw/KVGQZ86Ho/ERgoIqKqJ1U190HikYACqOoeVf3V9f4YsA5o6OM5GA8Gtq3Dxn3H2R/RzNlg4yTGmCLyKZGo6kxVvRHoDGwDvheRn0TkNhEJ93JYQ2CH2+ednP1LP6+MK3EcAaIBRKSHiKwBVgEj3RILrv3NgE7AIl/OwXg2sK1zl/usg7WcDTZOYowpoqJcqooGbgV+DywDXsNJLN97O8TDtvwX4L2WUdVFqpoAdAMeE5FIt1iigC+BB1T1qJd47xSRJSKyJCUlxet5VXTNY6rSIrYq327Jgsq1rEdijCkyX8dIJgE/AlWAy1V1qKpOUNU/AlFeDtsJNHb73AjY7a2MawykBpDqXkBV1wEngPaucuE4SeQTVZ3kLWZVfU9Vu6pq19jYWF9Os8IaGF+HRVsPkR0dZz0SY0yR+doj+UBV26nq86q6B5wZVwCq2tXLMYuB1iLSXEQqAdcBU/KVmQLc4no/HJilquo6JszVTlMgDtgmIgKMBdap6is+xm4KMbBtXTKyc9gd3gRS1tnMLWNMkfiaSJ71sG1hQQe4xjTuBabjDIp/rqprROQZERnqKjYWiBaRzcCDQO4U4d7AChFZDkwG7lHVA8D5wM3AALfpwZf4eA7Giy5Na1E9MoylaXUh7RCcOBDskIwx5UhYQTtFpB7OgHhlEenE6TGN6jiXuQqkqlOBqfm2PeX2Ph0Y4eG4j4CPPGyfj+dxFVMC4aEh9I+rw4xNNbkCnF5JlF0ONMb4psBEAgzGGWBvBLhfSjoGPB6gmEwQDGxbh3+sqOtMtk7ZAM37BjskY0w5UWAiUdUPgQ9F5GpV/bKUYjJB0K9NLA+G1OZUaFUibOaWMaYICru0dZOqfgw0E5EH8++3Ae9zR80qlejStDbJ+xvTzmZuGWOKoLBLW1VdX71N8TXnkAvb1mHljnrE7V9NaLCDMcaUG4Vd2nrX9fZtVbW7+s5xA+Lr8un0hoSenAMnDkLV6GCHZIwpB3yd/vuTiMwQkTtEpFZAIzJB0zK2KseiWjofDtjlLWOMb3xda6s18BcgAVgqIt+IyE0BjcyUOhGhQZtOAGTsWRvkaIwx5YXPa22p6i+q+iDOqr6pnF6115xDuid24LhGsmfz8mCHYowpJ3xda6u6iNwiIt8BPwF7cBKKOcd0bR7NVhqSsXddsEMxxpQThc3ayrUC+Ap4RlULXBrFlG+VwkI4Ub0VdY8tJidHCQmxhQSMMQXz9dJWC1X9kyWRiqFqowTqkMrarduDHYoxphwo7IbEV1X1AWCKiJy1JKyqDvVwmCnnmrXtAutg1fLFtG/ZNNjhGGPKuMIubeUunPhyoAMxZUe1RgkA7N+yAmd1f2OM8a7AS1uqutT1NklV57q/gKTAh2eComYTMkMiqHJkM3uPpAc7GmNMGefrGMktHrbd6sc4TFkSEkp27da0ll3MWr8/2NEYY8q4AhOJiFwvIl8DzUVkittrNnCwdEI0wRBRvx3xYbv5Yd2+YIdijCnjChsjyb1nJAb4p9v2Y8DKQAVlgk9i46inn7Ns8w7SMjpTuZIt42iM8aywRRt/A34DepVOOKbMiI0HoFH2Tn5KPsDAtnWDHJAxpqwq7NLWfNfXYyJy1O11TESOlk6IJihciaR9+B5+sHESY0wBCpu11dv1tZqqVnd7VVPV6qUTogmKWs0gtBL9aqcya91+VM+6jcgYYwDf19pqKSIRrvf9ReQ+EakZ2NBMUIWGQXRrOlTaw96j6azZbR1QY4xnvk7//RLIFpFWwFigOfDfgEVlyobYOOqe2oYI/LDOLm8ZYzzzNZHkqGoWcCXwqqr+CagfuLBMmRAbT+iR7XRvFMms9TYN2Bjjma+JJFNErse5MfEb17bwwg4SkSEiskFENovIKA/7I0Rkgmv/IhFp5treXUSWu14rRORKX+s0fhQbByhXNj7Jip1H2H/M7nI3xpzN10RyG84U4OdUdauINAc+LugAEQkF3gIuBtoB14tIu3zF7gAOqWor4F/AaNf21UBXVU0ChgDvikiYj3Uaf6nTFoA+NZx7T2fb7C1jjAe+Pmp3rarep6qfuj5vVdUXCjmsO7BZVbeoagbwGTAsX5lhnH7S4kRgoIiIqp50XUoDiARypwz5Uqfxl9otICSMBpm/0aBGpI2TGGM88nXW1vki8r2IbBSRLSKyVUS2FHJYQ2CH2+edrm0ey7gSxxEg2tVmDxFZA6wCRrr2+1Jnbsx3isgSEVmSkpLiy2ma/ELDIboVkrKBAW3r8OOmA6RnZgc7KmNMGePrpa2xwCtAb6Ab0NX1tSCeHq2X/2YEr2VUdZGqJrjaeUxEIn2sE9fx76lqV1XtGhsbW0ioxqvYOEhZz8C2dUnLzObnLbbEmjHmTL4mkiOq+p2q7lfVg7mvQo7ZCTR2+9wI2O2tjIiEATWAVPcCqroOOAG097FO40+x8XBoK72aVKVyeKhd3jLGnMXXRDJbRF4SkV4i0jn3Vcgxi4HWItJcRCoB1wFT8pWZwukl6ocDs1RVXceEAYhIUyAO2OZjncafYuNAc4g8spXerWOYtd7ucjfGnKmw1X9z9XB97eq2TYEB3g5Q1SwRuReYDoQC41R1jYg8AyxR1Sk4l8w+EpHNOD2R61yH9wZGiUgmkAPco6oHADzV6eM5mOJwrblFynoGxnfn+7X7WL/3GG3r2wo5xhiHT4lEVS8oTuWqOhWYmm/bU27v04ERHo77iNOP+S20ThNA0a1AQiBlAwO6XgbArPX7LZEYY/L4OmurroiMFZHvXJ/bicgdgQ3NlAlhEc404JT11KkeSWKjGsy0h10ZY9z4OkYyHudyUgPX543AA4EIyJRBsfGQsgGAgfF1Wb7jMAeOnwpyUMaYssLXRBKjqp/jjFfk3vNhNxRUFLFxkJoMWRkMbFsHVbvL3Rhzmq+J5ISIROO6Z0NEeuLcPGgqgth4yMmC1C0kNKhO3eoRzLJEYoxx8TWRPIgzzbaliCwA/gP8MWBRmbIlNs75mrIeEWFAfF3mbUzhVJZ1So0xhT9qt5uI1FPVX4F+wOPAKWAGzs2BpiKIbg1I3jjJhW3rcCIjm1+2phZ8nDGmQiisR/IukOF6fx7wBM7qu4eA9wIYlylLKlWBWk0hZT0A57WMISIsxO5yN8YAhSeSUFXN/bPzWuA9Vf1SVZ8EWgU2NFOmuM3cqlwplPNbxfDD+n12l7sxpvBEkrtUCTAQmOW2z9e74s25IDYODm6CbGd1/4Ft67AjNY3N+48HOTBjTLAVlkg+BeaKyP+ANOBHANez223WVkUSGw/ZGXBoGwAD4usAMNMubxlT4RWYSFT1OeAhnBsSe+vp6xgh2KytisVt5hZA/RqVSWhQ3Z7lbowpfPqvqv6sqpNV9YTbto2umVymoohp43x1JRKAgfF1WPrbIQ6dyPBykDGmIvD1PhJT0UVUgxqN8wbcAQa2rUuOwpyNdnnLmIrMEonxnetpibk6NKxBbLUIGycxpoKzRGJ8FxsPBzZCjnNHe0iIMCCuDvM2pJCZnRPk4IwxwWKJxPguNh6y0uHw9rxNA9rW4dipLBbbXe7GVFiWSIzv3J6WmKt3qxgqhYXwgy3iaEyFZYnE+C727JlbVSPC6NUimh/W2V3uxlRUlkiM7yJrQLUGZ8zcAmcRx20HT7LlwAkvBxpjzmWWSEzR5Ju5BXCB6y73H+wRvMZUSJZITNHExkPKRsg5PUurUa0qxNerZqsBG1NBWSIxRRMbB5kn4OiZj6MZ2LYOS347xJGTmUEKzBgTLJZITNHkzdw6c5xkYNu6ZOeo3eVuTAUU0EQiIkNEZIOIbBaRUR72R4jIBNf+RSLSzLX9IhFZKiKrXF8HuB1zvWv7ShGZJiIxgTwHk0++xRtzdWxUk+iqlexZ7sZUQAFLJCISivM0xYuBdsD1ItIuX7E7gEOq2gr4FzDatf0AcLmqdgBuAT5y1RkGvAZcoKqJwErg3kCdg/GgSm2oWuesRBIaIvSPq8OcDSlk2V3uxlQogeyRdAc2q+oWVc0APgOG5SszDPjQ9X4iMFBERFWXqepu1/Y1QKSIRADielUVEQGqA7sxpSs27qxLW+BMAz6SlsnS3w4FIShjTLAEMpE0BHa4fd7p2uaxjKpm4TwsKzpfmauBZap6SlUzgbuBVTgJpB0w1lPjInKniCwRkSUpKSklPRfjLvexu/luQOzdOobwULG73I2pYAKZSMTDtvy3PhdYRkQScC533eX6HI6TSDoBDXAubT3mqXFVfU9Vu6pq19jY2KJHb7yLjYNTR+HYnjM2V4sMp6frLndjTMURyESyE2js9rkRZ1+GyivjGv+oAaS6PjcCJgO/U9VkV/kkAFVNdj2t8XPgvECdgPHCw5pbuQbE1yE55QTb7C53YyqMQCaSxUBrEWkuIpWA64Ap+cpMwRlMBxgOzFJVFZGawLfAY6q6wK38LqCdiOR2MS4C1gXsDIxnXqYAAwyMrwtgl7eMqUAClkhcYx73AtNxftl/rqprROQZERnqKjYWiBaRzcCDQO4U4XuBVsCTIrLc9arjGoD/GzBPRFbi9FD+EahzMF5UjYHKtT32SJpEV6F1nSi7vGVMBRIWyMpVdSowNd+2p9zepwMjPBz3LPCslzrHAGP8G6kpEpHTA+4eDGxblw9+3MLR9EyqR4aXcnDGmNJmd7ab4omNg/3rzpq5BVApTMjKURKfnsH5L8ziq2W7ghCgMaa0WCIxxRMbD+mH4cSZU6u/WraL9+dtyfu863Aaj01aZcnEmHOYJRJTPF6WSnlp+gbSMs+8sz0tM5uXpnu+DGaMKf8skZji8TJza/fhNI/Fdx1O46fkA+Tk2FMUjTnXBHSw3ZzDqtWDiBpn9Uga1KzMLg/JRIAb3l9EveqRDEtqwLCkhrStXw1npRtjTHlmicQUj4jHNbceHhzHY5NWkZaZnbetcngofxuaQGSlUP63bBdj52/l3XlbaFM3imFJDRmW1IBGtaqU9hkYY/zEEokpvjrxsOG7MzZd0clZTu2l6RvYfTiNBjUr8/DguLztQzs2IPVEBt+u3M1Xy3fz0vQNvDR9A92a1eKKTg25tEN9alapVOqnYowpPlEP0zfPNV27dtUlS5YEO4xzz8K3YPrj8PAWqJp/rU3f7Eg9yf+W7+Kr5bvZvP844aFCvzZ1uKJTAy5sW5fI8FA/B22M8ZWILFXVroWVsx6JKT73mVtVzy9WFY1rV+HeAa35wwWtWLP7KF8t28WUFbuZuW4fURFhDGlfjyuSGtKrZTShITaeYkxZZInEFJ/74o3NipdIcokI7RvWoH3DGjx2SVt+3nKQr5btYtrqvUxcupM61SK4vGMDrkhqSPuG1fMG6b9atsvrZTRjTOmw6b8FSDmZwq3TbuVA2oFzoh2/t1W9IVSK8rhUSknaCQ0Rzm8Vw0sjOrL4Lxfy9o2dSWpck/8s3Mblb85n4CtzeeOHTYz9cQuPTVrFqeMb6dzscU4d2xSwmx+/WraLS0b/m6vHdOTi0eMDeoNlabVl51Sy2EqrrdI8p+KyHkkB3ln+Dr/u+5W3lr3FYz08PvbEL95e/naptBOQtmLbQMo6yM4ISDshIXBhu2gubBfNkZOZTFuzl/+t2MU/Z67NK9O5/kesisymQ93/sGDPI/z925XUqxlGiDhJKUQEcXsfgpy5L8T1WQQRcW2HENf+aav28PdvVtMl9gNWRmaTmPE+T06qQUZGOy7tUM9tlRhFXY/TUcXD8jGat0nRvP3OF+f99NX7GD19Pe1j32dFZDYdM97jb5Oqkn4inkHtnZWVxe2/ufJmUbtNp5azHvcjebu/W72X56auo6OrHeecqnEqvR0Xd6jn+z+QD75btZdnp64l6Yy2qpOR3o6LE/3X1ncr9/L3qWvpeMa/U8HtFHTMJYn188rln6U+deUenvn27OOyTyVwWVJ91zFnHvTNij389es1+Y6pQU5mApd1rI8nfjnm1Ps8Nsk5/0D11m2w3YMuH3chI98vRmOMKdc0jFW3LivSITbYXgLTrprGy0teZsZvM8jKySJMwmhZsyV9GvWhanhVv7VzIvMEP+78keTDyWRp4NoJaFvbFsDmmdDvEQivXKrn9M7cNdSrNJPdlZQsEcJUaZyRQ6uTUTStHoUCiqAIIKiI6zNAiNt+13sRVCXvmNz9Gw8d53DU3jPaqZ8hVD/WgDbR1Z1g8v749DYhwMt2OfPD+gOHOR61iz352oo61oj42Bqucl7++CvC34QbDh7hmId2qh9rSJuY6r5X5IONB49y1Gtb1fzYzrEit+PTMR6+rxtTvR/XunY118/OmZIPeT+mec2qeeVP91ph+9ETXo+pUzmSHFcvN0dzv8LRrIwzjonIyaHF8Ros2XufH77Lnlki8SC2SixVw6uSnZNNpdBKZGZn0rFOR+7vfL/f2zpy6ggbD20MeDsBayuiMSydBDHdoUmPwLXjQfvtHzJjy2QmVoqiUk4OmSJ0Tk+nR4//cnGvJL+1M/T5iXQNf4TJlarktdM9/QRLM+7kmVuv9ls7uW1189jW73nmd/5ry3s7/8czt5TWOfm3reK0U9zYCjru77cVva3n7yj6Ma8/5NsxGSK0072kRQXukQ422O5Fanoq18Rdw38v+S/XxF3DwbSD5bqdgLXlYfHGUjmntMP0XPEYqaGhXHPsOP/ds49rjh3nUFgIFx/8j1+berX+9xwKDTmjndSwEF6rP8Ov7ZRmW3ZOJYuttNoqzXMqCRsjMSWTkwP/aABdb4chpfSwSlWYeDusmeR5f70OMHK+/9ob0xv2rgp8O6XZlp1TyWIrrbZK85w88HWMxBKJKbl3+0KVGLjZyy92f1v+KXw1Egb8Bfo+XDptGlMB+ZpI7NKWKbkCHrvrd6lbYOqfocl50PvB0mnTGFMgSySm5GLj4OhOSD8a2Hays2DSnSChcNV7EGLrcBlTFlgiMSWXu1TKgU2BbWfei7BzMVz+L6jZOLBtGWN8ZonElJz7mluBsv1nmPcSdLwe2vt3eqoxpmQCmkhEZIiIbBCRzSIyysP+CBGZ4Nq/SESaubZfJCJLRWSV6+sAt2Mqich7IrJRRNaLiP1WCbaaTSE0InCJJP0IfPl/ULMJXPxiYNowxhRbwG5IFJFQ4C3gImAnsFhEpqjqWrdidwCHVLWViFwHjAauBQ4Al6vqbhFpD0wHcheJeQLYr6ptRCQEqB2oczA+Cg2DmNaBG3D/9s9wdBfcPg0i/XvXtTGm5ALZI+kObFbVLaqaAXwGDMtXZhjwoev9RGCgiIiqLlPV3a7ta4BIEYlwfb4deB5AVXNUNfBL5prCxcYFpkey8nNY9Tn0exQad/d//caYEgtkImkI7HD7vJPTvYqzyqhqFnAEyP+ovauBZap6SkRqurb9XUR+FZEvRKSu/0M3RRYbD4e3Q8YJ/9V5aBt8+xA07gl9HvJfvcYYvwpkIvG0Sl3+ux8LLCMiCTiXu+5ybQoDGgELVLUzsBB42WPjIneKyBIRWZKSklLU2E1RxcYB6r+ZW9lZMMn1z37Ve87lM2NMmRTIRLITcJ+j2QjY7a2MiIQBNYBU1+dGwGTgd6qa7Cp/EDjp2g7wBdDZU+Oq+p6qdlXVrrGxsSU/G1OwvJlbfhon+fGfsONnuPSfUKupf+o0xgREIBPJYqC1iDQXkUrAdcCUfGWmALe43g8HZqmqui5hfQs8pqoLcgurs57L10B/16aBgPvgvQmW2i0gJNw/4yQ7foG5o6HDNZB4TcnrM8YEVMASiWvM416cGVfrgM9VdY2IPCMiQ13FxgLRIrIZeBDInSJ8L9AKeFJElrtedVz7HgWeFpGVwM2AXTwvC0LDIbpVyXsk6Ufhy99DjYZwqcerlsaYMiagF55VdSowNd+2p9zepwMjPBz3LPCslzp/A/r6N1LjF7FxnlcdLYrvHoEjO+C27yCyRuHljTFBZ3e2G/+JjYdDWyEzvXjHr5oIKz51VvRt0tO/sRljAsYSifGf2DjQHDhYjJlbh7fDNw9Co27Q9xH/x2aMCRhLJMZ/ijtzKyfbmeqr2XDV+zbV15hyxv6PNf4T3dJZ4r2oM7fmvwLbf4IrxkDt5oGJzRgTMNYjMf4TFuFMAy5KItm5FGY/76zo2/G6wMVmjAkYSyTGv2LjfL+0deoYfHkHVG8Al74C4mmhA2NMWWeJxPhXbDwcTIasjMLLfjcKDv/mLIFSuWbh5Y0xZZIlEuNfsfHOoHlqcsHl1kyG5R87z11vel7pxGaMCQhLJMa/YuOcrwWNkxzZCV/fDw27QP+znndmjClnLJEY/4ppDYj3cZLcqb7ZWa6pvuGlGp4xxv9s+q/xr/DKUKuZ9x7Jgtfgt/kw7C1nurAxptyzHonxv9h4zz2SXb/C7Oeg3RWQdGPpx2WMCQhLJMb/YuOcB1xlZ53elnHCWdU3qi5c/qpN9TXmHGKJxPhfbDzkZDoLOOaaNgpSt8CV70LlWsGLzRjjd5ZIjP/ln7m1dgr8+h/o/QA07xO8uIwxAWGJxPhfTBvna8p6OLobvr4P6idB/8eDG5cxJiBs1pbxv4goqNYQfh4Dm2ZC1im4eiyEVQp2ZMaYALBEYgJDBE4ecF6Xvw4xrYIdkTEmQOzSlvG/Y3vh2B7nvYRAm8HBjccYE1CWSIz/zX3x9PTekFDnszHmnGWJxPjXsb2w/BPIcd1Dkp3pfD62L7hxGWMCxhKJ8a+5LzrPbXenOTB3dHDiMcYEnCUS4187f4HsfM8iyc5wthtjzkkBnbUlIkOA14BQ4ANVfSHf/gjgP0AX4CBwrapuE5GLgBeASkAG8LCqzsp37BSghaq2D+Q5mCIaOT/YERhjSlnAeiQiEgq8BVwMtAOuF5F2+YrdARxS1VbAv4Dc6x8HgMtVtQNwC/BRvrqvAo4HKnZjjDG+C+Slre7AZlXdoqoZwGfAsHxlhgEfut5PBAaKiKjqMlXd7dq+Boh09V4QkSjgQeDZAMZujDHGR4FMJA2BHW6fd7q2eSyjqlnAESA6X5mrgWWqesr1+e/AP4GTBTUuIneKyBIRWZKSklK8MzDGGFOoQCYST+uEa1HKiEgCzuWuu1yfk4BWqjq5sMZV9T1V7aqqXWNjY32P2hhjTJEEMpHsBBq7fW4E7PZWRkTCgBpAqutzI2Ay8DtVTXaV7wV0EZFtwHygjYjMCVD8xhhjfCCq+TsJfqrYSQwbgYHALmAxcIOqrnEr8wegg6qOFJHrgKtU9RoRqQnMBZ5R1S+91N8M+MaXWVsikgL8VsJTCrQYnEkG5xI7p/LBzql8KO1zOgCgqkMKKxiw6b+qmiUi9wLTcab/jlPVNSLyDLBEVacAY4GPRGQzTk/kOtfh9wKtgCdF5EnXtkGqur+YsZT5a1siskRVuwY7Dn+ycyof7JzKh7J8TgHrkZiiKcs/JMVl51Q+2DmVD2X5nOzOdmOMMSViiaTseC/YAQSAnVP5YOdUPpTZc7JLW8YYY0rEeiTGGGNKxBKJMcaYErFEYowxpkQskZRhItJCRMaKyES3bVVF5MP/b+9+Xqwq4ziOvz9CbdwEQqsIWpirWfQnBKItJH9CQVCaBNPCiBBSElyIiVslEgb0ripE/FVqIBi4cRH9WAyK0MLFbAqEQLRF6LfFOTPehpnxnnmec55zr58XzObMDz4f7j3Pd84MPI+kGUnvlcyXQtKrki5LOi3pQOk8OUhaI+mopJOSPiidJ5f6PfeLpC2ls9WpRYgAAAMESURBVOQgaVt9/1yStKl0ntXq01rgQdKSeoH8S9LsoutvSbor6Y9nLaD1zsl7F13eAZyLiI+AtzPHHkmObsDrwJWI+JDqmIGiMnXaSrUR6b9U2/8UlakTwOfA2XZSNpPpvrpY3z+7gXdajNtYw37F14J5HiTtGQD/21pguTNaJE1J+mHRx8vL/NxXeLqr8uOWsj/LgPRuvwHvSroB/NRx/qUMSO+0AbgVEZ8BH3ecfykDEjtJ2gjcBv7sOvwyBuS7rw7V39cnA0bsRz/WAqDlExKfZxFxU9V+YMMWzmgBkPQdsDUijgGj/tlgjuoN9DuFfhHI0U3SfuBw/bPOAWfaTb2yTJ3mqE70hMI3NmTr9CawlmoB+0fS1Yh40mrwFWTqJKoTWK9FxK/tJm6mST96sBbM8xNJt0Y5o2WBpHWSTgFvSDpYXz4P7JT0NfB9a0mba9QN+BH4pO53r8VcKZp2Og9slnQSuNlmsASNOkXEFxHxKfANMFNyiKyg6eu0D9gI7JI03WawTJbr15u1wE8k3RrljJann4i4D0wvuvYQ2JM5Vw5Nu80Cu9qLk0XTTo+ojo/us0adFr4gYpA/SjZNX6cTwIn24mS3ZL8+rQV+IunWKGe0jKtJ7OZO42ESOw3rfT8Pkm79DKyX9JqkF6m2zb9cOFMuk9jNncbDJHYa1vt+HiQtkfQtcAvYIGlO0t76XPr5M1ruAGeHD/oaF5PYzZ3GwyR2Gjau/bxpo5mZJfETiZmZJfEgMTOzJB4kZmaWxIPEzMySeJCYmVkSDxIzM0viQWJmZkk8SMzMLIkHiVkBks5I2iLpJUnXJG0vnclstTxIzMqYAv4GLgFHIuJC4Txmq+YtUsw6JmkN8AC4D3wVEccLRzJL4icSs+6tp9oGfDcwLemFsnHM0niQmHVvCrgeETeAWeD9wnnMkniQmHVvimqAAHwJHJTk00ptbPl/JGZmlsRPJGZmlsSDxMzMkniQmJlZEg8SMzNL4kFiZmZJPEjMzCyJB4mZmSXxIDEzsyT/Ad8I20aW/0sqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEcCAYAAADtODJSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd8U9X7wPHPSbooLZS9oWWV1VL23nuIiIAsFXEAypKfKCii8hUFRGQpiICIDNmobNmryJQ9C6WUWUaheyTn98dNS1rSNh3poOf9euVFc3Nz7pPQ5sk959znCCkliqIoipJWuqwOQFEURcnZVCJRFEVR0kUlEkVRFCVdVCJRFEVR0kUlEkVRFCVdVCJRFEVR0kUlEkVRFCVdVCJRFEVR0kUlEiVZQgh/IUTbrI4jswkhPIUQJ4UQIUKIEVkcyzkhRMu0Pp7WfdMjO71/iu2pRKIAIITYI4R4LIRwzOpYsomPgT1SSlcp5aysDERKWV1KuSfufuLknvhxa9uy8ZeEFN8/IcRAIcQZIUS4EOKuEOInIUT+pBrMrV9qcgKVSBSEEO5AM0AC3bI0mESEEHbWbEttG1YoB5xLw/MUTbLvnxDi/4ApwBggP9AQcAe2CyHsMyNAJQNJKdUtl9+ACcBBYDqwMdFj/sA44DzwGPgVcDI99glwCwgBLgFtrDxeSWAtEARcB0ZYOOYnwGkgCrCzsM0L2AMEo31gdUupDQtxVLXUBrALMACRQChQ2cJzLb52K1/bR6a4ngArU3o/Tc9pa/r5d8AIRJhi+zjR42OBNYmOOROYZd5WEu2MAdYmeu5sYEYS/49pev+AfKbtvRNtdwHuA29aONZz8SYXQxLx9gQOm34uCvgCo7P67+9FuGV5AOqW9TfgKvA+UAeIAYqZPeYPnAXKAAXREs7XgCdwEyhp2s8dqGDFsXTAcbTk5QCUB64BHRId8z/TMfNY2OZqivlTUxutTR++nsm1kSgO++TaMH04vZPEa7D42lPx2o6gJZyCwAVgSHLvJ2aJIqX7aGcC4UA+0309cAdoaGHfxO2UAMIAN9N9O7QP9joZ/P51BGKxnNx/A5Yl8bzE8SYbg4XnTwR+Bqqg/U53z+q/vRflprq2cjkhRFO0D59VUsrjgB/QL9Fuc6SUN6WUj4BJQF+0b5yOQDUhhL2U0l9K6WfFIesBRaSUE6WU0VLKa8AvQJ9E+80yHTMi8TbAB+3b62RTG7uAjaa4UmojTkMr27Akqdeemtd22/R+/m16PWl9PxOQUt4ATgDdTZtaA+FSysNWPPcOsA/oZdrUEXhg+r1ILD3vX2FTu7EWHrsDFLGijbTEUBMoABwC3pNSbrDyOEoKVCJR3gS2SykfmO4vN20zd9Ps5xto35qvAqOAL4H7Qog/hBAlrTheOaCkECI47ob2jbJYMsdMvK0kcFNKaUwUVykr2ohjbRvPSea1W/va7pr9HA64pOP9tGQ5zz5M+5nuW+s3YIDp5wFoXUqWpPn9Ax4AhZMYuyqB1i1ojdTGUBOohHbWUsXKYyhWUIkkFxNC5AF6Ay1Ms2buAh8CNYUQNc12LWP2c1ngNoCUcrmUMu6MRqINnqbkJnBdSulmdnOVUnZOtJ+lhXLitt0GygghzH9/y6KNL6TURhxr27Aoiddu7WtLTZsWd02hqdVASyFEaeAVkk4kltrZAHgLIWoAXYFlSTw3Pe+fL9q4VQ/zjUKIvEAnYK+V8Vodg2k2WFmgDfAt8FkaJ2EoFqhEkrt1R+tSqYbWveKDNni5H3jDbL8PhBClhRAF0b5hrzRdJ9DaNF04Em0Q1GDFMY8AT4UQnwgh8ggh9EKIGkKIeqmI+1+0vvyPhRD2pusiXgL+yIw2knntaX5tqXw/76GNv1gkpQxCG6P4FS2xXbC2HSllJLAGLfkckVIGJPHcNL9/UsonwFfAbCFER9Pz3dES4AOSTl6J401NDDXR3otHwCK0MZXEZ95KGqlEkru9CfwqpQyQUt6NuwFzgP5m39iWA9vRBo6voQ22OwKT0f7w76LNgvk0rmEhxBYhxKckIqU0oP2x+6DNanoALECbAmoVKWU02jTlTqbn/wS8IaW8mEltWHzt6Xxtyb6fiXwLjDd1n32UxD7L0WZnJdetlVQ7v6HNikuqWyvd/wdSyqlor28aWlfTdcAZbTA9zJp4UxlDTbSZcnGxT0adlWQYIaVaaldRlGeEEGWBi0BxKeXTTDrmILSzlCbJnAUp2ZTKxoqixDONN4wG/sisJAIgpVwkhIgBGgMqkeQw6oxEURQgfrD7HtrMp46mqdaKkiKVSBRFUZR0UYPtiqIoSrqoRKIoiqKkS64YbC9cuLB0d3fP6jAURVFylOPHjz+QUqZYsiZXJBJ3d3eOHTuW1WEoiqLkKEKIG9bsp7q2FEVRlHRRiURRFEVJF5VIFEVRlHRRiURRFEVJF5VIkrDh5C26fbuGfyc04KVv17LhpFXVxRVFUXIdlUgs2HDyFuPWnaFX2ArqiUv0ClvOuHVnbJZMVNJSFCUnyxXTf1Pru22XcIl5QC/HveiEpI9+D8diPNn+91k8Qiuj09mBTo/Q2SH0dgidHvR26HR6hM4edDqE3h69To/Q24HODr2dtp9Orz1HZ9p/27m7fLP5Ap/JFdTTxyWtvAB0r2XNYnOpFHIX1rwFPReDa+KF+xRFUVIvV9Taqlu3rkzNdSQeYzcx0W4RffS7sRfWrNWUdjFSjxGBA7EIAVLCDVmUMOFMHmcXpH1edA7O6B2dsXdywTGPC055XXFydkHn4Az2ppuDM9jneXbf3nTfIa/2r50TG/67jXHjaLrHbmO9XQf0XafbJlkpivJCEEIcl1LWTWk/dUZigVf+CHpF7k2QRKKkHeP0/8d77X3AaEAaDUhDDBiNSGMs0hgLxtj4+xgNYIxFmv4V0hh/XxhjQRrAaGDvpbu01P1HFRGAHRIjAongjrEABaJisQt/iIO8TR6i0Yso7IjCgSh0wpjMK3ieRNBe2pOHaISAl2O3M2X9FJzDB9C+Sf2MfgsVRclFVCKxYEaJfxDXnz9TG1b2BuUbJrVoXdp89u0a3oncjJ3QjqcXkuI8ZkSeb/l73KsAhEXFEhQSxY3QKIJCtNvDJ6E8efqEkJCnhIY+JTw0hKiIUBxkJM5EkYdonEQUzkThqo/BQUbRRhzDU9xEj0SPkc90v8E/v8HRsuDeDNybaje3shn6GhVFebGpRGJB+chzIGITbHMUsdr2DGYpaekwMrPEdkBLJHkd7cjraId74bzJtmU0SoIjYuKTTVBoZPzPf+w/wbuOf6M3JSwhIEraMzv2FT4qEQGXtsB/pqWy3colSixlMvx1K4ry4lCJxJIhBzLtUBmZtHQ6QcG8DhTM64BncdcEj1U98RUiNvFZlqSoeMQ7keMY3ncWNR3vwPX94L8fLm2C/5ZquxVwNyUVU3LJXzrVsSmK8uJSg+25RPD0Brg9vfjc9psOFegaM5knETE0r1yE4a0rUs+9IBiNcP88+B/QEov/AYgM1p5UwCNRYlED9oryIrJ2sF0lklxkw8lbfLftEreDIyjplocxHTzpXqsUoVGx/O57gwX7r/EwLJqG5QsyvHUlGlcohBBCe7LRCPfPaQnl+n64cQAin2iPFSyfMLHkK6ltV1ONFSVHU4nEjEok1gmPjmXFkZv8vNeP+yFR1C7rxvDWlWjpWeRZQoljNMC9c8/OWG4cNEssFbSE8vgGXN8LdQdB1+mZ/4IURUkXlUjMqESSOpExBlYfD2TeHj9uBUdQo1Q+hrWqRPtqxdDphOUnGQ1w7+yzMxb/AxAdoj1m5wQjT6uzEkXJYVQiMaMSSdpExxrZcPIWP+65yo2H4XgWc2VY64p09iqBPqmEEufvD+HkEu3aGoBar8PLc2wftKIoGcbaRJIja20JIaoKIeYJIdYIIYZmdTwvKgc7Hb3rlWHn6BbMeM0Hg5QMX3GSdj/sZe3xQGIMSVwUGXIXTi1/lkRAmwH28FrmBK4oSqbK9EQihFgkhLgvhDibaHtHIcQlIcRVIcTY5NqQUl6QUg4BegMpZkslfez0OrrXKsX2Uc35qX9tHO30/N/qU7T+fg/L/w0gKjZRGZm9U0EmSjJSwuLOEBOZeYEripIpsuKMZDHQ0XyDEEIP/Ah0AqoBfYUQ1YQQXkKIjYluRU3P6QYcAHZmbvi5l04n6OxVgs0jmrLgjboUdHbg0/VnaPndHn475E9kjCmhBB4BQ/TzDYTcgVVvQKyFxxRFybGyZIxECOEObJRS1jDdbwR8KaXsYLo/DkBK+a0VbW2SUnaxsP094D2AsmXL1rlxw6o17JVUkFKy/8oDZu+6wlH/xxRxdeS9ZuXp16As/5y/9/xUY8N22DgKqnWHVxeCXl0PqyjZWbYebLeQSHoCHaWU75juvw40kFIOS+L5LYEegCNwWkr5Y3LHU4Pttnf42kNm77rCwasPcbbXEW2QxBqf/W7lsdfzbQ8vukdugG2fQs2+8PJPoMuRw3SKkivktOq/lqYAJZnhpJR7gD22CkZJvYblC9GwfCGO33hMv18OJ0giABExBr7bdonuYz+A6HDY/bVW6r7L91rhL0VRcqzskkgCAfPKgKWB21kUi5IOdcoVIDrW8myu28ER2g/NP4LoUDg4Q1tHpd3/VDJRlBwsuySSo0AlIYQHcAvoA/TL2pCUtCrplodbcUnDTBFXR+0HIaDtlxATDodmg4MLtEx2op6iKNlYVkz/XQH4Ap5CiEAhxNtSylhgGLANuACsklJmfM12JVOM6eBJHnv9c9sfh0ez8bTpRFMI6DgFfAbAnm/h4KxMjlJRlIyS6WckUsq+SWzfDGzO5HAUG4hbvtd81ta7zT34+9Qdhi0/yYkbwYzrXAV7vQ66zdLOTP75XFsSuP67WRy9oiiplV26tpQXTPdapZ5bD75f/XJ8s/kCiw5e53RgMD/2r02xfE7QYz7ERMDmj7Q15n1Ur6ai5CRq7qWSaRzsdHzZrTqz+tbi/J2ndJl1gMPXHoLeHnothvKt4M8P4Nz6rA5VUZRUUIlEyXTdapZkwwdNyJfHjv4L/mX+Pj+knSP0WQZlGsDad+DS1qwOU1EUK6lEomSJysVc+fODJrSvVoxvNl9k6NIThBgdoN8qKO6llVK5tierw1QUxQoqkShZxtXJnp/61+azzlX558I9Xp5zkMtPBAxYB4Uqwoq+EHA4q8NUFCUFKpEoWUoIwbvNy7P8nQY8jYzl5TkH+fNyBLyxQVuyd1kvuH0yq8NUFCUZKpEo2UKD8oXYNKIpNUrlY+Qf//HlriCi+2+APG7w+ytw73xWh6goShJUIlGyjWL5nFj+bkPeburB4kP+9PkjgPuvrNaW6l3yMjz0y+oQFUWxQCUSJVux1+v4vGs15vSrxcW7IXRacpOTLX/TFsr6rRsEB2R1iIqiJKISiZItdfUuyV/DmuDmbM+rax6wqupsZHSIlkye3snq8BRFMaMSiZJtVSzqyp/DmtKpRgk+PiiZXHASMuw+/N4dwh5kdXiKkq1tOHmLJpN34TF2E00m72LDyVs2O5ZKJEq25uJox5x+tfi8azUW+hdilG4cxkfXtQH4iOCsDk9RsqUNJ28xbt0ZbgVHIIFbwRGMW3fGZslEJRIl2xNC8HZTD1a81xBfQ1Xei/4Q473zsKwnRIVmdXiKku1M3XqRiBhDgm1xi8vZgkokSo5Rz70gG0c05WnplgyNGo4h8ATG5X3gkT/82glC7mV1iIqSpUIiY5i7x4/bTyItPn7bwjpBGUElEiVHKerqxLJ3GlCu6Wt8GD0EbhzgzqzWGP19WTtjhE37gRUluwoOj+aHfy7TZPIupmy9iKOd5Y/2km55bHL8HJlIhBAthRD7hRDzhBAtszoeJXPZ63V82rkqhRr245uYfpTgIToh6RK7k+/X7VfJRMk1gkKi+HbLBZpM3sXMnVdoWL4Qfw1rwpRXvZ9bXC6PvZ4xHTxtEkemr0cihFgEdAXuSylrmG3vCMwE9MACKeXkZJqRQCjghLbeu5ILbT9/nyHiHrFSh50w4kAMH8tFTN5W9Lm1UBQlJ9tw8laCheLeaebBjYfhrDgSQIzBSFfvkrzfqgJViucDwLu0G5BwcbkxHTxt9nchpJQ2aTjJAwrRHC0JLIlLJEIIPXAZaIeWGI4CfdGSyreJmhgEPJBSGoUQxYDpUsr+yR2zbt268tixYxn7QpQsV3/sUvY5jsJJxMRvkxJ+iH2Vtz+bS35nxyyMTlEyRtwMrMSD5wLoVbc0Q1tWxKNwXpscWwhxXEpZN6X9smKp3X1CCPdEm+sDV6WU1wCEEH8AL0spv0U7e0nKY0B9WuRSY/P+jYhN+EXIiGC0/VoOTb2Kf7Np9GxRF4ck+osVJSf4btul55IIQNF8jkztWTMLInpedvkLKwXcNLsfaNpmkRCihxDiZ+B3YE4S+7wnhDgmhDgWFBSUocEq2UPrvP44itgE2/RCEuJYlDriIh329WDC1O/YePo2mX3mrSgZ4eytJ9xKYqbV/adRmRxN0rLLmu3CwrYk//KllOuAdck1KKWcD8wHrWsrXdEp2ZLb6H+f6zuO6weW9y/isHwgk4O/4fdVR3ht3/uM6VqLeu4FszpsRUnR8RuPmLPrKrsvBSGw/GFoqxlYaZFdEkkgUMbsfmngdhbFouQg3WuVsjiAKIpWwXXYXow7JvL64Tk0fXCR93/+gLLV6vNJxyqUL+KSBdEqStKklBzye8icXVfxvfaQgnkdGNPBk4J57Zn494UE3Vu2nIGVFtklkRwFKgkhPIBbQB+gX9aGpOR4do7oOk6CSm0ot34IG8UEpl3pS4cLHejbwJ2RbSpRyEUNsSlZS0rJrov3mbP7KicDginq6sj4LlXp16Aszg7aR3Qee7tMm4GVFinO2hJCHAN+BZZLKR+n+4BCrABaAoWBe8AXUsqFQojOwAy0mVqLpJST0nusOGrWlkLYQ/hrOFzaxBXX+rz+cCCh9oUZ2rICg5p4kMdBn3IbipJO5l2xJdycaFe1GEf9H3P+zlNKF8jDkBYV6FmnNE722eP30dpZW9YkkorAW8BrQFxS2S5z0OilSiQKoM0NPr4Yto7DYOfE3PyjmXajAiXyO/F/7T15pVYp9DpLw3WKkn5JTeMt4urAJx2r8rJPSez12WX+k8baRJJi1FLKq1LKz4DKwHJgERAghPhKCKFGLpWcQwio+xYM3oferQzD7n3Ov94bKe0i+Wj1KbrOPsD+K2qGn2IblgopglapoWed0tkuiaSGVZELIbyB74HvgLVAT+ApsMt2oSmKjRSpDO/sgMbDKXZ5OavEpyzu5ERIZAyvLzzCG4uOcOHO06yOUnlBREQbWHTgepKFFO8EW96ek6Q42C6EOA4EAwuBsVLKuMnL/wohmtgyOEWxGTtHaP81VGiD2DCUlntfY0/rz1ls7MLs3dfoPGs/PWuX5v/ae1I8v1OmhpbUlGYlZwmJjGHp4QAW7L/Gw7BoHPQ6og3G5/bLTtN408qaMZLycVecm23zkFJet2lkGUiNkSjJCn+kDcRf3AjlW/GkwyxmHw1lie8NdDp4t1l5BreowI7z92z+AR/Xj+4S84A5DrMZFj2CUPtCfNvDSyWTHCI4PJpfD/rz68HrPI2MpXnlIgxrVZHbpsWlEk/jzc7/txk52H5CSlnbQuN10hljplGJREmR2UA89nng5TncLNqKqdsu8fep27g46omKNRJjePb3ktyHQKzBSHiMgfAoA+HRsYRHGwiPNhAWHUtEtIGwqFgiYgyERRmIiI4lzPT4hpO3iIgx8D+7RfTX72SpoQ0TYgdRyi0PB8e2zsQ3RElJ4jPHwS3Kcys4gqW+NwiLNtC+WjGGta4YX0DR0nOy+9lmuhOJEKIKUB2YCowxeygfMEZKWT0jAs0MKpEoVgu6DOvegTunoM5b0OEbTt2Lptc8X4vdEg52OqoWd41PFOGmpBAd+/y+SRECnO315HGw40FoJFW4wZ+OE3AUsURIe5pHzeQBblyf3CUjX6mSDknNwALoVjNhJd6cLCOKNnqiFUx0A14y2x4CvJu+8BQlmypSGd7eAbu/hoOz4MZBar66gBgLSQQgOtaIm7MDJd30ODvY4eygx9lRj7O9HXkdzbY5mH521ONsL8gXHUTesADyhNzA/ok/4vE1eHSdiNgr5OFZDSUnYljh8DVTdYMwxrZHZ2efWe+EkowkCym6OjKrb60siChrWdO11UhK6ZtJ8diEOiNR0uTaXlg/GMIeMEfXj+9D21GYJ/FjF0G4Jd3lZIiFJzfh0TXT7To8vv7sZ4NZwT29AxTwgILluR7lQmn/ddiLZx9SUmpnLcG6Ajj5vIpTrT5Quq62Ucl0V+6F0O6HfRYfE/BCnTmm+4xECPGxlHIq0E8I0Tfx41LKEemMUVGyt/ItYOgh+Gs4wy7+Rm3HE9w2uFFPXGK43Tqmibf4sokjXN5uljBMt+AbYDSrTGyXBwqWh0IVoVJ77ee4W76SoNOuZPbYOBpDgA7ks0Ri1Nlzt2A9zgbF0vLEb3BiAbiVA6+eUKMnFKuW2e9MrnT21hPm7LrK1nN3c0QhxcyUXNfWBdO/6qu8kns5F4TXlsKJJdTf9DF6EYkABtjt4HV2InaafZw4uEJBDyjuBdVeTpgsXItbdwYReAS9jEmwSS9jKGUfwsP3NvLS7/vxCTvAKP0pShz4AbH/eyhaHbxehRqvQgH3DH35ilaJd/auq+y5FISrkx3DW1ekWD4nJm3K3oUUM5M1XVu1pJQnMykem1BdW0qGWPM2nF2L9l1UByW8ocGQZ8kib2GbdzcFh0czetUpdl28z4AaTnzucRnHi+vh5r/aDqXra2cq1V8Bl6I2jeVFkng21UftK1M0nxOzd13h8LVHFMzrwNtNPXi9UTnyOdlbfE52n4GVFhk5/Xc3UAJYDfwhpTyXMSFmHpVIlHQLuQsza0Ks2VXIdk4w8jS4FsvUUIxGydy9fny//RIVirgwd0BtKto/0pLc2bVw7ywIHXi0AK9eULUrOOXP1BhzEkszsITQxqaKujryXvPyCSrx5iYZlkhMjRUHeqMVbswHrJRSfp3uKDOJSiRKum0cDSd/B0P0s216B6j1OnSdniUhHbz6gBErThIRY2DKq968VLOk9sD9C3BmDZxZrY3V6B2hcnttPKVyB+06GSVek8m7LK5CmD+PPf9+2ibbVOLNChlWtBFASnlXSjkLGAL8B0xIZ3yKkrMEHkmYREC7H3gka+IBmlQszKYRzahWIh/DV5zky7/OadevFK0KbT6HkafgnZ1aocqAf2H1m/BdJVg/BK7s0GaWgXa29WsnCLmXZa8lq8QajEkuZfs0IiZXJ5HUsKZrqyramUhP4CHwB7BWSnnf9uFlDHVGorzIYgxGpmy5yIID1/Ep48ZP/Ws/P3vIaIDr++DsGjj/N0Q9AefCUL27lkgubdYuwMyis6vMFh1rZN2JQH7a40fAo3CL+6hqAhk7RnIYWAGsllLmyOVvVSJRcoPNZ+7w8ZrT2OsFM/vUonnlIpZ3jI2CK/9oSeXi5mfXtNg5wsgzmT7mk5kiYwysOnaTeXv8uP0kEq9S+anvXoDlRwKIiHl20Wl2r4GVWTJ0jCS7EUI0A/qjTV+uJqVsnNz+KpEoucW1oFCGLj3B5fshjGxTiRGtK6FLbrGuP4fBf8ufXbfi3hze/OuFu9gxPDqW5f8G8PO+awSFRFGnXAGGt65Ii8pFEELkihlYaZERtbZWSSl7CyHOkPDaGwFIKaV3GgNbhFZ65b6UsobZ9o7ATLSldhdIKSdb0VZ3oJiU8ufk9lOJRMlNwqNjGb/+LOtO3qJ55SLMeM2Hgnkdnt/R0kw0gIptofs8cEnijCabM08KxfM7UbusG77XHvEoLJrGFQoxrHVFGpUvhHjBkqUtZEQiKSGlvCOEKGfpcSnljTQG1hwIBZbEJRIhhB64DLQDAoGjQF+0pPJtoiYGxY3PCCFWAe9IKZNdhUglEiW3kVKy4shNvvzrHIVdHPixf21qlS2QcCdLM9GEXpv36lwQXpqpTR3OQZIqpli1uCtfv1KDOuXUoq6pke5ZW1LKO6Yf35dS3jC/Ae+nNTAp5T7gUaLN9YGrUsprUspotAH9l6WUZ6SUXRPd4pJIWeBJUklECPGeEOKYEOJYUJBaPlXJXYQQ9GtQlrVDG6PTCXr/7MsSX38SfHG0NBNNGqBQBa1sy8r+sH4oRD7J1NjT45vNFywWU3waGauSiA1Zc4VNO+CTRNs6WdiWHqWAm2b3A4EGKTznbeDXpB6UUs4H5oN2RpLeABUlJ/IqnZ+Nw5syetUpJvx5jmP+j/m2hxd5He1gyIGknxgbDfumwv7vwX8/dP8JPJpnXuCpYDRK9ly+z68H/bkfEmVxn9tJTPFVMkaSZyRCiKGm8RFPIcRps9t14HQGx2GpszLZD38p5RdSykMZHIeivHDcnB1Y8EZdxnTwZOPp27z840Gu3g9J/kl2DtB6PAzarl14+dtLsGUsxGSfD+SnkTEsOnCd1t/vYdDiY1y6G4Krk+Xvxrm1mGJmSe6MZDmwBW2MYqzZ9hApZeKuqfQKBMqY3S8N5MipxoqSHel0gg9aVaRWGTdG/HGSbnMOMvlVb4xGmfxspTL1tDOXHV/Av3PBbye8Mg9KZd0CqX5BoSw55M+a44GERRuoXdaN0e096Vi9OJvP3LG4nG1uLaYYFB7EmH1jmNZiGoXzFLbZcZIbbM8npXwqhLDYsZieZCKEcAc2mg2226ENtrcBbqENtvfLqLpearBdUZ65+ySSYctPcOzGY/Q6gcFo3fLB+O2CDR9A6D1oPgaafwT6zFloy2iU7L0cxOJD/uy9HISDXkfXmiUY2Ng9wVK2kDuKKVrrf4f/x+pLq+nl2YvPG36e6udnxKytjVLKrqauLEnC7icppSyf6qi0dlcALYHCwD3gCynlQiFEZ2AG2kytRVLKSWlp3xKVSBQloRiDkVoTtxMa9fzAdLJXdEcEw5ZP4PQfUMIHesyHIrb7th8SGcOa44Es8b3B9QdhFHV1ZEDDcvStX5Yiro42O25OV2dpHaITT6QAHPQOHB8e2AOpAAAgAElEQVRw3Op2XugLElNLJRJFeZ7H2E1JDkR+84oXPmXcqFzMBTu9haHU83/C36MgOgzafqmV09dZVbrPKteCQlnie4PVx24SFm2gVlk3BjZ2p1ONEjjYZdxxXlRB4UF8tPcjTtw/AYCT3ok2ZdvwUb2PUtXFlRFrtsc11AT4T0oZJoQYANQGZkgpA6yORlGUbKekWx6LBQt1Aj5dfwYAZwc93qXz41OmALXKulGrjBtF8zlpC3eVaQh/j4Rt47RaXd1/AreyaY7HaJTsu6J1X+25FIS9XtDVuyRvNnbHp4xbyg0o8Y7fOx6fRBx0DkQZosjrkNdm4yTWTP+dC9QUQtQEPgYWAr8DLWwSkaIomWJMB0+LA9PfvFKD2uUKcDIgmP9uBnMy4DELD1wjxqCdv5Ryy4OPKanUavwj3pU2Yv/Pp/BTY+g0GXz6p6rESmhULGuPB/LbIX+uPQijiKsjH7atTN8GZSjq6pThr/tF9+fVP5lwaAL5HfLTtlxb+lbpy+rLq3kQ8cBmx7SmaOMJKWVtIcQE4JZpPOOElLK2zaLKYKprS1Ess3ZgOjLGwLnbTzkZ8NiUXILjz2bs9YKWRcP5LGY27qH/Ee7Rnjw95iDMij9aOo5PGTd+8/Vn9bFAQqNi8SnjxltNVPdVeqy+vJqJvhNpWKIhM1vNxNneOV3tZWT1373AVuAtoDkQhNbV5ZWuCDORSiSKkvHuh0TyX0AwJ01nLWcCH9PHsJGP7VYRLvKwrMhoZJWuhEXHsviQP5Fm1XV1AoxSS0JdvErwZmP350u4KKmy9PxSphydQrNSzfih1Q846tM/GSEjE0lxoB9wVEq531SapKWUckm6o8wklhJJTEwMgYGBREZGJvEsJTdwcnKidOnS2NtnzjTWF5nBKLl8L4Rr545R89gnlI68zFpDM76MeZMQnv9m7Opkx87/a6G6rzLAgjMLmHliJm3KtuG75t9hHzctO+QurHkLei5O0/IAataWGUuJ5Pr167i6ulKokKoCmltJKXn48CEhISF4eHhkdTgvFkMM7PsOuW8at41ujIkZzBVjKeY4zGZY9AiCcEMA1yd3yepIczQpJXNPzWXuqbl0cu/EpGaTsNeZfSnaOBqO/5rmRcsybKldIUQPIcQVIcQTIcRTIUSIECLZars5QWRkpEoiuZwQgkKFCqmzUlvQ20OrTxFv/4NB58hyh29Y4jCZeuISw+3WAapsSXpJKfnhxA/MPTWXlyu8zLfNvk2YRJ7egeOLQRrhv2U2XUrZmhGtqUA3KWV+KWU+KaWrlDKfzSLKRCqJKOp3wMZK1+FU579ZbWxFVd1NdELymn4PZexDcm3ZkowgpWTK0Sn8evZXelfuzcQmE9HrzNaXNxpgycvPFiyTRtg7xWbxWJNI7kkpL9gsAkVRXmgv1atInfLFiEH7oHMUsWxx/R/dPYwpPFOxxCiNTDw8kWUXljGg6gDGNxyPTph9lBtiYdWb8OCS2bZom56VWJNIjgkhVgoh+pq6uXoIIXrYJJpcRq/X4+PjE3/z9/fP6pDYs2cPXbumvJhRy5YtiRt36ty5M8HBwQDMmjWLqlWr0r9/f6Kiomjbti0+Pj6sXLnSpnEr2VjIXcoHbsCeZ9eruIQHwpx6sH+6VrJeSVFQeBBvbnmTMXvHsObyGt71epeP632c8Kw6NlobXL/4t7ZImTkbnpVYc0FiPiAcaG8eErDOJhFlU7YoBJcnTx7++++/VD8vNjYWOztr/utSZjAY0Ov1Ke+YjM2bN8f//NNPP7FlyxY8PDw4fPgwMTExqXqNGfnalGxi71TtQ8yczh5cisLOr7Rvyp2mQsU2WRNfDjH31Nz4q9WH+QxjcM3BCXeIiYRVr8OV7eBaAkLuJHzcEK0tZmYDKf7FSinfssmRc5DEy3feCo5g3DqthERGVxWNjIxk6NChHDt2DDs7O6ZPn06rVq1YvHgxmzZtIjIykrCwMKpUqULHjh3p1q0br7zyCgUKFGDRokUsXLiQ69ev8/XXX9O9e3du3rxJZGQkI0eO5L333gPAxcWF0aNHs23bNr7//ntCQ0MZNWoUhQsXpnZty9eZRkRE8NZbb3H+/HmqVq1KRMSz0hru7u4cO3aM8ePHc+3aNbp168aAAQP45ZdfCAoKwsfHh7Vr1xIcHMzo0aMJDQ2lcOHCLF68mBIlStCyZUsaN27MwYMH6datG2+88QZDhgwhIECrwjNjxgyaNGnCl19+SUBAANeuXSMgIIBRo0YxYsQIAJYsWcK0adMQQuDt7c3vv/9OUFCQxXaUTGZpJUZjDDjlg/5rYMvHsLQHVO0GHb4BtzKW28mlLBVgnPPfHOafmf+sAGN0GKzoC9f3QdcZUDdzP7atqbVVGa1MSjEpZQ0hhDfa4PvXNo8uk3z19znO3056ItrJgGCiDQm/UUXEGPh4zWlWHLFccqxayXx88VL1ZI8bERGBj48PAB4eHqxfv54ff/wRgDNnznDx4kXat2/P5cuXAfD19eX06dMULFiQP/74g/3799OtWzdu3brFnTvat48DBw7Qp08fABYtWkTBggWJiIigXr16vPrqqxQqVIiwsDBq1KjBxIkTiYyMpFKlSuzatYuKFSvy2muvWYx17ty5ODs7c/r0aU6fPm0x4cybN4+tW7eye/duChcuTIMGDZg2bRobN24kJiaG119/nT///JMiRYqwcuVKPvvsMxYtWgRAcHAwe/fuBaBfv358+OGHNG3alICAADp06MCFC9ow3cWLF9m9ezchISF4enoydOhQLl++zKRJkzh48CCFCxfm0SNthYORI0cm2Y6SiZJbiRHA3Rd8Z8O+7+HqDq1EfaNh2uJaCn++/Cevb3k9vsSJeQFGQFsKeVlvLWG/Mg9q9sn0GK3pQ/gFGAP8DCClPC2EWA68MIkkJYmTSErbrWWpa+vAgQMMHz4cgCpVqlCuXLn4RNKuXTsKFtSWh2nWrBkzZszg/PnzVKtWjcePH3Pnzh18fX2ZNWsWoI1XrF+/HoCbN29y5coVChUqhF6v59VXXwW0D2YPDw8qVaoEwIABA5g/f/5zse7bty/+27+3tzfe3t6peq2XLl3i7NmztGvXDtC61EqUKBH/uHkC27FjB+fPn4+///TpU0JCtBX9unTpgqOjI46OjhQtWpR79+6xa9cuevbsSeHCWkG6uPcoqXZcXV1TFbtiY/ZOWvLwfg22jlPdXWbCY8L58tCXPIh4gEBgr7dPWIAx/JF2Nnf3DPRcBNVfyZI4rUkkzlLKI4mmScbaKJ4skdKZQ5PJuyxWSS3lloeVgxtlaCzJXSCaN2/eZ8cuVYrHjx+zdetWmjdvzqNHj1i1ahUuLi64urqyZ88eduzYga+vL87OzrRs2TL+egknJ6cE4yLWToFNz1RZKSXVq1fH19c3xddmNBrx9fUlT57nrzNwdHxW9kGv1xMbG4uU0mJsybWjZENuZaHPMrjyj+ruAkKjQ3l/5/ucCjpFtULV8CrsRa/KvZ4VYAwN0qb4PrwKry0Dz45ZFqs1s7YeCCEqYFpDXQjRE7iT/FNeLGM6eJLHPuGAtK2W72zevDnLli0D4PLlywQEBODpafk4jRo1YsaMGTRv3pxmzZoxbdo0mjVrBsCTJ08oUKAAzs7OXLx4kcOHD1tso0qVKly/fh0/Pz8AVqxYkWJcZ8+e5fTp06l6XZ6engQFBcUnkpiYGM6ds7wAZvv27ZkzZ078/ZQG69u0acOqVat4+PAhQHzXVmrbUbKJSu1gqK+2ZvyVf+DH+rludteTqCe8u/1dzgSdYWrzqazsupLxDcfjWdCT8Q3HM6POx/BrJ3h0DfqtzNIkAtYlkg/QurWqCCFuAaOAoTaNKgVCiGpCiFVCiLmmxGZT3WuV4tseXpRyy4NAOxNJcjnSdHr//fcxGAx4eXnx2muvsXjx4gTfws01a9aM2NhYKlasSO3atXn06FF8IunYsSOxsbF4e3vz+eef07BhQ4ttODk5MX/+fLp06ULTpk0pV66cxf2GDh1KaGgo3t7eTJ06lfr166fqdTk4OLBmzRo++eQTatasiY+PD4cOHbK476xZszh27Bje3t5Uq1aNefPmJdt29erV+eyzz2jRogU1a9Zk9OjRaWpHyUbiuruGHYEKrbXurrmN4OrOrI7M5h5FPuLtbW9z6fElprecTgf3Dgl3eHxDSyIhd+H1dVChVdYEasbqWltCiLyATkoZkq4DCrEI6Arcj1uz3bS9IzATbandBVLKycm08X/AEVMRyb+klN2SO6alWlsXLlygatWq6XglyotC/S7kAFd2wJYx2jfwF7i7Kyg8iHe3v0tgaCAzW82kSSmzWYYhd7WZWU9vQ2wEDFgPpevYNJ6MrLU1UggRdy3JD0KIE0KI9ik9LxmLgQTnYUIIPfAj0AmoBvQ1nXV4CSE2JroVRVtYq48Q4jugUDpiURQlJ6jUFt4//EJ3d90Nu8tb297idthtfmrzU8IkArB1LNw+ARGPYeAmmyeR1LCma2uQlPIp2gWJRdHWJUnybCElUsp9wKNEm+sDV6WU16SU0cAfwMtSyjNSyq6JbvdNtw+AsYDFZb+EEO8JIY4JIY4FBQWlNVxFUbILO8cXtrsrMCSQgVsH8iDiAT+3+5n6JRJ1HV/dAefWP7uft2jmBpgCaxJJ3HSYzsCvUspTZtsySingptn9QNM2ywEJ4S6EmA8sAb6ztI+Ucr6Usq6Usm6RIkUyNFhFUbJQ3Oyu/mu1K+aX9oCVr0PwzZSfmw35P/Fn4NaBhESHsKD9AmoVrZVwhxu+sNz82hBp0wKMaWFNIjkuhNiOlki2CSFcgYyutmYpMSU5eCOl9JdSviel7C+lTOFqJ0VRXkjJdXeF3DUNSNuudHpG8Av2461tbxFtiGZRh0XUKFwj0Q67YEl3MJpdcWHjAoxpYU0ieRutC6melDIccEDr3spIgYD5yFlp4HYGH0NRlBdNUt1df42AgMPZ7ps7aAPqA7cOxPe2L29t1T5Kf+34K54FE03zv7ARlr8GDs6gT3TJn43LwqdWiolESmmUUp6QUgab7j+UUqbuIoKUHQUqCSE8hBAOQB/grww+hqIoLyrz7i5DNFzZlikLOqXFvNPzOHHvBB/s/AAHvQOLOy6mgluFhDudWgmr3oDi3uBaXFtx0pwNCzCmhTVnJBlKCLEC8AU8hRCBQoi3pZSxwDBgG3ABWCWltHy12gvk7t279OnThwoVKlCtWjU6d+4cXw4ltRYvXszt2xl7Eufu7s6DB9pchsaNG8dvHzNmDNWrV2fMmDEEBQXRoEEDatWqxf79+zP0+IqSapXaamcmcR9txphs8829ztI6eP3mxapLq5BIYowx3Au/R4+/Eq3KcXQhrB8M5RrDGxvgfV/48snzt5RqmGWiJBOJEMImi1hLKftKKUtIKe2llKWllAtN2zdLKStLKStIKSfZ4tjploH9rlJKXnnlFVq2bImfnx/nz5/nm2++4d69tLWdlkQSG2t9pRvziwd//vlnTpw4wXfffcfOnTupUqUKJ0+ejL8YMiUGgyHlnRQlLULuwqk/iB/GNRrg5O/Z4qxka4+tNCjeIP6+o96RLh5d2Pbqtmc7HZwJm0ZDpfbQfzU45oy6cMmdkawBEELk/Ll1GWXv1Azrd929ezf29vYMGTIkfpuPj0/8h/F3331HvXr18Pb25osvvgDA39+fqlWr8u6771K9enXat29PREQEa9as4dixY/Tv3x8fHx8iIiI4fvw4LVq0oE6dOnTo0CG+OnDLli359NNPadGiBTNnzkwQ08OHD2nfvj21atVi8ODBCep+ubi4ANCtWzfCwsJo0KABU6ZM4eOPP2bz5s3xx92+fTuNGjWidu3a9OrVi9DQUEA7u5k4cSJNmzZl9erV+Pn50bFjR+rUqUOzZs24ePEiAAMHDmTEiBE0btyY8uXLs2bNmvgYpk6dipeXFzVr1mTs2LEASbaj5FKW1j4xRMM/n2dNPGYuPb7E0btHAXDQORBtiH5WfFFK2DUJ/pmgFV58bSnY55wacckVbdQJIb4AKgshRid+UEo53XZhZbItY7XqmcmJjYbbx7Rf0uO/avvrkylzXdwLOiV9uc3Zs2epU8fyBUXbt2/nypUrHDlyBCkl3bp1Y9++fZQtW5YrV66wYsUKfvnlF3r37s3atWsZMGAAc+bMYdq0adStW5eYmBiGDx9uVcl2c1999RVNmzZlwoQJbNq0yWIV4L/++gsXF5f4ulXFihXj2LFjzJkzhwcPHvD111+zY8cO8ubNy5QpU5g+fToTJkwAtHIsBw5op+Nt2rRh3rx5VKpUiX///Zf333+fXbt2AXDnzh0OHDjAxYsX6datGz179mTLli1s2LCBf//9F2dn5/h6Wu+9916S7Si5kKW1TwDObYDO07Q1ULLAroBdfLT3I5ztnWlbti0Dqg14VnxRStj2KRz+CWoNgJdmgS59i81ltuQSSR+gu2mfnHF+ZUtPArT/cND+DQ6AQhVtcqjt27ezfft2atXS5pOHhoZy5coVypYti4eHR/waJnXq1LG4PG9qSrab27dvH+vWaQtfdunShQIFCqQq7sOHD3P+/Pn4xaOio6Np1OhZdeS444aGhnLo0CF69eoV/1hUVFT8z927d0en01GtWrX4rr4dO3bw1ltv4ezsDGil4lNqR8mFLI0b+O2Gpa/C2rehz4rnZ0DZ2Db/bYzdN5aqhaoyt+1c8jvmB2B8w/Fa19tfw7XutwZDtdIvukwfuk63JN9RKeUlYIoQ4rSUcksmxpT5kjlzALR+15k1eXZpi4TIYK3+v2uxNB2yevXqCbptzEkpGTduHIMHJ1xK09/f/7ky6uYrFZo/39qS7Ymlt1R8u3btkqwgHHdco9GIm5tbktV4zV9jXPeapVLxKbWjKIBW1LDLNNj4IWz/DDpl3uD7335/M/7geHyK+PBjmx9xcXDRPk/WvAU9foHtn8O5ddoU5lafQTr+/rKSNanvkBBiely5ESHE90KI/DaPLDux1O+aznncrVu3Jioqil9++SV+29GjR9m7dy8dOnRg0aJF8eMLt27d4v79+8m25+rqGr/4U2pKtpszLxW/ZcsWHj9+nKrX1LBhQw4ePMjVq1cBCA8PtzgLLV++fHh4eLB69WpASxKnTp1Ktu327duzaNEiwsPDAa1UfFraUXKpuoOg4Qfw7zw48kvK+1srmQk4ay6v4bMDn1GvWD3mtp2rJREwjbX6wqKOWhJp+5V2UWUOTSJgXSJZBIQAvU23p8Cvtgwq27HU75rOedxCCNavX88///xDhQoVqF69Ol9++SUlS5akffv29OvXj0aNGuHl5UXPnj3jk0RSBg4cyJAhQ/Dx8cFgMFhdst3cF198wb59+6hduzbbt2+nbNmyqXpNRYoUYfHixfTt2xdvb28aNmyY5OD3smXLWLhwITVr1qR69er8+eefybYdtz593bp18fHxYdq0aWlqR8nF2v8PKneCLZ9otasywq6vtRImGz/U1kv32w1Xd7B87wS+8v2KJvkrMadUR5wvb4Mza+Df+XBiidY9/uQmtPkCmo7KmFiyUIpl5IUQ/0kpfVLalp2pMvJKctTvQi4SFaqdCQTfgLe3Q9E0/r9Hh8G+aXDg+TlHv+Z3ZXrBArQKC2fa/QckOSVHp4faA6Fr9p23ZG0ZeWtGnSKEEE3jaloJIZoAz3fMK4qiZHeOLtDvD/ilDSzvDe/sApdUFHWNiYRji7QEEhaEViZQgs4OWbEdP5eqwI/+f9GxaH2+qTEYeztHEHotaYQ/ht+7g8E0IcRo0K68b/FJmsdaswtrEskQYInZuMhj4E3bhaQoimJD+UtD3xXwa2f4ox+8+be2ImNyYqPh5BLY9z2E3IYyDSHySXyXtzTGMuvhERbEnKNbhW5MbDwRfeIpvBtH81wt2rix1mx8VmINa2ptnZJS1gS8AW8pZS0b1NpSFEXJPKVqwyvztHHOPz94NrU/MUMsnFwKc+rApv/Tanq9uRGKVQcgSK9jYPGifFWoAAvy5aWnYyn+1+R/zycRsMlYa3Zh9YRq0+JWiqIoL4bq3eHRBNg5EQpXhjpvatNyey6GvIXh7DrY8y088oOStaDrD1ChjTa7ats4MEQzt1ABjjs5cjyPE/2fhPCJQwhCJPH9PBvVxspomXtljqIoSnbSdDQ8uAp7voEbh7QSSBuGauuiB12AotWhz3Lw7Jxgem4dlwiiPRLOalyW35XV+kiOZ/ZryAaS7doSQuiEEI2T20dRFCXHEgJemgml6sH1PdqYhd9OiI3ULjgecgCqdHnuGo+/u/9NCedn1SKc9E7PF2DMRZJNJFJKI/B9JsWS62RkGXlr7Nmzh65du9qk7eDgYH766SebtK0oNmXnAEXMFpUSeijfCmq8arFcSbQhmqlHp3InXCuE6qB3IMoQ9awAYy5kzQWJ24UQr4r01M54QcStbPYg4kG628roMvJZTSUSJccKuQtnzcoVSQOcWm7xavXI2EhG7h7JzoCdVC5Qmdc8X2N55+X09uzNw4iHmRh09mJNIhkNrAaihRBPhRAhQohMG3gXQpQXQiwUQqxJbltmiFvZbO6pueluK7ky8lJKxowZQ40aNfDy8mLlypWAdkbRokULevfuTeXKlRk7dizLli2jfv36eHl54efnBzy7yr1Zs2ZUrlyZjRs3Pnf8sLAwBg0aRL169ahVq1b8FeHTp09n0KBBAJw5c4YaNWrElyWJc+7cOerXr4+Pjw/e3t5cuXKFsWPH4ufnh4+PD2PGjAGSLoVfpUoV3nzzTby9venZs+dz7StKprKyBFJ4TDjDdg7j4K2DfNHoC9Z2W8v4huPxLOjJ+IbjmdFqRiYGnb2kONgupUxz5V8hxCKgK3BfSlnDbHtHYCagBxZIKZOsmiilvAa8bZ40LG1LjylHpnDxUdLrWBy/dxxpNv971aVVrLq0CoGgTjHLpeCrFKzCJ/U/SbLN5MrIr1u3jv/++49Tp07x4MED6tWrR/PmzQE4deoUFy5coGDBgpQvX5533nmHI0eOMHPmTGbPns2MGdovs7+/P3v37sXPz49WrVrF17+KM2nSJFq3bs2iRYsIDg6mfv36tG3bllGjRtGyZUvWr1/PpEmT+Pnnn+Mr7saZN28eI0eOpH///kRHR2MwGJg8eTJnz56NL6CYXCn8S5cusXDhQpo0acKgQYP46aef+Oijj5J8rxTFpqyYlhsaHcr7O9/nVNApJjWdxEsVXsrkILM3q2ZtCSG6Ac1Nd/dIKZ//imvZYmAOsMSsLT3wI9AOCASOCiH+Qksq3yZ6/iApZfLVCjOBV2EvAkMCeRz1GIlEICjgVIAyLmVscrwDBw7Qt29f9Ho9xYoVo0WLFhw9epR8+fJRr169+JLwFSpUoH379lqMXl7s3r07vo3evXuj0+moVKkS5cuXf67m1fbt2/nrr7/ia1ZFRkYSEBBA1apVWbx4Md7e3gwePDi+JLy5Ro0aMWnSJAIDA+nRoweVKlV6bp/kSuGXKVMmvt0BAwYwa9YslUiUrJPCtNwnUU8YumMoFx5eYGrzqXRw75BJgeUcKSYSIcRkoB6wzLRppKlkytiUniul3CeEcE+0uT5w1XRWgRDiD+BlKeW3aGcvmS65M4c4E30nsubyGhz0DsQYYmhbri2fN0z7qmsplZFPinmJdZ1OF39fp9MlWDo38ZBW4vtSStauXYunpyeJXblyBRcXlySX7u3Xrx8NGjRg06ZNdOjQgQULFlC+fPnn2k+qFH5KsSlKdvEo8hGD/xmMX7Af01tOp1XZVlkdUrZkzRhJZ6CdlHKRlHIR0NG0La1KATfN7geatlkkhCgkhJgH1BJCjEtqm4XnvRdX+j4oKCgd4WoeRT6it2fvDBtYS66MfPPmzVm5ciUGg4GgoCD27dtH/fr1U9X+6tWrMRqN+Pn5ce3atecSRocOHZg9e3Z80jp58iQAT548YeTIkezbt4+HDx9aTHbXrl2jfPnyjBgxgm7dunH69OkEZezj2k+qFH5AQEB8ifsVK1bQtGnTVL02RckMQeFBDNo6iOtPrjO79WyVRJJh7QWJbsAj08/pXYvE0tfPJL+CSykfotX7SnabhefNB+aDVv039WEmZD6QNr7h+PQ2F19GftSoUUyePBknJyfc3d2ZMWMGzZs3x9fXl5o1ayKEYOrUqRQvXjxV65F7enrSokUL7t27x7x583BySlhL6PPPP2fUqFF4e3sjpcTd3Z2NGzfy4Ycf8v7771O5cmUWLlxIq1ataN68OUWLFo1/7sqVK1m6dCn29vYUL16cCRMmULBgQZo0aUKNGjXo1KkT3333HRcuXIhfIdHFxYWlS5ei1+upWrUqv/32G4MHD6ZSpUoMHTo03e+nomSku2F3eWf7O9wPv8/ctnOpV7xeVoeUrVlTRr4vMBnYjZYEmgPjpJR/WHUArWtrY9xguxCiEfCllLKD6f44AFPXlk3ktjLyAwcOpGvXrvTs2TOrQ3mOv78/Xbt25ezZs1kdSrwX+XdBSb3AkEDe2f4OT6KeMLftXHyK5pgVMzJchpSRN107cgBoiDZOIoBPpJR30xHbUaCSEMIDuIW2Nny/dLSnKIqSbkHhQQzfNZy74XeJNcayoMMCqheqntVh5QjJJhIppRRCbJBS1gH+Sm3jQogVQEugsBAiEPhCSrlQCDEM2IY2U2uRlDLldWAVqy1evDirQ0iSu7t7tjobUZQ4U45O4dzDczjqHVnWeRmeBZ+fiKJYZs0YyWEhRD0p5dHUNi6l7JvE9s3A5tS2l9GklGrGUC6XUteu8uKrs7QO0WbXkUQZouj5d08c9A4cH5AbSzCmnjWztloBvkIIPyHEaSHEGSFEjl+PxMnJiYcPH6oPklxMSsnDhw+fm4ig5C7TW0zHTjz7Tp3bCzCmhTVnJJ1sHkUWKBg5L+UAABdcSURBVF26NIGBgWTE1GAl53JycqJ06dJZHYaSRU7cO8En+z/B0c4RQ4wBe719ri/AmBYpDbbrgE3m5U1eFPb29nh4eGR1GIqiZJHDdw4zYtcIiuctTqm8pSjlWopelXux+vLqDCnMmpukNNhuFEKcEkKUlVIGZFZQiqIotrQ/cD+jdo+iXP5yzG83P8HZR0ZcJ5bbWNO1VQI4J4Q4AoTFbZRSdrNZVIqiKDayM2AnH+39iEpulZjfbj5uTm5ZHVKOZ00i+crmUSiKomSCrde3Mnb/WKoXrs7ctnPJ55Avq0N6IVhTRn6vEKIcUElKuUMI4Yx2/YeiKEqO8efVP5lwaAK1itbixzY/ktc+b1aH9MJIcfqvEOJdYA3ws2lTKWCDLYNSFEXJSKsvr2b8wf9v797jqqrSBo7/Ho7cRFEBLwli2AWhzPJCF9KxwqYyp5kuZpeZt8mmcd56p8snM199x6YyzXGczMxbmTWVldbMmF2cqTHRdBJwLEvEiMbEShEINUVuz/vHASIEPHDOYR/g+X4+fGKvzt77WXLOfs5aa++1ppHSJ4WFaQstifiYJ8+R3AGkAgcBVPUzoFeTexhjTIB4MftFHtr8ECPjRvLkJU8S3inc6ZDaHU/GSI6palnNE+Ai0okmZus1xhinFRwpYFL6JIb0GsLS7UtJi09j9sjZBLuCnQ6tXfIkkawXkf8FwkVkNPDfwBv+DcsYY1pu0UeLyNqXRda+LC5PuJxHL3yUTkGerpphmsuTaeSDgAnApbhn/12Le531NtMqaWgaeWNM+1N/3qwaNm9Wy3g6jfwJx0hUtUpVl6rqdap6bfXvbSaJGGM6jrd/9jYDun2/7LPNm9U6PBlsN8aYgKeqrMhZQV5JHuBuhdi8Wa3DOg2NMe3Cgm0LeHr70/Tt0pcLYy9k3OnjbN6sVmKJxBjT5i3ctpDFHy/m6tOuZvr50wkSd2eLzZvVOlrUtSUit/s6kCbONUBEnhGRVXXKkkRkkYisEpHftFYsxpjAs/ijxTz10VNcdcpVP0gipvW09F/co2UFRWSZiOwXkU/qlV8mIjkikisiDzR1DFXNU9UJ9cqyVXUiMA444R0Fxpj26entT/PkticZO2Asv7/g95ZEHNKif3VVXXziVwGwHLisboGIuIAFuBfMSgZuEJFkERkkImvq/TT6BL2I/ATYCLzXkjoYY9q2Zz95lnlb53FFwhU8nPowriCbAtApnsy1dZeIRIrbMyKyVUQu9eTgqpoOFNUrTgFyq1saZcDLwFWqul1Vr6z3s7+JY69W1QuAmxqJ+3YRyRSRTFsF0Zj25blPn2Nu1lwuP/lyZlw4w5KIwzxpkdyqqgdxP5DYE/glMMuLc8YCe+ps51eXNUhEokVkEXCOiEypLhslIk+IyGLgrYb2U9UlqjpMVYf17NnTi3CNMYHkhR0vMCdzDpf2v5RHR9gT64HAk79AzXjIFcCzqvqR1Ey81TIN7dvoA46qWghMrFf2PvC+FzEYY9qgl7Jf4rGMx0iLT2PWyFmWRAKEJy2SLBH5O+5EslZEugJVXpwzH+hXZzsO+MqL4xljOoBXdr7CzC0zuajfRe4JGINsAsZA4Uk6nwCcDeSp6hERicbdvdVSGcBpIpIA7AXGAzd6cTxjTDu3ctdKHvnwEUbFjeKPP/qjzeIbYDxZIbFKRE4GbhYRBTaq6l88ObiIrABGATEikg9MV9VnRORO3JM/uoBlqvppC+M3xrRzr3/2eu16In8cZUkkEJ0wkYjIU8CpwIrqol+LSJqq3nGifVX1hkbK36KRQXJjjKnx19y/8uCmB0mNTWXuqLmEuEKcDsk0wJOurR8BZ9bM+CsizwHb/RqVMabDW/35an73we8476TzmHfRPEJdoU6HZBrhyWB7DhBfZ7sf8LF/wjHGGFiTt4ZpG6eRclIKT1z8hCWRANdoi0RE3sB9W243IFtEtlRvnwtsap3wjDEdzdtfvM3UjVMZ1mcY8y+eT1inMKdDMifQVNfWnCb+ny1sZYzxubX/WcuUDVM4p9c5PHnxk4R3Cnc6JOOBRhOJqq5vqFxEUnHfrpvur6CMMR3PP3b/g8npkxncczBPXfIUnYM7Ox2S8ZBHj4WKyNm4k8c44AvgNX8GZYzpWN778j3uX38/g2IG8VSaJZG2pqkxktNxPyx4A1AIvAKIql7USrEZYzqA9/e8z33r7yM5OpmFaQuJCI5wOiTTTE21SHYCG4CxqpoLICL3tEpUxpgOIT0/nXvev4eBPQayaPQiuoR0cTok0wJN3f57DfANsE5ElorIJXi4oJUxxpzIxr0buXvd3Zze43QWX7qYriFdnQ7JtFCjiURV/6Kq1wMDcc+0ew/QW0QWeroeiTHGNGTT3k3c9c+7OLX7qSwZvYTIkEinQzJeOOEDiar6naq+qKpX4p6pdxvQ5PK4xhhTX8GRAm555xbe+eIdfrvutyR0S2DJ6CV0C+3mdGjGS82azF9Vi4DF1T/GGOOxRR8vYuu+rWzbt40BPQaw9NKldA/r7nRYxgdsVRhjjF8NfWEoZZVltduVVPJZ8WekrUoj6+YsByMzvuLJXFvGGNNi71z9DmnxabXboa5QxiSMYe01ax2MyviSJRJjjF/FhMeQXZQNQHBQMGWVZUSERBATHuNwZMZXAj6RiMgAEXlGRFbVKRslIhtEZJGIjHIwPGPMCbz22WvsPbyXwTGDWTFmBeMSx1F4tNDpsIwP+XWMRESWAVcC+1X1zDrllwHzcK+Q+LSqzmrsGKqaB0yom0hwTxp5GAjDvQa8MSYAfVHyBbMzZnPuSeeyZPQSgiSIaedNczos42P+HmxfDjwJPF9TICIuYAEwGncSyBCR1biTysx6+9+qqvsbOO4GVV0vIr2BucBNfojdGOOF8spyJqdPJtQVyozUGQRJwHeAmBbyayJR1fTq9d7rSgFyq1saiMjLwFWqOhN368WT41ZV/1oM2Io3xgSg+f+eT3ZRNvMumkfviN5Oh2P8yImvCLHAnjrb+dVlDRKRaBFZBJwjIlOqy64WkcXAn3G3eBra73YRyRSRzIKCAt9Fb4w5oc1fbebZT5/lutOv4+L4i50Ox/iZE8+RNDRfV6MLZalqITCxXtnrwOtNnURVlwBLAIYNG2YLcRnTSopLi5m6cSoJ3RKYNHyS0+GYVuBEIsnHve57jTjgKwfiMMb4mKoyfdN0vj32LQsuWWArHHYQTnRtZQCniUiCiITgXvNktQNxGGN8bOWulazbs467htxFUnSS0+GYVuLXRCIiK4DNQKKI5IvIBFWtAO4E1gLZwKuq+qk/4zDG+F/et3n8IeMPXND3An6e/HOnwzGtyN93bd3QSPlbwFv+PLcxpvWUVZZxf/r9hHcK55HUR+xW3w7GJm00xnht3tZ55BTnMP/i+fTs3NPpcEwrs68NxhivbNq7ied3PM/1idczqt8op8MxDrBEYoxpsaLSIqZ+MJVTup3CfcPuczoc4xDr2jLGtIiq8rsPfkfJsRIWpS0irFOY0yEZh1iLxBjTIq/kvML6/PXcO/ReEqMSnQ7HOMgSiTGm2XKLc5mTOYfU2FRuSrI5Uzs6SyTGmGY5VnmM+zfcT0RwBI+kPoJIQ7MemY7ExkiMMc3yeNbjfFb8GQsuWWCrHBrAWiTGmGbYuHcjL2S/wI0Db2Rk3EinwzEBwhKJ8ZuCIwXc8s4tHDh6oF2cp6M7cPQAUzdO5dTup3LvsHudDscEEEskxm8WfbyIrfu2svCjhe3iPB1Zza2+h8sOM3vkbEJdtp6c+Z6otv+lOoYNG6aZmZlOh9FhDH1hKGWVZceVd5JOTE6Z7LPzPLblMSq04rjyEFcIWTdn+ew8Bl7MfpFZW2bxQMoDdpdWByIiWao67ESvs8F24xOFRwvJ3JdJxjcZ9O7cmz2H9hz3mgqtYMaHM/weS9+IvkzZMIWBUQNJjk4mMSqRyJBIv5+3vdpVvIu5mXMZETuCGwfe6HQ4JgBZIjEtUlxaXJs4Mr7JIPfbXAA6d+rMkN5D6BHag+0HthMcFEx5VTljTxnLvUN9368+N3Mub+S9UXuegVED6dm5Jx9+/SFr8tbUvi62SyxJUUkMjBpIUnQSSVFJNrmgB0orSpmcPpmuIV15OPVhu9XXNMgSifFIybGSHySOXcW7AAjvFM45vc5hzIAxpPRJISk6ieCgYO5edzfjEsdx3enXsXLXSg4cPUB0eLTP4/qu4rvjzvP4RY8D7sHhnUU72Vm0k+zCbLKLsnn3y3dr940Oi2Zg9MDaBJMclUxc17gmL5YFRwqYlD6JOT+a0yFufZ2bNZfcb3NZmLbQL38/0z7YGIlp0MGyg2R9k0XGPnfiyCnKQVHCXGGc3etshvcZTkqfFM6IOYPgoGCnw/XYobJD5BTluJNLkTu55H2bR6VWAtAluAuJUYkkRSWRFO1OMAndEmrr+PC/HmZlzkquS7yO/zvv/5ysit+l56dzx3t3cHPSzT4d2zJth6djJAGfSERkADAV6Kaq11aXjQBuwt2iSlbVC5o6hiWS7zX2jfpw2WG27t/Klq+3kLEvg51FO6nSKkKCQn6QOM6MOZMQV4iDNfC9Y5XHyC3OJbsou7b1sqt4F6WVpQCEBIVQXlWOcvxnpb0O7B84eoBrVl9DTHgML415ye7S6qACYrBdRJYBVwL7VfXMOuWXAfMAF/C0qs5q7BiqmgdMEJFVdco2ABtE5Ke414A3Hqq5VXb+1vmMPnk0W77ZQsbXGewo2kGVVhEcFMzgnoP59Vm/Znif4ZzV86x2fxEJdYVyRswZnBFzRm1ZRVUFuw/udieXwp18XPAx2w9s/8FdYvFd45k4eCJHK44S3incidD9okqrmPbBNL4r/45lP17W7v/+xnt+bZGIyEjgMPB8TSIRERewCxgN5ONOBDfgTioz6x3iVlXdX73fqpoWSZ3jvwrcpqoHm4qjLbRIfNX3rqocKj9E0dEiikqLKC4tprC0kEc/fLS2+6a+Ib2GMLzPcIb3Gc7gnoNtOvBGPLT5IVbtWoVLXFRoBS5xUamVhLpCGd5nOCNiRzAibgT9uvZzOlSv/HnHn5mdMZup505l/MDxTodjHBQQLRJVTReRk+sVpwC51S0NRORl4CpVnYm79eIREYkHShpLIiJyO3A7QHx8fPODb2V1H6qr3/d+pPxIbVIoKnUniMLSwh9sF5UWuZPHsSIqqo5/tgLcz3FUaiWK4hIXQ3sP5cHzH6RfZNu+8LWWotKiHwzs7z+yn/GJ49mwdwPp+enM3DuTmVtmktAtgRGxIxgZN5IhvYYQ7Go7Y0g5RTn8KetPjIobxfWJ1zsdjmkj/D5GUp1I1tRpkVwLXKaqt1Vv/xw4V1XvbGT/aGAG7hbM09UJBxH5PbBWVTedKIZAbZFUVFWQ8mIK5VXlx/0/QejbpS9FpUUcrTja4P7hncKJCotq8KdHWA+iw6KJCq/eDu3BzC0zWbVrFcGuYMoryzvEgHFr2n1wNxvy3Uklc18m5VXlRARHcP5J5zMibgQjYkcE9C3HpRWljF8znpKyEl77yWtEhUU5HZJxWEC0SBrR0L2VjWYzVS0EJjZQPt2XQflLlVbxzXffsPvgbr48+CW7D1X/9+Bu8g/nN9h6iAyJJLFHIn0i+tAjrEdtcogOj65NEj1Ce9A5uHOzYqn/jdrmpvKt/pH96Z/cn5uTb+ZI+RH+9fW/2LB3AxvyN9TedpwUlVSbVAbFDMIV5HI46u/NyZzD5yWfs3j0YksiplmcaJGcDzyoqj+u3p4CUNPS8IeWtkg8HbdQVQqOFnyfLA7udv9+6Ev2HNrDscpjta8Nc4XRL7If/bv2Jz4ynpMjT2bdnnW8v+d9aym0U6rKruJdtUllW8E2qrSK7qHdSY1NZWTsSFJjU+kW2q12n9Z8XqXgSAG/+vuv+Lzkc36R/AsmDZ/k1/OZtiOQWyQZwGkikgDsBcYDATnvQt1xi2nnTqOotIgvD33ZYMKo2/0UHBRMv679iI+MJ7Vvam3CiI+Mp1fnXgTJD+fKXJ+/3loK7ZiIkBiVSGJUIrcNuo2SYyVs/moz6fnpbNy7kTfz3iRIgjgr5ixGxo1kRNwIVuasbHTMzNcez3qcz0s+p3tod+4acpdfz2XaJ3/ftbUCGAXEAPuA6ar6jIhcATyO+06tZarq1wmYmtsiaWzSwbpc4iK2S6y7OyPS3bqoaWWcFHFSQHVZmMBVWVXJp4Wf1g7Y7yjc0eDrgiSIMQljfHruN794kyqtOq68vT4bY5qv3TyQ6AvNTSQFRwqYkzmHd3e/S1lVGUEEEdc1jrGnjCU5Opn+kf3p26Vvm3qi27QNOUU5PLT5IT4p/KT2Ih/eKZxuId18/uWksqqSkrISSitKa2ctuCT+Eu4bfl+HmP7FnFggd20FvJ6dexIRHEF5VTkhrhDKK8s5r+95TBx83Ji/MT5V0wW2/cD22vfe2FPG+q17q+bZmBBXCMcqjxEREmFJxDSbJZJG2B1Oximt+d6z97nxBevaMsYY0yBPu7ZsqV1jjDFesURijDHGK5ZIjDHGeMUSiTHGGK9YIjHGGOMVSyTGGGO80iFu/xWRAmC303F4IAZobzfyW53aBqtT29Dadeqvqidc+6BDJJK2QkQyPblnuy2xOrUNVqe2IVDrZF1bxhhjvGKJxBhjjFcskQSWJU4H4AdWp7bB6tQ2BGSdbIzEGGOMV6xFYowxxiuWSIwxxnjFEokxxhivWCIJYCIyQESeEZFVdcoiROQ5EVkqIjc5GZ+3RCReRFaLyDIRecDpeLwlIkEiMkNE5ovIfzkdj69Uv+eyRORKp2PxBRH5afXn528icqnT8bRUIF0LLJH4SfXFcb+IfFKv/DIRyRGR3BNdPFU1T1Un1Cu+Glilqr8CfuLjsD3mi/oBpwNvquqtQLLfgvWAj+pzFRALlAP5/orVUz6qE8Bk4FX/RNk8Pvpc/bX683MLcL0fw222ZtYvIK4FYInEn5YDl9UtEBEXsAC4HPeF8wYRSRaRQSKypt5Pr0aOGwfsqf690k+xe2I53tfv38B4EfknsK6V469vOd7XJxHYrKr3Ar9p5fgbshwv6yQiacAOYF9rB9+I5fjuczWter9AshwP60fgXAtszXZ/UdV0ETm5XnEKkKuqeQAi8jJwlarOBDztNsjH/QbahoNfBHxRPxG5D5hefaxVwLP+jbpxPqpPPlBWvenoBxt8VqeLgAjcF7CjIvKWqlb5NfAm+KhOAswC3lbVrf6NuHmaUz8C5FqA0yfvgGL5/hsEuN8IsY29WESiRWQRcI6ITKkufh24RkQWAm/4LdKWaVb9gHeA31bX8T9+jKulmluf14Efi8h8IN2fgXmhWXVS1amqejfwErDUySTShOb+nf4HSAOuFZGJ/gzMRxqrX8BcC6xF0rqkgbJGnwhV1UJgYr2y74Bf+jguX2lu/T4BrvVfOF5rbn2OAPXHtAJNs+pU+wLV5b4PxWea+3d6AnjCf+H4XIP1C6RrgbVIWlc+0K/OdhzwlUOx+EN7q197qw9YndqigK+fJZLWlQGcJiIJIhICjAdWOxyTL7W3+rW3+oDVqS0K+PpZIvETEVkBbAYSRSRfRCaoagVwJ7AWyAZeVdVPnYyzpdpb/dpbfcDq5GScLdVW62eTNhpjjPGKtUiMMcZ4xRKJMcYYr1giMcYY4xVLJMYYY7xiicQYY4xXLJEYY4zxiiUSY4wxXrFEYowxxiuWSIxxiIg8KyJXikh3EXlbRH7mdEzGtIQlEmOcMwj4Fvgb8LCq/sXheIxpEZsixRgHiEgQcAgoBBao6mMOh2RMi1mLxBhnnIZ7KvBbgIkiEuxsOMa0nCUSY5wxCPiHqv4T+AT4hcPxGNNilkiMccYg3AkE4FFgiojYiqWmTbIxEmOMMV6xFokxxhivWCIxxhjjFUskxhhjvGKJxBhjjFcskRhjjPGKJRJjjDFesURijDHGK5ZIjDHGeOX/Ab1IrozjmC+JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "deltas = 10.**np.array([1.5,1,0.5,0,-0.5,-1,-2,-3,-4,-5,-6,-7,-8,-9,-10])\n",
    "\n",
    "# store forward difference results\n",
    "sens_FD = np.zeros(deltas.size)\n",
    "\n",
    "# store central difference results\n",
    "sens_CD = np.zeros(deltas.size)\n",
    "\n",
    "# store complex step approximation results\n",
    "sens_CS = np.zeros(deltas.size)\n",
    "\n",
    "# Q calculated from the mean value of parameters\n",
    "Q_nom = Q\n",
    "\n",
    "# kappa value in the x scale \n",
    "mean_kappa = kappa_func(xs, kappal_nom, kappah_nom)\n",
    "\n",
    "# Calculate the sensitivity at each delta\n",
    "for count, delta in enumerate(deltas):\n",
    "    i = 1250\n",
    "    \n",
    "    # Get sensitivity from forward difference \n",
    "    pert = np.ones(Nx)\n",
    "    pert[i] += (delta)\n",
    "    sol, Qnew = ADRSource(Lx, Nx, Source_func(xs, q_nom), omega_func(xs, omega_nom), v_func(xs, v_nom),\n",
    "                              mean_kappa*pert)\n",
    "    sens_FD[count] = (Qnew - Q_nom)/(delta*mean_kappa[i])\n",
    "    \n",
    "    # Get sensitivity from central difference\n",
    "    pert = np.ones(Nx)\n",
    "    pert[i] += (delta/2)\n",
    "    sol, Qfor = ADRSource(Lx, Nx, Source_func(xs, q_nom), omega_func(xs, omega_nom), v_func(xs, v_nom),\n",
    "                              mean_kappa*pert)\n",
    "    \n",
    "    pert = np.ones(Nx)\n",
    "    pert[i] += (-delta/2)\n",
    "    sol, Qback = ADRSource(Lx, Nx, Source_func(xs, q_nom), omega_func(xs, omega_nom), v_func(xs, v_nom),\n",
    "                              mean_kappa*pert)\n",
    "    sens_CD[count] = (Qfor - Qback)/(delta*mean_kappa[i])\n",
    "    \n",
    "    # Get sensitivity from the complex step approximation\n",
    "    pert = np.ones(Nx,dtype=\"complex\")\n",
    "    pert[i] += (delta*1j)\n",
    "    sol, Qnew = ADRSource(Lx, Nx, Source_func(xs, q_nom), omega_func(xs, omega_nom), v_func(xs, v_nom),\n",
    "                          mean_kappa*pert)\n",
    "    sens_CS[count] = np.imag(Qnew)/(delta*mean_kappa[i])\n",
    "\n",
    "\n",
    "# Calculate the sensitivity with complex step, delta=1E-12 as the exact sensitivity\n",
    "delta_c = 1e-12\n",
    "pert = np.ones(Nx,dtype=\"complex\")\n",
    "pert[i] += (delta_c*1j)\n",
    "sol, Qnew = ADRSource(Lx, Nx, Source_func(xs, q_nom), omega_func(xs, omega_nom), v_func(xs, v_nom),\n",
    "                      mean_kappa*pert)\n",
    "exact = np.imag(Qnew)/(delta_c*mean_kappa[i])\n",
    "\n",
    "# Sensitivity plot \n",
    "plt.semilogx(deltas*mean_kappa[i], sens_FD,'o-', label='Forward difference')\n",
    "plt.plot(deltas*mean_kappa[i], sens_CD,'^-', label='Center difference')\n",
    "plt.plot(deltas*mean_kappa[i], sens_CS,'*-', label='Complex step')\n",
    "plt.ylabel('Sensitivity')\n",
    "plt.xlabel('$\\kappa$')\n",
    "plt.title('Sensitivity of Q to $\\kappa$')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# Absolute error plot\n",
    "plt.loglog(deltas*mean_kappa[i], np.abs(exact-sens_FD),'o-', label='Forward difference')\n",
    "plt.loglog(deltas*mean_kappa[i], np.abs(exact-sens_CD),'^-', label='Center difference')\n",
    "plt.loglog(deltas*mean_kappa[i], np.abs(exact-sens_CS),\"*-\", label='Complex step')\n",
    "plt.ylabel('Abs. error of sensitivity')\n",
    "plt.xlabel('$\\kappa$')\n",
    "plt.title('Abs. error of sensitivity of Q to $\\kappa$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.3 Second-Derivative approximations](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.3-Second-Derivative-approximations)",
     "section": "4.1.3 Second-Derivative approximations"
    }
   },
   "source": [
    "## 4.1.3 Second-Derivative approximations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.3 Second-Derivative approximations](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.3-Second-Derivative-approximations)",
     "section": "4.1.3 Second-Derivative approximations"
    }
   },
   "source": [
    "The second derivative: \n",
    "\n",
    "$\\frac{\\partial^2 Q}{\\partial x_i^2} \\bigg|_{\\overline{x}} \\approx \\frac{Q(\\overline{x} + \\delta_i\\hat{e_i}) - 2Q(\\overline{x}) + Q(\\overline{x} - \\delta_i\\hat{e_i})}{\\delta_i^2}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.1.3 Second-Derivative approximations](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html#4.1.3-Second-Derivative-approximations)",
     "section": "4.1.3 Second-Derivative approximations"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The second derivative for the five parameters are [ 3.55503929e+00  9.36436493e+00  7.08268999e-02 -2.06039637e+01\n",
      " -7.10542736e-06]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:47: ComplexWarning: Casting complex values to real discards the imaginary part\n"
     ]
    }
   ],
   "source": [
    "# Get Q value with average value of parameters\n",
    "xs = np.linspace(dx/2,Lx-dx/2,Nx)\n",
    "source = Source_func(xs, 1)\n",
    "kappa = kappa_func(xs, 0.1, 2)\n",
    "\n",
    "omega_nom = 20\n",
    "omega_var = 0.3195214\n",
    "v_nom = 10\n",
    "v_var = 0.0723493\n",
    "kappal_nom = 0.1\n",
    "kappal_var = 8.511570e-6\n",
    "kappah_nom = 2\n",
    "kappah_var = 0.002778142\n",
    "q_nom = 1\n",
    "q_var = 7.062353e-4\n",
    "vs = v_func(xs, v_nom)\n",
    "\n",
    "# Q(x) needs to be evaluated only once \n",
    "sol,q_center = ADRSource(Lx, Nx, Source_func(xs, q_nom), omega_func(xs, omega_nom), \n",
    "                                 vs,                       \n",
    "                                 kappa_func(xs, kappal_nom, kappah_nom))  \n",
    "\n",
    "# define perturbation\n",
    "delta = 1e-4\n",
    "\n",
    "sens_q2 = np.zeros(5)\n",
    "\n",
    "# loop over each parameter \n",
    "for i in range(5):\n",
    "    deltas = np.zeros(5)\n",
    "    # only one parameter is perturbed\n",
    "    deltas[i] = 1.0*delta\n",
    "        \n",
    "    # get the forward term \n",
    "    sol,forward_pert = ADRSource(Lx, Nx, Source_func(xs, q_nom*(1+deltas[4])), omega_func(xs, omega_nom*(1+deltas[1])), \n",
    "                                 vs*(1+deltas[0]),                       \n",
    "                                 kappa_func(xs, kappal_nom*(1+deltas[2]), kappah_nom*(1+deltas[3])))  \n",
    "    \n",
    "    # get the backward term\n",
    "    sol,backward_pert = ADRSource(Lx, Nx, Source_func(xs, q_nom*(1-deltas[4])), omega_func(xs, omega_nom*(1-deltas[1])), \n",
    "                                 vs*(1-deltas[0]),                       \n",
    "                                 kappa_func(xs, kappal_nom*(1-deltas[2]), kappah_nom*(1-deltas[3])))\n",
    "    \n",
    "    noms = np.array((v_nom,omega_nom,kappal_nom, kappah_nom, q_nom))\n",
    "    \n",
    "    # the second derivative for the perturbed parameter \n",
    "    sens_q2[i] = (forward_pert - 2*q_center + backward_pert)/(np.sum(noms*deltas))**2\n",
    "    \n",
    "print('The second derivative for the five parameters are', sens_q2*noms*noms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [4.0 Local Sensitivity Analysis Based on Derivative Approximations](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.00-Local-Sensitivity-Analysis-Based-on-Derivative-Approximations.html) | [Contents](toc.html) | [4.2 Advection-Diffusion-Reaction (ADR) Example](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html)<p><a href=\"https://colab.research.google.com/github/ndcbe/cbe67701-uncertainty-quantification/blob/master/docs/04.01-Local-sensitivity-analysis.ipynb\"> <img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.ipynb\"> <img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}