{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import *\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cubic splines with Python\n",
    "\n",
    "In this notebook, we show how to use the `interpolation` library to interpolate a function using cubic splines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# we need to import the library first\n",
    "from interpolation import splines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolating a function\n",
    "\n",
    "Our goal here, consists in approximate a function $R^d \\rightarrow R$ given its values on \n",
    "a regularly spaced rectangular grid.\n",
    "\n",
    "We consider the two dimensional function defined by $f(x,y)=sin(\\sqrt{x^2+y^2}/\\pi)/(\\sqrt{x^2+y^2}/\\pi)$ if $(x,y) \\neq (0,0)$ and by $f(0,0)=1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numba import vectorize\n",
    "@vectorize(['double(double, double)'])\n",
    "def f(x,y):\n",
    "    if x==0 and y==0:\n",
    "        out = 1.0\n",
    "    else:\n",
    "        t = sqrt( x**2+y**2 )\n",
    "        out = sin(pi*t)/(pi*t)\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Value of f at (0.5,0.5) : 0.358187786013244\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7faff97d9f98>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcjeX/x/HXZ4yxk72ykyVbWSPFse+GhLEkW5aM8i2h\n1UgqWaJFRETRRDLWsXfsW4QwgyiyJGtlyyzX7w/z7eerwSxnznWWz/PxmEdzjmvO/Z7TmfM+93Vv\nYoxBKaWUfwqwHUAppZQ9WgJKKeXHtASUUsqPaQkopZQf0xJQSik/piWglFJ+zCUlICKfichpEdlz\nhzEfiMghEdklIg+7YrlKKaVSx1VrAtOBxrf7RxFpCpQwxpQE+gCTXLRcpZRSqeCSEjDGbAAu3GFI\nMDAzYexWIIeI5HfFspVSSqWcu7YJFAB+ven2iYT7lFJKWaQbhpVSyo8Fumk5J4BCN90umHDfv4iI\nnsxIKaWSyRgjKfk5V64JSMJXYhYCXQFEpAZw0Rhz+nYPZIzRLxd8DRs2zHoGT/26HnudqTumUnxC\ncQq/X5in5z/NjF0zOHbxWLKfzyPnj/DZzs/o8m0XCowtQOkPS/PF7i+IjYu1/nt68pe+Pl33lRou\nWRMQkdmAA8gtIseAYUAQYIwxnxpjlopIMxH5CbgMdHfFcpVKrtj4WL7Y/QVvrX+LovcUZUbrGdQq\nVAuRFH2IAqBYzmIUy1mMHpV6YIxhzc9rGOYcxlvr3uKNOm/QoVwH0gWkc+FvoZTruKQEjDGdkjAm\n1BXLUiolYuNjmbVnFiPWjaBQjkJMD55O7SK1Xb4cEaF+8frUK1aPVUdW/U8ZtC/XngDRzXDKs+gr\n0oc5HA7bETxCRHQEZT8uy7Rd05jaairfPf1digogOc+niNCwREM29tjI+43fZ/yW8VT4pALLflqW\n7OX6Kn19egZJ7XySq4mI8bRMyjtdi73GC8tfYPnh5UxuMZn6xeqnatonNYwxLD20lH5L+tGhXAfe\nrv826dOlt5JF+R4RwaRww7CWgPJJB88dpP3c9pTMXZKpLaeSI2MO25EAOHvlLN0iunH2ylnCnwyn\n6D1FbUdSPiA1JaDTQcrnzNozi1rTatGnSh/mPDnHYwoAIE/mPCzquIj25dpTfUp15kfNtx1J+Tld\nE1A+40rMFZ6LfI51R9cxp90cHr7Xs89TuPX4VkLmhdCyVEtGNxxNhsAMtiMpL6VrAsrvHTh7gOpT\nqnM19io7eu/w+AIAeKTgI+zsvZPjfx7n0WmP8svFX2xHUn5IS0B5vT2n9+CY4SC0eihftvmSbBmy\n2Y6UZDkz5WRe+3l0rtCZ2tNrc/DcQduRlJ9x12kjlEoTO0/tpNmsZoxvMp6Q8iG246SIiPBCzRfI\nkSEHdWfUZeVTKymbt6ztWMpPaAkor7X1+FZaftWSyS0m0+bBNrbjpFrPyj3JEJiB+jPrs6zzMh66\n9yHbkZQf0BJQXmn90fW0ndOW6cHTaV6que04LtOlYheC0gXR6MtGLOm0hKr3V7UdSfk4LQHlddb8\nvIYO33Rg1hOzaFSike04Lte+XHuC0gXRbFYzFoQsoGahmrYjKR+mJaC8yrKflvHU/KeY224ujqIO\n23HSTOsyrQlKF0Sr8FbMaz8vTc5zpBTocQLKi6w6sopO8zoRERLBo4UetR3HLVYdWUXHeR1ZELLA\nb35nlXx62gjl8/af2Y/jcwdz282lTtE6tuO4VeShSHos7MHGHhspnrO47TjKA+nBYsqn/X75d1rM\nbsHohqP9rgAAmpZsyquPv0qL2S24eO2i7TjKx+iagPJo12KvUW9GPeoVq8db9d6yHceq5yKfI+ps\nFEs7LdUzkKr/odNByicZY+j8bWfiTBxftf3K7y/IEhsfS3B4MIWyF+KT5p9YOy228jw6HaR80vC1\nwzly4QifB3/u9wUAEBgQSHjbcDb9uon3t7xvO47yEbqLqPJIs/bM4vNdn7O111Yypc9kO47HyJYh\nG4s7LabmZzUpkbMEwWWCbUdSXk4/XimPs+HYBv6z/D8s7rSY/Fnz247jcQrnKExEhwh6LerFzlM7\nbcdRXk5LQHmUY38co93cdnzR5gvK5ytvO47HqlagGpNbTCY4PJjTl07bjqO8mG4YVh4jNj4Wx+cO\nWpZqyZDHhtiO4xVeW/Ma209uJ7JzpG438WO6YVj5hDBnGFmCsvBSrZdsR/EaYY4wLl+/zJhNY2xH\nUV7KJSUgIk1EJFpEDorIvz7CiUh2EVkoIrtE5EcR6eaK5SrfsfrIaqb9MI2ZrWfqJ9pkCAwIZHbb\n2YzdPJYtx7fYjqO8UKr/2kQkAPgIaAyUAzqKSJlbhvUH9hljHgbqAmNFRPdMUgCcvnSarhFdmdlm\npm4IToHCOQozucVkOs7rqEcUq2RzxUeu6sAhY8xRY0wMEA7cut+aAf57zb9swDljTKwLlq28XLyJ\n5+mIp+n2UDcaFG9gO47Xal2mNS1KtqDXwl7oNjWVHK4ogQLArzfdPp5w380+AsqKyElgN/C8C5ar\nfMCYTWP46/pfhDnCbEfxeqMbjebwhcNM3jHZdhTlRdw1JdMY+MEYU09ESgArRaSiMeZSYoPDwsL+\n+d7hcOBwONwSUrnXluNbGLt5LNt6bdNz4bhAxsCMhLcN57Hpj/FooUepmL+i7UgqjTidTpxOp0se\nK9W7iIpIDSDMGNMk4fZQwBhjRt00ZjHwjjFmY8Lt1cAQY8z3iTye7iLqBy5eu0ilyZUY12icT1wf\n2JPM3D2Tdze8y/ZntpMlKIvtOMoNbO8iuh14QESKiEgQEAIsvGXMUaABgIjkB0oBR1ywbOWFjDH0\nXtSb5iWbawGkga4PdaXq/VV5fpnOuqq7S3UJGGPigFBgBbAPCDfGRIlIHxHpnTDsLeBREdkDrAQG\nG2POp3bZyjt9ve9r9p3Zx5hGum97Wvm42ces+XkNSw4usR1FeTg9Yli51e+Xf6fiJxVZ2HEh1QtU\ntx3Hp635eQ1PRzzN3n57yZExh+04Kg3p9QSU12g/tz1F7ynKew3fsx3FL/Rd3JfY+FimtppqO4pK\nQ7a3CSiVJPP2z2P36d0Mdwy3HcVvvNfwPVYeWcmKwytsR1EeSktAucW5K+cIjQxlevB0vT6AG2XP\nkJ1PW3xK70W9+evvv2zHUR5Ip4OUW3T5tgt5M+fl/SZ6RSwbei7oSYbADExsPtF2FJUGdDpIebRF\nBxax5fgWRtYfaTuK3xrbeCwLDyzku5+/sx1FeRgtAZWmLly9QL8l/ZjaaiqZ02e2Hcdv3ZPxHia1\nmESvRb24fP2y7TjKg+h0kEpTPRb0IFNgJj5u/rHtKAp4av5T5M6Um/FNxtuOolxIp4OUR1r20zK+\n++U7RjUcdffByi0mNJnAnH1z2Hhso+0oykNoCag0cfn6Zfou7sunLT4la1BW23FUglyZcvFRs4/o\nubAn1+Ou246jPICWgEoTI9ePpFbhWjQs0dB2FHWLNmXaUDJ3ScZuGms7ivIAuk1AuVzUmShqf16b\nPX33cF+2+2zHUYn4+cLPVJtSjR29d1DkniK246hU0m0CymMYY+i/tD+v135dC8CDFctZjIE1BuqZ\nRpWWgHKt8L3hXLh2gWerPWs7irqLlx59if1n9rP44GLbUZRFWgLKZf649geDVg5iYrOJBAa466J1\nKqUyBGbgo2Yf8Vzkc1yJuWI7jrJES0C5zDDnMJo90IyahWrajqKSqFGJRlQrUI131r9jO4qyRDcM\nK5fY9dsuGn/ZmH3P7iNP5jy246hkOPHnCR6a9BCbem6iVO5StuOoFNANw8qqeBNPvyX9GFlvpBaA\nFyqQvQCvPP4K/Zf2Rz+A+R8tAZVq03+YDkCPSj0sJ1EpNaD6AH679Btz98+1HUW5mU4HqVQ5d+Uc\nZSeWZVnnZVS6r5LtOCoVNhzbQMg3IUT1jyJbhmy246hk0MtLKmv6LOpDhsAMfND0A9tRlAt0i+hG\nnsx5GNNojO0oKhm0BJQV/90YHN0/mpyZctqOo1zgt0u/UW5iOTb33Kwbib2IbhhWbmeMYeCygYTV\nCdMC8CH3Zr2XwY8OZtCKQbajKDfRElApMj96PueunuOZKs/YjqJcbGCNgew7s4+Vh1fajqLcwCUl\nICJNRCRaRA6KyJDbjHGIyA8isldE9Bp3Xuxa7DUGrRjE+Mbj9chgH5QhMANjGo7hP8v/Q2x8rO04\nKo2lugREJAD4CGgMlAM6ikiZW8bkAD4GWhhjygPtUrtcZc/4LeOpmL8i9YvXtx1FpZHWZVqTL0s+\nJn8/2XYUlcZcsSZQHThkjDlqjIkBwoHgW8Z0AuYZY04AGGPOumC5yoJTf51izKYxuveIjxMRxjcZ\nz/C1wzl/9bztOCoNuaIECgC/3nT7eMJ9NysF5BKR70Rku4g85YLlKgteXfMqPSr14IFcD9iOotJY\nxfwVeeLBJ3hz7Zu2o6g05K4J3UCgMlAPyAJsFpHNxpifEhscFhb2z/cOhwOHw+GGiOpudpzcQeRP\nkRwIPWA7inKTEXVHUHZiWfpW7UuZPGXu/gPKLZxOJ06n0yWPlerjBESkBhBmjGmScHsoYIwxo24a\nMwTIaIwZnnB7KhBpjJmXyOPpcQIeyBjD49Mfp9vD3ehVuZftOMqNxm0ex6ojq1jaeantKOo2bB8n\nsB14QESKiEgQEAIsvGXMAuAxEUknIpmBR4AoFyxbucmcfXO4HHOZ7g93tx1FuVlo9VAOXzhM5KFI\n21FUGkh1CRhj4oBQYAWwDwg3xkSJSB8R6Z0wJhpYDuwBtgCfGmP2p3bZyj2uxlxl8KrBjG88nnQB\n6WzHUW4WlC6IsY3G8sKKF4iJi7EdR7mYnjZC3dXIdSP54bcf+Kb9N7ajKEuMMTT+sjEtSrXguUee\nsx1H3ULPHaTSzOlLpyk3sRxbe22lRK4StuMoi348/SMNvmjAgdAD3JPxHttx1E1sbxNQPmz42uE8\nVfEpLQBFhfwVaFGyBe9ueNd2FOVCuiagbivqTBS1P69NdP9ocmfObTuO8gAn/jxBxUkV2dl7J0Xu\nKWI7jkqgawIqTQxZNYQhtYZoAah/FMhegNBqoby65lXbUZSLaAmoRDl/cfLj7z8SWj3UdhTlYV6q\n9RJrfl7DjpM7bEdRLqAloP4l3sQzaMUg3q73NhkDM9qOozxM1qCshDnCGLRykF6Y3gdoCah/Cd8b\nToAE0KF8B9tRlIfqUakHpy+dZsmhJbajqFTSElD/41rsNV5Z/QpjGo0hQPTloRIXGBDI6IajeWnl\nS3rNAS+nf+Xqf3y49UMevvdhahepbTuK8nDNSjbj/mz389nOz2xHUamgu4iqf5y7co7SH5VmY4+N\nlM5T2nYc5QV2ntpJ89nNORh6kGwZstmO47d0F1HlEiPWjaBDuQ5aACrJKt9XmYbFGzJ602jbUVQK\n6ZqAAuCn8z9RY2oN9vffT74s+WzHUV7k2B/HqDS5Env67qFA9luvJ6XcQdcEVKq9svoV/lPjP1oA\nKtkK5yjMM5WfYZhzmO0oKgW0BBTbTmxj468bGVhjoO0oyksNfWwoCw8sZP8ZPUO8t9ES8HPGGAav\nHExYnTCyBGWxHUd5qXsy3sPQx4YydNVQ21FUMmkJ+Lmlh5Zy+vJpulfSK4ap1OlfrT97Tu9h/dH1\ntqOoZNAS8GNx8XEMXT2Ud+u/S2BAoO04ystlCMzAW/XeYvCqwXo6CS+iJeDHZu6eSY4MOWhVupXt\nKMpHdKrQiWux1/g26lvbUVQSaQn4qasxV3nD+QajG45GJEV7lin1LwESwKgGo3h59ct6PWIvoSXg\npz7Y+gHVC1SnZqGatqMoH9OoRCOK3FOEqTun2o6ikkAPFvND566co8zHZdjQfYMeHazShJ5Owr30\nYDGVLCPXj+TJB5/UAlBppvJ9lalXrB5jN4+1HUXdha4J+JlfLv5ClU+rsO/Zfdyb9V7bcZQP+/nC\nz1SdUlVfa25gfU1ARJqISLSIHBSRIXcYV01EYkTkCVcsVyXfa2teI7RaqP5RqjRXLGcxulbsyptr\n37QdRd1BqtcERCQAOAjUB04C24EQY0x0IuNWAleBacaYRPch0zWBtPPDqR9oOqsphwYc0nla5Rb/\nPT35pp6bKJW7lO04Psv2mkB14JAx5qgxJgYIB4ITGTcA+Ab43QXLVCkwZNUQXq/9uhaAcpvcmXPz\nYs0XeWX1K7ajqNtwRQkUAH696fbxhPv+ISL3A62NMZ8AulO6BSsPr+SXi7/Qu0pv21GUn3m+xvNs\nPbGVLce32I6iEuGucwWMB27eVnDHIggLC/vne4fDgcPhSJNQ/iLexDNk1RDerv826dOltx1H+ZnM\n6TMz3DGcwSsHs7bbWj040QWcTidOp9Mlj+WKbQI1gDBjTJOE20MBY4wZddOYI//9FsgDXAZ6G2MW\nJvJ4uk3AxWbtmcWH2z5kc8/N+georIiLj+OhSQ/xTv13aFm6pe04Pic12wRcUQLpgAPc2DB8CtgG\ndDTGRN1m/HRgkW4Ydo+/Y/+mzMdlmNF6hl48Xlm1+OBihqwawu6+u/WEhS5mdcOwMSYOCAVWAPuA\ncGNMlIj0EZHEJqD1Hd6NPt7+MeXzldcCUNY1L9mcPJnz8Pmuz21HUTfRg8V82MVrFyn1YSm+e/o7\nyuUrZzuOUmw7sY3g2cHU2dOD08cDKVAggBEjulGsWBHb0bxaatYEdJ3Mh7274V1alW6lBaA8Rt7r\n+flrf3a+Ppoe1ocBl9myZRgrVw7QIrBEzx3ko37941em7JzCcMdw21GU+sfrr3/O5YVzoeZHkPks\nkIXDh4fz+uuf247mt7QEfNQw5zD6VOlDgewF7j5YKTc5cSIezleEvSFQ+62Ee7Nw8mS81Vz+TEvA\nB+39fS9LDi1hSK3bnsZJKSsKFAgALsPaN6Dil5DzCHCZ++/XtyJbdMOwD2oxuwUNizfk+RrP246i\n1P/4+eejNGz4IYcPD4fa4yDvHkrsKqLbBFLJ9rmDlAdZ8/Maos5G0bdqX9tRlPqXYsVuvOF37jyG\n2umvkOnB5Yz5qrYWgEW6JuBD4k08VT+tytDHhtK+XHvbcZS6q6k7p/LFni9wPu3Uo9lTQdcEFACz\nf5xNULog2pVtZzuKUknS/eHunL96nkUHF9mO4re0BHzE1ZirvLL6FcY2GqufqJTXSBeQjtENRzN4\n5WBi4mJsx/FLWgI+YsLWCVQrUI1ahWvZjqJUsjQu0ZhCOQoxZecU21H8km4T8AFnLp/hwY8fZHPP\nzZTMXdJ2HKWSbddvu2jyZRMODjhI9gzZbcfxOrpNwM+9ufZNOlXopAWgvNbD9z5MkweaMGrDqLsP\nVi6lawJe7sDZA9SaVovo0GjyZM5jO45SKXb8z+M8NOkhdvXZRaEchWzH8Sq6JuDHhq4eyuBag7UA\nlNcrmL0g/ar247XvXrMdxa9oCXix9UfXs/PUTp575DnbUZRyicG1BrPi8Ap+OPWD7Sh+Q0vAS8Wb\neF5c8SJv13ubjIEZbcdRyiWyZ8jOG7XfYNDKQei0sHtoCXipOfvmEG/i6Viho+0oSrlUr8q9OPnX\nSSJ/irQdxS9oCXiha7HXeHn1y4xpNIYA0f+FyrekT5ee9xq8x6AVg4iNj7Udx+fpO4gXGrd5HJXu\nrYSjqMN2FKXSRItSLbg/2/1M+n6S7Sg+T3cR9TIn/zpJxU8qsu2ZbRTPWdx2HKXSzI+nf6T+zPpE\nh0aTK1Mu23E8Wmp2EdUS8DLdF3Qnf5b8vNvgXdtRlEpzzy55lsCAQD5o+oHtKB5NS8BPfH/ye1p9\n1Yro0Gg9tF75hTOXz1B2YlnWdltL2bxlbcfxWHqwmB8wxjBw2UBG1B2hBaD8Rt4seXnlsVd4ccWL\ntqP4LJeUgIg0EZFoETkoIv+6sK2IdBKR3QlfG0SkgiuW60++3vc1V2Ku0O3hbrajKOVW/av358iF\nIyw9tNR2FJ+U6hIQkQDgI6AxUA7oKCJlbhl2BKhtjHkIeAvQc8Ymw9WYqwxZNYQJTSaQLiCd7ThK\nuVVQuiDGNRrHC8tf0GsOpAFXrAlUBw4ZY44aY2KAcCD45gHGmC3GmD8Sbm4BCrhguX5jzKYxPFLg\nER4v8rjtKEpZ0axkM4reU5SPt39sO4rPcUUJFAB+ven2ce78Jt8L0EMBk+jEnyeYsHUC7zV8z3YU\npawREcY1HsfI9SM5e+Ws7Tg+JdCdCxORukB34LE7jQsLC/vne4fDgcPhSNNcnmzo6qH0qdKHovcU\ntR1FKavK5i1LSLkQ3vjuDSY2n2g7jlVOpxOn0+mSx0r1LqIiUgMIM8Y0Sbg9FDDGmFG3jKsIzAOa\nGGMO3+HxdBfRBFuPb+WJOU9wIPQAWYOy2o6jlHXnrpzjwY8fZM3Tayifr7ztOB7D9i6i24EHRKSI\niAQBIcDCWwIW5kYBPHWnAlD/Ly4+jgGRA3i73ttaAEolyJ05N2/UeYMBkQP0LKMukuoSMMbEAaHA\nCmAfEG6MiRKRPiLSO2HY60AuYKKI/CAi21K7XF83dedUgtIF8dRDT9mOopRH6Vu1LxevXeSrvV/Z\njuIT9IhhD3Tm8hnKTSzHqq6rqJi/ou04Snmczb9upu2ctkT1jyJHxhy241inp43wMT0X9CR7huy8\n3+R921GU8lj6d/L/tAR8yKZfN9Fubjui+kfp6SGUugNdY/5/tjcMKxeJjY/l2SXPMqbhGC0Ape4i\nb5a8jKg7gmeXPEu8ibcdx2tpCXiQidsnkitTLkLKh9iOopRX6FW5F9fjrjNz90zbUbyWTgd5iFN/\nnaLipIqs67aOB/M+aDuOUl7j+5Pf02J2C6L6R5EzU07bcazQbQI+oMu3XSiYvaBeLEapFHh2ybMA\nfnsksW4T8HLOX5ysP7ae12u/bjuKUl5pZL2RzI+ez/cnv7cdxetoCVgWExdD/6X9eb/x+2QJymI7\njlJeKWemnLxT/x2eXfIscfFxtuN4FS0By8ZtHkfhHIVpU6aN7ShKebWuD3UlKF0Qk3dMth3Fq+g2\nAYsOnjvIo589yvZntlMsZzHbcZTyevvP7KfO53XY0XsHhXMUth3HbXSbgBeKN/H0XNiTN+q8oQWg\nlIuUzVuW5x95nr6L++oJ5pJIS8CST7Z/Qlx8HP2r9bcdRSmfMqTWEE78dYIv93xpO4pX0OkgC45e\nPEqVT6uwvvt6PSZAqTSw4+QOms1uxp6+e8ifNb/tOGlOp4O8iDGGPov78ELNF7QAlEojVe6vQveH\nuzMgcoDtKB5PS8DNZu6eyenLp3np0ZdsR1HKpw2rM4zdp3czP2q+7SgeTaeD3Oi3S79R8ZOKLO+y\nnEr3VbIdRymft/7oekLmhbC3316fPqWEnjbCSzw550lK5S7F2/Xfth1FKb8RujSUKzFXmBY8zXaU\nNKPbBLzAvP3z2HdmH2/UecN2FKX8yjv132HNz2tYcXiF7SgeSUvADc5fPc+AyAFMbTmVjIEZbcdR\nyq9ky5CNyS0m03tRb/76+y/bcTyOTgelMWMM7ea2o1D2QnoZPKUs6rmgJwbjk9NCOh3kwab9MI1D\n5w/xToN3bEdRyq9NaDqBDcc2MGffHNtRPIquCaShA2cPUGtaLdZ1X0fZvGVtx1HK720/sZ3ms5vz\nfe/vfercQrom4IGux12n07edGFF3hBaAUh6iWoFqvFjzRbp820VPOZ3AJSUgIk1EJFpEDorIkNuM\n+UBEDonILhF52BXL9WSvrXmNgtkL0rdqX9tRlFI3eanWSwQGBPLOBp2iBReUgIgEAB8BjYFyQEcR\nKXPLmKZACWNMSaAPMCm1y/Vkq46sYvaPs/ms1WeIpGgNTSmVRgIkgJltZvLhtg/ZcnyL7TjWuWJN\noDpwyBhz1BgTA4QDwbeMCQZmAhhjtgI5RMQnz+p09spZukV0Y3rwdPJkzmM7jlIqEQWzF2RS80l0\n/rYzf/79p+04VrmiBAoAv950+3jCfXcacyKRMV7PGEOvhb3oWL4jDUs0tB1HKXUHbR5sQ4NiDQhd\nGmo7ilWBtgMkJiws7J/vHQ4HDofDWpbkmLxjMsf+OMacdroLmlLeYFzjcVT5tAqzf5xNpwqdbMdJ\nMqfTidPpdMljpXoXURGpAYQZY5ok3B4KGGPMqJvGTAK+M8Z8nXA7GqhjjDmdyON55S6iu3/bTYMv\nGrC++3rK5Clz9x9QSnmEH079QKMvG7Gxx0ZK5S5lO06K2N5FdDvwgIgUEZEgIARYeMuYhUBX+Kc0\nLiZWAN7q7JWztPm6DR80+UALQCkvU+m+SoysN5Lg8GC/3D7gkoPFRKQJMIEbpfKZMeZdEenDjTWC\nTxPGfAQ0AS4D3Y0xO2/zWF61JhAbH0vjLxtT9b6qjGo46u4/oJTySP0W9+PEXyeICIkgQLzrECo9\nlbRFA5cN5MC5AyzuuJh0Aelsx1FKpdD1uOvUn1mfukXr8mbdN23HSRbb00F+a8auGSw9tJTZT8zW\nAlDKywWlC+Kbdt8wY/cMvo361nYct9E1gRTadmIbLWa3wNnNqaeFUMqH7Di5gyazmrCm6xoq5K9g\nO06S6JqAm/126TfazmnLlJZTtACU8jFV7q/C+43fp/XXrTl/9bztOGlOSyCZ/o79m7Zz2tKrUi+C\ny9x6YLRSyhd0qdiF1qVb0+GbDsTGx9qOk6Z0OigZjDH0WdyHM1fOMK/9PK/bg0AplXSx8bE0m9WM\n8vnKM67xONtx7king9wkzBnG9pPbmdl6phaAUj4uMCCQ8CfDifwpktEbR9uOk2Y88rQRnmjc5nF8\nve9r1nVfR7YM2WzHUUq5Qa5MuVj51Eoen/44OTLmoHeV3rYjuZyWQBJM3TmVD7Z+wPru68mXJZ/t\nOEopNyqYvSArn1pJnc/rkD1DdkLKh9iO5FJaAncxZ98chjmH4XzaSaEchWzHUUpZ8ECuB1jWeRkN\nvmhAtqBsNC/V3HYkl9GJ7TtYemgpAyIHENk5kpK5S9qOo5SyqEL+CiwMWUj3Bd1Z+8ta23FcRkvg\nNtYdXUe3iG4sCFlAxfwVbcdRSnmARwo+QviT4bSb247tJ7bbjuMSWgKJ2HFyB0/OeZKv2n5FjYI1\nbMdRSnkyEJjIAAAJbElEQVSQesXqMbXVVFp+1ZJ9v++zHSfVdJvALZy/OOnwTQemtJxC/eL1bcdR\nSnmgVqVbcen6JRp80YCIDhE8UvAR25FSTNcEbhK+N5z2c9sz+4nZejSwUuqOOlXoxJSWU2j5VUsW\nHrj1EireQ48Y5saRwGM2jeHDbR+ypNMSrzlplFLKvu0nthMcHszrtV+nX7V+VjLo9QRSIS4+joHL\nBuI86iSycyQFsxd027KVUr7h8PnDNJ3VlLYPtmVk/ZFuP6OAlkAKXYm5QudvO/Pn33/ybftvyZEx\nh1uWq5TyPWevnKXVV60onrM404KnEZQuyG3L1nMHpcDZK2epP7M+WdJnIbJzpBaAUipV8mTOw+qu\nq7kSc4UmXzbh4rWLtiMliV+WwOojq6k8uTJ1i9blizZfuLWxlVK+K1P6TMxtN5cK+SpQ5dMqbDy2\n0Xaku/Kr6aCrMVd5efXLfLP/Gz5r9RmNH2icJstRSqmI6Aj6Lu5Lj0o9CHOEpemHTZ0OSoIdJ3dQ\n5dMqnLp0it19d2sBKKXSVOsyrdnddzd7f9/LI1MfYe/ve21HSpTPl0BsfCxvrXuLprOa8lrt1whv\nG07uzLltx1JK+YH8WfOzIGQBodVCqTujLuM2jyPexNuO9T98ejpo3+/76LWoF1mDsjI9eLru/qmU\nsubIhSN0nd+V9OnSM7nFZErlLuWyx7Y2HSQiOUVkhYgcEJHlIvKvXWxEpKCIrBGRfSLyo4g8l5pl\nJsX+M/sJ+SaEejPr0blCZ5Z3Wa4FoJSyqnjO4qzttpbmJZtTa1otno54mp/O/2Q7Vqqng4YCq4wx\npYE1wMuJjIkFXjDGlANqAv1FpEwql5uoqDNRdJrXiboz6lLp3kocfu4wodVD9VKQSimPkC4gHYMe\nHcShAYcokbMENT+rSfcF3Tl8/rC1TKmaDhKRaKCOMea0iNwLOI0xd3yDF5EI4ENjzOrb/Huyp4Oi\nz0YzYt0IVh5eyQs1X6B/tf56CUillMe7eO0iE7ZM4MNtH9KqdCteffxVSuQqkezHsXbEsIicN8bk\nut3tRMYXBZxAeWPMpduMSVIJHPvjGBHREURER7D3970MrDGQAdUH6Ju/UsrrXLx2kfFbxvPRto+o\nfF9lWpdpTXDpYApkL5Ckn0/TEhCRlUD+m+8CDPAa8PktJXDOGJPorjcikpUbBTDCGLPgDsszw4YN\n++e2w+HA4XBgjGHfmX1EREcwP3o+Ry8epWXplrQu3ZqGJRqSOX3mu/+2SinlwS5dv8Tyn5YTcSCC\nJQeXUCp3KVqXaU2bMm0onaf0P+OcTidOp/Of28OHD7e2JhAFOG6aDvrOGPNgIuMCgcVApDFmwl0e\n0xw4e4AjF45w+Pxhjlw4wpGLR9hzeg8xcTG0KdOGNg+24bHCjxEYoJdDUEr5ppi4GJy/OJkfPZ+I\n6AiyBmWlfL7ylMhZguI5i1Mi143/Fs5RmAyBGayVwCjgvDFmlIgMAXIaY4YmMm4mcNYY80ISHtMU\nn1D8xi/53182ZwlK5ylNubzlEEnR76mUUl4r3sTz4+kfOXju4I0PyBcSPiBfOMKJv05w/fXr1kog\nFzAHKAQcBdobYy6KyH3AFGNMCxGpBawDfuTGNJIBXjHGLLvNY7r9egJKKeWtYuJiCAoM0lNJK6WU\nv9JzBymllEoRLQGllPJjWgJKKeXHtASUUsqPaQkopZQf0xJQSik/piWglFJ+TEtAKaX8mJaAUkr5\nMS0BpZTyY1oCSinlx7QElFLKj2kJKKWUH9MSUEopP6YloJRSfkxLQCml/JiWgFJK+TEtAaWU8mNa\nAkop5ce0BJRSyo9pCSillB/TElBKKT+WqhIQkZwiskJEDojIchHJcYexASKyU0QWpmaZSimlXCe1\nawJDgVXGmNLAGuDlO4x9HtifyuWpZHA6nbYj+BR9Pl1Ln0/PkNoSCAZmJHw/A2id2CARKQg0A6am\ncnkqGfSPzLX0+XQtfT49Q2pLIJ8x5jSAMeY3IN9txr0PvASYVC5PKaWUCwXebYCIrATy33wXN97M\nX0tk+L/e5EWkOXDaGLNLRBwJP6+UUsoDiDEp/3AuIlGAwxhzWkTuBb4zxjx4y5i3gS5ALJAJyAZ8\na4zpepvH1LUFpZRKJmNMij5gp7YERgHnjTGjRGQIkNMYM/QO4+sALxpjWqV4oUoppVwmtdsERgEN\nReQAUB94F0BE7hORxakNp5RSKm2lak1AKaWUd7N6xLCIPCkie0UkTkQq32FcExGJFpGDCdNOKhFJ\nPXhPRH4Rkd0i8oOIbHN3Tk+XlNebiHwgIodEZJeIPOzujN7ibs+liNQRkYsJB5LuFJHEdjhRCUTk\nMxE5LSJ77jAmWa9N26eN+BFoA6y93QARCQA+AhoD5YCOIlLGPfG8TlIP3ovnxgb9SsaY6m5L5wWS\n8noTkaZACWNMSaAPMMntQb1AMv521xljKid8veXWkN5nOjeez0Sl5LVptQSMMQeMMYe4826j1YFD\nxpijxpgYIJwbB6mpf0vSwXvceL5tfwDwVEl5vQUDMwGMMVuBHCKSH3WrpP7t6m7jSWSM2QBcuMOQ\nZL82veGNoADw6023jyfcp/4tqQfvGWCliGwXkWfcls47JOX1duuYE4mMUUn/262ZMHWxRETKuiea\nz0r2a/OuB4ul1h0ONnvVGLMorZfva1J78F6CWsaYUyKSlxtlEJXwCUMpd9sBFDbGXEmYyogASlnO\n5FfSvASMMQ1T+RAngMI33S6YcJ9futPzmbDBKP9NB+/9fpvHOJXw3zMiMp8bq+1aAjck5fV2Aih0\nlzEqCc+lMebSTd9HishEEclljDnvpoy+JtmvTU+aDrrdvOB24AERKSIiQUAIoKejTtxCoFvC908D\nC24dICKZRSRrwvdZgEbAXncF9AJJeb0tBLoCiEgN4OJ/p+HU/7jrc3nzfLWIVOfGbutaAHcm3P79\nMtmvzTRfE7gTEWkNfAjkARaLyC5jTFMRuQ+YYoxpYYyJE5FQYAU3SuszY0yUxdiebBQwR0R6AEeB\n9nDj4D0Snk9uTCXNTzg9RyAwyxizwlZgT3O715uI9Lnxz+ZTY8xSEWkmIj8Bl4HuNjN7qqQ8l8CT\nItIPiAGuAh3sJfZ8IjIbcAC5ReQYMAwIIhWvTT1YTCml/JgnTQcppZRyMy0BpZTyY1oCSinlx7QE\nlFLKj2kJKKWUH9MSUEopP6YloJRSfkxLQCml/Nj/AZX/COSBywoZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7faff97d9dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's try it:\n",
    "\n",
    "# The function f can be evaluated at any point\n",
    "x0 = 0.5\n",
    "y0 = 0.5\n",
    "val = f(x0,y0)\n",
    "print('Value of f at ({},{}) : {}'.format(x0,x0,val))\n",
    "\n",
    "# Or at any list of points\n",
    "x = linspace(-1,1,50)\n",
    "y = linspace(-1,1,50)\n",
    "vals = f(x,y)\n",
    "\n",
    "# Plot\n",
    "plot(x0, val, 'o')\n",
    "plot(x, f(x,y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The grid over which we approximate the function is defined by its boundaries $[a_i, b_i]$ and the number of points $n_i$ in each dimension $i \\in [1,d]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In each dimension $i\\in[1,d]$ the $n_i$ points $(a_i=s^0_i, s^2_i, ..., s^{n_j-1}_i=b_i)$\n",
    "are regularly spaced, meaning that for any $j$, $s^j_i = a_i + j\\frac{b_i-a_i}{n_i-1}$. These points are called *nodes*.\n",
    "\n",
    "For the two dimensional function defined above, we define $[a_1,b_1]=[a_2,b_2]=[-5,5]$ and we use 50 points in each dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define the approximation space\n",
    "a = [-5.0, -5.0]  # a_1, a_2\n",
    "b = [5.0, 5.0]    # b_1, b_2\n",
    "orders = [15,15]     # n_2, n_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cartesian product of these nodes is used to approximate the compact state $\\mathcal{S}=[a_1, b_1], ... [a_d, b_d]$, by a finite set of points $S=(s_{i_1,...i_d})_{i_1\\in[1,n_1], ..., i_d\\in[1,n_d]}$. By definition, the point indexed by i_1,...i_d is simply the tuple $(s_{i_1},...s_{i_d})$.\n",
    "\n",
    "In general, this multivariate index is sufficient and it is not needed to list these points one after all. By convention, we list the points as a 2d arrays, where each line corresponds to a point, and where *the last index varies faster* meaning that points are listed like in $s_{00}, s_{01}, s_{02}, s_{10}, s_{11}, s_{12}$.\n",
    "\n",
    "This grid can be constructed using the mlinspace function, which extends the linspace function to multidimensional spaces:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First points:\n",
      "[[-5.         -5.        ]\n",
      " [-5.         -4.28571429]\n",
      " [-5.         -3.57142857]\n",
      " [-5.         -2.85714286]\n",
      " [-5.         -2.14285714]\n",
      " [-5.         -1.42857143]\n",
      " [-5.         -0.71428571]\n",
      " [-5.          0.        ]\n",
      " [-5.          0.71428571]\n",
      " [-5.          1.42857143]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fafeaed2240>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAElNJREFUeJzt3X+s3fVdx/HXqxYMjkF0Jhja8FPRhS0iOpgh6skIsWGx\n/mFMZjAL21/qlEXIAoMtvQ0xwWkyceofZFsNusUIM44ZWOiCZyYmGH6spBYILAVSSsaCmhhmsoB9\n+8c9XK7d/fH59ny+38+nn8/zkZyktz332fe5HZ+dfu857zoiBAA4ve0oPQAAYHkc5gDQAA5zAGgA\nhzkANIDDHAAawGEOAA3IcpjbPtf2fbafsX3E9tU5ugCANDszde6W9GBE/KbtnZJ+JFMXAJDAy75p\nyPY5kr4VEZfmGQkAMFSOyywXS3rN9gHbT9q+x/ZZGboAgEQ5DvOdkq6U9JcRcaWk/5F0W4YuACBR\njmvmL0s6FhGPLz6+X9KtJ9/JNktgAOAURIS3u8/Sz8wj4lVJx2xftvipayU9vcl9m73t27ev+Aw8\nPh4bj6+9W6pcr2a5SdKXbJ8h6aikj2TqAgASZDnMI+IpSe/L0QIADMc7QDOZzWalRxhVy4+v5ccm\n8fh6sfTrzJN/Izum+r0AoBW2FVN8AxQAUB6HOQA0gMMcABrAYQ4ADeAwB4AGcJgDQAM4zAGgARzm\nANAADnMAaECuRVvVeOGFl/TpT/+1jh8/oV27dujOO2/UxRdfSK9Ar+bZ6NEr3ctuwjWOMbajR1+M\nSy+9JaTXQ4qQXo9LL70ljh59kd7EvZpno0evdG+Ixdm5/RmbcqcctykO8xtuWFn3xY61L/oNN6zQ\nm7hX82z06JXuDZF6mDd1zfz48ROS3nHSz75Dr7xygt7EvZpno0evdG8MTR3mu3btkPS9k372ezr/\n/FN7mPROvVfzbPTole6NIuXpe46buGbeVa/m2ejRK90bQomXWZrbZ/7Wd5xfeeWEzj8/33ew6bU1\nGz16pXupUveZN3eYA0BL+McpAKAjHOYA0AAOcwBoAIc5ADSAwxwAGsCiLXqj9WqejR690r3sUl6M\nnuMm3jTUVa/m2ejRK90bQizairUvei3LdXrq1TwbPXqle0OkHuZNXTOvfblOT72aZ6NHr3RvDNkO\nc9s7bD9p+4FczaFqX67TU6/m2ejRK90bRcrT95SbpD+U9LeSHtjk10f/60jt18l66tU8Gz16pXtD\naMpFW7Z3Szog6Y8k3RwReze4T+T4vbZT+3Kdnno1z0aPXuleqkkXbdm+T6sH+bmSbil5mANASyZb\ntGX7g5JejYhDkry4AQAmlONNQ9dI2mv7eklnSXqn7Xsj4sMn33FlZWXtx7PZTLPZLMNvDwDtmM/n\nms/ngz8v6z5z278iLrMAQDbsMweAjvAvDQFAxXhmDgAdYWsivdF6Nc9Gj17pXnYp7yzKcRPvAO2q\nV/Ns9OiV7g0htibG2he9lk1pPfVqno0evdK9IVIP86aumde+Ka2nXs2z0aNXujeGpg7z2jel9dSr\neTZ69Er3RpHy9D3HTVwz76pX82z06JXuDaEptyamYGtif72aZ6NHr3Qv1aRbE1PwpiEAGI43DQFA\nRzjMAaABHOYA0AAOcwBoAIc5ADSARVv0RuvVPBs9eqV72aW8GD3HTbxpqKtezbPRo1e6N4RYtBVr\nX/Raluv01Kt5Nnr0SveGSD3Mm7pmXvtynZ56Nc9Gj17p3hiaOsxrX67TU6/m2ejRK90bRcrT9xw3\ncc28q17Ns9GjV7o3hFi0VedynZ56Nc9Gj17pXioWbQFAA1i0BQAd4TAHgAZwmANAAzjMAaABHOYA\n0AAWbdEbrVfzbPTole5ll/Ji9Bw38aahrno1z0aPXuneEJpq0Zak3ZIekXRE0mFJN21yv9EfdO3L\ndXrq1TwbPXqle0OkHuY5LrO8KenmiDhk+2xJT9h+OCKezdAepPblOj31ap6NHr3SvTEs/Q3QiPhO\nRBxa/Ph1Sc9I2rVs91TUvlynp17Ns9GjV7o3ipSn76k3SRdJelHS2Rv82ph/E4mI+q+T9dSreTZ6\n9Er3htDUi7YWl1jmku6MiK9u8Ouxb9++tY9ns5lms1mW33u92pfr9NSreTZ69Er3NjOfzzWfz9c+\n3r9/v2KqRVu2d0r6J0kPRcTdm9wncv0fBwD0YtKtibbvlfRaRNy8xX04zAFgoMkOc9vXSPoXrb4s\nMRa32yPi6yfdj8McAAZinzkANIB95gDQEQ5zAGgAhzkANICtifRG69U8Gz16pXvZpbyzKMdNvAO0\nq17Ns9GjV7o3hKbamph6m+Iwr31TWk+9mmejR690b4jUw7ypa+a1b0rrqVfzbPTole6NoanDvPZN\naT31ap6NHr3SvVGkPH3PcRPXzLvq1TwbPXqle0No6q2J25nqHaC1b0rrqVfzbPTole6l4u38ANAA\n3s4PAB3hMAeABnCYA0ADOMwBoAEc5gDQABZt0RutV/Ns9OiV7mWX8mL0HDfxpqGuejXPRo9e6d4Q\nYtFWrH3Ra1mu01Ov5tno0SvdGyL1MG/qmnnty3V66tU8Gz16pXtjaOowr325Tk+9mmejR690bxQp\nT99z3MQ18656Nc9Gj17p3hBi0Vady3V66tU8Gz16pXupWLQFAA1g0RYAdITDHAAawGEOAA3gMAeA\nBmQ5zG3vsf2s7eds35qjCQBIt/SrWWzvkPScpGslvSLpMUkfiohnT7rfpC9NrHW5Tk+9mmejR690\nL1Xqq1lyvBno/ZIeWvfxbZJu3eB+o7ygfr3a3yjQU6/m2ejRK90bQlMt2pL0G5LuWffxb0v68w3u\nN/qDrn25Tk+9mmejR690b4jUw3zSfeYrKytrP57NZprNZln7tS/X6alX82z06JXubWU+n2s+nw/+\nvByH+XFJF6z7ePfi537A+sN8DG8vw1n/Rc+xXIdeS7PRo1e6t5WTn+ju378/7RNTnr5vdZP0Q5K+\nLelCSWdKOiTp3Rvcb/S/jtR+naynXs2z0aNXujeEply0ZXuPpLu1+lLHL0TEXRvcJ3L8XtupfblO\nT72aZ6NHr3QvFYu2AKABLNoCgI5wmANAAzjMAaABHOYA0AAOcwBowKTvAJ1C7ct1eurVPBs9eqV7\n2aW8GD3HTbxpqKtezbPRo1e6N4SmWrSVepviMK99uU5PvZpno0evdG+I1MO8qWvmtS/X6alX82z0\n6JXujaGpw/ztZTjr5ViuQ6+l2ejRK90bRcrT9xw3cc28q17Ns9GjV7o3hKZctJWCRVv99WqejR69\n0r1ULNoCgAawaAsAOsJhDgAN4DAHgAZwmANAAzjMAaABHOYA0AC2JtIbrVfzbPTole5ll/LOohw3\n8Q7Qrno1z0aPXuneEGJrYqx90WvZlNZTr+bZ6NEr3Rsi9TBv6pp57ZvSeurVPBs9eqV7Y2jqMK99\nU1pPvZpno0evdG8UKU/fc9zENfOuejXPRo9e6d4QYmtinZvSeurVPBs9eqV7qdiaCAANYGsiAHRk\nqcPc9mdsP2P7kO2v2D4n12AAgHTLPjN/WNLlEXGFpOclfXL5kQAAQy11mEfENyLirRdaPipp9/Ij\nAQCGynnN/KOSHsrYAwAk2nbRlu2Dks5b/1OSQtIdEfG1xX3ukPRGRHx5q9bKysraj2ezmWaz2fCJ\nt1H7cp2eejXPRo9e6d5m5vO55vP58E9MeTH6VjdJN0r6V0k/vM39xng9/f9T+xsFeurVPBs9eqV7\nQ2iKRVuS9kg6IuldCfcd/UHXvlynp17Ns9GjV7o3ROphvuw1889JOlvSQdtP2v6rJXtLqX25Tk+9\nmmejR690bwxL/eMUEfFTuQbJ4e1lOOu/6DmW69BraTZ69Er3RpHy9D3HTVwz76pX82z06JXuDSEW\nbdW5XKenXs2z0aNXupeKRVsA0AAWbQFARzjMAaABHOYA0AAOcwBoAIc5ADRgqTcN1aj25To99Wqe\njR690r3sUl6MnuMm3jTUVa/m2ejRK90bQlMs2hpym+Iwr325Tk+9mmejR690b4jUw7ypa+a1L9fp\nqVfzbPTole6NoanD/O1lOOvlWK5Dr6XZ6NEr3RtFytP3HDdxzbyrXs2z0aNXujeEWLRV53Kdnno1\nz0aPXuleKhZtAUADWLQFAB3hMAeABnCYA0ADOMwBoAEc5gDQAA5zAGgAWxPpjdareTZ69Er3skt5\nZ1GOm3gHaFe9mmejR690bwixNTHWvui1bErrqVfzbPTole4NkXqYN3XNvPZNaT31ap6NHr3SvTE0\ndZjXvimtp17Ns9GjV7o3ipSn7zlu4pp5V72aZ6NHr3RvCE25NdH2LZL+RNKPR8R/bnKfyPF7baf2\nTWk99WqejR690r1Uk21NtL1b0ucl/bSkny99mANAS6bcmvhZSZ/I0AEAnKKlDnPbeyUdi4jDmeYB\nAJyCbd8BavugpPPW/5SkkPQpSbdLuu6kXwMATGzbwzwirtvo522/R9JFkp6ybUm7JT1h+6qI+O5G\nn7OysrL249lsptlsNnxiAGjYfD7XfD4f/HnZ/tk42y9IujIi/muTX+cboAAwUOo3QHMu2gpVcJml\n9uU6PfVqno0evdK97FJejJ7jJt401FWv5tno0SvdG0Is2oq1L3oty3V66tU8Gz16pXtDpB7mFS0W\nWF7ty3V66tU8Gz16pXtjaOowr325Tk+9mmejR690bxQpT99z3MQ18656Nc9Gj17p3hCactFWChZt\n9dereTZ69Er3Uk22aCsVrzMHgOGmXLQFACiMwxwAGsBhDgAN4DAHgAZwmANAA3Iu2qpC7ct1eurV\nPBs9eqV72aW8GD3HTbxpqKtezbPRo1e6N4RYtBVrX/Raluv01Kt5Nnr0SveGSD3Mm7pmXvtynZ56\nNc9Gj17p3hiaOsxrX67TU6/m2ejRK90bRcrT9xw3cc28q17Ns9GjV7o3hFi0VedynZ56Nc9Gj17p\nXioWbQFAA1i0BQAd4TAHgAZwmANAAzjMAaABHOYA0AAOcwBoAIc5ADSAwxwAGsBhDgANWPowt/0H\ntp+xfdj2XTmGAgAMs9Rhbnsm6dckvTci3ivpT3MMdTqaz+elRxhVy4+v5ccm8fh6sewz89+VdFdE\nvClJEfHa8iOdnlr/H1TLj6/lxybx+Hqx7GF+maRftv2o7X+2/Qs5hgIADLPtP+hs+6Ck89b/lKSQ\n9KnF5/9oRLzf9vsk/b2kS8YYFACwuaVW4Np+UNIfR8Q3Fx9/W9LVEfEfG9yX/bcAcApSVuBu+8x8\nG/8o6QOSvmn7MklnbHSQpw4DADg1yx7mByR90fZhSd+X9OHlRwIADDXZvzQEABjP5O8Abf1NRrZv\nsX3C9o+VniUn259Z/Lkdsv0V2+eUnikH23tsP2v7Odu3lp4nJ9u7bT9i+8jiv7ebSs+Um+0dtp+0\n/UDpWXKzfa7t+xb/3R2xffVW95/0MG/9TUa2d0u6TtJLpWcZwcOSLo+IKyQ9L+mThedZmu0dkv5C\n0q9KulzSb9n+mbJTZfWmpJsj4nJJvyjpY409Pkn6uKSnSw8xkrslPRgR75b0s5Ke2erOUz8zb/1N\nRp+V9InSQ4whIr4REScWHz4qaXfJeTK5StLzEfFSRLwh6e8k/XrhmbKJiO9ExKHFj1/X6mGwq+xU\n+SyePF0v6fOlZ8lt8TffX4qIA5IUEW9GxH9v9TlTH+bNvsnI9l5JxyLicOlZJvBRSQ+VHiKDXZKO\nrfv4ZTV02K1n+yJJV0j6t7KTZPXWk6cWv/F3saTXbB9YXEa6x/ZZW33Csq9m+QEtv8lom8d2u1Yv\nsaz/tdPKFo/vjoj42uI+d0h6IyK+XGBEnALbZ0u6X9LHF8/QT3u2Pyjp1Yg4tLh8e9r997aNnZKu\nlPSxiHjc9p9Juk3Svq0+IauIuG6zX7P9O5L+YXG/xxbfKHzXZq9Nr81mj832eyRdJOkp29bqJYgn\nbF8VEd+dcMSlbPVnJ0m2b9TqX2s/MMlA4zsu6YJ1H+9e/FwzbO/U6kH+NxHx1dLzZHSNpL22r5d0\nlqR32r43Ilp5efTLWv2b/uOLj++XtOU36Ke+zPLWm4y03ZuMTicR8e8R8RMRcUlEXKzVP4ifO50O\n8u3Y3qPVv9LujYjvl54nk8ck/aTtC22fKelDklp7VcQXJT0dEXeXHiSniLg9Ii6IiEu0+uf2SEMH\nuSLiVUnHFuekJF2rbb7Rm/2Z+TZ6eZNRqL2/9n1O0pmSDq7+5UOPRsTvlR1pORHxv7Z/X6uv1Nkh\n6QsRseUrBk4ntq+RdIOkw7a/pdX/Xd4eEV8vOxkS3STpS7bPkHRU0ke2ujNvGgKABvDPxgFAAzjM\nAaABHOYA0AAOcwBoAIc5ADSAwxwAGsBhDgAN4DAHgAb8H5/rHJ/nYRzTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fafecfb94e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from dolo.numeric.misc import mlinspace\n",
    "\n",
    "list_of_points = mlinspace(a,b,orders)\n",
    "print('First points:')\n",
    "print(list_of_points[:10,:])\n",
    "\n",
    "# optional plot the grid\n",
    "plot(list_of_points[:,0], list_of_points[:,1],'o')\n",
    "\n",
    "\n",
    "# vectorized version\n",
    "#all_points =  [ [x y] for x=linspace(smin[1],smax[1],orders[1]),  y=linspace(smin[2],smax[2],orders[2])] \n",
    "#all_points = vcat(all_points...)\n",
    "#PyPlot.plot(all_points[:,1], all_points[:,2], \"x\", color=\"blue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpolating data is supposed to fit the exact values $v_{i_1, ..., i_d}=f(s_{i_1, ..., i_d})$ of the function on the grid. We compute them below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_array = list_of_points[:,0]\n",
    "y_array = list_of_points[:,1]\n",
    "vals = f(x_array, y_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To intepolate this function, we construct an interpolating object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from interpolation.splines import CubicSpline\n",
    "interp = CubicSpline(a,b,orders,vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.980428387366858"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# we can evaluate this interpolator at a given point:\n",
    "x0 = array([0.0,0.1])\n",
    "interp(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Array of points'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "array([[-0.        , -6.        ],\n",
       "       [-0.        , -5.87878788],\n",
       "       [-0.        , -5.75757576],\n",
       "       [-0.        , -5.63636364],\n",
       "       [-0.        , -5.51515152],\n",
       "       [-0.        , -5.39393939],\n",
       "       [-0.        , -5.27272727],\n",
       "       [-0.        , -5.15151515],\n",
       "       [-0.        , -5.03030303],\n",
       "       [-0.        , -4.90909091]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Array of interpolated values'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "array([-0.1852616 , -0.16280564, -0.14034969, -0.11789374, -0.09543779,\n",
       "       -0.07298184, -0.05052589, -0.02806994, -0.00561399,  0.01668887])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# to evaluate it many times, it is more efficient to supply an array of points\n",
    "# where each line containt the coordinates of one point\n",
    "vec = linspace(-6.0, 6.0, 100)\n",
    "points = column_stack([ vec*0, vec ])\n",
    "interpolated_values = interp(points)\n",
    "display(\"Array of points\")\n",
    "display((points[:10,:])) # first 10 points\n",
    "display(\"Array of interpolated values\")\n",
    "display(interpolated_values[:10]) # first 10 interpolated values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Setting\n",
      "[15 15]\n",
      "(15, 15)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fafeac59240>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lFX68PHvCZ0okICEHpp0UAQhIEKCsoCKolhABFFW\nxVewrquoCC4WXHEVdVWw+0MF22IDEYWoqFGQGkgASQiYAFIFQksy9/vHGWJ6fabfn+uai8zMM+fc\nM8zcc+Y8pxgRQSmlVPAL83UASimlvEMTvlJKhQhN+EopFSI04SulVIjQhK+UUiFCE75SSoUIRxK+\nMeY1Y8xuY8y6Yu6/1hiz1n1Zbozp6kS9Simlys6pFv4bwOAS7k8B+ovIWcCjwCsO1auUUqqMqjpR\niIgsN8ZEl3B/Qp6rCUBTJ+pVSilVdr7ow/87sMgH9SqlVEhzpIVfVsaYOOAGoJ8361VKKeXFhG+M\n6QbMAYaIyIESjtPFfZRSqpxExJR2jJNdOsZ9KXyHMS2Aj4AxIrK1tIJEJCgvU6dO9XkM+vz0+enz\nC75LWTnSwjfGvAvEAvWNMduBqUB1m7tlDjAFiAReNMYYIEtEejlRt1JKqbJxapTOtaXcfxNwkxN1\nKaWUqhidaetFsbGxvg7Bo/T5BTZ9fsHPlKf/xxuMMeJvMSmllD8zxiBlOGnr1WGZSgW6li1bkpaW\n5uswVIiKjo5m27ZtFX68tvCVKgd3S8rXYagQVdz7r6wtfO3DV0qpEKEJXymlQoQmfKWUChGa8JUK\nIl26dOG7777zdRhl0qpVK5YuXeqVum644QYefvhhx8uNi4vj9ddfd7xcT9GEr1QQSUxMpH///mU6\n1psJt7LCwsJISUnxdRgBTxO+UqpCvDlaya7IoipLE75SQSRvq/2RRx7hmmuu4frrr6dOnTp07dqV\nVatWATB27Fi2b9/OsGHDqFOnDjNnzgQgISGB8847j4iICLp37863336bW3ZcXBwPPfQQ/fr1Izw8\nnNTUVOLi4njggQfo3bs3devW5fLLL+fgwYO5j/n000/p0qULkZGRDBw4kOTk5CLjXrFiBX379iUi\nIoKmTZsyadIksrOzARgwYAAiQrdu3ahTpw4ffPABAJ9//jndu3cnIiKCfv36sX79+tzyVq9eTY8e\nPahbty4jR47k+PHjRdZ78uRJIiIi2LhxY+5te/fupXbt2uzdu5eDBw8ybNgwGjZsSP369Rk2bBjp\n6elFlvXII48wZsyY3OtpaWmEhYXhcrkAOHToEH//+99p0qQJzZs3Z8qUKblfmlu3biU2NpZ69erR\nsGFDRo0aVWQdlebrVd6KWPVNlPJX/v7+bNmypXzzzTciIjJt2jSpVauWfPnll+JyuWTy5MkSExOT\n79ilS5fmXk9PT5f69evLl19+KSIiX3/9tdSvX1/27t0rIiKxsbESHR0tSUlJkpOTI1lZWRIbGyvN\nmjWTjRs3ytGjR2XEiBFy3XXXiYjIpk2bJDw8XL755hvJzs6Wf//739K2bVvJysoqFOuvv/4qP//8\ns7hcLklLS5NOnTrJrFmzcmMzxkhKSkru9VWrVknDhg1lxYoV4nK55O2335aWLVvKyZMn5eTJkxId\nHS2zZs2S7Oxs+fDDD6VatWoyZcqUIl+z8ePHy0MPPZR7/b///a8MHTpURET27dsnH3/8sRw/flyO\nHDkiV199tQwfPjz32NjYWHnttddyX+8xY8bk3rdt2zYJCwuTnJwcEREZPny43HrrrXLs2DHZs2eP\n9O7dW+bMmSMiIqNGjZLHH39cREROnDghP/zwQ5GxFvf+c99ean7VFr5SDjPGmYsT+vXrx+DBgzHG\nMGbMGNatW5fvfsnTLTN37lwuvvhiBg+221NfcMEF9OzZk4ULF+YeM27cODp06EBYWBhVq9qJ+mPG\njKFjx47UqlWL6dOn88EHHyAivP/++1xyySUMHDiQKlWq8I9//INjx47x448/ForznHPOoVevXhhj\naNGiBTfffHO+XxcFY33llVeYMGECPXv2zH1uNWrUICEhgYSEBLKzs7n99tupUqUKI0aM4Nxzzy32\nNRo1ahTvvfde7vV3332Xa6+160FGRkZy+eWXU6NGDcLDw5k8eXKFTorv3r2bRYsW8cwzz1CzZk0a\nNGjAnXfeybx58wCoVq0aaWlppKenU716dfr27VvuOspCE75SDhNx5uKERo0a5f5du3Ztjh8/ntvF\nUFBaWhrvv/8+kZGRREZGEhERwQ8//MCuXbtyj2nevHmhx+W9LTo6mqysLPbu3UtGRgbR0X9tdW2M\noXnz5kV2iWzZsoVhw4bRuHFj6tWrx4MPPsjevXuLfV5paWk8/fTT+WL9/fffycjIICMjg6ZN82+b\nnTeOguLi4jh27BgrVqwgLS2NtWvXcvnllwNw7NgxbrnlFlq2bEm9evUYMGAABw8eLPf5i+3bt5OV\nlUXjxo1z450wYQJ79uwB4KmnnsLlctGrVy+6du3KG2+8Ua7yy0rX0lEqRBU8Edq8eXPGjh3L7Nmz\ny/wYgB07duT+nZaWRrVq1WjQoAFNmjQhMTGx0LHNmjUrVMatt97KOeecw/z586lduzazZs3io48+\nKjaO5s2b8+CDDzJ58uRC93333XeFvlS2b99O27ZtiywrLCyMq6++mnfffZeoqCguueQSwsPDAXj6\n6afZsmULK1as4IwzzmDt2rWcc845iEih1yI8PJyjR4/mXt+5c2e+eGvWrMm+ffuKfA0bNmzInDlz\nAPjhhx+48MILGTBgAK1bty72NagIbeErFULytkwbNWqUb6jjddddx2effcZXX32Fy+Xi+PHjfPvt\nt2RkZJRY5ty5c0lOTubo0aNMnTqVq666CmMMV199NV988QXLli0jOzubmTNnUrNmTfr06VOojMOH\nD1OnTh1q165NcnIyL730Ur77C8Z600038fLLL/PLL78AkJmZycKFC8nMzKRPnz5UrVqV559/nuzs\nbD7++OPc44ozatQo5s+fn68751RctWrVok6dOuzfv59p06YVW8bZZ5/Nd999x44dO/jzzz+ZMWNG\nvvj/9re/cdddd3H48GFEhJSUlNzuoQ8//DD3S6pevXqEhYURFuZ8etaEr1QQKW34Yt7777//fqZP\nn05kZCT/+c9/aNasGZ988gmPP/44Z5xxBtHR0cycOTO3C6i4sseMGcP1119PkyZNOHnyJLNmzQKg\nXbt2zJ07l4kTJ3LGGWfwxRdf8Nlnn+X2/ectb+bMmbzzzjvUqVOHW265hZEjR+arY9q0aYwdO5bI\nyEg+/PBDevTowSuvvMLEiROJjIykXbt2vPXWW4DtD//444954403qF+/Ph988AEjRowo8XXp1asX\n4eHh7Ny5k6FDh+befuedd3L06FEaNGhA3759ueiii4p9PS+88EKuueYaunXrxrnnnsuwYcPyHfv2\n229z8uRJOnXqRGRkJFdddVVud9mKFSvo3bs3derUYfjw4Tz33HO0bNmyxJgrQlfLVKocdLXM/OLi\n4hgzZgw33nijr0MJCbpaplJKqTJxJOEbY14zxuw2xqwr4ZjnjDFbjDFrjDFnO1GvUsq3dAZsYHGk\nS8cY0w84ArwtIt2KuH8oMFFELjbG9AZmiUhMMWVpl47yW9qlo3ypsl06jgzLFJHlxpjiB7rCZcDb\n7mN/NsbUNcZEichuJ+pXypPSUlN5c8oUXMVMqVcqUHhrHH5TYEee6+nu2zThK7+WlprK84MG8cjW\nrYQD//J1QEpVgp60VaoEb06ZkpvslQp03mrhpwN552Q3c99WpLyTG2JjY4mNjfVUXEqVaP/GdE32\nyu/Ex8cTHx9f7sc5Ng7fGNMS+ExEuhZx30XAbe6TtjHAs3rSVgWC2KbX8UXGO7lJ3+DddeCVyssv\nxuEbY94FfgTaGWO2G2NuMMbcYoy5GUBEFgKpxpjfgNnA/3OiXqU8KTkZNhyfzpRWbcj0dTBeFmhb\n96mycWqUzrVlOGaiE3Up5S3PPgv/b2Irbhy3hJlTpuDKyIBly3wdllIVpqtlKlWEvXth/nzbyo+K\nasXUuXMB+FcJE43yDt8Ma9qUcdOnE92qVbnqdaIMpYpVll1SvHnBz3cUUqFh+nSR8eML317c+3Nb\nSorc06aNHHEvZ38E5J42bWRbnl2aSuNEGS1btpSZM2dKt27dpF69ejJy5Eg5ceKEiIjMmTNH2rZt\nK/Xr15fLLrtMMjIych/31VdfSYcOHaRevXoyceJEGTBgQO5OTiIir732mnTs2FEiIyNlyJAhkpaW\nlnvfnXfeKQ0bNpQ6depIt27dZMOGDWWOV5VPce8/yrjjlc8TfKGANOErHzt+XKRRI5HExML3Fff+\nnDZ6dG6iljwJe9ro0WWu14kyWrZsKb1795Zdu3bJgQMHpGPHjjJ79mxZunSpNGjQQNasWSMnT56U\nSZMmSf/+/UVEZM+ePXL66afLxx9/LNnZ2fLMM89I1apVcxP+ggUL5Mwzz5RNmzZJTk6OPPbYY9K3\nb18REVm8eLH07NlTDh06JCIiycnJsmvXrjLHq8qnsglfx+ErVcB778FZZ0HnzmV/jCu98PDNcLD9\n/l4sA+COO+4gKiqKevXqMWzYMFavXs0777zD+PHjOeuss6hWrRpPPPEECQkJbN++nUWLFtGlSxcu\nv/xyqlSpwp133plvp6zZs2czefJk2rVrR1hYGPfffz9r1qxhx44dVKtWjcOHD7Nx40ZEhPbt2xMV\nFVWueJX3aMJXqoD334ebbnJfeest+PrrUh8T1rRpoZE8mUBYkyZlrteJMoB8Cbd27docOXKEnTt3\n5tvmLzw8nMjISNLT08nIyCi0dWHe62lpadxxxx252wnWr18fYwzp6enExcUxceJEbrvtNqKiopgw\nYQJHjhwpV7zKezThK5VHdjb88APkzvWbOxdOnCj1ceOmT2dqm7+Gb2YCU9u0Ydz06WWu24kyimKM\noUmTJmzbti33tszMTPbt20fTpk1p3Lgx27dvz/eYvNsWNm/enNmzZ7N//37279/PgQMHOHLkCDEx\ndirNxIkTWblyJRs3bmTTpk089dRTlYpXeY4mfKXyWLUKoqOhfn0gJwd++QV69y71cdGtWjFpyRJm\njh7N1Lg4Zo4ezaQlS8o1wsaJMoozatQo3nzzTdatW8eJEyd44IEHiImJoUWLFlx88cVs3LiRBQsW\nkJOTw6xZs/JtXD5hwgQef/xxNm7cCMCff/7Jhx9+CMDKlSv55ZdfyM7OplatWtSsWdMjW/MpZ+iw\nTKXyiI/P07pPSoKGDaFBgzI9NrrVX8M3K6qyZRS3Pv3AgQOZPn06V1xxBQcPHqRv377MmzcPIHcb\nwEmTJnHDDTcwZswY+vXrl/vY4cOHk5mZyciRI9m+fTt169Zl0KBBXHnllRw6dIi77rqL1NRUatas\nyeDBg7n33nsrHL/yLN3iUKk8Lr4Yxo+HK64AXn0VvvsO3n47935dD1/5kl8sraBUMMjOhuXLoX9/\n9w0JCRBT5JJPSgUk7dJRym31amjRIk8Pzt13uzvzlQoOmvCVcsvXfw/QqZOPIlHKM7RLRym3b78t\nkPCVCjJ60lYpbP99gwawZQuccUbxx+lJW+VLetJWKQesWQPNmpWc7JUKdNqHrxQF+u9FoJjx7NHR\n0cWOdVfK0/Iuj1ER2sJXigL997/9Bn36FHnctm3b/lp9cPx4ZNcuRITXXhOuvdb3q83qJbgveZfH\nqAhN+CrkuVzw/fcFxt+3aFH6A9PT7bHAgAH2S0Mpf6YJX4W8rVuhXj27igJgk3gZ1s8hJgZ++gmA\n1q3h6FHYvdtzcSpVWU5tYj7EGJNsjNlsjLmviPvrGGM+NcasMcasN8aMc6JepZywZg10757nhoSE\nYrt08omJgZ9/BmyXf/futiyl/FWlE74xJgx4ARgMdAZGGWM6FDjsNmCDiJwNxAFPG2P0hLHyC2vW\nwNlnu68cPWo3ss33DVCMXr1g5Uo7phNbhiZ85c+caOH3AraISJqIZAHzgMsKHCPA6e6/Twf2iUi2\nA3UrVWmrV+dJ+Bs32u2uatYs/YEREXYs54YNgC1j9WrPxalUZTmR8JsCO/Jc/919W14vAJ2MMRnA\nWuAOB+pVyhH5Wvg9e9ozuGW1YAG0awdoC1/5P291qwwGVovIQGNMG2CJMaabiBS5F9q0adNy/46N\njSVW57srD9m9G44fLzAop0qVshfQvn3unx06wPbtkJkJ4QU3p1XKQfHx8cTHx5f7cZVeWsEYEwNM\nE5Eh7uv3Y3dQfzLPMZ8DT4jID+7r3wD3icjKIsqTysakVFktXgxPPglLlzpTXs+e8PzzZTvnq5RT\nvLm0wgqgrTEm2hhTHRgJfFrgmDTgQndgUUA7IMWBupWqlHzdOQ7Qbh3lzyrdpSMiOcaYicBX2C+Q\n10QkyRhzi71b5gCPAm8aY9a5H/ZPEdlf2bqVqqw1a2DoUOfK06GZyp/papkqpHXsCPPnQ7du2CUV\nGjeuWAd8Tg5UqcIPP9h9U9zD85XyCl0tU6lSZGZCWppN+gCMHQsrVpS/oG++gWHDAPvFkZiYOzRf\nKb+iCV+FrPXrbbKvVg04cQLWrrVnXcurSxe7xILLxemnQ5MmsHmz4+EqVWma8FXIynfCds0aO57+\ntNPKX1BUlJ2E5c7y2o+v/JUmfBWy8s2wTUiwa+NUVJ6F1HTGrfJXmvBVyMrXwq9swu/TJ3epZB2a\nqfyVJnwVkrKz7cnVs85y3xAVBeedV/ECY2LsGWD+Svg62Ez5Gx2WqUJSUpIdWPPbbw4VKH9tiygC\njRrBqlXQtOCqUkp5gA7LVKoEa9bkad07Ic8+t8ZoP77yT5rwVUjasMGOpvSUzp3tSstK+RNN+Cok\nbdwInTp5rvxOnTThK/+jCV+FJE34KhRpwlch58QJ2LbNvW9JaqpdTMcp69fD/v107GhPDOv4A+VP\nNOGrkLNlC0RHQ40a2AXxFy92rvBHHoEvvyQiAk4/HXbsKP0hSnmLJnwVcvJ151R2wlVBMTG5E7C0\nW0f5G034KuRowlehShO+CjlJSe6Ev38/ZGTYMZROOeccO+bz2DE6dbJ1KeUvNOGrkJPbwv/5Z7sc\ncnk2LS9N7dp2zeXVq7WFr/xOpbc4VCqQZGfb5RTatwdqt4bJk52v5MYbQSQ34edZdUEpn9K1dFRI\n2bQJLroItm71Tn0NG9p9VRo39k59KjR5dS0dY8wQY0yyMWazMea+Yo6JNcasNsYkGmOWOVGvUuXl\n6QlXBWm3jvInlU74xpgw4AVgMNAZGGWM6VDgmLrAf4FLRKQLcFVl61WqIjThq1DmRAu/F7BFRNJE\nJAuYB1xW4JhrgY9EJB1ARPY6UK9S5ebthN+xoyZ85T+cSPhNgbzzCX9335ZXOyDSGLPMGLPCGDPG\ngXqVKjdt4atQ5q1ROlWBc4CBQDjwkzHmJxEpcvuJadOm5f4dGxtLbGysF0JUwS4nx5607dABGD0a\nZsyA5s09U9kff8A779Bp5F2a8JXj4uPjiY+PL/fjKj1KxxgTA0wTkSHu6/cDIiJP5jnmPqCmiDzi\nvv4qsEhEPiqiPB2lozxi61YYOBDSEg/bLakOHIDq1T1T2aFD0KQJsv8AkVHV2LwZzjjDM1Up5c1R\nOiuAtsaYaGNMdWAk8GmBYz4B+hljqhhjagO9AZ2DqLwqtztn5Uq7JZWnkj1AnTrQqhVm3Vrt1lF+\no9IJX0RygInAV8AGYJ6IJBljbjHG3Ow+JhlYDKwDEoA5IqIfAeVVuQnf6fVzitOnDyQkaMJXfsOR\nPnwR+RJoX+C22QWuzwRmOlGfUhWxcSP07w98mgDXXef5CmNiYOlSOvWYqAlf+QVdS0eFjKQk6NhB\nvNfCd6+cqUMzlb/QpRVUSBCBunVhW6oQ+WcqtGrl+QVuXC6YP5+0PiM5r5/h9989W50KXWU9aasJ\nX4WEjAzo3h127/Z+3S6X3f1q1y77r1JO8+paOkr5u+Rk9/h7HwgLs/vnbtrkm/qVOkUTvgoJycnu\nJZF9pH17G4NSvqQJX4UEX7bwwdatCV/5miZ8FRI2bYJOrY/bDnUf0ISv/IEmfBUSkpPh3BUvwj33\neL/yp56iT8o7mvCVz2nCV0EvM9OuZRaxOQF69PB+AOHhNN20lN9+s1ssKuUrmvBV0Nu8Gc48E8IS\nfvLOhKuCYmKoujKBRo1g2zbvV6/UKZrwVdBLToa+LX6H48ehTRvvB9CtG6Sl0aPNQR2aqXxKE74K\nesnJEFvrZ9u69/Ts2qJUrQo9enBBnRXaj698ShO+CnrJydA2fCf4ciOdmBh65iRowlc+5a0dr5Ty\nmeRkCHtjot1zzVcefJBjv9Qi+REfxqBCnq6lo4KaywWnnWZH6Zx2mm9j2bULunaFPXt8G4cKPrqW\njlLA9u3QoIHvkz1AVBRkZcHevb6ORIUqTfgqqPl6SYW8jLGx6Egd5Sua8FVQ86eED7rEgvItTfgq\nqCUnw6CcL/1miuvZzfdpwlc+40jCN8YMMcYkG2M2G2PuK+G4c40xWcaYK5yoV6nS7Fy3hyFvjbSL\n0vvayZNMfKoFqYmZvo5EhahKfwqMMWHAC8BgoDMwyhhT6Ee0+7gZwOLK1qlUWdVJ+pmcHr38I+FX\nr05W+67UWL/S15GoEOXEp6AXsEVE0kQkC5gHXFbEcZOAD4E/HKhTqVIdPAhdMxOo0d8H6+cUo3r/\nGFruSuDECV9HokKREwm/KbAjz/Xf3bflMsY0AYaLyEuAD+a2q1CUlARxNX/C9PGfhF/lvBhiayaw\nZYuvI1GhyFszbZ8F8vbtl5j0p02blvt3bGwssb6cEq8CVvKGHEYeXwG9e/s6lL/ExNAj+y6+2Sh0\n6aJtH1Ux8fHxxMfHl/txlZ5pa4yJAaaJyBD39fsBEZEn8xyTcupPoAGQCdwsIp8WUZ7OtFWOmHL7\nn1y09nH6fPtk6Qd7iwgp0XF8OPp//POJCF9Ho4JEWWfaOtHCXwG0NcZEAzuBkcCovAeISOs8gb0B\nfFZUslfKSatT6tLjLj9K9gDGsPzReFYv8nUgKhRVug9fRHKAicBXwAZgnogkGWNuMcbcXNRDKlun\nUmWRlAQdO/o6isI6drSxKeVtuniaCjppqam8OnkK372fTuzIptz42HSiW7XydVi5NqxP5cruU7iq\nXzpVmjVl3HT/ik8FnrJ26WjCV0ElLTWV5wcN4pGtWwnHniya2qYNk5Ys8Yuk6u/xqcCkq2WqkPTm\nlCm5yRQgHHhk61benDLFl2Hl8vf4VHDThK+Ciis9PTeZnhIOuDIyfBFOIf4enwpumvBVUAlr2pSC\nK9VkAmFNmvginEL8PT4V3LQPXwWVtNRU3oyL4+G0NAz+10euffjKE/SkrQpZu57/L8m3P87X/dtT\nrXkTvxsFk5aaysv/mMKojz9nTVwPBrz2ql/FpwKPJnwVsg5cfwezP2vC/fuLXanb50Tg2doPcNOt\nVTntP//ydTgqwOkoHRWyzPLv2dP+fF+HUSJjIL31+WQt+97XoagQ4q3F05Tyms8Hv0CV6j19HUap\nTvQ8j8TjI/DvryYVTLSFr4LO10f70q5LdV+HUapWZ9Xhg6iJvg5DhRBN+Cro+OsaOgXpmjrK2zTh\nq6AiYpNoh0KbbPofTfjK2zThq6CSkQE1a0L9+r6OpHQtWsCBA3DokK8jUaFCE74KHi5XwHTngN1X\nvX17SE72dSQqVGjCV8HjwQep/tKsgEn4YL+cjr0xD+bM8XUoKgRowlfBY9kyVmadRadOvg6k7Dp3\nhtTdteHDD30digoBmvBVcDh8GBIT+XxvDF27+jqYsuvSBRYePh9++glOnvR1OCrIacJXwWH5cuTc\nc1m1sSZduvg6mLLr2hV+So6Adu1gxQpfh6OCnCZ8FRyWLePQOXHUrAlnnOHrYMouOhoOHoTjfeNg\n2TJfh6OCnCMJ3xgzxBiTbIzZbIwptGKVMeZaY8xa92W5MSaAfnSrgLB1KxujYgOqOwfsSJ3OnWFr\ns1hN+MrjKp3wjTFhwAvAYKAzMMoYU3DaSwrQX0TOAh4FXqlsvUrl89FHfC/nB1zCB9uP/2PtC2H+\nfF+HooKcEy38XsAWEUkTkSxgHnBZ3gNEJEFE/nRfTQCaOlCvUvmsTzQB1X9/SteusCa5JjRo4OtQ\nVJBzIuE3BXbkuf47JSf0vwOLHKhXqXwSEwnYFv769b6OQoUCry6PbIyJA24A+pV03LRp03L/jo2N\nJTY21qNxqcCXnQ2bNhFQY/BP6drVflmJ2HXylSpNfHw88fHx5X5cpXe8MsbEANNEZIj7+v2AiMiT\nBY7rBnwEDBGRrSWUpzteqXJLToZLLoHffvN1JBXTsCGsXg1NtbNTVYA3d7xaAbQ1xkQbY6oDI4FP\nCwTTApvsx5SU7JUqt82bISmJ9esJyP77U0618jl+HP78s9TjlaqISnfpiEiOMWYi8BX2C+Q1EUky\nxtxi75Y5wBQgEnjRGGOALBHpVdm6leKll6BhQxKPdwzI/vtTuna1/fiDv/0XVKsGjzzi65BUEHJk\nHL6IfCki7UXkTBGZ4b5ttjvZIyI3iUh9ETlHRLprsleOWbYMYmNZvz4wT9ie0qWLu4U/YICOx1ce\nozNtVeDavx9SUqBnTxITA79LZ/164LzzYNUqOHrU1yGpIKQJXwWub7+F887jaFY1duyAM8/0dUAV\n16mT3f0qp9ZpcNZZdjE1pRymCV8FLnd3TlKSXXusWjVfB1Rxp58OjRrB1q1AbKx26yiP0ISvAlf3\n7jBsWMD335+SOwFryBBwuXwdjgpCXp14pZSjbrgBgMTXgyPhnxqaOWLq+XD++b4ORwUhbeGrgBfo\nY/BP0SUWlKdpC18VsnEjZGVBvXr2UqeOf0/5D5Yuna5dYepUX0dRMhE7L+zgQXupVctuxK4CgyZ8\nBdgP8jffwGOPwZYtUL8+HDhgRz526wbvvgstW/o6ysJ27rQ7AzZv7utIKq9DB0hPh0OH7Jesv9m8\nGUaNsu+PiAh7+eMP+/548EHthQoE2qWj2LAB+vSBSZNst3hqKqxdC9u32+Rz5ZXQqxd89JGvIy3s\n11+hRw+0ADH4AAAb/ElEQVT//gVSVlWr2uS5erWvIyns//7PThG46Sbbwk9LgzVr7HtlxAgYNw76\n9w/ctYxChSb8ELd8OQwcCOPH2xOGY8fmH94YFgZ33w1ffAH//Cfcfrv9NeBTixfbJRWwCf+cc3wc\nj4POOcc+J8A+z7VrfRqPywU33giPP25/AU6YkP/LtUYN+yWwaZNtGPTvDytX+i5eVTJN+CFswQK4\n4gqYO9d+aKtUKf7Yc8+1E0Dj4+Gdd7wWYtEWLIBjx4C/WvjBokcP+zoDdvKVj1/s2bPtOZ2VK+2v\nj+JUrWobAy+9BBddZL+rlB8SEb+62JCUp73+ukjjxiK//lq+x61YIdIgMkXuv2K0PBwbK9NGj5Zt\nKSmeCbI4HTqIrFolIiJNmoh4u3pPWrPGPj0REYmPFzn3XK/HsC0lRaaNHi3/jImVLjVGy5KvyvcC\nL18u0rChyPz5HgpQFeLOm6Xn17Ic5M1LICb8Ux8QnyXAclq8WKRRI5HNm8v/2G0pKXJ9vTZyxPbs\nyBGQe9q08d5zzsgQiYgQycmRnTtFIiNFXC7vVO0NJ0+K1K4tcuiQiBw/LnLaaSIHD3qt/m0pKXJP\nm8r//65bJ3LGGSLff++hQB0UaJ/fomjCd1hWlm2xjBsncsEFIu3aiTRoIDJkcIrc3NCHCbAM8r6h\n77p4tNSPSJHvvitnIS6XyB9/yKxLL819rpLnOU8bPdojsRfy7rsil10mIiJffCFy4YXeqdabevWS\nv/5/Bg4U+ewzr9U9bfToQv+/mSBPjhhR7m/WRYtsw+K7eP9NqEV9wV1bu41E1E2RDh1EBg0SufFG\nkf/9TyQ729fRFk8TvkMOHRL5z39EoqNF+vUTefll20LesEFkxw6RsTGFPyBeTYClKOoNfVPDAl9I\nLpdITk7RBdx8s0jbtiI1a4pERsrBGjXyPddTl4fj4rzzhG6+WeTZZ0VE5F//ErnvPu9U60233iry\nzDPuK9Oni9xzj9fqfjg2tsj/36NVq4rUqiVy5pkiQ4YUnfyLuO2RqSlyVTX/bRAV9QV3BOSfw0dL\nYqL90vrvf+2XcNu29u8jR3wddWGa8B3w+ee2n/vqq0V+/rmIA44fl69at/ZtAixFUW/oEyCr2rSx\nzZf27UXCw0W+/rroAlavFklKyn2XTxs9WjJ9+QW3c6fIvn0iYhv677/vnWq96dVXRcaMcV/ZssX2\n5XtJcQlw2ujRIocPi2zcKPLtt0U/eOdOkdNPF+nYUeRvfxMZP16Wde0qx/y4QXRfn6K/4B47//x8\nx7lctntq+HCR5s2L/7j4iib8SjhyRGTCBNuqX7ashANzcuSHjh3laMHWEMikIf7xhi6uxfZJ+/a2\n+ZKYKPLnn2Uur6hfDHe39k2LrVkzkd9+83q1HrdqlUinTr6pe2NiilwRVsEWucslcuCAyPr1IgsX\nisyeLfHR0X7bIFq+XGRojUslp0BsJ0H+PWJEsY9btEikaVORu+4SOXbMiwGXQBN+BSUm2l+tY8eW\n7VxZwQSYCTKhYWupH5Ei69Z5Pt7SlNhiq6DccwJxcdIvarQ8+x/vJ/vdu0Xq1QuuE7annDhhe098\n0XXw4osigwf99f9b2T53T7z/nJCQYM/BvfFaitzXqlW5v+D27hW58kqRzp0rNvjBaZrwK+DLL+3I\ngrffdt/gctm+nBtusB35xcibAE99QN59V6RFC5H0dO/EXpyff0qRK6p4rg/1k09EYnpmietIpiPl\nldXChfZ8ZrDq2dO2QL0pJ2OXxLbZXvKv2nIq6hfhnS1924efkmJPJn/66V8xlukLbskSew7JPSTY\n5RJ56SWRqCgp/yAIh3k14QNDgGRgM3BfMcc8B2wB1gBnl1CWJ1+XYuX7jzt8WGT2bJHu3UVatRKZ\nMcM2KcvpscdsEYcPOx9vqXbvlpw335YBA0TuvtO5FltB2dkiz9d7SH6/zrtnT6dPF7n3Xq9W6VW3\n3CIya5Z360y75FZ5/Yx/Ov6rKW9CvbTjaBkYlyKuX1aI7N/vbEVlsH+/nefw/PMVePDOnXakQPPm\nIr17i7zxhkhmpixeXKCh6ANeS/jY2bq/AdFANXdC71DgmKHAF+6/ewMJJZTn2VdG8g9TnHrtaLnl\nphRp186eH5P0dDvOe/hw2+QvbvRKGbhcIuPHi4wctFeyT1a8nHLbvl2kfXuJH/CwxMV5fjjZKw9v\nl0PVIyv0pVhmx47ZQepul18uMm+e56rztTlzbLei12zbJn9Wi5S5z+7xaDUnTthc+dOFD4l06yay\na5dH68tXd2aWjD93rdx1VyULys62Q2Uvvlikfn2RhARJTBRp2VLkrjtSZOq13h+C6s2EHwMsynP9\n/oKtfOBl4Jo815OAqGLK8+TrUuRPzGtqtJHVv+b5j3HwTXjyhEsWNR4nK87+e6W+PMpsyxaRli1l\n661PSVSUd7qUDh4UmV1johy66W7PVfLWWyKjRuVebdHC/QUdpH79VaRLlzw33HvvX30QHnDgqpvk\n2dqTvXISMiVFpEF9l6TfNNVOaNm+3fOVZmXJyjNHyg+Nr3C2AZSaKpJpuzN/+TlFrqrumyGo3kz4\nI4A5ea5fBzxX4JjPgL55rn8NnFNMeZ58XXxyEumPlMPyc/XzZNuQWzx7ljExUaRpUzny9MvSooUd\nVuotD96QLpk1I+xMWE8YN84OghaRPXtE6tb1zvenrxw/bk/cZp46NTJjhsikSZ6pbOtWOVIjUh69\na69nyi/Ce+/Zce3HHptpm8ae/PbOzpbU80bL97UGyYGdnvtG8+UJ6rImfL9cD3/atGm5f8fGxhIb\nG+tMwSdP0mH1asIL3BwOuDIynKmjCGe0Oo2tny/ij4sGE3H9JOq89bzz6/m6XDB2LDLjScZ8PJor\nroCLL3a2ipKMvrcJc9+/nr8/8SRhzz3rfAXLltnlOoEVK+yqkmFBvPRfjRrQqZNdKvm887Abm990\nk0fqyv7XY7xobmPMnfU9Un5RRo6EJUtgwuZ7eHPy6XDhhZCUZHdUcVJODgdH3Mj2X3ZS97vPqdeo\nprPl5+FKTy8yt7Rbuxays+0Kcw6Jj48nPj6+/A8sy7dCSRdsl86Xea6XpUsnGW926ezeLfLAAyKN\nGslvDaJ8NhHk+UcPyvrwXpJ12+2eaekfPSovvihy9tm2hehtQ87eKasne2DFrFPDKtyv2eTJIlOm\nOF+Nv7njDpEnnnBfycoSqVPH/rxx2MJZm+XS/gccL7c0R47YeX9z54rHunVO3niL/Fw7Vl5/wfOj\nyIpq4R8H2Va/gR24P3Wqx/pYKWML34k20gqgrTEm2hhTHRgJfFrgmE+BsQDGmBjgoIjsdqDusjl6\nFI4eJWPuUv522k/cGtGGTPddmcDUNm0YN326x8O47YG6PBm3mA/XnumR8n9ZX4uHH4Z582wL0dsu\n+Xsjnky92vmC4+NtC9f9q+jbb2HAAOer8TcDBtjnCtjWYb9+eW5wzn+/OpMr/17P8XJLEx5u36t3\n3gmJfzq/ZZkIPJk2klcv+4xx/6+24+UXNG76dKa2yZ9bbq7bhqENfmHv/y2CPXtg3TqPx1Gisnwr\nlHbBDsvchB12eb/7tluAm/Mc8wJ2NM9aium/Fw/24Scl2RN9zz5bjnG3HnDwoG3VvPqqs+Xu2mVH\ni/3vf86WWx5799q+dccXd5w50560FdunHR7un+uZOG3PHtuoz8py3/Dvf4vc7eyJ8V277AQ2nwwd\ndnv7bduf7/QozWeesb92M704RaRgbkndmiKPPirSpo3I1q2eq5eQmniVkyPy1VciV1wh8uOPhe5e\nscL2CLhzhs8lJdlxu0Wuz1NWeabxnjwp0r+/yEMPVT62yrriCpFXXvFc+d98I9Knj+fK9zedO9v3\nr4jYoakOn6l++mmR6693tMgKuf12kaFDCwwh3lvxk8hLl9p5NamplQ7NES++aHt11q8v4s59+0RG\nj7ZrJlWwqzc0Ev6ePSJPPWWbB2edZWdPHTqU75BvvrHJ9ZNPyl6sNyxYYFvkFRoB+uKLdiEZ9yJi\nd9xRxIfFRxYssKuKesrUqcG5QmZxbr3VJmVPOessmxx9rVCj5cQJ+7muQCstLc028JYscTbGynr3\nXfslVKhNeviwyHPP2QWUOnSw3RDl/LkT/An/iy/sb9Hrr7evYBHfjP/7n032Tk4Vd9LDD9vkeGzz\ndvsfXhYzZtjZv+7fh08/bdf+8cGkxSKdOGFf862JR+060g6LjbXLKoSKefNELr3UAwU/+6wkLsmQ\nFi38Z3jr7t12wcLZs903JCXZhs0LL5T+YJdLZPJkyfzyO+nRw7YD/dHChfbzUeRHw+WyK5GOGmX7\nRsvxUzn4E/6hQ7kt3KK8/rr9ll+5smzF+UJOjv0ld+X5uySnfQeRRx8t/mCXy4406thR5PffRcQu\n79OmjXfmrZTH7beLzPjnPrsdlYPLWR4/bvvvy7G4Z8DLs8GXc9atE4mKkvtuOywPPuhguQ7YssX+\n8s1N+ikpIq1b5xmuVASXS+Sf/5TsbmfL4HP3yS0enu5SWWXaAnL3bvtTpYyCJ+GvXJnnrFXZPP20\nPUGbnFyuh/lEdrZN+lefnyE5Z7Yr/o396KN2YZ4//hARm+xbt/a/ZC9iZ4m2aiXimvKwnTDlkO+/\nF+nRw7HiAsaZZ4qsXetggSNGSPa/Z0pUlMimTQ6W65BTSX/OHPcNv/9uGzoPP1z4YHdDKLtLNxly\n7l6ZMMF/frGUZM0aux9z7hdbeRQxyCSwE/6RI/bnTI8edhZeGVuJpxrBHTr4ZyIsTm7SP+93yWnT\nVvZPnlx4S7itW0UOHpSjR+049Naty9UA8CqXyy4LsPzzA3atkYquH/v77/lOvjz2mFR+HZQANH58\ngcW+UlIq/jNn9WqRxo3liw8y/frk9+bNtjdn+nTbTSh79tiBGZJ/Laz4rl3lWNt2ctG5fwRMsj9l\nyxbbMHr88XL8IjlyxLZm+/Sx5zeOHhWRQE/4ERG243LhwjL/D2Zn2xUGe/TIbQQHlFPxnx35g+wO\nq5q7qUre9Ti++MIm+hEjfL/scmlmzrSrSsu//iVy3XUVK2T27DxbP4kMHmxPCoeat98WueqqPDdc\ndpldm6AiLrtM5JlnZMSIPC1oP7Vtm12frEMHO/hCpPBaWNkg48Ki5e83pgRUsj8lPd2OxLrnnnIk\n/aws+0EYOtQ2qO66K8ATfjmb5ydO2G0I4+IKDdIJOJMGF70ex/mNRkvbtnYBz0CQO747/U+700RF\nBiGPHGlPxshfE00rMVIvYG3bZvt8cxPCM8/YddnLa+dOkTPPlL07jnpmvoQHuFw2t0VH2610L27n\nnxuqVMa+fSIxMXaz9HL2Xttfe5Mne3WmrfOal33WXWYmXHopnDwJCxfC6ad7MC4viDhR9Hocnepn\nsH49DB7si6jKLyoKzj8fPlpSB379FVq1Kl8BIn/NsMWuKRMdDfW9t9yL34iOtkvMbN7sviEuzq4t\nVF6NGsHGjby3oBYXXwx16zoapkcYA5ddBhs3wtix0NBV9OfDk2theVpkpF1XaMcOuPpqOH68HA9u\n1Qoef7zMh/tnwi+j/fth0CBo3Bg++ABqem5dJK8Ja9o0d2r2KZlA47ObBNzzGzcO3nwTaNGi/IvF\nbdpk14dwf1F8+y307+90hIGjf/88qyp07Wrf/Onp5S+oalXefNP+3wSS2rVh1CiI7l305yOsSRNf\nhOWY006Dzz6zK2hcfDEcPuyZegI24e/cadca6dMHXnvN0YXofKqo9Ti8tdaP0y65BBITITW1Ag9e\ntsy2ZN2WLMl3NeTExcFXX7mvhIXZb4AKrJa4fj3s3g0DBzoantcE0+ejoBo14L33oE0buOAC2LvX\n+TqM2H5zv2GMkdJi2roV/vY3GD8eJk92fqVhX0tLTeXNKVNwZWQQ1qQJ46ZPJ7q8XSJ+4vbbbTfM\n1KnlfGBCgv03Joa9e+2HICPDLrgVivbtg9atbaP+tNOwP2kBrrqqXOXcc4/9JfzYY87H6C3B9Pko\nigg88AB88on9km/WrPTHGGMQkVIzYcAl/HXrYOhQmDIFJkzwYmCqQlatghEj7Jd0RdevnzMHvvkG\n5s93NrZAM3QoXH+9XUu+IrKy7Omx77+HMz2zYKty0MyZ8MILNum3a1fysWVN+AHVpfPjj7bP/j//\n0WQfKLp3tyfSc/ufn34a1qwpVxnz58M11zgfW6C55poKfOl9/31uh/2XX0LbtprsA8U//gEPP2zH\nLaxa5UyZAZPwv/zSnq1/6y398AcSY+xGTS+84L6hWjX7Li6j3bvtm33oUM/EF0iGD4elS+HQoXI8\naOrU3M0Dnn/eY5tmKQ+58Ub72RkyBL77rvLlBUTCnz/f/pT95BP7xFVgueEG+2b97Tfg5pttBl+5\nskyP/egjO2rB6Z3vAlG9ejZ3f1pwe6HiLFsG27fDmDGsWWOHNo4a5dEQlQdccYU9mXvllXYkT2X4\nfcJ/+WW4+247SqNvX19HoyritNPglltsVxw1a9oz7Xn2LS6JdufkV+ZuHRHbup86FapW5emn7Qn0\n6tU9HqLygAsugC++sL/Q5s6tREFlmZ3lzQvuxdNcLrt2SqtWji64qHxk1y67YsYff4hd9rJ5c5GE\nhMIHJia612Sw084jInyzP6+/OnTIzjjOXQ575syiZ6YvWWK3VsvKku3b7cKlB7y/ba1y2IYN9qMz\na1b+2wnkmbYulz1h8d57sHy5HZKnAltUlP1J+uKL2AHHDz5of7YV9M03uZMqPvjAzqL2xf68/ur0\n021rb8EC9w0rVtjXrKC9e2HGDKhalWeftedt63l/21rlsE6d7Hn4F16wP97KO8jSL4dljhsnbNoE\nn39upx2r4JCcbPugU1Ohdi0pegLF5ZfbseXXXkvfvnb4rZ6wzW/+fHjjDTuQgZdfhp9+sqMZinDw\noB2/v3ZtuVYsUX5u9257PrNfP5g1C6pU8cKwTGNMhDHmK2PMJmPMYmNModU5jDHNjDFLjTEbjDHr\njTG3l1buzp228afJPrh06AAxMe7cVFSyd7ns+M3YWH76CbZtgwsv9HaU/u+SS+x573Xr+GtdnWIa\nbrNnw0UXabIPNlFRdqL12rV2jaGyqmyXzv3A1yLSHlgKTC7imGzgbhHpDPQBbjPGdCip0E8/Dd0Z\nlcHuvvvgiSfszNFCfvwR6tfnZIMm3HwzPPOMHcWp8gsPt+tl3Xwz5LRpB9nZdnW5AnbuhGefhXvv\n9UGQyuPq1oXFi+HPP8v+mEp16RhjkoEBIrLbGNMIiBeREpO5MWYB8LyIFNHxWLalFVRgu/deu8bO\n559DlSp57rjnHgAer/80P/xg7w+2ZTOc4nLZCTlXXw0Tk26zSX/27Nz7s7Lsejl/+5vtFlPBKysL\nqlf3wtIKxpj9IhJZ3PUijm8JxANdRORIMcdowg9y2dn2xOPAgfbEU+7aKOnpHDm9KW98N51Va1rR\nsqWvI/VvSUl2DbXPPkll8Yv29Qtr2pRx06fz3AutSE6247YruqSFChxlXVqh1DUmjTFLgKi8NwEC\nPFTE4cVmamPMacCHwB3FJftTpuUZox0bG0use010FRyqVrUnHnv2hJbRqax/dBCPbN1KOHb1wz/r\nJ2BkCRA8C2J5QseOcN21qTw/aBBzjv71+t3zdQILqy1hzdpWmuyDVHx8PPEVWC21si38JCA2T5fO\nMhHpWMRxVYHPgUUiMquUMrWFHyKWL4dbL7yOhBPv5NvUIhOYOXo0Uys1wyQ0PDzqOu6bV/j1mzxk\nNM8t0tcvVHhr8bRPgXHuv68HPinmuNeBjaUlexVa+vWDCzoH3w5G3mR2Ff36RZzQ108VVtmE/yQw\nyBizCbgAmAFgjGlsjPnc/fd5wGhgoDFmtTFmlTFGV8RRAER0DM4djLyluB3S9PVTRfHLiVf+FpPy\nnLRU2wedtw9/aps2TFqyJKg2tfAUff0UBPEGKCr4BPsORp6mr5/ShK+UUiEiKHe8UkopVXGa8JVS\nKkRowldKqRChCV8ppUKEJnyllAoRmvCVUipEaMJXSqkQoQlfKaVChCZ8pZQKEZrwlVIqRGjCV0qp\nEKEJXymlQoQmfKWUChGa8JVSKkRowldKqRChCV8ppUJEpRK+MSbCGPOVMWaTMWaxMaZuCceGufez\n/bQydSqllKqYyrbw7we+FpH2wFJgcgnH3gFsrGR9AS0+Pt7XIXiUPr/Aps8v+FU24V8GvOX++y1g\neFEHGWOaARcBr1ayvoAW7G84fX6BTZ9f8Ktswm8oIrsBRGQX0LCY454B7gV0s1qllPKRqqUdYIxZ\nAkTlvQmbuB8q4vBCCd0YczGwW0TWGGNi3Y9XSinlZUak4o1uY0wSECsiu40xjYBlItKxwDGPA9cB\n2UAt4HTgYxEZW0yZ+itAKaXKSURKbUxXNuE/CewXkSeNMfcBESJyfwnHDwDuEZFLK1ypUkqpCqls\nH/6TwCBjzCbgAmAGgDGmsTHm88oGp5RSyjmVauErpZQKHH4509YYM8kYk2SMWW+MmeHreDzBGHOP\nMcZljIn0dSxOMsb82/1/t8YY85Expo6vY6osY8wQY0yyMWazu+syaBhjmhljlhpjNrg/b7f7OiZP\nCOaJn8aYusaYD9yfuw3GmN7FHet3Cd89kmcY0FVEugIzfRuR89zzEgYBab6OxQO+AjqLyNnAFkqe\njOf3jDFhwAvAYKAzMMoY08G3UTkqG7hbRDoDfYDbguz5nRLMEz9nAQvdA2bOApKKO9DvEj5wKzBD\nRLIBRGSvj+PxhFPzEoKOiHwtIi731QSgmS/jcUAvYIuIpIlIFjAPO+EwKIjILhFZ4/77CDZZNPVt\nVM4K5omf7l/Q54vIGwAiki0ih4o73h8TfjugvzEmwRizzBjT09cBOckYcymwQ0TW+zoWL7gRWOTr\nICqpKbAjz/XfCbKEeIoxpiVwNvCzbyNxXDBP/GwF7DXGvOHusppjjKlV3MGlTrzyhFImc1XFDu+M\nMcacC7wPtPZ+lBVXyvN7ANudk/e+gFLC83tQRD5zH/MgkCUi7/ogRFVOxpjTgA+BO9wt/aAQAhM/\nqwLnALeJyEpjzLPYNc6mFnew14nIoOLuM8ZMAD52H7fCfWKzvojs81qAlVTc8zPGdAFaAmuNMQbb\n3fGrMaaXiPzhxRArpaT/PwBjzDjsT+iBXgnIs9KBFnmuN3PfFjSMMVWxyf7/ROQTX8fjsPOAS40x\nF+Ge+GmMebu4iZ8B6Hdsj8FK9/UPgWIHFvhjl84C3InCGNMOqBZIyb4kIpIoIo1EpLWItML+Z3UP\npGRfGmPMEOzP50tF5ISv43HACqCtMSbaGFMdGAkE20iP14GNIjLL14E4TUQeEJEWItIa+3+3NIiS\nPe61zHa4cyXY+VDFnpz2SQu/FG8Arxtj1gMngKD5zymCEHw/MZ8HqgNL7I8YEkTk//k2pIoTkRxj\nzETs6KMw4DURKXYURKAxxpwHjAbWG2NWY9+TD4jIl76NTJXD7cA7xphqQApwQ3EH6sQrpZQKEf7Y\npaOUUsoDNOErpVSI0ISvlFIhQhO+UkqFCE34SikVIjThK6VUiNCEr5RSIUITvlJKhYj/D9dx35l1\nq8JsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fafea137a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure()\n",
    "plot(points[:,1], interpolated_values,\"-\",label=\"interpolated values\", color='blue')\n",
    "# plot(points[:,1], f(points[:,0], points[:,1]),\":\", color=\"black\",label=\"true values\")\n",
    "\n",
    "\n",
    "\n",
    "nodes = linspace(a[1], b[1], orders[1])\n",
    "plot(nodes, f(0*nodes,nodes),\"o\", color=\"red\", label=\"nodes\" )\n",
    "\n",
    "# compare with linear\n",
    "from interpolation.splines import LinearSpline\n",
    "plot(points[:,1], LinearSpline(a,b,orders,vals)(points), color='red',linestyle='--' )\n",
    "legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi splines\n",
    "\n",
    "Sometimes, it is useful, to interpolate several functions at the same points. In that case, the MultiCubicSpline\n",
    "can be used instead of CubicSpline.\n",
    "\n",
    "Let's define $g(x,y) = exp(-x^2-y^2)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We define the new function and evaluate it on the same set of points\n",
    "g = lambda x,y: exp(-x**2-y**2)  # ne need to vectorize it manually since we use only ufuncs here\n",
    "vals2 = g(x_array, y_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By convention, the values are given as an array whose first index corresponds to the number of the spline.\n",
    "For this reason, we concatenate the two set of values `vals` and `vals2` by creating a new last dimension.\n",
    "Note that in the resulting array `mvals`, vector values are contiguous in memory (if C-ordered)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mvals = concatenate([vals[:,None], vals2[:,None]],axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ORders\n",
      "[15 15]\n",
      "(15, 15, 2)\n"
     ]
    }
   ],
   "source": [
    "from interpolation.splines import CubicSplines\n",
    "interp = CubicSplines(a,b,orders,mvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpolation is made exactly as before but the returned objects have one more dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.98042839,  0.9884649 ])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "interp(x0)  # returns a one dimensional array with two elements instead of a float"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Array of points'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "array([[-0.        , -6.        ],\n",
       "       [-0.        , -5.87878788],\n",
       "       [-0.        , -5.75757576],\n",
       "       [-0.        , -5.63636364],\n",
       "       [-0.        , -5.51515152],\n",
       "       [-0.        , -5.39393939],\n",
       "       [-0.        , -5.27272727],\n",
       "       [-0.        , -5.15151515],\n",
       "       [-0.        , -5.03030303],\n",
       "       [-0.        , -4.90909091]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Array of interpolated values'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "array([[ -1.85261596e-01,   4.69322844e-05],\n",
       "       [ -1.62805645e-01,   4.12435244e-05],\n",
       "       [ -1.40349694e-01,   3.55547643e-05],\n",
       "       [ -1.17893743e-01,   2.98660042e-05],\n",
       "       [ -9.54377916e-02,   2.41772442e-05],\n",
       "       [ -7.29818407e-02,   1.84884841e-05],\n",
       "       [ -5.05258897e-02,   1.27997240e-05],\n",
       "       [ -2.80699387e-02,   7.11096397e-06],\n",
       "       [ -5.61398774e-03,   1.42220390e-06],\n",
       "       [  1.66888656e-02,  -4.19742306e-06]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# to evaluate it many times, it is more efficient to supply an array of points\n",
    "# where each line containt the coordinates of one point\n",
    "vec = linspace(-6.0, 6.0, 100)\n",
    "points = column_stack([ vec*0, vec ])\n",
    "interpolated_values = interp(points)\n",
    "display(\"Array of points\")\n",
    "display((points[:10,:])) # first 10 points\n",
    "display(points.shape)\n",
    "display(\"Array of interpolated values\")\n",
    "display(interpolated_values[:10]) # first 10 interpolated values\n",
    "display(interpolated_values.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Low level interface (TODO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: explain how basis splines are defined on *knots* and how to directly deal with the spline coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(17, 17)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# fit coefficients of the spline (note the size of the matrix ( 15+2 x 15+2 )\n",
    "from interpolation.splines.filter_cubic import filter_coeffs\n",
    "coeffs = filter_coeffs(a,b,orders,vals)\n",
    "display(coeffs.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# cartesian product with 50 points along each dimension\n",
    "fine_grid = mlinspace( [0,0], [3,3], [10,10] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from interpolation.splines.eval_cubic import vec_eval_cubic_spline, vec_eval_cubic_splines_G\n",
    "\n",
    "# evaluate the spline\n",
    "V = vec_eval_cubic_spline( a, b, orders, coeffs, fine_grid );\n",
    "\n",
    "# evaluate with gradient\n",
    "V, dV = vec_eval_cubic_splines_G( a, b, orders, coeffs[...,None], fine_grid );\n",
    "V = V[:,0]\n",
    "dV = dV[:,:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 10)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dV_m_1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.quiver.Quiver at 0x7fafea97fb00>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXmYHUX5/T/vzCST2bKRFbJBwhaQnRBEMRGEAAqIIqKA\ngAKiKJsLIn4FxAX8iYqIggKCiiyKAsqmQIigQAiExSQQCIQsZJKQbbZMZubW74/qzvT09FK93juT\ne56nn7m3u7qq7p3u0+eeeustUUpRRhlllFFG30ZFsTtQRhlllFFGcpTJvIwyyiijH6BM5mWUUUYZ\n/QBlMi+jjDLK6Acok3kZZZRRRj9AmczLKKOMMvoBQslcRKpF5FkReVFEXhGR7/qUu05EFovIfBHZ\nJ/2ullFGGWWU4YeqsAJKqXYRmamUahWRSuBpEXlIKfWcXUZEjgImK6V2FpGDgF8D07PrdhlllFFG\nGU4Y2SxKqVbrZTX6AeCeaXQccLtV9llgiIiMTquTZZRRRhllBMOIzEWkQkReBFYB/1RKzXUV2QFY\n5ni/wtpXRhlllFFGDjBV5gWl1L7AOOAgEZmabbfKKKOMMsqIglDP3Aml1CYReQKYBSxwHFoBjHe8\nH2ft6wERKSeCKaOMMoyhlJIk5w8TURvMiy9VSk1K0l4xEUrmIjIC6FBKbRSRGuAjwI9cxe4Hvgzc\nJSLTgQ1KqUbPCmdG5PMmj30da6DpCdi8yLG9BoVWqKiBqfNh0C7R2vHDisthh8t77lv9JKw6G7a8\n3nN/7aGww71QtZ15/c0x+lRogXWfhbb7eu4feDbIKKj9nlk9Xt+tCZSCwvmgbkCPcx+Jfr7vT7xo\n18utrS1mh4qBq4DLfI7VRK8ukqyy0BByvN7xeuPlMOTy4DJeqEupLzaGGpSx+/T7RDwOwAb0f8oE\nl8HExA0WESaX0FjgNhGpQN+pdymlHhSRcwCllLrJen+0iLwBtABnpNZDz4tkJAz/VM9dqgDty6F1\nEVSsggaLzOMSlh9a5kHbnTD4UFDvB9UJdEJHp37d8VY0MnffTGHk3rUK1nwMOuZB5Xio2hEqd9J/\nq3aElvnmbbu/W+PvainIDKj4HsiQ7t2d5k17w0mCfYnY3WgjMqHb310UUm8imESbCSfrsDIt1t8w\nUg/ri40NhBN6HIFThlFo4ivAfh77b3S9Py/FfkWHVMCgCXpzIugCi0P0dfvrzRQt4UV6IIzcC+th\nuzugagJIde/zO9/sriPqTeH8roK+G5mkNzfsqykxqUPfJ3a7zzFIPW1CD4MJ6beQL6GniMH5NVVU\nxPlx13/gd+E5iaxhRrI23DdAYnLfPbh89Qzvc7Midi84ryojYp8RcrwUif1Qw3IxVXqahA5QmBF8\n3L4+wlR6HyX0bQGS5+IUIqL4uKO9CCMTJYHULZuE58f9OZrkZ2zc7yAVte5GqRC7CWL46BBdbkXx\n0eOWSdNHDyL0v0riAVARUdcZlv0qyQdci4niKvO0nsx5PRRMlHwUpK7aY5yXl2JP1YKxYRNkXyD1\nErNdggjbxEc3IXQTlW7fuxmqdNOxWT+IyCzgZ+gxw5uVUle7jg8FbgEmo//RZyqlFjiOVwDPA8uV\nUsda+74LnAWstopdqpR6OEk/+4fNEvVCSJv80yL5NMm9VIk9sgVjgr5G6iVgu4QRdikOjBYBFhFf\nDxwGrATmish9SqlFjmKXAi8qpU4QkV2BXwKHO46fjw7ldtv31yqlrk2rr9tm1sShPlvaaHBtUVHn\n2qKg3rFleY6NOJ+xipTlREwrI3fEeOhYQVPGaCL8ARv28G42KGMiOEwf9Bn9wt7OcPPBNGCxUmqp\nUqoDuBOdvsSJqcDjAEqp14BJIjISQETGAUcDv/WoO1VLZ9skcz9kTfJJiB36BrEnIfVUiL0fEzpE\n/zWTlNBNyvQBQk8Ad6qS5fROVfIScAKAiEwDJqAnTgL8FPg6vfNZAZxnZZn9rYgzxjce8ifz+ohb\nKSALki+Was+L2ON+rlRIvYa+QeptxFbpUWBC6CYqPQimhG5C6qVH6GH4ETBMRF5AT558EegSkWOA\nRqXUfLQKdyrxG4CdlFL7oHNeJbZbSt8zT4vQs5iI4EfocS7G2BN4LMTx2+PEo0f12O3PFTe8MZGv\nXkPZR7eQl4+e9sBoCvBrar61hWAFWmnb6JWqRCnVBJxpvxeRJcAS4NPAsSJyNPof3CAityulTlNK\nrXFU8RvggfCuBCP/0MRTSzA9Sx4zzpJenHEjZqIOomYd7pgkvDM2sfcFQodtLnwxrB9PpBOa+IRh\n2Zn0Dk201nB4DT0A+i7wHHCyUmqho8wQoFUp1SEiZwGHKKVOd9XzIeBiRzTLGKXUKuv1hcCBSqnP\nxPmMNkpfmecBvwszTZJ3qvikyj0KITpvqihqHeIp9rBzGohP6LHVel+JdtnGwhdNI12KCKVUl4ic\nBzxKd2jiQmc6E2B3dMqTAvA/4PMGVV9jrchWAN4Gzkna1/yV+blWe305/0Lafc9bteeh1k3OSToJ\nK5ZSL3VCt5FDsi4TIg1T4Fkq9JSU+TzDsvtTnjQUD1kMbub1gAhTLFHh9t6jkntUvz0PtW6i1OP6\n6TZiKfV+7qOD+V1dKj66/f8vcZVe6uhfNkuUB0RWxO/Vh6htpWXJmJCkfaOZqvWoxF6SpF62XbbC\nhEhLaWC0DF/0LzKPgrTVdZS2otSfRLVH8dmjqnWIRuz1BmXSIPWyStfIwkc3uWf8yqQ9YzQCtpXn\nQ3nSkBeyjntPUl/cWPcocd9JJiYlLQPJJlZFlid9IR4dcps1GgRnPHrHAih49CnPCUZl9ED+yjwq\nSUDy7IJpIu3IlySqHaJbMnmp9aDPkXXkS9l26Yk0bZfCMmj8LlRt1nn1vZCmj54CyvnMSwlxHgBe\nyPKhkBbJJ0l8FZfYs/LWTawXQsolsV7Ktks3ktouhXWw+UfQfh3QBYMW+JxoIS0fvQxj9A0yTwtB\nF09WRJ9EeafltYcRe1Zq3VSBZ0nqZULvRpxol/o2TeCbfwT20sjVX4LKnbMdGO18z7CTZdjIP878\nazm1t2oevHwjzPw5DEjoi5oSfUczDIhghCcdaI16fpTB0yikafr9pJHUyUYcpR45Lr3USF0B/wSO\ncO3PKCZddUHhElC/oscI5pA3oGJMd7msZow+n06ceavhL4Dalr4dZ96/BkALXbD4b3DnofCHA2D1\nS7D00eT1upNauRNcta6Gxz4Hb/5Fv29eBu+9HF5vksFVVeh9btsKaF/tf06UgdMoScAq34BaFf6z\nOe1B0rp3DQo6kMrg6LV4J8DLGl3AhcAcj2MZJeuSSqj4JnpNdwuDLu5J5JBNoq6mJw06WIYT+ZO5\nO1tglM0PW5rghevg5l3gvo/D8n/r/aueg79/GrL69aEK8MaNcNdu8PrtsPTPcOeu8PsJ8NT50esz\nJfa2VfDUZ/XrzjZ46054dhY8MgHeuc6MDOMQe+c6/ZltND8Lb5wAr+4Cm6xFUkwiYeKS+mYrsVxh\nJbR8FjbuAnXLo0W+RM7K6MzAeAtwGXAk+RN6JfBReqfStpFVtMtwkJkg3wJGajL3Q5qEPmiPsI6V\n4ULf8sy9btqWNfDs9+Cdp6HFpdSqauFjd2kCksp0+9L4Ivzzi/qBYWPJ32HkXrD/RbDjkdEHDZ3w\nGwhd9QT8+2SoqIJnzoG374KOjVBVBzudBuM+2vN8E9vCJvQgG2bzUnj9aJj2Cqx7CJZcDc1PglTB\ndqdC9eSe5U0+e5TIl6Zroc1KC912uT5p4CnAQH08auRL5IiXp9ALxgC8DrwM7B2hwTRweHiROAga\nHJUKqLxJCyI1HZoH5zPBaMCIwC5HQY1ppvBSipqLgfw98ysybK/QBRvehtX/09ua/8GUWbD3Kfp4\nGuFO7Zvg6e/Ai9f3VKkADRPglGehboz3uRD/glEFeOGHMPf/erY7eiZM/hxM+ES4Xx/FY3cSe/NL\n8NJRsGUV1L0PWl6GyjoYexaMvxA6JvhWsxVJfPWW38O607rfV74Pan4JAz7oXUfqfvr/gKOAA62/\nswhcm6ZoSDg2FEXaJc28aOKhp+SZq+0Ny67s2555/yLzqIhz07e9By2N2t7obIOuzdZr6+/QybDD\nIeb1mZBc21p47FRY5lrvddzhcOjfobI60kcAzIl99WPwzMehy/qyKgfD+Ith3JdhgIvQTL7PqKTe\n9g9YexzaMwaogOG3QO1p0BJy30X9//oS+pvAeLb+CtiKUhogTWnyUykR+sIymUdB37JZ0obXxRhG\nADXb6S0thIX6vfcKzDlXR8psPwOqh/XctvwPRu6ny6a9yMSyO2De6aA6uvcVWqBqCFQN713eJMTR\nNLSxHmj9D7x3IluJvGIUVM8A1QaqCeoHB/c/ajijr/Uy2b3DQl+ZbBQBUcIX00gBkMdqYttIPPu2\nrcxNUKypxUn9uyRpa5WCxf8PXrsS6iZD3c5Qb23267ZRIAYiJq5a3/wqrPgEVO8FtTP0IFzV7sFt\nhn3m1JS6H4pF6hmlJDCVelml0k1Lme9sWHZx31bmZTKPg7wJPm9i79oC69bCoLFmhG0Svx6V1Dsb\noXKkHoBzI2nURL8j9QzzyxTTdkmLzPcxLDu/b5N5/4ozzwtRwibTQNzFm21EjWWvHAgjt4cGSS/M\n0eR7cn6+qtHeRA7hnyWsz1H/Z7HCGfsJ0k7UleR4kSAis0RkkYi8LiLf9Dg+WETuF5H5IvKKiJzu\nOHaziDSKyMuuc4aJyKMi8pqIPGItPZcIoWQuIuNE5HER+Z/V0a96lPmQiGwQkRes7bKkHetzKBa5\nR0WcSUpRzjEl9aDvyPTzBfXJpL+ZkrozPj1L5NDGNkzoIlIBXI+eXLAHcLKI7OYq9mXgf0r/BpgJ\n/ERE7CvlVutcNy4B/qWU2hV4HPhW0r6aXJqdwEVKqfkiUg/ME5FHlVKLXOXm2IuVlkG8wdU4iJPZ\n0EaSFYTCzjGJXTcZoDSNV/frSxZ5XyLFqPeTQdKomRdLaWA02QDoNGCxUmopgIjciZ655eQ/Rfcn\nbgDeU0p1AiilnhKRiR71Hgd8yHp9GzAbTfCxEarMlVKrlFLzrdfNwEJgB4+ifdZryg1Z2zPFUuxB\nMJllGtWCidOXLPKoF916ydnO2TYV+g7AMsf75fTmv+uBqSKyEniJ7tllQRillGoEzbHAqKQdjXQ5\nisgkYB/gWY/DB4vIfGAF8HWlVEiOzDKA7BS8m/iiqPa4S8MFlc9LqQcp8Sh51MP6YWNbU+l5KvSM\nMXuj3lLAkcCLSqkPi8hk4J8ispclfk2RODLEmMwti+XPwPkenZwHTFBKtYrIUcDfgF2Sdm6bRdQF\nmk0Q145Je81Pk9S8UUgdvD9PUuvFtB82IqXaTYPUizjImmcsehrwqX9GPcxwaOwrlnkWWwE4pziP\ns/Y5cQbwQwCl1Jsi8hawG/B8QK8aRWS0UqpRRMYAARnyzGBE5paZ/2fg90qp+9zHneSulHpIRG4Q\nkeFKqXW9Knvi8u7Xk2bAjjOi9nnbQ9rqPa5qj0LsJoQZptZN86zX4U/oQX0wWUDD2Y+w77zfrnDk\nA1OVbkLoAFWzoX12wk6ljrnAFMv3fhf4NHCyq8xSdOKcp0VkNFrILnEcF3rb0PcDpwNXA58DevFq\nVBjFmYvI7cBapdRFPsdH2/6PiEwD7lZKTfIo1z/izEsVaSj4ODHtaeUpTyNe3a//aXq0Jt9zprHp\nJRb6mFUs+rKU4sy9Ykm8yj7iHWcuIrOAn6PHGG9WSv1IRM4BlFLqJhEZC/yO7lzBP1RK/ck69w5g\nBjqZTyPwXaXUrSIyHLgbnSdiKfAppZTJHeDf/zAyF5FD0EmUX0H7Ogq4FJjo+DBfBs4FOtBX5YVK\nqV6+epnMc0ZScs+C2PMg9aB+p0XqmRA6mJF6iZE5ZEPoaZH5CYZl7+3bk4bKM0BNEdfXK51R+XzJ\nPQ21XiZ1H5QgmUP6hF4m80jo+4m2ijgaboSw/uVJ9kkHVtNe/9N0sDRuBExQ5EuYXx7FTzfx0iMT\nup+fXqJEDulHuqSFbSTRVumSeamTdFow+ZxZEX6UhZzdMB1ETZPUIV4ETJIB0qDjJm3biDw4aqOv\nLDBtIU1CLyMSikPm2wpRp4U8CD8t1R5G6pCMPE0iYLJQ6WH9CmrbiUQqPU8k+EVQJvSiIH8yLxN5\nNvD6XpMQfFxyT2O5uKSknoVKt8ukpdJj+ehZI+gXQBt9ltATp7DqGyhdm6WM5EiT4KOSe6mQehyV\nnrRPQW3biG27pImo9k2OhK62mLdTBlAm820PboLPmtxNBk3TIvU4Kt2vX2kMkGZmuyRBUv89Y0IH\n/b11PBCpV2WUybx/YfMmGFADlQPMz0mb3LPKq+I87lcmb+slrYyMRipdET+XXRoDqMqqp5ZYhA7R\nSL3eMJ7QBNtINEv+i1Ncux88/F14Zy4UCuHlk6KpEX56APx4T/jbBdCedNmeENhZAuf9AjqW9Fy4\nIcoWBYUCzL1Nf7dNjcn6X+/avKAUvPR7q+0uWHAvdFo/i/PIgAhQp2Dzc97Hgr5Dv74F9cnZH9UO\nHa/3Ph4GE0/Yl+wagRsxD/Rvc21JUAC+AXwMcFofMeqN8gtk053R69/Gkb8yX/EiNK2Cpnf1+wkH\nZtte7XCoGgTTz4bpZ5ktgxYHbvKY81sYsSOM3Cmd+rywAVj2PPz1K7D0Gb3v5T/DoRfEa9MLbvVZ\n6IIHvwJvPQEIzPk+rF0Ex98C+57RfZ6pUo/jqW9eC/85A1Y+CLP+C4OmeZ8f13rx69OgtbD049D5\nBox5WS9rF9ZXd7uxbJfvATcBjwL3+pyYVfhiBbAn8BV6X5R2mxnYLvJJ4DPm9ZZRhBmgF8yFHfaD\nihx/FHS0afuhv2HjSpj7O9i8UVss7U2wy0fgwM9l096GLfDX0+DVu7r31Y2C938NDjwXqgMkatyc\nKqDXJG0bqF+vmg1PfRbaVsJuF8B+P4LKan3Mj0iDJh2FzSBtfgTqPgJbFsOyY6DjTRjyA2i4xF8Y\nmFhVxjNH3wB2t94MBe5B53TyQjHj0SPeX0bL4qQ0A/RKw7L/17dngPaf6fzlkMds0d4Ct30CXnuk\ne9+wneDsufrXjwniEHpXBzx+Ghz2e3j+Sph3FVQPh/f/DsZ9tHf5ICL1I3W/fnWug1d2hxGXwpor\nQLXB9rfD4BPD20qN0D8N3IVemOYHwNSQE/oRoZfJPBL6zwBo1IG7Mvmbo3U9/PYYWPrf7n11I2HI\neFh4L+z/BbN6wmwGrwHSZy+FN+6EpqXQ+F/YfgYc/Aeo9VrsiuAok6i2y7KLoGs1NF4AlaNgwkNQ\nc5BZW2lEu6jn0amznwIOCanMRjFnjGYY6ZIE28gAaP8h86iIuoLOtor2Fnj0Chi3PxxwKoyeqrd6\nyy+O+hCN4qW/dT+89P/0vsb/wn7fhgOvgIrK4LaD/Gs/Qrf7Zvdr46Pw3m3dxwaOhY6VvbkqU0If\nC5VztJ0TKXyxTOjbIvL/Gk1+WpbK9N4y4UN1HRz/M//jpgmp3AhTpZ1vwxMu7/+NO2Hnz8DwqWZt\n+x0PC2Hc0AxLz+neVzkMhhwF2x2okzybtmMfI6Sfft+FOH59lOysUS+UCb0YKM2vMM1V7LN+MJjk\n9OjvSELo0Pv/3dkOD3wK2i22HXOIJvHJJ0LNyJ5lk6wm5KfS11wGW96GQbvCqAtgu1Oh0vqtPpBs\nJhmlnqyr2Am6SojQyzZLP4HfzZGH+nfevP2d2E2nuXvBrUyf/Dp0tMAHvg+7nQxDd9T7w8IY07Bd\nNj4Lra/BXg/C8COh2SPqym/yUxrL1KU6a7RM6GnAWmnoZ3SvNHS1R5kZwE+BAcAapdRMEdkFPXpt\nz/jaCfiOUuo6EfkucBbda39eqpR6OFE/c49mOVeV9pMyL4unv5N73Jmk7zXDhjdh5F7eoX9h82aS\nLDqxAehsgirXRRBnAYxcIl0MygDFT6EbM8olrWiWXxuW/WLvaBYRqQBeBw4DVqLXBP20UmqRo8wQ\n4D/AEUqpFSIyQim11qOe5cA0pdRyi8yblFLXxv90PVGc52CSSZhZPwiCbqI0ib6/q/a4Sn27ehi4\nt//xNFICBNouHv/kuKkAygrdQpEVejLOmAYsVkotBRCRO9FxooscZT4D/EUptQLATeQWDgfeVEot\nd+xLNQwy/+n8SdESY0sLTT5bUjQ7tv4Gk6n5bpimBAhrN+iY3/GgmbdBqQDS7ENYe04Yk16xJ81F\nfJiUzmDvDsAyx/vl1j4ndgGGi8gTIjJXRE71qOck4E+ufeeJyHwR+a2l7hOhBB2qDOBF6GkqfC9C\nj6vi3YTeX1R7nEHSOHHpUdqMG+2SlkI3Od6n86O7EWP6f4aYvVBvKaAK2A/4MPpK+K+I/Fcp9QaA\niAwAjgUucZxzA3ClUkqJyFXAtcDnk3YiX4Td0HmRV14En9Sa6U92TBzrJWmOlyyiXfz6VMyBUSNC\nL7bdYiOi7ZIUPpp3xnS92bjir57FVgATHO/HWfucWA6sVUptBjaLyBxgb3QuBoCjgHlKqTX2Cc7X\nwG+AxDl/S89maY6xpYUsLJq0rBjoP3ZMXOslCP3BdglD2HdgJM1KQxWXxkPFCHOBKSIyUUQGAp8G\n7neVuQ/4gIhUikgtcBDg1Pwn47JYRGSM4+0JwKtJO9o/bJYsJ/e4CT2uek9LqdvoD4o9qlIPU+lZ\nDY7mZbuUFXo2SPCLWynVJSLnoVNW2qGJC0XkHH1Y3aSUWiQijwAvA13ATUqpBQAWuR8OnO2q+hoR\n2QedY/ht4BwSIv/QxI9b7RWbgNJsP87FklUIZLG/1ySI8osjSRbGsLbSTNYVJ3TR5HjY5zfy0EuB\n0MGf0FMKTXzMsOxh5URb8RDlxs2CoNIcaIzjv6et1G0003cJPYpSNxkczWOSUVh/shwYLSv0Mhzo\nGzZLHhkR044isW/gYpB6XyZ0MI98KTXbJe+B0VQIvVSQIaGX8iTFFFF6A6BpII2B0rQGWU0HUtMc\nKIVta5C01AZH8xwYTTwoWkqKuBR+JfRd5O+ZzwxoL+ral0mQRLkmOdd9Qy//N6x+CZ2+Qen1NVEw\n0Pq72/EwPObSc+Dd1/YWeO8NaFoNOx+W76pPcWH6cEripcf1sfPw0TP30EuJSO0HTEqe+TOGZaeX\nPfP0ELS8lxtJiT/IVtm8Btb8BxqmQP1OUOVSL0ksGbf9st1UeOTzsH5xz3J1Y+DYG5MReVcHzP0H\nNC+GtYthzWJY8zpsWgmDx8JZD/UNIgdzPz2Jl57Edskj0iWob4ktl1LxzyF1y2UbsVlCyVxExgG3\nA6PRYTS/UUpd51HuOnRwfAtwulJqfsp97QlT4jclfedNUjcCFt8EKx7U72u2h4bJequfDKM/CKM/\n5H0uhJO7KsBb8+Dt++Gd+3sT+dTTYOZPoWa4vkHj+umVA6CzDR66TK+jaWPUbnDWwzB8YsyKXXB+\nx1EeyHGQxsLJWQyOBvnoeQ2M9jtCLyMKjJZVBS5SSs0XkXpgnog86soadhQwWSm1s4gcBPwamO5Z\nW9QBP1WATc/C4GkglYYnOeC+udzkXuiAliXQtAiaF+m/TYugaUF3mbaVemtZBntMgxEHEQivGPDO\nNlj+mCbwpQ9A6yq9v2FHTd4Lbof67eEjN8Jk19qWcQZJCwVY/CA8d0NPIp90CHz+fvN1O8Pg/j7z\nInYTpQrFGRwt5sBovyH0LeFFyuiBUDJXSq0CVlmvm0VkITrRjDNr2HFo9Y5S6lkRGSIio5VSjb4V\nmw72NVTAil/Cyx+F7Y6GER/TeaarIual6dyo81SvWgQtC6F1EWxeBC1vgHJc3VIJtTvB0P1g7Wy9\nr3432PVbMO5kGDzAvM22RnjjH/Du/bD8UU3oCIw+CN73VZh0LAybCu/Oga4K+Mi1MGiYf30mpN7V\nAa/eCU9dDav/BwPrYe/T4KXbYfcT4HN/gAEp/YS1ibvQBRtWQtNaaF4LQ8bAuPcFT7ZJis52qKpO\nHvGRVKUX03bJjNCLjQI61XdKKNssvSEik4B9gGddh9yZxVZY+3qT+YL9ojSp0bler5Te+Ae9SRXU\nfwiGfgyGfBRGTvY/d8tqmLsPbHm35/7KBqjdDUadrP/W7Qajd4O6yVBZDcv+CB0bYddvw/YfB7G8\nZdOZl2/dAU+dAiiorIGxR8CUY2HiMVA7umfZ0dNh+w/puWMtJAtnvPdUePUuqBsFh30fDjwXWtdC\ndQMc9XNor9Tp85OihyIXuO5YWDYfDjkdTrnBu2xapK4U3PV5mPkN2H6v5H5yFEJvXwfrX4IxM3sf\ncyJtQu9YCFIHVRN6H/NCIkK31fl8YDOwL1AdUFnaWAfsnGN7/QPGZG5ZLH8GzldKxQ98c5JqZb0m\n1ciQ7k0kZIXzEVC/L9R8UhN2rbUNHNt78YMC3XWNmAnjPuO9QIKNIK98xMEw5Qsw/lgYc1j3IGrB\no55K142SJEb9gC/CxA/Bvqd3K/CBDXD0L7o/S9px6BUVMPNLgIJDzwr+ztLAkn/DC3/U29mPwK5H\n6P1h5AbJbJfN6+Afh8OGBXD8G1A3LrjdtAh9UCOsPEqLmDELQRxP4ySEHoga4B/Aj9CpQ64GxgSe\nkRxzrC1l9OU5FxFgFJooIlXA34GHlFI/9zj+a+AJpdRd1vtFwIfcNouIKHaPEQrZeApsehiGHA1D\nPgZDjoTKwdHrcSPOcyRKFI3pRWRSzvSnYpzPlORit78PpbpJPOtB0NtOhJf/rF+P2QO+8BAMG999\nPGkYn5NQN70N9eNhy0Z44HBY+yIc9GvYxSOVRpqhi2oLyEAotMHSmbD5ORjxV6g5zrzdsHYgxG6Z\njv4RvifwF2B8UOEMUJtOaOLq8HIAMqpvhyaakvnt6BSPF/kcPxr4slLqGBGZDvxMKdVrADQWmauC\nvpAHHRhvADSqX5YlwZuQpimxmnyuvIk9D6xfBj/YUV8XH/gqHPV9qPb4MtIi9Ec+CbucAs9fqYn8\n0F/DHuf4158WoS/9ITQcC2uugKZ7YNRPYLuLorcb1IYNT0JfC4wCjgBuA1IQT5FRJvMoCCVzETkE\n/dvnFbatnPpgAAAgAElEQVTObOFSYCJW1jCr3PXALPRtcIZS6gWPuhTjI5B5lsSSNcmHEXyexB43\ntLEUif3Bb8NL98BJt8BOHwgum5TQlzwL907XIkJ1dRO5Sf1ex0wJvbAFXtkRVAV0Loeh58CYX/W0\nyUzb9GvDjV6Efgf6tr8e6Ag5OSukQ+bNm83mUtQPKvRpMjeJZnkaCJXESqnzUumRE0mmpIcRUdTk\nWFFXE3LeuF7EbjKQahq/7vwsXp/B2fcoxG63Xyqk3tWhB3Ivng8Da8PLJxkYVQqe+6b1ugtqxwAK\nOlpgQF13/VGiTUxj0dffAx0r9WsZBDXToWstVI0MbjdV/3x/dBpuQdNEKYQrlhGE/rvSUNjECy/4\nDX75kbzfzeEmzLBYd9MIGZNyYYOncWLWSyV3euUAOOyS8HJuxBkYfethWP5k9/tBI7R3XuV6iKQV\numj3Y5OCxp9171dboP0VGHyiWd1xCb1XdMuurgKlEn9ehh/yz80yLMP28sy3ktaAZB52zLZkwwTB\n1HYpdMHt+8LaV6BhPBzyPZh6CrQF/ECNk1fFy3bZ8BS8+EH9umZPmHgz1E/T773ERlTLpU/lcEnH\nZlmjzC7UkdLcv22W1GGqZuMgyYpDfirHD6YqPsziSNuO8Spjqtahbyp2U5hOMFp4BzSvgBk/gX2+\nBFWD9PGwvC5pKPTlP9Ohh2O/A2O+CRUDu495tW+aLthGv5khmh9EZBbwM7pXGrradfxD6KXjlli7\n7lVKXeU4XgE8DyxXSh1r7RsG3IUee3wb+JRSamOSfpZOoq0k6V/jEpCNKAQfVB6Cl5nLk9jjkDrE\nz6/eV3Koh5FfXUHP1P3CmzDI45+QJaG3vQVbGuHA+VA31TwW3avuzOLPYVsidIuIrwcOA1YCc0Xk\nPmc6EwtzbKL2wPnAAnqGBF0C/EspdY2IfBP4lrUvNvInc+dTP63WjVMD+OyPqsqjlPcj0DSJ3av9\ntEgdzIm91AZM/RA0MFpRAR84O1nmxbiEXjEQ9p3dHYKblHRT88+90HcIvRWDwXLA58uaBixWSi0F\nEJE70elL3GTuac9YiQqPBr4POEO7jwPsbH23AbNJSObFzX/amXCLiiaPzQ9RF7cIK9/i2IL65YUN\nji2sfdP9zj6FIerCGUmikPJE0EOngeCHWNCDMM5iE0OB6h16z6Xw6oNX2371xn2wGgmtGkprcYtM\n4E5Vstza58bBIjJfRP4hIlMd+38KfB0d0u3EKHtSpZX/alTSjpaOzRIHUQjd75NGCTmMarv4KdWg\nUMK0FLuXUo8bAWPaN5N+lBqSJOrKw0P360MU/9xvf2oJufqOSnfiv7O38MzsVLIzzgMmKKVarQyy\nfwN2EZGPAo1WxtkZ+Kh3C4kjQ/KPZqE1YS0pKoEoj7KoHnJUpZQ04sQvKsarvaTRL06YfC+lTuiQ\nbgoA03rTmC2adYRLpF/AaRN6OtEsC5VZ3v7dZWmv9qwZ7ZcrpWZZ7y9BT5a82qsOq8wS4ADga8Ap\n6G+xBn233KuUOs3KPjtDKdUoImPQ6VB2j/4Ju9FHlplxoi3G5oMo1k0Uiwai2x5BNoy7fS/4WTFe\n7fn1oaMV3v47/PNrsPpdjwIB/QqCiU1VbIQ9cOKuMxrXcjHtQxqWS9Bni/TbvV9aLnOBKSIyUUQG\nAp8G7ncWEJHRjtfTgAql1Dql1KVKqQlKqZ2s8x5XSp1mFb0fON16/Tl0NEwiFMFm2RTjnKR5IbwI\n3efC8yL0JBaNXwRKUGSK6YxOv5twA94Tk7ysl6a3YO0/YOk/YOUTeiGLWfdC3Vhz+8WkT359KCVk\nlUo3qN68LZc4iJT/vPQslzbjAdDeUEp1ich5wKN0hyYuFJFz6E5n8kkRORed96ANOMmg6quBu0Xk\nTGAp8KnYnbRQBJtlVW7tdcP0YRBBWZg+Br3ILc5gVVwbxkvl1QMbX4N/n6Rzcztx8I9hn6/F60NY\nXwDWLIR5v4Xpn4VxMXLb54UktkucRaPztFwym1DkRBqEno7N8oKhe7GfLOzTk4aKYLNsirBl1aYf\nMrBovGwIv4gX536loGsTtC+Clidg4+p4Noyf9VK5K+z42Z77p3wB9r7YpwHCrSB3X2x0tMH82+Hm\nD8L1U6HQqVdyKmUksV3iRLrkabnEsVugbLmUOIqgzF/Prb1o9kyUsoYXqt/F73XTVNwDrX+BwrvQ\ntVJvqhWogFE/huEXei/64EccXm3YxLC5EV7+Kqy4GwY0QEcTjPkwHPYwVFgLH6SRsXHNK7DoN/DS\n72Gz9UQZuTucM6/n0nWlbL2UFXpv5KbQ01HmT6n9jcp+QOb1aWVeBM887EpJY16/DT8V7kXcXmX9\nCN7Qg7cveve37OUvdx0J7ZdBwfGwk6Ew4i6oPkLfoEETgMB7tqmzjfVKL7v3xgXQsQ4mfgF2+Qb8\n92Nw6J+7iRzMQwv9fPVCJyyfA6/c3U3kFVVwgscapKXspecduujnr2fhocf11fu4h95fUYJx5onm\nGWP2MHATtx9pJyF4B2H5zXptAtQ7MOBGaP8tOLPoV+wOo+6HAVO698XJwWJ/nQPegdfOgXUPw6Cd\nYOrdsNNh2s55/4PQMUwP30SZSepu29luRRVUDoT29d37DrkStvexV0o5Lj3JwGiaseh+hJ4EXu2Y\nzDwtE3rJoQg2yxM+R/NaycSP7BXwV+v1GGAyMJZwdjHpt1u1F6DiUVC/AvV3/Z59oOJLULgYBsyE\nut+DWHXHCW+ziVUVYM2vYPklUGiF8RfAjldCpYN5/QZJo7TnRHU7PP4VePk3MGJPGDMN1r8OJ82G\nisrw520pErqNICWbtuXid8yL0Itht0DGlks6NstDaoZR2aNkdtlmSQdJBzxNHwbuK9RmFgFGA2e4\njtdb+w8Hvgm9wpxM1Lt9EbcAtwI3QuFNYCDwWaj8EnAQoEBWQOH/QBxj00GK3O9YC9D+Gqz+AjQ/\nBYP2gEk3Q/1BvZcaMQ1lDNq/9fhy+MsnYfWzsOtJcOTNsOYlHepYYTUcFsYYR6UvfxEqq2Ds+yKc\nFANx851krdCzsFtMFTr065mifQVFUObOeHuvu/k5YBCwBwYLHEWEk2SbgMXAa+hB2dfQGSzdS2Tt\nik56dhi9Z+OaPEAa0PMObkIvitsOTAK+hH5wjOguGmXANDCMsQPe+3+w9gqtzEd8GyZ8q2c6Vb96\nk6j0lU/Co5+CzWth+jWw90VQHyJ0kqr095bAQ9+B5c/D11/VC1jkgW1FoUdxPY1Vuimhl5V5FBRB\nmTsfHl6qdgCaPBuA/dBJy/bDO2zCxB9XwGrgDWt70/rrnOFYgSbYw4BX0bl0JgJfRCvyCrzTKgSp\n8ha0bXMbevlUAY4CznbU6UKUAVM/Rb7lRVh1JnTMh0HTYOzNMGjP7vsnbJDUJggnqQflewGoU/DK\ndfCfi6F6GHz0nzDuw/qY20t3I65KVwoe+wE8eoVeTu60u/Mh8qXPwMTp245Cj5K5MVIuF8hLpbdu\nI2GSRVDmcVCJVuoHo1cLj3LTrqN71izoC2mytU2xth2BarR3fT7wSTTxBj3rwh4kdwL/h1ben0Wn\naBhvHfNS9D4XnFcX/JquB9Z8DNofgyFXQf350ODx68b0uWia76V9Hfx9D6gbB0f+BRomeJ+XNP+M\nu921b8BP9obRe8D5z3qHbqaJRQ/DfRfCoRfCwWfrfXkp9FWz9aByrcci1lkq9C1/har9ocXnf+qF\nVHz07wA/SUWZ/0UdZVT2E/JQn1bmRSDzQ0NKdQL/cbyvA/YCPgjsS7QsUKCV+d3AODRpj6FbFbvZ\nowtN6O6Hhamd4kQzOkXxEWhv3KsOv3o9iN3Ugiksg9otUDW5534vuyIJqbvr27gQ6neEIYN8OhrQ\nZlj7QW1vXAlrXocpMwxOTIC2jfDjPWHjcqgZCt9+G2qG6GNZE3pbI/x9b50e9+NvQpvHd5wmodvv\nu5bCpt2h6iBoeMJcoadC5nOAWamQ+R/UJ4zKniJ/6dNkXgSbxT3A6MaT6FWU9gb2QStoW11utjYb\n2xm0J/inSnD/xq/E26e37ZQgUnfnh60HPuqqw32+X71t9CL0IAvGSYAV4/VXZJIC18sCiWu9DNnd\nv52wNsPad8PZxpDt9ZY1HviaJnLQIZcv3gEHf1H/GsjSclEFePo0PdHrsEegMuRhmQaq34D2KdB2\nEbAZaq6Jdn6ksEU/hIm+MtwooWgWGwfQvQBHGN7z2GdC8G6YJup2euQmxO6uy4+8vfbbisWA1IM8\nda/kXu6oF/AmdXf3oyTwCiN0rzbD2o/SRpp47VF49rcwYRp84Cuw94lQVd2zTFaE/uxP4N1HYerX\nYfsj/NtK0z9fdzrUHAcd98LAs6DqwPDPERvlCJe0UIJkHtVGccNN8FHJ3XQRzChqPQqp+4U1epC6\n6UBpEpXuRejQW6XHmWxkotKLTeid7bBkjvbkJ0yLX08UQm9bA01v60He+ZfCdtNgn6t6npOE0INQ\naIItz8CWp4EBMPCzUFgPFcOi1ZeKOk8H5svG9W2UIJmnDSe5RyH2NNV6FFIPsl7Ac2ZpmPWSRKX7\ndd2t0oMiXvoyoVdVw1FXhZeD9Kb+v3kPvHk3NL0DVbUw60+9w0qTIEidt/8bPXYE0AGtX4D6R4Bh\n3WVTJ/SyOk8DfXBxiiR4z7FFgekimGFZGf3q8TrHry6Pi94ra6NftkY3/BbKcMOr216DblFmGAa1\nF9Z2lPrzRBrZFhf/UcfsN70Fe3wJakaZZz/0GqyOkmGxfXb3+wH7wpinoHInszrLKCq2MTJ3Ig6x\nmy4zZELqpud47fdJz+tH6k4Epd51wivVbZQUu26kQehBX3t/IfSOt2CVI5rrlZ/Dknv96zUldFNs\nttJtVM+EUbOhcnRg8VAY//bPLha8lRqjzQ8iMktEFonI6yLyzYByB4pIh4icEHauiHxXRJaLyAvW\nNivp59yGydyJJMQehCBSD1LpfqTuRgCpu5tyI2uVHpfQk6j0/kDoi/7U/XrobnDCc7Drad374lpK\nJuq8ayN0vAA1J8DIB6EiJDfQNqDORaQCuB44Ej3Z5WQR2c2n3I+ARyKce61Saj9rezhpX0PJXERu\nFpFGEXnZ5/iHRGSD4wlzWXCNqQ+Hp4yoxG6i1tMi9SDrxUXqprZL2irdXZcbXvW70V8IPQxuMlQK\nFv5Rv979M/DJubDdnuH1pGW3tP4bhn4BJt4N4gqBTELoJaDOE2AasFgptVQp1YGeDXicR7mvAH9G\nTzc3PTfVmHaTr/lW4BfA7QFl5iiljjVvNiqhF0sCRB08DYuECYqACRokjRmf7hfG6DW/Ka2IF3e0\nSzEGRkslP7pJvnDnYOKal2DDm/CRG2Gvs3QMu9+aoiZT8aOmzK3eHeqP8Z9JmzQhV5GiW5KsAQrs\nACxzvF+OJumtEJHtgeOVUjOtBZ1Nzz1PRE4FngcuVkptTNLRUGWulHoKWB9SLONZU00RtywQRbEn\nVeqm5Q399FJV6es26gWk/dAfFLrJQ8V+KK2aC599BvY+u5tQ/R5oWdgtAyd3txu1/tT0Vkmq8zD8\nDJ1SNQpuAHZSSu2DXhj52qSdSCs08WARmQ+sAL6ulFqQUr0x4XeXp3XF2YQeptbjKvW04tNLQKV7\nxaN3tcPrv4K1z8EH7ygrdNCf4X2f75n62EbQxKKwdvJKl2uCIqlzvzjzt2cvZenspWGnrwCciWnG\nWfucOAC4U0QEnYzpKBHpDDpXKbXGsf83wANhHQlDGmQ+D5iglGoVkaOAvwG7+Bd/0PF6Z2vLC14k\nn4TgTW2YPEjdcBape7JR2rNHgwhdFWDBHbDoMmhZCkfP6647yYzRIEJf3waNc2C3IwMaKBEMrog2\nUzQruyUMSeyWQMy2Nuidijp9TJoxkUkzJm59P+eKp7yKzQWmiMhEdKrVTwMnOwsopbbGborIrcAD\nSqn7RaTS71wRGaOUWmWddgI6XWsiJCZzpVSz4/VDInKDiAxXSq3zPuPopE2mDPfVF5fcTYjdhNSj\n+ulgRuo5qnR3XesVbHkE/ncJbHxJ75t0MmznWEIu6YxRr/6vfAHuPQWO+H8BleYEU0Ubdep/XKWc\nRJ0H7Q8j9EB1PsPaQF+zPwioKB8opbpE5DzgUbQtfbNSaqGInKMPq5vcp4Sdax2+RkT2QWf2exs4\nJ2lfTclc8PHFRWS0UqrRej0NnYnRh8j7AtIg9zBiDyL1NAdJ/VT6HOBd6DwpXZXuReiFDr326Kpb\nu4/JANjlqvh5XcIIvdAFT10DT/wfNIyFKUeWhuWSFaGbtJPGVH+TdkzqNbJb0vPOmxLaq1bY4K6u\nfTf6lD0z7Fxr/2nufUkRSuYicgf6cbmdiLwDfBed09V+Kn1SRM5F/y5qwz9FYR9FUnIP8teTkHoc\nlf46Ok/0bPSqSvjneImj0r0InQGw83Ww6RlotUTJjudA3U7m9boRROjLlsCjp8E7T+v3+36+e7m6\n/kLoJnXGVexpqf9MknKVEYRQMldKfSbk+C+BX6bWo5JHXN89bVKPotKXoKNLf4vOu3EpPZari5KJ\nMSyhltvf7mqGF4/RRD58Fmx8CsZ+p3cdSQldKXj1d/D4V6HD6pRUwH5n9jyvPxB6XMJNezDUtF0v\nlFAirv6C8gzQVBAlNDIovDHofL/wRK9z7LJtaBKfDtyIJvLhwFcxCmG0q3fCb9KPVwhjVzMsPgaa\n58D2V8Ge98KkK2DgqGxSAOzwfpjw4e73Ox8NQ8bTC0nDFpNMl48KP51gklw0yWQi0/qC9ofVm1Oa\nv6TT+fsK+jmZb3JtecGE2PMg9eXo8RgnaZ+Ht5duwSS/C4STb6EZXnMQ+fbfhtYaGH9hd5k0Zoza\nhC4Cm9fBWw/BiD1hQD1MPcv/vGITepRfB0nINu6vkKjx7dvwlP9SQR9PgRuVoE3KmywRFwVhESxZ\n2C92+RHAS2hFPhCtys8geHDUQpIQxtpmWHYMtM6BkRaRbz0u6c8YbQHUu3DfJ2BgPRx/H7xxH+x0\ndHDYYlIkDfmLYlEktUOcSGK3BCHOgGjZbkkNfYDM81TUJu3FJfs8Sb0JTeDfB/4BnAhMQpN7TcB5\nKYQwFprh7WOg3SLyEd/2GRilN6mbLHjhRehdW+D+E6FlFXziQRi6E+x/QfdsRj9Cj+Ofp22xZEHo\nJm0kiW5JMnGoCEg4nb/PoATIPG+yToooZO+1SlAepF5AE/l96PkI30NnZGigtypPUaXXNsNai8iH\nXAWDHIo86gQju04TQv/PxbDqaZj2fdjRyiTqzi+SJqEDLH4axuwKQ0ckn5CTlNBN6kszuiWovrI6\nLxqK4Jnn5WO3A4szrN8P9udah56l6/e8zMpT34jOxHkfMAu4BH03boe2Wuw+evXZiYheumqGRgeR\nD7aI3Gtg1K8OiD4w+trt8Or1sOPxsN8lwSrV7yuLSnIbV8GvPwU1Q/T7PAdEveC30EQYkixkEdRG\nHP88Q1nZSq3R1tfRDwdAO4H/ovPWrKQ4wa4K+CuwFk2KcVYfshGF1BU6euUv6EWxrwAqHWWd8Eux\n64Rh0q5NzdB8DHTOgUFXQcW3e5bJitCXvgBPnqPzfn/4tu68JlkSemcH/PokaG+BqgGGJxkg6YBo\nmtEtUZDJgGgUo74MG/2IzLvQaWJ+ivaJW4DHgZ/jmGGbY192BN7n2Bd3STkbJqSugM3AB4DLgFaP\ncm6EETqEhzB2QGeLJvIaH0Ueh9CdpO4V1VLYAg0TYdZfYaBrLCMLQh8KdHXAIafDtJN6H0uKIEJf\n/BA8d4NeYBrMSDFudIuz7tU3wIYHoDaHe2irOv8GcHn27fUziFL5EZ2IKLgupdps0ikAC4B/oZWw\nEw3A/miVmqKKMoaiZxaEsMFT07CLoKReBaCWnp/Xq113W15l3Ps8YnHtG1Bt1gsauKt1k4f7vZei\ndNfhJkpnHYUuGFyJL/wUa9BXHUR4W5OGqZ6+fJqJrHqEeBb031/vC5uWw4VvQ7XVea+HkvshZhq/\n79X/9e/ByxOgbhrs+oR3/UF1Bu33e6B2KnRywZUAKKUSpdcWEXWU+otR2YfkE4nbKyaKoMwbPLY4\nGGxtXWjCPAitSN8HTETfdZuBPdHkN9i15QH7ughrM+r3EJRTvQL9uZ3w+lWQsu1ir0zjNcko6L3J\nL+og26Wi0iwO3Y2gH0EmlktWRA49HyYPXwCv3Q+NL8PBF3YTeZy6gvZ5YfUvodAKY6Km6k6Ayv+h\nibwY4qtvowSiWSA+oTehifpgn+MF/FNpRiH0JAO1aalxN94jPO1u2FR/kzKGudKDwhfdESMmOV3c\n3QqLdIkTlRInBt3dj7SJ3EY9sGYDPPsLmH8bDBoCB32lZ5k085+4QxW7WmDNdVCzNwx2pA/OOrpF\nPWJ15m90Z1AswwRFIPMgAoJoiypHvRPjXvlRlXxQkiwbeUyNKxO6Z90mCHs4ZEXiPfpgrfHSbl1P\n1+0CZzwJI3utJ9yNpItMFLZAxUB49xboeA/GXee/jFwWUK8B/wF2T63Kcpx50RBG9k5EIX7wJtAs\nol2yUuNuhKlz8A4K9yL0sDKGk4zyJnQ3ohJ6XHWeBxodC3YNrIfP/j2YyJNiKPDqb6B2Z1j2Exg0\nCSZ8qvfDIUt1XnENyNBy3HkMlIhnHtdH3y7CFqU/WSLt+k0faF4eeZwyBvHoQVUmjXKB3lEuWaHY\nsxxtMq8cCCf/DXY4sHcZ0zBFU++8cwG8fAxsXgqjToK2xXqwNy9IsYP2+y6KQOamg5BpEb4NE2IP\najspsnxQJCH0tAZGHSg2oXuRsFKwfJHHAZ/6SwGNC3Ts/Cl/gr0Oy6fNlsWgrH/gO1fDukdgsIfN\nkvVEohQ9gyYajDY/iMgsEVkkIq+LSK/RYBH5jIi8ZG1Picj7rP27iMiLIvKC9XejiHzVOjZMRB4V\nkddE5BERGZL0cxY5ztwdYRI14iSuyo9C7EFtRTk3a8QldAgndK8yJUDopvV1dcCcL8KaF6LPRymm\nOm9cACfeBHudoN9HIcS46rzZMWt6/MUw/gL//kUdh+iDEJEK4HrgSGAP4GQRcXtdS4BDlVJ7A1eh\np36jlHpdKbWvUmo/dIx0C3Cvdc4lwL+UUruiJ8R8K2lfi+CZR/GtTQjdJNLE2aa7LSehp+3B50Hi\ncZDVwGjOHrrJgOiADfDIJ2HlE3DQD/FFltkV42DzJvjAV+Ggz+fXZlc7tL6jX4/6NEy+pvtYlMiZ\ntNcNTYiEU/WnAYuVUksBRORO4Dhg6888pdQzjvLPADt41HM48KZSarn1/jj0BBiA29DLf12SpKMl\nMgM0iZUSVd0H1WvqswchL+/dC1EeRl4zTr0UelgZtw2Ts0IPikFvehP+fDCseAxGHwyDhnvX6Ve3\nX5/ywMA6mPm13vuTrpIUpM5blgAFGDETpv+uOz1CEPq/Ot8BWOZ4vxxvsrbxBeAhj/0nAX9yvB9l\nr52slFoFjErYz2J75ia2SlLf3K9+k/OTEnsx8B7RSd2JNAZGS4DQVz8ND02HTZaAmnhMcJ2lhoqA\nWa1eMLVagtCyGAa/Dw66Fyqrzdrwwza4iIWIzEQvGPBN1/4BwLHAPQGnJx5lLoLNYpOqFyH4EXqQ\nlRJ0Fbjveq+2gywYJ5LYMcWASdiiDbfH4JcjPayM03YpouWy5nF49njocPw/JxyNEfzsllJYOzQt\n+MWdF9rh/Q/BQOuLTLqIRYnkPfdbEm7z7OfYPPu5sNNXABMc78dZ+3pARPYCbgJmKaXWuw4fBcxT\nSq1x7GsUkdFKqUYRGQOsDutIGIqQm8Wd/MmNLFYPsmGiML3KBaHUiT3qIK8bXr9owso437tuJLd8\ncFcXNZeL+/ythD4bnpqpX9eOgxPegQbDqAw/fRBE5kpBoRMqc5qG7kWSfpetF9F6nq96TxDyiqk3\nyQkT1I5v+673nZJKbpaJaqFR2aWye6/2RKQSeA04DHgXeA44WanuSkVkAvAYcKrLP7eP/wl4WCl1\nm2Pf1cA6pdTVVoTMMKVUIs+8CGTubq/Ns2xvxJlSb0rWYdkMTVGqxL6NEXp9Kzy+N2xZCwf9Ehqf\nhOk3etflVZ9XnX59sfGv78P7vwS1w3wKpIwoCayiEK1pEq4s2nHXWwJkbtUxC51+tQK4WSn1IxE5\nB1BKqZtE5DfoVWCWohMydSilplnn1lr7d1JKNTnqHA7cDYy3jn9KKZVoOloJkLkXTAkeopG8SUy1\nH+IMt5cSuUf1/d1sFifzopvQnwIagU9kS+jLvwGLfwz73QITz4CKd6F2rH9dYfUF9QPgud/B378O\nV67xOJgRomYjNFXnScg8Sjt++zMg81E6ECUUq2ViOWtiJFR5bL1Q47N5IcpgqvuYexDUZCA2Ckpp\n8LQYA6P2+/noMN2PsTV5UlaDopueh8U/gZGHw4TT9b7C2J5lTXxc08iWRQ/DPV+AEVMMKk0RWQ0k\nmi5gkcZgqBf6wUBosZA/mXv9s4wIHswJHoLJPQ1Sj3LVlQqhQ3JCj/LrZjE6Iuv9wKPoRQcc30Xa\nhF7YAovOhIpBMPmmnv5vGIFHiW6x61o2D277pM6pPmLnCBWUCLIY0I0SQdNfBpRLBMVJtGUyb8ir\nZ57Jd/wI3W3VeEWyuPfZHWvyOe6Gu3wQbBIrBeslSaQLhE8wWgjciLYErcUVmAB8kUyjXBZfAy2v\nwJSfQc2OwUm5vKJTTBJ92XhnCdx8NGyxngIji0DmUaJF/CJOTOs1jWyJUqcfUp5E1Nq6bWRNzF+Z\n1zs2J0xCyI0VPPgr+LSVulf5IJSK9ZJEoUOwQq9Hk3fBse9bgLWARRZx6G0L4N3vQd10GHde97Go\nOVzc8BzoWwt/mAXNjmiyEbsYVJYTsowHN2mrrM6LgvwHQMcHtBd2c5k+rQPTZ7oVu0lmwCgDpX7n\neKEUVHraA6MN6M91OjqiqwbYFXgQrR0yiHJRXbDsg9A6D6a+CDVT/UMWg+qxETYY2tUJqgC/2gua\n3vyrOREAACAASURBVIUtTXDB8zBuP3JHVgOhfvvzjGxZn84AaMUqs58EhTH1/XsAVERuFpFGEXk5\noMx1IrJYROaLyD6BFdbjr87dx+KodwhR8Hkoda9zvJCFSveagh+EtAdG3wJOBV5H5xw6Afgu3Zda\nBjNF118PLf+FsZdpIvc616QeG2EzTiurYPGDsPY1+OClMONyqCmSZ56Gso2izk0HQ8vqPHeEKnMR\n+QD6cr9dKbWXx/GjgPOUUseIyEHAz5VS033qUhygwr27oAdpkkgE8FDteSh1U3JNqtTd7cSJvjGF\nl0JfjfbF3wF+iCZy25/PKA698y1YtScMnAI7zoX6gf7nBi0MHUWdKwW/PRjWLISL3oHqwT0HW/Mm\np1JV51Ha8tpXVuaREKrMlVJPAe7pqU4cB9xulX0WGCIiowMrrfPYnEii3iFYvRdFqZt66klUutfd\nFHUUKYlCXwychc5J9D3gI9Z++zMFpdCNqdCVgnVng9oMY28GGRisqoOSckVR50vnwPJn4cBz9dqc\n7lmTJTCFPTKyUOdptJ8CChvrjLa+jjQGQN1ZxVYQlFXMj2SjkHsYwdtQndD1ShFIvQBspDeysl5s\nplFon7rgccwUcQh9HXABerbzD9CpKLzaTpnQW34H7f+Chouh5oDuY0G//ILm2IURuo2nrtar/0w/\nP7iuvEi9VJJXmdotZVslE+Qfmvju5bDZej10Bgyb0btME/6em32DeV0Q9s1T/Tq03Aqtt0NhI2y/\nClrrvXNq2d/AVvKwCd0ml6ghjW+hlWkB+Cnglf3OJPZqO8yI1VnPCuAG4DJgpKtMlDs7auhiBbAX\n8EHgII/j7gRdKSXn6poJdZ+HwZcHhywGJeUyXRTaWc/082HykdAw1qegYX1RcP/FsKUVPnIZDN4+\n2SLLURJjvXgpSCVMvrJnm16hiklRNRs23gUyErbckXLl/R9pkPkKdH4BG55ZxbbiwMv1Xz/VsgF/\n3gki+U1N0Hk3bLgV2p7W+6r3gqHf0Bep/RvEbtcdIm5M6tA7Y+AmoAP4A3pREgV8Ht2wX2oCU0KH\nYFJ3fpBxaGVc53E8Kuw2TUi9EjjP0VZYHHpKhF45CYb8tvuYCRF7Ieg8r4WgpxypN1OkQejP365D\nIlWXXn0oD3Q0w6LrYOQhMCXlh4fXg2PQDFj/PahoA/Vm/Pbc6Iu2VwyYkrlYmxfuB74M3CUi04EN\ndtJ1T8S9qL1IXinY+G949xZYfQ8UWqFyGIz8Mow4E2r31WqixaPdyKQO/mp9CXARsADYG02oO1uV\nDnCVdcJ00pGJSrcfDmkQuRNRSN2JqITuRARCj0rENqJOKPJqO09UN+isjEd9P536TNT5O/dCZwtM\n/px3HUknEnmhsAQ6H09QwbaLUDIXkTvQCTW2E5F30HFmA7EyhimlHhSRo0XkDfTlcUZghWmMM7Qu\ng3dug3d+By1vAgLDj4AxZ8CI46ByUHdZJ7/FIXUIUOtVaEvlF0C19fpMuu8S55U9mOCZpGkSehZs\nY0LqYYQd1LeIS9A54STfvNR5VCRV59UNcOiFUD8yvKwTSch1yW0wYDCMP15/91krXNUBhXcybqT/\nIpTMlVKfMShzXliZrYjLM12bYdl98MYt8O4/AQUNk2HqVTD+NKgd3/scP8vGi9RBX6xeYtlTrT+O\nXiHqLXSq418A1rJkPZS7s8K8CD1LhPnpURR4xDVFg5pxwkmc/UWd73gIvP9c/+NxFoLwU+dr/gY1\n+8Kqx2HKWVAVcTq816VsYrV0LWPr4H3FJCi8Ha1dP5RtlowQpE7cLoFS0PgCzL8VFt8B7ev1hbXr\n52C3M2DsB6HFx/1pxv9GteEkdbe6cw+WqvlQtQ+oDdD1NeBmYBjwO+BEerpQXp66fYXnQehZI671\nAtH8c+hB6FHtlq4N9Jp1mpY6j0PoSdT50T+A7arSH3T0wtzLoWDNP5n8Oehs1fedac6WuOh8S/+V\noVD/IGyamlLF2wbyJ/MwggVoWQMv/xFevAUaX9H7djgE9jwTdj0RBlp3kZ/C9oJN7l5qPYjUAQb+\nBjr+AV2fg8KXgXdBTgT1C8AZUu/nqdsEniahQ+mSelS7JYzQHQiyW3p0oQvWfwbG3ZONOod8FXqN\nyY3jg6hWS/t6aLbsjn99BN5/K0w6KX77QXA+IDrfAgbAiHuhc/f02igr84zgd1N1dcKCh+G5W2HB\nA9DVoUOwDvsWTD09WiIjE5K3Sd3PU7dv5M3/gvXnopX3fSBjQf4KFcd31+U7WOokdS9Ct4+5kebA\naB7wsl6SEroTAf65nzrf+G0d1VRh/XODQhXTwOZN0L4RhnjYfaUML/uj3TFHcPIZ4USedK1QG51v\nwfDfwKCZJUXA1kpDP6N7paGrXcd3BW4F9gMuVUpda+2vBuagxxirgD8rpa6wjn0XPcvOztZ2qVLq\n4ST9zJ/Mazvhjivh2MuhogJWLILZt8K/b4cNq/QaivscBx88E6Z+ROfB8PoZ57ZCbASFLwbBi9Tb\nF8DyTwJder+MgIb/6JA4+xwwGCy1idxJ6HYF/cF2gfQJPcKAqLualj9B09VQtatZ14Pizp3wUvh2\n249+HQ44O5zMS31h6EJn90LYw/eHA67teTyON2+K2hNh4D7d7QTNO88JIlKBjjc+DFgJzBWR+5RS\nixzF3gO+AhzvPFcp1S4iM5VSrdZaok+LyENKKXsV6Wtt4k8D+ZP5d4+AJa/AyPHwz1th0X/1/p32\ngeMvhUM+Aw0OUmjBzJoxgU3SzY6/7rptUu9YDW98FJRjJqdUgLoB1A9AqiJMQvIidOhfPjpEm2xk\nghgRLhvmQfOZ+nXl9j2PxVHnJuT76r9g3k1wwDkGFRYJQZeRUzG3W8pp4BA46h6oqE63H0EDoTaR\np40koZIwDVislF57TkTuRKcw2UrmSqm1wFoR+aj7ZKWUvYJ9NfpqdSbDSjUPTP5k/uIT+u8vzobB\nw+GEr8BRZ8Au+5o/8f0IfgPhUSpecPvpXZth3vHQ8TbUHQQ1R0P90TBoP2ip6H0u9HZGepC6F6FD\n/xsY9ULa/rlPhEsTUNcIzcfrXC0AlWPNlXBcdb6lCR7+vH5tOiszK3Wehmq2LZYP3waDd/SuL4vF\nK0oX7nQly9EEbwRL2c8DJgO/VErNdRw+T0ROBZ4HLlZKeeUAMUZxVhoC2ONAuOlJGBS09FsK8CP4\nICx+BCacC8Pvg3ZHXK/bU8fjvRepexI69M+B0Sz8c4MBUdUFrefTI6mLW5lD+ur8yW9Akx0b3WcT\n7nUr5vb1sPfFsONxen+WtooTebXjxP9mw4LZmTahlCoA+4rIYOBvIjJVKbUAnXfjSqWUEpGrgGvR\n08ZjI38y33MvGDYchg6DRf+GDxyRTr1xvXIv62Vn60LegB66cNcfldQ9CR3678Coid2S4pR/0PlD\n1J0w4FSd16Nyb63MITt1vvQxeOnX3ftbI5B5MbxzEz1QPw4O+mG8+tMaCE0bfm1NmqE3G3+5wqvU\nCvSyWTaC05X4QCm1SUSeAGYBC5RSaxyHfwM8ELVON3In8/pnH9v6urWpBuc3XTBh4yg3gRfB20rd\n1Hqxz3HyaFRS70XokP3AaKnBRAbHjHDZqs7Xwpa7YcAxUHsLDFzcXdxvIlFcdb6+CR5xCSmR4k33\nTwN1AB6/ZrxUczGUdHEwF5giIhPRaUE/DZwcUH7rE11ERgAdSqmNIlKDzgv9I+vYGKXUKqvoCcCr\nSTuaO5nX1rb6HnMfKWysMyNuL4L3mvgTR7l7kbrNo6ak3ovQIfuB0VJU51HI2gshE4oqbgG2QPWX\noGIEdIzQw05xESQcujbDSU/CvUfrkMRCJ1vv475M6HnBNPlWkaGU6hKR84BH6Q5NXCgi52ClNLHW\nb3ge/V8viMj5wFRgLHCb5ZtXAHcppR60qr7GWpWtALwNJB49z30N0Ilq4db3bfScJuxeRVsr957w\nTCJvsnKJ/d7mP/tC2uA47j7mtXiB82ekm0udF6fXuc7yW0ndnmjkJG37dZPHMTfCFHoehG73wYvB\nvOwWdzmvdUSDjnutUlQApgAClYthsGOg2m9FIb+cZH7rhbqJfeObcMcU2P8i2OeLepZkww7eH8EP\nUawWk9mWYWQYdrkE2R+mKwQlWSfUWd+ydFYa4ieGHHdx8vaKifyVuW9KWMAgBYSReje1YuyomLBB\nUq9QRrf1AsGTj9x5X3IbGM1aoTe5XrtZLC///BHgLai4RoeQZgH3dfW2ZXNOORaGudYALZY6D1O3\nYZdL3n52GamhCGTub7N4FO6B1tZaaht6Pgy8avP13r2INQ6ikrqzXfDx0bMcGM2K0L3ai+ONm5QJ\nI/RfAdUgZ/Q+3S+rop93bjoQ+vb9UD0Mhh5i9hG8UOqTiNJC3ORbaaDErJuskDuZNyQZuDNR7k01\nVAzpeYUYe+9esL12NzGDGan7nZf7wGjahB7l/xjXPzf12JcCfwdOga4R/jND04D9f2xfD+/OgZ1P\nhoqA26gv+udRViLKOvlWGcbIncxroijzELRRGzigasMuYRQtE8d6gZ4XsDvyJch2yXVgNC1CDyNy\nU7slzoCo15T/m9ET676kd3vNDE1bnb/zsI5tn3Rs73qiIm91XurBUCU4ENoXUFqeeUJ4kbvTmvEk\n9STWS5BKh56JvJxrl/oROuQ4YzQuTOuNK0mjLGgBsAa4Hb0GqXv9UczslqhoRlssFQNgvGP5OD9C\nT1Od56V6o3jnSWaE5uHRl/KDK0WUtmeeEJ7kbtkwsawXJ/nWufYHWS9eselFGxhNos797gp3/nY/\nmKhzuz7TAdEH0Ynn/g+9UrhH7vMwMo2qzru2wNKHYIeZMND1meMSel/xzsuquWTRp22WKHATu5H1\n4lTVJsm+/KwXP9sll4HRtAg9jMj9zoljt5jUZX8Xt1n77TStAasTgZndEobV/4aOjd0Wixt5KPQg\nmBCuyY+3LFRzlIHQMiKhqDZLLW20Bt14EdBEQw/V30qtt6XjjpCx/hao60muSX6Cu1W6n+2S28Co\nF6IQuimRb0DHe9cC9jqsaYUretU1F3gGndLCgz2jqHO/ZrzU+fL79fuRHzOssJ+hrw2EbiOhln0r\nmiUANbT2moRU2NxO81MvM/jwA3ufUNvtp2+1XuKMYHlZL13tsHY11I8Ptl1SHxitADqASkcHkw6I\nmhB5F/B79Hqo7wc+4VFHnHDFsAHRe6y/n3OVDVhqDpKpc6Vg+QMwbB+om+D/4I+jzoNEROcW+MdP\n4PDzMQrrSgtKmWeDTLvdMiKh33jmrdRstXCUUjTeOZsl37qVLSveY9qim6idvOtWtb61bBSV7iZt\nP+ulsx3+eiKsfRGOfgHqR/ZWLF4DozZiEzrAlegBwR/TnSHM2WBUmCrySnQuot3Rg5DN5OOfX4DO\nRrqLT9mA6sN+ffmpcxH4yBPQtqq7bJqE7tufRpj9azjgRKieEvHkBPjPObDubfjow/EmY8UZCFUK\n3to7elvbOPqsZ95GreeDYe3Tr/PyRX9k/XNvUj12GHve9EWGTWpgs3X1uAnd9tKNVLofodthjBva\n4aETYekDsP+FsN2InnaKieXihCehu2ET+n+Ax9A5ewZ6lMtiQNSJGehfBn6Ia7cE1VUP+GXdDFHn\nTkSx1Zr/f3tnHiZFde/9z2+GZRYWAWUTAwoC7mjU65Ko4ZWIW9T3OjHiva5JJJGYZ0QjGrfovXGd\n1yQ3xCVXjd6rVyUhyrziHnFNEFmCiKCiAdlE9ll6tp5z/zhV0zXVtXZX99SM9X2efrq66tSp06e7\nv/Xr7/ktQL/R+lEIuI1l0D5w91q9HZV04XePVwo+fQ76jckm8kIm32pZCc3vR9CRga/Igm3RyXxP\nttGYw99ErYm7uzXWrdnCkllzWffHJZRW9OGAm89h/DWn06uyDGhBHIjfl9BNjxenH5id0Lc2w5+q\nYHUtHFsNx9VoK07obN07WSpei6ImPL1cdqGt8qHAlehfaRj9PCyc9PmoQuhzkVtCJOwK49lSSOS6\nGFosDXrFozB0EtRvhEO+X5xrtjdC0zJI/bU41+th6LaWuYmmHSne/ffXWPabd2hva2f8JUdz1G2n\nUbL3SPRiXGPWYqh9sbQDTjq6ndCdiL21GZ6o0oWoT6iGyTWZdKgmMTgRupdlZJdbHGGS6r1oeeU+\nMgzhRLhhrfNiRJdE5d1iIgSxh5Facu2rWDeHqPHJn+H1q/V2eyusfgwmXFjgi7bDupOh10goGQDt\nERkkiWVeGEQVNJRuTbP4/iW8cctbpLan2GfyWL5ZcxpDJ400PGS8bxomoXvp6L6EXg+UNcPdVfD3\nWphSDefVwC6xXkjDfp6d0EPJLaZ1DvAm8ATwHXTNWetFw3q4FBr5+OYF8W6xw8dN0Qo/YnfzOd+2\nGMqGQeWoYNfJB4Wwyt3u16VlkNqqtxfeDme+6tAo4utW9gKVgtY10V7rK4JYebM4eaTYoZRixby1\nPPuzRWz9aAd7TRzEdx49g/1PH4eIkIsl6aajZ45rZC2M9mmGX1TB4lo4rRr+tUZXm3FKB+BE1H6E\n7qlBlqOt2p8ZHf275ZiXRV1o7dxEmFzluVrnIa7np5tbEcaaXnE7TPyJJvPuEvgTBL0sN8EDLoBR\nk4Odl4+LonRdFcuegC7xZjEtYtDkXhfwR7t2yTbmzFzMRwu+oHLPMv75t8dz3A8PoKV3P8jT4g+9\nMNrSDDdWwbu1cHY1fL8G6iV7YVS/SQ2nIhZu8Ar771gMvQ3tDvh7wKxVGiSPSxj5JE6JPPxuECYL\nh/Bq8YNX210fwrq5sPdpMOxE/2HlM45io5cRJ9B3IJxU49wm8mhQi0ttxWRo/Es03cbB170I6Baa\n+bb1KR77+SLe/q91lPYu4fRrxjPl+sOp2KMP0ExjgMU3vwVUK3wJ/ctecH0VvFELF1TDDyyLnU6u\ni9aQfXuOFjfr3Au9gLZFwH8AJwPT8JZWgsotUVvnTnBjsHxdFXcAz6IXgK0IERHqB7vUsuJ2QEFd\nxLJAHBJvmZb5N34JlcP0dqGDb0SAUpDeMOIBWLO/7ynFgIhMBX5FptLQnQ5tfgOcip6li5VSy7zO\nFZFBwFPAaHSloe8qpXblM85AZO73ZkTkRPQv6VNj11yl1L859RVGM0/VtzHnrs+Ze8/nNKfaOf68\nEVxw+3iG7VtBijY8/PUyfbi4MEJ2lKjTwqid0Bu2lsDM78Jr8+HiaphhELnTD9CN0M3XQfVzu3Wu\nWtBRj+Voq1zwt8Dtx+NkcUeB/0QTugeidFFs+Aw+e0Jv132SWx9xRq9yGHYkHJZ3NTN3OH0FpRfs\neQv0KaIvvQeMkm+/RS9IbQQWicizSqlVljanAmOVUvuLyD8B9wPH+Jw7C3hFKXWXiFwLXGfsyxm+\nZB7kzRh4QynlkqwigyBBQ+m0ovYPu5h9w5ds29zGQcdUctn/G88Bxw40WhQu86J9YbQDFaCam0n9\n6CekX3oRLp8B1TYid0qeZSd0CObh4kTodU3Qvwx23QW8D/waGENmPryiQ7sDcrXOv0RHoE4Kfql8\nXRQ/vkunwIVglnl382rpXQlT7oOSUv+2QRBUN688CoZcFc01TeT39T8a+FgptRZARJ4EzgKs/HcW\nOnUnSqmFIjLQqAu6r8e5ZwGmNvcosIBCk3nANwOWqtReKG82yNyl0O7bL6e4++odrF7eysgxvbnj\nqb2ZUtUfkaZAckqUsOroZc072Vw1k/T8V+g9Yzp97riJVL1+L54Ro5AdLZqry2Lj96HvTGi/DTgG\nSq/QkfSdvFtMBJFbnC7kJLU4tcvHUyasWOzX/kF0xkRLVGboYtEhLp/aCGsfzryu7wHeF/aP+OBL\noN+I3PrKR0sf/Xto753jyQXB3sDnltfr0Zzo12Zvn3OHKaW+AFBKbRaRofkONAiZB3kzAMeKyDJg\nA3CNUmqlU2eVu9v1M/U0mEV3+8Lqle3cck0zr8xP03+AcN1dlZz3kyH0LRPCWOLWsP4wx9z7q6Cs\neRfrqq6jpfY1yqsvo1/NLFIpCR4xCp21cnAndCus1vmubdDyP9DyjD5W+hCI1WoKkjIXiiu35EKo\nQaNCzb43o10zAb5AF6mw2xUR6+ZqIxz/Kvz1VOh/INAMzduh7+Dw/cV1EdSJyItR5q18YvESY+1e\nAHULCtFzLsls8k5GE9UC6GLga0qpRkM/eoZMwozOsHzwlbSz5Uu4+a56fv/fet+lV/Tmmpv7sOde\nQiNNEQ3PHX4Lo+3NLXxadSOp2rfpX30Rg2pm0iQSPGLUrpVbXRbtPujgsSD6IjoIqgEoB7kS1EPQ\na7TD0kG+HizFWAgNCzfWewidXAw0aW8APHy+o9DNBx0JTZugrR72mgxH3gLt/us33U5q6SpEnXrX\nra/Sk2CPkzKvN/3CqdUGdOIhE6OMffY2+zi06eNx7mYRGaaU+kJEhqOT8ueFIGTu+2aUUvWW7edF\n5HciMlgptd3e2S136+e2NtiyG558Durq4cxT4K6bYOL+rTQMSEOz4Xfet4gZ4gyY8kq6uZV1Vfey\nq3YhQ6vPY3DNlYg0IUi4FAB2rdzNB91s60ToG+Zn2ksFlN8CKSM/iGMwkYkgZeeCEH3UFryXSRrU\nOl8FrDTabkcbN5vJkHmeLope5FtnqIz9JkJpX/2wojsuhBZ6ecUv6ZbVUm4lLlgEjBOR0cAm4HvA\n+bY284ArgKdE5Bhgp0HSWz3OnQdcDNyJTvv5bL4DDULmvm/GvMMY20cD4kTkADdfCE++DLN+B+s2\nw6SJUHM1TDZjEuq1xd6BARZZxEVnd0MjFTmn3E03t/JO1X+yrXYpo6rPYWjNDERSjgujHdfzy+li\nlVb8fNCthL47DQ0v6GO9JsBez0HTWJcfn5fcYoUfoXcH63wY8DhwNjACuBT3//x5Si32S5tk3n9i\nYYi7O94MrMhFfhlwkn6A/j1sdbSUwyMPGUgplRaRGcBLZLz5PhSRy/Vh9aBSar6InCYin6BHfonX\nuUbXdwJPi8il6Irk3819lBq+ZB7kzQDnisiP0PfTFJmyL1k45iJ49yMYMQQevgEuPBVKS3H9KxQV\nsZsIopunm9t4ueoRNtV+wLjqqYyruYAmyXwjQkeMgrNWDsF80BsWQXobVHwLRv0JUoN8golMhHVX\njBOCWOfmF+dz9DLOGcb+PKNBg4T21xm/yX4TfMYYAHHVzYOih5eSU0q9AEyw7XvA9npG0HON/dvR\nQSKRIZBm7vdmlFKzgdlB+lrxD7jlfLj6HKgsgzCyeNTEbsJK8G3NbcyrqmVd7YccUn0ih9WcQ8qQ\nVjqfEzJi1EoQJpcG9UH/dD4MuQRG3w+pPu4/Hke5JYy7YiGscy9izSmxt+2cZcaz8xJNpNGgVtSv\ngr7DoY/xITrdANxiD7qzbl6MAszmdRKEQtEjQD++FUbuDTQbjxxRCGJva07zTNVrfFr7KYdXH88J\nNaeQklQWkdvhReiAc1UbCO6DPvA42PeUzukCINs674Rc3RW9UGyfdTfr3MrAZsDOBPyZOUTiLRNu\n5LtrlZZYehq6Iiyh0Nfswf8arCg6mY/sjfOdfRcw0GF/AERB7K3NaR6reo8Pa9dzfPXhfKvmBETq\nOxF5LhGjjZDR0CHzxQrjgz52arb/eaDcLSbCuCuG+WVFIdPkqzGsQXuC2SMGc3CLDKqbt9VD8+cw\n8vRw5yeIT13QHojipykr8F+0IMRuD/FvbU7zq6plLK/dwsnVB3B6zZE0STDfdq+IUVdCh3A+6NB5\nQdQ8Jye5xYogcktXL4T6WeefoL1XQlrc+ejmqY/0c78AlnlQqaW76+ZO6OFaetzQNWTupxtGpCtm\nEbuDpd7S3M69Vet4r3YbZ1SP4byaAxw18iBwyo2eRegQ3gfdbGu1aiKTW7yyKzohTukB0uiskcdb\n9gUJ/7chrItio+HJIgWWWbrK0s/lI44zccfl61pgdL1lHiWxeyzOVNKeRegtzYpbqj5nYW0951SP\n4gc1oztp5LlGjHoSujX0H3LzQQf3CFJHuSUXd8ViWud+bOpmna8CWsiWWLwQQTRog0HmlQ5kHnQh\nNFc07YayiNIU5IpCLILGyUbopig+meeV5DEErF+4Sjr811X/BiiD5mbFzKqdvFHbwgXVg5leMzRL\nI3dCLql0TUIHMh4u+figm8ecIkitiI27YoR5UjrhM+N5nMs1IixYYaIObZmXlEPffTqnxQ2DXP59\n7lwPix6BKTfmcMEYIdHNC4Kut8ztx/Kx0gP8QO6tUVxxfQOX/UsJr9Smuay6nJk1/QJr5E7wWhi1\nl6Pr5LKYqw86ZOvn5jG73JKXu2LctXMz4/JY2/6IpBY3i7pxFVRMyK5Y73eeH7z+oCgFf7wchh6Q\nQ8c5ICpLOQ7yS4BMCz0BXUfmhSbuBmgrg+f+Bmd9W+96+R24+S54/R3Fq2+mmXEl3FRT0qGRe+U+\nD4N0YxPbNzRRsb9O19Bo/K03A4si8UGHzvo5wK4voaIcGvtF6K7oR/KggwX6QF5ZLXNZAdyCdoEa\nHPK8PFwUVTu0bIRBUzLHgljn+UotX3wIm1fAyMPy6CQHtDbqqkNuN65CId3Vd4DuhyJ/QnS+SzfY\nnqGzDON03Gnb4fy2NEz7BZx9E7y4ANrb4dpfaQPn1Tfhe+fA7XeKUTdUF82w6+NOcopTG/u+zfc+\nzcoD/4WG91Y5nlfR39av/UfudZOyt7Xy3+7bYOMIaP+ycxvHW7aTlZoLmW4E5gJ3E0HiNx/Y/xlc\nC9xHJkmdkzyUg2TkxSNSAodugQm/z70PE2F05+EHwnWfwOS8Ul4HR8NWWPk4vFYNv66EtmZo8THV\no1ys/fDICDtLBXx0bxSfzMGd0M1tN0L3InfL8bbtMO0emPM2TJ8KU74OT70CSz/MnPb0s/DgAzq/\nut0ar6AxK6eLrl3a+QN3yvtSsnkTn97+Z/pN2o+KI8Yb52bOs4b9A97k7HTcDe3NkHoceh8G5YGP\n5gAAEKFJREFUJXv5t48MD6J9vUuA5+mcU7wYyJNBcvkLLiXQq0huJtbfSmnv4i1+lg+Gl34IH/xB\nZ4R8/Gj4/PXCXtP87qs0NH/m2TRBNopP5ib51pP5ojqRsxOhW7ddrPW23TDtdxkinz1d7/v5fZk2\nZ38bVrwAV03rbEm6Wef2feVk3wDMfR/c9CfSDU2MrfkBA0oaOrWz91My0Gaa2fnB/tfdbrFbj9fP\ng/R2qLwks8/R0HaSF5wIwmmfXbMeiGbDZvQHth6d/CpBlyIKhaKkBPY8CNItmszbmmC/0yLoOABa\nNhilEROEQddY5lYO8yP0XR7HbdZ6WxqmPQRz3oXpk2H2+VBSBw+8CJ9tgm8cBm8/CH+ugQMs3mxO\n1jk4W95O0ou5b8f76/nHQwsYec6RjDhhrOd5oaUWP4u97mGdGrfiu0X0Te5LZ+f9KTjn5Y/Smsxn\nETYiyWX3K1qv6+kYflBm+4ifFkc3//IhSC2PuNPdAR/dG11nmVuJuN5hn5OO7iG1tKVh2oMwZwlM\nPwFmn6uNi7oU/PlNmHcDvFEDxx1qOa8+U/nICi/L2wq7HLP46j8iJcLBd34vq12n11FLLb3Xw+6X\nYEAVlLjo3oF187AwczDsjy55+BXAlvugfmnnfT3R1W6oQeZlg+CgizofK1QirPo3YM25Beq8Z6Pr\nLXM7oVv3BVwYbUvDtEdhzt9h+rEw+0xN5DRAuh1engVnHq1rL3ec57D45KaLZ1nVDtLLlheWseml\nlUyY8S367z8c0DcAu3XvK7UEgZ3cNz8GtMPAS7Lbhg4RdzrBqxPzhhBpNs8IkINfXZBT2pv1jXPr\n/w/ff3dz0DDJ/NAfQp+I2dvN+6dsIqg8MvB9hREPy9yP0D109LY6mPYkzPkAph8Fs/9vhshBf2dK\nU8a5LiReubs9U2jagJMlDtlWNkBl2y7emDmfvoMqOPTGM4x27nJMx2u71ALhdfOBCjY/An3HQsUJ\nmfM9rfmQbnmeGAgchC4Vmy+6QQhg3QJor4dttV09ksJjr4OgpBcc7piquzAoK0SKhERmKRy8FjRD\nLIy21cG0eTBnNUyfBLNPh5KUwzluY2ggy1pyss71fveF0KUPLWP7yi0cfuO32WOwBLoJZEktkJtu\nvu1tSH0CQy6Gfj75ZDyjCsJKLuZgBgH/J+S5+cLUzdPGI5ebQA6uaA3AToPE696D5o05XNcHcbqf\nDdwHDrkM+nvUVI0aBSHzwkBEBonISyKyWkReFBHHvK8i8pCIfCEiy237bxWRv4vIUhF5wagFioiM\nFpFGEVliPH4XZDxFJ/OUl/4dYmG0bQdMe8lC5N+2ELlb3zla515uik27m3n9pjfYY9wQjrri61n9\nBpZactXN1z4MCOx5kUsDP3iReBCC/yYwNM8+csUqYHUB+7dBKdhlkVe2PefdvrvJKnaIwNRfFvea\nfccCpVB+SHGvmxtmAa8opSYAfwGuc2n3CHCKw/67lFKHKaUOB54DbrYc+0QpdYTx+HGQwXSJZZ7y\nCgwKsDDa1g7T3oI5n8D0g2H28VDSaOvHqW87IrDOX7/jPRq2NHLynSdS2qeXpa17H51e5yK1mG1a\n62HD0zoaccg+nY9ZnwuaWnVQITv3wbvGww8RebE0r4CWtZnXW2MotUR9A6kIG13rgDDeVSV9oO9+\nMPya/K/bgbqAj9A4C3jU2H4UXYw2C0qpt4AdDvutn1YlWNO8hk/bWnQy320Qa2qXzUq3Plu3bTp6\nWztMWwBz1sL08TB7EpSI7RyviNGQ1jk4BxEBpNZt5e17lzDmm3sz8ZwJRttswgdzITUCqcW6b+0c\nSDfACIeFT1946eZe7O9Xl7NY2IauNb6oeJdsfB2GXaW3+0/WvtDp7h85GDsMOBkGuZYRjhOGmoXs\nlVKb8f6L6ggR+TcRWQdMA26yHBpjSCyvicg3gvTVJZb5bguROhK6y8JoWztM+yvM2QDT94XZBxlE\nHiRiNKR1Du7EbMovz1y/lNamNKfVfBMRcfVLd+sjC35Si33fmkeg9x4wztEgcEakunlX4gt0kNLn\nwCbL/hys8KBRoIN+DCNv1dvlh8Chz0OJxc8+F/fEHJyZejz2vk1b6AXHe8ADlkc2RORlEVluebxv\nPH/HoXno4AOl1A1Kqa8BjwM/MXZvAr6mlDoCmAk8ISK+/2+KnmirI1ef8SUeUKkJvXwgnZNnmdvG\nc1sdTFsNc7bB9H0sFrnZzlp2zqkfLK+tsFm+lbszec8bzVSHBkyybqScNe9uZ+Hjn/H1C8Yx/qgB\npIxqQyZMsrb3AXTkO+94bS9eAc7Z5qz7dn4MW96E8T+G0rLO78eJIPLOghe3hNPLLNvv4l7QOULY\ng2ZEyOHfcAI/9Ir635/bDX48nb83D2a1UEpNydppwFjUHKaU+sJYvNySxyCfAOYDtyilWtCJ+lFK\nLRGRNcZAl3h10DWWOZnptcougKNl3lZvIfLhMPsAi0ZuaecZMZqndW6VSJRSPD5zOb3LSjj3l4c6\nnmOHuRCak9QC2br56j/o53GXOjQmoIuiiUJZ5IW09P9u2X4vxHkRuqDF6d7WU9C9SufNAy42ti8C\nnvVom3XnFxFrVZWzgQ+N/XuKaMtBRPZDJ+z/FB8Uncyt339XQofORL7WIPIhMHs/KPGSZMA5YtSK\ngNq5W4j/0rmf89Fb2zj1qv0Z8rWMNR00QZeJQFILZOvm7Wn46FEYfAgMPsI6uGDo+D8WxN88CCFH\nSdpBGLIR+MjyegUFdR2Jsuvu6OEShmBjWdC6YAugdwJTRGQ12j/3DgARGSEiHW5PIvIE8A4wXkTW\niYi5yHWHIdksQ0fd/dTYfwKwXESWAE8DlyulfEW84uczJzNt/eksu5iSC2jZpa0eptXBnHqDyEcZ\nFrlROajji2OTZIBs2cUKp0VFaxuHL6RJvLtaynj42jUMHNqHM2dNQJGxxFMWScVpX+ZYpnBFxz4n\nqQWc1Y0tL0P9Bjju6kxYa5DqLaGUkrjJKlZsRcuLc4FS4F/RGnosmaR46I43im4MpdR2HEKflVKb\ngDMsr6e5nO+Yt0ApNRf95Q6F4nuzWLbrbPusC6N1O2HaDoPIB8LsQQ7uh2EiRu3wsc4hk4DLalnP\nnf0lm9Y0ceFtYxjcvzVwgi5w9jkHD6nFiZv6AUsf0ZF54y/I7LOiks77Q3Fcd1gE/RpgFmkQ4Eiy\nqw1BwUL6i42EpBMEQJdY5taSwmZ9GauF3gZcWQpz0zC9AmZX2NwPrVa4tfSa08KnaaEHsc7N89D1\nQhsGdL7XNW2r4w+3bmTfg8o45dLh2DNIuC16Olnier+L9e60+GmOd+N2WPUMjD0TyvfS+/L6sTsV\new6KOFvvMUGUxZx7GopWC7Qryx0WD7HQzK0W+nbgR2gi/34vuAeL10oeEaOOsFvnDha61Tq/77Zd\n1O1s54p7RlHaS0+duTgaNMgInH3OQUstjom3rJrlkid0jumDL3W+Ifnpm+bxvG7jQbwNim3hFymk\nP/+qggkSFARd4s1iXW6wEnob2qmyFrgQuLVNE3m+EaOdzvXzbDHb2dLjfvpRG/8zu45jT6lk8tTe\njtkUTbgFGbmRO7hILZBt1b37MPQbDodMzW5r9XjJS2pxQ/dyNegS9MRUuAm6BQKRuYhMFZFVIvKR\niFzr0uY3IvKxiCwTkUlufVn/8FgJvQ24El147HzgemNw+UaMdjznaZ3ffe0u2tth1j3OFqdpnXsF\nGTmf50LwTpr5umWwYSkcdSGU9spuHzlM8u4OOrofIq4FaqIQKpNbn2sXFuBiXwUUzJslVvAlc8Pf\n8bfoRDEHAeeLyERbm1OBsUqp/YHLgfu9+txGhtTr0Hnvfga8AvwzcJUxMKeF0TARo532WV+HsM4X\nPd/Om28o5j+T5rzLyhh/cCYyzbTO3axtL+s9tNQCsPAR/XyULXzfSuQ7Fri8IS+YWv4A23Ox8HGR\nrxcBWtYX93pKwdwrtFtqLvhsQaTDAQpXoCJBTghimR8NfKyUWquUagWeRCeYseIs4DEApdRCYKCI\nOBaDtN7/tqGJ/GbgdeA7wI+NQTl5ujgGGGHbthK6XXbJwTp/7S24fpaish/8/NYSxyAiK7yO+1rv\nXlJLWTO8/TiMPRaGWe6l9h/UzgUZg9pNaolEN48SdjLP10oqYEi/iSjJPEhI/6blsH4xrHNJLOb3\nL+IfC0IOKkF3QxAy3xudAMPEerIrEdjbbHBo0wHzT00a+CXwNjAV7TFvtcithB4mYtRxYdSKENb5\n8vdh6RKYebUwbHj2dFmt89BSCt5BRZ2kloW1ULcNjr8kcyyRsL86+MDI0Lgyhpkag6BLPXqS4hQF\ngTllabR2swg4CZhOZ161E3qnPjwiRu3PKxpg9Rb8rXPIss4bU/DqOzBqJMy4MtPMj7zTre1sWPwF\nrRu2uPqcg7s/OjiUk3v5EehTDkef515yq2UHqI3Q5vUmC4Ed6Pu3uWBcLJlmJDkkqiM//bwEKo6A\n0jzeYy6upBWDYcxx0H947tcNA6WgNQV1m2DjYi3vtBYpQ2R/oKmIeep7CET5VBkXkWPQyV+mGq9n\nAUopdaelzf3Aa0qpp4zXq4ATzfSQlnZfgZLmCRIkiApKqbwymYnIP4DRAZuvVUqNyed6XYkgquki\nYJyIjEanZvwe2uHEinnAFcBTBvnvtBM55P/BJEiQIEEYdGdyDgtfMldKpUVkBvASWpZ5SCn1oYhc\nrg+rB5VS80XkNBH5BC1S5FItIUGCBAkS5AhfmSVBggQJEsQfBVkAjTLIqBjwG6+InCgiOy3Vsm/o\ninHaxuRY8dvWJk5z7DnemM7xKBH5i4h8YFSYudKlXSzmOch44zbPItJXRBYaFerfF5GbXdrFYo5j\nDaVUpA/0DeIT9KJDb3RJmIm2NqcCzxnb/wT8LepxRDzeE4F5XTVGl3F/A5gELHc5Hps5DjjeOM7x\ncGCSsd0PWB3z73KQ8cZxniuM51Lgb8DRcZ3jOD8KYZlHGmRUBAQZL8SsPphyqfhtQZzmOMh4IX5z\nvFkptczYrkdXgrHHT8RmngOOF+I3z6aPb1/0Op5d+43NHMcZhSDzyIOMCowg4wU41viL95yIHFic\noeWFOM1xUMR2jkVkDPqfhT1BSizn2WO8ELN5FpESEVkKbAZeVkotsjWJ5RzHDbEJ6I45FqOrZTca\neWieoSgVhL9SiO0cG5XR/wj81LB4Yw2f8cZunpVS7cDhIjIAeEZEDlRKrezKMXVHFMIy34AuBWNi\nlLHP3mYfnzbFgu94lVL15l9BpdTzQG8RGVy8IeaEOM2xL+I6xyLSC02M/6WUcirYG6t59htvXOcZ\nQCm1G3gNnd3DiljNcVxRCDLvCDISkT7oIKN5tjbz0CnLzQhTxyCjIsF3vFZ9TkSORrt0bi/uMB2R\nVfHbgjjNsQnX8cZ4jh8GViqlfu1yPG7z7DneuM2z6Er0A43tcmAKsMrWLG5zHEtELrOobhZkFGS8\nwLki8iOgFV2e5ryuGq8J0RW/TwKGiMg6dPLJPsRwjsF/vMRzjo8HLgDeNzRdhU61P5oYznOQ8RK/\neR4BPCo61XYJ8JQxp7HkizgjCRpKkCBBgh6ALikblyBBggQJokVC5gkSJEjQA5CQeYIECRL0ACRk\nniBBggQ9AAmZJ0iQIEEPQELmCRIkSNADkJB5ggQJEvQAJGSeIEGCBD0A/wsqLaJgqn+HNQAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fafe8f62c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xm = fine_grid[:,0].reshape((10,10))\n",
    "ym = fine_grid[:,1].reshape((10,10))\n",
    "V_mesh = V.reshape((10,10))\n",
    "contourf(xm,ym,V_mesh,100)\n",
    "colorbar()\n",
    "# add gradient field\n",
    "# see api of quiver at http://matplotlib.org/api/pyplot_api.html\n",
    "dV_m_1 = dV[:,0].reshape((10,10))\n",
    "dV_m_2 = dV[:,1].reshape((10,10))\n",
    "streamplot( linspace(a[0], a[0], 10), linspace(b[1], b[1], 10), dV_m_1, dV_m_2)\n",
    "quiver( linspace(0,3, 10), linspace(0,3, 10), dV_m_1, dV_m_2)\n"
   ]
  }
 ],
 "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.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}