{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adaptive PDE discretizations on Cartesian grids\n",
    "## Volume : Algorithmic tools\n",
    "## Part : Domain representation\n",
    "## Chapter : Finite differences, interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook presents the basic types of finite differences and interpolation methods that can be considered on a cartesian grid. The tools are presented in two dimensions, but apply in arbitrary dimension. They also apply in the context of:\n",
    "* non-square domains, [example](SubsetRd.ipynb)\n",
    "* point dependent offsets, [example](../Notebooks_NonDiv/LinearMonotoneSchemes2D.ipynb)\n",
    "* vector valued functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "[**Summary**](Summary.ipynb) of volume Algorithmic tools, this series of notebooks.\n",
    "\n",
    "[**Main summary**](../Summary.ipynb) of the Adaptive Grid Discretizations \n",
    "\tbook of notebooks, including the other volumes.\n",
    "\n",
    "# Table of contents\n",
    "  * [1. Degenerate elliptic finite differences](#1.-Degenerate-elliptic-finite-differences)\n",
    "    * [1.1 Testing data](#1.1-Testing-data)\n",
    "    * [1.2 Upwind finite difference](#1.2-Upwind-finite-difference)\n",
    "    * [1.3 Second order finite difference](#1.3-Second-order-finite-difference)\n",
    "    * [1.4 Centered finite difference](#1.4-Centered-finite-difference)\n",
    "  * [2. Higher order finite differences](#2.-Higher-order-finite-differences)\n",
    "    * [2.1 Upwind finite difference](#2.1-Upwind-finite-difference)\n",
    "    * [2.2 Second order finite difference](#2.2-Second-order-finite-difference)\n",
    "    * [2.3 Centered finite difference](#2.3-Centered-finite-difference)\n",
    "  * [3. Composite finite differences](#3.-Composite-finite-differences)\n",
    "    * [3.1 Gradient](#3.1-Gradient)\n",
    "    * [3.2 Hessian](#3.2-Hessian)\n",
    "  * [4. Interpolation](#4.-Interpolation)\n",
    "    * [4.1 Linear splines](#4.1-Linear-splines)\n",
    "    * [4.2 Quadratic splines](#4.2-Quadratic-splines)\n",
    "    * [4.3 Cubic splines](#4.3-Cubic-splines)\n",
    "    * [4.4 Vector data](#4.4-Vector-data)\n",
    "    * [4.5 Comparison with scipy ndimage map_coordinates](#4.5-Comparison-with-scipy-ndimage-map_coordinates)\n",
    "  * [5. Functions associated with an AD variable](#5.-Functions-associated-with-an-AD-variable)\n",
    "    * [5.1 Taylor expansions](#5.1-Taylor-expansions)\n",
    "    * [5.2 Sparse differentiation](#5.2-Sparse-differentiation)\n",
    "\n",
    "\n",
    "\n",
    "**Acknowledgement.** Some of the experiments presented in these notebooks are part of \n",
    "ongoing research with Ludovic Métivier and Da Chen.\n",
    "\n",
    "Copyright Jean-Marie Mirebeau, Centre Borelli, ENS Paris-Saclay, CNRS, University Paris-Saclay"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Import the required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.347477Z",
     "iopub.status.busy": "2024-04-30T09:22:36.347217Z",
     "iopub.status.idle": "2024-04-30T09:22:36.357104Z",
     "shell.execute_reply": "2024-04-30T09:22:36.356628Z"
    }
   },
   "outputs": [],
   "source": [
    "import sys; sys.path.insert(0,\"..\") # Allow imports from parent directory\n",
    "#from Miscellaneous import TocTools; print(TocTools.displayTOC('FiniteDifferences','Algo'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.359779Z",
     "iopub.status.busy": "2024-04-30T09:22:36.359580Z",
     "iopub.status.idle": "2024-04-30T09:22:36.493755Z",
     "shell.execute_reply": "2024-04-30T09:22:36.493470Z"
    }
   },
   "outputs": [],
   "source": [
    "from agd import FiniteDifferences as fd\n",
    "from agd import AutomaticDifferentiation as ad\n",
    "from agd import Interpolation\n",
    "from agd import LinearParallel as lp\n",
    "norm_infinity = ad.Optimization.norm_infinity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.495569Z",
     "iopub.status.busy": "2024-04-30T09:22:36.495420Z",
     "iopub.status.idle": "2024-04-30T09:22:36.679785Z",
     "shell.execute_reply": "2024-04-30T09:22:36.679508Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Degenerate elliptic finite differences\n",
    "\n",
    "The finite difference presented in this section are a typical ingredient of monotone numerical schemes, see the corresponding [volume](../Notebooks_NonDiv/Summary.ipynb). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Testing data\n",
    "\n",
    "In order to test the finite difference and interpolation methods, we need some polynomial functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.681489Z",
     "iopub.status.busy": "2024-04-30T09:22:36.681347Z",
     "iopub.status.idle": "2024-04-30T09:22:36.683747Z",
     "shell.execute_reply": "2024-04-30T09:22:36.683532Z"
    }
   },
   "outputs": [],
   "source": [
    "def u1(X): return X[0]+2*X[1]\n",
    "def u2(X): return X[0]**2+2*(2*X[0]*X[1])+3*X[1]**2\n",
    "def u3(X): return X[0]**3+X[0]*X[1]**2\n",
    "\n",
    "def u123(X): return ad.array( (u1(X),u2(X),u3(X)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to choose a direction, with *integer coordinates*, for the finite differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.685132Z",
     "iopub.status.busy": "2024-04-30T09:22:36.685045Z",
     "iopub.status.idle": "2024-04-30T09:22:36.686670Z",
     "shell.execute_reply": "2024-04-30T09:22:36.686469Z"
    }
   },
   "outputs": [],
   "source": [
    "e = (1,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us also define a domain, here a square."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.688013Z",
     "iopub.status.busy": "2024-04-30T09:22:36.687928Z",
     "iopub.status.idle": "2024-04-30T09:22:36.691026Z",
     "shell.execute_reply": "2024-04-30T09:22:36.690795Z"
    }
   },
   "outputs": [],
   "source": [
    "aX,h = np.linspace(-1,1,retstep=True)\n",
    "X=np.array(np.meshgrid(aX,aX,indexing='ij'))\n",
    "shape = X.shape[1:]\n",
    "\n",
    "def interior(shape,k):\n",
    "    \"\"\"Boolean array excluding k boundary layers\"\"\"\n",
    "    interior = np.full(shape,True)\n",
    "    for i in range(len(shape)): \n",
    "        #interior[*(slice(None),)*i,:k]=False\n",
    "        #interior[*(slice(None),)*i,-k:]=False\n",
    "        interior.__setitem__((slice(None),)*i+(slice(None,k),),False)\n",
    "        interior.__setitem__((slice(None),)*i+(slice(-k,None),),False)\n",
    "    return interior\n",
    "    \n",
    "def close(u,v,k,ndim=2,**kwargs):\n",
    "    \"\"\"\n",
    "    Wether u and v are close in the domain minus k boundary layers.\n",
    "    - **kwargs : passed to np.allclose\n",
    "    \"\"\"\n",
    "    dom = interior(u.shape[-ndim:],k)\n",
    "#    print(norm_infinity(u[...,dom]-v[...,dom]),kwargs)\n",
    "    return np.allclose(u[...,dom],v[...,dom],**kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.692374Z",
     "iopub.status.busy": "2024-04-30T09:22:36.692289Z",
     "iopub.status.idle": "2024-04-30T09:22:36.695361Z",
     "shell.execute_reply": "2024-04-30T09:22:36.695102Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False, False, False, False, False],\n",
       "       [False,  True,  True,  True, False],\n",
       "       [False,  True,  True,  True, False],\n",
       "       [False,  True,  True,  True, False],\n",
       "       [False, False, False, False, False]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "interior((5,5),1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following variables are used for validation, by comparison with automatic differentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.714363Z",
     "iopub.status.busy": "2024-04-30T09:22:36.714238Z",
     "iopub.status.idle": "2024-04-30T09:22:36.718854Z",
     "shell.execute_reply": "2024-04-30T09:22:36.718600Z"
    }
   },
   "outputs": [],
   "source": [
    "X_ad = ad.Dense.identity(constant=X,shape_free=(2,))\n",
    "X_ad2 = ad.Dense2.identity(constant=X,shape_free=(2,))\n",
    "\n",
    "du1  = u1(X_ad).gradient()\n",
    "du2  = u2(X_ad).gradient()\n",
    "ddu2 = u2(X_ad2).hessian()\n",
    "du3  = u3(X_ad).gradient()\n",
    "ddu3 = u3(X_ad2).hessian()\n",
    "\n",
    "_e = fd.as_field(e,shape)\n",
    "du1_e  = lp.dot_VV(du1,_e)\n",
    "du2_e  = lp.dot_VV(du2,_e)\n",
    "ddu2_e = lp.dot_VAV(_e,ddu2,_e)\n",
    "du3_e  = lp.dot_VV(du3,_e)\n",
    "ddu3_e = lp.dot_VAV(_e,ddu3,_e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Upwind finite difference\n",
    "$$\n",
    "    \\frac{u(x+he)-u(x)} h = <\\nabla u(x),e> + O(h).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.720249Z",
     "iopub.status.busy": "2024-04-30T09:22:36.720165Z",
     "iopub.status.idle": "2024-04-30T09:22:36.722334Z",
     "shell.execute_reply": "2024-04-30T09:22:36.722113Z"
    }
   },
   "outputs": [],
   "source": [
    "Du1_e = fd.DiffUpwind(u1(X),e,h)\n",
    "assert close(Du1_e,du1_e,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Second order finite difference\n",
    "$$\n",
    "    \\frac{u(x+he)-2u(x)+u(x-he)}{h^2} = <e,\\nabla^2 u(x),e> + O(h^2).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.723697Z",
     "iopub.status.busy": "2024-04-30T09:22:36.723612Z",
     "iopub.status.idle": "2024-04-30T09:22:36.725709Z",
     "shell.execute_reply": "2024-04-30T09:22:36.725502Z"
    }
   },
   "outputs": [],
   "source": [
    "DDu2_e = fd.Diff2(u2(X),e,h)\n",
    "assert close(DDu2_e,ddu2_e,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Centered finite difference\n",
    "\n",
    "Centered finite differences *are not* degenerate elliptic, but they can be used within degenerate elliptic schemes if they are suitably compensated by second order finite differences.\n",
    "$$\n",
    "    \\frac{u(x+h e)-u(x-h e)} {2 h} = <\\nabla u(x),e> + O(h^2). \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.727008Z",
     "iopub.status.busy": "2024-04-30T09:22:36.726931Z",
     "iopub.status.idle": "2024-04-30T09:22:36.729020Z",
     "shell.execute_reply": "2024-04-30T09:22:36.728790Z"
    }
   },
   "outputs": [],
   "source": [
    "Du2_e = fd.DiffCentered(u2(X),e,h)\n",
    "assert close(Du2_e,du2_e,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Higher order finite differences\n",
    "\n",
    "High order finite differences are not degenerate elliptic. They may be used within filtered schemes, which combine a stable degenerate elliptic scheme with a higher order scheme.\n",
    "By construction, they are exact on higher order polynomials."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Upwind finite difference\n",
    "$$\n",
    "    \\frac{-u(x+2he)+4u(x+he)-3u(x)}{2h} = <\\nabla u(x),e> + O(h^2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.730346Z",
     "iopub.status.busy": "2024-04-30T09:22:36.730272Z",
     "iopub.status.idle": "2024-04-30T09:22:36.732382Z",
     "shell.execute_reply": "2024-04-30T09:22:36.732164Z"
    }
   },
   "outputs": [],
   "source": [
    "Du2_e = fd.DiffUpwind(u2(X),e,h,order=2)\n",
    "assert close(Du2_e,du2_e,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "    \\frac{2 u(x+3he)-9 u(x+2he)+18 u(x+he)-11 u(x)}{6h} = <\\nabla u(x),e> + O(h^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.733740Z",
     "iopub.status.busy": "2024-04-30T09:22:36.733660Z",
     "iopub.status.idle": "2024-04-30T09:22:36.735836Z",
     "shell.execute_reply": "2024-04-30T09:22:36.735610Z"
    }
   },
   "outputs": [],
   "source": [
    "Du3_e = fd.DiffUpwind(u3(X),e,h,order=3)\n",
    "assert close(Du3_e,du3_e,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Second order finite difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.737249Z",
     "iopub.status.busy": "2024-04-30T09:22:36.737172Z",
     "iopub.status.idle": "2024-04-30T09:22:36.739553Z",
     "shell.execute_reply": "2024-04-30T09:22:36.739340Z"
    }
   },
   "outputs": [],
   "source": [
    "DDu3_e = fd.Diff2(u3(X),e,h,order=4)\n",
    "assert close(DDu3_e,ddu3_e,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Centered finite difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.740945Z",
     "iopub.status.busy": "2024-04-30T09:22:36.740870Z",
     "iopub.status.idle": "2024-04-30T09:22:36.743011Z",
     "shell.execute_reply": "2024-04-30T09:22:36.742792Z"
    }
   },
   "outputs": [],
   "source": [
    "Du3_e = fd.DiffCentered(u3(X),e,h,order=4)\n",
    "assert close(Du3_e,du3_e,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Composite finite differences\n",
    "The following finite differences can be used to estimate numerically the derivatives of a function, but they are rarely adequate for building numerical schemes. We denote by $e_i$ the $i$-th element of the canonical basis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Gradient\n",
    "$$\n",
    "    \\frac{u(x+h e_i)-u(x-he_i)}{2h} = \\frac {\\partial u} {\\partial x_i} + O(h^2), \\qquad 0 \\leq i < d.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.744406Z",
     "iopub.status.busy": "2024-04-30T09:22:36.744333Z",
     "iopub.status.idle": "2024-04-30T09:22:36.746762Z",
     "shell.execute_reply": "2024-04-30T09:22:36.746514Z"
    }
   },
   "outputs": [],
   "source": [
    "Du2 = fd.DiffGradient(u2(X),gridScale=h)\n",
    "assert close(Du2,du2,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.748137Z",
     "iopub.status.busy": "2024-04-30T09:22:36.748064Z",
     "iopub.status.idle": "2024-04-30T09:22:36.750138Z",
     "shell.execute_reply": "2024-04-30T09:22:36.749908Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 50, 50)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Du2.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Hessian\n",
    "$$\n",
    "    \\frac{u(x+h e_i)-2u(x)+u(x-he_i)}{h^2} = \\frac {\\partial^2 u} {\\partial^2 x_i} + O(h^2),\n",
    "$$\n",
    "for all $0\\leq i < d$, and \n",
    "$$\n",
    "    \\frac{u(x+h e_i+h e_j)+u(x-he_i-h e_j)-u(x+h e_i -h e_j) - u(x-h e_i+he_j)}{4h^2} = \\frac {\\partial^2 u} {\\partial x_i \\partial x_j} + O(h^2), \n",
    "$$\n",
    "for all distinct $i,j$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.751559Z",
     "iopub.status.busy": "2024-04-30T09:22:36.751479Z",
     "iopub.status.idle": "2024-04-30T09:22:36.757860Z",
     "shell.execute_reply": "2024-04-30T09:22:36.757615Z"
    }
   },
   "outputs": [],
   "source": [
    "DDu2 = fd.DiffHessian(u2(X),gridScale=h)\n",
    "assert close(DDu2,ddu2,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.759294Z",
     "iopub.status.busy": "2024-04-30T09:22:36.759191Z",
     "iopub.status.idle": "2024-04-30T09:22:36.761308Z",
     "shell.execute_reply": "2024-04-30T09:22:36.761090Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2, 50, 50)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DDu2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.762576Z",
     "iopub.status.busy": "2024-04-30T09:22:36.762479Z",
     "iopub.status.idle": "2024-04-30T09:22:36.765570Z",
     "shell.execute_reply": "2024-04-30T09:22:36.765362Z"
    }
   },
   "outputs": [],
   "source": [
    "DDu3 = fd.DiffHessian(u3(X),gridScale=h,order=4)\n",
    "assert close(DDu3,ddu3,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## 4. Interpolation\n",
    "\n",
    "The agd library contains a partial reimplementation of the `scipy.ndimage.map_coordinates` function, which allows AD types, both for the coordinates and the input values. This is another avenue for numerically differentiating a function defined on a grid.\n",
    "\n",
    "**Boundary conditions.**\n",
    "We only support reflected (default) and periodic boundary conditions (`grid-mirror`, `grid-wrap`). \n",
    "As a result, the spline interpolation methods *does not* exactly reproduce polynomials. In the tests below, we exclude some boundary layers, and choose a sufficiently high tolerance, to account for this non-exactness.\n",
    "\n",
    "Exact reproduction of low degree polynomials could be achieved with the `not a knot` boundary conditions, but these are substantially more difficult to implement, and numerically coslty. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define a finer grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.766863Z",
     "iopub.status.busy": "2024-04-30T09:22:36.766783Z",
     "iopub.status.idle": "2024-04-30T09:22:36.768663Z",
     "shell.execute_reply": "2024-04-30T09:22:36.768438Z"
    }
   },
   "outputs": [],
   "source": [
    "aX_ = np.linspace(-1,1,80)\n",
    "X_=np.array(np.meshgrid(aX_,aX_,indexing='ij'))\n",
    "shape_ = X_[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.769961Z",
     "iopub.status.busy": "2024-04-30T09:22:36.769886Z",
     "iopub.status.idle": "2024-04-30T09:22:36.771442Z",
     "shell.execute_reply": "2024-04-30T09:22:36.771191Z"
    }
   },
   "outputs": [],
   "source": [
    "Interp = Interpolation.UniformGridInterpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Linear splines\n",
    "\n",
    "(Piecewise) Linear splines are continuous, and reproduce linear functions.\n",
    "They are second order consistent\n",
    "$$\n",
    "    U_1^h(x) = u(x)+O(h^2).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.772746Z",
     "iopub.status.busy": "2024-04-30T09:22:36.772667Z",
     "iopub.status.idle": "2024-04-30T09:22:36.774895Z",
     "shell.execute_reply": "2024-04-30T09:22:36.774642Z"
    }
   },
   "outputs": [],
   "source": [
    "U1 = Interp(X,u1(X),order=1)\n",
    "assert np.allclose(U1(X_),u1(X_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spline can be differentiated, and yields the a first order consistent estimation of the gradient, *except possibly at boundary points*.\n",
    "$$\n",
    "    \\nabla U_1^h(x) = \\nabla u(x) + O(h)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.776275Z",
     "iopub.status.busy": "2024-04-30T09:22:36.776199Z",
     "iopub.status.idle": "2024-04-30T09:22:36.778673Z",
     "shell.execute_reply": "2024-04-30T09:22:36.778425Z"
    }
   },
   "outputs": [],
   "source": [
    "dU1 = U1(X_ad).gradient()\n",
    "assert close(dU1,du1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Quadratic splines\n",
    "\n",
    "Quadratic splines are continuously differentiable, and reproduce quadratic functions.\n",
    "We use a not-a-knot boundary condition : in one dimension, the second derivative is continuous accross the second node from the left.\n",
    "$$\n",
    "    U_2^h(x) = u(x)+O(h^3).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.779984Z",
     "iopub.status.busy": "2024-04-30T09:22:36.779906Z",
     "iopub.status.idle": "2024-04-30T09:22:36.782829Z",
     "shell.execute_reply": "2024-04-30T09:22:36.782610Z"
    }
   },
   "outputs": [],
   "source": [
    "U2 = Interp(X,u2(X),order=2)\n",
    "assert close(U2(X_),u2(X_),k=10,atol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spline can be differentiated, one time or two times, and yields a second order consistent estimate of the gradient, and a first order consistent estimate of the hessian, *except possibly at boundary points*.\n",
    "$$\n",
    "    \\nabla U_2^h(x) = \\nabla u(x) + O(h^2),\n",
    "$$\n",
    "$$\n",
    "    \\nabla^2 U_2^h(x) = \\nabla^2 u(x) + O(h).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.784277Z",
     "iopub.status.busy": "2024-04-30T09:22:36.784198Z",
     "iopub.status.idle": "2024-04-30T09:22:36.787079Z",
     "shell.execute_reply": "2024-04-30T09:22:36.786853Z"
    }
   },
   "outputs": [],
   "source": [
    "dU2 = U2(X_ad).gradient()\n",
    "assert close(dU2,du2,k=10,atol=1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.788533Z",
     "iopub.status.busy": "2024-04-30T09:22:36.788448Z",
     "iopub.status.idle": "2024-04-30T09:22:36.792107Z",
     "shell.execute_reply": "2024-04-30T09:22:36.791863Z"
    }
   },
   "outputs": [],
   "source": [
    "ddU2 = U2(X_ad2).hessian()\n",
    "assert close(ddU2,ddu2,k=10,atol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Cubic splines\n",
    "\n",
    "Cubic splines are twice continuously differentiable, and reproduce cubic functions.\n",
    "We use a not-a-knot boundary condition : in one dimension, the third derivative is continuous accross the second node from the left, and likewise from the right.\n",
    "$$\n",
    "    U_3^h(x) = u(x)+O(h^4).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.793659Z",
     "iopub.status.busy": "2024-04-30T09:22:36.793576Z",
     "iopub.status.idle": "2024-04-30T09:22:36.796966Z",
     "shell.execute_reply": "2024-04-30T09:22:36.796739Z"
    }
   },
   "outputs": [],
   "source": [
    "U3 = Interp(X,u3(X),order=3)\n",
    "assert close(U3(X_),u3(X_),k=10,atol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spline can be differentiated two times, and yields consistent estimates of the gradient and hessian, *except possibly at boundary points*.\n",
    "$$\n",
    "    \\nabla U_3^h(x) = \\nabla u(x) + O(h^3),\n",
    "$$\n",
    "$$\n",
    "    \\nabla^2 U_3^h(x) = \\nabla^2 u(x) + O(h^2).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.798511Z",
     "iopub.status.busy": "2024-04-30T09:22:36.798431Z",
     "iopub.status.idle": "2024-04-30T09:22:36.801525Z",
     "shell.execute_reply": "2024-04-30T09:22:36.801303Z"
    }
   },
   "outputs": [],
   "source": [
    "dU3 = U3(X_ad).gradient()\n",
    "assert close(dU3,du3,k=10,atol=1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.802932Z",
     "iopub.status.busy": "2024-04-30T09:22:36.802850Z",
     "iopub.status.idle": "2024-04-30T09:22:36.807913Z",
     "shell.execute_reply": "2024-04-30T09:22:36.807652Z"
    }
   },
   "outputs": [],
   "source": [
    "ddU3 = U3(X_ad2).hessian()\n",
    "assert close(ddU3,ddu3,k=12,atol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.4 Vector data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.809566Z",
     "iopub.status.busy": "2024-04-30T09:22:36.809483Z",
     "iopub.status.idle": "2024-04-30T09:22:36.815029Z",
     "shell.execute_reply": "2024-04-30T09:22:36.814795Z"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "U123 = Interp(X,u123(X),order=3)\n",
    "assert close(U123(X_),u123(X_),k=12,atol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### 4.5 Comparison with scipy ndimage map_coordinates\n",
    "\n",
    "As announced above, our spline interpolation function (partially) reproduces `scipy.ndimage.map_coordinate`, and extends it to AD types. We make this transparent by providing a map_coordinates function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.816441Z",
     "iopub.status.busy": "2024-04-30T09:22:36.816361Z",
     "iopub.status.idle": "2024-04-30T09:22:36.818082Z",
     "shell.execute_reply": "2024-04-30T09:22:36.817840Z"
    }
   },
   "outputs": [],
   "source": [
    "agd_map_coordinates = Interpolation.map_coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also provide a small wrapper over ndimage.map_coordinates, which allows AD types for the input values (not the coordinates), as well as the interpolation of tensor data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.819422Z",
     "iopub.status.busy": "2024-04-30T09:22:36.819347Z",
     "iopub.status.idle": "2024-04-30T09:22:36.820897Z",
     "shell.execute_reply": "2024-04-30T09:22:36.820655Z"
    }
   },
   "outputs": [],
   "source": [
    "ndimage_map_coordinates = Interpolation.ndimage_map_coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The implementations works both on the cpu and gpu, select by commenting/uncommenting the cell below. In the latter case, `scipy.ndimage.map_coordinates` is replaced with `cupyx.scipy.ndimage.map_coordinates`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.822317Z",
     "iopub.status.busy": "2024-04-30T09:22:36.822229Z",
     "iopub.status.idle": "2024-04-30T09:22:36.823785Z",
     "shell.execute_reply": "2024-04-30T09:22:36.823592Z"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "caster = lambda x:x # CPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.825077Z",
     "iopub.status.busy": "2024-04-30T09:22:36.824997Z",
     "iopub.status.idle": "2024-04-30T09:22:36.826562Z",
     "shell.execute_reply": "2024-04-30T09:22:36.826349Z"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "GPU_config"
    ]
   },
   "outputs": [],
   "source": [
    "#caster = ad.cupy_generic.cupy_set # GPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.827865Z",
     "iopub.status.busy": "2024-04-30T09:22:36.827776Z",
     "iopub.status.idle": "2024-04-30T09:22:36.868114Z",
     "shell.execute_reply": "2024-04-30T09:22:36.867856Z"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- order=1, dom_shape=(5,), x_shape=(8,), c_shape=(2, 3), size_ad=1 ---\n",
      "Absolute error :  denseAD(0.0,[0.])\n",
      "Reproduction error. impl: 0.0  scpy: 0.0\n",
      "Reproduction max error. impl: 0.0  scpy: 0.0\n",
      "--- order=2, dom_shape=(6,), x_shape=(4,), c_shape=(3,), size_ad=1 ---\n",
      "Absolute error :  denseAD(7.083111874806036e-12,[4.10549372e-12])\n",
      "Reproduction error. impl: 2.220446049250313e-16  scpy: 1.1591035076197187e-10\n",
      "Reproduction max error. impl: 2.220446049250313e-16  scpy: 6.718436917907411e-11\n",
      "--- order=3, dom_shape=(4,), x_shape=(2,), c_shape=(), size_ad=2 ---\n",
      "Absolute error :  denseAD(1.4526723157715082e-09,[-2.47407976e-07 -4.96368829e-07])\n",
      "Reproduction error. impl: 1.3877787807814457e-17  scpy: 2.307997604145129e-08\n",
      "Reproduction max error. impl: 2.220446049250313e-16  scpy: 7.886280012692204e-06\n",
      "--- order=4, dom_shape=(10,), x_shape=(3,), c_shape=(1,), size_ad=-1 ---\n",
      "Absolute error :  5.551115123125783e-17\n",
      "Reproduction error. impl: 3.3306690738754696e-16  scpy: 2.220446049250313e-16\n",
      "--- order=5, dom_shape=(10,), x_shape=(3, 2), c_shape=(1,), size_ad=-1 ---\n",
      "Absolute error :  6.661338147750939e-16\n",
      "Reproduction error. impl: 2.220446049250313e-16  scpy: 4.440892098500626e-16\n",
      "--- order=1, dom_shape=(4, 5), x_shape=(3,), c_shape=(), size_ad=-1 ---\n",
      "Absolute error :  1.1102230246251565e-16\n",
      "Reproduction error. impl: 0.0  scpy: 0.0\n",
      "--- order=2, dom_shape=(4, 6), x_shape=(3, 4), c_shape=(5,), size_ad=2 ---\n",
      "Absolute error :  denseAD(4.440892098500626e-16,[4.4408921e-16 1.2490009e-16])\n",
      "Reproduction error. impl: 3.3306690738754696e-16  scpy: 2.220446049250313e-16\n",
      "Reproduction max error. impl: 3.3306690738754696e-16  scpy: 3.3306690738754696e-16\n",
      "--- order=3, dom_shape=(4, 6), x_shape=(3, 4), c_shape=(5,), size_ad=2 ---\n",
      "Absolute error :  denseAD(8.41014820128494e-06,[3.36075072e-06 8.63728999e-06])\n",
      "Reproduction error. impl: 5.551115123125783e-16  scpy: 8.262379805445974e-06\n",
      "Reproduction max error. impl: 4.440892098500626e-16  scpy: 9.46352712638543e-06\n",
      "--- order=5, dom_shape=(4, 3), x_shape=(2, 1), c_shape=(2, 1), size_ad=1 ---\n",
      "Absolute error :  denseAD(0.0017610646499264804,[0.00317475])\n",
      "Reproduction error. impl: 3.3306690738754696e-16  scpy: 0.0031852717644261785\n",
      "Reproduction max error. impl: 3.3306690738754696e-16  scpy: 0.0042962930820822365\n",
      "--- order=1, dom_shape=(4, 5, 6), x_shape=(2, 5), c_shape=(2, 3), size_ad=1 ---\n",
      "Absolute error :  denseAD(2.220446049250313e-16,[1.11022302e-16])\n",
      "Reproduction error. impl: 0.0  scpy: 0.0\n",
      "Reproduction max error. impl: 0.0  scpy: 0.0\n",
      "--- order=3, dom_shape=(6, 2, 3), x_shape=(1, 4), c_shape=(), size_ad=2 ---\n",
      "Absolute error :  denseAD(0.0006403965303843351,[0.00030922 0.00044751])\n",
      "Reproduction error. impl: 5.551115123125783e-16  scpy: 0.000487351300054617\n",
      "Reproduction max error. impl: 5.551115123125783e-16  scpy: 0.0004580943068471788\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "for i,(order,dom_shape,x_shape,c_shape,size_ad,mode) in enumerate([\n",
    "    # One dimensional tests\n",
    "    (1,(5,),(8,),(2,3),1,'reflect'),\n",
    "    (2,(6,),(4,),(3,),1,'reflect'),\n",
    "    (3,(4,),(2,),tuple(),2,'grid-mirror'),\n",
    "    (4,(10,),(3,),(1,),-1,'grid-wrap'),\n",
    "    (5,(10,),(3,2),(1,),-1,'grid-wrap'),\n",
    "\n",
    "    # Two dimensional tests\n",
    "    (1,(4,5),(3,),tuple(),-1,'grid-wrap'),\n",
    "    (2,(4,6),(3,4),(5,),2,'grid-wrap'),\n",
    "    (3,(4,6),(3,4),(5,),2,'reflect'),\n",
    "    (5,(4,3),(2,1),(2,1),1,'reflect'),\n",
    "\n",
    "    # Three dimensional tests\n",
    "    (1,(4,5,6),(2,5),(2,3),1,'reflect'),\n",
    "    (3,(6,2,3),(1,4),tuple(),2,'grid-mirror'),\n",
    "]):\n",
    "    print(f\"--- {order=}, {dom_shape=}, {x_shape=}, {c_shape=}, {size_ad=} ---\")\n",
    "    vals = np.random.rand(*c_shape,*dom_shape)\n",
    "    if size_ad>0: vals=ad.Dense.denseAD(vals,np.random.rand(*vals.shape,size_ad))\n",
    "    pos = np.random.rand(len(dom_shape),*x_shape)*3*np.max(dom_shape)\n",
    "    vals,pos = map(caster,(vals,pos))\n",
    "    impl = agd_map_coordinates(vals,pos,order=order,mode=mode)\n",
    "    scpy = ndimage_map_coordinates(vals,pos,order=order,mode=mode) #scipy.ndimage\n",
    "    print(\"Absolute error : \",norm_infinity(impl-scpy))\n",
    "    assert np.allclose(impl,scpy,atol=1e-2)\n",
    "\n",
    "    # Note : the accuracy of scipy.ndimage uses float32 internally, hence its rather low accuracy \n",
    "    pos = np.array(np.meshgrid(*[np.arange(s,dtype=np.float64) for s in dom_shape], indexing='ij'))\n",
    "    pos = caster(pos)\n",
    "    impl = agd_map_coordinates(vals,pos,order=order,mode=mode)\n",
    "    scpy = ndimage_map_coordinates(vals,pos,order=order,mode=mode)\n",
    "    print(\"Reproduction error. impl:\",norm_infinity(ad.remove_ad(impl-vals)),\" scpy:\",norm_infinity(ad.remove_ad(scpy-vals)))\n",
    "    if ad.is_ad(impl) and impl.size_ad>0:\n",
    "        print(\"Reproduction max error. impl:\",norm_infinity((impl-vals).coef),\" scpy:\",norm_infinity((scpy-vals).coef))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Functions associated with an AD variable\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1 Taylor expansions\n",
    "\n",
    "The automatic differentiation classes store a Taylor expansion of the approximated function, which can be evaluated directly. The tangent, adjoint, and hessian operators may also be extracted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.869646Z",
     "iopub.status.busy": "2024-04-30T09:22:36.869540Z",
     "iopub.status.idle": "2024-04-30T09:22:36.871610Z",
     "shell.execute_reply": "2024-04-30T09:22:36.871370Z"
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([0.2,0.5])\n",
    "H = X-fd.as_field(x,shape)\n",
    "x_ad  = ad.Dense.identity(constant=x)\n",
    "x_ad2 = ad.Dense2.identity(constant=x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.872919Z",
     "iopub.status.busy": "2024-04-30T09:22:36.872842Z",
     "iopub.status.idle": "2024-04-30T09:22:36.995027Z",
     "shell.execute_reply": "2024-04-30T09:22:36.994743Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAGxCAYAAABslcJTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB/Z0lEQVR4nO3deVwU9f8H8Ndy7C4IrCByKQKa4YEnpIJ5K2pS+f1mahqiqUVmZmYlWXl8S7PDnx0emRqVF5VamWSiqVlgeaBlml/7poIK4slhcu7n9wftxrL3zuzuHO9nj3k8cvjszGdmZ2fe8/58PjMKxhgDIYQQQojMebi7AoQQQgghQkBBESGEEEIIKCgihBBCCAFAQREhhBBCCAAKigghhBBCAFBQRAghhBACgIIiQgghhBAAFBQRQgghhACgoIgQQgghBACHoOidd96BQqFAXFycyb+fO3cOCoUCb775psOVc8T8+fOhUCgc+uzJkycxf/58nDt3jt9KmbFv3z4oFArs27dPP2/ixImIjo42KHf9+nWMHTsWISEhUCgUGDlyJID6fTxixAgEBQVBoVBg5syZLqm3EDnjeLt06RLmz5+PY8eO8bZMoenfv7/Z3/DVq1ehUCgwf/58/bxjx45hxIgRaNWqFXx8fBAUFITExESsX7/eRTW2XXR0NCZOnGi1XONtdNZ54MUXX0SrVq3g5eWFpk2bAqjf//379+dtHbm5uZg/fz5u3rzJ2zIdZev+58OiRYvwxRdfOH09uvNMZmam09clN0LZt16OfnDdunUAgN9++w0//fQTevbsyVul3OXkyZNYsGAB+vfvbxSYuMpLL72Ep556ymDef/7zH2zbtg3r1q1DmzZtEBQUBAB4+umn8dNPP2HdunUICwtDeHi4O6osWZcuXcKCBQsQHR2Nrl27urs6gnDz5k1ERkbioYceQosWLXDr1i1s2LABqampOHfuHF588UV3V9FueXl5aNmypf7fzjgPfPnll3j11Vcxd+5cDB8+HCqVCgCwYsUKXpavk5ubiwULFmDixIn6wEsOFi1ahFGjRulvGIn4hIeHIy8vD23atHFrPRwKig4fPozjx49jxIgR2LFjB9auXSuJoEgITB0QJ06cQJs2bTB+/Hij+T169ODtRMAYQ2VlJXx8fHhZHpEeU5mNlJQUnD17FqtXrxZlUNSrVy+nr+PEiRMAgBkzZiAkJEQ/v0OHDlY/W1dXh9raWn0gJRR//fUXfH19XbIuoe4Drm7fvg21Wu1w64aUqFQql/wWrXGo+Wzt2rUAgNdeew1JSUnYvHkz/vrrL5NltVotXn31VbRq1QpqtRoJCQnYs2ePQZkrV67g0UcfRWRkJFQqFZo3b47evXtj9+7dBuXWrVuHLl26QK1WIygoCP/6179w6tQpq/VtnB7XaZjezczMxIMPPggAGDBgABQKhVEqb/fu3Rg0aBACAgLg6+uL3r17G22LOb///juGDRsGX19fBAcHIz09HeXl5UblGjaf6dKJu3fvxqlTp/R10jW7/fHHH/jmm2/083Xp/rKyMsyePRsxMTFQKpVo0aIFZs6ciVu3bhntl+nTp2PVqlVo3749VCoVPvroIwDAmTNnMG7cOISEhEClUqF9+/ZYvny5wed19di0aRPmzp2LiIgIBAQEYPDgwTh9+rTRtu3cuRODBg2CRqOBr68v2rdvj8WLFxuUOXz4MO677z4EBQVBrVajW7du+PTTT23ax4Btx5st27dv3z7cddddAIBJkybp9/H8+fOxY8cOKBQKHDp0SF9+y5YtUCgUGDFihMF6OnfujAceeED/b8YYVqxYga5du8LHxweBgYEYNWoU/vzzT6M62nK86ZqLf/vtNzz00EPQaDQIDQ3FI488gtLSUpv3G1fBwcHw8rJ+j3X48GGMHTsW0dHR8PHxQXR0NB566CGcP3/eoFxmZiYUCgX27t2Lxx9/HMHBwWjWrBn+/e9/49KlSwZla2pq8NxzzyEsLAy+vr64++678fPPP9tc94bnB2vngfz8fKSkpOiPm4iICIwYMQIXLlwwu/zo6Gh9sBgaGmqwvsZBpu43//rrr+OVV15BTEwMVCoV9u7dC61Wi1deeQWxsbHw8fFB06ZN0blzZ7z99tsA6o+FZ599FgAQExNjcL6w5KuvvkJiYiJ8fX3h7++PIUOGIC8vz6CM7jg7evQoRo0ahcDAQP0NnD37v7i4GI899hhatmwJpVKJmJgYLFiwALW1tTbtA1MUCgVu3bqFjz76SL/Nun165coVTJs2DR06dICfnx9CQkIwcOBAHDhwQP95xhjatm2LoUOHGi27oqICGo0GTzzxhMV9+MMPP2DQoEHw9/eHr68vkpKSsGPHDoMyumN6165deOSRR9C8eXP4+vqiqqrK7HJtOZenp6dDrVbjyJEj+nlarRaDBg1CaGgoioqKDNafk5ODSZMmISgoCE2aNMG9995rdP7JycnB/fffj5YtW0KtVuOOO+7AY489hqtXrxqUs+f889lnn6Fnz57683/r1q3xyCOP6P9urvnMnn1ry/nCKmanv/76i2k0GnbXXXcxxhhbs2YNA8AyMzMNyp09e5YBYJGRkezuu+9mW7ZsYZ999hm76667mLe3N8vNzdWXHTp0KGvevDlbvXo127dvH/viiy/Yyy+/zDZv3qwvs2jRIgaAPfTQQ2zHjh3s448/Zq1bt2YajYb997//1ZebN28ea7xZANi8efOMtiUqKoqlpaUxxhgrKSnRr2P58uUsLy+P5eXlsZKSEsYYY5988glTKBRs5MiRbOvWrWz79u0sJSWFeXp6st27d1vcZ8XFxSwkJIS1aNGCffjhhyw7O5uNHz+etWrVigFge/fu1ZdNS0tjUVFRjDHGKisrWV5eHuvWrRtr3bq1vk6lpaUsLy+PhYWFsd69e+vnV1ZWslu3brGuXbuy4OBgtnTpUrZ792729ttvM41GwwYOHMi0Wq3BfmnRogXr3Lkz27hxI/vuu+/YiRMn2G+//cY0Gg3r1KkT+/jjj9muXbvYM888wzw8PNj8+fP1n9+7dy8DwKKjo9n48ePZjh072KZNm1irVq1Y27ZtWW1trb7smjVrmEKhYP3792cbN25ku3fvZitWrGDTpk3Tl/nuu++YUqlkffr0YVlZWWznzp1s4sSJDAD78MMPLe5je443W7avtLSUffjhhwwAe/HFF/X7uLCwkJWXlzNvb2+2aNEi/TLT09OZj48Pa9KkCauurmaMMXb58mWmUCjYihUr9OWmTp3KvL292TPPPMN27tzJNm7cyNq1a8dCQ0NZcXGxvpytx5vueI+NjWUvv/wyy8nJYUuXLmUqlYpNmjTJ4j5jjLF+/fqxjh07mvzblStXzP526urqWE1NDSspKWHLly9nXl5ebNWqVVbX99lnn7GXX36Zbdu2je3fv59t3ryZ9evXjzVv3pxduXJFX06371u3bs2efPJJ9u2337I1a9awwMBANmDAAINlpqWlMYVCwZ599lm2a9cutnTpUtaiRQsWEBCg/31b0nAbLZ0HKioqWLNmzVhCQgL79NNP2f79+1lWVhZLT09nJ0+eNLv8o0ePssmTJzMAbOfOnfrjiLH6/d+vXz99Wd1x3KJFCzZgwAD2+eefs127drGzZ8+yxYsXM09PTzZv3jy2Z88etnPnTrZs2TL9MVtYWMiefPJJBoBt3brV4HxhzoYNGxgAlpyczL744guWlZXF4uPjmVKpZAcOHNCX0x1nUVFR7Pnnn2c5OTnsiy++sGv/FxUVscjISBYVFcXef/99tnv3bvaf//yHqVQqNnHiRJv2gSl5eXnMx8eH3XPPPfpt/u233xhjjP3+++/s8ccfZ5s3b2b79u1jX3/9NZs8eTLz8PAwOO++/fbbTKFQGFxLGGNs+fLlDIB+ebq6NTwf7du3j3l7e7P4+HiWlZXFvvjiC5acnMwUCoXBNUx3TLdo0YI9+uij7JtvvmGff/65wXmyIVvP5bdv32Zdu3ZlrVu3Zjdu3GCMMfbyyy8zDw8PtmvXLqP1R0ZGskceeYR98803bPXq1SwkJIRFRkbqP8sYYytXrmSLFy9mX331Fdu/fz/76KOPWJcuXVhsbKz+/NbwuLB2/snNzWUKhYKNHTuWZWdns++++459+OGHLDU11eh757JvbTlfWGN3UPTxxx8zAPoTYHl5OfPz82N9+vQxKKfbwIiICHb79m39/LKyMhYUFMQGDx6sn+fn58dmzpxpdp03btzQH/QNFRQUMJVKxcaNG6ef52hQxFj9CbtxkMJY/cEZFBTE7r33XoP5dXV1rEuXLqxHjx5m684YY88//zxTKBTs2LFjBvOHDBliMSjSMXfhioqKYiNGjDCYt3jxYubh4cEOHTpkMP/zzz9nAFh2drZ+HgCm0WjY9evXDcoOHTqUtWzZ0uhkOn36dKZWq/XldUFR4+/l008/ZQBYXl4eY6z+GAkICGB33323QVDWWLt27Vi3bt1YTU2NwfyUlBQWHh7O6urqzH7WnuPN1u07dOiQ2YDs7rvvZgMHDtT/+4477mDPPvss8/DwYPv372eM/XPB0Z1o8/LyGAD21ltvGSyrsLCQ+fj4sOeee44xZt/xpjveX3/9dYOy06ZNY2q12uL+ZszxoOixxx5jABgAplQqDQI/e9TW1rKKigrWpEkT9vbbb+vn605yDYNmxhh7/fXXGQBWVFTEGGPs1KlTDAB7+umnDcrp9r29QRFj5s8Dhw8fZgD0wYA9dN9Tw8CPMfNBUZs2bQwuPozV/w66du1qcT1vvPEGA2A2gGiorq6ORUREsE6dOhn8tsrLy1lISAhLSkoyqv/LL79ssAx79v9jjz3G/Pz82Pnz5w3KvvnmmyYDD1P7wJwmTZrY9F3X1taympoaNmjQIPavf/1LP7+srIz5+/uzp556yqB8hw4dDC6qpi7cvXr1YiEhIay8vNxgPXFxcaxly5b636DumJ4wYYJN22TPufzMmTMsICCAjRw5ku3evZt5eHiwF1980eBzuvU33G7GGPvxxx8ZAPbKK6+YrIdWq2U1NTXs/PnzDAD78ssv9X+z9fyj+45v3rxpdnv52LfWzhe2sLv5bO3atfDx8cHYsWMBAH5+fnjwwQdx4MABnDlzxqj8v//9b6jVav2//f39ce+99+L7779HXV0dAKBHjx7IzMzEK6+8goMHD6KmpsZgGXl5ebh9+7bRSIbIyEgMHDjQ5iYsR+Xm5uL69etIS0tDbW2tftJqtRg2bBgOHTpk1DTV0N69e9GxY0d06dLFYP64ceN4r+vXX3+NuLg4dO3a1aCuQ4cONZlKHzhwIAIDA/X/rqysxJ49e/Cvf/0Lvr6+Bsu45557UFlZiYMHDxos47777jP4d+fOnQFA3ySSm5uLsrIyTJs2zWzb+R9//IHff/9d32+q8XqLiopMNsk1Zu14c2T7TBk0aBB+/PFH3L59G+fPn8cff/yBsWPHomvXrsjJyQFQ3/zVqlUrtG3bFkD9d6NQKPDwww8brDcsLAxdunTRfzeOHG+mvoPKykqUlJRY3RZHvPDCCzh06BB27NiBRx55BNOnT7dp5F9FRQWef/553HHHHfDy8oKXlxf8/Pxw69Ytk03h1o4tXZNK4/52o0ePtqk5zx533HEHAgMD8fzzz2PVqlU4efIkr8tv6L777oO3t7fBvB49euD48eOYNm0avv32W5SVlXFax+nTp3Hp0iWkpqbCw+OfS4Gfnx8eeOABHDx40KhbRMOmYMC+/f/1119jwIABiIiIMDiuhw8fDgDYv3+/QXlT+8ARq1atQvfu3aFWq+Hl5QVvb2/s2bPH4Hjz9/fHpEmTkJmZqf9tfffddzh58iSmT59udtm3bt3CTz/9hFGjRsHPz08/39PTE6mpqbhw4YLReavxPjTHnnP5HXfcgQ8++ABffPEFUlJS0KdPH5NdRgDj7yopKQlRUVEGzZMlJSVIT09HZGSkfp9FRUUBgM2/04bnH11XhNGjR+PTTz/FxYsXrW6/I/vW2vnCFnYFRX/88Qe+//57jBgxAowx3Lx5Ezdv3sSoUaMA/DMiraGwsDCT86qrq1FRUQEAyMrKQlpaGtasWYPExEQEBQVhwoQJKC4uBgBcu3YNAEyOroqIiND/3VkuX74MABg1ahS8vb0NpiVLloAxhuvXr5v9/LVr18zuB2fU9ZdffjGqp7+/PxhjRm3CjffptWvXUFtbi3fffddoGffccw8AGC2jWbNmBv/WdYa8ffs2gPp2fQAGI3xM1RsAZs+ebbTeadOmmVyvKdaON0e2z5TBgwejqqoKP/zwA3JychAcHIxu3bph8ODB+r5we/bsweDBgw22kTGG0NBQo3UfPHhQv15Hjjdr34E5Xl5e+puTxnT9PExdmFq1aoWEhATcc889WLlyJR599FFkZGTov2tzxo0bh/feew9TpkzBt99+i59//hmHDh1C8+bNTdbV2nbpfvuNv3cvLy+jz3Kl0Wiwf/9+dO3aFS+88AI6duyIiIgIzJs3z+hGjitT57qMjAy8+eabOHjwIIYPH45mzZph0KBBOHz4sEPrsHZe1Wq1uHHjhsV62bP/L1++jO3btxsd0x07dgRg/LvjYzTt0qVL8fjjj6Nnz57YsmULDh48iEOHDmHYsGFGx9uTTz6J8vJybNiwAQDw3nvvoWXLlrj//vvNLv/GjRtgjJndhwCMrk+2bpe95/IRI0YgNDQUlZWVmDVrFjw9PU0u19w5UldPrVaL5ORkbN26Fc899xz27NmDn3/+WX+z6MjvtG/fvvjiiy9QW1uLCRMmoGXLloiLi8OmTZvMbr8j+9bR82BDdt1KrVu3DowxfP755/j888+N/v7RRx/hlVdeMfgydIFNQ8XFxVAqlfroLzg4GMuWLcOyZctQUFCAr776CnPmzEFJSQl27typ31Bdh7GGLl26hODgYIv1VqlUJjuz2RpM6Zb/7rvvmu0dHxoaavbzzZo1M7sf+BYcHAwfHx+TAaru7w01ztwEBgbqI3FznQtjYmLsqlPz5s0BwGJnVF29MjIy8O9//9tkmdjYWKvrsna8eXt787J9PXv2hJ+fH3bv3o1z585h0KBBUCgUGDRoEN566y0cOnQIBQUFBkFRcHAwFAoFDhw4YHIUjW4e1+PNHqGhoTh06BAYY0bHgu5uzpZ19ejRA6tWrcKff/6p/74bKy0txddff4158+Zhzpw5+vlVVVUWbyos0Z0biouL0aJFC/382tpap9wsderUCZs3bwZjDL/88gsyMzOxcOFC+Pj4GGwTV6Yyql5eXpg1axZmzZqFmzdvYvfu3XjhhRcwdOhQFBYW2j0SzNp51cPDwyCLbKpe9uz/4OBgdO7cGa+++qrJ+ugudObW5Yj169ejf//+WLlypcF8U4Nc7rjjDgwfPhzLly/H8OHD8dVXX2HBggVmgwug/nzp4eFhdh8C1s+55th7LtcN3unYsSNmzJiBPn36GH1/gPlz5B133AGgfqTk8ePHkZmZibS0NH2ZP/74w6Z6m3P//ffj/vvvR1VVFQ4ePIjFixdj3LhxiI6ORmJiolF5R/YtH2wOiurq6vDRRx+hTZs2WLNmjdHfv/76a7z11lv45ptvkJKSop+/detWvPHGG/omjfLycmzfvh19+vQxebC1atUK06dPx549e/Djjz8CABITE+Hj44P169frR4YA9RfZ7777Tp+pMic6Ohq//PKLwbzvvvtOn6nSMRdV9u7dG02bNrWaSjVnwIABeP3113H8+HGDJrSNGzfavSxrUlJSsGjRIjRr1szu4AUAfH19MWDAAOTn56Nz585QKpWc65SUlASNRoNVq1Zh7NixJk8KsbGxaNu2LY4fP45FixY5vC5rx5s922fpLsPb2xt9+/ZFTk4OCgsL8dprrwEA+vTpAy8vL7z44ov6IEknJSUFr732Gi5evIjRo0ebXS/X480egwcPxsaNG7Fz5059M4bOp59+Cg8PDwwcONDqcvbu3QsPDw+0bt3abBmFQgHGmFFAuGbNGrPZKmt0o4w2bNiA+Ph4g7o3HNFkD1vuLhUKBbp06YL/+7//Q2ZmJo4ePerQuhzVtGlTjBo1ChcvXsTMmTNx7tw5dOjQwa4749jYWLRo0QIbN27E7Nmz9b/LW7duYcuWLfoRaZbYs/9TUlKQnZ2NNm3amLxYc6FSqUxus0KhMDrefvnlF+Tl5SEyMtKo/FNPPYXk5GSkpaXB09MTU6dOtbjeJk2aoGfPnti6dSvefPNN/eNMtFot1q9fj5YtW+LOO+90aJvsOZevWbMG69evx7p169CvXz90794dkyZNMvlAyw0bNhg04eXm5uL8+fOYMmUKgH+Ctsb77f3333doOxpTqVTo168fmjZtim+//Rb5+fkmgyJn7ltLbA6KvvnmG1y6dAlLliwx+QTWuLg4vPfee1i7dq1BUOTp6YkhQ4Zg1qxZ0Gq1WLJkCcrKyrBgwQIA9XePAwYMwLhx49CuXTv4+/vj0KFD2Llzpz5j0LRpU7z00kt44YUXMGHCBDz00EO4du0aFixYALVajXnz5lmse2pqKl566SW8/PLL6NevH06ePIn33nsPGo3GaBsAYPXq1fD394darUZMTAyaNWuGd999F2lpabh+/TpGjRqFkJAQXLlyBcePH8eVK1eM7kQamjlzJtatW4cRI0bglVdeQWhoKDZs2IDff//dpn1vj5kzZ2LLli3o27cvnn76aXTu3BlarRYFBQXYtWsXnnnmGavPlHr77bdx9913o0+fPnj88ccRHR2N8vJy/PHHH9i+fTu+++47u+rk5+eHt956C1OmTMHgwYMxdepUhIaG4o8//sDx48fx3nvvAaj/0Q0fPhxDhw7FxIkT0aJFC1y/fh2nTp3C0aNH8dlnn1ldl7XjzZ7ta9OmDXx8fLBhwwa0b98efn5+iIiI0N/RDho0CM888wwA6DNCPj4+SEpKwq5du9C5c2eDZ9L07t0bjz76KCZNmoTDhw+jb9++aNKkCYqKivDDDz+gU6dOePzxx+Hn58fpeLPH+PHjsWLFCowePRpz5szBXXfdhdu3byM7OxsffPABnnzySYNA59FHH0VAQAB69OiB0NBQXL16FZ999hmysrLw7LPPms0SAUBAQAD69u2LN954A8HBwYiOjsb+/fuxdu1ahx802L59ezz88MNYtmwZvL29MXjwYJw4cQJvvvkmAgICHFqmufNAXl4eVqxYgZEjR6J169ZgjGHr1q24efMmhgwZ4tC67HHvvfciLi4OCQkJaN68Oc6fP49ly5YhKipK32+tU6dOAOqP8bS0NHh7eyM2Nhb+/v5Gy/Pw8MDrr7+O8ePHIyUlBY899hiqqqrwxhtv4ObNm/pA3xJ79v/ChQuRk5ODpKQkzJgxA7GxsaisrMS5c+eQnZ2NVatWWWxit6RTp07Yt28ftm/fjvDwcPj7+yM2NhYpKSn4z3/+g3nz5qFfv344ffo0Fi5ciJiYGJNB85AhQ9ChQwfs3bsXDz/8sMHv15zFixdjyJAhGDBgAGbPng2lUokVK1bgxIkT2LRpk8MZL1vP5b/++itmzJiBtLQ0TJo0CUB9399Ro0Zh2bJlRm86OHz4MKZMmYIHH3wQhYWFmDt3Llq0aKHvptCuXTu0adMGc+bMAWMMQUFB2L59u76vpCNefvllXLhwAYMGDULLli1x8+ZNvP322/D29ka/fv3Mfs5Z+9YiW3tkjxw5kimVSv0QdVPGjh3LvLy8WHFxsb4n+ZIlS9iCBQtYy5YtmVKpZN26dWPffvut/jOVlZUsPT2dde7cmQUEBDAfHx8WGxvL5s2bx27dumWw/DVr1rDOnTszpVLJNBoNu//++/UjFnRMjT6rqqpizz33HIuMjGQ+Pj6sX79+7NixY0ajzxhjbNmyZSwmJoZ5enoa9YTfv38/GzFiBAsKCmLe3t6sRYsWbMSIEeyzzz6zuv9OnjzJhgwZwtRqNQsKCmKTJ09mX375Je+jzxhjrKKigr344ossNjZWv686derEnn76aYNh3wDYE088YbK+Z8+eZY888ghr0aIF8/b2Zs2bN2dJSUkGIxR0o88ab7+pUQSMMZadnc369evHmjRpwnx9fVmHDh3YkiVLDMocP36cjR49moWEhDBvb28WFhbGBg4caHW4t63Hmz3bxxhjmzZtYu3atWPe3t5Go5SOHz/OALC2bdsafObVV19lANisWbNM1nXdunWsZ8+erEmTJszHx4e1adOGTZgwgR0+fNignC3Hm7lRTbrRGLaMQiorK2PPPfcca9u2LVMqlczX15clJCSwVatWGY1eW7duHevTpw8LDg5mXl5erGnTpqxfv37sk08+sboexhi7cOECe+CBB1hgYCDz9/dnw4YNYydOnDD6Lerq33jkje6Ya/ibqaqqYs888wwLCQlharWa9erVi+Xl5Zn8fZvS+HtlzPR54Pfff2cPPfQQa9OmDfPx8WEajYb16NHD6HEkptg7+uyNN94wWsZbb73FkpKSWHBwMFMqlaxVq1Zs8uTJ7Ny5cwblMjIyWEREBPPw8DA5iq6xL774gvXs2ZOp1WrWpEkTNmjQIPbjjz/aVH/G7Nv/V65cYTNmzGAxMTHM29ubBQUFsfj4eDZ37lxWUVFhdR+Yc+zYMda7d2/m6+vLAOj3aVVVFZs9ezZr0aIFU6vVrHv37uyLL74weZ7VmT9/PgPADh48aPQ3c+e2AwcOsIEDB+p/07169WLbt283KGPumLbE2rm8oqKCtWvXjnXo0MHoevnEE08wb29v9tNPPxmsf9euXSw1NZU1bdpUP6r7zJkzBp/VXa/8/f1ZYGAge/DBB1lBQYHRb8XW88/XX3/Nhg8fzlq0aMGUSiULCQlh99xzj8FjH5yxb02dL6xRMMYY/6EWIYQQIj4JCQlGD2eVgszMTEyaNAmHDh1CQkKCu6sjWPyOWSWEEEJEpqysDCdOnMDXX3+NI0eOYNu2be6uEnETCooIIYTI2tGjRzFgwAA0a9YM8+bNoxfLyhg1nxFCCCGEwMEXwgrd999/j3vvvRcRERFQKBQmhyU2tn//fsTHx0OtVqN169ZYtWqVUZktW7boh7126NCBUqyEEEKIg5x1reZCkkHRrVu30KVLF/1Qb2vOnj2Le+65B3369EF+fj5eeOEFzJgxA1u2bNGXycvLw5gxY5Camorjx48jNTUVo0ePxk8//eSszSCEEEIkyxnXaq4k33ymUCiwbds2i23Ezz//PL766iuDd7qkp6fj+PHjyMvLAwCMGTMGZWVl+Oabb/Rlhg0bhsDAQIuPKieEEEKIZXxdq7mijtaozwIlJycbzBs6dCjWrl2LmpoaeHt7Iy8vD08//bRRmWXLlpldblVVlcHrRbRaLa5fv45mzZo556FThBBCJIMxhvLyckRERBi8tJdvlZWVqK6u5rwcZuJ1QSqVyuRrjRxhy7WaKwqKUP/el8bvdwoNDUVtbS2uXr2K8PBws2Usvb9s8eLFBk9SJoQQQuxVWFjo8NO+ramsrERMlB+KSxx71U5Dfn5+Rq/PmjdvHubPn8952YBt12quKCj6W+PoVteq2HC+qTKWMj4ZGRmYNWuW/t+lpaVo1aoVCgsL7X4FwTun/2VyfuFt4zeBX6zUmCgJXL7lZ3J+6V9qs+utvGUhwq8wf/h43bJ8V+N5y3qmzLvCapH6crdsKwcAygqt7YV1nyl1/GShLOPv7eme5ZW8LasxxU3jF2Q6k/bGTZeuzxU8ApvyujzW1PjVHPaq8zf/2zanOsC+u+1qjfkXphqU87Oe6ahpYn05NaZPYwCAuiaWe4PUNjHz+/cz/548dRPjl4lrfE3/FkObGJ+0WqhL9f//Rlfj94ZaUlZWhsjISJOvaeFLdXU1ikvqcPZIFAL8Hc9GlZVrERN/3uj6xleWSMeWazUXFBQBCAsLM8r4lJSUwMvLS/8WaHNlLL1B3FzaMCAgwO6gSO1n/FWdvx0MpYkThLen8UtOi2/5w8vECefmLR94mnnn4+0KFTx8zFSowgsw8zevCg/AwrnYq0IBWPmdeFfAahl9ORvfWasq1wJ2ZleVpXWAgylZZWkN4GXbBcMaz7LbgCe/JxcdxY1ywMM5yzZHq+D+omHBufkXPIJ4fNFpWTVYILeLocLL/qDI6y+gWmP7Ma/1tu0Y96oCqqxcdD1rLAc9AKC1cKh61gK1fuYDIw8fM0FRHcwGRp6+xhfbCqjRtInxC2ivQYWwJoY3GCVojpY+NwHA4ffxuaK7RYC/B6egSL8cB65vtrLlWs2VJEef2SsxMdHoZXe7du1CQkKCvo3SXJmkpCSX1bOh87eDTc6/cLspL8u/XeFYhsiVbM0kOYpThqiUxwxRmfU3njtKccO1GSIA0F6/4fJ1ypUzjx0dLr8TU6z9rq393avCfADhVWH/Jc/iuVBC6piW8+RstlyruZJkUFRRUYFjx47h2LFjAOqH8R07dgwFBQUA6pu1JkyYoC+fnp6O8+fPY9asWTh16hTWrVuHtWvXYvbs2foyTz31FHbt2oUlS5bg999/x5IlS7B7926jNxC7k7mAqPiW6bvNm7fMpYEssBIQWTvpWDph6Tgj2FGV2/eDpYCIOILvgM9d35G9x7Ctvxd7f4cuZefNnrnzp6nzLV83q86kBeM82csZ12quJBkUHT58GN26dUO3bt0AALNmzUK3bt3w8ssvAwCKior0Ox0AYmJikJ2djX379qFr1674z3/+g3feeQcPPPCAvkxSUhI2b96MDz/8EJ07d0ZmZiaysrLQs2dP124czGeJTHEkIHL0zsiVAZE9gZOgT8QWuOIu39UoS2Q/roGRo8cRn8G9vaSYLRJ6YKTl4T97OeNazZXkn1MkJGVlZdBoNCgtLbW7zfXNU0P1/28qKOIrS+Ros5ktJxprQZFQAiIhZImcHRC5KwMhh6CI135Ff+Pat6guwIGsMOzrW1Rf3rb+Rdb6FgHW+xZZK2Opb1Gtn4Xzgpm+RT5+xh2uAZjsWwTAqG8RAGzrvdz8ek3gcs2wdx2XTrfk3NE6IvaCU+vqCpLMFEmZMwMiizj2I7IlS2QLCoi4o4DIuZyxnVLLFtny2+TajC7EbJGQ1THGeZICCopExJ5mM0cIvdmMAiLuqB+RfImx07U1ThlswVPfIrFxR58iIaKgSOSE3mzmjo7VFBCZ5s6ASC5ZIh0hZoscRdmif9hz42juHEyEjYIikXB252qznDz83hn9iOxFARGRC8oW2UiG2SItGOo4TJQpIm7nkmcSWcA1SySEjtUUEDmf3LJEOnLOFvE5RJ+yRa5BzWf1KCgSKSk0m9mCAiJu3B0QEeERWqdrW3Adou8QGWaLCAVFkiKUZjMxdKymgMg15Jol0hFqtkhIzWiueI4YZYuso9Fn9SgoEiGhN5tZI4SO1Y6igMh2cg+IpEjIna4pW8SNlodJCigokgg5N5vZy9EsEQVERCikli1yBcoWEVtQUCQy9mSJqNnMGAVErkFZon84a19Qp2v7/+7QTZsT34kmJFxGnukmKaCgSETs7VxtjpCbzSggchwFRMReUny/nqNZZrk/5bqOcZ+kgIIikRNLs5m7+xFRQOQaFBCZRtkikQ/Rl0HfIupTVI+CIpFw+huWRdSPyJUjzfggl4CIiJMYh+hb48w+iY1JKVtEKCgSNV6zRA5ydT8isQ29l1NARFkiy4ScLRJSp2t3Z4vMkni2SAsF6jhMWvBzc+xuFBRJDN+dq4XWj8hVKCCyDwVE8iTFbJFc+xZpGfdJCigoEineOldTPyITn6OAyB4UENmOskWULSLCRkGRCIml2cwWQms2o4DIPhQQCYdYOl3bio+HsroyWyR2XJrOdJMUyO+bJ/8QSbMZBUT/oIBI/IS831zR6ZrPgQ9CyhaZu/kUS7aIgqJ6FBSJjBiyRBQQOQcFRMQaakaz/+/myDFbRCgokgRXPpOIr2YzZ6CAyDUoIOJOivtQyJ2uLaFsUT0tU3CepICCIhHh7THxEmw2o4DINaR4MXcX6nRN2SIhoeazevL5xiWKms0cJ/SASHGjnAIi4hCpdbp2NsoWER0KikTC3VkiIQdEjmSJxBAQCQkFRM4h5P1K2aJ6cskW1cGD8yQF0tgKmbI7S+TEZjN3oYDI+YR84SbmuasZjbJF4nyYI+PYn4hRnyIiF0LNErnrnWYUEBE+SXEfO2OIPmWLnIv6FNWT/jctUa7KEkktIOJ6F0sBERETsXS6FjUJZYsIBUWixGfnamffAVFAZBsKiOTNmftbDJ2uhZItsnQTKPVsUR3z4DxJgX0vciHiZOf7enRc9W4zMYw0o4CIyJln2W3UBdDoKbMqvAC/WnfXghMtFNByyJNoIY03wkojtJMRKTab2cMdHaudFRAJccg9BUTuQ9ki6WWLqAlNfChTRJxCKM1mQg6IhISCIWKNK7JFytI6VGs8nboOYhrXztLU0Zq4nJiyRLaggEgYKCASDrlni2wlyGyRyDtcU5+ietLYChlzdedqvprNKCByP2ouEyYhfydCeqAjIc5AQZFI2P2IeCd1rrbGGf2I7EUBkXVCvvAS55HSAx35yBZZIreHOdZ3tOY2SQEFRSLGd5ZIKs1mXFBARIRAis1o9nDV754e5vgPLcdXfHAZuSYk0tgKYsjBLJElcmg2c0ZAJKQRZtRcRvhC2SILRJwtIhIPilasWIGYmBio1WrEx8fjwIEDZstOnDgRCoXCaOrYsaO+TGZmpskylZWVrtgcA67uXG0N381mUgmIhIKCIfGRYrbIGUP0uaJsUT3qaF1PGlthQlZWFmbOnIm5c+ciPz8fffr0wfDhw1FQUGCy/Ntvv42ioiL9VFhYiKCgIDz44IMG5QICAgzKFRUVQa1Wu2KTOHFms5kzht/bS2gvoRRKQETZIXET8ncnlNd/OPvhr1xvGMVC+3cTGJdJCqSxFSYsXboUkydPxpQpU9C+fXssW7YMkZGRWLlypcnyGo0GYWFh+unw4cO4ceMGJk2aZFBOoVAYlAsLC3PF5hhwJEvkKDH0IxJax2ohBUSEmCOnbJHLXxTrhHOxs9UxBedJCiQZFFVXV+PIkSNITk42mJ+cnIzc3FyblrF27VoMHjwYUVFRBvMrKioQFRWFli1bIiUlBfn5+WaXUVVVhbKyMoPJHdzdbEYBkWtRdkhahNyMJpdsEZEPSQZFV69eRV1dHUJDQw3mh4aGori42Orni4qK8M0332DKlCkG89u1a4fMzEx89dVX2LRpE9RqNXr37o0zZ86YXM7ixYuh0Wj0U2RkpOMb9TchZYmE0I+IC6kGRIS4krM7XQshW8Rnh2uh4jLyTDdJgTS2wgyFwvBAZowZzTMlMzMTTZs2xciRIw3m9+rVCw8//DC6dOmCPn364NNPP8Wdd96Jd9991+RyMjIyUFpaqp8KCwsd3hZHOTNLZAuh9iOSWkBE2SFpE3K2SCjclS2SSodrLfPgPEmBuEJZGwUHB8PT09MoK1RSUmKUPWqMMYZ169YhNTUVSqXSYlkPDw/cddddZjNFKpUKKhV/wzDtzRK5u3O1UJvNpBgQEeJOjrwXTVlag2qNt41l+XknmncFUOPn2N+9KhSo9ZPGm+CJedII7RpRKpWIj49HTk6OwfycnBwkJSVZ/Oz+/fvxxx9/YPLkyVbXwxjDsWPHEB4ezqm+tqBms4blKSACKDskN5Qtsk5w2SIRNaFR81k98Xxjdpo1axZSU1ORkJCAxMRErF69GgUFBUhPTwdQ37R18eJFfPzxxwafW7t2LXr27Im4uDijZS5YsAC9evVC27ZtUVZWhnfeeQfHjh3D8uXLXbJN9hBLs5lYAyJ3X0QoGJIn7fUb8AgKdHc1TJJKtsgSKWeLtACnEWRS6eou2aBozJgxuHbtGhYuXIiioiLExcUhOztbP5qsqKjI6JlFpaWl2LJlC95++22Ty7x58yYeffRRFBcXQ6PRoFu3bvj+++/Ro0cPp2+PSQLNEjkrIOJCKgERBUPEWRQ3ysEC/d1dDc5U5VpU+XN76bUjQZNXhQdq/Uyc+yq8AL9ah+tDXEvBGJNm2CtAZWVl0Gg0KC0tRUBAgF2fjf7kNeOZPPclcvWrPFyVJaKAiEiNM7NFXAMje7NFAGzOFtWXtZ4tsiUoshT4WPqbpUyRyaAI0AdF51LnWK1XQ1yuGfauY+XRu+Dj5/iN9u2KWjze/ZBT6+oKks0USZ4DWSIhPZOIAiLbUTBEiH2cmS1yqAlNBNkirq/qoNd8EEFydHiokJvNKCAixJCQO12L5blFziCV4flyRt+gGLkhS+Qucg2IaGQZsYaOD27oYY6GtFBwnqSAgiIJkWKWyN3cFRAR4k5SyBYJbni+wJl66729kxQIO3QlxlzcuVrOzWauDogoGCL2EvIQfTHgMjxfarg+a0gqzymSxlYQyaGAiBD3Enq2yBZcs0W8N6ERwaOgSExkkiWigIgQ2wm507UzuepFsY4QYxOalik4T1Igvm+OCAYFRNxQZ2rCF6EeR5QtEg8tx1d8aCUSTkhjK+RAgFkiW7iiY7VYAyJCxEAO2SJCdCgoIkb4bDazB993iPZw1YmfskPEWYR6XIkhW0RNaICWeXCepEAaWyFTznydhyVSazZzZUBEiBjJPVskhya0Oig4T1JAQZEEcfmRyq3ZjAIiIiVCPc74fJiqs1C2iAAUFImWsx7UaI1Qms3EEhBRcxlxNWcdb67OFlGHa9ei5rN69PBGiRFb52qpB0SEkHqeZbdRF+DjlGUrS+tQrfF0yrLlog7g1AQmlS7v0gjtZMZdKVlb7rQoIKpHARFxJ8oWmUYdrok1lCmSEGdmidz1HiG+OfukTsEQIeaJPVtk6bUgXhUK1Poxp67fmbg2gUml+UwaWyEjzrjrkEuzGQVERE4oW+QYZ2SLxIBeCFtPGltBBDEE31YUEBEiX84ciWbLjZk7st5iaEJjUEDLYWIO9kdasWIFYmJioFarER8fjwMHDlgsv2HDBnTp0gW+vr4IDw/HpEmTcO3aNYfWbYrwvymiJ5Uskb0oICLEMZQt4p/UR6G5UlZWFmbOnIm5c+ciPz8fffr0wfDhw1FQUGCy/A8//IAJEyZg8uTJ+O233/DZZ5/h0KFDmDJlCm91oqBIAtydJXLFQxod5cyTNw23J2IgxGPU3c8tog7XxtzRfLZ06VJMnjwZU6ZMQfv27bFs2TJERkZi5cqVJssfPHgQ0dHRmDFjBmJiYnD33Xfjsccew+HDh7luvp6wvyWi544skdibzZwdEBEiZ0J9yrW7n3AtVqbeem/vBABlZWUGU1VVlcn1VVdX48iRI0hOTjaYn5ycjNzcXJOfSUpKwoULF5CdnQ3GGC5fvozPP/8cI0aM4G0/UFAkcs7MEtlCyM1mzkIBEREbIR6z9v6+hdSEZoncm9AiIyOh0Wj00+LFi02Wu3r1Kurq6hAaGmowPzQ0FMXFxSY/k5SUhA0bNmDMmDFQKpUICwtD06ZN8e677/JWfwqKZIqPLJEzm82E2o9IiBcXQtxFzNkidzShCVkdPDhPAFBYWIjS0lL9lJGRYXG9CoVhEMkYM5qnc/LkScyYMQMvv/wyjhw5gp07d+Ls2bNIT0/nZyeAnlMkau7OEtlKKv2IKCAiYqa9fgMeQYHurgYnytIaVGu83V0NAJafWSRGDZvAHP08AAQEBCAgIMBq+eDgYHh6ehplhUpKSoyyRzqLFy9G79698eyzzwIAOnfujCZNmqBPnz545ZVXEB4e7nD9dShTJEOuzhLZS4j9iCggIsQ0Lr85oTaR6zh6cyn3JjRbKJVKxMfHIycnx2B+Tk4OkpKSTH7mr7/+goeHYdji6Vn/wE7G+HlwJgVFEiTmLJHQAiIaYUakRArHsq3nFKE+s0iotPDgPNlr1qxZWLNmDdatW4dTp07h6aefRkFBgb45LCMjAxMmTNCXv/fee7F161asXLkSf/75J3788UfMmDEDPXr0QEREBC/7gZrPRMrROxFXZomkEBARQqxT3CgHC/R36LPOfPWHs0mpCa2OKVDHofnMkc+OGTMG165dw8KFC1FUVIS4uDhkZ2cjKioKAFBUVGTwzKKJEyeivLwc7733Hp555hk0bdoUAwcOxJIlSxyud2MUFEmMkLJEYkYBEZEqKfQtspUt70NTlWtR5W8+y+Fo4CP2d6G5yrRp0zBt2jSTf8vMzDSa9+STT+LJJ590Wn2o+UyEKEtkiO8sEQVEhNjPlX2LxDI8X0z4ek6R2FGmSEKcnSWigIgQaZBTtsjZpNKExpgHpzfdM3ohLHEHd2aJhIYCIkKINe58ZpGYRqHVQcF5kgIKiiSCS5ZIjM1mFBARwo0zjnlqQiNiR81nIuKsLJG7COUZJRQQEUK4kEITmpaB48MbeayMG1GmSALEmiVyFJ9ZIgqIiJzJJVtETWjWaf/uU8RlkgLKFImE0LNEYmw2o4CImOPBtIjTliCI3cZ1hQ9OeIRAq5DGSZ8QYp6kf+UrVqxATEwM1Go14uPjceDAAbNl9+3bB4VCYTT9/vvvBuW2bNmCDh06QKVSoUOHDti2bZuzN8Np3PU6DwqIiJD1ri3Ax5Xb8EZVDjKqf8AbVTn4uHIbetcWWP+wCAktWyRWQu2mYCstFJwnKZBsUJSVlYWZM2di7ty5yM/PR58+fTB8+HCDp2Oacvr0aRQVFemntm3b6v+Wl5eHMWPGIDU1FcePH0dqaipGjx6Nn376ydmbY5alH6JQs0TuRgERMad3bQFeqt6PYPaXwfxm7C+8VL1fsoGRkEixCU0MdE+05jJJgWSDoqVLl2Ly5MmYMmUK2rdvj2XLliEyMhIrV660+LmQkBCEhYXpJ93L5gBg2bJlGDJkCDIyMtCuXTtkZGRg0KBBWLZsmZO3hn9yzRJRQETM8WBaPF5zCACM7nk9ADAA6TWH4MHE9/gKayhb5Fxi6VdEJBoUVVdX48iRI0hOTjaYn5ycjNzcXIuf7datG8LDwzFo0CDs3bvX4G95eXlGyxw6dKjZZVZVVaGsrMxg4pPYskQUEBEhi9OWoDn7y2wjgAeAEPYX4rQlrqyWLAllZKqcUEfretLYikauXr2Kuro6hIaGGswPDQ1FcXGxyc+Eh4dj9erV2LJlC7Zu3YrY2FgMGjQI33//vb5McXGxXctcvHgxNBqNfoqMjOS4Za7jjCyRIyggIq4SxGy7ENtaTmzE/BtxZROaJZZuRoWeLdKC42s+JNKnSNKjzxQKwy+JMWY0Tyc2NhaxsbH6fycmJqKwsBBvvvkm+vbt69AyMzIyMGvWLP2/y8rKeAuMuGSJ+Hp6tSuyRHwQ88meuM51hW1vare1HKm/qWGB/u6uBq+k8EwiYp4kM0XBwcHw9PQ0yuCUlJQYZXos6dWrF86cOaP/d1hYmF3LVKlUCAgIMJjEQGpZIkJsccIjBFcUvjB3y6AFUKLwxQmPEFdWy6WEdANBTWiuxTiOPGMSyRRJMihSKpWIj49HTk6OwfycnBwkJSXZvJz8/HyEh4fr/52YmGi0zF27dtm1THcTS5aIms2Iq2kVHljpfRcUgFFgpEV95+tV3nfR84oEis8mNDni1HT29yQFkm0+mzVrFlJTU5GQkIDExESsXr0aBQUFSE9PB1DftHXx4kV8/PHHAOpHlkVHR6Njx46orq7G+vXrsWXLFmzZskW/zKeeegp9+/bFkiVLcP/99+PLL7/E7t278cMPP7h025zZwdqWEwY1mxGp+tGrFf6Dfni85hCaNxiWf1Xhi1Xed+FHr1ZurJ1raK/fgEdQIG/LE2MTmqpciyp/88GvpSY0sTavce0sLZWO1pINisaMGYNr165h4cKFKCoqQlxcHLKzsxEVFQUAKCoqMnhmUXV1NWbPno2LFy/Cx8cHHTt2xI4dO3DPPffoyyQlJWHz5s148cUX8dJLL6FNmzbIyspCz549Xb59juArS+RsfGSJKCAijvrRqxXyPFvSE60FwLPsNuoCqA8XcR0FY0wir3ETvrKyMmg0GpSWltrdv6jdvP8D4NwO1kLIElFARIj78ZkpAsApU2RvUFSt8baxnKfFv1vKFAGWs0Hm/vb7gqetVcsAl2uGveu4f9cj8G6idHg5Nbeq8WXyOqfW1RXo1kcmxJIl4ooCIkK44/t3JMZBE84ami9U9JqPehQUSYQU+hKJ8cRJCHEuoY5CE2PgQ6yjoEhEHP0RujpLRM1mhIifWLNFYnsXo1DQ6LN6ku1oTWzHd5bIHSggIoTYQ1laZ7VfkZxwDWykEhRRpkgkhPCeM1u4I0tEAREh0uauJjS59SsiFBRJnqtHnNmLAiJChEusTWiuIqXAh5rP6lHzmcgJ6Ucp1A6RhBB5UZbW2DQ0n5rQ/kHNZ/UoUyRhlCUihHBFvzMiJxQUEV7YmyWigIgQeXL0ty/GfkViwsDtWUVSeQo0BUUiZqnpjI8fslBHnFFARAixxlXnL7EMgrGG+hTVo6BIpvh8U7Srs0SEENeS8o0In+dCMaOgqB4FRSIlpjsQPkn55EyIXIitCY3IB40+kyBXdrB2ZZaIAiJCiDuoyrVWXxArdjT6rJ60v2VCCCG8EONNCfUrsh01n9WjoEiEuPzIKEtECCGGqF8R0aGgSGKkOHyUAiJCpIf6FQkLYwrOkxRQUEQMCDFLRAgRBjnfoEjxhrMhLs8o0k1SQEGRyHB5NpEYU8RyPgkTQghxLQqKiN1clSWigIgQwhV1trYNdbSuR0GRTLjrPWfUbEYIMUdI/YrEmEnnE/UpqkdBkYg4+7UetnBVJ0fKEhEiTPTbJFJGD28kAChLRAgh1nB5iKPQm9Do4Y31KCiSAbGlhelOlBDCJ2VpDao13u6uhqBxbQKj5jPiUs5sOnPGMHzKEvHLIyjQ3VUgRI/PGxch9SviSujZIEsYx07WUgmKKFNEBEXOWSJrgY89gZGc9yMhjlKW1qFa4+nuahA3oqBI4vhqOnPFXZncLuTOzP40Xrbc9i0hxD4MAGPcPi8FFBSJnKuazuxBTWfmuasZTLdeCo4IIaZooYCCw1Op6YnWRDYoS8SdR1CgIPoFCaUeRPzE9pvl6wZQ6q/7kDvKFEmYO0adUZboH0IOPhrWTWwXNyI9ihvlYIH+7q6GrNHos3qUKRIxITadOUKKF2UhB0SNiamuhDQkxCdbi3UEGr3mox5liohFNAzfPmINMKjPESGEUFAkWWJ5YKNULsJiDYYa8wgKlMx3QgixHWMcR59JZPgZBUUiJbSmM7lmiaQSDDVEWSNiK+31G5L8DcgR9SmqR32KiFnOHnUm9ouu1C8GUt8+IiyuuLGiEWjEGsoUSZBYms7ETC4BAzWnESIPlCmqJ+lM0YoVKxATEwO1Wo34+HgcOHDAbNmtW7diyJAhaN68OQICApCYmIhvv/3WoExmZiYUCoXRVFlZ6exN4ZUtd0vO7mAt1gutHJ/zI7ftJeIixHegiRGNPqsn2aAoKysLM2fOxNy5c5Gfn48+ffpg+PDhKCgoMFn++++/x5AhQ5CdnY0jR45gwIABuPfee5Gfn29QLiAgAEVFRQaTWq12xSbpUerWPeQcHMgxGCTyJcdh+bqO1lwmKZBs89nSpUsxefJkTJkyBQCwbNkyfPvtt1i5ciUWL15sVH7ZsmUG/160aBG+/PJLbN++Hd26ddPPVygUCAsLc2rduXB105lcskQUEBBCiPRJMlNUXV2NI0eOIDk52WB+cnIycnNzbVqGVqtFeXk5goKCDOZXVFQgKioKLVu2REpKilEmqaGqqiqUlZUZTO7Gd9OZHFBA9A/aF4RIU322R8FhcvcW8EOSQdHVq1dRV1eH0NBQg/mhoaEoLi62aRlvvfUWbt26hdGjR+vntWvXDpmZmfjqq6+wadMmqNVq9O7dG2fOnDG5jMWLF0Oj0einyMhIxzdKIsSWJaIgwBjtE9KQ2H7TxDRuARG3TtpCIsmgSEehMPySGGNG80zZtGkT5s+fj6ysLISEhOjn9+rVCw8//DC6dOmCPn364NNPP8Wdd96Jd9991+RyMjIyUFpaqp8KCwu5bRCE1Z9I6s8moou/ebRviDOIaVg+kSZJ9ikKDg6Gp6enUVaopKTEKHvUWFZWFiZPnozPPvsMgwcPtljWw8MDd911l9lMkUqlgkqlsq/yHFjrT0RNZ7aji751NFyfyJmqXIsqf+nkFdjfE5fPS4F0vtEGlEol4uPjkZOTYzA/JycHSUlJZj+3adMmTJw4ERs3bsSIESOsrocxhmPHjiE8PJxzneVADBdQGmVlH9pXRAjoZo47aj6rJ8lMEQDMmjULqampSEhIQGJiIlavXo2CggKkp6cDqG/aunjxIj7++GMA9QHRhAkT8Pbbb6NXr176LJOPjw80Gg0AYMGCBejVqxfatm2LsrIyvPPOOzh27BiWL1/ukm2ipjMiRJQxIoRIhWSDojFjxuDatWtYuHAhioqKEBcXh+zsbERFRQEAioqKDJ5Z9P7776O2thZPPPEEnnjiCf38tLQ0ZGZmAgBu3ryJRx99FMXFxdBoNOjWrRu+//579OjRw6XbZgofTWfOJIaLJmU9HEeBESEiR+1nACTafKYzbdo0nDt3DlVVVThy5Aj69u2r/1tmZib27dun//e+ffvAGDOadAERAPzf//0fzp8/j6qqKpSUlODbb79FYmKiC7fIueScgqaAiDvah0RKZPe6JK5NZw42n9nz5gmg/lE3c+fORVRUFFQqFdq0aYN169Y5tG5TJJspIs4jtaYzupgTQvjmXQHU+Lm7Frbj+lRqRz6re/PEihUr0Lt3b7z//vsYPnw4Tp48iVatWpn8zOjRo3H58mWsXbsWd9xxB0pKSlBbW+t4xRuhoEgkhNSfyF5CblahgIhf1IxGCLGVvW+e2LlzJ/bv348///xT/2Dl6OhoXusk6eYzuaCh+ERIKNAkRHz4Gn3W+C0OVVVVJtfnyJsnvvrqKyQkJOD1119HixYtcOedd2L27Nm4fZu/6xcFRcQu9jadCTlrQBdv56F9SxwlteZ50dD1C+IyAYiMjDR4k4OpjA/g2Jsn/vzzT/zwww84ceIEtm3bhmXLluHzzz83GBzFFTWfEVmiizYh/NFev0G/KQIAKCwsREBAgP7f1h5gbM+bJ7RaLRQKBTZs2KB/VM7SpUsxatQoLF++HD4+PhxrT5kiAvk1ndHJ2zVoPxOhcvcjSoRI19GaywQAAQEBBpO5oMiRN0+Eh4ejRYsW+oAIANq3bw/GGC5cuMDLfqCgSOL4/PFLoemMLtSuRfubEJFgPEx2cOTNE71798alS5dQUVGhn/ff//4XHh4eaNmypX0VMIOCIpGT3bM0CCGESMKsWbOwZs0arFu3DqdOncLTTz9t9OaJCRMm6MuPGzcOzZo1w6RJk3Dy5El8//33ePbZZ/HII4/w0nQGUJ8iIiOUtXAPGqZPXMGz7DbqAvi5MNpKSi+F5fr+Mkc+a++bJ/z8/JCTk4Mnn3wSCQkJaNasGUaPHo1XXnnF4Xo3RkGRzNnan0gKTWfEfSgwIkQE3PCqjmnTpmHatGkm/9bwjRI67dq1M2py45M0QlxCrKAsESGEEGsoUyRhNMKiHgVEwkDZIkKEyx3NZ0JEmSIRE2ona7rwEXMoQCVEoFw8+kyoKFNErBLzE2bpIkwIIbZQ/D1x+bz4UaZIxuT20EYiDBSoSg99p0QqKFNEJItO1IQQYiOuTWASaT6jTJFEuauTNfUnIraggJUQgaE+RQAoKBItV3WyFmt/IrroCh99R4QQoaHmM5mi/kSEEEL0mKJ+4vJ5CaBMEZEcykCIB31XhAiDqbfe2ztJAQVFhDfUn4gQQoiYUVAkQXJ+kjVlHsSHvjNCBIA6WgOgPkXEArF2siaEEGIn6lMEgDJFoiTU13u4G2UcxIu+O8JVXYCPu6tAJIAyRTLkjJFn1J+IEHniO6Blgf68Lo/YRsHqJy6flwLKFBFCBIGyRYS4EfUpAkBBEZEIuqBKA32PhLiJrk8Rl0kCKCiSGL5GnlEna0IIIXJDQRHhzN39iSi7IC30fRLiBtR8BoA6WhNCCCGEa2AjkaCIMkUyQ+88I2JA2SLiTNUab3dXgQgUZYpEhp5RZIgunoS4D/3+JIQyRQAoU0QIESi64BKhq/KX0CWURp8BoKCIcOTuTtaEEOmgBzcSd6OgiBgRy3B8yiRIH33HhLiG7onWXCYpoD5FEsLXM4oIIcQaClglhvoUAZB4pmjFihWIiYmBWq1GfHw8Dhw4YLH8/v37ER8fD7VajdatW2PVqlVGZbZs2YIOHTpApVKhQ4cO2LZtm7OqTwgBXXyJdfQyWMIXyQZFWVlZmDlzJubOnYv8/Hz06dMHw4cPR0FBgcnyZ8+exT333IM+ffogPz8fL7zwAmbMmIEtW7boy+Tl5WHMmDFITU3F8ePHkZqaitGjR+Onn35y1WYRiWCB/hYnYogCI0KIKygYYxJJehnq2bMnunfvjpUrV+rntW/fHiNHjsTixYuNyj///PP46quvcOrUKf289PR0HD9+HHl5eQCAMWPGoKysDN98842+zLBhwxAYGIhNmzYZLbOqqgpVVVX6f5eVlSEyMhKlpaUICAiwa3sSx70FwPKQfGvNZ7Y+o8iePkXu6mgttosk10BHLP28nIk69QuHs35/jv5O7M0U2fKcomqNp9Uy1kaf1fgZ/vvXt562usyGysrKoNFoHLpm2LuOqCWvwEOtdng52spKnH/+RafW1RUkmSmqrq7GkSNHkJycbDA/OTkZubm5Jj+Tl5dnVH7o0KE4fPgwampqLJYxt8zFixdDo9Hop8jISEc3iYgUX5kfyiCJLxAm9pH78e12NCQfgESDoqtXr6Kurg6hoaEG80NDQ1FcXGzyM8XFxSbL19bW4urVqxbLmFtmRkYGSktL9VNhYaGjmyRIdOdumTNO8hQcEUKI80h69JlCYRi5MsaM5lkr33i+PctUqVRQqVR21dkSepq1OLgiaGGB/rJsUvMICqRg3M3EnrFz1Ss+GjedCR6NPgMg0aAoODgYnp6eRhmckpISo0yPTlhYmMnyXl5eaNasmcUy5pYpRmK40Ar5pOzKLI5uXWL4zgghAkdBEQCJNp8plUrEx8cjJyfHYH5OTg6SkpJMfiYxMdGo/K5du5CQkABvb2+LZcwtk8iLu5q15NacJuSgmLgeDccnfJJkpggAZs2ahdTUVCQkJCAxMRGrV69GQUEB0tPTAdT397l48SI+/vhjAPUjzd577z3MmjULU6dORV5eHtauXWswquypp55C3759sWTJEtx///348ssvsXv3bvzwww9u2UYiHO4OTOTanEZcS2ijztxJUu89A/enUtMTrQVuzJgxuHbtGhYuXIiioiLExcUhOzsbUVFRAICioiKDZxbFxMQgOzsbTz/9NJYvX46IiAi88847eOCBB/RlkpKSsHnzZrz44ot46aWX0KZNG2RlZaFnz54u3z4iHEI5ocspMKK+RYTwjJrPAEg4KAKAadOmYdq0aSb/lpmZaTSvX79+OHr0qMVljho1CqNGjeKjeoTwjgIjQghxnLTyfzIml/eeCa0/iVCyRA0JsU5E/IT223MmWx7cKDmMh0kCKCgixEFCDj6EXDc+yelCLVVcjlVnPMlarky99d7eSQooKCJEouQSGBHno+CTyAUFRYQ4QCwBh1jqyQVdsAnhAb3mAwAFRbJh68tgbSXnDq5iCzTEVl9HUGDkPLRv7Se6p1kD1KfobxQUESIDcgiMiPiI9biU2jOKAOpTpCO9b5YQJxLrSVwOKKMhL9TJmjgDBUWEyAQFdMRecgw0ZTkcH6Dms79RUEREw90naCkEFVLYBkvcfYxIibP3pdSPRdHh2nRGQREhhAgPBUaEEEdRUESIzNAdOrGGAksZouYzABQUESJLUg+M6KIube7sZG1t5Jkoh+MDFBT9jYIiQmwgxSBCitvUEAVGjnHFfhPqsSfbTtZEj4IiQohkUWBkH9pf8kXPKapHQREhMibUO3YiXVyPOXubzgixBwVFhBBJo+yHbaS6n+ihjcQeFBQRYoXUsylS3z5Auhd8vtD+4YdoO1kD1NH6bxQUEUIoMCIuIeTjzJZO1lJ855kO9Smq5+XuChBCCHEfMQWL1J/IySQS2HAh3bCXEGIXId/F80VMAYAruHJ/uOP4ov5ExF4UFBFCZIUCo3q0H1xL8P2NqE8RAAqKZIPSzsQWcsgWARQQuHr7+TiunHkOk3t/IoD6FOlI+1smhBAz5BoYyXW7nU3wmSBiEwqKJILazp1DLpmThuS0zXILEMS6vY5kiYR0ThRFwETNZwAoKCKEmECBkfS4azvldCyJGTWf1aOgiBAie1IPjKS+fVxQfyLSEH3ThBCT5HaHL9XAwZ3b5a4O1q5uOhNF85g1bmo+W7FiBWJiYqBWqxEfH48DBw7Y9Lkff/wRXl5e6Nq1q2MrNoOCIuIQqV5AiLxJ7biW2vYQJ3JDUJSVlYWZM2di7ty5yM/PR58+fTB8+HAUFBRY/FxpaSkmTJiAQYMG2b9SKygoEhFb0ryE8Elu2SJAOoGEu7dDDMeOK86pksgi2aGsrMxgqqqqMlt26dKlmDx5MqZMmYL27dtj2bJliIyMxMqVKy2u47HHHsO4ceOQmJjId/UpKCKEkMY8ggLdHlQ4Ssx1b0wITWdy6U/EV0fryMhIaDQa/bR48WKT66uursaRI0eQnJxsMD85ORm5ublm6/nhhx/if//7H+bNm8fbtjdE7z4jhFjEAv2huFHu7mq4hUdQILTXb7i7GjYTSjAkhiwRXySTCeI6rP7vzxYWFiIgIEA/W6VSmSx+9epV1NXVITQ01GB+aGgoiouLTX7mzJkzmDNnDg4cOAAvL+eELxQUEUKIBWIIjIQSDPGJnsLvYjwFRQEBAQZBkTUKhcJwMYwZzQOAuro6jBs3DgsWLMCdd97JoaKWUVBECLFKztki4J+gQ4jBkdACIndmiexpOqP+RO4VHBwMT09Po6xQSUmJUfYIAMrLy3H48GHk5+dj+vTpAACtVgvGGLy8vLBr1y4MHDiQc73k0VhKbCantDch9hJSACLEvkNSO3/IpT8R4PqHNyqVSsTHxyMnJ8dgfk5ODpKSkozKBwQE4Ndff8WxY8f0U3p6OmJjY3Hs2DH07NmTy+brUaaIiIb2+g3BXQTkRO7ZIp2Gx6A7Mkdy+A2IpelMUpkgnprP7DFr1iykpqYiISEBiYmJWL16NQoKCpCeng4AyMjIwMWLF/Hxxx/Dw8MDcXFxBp8PCQmBWq02ms+FJMPgGzduIDU1Vd/7PTU1FTdv3jRbvqamBs8//zw6deqEJk2aICIiAhMmTMClS5cMyvXv3x8KhcJgGjt2rJO3hj9iOdEQIhauytbo1iPkgMjdWSKhNZ0R68aMGYNly5Zh4cKF6Nq1K77//ntkZ2cjKioKAFBUVGT1mUV8UzDGJPLGkn8MHz4cFy5cwOrVqwEAjz76KKKjo7F9+3aT5UtLSzFq1ChMnToVXbp0wY0bNzBz5kzU1tbi8OHD+nL9+/fHnXfeiYULF+rn+fj4QKPR2FSvsrIyaDQalJaW2tURDQASx70FAFCW1lkspyytsfh3z7LbVtdlazZADnfJ7j7RCxFliyzj63ch5ADIFL5+K47evPEdFNnSdGYpU2Tqb78veNrqMhvics2wdx3tpy+Cp0rt8HLqqipx6r0XnFpXV5Bc89mpU6ewc+dOHDx4UN/G+MEHHyAxMRGnT59GbGys0Wc0Go1Ru+a7776LHj16oKCgAK1atdLP9/X1RVhYmHM3ghABo2Y0y0wFM7YESmILghpy982Dq1/rIUluaD4TIsk1n+Xl5UGj0Rh0uurVqxc0Go3FB0I1VlpaCoVCgaZNmxrM37BhA4KDg9GxY0fMnj0b5eXmLw5VVVVGT/eUEjGfxAlxpYbNX+YmseIzIHJFEz9fTWeS6k9E9CSXKSouLkZISIjR/JCQELMPhGqssrISc+bMwbhx4wzSgOPHj0dMTAzCwsJw4sQJZGRk4Pjx40ZZJp3FixdjwYIFjm2IGdUaT6tNaIQ/ihvlbr8LFiLKFhEiMZQpAiCiTNH8+fONOjk3nnT9f0w9+MncA6Eaq6mpwdixY6HVarFixQqDv02dOhWDBw9GXFwcxo4di88//xy7d+/G0aNHTS4rIyMDpaWl+qmwsNCBLSeEEGESQpbIGU1nXIfiizGLpOBhkgLRZIqmT59udaRXdHQ0fvnlF1y+fNnob1euXDH5QKiGampqMHr0aJw9exbfffed1c5i3bt3h7e3N86cOYPu3bsb/V2lUpl9xDkhYkfZIiI21HRGrBFNUBQcHIzg4GCr5RITE1FaWoqff/4ZPXr0AAD89NNPKC0tNflAKB1dQHTmzBns3bsXzZo1s7qu3377DTU1NQgPD7d9Q0SALnaEEGukmiWSLWo+AyCi5jNbtW/fHsOGDcPUqVNx8OBBHDx4EFOnTkVKSorByLN27dph27ZtAIDa2lqMGjUKhw8fxoYNG1BXV4fi4mIUFxejuroaAPC///0PCxcuxOHDh3Hu3DlkZ2fjwQcfRLdu3dC7d2+3bKsjxP6sIiG+ZkHOqL+VPEn5e5dj0xng+idaC5XkgiKgfoRYp06dkJycjOTkZHTu3BmffPKJQZnTp0+jtLQUAHDhwgV89dVXuHDhArp27Yrw8HD9pBuxplQqsWfPHgwdOhSxsbGYMWMGkpOTsXv3bnh6CudBYHTnRAhxJr4DIlfdqFHTmRWMh0kCRNN8Zo+goCCsX7/eYpmGz6yMjo6GtWdYRkZGYv/+/bzUT0rE8AZx4lzU3ErcgW4AiTNIMlMkdfSIeteiCz4h9aSeJZLTC2BNknmWCKCgiBDCAyn3MSH1hPQdCzVLJOamNepTVI+CImKSkE6ARBzomJEuZ3y3Qhv0wfVdZ9bU+kkkapA4CopkSGgnI3tRHyZCXEdowa69WSLqbmAj6mgNgIIiSRJqaplIn9AuoESYxHhjZi1LZOnvYsgSUfNZPQqKCGdifpmlraizNZEjoQW5zsoSyb6DNdGjI4EQwiuhXUiJY5z1PYoxSyQL1HwGgIIi0ZJ7Ozn1KxI2CozETYgBkTv7Ekm96Qyg5jMdCoqIWXRhI0R+5Pa7p6Yz0hAdDTJFKWz7Ub8i+8jt4ioFzvzOXJklIg6g5jMAFBRJlqtPInLobE3sR4GReAg1IHIEnx2s5dB0BoCCor9RUEREi/oVEcIPIQevlCVyDepTVI+CIhGTe2drd6AmNPsJ+YJLnP/9yDVLRMSJgiJiEV3QCB/oOBImoQdEYs4SiarpDKDms79RUCRjfN/BUb8iYgkFRsIixe+DsueOUzDGeZICCookTMx3WbZyR78iakIjYueKgEjIWSJqOiPmUFBECHEZKWYnxEYMAZEjhJQlEl3TGUDNZ3+joEjkXHEioAsZ4RMdT+7BAv1Fs+/dnSWSIxp9Vo+ODsIrd/QroiY08RHLxVkqXLm/pZAloqYz+aKgSOboydaESJvYAiKxZ4lE2XQGUPPZ3ygokjg5dLZ2F8oWcUPZIueTQ0BEWSJ+UPNZPQqKiE2EfgGjp1uLk9CPK7ESU/8hQoSEgiKRsJT2FdKoC4CeV0TsQxdvfrljf4ohS8THMHxrLDWd1fppuS3c2aj5DAAFRQTUr4gLakLjBwVG/BBrQCQWUm06A6j5TIeCIhmgfkVEDCgwcpy7msv4CoiEkCXiStRZIoAyRX+joIjYzJ6TrlyG5gOULeITBUb2E/s+E8pNm1w7WBNDXu6uALFdlb8HVOWm7ziqNZ5Qlta5uEaE8E93kadg0zJ3B0PubDYTWpZIKqTSBMYFBUUEQP0JzrPstrurQYgeC/SnwMgEdwdDgHubzZxB1h2sdRirn7h8XgIohJYJvk4+QjghW0JNaNIi9OPNlYQyzN7dAZE7skTUdCYfFBQRp6Gh+YQPQgkG3Eko2+/ukWbOePwIZYnq0eizehQUiYwzn1fk7hMeXyhbJE1CCQxcRRcMCmW7+Tw/uKLZjLJEdqLRZwCoTxEhRETk0AlbKEFQQ0IIiMSWJSLiRJkiGXFHvyK5NaFJ+WItJELKoPBFitvkLkLIEomp6QwAFFrukxRQpohIkvb6DdkFZHIk9syRGIIgsWWJbA2IZNMsZiuuTWASSZpRpkiEqF+RsIn1Ai1mYsuyiKW+YguI+GQtaJJKB2tiSJJB0Y0bN5CamgqNRgONRoPU1FTcvHnT4mcmTpwIhUJhMPXq1cugTFVVFZ588kkEBwejSZMmuO+++3DhwgUnbgn/5NSE5q4O1wAFRu4itM7JOg3rJbS6mSOEgMhelCVyHI0+qyfJoGjcuHE4duwYdu7ciZ07d+LYsWNITU21+rlhw4ahqKhIP2VnZxv8febMmdi2bRs2b96MH374ARUVFUhJSUFdHT1JmhChcWcQIsYgqCGhZIwpS+RCuoc3cpkkQHJ9ik6dOoWdO3fi4MGD6NmzJwDggw8+QGJiIk6fPo3Y2Fizn1WpVAgLCzP5t9LSUqxduxaffPIJBg8eDABYv349IiMjsXv3bgwdOpT/jbHAma/8kNLTrd3Zt0hxo1yUF0QpMvU98JXNk9p3zHdA5KpmM8oSccM12yOVTJHkgqK8vDxoNBp9QAQAvXr1gkajQW5ursWgaN++fQgJCUHTpk3Rr18/vPrqqwgJCQEAHDlyBDU1NUhOTtaXj4iIQFxcHHJzc00GRVVVVaiqqtL/u6ysjI9N5Kxa4w1laQ3n5djzGgaPoEC3NmcR0pjUghk+CCUgcicuWSIifpJrPisuLtYHMg2FhISguLjY7OeGDx+ODRs24LvvvsNbb72FQ4cOYeDAgfqgpri4GEqlEoGBhlmH0NBQs8tdvHixvl+TRqNBZGQkhy0jYkR9i4hYCCkgEmuWyGLTmV+tc1fOFT28EYCIgqL58+cbdYRuPB0+fBgAoFAojD7PGDM5X2fMmDEYMWIE4uLicO+99+Kbb77Bf//7X+zYscNivSwtNyMjA6WlpfqpsLDQji12L6H0KeADZagIsUzMv3e+nknkVEIPiEAdrXVE03w2ffp0jB071mKZ6Oho/PLLL7h8+bLR365cuYLQ0FCb1xceHo6oqCicOXMGABAWFobq6mrcuHHDIFtUUlKCpKQkk8tQqVRQqVQ2r9NezuxXZCtqQrOO+hYRIXNGQOTKLJGtbMkSybKDNTEgmqAoODgYwcHBVsslJiaitLQUP//8M3r06AEA+Omnn1BaWmo2eDHl2rVrKCwsRHh4OAAgPj4e3t7eyMnJwejRowEARUVFOHHiBF5//XUHtsi9+OpXJBb0MEdCjIk9IBJFlkgsuI4gk8joM8kdUe3bt8ewYcMwdepUHDx4EAcPHsTUqVORkpJi0Mm6Xbt22LZtGwCgoqICs2fPRl5eHs6dO4d9+/bh3nvvRXBwMP71r38BADQaDSZPnoxnnnkGe/bsQX5+Ph5++GF06tRJPxpNapxxwpRrYEJ9i4iQ1AX4CC4gcia3ZolE0HQGUPOZjuSCIgDYsGEDOnXqhOTkZCQnJ6Nz58745JNPDMqcPn0apaWlAABPT0/8+uuvuP/++3HnnXciLS0Nd955J/Ly8uDv/0+zx//93/9h5MiRGD16NHr37g1fX19s374dnp7ueZYG4NynW9tKLE1Dcmy6I6QxZ/Uf4hoQubNzNQ3TJzqiaT6zR1BQENavX2+xDGuQ6vPx8cG3335rdblqtRrvvvsu3n33Xc51FAK5NaG5G/UtIu4m1A7VQm82k3qWCAC9++xvkswUSZG77mSoCY1f1IxG3MWZAZEQO1YDlCWyBzWf1aOgSOKoCc0QNaERORJqQOQIyhIRZ6KgSAK4nCSE2jFSyihbRFxJyAERZYkERMu4TxJAQZGICL0JzZ5skTub0ISQLaLAiDibs0aY6bgjIBJSlsgiMWaJ6InWACTa0ZoYctWDHAkhwuDsDtXuyDDbExC5IksktYc1KsDxhbC81cS9KFMkEWJsQqNsEWWLCP/EEBC5u9nMFrLKEhE9CopERkpNaIQQfgl1yH1DQmg24/Kgxvq/SytLBOCfJ1pzmSSAgiKZcNUoNHtRtoiyRYQ7Z/cf0hF6PyIhd5z28atydxUsoiH59SgokhBnN6GJ4S6UELlx1e9SKiNVnZoloqYzu61YsQIxMTFQq9WIj4/HgQMHzJbdunUrhgwZgubNmyMgIACJiYk2PXjZHhQUiZCQ74YAcTWhUbaIiJmYAiLKEgmcG0afZWVlYebMmZg7dy7y8/PRp08fDB8+HAUFBSbLf//99xgyZAiys7Nx5MgRDBgwAPfeey/y8/PtX7kZFBRJDJd3ocmxwzVAgRERH1c1lwHSCogoS2SegjHOEwCUlZUZTFVV5gPCpUuXYvLkyZgyZQrat2+PZcuWITIyEitXrjRZftmyZXjuuedw1113oW3btli0aBHatm2L7du387YfKCgidqEO14S4lyubsYU+0swe7sokiSJLxKPIyEhoNBr9tHjxYpPlqqurceTIESQnJxvMT05ORm5urk3r0mq1KC8vR1BQEOd669BzikSqxg/wrnB3LfjjERTo1oyN9voNt2es6IWxxBJX9+lzZ0DkjmYzOWeJAADavycunwdQWFiIgIAA/WyVSmWy+NWrV1FXV4fQ0FCD+aGhoSguLrZplW+99RZu3bqF0aNHO1ZnEygokqAqfw+oyk0f3dYe5Fit8YaytMbi8usCfOBZdttqPVigPzUL2YkCI2KKGAMiR7mj2cwaR4fgiylL1LAJzNHPA0BAQIBBUGT1cwrDxz4yxozmmbJp0ybMnz8fX375JUJCQuyrrAXUfEYEw92ZGiH0LSKkIVf2HdLhKyASwvOIANsCIocf1AhII0vkBsHBwfD09DTKCpWUlBhljxrLysrC5MmT8emnn2Lw4MG81ouCIhFz9O5HqB2uST3KrhHAPY/AcGdAZC9X9RGSQ5YIgMtHnymVSsTHxyMnJ8dgfk5ODpKSksx+btOmTZg4cSI2btyIESNG2LdSG1DzmURZakLjg7Oa0KhvEZE7dz0PzN0BEWWJ3IzrU6kd+OysWbOQmpqKhIQEJCYmYvXq1SgoKEB6ejoAICMjAxcvXsTHH38MoD4gmjBhAt5++2306tVLn2Xy8fGBRqNxvO4NUFBEiABR3yJ5ooDINoLoXG2B6LJE4P5Uakc+O2bMGFy7dg0LFy5EUVER4uLikJ2djaioKABAUVGRwTOL3n//fdTW1uKJJ57AE088oZ+flpaGzMxMxyvfAAVFIufoKDTqcG2eULJFFBjJhzufFu/ugMherupcbZWUskRuNG3aNEybNs3k3xoHOvv27XN6fahPkUg48gN3RjraFYQQkBDiKlIIiLhw12gzZw3BF2OWCAC9EPZv4rxqEl64ssO12DIeQhmJJqYMG7GPO0aWNcTn71tI/YiEqmkT61l1d1JouU9SIJ8jUsIs3RU5+6TjrJM6ZYuIVLk7GALEGRCJOUsk9ICI/IOCIhER8gsPbUHZIsdQtkg63B0MAcIIiOwl9nOfKFDzGQAKimSPjyY0KWeLKDAifBBCdggQTkDkrAw2ZYk4cPFzioSKgiKRMfejd2cTGiHENKEEQ4B4AyLBNJsRWaCrJbGKz2yRvU1olC36B2WLxENowZDUAyKXkHKWCP+8+4zLJAUUFEmIo9kiV/ULIEQOhBIMAcIYcq/jzIw1ZYl4QH2KAFBQJEruuHuibBFli4hlQsoOAc4JiFx5A+XKZjOrJJ4lIv+gJ1pLjLOecE0IMU1IgZCO0AIioTebyel1HmYxAFySZdJIFFGmSKyE+oRryhY5H2WLhEFomSEdOQVELmk2c+B1HmLMElGfonoUFBE9Vz7hWqwoMCJCDYYA8QdE9nBJNkmKr/Mwh4FjnyJ3bwA/KCiSILEMzxdjtojIm1CDIUAaARHfgY47OleLMUtE/iGcKySxmzPulIT8MEehoGyR/Ag5OwQILyByhKubzaziOUsU1kTgv1cafQaAgiLJEtTzPSygbBERMjEEQ0IMiNzdj8gWrswSCT4gAuo7WXOdJICCIhni2oRG2SLKFkmd0IMhwHl9/FwdEDmDMztXS64vETHg/qPXCW7cuIHU1FRoNBpoNBqkpqbi5s2bFj+jUChMTm+88Ya+TP/+/Y3+PnbsWCdvTT1zP3J3NKHxTazZIqEERoQ/YgiGAGkFRGJqNrNE1Fki0OgzHUkGRePGjcOxY8ewc+dO7Ny5E8eOHUNqaqrFzxQVFRlM69atg0KhwAMPPGBQburUqQbl3n//fWduCifO7HDt7myRUAIjIaBsET/EEAwBFBBxxaXZTNJZIupTBECCD288deoUdu7ciYMHD6Jnz54AgA8++ACJiYk4ffo0YmNjTX4uLCzM4N9ffvklBgwYgNatWxvM9/X1NSrrKrV+DF4VCqP5jj6w0RJXP8yRBfqL8uKuvX6DAjSRE0swBFBAZI27OleLPUtE/iG5TFFeXh40Go0+IAKAXr16QaPRIDc316ZlXL58GTt27MDkyZON/rZhwwYEBwejY8eOmD17NsrLzR/0VVVVKCsrM5hczd3D8ylb5BpiDCjdTSxNZYDzOlTXL1sa7z60JSCiIfgWUKYIgAQzRcXFxQgJCTGaHxISguLiYpuW8dFHH8Hf3x///ve/DeaPHz8eMTExCAsLw4kTJ5CRkYHjx48jJyfH5HIWL16MBQsW2L8RDnBHtqha4w1laQ1v66NsEXEVsQRDgHMfmspHQCSULJE1ru5cLbosEdfARiJBkWgyRfPnzzfbGVo3HT58GEB9p+nGGGMm55uybt06jB8/Hmq12mD+1KlTMXjwYMTFxWHs2LH4/PPPsXv3bhw9etTkcjIyMlBaWqqfCgsL7dxqY5zTw24gh2yREDpdizGgdDUxZYcACoiE0GxmiWSyRERPNJmi6dOnWx3pFR0djV9++QWXL182+tuVK1cQGhpqdT0HDhzA6dOnkZWVZbVs9+7d4e3tjTNnzqB79+5Gf1epVFCpVFaX42yWskhV/h5QlZu/g6JsEZEKMQVDAAVEfAZEru5cLbosEVD/nCHb8gbmPy8BogmKgoODERwcbLVcYmIiSktL8fPPP6NHjx4AgJ9++gmlpaVISkqy+vm1a9ciPj4eXbp0sVr2t99+Q01NDcLDw61vAI9c2eGaL3UBPvAsc85dlUdQoCAyNdSMJkwUDDVevvADIj45q9lMap2ruQ6rpyH5AtW+fXsMGzYMU6dOxcGDB3Hw4EFMnToVKSkpBiPP2rVrh23bthl8tqysDJ999hmmTJlitNz//e9/WLhwIQ4fPoxz584hOzsbDz74ILp164bevXs7fbu4cvfwfHvY+9wi8g/KshmigKjx8t0TENnLZc1mLtTS56a7q2AZdbQGIMGgCKgfIdapUyckJycjOTkZnTt3xieffGJQ5vTp0ygtLTWYt3nzZjDG8NBDDxktU6lUYs+ePRg6dChiY2MxY8YMJCcnY/fu3fD0dP3oDTE+zJH6FhFXooCo8fLdFxAJttnMRVkiwQdERE80zWf2CAoKwvr16y2WYSai2kcffRSPPvqoyfKRkZHYv38/L/VzFy59i6yhvkX/cHczmuJGuayzbWILhgAKiOzl7IDIEsl2rtYyQMEh26OlTBERKMoWEbkSW0DkzOcP/bMO8QREQnmRNZ+dq0WTJaLmMwAUFImaI+3pUu9bJJTAyN3NaGLNsnEhxoDI+etw34MZnRUQCbHZjEgHBUUSJZVskZybgYjtKCAytQ5+fq9i61jtrGYzS0SfJQIAcM0SUaaIuJC5H7rUs0WOoGyRvFBAZGod7g2I3PHEaj7IunM1NZ8BoKCI2ImyReIh9SY0sT2dGqCAiAvqXE1cgYIiEbE3W+ToyYmyRfyhbJFziDEYooDI8fKuCIhk2bm6IS3jPkkABUUyxeVujrJF9nFnYCTFbJEYAyLXrEe+ARFX1LkaANNynySAgiKRoWyR7YSSLSLyRQERP+WtEUPn6iifq7zXgfCPgiIZk0O2SCiBETWj8UNMWSK5BETO5O5mM76yRKIIiKijNQAKikTJ3jc+Cz1bJKYLnRhJpQlNTMeJnAIiQfcjssLRgEhSfYl0qE8RAAqKJIXv4fnWuOPBcJQtkicKiEytR94BkU2o2cx2lCkCQEGRaMk5WyTmTtdEulw1wqx+XdINiGwlhmYzIj4UFEmMELNF1On6H+7KFom5CU0MWSJXHuNSD4jc3WxmiWSzRED9A6k5ZYrcvQH8oKBIRtyVLbKVXLJF1IxmOwqIGq9LfAGRPYTQbGbvk6vNEVVABFDz2d8oKBIxV736wxrKFhFnoICo8brEGRDx3Y9IaM1mou5cTYxQUCQWPHUYlHu2SCiBkTuyRWJuQhMiMQZEXDgzILIVNZs5kVbLfZIACopEjrJFRIqEniVyZYdqPgMiVz6LyC39iJzQbCYb1HwGgIIicRF4toge6Ggf6lskTmIbYaYj5o7VNnFSsxlfWaLWyhLzlSOCQUGRBAglW2QLZ1xQxNzp2tXE0IQm5CwRBUT84q0fkRV8N5vZSxQBEWWKAFBQJD6ULeIFZYuIvSggct5nLBFis5mk+hLp0BOtAVBQJBlSzBbJpRmN/EOoWSIKiPj/jKv6EVGzGbEHBUViJJNskVxQtkjYKCDi/zNC6Edkiayazf7GmJbzJAUUFEkIZYvEmy1yZWAk1H5FQswSUUDE/2dc1Y/IGmo2a4RxbDqjPkVECihbRIhpFBDx/xleAyIXN5vZS0xZIgDU0fpvFBSJlZkTAmWLKFskRkLLEskxIBIKZwdEjpJ0lojoif8XRDhzZ7ZIKEP0hRIYuYpQm9CEQK4BkVCyRM5GnavNoCdaA6CgSNwkkC2yldAyCc4ix2yRkL5bCoic8xlqNgOiva/xsnynoeYzABQUiYYz0sENiSFbJJdmNOIeFBA55zNCCYgcxUezmeADIqJHQZHYyShbJBdyzBa5GwVEzvmMkAIiajazjGm1nCcpoKBIRNyZLbL0N8oWiZMQ+hUJoemMAiLnfkYIXNFsZo5oskTUfAaAgiJpcEG2iCt3DdEX63vRKFvkGhQQOe8zQsoSOYKPLJFoAiKiR0GRyIg1W2QLZ2SLHCGnbBERH3cERPYSUkDEZ7OZOWJuNtOjd58BoKBIOgSQLaIh+vySerbI3U1nYswSuSsgEmo/Imv4bjaTdOdqxgCm5TBRUETcRKjZIldy9wWViJsYAyJ3EfKziBwdfm+JnDpXE2MUFEmJzLJF1OmaOyF0tnY1sQZEYuhHZCt39yNydrOZ6LJEAJiWcZ6kQJJB0auvvoqkpCT4+vqiadOmNn2GMYb58+cjIiICPj4+6N+/P3777TeDMlVVVXjyyScRHByMJk2a4L777sOFCxecsAXWUbZIHqTahCb1TJ8cAyKx9yOyhI9ms9ZeAg8aODWd/T05YMWKFYiJiYFarUZ8fDwOHDhgsfz+/fsRHx8PtVqN1q1bY9WqVQ6t1xxJBkXV1dV48MEH8fjjj9v8mddffx1Lly7Fe++9h0OHDiEsLAxDhgxBefk/P6CZM2di27Zt2Lx5M3744QdUVFQgJSUFdXV1ztgMx1C2yCzKFhFXZIkoILJUTrj9iJzZbCb4gAjuyRRlZWVh5syZmDt3LvLz89GnTx8MHz4cBQUFJsufPXsW99xzD/r06YP8/Hy88MILmDFjBrZs2cJ18/UUjEmkd5QJmZmZmDlzJm7evGmxHGMMERERmDlzJp5//nkA9Vmh0NBQLFmyBI899hhKS0vRvHlzfPLJJxgzZgwA4NKlS4iMjER2djaGDh1qtT5lZWXQaDQoLS1FQECAXdvSftsCk/NvV6iMZ1Z4mSzrVWH65OpVoTC7Xu8K83Wy9DdVueWTn7LUeiCpLK2xWgYAPMvs6zvgSJORuzM2zgzO3PHYAndkiiggcs5neA2IAMGNNuM6BL9hQNQ0otDsekzhcs2wdx39Ff+Cl8Lx30gtq8E+ts2uuvbs2RPdu3fHypUr9fPat2+PkSNHYvHixUbln3/+eXz11Vc4deqUfl56ejqOHz+OvLw8h+vekOmrp8ycPXsWxcXFSE5O1s9TqVTo168fcnNz8dhjj+HIkSOoqakxKBMREYG4uDjk5uaaDIqqqqpQVfXPD7i0tBRA/UFoLz/cROlfaqP52tsmTkieMBkYVXsCXreMT7LVXoDnLdOBkYeF849Htfm//aUClBXmT4K1vtYDo1pfQFlmPTCq9VXAs7zSajkdhdb+pkcts7CxLuDhQJ1tdq0KrKlrA6O6WvOBuDNUB3gDNc7N6FZrPIEa2wJ5m5bn5wFwOOxqmgBw4LCps+O6WNeEATb+9LSeNmaJLNzjqJtUoe4v83/X+Fai9pbpv4U2qUCNmb9V15n+3iJ9rqHSxM1ftPcV/GVi31Z4G29jWYOgyMPOc7/uWuGK3EUtq3K4CQwAalG/Dxtf31QqFVQq45v36upqHDlyBHPmzDGYn5ycjNzcXJPryMvLM7gGA8DQoUOxdu1a1NTUwNub+40PBUUAiouLAQChoaEG80NDQ3H+/Hl9GaVSicDAQKMyus83tnjxYixYYJzhiYyM5KPaRE4cG0hju0tOXj4hBIDGoU+Vl5dDo3Hss9YolUqEhYXhh+Jszsvy8/Mzur7NmzcP8+fPNyp79epV1NXVmbzumrumFhcXmyxfW1uLq1evIjw8nNsGQERB0fz5800GGA0dOnQICQkJDq9DoTC8e2WMGc1rzFKZjIwMzJo1S/9vrVaL69evo1mzZlaX21BZWRkiIyNRWFjotBSqM1C9XU+sdad6uxbV2/UcqTtjDOXl5YiIiHBavdRqNc6ePYvqau7ZcFPXQ1NZoobsve6aKm9qvqNEExRNnz4dY8eOtVgmOjraoWWHhYUBqI9CG0aaJSUl+qg0LCwM1dXVuHHjhkG2qKSkBElJSSaXayptaOtoOFMCAgJEdyIAqN7uINa6U71di+rtevbW3VkZoobUajXUauPuGc4UHBwMT09Po6xQw+tuY2FhYSbLe3l5oVmzZrzUSzSjz4KDg9GuXTuLk6NfakxMDMLCwpCTk6OfV11djf379+sDnvj4eHh7exuUKSoqwokTJ8wGRYQQQggxplQqER8fb3BNBYCcnByz19TExESj8rt27UJCQgIv/YkAEQVF9igoKMCxY8dQUFCAuro6HDt2DMeOHUNFxT895tq1a4dt27YBqE+7zZw5E4sWLcK2bdtw4sQJTJw4Eb6+vhg3bhyA+mh98uTJeOaZZ7Bnzx7k5+fj4YcfRqdOnTB48GC3bCchhBAiVrNmzcKaNWuwbt06nDp1Ck8//TQKCgqQnp4OoL4LyoQJE/Tl09PTcf78ecyaNQunTp3CunXrsHbtWsyePZu/SjEJSktLYwCMpr179+rLAGAffvih/t9arZbNmzePhYWFMZVKxfr27ct+/fVXg+Xevn2bTZ8+nQUFBTEfHx+WkpLCCgoKnL49lZWVbN68eayystLp6+IT1dv1xFp3qrdrUb1dT8x1d6bly5ezqKgoplQqWffu3dn+/fv1f0tLS2P9+vUzKL9v3z7WrVs3plQqWXR0NFu5ciWv9ZH0c4oIIYQQQmwlyeYzQgghhBB7UVBECCGEEAIKigghhBBCAFBQRAghhBACgIIiQgghhBAAFBQJwquvvoqkpCT4+vra/MRrxhjmz5+PiIgI+Pj4oH///vjtt98MylRVVeHJJ59EcHAwmjRpgvvuuw8XLlzgrd43btxAamoqNBoNNBoNUlNTcfPmTYufUSgUJqc33nhDX6Z///5Gf7f2NHNX1H3ixIlG9erVq5dBGaHt85qaGjz//PPo1KkTmjRpgoiICEyYMAGXLhm+7Izvfb5ixQrExMRArVYjPj4eBw4csFh+//79iI+Ph1qtRuvWrbFq1SqjMlu2bEGHDh2gUqnQoUMH/XPG+GRPvbdu3YohQ4agefPmCAgIQGJiIr799luDMpmZmSaP98pK219g7Iy679u3z2S9fv/9d4NyQtvnpn6DCoUCHTt21JdxxT7//vvvce+99yIiIgIKhQJffPGF1c8I5RgnVvA6wJ845OWXX2ZLly5ls2bNYhqNxqbPvPbaa8zf359t2bKF/frrr2zMmDEsPDyclZWV6cukp6ezFi1asJycHHb06FE2YMAA1qVLF1ZbW8tLvYcNG8bi4uJYbm4uy83NZXFxcSwlJcXiZ4qKigymdevWMYVCwf73v//py/Tr149NnTrVoNzNmzd5qTOXuqelpbFhw4YZ1OvatWsGZYS2z2/evMkGDx7MsrKy2O+//87y8vJYz549WXx8vEE5Pvf55s2bmbe3N/vggw/YyZMn2VNPPcWaNGnCzp8/b7L8n3/+yXx9fdlTTz3FTp48yT744APm7e3NPv/8c32Z3Nxc5unpyRYtWsROnTrFFi1axLy8vNjBgwcdqiMf9X7qqafYkiVL2M8//8z++9//soyMDObt7c2OHj2qL/Phhx+ygIAAo+Oeb/bWfe/evQwAO336tEG9Gh6nQtznN2/eNKhvYWEhCwoKYvPmzdOXccU+z87OZnPnzmVbtmxhANi2bdsslhfKMU6so6BIQD788EObgiKtVsvCwsLYa6+9pp9XWVnJNBoNW7VqFWOs/uTh7e3NNm/erC9z8eJF5uHhwXbu3Mm5ridPnmQADH6weXl5DAD7/fffbV7O/fffzwYOHGgwr1+/fuypp57iXEdzHK17Wloau//++83+XSz7/Oeff2YADC48fO7zHj16sPT0dIN57dq1Y3PmzDFZ/rnnnmPt2rUzmPfYY4+xXr166f89evRoNmzYMIMyQ4cOZWPHjuWlzozZX29TOnTowBYsWKD/t62/aa7srbsuKLpx44bZZYphn2/bto0pFAp27tw5/TxX7XMdW4IioRzjxDpqPhOhs2fPori4GMnJyfp5KpUK/fr1Q25uLgDgyJEjqKmpMSgTERGBuLg4fRku8vLyoNFo0LNnT/28Xr16QaPR2Lz8y5cvY8eOHZg8ebLR3zZs2IDg4GB07NgRs2fPRnl5Oec681H3ffv2ISQkBHfeeSemTp2KkpIS/d/EsM8BoLS0FAqFwqiplo99Xl1djSNHjhjsAwBITk42W8e8vDyj8kOHDsXhw4dRU1NjsQwf+9XRejem1WpRXl6OoKAgg/kVFRWIiopCy5YtkZKSgvz8fF7qrMOl7t26dUN4eDgGDRqEvXv3GvxNDPt87dq1GDx4MKKiogzmO3uf20sIxzixjZe7K0Dsp3tLcOM3CYeGhuL8+fP6MkqlEoGBgUZlGr9l2NE6hISEGM0PCQmxefkfffQR/P398e9//9tg/vjx4/Uv6T1x4gQyMjJw/PhxoxcBurruw4cPx4MPPoioqCicPXsWL730EgYOHIgjR45ApVKJYp9XVlZizpw5GDdunMGbuvna51evXkVdXZ3JY9NcHYuLi02Wr62txdWrVxEeHm62DB/71dF6N/bWW2/h1q1bGD16tH5eu3btkJmZiU6dOqGsrAxvv/02evfujePHj6Nt27Zuq3t4eDhWr16N+Ph4VFVV4ZNPPsGgQYOwb98+9O3bF4D570Uo+7yoqAjffPMNNm7caDDfFfvcXkI4xoltKChykvnz52PBggUWyxw6dAgJCQkOr0OhUBj8mzFmNK8xa2Vsrbep9dtaB51169Zh/PjxUKvVBvOnTp2q//+4uDi0bdsWCQkJOHr0KLp37+62uo8ZM8agXgkJCYiKisKOHTuMAjt7luuqfV5TU4OxY8dCq9VixYoVBn9zdJ+bY++xaap84/mOHO/2cnQdmzZtwvz58/Hll18aBK69evUy6Izfu3dvdO/eHe+++y7eeecd/ioO++oeGxuL2NhY/b8TExNRWFiIN998Ux8U2btMRzm6jszMTDRt2hQjR440mO/KfW4PoRzjxDIKipxk+vTpVkfvREdHO7TssLAwAPV3H+Hh4fr5JSUl+juNsLAwVFdX48aNGwaZi5KSEiQlJXGu9y+//ILLly8b/e3KlStGdzumHDhwAKdPn0ZWVpbVst27d4e3tzfOnDlj8QLtqrrrhIeHIyoqCmfOnAEg7H1eU1OD0aNH4+zZs/juu+8MskSm2LrPGwsODoanp6fR3W3DY7OxsLAwk+W9vLzQrFkzi2Xs+b74rrdOVlYWJk+ejM8++wyDBw+2WNbDwwN33XWX/pjhA5e6N9SrVy+sX79e/28h73PGGNatW4fU1FQolUqLZZ2xz+0lhGOc2Mj13ZiIOfZ2tF6yZIl+XlVVlcmO1llZWfoyly5d4r3T708//aSfd/DgQZs7/aalpRmNgDLn119/ZQAM3p7MBde661y9epWpVCr20UcfMcaEu8+rq6vZyJEjWceOHVlJSYlN6+Kyz3v06MEef/xxg3nt27e32NG6ffv2BvPS09ONOqEOHz7coMywYcN47/RrT70ZY2zjxo1MrVZb7Wiro9VqWUJCAps0aRKXqhpxpO6NPfDAA2zAgAH6fwt1nzP2T0fxX3/91eo6nLXPdWBjR2shHOPEOgqKBOD8+fMsPz+fLViwgPn5+bH8/HyWn5/PysvL9WViY2PZ1q1b9f9+7bXXmEajYVu3bmW//vore+ihh0wOyW/ZsiXbvXs3O3r0KBs4cCDvw8M7d+7M8vLyWF5eHuvUqZPR8PDG9WaMsdLSUubr68tWrlxptMw//viDLViwgB06dIidPXuW7dixg7Vr145169aNt3o7Uvfy8nL2zDPPsNzcXHb27Fm2d+9elpiYyFq0aCHofV5TU8Puu+8+1rJlS3bs2DGDIcpVVVWMMf73uW6Y9dq1a9nJkyfZzJkzWZMmTfQjhObMmcNSU1P15XXDlZ9++ml28uRJtnbtWqPhyj/++CPz9PRkr732Gjt16hR77bXXnDY83NZ6b9y4kXl5ebHly5ebfZTB/Pnz2c6dO9n//vc/lp+fzyZNmsS8vLwMAlt31P3//u//2LZt29h///tfduLECTZnzhwGgG3ZskVfRoj7XOfhhx9mPXv2NLlMV+zz8vJy/XkaAFu6dCnLz8/Xj+gU6jFOrKOgSADS0tIYAKNp7969+jIA2Icffqj/t1arZfPmzWNhYWFMpVKxvn37Gt013b59m02fPp0FBQUxHx8flpKSwgoKCnir97Vr19j48eOZv78/8/f3Z+PHjzca4tu43owx9v777zMfHx+Tz8EpKChgffv2ZUFBQUypVLI2bdqwGTNmGD0PyNV1/+uvv1hycjJr3rw58/b2Zq1atWJpaWlG+1No+/zs2bMmj62Gx5cz9vny5ctZVFQUUyqVrHv37gYZp7S0NNavXz+D8vv27WPdunVjSqWSRUdHmwyYP/vsMxYbG8u8vb1Zu3btDC7gfLGn3v369TO5X9PS0vRlZs6cyVq1asWUSiVr3rw5S05OZrm5ubzX2966L1myhLVp04ap1WoWGBjI7r77brZjxw6jZQptnzNWn5H18fFhq1evNrk8V+xzXabK3Hcv5GOcWKZg7O/eXoQQQgghMkbPKSKEEEIIAQVFhBBCCCEAKCgihBBCCAFAQREhhBBCCAAKigghhBBCAFBQRAghhBACgIIiQgghhBAAFBQRQgghhACgoIgQQgghBAAFRYQQQgghACgoIoQQQggBAPw/SCMooMZ0RAQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Absolute difference between U3 and its first order taylor expansion\")\n",
    "plt.contourf(*X, np.abs(U3(x_ad).as_func(H) - U3(X)), levels=20);\n",
    "plt.axis('equal'); plt.scatter(*x,color='red'); plt.colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:36.996575Z",
     "iopub.status.busy": "2024-04-30T09:22:36.996462Z",
     "iopub.status.idle": "2024-04-30T09:22:37.121135Z",
     "shell.execute_reply": "2024-04-30T09:22:37.120836Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAGxCAYAAADF4QrHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB2AUlEQVR4nO3deXgTVdsG8Dvd0o0GSukClhaQfacItMgOBQR3BQQLKKIVUcvigqgsr6+IviKKLKJgRRGKsrhQkbKjFAUsKLKIChSkpRZKW9CuOd8f0HxNk7ZJZpJMMvfPa65LpmcmZyaTmWeec+aMRgghQEREREQO4+HsChARERGpDQMwIiIiIgdjAEZERETkYAzAiIiIiByMARgRERGRgzEAIyIiInIwBmBEREREDsYAjIiIiMjBGIAREREROZjNAdg777wDjUaDdu3amf37mTNnoNFo8L///c/mytli9uzZ0Gg0Ni177NgxzJ49G2fOnJG3UtXYtWsXNBoNdu3aZZg3fvx4REdHG5W7fPkyRo0ahdDQUGg0Gtx1110Aru/jYcOGITg4GBqNBklJSQ6ptxLZ43i7cOECZs+ejcOHD8u2TqXp27dvtb/h3NxcaDQazJ492zDv8OHDGDZsGBo3bgw/Pz8EBwcjNjYWn3zyiYNqbLno6GiMHz++1nJVt9HR5wFXV3X/OZqUc761UlNTHbatffv2Rd++fR3yWWqjlH3rZeuCK1euBAD8+uuv+OGHH9C9e3fZKuUsx44dw5w5c9C3b1+TIMhRXnrpJTz99NNG8/7zn/9g48aNWLlyJZo1a4bg4GAAwJQpU/DDDz9g5cqVCA8PR0REhDOq7LYuXLiAOXPmIDo6Gp06dXJ2dRThypUriIyMxAMPPIBGjRrh2rVrWL16NRISEnDmzBm8+OKLzq6i1dLT03HTTTcZ/q2E8wApU2pqKhYvXuzUgJOkW7JkibOrAMDGAOzgwYM4cuQIhg0bhs2bN2PFihVuEYApQbNmzUzmHT16FM2aNcOYMWNM5nfr1s2QEZNKCIGioiL4+fnJsj5yP+buHIcPH47Tp09j+fLlLhmA9ejRw9lVoBqUl5ejrKwMWq3WIZ/3zz//wN/f3yGf5Sg8txtr06aNs6sAwMYmyBUrVgAAXnvtNcTFxWHt2rX4559/zJbV6/X473//i8aNG8PX1xddu3bF9u3bjcr8/fffePTRRxEZGQmtVosGDRqgZ8+e2LZtm1G5lStXomPHjvD19UVwcDDuvvtuHD9+vNb6Vpcir9xEkZycjPvvvx8A0K9fP2g0Gmg0GiQnJxvKb9u2DQMGDEBQUBD8/f3Rs2dPk22pzokTJzBkyBD4+/sjJCQEiYmJKCwsNClXuQmyollt27ZtOH78uKFOFU2Xv//+O7755hvD/Iomk4KCAkyfPh1NmjSBj48PGjVqhKSkJFy7ds1kv0yePBnLli1D69atodVq8dFHHwEATp06hdGjRyM0NBRarRatW7fG4sWLjZavqMeaNWswc+ZMNGzYEEFBQRg4cCBOnjxpsm1btmzBgAEDoNPp4O/vj9atW2PevHlGZQ4ePIg77rgDwcHB8PX1RefOnbFu3TqL9jFg2fFmyfbt2rULt9xyCwDgoYceMuzj2bNnY/PmzdBoNDhw4ICh/Pr166HRaDBs2DCjz+nQoQPuvfdew7+FEFiyZAk6deoEPz8/1KtXD/fddx/+/PNPkzpacrxVNL/8+uuveOCBB6DT6RAWFoaHH34Y+fn5Fu83qUJCQuDlVfv93MGDBzFq1ChER0fDz88P0dHReOCBB3D27FmjcsnJydBoNNi5cycef/xxhISEoH79+rjnnntw4cIFo7KlpaV49tlnER4eDn9/f9x666348ccfLa575fNDbeeBjIwMDB8+3HDcNGzYEMOGDcP58+dr/AxLlrPm2LDkt/Tll18iNjYW/v7+qFOnDgYNGoT09HSjMtYcPwUFBZg4cSLq16+PwMBADBkyBL/99ptF+xgAMjMz8eCDDxr95t58803o9XpDmYpz3uuvv45XXnkFTZo0gVarxc6dOwEAmzdvRqdOnaDVatGkSZNquxxYui8rmuD37NmDuLg4+Pv74+GHHza7zvHjxxvOERXHReXz7uLFi9G7d2+EhoYiICAA7du3x+uvv47S0lLDOv7zn//Ay8sL586dM1n/ww8/jPr166OoqKjafXj58mVMmjQJjRo1go+PD5o2bYqZM2eiuLjYqFxN5/bqpKSkIDY2FgEBAQgMDMTgwYORkZFh+Pt3330Hb29vTJ8+3Wi5it9qRVxQ+fPfe+89tGjRAlqtFm3atMHatWuNlv37778xadIktGnTBoGBgQgNDUX//v2xd+9eo3KVu5gsWLAATZo0QWBgIGJjY7F//36jsn/++SdGjRqFhg0bQqvVIiwsDAMGDDDqSmLuRtLaffvxxx+jdevW8Pf3R8eOHfH111/XuH/NElb6559/hE6nE7fccosQQogPPvhAABDJyclG5U6fPi0AiMjISHHrrbeK9evXi88++0zccsstwtvbW+zbt89QdvDgwaJBgwZi+fLlYteuXWLTpk3i5ZdfFmvXrjWUefXVVwUA8cADD4jNmzeLVatWiaZNmwqdTid+++03Q7lZs2aJqpsFQMyaNctkW6KiosS4ceOEEELk5OQYPmPx4sUiPT1dpKeni5ycHCGEEB9//LHQaDTirrvuEhs2bBBfffWVGD58uPD09BTbtm2rcZ9lZ2eL0NBQ0ahRI/Hhhx+K1NRUMWbMGNG4cWMBQOzcudNQdty4cSIqKkoIIURRUZFIT08XnTt3Fk2bNjXUKT8/X6Snp4vw8HDRs2dPw/yioiJx7do10alTJxESEiIWLFggtm3bJt5++22h0+lE//79hV6vN9ovjRo1Eh06dBCffvqp2LFjhzh69Kj49ddfhU6nE+3btxerVq0SW7duFdOmTRMeHh5i9uzZhuV37twpAIjo6GgxZswYsXnzZrFmzRrRuHFj0bx5c1FWVmYo+8EHHwiNRiP69u0rPv30U7Ft2zaxZMkSMWnSJEOZHTt2CB8fH9GrVy+RkpIitmzZIsaPHy8AiA8//LDGfWzN8WbJ9uXn54sPP/xQABAvvviiYR+fO3dOFBYWCm9vb/Hqq68a1pmYmCj8/PxEQECAKCkpEUIIcfHiRaHRaMSSJUsM5SZOnCi8vb3FtGnTxJYtW8Snn34qWrVqJcLCwkR2drahnKXHW8Xx3rJlS/Hyyy+LtLQ0sWDBAqHVasVDDz1U4z4TQog+ffqItm3bmv3b33//Xe1vp7y8XJSWloqcnByxePFi4eXlJZYtW1br53322Wfi5ZdfFhs3bhS7d+8Wa9euFX369BENGjQQf//9t6Fcxb5v2rSpePLJJ8W3334rPvjgA1GvXj3Rr18/o3WOGzdOaDQa8cwzz4itW7eKBQsWiEaNGomgoCDD77smlbexpvPA1atXRf369UXXrl3FunXrxO7du0VKSopITEwUx44dq3b9li5n6bFhyW9p9erVAoCIj48XmzZtEikpKSImJkb4+PiIvXv3GspZevzo9XrRr18/odVqxX//+1+xdetWMWvWLNG0adNqj5HKcnJyRKNGjUSDBg3EsmXLxJYtW8TkyZMFAPH4448bylX8jhs1aiT69esnPv/8c7F161Zx+vRpsW3bNuHp6SluvfVWsWHDBsPvu+I8Wpml+7JPnz4iODhYREZGikWLFomdO3eK3bt3m92G33//Xdx3330CgOG4qDjvCiHElClTxNKlS8WWLVvEjh07xFtvvSVCQkKM9uPFixeFVqsVM2fONFr3pUuXhJ+fn3jmmWeM6tanTx/Dv//991/RoUMHERAQIP73v/+JrVu3ipdeekl4eXmJ2267zWh91Z3bq/Pf//5XaDQa8fDDD4uvv/5abNiwQcTGxoqAgADx66+/Gsq99tprAoD44osvhBBCHD16VPj7+4sHH3zQ5PMjIyNFmzZtxJo1a8SXX34phgwZIgCIzz77zFDuxIkT4vHHHxdr164Vu3btEl9//bWYMGGC8PDwMLouVhwX0dHRYsiQIWLTpk1i06ZNon379qJevXriypUrhrItW7YUN998s/j444/F7t27xfr168W0adOM1id130ZHR4tu3bqJdevWidTUVNG3b1/h5eUl/vjjj2r3sTlWB2CrVq0SAAwn28LCQhEYGCh69eplVK5ihzVs2FD8+++/hvkFBQUiODhYDBw40DAvMDBQJCUlVfuZeXl5ws/Pz2RHZGZmCq1WK0aPHm2YZ2sAJsT1i0PVgEgIIa5duyaCg4PF7bffbjS/vLxcdOzYUXTr1q3augshxHPPPSc0Go04fPiw0fxBgwbVGIBVqO4iGRUVJYYNG2Y0b968ecLDw0McOHDAaP7nn38uAIjU1FTDPABCp9OJy5cvG5UdPHiwuOmmm0R+fr7R/MmTJwtfX19D+YoArOr3sm7dOsNJSojrx0hQUJC49dZbjQLAqlq1aiU6d+4sSktLjeYPHz5cREREiPLy8mqXteZ4s3T7Dhw4UG3wd+utt4r+/fsb/n3zzTeLZ555Rnh4eBhO4BUXwYobhPT0dAFAvPnmm0brOnfunPDz8xPPPvusEMK6463ieH/99deNyk6aNEn4+vrWuL+FsD0Ae+yxxwQAAUD4+PgYBZnWKCsrE1evXhUBAQHi7bffNsyvCMAqBxVCCPH6668LACIrK0sIIcTx48cFADFlyhSjchX73toATIjqzwMHDx4UAMSmTZus2kZLlrP02LDkt1ReXi4aNmwo2rdvb/SbKSwsFKGhoSIuLs4wz9Lj55tvvhEAjL4jIa5fuC0JwJ5//nkBQPzwww9G8x9//HGh0WjEyZMnhRD//ztu1qyZ4UamQvfu3av9fVc+51u6L4W4fvwDENu3b6+x/hWeeOIJk+uLORU3KKtWrRKenp5G59hx48aJ0NBQUVxcbJg3f/584eHhIU6fPm1Ut8pBwrJlywQAsW7dOqPPmj9/vgAgtm7daphX3bndnMzMTOHl5SWefPJJo/mFhYUiPDxcjBgxwjBPr9eL2267TdStW1ccPXpUtGnTRrRq1UpcvXrVaFkAws/PzyjYLSsrE61atRI333xztXUpKysTpaWlYsCAAeLuu+82zK84Ltq3b290Y//jjz8KAGLNmjVCCCFyc3MFALFw4cIat1nqvg0LCxMFBQWGednZ2cLDw0PMmzevxs+tyuomyBUrVsDPzw+jRo0CAAQGBuL+++/H3r17cerUKZPy99xzD3x9fQ3/rlOnDm6//Xbs2bMH5eXlAIBu3bohOTkZr7zyCvbv32+UsgWud5L9999/TZ5oioyMRP/+/S1uBrTVvn37cPnyZYwbNw5lZWWGSa/XY8iQIThw4IBJ815lO3fuRNu2bdGxY0ej+aNHj5a9rl9//TXatWuHTp06GdV18ODBJk9cAkD//v1Rr149w7+Lioqwfft23H333fD39zdax2233YaioiKTlO8dd9xh9O8OHToAgKFZad++fSgoKMCkSZOqfVrp999/x4kTJwz93Kp+blZWltlmzapqO95s2T5zBgwYgO+//x7//vsvzp49i99//x2jRo1Cp06dkJaWBuB6E2Ljxo3RvHlzANe/G41GgwcffNDoc8PDw9GxY0fDd2PL8WbuOygqKkJOTk6t22KLF154AQcOHMDmzZvx8MMPY/LkyRY9gXr16lU899xzuPnmm+Hl5QUvLy8EBgbi2rVrZrsT1HZsVTRNVe0fOWLECIuaRK1x8803o169enjuueewbNkyHDt2TLblrDk2avstnTx5EhcuXEBCQgI8PP7/FB8YGIh7770X+/fvN+kyUtvxU91+tvQctmPHDrRp0wbdunUzmj9+/HgIIbBjxw6T+nh7exv+fe3aNRw4cKDa33dllu7LCvXq1UP//v0t2o6aZGRk4I477kD9+vXh6ekJb29vjB07FuXl5UZNtU8//TRycnLw2WefAbjebWLp0qUYNmxYjQ997NixAwEBAbjvvvuM5ldcF6teB6ue26vz7bffoqysDGPHjjXaX76+vujTp4/R/tJoNFi1ahXq1KmDrl274vTp01i3bh0CAgJM1jtgwACEhYUZ/u3p6YmRI0fi999/N2p6X7ZsGbp06QJfX194eXnB29sb27dvN3s+GDZsGDw9PQ3/rno+CA4ORrNmzfDGG29gwYIFyMjIMGriro61+7Zfv36oU6eO4d9hYWEIDQ016UpRG6sCsN9//x179uzBsGHDIITAlStXcOXKFUOlK56MrCw8PNzsvJKSEly9ehXA9bbncePG4YMPPkBsbCyCg4MxduxYZGdnAwAuXboEAGaf8mvYsKHh7/Zy8eJFAMB9990Hb29vo2n+/PkQQuDy5cvVLn/p0qVq94M96vrzzz+b1LNOnToQQiA3N9eofNV9eunSJZSVlWHRokUm67jtttsAwGQd9evXN/p3RWfZf//9F8D1dn4ARk+amas3AEyfPt3kcydNmmT2c82p7XizZfvMGThwIIqLi/Hdd98hLS0NISEh6Ny5MwYOHGjou7h9+3YMHDjQaBuFEAgLCzP57P379xs+15bjrbbvoDpeXl6GG6GqysrKAMDoQlihcePG6Nq1K2677TYsXboUjz76KGbMmGH4rqszevRovPvuu3jkkUfw7bff4scff8SBAwfQoEEDs3WtbbsqfvtVv3cvLy+TZaXS6XTYvXs3OnXqhBdeeAFt27ZFw4YNMWvWLJObRmuXs/TYsOS3VNv5Uq/XIy8vz2i+JfvZ3D619Bx26dKlautTuc4VqpbNy8uDXq+36Dxq6b6s7rNskZmZiV69euGvv/7C22+/jb179+LAgQOGPmOVj+3OnTujV69ehr99/fXXOHPmDCZPnlzjZ1RcR6oG3qGhofDy8qp1H1an4nxzyy23mOyvlJQUs+f7O+64A0VFRRgyZAjat29vdr01fVcVdV2wYAEef/xxdO/eHevXr8f+/ftx4MABDBkyxKbzgUajwfbt2zF48GC8/vrr6NKlCxo0aICnnnrKbJ/rCtbuW3PnFq1WW+v5tiqrbhFXrlwJIQQ+//xzfP755yZ//+ijj/DKK68YRagVQVRl2dnZ8PHxQWBgIIDrHXgXLlyIhQsXIjMzE19++SWef/555OTkYMuWLYaNzcrKMlnXhQsXEBISUmO9tVqtSUc6wPRHX52K9S9atKjaJ6YqR/pV1a9fv9r9ILeQkBD4+fmZDYYr/l5Z1QOuXr168PT0REJCAp544gmz62jSpIlVdWrQoAEA1NhRuaJeM2bMwD333GO2TMuWLWv9rNqON29vb1m2r3v37ggMDMS2bdtw5swZDBgwABqNBgMGDMCbb76JAwcOIDMz0ygACwkJgUajwd69e80+0VUxT+rxZo2wsDAcOHAAQgiTY+Gvv/6y+LO6deuGZcuW4c8//zR831Xl5+fj66+/xqxZs/D8888b5hcXF9d4A1OTinNDdnY2GjVqZJhfVlZmlxuz9u3bY+3atRBC4Oeff0ZycjLmzp0LPz8/o22ydjlLjw1Lfku1nS89PDwsyoxUXWfFPq188bH0HFa/fv1q6wNYdl7SaDQWnUct3ZfVfZYtNm3ahGvXrmHDhg2IiooyzK9uDMGnnnoK999/P3766Se8++67aNGiBQYNGlTjZ9SvXx8//PCDyW81JycHZWVlte7D6lQs9/nnnxvVvTppaWlYunQpunXrho0bN2L9+vVGDxpVqOm7qjiGPvnkE/Tt2xdLly41KldTsFSbqKgowwMBv/32G9atW4fZs2ejpKQEy5YtM7uMtftWLhZnwMrLy/HRRx+hWbNm2Llzp8k0bdo0ZGVl4ZtvvjFabsOGDUZPdRQWFuKrr75Cr169jAK1Co0bN8bkyZMxaNAg/PTTTwCA2NhY+Pn5mQz2eP78eezYsQMDBgyose7R0dH4+eefjebt2LHDkIGrUF3WoGfPnqhbty6OHTuGrl27mp18fHyq/fx+/frh119/xZEjR4zmf/rppzXW2xbDhw/HH3/8gfr165utZ23jGvn7+6Nfv37IyMhAhw4dzK7D2sxCXFwcdDodli1bhutN6KZatmyJ5s2b48iRI9Xu48op3+rUdrxZs301ZZG8vb3Ru3dvpKWlYceOHYaTZ69eveDl5YUXX3zREJBVGD58OIQQ+Ouvv8x+bsWdpNTjzRoDBw5EQUEBtmzZYvK3devWwcPDw6LmmZ07d8LDwwNNmzattoxGo4EQwuQC+MEHH1SbhatNxZNMq1evNpq/bt06QwbPWpZkDzUaDTp27Ii33noLdevWNZyralPdcpYeG5b+lho1aoRPP/3UqMy1a9ewfv16w5OR1ujXrx8A0/1s6TlswIABOHbsmMl+WrVqFTQajWH91QkICEC3bt2q/X1XZum+tEV1x0bFRbvysS2EwPvvv292PXfffTcaN26MadOmYdu2bTU2KVcYMGAArl69ik2bNhnNX7VqleHvthg8eDC8vLzwxx9/VHu+qZCVlYUHH3wQffr0wb59+3DHHXdgwoQJOH36tMl6t2/fbsiuAddjiJSUFDRr1syQwdVoNCbng59//tnkaV1btWjRAi+++CLat29f42/UXvu2NhZnwL755htcuHAB8+fPNzuCbLt27fDuu+9ixYoVGD58uGG+p6cnBg0ahKlTp0Kv12P+/PkoKCjAnDlzAFy/K+7Xrx9Gjx6NVq1aoU6dOjhw4AC2bNliyITUrVsXL730El544QWMHTsWDzzwAC5duoQ5c+bA19cXs2bNqrHuCQkJeOmll/Dyyy+jT58+OHbsGN59913odDqTbQCA5cuXo06dOvD19UWTJk1Qv359LFq0COPGjcPly5dx3333ITQ0FH///TeOHDmCv//+2ySCrywpKQkrV67EsGHD8MorryAsLAyrV6/GiRMnLNr31khKSsL69evRu3dvTJkyBR06dIBer0dmZia2bt2KadOm1Tpm29tvv41bb70VvXr1wuOPP47o6GgUFhbi999/x1dffWXSX6M2gYGBePPNN/HII49g4MCBmDhxIsLCwvD777/jyJEjePfddwEA7733HoYOHYrBgwdj/PjxaNSoES5fvozjx4/jp59+MvSZqEltx5s129esWTP4+flh9erVaN26NQIDA9GwYUNDs8mAAQMwbdo0ADBkuvz8/BAXF4etW7eiQ4cOCA0NNXxuz5498eijj+Khhx7CwYMH0bt3bwQEBCArKwvfffcd2rdvj8cffxyBgYGSjjdrjBkzBkuWLMGIESPw/PPP45ZbbsG///6L1NRUvP/++3jyySeNgqpHH30UQUFB6NatG8LCwpCbm4vPPvsMKSkpeOaZZ6rNfgFAUFAQevfujTfeeAMhISGIjo7G7t27sWLFCtStW9em+rdu3RoPPvggFi5cCG9vbwwcOBBHjx7F//73PwQFBdm0zurOA+np6ViyZAnuuusuNG3aFEIIbNiwAVeuXKkxe/H111/Xupw1x0ZtvyUPDw+8/vrrGDNmDIYPH47HHnsMxcXFeOONN3DlyhW89tprVu+T+Ph49O7dG88++yyuXbuGrl274vvvv8fHH39s0fJTpkzBqlWrMGzYMMydOxdRUVHYvHkzlixZgscffxwtWrSodR3/+c9/MGTIEAwaNAjTpk1DeXk55s+fj4CAAKMMqqX70hYVwdv8+fMxdOhQeHp6okOHDhg0aBB8fHzwwAMP4Nlnn0VRURGWLl1q0tRbwdPTE0888QSee+45BAQEWPTGhrFjx2Lx4sUYN24czpw5g/bt2+O7777Dq6++ittuu80o226N6OhozJ07FzNnzsSff/6JIUOGoF69erh48SJ+/PFHBAQEYM6cOSgvL8cDDzwAjUaDTz/9FJ6enkhOTkanTp0wcuRIfPfdd0Y3hiEhIejfvz9eeuklBAQEYMmSJThx4oTRUBTDhw/Hf/7zH8yaNQt9+vTByZMnMXfuXDRp0sSmG6iff/4ZkydPxv3334/mzZvDx8cHO3bswM8//1xjhtpe+7ZWlvbWv+uuu4SPj49hWAZzRo0aJby8vER2drbhqYX58+eLOXPmiJtuukn4+PiIzp07i2+//dawTFFRkUhMTBQdOnQQQUFBws/PT7Rs2VLMmjVLXLt2zWj9H3zwgejQoYPw8fEROp1O3HnnnUaPyAph/inI4uJi8eyzz4rIyEjh5+cn+vTpIw4fPmzyFKQQQixcuFA0adJEeHp6mjwBt3v3bjFs2DARHBwsvL29RaNGjcSwYcOMHqutzrFjx8SgQYOEr6+vCA4OFhMmTBBffPGF7E9BCnH9sfcXX3xRtGzZ0rCv2rdvL6ZMmWL0VAoA8cQTT5it7+nTp8XDDz8sGjVqJLy9vUWDBg1EXFyceOWVVwxlKp6CrLr9Fd991acHU1NTRZ8+fURAQIDw9/cXbdq0EfPnzzcqc+TIETFixAgRGhoqvL29RXh4uOjfv3+tQxxYerxZs31CCLFmzRrRqlUr4e3tbfK015EjRwQA0bx5c6NlKp4Mmzp1qtm6rly5UnTv3l0EBAQIPz8/0axZMzF27Fhx8OBBo3KWHG8Vx3vlIRyE+P+nCCs/VVWdgoIC8eyzz4rmzZsLHx8f4e/vL7p27SqWLVtm8qTdypUrRa9evURISIjw8vISdevWFX369BEff/xxrZ8jhBDnz58X9957r6hXr56oU6eOGDJkiDh69KjJb7Gi/lWf5q045ir/ZoqLi8W0adNEaGio8PX1FT169BDp6elmf9/mVP1ehTB/Hjhx4oR44IEHRLNmzYSfn5/Q6XSiW7duJkPwVGXNcpYeG5b8ljZt2iS6d+8ufH19RUBAgBgwYID4/vvvjcpYc/xcuXJFPPzww6Ju3brC399fDBo0SJw4ccKipyCFEOLs2bNi9OjRon79+sLb21u0bNlSvPHGG0ZPalb8jt944w2z6/jyyy8N14DGjRuL1157zew539J9WdNTwOYUFxeLRx55RDRo0EBoNBqjffTVV1+Jjh07Cl9fX9GoUSPxzDPPGJ4erfpErRBCnDlzRgAQiYmJZj+r6pN6QlwfriIxMVFEREQILy8vERUVJWbMmGEYCqNCTef26mzatEn069dPBAUFCa1WK6KiosR9991nGPZm5syZwsPDw+SJ0X379gkvLy/x9NNPm3z+kiVLRLNmzYS3t7do1aqVWL16tdGyxcXFYvr06aJRo0bC19dXdOnSRWzatMnkOljTcVH5+Lt48aIYP368aNWqlQgICBCBgYGiQ4cO4q233jJ6etIe+9bS801lmhsrJCIiIgdZtGgRnnrqKRw9ehRt27Z1dnVkpdFo8MQTTxhaN8g8eZ/TJiIiomplZGTg9OnTmDt3Lu688063C77IcgzAiIiIHOTuu+9GdnY2evXqVe1TeaQObIIkIiIicjCbXsbtjvbs2YPbb78dDRs2hEajMXkc1Zzdu3cjJiYGvr6+aNq0qdm7mfXr16NNmzaGl5Fu3LjRDrUnIiJyf/a6VjsDA7Abrl27ho4dO1rcafD06dO47bbb0KtXL2RkZOCFF17AU089hfXr1xvKpKenY+TIkUhISMCRI0eQkJCAESNG4IcffrDXZhAREbkte1yrnYVNkGZoNBps3LgRd911V7VlnnvuOXz55ZdG76tKTEzEkSNHDIPIjRw5EgUFBUaD01aMsbJmzRq71Z+IiMjdyXWtdhZ2wrdReno64uPjjeYNHjwYK1asQGlpKby9vZGeno4pU6aYlFm4cGGN6y4uLjZ6dZJer8fly5dRv359WV6bQURE7kkIgcLCQjRs2NDoZexyKyoqQklJiSzrEmZehabVas2+SspallyrnYUBmI2ys7NN3pEXFhaGsrIy5ObmIiIiotoytb0/bd68eUYjtxMREVnj3LlzNb60XYqioiI0iQpEdo5trxCrKjAw0OTVgLNmzcLs2bMlr9uSa7WzMACToGrEXtGaW3m+uTK1ZbFmzJiBqVOnGv6dn5+Pxo0b49y5c1a9XqXTwsUAAK+rtRS8wcfK95/6FFrfeu1ToLd6mcq0V0olLV/B64rpy9nl5HGlwK7rt5T+kvlXoRC5C4/61r1Y3J70dW17/VVNNn73glXlCwoKEBkZadG7c21VUlKC7JxynD4UhaA60rJsBYV6NIk5a3J9kyP7VcGSa7UzMACzUXh4uEkmKycnB15eXoaXOVdXpmo0XlV1qdegoCCrAjAPX194FwKw8Dj2tDKb7FlsfQDm5W17AKbNKwW8TF/gbguZVmOWR14B4CHfyUMKvUael3YTKdbla/AICXZ2La4rKIa+nrxBmK3vNHVEcBFUx0NyAGZYl5XXN0tZcq12Fj4FaaPY2FikpaUZzdu6dSu6du1qaFOurkxcXJzD6mkvPgWu++yGV16Rs6vgEPrcy7UXIiKyUbnQyzLZkyXXamdhAHbD1atXcfjwYRw+fBjA9UdXDx8+jMzMTADXmwXHjh1rKJ+YmIizZ89i6tSpOH78OFauXIkVK1Zg+vTphjJPP/00tm7divnz5+PEiROYP38+tm3bhqSkJEdumkWsbX60hTbfvj80JfDIU0bTI5Ga8GbDOfQQskzWsMe12lkYgN1w8OBBdO7cGZ07dwYATJ06FZ07d8bLL78MAMjKyjJ8wQDQpEkTpKamYteuXejUqRP+85//4J133sG9995rKBMXF4e1a9fiww8/RIcOHZCcnIyUlBR0797dIdvk7YCgylG0efL0/SIisic13YTpZfrPGva4VjsLxwFzAQUFBdDpdMjPz7eqjbz1S29ZXNbqDvg2NEFKyYDJFYDZu/lRKSdfZgRIjZTSF0zOfmDfZlj3RLyt1wtbPuPCyZtk6YTfsOV5u9ZXqdgJn6zmyv2/7EkpwRcROZdHXoHsnfGVqFwIlEvM4Uhd3pUxACPF9/9i86N1mP0itdLnXlZMFkwNbOnDZW4dasU+YKQa9mx+ZPaLiIiswQCMyI0w+0Vqp5TfgBpuyvQQKJc4qTkDxiZIsoqjO98TEZEysQlSGmbAVM4R/b+kcIWnH5Vyp6uUO38iZ1PKb0Ep5wZSJmbAiIiIyGp8ClIaBmBEbkApd/xEpB76G5PUdagVmyDJYhx81RSbGIiUSSk3JTxHUHWYAVMxpff/IiIi5ap4klHqOtSKARiRi1PKnT6R0nBgVvsqF9cnqetQKzZBEtmITQtEZAl3PVfoZZrUigEYWYT9v5SJ2S+imvE3QkrFJkgiIiKymh4alEMjeR1qxQBMpdgBXxolNCnwzp7IdXjkFUBfL8jZ1ZCVXlyfpK5DrdgESW7L3ZsficgyvFkhJWIGjGrlqv2/iIjIfsplaIKUurwrYwaMyEpsfiRyPUr4zSjh3CGnigBM6qRWDMCIiIiIHIxNkCqkhg747P9FRGRfeqGBXkh8ClLi8q6MARgpitL7fymhCUEJTSlErogj48uLfcCkYRMk1cjRHfCJiNyZEm7iSBmYASMiIiKrlcMD5RLzOOUy1cUVMQAjt+PO/b/Y/EgkDZsh5SNk6AMm2AeM1ELJHfDZ/4uI1MBdRsVnHzBp2AeMZMX+X/bD7BcRkftgBoyqZUsHfCIipWMzpDzKhQfKhcQ+YCq+zDAAI7dir/5fbH4kIjKmhwZ6iQ1peqg3AmMTJJELYPMjkXvhTR0xA0aKoPQO+ETkXtgMKR074UvDAExFlPwEJBERuRZ5+oCxCZJIMmc/Aemu/b/Y/EhE5H4YgJFZfAKSiNyds29unH1zJ9X1TvjSJ7ViEyQRERFZTS/Dq4j4FCSRE7EDPhERqQ0DsCqWLFmCJk2awNfXFzExMdi7d2+1ZcePHw+NRmMytW3b1lAmOTnZbJmiIvd9XyHJx9lNJETujr8x21V0wpc6qZV6t9yMlJQUJCUlYebMmcjIyECvXr0wdOhQZGZmmi3/9ttvIysryzCdO3cOwcHBuP/++43KBQUFGZXLysqCr6+vIzbJwN2fgHTXDvhE5N5c+Ryjh4csk1qxD1glCxYswIQJE/DII48AABYuXIhvv/0WS5cuxbx580zK63Q66HQ6w783bdqEvLw8PPTQQ0blNBoNwsPD7Vt5J3P2E5BERORY5UKDciFxHDCJy7sy9YaeVZSUlODQoUOIj483mh8fH499+/ZZtI4VK1Zg4MCBiIqKMpp/9epVREVF4aabbsLw4cORkZFR43qKi4tRUFBgNDkSn4BUBjaNEDkGf2vkDAzAbsjNzUV5eTnCwsKM5oeFhSE7O7vW5bOysvDNN98YsmcVWrVqheTkZHz55ZdYs2YNfH190bNnT5w6daradc2bN8+QXdPpdIiMjLRto4iIiOyk/MZTkFIntVLvlldDozFOhwohTOaZk5ycjLp16+Kuu+4ymt+jRw88+OCD6NixI3r16oV169ahRYsWWLRoUbXrmjFjBvLz8w3TuXPnbNoWV8AnIImIXJNeeMgyqRX7gN0QEhICT09Pk2xXTk6OSVasKiEEVq5ciYSEBPj4+NRY1sPDA7fcckuNGTCtVgutVmt55ckuXLlzLBG5Do+8AujrBTm7GuRg6g09q/Dx8UFMTAzS0tKM5qelpSEuLq7GZXfv3o3ff/8dEyZMqPVzhBA4fPgwIiIiJNWX/p+9noB0JvZJISKlYxOkNMyAVTJ16lQkJCSga9euiI2NxfLly5GZmYnExEQA15sG//rrL6xatcpouRUrVqB79+5o166dyTrnzJmDHj16oHnz5igoKMA777yDw4cPY/HixQ7ZJsD+Q1DwCUgicnX63MvwCAl2djVcih7Sn2JU89WDAVglI0eOxKVLlzB37lxkZWWhXbt2SE1NNTzVmJWVZTImWH5+PtavX4+3337b7DqvXLmCRx99FNnZ2dDpdOjcuTP27NmDbt262X17iIiISJkYgFUxadIkTJo0yezfkpOTTebpdDr8888/1a7vrbfewltvvSVX9YiIiBRBjoFUORAr0Q2OHANMyU9AsgM+ETmSK3bEl+NVQnwVEREpBjvgEzkHf3vkSMyAkUtzxycgiYhcgR4a6CG1E756X0XEAIyIiIisxiZIaRiAERERkdXkGMdLzeOAqXfLSRYcA4yIiMh6DMDcnL0HYXVHznwCkp2AiZzLmb9BV3v6Wi80skxqxSZIIiIisppehiZINY8Dpt4tJ6dS8hhgRERE9sYMGBk4chBWOXAICiIi59ELD+glPsUodXlXxgCMiIiIrFYODcoljuMldXlXpt7Qk0hh2AGfiEg9mAEjIiKqRJ97GR4hwc6uhuKxCVIaBmBElbjaY+BERM5SDulNiOXyVMUlqTf0JMk4CCsRkbx4E6gezIC5MQ7CSkRE9sImSGkYgBEREZHV+DJuadS75URERGQzAQ30EidhQx+yJUuWoEmTJvD19UVMTAz27t1bY/nVq1ejY8eO8Pf3R0REBB566CFcunTJ1s2WDQMwcjiOgm+KQ1AQEdUuJSUFSUlJmDlzJjIyMtCrVy8MHToUmZmZZst/9913GDt2LCZMmIBff/0Vn332GQ4cOIBHHnnEwTU3xQCMXBJHwScie+JNUe0qmiClTtZYsGABJkyYgEceeQStW7fGwoULERkZiaVLl5otv3//fkRHR+Opp55CkyZNcOutt+Kxxx7DwYMH5dgFkjAAIyIiIqvphUaWCQAKCgqMpuLiYpPPKykpwaFDhxAfH280Pz4+Hvv27TNbx7i4OJw/fx6pqakQQuDixYv4/PPPMWzYMPl3iJUYgBEREZFTRUZGQqfTGaZ58+aZlMnNzUV5eTnCwsKM5oeFhSE7O9vseuPi4rB69WqMHDkSPj4+CA8PR926dbFo0SK7bIc1+BQkAXC9F3HbA8ffISKyXDk8UC4xj1Ox/Llz5xAUFGSYr9Vqq11GozHuuC+EMJlX4dixY3jqqafw8ssvY/DgwcjKysIzzzyDxMRErFixQlLdpWIARkRERFar3IQoZR0AEBQUZBSAmRMSEgJPT0+TbFdOTo5JVqzCvHnz0LNnTzzzzDMAgA4dOiAgIAC9evXCK6+8goiICEn1l4JNkERERArCbLx5Pj4+iImJQVpamtH8tLQ0xMXFmV3mn3/+gYeHcajj6ekJ4HrmzJmYASMiIiKr6eEBvcQ8jrXLT506FQkJCejatStiY2OxfPlyZGZmIjExEQAwY8YM/PXXX1i1ahUA4Pbbb8fEiROxdOlSQxNkUlISunXrhoYNG0qqu1QMwIiIiMhq5UKDcolNkNYuP3LkSFy6dAlz585FVlYW2rVrh9TUVERFRQEAsrKyjMYEGz9+PAoLC/Huu+9i2rRpqFu3Lvr374/58+dLqrccGIARERGRy5g0aRImTZpk9m/Jyckm85588kk8+eSTdq6V9RiAETkZB3wkIlckZyd8NWIARkRERFYTwgN6iS/TFip+GTcDMCIiIjP0uZfhERLs7GooVjk0KLfhZdpV16FW6g09iYiIiJyEGTCyiTZf7+wqEBGRE+mF9D5cehW/hIUBGBEREVlNL0MfMKnLuzIGYEREMvAQerTT5yBY/IvLGj8c9QiFXqPeiwsR1YxnhyqWLFmCJk2awNfXFzExMdi7d2+1ZXft2gWNRmMynThxwqjc+vXr0aZNG2i1WrRp0wYbN26092YQkQP1LDuLVf9uwBtFWzGjeC/eKNqKVf9uQM+ys86uGpHd6KGRZVIrBmCVpKSkICkpCTNnzkRGRgZ69eqFoUOHGo2qa87JkyeRlZVlmJo3b274W3p6OkaOHImEhAQcOXIECQkJGDFiBH744Qd7bw4ROUDPsrN4qXg3QsQ/RvPri3/wUvFuBmHktipGwpc6qRUDsEoWLFiACRMm4JFHHkHr1q2xcOFCREZGYunSpTUuFxoaivDwcMNU8aJPAFi4cCEGDRqEGTNmoFWrVpgxYwYGDBiAhQsX2nlriMjePIQej5ccAACT+3gPAAJAYskBeAg+tEJExhiA3VBSUoJDhw4hPj7eaH58fDz27dtX47KdO3dGREQEBgwYgJ07dxr9LT093WSdgwcPrnGdxcXFKCgoMJqISHna6XPQQPxTbSOKB4BQ8Q/a6XMcWS0ih6johC91Uiv1bnkVubm5KC8vR1hYmNH8sLAwZGdnm10mIiICy5cvx/r167Fhwwa0bNkSAwYMwJ49ewxlsrOzrVonAMybNw86nc4wRUZGStgyIrKXYPGvrOWIXIkeGsPriGyeVNwHjE9BVqHRGB8MQgiTeRVatmyJli1bGv4dGxuLc+fO4X//+x969+5t0zoBYMaMGZg6darh3wUFBQzCiBTossZP1nJEpB7MgN0QEhICT09Pk8xUTk6OSQarJj169MCpU6cM/w4PD7d6nVqtFkFBQUYTESnPUY9Q/K3xR3U9vPQAcjT+OOoR6shqETmEkOEJSKHiDBgDsBt8fHwQExODtLQ0o/lpaWmIi4uzeD0ZGRmIiIgw/Ds2NtZknVu3brVqnUSkTHqNB5b63AINYBKE6XG9Y/4yn1s4Hhi5JcnNjzcmtWITZCVTp05FQkICunbtitjYWCxfvhyZmZlITEwEcL1p8K+//sKqVasAXH/CMTo6Gm3btkVJSQk++eQTrF+/HuvXrzes8+mnn0bv3r0xf/583Hnnnfjiiy+wbds2fPfdd07ZRiKS1/deUfgP+uDxkgNoUGkoilyNP5b53ILvvaKcWDsi++FI+NIwAKtk5MiRuHTpEubOnYusrCy0a9cOqampiIq6fgLNysoyGhOspKQE06dPx19//QU/Pz+0bdsWmzdvxm233WYoExcXh7Vr1+LFF1/ESy+9hGbNmiElJQXdu3d3+PYRkX187xWFdM9IjoRPRBbTCCFU/CpM11BQUACdTof8/Hyr+oN1THrL4rI+BdYdBlJexq3NK7V52QpeeUWS11GVR55zhvvQ5152yucSUe08QoKd8rn6etfP9d9mzLFqOVuvF7Z8xp1bH4Z3gI+kdZVeK8EX8SvtWl+lYgaMbFKs85AUhBERkWuT41VCah6GgvlxIiIiIgdjBoyIiMgMZzU/ugo5nmLkU5BEREREVmAAJg2bIImIiIgcjBkwIifzCAnmk5BE5HKYAZOGARgRERFZjQGYNGyCJCIiInIwZsCIiIgUpGIQVqUTkD6Ol5pHgmcARkRERFZjE6Q0DMAIAFASpLH6dUTuRl8vyGmvIyIicjUMwKRhHzAiIiIiB2MGjIiIiKzGDJg0DMCIiIjIagzApGETJDlccT1vyesoq+crQ02Ug++cI1IW/ibJ3pgBIyIiIqsJoYGQmMGSurwrYwBGREREVtNDI3kcMKnLuzI2QbqxkjrOrgERERGZwwwY2axY5wFtvt7Z1SAichuuMgo+wE74UjEAI6qEg7ESEVmGfcCkYRMkERERkYMxA0akEB4hwdDnXnZ2NYhUj0NQWIZNkNIwACMiIiKrsQlSGgZgZOBqL+Quq+cLr7wiZ1eDiEiVhAwZMDUHYOwDRk4hx2j4RERErooZMCIiIrKaACAkNpq4TpuL/JgBc3McjNV6rjQODxHJy5kd8F3t3FMxEr7USa0YgBEpCJ++IiJSBzZBkiQcDZ+ISJ34FKQ0DMCIiIjIanqhgYbjgNmMTZDk0srq+Tq7CkRERFZjAEakMOwHRuQc/O1ZRwh5JrViAEZGSoIclw5W8lhgrvY0EhG5Nlc851T0AZM6qRUDMCIiIiIHYyd8IiIishqfgpSGGbAqlixZgiZNmsDX1xcxMTHYu3dvtWU3bNiAQYMGoUGDBggKCkJsbCy+/fZbozLJycnQaDQmU1GR495haO/BWIt1PIyIiNRGf+NdkFInteKVs5KUlBQkJSVh5syZyMjIQK9evTB06FBkZmaaLb9nzx4MGjQIqampOHToEPr164fbb78dGRkZRuWCgoKQlZVlNPn68uk9qh47AxM5Fn9z1mMnfGnYBFnJggULMGHCBDzyyCMAgIULF+Lbb7/F0qVLMW/ePJPyCxcuNPr3q6++ii+++AJfffUVOnfubJiv0WgQHh5u17qrWVk9X3jlOS6jSEQkJ1fsgE/SMQN2Q0lJCQ4dOoT4+Hij+fHx8di3b59F69Dr9SgsLERwsPGd1NWrVxEVFYWbbroJw4cPN8mQVVVcXIyCggKjyV3xSUgiItd0PYMl9SlIZ2+F8zAAuyE3Nxfl5eUICwszmh8WFobs7GyL1vHmm2/i2rVrGDFihGFeq1atkJycjC+//BJr1qyBr68vevbsiVOnTlW7nnnz5kGn0xmmyMhI2zaKiIjITjgMhTQMwKrQaIwPBiGEyTxz1qxZg9mzZyMlJQWhoaGG+T169MCDDz6Ijh07olevXli3bh1atGiBRYsWVbuuGTNmID8/3zCdO3fO9g2ygSPHAqPqsU8KkWPwt0bOwD5gN4SEhMDT09Mk25WTk2OSFasqJSUFEyZMwGeffYaBAwfWWNbDwwO33HJLjRkwrVYLrVZreeWJiIgcTNyYpK5DrZgBu8HHxwcxMTFIS0szmp+Wloa4uLhql1uzZg3Gjx+PTz/9FMOGDav1c4QQOHz4MCIiIiTX2RruPhSFvd4JyX5gRGRPrnyOYROkNMyAVTJ16lQkJCSga9euiI2NxfLly5GZmYnExEQA15sG//rrL6xatQrA9eBr7NixePvtt9GjRw9D9szPzw86nQ4AMGfOHPTo0QPNmzdHQUEB3nnnHRw+fBiLFy92zkYSEZEBmx/JWRiAVTJy5EhcunQJc+fORVZWFtq1a4fU1FRERUUBALKysozGBHvvvfdQVlaGJ554Ak888YRh/rhx45CcnAwAuHLlCh599FFkZ2dDp9Ohc+fO2LNnD7p16+bQbVOy4nre0OaVOrsaiuQREgx97mVnV4OIyBTbICVhE2QVkyZNwpkzZ1BcXIxDhw6hd+/ehr8lJydj165dhn/v2rULQgiTqSL4AoC33noLZ8+eRXFxMXJycvDtt98iNjbWgVtERERkB3I0P9rQBGnNG2uA60M7zZw5E1FRUdBqtWjWrBlWrlxp61bLhhkwMqskSAOfAhXfmhCR23N286Mr9/8C5BnJ3trlK95Ys2TJEvTs2RPvvfcehg4dimPHjqFx48ZmlxkxYgQuXryIFStW4Oabb0ZOTg7KysqkVVwGDMBINsU6D2jz9U77fHuNiK+vFwSPPOcNhstmSCKi66x9Y82WLVuwe/du/Pnnn4ZB0qOjox1Z5WqxCVJF7P0kJBERqYecT0FWfftLcXGxyefZ8saaL7/8El27dsXrr7+ORo0aoUWLFpg+fTr+/fdf+XeIlZgBI0VgR3wiciRnNz+6BRv7cJmsAzB548usWbMwe/Zso3m2vLHmzz//xHfffQdfX19s3LgRubm5mDRpEi5fvuz0fmAMwIhcAJshidyLq/f/ktu5c+cQFPT/+6SmwciteWONXq+HRqPB6tWrDcNDLViwAPfddx8WL14MPz8/GWpvGzZBElmAJ0siImMVnfClTgAQFBRkNJkLwGx5Y01ERAQaNWpkCL4AoHXr1hBC4Pz58/LtDBswAKNqueI7Ie01Ij4RuQ82P8pEyDRZyJY31vTs2RMXLlzA1atXDfN+++03eHh44KabbrL8w+2AARjJytmvJHJnvGgQkdpNnToVH3zwAVauXInjx49jypQpJm+sGTt2rKH86NGjUb9+fTz00EM4duwY9uzZg2eeeQYPP/ywU5sfAfYBU52SOoBPobNrYZ7SO+I7ezgKInIP7tKlQY53OVq7vLVvrAkMDERaWhqefPJJdO3aFfXr18eIESPwyiuvSKq3HBiAERGRajCTLDMnjNc9adIkTJo0yezfKr+JpkKrVq1Mmi2VgO1F5HbcuR8YLx5ERO6BGTCqEV9JREQkH3dpfgSc0wTpTpgBI9lJ6YhfXM9bxprITwknT2bBiGzD347MHPwUpLthAKZCanglkTs3QxIRKYNGpkmdGIARERERORgDMCIrsRmSyPUo4TejhHOHrNgEKQkDMKqVLSPiu3M/MCIiAgMwiRiAkdtiPzAiApSR/SKqigGYSqmhI749KaEpgRcVItehhHOG7IRGnkmlOA4YERERWU2I65PUdagVM2BkEfYDUyZmwYhqxt8IKRUDMHJr9uwH5pZNCkQkO7c9V7ATviRsgiQiIrfE7JedydGHS8V9wJgBUzF2xHcPvMgQEbkeBmBkMVftB8ZmSCJyFnc+R2iEPJNaMQAjcgPMghEZ42/CAdgHTBL2ASMiIiLrsQ+YJMyAqZzS+4GxGdJyvOMnuk4pvwWlnBtImRiAkVUc3Q+MiIgUik2QkvDKSORGlHLnT+QsSvkNqCL7xQBMEgZgpBpqaIYkIiLXwACMHNIPTAnDUaiFUjIARI7GY9/BmAGThAEYWc2WfmBqwCwYEQEqOhdUPAUpdVIpBmCkKvZshiQiIrIUAzByCWyGtA6bYkhteMw7HkfCl4YBmBsrtaJvl9L7gbkK1TQ9EJFZqjoHsA+YJO5/RbTSkiVL0KRJE/j6+iImJgZ79+6tsfzu3bsRExMDX19fNG3aFMuWLTMps379erRp0wZarRZt2rTBxo0b7VV9h3HlfmBqaYZkRoDUgsc6uSIGYJWkpKQgKSkJM2fOREZGBnr16oWhQ4ciMzPTbPnTp0/jtttuQ69evZCRkYEXXngBTz31FNavX28ok56ejpEjRyIhIQFHjhxBQkICRowYgR9++MFRm+U2XKUZUlV3wERkwN8+WUMjhFBxAtBY9+7d0aVLFyxdutQwr3Xr1rjrrrswb948k/LPPfccvvzySxw/ftwwLzExEUeOHEF6ejoAYOTIkSgoKMA333xjKDNkyBDUq1cPa9asMVuP4uJiFBcXG/5dUFCAyMhI5OfnIyjI8h/4za+/BQDwLrSsvI+F5QzlC2w7dLT5epuWAwBtXqnNy1bmlVcky3qq45FXYNf1W0Ofe9nZVSCyGyVlv+wRgH2bMceq8gUFBdDpdFZfL2z5jKj5r8DDV1qLgr6oCGefe9Gu9VUqZsBuKCkpwaFDhxAfH280Pz4+Hvv27TO7THp6ukn5wYMH4+DBgygtLa2xTHXrBIB58+ZBp9MZpsjISFs2yWrW9gNjMyQROZOSgi9V4jAUkjAAuyE3Nxfl5eUICwszmh8WFobs7Gyzy2RnZ5stX1ZWhtzc3BrLVLdOAJgxYwby8/MN07lz52zZJLfEZkjr8SJFZH9K+s2Ta/BydgWURqMxjsaFECbzaitfdb6169RqtdBqtRbXuTaldSxvhiQicgW8sVAAOZ5iVHEnKGbAbggJCYGnp6dJZionJ8ckg1UhPDzcbHkvLy/Ur1+/xjLVrdPZ1DQchb2bIZV0R8yLFZH9KOm37lAchkISZVwJFcDHxwcxMTFIS0szmp+Wloa4uDizy8TGxpqU37p1K7p27Qpvb+8ay1S3TlfjjH5grtIMSUT2wRsKcgcMwCqZOnUqPvjgA6xcuRLHjx/HlClTkJmZicTERADX+2aNHTvWUD4xMRFnz57F1KlTcfz4caxcuRIrVqzA9OnTDWWefvppbN26FfPnz8eJEycwf/58bNu2DUlJSY7ePDJDTZ3xedEiIjlxJHxp2AeskpEjR+LSpUuYO3cusrKy0K5dO6SmpiIqKgoAkJWVZTQmWJMmTZCamoopU6Zg8eLFaNiwId555x3ce++9hjJxcXFYu3YtXnzxRbz00kto1qwZUlJS0L17d4dum9L6gRXrPCQNR+Eq9PWCFDUkhUdIMIelIJemtBsJ1TY/AuwDJhEDsComTZqESZMmmf1bcnKyybw+ffrgp59+qnGd9913H+677z45qucQJXWsGxOsJEhj85hgtiqu5y3bmGBERESOxiZIUj01dcYHlJdBILKU0o5de/+2Fd9Fgp3wJWEApiLWvJzbEZTyNCQREVmPfcCk4RWQzHKFUfFd6WlIZsGIpFHaMav67BdJxj5gRLh+srP3+yGJiNyKHK8S4quIiJxDajMks2C2U1pGgag6SjtWlfZbdhr2AZOEAZjKWNMPzBWaIeXElD8RKYGrnIvYB0waBmDkdGrqjK+0O2elZRaIqlLaMaq03zC5LvVc+chtydkM6Sp3nkTknlzqHMQmSEkYgKkQmyGpMqVlGIgqKO3YZParCjmaHxmAETkXmyGJqDKlBV+O4FLZL5JMPVc9cmtshpRGjRc7IpKITZCSMABTKTZDOpcSs2AMwkgplHgsKvE363QMwCRhAEaKoaQxwdSYBSMi5+E5R30YgBE5iRLvqJWYeSB1UeIxqMTfqhJwHDBpGICpGJsha8Y7UiJyBJ5r1IkBGCmKkpohHUGJd9ZKzECQOijx2FPib5TcA1/GTURE5CQunf2SoxM9myBJrZTYDKmkLJgjTo5KvMNWYiaC3JsSjzkl/jaVhH3ApGEGjIiInEqJwRdZSMUBlFTMgJFdsTO+ZZR4p82LIqmZI36TLt38SJIxACO7NkPaSknNkGrGIIzsjceYC+NArJIwAHNjZUHlzq6C21BrFoxIjZj9sgz7gEnDAIzszh0646sZMxRkL0o8tngjRI7CAIwAKLMZUmmYBSMiObhD9gsAmyAlYgDm5tTeDMksmDyUmKkg16bEY4o3QNZhE6Q0DMDIIZzVDCk3NWfBlHjBJNek5mPJbbJfJJmyrm7kVO7aDOmKWTClBmFE7oq/ORs4qQlyyZIlaNKkCXx9fRETE4O9e/datNz3338PLy8vdOrUyfoPtQMGYCqglGZIZsFcn5ozFyQPNR9DbnfucEIAlpKSgqSkJMycORMZGRno1asXhg4diszMzBqXy8/Px9ixYzFgwADrPtCOlHVlI5fiSlkwV8Q7ciLH4G/N+QoKCoym4uJis+UWLFiACRMm4JFHHkHr1q2xcOFCREZGYunSpTWu/7HHHsPo0aMRGxtrj+rbhAGYSliaBbOmGdKVyN0M6XZ3slZQcwaDpFHzseOO5ww5O+FHRkZCp9MZpnnz5pl8XklJCQ4dOoT4+Hij+fHx8di3b1+19fzwww/xxx9/YNasWbJuv1R8FyQ5VEmQBj4F1jf6F+s8oM3X26FGyqavFwSPvAJnV4NIMqUGX8x+SSDHMBI3lj937hyCgv7/u9BqtSZFc3NzUV5ejrCwMKP5YWFhyM7ONrv6U6dO4fnnn8fevXvh5aWskEdZtSGXU1IH8Cl0di0sU1zPG9q8UtnWV1bPF155RbKtz5V4hARDn3vZ2dUgcgnumP0CIGsAFhQUZBSA1USjMe5PLIQwmQcA5eXlGD16NObMmYMWLVpIrKj8GICpSFlQObwKPGstV1oH8LZjUMUsmHWUmgVjEEaWYvaL5BASEgJPT0+TbFdOTo5JVgwACgsLcfDgQWRkZGDy5MkAAL1eDyEEvLy8sHXrVvTv398hdTeHfcBIMlfqjM++YETkSO58jnD0QKw+Pj6IiYlBWlqa0fy0tDTExcWZlA8KCsIvv/yCw4cPG6bExES0bNkShw8fRvfu3aXuAkmYASOz7J0FsxWzYMrCLBjVhtkvNyZjE6Slpk6dioSEBHTt2hWxsbFYvnw5MjMzkZiYCACYMWMG/vrrL6xatQoeHh5o166d0fKhoaHw9fU1me8MzIDdkJeXh4SEBMMTGAkJCbhy5Uq15UtLS/Hcc8+hffv2CAgIQMOGDTF27FhcuHDBqFzfvn2h0WiMplGjRtl5a6rn6mOCyYFZMCJyBJ4b5Ddy5EgsXLgQc+fORadOnbBnzx6kpqYiKioKAJCVlVXrmGBKoRFCqPhNTP9v6NChOH/+PJYvXw4AePTRRxEdHY2vvvrKbPn8/Hzcd999mDhxIjp27Ii8vDwkJSWhrKwMBw8eNJTr27cvWrRogblz5xrm+fn5QafTWVy3goIC6HQ65OfnW9xJEQCil/3P7HxL+oEB1mfArO2Mb0s/sApSs2BydsYH4LDO+ErMggFgFozMUnv2S44AbPuOGVaVt/V6YctntJ78Kjy10raxvLgIx999wa71VSo2QQI4fvw4tmzZgv379xvahN9//33Exsbi5MmTaNmypckyOp3OpB160aJF6NatGzIzM9G4cWPDfH9/f4SHh9t3I+xAqZ3x5eCqT0QqtSmSyFW4UvCleE5ognQnbIIEkJ6eDp1OZ9Qhr0ePHtDpdDUO7lZVfn4+NBoN6tatazR/9erVCAkJQdu2bTF9+nQUFtYc1RQXF5uMCiwnezVDulJnfJKXUjMd5Dw8JohqxgAMQHZ2NkJDQ03mh4aGVju4W1VFRUV4/vnnMXr0aKM06pgxY7BmzRrs2rULL730EtavX4977rmnxnXNmzfPaETgyMhI6zZIRvYeGd+Z74d01b5gSu08zAsuVVDqscDsl8yc9DJud+HWAdjs2bNNOsBXnSr6a5kbxK26wd2qKi0txahRo6DX67FkyRKjv02cOBEDBw5Eu3btMGrUKHz++efYtm0bfvrpp2rXN2PGDOTn5xumc+fOWbnltVNKZ3wiInJNGpkmtXLrPmCTJ0+u9YnD6Oho/Pzzz7h48aLJ3/7++2+zg7tVVlpaihEjRuD06dPYsWNHrZ0Iu3TpAm9vb5w6dQpdunQxW0ar1Zp9DYMrcOTI+HIMScG+YPLisBTE7JdKsl8kmVsHYCEhIQgJCam1XGxsLPLz8/Hjjz+iW7duAIAffvgB+fn5Zgd3q1ARfJ06dQo7d+5E/fr1a/2sX3/9FaWlpYiIiLB8Q5zMnTvjExGRjdgJXxK3boK0VOvWrTFkyBBMnDgR+/fvx/79+zFx4kQMHz7c6AnIVq1aYePGjQCAsrIy3HfffTh48CBWr16N8vJyZGdnIzs7GyUlJQCAP/74A3PnzsXBgwdx5swZpKam4v7770fnzp3Rs2dPp2xrZe7QGV+OvmByY18wZWZAyP6U+t0z+2Ufjh4J390o7+rlJKtXr0b79u0RHx+P+Ph4dOjQAR9//LFRmZMnTyI/Px8AcP78eXz55Zc4f/48OnXqhIiICMNU8eSkj48Ptm/fjsGDB6Nly5Z46qmnEB8fj23btsHT07KxuJRCqZ3x5SB3Z3wiUg6l3qhYQvHnJnbCl8StmyCtERwcjE8++aTGMpXHrI2OjkZtY9hGRkZi9+7dstSPqqfE1xOxLxj7gqmNUrNfjiJ39kvxwRdJxgyYyimpGZJZMNu48h0+uQelBl/8bTgAs182YwBGFrN3M6QUau4LplRKvSgTyU2t2S/2AZNGeVctko2mjmXDKzALdp2rnPTMUeqdPoMw96fU71ipvwmiCgzAyCrMgllH7VkwInen1uwXAHbCl0h5VyySlStmwZzJpU5+VSj1jl+pGRKSTqnfLYedcAw2QUrDAIwUx5nNkPag9pM0kSMp9UbEEq58A0jWYwCmApZmwSzl7s2Q9jgJcnBWZWZKyHZq/055YwU2QUrEAIwMlNQM6ewsmCvfiSo1CCOyN1c+9l3xnMMmSGkYgJFN3D0LZg9qv2NWe8bEnaj9u1T7b5nkocwrFcnOFTvjMwtmO1fOBBDZwpWPeZc917AJUhIGYGQzZsGsp/Y7Z7VnTtyB2r9Dtf+GjTAAk0SZVymyC7k741uLWTDHcuWMAJE1XPlYd+VzDPuAScMAjEzYqxnS0ZgFUya1Z1BcmRK/O0cGX2r/7ZK8lHmFIpdhbTMks2CO5cqZASIlsUfw5crnFgBsgpSIAZjKOLszvqMxC6bMIEyJmRSqmRK/MyUe22qiEUKWSa2UeXUil8IsmG3YnEHkGpj9IntgAKZCzIKpjxIzBUrMqJB5SvyulHhMqw6bICXhlYlk4YghKZgFIyJHY/arenwKUhoGYCrlikNSSMEsmDIzBkrMrJAxJX5HSjyWiazFqxLVyJpmSKVnwdT+om6AFy5yfa4+7IS7ZL8AsAlSIgZg5DSOzoKRMikxw0LX8buRl1sFX2ATpFQMwFTMHp3xmQWzDbNgRLVz9ewXUWUMwMipXDUL5upBmNIw06I8av5O2PRoITZBSsIATOVcdUgKZ2fBXB2zYORqeMwqD5sgpeGVyI0F1ilyyuc6YmBWJWAWTF5qzrgojZq/C2a/rMAMmCQMwIhZMJViRoFcBY9Vcke8Crk5ZsHsi1kweak586IUSvsOXL3jvdtmv25g86PtGIARAGbB1IqZBSKymRDyTCrFK5AKMAtmX66eBVNaEKa0DIyaKG3fM/tF7owBGLk8d82CqbkpksiR+FuzDZ+ClEa5Vx+SlSVZMCUMzMosmHMoLQtG5OrHpKufEyzCpyAlYQBGbkEJWTBXb4pUEqU1hamBWve5kpselZyhJ+n47ZIRZsHUy9UzDuQ+eCy6RvCl0cszqZXyv2GSjbM64zsKs2DuRa0ZGWdQ0r5mx3sXwiZISRiAkQlmwdSLmQdSCyXf2LhC9ouk47d8Q15eHhISEqDT6aDT6ZCQkIArV67UuMz48eOh0WiMph49ehiVKS4uxpNPPomQkBAEBATgjjvuwPnz5+24JTVjFqx6zIJdp6QgTEmZGXelpH2spGPPFqrKfoFPQUrFAOyG0aNH4/Dhw9iyZQu2bNmCw4cPIyEhodblhgwZgqysLMOUmppq9PekpCRs3LgRa9euxXfffYerV69i+PDhKC9X1oCmVak1C6bkO08l37ETuRIl/5aUfA4ywYFYJfFydgWU4Pjx49iyZQv279+P7t27AwDef/99xMbG4uTJk2jZsmW1y2q1WoSHh5v9W35+PlasWIGPP/4YAwcOBAB88skniIyMxLZt2zB48GD5N8YCgXWKcLVQuScgqUqCNPApcO6PurieN7R5lgWxSqSvFwSPvAJnV4NUhNkv1yNHBosZMJVLT0+HTqczBF8A0KNHD+h0Ouzbt6/GZXft2oXQ0FC0aNECEydORE5OjuFvhw4dQmlpKeLj4w3zGjZsiHbt2tW43uLiYhQUFBhNzsAsmPIo+c7dXpTUROZu1LhvldzxXsnnHpIfv20A2dnZCA0NNZkfGhqK7OzsapcbOnQoVq9ejR07duDNN9/EgQMH0L9/fxQXFxvW6+Pjg3r16hktFxYWVuN6582bZ+iLptPpEBkZaeOWVY99wezP1e+IXT0jQa6Dx5qL4lOQkrh1ADZ79myTTvJVp4MHDwIANBrTC7YQwuz8CiNHjsSwYcPQrl073H777fjmm2/w22+/YfPmzTXWq7b1zpgxA/n5+Ybp3LlzFm6x8zALVs16XLxDvlKoMVNjb2rcp8x+yYud8KVx6z5gkydPxqhRo2osEx0djZ9//hkXL140+dvff/+NsLAwiz8vIiICUVFROHXqFAAgPDwcJSUlyMvLM8qC5eTkIC4urtr1aLVaaLVaiz/XnjR1SiEKXTOTo4S+YK6OfcHI3hyV/VLjjQspm1sHYCEhIQgJCam1XGxsLPLz8/Hjjz+iW7duAIAffvgB+fn5NQZKVV26dAnnzp1DREQEACAmJgbe3t5IS0vDiBEjAABZWVk4evQoXn/9dRu2SF5yd8YvCyqHV4GnRWVL6wDehZavu6QO4GNFecNyEoKwYp0HtPnSh2m2R4f8snq+8Mpz72ZkIqVTc/YLgDxPMar4KUgX/dbl1bp1awwZMgQTJ07E/v37sX//fkycOBHDhw83egKyVatW2LhxIwDg6tWrmD59OtLT03HmzBns2rULt99+O0JCQnD33XcDAHQ6HSZMmIBp06Zh+/btyMjIwIMPPoj27dsbnop0BZZ2xndHSj4xOuqOXin9c9TYZGYvStmXzH4p+xxTGzZBSuO637zMVq9ejfbt2yM+Ph7x8fHo0KEDPv74Y6MyJ0+eRH5+PgDA09MTv/zyC+688060aNEC48aNQ4sWLZCeno46df6/09Jbb72Fu+66CyNGjEDPnj3h7++Pr776Cp6elmWKpAgJvFprGbk74yuxLxg75BMpj1ICe1vxN01SuXUTpDWCg4PxySef1FhGVEqV+vn54dtvv611vb6+vli0aBEWLVokuY7OpJS+YLY2RUohV1OkPTiqKZJ9wdyHUrJfjsKO93Ykx1OMzICRu3K3LJitmAWTTgkZC7UFD+5KCccSSccmSGkYgJHFlNIXjMNSGFNy/xYiZ2L2i5SMR4AKMAt2nRKyYK5OCZkLZsFsp4R9x473bkQv5JlUigEYWYVZMBnWwywYkUti9qsKjoQviZscBVQbZsGuc+csmNqGpSDXw+yXGwVfADSQoQ+YszfCidznSCCHYRZMhvW4eId8Z1NCU5qr4T6Thr9ZkhsDMBVhFuw6qVkwJd/BMgtGSuXK2S82PVajYiR8qZNKudnRQI7i6lkwJTRF8o6aiFwZh6GQhgGYyjALJg8lN0WqJQvGJjXLOXtfMfvlhtkvJ1qyZAmaNGkCX19fxMTEYO/evdWW3bBhAwYNGoQGDRogKCgIsbGxFg2i7gg8Ishm9sqCqekVRUQkDyV3vHdbTngKMiUlBUlJSZg5cyYyMjLQq1cvDB06FJmZmWbL79mzB4MGDUJqaioOHTqEfv364fbbb0dGRob12yszBmAqpPQsmKtgFsz5WTBSPlc+Rpj9qplGCFkmACgoKDCaiouLzX7mggULMGHCBDzyyCNo3bo1Fi5ciMjISCxdutRs+YULF+LZZ5/FLbfcgubNm+PVV19F8+bN8dVXX9ltv1jKPY8KchhmwZRLDRkBZzetuQI17CMlH+vuGnzJLTIyEjqdzjDNmzfPpExJSQkOHTqE+Ph4o/nx8fHYt2+fRZ+j1+tRWFiI4GDn/y74Mm6VCgm8ityrgTWWCaxThKuF8p3YyoLK4VXgKdv6qnLlF3UX1/OGNk8ZDzZYiy/qpuow++Xm9DcmqesAcO7cOQQF/f/xotVqTYrm5uaivLwcYWFhRvPDwsKQnZ1t0ce9+eabuHbtGkaMGGF7nWXCAIwk09QphSiU/2RVWgfwdkBAVRKkgU+B7Y/iKDkIK6vnC688eZuTiZTEnTveKz1DX7kJUco6ACAoKMgoAKtxGY3xfhFCmMwzZ82aNZg9eza++OILhIaGWl9ZmTE3qmLu2BeMTZHO4cxMhxqa2GzlzH3jytkvJeA5yVRISAg8PT1Nsl05OTkmWbGqUlJSMGHCBKxbtw4DBw60ZzUtxgCMZKGUvmDOwg75RI7nztkvl+DgpyB9fHwQExODtLQ0o/lpaWmIi4urdrk1a9Zg/Pjx+PTTTzFs2DDLP9DOVHCEqFfjwCu1lmEWrNJyvOOUhBkPqsBjQRqXORc5YST8qVOn4oMPPsDKlStx/PhxTJkyBZmZmUhMTAQAzJgxA2PHjjWUX7NmDcaOHYs333wTPXr0QHZ2NrKzs5Gfny/rrrAFAzA3Z0kQJhdmwZgFcxY2Q5py933C7JfzOWMk/JEjR2LhwoWYO3cuOnXqhD179iA1NRVRUVEAgKysLKMxwd577z2UlZXhiSeeQEREhGF6+umn5dwVNmEnfOITkZWXk9ghX+34RCQ5IvvlzsGXy2S/nGjSpEmYNGmS2b8lJycb/XvXrl32r5CN1BGmq5xSs2B8RVEN62EWjIiUji/jloQBGAGQty+Y2l/UrfYgzFn9f9y9yc0aztoXzH6pK/ul0cszqRUDMJVwZBbMGvbOgrFDPhG5Ap5z1IcBGBm4QhaMTZHSuXMWjJyH2S8VXk7ZBCmJCo8Y9XKHLJgtmAUjIiVz2XONg8cBczcMwMgIs2DyUXsWzBnYD8w5+4DZL15KyXo8alSGWTArl3PVO1MFYDMkKZlSgi9XPsdUvAtS6qRWDMDIBLNg8mEWjEg6HqsKxT5gkjAAUyFmwaxczo2HpbA3ZsEcy12bH+XG7BcpAQMwMstds2BsijTmjpkF9gNzL+54jAJuck4RAPQSJ/UmwBiAqZVas2DOwiwYuSN7f7fseK9s7AMmDY8gqhazYFWWc4c7VjPcNcNApFRucy4RkKEPmLM3wnkYgKmYXFkwS4MwS1mbBWOHfOnsHYQ5OgumxmZIR28zs19E0vAocmNN/HMlr8OSLJil7PWiblu4eod8InIOdryvhE9BSsKrgcopNQtmLVfJgsnFFbNg5D7UnP2iSqR2wK+YVIoBmJtjFqx6rp4Fc7ULCjvjk7MopenRrbJfJBkDsBvy8vKQkJAAnU4HnU6HhIQEXLlypcZlNBqN2emNN94wlOnbt6/J30eNGmXnrbGOo7Ng7JCvXO6UBVNTPzB32lZ3OgYrc8dzBp+ClIYB2A2jR4/G4cOHsWXLFmzZsgWHDx9GQkJCjctkZWUZTStXroRGo8G9995rVG7ixIlG5d577z17booJpWXBrMFhKWpZD7Ng5GCu9h0qJfvlltgHTBIvZ1dACY4fP44tW7Zg//796N69OwDg/fffR2xsLE6ePImWLVuaXS48PNzo31988QX69euHpk2bGs339/c3Kas0jQOvIPNqXcnrCaxThKuFtd/BauqUQhTKHzyU1gG8C61bpqQO4GPlMsD1O1qfAvc7eZTV84VXnnP79JE6MftVZTmV9W1VG4b0ANLT06HT6QzBFwD06NEDOp0O+/bts2gdFy9exObNmzFhwgSTv61evRohISFo27Ytpk+fjsLCmq/2xcXFKCgoMJqkUlMWzFU65DMLRq7G1b47Zr/sjBkwSXhUAcjOzkZoaKjJ/NDQUGRnZ1u0jo8++gh16tTBPffcYzR/zJgxWLNmDXbt2oWXXnoJ69evNylT1bx58wx90XQ6HSIjIy3fGAncpS+YLVy9Q77c3CUT4U59o6rjLtuo1GPOWR3vXSL7xQBMEmWe/WUye/bsajvKV0wHDx4EcL1DfVVCCLPzzVm5ciXGjBkDX1/jk8jEiRMxcOBAtGvXDqNGjcLnn3+Obdu24aeffqp2XTNmzEB+fr5hOnfunBVbXT1mwdyXq2XBiOxNCb8Jd+x4T/Jx6z5gkydPrvWJw+joaPz888+4ePGiyd/+/vtvhIWF1fo5e/fuxcmTJ5GSklJr2S5dusDb2xunTp1Cly5dzJbRarXQarW1rssWTfxzcfqfkGr/bklfsJDAq8i9GlhjGXv0BSsLKodXgadFZW3hrL5gxToPaPOVNxiOPfuC6esFwSNPetM6OY49mx/dNftlK5fIfgHXx/CSGmMq79TnMG4dgIWEhCAkpPpgo0JsbCzy8/Px448/olu3bgCAH374Afn5+YiLi6t1+RUrViAmJgYdO3asteyvv/6K0tJSRERE1L4BJInaOuQX1/OGNk85zbpEzuLK2S+XCb4AWYaR4DAUKte6dWsMGTIEEydOxP79+7F//35MnDgRw4cPN3oCslWrVti4caPRsgUFBfjss8/wyCOPmKz3jz/+wNy5c3Hw4EGcOXMGqampuP/++9G5c2f07NnT7ttVndqaIi3pC+asF3VzWArHsmdmwlEdut2lj5Q5jto2V8p+seO9A7EPmCQ8wm5YvXo12rdvj/j4eMTHx6NDhw74+OOPjcqcPHkS+fn5RvPWrl0LIQQeeOABk3X6+Phg+/btGDx4MFq2bImnnnoK8fHx2LZtGzw97deUpiTO7pCvtsFZlXDnT6R2ash+kXRu3QRpjeDgYHzyySc1lhFmIvVHH30Ujz76qNnykZGR2L17tyz1k5uj+oLZg7V9wWxpinQGufqCyd0UyXHByF6Y/XJxegFoJGaw9MyAEdmNs7NgtuCwFI7hauNKqZHaviMOO2EFNkFKoq6zvco005o+2VmZo/qC2QOHpaiZ3E2RSn1KjVyXUrNfUqgq+CLJGIBRjZQ8OKu9O+QzC+YYjsiwuGNHfEdsE7NfVDM5sl/MgJGbkpoFs4SzsmDWYod8aZgFI7kw+1VpOVfOfrEJUhIGYFQrNWfBnIV34kSOxd8cORqPOBVgFuz/MQumTGpr6nIF9vpOmP2qtJwrZ7+A608wyjGpFAMwsog7Dc7qKh3ylXhHzmZIZXDHPm3OpMTfmksQenkmleJRpxKOyILJjcNSyEMJGQIlYNDifuQ4tjnsBDkLAzCymKOzYNZgFsxx7JUFYzOkcrhK86OzqT74Yid8SZR3die7UVIWTAlNkdZy5WEpmAUjd6OE7JfqsQ+YJDz6yCpKHpzVWmrrkC8nd8tk0P9zheyXEm4oVJ/9ApgBk4gBmMowC+aa1JIFYzNk7diXTR7MfpGz8QhUIXd+RZG1mAWznStmwRi8OAezXzeWc6fsF3B9EHvJGTBnb4TzMAAjp1JCFkxNHfKVcPEi5VJLBpLZL5mwCVISHoUq5YpZMA5LoTz2yIKpJQhQC2a/bixnwznDVW4OyTYMwMhmah+WwhnNCcyCEUkj5TfE4KsKvV6eSaUYgLmxFj41Z7nYId853C0LZg/MgjmePfa50rJfbHqUGZsgJeHR6OZqC8KkYod82z5LShCmtCyYq3XGd+WO+K5cd1fH7BfJjQGYyjEL5pp4J09qoqbsl0sFX8yASeIaRyRJwiyY5VwlCyYHNWfByJjSmx+djcNOVIMj4UvCAIyYBavCVe5AXeWO3lbsB0YAs1/kvlzjqCTJ3DELxmEppGMWjOzBnY4FZr+qJ4RelkmtGIARANfMglmDw1KQq7NXB3wlZxpdNfulmo73QobmR/YBIzWw97AUcmbBlNAU6SjulAWTm72CAz5NaD/MfqkIO+FLwgCMZMUO+bZ9lrOHpZCLO118ybmcnf3isBNkb8o5c5NDcHBW9+yQL5WSs2BE7sylzzEcCV8SBmAkO3fKgtmCWTB5s2BK7qPk6uTet0rKgLpK9sulsQlSEuWctclhmAVjFswcZsFIKdRyLKrl3ELmMQBTKQ5LYV/MgiknE+IOXOGhASV958x+OYbQ62WZ1Eo5Z2yS3c3etgckrpgFswazYOYpNfNgj2ZIVwhq7EmpTbtKPQZrotqO92yClIQBmIq5UhZMCU2R1lJ7FozIGZyR/SKyBc/Wbs7ZWTB36pDvKoOzSiVXBkJJTVJkX3J918x+uRi+C1ISBmAqJzULpsSmSHfokM8sGNmTUpsfpXKF7JfbBF/AjSZEvcSJARi5MXtmwSzhTlkwWzALJg93DRpcnTtkv9jxnpyBARhxWApmwYhcnisc89aeO5T6+rQKQi9kmdRK+Uesg/z3v/9FXFwc/P39UbduXYuWEUJg9uzZaNiwIfz8/NC3b1/8+uuvRmWKi4vx5JNPIiQkBAEBAbjjjjtw/vx5O2xBzaRkweSg9GEp3LFDvlSu2B/HWq7wJKQr1FEqqceaOw47ofTgC4AMzY83JistWbIETZo0ga+vL2JiYrB3794ay+/evRsxMTHw9fVF06ZNsWzZMlu3WFYMwG4oKSnB/fffj8cff9ziZV5//XUsWLAA7777Lg4cOIDw8HAMGjQIhYWFhjJJSUnYuHEj1q5di++++w5Xr17F8OHDUV7u+B9XTUGYO2bB7MkV+nEoJSPAZkhlkXMf8kELy7nCOcNazsiApaSkICkpCTNnzkRGRgZ69eqFoUOHIjMz02z506dP47bbbkOvXr2QkZGBF154AU899RTWr18vxy6QRCOEinvAmZGcnIykpCRcuXKlxnJCCDRs2BBJSUl47rnnAFzPdoWFhWH+/Pl47LHHkJ+fjwYNGuDjjz/GyJEjAQAXLlxAZGQkUlNTMXjwYIvqVFBQAJ1Oh/z8fAQFWX7yvHThJpN5v5fWfKf5W0lYtX/7o7j6vwHA6X9Caq1T5tW6tZbJvRpYaxkAuFpo2clfFFp+d+1V4GlxWe/C2stU5WPDMgDgU2Dbz1SbL32QQ22e9OypV558QbNHXoFs6wIAfe5lWdcnN7kzYEoLwNSS/bImAKuc/TqTON2qz7H1emHLZ/TV3A0vjbTvr0yUYpfYaHF9u3fvji5dumDp0qWGea1bt8Zdd92FefPmmZR/7rnn8OWXX+L48eOGeYmJiThy5AjS09Ml1V0qL6d+ugs7ffo0srOzER8fb5in1WrRp08f7Nu3D4899hgOHTqE0tJSozINGzZEu3btsG/fvmoDsOLiYhQXFxv+nZ+fD+D6QW+NwkLTi28YivFnDUHYv6Vl1f6tIf7C6eLQGv6ehTP/1K+xTuHIwflruhrL1NUU49K12oMw/T+1FgEAiH8tzzbqiywPwIq9AS8rW03/9bEtCCsvsS0A+8cP8CmQFoR5lsnQfF1eXHsZC3no5VsXAOhFiazrk5vs2yvTd1FWVwuUSQ+sy0qltQaUl9gegJUXWx+AldQBYOUuLAsEYMWu0nv//z6x9rxfUd4RuZUyUWxTE6LROnD9/FJ1O7VaLbRardG8kpISHDp0CM8//7zR/Pj4eOzbt8/s+tPT042uwQAwePBgrFixAqWlpfD2dl5XCwZgNsrOzgYAhIUZZ4XCwsJw9uxZQxkfHx/Uq1fPpEzF8ubMmzcPc+bMMZkfGRkptdoWMJ/GJSInsfBGw2LnZF4f2ZVu6ks2LVdYWAidruabXVv5+PggPDwc32WnyrK+wMBAk+vbrFmzMHv2bKN5ubm5KC8vN3vdre6amp2dbbZ8WVkZcnNzERERIX0DbOTWAdjs2bPNBjKVHThwAF27drX5MzQa4zsoIYTJvKpqKzNjxgxMnTrV8G+9Xo/Lly+jfv36ta67QkFBASIjI3Hu3Dm7paHtwVXrDbhu3Vlvx2K9Hc9V625LvYUQKCwsRMOGDe1WL19fX5w+fRolJfJkj81dE6tmvyqz9rprrry5+Y7m1gHY5MmTMWrUqBrLREdH27Tu8PBwANej68oRdE5OjiHaDg8PR0lJCfLy8oyyYDk5OYiLi6t23eZSr5Y+mVlVUFCQS51wKrhqvQHXrTvr7Vist+O5at2trbe9Ml+V+fr6wtfXsQ9hhISEwNPT0yTbVfm6W1V4eLjZ8l5eXqhfv+YuM/amjMek7CQkJAStWrWqcbL1AGrSpAnCw8ORlpZmmFdSUoLdu3cbgquYmBh4e3sblcnKysLRo0drDMCIiIjImI+PD2JiYoyuqQCQlpZW7TU1NjbWpPzWrVvRtWtXp/b/Atw8ALNGZmYmDh8+jMzMTJSXl+Pw4cM4fPgwrl79/17WrVq1wsaNGwFcT10mJSXh1VdfxcaNG3H06FGMHz8e/v7+GD16NIDrdyETJkzAtGnTsH37dmRkZODBBx9E+/btMXDgQKdsJxERkauaOnUqPvjgA6xcuRLHjx/HlClTkJmZicTERADXu/CMHTvWUD4xMRFnz57F1KlTcfz4caxcuRIrVqzA9OnWPV1qF4KEEEKMGzdOADCZdu7caSgDQHz44YeGf+v1ejFr1iwRHh4utFqt6N27t/jll1+M1vvvv/+KyZMni+DgYOHn5yeGDx8uMjMz7b49RUVFYtasWaKoqMjunyUnV623EK5bd9bbsVhvx3PVurtqve1t8eLFIioqSvj4+IguXbqI3bt3G/42btw40adPH6Pyu3btEp07dxY+Pj4iOjpaLF261ME1No/jgBERERE5GJsgiYiIiByMARgRERGRgzEAIyIiInIwBmBEREREDsYAjIiIiMjBGIC5qP/+97+Ii4uDv7+/xaPkCyEwe/ZsNGzYEH5+fujbty9+/fVXozLFxcV48sknERISgoCAANxxxx04f/68rHXPy8tDQkICdDoddDodEhIScOXKlRqX0Wg0Zqc33njDUKZv374mf6/tTQj2rvf48eNN6tSjRw+jMvbe59bWu7S0FM899xzat2+PgIAANGzYEGPHjsWFCxeMysm9v5csWYImTZrA19cXMTEx2Lt3b43ld+/ejZiYGPj6+qJp06ZYtmyZSZn169ejTZs20Gq1aNOmjWEcP7lZU/cNGzZg0KBBaNCgAYKCghAbG4tvv/3WqExycrLZ472oSPrLr22t965du8zW6cSJE0blHLHPram3ud+gRqNB27ZtDWUcsb/37NmD22+/HQ0bNoRGo8GmTZtqXUZJxzjZgZOHwSAbvfzyy2LBggVi6tSpQqfTWbTMa6+9JurUqSPWr18vfvnlFzFy5EgREREhCgoKDGUSExNFo0aNRFpamvjpp59Ev379RMeOHUVZWZlsdR8yZIho166d2Ldvn9i3b59o166dGD58eI3LZGVlGU0rV64UGo1G/PHHH4Yyffr0ERMnTjQqd+XKFafWe9y4cWLIkCFGdbp06ZJRGXvvc2vrfeXKFTFw4ECRkpIiTpw4IdLT00X37t1FTEyMUTk59/fatWuFt7e3eP/998WxY8fE008/LQICAsTZs2fNlv/zzz+Fv7+/ePrpp8WxY8fE+++/L7y9vcXnn39uKLNv3z7h6ekpXn31VXH8+HHx6quvCi8vL7F//36b6ihX3Z9++mkxf/588eOPP4rffvtNzJgxQ3h7e4uffvrJUObDDz8UQUFBJse9M+u9c+dOAUCcPHnSqE6Vj1NH7HNr633lyhWj+p47d04EBweLWbNmGco4Yn+npqaKmTNnivXr1wsAYuPGjTWWV9IxTvbBAMzFffjhhxYFYHq9XoSHh4vXXnvNMK+oqEjodDqxbNkyIcT1E5W3t7dYu3atocxff/0lPDw8xJYtW2Sp77FjxwQAoxNEenq6ACBOnDhh8XruvPNO0b9/f6N5ffr0EU8//bQs9azK1nqPGzdO3HnnndX+3d77XK79/eOPPwoARhc5Ofd3t27dRGJiotG8Vq1aieeff95s+WeffVa0atXKaN5jjz0mevToYfj3iBEjxJAhQ4zKDB48WIwaNUqWOlewtu7mtGnTRsyZM8fwb0t/11JYW++KACwvL6/adTpin0vd3xs3bhQajUacOXPGMM8R+7sySwIwJR3jZB9sglSJ06dPIzs7G/Hx8YZ5Wq0Wffr0wb59+wAAhw4dQmlpqVGZhg0bol27doYyUqWnp0On06F79+6GeT169IBOp7P4My5evIjNmzdjwoQJJn9bvXo1QkJC0LZtW0yfPh2FhYVOr/euXbsQGhqKFi1aYOLEicjJyTH8zd77XI79DQD5+fnQaDQmzd1y7O+SkhIcOnTIaB8AQHx8fLV1TE9PNyk/ePBgHDx4EKWlpTWWketYtrXuVen1ehQWFiI4ONho/tWrVxEVFYWbbroJw4cPR0ZGhiLq3blzZ0RERGDAgAHYuXOn0d/svc/l2N8rVqzAwIEDERUVZTTfnvvbFko5xsl+vJxdAXKMirfBV31jfFhYGM6ePWso4+Pjg3r16pmUqfo2eSn1CA0NNZkfGhpq8Wd89NFHqFOnDu655x6j+WPGjDG8JP3o0aOYMWMGjhw5YvIiVkfWe+jQobj//vsRFRWF06dP46WXXkL//v1x6NAhaLVau+9zOfZ3UVERnn/+eYwePRpBQUGG+XLt79zcXJSXl5s9NqurY3Z2ttnyZWVlyM3NRURERLVl5DqWba17VW+++SauXbuGESNGGOa1atUKycnJaN++PQoKCvD222+jZ8+eOHLkCJo3b+6UekdERGD58uWIiYlBcXExPv74YwwYMAC7du1C7969AVT/vci1z6Xu76ysLHzzzTf49NNPjebbe3/bQinHONkPAzAFmT17NubMmVNjmQMHDqBr1642f4ZGozH6txDCZF5VlpSxtO7m6mDpZ1RYuXIlxowZA19fX6P5EydONPx/u3bt0Lx5c3Tt2hU//fQTunTp4pR6jxw50qhOXbt2RVRUFDZv3mwSQFqzXkft79LSUowaNQp6vR5Lliwx+pst+7sm1h6b5spXnW/L8W4LWz9nzZo1mD17Nr744gujQLlHjx5GD2v07NkTXbp0waJFi/DOO+84pd4tW7ZEy5YtDf+OjY3FuXPn8L///c8QgFm7TlvZ+hnJycmoW7cu7rrrLqP5jtrf1lLSMU7yYwCmIJMnT671KbLo6Gib1h0eHg7g+l1VRESEYX5OTo7hDio8PBwlJSXIy8szysjk5OQgLi5Olrr//PPPuHjxosnf/v77b5M7OXP27t2LkydPIiUlpdayXbp0gbe3N06dOlVtQOCoeleIiIhAVFQUTp06BcD2fe6IepeWlmLEiBE4ffo0duzYYZT9MseS/W1OSEgIPD09Te7aKx+bVYWHh5st7+Xlhfr169dYxprvyx51r5CSkoIJEybgs88+w8CBA2ss6+HhgVtuucVw3Eglpd6V9ejRA5988onh3/be51LqLYTAypUrkZCQAB8fnxrLyr2/baGUY5zsyPHdzkhO1nbCnz9/vmFecXGx2U74KSkphjIXLlywSyf8H374wTBv//79FncKHzdunMnTeNX55ZdfBACxe/dum+tbQWq9K+Tm5gqtVis++ugjIYT997mt9S4pKRF33XWXaNu2rcjJybHos6Ts727duonHH3/caF7r1q1r7ITfunVro3mJiYkmHZSHDh1qVGbIkCF26YRvTd2FEOLTTz8Vvr6+tXbErqDX60XXrl3FQw89JKWqRmypd1X33nuv6Nevn+Hfjtjntta74iGCX375pdbPsMf+rgwWdsJXyjFO9sEAzEWdPXtWZGRkiDlz5ojAwECRkZEhMjIyRGFhoaFMy5YtxYYNGwz/fu2114ROpxMbNmwQv/zyi3jggQfMDkNx0003iW3btomffvpJ9O/f3y7DUHTo0EGkp6eL9PR00b59e5NhEarWXQgh8vPzhb+/v1i6dKnJOn///XcxZ84cceDAAXH69GmxefNm0apVK9G5c2dZh3Owpt6FhYVi2rRpYt++feL06dNi586dIjY2VjRq1Mih+9zaepeWloo77rhD3HTTTeLw4cNGj+UXFxcLIeTf3xVDC6xYsUIcO3ZMJCUliYCAAMOTas8//7xISEgwlK94RH/KlCni2LFjYsWKFSaP6H///ffC09NTvPbaa+L48ePitddes+swFJbW/dNPPxVeXl5i8eLF1Q7hMXv2bLFlyxbxxx9/iIyMDPHQQw8JLy8vo0Da0fV+6623xMaNG8Vvv/0mjh49Kp5//nkBQKxfv95QxhH73Np6V3jwwQdF9+7dza7TEfu7sLDQcJ4GIBYsWCAyMjIMTxYr+Rgn+2AA5qLGjRsnAJhMO3fuNJQBID788EPDv/V6vZg1a5YIDw8XWq1W9O7d2+Ru8N9//xWTJ08WwcHBws/PTwwfPlxkZmbKWvdLly6JMWPGiDp16og6deqIMWPGmDzaXrXuQgjx3nvvCT8/P7NjTWVmZorevXuL4OBg4ePjI5o1ayaeeuopkzG3HFnvf/75R8THx4sGDRoIb29v0bhxYzFu3DiT/WnvfW5tvU+fPm322Kp8fNljfy9evFhERUUJHx8f0aVLF6NM2rhx40SfPn2Myu/atUt07txZ+Pj4iOjoaLOB+WeffSZatmwpvL29RatWrYyCBTlZU/c+ffqY3bfjxo0zlElKShKNGzcWPj4+okGDBiI+Pl7s27fPqfWeP3++aNasmfD19RX16tUTt956q9i8ebPJOh2xz609Vq5cuSL8/PzE8uXLza7PEfu7IgNX3feu9GOc5KcR4kavPiIiIiJyCI4DRkRERORgDMCIiIiIHIwBGBEREZGDMQAjIiIicjAGYEREREQOxgCMiIiIyMEYgBERERE5GAMwIiIiIgdjAEZERETkYAzAiIiIiByMARgRERGRg/0frqXPPb6J4kMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Absolute difference between U3 and its second order taylor expansion\")\n",
    "plt.contourf(*X, np.abs(U3(x_ad2).as_func(H) - U3(X)), levels=20);\n",
    "plt.axis('equal'); plt.scatter(*x,color='red'); plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Sparse differentiation\n",
    "\n",
    "When differentiating a high dimensional function, for instance in the context of a PDE discretization scheme, sparse automatic differentiation becomes for reasons of memory and computation cost. We check here the consistency of the related Taylor expansions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:37.122665Z",
     "iopub.status.busy": "2024-04-30T09:22:37.122559Z",
     "iopub.status.idle": "2024-04-30T09:22:37.124402Z",
     "shell.execute_reply": "2024-04-30T09:22:37.124168Z"
    }
   },
   "outputs": [],
   "source": [
    "x_sp  = ad.Sparse.identity(constant=x)\n",
    "x_sp2 = ad.Sparse2.identity(constant=x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:37.125773Z",
     "iopub.status.busy": "2024-04-30T09:22:37.125679Z",
     "iopub.status.idle": "2024-04-30T09:22:37.151366Z",
     "shell.execute_reply": "2024-04-30T09:22:37.150998Z"
    }
   },
   "outputs": [],
   "source": [
    "assert np.allclose(U3(x_sp ).as_func(H), U3(x_ad ).as_func(H))\n",
    "assert np.allclose(U3(x_sp2).as_func(H), U3(x_ad2).as_func(H))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sparse classes also provide tangent, adjoint, and hessian linear operators, stored as opaque sparse matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:37.153392Z",
     "iopub.status.busy": "2024-04-30T09:22:37.153248Z",
     "iopub.status.idle": "2024-04-30T09:22:37.159129Z",
     "shell.execute_reply": "2024-04-30T09:22:37.158898Z"
    }
   },
   "outputs": [],
   "source": [
    "tangent_op = U3(x_sp ).tangent_operator()\n",
    "hessian_op = U3(x_sp2).hessian_operator()\n",
    "H_ = H.reshape(2,-1) # Depth must be at most two"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:37.160552Z",
     "iopub.status.busy": "2024-04-30T09:22:37.160463Z",
     "iopub.status.idle": "2024-04-30T09:22:37.172981Z",
     "shell.execute_reply": "2024-04-30T09:22:37.172672Z"
    }
   },
   "outputs": [],
   "source": [
    "assert np.allclose(U3(x_sp ).as_func(H_), U3(x) + tangent_op*H_)\n",
    "assert np.allclose(U3(x_sp2).as_func(H_), U3(x) + tangent_op*H_ + 0.5*lp.dot_VV(H_,hessian_op*H_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we check the adjoint operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:37.174634Z",
     "iopub.status.busy": "2024-04-30T09:22:37.174547Z",
     "iopub.status.idle": "2024-04-30T09:22:37.177857Z",
     "shell.execute_reply": "2024-04-30T09:22:37.177584Z"
    }
   },
   "outputs": [],
   "source": [
    "adjoint_op = U3(x_sp).adjoint_operator()\n",
    "\n",
    "np.random.seed(42)\n",
    "R_ = np.random.rand(1,H_.shape[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T09:22:37.179222Z",
     "iopub.status.busy": "2024-04-30T09:22:37.179139Z",
     "iopub.status.idle": "2024-04-30T09:22:37.181167Z",
     "shell.execute_reply": "2024-04-30T09:22:37.180924Z"
    }
   },
   "outputs": [],
   "source": [
    "assert np.allclose(lp.dot_VV(R_,tangent_op*H_), lp.dot_VV(H_,adjoint_op*R_))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}