BayesicFitting/examples/simplemap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Purpose: simple map example\n",
    "\n",
    "Artificial data.\n",
    "\n",
    "Demonstration: \n",
    "1. Data in the form of a map\n",
    "2. Parameter fit\n",
    "3. Results in the form of a map\n",
    "\n",
    "Author: Do Kester\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T05:00:50.086340Z",
     "iopub.status.busy": "2023-07-05T05:00:50.085608Z",
     "iopub.status.idle": "2023-07-05T05:00:50.885937Z",
     "shell.execute_reply": "2023-07-05T05:00:50.886129Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as numpy\n",
    "import math\n",
    "\n",
    "from BayesicFitting import GaussModel\n",
    "from BayesicFitting import PolySurfaceModel\n",
    "from BayesicFitting import Kernel2dModel\n",
    "from BayesicFitting import Gauss\n",
    "from BayesicFitting import AmoebaFitter\n",
    "from BayesicFitting import LevenbergMarquardtFitter\n",
    "from BayesicFitting import ImageAssistant\n",
    "from BayesicFitting import formatter as fmt\n",
    "from BayesicFitting import plotFit\n",
    "from BayesicFitting import Tools\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define some constants to make a gaussian source on a constant background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T05:00:50.888168Z",
     "iopub.status.busy": "2023-07-05T05:00:50.887872Z",
     "iopub.status.idle": "2023-07-05T05:00:50.889281Z",
     "shell.execute_reply": "2023-07-05T05:00:50.889446Z"
    }
   },
   "outputs": [],
   "source": [
    "# define some constants\n",
    "N = 71\n",
    "b0 = 3.5            # amplitude of source\n",
    "b1 = 30             # position\n",
    "b2 = 45             # position\n",
    "b3 = 13             # size of source\n",
    "b4 = 0.3            # background\n",
    "\n",
    "sig = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a map as a direct product of 2 one dimensional Gaussians. Add the background and noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T05:00:50.892258Z",
     "iopub.status.busy": "2023-07-05T05:00:50.891917Z",
     "iopub.status.idle": "2023-07-05T05:00:50.935668Z",
     "shell.execute_reply": "2023-07-05T05:00:50.935855Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD7CAYAAACscuKmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABO6klEQVR4nO29aYxk2XXfeW68eLFHZuReuVTW0lVdvW9qtdgmR6Al0RY9gjQfNIRke6AxBOiLxpAwGliiv4wHGGPoL7YEjEcAIS8ERmNJI8tjWZAlcSjStjwy2WQ3m+xmdXXtS1bumZERGeuLiDsfIjLO7wYju7J6yarquH+g0S+j3nLfFvE/5/7P/xhrrXh4eHz8EXvQA/Dw8Dge+Jfdw2NE4F92D48RgX/ZPTxGBP5l9/AYEfiX3cNjRHDPl90Yc8EY8238VzLG/LIxZtIY82VjzOXe/yeOY8AeHh7vD+Z+5tmNMYGIrIjID4nIL4rIjrX2C8aYXxORCWvtr340w/Tw8PiguN+X/a+JyP9srf2kMeaSiHzaWrtqjJkXka9Zay+81/ZhImtT6S4BiEWd/uedhBIM0xo+Hhsz+gf4iInc9W1c1zNt/TceI6i1sL5+3sG2sWj4OEzUxh+62EnG9ePOIdu2Os7fnTAYvg3PFZ/zGtgAY23qmAx2Yw3Wjw9fn2N4L57Ha+mMKdCNYo2ov9zOJvTz1iHbxg67VxyrXjNeI963wetqYzGsh2M448Y15u3i9eOlwTnwusZwbOf+8Ho32lhn+NgGj82xmjaOgfN23one6vX6rkTNysCOu4gP+/A98DMi8q96y3PW2tXe8pqIzN1r41R6Qn7g1b8rIiLJu/v9z6unx/rLyZ3G0G3baR1qO6UnnNqoO+s1plL95bDU1GPM6+dj7xT7y61Cur9cn9EHNLU+fBzh6m5/2YY6ptpjU/3l+L4+9Hz54jsVd6wLet78Auok9SmL1fVBaWf0eM0xXc7e0mvJLyObDPVYuC7p26X+cn0+p/tP48un7X5hhRUdX4DrGk3qflNXNvrLxVcW9Hibej2Cqi63cnq9wzLu1YLek8xKtb8cK+u9jmbzuu1G2RlrZwz3dBrPQ1XPIcrp9eMXO39sorxej2RRx93BF1xiQ+9pB+fTmEzqOVwvYh0dD8cmIhI09KWO4zoFe7X+cmsy219ujun9Pfjx/NZ/+d/lMBw5QWeMSYjIT4rI/z34b7ZLD4b+nBljfsEY801jzDejZmXYKh4eHseA+/ll/6yIvG6tXe/9vW6MmQeN3xi2kbX2iyLyRRGR/NiSPfjFaM7qN1RyW39F+QselvTzKK/fYvz25a+diEhY1m/EcG1Pj5HS9Tpp/OJN6rdx7uKOroNfh1hFx1G9MIvx6bHi+GWO7+uvlJDmBeCF4jKP2hyOhxAnqOAXJdTvZv6al87pr1xqW9dP3dVfvBD0r3ZSGQWpau61m/3l5gX9ZRYRMRiTBfNI3gLTyYA94Nc8LOovcm1R73v2rTU93kkwo4qyk+pSRs9nQ48bVPXa7T8x6Yw1saf3InNlW/8hrtvHi/o8kBq38vqLnHtHnweb1uekk9D9kJUGYGGJPdy3rO6TzItMo3tO+ncbz2vjdAFjVYYelnX9g9DssBBS5P6m3n5WlMKLiPyhiPxcb/nnROTf3se+PDw8jhlHetmNMVkR+YyI/AE+/oKIfMYYc1lEfqz3t4eHx0OKI9F4a21FRKYGPtsWkR+9n4MZUZrBBEmAbGUMmcfKslK+RFEpS7intLAxrfRXRCQOKtUBrWyn9XstXsUyKGP95LgeA0k20v54FevPKj0Ly/p5cwqUHOfm0HsRMU3OCujnHaPjayCJw2NH40g4vqt0nRS7Na7jsInh3+sh6GZ7aUb/YSBBJ2b4TEB0Qq9ZrDE8AWasXqf0TR1ra64wdEwE7zuTXjbQ5ez1krNNO6f/1prWEIfUv/PmRR3fy89grHre9VM6vqCm176yoPsvfHtLD4wwob6I44J6M7nXmHRfv9x1HR+vceaqhhONkzqmdkqPd3C9vy/DD3gFnYfHiMC/7B4eI4L7nWf/0MCMeoD50/3HlBamV0nXlTqJKIWNDYhwOPfYKOhy9qZmr/fPKMVKlJQmJoqadWd4YPGVyAx85q7OfzYndHxuVlWplg3dedVwW+eQmUFm5pvz3RRwtEHX6wuasU5gPjhW1+X9RZ1PT5SQEV5XCmwzeg6xuPs70MHfjSnOj0MHkNbtU1t1fI7MN0KiVh7CG4Q7UR4agts6XRvu6vUmVY+m9PxFBoU4uhxNIuT79Ev9RYZXEWZ30qt6bApp4g09h+KL0/1l57oyU87wBnqH5K6bja+f0PPg3H/pWT1GWKGIR7dNr3THymdnEP6X3cNjROBfdg+PEcHx0viO7dO1GAQipQtK3VPbSm1aWR1eYlepVmJFhRzRfME5BCl0dUEpUxMULr0O2SXEOsxAJ9eVYncyug7pc7ij67RPYP8XNUNbu6Aq4kHdPyl6uIsxFUD3A2wDGh+HUIWS3NIZ3Ta7irGCbjt1CQWd8bBI5A7qzUl1g5qGB42ZzNBtDGdYQIHrMzg3HC9eoUAJx8Y1aoKuM5Od3nRnOcRS/qr7YqjgnB9DJYqHWDeA9VMbGu5RxBLDNWqPaZjRoVBsH9R94KfWVBDa1hFebmOmZ0KvX4DnqVXoHS84/Pfb/7J7eIwI/Mvu4TEiOFYab1ptiW92RRXxbf2eiZ5ixZhSIZZiskKqcUrXp+5axKXrzFxSuBPf0Ax0O6W6aovQgiWrFIikNpS629ur/eUMS2UnNPPtlCQOCFWcMlBkbJvjPJ5SxtieHjuaU002M+oFXrOshh8hxEYsmRRkq9vQhctg9SVoKceXvaNjMqD3dWb/9xAC1PV4FEcJhsRwIravGfgE6LmN4T5v6RgGzyPCuOtTOu78dc20MwueuaYhIkU1LM1ljUIcYpsAVY1RXs+f5dUMR61xf2sPKz0OGEKg3qG+hGfgoM7iQ9LGe3h4PMLwL7uHx4jgWGm8DWLSyfXoF+gthS0xFO2bCNR2SWkRS2JZktjdSBeZWWXWPdiDAIZZ4/pwY4W4Q4FBk84u9RdJyWnQwNJXCjNERPZPQQxT1rEmdzAOzEh0UjprUZ/S87Fx/Xy4q4BbWslr7JSr3tRy0PoZFXKIuEYO2RsqUOL20Zxm9ikKid9BmemShmAxZJxpouG47WQ1+0y9ePrtlf5y+ZVlZ6wUUCXqLAPVe8QMfAJCqTbqCRLbGkLQOKO+XOgvx8soz4axBM+BGfjyKcyWrOlxRVzTEsdcA89AvKb3vTmm6zcK3Wep847Pxnt4jDz8y+7hMSI4XlGNMf2MdxPeZclt6KhRukm6mb083DWkDcGLiEhQp4cYstHICDcWlDLSLYYloe0ksuigj4kVpW10s4kXlfLVlqHvv6oUlgIWEZEMsqx0P+FXcBwioRZ029RV03iQ/nwEQ5rWBGghaChnOQYdTxhqNaHvTq0qZaa7SuKOZrVboO4ca2Nu+H1wDBaRlaYrkZkp9JeDxkB9BGZkOiGo/x3NZJcvcHuIlSp6/VjXEIemn440lWU9hzycbRxjTZzD2HV9TgaFS0S4iZkXCJc480ABWSvXHd9grQjhf9k9PEYE/mX38BgR+Jfdw2NEcMwxu8am7RRM9GEXVJ+FWydipoC+8WnGnG4RRG1e47XsDfiJs+jkthaqMJ6kb3oM68cR7wu84qmyo1KOMV35GbV6otqsuyJqpKusQ4eHPGL5BHMbVMdBQRctaL6A9fws4GllEIuWdJ9Ueg0qsRhTx1BsUTmjx0uUh0+bOo7BuzweavJnMeVV1GMxp2BxztVlfU4ydzRvICKuoy/OmwUy6TUUEuFca3PIR2zi2E4DCD1/+slHsxq/11HzP/b12/3lxuMn+suV066l2vhF2IsxR4AiIboy89lv9+6dPdyVyv+ye3iMCvzL7uExIjhWGt+Jx6Q+01Wvjb1+t/95c1nVWgf2OiLulAXroEnznCkrce2dbGL46bXDQn/ZqcEGdecUVrCP6TZML5EyBSzgAQ3nVAkbCoi41k1RDg62dHJlH7z4cEVgJ6kNdJN3tDFGfFvHuv84PANQ/93OQPkHNWEzPzClGWFaDUq7WKTn5Fh5ndHwxWmHVUNIxNrs5vD70IAqr41rwTZhg89Ac06nqjhtmoJHAWl5s6DXgE61tLuKw7IsKGrYUDuv58lpProWN84rdWe7qEF/gzbbYcEroYXpaIZjibsIieZ7Y/2gNN4YUzDG/L4x5h1jzEVjzKu+ZbOHx6OFo9L43xCRP7HWPiEiz4vIRRH5NRH5irX2vIh8pfe3h4fHQ4p70nhjzLiI/LCI/PciItbapog0jTE/JSKf7q32JRH5mojcuz97j2bYBAs5kJlfU4rYXlR6zzpo2k0Nqqeohuo08V3G4gooscJVpW0GNCq1rkoniQ3/TozoZDsNOgxbJSqaBlVpDCFYI07lFosoWFxBRRwtnVrTw+uoc+8qvWcREmluB49DbMCllJSzdaLQXw62ldIGuH5OwRCy4O0JqMF29Bo7Sjk25UAIEOL6VTnrMjDLwRCMXWB5vanGy9zUZ4B98Kh+bC4W9AB4HJJrKApC2MjZEjrz8lrQempwPaoLGR5kLqGlIsK9TC/8oLPuII7yy35GRDZF5F8YY94wxvxWrx3Ufbds9vDweHA4ysseF5GXROQ3rbUvikhFBij7kVs2N/aHreLh4XEMMNYeLpwXETHGnBCR/2KtPd37+7+S7st+TkQ+jZbNX7PWXnivfY2n5u2rp7uNX2l8n0FWcf/UcFsp9jlzmgC0B2gLTqfDQga0f6ZDZ21G6dbEm1rIwCYEwZ7SzdYUillYn47rWFnSc6D9EW2lRETKz6DBAHqasXa8hR51CTjEMtvbLCAswToJzAS0xuDqSnsnUEdafLEPnYhIB1l7McNTvsG2ikLo+st7xLp6OtDun9fZAtbLBzuwEJt2ZzP6Y0u5MwdOH0E8K3SX7ThZbT0Gmz6w2IiZfHod1E7oczL+ms4wlV+c7y9nVvT5oaXV4HVsQdTFRhnMsDM8ytzRZ+vA1+Gb3/ynUirdGXqD7vnLbq1dE5HbxpiDF/lHReR74ls2e3g8UjjqPPvfFZHfNsYkROSaiPwd6X5R/J4x5udF5KaIfO6jGaKHh8eHgaO2bP62iLw85J/uq2VzJxVI9VzXzTW91XQ+P0AW9eIO/UMm1tFa77justRM09W0g0xpvKL0LI+aYFODUGMGYo6s7pMhhHVa8eo5OP3ZKrpP2Sk6Yw1q6mybXFMK3E4qpR2/oVl0ZpOpSa9O6bFnriklJXWPIxRx2i9jpoEUm1lpEdeZNXdreIOK2lOzej4QvdDltg6hSu6G0lCKoVgv3xeLiOvkSnFJJz5Ih/Xv3G19BsqnYQOGvmx7f0WFMYUrem6sG2BtRYsCKISa0aLez8xNPTeiOQtx2K5b18HZE84iODMV7PEH8RDD3MPg5bIeHiMC/7J7eIwIjtddNmYkyna/X0wbrqaHUD6WG7YTSs3ysPahq6uI28SBWdomyj3ZTtcB6K2jW8ZVopF/O45yV1Cw8W9qVpYZZHNC6aKIG5rsPQ19O7LALfR9S1xZ7y/TUXb2Lzb1c5TgMtPOa8EW2ZKE+ANilkF7o7FruK7YbxyNESoLaPmcJ3VH0wv08qPtWACqSt06Qy4eq7Gk55/cHqDDGB+tsvLXoEuHaChzQ8dRRXjktAZHE48EQsdgV/fZnlIhDEU19Vm9LrlrGq45dQLiNi8pP6chEUMcp9S4iVmbngOtL3H18PDwL7uHx6jgWGl8LOpI9m6XstehJadTDaTaEq8qTcmssCwTLqMD/dMoRuC/URTh9EyD5rk9qTSMTi0hRBf1U0q3HWfaST0fm4VLLahxe8p1JkluaDgSNJTq0dk2vAvqCpcTp1EDqPth/drYAKN5Ii/DQOcTuu2IuHS4NouSULa5RviRfmdNj41xUw/uNJgA7SdNFlD6aAZutHvogVd16bDBzAvFM6TWhI2Rrh+uLR82ptJzKOVF/Ubqrj5j2etoBkFBU+i+fg30los14ZRERxrcBzorpW8faOMPz8r7X3YPjxGBf9k9PEYE99TGf5gYT52wr57qKmzLz4L+VEF/0HQgtqOZy91Pal+1xB404jWXtiRQskrKxP5cpKg082PpZ5RHrzen/xz01aCL4XWlrTanYod9GE4md93eXlXQYYo8KJCgcIfaaWrmKVDaO62UdOpNFeREyOqzLTZdXqIxnHPNnbGI4e/4FpotPDuLbdAiu8brBOEOaL/bCltDGjYQOZi9ERFJoQdeB/He4MxBdQ5hxl2YRvJa0swU1J2lvY1JzuDo+RzWappiIIJ6eM4W1RbcpiHtJJqRoFaig7AuAYNLis4OZiBee+P/kFJ55f1p4z08PD4e8C+7h8eI4Hh942Oxfgtemh7GouEZxPYseqZtwledJpEDjir7T6g+OUX9PahU+uJqf9kWNDNtLYQ+6xpOdA5xIKEpYP1pDTOon09toKVvyr3c+Svok4b90vgxuK3OJLEzmtVOY3Yitq2hy0RDffBJE0n7x1aKOohxDTlIYSls6f4j+slN6H4db3VcY1J3uvAwlEnswXO9oLMIzljf1POvn9Z7y7oH6sVFRFLbeFYghqmjl1/mtoYi+2dV+MQQyjEvpW88BECd0D12fz8VOgzBlx73JHPb9XfYeU6f98wKfPctZphYIsy20L2Qd/B9IPwvu4fHiMC/7B4eI4JjpfHtZCD7Z7q0OYTTCkUa+asodYSjDDOuzOIOiiByVzQDTXeV1JpSO1J3ttMhVavOKq3Oo6QzDjeX8gUV2DBkoNghQIVrYsste4wmh/ubUzMfG1MhCbXuzWndtv640tDcTQ0tGuMoHV4b3sLKGQ8cchhmibh1CsyQ52+iJVVyuL49Be06/fVbY3ADYvth0PjS85rtz1/RcKV0Ac42t1G+K24IEVTQCrqj57DxiUJ/ObkH95dVXZ9e/gTDMd63ABSaZdsWNRfZq8X+cmPeFTdNvq5trjtww+Ex6nMainDGI96bVbErh/9++192D48RgX/ZPTxGBMdK44NGW3JXuzSbWcUkyixr86DV9DSn0d4VpTsykDXe/kE1DMyCutPf22nhhO0ptkkWh1NBlpxSRELBBoUjbP9UW3JpG7XxEYQkFl1SG0sFPR47veL6BaDfDbTJimGSg6EBww9mx1vwW3eMEUUkUYKopoLMNGsZ4PHeHFcaypLTFlpEpdY0rCkhFInXday5a3rfaIJJlxteLxG3xwD919N3NAPfSmsYYCiMwawAax+o3Td4TtooHa5N63LhLX1G2farifHQEFRExDQQ7hzSBZfhTm1OxzooLhsG/8vu4TEi8C+7h8eI4HhFNUAE/3Ua+IWgKeGGUrjWpK5feZzGftDCi0hmHbQ0A00yxBUtiCuSq0rt6ktKtxwdNcpG6XiT3GHmlq2TlIZyRoBCCxGRGGgpjTMbMBJMoHtoY2J4+SrLXZsoo3VmKrCY/562NeK1HH+rqLsfqJmontJrw3ZEpLrVRR03nYWSRZqFovUUBCkBrllqEzMe5zX0SW3pvW1MgNoOtADj2BmmdUKl0KzH4L1ujumYYuhc67TearI2Y/gsShOGpclt3MMpuPnUXZclYa8ChI70u2fH2fxlnXniM30YjvSyG2NuiEhZRNoi0rLWvmyMmRSR3xWR0yJyQ0Q+Z63dPWwfHh4eDxb3Q+P/qrX2BWvtgaW07+Lq4fEI4YPQ+Pvu4tqJx6TRE7E4lBa0i6WUNce5Q2lXojS85FTEzXBTq116RjXjuUM8vWvTejnGL8EYEFr8cEM/byzC9BDH7SB7WoMIIlF0hSo0y2RZJ+kgDTg7KIHMv6H6/k5eKSPrBkj/HIHILOoBIO2unoE55kDH2RpNI3d0v1svKDWeeUPDif1lUHpk+Vso1yw+AVpdH37O1KeTunMmhNR7EOE+W4XB4POMHiO3queTvYMZkjzrFUD7kf2Pb8P8dEHPhy25TFuXWaI6CLo3sVi2ARckhogUFo1/pxua8VkdxFF/2a2I/Jkx5lvGmF/ofea7uHp4PEI46i/7p6y1K8aYWRH5sjHmHf6jtdYaY+ywDXtfDr8gIpJMFT7IWD08PD4Ajtr+aaX3/w1jzL8RkVdEZN0YM48urhuHbPtFEfmiiEh+bMkeCBgoFODXRCepn6evaVfVzpjSYcctJuNmNGP7aOGUVzLEjG0LLjRJdB7lOtE4/Np3oL2OI4OMkINZ4/G3NU9JcUlrQGsd1JVYJbc1A107gXNFTUAIN5v6OdWMU3zEMlUKREh1syt6rMxtDWnYTogZcRGRJkKF6qzua/ISOvCCusdoRAnqaid5v5jJ1k9J1w2itGYOYqMxXY4N6ElYR8FwhDMbU2/rAavzep3o3MNy150n9XmY+8/6XLINV+67GlpVntJy5NSaimfoZhMbaFtFNxzOImSu64xT7aQ+Z/kb2O98NwSzK4eHNPek8caYrDEmf7AsIn9NRN4S38XVw+ORwlF+2edE5N+YbuVOXET+L2vtnxhjXhPfxdXD45HBPV92a+01EXl+yOfbcp9dXD08PB4cjlVBZ0TVTQm0KGbc7TR9QN11hNpnTuUEA/XsCcbzUGilthgTa9wU5jU+PkyRxHpky2YJ87ofxsEsfAhhvdROufFUKzu8/pvTK7SoSl7WXm/llxf7y1SD1eGImr2N2vYCCowQK5bP6nRREq691UVXkcUxxVrY/iTiXew3QJODJpVhzItkYHXVxjKEguPX4BOQwFQY9l+dcaPReB3qR6goo3HdMfMCSbjW1g5x/A2rGDe8BBjXF39I70kzC6fYPUxb3kAeCpZgIq6ikFOlLRRQseDFKVZ6jx5vB/DaeA+PEYF/2T08RgTH27JZlBKzsIV2RuwFJoKikzvF/jIbPgwqvYiDNrYiIslN9N66AcfRedQXw7GUVHfvSaVR6Q2lfJm7OiXHaaf0BopDMK1DNZeISBzWXLRScpxJcX6lH1QHW05nNSZgoQXrptoiKDr6sNWg1HKmPcPhoYSIq3xjsdFhDSDqE3pPOa3Ggpf6FKhuGb3NUOa98SJDAP083Nf1GxMDU1h39ICVOR0H6f7EJb1HbCnNsLA+qduOXdcwrXRmeCvxwrv6PKTwTFMV2prREG9/2Q2VEvuYekPoyAYabAxB74fMre59/zAUdB4eHo84/Mvu4TEiON5svLV9hVcApVttVqlNfgX9ydBauDWFdspAayDDHYNqjhSVRSdsHECrItpVGVgpFVDn3YAbbXxFM6s50G0WwjBLP9gkIrGp6jXWi9vs8NviWEPRAmpM12dYQiUZ7aPYG660DJUhilFqU+51TZaUHm49p9dy7KZ+HqVRw17Wz8uLuq9mDqo53h/QYaoG4xAvdiC+y66DGmdcZWJ5We9d7g5mNpqk/qDuCC2ccYARVxY0nEggHGNPtuGCcZH6goZ4yU268brhR+oqirjw/NFRlspGhk3tbHd839fcA/C/7B4eIwL/snt4jAiOl8ZHbQnXujS9fkodTjMrmn5lwQtdSROoG6ajbLDh1qZbONUyG88ilM6E8sE2ssmkvbQS6sBBlDZMzbPajpkZV4Pa+2gMrqRtl+e1YXflNFXY1HOtLMGiCkUu+6f0c9JNG0DMAnrPGYLSKRT5gJKipZi0BlyO2LghszY8PIjyEB/FIGZBBGYgnglQa9OCvqRZIKXXzxMacUlpWQfb1Ajo+1CbhQfALRS/IDPP60dhULym55ld10w5rwXtyChsYd16GzMZMTgET7ztWqqx+Iq9BhO07yrpRbPTcCQ+eI7fQ1zjf9k9PEYE/mX38BgRHCuN7yTjUj/TdTMlXWJ2vDmFmugG69ZBw9kaeMrVFydvqHOqnVSaHe7rMSgE4eeOjdMYNNKbw62uQrQ+Lj+r9eWtFIQWF5Wqlc+5TSKYcc1e1sx+9TE0cQAd5HlnVjGbgWYBtUnosItwJUXb5A4S7aThHSekcYYqFj8LpOU1UFrS3ha04aSWrHMvPYbP4dg1+baus/MMnpM2nxnUe6+54dH+MjX6+vneWT3xzDrGipClg8R+DB2VqbFnuHcYdS+dBiVHnXoZVlJmQP9CCy6609KjgK7M6RV4EfQy/vQ2GIT/ZffwGBH4l93DY0RwvNr4mGYymYEnnQ13NdtYPwG9OVrpUmTArLmISHNRKTD13eGtLV3ntNJ7tlduQQxD6mQ6Hayv9D5I6nIW9k4sG21O6lgzK67VkwVtrp9UehdCGx8g+0/rK45PwCpJt3n+zTxChg3QRfSJq8yjqca2S40phumgBLUNG1SDmKANAQyz7rUZ0kwIXtBsYecZrMI+bGO4n2Vdv7goDtJaCeyUy06+o3R6+2k970RR16HFFem96ZihywxHo8zw384YxDwp1F/IQCMOutbuPaYXduIdfbb2MTuTXdNt03f2e8fy2ngPj5GHf9k9PEYEx9uyeb8hub+4IiIirSeW+59TaNDOIQt+COVxqIpxs4+Ju8X+cnRCqXG0rK2c+RVHgUR6VbPubH1cOU36DCq5ppzPNNnbTGlyYksp2KATbgWuKOlNOLBCYFM5fZjWnRpzOKhW9Nowg0y6WS/AIYfOtLjcLXYpEJe612ZQB4DsdQRxSwC3mCjPXm/DBTn1eYRNuJadFO41ZwQgvAl33d+s+jTHp8tFZOPbKVwzaFt4jTkjQd17Y4ztufXzNEpUSel5Hxh+UHA1uK/pv9zsL+8/qf34Jv7iVn+59vSCbtBrJe6z8R4eHv5l9/AYFRyvqCaTlOorXSUF9d+1JRUK0C0luYuebi2IQiiKQZmpiEh7Rqk7W+vS/D9883p/uf5Xn9BjgKplLmn2nuFAfF+5bgdGivVZaNgRfkST7Hnmlo2OXUI57yFiorCi57q/QDNJPcbeY0r1O7ijU29pGnzreRpo6jrZdbi6zIPSt1w66GTXIWhJ6Ck4Jaj1afScAx1up0HLk7qcnUA9wA7S4E2Yi+bRJ+2uxhnNGbdLxMSCDqp0WWdnmuN6g9NrCJXAhtmCuTaHUtsKahe2EMoh9AnQf441FyEciRpT4dB1RERyt3SGqvSs9iZM7qDf3XltPhGgvXntRPdZdEwoB3DkX3ZjTGCMecMY80e9v88YY75ujLlijPldY0ziXvvw8PB4cLgfGv9LInIRf/8jEfkn1tpzIrIrIj//YQ7Mw8Pjw8WRaLwxZklE/msR+Yci8j+abnuYHxGRv9lb5Usi8g9E5Dfvua+eowt94Jldp6kg6Swz9plbEIWcUrojIhLCj774ombgx9CCuf7K+f5yCq4t7Qwuh4UhJuh3rMll+ooPN4lkZn3QDLBGP/CqS0UP4FD3NdI2pbq5uzg2ssnlZZhdwtCRpajMLFP8wuy7iEgrO9yGpTGDEtkMzqGhvyPMfAc5HavdgKYfAqV4BnUMcHOJxZCxX9ZZjlToXrvdFQ27JENBlH7cnEU58xqMTXfZG44GnDIUDH1SuyhFxTPNdtcUPaW33XEHexrKJNG2OgaTyfieviCmoYM6CC+D+vDnSOTov+y/LiJ/T1TPNCUiRWvtwV25IyKLQ7bz8PB4SHCUxo4/ISIb1tpvvZ8DGGN+wRjzTWPMN6Nm5d4beHh4fCQ4Co3/pIj8pDHmb4hISkTGROQ3RKRgjIn3ft2XRGRl2MZs2TyeWbCpla6Gt1WgqFrpT/KWtjsWtPuhqwvbPQ+iuaDKDnq8t9D6h7ryoMYSUh1HbT6Dz6H5TlFTrSEA6RUz/6T6rAEQEYmgVyfdt3H6mIMCg6I183pt6ArjeI9DYNEkXWf7arRgqoPaxvfdbHwrg9AEQhUb1+XMOERJDR1fPK5jOjlV7C+vZ3UWYbmgn+/WlW/v1/W+lUv6+fws9rMN2i4iZ8+raPzGW5pqp4uPhVgnGoNXfIsuNLp+7jZaSqF8N4R7TpSFiIllxLgPiTLpvTs70wl15oDhWG1Gr2VhQ+tu9y/o+gflsYPPGHHPX3Zr7eettUvW2tMi8jMi8ufW2r8lIl8VkZ/ureZbNnt4POT4IKKaX5Vusu6KdGP4f/bhDMnDw+OjwH2Jaqy1XxORr/WWr4nIK/ezfSeMSb3XkqhOR5VdZBDRubUFAz5SdxopDrYpytxWXmUT8FOfHhB7H2zfRKdYUGOWzsa34e8OFxl2WG2n0WG0ClPKST3uYPnh2LeVbranNfxgpj0FQQVDggbDEohWYqChIbq7UqtdmwNFx+UjnR0UfHSyuEegwwEy5+02whKUgWbTei3nM1DhALUWjDnxeTLU/Z9Y3Bi67WbgZqDXSzrLMX5Ww8LijYIeo4znbxPXZhFdczeUZlOsRIPL2as6vt1zulJmY7gTjlOLse0+D3w++H4wjGwgTE1t6HU9cFaybvTl7v/wf/Lw8Pg4wb/sHh4jgmPVxotRXTupe/qqmkTSN54Uk9SdJoyttJvRbCFrT5pNoUv6murpa49Bg4xOr3QNoU995npRx4T2VAYZfrZgcow1B4QZ7OhJg8vkltKzxpSeTwuinzhaNUWZ4dwtAmUkN3bEM3gC0qt6LVmWKiLSBj8Mc6CPuxqmTC4rZY4hk52Kw3kHBQhRW4/3qdmr/eXXd0/2l+fSGpa9W1SHobmMfp7Pug5AaVB/nkVlDsah+Id6B6alreH3q4mEP69feUEvYO4u6gxOoAXVCmZ/4FpTn3B/a/dO6/WYfUOvcYCeBAw12+h/cFDLQfPWQfhfdg+PEYF/2T08RgTHazhpVBvMrGFzsdBfDneVarXQaZPa8TYa3afvuqo8U0F7HJhDplbghjMB7+3balNSO6m0mgKWqKBUOl7W/cSwDnX13DaF86HjjYhIUKN+HMaXOTjJ0JGmw0w7NNxo4dRCq6HyqeGUlKW8zD7HMJ7UthsadE7pDEMILXpiTq9/PqnUcyqln98q6QxGBLVJKq6Diqx+PpnUcOqnp1/rL38jrUbz39pVp6OXZl0915tbENLY4bMC+xW0TpqAxnxNnxmKalJbEMaUQMWnKFbS5fEbeo3q48NnSArvoqWZiFQX0JZrXa/f9st6/cavowssy1kPZrHM4el4/8vu4TEi8C+7h8eI4FhpfKxl+yWl8S3V+LamNXu983yhv5y9q9SR7i/MUMtg2ehZNedLbivlqSwrdc9dUWEH2+xQ25woKuVjqyG3DFYX2V4pLKFEE2IbeuWLdNth9ZcTNITUcUxc1G12nzykc+shzM1pcYSWVJUllHduKn2un9ZzriZcM0TZ0mt+6lkVt0yn9D5e3dOS4mxc79ezU6v95Rm0Yn0ip6KiHFLc8zByL8T0/M8k1YTxakIz89m42zvgr5xQJ6Kvb5zqL9cjCJEClMtiOZqCiCmty+VJvR6pNb1v+Rt6LWkY2Rhj3YN+zmesNum2ynXMQqf035yZF3SNDfc0/GjOdUNEe9udnSL8L7uHx4jAv+weHiMC/7J7eIwIjldBJxpfRrNQn8F2Z+Li/vetKyLSyuuUCNss25R7CmFpeH06nT+d8Rxiql9Z1JgpvQHrKkxtsXCB9cvML/AcTN1VNwVNnMcJjSfHr7HAATH/BqbJoNirzLOdte6fsSIVYJyGY9MGs49rOQHbVBExO3r9L72h016bj6sa8YUZnQLbber1e3bsbn85RGD7UuZGfzkT03P+0san+suTgT4PJ+Kaa/nhiXf7y7NxdHkQkRtNzR1cTKgbay6h57SyqxekUxnIT/RgVnUqjMq/OFIvtVnYUm0Ob54RuCmFQ9GCBVcVDUTEDllZ3MIoU+nd1Pbw51zE/7J7eIwM/Mvu4TEiOGYFnekrxagyIyVtwVWTqrLkunIntnKmm6yISFDClN5JnZ5JbaDIBXXurNuO13UcnD6Ll3RaqDqv9C8LVV6sjWkunENY1G07Oddan1OIDaisAtTYF8/qcnoL00Xg6BHaKQfONA3VXTgwVFY1TDWFcEedOufWnTcKOo7Hp3QKLAaOySmwckvP7bNjb/aXv149119+vXq6v/yp3KX+8mcm3u4vP5HQ6bk/rzzZX26jsqfaca/rlepcf3k5p8U5LVyzDuKr2/AASCT0uaw2EGqyvaDjH6DL5TO6nNfZP4nTBizJuM4ZtrOvAD4NAZ7LNrZ3Gp+U671xHsL5xf+ye3iMDPzL7uExIjj+evaeeJ+tcQVZ6sYEDfvRxjhkMQGsgH4ArZhFJFlEb68xhAFFZC5bDBt0v/nLKIpZVFVflNfl1Bay/QgBkteU2tJiqgklVCvjqpuCmo6JRRE814krep2aueE2SelNPZ/ySd3WHMLoWsjAB+gl1z6pIcfaDbf5Rn5ew6U4OO0+VGm3KnrtXyrc7i//+9Lz/eVfmXq9v/zrOy/0l38YrmHXA932bltDNhbLzIUaZryU0jbGIiKrzUJ/+Y11bWcwldZQbqOs9/SVZd3+P1/UMCNEfzfOWrAHHF1n2Y6ZVlRhDTMeuCem7d6gzDr6CMIHgY6xpPSOc3HCF8J4eHj04F92D48RwbHSeNO2/cYKcfS1qpwt9Jdz15RKN9HGuD6rvChzS9fJBO73FbORFMAw+0+KlLuhdcONOS2KSTH7j3bMAUKOAP3nLHqVNaZ1rCHq3wcLVhgGNFC7n1nVrHa9ALuqveENBhgSUfBRn8WxwuGCj1ikg4rKcMgtuKKayYzeryYURBkUvFRbsEmK6efPpJSW/58lrUn/3Lg2GbqK7DMz7eUOGkOEmlkfQ+FMteOKYp5O3+kvr4wV+suL6WJ/uYhGFN+4pSKhRA7OwNinqeg5x1H3X0cUCV2QFK7q1jtP6HVhLfz4Lfcas8iFjrLM2lO81U7pOXwoTSKMMSljzDeMMW8aY942xvwvvc99y2YPj0cIR6HxDRH5EWvt8yLygoj8uDHmE+JbNnt4PFK4J4231loROVCqhL3/rLyPls2mY/sNFNpjmn5Nr0Hwgr5lDTaS2GZmHpr3Ab15cxyUaU//rT6tVC+9Ce06M51wo+1AeJPC+KIJHXdjUjPR6TWY/9/STHFrUkODTuh+tzr186XhFleZTTRhQFhC3XujoDyPtc/s18aebs0JCIY20XCjCYER+rOJiMzCzfUTBVWMrDQKOg7Q6QVQ7huRipteQOb8dgvdFoDnE9C6h6q9p36+3NH7QL29iMg7NbWl+tuzf9lf/sK1z+r4cnqPVjb0HAxvUQy94RIMg3SZGXiivIT6d7SBDlGzXllwyXATrbTHEQa04aDMuob0GnrrTXWvxwduEmGMCYwx3xaRDRH5sohclSO2bGYX12arOmwVDw+PY8CRXnZrbdta+4J0u7W+IiJPHPUA1tovWmtftta+nIhn7r2Bh4fHR4L77fVWNMZ8VURelSO2bHa2N0rB6ZrKdsoUnlDAUp0DPQetGSxRZfkrG0iEcGANIMpp5ZRutaDFT28rbY1m9EuqndLvx+xVpapsGBGnsMGS8rldIsIq3F/RCKCd0TEVH9Nr4+izoccgtSO9r0/DgXZPj5Xc0vNszGNj9Gdr1t0M97dvL8kwnEgp5T6V3uovv1XT9c+l1vvLzLR/KqUZ9T+tas1BCC59Jq73fQYWVf+uomnwk/GiM6ZP5K70l/+4+Fx/ebsCWy9cwERKr30mpaHC3jsqLOqkoU9P4eJjkXp4R+feHC6qsQM/tWM3YYmFGRZTgXiG2+DZOrBRGxTqEEfJxs8YYwq95bSIfEZELopv2ezh8UjhKL/s8yLyJWNMIN0vh9+z1v6RMeZ7IvI7xpj/VUTeEN+y2cPjocZRsvHfEZEXh3x+3y2bbRCTqEfZSTdYPhmHiT4z7aldlCdCcJC9qJp0EZEIDScyVzSTu/+UUjJm7OkKm7mLclRQaTrSpMoQQoCuO2EJSlwpiAh33KYAlTNK/aNFpevJPaVtuVW9BtT6R2m08Z1AL7p1lOnu6zVroeEBJOYS7OERQPa5Hbg6/sKUhjXXi3otx2f0mv3Z+lP95f92QQUz65FS9JjgWhp1sHkioffxd8vn+8vPJVWQM4lW0T+IrP5eB27DInKzqdn/qVBFUy/P675u7auOv93W67ezUugvm5SO1eLYsQjXNUvaDIfgaTSVQLVwHI9AdtUN6+pTcCWqDK8XoTa+xlLvg/Vjh6fjvVzWw2NE4F92D48RwfGWuMY0m00DyHAfvcN2NRvKNsjxgUz2AaL5gvM36XTpWc3YMiOa3EI5aaBUqINsfLx4SHb8nFLvwhXdT7iHHlxwpAlgPmkiV/yRorhnTqlovUDdu65PJ53KCeqodZ0OZifoThNCYEPzyegEMsA59DwL3Kxusaw67BMTSulrbT34CxOqSf+Pu4/3lz+/8O/7y2ttLS1daRV03PjdGYvpdV0A7/0anG3WW3oSM3HXreg/7GgYsFHV+5VL6AVZL6HvHpQoEwvKuXfXIfphq2lcv+zN+NB10htYX3VVTrlrlHdDJafsGz/DQQNlrQh/A7wTB0KsDyyq8fDwePThX3YPjxHB8Za4tqwkN7t0t42MNaXNsQa81KFPp386M/kUmoi41CgJKh7uoY/ZklLJNnqgBTB0FGjmaYhJ6h5Dievus0r5cndA3dN6Djbp0rZmAXr9NQ0Dmuc0tIiBwpWXdF/T34X447SGDZ1D7iidUyjCYd+yBtoP2ym3/HKuoEaetUjHXW/r9pfLWlPLdsyXoY1nNr5udT8JPASv7atz4xroegAFy3ig9+E/FC84Y/2xqYv95f8vriW1TVycyjbUnMjG61mKhNu6fhItrOPoEt7GdU3t6rlVZ1ErsQ3qjfu584QrXJr7hp5TlMdzw9CsjPASYefBjNFh7kQi/pfdw2Nk4F92D48RwfHS+I7tU9/mhGafnXbM0KHz83hVaV741o3+cuecq9kunlOKPg5NcYsZcopktmjyB7PGBowlGXJAj9yY07GybVVtTrldO637TG0oVRcRSW4rFaeH/NhNuLAgS5+/gxAHWdcAtDyEKKl0Djr+OzBPhLkjaX+c5pMtV6iyC1eUVEIP+PbafH85ndRr+cOLV/vLf1Z8Ro8xUI56gB8ZV+pNY0m2fyK2UB5bjtyxvlbSMOC7G1ru2sZFm1ko9pc3N3RfURXUGqXAQQ2hyxQ06SW46pyEISiqV2HaI3UIyJI7LueuoeVT9rbWATDkPSgRFxFpZfSemAMnJvsBtPEeHh4fD/iX3cNjRHCsNL4TxqS+1BUzMMMdlJXOxvZBqyFmoUtLOK/ZXUeHLiJzX9voL1fPT/aX46A/8TKFO/CTD1g6Cx06MqPpVRhwYFbACQFgenlgsNndp6t4iFBeayCoYLlralv3W5vR9aszOj62I9o7ixJhlA2EFd1/+bR+Tq12NKHXIl50Zw4ySYQ1YIosCSUuFrV7KtsuzWdKQz//f4tP95c/O6Htoi43dD+vl9QYsgXnxe+sKFUXEXnppIp72OapXlOavH9bqXsMHWvttoYEmVUdX3MMoqxdhESIIMJ9OtLocn0CoVwRIeS267JkoWvfO69KHIZp8SxnCPTal852KX377Q9gOOnh4fHxgH/ZPTxGBMerjbc6+c8sdTsLUQg6U+6f1LRx4e2irg+zysFOmNFsXoYC1JPaddJvOthQxJPcUrpEqsVwIFYGH44V+ov1WeV5YcnNRNOAkmKIOrL5LYh+2MKKdarU/ccboK0F3X/1hH6eUkMZRz8fr0LElHezuts3tCQ0XsK9W9aZgyBQijqX0yx6paHXe7OmsyVhoOfz4piWn/4/2z/QX95t6rV4PKchWg4DfzulVF9E5PKO1kQ8Nqkne31XS3OLARyAtsDFQenr9MG/CzPOKlpyjcF3H11zWdOQ3tbrwvLl4jnXcDIG7/zpb+v1Y5un2qLS+zYEW5Ovb/fG5oYGzv4P/RcPD4+PFfzL7uExIjheUY21fUpCGt+cTA1dP7PO9kpwjkFmvTbnbhtH6WwcoppYU+kNDfyos2/OIgMK4U1QUcrYGgfFzisNs2gdRf/57A2IQgYED/X5nAwDM/PJXd1XeVmPx2wv6V99HsIOsNPUFsRAE8Ozya0MxjdQc8BwKZrT+zKWUxpf2tDzudSe003hgHNqVt2DOtCkX6mqrv7yns62JAO9b2+VNOvebOt55gZmBJiBvwpKX9rRexSWdPs2POGDNYi9inQi0v3ThSazigw86H0rie6uEDqxFmPyojtuzrbw2TqsdVkCszb998N3cfXw8PAvu4fHiOBYabyNGWn1RAHZ72lmNTpR6C93UAZK4Q315kkITfLfUONBEZHyD6nwog0dcr6i1DNWVmpdelEzubmrEHwg49/JKKWqTyNDu6oZeNPSsR604hER6QQ67mjMFaqQkqXXNTTZX4KPPjKuWZhPWmzLjq7ZDaV8jTwyyA0KgHQMTbSOSq8hk48OsCIiBll+M67HKK1i9iMOU0tQ9Bho/PUVpdXxhO5na18pdgSKPjemLjSlqq6zsaWimPFxt9OQBecu39H1SHCjKb2W2augzwhlOCOR2sJ12hh+LQNc49SOnlt9YvjMCYViIiL561o7u/e4nmt6c3hH1/g7+uyXP911Bmpfd58x4ii+8SeNMV81xnyv18X1l3qfTxpjvmyMudz7/8S99uXh4fHgcBQa3xKRX7HWPiUinxCRXzTGPCUivyYiX7HWnheRr/T+9vDweEhxz5fdWrtqrX29t1yWbjeYRRH5Kel2b5Xe//+bj2iMHh4eHwLuK2Y3xpyWbsOIr4vInLV2tfdPayIyd9h2/e07tl+XXj+jsRun4RJFTLfFh08jmDbM+8fc6avMHY3fWMPOgpnGlAak+e9qtYjNaqwd1HQcjWmNu7O3ULc+r3FVchfFFFDZVeECm7/uNolga2bHTRRTaVRccTpQMGsDZ6RDnWmpxMvf0X3u5OBSC0HXYNn5+GXdvphE370iplCn0JOshUYKdcSseb2uUV2PvRehCAlKvjXMBtbuQn03q9fScYEVEYEzbjit63Vuw+a1jiYbk7p+7iabb+jq9ABIaos/5xqXl/Uc8ioIlDam4bjMZiciIgGUhkkWzKzjmUZfxOi50/3lg/wR81yDOHI23hiTE5F/LSK/bK0t8d96PdyHVs07LZujyrBVPDw8jgFH7c8eSvdF/21r7R/0Pl43xsz3/n1eur3bvw9Oy+YwO2wVDw+PY8A9abwxxki3aeNFa+0/xj/9oXS7t35BjtjF1bRtv3jEqe02oM/7cIFd1mmd5I7SZCruEgO9rVpOjzbQShSdNEB1ExP6BdTBNBdrytO3lchEE5gC3IGyDpScRS2cghmsZ29MDO9xl9yF2q89fKqGraNJuTObuu3eaYwJ7I618Gko62ozw6eXRETqcJ7NXdXtmwUo+arox0dhI6bkZAVWSmiDbLM45xbCoE12WIBCbQXusGNu8Uewrc+A3YBnAFfCc0Pn2CYiguwKzg2HoCKuNoVp4Nt6I5o5/Xz8OnoIovAqse4yXSo4aX9WelwHlb2NsCRFtWRv+T2aRBwlZv+kiPx3IvJdY8y3e5/9fem+5L9njPl5EbkpIp87wr48PDweEI7SxfUv5PDvix/9cIfj4eHxUeF4baniMWlOdmlciNbHYVFpTiuP1sWoI2cxSqyuGd3WmOssGpZ0vU6ILHBdjzfx7lp/ufrsoq6DAhZS8VZBqWe4rpS+uVToLyeKetzKSaVjqU1YYLVdakzqTopfm1bqWfiuFo40ZzUbHWdLX6j6ysj+03Iqt4LiIRRcNPOg7pvDKb2ICByknJbPpsM6eWaadYMoj/MEFU+tIwP/nGacOyt6PtEE4g/Qe8vilaL7GJNy02pr7B1dj9n1LIpZ6pPD1Wq8HhHCFTbccBSLdVxL3B/S+9xAWJfYRgOStm7DwiiLMCC+D1VorziMSs5BeG28h8eIwL/sHh4jguMthAmMNMe7hwzqMOAHJQ3fvqkbnFDhTTSFrHmStcju91WIflvtAYrfH8ei2hNR5BCgh0NsD7RyWulz7aw61kY5iCg2df2xN9f7y6XnVWsUVlylCsUSHfS1i8ESq76E9DBsXemEG0fL6xyKMVoZCFswVgp1bIzXcritksgAxZ8FnY4xNEGtdRHDxn6TdylaQbb7bZ15IUVPr+j4QnRmJg0ftNByQot39bqy3x3bXFdxPnSI5TVgpj3G3oQobGnABiwJC7JwD62w63pgPj8irqgmjtC2eVKvTStrsKznFvbs3EjzB+F/2T08RgT+ZffwGBEcK40PGm3JXetms0k3oknlV/E4enOhjpxZ8HZBBRVhx6Vw1MMn15T3NU6AJsYpqNAsNUU40Rz6f0E/nr2oFD1aVErfxixCawrCEYgjUncclbFUT2k74jish5JbmpXlNSDY2rpyUq9HvI5sLJK9FnZFjlYbFJu17ePXXB7P9aJ9pZ+tDCm9rt/UU5N4ZXjGvzHJdXQ5qA+f6Y1QBkGKbQdKuBlCMFtO6s4MfBsa9QB92djsozGGLDp8BdhGO7eq9zC9oifEGgp6GNBOTMT1LqBfWLgPwRGep/A7N3TbJ3s+Dr5ls4eHh3/ZPTxGBMfe66220OViFK2kr0E4sghqC4EMmz+QXgXbrr64g2IbG0N/OIh42I65nVW6RLGODaBhRuvozhg02UBzAmWfyL4yW0txziBYFttEGBCWUBMwDrEJ2vtSkBPAXdehjOwlB/rIMkvSU1opibja+rpOkkjmLhpUQMRTnWdIoPS2MotZAZSKZtb0ANU5NO5AKa9z3zFzknCjI7F4uBJlXP8qrgcz6rgetI2KEKJQJJPc0oOzn15jUp+l+gl9DltZng/3M9Anj8/AuMYcmcvbut9TWndb/bHHddPeraOYaRD+l93DY0TgX3YPjxHBsdL4WNSR9Eo3Q97OqyrCZpT+BCgtbaeVysRLSp2iSaXS5adVICMiEiAbnUI7ZrrCxEB7nTLBlO63lWbrYz22qSr1ikOHHOU0e9+B0CcBBxvHmUZEUqCDzr8hI0whTXIb66d1rLHacD10ZU7PYept5di1OTRCgPgjgm6btF/ELYuduAThkuO8outnYDFTOgVhDEQrBhEY+6Q56+DUKotozrDGsMQZqpOp5ywCGzSEO7rjENePIq3U7vBS1s0f0GmBzCaaiYCiZy8qv4+dwrQDQFGMiDtDFZbxHkzp8UjT89fUNSkqdN+nwfoLwv+ye3iMCPzL7uExIjh2bXyr5/QSrkIkMwmHDtCQEG2QmzPopQZdffa6m9Esn89jPaWrdHax0J6zHxz3G4LGsumDIOOa2EX5KlsuTyKTiixzuOeONVZGShm2ke0UWicX9HjM2Ce3lOsnrqsjWO3J+f5y4SpMMJGZz1/U2Y/KOc3ukrrbgZ5h3FcDzTfSW3rNatPDH6fk7vAWx3TYYe86lpNm1vXzyYuk27ptM+eOlbQ+gf3aQ37aOIPBttrM/id3Ea6A6qev65RC5TzqJuYL/WWWFI9dhsB/IPqK7evzXj+j+0rc0JCgMTWPDYY4F73Hz7f/ZffwGBH4l93DY0RwrDReOiJBtUv7Wigbrc8qVaV7SxyZx8RKsb9Mf/fmpCtUKfynG/oHdMSdaaWr0TS05OgBR9EKs8CpNeicFzScYD+uydeVGicTpOE6VhO5dLP6mI4pLOk4gn2lzK2Mnl/qjlLA9jiuwRn1wY8hDGoUkB7HobefUlXM1FsaStRmh5+/iFuOmd7Sse4vwJsf0Q7FMBSqhOhgPXFZj108qxtPf2cgvX6wf5TsUiwTNN3frMyqHrx8Chrz6nAxUX0KpaKV4bUFFN4kYH66/yTqI7DPzC2l5BnM/lSW9bnPv1t0xl18Se9j9i4EZUvujFN/3LN6zQ76BXQG3G8I/8vu4TEi8C+7h8eI4Hi18Qkj1YUuLc1dVL1vBi1rqovUhSt1qp5XKsOsMU36RETqT6mBZPKuZvw78JMnmPlmRp3GlQwHDtpXiYiM0zd+QkMDZs0d3/esOwZS5U4iwOcQlTDLvwjDSZTEkq7T4DJeo5uNUuPsTV2uz3Lceqz4gFCH1zzK635zd1E3gLCrModHC+eQxD0tPqY0NL0N9xzMHHBMEWh8dh31DQOtj0ndWR9gsBrpehozG/UpuMUgpCTtpwCK2XuaxBSfLuj+NxGiQfTFtuAi7n0Pd9TFqD4PEc/1Pd2gBTecXkm20yJsAEdp2fzPjTEbxpi38Jlv1+zh8YjhKDT+X4rIjw985ts1e3g8YjhKk4j/2OveSvyUiHy6t/wlEfmaiPzqvfZlOkpj9p9UWp5e16ysI2qA9txttYSOrDNuySn926m/pwsNRTKpDaVLUV731VzWTHtiD2KbIl0p0YKppes0IYQJQLetcb9bW+hem4LpZpTX8yM1ZqdYCjvy31H3HHbHTa1q6ttihqA+h2uG5G1mDSKcgazugVGoiCtCIaWn0CeHe8cQYPVVvfZjN+Enj9Angqli9i6tY3SR9JwacRGR0jmEO7Xh2XWGB60Uy26Vcoe7eq/ZyZczOGznxPuWLILe1yDcwsxO9bTbfTa7guOd0HNIYBylJ5VEU4t/QN9t7MPPxt93u2YPD48Hiw+cjX+vds0ibsvmqOlbNnt4PCi832z8ujFm3lq7+l7tmkW6LZtF5IsiImP5RXvQtih9Q6kXNcUs+6zPgDqBDic29EuDnS9FRFrIeMfQzondYTM3iv3lxrxSKdItCjiSazrW2pJq76m3J9WP8uzOqp8zlOhuo8vhhgpmWqcL/eX0us42kFonVnT96gWIajCzQerODDJnAeIoKa4sMovtDFUSZYRUKLWtz0LUBD09w5cItQiTl+hZL0OXWaLaQJ1BCH/8KKvUu3gBTpQiktgfHkI42XyIVuhQFKvpMikzr30SPQIMWoxxRmV/SbP66bv6/FAMxedNRCS8pm3Jmi8u65hQ6k2BEo02D2ZbPooS14N2zSJHbNfs4eHxYHGUqbd/JSJ/KSIXjDF3ei2avyAinzHGXBaRH+v97eHh8RDjKNn4nz3kn95fu+Ye1aMuPHMHWWNQQdKfGgQlZgL0vurqqKlvj0EXQ994J4NchMBkUrfN3mHTe9JyeIOjvNGmYTi5AQcRlOYOlo1SqMFaAVLG/WXNnGdXdEyVx3U2g1l6i0NQO009d4DaUp7z2BUdN0MoEZFOkvpunpOuU5vW4zGLzvZPpJnMZNMxx8mUpznjgW6rEygz3XHpMEtnUwiDEqDZNIcMmnp/k+saAlDQlLmt++EMUG0W5cwU+gzXcDnZcpqJioh0Tug9DcvI+KPdGUuyKcI5EId9FNl4Dw+PRwz+ZffwGBEcr1NNzPRN9tK3NJtsQVNaaKMULyNjCgoX34Ojx7ybiSX1isaHO9Kw5RNN/lKgunVQUgp66jNsx4QWUcjAx5rDs8FtiDdE3Ex4uKIlsiZS7/wUzqedpPgD5wDqxmtmIciJwYXHMTo8hPUNZoptE5l2zDC0ICTJX9V7un9Ww6bcFZ12cGoFcK+jHGZRWKMAsVGbZot3lObSMUjEDX2iMQqckOZ3WmPpMsNFx4Mf1z4aQzswhCuNCT2H7AZCJXjIc8Yi3HLFQOUndFYqe0cz/k5H1zndF/XzpmeuGovc+0b4X3YPjxGBf9k9PEYE/mX38BgRHGvMbtq2P91gMZ3FmJPTVmxp7MRYUIaZAcFQHDFe8ZxOH7FcJrWGYhmok1jIwJg1htg6w1jqrsbZAmso1mM7LZtX3RhNEL9WntLyAtZnJ0qIj9G4oplH/HpDcxixEuK4ccSre5AqoxCIyj9OqWVvudJm1uJHmN5knTyLRdjjjkUdqdsav7fRopiFPU5fuUnkQhCOstY8yrttrZNF5EL2NN6tnhxe3FSdg9oNteeO0zHi43hFr4WB4s7GdP+xNvI2bHyCIq7qmYIz7uxtKEPpKoxpQrbzjqb0qaaq9DD4X3YPjxGBf9k9PEYEx+suC0SgmKzZjoGSZa6odVVrFtNcmE5JX9JabhGRxpmZ/vL4VaU8nHaJ7Skla0yraomFLalN3ZZUnIU2jecW+sts05y4ruNunNb6cmcKRlyVFMOGzE2106ou63mzuCSN6T3S4dbZ4fZdtXM6DoY6nPKz5/Ta8RqLiDQmhj8qsQjqLkzJBSjISVzVOqlOQafkYjgf7p0Oqaaj92TsG3f6y/ULJ/RYoMndjTA+NuKwaDSCbXK3UGCDUC6NAigWXLEfXHMc+2yyWEiX908p3c5fV7odVN15z2C9qMfeRsi2pFNydDd2il56Sk17053eJfwvu4fHiMC/7B4eI4IHRuNZBEAFkw30+6dTAO0FlSattkk3E0u6zowm1WuxnG7D5gxUNzFTHk2gBhmusyZCrTXocOlF7ceVu4aedgNqvxZUYwFq70tPFHT7m0r7YiUU54wrNTRsTY3ZBV4n3mqDTHG0qBSRDSYSKyi2F5GgrmNvoTa+PqXHy93Qexrf1VDJjul9LD5X6C+zTxybUCR39J7k4Htgs5rtT94u6nimBlSUsMRqLuqMTgYzDHQb5iwHa/XZGryVQQvvu3puiV3YUo3pc9XGPjMbKELa1OtaO68zOF3AtxUzVDx2EipPjrtvleVbNnt4ePiX3cNjRHC8opqO7dtDNZDdpCNsgLa1HVBSGwwvNGkXXHfZYEuLMVJVpTztLCg9aCgz2eylJsgIswDFwlaIFlDMmLKIgZQy8447c9Ce0sw0a6RZUEKhShg/pLAFMxtxtIVuTrHhhl6LWFHp7O4Psr8YKOKMjk3EvQYdiH7GLkEkk1Ua28nocoRrNnZFj01LK4pwGHIZ9GYuPqFjYggwKFZqTbi1+P19IaxpoI6/jXr7oIF6dt67S1v9ZZvR671/FoUzcHtlAU9tWo+1v6hNTHIrbgFPgJkDg/dARJ+tDmywOFbpHCx7Gu/hMfLwL7uHx4jgeOvZA6P9rZBtZNY8DusgZsfju6D3yKYbUGkRkeiEZl9Zz00RBXvANRZUtMIMbXMCwp2bxf4yG0/QTZQ14jxu7TT0yyWXGlNMlLmlY6LmmfX2tNxK1Vw7Lj32kAytiMQQ0uw/rQIb0mE2wGjMuuFRfWJ4PUIi54pv+mPFDEbqVrG/XDurlDTch9AJsxF8NuJFtD5G9pkzODLQ3izcRKiwADfgLB1fsc4JPdcGZhcMtDpBVc+HNe85NKgItjSkiZbRrGNb70MLrriD445mdb9N2LYlIL5izXxym9ZVvftjvC2Vh8fIw7/sHh4jguMV1XRc8Uj/Y9BhijFYApiEEKQ2p5QqLA80XkCpX+0kdNg4rqkqXQ3q0Dbv6rapBjgce7phHG1kZUmfnHJD6LFJwURE0huajWX2mqFJh+EBZghI18M1NIwA/WNGuD6t4U36LmoGduCEu1joL8cGwiNqugmKoGLIXqdguUTRS3IdMxXIjnPmgLZXFLY47blxTwZ/shpzet8ZQnCGgbMfiV3X5bU/jgpmMMpwqd3DzA5CSonrtXd69N3UWokQwqDashvWUdNPF+MIs0f5N3VGp/aYhgrBAaM/nMV/sF92Y8yPG2MuGWOuGGN8J1cPj4cY7/tlN8YEIvJPReSzIvKUiPysMeapD2tgHh4eHy4+CI1/RUSuWGuviYgYY35Huq2cv3foFoHpl4gmVjX73AEdJp1Nou8bqT5bPFcXXAFFuMdsvi5y+ybKTivzeuxx0DZT02VS4xRoKEEHkeopzfDTmTa15QoeGBKQDtMVN3tVG1G0JpClj5HSU/cOaoxyzdx3V/vL0ZKWwTaWNTuevFPUc0D/PRG3HNXJln/vpu736VM6jiTKbhmKoO8b3Wmy1zST3SroPaWwheWdCdQ01JfcmYMUWkdXLqCEGeIZ9o3jMcpn9BjJPTQcmUYtApozMDxk0ws+e+zFx6YXY2+6Iqv6GTSJQJtm1nJ0UFKbevNWf7n5zMnecT+abPyiiNzG33d6nzlgF9em7+Lq4fHA8JFn4621X7TWvmytfTmRyN57Aw8Pj48ExtrDtbTvuaExr4rIP7DW/vXe358XEbHW/m/vsc2miFREZOuwdT7GmJbRO29/zsePU9bamWH/8EFe9riIvCvdBo8rIvKaiPxNa+3b99jum9bal9/XQR9hjOJ5+3N+uPC+E3TW2pYx5n8QkT8VkUBE/vm9XnQPD48Hhw8kqrHW/rGI/PGHNBYPD4+PEA9CLvvFB3DMhwGjeN7+nB8ivO+Y3cPD49GCL4Tx8BgRHOvLPgpaemPMSWPMV40x3zPGvG2M+aXe55PGmC8bYy73/j9xr309ajDGBMaYN4wxf9T7+4wx5uu9+/27xpjEvfbxqMEYUzDG/L4x5h1jzEVjzKsP670+tpd9hLT0LRH5FWvtUyLyCRH5xd55/pqIfMVae15EvtL7++OGXxKRi/j7H4nIP7HWnhORXRH5+Qcyqo8WvyEif2KtfUJEnpfu+T+c99paeyz/icirIvKn+PvzIvL54zr+g/pPRP6tiHxGRC6JyHzvs3kRufSgx/Yhn+eSdB/sHxGRP5KuOnxLROLD7v/H4T8RGReR69LLfeHzh/JeHyeNP5KW/uMEY8xpEXlRRL4uInPW2oNqlDURmTtsu0cUvy4if0/UbGlKRIrW2oPi9I/j/T4jIpsi8i964ctvGWOy8pDea5+g+4hgjMmJyL8WkV+21pb4b7b7lf+xmQYxxvyEiGxYa7/1oMdyzIiLyEsi8pvW2helKwV3KPvDdK+P82VfEZGT+Hup99nHDsaYULov+m9ba/+g9/G6MWa+9+/zIrJx2PaPID4pIj9pjLkhIr8jXSr/GyJS6MmqRT6e9/uOiNyx1n699/fvS/flfyjv9XG+7K+JyPlehjYhIj8jIn94jMc/FhhjjIj8MxG5aK39x/inPxSRn+st/5x0Y/mPBay1n7fWLllrT0v3vv65tfZvichXReSne6t9rM5ZRMRauyYit40xF3of/ah0/Rweynt9rKIaY8zfkG5sd6Cl/4fHdvBjgjHmUyLyn0Tku6Lx69+Xbtz+eyKyLCI3ReRz1tqdoTt5hGGM+bSI/E/W2p8wxpyV7i/9pIi8ISJ/21rbeI/NHzkYY14Qkd8SkYSIXBORvyPdH9GH7l57BZ2Hx4jAJ+g8PEYE/mX38BgR+Jfdw2NE4F92D48RgX/ZPTxGBP5l9/AYEfiX3cNjROBfdg+PEcH/D8aT/TzyIFmZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = numpy.arange( N, dtype=float )\n",
    "x1 = ( x - b1 ) / b3\n",
    "y1 = numpy.exp( - x1 * x1 )\n",
    "x2 = ( x - b2 ) / b3\n",
    "y2 = numpy.exp( - x2 * x2 )\n",
    "\n",
    "y = b0 * numpy.outer( y1, y2 ) + b4\n",
    "numpy.random.seed( 130105 )\n",
    "y += sig * numpy.random.randn( N, N )\n",
    "plt.imshow( y, origin=\"lower\" )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the model as a 2-d KernelModel using a circular Gaussian kernel.\n",
    "Add a constant 2d polynomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T05:00:50.938481Z",
     "iopub.status.busy": "2023-07-05T05:00:50.938196Z",
     "iopub.status.idle": "2023-07-05T05:00:50.939568Z",
     "shell.execute_reply": "2023-07-05T05:00:50.939739Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2-d-Circular-Gauss: exp( -0.5 * x^2 ) +\n",
      "PolySurface: f( x,y:p ) = p_4\n"
     ]
    }
   ],
   "source": [
    "mdl = Kernel2dModel( kernel=Gauss(), shape='circular' )\n",
    "mdl += PolySurfaceModel( 0 )\n",
    "print( mdl )\n",
    "lo = [0, 0, 0, 0.01, 0 ]\n",
    "hi = [10, 70, 70, 70, 10 ]\n",
    "#mdl.setLimits( lo, hi )\n",
    "par = [5.0,30.0,30.0,30.0,5.0]\n",
    "mdl.parameters = par\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the fitter: AmoebaFitter.\n",
    "The map, y, defines both the xdata and the ydata. By setting map=True the fitter transforms the map into a (N*N,2) array of pixel positions, one for every pixel in the map, which is used as xdata.\n",
    "The fitter attribute 'yfit' returns the result as a map, while the result of the model is just a 1d array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T05:00:50.943764Z",
     "iopub.status.busy": "2023-07-05T05:00:50.943431Z",
     "iopub.status.idle": "2023-07-05T05:00:52.417083Z",
     "shell.execute_reply": "2023-07-05T05:00:52.417262Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y       (71, 71)\n",
      "xdata   (5041, 2)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yfit    (71, 71)\n",
      "mdlfit  (5041,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD7CAYAAACscuKmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfvUlEQVR4nO2dW6wkV3WG/9XV5zJzjG3GCZMRJrEjLCO/gNGIi4wiBzACgiAPyOKSiCAkXkgEChEYXpJIQYIXwA8RksUlfiCAY0BYFoI4xlYSKXIwmATwQEwcELZ8AXvMjGfmXLp75aF7pv7dp1b37tPd5/Q5+/+k0eyu3lW1q6r32f9etfZa5u4QQhx8WnvdACHE7qDOLkQhqLMLUQjq7EIUgjq7EIWgzi5EIYzt7GZ2tZn9gP6dMrMPmNkRM7vLzB4a/P/c3WiwEGJn2CTv2c2sAvAogJcDeB+Ap93942Z2E4DnuvuH59NMIcS0TNrZXwfgr939OjP7KYDr3f0xMzsG4F53v3rU/su24qtYm67FQoiQdZzBpm9Y03ftCY/1NgBfGpSPuvtjg/LjAI6O23kVa3i5vWbCUwohcrnP7w6/yzbQmdkygDcD+Kfh77wvDxolgpm918zuN7P7t7CRezohxIyZxBr/BgDfd/cnBp+fGMh3DP5/smknd7/F3Y+7+/ElrEzXWiHEjpmks78dtYQHgDsAvGtQfheAb8yqUUKI2ZPV2c1sDcANAL5Gmz8O4AYzewjAawefhRALSpaBzt3PALhsaNtTAGRtE2KfIA86IQpBnV2IQlBnF6IQ1NmFKAR1diEKQZ1diEJQZxeiENTZhSgEdXYhCkGdXYhCUGcXohDU2YUoBHV2IQpBnV2IQlBnF6IQ1NmFKAR1diEKQZ1diEJQZxeiENTZhSgEdXYhCkGdXYhCyI0bf6mZ3W5mPzGzE2b2SqVsFmJ/kTuy3wzgW+7+IgAvBnACwE0A7nb3qwDcPfgshFhQxiaJMLNLAPwBgD8DAHffBLBpZm8BcP2g2q0A7gWg/OwlY42ZgheHCdKTH0RyRvYrAfwKwBfM7AEz++wgHdTEKZuFEHtHTmdvA3gpgM+4+7UAzmBIsitlsxCLT05nfwTAI+5+3+Dz7eh3fqVsPsiYTf5v0Tlo1zMhYzu7uz8O4JdmdvVg02sAPAilbBZiX5GVxRXAXwD4opktA3gYwLvR/0Nxm5m9B8AvANw4nyYKIWZBbsrmHwA43vCVUjbvd2YpV23BfbS8N1n96N7sU6v+gj8dIcSsUGcXohBy5+xivzOtXJ9QoltrMazZ3iPJnXMNOVJ/n8p7jexCFII6uxCFIBl/0JhUrmfK84ll+YJY5q01XpbPTOovuLxfjCcihJg76uxCFIJk/EEgR7pnyNORUj1H3mZIfZuDz7mPlMlVXew11zOqEkn0LKmfI+/3UNJrZBeiENTZhSgEyfj9xIws7aFcHyXVaZ9Qircyxo45yHjLlcZWy+xQ+veiezahvJ/UYt9vVPzdDNDILkQhqLMLUQiS8YvOFJb2ieX6KKkeSXSuF9SJZf8Ukj6wrAMjJDq1I5T+kdSfWN5PaLEfat88JL1GdiEKQZ1diEKQjF9EZiXdJ5XrLMNHyPh0n8nON5dAjkOS1yKJzxKa6iRyPZL6OfKejs/PYWKL/VA7ZiXpNbILUQjq7EIUgmT8ojAP6T6pXI+kOgBU5EAeyfKcc0/jkNMLZO+QzE1lNpVZNlMd63Yb2ze5vKdro/OyxX6kj/2cfeuzOruZ/RzAaQBdAB13P25mRwB8BcAVAH4O4EZ3P7njlggh5sokMv4P3f0l7n4+pLSyuAqxj5hGxiuL67TslXQnSZ7UiaQ6ALSavzM+bljOkPE59yKSsKOs8Sz9eyynaTtfW6+W9FnyPjkvHZ8fG1nsQ0nf/xL0ZeM5ppH0uSO7A/hnM/uemb13sE1ZXIXYR+SO7K9y90fN7HkA7jKzn/CX7u5mFmZxBfBeAFjF4akaK4TYObnpnx4d/P+kmX0dwMswyOLq7o+Ny+IK4BYAuNiOLEbkvb0i16FkVtI9keUkJatmSc9yNqkz4liomi37XgWSPnLImZRIqgNDlvbIAt8s79Gt2+oZ8j75QbPFnutkSHpgyuWyGZJ+rIw3szUze875MoDXAfgRlMVViH1Fzsh+FMDXByNGG8A/uvu3zOy7UBZXIfYNYzu7uz8M4MUN25+CsrgKsW+QB92iMM08nefTFmxPysFrOJ5bt4d+GrSPt4P5O2335NUbGut7Mn/HeGjqmnixDc3ZrUuLXJL5ePM83Tr0ii3xIuzQcZpfH5rx8Zvn8jnzdyDT027StNOEfOOFKAR1diEKQTJ+3ox63bZH0t1YoidlkuEs1YePldRrNZYTuV7RK7nAm47XkESvlBJPDt4+5InGMh4dkuuJdCfJXNFrta1auif3m6cA3G5QfW4DNzVpWyDpgTxPu0jSn2/TiDdwGtmFKAR1diEKQTJ+Huwg9NLE3nE50p0kurWbtyeSfKl5e/+7qrlM0r0XSPpExlfjpTtvD6U7K9tOvBCG5Xqr0yzjbSvwCCRJn1js+VzcPGu2lOdIemCHi2cmQCO7EIWgzi5EIUjG7yY27ESx84UtU0l3kuselVdiGd9LyiTjl0iKt+tyrz1exm+zTDfAKjktp9K2tcULXlqN21ubtYRutVkyN5fRIkkfTafYCSe6BioP1wkt9bysPolaO5mzjUZ2IQpBnV2IQpCMnxVhuKXJ/56G4aSC0FKpw0yGdF9ZqsvL7cZybzmV8b3l+hxdKvN2lvHdJZbx9XFSBxvanljmqcxWdy6zO3t3SMazX0yHpTuVqX3VJlnsSSa3Kn4OaCR678K+++6B3B6euySOQmSBz1rbPv53ppFdiEJQZxeiECTjd5FtKZQDC3wYCTZZmprh654h3XsrzdK9uzok41dIxnN52ZrL9eng1KTEMs/qlIedSMYH1vjWkHt6a4vLJNeX64NVm81vCCq2dlO5ynCUSmqwJOfNSf204WnEW745zbnlJnW20cguRCGoswtRCJLx0zAvC3ziqBFEiw0caSJf91C6U7l7iGT8SnoNnVWS6Ct1ubPSvL23XO/bI0mfWOaTCDa0PfKND2V8+hxam3W5Igt8l6Q737IeR+FJpDsaqZo3JySRdHLKSJ2DnBxsLFlSK994IcQY1NmFKATJ+DkT+r8DsQU+cp7hJA5RoMfI1z1wmGHp3jlUH7OzOkLGr47fnsp4p3K93enSvMVmdzpxYo3nyDH1dra+A0CL5HqPyq0NOmwynQjmDRkrlavEUk5f9IIIOYmzzHA+vWCxQCT9myzzI1zks0d2M6vM7AEzu3Pw+Uozu8/MfmZmXzGz5XHHEELsHZPI+PcDOEGfPwHgU+7+QgAnAbxnlg0TQsyWLBlvZpcD+CMAHwPwl9Y3D74awDsGVW4F8DcAPjOHNi4WO4hCEx8qwwKfBJlsjuueBIdkSU/LVCOHGba6s3TvHEqvs3OIyofJGk/SvbtKcn2Ftq+QlXmJdGZFkjRxqmHZSu1gidoleb6Zjlks17tkmW9HS21pCpFKa5bSzc1jCzr76Bu/duBgle3MHHW95jZNY5nPHdk/DeBDqG/3ZQCecffzLkCPAHj+RGcWQuwqOYkd3wTgSXf/3k5OYGbvNbP7zez+LWyM30EIMRdyZPx1AN5sZm8EsArgYgA3A7jUzNqD0f1yAI827VxMymZjy+j4CDT9z5NZ4HNSOEWBIZOlqIl0t+YyyXYA2Fojq3si6clp5RDJ+FWSqCu16bxapuWkFK+9IklvpJOdZHyXpHu3Q9e2OeTHT7K+dY6DYNZ1ovj1rNeNzp28CWAreJfekLD0Zone4TpkmW+lY22SIrrV5S+oHFjmm+NhJowd2d39I+5+ubtfAeBtAL7j7u8EcA+Atw6qKWWzEAvONE41H0bfWPcz9Ofwn5tNk4QQ82Aipxp3vxfAvYPywwBeNvsmHVxshCU/ssBHjjehBT4IDNldabbAsz974ixzOG0rS/etNZLua+Qkcohk+eF6+ebKSu31srpM5aW6zhIt12QZ3yXJvEXlja36p7u+SZ46ADY36s9dvjdJyqxgTS1L98S6zttbzXU6vOw4iB7E8eerdPqRWOerYHqAgAvXE//G5C4rRCGoswtRCPKN30si6c5OFImVPsiSyimYMmK6c0QZXqKaWubTpiZWd5bua7UsXV6rPVjWDtWvWS9ercsXLVN5qS4vU7iZimU8yep1Cn9ztlN7Z5/aSBt7eqn26DlT1fW2yKM7iRHD1nWS67x0NglimUTCoedA99vJAu9bJN1J0g+nf+Lpgef8NiZ0sNHILkQhqLMLUQiS8TlME5EmWsY6/NmCcuR4w1bmIGMqS/dIxkfWePZzB1KHGba6s3S/eG39QvnIobMXypetnqm3L9fb19q1jD9M4WUqCkPTJYeSs7Ru9jTNM55eOpy09an22oVyi6z8p6nOFjnvs7NOGnOetkdBLKmcOPNs8jPhoKG8Njfz9xDV4WnAsMNW065jawghDgTq7EIUgmT8jBgZkWZSIkeaJKb5eMs8Z1LtJqmZ6l2TiDJcpiWqQOrrzg4zbHVn6X708Km6vFIL6Oct19svqc7Vx6F1qRWtZd2kcDZnqVG/adfSfa1KF1i1rTlcS4+s7qfJp71Lcr1LfvYcuLJFU5+KI+9E2WqpjHbwrKr0dxI6ViW/pyC2fHe8c7xGdiEKQZ1diEJQZxeiEDRn30PC8FPhohj2pgtCLFXNc8gkOQMnSKD5e5dyoXEoKQDJmnRe2MLecfyKjefpz185eaH8O+3fXCgfqZ69UE7n7PW5N9E8Z3+qe9GF8pKl89Ueed116NXdRqe+cF5Ic24lWPefvLqk47eby8m9p/l4K3g+2xZGRd5xwW9j0uAQGtmFKAR1diEKQTJ+zkQRZBsq7ricvJKbkbxPosAiDSfFa9J5YQt7x/ErNpbuqYwnbzqjte10mVukVc+0ag+9VRvKDEFs0eu6c/Q+8dRy7XX37HK9fWOZ01nzQiLyGgxesfVoSXqcjjrjeQJZXnAJDbJ/1BE0sgtRCOrsQhSCZPxekiPbAsu8RzI+ShfG5Wr89iSBA4CqXVu8OZwUr0m/uF17xLF3HFvdWbofoVUna+QltkTlLa+nDyte1+95ffz1Kg1Ldbqq5fpaVVvtD7drl7gVuoaqqs+xRdcd3SdE9zW5lxNOy3Lr8W8mI6Iso5FdiEJQZxeiECTjF5Hhdc6NdYLtFkj6oJxKUpLuQzK+Rd9xJFgOJ7XCspycZJKyNUv3i1q1w0yLbMqJw0yvtsBv0ALzU54uhDlM5ztc1dJ9tWqOZtuqeHEJyfgMuZ6Yv5Ny9Bw45NiE1vcpyUn/tGpm/2lm/2VmPzazvx1sV8pmIfYROTJ+A8Cr3f3FAF4C4PVm9gooZbMQ+4qxMt7dHcB5c+rS4J+j1JTN82KaNfCBZEzrNJc9kqHbTtGci61KyrTmHVymKQCdg63uLN0rvhec1ZmmAHy3Kgw5AFGbWnSAVrDO3XLuQc6UaFpmmA68iaxfmJlVZvYDAE8CuAvA/yIzZbOyuAqxGGR1dnfvuvtL0M/W+jIAL8o9gbvf4u7H3f34ElbG7yCEmAuT5np7xszuAfBKZKZsFpl4s8TM25fTDEd1mssWbN9+Cmssd5My5ZADl+s67OueOMxwQ2h7hzxHenSdW3TM7tCYxW3qUb0eta+XXA/tnHH/LLp/Sf0dZCffyT4TkGON/20zu3RQPgTgBgAnoJTNQuwrckb2YwBuNbMK/T8Ot7n7nWb2IIAvm9nfAXgAStksxEKTY43/bwDXNmxXyuZ50cuQ9FGVQNKH0pOOw3nO0E0tw5xIgVMnJ/nXaDnpGYoqc7pX53t+DoVsZV939GrjbctYeteNPU3y/iyFjjndS3O9cUSbs926Tevd+ufeoWvoUY41JKmZ683Way7nTI/SMn3ozVe2DyN3WSEKQZ1diEKQb/xekiPjAtnHctB6gXQPpGcoT3nJ5JCM71FShXUK1sipkzn/2slOnW/t4la9HDWNMFMvd90IHGZYMZ8m6f5Ur04S8Ux3DczTnXpZ6ylqE7eVA052O5Syma47lO7RfU3qBHLdg/Ko72Yk/TWyC1EI6uxCFIJk/Jxxb5beDRUnK5PFPpHx3fHlrLTEm0OOKvR5fbOW089u1pZvTp3M+ddWWs3BIde9Ps7hINcbO8xEceOf2LokOe5TW7WsP7lZt+nURi3pWcb3kvxufD9I0if3zMeWE3mf8zyBySU67X/+dzbqCBrZhSgEdXYhCkEyfg/xwKIebU+cbRIreiTXm6V7IkkT6U4SltIVA0Bvo5a6Gxu1/P5Nu5bGKxQJJkqbvEEW9Wco+GQi4+mVAvu5s4xnizvLdgB4YuPiC+WTG7WMP71e77+xTkEqaYrS2mi+B8l94vvHbzb43ndpmtXpNdYZZY33wALvU/jPa2QXohDU2YUoBMn4GeFsEW+xhK22Vx5HIN2t22vc7r1AMpKzSEVrSyuyMncjqToUZ6R1jpastuufzZlW7ajy61YtmXlpKWdSPdOtpfRz2rX8XqWTc3QZPg77ubOzDFvcgVS6P3WmLp87V+/fO1dfA18bxaccKpNE5/u0FVjjqYxuxvMEYueZaPlzzhoKQiO7EIWgzi5EIUjG58DyioMCsryKAkayHBu2UEfH8ma5nrSjU5uBrUNOISTjE4lJMpQlKfm+oEfRIJP0RRjK8FrV59siGX+a6vdoCelGp96ZM6ly6qhlMmtzYEiOLsNLVNnPnZ1lgNTqztJ962xtgbez9TVU6xTs8hyV6d4k5UTSc7lZoluXTPbBVGzb5+Q3ELyRSfYdb6XXyC5EIaizC1EIkvG7yLBDRBi1hJ0rMqzxqQWeZTyXSZ6S80yb5Xki1YdkfJKVtB4jyL8EWxSp5jQtiWU/9GeX6zorS3UEmyQdEznVcGDIJEIO+efz8YHUYYat7izd22etubxeH6da5+lOMA3apHvPMn6Lp1n8fAJJj/RZD39XV9q5g41GdiEKQZ1diEKQjF8QEn94tt62yCmn1ywHvSLJuEV+3m12FiFLe5sCOlYUeJGkepLRFYhzJHGseIrywtFfzrFf/XItsas2WeDpfMlLCp7pdPlc1O7NIcelDboH7DCzHkj3OpAOqnP1CVnSt1nGb9CUY5PL/BzoGbJ0T6Q6hwYa+kzWeA+WNidcqB9L+5y48S8ws3vM7MFBFtf3D7YfMbO7zOyhwf/PHXcsIcTekSPjOwA+6O7XAHgFgPeZ2TUAbgJwt7tfBeDuwWchxIIytrO7+2Pu/v1B+TT62WCeD+At6GdvxeD/P55TG4UQM2CiObuZXYF+woj7ABx198cGXz0O4Ohsm7bP4DkWvSJKFsX0hv62JrGLrLnc41c1tH+Lyt3mObvRnL1Fc/OqxXNzPm/gKdivSUVOpNAcuonn0b2V+ty+RN53Fedj5nM3nzZN4MDrztO28ufE8+1c8IqN5+k0f2/zq7f1aM7O9z56Hdo8f/fh+Xf0KnYKrzkm2xpvZhcB+CqAD7j7Kf5ukMO98cxK2SzEYpCbn30J/Y7+RXf/2mDzE2Z2bPD9MfRzt29DKZuFWAzGyngzM/STNp5w90/SV3egn7314ygpi2u0KGZHh2qOPBuGq+LFFaSZjSS6k7w3fgVFbSUVjqqKrmFoHGDpnkSzJTlNkWpp6TmncUOPXvs5vTFzPp0FU4YgoQV7BwJpOKmK1p5XgXdc8oqNt5/rUX16TbhB5eh12xbNaZKFMEF56LMHUWinCUuVM2e/DsCfAvihmf1gsO2j6Hfy28zsPQB+AeDGHbdCCDF3crK4/jtSkwnzmtk2RwgxL+RBN2ficFVIrfORZT7xpiMJzJZpXtseWfXpVFVo+R+VxILkemiBp9Mtk5cer5PnxTat5nLqQkeboxxrvBoHadioaO0+L2xJvOMSC3y3uUwy3jZpOrVZnzjLAj+82CWMJBysbU/eAA22j3iE8o0XohDU2YUoBMn4aZgmXNW2Q2VY5lnSR1KcyomkbzXXSZaQkHnchuSgkUNPi8skmbsk3aslWpNO+RiSxTZ8cr59gZ9PIul5dtNJG9sKkmBUYWiuXmOZre6JdN8gub5BMp4t8JE1vkPbh6LGemCpn8YCz2hkF6IQ1NmFKATJ+F3Eh6zdod98kCeN9zaq4xQcygJHn2RrEAKrinLMAbAeSXzOabbEiSjIAs9r5ulX1gscaTwYdkIZT7eo1R0l48dHgo183ROHGba6R9K90yzXvdMs6X1Ekogo11ujBT4TjexCFII6uxCFIBk/K6a1zCdhiNgfPMPZhg9DZurI7TH06B8hjRN/eFpG21umpaxbbGlvLie+8YkjTdQoIpLxQ9b4MG11Ev01Y5kqR4hlh5lIum/VdVLpTs82N0kEkyPdo3xwhEZ2IQpBnV2IQpCM302GnSiSiDbNPuosddmpJrHMR6fjOoFjhgWRS62X/jSSZBUcAYcisjgvtQ0i2CYW+CpwDIpIHIz4zUFarcVLgTmNcjeQ6EFCB3aMSXzdExkfWd2DKLIjosuGFviIDOnOaGQXohDU2YUoBMn4ebCDaDbhUtheswSeWNLTzpZ+UW8ftfyS0iVbu7bAO8t4DoKZBLuMAlw2R8+JSB1sRjgA8ZsEDvzYDeR6kv46iCQTbM+R7tk+7zNynonQyC5EIaizC1EIkvHzZpSkDxxu5i7pOdFyYuFmGT7sVMOOIeQnz3HqyerOce0jq3si+zOwIAgjfISMD6YmFuVfi+R6UidYvpoEB82Q7sPyfBrnmYxlsBrZhSgEdXYhCkEyflGYs6RP5Hri2cIOPM3pofvnJrneoXqhdKf6k0bSiQiW5m6TsMnUJPBFjyznkdML12e5HR4/Q7pvc7Laud97Djkpmz9vZk+a2Y9om9I1C7HPyJHx/wDg9UPblK5ZiH1GTpKIfx1kb2XeAuD6QflWAPcC+PAsG3YgGZabkcPNHCR9Ite5CSwdW1S/SkJRhhLfgsyyiUSvgjElCpoZEUn3bVbtqB5buwO5zvVzLO3R8aeNNJMj3ScMRLlTA53SNQuxz5jaGj8qXTMApWwWYkHYqTX+CTM75u6PjUrXDADufguAWwDgYjsyGyffg0KOD/2sJH0QBcWTqDokw4clZiDxk/05Eg9Z160zoVxnZ5soegszJGezrN+BdT2R64FEz4rpviDSndnpyH4+XTNQUrpmIfYxOa/evgTgPwBcbWaPDFI0fxzADWb2EIDXDj4LIRaYHGv824OvlK55lsxD0rPR2JodaRKrOVvst/nxBxI/cobh9gXLWqdilJwd4bjSVCeU61NY2pOmTusso/RPQohJUGcXohDkG7+IzEzSB/Hr2WLPb00jeT/0XSIqE/92HjvI8Sa6hhx/+IgRVu3YQj5eoofHyZDiWRFldlm6MxrZhSgEdXYhCkGdXYhC0Jx90Zl0/p5sbv5bnryeY6K5/HA7+FhR/jle257uQPWbmzExI1/DZczHk/qTvSab6Rr0OczTGY3sQhSCOrsQhSAZv5+IZN4Ur+eSw0TyHkgkfnIKFumtQNLzOWblQcdtyJW/U0jumb1Wm7NUH4VGdiEKQZ1diEKQjD8ITGGxz5H3fSgfWuT5xsrdMmT/bjAriT7hMdP6ixHGQSO7EIWgzi5EIUjGHzSmsdiPIlvun68+m8QG0zJ3iZ7suxhyPUIjuxCFoM4uRCFIxpfCpPJ+2/4Z8nZCqb/rzChn2qLL9QiN7EIUgjq7EIUgGV8608r75FiLYYGfin0q0XOYamQ3s9eb2U/N7GdmpkyuQiwwO+7sZlYB+HsAbwBwDYC3m9k1s2qYEGK2TCPjXwbgZ+7+MACY2ZfRT+X84CwaJvaYncjZOSxfnSkHWKLnMI2Mfz6AX9LnRwbbEpTFVYjFYO7WeHe/xd2Pu/vxJazM+3RCiIBpZPyjAF5Any8fbAs5jZO//he//QyAX09x3v3Kb+GgX/d2lXzwr3k7e33Nvxd9YdkhfYZ3NGsD+B/0Ezw+CuC7AN7h7j8es9/97n58Ryfdx5R43brmxWLHI7u7d8zszwF8G0AF4PPjOroQYu+YyqnG3b8J4JszaosQYo7shbvsLXtwzkWgxOvWNS8QO56zCyH2F1oII0Qh7GpnL8GX3sxeYGb3mNmDZvZjM3v/YPsRM7vLzB4a/P/cvW7rrDGzysweMLM7B5+vNLP7Bs/7K2a2vNdtnDVmdqmZ3W5mPzGzE2b2ykV91rvW2Qvype8A+KC7XwPgFQDeN7jOmwDc7e5XAbh78Pmg8X4AJ+jzJwB8yt1fCOAkgPfsSavmy80AvuXuLwLwYvSvfzGftbvvyj8ArwTwbfr8EQAf2a3z79U/AN8AcAOAnwI4Nth2DMBP97ptM77Oy9H/Yb8awJ0ADH3nknbT8z8I/wBcAuD/MLB90faFfNa7KeOzfOkPEmZ2BYBrAdwH4Ki7Pzb46nEAR/eqXXPi0wA+BOD8ovbLADzj7p3B54P4vK8E8CsAXxhMXz5rZmtY0GctA92cMLOLAHwVwAfc/RR/5/0/+QfmNYiZvQnAk+7+vb1uyy7TBvBSAJ9x92sBnMGQZF+kZ72bnX1iX/r9ipktod/Rv+juXxtsfsLMjg2+Pwbgyb1q3xy4DsCbzeznAL6MvpS/GcClA7dq4GA+70cAPOLu9w0+345+51/IZ72bnf27AK4aWGiXAbwNwB27eP5dwfo5iT8H4IS7f5K+ugPAuwbld6E/lz8QuPtH3P1yd78C/ef6HXd/J4B7ALx1UO1AXTMAuPvjAH5pZlcPNr0G/XgOC/msd9WpxszeiP7c7rwv/cd27eS7hJm9CsC/Afgh6vnrR9Gft98G4HcB/ALAje7+9J40co6Y2fUA/srd32Rmv4/+SH8EwAMA/sTdD1RQAzN7CYDPAlgG8DCAd6M/iC7cs5YHnRCFIAOdEIWgzi5EIaizC1EI6uxCFII6uxCFoM4uRCGoswtRCOrsQhTC/wNT41p7YzNPRwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fitter = AmoebaFitter( y, mdl, map=True )\n",
    "print( \"y      \", y.shape )\n",
    "print( \"xdata  \", fitter.xdata.shape )\n",
    " \n",
    "# find the parameters\n",
    "param = fitter.fit( y, temp=10, maxiter=100000, size=5 )\n",
    "\n",
    "print( \"yfit   \", fitter.yfit.shape )\n",
    "print( \"mdlfit \", mdl( fitter.xdata ).shape )\n",
    "\n",
    "plt.imshow( fitter.yfit, origin=\"lower\" )\n",
    "plt.show()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T05:00:52.419723Z",
     "iopub.status.busy": "2023-07-05T05:00:52.419424Z",
     "iopub.status.idle": "2023-07-05T05:00:52.427393Z",
     "shell.execute_reply": "2023-07-05T05:00:52.427890Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters : [    3.501   30.025   44.957    9.220    0.300]\n",
      "StDevs     : [    0.154    0.404    0.404    0.331    0.033]\n",
      "Chisq      :   207.455\n",
      "Scale      :     1.773\n",
      "Evidence   :   -87.920\n",
      "Covar matrix :\n",
      "[[    0.024    0.000   -0.000   -0.031   -0.000]\n",
      " [    0.000    0.163   -0.000    0.000   -0.000]\n",
      " [   -0.000   -0.000    0.163   -0.000    0.000]\n",
      " [   -0.031    0.000   -0.000    0.109   -0.005]\n",
      " [   -0.000   -0.000    0.000   -0.005    0.001]]\n"
     ]
    }
   ],
   "source": [
    "print( \"Parameters :\", fmt( param, max=None ) )\n",
    "print( \"StDevs     :\", fmt( fitter.stdevs, max=None ) )\n",
    "print( \"Chisq      :\", fmt( fitter.chisq ) )\n",
    "print( \"Scale      :\", fmt( fitter.scale ) )\n",
    "print( \"Evidence   :\", fmt( fitter.getEvidence( limits=[-100,100] ) ) )\n",
    "print( \"Covar matrix :\" )\n",
    "print( fmt( fitter.covariance))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See also randommap.ipynb"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}