examples/fitting_with_jacobian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting with automatically obtained Jacobian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitting functions works in python primarily with scipy's [curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html), which is essentially a wrapper around several non-linear least squares minimizers. By defaut, it engages the [leastsq](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html) function, which it itself a wrapper around [minpack](http://www.netlib.org/minpack/).\n",
    "\n",
    "For the minimization, the derivative of the residuals with respect to the fitting parameters (Jacobian) has to be estimated. However, fits converge better if the Jacobian is analytically known. Of course you can derive it by hand or ask Wolfram Alpha, but it would be to have automatic differentiation in the fitting framework.\n",
    "\n",
    "Automatic computation of derivatives (autograd) is widely used in the deep learning community, and we can use the tools made for neural networks for our fitting purposes. Some libraries like [autograd](https://github.com/HIPS/autograd) focus solely on the the automatic differentiation functionality, while other solutions like [pytorch](https://github.com/pytorch/pytorch) are much larger frameworks.\n",
    "\n",
    "In geeksw, it is possible to exploit the autograd feature of pytorch (the library I am more familiar with) to easily get the Jacobian of any fitting function. A function `geeksw.fitting.utils` is provided, which takes any function that can be evaluated with both numpy arrays or pytorch tensors. Let's try it out!\n",
    "\n",
    "Start with the imports needed for random number generation, plotting and fitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import the [wrap_jac](https://github.com/guitargeek/geeksw/blob/master/geeksw/fitting/utils.py#L4) function from geeksw. In order to build fitting functions which are both compatible with pytorch and numpy, we import a little wrapper around the two of them: [geeksw.nptorch](https://github.com/guitargeek/geeksw/blob/master/geeksw/nptorch.py). Is provies all the mathematical functions which are included in both numpy and pytorch (`sin`, `cos`, `exp`, `sqrt`, etc.) and passes the inputs to either the numpy or pytorch version depending on the input type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geeksw.fitting.utils import wrap_jac\n",
    "import geeksw.nptorch as nt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the true parameters of our toy model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = np.array([90., 5, 100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the function according to which the data is distributed. We will also use it for fitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x, mu, sigma, tau):\n",
    "    sig = nt.exp(-(x - mu)**2/(2*sigma**2))\n",
    "    bkg = nt.exp(-(x - mu)/tau)\n",
    "    return sig + bkg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the data, include some Gaussian noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(60, 120, 100)\n",
    "y = f(x, *p) + np.random.normal(0, 0.1, size=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting parameters for the fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = [85., 1., 10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try it out without the analytic Jacobian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4.06 ms, sys: 94 µs, total: 4.16 ms\n",
      "Wall time: 4.74 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "popt_scipy, pcov_scipy = curve_fit(f, x, y, p0=p0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now with the Jacobian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 11.7 ms, sys: 655 µs, total: 12.4 ms\n",
      "Wall time: 25 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "popt_jac, pcov_jac = curve_fit(f, x, y, p0=p0,  jac=wrap_jac(f))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare the fits to the data and to the true function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD1CAYAAABJE67gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl81MX9+PHXbO47gdzhCrdALu5DbhWsXKJWxHq2UGu1Vvu1Qv1WsdWvKFZbfy1iUIs3CCiC3Bq5FeQIdwgBAuQigSSbO9nszu+PHObYTTbJJntkno8HD81+Zj87E5b3zs68Z0ZIKVEURVHsj8baFVAURVFaRwVwRVEUO6UCuKIoip1SAVxRFMVOqQCuKIpip1QAVxRFsVPOHfligYGBslevXvUe02q1+Pn5NfvcnJwcgoKCmi1n7v0sXa4lZW29LZZuhzXraOt/J+1xT0d5f4HjtKWt7Thy5Mh1KWXjG0gpO+zPsGHDZEMLFixo9Jgxxp5rjLn3s3S5lpS19bZYuh3t8drWaos13ze23pb2+N04Slva2g7gsDQSUzu0B67Valm4cCEAM2fOZObMmR358oqiKHZj06ZNbNq0qeZHo938Dg3gfn5+xMfH13vM0kHc3PtZulxLy1ryftZss7kcpS3WfN/YelvU+6vt6t6vbid35cqVWmPlhezApfTDhw+Xhw8fbu1zae1zbY2jtMVR2gGqLbbKUdrS1nYIIY5IKYc3fNxuslDqfJWwe47SFkdph6LYK7sJ4M7OHTra064cpS2O0g6gdm7GEai22J72aofdDKEoiqJ0VnY/hKIoiqLU1+nTCPPz8/nss894/PHHrV0VpR1sOJbOsu3nyMgvJdzfg2enDWBOXIS1q6UozTInjbDTD6GkpqYyY8YMTp06Ve9xvV6Pk5OTlWqlWMKGY+ks/vIkpTp97WMeLk68OjdKBXHFrqghFBMWLVrEhQsXiI2NZcSIEUyePJn58+cTFRVFamoqQ4YMqS37xhtvsGTJEgAuXLjA9OnTGTZsGOPHjycpKclKLVBMWbb9XL3gDVCq07Ns+zkr1UhRLMum0ghe2nSaMxkFFr3noHBfXpw52OT1pUuXcurUKRITE9m1axd33HEHp06dIjIyktTUVJPPW7hwIStWrKBfv34cPHiQxx9/nISEBIvWXWmbjPzSFj2uKPbGpgK4LRg5ciSRkZFNlikqKuLAgQPcc889tY+Vl5e3d9WUFgr39yDdSLAO9/ewQm0UxfJsKoA31VPuKF5eXrX/7+zsjMFgqP25rKwMAIPBgL+/P4mJiR1eP8V8z04bYHQM/NlpA6xYK0WxnE4/Bu7j40NhYaHRayEhIWRnZ3Pjxg3Ky8v55ptvAPD19SUyMpK1a9cCVTs6Hj9+vMPqrJhnTlwEr86NIsLfAwFE+HuoCUzFodhUD9waunbtyrhx4xgyZAgeHh6EhITUXnNxceGFF15g1KhRREZGMnDgwNprn376Kb/73e94+eWX0el0zJs3j5iYGGs0QWnCnLgIFbAVh9WhaYT9+vWTkydPBmwnD1xRFMUW1c0DX7lyZYqUsl/DMp0+D1xRFMXWqTxwRVEUB9NsABdCdBdCfC+EOCuEOC2EeMpIGSGEeFsIkSKEOCGEGNo+1VUURVFqmDOJWQn8SUp5VAjhAxwRQuyUUp6pU+Z2oF/1n1HAO9X/VRRFUdpJswFcSpkJZFb/f6EQ4iwQAdQN4LOBj6oP3/xRCOEvhAirfq6i2KS1PyTz7c6NHCwPwMsnUm10pdidFqURCiF6AXHAwQaXIoCrdX5Oq35MBXDF5mw5sJrMH9dyt/YAMa467ukXSmiZE+u/7UOu9q88OinO6PPUzoaKrTE7gAshvIH1wB+llA03LBFGntIovSUnJ4fhw+tPpG7atAlnZ2eCgoLMrYqitNrfPnqItfIoE93L6JoXy4byKGJyj3HdJ5Pk4POUJj3Czf020T+ie73nNdzZMD2/lMVfngRQQVyxmPj4+EYHv1cLNPagWWmEQggX4Btgu5TyTSPX3wV2SSk/r/75HDCp4RCKSiM07u233+add95h6NCh3HvvvZw5c4ZFixaxYcMG+vfvz6BBg6xdRYew9NMn+LRyN8PK3Tma+nsKDPU7DaP8VnM27BjRpa68MHM7/bv9fH3c0gSj+6pE+Huwf9GUdq+70rm1Oo1QCCGA94GzxoJ3tY3Ag9XZKKMBraOPf+v1+uYLmWn58uVs2bKFTz/9lFmzZrFo0SIANmzYwJkzZ5p5tmKON9f8hdW6XQwpd2H5Azvx8e3RqMxB7TwGZ8Xy+/x0rq28m6S0nNpramdDxRaZM4QyDngAOCmEqNm96S9ADwAp5QpgC/ALIAUoAR5pVW22LoKsk616qkmhUXD70iaLfPTRR7zxxhsIIYiOjsbJyYkZM2Zw9913A+Dt7U1RURG7du3ipZdeIiwsjMTERGbOnEnPnj1rT/NZsmQJPj4+/OlPf2LZsmV88cUXlJeXc+edd/LSSy8Zfe3HHnuMixcvMmvWLB599FECAgI4fPgw8+fPZ+PGjezevZuXX36Z9evX06dPH8v+bjqJ/bu3czxvHX2cPFhx72Y8PfyNbnQF8EP+ffRwCuU1l5W8+d5vedTt9/x5+kC1s6Fik8zJQtmH8THuumUk8HtLVaojnT59mldeeYX9+/cTGBhIbm4uzzzzjMnyhw4dqt0v/NixY/zxj3+sDeBffPEF27ZtY8eOHZw/f55Dhw4hpWTWrFns2bOHCRMmNLrfihUr2LZtG99//z2BgYGsWrUKgLFjxzJr1qx6HySKeUrKK7n1rT1Vk42+LsRXLOJljQ6PJ7bj5xcG/DxuvWz7uUaBeY1+MsLvKNuDk/G7dJDFX+q4a1gE64+kq50NFZtiW5tZNdNTbg8JCQncfffdBAZWzRF06dKlyfJ19wuPi4sjOzubjIwMcnJyCAgIoEePHrz99tvs2LGDuLiqbIaioiLOnz9vNIArlpVXUkFJhb42KA/Wr6abJpUf+v2N6V3q7/Nes9FV5KLNjWbct2nvxzPwTQJDvyD1chTfJ+Xw6twolYWi2BTbCuBWIKWkapj/Z3X3AZdSUlFRUXut7n7hAHfffTfr1q0jKyuLefPm1T5n8eLF/Pa3v23n2isNXdOW4e1e9bb21eRwJuIgT5b25ELqIKabeI6x4ZF8fSiDrw/kREgyI3w38FP+XLWzoWJzOnQvlJpT6RcuXFj3tGWrmjp1Kl988QU3btwAIDc3l169enHkyBEAvv76a3Q6ncnnz5s3j9WrV7Nu3braoY5p06bxwQcfUFRUBEB6ejrZ2dktrltTe5UrxlXofz6AIybkfYo0gvScu8jQlpl8zrPTBuDh0vgA6x9zf0XPcrgR8iOhnqafryjtYdOmTbXxEhOn0ndoD9zPz89UjqPVDB48mOeff56JEyfi5OREXFwcr732GrNnz2bkyJFMnTq1Ua+74fMLCwuJiIggLKxqfPW2227j7NmzjBkzBqiaBP3kk08IDg5uUd3mzZvHggULePvtt1m3bp2axDSDq1NVn6SP+1FO+OURre3C3rKhRDQx2dhwPFxQtYhBjytkTae0+xamuq/FYJiLRtPkdJCiWEzdLbdXrlypNVZGbSerOJS8kgpKyvU8++kELnmWUJTyLJVOoS06iafhist3XZcSok0i8a593BrdOP1QUdqbqTzwTj8GrjiWAE9XNLKUYE0uLtpepPj2qJ1sNHcpfMOxbn2yAc1nd3H02xXcEvVKozkTRbEWFcA7yI0bN5g6dWqjx7/77ju6du1qhRo5JiklzpVFvJGdh+aPu9H4dwPathTeqd9UngjryWWxicmXnmFkb6OrmhWlw6kA3kG6du2qTrHvANqSErxkGZnd76B7dfCGqvHthot2SnV6lm0/1/zQihDc3HMKu699x+pvlzNy4QvtUXVFaTF1Io/iMKSUaIszyXF2Qk56uN61ti6Fnz3hBfz0BjJ1GzmXpTKDFNvQoT3wmjRCUIcaK5ZXVFZBmSjDS3oS3qf+cFVbl8J7eHZhtnd/PtGc5+PvvuHl+++zSJ0VxZS6hxpjIo1QZaEoDuNqTioFopg+Pt1xd/Otd63hGDhULYVvSXZKVuYxpm9/gKH5ofz78W14uqoRSKVjqEONFYemqzRQLotxkeDm6tPo+py4CF6dG0WEvweCqm1gWxK8AULD4nhb34MleSlsO5luwdorSuuoAG4j3n77bW666Sbuv/9+Nm7cyNKlVfvCtGVL2abuM2nSJJr7NpSRkWHxjbRMvW7durZGXlEBFRrwc/YwmeY3Jy6C/YumcGnpHexfNKVVy+LHj/oNPUQeZ37Y2uq6KoqlqO+AbaDX63FyarwEuzWWL1/O1q1bazfKmjVrFlAVeGfMmNGqQx1mzZrVpvuEh4ezbt26Fr9ua9Sta0tJKREl+fSjghRDMAZNKQcv5rbLviViwO2879+VjIpPSct7kG4BnhZ/DUUxl831wB/Z9kijP6uTVgNQWllq9PqGlA0A5JXlNbpmjo8++ojo6GhiYmJ44IEHAHj44YfrBS9vb28Adu3axeTJk5k/fz5RUVE899xzLF++vLbckiVL+Mc//gHAsmXLGDFiBNHR0bz44osmX7/unuBvvfUWq1at4oknnuDAgQNs3LiRZ599ltjYWC5cuFD7HL1eT+/evZFSkp+fj0ajYc+ePQCMHz+elJSUZu+zdu1aRo4cSf/+/dm7d2+jeqWmpjJkyJDa/x8/fjxDhw5l6NChHDhwoLbc66+/TlRUFDExMbWHUSQmJjJ69Giio6O58847ycvLqy3/ySefMHbsWIYMGcKhQ4cAausKVZM3o0aNIi4ujltuuYVr167V/m4fffRRJk2aRO/evXn77bcByNKW4ScLKZMe6HBBb5As/vIkG461wzCHiwfJXXtyNCCX9T8lWf7+itICNhfAO1rNfuAJCQkcP36cf/3rX80+59ChQ7zyyiucOXOGefPmsWbNmtprX3zxBffcc0+9PcETExM5cuRIbYBtaMWKFYSHh/P999/z9NNP1z5esyf4smXLSExMrLcXipOTE/379+fMmTPs27ePYcOGsXfvXsrLy0lLS6Nv377N3qeyspJDhw7xz3/+0+SBEzWCg4PZuXMnR48eZc2aNfzhD38AYOvWrWzYsIGDBw9y/Phx/vznPwPw4IMP8tprr3HixAmioqLq3b+4uJgDBw6wfPlyHn300UavdfPNN/Pjjz9y7Ngx5s2bx+uvv157LSkpie3bt3Po0CFeeukldDodxaXXSXeFG/zcG67J8W4Ps2+6h2KNhn1HltNr0WbGLU1onw8LRWmGzaUR/nf6f00+38PZo8nrAe4BTV43pqX7gYPt7Ak+fvx49uzZw6VLl1i8eDErV65k4sSJjBgxwqznz507F4Bhw4aRmpraZFmdTscTTzxBYmIiTk5OJCcnA/Dtt9/yyCOP4OlZFTy7dOmCVqslPz+fiRMnAvDQQw9xzz331N7rvvuqUvAmTJhAQUEB+fn59V4rLS2Ne++9l8zMTCoqKmp/1wB33HEHbm5uuLm5ERwcTFpGJhqPYsqEQC+9692nvY47yxK306XyPzj5JoL2TnXAsdIuzEkj7PS7ERrbDxzsY0/w8ePHs2LFCjIyMvjb3/7GsmXL2LVrl9kfEm5ubkBVb76ysrLJsm+99RYhISEcP34cg8GAu7s7YPr315SG5Rv+/OSTT/LMM88wa9Ysdu3axZIlSxrVuabe17WFeHmBu96ZogZfKNvruLN/7LxEH/cgTvlfx0eTQ6EhyPxVnYpiJnN2I+z0QyjG9gMH7GJP8FGjRnHgwAE0Gg3u7u7Exsby7rvvMn78+BbdxxxarZawsDA0Gg0ff/xx7aHOt912Gx988AElJSVA1e/Pz8+PgICA2nH1jz/+uLY3DtQOOe3btw8/Pz/8/PwavVZERFUg/PDDD5usV3n5jaqtX/X1Uwfb87izjPxSMvMncnNJKZNdD9Z7XFE6kjmn0n8ghMgWQpwycd1PCLFJCHFcCHFaCNG6A42tpO5+4DExMbXnYS5YsIDdu3czcuRIDh482Ko9wefPn8+YMWOIiori7rvvblUAnTdvHsuWLSMuLq7eJCZU9Ua7d+/O6NGjgaoeeWFhIVFRUS26T1NqesePP/44H374IaNHjyY5Obn29zF9+nRmzZrF8OHDiY2N5Y033gCqAu+zzz5LdHQ0iYmJvPDCz/uHBAQEMHbsWB577DHef//9Rq+5ZMkS7rnnHsaPH187tGWMQYJelOEqIdAvqHYvcCeNaHGOd0uE+3twvmwET2W58Ev98XqPK0pHanYlphBiAlAEfCSlHGLk+l8APynlc0KIIOAcECqlrGhYVq3EtC9HjhzhmWeeYffu3dauilHX8gpxLr+IcPOjS0CvDnvdmlWdT8hPmem2lTtK3qLSOahdPzSUzq3V+4FLKfcIIXo1VQTwEVVdNW8gF2h6QFWxeYcPH2b+/PltWlxjKXklFVzTllGhN+DqpCHEzx1/DxdkaT5dMYB3SIfWpyZIf7FzACtCEhmcuYNf3vqqCt5Kh7PEJOa/gY1ABuAD3CulNDT9lM7JnvYEHz58eG2miTXllVSQnleKofqbYoXeQHpeKWU6Pc6ikAqNG64uHT90MScuglnRi5j64Wo0PmcJU8MnihVYIoBPAxKBKUAfYKcQYq+UsqBhwZycHIYPr/8tYNOmTTg7OxMUFGSBqtg2tSd4y13TltUG7xoGKdEWF6J3lUgnD1p20qjlaJycmegRwRaRxrbjyYzpM8ZKNVEcRXx8vKlMPaOTQZYI4I8AS2XVYHqKEOISMBA41LBgUFBQs/tvKEpddU+Zr8tVU0Ap4Odp3Q/+W/rMYP3ZeC5e+AyDYbQ69Fhpkzqn0NcjhLhurLwl0givAFOrXyQEGABctMB9FaU2s6Qhg1M5bhLc3BrvPNiRRsY8jJfBgHQ7wtErec0/QVEsqNkeuBDic2ASECiESANeBFwApJQrgL8Dq4QQJwEBPCelNPppoSgt1aOLJ3mlOupmSzmJSqRzAD4aVyvWrIqrmw//cBtA8JUjrD6RzvBeza/kVRRLMScLpcmjR6SUGcBtFquRotTh6ebMjjPXak+TD/VzZ7DTfzgYcp4tU+Lx8bT+AcPjouZD8ndknN6LnDlEnVqvdBi1naxi8+bERdSm6H2dmE6PLy+zyNOb7t1tZNKw36184eOLi34TKdkP0C/EusM6SufR6ZfSK/Zl39ETxGlS6B11l7Wr8jN3P3Z0DSU54CrfJbV8uwRFaS2b241QUUwpLNNx/fp/eSY4kFcG3Mb2Y+m1Qyvh/h48O22AxRfTbDDzNSYGD+X1nAP8dGo3j03sY+ROitIy6lBjxaFsOJbOJz/eQbFzJWcvvIagahlwjZYeUmzO65l7EPKl1F3M2v0kMdcG8P9+/xkBXtafYFUchzrUWLF7n+45QIpHJUFFVUvnG3Y9LH2Iw7Lt5+oF76Zeo1ePCUToJZXeF9mdnGOxOihKU1QAV+xCYZkO57Kv0AtBVsE4k+UsuaWrqXsZe1xoNEz06omnpohdp69YrA6K0hQVwBW78O3Za+B9jqBKA8mlpk8csuSWrqbu1fDxDcfSGbc0geTjo/kg6xqF5/aiM7GCVFEsSQVwxS58m3iJcRU3GJgfiqm3raUPcXh22gA8XJyafI2acfL0/FJ+MAymXLowWh7hX9+et1g9FMUUFcAVm1dUXonm4rc8oc2jKP8X9a7VLJmJ8Pew+H7cc+IieHVuFBH+HggTr1F3nLwMN57378WmHqdZdSDVYvVQFFNUGqFi8xKSshniup8S5y7cddvdXNh5oV1TB+uqu4jImIbj4VmV3bjkdgFfeZaqjToVpXVUGqHiEB7/aB/HKn/LdI9wXrxvp7WrU8+4pQmk1wniPdxOk9f7Y6KvDWDpbz6iexdPK9ZOcRQqjVCxS6UVerLT1lHkpGFCL9vr0TYcJ79SfhOhOgOV3heZ+f/2seFYuhVrpzg6FcAVm7Y7ORsPryO4GyRjYn9t7eo0UnecHECgIbS4C5c8yikpLeDpNYn0WrSZcUsTVDBXLE4FcMWmbT6eQZp3LmOd/HD3CLB2dYyaExfB/kVTiPD3QAK6glgW5BcQ45Rcu9goPb+UxV+eVEFcsSgVwBWbVabTczl1E9edNUztNtHa1WlWzYTmqeIJPJJfxCRxtt51S68UVRQVwBWbte/8daaWXuSjzBwmDf2ttavTrJoFPkV4sp++dPE60aiMJVeKKopKI1Rs1tZTWfzO+Si9Q4aj8e9p7eo069lpA2o3v3rftyuJXdMIKMoiTx9aW8aSK0UVx2ZOGmGHBnA/Pz9TJy4rSj0VlQZOJu9ibZciHuo7gXBrV8gMNfniy7afI7toGDIwnf7eezmovQew/EpRxbHV7eSuXLlSa6yMOpFHsUkHLlwn2P07PvPz4cHek6xdHbPVLPzRV97M+I834OKVDNqqVZztvehI6XyaHQMXQnwghMgWQpxqoswkIUSiEOK0EGK3ZauodEZbT2ZR4nOZAQYNEREjrV2dFnNydmW0cxcyvLT8YkgIz04bwLLt54hUKYWKBZkzibkKmG7qohDCH1gOzJJSDgbusUzVlM5Kpzdw8OxhUtz0TA0YbO3qtNrY0BHkOGtIuZjAovUnSM8vRaJSChXLaTaASyn3ALlNFJkPfCmlvFJdXh0KqLRKzbas/Z7fSrjbNqQQTBl8v7Wr1Wq3Rj/Kl2mZjC2/TFll/e1lVUqhYgmWSCPsDwQIIXYJIY4IIR60wD2VTqbutqwA/TSp9C8zcKYg2so1az2/kCH09u7BBI3x0UeVUqi0lSUmMZ2BYcBUwAP4QQjxo5QyuWHBnJwchg9vtB8LCxcurE0vVDqnutuyelPCXwov8GH+NN7IO8+dw7pbuXatl9RrGAmaBFwzSqmQ9VMIVUqh0lB8fLypTL1AYw9aIoCnAdellMVAsRBiDxADNArgQUFBqN0IlbpqTn2vu6PfeKfDuIpKtuuH230vNTu4P1/n/8iA/B84WTyl9nGVUqgYY6ozK4S4bqy8JYZQvgbGCyGchRCewCjgbDPPUZRGwyY1tGEJzAsL55jsZ/e91BFDfoWTlPh7H6eLp6vJgyEUpTWa7YELIT4HJgGBQog04EXABUBKuUJKeVYIsQ04ARiA96SUJlMOlc6ppqdd9yAGY6e+u4liLniV0L8wADcXF7vvpXr7hBGFKwWe2UzvFcr/3Rll7SopDqTZAC6lvM+MMsuAZRapkeJwanraNcG6Jo2uYfAGGOS9ixSNBqEb6TC91LH+A3kn/wRlySeRcghCiOafpChmUJtZKe3OWE+7VKfHyUgg8/JJxNNg4L0n/u4QwRtgTJ87CK/UE1z2A6k3SqxdHcWBqACutDtTE5F6KeudZuNEBaneWkbIAFzdfDqqeu0u5qZfsiW7kNt1V9mTnGPt6igOpEMDeM1uhAsXLqy7y5bi4ExNRNZM5tWcZjNSk8SL12+wYMi8jqxeuxPOLmgiJzDR5TR7kq9ZuzqKndi0aVPdrBR1qLFiHQ3HwKEqja7uGPf/rD3OiNN/55cu+xHPXgBXxzoM+EDC/7I49Utc0haw5S9P4eqsvvwq5lOHGitWU/fcSGNpdOWVenacTiPX/wSZfSY6XPAGCO99K7lOToS47+fI5TxrV0dxEGo7WaVD1Gyzasze5Ov0dEpgeYAbvSMGE47xtEN7ntTs2WM8YXqJ9Eplz/kcxvTpau0qKQ5A9cAVq9t0IoOu/j/hbpCMH/a7egt8HGX3PqHRMNYjnIuepexJst92KLZFBXDFqkoqKvnuTDqpXjcY7+yPp2egybRDe9+9b0y38RRrNBi035JTWG7t6igOQAVwxap2nL5GL+fd3HDScFuPWwDTaYf2vi/K6KgHuF9byChS2XtepRMqbacONVasav3RNIb4nOCGwcCEYY8DVWmHDfdHqXncnvn59+I5996ckCn8NzmHuUO7WbtKig1ThxorNi27oIwfUrJJ9E7hpdCRuHkHA/VPd6/hKLv36ftMRn/j3xw9n4zBEItGo5bVK8apQ40Vm/Z1YgbDxWm8ddch6ueT+Oqe7u4oWSg1Tgb14qHwYAak7+RM5nRSsoscsp1Kx1ALeRSruf1fe4lwXYyzJod3Hj6KcPOydpXanU5XwvhPRtCv0J+9GX9BAHX/BTZc4KQooBbyKDbmbGYBF7MyOeWeR1efiE4RvAE2n8qjd6k7GZ75gIGG3SdHyLZROo4K4IpVfHUsnVjfHRRrNNzR7y5rV6fDLNt+DlEUSbaLhh6uSUbL2Hu2jdJxVABXOlyl3sCGY+l4djlBV71kZOyvrV2lDpORX0pq0VgAuvv8YLSMvWfbKB1HpREqHW7XuRxKi9M5HlHCvZ6ROLu4W7tKHaYqRXIA/0kvR1Oay/4G1x0l20ZpO3PSCNUkptLhfr3qJ7qlreEmj88ZNmM5/QbMtnaVOkzNNgHPy3judNpHbPlKdNX9qAiVhaKYYGoSU6URKh0qI7+U789lsyvoKD2cQ6D/LGtXqUPVBOfd26LJ8TzMIO1+po1fyOOT+1q5Zoo9anYMXAjxgRAiWwjR5EHFQogRQgi9EOJuy1VPcTRrfrpKqPMFDspzlAy5Ezrh+ZBz4iJ45cnf8ZGfL36+P6kzMpVWM2cScxUwvakCQggn4DVguwXqpDioSr2BLw5fZVj4dv7eNQDtwNutXSWr8fQOYZjw4IbXdXady7Z2dRQ71WwAl1LuAXKbKfYksB5Q70TFpN3JOVzTFnHG7TJjhRdh4Y2G9DqV8YGxXHGFq+lHKCjTWbs6ih1qcxqhECICuBNY0fbqKI7s80NXGNP1O7KdBHf16TwTl6aMH1R19mcv7z3sP3/dyrVR7JElJjH/CTwnpdQ3N5aXk5PD8OGNe111Du5UHNTV3BISkrKZ0ucgXQySSSOesnaVrC6y52QGVEoiXa7w/blsbo8Ks3aVFCuLj483teFfoLEHLRHAhwOrq4N3IPALIUSllHJDw4JBQUFYKo3Q0Y7ccnT/3Z9KkMijQp/HLP/BuHSSpfNNERoNa4NvoTxtNbecy0DKaDWh2cmZ6sysN5mGAAAgAElEQVQKIYx+RWvzEIqUMlJK2UtK2QtYBzxuLHhbkrEjt55ek0ivRZsZtzTBro/eckTaUh1rfrrC/0ac4KPMazw14VVrV8lmiP7TcJeldC9O5HRGgbWro9iZZnvgQojPgUlAoBAiDXgRcAGQUlpl3NvYkVs1y5Fqzk8EVI/cCox9M7pWUEZxhY7RFdugx1icgwdau5o2w9DrZu6NCMOjeAffnb2bIRFGF9wpilHNBnAp5X3m3kxK+XCbamOm5jb7qdnRTQXwjlXzzajmwzU9v5RF60/g7uLErB57me5h4L3+U7mshr9qadx8CHD1JU1c59uz13jqln7WrpJiR+xyMytzNvtRO7p1PGPfjMoqDeSX6ih0T8BPwkXXWQ534nxb3RwUy1UXuJFzhCxtmbWro9gRuwjgS7cm8dq2JFKyi4CqI7c8XJyafI7a0a3jmfrQ7O1+jMNO5cwPHMab36U55InzbTFx8AMA9PbZxXdJ16xcG8WedGgAr9mNcOHChXV32WpWen4p8Xsucsubu5nzn/0Ullfy1xk3EVEdpBvO26sd3azD1Idmt8AteBgk90x4yWFPnG+Lnj3H01sv0Pmk8u0ZFcCVKps2baqblWLfuxFmF5bx9bEM1h1J49y1QlydNNwyKJi7hnYjv0THmzuTLT6mqlIVW6bhGDhAoOYa9HuTOR7d+d95Wxm3NMHoifMR/h7sXzSlI6trUzZvfBR5bjsvav/F3hfm4Omq9plTfmb3uxEG+7izYEJvfjM+ktMZBaw7ksbG4xlsOZlFoLcrs2MjmDs0gkFhvhbJpTU2IdeZs1vM+TCrexhxTZB+rdtx+qZn4vbQe4BjnzjfFnfEPQZH17PDcIy9529m2uBQa1dJsQN20wM3pqLSwO7kHNYfSeO7pGvo9JKBoT7MHRrB7NgIQnxbf1BAW3uK7dV7t8a3AmM966YO39UbJLf/aw9CX8o2w+OIbiNh/mqrtsHmGfRcWNafrwzd+Lf2z2pvcKUeu++BG+PqrOHWQSHcOiiEvOIKvjmZyZdH0/i/LUks3ZrEzf2CmBsXwW2DQ+p9JTUngLRlrNbc3ntLA5m1vhUYyy5pKlXzmxMZJF8r4vGYNTydJ1g66rfU/SidExehAlMDG45n8aVnCMne13HRlpGeT6f+xqeYxy6yUNizjDWbF3ApdZfJIgFerjwwuidfPT6OhD9N5PeT+3Ixp4g/rklkxMvf8swXiew7f50vj6SZlcZmakLOnOyWpgJeDWOrSZtLpzPnvu2hJR9mOr2Bt3YmExtayTdlRynx6op7n8ntWj9HsGz7OYoLoinWaLjJew+gsnOU5tlFAE+/so+Xr//IrN1PcucH0fznq/s4l/wN0mBoVHbDsXQeeP8Q/05IwWCQPDmlLzNjwtl55hq/ev8g/7PuuFlB0FiqorljteYEvNYEY2tlcLTkw+zDA6mk3ihhVMhn5GkETw5/tl3r5igy8ks5XTQJN4PEz+dovccVxRS7SCOM+NXX7Jz2MYuCx+OvcSVee5K7f1jMzhVxsOOvVFzej0Ff2ahXm6Et4729lxjduys/PX8L/5k/FIOJIf+G/1DmxEXw6twoIvw9EFSNfZsa823InIDXmmDclm8FbWHsw0xQ9a2h7t4zmdpS3tqZzG19DWwqO80k4UPU4HvatW6OItzfgzLpQ98SNzK8cgFD7eNK5+RQaYR13biezK5j73Jr1gV8L+3jQ293PvT3p1dJF67nx3GmaCKVuNWWrzvxaGpy0tVJw0uzB/OLIWH4ebq0qX7mTPq1ZpK0pZOJllQzXp+eX4rg571n6tZhx5ksvjubzW+jP+C/pUmsG7uUAf3uaNd6OYqav9tY71WcCTlFyMX7uKwfymt3R6sxcMXkJKZdBvB6yrQcOPwOay99w77KfMo0Ah+9gT7FPrhnTmOvIYYy3Gpn9YFGQVBTnXVY0zuPivDjd5P6MGVgMO7NrPg0pbkJytYGY2tncJj64Onq5cqN4gr+Oqkrdx27i329hjHjvo0dVi9HsOFYOiu27+Oj8sdZq/sFR/r+gQ8eHmHtaik2wHEDeB2Tl35NV8M3uPucxN8pl39nX6NUuvIX/z5klvfhUul0Fs+p6t3WBEE/DxeKKyrR6Rv/HnzcnJk2JJTZseGM6d0VZyfLjjhZMxi39rUjF23G1Dumd5AX23t8gsvZDfD4jxCoTlpvlY/mcO3KOX7p+h92PTtZ7RGuOGYaYUNPTRvO4i/dKC2YgzOVzNckMdXpICd9znGtywmc5HE+O/gmc3qN5cuHFxASGsu4pQnklzY+j7CrlytTBgaz7VQW646kEejtxozoMGbFhhPX3d8i/6islU7XlnTEcH8Poz1wgIVRR7kvdS//HL2Abip4t9rFvhP4v5LTuGbt5UzmMAaHqy1mFeMcqgcO9cdqawgq6e9xmBDfA2R7XyPdVfA/N/J4yKsPS6/2ZrtTXy6VR1N3TlcAl5beQZlOz65z2XydmMF3SdlUVBro0cWTmTFhzI6NoH+Ij9l1spWFK21ZpGRs6Adgan9PrvEHDMBX8/fi5u5vySp3Knm5F5i8cTYxeSFED17On27r3KtUFQfpgZu7nHtOXES9ICVx5lzpaM6Vjiai3ION8zT4Xf4Bzn9HlN82Pg0OpL/uc0KLAsktGM7ZknGE+lf1etxdnJg+JIzpQ8IoKNOx/VQWG49n8M6uC/zn+wsMCPFhVmw4M6LD6Nm18TFhtrgkvy3piA2XywugR1dPhgasZEUBvDv4dyp4t1FAlz6MFJ5c9rnGlRPpPHNrfzWMohjVoT3wfv36ycmTqxZ1zJw5k5kzZ5r93JZO+plbfu3e3Wz8cTkVXsmc99ShEwI/vYF3ZAxRY+Yh+0xBuHk3un9OYTlbT2WyMTGDw5fzAIjp5sfMmHDuiA4jzK8q/csWN2+yRJ3KK/Xc9c4BruaWEj9Xx++OPMctzl15/Ve7LV3dTmn9zj+xJGMHoZfm8u/H/ocBoc1/01Mcy6ZNm2rTrVeuXJkipWx02ofdDKG0Nu3OnKGLmnLagkxGBO6hh38yL2VfgDItfw0KItcniD5u0XybOpqzecH4ebggBOSX6Aj39+A3N0eiMxj4+IfLXM2rqqOrk4aZMWGsP2p8ZWXNEE1b6t1abU1H1Bskf/j8GJtPZvL+/EEkHLqTfU6VrJ+xlsCgmyxWz84sL+8ik7+eRUx+CENu+g/PTlPH0HVmdj+E0pqv/eZOEtYvd3/Vf/Q6uPIDIYff4mDhefbIAxB6gGh/QbC2F1fybieP7qTnl/L69nPcNSyC60UVtfes0BtMBm8wvUCjI4Zc6g6DtPRDwmCQLP7yBJtPZvK/vxjI1JRXmZCWyuW5y1XwtqCAgN7MqfSlhy6Nv32fwldH0/nz9IEqJ1ypx24CuKnsh3ZbqebkApETeCJyAp+9+i2BZQfp5rufQu80bnI7zkq33aTKQJYEdCOvKJZ1BydSJt0a3cbNWUNFpaFe6p2LRvDYxN5GX7alG0e1VmsyYKSUvLz5LF8cTuMPU/sxUP8h2lNr8Zu0mL7R8xuVt7XJW3uy4Vg6lVdu5lHNu3wlLnFS29vqcyeK7Wk2sVkI8YEQIlsIccrE9fuFECeq/xwQQsRYvppt25ukrTK15Vwqj2Vvzu9JvPQqn2W+zHO6Bex2juCU3w0udk/Ar/9fGdv9RUb7f4a/5udTVSoqDbz5yxiCfaqCu7NGoDNIXth4ml+++wP/3X+JTO3PH0y2emKNTm/gxY2n+WD/JR4Z14spYUd46vKX/KPnTTCh8X4nrdmsS/nZsu3n2FgxjBu4MtprB6A2t1IaM6cHvgr4N/CRieuXgIlSyjwhxO1APDDKMtX7WVu+9reEsV5jw97/dbqwRj8ZiifjnlzIIO/dePmc4KpXPje8T/C+/BZNaQ82ips44zWEObG3c+fQbkBVLzb5WhFbT2Wy9WQWL206w0ubzhDTzY/bBocS5ONGdmF5o3pphCBy0War9GRziyt4/NMj/HgxlwXjI5nT+xwLfnyFYCl4es7HoGncD+iobxKOKiO/FIk3TwX2ItXzIi7JZehwt/oHuWJbzJrEFEL0Ar6RUg5pplwAcEpKafRfaEfkgbeFqcm9u4ZFsP5IeqOA1LDMl0cu0835ENN1V7hFk8j2wDw+8vMlTC+Z4NmdCZG3MSLqITw8u9Q+90JOEdtOZbHjdBbH07QAjfYaMfZ6TU04WnLo4tiVPJ78/BjZheUsnRvFYJ/j/Hrfc3hKWDV9FeHhjeZVANMrNpuavFV+VjNpH+ezmZRuexmQPprDBXM6/dFznVVHTWL+Gthq6mJOTg7Dhzf+B19nxy2rMtVr/D4ph1fnRtVbfl83C6UmQA7v2YVl2z35f6WlrPd/kD/EwpLyzezJ+pGNpVdZk/QBwadW8q3bIET/aRT2HEOfkCH8fnLfqt0Ot57lWkE5ThpBpaltEzHekzW12VRrJ0FzCst5fVsSa4+kEebnztrfjiHa+QoPffMsbk4a3r91pcngDVaYs3AwNUfPnSi8lYjK3bj4HUVTOKfTHz3n6OLj44mPjzd2KdDYgxbrgQshJgPLgZullDeMlbFmD9ycXml79horygs5fOIjTp9J4PbUU3STWczoFoYULsQ4dyM5ox8nC26mXFYtBnJ31jBvZA9WHUg1ec/EF27F39PV5OrIusztuWVqS1l96Cof7L9EmU7PozdH8uSUfrif/QLnb54hx9Ofkrvi6dlzQpP3sebOiY6i5j3bx/3vnPTPwenC03y3+ME275ap2J927YELIaKB94DbTQVvazI3Na89e42ubj5kO8/hzZQ+vKLT00uk0y1vCzrvy2zXXKUiPA2f0AR+meuBa+4I9uij2HHahYgm9h6J/dtOBob6cCW3pMngDU1Pgmbkl/JTai6bjmeQkJSNBKYODGbxL26iu6/k9a9mkpubzD+6DSPo7v+Cd3Cz7e2oOQtHVpMplJxSzl37/0yc1zY2nriFB0b3tHbVFBvR5gAuhOgBfAk8IKVMbnuVLM/cCbX2PjG9bj1SZQSpNxbADXAXhQz02oev90mG6LK51+VTLro482hwKNGaLkQ6deeCdjSZlX2Aqt75gvG9cXXW8NPlPJKyCpt9bV8PFz7+8TIABaU6srRlZBWUcSajoPYDIsjHjd9N6sO8ET3o5u/OtweW8mTy51x1gkdCYjDM+hwnl8apkqaosy8to3/f29mY8DKVugs8dfAKvxrVQy2tVwAzArgQ4nNgEhAohEgDXgRcAKSUK4AXgK7A8uo3VaWxrr41mZua1969RlP1KJM+JBbdDkW3swf4F9eJFj/Qo/wYxz3yyA3Oh+CT9NEZmHe9O8MH3s5N0eFsyfBn9U9FZr22tlTHXzf8nAnq6+5MmJ8HMd39+M34SEb06sJNYb44YSDr7Fc8tPFVjokK+iJ4d9BvGTvi95b4FSitFBn9AGx7Dn3WKY6nRRHbXe03o9jRUvq2sJX9SEzVw9/DhfJKg9Hx4tkxYaRc3MGh819z5PopXriWiX9RDu/5+fKevy/dy1zwKOtCWVkPrpUNJL18IJW41E5khvu588SUvtw6KBQAicTbzRlP158/u8vLtJxN2YxM3U/cue8oKsrinu4R3Oo1nq8uziFNq1dDIFYmi66z5MNRXK3sRmDE27x+d7sst1BslN0vpW+L9h4aaWs9lswaDJju+ffrO51+fadXLfKXEm6kcHzVG/QrOE2+h5aUgGwqRQ6uhsOcSU3nmiaEzWHBZPm4E+DqS+WNLnz/gyc+Tu7c7n8TVBTxadYPnC3N4kL5DZKEjkohiCqv4E3Pgfy/sl9xIXkQSbghsZ1dFDsz4R1IZUAvTpVfQ3/yNAUzBuHrriYzOzu72Y2wrWxlWbel6lE3Y8aFMnq4nybM9QKjS+CJGMHfCk6yS5STp4HK6vHSbjodW9MyAXgiJJgkNze6a9yJ9ulJdOgIbsjxvLi90CLZLIrlnT67nnmHlhCT3Y/bbn6TB8f0snaVlHbkULsRKvWZOywkDQYKCq5QXl6ERkCgqx+4eoFHADSYCDN1z7oaplTaygdjZ3H/f+O4oddB4dts/aM6bs1c9v4+NTWEYtlDHtvBhmPpjFuaQOSizYxbmqD20qhm7t4wQqPBz78XwSFDCAweAv7dwbNLo+AN5u23UjelUu130vHm9bqDdBeBW/HXHL2SZ+3q2AVHfp/adAB35F98W82Ji+DVuVFE+HsgqOp5t3WRTHP57g0/IJpKz1Tax7Sxz/GbIh0PyuN8sC/V2tWxC478PrXpSUy1IVLTLJ1nbWyStSabJcLI105b3TnRkbm6+fDUwPuRe//BP0+d4MqNgfTo6mntatk0R36f2nQAd+RfvC1qaR682u/ESoY/yk+HlzPC61Oe/+omLl4vttux3Y7gyO9Tmw7g9vCLt/fJkYZa0qu3lfTMTscvgk3dB7OnPBNSErmur9qqWKV6GufI79MODeBarbZ210Fz0ght/RdviyfOdyS130nHqdtR8PNwIcx5BrqI94gJWs2erP+pLaeGGBuz1/dp3TRCwM9YGZtPI7TlHq6trPBUHJup3SbHdfsb572KqEh5hnx9aO3jas91x2O3KzFteUMkNUavdARjk/kA6TlzKPP+lNjAz9lz7enax21piFFpXzadRmjrTP1DUf+AFEsy1SFILY9mlNaLKYYUfKna1MyWhhiV9qcCeBtY86BlpfNoqkOQlvUQDxfmscB5C77uzurAjE5GBfA2aI/FNIrSkLGOQo2zsicbDWMJ9f+eYKcLTB7Y/GEbiuOw+UlMRVEaZ6HUPZP1qTGVLLv4DP2LPYnr9zH/o74BOhybmMRsaRqhoihVmpvMz8kdwgrNaYoOf8hDY18kyMf8k5MU2+QQaYSKojSvpOQ6M1ZPwkenYXDgh4zsHWyz6bdKy9lED1yxTbaca6+Yx9MzkKcjZ/OXKxvxTFnG+qMPUK43AJ1vgVlnoiYxOzm146PjuGPCS0zTu/Nrwz789DfqXXOU3feU+poN4EKID4QQ2UKIUyauCyHE20KIFCHECSHEUMtXU2kvjrzVZmejcXLmjVmfM66slKUuK4H6w6NqgZnjMacHvgqY3sT124F+1X8WAu+0vVpKR1GrSR1MYF+WO9/PqcBLjAn4tN4ltcDM8TQbwKWUe4DcJorMBj6SVX4E/IUQYZaqoNK+1GpSxxP5i6fZ49aFi8EnCHO+ANj2AjN16lbrWWISMwK4WufntOrHMhsWzMnJYfjwRhOpLFy4sDa9UOlYtr7jo9Jyc4f1pLzgZf6Z+hyh3d6jIuvv/HV2rE1OYHb2HT0bio+PJz4+3tilQGMPWiKAGztV1WhuYlBQECqN0LbY61abSmP1s4k8ebDvHazQbCWu6+v0CfzG2tUzSp26VZ+pzqwQ4rqx8pYI4GlA9zo/dwMyLHBfpYPY8o6PinmM9WT/fXwKd/dNYov3RdZ//jy9n3obLzfbyhxWczBtY4k0wo3Ag9XZKKMBrZSy0fCJYn1qrNFxmerJJmQ8wcrKnjxfvJoPPvmIjly4Zw41B9M25qQRfg78AAwQQqQJIX4thHhMCPFYdZEtwEUgBVgJPN5utVVaTeV7OzaTPVmtjqEPfE6hVw+Ccl/ijbVvd3DNmqZ29GybZr9PSSnva+a6BH5vsRop7UKNNTq2Js+PdffF8zdfsurLWWiL4um62YNH77CNpAE1B9M2ai8UB1czsWXsHzeo47cchbFj1zxcnHh1bhRQFSD1Rafx7rWCYg14X7mHEvcJKljaCZvYC0XtRtixTJ2lWJcaa3QMpnqyQJ33QG/CryzEs0c8hT3XEpaRw+IvZ9Z7vmI71G6EnZypQ5dr1PTQ1D9ex2XsPRDifIngHu9yV1Eel2/MYYf3nexfPNVKNVTMYRM9cKVjNZWKFaHGGjsFY++Ba5WRFFx6nlDneB52+Rhv/XlKSmLw9DS6VkQxg7V29FS7ETowU8MjEf4e7F80RQXvTsDUe6BU+vA73dP8Qz+TL8OvMvfzSRw78VkH184xWDPDSwVwB6ZStJSmztOUaFgp7yem9C50UvLw0f/jH2vvoqQou4Nrad+suaOnGkJxYCpFS2n4Hmh4nuaz0wYwO3Y6q3bPJOHkH1klktn8xRS+in4Gv6EPg8Z6fTx7OWjEmqtJ1SSmoigAHLqUyz/X/5O+Tl/zN20qhtAYzgy7j0FDF6Bx6ti+XlNpkbYWxE0lC9QMVVqCqUnMDg3g/fr1k5MnTwZUGqGiWJuxHu4tg0J4bctpCn9azf1e6/lNhCuRBg2P9prB9HGLcXXz6ZC6dURQtJT2+rCpm0a4cuXKFCllv4ZlVA9cUTqh5oLOgZTrvLDhKF11n1IceJTLruCvNzCgOJSUa7PReQ8xa0ijtcMgkYs2G93S1FYXnrX3cI9N9MBVAFcU22BOD1enN7D60BXe2nmW7ppteAXs57xnKduuZnC1MpLPXWMYOvEh7h0/yehrtKVnak898I5gKoCrLBRF6YTMmXhzcdLwwJhe7PrzLVys/AUHrr6EPP8M8eX34kIl2YH7eOXCEzywaijvb3qEi5cSkAZD7fPbkp2hMqjMo7JQFKUTanLzqwZ83V0oKqsE4LohlHeZybv6mfTMPEms7y5KAzL5Z+5h/rnnMLdXwOuhU6HvVAq1+UCXRvdrKjuj7lCEn4cL7i6aehkztjaBaW0qgCtKJ9TSo/SMBfzL5VFczomCHJjQ7TqDAvYQIzLg5DoKj67CaUA3oioEASWBFBUP5HzxKAoNQSYXFzUccskv1eHh4sRb99Y/Ds5e0gs7ggrgitIJtXSNgKmA/+fpAyivNPDV0XT2nAxECBjVw4c7wy4xQ/s1Zw0XOOWfgy7gOkLu5fkcwdSgsRQdSyc/qD8R4SMQ1bnm5mx5rM7QrE+lESqKYpbmer7nsgrZdiqLracyScoqBCDYx41KXQHhmgOE+p3hcZcKhuYnscVFz3PBgXSpNNC9woPBvgPYfyGClNLhlEmveq9bN/OkM01uqjRCRVFazBJDFJdvFJOQlE1CUjYHL+ZSoTfg5qxhVO+uBLhrOJ+cQIDnQfBII8ejmEyXqrPRt1zORFsZwdduYZxyDiCzJIpK92HsWzwNMJ1eCOBvZJWpo/TKVRqhoijNao9FKcXllRy6lMue8znsSc7hQk5xozJdnDIY6Hucm4uLiJLn2RySzWYfdwA8DQYGCQ9ivHuQcy6W3UU9SJOBVPXNTbPVlZutoQK4oijNamoPeUttQdxr0WaT1waF+ZKWV0JhWRmxAZeY3PsShYbzHC/KQOhLWZORBcBTQaFkCU/cyoIoKu3NlZI4ruu7Ga1zU0Mr9jIh2qb9wIUQ04F/AU7Ae1LKpQ2u9wA+BPyryyySUm5pc60VRelQTaX4WWrCMMJECqOnqxNllXoKyioBZ84UDsAldxR+7i4kXswBfRkzRRoxmgsU636i0FNLUpcMDCIT2M+MIh2zssI5YejNAXdvLpUNISM/zGQ9HGFCtNkeuBDCCUgGbgXSgJ+A+6SUZ+qUiQeOSSnfEUIMArZIKXs1vJfqgSuKbWvuFCdo+4Rhc8M0OYXlHLmcy0+pefyUmsuJNK3R+2gAN1FApMcxunicY4BBy/ziLMKcshjTsxsGIQiuNBDtGsAgv95M6HkrA/rdAZ5dmmyrLU6ItqUHPhJIkVJerL7RamA2cKZOGQn4Vv+/H5DRtuoqimINxtIFG2rrNqnNpTAG+bgxfUgY04dU9Z5NDbkYgFLpy5mSiVAykX3A+4APhfS6+hNdPFPwDrhBsk7Lt3nH8EhJYMBXT5Dl34PXuvrRx8WXIO/eXC6NJU//c0+9YfsaLi6ypYlScwJ4BHC1zs9pwKgGZZYAO4QQTwJewC3GbpSTk8Pw4Y0+RFi4cGHtYceKonQ8Yysg80p0RssaW4jT0rHkOXERZgc+U0Mu/h4u9A/x4eiVPCoNP48kFOLDubJbcNbdSkmOnnA/d16c5Me0wBuQfY6ctH0kF5/lSqCWqtC2m96VBn6ZFUBxyQCuePXiek5vAoNuMrq4qEZ7DLnEx8cTHx9v7JLR8+7MGUK5B5gmpfxN9c8PACOllE/WKfNM9b3+IYQYQ9UH4RAppaHuvdQQiqLYHlNDGncNi2D9kfRmM1LMzVxp7YShOffPK67g2NU8jl7OZ/vpLM5nF9W7h0bA9CGh3D+qJ9Hd/PBxd2HdD4f57PvP8XJNxuB+jd/nFjDWkMkaH29eDuxCkF4SWu6Gc2lXissiySgcQaYMpWH2S0cMubQ6C6U6IC+RUk6r/nkxgJTy1TplTgPTpZRXq3++CIyWUtY7m0kFcEWxPU2NBT87bUCzQdeczBWgTemJLQn+zY3jCwH9gr2J6x4AwPfnsskuLCfC34PFUyIY7HOafZe2cyYvmeMlOaS7gEEIfki9SoXBi/e9wjnp4kNZWQ+ySgeTVtGfi0tnt6qu5mpLAHemahJzKpBO1STmfCnl6TpltgJrpJSrhBA3Ad8BEbLBzVUAVxTb09a9t5taXANVgdrUkEx79F6bqs9Hj47k2JV8jl7JI/FqPtrqIREvVyeiuvkR2z2A2O7+xPXwJ8TXnXFLE7ihzaa320mG6ooZLFI5EnSBBD+JTlT1xD0NBoYKD94JnsLRiu68mFhOUlkkOtxr29/WfPRWT2JKKSuFEE8A26lKEfxASnlaCPE34LCUciPwJ2ClEOJpqiY0H24YvBVFsU0t2ZmwJc+vUarTm5wUbY9zI03VJ8Lfgwn9g8gtruCLw1fRluoI8XHj1kEhOGkEiVfzeX/fRXT6qtAV6utOsK8b1wp8OVM29uesjUxwySyjp/spwrySiAi+gbdBC8c+YaiuGI/wULxdXehWIfAr84byMDZtH8Wc6D+BhY+mUwt5FKWTa83qy4aTnsUVlbWBryXaowfeVHug6aGcMp2eM5kFfPzDZbadymr0wTucxZEAAAe7SURBVKMR4KwRVOglQd5uLLp9IHcNq15AZNAz9fkPCPI5gKvHZcrc88h0qyDfScP76TmM/PMVcHZtVZvatJBHURTH1dKdCY1lZrhoBAGeLiYzV/w9XCivNJi9fW1bNNWecUsTmtzx0N3FiSs3ShoFb1cnDVMGBlNeqed4mpbc4gpyisr5y1cn+fTgZaK7+RPb3Z9C70guFIRDQc0zDQQ5X+W/nuWMbGXwbooK4IqitCitz9i2rzqDxNPVmRdnDjbaw10ya3Dtczti2bqp9phzEpGx9lXoDZxM17J/0RSklKTllZJ4NZ/Eq/mcSMtnzU9XWXUg1cidNRQQyW3To9rSHJM6NIBrtdrafG+1nayi2KemgmBzvXlrL1Fvary/ZljI1Hh+TbuFEHTv4kn3Lp7MjAkHoFJv4Hx2ESfS8vk6MYOfUnNrh5TKKw28siWJ6UPC8HB1MnpvY+puJ0vVAslG1Bi4oigtYk9L0BtqSc57Qy1tX5lOz9nMAk6kabmaW8L/zhjU6nqrQ40VRbEIez5weE5cBK/OjSLC3wNBVVB+dW4U3yflNBm8jbVvw7F0xi1NIHLRZsYtTWDDsfR6191dnLh8o4T4PRd5f98lo2XaSo2BK4rSIi2d9LQ1xsbHn16TaLK8sW10zdnJsCN2O1QBXFGUFmvJpKc9aCp33NiwiTnnd5pTpq3UEIqiKJ1eS4eFzMlmMadMW6kArihKp2dqbNxUT9nUKtW6j5tTpq1UGqGiKAotGxYytm96wx67OWWaotIIFUVR2ok5uw5aamdCdaixoiiKnVJ54IqiKA7GbgK4iWOG7JKjtMVR2gGqLbbKUdrSXu1QAdwKHKUtjtIOUG2xVY7Slk4fwBVFUZT6OjSA16QRLly4sDY9pk6ajEWYez9Ll2tpWUvez5ptNpejtMWa7xtbb4t6f7Vd3ftt2rSpNl5iIo2wQwO4n58f8fHxxMfH1+aAO8qbsqVlLXk/W39Ttsdr28I/MEu9rqO0Rb2/2q7u/WbOnFkbLwGtsfJqCEVRFMVOdWgeuBAiB7jc4GE/THy6NBAIXDejnLn3s3S5lpS19bZYuh3t8drWaos13ze23pb2+N04Slva2o6eUsqghg92aABXFEVRLEcNoSiKotgpFcAVRVHslE0GcCGEvxBinRAiSQhxVggxRgjRRQixUwhxvvq/AdauZ3OEEAOEEIl1/hQIIf5oj20BEEI8LYQ4LYQ4JYT4XAjhLoSIFEIcrG7LGiGEq7Xr2RwhxFPVbTgthPhj9WN28XcihPhACJEthDhV5zGjdRdV3hZCpAghTgghhlqv5o2ZaMs91X8vBiHE8AblF1e35ZwQYlrH19g0E21ZVh3DTgghvhJC+Ne5ZpG22GQAB/4FbJNSDgRigLPAIuA7KWU/4Lvqn22alPKclDJWShkLDANKgK+ww7YIISKAPwDDpZRDACdgHvAa8FZ1W/KAX1uvls0TQgwBFgAjqXpvzRBC9MN+/k5WAdMbPGaq7rcD/ar/LATe6aA6mmsVjdtyCpgL7Kn7oBBiEFXvt8HVz1kuhDD/iPf2t4rGbdkJDJFSRgPJwGKwbFtsLoALIXyBCcD7AFLKCillPjAb+LC62IfAHOvUsNWmAheklJex37Y4Ax5CCGfAE8gEpgDrqq/bQ1tuAn6UUpZIKSuB3cCd2MnfiZRyD5Db4GFTdZ8NfCSr/Aj4CyHCOqamzTPWFinlWSnlOSPFZwOrpZTlUspLQApVH8I2wURbdlS/x/j/7Z2/axRREMc/A8LBpRANRNBEiH+BpgoWaWJhRBIEi4jggf4RFqnsbUULxUIkRUT0GhULCxsN/oqKP4giJKfxDIIpFCTCWLy33CLZc1kWdgfmA8t7++4dzJdZZnfmvWWBR8Bw7JempXYBHNgHrANXReS5iFwWkQFgl6quAcR2qEojCzALzMe+OS2q+hk4D6wQAvcG8BT4kbpIO0DdP5T4GpgQkUERaQJHgBEM+iRFlu17gNXUPAv+ycK6ltPAndgvTUsdA/g2YAy4qKoHgJ/UN53NRawLTwMLVdtSlFhXnQFGgd3AACFF/5da70tV1beEss994C6wBPzp+ye7yBZjtfZPH8xqEZE5wjV2PRnaYlohLXUM4B2go6qP4/kNQkDvJulfbL9VZF8RpoBnqtqN5xa1HAI+qeq6qm4CN4GDhLQ8+TTfMPClKgPzoqpXVHVMVScIae8yNn2SkGV7h5BdJJjwTwYmtYhICzgKnNTeSzelaaldAFfVr8CqiCQfjpsE3gBtoBXHWsDtCswrygl65ROwqWUFGBeRpogIPb88AI7HOSa0iMhQbPcSFszmsemThCzb28CpuBtlHNhISi0GaQOzItIQkVHCwuxixTb1RUQOA2eBaVX9lfqpPC2qWrsD2A88AV4Ct4AdwCBhhX05tjurtjOnlibwHdieGrOq5RzwjlBHvgY0CGsWi4SFmAWgUbWdOXQ8JNx8loBJSz4h3GzWgE3Ck9yZLNsJqfoF4CPwirCDqHIN/9FyLPZ/A13gXmr+XNTyHpiq2v4cWj4Qat0v4nGpbC3+Kr3jOI5RaldCcRzHcfLhAdxxHMcoHsAdx3GM4gHccRzHKB7AHcdxjOIB3HEcxygewB3HcYziAdxxHMcofwEmsmwZu/fNZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y)\n",
    "plt.plot(x, f(x, *p), label=\"true\")\n",
    "plt.plot(x, f(x, *popt_scipy), '-', label=\"curve_fit\")\n",
    "plt.plot(x, f(x, *popt_jac), '--', label=\"curve_fit with jacobian\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In out easy example, both fits converged well. The Jacobian version is just a bit slower due to the overhead from pytorch and the autograd.\n",
    "\n",
    "It would be interesting to know from more complicated examples, if including the Jacobian can make the difference between convergence and non-convergence! Please let me know if you have data that you can't fit with `curve_fit` without the Jacobian, and tell me if including the Jacobian did help.\n",
    "\n",
    "The open question is of course: did including the Jacobian make actually any difference in the example above? To find out what's happening, to hack around a little bit such that everytime the minimizer calls the numpy version of the function (that is not the Jacobian), the sum of squared residuals will be appended to an array which we will inspect later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_of_squares = []\n",
    "\n",
    "def f(x, mu, sigma, tau):\n",
    "    sig = nt.exp(-(x - mu)**2/(2*sigma**2))\n",
    "    bkg = nt.exp(-(x - mu)/tau)\n",
    "    res = sig + bkg\n",
    "    if isinstance(x, np.ndarray):\n",
    "        sum_of_squares.append(np.sum((res - y)**2))\n",
    "    return res\n",
    "\n",
    "curve_fit(f, x, y, p0=p0)\n",
    "sum_of_squares_nojac = np.array(sum_of_squares)\n",
    "\n",
    "sum_of_squares = []\n",
    "curve_fit(f, x, y, p0=p0,  jac=wrap_jac(f))\n",
    "sum_of_squares_jac = np.array(sum_of_squares)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Function calls without using Jacobian: 39\n",
      "Function calls with Jacobian: 12\n"
     ]
    }
   ],
   "source": [
    "print(\"Function calls without using Jacobian: \" + str(len(sum_of_squares_nojac)))\n",
    "print(\"Function calls with Jacobian: \" + str(len(sum_of_squares_jac)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, there was a difference and with the Jacobian, `curve_fit` needed about 3 times less calls to come to result it was happy with. Let's compare how the residues evolve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEDCAYAAADTIbj3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOX1+PHPyUISyMYSEQIVtLgghB3xi0UQCrig4FpKrUsVfVGXb/0Vi60iaPsV0dZWv35VsK7VqqUUcN8Qqbthd9cKCgEhCGFNQpbz++PeCUMyk9wkc2dJzvv1ymtm7nPnzslA5sx9nueeR1QVY4wxrVNSrAMwxhgTO5YEjDGmFbMkYIwxrZglAWOMacUsCRhjTCtmScAYY1qxlFgH0FidOnXSHj16hGzbtWsXOTk5YZ9bXFxMXl5e2PaGnu93u8Vn8Vl8Fp8f7StWrNiuqqGDU9WE+hk0aJCGc/nll4dtU1Wt77lenu93u8Vn8cWy3eJrufEBhRrmM9W6g4wxphVrUUlgwoQJvj7f7/aGWHwWn5/tDbH4Wmh84U4R4vWnoVOu+jTnudFg8TWPxdc8Fl/zxHN8WHeQMaYhzz77bKxDqJfF54+Emx20a9cupk6dCjinP405BQo8L15ZfM1j8TVPSkp8fxxYfI3z7LPPBiemsNOKRBOsiujgwYO1sLCwUc/5cMkDdF95B4dpMdskj40DpzPkzCsi1r5oVRF3vPw5m0tK6ZqbwfRxxzBxQH7E2o0xpjlEZIWqDg7VFl+pywcfLnmAPituJEMOgMDhFJOz4kY+rK5m0ISprHh2Hn1WzWx8uyqDJlzB4lVF/HbROkorqgEoKinlhoVr0WrlrAH5TWxfB2CJwBjjuxZ/JvDdrB9yOMURj2NTdSdOOnB3xI8bkJ+bwdszTvHt+MbE0iOPPMLYsWPp2rUrAJdddhnXXXcdvXv3btZxN2zYwDvvvMNPf/rTRj3v4osv5owzzuDcc8/1tL2pli1bxp133slzzz1Xp+20007jySefJDc3NyKvFay+M4EWPzB8mIZOAKrw7g+uIFwObKi9q3zPr8YcXe9rN6d9c0lpvc81pjkWrSpi+Jyl9JzxPMPnLGXRqqKovv4jjzzC5s2bax4/+OCDzU4A4CSBJ598stnHiYUXXnjBlwTQkBafBLZJ6Cult0oeJ146l61NbN8mnbh2TC/yczNCtufnZjSrvWuY7cY016JVRdywcB1FJaUoB7sgm5sI/va3vzF06FD69+/PFVdcQVVVFVVVVVx88cX06dOHvn37ctddd7FgwQIKCwuZMmUK/fv3p7S0lJEjRxI4w8/MzOQ3v/kNgwYNYsyYMXzwwQeMHDmSI488kiVLlgDOh/2PfvQjBg4cyMCBA3nnnXcAmDFjBv/+97/p378/d911F1VVVUyfPp0hQ4ZQUFDAAw88ADhT46+66ip69+7N6aefzrZt2xr8/W655RaGDBlCnz59mDp1KoFelK+++ooxY8bQr18/Bg4cyH/+8x9UlenTp9f83k8//XTNcXbv3s2kSZPo3bs3V155JdXVTldwjx492L59OwATJ05k0KBBHH/88cybN6/muZmZmfzud7+jX79+DBs2jK1btzbr3wxawZjAxoHTyQmMCbhKtQ0bB03n8Ai0Tx93DDcsXEdpRVVNe0ZqMtPHHQMRaDemsWY/+zGfbN4dtn3VtyUcqKo+ZFtpRRXXL1jL3z/4NuRzenfN5uYJx4c95qeffsrTTz/N22+/TWpqKtOmTeOJJ57g+OOPp6ioiI8++giAkpIScnNz+d///V/uvPNOBg+u20Oxb98+Ro4cye23386kSZO48cYbefXVV/nkk0+46KKLOPPMMznssMN49dVXSU9P58svv2Ty5MkUFhYyZ86cQ7pb5s2bR05ODh9++CHl5eUMHz6csWPHsmrVKj7//HPWrVvH1q1b6d27N5deemm97+tVV13FzJkzAbjwwgt57rnnmDBhAlOmTGHGjBlMmjSJsrIyqqurWbhwIatXr2bNmjVs376dIUOGMGLECAA++OADPvnkE4444gjGjx/PwoUL63Q3PfTQQ3To0IHS0lKGDBnCOeecQ8eOHdm3bx/Dhg3jD3/4A9dffz3z58/nxhtvrDfuhrT4JDDkzCv4ENzZPdvZJp3YOOjg7B6v7QUrf0cbrWCr5B3SHhi8DTe7x2v7b/65lvLKavJtdpDxWe0E0NB2L15//XVWrFjBkCFDACgtLeWwww5jwoQJfP3111x99dWcfvrpjB07tsFjtWnThvHjxwPQt29f0tLSSE1NpW/fvmzYsAGAiooKrrrqKlavXk1ycjJffPFFyGO98sorrF27lgULFgDOFPMvv/yS5cuXM3nyZJKTk+natSunnNLw+Nsbb7zB3Llz2b9/Pzt27OD4449n5MiRFBUVMWnSJADS09MBeOutt2qO37lzZ04++WQ+/PBDsrOzGTp0KEceeSQAkydP5q233qqTBO6++27+9a9/AbBx40a+/PJLOnbsSJs2bTjjjDMAGDRoEK+++mqDcTekxScBcD7IcT+0D3d/Gt2+fznsWM/h096p0z5xQH69H9pe2l/7dCufbN7N0l+P9PIrGRNWfd/YAYbPWUpRiDGn/NwMnr7ixCa9pqpy0UUXcdttt9VpW7NmDS+//DL33nsvzzzzDA899FC9x0pNTUVEAEhKSiItLa3mfmVlJQB33XUXnTt3Zs2aNVRXV9d8+IaK65577mHcuHGHbH/hhRdqXsOLsrIypk2bRmFhId27d2fWrFmUlZURbmJNfRNuar9u7cfLli3jtdde491336Vt27aMHDmSsrIy4ND3Jjk5ueb9aI4WPyYQMWlZUL7Ht8Nnpaewp7z5/6DGNGT6uGPISE0+ZFtzuyBHjx7NggULavrWd+zYwTfffMP27duprq7mnHPO4dZbb2XlypUAZGVlsWdP0/+edu3aRZcuXUhKSuLxxx+nqqoq5HHHjRvHfffdR0VFBQBffPEF+/btY8SIETz11FNUVVWxZcsW3njjjXpfL/Ah3KlTJ/bu3VtzZpGdnU23bt1YtGgRAOXl5ezfv58RI0bw9NNPU1VVRXFxMcuXL2fo0KGA0x20fv16qqurefrppznppJPq/G7t27enbdu2fPbZZ7z33ntNfp+8aBVnAhGRlgXl4ftZmyszLYW9ZZYEjP8a6qJsit69e/P73/+esWPHUl1dTWpqKvfeey8ZGRlccsklNYOfgTOFiy++mCuvvJKMjAzefffdRr/etGnTOOecc/jHP/7BqFGjaNeuHQAFBQWkpKTQr18/Lr74Yq699lo2bNjAwIEDUVXy8vJYtGgRkyZNYunSpfTt25ejjz6ak08+OeTrVFZWkpaWRm5uLpdffjl9+/alR48eNd1eAI8//jhXXHEFM2fOJDU1lX/84x9MmjSJd999l379+iEizJ07l8MPP5zPPvuME088kRkzZrBu3TpGjBhR05UUMH78eO6//34KCgo45phjGDZsWKPfn8ZIuOsEevXqpaNGjQIaXzaiWV6bDe/cDTdth0acRnr159e+4M+vfcl//uc0kpMif3xjTONUV1czZMgQHnvsMY4/vv4utngUXDZi/vz5X6lqr1D7JdyZQE5OziFTpqImLQuqK6GyDFIjP30zM835p9hbXklORmrEj2+M8W7z5s2MGTOGUaNGJWQCgEO/JM+fP39XuP0SLgnETFqWc1u+x5ckkJVuScCYeNG1a1c++eSTWIcRFTYw7FVatnPr0+BwZprzwW/jAsaYaLIk4FXNmYA/g8OZNWcCFb4c3xhjQrEk4FVwd5APAmMCe+xMwBgTRZYEvErLdG7L9/py+EAS2Fde1cCexhgTOZYEvPJ7TMC6g4wBnAqjV111VaOft2TJEubMmeN5/+CCbZEwa9Ys7rzzzjrbN2/eHLFS1H6w2UFe+T0mYN1BJprWPgOv3wK7NkFONxg9EwrOj3VUzXLmmWdy5plnxjqMOrp27VpzhXE8sjMBr6I0JrDXSkcYv619Bp69BnZtBNS5ffYaZ3szNLb88bPPPssJJ5zAgAEDGDNmTJ2yyHv27KFnz541JR92795Njx49qKio4O6776Z3794UFBTwk5/8BDj0DGLr1q1MmjSJfv360a9fv5pS042N/aWXXmLgwIH069eP0aNHA05JjIkTJ1JQUMCwYcNYu3Ztzf5r1qzhlFNOoVevXsyfPx9wyl736dOn5n6oEtjLli1j5MiRnHvuuRx77LFMmTKl3vpDkWRnAl6lpENSim9JIDlJaNsm2aaImuZ7cQZ8ty58+6YPoar80G0VpbD4KljxaOjnHN4XTq2/q6Wx5Y9POukk3nvvPUSEBx98kLlz5/LHP/6x5nhZWVmMHDmS559/nokTJ/LUU09xzjnnkJqaypw5c1i/fj1paWmUlJTUieWaa67h5JNP5l//+hdVVVXs3Vv/WF6o2Kurq7n88stZvnw5PXv2ZMeOHQDcfPPNDBgwgEWLFrF06VJ+/vOfs3r1agDWrl3Le++9x759+xgwYACnn376Ia8TrgQ2wKpVq/j444/p2rUrw4cP5+23365TV8gPCZcEdu3axdSpU4Eol40Q8b2IXGZaip0JGP/VTgANbfeoseWPN23axAUXXMCWLVs4cOAAPXv2rHPMyy67jLlz5zJx4kQefvjhmm/XBQUFTJkyhYkTJzJx4sQ6z1u6dCmPPfYY4FTbzMnJaXTsxcXFjBgxoiauDh06AE6Z6H/+858AnHLKKXz//ffs2uVckHvWWWeRkZFBRkYGo0aN4oMPPqB///41r1NfCeyhQ4fSrVs3APr378+GDRualQSCy0YAYd+AhEsCMSsbAdAmCw74MzsInCRglURNszXwjZ27+rhdQbXkdIdLnm/SSzal/PHVV1/Nddddx5lnnsmyZcuYNWtWneMOHz6cDRs28Oabb1JVVVXTrfL888+zfPlylixZwq233srHH3/cpLjri11VQ5abDtVNE9ivoTLR9ZXADpTMhsiUifZaNsLGBBrD7zOB9BT2WRIwfhs9s27pk9QMZ3sTNaX88a5du8jPdyqXPvpomG4o4Oc//zmTJ0/mkksuAZzCbhs3bmTUqFHMnTuXkpKSOt09o0eP5r777gOgqqqK3bvDT+gIF/uJJ57Im2++yfr16wFquoNGjBjBE088ATgJpFOnTmRnO7MHFy9eTFlZGd9//z3Lli07pNpo4LVClcCOJUsCjWHlpE1LUHA+TLjb+eaPOLcT7m7W7KDx48dTWVlJQUEBN910k6fyx7NmzeK8887jRz/6EZ06dQq735QpU9i5cyeTJ08GnA/1n/3sZ/Tt25cBAwbwq1/9qs4C7X/5y19444036Nu3L4MGDQp5phAoEx0u9ry8PObNm8fZZ59Nv379uOCCC2riLiwspKCggBkzZhySwIYOHcrpp5/OsGHDuOmmm+jateshrzlt2jQeffRRhg0bxhdffFFTAjuWEq6U9ODBgzUwkBJ1fzsX9m+Hqct8OfzUxwr5dsd+XvrvEb4c35hEtGDBAhYvXszjjz8esWMWFxfTv39/ioqKInbMeCYiK1S17oLOJOCYQEylZcHO9b4dPjM9xa4TMCbI1VdfzYsvvsgLL7wQsWMuWbKE66+/PuRSmK2Rr0lARH4FXAYosA64BOgCPAV0AFYCF6rqARFJAx4DBgHfAxeo6gY/42s0v5eYtNlBxhzinnvuifgx4/WisljxbUxARPKBa4DBqtoHSAZ+AtwO3OWucrMT+IX7lF8AO1X1h8Bd7n7xJS3Lt9pBAO3cJJBoXXTGmMTl98BwCpAhIilAW2ALcAoQuIb6USAwyfcs9zFu+2gJNT8rltKyoGIfVPszop+ZnkJVtVJeWe3L8Y0xpjbfkoCqFgF3At/ifPjvAlYAJaoa6PPYBARWt84HNrrPrXT37+hXfE3ic+mILKsfZIyJMt/GBESkPc63+55ACfAP4NQQuwb6PkJ966/TL1JcXMzgwXUHuadOnVpzJbFvgpNARm79+zZBZtASk3lZaQ3sbYwxh9qwYQOXX345O3furN0Udg6unwPDY4D1qloMICILgf8CckUkxf223w3Y7O6/CegObHK7j3KAHbUPmpeXR8ymiPpeRM6WmDTGNF2PHj1qynIEE5GwNbP9HBP4FhgmIm3dvv3RwCfAG0CguPZFwGL3/hL3MW77Uo23EdJorS5mawoYY6LEzzGB93EGeFfiTA9NAuYBvwGuE5GvcPr8/+o+5a9AR3f7dcAMv2JrsjZuEjjgczlpOxMwxkSJr9cJqOrNwM21Nn8NDA2xbxlwnp/xNJvfZwLumMC+A5YEjDHRYbWDGiNaC8vYmYAxJkosCTSG31NE0wNjApYEjDHRYUmgMXxOAmkpSaQkiZ0JGGOixpJAYyQlQ2o735KAiJCZbvWDjDHR02ASEJFrRSRbHH8VkZUiMjYawcWltExfi8i1a2NrChhjosfLmcClqrobGAvk4VQCbWD9uhbM70qidiZgjIkiL1NEA+UcTgMeVtU1sSzsFrOF5gNssXljTAKI5ELzK0TkFZwaQDeISBYQszKXMV1oHqKyzvCOfQd8O74xpnXwutC8lyTwC6A/8LWq7heRjjhdQq1TWjbs83F1sbQUvv1+v2/HN8aYYGGTgIgMrLXpyHgr7x8TaVm+lY0AZ0zArhMwxkRLfWcCf6ynTXEWh2l92tjsIGNMyxE2CajqqGgGkjACYwKq4MOZUWZ6CqUVVVRVK8lJduZljPGXpwJyItIH6A2kB7ap6mN+BRXX0rKguhIqyyA1I+KHr6kfVF5JTkZqxI9vjDHBGkwCInIzMBInCbyAszrYW0DrTQLgnA34kASy0i0JGGOix8vFYufiLAjznapeAvQDWu/ah2nZzq2tLmaMaQG8JIFSVa0GKkUkG9gGHOlvWHEsSmsK7LXVxYwxUeBlTKBQRHKB+cAKYC/wga9R1SMurhgG/5eYtDMBY0wzROyKYVWd5t69X0ReArJVdW3zQ2ya2F8xnOnc+pwE9pVX+XJ8Y0zrELErhkVkRKhtqrq8WREmKr/HBKw7yBgTRV66g6YH3U/HWR94Ba31YrGa7qDdvhzeuoOMMdHkpTvokE53EekOzPUtongXrXWGrXSEMSYKmrKy2CagT6QDSRgp6ZCUAgf2+nL45CShbZtkmyJqjIkKL2MC9+DUCgInafQH1vgZVFwTsTUFjDEthqcpokH3K4G/q+rbPsWTGNpYEjDGtAxexgQejUYgCSUKC8tYEjDGREN96wms42A3UB2qWuBLRIkgLcu32UHgngnYmIAxJgrqOxM4w739pXv7uHs7BWjdS1+lZcH+7b4dPjMthW/3te632BgTHfWtJ/ANgIgMV9XhQU0zRORt4Ba/gwsl5mUjwEkCOzf4dvjM9BS7TsAY0yyRXGi+nYicpKpvAYjIfwHtmh9i08S8bAT4PiaQZQPDxphmivRC8w+JSCCTlACXNjvCROZzEmiXlsK+8kpUFVvX2RjjJy+zg1YA/dwy0qKqYTNKq5GWBRX7oLoKkpIjfvjM9BQqq5XyymrSUyN/fGOMCahvdtDPVPVvInJdre0AqOqffI4tfgWXjsjIjfjhs4LqB1kSMMb4qb6yEYF+/6wwP61X1BaWsXEBY4y/6psd9IB7Ozt64SSIQBLwqX6QLTFpjImWBgvIichcEckWkVQReV1EtovIz7wcXERyRWSBiHwmIp+KyIki0kFEXhWRL93b9u6+IiJ3i8hXIrJWRAY295fzTbRWF7M1BYwxPvNSRXSsqu7GuXhsE3A0h64xUJ+/AC+p6rE4C9R/CswAXlfVXsDr7mOAU4Fe7s9U4D6vv0TUtYnOmgJ2JmCM8ZuXJJDq3p6GUzxuh5cDu7OJRgB/BVDVA6paApwFBOoRPQpMdO+fBTymjveAXBHp4u3XiLIojQnsO2BJwBjjLy/XCTwrIp8BpcA0EckDyjw870igGHhYRPrhrEZ2LdBZVbcAqOoWETnM3T8f2Bj0/E3uti3BBy0uLmbw4MF1Xmzq1Kk1VxL7LloLy9iZgDGmEebNmxfuYtpO4Z7j5TqBGSJyO7BbVatEZD/Ot/aGpAADgatV9X0R+QsHu35CCXVVVJ0Cdnl5eRQWFobYNYp8TgJZ6YExAUsCxhjvwn0ZFpGwxc68DAy3xSkiF+ij7wrU/Spe1yZgk6q+7z5egJMUtga6edzbbUH7dw96fjdgs4fXib6aJODP7KC0lCRSksTOBIwxvvMyJvAwcAD4L/fxJuD3DT1JVb8DNorIMe6m0cAnwBLgInfbRcBi9/4S4OfuLKFhwK5At1HcSUqG1Ha+DQyLiK0pYIyJCi9jAkep6gUiMhlAVUvFe0Gbq4EnRKQN8DVwCU7ieUZEfgF8C5zn7vsCzuDzVzilqi/x/mvEQFqmv/WD2tiaAsYY/3lJAgdEJAO3f15EjgLKvRxcVVcTuutodIh9lYNrF8Q/vyuJ2pmAMSYKvCSBm4GXgO4i8gQwHLjYz6ASgi02b4xpAepNAm63z2fA2cAwnBk816qqf8tqJYoorDO8Y98B345vjDHQQBJQVRWRRao6CHg+SjElhrRs2L/Bt8NnpqXw7fe2xKQxxl9eZge9JyJDfI8k0fi82HxWeopdJ2CM8Z2XMYFRwBUi8g2wD6dLSFW1wNfI4l0bmx1kjEl8XpLAqb5H0QhxsdA8HBwTUAUfloDMTE+htKKKqmolOcmWmDTGNE7EFppX1W8iFVQkxMVC8+AkgepKqCyD1IyIH76mflB5JTkZqQ3sbYwxh/K60LyXMQETSpTqB9k0UWOMnywJNFVatnPrWyVRW13MGOM/SwJNFbV1hm11MWOMf8KOCYjIHkKUcg5Q1WxfIkoUaZnOrW9nAskA7LEzAWOMj+pbaD4LQERuAb4DHseZHjoFyIpKdPHM94VlnO6gfeVVvhzfGGPAW3fQOFX9P1Xdo6q7VfU+4By/A4t7fo8JWHeQMSYKvCSBKhGZIiLJIpIkIlMA+3qaFp3F5q07yBjjJy9J4KfA+cBW9+c8d1vrFkgCB/xZXSz4OgFjjPGLl4vFNuBtTeHWJSUdklJ86w5KThLatkm2KaLGGF81mARE5Gic9YU7q2ofESkAzlTVBpeY9EPclI0Q8b9+kK0pYIxpooiVjQDmA9OBBwBUda2IPImHdYb9EDdlI8AZHPZzdTFLAsaYJopk2Yi2qvpBrW32yQRRWVjGkoAxxk9eksB2d13hwBrD5wJbfI0qUURjiUkbEzDG+MhLd9AvgXnAsSJSBKzHuWDMpGXB/u99O3xmWgrf7rPVxYwx/mlojeEkYLCqjhGRdkCSqvr31TfRpGXBzg2+HT4zPcWuEzDG+Kre7iBVrQaucu/vswRQi8/dQTYwbIzxm5cxgVdF5Nci0l1EOgR+fI8sEficBNqlpbCvvBLVsHX8jDGmWbyMCVzq3v4yaJsCR0Y+nASTlgUV+6C6CpKSI374zPQUKquV8spq0lMjf3xjjPFyxXDPaASSkIJLR6SHvRajybKC6gdZEjDG+MHLmQAi0gfoDaQHtqnqY34FVZ+4uWIYDi0n7UMSyAxaYjIvKy3ixzfGtFwRu2JYRG4GRuIkgReAU4G3gJgkgfi6Yjg6awrYtQLGmMaK5BXD5wKjge9U9RKgH2BfSyEKScDtDrI1BYwxPvGSBErdqaKVIpINbMMGhR1torOmgK0uZozxi5cxgUIRycUpJLcC2AvUriXUOtli88aYBOdldtA09+79IvISkK2qa/0NK0HUJAGfF5axMQFjjE+8DAyPCLVNVZf7E1IC8flMICs9MCZgScAY4w8v3UHTg+6nA0NxuoVO8fICIpIMFAJFqnqGiPQEngI6ACuBC1X1gIik4cw4GgR8D1zgrmoWv3xOAmkpSaQkiZ0JGGN80+DAsKpOCPr5MdAHZ61hr64FPg16fDtwl6r2AnYCv3C3/wLYqao/BO5y94tvScmQ2s63gWERsTUFjDG+8jI7qLZNOImgQSLSDTgdeNB9LDhnEAvcXR4FJrr3z3If47aPdvePb2k+LzHZxtYUMMb4x8uYwD24C8rgJI3+wBqPx/8zcD3g9pvQEShR1cCn2iYg372fD2wEUNVKEdnl7r/d42vFht+VRO1MwBjjI09TRIPuVwJ/V9W3G3qSiJwBbFPVFSIyMrA5xK7qoa1GcXExgwcPrrPj1KlTa8pJRFVallM7yCeZVk7aGOPRvHnzwlVU6BTuOV6miD7a0D5hDAfOFJHTcAaUs3HODHJFJMU9G+gGbHb33wR0BzaJSApOrYsdtQ+al5dHYWFh7c2xE4V1hnfsO+Db8Y0xLUe4L8MiErZHpcExARFZJyJrQ/ysE5Gw1wuo6g2q2k1VewA/AZaq6hTgDZxSFAAXAYvd+0vcx7jtSzURCumnZds6w8aYhOWlO+hF9/Zx93YKsJ+Dg7iN9RvgKRH5PbAK+Ku7/a/A4yLyFc4ZwE+aePzoSsvybXYQOGMCdp2AMcYvXpLAcFUdHvR4hoi8raq3eH0RVV0GLHPvf41zrUHtfcqA87weM260sdlBxpjE5WWKaDsROSnwQET+C2jnX0gJJjAm4FPPVWZ6CqUVVVRVx3/PmDEm8Xg5E/gF8JCIBBYlKOHgkpMmLQuqK6GyHFLTG96/kWrqB5VXkpORGvHjG2NaNy+zg1YA/dwy0qKqYRcnaJWCS0f4kASy0i0JGGP842V20LVuAtgD/FFEVorIWP9DSxBp2c6tb2sK2Opixhj/eBkTuFRVdwNjgcOAS4A5vkaVSGxNAWNMAvMyJhC4kvc04GFVXRPLmj5xtdA8OLWDwMclJpMB2GNnAsaYRojYQvPAChF5BegJ3CAiWUB180NsmrhaaB6itti8LTFpjGkMrwvNe50d1B/4WlX3i0hHnC4hAwfHBHyqH2TdQcYYP3mZHVSNs/hL4PH3OIu+GAg6E/B3sXnrDjLG+KEp6wmYYL53Bx2cImqMMZEWNgm4y0CahqSkQ1KKb0kgOUlo2ybZpogaY3xR35nAAgAReT1KsSQmEf/rB9maAsYYn9Q3JpAkIjcDR4vIdbUbVfVP/oWVYHwuJ51lScAY45P6zgR+ApRaDco5AAAVwklEQVThJIqsED8mIAoLy1gSMMb4IeyZgKp+DtwuImtV9cVw+xn8TwK2sIwxxideZge9IyJ/EpFC9+ePQRVFDUQnCdiZgDHGB14uFnsI+Ag43318IfAwcLZfQdUn7spGgJMEdm7w7fCZ6Sl2nYAxplEiWTbiKFU9J+jxbBFZ3ZzgmiPuykaAUz/IzgSMMXHEa9kIL91BpbVWFhsOlDY7wpYkCovN7yuvRH1avcwY03p5ORO4EngsaBxgJ3CRfyEloLQsqNgH1VWQlBzxw2emp1BZrZRXVpOeGvnjG2NaLy+1g9ZwcGUx3LUFTLBA6YgDeyE98mPmWUH1gywJGGMiyXPtIFXdbQkgjKgtLGPjAsaYyLICcpEQpTUF7FoBY0ykWRKIhDb+JoF2gdXFbE0BY0yENTgmICLJwOlAj+D9rXZQEJ/XFMiy1cWMMT7xMjvoWZwaQuuI4bKSca0mCdjqYsaYxOIlCXRT1QLfI0lk0VpYxsYEjDER5iUJvCgiY1X1Fd+j8SAuy0Z85S65sOQqePN2GD0TCs6v/zmNsOzzbQDctPhj7n/za6aPO4aJA/Jr2hetKuKOlz9nc0kpXXMzIt5ujEk8XstGSENXoYrIJOBvOIPIFYAAqqrZkQm1cQYPHqyFhYWxeOnQ1j4Dz14DFUEXUadmwIS7I5IIFq0q4oaFaymtONgTl5GazG1n92XigHy3fR2lFVW+tBtjEp+IrFDVwSHbPCSBr4GJwDqNg7oFcZcE7uoDuzbW3Z7THX71UbMPP3zOUopK6lbpSE0WjuuSzadbdlNRVfefpbnt+bkZvD3jlGbHb4yJvfqSgJcpol8CH8VDAohLuzY1bnsjbQ6RAAAqqpSO7dqE/ACPRHu41zXGtCxeksAWYJmI3CAi1wV+/A4sYeR0a9z2RuqamxFye35uBg9fMpR8n9rDva4xpmXxkgTWA68DbbDlJesaPdMZAwiWmuFsj4Dp444ho1a9oIzUZKaPOyYq7caYls1LAbnZ0QgkYQUGf1+9CfZ8Bxnt4dS5EZsdFBicDTd7J1Ltc178jO92l5GdnsItZ/WxQWFjWgkvA8NvAHV2UtV6Rw1FpDvwGHA4zkVm81T1LyLSAXga5wrkDcD5qrpTRAT4C3AasB+4WFVX1j5u3A0MB6jCHUfB0afCxHtjHU2TnHzHGxx3eDb3Xzgo1qEYYyKovoFhL9cJ/DrofjpwDuDlqqVK4P+p6koRyQJWiMirwMXA66o6R0RmADOA3wCnAr3cnxOA+9zbxCAC3U+Aje/HOpImG/SD9vz7q+2oKk5ONsa0dA2OCajqiqCft1X1Ojx8OKvqlsA3eVXdA3wK5ANnAY+6uz2KM/0Ud/tj6ngPyBWRLo3/lWKo+1D4/kvY932sI2mSgUe0p3hPOZt22swgY1oLLwXkOgQ9TAIG4XTxeCYiPYABwPtAZ1XdAk6iEJHD3N3ygeAJ95vcbVuCj1VcXMzgwXXPaqZOnVpzJXHMdHdz46YP4ZjxsY2lCQYd0R6AFd/spHuHtjGOxhjTWPPmzQu3BnuncM/x0h20AmdMQHC6eNYDv/AalIhkAv8E/ltVd9fTzRCqoc5YRF5eHnE5JgDQdQAkpThdQgmYBI7unEVmWgorvtlpA8PGJKBwX4ZFZHu453iZHdSzqQGJSCpOAnhCVRe6m7eKSBf3LKALsM3dvgnoHvT0bsDmpr52TKRmQJd+sPGDWEfSJMlJwoAf5LLim52xDsUYEyUNjgmIyHnuwC4icqOILBSRgR6eJ8BfgU9rrT2whIML1V8ELA7a/nNxDAN2BbqNEkr3E6BoBVQlZtnngT9oz2ff7balLI1pJbxcLHaTqu4RkZOAcTiDufd5eN5w4ELgFBFZ7f6cBswBfiwiXwI/dh8DvAB8DXwFzAemNe5XiRPdh0JlKXy3NtaRNMmgI9pTrbBmY0msQzHGRIGXMYFAecnTgftUdbGIzGroSar6FqH7+QFGh9hfgV96iCe+dRvq3G78APITb759/x/kIuIMDg//YdixJGNMC+HlTKBIRB4AzgdeEJE0j89rnXLynQqiCXq9QHZ6Ksd0zrJxAWNaCS8f5ucDLwPjVbUE6ABM9zWqRNd9aMIODoNzvcDKb3dSXW2FY41p6bxcLLZfVReq6pfu4y3xsspY3Op+Auwuilg56Wgb9IP27Cmr5Ktif9ZMNsbED+vW8UP3wLhAYnYJBV80Zoxp2SwJ+KFzH0htm7BdQkd0bEvHdm0o3GBJwJiWzsvsoLgSlwvN15ac6swMStAzARGpGRcwxiSmiC00H2/itpR0ba/fAm/9GW7YCG3axTqaRrv/zf8w58XPWHHjGDpmpsU6HGNMMzR3jWHTFN1PAK2CzatiHUmTBMYFVn5rF40Z05JZEvBLtyHObYJ2CfXNzyE1WWxw2JgWzpKAX9p2gE7HJOzgcHpqMsd3zWGlJQFjWjRLAn7qPtQ5E0iwcZeAQUe0Z82mEg5UVsc6FGOMTywJ+Kn7CVC6E77/KtaRNMmgI9pTXlnNJ1t2xzoUY4xPLAn4KbDSWIKOC9hFY8a0fJYE/NTxh5DRPmGTQOfsdPJzM2xcwJgWzJKAn5KSnNLSCTo4DM7ZQOE3O0i060mMMd7YFcN+6z4UvnzZGRvIaB/raBpt0BHtWbJmM5t3lZGfmxHrcIwxHnm9YjjhkkBOTg7z5s2LdRjeBcYFNhVCrx/HNpYmCB4XsCRgTOII/pI8f/78XeH2s+4gv+UPBEmGb9+LdSRNcuzhWWSkJtu4gDEtlCUBv7VpB4f3TdjB4ZTkJPp3z7UZQsa0UAnXHZSQMjrA12/ArFzI6QajZ0LB+bGOyrPMtGTe/fp7es54nq65GUwfdwwTB+TXtC9aVcQdL3/O5pJSX9qNMf6xJOC3tc/AN28B7uyaXRvh2Wuc+wmQCBatKuLNL7YDzm9QVFLKDQvXATBxQD6LVhVxw8J1lFZUgQ/txhh/WRLw2+u3QNWBQ7dVlDrbEyAJ3PHy5xyoOrRsRGlFFTMXf8SWXWXct+yrmg/wSLbf8fLnlgSMiQJLAn4Lt85wgqw/vLmkNOT23WWV3P7SZ2Gf19z2cK9rjIksGxj2W063xm2PM13DTAvtmpPOZ7eOp2tOuj/tNh3VmKiwJOC30TMhtdYHWlKqsz0BTB93DBmpyYdsy0hN5vrxx5Lu3ka6PSVJmD7uGH9+IWPMIaw7yG+Bfv/Xb3G6gFLSoLIcOhwZ27g8CvTLh5u9E+n2jDbJ7D9QRW7b1Kj+nsa0Vgm3xnCvXr101KhRQIKUjaitdCfcPwJQuGK5s/iMqVFWUcXEe99m6+4yXrj2R3TJsW4hY5oiuGzE/Pnzv1LVXqH2S7gkkDALzddn0wp4aBz8cAxM/juIxDqiuPKf4r2cec9bHNclm79PHUZqsvVaGtMcttB8vOk2CMbeCl+8CO/+b6yjiTtH5WXyP2f3pfCbndz5yuexDseYFs2SQKyccCUcNwFem5XQpab9clb/fH56wg944M2vef3TrbEOx5gWy7qDYqm0BB4YAWW7oU1b2L05IctK+KWsooqz/+8d1m/fS3ZGKtt2l1vZCmOawLqD4lVGLgz4GZTthN1FgB4sK7H2mVhHF3PpqclMGpBPaUU1W3eXH1K2YtGqIoCashNFJaW+tBvT0tmZQKzd1cf54K8tpzv86qPoxxNnhs9ZSlGIq4dTkoQj89rxdfE+Kqvr/h9ubnt+bgZvzzglMr+EMTFW35lAXF0nICLjgb8AycCDqjonxiH5L2xZiY2w7TM47FjnrCBwnUGo7qJ4b2+GcOUjKquVo/Iy+WLrXl/ai0pKWbOxhIJuOSxevdnX7iZrt/ZYdlcmz5o1K2IHaw4RSQZeAsYBtwF3z549e/msWbOKg/ebN2/erMDyki3Cqr9B+e7QbR8+CKuehLVPO11G4Oz71WuQ+wPofLzzAfzsNbD/+/hsb6Z/FG5iT1llne35uRksnDbct3aApz7cyN8/+JbFqzezq7QCgD1llbz5RTHd2mdwbJfsmu6kHfsPWLu1R73dq9mzZ2+ZNWtWyCUZ46Y7SEROBGap6jj38Q0Aqnpb8H4trjso8CFaEfSNNzUDfvwHqCqHV2dCdUXd56VmQN/z4eOFUL6nbntGBxh/G7x0A5TuqNvetiOcdie88OuDH+C120//Izz//5rWHqHurNqlpsEpO3Hb2X1DlqKOVPtNE45DEG5e/HGdKqoAySIcnpPOd7vKqArxN9Tc9vTUJM7ql8/z6zazt7yqTntuRiozTj2WOS9+Rklp3f8fzW5vm8rMM3pz63OfsHN/3fb2bVO5dWIfblr0Uczafz+xLzcuWtcq2xvbXVlfd1A8JYFzgfGqepn7+ELgBFW9Kni/FpcEoP7ulFm51KxFUFtmZ9gbr9MnBWaVRORIsTzd7jnj+XDvPucO6saCFeGrwTa3vXN2Glt3l4dtN62XAOvnnO59/wSZHRTqstk6f3/FxcUMHjy4zk9CLT5fW8H5zrfmWSXObXB/etgqpN3h1184t6FkdYGrV0LW4aHbMw+Hae87t2Hb32t6ewSrpE4ckM/bM05h/ZzTeXvGKXX6Q/1sD1fNND83gzvP60e+j+3v/3ZM2PbDs9N494ZTODw7zZf2ztlpvDl9JJ3DtB+WlcarvxrBYVmxa3+lFbeH+385b968kJ+PQKeQTyC+ksAmIPgTrRuwufZOeXl5FBYW1vlpUeMEwUJVIU3NOFiFNFz7j2+BjkfBj28N3T72VmfQeWx97cc1vT1BqqQ2JFwV1UCV01i1zzj1OLrkZDDj1ON8ab/h1OM4omM7bgjT/tvTjqNX5yx+e1rs2o9uxe3hquxOnTo15OcjsD3kE4ijgeHZs2d/B8yaPXv2ktmzZ+8H7gb+p8UPDDek8/HOIOvm1U7ff053GD/n4NlCvLcnuGO7ZNOtfQbrinaxt6yS/NwMZk7oXXO2YO3WHst2r+obGEZV4+YHOA34AvgP8LtQ+wwaNEib6oEHHmjyc6PB4msei695LL7mief4gEIN87kbT91BqOoLqnq0qh6lqn+I9PHjfdzA4msei695LL7miff4womrJGCMMSa6WlQSCCyg4Nfz/W5viMVn8fnZ3hCLr2XGZ0mgEc+3/0T1s/gsPj/bG2LxWRIwxhjTSHFzxbBXIlIMfBOmOQfYVc/TO1HPfFkPz/e73eKz+Cy+8Cy+prcfoap5oRoSLgkYY4yJHOsOMsaYVsySgDHGtGKtIgmIyHgR+VxEvhKRGbGOpzYR2SAi60RktYjERYlUEXlIRLaJyEdB2zqIyKsi8qV72z7O4pslIkXu+7haRE6LUWzdReQNEflURD4WkWvd7XHx/tUTX1y8f24s6SLygYiscWOc7W7vKSLvu+/h0yLSJs7ie0RE1ge9h/1jEV9jtPgxAXexmi+AH+MUqfsQmKyqn8Q0sCAisgEYrKr1DSpFlYiMAPYCj6lqH3fbXGCHqs5xk2l7Vf1NHMU3C9irqnfGIqag2LoAXVR1pYhkASuAicDFxMH7V0985xMH7x+AiAjQTlX3ikgq8BZwLXAdsFBVnxKR+4E1qnpfHMV3JfCcqi6IdkxN1RrOBIYCX6nq16p6AHgKOCvGMcU9VV0O1F6N5izgUff+ozgfHDERJr64oKpbVHWle38P8CmQT5y8f/XEFzfckjeBtT9T3R8FTgECH7CxfA/DxZdwWkMSyAeCV3LfRJz9h8f5z/OKiKwQkXgukdpZVbeA80ECHBbjeEK5SkTWut1FMeuuChCRHsAA4H3i8P2rFR/E0fsnIskishrYBryKU1iyRFUD64HG9G+5dnyqGngP/+C+h3eJSOgFAeJIa0gCnharibHhqjoQOBX4pdvVYRrvPuAooD+wBfhjLIMRkUzgn8B/q2qYhaRjJ0R8cfX+qWqVqvbHWVtkKHBcqN2iG1XQC9eKT0T6ADcAxwJDgA5ATLpLG6M1JAFPi9XEkqpudm+3Af/C+Q8fj7a6/cmBfuVtMY7nEKq61f3DrAbmE8P30e0n/ifwhKoudDfHzfsXKr54ev+CqWoJsAwYBuSKSIrbFBd/y0HxjXe72lRVy4GHiZP3sD6tIQl8CPRyZxW0AX4CLIlxTDVEpJ07OIeItAPGAs1fod0fS4CL3PsXAYtjGEsdgQ9Y1yRi9D66g4Z/BT5V1T8FNcXF+xcuvnh5/9xY8kQk172fAYzBGbt4AzjX3S2W72Go+D4LSvKCM14Rr3/LNVr87CAAd6rbn4Fk4CE/1ipoKhE5EufbP0AK8GQ8xCcifwdG4lwKvxW4GVgEPAP8APgWOE9VYzI4Gya+kThdGQpsAK4I9MFHObaTgH8D64Bqd/NvcfrdY/7+1RPfZOLg/XNjLMAZ+E3G+bL6jKre4v69PIXT1bIK+Jn7rTte4lsK5OF0Q68GrgwaQI5LrSIJGGOMCa01dAcZY4wJw5KAMca0YpYEjDGmFbMkYIwxrZglAWOMacUsCZiEJyK3ichIEZkoYarEuvO63xeRVSLyowi+9sUi0jXo8YMi0jtSx29EHMtEZLB7f4OIdIp2DCYxWRIwLcEJOHPwT8aZ/x7KaOAzVR2gquH2aYqLgZokoKqXxVOFWmMaYknAJCwRuUNE1uLUaXkXuAy4T0Rm1tqvPzAXOM2t8Z4hInuD2s8VkUfc+4+IyN0i8o6IfC0i5wbtd7046z6sEZE5bttg4Img4wZ/I5/s7v+RiNwedJy9IvIH9zjviUjnEL9bpog87D5/rYic426/T0QKJaiGvTHNYUnAJCxVnY7zwf8ITiJYq6oFqnpLrf1WAzOBp1W1v6qWNnDoLsBJwBnAHAARORWnDMAJqtoPmOvWjC8EptQ+rttFdDtO6eP+wBARCZQ9bge85x5nOXB5iBhuAnapal9VLQCWutt/p6qDgQLgZPfKVWOazJKASXQDcC7PPxaIVDfMIlWtdrt1At/SxwAPq+p+AA/lHoYAy1S12C19/AQQqA57AHjOvb8C6BHi+WOAewMPVHWne/d8EVmJUzLheCDq4w+mZUlpeBdj4o/bxfMITiXJ7UBbZ7OsBk708G0/uF5Keq224Fo0EnTbmBoroUqYB1TowXotVYT+O6zzeiLSE/g1MERVd7pdWLVjN6ZR7EzAJCRVXe3Wcv8C59vwUmCcx+4ecMo6HyciSTgVMxvyCnCpiLQFZ71gd/seICvE/u/jdNd0EmeJ08nAmx5eJ/j1rgo8cBd4yQb2AbvccYRTG3E8Y0KyJGASlojkATvd+vfHNnJWzgycLpmlOAuo1EtVX8IpBV3onm382m16BLg/MDActP8WnAVG3gDWACtVtTFlj38PtHcHldcAo1R1DU430MfAQ8DbjTieMSFZFVFjjGnF7EzAGGNaMUsCxhjTilkSMMaYVsySgDHGtGKWBIwxphWzJGCMMa2YJQFjjGnFLAkYY0wr9v8BXAzpJHr8KdoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(sum_of_squares_nojac, 'o-', label=\"estimated Jacobian\")\n",
    "plt.plot(sum_of_squares_jac, 'o-', label=\"analysic Jacobian\")\n",
    "plt.ylabel(\"sum of squared residuals\")\n",
    "plt.xlabel(\"# function call\")\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is about what was expected. To estimate the Jacobian, one needs one additional call per parameter, because the parameters will be slightly varied one at a time to estimate the derivatives. With the analysic Jacobian available, this is not necessary. The two extra calls in the beginning give also a little plateau to the residue curve where the Jacobian was used, which must be somehow related to the details of the algorithm in the beginning."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}