LogGabor_fit_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LogGabor user guide\n",
    "\n",
    "\n",
    "Table of content\n",
    "----------------\n",
    "\n",
    "* [What is the LogGabor package?](LogGabor.ipynb#What-is-the-LogGabor-package?) \n",
    "\n",
    "* [Installing](LogGabor.ipynb#Install) \n",
    "\n",
    "* [Importing the library](LogGabor.ipynb#Importing-the-library)\n",
    "\n",
    "* [Properties of log-Gabor filters](LogGabor.ipynb#Properties-of-log-Gabor-filters)\n",
    "\n",
    "* [Testing filter generation](LogGabor_testing-filter-generation.ipynb) \n",
    "\n",
    "* [Testing on a sample image](LogGabor.ipynb#Testing-on-a-sample-image) \n",
    "\n",
    "* [Building a pyramid](LogGabor_making-a-pyramid.ipynb) \n",
    "\n",
    "* [An example of fitting images with log-Gabor filters](LogGabor_fit-example.ipynb) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing the library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:43.796812Z",
     "start_time": "2018-06-20T15:16:42.763544Z"
    }
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "from LogGabor import LogGabor\n",
    "parameterfile = 'https://raw.githubusercontent.com/bicv/LogGabor/master/default_param.py'\n",
    "lg = LogGabor(parameterfile)\n",
    "lg.set_size((32, 32))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To install the dependencies related to running this notebook, see [Installing notebook dependencies](#Installing-notebook-dependencies)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to top](#LogGabor-user-guide)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:43.819100Z",
     "start_time": "2018-06-20T15:16:43.798623Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "np.set_printoptions(formatter={'float': '{: 0.3f}'.format})\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "fig_width = 12\n",
    "figsize=(fig_width, .618*fig_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perspectives: Better fits of the filters\n",
    "\n",
    "Basically, it is possible to infer the best possible log-Gabor function, even if it's parameters do not fall on the grid\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining a reference log-gabor  (look in the corners!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:44.504699Z",
     "start_time": "2018-06-20T15:16:43.821115Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial_guess : [15, 9, 0.26179938779914935, 0.12435117250106024] , idx : (15, 9, 13, 3)\n",
      "popt : [ 14.619  8.855  0.212  0.102] , true :  14.6 8.5 0.20943951023931953 0.1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAEvCAYAAADB37lNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dabBkZ33n+d9z97p7lW7tVVIJSWhBHguoYWloW4bGbbtpY0cQjGkb4xi75RdNhB12zJhgYsJ090wEM2Gb8Yw76CkPBBBjGxPItAkMHmMbt4wXjCQktJQEKqkK1a5abtW9VXe/z7y4mVLW1X3+J5ezPOec7ydCIameOplPbifPN0/mOc57LwAAAABAsfqKngAAAAAAgDgDAAAAgCgQZwAAAAAQAeIMAAAAACJAnAEAAABABIgzAAAAAIjAQJ5XNjMz4w8dOpTnVQLI2COPPHLBe7+z6Hn0gnUTUE2snwDEyFo35Rpnhw4d0sMPP5znVQLImHPuRNFz6BXrJqCaWD8BiJG1buJrjQAAAAAQAeIMAAAAACJAnAEAAABABIgzAAAAAIgAcQYAAAAAESDOAAAAACACxBkAAAAARIA4AwAAAIAIEGcAAAAAEAHiDAAAAAAiMFDElR7+pd8p4mq7NneLK3oKADL01Mlzuvd/+HhwfGm7D44NXLPXDzcdXTXH3Zo9t5Wx8Gdoy2P2dV95rX3Za3uXwoOzg+aye/7Ovu5zb7Gve/K58O3qXwrf35J0bb993SMXw2Oj5+w7fPb2fnN8xzP28qd+xFj2Cfvz0KVp+3Yt7LLvl5EL4eVHz9rLKuFt7tre7i/70g/a42Mv2vdLv/E0vfqadXPZKnjmhXN6y8/+dtHTKMT8fvYhoH5yjbMnz5zTnf/x49ItThMnEt4oCkSMAQAA5IsYA3KOs/5lRRllxBgAAEC+iDHg1Qr5WmMMCDIAAID8EGOos/FT7X0Nu1ZxRpAB6YpxTzgAIB4EGeqs3SBrVYs4I8qA3hFiAIB2EGSoq25ibLNKxxlRBnSPGAMAdIIoQ92kEWObVTbOCDOgfYQYAKBbRBnqJIsga1W5OCPKgGTE2BaMu8Q6l9nwrH1frvfb6yQ/ZI+vG2vppR32sjOP228g/omh4NjY6RVz2SST3xs2x4fmwvfbWnhakqS+ZXt8eTI85lbtjch1+/RuOvVOe3zmYeO8dBP245X0XEo6p97yVHhsYZe97Pyt9vn49n09vPziDvs+3fdf7dt19i32+P6HwnMbnrXPS3fcHEUvCDLURdYxtllinDnnRiQ9JGm48fe/4L3/TefcrZI+J+kmSY9I+oD3PuEtM33EGGAjxAAgX7FvO3WDGEOd5B1krdrZc7Yk6R3e+3nn3KCkbzjnvirp1yR93Hv/Oefcf5b0i5I+keFcX0aQATaCDAAKFd22U6eIMdRJkTG2WWKcee+9pPnG/w42/vGS3iHp3zT+/DOSPqoMVzAEGRBGjAFAPGLZduoGUYa6iCnIWrX1mzPnXL82dr/fLuk/STomadZ73/wi+ElJ+zOZoQgzYDNibEMVvzoEoBqK3nbqBEGGuog1yFq1FWfe+zVJ9znnpiV9UdJd7V6Bc+4BSQ9I0tDY9o4mR5QBG4ixoNJ/dQhANaW27TTa2bZTuwgy1EEZYmyzjl6Z3vtZSV+X9FZJ0865ZtwdkHQqsMwR7/1h7/3hgW1jbV8XYYa6mzjhX/4HW/MbQl8d+kLjzz8j6acKmB4A9LztNDjS/rZTO+b39xFmqLTxU+sv/1NG7RytcaekFe/9rHNum6R3SfrftLGiea82vjr0QUl/msaEiDLUFRHWnTS+OrQ+4jV3d/jQ8ZNPhY+vvrDTXmdd32Mf6ntt2H7cV3eEDyO+7fv2BtbAon3Z594UXv7Kbfah8Jen7Te9befMYS1tD99vSzP2Za+P2od9HzkVfrxWJhMOR7/dvu6hS/Z97ozFx87Zl33+jfbcpp+xH0/r1AoDc+aiGpi3n6cLM8b1TtvzXhm3L3vbeXNYi9Ph5a/tiTMy8t52agdBhiora4htpZ2vNe6V9JnGBlCfpM9777/snHta0uecc/+LpG9L+mQvEyHKUEcEWe+6/epQ69eG+ndMZzdBAHWUy7ZTO4gyVFWVgqxVO0dr/I6k12/x589LelOvEyDKUDcEWTYan1Df8NWhxt6zLb865L0/IumIJA0fOsCDAiA1WW87JSHIUFVVDbJWbR0QJG0EGeqEGMtOjF8dAoAiEGSoqjoEWatc42xtiDBD9RFjuYrmq0MAkDeCDFVVtyBrVcieM6BqCLJiFP3VIQAoCmGGKqlzjG1GnAFdIsgAAHkjylAVBNnWiDOgTcQYAKAoRBmqgCBLRpwBCYiy6huYd9r1UHh1OLAUfjO5vmJvME2css/JdfWgvRqe+kb4d7qX7jYX1fJ4wnnQroUvuy982jdJ0uRz9mVf32u/boauhK976LJ92bv+KuEca2evBsdOvHsyYVn7ulfH7Ns1v9+4XXP2b66HZs1hDc+tmePLxrnKlu2bLd9n367+pfDc14fsy975HfvJtDRpnwdt/kD4MRk7Xf0NvbVBwgzlR5S1jzgDtkCQAQAAdI8g6w5xBjQQZJ3j6KsAAKCJIOsdcYbaIsa6R5QBAACJIEsbcYZaIch6Q5QBAACCLDvEGSqPIEsHYQYAQD0RY+lo5+A+xBkqiSBLD1EGAED9EGTpaveoq8QZKoEYywZhBgBAvRBl6er0VBjEGUqNKMtOncJsbXpNc++ZC45Pf2E8OHbtgP0cXB+0V7NJ55/6vnEus5Gz9rLWObckyRvvF0vb7dvlVu3L7l+wx63LHzttL3v+sH2f7nhqLDyvRXNRLeyxN0qSbte2K+GxsfP2ecqu7bPfwC++zj4f2OTz4fv0pdvt23Xzn9njizvC173zMft2XfiBQXN8PWFLZPRs+HZdvrs+6ykgSwRZNro5RyFxhtIhyLJVpygDAKDOiLJs9HLieOIMpUCQ5YMwAwCg2giybPUSZhJxhogRZPkhygAAqCZiLB+9RlkTcYZoEGPFIMwAAKgWgixfaYWZRJyhYARZsQgzAACqgSDLX5pR1kScIXcEWfGIMgAAyo8gK04WYSYRZ8gJQRYPwgwAgPIiyIqVVZQ1EWfIFFEWF8Jsa31X+jX65xPB8ZVt4WVHE87JdX2f/RoYumIvP3IuvJrecdQ+v9TyuP0GsrQjfN1jp8xFtTaUMG7cZ5I0cD183YNz9n02cN2+7Iv3hs/JtXyrfaKz0adHzPHrB+z7fO628O3qW7Mfj+HL9nNh9Jx9v2y7sBocm37KfrtfH7Yve3U0PLfpZ+0H5PpM+LUlSZfvTXqNhO+3IeO8ckDVEGXFyzrMJOIMGSDI4kSYAQBQLgRZPPIIM4k4Q0oIsrgRZgAAlANBFp+8wkwiztAlYqwciDIAAOJHkMUpzyhrIs7QNoKsXAgzAADiRIzFr4gwk4gzGIix8iLMAACIC0FWDkVFWRNxhhsQZOVHmAEAEAeCrFyKDjOJOIMIsiohzLrUJ61uC99346fDh09f2h4+bLskuYT35YFr9vjIxfDrc2XUfhNZ2G0/H9aNd4DlKXNR9S/Y44NX7fG5O8L36fK0Pe+bHrfHD/5leHJn32of43/H0+HD0UvS1PP2fX71lvD49HNL5rIX7xk2x/tW7HX1S/cNBsdmnrBv1/w++3l87UB47Mod9qHyhy/aj9eubyXcrneHD9U/+MyouSxQBIKsnGIIM4k4qzWirFoIMwAAikGQlVssYSYRZ7VDkFUPUQYAQDGIsvKLKcwkKXE2zrmDzrmvO+eeds495Zz7lcaff9Q5d8o591jjn5/Ifrro1MQJf8M/qBbCDADiw7ZTdY2fWr/hH5TX/P6+6MJMam/P2aqkX/feP+qcm5D0iHPua42xj3vvfyu76aFTBFh9EGYAEC22nSqCAKumGKOsKTHOvPdnJJ1p/Pecc+6opP1ZTwztIcbqiTDb+GRa0mcl7ZbkJR3x3v+uc+6jkv6tpJcaf/Uj3vuvFDNLAHXEtlO5EWTVFnOYSR3+5sw5d0jS6yV9U9LbJH3IOffzkh7WxidEl7dY5gFJD0jSwNT2HqcLYqzeiLIb8Mk0gOj1uu00OMm2U9aIsXqIPcqa2p6lc25c0oOSftV7f1XSJyTdJuk+bXw69NtbLee9P+K9P+y9P9w/NpbClOuH34xBIsw2896f8d4/2vjvOUl8Mg0gKqlsO21j2ykL/G6sXsoSZlKbe86cc4PaWLn8gff+TyTJe3+uZfz3JX05kxnWFCGGVoSZrZtPplutD0rX94Vfc9eMlfrqzLI9twX7/FGTz9lvGBfeHD4f2NRT9ip8+rnwspJ07nD4um960l4HXb7bHJZbt5+zt3w5vEF09k3h83VJ0sTJRXP8+L8eCY7tetjeEJu9w75PR8/ay4+dDo9futs+j5lPeEcevG5f9/Bs+D4fWLCfC5fvtR+vgevh50q//XCY5+qTpKXJhPPa/blxbjof73sl205xIcLqq0xhJrV3tEYn6ZOSjnrvf6flz/e2/LWflvRk+tOrD46qiBDCzNbtJ9POuQeccw875x5eu5ZwJmgA6ADbTnFg7xjKFmZSe3vO3ibpA5KecM491vizj0h6v3PuPm38EP+4pF/OZIYVVrcImzq21NHfv3Kb/Slz1RFlyXr5ZNp7f0TSEUkaOXCwXi9GAFlj26kghBikckZZUztHa/yGpK22Ejn6WYeqHmOdxlenl1enWCPMklmfTDeOlCbxyTSAArDtlB9iDK3KHGVNHR2tEe2rcoilHWF4BVHWET6ZBoAaqUuITbyw0PNlzN1q/FazoqoQZhJxloqqhlhsETZ1bKnSe88Is87wyTQAVFuVYyyNAEu6/LoEWlWirIk461AVQyy2CLNUMdCIMgBA3VUxxLIOsHavv8qRVrUwk4gzEyGGLBFlAIC6qlqMFR1iliruRatilDURZy2IsXKowt4zwiwublUavhR+TOZvXQ2O7f2avRpduCnhDcTZ652bHg6fJ23H0evmssc/ZF/26D+GT267cJP9HN3xpL1htRo+1ZgkaW5/+H6becI+J9fcAfv1v+9vw8uvD9i3a+Zxe515/rB93Xv/LvyYzK/Zd4pbsx8v6z6TJG+cUu/6LvvccQf+2n48z/yz8Fj/gn2f+n77dl25yx7fdjb8Glq3TyOIiFQpxmIOsZCqBFqVo6yp1nFWtRirYoiFlDXQiDIAQB1UJcbKGGIhZQ+0OoSZVLM4I8ZQFKIMAFBlxFg5lDHQ6hJlTZWMs6pFmESItSrTHjOiDABQJVWJsKaqx9hmZQqzukVZU6njrIoR1kSMba0MYUaQAQDKrmoRJtUvxDYrQ5jVNchalSLOqhxhEiHWrtjDjCgDAJRNFSNMIsQ2iz3MiLJXRBVnVY8wiRDrRsxRRpABAMqgqhHWRIyFxRpmBNnWCouzOoSYRIz1KtYwI8qqxfdLK+Ph8d1/F34DmTtov7lMnkg4LPxB+1jg1gbVwu6Ew7r/ob2eHbgWXj+dfbN92UPz9obeyqh9u7Y/txgcG3zpmrns1bu3m+NDsyvBsVM/bG+kHPi6vc7uD09bkvTiu0aDY/v/q73w2lBvGyoTz80Hx469f9Jcdubb9nXPPBZ+Ls0dtNeHS9vt8aT7dOr58HPt8mvZuNsKIYamGMOMKLPlGmf9y9WOMkIsXbGFGUEGAIhN1UNMIsa6EVuUEWTti+prjWVDjFUfQQYAiAkxhjIhyjpHnHWAGMteDHvLCLL2VXlPOADEoH+l+kFGjKWn6D1mxFjviDNDGWNs6OiLXS23fPfBlGfSmaKjjCBrDzEGAOgFIZY+gqxaiLOGsoVYtxEWm6KijBhrDzEGAOgWIZatIqOMIOte0p7w2sZZWWKsKhG2Wd5RRoy1hxgDAHSrjDE2cOxU23939bb9Gc6kfUVEGTHWnW6+klyLOCtDiFU1wloRZHEiyAAAnSpTiHUSYJ1cTt6xlneUEWTd6fU3opWMs9hjLLYQy/r3ZnlFGTHWGaKshZPWh8P3x+BCeGzXt+31zfDZ8LmnJOnynTeZ40NXw+dJ619O+GrEXz1ijrs3vi44Nv2c/faw7aR9LrKB6yPm+OAz4Y2lpdfZ66TBefvccaffHt6AmTxu32erY/btvnLPqjk+8d3w8sc+YG/ojH13yBx3Ce/3fWvhk/UNzNvrx6u32uPrg+GxpHVJ34o9vjxrX/f1mfD9NnzZXBQpij3G0gqwWOURZsRY99I8aE8l4izmGIstxPKWdZgRZJ0hyAAA7Yg5xmIJsYFjpzLfe5Z1lBFk3cvqKKqljLNYY6zuIdYqyygjyDpHlAEALLHGWCwhlrcso4wg603Wp7YoRZzFGGNVCLG0v86YVZARY90jygAAW4kxxsoUYlnsMcsqyIixdOR1vsFo44wgy0YZgowYSwdhBgBoIsbSkXaUZRFkxFh6ijgBfFRxRpBlI/YgI8bSRZQBACSCLC2xBxkxlr4ioqyp8DgjyLKTZpSlGWTEWDaIMgCotxhjTCLIJIKsLIqMsqZC4iy2IKtKjDWlFWVpBRkxli2iDADqK8YgK2OMtUorzNIMMmIsezGEmZRznPUv+mjCjCALSyPKCLJ8EGbpGLqyrlu+HN7AuXxX+A125CV7Zb4+ap+7avp79vJjz75kjluu/cR/a1/20+eCY1MLK+ayi/sm7PEZ48RYkgYP7gqOrY3YGyFXXmNf9q5Hl4Nj5w7bj8f8AXt84rvmsHx/eGz6Efuyl7bbl+3s07vp9E+Hb3f/Sfu8c31X7XV2/2J4bHnSXnb/F+3326f/573m+I5HwpsqizfV870mtiAre4w1xRRlxFh+YomypsK/1pinqgVZUxph1muQEWP5I8wAoD4IsmylEWa9RhlBlq/Yoqyp0nFW1RhrKjrKigiypZvDnwy3Y/j79ifXZUCUAUB99C2tRxFmVYuxpqKjLLYgu37gxmAZPRnX/NISa5hJbcSZc+6gpM9K2i3JSzrivf9d59wOSX8s6ZCk45Le572/nN1U21PlIKt6jPUaXp1eRxlDLbYwi+VrygAQk7JtO4VUMcjqEGObAyvtyyp7sMUcZlJ7e85WJf269/5R59yEpEecc1+T9AuS/sp7/zHn3IclfVjSb2Q31WRVCrMYDuqRRYzlEWDtas6lLJFGmN2oKhs/ACqpNNtOW6lSlFU1xtIMsG6vu2yRFnuUNSXGmff+jKQzjf+ec84dlbRf0nsk3d/4a5+R9DcqcAVT9jCrWozFFGFJyhBpMYVZ0VHWotQbPwCqqyzbTlspc5jFcECPtGOsyAhLcv3AemkCrSxhJnX4mzPn3CFJr5f0TUm7GysfSTqrjU+v0YYYjqyYVoiVKcLKKKYok6IKs1Jv/ACojzJtO5UtzKoUYzFHmKUMgVamMJM6iDPn3LikByX9qvf+qnOvbNx7771zbsutSOfcA5IekKTh4eneZhsQ+16zKsRYHSJs6eblqPaexRRmMUXZVsq08QOgPtLYdhoZmspjqqUIszRPDN1tkPUaY2WNMEusgRZjlLVzcJ+24sw5N6iNlcsfeO//pPHH55xze733Z5xzeyWd32pZ7/0RSUckaXLiQOpbm7GGWdmDrA4xFjPCrH3pfHBkb/ysjIVfS+tD9hvS7B3j5vhNT8yb41d/MNyWQ1dWzWVHT1w1xy/88/CGztTz9hvItX32BxmD1+w3xeM/GT5P2o6n7ef/7Ovt9dPqO68HxxYumYtqeMJ+vq+u2G+ba9fC4/v/wn6uzB2yxydOJJyL7Nvhjc25u+zz1o1csM8d12csPnrOfqxPvvdmc/zQg/bjOXtH+HbvetS+Xc+Yo9lKbdtpfH/mbwgxh1nRQdZLjFUxxEJiC7TYwqyTI662c7RGJ+mTko5673+nZehLkj4o6WONf/9pZ9PsXUxhlmaMSd0FWbcxRoi9Ioa9Z7GEWexRJqX5wVH2Gz8A6iPmbafNYgyzMgZZnUIsdmUOM6m9PWdvk/QBSU845x5r/NlHtLFi+bxz7hclnZD0vo6uuUdFh1naMSblF2RZxNihAy+lcjnHT+5M5XJ6UWSgEWbtK9PGD4DaiXLbabOYwqzIICPG0hHD3rOYwqzb8xO2c7TGb0gKFcA7u7rWFCTFUTfxlkVwhZQxxNIKsKTriCHQikCYdawUGz8A6ifWbafNNgdRFrGWZnRtpax7xgb2hb9+3anV06OpXVaZxRJmvZ40vqOjNZZJnqFlyfP3YmmEWB4B1s4cigy0IvaaxRJmZVKWjR8AKIusQ6oXeR7AI60QSzPA2rmeOkdaLGGWhsrGWd56Ob+Y1HmM9RpiMUSYpahAI8wAAChWmUIsrwBrR9GRVsRXGqsUZU3EWYd6jTAp3xCLPcIseQcaYQYAQH56Ob+Y1HmM9RJiMUVYkiIijTBLD3EWkEaESfl9PbHMEWbJK9AIs3pbHe3ThfvCb2KjL4XfAJYn7dXotb32OmBtyD7U/sp4ePk937TXFc/8xpg5PvRC+LIvvHHEXNbZRzDX0JV+c3xpX3ju529bM5cdPGHPbW3K2EhYsjcg1i7Yj4cfTHjdjoSfK6fvtxfddtp+riwkrArHXwzPbeff28/ToTn7Pp87GH48p7/6tLns5G32zwwuvGHSXv5E+JQR5w7bpwDQV+1h5KcsIVamCEsysO96LoFGmKWLOBMhttm79mx9Zpivnb0r0+stSt5hFnuUXbltuEwHBQEARKTXCGvKOsaqFGFFIszSV6s4SyvCmvKIsV5DLBRaWV1WFgGX5d4zwgwAgO4UFWJSZzGWRojdvedcz5cRcvTs7swuO0t5h1nVo6ypknGWdoQ1ZR1j3YRYmvGVhtb5pBlqRR/BMQ2EGQCgrNIKMSnOGMsyvjq97jRjLauvNhJm2Yk+zrIKLUs3ESbls1csthgrszz3mpUtzKzXHV95BIBqSDO4QmIMManYGCs7wsy21euqk3Of5RpnayOukNgK6TbCJEIsSZp7zdLeY0aU9eZVr+G/KWQaAFAL68N9uURUr7qJMCmf34qVKcTS3GuW9h4zoqx7r3oN/0P470a/5ywNvUSYlM9BO8ocYpulFWZljjKpmmEGAKi3biNMIsSSpBVmRFm5VSbOeg2wJkKsN2mEGVEGAECxeokwiRDrVIxhxpEYi1G6OEsrwqRyhdh7Jx9N5XK+cPUNqVzOZjHuLSPK0K71Qen63vDjd/X28LITx+zzeS3usc8ftbzdXqdNvvZScOy526fMZfv77ZORrd85Hxz7oUMvmMv2yX6+v278tDn+azueD479+5fuMZf93NAbzfGBgfB9Pv6C/ba39IZr5nj/cXvDZ+B6+PmwvDt8vi5J8n32c2k4/FSQtPE8DhlYsB+v2dvt6x47E17+9C/cay47cN2+7j77btG1XeG5DV2xl0WyXiNM6u6cYp3GWJoh9o6ZzrbN/vpCtqcSSiPMiLLqiDrO6hRiacVXu9eTVaR1q8xRJhFmAIC4pRFhUvwh1ml4dXqZaYdar2FGlFVPFHGWZoQ1xRhjeQVYO9KMtF72mhFlAACkK60Qk7KPsW5CLIsA6/a6e4m1WMKMKItL7nEWS4hJncVYmUPM0mukxRBmfH0RAFBXZQoxqbMYKzLC2vWOmWcy/9rjZmWNMoKsPfkeSj+lbeg8Ykzq7muKZYmyoqURZkQZAKAu0oywpm5iTMr+a4pliLJedbvXLI0wI8riFsXXGpN0G2NSfr8bK2OU5bm3jL1k1ZPFXnAAwI3WBov9vVgTQfaKbveUEWRoR5RxlmeMSfUJsry/uljGvWMSQbYZEQYA5ZRnjEnl+/1YJwiyZARZsnY+aIkiznqJMYkgsxBk7SHINhBiAFBevcSYlN9RFQmyrfUaZMRYXLrd411InBURY1J9fkNWtiAjxopBiL2if0Ha8WT4ObH8b8InmJq/dpN94Ql382vf+H1zfO+2q8Gx7w3b69JtA/Z5zv7HQ18Njv3qd/47c9ldE+FzpEnS05ftjZNP/Nm/DI65hPf79ZsXzfGls+ENnL0n7Au/PDxmjvuEd809/xB+TOb3Gycik7Q23Nt6afR8+LYtTtsbCX32U0UTLy4Fx9YHhs1lF3bZL4JR4xxqkjR6IXzeupd+bsFcVv+nPVx2eceYVO29Y03dRFk3QVa2vWPEmC2trx/nGmd+yOd2MA+pPnvHmrqJsjoEGSG2gRgDgHIrIsYkgiykiCAjxuKQxQF6mqL4WuNW2DvWvjyirCwxRojdiCADgHJbH+K3Y1nKI8p6CbK8YowQS5ZlkLWKKs4Iss50GmV5BlkeMUaIvRoxBgCoSpD9xPjT+sr8PZlcdqdRlleQZR1jRFj78oqxzaKJszy/tlh2WUZZt0GWZYwRYcnyjrJefzcKAEhfnl9bzEMWYZZllMUWZIRY54oKslaFx1mee8tafeHqG0qx9yzWvWNZxRgh1r4i9pIRZQAQn26jTIo3zNIQ494xQiw+MQRZq8LirNsok9LbYxZToOV1hMWYYowQ6w5RBgCQeosyqfcwa8ZP0b83y+uQ97HEGCHWuyKCrN3fjuYeZ71EmVT+rzJ2G2FSeUOMCOsdQZat9e1ruvreueD49fOTwbF9b7Q3bkYH7WOUD/TZK+vh/tXg2I/s+a657P/79X9ujv/Ss78Ynte0/fhfe3CfOX59t/3GN2QcfX0gYXtzeWGbOb774fB97vvt19LECXs88bXowuNTz4cPRy9Jp98+Yl90+IjyjfHwfX7ltfZ6eG0i/DyTpL7V8NwWdtmXnTTvC2+yXwMT3w1vqixetJ8LVVV0lG2WV6R1G2FSOUOMCEtfzFHWlGucDQ3ZK/8kWYRZFnvPegmwpm4Ocd9NjKUZYkRY+ogyAEBTbFG2WRqR1kuANeV1iPu0YowIy14Zoqyp8N+ctSPrvWWdBloa8bVZHnvF0goxIiwfHOQDACD1HmVSvr8t2yrS0oiukLLsFSPC8pd3lPV6LkIp8jjL8yuMWQRXSB57xXoNMQKsOHlGGUEGAHFzg71t7BV5wI8sgiz2ECPA4lDGKGuKNs7K/tuyVlnvFeslxIiweBBlAIC0VOUojFnHWLchRoTFqcxR1pQYZ865T0l6t6Tz3uiNJ68AAB/dSURBVPt7G3/2UUn/VlLz6B4f8d5/JfXZlVSMMVZkhE0ds38E36krtxlHEigpwqw7rJ8AxCiGddPRs7tLGWhZxlgZQmzihYXcrkuS5m6tzkF1qhBmUnt7zj4t6fckfXbTn3/ce/9bqc+ooRk4Me9B6+briVJnMdbtXrG8Yizt8OrmOssea4RZTz6tAtZPAJDg04pg3dQMnRgjrZuDdkjZx1jWIZZ3fCWZeGGhEoFWlTCT2ogz7/1DzrlDmc0gQdGR1m2ANZV9r1gR8dWpqWNLpQw0DvjRu6LXTwCwldjWTZtDKK9Y6zbAmvL4imIWMRZbgCUpc6CV6SiM7erlN2cfcs79vKSHJf269/5ySnPakhVJnYZbr8EVUtYQK0OAJSlboBFmmeto/bS+0q+Fs+PB8b6V8ON15pj9uvcD9krcjdgngTo2flN4Xt+csq/7gH3ZQxf7g2PT/2C/Ua8mvI9PnrCve+pbp4Nj69Phx0KSfJ/9Zrw2Hl63rkzab3sjs/a8Z752whw/+TO3BccOfCF8myXp4GL4sZak82+w75ex8+G5X77XXueMnBk0xxd2ht97lvcnrF+W7cdr28mE637LteCYO1+6Dcpct51CrGiywq3X2LKULcTKFmBJyhhoVdpb1qrbOPuEpP8oyTf+/duS/vut/qJz7gFJD0jS8K6JLq/OllVsheR1GPs0Q6wKAZakLIFGmGWurfVT67qpf8d0nvMDUE9dbTsNzNgfwqQtywCT8jt6YlohVrUIs5Ql0Kq4t6xVV3HmvX/5YxXn3O9L+rLxd49IOiJJE3fuKd2hAcsUYnUIsCSxBxphlr1210+t66bhWw6Wbt0EoFy63XYauW1/addPZQmxOgVYktgDrephJnUZZ865vd77M43//WlJT6Y3pWJ0E2FNncZYLyFGgCWLNdAIs3xUcf0EoPyqvG7qJsKaOo2xXkKMCGtPrIFWhzCT2juU/h9Jul/SjHPupKTflHS/c+4+beyaPy7plzOcY+rKEGJEWG9iCzTCLBtVXD8BKL8qr5tiDzECLB2xBVpdwkxq72iN79/ijz+ZwVxS10uESdULsaGjL3a97PLdB1OcSfbqHGZ1Uub1E4DqqsK6qZcIk6oRYgPHTqV2Wau37U/tsrIWU5RJxYRZkXo5WmPheo2vzfL6rVgeMdZLiJVdLGFGlAEAeuVX+l4OpYF911O5zF7Dayt5/VYsyyBLM8bKKqYwq1uUNeUaZ8vLA6kHVS/yPHBHlkFW5xDbjDDbMPz9odp8tREA6iKLqOpWngfuyCrIiLFXxBRlUn3DTCr5nrNOdBNiTTHtHSsqxIaOvhj9VxsJM3Srf1GafDZ8zq/h2fA64Poe+w1k9Kw9fvnH7YiefDB8CpK+VXtDZ+/f2euh5anw+aX6lu3LHjlhn57Jraya44u37wqOXb3FXl9vP2rvPVjcGV5+ccp+PAaW7PX9qfeFz2MmSUNXwssf//lbzGVv/r+eMMfHDtxjjr/4k+HHbOJJ+z6df419frfhl8Kvj5ET9rrXJ6wSl6fs55o/F95o3Pv39mXbZ6VDu7oJsaZY9o4RYmGE2dZGT/bFeUCQMso7xJrSDjL2iLUnliiTCDMAQLn1EmJSHDFGiLWPMItPqeOslwhriuFrikRY9wgzAAC6U0SISenGGCHWndiiTCLMmkoTZ2mEmFR8jJU5xGL7aiNhZuN3ZwAAqfcIayo6xgixdMQWZkTZjaKLs7QirKnIGCtziMWOMAMA4NXSCjGpuBirWoTFchj92KJMij/MivjdWWFxlnaENRUVY1UOsVj2lsUUZFI5ooy9ZwBQPWkGWCtiLD2xBJkUZ5RJ8YdZU96BlmucuWUXzZ4xib1jlliCTIovyqRyhFnTVq85gg0AyqFvObsYk7oPMimdrytWKcxiCjIp3iiTyhNmTZtfg1nGWnRfa2xHL0EmEWUhBFl7yhRmIVntuQYAxK/oIJOqE2UEWWfKFmUhWX5gUoo46zXGJA7oERJTkElxR5lUjTDDq60NSfMHw+uZteHw4748Ya+f+hcSnjMn7JPK7vj78AbM5bfsM5ftW7bPXTVyPnwusrP/bNJcdv+z9t7XpdvC5zGTpLWR8Bvb4DX7Pl0btd+6Lt4dPifXzsft86/NHQwvK0m7HrHPsfb9fxneMOpP2KY9+cs/YI6vhE95t8E4791qwvbattP27R66Eh6bu9Xe0B+4br8GfMKWyMC18PJrQ71vH9RBLzHWRJS9giDrTlXCLGtRxlkaMSZxuPsQgqxzRBkAoCzSiDGJIGtFkHWPKOtMNHFGkGUrtiCTyhFlEmEGAIhbWjEmcQ6yVgRZb4iy7hQSZ2mFWBNBFhZblJUlyCSiDAAQpzRjTCLIWhFk6SDMupdrnPUvs4csDwRZb4gyAEBMYo4xiSBLG0FWb9F8rbEdBNmrxRZiUvlirBVhBgCISf9KOpdDkG2ILcSk8sZYK8IsPaWIszSjrOxijDGp3EEmEWUAgOpKO8zKJMYYk6oRZBJRloWo4yyrKCvLXrNYQ0wqf4w1EWWdS/s3ozHoX5SmnwmPrxhHlb/lq/ah1Zen7XPKbbtoH8J86daZ4NjAov1Vp5deP26O7/mzE8Ex32cfSv/C/fb6acd3Zs3x2dung2NjZ+1TAFzfPWiOrxinN7jyGvttb/5m+z4dO2ev+6zrHrxqr2/WRsxhDSRsX7uT4ftlaca+XXf8oP2+eOKhW4Jj+/7Wvuxru+3n+OqYvXHXZ+w5Gpq3nyt1l1WUxbrXLNYQk6oTY01EWfeSvqYcZZxluacs1jCLOcSk6sRYE1HWmSoGGQBUVR32lMUcYlL1YkwiyHrRye9Go4qzunx9MfYQk6oXYxJB1gliDADKJ48oy3uvWewRJlUzxJoIsu70chCfwuMszyArYq8ZIVYsgqw9xBgAlFdV9pSVIcQkYgxbS+uoqoXFWd57yfIKM2KsWMRY+wgyACivIoIszb1mhFgciLHupX2Ki6Z8z3O26Av56mIWYVaGCJOqHWISMdYJYgwAyq3oPWTNoGo30soSYK2IMYRkFWObFf61xjwkhVRrvJUlukKqHmMSQdYuYgwAqqHoKNusjNEVUvUYkwiyXuQVZK1qEWdJyhxkdYgxiSDrBFEGAMCr1SHEJGIsDUVEWRNxVjLEGLZCkPXGeal/OTw++VR48OIPjJqXveMp+xPv1VH7TfTSXeHX/Ng5+83j8g/Y43v+MnwOtj3/OG8uO/vaMXP8xE9uN8e98RK/cod9XqyBa/b6YfL58Ovh0htWzWX75+zrPvmv7OWnHwnfp9d/2L5PZx60n0uL2+3bffXW8Fj/or3s8/94szk+bJzOb3Havs9Wxu3rHj9tP0+HL4fv8xfflbAZ80V7uAz6lorbSCy7OsQYIZaeIoOsFXEWMUIMFoIMAIBXVD3GCLH0xRJkrYiziNQhxgix3hBkAAAQYuhejEHWijgrSNVDjAhLD0EGAKizKocYEZaf2KOsiTjLQVVDjADLVuxRVsRpMTZzzn1K0rslnffe39v4sx2S/ljSIUnHJb3Pe3+5qDkCADpTtRgjwIpTliBrRZxloIoxRojlI/Yga4ohzBo+Len3JH225c8+LOmvvPcfc859uPH/v1HA3AAACQgxpK2MQdYqMc74ZDoZMYZelSXKpKjCTN77h5xzhzb98Xsk3d/4789I+hsRZwByxLZTGDGGrJQ9ypra2XP2afHJ9MsIMaSpTFFWIru992ca/31W0u4iJwOglj4ttp0kVSvGCLE4VSXKmhLjrO6fTBNjyAJRlg/vvXfObXlnO+cekPSAJA2Ob9ea8VJ/4afCb8i3/pcVcw5n32Kfu2rq+Jo5PvN4+ARTJ99hn2ts5Kw5rKf/p5ng2OC5QXPZPf9kz/ulGXsj5qbHw+uh2deai8olvA8vT4Yve3DWPifX1LP2ZU+9YF/58XeHxwePjpvL9i/Z51Db8bT9XLt0X/gx2/mP9u0emk8419hseG5XXhM+t5skXd9rr+/2/e01c3xxZiQ4tvfv7ct+wRzNTl23naoUYhIxVgZVCzOp+9+cVfKTaUIMeSDMMnfOObfXe3/GObdX0vmt/pL3/oikI5I0uvMgDwqArFVq24kQQ5GqGGVNPR8QxPpkWrrx0+nh4elery41hBjyRpTl5kuSPijpY41//2mx0wGAG3Wy7TQyNJXbvEIIMcSiylHW1G2ctfXJtHTjp9OTEwcK2TolxFAkoiw7zrk/0sbXhGaccycl/aY2ouzzzrlflHRC0vuKmyEAvKy7bafx/bm9iVQtwiRCrCrqEGVN3cZZlJ9ME2GICVGWPe/9+wND78x1IgCQLJptpypGmESIVVWdwkxq71D60X0yXcUIayLGqqHKYXbltuGoDqcPALGJYdupqgHWihirtrpFWVM7R2ss5JPpKgfYZgRZdVQ5ylrd8Pr8m8KmAQBRSnPbaX24rxah1S6CrPrqGmVNPR8QpBNrI65W0RVCjFVTXcIMAIA8EWT1Ufcwk3KOszojyKqNMCs3J/vcWXd+MnweprVx+xxPi7vs58aeby2b4yffGT6X2cx37PNiOftUZLp0OHwOtbG/njSXvbbbPm/Wjifs6/a9bGslrE4P/Pml4NipH91hLjtxyn48jv2sfbv7xxaDY+N/ae/9OPO2hHOwfc8eH38hfMf0rdkbPJfuTng8nwmPrQ3ZD8j0s/Zr4Pi/ts//tu1s+PLnDyWsex+0hxEnYqyeCLMNxFlGiLH6IMwAAOgNQVZfRNmNiLOUEGP1Q5QBANAdYgwSYbYV4qxLxFi9EWYAALSPGMNmhNnWiLM2EWNoIswAALARYwghymzEWQAxhs2IMgAAtkaMoR2EWTLirAVBhhDCDACAVyPK0C7CrD3EmYgy2Aiz6vOSvHEk8dm7JoJj4wmHXp88Zl/3yfvtw6svT4XfzC7eY6/C17bZz939vxe+7tNvsze4dhy1j9N/9q32enX0bPjyB8JH+JckbTtn365n/l34NACTxiHhJen4v7Lv05u+ad+u4avhJ9Lsbfay60P2hsvsa+3l+xfD94vbZT+ee/7Jfh7P7R8Mjq2OmotK3p73zOP2c+nyneH7dPQM799FIMjQKcKsfbWMM2IM7SLMAAAgyNAdoqxztYkzggydIswAAHVGkKEXhFl3Kh1nBBm6RZgBAOqIIEMaCLPuVSrOiDH0iigDANQJMYa0EWa9KXWcEWNIE2EGAKg6YgxZIcrSUbo4I8iQBcIMAFBVBBmyRpilpzRxRpQhK4QZAKCKiDLkgTBLV9RxRpAhS0QZmvpWvUZfCp9rab0/vC568V8Mm5eddB6mkQv23OR72LhKOL/UqR8Kz71vxb7olVF7XsOX7OteNU7vlnR+tvHT9nmxLt0XHlsJnwJNkjT1rH27Zh67ao7P3jkeHBuetW/XyoR9n02+YA5LxvaRH7Cve30g4Xl6JXzh/StJj7U9fm2PcZJBSX3GKdjYKHw1ggx54fWXjejijCBDHggzAECVEGXIE2GWnSjijCBDnggzAEAVEGQoAmGWrcLijCBDEQgzAEBZEWMoElGWj1zjbG2IKENxCDMAQNmsDRJlKB5hlp8ovtYIZIkoAwAA6BxRlj8+ikGlEWYAAACdI8yKwZ4zVBJRlg++pgwAQLUQZcUizlA5hFm2qhhkKxNOp344/EWCvuXwbZ5+1r7sa/vt8eFL9vj6cPj5vPuRVXPZtWH7yxEL28Pjl19nv45GLtrPA2+fukojF8NjO542Tmwl6cprhszxIeOy+5fMRTX1gn2Ct+/9bPg8ZpL0mgcXg2OX7jFO7iZp5IJ9ny7MmMMaPRd+zEYu2htbp+63Nwdu/8MrwbGL903ZE0vYzluatsfHT4Vv1+pI9dZHQJEIs+y0+9tR4gyVQphlo4pBBgAANhBl2ejmYD7EGSqDMEsXQQYAQPURZunq9eiqxBlKjyhLD0EGAEA9EGXpSfN0F8QZSo0w602VY8w5d1zSnKQ1Save+8PFzggAgOIRZenI6vyDxBlKiSjrXpWDbAs/4r2/UPQkAAAoGlHWm7xOBt9TnPHJNPJGlHWnZkEGANFi2wl5I8p6k1eUNaWx54xPppELwqwzBJm8pL9wznlJ/7f3/oj1N91a+IKmvhseu3if/aY3ct4+pvz8wYRD1huHV79056C5rE9YwzvjSPyTz9vPn/4le94DC/Z1e+O97vhP2rdr8jl7btZh+kcu2vO+fJd93TuesJc/+S9Gw9ed8E65MmGPTx2zn2vn3xqe2/bH7Y2L6aP2dZ9/U/hw+UmnTVgbsh+v4Vl7+cUd4eVHz5V2o5NtJ+SCMOtO3kHWiq81InpEWfsIshu83Xt/yjm3S9LXnHPPeO8fag465x6Q9IAk9W/fXtQcAQBIHVHWuSKDrFWvcdb+J9NAFwizZATZ1rz3pxr/Pu+c+6KkN0l6qGX8iKQjkjR8MGH3FQCkh20nZIYo60wsQdaq1zgzP5mWbvx0emCKT6fRHqLMRpDZnHNjkvq893ON//5RSf+h4GkBgNThttPgJNtOSEaUtS/GIGvVU5wlfTLdGHv50+mR/Xw6DRtRFkaQdWS3pC8656SN9dwfeu//vNgpAUDn207b9rDthDCirD2xB1mrruOMT6aRNsLsRsRY97z3z0v6waLnAQCt2HZCWoiyZGUKsla97Dnjk2n0jCB7BTEGAJXHthO6RpAlK2uQteo6zvhkGr0gyjYQZABQH2w7oVMEWbIqBFkrDqWPXBFlBFmM+pekyWPhx2XyxZXg2MU32m8Kq6P2c37ihD234dnwG/PKmP1cWpq2x61zuyWdf+3AX4fvE0ma+7Elc3z8i2PBsf5Fe94LO81hbbsQnrtbt2/X8lTCObku2dc9+Xz48pMej8E5+7LPvd3eSBv9fvgtfXnavuyx0/ZlL0+Hn+dJ545za/b49X32/TIwHx6bu7laG2VAE1Fmq1qQtSLOkJs6hxlBBgAAkhBlYVUOslbEGTJXxygjxgAAQLuIsq3VJchaEWfITN2ijCADAACdIMperY5B1oo4Q+rqEmXEGAAA6AZR9oq6x9hmxBlSU/UoI8YAAEC3CLJXEGRhxBl6UuUgI8YAAEC3iLFXEGPtI87QsaoGGTEGAAB6QZC9giDrDnGGtlUtyogxNPUvek1/bzk4fm3vUHBs2yn7ebQ6Zr9utl2w38hnb+8Pjg0lnBdr6Sb7useMuY+8ZN+ul+4L3yeSpGcGzeFre8OXvzxlX/TOx4wTtEm6eG/4PlvYZd+u6Wftx2Mp4TxoAwvhsdHz9mWvDtuX3f8d+y17cSY8lnQOtZHL9n06dzC8kbU8nrAuTdg+s863J0lDV8PP4/mb7WWBPBBkryDIekecIVFVoowYAwAAaSHKiLEsEGcIKnOUEWIAACBtBBlBljXiDK9SxigjxgAAQFbqHGXEWL6IM7ysTFFGjAEAgCwRZCgCcVZzZQkyYgwAAGSJGEMMiLOaIcYAAADqHWMSQRYr4qziyhBjhBgAAMgaMUaMlQFxVlExRhkRhlgt75BeeH/4+TnxVHhseNa+7LURe/z0j9snedr2QvicXda52SRpaco+F9mVO8PX3X/NfhPflnAetJ2P2RtBl39uPjg29A37RGcLO+y57Xp4JTh2+ofst72k85itTNjjl+9bDY5NP5nwlpuw3Tj7Bvvx7psLX/7UMfuy5w7Yc1uett5T7Ptk/KT9frS6zV5+6nj48VzannC+PdRenYOMGCsn4qyCYggzQgwAABSljlFGjFUDcVYhRUQZEQYAAGJQtyAjxqqJOKuAPKOMGAMAALGoU5ARY/VAnJUQMQYAAOqIGEPVEWclkfdXFokyAABQtDrFmESQgTiLUlEH9CDIAABAkeoWYxJBhhsRZwUiwoA49C06TTwVPiT39b3h1+r0swmXvWK/3nZ9fdAcX9xhjN1kr8IXd9kbOdufCG8QLE/a8x65YK+/BuftUwSsPxo+XP6gfcR4rY7Zc1vbFr5d64MJh3UftTeSlrbby297MfyYzN4bPsy+JE181348x5+xDxvfvxQeO/9m+/GY/J59u9dGws+l1TFzUV3fbT9ePmG79Nzh8O3uCx9lHyVQxxCTiDHYiLOCZB1mBBgAAIhN3YKMEEOniLOcZRFlhBgAAIhZXaKMGEOviLMcpRFmhBgAACiDOgQZMYa0EWcZSmsvGUEGdM4592OSfldSv6T/x3v/sYKnBACVVocYkwgyZIs4S1HaX1kkyoDuOOf6Jf0nSe+SdFLSt5xzX/LeP13szACgOuoSYxJBhvwQZ13gYB5A9N4k6Tnv/fOS5Jz7nKT3SCLOAKBLxBiQPeLMUMSh7gkzIBX7Jb3Y8v8nJb25oLkAQCnVKcaaiDIUzXnffYB0+puOw4cP+4cffrjr6wMQH+fcI977w0XPo5Vz7r2Sfsx7/0uN//+ApDd77z/U8ncekPSAJN18881vPHHiRCFzBZCdSNdPbDsBNWetm7r+eKDlNx0/LukeSe93zt3T7eUBQIpOSTrY8v8HGn/2Mu/9Ee/9Ye/94Z07d+Y6OQD1xLYTgCS97Lt9+Tcd3vtlSc3fdABA0b4l6Q7n3K3OuSFJPyPpSwXPCQDYdgJg6iXOtvpNx/7epgMAvfPer0r6kKT/T9JRSZ/33j9V7KwAgG0nALbMDwiy6XcdWV8dAEiSvPdfkfSVoucBAJ1i2wmor172nCX+pkPidx0AAAANbDsBMPUSZ/ymAwAAoH1sOwEwdf21Ru/9qnOu+ZuOfkmf4jcdAAAAW2PbCUCSnn5zxm86AAAA2se2EwALp0EHAAAAgAgQZwAAAAAQAeIMAAAAACJAnAEAAABABJz3Pr8rc+4lSSdyu8JkM5IuFD2JgJjnJjG/XsU8v07ndov3vtQn4nHOzUl6tuh5BMT6XIl1XlK8c4t1XlK8c+t1XlVYP7Ht1L6Y5yYxv17EPDcpxW2nXOMsNs65h733h4uex1ZinpvE/HoV8/xinltWYr7Nsc4t1nlJ8c4t1nlJ8c4t1nnVWcyPScxzk5hfL2Kem5Tu/PhaIwAAAABEgDgDAAAAgAjUPc6OFD0BQ8xzk5hfr2KeX8xzy0rMtznWucU6LyneucU6LyneucU6rzqL+TGJeW4S8+tFzHOTUpxfrX9zBgAAAACxqPueMwAAAACIQi3jzDl33Dn3hHPuMefcwxHM51POufPOuSdb/myHc+5rzrnvNf69PbL5fdQ5d6pxHz7mnPuJguZ20Dn3defc0865p5xzv9L48yjuP2N+sdx/I865f3LOPd6Y379v/PmtzrlvOueec879sXNuqIj55cE592POuWcbt/XDRc+nKab1VMzrqFjXT7Gum2JeJ7E+iltM66TGfKJdLxnzK/x11phHlOunNuYXy/2X6bqqll9rdM4dl3TYex/F+RKccz8kaV7SZ7339zb+7H+XdMl7/7HGBuN27/1vRDS/j0qa997/VhFzapnbXkl7vfePOucmJD0i6ack/YIiuP+M+b1Pcdx/TtKY937eOTco6RuSfkXSr0n6E+/955xz/1nS4977TxQ51yw45/olfVfSuySdlPQtSe/33j9d6MQU13oq5nVUrOunWNdNMa+T6r4+il1M6yQp7vWSMb+PKo73/ijXT23Mr/D1VGN+ma6rarnnLDbe+4ckXdr0x++R9JnGf39GG0/KQgTmFwXv/Rnv/aON/56TdFTSfkVy/xnzi4LfMN/438HGP17SOyR9ofHnhT7/MvYmSc9575/33i9L+pw2njtoEfM6Ktb1U6zrppjXSayP0ImY10tSvOsmKd71Uxvzi0LW66q6xpmX9BfOuUeccw8UPZmA3d77M43/Pitpd5GTCfiQc+47jV33hX11oMk5d0jS6yV9UxHef5vmJ0Vy/znn+p1zj0k6L+lrko5JmvXerzb+yklFtFJM2X5JL7b8f0y3Nfb1VHSvsU2ieH1J8a6bYlwn1Xx9FLvY10lSRK8vQ+Gvs1axrp+aYlxPNeaV2bqqrnH2du/9GyT9uKR/19j1HC2/8d3T2L5/+glJt0m6T9IZSb9d5GScc+OSHpT0q977q61jMdx/W8wvmvvPe7/mvb9P0gFt7Em6q6i54AalWU/F8BrbJJrXV6zrpljXSayPolaadZIU5XpJiuR11hTr+qkp1vWUlO26qpZx5r0/1fj3eUlf1MadGptzje/cNr97e77g+dzAe3+u8cRcl/T7KvA+bHzf90FJf+C9/5PGH0dz/201v5juvybv/aykr0t6q6Rp59xAY+iApFOFTSxbpyQdbPn/aG5rCdZT0bzGNovl9RXruqkM66Saro+iVoJ1khTB68sS0+ss1vWTNb+Y7r+mLNZVtYsz59xY48eFcs6NSfpRSU/aSxXiS5I+2PjvD0r60wLn8irNF2/DT6ug+7Dxo8xPSjrqvf+dlqEo7r/Q/CK6/3Y656Yb/71NGwfGOKqNFc17G38tuudfir4l6Q63cYSlIUk/o43nTqFKsp6K4jW2lRheX7Gum2JeJ7E+ildJ1klSxOslKY7XWWMeUa6fmmJeTzXmkem6qnZHa3TOvUYbn/hI0oCkP/Te/68FTknOuT+SdL+kGUnnJP2mpP8i6fOSbpZ0QtL7vPeF/LA0ML/7tbFb2Us6LumXW76nnOfc3i7pbyU9IWm98ccf0cZ3kwu//4z5vV9x3H//jTZ+tNqvjQ9rPu+9/w+N18nnJO2Q9G1JP+e9X8p7fnlwG4fi/T+0cR98quj1gRTfeirmdVSs66dY100xr5NYH8UrtnWSFPd6yZjf/YrjvT/K9VMb8yt8PdWYX6brqtrFGQAAAADEqHZfawQAAACAGBFnAAAAABAB4gwAAAAAIkCcAQAAAEAEiDMAAAAAiABxBgAAAAARIM4AAAAAIALEGQAAAABE4P8HLoXYar5pziIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def twoD_Gaussian(xy, x_pos, y_pos, theta, sf_0):\n",
    "    FT_lg = lg.loggabor(x_pos, y_pos, sf_0=np.absolute(sf_0), B_sf=lg.pe.B_sf, theta=theta, B_theta=lg.pe.B_theta)\n",
    "    return lg.invert(FT_lg).ravel()\n",
    "\n",
    "# Create x and y indices\n",
    "x = np.arange(lg.pe.N_X)\n",
    "y = np.arange(lg.pe.N_Y)\n",
    "x, y = xy =  np.meshgrid(x, y)\n",
    "\n",
    "#create data\n",
    "x_pos, y_pos, theta, sf_0  = 14.6, 8.5, 12 * np.pi / 180., .1\n",
    "data = twoD_Gaussian(xy, x_pos, y_pos, theta=theta, sf_0=sf_0)\n",
    "\n",
    "\n",
    "# plot twoD_Gaussian data generated above\n",
    "#plt.figure()\n",
    "#plt.imshow(data.reshape(lg.pe.N_X, lg.pe.N_Y))\n",
    "#plt.colorbar()\n",
    "\n",
    "# add some noise to the data and try to fit the data generated beforehand\n",
    "data /= np.abs(data).max()\n",
    "data_noisy = data + .25*np.random.normal(size=data.shape)\n",
    "# getting best match\n",
    "C = lg.linear_pyramid(data_noisy.reshape(lg.pe.N_X, lg.pe.N_Y))\n",
    "idx = lg.argmax(C)\n",
    "initial_guess = [idx[0], idx[1], lg.theta[idx[2]], lg.sf_0[idx[3]]]\n",
    "print ('initial_guess :', initial_guess, ', idx :', idx)\n",
    "\n",
    "import scipy.optimize as opt\n",
    "\n",
    "popt, pcov = opt.curve_fit(twoD_Gaussian, xy, data_noisy, p0=initial_guess)\n",
    "\n",
    "data_fitted = twoD_Gaussian(xy, *popt)\n",
    "\n",
    "extent = (0, lg.pe.N_X, 0, lg.pe.N_Y)\n",
    "print ('popt :', popt, ', true : ', x_pos, y_pos, theta, sf_0)\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5))\n",
    "_ = axs[0].contourf(data.reshape(lg.pe.N_X, lg.pe.N_Y), 8, extent=extent, cmap=plt.cm.viridis, origin='upper')\n",
    "_ = axs[1].imshow(data_noisy.reshape(lg.pe.N_X, lg.pe.N_Y), cmap=plt.cm.viridis, extent=extent)\n",
    "_ = axs[2].contourf(data_fitted.reshape(lg.pe.N_X, lg.pe.N_Y), 8, extent=extent, cmap=plt.cm.viridis, origin='upper')\n",
    "for ax in axs: ax.axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to top](#LogGabor-user-guide)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## performing a fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:45.394690Z",
     "start_time": "2018-06-20T15:16:44.506687Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from LogGabor import LogGaborFit\n",
    "lg = LogGaborFit(parameterfile)\n",
    "lg.set_size((32, 32))\n",
    "\n",
    "x_pos, y_pos, theta, sf_0  = 14.6, 8.5, 12 * np.pi / 180., .1\n",
    "data = lg.invert(lg.loggabor(x_pos, y_pos, sf_0=np.absolute(sf_0), B_sf=lg.pe.B_sf, theta=theta, B_theta=lg.pe.B_theta))\n",
    "data /= np.abs(data).max()\n",
    "data_noisy = data + .25*np.random.normal(size=data.shape)\n",
    "\n",
    "\n",
    "data_fitted, params = lg.LogGaborFit(data_noisy.reshape(lg.pe.N_X, lg.pe.N_Y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:45.418079Z",
     "start_time": "2018-06-20T15:16:45.396784Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(32, 32)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_fitted.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:45.442571Z",
     "start_time": "2018-06-20T15:16:45.421098Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Name        Value      Min      Max   Stderr     Vary     Expr Brute_Step\n",
      "B_sf       0.3755    0.001      inf  0.01536     True     None     None\n",
      "B_theta    0.1684    0.001      inf 0.007364     True     None     None\n",
      "phase     -0.2989     -inf      inf   0.1059     True     None     None\n",
      "sf_0       0.1011    0.001      inf 0.001795     True     None     None\n",
      "theta      0.2045   -1.571    1.571 0.005312     True     None     None\n",
      "x_pos       14.31        0       32   0.1593     True     None     None\n",
      "y_pos       7.726        0       32   0.2705     True     None     None\n"
     ]
    }
   ],
   "source": [
    "params.pretty_print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:45.865095Z",
     "start_time": "2018-06-20T15:16:45.445028Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "params : Parameters([('x_pos', <Parameter 'x_pos', value=14.305894604633949 +/- 0.159, bounds=[0:32]>), ('y_pos', <Parameter 'y_pos', value=7.725755384263746 +/- 0.271, bounds=[0:32]>), ('theta', <Parameter 'theta', value=0.20451580315562445 +/- 0.00531, bounds=[-1.5707963267948966:1.5707963267948966]>), ('sf_0', <Parameter 'sf_0', value=0.10113895438265275 +/- 0.00179, bounds=[0.001:inf]>), ('phase', <Parameter 'phase', value=-0.29890845891507284 +/- 0.106, bounds=[-inf:inf]>), ('B_sf', <Parameter 'B_sf', value=0.37545690182138414 +/- 0.0154, bounds=[0.001:inf]>), ('B_theta', <Parameter 'B_theta', value=0.16840243094402985 +/- 0.00736, bounds=[0.001:inf]>)]) , true :  14.6 8.5 0.20943951023931953 0.1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAEvCAYAAADB37lNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dabBtZ33f+d9z5vnOuvOgERACCSELE2EbC9sFMQG7QtEmHULKdkSqQmJX/MIU3dUmcbuKdNmmnU6Mcx0IkDgMzRAUl9yxcEvNYDNcyZoFmtCd56t7zzw//eLsjfY9Os9/T2t41lrfT5VK0n3u2vs5++yzzvrutfd6nPdeAAAAAIB89eQ9AQAAAAAAcQYAAAAAUSDOAAAAACACxBkAAAAARIA4AwAAAIAIEGcAAAAAEIG+LO9s+/bt/tChQ1neJYCUPfTQQxe89zvynkc32DcB5cT+CUCMrH1TpnF26NAhHTlyJMu7BJAy59zRvOfQLfZNQDmxfwIQI2vfxNsaAQAAACACxBkAAAAARIA4AwAAAIAIEGcAAAAAEAHiDAAAAAAiQJwBAAAAQASIMwAAAACIAHEGAAAAABEgzgAAAAAgAsQZAAAAAESgL487vePX/zCPu+3Y1EGX9xQApOjZJ0/q7a/7X8Pj/8tweOPjxpik7Y96c/zM3cvm+PhTA8GxHnvTpqb3h+e2smXJ3NbN95rjw6fscRkPS9+cvenAFfsxndkb3mePnLW3fek19riuWTCHRx4JPx8WN9m3vbhjxRwfPmn/yu6fDo/N7rbve+SU/Xtu6vb54FjvmUFz2x77IVPPsn3fK0NNvicl99Sxs7rtn38872m8wuyuvGcAZGPkTLb3l2mcPXH6rF71ux+XDjqNH413Z0uMAQAArCHEgOxkGme9i4oyyogxAABQdUQYsCbrs2WNcnlbYwwIMgAAUEVEGLCxPKOsrlJxRpAByYrxTDgA4GWEGNBcDFFWV4k4I8qA7hFiABA3QgxoXUxB1qjUcUaUAZ0jxgAgTkQY0LlYo6yutHFGmAGtI8QAIE6r/cQYkITYo6yudHFGlAHNEWNXm9/Rpx98cHNw/MadJ4NjLxzbb972mZ+x164aPBlex0ySpm4IL2a243v2WmKv+6ePm+P/3zdfFxz7iVf/yNz22J/caI5P7zWHzTW5rrzGfsw2P9ljjjtj80uvXzW3Hd8/aY4P9Nlzm+0Pr3PWs2L/fuq7Yn8/x47bP7fzW8O3/+afetLc9m+/eos57hfCc+tZNDfVyGn7657/2SlzfOli+DHd9GTpDmMAJKwoUVbXdK/mnBuS9A1Jg7W//yXv/e84566V9HlJ2yQ9JOn93vsmu+jkEWOAjRADgGzFfuwEVEHRoqzOfvlxzYKku733t0q6TdLbnXM/KenfSPq49/4GSS9J+rX0pnm1qYPux/8AeKXxo/7H/wAAMhfdsRNQBSNnXv6nqJrGmV9TfwNKf+0fL+luSV+q/flnJP1SKjOsIciAsMYYI8gAIF+xHDsBVVH0IGvU0pu1nXO9Wjv9foOkfy/peUmXvff1D0OckNTkEwadI8iAqxFga3jrEIBY5X3sBJRdWWJsvZbizHu/Iuk259xmSV+V9OpW78A5d4+keyRpYHRLW5MjyoA1xFhQ/a1D0865fknfcs79haR/qbW3Dn3eOfcnWnvr0CfynCiAaknq2KlvU3vHTkDZlTXK6lr5zNmPee8vS3pA0pslbXbO1eNun6QNL2fmvT/svb/De39H3/Boy/dFmKHqeJtic7x1CEDsuj126h1t/dgJKLMyvXXR0srVGndIWvLeX3bODUv6ea19oPUBSe/R2luHPiDpa0lMiChDVRFhnUnsrUM+vO85cXlTcGxg0t5n7f2mfbcLm+xLu58fC7+Gtjhh3/cTfxy+VL4krd4dfqfnkRcOmtsOXm+/trew1f66dv21se1p+5Ly0z89Y44PD4e/Ln9uzNx2/L9MmOP9H7SPDC4cCC990Ddpf13LE/Zl+qcO2b+yx3/ifHDsu//va81te/vNYblZ41L6rzLWRZDU/5z9mA9/2Y6PVeNhu3yTuWlusj52AsqqCjG2Xitva9wt6TO1A6AeSV/03v+5c+4pSZ93zv3vkv5W0ie7mQhRhioiyLrX6VuHGt821LuFtw0BSFQmx05AWVUxyuqaxpn3/jFJb9jgz1+QdGe3EyDKUDUEWTpqr1Bf9dah2tmzDd865L0/LOmwJA0e2M83BUBi0j52AsqqylFW19IFQZJGkKFKiLH08NYhAACKjSC7WqZxtjJAmKH8iLFM8dYhAAAKhiALy+XMGVA2BFk+eOsQAADFQJC1hjgDOkSQAQAA2Iiy9hBnQIuIMQAAgOYIss4RZ0ATRFkFeKknvDyVZs+H12EaX7Bv+tgv2GtbbXvU/hzuxLPhMbdqPzfP3WV8UZJcj7H9ZXvhq4Vr7S+8f3jJHJ/bHl77auSM/XVdumbYHJ8cHAzP67L9/bh8g/39WPrOHnN87Ep4+4FJ++u6/Bp77bixJvuiKb8jOLYybq8713OTvVbZ+LfHg2PLZ8JjkrTaZ8/7zJvNYfXtnA3f9wKHMUBMiLLusVcDNkCQAQAAtIYoSw5xBtQQZO3j6qsAAFQTQZYO4gyVRYx1jigDAKB6CLL0EWeoFIKsO0QZAADVQpBlizhD6RFkySDMAACoBoIsP8QZSokgSw5RBgBANRBl+SPOUArEWDoIMwAAyo0giwtxhkIjytJTqTDrkVaHws+lLY+G18ZaGbJvevPT9uM4eZ29/fbHwvM6dbe9dpV8k+/hpLGW2WZ7nbLr9543x099fb991zeG5z58xl7v6/pbTpjjF//vfeH7vd7ctOlD1j9t/4XZveGva+ykfdtqsju7eLv9/d5+JPy4LWyxH9PJTfYTefO58H1P7bdve+T8ijl+MbwsnSRp5cxIcGziefu+AYQRZXEizlA4BFm6KhVlAABUCEEWP+IMhUCQZYMwAwCgfIiy4iDOEC2CLDtEGQAA5UGMFRdxhmgQY/kgzAAAKD6CrByIM+SKIMsXYQYAQHERZOVDnCFzBFn+iDIAAIqJICs34gyZIMjiQZgBAFA8RFk1EGdIFVEWF8JsY30z0vbvh9dLuvyq8La9c/Ztz99irxe29XvGWmOS5reEv2dDp+xdeO9tV8xxfXdTcGixybpYzy3vNMevOWH/7M/uD6+b5Vbt+577xB5zfP7V4cds+Jy5qRY32/Pe8rQ9vrg5fN8XX2/f9/BZ++dzedS+77lrwtsvbrK37b1sP5fmtodve/y4vf7a8XfZ42M/HDDH+2bDY9MH+B2DciPIqoc4Q+IIsjgRZgAAFANRVl3EGRJBkMWNMAMAIG4EGSTiDB0ixoqBKAMAIF4EGdYjztAygqxYCDMAAOJCjKEZ4gxBxFhxEWYAAMSBIEM7iDNchSArPsIMAIB8EWToFHEGgqxECLPOeCet9obHl7YuB8eWl+zLvo8+a18m/NIbwrctSYNnO99Nz8/Z9726byU45pbs59LwC/Ztz/yyfRn/kSPhy/jPvHbe3HbowqA5Png5PDZ8wb6s+/yU/XX32CsjaPDV4a97/tnw1yxJ8zvsuXn7qaZV46kyeMn+ulyTXwNjp8LPldPvXjS37Wny/Zo5GL5tSXIL4bkPTDZ5UIAMEWRIAnFWYURZuRBmAABkjyhDq2Z3Nf87xFnFEGTlQ5QBAJAtggztaCXK6pq+H8A5t98594Bz7inn3JPOud+o/flHnXMnnXOP1P75u51PGWkZP+qv+gflQpgBQHw4diqnkTMv/wO0YnZXe2EmtXbmbFnSb3nvH3bOjUt6yDl3f23s497732/vLpEmAqw6CDMAiBbHTiVAhKFT7QZZo6Zx5r0/Lel07b+nnHNPS9rb+V0iScRYNRFma69MS/qspJ2SvKTD3vs/cs59VNI/kXS+9lc/4r2/L59ZAqgijp2KiRhDN7oJskZtfebMOXdI0hskfVfSXZI+5Jz7R5KOaO0Vopc22OYeSfdIUt+mLV1OF8RYtRFlV+GVaQDR49gpXsQYkpBUlNW1HGfOuTFJX5b0m977SefcJyT9rtZesf5dSX8g6VfXb+e9PyzpsCQN7d1PWXSAIINEmK3HK9MAYsexU3zKEmTW8hatmt5jrCGDppKOsrqW4sw516+1ncufee+/Ikne+7MN438q6c9TmWFFEWRoRJjZOnllupHvk+a3hx/jkRf7O55b74I93jNrX5dp/I0XgmMXL46Z2+77oj3v4+8O/3K/9sA5c9sfPbnHHN9+74Q5vjwc3seNPtBkXazd9s/D0MXwbQ9darKunPlMkS683l7fbfhr4bXMBpseB9lf18qQPb6wOTw2ctb+nXLlBvu2z/ykMX7e/n75Jl/30Gn7Lzhj+TfreZQ3jp3iUcQgSyK+2r0PYq25tIKsUdM4c845SZ+U9LT3/g8b/nx37ZVrSfplSU+kM8VqIMYQQpjZOn1l+qq3DU3wtiEAyeHYKV9FirEsIqxVjXMh1K6WRZTVtXLm7C5J75f0uHPukdqffUTS+5xzt2ntAOhFSR9MZYYlVrUg2/R8k1MI61y53n41tuyIsua6eWW68W1Dw7t52xCARHHslIPYoyymEGumPteqR1qWUVbXytUav6WN32vB1c/aVPYYaze+2r29KsUaYdYcr0wDiBXHTukjxNJX1TDLI8gatXW1RrSuzCGWdIThZURZW3hlGgAqIuYYK0OIVV3eQdaIOEtAWUMstgjb9PxCqc+eEWbt4ZVpACgnQix/VThrFlOQNSLO2lTGEIstwixlDDSiDABQVYQYshRrkDUizgyEGNJElAEAqiTmEJOIsbqynTUrQpA1Is4aEGPFUIazZ4RZXLyTVo0lwZbHjH1Dk93G2El7fGnEXuds9tvbg2ODt0+Z207tHzLHRzdfCY6d+I69nvfwjP0cntltDmvrD8IHQVeusw8M+mfs27bWA1sZtB/vK9favxYnXrQP3ua3hm9/bof9mO359rw5PnnQ3u/1zRqDTXY5rsnzeGVT+Oved5/9mJ68277thVfNmeP7vhD+4Tx/G4cxMYs5xgixcitakDWq9F6tbDFWxhALKWqgEWUAgLIixpC3IkdZXaXijBhDXogyAECZxBxiEjHWiSK+nbEMMbZeKeOsbBEmEWKNinTGjCgDABRZ7BEmEWLdKlqUlTHIGhU6zsoYYXXE2MaKEGYEGQCgSIoQYI2IsWQUJcrKHmPrFSLOyhxhEiHWqtjDjCgDAMSoaPHViBBLR+xhVrUgaxRVnJU9wiRCrBMxRxlBBgDIW5Hjq1FRQmz0mH3J2JkDoxnNpH2xRlmVY2y93OKsCiEmEWPdijXMiLJy8UOrWnxN+HLe/c8MB8cWtq6at33up+yDjZ0P2r8ofW/4uTb8mH2p/DNvsvezO/9T+ABicdzednqP/TOw7aklc3xxIvx17//qKXPby2/caY5bl7Nvdqn8wcv2133mXYv29j8MP1dWBu3bfuEDTfYrM/ZzTWPhx3ziiP1cGTltz23HI+G5zVxjX0q/b9oc1up8s+dxeGxxr/39KIuyxFejGEOsWXR1ejuxxFpsYUaQbSzTOOtdLHeUEWLJii3MCDIAQNZ67Nc5CiGmEEsqwDq5z7wiLZYoI8ZaE9XbGouGGCs/ggwAgPbkHWN5BFgr8oi0vMOMIGsfcdYGYix9MZwtI8haV+Yz4QCA5gix9o0em0k90PKKMmKse8SZoYgxNvD08Y62W3zN/oRn0p68o4wgaw0xBgDVRYh1L80oyyPIiLHkEWc1RQuxTiMsNnlFGTHWGmIMAKop7xCTyhFjdWlFWdZBRoylr7JxVpQYK0uErZd1lBFjrSHGAKBaYogwqVwhtl4aYZZllBFk2apEnBUhxMoaYY0IsjgRZABQDYRYtpKOsqyCjBjLVynjLPYYiy3E0v68WVZRRoy1hyh72cAlp31/Ft4drgyGD2jO3W6v8TRwvN8cv3yTPbdNz4W/T9O77V34wfvsA6DLN4UPHNyq/fzY/zV74aXJW3eY49542Fa2jZvbnv7FZXN89KnwPmf1TVfMbfvumzDHe0/Ya3LNXR/+/dNz2X4u9FxsMr5o7+N6z4bnduVWez2w3iv2c6lnMXxQOHWdualGX/OSOT7z1BZzfPOz4bFze+37rjJCLF9JhlnaUUaMxaUUcRZzjMUWYllLO8wIsvYQZABQPoRYXJIKs7SijBiLWyHjLNYYq3qINUozygiy9hFlAFAOhFjckgizNKKMICuOQsRZjDFWhhBL+u2MaQUZMdY5ogwAiq1nKY4gizHG3I9OBsf8tdm95zXGICPGiivaOCPI0lGEICPGkkGYAQA6VbQYq8sqymILMmKsPKKKM4IsHbEHGTGWLKIMANCuGGNMai3I6tIOs26DjBhDK3KPM4IsPUlGWZJBRoylgygDALQqxhhrJ8QapRllsQQZMVYducRZbEFWlhirSyrKkgoyYixdRBkAoJkyxVhdrFGWRJARY9WVaZz1zvtowowgC0siygiybBBmyfA9Tovj4V+mffOrwbHNz9i3PThpr8k1fHbevoHvPBYcmn/nneamfccvmOOjEwPBsf7L9r566hZ7HTM1eWr2zYf/wtS19kGRX7YfU/+T4bXMVh/bZG576Vb7wgs92+3H5VW7zgfHjj9z0Nx2edR+0Ba3NbsohHFAuGLvkydutNciu9i/OTi269v2bftHw9tK0tZz9hpsszvD67/tvc8+CD5mjpZfLEHWbYRtJK0w6zTKug2ymGJsYddSR9sNnrHXakRrcn9bY5bKFmR1SYRZt0FGjGWPMAMANIohxtIIsSxkHWVZxlinsZX0/RBvrSl1nJU1xuryjrI8gmzhgP0KazODx8JnCoqCKAMASMRYEjqJshiDLKsA68ZGcyTYXqlpnDnn9kv6rKSdWnujymHv/R8557ZK+oKkQ5JelPRe7739vogMlDnIyh5j3YZXu/dRxFCLLcxieZsyAMQkzWOnvIMs6xhL+u2LRY6xIgRYuxZ2LRFo67Ry5mxZ0m957x92zo1Lesg5d7+kfyzpr7z3H3POfVjShyX9dnpTba5MYRbDRT3SiLEsAqxV9bkUJdIIs6sV7YUjAJWS+LFT3lFWV4+lJCMttgt75B1jZYwwtK5pnHnvT0s6XfvvKefc05L2Snq3pLfW/tpnJD2oHOOs6GFWthiLKcKaKUKkxRRmeUdZg8K8cASgWpI8dupdXIkmzBq1GmlZLQotZfu5MUIMaWnrM2fOuUOS3iDpu5J21nY+knRGa69eowUxXFkxqRArUoQVUUxRJkUVZoV54QhAtZX92CnL+FqvSDFGhIXx1sartRxnzrkxSV+W9Jve+0nnXj64995759yGR5HOuXsk3SNJg4P25XQ7FftZszLEWBUibOHAYlRnz2IKs5iibCNlP/gBUExJHDsNDdjLP1RNVm9T7CbGCDF0o6U4c871a23n8mfe+6/U/visc2639/60c263pHMbbeu9PyzpsCRNjO9L/Ggz1jArepBVIcZiRpi1LomDn4GRLdr4b62Z3R7+xT5x1P5ZWRm2Dwr6zobX5JKkk//i7wTHdv31pLnthbvtdbV6F8Nf9OwO+9fDan+T/UuTp/CVG8JjK0P2tpsetV9E2fKu8FpjR6+z96fX7w1vK0mHxi6Z49/58q3BsR0/d8rc9qW/2GOOL+4Jr7cnSSuDPcGx4RP2q9Jzp7eZ46Oz4bEzb7HXX9v6WHhekuQv2s+lmV3h7Xc8Gu/+Kaljp02je+L5hZCDrM6OdRJjRBjS0MrVGp2kT0p62nv/hw1D90r6gKSP1f79tVRmaIgpzJKMMamzIOs0xgixl8Vw9iyWMIs9yqTkDn7Gtu6P40EHUAoxHzsVQRZBFluMbdk1lejtvXRmPNHbQ3ZaOXN2l6T3S3rcOfdI7c8+orUdyxedc78m6aik96YzxY3lHWZJx5iUXZClEWOH9tmvNrfqxRM7ErmdbuQZaIRZ6zj4ARCxKI+dYpb22xXbjbGkQizp6Or2fmONNj539rJWrtb4LUmhAnhbstNpXbM46iTe0giukCKGWFIB1uw+Ygi0PBBmbePgB0CUYj12ikGMb1NMIsTyirB2bdk1FW2gYU1bV2sskixDy5Ll58WSCLEsAqyVOeQZaHmcNYslzIqEgx8AiE+n8dUo7c+LdRNjRYmwIuGM2dVKG2dZ62Z9Man9GOs2xGKIMEtegUaYAQAQlkR8NYo1xIiwbFQxzEbO2OPEWZu6jTAp2xCLPcIsWQcaYQYAqLKkw6tRFpez7yTEiLB8EGVhxFlAEhEmZff2xCJHmCWrQCPMqm21T5rdHr5c97an54NjJ3/Gvu77NQ8tm+M/+of2Aq4DxpX2n/0N+5fb6pJ9oOJmjYOlcXvbnn77su7/82u/b47fOfp8cOxf/fDvmdtuuX3OHD83PRYc6+mz5/2Lux43x//4a+8wx5dvCT9Xpk/al6sfavIxkPEn7P1U31x4nzJ1rb2/2fSMfd9L4+HfZVsftS+VP7fD/j3olu3ncd9seO79l8OPd1msDPSmGk1piDHEsoywW3e8vGzGo+ftJTKqqEph1mqQNSLORIit9/O7frDhn99/5tWp3m9esg6z2KPsyvWDRbooCAAgR52EmJTeBTuSiLDGuMritqoScERZayoVZ0lFWF0WMdZtiIVCK63bSiPg0jx7RpgBANCasoRYkvGVhPXzKWOsVSHMugmyRqWMs6QjrC7tGOskxJKMryQ0zifJUMv7Co5JIMwAAEXRaYhJ6cRYuyEWW4C1o2xviyx7mCUVZXXRx1laoWXpJMKkbM6KxRZjRZblWbOihZn1c8dbHgGg+LqJr0axnBUrcoyVVZmjLOkga5RpnK0MuVxiK6TTCJMIsWaSPGuW9Bkzoqw7r/gZfjCXaQAAApIKr0ZpXrCDENtY0mfNslp8uqxRlmaQNYr+zFkSuokwKZuLdhQ5xNZLKsyKHGVSOcMMAJCfNKJrvXYjTEr/yolJx9hPbW5yedIG37x8U6L33aokwyyLKCPIklOaOOs2wOoIse4kEWZEGQAAa1b74zgTVhd7iLUTXu3cXlaRRpTlL48ga1S4OEsqwqRihdh7Jh5O5Ha+NHl7IrezXoxny4gytGplUJq6Prz+1ZXXhn/59G2bNm/72C57HbTBrfaBy0r/SnDsp3cfN7cd7rUPot449mJwbG//S+a2//b428zxpVX7YPKf//U/CI79b3f+ubnt507eaY4vLIV/tU2M2Wuk/fG99jpmunbWHO49OhIc2/xD+6ZXB+x9yEuvt9do88bac6PP2wdQV26073vihfD40qj9e3nokn3bffP2+Mye8Dpqz/9PE+a2OmIPV0GnESalG2LdRFjSAdbOfaYZaUmFWdpRVrYgyzvG1os6zqoUYknFV6v3k1akdarIUSYRZgCAfHUTYVJ8IZZHgLUirUhLIszSjLIyBVlsMbZeFHGWZITVxRhjWQVYK5KMtG7OmhFlAAC0p9sQk9KLsXZDLNYIayapSCPK0hd7jK2XeZzFEmJSezFW5BCzdBtpMYQZb18EAJTVan95zoolEWI/N/Jc17fx9dkbur6NJHQbZmlFWdGDrGgxtl62l9JP6Bg6ixiTOnubYlGiLG9JhBlRBgDAK7UbY1LcQVZGMYZZEaOs6CG2kSje1thMpzEmZfe5sSJGWZZnyzhLVj5pnAUHALSnkxCT0r+KYlpR9vXZG9o6e5bGWbJO38rYTZARY+UMsY1EGWdZxphUnSDL+q2LRTw7JhFk6xFhABCPWGNMyu4sWSjQ0nq7YjefK4spyIoUY1UJsY1EEWfdxJhEkFkIstYQZGsIMQCIi+/3qb9NsS6Ws2PSy58tCwVX2p8bK0OQFSXGqhxiG8klzvKIMak6nyErWpARY/kgxBr0eK2MhNeI6hkLHxj9vRufMG+62Vpj3zp3vTk+1Bfe/tTMJnPb/+uGL5jjs6vhXwH/+dKbzW3ve9V95vi7nn27Of7Omx8Pjv3evX/f3Nb32D+/K5uXg2OrJ+11sQbtZes0Ozhsjq9uDd/34oR9oOSbrDM89rz9FzYba5FdarJ7b7bG2tJoeK2x6YNN1l8btMcnnrYPRRZvDq9Nt7IYnlcVFfHsWLO3KP7cyHOZXcCj6G9XjD3GiLDWZBpnfsBndjEPqTpnx+o6ibIqBBkhtoYYA4ByKdLZsW6usph2oHUSZTEEWYwxRoB1L4q3NW6Es2OtyyLKihJjhNjVCDIAKI8ixFgSl7rPQpZBVrYYI8DSFVWcEWTtaTfKsgyyLGKMEHslYgwAyqW3fzW1y9zXdRNkWcVYUmfPsoqyJIIszxgjwPITTZxl+bbFokszyjoNsjRjjAhrLuso6/ZzowCA5MXwGbJYtRtlWQdZHiFGgMUp9zjL8mxZoy9N3l6Is2exnh1LK8YIsdblcZaMKAOA+HQSZUnZ6GxWJ2fTkvxMWewhJmUbY0RYseQWZ51GmZTcGbOYAi2rKyzGFGOEWGeIMgBAXTdhluZC0VkpwtsUswoxIqwcMo+zbqJMKv5bGTuNMKm4IUaEdY8gS1dP/6rGd4c/09HfuxIcu3nEPjD67+duNcePP7HLHL/h1hPBsXv2f8Pc9iPH3m2OP3U2fN/uIfuS8//tZvvr6js2ZI4/s3xdcKynyaXyh8/bl093PwrvwxY2m5tq8JJ937O3Lpjj4w+HL7XfP23f9ubn7ds+dZf9mE7tDf9Kn91tX85+4HKTx9TYfHDPjLltM/5nw8sPSFLfonGocma0q/susm7PlhXx7YyEGAFWBZnG2cCAvQNuJo0wS+PsWTcBVtfJJe47ibEkQ4wISx5RBgCoS+Lti0WIsiwXgI4xxAiwasv9M2etSPtsWbuBlkR8rZfFWbGkQowIywYX+QAASMl9pizGMMsyxKTOYyzpECPAEBJ1nGX5FsY0giski7Ni3YYYAZafLKOMIAOAuI30JbOfjiHMugkxKZuzYklGGAGGTkQbZ0X/bFmjtM+KdRNiRFg8iDIAQBryCrPY356YRIgRYEha0zhzzn1K0jslnfPe31L7s49K+ieS6lf3+Ij3/r60Jlk0McZYnhG2qcmH3Nt15bfhLgkAAB7ZSURBVPrBRG8vBoRZZ9g/AYhRbPumb16+KfVAi/2sWLchRoQhK62cOfu0pH8n6bPr/vzj3vvfT3xGNfXAifkMWidvT5Tai7FOz4plFWNJh1cn91n0WCPMuvJp5bB/AoAmPq3I9k3dxlOSsvisWDcxFnuIjZ0KX0G40fSe3pRngjQ0jTPv/Tecc4fSn8rG8o60TgOsruhnxfKIr3Zten6hkIHGBT+6l/f+CQA2wr5pTScRJlUjxFoNrCzug4iLSzefOfuQc+4fSToi6be89y8lNKcNWZHUbrh1G1whRQ2xIgRYM0ULNMIsdW3tn1aXezR1biw43jsaXgbk9/7KXktsaJe9BlT/Pnv8xNcPBMf+7eO/Ym/7NnvtKt8b3p8MNtk97fgf9s/b7E77Ob7tifBzsmfJ3s/N7LEPxIYvhA9IfJMfvbkd9q/FoSfD65hJ9npgsr8deuE99tc1+qK9/eDl8J1f91V7H3DyrfYaapOvM7a/0mTfazzPJKl3k70Gm7Xu3U23HzO3/ZE5motMj52y0mmESeULsSyiK2kbzZlgy0+ncfYJSb8rydf+/QeSfnWjv+icu0fSPZI0eE3na0lY0oqtkKwuY59kiJUhwJopSqARZqlraf/UuG/q3dpkZWIA6F5Hx079Oyb06Pk9iV1Ov1PdBFhdViGWVoQVMbw61fi1EmrZ6ijOvPdn6//tnPtTSX9u/N3Dkg5L0virdhXu0oBFCrEqBFgzsQcaYZa+VvdPjfumwUP7CrdvAlAsnR47jdy4x0tXx1EaoZZEfNVleeXEJEOsSvHVDkItWx3FmXNut/f+dO1/f1nSE8lNKR+dRFhduzHWTYgRYM3FGmiEWTbKuH8CUHxJ7puSDKluZRViSUUYAdad+uNHpKWnlUvpf07SWyVtd86dkPQ7kt7qnLtNa6fmX5T0wRTnmLgihBgR1p3YAo0wS0cZ908Aiq+s+6aihBgBlj4iLT2tXK3xfRv88SdTmEviuokwqXwhNvD08Y63XXzN/gRnkr4qh1mVFHn/BKC8yrBv6iTEpPZjrJsQI8JQRt1crTF33cbXell9ViyLGOsmxIouljAjygAARdBpiEnZnhkjxuLCWbN0ZBpni4t9iQdVN7K8cEeaQVblEFuPMFszeGygMm9tBAC0ppsIq6tSjI0es5c6mTkwmtFM4kOYpafQZ87a0UmI1cV0diyvEBt4+nj0b20kzNCpnnmn8R+GDzimbg7vA4bP2L+gJq6bN8dnHrzGHF81bv78bfYufP9fLpnjS2PWwlv2QdHIuWbxb/88Dp2cDo4tb7bX3Jo6YN92/2z4+zWz015sbPCyvb93TRZKW3rLZHBs+eEJc9uBi/ZzabXJr7G5HeGv7dyb7QPqni1z9rixFNnI4/bab83Wljv0DvvKg8+cCb+o++z3D9o3XjFJxFejrK+gmGWMNQuvTm+jCsFGmKWrlHGWdYjVJR1knBFrTSxRJhFmAIDkJR1dG8ljcee0YyyJAOvmPssYaoRZ+godZ91EWF0Mb1MkwjpHmAEAymxlyT7b24k8QkxKJ8byCLBWjR6bKVWgEWbZKEycJRFiUv4xVuQQi+2tjYSZjc+dAQC6CTEpnrNiMUeYpSyBRphlJ7o4SyrC6vKMsSKHWOwIMwAAXqnIb08saoA1U/+6ihpphFm2couzpCOsLq8YK3OIxXK2LKYgk4oRZZw9A4DyKvJZsbKG2HoEGdqVaZy5RRfNmTGJs2OWWIJMii/KpGKEWd1GP3MEGwAUS7chVpf3Z8aqEGUEGboR3dsaW9FNkElEWQhB1poihVlIWmeuAQDJcksutzNkEm9ZbBVBhmZmd7X29woRZ93GmMQFPUJiCjIp7iiTyhFmeCXfKy1aS1AZ3/blEXv/dOHp7eb4oe/b66BN7w2H9LYH7f2Rn7APFiZv3hIcu/ha+xd276J9sDh5rX2FufGj4bWxrlxnr3M2ctp+zBdHw/c9t8P+GZ7dbQ6rP7w8myRp+C/DlzxftpcD0+y+ZXO8Z9F+THvnwl/b+PP293PuGvsxX9ocPkDvsZfT09QhY5E0Sc8/eK057o2fsaHz7JM30k2MSclezKPMUUaQoRWtRlldlHGWRIxJXO4+hCBrH1EGAIhVtzEmEWTtKGKUEWTZazfK6qKJM4IsXbEFmVSMKJMIMwBAfAiybBFkaEWnQdYolzhLKsTqCLKw2KKsKEEmEWUAgHgQY/koWpQRZPlIIsrqMo2z3kXOkGWBIOsOUQYAiEXPEhf0yBpBhlYkGWSNonlbYysIsleKLcSk4sVYI8IMAFBknB1rT9FCTCLG8pZWlNUVIs6SjLKiizHGpGIHmUSUAQCKLckok8obZkWLMUIsDmkHWaOo4yytKCvKWbNYQ0wqfozVEWXtS/ozozHoWZKGz4bH998fvsS577cvE977wMPm+Nwv3WmOD10KH3At3LjT3PbMm+zLow9cCX8vNz9jf13zm+0DhqHz9vPk9F3h68o3u+S8nH3b1zwUfsyWx+x5rw40uUz/Ifv3Ut/3wo95s0vO3/if7ds+fZd9ULmwLTz3pSbHo0tb7cv4j/8wvHTC7C77MeszLvEvSVufavZcC2+/Yj/FK4Mw21jRQkwixmKSZZA1ijLO0jxTFmuYxRxiUnlirI4oa08ZgwwAiizpIKsrapgRYkhCXkHWKKo4q8rbF2MPMal8MSYRZO0gxgAgTmlFmVScMCtiiEnEWIxiiLH1co+zLIMsj7NmhFi+CLLWEGMAEK80g0yKP8qKGGOEWJxijLH1couzrM+SZRVmxFi+iLHWEWQAEL9mn1XsVoxhVrQYI8TiVYQYWy/bdc7mfS5vXUwjzIoQYVK5Q0wixtpBjAEA6mKJMkIMSSliiG0k97c1ZqFZSDXGW1GiK6TsMSYRZK0ixgAAITMHRjMLtKIFWB0hFq+yhNhGKhFnzRQ5yKoQYxJB1g6iDADQivXR1GmsFTW+GhFi8SpziG2EOCsYYgwbIci607voNXEsvM7TwraB4NjUfns3OrzzJ83xsePz5vjkofBCTnOvCq89JUmDL9nPi8nrwmNb/2LR3PbY2+190dKEfQGDzU+FH7fhc/a8Z3fa+4cr14Vvu/+6SXPbpefGzXGdtr/uZWPdreUmx68z++xFuxa224/Lan94fLHJ+m39l+zncd9cePttj9m3/dK7Zs3xM0Mj5vjub4bXQbt8IwfUSSpDZLWCEItP1QLMQpxFjBCDhSADAMBGiMWB+GodcRaRKsQYIdYdggwAgI0RYvkgvJJFnOWk7CFGhCWHIAMA4GqEWHqIrXwRZxkoa4gRYOmKPcryWBZjPefcpyS9U9I57/0ttT/bKukLkg5JelHSe733L+U1RwBAdwixzhFaxUOcpaCMMUaIZSP2IKuLIcxqPi3p30n6bMOffVjSX3nvP+ac+3Dt/387h7kBADpAjNkIrnJrGme8Mt0cMYZuFSXKpKjCTN77bzjnDq3743dLemvtvz8j6UERZwAyxLFT6wixjRFg1dXKmbNPi1emf4wQQ5KKFGUFstN7f7r232ck7cxzMgAq6dPi2GlDxNga4gshTeOs6q9ME2NIA1GWDe+9d85t+GA75+6RdI8kDQ5tVt9ceF2uue3h9cSmDtrfy933XzDHf/DPtpnjI6d6gmNbnrHXEjv9d+yf9d658LjvsbddbXJ8tevai+b4halrgmMDl+37XpqwH3Pnw9svPW+vY7b9Ufu2z/2EOawdjy2Ft709vF6eJE3va/KgHpo2h1cuh9dJGzpl/7qfP2Cvazf4WHj7xXH7+zX+dXvtrKmD5rAmHjsXHLtybZxHuFU/dqqreogRYOhEp585K+Ur04QYskCYpe6sc2639/60c263pA2P7Lz3hyUdlqTxiX18UwCkrZTHTnVVDDHiC2no+oIg1ivT0rpXpwc3d3t3iSHEkDWiLDP3SvqApI/V/v21fKcDAFdr69hpOJ5jp0ZViTECDFnrNM5aemVauvrV6YnxfF6dJsSQJ6IsPc65z2ntbULbnXMnJP2O1qLsi865X5N0VNJ785shAPxYR8dOY1v25/5LpOwhRoAhJp3GWZSvTBNhiAlRlj7v/fsCQ2/LdCIA0FyUx07rlTnEiDAUQSuX0o/ulekyRlgdMVYOZQ6zK9cPRnU5fQCITYzHThspa4gRYSiyVq7WmMsr02UOsPUIsvIoc5Q1uurn88HcpgEAUYrxrH5ZQ0wixlAuXV8QpB0rQ65S0RVCjJVTVcIMABC3ModYHUGGsso0zqqMICs3wqzYVvudZq8Jr2U2fD68dtXWJ+y1q47+/fB6XpK0829WzfHzbwyPX26y2JhbtZ+Xi3vCX9eZN9kvpB24314X69TBEXN85FR4nzh1g71+29iLzb7u8JgPLxsnSepdsh+znd+z9+VH3x+e+6577a/rwuvtya0s2V93z2x4+5Gz9tfVP2M/jycPhccGrpibaq7JBeMHJu3xU28PH4X32A9pqVUhwOoIMVQJcZYSYqw6CDMAQFpW+6sVYhIxhmojzhJCjFUPUQYAQPeIMeBlxFmHiLFqI8wAAGgfIQbYiLMWEWOoI8wAAGgNMQa0hzgLIMawHlEGAICNGAO6Q5w1IMgQQpgBALAxggxIDnEmogw2wqz8/NZlLf6DS8Hx8z/YFhzrta8or8Uty+b43G57/zNwKXyVNtfkMuLW5eolSSfDl0+f3W0/78+9wb70+uJp+9fL2Ez49nffeN7c9tTwVnN8y5HwsgibnjM31fKvXjTHV760wxyf+Jvh4NhLr7bve3FXeGkDSRp8IXzbkjRyOvyYzu2wnwvLo/b3e/hsePsF+9uhvll7fGGzfd9jU8ag/ZAhRQQZkI5KxhkxhlYRZgAArCHIgPRVJs4IMrSLMAMAVB1BBmSr1HFGkKFThBkAoIqIMSBfpYozYgzdIsoAAFVCjAFxKXScEWNIEmEGAKgCggyIV+HijCBDGggzAECZEWRAMRQmzogypIUwAwCUEUEGFE/UcUaQIU1EGepWV3s0PTsUHB85E94XjZ5aNW/74uvC65RJ0uAlez8387r54NjEdmsBKGn6O/aaXM5egs00c5O9wNu2vw6vNSZJ/bPhn7/zD+80t504Zz9mbjV824ub7G0nn2jymN1kDmtpW3jhrYkd0+a2fsn+lRx+JqxZMNYyO/jf7f3d+V+3FyNbntwUHOubs+c1fcD+GdnypP09WZwIjw9MsR/fCFEGFFd0cUaQIQuEGQCgTAgyoByiiDOCDFkizAAAZUCQAeWTW5wRZMgDYQYAKDKCDCi3TONsZYAoQ34IMwBAERFkQHVE8bZGIE1EGQCgqFbta+sAKBniDKVGmAEAACAGI2ea/x3iDKVElGWDtykDAADYWomyOuIMpUOYpauMQeZmejRwZCw4Prsr/Jyavs5+vvVNNlmTq8nT1S/1BMcuvLjV3na3vZDZyPHwr4ChJuuvDT1jv9fqwl3h9b4kqf+8sf2hGXPbXfeaw5rePxwcWx2w153b8qR9297eXHPvDK8XtvqtLea2y7fba43tesC+84u3hL9n0/vMTTU3PWiO+xuMde2W7efK9u/ahxoXb7fXQetZCN/+8M9eNLfVn9jDAJCmdqKsjjhDqRBm6ShjkAEAAKSlkzCTiDOUCGGWLIIMAACgPZ1GWR1xhsIjypJDkAEAALSv2yirI85QaIRZd8ocY865FyVNSVqRtOy9vyPfGQEAgLJJKsrqiDMUElHWuTIH2QZ+1nt/Ie9JAACAckk6yuq6ijNemUbWiLLOVCzIACBaHDsBxZVWkDVK4swZr0wjE4RZewgyeUl/6Zzzkv6D9/5w6C+u9ktzO8OX897+SPhOLt/Y5NrqTa6Vv+MdJ8zxqed2Bce2PRS+zL4kTV5nDmtxS3hugxft58+lW+yva+jogDk+eDk8tno+vKyBJF18rX3fM3vDcx+YNDfVlRvs8dFT9vjsqfDce6+x5731L8JLAEjS2TfZ2/fOh8cGpuxtR58YMscXJ8Lb98/Yz5WLP2NMTJK7aD9XevYayxN8Zbu5bcQ4dgIKJIsoq+NtjYgeUdY6guwqb/Hen3TOXSPpfufcD7z336gPOufukXSPJPVusdefAgAA1ZNllNXZL7s2V39l+qHagQ6QKMKsuamD7sf/4GXe+5O1f5+T9FVJd64bP+y9v8N7f0fv6GgeUwRQTRw7AZEbOZNPmEndnzkzX5mWrn51um8Tr06jNUSZjRCzOedGJfV476dq//0Lkv51ztMCAIljJyBaeQVZo67irPGVaedc/ZXpb6z7O4clHZakob37OeKGiSgLI8jaslPSV51z0tp+7r967/+ffKcEABw7ATGKIcrqOo4zXplG0gizqxFjnfPevyDp1rznAQCNOHYC4hFTkDXq5swZr0yjawTZy4gxACg9jp2AHMUaZI06jjNemUY3iLI1BBkAVAfHTkA+ihBldVxKH5kiygiyGPUuSmPHwhevHZheDo4tXrto3naz9aOufG6vPbe3zQXHJq+z18Xa+qT983b+9vBzceZG++uaeMJem2pmr33fC9vC435kxdx2x7ftX107Hgt/v1660d62b9b++Wz2dckY3vSsven0fvu+9z1gPy4LE+E196b22RdnduFl/iRJK8PhL6zZY6bpfnN4ddT+ul63+2xw7Plt19v3DaCyihRkjYgzZKbKYUaQAQAApK+oUVZHnCF1VYwyYgwAACAbRQ+yRsQZUlO1KCPIAAAAslGmIGtEnCFxVYkyYgwAACA7ZQ2yRsQZElP2KCPGAAAAslOFGFuPOENXyhxkxBgAAEB2qhhj6xFnaFtZg4wYAwAAyA4x9krEGVpWtigjxlDnnbRq7A2P/2L4uT/0vL2O2czNC+b47CV7DajBJ0fCg/bSVZq81n6Oj54ybvqovY7ZsjEtSXJNdhdbnjTWWNtrPyaLE/ZtD0yFb3tg0p6Yb/KYjp2wx2d3h9cam7x7xty278lRc/zUXeHblqTxH4W/7uHz9tc9++5Jc9w9E37Q53bZi6T1zNoP6tBF++t68sp1wbGJn7pgbqt/Yw8DyB5BZiPO0FRZoowYAwAAyBYx1h7iDEFFjjJCDAAAIHvEWHeIM7xCEaOMGAMAAMgeMZYs4gw/VqQoI8YAAACyRYiljziruKIEGTEGAACQHUIsH8RZxRBjAAAAqCPC4kKclVwRYowQAwAASBcRVgzEWUnFGGVEGGK1OiBNH1oJjm9+NLzu1tQhe42nvV+zd7O98+H7laQT758PD56211gbe9FeX6pnKbyfWNxs/7xOvNhkbasl+7775sP3PfGCvf+aarJ+266/mQ1vu3/c3HbytUvm+PgP7DXY+oylzHq+Z69jNrPXfkzHjtmP6exu4zF9w0Vz27mX7Ln1GHc9eNGe12q//f1cGm2y9pzx7Xb3bjO3BaqKECsu4qyEYggzQgwAACAbxFh5EGclkkeUEWEAAADZIcTKjTgrgSyjjBgDAADIBiFWPcRZARFjAAAA5UGEoY44K4is37JIlAEAACSPEIOFOItQXhf0IMgAAACSQYShE8RZjogwIA5uReqfDF8OvH86/LO66zv2z/HsDvsy4966Rrmk8W+HL91uzUuSloft8T7jKv0zY02+rmvseS/ZV2bXwOXw2OXX2NsOn7XHL94yFhybPmhfrn7wtP1rceqWBXN84tHB4NjMHvsxXd1uX8Z/8pomvzMmw8+Viz/aYm669TH7+3npdeHHbXXA/p2y90H7MT/+C+awxp8Pf0+m99nbAlkgwpAk4iwnaYcZAQYAAJA8YgxpIs4ylkaUEWIAAADpIMaQJeIsQ0mEGSEGAACQHmIMeSLOUpTUWTKCDGifc+7tkv5IUq+k/+i9/1jOUwIARIYQQ2yIswQl/ZZFogzojHOuV9K/l/Tzkk5I+r5z7l7v/VP5zgwAkCdiDLEjzjrAxTyA6N0p6Tnv/QuS5Jz7vKR3SyLOAKBCiDEUDXFmyONS94QZkIi9ko43/P8JSW/KaS4AgAwQYigD533nAdLuZzruuOMOf+TIkY7vD0B8nHMPee/vyHsejZxz75H0du/9r9f+//2S3uS9/1DD37lH0j2SdODAgTcePXo0l7kCSE+k+yeOnYCKs/ZN9qqT9o3WP9PxDkk3S3qfc+7mTm8PABJ0UtL+hv/fV/uzH/PeH/be3+G9v2PHjh2ZTg5ANXHsBKCZjuNMDZ/p8N4vSqp/pgMA8vZ9STc65651zg1I+hVJ9+Y8JwDg2AmAqZs42+gzHXu7mw4AdM97vyzpQ5L+h6SnJX3Re/9kvrMCAI6dANhSvyDIus91pH13ACBJ8t7fJ+m+vOcBAO3i2Amorm7OnDX9TIfE5zoAAABqOHYCYOomzvhMBwAAQOs4dgJg6vhtjd77Zedc/TMdvZI+xWc6AAAANsaxE4BmuvrMGZ/pAAAAaB3HTgAs3bytEQAAAACQEOIMAAAAACJAnAEAAABABIgzAAAAAIiA895nd2fOnZd0NLM7bG67pAt5TyIg5rlJzK9bMc+v3bkd9N4XeiEe59yUpB/mPY+AWJ8rsc5Lindusc5Lindu3c6rDPsnjp1aF/PcJObXjZjnJiV47JRpnMXGOXfEe39H3vPYSMxzk5hft2KeX8xzS0vMX3Osc4t1XlK8c4t1XlK8c4t1XlUW8/ck5rlJzK8bMc9NSnZ+vK0RAAAAACJAnAEAAABABKoeZ4fznoAh5rlJzK9bMc8v5rmlJeavOda5xTovKd65xTovKd65xTqvKov5exLz3CTm142Y5yYlOL9Kf+YMAAAAAGJR9TNnAAAAABCFSsaZc+5F59zjzrlHnHNHIpjPp5xz55xzTzT82Vbn3P3OuWdr/94S2fw+6pw7WXsMH3HO/d2c5rbfOfeAc+4p59yTzrnfqP15FI+fMb9YHr8h59z3nHOP1ub3r2p/fq1z7rvOueecc19wzg3kMb8sOOfe7pz7Ye1r/XDe86mLaT8V8z4q1v1TrPummPdJ7I/iFtM+qTafaPdLxvxy/zmrzSPK/VML84vl8Ut1X1XJtzU6516UdIf3Por1EpxzPy1pWtJnvfe31P7s/5B0yXv/sdoB4xbv/W9HNL+PSpr23v9+HnNqmNtuSbu99w8758YlPSTplyT9Y0Xw+Bnze6/iePycpFHv/bRzrl/StyT9hqR/Kekr3vvPO+f+RNKj3vtP5DnXNDjneiU9I+nnJZ2Q9H1J7/PeP5XrxBTXfirmfVSs+6dY900x75Oqvj+KXUz7JCnu/ZIxv48qjt/9Ue6fWphf7vup2vxS3VdV8sxZbLz335B0ad0fv1vSZ2r//RmtPSlzEZhfFLz3p733D9f+e0rS05L2KpLHz5hfFPya6dr/9tf+8ZLulvSl2p/n+vxL2Z2SnvPev+C9X5T0ea09d9Ag5n1UrPunWPdNMe+T2B+hHTHvl6R4901SvPunFuYXhbT3VVWNMy/pL51zDznn7sl7MgE7vfena/99RtLOPCcT8CHn3GO1U/e5vXWgzjl3SNIbJH1XET5+6+YnRfL4Oed6nXOPSDon6X5Jz0u67L1frv2VE4pop5iwvZKON/x/TF9r7Pup6H7G1oni50uKd98U4z6p4vuj2MW+T5Ii+vky5P5z1ijW/VNdjPup2rxS21dVNc7e4r2/XdI7JP2z2qnnaPm1957G9v7TT0i6XtJtkk5L+oM8J+OcG5P0ZUm/6b2fbByL4fHbYH7RPH7e+xXv/W2S9mntTNKr85oLrlKY/VQMP2PrRPPzFeu+KdZ9EvujqBVmnyRFuV+SIvk5q4t1/1QX635KSndfVck4896frP37nKSvau1Bjc3Z2ntu6++9PZfzfK7ivT9be2KuSvpT5fgY1t7v+2VJf+a9/0rtj6N5/DaaX0yPX533/rKkByS9WdJm51xfbWifpJO5TSxdJyXtb/j/aL7WAuynovkZWy+Wn69Y901F2CdVdH8UtQLsk6QIfr4sMf2cxbp/suYX0+NXl8a+qnJx5pwbrX24UM65UUm/IOkJe6tc3CvpA7X//oCkr+U4l1eo//DW/LJyegxrH8r8pKSnvfd/2DAUxeMXml9Ej98O59zm2n8Pa+3CGE9rbUfzntpfi+75l6DvS7rRrV1haUDSr2jtuZOrguynovgZ20gMP1+x7pti3iexP4pXQfZJUsT7JSmOn7PaPKLcP9XFvJ+qzSPVfVXlrtbonLtOa6/4SFKfpP/qvf+9HKck59znJL1V0nZJZyX9jqT/JumLkg5IOirpvd77XD5YGpjfW7V2WtlLelHSBxvep5zl3N4i6ZuSHpe0Wvvjj2jtvcm5P37G/N6nOB6/12vtQ6u9Wnux5ove+39d+zn5vKStkv5W0j/03i9kPb8suLVL8f6fWnsMPpX3/kCKbz8V8z4q1v1TrPummPdJ7I/iFds+SYp7v2TM762K43d/lPunFuaX+36qNr9U91WVizMAAAAAiFHl3tYIAAAAADEizgAAAAAgAsQZAAAAAESAOAMAAACACBBnAAAAABAB4gwAAAAAIkCcAQAAAEAEiDMAAAAAiMD/D5XQW+B3pjkaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "extent = (0, lg.pe.N_X, 0, lg.pe.N_Y)\n",
    "print ('params :', params, ', true : ', x_pos, y_pos, theta, sf_0)\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5))\n",
    "_ = axs[0].contourf(data.reshape(lg.pe.N_X, lg.pe.N_Y), 8, extent=extent, cmap=plt.cm.viridis, origin='upper')\n",
    "_ = axs[1].imshow(data_noisy.reshape(lg.pe.N_X, lg.pe.N_Y), cmap=plt.cm.viridis, extent=extent)\n",
    "_ = axs[2].contourf(data_fitted.reshape(lg.pe.N_X, lg.pe.N_Y), 8, extent=extent, cmap=plt.cm.viridis, origin='upper')\n",
    "for ax in axs: ax.axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With periodic boundaries, check that the filter \"re-enters\" the image from  the other border:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:46.560491Z",
     "start_time": "2018-06-20T15:16:45.867344Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "params : Parameters([('x_pos', <Parameter 'x_pos', value=14.305893862225966 +/- 0.159, bounds=[0:32]>), ('y_pos', <Parameter 'y_pos', value=7.725755287883015 +/- 0.271, bounds=[0:32]>), ('theta', <Parameter 'theta', value=0.20451580451903784 +/- 0.00531, bounds=[-1.5707963267948966:1.5707963267948966]>), ('sf_0', <Parameter 'sf_0', value=0.10113895477158497 +/- 0.00179, bounds=[0.001:inf]>), ('phase', <Parameter 'phase', value=-0.29890899924549347 +/- 0.106, bounds=[-inf:inf]>), ('B_sf', <Parameter 'B_sf', value=0.37545689608736066 +/- 0.0154, bounds=[0.001:inf]>), ('B_theta', <Parameter 'B_theta', value=0.16840243358123808 +/- 0.00736, bounds=[0.001:inf]>)]) , true :  14.6 8.5 0.20943951023931953 0.1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAEvCAYAAADB37lNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZBlZ3nn+d+be1YutauqVFVSaQNJCBBCFqaFbSzaHmjTYEcTjHE3TUfbLTqimbGj/YcJZmJMd48j6Akv7V6MpzwwQI+bZVgatUMeW/agwWCzlGTtJbShUu2rqnKv3N75I2+irFS+z93O8r7nfD8RCkn11rn3vSfvffL5nXPueZ33XgAAAACAcvWUPQEAAAAAAOEMAAAAAKJAOAMAAACACBDOAAAAACAChDMAAAAAiADhDAAAAAAi0Ffkk+3YscMfOHCgyKcEkLOHHnronPd+Z9nz6Aa1Cagm6hOAGFm1qdBwduDAAR06dKjIpwSQM+fckbLn0C1qE1BN1CcAMbJqE5c1AgAAAEAECGcAAAAAEAHCGQAAAABEgHAGAAAAABEgnAEAAABABAhnAAAAABABwhkAAAAARIBwBgAAAAARIJwBAAAAQAQIZwAAAAAQgb4ynvTOX/ndMp62Y5PXurKnACBHzz55XO98/f8cHv+fhsMbHzXGJO141Jvjp+5ZNMfHnhoIjvXYmzY1tT88t6WtC+a2bq7XHB8+YY/L2C19s/amA5fsfTq9N1yzN522t335FntcV102hzc9En4/zG+2H3t+55I5Pnzc/pXdPxUem9ljP/emE/bvuck75oJjvacGzW177F2mnkX7uZeGmvxMKu7wkdPJ9U2xoH9DGcaOdFezCg1nT5w8rdf+m9+TrnVdTzxPfJgBAADSQv+GKig0nPXOd58m88CHGQAAIC30b4hFlvmmlMsaY8AHGgAAIC30b4hJHiedahXO+EAD2YrxTDgAoDro3RCbvHufWoQzPthA9whiAIAi0LchRkX1QZUOZ3y4gc4RxgAARaFnQ4zK6IUqG874kAOtI4gBAIpEn4aYldkXVS6c8WEHmiOMXWluZ5+e/vCW4PhNu44Hx154ab/52Kd+yl67avB4eB0zSZq8MbyY2c7v2WuJvf6fP26O/39/9frg2I/d/ENz25f+8CZzfGqvOWyuyXXpFnufbXmyxxx3xuYX3rBsbju2f8IcH+iz5zbTH17nrGfJ/v3Ud8n+eY4etT+3c9vCj//Wn3jS3PZvv3abOe4vh+fWM29uqk0n7dc999OT5vjC+fA+3fxk5dqYSqNHQ+xi6I+aVjXn3JCkb0oabPz9L3vvf9M5d52kL0jaLukhSR/03jcp0dnjgw7YYig0AFAnsfdORaFHQwpi65NaOeR0WdI93vsp51y/pG855/5U0r+U9Hve+y845/5Q0i9L+mSOc/0RPuyALbZCAwA1E13vVAT6M6Qk1l6paTjz3ntJqxeg9Df+8ZLukfRLjT//rKSPK8cCwwceCIu1wABAHcXSO+WN3gypSaFfaulibedcr1ZOv98o6T9Jel7SRe/96pchjklq8g2DzvHhB66UQnEpApcOAYhV2b1TXujJkKKU+qaWwpn3fknS7c65LZK+JunmVp/AOXevpHslaWBka1uTowAAK1IqKgWr5aVDAOJXVu+UNXoxpCrV3sm+5dU63vuLkr4h6a2StjjnVsPdPkkb3s7Me3/Qe3+n9/7OvuGRlp+LYoC6Gzvif/QPNuZXhC4d+nLjzz8r6edLmB4AFNo7ZWXyWvejf4DUpN47tXK3xp2SFrz3F51zw5J+RtK/1UqheZ9WLh36kKSvZzEhCgHqKuVCUqbMLh3y4dpz7OLm4NjAhF2z9v6V/bSXN9u3dj87Gj6GNj9uP/cTfxC+Vb4kLd8TvtLz0AvXmtsO3mAf27u8zX5du//a2PakfUv5qZ+cNseHh8Ovy58ZNbcd+7/GzfH+D58yx89dE176oG/Cfl2L4/Zt+icP2L+yx37sbHDsu//v68xte/vNYbkZ41b6rzXWRZDU/5y9z4e/YoePZWO3XXyNuWlpiu6dskIPhpRVpY9q5bLGPZI+22iAeiR9yXv/J865pyR9wTn3v0r6W0mf6mYiFATUUVUKSZk6vXRo7WVDvVvLvWwIQOUU0jtlgf4LKatiH9XK3Rofk/SmDf78BUl3dTsBigLqpoqFJAaNI9RXXDrUOHu24aVD3vuDkg5K0uA1+/mhAMhM3r1TN+i7UAVV7qVauiFI1igMqJMqF5CypXrpEAAUhZ4LVVGXfqrQcLY0QJFA9dWleEQimUuHAKAI9Fmomrr1VaWcOQOqpm6FIxYxXzoEAEUhkKFq6txXEc6ADtW5cAAAykUgQ9XQV60gnAEtomgAAMpGKEPV0F9diXAGNEHRqAEv9YSXp9LM2fA6TGOX7Yd+6Wftta22P2o3WuPPhsfcsv3ePHO38aIkuR5j+4v2wleXr7NfeP/wgjk+uyO89tWmU/brunDVsDk+MTgYntdF++dx8Ub757HwnavN8dFL4e0HJuzXdfEWe+240Sa1aNLvDI4tjdnrzvW8xl6rbOzbY8GxxVPhMUla7rPnfeqt5rD6ds2En/tyfdoYQhmqhN4qrD5VDWgDRQMAEANupoYqob9qjnAGNFAw2kfDAAAALPRX7SGcobYoFp0jlAEAgBB6rM4RzlArFIvuEMoAAMBG6LFa06yXIpyh8igW2SCYAQCAteixWtdqH0U4QyVRLLJDKAMAAKvosVrXSQ9FOEMlUCjyQTADAAD0We3ppn8inCFpFIv81CqY9UjLQ+H30tZHw2tjLQ3ZD73lsL0fJ663t9/xWHheJ+6x166Sb/IznDDWMttir1N2w96z5viJv9hvP/VN4bkPn7LX+7rhtmPm+Pn/e1/4eW8wN226y/qn7L8wszf8ukaP24+tJuXs/B32z3vHofB+u7zV3qcTm+038pYz4eee3G8/9qazS+b4+fCydJKkpVObgmPjz9vPDaBz9Fid6bZ/IpwhORSLfNUqlAEAgCvQZ3Umq/6JcIYkUCiKQTADAKB+6LM6l3XvRDhDtCgUxSGUAQBQH/RY3cmzbyKcIRoUinIQzAAAqD76rGzk3TcRzlAqCkW5CGYAAFQXfVZ2iuqZCGcoHIWifIQyAACqiT4re0X2TYQzFIJCEQ+CGQAA1UKflY8yeibCGXJFsYgLwWxjfdPSju+H10u6+Nrwtr2z9mPP3WavF7bte8ZaY5LmtoZ/ZkMn7BLee/slc1zf3Rwcmm+yLtZzi7vM8auO2Z/9mf3hdbPcsv3cs5+82hyfuzm8z4bPmJtqfos9762H7fH5LeHnPv8G+7mHT9ufz8UR+7lnrwpvP7/Z3rb3ov1emt0Rfuyxo/b6a0ffY4+P/mDAHO+bCY9NXcPvGNQXPVa+yuqZCGfIHMUiTgQzAADSRo+Vv7L7JcIZMkGxiFvZhQYAAHSGHqs4MfRLhDN0hEKRhhiKDAAAaB09VvFi6pcIZ2gZxSItMRUaAAAQRo9Vntj6JcIZgigU6Yqt0AAAgFfQY5Uv1l6JcIYrUCzSF2uxAQCgzuix4hFzr0Q4A8WiQmIuNjHzTlruDY8vbFsMji0u2Ld9H3nWvk34hTeFH1uSBk93XqbnZu3nXt63FBxzC/Z7afgF+7Gnf8G+jf+mQ+Hb+E+/bs7cdujcoDk+eDE8NnzOvq373KT9unvslRE0eHP4dc89G37NkjS3056bt99qWjbeKoMX7NflmvwaGD0Rfq+cfO+8uW1Pk5/X9LXhx5Ykdzk894GJJjsFKBk9Vnxi75UIZzVGwaiW2IsNAAB1QH8VrxR6JcJZzVAwqieFQgMAQNXRY8UrpV6p6fUAzrn9zrlvOOeecs496Zz71caff9w5d9w590jjn7+X/3TRrrEj/op/UC0pFRsAqAt6p3qgx0pDar1SK2fOFiX9uvf+YefcmKSHnHMPNMZ+z3v/2/lND+2iONRHasUGAGqE3qmC6LHSk2Kv1DScee9PSjrZ+O9J59xhSXvznhhaQ6GopxSLTdacc/slfU7SLkle0kHv/e875z4u6Z9JOtv4qx/z3t9fziwB1BG9UzXQY6Ur5T6pre+cOecOSHqTpO9KulvSR5xz/1jSIa0cIXp5g23ulXSvJPVt3trldEGhqLeUi00OODINIHr0Tumgx6qG1Hullu9B65wblfQVSb/mvZ+Q9ElJN0i6XStHh35no+289we993d67+/sHRnJYMr1w/XMkNIvNlnz3p/03j/c+O9JSRyZBhAVeqf40WNVSxV6pZbOnDnn+rVSXP7Ye/9VSfLen14z/keS/iSXGdYURQJrVaHY5KmTI9Nr+T5pbkd4H296sb/jufVetsd7ZuxjZGNvPhccO39+1Nx235fseR99b3h9qeuuOWNu+8MnrzbHd9w3bo4vDodr3Mg3mqyLtcf+PAydDz/20IUm68qZ7xTp3Bvs9d2Gvx5ey2zQWEtvhf26lobs8ctbwmObTtu/Uy7daD/2qR83xs/aPy/f5HUPnbT/gjOWf7PeR2Wjd4oXPVY1VaVXauVujU7SpyQd9t7/7po/37Pmr/2CpCeyn159cMcfhFSl2OSl0yPTzrl7nXOHnHOHFmemC5svgOqjd4oLPVb1ValXauXM2d2SPijpcefcI40/+5ikDzjnbtfKF/FflPThXGZYYXUrEJufb3IKYZ1LN9hHY6uuSoUmL90cmfbeH5R0UJKG9+yv14cRQN7onUpWtx6rzqrWL7Vyt8ZvaeNrLbj7WZuqXijaDV/tPl6dwlrVCk0erCPTjTulSRyZBlACeqfiVa3HaqenqlN/tFZVe6W27taI1lWtSKyVdQjDK6paaHLCkWkAqKnU+6wse6mNHqvqga3K/RLhLAOpF4iQ2ELY5ucvV7rYVLnQ5IEj0wBQD6n2WWX2UVXtmerQKxHO2pRqgbDEFsIsVSw2dSg0AAC0KqVeK+Yeqmo9U136JcKZIaXi0KqYi0jd1KXIAAAQklKvlWIPVYWAVrd+iXC2RkoFolUpFpJmKDTImnfSsrEk2OKoURualI3R4/b4wiZ7RZOZb+8Ijg3eMWluO7l/yBwf2XIpOHbsO/Z63sPT9nt4eo85rG1Ph9dYu3S9ve5Vf5OVD6z1wJYG7f196Tr71+L4i+F5S9LctvDjz+6099nV354zxyeutete34wx2KTkuCbv46XN4de97357nx6/x37sy6+dNcf3fTH84Tx7O21MalLptarUP6XcN9WxX6p1VUulQLSqSoWkmVQLTR2LDACgvlLoterQP6XYN9W1Z6pVOEuhQLSjDsWkKupaYAAA9RN7v0X/FLe690yVDGexF4VOUEhekdKRn7oXGABAdaXQb9E/rUihd6JnWpF0OEuhKHSKYrIxigsAAMVKpd+id9pY7L0TfdOVkghnqRSFTlFMWkNxAQAgPyn1W/ROrYm5d6Jv2lhU4SylotApikn7KCwAAGQnxX6L/ql9sfZP9E620sJZioWhExST7lBYUAQ/tKz5W8K38+5/Zjg4dnnbsvnYZ37CvvX6rgft28b73vB7bfgx+1b5p95i19ld/+dIcGx+zN526mr7M7D9qQVzfH48/Lr3f+2Eue3FN+8yx63b2Te7Vf7gRft1n3rPvL39D8LvlaVB+7Ff+FCTujJtv9c0Gt7n44fs98qmk/bcdj4Sntv0Vfat9PumzGEtzzV7H4fH5vfaPw+k2W/RO3Uvxv6J3qk1hYaz3vk0i0SrKCbZiq2wUFQAADFLsceid8oe/VPaorqsMTUUlOqjoAAAYpRiEJPonfIWSzCjf+oc4awNFJT8xVBUKCitS7U5AIDUpFpvY+ydBp4+lsvjzt+8L5fHbQX9U3UQzgwxFpRmBg4f7Wi7+Vv2ZzyT9pRdVCgorUm1OQCAVKX0lZDY+qa8Qpj1fEUHNPqn6iGcNcRWUJrpNITFpqyiQjFpTSoNAQCgWLH1TUUHsY0UFcwIZOlqpa+qbTiLraiEVCWErVd0YaGQtIYwBgDYSCx9UwwhbCNFBLMyQxl9VPs67alqEc5iKSiWqoawtQhkcSKQAQDWiqVvijWIFa2MUEYP1Z4se6lKhrNYikpIbEEs7++bFVVUKCTtIZS9YuCC074/DpfDpcHwWmVn7rDXeBo42m+OX3yNPbfNz4V/TlN77BJ+7f3TTZ47vM6ZW7bfH/u/fsocn3jjTnPcG7ttafuYue3Jn1s0x0eeCtec5bdcMrftu3/cHO89Zq/JNXtD+PdPz0X7vdBzvsn4vF3jek+H53bpjfZ6YL2X7PdSz3x4XbrJ681NNXLLy+b49FNbzfEtz4bHzuy1nxvNxdIzpRzE8jprxkHteOXZQ1UinMVSWDYSWxArWt6FhULSHgIZANRbDD1TykFsvTyCGQe141RUD5VkOIuhsGyk7kFsrTwLC8WkfYQyAKinsnumKgWxInBQOy5l9E9JhLOyC8tGqhDEsr6cMa+CQiHpHKEMAOqlzJ6pLkEs67NlHNCOS9m9U7ThjECWjxQCGYUkG2UXFwBA/srulwhk7SOMxSmWvimqcFZ2gdkIgezVsi4qFJJsxVJcAAD54OxYMVIIZPRQ3YmxZyo9nBHI8pNlKMuyqFBI8hFjgQEAdK/sXqlOgUzKLpRxdVG8Yu6ZSglnZReZ9aoSxlZlFcqyKioUknzFXGAAAN3pnSu+xtctjK3KIpRxdVHcUuiZCg1nvXM+mmBGIAvLorBQTIqRQpFJge9xmh8Lr+PUN7ccHNvyjP3YgxP2mlzDp+fsB/jOY8GhuXffZW7ad/ScOT4yPhAc679o1+rJ2+x1zNTkrdlnNJyT14XXX5Mkv2jvU//j4bXMlh/bbG574Y3hNe0kqWeHvV9eu/tscOzoM9ea2y6O2Dttfrs9Nyn8HtaSXZPHb7LXIjvfvyU4tvvb9mP7R8PbStK2M/YabDO7wuu/7b3feM2SXjJHsVZdA5nUfSjj6qL4pdQvlX5ZY5GqFshWZRHMui0sFJPipVRoAABxqXMYWxVDKKN/yleMvVKzE1WVDmdVDWOryg5lZRSUy9fYR1ibGXwpfKYgFTEWGgBA/AhkK7oJZQSyNMTYK7V69WDTcOac2y/pc5J2aeVClYPe+993zm2T9EVJByS9KOn93nv7uogCVDmQVT2MdRu82n2OFINabMUmlsuUASAmsfROhLEVZZ4hI4gVL6ZeqZM+qZUzZ4uSft17/7BzbkzSQ865ByT9E0l/6b3/hHPuo5I+Kuk32p5BhqoUzGK4qUceBaWIANaq1bmkEtJiKjZS+cEsluYHADZQeu9U12BW9qWKZYWxbvqrVPqgVsTSK3XTIzUNZ977k5JONv570jl3WNJeSe+V9PbGX/uspAdVYjhLPZhVLYzFFMKaSSGkxVJspPJD2RqlNz8AsJEye6e6hbKqnxUr+qqiVTH3RBuJpU/Kokdq6ztnzrkDkt4k6buSdjWKjySd0srRa7QghjsrZlVQUgphKYql2KyKKJglc+AIQL0V2TvVIZhVLYzF2kelcOB6VQy9Upb9UcvhzDk3Kukrkn7Nez/h3CtvUO+9d85tuGecc/dKuleSBgft2+l2KvazZlUIY7EWjyxdvmY+qiIUQ7FZFVMo2wgHjgDEqMjeqarBrCphLNU+KrbeaK0Y+qQ8+qOWwplzrl8rxeWPvfdfbfzxaefcHu/9SefcHklnNtrWe39Q0kFJGh/bl/lejDWYpR7IUi0iVRFDwVmVQDDruvkZ2LRVG/+tFTM7wmspjR+xPytLw/Y6TH2nw2tySdLx//HvBMd2//WEue25e+x1tXrnwy96Zqf962G5v0l9afIWvnRjeGxpyN5286N2o7D1PeG1xo5cb9fTG/aGt5WkA6MXzPHvfOWNwbGdf/eEue3Lf3q1OT5/dXi9PUlaGuwJjg0fC68VJkmzJ7eb4yMz4bFTb7PXX9v2WHhekuTP2++l6d3h7Xc+Gm99Kqp3qmIoKyuQcSB7YzGeRYuhT8qrP2rlbo1O0qckHfbe/+6aofskfUjSJxr//nouMzTEFMyyDGNSZ4Wl06JS1WLSiRiOEMVQcKT4Q5mUXfMzum1/HDsdQCUU1TtVKZiVEcg4iN2eGHokqfw+Ke/+qJUzZ3dL+qCkx51zjzT+7GNaKSxfcs79sqQjkt6fzxQ3VnYwyzqMScUVljwKyoF99tHmVr14bGcmj9ONMotP2QVnVSLBLNoDRwBqL9feqQqhrNswJrXfN3EQu3tlB7Qy+6SieqNW7tb4LUmhd/M7sp1O65qFo07CWx6BKyTFIJZVAGv2HDEEtDIQzNoW5YEjAMi7d5q/eV8SAS2LALYqxb4J2SqrTyq6L2rrbo0pKTJoWYq87jmLglJEAGtlDmUGtDKOCMUSzFIS64EjACjC+uBTdFjLMnitlVrfZLlm37muH+OlYzsymEl26n5lUREqG86KVvSihd0WlBhCmKWsgEYwAwCkKK+wlKeUg1gWwaud54ktpBWpbn0S4axN3YYwqdggFnsIsxQd0AhmAABkr+jb2WcRxIoKX61aO5+yglrRfVIMPVIZX/UgnAVkEcKk4opKyiHMUlRAI5jV23KfNLMjfLvu7YfngmPHf8q+7/tVDy2a4z/8R3vN8QHjTvvP/qp9e/TlhQVz3M0Yt/kfs7ft6bdv6/4PX/d9c/yukeeDY//qB3/f3HbrHbPm+Jmp0eBYT58975/b/bg5/gdff5c5vnhb+L0yddy+Xf3QmDmssSfsOtU3G64pk9fZ9WbzM/ZzL4yFf5dte9S+Vf7sTvv3oFu038d9M+G5918M72+UI6UgFlsAa1UZZ9PqGMzKQjgTQWy9n9n99IZ//sCpm3N93rJQcK506YbBlG4KAgAoQdFf55A6D2KphrBmrtl3rpCAVtc+qaxeqFbhLKsQtqqIwtJtEAsFrbweK48Al+fZs7oWHAAAWpXKVzqqGsIsRQW0otAnVTScZR3CVuUdxjoJYlmGryysnU+WQa3sOzhmgYIDAIhdWVcTtRvGug1i/93uw11t364/O3VLbo+dZ0Ar8iB2TH1SmVcQRR/O8gpalqIWKaxCGEtZXQtOK6zPHZc8AkC6UryKSEovjFVBnfuktZ+TSq9ztjTkSglbIZ2GMIkg1kyWZ82yPmNW52KThVd9hh8sZRoAgDWK6K+KOnhNEGtdHmfM6JOuVHRQi/7MWRa6CWFSMTftSDmIrZdVMEs5lElpFBwAQNzKOLBd5MHrboJYKiEsr0saCWbFa/ZZzCK8VSacdRvAVhHEupNFMCOUAQCQv6IPXtchiBWBUBavlg+kPBgeSi6cZRXCpLSC2PvGH87kcb48cUcmj7NejGfLCGVo1dKgNHlDeP2rS68Lr8PUt33KfOyXdtvroA1um7Tn1r8UHPvJPUfNbYd77bXK3jz6YnBsb//L5rb//ug7zPGFZWMNNUn/w1//UnDsf7nrT8xtP3/8LnP88kL4V9v4qL1G2h/cZ69jputmzOHeI5uCY1t+YD/08oBdQ15+g71GmzfWnht53l5L7NJN9nOPvxAeXxixfy8PXbAfu2/OHp++OryO2vP//bi5rQ7Zw3WRRe9U1PfEqhbCsj5rlnUwo1eKT9ThrE5BLKvw1erz5BXSOpVyKJMoNgCAOKQSxPIMYRv1VLH1Pe0ilNVHFOEsyxC2KsYwVlQAa0WWIa2bs2aEMgAA2lfm1znKCGLd9lCtbJ91gMvqrFmWwYxQFr/Cw1ksQUxqL4ylHMQs3Ya0GIIZhQYAUGVlX0lUVBgru3da+/wxnGkjlNVTsbfSz+h9UUQYkzq7TLHswpKKLIIZhQYAUHXd9k6d9kxFnh2jd3q1rIIZvVJ6orissZlOC4tU3PfGUiwsRZ4t4yxZ9eRxFhwA0J2iwphUvUCWxdmybi5lTDWQSfRLrWqld4oynBUZxqT6BLKiL11M8eyYRIFZjxAGAPHqpmeqeyArO4xJBLIq67R/iiKcdVNYJAKZhUDWGgrMCoIYAMQtpTAmxdc3Eci6Q7/0aln3TqWEszLCmFSf75ClFsgoLuUgiK3R47W0KbxGVM9oeL2wv3/TE+ZDN1tr7FtnbjDHh/rC25+Y3mxu+x9u/KI5PrMc/hXwny+81dz2/tfeb46/59l3muPvvvXx4Nhv3fcPzG19j/35XdqyGBxbPm6vizVoL1unmcFhc3x5W/i558fttca8vTScRp+3/8IWYy2yC03Ke7M11hZGwmuNTV3bZP21QXt8/LDdiszfGl6bbmk+PK+q8AO+0EsVJQLZRghkkIrpnQoNZ90UGM6ONddJ8alDIKOorCCMAUD1cXbsFWWHsiwCGWGsXGX0TlFc1rgRzo61rohQlkoYo6BciUAGANVWxtmx2PumbkNZp4Es1bNj9E6viKFviiqcEcja027xKTKQFVFYKCavFkNRAQDkq+438rB0E8w6CWWpnR2jd7pSkX1Tq1cPRhPOirxsMXV5hrJOA1mehYVC0lzRoazb740CANpX1nfIUtBpKCsjkHE1UbliDGRrlR7OijxbttaXJ+5I4khQrGfH8iosFJLWlXGWjFAGAMXrNJRJ2QSz1V4khr6pjEsWUwhj9E+2lA5ilxbOOg1lUnZnzGIKaEXdYTGmMEYh6QyhDADqY2AgfAfQZrI+Y1Zk31TmzTxiD2P0T61JKZCtVXg46yaUSelfythNsUk1iFFEukcgy1dP/7LG9kwGx/t7l4Jjt246YT72fzvzRnP86BO7zfEb33gsOHbv/m+a237spfea40+dDj+3e8i+5fx/vdV+XX0vDZnjzyxeHxzraXKr/OGz9u3T3Q/DNezyFnNTDV6wn3vmjZfN8bGHw7fa75+yH3vL8/Zjn7jb3qeTe8O/0mf22LezH7jYZJ8amw9ePW1u24z/aTt89M0brcqpka6eu6ryvIwxq7NoWYSvtcq6qyJBLC5V6JcKDWfdHP2R8glmeRwFyqLgdHKL+07CWJZFhQKSvSoUGQBAcYr6flnW4apdZZwV4ysd8Ur1LNlGSv/OWSvyPlvWbkDLoyAVcVYsq6JCESlGlQoNACBfVb7pR1WCGP1T9qrYK0Udzoq8hLHII0BFnBXrtqhQQMoT+12EAADxqFooK+PyRL7OkRDc7UYAACAASURBVJ5Ue6XVm/scMf5OtOEs9e+WrZX3WbFuigoFJB6pFhoAQDmqEMyKPivGVURpS7VXaueOq03DmXPu05LeLemM9/62xp99XNI/k7R6d4+Pee/vb3umFRVjGCuziGxu8iX3dl26YTDTx4tBqsWmbNQnADGiNoUVGcayCGKEsDhU8fLFkFbOnH1G0n+U9Ll1f/573vvfznxGDasBJ+YzaJ1cnii1F8Y6LSxFFZOsg1cnz5l6WCOYdeUzKqE+AUATn1EBtWlt0IntLFqKlycSxOJUp2AmtRDOvPffdM4dyH8qGys7pHUawFalflasjPDVrs3PX04yoNWt2OSh7PoEABspozYVHdS6CV/rEcawkbr2Sd185+wjzrl/LOmQpF/33r+c0Zw2ZIWkdoNbt4ErJNUglkIAaya1gFbXglOgturT8mKPJs+MBsd7R8LLgPzWX9priQ3ttteA6t9njx/7i2uCY//+8V+0t32HvXaV7w3Xk8Em5Wnnn9mft5ld9nt8+xPh92TPgl3npq/uN8eHz4XXpfNNPnqzO+1fi0NPhtcxk+z1wGT/OPTC++zXNfKivf3gxfCTX/81uwYcf7u9htrE643tLzWpvcb7TJJ6N9trsFnr3r3mjpfMbX9ojpaikN4py+CUpbLunlilr3ZYUuqDmqlSn9TO982kzsPZJyX9G0m+8e/fkfRPN/qLzrl7Jd0rSYNXjXX4dLa8wlZIUbexz7KYVCGANZNKQKtSwYlUS/VpbW3q3dZkZWIA6F5UvVPeqh7EYuyrqvK1j7r3SR2FM+/96dX/ds79kaQ/Mf7uQUkHJWnstbuTO3+cUhCLsVAULfaAVveCU4RW69Pa2jR4YF9ytQlAWqraO3UTwqQ4g1hV+qnYe6L1iu6RpDj7pI7CmXNuj/f+ZON/f0HSE9lNqRydhLBV7RaWbopJVQpGnmItRgSzYlSxPgFIX+q1qdsQJsVzk7M69VKrrzXGvmgteqRXtHIr/c9LerukHc65Y5J+U9LbnXO3a+XU/IuSPpzjHDOXQhCrU+HIQ2wBjaKTjyrWJwDpS702VSGI0UddKeaQVuUeqd3vm0mt3a3xAxv88afafqYSdBPCpOoFsYHDRzvedv6W/RnOJH8xFZ8yTtPXRcr1CUB1pVKbsghhUnlBjACWvqr3SC8d21HYDUGi0G34Wq+o74oVUUy6CWKpiyWYVb3gAADillX4Wq/Ms2IEss7E0hutVZc+af3nsFlYKzSczc/3ZR6oulHkjTvyLCZ1DmLrxVJ8yi44gy8N1ObSRgCosvn5vtxCVivKXkuMMNa9WHqjtcruk8rU7POc9JmzdhR9N6C8iklZQWzg8NHoL22MpfjUueCkqmfOaewH4TWmJm8N14DhU73mY49fP2eOTz94lTm+bDz82dvtEr7/zxfM8YVRa+Gt8FphkrTpTLPwb38eh45PBccWt9hrbk1eYz92/0z45zW9y15sbPCiXe9dk4XSFt42ERxbfHjc3HbgvP1eWm7ya2x2Z/i1nXmrvYZaz9ZZe9xYimzT4/bab83WljvwrhPm+DOnwgd1n/3+tfaDoy3dBjGJMBabWHqjVfRIzVUynJV1W9asCwpnxFoTU+Gh6AAAUhBDEJMIY3mKqT+S6JFalXQ4K7OwZFlMCGGdi6nwUHQAADHKol+SOCuWkpj6IymuHin2r34kE87KLixZFZSUg1hslzbGVHhiKjqrYi8+AIBsZdUrSemdFRt4+lgujzt/875cHjdPMfVHUpw9UsyiC2dZFhap3DCWchCLXUyFh6IDAChaLP3SqlxvfJZT8Gr3uVMIajH1RxI9UidKC2dZF5VVZYWxKgexWM6WUXDax9kzAEibm3e59EyxhrEyg5gl5mAWW3+0KuY+Keb+qNBwlkeBKfsGHlUNZbEEMinOohNzwVlvo89crAUJAJCfWAOZFGcoizmQSXH2R1I6PdL6/iiW3ii6yxpbEUNxqWIoI5C1JpWiY8nrzDUAIB5ZfG9Mqlcoiz2QSfRIebF6o1aDWxb9VRLhLKYvpVYtlMUUyKS4C46UdtFBmO+V5q0lqIwf++Imuz6dO2wvNnng+/Y6aFN7w4V++4N2PfLjI+b4xK1bg2PnX2evudU7b6+bNXGdvZ7Y2JHw2liXrrfXOdt00t7n8yPh557daX+GZ/aYw+oPL88mSRr+87Hg2KK9HJhm9i2a4z3z9j7tnQ2/trHn7Z/n7FX2Pl/YEl73rsdeTk+TB4xF0iQ9/+B15rg3PmNDZ6nJa2UVxiQCWWxi74+kavdIRR7UjjKcxXikp0qhjEDWvioXHABAmrIMY1L+d1cklLWPHql+oglnBLJ8xRbIpDQKjkTRAQDEIeswJtUvkEnxh7JU+iOJHikPpYSzmI/0VCmQSfGFMgoOAACt651P53JFKc4wtopQlh16pPwUGs6yLDAEsjACWXcoOACAqqhzGJMIZFmjR8pfNJc1toJA9mqxBTEpvUKzFkUHAJC6vAOZFGcoiz2ISfRIaH6iKolwVkSRSUWMYUxKu9hIFBwAQPrq1C+lEMRW0SOhnSsHow5nua0+n8hZs1iDmJR+oVlFwWlfHl9IL1vPgjR8Ojy+/4HwLc59v32b8N5vPGyOz/78Xeb40IXwLcwv37TL3PbUW+zbow9cCv8stzxjv665Lfat2YfO2u+Tk3eH7yvf7JbzcvZjX/VQeJ8tjtrzXh5ocpv+A/bvpb7vhfd5s1vO3/Sf7cc+ebe9NMLl7eG5L9ibamGbfRv/sR+El06Y2d3kKLBxi39J2vZUs/daePsl+y1eG2WEsvmb9xVy9iylELaKHglS5/1SlOEs17U1Ig1mMQcxqTqFZhUFpz1VDGQAkLoqnSlLMYStokeClF2vFFU4q1KRscQexKTqFRqJYtMOwhgAxCuWfqnds2cpB7C16JEg5dcrlR7OiiwwZZw1I4iVi2LTGsIYAMQtlkC2XlUCl6WqfRI9UvuK6JdKC2dFF5mighlhrFwUmtYRyAAgfr1zPtpgVlX0SVhVRq9U7DpnJRWYPIJZCiFMqnaBkSgy7SCMAQDwavRKWBVDr1T6ZY1FaBak1oa3VEJXSNULjESRaVUMBQYAgJjUoU+S6JXaEVu/VItw1kzKgYwig/ViKzIAAJSFPgkbiblXIpwlhiKDjcRcZFLQO+81/lJ4nafL2weCY5P77TI6vOvHzfHRo3Pm+MSB8EJOs68Nrz0lSYMv2++LievDY9v+dN7c9qV32rVoYTy81pgkbXkqvN+Gz9jzntll14dL14cfu//6CXPbhefGzHGdtF/3orHu1mKTtcam99mLdl3eYe+X5f7w+HyT9dv6L9jv477Z8PbbH7Mf++X3zJjjp4Y2meN7/iq8DtrFm+x161A/demTJHqldqXSKxHOIlaXAkNx6UwqRQYAgDzUpU+S6JU6lWKvRDiLSB2KDMWlOykWGQAAulWHHmkVvVJ3Yu6VWrkxIuGsJFUvMhSW7MRcZAAAyFrVe6S16JeyEXuv1M7d6glnBahqkaGg5KtKhSYvzrlPS3q3pDPe+9saf7ZN0hclHZD0oqT3e+9fLmuOAABbVfukteiZshd7n7Sq3X6JcJaDKhYZikoxqlpocvQZSf9R0ufW/NlHJf2l9/4TzrmPNv7/N0qYGwBgnSr2SGvRL+Wv6r1S03DGkenmqlhoKC7FSqXQSFEFM3nvv+mcO7Duj98r6e2N//6spAdFOANQIHqnV1SxR1qLfqkYKfVJUne9Uitnzj4jjkz/SBWLDIWlPKkVm0Ts8t6fbPz3KUm7ypwMgFr6jGrYO1WxR1qPnqlYdeyTmoazuh+ZrmKhobCUr47Fpgzee++c23BnO+fulXSvJA0ObVHfbHhdrtkd4fXEJq+1f5Z7Hjhnjj/9L7ab45tO9ATHtj5jryV28u/Yn/Xe2fC477G3XW6yvNTu686b4+cmrwqODVy0n3th3N7nzoe3X3jeXsdsx6P2Y5/5MXNYOx9bCG97R3i9PEma2tdkpx6YMoeXLobXSRs6Yf+6n7vGXtdu8LHw9vNj9s9r7C/sBd4mrzWHNf7YmeDYpet22xuXpA69UxX7o/Xol8pT5z6p0++cVfLIdBULDYUlPnUuOAU57Zzb470/6ZzbI2nDzs57f1DSQUkaG9/HDwVA3pLtnarYH61HvxSPKvRJl24YzO87Z81YR6aldUenB7d0+3SZqWKhobDErQrFJhH3SfqQpE80/v31cqcDAFeKuXeqYn+0Hv1SnOiTVnQazlo6Mi1deXR6fKyco9NVLDQUlnRQbPLjnPu8Vi4T2uGcOybpN7USyr7knPtlSUckvb+8GQLAj0TVO1WxN9oI/VL86JOu1Gk4i/LIdBULDUUlXRSb/HnvPxAYekehEwGA5krpnarYG22Efik9Ve+TOr20sZVb6Ud3ZLrKhYbiUg1VLjjdXEcNAHVQdO9U5b4ohH4pXVXukdZb+9lstXdq5W6NpRyZrlOhocBUR10KzhWfzwdLmwYARCnL3mlpyNWqJwqhV6qGuvRJG2m1d+r6hiDtoMCsoMBUU50LDgAAWaNfqg56pNYVGs7qjAJTbRSdtC33O81cFV7LbPhseO2qbU/Ya1cd+Qfh9bwkadffLJvjZ98cHr/YZLExt2y/L+evDr+uU2+xD6Rd84C9LtaJazeZ45tOhGvi5I32+m2jLzZ73eExH142TpLUu2Dvs13fs2v5kQ+G5777Pvt1nXuDPbmlBft198yEt9902n5d/dP2+3jiQHhs4JK5qWab3DB+YMIeP/HO8FpmPfYuRWLolaqJHqk9hLOcUGDqg6IDAED76JWqjx6pfYSzjFBg6oeCAwBA6+iV6oMeqXOEsw5RYOqNogMAgI1eqZ7okbpDOGsRBQarKDoAAFyJPgn0R9kgnAVQZLAeRQcAgFfQK2EVPVJ2CGdrUGQQQtEBAIBeCa9Gj5QtwpkoNLBRdKrPb1vU/C9dCI6ffXp7cKzXvqO85rcumuOze+z6M3AhfPt01+Q24tbt6iVJx8O3T5/ZY7/vz7zJvvX6/En718vodPjx99x01tz2xPA2c3zrofCyCJufMzfV4j89b44vfXmnOT7+N8PBsZdvtp97fnd4aQNJGnwh/NiStOlkeJ/O7rTfC4sj9s97+HR4+8v2j0N9M/b45S32c49OGoP2LkNG6JMQQo+UvVqGM4oMWkXRAQDUDX0SmqE/yk9twhmFBu2i8AAA6oI+Ca2iP8pXpcMZhQadovAAAKqOPgntoj/KX6XCGUUG3aLoAACqij4JnaI/Kk7S4YwigyxReAAAVUKfhCzQHxUruXBGoUEeKDwAgCqgT0KW6I+Kl0w4o9ggLxQeAEDK6JGQNXqj8kQdzig2yBOFB6uWl3s0NTMUHN90KlyLRk4sm499/vXhdcokafCCXeemXz8XHBvfYS0AJU19x16Ty9lLsJmmX2Mv8Lb9r8NrjUlS/0z483f24V3mtuNn7H3mlsOPPb/Z3nbiiSb77DXmsBa2hxfeGt85ZW7rF+xfyeF3worLxlpm1/43u96d/RV7MbLFic3Bsb5Ze15T19ifka1P2j+T+fHw+MBkfes4PRLyQn9UrujCGcUGRaDwAABSQ4+EPNEbxSGKcEaxQZEoPgCAVCwN0Cchf/RG8SgtnFFoUAaKDwAAwCvojeJSaDjj6A/KRPEBAABYQV8UpyguawTyRPEBAABYQV9UnlZOUhHOUGkUIAAAgBX0RcVr96pBwhkqieJTDC5TBgAgfvRFxeqmPyKcoXIoQPmqYiBz0z0aODQaHJ/ZHX5PTV1vv9/6JpqsydXk7eoXeoJj517cZm+7x17IbNPR8K+AoSbrrw09Y69jdu7u8HpfktR/1tj+wLS57e77zGFN7R8Oji0P2OvObX3Sfmxvb67Zd4fXC1v+1lZz28U77LXGdn/DfvLzt4V/ZlP7zE01OzVojvsbjXXtFu33yo7v2q3G+TvsddB6Locff/inz5vb6g/tYaAO6IuKkVV/RDhDpVCA8lHFQAYAQJXRE+Uvj/6IcIbKoAhli0AGAEB66Ifyl2ePRDhD8ihC2SGQAQCQJvqh/BTZHxHOkDQKUXeqHMaccy9KmpS0JGnRe39nuTMCACB79EL5KKtHIpwhSRSizlU5kG3gp73358qeBAAAWaMXylYs/VFX4Ywj0ygahagzsRQcAKg7eid0i14oW7H1SFmcOePINApBMWpPbMWmBF7SnzvnvKT/3Xt/MPQXl/ul2V3h23nveCT8JBdvanJv9Sb3yt/5rmPm+ORzu4Nj2x8K32ZfkiauN4c1vzU8t8Hz9vvnwm326xo6MmCOD14Mjy2fDS9rIEnnX2c/9/Te8NwHJsxNdelGe3zkhD0+cyI8996r7Hlv+9PwEgCSdPot9va9c+GxgUl725Enhszx+fHw9v3T9nvl/E8ZE5PkztvvlZ69xvIEX91hbhsxeie0hR4oWzH3SFzWiOhRkFoXc7Epwdu898edc1dJesA597T3/purg865eyXdK0m9W+31pwAAKAt9UDZS6ZG6DWctH5kGOkFBai6VYlM07/3xxr/POOe+JukuSd9cM35Q0kFJGty/nzcagKLQO6Ep+p/spNYndRvOzCPT0pVHp/s2c3QaraEo2VIrNEVzzo1I6vHeTzb++2cl/euSpwUAEr0TDPQ/2Ui5T+oqnDU7Mt0Y+9HR6aG9HJ2GjaIUlnKhKcEuSV9zzkkrde6/eO//n3KnBAD0Tng1ep9sVKVP6jiccWQaWaM4XakqRaYM3vsXJL2x7HkAwFr0TliLvqd7VeyVujlzxpFpdI3C9IoqFhgAwBXonWqOvqd7Ve+XOg5nHJlGNyhOK6peYAAAr6B3qh/6nWzUqV/iVvooFEWqXgUmFb3z0uhL4TXDBqYWg2Pz182bj91s/ahLn99rz+0ds8GxievtdbG2PWl/3s7eEX4vTt9kv67xJ+y1qab32s99eXt43G9aMrfd+W37V9fOx8I/r5dvsrftm7E/n81el4zhzc/am07tt5973zfs/XJ5PLzm3uQ+e008F17mT5K0NBx+Yc32mab6zeHlEft1vX7P6eDY89tvsJ8bKAn9Tvfq2i8RzlCYOhequhYYAADqos59Tpbq3jMRzpC7OharuhcWAADqoo59Ttbom15BOENu6lasKCwAANRD3XqcPNA3bYxwhszVpWBRVAAAqIe69DZFoH+yEc6QmaoXLooJAAD1UfW+pkj0UK0jnKErVS5cFBIAAOqjyj1NGeijOkM4Q9uqWrwoIgAA1EdV+5my0U91h3CGllWtiFE8sMo7admohkd/LvzeH3reXsds+tbL5vjMBXsNqMEnN4UH7aWrNHGd/R4fOWE89BF7HbNFY1qS5JqUi61PGmus7bX3yfy4/dgDk+HHHpiwJ+ab7NPRY/b4zJ7wWmMT90yb2/Y9OWKOn7g7/NiSNPbD8OsePmu/7pn3Tpjj7pnwTp/dbS+S1jNj79Sh8/brevLS9cGx8Z84Z26rf2sPo36q1svEgp4qO4QzNFWVQkbhAACgfqrSx8SGviofhDMEpVzMKBgAANRXyj1M7Oix8kU4w6ukWNAoFAAA1FuK/Usq6LOKQzjDj6RU1CgSAADUW0p9S4rotcpBOKu5VAobBQIAgHpLpWdJHT1XuQhnNZNKYaMwAACAVPqW1NF3xYNwVnEpFDUKAgAASKFnqRL6rzgRzioqxgJHEUCslgekqQNLwfEtj4bX3Zo8YK/xtPfrdpntnQs/ryQd++BcePCkvcba6Iv2+lI9C+E6Mb/F/ryOv9hkbasF+7n75sLPPf6CXb8mm6zftvtvZsLb7h8zt5143YI5Pva0vQZbn7GUWc/37HXMpvfa+3T0JXufzuwx9umbzpvbzr5sz63HeOrB8/a8lvvtn+fCSJO154wft7tvu7kt4hdjv1Jl9GLxI5xVUAyFjg8/AADYSAx9Sp3Qk6WFcFYhZRQ7PvAAAKAZAlmx6M/SRTirgCILHh92AADQCgJZ8ejT0kc4SxBhDAAAxIYwVg56tWohnCWi6ILHBx0AAFgIY+WhT6suwlmEyip2fNABAEAIYaxc9Gn1QDgrESEMiINbkvonwrcD758Kf1Z3f8f+HM/stG8z7q17lEsa+3b41u3WvCRpcdge7zPu0j892uR1XWXPe8G+M7sGLobHLt5ibzt82h4/f9tocGzqWvt29YMn7V+Lk7ddNsfHHx0Mjk1fbe/T5R32bfwnrmryO2Mi/F45/8Ot5qbbHrN/nhdeH95vywP275S9D9r7/OjPmsMaez78M5naZ2+L7hHIykXPVj+Es5LkXez4MAMAgE4QyMpFD1dvhLOC5VHw+BADAIBuEMjKQx+HtQhnBcqi8PEBBgAAWSCQlYd+DiGEsxxlVfT4AAPtc869U9LvS+qV9H947z9R8pQAoFSEsXLRz6EVhLMMZV30+BADnXHO9Ur6T5J+RtIxSd93zt3nvX+q3JkBQLEIZOWil0O7CGcd4GYeQPTukvSc9/4FSXLOfUHSeyURzgBUGmGsfPRx6AbhzFBGgeMDDWRir6Sja/7/mKS3lDQXAMgVgSwO9HDIgvO+8w90u9/puPPOO/2hQ4c6fj4A8XHOPeS9v7PseazlnHufpHd673+l8f8flPQW7/1H1vydeyXdK0nXXHPNm48cOVLKXAHkJ9L6RO8E1JxVm+xVJ+0HXf1Ox7sk3SrpA865Wzt9PADI0HFJ+9f8/77Gn/2I9/6g9/5O7/2dO3fuLHRyAOqJ3glAMx2HM635Tof3fl7S6nc6AKBs35d0k3PuOufcgKRflHRfyXMCAHonAKZuwtlG3+nY2910AKB73vtFSR+R9GeSDkv6kvf+yXJnBQD0TgBsud8QZN33OvJ+OgCQJHnv75d0f9nzAIB20TsB9dXNmbOm3+mQ+F4HAABAA70TAFM34YzvdAAAALSO3gmAqePLGr33i8651e909Er6NN/pAAAA2Bi9E4BmuvrOGd/pAAAAaB29EwBLN5c1AgAAAAAyQjgDAAAAgAgQzgAAAAAgAoQzAAAAAIiA894X92TOnZV0pLAnbG6HpHNlTyIg5rlJzK9bMc+v3bld671PeiEe59ykpB+UPY+AWN8rsc5Lindusc5Lindu3c6rCvWJ3ql1Mc9NYn7diHluUoa9U6HhLDbOuUPe+zvLnsdGYp6bxPy6FfP8Yp5bXmJ+zbHOLdZ5SfHOLdZ5SfHOLdZ51VnMP5OY5yYxv27EPDcp2/lxWSMAAAAARIBwBgAAAAARqHs4O1j2BAwxz01ift2KeX4xzy0vMb/mWOcW67ykeOcW67ykeOcW67zqLOafScxzk5hfN2Kem5Th/Gr9nTMAAAAAiEXdz5wBAAAAQBRqGc6ccy865x53zj3inDsUwXw+7Zw745x7Ys2fbXPOPeCce7bx762Rze/jzrnjjX34iHPu75U0t/3OuW84555yzj3pnPvVxp9Hsf+M+cWy/4acc99zzj3amN+/avz5dc657zrnnnPOfdE5N1DG/IrgnHunc+4Hjdf60bLnsyqmOhVzjYq1PsVam2KuSdSjuMVUkxrzibYuGfMr/XPWmEeU9amF+cWy/3KtVbW8rNE596KkO733UayX4Jz7SUlTkj7nvb+t8Wf/m6QL3vtPNBrGrd7734hofh+XNOW9/+0y5rRmbnsk7fHeP+ycG5P0kKSfl/RPFMH+M+b3fsWx/5ykEe/9lHOuX9K3JP2qpH8p6ave+y845/5Q0qPe+0+WOdc8OOd6JT0j6WckHZP0fUkf8N4/VerEFFedirlGxVqfYq1NMdekutej2MVUk6S465Ixv48rjt/9UdanFuZXep1qzC/XWlXLM2ex8d5/U9KFdX/8Xkmfbfz3Z7XypixFYH5R8N6f9N4/3PjvSUmHJe1VJPvPmF8U/Iqpxv/2N/7xku6R9OXGn5f6/svZXZKe896/4L2fl/QFrbx3sEbMNSrW+hRrbYq5JlGP0I6Y65IUb22S4q1PLcwvCnnXqrqGMy/pz51zDznn7i17MgG7vPcnG/99StKuMicT8BHn3GONU/elXTqwyjl3QNKbJH1XEe6/dfOTItl/zrle59wjks5IekDS85Iueu8XG3/lmCIqihnbK+nomv+P6bXGXqei+4ytE8XnS4q3NsVYk2pej2IXe02SIvp8GUr/nK0Va31aFWOdaswrt1pV13D2Nu/9HZLeJelfNE49R8uvXHsa2/Wnn5R0g6TbJZ2U9DtlTsY5NyrpK5J+zXs/sXYshv23wfyi2X/e+yXv/e2S9mnlTNLNZc0FV0imTsXwGVsnms9XrLUp1ppEPYpaMjVJirIuSZF8zlbFWp9WxVqnpHxrVS3Dmff+eOPfZyR9TSs7NTanG9fcrl57e6bk+VzBe3+68cZclvRHKnEfNq73/YqkP/bef7Xxx9Hsv43mF9P+W+W9vyjpG5LeKmmLc66vMbRP0vHSJpav45L2r/n/aF5rAnUqms/YerF8vmKtTSnUpJrWo6glUJOkCD5flpg+Z7HWJ2t+Me2/VXnUqtqFM+fcSOPLhXLOjUj6WUlP2FuV4j5JH2r894ckfb3EubzK6oe34RdU0j5sfCnzU5IOe+9/d81QFPsvNL+I9t9O59yWxn8Pa+XGGIe1Umje1/hr0b3/MvR9STe5lTssDUj6Ra28d0qVSJ2K4jO2kRg+X7HWpphrEvUoXonUJCniuiTF8TlrzCPK+rQq5jrVmEeutap2d2t0zl2vlSM+ktQn6b9473+rxCnJOfd5SW+XtEPSaUm/Kem/SvqSpGskHZH0fu99KV8sDczv7Vo5rewlvSjpw2uuUy5ybm+T9FeSHpe03Pjjj2nl2uTS958xvw8ojv33Bq18abVXKwdrvuS9/9eNz8kXJG2T9LeS/pH3/nLR8yuCW7kV77/Tyj74dNn1QIqvTsVco2KtT7HWpphrEvUoXrHVJCnuumTM5IcVlQAAAGJJREFU7+2K43d/lPWphfmVXqca88u1VtUunAEAAABAjGp3WSMAAAAAxIhwBgAAAAARIJwBAAAAQAQIZwAAAAAQAcIZAAAAAESAcAYAAAAAESCcAQAAAEAECGcAAAAAEIH/HwXz6usIuiShAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_fitted, params = lg.LogGaborFit(data_noisy.reshape(lg.pe.N_X, lg.pe.N_Y), do_border=False)\n",
    "extent = (0, lg.pe.N_X, 0, lg.pe.N_Y)\n",
    "print ('params :', params, ', true : ', x_pos, y_pos, theta, sf_0)\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5))\n",
    "_ = axs[0].contourf(data.reshape(lg.pe.N_X, lg.pe.N_Y), 8, extent=extent, cmap=plt.cm.viridis, origin='upper')\n",
    "_ = axs[1].imshow(data_noisy.reshape(lg.pe.N_X, lg.pe.N_Y), cmap=plt.cm.viridis, extent=extent)\n",
    "_ = axs[2].contourf(data_fitted.reshape(lg.pe.N_X, lg.pe.N_Y), 8, extent=extent, cmap=plt.cm.viridis, origin='upper')\n",
    "for ax in axs: ax.axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to top](#LogGabor-user-guide)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TODO: validation of fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to top](#LogGabor-user-guide)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## more book keeping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-20T15:16:46.610433Z",
     "start_time": "2018-06-20T15:16:46.562532Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The watermark extension is already loaded. To reload it, use:\n",
      "  %reload_ext watermark\n",
      "2020-09-08T16:56:12+02:00\n",
      "\n",
      "CPython 3.8.5\n",
      "IPython 7.16.1\n",
      "\n",
      "numpy 1.20.0.dev0+7d04e22\n",
      "matplotlib 3.2.2\n",
      "scipy 1.5.2\n",
      "imageio 2.9.0\n",
      "SLIP 20191113\n",
      "LogGabor 20191113\n",
      "\n",
      "compiler   : Clang 11.0.3 (clang-1103.0.32.62)\n",
      "system     : Darwin\n",
      "release    : 19.6.0\n",
      "machine    : x86_64\n",
      "processor  : i386\n",
      "CPU cores  : 36\n",
      "interpreter: 64bit\n",
      "host name  : fortytwo\n",
      "Git hash   : 8610d956aee6574a83ae83666bd5d8ebbd88a273\n",
      "Git repo   : https://github.com/bicv/LogGabor\n",
      "Git branch : master\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -i -h -m -v -p numpy,matplotlib,scipy,imageio,SLIP,LogGabor  -r -g -b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to top](#LogGabor-user-guide)\n",
    "[Back to the LogGabor user guide](LogGabor.ipynb#LogGabor-user-guide)"
   ]
  }
 ],
 "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.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "279px",
    "width": "252px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}