cubic_spline_interpolation.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "cubic_spline_interpolation.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "Source : \n",
        "https://www.amazon.com/Numerical-Analysis-Richard-L-Burden/dp/1305253663"
      ],
      "metadata": {
        "id": "QxKtk8At8cEU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Cubic Spline Interpolation\n",
        "\n",
        "Polynomial regression was concerned with the approximation of arbitrary functions on closed intervals using a single polynomial. However, high-degree polynomials can osciallate eratically.\n",
        "\n",
        "An aternative approach is to divide the approximation interval into a collection of subintervals and construct a (generally) different approximating polynomial on each sub-interval. This is called **piecewise-polynomial approximation**.\n",
        "\n",
        "\n",
        "The most common piecewise-polynomial approximation uses cubic polynomials between each successive pair of nodes and is called **cubic spline interpolation**. \n",
        "\n",
        "#### <font color='blue'>**Definition**</font>\n",
        "\n",
        "Given a function $f$ defined on $[a,b]$ and a set of nodes $a = x_0<x_1<\\ldots < x_n = b$ , a cubic interpolant $S$ for $f$ is a function that satisfies the following conditions : \n",
        "\n",
        "\\begin{align*}\n",
        "\n",
        "$(a)$ $S(x)$ us a cubic polynomial,deonted $S_j(x)$, on the subinterval $[x_j,x_{j+1}]$ for each $j=0,1,\\ldots, n-1$;\n",
        "\n",
        "$(b)$ $S_j(x_j) \\text{  and } S_j(x_{j+1})$ = $f(x_{j+1})$ for each $j=0,1,\\ldots, n-1$.\n",
        "\n",
        "$(c)$ $S_{j+1}(x_{j+1}) = S_j(x_{j+1}) \\text{ for each }j = 0,1, \\ldots, n-2$ (Implied by (b)).\n",
        "\n",
        "$(d)$ $S'_{j+1}(x_{j+1}) = S'_j(x_{j+1})$ for each $j = 0,1,2,\\ldots, n-2$.\n",
        "\n",
        "$(e)$ $S_{j+1}^\"(x_{j+1}) = S_j^\"(x_{j+1})$ for each $j = 0,1,2,\\ldots, n-2$.\n",
        "\n",
        "$(f)$ One of the following sets of boundary conditions is satisfied: \n",
        "\n",
        "  - (i) $S^\"(x_0) = S^\"(x_n) = 0$ ( **natural**( or free) **boundary**)\n",
        "  - (ii) $S'(x_0) = f'(x_0)$ and  $S'(x_n) = f'(x_n)$ (clamped boundary).\n",
        "\\end{align*}"
      ],
      "metadata": {
        "id": "mufn05bE9wjS"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### <font color='blue'> Construction of a Cubic Spline\n",
        "\n",
        "A spline defined on an interval that is divided into $n$ subintervals will required determining $4n$ constants. To construct the cubic spline interpolant for a given function $f$, the conditions in the definition are applied to the cubic polynomials\n",
        "\n",
        "$S_j(x) = a_j + b_j(x-x_j)+c_j(x-x_j)^2+d_j(x-x_j)^3$\n",
        "\n",
        "for each $j=0,1,ldots,n-1$. Since $S_j(x_j) = a_j = f(x_j)$, condition **(c)** can be applied to obtain. \n",
        "\n",
        "\\begin{align*}\n",
        "a_{j+1} = S_{j+1}(x_{j+1}) = S_j(x_{j+1}) = a_j +b_j(x_{j+1}-x_j) +c_j(x_{j+1}-x_j)^2 + d_j(x_{j+1}-x_j)^3.,\n",
        "\\end{align*}\n",
        "\n",
        "For deach $j = 0,1,2,\\ldots,n-2$.\n",
        "\n",
        "The terms $x_{j+1}-x_j$ are used repeatedly in this development,so it is convenient to introduce the simpler notation.\n",
        "\n",
        "\\begin{align*}\n",
        "h_j = x_{j+1}-x_j\n",
        "\\end{align*}\n",
        "\n",
        "for each $j = 0,1,\\ldots,n-1$. If we also define $a_n = f(x_n)$, then the equation:\n",
        "\n",
        "\\begin{align*}\n",
        "a_{j+1} = a_j + b_jh_j + c_jh_j^2+d_jh_j^3\n",
        "\\end{align*}\n",
        "\n",
        "holds for each $j=0,1,\n",
        "\\ldots,n-1$.\n",
        "\n",
        "In a similar manner, define $b_n = S'(x_n)$ and observed that \n",
        "\n",
        "\\begin{align*}\n",
        "S'_j(x) =b_j +2c_j(x-x_j) +3d_j(x-x)j)^2\n",
        "\\end{align*}\n",
        "\n",
        "implies $S'_j(x_j)=b_j$, for each $j=0,1,\\ldots,n-1$. \n",
        "\n",
        "Applying condition **(d)**  gives\n",
        "\n",
        "\\begin{align*}\n",
        "b_{j+1} = b_j +2c_jh_j + 3d_jh_j^2\n",
        "\\end{align*}\n",
        "for each $j=0,1,\\ldots,n-1$\n",
        "\n",
        "Anoter relationsip between the coefficients of $S_j$ is obtained by defining $c_n = S^\"(x_n)/2$ and applying condition **(e)**. \n",
        "\n",
        "Then , for each $j=0,1,\\ldots,n-1$\n",
        "\n",
        "\\begin{align*}\n",
        "c_{j+1} = c_j + 3d_jh_j\n",
        "\\end{align*}\n",
        "\n",
        "Now, solving the above equation for $d_j$ and substituting in the equation of \n",
        "\n",
        "\\begin{align*}\n",
        "&a_{j+1} \\text{ and } b_{j+1}, \\text{ we get }.\\\\ \n",
        "&d_j = \\frac{1}{3h_j}(c_{j+1} - c{j})\\\\\n",
        "&a_{j+1} = a_j + b_jh_j + \\frac{h_j^2}{3}(2c_j+c_{j+1})\\\\\n",
        "&b_{j+1} = b_j +h_j(2c_j + c_{j+1}) \\tag{1}\n",
        "\\end{align*}\n",
        "\n",
        "\n",
        "The final relationship involving the coeffecients is obtained by solving the appropriate equation. \n",
        "\n",
        "First for $b_j$\n",
        "\n",
        "\\begin{align*}\n",
        "b_j = \\frac{1}{h_j}(a_{j+1}-a_j) - \\frac{h_j}{3}(2c_j + c_{j+1})\n",
        "\\end{align*}\n",
        "\n",
        "amd then, with a reduction of the index,for $b_{j-1}$. This gives\n",
        "\n",
        "\\begin{align*}\n",
        "b_{j-1} = \\frac{1}{h_{j-1}}(a_{j}-a_{j-1}) - \\frac{h_{j-1}}{3}(2c_{j-1} + c_{j})\n",
        "\\end{align*}\n",
        "\n",
        "Substituting these values into the equation (1), with index reduced by (1), gives the **linear system of equations**\n",
        "\n",
        "\\begin{align*}\n",
        "h_{j-1}c_{j-1} + 2(h_{j-1}+h_j)c_j + h_j c_{j+1} = \\frac{3}{h_j}(a_{j+1} -a_j) - \\frac{3}{h_{j-1}}(a_j - a_{j-1}) ----(2)\n",
        "\\end{align*}\n",
        "\n",
        "**for each $j= 1,2,\\ldots,n-1$.**\n",
        "\n",
        "- This system involves only the $\\{c_j\\}_{j=0}^n$ as **unknowns**. \n",
        "- The values of $\\{h_j\\}_{j=0}^{n-1}$ and $\\{a_j\\}_{j=0}^n$ are given, respectively, by the spacing of the nodes $\\{x_j\\}_{j=0}^n$ and the values of $f$ at the nodes.\n",
        "\n",
        "- So once the values of $\\{c_j\\}_{j=0}^n$ are determined, it is a siple matter to find the remainder of the constants $\\{b_j\\}_{j=0}^{n-1}$ and $\\{d_j\\}_{j=0}^{n-1}$\n",
        "\n",
        "Then we can construct the cubic polynomials $\\{S_j(x)\\}_{j=0}^{n-1}$\n",
        "\n",
        "\n",
        "## <font color='blue'>Natural Splines</font>.\n",
        "\n",
        "**Theorem : If $f$  is defined at $a= x_0<x_1,\\ldots <x_n =b$, then $f$ has a unique natural spline interpolant $S$ on the nodes $x_0,x_1,\\ldots,x_n$; that is, a spline interpolant that satisfies the natural boundary conditions $S^\"(a)=0 \\text{ and } S^\"(b)=0$.**\n",
        "\\\n",
        "\\\n",
        "Proof: The boundary conditions in this case imply that $c_n = S^\"(x_n)/2 =0$ and that $0= S^\"(x_0) = 2c_0 + 6d_0(x_0-x_0)$\n",
        "\n",
        "So $c_0 = 0$. The two equations $c_0=0$ and $c_n=0$ together with the equation (2) to produce a linear system described by the vecotr equation $Ax=b$ where $A$ is the $(n+1) × (n+1)$ matrix.\n",
        "\n",
        "\\begin{align*}\n",
        "A = \\left (\\begin{array}{cccc} \n",
        "    1 & 0 & 0 & \\dots  & 0 \\\\\n",
        "    h_0 & 2(h_0+h_1)& h_1 & \\ddots  & \\\\\n",
        "    0 & h_1 & 2(h_1+h_2)& h_2 & \\ddots  & \\\\\n",
        "      \\ddots & \\ddots  & \\ddots & \\ddots  & 0 \\\\\n",
        "      \\ddots & \\ddots  & \\ddots & \\ddots  & 0 \\\\\n",
        "    0 & \\dots & 0 & 0 & 1\n",
        "\\end{array}\\right)\n",
        "\\end{align*}\n",
        "\n",
        "and $b$ and $x$ are the vectors. \n",
        "\n",
        "\n",
        "\\begin{align*}\n",
        "b = \\left (\\begin{array}{c}\n",
        "  0 \\\\\n",
        "   \\frac{3}{h_1}(a_2-a_1) -\\frac{3}{h_0}(a_1-a_0) \\\\\n",
        "   \\vdots \\\\\n",
        "   \\frac{3}{h_{n-1}}(a_n-a_{n-1}) -\\frac{3}{h_{n-2}}(a_{n-1}-a_{n-2}) \\\\\n",
        "  0\n",
        "\\end{array} \\right)\n",
        "\\end{align*}\n",
        "\n",
        "\\begin{align*}\n",
        "x =\\left (\\begin{array}{c}\n",
        "  c_0\\\\\n",
        "  c_1 \\\\\n",
        "   \\vdots \\\\\n",
        "  c_n\n",
        "\\end{array} \\right)\n",
        "\\end{align*}\n",
        "\n",
        "The matrix $A$ is stricly diagonally dominant, that is, in each row the magnitude of the diagonal entry exceeds the sum of the magnitudes of all the other entries in the row. \n"
      ],
      "metadata": {
        "id": "rtJYYCtgJYYS"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## <font color='blue'>Example </font>\n",
        "**Blueprint of the Implementation**\n",
        "\n",
        "Use the data points  $(0,1),(1,e),(2,e^2)$ and $(3,e^3)$ to form a natural spline $S(x)$ that approximates the function $f(x)=e^x$.\n",
        "\n",
        "**Solution:** :\n",
        " We have \n",
        " 1. $n=3$\n",
        "\n",
        " 2. $h_0= h_1 = h_2 = 1$,\n",
        "\n",
        " 3. $a_0 = 1, a_1 = e, a_2 = e^2, a_3 = e^3$\n",
        "\n",
        " 4. Calculate the matrix $A$. \n",
        "\n",
        "From the above derivation, the below condition is true for $j=1,2,\\ldots,n-1$.  In this example . $j = 1,2$ \n",
        "\n",
        "$h_{j-1}c_{j-1} + 2(h_{j-1}+h_j)c_j + h_j c_{j+1}$\n",
        "\n",
        "From the boundary condition for natural spline : $c_0=c_3 =0$\n",
        "\\begin{align*}\n",
        "&j:1\\\\\n",
        "&=h_{0}c_{0} + 2(h_{0}+h_1)c_1 + h_1 c_{2}\\\\\n",
        "&=(1)*c_0 + 2(1+1)*c_1+1*c_2 = 1.c_0 + 4.c_1 +1.c_2\\\\\n",
        "\\\\\n",
        "&j:2\\\\\n",
        "&=h_{1}c_{1} + 2(h_{1}+h_2)c_2 + h_2 c_{3}\\\\\n",
        "&=(1)*c_1+ 2(1+1)*c_2+1*c_3=1.c_1+4.c_2+1.c_3\\\\\n",
        "\\end{align*}\n",
        "\\\n",
        "\\\n",
        "We got 2 of the matrices, ie. $A$ and $x$.\n",
        "\n",
        "\\begin{align*}\n",
        "A = \\left (\\begin{array}{cccc} \n",
        "    1 & 0 & 0 & 0 \\\\\n",
        "    1 & 4& 1 & 0 \\\\\n",
        "    0 & 1 & 4& 1  \\\\\n",
        "    0 & 0 & 0 & 1\n",
        "\\end{array}\\right)\n",
        "\\end{align*}\n",
        "\n",
        "\\begin{align*}\n",
        "x = \\left (\\begin{array}{c} \n",
        "    c_0 \\\\\n",
        "    c_1 \\\\\n",
        "    c_2  \\\\\n",
        "    c_3\n",
        "\\end{array}\\right)\n",
        "\\end{align*}\n",
        "\n",
        "Now calculating $b$.\n",
        "\n",
        "$\\frac{3}{h_j}(a_{j+1} - a_j) - \\frac{3}{h_{j-1}}(a_j - a_{j-1})$\n",
        "\n",
        "Since $c_0 = c_3 = 0$ We will need to solve the linear systems of equations $Ax=b$ for the remaining 2 values of $c: c_1, c_2$. \n",
        "\n",
        "\\begin{align*}\n",
        "j:1: \n",
        "&=\\frac{3}{h_1}(a_{2} -a_1) - \\frac{3}{h_0}(a_1 - a_0)\\\\\n",
        "&= 3.1.(e^2-e)-3(e-1)\\\\\n",
        "&= 3(e^2 -2e +1)\\\\\n",
        "\\\\\n",
        "\\\\\n",
        "j:2:\n",
        "&=\\frac{3}{h_2}(a_3 -a_2) - \\frac{3}{h_2}(a_2- a_1)\\\\\n",
        "&= 3.1.(e^3-e^2)-3(e^2-e)\\\\\n",
        "&= 3(e^3 -2e^2 +e)\\\\\n",
        "\\end{align*}\n",
        "\\\n",
        "\\\n",
        "\\begin{align*}\n",
        "\n",
        "\n",
        "\\end{align*}\n",
        "\n",
        "\\begin{align*}\n",
        "Ax &= B \\\\\\\\\n",
        "\\begin{bmatrix}\n",
        "    1 & 0 & 0 & 0 \\\\\n",
        "    1 & 4& 1 & 0 \\\\\n",
        "    0 & 1 & 4& 1  \\\\\n",
        "    0 & 0 & 0 & 1\n",
        "\\end{bmatrix}\n",
        "\\begin{bmatrix}c_0 \\\\ c_1 \\\\ c_2  \\\\c_3\\end{bmatrix} &= \n",
        "\\begin{bmatrix}\n",
        "0 \\\\\n",
        "    3(e^2 -2e +1) \\\\\n",
        "    3(e^3 -2e^2 +e)  \\\\\n",
        "   0\n",
        "\\end{bmatrix}\n",
        "\\end{align*}"
      ],
      "metadata": {
        "id": "vX9Ke4pa8Do3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "1. calcualte the matrices. A and b\\\n",
        "    *function to calculate the exp \\\n",
        "    *function to calculate for each entry of A\\\n",
        "    *function to calculate for each entry of b\n",
        "2. Solve the linear systems of equations, using the Gaussian Elimination.\n",
        "3. By Substitution calculate the remaining constants of the cubic splines, i.e co-ordinates of the quadratic and cubic terms. (b & d)\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "a6CC_vjVSDhv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "###########################################################\n",
        "# function definition\n",
        "###########################################################\n",
        "def f(x):\n",
        "  return np.exp(x)\n",
        "\n",
        "###########################################################\n",
        "# Evaluate f(x) for given x's and those\n",
        "# will be the values of constants a, \n",
        "# in the cubic spline equation\n",
        "###########################################################\n",
        "\n",
        "def vec_a(x):\n",
        "  m = len(x)\n",
        "  a = np.zeros([m])\n",
        "  for i in range(0,m):\n",
        "    a[i] = f(x[i])\n",
        "  \n",
        "  return a\n",
        "\n",
        "###########################################################\n",
        "# Formula Evaluation to get matrix A.\n",
        "###########################################################\n",
        "\n",
        "def matrix_A(h,n):\n",
        "  A = np.zeros([n,n])\n",
        "\n",
        "  A[0,0], A[n-1,n-1] = 1,1\n",
        "  \n",
        "  for i in range(1,n-1):\n",
        "    # hardcoded for this example, will need to check if same is applicable for others\n",
        "    # except changing the value of h, anything else also needs to be changed?\n",
        "    A[i,i-1] = h[0]\n",
        "    A[i,i] = 2*(h[0]+h[1])\n",
        "    A[i,i+1] = h[1]\n",
        "  return A\n",
        "\n",
        "###########################################################\n",
        "# Formula Evaluation to get value vector B.\n",
        "###########################################################\n",
        "\n",
        "def matrix_B(a,h,n):\n",
        "\n",
        "  B = np.zeros([n])\n",
        "  for i in range(1,n-1):\n",
        "    B[i] = ((3/h[i])*(a[i+1]-a[i])) - ((3/h[i-1])*(a[i]-a[i-1]))\n",
        "  return B\n",
        "\n",
        "###########################################################\n",
        "# Formula Evaluation to get vector b. \n",
        " #(coeffecients of the quadractic term of the cubic splines)\n",
        "###########################################################\n",
        "\n",
        "def cal_b(a,h,c,n):\n",
        "  b  = np.zeros([n-1])\n",
        "\n",
        "  for i in range(0,n-1):\n",
        "    b[i] = ((1/h[i]) * (a[i+1]-a[i])) - ((h[i]/3)*(2*c[i]+c[i+1]))\n",
        "  return b\n",
        "\n",
        "###########################################################\n",
        "# Formula Evaluation to get vector d.\n",
        "#(coeffecients of the cubic term of the cubic splines)\n",
        "###########################################################\n",
        "\n",
        "def cal_d(h,c,n):\n",
        "  d  = np.zeros([n-1])\n",
        "\n",
        "  for i in range(0,n-1):\n",
        "    d[i] = (1/3*h[i])*(c[i+1]- c[i])\n",
        "  return d\n"
      ],
      "metadata": {
        "id": "RrJPMOS37PgS"
      },
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "\"\"\" \n",
        "Implementing the Gaussian Elimination to find unique solutions to the given system of linear equations.\n",
        "The algorithm implementation is largely divide into 2 parts \n",
        "1. Convert the original matrix to a upper triangular matrix, by matrix elimination.\n",
        "2. Calcuate value of x's in reverse order, by Backsubstituion\n",
        "\"\"\"\n",
        "def gaussian_elimination(a,b):\n",
        "\n",
        "  # Create the augmented matrix.\n",
        "  ab = np.hstack([a,b])\n",
        "  \n",
        "  #First loop for the Pivots.\n",
        "  n = len(a)\n",
        "  ##############################################################################\n",
        "  \"\"\"\n",
        "  Part 1 of the Implementation : converting the original matrix to a upper \n",
        "  triangular matrix, by matrix elimination.\n",
        "  STEPS:\n",
        "  (i). Outer for loop to loop across each of the pivots.\n",
        "  (ii). logic for row exchange, to find the smallest pivot for each loop around a new pivot \n",
        "  - Note: Point (ii) is not needed for tridiagonal matrices. like this example.\n",
        "  (iii). Elimination step to arrive at the upper triangular matrix reduced matrix.\n",
        "  \"\"\"\n",
        "  ##############################################################################\n",
        "  for i in range(0,n):\n",
        "    '''\n",
        "    # find the minimum row value for the respective column i \n",
        "    min = 999999999\n",
        "    for j in range(i+1,n):\n",
        "      if  ab[j,i] != 0 and ab[j,i] < min:\n",
        "        min = ab[j,i]\n",
        "        row = j\n",
        "    \n",
        "    # Exception Handling in case, min = 0\n",
        "    if min == 0:\n",
        "      print(\"No integer can be found, hence no unique solution\")\n",
        "      return \n",
        "\n",
        "    # swapping the rows \n",
        "    temp = np.copy(ab[i,:])\n",
        "    ab[i,:] = ab[row,:]\n",
        "    ab[row,:] = temp\n",
        "    '''\n",
        "    # Once,the swapping is done, we start the elimination loop for rows below the Pivot\n",
        "    for j in range(i+1,n):\n",
        "      m = (ab[j,i]/ab[i,i]) * ab[i,:]\n",
        "      ab[j,:] = ab[j,:] - m\n",
        "\n",
        "  ##############################################################################\n",
        "  \"\"\"\n",
        "  Part 2 of the Implementation :  Calcuate value of x's in reverse order, \n",
        "  by Backsubstituion\n",
        "  \"\"\"\n",
        "  ##############################################################################\n",
        "  x = np.zeros(a.shape[0])\n",
        "  if ab[n-1,n-1] == 0:\n",
        "    print(\"No unique solution\")\n",
        "    return  \n",
        "  \n",
        "  # Get the last value of x, before starting the loop, to calculate rest of the\n",
        "  # elements.\n",
        "  x[n-1] = ab[n-1,n]/ab[n-1,n-1]\n",
        "\n",
        "  for i in range(n-2,-1,-1):\n",
        "    sum = 0\n",
        "    for j in range(n-1,i,-1):\n",
        "      sum += ab[i,j]*x[j] \n",
        "    x[i] = (ab[i,n] - sum)/ab[i,i]\n",
        "\n",
        "  return x"
      ],
      "metadata": {
        "id": "i5jYyK6VvvzF"
      },
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# function to create the splines \n",
        "def splines(a,b,c,d,point_x,x):\n",
        "\n",
        "  spline = []\n",
        "\n",
        "  for i in range(0,len(x)):\n",
        "     val = a + b*(x[i]-point_x) +c*(x[i]-point_x)**2 + d*(x[i]-point_x)**3\n",
        "     spline.append(val)\n",
        "    \n",
        "  return spline"
      ],
      "metadata": {
        "id": "BmyDT1Lgu9Ex"
      },
      "execution_count": 24,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "##############################################################################\n",
        "\"\"\"\n",
        "Function for calculating the \n",
        "1. Unknown coeffecients of the spline equation\n",
        "2. Corresponding function values \n",
        "\"\"\"\n",
        "##############################################################################\n",
        "def cubic_pline(h,x,n):\n",
        "\n",
        "  #########################################################\n",
        "  \"\"\"\n",
        "  Part-1 of the Implementation\n",
        "  - Calculating all the needed constants to build our spline equations\n",
        "  \"\"\"\n",
        "  #########################################################\n",
        "\n",
        "  A = matrix_A(h,n)\n",
        "  a = vec_a(point_x)\n",
        "  B = matrix_B(a,h,n)\n",
        "  B = B.reshape(n,1)\n",
        "  c = gaussian_elimination(A,B)\n",
        "  b = cal_b(a,h,c,n)\n",
        "  d = cal_d(h,c,n)\n",
        "\n",
        "   #########################################################\n",
        "  \"\"\"\n",
        "  Part-2 of the Implementation\n",
        "  - Generating y-values for corresponding intervals\n",
        "  For example, in this example we have 4 knot points\n",
        "  @[0,1,2,3] . Hence, we will get 3 cubic spline equations.\n",
        "  and for each we get corresponding y values. \n",
        "  \"\"\"\n",
        "  #########################################################\n",
        "\n",
        "  # Create equally spaced points between the knot points\n",
        "\n",
        "  p1 = np.linspace(0, 1, num=10)\n",
        "  p2 = np.linspace(1,2,num=10)\n",
        "  p3 = np.linspace(2,3,num=10)\n",
        "\n",
        "  p = [p1,p2,p3]\n",
        "\n",
        "  y = []\n",
        "  for i in range(0,len(p)):\n",
        "    _a = a[i]\n",
        "    _b = b[i]\n",
        "    _c = c[i]\n",
        "    _d = d[i]\n",
        "    _point_x = point_x[i]\n",
        "    ys = splines(_a,_b,_c,_d,_point_x,p[i])\n",
        "    y.extend(ys)\n",
        "\n",
        "  p = np.concatenate([p1,p2,p3])\n",
        "  return p,y\n"
      ],
      "metadata": {
        "id": "MKzYxB9ov4ZH"
      },
      "execution_count": 19,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "h = np.array([1,1,1])\n",
        "point_x = [0,1,2,3]\n",
        "n=4\n",
        "cubic_x,cubic_y= cubic_pline(h,point_x,n)"
      ],
      "metadata": {
        "id": "4pJ2VRY-o0zg"
      },
      "execution_count": 25,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Continuing our solution** \n",
        "\\begin{align*}\n",
        "\\begin{bmatrix}\n",
        "    c_0 \\\\\n",
        "    c_1 \\\\\n",
        "    c_2  \\\\\n",
        "    c_3\n",
        "\\end{bmatrix}\n",
        "&= \n",
        "\\begin{bmatrix}\n",
        "0 \\\\\n",
        "    0.75685264 \\\\\n",
        "    5.83006675  \\\\\n",
        "   0\n",
        "\\end{bmatrix}\n",
        "\\end{align*}\n",
        "\n",
        "Now calculating $b_0,b_1,b_2$ and $d_0,d_1,d_2$ which are the coeffecients of the quadratic and the cubic terms.\n",
        "\n",
        "$b_j = \\frac{1}{h_j}(a_{j+1} -a_j) - \\frac{h_j}{3}(2c_j + c_{j+1})$\n",
        "\n",
        "\\begin{align*}\n",
        "&b_0 = \\frac{1}{h_0}(a_1-a_0) - \\frac{h_0}{3}(2c_0 + c_1)\\\\\n",
        "&= (e-1) -(1/3)(0+0.7568)\\\\\n",
        "&=1.718-0.252\\\\\n",
        "&=1.466\n",
        "\\\\\n",
        "\\\\\n",
        "&b_1 = \\frac{1}{h_1}(a_2-a_1) - \\frac{h_1}{3}(2c_1 + c_2)\\\\\n",
        "&= (e^2-e) -(1/3)(2*0.7568 + 5.830)\\\\\n",
        "&=4.670-2.4478\\\\\n",
        "&=2.222\n",
        "\\\\\n",
        "\\\\\n",
        "&b_2 = \\frac{1}{h_2}(a_3-a_2) - \\frac{h_2}{3}(2c_2 + c_3)\\\\\n",
        "&= (e^3-e^2) -(1/3)(2*5.830+ 0)\\\\\n",
        "&=12.696-3.8866\\\\\n",
        "&=8.80933\n",
        "\\end{align*}\n",
        "\n",
        "$d_j = \\frac{1}{3h_j}(c_{j+1} - c_j)$\n",
        "\n",
        "\\begin{align*}\n",
        "&d_0 = \\frac{1}{3h_0}(c_1-c_0)\\\\\n",
        "&= (1/3)(0.7568)\\\\\n",
        "&=0.25226\\\\\n",
        "\\\\\n",
        "&d_1 = \\frac{1}{3h_1}(c_2 - c_1)\\\\\n",
        "&= (1/3)(5.830 - 0.7568)\\\\\n",
        "&=1.6910\\\\\n",
        "\\\\\n",
        "&d_2 = \\frac{1}{3h_2}(c_3 - c_2)\\\\\n",
        "&= (1/3)(0-5.830)\\\\\n",
        "&=-1.9433\\\\\n",
        "\\end{align*}\n",
        "\n",
        "The natural cubic spline is described piecewise by :  \n",
        "\n",
        "  \n",
        "\\begin{align*}\n",
        "S(x) =\n",
        "\\begin{cases}\n",
        "  1+1.46600x+0.2522x^3 & \\text{ for  }x\\in [0,1]\\\\\n",
        "2.71828 + 2.222(x-1)+0.75685(x-1)^2 +1.69107(x-1)^3 &\\text{ for  }x\\in [1,2] \\\\\n",
        "      7.38906 + 8.80933(x-2)+5.83(x-2)^2 -1.9433(x-2)^3 &\\text{ for  }x\\in [2,3]\n",
        "\\end{cases} \n",
        "\\end{align*}"
      ],
      "metadata": {
        "id": "3wxL5qn_txFj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "y0 = 1 +1.46600*x +0.2522*x**3\n",
        "y1 = 2.71828 + 2.222*(x−1)+ 0.75685*(x−1)**2 +1.69107*(x−1)**3\n",
        "y2 = 7.38906 + 8.80933*(𝑥−2)+ 5.83(𝑥−2)**2 − 1.9433*(𝑥−2)**3 "
      ],
      "metadata": {
        "id": "eZn_bkPnDYEm"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from numpy.core.function_base import linspace\n",
        "# Try to plot the splines against the function \n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.figure(figsize=(8,6))\n",
        "x = np.linspace(0, 3.0, num=30)\n",
        "y = np.exp(x)\n",
        "plt.plot(x,y,label='y=exp(x)',color='blue',linewidth=1 )\n",
        "plt.plot(cubic_x,cubic_y,label ='Cubic Spline Approx',linewidth=2,color='green')\n",
        "#knot Points Scatter\n",
        "x_knot = [0,1,2,3]\n",
        "y_knot = [1,2.718,7.38,20.08]\n",
        "plt.scatter(x_knot,y_knot,color='red',s=100,label='Knot points')\n",
        "leg = plt.legend(loc='lower right')\n",
        "plt.title(\"Cubic Spline Interpolation for f(x)=e^x\")\n",
        "plt.grid()\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 390
        },
        "id": "wEq3PzA4DV72",
        "outputId": "ea64b69d-346e-4889-9bb9-112d6a2f3d3d"
      },
      "execution_count": 43,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAF1CAYAAAAnXamsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZyN5f/H8dc1Y6xjXyYh4osWRdFOKHvEj1LKkrIUlSXZs7ZIJNlC1oRQKrJmSaQwouwSMvbdjHXMXL8/7oNZzjDMmblneT8fj/Nwzn1d931/znXG+Zzruq/7vo21FhEREUme/NwOQEREROKmRC0iIpKMKVGLiIgkY0rUIiIiyZgStYiISDKmRC0iIpKMKVGLq4wxe4wxVeIoq2CM2Z7E8bxijFkZ5XWYMaZoUsaQ3BhjJhpj3k/A+onShsaYksaYDcaYUGPM2z7YnjHGTDDGnDTGrPEsy2CM2WKMyR+P9d8yxnyc0DhusI+yxpjexpjbE3M/krwoUUuCGWNeMsas83whHzTGzDfGlE/odq21v1prS95CPDmMMeONMYc8X+I7jDFdbzGGQGvtv7ey7g1iXG6MaRHPuglKlEnJ2/tKrDYEOgPLrLVZrbWf+2B75YGqQEFr7cOeZa2AFdbag/FYfyzwsjEmnw9iicUYcxewAHgaWGiMyZEY+5HkR4laEsQY0xH4DPgQCALuAEYCdV0MawgQCNwNZAeeBf5xMR5XGWPSuR1DIikMbL6VFeNok8LAHmvt2SjLXge+is82rbUXgPlA01uJ6XqMMQU92+4GVASWAT8aYzL6el+SDFlr9dDjlh44STAMeP46dSYC70d5XQkIifJ6D86XzxbgJDAByBhH3ULAd8BR4DgwPI59bgLqXScmC7wN/AscAz4B/DxlrwArY9T9X5T3MgL4CQgF/gCKRal7F7AYOAFsBxpeJ4blQIuo7xN4BzgCHASae8paAeHAJU9bz/Esvx341tMWu4G3o2y7DzALmAKcAVpEWfaNJ/b1QOko69ztiekUTvJ71ttnCOQE5nr2e9LzvKCn7AMgArjgiXW4lzbMDkz2rL8X6Bmz7YFBnm3vBmrG0X5LY+yrRDy2vQrnR9xxovxNespf82wrwrO9vjg/Os8D6Tx10gMbgLc8r/092+wVZTsv4/Tyb/b/0vU+z1zARqBxjHUGAN8D/p7X84DBUcqnA+Pd/p7QI+EP1wPQI+U+gBrA5StfZHHUufol73ldidiJehNOEs7l+eJ7P2Zdz5fiRs8XbRYgI1A+jn1+iZNsmgPFvZRbnB5JLs+X8Q6uJc1XuH6iPg48DKQDvgame8qyAPs8+0wHPIDzI+CeOGJcTvREfRnoBwQAtYBzQM442tAPCAZ6eZJHUZwfHdU95X1wkns9T91MUZY959lHJ09CCPA8/gG6e7b3FE4yLxlz/0BuoAGQGcgKzAS+9/a+4mjDycAPnnWLeNr+tShtHw609HzebwAHAHOjNoznti8Db3k+n0xethfzs38G2ByjTimcHxF3Az2A3/EkSk/5g8CJKK/n4vz48faYG5/P8yb+P96G80PvKZwfDP8CWd3+ntAj4Q8NfUtC5AaOWWsvJ3A7w621+6y1J3B6ZY281HkYp9fxrrX2rLX2grV2pZd64HwZfw28CWwxxvxjjKkZo87H1toT1tr/cIbuve3Tm9nW2jWe9/w1UMazvDbOsOkEa+1la+2fOD2k5+O53XCgn7U23Fo7D6dXF9fx+YeAvNbaftbaS9Y5/jsWeDFKndXW2u+ttZHW2vOeZcHW2lnW2nDgU5wfO496HoHAAM/2luIkmFhtYq09bq391lp7zlobivN5VYzPGzTG+Hti7GatDbXW7gEGA02iVNtrrR1rrY0AJgH5cQ6p+GLbB6y1wzyfz3lv24khB84PlqustZuA93F6sp2AJp5YrwjF6dlfqV/bWpsjjkdtT7X4fJ43ZK09hPPjZhIwFGjq+YwkhVOiloQ4DuTxwTHQfVGe78VJyDEVwvkSv+GPAmvteWvth9basjg/JmYAM40xuW5yn94civL8HE6CA+f45iPGmFNXHji9mtviud3jMd5b1G3HVBi4Pca+uhM9oe3zst7VZdbaSJzh9ts9j32eZVfsBQrE3IAxJrMxZrQxZq8x5gywAsjhSZQ3kgen9773Ovu52r7W2nOep3G1w81u21ubXM9JnN55TJNwPoN51tqdMcqyAqdvcj/x+Tzjaw7OaMT26/yQlRRGiVoSYjVwEWeINS5ncYZJr/CWuApFeX4HznBnTPuAO272R4G19gzORLcswJ03uc+bsQ/4JUaPKdBa+0YCtwvO0HHMfe2Osa+s1tpa11kHorxnY4wfUBDnfR8ACnmWXXEHsN/LNt7B6ek/Yq3NBjx5ZZPX2e8Vx3BGDgrHYz83Kz7bvtlbBf4F3Onlb24kzohDdS9nN9yNc4gGAM8ZEGFxPOZ7qsXn84yvD4CtQH5jTHxHiSSZU6KWW2atPY1zXG2EMaaep7cVYIypaYwZ6Km2AahljMlljLkNaO9lU22NMQU9Pd4eOBOeYlqDM8lqgDEmizEmozHmCW9xGWPeM8Y8ZIxJ75kV2w7nmGDUc7LfNcbkNMYU8pR72+fNmAuUMMY08bRBgCeGuxO4XYDDOMctr1gDhBpjuhhjMhlj/I0xpYwxD91gO2WNMfU9iac9zo+s33EmxZ0DOnvirgTUwZmMFFNWnAlWpzyfV+8bxHqVZ4h4BvCBMSarMaYw0BFn0luCJMa2rbUhOMfur5yqhTGmCVAW53j228AkY0zUHn9FnNnZV7ZR0/ODzdvjyuGYW/08ozHGPIkzR6Ip0AwYZoyJNSoiKY8StSSItXYwzhdiT5wZq/twjg1/76nyFU4PYw+wCO8Jcaqn7F9gF84xwJj7icBJHv8D/sMZtn0hrrBwZo8fw+ktVgWesdaGRanzA84Eng04s7jHxePtxslzLLAaznHFAzhDuB8DGRKyXY9xwD2eYdHvPW1RG+f4+G6c9/klUY6NxuEHnDY7iXPstr7nmPglnLat6dnWSJzjm9u8bOMznMlpx3CS/IIY5UOB5zwXDfF2bvNbOKMs/+LM8J4KjL9B3PGVGNsejec4tzHmDpz339RaG2atnQqsw5ngiOdHYS2cofF4S8DneZUxJhvOZLo3rbX7rbW/4vzdTDDGmOuvLcmdsfZmR4NEUjZjjMWZDZ5mzq02xvTBmXnd2O1YUhJjTAbgT+Bpe4OLnhhj3gIKWWs7J0lwkmak1gshiIgkmLX2InBPPOsOS+RwJI3S0LeIiEgypqFvERGRZEw9ahERkWRMiVpERCQZS5aTyfLkyWOLFCnis+2dPXuWLFmy+Gx7KZ3a4xq1RXRqj2vUFtGpPaLzdXsEBwcfs9bm9VaWLBN1kSJFWLdunc+2t3z5cipVquSz7aV0ao9r1BbRqT2uUVtEp/aIztftYYzZG1eZhr5FRESSMSVqERGRZEyJWkREJBlTohYREUnGlKhFRESSMSVqERGRZEyJWkREJBlTohYREUnGlKhFRETiY9cuaNMGsmWD4GDn3zZtnOWJSIlaRETkRubPh/vvhy+/hNBQZ1loqPP6/vud8kRyw0RtjClkjFlmjNlijNlsjGnnWZ7LGLPYGLPT82/OONZv5qmz0xjTzNdvQEREJFHt2gXPPQfnzkF4ePSy8HBn+XPPJVrPOj496svAO9bae4BHgbbGmHuArsASa21xYInndTTGmFxAb+AR4GGgd1wJXUREJFkaPPhqgj4XAHVe8GdzxJ7odcLDYciQRNn9DRO1tfagtXa953kosBUoANQFJnmqTQLqeVm9OrDYWnvCWnsSWAzU8EXgIiIiSWLKlKuJ+qMnDHPvjmBg2I/YqHXCw+GrrxJl9zd1jNoYUwR4APgDCLLWHvQUHQKCvKxSANgX5XWIZ5mIiEjKEBYGwK6c8HF5A0CnwDqYOOr5Wrxvc2mMCQS+Bdpba88Ycy1Ea601xtg4V47f9lsBrQCCgoJYvnx5QjYXTVhYmE+3l9KpPa5RW0Sn9rhGbRFdmm6PwYMhMpLu58cTHrGFMmGVuPOuCiwfdGf0ev7+kBhtZK294QMIABYCHaMs2w7k9zzPD2z3sl4jYHSU16OBRjfaX9myZa0vLVu2zKfbS+nUHteoLaJTe1yjtoguTbfHG2/YuXf7W/pg/bplsSGB2GWDBlkL1x4BAda2bXvLuwDW2ThyYnxmfRtgHLDVWvtplKIfgSuzuJsBP3hZfSFQzRiT0zOJrJpnmYiISIpwof2btKsWCUCDZaUo4G2EOyAAOnRIlP3H5xj1E0AT4CljzAbPoxYwAKhqjNkJVPG8xhhTzhjzJYC19gTQH1jrefTzLBMREUkRBh36jl05LemOlGTc+k3RCwMCIHNmmDULihVLlP3f8Bi1tXYlxD5m7vG0l/rrgBZRXo8Hxt9qgCIiIm7Ze2ovH/76IQBNswwga8ufr83uzpYNmjRxetKJlKThJiaTiYiIpDXvLHqH85fPE7D9BT4ZXg9y1YPhw51JY6dPJ0kMuoSoiIiIF4t3Lebbrd8SYLPQON8gcuVyJw71qEVERGK4FHGJt+a/BUC6397jvdEFXYtFiVpERCSGz37/jO3Ht5PXlKBirg7ceeeN10ksStQiIiJR7D+zn36/9APAzhtG1xHpXY1HiVpERCSKTos7cTb8LOWy/B9ZM1ejbFl349FkMhEREY/le5YzfdN0MqbLSNisIbz7rtsRKVGLiIgAEB4Rzpvz3gSgYf7upDtbmBrJ4H6PStQiIiLAiLUj2Hx0M0VzFmXf9Hd55x0wcV3uKwkpUYuISJp3KOwQvZf3BuDt4kPZsSUjL73kclAeStQiIpLmdfm5C2cunqF2idqsmVKbt9+G9O5O9r5KiVpERNK0Vf+tYvLGyWTwz0Dn+z5jwQJo3drtqK5RohYRkTQrIjKCtvPaAtD5ic7MHl+M5s0he3aXA4tC51GLiEia9cW6L9h4eCOFsxfm9VJdKdUANm50O6ro1KMWEZE06ejZo/Rc1hOAIdWHMHlcZmrXhkKFXA4sBvWoRUQkTeq2pBunLpyiWrFq1LyzHm0/h/nz3Y4qNvWoRUQkzVmzfw3j/hxHgF8An9f4nKlTDffdB6VLux1ZbOpRi4hImhJ1AlnHxzpSPFdJ/m8QfP65y4HFQT1qERFJU8b/OZ51B9ZRIGsBej7Zk/nzIUMGePpptyPzTolaRETSjOPnjtN1SVcABlcbTGD6QD75BN59N3lcLtQbJWoREUkzui/pzonzJ3j6zqdpeG9D1q6F3bvh+efdjixuStQiIpImrNm/hrHrxxLgF8DwWsMxxvDJJ9C+PQQEuB1d3DSZTEREUr0rE8gslo6PdeSuPHfx77+wdCmMG+d2dNenHrWIiKR6Y9ePZd2BdRTMVpCeTzoXOfn0U2jZErJmdTm4G1CPWkREUrWjZ4/SfUl3wLkCWWD6QA4ehKlTYfNml4OLB/WoRUQkVeu2pBsnL5ykatGqNLi7AQADBsArr0D+/O7GFh/qUYuISKr1e8jvV69ANqzmMIwx7N8PX30FW7e6HV38qEctIiKpUkRkBG1+agNAp8c7UTJPSQA++gheew2CgtyMLv7UoxYRkVRpdPBo/jz0J3dkv4MeFXoAsG+fc2x62zaXg7sJ6lGLiEiqc+TsEXosdZLzZ9U/I0v6LAB8+CG0agX58rkZ3c25YY/aGDMeqA0csdaW8iz7BijpqZIDOGWtLeNl3T1AKBABXLbWlvNR3CIiInHq+nNXTl04RY3/1aDeXfUA2LsXZsyA7dtdDu4mxWfoeyIwHJh8ZYG19oUrz40xg4HT11m/srX22K0GKCIicjN+2/cbEzZMIL1/ej6v8TnGcxHvDz6A11+HPHlcDvAm3TBRW2tXGGOKeCszzrtvCDzl27BERERu3uXIy1dvYdn58c4Uz10ccK7n/e23sHOnm9HdmoQeo64AHLbWxvXWLbDIGBNsjGmVwH2JiIhc1xfrvmDDoQ0Uzl6YbhW6XV3+/vvQti3kyuVicLfIWGtvXMnpUc+9cow6yvJRwD/W2sFxrFfAWrvfGJMPWAy8Za1dEUfdVkArgKCgoLLTp0+/mfdxXWFhYQQGBvpseymd2uMatUV0ao9r1BbRpYT2OHHpBE3XNOVsxFn639uf8nnKA7B/fybatHmQKVP+IGvWyz7Zl6/bo3LlysFxzuOy1t7wARQBNsVYlg44DBSM5zb6AJ3iU7ds2bLWl5YtW+bT7aV0ao9r1BbRqT2uUVtElxLao9nsZpY+2Fpf17KRkZFXlzdtam2fPr7dl6/bA1hn48iJCTmPugqwzVob4q3QGJMF8LPWhnqeVwP6JWB/IiIiXq38byWTNk4ig3+GaBPIduyAefPgn39cDjABbniM2hgzDVgNlDTGhBhjXvMUvQhMi1H3dmPMPM/LIGClMWYjsAb4yVq7wHehi4iIRJ9A1uWJLhTLVexqWb9+0K4dZM/uVnQJF59Z343iWP6Kl2UHgFqe5/8CpRMYn4iIyHWNXDuSvw7/RZEcRehavuvV5Vu3wqJFMHKki8H5gK5MJiIiKdahsEO8t+w9AD6v8TmZAjJdLevXDzp0gGzZ3IrON3StbxERSbE6L+7MmYtnqF2iNnVK1rm6fPNmWLoUxoxxMTgfUY9aRERSpBV7V/DVX1+RwT8DQ2sMjVbWty+88w5kzepScD6kHrWIiKQ44RHhVyeQdSvfjaI5i14t+/tvWLECJkxwKzrfUo9aRERSnGFrhrHpyCaK5ixK5yc6Ryvr0wc6d4YsWdyJzdfUoxYRkRQl5EwIvZf3BmBYzWHRJpBt2ACrV8NXX7kVne+pRy0iIilK+wXtCbsURv2761OreK1oZX36QJcukDmzO7ElBvWoRUQkxZi/cz7fbv2WLAFZ+Kz6Z9HKgoNh7VqYNi2OlVMo9ahFRCRFOB9+/uoEsr6V+lIoe6Fo5X36QLdukCmTl5VTMPWoRUQkRfjg1w/YfWo39+W7j7cfeTta2Zo1zvHpmTNdCi4RqUctIiLJ3rZj2xi4aiAAX9T+ggD/gGjlffpA9+6QMaMLwSUy9ahFRCRZs9bS5qc2hEeG0+KBFjxe6PFo5atXO1cimz3bpQATmXrUIiKSrE39eyrL9iwjT+Y8DKgyIFZ5r17QowdkyOBCcElAiVpERJKtUxdO0XFRRwAGVhlI7sy5o5UvWAD//QfNm7sRXdJQohYRkWSrx5IeHDl7hPJ3lKdZmWbRyi5fdq7n/cknEBAQxwZSASVqERFJltbuX8uodaNI55eOUc+Mws9ET1njxkFQENSpE8cGUglNJhMRkWQnIjKC1396HYulw6MdKJWvVLTyM2egd2+YPx+McSnIJKIetYiIJDsj145k/cH1FMpWiF4Ve8Uq/+gjqFkTHnjAheCSmHrUIiKSrBwMPUjPZT0B56YbgekDo5Xv3QtjxsBff7kRXdJTj1pERJKVjos6cubiGeqUqEPdu+rGKu/WDd56CwoUcCE4F6hHLSIiycbiXYuZvmk6mdJl4vOan8cq/+MPWLECxo51ITiXqEctIiLJwoXLF67edKNXxV4UyVEkWrm10LEjvP8+ZMniQoAuUaIWEZFkYeCqgew8sZN78t5Dx8c6xiqfNQvOnYOmTV0IzkUa+hYREdf9c+IfPvz1QwBG1hpJev/00covXoQuXeDLL8EvjXUx09jbFRGR5MZaS9t5bbkYcZGmpZtSsUjFWHWGDYNSpeCpp1wI0GXqUYuIiKtmbpnJol2LyJkxJ59U/SRW+dGjMGAArFrlQnDJgHrUIiLimjMXz9B+QXsABlQZQL4s+WLV6dsXXnoJSpZM6uiSB/WoRUTENb2W9eJg2EEeLfgoLR5sEat82zb45hvYutWF4JIJ9ahFRMQVfx78k2FrhuFn/LzedAPg3Xeha1fIk8eFAJMJ9ahFRCTJRURG0HJOSyJtJO0faU+Z28rEqrNkCWzZ4pyWlZbdsEdtjBlvjDlijNkUZVkfY8x+Y8wGz6NWHOvWMMZsN8b8Y4zp6svARUQk5Rq2ZhjBB4O5I/sd9H+qf6zyiAjnXtMffwwZMrgQYDISn6HviUANL8uHWGvLeB7zYhYaY/yBEUBN4B6gkTHmnoQEKyIiKd/eU3vpudS56cbIWiNj3XQDYNIkCAyEBg2SOrrk54ZD39baFcaYIrew7YeBf6y1/wIYY6YDdYEtt7AtERFJBay1tJnXhrPhZ2l4b0OeKfFMrDphYfDeezB7duq/13R8GGvtjSs5iXqutbaU53Uf4BXgDLAOeMdaezLGOs8BNay1LTyvmwCPWGvfjGMfrYBWAEFBQWWnT59+S2/Im7CwMAIDY/9iS6vUHteoLaJTe1yjtojOV+2x9MhS+m/tTxb/LEx+eDK50ueKVWfChCLs35+Jnj2T71RvX/99VK5cOdhaW85robX2hg+gCLApyusgwB9n6PwDYLyXdZ4DvozyugkwPD77K1u2rPWlZcuW+XR7KZ3a4xq1RXRqj2vUFtH5oj1OnDth832Sz9IHO2bdGK919u2zNlcua/fuTfDuEpWv/z6AdTaOnHhLp2dZaw9bayOstZHAWJxh7pj2A4WivC7oWSYiImlQ58WdOXL2CBXuqMBrD77mtU7PnvD663DHHUkcXDJ2S6dnGWPyW2sPel7+H7DJS7W1QHFjzJ04CfpF4KVbilJERFK0X/b8wpd/fkl6//SMrj3a6znTwcGwcCHs2OFCgMnYDRO1MWYaUAnIY4wJAXoDlYwxZQAL7AFae+rejjPcXctae9kY8yawEGeYfLy1dnOivAsREUm2Lly+QOu5rQHoXr47d+e9O1adK/ea7tMHsmZN4gCTufjM+m7kZfG4OOoeAGpFeT0PiHXqloiIpB0f/foR249v5648d9G1vPdLakyaBGfPQovYVxFN83RlMhERSTRbjm7ho5UfATC2zlgypIt99ZJjx5x7Tc+fD/7+SR1h8qdrfYuISKKItJG0nNOS8MhwWpdtTfk7ynut9+670KgRPPhgEgeYQqhHLSIiiWJM8Bh+2/cbtwXexoAqA7zW+eUX+Pln55re4p161CIi4nMHQg/Q5ecuAAyrOYwcGXPEqnPxonMq1uefawLZ9ShRi4iIz709/23OXDxDnRJ1aHC39wt2DxwIxYtDvXpJHFwKo6FvERHxqR+2/cC3W78lMH0gI2qNwHi5YPfOnTB0qHPutK7nfX3qUYuIiM+cuXiGtvPaAvDBUx9QKHuhWHWshTZtoFs3KFw4qSNMeZSoRUTEZ3ou7cn+0P08dPtDtH2ordc6U6fC0aPQrl0SB5dCaehbRER84o+QPxi+Zjj+xp+xdcbi7xf7pOiTJ6FTJ/j+e0inDBQv6lGLiEiChUeE03JOSyyWTo93ovRtpb3W69oV6teHRx5J4gBTMP2eERGRBBv02yD+PvI3xXIWo3fF3l7rrFoFc+fqnOmbpR61iIgkyM7jO+n7S18Avqj9BZkCMsWqEx7unDP96aeQPXtSR5iyKVGLiMgts9by+k+vczHiIk1LN6VK0Spe6336KRQsCA0bJnGAqYCGvkVE5JZN3DCRpbuXkjtTbgZXG+y1zu7d8MknsGaNzpm+FepRi4jILdl/Zj8dFnYAYEj1IeTJnCdWnSvnTL/zDhQtmtQRpg5K1CIictOuDHmfvniaZ4o/Q+P7G3utN3Mm/Pefk6jl1mjoW0REbtrXf3/N3B1zyZ4hO6Nrj/Z6mdDTp6FDB/jmG0if3oUgUwn1qEVE5KYcCjvE2/PfBuDT6p9SIFsBr/V69IBataC899tQSzypRy0iIvFmraXNT204eeEk1YtVp3mZ5l7rrVkDs2bpnGlfUI9aRETibcbmGczeNpus6bMyts5Yr0Pely9D69YwaBDkyuVCkKmMErWIiMTLkbNHeHP+mwAMqjbI652xAD7/HHLnhpdfTsroUi8NfYuISLy8Oe9Njp07xtN3Pk3LB1t6rbN9O3z4Ifz2m86Z9hX1qEVE5IZ+OfoLM7fMJEtAFr589ss4h7ybNoXevaFECReCTKWUqEVE5LqOnTvG0J1DAfi4yscUyVHEa72PPoJs2aCt99tQyy3S0LeIiFxXuwXtOBl+koqFK/LGQ294rRMcDMOGwfr14KcuoE+pOUVEJE4/bv+RqX9PJaNfRsY9Ow4/EzttnD8PTZrAkCHOjTfEt5SoRUTEq5PnT/L63NcBeO3O1yiWq5jXet27Q6lS8NJLSRld2qGhbxER8arDwg4cDDvIE4WeoH6B+l7rLFsGM2bAX39plndiUY9aRERimbdzHpM2TiJjuoyMrzve65D36dPwyivw5ZfOedOSOJSoRUQkmtMXTtNqTisA+lfuT4nc3s+1evttqFnTeUjiueHQtzFmPFAbOGKtLeVZ9glQB7gE7AKaW2tPeVl3DxAKRACXrbXlfBe6iIgkhncWvcP+0P08WvBROjzawWud776DVatgw4YkDi4Nik+PeiJQI8ayxUApa+39wA6g23XWr2ytLaMkLSKS/C3atYhxf44jg38Gxj87Hn8//1h1Dh2CNm1g8mQIDHQhyDTmhonaWrsCOBFj2SJr7WXPy98BTcgXEUnhQi+G0nKOc2nQPpX6cHfeu2PVsRZatYJXX4XHH0/qCNMmY629cSVjigBzrwx9xyibA3xjrZ3ipWw3cBKwwGhr7Zjr7KMV0AogKCio7PTp0+P5Fm4sLCyMQP3su0rtcY3aIjq1xzVpsS2G7BjCjwd/pGTWkox4YAT+5lpv+kp7/PTTbXz/fQFGjlxPQMCN80dq5eu/j8qVKwfHOfJsrb3hAygCbPKyvAcwG0/C91JewPNvPmAj8GR89le2bFnrS8uWLfPp9lI6tcc1aovo1B7XpLW2WPLvEksfbEC/APv34b9jlS9btszu2mVtnjzW/h27OM3x9d8HsM7GkRNveda3MeYVnElmL3t24u1HwH7Pv0c8Cf3hW92fiIgkjrBLYbT4sQUA7z35HqXyxRo8JSLCORWra1fn4iaSdG4pURtjagCdgWettefiqJPFGJP1ynOgGrDpVgMVEZHE8e6id9l9ajdlbitD1/JdvdaZObMQxkD79kkcnH6xkdsAACAASURBVMTr9KxpQCUgjzEmBOiNM8s7A7DYc6uz3621rxtjbge+tNbWAoKA2Z7ydMBUa+2CRHkXIiJyS37a8RNfBH9Bev/0TKw7kQD/gFh1/voLpk8vxIYN4B97ErgkshsmamttIy+Lx8VR9wBQy/P8X6B0gqITEZFEc/TsUV778TUA3q/8PqVvi/2VffGic8ON1q3/pUiRu5I6REFXJhMRSZOstbSc05LDZw9TsXBFOj7W0Wu9Pn2gSBGoUeNQksYn1+imHCIiadD4P8fzw/YfyJYhG5PqTfJ6YZNVq2DiRNi4EbZsSfoYxaEetYhIGrPrxC7aLWgHwIhaIyico3CsOidPOkPeX3wB+fIldYQSlRK1iEgacjnyMk1mN+Fs+Fka3tuQl+97OVadyEho1gyefRbq1nUhSIlGQ98iImnIxys/ZnXIagpkLcCoZ0ZhvNxEetAgOHoUZs1yIUCJRYlaRCSNWHdgHX1+6QPAxHoTyZUpV6w6v/wCn34Ka9dC+vRJHKB4paFvEZE04Fz4ORp/15jLkZdp90g7qhStEqvOoUPw0kswaRIUKuRCkOKVErWISBrQeXFnth/fzj157+Gjpz+KVX75MjRqBC1bQvXqLgQocVKiFhFJ5Rb8s4ARa0cQ4BfA1/W/JlNAplh1evWCgAB47z0XApTr0jFqEZFU7Ni5YzT/oTkA/Sv3p8xtZWLVmTsXpkyB4GBdIjQ5UqIWEUmlrLW0ntuaQ2GHqHBHBTo93ilWnd274bXXYPZsyJvXhSDlhjT0LSKSSk3aOInvtn5H1vRZmfx/k2NdfeziRXj+eejWDR5/3KUg5YaUqEVEUqHdJ3fz9vy3ARheazhFchSJVadDB+c63u3aJW1scnM09C0ikspEREbQZHYTQi+F0uDuBjS5v0msOl9/DYsXw7p14OWaJ5KMKFGLiKQyA1cNZNW+VeQPzM/o2qNjXX1syxZo3x5+/hmyZ3cpSIk3DX2LiKQi6w+up9fyXgBMqDuB3JlzRysPC4MGDWDgQCgd+/bTkgwpUYuIpBLnw89fvfrYmw+9SfX/Rb9yibXOBU0efxyaN3cpSLlpGvoWEUkluvzcha3HtnJXnrv4uOrHscpHjXKGvX//3YXg5JYpUYuIpAIL/1nIsDXDSOeXjq/rf03mgMzRyteuhT59YNUqyBT7wmSSjGnoW0Qkhdt/Zj+NZzcGoG+lvjyY/8Fo5SdOQMOG8MUXULy4GxFKQihRi4ikYJcjL/PSdy9x7NwxqhatSpcnukQrj4iAxo2hfn3nISmPErWISArWZ3kfVuxdQf7A/EypPyXW1cfeeQcuXYIBA1wKUBJMx6hFRFKohf8s5MNfP8TP+DG1wVTyZckXrXzkSFi4EFavdu6MJSmTErWISAp05bi0xdKvUj8qFakUrXzhQujXz5k8liOHOzGKbyhRi4ikMJcjL9Po20ZXj0t3r9A9WvnmzdCkCXz3HRQr5lKQ4jM6Ri0iksL0WtaLX//79epxaT9z7av88GGoXRuGDIHy5V0MUnxGiVpEJAVZ8M8CPlr5EX7Gj2kNpkU7Ln3+PNSr5/SmX37ZxSDFp5SoRURSiJAzITSZ7dwJq3/l/lQsUvFqmbXw6qtQuDD07etWhJIYdIxaRCQFiHpculqxanQt3zVaeZ8+sGcPLF2q21amNkrUIiIpwHtL32Plfyu5PevtfPV/X0U7Lj1lCkye7FzDW5cHTX3iNfRtjBlvjDlijNkUZVkuY8xiY8xOz78541i3mafOTmNMM18FLiKSVszfOZ8BqwZ4PS69ciV07Ahz5kBQkItBSqKJ7zHqiUCNGMu6AkustcWBJZ7X0RhjcgG9gUeAh4HecSV0ERGJLepx6fcrv8+ThZ+8WrZrFzz/PHz1FZQq5VaEktjilaittSuAEzEW1wUmeZ5PAup5WbU6sNhae8JaexJYTOyELyIiXlyOvMyLs17k+PnjVC9WnS7lr13H+9Qp5zSs996D6tWvsxFJ8Yy1Nn4VjSkCzLXWlvK8PmWtzeF5boCTV15HWacTkNFa+77n9XvAeWvtIC/bbwW0AggKCio7ffr0W31PsYSFhREYGOiz7aV0ao9r1BbRqT2uSQ5tMebfMUzbN4086fMwtuxYcqR3vmIvXzZ06XI/RYqc5a23/kmSWJJDeyQnvm6PypUrB1try3kr88lkMmutNcbEL+PHvY0xwBiAcuXK2UqVKvkiNACWL1+OL7eX0qk9rlFbRKf2uMbttpi3cx7TfpmGv/Fn9suzKX+Hc/USa6F1a7jtNpgxIyf+/gWTJB632yO5Scr2SMh51IeNMfkBPP8e8VJnP1AoyuuCnmUiIhKHfaf30XR2UwDef+r9q0ka4NNPndnd06eDv39cW5DUJCGJ+kfgyizuZsAPXuosBKoZY3J6JpFV8ywTEREvwiPCefFb57h0jf/VoPMTna+W/fijk6jnzoWsWV0MUpJUfE/PmgasBkoaY0KMMa8BA4CqxpidQBXPa4wx5YwxXwJYa08A/YG1nkc/zzIREfGi59Ke/LbvNwpkLRDtfOlff4XXXoPZs+GOO1wOUpJUvI5RW2sbxVH0tJe664AWUV6PB8bfUnQiImnInO1zGPjbQPyNP9Ofm06ezHkACA6GBg1g2jR4+GGXg5Qkp2t9i4gkA1uObuHl75w7aXzw1AdXj0tv2QLPPANjxkCVKm5GKG5RohYRcdnxc8d5dtqzhF4KpeG9Da8el9692zlH+pNPnLtiSdqkRC0i4qLwiHAazmrIrpO7eDD/g0yoOwFjDAcOOD3obt2c21ZK2qVELSLioo4LO7J091KCsgTx/QvfkzkgM8eOQdWq0KIFtGnjdoTiNiVqERGXjAkew/C1w0nvn57ZL8ymUPZCnDkDNWpAnTpOb1pEiVpExAUr9q6g7by2AIypPYbHCj3GuXNOgn74YfjoI5cDlGRDiVpEJIntObWHBjMacDnyMh0f7UizMs24dAmeew4KFYLhw8EYt6OU5EKJWkQkCYVdCuPZac9y7NwxavyvBgOrDiQiAho3hvTpYcIE8NM3s0Thk5tyiIjIjUXaSJrMbsLfR/6mZO6STGswDT/jT8tWcOKEc2nQgAC3o5TkRolaRCSJ9Fneh++3fU/2DNn5sdGPZM+Qg3fegc2bYfFiyJjR7QglOVKiFhFJAjM2z6D/iv74GT9mPD+DErlL0LcvLFkCy5eDbvUscVGiFhFJZOsPrueV718BYHC1wVQrVo3PPoOvv3ZutpEzp7vxSfKmRC0ikogOhR2i7vS6nL98nlfLvEq7R9oxbhwMGeIk6aAgtyOU5E6JWkQkkVy8fJH639Qn5EwIjxd6nJHPjGTMGEP//rB0qW5XKfGjRC0ikgistbz+0+usDllNoWyF+K7hd4wcloHPP4dffoFixdyOUFIKJWoRkUQw5PchTNwwkcwBmfn+hR8Y+1kQkyfDihXORU1E4kuJWkTEx+bvnM+7i98FYELdicwc9gBz5jhJ+rbbXA5OUhwlahERH9p2bBuNvm1EpI3kvQq9WDXmeVaudE7BypPH7egkJVKiFhHxkZAzIVSfUp3TF09T/64G7J/am61bnHOlc+RwOzpJqZSoRUR84Ni5Y1T7qhr/nf6PRws8Rro5k9h9wI9Fi3QxE0kYJWoRkQQKvRhKra9rsfXYVu7NW4rcC+cSejYLP/0EmTK5HZ2kdLpHi4hIAly8fJH6M+qz9sBaCmcvQt4FC0kfkYvZs5WkxTeUqEVEblFEZASNZzfm539/Jl/mIIIWLeb2rLczYwZkyOB2dJJaKFGLiNwCay1v/PQGs7bMIlv67AQtXsj9Bf7H5MmQTgcVxYeUqEVEbkGPpT0Yu34sGf0zErRkDpXuKs3o0eDv73ZkktooUYuI3KTBvw3mo5Uf4W/8ybt8Jg0eqsDQoeCnb1RJBPqzEhG5CRP+nECnxZ0AyP3rRFpVqs1HH4ExLgcmqZaOpIiIxNP3276nxZwWAAT+OpQPXmxMixYuByWpnhK1iEg8LN+znBdnvUikjSTz2vf4rsvbVK3qdlSSFihRi4jcQPCBYJ6d9iwXIy4SuLUNqz7sy/33ux2VpBW3fIzaGFPSGLMhyuOMMaZ9jDqVjDGno9TplfCQRUQSya5d0KYNZMsGwcGQLRvb32pEjclVCb0USs6QF9k6eBj3368D0pJ0brlHba3dDpQBMMb4A/uB2V6q/mqtrX2r+xERSRLz58Nzz0F4uPMAQkwoVQJmcOxiJHmPVmbboEnkyq45uJK0fPUX9zSwy1q710fbExFJOrt2OUn63LmrSfq0PUvlJgGEZI/ktn1F+OerNeQ6ts/lQCUtMtbahG/EmPHAemvt8BjLKwHfAiHAAaCTtXZzHNtoBbQCCAoKKjt9+vQEx3VFWFgYgbp9zVVqj2vUFtGl2fb47z84dgw834fn7AXeOPcV/9nt5Dp/JxPyNCebXxbImxcKFXI5WHek2b+NOPi6PSpXrhxsrS3nrSzBidoYkx4nCd9rrT0coywbEGmtDTPG1AKGWmuL32ib5cqVs+vWrUtQXFEtX76cSpUq+Wx7KZ3a4xq1RXRptj2yZYPQUABOZYTHXsrFtjtOEHQyM+vHn+P20Cj1Tp92L04Xpdm/jTj4uj2MMXEmal8MfdfE6U0fjllgrT1jrQ3zPJ8HBBhj8vhgnyIivhMWBsDhLFCqWV623XGCnJG5WflVlCQdpZ5IUvJFom4ETPNWYIy5zRjnej3GmIc9+zvug32KiPhOYCB7shnuap6P/fmPcufxdIwKfJ3/nYhdTySpJShRG2OyAFWB76Ise90Y87rn5XPAJmPMRuBz4EXri4PiIiI+tOHl+tzzah5O5TnCvYf8WD3+MkF+OaNXCgiAJk3cCVDStARd8MRaexbIHWPZF1GeDweGx1xPRCS5mB/8N89mX8DlTEd5ZB/M/zqSnBdga8yKAQHQoYMbIUoapxMCRSTN+nTG7zwzqyKXMx3m6Wxl+HlWJnJGBESvFBAAmTPDrFlQrJg7gUqapkQtImlORAQ07b2UTn9VwWY8Sd2SdZn71moCg/+GVq2c2d3g/NuqFfz1F9Ss6W7QkmbpWt8ikqYcPw5V3/qBDf97ARtwkSb3N2F83fGk80vn9JiHD3cey5en2VOxJHlRj1pE0ozgYLir4ddsKNEA63+Rtg+1ZWK9iU6SFkmmlKhFJE0YNw4qdRrF8SebYE0E3ct3Z1jNYfgZfQ1K8qafkSKSql24AG++CT+eGEBYpW4AfFzlYzo/0dnlyETiRz8lRSTV2rMHnihvWZWpK0dLd8Ng+OKZL5SkJUVRohaRVGnBAnjk0UgyNmjDtjwf42/8+br+17Qu19rt0ERuioa+RSRViYyEDz6AUaPDub/PK/x8eCoZ/DMw8/mZ1ClZx+3wRG6aErWIpBonTzpX+Tx+7gTF+7zIz/sXE5g+kB9f/JHKd1Z2OzyRW6KhbxFJFX79FR58EHKW/JsjdR9ixf7F5M2clyVNlyhJS4qmRC0iKdqlS9C9O7zwArzYbxazcz/Gv6f+5cH8DxLcKpiHCzzsdogiCaJELSIp1tat8Nhj8PemSBqO7sGAf5/nbPhZGt/fmJXNV1IoeyG3QxRJMCVqEUlxrIURI+DJJ6Fpy9PYF59l6PoP8TN+fFrtUybXm0ymgExuhyniE5pMJiIpyqFD0Lw5nDgBk+dvo/3vddmxcwe5MuXim+e+oUrRKm6HKOJT6lGLSIoxezaUKQMPPwxdJ8zhhcUPs+P4Du7Ldx9rW65VkpZUST1qEUn2QkOhQwdYtgy+/S6SpeEfUH9mLwCev+d5JtSdQJb0WVyOUiRxqEctIsna6tXwwAPO81/XhPJpyPP0Wt4Lg+Gjpz/im+e+UZKWVE09ahFJlsLD4f33YfRoGDUK7qv4D9Wm12Pz0c1kz5CdaQ2mUbN4TbfDFEl0StQikuzs3AmNG0POnPDnn/DX2YU8NPZFTl04xd157uaHF3+geO7ibocpkiQ09C0iyUZ4OHz8sXNudOPGMG+e5at/B1Jrai1OXTjFsyWf5fcWvytJS5qiHrWIJAtr1kDLlpA/P6xdC1nyHaH+jFb8sP0HAHpX7E2vir3wM+pfSNqiRC0irgoNhffeg2++gcGDoVEj+H7bbFqPbM3Rc0fJliEbk+pNot5d9dwOVcQVStQi4pq5c6FtW3jqKdi0CfyznOKVH9oxeeNkAJ668ykm1J3AHdnvcDlSEfcoUYtIkjt4ENq1cyaKTZjgJOqf//2Z5pObE3ImhIzpMjKwykDaPtxWQ92S5ul/gIgkmchIGDsWSpeG4sXhr7/g0QrneGveW1T9qiohZ0J4uMDDbGi9gbceeUtJWgT1qEUkiWzbBq1aObelXLIE7rsPfg/5nWbfN2PH8R2k80tHn4p96FK+C+n89NUkcoV+ropIorp4Efr2hQoVnHtGr1oFJe+5RM+lPXli/BPsOL6DUvlKsabFGno82UNJWiQG/Y8QkUTz669OL7pECed4dMGC8Pfhv2n6fVM2HNqAwfDu4+/Sr3I/MqbL6Ha4IslSghO1MWYPEApEAJetteVilBtgKFALOAe8Yq1dn9D9ikjytXcvdO7sXKd7yBCoXx8ibQQfrxxEr+W9uBRxiaI5izKx7kQqFK7gdrgiyZqvetSVrbXH4iirCRT3PB4BRnn+FZFUJiwMBgxwrs3drp0zoztzZth1YhfNvm/Gqn2rAGhdtjWDqg0iMH2gyxGLJH9JMfRdF5hsrbXA78aYHMaY/Nbag0mwbxFJApGRMHky9OjhnGq1caMzzB0eEc5nv4+g59KenA0/S/7A/Ix7dpxupiFyE3yRqC2wyBhjgdHW2jExygsA+6K8DvEsU6IWSQV+/dW5V3RAAHz3HTziGS9b+M9C2i9sz7Zj2wBoVKoRw2sNJ1emXC5GK5LyGKejm4ANGFPAWrvfGJMPWAy8Za1dEaV8LjDAWrvS83oJ0MVauy7GdloBrQCCgoLKTp8+PUFxRRUWFkZgoIbYrlB7XKO2iO5m2uPQoYx88UVRtm7NRqtW//LUU0cwBkLOhTDq31H8dvw3AApkKkCbom14PM/jiRm6z+lvIzq1R3S+bo/KlSsHx5zjdZW11mcPoA/QKcay0UCjKK+3A/mvt52yZctaX1q2bJlPt5fSqT2uUVtEF5/2OHPG2u7drc2Vy9q+fa09e9ZZfvrCafvuondtQL8ASx9s1g+z2oErB9oL4RcSN+hEor+N6NQe0fm6PYB1No6cmKChb2NMFsDPWhvqeV4N6Bej2o/Am8aY6TiTyE5bHZ8WSXEiI2HSJOc49NNPXzsOHWkjmfDnJLot6cbhs4cBaF6mOR8+/SG3Bd7mctQiKV9Cj1EHAbOdM7BIB0y11i4wxrwOYK39ApiHc2rWPzinZzVP4D5FJIn98gt07Ajp08Ps2deOQ6/et5q3F7zNugPOkazHCj7G0BpDeajAQy5GK5K6JChRW2v/BUp7Wf5FlOcWaJuQ/YiIO9ascW5BuWMHfPCBcwtKY2D/mf10XdKVKX9NAeD2rLczsMpAXrrvJTw/3EXER3RlMhGJZcMG6NUL1q+Hnj3h1Ved3vSFyxcY/NtgPlz5IefCz5HBPwOdHu9E1/JddU60SCJRohaRq/buzUzDhs4pV127wowZkDGjM+n0u62zeWfRO+w5tQeA+nfXZ1DVQdyZ8053gxZJ5ZSoRYRdu5wbZ8yZU4auXZ0rimXJ4iToBf8spP+K/vy2zzndqlS+UgytMZSn7nzK5ahF0gYlapE07L//oH9/Z4LY22/DlCl/8MwzFbDW8uP2OfRf0f/qRLHcmXLTr3I/WpVtpTtciSQh/W8TSYMOHoQPP4SpU6F1a2eyWK5csHRZOLO2zOL9Fe+z8fBGAPJlycc7j73DG+XeIGuGrC5HLpL2KFGLpCFHj8LAgTBuHDRvDlu3Qr58EBEZwdS/v6H7uu7sPbcXgPyB+en8RGdalW1F5oDMLkcuknYpUYukAXv3wqefwldfOadY/f03FCjg3DRj4oav+fDXD9l5YicAhbIVolv5bjR/oLnuES2SDChRi6Rif/3l9KDnz4cWLWDTJrj9drgUcYkxwRMZsHIAu0/tBqBozqLUz1ufDxp+QHr/9C5HLiJXKFGLpDLWwooV8PHHzvnQ7drBiBGQPbtzHvSINeMYsGoAIWdCACiRuwQ9KvTgpfteYuWKlUrSIsmMErVIKhEZCT/84CToEyfg3Xed205mzAhHzh5hwMrxfP7H5xwMcy61f2/ee+n5ZE+ev+d5/P38XY5eROKiRC2Swl286Bx7/uQTyJEDunSBunXBz8+yat8qRq0bxczNMwmPDAegzG1leO/J96h3Vz38jJ/L0YvIjShRi6RQp0/D6NEwdCjcf7/zvGJFCLsUypj1Uxi5biSbjmwCwM/4UadEHdo81IbqxarretwiKYgStUgKs3s3jBoF48dD9eowbx6ULg1/Hf6LNj+NYsrfUwi7FAY450C3eKAFrcq2onCOwi5HLiK3QolaJAWIjIRFi5xJYatXQ7NmsG4d5C94kW+3fkvb8SNZtW/V1fpPFn6SN8q9Qf2762tymEgKp0QtkoydPAkTJ8LIkRAYCG++Cd98A0cu7WH0utGMmzmOo+eOApA1fVaalm7K6+Vep1S+Uu4GLiI+o0Qtkgxt2OD0nmfNglq1YNIkeOChCyz4Zz7Pf/8l83fOx2IBKB1UmjfKvcHL97+sW02KpEJK1CLJxKVLTmIeMcK5Wcbrr8OmLZfZcm4ZYzdN5bvB33Hm4hkA0vunp+G9DWlTrg2PFnxUk8NEUjElahGXhYQ4M7a//BLuuQc6drTkK/s7M7ZM5cGvZ3Dk7JGrdR+47QFeuu8lmpVuRt4seV2MWkSSihK1iAvCw2HBAue+z8uXw8svwxff/c3vZ6fSafN09mzac7Vu8VzFeem+l2hUqhEl85R0LWYRcYcStUgS2rzZSc5TpkCxYvBM438p1WYa3+2YxvBFm6/WK5C1AC+WepFGpRrxYP4HNbSdwoSHhxMSEsKFCxfcDsVnsmfPztatW90OI9m41fbImDEjBQsWJCAgIN7rKFGLJLKTJ2HaNGf29v798H9N9/Pa2G9ZcngqPfb/AZ6R7VyZcvH8Pc/TqFQjKhSuoKuGpWAhISFkzZqVIkWKpJofWaGhoWTNqvuRX3Er7WGt5fjx44SEhHDnnXfGez0lapFEEBEBixc7yXn+AsvDdddTotUcIiLmMOLQeljv1MsSkIV6d9WjUalGVC1WVec8pxIXLlxIVUlafMMYQ+7cuTl69OhNradELeJDO3Y4yXnS1+fJct9SgirMIXO5Ofx89gDsd+pkSpeJasWq0ahUI2qXqE2W9FlcjVkSh5K0eHMrfxdK1CIJdOyYc5eqsdMOsp2fyPfEHE60XMyBiPPs9ByiLJC1ALVL1KZOiTo8dedTZArI5G7QIkmgffv21K9fnyeffDLOOlWqVGHmzJnkzJkzCSNLWZSoRW7BqVMwe7bly7kbCQ6dQ+YH5nCy0loAQgEioGz+stQpUYc6JevwwG0PqIclacrx48f5/fff+eyzz65br0mTJowcOZIePXokUWQpj2arSNqzaxe0aQPZskFwsPNvmzbO8usIDYUvJh+nXNNZ5HutNW3/Kcxv9z/AxSd6cTLzWjKmy0jtErUZXXs0+zvuZ12rdfSu1FuztsUVvXr1ipYke/TowdChQ2+4XnBwMBUrVqRs2bJUr16dgwcPcvr0aR588EG2b98OQKNGjRg7diwAgYGBdOjQgXvvvZenn3766vHXb7/9lho1agBw+vRpSpYs6XX9Z599lmnTpvnujadC6lFL2jJ/Pjz3nHMic7hzf2ZCQ52rjUya5FwarGbNq9VPnLnAZ9+uYsa6n9kZuZjIoPVQzLl0ZziQPzD/1SHtp4s+TeaAzC68KZHYXn31VerXr0/79u2JjIxk+vTpLF26lDJlynitP3XqVIoXL85bb73FDz/8QN68efnmm2/o0aMH48ePZ9CgQbzyyiu0a9eOkydP0rJlSwDOnj1LuXLlGDJkCP369aNv374MHz6cVatW8dxzzwHOqUzDhw/3un7OnDm5ePEix48fJ3fu3EnTOCmMErWkHbt2OUn63LnYZZ7EHflcA9Ys+oZRm7ex8J/FHM7wK6S7APmcaun90/NEoSeoWrQqVYtV5cH8D+o0KrmhxBhQsfb65UWKFCF37tz8+eefHD58mAceeIDChQuzYcOGONfZtGkTmzZtomrVqgBERESQP39+AJ566il++ukn2rZty8aNG6+u4+fnxwsvvABA48aNqV+/PgAHDx4kb95rV8+rWrUqM2fOjLU+QL58+Thw4IASdRyUqCXtGDz4Wi86ir3ZYW7RDEwtlpfgO09z8ednnQLPZOzSQaWpWrQqVYpWoULhCuo1y027UVJNLC1atGDixIkcOnSIV199ldDQUCpUqOC17tSpU7HWcu+997J69epY5ZGRkWzdupXMmTNz8uRJChYs6HU7Vw7zZMqUKdoFX663/oULF8iUSRMs46JELWnHlCkQHs7+rLCiMMwrkpm5x4dwqgPARSAEgNvPpKPGk035//buParKKn/8+HsDRxDBK4mNmpiKcpNQxAteMB2n0kkovqlLa9Rpyubr9KU1rXS+a8bSn/qrlqWV9asm+2mWaF7za1rjr5HUJBEQFBQUFU25oyCkct2/Px5ADhc9yuVw+bzWOotzzrOf5+zzccuH/Tz72Xvyw5OZ9PAkenbqac1aC3HfQkJCWLJkCSUlJWzatAlbW9s79qiLi4vJzs4mMjKS0aNHU1JSwpkzZ/Dy8uLDDz/Ew8ODlStXMm/ePCIjIzGZTJSXl7Nt2zZmzpzJpk2bGDt2LAAeHh6kpKQQFBQEwOrVq+vcX2tNRkYGbm5uzRCR1um+E7VSqi/wBeAKaOBTrfV7NcoEAd8ARW/v1QAAGSlJREFUFyre2qG1Xna/nynE/biYd5EfL/5IxMQC/p9bB37pXlyx5QZwg4637Jh0AR4/X8rk8zDoWhnqnXXWrLIQjaJDhw5MnDiRrl27Ymtra1H5bdu28fLLL5Ofn09paSlhYWHY2dmxYcMGoqOjcXZ2Zvz48SxfvpylS5fSqVMnoqKiWL58OT179mTLli0ATJ06lU8++YTnn3+e5ORkPvvsM6KiomrtHxMTw6hRo7Czk35jfRoSmVLgr1rrWKWUMxCjlNqvtT5Vo9whrfW0BnyOEBbTWnP+2nl+vPij8Uj9kYv5F42NwwCKMRXZ43epE09ezKPb7/+TF5Z9gF15tYN0lmkSRdtQXl7Ozz//zNatWy3e55FHHuHgwYO13q9M0gDvvvuu2baarwHGjRvH3/72N/Ly8hg8eLDZvNjVy2/cuJE///nPFtevPbrvRK21TgfSK54XKKVOA72BmolaiCajteZM7hmzxHyl4IpZGduSrgzqMI6p6SU88+1+hl0pwq68CICI4H7mSdpkgmefbcZvIETTOHXqFNOmTSMkJIRBgwZZpQ7vvPMOly5domvXrvWW8fb2ZtKkSc1Yq9anUc41KKXcAD/gaB2bRyul4oE04FWtdWIdZYSwSFFpETHpMRz55Qg//fITR345YrZeM4C62QPXW+MZ99AEnp88gUnePtja2Bqjvr8YCuV1jPquZDLBK6808bcQoul5enpy/vz5Jv+cwsLCereNHDnyrvtX3qYl6qd0A4cjKqWcgB+BFVrrHTW2dQbKtdaFSqkngPe01nX+aaeUegF4AcDV1XX45s2bG1Sv6goLC3Fycmq047V2rSke14qvkXg9kYT8BBKuJ3Cm4AwlusbI7V974pQTyBAHPx4d4MGjPj2w71DP/TDXrxsJW2vQmsI+fXC6fNm4f0YpY+3Jzp2b/ou1UK2pbTS1hsSiS5cuDBw4sJFrZF1lZWUWXeduLxoSj5SUFPLz883emzhxYozW2r+u8g1K1EopE7AH+F5rXfsiRe3yqYC/1jrnTuX8/f11dHT0fderpoiIiKqRh6LlxqNcl5OUk8RPl36q6i2fvXq2VrkOeV7YXAnEz2UMT48IZPbjA+jV6x5uVD13Dlavho0biViyhKBly4zT3a+8YiTqdqyltg1raEgsTp8+jYeHR+NWyMpkmUtzDYlHXe1DKVVvom7IqG8FrANO15eklVK9gEyttVZKBWBMWZp7v58p2pa8W3kcu3KMo1eOEnk5kshfIrl265pZGROOdLo2kl+TxuDhFEiw/yiCn+uGry/Y3O88IwMGwNq1xiMiAmr8ZSuEEC1JQ65RBwLPAieVUpU35v038BCA1vpjIBR4SSlVCtwEZuqGnmsXrVJxWTEnMk8QdSWKo1eOcvTyUZJzk2uV60xvOuYEkncikIfUGKb4+vLYb00ELQc5IyuEaI8aMur7MHDH841a67XA2vv9DNE6aa25kHeBo5ePcvTKUaKuRBGbHktRWZFZuQ429rhqP9SVkWTGjmSAKZApIx8iaDKMWwrdu1vpCwjRBmRkZBAWFsaxY8fo2rUrrq6urFmzBnd393r3cXNzIzo6GhcXF7P3P/74Y5RSvPjiixZ99ooVK6omWLGxseGTTz6548CyN954AycnJ1599VWWLFnC+PHjmTx5smVf9C7WrFnD4sWLyczMpEuXLo1yzOYmd5iLBsv+NZuY9Jiq3nLUlShybtQehvCQ42C63wzgVspILh4ZiXuXoTw6oQMTnoBx/1sSsxCNRWtNSEgIf/jDH6gcmBsfH09mZuYdE3V9FixYQEFBgUVlIyMj2bNnD7Gxsdjb25OTk0NxcfHdd6ywbFnjzokVHh7OiBEj2LFjB/PmzWuUY5aVlTXKcSwlqwmIe5JzI4fvU75nxcEVPLXlKfqt6UfPVT15/KvHeT3idfae3UvOjRx6ODyAf+dpjC3+X3gf/x7H96/SIzyJiXlf8Fbof3LlmD9xMR14912YPl2StBCN6cCBA5hMJhYsWFD1nq+vL+PGjSMiIoJp027PQbVw4ULWr19f9frtt9/Gx8eHgIAAUlJSAKPH+/777wPGiOXJkyfj6+vLsGHDOFdjedj09HRcXFywt7cHwMXFhd/85jeA0WN/7bXXah2/urlz57Jt27aq8q+//jrDhg3Dx8eHpKQkwFixa/78+QQEBODn58c333xTZxzOnTtHYWEhy5cvN1tKc/369UyfPp2goCAGDRrE0qVLAUhNTWXIkCHMnj0bDw8PQkNDuVGxiI+bmxuLFi1i2LBhbN26la1bt+Lj44O3tzeLFi0CYOfOnUyaNAmtNenp6bi7u5ORkXHHfytLSI9a1Cv3Ri4x6TFEp0UTkx5DTFrM7Vm+qulk6sTgzsPodtOf4vMjuRgZQP5FN3qMVIwZA2MWQsAX7fquJ9GOqaVNsxa5fr3+4T4JCQkMHz78vo7bpUsXTp48yRdffEFYWBh79uwx2z579mwWL15MSEgIt27dory83Gz7lClTWLZsGe7u7kyePJkZM2YwYcIEi49fk4uLC7GxsXz00UesWrWKzz77jBUrVvDoo4/y+eefk5eXR0BAAJMnT6ZTp05m+27evJmZM2cybtw4kpOTyczMxNXVFYCoqCgSEhJwdHRkxIgRTJ06FRcXF5KTk1m3bh2BgYHMnz+fjz76iFdffRWAHj16EBsbS1paGiNHjiQ2NpZu3boxZcoUdu3aRUhICNu3b+fDDz/ku+++Y+nSpfTq1eu+/h2qk0QtAOP09fGM48SkxRhJOT2G1LzUWuUcTY74uPjRS/tD2nAyY/1JPOTOr71s8R0Do0cbidnDowGjsoUQVjNr1qyqn6/UmPynoKCAK1euEBISAoCDg0Ot/Z2cnIiJieHQoUMcOHCAGTNm8OabbzJ37ty7Hr8ulctmDh8+nB07jKk6/vWvf7F7925WrVoFGKtvXbp0qdYtT+Hh4ezcuRMbGxuefvpptm7dysKFCwFj2c3KZTWfeuopDh8+THBwMH379iUwMBAwlu18//33qxJ15XKex44dY+zYsVXLeM6ePZuDBw8SHBzMBx98gLe3N6NGjar6rg0libqd0VqTdjON7ae2czzjOHEZcRzPOE5aQVqtsh3tOuLb04++tv50yB3O9aThnP15CImXbHEcYSTlPy2AUetBlpEVom536vk2FS8vr6rTxzXZ2dmZ9YKrL0UJt5eprPn8Xtja2hIUFERQUBA+Pj5s2LChKlHf6/ErT6Hb2tpSWloKGL/Htm/fzuDBg+vd7+TJk5w9e7Zqbe3i4mL69+9flahrfnbl6/reB2r12Oty+fJlbGxsyMzMpLy8HJtG6LFIn6cNKykrIT4jng1xGwj7LowJ6yfQ9a2uzI6aTejWUFYcWsG3Z78lrSANpw5OjO4dyH/0W8izzusJTjvJgK+vE/9fP3Hh/7yH8/nnmD7Gi/CvbLl2Df79b1ixAqZOlSQtREvz6KOPUlRUxKefflr13okTJzh06BD9+vXj1KlTFBUVkZeXxw8//GC2b+XqV1u2bGH06NFm25ydnenTpw+7du0CoKioqOoabqXk5GTOnr09UVFcXBz9+vWz6PiW+t3vfscHH3xA5d2+x48fr1UmPDycN954g9TUVFJTU0lLSyMtLY2LF43Ld/v37+fq1avcvHmTXbt2VfWiL126VLUed/VlO6sLCAjgp59+Iicnh7KyMsLDw5kwYQKlpaXMnz+f8PBwPDw86lys5H5Ij7qNyLuVx4nME8RnxBOXEUdcZhwJWQkUl9UebdnN1A3/PgH0tvXDlONHXvIjnP15IHFJNtwcDP7+EOgP/i+AtzdU/EErhGgllFLs3LmTsLAw3nrrLRwcHHBzc2PNmjX07duXZ555Bm9vb/r374+fn5/ZvteuXWPo0KHY29ubDcCqtHHjRl588UWWLFmCyWRi69atPPzww1XbCwsL+ctf/kJeXh52dnYMHDjQ7A+Gux3fEv/4xz8ICwtj6NChlJeX079//1rXujdv3szevXvN3gsJCWHz5s24uroSEBDA008/zeXLl5kzZw7+/v6kpqYyePBgPvzwQ+bPn4+npycvvfRSrc9/8MEHeeONN5g4cSJaa6ZOncr06dNZtmwZ48aNY+zYsfj6+lZd+27oLHUNnuu7KcgUovUrKy8j5WqKkZQz44nPjOdE5gku5V+qs/yAbgPw7OaHS6kfNlmPkJfkx/GDzqSlOTFwoJGUKx9Dh0LHjs38haysLbWNxiDxuE2mEDXXGFOI1nefdnNbv3490dHRrF1rPs1Hamoq06ZNIyEh4a7HaBVTiIqmd+3mNU5mnSQ+I74qMSdkJXCz9Gatsg52Dng94M1D9kNx+tWX0l/8yDrhS2JsZw7dAl9fIxGPHgu/HRXDc88Nb3dJWQghWiNJ1C1AcVkxyTnJnMw6SUJWAiezTt6xl9y3c1/cuwzlgTJfTNeGUpDiS2rMQE6dtqOgr5GQfX1h1gLjZ9++xsJQlSIiCiRJCyGaXWpqqrWrABj3alcObqvOzc3Not50c5NE3Yy01lzMv2gk48yTnMwyHsk5yZSUl9Qq72DngHsXb3rZ+NIxfyjFl3zJSvDhTHx3EjuBp6fxGOsPvn80ridbMChRCCFEKyKJuolcvXn1djKu+JmQlUBBce1p+BSKh5wH8BtbHzoV+lB6xYfc0z5ciBlEjrMtPT3BzQs8/cBztpGcZSYvIYRoHyRRN1BBUQGJ2YkkZCWQmJVIQnYCCVkJZBTWPW1cN1NPetn44FjoQ1maD3lnfEiL90R360SXih6y5zDwetaYNKRr12b+QkIIIVqUtpuoz52Dd96BL7+E11+HJ5+EOXPgr3811iO+RzdLbpKUk0RClpGIE7KNxFzXlJoA9qoTLtoLh3wfii/7kHvKB4fr3gzp2xN3d4zHVHB/BQYObH+jrYUQQlimbSbqffsgNBRKSowHQEEBfPYZbNgA27bB44/XuWvlwK7E7EQSsxKresvnrp2jXJfXKm+r7XEuGoJtrjc3L3pTlu7Nw07eePd9iMHuNrj7gfsMGDRITlcLIepQvVNRWGgsvN6ATkUlJycnCgsLAdi7dy9hYWHs37/fbPIRS6SmpnLkyBF+//vf33ddqnviiSfYtGkTXe9wunD9+vVMmTKlajGP9q7tJepz54wkXWO2HOB24g4NpSQuhpSu2jhlnZ1YlZjPXj1LaXlprV2VtqXDdQ9K07zp9KsXfTp4M7i7N0P7DMDd246HHzb+Tz3wgPkIayGEqFcDOhWW+uGHH3j55Zf5/vvv7zlJg5GoN23a1GiJuuYkJHVZv3493t7ekqgrtL1E/c47txs8UKbgl/Jsdg6BhJ6Q2BPiexaT8qU3pTZ1rCmqFSpvAE43vHG18WKAsxe+D3ozauBghgyyx81NTlMLIRqBhZ0KTpy47571wYMH+dOf/sTevXsZUHGMuXPn0rlzZ6Kjo8nIyODtt98mNDQUrTWvvfYa+/btQynF3//+d2bMmMHixYs5ffo0gYGBzJs3z2whjYiICJYsWYKzszMpKSlMnDiRjz76CBsbG8LDw1m5cmXVzF1vvfUWcHvSk8LCQh5//HHGjh3LkSNH6N27N9988w3ffvst0dHRzJ49m44dOxIZGcnSpUvZvXs3dnZ2TJkypWoxjvai7SXqL7+sStTHe8HI520oufEWzKxeyOgxOxa74aq8jFPVPb0I6O/FeM8h9O7pKL1iIUTTqtGpqFNJCaxeDTVm0LJEUVERwcHBREREMGTIELNt6enpHD58mKSkJJ588klCQ0PZsWMHcXFxxMfHk5OTw4gRIxg/fjxvvvkmq1atIjw8vM6ZuKKiojh16hT9+vXjscceY8eOHYwZM4ZFixYRExNjtgxkcHCw2b5nz54lPDycf/7znzzzzDNs376dOXPmsHbtWlatWoW/vz+5ubns3LmTpKQklFLk5eXdcyxau7aXqCuuyQD0z4MSu3Kcix9g8EWNd3Y5o7IKGZZVjEeuwunWBStWVAjRrlXrVNSrpAQ2bryvRG0ymRgzZgzr1q3jvffeM9sWHByMjY0Nnp6eZGZmAnD48GFmzZqFra0trq6uTJgwgWPHjtH5LgvJBwQEVM31PWvWLA4fPozJZCIoKKjOZSCr69+/P4888ghgLGNZ14QoXbp0wcHBgT/+8Y9MmzaNadOm3XMsWru2t3qWk1PV06634PpK2N19Ece+yuH//usqL8YVMyINnOwbNmetEEI0SLVORaOUq8HGxoavv/6aqKgoVq5cabbNvtpKOw1d7+FOy0LeTfV6VF/Gsjo7OzuioqIIDQ1lz549PPbYY/df2Vaq7SXqOXPAZKp66Vx78Shj+7PPNl+dhBCipmqdikYpVwdHR0e+/fZbvvrqK9atW3fHsuPGjWPLli2UlZWRnZ3NwYMHCQgIwNnZmYKC2hM1VYqKiuLChQuUl5ezZcsWxo4dS0BAAD/++GOtZSAtVf0zCwsLyc/P54knnmD16tXEx8dbfJy2ou2d+v7rX43Rknc6pWQyQbUBEUII0ezmzDFGd9/td1UDOxXdu3fnu+++Y/z48VWnousSEhJCZGQkvr6+KKV4++236dWrFz169MDW1pYxY8Ywf/58s8FkACNGjGDhwoVVg8lCQkKwsbHhzTffrLUMpKXmzp3LggUL6NixI/v27WP69OncunULrXWjrfHcqmitW9xj+PDhukH27tXa0VFrk0lr0AdWrdIajNeOjsb2duzAgQPWrkKLIbEwJ/G4rSGxOHXq1N0LpaQYv4+g/oejo1GuBbh+/Xqt9w4cOKCnTp1qhdpYX13xsFRd7QOI1vXkxLZ36huM+w5PnIAXXoDKgRCdOxuvT5xo8H2JQgjRYAMGGPdJOzqaXa4DjNeOjsb2Bkx6ItqGtpmowWjca9dCfj4MH278XLtWGr0QouWo2amwsWlVnYqgoCD27Nlj7Wq0eW3vGrUQQrQmlZ2K+7gFS7QPbbdHLYQQVqQbeNuTaJvup11IohZCiEbm4OBAbm6uJGthRmtNbm4uDg4O97SfnPoWQohG1qdPHy5fvkx2dra1q9Jobt26dc8Jpi2733g4ODjQp0+fe9pHErUQQjQyk8lE//79rV2NRhUREYGfn5+1q9FiNGc85NS3EEII0YJJohZCCCFaMEnUQgghRAumWuKoRKVUNnCxEQ/pAuQ04vFaO4nHbRILcxKP2yQW5iQe5ho7Hv201nVOxt4iE3VjU0pFa639rV2PlkLicZvEwpzE4zaJhTmJh7nmjIec+hZCCCFaMEnUQgghRAvWXhL1p9auQAsj8bhNYmFO4nGbxMKcxMNcs8WjXVyjFkIIIVqr9tKjFkIIIVqlNpWolVKPKaWSlVIpSqnFdWy3V0ptqdh+VCnl1vy1bB4WxGKuUipbKRVX8XjeGvVsDkqpz5VSWUqphHq2K6XU+xWxOqGUGtbcdWxOFsQjSCmVX61tLGnuOjYXpVRfpdQBpdQppVSiUuq/6ijTbtqHhfFoT+3DQSkVpZSKr4jH0jrKNH1e0Vq3iQdgC5wDHgY6APGAZ40yfwY+rng+E9hi7XpbMRZzgbXWrmszxWM8MAxIqGf7E8A+QAGjgKPWrrOV4xEE7LF2PZspFg8CwyqeOwNn6vi/0m7ah4XxaE/tQwFOFc9NwFFgVI0yTZ5X2lKPOgBI0Vqf11oXA5uB6TXKTAc2VDzfBkxSSqlmrGNzsSQW7YbW+iBw9Q5FpgNfaMPPQFel1IPNU7vmZ0E82g2tdbrWOrbieQFwGuhdo1i7aR8WxqPdqPg3L6x4aap41BzY1eR5pS0l6t7AL9VeX6Z2A6sqo7UuBfKBHs1Su+ZlSSwAnq44lbdNKdW3earWIlkar/ZkdMXpvn1KKS9rV6Y5VJyy9MPoNVXXLtvHHeIB7ah9KKVslVJxQBawX2tdb/toqrzSlhK1uDf/A7hprYcC+7n9F6EQsRjTGfoCHwC7rFyfJqeUcgK2A2Fa6+vWro+13SUe7ap9aK3LtNaPAH2AAKWUd3PXoS0l6itA9V5hn4r36iyjlLIDugC5zVK75nXXWGitc7XWRRUvPwOGN1PdWiJL2k67obW+Xnm6T2u9FzAppVysXK0mo5QyYSSlr7TWO+oo0q7ax93i0d7aRyWtdR5wAHisxqYmzyttKVEfAwYppforpTpgXNTfXaPMbuAPFc9DgX/rihEAbcxdY1HjGtuTGNei2qvdwHMVo3tHAfla63RrV8palFK9Kq+xKaUCMH5PtMU/aKn4nuuA01rrd+sp1m7ahyXxaGft4wGlVNeK5x2B3wJJNYo1eV6xa8yDWZPWulQptRD4HmPU8+da60Sl1DIgWmu9G6MBblRKpWAMpplpvRo3HQtj8bJS6kmgFCMWc61W4SamlArHGKnqopS6DLyOMSgErfXHwF6Mkb0pwA1gnnVq2jwsiEco8JJSqhS4Ccxso3/QAgQCzwInK65DAvw38BC0y/ZhSTzaU/t4ENiglLLF+IPka631nubOKzIzmRBCCNGCtaVT30IIIUSbI4laCCGEaMEkUQshhBAtmCRqIYQQogWTRC2EEEK0YJKohRBCiBZMErUQQgjRgkmiFkIIIVqw/w+ep4SRcWKGTAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 576x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}