doc/getting_started.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Getting started"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Fitting a data to a user defined model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "*refnx* can examine most curve fitting problems. Here we demonstrate a fit to a simple user defined model. This line example is taken from the [emcee documentation](http://emcee.readthedocs.io/en/stable/user/line.html) and the reader is referred to that link for more detailed explanation. The errorbars are underestimated, and the modelling will account for that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "To use *refnx* we need first need to create a dataset. We create a synthetic dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "np.random.seed(123)\n",
    "\n",
    "# Choose the \"true\" parameters.\n",
    "m_true = -0.9594\n",
    "b_true = 4.294\n",
    "f_true = 0.534\n",
    "\n",
    "N = 50\n",
    "x = np.sort(10 * np.random.rand(N))\n",
    "yerr = 0.1 + 0.5*np.random.rand(N)\n",
    "y = m_true*x + b_true\n",
    "y += np.abs(f_true * y) * np.random.randn(N)\n",
    "y += yerr * np.random.randn(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We create a `Data1D` object from this synthetic data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from refnx.dataset import Data1D\n",
    "data = Data1D(data=(x, y, yerr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Then we need to set up a generative model. Firstly we write a fit-function that returns our straight line model. The `Parameter` objects describe the parameters we're going to use in the fit. We give the parameters values, names, and specify their limits. The parameters are combined into a `Parameters` set, `p`, using the or operator.\n",
    "Then we create a `Model` object from our parameter set and the fit-function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from refnx.analysis import Parameter, Model\n",
    "\n",
    "def line(x, params, *args, **kwds):\n",
    "    p_arr = np.array(params)\n",
    "    return p_arr[0] + x * p_arr[1]\n",
    "\n",
    "# the model needs parameters\n",
    "p = Parameter(1, 'b', vary=True, bounds=(0, 10))\n",
    "p |= Parameter(-2, 'm', vary=True, bounds=(-5, 0.5))\n",
    "\n",
    "model = Model(p, fitfunc=line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We set lower and upper limits on each of the parameters. This means that the log-prior probability from those parameters is described by a uniform distribution. Only solutions which have finite probability (i.e. lie between the limits) will be considered by the fit/sampler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2.3025850929940455 -inf\n"
     ]
    }
   ],
   "source": [
    "# ln(1 / 10)\n",
    "# a value lying outside the limits is not possible\n",
    "print(p[0].logp(1), p[0].logp(-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "It's not required to give each parameter bounds unless the fit method requires it. The bounds do not have to be from a uniform distribution, any of the `scipy.stats.rv_continuous` distributions can be used:\n",
    "\n",
    "```\n",
    "import scipy.stats as stats\n",
    "# a normal distribution of mean 5 and standard deviation 1.\n",
    "p[0].bounds = stats.norm(5, 1)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now we create an `Objective` from the model and the data. We use an extra parameter, `lnsigma`, to describe the underestimated error bars. Objectives use the model and data to calculate statistics about the curve fitting system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1185.7110246179552 -6.405228458030841 -738.614995170284 -745.0202236283149\n"
     ]
    }
   ],
   "source": [
    "from refnx.analysis import Objective\n",
    "\n",
    "lnf = Parameter(0, 'lnf', vary=True, bounds=(-10, 1))\n",
    "objective = Objective(model, data, lnsigma=lnf)\n",
    "print(objective.chisqr(), objective.logp(), objective.logl(), objective.logpost())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Then a `CurveFitter` is created from the `Objective`. This is responsible for doing all the curvefitting/Bayesian sampling. Let's do a quick fit using Differential Evolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "22it [00:00, 127.49it/s]\n"
     ]
    }
   ],
   "source": [
    "from refnx.analysis import CurveFitter\n",
    "\n",
    "fitter = CurveFitter(objective)\n",
    "fitter.fit('differential_evolution');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In the final output of the sampling each varying parameter is given a set of statistics. `Parameter.value` is the median of the chain samples. `Parameter.stderr` is half the [15, 85] percentile, representing a standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________________________________________________\n",
      "Objective - 4829542384\n",
      "Dataset = <None>, 50 points\n",
      "datapoints = 50\n",
      "chi2 = 46.739546242716216\n",
      "Weighted = True\n",
      "Transform = None\n",
      "________________________________________________________________________________\n",
      "Parameters:      None      \n",
      "<Parameter:     'lnf'     , value=-0.790441 +/- 0.205, bounds=[-10.0, 1.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters:      None      \n",
      "<Parameter:      'b'      , value=4.52831 +/- 0.358, bounds=[0.0, 10.0]>\n",
      "<Parameter:      'm'      , value=-1.00301 +/- 0.078, bounds=[-5.0, 0.5]>\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAfu0lEQVR4nO3deXxU1d3H8c+PREBFH1FQgYDBHVQUjEpwKYuCFgVxQawbosQNRcVqsajl0SruolJsAFFcqIj4oJayFEGpRDSAC4siguyWaN2rIHKeP04iSUjITOZO7tzM9/165dWZITn3MJZvzpx7zu+Ycw4REYmuOmF3QEREEqMgFxGJOAW5iEjEKchFRCJOQS4iEnGZYVy0UaNGLjs7O4xLi4hE1vz5879wzjUu/3ooQZ6dnU1hYWEYlxYRiSwzW1XR65paERGJOAW5iEjEKchFRCJOQS4iEnEKchGRiFOQi4hEnIJcRCTiFOQiIhEXqSA3818FBWH3REQkdQQS5GZ2g5ktNrNFZjbezOoH0W5ppcO7SxeFuYhIiYSD3MyaAdcBOc65w4EMoE+i7ZY3e/a2x5s3l30uIpLOgppayQR2NrNMYBdgfUDt/qpjx22P69Yt+1xEJJ0lXDTLObfOzB4AVgM/AtOdc9PLf5+Z5QF5AC1atIj7Orm5MHeuH4l37Oifi4gIWKKHL5tZQ+Al4Dzga+BFYKJz7tnKfiYnJ8ep+qGISHzMbL5zLqf860FMrZwMrHTOFTnnfgYmAR0CaFdERGIQRJCvBtqb2S5mZkAXYGkA7YqISAwSDnLn3DxgIrAA+LC4zfxE2xURkdgEckKQc+4O4I4g2hIRkfhEameniIhsT0EuIhJxCnIRkYhTkIuIRJyCXEQk4hTkIiIRpyAXEYk4BbmISMQpyEVEIk5BLiIScQpyEZGIU5CLiEScglxEJOIU5CIiEacgFxGJOAW5iEjEKchFRCIukCA3sz3MbKKZfWRmS80sN4h2RUSkaoEc9QYMB6Y6584xs7rALgG1KyIiVUg4yM3sf4CTgL4AzrnNwOZE2xURkdgEMbXSEigCxprZQjMbbWa7lv8mM8szs0IzKywqKgrgsiIiAsEEeSbQDhjpnGsL/AD8ofw3OefynXM5zrmcxo0bB3BZERGBYIJ8LbDWOTev+PlEfLCLiEgNSDjInXOfA2vM7JDil7oASxJtV0REYhPUqpVrgeeKV6ysAC4NqF0REalCIEHunHsPyAmiLRERiY92doqIRJyCXEQk4hTkIiIRpyAXEYm4tA1yM/8lIhJ1aRvkIiK1hYJcRCTi0j7ICwrC7oGISGLSMshLh3eXLgpzEYm2tAzy2bO3Pd68uexzEZGoScsg79hx2+O6dcs+FxGJmrQM8txSJ4rOnFn2uYhI1KRlkJemEBeRqEv7IBcRibqg6pFHjnNh90BEJBgakYuIRJyCXEQk4gILcjPLMLOFZvZaUG2KiEjVghyRDwSWBtieiIjEIJAgN7MsoDswOoj2KrVuHWzcGGOfVKZWRNJDUCPyR4Cbga2VfYOZ5ZlZoZkVFhUVVe8qd94JWVnQuzfMmAFbK72ciEjaSDjIzex0YKNzbv6Ovs85l++cy3HO5TRu3Lh6Fxs4EAYM8Nsxu3aFAw+Eu++GDRuq156ISC0QxIj8eKCHmX0G/A3obGbPBtDu9lq1goce8lMszz8P2dnwxz9C8+bQqxdMmQK//JKUSyebpoJEpLoSDnLn3GDnXJZzLhvoA7zunLsw4Z7tSP36cP758PrrsGwZDBoEc+dC9+7QsiUMHUoWa5LaBRGRVBH9deQHHQT33gtr1sCLL/pR+5/+xGdk8xrd+ejeybBlS9i9jJlqo4tIvAINcufcbOfc6UG2GbO6deGcc2DaNBZMXME9DKYtCzn0D2eyuUkLGDIEVq4MpWtV0UEXIpKI6I/IKzBtWUtu4y5asJpedf6P1Xu1g3vugQMOgG7dYOJEf6JEitBBFyKSiFoZ5CUHRfxCJtPq9aRo7Gvw2Wdwxx2wdCmce66/QXrLLfDJJ4Feuzo3LXXQhYgkwlwIZQBzcnJcYWFhUq9REqZz55arOf7LLzBtGowaBa++6p936gT9+8NZZ0G9eoFcN963tdL+iogUM7P5zrmc8q/XyhF5aduFYkYG/Pa38PLL/gbpn//sR+u/+x00awY33uhH7QmKd57bOf9VWYhreaKIVKbWB/kONWkCt94Ky5fD9OnQuTM8/ji0bg0nngjjxsGPP8bcnG5aikgY0jvIS9SpA6ecAhMmwNq1cN998O9/wyWX+LC/9lr44IMqm9FNSxEJQ60N8pKpirjtvTf8/vfw8ccwa5bfZDRqFBx5JBx3HIwZA99/X+GP6qaliISh1gZ5wsx8Ej/3nC8J8PDDPsAvvxyaNoUrr4T5ZcvLlJ7fnjkzOTctNV0jIuUpyGOx115w/fWwaBG89ZZf3TJuHOTkQLt2MHIkfPMNUPVNy+rQ3LuI7IiCPB5m0KEDPPUUrF/vb4xu3QpXX+1H6f36+ZQNeEmn5t5FZEcU5NW1xx5wzTWwcCHMm+eXL06Y4IO+TRt49FH46qtALqW5dxHZEQV5oszg2GP9DdENG+Cvf/XVGQcO9KP0iy6COXMSGqXXxNy7iESXgjxIu+0GeXnw7ruwYAFceim88gqcdJJfm/7gg/DFFwldQiEuIuUpyJOlbVv4y1/8XPrYsdCwIdx0k9892qePH1rrqDoRCYCCPNl23RX69vVFVD780C9bnD4dTj4ZDj4Yhg2Dzz+vsplqr4sXkVpPQV6TDj8chg/369KffdaPzgcP9pUYzz7bF/PSKF1E4qQgD8POO8MFF8Abb/gCXQMHwptvwqmnwv77w513+rAXEYlBwkFuZs3NbJaZLTGzxWY2MIiOpY1DD4UHHoC1a+nNC8xYdRDcfju0aAE9esBrr0XqqDoRqXlBjMi3AIOcc62B9sA1ZtY6gHbTS716vEhvujLDV2O8+WZ45x044wzIzvbhvmpV0i6vMrki0ZVwkDvnNjjnFhQ//g5YCjRLtN10VrDxAH803Zo1MGkSHHEE3HUXtGwJp53ma6n//HPY3RSRFBHoHLmZZQNtgXkV/FmemRWaWWFRUVGQl60VKqynstNO0KsX/OMf/uDo227zK1/OOstPvQweDJ9+GlqfRSQ1BBbkZtYAeAm43jn3bfk/d87lO+dynHM5jRs3DuqytUaV9VT22w+GDvWnGb36KhxzjK+bfuCBfinjhAmwaVPM16tsKkUFuUSiJ5AgN7Od8CH+nHNuUhBtppuY66lkZsLpp/sdo6tX+xUuy5fDeedBVpavpb5sWVzXVnVFkWgLYtWKAWOApc65hxLvUnqqVj2VZs1gyBA/vTJ1qi8F8MgjcMgh22qpx3BUnaorikRbECPy44GLgM5m9l7x128DaDdtxV1PJSMDunWDl17yN0iHDfNH1l14oQ/7gQN9LfVKqLqiSLSZC2Hfd05OjissLKzx66aVrVv90HrUKL/yZfNmaN/eF/Xq3RtrsCvgKwfk5m6bLy95LiKpx8zmO+dyyr+unZ21VZ060LkzjB/vd4k++CB8/TX068eWfZoygqs5ioXbzYkrxEWiR0FeS5VZldKoEdx4IyxZAnPmsPTgnlzKWBbSjjd+PIaxHfJpwHeh9ldEqk9Bnk7M4IQT+H7EOJqyngE8Rj02kc8VbKAJ+fT3u0lVZlEkUhTkAUjl7e0V9S03F76mISMYwJG8T3sKaNDvPPrv8jwcd5yvpT5ihJ+KEZGUpyBPe8Y82lNw+Rh/VN3IkX4VzIAB/qi6vn3hrbc0ShdJYQryNFHVJp8uXaBg8e7+4Iv586GwEC6+2K94OeEEX0v9kUfgyy/jum7JJwJtMhJJHgV5gFI5rMqvTnEO7r572/PtNgIdfTQ88YQ/qm70aH8e6Q03+HXpF1zgv7mKUbp2jIrUDAV5glIxrCrqQ0U7NmPaCNSgAVx2Gbz9Nrz/PvTvD3//O3Tq5Gup338/bNxYYT+0Y1SkZijIE5SKYVVRHyoK6rjLArRpA489BuvXczFP869ljX3d9Kws6N0bZswoc1Rdx47+MKSMDO0YFUkmBXmCUnF7e0V9qCqo49oItMsuPMPFnMi/YPFif2P09deha1dfjfHuu2HDBnJz/XXvvDOO+jEiEjdt0Q9AKm5vL7/ksLL/zCXfF+//Dbb7uZ9+8gdejBoFs2b5YfgZZ/ipmG7d/HMRSYi26NeAVAnxeDgX0MrC+vXh/PP9yHzZMhg0yP9m697dn2w0dKgv6CUigVOQ13KBBXUlKry5e9BBcO+9PrhffBFatfJBnp3tg33y5F8PlC69YanksZYrisRHQR6AZIdlmCraGRrzSp26deGcc2DaNFixAm69FRYuhDPP9EfVDRlCNiu3axOgQ4eywS4ilVOQS9yqs1LHWmZjd93pTzWaPBnatcPdcw+fcgBT6cbwkyayE5sD62P5XwD6hSC1mYK8lkrmp4SEVupkZkKPHvDaa4y46TOGcgetWMrftpzLGpozjFs4iPiOqhNJd0Gd2XmqmX1sZsvN7A9BtCmpq1rH0lXg6DOb87/cQUtW0iNzCnPpwCAeZBmH8DqdOJ/nqcdPwXQ6ThrBS5QEcWZnBjACOA1oDZxvZq0TbVdSS2Xz4Ims1Cn52a1kMPjN09h37ss0Zw238mf2YxXPcwHraMZD3MB7zy+p1jXK9zvem6i66SpREMSI/FhguXNuhXNuM/A3oGcA7UrIklF+oHwbJVNAubn+63OacA+3ciDLOYXpzKQL1zCCoy44jG+PPAHGjYP//jeufufnx/f3SKWyC/pkILEIIsibAaUXCK8tfk0iLqjyA9UJRkcd/skpnMcEsljLLXYfW9ZthEsu8eV1r70WPvggpn6/9FLZ51X9PVKx7II+GciO1NjNTjPLM7NCMyssKiqqqctKAoIqP5BoMBaxN4/V/z0fv/Kx3zXavbvfQXrkkf4gjDFj4PvvK+332WeXfV7V3yNVyi6k0icDSW1BBPk6oHmp51nFr5XhnMt3zuU453IaN24cwGUl2XZ0UzOeVTGxBuOOgmrmTMjtYP6Hn3vOHyj98MM+wC+/HJo0gSuugPnzt+t3Xl7lf4+KBHUzN1Gp+MlAUlMQQf4ucJCZtTSzukAf4JUA2pUUEsRNTdhxMO4oqLb7mb32guuvh0WL/AlGZ58NzzwDOTnQrh1XMpLd+Wa7n4v171F67j4sqfLJQIKTrHseCQe5c24LMACYBiwFJjjnFifartROOwrGagWVmd8G+tRT/hCMESNg61ZGcjXracrG0/v5WupEb+ttqnwykOAFPU2m6odSI2Ktslh6tOJc9aozFsx1XHd8IXnkcz7jacAP/qi6vDy48EJo2DC+zoeoutUpJfUUFPgxB/g6/dX55azqhxJJ1dmhOvsNo5BjyGMUWXU2MKVXvv+Xc9110LQp4+xiTrQ5FMxN/XSszXV80k0y73koyKVG1GQglZ6i2VxvNxr+vj+88w4sWMDnp11KTyYzh5NoeEJrPrvuIfjii5rpWJy0hrx2SeY9DwW5pJSSwE8k9HNzfSn0u+8u9/G1bVvGHvMXmrKevozlK9eQ7McG+QOl+/Tx31zqqDqRICXznofmyCWtlJ+nnJu/iKPeHeVXvHz1FRxwgF/O2Lcv7LtvqH3V/Hjtk+h/U82Ri7D9qOioCw+H4cPZ+at1XMgzfnQ+eDA0bw5nnQVTp8Ivv4TXYalVkjXFqCCXtFPRGvGf2JnnuJCCYW/A0qUwcCDMmQOnnQb77+9PkF633T43kZSgIJe0t91W+K8OhQcegLVr4YUX4OCD4fbb/alGPXrAq6/+elRdTfdPpCIKckl7FS0LMwOrXw9694YZM2D5crjlFr/6pUcPf/7o7bfDqlVJ6ZPqrEg8FOSS9mJaFnbAAX4ZzJo1MGkStGkDd90FLVv66ZdJk+DnnwPrU7LWHGtJY+2kIJe0V+lyxYrstBP06gVTpsDKlXDbbfDhh77WS/Pm/kbpp58m3Kdk11nRCL920fJDkQrEtUxsyxa/uiU/H/7+d78WvUsX6N8fzjwT6tVLqA9z5waz5jiILeISLi0/FKmGmEaumZlw+unwyiuwerVf4bJ8ud9klJUFN90EH39c7T4EFbYqi1t7KchFyknoRmOzZjBkiJ9emToVfvMbGD4cDj3UP372Wfjxx8D7HAuVxa29FOQi5QQycs3IgG7dYOJEv4xx2DC/Dv2ii3zYDxzoa6lXIFk3JFUWt/ZSkIuUE/jIdZ99/NLFZct8gnbrBk88AUcc4Setx46FH35I8CKxSYUDMyR4CnKRcpI2cq1TBzp3hvHj/ej8wQd9fZd+/fyB0ldfDQsX/vrtKmErsVKQi+xA0kaujRrBjTfCkiXw5pvQs6cfmbdrxzscQ3/y4bvvknRxqW0SCnIzu9/MPjKzD8zsZTPbI6iOiaSqQOewzeDEE2HcOH9U3fDh1Ocn8rnCHyh9+eV+N2kNDs21aSh6Eh2RzwAOd861AZYBgxPvkkj6KAlNM/wRdNddRxs+oD0FvjzA+PFw3HHQtq0/j/Trr5PeH4mehILcOTe9+PBlgLeBrMS7JBK+cOenjXm0p6D/k7BhA4wc6VfBDBjg59L79oW33tIEuvwqyDnyfsA/AmxPJK0UFFSwhn3x7nDllTB/PhQWwsUX+7ouJ5zgD5R+5BH48suk9UeiocogN7N/mtmiCr56lvqePwJbgOd20E6emRWaWWFRUVEwvRcJUdBB16WLnyovsd0a9qOP9ssW16+H0aNht93ghhv8uvQLLvDfnMAovfzfp1MnhXlUVBnkzrmTnXOHV/A1GcDM+gKnAxe4HRRucc7lO+dynHM5jRs3DuwvIFKTgiwvW/5nN20q+7zSNewNGsBll8Hbb8P77/uaLlOm+OQ99FC4/37YuDHu/pTf+KRt/NGR6KqVU4GbgR7Ouf8G0yWR1BVkvZLyP5uR4WdOSsS0hr1NG3jsMT9KHzcO9t4bbr7Z13g591xfSz3GA6XL/9LQNv7oSHSO/HFgN2CGmb1nZk8E0CeRlBXkrs/yP/v442WDO6417Dvv7Lf/z5nj16Zfey3MmgVdu/JT8wN5o+ufKXxl/Q6bKH+9WbO0AzQqVMZWJE5BlpetqK1ET1r/1aZNLLv3ZdbcMYouvM4WMvj2pDPY85b+vkxARkal/Qnk+hI4lbEVCViQo9WkjHzr1eOlnfpwMjM5iGU8ZDdRb/5c6N7dn2w0dKg/8agUhXc0KchFarGS6ZvlHMSf6g/jwylrfEXGVq18kGdn+1rqkyfX6IHSYSrZgFWbVuQoyEVCVNHGoyA3I5U/xq79SXX9sXTTpsGKFXDrrbBggT/JqEULGDIEt2JlrR2Z19ZDrRXkInGKWlXC3Fx/lOh20zfZ2f40o9Wr/Yi8XTu45x5/0HTXrn7kvnlzGF1Omtp6SpKCXCTdZWZCjx7w2mvw2Wdwxx3w0Ud++WJWll/OuGxZ2L0MRG09JUmrVkRkOxn2C12Zzj96jfJnkf7yi0+9/v3hrLOgfv2wu1htBQV+JN6xY/SWV1a2akVBLiLbKbMEcsMGeOopXxZgxQrYc0+/c6l/f2jdOsxuph0tPxSR6mnSxE+yf/KJ3yl68sm+pO5hh/niXU8/Df/Vxu4wKchFJDZ16vgQf+EFf1Td/fdDUZEvq9u0qS+z+/77oXYxXQ/FUJCLSPwaN4abbvI3RWfP9mvRR4+Go47yB2GMHg3ffx92L9OGglxEKlXlOmsz+M1v4NlnfeGuRx6BH37w8+dNmsAVV/g66lFarxlBCnIRKaOiTTMxTVnsuScMHAgffuhPMDrnHHjmGTjmGF9LfeRI+OabpPY9XSnIRaSMhDfNmEGHDjB2rB+ljxjhR+RXX+3n0vv1878dkjhKry07NmOlIBeRMgLdNLPHHj7AFyyAd9/1Jxm9+KIP+jZt4NFH4T//SbDHXm3dfh8LBbmIlFF6k0xMh1vEwgxyciA/34/S8/N9DfWBA/0o/aKL4M03Exql19bt97FQkItIpcqHeCCj3N128zdD33kHFi70x9a98oq/adqqFTz4oF/WGKfauv0+FgpyEdmhpE5ZHHWUn0Nfv97Pqe+1l1/W2KwZ9OnjPxLEeFRdUj5JREQgQW5mg8zMmVmjINoTkdRRI1MWu+7qNxa99RYsWuTn1adP9xuQDj4Yhg2Dzz+Publ0CnEIIMjNrDnQFVideHdEJJWYhTBlcdhhfj36+vV+fXpWli8R0Ly5r6U+daov4iW/CmJE/jBwM6AV/yK1UGhTFvXr+1Uus2f7HaTXX+9viJ52Guy/v6+lvm5dDXUmtSUU5GbWE1jnnAu3wIKIBKr04pHSc+KhTVkccoiv7bJuHUyY4Kdbbr/dn2rUowe8+ips2RK5Qz+CUmUZWzP7J7BvBX/0R+BWoKtz7hsz+wzIcc59UUk7eUAeQIsWLY5etWpVIv0WkSQqKPBLvcGvEvzxR/84pUJyxQoYMwaefNLPnzdr5jcbXXYZ7Ldf2L1LisDrkZvZEcBMoKR+ZRawHjjWObfDuxKqRy6S2u65xx/nCZCRsW1KOqWCvMTPP/vTjUaN8vPnAN26+SWOZ5wBO+0Ubv8CFHg9cufch865vZ1z2c65bGAt0K6qEBeR1Ff+BufcuSka4uCDulcvmDIFVq6E227z9V7OPtvfIB08GD79NOxeJpXWkYvIdiK7Jnu//WDoUH/26Kuv+pK6998PBx64rZb6pk1h9zJwgQV58ci8wvlxEYmuyIR4aZmZvkb65MmwapVf4bJ8ud9klJUFgwb5lTC1hEbkIlKhWrMCpFkzGDLE3xydNs2XAnj0UV8O4KST/Fr1kru5MUq1k4gU5CISSXGHaZ060LUrTJwIa9f63aIbNviCXU2bwnXX+bn1CFKQi0j62WcfuOUW+PhjeP11OPVU+OtffWnd3Fxf9+WHH6psJlVK5SrIRSR91akDnTrB+PF+s9FDD8HXX/v16E2bwlVX+VrqpaRi3XMFuYhEWmBB2qgR3HADLFkCc+bAmWfCU0/5Y+pycvyI/dtvU7LuuYJcRCInqaNiMzjhBHj6aV+467HHfGJfeSU0bUq/gss5lnmAS5m65wpyEYmcGhsVN2wIAwbA++/D22/Deeexz8zxzKM973EUi696nNxWXyfp4rFTkItI5NR4aV0zv7lozBjYsIErGckWMmn50LV+Lv2SS3wt9ZDWayrIRSRyQt15uvvuPOGu5Gg3HwoL4eKL4eWX/XTMYYfBww/Dl1/WYIcU5CIScaHuPD36aHjiCT+XPmYM7L473HijH6X/7ncwa1aNjNIV5CIiiWrQwC9ZfPttP5+el+eLeHXu7Gup33cfbNyYtMsryEUkklK2hECbNn6ly/r1fuXL3nv7zUdZWXDuuf5M0oApyEVEkmGXXfz8+b/+BYsX+9Uvs2bBTz8FfqlqHyyRCB0sISJpadMmv8ymmhW3KjtYIjPhjomISGzq1UtKs5paERGJOAW5iEjEJRzkZnatmX1kZovN7L4gOiUiIrFLaI7czDoBPYEjnXObzGzvYLolIiKxSnREfhUwzDm3CcA5l7wV7yIiUqFEg/xg4EQzm2dmb5jZMUF0SkREYlfl1IqZ/RPYt4I/+mPxz+8JtAeOASaY2f6ugsXpZpYH5AG0aNEikT6LiEgpVQa5c+7kyv7MzK4CJhUH9ztmthVoBBRV0E4+kA9+Q1C1eywiImUkOrXyf0AnADM7GKgLfJFop0REJHaJ7ux8EnjSzBYBm4FLKppWERGR5EkoyJ1zm4ELA+qLiIhUg3Z2iohEnIJcRCTJzKpd8DAmCnIRkRpSUJCcdhXkIiJJVDq8u3RJTpgryEVEkmj27G2PN28u+zwoCnIRkSTq2HHb47p1yz4PioJcRCSJcnO3PZ45s+zzoOioNxGRJEv2NkmNyEVEIk5BLiIScQpyEZGIU5CLiEScglxEJOIU5CIiEacgFxGJOAW5iEjEKchFRCLOwjiZzcyKgFU1fuFwNULnmYLehxJ6Hzy9D/G9B/s55xqXfzGUIE9HZlbonMsJux9h0/vg6X3w9D4E8x5oakVEJOIU5CIiEacgrzn5YXcgReh98PQ+eHofAngPNEcuIhJxGpGLiEScglxEJOIU5ElmZs3NbJaZLTGzxWY2MOw+hcXMMsxsoZm9FnZfwmJme5jZRDP7yMyWmlkSDv5KfWZ2Q/G/h0VmNt7M6ofdp5pgZk+a2UYzW1TqtT3NbIaZfVL8vw3jbVdBnnxbgEHOudZAe+AaM2sdcp/CMhBYGnYnQjYcmOqcOxQ4kjR8P8ysGXAdkOOcOxzIAPqE26sa8xRwarnX/gDMdM4dBMwsfh4XBXmSOec2OOcWFD/+Dv8Pt1m4vap5ZpYFdAdGh92XsJjZ/wAnAWMAnHObnXNfh9ur0GQCO5tZJrALsD7k/tQI59ybwH/KvdwTeLr48dPAmfG2qyCvQWaWDbQF5oXbk1A8AtwMbA27IyFqCRQBY4unmEab2a5hd6qmOefWAQ8Aq4ENwDfOuenh9ipU+zjnNhQ//hzYJ94GFOQ1xMwaAC8B1zvnvg27PzXJzE4HNjrn5ofdl5BlAu2Akc65tsAPVONjdNQVzwH3xP9iawrsamYXhtur1OD8evC414QryGuAme2ED/HnnHOTwu5PCI4HepjZZ8DfgM5m9my4XQrFWmCtc67kE9lEfLCnm5OBlc65Iufcz8AkoEPIfQrTv82sCUDx/26MtwEFeZKZmeHnRJc65x4Kuz9hcM4Nds5lOeey8Te1XnfOpd0IzDn3ObDGzA4pfqkLsCTELoVlNdDezHYp/vfRhTS86VvKK8AlxY8vASbH24CCPPmOBy7Cj0LfK/76bdidktBcCzxnZh8ARwF3h9yfGlf8iWQisAD4EJ9DabFV38zGAwXAIWa21swuA4YBp5jZJ/hPK8Pibldb9EVEok0jchGRiFOQi4hEnIJcRCTiFOQiIhGnIBcRiTgFuYhIxCnIRUQi7v8BSnm79uN1MLAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(objective)\n",
    "objective.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The trouble with a single fit that minimises $\\chi^2$ is that it doesn't reveal the range of solutions that are consistent with the data. It also assumes that the parameter uncertainties will be normally distributed and uni-modal. To investigate the parameter probability distributions we need to use MCMC to sample the posterior probability distribution of the system. \n",
    "\n",
    "Note: `pool=` specifies that no parallelisation is done during sampling. On platforms that use `spawn` for multiprocessing special precautions [must be used](https://docs.python.org/3/library/multiprocessing.html#the-spawn-and-forkserver-start-methods) when using `pool > 1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [00:20<00:00, 48.70it/s]\n"
     ]
    }
   ],
   "source": [
    "fitter.sample(1000, pool=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Once the sampling is done we burn/discard some of the initial steps because the initial locations of the walkers won't be around their 'equilibrium' position. We thin out the chain to reduce auto-correlation between successive steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________________________________________________\n",
      "Parameters:      None      \n",
      "<Parameter:     'lnf'     , value=-0.761916 +/- 0.153, bounds=[-10.0, 1.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters:      None      \n",
      "<Parameter:      'b'      , value=4.53947 +/- 0.355, bounds=[0.0, 10.0]>\n",
      "<Parameter:      'm'      , value=-1.0064 +/- 0.0775, bounds=[-5.0, 0.5]>\n"
     ]
    }
   ],
   "source": [
    "from refnx.analysis import process_chain\n",
    "\n",
    "process_chain(objective, fitter.chain, nburn=300, nthin=100, flatchain=True)\n",
    "print(objective.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now we can see the range of solutions that are consistent with the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3Sc13nv+++e3jHojQB770WsEkVZki3bkpPIjmtck6OcJDfNTnJc4pz4JOfaXonPshOv2NGJ5CqXa0u2FVtRrxR7J0WKnSBBkAQIYIDBYPrs+8dgvxxAYAExIDDk81mLi4P24p2R+MPGs5+9t9JaI4QQonTZxvsGhBBCjI4EuRBClDgJciGEKHES5EIIUeIkyIUQosQ5xuObVlVV6SlTpozHtxZCiJK1c+fOi1rr6qHvH5cgnzJlCjt27BiPby2EECVLKdUy3PultCKEECVOglwIIUqcBLkQQpQ4CXIhhChxEuRCCFHiJMiFEKLESZALIUSJkyAXQogSV1JBrlT+z+bN430nQggxcRQlyJVSf6mUekMpdUAp9WOllKcY1y1UGN533y1hLoQQxqiDXCnVCPwZsEJrvQCwAx8c7XWHevnlS49TqcFvCyHEraxYpRUH4FVKOQAf0Fak61o2bLj02OUa/LYQQtzKRr1pltb6rFLqn4HTQBx4Vmv97NDPU0o9BDwE0NzcPOLvs2YNbNqUH4lv2JB/WwghBKjRHr6slCoHHgc+AESAnwE/11r/8HJfs2LFCi27HwohxMgopXZqrVcMfX8xSiv3ACe11h1a6zTwBLC2CNcVQghxDYoR5KeB1Uopn1JKAXcDh4pwXSGEENdg1EGutd4K/BzYBewfuObDo72uEEKIa1OUE4K01v8T+J/FuJYQQoiRKamVnUIIId5KglwIIUqcBLkQQpQ4CXIhhChxEuRCCFHiJMiFEKLESZALIUSJkyAXQogSJ0EuhBAlToJcCCFKnAS5EEKUOAlyIYQocRLkQghR4iTIhRCixEmQCyFEiZMgF0KIEidBLoQQJa4oQa6UCiulfq6UelMpdUgptaYY1xVCCHF1RTnqDfgG8LTW+n1KKRfgK9J1hRBCXMWog1wpVQasBz4BoLVOAanRXlcIIcS1KUZpZSrQAXxHKbVbKfUfSin/0E9SSj2klNqhlNrR0dFRhG8rhBACihPkDmAZ8C2t9VIgBnx26CdprR/WWq/QWq+orq4uwrcVQggBxQnyVqBVa7114O2fkw92IYQQN8Cog1xrfR44o5SaPfCuu4GDo72uEEKIa1OsrpU/BR4b6Fg5AXyySNcVQghxFUUJcq31HmBFMa4lhBBiZGRlpxBClDgJciGEKHES5EIIUeIkyIUQosTdskGuVP6PEEKUumK1H94QuVwOrTVKKWy2W/ZnkBBCDFJSaaiUQimF1ppsNksulxvvWxJCiHFXckFus9mw2WxFC/TNm4t4g0IIMQ5KKsiN0QZ6YXjffbeEuRCitJVUkOdyObLZLFpr4PoD/eWXLz1OpQa/LYQQpaakgjybzZJMJslms9cU6ObjQ23YcOmxyzX4bSGEKDUl1bWSSqWIxWLY7XZcLhdutxu73W4FeOFkqNbaGpmbjxtrCk4UfeGFwW8LIUSpKakgz2QypNNpa0SeSqUGBboZmV9roIOEuBCi9JVUkKfTafr7+3E4HLhcLjKZDJlMhlQqhdPpxOVy4XQ6rynQQVYDCSFuDiUV5PF4nGg0itvtJpVKYbfbcbvdpNNp0uk0yWTSGqFfLdAzGQl0IcTNoaSCvL+/n0gkgs1mw+12EwgEBoW3w+Egm82STqdxOBx4PB4cDseIa+hCCFFKSirI0+k0XV1dBAIBMpkMiUQCj8eD1+slkUjgdrut0brD4SCdTmO32/F6vTidTgl0IcRNqWhBrpSyAzuAs1rr+4t13UI9PT10dXXR19eH1+u1Ary/vx+Px2OVVFwuFy6Xi3Q6jdPptCZFPR4PLpdLAl0IcVMp5oj8z4FDQKiI1xzEjMLT6TSJRIKenh68Xq9VYjGToCbkHQ6H9T7T5eJ0Oq3Ql0AXQtwMihLkSqlJwLuB/w18uhjXHE5/fz99fX3YbDa01ni9XnK5HLFYzBpxBwIB0uk0q1fHACc7d/rxeDzYbDY8Ho81MZpIJKyvKVxMJIEuhCg1xRqRfx34GyB4uU9QSj0EPATQ3Nx8Xd/E4XCQyWSsbpT+/n5rsjObzVpB7/V6AQ9QRldXEp/Ph8PhIJFI4PV6cblcVpAP7UWXQBdClJpRB7lS6n6gXWu9Uym14XKfp7V+GHgYYMWKFcOvnb8Ke0sL5V1dXCwrs4LWbreTSqXQWlu18Xg8Tj7Ie+jo8OPz+fB6vbjdbpLJJHa7HZ/PZ5Vq3G63FeamF10CXQhRKooxIl8HvEcp9S7y6RlSSv1Qa/17Rbj2IPN+8Qve9tRTnJw8mZ2LF3N4wQKSAyFrs9mIx+NWxwqkAA+dnQkikQjBYNDqcPF4PPT391uPTbui0+m0RuzD9aKPZaCby1xmexghhLisUQe51vpzwOcABkbkfzUWIQ7Q8qlPcUIp5m3dyvuefJLEM8+wb948di9bRltDA5ls1lrtCUkgRk+PE4fDYXW2OJ1OQqEQgUCARCJhhbcZsZuWRTMhaiZNh27KJSN0IcREUVJ95K5p0zj4oQ9x5MEHKdu3jzmbNrF0zx5W7t7Nhepqdi5ezN6FC4n5fEADMIMjR04weXIbDoeDVCqFz+ejv78fv9+Py+XC7/fj9/vp6+vD5/NZXS+m7BKPx62RvCm5jGWgb94s+78IIUZGXW6r17G0YsUKvWPHjhF/3f79+zlx4gTpdNoaeee6u5n0+uvM3byZ+tOnydrt7GpcxN+d/nue5T5sDvjQhx5h8uQ2qwRjdk80JRTTwmhG5qbf3AS4qZu73W4r7IeroZvXcqSBvnkzrF2bf+z1yo6MQojhKaV2aq1XDH1/SY3IQ6EQ4XCYeDxu7Uuedjo598ADtL7rXTgPH2bma68xa8s+/ovfopVGvpf5GIeOzibdkLaC1nS/JJNJa/+WaDSKw+EgGAwSCASsEbsJc7fbjcfjIZlMWitIPR7PoF0Xr3eEPtxBFxLkQohrVVJBXllZid1up6+vj76+PmKxGH19fSSTSQCSCxZweN48frnaw5v/lOaTfJ/P8lXsW3IcP9vE5rlzOTRvHqlcDqWUVW6x2+0kEgmrRbGnpwefz0cwGLRG6IWjdVNX9/v91mMT6Ha7DdBks9ce6HLQhRBiNEqqtJJIJEgmkwO7F+Zr2CbMe3t76evrs/Yr/8AHwkAl//CHu7jz9C9YsHUr5V1dJNxu9syZw7aFC2mtqyM10LFit9utFkan02ltuuVyuQiHwwSDQetjgUCAYDCIw+Gwyi2XRu1OQKH1yEou5t2bNsloXAgxvJuitGLq2/mRr90K1YqKCmKxGNFolN7eXnp6eoBeoI+Zd4U4lfoQR37nd6wJ0uX797N6717OVVWxee5c9ixYQJ/HY32fdDptLThyOBz09PQQCASsAE+lUvT29uL1egkGg9aEqMfjIb8mysPGjTbWrbNdc8nlaj9PpT1RCHE5JTUiB6xANBOe5mzOXC5HLpcjlUoRj8eZMSMKRHjyySjJZNLarzyTyUBPD82bNjF/2zYaW1vJ2GwcmD6drfPn80ZDA1mwesidTidOpxOtNT6fD5/Ph9vtpqysjEAgMKjscuxYNR/5yGRA4fXCU0+luP12u/XDZzSTohLkQoibYkQOWKNxUwoZGuqmJJLfu6uKWbNi9Pb2EolEiEbzoZ5xu2l997s5+fa34ztxgtmvv8783btZcvQo3YEAm2bNYseCBXQEg1YZRylljcQdDgfd3d34/X5rcjQQCPDCCxXkD6pQpFKap59OsHhxzqqrmy0ApA9dCFFMJRfkha4U6pFIFq2dZLNeKioqiMfjxGIxIpEIPT09RKNRUqkU2fnz2T9rFns++EFqt21j3ubNvHP3bt69axdHGxvZOHs2+2bOJDawHYDdbsfhcJBMJolGo0QiEcrKyvB6vTQ3u4ApgA2HQzN37gXOn8/i8/nw+/2DWhcl0IUQxVJypZVrcbnyi3mcSCSsCdKuri5isRjJZNIqwbja25m+cSOLduygIhKh3+Vi+/TpbJw9m9PV1TgKFga5XC5rz5dQKMTf/u0UYA7/8A8RFi2KWd0vZouAwn71oYF+pZKLTIYKIS5XWimpIDdhZxbiXOvXXC7UCztfurq66O7uttoZE4kEmVSK6kOHmLt5M3MPHsSVyXC2spLXZsxgx5w5dNlsA/u6MFDOwVpkFAqFKCsro6ysDI/HY43IQ6GQ1dZoul2GllyGBrosGBJCwE1UIy8MOhPmVwr2q9XUzSjZlF96e3u5ePEiPT09+SX6q1bx2oIFbOztZfLmzSzetYsPbt3K+7ZvZ09zMy9OncqhSZPQWlsLjXK5nHWdcDhMIBCwJkcjkYhVUw+Hw4PKLmZS1EyMmpLLSy9dOiRaFgwJIYYqqSA3oWyCfOgItvC0nyt9/XChbnrCA4EA1dXVxGL5SdLOzk56e3vpDwZpe+ABjtx1F2VnzjB382YW793LilOn6PL72Th9Oq/PmkWbx0M8HrfODE0kEgBWcIfDYWvTrq6uLoLBIOXl5VZro5kYLQz09etzgA0Al0vJgiEhxCAlVVq5nKGBDoNH6VcrwwwNdTOqNsv4zaRmd3e3tfCov7+fZDTKlAMHWLxjBzNPnMCmNQdra3lh8mR2TpmCze8nl8tZC4lyuRxer5fy8nL8fr+10CgQCBAKhaw/wWDQmhg1926z5Z/ba69p1q2TSVEhbkU3RY38WgwdrRtXG60Xfv3QUM9ms9ZkqGll7O7uJhqNkkgkiEaj2NraWLBzJ0t276aqp4eY08nmKVN4edo0TlZUkBs4ms78VmB2XqysrCQcDuPxeAgGg4TDYSorK61wN3u+OBz5EXk2e/2bcwkhStstE+RDjWa0br7OnPNpgt1sttXX10ckEqGrq4tIJEIsFqO/r4/aw4dZuH07Cw8fxpXN0lJWxguTJ/NqUxP9AwuKcrkcgUCATCaDz+ejqqoKv99PTU0NHo+HUChERUUFVVVVBAIBamv9gJtczj7o3kACXYhbxS0b5IUKw+96RutmpF4Y6olEwjor9OLFi3R2dlqhnu3qYuaOHSzft4/JFy6QttnY0dDAc01NvNHYSCaXsyY5s9kswWAQv99PVVUVZWVlVudLMBiktraWcDiM3++3dl0c+pwk0IW4uUmQD2M0o/WhoZ5KpUgmk1Z/uul8iUQiRCIRvMePs3TXLlYcOkQgmeSi18uLTU0829hIZ1kZLpeLXC5HcGA1aVlZGeFwmIqKCsLhsPV2ZWWlNXo3pxcVPheQQBfiZiVBfhWjGa2bUE+lUla4JxL5s0JN54s1WdrRweR9+7ht3z7mnjmDDdhXVcWzkyaxvamJxMDe5oFAgFQqRSgUorKykmAwSEVFBeXl5ZSVldHQ0EB5eTmhUAifz4fD4Ri0nwtIoAtxsxmzIFdKNQHfB2oBDTystf7Glb5mIgb5UFcarY8k1FOplLUzY3d3N11dXbS3txOJRLC3tTFvxw5WHThAdV8fEZz8iAbeXNPAkWCQwMDo3Jwh6vf7aWpqstoXKysrqa2tpaqqioqKCnw+n9W6KIEuxM1nLIO8HqjXWu9SSgWBncBva60PXu5rSiHIC13vaD2bzVoj9GQySSqVIhqN0tfXR2dnJ11dXVy8eJHuzk6q3niDwM928z5O4iXLiWCQZxsbebWpiZ6BHRjNqNuMzk3ppb6+nvLychoaGqwSzEgDXXZXFGLiu2GlFaXUr4Bvaq2fu9znlFqQDzXS0XphS2NhqBe2MnZ2dvLRj9oow8WXFz3GfWdfY2pnJ2ml2FRdzTONjWwNhfD4/QDWcXNVVVVW7byyspL6+nrq6+uprKwkFArhdrute7qWvVwkyIWYuG5IkCulpgCvAgu01r1DPvYQ8BBAc3Pz8paWlqJ93/E0XKjD5YN9uFCPxWK8+mqaj31sKqBwOrN87GPfZ2b8RZbu3s3q48cJpVK0u908U1/Psw0NtAcC1ipVs2FXdXU1DQ0NBINB6urqmDRpEnV1ddbeLjabzboHGBzoEuRCTHxjvteKUioAPA78xdAQB9BaPww8DPkRebG+73gbGtaFwT7c1gGF2wSYY+ICgQDHjkF+Gb4il1NofSdT3hPgyJo1bDxzhsZdu1h54AAfOXWKj546xe5wmN/U1fFCWRk2v5/+/n6i0Sitra2EQiHq6+s5deoUtbW1TJ48mYaGBmtvl/xJSwA5MhkGAj7/HDZvln1chCg1RRmRK6WcwK+BZ7TW/+dqn1/qpZVrNZLR+qZNmnXr8o+9XvjZz7qZOzffutjR0UFHRwdnz54le+oUc7ZuZe3hw9T19xO123muqoqn6uo4HAhgdzhwOp3WDoz19fU0NDRQXV1Nc3Mzzc3NlJeXU1ERAOxks/nvfccdtoHvrWR3RSEmqDEbkat8Gj0CHLqWEL+VjGS0vnbtpc97/nnNqlVlpFJeazLTdL20t7dzbulSvt3aSnDXLpbv28e7z57lwQsXOOrz8avKSv6rooLeUMjaI+bEiRPU1NRw9OhRGhoaaGpqAqYDtaTTfl591WnumFRK89JLsGaNdLkIUSqKUVpZB3wU2K+U2jPwvs9rrZ8qwrVvKoXBPlyom/LGmjUKpex4vV5rD5by8nJqa2uZNGmSdSDGhSVL2H/2LC8dO8b07dtZ/cYb/NWZM/xZayuvlJfzn1VVbPR48Pj9dHd309LSQllZ2UCQHwMaOHBgOkuXTgfKAI3LBevX58hmpW1RiFIx6iDXWm/EJJC4ZsON1s2GWLnc4NF6YU09GAxSWVlJXV0dzc3NzJo1i4sLF3J25Up+3tKC/Y03WLBtG7e3tHBvVxfnnE6erKjgV+XlnPN4iEajXLhwgfx/sia+/e293HHHG8BSoImnnvKyZk2+zVGOoBOiNMjKzgnocu2N5u/C/vBUKmUdiNHb20tbWxtnzpyh9fhxqrdsYeX+/Szv7MQGbPH5eCIc5knHOg6f/h7gxGbLsH79/+Lll5PAFJ5/fh7Tpk2juroaj8djnVwEEuhCjDdZol+irmVbXvOxwr1eZszoBNr55jdb6Nq7l5mvv876EydoSCToUR4e0x/nUf6AnSwGvgX8EKjjv/23OmbOnMns2bOZPXs2NTU1+P1+HA6HBLoQ40yC/CZxtdE65LcJcDiSQB9HjvRw8eJFTp48ybEjR/Bv3878rSe4s/MUXhLsZSGPoHiMDroIsmRJ/ui52tpa5syZw9y5c5kzZw5Tp04lGAxauy6CBLoQN5oE+Rgar8U0Vxqt2+0mYHNAgnPn8h0sFy5c4MiRI/z1Q1E+xH5+n1+zgoskgafdbr7ncLAtGKS8qgqHw0FDQwOTJk1i8eLFLFq0iBkzZhAOh60Do0ECXYgbRYJ8DE2UVZGFo3WHNY2tAGVNpJoDMerr24HzQIRFtPG39T/g3va9hLNZzijFj1wuHnM66Rw44MLn8zF9+nQWLVrEwoULmT9/PjU1NTidTivIr2VPdyHE9ZMgH0MTJcgLKaXJb0aZ99prmjVrLpVf7HYNxAEfoHC5snzqI48w+/D3WXPwILdFItiAl2w2/m8ux1MuFxWNjfh8Ppqamli4cCELFy5kyZIlTJs2DafTid1uf0ugm9dm0yZZZCTEaEmQj6GJGFZDB8Zer+a55zSrV1/67/3Vryq+8IX8iN1u13zmMz387u8e49ChQ7Ru3kzziy9y54kTTEqn6QZ+RH7lV0t5OcFQiNraWubOncvq1atZvHgxs2fPtvZGt9lsbNmiWLdODXx/ZMWoEKM05nut3Ko2b770+O67J0ZYFd6TkUopXn1VsW7dpRLMnXeaj+YXAj3wQJCFCxcydepUImvXcu7DH+aJHTvQL73E3C1b+P2ODv5Ea/Z0d/NIdzc/amnh+PHjbNmyhenTp7N69WpWrVrFwoULKS8v56WXnJj9Y1IpePnl8X9thLgZSZCP0ssvX3o8UcKq8J4Mlws2bMg/NqUPs7cLmNE6gMM6Wq6+vp5Zs2bR/c53curUKf7t9dcJ/Od/su7wYf41HuefgV92dvJIZycvHD/O3r17+c1vfsOKFSu4/fbbmTFjLV5vE6lU/gdF/geH1NCFKDYprYzS5s2wdm3+8UQpHxTek3G5ss/Q+v5wh2jkcjni8TixWIyuri727dvHVz78FJ/iZT5ma6M8l6MF+O7AnzN2O1OmTGHx4sVMmfIhEonVvOc9Ie66y2MdjiGTokKMnNTIx1Ap1Mgv95/5ahO1Q/eEyWazJJNJyspiwDke/tdNpH7+c5bt2cOqnh5swPPAo8AvAHdZGYsWLWL9+vWsX7+ehQsXUlFRYdXRJdCFuHYS5GNoYnatDH67GPc2uL1RAxk6O/vp7OykpaWFE6+8gv9nP2P98eM0ZTKDJkh3A7NmzWLFihU88MADrFy5kurqamt/dAl0Ia5OgnwMTeQgH6t7yl9fs3GjZtWqHMlk0tpq9+CBA5z/6U+Z8fLL3NnZiYd8kD9KPthTgQDLly/n/vvv56//+h5gCul0AKfTjqmhT6TfboSYKCTIbzHFCvLhrjPcvMDq1ZdKL4lEgq6uLjo6Oji0aROZH/yA2/bvZ0EySZJ8yeVR8iUYze8CH+NrX6vkM5+ZA4TNdx50HxPph6QQ40XaD0XRDN+pk5/AtNlsOJ1OAoEA9fX1TJ8+ncgDDzBtWguL+ClfbHiSt587xwe1pgUX32UG32EOn/nMJC7976jIL2bS5NsXR27oD6CJ+FuTEMVyff9KxISn9diFlmljhMFtjYWUUrhcLsrLy5kyZQqwln38PfOff57Xf/pTHn372znhnMoX+QqnmMlzvJsP8RPcJMkHeA44AHQC6bF5IkLcJIp1Zud9wDcAO/AfWuuvXOnzpbRSOi43kh1pp475/FxOk0gkiEQiPP10D3//KQcf54d8ku8ylRa6CfMYH+QRHOyhDMgAS2lpWUNVVRU+n++67nukI3IZwYuJaMxq5EopO3AEuBdoBbYDH9JaH7zc10iQl47LBXYxgjGbzeJw9ANpvvZPh0j810+Y/OJm3ssePGTZzWIe5WM8Rg5vwxY+8YnZvPvd72bBggX4/f5BW+pe7b6v9wePTLqKieRyQf6WPuGR/gHWAM8UvP054HNX+prly5drMfFt2mQKNFp7vfm3DfP+a2U+v/AahTKZjO7o6NDwhg6zVf8x/6p3sExr0Alc+se8T9/DHK1AP/jgg/qnP/2pPnPmjE4kEle973//98s/j5E+7xttpK+zuLkBO/QwmVqMGnkjcKbg7daB94kSN9yk5vUYuh/NcHvB2O12qqqqgDlEuI1/409YwQ6WsItv89+5lxd5jjc5CSx84gn+5gMf4K677uJLX/oS27ZtIxqNWmeMDr3vxx8f2fMo1vMupuFeMyGMGzbZqZR6SCm1Qym1o6Oj40Z9WzEK1zKpeS1GFoz5TbZM++FelvAXfJ2pzha+uuyDnHC5+DvgBPBvx45x8stfZv2qVbz//e/n+9//PidOnGDNGjNhmr/v9753ZM+jWM97tK7lB6AQUJwgPws0Fbw9aeB9g2itH9Zar9Bar6iuri7CtxVjrbA2PHQPmZF0xVxrML41qC4F+n8+5+Kj//k13C+9xLc/+1kebWpiBvBjoA1419NP841PfpL77ruPX/7yfwC7gBjPPpvloYcu/zyGc6XnfSNNxN8MxMRUjMlOB/nJzrvJB/h24MNa6zcu9zUy2Vk6ir2w6EqTh1/+Mnz+88N/zHz/XC5HNBqlvb2d3Tt3cvLRR5n52mu8K5HAQz6+HwUeAyLcy7/8ywOsX7+eJUtmAF60Lp2O24m4IZsYndH+exrTlZ1KqXcBXyfffvio1vp/X+nzJchLx1iuEB1quF0bjaFfp7UmnU7T1dXFiRMn2PbMM2Qfe4wNx4+zHEig+AVzeBRN2xw4+OY64J2cPLmcuro63G53SezvIt0zN5fR/veUJfpiXF3rD4TCbNX62r4ul8vR3d3NhQsX2L17N6984yUWbHfxe/yECrppwcV3SPFdYPLALoz33nsvs2fPJhwO43K5JmyoSz/7zaMYv2FdLshL5/dMcUu6llq8zWajsrKSuXPn8t73vpeKt32NP+ebNNDGB3iMwzTwd8Ap4Iuvvsqxf/wJb7/zVT70oX/hkUceYc+ePXR1dZFMJhmPgc2VjOUKXXFjjeWch4zIxYQydER+PfIjn/wXu92aT3/6N3Tt/jemvPwyH0zkmEJqYAXpB/iubSuuVV5WrVrF7bffzuLFi6mqqsLv9+N0OovwjK6fjMZvLmM5IpcgFzelzZvzI5477siybFmStrY2Dh48yD99VeHe5OJTfI8HeQIPSXaT3zP91cZGpt92G6tWrWL58uVMnz6d6upqPB7PuIS6BPnNR2rkQoyC1ppkMskLL/Rz//1lAFSqdj7puYcPxg8OTJDmt9j9ns1G7223sWjJEpYM/GloaKC8vBy3243T6bwhNXUJ8pvPhO5aGSkJcjFetNbYbFkgzQ9+0IZSW3jttdfY8u/P8/sc5yNABfl6+neB5xobKVu0iGXLlrFo0SLmzJlDTU0Nfr8fj+fSGaRjQYJcDCVBLkQBPXAIRjqd5sKFC0ydeg7oZe2yX7Ho5E94sLubuwc+93ngR243R+bNY9q8eSxevJiFCxcybdo0wuEwPp8Pl8tV9IOlJcjFUBLkQgxDa82mTZrbb8+nptud44/+6HE6Op7k4o4drD95kg+nUkwBusgfVfd0QwPpBQuYN28eS5Ysobm5mebmZoLBYFFDXXrIxVByQpAQw1BK8eqr5kQiRSZjIxC4n69/fR6wk3l/uZeH9u2j9o03uP/CBf5Aa/6ftjZ2tbXxo1df5TuzZlE1cyYLFy5kzpw5NDc3U1dXh9/vtzpfbDabdcD0tQb70H1WZFWnuBIJcnHLy+/9kg9Ylwvuu8/DP/7jXGAKf/zH6zhw4AB79uzh8YMH+b/HjrH6xAke7OnhnxMJEvv28eSBAzy1ZQsbZ89m5uzZzJ49m6RzGfkAACAASURBVOnTpzN58mTKysrw+/3WJKndbr+mUB/+OL3RP1cp19ycpLQiBJfaFTdsMAdR5M8MzWQ0mUyGixcvcvLkSbZv386bb77JsWPHCBw9ygMdHTyYSFABtCjFTzweNs6YgWP6dObOncvs2bOZPHkydXV1BINBPB4Pbrcbl8tljdSBtwT7WO2zIuWa0iY1ciFGoPBoOrN5fzabJRaLcfToUQ4ePMiRI0c4efIkbSdOsOTUKd5z8SJ3DeyJ/prLxa9rajgwYwa1kyczb948mpubaWpqoqamBp/Ph8/nw+1243A4Bo3SL/3J30OxQlc24Sp9UiMX4jps2aJYsyafqEopQqEQS5cuZe7cubS2tnLgwAEOHz7MyZMn+eqZM3yttZW7z5zhd2Mx/qm1le6zZ/lVIMCz27bx+syZzJgxg2nTptHc3Ex9fT3hcBi/34/X67UmSU2om3LP6tWXHo/GWJVrxPiTIBdiiMtNNBaWQXw+H7NmzaK5uZlly5Zx/Phxjh07xpEjR9jV0cEvT51i1pkzPHDxIh+MRvnEoUMcOHaMJ6ureWbaNAJNTcyYMYPJkyfT3NxMdXU14XAYj8eD1+vF7XaT30xUkctd+kEymk6YiXJghig+CXIhhrjayNUEei6Xw+PxMGXKFBobG5k3bx4tLS0cPnyY1tZWzp49y7dOneJbXV2sPXWKByMRPt/WRqKtjedDIf6rvp7d06czqbmZadOm0djYSHNzM3fcUQn4AA/gI53OT5IWhvjQx9diohyYIYpPglyIIa515GoCXWuN0+mkvr6e6upqZs2axenTpzly5Ajnzp3j/PnzHG5t5S8uXKB8oJb+nmiU+3t7aT1+nKdqanh1+nT2NjczefJkoA5oBhqAMi5ezNfSvV4vDocDh8NhfV/jWkfr0q1yc5LJTiGGcT0TjVpr6wDobDZLX18f586d4/Tp0xw/fpz29na6u7s5duwY6WiU5WfO8J6ODtbF4+SAbcEgzzU3849vTCLFZL70pUZrtB8OhwmFQlbZxev1Wu2MQ8P7ekbrojRI14oQIzCafuvCQM/lcvT399Pe3k5LSwtnz57lwoULdHR0cPr0abq7uwl2dvK206f57UiExkyGTuw8RjVHb59BR2MjdXV1NDU10dDQYIW66U0vDPXChUfXM1oXE9+YdK0opf4JeABIAceBT2qtI6O5phAT3dVCXimF3W5Ha41SikAggN/vp76+ns7OTs6ePcuZM2dobm6ms7OTzs5Ofl1Twy9TKWa1tnL7kfP8IRdwbzzPYb+fZyZNYsvUqRyoq7NCvb6+nrq6OsrKyqxRuulPdzqdVvdL/j61FezXMlqXRUOlZ1QjcqXU24EXtdYZpdRXAbTW/+NqXycjcjHRXSnMRhp0JkjNv7V4PE5PTw9tbW2sXHkeuAB08ulPt9Pe3s4Pf9hJOd18rrmN+9vbmZtIkFCKjdXVvDJ9Oi3TplFZXU1VVRUNDQ3U1tZSXV1NMBi0Rufm78LVpOYHy5VG68U42EOMnTEZkWutny14cwvwvtFcT4iJopghZoLShLnX68Xr9VJZWQk0Ap3ARZYta6Orq4sf/tBDN3X8qukVXl5wiKozZ3j3+fPc19nJPe3tnNu9mxeamtg4YwZnmpsJBALU1tbS1NREbW0t5eXleL1e/H4/DofD2m7XlF/sdnvB8xw8Ws/3q0sJptQUs2vlU8BPi3g9IW4qQwM9f+pQLVANJFGqA7c7CswDbGzbdh/vfe+3UMt28eOeHr4TibCqrY13njvH7x09yoePHmVXZSUvTZnCm3PmcPr0aUKhEI2NjTQ2NlJRUUE4HMbr9eLz+axQt9vtVumlMNTzcpggf/11WLtWauul4KpBrpR6nnw/1FBf0Fr/auBzvgBkgMeucJ2HgIcAmpubr+tmhZhINm++7uO6Ckocmvw/Qwd/8Ad+PvCBJGYhUDYLWt/JbbfZuXDhAu3t7RyorGTHzJmEurq4p7WVd7W389c7d9K7bx+vTJrEazNnsqejg5MnTxIOh6mtraWxsZHKykoCgQBer9caobvd7kGhvm2bncLz2O++G55/Pjew94xMmE5kVw1yrfU9V/q4UuoTwP3A3foKBXet9cPAw5CvkY/sNoWYGIq5veyWLYPLGMmkxuVyDbylrZ0YGxtnUV1dTSQS4eLFi0QiEbq7u/lVdTWPJ5PMO3+ed7S28o5Tp3jg5EmOhUK8OHUqe+bOpb29nVOnTlFRUWEFenl5OYFAwJogNaH+/PMewGndUyqlefVVG2vXDq7xgwT7RDParpX7gL8B7tRa9xfnloSYuIq5X0nhtQDsdsXHPqZ5+OF8eeOZZ7IsWlRLNOqjs7OTWCxGX18f58+fp6enh/b29ny419fzzSlT+HZPD3eePcs729p4aO9eUvv3s6WhgU1z5nB02jTa29vxer1UVVUxadIkKioqCAQC+Hw+PB4Py5b5gUryvyXkF0PdcUcWGNzWOHTyVkJ9/I22Rv5NwA08N/AfcYvW+r+P+q6EmKCKuV/J0K/95jdh3TpT2sgNPM6fOlReXk5PTw+RSIRAIEA0GqWmpoZoNMr58+eJRCJEIhFeqazk2VmzmNzVxT1nznD3+fOsb22l1VnFC03zeeO2Ojprajh37hx+v5/q6mrq6+spLy9nypQgEAfCgI0nnkixeLGdeHxwTV0pZXXBDA11kGAfD7IgSIgRKub2ssNdq7C9sXBxkdaaZDJJb28v0WiUvr4++vr66OzstEK+o6ODWCxGb29vPtzPNVG3ZSmf4gfcw4vkgH01Nbw+ezZH58/HU1aGx+Ox6ukf/nANEAL8HD/usXrSnU4nLpfrLROlhYEtoT72ZGWnEEVSzAUzw11ruPcNDfRUKkVfXx/RaJRYLGaVXUwN/eLFi8TjcV544Ta2bXsPYGcKR/l0+Wf4QPw5ahIJoi4XW2fMYPvChXQ0NBAIBPB4PHzhC2Gghpdfrh5UejGhbhYd2e1268/QU48uF+rmbwn26yP7kQtRwgpXi+ZyOWsVZygUIhqNEo1GCQQCBAIBUqkU3d3ddHd309eXYNu2LABn7JN5YvECXgzCwvZ27jlzhjvffJN7Dh7kREUFm+fMYe/8+UAZcIFt2wJUVlZSUVFBbW0twWAQp9NpHVtn/pjTjkygm1A3K0uBCVVbvxlPSZIRuRAjNBGWsA/dz8Xs6dLX10dPTw/pdNp6vG2bnd27y2hoOEJFxWGrJNPb24s/kWBdSwv3njnDtN5eUnY7e6ZNY9v8+ZyZMQOn222NyisrK6mtraWiosIaobtcLuuP2ZmxMNCHO590PEfrpX5KkozIhSiSibB0feh+LjabzQrcYDBIX1+f9bimJsUdd/TQ1eUlGp1OJBIhGo2SSqXo6upiY0UFz8+ZQ9PFi7yjtZXbT51i5dGjdAQC7FiwgF2LFtFRU0MkEuHcuXMEAgHq6uqoqKggFArhcrmsLXYLR+qmlj401GHwFsDm7xsxWr9ZT0mSEbkQN4HCIDQj9HQ6bdXQk8kk8Xic/v5+ent76erqIhKJDHpfb28v/f392FIpVpw5w9tOnWJxRwc54FBTE9vmz+fEwoUw0HseDocpLy+nvLyccDhMWVkZDocDn89n1dAvF+pDzycd7nkYxRyt36wjcglyIW4iwwV6Npu16uiZTIb+/n6rjh6NRunp6aG3t5dYLGZt6PWtb/UD/fzJu2JsaGlhw8mTVPX30+d2s3v+fHYvXUrXpEnWsXfBYJBwOGztm24WHJmyy/WGeuFzMkY7Wt+8OT8S37ChtEIcJMiFuKUMF+i5XI5oNEoikSAWi5HJZKwReUdHB6lUiosXL5JIJPjsZxNAHw89lG91VLkc88+f566TJ1nR2oozl+NMTQ3bFy3i4OLF6LIybDYbXq+XQCBARUWFtddL4fa6ZmJ0pKE+9DkZt9ohGlIjF+IWMjTgbDYbuVyOcDhMNpvF7/cTj8dxOp34fD5CoRCJRIKysjKi0SjQD/QwfXo13d3d9Pb2ciIQ4I3GRkLpNGtPneLO48d58Pnnuf+llzgwaxZ7ly+ndcYM4vG41dPu9XoHlV5MoBeGeuFE6ZVCfbj2RsCa9C183rdCqBeSIBfiJmcmFs0Se6UUoVAIn89nBXo8HsftduP3+0kkEkAv0ENNTR9er5eKigoSiYS1GOl5r5eX5s+n6eJF1h89yprjx1l26BCdZWXsXryYQ6tW0as1vb29dHd34/f7rfJLeXm5tcWuy+UimUxaG3iZ2rrD4bB+AF0t1OFSsNvt+Y3IMplba7QupRUhbjG5XG5Q+2ImkyGVSlmBHo/HmT9fA3FefjlmTYr29PTQ399PNpslEolYi5DS6TSubJaVZ8+y7sgR5p49Sw44NnUqu5Yu5czSpSiPh2w2a5VeysrKKCsrIxwO4/P5BtXQC0fqZsHRlUK9UP5NTTZ7qb5+6WOlP1qXGrkQYpChgZ7NZq1Ar6vLAQ5+9KMYc+dGrL7zRCJh9aGnUilisZi1XUAqlUJrTW1/P6vffJN1R45Q3tdHzOtl38KFHFi1ikhTE4DV3eL3+ykrK6O8vJxgMDhosVFhX3rhoRhXHqnnn9vQVbE3y9YBEuRCiGEVTohCvqvjzjvzVVePR/P977fx/vf3ASm2bIkTjUbp7e0lHo8PmjSNRCL09/fT399PJpPBabMxt62NtYcOsfjUKRzZLGcbGti7bBnHbruNTCBgLR7yer0Eg0HKysoIhUIEg0Grdm7+djqdKKWuGOoOR35rYK2HD+jCQC/FYJfJTiHEsMziIjMh+uqr1kdIp2Hv3nLyW9s6CYWc+P1+KioqrB71vr4++vv78Xg85HI54vE4XV1d9Pf3c6ChgYOTJhFKp1l59ChrDh3iXb/+Nemnn+bNefM4sGoV52bNIpfLEYvFuHjxojVKNyUYc7i0GaGbEbvNZnvLBl75QzlsbNx46XSjwnC+Ut/6eG8dMBoyIhdCDLJpk2bduvxjrxeefjrDnXdmgDhnzyZIp9Mkk0mrFBOLxejp6bEex+NxUqmUtdAomUwOdMKA0+Gg8cIFVh88yIojR/Amk3SVl7Nv2TKOrFlDvKrKCmWn02lNkgaDQUKhkFV6GbqK1GazsWuXm3vv9Vn3/cwz2becbnSlcB6uBAMTK9iltCKEuGZm0nDjRs3KlVlcrnxOxOP5PV0SiXygZ7NZ0uk06XSaeDxOX18fkUjEqrVns1mSySQ9PT309fURi8VIp9MopfDbbCw4epTVBw8y8/RpNHBi1iwOrlpFy+LFOPx+crmcNRI3o/NAIDConm66XL797TBf/nIQUNjtmr/92xSf/eylLQyGC/SRBvt4h7oEuRDimhVOGmqtsdnybX2vvJI/w9MEdDweJ5lMksvlyGQyJJNJUqmUtW96X18fiUSCZDJJOp229kpPp9P09/eTy+Ww2+1U9/WxYv9+bnvjDcLRKP0+HweXLePwunVEmpoG7X3uHtjIKxQKEQqF8Pv9uFwu9u718f73NwAKjwd+8YteVq/WVumlsKY+NMxLZbQuQS6EuGaFQX5pfxJtlSxWrcphs9kGtS6mUinS6TSAFeim1GImRuPxOIDV+RKNRonH42QymfzZoQ4HM1taWLF/P/OOHMGRzXKusZEDq1ZxYuVK0gOhbQK5cJLU7/ezYkUlEODxx9MsX54atOeL6U0v3D/9ekI9/7qMz2h9TCc7lVKfAf4ZqNZaXyzGNYUQE8OlHQMVqRS89pqddetsZLNZK0zdbrdVXkmlUtjtdmtXxMKyiymtFC7lN9vvptNpenp6eGPSJI5MnYo/kWDxgQOs3L+fe594gvSTT3J08WLeXLOGC3Pm4PJ4rDJPV1cXfr8faAdCTJ8eIJHwkEqlcLvdpFKpQa2M5pDrwiAfSagPt8p06ITp0M8bS6MOcqVUE/B24PTob0cIMZEolT+AwXC54K67FDabslaKZrNZq+RhQtMEejabtWrYZrm+mRz1er3E43GCwSDl5eXW6tJEIpHf1Mtm4/Xly9m2ahUN58+zYu9eFh04wLydO4lUVHBozRqOrVtHvLKSXC5HT08PEAU6OX7cZ02Umg28TCujKf0MHakXI9RhfDphRl1aUUr9HPgH4FfAimsZkUtpRYiJrTBvtL76qTpml0WzBQBgjcRNDT2Xy1nhburmsViM/v5+kskkmUzG2p3RLD5KJBL09/fjcrnQWuMF5h0+zLI9e5h28iQ5pWiZNYs3163jzJIl/OGf+QA73/2uIpfL4XQ6ra6X8vJyazdGs3LUtDCaHzajDfWhhivBFO7LPvL/LmNQI1dK/RbwNq31nyulTiFBLsRNozC8zR7eV4uLoYFuFhuZOrlZTWo6XQpH7729vYO6YBKJBPF43Fp8lEgkyGQy1gKhcFcXy/btY+nevYR6eoj7fBxZuZI3166lZ/JkK5hzuRwejweXy2UdvuH3+6332e12gEGll7EKdfM11+u6g1wp9TxQN8yHvgB8Hni71rrnakGulHoIeAigubl5eUtLy8iegRDihhl6AMPAHOU1n4403Ai9MNBNwGcyGasDxuybXjhCT6fTaK3p6emhp6fHKruYco7WGodSzG5tZfGOHcw+dAh7Nsv5piYOr1vHmXXrSPv9VkBns1mcTiderxev12vtpe7z+awVosBb9k8fruNlNKF+vYo+IldKLQReIL/fJcAkoA1YqbU+f6WvlRG5EBPbl78Mn/98/rHdDtn8+c0jPuZOa00mk3lLoJvWxUwmg1LKCmfzdzKZJBKJWCNzyHfCJJNJ+vv7icVixGIxa38Xs0DIHY2yaP9+lu7eTc3582ScTk4sWcLxDRtonzsXXRDAuVzOGqWb7hefz4fH4xnU2WIWHE2EUB/z9kMprQhx8yj2kWhDA10pNagX3ZRcTGnFtDJms9lBuyzabDar5FLY7RKLxayJVa01ToeD5o4OFu7Ywbw9e/AkEvRUVvLm2rUcv/120rW1QL6ckt/+1o7b7cbn81mjdJfLZW2va4LdhLh5/3BtjIXbBg+3CGk0JMiFECNytQnO6zE00M3+LmZiM5PJAPkFR2aEbhYaJZNJ+vr6rJJL4alHhbsw2u12crmcNYq2JZPMP3yYhTt2MPnYMXJK0Tp3LkfvuINzt91GbshI20x8mpF6YdfL0P1dzPe4XKhfKdivhywIEkKMSOGioGIrDHTACt9kMmkFuun2MOUWM0KPx/M7MKZSKWtkb7pbzGjdjNiBQXXvmr4+5m/bxoKdOwlGIsT9fo6tXs3Ju+6ib+pUK2TNfZme82AwiN/vx+/3v2WUXrjIyAT+0BJLYbg7HNff9S1BLoSYcEw5xWyha0ojZgRuwtoEfyqVsuri6XTaGo2br7fZbNaEqdnIKx6PW9fNbzdgg2yWmS0tLNy+nen792PPZmmfMoVjd9zBqTVryAycYGSuaRY+OZ1OPB4Pd90VBuwcPeq12hjN55qQNx0xQ0PdjOCvhwS5EGLCGi7Q4dJSfzPhaSYpk8kkU6emgTQHDui3TH6asktfX5+1qtSM1k3YQj54g8kk83bvZv7WrVSfO0fG6eTksmUc37CBrgULcAwEtSnXOBwOfud3nICXjRvd1ijdTJKawDaBbiZizajc6/UWvbQi+5ELIcadKTmYQM8OtMmYCUePx2PV0e12Oz6fD8gAKdzujLXlbSKRIBqNWpOnLpeLdDpNIBCwDr0wPxzMCL3X5WL72rXsWLuWxvPnWbBtG7N27mTm1q30VlVxbP16Tt15J9naWuu3BQgBNbz2WicLFrRadXSfz2dtTWB+izD3bLpeTB98UV8/GZELISaaoSN0M5o1wZhMJqmtzQE2zp3LWYuITGeLWUiUTCYHdcoULjAyo3Uz0i6sY3tyOabt28fCbduYdORIfoJ0/nxObNjAa+H7+Lt/fAdgw+XK8cUvvsKsWZ34fD6r9OLxePD5fLjd7reUXmpqaqy9XkZKRuRCiJJROCmYyWQGlVxMNwmkAc3OnSlWrry04CeZTOJyuaiqqiKTyRCNRq39003Ymm4XcwapWYBkNgLry+U4tHQpR1asoLy7m3lbtzJn61aav/lNlrp+SICP8wh/wKHMPN58s5Z58y4tVirssDF96h6Px5okNc+lqK+XjMiFEKWgMNC3brWxfn2+zu31wuOP97B0acJakm9WkZrVo1prK7jNBKrdbrdWk5pOGLPgyNTCTcdJOp3GZbcz+cgRmp/fzdwjO3CRZrtaQdt9q8n+7jIoK7Nq6eYHkamTmyAPBoNMnjyZQCBwXa+BjMiFECWtcIT+yivmvYpUSrNzZ4D1653WpKjNZiMUClmlmGw2SzgcRmtt1dFNoAeDQSorKwcdTdfT02OVZUy5JpnJcHzmTE7NmcPdf/YQH2Efn616gtv+65tknndyasUKTt51F5HFi3G53VbLpFKKWCyGz+cjFovR0NBQ/Nem6FcUQogx5HA4uPtu+OIX86cWuVxw551YZROzrN/s9RIIBKydFZPJpLXPimlfNKtLTXthKpUiGAxao3Qzii+soXfSx78wmVl/91fUt7UxZ9MmZmzbxozNm+mtqeHY7bdz5m1vI1dbi9PpRCk1qD++2KS0IoQoSaaD77XXsqxcme9yMa2FplZuJkDNMnyzla5ZMWpaFc1y/2w2O2hBj9keoL+/31pkZEonptXQnCvq1Zppe/Ywe9Mm6g8dIqcUbQsWcOruu7m4Zg12rxe73c7tt99OKBS6zucsfeRCiJvI0JWnhW2LJtBNz3k6nSaZTFoBbereZpKzcLvd/v7+QQuRnE6ntVGXCfz+/n6rG8as6jTnjyqlKO/uZvamTcx8/XX83d0kgkFOrl3LibvuYsUnPkHtwF4vI3/OEuRCiJvI5bYQKAx0M+lYOBI3I3QzkZnJZKyJUdO1YtoUTci7B2repiOl8Oi6VCplfS9TfjEje4/TScMbbzD79ddp2rULezZL9DvfIfiJT1znc5YgF0LcQkwfugltE+hmbxazKKhwW1uzta4J9EwmYx0UDVhL7M1I32wH0Nvba/2AAKyRubm+1hp/PM7s7dtp/tKXqJs587qek3StCCFuKab1zwS66SAxe514PB5rhG0WBeUPcMZaRepyuSgrKyMQCFg954U199raWhKJBJWVlUSjUXp6eqzNvcwh1JAv5fS6XOxav57mYLDoz1WCXAhxU7tcoJs+b7fbbZ0vagLdrM405RhzxqdZ6m/2P9daD1rFWVZW9pYzR81kqKmljwUJciHELWFooJvTicz5nKZ10Rw7Z8LbBL3pUff7/QQCAeucUdNvbvaE8fv9lJeXE4lESCQSxGIx63g7U08vNglyIcQt5UqB7nK5rEA3f8zHgsGgtR+6OY2osrKSTCZj7b5onSPqcFBbW2tdw5RdzOcUmwS5EOKWZEodZoOuoYFujpwrPATaZrNRVlZmlVxSqRROp5NQKEQoFLIC3eysaEb0DoeDcDhMf38/Ho+n6M9l1EGulPpT4E+ALPAbrfXfjPquhBDiBig8em1ooJvtcgvPETWHWyil8Pv9eL1ea0MugFAoRDAYtJb5myPpXC4XSqkx2cIWRhnkSqm7gN8CFmutk0qpmuLclhBC3DhXC3SPx2ON0k0N3Xzc1MXN5CaAx+PB6/WSSqWsUbppcxwLo/3R8EfAV7TWSQCtdfvob0kIIcbH0EAvrKGbDpfCGrqZvDSdKeFw2DrRKJPJWO8rPNjCHDBdTKPthZkF3KGU2qqUekUpdVsxbkoIIcaTWaVpjmozh1qYUDenAZmTgcwKUBPefr+fYDBo1cNNoNfV1Y1PjVwp9TxQN8yHvjDw9RXAauA24P9TSk3TwywXVUo9BDwE0NzcPJp7FkKIG6JwhG72UzFhbtoTh2tdNJOmTqcTt9ttrRg1K0iL7apBrrW+53IfU0r9EfDEQHBvU0rlgCqgY5jrPAw8DPkl+td9x0IIcYMNF+iFy//NeZwm0AtH72bP81wuRyKRwOl0Fv3+Rlta+SVwF4BSahbgAi6O9qaEEGIiKiy5mOX3Zv8Ws1Oiz+fD6/Xi8XisfnWzbN/r9U7IPvJHgUeVUgeAFPDx4coqQghxMxluUnToCN3pdFoj81QqZXW6mO1ui2lUQa61TgG/V6R7EUKIknK1QLfb7VbZxSwsGov9VmRlpxBCjNLlAt28bQLd7LxYbBLkQghRJEMD3YS6y6UATTY7Nrsfjs1VhRDiFmYmRc3EKOQnOF9/PTcmhy9LkAshxBhRSrF1q4181Np4xzvsbN5c/O8jQS6EEGPo5ZchPyJXpFLwyivFbz+UIBdCiDG0YcOlxy6XGvR2sUiQCyHEGFqz5tLjF14Y/HaxSNeKEEKMsbFeJikjciGEKHES5EIIUeIkyIUQosRJkAshRImTIBdCiBInQS6EECVOglwIIUqcBLkQQpQ4CXIhhChxajxOZlNKdQAtN/wbj68q5DxTkNfBkNchT16Hkb0Gk7XW1UPfOS5BfitSSu3QWq8Y7/sYb/I65MnrkCevQ3FeAymtCCFEiZMgF0KIEidBfuM8PN43MEHI65Anr0OevA5FeA2kRi6EECVORuRCCFHiJMiFEKLESZCPMaVUk1LqJaXUQaXUG0qpPx/vexovSim7Umq3UurX430v40UpFVZK/Vwp9aZS6pBSagwO/pr4lFJ/OfDv4YBS6sdKKc9439ONoJR6VCnVrpQ6UPC+CqXUc0qpowN/l4/0uhLkYy8DfEZrPQ9YDfyJUmreON/TePlz4NB438Q4+wbwtNZ6DrCYW/D1UEo1An8GrNBaLwDswAfH965umO8C9w1532eBF7TWM4EXBt4eEQnyMaa1Pqe13jXwOEr+H27j+N7VjaeUmgS8G/iP8b6X8aKUKgPWA48AaK1TWuvI+N7VuHEAXqWUA/ABbeN8PzeE1vpVoGvIu38L+N7A4+8Bvz3S60qQ30BKqSnAUmDr+N7JuPg68DdAKULuOQAAAaNJREFUbrxvZBxNBTqA7wyUmP5DKeUf75u60bTWZ4F/Bk4D54AerfWz43tX46pWa31u4PF5oHakF5Agv0GUUgHgceAvtNa9430/N5JS6n6gXWu9c7zvZZw5gGXAt7TWS4EY1/FrdKkbqAH/FvkfbA2AXyn1e+N7VxODzveDj7gnXIL8BlBKOcmH+GNa6yfG+37GwTrgPUqpU8BPgLcppX44vrc0LlqBVq21+Y3s5+SD/VZzD3BSa92htU4DTwBrx/mextMFpVQ9wMDf7SO9gAT5GFNKKfI10UNa6/8z3vczHrTWn9NaT9JaTyE/qfWi1vqWG4Fprc8DZ5RSswfedTdwcBxvabycBlYrpXwD/z7u5hac9C3wJPDxgccfB3410gtIkI+9dcBHyY9C9wz8edd435QYN38KPKaU2gcsAf7fcb6fG27gN5KfA7uA/eRz6JZYqq+U+jGwGZitlGpVSv0+8BXgXqXUUfK/rXxlxNeVJfpCCFHaZEQuhBAlToJcCCFKnAS5EEKUOAlyIYQocRLkQghR4iTIhRCixEmQCyFEifv/AeOF07gOzo1EAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "objective.plot(samples=300);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Fitting a neutron reflectometry dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We start off with all the relevant imports we'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import os.path\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "\n",
    "import refnx\n",
    "from refnx.dataset import ReflectDataset, Data1D\n",
    "from refnx.analysis import Transform, CurveFitter, Objective, Model, Parameter\n",
    "from refnx.reflect import SLD, Slab, ReflectModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "It's important to note down the versions of the software that you're using, in order for the analysis to be reproducible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "refnx: 0.1.18.dev0+0fbfe9f\n",
      "scipy: 1.5.3\n",
      "numpy: 1.19.4\n"
     ]
    }
   ],
   "source": [
    "print(f\"refnx: {refnx.version.version}\\n\"\n",
    "      f\"scipy: {scipy.version.version}\\n\"\n",
    "      f\"numpy: {np.version.version}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Loading/Creating a dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "*refnx* reads 2, 3, or 4 column plain-text files using `Data1D` or `ReflectDataset`.\n",
    "\n",
    "\n",
    "| columns |   data                  |\n",
    "|---------|-------------------------|\n",
    "|  2      |         $x, y$          |\n",
    "|  3      |    $x, y, y_{err}$      |\n",
    "|  4      |$x, y, y_{err}, x_{err}$ |\n",
    "\n",
    "\n",
    "$y_{err}$ being the standard deviation of the measured $y$ data, $x_{err}$ being the uncertainty in $x$.\n",
    "\n",
    "In a reflectometry context $x$ is the momentum transfer $Q$ ($A^{-1}$), $y_{err}$ is the uncertainty in the reflectivity, and $x_{err}$ is the full width at half maximum (FWHM) of the Gaussian approximation to the resolution function, $dQ$.\n",
    "\n",
    "The dataset we're going to use as an example is distributed with every install. The following cell determines its location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "pth = os.path.dirname(refnx.__file__)\n",
    "DATASET_NAME = 'c_PLP0011859_q.txt'\n",
    "file_path = os.path.join(pth, 'analysis', 'test', DATASET_NAME)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`ReflectDataset` uses a file path to load the data. However, you can also make a dataset directly from numerical arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "data = ReflectDataset(file_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Creating an interfacial Structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`Structure` objects describe the interface of interest. They are made by assembling a series of `Components` (the simplest `Component` being a `Slab`. However, the first step is to create `SLD` objects that represent each of the materials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "si = SLD(2.07, name='Si')\n",
    "sio2 = SLD(3.47, name='SiO2')\n",
    "film = SLD(2.0, name='film')\n",
    "d2o = SLD(6.36, name='d2o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`Slab`s are created from these `SLD`s to represent each layer in the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# first number is thickness, second number is roughness\n",
    "# a native oxide layer\n",
    "sio2_layer = sio2(30, 3)\n",
    "\n",
    "# the film of interest\n",
    "film_layer = film(250, 3)\n",
    "\n",
    "# layer for the solvent\n",
    "d2o_layer = d2o(0, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "A `Slab` has the following parameters, which are all accessible as attributes:\n",
    "\n",
    " - `Slab.thick`\n",
    " - `Slab.sld.real`\n",
    " - `Slab.sld.imag`\n",
    " - `Slab.rough`\n",
    " - `Slab.vfsolv`\n",
    " \n",
    "We need to specify which parameters are going to vary in a fit, and what the limits are on those parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sio2_layer.thick.setp(bounds=(15, 50), vary=True)\n",
    "sio2_layer.rough.setp(bounds=(1, 15), vary=True)\n",
    "\n",
    "film_layer.thick.setp(bounds=(200, 300), vary=True)\n",
    "film_layer.sld.real.setp(bounds=(0.1, 3), vary=True)\n",
    "film_layer.rough.setp(bounds=(1, 15), vary=True)\n",
    "\n",
    "d2o_layer.rough.setp(vary=True, bounds=(1, 15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now we assemble the `Structure` from the `Components`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "structure = si | sio2_layer | film_layer | d2o_layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________________________________________________\n",
      "Parameters:     'SiO2'     \n",
      "<Parameter:'SiO2 - thick' , value=30          , bounds=[15.0, 50.0]>\n",
      "<Parameter: 'SiO2 - sld'  , value=3.47 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter: 'SiO2 - isld' , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'SiO2 - rough' , value=3          , bounds=[1.0, 15.0]>\n",
      "<Parameter:'SiO2 - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n"
     ]
    }
   ],
   "source": [
    "print(sio2_layer.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`Structure` has a `sld_profile` method to return the SLD profile. Let's also plot that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAELCAYAAAA1AlaNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAewklEQVR4nO3de5QcZ3nn8e8z96tmJM1YlqzLWLaxsA2+yeAbF5sQbCBAElgIi7kGEWBJAkuCWdgs5IQcWBZ27RyCMbcAIVxiMHAIeE0cgw0sBknItmTJsrEl666RNCPNpaevz/5R1ePReEbqlru7qrp/n3P6TE9VTdVTqtHT7zz11vuauyMiIo2hKeoARESkdpT0RUQaiJK+iEgDUdIXEWkgSvoiIg2kJeoATmZgYMCHhoaiDkNEJDE2bNhwyN0H51oX+6Q/NDTE+vXrow5DRCQxzGznfOtU3hERaSBK+iIiDURJX0SkgSjpi4g0ECV9EZEGoqQvItJAlPRFRBpI7Pvpi4jUm3Quz4adIzxyYJzJTJ50Lk+hcPww913tLfzZC86q+LGV9EVEamjznqO842sb2DOaOm652fHbDfS0K+mLiCTZ0cksb/2n39Da3MTnbriUi1f009vRSltLE81NdvIdVICSvohIjXzpF49zcCzND99zNRec0RdJDLqRKyJSA4WC883fPMELnjEYWcIHJX0RkZp4aN8xDhxL84oLl0Uah5K+iEgN/OLRQwBcfc5ApHEo6YuI1MCmXaOsWtzFkgUdkcahpC8iUgPb9o9x3tIFUYehpC8iUm2TmRw7Dk+w5nQlfRGRuve7gxO4w7mn90QdipK+iEi1PXFkEoCVi7ojjkRJX0Sk6naPBEl/+aLOiCNR0hcRqbpdI5P0dbayoKM16lBqn/TNrN/MbjOzbWa21cyuqHUMIiK1tOtIihUxaOVDNGPv3ATc4e6vNrM2oCuCGEREambf0RRDi6Ov50ONW/pm1gc8H/gigLtn3H20ljGIiNTa8Fia0xa0Rx0GUPvyzpnAMPBlM/utmX3BzOLx8SciUgWZXIGRySyDPdE+iVtU66TfAlwCfNbdLwYmgBtnb2Rm68xsvZmtHx4ernGIIiKVc3giDcBgb2O29HcDu939vvD72wg+BI7j7re6+1p3Xzs4OFjTAEVEKunQWAZo0KTv7vuBXWZ2brjoRcBDtYxBRKSWhsengPgk/Sh677wH+HrYc+cx4C0RxCAiUhPDY0F5Z6CnLeJIAjVP+u6+CVhb6+OKiEThyEQWgMXd8Wjp64lcEZEqGk1laGtpoqM1Huk2HlGIiNSp0Yks/Z2tmFnUoQBK+iIiVTWayrCwKx71fFDSFxGpqpHJLH1d0Q+0VqSkLyJSRUcng/JOXCjpi4hU0cikyjsiIg3B3RlNZelXeUdEpP6lcwUyuQILVN4REal/4+kcAD3tUQx+MDclfRGRKhmfUtIXEWkY0y39DiV9EZG6p/KOiEgDmVDSFxFpHMWWfreSvohI/Ssm/V7V9EVE6t+EWvoiIo1jfCqHGXS1NkcdyjQlfRGRKhlP5+lua6GpKR5j6YOSvohI1Yyns7HquQNK+iIiVTORztPdHp/SDijpi4hUzVg6R09HfAZbAyV9EZGqmUjn6FFLX0SkMYxP5VTTFxFpFOPpXKz66IOSvohI1Yync/Qq6YuI1D93Z0ItfRGRxpDOFcgVPFZj6YOSvohIVUwPtqaWvohI/StOlajyjohIA4jjrFmgpC8iUhVxHFYZlPRFRKpiKlcAoLNNT+SKiNS9VCYPQGeMxtIHJX0RkaqYyirpi4g0jFSY9DuU9EVE6p/KOyIiDWS6pd8WrzRb875EZrYDGAPyQM7d19Y6BhGRapvK5mkyaGtu8KQfusbdD0V0bBGRqktl8nS2NmMWn0nRQeUdEZGqSGXzseujD9EkfQfuNLMNZrZurg3MbJ2ZrTez9cPDwzUOT0Tk6Utl87HruQPRJP2r3f0S4Hrg3Wb2/NkbuPut7r7W3dcODg7WPkIRkadpKpuPXc8diCDpu/ue8OtB4HbgObWOQUSk2lIZlXcws24z6y2+B34f2FzLGEREaiGu5Z1a995ZAtwe3s1uAf7F3e+ocQwiIlWXyhbo72yNOoynOGnSN7MXA/8J+Iy7bzKzde5+66kczN0fAy48lZ8VEUmSqUyezgUdUYfxFKW09N8KvBP4sJktAi6qbkgiIsk3lUtuTX/M3Ufd/f0ENfjLqhyTiEjipTJ5Olrj9yhUKRH9W/GNu98IfLV64YiI1Ie43sg9adJ39+8DmNlA+P0/VDsoEZGkq4d++l+qWhQiInUkmy+QzXvik368Rg0SEYmp6VmzEnojt8irFoWISB2J66xZoJa+iEjFTWUKQPxmzYLykv4HqxaFiEgdSdVDecfdNUaOiEgJppN+wlv68zKz/krsR0SkHhQnRY9jTb+UsXcuBV4O3AzkgPNnvS4AugElfhER4t17p5Sxdz4HvAN4gmBC8y3ANmAr8DrgonBsfBERId7lnVKS/i+BvwI2Al3A59392wBm9ldK+CIix3uyvBO/sXdOmvTd/c/NrMvdJ8NRNj9sZu8F/hb13RcReYqpXHxb+iV9DLn7ZPj1iLu/j6Cs83pgiZldU8X4REQSZ7qlH8Oa/in97eHuO939BuAq4ANm9rPKhiUiklxTCa/pz8vdNwHXqbUvIvKkVDZPS5PR2pzAmr6Z3QDsBt5F0GXzHnf/7Mxt3P3u6oQnIpI8qUwhlq18KK2lfxlwobu/BsDMbqpuSCIiyZbK5mNZz4fSkv4xYLmZvR0YIXgQS0RE5hHXCVSgtBu5/x34HrAIaAPeU9WIREQSLpWJb9IvpZ++A98zswF3P1SDmEREEi3O5R1NlygiUmGpbJ7OGD6NC5pERUSk4pJe0y/SkAsiIiVIZfKxHFYZ1NIXEam4VJ209DVdoohICaZifCO35GEY3H2zma0BXgmcES7eA/zA3bdWIzgRkSSaysb3idySW/pm9gHgmwRlnl+HLwO+YWY3Vic8EZFkcfdYl3fKGXDtbcD57p6dudDMPk0wm9bHKxmYiEgSZfNOvuCxnCoRyqvpF4BlcyxfGq4TEWl4xakS49p7p5yW/l8Cd5nZI8CucNlK4Gw0NIOICBDvsfShvBu5d5jZM4DncPyN3N+4e74awYmIJE1x1qzOtng+kVvWJCruXgB+NXu5mb3F3b9csahERBIqFfOWfqU+ij5aof2IiCRa3dT0zeyB+VYBS8o5qJk1A+uBPe7+8nJ+VkQkzqYy8W7pl1PeWQK8hGAilZkM+GWZx/0LYCuwoMyfExGJtbi39Msp7/wQ6HH3nbNeO4CflroTM1sOvAz4QlmRiogkwHRNP6b99MvpvfO2E6x7fRnH/D/AXwO9821gZuuAdQArV64sY9ciItFKxby8U84wDFeY2dMaadPMXg4cdPcNJ9rO3W9197XuvnZwcPDpHFJEpKam6qi880Zgg5l908zebGann8LxrgJeYWY7CMbxudbM/vkU9iMiEkv1VN55J0A40ub1wD+ZWR9wN3AH8IuTPaTl7h8kHKLZzF4IvN/d33BqoYuIxE8qE4xK09ESz4ezyo7K3be5+/929+uAa4GfA68B7qt0cCIiSTOVy9PW3ERLczyTfllP5M7m7ingR+Gr3J/9KWX0+hERSYJgqsR4Jnyo3BO5IiJCOCl6TOv5oKQvIlJRcZ5ABU4x6ZvZoJmpL6WIyCxBeacOkr4FPmJmh4CHge1mNmxmf1O98EREkiVVR+Wd9xL0s7/M3Re5+0LgucBVZvbeqkQnIpIwU9k8HS31kfRvAP7E3R8vLnD3x4A3EDy4JSLS8Oqppd/q7odmL3T3YaC1ciGJiCRXKlM/N3Izp7hORKRhTGULsb6RW87DWRea2bE5lhvQUaF4REQSLSjvxLc3fDlj78z70aUbuSIigXoq75yIkr6INDx3r8+Hs+bwtMbZFxGpB+lcOMJmnfTeORGv0H5ERBKrOIFKnFv6Jdf0zWyM45P7zNZ9Z8UiEhFJqFQ9JX13n3dOWxERmTE/bj2Ud8zssplTJJrZG83s+2Z2k5ktqk54IiLJkYr5/LhQXk3/c4QPYZnZ84GPA18FjgG3Vj40EZFkqauaPtDs7kfC968FbnX37wDfMbNNlQ8t3g6Np1m/4wgHjqXJ5ArkCk6+UMDMaDKjyaC5ycLvocmM9pYmhga6uXhlP+0xHpBJRE7N9Py49ZL0zazF3XPAi4B1p7ifREtl8nzkB1v4zsbd5Aqn1mlpoKedj//Rs/i985ZUODoRidJkJgfUT0v/G8DPwvH0U8C9AGZ2NnC0CrHFjrvzvm9v4o4t+3nTFUO86uIzWLGwk/bWZlqagha+47hDvuAU3Cl48HP5gjOZybNt/xg33bWddV9bz7+8/XIuX7046tMSkQqZDG/kdrfXQdJ394+Z2V3AUuBOdy82c5uA91QjuLi5a+tBfrx5Px+4bg3vfOFZZf/8YmDFoi6uOGsxL7/5Xv7bdx/kJ+97Ac1NerZNpB4Uk35XW3yLH2U9nOXuv3L32919Ysay7e6+sfKhxc8tP/sdqxZ38afPO/Np7aenvYW/vm4Njx2a4M4t+ysUnYhErVje6YpxSz++Q8HFzGPD46zfOcLrn7OS1uan/8/2kvNP57Tedm7/7Z4KRCcicTDd0o9xTV9Jv0R3bT0IwB9cuKwi+2tuMl76rKX8dPvw9AMdIpJsE5kcbS1NtFSgYVgt8Y0sZu5++CBrTu9lWX/lRpx4wbmDZHIFNuwcqdg+RSQ6qUyerhg/jQtK+iXJ5gtsfGKEK88aqOh+LxtaREuT8cvfPWUWShFJoIl0nu4Y38QFJf2SbN13jKlsgUtW9Vd0vz3tLZy3bAGbdo1WdL8iEo1UNhfrcXdASb8kv30iSMoXr1xY8X1fcEYfm/cc5ckesCKSVEFLX0k/8bbuO8bCrlaW9VV+KuALlvVxbCrHE0cmK75vEamtoKav8k7ibds/xprTF2BW+Yeonrk0GLF6+4Hxiu9bRGprIpPTjdykKxSc7QfGWLO0OtMJrB7oAeDxQ0r6IkmXyuTpaldLP9EOjE0xmclz9mk9Vdl/X1cri7rbePzQxMk3FpFYm8jkYv1gFijpn9TukRQAyxd2Ve0YZw5089iwkr5I0k2m87EeggGU9E9q90hwg3X5wupNA3zmQLda+iIJ5+5MZvVwVuLtPhK09M+o4JO4s5050M3BsTTj6VzVjiEi1ZXOFcgXXL13km73SIrB3vaqzoSzeqAbgB1q7YskVmp6WGW19KeZWYeZ/drM7jezLWb20Voe/1TsHp2samkH4MzBIOmrxCOSXBPhsMoahuF4aeBad78QuAi4zswur3EMZdk9kqrqTVyAVYvU0hdJuumWvm7kPskDxQ7preErtuMP5AvO3tFU1Vv6nW3NLO5uY+/RqaoeR0SqZ0LlnbmZWbOZbQIOAj9x9/tqHUOpDo5Nkc171ZM+wLL+TvaOpqp+HBGpjulZs1TeOZ675939ImA58Bwzu2D2Nma2zszWm9n64eHhWoc4rRZ99IuW9Xco6Ysk2GRaLf0TcvdR4G7gujnW3erua9197eDgYO2DC+2ZTvq1a+lrtE2RZJrMxn9SdKh9751BM+sP33cCLwa21TKGchQfzKpmH/2iM/o7mcjkOZZSX32RJJpMF8s78W7p1/ojaSnwFTNrJvjA+ba7/7DGMZSsFn30i4rTMO4ZTdHX1Vr144lIZRUnRY97l82aRufuDwAX1/KYT0fQXbP6rXx4MunvHU1x3rIFNTmmiFRO8UauZs5KsN0jkzW5iQvBjVyAvUd1M1ckicbTedqam2hriXdajXd0ESoUnD2jqZrU8wEGuttpa25ij3rwiCTS2FSW3o54l3ZASX9eB8fSNeujD9DUZCzt75juMSQiyTKezinpJ1kthlSebWlfB/v0VK5IIo1N5ehR0k+uWj6YVbSsr5N9Ku+IJNLYVJbe9vj3vFPSn0ckLf3+Dg6MpckX9ICWSNKMTam8k2i7R1IM9NSmj37R0r5O8gVneCxds2OKSGWovJNwteyjX6RumyLJNTaVZUGHyjuJFfTRr23SX9oXHG/fqG7miiSJu6v3TpIV++jX8iYuBL13APappS+SKJOZPAWHnnYl/USqdR/9or7OVjpbm9mrlr5IooxNBUMw9Kq8k0xR9NwBMAse0FJLXyRZRlMZIGi4xZ2S/hyi6KNftKyvU9MmiiTM6GQWgIUJGCFXSX8OUbX0Iajr71dLXyRRikk/CcOiK+nPIYo++kVL+zvDewqFmh9bRE7N0bC809/VFnEkJ6ekP4ddEXTXLFrW14E7HDimEo9IUoyELf1+1fSTacehSYYW176eD0FLH9DAayIJMjqZpa25KfZTJYKS/lOkc3n2Hk2xanF3JMcv9tXfq4HXRBLjaCpDX1crZhZ1KCelpD/LriMp3GFoIKKW/vQDWmrpiyTFyEQ2EaUdUNJ/ip2HJwAia+n3drTS19k63YNIROLvyGSGhQm4iQtK+k+x43CQbIciSvrBsbvYeVhJXyQpDo2lGextjzqMkijpz7Lz8AS9HS2RPmSxanE3O8K/OEQk/obHlfQTa8fhSYYWd0d6Q2ZocRd7RlJkcuqrLxJ3U9k8Y1M5Jf2k2nl4glURddcsWrm4m4LDHvXgEYm94qRHgz1K+okzlc2z68gkqwd7Io2j+IyASjwi8Tc8HiZ9tfST55ED4xQcnnl6b6RxFHsO7TykpC8Sd8WW/oBa+smzdf8xANYsXRBpHAM9bXS3NbPziHrwiMRd8UHK4nSncaekP8NDe4/R2drMykXR1vTNjFWLu3lcLX2R2Ns9kqKztZlF3eqnnzgbnxjhwhV9NDdF/yj1OUt62L5/LOowROQkdh0JBmhMwhAMoKQ/LZXJ89DeY1y6amHUoQDwzKUL2Ht0itHJTNShiMgJ7B5JRTYq76lQ0g/9escRcgVn7dCiqEMBYE14M3nrPrX2ReLK3dl1ZJIVEZeEy6GkH7p720E6Wpu4YvXiqEMB4IIz+gB4YPdoxJGIyHz2HZ1iLJ3jnNOi7eZdDiV9IJcv8KMH93H12YORzJY1l4GedoYWd7F+50jUoYjIPLbFpMdfOZT0gf+75QAHx9K89rIVUYdynEtXLWLDzhEKBY86FBGZQ7H8em7Ez/aUo+GT/rGpLJ+4Yxtnn9bDNecORh3OcZ53zgBHJjLcrxKPSCxt2DnC6oFuFnQkYyx9gJaoA6iWQsHJu5MvhC93CgUnVwi+pnMFHtp3jJvveoS9oym+/qfPpaU5Xp+BLzx3kOYm444t+7l4ZTx6FYlIIJsv8KvHDvPHlyyPOpSy1G3SP+9/3MFU9uSjVA70tHPLGy7luTG5gTtTf1cb15x7Gret3837XvwM2lvicb9BROCe7cNMZvI875yBqEMpS02TvpmtAL4KLAEcuNXdb6rGsd5z7TkUCk5zs9FsRnNT8GppMprCr6sHe3j28r5YJ9O3Xj3E6z9/gC/c+zjvvubsqMMREYKuml+493EGe9u5Zs1pUYdTllq39HPAf3X3jWbWC2wws5+4+0OVPlC9JMgrzxrgZc9eyqfufJjWZuO1a1fSF+EELyKN7thUln+8+3f8v8cO89FXnE9rzMrCJ1PTpO/u+4B94fsxM9sKnAFUPOnXk0+++tlMpHP8/Y+28fc/2kZfZzCPbpMF4/QYgIFBWY+Cu5fWK6jkvkNldDIqddNKx1ji7sJ9lnjsEvdZzrFLVVf/PqUfuow4Kx/jkYngKfnXXLqcGy5fVWogsRFZTd/MhoCLgfvmWLcOWAewcuXKmsYVR11tLXz5zZfxwO6j3LN9mOHxNEdTWdyDX0L38L9f+IvrOOFHwclVdrOyPnRK32el91f5GEv/dyzj2BU/71L3V4UYK/0PWc6xK7y/Jb0dXHn2QGyGbCmXldpSqOhBzXqAnwEfc/fvnmjbtWvX+vr162sTmIhIHTCzDe6+dq51NS9GmVkr8B3g6ydL+CIiUlk1TfoW/F39RWCru3+6lscWEZHat/SvAm4ArjWzTeHrpTWOQUSkYdW6987PKedOjYiIVFSyOpiKiMjToqQvItJAlPRFRBqIkr6ISAOJ5OGscpjZMLAzosMPAIciOnalJP0ckh4/JP8ckh4/JP8cyo1/lbvPOUFI7JN+lMxs/XxPtSVF0s8h6fFD8s8h6fFD8s+hkvGrvCMi0kCU9EVEGoiS/ondGnUAFZD0c0h6/JD8c0h6/JD8c6hY/Krpi4g0ELX0RUQaiJK+iEgDUdKfh5ldZ2YPm9mjZnZj1PGUwsx2mNmD4eil68Nli8zsJ2b2SPg1VtP9mNmXzOygmW2esWzOmC1wc3hNHjCzS6KLfDrWueL/iJntmWskWTP7YBj/w2b2kmiifpKZrTCzu83sITPbYmZ/ES5P0jWY7xwScR3MrMPMfm1m94fxfzRcfqaZ3RfG+S0zawuXt4ffPxquHyrrgO6u16wX0Az8DlgNtAH3A+dFHVcJce8ABmYt+5/AjeH7G4FPRB3nrPieD1wCbD5ZzMBLgR8TjNR6OXBfTOP/CPD+ObY9L/xdagfODH/HmiOOfylwSfi+F9gexpmkazDfOSTiOoT/lj3h+1aCKWQvB74NvC5cfgvwzvD9u4BbwvevA75VzvHU0p/bc4BH3f0xd88A3wReGXFMp+qVwFfC918BXhVhLE/h7vcAR2Ytni/mVwJf9cCvgH4zW1qbSOc2T/zzeSXwTXdPu/vjwKMEv2uRcfd97r4xfD8GbAXOIFnXYL5zmE+srkP4bzkeftsavhy4FrgtXD77GhSvzW3Ai6yMiZ+V9Od2BrBrxve7OfEvUVw4cKeZbQgnlwdY4u77wvf7gSXRhFaW+WJO0nX5L2H540szSmqxjj8sE1xM0NJM5DWYdQ6QkOtgZs1mtgk4CPyE4K+PUXfPhZvMjHE6/nD9UWBxqcdS0q8vV7v7JcD1wLvN7PkzV3rw92Ci+ugmMWbgs8BZwEXAPuBT0YZzcmbWQzB39V+6+7GZ65JyDeY4h8RcB3fPu/tFwHKCvzrWVOtYSvpz2wOsmPH98nBZrLn7nvDrQeB2gl+eA8U/v8OvB6OLsGTzxZyI6+LuB8L/xAXg8zxZOohl/GbWSpAsv+7u3w0XJ+oazHUOSbsOAO4+CtwNXEFQOivObjgzxun4w/V9wOFSj6GkP7ffAOeEd8/bCG6W/CDimE7IzLrNrLf4Hvh9YDNB3G8KN3sT8P1oIizLfDH/AHhj2IPkcuDojBJEbMyqcf8hwXWAIP7Xhb0vzgTOAX5d6/hmCmvBXwS2uvunZ6xKzDWY7xySch3MbNDM+sP3ncCLCe5L3A28Otxs9jUoXptXA/8R/jVWmqjuWMf9RdBLYTtBbe1DUcdTQryrCXok3A9sKcZMUOu7C3gE+HdgUdSxzor7GwR/emcJ6pZvmy9mgl4OnwmvyYPA2pjG/7UwvgfC/6BLZ2z/oTD+h4HrYxD/1QSlmweATeHrpQm7BvOdQyKuA/Bs4LdhnJuBvwmXryb4MHoU+FegPVzeEX7/aLh+dTnH0zAMIiINROUdEZEGoqQvItJAlPRFRBqIkr6ISANR0hcRaSBK+iIiDURJX0SkgSjpi8Scmf2DmW00s8uijkWST0lfJMbCITVOA94BvDzicKQOKOlLXQpnTXp/+P6XJ9iu38zeVbvI5ozhFjO7aq517j5BMEnIT4GbaxmX1Cclfal77n7lCVb3E8xEFKXLgV/NtcLMFgNdwBiQm2sbkXIo6UvdMLMPmdl2M/s5cO6M5ePh124z+7dwLtLNZvZa4OPAWeEcqp8Mt/teOBHNluJkNGY2ZGZbzezz4fI7wxERMbM3hhN13G9mX5tx3DeEc59uMrPPmVnzHDE/E9ju7vl5TuvDwP8iGETv/Er8O0ljazn5JiLxZ2aXEgyBfRHB7/VGYMOsza4D9rr7y8Kf6SOYYekCDyawKHqrux8Jk/pvzOw74fJzgD9x97eb2beBPzaz3xIk5ivd/ZCZLQr3/UzgtcBV7p41s38E/jPw1VkxXQ/cMc85DQFXAu8jGEnyfGDeUpVIKZT0pV48D7jd3ScBzGyu+Q8eBD5lZp8Afuju986YQm+mPzezPwzfryBI9vuBx919U7h8AzAELAT+1d0PAbh7cb7cFwGXEnxoAHQy9wQ2LwHeMs85/R3wt+7uZrYVtfSlApT0pWG4+3Yzu4RgrPW/M7O7mNXyNrMXAr8HXOHuk2b2U4LxywHSMzbNEyTy+RjwFXf/4LwbmHUB/e6+d451FwF/BFxtZp8JY3jwxGcocnKq6Uu9uAd4lZl1hjOI/cHsDcxsGTDp7v8MfBK4hOAGae+MzfqAkTDhryG4yXoi/wG8JrzhSrG8QzAByavN7LTicjNbNetnryGYHWkunwBe4e5D7j4EXIha+lIBaulLXXD3jWb2LYKZww4STHk527OAT5pZgWCmq3e6+2Ez+4WZbQZ+TFCf/7OwnPIw8/SqmXHcLWb2MeBnZpYnmAHpze7+kJl9GLjTzJrC470b2Dnjx68Hbpu9TzO7Fuhy93+fcZwDZtZjZotmlJBEyqaZs0QiYmYbgee6ezbqWKRxKOmLiDQQ1fRFRBqIkr6ISANR0hcRaSBK+iIiDURJX0SkgSjpi4g0ECV9EZEG8v8BgN+WOubmFcAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(*structure.sld_profile())\n",
    "plt.ylabel('SLD /$10^{-6} \\AA^{-2}$')\n",
    "plt.xlabel('distance / $\\AA$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### `ReflectModel` calculates the generative model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "A `ReflectModel` is made from the `Structure` and is responsible for calculating the reflectivity of the system. `ReflectModel` performs resolution smearing, applies scaling factor and adds a Q-independent constant background. It can use constant `dq/q`, point-by-point, and full resolution kernel smearing. The resolution parameter, `dq`, can be fitted, but this will only be valid if your dataset didn't supply the instrument resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "model = ReflectModel(structure, bkg=3e-6, dq=5.0)\n",
    "model.scale.setp(bounds=(0.6, 1.2), vary=True)\n",
    "model.bkg.setp(bounds=(1e-9, 9e-6), vary=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Let's quickly have a look at the model generated by the structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3hc9ZXw8e+ZGfUuq9hWseRecEUIMCQQDAESjFMgwWQXEtiwpO4mu5uXPMnum21hN9nk3SSQJSawEEIoISwQ1gmhGTBgcMdV2JaLJBdZXbK65rx/zMgIIY1GZdrV+TzPPJr5zZ17z2WMjn5dVBVjjDFmOK5IB2CMMSa6WaIwxhgTkCUKY4wxAVmiMMYYE5AlCmOMMQF5Ih1AKOTk5GhJSUmkwzDGmJiydevWOlXNHVzuyERRUlLCli1bIh2GMcbEFBE5OlS5NT0ZY4wJyBKFMcaYgCxRGGOMCcgShTHGmIAsURhjjAko6kc9iUgK8HOgG9igqg9HOCRjjJlUIlKjEJH7RaRWRHYPKr9KRCpE5KCI3OEv/hTwhKp+Ebg27MEaY8wkF6kaxQPAXcCv+gtExA3cDVwBVAObReQZoBDY5T+sL5RB/XH3SfYeb36vQOS9pwOOy01LYG15MW7XwFJjjHGmiCQKVX1VREoGFZcDB1W1EkBEHgXW4EsahcAOAtSAROQ24DaA4uLiMcX14r5TPLGt2h9j4GPTEj2sWVYwpusYY0wsiaY+igKgasDrauB84KfAXSLyceD3w31YVdcB6wDKysrGtBvTD69fyg+vXxrwGK9XWf7Pz/PmoXpLFMaYSSGaEsWQVPUM8IVgjhWR1cDq2bNnhywel0soyk7iZEtnyK5hjDHRJJqGx9YARQNeF/rLgqaqv1fV2zIyMiY0sMFyUhOoa+sK6TWMMSZaRFOi2AzMEZFSEYkHbgCeGc0JRGS1iKxrbm4e+eBxyE1NoK61O6TXMMaYaBGp4bGPAG8C80SkWkRuVdVe4KvAc8A+4HFV3TOa84arRpGVEk9ThyUKY8zkEKlRT2uHKV8PrB/recPRRwGQluChs8dLd6+XeE80VcqMMWbiOeq3XLhqFOlJcQC0dvaE9DrGGBMNHJUowtVHkZboq4i1dvaG9DrGGBMNHJUowlajSPTVKFqsRmGMmQQclSjCxWoUxpjJxFGJInxNT9ZHYYyZPByVKMLXme2rUbR0WI3CGON8jkoU4ZJmfRTGmEnEEsUYpCV4EIEW66MwxkwCjkoU4eqjcLmE1HiP9VEYYyYFRyWKcPVRgG/SnY16MsZMBo5KFOGUluihpcNqFMYY57NEMUbpiVajMMZMDo5KFOHqowB/jcL6KIwxk4CjEkU4+yiyUuJpOGNLjRtjnM9RiSKc8tMTqG3twusd0/bcxhgTMyxRjNHU9ET6vEq91SqMMQ5niWKM8tITATjV0hnhSIwxJrQsUYxRviUKY8wk4ahEEc5RT/npCQCctERhjHE4RyWKcI56yktLJM4t1DR2hPxaxhgTSY5KFOHkdgmFWckcrW+PdCjGGBNSlijGYcaUZI42nIl0GMYYE1KWKMZhRnYyR+vaUbW5FMYY57JEMQ7FU1Jo7eqlsd2W8jDGOJclinEomZIMwJF6a34yxjiXJYpxmDElBYAjdZYojDHOFfWJQkRmish9IvJEpGMZbMaUZOLcwrun2iIdijHGhExIE4WI3C8itSKye1D5VSJSISIHReSOQOdQ1UpVvTWUcY5VnNvFrNxUKk62RDoUY4wJGU+Iz/8AcBfwq/4CEXEDdwNXANXAZhF5BnADdw76/C2qWhviGMdl3tQ0thxpjHQYxhgTMiFNFKr6qoiUDCouBw6qaiWAiDwKrFHVO4FrxnotEbkNuA2guLh4rKcZtXlT03h6x3FaO3tIS4wL23WNMSZcItFHUQBUDXhd7S8bkohMEZF7gOUi8u3hjlPVdapapqplubm5ExftCOblpwHw7qnWsF3TGGPCKdRNT+OmqvXA7cEcKyKrgdWzZ88ObVADzJvqSxQVJ9s4d0Z22K5rjDHhEokaRQ1QNOB1ob9s3MK5KGC/gswkUhM87LcObWOMQ0UiUWwG5ohIqYjEAzcAz0zEicO5zPiAa7JoejrvVIfvmsYYE06hHh77CPAmME9EqkXkVlXtBb4KPAfsAx5X1T0Tcb1I1CgAlhRmsPdECz193rBe1xhjwiHUo57WDlO+Hlg/0deLRB8FwOLCTLp7D/PuqVYWTQ9vkjLGmFCL+pnZoxGxGkWB73rW/GSMcSJHJYpI9FGAbymP9ESPJQpjjCM5KlFEqkYhIiwpzGRXTVNYr2uMMeHgqEQRSYsLM6g42UpnT1+kQzHGmAnlqEQRqaYn8PVT9PQp+0/aDG1jjLM4KlFEqukJfDUKgHeqrfnJGOMsjkoUkVSQmURuWgLbj1miMMY4iyWKCSIiLC/KZPsxW3LcGOMsjkoUkeyjAFhenMWR+nYaznRH5PrGGBMKjkoUkeyjAFhRnAnAjiqrVRhjnMNRiSLSFhdm4HaJ9VMYYxzFEsUESo73MH9qGtusn8IY4yCOShSR7qMAWF6cyc6qZvq8GrEYjDFmIjkqUUS6jwJgeVEWbV29HKxti1gMxhgzkRyVKKLBcn+Htg2TNcY4hSWKCVaak0Jmcpx1aBtjHMMSxQQ7O/HOhsgaYxzCEkUILC/O4kBtGy2dPZEOxRhjxs1RiSIaRj2Br59CFXZWWfOTMSb2OSpRRMOoJ4ClRZmIYP0UxhhHcFSiiBbpiXHMyUu1kU/GGEewRBEiy4uy2F7VhKpNvDPGxDZLFCGyvDiTpvYejtS3RzoUY4wZF0sUIbK8OAuAbUet+ckYE9ssUYTI7LxUUhM8Np/CGBPzLFGEiNslLCvKtJFPxpiYFxOJQkQ+ISL3ishjIvLRSMcTrOXFmew/2Up7d2+kQzHGmDELeaIQkftFpFZEdg8qv0pEKkTkoIjcEegcqvqUqn4RuB34bCjjnUjLizPp8yq7qiM7AdAYY8YjHDWKB4CrBhaIiBu4G7gaWAisFZGFIrJYRJ4d9Mgb8NHv+j8XE5YV+Tq0t9sMbWNMDPOE+gKq+qqIlAwqLgcOqmolgIg8CqxR1TuBawafQ0QE+DfgD6q6bajriMhtwG0AxcXFExb/eGSnxFMyJdlGPhljYlqk+igKgKoBr6v9ZcP5GnA5cJ2I3D7UAaq6TlXLVLUsNzd34iIdp+XFNvHOGBPbgkoUIvI1EckKdTDDUdWfquq5qnq7qt4z3HHRsijgQCuKMznd2kVNU0ekQzHGmDEJtkaRD2wWkcf9ndAyzuvWAEUDXhf6y8YlWhYFHKh/4p0NkzXGxKqgEoWqfheYA9wHfB44ICLfF5FZY7zuZmCOiJSKSDxwA/DMGM91VjTWKOZNTSMxzmWJwhgTs4Luo1BfI/tJ/6MXyAKeEJEfBPqciDwCvAnME5FqEblVVXuBrwLPAfuAx1V1zxjvYWCMUVejiHO7WFKQyTZbSdYYE6OCGvUkIn8F3ATUAb8E/k5Ve0TEBRwAvjXcZ1V17TDl64H1o444cJyrgdWzZ8+eyNOO2/LiTP779SN09faR4HFHOhxjjBmVYGsU2cCnVPVKVf2tqvYAqKqXIYazRko01ijA10/R3edlz/GWSIdijDGjFmyimKmqRwcWiMhDAKq6b8KjGqNo7KMAX40CrEPbGBObgk0Uiwa+8M+sPnfiwxmfaK1R5KcnUpCZZDveGWNiUsBEISLfFpFWYImItPgfrUAt8HRYInSIZcW2kqwxJjYFTBSqeqeqpgE/VNV0/yNNVaeo6rfDFGPQorXpCWB5USY1TR3UtnRGOhRjjBmVkWoU8/1PfysiKwY/whDfqERr0xPAihm+iXdbbN0nY0yMGWl47DfxLbT3oyHeU+CyCY/IoRYXZJAS7+aNQ3V8bPG0SIdjjDFBC5goVPU2/8+PhCec8YnWeRTgm3h3wcwpvH6wPtKhGGPMqAS7KOA7/o7tsS7ZERbR3PQEcNHsHA7XnaG6sT3SoRhjTNCCHR67GugDHheRzSLytyISHZs+xJCL5+QA8PrBughHYowxwQt2UcCjqvoDVT0XuBFYAhwOaWQONCcvlWkZiTy/91SkQzHGmKAFvSigiMwQkW8BjwLzCbC+kxmaiPDxxdN45d3TNLf3RDocY4wJSrB9FG8B/wO4getVtVxVhxoJFVHRPI+i3+ql0+npU9bvPhHpUIwxJijB1ihuUtUV/gl4lSGNaByivTMbYElhBgumpXPfxsP0eT+4PWpHdx+nW7vwDvGeMcZEQsDhsSLyZ6r6a+DjIvLxwe+r6o9DFplDiQhfu2w2X354G/e8coivfGQ23b1eXq6o5clt1by0v5aePiUjKY615cX81ao5JMXb0uTGmMgZacJdiv9n2hDv2Z+8Y3T1OVP5+OJp/PC5Cp595wTVje20dvaSk5rAn19QQnF2EpuPNPKLVw+xoaKWX//F+eSkJkQ6bGPMJCW+jetGOEjkIlV9faSyaFFWVqZbtmyJdBgBdfd6ufe1St48VE9RdhJXLprKxbNz8Ljfaw3cUFHL7b/eyqzcVH73pZUkxlnNwhgTOiKyVVXLPlAeZKLYpqorRiqLFrGQKIL14r5T3PrgFm6+cAb/uOacSIdjjHGw4RLFSH0UFwIrgVwR+eaAt9LxjYCKKtG8hMdYrVqQzxcuKuG/Xz/Cp88tZElhZqRDMsZMMiONeooHUvEllLQBjxbgutCGNnqxMOppLL55xVympMTzL8/uI5gaoDHGTKSRFgV8BXhFRB4YvBWqCZ+0xDi+ccVcvvvUbjYerONDc3IjHZIxZhIJdh7FL0XkbJuHiGSJyHMhiskM4fqyQnJSE7hvo62cYowJr2ATRY6qnt3HU1UbgbzQhGSGkuBxc/OFM9hQcZoDp1ojHY4xZhIJNlF4B64WKyIzsHkUYXfj+cXEuYXHNldFOhRjzCQSbKL4DrBRRB4SkV8DrwJRt2e2001JTeAj8/J4eudxevu8kQ7HGDNJBLvM+B+BFcBj+FaPPVdVrY8iAj61ooDTrV1stD0tjDFhEuzqsQJcBaxQ1WeBZBEpD2lk7117gYjcIyJPiMiXwnHNaPaR+XmkJ3p4ZufxSIdijJkkgm16+jlwIbDW/7oVuHukD4nI/SJSKyK7B5VfJSIVInJQRO4IdA5V3aeqtwOfAS4KMl7HSvC4WbUgn5f311rzkzEmLIJNFOer6leATjg76ik+iM89gK8mcpaIuPElmauBhcBaEVkoIotF5NlBjzz/Z64F/hdYH2S8jnb5gnwa23vYdqxp5IONMWacgk0UPf5f8AogIrnAiH/OquqrQMOg4nLgoKpWqmo3vj6PNaq6S1WvGfSo9Z/nGVW9GvjccNcSkdtEZIuIbDl9+nSQtxWbPjw3hzi38MI+21LVGBN6wSaKn+Lb4S5PRP4V2Ah8f4zXLAAGju+s9pcNSUQuFZGfisgvCFCjUNV1qlqmqmW5uc6euZyWGMcFM6fw0v7aSIdijJkERtqPAgBVfVhEtgKrAAE+oar7QhrZe9feAGwI5lgnLgo4nItn53DnH/ZT29JJXnpipMMxxjhYwBqFiGT3P4Ba4BHgN8Apf9lY1ABFA14X+svGzamLAg5l5awcAN6srI9wJMYYpxupRrEVX7+E+F/3z8YW//OZY7jmZmCOiJTiSxA3ADeO4TwfMJlqFAunp5Oe6OH1g3WsWTZsy50xxozbSH0Uf66qM4EFqlqqqjP9j1J/eUAi8gjwJjBPRKpF5FZV7QW+CjwH7AMeV9U9470RmFw1CrdLuHDWFF4/aDUKY0xojVSj+AlwLvAGvpnZo6Kqa4cpX08IhrpOphoFwAUzp/DcnlOcaO5gWkZSpMMxxjjUSImiR0TWAYUi8tPBb6rq10MT1tio6u+B35eVlX0x0rGEw4riLAC2H2ti2mJLFMaY0Bip6eka4CWgA19/xeBHVBGR1SKyrrm5OdKhhMWCaenEe1xsO9oY6VCMMQ420g53dcCjIrJPVXeGKaYxm2w1iniPi8UFGWyvshnaxpjQCXbCXYeIvNi/ZpOILBGR74YwLhOkFcWZ7KppprvX1n0yxoRGsIniXnz7T/QAqOo7+Ia1RpXJ1vQEsLw4i+5eL3tPtEQ6FGOMQwWbKJJV9e1BZb0THcx4Tabhsf36O7Stn8IYEypBLeEB1InILN5bFPA64ETIojJBm5qRyNT0RN6pHl0/RXVjOw++cYTa1i7KS7P5TFkRce5g/24wxkwmwSaKrwDrgPkiUgMcJsBKrpEy2eZR9FtSmME71cE3t22oqOX2X2/F64Wc1Hie3nGcxzdX8eAt5WQmB7N6vDFmMgl2K9RKVb0cyAXmA5cAF4cysLGYjE1P4EsUlXVnaO7oGfHY/SdbuO2hrczKTWXD313KG99exV03LmffyVa+8MBmemwzJGPMICMtCpguIt8WkbtE5AqgHbgZOIhvxzkTBZYUZgKwpyZwraLPq3zzsZ2kJ3p48JZypmf6Julds2Q6P7p+KduPNXHXSwdDHq8xJraMVKN4CJgH7AK+CLwMXA98UlXXhDg2E6Qlhb4a1M4Rmp/+Z3sNe0+08L1rF5GTmvC+91Yvnc6aZdP5r1cOUdPUEbJYjTGxZ6REMVNVP6+qv8C3X/ZC4EpV3RH60EywMpPjKc5ODtih7fUqP3vpAIsLMvj44mlDHvOtq+YD8JMX3g1JnMaY2DRSojjb6K2qfUC1qnaGNqSxm4zzKPqN1KG97VgjR+vbueXiEkRkyGMKMpNYe14RT20/Tm1r1H7NxpgwGylRLBWRFv+jFVjS/1xEom6G12TtzAZYWphJTVMH9W1dQ77/9I7jJMa5+OjCqQHP8/mLSunu8/LwpmOhCNMYE4MCJgpVdatquv+RpqqeAc/TwxWkGVl/P8VQtYqePi//u+sEly/IJyUh8Ijo0pwUPjQnhye3V6OqAY81xkwONsPKIRYVZCAydKLYeLCOhjPdQe+Et2ZZAVUNHWw7ZosNGmMsUThGaoKH2bmpQ3ZoP7PjOBlJcVwyNzeoc125KJ8Ej4und0zIVubGmBhnicJBlhRmsrO6+X1NRh3dfTy35yQfWzyVeE9wX3daYhyrFuSxftdJvF5rfjJmsnNUopjMo57A109R19bFyZb3Riy9sO8U7d19XLs0uGanfqvm51PX1mWr0hpjnJUoJvOoJxgw8a7qvUT59I7j5KcnUF6aPapzfdjfTLWhonbiAjTGxCRHJYrJbsG0dDwuOdtP0dzewyvv1rJ6yXTcrqHnTgwnNy2BxQUZbKg4HYpQjTExxBKFgyTGuZk3Ne3syKf1u0/Q06dBj3Ya7NJ5uWw71khL58iLDRpjnMsShcOcOyOLbccaae/u5bdbqpiZm8I5BWOb8nLhzCl41TZFMmays0ThMB9bPI327j7uXL+fbceauLG8eNglO0ayrDgTt0vYfKRhgqM0xsQSSxQOU16SzZy8VB7adJRpGYncUF485nMlx3s4Z3o6m49YjcKYySzYHe5MjHC5hHtvKuORt4+xZlkBqSMs2TGSspJsHtp0lK7ePhI87lF99q3Keta9WsnJlk7OK8nm66vmkJ1iO+gZE2tiokYhIikiskVErol0LLGgJCeFb39sAQunj385rvNKsunu9bJ7hE2RBntiazVr793E7uPNZKfE8+tNR7n6J69yuO7MuGMyxoRXSBOFiNwvIrUisntQ+VUiUiEiB0XkjiBO9X+Ax0MTpQlkWZFv97zR7Mn9TnUTd/zuHVbOyuGlv7mUh249n6e/ehG9fcpN979Fq42iMiamhLpG8QBw1cACEXEDdwNX49sIaa2ILBSRxSLy7KBHnn8L1r2AzfyKgPz0BHJSE9gVZI2iz6t864l3yEtL4O4bV5xdrXbR9AzW3XQuNY0dfH/9vlCGbIyZYCHto1DVV0WkZFBxOXBQVSsBRORRYI2q3gl8oGlJRC4FUvAllQ4RWa+q3iGOuw24DaC4eOwduOb9RITFBelBNz09s7OG/Sdb+dna5WQkx73vvXNnZHPLRaXc9/phbl5ZwvyptlK9MbEgEn0UBUDVgNfV/rIhqep3VPWvgd8A9w6VJPzHrVPVMlUty80NbpVUE5zFBRkcrG2jvbs34HGqyi9eqWT+1LRht1v92mVzSEvw8OM/2XarxsSKmOjMBlDVB1T12UDHTPZFAUPlnIIMvAr7RlggcOvRRvafbOXmlSW4hlkyJCM5jpsuLOH5fac4Wm8d28bEgkgkihqgaMDrQn/ZuE32RQFDZUmhr0N71wgd2o9vqSI1wcO1S6cHPO6mC2fgcQkPvXl0wmI0xoROJBLFZmCOiJSKSDxwA/DMRJzYahSh4evQjmfP8eFrFL19Xp7fe4rLF+SNuN1qXnoil83P4+mdx+ntG7Il0RgTRUI9PPYR4E1gnohUi8itqtoLfBV4DtgHPK6qeybielajCA0RYcG09ICJ4u0jDTS293DVOVODOucnlhVwurWLNw7VT1SYxpgQCfWop7XDlK8H1k/09URkNbB69uzZE33qSW/R9Azu21hJd693yJ3yntt9ksQ419l9LEbykfl5pCV6eGpHTdCfMcZERsx0ZgfDahShs3B6Oj19yoHa1iHff7niNBfPziU5Pri/PRLj3Hx04VRe2HvKmp+MiXKOShTWRxE6i/zLgewdovnpeFMHxxraWTlryqjOedn8PFo6e9le1TQhMRpjQsNRicJqFKFTMiWFpDj3kP0Ubx329TOcP3N0261+aG4Obpfw8n6bdG9MNHNUojCh43YJC6alDVmj2HSogYykOBaMcqZ1emIcZTOyeNm2WzUmqjkqUVjTU2gtmp7B3hMteL36vvK3DtdzXkn2sJPsArl0Xh77TrRwurVrTDHtOd7Mj/5UwQ/+uH/UK9waY4LjqERhTU+htXB6Om1dvVQ1tp8tO9ncyZH6di4YZbNTv/7PjWUXvZ+9eIBrfraRn284xLpXK7nmZxv5tz/sR1VH/rAxJmiOShQmtPo7tAf2U/T3T1wwc3Qd2f3OKcggOd7NW5Wjm0/x4BtH+NHz77Jm6XS2ffcKtv3DFawtL+KeVw7xw+cqxhSLMWZojkoU1vQUWnPz03C75H39FJsqG0hL9LBg2thWgo1zuzh3RhZvHQ6+RlF5uo1/Xb+PVfPz+NFnlpGRHEd6Yhzf/+Ri1pYX8fMNh3jlXev3MGaiOCpRWNNTaCXGuZmdm8qe4+8l4rcq6ykvycY9hv6JfuUl2VScaqWpvTuo4//52b0keFzc+enF77uuiPB/Vy9iVm4Kf//Ubjp7+sYckzHmPY5KFCb0Fk5PZ69/FdmTzZ1U1p0Zc7NTv/NKs1GFbccaRzx2V3UzL1ec5vZLZpGXlviB9xPj3Hzv2kUca2jnN28dG1dcxhgfSxRmVBZNT+dUSxd1bV28WVkHwMrZ40sUiwsycAnsrBq5yfDe1ypJS/Tw5xfOGPaYD83Jpbw0m3tf8y05YowZH0sUZlQWDujQfuNgPZnJo58/MVhKgoe5+WnsrA48Q7u5vYc/7jnJp5YXkJ4YF/DYr3xkNieaO3lq+4SsYG/MpOaoRGGd2aG30N9pvbOqiVfePc2FM6eMaf7EYEsLM9lZ1RRwaOsz7xynu9fLdecWDXtMvw/PyWFefhoPv2V7XhgzXo5KFNaZHXqZyfEsKczgx8+/S21rF2uWBd6kKFhLizJpbO+hqqFj2GOe3FbN/KlpnFMwcg1GRPjMeUXsrG5m/8nAO/MZYwJzVKIw4XHrxaUAFGcnc9n8/Ak559IiX3LfMUzzU21LJ9uPNXHNkmmIBFeD+eTyAuLcwuObqyckRmMmK0sUZtSuXTqdJ7+8kie/vHLIvSnGYm5+GolxLnYOs5Lsi/6FAy9fGHxiyk6JZ9X8fJ595/gHlh0xxgTPEoUZNRFhRXEWOakJE3bOOLeLc6ZnDJ8o9p2iIDOJeflpozrv1YunUtvaxfaqkYfeGmOGZonCRI2lRZnsqmmmZ9BGRh3dfWw8WMcVC/ODbnbqd9n8POLdLv6w6+SYYlJVdtc08/SOGrYcabCaiZmUQroVarjZVqixbVlRJvdtPEzFyVbOKXhvQMIbh+ro7PGyakHeqM+ZlhjHRbOn8Mc9J/nOxxeMKtGcaunkG4/teN++3vOnpvHjzyw7O0zYmMnAUTUKG/UU25YXZwIfnKH9ckUtyfFuykvHtkLt1edMo7qxY8hNl4ZT29LJp//rDXZWNfEP1yzkT9/4MD+6fikNZ7q5/p432Hp09KvdGhOrHJUoTGwryEwiNy2B7cfe66dQVTZUnGblrBwSPO4xnXfVgjxE4IV9p4I6vrfPy5cf3kZ9Wze/+eIF3HJxKXPz0/j0uYX8/msXk5+eyBd/tZWTzZ1jiseYWGOJwkQNXyd5JtsH1Cgq685Q3djBJfNyx3zeKakJLC/K5MV9wW25+vBbx9hytJF//eQ5LC3KfN97+emJ3HtzGZ09fXz90e3WZ2EmBUsUJqosL87iSH079W2+He82+LdJvXTu2BMFwKoF+eyqaeZUS+BaQH1bF//xpwounp3DJ5cXDHnMrNxUvrd6EW8fbuB322yOhnE+SxQmqqwozgJg8xFfreKVd08zMzeFouzkcZ23vyP8pf2BaxUPvnGEtq5evnftwoAd39edW8iK4kz+/Y/7aensGVdsxkQ7SxQmqiwryiQ1wcMr756mqb2bNw/Vcdm80Y92GmxefhoFmUm8GKCforOnj1+/dYxV8/OZnRd4vobLJXzv2kXUtXXz4OtHxh2fMdHMEoWJKvEeFx+ak8Pze0/xxNZqevqUaydgPSkR4fIFeWw8WDfshkZPba+h4Uz32SVKRrKkMJPLF+Txy42HaevqHXeMxkSrqE8UInKpiLwmIveIyKWRjseE3ufOn0FdWxf/8r/7WFaUyeKCiRnufNmCfDp7vLxxqO4D76kq9208zKLp6VwwM/hhuF9fNYfmjh4efOPIhMRoTDQKaaIQkftFpFZEdg8qv0pEKkTkoIjcMcJpFGgDEgHrOZwELpo9hc+WFVGak8K/fOKcUc/GHs4FM7NJjncPOfrp1QN1HKht49aLS0d1vSWFmVw6L5f7N/zQy3YAAA3hSURBVB62rVeNY4V6ZvYDwF3Ar/oLRMQN3A1cge8X/2YReQZwA3cO+vwtwGuq+oqI5AM/Bj4X4phNhIkI/37dkgk/b4LHzYfm5PDS/lpU9X0J4b6Nh8lLS+CaJaNv5vrSJbP47LpN/HZLFX9+YcmoP9/d6+Vw3RkASnNSJmyhRWMmSkgThaq+KiIlg4rLgYOqWgkgIo8Ca1T1TuCaAKdrBIZdhU5EbgNuAyguLh5H1MbJVi3I57k9p9h7ooVF031NWgdOtfLqu6f524/OHdMv6fLSbJYXZ7LutUrWlhfjcQd3jo7uPn760gF+vekorZ2+Po6UeDfXlxXxjcvnkpEceBc/Y8IlEn+6FABVA15X+8uGJCKfEpFfAA/hq50MSVXXqWqZqpbl5o5vzL1xro/M88/S3vte89O9r1WS4HFx4/nD78MdiIhw+yWzqGroYP3u4BYfrGpo5+M/e43/2nCIS+bm8pMblvGTG5Zx5TlTeWjTUa78z1fZetRWvDXRIerruKr6pKr+pap+VlU3BDrWtkI1I8lNS6C8JJsntlXR51WqGtp5clsNa8uLyU6JH/N5r1iQz6zcFO7ZcCjgdq4Ax+rbuf6eN6lv6+bhvzifu25cwZplBaxZVsCPP7OMp758EYlxLj73y01sqAhuNrkxoRSJRFEDDNz0uNBfNm62KKAJxudXllDV0MFv3jrKPz27F5fLVyMYD5dL+MtLZrH3RAuvHvjgqKp+rZ09/MWvNtPR08ejt13ARbNzPnDM4sIMnvjSSmblpvKXD221BQhNxEUiUWwG5ohIqYjEAzcAz0zEia1GYYJx5aKpXDhzCn//9B6e33uKb105j6kZieM+7yeWFTA1PZF7Nhwa8n2vV/nGYzs4dPoMP//cChZMG36p8pzUBB68pZzpmUnc8sAWDta2jTs+Y8Yq1MNjHwHeBOaJSLWI3KqqvcBXgeeAfcDjqrpnIq5nNQoTDJdL+OXNZXzrqnncfeOKoCfYjSTe4+LWi0t5s7J+yJ36fvR8BS/sq+Ufrlk4ZE1isJzUBH51SzlxbuELD7x9dv0rY8JNRmpPjSUDNi764oEDByIdjpmE2rp6WXnniywpzOShW8vPDsH9/c7jfO2R7awtL+L7n1w8qrka2481csO6TZxTkMHDf3E+iXFjW27dmJGIyFZVLRtcHvWd2aNhNQoTaakJHr511Xw2Hqzj5/4mqOf3nuJvfruTshlZ/OO1o59AuLw4i//32WVsPdrI3z3xji1tbsLOtkI1ZoJ97vxiNlXW88PnKvjNW8eoaepgcUEG995UNubJdB9bPI1vXTWPH/yxgpIpyfzNR+eN+Jmj9Wf43bYath9rpKqhna5eLzmpCczJS+XiOTlcvjCf9ESbq2FG5qimp35lZWW6ZcuWSIdhJrHePi8PbTrK24cbWFyYwRdWlpIUP74mI1Xljt/t4rEtVdxx9Xz+8sMzP1A76fMqL++v5aFNR3nl3dO4BOZPTWdWXioJHhe1rV3sqWmm/kw3SXFu1iybzp9dMON9e5SbyWu4pidLFMbEkJ4+L3/92A7+950TXL4gjy9dOovSnFSON3Xwpz0n+d22GmqaOshPT2BteTE3nFf8gRFdXq+yo7qJxzdX8fSO43T09HF+aTa3XFzK5QvycbsmZm0tE3smRaKwzmwzGagqv3ztMP/5wruc6X5vIUKXwEWzc7ixvJjLF+YTF8RSIs0dPfx2SxX//foRapo6KMpO4vMrS/lMWSFp1iw16UyKRNHPahRmMmjp7GHjgTpONneSl57AeSXZ5KePbT5Ib5+XF/ad4v6NR3j7SAMJHhfnz5zC+aXZzM5LZXpGEmmJHlITPSR4XMR7XMS7XRO2sq+JDpYojDFB2VXdzJPbq3ntQN2IE/3iPS4S3L7E0Z9AEjxuUhM9ZCXHkZkcT3ZKPJnJcWQlx3+gLDMp3lbLjSLDJQob9WSMeZ/FhRksLvR1brd29nDo9BlOt3bR2tlDW1cv3b1euvwP3/O+s2X9r1s7e6lu7GB3TQuN7d109XqHvV5qgudsIslOiacoO4mSKSkUZyczw/9zNAMBVJUz3X00d/j2Mo9zCamJHpLjHfXrLqysRmGMCbmO7j4a27tpbO+mqb2HhjPdNLV309jec7assb2burYujtW309L5/q1l89ISfMljSjJZyXG4RPCq0tLRS1OH7/NNA87V3ffBxJSW6GFqeiIlOSnMzU9lbn4ac/PTmJmbQoIn+ETk9SpNHT00nOmi4UwPvV4vCR43CR4XOakJ5KTGB73UfLSZFDUKY0x0Sop3kxSfxPTMpKCOb2rv5mh9O0cb2jlad4ajDe0cq2/ntQOnaenoxauKCGQk+ZqvMpLimDElmeXFmWT6m7gykuIQgV6vL6GcaunkRHMHlafP8PL+Wnr9ExfdLmFqeiIFWUkUZCaRmuAh3uPCq0prZy+tnT00+pNb4xlfsgs059ElvlWKp2UkMT0zkekZSUzLTKIgM5FpGUnkpScQ73bhcbvwepW2rl7auno509VLa1cvbZ0DXnf20tnbR5zLhcctxHtcZCTFkZ0cT1aKrwaWnRJPVnJ8SEerWaIwxkSdzOR4MpPjWVqUGZLz9+8qWHGqlYOnWqlq7KCmsYO3DzfQ3u1rXhMR0hI9pCV6yEyOZ05eKlkp8UwZ9Ava4xa6e7109nipa+uitqWTky2dnGjupOJkKy/vP03HOLbJjfe46PMqfQGykwhkJsWRnRLPf1y/lOXFWWO+3lAsURhjJp14j4t5U9OYNzUt5NdSVZraezje3MHxpk5Ot3bR0+el16u4BFISPKQmeM7+TPePLktN8JAS78Hlryl4vUp3n5fmDl/tpuFMN/X+Wk79mW4aznRR39ZNetLED2t2VKKwzmxjTLQREbJSfE1F/dvvjoXLJSS63CTGucc8DHrM1w7r1ULMFgU0xpiJ56hEYYwxZuJZojDGGBOQJQpjjDEBWaIwxhgTkKMShYisFpF1zc3NkQ7FGGMcw1GJwkY9GWPMxHNUojDGGDPxHLkooIicBo4GcWgOUBficMLJSffjpHsBZ92Pk+4FnHU/472XGaqaO7jQkYkiWCKyZaiVEmOVk+7HSfcCzrofJ90LOOt+QnUv1vRkjDEmIEsUxhhjAprsiWJdpAOYYE66HyfdCzjrfpx0L+Cs+wnJvUzqPgpjjDEjm+w1CmOMMSOwRGGMMSYgxyYKEblKRCpE5KCI3DHE+wki8pj//bdEpGTAe9/2l1eIyJXhjHsoY70XESkRkQ4R2eF/3BPu2IcSxP18WES2iUiviFw36L2bReSA/3Fz+KIe2jjvpW/Ad/NM+KIeXhD3800R2Ssi74jIiyIyY8B7sfbdBLqXWPxubheRXf6YN4rIwgHvje93mqo67gG4gUPATCAe2AksHHTMl4F7/M9vAB7zP1/oPz4BKPWfxx2j91IC7I709zGG+ykBlgC/Aq4bUJ4NVPp/ZvmfZ8Xivfjfa4v09zGG+/kIkOx//qUB/9Zi8bsZ8l5i+LtJH/D8WuCP/ufj/p3m1BpFOXBQVStVtRt4FFgz6Jg1wIP+508Aq0RE/OWPqmqXqh4GDvrPFynjuZdoNOL9qOoRVX0H8A767JXA86raoKqNwPPAVeEIehjjuZdoFMz9vKyq7f6Xm4BC//NY/G6Gu5doFMz9tAx4mQL0j1Qa9+80pyaKAqBqwOtqf9mQx6hqL9AMTAnys+E0nnsBKBWR7SLyioh8KNTBBmE8/31j8bsJJFFEtojIJhH5xMSGNiajvZ9bgT+M8bOhNp57gRj9bkTkKyJyCPgB8PXRfDYQz6hCNbHmBFCsqvUici7wlIgsGvSXh4mcGapaIyIzgZdEZJeqHop0UMEQkT8DyoBLIh3LeA1zLzH53ajq3cDdInIj8F1gQvqKnFqjqAGKBrwu9JcNeYyIeIAMoD7Iz4bTmO/FX9WsB1DVrfjaJueGPOLAxvPfNxa/m2Gpao3/ZyWwAVg+kcGNQVD3IyKXA98BrlXVrtF8NozGcy8x+90M8CjQXxMa/3cT6U6aEHX8ePB1ppXyXsfPokHHfIX3dwA/7n++iPd3/FQS2c7s8dxLbn/s+DrBaoDsaP9uBhz7AB/szD6Mr7M0y/88YvczznvJAhL8z3OAAwzqnIzG+8H3C/MQMGdQecx9NwHuJVa/mzkDnq8Gtvifj/t3WsRuPAz/YT8GvOv/h/Adf9k/4fvLASAR+C2+jp23gZkDPvsd/+cqgKtj9V6ATwN7gB3ANmB1pO8lyPs5D1876hl8tbw9Az57i/8+DwJfiNV7AVYCu/z/A+8Cbo30vQR5Py8Ap/z/pnYAz8TwdzPkvcTwd/OTAf+/v8yARDLe32m2hIcxxpiAnNpHYYwxZoJYojDGGBOQJQpjjDEBWaIwxhgTkCUKY4wxAVmiMCYMRKRQRJ72r6xaKSJ3iUhCpOMyJhiWKIwJMf8CjU8CT6nqHGAOkIRvPR5jop7NozAmxERkFfB/VfXDA8rSgaNAkaq2RSw4Y4JgNQpjQm8RsHVggfoWZjwCzI5EQMaMhiUKY4wxAVmiMCb09gLnDizwNz1Nxbf2jjFRzRKFMaH3IpAsIjcBiIgb+BFwl6p2RDQyY4JgicKYEFPfiJFPAteJyAF8q8h6VfVfIxuZMcGxUU/GhJmIrAQeAT6pqtsiHY8xI7FEYYwxJiBrejLGGBOQJQpjjDEBWaIwxhgTkCUKY4wxAVmiMMYYE5AlCmOMMQH9f2Y/fO2uDz68AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "q = np.linspace(0.005, 0.3, 1001)\n",
    "plt.plot(q, model(q))\n",
    "plt.xlabel('Q')\n",
    "plt.ylabel('Reflectivity')\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### `Objective` combines the model and data, calculating statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "An `Objective` is made from a model and dataset. Here we use a `Transform` to fit as logY vs X."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "objective = Objective(model, data, transform=Transform('logY'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The `Objective` can calculate statistics for the fitting system. Note how the log-probability is the sum of the log-prior and log-likelihood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "34376.69092305155 -5.013178251637257 -16157.428065740693 -16162.441243992329\n"
     ]
    }
   ],
   "source": [
    "print(objective.chisqr(), objective.logp(), objective.logl(), objective.logpost())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### `CurveFitter` does the fitting/sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The final setup step is to create a `CurveFitter` from the `Objective`. These objects do the fitting/sampling. Let's do an initial fit with differential evolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "51it [00:06,  7.53it/s]\n"
     ]
    }
   ],
   "source": [
    "fitter = CurveFitter(objective)\n",
    "fitter.fit('differential_evolution');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "An `Objective` has a plot method, which is a quick visualisation. You need matplotlib installed to create a graph. You can see that the fit looks good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXgUVdbA4d/p7GwCAUVZBpB9TSCKjYIRFFCRRUBcQUGDuOsoisznOCOCIOPOCFHUwQUVUdkEZIugBDBI2FFAQEFRCLIIhGz3+6M6SafTSToxnUrS532eetJVXV11utA6fZe6V4wxKKWUCjwOuwNQSillD00ASikVoDQBKKVUgNIEoJRSAUoTgFJKBahguwMojjp16pjGjRvbHYZSSlUoGzZsOGKMqeu5vUIlgMaNG5OUlGR3GEopVaGIyH5v27UKSCmlApQmAKWUClCaAJRSKkBVqDYApdRfl56ezoEDB0hNTbU7FFXKwsPDadCgASEhIT7trwlAqQBz4MABqlevTuPGjRERu8NRpcQYQ0pKCgcOHKBJkyY+fUargJQKMKmpqURGRurNv5IRESIjI4tVstMEoFQA0pt/5VTcf1dbE4CI9BGR70Vkt4g84c9zxcdD797WX6VU8YhYi6pcbEsAIhIETAWuBtoAN4lIG3+c69ZbYdQo+PJL6+/AgZCY6I8zKaVUxWFnCeBiYLcx5kdjTBrwIdC/tE8SHw/vv5932+efw6WXamlAqfLu6aefpn79+kRFRdGuXTvmzZuXs33KlCn59g8KCsrZd8iQIZw+fRqAxYsX07JlS5o1a8Zzzz2Xs//evXvp0qULzZo1Y+jQoaSlpQGwatUqOnXqRHBwMJ988kmec/Tp04eaNWvSt2/fPNuXL19Op06diIqK4rLLLmP37t0AvPPOO9StW5eoqCiioqJ48803cz7z+OOP065dO9q1a8dHH31UCleseOxMAPWBn93WD7i2lao5c7xvN8YqDWhJQCnf2fH/y8MPP0xycjKzZ89mxIgRZGVlFbhvREQEycnJbN26ldDQUKZNm0ZmZib33nsvixYtYvv27cyaNYvt27cD1g344YcfZvfu3dSqVYsZM2YA0KhRI9555x1uvvnmfOd47LHHePfdd/NtHz16NO+//z7JycncfPPNjB8/Pue9oUOHkpycTHJyMnfeeScACxcu5LvvviM5OZl169YxZcoUTpw48ZeuVXGV+0ZgEYkTkSQRSTp8+HCxP1+3LjzJsyTTkfVcxDd0ZSWxTCcOMFx9denHrFRl4n7T79mz9JLAzJkz6dChAx07duS2224rcv/WrVsTHBzMkSNHfDp+t27d2L17N+vXr6dZs2Y0bdqU0NBQbrzxRubOnYsxhhUrVjB48GAAhg8fzueffw5Y44516NABhyP/LbJnz55Ur14933YRybmBHz9+nAsuuKDQ+LZv30737t0JDg6matWqdOjQgcWLFxe4/+LFi2nVqhWdOnXigQceyFcCKQk7E8BBoKHbegPXtjyMMfHGmBhjTEzduvkGsyvSrl1wmLrspQm/cy6nqEoX1hHHGzRjN8ePl/wLKBUIEhJyX6el5V0vqW3btjF+/HhWrFjBpk2bePnll4v8zLp163A4HPhyH8jIyGDRokW0b9+egwcP0rBh7q2mQYMGHDx4kJSUFGrWrElwcHCe7SX15ptvcs0119CgQQPeffddnngit1/LnDlz6NChA4MHD+bnn62Kj44dO7J48WJOnz7NkSNHWLlyZc57nlJTU7nrrruYP38+GzZs4NChQyWO052dCeBboLmINBGRUOBGYF5pnyQ8HN4gjoF8Tl8W0oulXMMXANTnIBERpX1GpSqX2Njc16GheddLasWKFQwZMoQ6deoAULt27QL3ffHFF4mKiuLRRx/lo48+KrSr45kzZ4iKiiImJoZGjRoxcuTIvx6sj1588UW++OILDhw4wB133MEjjzwCwHXXXce+ffvYvHkzV111FcOHDwegV69eXHPNNXTt2pWbbroJp9NJUFCQ12Pv3LmTJk2a0Lx5c0SEW2+9tVRiti0BGGMygPuAJcAO4GNjzLbSPk8bL/2KjnMOAOdwnLNntR1AqcI4nbmvly/Pu14WstsAVq9eTbdu3QrdN7sNIDk5mVdffZXQ0FDq16+f55f1gQMHqF+/PpGRkRw7doyMjIw820vi8OHDbNq0iS5dugBWnf+aNWsAiIyMJCwsDIA777yTDRs25Hxu3LhxJCcns3TpUowxtGjRokTnLylb2wCMMV8YY1oYYy40xjzrj3MMG5Z/m3sCyMqCmTP9cWalKp/Suvn36NGD2bNnk5KSAsDRo0dL58BeXHTRRezatYu9e/eSlpbGhx9+SL9+/RARrrjiipxePv/73//o379kHRFr1arF8ePH+eGHHwBYunQprVu3BuDXX3/N2W/evHk52zMzM3O+/+bNm9m8eTO9evXyevxWrVqxb98+9uzZA8CsWbNKFKenSj8WkNMJa9bA5Mnw/ffQsiVkHqoBa6EGZdvirpSytG3blnHjxnH55ZcTFBREdHQ077zzTrGOMX78eF566aWc9QMHDnjdLzg4mNdee43evXuTmZnJiBEjaNu2LQCTJk3ixhtv5B//+AfR0dE5VUbffvstAwcO5I8//mD+/Pn885//ZNs2q4KiW7du7Ny5kz///JMGDRowY8YMevfuzRtvvMGgQYNwOBzUqlWLt956C4BXXnmFefPmERwcTO3atXO+Z3p6ek6JpkaNGrz33ns57RGewsPDiY+P59prr6VKlSp069aNkydPFut6eSPGmL98kLISExNjSmNGsLUJqVxyRQRjmcB/Qsby1VdlX6xVyi47duzI+RWqKqaEhASmTJnCggUL8r3n7d9XRDYYY2I89y333UD9wYSGkU4w1fkTL728lFIqIFT6KiBvEr4SWlGNapwkPd3q1qYlAKXs8+yzzzJ79uw824YMGcK4ceNsisg+AwcOZO/evXm2TZo0id69e+esx8bGElsK3bECMgFERsJJqlOdk2RlWetKBRJjTLkaEXTcuHEBebP35rPPPivxZ4tbpR+QFSApKbkJQMRaVypQhIeHk5KSUuybhSrfsieECQ8P9/kzAV8CMEZLACqwNGjQgAMHDlCSoVVU+ZY9JaSvAjIBbNwITalGNf7MWVcqUISEhPg8ZaCq3AKyCghySwBKKRWoAjIBREfnTQDR0TYHpJRSNgjIBLBxI/xBLSJJyVlXSqlAE5AJ4NAh+JXzqcFJIjhNKY2sqpRSFUpAJgCAQ9QDoB6H8OM4VEopVW4FZAKoV88qAQCcz698/bUOCa2UCjwBmQCGDYPfJbcEoENCK6UCUUAmAKcTGnax5uus75qFUtsBlFKBxpYEICJDRGSbiGSJSL4hSsvC8ZA6/ElVmvIjgLYDKKUCjl0lgK3A9cAqm85P6llhDxdyIdYMO6mpdkWilFL2sCUBGGN2GGO+t+Pc2Zo3J08CaN7czmiUUqrsBWQbAMCuXVYCaMqPCFl8953dESmlVNnyWwIQkWUistXLUqxZl0UkTkSSRCSpNEcvvOAC2EVzwjlLI35i507tCqqUCix+SwDGmCuNMe28LHOLeZx4Y0yMMSambt26pRbfmDGwlfYAtGcLxmhXUKVUYAnYKiCnE6R9OwA6sBmA7dvtjEgppcqWXd1AB4rIAcAJLBSRJXbE8UdGdfbSmPZsAUDnx1BKBRJbJoQxxnwGlHziy1ISFgab6ZCTAMLCbA5IKaXKUMBWAQGcPQtbaE9LvieUs/o0sFIqoAR0AmjZ0ioBBJNJW7Zx6BDEx9sdlVJKlY2ATgBjxkAyUQBEkQzAjBl2RqSUUmUnoBOA0wnSrBknqUY01rRgaWk2B6WUUmUkoBMAQM3aDjbRMScBhIbaHJBSSpWRgE8AzZvDRqLpyCaELGrWtDsipZQqGwGfAHbtshJAdf6kGbtZulSHhFBKBYaATwAXXGAlAIBoNuqQEEqpgBHwCWDMGNhGW9IIyWkH0OcBlFKBIOATgNMJbaNC2UbbnASgs4MppQJBwCcAsJ4I3ki0KwEY9u+3OyKllPI/TQBA3bpWAjiXw1zAL+zfrw3BSqnKTxMA0KZN3oZggMmT7YxIKaX8TxMAMGwYbKYjWUhOAvjlF5uDUkopP9MEgNUQ3O+W6uymWU4CiI21NyallPI3TQAup065NwTD2rU2B6SUUn5m14xgz4vIThHZLCKfiYjtAzB8/701MmgT9lGTP1i9WhuClVKVm10lgKVAO2NMB+AHYKxNceRo2TK3ITiKZIyBhAR7Y1JKKX+yJQEYY740xmS4VtcCDeyIw501N0DenkDHjtkZkVJK+Vd5aAMYASyyOwinE2q3Po9fOD8nAWgJQClVmfltUngRWQbU8/LWOGPMXNc+44AM4P1CjhMHxAE0atTID5HmatkSNu7IbQi+4AK/nk4ppWzltwRgjLmysPdF5HagL9DTGGMKOU48EA8QExNT4H6loUULqx2gN0sI5wxVq0b483RKKWUru3oB9QHGAP2MMaftiMGb5GQrAQSTSXu28MEH2hNIKVV52dUG8BpQHVgqIskiMs2mOPIYNEjnBlBKBQ6/VQEVxhjTzI7zFiUuDl7/bxOObTpH5wZQSlV65aEXULnSuImQTFROAlBKqcpKE4AXyUTRni04yLQ7FKWU8htNAF5soiNVOEMzdrNvn93RKKWUf2gC8FCvnpUAADqwmeRkiI+3OSillPIDTQAehg2D7bQhgyA6sgmAGTNsDkoppfxAE4AHpxNadQxnJ61yEkB4uM1BKaWUH2gC8MLptKqBshNA7do2B6SUUn6gCcCL6GjYTAca8TO1OMr8+fpEsFKq8tEE4EVKSt6G4MzM3CeCExNh4kRNCEqpik8TgBexsdazAACd+A6A5cth9Gi4/HJ48km49FJ4/HEbg1RKqb9IE4AXTieYc+uxmwuJJQGAXbtg2jRokP4jg/iE+uZnJk+GW2+1N1allCopTQAFuOoqWE5PLucrgrAmL3uYF/iBFnzCEHbRnCF8zPvv63MCSqmKSRNAAdq2hWVcyTmcoDuruIepvMDfmUt/uvIN33IRMxlGG7bx/PPaJqCUqnikkLlYyp2YmBiTlJRUJudKTIQeXc/wMw0RDJEc5XP6M5hPyCSYc/mNrbRjEx25imVERFjtBE5nmYSnlFI+E5ENxpgYz+1aAiiA0wn1GkdwB2/zB7V4gzsZykdkukbQ/p3zmMhYrmQ5TtaQmqpzCCulKhZNAIUYOxYWcB3N2U0cb5BGGAAO11WbziiOcQ53Mw1jIDLSxmCVUqqY7JoS8hkR2eyaDexLESmX06/HxcH06dYAce5ef93adpqqfMRQBjGHapwkJcWeOJVSqiTsKgE8b4zpYIyJAhYAT9kUR5Hi4uDXX61E0KuX9TcuDv71L+v9mQyjKqe5loVaBaSUqlBsSQDGmBNuq1WBct8SHRcHS5ZYf7PXBwyAtVzCYepwLQv58kt9OEwpVXHY1gYgIs+KyM/ALZTjEkBhxowBI0Es4mquZhEOMvn0U7ujUkop3/gtAYjIMhHZ6mXpD2CMGWeMaQi8D9xXyHHiRCRJRJIOHz7sr3BLxOm0HhhbxNXUIYUokmna1O6olFLKN35LAMaYK40x7bwscz12fR8YVMhx4o0xMcaYmLp16/or3BKLjYXVdAPgMr5m5Up9KEwpVTHY1Quoudtqf2CnHXGUhthY+C24AXtpTDdWk5GhzwMopSoGu9oAnnNVB20GegEP2hTHX+Z0wiOPwNdcxmV8jTFGnwdQSlUIwXac1BhTYJVPRXTiBCTi5DbeoyE/s3FjI7tDUkqpIumTwKVkI9EARLORQ4dsDkYppXygCaAUDBsG24M6kIUQzUYWLdKGYKVU+acJoBQ4nXDzXVX5npZEs5H0dG0IVkqVf5oASkl0tFUNFM1GsrIKHhhO5xRWSpUXtjQCV0YbN0JN2nMzs6jGSTZurJ5vn/h4a17hrCwICYGvvtL5A5RS9tESQCnaQWsAWrEzX0NwfDz8fdRJ/pX1D97mdqLS19O7t04nqZSyjyaAUjJsGOwJbgVAa3bkaQhOTIT7R2fwKdfzJBPoz1xW0Z1WJ9czapQOIKeUsocmgFLidEL3Oy4knWBasTNPQ3BCAtyeNYOrWEYc8TRnF79yPh9wM8GkM2WKtgkopcpeiRKAiLQQkTdKO5iKrmNMCLtoTmt25GkI3rElgyd4jrV0YQYjSaEO9/MqzdjD7bxDVhbMnGlv7EqpwFNoAhCRDq4Zu7aKyHgROV9E5gArgO1lE2LFkZICO2lFa3bgcFjr8fHw66yVNGEfU3gUEM4/HxZyLRvoxH28BsCMGVoKUEqVraJKAG8AH2CN1nkYSAb2AM2MMS/6ObYKJzYWvpfWXMgeQiWd2FiYMwduYhbHqcEC+gLQpQsEBwszGElHNtORZB1ETilV5opKAGHGmHeMMd8bY14GThljxhhjUssiuIpmyxbYZloTQgaNM3ezZQsMut7QlwUsoC9nCScoyJpIZupU+MQxlDRCuIlZBAdbCUQppcpKUQkgXESiRaSTiHQCznqsKzdz5uR2BW3DdubMgYurbuNcDrOMK3E44L//tRqM4+Jg/OuRfMOl9GExWVlWAlFKqbJSVAL4FXgB+I9rOeS2PsW/oVU8gwZZCSALoS3bGDQIfpi2AoAV9EDEahfIlpICS6QPHdnMuZm/cN992g6glCo7hT4JbIy5oqwCqQysCeOrcOihJgxvu53f28Nv36xgD035ib8R5lHNExsL9zn68FzmE1zFUt7LHE5Cgj4drJQqGz51AxWR670sPUXkXH8HWNHExUF4pzZU/3kbUyZlEksCK+gBwNVX5725O50w6rX2HOMcupKo7QBKqTLl63MAI4E3gVtcyxvA48A3InJbSU8uIn8XESMidUp6jPImMRHiv2lLzd++59C89dTkOCspuCDVvqODdXShC2vJzNR2AKVU2fE1AQQDrY0xg1yzebUBDNAFKxEUm4g0xJoO8qeSfL68SkiAbbQllHTuM68AsIwrAahXz/v+67iE9mwhPPNPbQdQSpUZXxNAQ2PMb27rv7u2HQXSS3juF4ExWImk0oiMhHV0AeAmPuRbLuKInEtoqDVekKfYWFjvuIQgsoghiczMgp8HiI9HB5BTSpUaXxNAgogsEJHhIjIcmOfaVhU4VtyTikh/4KAxZlNxP1vebdwIu2jOPv4GwNd1BzJqFAU27jqdMPQ5azrJKNlEWJj3doDHH4epozYR8uUCHhx1RgeQU0r9Zb7OB3AvcD1wmWv9f8AcY4wB7xXcIrIM8FLpwTjgSazqnyKJSBwQB9CoUUWZbF0YwOcM5SPGH36AzLe9//rPdtuj53Hmmbr0rrGFNk/lTxSJiSCTJ7GJJwDYQSuumJzAgAHnaY8hpVSJiXUP92FHkfOAi7GqbNYbY34v0QlF2gPLgdOuTQ2AX4CLjTGFTqceExNjkpKSSnLaMpOYCN26QWZm3u133w2vv17wZ0537Uk1/uSKiHUsX543CTx++VomrurKbIYwmyH8j+F8w6U8220JX60S/30ZpVSlICIbjDExntt97QZ6A7AeGAzcAKwTkcElCcQYs8UYc64xprExpjFwAOhU1M2/onA6oUWL4n0mIQG20p52bCUtNStPG0B8PFy56v84TF1GMoM5DOYJnqMXS4lYvUTbA5RSJeZrG8A44CJjzHBjzDCsksD/+S+siq1ly7zrIoVXAUVGwmbaU5XTNDY/5gwjnZgIk0ft5iqW8TIPcopqAExnFD/RkEeZwowZfvoSSqlKz9cE4PCo8kkpxmcL5SoJHCmNY5UXY8ZYc/4COBwwbVrhT/empMBW2gHQTrbnDBeRkAC38j5ZCP9jOA4H9OoF6YTyNnfQgxUcWP+LNggrpUrE15v4YhFZIiK3i8jtwELgC/+FVbE5ndaE7xMmwNdfZw8RUbDYWNgf7ppOUnbklAAiI+E65pGIk1+oz6OPwpIlMGAAfMDNODAM5UMmT9auoUqp4vMpARhjHgPigQ6uJd4Yo787C+F0wtixvo3r43TCv18+h4NcQMusHTz0kFX9833iUTrzHYu4GocData09h8zBnY7WpJEZ27kQ8AaiVQppYrD52ocY8wcY8wjruUzfwYViFJSrJFEW7ODtDRrisjd760F4GsuIyQk9/kApxMefRTmcx0xJFGTP4iKsi92pVTFVNSUkCdF5ISX5aSInCirIANBbCzsCrISQGiI1TX3oow1ZBBEEhdxxx15SxOTJkG9m3rgwHA5q3j1VR1CQilVPIUmAGNMdWNMDS9LdWNMjbIKMhA4ndD7odbU4CSrP/qF6GhwksgmOnKKqkRH5//MydYXc5oIrmAFZ8/qlJJKqeIplZ48qnScaWzNJhaxdztHf8+gC+tYQ9ecCeY91TwvjG+4lCtYSVYWOY3HSinlC00A5URiIvS830oAMx7bQfOzW6nGKdaJs8DxgVJSYCVX0IEtRMpRr0kCrB5CXbrA5ZfD6NFaVaSUsmgCKCcSEuCI4zz+oCbNM3aQ9pV1l/7bTV3zDQ2RLTYWkkOtkUdjZIPXEkB8PCSM+oBZ65vyxqoWHJ72CZdfrklAKaUJoNyIjYXQMGGnqydQxuo1/Eo9Xvz0bwV+xumEG57rBEC02ZDTfdTd9v8s4gNu4Si1OUVVPmEIV6Z/weTJfvwySqkKQRNAOeF0wvLlkNWqNe2DdtCVNayhK2fTpNDG3V9Ta7GHpnQyG0hLy9sQvHZ1Og/+cA9baUs3VtOVNWyiAzMYydLPT+nDY0oFOE0A5cycnW2pnfE7F/IjX3F5kY27sbGQHNSZzmxAJO++6x6aRRP2MYbJpBJBKhHcw385n0OM5r86+5hSAU4TQDmSkADz6JezPo9+BfYAyuZ0QoN+nWnKXqpnHM2pBoqPh+jv3mQnLVnE1YSEWE8Qr3VcyjJ6cj+vkpWeycyZ/v9eSqnySRNAORIZCXtoxiP8h6F8yH4a53kCuCA/1uoMQGdyq4FWv7ef7qzmHW4HhOho6+Gx11+H6Y7RNOJnerOYt9/WUoBSgUoTQDmS/Uv/RR7hY4YC5HsC2JtmN1gNwZ35jqAgK2FcH74QgDkMAmDkSGvfuDg4d2Q/jlKLIcwmPV0fIFMqUGkCKEdiYyHYY5JOb08Ae8qoUZu9NCaa7xCBLVsgdOlCdtGMHx3NGTMm74ikHWNCWEBfrmM+kpWhD5ApFaA0AZQjTifceWfuelH1/9kSEmADnXOqgN574ww9WMFCrkUkdxTRbCkpME8GEMlRusvXPp1DKVX52JIARORpETkoIsmu5Ro74iiPhg2DiAgICqLAJ4A9xcbCpqBONGMPNcwxzvluJRGksliuITQ0/zFiYyEhrDdnCGcAnxdaAkhMhIkTtZ1AqcrIzhLAi8aYKNeik8u4ZD8P8MwzFPgEsLfP1O1jNQRHs5E+WQs5RRVWyeW89FL+YzidMOHlqizlKq4zc7n3HuP1mYD4eGv4iHHjrInuBw7URKBUZaJVQOVQcSaTyfnMPVZD8CWs5Rq+YBlXciYrjI0bve+fkgJfcA1N2EfTzB/yPROQmAhTR21icXoPfjINmJo5ipWfH9NhJJSqROxMAPeJyGYReUtEatkYR6WQUasu3xLDRJ6kCftYyLWF7h8bC8uDegPQmyVkZubtDbTl4x0kEEsbtvMNlzKCt1hCb4LTT+swEkpVEn5LACKyTES2eln6A68DFwJRwK/Afwo5TpyIJIlI0uHDh/0VboWXkAAfcwMAp6jCLG4iJMRqU/DG6YTH/tuEXTSnN0twONyeIs7MZOCnt5FGKJewlhv5iCHMpgvreZZxzJ+vpQClKgVjjK0L0BjY6su+nTt3Nsq7NWuMiQg6ax7gJdOHL4yIMXffXfTntsTeZ04RYcI4YyIirOMkDJthDJgb+cAEBxtzyy3GiBjzGveYTMR0kM1mwgT/fyelVOkAkoyXe6pdvYDOd1sdCGy1I47KxOmEl/4byushD/Kl42rCwwv+9e9uR6PeVOEM3VhtdSF9O51GM5/hW2L4kBvJyoK2bWHaNPh38DOcoAbP8A99dkCpSiC46F38YrKIRAEG2AeMsimOSiUuDtq3t6qDYmN9a0RuNPwKTs2swhD5hG9Cr6JV0ns0YR/38yogiLgfqzYvjf47T2c9Rad7NgMd8jxgppSqWMQqHVQMMTExJikpye4wKp3DvW+h2upFbJm3j4YDOnHwVE0u4ltAGDAAPvvM2m/iRJj85B/8RCPm0p8RIe/x1VfF662klCp7IrLBGBPjuV27gSp+7zuSiDN/cPFV53D+qT2MZSIgOSOIZouNhT+DaxFPHDfyIfUz9hc4jpA+QKZU+acJIMAlJkK7B3ownTjOEsoExrKMqxCxBpBz/3XvdMLUqTA1+CEMwsO8wPr1eW/y8fHQpg3EXprO/z2ZqQ+QKVWOaQIIcNm/4O9mOhGcYRwTAGtQOm+NyHFx8MTUhrzPLYw0b7Lq8xS6dbNu/PHxMGXUD0zZcQ1nTBgnqcZrmXez/PMTOfsopcoPTQABLjYWQkOt19mdwoKC4LXXCq7bT0mByYyhKqd5jOfJzIS774aEMV+wnovpwjqm8Cjvcwt38iar6cY5mSmMHq1JQKnyxK5eQKqccDqtUsDMmXDoENSrZ/3yL6xhNzYWng5twztpw3mM59lINM3NLv59/Ck20ZGBfMZ+GgMwmyHMox+fMJheWV8yenQIgPYeUqoc0F5AqkQSE+HVZ08wduGltHc9xvEhQxnBW5yhCt27wy23wKJFUP3zd5nJMMYzjv9jPCEhFNh7KDGxeN1YlVJF015AqlQ5nfDBghqsm7qBwY45dGEtNzGLM1TB4YA+faxf+Z99BpdNv423uYMnmUA3VuUbdyhbfDx07QpPPgndu2t1kVL+pglA/SV33hPK37++ngsGdCEoSHA48s9jEBcH5uVX2MOFvMet1JUjOU8SZ3cXfWNaJs+O2k9trNlpMjLIN0KpUqp0aQJQf5nTaf3SX70axo/3Po/BiAeqsWXsLM7ldz7N7MfTD/7B44/DgK6/cfbJp+k/+nz205gU6rCay4jmO87T2OkAABKfSURBVNLT4aGHNAko5S/aBqDKzMSJkPTkHD5iKH9Qix9owcWsJ4QM5tOXefTjXH7nfl6lFn9wAx8zj/6EhcHKldomoFRJaRuAsl1sLCyKGEQX+ZavuByD8AKP0JKd9GM+b3IXExjH9S22sZFo5jCInizj7Fl4+unCSwL65LFSxaclAFWmEhOtm/mXX3p/PywMXnkFnrzvBCvTL6UhP+MkkZ20JiLCe/VSYqLVeAwUuI9SgUxLAKpccDqtBJAtONgab+juu61l5Uqr0Xj+VzX48JYFZAWHsYC+nMtvpKZazyvk2LePg7c9QcaAwbzGvVzK1z6VFpRSFi0BKFv42t9/y5vruPCuK/iZhvRnLvvCWvH+o9/Rev7ztNo6m8wsYRfN+Rv7qcpp4rmLB3mZjOAIpk7VB86UgoJLAJoAVLk3ZcDX3Dm3LzU4wWHqch6/c4LqvCGjeNE8yEEaEMFpnuIZxjCJpVxFf+aSJuH072+VMLRKSAUyTQCqwkpMhJt7HGJY6nTqc5DttOFt7uAE5+TZTwTu4G1mmBF8wdUM5DPSCMvpRQRWqSMy0hrPyL30IWL9XbNGk4WqfApKADoWkCr3nE4Y+3I9Ro36Z6H7hYdDnfvv4M7JmbzJXXzCYAYxh7S0UGbOtKa1BAgjlRv4mNqSiIlyMGVjT8K4hrOE07OnNiKrwGFbI7CI3C8iO0Vkm4hMtisOVTGkpORdz/7FDuBwQK9e1o170iS4ePqdjJZpXMcCZjOE80KOMm0aNOQnnuVJfqYhMxnODeZD2m18l08ZxFba0Y1VpKZaD5+NHp2/IVnEWhIT8772lXZVVeWOt5ni/b0AVwDLgDDX+rm+fK5z586FTXyvKrE1a4wBawkONmbMGGMiIowJCrL+rlmTf/9Zl71mMhFzgmpmG61NBg6TgcN8ygDTg2UGskwQ6eZa5pvdNDWZiHmUyQayDBgTFpZ7XPfzh4bmvvZ27qLi9/UzSpUWIMl4uafaVQIYDTxnjDnrSkK/2xSHqiCcTqt+fsIEWLXK+qW/fDk884z3KhunE25cfS+b393MRwxlG22ZyFgu5Ede6v4Zq0N6IiJkEsxC+tKBzXzCYJ5nDG9zB6GcJS0td9A698HrwtJOcg7HCCKDM2esZxCKKg24f979uErZyZZGYBFJBuYCfYBU4FFjzLcF7BsHxAE0atSo8/79+8ssTlXxTZxojS6aLXsoashtEN640VqP7pjFwdHP8C+eZj0X8UTICzybcCnOzmnEhicymE/oxZe0YBcAf1KVL+nFNO5mKVcRFiYFDlnh/rBatrJscM6uMqtAfT5UKSqoEdif1TzLgK1elv6uv68CAlwM7MWVjApbtApIFZdn1dH06UXvP/2q2eZYlXrGgDlNuEnFqvM5RYSZL33N+v7jzd4HXjALGowyv2Dt9zVdTU+WmQnPZhV67LvvtqcqKPucWvUUmCigCshvvYCMMVcW9J6IjAY+dQW2XkSygDrAYX/FowJTdtWRr5PMOJ3g/HIwnL6GPf96j08n7yILB2voynJ6kuqoRuZca9/p0+H6US8zgrcYx7Ms40pWjevG9mpPEfVgd9Kx5tqsyR/8kfQjzgM/cuTIGX6jGttpw56zLUhIcBQZ01/poureWA5oLycblcdSmF1VQHcDFxhjnhKRFsByoJEpIhh9DkCVJc+bJ1jzJ6elWa8jIuCll6wqpJSDqZw3/w3GMpEL+JWzhJJCJNU5SXX+9Hr8o9Qipcu1NH9sgDWDTtWqhZ6/sHGOPPc1xnu1U1CQ1W4ydmxh39xe7jfKwm6a5fGGCgXHb+ezJuXtOYC3gLdEZCuQBgwv6uavVFnybNANCYGRI63X2c8TpKVZ3VNffx0mTgznyfn38yZ3ch3z6cwGanOUk1TnAA34kaYcDGnC9PersXfTcRY8u4lYEui7bgEMfo8zhPMtF9H9wU5w3nlQrRoPkkkQmQSTQRCZBKUaUl5oAA83p+6lzTlCHYzxkqXwnrzASmDuk/UUdBP1llBKQ2nctD1jc090xTmuvxOI+39D7q8LKoXZkdBsSQDGmDTgVjvOrZQv3HvpiFg3/9dft/5Hzk4A7jfT7L+pRDCbG5jNDfmOKRmwZDdQFd6hE+9wB0Fk0I3V9GMel7CWrNen40hLBeAlzwMY4BNrOQwc4xw2yIX8SFMmciGHqMcpqnKaKpymSp7X6SFVuW5oFaa/V4XLutYki6A8h05M9O1XqS83qdLax/MGWlB8PXsW77ilxdu5PG/02dwHMczuOVZYCcdzu7++lw4FoZQXhQ0xXdBAdgX96s7mPiSFZ9VMLkMEZ6jKKTIIJpMgsiSY2J5BLFtmaMRPNGcXzdlFM3ZzIXtoyo80YS+hpPv03bIQjlKbI9ThMHX5k2qkO8Lp2iOcOg3COXQ8nFmfhZNK7nKWsDzrT08Mp1VUOISHc8kVeffNXo6dCbe+tNuF8bxGntUh7tfdvbot5+qY/McICoLMzPz7eZ6zoCoZb59x58tNOvv1hAm5vc4cDsjK8u27FPa6sBh8Vd6qgJQq19xvSp7FdafT+6/RNWvy3rxGjIDo6NxupsOG5X5uzZqC5kUQzlCFM1TJ2RIRDv2GwMJlsIsW7KJFvnM7yKQGJ1y/909TlVNU4TTVHae56+ZTDL7mNHu2nOaViX9Sm6PU5TB1OUwdjhBJCuFZqRxdlkqdhqlUP3qWEa7beBhp+c4FgFsbwlrve0CE629YGISFkS4h/EIIGQST7vqbcWkIpy4Mpuo5IRAcTINDIaxwvZ+eFpKTBLOXI70czHRbz8KBEMRZ1+vs7YwJsu7AQUH8y/UezwQx1vV6/z3CYwjGbeEF4aG/OzAIL78s3P+gtf2e7PdfF+uYrsfAR2Zvf0u4HSELB4PPCNtd2x0IGa7Xki5k4Mhzvu8nCgOyjzFX6Ot6f8cU4ers7YutO38vhA10xuorU3q0BKBUKSpOQ5+3RlpvsnsxuT/P4Cv30ovnMxG+ELII46yX3/fWkrgilat7pBa6z9OPp/LSpFRCSCeYjDx/Q0mnfesMWjRJh4wMjh9JZ9N31nshkkGwSc+53TvI4tzITI6l5K4HkUnN6pmcOpnllhYyiQjLIuNs7n4OKs59riB9WMQ/1/QpUQOylgCUKgPF+T3l3kX12DGY7DEiVnFmN/Ps3ZOdiNw/79746yuDg1QiSM35OZ9XYjgsLuIYTz8HD0/y/l5EBJzZAezI/96ab+AijwQ5fQKMGpV324Sx+RNbhAPO5NliMOmZhIVYCSFhhaFHD5OnDFAl3JCaanCQxdEjhjp18r6/cL7huuus9+d9bhgwwNoeEWY4e9Z6vXd3VvYjF2AMLVta2z/8wBDV0fDG9CxeecXaFuwwZGVZr4cMMsyZYx07SKw++oIhLAweftAwaTLspFXpd+P19nBAeV30QTBVWU2Y4HbXKOChLff3IiKsh9qK+1CZ5zkKWnzd3/3BtsKOVdB77g/qeZ47m+c4Sp7Xx9sxgoIKj8NbTO6fcb+23s7t/u/l/rmC/s28fRf3MaXcXwcH5z12r1551ydMKPrfOf+/e/kaC0gp5cb913lERNFVSNldULP5+qsw+1aSLXt8JU++jlj61lu+7ed+Pvfz+xKz5zhKnp/1dozQ0OLFBVYpKtu99+Z/3/3c7v9e7p/r2TPvtfO83u6xjhiR+9q9EXvq1NzXoaEwaFDe9ZKU5AqiCUCpcqCwRudsa9bkvs6+EWTfYIpbJeD+OW83lOwbWVEJw7P3TfY+7p/x7Brp7biFcY8vNDT/8b1Zvrzo43rG4H5Dzu694849qRR0Iy/OQH/Dhnk/tvs0psuX518vzYfINAEoVU4UdTP3JUmUhLfjFHQj80wY3n5pex7PvQ98SUZC9eV7F/ZL21fuN+SwsPzvuycV96TmeSMv6he6t3/nghKWt1FuS5MmAKUqIH8PJeB5I3O/wRZ043Lfx/0G6V5N5Hlcb1VNhVU/FfW9s2MoyaQ7vtyQs/XsmVsSKY3E7HQWr1RUWjQBKFWB+Osm4XlcX29kBd243H/lu1cTeR7XW2nAsx7dW3xFKU4pw9uxvX13X+Z0KK3E7BmTv/7dNQEopfL5qzeygqqJPI/rrbqkNCbMKY2GUs9+QJ7fqTQbY+2iCUApVWxF/SL1tTqloB48f/Xm6o8qMn+0wdhR7eNOHwRTSvlVdjWRr+ycr8DXOD3js/Mm/ldoAlBK5SgPNzKdrKbsaBWQUkoFKC0BKKX8wtfSRPZ+RQ2nXdLzl/ZxKxNbSgAi8pGIJLuWfSKSbEccSqnAVJLnBMD+RtvSZksCMMYMNcZEGWOigDnAp3bEoZQKHN6GpAh0trYBiIgANwCz7IxDKWU/f/+6/qtDUlRGdjcCdwN+M8bsKmgHEYkTkSQRSTp8+HAZhqaUquh8HZIiUPmtEVhElgH1vLw1zhgz1/X6Jor49W+MiQfiwZoRrFSDVEpVar4OSRGo/JYAjDFXFva+iAQD1wOd/RWDUiqweQ7fcMY1TZje/C12VgFdCew0xhywMQalVCVWnBE+A5GdCeBGtPFXKVVG9Fd/frY9CGaMud2ucyullLK/F5BSSimbaAJQSqkApQlAKaUClCYApZQKUJoAlFIqQGkCUEqpAKUJQCmlApQmAKWUClCaAJRSKkBpAlBKqQClcwIrpQJGZZrOsTRoCUAppQKUJgCllApQmgCUUipAaRuAUqpS03r/gmkJQCmlApQtCUBEokRkrYgki0iSiFxsRxxKKRXI7CoBTAb+ZYyJAp5yrSullCpDdiUAA9RwvT4H+MWmOJRSKmDZ1Qj8ELBERKZgJaGuBe0oInFAHECjRo3KJjqllAoAfksAIrIMqOflrXFAT+BhY8wcEbkBmAFc6e04xph4IB4gJiZG2/OVUqqU+C0BGGO83tABRGQm8KBrdTbwpr/iUEop5Z1dbQC/AJe7XvcAdtkUh1JKBSy72gDuAl4WkWAgFVcdv1JKqbIjpgI9Jicih4H9PuxaBzji53AqAr0OufRaWPQ65Aqka/E3Y0xdz40VKgH4SkSSjDExdsdhN70OufRaWPQ65NJroUNBKKVUwNIEoJRSAaqyJoB4uwMoJ/Q65NJrYdHrkCvgr0WlbANQSilVtMpaAlBKKVUETQBKKRWgKnQCEJE+IvK9iOwWkSe8vB8mIh+53l8nIo3LPkr/8+E6dBeR70QkQ0QG2xFjWfHhWjwiIttFZLOILBeRv9kRp7/5cB3uFpEtrjk5vhaRNnbEWRaKuhZu+w0SESMigdM11BhTIRcgCNgDNAVCgU1AG4997gGmuV7fCHxkd9w2XYfGQAdgJjDY7phtvhZXAFVcr0cH8H8TNdxe9wMW2x23XdfCtV91YBWwFoixO+6yWipyCeBiYLcx5kdjTBrwIdDfY5/+wP9crz8BeoqIlGGMZaHI62CM2WeM2Qxk2RFgGfLlWqw0xpx2ra4FGpRxjGXBl+twwm21KtYcHZWRL/cJgGeASVhD0wSMipwA6gM/u60fcG3zuo8xJgM4DkSWSXRlx5frECiKey1GAov8GpE9fLoOInKviOzBmpHvgTKKrawVeS1EpBPQ0BizsCwDKw8qcgJQqsRE5FYgBnje7ljsYoyZaoy5EHgc+Ifd8dhBRBzAC8Df7Y7FDhU5ARwEGrqtN3Bt87qPa+TRc4CUMomu7PhyHQKFT9dCRK7EmpionzHmbBnFVpaK+9/Eh8AAv0Zkn6KuRXWgHZAgIvuAS4B5gdIQXJETwLdAcxFpIiKhWI288zz2mQcMd70eDKwwrhafSsSX6xAoirwWIhINTMe6+f9uQ4xlwZfr0Nxt9Voq75wchV4LY8xxY0wdY0xjY0xjrHahfsaYJHvCLVsVNgG46vTvA5YAO4CPjTHbROTfItLPtdsMIFJEdgOPAAV2AauofLkOInKRiBwAhgDTRWSbfRH7j4//TTwPVANmu7pAVrpk6eN1uE9EtolIMtb/G8MLOFyF5uO1CFg6FIRSSgWoClsCUEop9ddoAlBKqQClCUAppQKUJgCllApQmgCUUipAaQJQ6i8QkQYiMldEdonIjyLymoiE2R2XUr7QBKBUCbkGFvwU+NwY0xxoDkRgja2jVLmnzwEoVUIi0hP4pzGmu9u2GsB+rMHF/rQtOKV8oCUApUquLbDBfYNrmOV9QDM7AlKqODQBKKVUgNIEoFTJbQc6u29wVQHVA763JSKlikETgFIltxyoIiLDAEQkCPgP8Jox5oytkSnlA00ASpWQa2jxgcBgEdmFNddEljHmWXsjU8o32gtIqVIiIl2BWcBAY8x3dsejVFE0ASilVIDSKiCllApQmgCUUipAaQJQSqkApQlAKaUClCYApZQKUJoAlFIqQP0/VxYW3sRw7oUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "objective.plot()\n",
    "plt.legend()\n",
    "plt.xlabel('Q')\n",
    "plt.ylabel('logR')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Let's see the results of the fit. For the case of DifferentialEvolution uncertainties are estimated by estimating the Hessian/Covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________________________________________________\n",
      "Objective - 4829712976\n",
      "Dataset = c_PLP0011859_q\n",
      "datapoints = 408\n",
      "chi2 = 921.484228321403\n",
      "Weighted = True\n",
      "Transform = Transform('logY')\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "________________________________________________________________________________\n",
      "Parameters: 'instrument parameters'\n",
      "<Parameter:    'scale'    , value=0.8797 +/- 0.00305, bounds=[0.6, 1.2]>\n",
      "<Parameter:     'bkg'     , value=4.42652e-07 +/- 2.13e-08, bounds=[1e-09, 9e-06]>\n",
      "<Parameter:'dq - resolution', value=5 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters: 'Structure - ' \n",
      "________________________________________________________________________________\n",
      "Parameters:      'Si'      \n",
      "<Parameter: 'Si - thick'  , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  'Si - sld'   , value=2.07 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  'Si - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter: 'Si - rough'  , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'Si - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters:     'SiO2'     \n",
      "<Parameter:'SiO2 - thick' , value=38.9178 +/- 0.377, bounds=[15.0, 50.0]>\n",
      "<Parameter: 'SiO2 - sld'  , value=3.47 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter: 'SiO2 - isld' , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'SiO2 - rough' , value=5.85997 +/- 0.283, bounds=[1.0, 15.0]>\n",
      "<Parameter:'SiO2 - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters:     'film'     \n",
      "<Parameter:'film - thick' , value=258.835 +/- 0.252, bounds=[200.0, 300.0]>\n",
      "<Parameter: 'film - sld'  , value=2.40753 +/- 0.0129, bounds=[0.1, 3.0]>\n",
      "<Parameter: 'film - isld' , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'film - rough' , value=8.97212 +/- 0.37 , bounds=[1.0, 15.0]>\n",
      "<Parameter:'film - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters:      'd2o'     \n",
      "<Parameter: 'd2o - thick' , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  'd2o - sld'  , value=6.36 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter: 'd2o - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter: 'd2o - rough' , value=3.69552 +/- 0.112, bounds=[1.0, 15.0]>\n",
      "<Parameter:'d2o - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n"
     ]
    }
   ],
   "source": [
    "print(objective)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now lets do a MCMC sampling of the curvefitting system. First we do 400 samples which we then discard (burn). These samples are discarded because the initial chain might not be representative of an equilibrated system (i.e. distributed around the mean with the correct covariance)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 400/400 [01:15<00:00,  5.30it/s]\n"
     ]
    }
   ],
   "source": [
    "fitter.sample(400, pool=1)\n",
    "fitter.reset()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We then follow up with a production run, only saving 1 in 100 samples. This is to remove autocorrelation. We save 15 steps, giving a total of 15 * 200 samples (200 walkers is the default)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 90%|█████████ | 1356/1500 [04:12<00:26,  5.52it/s]"
     ]
    }
   ],
   "source": [
    "res = fitter.sample(15, nthin=100, pool=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In the final output of the sampling each varying parameter is given a set of statistics. `Parameter.value` is the median of the chain samples. `Parameter.stderr` is half the [15, 85] percentile, representing a standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "print(objective)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "A corner plot shows the covariance between parameters. You need to install the *matplotlib* and *corner* packages to create these graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "objective.corner();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Once we've done the sampling we can look at the variation in the model at describing the data. In this example there isn't much spread."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "objective.plot(samples=300);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In a similar manner we can look at the spread in SLD profiles consistent with the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "structure.plot(samples=300)\n",
    "plt.ylim(2.2, 6);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Fitting the BornAgain example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "[BornAgain](https://www.bornagainproject.org/) is another program for fitting reflectometry and GISAS data. The following cells repeat the analysis in their [specular fitting example]( https://github.com/scgmlz/BornAgain/blob/master/Examples/python/fitting/ex03_ExtendedExamples/specular/FitSpecularBasics.py). The simulated thickness of the Titanium layer is 30 Angstrom, but we'll start the fit with a value of 50."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# necessary imports\n",
    "import numpy as np\n",
    "from refnx.util import q\n",
    "from refnx.analysis import Objective, CurveFitter, Transform\n",
    "from refnx.dataset import Data1D\n",
    "from refnx.reflect import SLD, Stack, ReflectModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# first grab the data from the BornAgain repository.\n",
    "# The data was originally created in genx. \n",
    "import requests as req\n",
    "import io\n",
    "import gzip\n",
    "url = (\"https://github.com/scgmlz/BornAgain/raw/master\"\n",
    "       \"/Examples/python/fitting/ex03_ExtendedExamples\"\n",
    "       \"/specular/genx_interchanging_layers.dat.gz\")\n",
    "f = gzip.open(io.BytesIO(req.get(url).content))\n",
    "dataset = np.loadtxt(f, usecols=(0, 1), skiprows=3).T\n",
    "\n",
    "# data is saved as two_theta, convert to Q.\n",
    "dataset[0] = q(dataset[0] / 2, 1.54)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# make a dataset\n",
    "data = Data1D(data=dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# make the structure\n",
    "air = SLD(0)\n",
    "si = SLD(2.0704) # silicon substrate\n",
    "ni = SLD(9.4245) # nickel\n",
    "ti = SLD(-1.9493) # titanium\n",
    "\n",
    "# make the layers\n",
    "ti_layer = ti(50)\n",
    "ni_layer = ni(70)\n",
    "\n",
    "# Make a multilayer by using a Stack Component\n",
    "stack = Stack()\n",
    "stack |= ti_layer\n",
    "stack |= ni_layer\n",
    "stack.repeats.value = 10\n",
    "\n",
    "structure = air | stack | si(0, 0)\n",
    "\n",
    "# put the Structure in a Model\n",
    "model = ReflectModel(structure, bkg=0, dq=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# we're only going to fit the Titanium thickness\n",
    "ti_layer.thick.setp(vary=True, bounds=(10, 60))\n",
    "\n",
    "# now do the fit\n",
    "objective = Objective(model, data, transform=Transform('logY'))\n",
    "fitter = CurveFitter(objective)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# now fit and plot\n",
    "fitter.fit('differential_evolution');\n",
    "fig, ax = objective.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "print(objective.chisqr())\n",
    "print(ti_layer)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}