{
  "cells": [
  {
   "cell_type": "markdown",
   "id": "ffdfe095",
   "metadata": {},
   "source": [
    "# SciPy minimizers and constraints\n",
    "\n",
    "The `Minuit` class can call SciPy minimizers implemented in `scipy.optimize.minimize` as alternatives to the standard Migrad minimizer to minimize the cost function. The SciPy minimizers may perform better or worse on some problems. You can give them a try when you are not satisfied with Migrad.\n",
    "\n",
    "More importantly, the SciPy minimizers support additional features that Migrad lacks.\n",
    "\n",
    "* Migrad does not allow one to use an externally computed hessian matrix.\n",
    "* Migrad does not allow one to use additional constraints of the form $\\vec a \\leq f(\\vec x) \\leq \\vec b$ in the minimization, where $\\vec x$ is the parameter vector of length $m$, $f$ is an arbitrary function $\\mathcal{R}^m \\rightarrow \\mathcal{R}^k$ and $\\vec a, \\vec b$ are vector bounds with length $k$.\n",
    "\n",
    "SciPy comes with a variety of minimization algorithms and some of them support these features. The ability to use constraints is interesting for HEP applications. In particular, it allows us to ensure that a pdf as a function of the parameters is always positive. This can be ensured sometimes with suitable limits on the parameters, but not always.\n",
    "\n",
    "We demonstrate this on a common example of fit of an additive model with a signal and background pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "unnecessary-vermont",
   "metadata": {},
   "outputs": [],
   "source": [
    "from iminuit import Minuit\n",
    "from iminuit.cost import ExtendedUnbinnedNLL\n",
    "import numpy as np\n",
    "from numba_stats import norm, bernstein\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display\n",
    "import joblib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bac1095",
   "metadata": {},
   "source": [
    "The signal pdf is a Gaussian, the background is modelled with second degree Bernstein polynomials. We perform an extended maximum likelihood fit, where the full density model is given by the sum of the signal and background component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "equivalent-minnesota",
   "metadata": {},
   "outputs": [],
   "source": [
    "xrange = (0, 1)\n",
    "\n",
    "\n",
    "def model(x, b0, b1, b2, sig, mu, sigma):\n",
    "    beta = [b0, b1, b2]\n",
    "    bint = np.diff(bernstein.integral(xrange, beta, *xrange))\n",
    "    sint = sig * np.diff(norm.cdf(xrange, mu, sigma))[0]\n",
    "    return bint + sint, bernstein.density(x, beta, *xrange) + sig * norm.pdf(x, mu, sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58a91a58",
   "metadata": {},
   "source": [
    "In searches for rare decays, it is common to fit models like this to small simulated samples that contain only background, to calculate the distribution of some test statistic (usually the likelihood ratio of S+B and B-only hypotheses). Here, for simplicity, we use the signal amplitude itself as the test statistic.\n",
    "\n",
    "We run one such fit. The mean and width of the Gaussian are fixed, only the signal amplitude and the background parameters are varied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "abroad-census",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <th colspan=\"5\" style=\"text-align:center\" title=\"Minimizer\"> Migrad </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:left\" title=\"Minimum value of function\"> FCN = -181.3 </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center\" title=\"Total number of function and (optional) gradient evaluations\"> Nfcn = 96 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:left\" title=\"Estimated distance to minimum and goal\"> EDM = 1.83e-06 (Goal: 0.0002) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center\" title=\"Total run time of algorithms\">  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Valid Minimum </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center;background-color:#92CCA6;color:black\"> No Parameters at limit </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Below EDM threshold (goal x 10) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Below call limit </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\"> Covariance </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\"> Hesse ok </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Is covariance matrix accurate?\"> Accurate </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Is covariance matrix positive definite?\"> Pos. def. </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Was positive definiteness enforced by Minuit?\"> Not forced </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th title=\"Variable name\"> Name </th>\n",
       "        <th title=\"Value of parameter\"> Value </th>\n",
       "        <th title=\"Hesse error\"> Hesse Error </th>\n",
       "        <th title=\"Minos lower error\"> Minos Error- </th>\n",
       "        <th title=\"Minos upper error\"> Minos Error+ </th>\n",
       "        <th title=\"Lower limit of the parameter\"> Limit- </th>\n",
       "        <th title=\"Upper limit of the parameter\"> Limit+ </th>\n",
       "        <th title=\"Is the parameter fixed in the fit\"> Fixed </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 0 </th>\n",
       "        <td> b0 </td>\n",
       "        <td> 16 </td>\n",
       "        <td> 17 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td> 0 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 1 </th>\n",
       "        <td> b1 </td>\n",
       "        <td> 85 </td>\n",
       "        <td> 35 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td> 0 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 2 </th>\n",
       "        <td> b2 </td>\n",
       "        <td> 16 </td>\n",
       "        <td> 16 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td> 0 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 3 </th>\n",
       "        <td> sig </td>\n",
       "        <td> -4.0 </td>\n",
       "        <td> 2.7 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 4 </th>\n",
       "        <td> mu </td>\n",
       "        <td> 0.500 </td>\n",
       "        <td> 0.005 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td> yes </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 5 </th>\n",
       "        <td> sigma </td>\n",
       "        <td> 50.0e-3 </td>\n",
       "        <td> 0.5e-3 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td> yes </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th> b0 </th>\n",
       "        <th> b1 </th>\n",
       "        <th> b2 </th>\n",
       "        <th> sig </th>\n",
       "        <th> mu </th>\n",
       "        <th> sigma </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> b0 </th>\n",
       "        <td> 306 </td>\n",
       "        <td style=\"background-color:rgb(162,162,250);color:black\"> -412 <strong>(-0.676)</strong> </td>\n",
       "        <td style=\"background-color:rgb(250,192,192);color:black\"> 110 <strong>(0.387)</strong> </td>\n",
       "        <td style=\"background-color:rgb(250,202,202);color:black\"> 15.1 <strong>(0.319)</strong> </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> b1 </th>\n",
       "        <td style=\"background-color:rgb(162,162,250);color:black\"> -412 <strong>(-0.676)</strong> </td>\n",
       "        <td> 1.21e+03 </td>\n",
       "        <td style=\"background-color:rgb(166,166,250);color:black\"> -366 <strong>(-0.645)</strong> </td>\n",
       "        <td style=\"background-color:rgb(164,164,250);color:black\"> -62 <strong>(-0.659)</strong> </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> b2 </th>\n",
       "        <td style=\"background-color:rgb(250,192,192);color:black\"> 110 <strong>(0.387)</strong> </td>\n",
       "        <td style=\"background-color:rgb(166,166,250);color:black\"> -366 <strong>(-0.645)</strong> </td>\n",
       "        <td> 266 </td>\n",
       "        <td style=\"background-color:rgb(250,204,204);color:black\"> 13.4 <strong>(0.305)</strong> </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> sig </th>\n",
       "        <td style=\"background-color:rgb(250,202,202);color:black\"> 15.1 <strong>(0.319)</strong> </td>\n",
       "        <td style=\"background-color:rgb(164,164,250);color:black\"> -62 <strong>(-0.659)</strong> </td>\n",
       "        <td style=\"background-color:rgb(250,204,204);color:black\"> 13.4 <strong>(0.305)</strong> </td>\n",
       "        <td> 7.28 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> mu </th>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> sigma </th>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td style=\"background-color:rgb(250,250,250);color:black\"> 0 </td>\n",
       "        <td> 0 </td>\n",
       "    </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "┌─────────────────────────────────────────────────────────────────────────┐\n",
       "│                                Migrad                                   │\n",
       "├──────────────────────────────────┬──────────────────────────────────────┤\n",
       "│ FCN = -181.3                     │              Nfcn = 96               │\n",
       "│ EDM = 1.83e-06 (Goal: 0.0002)    │                                      │\n",
       "├──────────────────────────────────┼──────────────────────────────────────┤\n",
       "│          Valid Minimum           │        No Parameters at limit        │\n",
       "├──────────────────────────────────┼──────────────────────────────────────┤\n",
       "│ Below EDM threshold (goal x 10)  │           Below call limit           │\n",
       "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n",
       "│  Covariance   │     Hesse ok     │ Accurate  │  Pos. def.  │ Not forced │\n",
       "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n",
       "┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n",
       "│   │ Name  │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │\n",
       "├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n",
       "│ 0 │ b0    │    16     │    17     │            │            │    0    │         │       │\n",
       "│ 1 │ b1    │    85     │    35     │            │            │    0    │         │       │\n",
       "│ 2 │ b2    │    16     │    16     │            │            │    0    │         │       │\n",
       "│ 3 │ sig   │   -4.0    │    2.7    │            │            │         │         │       │\n",
       "│ 4 │ mu    │   0.500   │   0.005   │            │            │         │         │  yes  │\n",
       "│ 5 │ sigma │  50.0e-3  │  0.5e-3   │            │            │         │         │  yes  │\n",
       "└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n",
       "┌───────┬───────────────────────────────────────────────────────┐\n",
       "│       │       b0       b1       b2      sig       mu    sigma │\n",
       "├───────┼───────────────────────────────────────────────────────┤\n",
       "│    b0 │      306     -412      110     15.1        0        0 │\n",
       "│    b1 │     -412 1.21e+03     -366      -62        0        0 │\n",
       "│    b2 │      110     -366      266     13.4        0        0 │\n",
       "│   sig │     15.1      -62     13.4     7.28        0        0 │\n",
       "│    mu │        0        0        0        0        0        0 │\n",
       "│ sigma │        0        0        0        0        0        0 │\n",
       "└───────┴───────────────────────────────────────────────────────┘"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyUUlEQVR4nO3deXhTVf7H8c9N0gaQtsgOUkR2WWURBERWF0QWl0HQQcQFddBRGRURRxQXGH8uqKPoIIqjbLKqgCKyiCwqq7IrqyC7YIsF0iY5vz+uVhlFSDlJmvJ+Pc992tzc9n5zmyafnHvuOY4xxggAAMACT7wLAAAABQfBAgAAWEOwAAAA1hAsAACANQQLAABgDcECAABYQ7AAAADWECwAAIA1vljvMBwOa+fOnUpJSZHjOLHePQAAyANjjA4dOqTy5cvL4zl+u0TMg8XOnTuVnp4e690CAAALtm/frgoVKhz3/pgHi5SUFEluYampqbHePQAAyIPMzEylp6fnvo8fT8yDxS+nP1JTUwkWAAAkmBN1Y6DzJgAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArIn5tOmnu0oPTj/hNluHdoxBJQAA2EeLBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGtOKVgMHTpUjuPonnvusVQOAABIZHkOFkuWLNFrr72mevXq2awHAAAksDwFi59++knXX3+9RowYoTPPPNN2TQAAIEHlKVj07dtXHTt2VPv27W3XAwAAEpgv0h8YN26cli9friVLlpzU9oFAQIFAIPd2ZmZmpLsEAAAJIqIWi+3bt+vuu+/W6NGjVahQoZP6mSFDhigtLS13SU9Pz1OhAAAg/3OMMeZkN546daquvPJKeb3e3HWhUEiO48jj8SgQCBxzn/THLRbp6enKyMhQamqqhYeQWCo9OP2E22wd2jEGlQAAcPIyMzOVlpZ2wvfviE6FtGvXTqtWrTpmXe/evVWzZk3179//d6FCkvx+v/x+fyS7AQAACSqiYJGSkqI6deocs+6MM85QiRIlfrceAACcfhh5EwAAWBPxVSH/a968eRbKAAAABQEtFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKzxxbsA5E+VHpz+p/dvHdoxRpXgRPhbJY4T/a1OhL+lHfzPRBctFgAAwBqCBQAAsIZgAQAArCFYAAAAayIKFsOHD1e9evWUmpqq1NRUNWvWTB9++GG0agMAAAkmomBRoUIFDR06VMuWLdPSpUvVtm1bdenSRWvWrIlWfQAAIIFEdLlpp06djrn95JNPavjw4fr8889Vu3Ztq4UBAIDEk+dxLEKhkCZMmKCsrCw1a9bsuNsFAgEFAoHc25mZmXndJQAAyOci7ry5atUqFS1aVH6/X7fffrumTJmiWrVqHXf7IUOGKC0tLXdJT08/pYIBAED+FXGwqFGjhlauXKkvvvhCd9xxh3r16qW1a9ced/sBAwYoIyMjd9m+ffspFQwAAPKviE+FJCcnq2rVqpKkRo0aacmSJXrhhRf02muv/eH2fr9ffr//1KoEAAAJ4ZTHsQiHw8f0oQAAAKeviFosBgwYoA4dOqhixYo6dOiQxowZo3nz5mnmzJnRqg8AACSQiILF3r17dcMNN2jXrl1KS0tTvXr1NHPmTF188cXRqg8AACSQiILFyJEjo1UHAAAoAJgrBAAAWEOwAAAA1hAsAACANQQLAABgDcECAABYQ7AAAADWECwAAIA1BAsAAGANwQIAAFhDsAAAANYQLAAAgDUECwAAYA3BAgAAWEOwAAAA1hAsAACANQQLAABgDcECAABYQ7AAAADWECwAAIA1BAsAAGANwQIAAFhDsAAAANYQLAAAgDUECwAAYA3BAgAAWEOwAAAA1hAsAACANQQLAABgDcECAABYQ7AAAADWECwAAIA1BAsAAGANwQIAAFhDsAAAANYQLAAAgDUECwAAYA3BAgAAWEOwAAAA1hAsAACANQQLAABgDcECAABYQ7AAAADWECwAAIA1BAsAAGANwQIAAFhDsAAAANYQLAAAgDUECwAAYA3BAgAAWEOwAAAA1hAsAACANQQLAABgDcECAABYQ7AAAADWECwAAIA1BAsAAGANwQIAAFhDsAAAANYQLAAAgDUECwAAYA3BAgAAWEOwAAAA1hAsAACANQQLAABgTUTBYsiQITr//POVkpKi0qVLq2vXrtqwYUO0agMAAAkmomDx6aefqm/fvvr88881a9Ys5eTk6JJLLlFWVla06gMAAAnEF8nGH3300TG3R40apdKlS2vZsmW66KKLrBYGAAAST0TB4n9lZGRIkooXL37cbQKBgAKBQO7tzMzMU9klAADIx/LceTMcDuuee+5RixYtVKdOneNuN2TIEKWlpeUu6enped0lAADI5/IcLPr27avVq1dr3Lhxf7rdgAEDlJGRkbts3749r7sEAAD5XJ5Ohdx5552aNm2a5s+frwoVKvzptn6/X36/P0/FAQCAxBJRsDDG6K677tKUKVM0b948nXPOOdGqCwAAJKCIgkXfvn01ZswYvffee0pJSdHu3bslSWlpaSpcuHBUCgQAAIkjoj4Ww4cPV0ZGhlq3bq1y5crlLuPHj49WfQAAIIFEfCoEAADgeJgrBAAAWEOwAAAA1hAsAACANQQLAABgDcECAABYQ7AAAADWECwAAIA1BAsAAGANwQIAAFhDsAAAANYQLAAAgDUECwAAYA3BAgAAWEOwAAAA1hAsAACANQQLAABgDcECAABYQ7AAAADWECwAAIA1BAsAAGANwQIAAFhDsAAAANYQLAAAgDUECwAAYA3BAgAAWEOwAAAA1vjiXUCiqfTg9D+9f+vQjjGq5PhiUWN+2Ed+cKqPMz88n061hpP5O+WH4xSLxxltp1pDfvjft/F8OR3+b/LDe0le0WIBAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrmDYdyIMUHVYlZ7e05TMpeFTKOfKbrwEpeETKOequK5YuVThfKl1L8njjXTpOkUdhVXd2qIHnW1V1dipHXgWUrKMmWQEl6aiSddQkueuUpCPya7sppR2mlAyf5XAaIFgAf6K4MlXN+V7VPDtUxdmpas4OVfXsVFnnoLvBWxH8sqQzpLMauiGjwvlShcZRqRmWZe2XdizRfb7xauBsVH3PJhV1jkb8a46aJG0y5bXRnKVvw2dpozlLG015bTNllcNLMQoQns3Az0rroBp5vlFDz7eq59msas4OFXd+Ou72e00xlS5VRkoqJPkK//z15yWpsPvVmyTt/0basUzKPiRt/cxdfjY/uZRWmGqaHWqo98PNJDkxeKQ4oY2fSF+Nl3YskQ5ukSTd+ZtXy59MIa0MV9E6c7aMHBVStgopW34n53ffn6GjqujsVSEnR7WdbaqtbdJvGq5yjFfbTBmtMZW0LFxNy8PVtN5UVJCXZyQonrk4LfkUVE3nOzXyfJsbJio4+3+3Xdg42mFK6ltTQRtzP21W0CZTXodURFvv7HhyOwyHpH0b3DeqHUukHUulfetV0bNPFbVPXbyLdG1orgYEb9F3pozlR4uTlrVf+rC/tHrisetL1dT4XWW1wlTVinBVfWsqKBzBaQ2Pwkp39qqa872qOt+rqmenqjo7VNXZqaLOUVV1dqqqdqqLd5Ek6bDx62tTOTdoLA9X00Gl2nykQNQQLHBaSFaOGjgb1dy7Wk0961XP2awiTuCYbULG0XpTUcvD1bQiXFUbTEVtMuV0VP5TL8DjlcrUcpdGvdx1RzN0/eDhauFZo5u8H6qFd41mevrr2eBf9EaoQ0RvXDhVRvr6XTdUHDkgOR7p/Fuk6pdJZzWSChdT/wen5/m3h+XRNlNW20xZfaJGUujX/ZbTAVX37FB9Z5Maer5VQ8+3SnUO6wJnnS7wrMv9HZvDZbU0XEMLw7W1KFxH+1TslB4xEC0ECxRIjsKq5WxTC89qtfCsURPPehV2so/ZJsMU0fJwNS0LV9dyU01fhasoS4VjV2ShNC0M19XCcF2ND7XWUN/rauZdq4eTRusK72L1z+kTu1pOY+X0g55MGilNXumuKF1b6vJvtz9M1DnapRLaFS6hT1VfCrnP3SrOTjXyfKuGjtuiVtWzU5U9u1XZs1vd9Kkk6ZvwWVoYrqOF4Tr6InyuDqlIDOoFToxggQKjorNHLT2r1NyzWs08a3/XP2KfSdWicB0tDtfSknANbTbl8k0v/W2mrHrkDFT38Fw95But8zybNS15oDTngHTRfZLPQqsJjuEorOu9s9XfN04pzhHJmyxd9IDU4m7Jlxy3uow82mgqaGOogsarjSQpTT+poedbXeBZqxaeNarlbFN1z/eq7vlevTVTIePoa1NFC8O1tTBcR0vDNegQirjhmYeElaSgzvesV1vPCrXxrFQVz65j7j9kCuuLcE0tCtfRgnAdfWMqKH93jnQ0LtRWc0IN9HjSm7rUu1Sa/7S09j2p80tSxabxLrDAqOzs1JCk19XUs16StDRcXY37viOVqhHnyv5YhopqbriB5oYbSJKK6ZCaedaqhWe1mnvWqLJntxo4G9XAs1F36j0dMoW1IFxHc8INNC9UX/t0ZpwfAU4nBAsklFI6qDbelWrrWakLPauOuewvx3i13FTTgpDbPPy1qZyQPev36kzdlnOvOoS+1PAzx0j7N0hvXCo1vV265AnJm3iPKT+52TtDD/jGy+/kKMv49a9gd70dulhb8mmo+CM/KkUfhpvqw7AbNstrv5p716i5Z41aer5WKSdTHbxL1MG7REqSVoUraU64gbS9tHuKh/FUEEW8QiGfM6rtbNMl3qVq61muup6tx9y7z6RpXqi+5oQbaEG4bgE6z+y4bxp975E+flhaOVr6YrgkI3X4V7yLS1g9vLP1z6R3JEmfhurpoZyb9b1KxbmqU7dTJTUx1EoTQ63kKKy6zha19a5Qa89KnefZrLqere7/zsgpUpESUtWLpZqXS1XbS8lnxLt8FDAEC+Q/4ZC0/Us97Htbl3qWKt2z75i7V4Yra26ogeaEG2i1qZRv+klERZHiUtdXpCptpUk3S1+86jbXN74p3pUlnGaeNRrsGyVJeiF4lZ4PXq38fWosb4w8+tpU0dfBKhqma1RSGWrl+UptvCt0RZF10uEfpK/HuYuvkFSlnXRuJ6n6pe7zDThFBAvkD8Fsaet8ad0H0voZUtZe3fLzs/OISdan4fqaFWqkT8P1tV9p8a01Hupe4w7UNOcJacb9Uomq0jkXxbuqxPHDJg1PGqYkJ6T3Qs0LbKj4I/uVpknhizQpfJGueOASafsX0oYP3f+1H7dJG6a7i+OVzmnphowaHaXUcvEuHQmKYIH4CWZLm2ZLa6ZIGz6SAhm/3lcoTZOy6urj0Pn6NFzPzlgSia7lfe4gW6smSON7SrfOkUpUiXdV+d+RH6Wx3VXMydLKcBU9kNNHp0uo+B1vklTpQne55Alpz2o3YKybJu1dI22e5y7T/+EOO1/7SndJLR/vypFACBaIrXBI2jJfWj1JWve+dPQ3YaJoGalmR6nmFVKllvrHw7PiV2d+5Dju1SEHNkvfL5PGXCvd8km8q8rfQkFp4k3S/m+00xTXrdn9FFD8LiXNVxxHKlvXXdo8JP2w6ecWw2m/GSF2iTRzoHR2C6nOVVKtLvGuGgmAYIHoM0ba/qUbJtZMkbL2/npf0bK/fiqqcL7kKcD9JWxIKix1HyONaCv98K00sbe86q2Q6OX/hz5+2G0V8xXWrVn/4LLLP1OiinThPe6SudMNGasnS9s/l7YtcJcZ92tUUh19EGqmj8ONC1BnadhEsED07FkjfT1eWj1Fyvju1/WFz3Q/+dS5Rjq7OZe+RSqlrNRjrPTGZdKmOXrYl6THgr3iXVX+s/TNn6+kkXTVa1rzX17uTlpqeanpbe7y43ZpzWT3g8Gur9Ta6y4Bk6S54fM0NdRCc8INlK2keFeNfIL/NFhVXJnq7F0kvTpU2v31r3ckF3VPcdS5WqrSxj3Xi7wrV1+68jXp3Z7q7ZupjeYsjQ61j3dV+ceW+dKM+9zv2zz8cxN+3uf6OK0VS3dHI21xt7T/Wz33/FB19i5SVc9OXeZdosu8S/SjOUPvh5prUqil20LpnKZ9WCCJYAELkhRUG88KXeOdrzaelUpyQtJuSZ4k9xK2un9xvybFcB6O00GtzlLbh6U5T+gx3yhtNuW0OFw73lXF3w+bpHdvkMJBt1XsovviXVHBUbKaXgxdpRdDV6qms11dvAvV1btQ5ZwDusE3Szf4ZkkvvyPV7yHV706nz9MUwQJ5ZFTH2aKrvZ+pi3fhMfNyfBWurPod73BbJ84oEccaTwMt79PUWXPU1btIw5OGqWv2YG01p/FlgkczpLHdpSMH3VlJu/ybT89R4c4EvD5YUf8XvFbNPWt0jfdTXepZqsL7v5FmPybNHixVbi2dd53bWplMf4zTBcECEUlVlrp4F+o67xyd6/m138QeU0xTQhdqUugifWsqaGvTjnGs8jTiOOqf00dnO3vVwLNRI5Oe0ZXZg5Wp03A0xXBImtBb2v+NlHqW28mVVrKoC8ujBeG6WhCuq6I6rNXdjkgrx0rfLZI2z3UXf6pU71qpEX2BTgcEC5wEo4bOt+rhnaMrvJ/nTj8eMEn6ONxIk0IX6bNwXa5MiJOAktUnu5/e8z+sKp5detA3Vg8Fb4l3WbG39I3cK0DUfYzbyRUx9ZOKSA3/IjW8QTqwRfpqnPTVWHcgriUjpCUjNCW5qsaE2mpa6AIdUaF4l4woIFjguFKVpa7eBerhnaNzPdtz128IV9CYUDtNCbVQporGsUL8Yp+K6e7sOzXBP1jdvXM1PtRaX5mq8S4rdn7aJ81+3P3+ksel8ufFtRxIKn6O1GaA1Kq/tOVTadkoaf00NZA7C+sjvrc1NdRCY0NttdZUine1sIhggWMZI+1YomeSXlVHz6+tE0dMsqaFLtDYUFstN9V02o5cmI8tMTU1KdRSV3s/0+CkUboye7DCBXkeld/6ZJA7cmu5+syjkt94PO6VYFXaSD/t1ZAhj6i7d47O8exRT98n6un7RF+FK2tsqK3eCzWnFaMAiPhVZ/78+erUqZPKly8vx3E0derUKJSFmMs5Kq0YLf2nlTTyYl3jna/CTrbWhdP1SE4vNQ28rPuDt2u5qS5CRf41JOc6ZZrCqu/ZrO7eufEuJyYaORvc2V8lqeNzjIuSnxUtrddCndQ2+1n1yB6oD0IXKNt4Vd+zWUOTXtfn/jv1kG+0exoFCSviFousrCzVr19fN910k6666qpo1IRYytghLRkpLX/LnfVQknyFNOFoE40JtdMKU1UEicSxX2l6NthNjyW9pQd84/RR6Px4lxRVXoX0eNIo90bDG6QKjeNaD06OkUeLw7W1OFxbxZWpq73z9VfvJzrbs1d9fNOlF2e4l6g36SNVbsOIvAkm4mDRoUMHdejQIRq1IFaMkbYtlL54TVo/XTIhd31aunT+zVLDXrp/8OL41og8eyfUXtd656mWZ5se8I2T1CPeJUXNX72fqJZnm1SomNTu0XiXgzw4oFSNCF2h10OXq7XnK93onalW3q+lbz5ylxLVpCa3umNjFEqNd7k4CVHvYxEIBBQIBHJvZ2ZmRnuXOA6/stXVu1Aa/qQ7k+EvKrV0h+6t3kHy0u0m0YXk1T9zbtQk/2Pq7psnbV8ipRe8lotS+lH/8L3r3mg/iDFTEpyRR3PDDTQ33EBb764ufTlCWjnGnRPnwwfczrnn9ZCa3h7vUnECUW9fGjJkiNLS0nKX9PT0aO8S/6O4MnW3d5IW+v+ufyWNcENFUhGpUW/pjsXSjdOkczsRKgqQZaaGJgQvcm9M7+eO8VDAPJg0RqnOEa0MV5YaMj5CgVKymnT509I/1kkd/s9ttcg+JH35H+mlRno16Xk1dL6Jd5U4jqgHiwEDBigjIyN32b59+4l/CFac4+zSk76RWuS/S/cmTVJJJ1M7TEnpkielfmulTsOkMrXiXSaiZGiwhzJNEXfOlqVvxLscq5o463S1d4HCxtEjOb3psFlQ+VOkpn2kO5dIPadI1S6RZHSZd4km+x/V5ORHdJnnS3kUjnel+I2of0T1+/3y+/3R3g1+YYwaO+vVxzdd7T3L5XGMJHeY7RHBjvow3ESbmneOc5GIhR+Upv8LdnM7N855XKrVVSpaKt5lnTKfghr8c4fNsaG2+tpUiW9BiD7Hkaq0dZe96zXuxf660rtADT0b9WryMG0Ll9bIUAdNCLXictV8gK62BUU4JK2ZIr3eThP9g3WJd5k8jtGsUEN1C/xTXbIf17RwM0bHPM2MDrWXytZz59D4ZFC8y7HiBu8s1fRs1wFTVP8X7BbvchBrpWvqwWAfXRh4US8Gu+qgKaqzPXs1OOktLfbfpft846VDe+Jd5Wkt4haLn376SRs3bsy9vWXLFq1cuVLFixdXxYoVrRaHkxDMlr4eJy14XjqwWZI71PakUEuNDHXQJnNWnAtEPIXlkTo+K4282B3roeENUsUL4l1WnpXSQd3rmyhJ+lewh35USpwrQrzsUzE9F+ym4cHOutr7mW7xzlAlzx7d6XtPemGm1KCn1OLvUjHel2It4mCxdOlStWnTJvd2v379JEm9evXSqFGjrBWGE8g+LC3/r7ToRSnze3dd4TOlJn3UfGYl/aC0+NaH/CO9ifsiu+Jtafp9Up95CdtR96GkMUpxjmhluIreDbWKdznIB46okN4JXawxoXa62LNMt/s+UIPgRndukmVvupOfXXiv2yEUMRHxq0vr1q1ljIlGLTgZRzOkJa9Li1+RDu931xUtIzW/y73Kw19UP8ycHt8akf+0f1Ra94G0Z5W0dKR7eXGCaeqs05XehQobRw/n9JbhTC5+IyyPZobP18zsxtraJ0X67Fl3jpKVo93LVmt1kVr2c4d9R1Ql5seW01HWD9Lnr7jXdgcy3HXFKkot7pHOu15KosMS/sQZJaV2j7iXns55Qqp9pVS0dLyrOmluh803JUmjQ+202lSOc0XIvxypcit32bFU+uw5acN0ae1Ud6l6sRo5zbXM1Ih3oQUWwSK/+2mftHCYe7lgzmF3XckabvKuc7XkTYpreUggjW50T5/tWinNekS68tV4V3TSenlnqoZnh34wKXqGDps4WRUaSz3GSHvWuP3QVk+SNs7SJP8sLQ7V0rDg1frCnBvvKgsc2hLzq6wf3Bf/F+pJi//thopy9aVub0t/+1yq351Qgch4vG5HTkn6aqz03RfxredkHdqje3yTJUn/CnZXhorGuSAknDK1patfl+5aJjW6UQHjUzPvWo33P67RSU+qsbM+3hUWKASLfKaYDkmfPCYNqystfMENFOUbSNdNkPp8KtXqzIQ8yLsKjaXz/up+/+H9iTEi5yeP/txhs7Im0GETp6J4ZanTC2odeF5vB9sr23jVwrtGE/2D9XbSU4zmaQmnQvKJNP2kW3wzdKN3prTgiLuyXH2p9UPuLH8OM4zCkvaDpHXvS7u+ck+NKB/3tdj+pfTVGEnSozk30mETVuxSCf0zeJOGBzurr+89/cU7Ty29q9XSu1rzQ3Wl7aUL5Pw6scJ/aZylKkv3+ibqM//duss3VSnOEalsXan7GLeFosZlhArYVbS01HqA+/3swUrTT/Gt53jCIWnG/ZKkd4OttNJUjXNBKGh2qqQGBm9Wm8BzGhNsoxzj1UXeVdLI9tI710jfL4t3iQmJYBEnhRTQbd4PNN9/j+72TVaqc0TrwhV1W/a9Up/5Us2OBApET5NbpVI1pSMHcgecyndWvON2NPWn6ulg93hXgwLse5XSQ8Fb1Sb7WY0PtpYcr7RxljSirTT+r9K+DfEuMaEQLGLMq5B6eGdrnr+fBiSNVTEnSxvCFXRH9t26PPspzQyfTx8KRJ83SerwtCSpp3eWajrfxbmg/3HkoDT7Mff71gO0nwHfEAM7TGn1D/aR7loq1b9Ocjzu+C+vXCC911f6kUk0TwbvYLESDktrpujj5Ac0JGmkyjoHtcOUVL/s29Uhe6g+DDfl/DFiq3IrqVYXeR2jx5JGScpHA9/NHSId/sFtVWlya7yrwemmeGXpyuHSHYulmldIJuy2oL3USJo50L1qD8fFO1m0GSNtmiONaCNNuFFVPLu036Tq0Zwb1DbwrCaHL3LncwDi4ZIndcQkq6lnvTp5Fse7GteeNe7ospLU4V9cVo34KV1T6j5auvkTqVJLKRRwL/9/ob7u8k5WER2Nd4X5Eu9o0bRjmfTfztLbV7rnipOL6rmca9Qq8LxGhS5TtnjBRJwVS9crwc6S3Hk44v5CaYw04wHJhKRzO0uVW8e3HkByrxDp9YH018nubMHZh/SPpIma779HvbwzlaRgvCvMVwgW0XBwmzTxJun1ttKW+ZI3Wbrgb9LdX+nF0FXKUuF4Vwjk+k/oCn0XLqVyzgH19U2NbzFrJkvbFki+wtKlT8a3FuC3HEeq2s69Wu+aN7U5XFYlnUw9lvSWPk6+X5d6vlS+Op0YRwQLi1J02B0t89/nu0PHypHq93BHe7tsiDtfA5DPBJSsx4M9JUm3eGfobGd3fArJzpI+/qf7/YX3Mt018iePR6pzlS7JfloP5dysfSZN53j26LXkYRqf/LjqOZviXWHcESws8Cmont6PNc9/rztaZiggnXORdNun7nwMvEAin5sVbqRPQ/Xkd4J6xPd2fIr47Dkp8/ufJ9f7e3xqAE5SUD6NCbVT68BzeiF4ZW5fpff9/9SwpH9LP+azK61iiGBxSozaeZZpZnJ/PZ40SiWcQ1LJ6lKP8dIN7zM9LxKIo8eCNyjbeNXOu0JtPCtiu/sDm6VFL7rfX/qUlMTpQiSGLBXW88G/qE3gWU0MXaSwcdTVu0h6qbE0a5B0NCPeJcYcwSKPajtbNSbpSY1MflZVPLv0g0nRwzm9pTsWMVomEtJmU15vhDpIkh7x/VfJyondzj96SAplS5XbuJf3AQlmt0rovpzb1Sn7SS0K1XJbrhcOk15sIH05QgqdPh08CRYRKqEMDfGN0AfJA9Xcu1YBk6RXgp3VOvC83gldzKVxSGgvBa/UHlNM53j26Fbv9Njs9JuZ0jcfSh6fO2gXoRwJbI2ppOtyBrot1yWru+OxzLhPevVCafO8eJcXEwSLkxXKkRa/orn+f6iHb648jtF7oeZqG3hGTwe765CKxLtC4JRlqbCG5FwnSernm6BLPEuiu8Pdq6VJPw+A1fR2qVT16O4PiAnHbbm+Y5F0+TNS4eLSvnXSf7u4Q4Qf3BbvAqOKYHEyNs2VhreQZg5QqnNYq8KVdHVgkO7OuVPfq1S8qwOsmhpuoXeDreR1jF5K+re0dWF0dnRwq/TO1VIgQ0q/QGozMDr7AeLFm+SOHPv35VKT29w5SNZ9IL3cRPf6JqqQAvGuMCoIFn/m4FZp3PXS212l/RukIiX0YM4t6pL9hJaZGvGuDogSRwOCt2hWqJH8To40tru0e5XdXfy0zx047qfdUula0nXjpGRa/VBAFT5Tuvxp6fYF7hWDwaO62zdZs/33qaPncxW08S8IFn8kO0ua84T07ybS+mluymx6h3TXMo0LtWUIbhR4IXl1Z85d+iJcUwpkSm9f5V65YcPRTGn01e7vS6vojmZY+Ew7vxvIz8rUcq8Y7PZf7TAldZbzg15OflHjkp/QuU7BOT3CO+RvGSOtmeoOcDX//34dj+KOhVKHobz44bQSULJuyb5PKlNXytrrtjAc2nNqvzTnqDTuOmnXV1KRklLPKVJqOTsFA4nAcaRaXdQu8Iyey7lGR0yyLvCs07Tkh/Sob5RSlRXvCk8ZweIXP2ySRl8jTejlDtKTVlHq9rabLkufG+/qgLg4pCLSXydJZ1b6tU9EXq/LD4ekybdKWz+TkotKf50olaxqs1wgYQSUrBdDV6ld4BlNC10gr2N0o+9jzfbfpy6eBe4H3QRFsMg5Ks37l/RKM2njJ+68Hq36S3d+KdXqzKVvQEoZt2XhjNLSnlXS2B7u/00kjJGm/0Na9777P9Z9jFS+QXTqBRLITpXUnTl/V4/sgdoULqdSToZeSH5FequTtO+beJeXJ6d3sNg0RxreXJr3lHvao3Jr6Y7FUpuHGPkP+K3ild0WBn+qtG2hNOnmyAb8mfuktOxNSY501QipcquolQokosXh2uqQPVRP53TTUZPktuwNby7NHixlH453eRE5PYNF5i5pwo3uOeMDm6SiZaVr3pB6TqVpFjiecvWlHmMlr9/t1DztnpNqrr3R+5HbZ0mSrnhOqt01qmUCiSpbSXol1FXts/9Pqn6ZFM6RPntWermptOHDeJd30nzxLiCmQkHpy/9Ic5+Ssg9Jjse9trjNQ1Kh1HhXB+R/lS50Q/i7PaUVb0tZ+zXYd/yWi8IK6C+++e6NNg9LjW+KUaFA4tphSks9xkkbZkgf9pcyvnMv+65xudThX/l+YsvTJ1jsXCl98He3N7okVThf6vicVK5eXMsCEs65V0idXpDev0v65kPdcDKvIk1uky66L+qlAQWG40g1O7qn6D99Wlr8bzdobP5Uavuw1PQ2yeONd5V/qOAHi+zDbh+Kxa9IJiQVKia1f1Rq2EvynJ5ngoBT1vAGKfUsafuXGjb7zzuYbQ2X1bDLnqQjNJAXyWdIFz8m1e8hTbtX+m6RNHOAtGqC1PlFqWzdeFf4OwU7WGya654HPrjVvV37KrcZqWjpeFYFFAxV20lV22nYzBNPVjaMEA+cmtI1pRunSyv+K338iLRzufRaK6nF390rGfPRBQcF87/98AFpyh3uUNwHt7qfrHqMl/7yJqECAJCYPB6p0Y3ucAjndnZb4Rc87149svnTeFeXq2AFC2Okrye4I2d+NUaS457b7fuFO9McAACJLqWsdO3b7ngwKeXc4fH/21l6r6/7wTrOCk6w+PE7d+TMybdIh/dLpc6Vbp7lTvziT4l3dQAA2FWzo/vBufHN7u0V70gvN5FWT4rryJ0Fo49F5i7p5QuknCx3VL+LHpBa3C35kuNdGQAA0VMozR0fpl436f2/uzNxT7xJKlLCvaIkDgpGsEgt5w66c2CLexlcqerxrggAgNipeIF0+2dun4tdX0nnxG9024IRLCTp8mckXyEuIQUAnJ58fqn1g+5pkDhe3l1wgkVykXhXAABA/MV5zBg+3gMAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArCFYAAAAawgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsCZPweLll19WpUqVVKhQITVt2lRffvml7boAAEACijhYjB8/Xv369dOgQYO0fPly1a9fX5deeqn27t0bjfoAAEACiThYPPfcc7r11lvVu3dv1apVS6+++qqKFCmiN954Ixr1AQCABOKLZOPs7GwtW7ZMAwYMyF3n8XjUvn17LV68+A9/JhAIKBAI5N7OyMiQJGVmZual3rgLBw7/6f0nelwn+vmT+R0ncqo1nszvOBEbf99TrSEWEuFvdarPydPlOW3jceZ3sfi/jMXzJd7P+fxSQ6z9UpMx5s83NBH4/vvvjSSzaNGiY9bff//9pkmTJn/4M4MGDTKSWFhYWFhYWArAsn379j/NChG1WOTFgAED1K9fv9zb4XBYBw4cUIkSJeQ4jrX9ZGZmKj09Xdu3b1dqaqq134tjcZxjh2MdGxzn2OA4x0Y0j7MxRocOHVL58uX/dLuIgkXJkiXl9Xq1Z8+eY9bv2bNHZcuW/cOf8fv98vv9x6wrVqxYJLuNSGpqKk/aGOA4xw7HOjY4zrHBcY6NaB3ntLS0E24TUefN5ORkNWrUSLNnz85dFw6HNXv2bDVr1izyCgEAQIES8amQfv36qVevXmrcuLGaNGmiYcOGKSsrS717945GfQAAIIFEHCyuvfZa7du3T4888oh2796t8847Tx999JHKlCkTjfpOmt/v16BBg3532gV2cZxjh2MdGxzn2OA4x0Z+OM6OOeF1IwAAACeHuUIAAIA1BAsAAGANwQIAAFhDsAAAANYkVLCIdLr2CRMmqGbNmipUqJDq1q2rGTNmxKjSxBbJcR4xYoRatmypM888U2eeeabat29/wr8LXJE+n38xbtw4OY6jrl27RrfAAiTSY/3jjz+qb9++KleunPx+v6pXr87rx0mI9DgPGzZMNWrUUOHChZWenq57771XR48ejVG1iWn+/Pnq1KmTypcvL8dxNHXq1BP+zLx589SwYUP5/X5VrVpVo0aNim6RkcwVEk/jxo0zycnJ5o033jBr1qwxt956qylWrJjZs2fPH26/cOFC4/V6zdNPP23Wrl1rHn74YZOUlGRWrVoV48oTS6TH+brrrjMvv/yyWbFihVm3bp258cYbTVpamtmxY0eMK08skR7nX2zZssWcddZZpmXLlqZLly6xKTbBRXqsA4GAady4sbn88svNggULzJYtW8y8efPMypUrY1x5Yon0OI8ePdr4/X4zevRos2XLFjNz5kxTrlw5c++998a48sQyY8YMM3DgQDN58mQjyUyZMuVPt9+8ebMpUqSI6devn1m7dq156aWXjNfrNR999FHUakyYYNGkSRPTt2/f3NuhUMiUL1/eDBky5A+379atm+nYseMx65o2bWpuu+22qNaZ6CI9zv8rGAyalJQU89Zbb0WrxAIhL8c5GAya5s2bm9dff9306tWLYHGSIj3Ww4cPN5UrVzbZ2dmxKrFAiPQ49+3b17Rt2/aYdf369TMtWrSIap0FyckEiwceeMDUrl37mHXXXnutufTSS6NWV0KcCvlluvb27dvnrjvRdO2LFy8+ZntJuvTSS4+7PfJ2nP/X4cOHlZOTo+LFi0erzISX1+M8ePBglS5dWjfffHMsyiwQ8nKs33//fTVr1kx9+/ZVmTJlVKdOHT311FMKhUKxKjvh5OU4N2/eXMuWLcs9XbJ582bNmDFDl19+eUxqPl3E470w6rOb2rB//36FQqHfje5ZpkwZrV+//g9/Zvfu3X+4/e7du6NWZ6LLy3H+X/3791f58uV/90TGr/JynBcsWKCRI0dq5cqVMaiw4MjLsd68ebPmzJmj66+/XjNmzNDGjRv1t7/9TTk5ORo0aFAsyk44eTnO1113nfbv368LL7xQxhgFg0Hdfvvteuihh2JR8mnjeO+FmZmZOnLkiAoXLmx9nwnRYoHEMHToUI0bN05TpkxRoUKF4l1OgXHo0CH17NlTI0aMUMmSJeNdToEXDodVunRp/ec//1GjRo107bXXauDAgXr11VfjXVqBMm/ePD311FN65ZVXtHz5ck2ePFnTp0/X448/Hu/ScIoSosUiL9O1ly1bNqLtkbfj/ItnnnlGQ4cO1SeffKJ69epFs8yEF+lx3rRpk7Zu3apOnTrlrguHw5Ikn8+nDRs2qEqVKtEtOkHl5Tldrlw5JSUlyev15q4799xztXv3bmVnZys5OTmqNSeivBznf/7zn+rZs6duueUWSVLdunWVlZWlPn36aODAgfJ4+Nxrw/HeC1NTU6PSWiElSItFXqZrb9as2THbS9KsWbOY3v1P5OU4S9LTTz+txx9/XB999JEaN24ci1ITWqTHuWbNmlq1apVWrlyZu3Tu3Flt2rTRypUrlZ6eHsvyE0pentMtWrTQxo0bc8ObJH3zzTcqV64coeI48nKcDx8+/Lvw8EuYM0xhZU1c3guj1i3UsnHjxhm/329GjRpl1q5da/r06WOKFStmdu/ebYwxpmfPnubBBx/M3X7hwoXG5/OZZ555xqxbt84MGjSIy01PQqTHeejQoSY5OdlMnDjR7Nq1K3c5dOhQvB5CQoj0OP8vrgo5eZEe6++++86kpKSYO++802zYsMFMmzbNlC5d2jzxxBPxeggJIdLjPGjQIJOSkmLGjh1rNm/ebD7++GNTpUoV061bt3g9hIRw6NAhs2LFCrNixQojyTz33HNmxYoVZtu2bcYYYx588EHTs2fP3O1/udz0/vvvN+vWrTMvv/wyl5v+1ksvvWQqVqxokpOTTZMmTcznn3+ee1+rVq1Mr169jtn+3XffNdWrVzfJycmmdu3aZvr06TGuODFFcpzPPvtsI+l3y6BBg2JfeIKJ9Pn8WwSLyER6rBctWmSaNm1q/H6/qVy5snnyySdNMBiMcdWJJ5LjnJOTYx599FFTpUoVU6hQIZOenm7+9re/mYMHD8a+8AQyd+7cP3zN/eXY9urVy7Rq1ep3P3PeeeeZ5ORkU7lyZfPmm29GtUamTQcAANYkRB8LAACQGAgWAADAGoIFAACwhmABAACsIVgAAABrCBYAAMAaggUAALCGYAEAAKwhWAAAAGsIFgAAwBqCBQAAsIZgAQAArPl/mVfo6Nk3fY4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rng = np.random.default_rng(2)\n",
    "x = rng.uniform(0, 1, size=35)\n",
    "\n",
    "cost = ExtendedUnbinnedNLL(x, model)\n",
    "n = len(x)\n",
    "m = Minuit(cost, b0=n, b1=n, b2=n, sig=0, mu=0.5, sigma=0.05)\n",
    "m.print_level = 0\n",
    "m.limits[\"b0\", \"b1\", \"b2\"] = (0, None)\n",
    "m.fixed[\"mu\", \"sigma\"] = True\n",
    "display(m.migrad())\n",
    "\n",
    "plt.hist(x, bins=50, density=True)\n",
    "xm = np.linspace(0, 1)\n",
    "yint, ym = model(xm, *m.values)\n",
    "plt.plot(xm, ym / yint);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1fac22f",
   "metadata": {},
   "source": [
    "In this example, the signal amplitude came out negative. This happens if the background has an underfluctuation where the signal is expected. This is not an issue if the sum of signal and background density is still positive everywhere where we evaluate it. As long as the total density is positive, individual components are allowed to be negative.\n",
    "\n",
    "There are, however, no principle restrictions in this example that prevent the sum of signal and background from becoming negative for some toy data sets. When that happens, the fit will fail, since the total density cannot mathematically become negative.\n",
    "\n",
    "If this happens anyway, the fit will fail since taking logarithm of a negative number will cause havoc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ecf9073",
   "metadata": {},
   "source": [
    "## Migrad fit on toys\n",
    "\n",
    "We apply the fit many times on randomly sampled background-only data to observe this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "following-bruce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28 failed\n"
     ]
    }
   ],
   "source": [
    "@joblib.delayed\n",
    "def compute(itry):\n",
    "    rng = np.random.default_rng(itry)\n",
    "    x = rng.uniform(0, 1, size=35)\n",
    "    cost = ExtendedUnbinnedNLL(x, model)\n",
    "    m = Minuit(cost, b0=n, b1=n, b2=n, sig=0, mu=0.5, sigma=0.05)\n",
    "    m.limits[\"b0\", \"b1\", \"b2\"] = (0, None)\n",
    "    m.fixed[\"mu\", \"sigma\"] = True\n",
    "    m.migrad()\n",
    "    return m.values[\"sig\"] if m.valid else np.nan\n",
    "\n",
    "sigs_migrad = joblib.Parallel(-1)(compute(i) for i in range(200))\n",
    "\n",
    "print(np.sum(np.isnan(sigs_migrad)), \"failed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bd848895",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAHHCAYAAAAf2DoOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAw8klEQVR4nO3deXgUVb7G8TchJGFJd0jIQoAE3MIOghKCO2SITGRgiAxwGWVzncCwXUe4dzSidwgqyijDNg4GLw6joqIiAiIKLgSEAKMgoDiBwECCgukAkgVy7h8+9LXNQjqEkzR+P8/Tz0OdOnXqVymafqk6XfEzxhgBAABY4l/XBQAAgJ8XwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHUI+dPHlSd911l6Kjo+Xn56eJEydq//798vPz0+LFi63v2xt+fn565JFH3MuLFy+Wn5+f9u/fX2s1jho1Sm3atKlW37KyMnXq1El/+tOfam3/tixYsECxsbEqLi6u61KAWkH4ALy0ZcsWjRs3Th07dlSTJk0UGxur3/zmN/ryyy8r7P/KK6+oV69eCg0NVXh4uG666SatXLmyWvuaMWOGFi9erPvvv19LlizRHXfcUWG/d955x+ODvjZUd9++4h//+IcOHjyocePGudtOnjyp9PR03XrrrQoLC6t2qCstLVWHDh3k5+enWbNmnbe/MUbTp09Xy5YtFRkZqYkTJ6qkpMSjz8mTJ9WyZUstXbq03PajRo1SSUmJFi5ceP4DBXyBAeCV1NRUEx0dbcaPH2+ee+4589hjj5moqCjTpEkT8/nnn3v0ffbZZ40kk5KSYubPn29mz55tunbtaiSZ11577bz7SkhIMNddd51HW1lZmTl9+rQ5c+aMuy0tLc3U9tu5on174/Tp06a0tNS9nJmZaSSZnJycWqjuByNHjjRxcXHV6tu1a1dzzz33eLTl5OQYSSY2NtbcfPPNRpLJzMw871hPPfWUadKkiZFknnzyyfP2X7JkiQkMDDQPPfSQmTlzpgkJCTEzZszw6DN16lTTu3fvSsf4wx/+YOLi4kxZWdl59wfUd4QPwEuffPKJKS4u9mj78ssvTVBQkBkxYoRH+5VXXmmuvfZajw8Ml8tlmjZtan71q1+dd19t27Y1KSkp5+13McJHdfddXXUZPrZt22Ykmffee8+jvaioyBw5csQYY8yWLVuqFT7y8/ON0+k0jz76aLXDx9ChQ83o0aPdy+np6aZXr17u5X379plGjRqZLVu2VDrG1q1bjSSzbt268+4PqO+47QJ4qXfv3goMDPRou/LKK9WxY0ft3r3bo72wsFCRkZHy8/NztzkcDjVt2lSNGjWqdB/r16+Xn5+fcnJytHLlSvn5+bnnS/x0zseoUaM0d+5cSXL3+/H+XnrpJfXo0UMhISFyOBzq3LmznnnmmRrtu6SkRA8//LB69Oghp9OpJk2a6IYbbtAHH3xQbpyfzvmozKpVq3TDDTeoSZMmCgkJUUpKinbt2lWu3xtvvKFOnTopODhYnTp10vLly8879o+3DQwM1I033ujRHhQUpOjo6GqPI0lTp05VfHy8fvvb31Z7m9OnT6tZs2bu5bCwMH3//ffu5SlTpmjYsGG65pprKh2jR48eCgsL05tvvulVvUB9FFDXBQCXAmOM8vPz1bFjR4/2m2++Wa+++qrmzJmjAQMGqKioSHPmzJHL5dKECRMqHa99+/ZasmSJJk2apFatWmnKlCmSpIiICH3zzTcefe+9914dPnxYa9eu1ZIlSzzWrV27VsOHD1ffvn31+OOPS5J2796tTz75pNL9V7XvwsJC/e1vf9Pw4cN1991368SJE1q0aJGSk5P16aefqlu3bl793JYsWaKRI0cqOTlZjz/+uL7//nvNnz9f119/vbZv3+6eTPruu+8qNTVVHTp0UEZGho4dO6bRo0erVatW1drPxo0b1alTJzVs2NCr+n7q008/1QsvvKCPP/7YI+Cdz7XXXqt58+ZpyJAhatKkiRYuXKjevXtL+uEcvf/++5XOGfqx7t2765NPPqlx/UC9UdeXXoBLwZIlS4wks2jRIo/2/Px807dvXyPJ/WrevLnZuHFjtcaNi4srd+vj3DyFH98eqOy2y4QJE4zD4fCYH1JdFe37zJkz5W45fffddyYqKsqMGTPGo12SSU9Pdy//9LbLiRMnTGhoqLn77rs9tsvLyzNOp9OjvVu3bqZFixamoKDA3fbuu+8aSdW67dKqVSuTmppaZZ/z3XYpKyszPXv2NMOHDzfG/P95qM5tl8LCQnP99de7/w507NjRHDp0yJSWlpoOHTqYmTNnnncMY4y55557TKNGjarVF6jPuO0CXKA9e/YoLS1NiYmJGjlypMe6xo0bKz4+XiNHjtSyZcv0/PPPq0WLFho8eLD27dt30WsLDQ3VqVOntHbt2loZr0GDBu5bTmVlZTp+/LjOnDmja665Rtu2bfNqrLVr16qgoEDDhw/Xt99+6341aNBACQkJ7ls5R44c0Y4dOzRy5Eg5nU739r/4xS/UoUOHau3r2LFjHrc9amLx4sX6/PPP3VeQvBESEqINGzZo165d2rFjh3bs2KGWLVtq3rx5Ki4u1qRJk/TFF1/olltuUcuWLfXb3/5WhYWF5cZp1qyZTp8+7XHLBvBFhA/gAuTl5SklJUVOp1OvvvqqGjRo4LF+yJAhys3N1eLFi3X77bdr9OjRWr9+vUpKSvTf//3fF72+3/3ud7rqqqvUv39/tWrVSmPGjNHq1asvaMwXXnhBXbp0UXBwsMLDwxUREaGVK1fK5XJ5Nc5XX30lSerTp48iIiI8Xu+++66OHj0qSTpw4ICkH+bV/FR8fHy192eM8aq+HyssLNS0adP0wAMPqHXr1jUaw9/fXx06dFDXrl0VEBCgb7/9Vo888ohmzZolPz8/3XbbbercubPefPNN5ebmavz48ZUegze3fID6iDkfQA25XC71799fBQUF+uijjxQTE+Ox/l//+pdWr16tv/71rx7tYWFhuv76663cu4+MjNSOHTu0Zs0arVq1SqtWrVJmZqbuvPNOvfDCC16P9+KLL2rUqFEaNGiQHnjgAUVGRqpBgwbKyMjQ119/7dVYZWVlkn6Y91HRpM+AgNr75yk8PFzfffddjbefNWuWSkpKNHToUPdD0g4dOiRJ+u6777R//37FxMSUm4hclYceekjdu3fXoEGD9NFHH+nIkSN64oknFBwcrOnTp+vWW29VZmam/P3///+I3333nRo3blzlZGXAFxA+gBooKirSgAED9OWXX+q9996r8PJ/fn6+JOns2bPl1pWWlurMmTO1Vk9V/xMODAzUgAEDNGDAAJWVlel3v/udFi5cqIceekhXXHGFV/t59dVXddlll+n111/32Gd6errXNV9++eWSfghISUlJlfaLi4uT9P9XSn5s79691dpXu3btlJOT43WN5+Tm5uq7774rN6FY+uFhbDNmzND27durPeH2n//8p55//nllZ2dLkg4fPqxmzZopODhYkhQTE6OSkhJ98803ioqKcm+Xk5Oj9u3b1/g4gPqC2y6Al86ePauhQ4cqKytLy5YtU2JiYoX9rrjiCvn7++vll1/2uOR/6NAhffTRR7r66qtrraYmTZpIkgoKCjzajx075rHs7++vLl26SFKNHtV97rbSj49n8+bNysrK8nqs5ORkORwOzZgxQ6WlpeXWn/tWT4sWLdStWze98MILHrd21q5dqy+++KJa+0pMTNTOnTtr/Hjy3//+91q+fLnH69zTRkeNGqXly5erbdu21R5vwoQJuuuuu9SpUydJUlRUlL755hsdP35c0g/fSAoICFDz5s09ttu2bZv7WzKAL+PKB+ClKVOm6K233tKAAQN0/Phxvfjiix7rzz3/ISIiQmPGjNHf/vY39e3bV4MHD9aJEyc0b948nT59WtOmTau1mnr06CHphw/J5ORkNWjQQMOGDdNdd92l48ePq0+fPmrVqpUOHDigOXPmqFu3bjX6H/Rtt92m119/Xb/+9a+VkpKinJwcLViwQB06dNDJkye9GsvhcGj+/Pm644471L17dw0bNkwRERHKzc3VypUrdd111+kvf/mLJCkjI0MpKSm6/vrrNWbMGB0/flxz5sxRx44dq7XfgQMH6rHHHtOGDRvUr18/j3V/+ctfVFBQoMOHD0uSVqxY4b6lMn78eDmdTnXv3l3du3f32O7c7ZeOHTtq0KBB1T7uZcuW6bPPPtNrr73mbktMTFRUVJSGDBmiwYMHa9asWRo8eLDHHKLs7GwdP35cAwcOrPa+gHqrbr9sA/iem266yeOrsz99/VhpaamZM2eO6datm2natKlp2rSpueWWW8z7779frX1V96u2Z86cMePHjzcRERHGz8/PXcerr75q+vXrZyIjI01gYKCJjY019957r/upnt7uu6yszMyYMcPExcWZoKAgc/XVV5u33367wieN6jxftT3ngw8+MMnJycbpdJrg4GBz+eWXm1GjRpmtW7d69HvttddM+/btTVBQkOnQoYN5/fXXvXq8epcuXczYsWMrPM7KzmVVT2P15qu253z//fcmLi7OPPvss+XWbdmyxXTv3t2EhISYAQMGmKNHj3qsf/DBB01sbCyPV8clwc+YC5gCDgA+YsmSJUpLS1Nubq5CQ0PruhyvFBcXq02bNpo6dWqVD6cDfAVzPgD8LIwYMUKxsbHuR9H7kszMTDVs2FD33XdfXZcC1AqufAAAAKu48gEAAKwifAAAAKsIHwAAwCrCBwAAsKrePWSsrKxMhw8fVkhICL88CQAAH2GM0YkTJxQTE+PxO4kqUu/Cx+HDh2v8WyMBAEDdOnjwoFq1alVln3oXPkJCQiT9ULzD4ajjagAAQHUUFhaqdevW7s/xqtS78HHuVovD4SB8AADgY6ozZYIJpwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArAqo6wIA1FybqSvrugSv7Z+ZUtclAKhjXPkAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVXoWPRx55RH5+fh6vdu3audcXFRUpLS1N4eHhatq0qVJTU5Wfn1/rRQMAAN/l9ZWPjh076siRI+7Xxx9/7F43adIkrVixQsuWLdOGDRt0+PBhDR48uFYLBgAAvi3A6w0CAhQdHV2u3eVyadGiRVq6dKn69OkjScrMzFT79u21adMm9erV68KrBQAAPs/rKx9fffWVYmJidNlll2nEiBHKzc2VJGVnZ6u0tFRJSUnuvu3atVNsbKyysrJqr2IAAODTvLrykZCQoMWLFys+Pl5HjhzR9OnTdcMNN2jnzp3Ky8tTYGCgQkNDPbaJiopSXl5epWMWFxeruLjYvVxYWOjdEQAAAJ/iVfjo37+/+89dunRRQkKC4uLi9Morr6hRo0Y1KiAjI0PTp0+v0bYAAMD3XNBXbUNDQ3XVVVdp3759io6OVklJiQoKCjz65OfnVzhH5Jxp06bJ5XK5XwcPHryQkgAAQD13QeHj5MmT+vrrr9WiRQv16NFDDRs21Lp169zr9+7dq9zcXCUmJlY6RlBQkBwOh8cLAABcury67fKf//mfGjBggOLi4nT48GGlp6erQYMGGj58uJxOp8aOHavJkycrLCxMDodD48ePV2JiIt90AQAAbl6Fj0OHDmn48OE6duyYIiIidP3112vTpk2KiIiQJM2ePVv+/v5KTU1VcXGxkpOTNW/evItSOAAA8E1+xhhT10X8WGFhoZxOp1wuF7dggPNoM3VlXZfgtf0zU+q6BAAXgTef3/xuFwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWBdR1AQB+XtpMXVnXJXht/8yUui4BuKRw5QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYdUHhY+bMmfLz89PEiRPdbUVFRUpLS1N4eLiaNm2q1NRU5efnX2idAADgElHj8LFlyxYtXLhQXbp08WifNGmSVqxYoWXLlmnDhg06fPiwBg8efMGFAgCAS0ONwsfJkyc1YsQIPffcc2rWrJm73eVyadGiRXr66afVp08f9ejRQ5mZmdq4caM2bdpUa0UDAADfVaPwkZaWppSUFCUlJXm0Z2dnq7S01KO9Xbt2io2NVVZW1oVVCgAALgkB3m7w0ksvadu2bdqyZUu5dXl5eQoMDFRoaKhHe1RUlPLy8iocr7i4WMXFxe7lwsJCb0sCAAA+xKsrHwcPHtSECRP097//XcHBwbVSQEZGhpxOp/vVunXrWhkXAADUT16Fj+zsbB09elTdu3dXQECAAgICtGHDBj377LMKCAhQVFSUSkpKVFBQ4LFdfn6+oqOjKxxz2rRpcrlc7tfBgwdrfDAAAKD+8+q2S9++ffX55597tI0ePVrt2rXTgw8+qNatW6thw4Zat26dUlNTJUl79+5Vbm6uEhMTKxwzKChIQUFBNSwfAAD4Gq/CR0hIiDp16uTR1qRJE4WHh7vbx44dq8mTJyssLEwOh0Pjx49XYmKievXqVXtVAwAAn+X1hNPzmT17tvz9/ZWamqri4mIlJydr3rx5tb0bAADgo/yMMaaui/ixwsJCOZ1OuVwuORyOui4HqNfaTF1Z1yX8LOyfmVLXJQD1njef3/xuFwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABglVfhY/78+erSpYscDoccDocSExO1atUq9/qioiKlpaUpPDxcTZs2VWpqqvLz82u9aAAA4Lu8Ch+tWrXSzJkzlZ2dra1bt6pPnz4aOHCgdu3aJUmaNGmSVqxYoWXLlmnDhg06fPiwBg8efFEKBwAAvsnPGGMuZICwsDA9+eSTuv322xUREaGlS5fq9ttvlyTt2bNH7du3V1ZWlnr16lWt8QoLC+V0OuVyueRwOC6kNOCS12bqyrou4Wdh/8yUui4BqPe8+fyu8ZyPs2fP6qWXXtKpU6eUmJio7OxslZaWKikpyd2nXbt2io2NVVZWVk13AwAALjEB3m7w+eefKzExUUVFRWratKmWL1+uDh06aMeOHQoMDFRoaKhH/6ioKOXl5VU6XnFxsYqLi93LhYWF3pYEAAB8iNdXPuLj47Vjxw5t3rxZ999/v0aOHKkvvviixgVkZGTI6XS6X61bt67xWAAAoP7zOnwEBgbqiiuuUI8ePZSRkaGuXbvqmWeeUXR0tEpKSlRQUODRPz8/X9HR0ZWON23aNLlcLvfr4MGDXh8EAADwHRf8nI+ysjIVFxerR48eatiwodatW+det3fvXuXm5ioxMbHS7YOCgtxf3T33AgAAly6v5nxMmzZN/fv3V2xsrE6cOKGlS5dq/fr1WrNmjZxOp8aOHavJkycrLCxMDodD48ePV2JiYrW/6QIAAC59XoWPo0eP6s4779SRI0fkdDrVpUsXrVmzRr/4xS8kSbNnz5a/v79SU1NVXFys5ORkzZs376IUDgAAfNMFP+ejtvGcD6D6eM6HHTznAzg/K8/5AAAAqAnCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsCqjrAgCgvmszdWVdl1Aj+2em1HUJQIW48gEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACs8ip8ZGRk6Nprr1VISIgiIyM1aNAg7d2716NPUVGR0tLSFB4erqZNmyo1NVX5+fm1WjQAAPBdXoWPDRs2KC0tTZs2bdLatWtVWlqqfv366dSpU+4+kyZN0ooVK7Rs2TJt2LBBhw8f1uDBg2u9cAAA4JsCvOm8evVqj+XFixcrMjJS2dnZuvHGG+VyubRo0SItXbpUffr0kSRlZmaqffv22rRpk3r16lV7lQMAAJ90QXM+XC6XJCksLEySlJ2drdLSUiUlJbn7tGvXTrGxscrKyrqQXQEAgEuEV1c+fqysrEwTJ07Uddddp06dOkmS8vLyFBgYqNDQUI++UVFRysvLq3Cc4uJiFRcXu5cLCwtrWhIAAPABNQ4faWlp2rlzpz7++OMLKiAjI0PTp0+/oDGA2tBm6sq6LgEAfhZqdNtl3Lhxevvtt/XBBx+oVatW7vbo6GiVlJSooKDAo39+fr6io6MrHGvatGlyuVzu18GDB2tSEgAA8BFehQ9jjMaNG6fly5fr/fffV9u2bT3W9+jRQw0bNtS6devcbXv37lVubq4SExMrHDMoKEgOh8PjBQAALl1e3XZJS0vT0qVL9eabbyokJMQ9j8PpdKpRo0ZyOp0aO3asJk+erLCwMDkcDo0fP16JiYl80wUAAEjyMnzMnz9fknTzzTd7tGdmZmrUqFGSpNmzZ8vf31+pqakqLi5WcnKy5s2bVyvFAgAA3+dV+DDGnLdPcHCw5s6dq7lz59a4KAAAcOnid7sAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwLqugAAwMXRZurKui7Ba/tnptR1CbCAKx8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArPI6fHz44YcaMGCAYmJi5OfnpzfeeMNjvTFGDz/8sFq0aKFGjRopKSlJX331VW3VCwAAfJzX4ePUqVPq2rWr5s6dW+H6J554Qs8++6wWLFigzZs3q0mTJkpOTlZRUdEFFwsAAHyf10847d+/v/r371/hOmOM/vznP+uPf/yjBg4cKEn63//9X0VFRemNN97QsGHDLqxaAADg82p1zkdOTo7y8vKUlJTkbnM6nUpISFBWVlZt7goAAPioWv3dLnl5eZKkqKgoj/aoqCj3up8qLi5WcXGxe7mwsLA2SwIAAPVMnX/bJSMjQ06n0/1q3bp1XZcEAAAuoloNH9HR0ZKk/Px8j/b8/Hz3up+aNm2aXC6X+3Xw4MHaLAkAANQztRo+2rZtq+joaK1bt87dVlhYqM2bNysxMbHCbYKCguRwODxeAADg0uX1nI+TJ09q37597uWcnBzt2LFDYWFhio2N1cSJE/U///M/uvLKK9W2bVs99NBDiomJ0aBBg2qzbgAA4KO8Dh9bt27VLbfc4l6ePHmyJGnkyJFavHix/vCHP+jUqVO65557VFBQoOuvv16rV69WcHBw7VUNAAB8lp8xxtR1ET9WWFgop9Mpl8vFLRhY1WbqyrouAfjZ2z8zpa5LQA158/ld5992AQAAPy+EDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYVau/1RY4h2dmAAAqw5UPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVQF1XQAAAOe0mbqyrkvw2v6ZKXVdgs/hygcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArAqo6wJwfm2mrqzrEgAAqDVc+QAAAFZdtPAxd+5ctWnTRsHBwUpISNCnn356sXYFAAB8yEUJHy+//LImT56s9PR0bdu2TV27dlVycrKOHj16MXYHAAB8yEUJH08//bTuvvtujR49Wh06dNCCBQvUuHFjPf/88xdjdwAAwIfUevgoKSlRdna2kpKS/n8n/v5KSkpSVlZWbe8OAAD4mFr/tsu3336rs2fPKioqyqM9KipKe/bsKde/uLhYxcXF7mWXyyVJKiwsrO3SfFZZ8fd1XQIAoBJ8Xv3g3M/BGHPevnX+VduMjAxNnz69XHvr1q3roBoAALzj/HNdV1C/nDhxQk6ns8o+tR4+mjdvrgYNGig/P9+jPT8/X9HR0eX6T5s2TZMnT3Yvl5WV6fjx4woPD5efn1+t1lZYWKjWrVvr4MGDcjgctTp2fXCpH5906R8jx+f7LvVj5Ph838U6RmOMTpw4oZiYmPP2rfXwERgYqB49emjdunUaNGiQpB8Cxbp16zRu3Lhy/YOCghQUFOTRFhoaWttleXA4HJfsXyrp0j8+6dI/Ro7P913qx8jx+b6LcYznu+JxzkW57TJ58mSNHDlS11xzjXr27Kk///nPOnXqlEaPHn0xdgcAAHzIRQkfQ4cO1TfffKOHH35YeXl56tatm1avXl1uEioAAPj5uWgTTseNG1fhbZa6FBQUpPT09HK3eS4Vl/rxSZf+MXJ8vu9SP0aOz/fVh2P0M9X5TgwAAEAt4RfLAQAAqwgfAADAKsIHAACwivABAACsuqTCx5/+9Cf17t1bjRs3rvRBZbm5uUpJSVHjxo0VGRmpBx54QGfOnKly3OPHj2vEiBFyOBwKDQ3V2LFjdfLkyYtwBN5Zv369/Pz8Knxt2bKl0u1uvvnmcv3vu+8+i5VXX5s2bcrVOnPmzCq3KSoqUlpamsLDw9W0aVOlpqaWe+JufbF//36NHTtWbdu2VaNGjXT55ZcrPT1dJSUlVW5Xn8/h3Llz1aZNGwUHByshIUGffvpplf2XLVumdu3aKTg4WJ07d9Y777xjqVLvZWRk6Nprr1VISIgiIyM1aNAg7d27t8ptFi9eXO5cBQcHW6rYO4888ki5Wtu1a1flNr50/qSK/03x8/NTWlpahf3r+/n78MMPNWDAAMXExMjPz09vvPGGx3pjjB5++GG1aNFCjRo1UlJSkr766qvzjuvt+9hbl1T4KCkp0ZAhQ3T//fdXuP7s2bNKSUlRSUmJNm7cqBdeeEGLFy/Www8/XOW4I0aM0K5du7R27Vq9/fbb+vDDD3XPPfdcjEPwSu/evXXkyBGP11133aW2bdvqmmuuqXLbu+++22O7J554wlLV3nv00Uc9ah0/fnyV/SdNmqQVK1Zo2bJl2rBhgw4fPqzBgwdbqtY7e/bsUVlZmRYuXKhdu3Zp9uzZWrBggf7rv/7rvNvWx3P48ssva/LkyUpPT9e2bdvUtWtXJScn6+jRoxX237hxo4YPH66xY8dq+/btGjRokAYNGqSdO3darrx6NmzYoLS0NG3atElr165VaWmp+vXrp1OnTlW5ncPh8DhXBw4csFSx9zp27OhR68cff1xpX187f5K0ZcsWj+Nbu3atJGnIkCGVblOfz9+pU6fUtWtXzZ07t8L1TzzxhJ599lktWLBAmzdvVpMmTZScnKyioqJKx/T2fVwj5hKUmZlpnE5nufZ33nnH+Pv7m7y8PHfb/PnzjcPhMMXFxRWO9cUXXxhJZsuWLe62VatWGT8/P/Pvf/+71mu/ECUlJSYiIsI8+uijVfa76aabzIQJE+wUdYHi4uLM7Nmzq92/oKDANGzY0Cxbtszdtnv3biPJZGVlXYQKa98TTzxh2rZtW2Wf+noOe/bsadLS0tzLZ8+eNTExMSYjI6PC/r/5zW9MSkqKR1tCQoK59957L2qdteXo0aNGktmwYUOlfSr796g+Sk9PN127dq12f18/f8YYM2HCBHP55ZebsrKyCtf70vmTZJYvX+5eLisrM9HR0ebJJ590txUUFJigoCDzj3/8o9JxvH0f18QldeXjfLKystS5c2ePJ60mJyersLBQu3btqnSb0NBQjysJSUlJ8vf31+bNmy96zd546623dOzYsWo9xv7vf/+7mjdvrk6dOmnatGn6/vvvLVRYMzNnzlR4eLiuvvpqPfnkk1XeJsvOzlZpaamSkpLcbe3atVNsbKyysrJslHvBXC6XwsLCztuvvp3DkpISZWdne/zs/f39lZSUVOnPPisry6O/9MN70pfOlaTznq+TJ08qLi5OrVu31sCBAyv996Y++OqrrxQTE6PLLrtMI0aMUG5ubqV9ff38lZSU6MUXX9SYMWOq/EWmvnT+fiwnJ0d5eXke58jpdCohIaHSc1ST93FNXLQnnNZHeXl55R7xfm45Ly+v0m0iIyM92gICAhQWFlbpNnVl0aJFSk5OVqtWrars9x//8R+Ki4tTTEyMPvvsMz344IPau3evXn/9dUuVVt/vf/97de/eXWFhYdq4caOmTZumI0eO6Omnn66wf15engIDA8vN+YmKiqp356si+/bt05w5czRr1qwq+9XHc/jtt9/q7NmzFb7H9uzZU+E2lb0nfeFclZWVaeLEibruuuvUqVOnSvvFx8fr+eefV5cuXeRyuTRr1iz17t1bu3btOu971baEhAQtXrxY8fHxOnLkiKZPn64bbrhBO3fuVEhISLn+vnz+JOmNN95QQUGBRo0aVWkfXzp/P3XuPHhzjmryPq6Jeh8+pk6dqscff7zKPrt37z7vpChfUpNjPnTokNasWaNXXnnlvOP/eL5K586d1aJFC/Xt21dff/21Lr/88poXXk3eHN/kyZPdbV26dFFgYKDuvfdeZWRk1OvHH9fkHP773//WrbfeqiFDhujuu++uctu6PoeQ0tLStHPnzirnREhSYmKiEhMT3cu9e/dW+/bttXDhQj322GMXu0yv9O/f3/3nLl26KCEhQXFxcXrllVc0duzYOqzs4li0aJH69+9f5a+A96Xz50vqffiYMmVKlalUki677LJqjRUdHV1uxu65b0FER0dXus1PJ9mcOXNGx48fr3SbC1WTY87MzFR4eLh+9atfeb2/hIQEST/8r9vGB9eFnNOEhASdOXNG+/fvV3x8fLn10dHRKikpUUFBgcfVj/z8/It2viri7TEePnxYt9xyi3r37q2//vWvXu/P9jmsSPPmzdWgQYNy3yyq6mcfHR3tVf/6Yty4ce7J597+77dhw4a6+uqrtW/fvotUXe0JDQ3VVVddVWmtvnr+JOnAgQN67733vL5a6Evn79x5yM/PV4sWLdzt+fn56tatW4Xb1OR9XCO1NnukHjnfhNP8/Hx328KFC43D4TBFRUUVjnVuwunWrVvdbWvWrKlXE07LyspM27ZtzZQpU2q0/ccff2wkmX/+85+1XFnte/HFF42/v785fvx4hevPTTh99dVX3W179uyp1xNODx06ZK688kozbNgwc+bMmRqNUV/OYc+ePc24cePcy2fPnjUtW7ascsLpbbfd5tGWmJhYbycslpWVmbS0NBMTE2O+/PLLGo1x5swZEx8fbyZNmlTL1dW+EydOmGbNmplnnnmmwvW+dv5+LD093URHR5vS0lKvtqvP50+VTDidNWuWu83lclVrwqk37+Ma1VprI9UDBw4cMNu3bzfTp083TZs2Ndu3bzfbt283J06cMMb88JemU6dOpl+/fmbHjh1m9erVJiIiwkybNs09xubNm018fLw5dOiQu+3WW281V199tdm8ebP5+OOPzZVXXmmGDx9u/fgq89577xlJZvfu3eXWHTp0yMTHx5vNmzcbY4zZt2+fefTRR83WrVtNTk6OefPNN81ll11mbrzxRttln9fGjRvN7NmzzY4dO8zXX39tXnzxRRMREWHuvPNOd5+fHp8xxtx3330mNjbWvP/++2br1q0mMTHRJCYm1sUhnNehQ4fMFVdcYfr27WsOHTpkjhw54n79uI+vnMOXXnrJBAUFmcWLF5svvvjC3HPPPSY0NNT9DbM77rjDTJ061d3/k08+MQEBAWbWrFlm9+7dJj093TRs2NB8/vnndXUIVbr//vuN0+k069ev9zhX33//vbvPT49x+vTpZs2aNebrr7822dnZZtiwYSY4ONjs2rWrLg6hSlOmTDHr1683OTk55pNPPjFJSUmmefPm5ujRo8YY3z9/55w9e9bExsaaBx98sNw6Xzt/J06ccH/WSTJPP/202b59uzlw4IAxxpiZM2ea0NBQ8+abb5rPPvvMDBw40LRt29acPn3aPUafPn3MnDlz3Mvnex/XhksqfIwcOdJIKvf64IMP3H32799v+vfvbxo1amSaN29upkyZ4pF8P/jgAyPJ5OTkuNuOHTtmhg8fbpo2bWocDocZPXq0O9DUB8OHDze9e/eucF1OTo7HzyA3N9fceOONJiwszAQFBZkrrrjCPPDAA8blclmsuHqys7NNQkKCcTqdJjg42LRv397MmDHD4yrVT4/PGGNOnz5tfve735lmzZqZxo0bm1//+tceH+b1SWZmZoV/Z398UdLXzuGcOXNMbGysCQwMND179jSbNm1yr7vpppvMyJEjPfq/8sor5qqrrjKBgYGmY8eOZuXKlZYrrr7KzlVmZqa7z0+PceLEie6fR1RUlPnlL39ptm3bZr/4ahg6dKhp0aKFCQwMNC1btjRDhw41+/btc6/39fN3zpo1a4wks3fv3nLrfO38nfvM+unr3DGUlZWZhx56yERFRZmgoCDTt2/fcscdFxdn0tPTPdqqeh/XBj9jjKm9mzgAAABV+1k95wMAANQ9wgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8ALBi1KhRGjRoUF2XAaAe4CFjAKxwuVwyxnj8wj8AP0+EDwAAYBW3XQDUqldffVWdO3dWo0aNFB4erqSkJJ06darcbZcTJ05oxIgRatKkiVq0aKHZs2fr5ptv1sSJE+usdgB2ED4A1JojR45o+PDhGjNmjHbv3q3169dr8ODBqugC6+TJk/XJJ5/orbfe0tq1a/XRRx9p27ZtdVA1ANsC6roAAJeOI0eO6MyZMxo8eLDi4uIkSZ07dy7X78SJE3rhhRe0dOlS9e3bV5KUmZmpmJgYq/UCqBtc+QBQa7p27aq+ffuqc+fOGjJkiJ577jl999135fr961//UmlpqXr27Oluczqdio+Pt1kugDpC+ABQaxo0aKC1a9dq1apV6tChg+bMmaP4+Hjl5OTUdWkA6hHCB4Ba5efnp+uuu07Tp0/X9u3bFRgYqOXLl3v0ueyyy9SwYUNt2bLF3eZyufTll1/aLhdAHWDOB4Bas3nzZq1bt079+vVTZGSkNm/erG+++Ubt27fXZ5995u4XEhKikSNH6oEHHlBYWJgiIyOVnp4uf39/+fn51eERALCBKx8Aao3D4dCHH36oX/7yl7rqqqv0xz/+UU899ZT69+9fru/TTz+txMRE3XbbbUpKStJ1112n9u3bKzg4uA4qB2ATDxkDUC+cOnVKLVu21FNPPaWxY8fWdTkALiJuuwCoE9u3b9eePXvUs2dPuVwuPfroo5KkgQMH1nFlAC42wgeAOjNr1izt3btXgYGB6tGjhz766CM1b968rssCcJFx2wUAAFjFhFMAAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABg1f8BFGTB0DTjyXsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nfailed = np.sum(np.isnan(sigs_migrad))\n",
    "plt.title(f\"{nfailed} fits failed ({nfailed / len(sigs_migrad) * 100:.0f} %)\")\n",
    "plt.hist(sigs_migrad, bins=10, range=(-10, 10))\n",
    "plt.xlabel(\"sig\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "432b611e",
   "metadata": {},
   "source": [
    "The distribution of the signal amplitude looks fairly gaussian which is good, but the fit failed to converge in a few cases due to the problem just described. Simply discarding these cases is not acceptable, it would distort conclusions drawn from the distribution of the test statistic, which is commonly needed to set limits or to compute the p-value for an observed amplitude.\n",
    "\n",
    "We can repair this by placing a limit on the signal amplitude. This is a suitable solution, although it will bias the signal amplitude and change the shape of the distribution of the test statistic. \n",
    "\n",
    "An alternative is to perform a constrained minimization, which allows us to directly add a condition to the fit that the model density must be positive at each data point. We merely need to replace the call `m.migrad` with the call `m.scipy` and pass the (non-linear) constraint. An appropriate algorithm is automatically selected which performs a constrained minimization. The SciPy minimizers are fully integrated into Minuit, which means that Minuit computes an EDM value for the minimum and parameter uncertainties."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d364cb83",
   "metadata": {},
   "source": [
    "## SciPy constrained fit on toys \n",
    "\n",
    "We run SciPy with the constraint on the same simulated samples on which we ran Migrad before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "young-ocean",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 failed\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import NonlinearConstraint\n",
    "\n",
    "@joblib.delayed\n",
    "def compute(itry):\n",
    "    rng = np.random.default_rng(itry)\n",
    "    x = rng.uniform(0, 1, size=35)\n",
    "    cost = ExtendedUnbinnedNLL(x, model)\n",
    "    m = Minuit(cost, b0=n, b1=n, b2=n, sig=0, mu=0.5, sigma=0.05)\n",
    "    m.limits[\"b0\", \"b1\", \"b2\"] = (0, None)\n",
    "    m.fixed[\"mu\", \"sigma\"] = True\n",
    "    m.scipy(constraints=NonlinearConstraint(lambda *par: model(x, *par)[1], 0, np.inf))\n",
    "    return m.values[\"sig\"] if m.valid else np.nan\n",
    "\n",
    "sigs_constrained = joblib.Parallel(-1)(compute(i) for i in range(200))\n",
    "\n",
    "print(np.sum(np.isnan(sigs_constrained)), \"failed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0ce87a47",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAHHCAYAAAAf2DoOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9TUlEQVR4nO3deVxWdf7//yeILCKLoLIkKC6BpqZpEZpLSpGp4zZZjiUqk1bUuGbRTG5NUWZq+XPJxsAWv6alLTOl4+6oaOZSZmpquCVgaoAr6/v3RzevT1cseiEeQB/32+263bje533OeZ1zxOvJOe9zLidjjBEAAIBFnCu6AAAAcHMhfAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8ALiiwYMHq0GDBhWy7uTkZDk5Oenw4cNX7Lt8+XK1atVK7u7ucnJyUmZmpmW1F7fuqzVx4kQ5OTnZtTVo0ECDBw8ut/oOHz4sJycnJScnl9sygbIifOCmk5OTo+eee07BwcHy8PBQZGSkVq5cWdFlXZMTJ05o4sSJ2rVrV0WXUmFOnz6t/v37y8PDQ7NmzdL7778vT0/PIv0uXLigiRMnat26dZavG8BvXCq6AMBqgwcP1scff6yRI0eqSZMmSk5O1oMPPqi1a9fqnnvuqejyyuTEiROaNGmSGjRooFatWpX78t955x0VFhaW+3LL07Zt23T27Fm99NJLio6OtrX/sfYLFy5o0qRJkqTOnTtf13VfrX/84x96/vnny6UWoCogfOCm8vXXX2vRokV6/fXXNXbsWEnSoEGD1Lx5c40bN06bN2+u4AqtceHCBdWoUeOq+1evXv06VlM+Tp48KUny9fW1a7ei9pLWfbVcXFzk4sJ/x7h5cNkFN5WPP/5Y1apV07Bhw2xt7u7uiouLU0pKio4dO3bFZWzdulUPPvigatWqJU9PT7Vs2VJvvvmmXZ81a9aoQ4cO8vT0lK+vr3r16qW9e/fa9bl8nf/gwYMaPHiwfH195ePjoyFDhujChQt2fVeuXKl77rlHvr6+qlmzpsLDw/XCCy9IktatW6c777xTkjRkyBA5OTnZXdvv3Lmzmjdvru3bt6tjx46qUaOGbd7PPvtM3bt3V3BwsNzc3NSoUSO99NJLKigosFv/H8dNXB4/MHXqVM2bN0+NGjWSm5ub7rzzTm3btq3IPtu3b5/+/Oc/y8/PT+7u7mrbtq0+//zzIv327NmjLl26yMPDQ/Xq1dM///nPqzrj0rlzZ8XGxkqS7rzzTjk5OdnGS/y+9sOHD6tOnTqSpEmTJtn21cSJEyVJ6enpGjJkiOrVqyc3NzcFBQWpV69epY43KW3d//vf//TQQw8pNDRUbm5uCgkJ0ahRo3Tx4kW7ZRQ35qM4mZmZGjlypEJCQuTm5qbGjRvrtddeK7KPLo918fHxka+vr2JjYx0agwJcb0Rt3FR27typW2+9Vd7e3nbtd911lyRp165dCgkJKXH+lStXqkePHgoKCtKIESMUGBiovXv36t///rdGjBghSVq1apW6deumhg0bauLEibp48aJmzpyp9u3ba8eOHUUGP/bv319hYWFKTEzUjh079K9//Ut169bVa6+9Jum3D+QePXqoZcuWmjx5stzc3HTw4EFt2rRJktS0aVNNnjxZ48eP17Bhw9ShQwdJUrt27WzrOH36tLp166ZHHnlEjz76qAICAiT9NpizZs2aGj16tGrWrKk1a9Zo/Pjxys7O1uuvv37F/blw4UKdPXtWw4cPl5OTk6ZMmaK+ffvqp59+sp1x2LNnj9q3b69bbrlFzz//vDw9PbV48WL17t1bn3zyifr06SPptw/+e++9V/n5+bZ+8+bNk4eHxxXr+Pvf/67w8HDNmzdPkydPVlhYmBo1alSkX506dTRnzhw9+eST6tOnj/r27StJatmypSSpX79+2rNnj5555hk1aNBAJ0+e1MqVK3X06NESB62Wtu4lS5bowoULevLJJ+Xv76+vv/5aM2fO1PHjx7VkyZIrbtfvXbhwQZ06ddLPP/+s4cOHKzQ0VJs3b1ZCQoLS0tI0Y8YMSZIxRr169dLGjRv1xBNPqGnTplq2bJktIAGVggFuIrfddpvp0qVLkfY9e/YYSWbu3Lklzpufn2/CwsJM/fr1za+//mo3rbCw0PZzq1atTN26dc3p06dtbd9++61xdnY2gwYNsrVNmDDBSDJDhw61W1afPn2Mv7+/7f306dONJPPLL7+UWNu2bduMJJOUlFRkWqdOnUrctgsXLhRpGz58uKlRo4a5dOmSrS02NtbUr1/f9j41NdVIMv7+/ubMmTO29s8++8xIMl988YWtrWvXrqZFixZ2yyssLDTt2rUzTZo0sbWNHDnSSDJbt261tZ08edL4+PgYSSY1NbXE7TfGmKSkJCPJbNu2za79j7X/8ssvRpKZMGGCXb9ff/3VSDKvv/56qetxZN3F7d/ExETj5ORkjhw5Ymu7/G/h9+rXr29iY2Nt71966SXj6elpfvzxR7t+zz//vKlWrZo5evSoMcaYTz/91EgyU6ZMsfXJz883HTp0KPHfCGA1LrvgpnLx4kW5ubkVaXd3d7dNL8nOnTuVmpqqkSNHFrm2f/mUeVpamnbt2qXBgwfLz8/PNr1ly5a677779OWXXxZZ7hNPPGH3vkOHDjp9+rSys7Ml/d84gs8++6zMgz7d3Nw0ZMiQIu2/P6tw9uxZnTp1Sh06dNCFCxe0b9++Ky734YcfVq1atexql6SffvpJknTmzBmtWbNG/fv3ty3/1KlTOn36tGJiYnTgwAH9/PPPkqQvv/xSd999t+0slPTbmYqBAweWaZsd5eHhIVdXV61bt06//vpruS3zsvPnz+vUqVNq166djDHauXOnQ8tasmSJOnTooFq1atn246lTpxQdHa2CggJt2LBB0m/70cXFRU8++aRt3mrVqumZZ54pl20CygPhAzcVDw8P5eTkFGm/dOmSbXpJDh06JElq3rx5iX2OHDkiSQoPDy8yrWnTpjp16pTOnz9v1x4aGmr3/vKH+eUPwIcffljt27fXX//6VwUEBOiRRx7R4sWLHQoit9xyi1xdXYu079mzR3369JGPj4+8vb1Vp04dPfroo5KkrKysKy73SrUfPHhQxhi9+OKLqlOnjt1rwoQJkv5vsOaRI0fUpEmTIusobl9eD25ubnrttdf01VdfKSAgQB07dtSUKVOUnp5e5mUePXrUFkRr1qypOnXqqFOnTpKubv/+3oEDB7R8+fIi+/Hy3TW/349BQUGqWbOm3fxW7UfgajDmAzeVoKAg21/av5eWliZJCg4OtrokVatWrdh2Y4yk3wLRhg0btHbtWv3nP//R8uXL9dFHH6lLly7673//W+L8v1dcqMrMzFSnTp3k7e2tyZMnq1GjRnJ3d9eOHTv03HPPXVW4uVLtl5cxduxYxcTEFNu3cePGV1yPVUaOHKmePXvq008/1YoVK/Tiiy8qMTFRa9asUevWrR1aVkFBge677z6dOXNGzz33nCIiIuTp6amff/5ZgwcPdvgsVmFhoe677z6NGzeu2Om33nqrQ8sDKhLhAzeVVq1aae3atcrOzrYbdLp161bb9JJcHkT4/fffl/gsh/r160uS9u/fX2Tavn37VLt27TI9fMrZ2Vldu3ZV165dNW3aNL3yyiv6+9//rrVr1yo6Ovqq7pT4o3Xr1un06dNaunSpOnbsaGtPTU11eFkladiwoaTfbne90vMv6tevrwMHDhRpL25fXosr7atGjRppzJgxGjNmjA4cOKBWrVrpjTfe0AcffODQenbv3q0ff/xRCxYs0KBBg2ztZX2gXaNGjXTu3Lmr2o+rV6/WuXPn7M5+lPd+BK4Fl11wU/nzn/+sgoICzZs3z9aWk5OjpKQkRUZGlnqnyx133KGwsDDNmDGjyG2Ll//SDwoKUqtWrbRgwQK7Pt9//73++9//6sEHH3S45jNnzhRpuxySLl9CuhxoHLmd8vJZi8u1S1Jubq5mz57tcI0lqVu3rjp37qy3337bdnbp93755Rfbzw8++KC2bNmir7/+2m76hx9+WG71SLI93+SP++rChQu2y2+XNWrUSF5eXsVeqruS4vavMabIbdlXq3///kpJSdGKFSuKTMvMzFR+fr6k3/Zjfn6+5syZY5teUFCgmTNnlmm9wPXAmQ/cVCIjI/XQQw8pISFBJ0+eVOPGjbVgwQIdPnxY8+fPL3VeZ2dnzZkzRz179lSrVq00ZMgQBQUFad++fdqzZ4/tQ+H1119Xt27dFBUVpbi4ONuttj4+PrbnSThi8uTJ2rBhg7p376769evr5MmTmj17turVq2d7ImujRo3k6+uruXPnysvLS56enoqMjFRYWFiJy23Xrp1q1aql2NhY/e1vf5OTk5Pef/99uw/L8jBr1izdc889atGihR5//HE1bNhQGRkZSklJ0fHjx/Xtt99KksaNG6f3339fDzzwgEaMGGG71bZ+/fr67rvvyq0eDw8PNWvWTB999JFuvfVW+fn5qXnz5srPz1fXrl3Vv39/NWvWTC4uLlq2bJkyMjL0yCOPOLyeiIgINWrUSGPHjtXPP/8sb29vffLJJ2UezPrss8/q888/V48ePTR48GC1adNG58+f1+7du/Xxxx/r8OHDql27tnr27Kn27dvr+eef1+HDh9WsWTMtXbrU4TEmwHVVYffZABXk4sWLZuzYsSYwMNC4ubmZO++80yxfvvyq59+4caO57777jJeXl/H09DQtW7Y0M2fOtOuzatUq0759e+Ph4WG8vb1Nz549zQ8//GDX5/LtlX+8hfbybZuXby1dvXq16dWrlwkODjaurq4mODjYDBgwoMgtl5999plp1qyZcXFxsbulslOnTua2224rdls2bdpk7r77buPh4WGCg4PNuHHjzIoVK4wks3btWlu/km61Le62VBVzG+uhQ4fMoEGDTGBgoKlevbq55ZZbTI8ePczHH39s1++7774znTp1Mu7u7uaWW24xL730kpk/f3653mprjDGbN282bdq0Ma6urrZ6T506ZeLj401ERITx9PQ0Pj4+JjIy0ixevLjU9Za27h9++MFER0ebmjVrmtq1a5vHH3/cfPvtt0Vueb2aW22NMebs2bMmISHBNG7c2Li6upratWubdu3amalTp5rc3Fxbv9OnT5vHHnvMeHt7Gx8fH/PYY4+ZnTt3cqstKg0nY8r5zxwAAIBSMOYDAABYivABAAAsRfgAAACWInwAAABLET4AAIClCB8AAMBSle4hY4WFhTpx4oS8vLzK9MhoAABgPWOMzp49q+DgYDk7l35uo9KFjxMnTpT6iGsAAFB5HTt2TPXq1Su1T6ULH15eXpJ+K/73X/wFAAAqr+zsbIWEhNg+x0tT6cLH5Ust3t7ehA8AAKqYqxkywYBTAABgKcIHAACwFOEDAABYqtKN+QAAoDQFBQXKy8ur6DJuSq6urle8jfZqED4AAFWCMUbp6enKzMys6FJuWs7OzgoLC5Orq+s1LYfwAQCoEi4Hj7p166pGjRo8iNJilx8CmpaWptDQ0Gva/4QPAEClV1BQYAse/v7+FV3OTatOnTo6ceKE8vPzVb169TIvhwGnAIBK7/IYjxo1alRwJTe3y5dbCgoKrmk5hA8AQJXBpZaKVV77n/ABAAAs5XD4+Pnnn/Xoo4/K399fHh4eatGihb755hvbdGOMxo8fr6CgIHl4eCg6OloHDhwo16IBALgRde7cWSNHjqyw9Q8ePFi9e/e+7utxaMDpr7/+qvbt2+vee+/VV199pTp16ujAgQOqVauWrc+UKVP01ltvacGCBQoLC9OLL76omJgY/fDDD3J3dy/3DQAA3Nymr/zR0vWNuu9Wh/oPHjxYCxYs0PDhwzV37ly7afHx8Zo9e7ZiY2OVnJyspUuXXtNAzqrCofDx2muvKSQkRElJSba2sLAw28/GGM2YMUP/+Mc/1KtXL0nSe++9p4CAAH366ad65JFHyqlsAACqjpCQEC1atEjTp0+Xh4eHJOnSpUtauHChQkNDbf38/PyueV15eXmVPsA4dNnl888/V9u2bfXQQw+pbt26at26td555x3b9NTUVKWnpys6OtrW5uPjo8jISKWkpJRf1QAAVCF33HGHQkJCtHTpUlvb0qVLFRoaqtatW9va/njZJS0tTd27d5eHh4fCwsK0cOFCNWjQQDNmzLD1cXJy0pw5c/SnP/1Jnp6eevnll1VQUKC4uDiFhYXJw8ND4eHhevPNN+1qKigo0OjRo+Xr6yt/f3+NGzdOxpjrtg9+z6Hw8dNPP2nOnDlq0qSJVqxYoSeffFJ/+9vftGDBAkm/PQBGkgICAuzmCwgIsE37o5ycHGVnZ9u9AAC40QwdOtTuysG7776rIUOGlDrPoEGDdOLECa1bt06ffPKJ5s2bp5MnTxbpN3HiRPXp00e7d+/W0KFDVVhYqHr16mnJkiX64YcfNH78eL3wwgtavHixbZ433nhDycnJevfdd7Vx40adOXNGy5YtK78NLoVDl10KCwvVtm1bvfLKK5Kk1q1b6/vvv9fcuXMVGxtbpgISExM1adKkMs0L3PTWJlZ0BY67N6GiKwAqxKOPPqqEhAQdOXJEkrRp0yYtWrRI69atK7b/vn37tGrVKm3btk1t27aVJP3rX/9SkyZNivT9y1/+UiTI/P6zNSwsTCkpKVq8eLH69+8vSZoxY4YSEhLUt29fSdLcuXO1YsWKa97Oq+HQmY+goCA1a9bMrq1p06Y6evSoJCkwMFCSlJGRYdcnIyPDNu2PEhISlJWVZXsdO3bMkZIAAKgS6tSpo+7duys5OVlJSUnq3r27ateuXWL//fv3y8XFRXfccYetrXHjxnY3eVx2OZz83qxZs9SmTRvVqVNHNWvW1Lx582yf11lZWUpLS1NkZKStv4uLS7HLuR4cCh/t27fX/v377dp+/PFH1a9fX9JvySowMFCrV6+2Tc/OztbWrVsVFRVV7DLd3Nzk7e1t9wIA4EY0dOhQJScna8GCBRo6dGi5LdfT09Pu/aJFizR27FjFxcXpv//9r3bt2qUhQ4YoNze33NZ5LRwKH6NGjdKWLVv0yiuv6ODBg1q4cKHmzZun+Ph4Sb8Nehk5cqT++c9/6vPPP9fu3bs1aNAgBQcHW3LfMAAAldkDDzyg3Nxc5eXlKSYmptS+4eHhys/P186dO21tBw8e1K+//nrF9WzatEnt2rXTU089pdatW6tx48Y6dOiQbbqPj4+CgoK0detWW1t+fr62b99ehq1ynENjPu68804tW7ZMCQkJmjx5ssLCwjRjxgwNHDjQ1mfcuHE6f/68hg0bpszMTN1zzz1avnw5z/gAANz0qlWrpr1799p+Lk1ERISio6M1bNgwzZkzR9WrV9eYMWPk4eFxxcecN2nSRO+9955WrFihsLAwvf/++9q2bZvd4zFGjBihV199VU2aNFFERISmTZumzMzMa97Gq+Hwt9r26NFDPXr0KHG6k5OTJk+erMmTJ19TYQAA3IgcGV7w3nvvKS4uTh07dlRgYKASExO1Z8+eK/5BP3z4cO3cuVMPP/ywnJycNGDAAD311FP66quvbH3GjBmjtLQ0xcbGytnZWUOHDlWfPn2UlZVV5m27Wk7Gqpt6r1J2drZ8fHyUlZXF+A/gSrjbBTeJS5cuKTU1VWFhYTf1mfTjx48rJCREq1atUteuXS1ff2nHwZHPb4fPfAAAAGusWbNG586dU4sWLZSWlqZx48apQYMG6tixY0WXdk0IHwAAVFJ5eXl64YUX9NNPP8nLy0vt2rXThx9+WOkfn34lhA8AACqpmJiYK94VUxU5dKstAADAtSJ8AAAASxE+AACApQgfAADAUoQPAABgKcIHAACwFOEDAIBKonPnzho5cmSFrX/w4MGWfBEsz/kAAFRtVn/NgINfETB48GAtWLBAw4cP19y5c+2mxcfHa/bs2YqNjVVycrKWLl1a5R8gdjU48wEAwHUWEhKiRYsW6eLFi7a2S5cuaeHChQoNDbW1+fn5ycvL65rWlZeXd03zW4HwAQDAdXbHHXcoJCRES5cutbUtXbpUoaGhat26ta3tj5dd0tLS1L17d3l4eCgsLEwLFy5UgwYNNGPGDFsfJycnzZkzR3/605/k6empl19+WQUFBYqLi1NYWJg8PDwUHh6uN998066mgoICjR49Wr6+vvL399e4ceNk1XfNEj4AALDA0KFDlZSUZHv/7rvvasiQIaXOM2jQIJ04cULr1q3TJ598onnz5unkyZNF+k2cOFF9+vTR7t27NXToUBUWFqpevXpasmSJfvjhB40fP14vvPCCFi9ebJvnjTfeUHJyst59911t3LhRZ86c0bJly8pvg0vBmA8AACzw6KOPKiEhQUeOHJEkbdq0SYsWLdK6deuK7b9v3z6tWrVK27ZtU9u2bSVJ//rXv9SkSZMiff/yl78UCTKTJk2y/RwWFqaUlBQtXrxY/fv3lyTNmDFDCQkJ6tu3ryRp7ty5WrFixTVv59UgfAAAYIE6deqoe/fuSk5OljFG3bt3V+3atUvsv3//frm4uOiOO+6wtTVu3Fi1atUq0vdyOPm9WbNm6d1339XRo0d18eJF5ebmqlWrVpKkrKwspaWlKTIy0tbfxcVFbdu2teTSC+EDAACLDB06VE8//bSk38JBefH09LR7v2jRIo0dO1ZvvPGGoqKi5OXlpddff11bt24tt3VeC8Z8AABgkQceeEC5ubnKy8tTTExMqX3Dw8OVn5+vnTt32toOHjyoX3/99Yrr2bRpk9q1a6ennnpKrVu3VuPGjXXo0CHbdB8fHwUFBdmFkfz8fG3fvr0MW+U4znwAAGCRatWqae/evbafSxMREaHo6GgNGzZMc+bMUfXq1TVmzBh5eHjIycmp1HmbNGmi9957TytWrFBYWJjef/99bdu2TWFhYbY+I0aM0KuvvqomTZooIiJC06ZNU2Zm5jVv49XgzAcAABby9vaWt7f3VfV97733FBAQoI4dO6pPnz56/PHH5eXlJXd391LnGz58uPr27auHH35YkZGROn36tJ566im7PmPGjNFjjz2m2NhY26WZPn36lHm7HOFkrLqp9yplZ2fLx8dHWVlZV31wgJuW1U92LA8OPh0SkH57IFdqaqrCwsKu+MF7Izt+/LhCQkK0atUqde3a1fL1l3YcHPn85rILAACV1Jo1a3Tu3Dm1aNFCaWlpGjdunBo0aKCOHTtWdGnXhPABVGEpP52u6BIcFnVvRVcAVB15eXl64YUX9NNPP8nLy0vt2rXThx9+WOW//4XwAQBAJRUTE3PFu2KqIgacAgAASxE+AABVRiW7R+KmU177n/ABAKj0Lo9xuHDhQgVXcnPLzc2VdOVnlFwJYz4AAJVetWrV5Ovra/tG1xo1alzxQVsoX4WFhfrll19Uo0YNubhcW3wgfAAAqoTAwEBJKvYr5WENZ2dnhYaGXnPwI3wAAKoEJycnBQUFqW7dusrLy6vocm5Krq6ucna+9hEbhA8AQJVSrVq1ax5zgIrFgFMAAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUg6Fj4kTJ8rJycnuFRERYZt+6dIlxcfHy9/fXzVr1lS/fv2UkZFR7kUDAICqy+EzH7fddpvS0tJsr40bN9qmjRo1Sl988YWWLFmi9evX68SJE+rbt2+5FgwAAKo2F4dncHFRYGBgkfasrCzNnz9fCxcuVJcuXSRJSUlJatq0qbZs2aK777772qsFAABVnsNnPg4cOKDg4GA1bNhQAwcO1NGjRyVJ27dvV15enqKjo219IyIiFBoaqpSUlBKXl5OTo+zsbLsXAAC4cTkUPiIjI5WcnKzly5drzpw5Sk1NVYcOHXT27Fmlp6fL1dVVvr6+dvMEBAQoPT29xGUmJibKx8fH9goJCSnThgAAgKrBocsu3bp1s/3csmVLRUZGqn79+lq8eLE8PDzKVEBCQoJGjx5te5+dnU0AAQDgBnZNt9r6+vrq1ltv1cGDBxUYGKjc3FxlZmba9cnIyCh2jMhlbm5u8vb2tnsBAIAb1zWFj3PnzunQoUMKCgpSmzZtVL16da1evdo2ff/+/Tp69KiioqKuuVAAAHBjcOiyy9ixY9WzZ0/Vr19fJ06c0IQJE1StWjUNGDBAPj4+iouL0+jRo+Xn5ydvb28988wzioqK4k4XAABg41D4OH78uAYMGKDTp0+rTp06uueee7RlyxbVqVNHkjR9+nQ5OzurX79+ysnJUUxMjGbPnn1dCgcAAFWTkzHGVHQRv5ednS0fHx9lZWUx/gO4gpT5Yyu6BIdFxU2t6BIAXAeOfH7z3S4AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUi4VXQCAm0vK/LEVXYLDouKmVnQJwA2FMx8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUoQPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACWInwAAABLET4AAIClCB8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJa6pvDx6quvysnJSSNHjrS1Xbp0SfHx8fL391fNmjXVr18/ZWRkXGudAADgBlHm8LFt2za9/fbbatmypV37qFGj9MUXX2jJkiVav369Tpw4ob59+15zoQAA4MZQpvBx7tw5DRw4UO+8845q1apla8/KytL8+fM1bdo0denSRW3atFFSUpI2b96sLVu2lFvRAACg6ipT+IiPj1f37t0VHR1t1759+3bl5eXZtUdERCg0NFQpKSnFLisnJ0fZ2dl2LwAAcONycXSGRYsWaceOHdq2bVuRaenp6XJ1dZWvr69de0BAgNLT04tdXmJioiZNmuRoGQAAoIpy6MzHsWPHNGLECH344Ydyd3cvlwISEhKUlZVlex07dqxclgsAAConh8LH9u3bdfLkSd1xxx1ycXGRi4uL1q9fr7feeksuLi4KCAhQbm6uMjMz7ebLyMhQYGBgsct0c3OTt7e33QsAANy4HLrs0rVrV+3evduubciQIYqIiNBzzz2nkJAQVa9eXatXr1a/fv0kSfv379fRo0cVFRVVflUDAIAqy6Hw4eXlpebNm9u1eXp6yt/f39YeFxen0aNHy8/PT97e3nrmmWcUFRWlu+++u/yqBgAAVZbDA06vZPr06XJ2dla/fv2Uk5OjmJgYzZ49u7xXAwAAqqhrDh/r1q2ze+/u7q5Zs2Zp1qxZ17poAABwA+K7XQAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUoQPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACWInwAAABLET4AAIClCB8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUoQPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACWInwAAABLET4AAIClCB8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJZyKHzMmTNHLVu2lLe3t7y9vRUVFaWvvvrKNv3SpUuKj4+Xv7+/atasqX79+ikjI6PciwYAAFWXQ+GjXr16evXVV7V9+3Z988036tKli3r16qU9e/ZIkkaNGqUvvvhCS5Ys0fr163XixAn17dv3uhQOAACqJhdHOvfs2dPu/csvv6w5c+Zoy5YtqlevnubPn6+FCxeqS5cukqSkpCQ1bdpUW7Zs0d13311+VQMAgCqrzGM+CgoKtGjRIp0/f15RUVHavn278vLyFB0dbesTERGh0NBQpaSklEuxAACg6nPozIck7d69W1FRUbp06ZJq1qypZcuWqVmzZtq1a5dcXV3l6+tr1z8gIEDp6eklLi8nJ0c5OTm299nZ2Y6WBAAAqhCHw0d4eLh27dqlrKwsffzxx4qNjdX69evLXEBiYqImTZpU5vmBcrM2saIrAICbgsOXXVxdXdW4cWO1adNGiYmJuv322/Xmm28qMDBQubm5yszMtOufkZGhwMDAEpeXkJCgrKws2+vYsWMObwQAAKg6rvk5H4WFhcrJyVGbNm1UvXp1rV692jZt//79Onr0qKKiokqc383NzXbr7uUXAAC4cTl02SUhIUHdunVTaGiozp49q4ULF2rdunVasWKFfHx8FBcXp9GjR8vPz0/e3t565plnFBUVxZ0uAADAxqHwcfLkSQ0aNEhpaWny8fFRy5YttWLFCt13332SpOnTp8vZ2Vn9+vVTTk6OYmJiNHv27OtSOAAAqJqcjDGmoov4vezsbPn4+CgrK4tLMLBWFRxwmvLT6You4aYQFTe1oksAKj1HPr/5bhcAAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUoQPAABgKZeKLgAAKruU+WMruoQyiYqbWtElAMXizAcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUoQPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACWInwAAABLET4AAIClCB8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKUcCh+JiYm688475eXlpbp166p3797av3+/XZ9Lly4pPj5e/v7+qlmzpvr166eMjIxyLRoAAFRdDoWP9evXKz4+Xlu2bNHKlSuVl5en+++/X+fPn7f1GTVqlL744gstWbJE69ev14kTJ9S3b99yLxwAAFRNLo50Xr58ud375ORk1a1bV9u3b1fHjh2VlZWl+fPna+HCherSpYskKSkpSU2bNtWWLVt09913l1/lAACgSrqmMR9ZWVmSJD8/P0nS9u3blZeXp+joaFufiIgIhYaGKiUlpdhl5OTkKDs72+4FAABuXA6d+fi9wsJCjRw5Uu3bt1fz5s0lSenp6XJ1dZWvr69d34CAAKWnpxe7nMTERE2aNKmsZQDlJuWn0xVdAgDcFMp85iM+Pl7ff/+9Fi1adE0FJCQkKCsry/Y6duzYNS0PAABUbmU68/H000/r3//+tzZs2KB69erZ2gMDA5Wbm6vMzEy7sx8ZGRkKDAwsdllubm5yc3MrSxkAAKAKcujMhzFGTz/9tJYtW6Y1a9YoLCzMbnqbNm1UvXp1rV692ta2f/9+HT16VFFRUeVTMQAAqNIcOvMRHx+vhQsX6rPPPpOXl5dtHIePj488PDzk4+OjuLg4jR49Wn5+fvL29tYzzzyjqKgo7nQBAACSHAwfc+bMkSR17tzZrj0pKUmDBw+WJE2fPl3Ozs7q16+fcnJyFBMTo9mzZ5dLsQAAoOpzKHwYY67Yx93dXbNmzdKsWbPKXBQAALhx8d0uAADAUoQPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACWInwAAABLET4AAIClCB8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUoQPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACWInwAAABLET4AAIClCB8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEsRPgAAgKVcKroAAMD1kTJ/bEWX4LCouKkVXQIswJkPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACW4m4XXB9rEyu6AgBAJcWZDwAAYCmHw8eGDRvUs2dPBQcHy8nJSZ9++qnddGOMxo8fr6CgIHl4eCg6OloHDhwor3oBAEAV53D4OH/+vG6//XbNmjWr2OlTpkzRW2+9pblz52rr1q3y9PRUTEyMLl26dM3FAgCAqs/hMR/dunVTt27dip1mjNGMGTP0j3/8Q7169ZIkvffeewoICNCnn36qRx555NqqBQAAVV65jvlITU1Venq6oqOjbW0+Pj6KjIxUSkpKsfPk5OQoOzvb7gUAAG5c5Ro+0tPTJUkBAQF27QEBAbZpf5SYmCgfHx/bKyQkpDxLAgAAlUyF3+2SkJCgrKws2+vYsWMVXRIAALiOyjV8BAYGSpIyMjLs2jMyMmzT/sjNzU3e3t52LwAAcOMq1/ARFhamwMBArV692taWnZ2trVu3KioqqjxXBQAAqiiH73Y5d+6cDh48aHufmpqqXbt2yc/PT6GhoRo5cqT++c9/qkmTJgoLC9OLL76o4OBg9e7duzzrBgAAVZTD4eObb77Rvffea3s/evRoSVJsbKySk5M1btw4nT9/XsOGDVNmZqbuueceLV++XO7u7uVXNQAAqLKcjDGmoov4vezsbPn4+CgrK4vxH1VZFfxul5SfTld0CcBNLypuakWXgDJy5PO7wu92AQAANxfCBwAAsBThAwAAWIrwAQAALEX4AAAAlnL4VlvganDnCACgJJz5AAAAliJ8AAAASxE+AACApQgfAADAUoQPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACWInwAAABLET4AAIClCB8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACzlUtEF4CqsTazoCgAAKDec+QAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnudgEAVBop88dWdAkOi4qbWtElVDmc+QAAAJYifAAAAEsRPgAAgKUIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8AAAASxE+AACApQgfAADAUoQPAABgKcIHAACwFOEDAABYivABAAAsRfgAAACWcqnoAiy3NrGiK3BYyk+nK7oEAADKzXU78zFr1iw1aNBA7u7uioyM1Ndff329VgUAAKqQ6xI+PvroI40ePVoTJkzQjh07dPvttysmJkYnT568HqsDAABVyHUJH9OmTdPjjz+uIUOGqFmzZpo7d65q1Kihd99993qsDgAAVCHlHj5yc3O1fft2RUdH/99KnJ0VHR2tlJSU8l4dAACoYsp9wOmpU6dUUFCggIAAu/aAgADt27evSP+cnBzl5OTY3mdlZUmSsrOzy7u035y/dH2Wex2dv5hz5U4AgApx3T6vqpjL+8EYc8W+FX63S2JioiZNmlSkPSQkpAKqAQDAQc/8fxVdQaVy9uxZ+fj4lNqn3MNH7dq1Va1aNWVkZNi1Z2RkKDAwsEj/hIQEjR492va+sLBQZ86ckb+/v5ycnMq1tuzsbIWEhOjYsWPy9vYu12VXBjf69kk3/jayfVXfjb6NbF/Vd7220Rijs2fPKjg4+Ip9yz18uLq6qk2bNlq9erV69+4t6bdAsXr1aj399NNF+ru5ucnNzc2uzdfXt7zLsuPt7X3D/qOSbvztk278bWT7qr4bfRvZvqrvemzjlc54XHZdLruMHj1asbGxatu2re666y7NmDFD58+f15AhQ67H6gAAQBVyXcLHww8/rF9++UXjx49Xenq6WrVqpeXLlxcZhAoAAG4+123A6dNPP13sZZaK5ObmpgkTJhS5zHOjuNG3T7rxt5Htq/pu9G1k+6q+yrCNTuZq7okBAAAoJ3yrLQAAsBThAwAAWIrwAQAALEX4AAAAlrqhwsfLL7+sdu3aqUaNGiU+qOzo0aPq3r27atSoobp16+rZZ59Vfn5+qcs9c+aMBg4cKG9vb/n6+iouLk7nzp27DlvgmHXr1snJyanY17Zt20qcr3PnzkX6P/HEExZWfvUaNGhQpNZXX3211HkuXbqk+Ph4+fv7q2bNmurXr1+RJ+5WFocPH1ZcXJzCwsLk4eGhRo0aacKECcrNzS11vsp8DGfNmqUGDRrI3d1dkZGR+vrrr0vtv2TJEkVERMjd3V0tWrTQl19+aVGljktMTNSdd94pLy8v1a1bV71799b+/ftLnSc5ObnIsXJ3d7eoYsdMnDixSK0RERGlzlOVjp9U/P8pTk5Oio+PL7Z/ZT9+GzZsUM+ePRUcHCwnJyd9+umndtONMRo/fryCgoLk4eGh6OhoHThw4IrLdfT32FE3VPjIzc3VQw89pCeffLLY6QUFBerevbtyc3O1efNmLViwQMnJyRo/fnypyx04cKD27NmjlStX6t///rc2bNigYcOGXY9NcEi7du2UlpZm9/rrX/+qsLAwtW3bttR5H3/8cbv5pkyZYlHVjps8ebJdrc8880yp/UeNGqUvvvhCS5Ys0fr163XixAn17dvXomods2/fPhUWFurtt9/Wnj17NH36dM2dO1cvvPDCFeetjMfwo48+0ujRozVhwgTt2LFDt99+u2JiYnTy5Mli+2/evFkDBgxQXFycdu7cqd69e6t37976/vvvLa786qxfv17x8fHasmWLVq5cqby8PN1///06f/58qfN5e3vbHasjR45YVLHjbrvtNrtaN27cWGLfqnb8JGnbtm1227dy5UpJ0kMPPVTiPJX5+J0/f1633367Zs2aVez0KVOm6K233tLcuXO1detWeXp6KiYmRpculfwlq47+HpeJuQElJSUZHx+fIu1ffvmlcXZ2Nunp6ba2OXPmGG9vb5OTk1Pssn744QcjyWzbts3W9tVXXxknJyfz888/l3vt1yI3N9fUqVPHTJ48udR+nTp1MiNGjLCmqGtUv359M3369Kvun5mZaapXr26WLFlia9u7d6+RZFJSUq5DheVvypQpJiwsrNQ+lfUY3nXXXSY+Pt72vqCgwAQHB5vExMRi+/fv3990797dri0yMtIMHz78utZZXk6ePGkkmfXr15fYp6T/jyqjCRMmmNtvv/2q+1f142eMMSNGjDCNGjUyhYWFxU6vSsdPklm2bJntfWFhoQkMDDSvv/66rS0zM9O4ubmZ//f//l+Jy3H097gsbqgzH1eSkpKiFi1a2D1pNSYmRtnZ2dqzZ0+J8/j6+tqdSYiOjpazs7O2bt163Wt2xOeff67Tp09f1WPsP/zwQ9WuXVvNmzdXQkKCLly4YEGFZfPqq6/K399frVu31uuvv17qZbLt27crLy9P0dHRtraIiAiFhoYqJSXFinKvWVZWlvz8/K7Yr7Idw9zcXG3fvt1u3zs7Oys6OrrEfZ+SkmLXX/rtd7IqHStJVzxe586dU/369RUSEqJevXqV+P9NZXDgwAEFBwerYcOGGjhwoI4ePVpi36p+/HJzc/XBBx9o6NChpX6RaVU6fr+Xmpqq9PR0u2Pk4+OjyMjIEo9RWX6Py+K6PeG0MkpPTy/yiPfL79PT00ucp27dunZtLi4u8vPzK3GeijJ//nzFxMSoXr16pfb7y1/+ovr16ys4OFjfffednnvuOe3fv19Lly61qNKr97e//U133HGH/Pz8tHnzZiUkJCgtLU3Tpk0rtn96erpcXV2LjPkJCAiodMerOAcPHtTMmTM1derUUvtVxmN46tQpFRQUFPs7tm/fvmLnKel3siocq8LCQo0cOVLt27dX8+bNS+wXHh6ud999Vy1btlRWVpamTp2qdu3aac+ePVf8XbVaZGSkkpOTFR4errS0NE2aNEkdOnTQ999/Ly8vryL9q/Lxk6RPP/1UmZmZGjx4cIl9qtLx+6PLx8GRY1SW3+OyqPTh4/nnn9drr71Wap+9e/decVBUVVKWbT5+/LhWrFihxYsXX3H5vx+v0qJFCwUFBalr1646dOiQGjVqVPbCr5Ij2zd69GhbW8uWLeXq6qrhw4crMTGxUj/+uCzH8Oeff9YDDzyghx56SI8//nip81b0MYQUHx+v77//vtQxEZIUFRWlqKgo2/t27dqpadOmevvtt/XSSy9d7zId0q1bN9vPLVu2VGRkpOrXr6/FixcrLi6uAiu7PubPn69u3bqV+hXwVen4VSWVPnyMGTOm1FQqSQ0bNryqZQUGBhYZsXv5LojAwMAS5/njIJv8/HydOXOmxHmuVVm2OSkpSf7+/vrTn/7k8PoiIyMl/fZXtxUfXNdyTCMjI5Wfn6/Dhw8rPDy8yPTAwEDl5uYqMzPT7uxHRkbGdTtexXF0G0+cOKF7771X7dq107x58xxen9XHsDi1a9dWtWrVitxZVNq+DwwMdKh/ZfH000/bBp87+tdv9erV1bp1ax08ePA6VVd+fH19deutt5ZYa1U9fpJ05MgRrVq1yuGzhVXp+F0+DhkZGQoKCrK1Z2RkqFWrVsXOU5bf4zIpt9EjlciVBpxmZGTY2t5++23j7e1tLl26VOyyLg84/eabb2xtK1asqFQDTgsLC01YWJgZM2ZMmebfuHGjkWS+/fbbcq6s/H3wwQfG2dnZnDlzptjplwecfvzxx7a2ffv2VeoBp8ePHzdNmjQxjzzyiMnPzy/TMirLMbzrrrvM008/bXtfUFBgbrnlllIHnPbo0cOuLSoqqtIOWCwsLDTx8fEmODjY/Pjjj2VaRn5+vgkPDzejRo0q5+rK39mzZ02tWrXMm2++Wez0qnb8fm/ChAkmMDDQ5OXlOTRfZT5+KmHA6dSpU21tWVlZVzXg1JHf4zLVWm5LqgSOHDlidu7caSZNmmRq1qxpdu7caXbu3GnOnj1rjPntH03z5s3N/fffb3bt2mWWL19u6tSpYxISEmzL2Lp1qwkPDzfHjx+3tT3wwAOmdevWZuvWrWbjxo2mSZMmZsCAAZZvX0lWrVplJJm9e/cWmXb8+HETHh5utm7daowx5uDBg2by5Mnmm2++Mampqeazzz4zDRs2NB07drS67CvavHmzmT59utm1a5c5dOiQ+eCDD0ydOnXMoEGDbH3+uH3GGPPEE0+Y0NBQs2bNGvPNN9+YqKgoExUVVRGbcEXHjx83jRs3Nl27djXHjx83aWlpttfv+1SVY7ho0SLj5uZmkpOTzQ8//GCGDRtmfH19bXeYPfbYY+b555+39d+0aZNxcXExU6dONXv37jUTJkww1atXN7t3766oTSjVk08+aXx8fMy6devsjtWFCxdsff64jZMmTTIrVqwwhw4dMtu3bzePPPKIcXd3N3v27KmITSjVmDFjzLp160xqaqrZtGmTiY6ONrVr1zYnT540xlT943dZQUGBCQ0NNc8991yRaVXt+J09e9b2WSfJTJs2zezcudMcOXLEGGPMq6++anx9fc1nn31mvvvuO9OrVy8TFhZmLl68aFtGly5dzMyZM23vr/R7XB5uqPARGxtrJBV5rV271tbn8OHDplu3bsbDw8PUrl3bjBkzxi75rl271kgyqamptrbTp0+bAQMGmJo1axpvb28zZMgQW6CpDAYMGGDatWtX7LTU1FS7fXD06FHTsWNH4+fnZ9zc3Ezjxo3Ns88+a7Kysiys+Ops377dREZGGh8fH+Pu7m6aNm1qXnnlFbuzVH/cPmOMuXjxonnqqadMrVq1TI0aNUyfPn3sPswrk6SkpGL/zf7+pGRVO4YzZ840oaGhxtXV1dx1111my5YttmmdOnUysbGxdv0XL15sbr31VuPq6mpuu+0285///Mfiiq9eSccqKSnJ1ueP2zhy5Ejb/ggICDAPPvig2bFjh/XFX4WHH37YBAUFGVdXV3PLLbeYhx9+2Bw8eNA2vaofv8tWrFhhJJn9+/cXmVbVjt/lz6w/vi5vQ2FhoXnxxRdNQECAcXNzM127di2y3fXr1zcTJkywayvt97g8OBljTPldxAEAACjdTfWcDwAAUPEIHwAAwFKEDwAAYCnCBwAAsBThAwAAWIrwAQAALEX4AAAAliJ8ALDE4MGD1bt374ouA0AlwEPGAFgiKytLxhi7L/wDcHMifAAAAEtx2QVAufr444/VokULeXh4yN/fX9HR0Tp//nyRyy5nz57VwIED5enpqaCgIE2fPl2dO3fWyJEjK6x2ANYgfAAoN2lpaRowYICGDh2qvXv3at26derbt6+KO8E6evRobdq0SZ9//rlWrlyp//3vf9qxY0cFVA3Aai4VXQCAG0daWpry8/PVt29f1a9fX5LUokWLIv3Onj2rBQsWaOHCherataskKSkpScHBwZbWC6BicOYDQLm5/fbb1bVrV7Vo0UIPPfSQ3nnnHf36669F+v3000/Ky8vTXXfdZWvz8fFReHi4leUCqCCEDwDlplq1alq5cqW++uorNWvWTDNnzlR4eLhSU1MrujQAlQjhA0C5cnJyUvv27TVp0iTt3LlTrq6uWrZsmV2fhg0bqnr16tq2bZutLSsrSz/++KPV5QKoAIz5AFButm7dqtWrV+v+++9X3bp1tXXrVv3yyy9q2rSpvvvuO1s/Ly8vxcbG6tlnn5Wfn5/q1q2rCRMmyNnZWU5OThW4BQCswJkPAOXG29tbGzZs0IMPPqhbb71V//jHP/TGG2+oW7duRfpOmzZNUVFR6tGjh6Kjo9W+fXs1bdpU7u7uFVA5ACvxkDEAlcL58+d1yy236I033lBcXFxFlwPgOuKyC4AKsXPnTu3bt0933XWXsrKyNHnyZElSr169KrgyANcb4QNAhZk6dar2798vV1dXtWnTRv/73/9Uu3btii4LwHXGZRcAAGApBpwCAABLET4AAIClCB8AAMBShA8AAGApwgcAALAU4QMAAFiK8AEAACxF+AAAAJYifAAAAEv9/7elv1UxvFljAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(f\"{np.sum(np.isnan(sigs_constrained))} constrained fits failed\")\n",
    "plt.hist(sigs_migrad, alpha=0.5, bins=10, range=(-10, 10), label=\"Migrad\")\n",
    "plt.hist(sigs_constrained, alpha=0.5, bins=10, range=(-10, 10), label=m.fmin.algorithm)\n",
    "plt.xlabel(\"sig\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "881008cc",
   "metadata": {},
   "source": [
    "There are no failures this time. \n",
    "\n",
    "For sig > 0, the distributions are identical in this example, as theoretically expected. In practice, there can be small bin migration effects due to finite precision of numerical algorithms. These are not of concern.\n",
    "\n",
    "Important are the differences for sig < 0, where Migrad did not converge in a few cases and where therefore samples are missing. Those missing samples are recovered in the distribution produced by the constrained fit.\n",
    "\n",
    "This demonstrates that it is important to not discard failed fits, as this will in general distort the distribution of the test statistic."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f04e7c57",
   "metadata": {},
   "source": [
    "## Bonus: unconstrained SciPy fit\n",
    "\n",
    "The issues we describe here are of a principal mathematical nature. We should not expect that an unconstrained minimiser from SciPy does better than Migrad, but let's test this assumption. The minimiser that SciPy uses when only box constraints are used is the L-BFGS-B method which is roughly comparable to Migrad. Let us see how well this algorithm does on the same toy samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a60a75f7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17 failed\n"
     ]
    }
   ],
   "source": [
    "@joblib.delayed\n",
    "def compute(itry):\n",
    "    rng = np.random.default_rng(itry)\n",
    "    x = rng.uniform(0, 1, size=35)\n",
    "    cost = ExtendedUnbinnedNLL(x, model)\n",
    "    m = Minuit(cost, b0=n, b1=n, b2=n, sig=0, mu=0.5, sigma=0.05)\n",
    "    m.limits[\"b0\", \"b1\", \"b2\"] = (0, None)\n",
    "    m.fixed[\"mu\", \"sigma\"] = True\n",
    "    m.scipy()\n",
    "    return m.values[\"sig\"] if m.valid else np.nan\n",
    "\n",
    "sigs_bfgs = joblib.Parallel(-1)(compute(i) for i in range(200))\n",
    "\n",
    "print(np.sum(np.isnan(sigs_bfgs)), \"failed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "589bf6b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAHHCAYAAAAf2DoOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABKB0lEQVR4nO3deXwN1/8/8NdNIvsmiyxksySx1xqhRYmiqCUfNLUk1tLQxtJqFKm0FVRRmtJqJVpUUYK2KNqgthJSW8TSEJGF0NyQNPv5/eFnvr3Ndm/czE3i9Xw85sE9c2bmPXcS92XmzFyFEEKAiIiISCZ6ui6AiIiIni0MH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EpBP79u3Dc889B2NjYygUCmRlZSEoKAju7u462ba63n//fSgUCpU2d3d3BAUFaa2+mzdvQqFQIDo6WmvrJKpJGD6ItODRo0cICwtDv379YGNjU+EHh0KhKHfq06dPpdv67zJmZmZo0aIFPvzwQ+Tm5qr0DQoKKndb+/btU+mbnZ2Njz76CB07doSVlRWMjIzg5uaGkSNH4qeffipVx82bNzFu3Dg0adIExsbGcHR0RPfu3REWFlbpPty/fx8jRoyAiYkJIiMj8e2338LMzKxUv9zcXLz//vuIjY2tdJ3qUnfbRFR9DHRdAFFdkJmZifDwcLi6uqJt27YVflh+++23pdrOnDmDTz/9FC+99JJa2+vTpw/Gjh0L4HHwOXr0KObPn48///wT27ZtU+lrZGSEr776qtQ62rZtK/39+vXr6Nu3L27duoWhQ4di7NixMDc3x+3bt/Hzzz9j4MCB+OabbzBmzBipf6dOnWBiYoLx48fD3d0daWlpOHv2LJYsWYKFCxdWWP/p06fx8OFDfPDBB/Dz85Pa161bh5KSEul1bm6utK6ePXuq9d5Uprxtq2vevHl49913tVIL0bOK4YNIC5ycnJCWlgZHR0ecOXMGnTp1Krfv6NGjS7XFxsZCoVAgICBAre15enqqrGfKlCkoKCjAjh07kJeXB2NjY2megYFBmdt8oqioCEOHDkVGRgYOHz6Mbt26qcwPCwvDL7/8guLiYqltxYoVePToEeLj4+Hm5qbS/+7du5XW/6SPtbW1Snu9evUqXfZplbdtdRkYGMDAgP90Ej0NXnYh0gIjIyM4OjpWadn8/Hz88MMP6NGjBxo1alTlGhwdHaFQKDT+YNy2bRsuXryI+fPnlwoeT7z00kvo37+/9PrGjRto1KhRqeABAA0aNKhwez179kRgYCAAoFOnTlAoFNJ4iX+P+bh58ybs7e0BAAsXLpQuF73//vsAgPT0dIwbNw6NGjWCkZERnJycMHjwYNy8ebNK2z569CiGDx8OV1dXGBkZwcXFBTNmzMA///yjso6yxnyUJSsrCyEhIXBxcYGRkRGaNm2KJUuWqJzZedIvKCgIVlZWsLa2RmBgoEZjUIhqI8Z3Ih37+eefkZWVhVGjRqm9TF5eHjIzMwEAOTk5OHbsGDZs2IDXXnutzPDxpO8T9erVg5WVFQBgz549AMo+I1MeNzc3HDx4EL/++it69eql9nIA8N5778HLywtffvklwsPD4eHhgSZNmpTqZ29vjzVr1mDq1KkYOnQohg0bBgBo06YNAMDf3x+XLl3C9OnT4e7ujrt37+LAgQNITk4ud9BqRdvetm0bcnNzMXXqVNja2uKPP/7A6tWrkZKSUupSVmVyc3PRo0cP3LlzB6+//jpcXV1x/PhxhIaGIi0tDStXrgQACCEwePBg/P7775gyZQqaN2+OnTt3SgGJqM4SRKRVp0+fFgBEVFSUWv39/f2FkZGR+Pvvv9XqD6DMaciQISIvL0+lb2BgYJl9e/ToIfVp166dsLa2LrWdR48eiXv37kmTUqmU5l28eFGYmJgIAOK5554Tb731loiJiRE5OTlq7UNUVJQAIE6fPl2qXjc3N+n1vXv3BAARFham0u/vv/8WAMTHH3+s1vbU2XZubm6pvhEREUKhUIhbt25JbWFhYeK//3S6ubmJwMBA6fUHH3wgzMzMxNWrV1X6vfvuu0JfX18kJycLIYSIiYkRAMTSpUulPkVFReKFF17Q6GeIqLbhZRciHcrOzsZPP/2El19+WaMxCIMHD8aBAwdw4MAB7Nq1C6Ghodi3bx9ee+01CCFU+hobG0t9n0yffPKJSg3m5ualtvHee+/B3t5eml577TVpXsuWLREfH4/Ro0fj5s2b+PTTTzFkyBA4ODhg3bp1mr8RGjIxMYGhoSFiY2Px999/a22dT+Tk5CAzMxNdu3aFEALnzp3TaF3btm3DCy+8gPr16yMzM1Oa/Pz8UFxcjCNHjgB4fNbLwMAAU6dOlZbV19fH9OnTtbJPRDUVL7sQ6dAPP/yAvLw8jS65AECjRo1U7tR45ZVXYGtri9mzZ+PHH3/EoEGDpHn6+voV3tVhYWGB+/fvl2p/4403MHDgQABlX5Lx9PTEt99+i+LiYly+fBk//vgjli5dismTJ8PDw6NKd5Koy8jICEuWLMGsWbPg4OCALl26YODAgRg7dmyVx94kJydjwYIF2L17d6lAo1QqNVrXtWvXcP78eWnMyn89GfR669YtODk5lQp/Xl5eGm2PqLZh+CDSoU2bNsHKykr6kH8avXv3BgAcOXJEJXxUxtvbG/Hx8bhz5w4aNmwotXt6esLT0xMAVO6e+S99fX20bt0arVu3hq+vL1588UVs2rSpWsMHAISEhGDQoEGIiYnB/v37MX/+fERERODXX39Fu3btNFpXcXEx+vTpgwcPHmDOnDnw9vaGmZkZ7ty5g6CgoFKDRCtTUlKCPn364J133ilz/pP3lehZxfBBpCNpaWn47bffEBQUBCMjo6deX1FREYDHz/3QxMCBA7FlyxZs2rSp3A9LdXXs2BHA433ThsruKmnSpAlmzZqFWbNm4dq1a3juuefwySefYOPGjRpt58KFC7h69So2bNggPT8FAA4cOFClups0aYJHjx5VGsDc3Nxw6NAhPHr0SOXsR2JiYpW2S1RbcMwHkY5s2bIFJSUlGl9yKc+Tu1b+/fAwdYwYMQItWrTABx98gJMnT5bZ57/jSI4ePYrCwsJS/X7++WcA2rtsYGpqCgClbj3Nzc1FXl6eSluTJk1gYWGB/Px8jbejr68PQHU/hRD49NNPNV4X8Pg9PXHiBPbv319qXlZWlhQUX375ZRQVFWHNmjXS/OLiYqxevbpK2yWqLXjmg0hLPvvsM2RlZSE1NRXA4zCQkpICAJg+fbp0a+sTmzZtgrOzc5We3Hn16lXpf/e5ubk4efIkNmzYgKZNm0pPIVVXvXr1sHPnTvTt2xfPP/88hg0bhhdeeEG67LB7924kJydjwIAB0jJLlixBXFwchg0bJt36evbsWXzzzTewsbFBSEiIxvtUFhMTE7Ro0QLff/89PD09YWNjg1atWqGoqAi9e/eWgpOBgQF27tyJjIwMvPrqqxpvx9vbG02aNMHs2bNx584dWFpa4ocffqjyYNa3334bu3fvxsCBAxEUFIQOHTogJycHFy5cwPbt23Hz5k3Y2dlh0KBB6NatG959913cvHkTLVq0wI4dOzQeY0JU6+j0XhuiOsTNza3c22CTkpJU+l65ckUAEDNnztR4O/9dt76+vmjUqJGYPHmyyMjIUOkbGBgozMzM1FpvVlaWCA8PF+3atRPm5ubC0NBQuLi4iP/9739iz549Kn2PHTsmgoODRatWrYSVlZWoV6+ecHV1FUFBQeLGjRuVbkvdW22FEOL48eOiQ4cOwtDQULrtNjMzUwQHBwtvb29hZmYmrKyshI+Pj9i6dWuVt3358mXh5+cnzM3NhZ2dnZg0aZL4888/S93yqs6ttkII8fDhQxEaGiqaNm0qDA0NhZ2dnejatatYtmyZKCgokPrdv39fjBkzRlhaWgorKysxZswYce7cOd5qS3WaQoj/nE8lIiIiqkYc80FERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK4YPIiIiklWNe8hYSUkJUlNTYWFhUemjlYmIiKhmEELg4cOHcHZ2hp5exec2alz4SE1NhYuLi67LICIioiq4ffs2GjVqVGGfGhc+LCwsADwu3tLSUsfVEBERkTqys7Ph4uIifY5XpMaFjyeXWiwtLRk+iIiIahl1hkxwwCkRERHJiuGDiIiIZMXwQURERLKqcWM+iIiociUlJSgoKNB1GfSMMTQ0rPQ2WnUwfBAR1TIFBQVISkpCSUmJrkuhZ4yenh48PDxgaGj4VOth+CAiqkWEEEhLS4O+vj5cXFy08r9QInU8eQhoWloaXF1dn+pBoAwfRES1SFFREXJzc+Hs7AxTU1Ndl0PPGHt7e6SmpqKoqAj16tWr8noYmYmIapHi4mIAeOrT3kRV8eTn7snPYVUxfBAR1UL87ivSBW393DF8EBERkaw0Dh937tzB6NGjYWtrCxMTE7Ru3RpnzpyR5gshsGDBAjg5OcHExAR+fn64du2aVosmIqK6q2fPnggJCdHZ9oOCgjBkyBCdbf9ZoNGA07///hvdunXDiy++iL1798Le3h7Xrl1D/fr1pT5Lly7FqlWrsGHDBnh4eGD+/Pno27cvLl++DGNjY63vABERASsOXJV1ezP6eGrUPygoCBs2bMDrr7+OtWvXqswLDg7G559/jsDAQERHR2PHjh1PNZiRaj6NwseSJUvg4uKCqKgoqc3Dw0P6uxACK1euxLx58zB48GAAwDfffAMHBwfExMTg1Vdf1VLZRERU27i4uGDLli1YsWIFTExMAAB5eXnYvHkzXF1dpX42NjZPva3CwkIGmBpMo8suu3fvRseOHTF8+HA0aNAA7dq1w7p166T5SUlJSE9Ph5+fn9RmZWUFHx8fnDhxQntVExFRrdO+fXu4uLhgx44dUtuOHTvg6uqKdu3aSW3/veySlpaGAQMGwMTEBB4eHti8eTPc3d2xcuVKqY9CocCaNWvwyiuvwMzMDB999BGKi4sxYcIEeHh4wMTEBF5eXvj0009VaiouLsbMmTNhbW0NW1tbvPPOOxBCVNt7QI9pFD7++usvrFmzBs2aNcP+/fsxdepUvPnmm9iwYQMAID09HQDg4OCgspyDg4M077/y8/ORnZ2tMhERUd00fvx4lbPn69evx7hx4ypcZuzYsUhNTUVsbCx++OEHfPnll7h7926pfu+//z6GDh2KCxcuYPz48SgpKUGjRo2wbds2XL58GQsWLMDcuXOxdetWaZlPPvkE0dHRWL9+PX7//Xc8ePAAO3fu1N4OU5k0uuxSUlKCjh07YtGiRQCAdu3a4eLFi1i7di0CAwOrVEBERAQWLlxYpWWJnnm/Rei6As29GKrrCkiHRo8ejdDQUNy6dQsAcOzYMWzZsgWxsbFl9r9y5QoOHjyI06dPo2PHjgCAr776Cs2aNSvV97XXXisVZP79+eLh4YETJ05g69atGDFiBABg5cqVCA0NxbBhwwAAa9euxf79+596P6liGp35cHJyQosWLVTamjdvjuTkZACAo6MjACAjI0OlT0ZGhjTvv0JDQ6FUKqXp9u3bmpRERES1iL29PQYMGIDo6GhERUVhwIABsLOzK7d/YmIiDAwM0L59e6mtadOmKjc6PPEknPxbZGQkOnToAHt7e5ibm+PLL7+UPrOUSiXS0tLg4+Mj9TcwMChzPaRdGoWPbt26ITExUaXt6tWrcHNzA/A4VTo6OuLQoUPS/OzsbJw6dQq+vr5lrtPIyAiWlpYqExER1V3jx49HdHQ0NmzYgPHjx2ttvWZmZiqvt2zZgtmzZ2PChAn45ZdfEB8fj3HjxvHbgGsAjcLHjBkzcPLkSSxatAjXr1/H5s2b8eWXXyI4OBjA4wE/ISEh+PDDD7F7925cuHABY8eOhbOzM++ZJiIiAEC/fv1QUFCAwsJC9O3bt8K+Xl5eKCoqwrlz56S269ev4++//650O8eOHUPXrl3xxhtvoF27dmjatClu3LghzbeysoKTkxNOnToltRUVFSEuLq4Ke0Wa0GjMR6dOnbBz506EhoYiPDwcHh4eWLlyJUaNGiX1eeedd5CTk4PJkycjKysLzz//PPbt28dnfBAREQBAX18fCQkJ0t8r4u3tDT8/P0yePBlr1qxBvXr1MGvWLJiYmFT6qO9mzZrhm2++wf79++Hh4YFvv/0Wp0+fVnlExFtvvYXFixejWbNm8Pb2xvLly5GVlfXU+0gV0/hbbQcOHIiBAweWO1+hUCA8PBzh4eFPVRgRVS45IwuZyhxdl6ERuybJKs90oGeTJpfYv/nmG0yYMAHdu3eHo6MjIiIicOnSpUr/U/v666/j3LlzGDlyJBQKBQICAvDGG29g7969Up9Zs2YhLS0NgYGB0NPTw/jx4zF06FAolcoq7xtVTiFq2A3N2dnZsLKyglKp5PgPogokJyejuVdT5OYV6roUjZiamiIhIYEBpIry8vKQlJQEDw+PZ/aMckpKClxcXHDw4EH07t1b1+U8Uyr6+dPk81vjMx9EVDNkZmYiN68QG+eOQHNXe12Xo5aE5HsYvWgrMjMzGT5Ibb/++isePXqE1q1bIy0tDe+88w7c3d3RvXt3XZdGVcTwQVTLNXe1R3vPhroug6jaFBYWYu7cufjrr79gYWGBrl27YtOmTXx8ei3G8EFERDVa3759K70rhmoXjW61JSIiInpaDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxfBBREREsmL4ICKiGi86OhrW1tZaXae7uzsUCgUUCkWt/j6X6OhoaT9CQkJ0XY5a+JwPIqK64LcIebf3YqjGi9y7dw8LFizATz/9hIyMDNSvXx9t27bFggUL0K1btwqXHTlyJF5++WXpdXR0NMaNGwfg8XeKOTs7o0+fPliyZAkaNGigdk3h4eGYNGkSrKyspLZ169bhs88+w40bN2BgYAAPDw+MGDECoaGP9/n9999HTEwM4uPjy1xnUlIS3nvvPcTGxuLBgwews7NDhw4dsGTJEnh7e0v9fvzxR3z88cc4e/YsiouL0bJlSwQHByMoKEhlfTt37sSSJUuQkJCAkpISuLq6ok+fPli5cqX03vTr1w/Dhg1Te791jeGDiIhk4e/vj4KCAmzYsAGNGzdGRkYGDh06hPv371e6rImJCUxMTFTaLC0tkZiYiJKSEvz5558YN24cUlNTsX//frVrsrCwgKOjo/R6/fr1CAkJwapVq9CjRw/k5+fj/PnzuHjxolrrKywsRJ8+feDl5YUdO3bAyckJKSkp2Lt3r8rZldWrVyMkJARz5szBmjVrYGhoiF27dmHKlCm4ePEili1bBgA4dOgQRo4ciY8++givvPIKFAoFLl++jAMHDpR6bwwNDdXeb11j+CAiomqXlZWFo0ePIjY2Fj169AAAuLm5oXPnzip95syZg5iYGCiVSjRt2hSLFy/GwIEDER0djZCQEJUPcIVCIQUHZ2dnvPnmm5g/fz7++ecfDBgwAC1atMBnn30m9b937x4aNmyIvXv3lvuFdLt378aIESMwYcIEqa1ly5Zq7+elS5dw48YNHDp0CG5ubtJ+/vvMzu3btzFr1iyEhIRg0aJFUvusWbNgaGiIN998E8OHD4ePjw/27NmDbt264e2335b6eXp6YsiQIWrXVBNxzAcREVU7c3NzmJubIyYmBvn5+aXml5SUoH///jh27Bg2btyIy5cvY/HixdDX11d7GyYmJigpKUFRUREmTpyIzZs3q2xr48aNaNiwIXr16lXuOhwdHXHy5EncunVLsx38/+zt7aGnp4ft27ejuLi4zD7bt29HYWEhZs+eXWre66+/DnNzc3z33XdSPZcuXVL7zEttwfBBRETVzsDAANHR0diwYQOsra3RrVs3zJ07F+fPnwcAHDx4EH/88Qd27NiBPn36oHHjxhg4cCD69++v1vqvXbuGtWvXomPHjrCwsJDGP+zatUvqEx0djaCgICgUinLXExYWBmtra7i7u8PLywtBQUHYunUrSkpK1KqjYcOGWLVqFRYsWID69eujV69e+OCDD/DXX39Jfa5evQorKys4OTmVWt7Q0BCNGzfG1atXAQDTp09Hp06d0Lp1a7i7u+PVV1/F+vXrywxwtQnDBxERycLf3x+pqanYvXs3+vXrh9jYWLRv3x7R0dGIj49Ho0aN4Onpqfb6lEolzM3NYWpqCi8vLzg4OGDTpk0AAGNjY4wZMwbr168HAJw9exYXL14sNZjzv5ycnHDixAlcuHABb731FoqKihAYGIh+/fqpHUCCg4ORnp6OTZs2wdfXF9u2bUPLli1VxmlU5sn4DTMzM/z000+4fv065s2bB3Nzc8yaNQudO3dGbm6u2uuraRg+iIhINsbGxujTpw/mz5+P48ePIygoCGFhYaUGk6rDwsIC8fHxuHjxInJycnDkyBGV8DJx4kQcOHAAKSkpiIqKQq9evaRxGJVp1aoV3njjDWzcuBEHDhzAgQMHcPjwYY1qGzRoED766CP8+eefeOGFF/Dhhx8CAJo1awalUonU1NRSyxUUFODGjRulQliTJk0wceJEfPXVVzh79iwuX76M77//Xu16ahqGDyIi0pkWLVogJycHbdq0QUpKinS5QR16enpo2rQpGjduXGZ4ad26NTp27Ih169Zh8+bNGD9+fJVrBICcnJwqLa9QKODt7S0t/7///Q8GBgb45JNPSvVdu3YtcnNzMXbs2HLX5+7uDlNT0yrXUxPwbhciIqp29+/fx/DhwzF+/Hi0adMGFhYWOHPmDJYuXYrBgwejR48e6N69O/z9/bF8+XI0bdoUV65cgUKhQL9+/aq83YkTJ2LatGkwMzPD0KFDK+0/depUODs7o1evXmjUqBHS0tLw4Ycfwt7eHr6+vlK/f/75p9RzPiwsLPDw4UOEhYVhzJgxaNGiBQwNDXH48GGsX78ec+bMAQC4urpi6dKlmD17tnR5qF69eti1axfmzp2LDz/8EK1atQLw+Jkiubm5ePnll+Hm5oasrCysWrVKuqW3tmL4ICKiamdubg4fHx+sWLECN27cQGFhIVxcXDBp0iTMnTsXAPDDDz9g9uzZCAgIQE5OjnSr7dMICAhASEgIAgICYGxsXGl/Pz8/rF+/HmvWrMH9+/dhZ2cHX19fHDp0CLa2tlK/q1evol27dirL9u7dG1u2bIG7uzsWLlyImzdvQqFQSK9nzJgh9Z0xYwYaN26MTz75BJ9++ql0FuO7777Dq6++KvXr0aMHIiMjMXbsWOnBbO3atcMvv/wCLy+vp3pvdEkhhBC6LuLfsrOzYWVlBaVSCUtLS12XQ1RjnT17Fh06dEDc2mC092yo63LUcvbqHXSYEom4uDi0b99e1+XUSnl5eUhKSoKHh4daH6bPups3b6JJkyY4ffp0qZ85d3d3hISE1IhHkj948AC9e/eGpaUl9u7dC1NTU43X0bNnTzz33HPSk0+rQ0U/f5p8fnPMBxER1TmFhYVIT0/HvHnz0KVLl3LD7pw5c2Bubg6lUilzhapsbGxw8OBB9O7dGydOnNBo2U2bNsHc3BxHjx6tpuq0j5ddiGq583eUyDeoHf8DTryj23/g6dlx7NgxvPjii/D09MT27dvL7HP48GEUFhYCeDxeQ9dsbW2xYMECjZd75ZVX4OPjAwBa//K96sLwQUREdU7Pnj1R2agCdW+7reksLCxqRHjSBC+7EBERkawYPoiIiEhWDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxVttiYjqgOTkZGRmZsq2PTs7O7i6usq2vejoaISEhCArK0tr63R3d8etW7cAAH///XeteUZGdYmNjcWLL74IABg8eDBiYmKqbVsMH0REtVxycjKaN2+O3Nxc2bZpamqKhIQEjQLIvXv3sGDBAvz000/S95S0bdsWCxYsQLdu3SpcduTIkXj55Zel19HR0Rg3bhyAx98a6+zsjD59+mDJkiVo0KCB2jWFh4dj0qRJsLKyAvB/H8DqhpF/1wEAZmZm8PLywnvvvYdhw4ZJ7T179sThw4dLLV9YWAgDg8cfxdevX8eiRYtw8OBBZGRkwM7ODt7e3hg/fjxGjhwp9Tt8+DAWLlyI+Ph45OXloWHDhujatSvWrVsHQ0PDMuv8d9DS09ODg4MD+vfvj2XLlqF+/foAgK5duyItLQ1vvfUW8vPz1Xj3qo7hg4iolsvMzERubi42btyI5s2bV/v2EhISMHr0aGRmZmoUPvz9/VFQUIANGzagcePGyMjIwKFDh3D//v1KlzUxMYGJiYlKm6WlJRITE1FSUoI///wT48aNQ2pqKvbv3692TRYWFnB0dFS7f1me1AEADx8+RFRUFEaMGIFLly6pfPnbpEmTEB4errLsk0Dxxx9/wM/PDy1btkRkZCS8vb0BAGfOnEFkZCRatWqFtm3b4vLly+jXrx+mT5+OVatWwcTEBNeuXcMPP/yA4uLiCut8ErSKi4tx9epVTJ48GW+++Sa+/fZbAIChoSEcHR1hYmLC8EFEROpp3rx5jf3CvqysLBw9ehSxsbHo0aMHgMdPGO3cubNKnzlz5iAmJgZKpVL6VtuBAweWedlFoVBIwcHZ2Rlvvvkm5s+fj3/++QcDBgxAixYt8Nlnn0n97927h4YNG2Lv3r3o3bu31vbt33U4Ojriww8/xLJly3D+/HmV8GFqalpm0BFCICgoCJ6enjh27Bj09P5vOGazZs0QEBAgPa31l19+gaOjI5YuXSr1adKkCfr161dpnf8OWg0bNkRgYCC+++67qu30U+KAUyIiqnbm5uYwNzdHTExMmf+rLikpQf/+/XHs2DFs3LgRly9fxuLFi6Gvr6/2NkxMTFBSUoKioiJMnDgRmzdvVtnWxo0b0bBhQ/Tq1Usr+1SW4uJibNiwAQDUDoLx8fFISEjA7NmzVYLHvykUCgCPw01aWhqOHDnyVHXeuXMHe/bskb4TRm4MH0REVO0MDAwQHR2NDRs2wNraGt26dcPcuXNx/vx5AMDBgwfxxx9/YMeOHejTpw8aN26MgQMHon///mqt/9q1a1i7di06duwICwsLabzFrl27pD7R0dEICgqSPsi1RalUSuHK0NAQU6dOxZdffokmTZqo9Pv888+lfubm5pg1axYA4OrVqwCgcpbk7t27Kn0///xzAMDw4cMREBCAHj16wMnJCUOHDsVnn32G7OzsSut88g2+JiYmaNSoERQKBZYvX66tt0EjDB9ERCQLf39/pKamYvfu3ejXrx9iY2PRvn17REdHIz4+Ho0aNYKnp6fa63vyoW9qagovLy84ODhg06ZNAABjY2OMGTMG69evBwCcPXsWFy9eRFBQUJXrX7RokUogSE5OBvD4ckZ8fDzi4+Nx7tw5LFq0CFOmTMGePXtUlh81apTULz4+HqGhoeVuy9bWVupnbW2NgoICAIC+vj6ioqKQkpKCpUuXomHDhli0aBFatmyJtLQ0AFCpccqUKdI63377bcTHx+P8+fM4dOgQAGDAgAGVjhWpDhzzQUREsjE2NkafPn3Qp08fzJ8/HxMnTkRYWBhmz56t8bosLCxw9uxZ6OnpwcnJqdSA1IkTJ+K5555DSkoKoqKi0KtXr6f6JtspU6ZgxIgR0mtnZ2cAj+8eadq0qdTepk0b/PLLL1iyZAkGDRoktVtZWan0e6JZs2YAgMTERLRr1w7A45DxpO+TQan/1rBhQ4wZMwZjxozBBx98AE9PT6xdu1a6C+YJS0tL6e92dnbSOps1a4aVK1fC19cXv/32G/z8/DR+P54GwwcREelMixYtEBMTgzZt2iAlJQVXr15V++zHfz/0/6t169bo2LEj1q1bh82bN6sMPq0KGxsb2NjYqNVXX18f//zzj1p927VrB29vbyxbtgwjRowod9xHeerXrw8nJyfk5OQAQIXvyX9rBKB2ndrE8EFERNXu/v37GD58OMaPH482bdrAwsICZ86cwdKlSzF48GD06NED3bt3h7+/P5YvX46mTZviypUrUCgUat3JUZ6JEydi2rRpMDMzw9ChQ9Ve7sKFC7CwsJBeKxQKtG3btsy+Qgikp6cDePxBfuDAAezfvx8LFixQa1sKhQJRUVHo06cPunXrhtDQUDRv3hyFhYU4cuQI7t27JwWFL774AvHx8Rg6dCiaNGmCvLw8fPPNN7h06RJWr15d4XYePnyI9PR0CCFw+/ZtvPPOO7C3t0fXrl3VqlObGD6IiOqIhISEGrsdc3Nz+Pj4YMWKFbhx4wYKCwvh4uKCSZMmYe7cuQCAH374AbNnz0ZAQABycnKkW22fRkBAAEJCQhAQEABjY2O1l+vevbvKa319fRQVFZXZNzs7G05OTgAAIyMjuLm5ITw8HHPmzFF7e126dEFcXBwWLVqE4OBgpKenw8zMDG3btsWKFSswfvx4AEDnzp3x+++/Y8qUKUhNTYW5uTlatmyJmJgY6Rbm8ixYsEAKRPb29ujUqRN++eUX2Nraql2ntijEk5uHa4js7GxYWVlBqVSqXKsiIlVnz55Fhw4dEDV/NLzcHHRdjloSb2Vg3AcbERcXV2OfR1HT5eXlISkpCR4eHtKHaW15wqku3Lx5E02aNMHp06dL/cy5u7sjJCQEISEhuimuhgoKCkJWVlaZj1cv6+fvCU0+v3nmg4iolnN1dUVCQkKd/m4XTRUWFuL+/fuYN28eunTpUm7YnTNnDubNm4c7d+5Ij1h/Vh09ehT9+/dHfn4+BgwYUK3bYvggIqoDXF1da3QYkNuxY8fw4osvwtPTE9u3by+zz+HDh1FYWAgAKuM7nlUdO3aU7pQxNzev1m0xfBARUZ3Ts2dPVDaq4Gluu62LTExM1L5T5mnxIWNEREQkK4YPIqJaqIbdK0DPCG393GkUPt5//30oFAqV6cnX/gKPR8EGBwfD1tYW5ubm8Pf3R0ZGhlYKJSKi/3sw1JPHbRPJ6d+PeX8aGo/5aNmyJQ4ePPh/K/jXY19nzJiBn376Cdu2bYOVlRWmTZuGYcOG4dixY09VJBERPWZgYABTU1Pcu3cP9erV0/hpmERVVVJSgnv37sHU1LTMR75rQuOlDQwM4OjoWKpdqVTi66+/xubNm6WvK46KikLz5s1x8uRJdOnS5akKJSKix0/DdHJyQlJSEm7duqXrcugZo6enB1dX16f+ZmCNw8e1a9fg7OwMY2Nj+Pr6IiIiAq6uroiLi0NhYaHKl9N4e3vD1dUVJ06cKDd85OfnIz8/X3qtztcCExE9ywwNDdGsWTNeeiHZGRoaauVsm0bhw8fHB9HR0fDy8kJaWhoWLlyIF154ARcvXkR6ejoMDQ1hbW2tsoyDg4P0zPuyREREYOHChVUqnojoWaWnp6fR48KJahKNwkf//v2lv7dp0wY+Pj5wc3PD1q1bS32VsbpCQ0Mxc+ZM6XV2djZcXFyqtC4iIiKq+Z7q3Im1tTU8PT1x/fp1ODo6oqCgAFlZWSp9MjIyyhwj8oSRkREsLS1VJiIiIqq7nip8PHr0CDdu3ICTkxM6dOiAevXq4dChQ9L8xMREJCcnw9fX96kLJSIiorpBo8sus2fPxqBBg+Dm5obU1FSEhYVBX18fAQEBsLKywoQJEzBz5kzY2NjA0tIS06dPh6+vL+90ISIiIolG4SMlJQUBAQG4f/8+7O3t8fzzz+PkyZOwt7cHAKxYsQJ6enrw9/dHfn4++vbti88//7xaCiciIqLaSaPwsWXLlgrnGxsbIzIyEpGRkU9VFBEREdVdfDQeERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxfBBREREsmL4ICIiIlkxfBAREZGsGD6IiIhIVgwfREREJCuGDyIiIpKVga4LIKJnT0JCgq5L0IidnR1cXV11XQZRncHwQUSyua/MAQCMHj1ax5VoxtTUFAkJCQwgRFrC8EFEsnmYmwcAmDykG3xbe+i4GvXcTHuAhV/9jMzMTIYPIi1h+CAi2TnZWcLLzUHXZRCRjnDAKREREcmK4YOIiIhkxfBBREREsmL4ICIiIlkxfBAREZGsGD6IiIhIVgwfREREJCuGDyIiIpIVwwcRERHJiuGDiIiIZMXwQURERLJi+CAiIiJZMXwQERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxfBBREREsmL4ICIiIlk9VfhYvHgxFAoFQkJCpLa8vDwEBwfD1tYW5ubm8Pf3R0ZGxtPWSURERHVElcPH6dOn8cUXX6BNmzYq7TNmzMCePXuwbds2HD58GKmpqRg2bNhTF0pERER1Q5XCx6NHjzBq1CisW7cO9evXl9qVSiW+/vprLF++HL169UKHDh0QFRWF48eP4+TJk1ormoiIiGqvKoWP4OBgDBgwAH5+firtcXFxKCwsVGn39vaGq6srTpw4Uea68vPzkZ2drTIRERFR3WWg6QJbtmzB2bNncfr06VLz0tPTYWhoCGtra5V2BwcHpKenl7m+iIgILFy4UNMyiIiIqJbS6MzH7du38dZbb2HTpk0wNjbWSgGhoaFQKpXSdPv2ba2sl4iIiGomjcJHXFwc7t69i/bt28PAwAAGBgY4fPgwVq1aBQMDAzg4OKCgoABZWVkqy2VkZMDR0bHMdRoZGcHS0lJlIiIiorpLo8suvXv3xoULF1Taxo0bB29vb8yZMwcuLi6oV68eDh06BH9/fwBAYmIikpOT4evrq72qiYiIqNbSKHxYWFigVatWKm1mZmawtbWV2idMmICZM2fCxsYGlpaWmD59Onx9fdGlSxftVU1ERES1lsYDTiuzYsUK6Onpwd/fH/n5+ejbty8+//xzbW+GiIiIaqmnDh+xsbEqr42NjREZGYnIyMinXTURERHVQfxuFyIiIpIVwwcRERHJiuGDiIiIZMXwQURERLJi+CAiIiJZMXwQERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxfBBREREsmL4ICIiIlkxfBAREZGsGD6IiIhIVgwfREREJCuGDyIiIpIVwwcRERHJiuGDiIiIZMXwQURERLJi+CAiIiJZMXwQERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxfBBREREsmL4ICIiIlkxfBAREZGsGD6IiIhIVhqFjzVr1qBNmzawtLSEpaUlfH19sXfvXml+Xl4egoODYWtrC3Nzc/j7+yMjI0PrRRMREVHtpVH4aNSoERYvXoy4uDicOXMGvXr1wuDBg3Hp0iUAwIwZM7Bnzx5s27YNhw8fRmpqKoYNG1YthRMREVHtZKBJ50GDBqm8/uijj7BmzRqcPHkSjRo1wtdff43NmzejV69eAICoqCg0b94cJ0+eRJcuXbRXNREREdVaVR7zUVxcjC1btiAnJwe+vr6Ii4tDYWEh/Pz8pD7e3t5wdXXFiRMntFIsERER1X4anfkAgAsXLsDX1xd5eXkwNzfHzp070aJFC8THx8PQ0BDW1tYq/R0cHJCenl7u+vLz85Gfny+9zs7O1rQkIiIiqkU0Dh9eXl6Ij4+HUqnE9u3bERgYiMOHD1e5gIiICCxcuLDKyxNpzW8Ruq5AM1fv6LoCIqIq0fiyi6GhIZo2bYoOHTogIiICbdu2xaeffgpHR0cUFBQgKytLpX9GRgYcHR3LXV9oaCiUSqU03b59W+OdICIiotrjqZ/zUVJSgvz8fHTo0AH16tXDoUOHpHmJiYlITk6Gr69vucsbGRlJt+4+mYiIiKju0uiyS2hoKPr37w9XV1c8fPgQmzdvRmxsLPbv3w8rKytMmDABM2fOhI2NDSwtLTF9+nT4+vryThciIiKSaBQ+7t69i7FjxyItLQ1WVlZo06YN9u/fjz59+gAAVqxYAT09Pfj7+yM/Px99+/bF559/Xi2FExERUe2kUfj4+uuvK5xvbGyMyMhIREZGPlVRREREVHfxu12IiIhIVgwfREREJCuGDyIiIpIVwwcRERHJiuGDiIiIZMXwQURERLJi+CAiIiJZMXwQERGRrBg+iIiISFYMH0RERCQrjR6vTlSXJWdkIVOZo+sy1JaQfE/XJRARVQnDBxGA5ORkNB+3Arl5hbouRWOZyhx46boIIiINMHwQAcjMzERuXiE2zh2B5q72ui5HLT//cRXz1x/Ao9w8XZdCRKQRhg+if2nuao/2ng11XYZaeNmFiGorDjglIiIiWTF8EBERkax42YWISA0JCQm6LkFjdnZ2cHV11XUZRKUwfBARVeD+/7/9evTo0TquRHOmpqZISEhgAKEah+GDiKgCD///3USTh3SDb2sPHVejvptpD7Dwq5+RmZnJ8EE1DsMHEZEanOws4eXmoOsyiOoEDjglIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxfBBREREsmL4ICIiIlkxfBAREZGsGD6IiIhIVgwfREREJCuGDyIiIpIVwwcRERHJiuGDiIiIZMXwQURERLJi+CAiIiJZMXwQERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK43CR0REBDp16gQLCws0aNAAQ4YMQWJiokqfvLw8BAcHw9bWFubm5vD390dGRoZWiyYiIqLaS6PwcfjwYQQHB+PkyZM4cOAACgsL8dJLLyEnJ0fqM2PGDOzZswfbtm3D4cOHkZqaimHDhmm9cCIiIqqdDDTpvG/fPpXX0dHRaNCgAeLi4tC9e3colUp8/fXX2Lx5M3r16gUAiIqKQvPmzXHy5El06dJFe5UTERFRrfRUYz6USiUAwMbGBgAQFxeHwsJC+Pn5SX28vb3h6uqKEydOlLmO/Px8ZGdnq0xERERUd2l05uPfSkpKEBISgm7duqFVq1YAgPT0dBgaGsLa2lqlr4ODA9LT08tcT0REBBYuXFjVMoi06vwdJfINjHVdhlqu3X2o6xKIiKqkymc+goODcfHiRWzZsuWpCggNDYVSqZSm27dvP9X6iIiIqGar0pmPadOm4ccff8SRI0fQqFEjqd3R0REFBQXIyspSOfuRkZEBR0fHMtdlZGQEIyOjqpRBREREtZBGZz6EEJg2bRp27tyJX3/9FR4eHirzO3TogHr16uHQoUNSW2JiIpKTk+Hr66udiomIiKhW0+jMR3BwMDZv3oxdu3bBwsJCGsdhZWUFExMTWFlZYcKECZg5cyZsbGxgaWmJ6dOnw9fXl3e6EBEREQANw8eaNWsAAD179lRpj4qKQlBQEABgxYoV0NPTg7+/P/Lz89G3b198/vnnWimWiIiIaj+NwocQotI+xsbGiIyMRGRkZJWLIiIiorqL3+1CREREsmL4ICIiIlkxfBAREZGsGD6IiIhIVgwfREREJCuGDyIiIpIVwwcRERHJiuGDiIiIZMXwQURERLJi+CAiIiJZMXwQERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxfBBREREsmL4ICIiIlkxfBAREZGsGD6IiIhIVgwfREREJCuGDyIiIpIVwwcRERHJiuGDiIiIZMXwQURERLJi+CAiIiJZMXwQERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERycpA1wUQEVH1SUhI0HUJGrGzs4Orq6uuy6BqxvBBRFQH3VfmAABGjx6t40o0Y2pqioSEBAaQOo7hg4ioDnqYmwcAmDykG3xbe+i4GvXcTHuAhV/9jMzMTIaPOo7hg4ioDnOys4SXm4OuyyBSwQGnREREJCuGDyIiIpIVwwcRERHJimM+qHr8FqHrCjRz9Y6uKyAiembwzAcRERHJSuPwceTIEQwaNAjOzs5QKBSIiYlRmS+EwIIFC+Dk5AQTExP4+fnh2rVr2qqXiIiIajmNw0dOTg7atm2LyMjIMucvXboUq1atwtq1a3Hq1CmYmZmhb9++yMvLe+piiYiIqPbTeMxH//790b9//zLnCSGwcuVKzJs3D4MHDwYAfPPNN3BwcEBMTAxeffXVp6uWiIiIaj2tjvlISkpCeno6/Pz8pDYrKyv4+PjgxIkTZS6Tn5+P7OxslYmIiIjqLq2Gj/T0dACAg4Pq0/QcHBykef8VEREBKysraXJxcdFmSURERFTD6Pxul9DQUCiVSmm6ffu2rksiIiKiaqTV8OHo6AgAyMjIUGnPyMiQ5v2XkZERLC0tVSYiIiKqu7QaPjw8PODo6IhDhw5JbdnZ2Th16hR8fX21uSkiIiKqpTS+2+XRo0e4fv269DopKQnx8fGwsbGBq6srQkJC8OGHH6JZs2bw8PDA/Pnz4ezsjCFDhmizbiIiIqqlNA4fZ86cwYsvvii9njlzJgAgMDAQ0dHReOedd5CTk4PJkycjKysLzz//PPbt2wdjY2PtVU1ERES1lsbho2fPnhBClDtfoVAgPDwc4eHhT1UYERER1U06v9uFiIiIni0MH0RERCQrhg8iIiKSFcMHERERyUrjAadE6kjOyEKmMkfXZagtIfmerksgInpmMHyQ1iUnJ8MrcDnyCop0XYrGMpU58NJ1EUREdRzDB2ldZmYm8gqKEDbxZbg72ei6HLWcuJCEL2OO4VFunq5LISKq8xg+qNq4O9nAy82h8o41wM20+7ougYjomcEBp0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawYPoiIiEhWDB9EREQkK4YPIiIikhXDBxEREcmK4YOIiIhkxfBBREREsmL4ICIiIlkZ6LoAUsNvEbquQDNX7+i6AiIiqsF45oOIiIhkxfBBREREsmL4ICIiIlkxfBAREZGsOOCUiIhqlISEBF2XoBE7Ozu4urrquoxaheGDiIhqhPvKHADA6NGjdVyJZkxNTZGQkMAAogGGDyIiqhEe5uYBACYP6Qbf1h46rkY9N9MeYOFXPyMzM5PhQwMMH0REVKM42VnCy81B12VQNeKAUyIiIpIVwwcRERHJiuGDiIiIZMXwQURERLJi+CAiIiJZMXwQERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkawMdF0AVS45IwuZyhxdl6G2hOR7ui6BiIhqsGcvfPwWoesKNJKckQWvwOXIKyjSdSkay1TmwEvXRRARUY1TbeEjMjISH3/8MdLT09G2bVusXr0anTt3rq7N1VmZyhzkFRQhbOLLcHey0XU5ajlxIQlfxhzDo9w8XZdCREQ1ULWEj++//x4zZ87E2rVr4ePjg5UrV6Jv375ITExEgwYNqmOTdZ67kw283Bx0XYZabqbd13UJRERUg1XLgNPly5dj0qRJGDduHFq0aIG1a9fC1NQU69evr47NERERUS2i9fBRUFCAuLg4+Pn5/d9G9PTg5+eHEydOaHtzREREVMto/bJLZmYmiouL4eCgeonAwcEBV65cKdU/Pz8f+fn50mulUgkAyM7O1nZpAID02/eQ/uBRtay7OiSmZD7+81YG/skv0HE16rmV9kD6M/7qbR1Xox7WLA/WLJ/aWHdtrDk5/W8AQFxcHB49qj2fLY6OjnB0dNTqOp98bgshKu8stOzOnTsCgDh+/LhK+9tvvy06d+5cqn9YWJgAwIkTJ06cOHGqA9Pt27crzQpaP/NhZ2cHfX19ZGRkqLRnZGSUmbJCQ0Mxc+ZM6XVJSQkePHgAW1tbKBQKrdaWnZ0NFxcX3L59G5aWllpdd01Q1/cPqPv7yP2r/er6PnL/ar/q2kchBB4+fAhnZ+dK+2o9fBgaGqJDhw44dOgQhgwZAuBxoDh06BCmTZtWqr+RkRGMjIxU2qytrbVdlgpLS8s6+0MF1P39A+r+PnL/ar+6vo/cv9qvOvbRyspKrX7VcqvtzJkzERgYiI4dO6Jz585YuXIlcnJyMG7cuOrYHBEREdUi1RI+Ro4ciXv37mHBggVIT0/Hc889h3379pUahEpERETPnmp7wum0adPKvMyiS0ZGRggLCyt1maeuqOv7B9T9feT+1X51fR+5f7VfTdhHhRDq3BNDREREpB3V8oRTIiIiovIwfBAREZGsGD6IiIhIVgwfREREJKs6FT4++ugjdO3aFaampuU+qCw5ORkDBgyAqakpGjRogLfffhtFRUUVrvfBgwcYNWoULC0tYW1tjQkTJtSIZ/jHxsZCoVCUOZ0+fbrc5Xr27Fmq/5QpU2SsXH3u7u6lal28eHGFy+Tl5SE4OBi2trYwNzeHv79/qSfu1hQ3b97EhAkT4OHhARMTEzRp0gRhYWEoKKj4e3xq8jGMjIyEu7s7jI2N4ePjgz/++KPC/tu2bYO3tzeMjY3RunVr/PzzzzJVqrmIiAh06tQJFhYWaNCgAYYMGYLExMQKl4mOji51rIyNjWWqWDPvv/9+qVq9vb0rXKY2HT+g7H9TFAoFgoODy+xf04/fkSNHMGjQIDg7O0OhUCAmJkZlvhACCxYsgJOTE0xMTODn54dr165Vul5Nf481VafCR0FBAYYPH46pU6eWOb+4uBgDBgxAQUEBjh8/jg0bNiA6OhoLFiyocL2jRo3CpUuXcODAAfz44484cuQIJk+eXB27oJGuXbsiLS1NZZo4cSI8PDzQsWPHCpedNGmSynJLly6VqWrNhYeHq9Q6ffr0CvvPmDEDe/bswbZt23D48GGkpqZi2LBhMlWrmStXrqCkpARffPEFLl26hBUrVmDt2rWYO3dupcvWxGP4/fffY+bMmQgLC8PZs2fRtm1b9O3bF3fv3i2z//HjxxEQEIAJEybg3LlzGDJkCIYMGYKLFy/KXLl6Dh8+jODgYJw8eRIHDhxAYWEhXnrpJeTk5FS4nKWlpcqxunXrlkwVa65ly5Yqtf7+++/l9q1txw8ATp8+rbJ/Bw4cAAAMHz683GVq8vHLyclB27ZtERkZWeb8pUuXYtWqVVi7di1OnToFMzMz9O3bF3l5eeWuU9Pf4yrRyrfJ1TBRUVHCysqqVPvPP/8s9PT0RHp6utS2Zs0aYWlpKfLz88tc1+XLlwUAcfr0aalt7969QqFQiDt37mi99qdRUFAg7O3tRXh4eIX9evToId566y15inpKbm5uYsWKFWr3z8rKEvXq1RPbtm2T2hISEgQAceLEiWqoUPuWLl0qPDw8KuxTU49h586dRXBwsPS6uLhYODs7i4iIiDL7jxgxQgwYMEClzcfHR7z++uvVWqe23L17VwAQhw8fLrdPef8e1URhYWGibdu2avev7cdPCCHeeust0aRJE1FSUlLm/Np0/ACInTt3Sq9LSkqEo6Oj+Pjjj6W2rKwsYWRkJL777rty16Pp73FV1KkzH5U5ceIEWrdurfKk1b59+yI7OxuXLl0qdxlra2uVMwl+fn7Q09PDqVOnqr1mTezevRv3799X6zH2mzZtgp2dHVq1aoXQ0FDk5ubKUGHVLF68GLa2tmjXrh0+/vjjCi+TxcXFobCwEH5+flKbt7c3XF1dceLECTnKfWpKpRI2NjaV9qtpx7CgoABxcXEq772enh78/PzKfe9PnDih0h94/DtZm44VgEqP16NHj+Dm5gYXFxcMHjy43H9vaoJr167B2dkZjRs3xqhRo5CcnFxu39p+/AoKCrBx40aMHz++wi8yrU3H79+SkpKQnp6ucoysrKzg4+NT7jGqyu9xVVTbE05rovT09FKPeH/yOj09vdxlGjRooNJmYGAAGxubcpfRla+//hp9+/ZFo0aNKuz32muvwc3NDc7Ozjh//jzmzJmDxMRE7NixQ6ZK1ffmm2+iffv2sLGxwfHjxxEaGoq0tDQsX768zP7p6ekwNDQsNebHwcGhxh2vsly/fh2rV6/GsmXLKuxXE49hZmYmiouLy/wdu3LlSpnLlPc7WRuOVUlJCUJCQtCtWze0atWq3H5eXl5Yv3492rRpA6VSiWXLlqFr1664dOlSpb+rcvPx8UF0dDS8vLyQlpaGhQsX4oUXXsDFixdhYWFRqn9tPn4AEBMTg6ysLAQFBZXbpzYdv/96chw0OUZV+T2uihofPt59910sWbKkwj4JCQmVDoqqTaqyzykpKdi/fz+2bt1a6fr/PV6ldevWcHJyQu/evXHjxg00adKk6oWrSZP9mzlzptTWpk0bGBoa4vXXX0dERESNfvxxVY7hnTt30K9fPwwfPhyTJk2qcFldH0MCgoODcfHixQrHRACAr68vfH19pdddu3ZF8+bN8cUXX+CDDz6o7jI10r9/f+nvbdq0gY+PD9zc3LB161ZMmDBBh5VVj6+//hr9+/ev8Cvga9Pxq01qfPiYNWtWhakUABo3bqzWuhwdHUuN2H1yF4Sjo2O5y/x3kE1RUREePHhQ7jJPqyr7HBUVBVtbW7zyyisab8/HxwfA4/91y/HB9TTH1MfHB0VFRbh58ya8vLxKzXd0dERBQQGysrJUzn5kZGRU2/Eqi6b7mJqaihdffBFdu3bFl19+qfH25D6GZbGzs4O+vn6pO4sqeu8dHR016l9TTJs2TRp8run/fuvVq4d27drh+vXr1VSd9lhbW8PT07PcWmvr8QOAW7du4eDBgxqfLaxNx+/JccjIyICTk5PUnpGRgeeee67MZarye1wlWhs9UoNUNuA0IyNDavviiy+EpaWlyMvLK3NdTwacnjlzRmrbv39/jRpwWlJSIjw8PMSsWbOqtPzvv/8uAIg///xTy5Vp38aNG4Wenp548OBBmfOfDDjdvn271HblypUaPeA0JSVFNGvWTLz66quiqKioSuuoKcewc+fOYtq0adLr4uJi0bBhwwoHnA4cOFClzdfXt8YOWCwpKRHBwcHC2dlZXL16tUrrKCoqEl5eXmLGjBlark77Hj58KOrXry8+/fTTMufXtuP3b2FhYcLR0VEUFhZqtFxNPn4oZ8DpsmXLpDalUqnWgFNNfo+rVKvW1lQD3Lp1S5w7d04sXLhQmJubi3Pnzolz586Jhw8fCiEe/9C0atVKvPTSSyI+Pl7s27dP2Nvbi9DQUGkdp06dEl5eXiIlJUVq69evn2jXrp04deqU+P3330WzZs1EQECA7PtXnoMHDwoAIiEhodS8lJQU4eXlJU6dOiWEEOL69esiPDxcnDlzRiQlJYldu3aJxo0bi+7du8tddqWOHz8uVqxYIeLj48WNGzfExo0bhb29vRg7dqzU57/7J4QQU6ZMEa6uruLXX38VZ86cEb6+vsLX11cXu1CplJQU0bRpU9G7d2+RkpIi0tLSpOnffWrLMdyyZYswMjIS0dHR4vLly2Ly5MnC2tpausNszJgx4t1335X6Hzt2TBgYGIhly5aJhIQEERYWJurVqycuXLigq12o0NSpU4WVlZWIjY1VOVa5ublSn//u48KFC8X+/fvFjRs3RFxcnHj11VeFsbGxuHTpki52oUKzZs0SsbGxIikpSRw7dkz4+fkJOzs7cffuXSFE7T9+TxQXFwtXV1cxZ86cUvNq2/F7+PCh9FkHQCxfvlycO3dO3Lp1SwghxOLFi4W1tbXYtWuXOH/+vBg8eLDw8PAQ//zzj7SOXr16idWrV0uvK/s91oY6FT4CAwMFgFLTb7/9JvW5efOm6N+/vzAxMRF2dnZi1qxZKsn3t99+EwBEUlKS1Hb//n0REBAgzM3NhaWlpRg3bpwUaGqCgIAA0bVr1zLnJSUlqbwHycnJonv37sLGxkYYGRmJpk2birffflsolUoZK1ZPXFyc8PHxEVZWVsLY2Fg0b95cLFq0SOUs1X/3Twgh/vnnH/HGG2+I+vXrC1NTUzF06FCVD/OaJCoqqsyf2X+flKxtx3D16tXC1dVVGBoais6dO4uTJ09K83r06CECAwNV+m/dulV4enoKQ0ND0bJlS/HTTz/JXLH6yjtWUVFRUp//7mNISIj0fjg4OIiXX35ZnD17Vv7i1TBy5Ejh5OQkDA0NRcOGDcXIkSPF9evXpfm1/fg9sX//fgFAJCYmlppX247fk8+s/05P9qGkpETMnz9fODg4CCMjI9G7d+9S++3m5ibCwsJU2ir6PdYGhRBCaO8iDhEREVHFnqnnfBAREZHuMXwQERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFEsggKCsKQIUN0XQYR1QB8yBgRyUKpVEIIofKFf0T0bGL4ICIiIlnxsgsRadX27dvRunVrmJiYwNbWFn5+fsjJySl12eXhw4cYNWoUzMzM4OTkhBUrVqBnz54ICQnRWe1EJA+GDyLSmrS0NAQEBGD8+PFISEhAbGwshg0bhrJOsM6cORPHjh3D7t27ceDAARw9ehRnz57VQdVEJDcDXRdARHVHWloaioqKMGzYMLi5uQEAWrduXarfw4cPsWHDBmzevBm9e/cGAERFRcHZ2VnWeolIN3jmg4i0pm3btujduzdat26N4cOHY926dfj7779L9fvrr79QWFiIzp07S21WVlbw8vKSs1wi0hGGDyLSGn19fRw4cAB79+5FixYtsHr1anh5eSEpKUnXpRFRDcLwQURapVAo0K1bNyxcuBDnzp2DoaEhdu7cqdKncePGqFevHk6fPi21KZVKXL16Ve5yiUgHOOaDiLTm1KlTOHToEF566SU0aNAAp06dwr1799C8eXOcP39e6mdhYYHAwEC8/fbbsLGxQYMGDRAWFgY9PT0oFAod7gERyYFnPohIaywtLXHkyBG8/PLL8PT0xLx58/DJJ5+gf//+pfouX74cvr6+GDhwIPz8/NCtWzc0b94cxsbGOqiciOTEh4wRUY2Qk5ODhg0b4pNPPsGECRN0XQ4RVSNediEinTh37hyuXLmCzp07Q6lUIjw8HAAwePBgHVdGRNWN4YOIdGbZsmVITEyEoaEhOnTogKNHj8LOzk7XZRFRNeNlFyIiIpIVB5wSERGRrBg+iIiISFYMH0RERCQrhg8iIiKSFcMHERERyYrhg4iIiGTF8EFERESyYvggIiIiWTF8EBERkaz+H1xoeL0SYEekAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(f\"{np.sum(np.isnan(sigs_bfgs))} BFGS fits failed\")\n",
    "plt.hist(sigs_migrad, alpha=0.5, bins=10, range=(-10, 10), label=\"Migrad\")\n",
    "plt.hist(sigs_constrained, alpha=0.5, bins=10, range=(-10, 10), label=\"SciPy[SLSQS]\")\n",
    "plt.hist(sigs_bfgs, bins=10, range=(-10, 10), fill=False, label=\"SciPy[L-BFGS-B]\")\n",
    "plt.xlabel(\"sig\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "933fcadf",
   "metadata": {},
   "source": [
    "In this example, the BFGS method actually failed much less than Migrad, but it still fails in some cases, while the constrained fit did not fail at all."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a2de7c9",
   "metadata": {},
   "source": [
    "## Speed comparison\n",
    "\n",
    "Since constrained fits are so useful, should you use them all the time? Probably not.\n",
    "\n",
    "Constrained fits are more computationally expensive. Satisfying extra constrains generally slows down convergence. Let's compare the speed of the three minimisers tested here. We set the strategy to 0, to avoid computing the Hessian matrix automatically, since we want to measure only the time used by the minimiser."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "84dfd255",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.strategy = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6e6c7b2b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.37 ms ± 555 µs per loop (mean ± std. dev. of 7 runs, 3 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n3 m.reset(); m.migrad()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bafa5a40",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.67 ms ± 713 µs per loop (mean ± std. dev. of 7 runs, 3 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n3 m.reset(); m.scipy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "05a2b917",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.4 ms ± 283 µs per loop (mean ± std. dev. of 7 runs, 3 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n3 m.reset(); m.scipy(constraints=NonlinearConstraint(lambda *par: model(x, *par)[1], 0, np.inf))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "365a71bf",
   "metadata": {},
   "source": [
    "Migrad is the fastest, followed by the L-BFGS-B method. The constrained fit is much slower.\n",
    "\n",
    "The constrained fit is much slower, since it has to do more work. Why Migrad is faster than L-BFGS-B is not so obvious. There are some general reasons for that, but there may be cases where L-BFGS-B performs better.\n",
    "\n",
    "Migrad is comparably fast because of its smart stopping criterion. Migrad stops the fit as soon as the improvement of the fitted parameters become small compared to their uncertainties. Migrad was explicitly designed for statistical fits, where the cost function is a log-likelihood or least-squares function. Since it assumes that, it can stops the fit as soon as the parameter improvements become negligible compared to the parameter uncertainty, which is given by the inverse of its internal approximation of the Hessian matrix.\n",
    "\n",
    "The SciPy minimisers do not expect the cost function to be a log-likelihood or least-squares and thus cannot assume that the Hessian matrix has a special meaning. Instead they tend to optimise until they hit the limits of machine precision. This is the main reason why the L-BFGS-B method is slower. You can also see this in the benchmark section of the documentation.\n",
    "\n",
    "We can force Migrad to do something similar by setting the tolerance to a tiny value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "70327e76",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.tol = 1e-20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3cb6d132",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.89 ms ± 364 µs per loop (mean ± std. dev. of 7 runs, 3 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n3 m.reset(); m.migrad()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "178aaef6",
   "metadata": {},
   "source": [
    "Now the runtime of Migrad is closer to L-BFGS-B, but it is still faster in this case."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.14 ('venv': venv)",
   "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.10.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "bdbf20ff2e92a3ae3002db8b02bd1dd1b287e934c884beb29a73dced9dbd0fa3"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}