{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'm currently taking part in a MOOC called [Computers, Waves and Simulations](https://www.coursera.org/learn/computers-waves-simulations). It is about numerical methods for wave equations. So far, it seems really well designed. In week 2, the content focuses on finite difference methods. \n",
    "\n",
    "In particular, the course introduces a very general method for deriving high-order stencils of derivatives by solving a linear system. In this post, I try to reproduce the method using a symbolic approach, while the class example focused on a numerical approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will start by solving the system with numbers and then move on to solve it symbolically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deriving a fourth order accurate approximation of a second derivative "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the exercices in the online class is to derive a fourth order accurate expression for the second derivative, using the value of the function at points: $x-2dx, x-dx, x, x+dx, +2dx$. After writing the Taylor expansion at each of these points, they are combined to form an approximation for the derivative.\n",
    "\n",
    "Let's follow along and start with the matrix obtained in the exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 1\\\\2 & 1 & 0 & -1 & -2\\\\4 & 1 & 0 & 1 & 4\\\\8 & 1 & 0 & -1 & -8\\\\16 & 1 & 0 & 1 & 16\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[ 1, 1, 1,  1,  1],\n",
       "[ 2, 1, 0, -1, -2],\n",
       "[ 4, 1, 0,  1,  4],\n",
       "[ 8, 1, 0, -1, -8],\n",
       "[16, 1, 0,  1, 16]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import *\n",
    "\n",
    "weight_matrix = Matrix([[1,  1, 1,  1,  1],\n",
    "                 [2,  1, 0, -1, -2],\n",
    "                 [4,  1, 0,  1,  4],\n",
    "                 [8,  1, 0, -1, -8],\n",
    "                 [16, 1, 0,  1, 16]])\n",
    "\n",
    "weight_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does this matrix mean? Each column represents the Taylor expansion at one of the points as a function of $\\frac{dx^i}\n",
    "{i!}f^{(i)}(x)$ (each row coefficient).\n",
    "\n",
    "This matrix can then be inverted:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & - \\frac{1}{12} & - \\frac{1}{24} & \\frac{1}{12} & \\frac{1}{24}\\\\0 & \\frac{2}{3} & \\frac{2}{3} & - \\frac{1}{6} & - \\frac{1}{6}\\\\1 & 0 & - \\frac{5}{4} & 0 & \\frac{1}{4}\\\\0 & - \\frac{2}{3} & \\frac{2}{3} & \\frac{1}{6} & - \\frac{1}{6}\\\\0 & \\frac{1}{12} & - \\frac{1}{24} & - \\frac{1}{12} & \\frac{1}{24}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0, -1/12, -1/24,  1/12, 1/24],\n",
       "[0,   2/3,   2/3,  -1/6, -1/6],\n",
       "[1,     0,  -5/4,     0,  1/4],\n",
       "[0,  -2/3,   2/3,   1/6, -1/6],\n",
       "[0,  1/12, -1/24, -1/12, 1/24]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight_matrixinv = weight_matrix.inv()\n",
    "\n",
    "weight_matrixinv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can easily solve for interesting solutions. Setting the first coefficient to one in the solution vector should give us a fourth order accurate approximation to $f(x)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\\\1\\\\0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[0],\n",
       "[1],\n",
       "[0],\n",
       "[0]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight_matrixinv.multiply(Matrix([1, 0, 0, 0, 0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recalling that the above are coefficients for the sampled points, we can go further and multiply the inverted solution by the function values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & f{\\left(x \\right)} & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[0, 0, f(x), 0, 0]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, h = symbols('x, h')\n",
    "f = Function('f')\n",
    "\n",
    "approx_points = Matrix([f(x+2*h), f(x+h), f(x), f(x-h), f(x-2*h)])\n",
    "\n",
    "weight_matrixinv.multiply(Matrix([1, 0, 0, 0, 0])).multiply_elementwise(approx_points).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can apply the same reasoning to obtain fourth order versions of the other derivatives.\n",
    "\n",
    "First derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{f{\\left(2 h + x \\right)}}{12 h} & \\frac{2 f{\\left(h + x \\right)}}{3 h} & 0 & - \\frac{2 f{\\left(- h + x \\right)}}{3 h} & \\frac{f{\\left(- 2 h + x \\right)}}{12 h}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-f(2*h + x)/(12*h), 2*f(h + x)/(3*h), 0, -2*f(-h + x)/(3*h), f(-2*h + x)/(12*h)]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight_matrixinv.multiply(Matrix([0, 1/h, 0, 0, 0])).multiply_elementwise(approx_points).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{f{\\left(2 h + x \\right)}}{12 h^{2}} & \\frac{4 f{\\left(h + x \\right)}}{3 h^{2}} & - \\frac{5 f{\\left(x \\right)}}{2 h^{2}} & \\frac{4 f{\\left(- h + x \\right)}}{3 h^{2}} & - \\frac{f{\\left(- 2 h + x \\right)}}{12 h^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-f(2*h + x)/(12*h**2), 4*f(h + x)/(3*h**2), -5*f(x)/(2*h**2), 4*f(-h + x)/(3*h**2), -f(-2*h + x)/(12*h**2)]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight_matrixinv.multiply(Matrix([0, 0, factorial(2)/h**2, 0, 0])).multiply_elementwise(approx_points).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Third derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{f{\\left(2 h + x \\right)}}{2 h^{3}} & - \\frac{f{\\left(h + x \\right)}}{h^{3}} & 0 & \\frac{f{\\left(- h + x \\right)}}{h^{3}} & - \\frac{f{\\left(- 2 h + x \\right)}}{2 h^{3}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[f(2*h + x)/(2*h**3), -f(h + x)/h**3, 0, f(-h + x)/h**3, -f(-2*h + x)/(2*h**3)]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight_matrixinv.multiply(Matrix([0, 0, 0, factorial(3)/h**3, 0])).multiply_elementwise(approx_points).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fourth derivative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{f{\\left(2 h + x \\right)}}{h^{4}} & - \\frac{4 f{\\left(h + x \\right)}}{h^{4}} & \\frac{6 f{\\left(x \\right)}}{h^{4}} & - \\frac{4 f{\\left(- h + x \\right)}}{h^{4}} & \\frac{f{\\left(- 2 h + x \\right)}}{h^{4}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[f(2*h + x)/h**4, -4*f(h + x)/h**4, 6*f(x)/h**4, -4*f(-h + x)/h**4, f(-2*h + x)/h**4]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight_matrixinv.multiply(Matrix([0, 0, 0, 0, factorial(4)/h**4])).multiply_elementwise(approx_points).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deriving the stencil with a formal Taylor expansion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you may have noticed, we started with a matrix of integers and derived the stencils from there. Can we obtain the same matrix, but keep it symbolic? We can! Let's do this using `sympy`.\n",
    "\n",
    "(as a side note, I found this discussion helpful: https://stackoverflow.com/questions/16869587/how-to-do-a-symbolic-taylor-expansion-of-an-unknown-function-fx-using-sympy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first write a function that generates a Taylor expansion up to order $n$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, h = symbols('x, h')\n",
    "f = Function('f')\n",
    "\n",
    "def TaylorExpansion(point=h, order=4):\n",
    "    return sum(point**i/factorial(i) * f(x).diff(x, i) for i in range(order+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{h^{2} \\frac{d^{2}}{d x^{2}} f{\\left(x \\right)}}{2} + h \\frac{d}{d x} f{\\left(x \\right)} + f{\\left(x \\right)}$"
      ],
      "text/plain": [
       "h**2*Derivative(f(x), (x, 2))/2 + h*Derivative(f(x), x) + f(x)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TaylorExpansion(h, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{h^{5} \\frac{d^{5}}{d x^{5}} f{\\left(x \\right)}}{120} + \\frac{h^{4} \\frac{d^{4}}{d x^{4}} f{\\left(x \\right)}}{24} + \\frac{h^{3} \\frac{d^{3}}{d x^{3}} f{\\left(x \\right)}}{6} + \\frac{h^{2} \\frac{d^{2}}{d x^{2}} f{\\left(x \\right)}}{2} + h \\frac{d}{d x} f{\\left(x \\right)} + f{\\left(x \\right)}$"
      ],
      "text/plain": [
       "h**5*Derivative(f(x), (x, 5))/120 + h**4*Derivative(f(x), (x, 4))/24 + h**3*Derivative(f(x), (x, 3))/6 + h**2*Derivative(f(x), (x, 2))/2 + h*Derivative(f(x), x) + f(x)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TaylorExpansion(h, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nice thing about this is that we can use the same function to generate an expansion at a point away, let's say $x + 2h$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{2 h^{4} \\frac{d^{4}}{d x^{4}} f{\\left(x \\right)}}{3} + \\frac{4 h^{3} \\frac{d^{3}}{d x^{3}} f{\\left(x \\right)}}{3} + 2 h^{2} \\frac{d^{2}}{d x^{2}} f{\\left(x \\right)} + 2 h \\frac{d}{d x} f{\\left(x \\right)} + f{\\left(x \\right)}$"
      ],
      "text/plain": [
       "2*h**4*Derivative(f(x), (x, 4))/3 + 4*h**3*Derivative(f(x), (x, 3))/3 + 2*h**2*Derivative(f(x), (x, 2)) + 2*h*Derivative(f(x), x) + f(x)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expansion = TaylorExpansion(2*h, 4)\n",
    "expansion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, to fill the matrix, all we have to do is take the coefficient that's before a derivative. We can use the `.coeff` method for that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{4 h^{3}}{3}$"
      ],
      "text/plain": [
       "4*h**3/3"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expansion.coeff(f(x).diff(x, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the above ideas, we can now:\n",
    "\n",
    "- generate a Taylor expansion at each point of the grid\n",
    "- put the coefficient for each derivative of that expansion in a matrix\n",
    "- invert that matrix to find accurate stencils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's write a prototype for order 4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 1\\\\- 2 h & - h & 0 & h & 2 h\\\\2 h^{2} & \\frac{h^{2}}{2} & 0 & \\frac{h^{2}}{2} & 2 h^{2}\\\\- \\frac{4 h^{3}}{3} & - \\frac{h^{3}}{6} & 0 & \\frac{h^{3}}{6} & \\frac{4 h^{3}}{3}\\\\\\frac{2 h^{4}}{3} & \\frac{h^{4}}{24} & 0 & \\frac{h^{4}}{24} & \\frac{2 h^{4}}{3}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[        1,       1, 1,       1,        1],\n",
       "[     -2*h,      -h, 0,       h,      2*h],\n",
       "[   2*h**2,  h**2/2, 0,  h**2/2,   2*h**2],\n",
       "[-4*h**3/3, -h**3/6, 0,  h**3/6, 4*h**3/3],\n",
       "[ 2*h**4/3, h**4/24, 0, h**4/24, 2*h**4/3]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "n_points = 5\n",
    "grid_points = np.arange(-(n_points-1)/2, (n_points-1)/2 + 1).astype(int)\n",
    "order = 4\n",
    "\n",
    "coef_matrix = ZeroMatrix(n_points, n_points).as_mutable()\n",
    "\n",
    "for p, h_coef in zip(range(n_points), grid_points):\n",
    "    \n",
    "    expansion = TaylorExpansion(h_coef * h, order)\n",
    "    \n",
    "    for derivative in range(order + 1):\n",
    "        term =  f(x).diff(x, derivative)\n",
    "        coef_matrix[derivative, p] = expansion.coeff(term)\n",
    "\n",
    "coef_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can just solve for the second derivative by inverting that matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{1}{12 h^{2}} & \\frac{4}{3 h^{2}} & - \\frac{5}{2 h^{2}} & \\frac{4}{3 h^{2}} & - \\frac{1}{12 h^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-1/(12*h**2), 4/(3*h**2), -5/(2*h**2), 4/(3*h**2), -1/(12*h**2)]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(coef_matrix.inv() @ Matrix([0, 0, 1, 0, 0])).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is exactly the expression we had before. Our little algorithm works!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stencils for the second derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's extend that procedure by adding even higher order terms and let's focus on the second derivative stencil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_second_derivative_stencil(n_points, order=4):\n",
    "    \"\"\"Returns a `order` accurate stencil for the second derivative.\"\"\"\n",
    "\n",
    "    grid_points = np.arange(-(n_points-1)/2, (n_points-1)/2 + 1).astype(int)\n",
    "\n",
    "    coef_matrix = ZeroMatrix(n_points, n_points).as_mutable()\n",
    "\n",
    "    for p, h_coef in zip(range(n_points), grid_points):\n",
    "\n",
    "        expansion = TaylorExpansion(h_coef * h, order)\n",
    "\n",
    "        for derivative in range(order + 1):\n",
    "            term =  f(x).diff(x, derivative)\n",
    "            coef_matrix[derivative, p] = expansion.coeff(term)\n",
    "    \n",
    "    second_derivative_vector = ZeroMatrix(order + 1, 1).as_mutable()\n",
    "    second_derivative_vector[2, 0] = 1\n",
    "\n",
    "    return coef_matrix.inv() @ second_derivative_vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{1}{12 h^{2}} & \\frac{4}{3 h^{2}} & - \\frac{5}{2 h^{2}} & \\frac{4}{3 h^{2}} & - \\frac{1}{12 h^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-1/(12*h**2), 4/(3*h**2), -5/(2*h**2), 4/(3*h**2), -1/(12*h**2)]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_second_derivative_stencil(n_points=5).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{1}{90 h^{2}} & - \\frac{3}{20 h^{2}} & \\frac{3}{2 h^{2}} & - \\frac{49}{18 h^{2}} & \\frac{3}{2 h^{2}} & - \\frac{3}{20 h^{2}} & \\frac{1}{90 h^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[1/(90*h**2), -3/(20*h**2), 3/(2*h**2), -49/(18*h**2), 3/(2*h**2), -3/(20*h**2), 1/(90*h**2)]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_second_derivative_stencil(n_points=7, order=6).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now do the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "stencils = []\n",
    "for n_points in [3, 5, 7, 9, 11, 13]:\n",
    "    x_grid = np.arange(-(n_points-1)/2, (n_points-1)/2 + 1).astype(int)\n",
    "    stencil = compute_second_derivative_stencil(n_points, order=n_points-1)\n",
    "    stencils.append((x_grid, stencil))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('bmh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAFgCAYAAACmDI9oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOydeXxcVd24n29mJvvSbE3SvaWFsm8ti7IIBQFBUXBBXlAQBBVUXvHnLqLoKy6v4vu6vCJQRERc2BFFLFBAodCWrXvTJU2afWn2ZLbz+2Mm7WQyk8xyZ7np9/l87ieZuXOfc+65Z745Offcc8QYg6IoiqIoiqIoAXIynQFFURRFURRFySa0gawoiqIoiqIoIWgDWVEURVEURVFC0AayoiiKoiiKooSgDWRFURRFURRFCUEbyIqiKIqiKIoSgjaQlWmJiCwQESMip2U6L+lGRG4VkfqQ11eJiDeTeVIUZfqhcVbj7HRGG8hKyhCRK0VknYj0iMiwiGwWkZtFRNKQfCNQB6yJ5yAR8YrIVckkLCJzgn803pWMJwl+DJySobQVRUkjwYaaibAtTkPyGmeVaYsz0xlQpjXtwG3AVmAUOB34JeAFfpbKhI0xPqA1lWlkK8aYAWAg0/lQFCVt7AZODXuvI9WJapzVODud0R5kJWUYY542xjxqjNlsjNlpjPkt8A/gXZMdF+wV+LyIPCQigyLSLCJfCPtMnYg8KCL7gr3Tz4vIspD94279hbz+sIg8ISJDIrJTRK4MOWY34ABWjvXCTJLH00TkXyLSH9zeFJHzgrsbgz+fC3p2hxx3bvC4YRHZKyIrRaQyZP+9IvJPEblORBpEpE9EHhOR6rD0zxGRF4Pn0Ssiq0XkkOC+cbf+FEWZ9viMMa1hm2+yAzTOapxVJkcbyEpakAAnAe8EnovhkG8BzwPHAz8Afigil4y5gEeBpcBFwElAG/CMiFRN4b0d+B1wDPAnAkF6SXDfcsAH3ETgtmFdlHNxAI8TuK14QnC7FRgKfuSE4M9Lg47lwePOBh4DHgym/35gAfBI8JzGWA6cBVwInA8cR+B23lj65wBPA+sI9BqdDNwHuKY4d0VRpidzRKQpuP1NRN4R43EaZzXOKtEwxuimW8o2oIzAbSg3gaB4SwzHGOB3Ye89ALwU/H1F8DNHhOzPA1rG/AQCogFOC3v9hZBjnMG8XR/ynhe4aor8lQdd74qyf06k/QT+EN0e9t684GePC76+l8Ct0byQz3wFaAl5/SLw5CT5uxWoD3l9FeDNdF3QTTfdrN+AC4APE2gMnh6MlT7g3CmO0zircVa3STYdg6ykmn4C/5kXAu8Avi8izcaYu6Y47uWw1/8i8F8+wJFAlzFm09hOY8yoiKwJ7puMN0KO8YpIG1Az9WkcwBjTIyJ3AU+LyLPAauARY8zWKQ5dDpwiIjdG2LckJG+bjTGjIfv2huXxRALBXFGUgxxjzN/C3npRRGYD/w94ZorDNc4eQOOsMg4dYqGkFGOM3xhTb4x5yxjzf8APge8moAqf+SLSuDWJ8n4o7gieuL8HxphPEgigzwBnAhtE5PopDsshcBvzuLBtCRD6Ry5SHmM5f0VRFAg0fBckcJzG2YnvKQcp2kBW0k0Ogdt0UxE+fc6pwObg7xuBKhE5YmyniOQRGCO3Mcn8uQk8QDIlxpgNxpifGGMuAO4GrgtxEMGzFjgy+A9D+BbP09DrgPOm/JSiKAcrx3PgIbbJ0DgbHY2zBznaQFZShoh8O/gU8CIROUxEPgl8mcCDDlNxkYjcKCJLROSzwEeAnwb3PQu8CjwgIu8UkaOCznzgV0lmexdwlojMivYgiogsFpEfBJ+wni8ipxIY+zd2K7KTwJi7d4tIrYiUB9+/BbhYRH4qIseJyCEicr6I3C0iBXHk8TbgAhG5Q0SOCZbtVSJyWGKnrCiKXRGRn4jI2cE4e5yI/AI4F7gjhsM1zkZH4+xBjjaQlVRSCvwfgd6GV4BPAV8FvjDZQUG+A5wDvAl8DfiqMeYvAMYYQ+DJ5C3AX4HXgFoCD6V0Jpnnmwnc0ttF9HlEBwncrnsQ2AY8BPwbuDGYPz9wA4EHZxqB14PvPwecDRxN4AGQtwj8MeoHPLFm0BjzD+A9BJ6qXkPgj9jH43EoijJtqCPQcN1MYBrNw4BzjDFPxHCsxtkoaJxVJPAdUJTsQQLzYl5pjLk/03lRFEWZjmicVZTJ0R5kRVEURVEURQlBG8iKoiiKoiiKEoIOsVAURVEURVGUELQHWVEURVEURVFCsM1Kes8//7zJy4tl+tyJeL1enM7Unar61a/+1IYSu59Dov6hoaHOFStWVKcgSxFJNM5ma/lliz8daahf/epPzB8tztqmgZyXl8fSpUsTOrahoYH58+dbnCP1q1/96fCnI41s9a9fv74hBdmJSqJxNlvLL1v86UhD/epXf2L+aHH2oBhiUVMT1xLw6le/+rPIn4407O7PNHYvP63D6le/+sPJigayiMwVkedEZLOIbBSRz1vp7+iINg+5+tWv/mz3pyMNu/szjd3LT+uw+tWv/nCyZYiFF7jZGLNeREqAdSLyjDFm01QHxoKIWKFRv/rVnwF/OtKwuz/T2L38tA6rX/3qDycrepCNMS3GmPXB3/sJLJk52yp/RUWFVSr1q1/9afanIw27+zON3ctP67D61a/+cLKlB3k/IrIAOJ7A2uf7aW9v55prrsHpdOLz+bjkkku44YYbaG1tpaioCIfDQV9fH9XV1XR3d2OMobq6mra2NoaHh6murmZgYICamho6OjoQESoqKujo6KC0tBSfz8fg4CC1tbW0trbicrkoKyujs7OTsrIy3G43w8PD+/fn5uZSUlJCV1cXo6OjlJWVMTIysn9/fn4+BQUF9PT0UFlZSX9/P263e//+goICcnNz6e3tpaqqit7eXjwez/79oefk8XjIz88fd07FxcUAlpyTz+fD4XCMO6fy8nKGh4ctOacx/1TXKdFzGhwcJDc3d8rrlOg59fT0kJ+fP+V1SvScWltbKSwsTKjuxXJODQ0NFBcXJ1T3YjmnHTt2UF5ebtn3KdI5DQ4OMmfOHEu+T3aKEekm0ThbX1/P7NmzLYtJ4WW4c+dOFixYYFlMCq8XjY2NLF682LKYpHFW46zd4mx9fT01NTWWtVvsFCOikVULhYhIMbAa+J4x5uHQfS+//LJJdBaLnp4eysvLLcih+tWv/nT705FGtvrXr1+/bsWKFctSkKWIJBpns7X8ssWfjjTUr371J+aPFmezYogFgIi4gIeA34c3jpPF5/NZqVO/+tWfRn860rC7P9PYvfy0Dqtf/eoPJysayBIYXX03sNkY8xOr/YODg1Yr1a9+9afJn4407O7PNHYvP63D6le/+sPJljHI7wSuBN4WkTeC733NGPOUFfLa2lorNOpXv/oz4E9HGnb3Zxq7l5/W4eisqu9m5dpmOgY8VBf3c/WyWaxYbP0DUZp/9WebPyt6kI0xLxljxBhzjDHmuOBmSeMYoLW11SqV+tWv/jT705GG3f2Zxu7lp3U4Mqvqu7njxT20D3gwQPuAhzte3MOq+m7L09L8qz/b/FnRQE41LpdL/epXv0396UjD7v5MY/fy0zocmZVrmxn1jX+Qf9RnWLm22fK0NP/qzzb/QdFALisrU7/61W9TfzrSsLs/09i9/LQOR6Z9wBPx/Y4o7ydDKvIfLZ92yb/6M+vPmgayiNwjIu0issFqd2dnp9VK9WeBf1V9N1c8uIEP/GEbVzy4ISW3zUDzn2l/OtKwuz/T2L38tA6Pp3fEy389uyvq/upi63vrUlE+VUWR81mS58DqKW7tdH3VHxtZ00AG7gXOT4XYjv+5qH9yQseWQWrHlmn+M+tPRxp292cau5ef1uEAxhhW7+zh2r9s5vmd+3AIOHPGL+Gb6xCuXjYr6bTCSUX5HD6zMOL7faM+blu1i+4h63qS7XB91R8fMTeQReRDUd7/oBUZMca8AKSkC83tdqdCq/4M+tM5tkzzn1l/OtKwuz+UVMfqSNi9/LQOQ8+Qh9tW7eJ7z+6md8TLsXXF3P2hI7j5jHnMDOmJPaauOCWzQFhdPj1DHl5r6gdgRn5gwq6ZRS7OW1JBoSuHl3b38smHNvPP7d2W9CZn+/VVf/zEM83b3cCfI7x/J/AXa7ITnWSXmnY4HCldRtYYk9Klpt1ud0qXQB0eHk7pEqjDw8OWLoEabQxZ+4CHhoYGy5dAHR4etnQJ1Mny39zcbOkSqB0dHeOOt3oJ1La2Ntxud8qXms7Pz0/pUtPZGCMSJOFYnWicbW1txeFwRKzrVpThWHmlaqnp1tZWiouLU77UdDbG2fb2dta0unlwyyD9bj/5TuHyI2dwcjVU5sFiVz8/PLOCLW19/GDdMK/v7eetnc2UOTxZHWfv2zjAsMfPMVUuvn7m7HFLTZ9V4+P+LcNs6HTzw9UN/HNbO1cdVUYB7oM6zhpjUrbUdDbHiGhMudS0iCwK/voWcDQQer9lEXCfMcaS+y0isgB40hhzVPi+ZJaaHh0dJS8vL7nMqT+r/Fc8uCHiAyQzi13cf9mE6pMUqcj/f/xhAx2D9s1/Ov3pSCNb/fEsNW1FrE40zmZr+WWLPx1pJOLvHHTzs5caWdPYB8CJs0v4z9PnMbM4N6L/52taeXpbN6ctKOOWcxZN+EwyWFk+u7qH+fQjWwC489LDmTcjf4LfGMMz27v51St7GXT7KHTlcP0pczj/0AoCa5dlLv/qT68/maWm64HtQCGwI/h6bLsPuDWhHKURO86/p/7JuXrZLPIcE4PY+w6vtjytVOR/YXn+hPfyUjS2z47XN91p2N0fJGOx2u7ld7DVYWMMf9/axScf2sKaxj6Kch3cfMY8/uv8QyI2jsf8V504i3xnYHjCWy0DVmV9v98qfvPqXvwG3nt4FfNm5Ef0iwjvPrSSuy49nFPmlTLk8fPTF/fw1b/voK0//tv12XR91W8NUzaQjTE5xhgH8GLw99BtljHmTstzZTG5uZG/8Oq3r3/F4gpuOn0eM4NPU+c7A1X5uZ09uH1+S9OyOv/rmvp4NTg2rrwgODau2MVNp89Lydg+O17fdKdhdz9kNlbbvfwOpjrc1u/ma3/fwU9e3MOg28fJc0v5zaVLOe/Qykl7TnNzc6kscvHhY2YC8Os1TfgtnAnCqvJ5rbGPtU39FOU6uOKEuin9lUUuvn3uIr7yrvmU5DlYv7ef6x7ezBObOuI6v2y5vuq3jpjHIBtjzrQ89RBE5A/Au4AqEWkCvmWMudsKd0lJiRUa9WeZf8XiClYsrmBgYICc3AI+/cgWdnQNs/K1Zq4/ZY5l6ViZ/33DHn60ugGAq06s4/LjaxkYGNg/3iwV2PX6pjMNu/tDSXWsjoTdy+9gqMN+Y/jr5k7ueq2ZYY+fkjwHnzl1DmcfUh7TkIIx/6VHz+SpLV1s7xzm2foezllizT/1VpSPz2+489W9AFx+XA1l+QeaOJP5RYSzF1dw/KwS/vffTby0ex//++8mXti1j/88fR6zSqe+dZ/p66t+6/3xzGKxUEQeEJFNIrIndLMiI8aYjxpj6owxLmPMHKsaxwBdXV1WqdSfpf7CXAdfOWsBOQIPbehgbVOfpX4rMMbw3y/soXvYy9G1xXzk2BpL/dGwuz8dadjdH0qqY3Uk7F5+070ON/eN8uWn6vnffzcx7PFz2oIZ3HXp4axYHPt42zF/gcvB1csDPbP3rG1mxGvNHTsryudvW7to6BmhtiSXi48cP9wuFn95oYtbzlnIN1YsoCzfyZstA1z/8BYe2dCOzz95b7Ld66j6JxLPLBYPEBjXdjMwZHlOUkh5ebn6DwL/4TOL+PiJdaxc28KPVjfwf5cspbwg+Qntrcr/E5s7WdPYR3Gugy+/az6O4Pyi06X87ZyG3f1hpD1W2738pmsd9vkNj2/q4J61LYx6/ZTlO/nsO+dwxsL48xLqX7G4gkc2dFDfNcxDb7fzH8fXJpX3cH8iDLp9/HZdCwDXnjSLXMf4/r94/GcsLOfYuhJ++XITz+3o4Vev7GX1zn3cfMY85s6Y+PxIvP5EUH/6/fEsFHIk8DFjzN+MMatDN8tzZTHDw8PqP0j8Hz6mhmNqi+kZ9vKTF/ZYMr+lFfnf1T3MnWsCt/5uOn3uuAdhplP52zUNu/vDSHustnv5Tcc63LhvhJuf3M6vXtnLqNfPWYeUc9cHD0+ocRzuzxHh+pNnA/DHN9vosmDBjWTL58E32+gd8XJkTRGnL5iRtL8s38lXz1rAt89dREWhk03tg3zqkS386c22iL3Jdq+j6p9IPA3kF4DjLc9BEBE5X0S2iki9iHzFSvfIyIiVOvVnsd+RI3wp+LDFmsY+Ht+U/PKTyebf7fXz/ed24/YZzju0YsIfqOlU/nZNw+7+MFIaqyNh9/KbTnXY5zf86a02Pv3IFja1D1JR6OTWcxfy1bMWjBuTm6h/jGNnlfCO+WWMeP3cF+y5TYZkyqet383DG9oBuO7k2RGHjSTqP3V+Gb+59HDOO7QCj89w12vN3PTENnZ1j2+Q2b2Oqn8ik35bROQ7IS93A0+LyMPAuPk0jDG3JJMJEXEAvwDOBZqA10TkcWPMpmS8Y9TWJn/7R/328c8szuXzp83lu6t2c+erezmmrpiFFQkvupB0/u96rZndPSPMLs3jM6dOfHhwupW/HdOwuz9dsToadi8/O9fhVfXdrFzbTMeAh/KCPnKdQmt/oEf33UsquP6U2ZTkJd4wHiNS/q89aRZr9vTy961dvO+IKg6pjLy0c6L+WLlnbTMen+GsQ8o5fGaR5f6SPCc3nzGfMxaW89OX9rC1Y4gbHt3Kfxxfy8xiF79d10LHgIfq4n6uXjYrJbMR2f07YEf/VD3Ic0O2IuAJwBX2/lwL8nESUG+M2WmMcQMPAhdb4AXsOf+e+pPzn7GwnAsOq8TjM3z/ud2MJvEgSTL5f7Wxl0c3duAQ+OpZCyhwOSz1x4Ld/elIw+5+0herI2L38rNrHV5V380dL+6hfcCDAbqHvbT2eyjOzeF75x3CF8+cb0njGCLnf05ZPu87ohoD3LmmOakhbYmWz+b2QZ7b0YPLIXxiknnkrSj/5XNL+c2lh3Ph0kq8fsNv17Xw49UHyr99wMMdL+5hVX130mmFY/fvgB39k35zjDFXW55iZGYDjSGvm4CTQz+QzFLTXq+Xrq6ulC0ja4yhra0tZUtN5+Tk0NTUlLIlUJ1OJw0NDSlbAnXMb+VS06HnJCI0NDRMuE4XzvbzepOD3T0j/GTVFm44dU5C5+T3+2loaIh7ac3te1r44ZrAbBrvX1zAvJIcGhsbJ5yTxzN+eexklwsNP6fR0dFxfquXQB0aGqKpqSmlS6C63W76+vpSttR0tsaIWLEqVicaZwcHB/c/RZ6KMhwcHKSvry9lS00PDg4yNDSU0qWmUxFn73m1j1HfxEZpnkOY6e9haMiZ8jh7ZrWHf7iE15v7+eu6es44tCZtcbaqqor/Wb0TgPccUsJIdwsjrsjnZGWcff9cw9EzKvnhmi7ChyOP+gy/eaWRo0u9toqzg4ODtLW1pWyp6WyOEdGYcqnp/R88sIxpOKNAizEm4S46EfkQcJ4x5trg6yuBk4wxnx37TDJLTff19VFaWppo9tRvY/+2ziFuenwbXr/htncv4uR5ZZb6o+E3hm88vYO1Tf0cN6uY2y9YTE6U6ZSmc/nbJY1s9cez1PQYycTqRONstpZftvhTlcZ5d71OpL/gAjx9rbXD0CfL/8Mb2vm/V/YytyyPX196OM6c+JdqTqR8XtjZw3ef3c2MfCcrP3wERbkT79Al45+KbCl/9SfnT2ap6THGljHdHvb7HmBURB4SkZqEchfoMQ69/TcHaE7QNYGenh6rVOq3mf/QqkKuXhaYs/PHL+xJ6GnrRPL/6MYO1jb1U5Ln4Etnzo/aOE7UHw9296cjDbv7w0hlrI6I3cvPrnW4ujjyNJbR3k+GyfL/3sOrmFWaR2PvKE9tSezB6HjLx+31c9drgWbCx5fVTdo4TsQfC9lS/upPjT+eBvIngd8DhwL5wGHA/cBngKMJDNf4RYL5eA1YEpzgPhe4DHg8QdcEKisrrVKp34b+S4+eyfGzSugd8fLj1Q1xL48ab/53dA1x96uBwP2F0+dRVTT5EpiZLp9s96cjDbv7w0hlrI6I3cvPrnX46mWzyHOM/+c7zyFcPclY3ESZLP8uRw7XnhRI83frWxl0+yz1R+LRTR209ruZX57P+YdOfex0Ln/1p8YfTwP528B1xpgdxhi3MaYe+DTwTWPMFuAqAktFx40xxgvcCDwNbAb+ZIzZmIgrEv39/Vap1G9Df44IXzpzPqV5Dtbt7efhDR2W+kMZ8fr5/nMNePyGi5ZW8c4I83Em408Eu/vTkYbd/WGkLFZHw+7lZ9c6vGJxBTedPo+ZwR7LmcUubjp9XkpmUZgq/++cX8bRtcX0jnj5wxvxPzAVT/nsG/bwwOuBNK4/efb+RZes8sdKNpW/+q33x/N4aw6wANgS8t48YOy+xkCcvnEYY54Cnkr0+Mlwu92p0KrfRv7KIhc3nzGfbz2zk3tea+a4umIWV8U2JVE8+b/zlb3s2TfCvBn5XHfKbMv9iWB3fzrSsLs/jJTG6kjYvfzsXIdXLK5gxeIKGhoamD9/fkrSgKnzL8HFQ258bCuPbOjgwsOrqCvJs8wfyv2vtzLk8bNsTgnL5sQ27nS6l7/6rffH04N8B/CsiHxPRD4lIt8FVgXfB7gQeNnqDFqBHeffU7/1/lPnl3HR4VV4/YGp34Y9sd0GjNX/74Z9PLmlE1eO8NWz5pPvjO3rlS3lk63+dKRhd38YaY/Vdi8/rcPW+A+tLuScxeV4/IZ7Xo3vMaJY87+nZ4QnN3eSI4FFQaz2J4r6p58/5gayMeaHwCeAWgJzFM8CrjHG/CC4/1FjzAWW59AC7Dj/nvpT47/u5NnMn5FPY+8ovw4u/WyFv3PQzX+/sAeATyyfFdeE+dlUPtnoT0cadveHkolYbffy0zpsnf+q4Ljc1bv2sbFtwHL/b17di9/ABYdVsqA89qkQs6V81G8ffzw9yBhj/m6MucYYc4Ex5hPGmL8nmwER+ZCIbBQRv4jENZ1RrMQzn6j6p7c/35nDV89agMshPLWli5d270va7zeGH61uoH/Ux7I5JXzgqOqY8xOLP1ns7k9HGnb3h5OKWD0Zdi8/rcPW+WcW53Lp0TMB+PUre2NePCQW//q9faxp7KPQlcPHTqiLyRuPPxnUP/38Uy01/XVjzPeCv38n2ueSXL50A3AJ8OskHJOSmzv5LALqP7j8iyoLuHb5LH71yl5++uIeDqsupHqSmSam8v/l7XZebx6gLN/JF8+YfEq3RPzJYnd/OtKwuz9NsToqdi8/rcPW+j9ybA1/39rFlo4hnt+5j7MOKU/a7/Mb7gze9bvsuBrKC+ObSi2bykf99vBP1YM8J+T38CVLLVm+1Biz2RizNRnHVPT29qZSr34b+t9/ZDXL55TSP+rjh8834AtfDilG/7aOIVYG5+L84hnzqIgzaE/ltwK7+9ORht39pCFWT4bdy0/rsLX+ApeDjwenOrvntWbc3qnXEZvK/4/t3ezsHmFmsYtLjpwZc15i9SeL+qeff6qlpj8d8nu6lp22nKqqKvWrfxwiwhfPmMf1D2/hzZYB/vJ2Ox85NvLaCdH8wx4f339uNz4DFx9RndAqfZP5rcLu/nSkYXd/pmO13ctP67D1/ncvqeCxje3s7B7hkY0dUeNrLP5hj4/frg10RFyzfBa5MT4AHavfCtQ//fxxTfUjIocDHwRqjDE3ishhQJ4x5q0pjvsngQdGwvm6MeaxWNJub2/nmmuuwel04vP5uOSSS7jhhhtiWtPc7XZTUVGRkvW/u7q68Pl8FBUVxb3+d6xrmhtj6OnpmXSd9mTOSUTo7OycdJ32ZM5pzB/v2vOxnpPb7U5o7fm+jlauPaaUH7/azb1rmzl0hoPaXM+Ec+rv78fpdE64Tj9atYO9faPMLXVxXq2XoaGhhM6po6ODvLy8hOpeLNepqamJgoKChOpeLNdp9+7dlJWVWfZ9inROIyMj1NXVWfJ9slOMSJREY3WicXbXrl3U1dVZFpPCy3D37t3MmzfPspgUXi/27t3LokWLLItJGmcD53T1CTP55j/38PvXW3hHnQv/cH/Uc4oWZ1tbW3lqj5fuYS8LyxycVJtHU1OTxlmL4+yuXbuorq62rN1ipxgRNY7GOoBeRD4E/BJ4CLjcGFMafKjudmPMOTFJJvc/D3zRGLM20v6XX37ZLF26NCF3qucnVL+9/b96uYlHNnYwqzSPX33gMApc45csjeR/YWcP3312Ny6H8POLD2NhReKNmWwvn0z705FGtvrXr1+/bsWKFXE9vJxMrE40zmZr+WWLPx1pZKv/m0/vYE1jHxctreJzp0Uf5RPN3z7g5po/b2LUZ/jpe5dwZE1x3HmYzG8V6revP1qcjec+xXeAc40xnwLGJpB9Ezg2oRylETvOv6f+9PmvWT6LRRX5NPeN8suXm6b0tw+4ueOlRiCwilMyjeNIfquxuz8dadjdH0baY7Xdy0/rcOr8nzxpNjkCT23tZHfPcNz+e9c2M+oznLFwRsKN48n8VqH+6eePp4E8k0CQBTAhP2Prgo6CiHxARJqAU4G/isjTyfgiYcf599SfPn9ucOq3XIfw9LZunt/RE9Xv8xt+8HwDA24fJ88t5b2HJz/uKdvLJ9P+dKRhd38YKYnVk2H38tM6nDr/vPJ8Llxahd/Ab9ZEXzwkkn9bxxD/rO/BlSNcs3xWQulP5rcS9U8/fzwN5HXAlWHvXQa8mkwGjDGPGGPmGGPyjDE1xpjzkvFFoqioyGql+qeZf355AdcHV2X62b8aaes/sGxlqP+Pb7bxdusAFQVObj5jHhLnlG6RsEP5ZNKfjjTs7g8jJbF6MuxeflqHU+u/8oRaCl05vNbUx9qmvpj8xpj9izm9/8hq6kpjX7Y6Fr/VqH/6+eNpIH8W+K6IrAaKgj29twH/aXmuLMbhcEz9IfUf9P6LDq/i1HllDLp9/OD53funfhvzb24f5L71LY6oZgYAACAASURBVAB88cz5zCiIf0q3SNilfDLlT0cadveHkfZYbffy0zqcWv+MAheXHx+4BX7nmr0Rp9UM9/+roZe3WwPzy3/0uMlnwIiFbC4f9WenP54Gch6wFPgF8A1gJXC0MWa75bmymL6+yP+xql/9oYgIXzhjHhWFTja0DfKHN9v2+wfdPm5/bjd+A5ceVc2yOaWWpDnmTyV296cjDbv7w0h7rLZ7+WkdTr3//UdUU1uSy+6eEZ7e1jWp3+Pzc9ergeEYV55QS3FeXBNuRSTby0f92eePp4H8JNAIXA54gW3AoOU5SgHV1fEt/av+g9dflu/kS2cGnoS9f30Lm9oGqa6u5hf/bqSl380hlQVcneRYuHDsVD6Z8KcjDbv7w0h7rLZ7+WkdTr0/15mzfxzxvWtbGHL7xu0P9T+xuZPmvlHmluXxnqXWzG+b7eWj/uzzx9xANsbMA5YDjwLHAH8GekTkyWQyICI/EpEtIvKWiDwiIjOS8UWiu7vbaqX6p7H/hNmlfOjomfgNfOHJbbz/ga38s74HhxB8mC/+Seonw27lk25/OtKwuz+UVMXqybB7+WkdTo//jIUzOGJmEftGvPwxeIcu3N834uX3rwceuLru5Nk4c5J/ziPUnyrUP/38cf2lN8bsBP4NvAy8QmAKofjXfBzPM8BRxphjCPR0fDVJ3wRinetZ/eofY355HgKEDpUTEbZ3Dlmelh3LJ53+dKRhd3+E9FIRqydLL1XqaeFPRxp28IsI158SeBj6oQ3ttA8ceBh6zP/711vpH/Vx/KxiTppr3VA2O5SP+rPLH3MDWUQeFJFG4D5gEfB7YIEx5qRkMmCM+Ycxxht8+QowJxlfJOzYta/+zPrvW986YU4sr9+wcm30aYoSxY7lk05/OtKwuz+UVMXqybB7+WkdTp//8JlFvGvRDNw+wz2vHYin1dXVNPWO8PimDoRA77EVswSF+lOJ+qefP56R78sI9EK8GdzeMMb0W5yfTwB/jLQjmaWmh4eHqa6uTtkysqOjo5SVlaVsqWmPx0N+fn7KlkD1+Xw4HI6ULYE65k/VEqiDg4Pk5uZaslzo2Dl1DHgiVtD2AQ8NDQ0JLxca6ZxaW1spLCxM2RKoDQ0NFBcXp2wJ1O3bt1NeXp7SpaYHBweZM2dOypaaztYYkSAJx+pE4+y2bduYPXu2ZTEpvAx37tzJggULUrbUdGNjI4sXL07pUtMaZw+c00Xznby0G57d0cMplT6OqC2hp6eHu7d48Bk4Y04+zoEORos1zqYrzm7bto2ampqULTWdzTEiGjEvNQ0gIrXAmcAZwGlAAfCCMebaKY77JxBpmZOvG2MeC37m6wQC+yUmQqaSWWq6q6uLysrKhI5V/8Hpv+LBDbRHaCTPLHZx/2VHWZqWHcsnnf50pJGt/kSWmobEY3WicTZbyy9b/OlIw27+u19r5o9vtnFkTRE/uWgJL21r5rYX28l35rDyw0dQWWjNNJpj2K181J8+f7Q4G9fcKcaYVhHZCswiMBTiLOCCGI47Z7L9IvJx4CJgRaTGsaKkm6uXzeKOF/cw6jtQHfMcwtXLrJ3BQlFSQaKxWlHSxWXH1vD3rV1sbBvkw/e/Te9oYFaLZXNKLG8cK0oixDMG+XER6QYeA44HngBONMbMTiYDInI+8GXgfcYY65+AItCdn0rUP/38KxZXcNPp85hZHAjUM4td3HT6PFYsrrA8LTuWTzr96UjD7v5QUhWrJ8Pu5ad1OP3+olwHpwQfwhtrHAO81tjHqnrrZySwW/moP/P+eHqQHwY+b4zZZXEefk5gYvtnggPyXzHGfMrKBGpqkl+FR/0Hn3/F4gpWLK5gZGSE/Pz8lKQB9i2fdPnTkYbd/WGkKlZHxe7lp3U4M/7XmycOjR/1BR6Gtrozwo7lo/7M+uOZB/neVARcY8xiY8xcY8xxwc3SxjFAR0eH1Ur1q1/9afKnIw27+0NJVayeDLuXn9bhzPg7BiM/DB3tIemk0rJh+ag/s35rVzzIUqycKkb96ld/ev3pSMPu/kxj9/LTOpwZf3Vx5LHG0d5PBjuWj/oz6z8oGsgVFdaPG1W/+tWfHn860rC7P9PYvfy0DmfGf/WyWeQ5xjdsUvUwtB3LR/2Z9R8UDWQ7du2rX/3qT18advdnGruXn9bhzPjT+TC0HctH/Zn1xzXNm10pLbVuuUr1q1/96fWnIw27+zON3ctP63Dm/GMPQ/f09FBeXp6SNMC+5aP+zPkPih5kn8839YfUr371Z6U/HWnY3Z9p7F5+WofVr371h3NQNJAHBwfVr37129SfjjTs7s80di8/rcPqV7/6wzkoGsi1tZFWuVa/+tVvB3860rC7P9PYvfy0Dqtf/eoP56BoILe2tqpf/eq3qT8dadjdn2nsXn5ah9WvfvWHc1A0kB999FH1q1/9NvWnIw27+zON3ctP67D61a/+cA6KBvLDDz+sfvWr36b+dKRhd3+msXv5aR1Wv/rVH85B0UD2er3qV7/6bepPRxp292cau5ef1mH1q1/94YgxxnJpKli1alUH0JDIsd3d3VUVFRWdFmdJ/epXfxr86Ugji/3zV6xYUW15hqKQaJzN4vLLCn860lC/+tWfsD9inLVNA1lRFEVRFEVR0sFBMcRCURRFURRFUWJFG8iKoiiKoiiKEoI2kBVlEkTkXhH5bqbzoSiKMh3RGKtkK9pAVhRFURRFUZQQtIGsKIqiKIqiKCFoA1lRQhCR40VkvYj0i8gfgfzg+18WkVdExBl8/WkR2Sgi+RnNsKIoio3QGKvYBW0gK0oQEckFHgV+B1QAfwYuDe7+EeAGviEiS4D/Aq4wxoxkIq+Koih2Q2OsYid0HmRFCSIiZwAPArNN8IshIv8GnjXGfENEFgDrgTbgPmPM9zOVV0VRFLuhMVaxE9qDrCgHmAXsNeP/a9y/qpgxZjfwHLAA+EVac6YoimJ/NMYqtkEbyIpygBZgtohIyHvzxn4RkfcApwKrCNwOVBRFUWJHY6xiG7SBrCgHeBnwAp8TEaeIXAKcBCAiVcDdwLXAx4H3BoO5oiiKEhsaYxXboGOQFSUEEVkG/AZYDDwVfHs7cATQboz5VPBzFxAI5kcbY7oykVdFURS7oTFWsQvaQFYURVEURVGUEHSIhaIoiqIoiqKEoA1kRVEURVEURQlBG8iKoiiKoiiKEoI2kBVFURRFURQlBG0gK4qiKIqiKEoI2kBWFEVRFEVRlBC0gawoiqIoiqIoIWgDWVEURVEURVFC0AayoiiKoiiKooSgDWRFURRFURRFCUEbyIqiKIqiKIoSgjaQFUVRFEVRFCUEbSAriqIoiqIoSgjaQFamJSKyQESMiJyW6bykGxG5VUTqQ15fJSLeTOZJUZTph8ZZjbPTGW0gKylDRHYHg2f4tjENyTcCdcCaeA4SEa+IXJVMwiIyJ3ie70rGkwQ/Bk7JUNqKoqQRESkSkdtFZKeIjIjI2yLywTQlr3FWmbY4M50BZVqzHHCEvC4C3gYeTHXCxhgf0JrqdLIRY8wAMJDpfCiKkhbuJNBQux7YCbwHeFBELjTGPJ3KhDXOapydzmgPspIyjDEdxpjWsQ04G3ABd092XLBX4PMi8pCIDIpIs4h8IewzdSLyoIjsE5FhEXleRJaF7B936y/k9YdF5AkRGQr2uFwZcsxuAg36lWO93ZPk8TQR+ZeI9Ae3N0XkvODuxuDP54Ke3SHHnRs8blhE9orIShGpDNl/r4j8U0SuE5EGEekTkcdEpDos/XNE5MXgefSKyGoROSS4b9ytP0VRpicikg98GPi6MeYZY8wOY8z/An8FvjbFsRpnNc4qk6ANZCWdXA88YYxpjuGz3wKeB44HfgD8UEQuARARAR4FlgIXAScBbcAzIlI1hfd24HfAMcCfCATpJcF9ywEfcBOB24Z1kQQi4gAeJ3Bb8YTgdiswFPzICcGflwYdy4PHnQ08RqAH/Rjg/cAC4JHgOY2xHDgLuBA4HziOwO28sfTPAZ4G1gGnAicD9xH450NRlIMHF4HG5kjY+8PAKSIyVUzQOKtxVomGMUY33VK+AcsAA5wXw2cN8Luw9x4AXgr+viL4mSNC9ucBLcAtwdcLgp85Lez1F0KOcRK4RXZ9yHte4Kop8lcedL0ryv45kfYT+EN0e9h784KfPS74+l6gA8gL+cxXgJaQ1y8CT06Sv1uB+pDXVwHeTNcB3XTTzfotGA9eC8a4HOACAg1kA9RNcpzGWY2zuk2yaQ+yki6uB3YB/4jx8y+Hvf4XcETw9yOBLmPMprGdxphRAj0NR07hfSPkGC+BHpGaGPM0dlwPcBfwtIj8TUS+IiKHxXDocuAmERkY24Cxc1gS8rnNwfMZY29YHk8k9nJUFGV6cwXQS2D8sZtAL+hdwX2+KY7VOHsAjbPKOLSBrKQcESkFPgrcaYL/aieiCXsdySNR3g/FHcET9/fAGPNJAgH0GeBMYIOIXD/FYTkEbmMeF7YtAf42RR5jOX9FUQ4yjDENxphzgGJgnjHmSAI9yH1AZ5w6jbMT31MOUrSBrKSDK4BcYGUcx4RPn3MqsDn4+0agSkTGejoQkTwCY+SSnULOzfiZN6JijNlgjPmJMeYCAg8eXhfiIIJnLXCkMaY+whbP09DrgPOm/JSiKAcNxpghY0yziOQCHwQeNcb4pzhM42x0NM4e5GgDWUkH1xMI1m1xHHORiNwoIktE5LPAR4CfBvc9C7wKPCAi7xSRowg8PJEP/CrJvO4CzhKRWdEeRBGRxSLyg+AT1vNF5FTgdA7cxuskMObu3SJSKyLlwfdvAS4WkZ+KyHEicoiInC8id4tIQRx5vA24QETuEJFjROQwCUxSH8vtR0VRphHBGRsuFJFFInImgd7WAqaYxSKIxtnoaJw9yNEGspJSROQUAk8S/zrOQ78DnAO8SSDQf9UY8xeA4DCN9wNbCExn9BpQC5xrjIn3lmI4NxO4pbeLwEMckRgkcLvuQWAb8BDwb+DGYP78wA0Epl9qBF4Pvv8cganujibwAMhbBP4Y9QOeWDNojPkHgblOTyYwHvBV4OPxOBRFmTaUAncQ6Pl9mMBY2lOMMXtjOFbjbBQ0ziqS+JBQRUkNEpgX80pjzP2ZzouiKMp0ROOsokyO9iAriqIoiqIoSgjaQFYURVEURVGUEHSIhaIoiqIoiqKEoD3IiqIoiqIoihKCM9MZiJXnn3/e5OXlJXSs1+vF6Uzdqapf/epPbSix+zkk6h8aGupcsWJFdQqyFJFE42y2ll+2+NORhvrVr/7E/NHirG0ayHl5eSxdujShY0dGRsjPz7c4R+pXv/rT4U9HGtnqX79+fUMKshOVRONstpZftvjTkYb61a/+xPzR4uxBMcSioyPaNIvqV7/6s92fjjTs7s80di8/rcPqV7/6w8mKBrKIzBWR50Rks4hsFJHPW+y3Uqd+9as/jf50pGF3f6axe/lpHVa/+tUfTrYMsfACNxtj1otICbBORJ4xxmya6sBYqKiosEKjfvWrPwP+dKRhd3+msXv5aR1Wv/rVH05W9CAbY1qMMeuDv/cTWDJztlV+O3btq1/96k9fGnb3Zxq7l5/WYfWrX/3hOG699VbLpckgIgsIrAn/9VtvvXV07P3169ffevnll3Pfffdx9913s2/fPk444QT27t2L1+tldHSU9vZ2cnNzaWtro7e3l7y8PPbu3YuI4PV66ejoID8/n5aWFvr7+3E6nTQ3N2OMYXBwkI6ODgoKCti7dy9DQ0Pk5OTQ0tICQF9fH52dnfv3Dw8PIyK0tLTgdDoZHh6mq6tr//7R0VH8fj+tra04nU66urro7u7ev9/tduPxeGhra8PlctHZ2Tluf+g55eXlsW/fvnHn5PP5GB4etuycOjs7x51TTk4O+/bts+ycOjs7p7xOiZ6T0+mku7t7yuuU6DkZY+jp6ZnyOiV6TqOjo/T29iZU92I5p/7+fvr7+xOqe7GcU3d3N8PDw5Z9nyKdk8/nw+FwWPJ9slOM6Ovra1m0aNGd6Ym+icfZrq4uRMTSmBRahj09PeTm5loak0LPqbe3l8LCQstiksZZjbN2i7OdnZ34fD7L2i12ihHR4mxWLRQiIsXAauB7xpiHQ/e9/PLLJtFZLDo7O6mqqrIgh+pXv/rT7U9HGtnqX79+/boVK1YsS0GWIpJonM3W8ssWfzrSSJV/VX03K9c20z7gYWaxi6uXzWLFYutvZ2v+1Z8pf7Q4mxVDLABExAU8BPw+vHGcLIODg1bq1J8l/lX13Vzx4AYuf7SRKx7cwKr67pSko/nPrD8dadjdn2nsXn5ahyOzqr6bO17cQ/uAB4D2AQ93vLgnJbFK86/+bPNnRQNZAo8f3g1sNsb8xGp/bW2t1Ur1Z9ifzsCn+c+sPx1p2N2faexeflqHI7NybTOjvvF3mUd9hpVrmy1PS/Ov/mzzZ0UDGXgncCVwtoi8EdzeY5W8tbXVKpX6s8SfzsCn+c+sPx1p2N2faexeflqHIzP2D3w4HVHeT4ZU5D9aPu2Sf/Vn1p8V07wZY14CUjZJnsvlSpVa/RnypzPwaf4z609HGnb3Zxq7l5/W4Ym82tiLAJGeUqoutv5cUlE+1UUu2gcnxtSqInvkX/2Z9WdLD3JKKSsrU/808g+MenE5Iv8/lYrAnYryiZbPkjyH5WnZ7fpmIg27+zON3ctP6/ABjDE8sqGdW/6xEwPkhIXaPIdw9bJZlqQVSirK550LIjuNMbT1uy1Nyy7XV/2xc1A0kDs7O9U/Tfz1nUPc8OhW3L6J/RrOnNQE7lSUz9Un1kW8ZdI36uNnL+3B7fVblpadrm+m0rC7P9PYvfy0Dgfw+g0/+1cjv3plL34DVxxfyxfPmEd1SI/rBYdVpmwWCCsxxvB2a+DBreLcQMdDZaGTigIHnUNePvvYVja2DViWnh2ur/rjIyuGWACIyD3ARUC7MeYoK912/M9F/eMxxvD3rV38/OUmPD7D4soCViwp55ENHfvHyRW5cjhzUXnSaYWTivLJdzn29874DcwscnHCnBJW1ffw1y1dbO0Y4pvnLKSuJC/ptOxwfTOdht39mcbu5ad1GPpGvNy2ahdvtgzgcghfPGM+Zx0SiKfnLKnkj+v2cPfrXbzdNogxxvKlfa0un9ea+qjvGqa8wMl9HzmS4YE+ZsyYwcCol+89u5t1e/v50l/r+c/T53HOkuQb/Nl+fdUfP9nUg3wvcH4qxG63tbdS1J9e/4jXz49f2MNPX2rE4zNcuLSSO957KJceVcP9lx3Fby+axZyyPHpHfTyzrcuiXB/A6vLxG8Pv1gcmPb/+5Nnc/77Z3P/Ro/jC6fP52XsPpa4kl/quYW54ZCsvN/QmnV62X99sSMPu/kxj9/I72Otw474RPvf4Nt5sGaCiwMmPL1yyv3E8xjtn5VFR4GRH1zCvNvYlm90JWFk+xhgeeL0NgA8ePZM8Z85+f3Gek++edwgXH1GFx2/44eoG7nmtGX+Sa0Jk8/VVf2LE3EAWkQ9Fef+DVmTEGPMCkJKJYIeHh1OhVX8a/E29I3z+sa08s72bPIfwpTPn8/nT5pHrPFB13aMjXHlCYIqX37/Rittn3fAEsL58Xtq1j53dI1QVurhwadU4/+KqQn7x/sM4dX4ZA24f33pmJ3e/uhefP/Hgnc3XN1vSsLs/lFTH6kjYvfwO5jq8rqmPzz2+jea+UQ6pLOB/Lj6Mw2cWTficzz3CB4+pAeCBN1qxepExK8vnzZYBNrUPUpLn4KLDqyb4HTnCDe+Yy43vmEOOwINvtnHbP3cx7PElnGa2Xl/1J048QyzuBv4c4f07gb9Yk53otLe3c8011+B0OvH5fFxyySXccMMNtLa2UlRUhMPhoK+vj+rqarq7uzHGUF1dTVtbG/n5+XR1dTEwMEBNTQ0dHR2ICBUVFXR0dFBaWorP52NwcJDa2lpaW1txuVyUlZXR2dlJWVkZbreb4eHh/ftzc3MpKSmhq6uLoqIi2traGBkZ2b8/Pz+fgoICenp6qKyspL+/H7fbvX9/QUEBubm59Pb2UlVVRW9vLx6PZ//+0HMqKyujqalp3DkVFxcDWHJOFRUVNDQ0jDun8vJyhoeHLTmnMf9U1yn8nP72dhMrNwww4jPUFObwxXfUMrvYT0NDw7hzKi0tZUFvL3NKXTT1efjDy9u5bPmCCdcp0XMqLCykoaFhyusUyzm1tbdzz6v7ADhvnovB/l5yc3MnnNP/O3UmD+S6ebh+mD++1c7rjd18Y8UiRns74z6nnJyccf546l4s5+T3+2lqarLs+xTpnJxOJ319fZZ8n+wUIxIk4VidaJz1+Xx0dXVZFpPCy9Dn89HX12dZTAqvFz6fj6GhoSnr+nSLs+v2ubhzbRt+A6fMKeaKQ124PIN0dvZNOKfS0lKO8fVSkpvD5vYhnn5jB2cdMTcr4+zKVxoAeO9h5bQ3N1FaWhoxzq6YX4ZrpIQ73x7kXw293PjwRr597iJ8/V0HXZz1+Xy0tbVZ1m6xU4yIxpRLTYvIouCvbwFHM346tkXAfcYYS56MEpEFwJORxiAns9R0Q0MD8+fPTy5z6k+b3+Pzc9drzTyyoQOAMxfO4D9Pn0dhbuQZHsb8L+zq4burdlNV6OLeDx8xrpc5nfmfjGfru7n9+QZmFrtY+aEjcDlyJvW/1TLAfz27i+5hLxUFTr529gKOqSuJK81su77ZmEa2+uNZatqKWJ1onM3W8ssWfzrSiMfv8xt+9UoTj28KPNj00WNr+PiyOnImGVc85n/g9VbuXdfCsXXF/OjCJZbkPdSfLJvaBrnpiW0UunK4/7IjKc5zTulv3DfCLf/Yyd6+USoKnNx67iKWRuhFT0f+1Z9+fzJLTdcD24FCYEfw9dh2H3BrQjlKI7m5ueq3ib9j0M3/+2s9j2zowJkjfObUOXzt7AVRG8eh/tMWzGBRRT6dQx7+usW6J1qtKh+f33D/64HJzP/juFpcjpwp/cfUFfPLDyzl2Lpiuoe9fOmpev74Zltc4+Wy6fpmaxp29wfJWKy2e/kdTHV4YNTL15/eweObOnHlBIatXb181qSN41D/xUdWU5Tr4M2WATa2WjcLhFXl84c3AjH24iOq9zeOp/LPnZHPz9536P44+8W/bue5HT1xpZst11f91jFlA9kYk2OMcQAvBn8P3WYZY+60PFcWU1ISX4+b+jPjX9fUx2ce2cqm9kGqi1z890VLeP+R1VM+LT3mzxHhYyfWAYExZSMWTZVmVfk8u6Obpt5R6kpyOffQypj9FYUubr9gMZcdW4PfwN2vNfPtZ3bRP+qNKd1sub7ZnIbd/ZDZWG338jtY6vDe3lE+9/g21u/tpyzfyQ8vXBzzDA5j/qJcBxcfERjX+8AbbYlnOIo/Geo7h1jT2EeeM4cPHFUdl78038n3L1jMe5ZW4vYZvv/cbu5b1xJzZ0Q2XF/1W0vM96CNMWdannoIIvIH4GXgMBFpEpFrrHKPjXtJFepPzu83hvvXt/C1v++gd8TLibNL+OUHlkZ8UGQq/6nzylhSVUDPsJcnN3Ukle9I/kTx+g33rw/2Hh9fizNk9v1Y/I4c4RPLZ/Gddy+iONfBy3t6ueHRrWzvHJry2ExfXzukYXd/KKmO1ZGwe/kdDHX4jeZ+Pvf4Vpp6R1lYns/PLz6MI2uKE/J/4KiZ5DtzeK2pj20xxKB4/Yky1mC/aGklMwrGL8YUi9+ZI3z+nXP59CmzyRG4//VWvv/s7pg6WzJ9fdVvvT+eWSwWisgDIrJJRPaEblZkxBjzUWNMnTHGZYyZY4y52wovQHm59XPjqt8af++Il288vYP7go3Hj51Qy3fPO4Sy/NifHw31iwgfD/Yi//Gt9qSeSo7kT5RntnfT0u9mTlnehEn24/GfMq+MX3zgMJZUFdDa7+amJ7bx1y2dkz5Rbvf6k4407O4PJdWxOhJ2L7/pXoef2tLJV/9WT/+oj5PnlvLT9x5KTUl8t6RD/WX5zv2zQ/whOGwsWZItn909w7y0ex8uh/DBo2sS9osIHzhqJre9+xAKXTms3rWPLz65na4IS1Yn4k8U9affH89TTA8AfuBm4MqwLaux4/QiB4N/c/sgn35kC2ub+inNc/C98w/hihPqcISvbRqnf/mcUpZWF9I74uUxC3qRky0fj8/PA8E/IlccXzvh/OL115Xk8dOLDuWipVV4fIafvdTIj1Y3RP1nwO71Jx1p2N0fRtpjtd3Lb7rW4bGH8e54qRGfgQ8dPZNbz1006TMdsfo/ePRMXA7hXw297O5J/tySLZ8Hg73H5x9aSWWRa8L+eP3L55bys/cF5qXf1jnEjY9tnbS33O51VP0TiaeBfCTwMWPM34wxq0M3y3NlMSMjI+rPIr8xhkc3dnDzk9vpHPRw+MxCfvmBpSybU2qJP7QX+c9vtTPoTq4XOdnyeXpbN20DbubPyI+40l8i/lxnDp87bS5fOnM+ec4c/lnfw+cf30bjvokuu9efdKRhd38YaY/Vdi+/6ViHB90+bvnHzv0PPH/h9Hl88uTZcXdARPNXFLp4z2GBZyn+YMFY5GTKZ2/vKM/v7MEh8OFjJvYeJ+qfXx6YF/qo2iK6hjzc/MQ2XtgV+eE9u9dR9U8knnmQXwCOB9ZZngtARM4HfgY4gLuMMbdb5a6trbVKNa38q+q7Wbm2mY4BD9XF/Vy9bNaE2/9WEJr/IbePn764h9W7AnMBf+DIaq49adb+GR2S9Y9xwuwSjqopYkPbII9s7OCK4xMvw2TK3+090Ht85QkTe4+T9Z+zpIJDKgu4bdUudveMcONjW/nC6fPGNcTtXn/Avt+xdPnDSGmsjoRdyy8ddTgTcbalb5Rb/rGThn0jlOY5uOWcRRxTF/t446n8Y3zomBr+uqWL1Tt7+NgJtcwuy7fUHyt/eiswl/N5h1ZEHTqSqL8s38kPLljM//yrkae3dfPdIj2CmQAAIABJREFUVbv5+ImjXH5czbgHyO1aRzNRP1NBKvyTtkpE5DtjG7AbeFpE7gx9P7gvKUTEAfwCuAA4AvioiByRrHeM1lZrxkhNJ/+q+m7ueHEP7QMeDNA+4OGOF/ewqt76xQzH8r+re5gbH9vK6l37KHDl8I2zF/DpU+ck1TgO9YcS2ov80NvtMc/4EKs/Vp7a2kXnkIdFFfmctnCG5X6AhRUF/Pziwzhz4QyGPX6+9+xufvlyE57gioJ2rz9gz+9YOv3pitXRsGP5paMOr9relfY4+3brAJ97fBsN+0aYPyOf/734sKQbx6H+UGYW53LO4gr8JjBzkNX+WGgfcPPM9m5yBC47NnLvcTJ+AJcjhy+cPo/rTpqFAL9d18LtzzcwGvLwnh3raCbaAakiFf6pepDnhr1+AnBFeD9ZTgLqjTE7AUTkQeBiYJMV8vz8xP+rna7+lWubGfWNf7BrNDiedUfXMKX5DsrynJTkOynNc1Ka7wj+dI6bgWEyxv4zbR/wUJrXy6Dbh8/AgvJ8vrliIXNnWHNe0crn2FklHFtXzJstAzy8oWN/g9kq/1SMeP08+MZY73H0SfituL6FuQ6+dvYCjtrUya/X7OXRjR28umcfoz7oGvIwM46eAWMMwx4/faNe+kZ99I146R/10jty4Pe/b+uOWH9Wrm1OSe+DHb9j6fSTvlgdETuWX7QY+JMX9vD4pg78JjCG128MPgN+vwm8ZwLv+f0c2GdM8LPB/cHfIz06OxZn9+wLLDdfWeSiqjCXyiIXM/KdcQ+BCI2zJXm9DIz6MMCyOSV8/eyFFCUw3jgS0a7BR46t4R/bu/jn9m6uOL4u7of/pvJPxZ/fasPrN5x1SPmkPdjJ1iER4YPH1DBnRj7ff243z+3ooaVvlHMWV/Cnt9toH4gvzo4x7PHRNeShczCwdQ956Ay+fmVPL17/xDr6g+cb+MkLe8gRyMkRckRwSGCq05wcgq8FR/D3/ftzJHBMcP/2ziE8EfypiON2jBGTNpCNMVdbnmJkZgONIa+bgJNDP5DMUtMOhyOly8jm5uamdKnpgoICy5ea7hiI/ETuiNfPX95un/RiFTiF0jwHhU4odApVJQU4fKOUF+ZSVpCLwzdCl9fFE1v34Qn+g903GhgHfFhFLp89roDKXL9l5+RwOGhoaIh4nS5akMubLfDQW62cMy8fp98d93Xy+XwJLYH6u1d30z3sZeEMF7PZx8hIfsRzcrvdE5ZATXS50OOKh/j6Oyr57zVdNPcfuMbtAx5+8sIe1u9sZWFFAcP+HDr7hvA68ugeGKHf7WPUONg35GbQa0h0Cun2AQ89PT2WLzVtjEnpUtPZGiNixapYnWic7e3txekM/DlJRRn29vZSUFBg6VLT7VFioMdv2NxuzdRl0Rjx+iOO280RmJGXQ0WhixKnn9qyQopzfJQ4DYvqqvAP7aNuRiFlhfn09vaycSCXX65pxR38vvYH4+yxM/O59jAn4h2lqd2a5bOjxVnvvlZOrsvn5eYR7v5XPTe+Y05C1ymROOssLt+/KNQFC/JoaGiIek5Wxdnja2v50onF/OLNQbZ0DLGl40BdGeuB7ezs5MyFZQz5Hezp6MWfV8ze7n66hjyMSB4t+wbp80DPiI8hT+yLPoWyv2HrS+z4yWgf8DAwMGBpnM3mGBGNKZea3v/BA8uYhjMKtBhjEl6VQUQ+BJxnjLk2+PpK4CRjzGfHPqNLTVvrv+LBDRH/QJTkOfjwMTX0jXgDPYgjvuDPQG9i/6gXfxLfx5lFLu7/6ISVxJNiqvL56t/qWbe3n48cW8M1y+NfFT2R8h/2+PjYHzfRO+Lltncv4uR5ZZb6p+LyP2ygc4ppiSYjz5lDaZ6D0nxn4Gfw7sHY69+/3rr/n55QZha7uP8ya68v2PM7ZoU/nqWmx0gmVh9MS01Hi4HlBU5uOWdhhF64kN64kF66sV68Az14gfccInz8jxtpj/A9LMlz8P4jq+kc9OzvPewa8tA7EttQsEJXDpWFLloH3HgiNJBS8T2c7Bo09AzzyYe24HII933kSCoLJ84ikYw/Gneu2ctf3m7nnfPL+Na50ap94v7J6Bn28B9/2DihhxcC/+gAMf2tdDmEykJXyN0EF5WFLiqLcvm/V5roGZ5YJ2YWubjnQ0cE72YcuNPh3383g+DriXc49n/Ob/j2ql3si+RPc/3JtD9anI3nIb16DtwxEsbfPfKLyOPAZ4wxiQxEamL8rcA5QHMCnohUVlZO/aGDzH/1slnc8eKecbcY8xyBpZ0nu7XiN4ZBt29/w7k/2IjuDTao+4PvvxB8CC+cjiQabdGYqnw+dmId6/b289jGDi45qprygviCdyLl//imTnpHvCytLuSkuZPPzpGK6zvZnJ3nLC7fP3SmLNjgLRlrCAffz3NOPi68NN8Zsf5cvSz+f0BiwY7fsXT6w0hlrI6IHcsvWgy87uTZcS2gMWkay+OLs26fn+4hD12DgdvsoY3nsZ9dg26GPH6GekejphvtDmEyTHYN5pcXcNqCGby0ex8Pvd3OdSfPttQfid4RL09uDvQefzSGh7CtrkPlBS58UVrAY2+X5TupKorU+A0Mq6kqclGS54i6WqzfmMhxdvkscqeI0bFw/cmz0xbH7Rgj4mkgfxI4E/g2geEQ84BvAv8GVgM/IPCg3QcTyMdrwBIRWQjsBS4DLk/AE5H+/v79t5VSgR39Y8F5bOzazGJXTGOnckQoyXNSkudkNnlRPxetd6a6OP6ehamYqnwOn1nEyXNLWdPYx5/fij94x1v+g24ff3or0Pb42Il1Uy6VnYrrW13silj+M4tdfOldC5L2J1p/EsWO37F0+sNIZayOiB3LLx11ON40ch051JbkUVsSPbYaYxhw++gc9PDlp+rZF6HXORNx9vLjanhp9z6e2NzJR46tiWuxp1j84TyyoZ0Rr5/lc0o5tKrQcn8sRIuzlYUu7vvIEUk/gJ7qOprOOG7HGBFPDf42sNgYMzbZXL2IfBrYZoz5tYhcBWxPJBPGGK+I3Ag8TWCat3uMMRsTcUXC7XZbpZpW/hWLK1ixuCIltz6i9c78f/bOOzyu6kzc7zdFvVjdvdsY4wYYCMlCIKZDICFASIEsIdkkS7JLNtmUJZtO+iZkN/llNwshmwChxZTQiQndFGNs427ZlmxZVpfVpdHMfL8/ZiSPxjPSlDvlWud9nnk8Mxq959zvnvl8dO6556TiL9NY4nP9qdN4/WA3j25v5UPLq+O6BBhv/B/e1krPkI+Tago5dcbE+8On4vymI/6pbD/h2PU7li5/GCnL1dGwa/zS0YatLkNCBik+8670jQBOdA4WVhaMDkSs3doSdx3iOce9Q14e3hbYBOqjJ0dfuSJRf6xEy7PJLl0aSqrbaLryuB1zRDxn0AHMDXtvNoEOLUAv8XW4x6CqT6jqYlVdoKq3JuqJhB3X37O7f83Ccm4+azbVRW6EwMjlzWfNztj6iosqC3j3nFI8PuW+OJcjiic+vUNe/hy8yfETMYwex+uPlWyLf7aXYXd/GCnN1ZGwe/zs2oaz7Xv+kVWBzzyyrZXeOJfWjCc+j25vo3/Yz8ppRTFPh5kM8Td+a4mng3wb8JyI3CoinxWR7wPrgu8DXAqst7qCVmDH9feOB/+aheXcde0yfnt+4N9UXX6Ptf7XnxJY5u3xnW209sX+12Y88Vm7tZVej4+V04pYNX3i0eN4/fGQbfHP5jLs7g8j7bna7vGzcxvOpu/50ppCVk0von/YzyPb2yz3Q+AG6LVbA4MQH10Ve6doMsTf+K0l5g6yqv4E+CQwlcAaxdOBG1X1x8GfP6yqF1teQwuIZ7kk4z9+/fMr8jl73hSGfRrX1qix+rsHvaOJ+/o41lzOlvhkqz8dZdjdH0omcrXd42fasHX+kU7r2q0tDAwfu9JNsv7Hd7bTPeTjxOoCVk2Pfc5ptsTH+O3jj2uSjKo+pao3qurFqvpJVX0q2QqIyNUisk1E/CIS13JGsZKTk9jC5cZ//PmvO2UqAjy1q53mnthGkWP1P/hOC/3Dfk6ZUczyqbEn7myKTzb601GG3f3hpCJXj4fd42fasHX+ldOKWFpdSM+Qb3SVCav8Hq+fB4M3QH901dSYprDF408G4z/+/BNtNX1LyPPvRnskWYetwJXAi0l6otLV1ZUqtfHbzD+nLJ9zFpTh9Sv3bIrtkkws/iMDw6M3jcS7Y182xScb/ekow+7+NOXqqNg9fqYNW+cXkdEb5x58p2XMdszJ+p/a3U7HgJcFFfkTLp+ZiD8ZjP/48080gjwz5PmscR4Jo6o7VHVXMo6JqKysTKXe+G3mv+6UqTgEnt7dTmN39LVE4/HfvyWw5NDps0o4sbowrvpkW3yyzZ+OMuzuJw25ejzsHj/Thq31nzazhIUV+XQOeHl6d7sl/mGff3T5zHhHj2PxJ4vxH3/+ibaa/lzI83RtOx2RZLaa9ng8lJeXp2wbWZ/PR2FhYcq2mlZVOjs7Ld1qOvSYRIS2trYJt/tN9JhG/BOdp0SPyePxxLVdqG9wkHPmlfLcvi5uf3U/N797+rjH1NPTg8vlinqeeoaVR7cF5h5fNNPBwYMH4zqm1tZWcnNzLduWOfw8NTQ0kJ+fb9m2zOHHVFdXR2lpqWXfp0jHNDg4yLRp01K21XS25ohYsSpXJ5pn9+/fz7Rp0yzLSeExrKurY/bs2ZZuNR16TIcOHWL+/PmW5SSTZwc4f6aT2na4Z2Mjp5b5KC4sSCrPvniwn5beYWaW5DDH1cPBgz0mz1qYZ/fv309VVZVl/RY75YhoxLzVNICInEhgcfkaVf28iJwA5Krqlgl+768EbhgJ5xZVfST4meeBL6vqhkgOs9W08VvpP9w9xA0PbAfg9qtOZGZpXsL+37zWwENbWzlzTinfmWC700T8yWJ3fzrKyFZ/IltNQ+K5ejJtNZ1OfzrKyDa/X5XPrN1JfecgXzxrNhefMP5OZ+P5fX7lxgd30Ng9xFfPmZPQShHZFh/jzx5/tDwb8016InI1gXnCM4Drg28XAz+f6HdV9TxVXRbh8Uis5SeDHdffM/7U+qeV5HLh4gr8Cn/cOP5c5PH8bX2e0RtRrj8lsePMxvhkkz8dZdjdH0oyuTpR7B4/04at9ztE+MjKwFzk+zY3Rd2WORb/C/s6aeweYlpxDufML4urHrH4rcD4jz9/PKtYfBc4X1U/C4ys3bIZWGl5rSzGjuvvGX/q/R9dNRWXQ3h+byf1nQMJ+e/d3MywTzlr3hQWVEy83Wm8fiuwuz8dZdjdH0bac7Xd42facGr8751fxvSSXBq7PTy/rzMhv1+PLst57coanI745h5P5LcK4z/+/PF0kKsJJFkADfk39jkaERCRD4pIA3Am8LiIPJ2MLxKFhfHdNGX8k8NfU5zDRSdUoIw/ihzN39Lr4cmd7QiBG/8SJVvjky3+dJRhd38YKcnV42H3+Jk2nBq/0yFcGxxF/tOmZvzjTOmM5n+1vov6I4NUFro5b1Him3BkY3yMP7v98XSQ3wKuC3vvWuCNZCqgqg+p6kxVzVXVGlW9MBlfJJxO58QfMv5J6f/IqhrcTuHF/UfY1x55FDma/55NTQz7lXMWlDG3LPFFyrM5PtngT0cZdveHkZJcPR52j59pw6nzr1lYRnWRmwNHBnmlLvpSXJH8qso9bwcGLz68oga3M66tGyb0W4nxH3/+eFrbF4Dvi8gLQGFwpPd7wBctr5XFdHd3G7/xR6SqMIdLlwSWh/nDxsMx+w93D/H0rnYcAh8/Obm5T9kcn2zwp6MMu/vDSHuutnv8TBtOnd/tdHDNipFR5CaiLQwQyf9mQze17QOU5bu4aIKb/CYiW+Nj/Nnrj6eDnAssAX4NfAO4E1iuqnssr5XFVFVVGb/xR+XalTXkOoVX67vY3dYfk/+eTU34FN63sJxZU6KvgBEL2R6fTPvTUYbd/WGkPVfbPX6mDafWf9HiCsrzXdS2D/BmQ+SOTLg/MHocmHt81fJqcl2Jjx5H8luN8R9//nha3GPAQeCjgBfYDfQlWwER+amI7BSRLSLykIhMSdYZTkdHh9VK4z+O/OUFbt6/NPDl+uNbx44ih/sPdQ3y7J4OHAIfW5X8nbPZHp9M+9NRht39YaQkV4+H3eNn2nBq/TkuB1cFR5HvfjvyKHK4f/PhXra39FGc6+SyE5PfBCKb42P82emPuYOsqrOB04CHgRXAA0CniDyWZB2eBZap6goCifzrSfqOIZ61no1/cvqvWVFNnsvB6we72dEyti8R7r/r7Sb8ChcsqmBGaW7SZdshPpn0p6MMu/vDykpVrh6vzFSpjwt/OsrIdv+lSyooyXWyo6WfTYd7J/Tfsykw9/jKZdXku5OfX5rt8TH+7PPHdc1CVfcBrwLrgdcILCFUnUwFVPUZVfUGX77G2C1TLcGOQ/vGn17/lHw3V5wUHEUOm4sc6j/QOchztZ24HMJHT65Jutxwfyqwuz8dZdjdH04qcvV42D1+pg2n3p/vdnLlskATHLnxLpp/W3Mvmxp7KXA7uGKpNVsIZ3t8jD/7/PFsFHKviBwE/gDMB+4G5qrq6RbW55PAkxb6AGhubrZaafzHof/q5dUUuB1saOhhW9PREY5Q/x/fPowSmFM3tTj50eNwfyqwuz8dZdjdH0qacvUY7B4/04bT47/ipCoKc5xsPtw7JseG+0fWPb5iaRVFua6kyw33pwLjP/788bS81QRGITYHH5tUtSeWX4xxq+lbCMyXuzuSo6WlhRtvvBGXy4XP5+PKK6/kpptuimlPc1Wlvb09Jft/t7e343A4aG5ujnv/71j3NHe73TQ0NIy7T3syx5Sbm0t9fX3ce8/Hekwj/nj3no/1mFwuF/X19XHvPR/pmM6fk88jtX3cueEQX1ge+JmIUF9fT5+rmBf3HcElcOXSMsuOye/3U19fn1Dbi+WYhoeHx/jjaXuxHNPg4CANDQ2WfZ8iHZPX66W7u9uS75OdckSCJJyrE82zAwMDtLe3W5aTwmM4MDBAd3e3ZTkpvF0MDAzQ399vWU4yeTb6MV28sJQHt3dwx/o6fnjJotFjGsmzXY4i3jjYTY4DLl1cavJsmvLswMAAzc3NlvVb7JQjoiHxzNsQkanAe4Gzgb8D8oEXVfVTMUsiez8BfBZYo6rHLiMArF+/XpcsWZKQv729nYqK5JaIMf7J4e8d8nLdfdvp8/j46SULWTm9eNT/3b/u4+W6Lq5YWsVN77ZuJpCd4pMJfzrKyFb/xo0b31qzZs3qeH8v0VydaJ7N1vhliz8dZdjF3zXo5bp7tzHo9fOrD5zA4sqCMf7v/nU/L9cd4UPLqvjMu0yeNf7U+6Pl2XjnIDcBu4BaoI7AqPDFCdUoiIhcBHwVuDxa5zhZenuPvSHA+I0/EkW5Lj60PDBP7v82HkZV6e3tpbatn5frushxCteusmbu8Qh2ik8m/Okow+7+cFKRq8fD7vEzbTh9/tI81+iqFPduOjoXube3l7rOAV6uO4LbKVy13ORZ48+sP545yI+KSAfwCHAy8BfgVFWdkWQdfgUUA8+KyCYR+e8kfcdQU2PtF834j2//B0+qojjXydamPjYe6qGmpmZ0E5H3n1hJRYHb0vLsFp90+9NRht39oaQwV0fF7vEzbTi9/quWV+N2Ci/XdVHXOTDqvzc49/iixRVUFJo8a/yZ9cczgryWQJKdo6rXq+rtViw8r6oLVXWWqq4KPj6brDOc1tZWq5XGfxz7C3OcXL0iMIr8h42HeX1PI68d6CbX5eCaldZ/Ce0Wn3T701GG3f1hpCRXj4fd42facHr95QVuLgnujDdyQ97Wuiae39eJUxjdec9K7BQf488OfzzrIP9eVfdbXoM0ICLGb/xxccXSKvJdwo6Wfm59vQuAk6cVUpZv7agG2DM+6fSnowy7+0PJRK62e/xMG06//+oVNbgcwgv7OjnUNciTdQP4Fc5bVE5NcY6lZYH94mP8mfcnt3ejTSgvLzd+44+LV+u7GPaPfW9jYy/raq3frceO8UmnPx1l2N2faeweP9OG0++vLsrhvIXl+BU+s3YnLzUMAjBnSp6l5Yxgt/gYf+b9k6KDbMehfePPrP/ODY14/WNXePH4lDs3NFpelh3jk05/Osqwuz/T2D1+pg1nxj97SmAteY/vaK79v7cOp2Qgwo7xMf7M+idFB7mkpMT4jT8uWnuH43o/GewYn3T601GG3f2Zxu7xM204M/6Htx/bqRlK0UCEHeNj/Jn1T4oOss/nM37jj4uqoshzjaO9nwx2jE86/ekow+7+TGP3+Jk2nBl/Ogci7Bgf48+sf1J0kPv6+ozf+OPihtXTyXWOnfSf6xRuWD3d8rLsGJ90+tNRht39mcbu8TNtODP+dA5E2DE+xp9Z/6ToIE+dGmmXa+M3/uisWVjOzWfNprrIjQDVRW5uPms2axZafyOAHeOTTn86yrC7P9PYPX6mDWfGn86BCDvGx/gz658UHeSmpqaJP2T8xh/GmoXl3HXtMn57fuDfVHSOwb7xSZc/HWXY3Z9p7B4/04Yz4w8diIDUDkTYMT7Gn1n/pOggP/zww8Zv/MZvU386yrC7P9PYPX6mDWfOPzIQcWnfiykdiLBrfIw/c/5J0UFeu3at8Ru/8dvUn44y7O7PNHaPn2nDxm/8xh/OpOgge71e4zd+47epPx1l2N2faeweP9OGjd/4jT8cUdWJP5UFrFu3rhWoT+R3Ozo6KsvLy9ssrpLxG7/xp8GfjjKy2D9nzZo1VZZXKAqJ5tksjl9W+NNRhvEbv/En7I+YZ23TQTYYDAaDwWAwGNLBpJhiYTAYDAaDwWAwxIrpIBsMBoPBYDAYDCGYDrLBYDAYDAaDwRCC6SAbDOMgIr8Xke9nuh4Gg8FwPGJyrCFbMR1kg8FgMBgMBoMhBNNBNhgMBoPBYDAYQjAdZIMhBBE5WUQ2ikiPiNwH5AXf/6qIvCYiruDrz4nINhHJy2iFDQaDwUaYHGuwC6aDbDAEEZEc4GHgj0A58ADwoeCPfwp4gG+IyCLgB8DHVXUwE3U1GAwGu2FyrMFOmI1CDIYgInI2cC8wQ4NfDBF5FXhOVb8hInOBjUAz8AdV/WGm6mowGAx2w+RYg50wI8gGw1GmA4d07F+No9vuqmod8DdgLvDrtNbMYDAY7I/JsQbbYDrIBsNRDgMzRERC3ps98kRELgHOBNYRuBxoMBgMhtgxOdZgG0wH2WA4ynrAC/yTiLhE5ErgdAARqQTuAD4FfAJ4fzCZGwwGgyE2TI412AYzB9lgCEFEVgP/CywEngi+vQdYCrSo6meDn7uYQDJfrqrtmairwWAw2A2TYw12wXSQDQaDwWAwGAyGEMwUC4PBYDAYDAaDIQTTQTYYDAaDwWAwGEIwHWSDwWAwGAwGgyEE00E2GAwGg8FgMBhCMB1kg8FgMBgMBoMhBNNBNhgMBoPBYDAYQjAdZIPBYDAYDAaDIQTTQTYYDAaDwWAwGEIwHWSDwWAwGAwGgyEE00E2GAwGg8FgMBhCMB1kg8FgMBgMBoMhBNNBNhgMBoPBYDAYQjAdZIPBYDAYDAaDIQTTQTZMKkREReTjma5HuhGRvxcRb8jrc4KxmJnJehkMhuMLk2NHX5sca3NMB9lgGSJytog8IiL1wcTwjQifOUlEHhCRPSLiF5Hb01zNacCD8fyCiNSKyLeTLVhEvCLy98l6EuQ+YEaGyjYYDBYQY469UETWi0ibiAyKyF4R+b6I5KSpmibHGo4LXJmugOG4ogjYDtwD3BblMwXAAeBR4F/SVK9RVLUp3WVmA6o6AAxkuh4GgyEpYsmx3cAvga1AD3Ay8FugEPhiqitocqzheMGMIBssQ1WfUNWvq+p9wFCUz7ypql9S1T8CXbG6RaRORG4VkdtFpDs4OvJjEXGEfKZYRP5HRFqDIycbROSCMM+Yy3/B1/8oIn8UkR4ROSgiXwn5+fPAAuBbwc+qiMyNUseTRORpETkiIn0iskNErhupP+AE7hzxhPzeqSLyjIj0Buu+VkTmhPz828ERlitEZGfQ/TcRWRBW/qki8lQwPr0i8oaInBH82ZjLfwaDwX7EmGPXq+q9qrpVVetV9WHgbuCc8dwmx5ocaxiL6SAb7MQXgEbgNAIjIZ8Hbg75+e+AC4GPExg1eQV4TESWTOD9FvAisAr4KfBjETk3+LMrgTrgPwhcOpwGHIzi+RPQDrwbWE5ghLwz+LPTAF+wviMeRGQp8AKwHlgNvC/4uWdFJC/EPQ34HPCxoH9K8HgJek4KHkNn0HEy8AvMd9xgmNQE89/FwN9i+LjJsSbHGkZQVfMwD8sfBBLeNyb4zPPA7XH4Xgp77wdAQ/D5QkCBS8I+sxH4XchrBT4e9vo/w35nJ/DDkNe1wLdjqGMX8Pfj/Nwb/nPg98C9Ye/lAv3AB4Kvvx383aqQz1wL+IG84Os/ApsBR5Sy/x7whrw+J3jsMzPdVszDPMwj/sdEORZoIDDKrMD/AM4YfCbHHv2MybGT/GH+8jHYifVhr18BZohICbA0+N6LYZ95EThpAu+msNeHgJoE6vcz4HYReT54ye6UGH7nNOCDwct1vSLSS2CEJA9YFPK5RlVtDaujANXB16cC61TVn0C9DQbD8cdZwCnAdcBlwDdj+B2TY8fW0eTYSYzpIBvsjMT4GZ3gM56w10oC3w1V/R6wGLgfWAa8JiLfn+DXHARGJlaFPRYDoSt8RKrjyO+Hv2cwGCY5qrpfVbep6l3Al4FviEhhnBqTY02OnbSYDrLBTrwr7PWZBP7q7wa2Bd87O+wzZ4X8LFE8BG7+mBBV3aeq/09VryIwYvO5CTwbgBXAXlWtDXt0EjtvAeeF3lBjMBgMQRzBh3uCz5kcGx2TYycZ5kQbLENEikRklYisAnKAqcHXC0M+kxPymSKgPPh6aTRvCKuCl9UWi8iAruIfAAAgAElEQVRHgX8mcJMEqroXeAD4fxJYB3SJiPySwCjDT5M8tP3Ae0RktohURkqQwWP/tYi8T0TmicjJwEUElmQK9ZwrItNFpDL43g+AE4G7ROT04O+eKyK/FJH5cdTxJwQuF94tIqtFZIGIXC0iZyZ0xAaDIeuIMcd+SUQuFZFFIrJQRK4lkB8eVdUjExRhcmx0TI6dZJgOssFKVgNvBx/TgJuCz0MvY00P+cypwAeDz5+Iwf9fwBwCIwK/An5DMHkH+RTwNHAXgZsp3gNcpqo7Ez6iAN8CSoFdQCswO8JnvEAZcAewI1iPZuCjIZ/5EoFj3h/0oKo7CNwxXRT8ne3A/wL5wET/mY2iqu8QuCmkisAd25sIXFb1xeowGAxZTyw51k2gw7op+LNvAL8GPhKD3+TYKJgcO/kQVTOlxpD9SGCNy9tVdaL5ZgaDwWCIE5NjDYaxmBFkg8FgMBgMBoMhBNNBNhgMBoPBYDAYQjBTLAwGg8FgMBgMhhDMCLLBYDAYDAaDwRCCK9MViJXnn39ec3NzE/pdr9eLy5W6QzV+4zf+1KYSux9Dov7+/v62NWvWVKWgShFJNM9ma/yyxZ+OMozf+I0/MX+0PGubDnJubi5LlixJ6HcHBwfJy8uzuEbGb/zGnw5/OsrIVv/GjRvrU1CdqCSaZ7M1ftniT0cZxm/8xp+YP1qenRRTLFpbWyf+kPEbv/FnpT8dZdjdn2nsHj/Tho3f+I0/nKzoIIvILBH5m4jsEJFtIvLPFvut1Bm/8Rt/Gv3pKMPu/kxj9/iZNmz8xm/84WTLFAsv8CVV3SgixcBbIvKsqm6f6Bdjoby83AqN8Ru/8WfAn44y7O7PNHaPn2nDxm/8xh9OVowgq+phVd0YfN5DYBvJGVb57Ti0b/zGb/zpK8Pu/kxj9/iZNmz8xm/84WTLCPIoIjIXOBl4PfT9lpYWbrzxRlwuFz6fjyuvvJKbbrqJpqYmCgsLcTqddHd3U1VVRUdHB6pKVVUVzc3NiAjt7e309vZSU1NDa2srIkJ5eTmtra2UlJTg8/no6+tj6tSpNDU14Xa7KS0tpa2tjdLSUjweDwMDA6M/z8nJobi4mPb2dlwuF83NzQwODo7+PC8vj/z8fDo7O6moqKCnpwePxzP68/z8fHJycujq6qKyspKuri6Gh4dHfx56TLm5uTQ0NIw5pqKiIgBLjik/P5/6+voxx1RWVsbAwIAlxzTin+g8JXpMOTk51NfXT3ieEj0mp9NJfX39hOcp0WMCqK+vT6jtxXJMXq93jD+ethfLMQ0NDdHQ0GDZ9ynSMfl8Prq7uy35PtkpR6SbRPPs4OAg7e3tluWk8BgODg7S3d1tWU4KbxeDg4P09/dblpOOpzz78Nv1PFw7QPugn/K8I3xsRTlnTM21TZ79y5YGHtk7SNuAj/K8I9xw6lSWFAyaPBt2TIODgzQ3N1vWb7FTjohGVm0UIiJFwAvAraq6NvRn69ev10RXsWhra6OystKCGhp/NvnX1XZw54ZGWnqHqS5yc8Pq6axZaP1lFlP/zPrTUUa2+jdu3PjWmjVrVqegShFJNM9ma/yyxZ+OMlLhX1fbwW0vHWDId7SfkOsUbj5rtuW5ytTf+DPlj5Zns2KKBYCIuIE/A3eHd46Tpa+vz0qd8WeBfyTxtfQOA9DSO8xtLx1gXW2H5WWZ+mfWn44y7O7PNHaPn2nDkblzQ+OYziXAkE+5c0Oj5WWZ+ht/tvmzooMsgdsP7wB2qOrPrfZPnTrVaqXxZ9ifzsRn6p9ZfzrKsLs/09g9fqYNR2bkD/hwWqO8nwypqH+0etql/safWX9WdJCB9wDXAe8TkU3BxyVWyZuamqxSGX+W+NOZuFNR/2j1jHZcyWDH85vuMuzuzzR2j59pw8fy0NaWqD+rKnJbWhakJj6VhZHrWV5gj/obf2b9WdFBVtWXVVVUdYWqrgo+nrDK73Zb/2Uw/sz4/arc8cahqD9PReJORXyi1dMhsLut39Ky7HR+M1WG3f2Zxu7xM234KKrK7zc08pvXAnnW5Ri7vmyuU7hh9XRLygolFfFZOa0o4vser48DRwYtLcsu59f4YycrOsipprS01PiPA/+wz89Pnq/nvi2BkY3wxO0QUpK4UxGfy088Ztt3RMCv8OXH9vDmwW7LyrLL+c1kGXb3Zxq7x8+04QA+v/Kfrxzknk3NOAS+fPZsvnT2bKpDRmLPmV+WkpuJrY6Pz6+80xSYl1qa50QIjChPL86hx+PnX/6ym12t1s1btcP5Nf74mBQd5La2NuO3ub/P4+MbT+/lub2d5Lkc3HrhgkDiDhmJ9SvMmWL9Xu+piM/BrsDoRb4r8BWsLnLzpbNm8b4FZQx6/fz7M3t5ene7JWXZ4fxmugy7+zON3eNn2jB4fH5++Lc6Ht/ZjtspfOu8+VywuII1C8u56yPL+IflhQC809yHz2/96ldWx+fF/Udo7vUwoySXez+6nN+eX849H1nGb65cwmkzS+ge8vGVJ2p5+1CPJeVl+/k1/vjJmg6yiPxORFpEZKvVbjv+5WL8R2nr8/Clx3bzdmMvZfkufnbZIk6bVRJI3Ncu4/6r5nHlssCI7H1bmq2o8hisjk9Lr4d1tZ04BH5z5RLuv2oed127jAsWV/KVc+bw4ZU1+BX+48UD3LXxMMkuxZjt5zcbyrC7P9PYPX6TvQ0PDPv45jP7eHH/EQrcDn540QLOnDPWd+6iKmqKcmjsHmJ9fVey1T0GK+OjqjwQ/L/gQ8urcTpk1J/vdvKdC+Zz7oIyBob9fOPpvby0/0jSZWbz+TX+xMiaDjLwe+CiVIg9Hk8qtMafBn9d5wD/9Ohu9nUMMrM0l9suX8ziyoJj/B9aXo3LIby0/wiHuoaSrfIxfitZu7UFr185a94UppfkjvE7RLjxtOl8/t0zcQj8YWMTt718MKkRm2w+v9lSht39mcbu8ZvMbbh70MtXnqhl46EepuS5+Nmli1gxrfiYz/m8w3xoeTUA929pTvoP93CsjM+mxl5q2wcozXNx/qLyY/wuh/DVc+ZwxdJKhv3Krc/t54mdSY7AZ+n5Nf7EibmDLCJXR3n/KisqoqovAtYvAgsMDAykQmv8KfZvOdzDF/+yh7a+YZZWF3Lb+xczrTg3or+qMIc1C8vwKzz4jrWjyFbGp3vQyxM7A1MnPryiJqr/8qVVfPO8eeQ4hSd3tfOtZ/cxMOxLqMxsPb/ZVIbd/aGkOldHwu7xm6xtuLXPw788toddrf3UFOXwi/cvYmHYAESo/8LF5RTnOtnZ2s/WZmvXnbUyPg8E/w/4wElV5AansYX7HSL845kzuf6UqfgVbnv5IPdtTrzjn43n1/iTI56tpu8AHojw/m+BB62pTnSS2Wo6Ly8vpdvIFhYWpnSr6dLS0pRuNV1eXp7SLVBH/PFsrbl7IJ9fvHIIr8K7ZhZx3WI3voEe2nqPPaaSkhLq6+u5cG4+z+yGp3e1c81J5Xi62y05poKCAsu2QP3j+noGvX5OnlaAu6+VTmfJ6FbZ4edpBl188+yp/OiVZt442M2XHtvN507Ko7I4L65jcjgcKd0C1e/3p3yraZfLldKtprM1RyRIwrk60Tzr8/lSuo3syFbjqdpq2ufzpXyr6WzLs7XNXfzn5j5a+7zMKHLy3fNm4j3STKc/8jGVlJTQ0tjABfOL+fOOI/zhtTq+f/Eiy7aatirP9ueUsqGhhxwnrJlbSH19PSUl0fPsWRUe3CdX8Lu327njzUba+wa5ZLqSm5s7qfKsz+dL6VbT2ZwjojHhVtMiMj/4dAuwHAhdOmA+8AdVtWTpABGZCzymqsvCf5bMVtP19fXMmTMnucoZf1r8qsqft7by29cDSwxdsbSSz75rJs6wFSui+b/z7D5eqe/iwyuqufH0GclXHuviMzDs47p7t9E95OOnlyxk5fTimPwNXYP821N7aerxMK04h1svWsDM0thvRsym85utZWSrP56tpq3I1Ynm2WyNX7b401FGPP7dbf3c8tReuga9nFhdwPcuWEBJ3vjjZSP+zoFhrrt3Gx6fcvuHTmR2mTU3RlsVn588X8dfazu5YmkVN717Zsz+v+3t5CfP1+FTuGBROV88a/a4/++Ek03n1/jjI5mtpmuBPUABsDf4euTxB+DbCdUojeTk5Bi/Dfw+v/Lfrx0a7Rx/6vTp/OOZ43eOw/3XrAxMW/jLjjb6PIlNSRjPnwxP7+6ge8jHkqoCVoSszzmRf2ZpHr98/2IWVeZzuMfDzY/uZnsclzez5fxmcxl29wfJWK62e/wmUxve1NjDVx7fQ9egl9Uzi/nRxQsn7ByH+svy3aPzeh+wcDqbFfFp6fXwt72BG6CvXD52Kc2J/OcuKOO7Fywg1+XgmT0dfPev+xny+mMuO1vOr/Fbx4QdZFV1qKoTeCn4PPQxXVV/a3mtLKa4+NgbDow/u/wer59bn6vjoW2tuBzC186ZwzUragjsQh67/8TqQlZOK6J/2M9fdrQmVe9I/kTx+nV0bvSHV449rlj8ZQVufnbpIk6fNbI80R5erY/tzutsOL/ZXobd/ZDZXG33+E2WNvxy3RFueWov/cN+zpk/he+cP598tzNu/1XLqxHgudpO2vut2f3Tivg8tLUFn8LZ86Ycc79KLP7TZpXw44sXUpzrZP2BLm55am/MAy3ZcH6N31pivklPVd9reekhiMifgPXACSLSICI3WuUemfeSKow/OX/3oJevPVnLy3WBJYZ+cNEC3hfHQvTh/g8HR5Ef2tqKJ44RgFj9ifD83k5aeoeZVZp7zPJJsfrz3U6+c/58LlpcgcenfPev+/nL9on/CMj0+bVDGXb3h5LqXB0Ju8dvMrThp3a18/11+xn2K5cvreRr587F7Yx9IatQ/4zSPN4zt5Rhv/LwNmsGIpKNT++Qlyd2BRxXB2+ATsS/tKaQn126iPICF1uaevny43vojOGPgEyfX+O33h/PKhbzROQeEdkuIgdCH1ZURFU/oqrTVNWtqjNV9Q4rvABlZWVWqYzfYn9zT+Au6q3NfVQWuPnF+xezanp8fwmG+0+dUczCinw6B7w8syf5hVGSjY9fdXR95g+vrMERNioej9/pEL541iyuC955/V+vNnDHm43j3nlt9/aTjjLs7g8l1bk6EnaP3/Hehu/f3MzPXzqAX+HjJ0/lpjNnHpOH4vWPdEIf29FGvwXT2ZKNz2M72xgY9rNqehGLIqzEEY9/Xnk+v3j/YqaX5LK3fYAvPraHpp7xlw+1exs1/mOJZx3kewA/8CXgurBHVmPH5UUmg7+2rZ9/fnQXB44MMrcsj19esZh55fHfuR/uF5HRUeQHtjQnvetTsvF5/UA39Z2DVBa6OXfBsV/ieP0iwnWnTONfzpqNQ+C+zc385IV6hn2RR8vt3n7SUYbd/WGkPVfbPX7HaxtWVf739UPc/mYjADedOZPrT50W09S1ifwnVheybGohfR7f6MhtMiQTH4/Pz8NbAyPZ10QYPU7EP604l19ctogFFfk0dg/xxb/sYX9HdIfd26jxH0s8HeSTgOtV9UlVfSH0YXmtLGZwcND4s8y/oaGbLz2+h44BLyunFfHzyxZRVZjYJPtI/r+bO4XpJTkc7vEkvUtSMvFRVe7bHBg9vmp5dcRLmon6Lzqhgu9dsIA8l4N1tZ184+nI8+Xs3n7SUYbd/WGkPVfbPX7HYxv2+ZWfv3SAB95pwSnwtXPmcMVJVVF+O34/wNXLR6azBTY/SoZk4vNcbScdA17ml+dx6ozIVyAT8Y/c+7F8ahHt/cN8+fE9UW+QtnsbNf5jiWcd5BeBk4G3LK8FICIXAb8EnMDtqvojq9xTp061SjWGdbUd3LmhkdbeYaqKerhh9XTWxDF39nj3jxAe/2d2t/OLlw7gUzhn/hS+/N455MQxF24iPwSmIly9ooZfvnyQ+7Y08975UxIaNYnmj5WtzX1sb+mjONfJxSdUWO4/bVYJP7tsEf/+9F7ebuzlS4/t5vsXLqAy5I8Nu7b/dJSRqe9Aiklprg7F7ufneGjDI4S2MY/Xzw/+Vser9V3kOoV/P28ep89KbiveSG34jNklzCrN5WDXEM/v7eS8RYkfV6LfEb8qD77TAsBVy6Pf2J2ovzDHyQ8uWsAPnqtj/YEuvvpkLd9cM4/TZpVY4p8Iu38H7Owft4MsIt8NeVkHPC0ia4Gm0M+p6jeTqYSIOIFfA+cDDcCbIvKoqm5PxjtCU1OT5evvravt4LaXDjDkC/zV3NI7zG0vBab4WXFyEvH7VfH6FZ9fGfYF/vWq4vUF3g99vHGgi/u3tDDsT039QxmJv6ryp03N/P6twwBcvbyaG0+fHvdcuGj+cM5fWM4f3zrM3vYB3jrUw+qZJRF+O3F/LIyMHl+xtCrq3eLJts/FlQXcdvlibnlqL/s6BvnnR3dz60ULmFuWb4k/EpHa589fOkBD1yDvml2KyyEhDwcuh+B0gNvpwBl83ymM+0dLNn7HEiUV5yCUdOXqUOx+fuL1+zWQV4d9foZ9isenDPtHnvuDPwt57vezubGHZ/Z0jo6upqON9Xl8fPvZfWw+3EtRjpPvXTifk2qKJhbE6A/FIcJVK2r4xUsHeGBLM2sWliU8EJHod+T1A90cOBKYwnZOhClsyfoBcl0OvnnePH7+0gGe3dPBN5/Zy1fOmcO5C46eQ6u/46rK4zvb+M1rhxgOaaP/8eIBdrb0cfKMYnKcDtwOwe10kOMU3E4JvOcMfc8RNddm23cs2/zjbhQiInfGIlHVG5KqhMiZwLdV9cLg668HvT8c+UwyG4U0NzdTUxN5XlKifPzerbT0Hntnq1MCd/hCcJX+YJuUo0852k5l9Hlo0xWB/R2DES9ZOQXKC9yBzm/YI8krXKP+ldOLmVqcE3gU5Y4+L81zJZT8mpubqayq5levHuTxne0I8LkzZ/KBJC73hfujnd/7Nzdz+5uNrJxWxE8vXWS5fzz2tQ/w2Yd2kutycNe1J1EaZa1Rq9pn96CXbz6zj+0tfRTlOPn2+fNYMa04Yb/Xr7T0emjqGaKpx8PhnqPPd7f2Y0FzC3acBXfw39COdVPPEL4IhbgcwtyQzQkUOJrGAk9UGa2fHn07+NnAi8buyP7qIjd3XXvMXkVJkeg5iHWjEKtydTx5drwcOK3k2C3hx6/Xse9FO/8jOXYkv47kVpFAPg3Ns8LR90Zej/xsd2v/6ABBuL+mOAePVxn2BzrEnuAgg1U4BU6ZUcLM0tzgI48ZpblUFroTHjBobm4mp6ScW57aS237AOUFLn540cKE7uuI5o/Uhj0+P9ffu42OAS8/uGhBwgMRiX5H/uWx3Wxt6uMfzpjBVcurLfeH4lfl9jcaefCdFgS46d0zuXxpVVL+gWEfjd1DNHQFHoe6BgP/dg/RM2TNWv4Co53l0E704Z4hIt264nYIi6sKRr+XftWj+VRB0dEcq8EEO/p65DOqUXOsU2Bq8bE5It6mf9iiHB4tz447gpxsxzcOZgAHQ143AGeEfiCZraadTqfl28hG+o8BwKdw4Ejq5tr4FFr7oi8543IEGp/LEfir0SFKjssF6sPtEHJcLtTvpa478hfPp7DxUE/En+U6hZoiN2VuZWZ5IaUuH2Vu5YRZ1UhfJ2XFBWO2oXxyexMP7uqlY9CP29HIsD/wxfvUskIumFdo2fbZTqeT+vr6iOfp1DIX97gdbD7cy4b9Lcwq0Li3ofT5fAltgfq79fsAOHdOAUeaD5Eb5Zg8Hk/ELVAT2S70pmU53Lndx4amQf718VqKch30DPkpz2vihlOnsqRgcHS70M4jR3AWTmFf8xGaejwMugqoa+2ic0hoHfDS3u9LqBM8f0oOXr/i8fkQh5OhYR8+Bb8S/GMOfKohr5Xx7w8fi9ev1Lan7oaP1t5hDh48aOlW0zk5OSndatqqXB1Pnh0vBzZ0xXNG4yMdObax2xPxZy5htIPhQMl1OXA5wIFSmJeD3ztMjstBQa4b//AwbzZH9vgU3mzo5s2Gse/nOKCmyMWMklzK3T4WVJcyxTlMRY6fBbOmR9zuNzTPOqQRv0JNoYt/XlVITZ6mJc+umZ3HA7t6+eMb9SwumZPQdr+J5Nl2itja1EeBW3hXlVBfXx/1mKzKs9evrME/0MPa2gF+9WoDr+xp5mCf0tbvjZhnu7q6KCuvYG9TBw1dQ/Q7C6ltPkLbIDT1eWkfSKwTfMrUfIZ9foaGfeB00T80jE/BhwRyLhK8khHItUM+ZcgXW1nDfmVbHJtRxYtP4VB36nJES+8wQ0NDqd9qevSDR7cxDWcIOKyqCS84KyJXAxeq6qeCr68DTlfVL4x8Jtu2mo42elKe7+LHlywMG9U6+lxDuhxH3zs62kXwr7NvPrOPzgHvMf6KAhe3vf+EwCibU8aMwDkmuGQda/2/eNZsmkJGC5t6PTT1eCZcML0k18nU4sCI85DXx1uHeo8ZdfnYqho+sdqSnclHmej83vFmI/dtbubdc0r59vnRmnHi/kgc7h7ihge2I8DvrzmJmuLoNyBa3T59fuWWp2rZ2Ng75n2nQ1heU4jLKTT1eGju9YxeuouEAJWF7tFzOq04Z/T5D56roy3C2qDx/OUeOiXI6w9OBdKjz/81eBNnOGX5Lr53wYLR0cORuoY3fQm5QkPIzyX4s688EdmfihHkdGw1PUIyudqKEeTyAhc/uWTs1ZpYslL4+fvXx/fQ3n/s+SnPd/GjSxYezZ96NK/69egVg/BRLg28OZpvv7duP0cinP/yfBf/cdmiY0bb3A6J+wraeHn2C++ZxaGuIQ52DXIoOHp4ZPDY+oxQkuscHWkeGXU+1DXI3W83jV5ihkCsbzpzJpdbdIVuhPHacO+Ql4/du42BYT+//sAJEZdZS8Yfje/+dR8v13Xx4RXV3Hj6DMv94/HEzjZue/ngMe+7HcL7FpZRnOuiITgaHG2kEwJXxKYV5zBzSh4zSwLndkZpHjNLc/mnR3bREmFALN4c5fOHTv05+vwrUXLslHwX/75m3pgrL2PyZ7CvIYzNvaNXbIK5+atP1Eb0VxS4jrmiO1FXNNKPv/LEHjoi5Ii0jCCHURtSR2Fsff0i8ijwj6qayN6TDcCskNczgcYEPBGpqIh8c1Qy3LB6+pi5LxAYYf30GTOYU5b8Ja1/OGNGRP+nTp8xbmcrVsar/xmzI9/M0TPkDXacQzrPwefNvR66h3x0D/Wzu60/arnP1nZY3kGe6Px+8KQq1m5t4dX6Lg50DjI75PK8Ff5IPPhOC36F8xaVT3i+rG6fTodEHMHz+ZVNh8d2mkvzXEen0wQ7v1OLAs+ri9xRNxK48fTI7eeGOM6tQ4QcpwRuy43Ap6N8B/7hjBksror/P+BY/fEcQ6ykIgeNQypz9ShRc8jpM5g9Jb7vWCQ+dXrk8/PpM2aMzq9Phs9EOf+fPmPG6DS5ZBkvz75n7pRjPt875OVQyOX2hpDOc/eQj+0tgZt+x0OB+99ptryDPF4bLsp1cfEJFazd2soDW5r5t/fNs9QfiUNdg7xS14XLIXzgpOhTKxL1T8QlSyq5c0MjXYNjB46G/crTu49df7+q0D36h02gExx4XlOUg9MR+Q+vG05LPs9C4P+EfIeTfPfY96PlwM+cMYPlU5Oftx7N/6nTZzDTgu/Yp6PkCKtyeDwd5E8D7wW+Q2A6xGzg34FXgReAHxO40e6qBOrxJrBIROYBh4BrgY8m4IlIT0/P6GUlqxiZAH7nhkZaeoepLnJbevdkNvqLc10U57oijg74Vekc8I52nH/8fH1ER2uUy7LJMNH5LS9wc+GiCh7b2cb9W5r58nvjG0WIt/10Dgzz9O7AuqDXrJg4caeifY43Dec7588f7RTHus1sOKlun+koIx3HMEIqzvE4pDJXj2L385ONbbgo18UJVS5OqCoc876q0tHvDYxIdg/RcCQwMvn6we6Inkzk2SuXVfPItlZe3H+ET/YMRZxjmow/nD+/04oCaxaWUVHonvDzqfgOdg9Gv6p6w+ppgU5wSR7TS3PJc8W/UpPdvwN298fTQf4OsFBVRyZ/1YrI54Ddqvo/IvL3wJ5EKqGqXhH5PPA0gfGk36nqtkRckfB4Is8DS5Y1C8tZs7A8JVM47OZ3iFBR4KaiwM1JNUcbbDhVRRMnsniJ5fxetaKaJ3a1sa62g+tPnUZ1Ueyj8PG2n4e3teLxKWfOLo1ppCsV7bOqyB0x/tVF7mO2uk6UVLfPdJSRjmOA1OWgKKQsV4dj9/NjlzYsIlQUuqkodLMyZKfRaFM4MpFnq4tyOGdBGetqO1m7tZV/PHOmpf5QOgeGeWZPcFvp5bHdGJfuPPuRVdYs+2b374Cd/fH8SeMA5oa9N5ujF0h7ia/DPQZVfUJVF6vqAlW9NVFPJFK9BqnxH8sNq6eT6xx72ShVl69jqf/0klzOnjcFn8Kft7ZY7h+hz+PjL9vbAEZ387PSHyvZFv9sL8Pu/jBSmqsjYff42bUNZ9v3fGQViSd3tdM9znzqRP0jPLq9DY9POWNWScxT5iZD/I3fWuLpIN8GPCcit4rIZ0Xk+8C64PsAlwLrra6gFTQ1NU38IeO3lDULy7n5rNlUB0cyqovc3HzW7JRcvo61/iMd1id2xpe844nPEzvb6PX4WD61iKU1hRP/Qpz+WMnG+GdzGXb3h5H2XG33+Nm1DWfb93xBRQGrZxYz5PXzlx1tlvshsCzao9uD20rHOAgRjz8esi3+xm8tMY8iqOpPRGQLcDVwCnAYuFFVnwr+/GHgYctraAHxLJdk/NYxcumjpaWF6uqJ5+ImSqz1H0neGxp6eGR7K9edMs1Sv8fnHx2d/vDK2I93ssQ/m8uwuz+UTORqu8fPzm04277nVy+vCeTYba1ctbya3Bjn3sbqf2Z3Bz1DPpZUFbAsxkGIePzxkm3xN37riGvWuKo+pao3qurFqvrJkYSbDPyDGgMAABYmSURBVCJytYhsExG/iMS1nFGs5OQkv+qD8R8f/muDIw4Pb2tlYDi2NSFj9a/b00FHv5f55XmcFsdi+dkUn2z0p6MMu/vDSUWuHg+7x8+0Yev8q6YXsbAinyODXp7dc+xqDsn4ff6j20pfvSL6ttKJ+pPB+I8//7gdZBG5JeT5d6M9kqzDVuBK4MUkPVHp6upKldr4beZfPrWIE6sL6Bny8dSudsv8Pr9y/5aR0eP4Enc2xScb/ekow+7+NOXqqNg9fqYNW+cXEa4Ort7z53da8MW4A2Es/pf2H6G518P0klzeHefNxtkSH+O3j3+iEeTQ21BnjfNIGFXdoaq7knFMRGVlZSr1xm8jv4iMzkV+8J2WmLaPjcX/Sv0RDnUPMbU4h7PnlcVcn1j9yWB3fzrKsLufNOTq8bB7/EwbttZ/9rwyaopyONQ9xPoDsXVcJvKrKg+8E1i6+6rl1VHXDk7UnyzGf/z5J9pq+nMhz9O17XREktlq2uPxUF5ebulW06HbUPp8PgoLC+PeWjN0u9DxttZUVTo7Oy3ZLjTSMYkIbW1t425hnMwxjfgnOk+JHpPH44lru9BpOsCs0hwOdnl4ZGMd5y+uGPeYenp6cLlcUc+Tw+Hgj28E9ow9f1YOjYca4jqm1tZWcnNzLdlqOtJ5amhoID8/P6G2F8t5qquro7S01LLvU6RjGhwcZNq0aZZ8n+yUI2LFqlydaJ7dv38/06ZNsywnhcewrq6O2bNnW5aTwtvFoUOHmD9/vmU5yeTZAc6d4eLeXR7uevMAC/MGKSgoSCrP7mj3sKdtgNJcJ0sLBjh48KDJsxbm2f3791NVVWVZv8VOOSIaMW81DSAiJxJYXL5GVT8vIicAuaq6ZYLf+ysQaQ2OW1T1keBnnge+rKobIjmybatp47e3/5nd7fzsxQPMmZLH/3xoCY5xpkRM5N94qJuvPbmXKXku/njtSTHflBKrP1ns7k9HGdnqT2SraUg8VyeaZ7M1ftniT0cZ2eYfGPbx8Xu30TPk4+eXLWLZBDuzTeT/t6dq2dDQw/WnTOXjMd5gHY8/WYzfvv5oeTbm/8lF5GoC84RnANcH3y4Gfj7R76rqeaq6LMLjkVjLTwY7rr9n/Kn1n7ugjMpCN/VHBnn9QOTdqGL137c5cNnvg8uq4u4cx+JPFrv701GG3f2hJJOrE8Xu8TNt2Hp/vtvJ+08MXPZ+YMvEa8+P59/fMcCGhh5yXQ4uX5rYFtrZFh/jz35/PP+bfxc4X1U/C4zc/r8ZWGl5rSzGjuvvGX9q/W6nY3RR+3s3NzHelZTx/Lta+3i7sZcCt2P0P4N4ycb4ZJM/HWXY3R9G2nO13eNn2nBq/FcsrcLtFNYf6OLAkcFxPzue/4HgyhUXLS6nJC+xPW6yMT7Gn93+eDrI1QSSLICG/Bv7HI0IiMgHRaQBOBN4XESeTsYXicLC2NdKNP7J47/4hAqKc53saOnnnaa+hPz3bQ4k7kuXVFKUm1jiztb4ZIs/HWXY3R9GSnL1eNg9fqYNp8ZfVuDm/EWBTTMenGAUOZq/pdfD32o7cAhcuTzxdYazMT7Gn93+eDrIbwHXhb13LfBGMhVQ1YdUdaaq5qpqjapemIwvEk6nc+IPGf+k8+e7nVwRvFw3Mk0iHn9D1yCv1B3B7RCuXJZ44s7W+GSLPx1l2N0fRkpy9XjYPX6mDafOf9XyagRYV9tBe/9w3P6Ht7XiUzhr7hSmFecmVIfx/FZh/MefP54O8heA74vIC0BhcKT3e8AXLa+VxXR3jz/H1Pgnr/+KkwLzht9s6GZve39c/ge2tKDAeYvKqSh0J1yHbI5PNvjTUYbd/WGkPVfbPX6mDafOP7M0j3fPKWXYrzyyrTUuf5/HxxM7A1tWX70i9m2lY/VbifEff/54Osi5wBLg18A3gDuB5aq6x/JaWUxVVWKT+o3/+PeX5rm45IQKgNGNPmLxt/V5eHZPBwJcsyK57UWzOT7Z4E9HGXb3h5H2XG33+Jk2nFr/SOf2sR1t9Hsi72Aayf/4jjb6h/2snFbE4qqChMuP5rcS4z/+/PF0kB8DDgIfBbzAbiD6xM0YEZGfishOEdkiIg+JyJRkneF0dMS+3aXxTz7/h5ZX4xR4YV8nh7uHYvKv3dqK16+cNW8KM0rzkio/2+OTaX86yrC7P4yU5OrxsHv8TBtOrX9pTSEn1RTS6/HxZJQdTMP9Hp+ftdtGtpVObhAikt9qjP/488fcQVbV2cBpwMPACuABoFNEHkuyDs8Cy1R1BYFE/vUkfccQz1rPxj/5/NVFObxvYTl+PXq39Hj+niEvjwcv+12zMrnLfpH8VmN3fzrKsLs/rKxU5erxykyV+rjwp6OMbPdfExxFXrs18g6m4f6/7e2ko9/L3LI8TptZklTZkfxWY/zHnz+uRVtVdR/wKrAeeI3AEkJJ/Wmnqs+oqjf48jXGbplqCXYc2jf+9PpHpkk8vbudzrAbScL9f9nexsCwn1NmFLO4MrnLfpH8VmN3fzrKsLs/nFTk6vGwe/xMG069/4zZJcwqzaW1b5gX9nWO6/erjq56cfWKamScjZxiJdvjY/zZ549no5B7ReQg8AdgPnA3MFdVT7ewPp8EnrTQB0Bzc/QVCozf+AHmlOVz5pxShn3KQ2E3koT6B73+0Z9/2ILR43B/KrC7Px1l2N0fSppy9RjsHj/ThlPvd4iMrj3/wJbmY0b8Qv1vHOym/sgglQVuzplfllS5kfypwPiPP388C7euJjAKsTn42KSqPbH8YoxbTd9CYL7c3ZEcLS0t3HjjjbhcLnw+H1deeSU33XRTTHuaqyrt7e0p2f+7vb0dh8NBc3Nz3Pt/x7qnudvtpqGhYdx92pM5ptzcXOrr6+Peez7WYxrxx7v3fKzH5HK5qK+vj3vv+fBjunxRMevru3hkWwsfWDKFno5W8vPzERHq6+uprKzkz5sP0zXoZVFFHmWedtrahpI+Jr/fT319fUJtL5bzNDw8PMYfT9uL5ZgGBwdpaGiw7PsU6Zi8Xi/d3d2WfJ/slCMSJOFcnWieHRgYoL293bKcFB7DgYEBuru7LctJ4e1iYGCA/v5+y3KSybORj2lRTi9T8pzs6xjkybf3suakWaPHFJpn795wEID3LymjseGgJcdk8uz4xzQwMEBzc7Nl/RY75YhoSDzzNkRkKvBe4Gzg74B84EVV/VTMksjeTwCfBdaoasS1ttavX69LlixJyN/e3k5FRUUSNTT+yeL/8mN72NLUy6dOmz46v3jE7/UrN9y/neZeD99cM4+/m2fN/aR2ik8m/OkoI1v9GzdufGvNmjWr4/29RHN1onk2W+OXLf50lGEX/582NXHnhsOcPL2YH1+y8Bj/zpY+/unR3RS4Hdz9kWUU5lizvq1d4mP86fdHy7PxzkFuAnYBtUAdgVHhixOqURARuQj4KnB5tM5xsvT29qZCa/zHoX9k2sTarS14vP4x/hf2ddLc62FmaS7vnltqWZl2ik8m/Okow+7+cFKRq8fD7vEzbTh9/stOrCTP5eDtxh5q247+lz/iH7lR+rITKy3rHIf6U4XxH3/+eOYgPyoiHcAjwMnAX4BTVXVGknX4FVAMPCsim0Tkv5P0HUNNjTVzRY3/+PevnlnM/PJ8Oga8PFvbMepX1dHd9q5ZUYPDgptGRrBTfDLhT0cZdveHksJcHRW7x8+04fT5i3NdXLwkMNIXumpQTU0Nh7qGeHn/EVwO4YMnWXtPqV3iY/zZ449nBHktgSQ7R1WvV9XbrVh4XlUXquosVV0VfHw2WWc4ra3Rd+8xfuMPRUT48MqRG0la8PmV1tZW3jjYTV3nIBUFbt630JqbRkawU3wy4U9HGXb3h5GSXD0edo+facPp9X9oWTWO4NrzzT2eUf+ftwZ2J12zsCyp3UkjYaf4GH92+ONZB/n3qrrf8hqkASuWiDH+yeM/e14ZU4tzaOwe4pW6I4gI920JjB5/aFkVOc64ZiZNiN3ik25/Osqwuz+UTORqu8fPtOH0+quLcjhnfhl+DUxnA+gZVp7ZHbiJa2S1CyuxU3yMPzv81v5Pn6WUl5cbv/HHjNMhXB1M0PdubqbZl8fWpj6KcpxcsqTS0rLAfvFJtz8dZdjdn2nsHj/ThtPvH9kd74ld7XQPelnfAh6fcsasEuaUJbyCS1TsFh/jz7x/UnSQ7Ti0b/yZ9V+wuIJ8l1DbPsDXnw0sObR8aiEFFt40MoId45NOfzrKsLs/09g9fqYNp9+/oKKAU2YUM+T18+d3Wnh8V+CeDyu2lY6E3eJj/Jn3T4oOcklJ8ttUGv/k8r9cdwRP2Haobx3qYV2t9fu92zE+6fSnowy7+zON3eNn2nBm/CM7mP5pczO9w4rLIbT2eiwvB+wZH+PPrH9SdJB9Pp/xG39c3LmhEZ9/7Hsen3LnhkbLy7JjfNLpT0cZdvdnGrvHz7ThzPg7+4cJnTnq9Su3vXwwJQMRdoyP8WfWPyk6yH19fcZv/HHR2jsc1/vJYMf4pNOfjjLs7s80do+facOZ8d/51mHCtyobStFAhB3jY/yZ9U+KDvLUqZF2uTZ+449OVVHkJYaivZ8MdoxPOv3pKMPu/kxj9/iZNpwZfzoHIuwYH+PPrH9SdJCbmpqM3/jj4obV08l1jl02Jtcp3LB6uuVl2TE+6fSnowy7+zON3eNn2nBm/OkciLBjfIw/s/5J0UF++OGHjd/442LNwnJuPms21UVuUKW6yM3NZ81mzULrl5KxY3zS6U9HGXb3Zxq7x8+04cz40zkQYcf4GH9m/ZOig7x27VrjN/64WbOwnLuuXcbQXV/krmuXpaRzDPaNT7r86SjD7v5MY/f4mTacGX/oQISmeCDCjvEx/sz6XZYbsxCv12v8xm/8NvWnowy7+zON3eNn2nDm/GsWlrNmYTnvec97uOuVV1JSBtg3PsafOb+oht9Dmp2sW7euFahP5Hc7Ojoqy8vL2yyukvEbv/GnwZ+OMrLYP2fNmjVVllcoConm2SyOX1b401GG8Ru/8Sfsj5hnbdNBNhgMBoPBYDAY0sGkmINsMBgMBoPBYDDEiukgGwwGg8FgMBgMIZgOssFgMBgMBoPBEMKk6iCLyBdEZJeIbBORn1js/raIHBKRTcHHJVb6Q8r5soioiFRa7P2eiGwJ1v0ZEbF0IUoR+amI7AyW8ZCITLHYf3XwvPpFZLWF3ouCbaZWRL5mlTfo/p2ItIjIViu9If5ZIvI3EdkRjM0/W+zPE5E3RGRz0P8dK/0h5ThF5G0ReSwF7joReSfY7jekwD9FRB4Mtv0dInKm1WVkEybHjus1OTayN2U5NuhPWZ5NdY4NlpHyPJvKHBv02zPPquqkeADnAn8FcoOvqy32fxv4coqPYRbwNIG7zCstdpeEPP8n4L8t9l8AuILPfwz82GL/icAJwPPAaoucTmAvMB/IATYDSy2s89nAKcDWFLWXacApwefFwG6L6y9AUfC5G3gdeFcKjuNfgHuAx1LgrrP6uxTm/z/gU8HnOcCUVJWV6YfJsRO6TY491pnSHBssI2V5NtU5NuhNeZ5NZY4N+m2ZZyfTCPLngB+p6hCAqrZkuD6J8AvgK4DlS4+oanfIy0Kry1DVZ1R1ZKHC14CZFvt3qOouK53A6UCtqu5TVQ9wL3CFVXJVfRHosMoXwX9YVTcGn/cAO4AZFvpVVXuDL93Bh6XtRkRmApcCt1vpTQciUkLgP+c7AFTVo6pHMlurlGJy7DiYHBuRlOZYSG2eTXWODXpTmmftnGMhtXl2MnWQFwNnicjrIvKCiJyWgjI+H7y89TsRKbNSLCKXA4dUdbOV3rAybhWRg8DHgG+mqhzgk8CTKfRbxQzgYMjrBixOfulCROYCJxMYfbDS6xSRTUAL8KyqWuoHbiPQYfFb7B1BgWdE5C0R+QeL3fOBVuDO4OXL20Wk0OIysgmTYycuw+TYsZgcG5s7lXk21TkWbJpnj6ud9ETkr8DUCD+6hcCxlgHvAk4D7heR+Rock7fA/xvgewQawveA/yCQpKyq/78RuISWMOP5VfURVb0FuEVEvg58HviWlf7gZ24BvMDdcVU+Rr/FSIT3bLdwuIgUAX8Gbg4bxUoaVfUBq4LzHR8SkWWqaslcPxG5DGhR1bdE5BwrnBF4j6o2ikg18KyI7AyOOFmBi8Cl3S+o6usi8kvga8C/W+RPOybHJu43OTZykRHeMzk2jFTl2TTlWLBpnj2uOsiqel60n4nI54C1wWT9hoj4gUoCf3kk7Q8r63+BuCe7R/OLyHJgHrBZRCBw6WyjiJyuqk3J+iNwD/A4cSbvifwi8gngMmBNPP9pxupPAQ0E5iSOMBNoTHMdkkJE3AQS992qav1m9UFU9YiIPA9cBFh1M8x7gMslcDNWHlAiInep6sct8qOqjcF/W0TkIQKXfK1K3A1AQ8hoz4MEErdtMTk2MX8ETI4NYHJsHKQgz6Y8x4J98+xkmmLxMPA+ABFZTGAit2XbHorItJCXH8S6TgKq+o6qVqvqXFWdS6BBnBJP4p4IEVkU8vJyYKdV7qD/IuCrwOWq2m+lO4W8CSwSkXkikgNcCzya4TrFjAT+p78D2KGqP0+Bvyo4ooGI5APnYWG7UdWvq+rMYJu/FnjOysQtIoUiUjzynMDooZXf2ybgoIicEHxrDbDdKn8WYnLsOJgcGxGTYycuI2V5NtU5FuydZ4+rEeQJ+B3wOwks9eIBPpHIX9jj8BMRWUXg8lAd8BkL3engR8EG5idwB/dnLfb/CsglcHkF4DVVtawMEfkg8F9AFfC4iGxS1QuTcaqqV0Q+T+CudifwO1XdlnxtA4jIn4BzgEoRaQC+pap3WOUnMDpwHfBOcP4awL+p6hMW+acB/yciTgJ/bN+vqilZJihF1BC4XAmBXHiPqj5lcRlfAO4O/ue/D7jBYn82YXLs+JgcG0aqcyykPM+mOseCybOxkJI8K/+/vTuE3WmN4wD+/W6CbBKbCRSbYDOborMpVJsiUG5VRCPIktmMwhRJQ3Nv0iSFYOZuNymC7RH+x/6PgmBe5/p80jnPe8KvnO++7/s+2/Nz8wsAANbtT9piAQAA36UgAwDAREEGAICJggwAABMFGQAAJgoyfEPbO22vbnoOgP8jGcvvSkEGAICJggwAABMFGSZtj7R90fZD2wfZOp8+bS+3/aftjuX+UtuXbXdudGCAFZGxrIWCDIvlmMpHSe4l2ZXkYZKzy8c3snV87pW2B5NcS3JujPFxE7MCrI2MZU0cNQ2LtieS3E+ydywvRtvnSZ6OMa603Z/kRZL3Se6OMa5valaAtZGxrIlfkGHbniRvx9ffGt98uRhjvE7yLMn+JDd/6WQA6ydjWQ0FGba9S7K3bae1fV8u2p5KcjzJk2z9HQjAj5OxrIaCDNv+TvIpyV9td7Q9k+RYkrTdneR2kgtJzic5vYQ5AD9GxrIa9iDDpO3RJLeSHEjyeFl+leRQkn/HGBeX505mK8wPjzH+28SsAGsjY1kLBRkAACa2WAAAwERBBgCAiYIMAAATBRkAACYKMgAATBRkAACYKMgAADBRkAEAYPIZxBETfTuiFjsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(ncols=2, nrows=3, sharex=True, sharey=True, figsize=(10, 5))\n",
    "for ind, (x_grid, stencil) in enumerate(stencils):\n",
    "    current_ax = ax.flatten()[ind]\n",
    "    current_ax.plot(x_grid, np.array((stencil * h**2)).astype(float), '-o')\n",
    "    current_ax.set_title(f\"{len(x_grid)} point stencil\")\n",
    "    current_ax.set_xlabel('dx')\n",
    "    current_ax.set_ylabel('weight')\n",
    "    current_ax.set_xticks(np.arange(-6, 7))\n",
    "    current_ax.set_yticks(np.arange(-2, 3))\n",
    "    current_ax.set_ylim(-3.5, 2.5)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important feature of this graph is that the more points you add, the smaller the edge coefficients get. A relevant question then is: what is the right stencil size? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite difference stencils at the edge of a mesh "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What I like about the \"automatic\" method described in this post is that it can also help to derive accurate stencils for edges of meshes. Admittedly, this is primarily relevant in 1D, but I remember [working on the 1D heat equation](http://flothesof.github.io/heat-equation-cook-my-meat.html) and wondering how to have accurate stencils at the edges. \n",
    "\n",
    "The example we're going to look at comes from these course notes:\n",
    "http://iacs-courses.seas.harvard.edu/courses/am205/notes/am205_fd_stencil.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is to find a second order accurate stencil for the derivative at the left edge of a 1D mesh. We can build a solution reusing the previously written code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1\\\\0 & h & 2 h\\\\0 & \\frac{h^{2}}{2} & 2 h^{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,      1,      1],\n",
       "[0,      h,    2*h],\n",
       "[0, h**2/2, 2*h**2]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_points = [0, 1, 2]\n",
    "order = 2\n",
    "coef_matrix = ZeroMatrix(order + 1, len(grid_points)).as_mutable()\n",
    "\n",
    "for p, h_coef in enumerate(grid_points):\n",
    "    expansion = TaylorExpansion(h_coef * h, order)\n",
    "    for derivative in range(order + 1):\n",
    "        term =  f(x).diff(x, derivative)\n",
    "        coef_matrix[derivative, p] = expansion.coeff(term)\n",
    "\n",
    "coef_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now use this to find the solution that gives a vector of [0, 1, 0]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{3 f{\\left(x \\right)}}{2 h} + \\frac{2 f{\\left(h + x \\right)}}{h} - \\frac{f{\\left(2 h + x \\right)}}{2 h}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-3*f(x)/(2*h) + 2*f(h + x)/h - f(2*h + x)/(2*h)]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "approx_vector = Matrix([f(x + h_coef * h) for h_coef in grid_points])\n",
    "\n",
    "(coef_matrix.inv() @ Matrix([0, 1, 0])).T @ approx_vector "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! This is the same solution as the PDF."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can go even further by adding points and derive a third order accurate stencil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1\\\\0 & h & 2 h & 3 h\\\\0 & \\frac{h^{2}}{2} & 2 h^{2} & \\frac{9 h^{2}}{2}\\\\0 & \\frac{h^{3}}{6} & \\frac{4 h^{3}}{3} & \\frac{9 h^{3}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,      1,        1,        1],\n",
       "[0,      h,      2*h,      3*h],\n",
       "[0, h**2/2,   2*h**2, 9*h**2/2],\n",
       "[0, h**3/6, 4*h**3/3, 9*h**3/2]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_points = [0, 1, 2, 3]\n",
    "order = 3\n",
    "coef_matrix = ZeroMatrix(order + 1, len(grid_points)).as_mutable()\n",
    "\n",
    "for p, h_coef in enumerate(grid_points):\n",
    "    expansion = TaylorExpansion(h_coef * h, order)\n",
    "    for derivative in range(order + 1):\n",
    "        term =  f(x).diff(x, derivative)\n",
    "        coef_matrix[derivative, p] = expansion.coeff(term)\n",
    "        \n",
    "coef_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{11 f{\\left(x \\right)}}{6 h} + \\frac{3 f{\\left(h + x \\right)}}{h} - \\frac{3 f{\\left(2 h + x \\right)}}{2 h} + \\frac{f{\\left(3 h + x \\right)}}{3 h}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-11*f(x)/(6*h) + 3*f(h + x)/h - 3*f(2*h + x)/(2*h) + f(3*h + x)/(3*h)]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "approx_vector = Matrix([f(x + h_coef * h) for h_coef in grid_points])\n",
    "\n",
    "(coef_matrix.inv() @ Matrix([0, 1, 0, 0])).T @ approx_vector "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, we can derive a third order accurate stencil for the second derivative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{2 f{\\left(x \\right)}}{h^{2}} - \\frac{5 f{\\left(h + x \\right)}}{h^{2}} + \\frac{4 f{\\left(2 h + x \\right)}}{h^{2}} - \\frac{f{\\left(3 h + x \\right)}}{h^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[2*f(x)/h**2 - 5*f(h + x)/h**2 + 4*f(2*h + x)/h**2 - f(3*h + x)/h**2]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(coef_matrix.inv() @ Matrix([0, 0, 1, 0])).T @ approx_vector "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I hope this notebook sheds some light on the linear system method used to derive high order stencils for finite-difference computations. If this sort of thing is interesting for you, check out the [Coursera class](https://www.coursera.org/learn/computers-waves-simulations)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*This post was entirely written using the Jupyter Notebook. Its content is BSD-licensed. You can see a static view or download this notebook with the help of nbviewer at [20200403_StencilsForFiniteDifferences.ipynb](http://nbviewer.ipython.org/urls/raw.github.com/flothesof/posts/master/20200403_StencilsForFiniteDifferences.ipynb).*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}