{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performance of cost functions\n",
    "\n",
    "This is not really a tutorial, but more of a benchmark of the builtin cost functions.\n",
    "\n",
    "We test the performance of the cost functions shipped with iminuit. We check that they produce unbiased results with proper variance. To do that, we generate normal distributed data many times and fit a normal distribution to each independent data set. The bias is computed from the averages of these reconstructed parameters. We also compute the mean of the estimated variance for each data set, which should converge to 1.\n",
    "\n",
    "Since we do the fit many times, we do not use implementations of the pdf and cdf of a normal distribution from `scipy.stats`, but Numba-accelerated versions from the `numba-stats` package. For the binned fits, we compute histograms of the data with $3 + n/10$ equidistant bins, where $n$ is the sample size.\n",
    "\n",
    "Disclaimer: This tutorial is targeted at experts, please read the code to understand what is going on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum-likelihood fits\n",
    "\n",
    "Here we check that the different maximum-likelihood cost functions produce asymptotically unbiased results with the expected variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from iminuit import Minuit\n",
    "from iminuit.cost import (\n",
    "    UnbinnedNLL,\n",
    "    BinnedNLL,\n",
    "    ExtendedUnbinnedNLL,\n",
    "    ExtendedBinnedNLL,\n",
    "    LeastSquares,\n",
    ")\n",
    "from argparse import Namespace\n",
    "import numba as nb\n",
    "import math\n",
    "from numba_stats import norm\n",
    "import joblib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/tl/pv6mt7z17tz0stm1fjfg01cc0000gn/T/ipykernel_13411/675801181.py:31: RuntimeWarning: invalid value encountered in true_divide\n",
      "/var/folders/tl/pv6mt7z17tz0stm1fjfg01cc0000gn/T/ipykernel_13411/675801181.py:56: RuntimeWarning: divide by zero encountered in true_divide\n",
      "/Users/hdembinski/Extern/iminuit/venv/lib/python3.8/site-packages/numpy/lib/function_base.py:1292: RuntimeWarning: invalid value encountered in subtract\n",
      "  a = op(a[slice1], a[slice2])\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 20 need to re-try [(True, True), (True, True), (False, True), (False, False)]\n"
     ]
    }
   ],
   "source": [
    "n_tries = 100  # increase this to get less scattering\n",
    "\n",
    "n_pts = np.array((10, 30, 100, 300, 1000, 3000, 10000))\n",
    "\n",
    "truth = Namespace(mu=0, sigma=1)\n",
    "\n",
    "\n",
    "# function that runs random experiments with sample size n\n",
    "@joblib.delayed\n",
    "def compute(n):\n",
    "    rng = np.random.default_rng(n)\n",
    "    np.random.seed(n)\n",
    "    u_nll = []\n",
    "    b_nll = []\n",
    "    e_u_nll = []\n",
    "    e_b_nll = []\n",
    "    for i_try in range(n_tries):\n",
    "        while True:\n",
    "            k = 2 * rng.poisson(n)\n",
    "            x = rng.normal(truth.mu, truth.sigma, size=k)\n",
    "            x = x[np.abs(x) < 2]\n",
    "            x = x[:k]\n",
    "            xrange = np.array((-2.0, 2.0))\n",
    "            nh, xe = np.histogram(x, bins=3 + n // 10, range=xrange)\n",
    "            m = [\n",
    "                # model must be a normalized pdf\n",
    "                Minuit(\n",
    "                    UnbinnedNLL(\n",
    "                        x,\n",
    "                        lambda x, mu, sigma: (\n",
    "                            norm.pdf(x, mu, sigma) / np.diff(norm.cdf(xrange, mu, sigma))\n",
    "                        ),\n",
    "                    ),\n",
    "                    mu=truth.mu,\n",
    "                    sigma=truth.sigma,\n",
    "                ),\n",
    "                # model must be a function that returns the integral over the scaled pdf and the scaled pdf\n",
    "                Minuit(\n",
    "                    ExtendedUnbinnedNLL(\n",
    "                        x,\n",
    "                        lambda x, n, mu, sigma: (\n",
    "                            n * np.diff(norm.cdf(xrange, mu, sigma)),\n",
    "                            n * norm.pdf(x, mu, sigma),\n",
    "                        ),\n",
    "                    ),\n",
    "                    n=n,\n",
    "                    mu=truth.mu,\n",
    "                    sigma=truth.sigma,\n",
    "                ),\n",
    "                # model must be a normalized cdf up to an arbitrary additive constant (only differences are used)\n",
    "                Minuit(\n",
    "                    BinnedNLL(\n",
    "                        nh,\n",
    "                        xe,\n",
    "                        lambda x, mu, sigma: (\n",
    "                            norm.cdf(x, mu, sigma) / np.diff(norm.cdf(xrange, mu, sigma))\n",
    "                        ),\n",
    "                    ),\n",
    "                    mu=truth.mu,\n",
    "                    sigma=truth.sigma,\n",
    "                ),\n",
    "                # model must be a scaled cdf up to an arbitrary additive constant (only differences are used)\n",
    "                Minuit(\n",
    "                    ExtendedBinnedNLL(\n",
    "                        nh, xe, lambda x, n, mu, sigma: n * norm.cdf(x, mu, sigma)\n",
    "                    ),\n",
    "                    n=n,\n",
    "                    mu=truth.mu,\n",
    "                    sigma=truth.sigma,\n",
    "                ),\n",
    "            ]\n",
    "            for mi in m:\n",
    "                mi.limits[\"sigma\"] = (1e-3, None)\n",
    "                if \"n\" in mi.parameters:\n",
    "                    mi.limits[\"n\"] = (0, None)\n",
    "\n",
    "            # only accept a random data set when all fits converged ok\n",
    "            all_good = True\n",
    "            for mi in m:\n",
    "                mi.migrad()\n",
    "                mi.hesse()\n",
    "                if not mi.valid or not mi.accurate:\n",
    "                    all_good = False\n",
    "                    break\n",
    "            if all_good:\n",
    "                break\n",
    "            print(f\"{n} {i_try} need to re-try {[(mi.valid, mi.accurate) for mi in m]}\")\n",
    "\n",
    "        # store parameter deviations and estimated variances for each pseudo-experiment\n",
    "        u_nll.append(\n",
    "            (\n",
    "                m[0].values[\"mu\"] - truth.mu,\n",
    "                m[0].errors[\"mu\"] ** 2,\n",
    "                m[0].values[\"sigma\"] - truth.sigma,\n",
    "                m[0].errors[\"sigma\"] ** 2,\n",
    "            )\n",
    "        )\n",
    "        e_u_nll.append(\n",
    "            (\n",
    "                m[1].values[\"n\"] - n,\n",
    "                m[1].errors[\"n\"] ** 2,\n",
    "                m[1].values[\"mu\"] - truth.mu,\n",
    "                m[1].errors[\"mu\"] ** 2,\n",
    "                m[1].values[\"sigma\"] - truth.sigma,\n",
    "                m[1].errors[\"sigma\"] ** 2,\n",
    "            )\n",
    "        )\n",
    "        b_nll.append(\n",
    "            (\n",
    "                m[2].values[\"mu\"] - truth.mu,\n",
    "                m[2].errors[\"mu\"] ** 2,\n",
    "                m[2].values[\"sigma\"] - truth.sigma,\n",
    "                m[2].errors[\"sigma\"] ** 2,\n",
    "            )\n",
    "        )\n",
    "        e_b_nll.append(\n",
    "            (\n",
    "                m[3].values[\"n\"] - n,\n",
    "                m[3].errors[\"n\"] ** 2,\n",
    "                m[3].values[\"mu\"] - truth.mu,\n",
    "                m[3].errors[\"mu\"] ** 2,\n",
    "                m[3].values[\"sigma\"] - truth.sigma,\n",
    "                m[3].errors[\"sigma\"] ** 2,\n",
    "            )\n",
    "        )\n",
    "\n",
    "    # means over pseudo-experiments are computed here\n",
    "    return (\n",
    "        np.mean(u_nll, axis=0),\n",
    "        np.mean(e_u_nll, axis=0),\n",
    "        np.mean(b_nll, axis=0),\n",
    "        np.mean(e_b_nll, axis=0),\n",
    "    )\n",
    "\n",
    "\n",
    "unbinned_nll = []\n",
    "extended_unbinned_nll = []\n",
    "binned_nll = []\n",
    "extended_binned_nll = []\n",
    "\n",
    "result = joblib.Parallel(-1)(compute(n) for n in n_pts)\n",
    "\n",
    "for a,b,c,d in result:\n",
    "    unbinned_nll.append(a)\n",
    "    extended_unbinned_nll.append(b)\n",
    "    binned_nll.append(c)\n",
    "    extended_binned_nll.append(d)\n",
    "\n",
    "unbinned_nll = np.transpose(unbinned_nll)\n",
    "extended_unbinned_nll = np.transpose(extended_unbinned_nll)\n",
    "binned_nll = np.transpose(binned_nll)\n",
    "extended_binned_nll = np.transpose(extended_binned_nll)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the measured bias as a point and the mean variance as an error bar. The deviations go down with $n^{-{1/2}}$, where $n$ is the sample size. We undo this for the plots by multiplying deviations with $n^{1/2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzEAAAH2CAYAAABA5AKlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6cUlEQVR4nO3dfZzcdX3v/deHEAkCkjaEUyDGQKWccqfIVk+VWq5KBZEWDtcBbdFTbC31eDxar4oVsNabauWkl629amsp1p+x6UoMkZhGQ0xNommMJgshCazRnDQQVtIkSxKJyeoSv9cfMztMNrN3mdmd+c6+no/HPJLv7/Y7+52Zz77nd7ORUkKSJEmScnFCszsgSZIkSWNhiJEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSuGGEmSJElZMcSoLUVEERF/Nsz8gxFx3kT2qbzfVRHx1onerySptUTEpyPiT5qw31sjYs1E71dqNEOMWlJEpIh48aBpH4yIf2rE9lNKp6aUtjdiW41Sfn4pIm6umnZiedqccnvIcFbrZyZJao6I2BERh8tfmu2LiKUR8cKB+Smlt6WUPtLMPg4WEXPKteQrg6b/U0R8sPz/KyPiySHWH/YLRKmRDDFSa3ka+FBETGl2RyRJdfuNlNKpwFnAfwD/X5P7M1qviIhXNrsT0nAMMcrSwDdBEfFHEbE7Ip6KiLcMWuyMiPhaRDwTEasj4kVV61eOWpS/OfpU+VuyZyLi2xHx84OWfVtEfD8i9peXjar5vxsR3eVv2h4ctJ9fj4jvRsSBiPgboLLeEJYBPwHeVMePR5LUQlJKfcBC4MKBadVHLUaqaaOoU/+5XO+ejoitg47oz4iIL0fEDyPiO0BlvWH8b+CjDXjq0rgxxChnPwecDpwD/B7wqYj4mar5twAfAc4ANgLzh9nWG4EPAT8DbOPYD+/rgF8CLgVuBq4GiIjrgTuBG4GZwDeBzvK8M4BFwPvLffg/wKtGeE4J+BPgTyNi6gjLSpIyEBHPB94ArBtmsZFqWs06FRGnAF8D/hk4s7zc30bEQGD6FNBH6WjQ75YfI/lb4Bci4qrRPD+pGQwxylk/8OGUUn9K6SvAQeCCqvlLU0rfSCn9GLgL+OXq85EH+VJK6TsppWcphZ2XDpr/8ZTS/pTSE8DKqvlvA/48pdRdXvdjwEvLR2OuBR5NKS1MKfUDfwXsGulJpZS+DOwBvAGAJOXtgYjYDxwAfh2YO8yyI9W0oerUdcCOlNJnU0rPppQeBu4Hbiqfmvx/Ax9IKf0opbQF+Nwo+n2YUkjy+ha1LEOMWtURYPCRiKmUPuQH9JY/zAccAk6tau8c+E9K6SCl603OHmJ/1eFi8HaGm/8i4JPl08z2l/cRlL5JO3tQH1J1ewTvpxS8po1yeUlS67khpTSd0mf5O4DVEfFzQyw7Uk0brg69YqAOlWvRLZSO7MwETuTo2vP4KPt+L/CfIuI3Rrm8NKEMMWpVTwBzBk07l9F/+AJUjrpExKnAzwI/qLtnR9sJ/EFKaXrV4+SU0lrgqUF9iOr2cFJKX6N0usDbG9xfSdIESykdSSktovQF3RUN3vxOYPWgOnRqSul/UDqq/yxH157Zo+zzTyidvvYRRr6eU5pwhhi1qvuA90fErIg4oXxe7m9QujBytK6NiCsi4nmUPoTXpZRGeyRktD4N3BERFwFExOkRcVN53lLgooi4MSJOBN5J6Zux0boLeG+N6VMiYlrV43lV8543aJ53OZOkJouS6yldz9Ld4M3/C6XrV94cEVPLj1+KiF9MKR2hdG3mByPi+eXrZH5nDNv+PKWjSNcMnjGo1kyruuHNcDVKahhDjFrVh4G1wBpgH6U7pdxSPp93tP4Z+FNKp3hdzjjc8Sul9CXgbuALEfFDYAvwuvK8vcBNwMeBXuB84N/GsO1/A75TY9b7KJ2vPPD4etW8RwfNG3zHNknSxFkSEQeBH1K6xuR3UkqPNnIHKaVngNdSuqD/B5ROO7sbOKm8yDsonXq2CyiAz45h20eAD1A6k6HaORxdaw7z3F3PhqtRUsNE6TR9SZIkScqDR2IkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVk5sxk7POOOMNGfOnGbsWpJU1tXVtTelNLPZ/WhF1ilJar7h6lRTQsycOXPYsGFDM3YtSSqLiMeb3YdWZZ2SpOYbrk55OpkkSZKkrBhiJEmSJGXFECNJkiQpK4YYSZIkSVkxxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVhoWYiJiSkQ8HBH/0qhtSpIkSdJgjTwS8y6gu4HbkyRJkqRjNCTERMQs4PXAvY3YniRJkiQNpVFHYv4KeC/w06EWiIjbImJDRGzYs2dPg3YrSVJjWKckKR91h5iIuA7YnVLqGm65lNI9KaWOlFLHzJkz692tJEkNZZ2SpHw04kjMq4DfjIgdwBeAX4uIf2rAdiVJkiTpGHWHmJTSHSmlWSmlOcAbga+nlN5Ud88kSZIkqQb/TowkSZKkrJzYyI2llFYBqxq5TUmSJEmq5pEYSZIkSVkxxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVgwxkiRJkrJiiJEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSuGGEmSJElZMcRIkiRJyoohRpIkSVJWDDGSJEmSsmKIkSRJkpQVQ4wkSZKkrBhiJEmSJGXFECNJkiQpK4YYSZIkSVkxxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZaXuEBMRL4yIlRHxWEQ8GhHvakTHJEmSJKmWExuwjWeBP0opPRQRpwFdEfG1lNJjDdi2JEmSJB2l7iMxKaWnUkoPlf//DNANnFPvdiVJkiSploZeExMRc4DLgG/XmHdbRGyIiA179uxp5G4lSapbw+rUZ19fekiSxk3DQkxEnArcD/xhSumHg+enlO5JKXWklDpmzpzZqN1KktQQ1ilJykdDQkxETKUUYOanlBY1YpuSJEmSVEsj7k4WwGeA7pTSJ+rvkiRJkiQNrRFHYl4FvBn4tYjYWH5c24DtSpIkSdIx6r7FckppDRAN6IskSZIkjaihdyeTJEmSpPFmiJEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSuGGEmSJElZqfsWy5Ik6TmPPnUAgIua3A9JamceiZEkSZKUFUOMJEmSpKwYYiRJkiRlJb8Q89nXlx6SJEmSJqX8QowkSZKkSS27u5N51xdJkiRpcmvKkZje3l42btwIwJEjRyiKgk2bNgHQ399PURRs2bIFgL6+PoqioLu7G4Afp6msevZytm7dCsDBgwcpioJt27YBcODAAYqiYPv27QDs27ePoijYsWMHAHv37qUoCnbu3AnA7t27KYqCnp4eAHbt2kVRFOzatQuAnp4eiqJg9+7dAOzcuZOiKNi7dy8AO3bsoCgK9u3bB8D27dspioIDB0pha9u2bRRFwcGDBwHYunUrRVFw6NAhALq7uymKgr6+PgC2bNlCURT09/cDsGnTJoqi4MiRIwBs3LiRoigqP8uuri7mzZtXaa9fv5758+dX2uvWraOzs7PSXrt2LQsWLKi016xZw8KFCyvt1atXs2jRokp75cqVLF68uNJesWIFS5YsqbSXL1/O0qVLK+1ly5axbNmySvtdc/+Rd859rr9LlixhxYoVlfbixYtZuXJlpb1o0SJWr15daS9cuJA1a9ZU2gsWLGDt2rWVdmdnJ+vWrau058+fz/r16yvtefPm0dXVVWkXRXHcr71Dhw5RFIWvvbJWf+0tXbqU5cuXV9q+9o597Wlo1qnJ81lhnfK1N8A61XqvveF4Opk0Sm/93Hoe+8EPm90NSZpYX/8o7Nrc7F5oFKxTmkwipTThO+3o6EgbNmw4rnUf/dgVAFx055oRllQreMPffwuA+/7gl5vck/q103ORACKiK6XU0ex+tCLrVJWBm+m8Zenwy2WqnT7b2+m5SDB8ncrumph20naFTnlr819UJB0fr0VVy7BOqYqnk0lSht7w99+qfOsqSVKrGe865ZEYjasP9N5e/p9HmzSxPK1C0mhYp9Qs1qn6GGKkUbLQSZJamXVKk4khRhLQfue9W8wlqb1Yp1TNECNJGbL4SZJa2XjXKS/slyRJkpQVQ4wkSZKkrBhiJEmSJGXFECNJkiQpK4YYjZsHHu7h1mfeznU/fB+v+vjXeeDhnmZ3SZIkSW3AEKNx8cDDPdyxaDN70ukkgp79h7lj0eZsg4yBTJLaTzt9trfTc5FGwxCjcTH3wa0c7j9y1LTD/UeY++DWJvXo+LVbIJMktddnezs9F2m0/DsxGhc/2H94TNNb2XCB7IbLzmlSrxpo7vnwo93P/fGwD55e+veUM+H27zerV5I0rtrps72dnktN1inV0JAQExHXAJ8EpgD3ppQ+3ojtKl9nTz+ZnhqB5ezpJzehN/Vpp0BW0492j226mu6Bh3v46DNvZ296AWd//OvcfvUF7fGLSu78RSsr7fTZ3k7PpSbrVHYmok7VHWIiYgrwKeDXgSeB9RHx5ZTSY/Vu+yjtVBza6bkM4V/TW5k2rfeY6X1pBrB94jtUh3YKZJNFO/+SP3DayOFU+twYOG0EaJvnmK12+0WrzWvVhmlvZwb7j5ney3Tg8YnuTl2sU/mxTtWvEUdiXg5sSyltB4iILwDXA40NMe1UHNrpuQxh2o+PDTDDTW9l7RTIJoN2/yW/7U8bUeto81pVK8AMN72VWafyYp1qjEaEmHOAnVXtJ4FXNGC7Uktop0A2GfzK4lfSPWV/6eTWKr2Lp8NleX27WssDh29l5rQDx0zfc/h04ImJ75CkprNO5cU61RgTdneyiLgtIjZExIY9e/ZM1G4lTTLt9O1qLTPj2MIw3HSNnnVK0kSwTjVGI0JMD/DCqvas8rSjpJTuSSl1pJQ6Zs6c2YDdSmqEvpNmjGm61K6sU1Jrsk6plkacTrYeOD8izqUUXt4I/HYDtitpArwm7qWn79gLQs+ZdjL/1oT+SLnqZfqQF4r7q5Z0/KxTqqXuEJNSejYi3gE8SOnsvn9MKT1ad88G6TtpRs1zO/tOmsG0Ru9svJ1yZu0LI085c+L7Ml4mw3NsE21/a05pgnzz+rWli3WrLmg9eeoU/vzGS7ihed06fu3+Od7uz6+NWKdUS0P+TkxK6SvAVxqxraFMu2N76XZ0C1aXbkc3/fn53o6ufGvKRz92BQAX3bmmmb0ZH+30HNu80HlrTqkxBupRW9QpaK/P8Vra6flZpzQJNSTETJQbLjuH87/6twBc9L6MP2yUl3YqdDXcfvUFNb89vv3qC5rYqzq0eTFv++eXOeuUmsI6lZd2/xyfoOeXVYiR1Hh+e5yZdn9+kjSIdSozE/T8DDGS/PZYktTSrFMabML+TowkSZIkNYIhRpIkSVJWDDGSJEmSsmKIkSRJkpQVL+zXuLrorNOb3QVJkoZknZLyZIiRRslCJ0lqZdYpTSaGGEkAfHjGXADua3I/JEmqxTqlaoYYScqQ37hKklrZeNcpQ4zG11uWNrsHmqT8JV+S1MqsU/UxxEijZSDLi+MlaTTa6bOinZ7LZOB41aUpt1ju7e1l48aNABw5coSiKNi0aRMA/f39FEXBli1bAOjr66MoCrq7uwH4cZrKqmcvZ+vWrQAcPHiQoijYtm0bAAcOHKAoCrZv3w7Avn37KIqCHTt2ALB3716KomDnzp0A7N69m6Io6OnpAWDXrl0URcGuXbsA6OnpoSgKdu/eDcDOnTspioK9e/cCsGPHDoqiYN++fQBs376doig4cOAAANu2baMoCg4ePAjA1q1bKYqCQ4cOcdGdazjhv/4DRVHQ19cHwJYtWyiKgv7+fgA2bdpEURQcOXIEgI0bN1IUReVn2dXVxbx58yrt9evXM3/+/Ep73bp1dHZ2Vtpr165lwYIFlfaaNWtYuHBhpb169WoWLVpUaa9cuZLFixdX2itWrGDJkiWV9vLly1m69Lk34bJly1i2bFmlvXTpUpYvX15pL1myhBUrVlTaixcvZuXKlZX2okWLWL16daW9cOFC1qxZU2kvWLCAtWvXVtqdnZ2sW7eu0p4/fz7r16+vtOfNm0dXV1elXRTFcb/2Dh06RFEUbfHaA+ju7j7qtfeCw7t4Ue96X3tlLf/ai5vY9iufBI7/taehWaee+6zo+emZrHr2cutUWXafFRm/9qxTmb/2xrlO+XdiJAHwrqvO58KzX9DsbkhqMS/8medzyvM8cUPNZ51StUgpTfhOOzo60oYNG45r3Uc/dgUAF925ZoQlJUnDiYiulFJHs/vRiqxTVT77+tK/nvoiaYINV6f8akWSJA3N8CKpBXk6mSRJkqSsGGIkSZIkZcUQI0mSJCkr2V0T0zYXSkqSJEk6Lh6JkSRJkpQVQ4wkSZKkrBhiJEmSJGUlu2tiJElqZReddXqzuyBJbc8jMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWakrxETE3Ij4bkRsiogvRcT0BvVLkiRJkmqq90jM14CLU0qXAt8D7qi/S5IkSZI0tLr+2GVKaXlVcx3w3+rrjiRJmXvL0mb3QJLaXiOvifld4KtDzYyI2yJiQ0Rs2LNnTwN3K0lS/axTkpSPEUNMRKyIiC01HtdXLXMX8Cwwf6jtpJTuSSl1pJQ6Zs6c2ZjeS5LUINYpScrHiKeTpZSuGm5+RNwKXAe8JqWUGtQvSZIkSaqprmtiIuIa4L3Ar6aUDjWmS5IkSZI0tHqvifkb4DTgaxGxMSI+3YA+SZIkSdKQ6r072Ysb1RFJkiRJGo1G3p1MkiRJksadIUaSJElSVgwxkiRJkrJiiJEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSuGGEmSJElZMcRIkiRJyoohRpIkSVJWDDGSJEmSsmKIkSRJkpQVQ4wkSZKkrBhiJEmSJGXFECNJkiQpK4YYSZIkSVkxxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVgwxkiRJkrJiiJEkSZKUFUOMJEmSpKwYYiRJkiRlpSEhJiL+KCJSRJzRiO1JkiRJ0lDqDjER8ULgtcAT9XdHkiRJkobXiCMxfwm8F0gN2JYkSZIkDauuEBMR1wM9KaVHGtQfSZIkSRrWiSMtEBErgJ+rMesu4E5Kp5KNKCJuA24DmD179hi6KEnS+LNOSVI+RgwxKaWrak2PiEuAc4FHIgJgFvBQRLw8pbSrxnbuAe4B6Ojo8NQzSVJLsU5JUj5GDDFDSSltBs4caEfEDqAjpbS3Af2SJEmSpJr8OzGSJEmSsnLcR2IGSynNadS2JEmSJGkoHomRJEmSlBVDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVgwxkiRJkrJiiJEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSuGGEmSJElZMcRIkiRJykqklCZ+pxF7gMdrzDodODCKaWcAe8eha6NVq08Tva3Rrjea5YZbZqh5Y5k+2cdrLOuMtOzxzs/lvQWTZ7xa4b31opTSzDrWb1ttUKcgn/eSdar5YzWW9SZyvFpxrKD549XsOjXUvImtUymllnkA94xy2oZW6+dEb2u0641mueGWGWreWKZP9vEayzojLXu883N5b02m8WqH99ZkfPheavw61qnmj1WrjlcrjlUrjFez61SrjFernU62ZJTTmq2RfTrebY12vdEsN9wyQ80b6/RmavZ4jWWdkZY93vm5vLdg8oxXO7y3JiPfS41fxzrV/LEay3oTOV6tOFbQ/PFqdp0aat6EjldTTierV0RsSCl1NLsfGh3HKx+OVV4cr9bl2OTF8cqHY5WX8RyvVjsSM1r3NLsDGhPHKx+OVV4cr9bl2OTF8cqHY5WXcRuvLI/ESJIkSZq8cj0SI0mSJGmSMsRIkiRJyoohRpIkSVJWDDGSJEmSsmKIkSRJkpQVQ4wkSZKkrBhiJEmSJGXFECNJkiQpK4YYSZIkSVkxxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUbZiYgrI+LJiV53lNtfFRFvHWLe7Ig4GBFTxmv/w/RrR0RcNdH7lSTVNp61LCJSRLx4iHm3RMTy49lvPSJiTrlfJ070vtWeDDEaUvkX38PlX7wHHn8zivXGNSiMl6E+YCOiiIg/q3f7KaUnUkqnppSO1LutRio/vxQRL6+a9uKISFXtmuHMoiSp1U22WjaSlNL8lNJrm92PwcrjtDsiTqma9taIWFXVrhnOIuLWiFgzQV1VizDEaCS/Uf7Fe+DxjmZ3SOPiaaDuoCZJLcpalocpwLua3QnlwRCj4xIRfxcR91e1746Ify1/g/JV4Oyqb7zOjogTIuJ9EfF/IqI3IhZExM+W1x34Nv93IuKJiNgbEXdVbfvk8tGCfRHxGPBLg/pydkTcHxF7IuLfI+Kdo133OJ73rRGxJiL+orzNf4+I1w1a7Ocj4jsR8cOIWFzjeZ5Ybq+KiI9ExL9FxDMRsTwizhjlz2TIn2d5/psj4vHyvLsY2eeASyPiV+v5+UhSTtq8ll0bEdvL/ZgbESeUt3XUUYtyn98WEd+PiP0R8amIiOplh6p5EXF6RHwmIp6KiJ6I+LMonzIdEVPK6+2NiO3A60fR57nAeyJi+iiW1SRniNHx+iPgkvIH3K8Avwf8TkrpR8DrgB9UfeP1A+B/ATcAvwqcDewDPjVom1cAFwCvAT4QEb9Ynv6nwM+XH1cDvzOwQvlDeQnwCHBOed0/jIirR1q3Dq8AtgJnAP8b+MzAB37Zfwd+FzgLeBb462G29dvAW4AzgecB7xk0f6ifyZA/z4i4EPg74M3leTOAWSM8p0PAx4CPjrCcJLWTdq5l/xXoAF4GXE+pLg3lOkrB6FLg5vI+BgxX8wpKde7FwGXAa4GBU49/v7zdy8r9+G+j6PMGYBXH1kLpWCklHz5qPoAdwEFgf9Xj96vmv4LSaUiPA79VNf1K4MlB2+oGXlPVPgvoB04E5gAJmFU1/zvAG8v/3w5cUzXvtoHtl/vwxKB93QF8dqR1azzfgX6cOGh6AfxZ+f+3Atuq5j2/vM7PldurgI9Xzb8Q+AmlQ+RHbb+87Purln07sGxQX4b6mQz38/wA8IWqeaeU+3DVEM+7oHQq2UnAE5QK94tLHw+VZVYBbx3tz8yHDx8+WuUx2WpZeX4atPzbgX8t//9WYM2gZa+oai8A3le1bM2aB/wn4MfAyVXzfwtYWf7/14G3Vc177XD1ojxOVwEXAweAmZQC0apBfX1xjXWPek4+JsfDi3E1khtSSitqzUgpfbt8iPhMSh96w3kR8KWI+GnVtCOUPgQH7Kr6/yHg1PL/zwZ2Vs17fNB2z46I/VXTpgDfHMW6gz1b/ndq1f8H2v21+plSOlT+QurUqvmD9zeV0jdYtQz1nEeaP9zP86jnnFL6UUT0DrF/qpb7cUR8BPgI8MaRlpekjEymWjZg8PJnD7PscLVoqJr3s5Tq21NVJyOcULXf4+kzKaUtEfEvwPsohUapJk8n03GLiP9J6dv7HwDvrZqVaiy+E3hdSml61WNaSqlnFLt6CnhhVXv2oO3++6DtnpZSunYU69baTz+lb9OqncsoP3zLBu+vH9g7hvVHY7if51HPOSKeT+mUstH4LDAduLHB/ZWkltSGtWzA4OV/MIp1xmInpSMxZ1T1+QUppYvK84+nzwP+lNLpaOc0pqtqR4YYHZeI+AVKpyC9idK1F++NiJeWZ/8HMCMiTq9a5dPARyPiReX1Z0bE9aPc3QLgjoj4mYiYRemc5AHfAZ6JiD8uX/g4JSIujohfGsW6R0mlWx/fX+7njIiYGhG/RemUsK+Osq8Ab4qIC8vh4cPAwtT42yoP9/NcCFwXEVdExPPKfRjVez2l9Cyl4vHHNWafGBHTqh5Tq+adNGieny2SWl471rIqt5eXfyGlO37dN8p+jkpK6SlgOfD/RsQLonTTg5+P524QswB4Z0TMioifoXRkZbTb3lbu7ztrzH7eoHoz8LfXYtD0aXU8PWXAXzQ0kiVx9L31vxSlu2v9E3B3SumRlNL3gTuBz0fESSml7wKdwPYo3enkbOCTwJeB5RHxDLCO0jnAo/EhSkdC/p3SB+bnB2aUw8F1wEvL8/cC9wKnj7TuEN5O6dzoTcBu4B3A61NK/zHKvlLeR0HpEPw0an8I12vIn2dK6VHgfwL/TOmbsH3AWP7WQWd5vcH+Djhc9fhs1byDg+b92hj2J0njbbLVMoDFQBewEVgKfGaU/RyL/07ppjSPUao1CyldJwTwD8CDlG5W8BCwaIzb/jClazoHe5Sj681bytNfOWj64fBvmLW1SKnW0VJJkiRJak0eiZEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSuGGEmSJElZacqt584444w0Z86cZuxaklTW1dW1N6U0s9n9aEXWKUlqvuHqVFNCzJw5c9iwYUMzdi1JKouIx5vdh1ZlnZKk5huuTnk6mSRJkqSsGGIkSZIkZcUQI0mSJCkrTbkmRpIkSZrM+vv7efLJJ+nr62t2V5pu2rRpzJo1i6lTp456HUOMJEmSNMGefPJJTjvtNObMmUNENLs7TZNSore3lyeffJJzzz131Ot5OpkkSZI0wfr6+pgxY8akDjAAEcGMGTPGfETKECNJkiQ1wWQPMAOO5+dgiJEkSZKUFUOMJEmSpKwYYiRJkqQMvOHvv8Ub/v5b47qPN7/5zUTEqB8ADzzwABHBd7/73XHtWzVDjCRJkiSeeuopLrjgAlJKo34AdHZ20tHRQWdn54T11RAjSZIktbgHHu7h4Sf28+1/f5pXffzrPPBwT8P3cd9993HLLbdU2jfeeCPvf//7efWrX83s2bNZsWLFMescPHiQVatWce+99x4TYq688srK0Zne3l4uvvjihvXVECNJkiS1sAce7uGORZv5yZGfAtCz/zB3LNpcd5Dp7e3lwQcfrLQff/zxo/5Wy+bNm5k+fTrf+MY3+OQnP8n8+fOP2cbixYu56qqreMlLXsKpp55KV1dXZd62bdv4hV/4BQA2bdrEJZdcUld/qxliJEmSpBY298GtHO4/ctS0w/1HmPvg1rq2u2PHDt7znvewf/9+uru7jwoZhw4d4sCBA7z73e8GoL+/n+nTpx+zjc7OTm6++WYAbr755srRmMcff5xzzjmHE04oxY1NmzZx6aWX1tXfag0LMRExJSIejoh/adQ2JUmSpMnuB/sPj2n6aF1++eXcdNNNfPGLX+T+++/npptuqsx77LHHuPzyy5kyZQpQCiGDTwd7+umn+fa3v80111wDlELMfffdR0qJRx555KjQ0tXV1ZohBngX0N3A7UmSJEmT3tnTTx7T9LF405vexOc//3kOHTrEaaedVpm+efNmXvrSl1batY6kLFy4kGuvvZaTTjoJgPPOO4+zzjqLb37zm2zcuJG+vj4Avv/977N48eLWO50sImYBrwfubcT2JEmSJJXcfvUFnDx1ylHTTp46hduvvqDubZ933nkcOXKEK6+88qjpg0PMli1bjjkS09nZyZIlS5gzZ07l0d3dTWdnJ4888gg//elPeclLXsKHP/xhLrzwQj73uc/V3d8BJzZoO38FvBc4bYTlJEmSJI3BDZedA8B7F27iJ0d+yjnTT+b2qy+oTK/X3XffzStf+cqjpn3iE584qr19+/Zj1lu5cuWQ2zz//PN56KGHjjq600h1h5iIuA7YnVLqiogrh1nuNuA2gNmzZ9e7W0mSGso6JamV3XDZOXR+5wkA7vuDX27otq+44oqGbu+ZZ54hIsYtwEBjjsS8CvjNiLgWmAa8ICL+KaX0puqFUkr3APcAdHR0pAbsV5KkhrFOSWp1jQ4v4+W0007je9/73rjuo+5rYlJKd6SUZqWU5gBvBL4+OMBIkiRJUqP4d2IkSZIkZaVRF/YDkFJaBaxq5DYlSZIkqZpHYiRJkiRlxRAjSZIkKSuGGEmSJElZMcRIkiRJyoohRpIkSVJWDDGSJElSDj77+tJDhhhJkiRJJW9+85uJiFE/AB544AEigu9+97sT1s+G/p0YSZIkSQ0293z40e7n2h88vfTvKWfC7d9v2G6eeuopLrjgAlJKY1qvs7OTjo4OOjs7+dCHPtSw/gzHIzGSJElSK6sOMKOZfpzuu+8+brnllkr7xhtv5P3vfz+vfvWrmT17NitWrDhmnYMHD7Jq1SruvfdeOjs7j5r3yCOP8OpXv5oLL7yQE044gYjgAx/4QEP6aoiRJEmSJqHe3l4efPDBSvvxxx/n3HPPrbQ3b97M9OnT+cY3vsEnP/lJ5s+ff8w2Fi9ezFVXXcVLXvISTj31VLq6ugDo6+vjDW94A3/xF3/BY489xl133cV73vOehh2pMcRIkiRJk9COHTt4z3vew/79++nu7uaSSy6pzDt06BAHDhzg3e9+NwD9/f1Mnz79mG10dnZy8803A3DzzTdXjsasWLGCl73sZbz85S8H4NJLL+Xpp5+uXEdTL0OMJEmSNAldfvnl3HTTTXzxi1/k/vvv56abbqrMe+yxx7j88suZMmUKAJs2beLiiy8+av2nn36ab3/721xzzTVAKcTcd999pJTYsmXLUaHooYce4mUve1nD+m6IkSRJkiapN73pTXz+85/n0KFDnHbaaZXpmzdv5qUvfWmlvWnTJi699NKj1l24cCHXXnstJ510EgDnnXceZ511Ft/85jeZMWMGmzZtAuB73/seixYt4o1vfGPD+u3dySRJkqRWdsqZtS/iP+XMujd93nnnceTIEa688sqjpm/evJlXvOIVlfaWLVuOORLT2dnJI488wpw5cyrTent76ezsZO7cuXz5y1/m4osv5owzzqCzs5MZM2bU3d8BhhhJkiSplQ3cRnngD12+ZWlDN3/33Xfzyle+8qhpn/jEJ45qb9++/Zj1Vq5cOex2lyxZUn/nhmCIkSRJknLQ4PAy4IorrhiX7Y4nr4mRJEmSlBVDjCRJkqSsGGIkSZKkJkgpNbsLLeF4fg6GGEmSJGmCTZs2jd7e3kkfZFJK9Pb2Mm3atDGt54X9kiRJ0gSbNWsWTz75JHv27Gl2V5pu2rRpzJo1a0zrGGIkSZKkCTZ16lTOPffcZncjW55OJkmSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVgwxkiRJkrJiiJEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSt1h5iIeGFErIyIxyLi0Yh4VyM6JkmSJEm1nNiAbTwL/FFK6aGIOA3oioivpZQea8C2JUmSJOkodR+JSSk9lVJ6qPz/Z4Bu4Jx6tytJkiRJtTT0mpiImANcBny7xrzbImJDRGzYs2dPI3crSVLdGlanPvv60kOSNG4aFmIi4lTgfuAPU0o/HDw/pXRPSqkjpdQxc+bMRu1WkqSGsE5JUj4aEmIiYiqlADM/pbSoEduUJEmSpFoacXeyAD4DdKeUPlF/lyRJkiRpaI04EvMq4M3Ar0XExvLj2gZsV5IkSZKOUfctllNKa4BoQF8kSZIkaUQNvTuZJEmSJI03Q4wkSZKkrBhiJEmSJGXFECNJkiQpK4YYSZIkSVkxxEiSJEnKiiFGkiRJUlbq/jsxkiTpOY8+dQCAi5rcD0lqZx6JkSRJkpQVQ4wkSZKkrOQXYj77+tJDkiRJ0qSUX4iRJEmSNKlld2G/F0xKkiRJk1tTjsT09vayceNGAI4cOUJRFGzatAmA/v5+iqJgy5YtAPT19VEUBd3d3QD8OE1l1bOXs3XrVgAOHjxIURRs27YNgAMHDlAUBdu3bwdg3759FEXBjh07ANi7dy9FUbBz504Adu/eTVEU9PT0ALBr1y6KomDXrl0A9PT0UBQFu3fvBmDnzp0URcHevXsB2LFjB0VRsG/fPgC2b99OURQcOFAKW9u2baMoCg4ePAjA1q1bKYqCQ4cOAdDd3U1RFPT19QGwZcsWiqKgv78fgE2bNlEUBUeOHAFg48aNFEVR+Vl2dXUxb968Snv9+vXMnz+/0l63bh2dnZ2V9tq1a1mwYEGlvWbNGhYuXFhpr169mkWLFlXaK1euZPHixZX2ihUrWLJkSaW9fPlyli5dWmkvW7aMZcuWVdrvmvuPvHPuc/1dsmQJK1asqLQXL17MypUrK+1FixaxevXqSnvhwoWsWbOm0l6wYAFr166ttDs7O1m3bl2lPX/+fNavX19pz5s3j66urkq7KIrjfu0dOnSIoih87ZW1+mtv6dKlLF++vNL2tXfsa09Ds05Nns8K65SvvQHWqdZ77Q3H08mkUXrr59bz2A9+2OxuSNLE+vpHYdfmZvdCo2Cd0mQSKaUJ32lHR0fasGHDca376MeuAOCiO9eMsKRawRv+/lsA3PcHv9zkntSvnZ6LBBARXSmljmb3oxVZp6oM3EznLUuHXy5T7fTZ3k7PRYLh61R218S0lTYvDMqMr0dJNXgtqlqGdUpVDDFNZGGQdLz8xlWS1MrGu04ZYjSuPtB7e/l/bXJahbLhL/mSRsM6pWaxTtXHECONkoVOktTKrFOaTAwxktqSxVyS1MqsU/UxxEgCvEYrNxY/SZONdSov412n/DsxkiRJkrJiiJEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSuGGI2bBx7u4dZn3s51P3wfr/r413ng4Z5md0mSpArrlJQvQ4zGxQMP93DHos3sSaeTCHr2H+aORZuzLRAWOklqL9YpKW/+nRiNi7kPbuVw/5Gjph3uP8LcB7dyw2XnNKlXx2eg0B1OpwNUCh2Q3XOpae758KPdz913/4Ol58kpZ8Lt329WryRpXFmnMmKdUg0NORITEddExNaI2BYR72vENpW3H+w/PKbprWy4QtcWfrR7bNPVdH7j2qLmng8fPJ2LfrKZi36yufSL1gdPL01Xy7FOZcQ6lZ2JqFN1H4mJiCnAp4BfB54E1kfEl1NKj9W77Xb2wMM9fPSZt7M3vYCzP/51br/6gvb4tqRsw7S3M4P9x0zvZTrw+ER3py7tVOgmi3Z+f7X9N645a7dftNr822/rlJrJOlW/RpxO9nJgW0ppO0BEfAG4HmhsiGmjD9O+Pz+PG37cyw0nDUwAFkPfshlMu2N7M7vWMLUKw3DTW1k7FbrJoN1/yW+nU2DU4totlA1inVKzWKcaoxEh5hxgZ1X7SeAVDdju0drow3Taj3vHNF3N1U6FbjL4lcWvpHvKfphy9PTexdPhsvyLud+4ShrMOpUX61RjTNiF/RFxG3AbwOzZsydqt5ImmXYv5n7jOn6sU5ImgnWqMRpxYX8P8MKq9qzytKOklO5JKXWklDpmzpzZgN1KaoS+k2aMabqaq92LXzNZp6TWZJ3Ky0TVqUYciVkPnB8R51IKL28EfrsB25U0AV4T99LTd+wh3nOmncy/NaE/Uq56mT7kt4/+qiUdP+uUaqk7xKSUno2IdwAPUjq77x9TSo/W3TPl7ZQza1+vdMqZE98XDctrLKTG+Ob1a0sX61Zd0Hry1Cn8+Y2XcEPzunXc+k6aUfNazb6TZjCtCf1pOOtUNqxTqqUh18SklL4CfKUR2xpKW32YToYPzvId4x792BUAXHTnmmb2pj5tPl5nTz+ZnhqF4OzpJzehN1K+Bu6689EFq0u3TZ3+/Kxvmzrtju2l28C2yfM5hnUqG9Yp1TJhF/bXq60+TNvpg3MyaPPxuv3qC2p+e3z71Rc0sVd1aPNirtZ2w2XncP5X/xaAi96X/2dFuz2ftmWdyot1qiGyCTHgh6k0Htrt2+N2L+YWP0mTjXUqMxNUp7IKMZLGh18QZKTdi58k1WCdysgE1alG3GJZkiRJkiaMIUaSJElSVgwxkiRJkrJiiJEkSZKUFS/s17i66KzTm90FSZIktRlDjDRKBjJJaj/t9NneTs9FGokhRhIAH54xF4D7mtwPSZJqsU6pmtfESJIkScqKR2I0vt6ytNk90CTV7qdVtPvzk6R21+6f4+P9/Awx0mgZyCSp/bTTZ3s7PRdpBIaYJmr3BC41VbsX83Z/fpLU7tr9c3ycn19Tronp7e1l48aNABw5coSiKNi0aRMA/f39FEXBli1bAOjr66MoCrq7uwH4cZrKqmcvZ+vWrQAcPHiQoijYtm0bAAcOHKAoCrZv3w7Avn37KIqCHTt2ALB3716KomDnzp0A7N69m6Io6OnpAWDXrl0URcGuXbsA6OnpoSgKdu/eDcDOnTspioK9e/cCsGPHDoqiYN++fQBs376doig4cOAAANu2baMoCg4ePAjA1q1bKYqCQ4cOwVuW0v1f/oKiKOjr6wNgy5YtFEVBf38/AJs2baIoCo4cOQLAxo0bKYqi8rPs6upi3rx5lfb69euZP39+pb1u3To6Ozsr7bVr17JgwYJKe82aNSxcuLDSXr16NYsWLaq0V65cyeLFiyvtFStWsGTJkkp7+fLlLF363It02bJlLFu2rNJeunQpy5cvr7SXLFnCihUrKu3FixezcuXKSnvRokWsXr260l64cCFr1qyptBcsWMDatWsr7c7OTtatW1dpz58/n/Xr11fa8+bNo6urq9IuiuK4X3uHDh2iKIr2eO0B3d3dR732/uRVp/G6k7b62iubDK89Dc069dxnxQmnz2H9Ca+wTpVNxs8K61SJr73WqlNe2C9JkiQpK5FSmvCddnR0pA0bNhzXuo9+7AoALrpzzQhLSpKGExFdKaWOZvejFVmnJKn5hqtTHomRJEmSlBVDjCRJkqSsGGIkSZIkZSW7Wyx7jrEkSZI0uXkkRpIkSVJWDDGSJEmSsmKIkSRJkpSV7K6JkSSplV101unN7oIktT2PxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZaWuEBMRcyPiuxGxKSK+FBHTG9QvSZIkSaqp3iMxXwMuTildCnwPuKP+LkmSJEnS0OoKMSml5SmlZ8vNdcCs+rskSZIkSUM7sYHb+l3gvqFmRsRtwG0As2fPbuBuJUmqX8Pq1FuWNqhHkqShjHgkJiJWRMSWGo/rq5a5C3gWmD/UdlJK96SUOlJKHTNnzmxM7yVJahDrlCTlY8QjMSmlq4abHxG3AtcBr0kppQb1S5IkSZJqqut0soi4Bngv8KsppUON6ZIkSZIkDa3eu5P9DXAa8LWI2BgRn25AnyRJkiRpSHUdiUkpvbhRHZEkSZKk0aj3SIwkSZIkTShDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVgwxkiRJkrJiiJEkSZKUFUOMJEmSpKwYYiRJkiRlxRAjSZIkKSuGGEmSJElZMcRIkiRJyoohRpIkSVJWDDGSJEmSsmKIkSRJkpQVQ4wkSZKkrBhiJEmSJGXFECNJkiQpK4YYSZIkSVkxxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScpKQ0JMRPxRRKSIOKMR25MkSZKkodQdYiLihcBrgSfq744kSZIkDa8RR2L+EngvkBqwLUmSJEkaVl0hJiKuB3pSSo80qD+SJEmSNKwTR1ogIlYAP1dj1l3AnZROJRtRRNwG3AYwe/bsMXRRkqTxZ52SpHyMGGJSSlfVmh4RlwDnAo9EBMAs4KGIeHlKaVeN7dwD3APQ0dHhqWeSpJZinZKkfIwYYoaSUtoMnDnQjogdQEdKaW8D+iVJkiRJNfl3YiRJkiRl5biPxAyWUprTqG1JkiRJ0lA8EiNJkiQpK4YYSZIkSVkxxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlBVDjCRJkqSsGGIkSZIkZcUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVgwxkiRJkrJiiJEkSZKUFUOMJEmSpKxESmnidxqxB3i8xqzTgQOjmHYGsHccujZatfo00dsa7XqjWW64ZYaaN5bpk328xrLOSMse7/xc3lswecarFd5bL0opzaxj/bbVBnUK8nkvWaeaP1ZjWW8ix6sVxwqaP17NrlNDzZvYOpVSapkHcM8op21otX5O9LZGu95olhtumaHmjWX6ZB+vsawz0rLHOz+X99ZkGq92eG9NxofvpcavY51q/li16ni14li1wng1u061yni12ulkS0Y5rdka2afj3dZo1xvNcsMtM9S8sU5vpmaP11jWGWnZ452fy3sLJs94tcN7azLyvdT4daxTzR+rsaw3kePVimMFzR+vZtepoeZN6Hg15XSyekXEhpRSR7P7odFxvPLhWOXF8Wpdjk1eHK98OFZ5Gc/xarUjMaN1T7M7oDFxvPLhWOXF8Wpdjk1eHK98OFZ5GbfxyvJIjCRJkqTJK9cjMZIkSZImKUOMJEmSpKwYYiRJkiRlpS1CTEScFxGfiYiFze6LhhcRN0TEP0TEfRHx2mb3R8OLiF+MiE9HxMKI+B/N7o+GFxGnRMSGiLiu2X3R0axTebFW5cM6lZ9G1aqWDTER8Y8RsTsitgyafk1EbI2IbRHxPoCU0vaU0u81p6ca41g9kFL6feBtwBua0d/Jbozj1Z1SehtwM/CqZvR3MhvLWJX9MbBgYns5eVmn8mKtyod1Ki/NqlUtG2KAArimekJETAE+BbwOuBD4rYi4cOK7pkEKxj5W7y/P18QrGMN4RcRvAkuBr0xsN8UYxioifh14DNg90Z2cxAqsUzkpsFblosA6lZOCJtSqlg0xKaVvAE8PmvxyYFv5G62fAF8Arp/wzukoYxmrKLkb+GpK6aGJ7qvG/t5KKX05pfQ64JaJ7anGOFZXAv8F+G3g9yOiZT/f24V1Ki/WqnxYp/LSrFp14vGu2CTnADur2k8Cr4iIGcBHgcsi4o6U0p83pXeqVnOsgP8FXAWcHhEvTil9uhmd0zGGem9dCdwInITfcLWKmmOVUnoHQETcCuxNKf20CX2TdSo31qp8WKfyMu61KrcQU1NKqZfSeatqcSmlvwb+utn90OiklFYBq5rcDY1BSqlodh90LOtUXqxV+bBO5akRtSq30w16gBdWtWeVp6n1OFZ5cbzy4Vi1NscnL45XPhyrvIz7eOUWYtYD50fEuRHxPOCNwJeb3CfV5ljlxfHKh2PV2hyfvDhe+XCs8jLu49WyISYiOoFvARdExJMR8XsppWeBdwAPAt3AgpTSo83spxyr3Dhe+XCsWpvjkxfHKx+OVV6aNV6RUmrk9iRJkiRpXLXskRhJkiRJqsUQI0mSJCkrhhhJkiRJWTHESJIkScqKIUaSJElSVgwxkiRJkrJiiJEkSZKUFUOMNISIWBQRfxYR34iIJyLiqmb3SZKkAdYpTWaGGGlolwD7U0qvBt4F3NLk/kiSVM06pUnLECPVEBHPB04H/rI8aSqwv8Zyt0bEdRPYNUmSrFOa9E5sdgekFnUh0JVSOlJuXwpsiYhbgf8LeAw4Avxn4PkRAfAC4ErgGeDOlNKPJ7jPkqTJwzqlSc0jMVJtlwAbq9qXApvK//9aSuluoANYC/xzSulfgJ8vL/NXFgZJ0jizTmlSM8RItQ0uDhcDW8r/HziCORV4dmCBlNJHgG8CcyPi/AnooyRp8rJOaVLzdDKphpTS/zOofR5A+XD8ayPiUuA7wCPAXRFxInAmcD7wU6B3QjssSZpUrFOa7CKl1Ow+SNkon2u8t3xYXpKklmKd0mRhiJEkSZKUFa+JkSRJkpQVQ4wkSZKkrBhiJEmSJGXFECNJkiQpK4YYSZIkSVkxxEiSJEnKiiFGkiRJUlYMMZIkSZKyYoiRJEmSlJX/HxtIxs2S/r4hAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1008x576 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2, 2, figsize=(14, 8), sharex=True, sharey=True)\n",
    "\n",
    "plt.sca(ax[0, 0])\n",
    "plt.title(\"Unbinned NLL\")\n",
    "plt.errorbar(\n",
    "    n_pts,\n",
    "    n_pts ** 0.5 * unbinned_nll[0],\n",
    "    np.sqrt(n_pts * unbinned_nll[1]),\n",
    "    fmt=\"o\",\n",
    "    label=r\"$\\sqrt{n}\\,\\Delta\\mu$\",\n",
    ")\n",
    "plt.errorbar(\n",
    "    n_pts,\n",
    "    n_pts ** 0.5 * unbinned_nll[2],\n",
    "    np.sqrt(n_pts * unbinned_nll[3]),\n",
    "    fmt=\"s\",\n",
    "    label=r\"$\\sqrt{n}\\,\\Delta\\sigma$\",\n",
    ")\n",
    "\n",
    "plt.sca(ax[0, 1])\n",
    "plt.title(\"Binned NLL\")\n",
    "plt.errorbar(\n",
    "    n_pts,\n",
    "    n_pts ** 0.5 * binned_nll[0],\n",
    "    np.sqrt(n_pts * binned_nll[1]),\n",
    "    fmt=\"o\",\n",
    "    label=r\"$\\sqrt{n}\\,\\Delta\\mu$\",\n",
    ")\n",
    "plt.errorbar(\n",
    "    n_pts,\n",
    "    n_pts ** 0.5 * binned_nll[2],\n",
    "    np.sqrt(n_pts * binned_nll[3]),\n",
    "    fmt=\"s\",\n",
    "    label=r\"$\\sqrt{n}\\,\\Delta\\sigma$\",\n",
    ")\n",
    "\n",
    "plt.sca(ax[1, 0])\n",
    "plt.title(\"Extended Unbinned NLL\")\n",
    "plt.errorbar(\n",
    "    n_pts,\n",
    "    n_pts ** 0.5 * extended_unbinned_nll[2],\n",
    "    np.sqrt(n_pts * extended_unbinned_nll[3]),\n",
    "    fmt=\"o\",\n",
    "    label=r\"$\\sqrt{n}\\,\\Delta\\mu$\",\n",
    ")\n",
    "plt.errorbar(\n",
    "    n_pts,\n",
    "    n_pts ** 0.5 * extended_unbinned_nll[4],\n",
    "    np.sqrt(n_pts * extended_unbinned_nll[5]),\n",
    "    fmt=\"s\",\n",
    "    label=r\"$\\sqrt{n}\\,\\Delta\\sigma$\",\n",
    ")\n",
    "\n",
    "plt.sca(ax[1, 1])\n",
    "plt.title(\"Extended binned NLL\")\n",
    "plt.errorbar(\n",
    "    n_pts,\n",
    "    n_pts ** 0.5 * extended_binned_nll[2],\n",
    "    np.sqrt(n_pts * extended_binned_nll[3]),\n",
    "    fmt=\"o\",\n",
    "    label=r\"$\\sqrt{n}\\,\\Delta\\mu$\",\n",
    ")\n",
    "plt.errorbar(\n",
    "    n_pts,\n",
    "    n_pts ** 0.5 * extended_binned_nll[4],\n",
    "    np.sqrt(n_pts * extended_binned_nll[5]),\n",
    "    fmt=\"s\",\n",
    "    label=r\"$\\sqrt{n}\\,\\Delta\\sigma$\",\n",
    ")\n",
    "\n",
    "plt.ylim(-5, 5)\n",
    "plt.legend()\n",
    "plt.semilogx();\n",
    "for i in (0, 1):\n",
    "    ax[1, i].set_xlabel(r\"$n_\\mathrm{pts}$\")\n",
    "for axi in ax.flat:\n",
    "    for y in (-1, 1):\n",
    "        axi.axhline(y, ls=\":\", color=\"0.5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Least-squares fits\n",
    "\n",
    "We do the same as before, but this time we use a least-squares fit of $x,y$ scattered data and vary the residual function. Other functions than the identity can be used to reduce the pull of large outliers, turning the ordinary least-squares fit into a robust fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 17 need to re-try [(True, True), (True, True), (False, False)]\n",
      "10000 69 need to re-try [(True, True), (True, True), (False, False)]\n"
     ]
    }
   ],
   "source": [
    "n_tries = 100  # increase this to 500 to get less scattering\n",
    "\n",
    "truth = Namespace(a=1, b=2)\n",
    "\n",
    "n_pts = np.array((10, 30, 100, 300, 1000, 3000, 10000))\n",
    "\n",
    "\n",
    "@joblib.delayed\n",
    "def compute(n):\n",
    "    rng = np.random.default_rng(n)\n",
    "    x = np.linspace(0, 1, n)\n",
    "\n",
    "    linear = []\n",
    "    soft_l1 = []\n",
    "    arctan = []\n",
    "    for i_try in range(n_tries):\n",
    "\n",
    "        def model(x, a, b):\n",
    "            return a + b * x\n",
    "\n",
    "        while True:\n",
    "            y = model(x, 1, 2)\n",
    "            ye = 0.1\n",
    "            y += rng.normal(0, ye, len(y))\n",
    "\n",
    "            m = [\n",
    "                Minuit(LeastSquares(x, y, ye, model), a=0, b=0),\n",
    "                Minuit(LeastSquares(x, y, ye, model, loss=\"soft_l1\"), a=0, b=0),\n",
    "                Minuit(LeastSquares(x, y, ye, model, loss=np.arctan), a=0, b=0),\n",
    "            ]\n",
    "\n",
    "            all_good = True\n",
    "            for mi in m:\n",
    "                mi.migrad()\n",
    "                mi.hesse()\n",
    "                if not mi.valid or not mi.accurate:\n",
    "                    all_good = False\n",
    "                    break\n",
    "            if all_good:\n",
    "                break\n",
    "            print(f\"{n} {i_try} need to re-try {[(mi.valid, mi.accurate) for mi in m]}\")\n",
    "\n",
    "        linear.append(\n",
    "            (\n",
    "                m[0].values[\"a\"] - truth.a,\n",
    "                m[0].values[\"b\"] - truth.b,\n",
    "                m[0].errors[\"a\"] ** 2,\n",
    "                m[0].errors[\"b\"] ** 2,\n",
    "            )\n",
    "        )\n",
    "        soft_l1.append(\n",
    "            (\n",
    "                m[1].values[\"a\"] - truth.a,\n",
    "                m[1].values[\"b\"] - truth.b,\n",
    "                m[1].errors[\"a\"] ** 2,\n",
    "                m[1].errors[\"b\"] ** 2,\n",
    "            )\n",
    "        )\n",
    "        arctan.append(\n",
    "            (\n",
    "                m[2].values[\"a\"] - truth.a,\n",
    "                m[2].values[\"b\"] - truth.b,\n",
    "                m[2].errors[\"a\"] ** 2,\n",
    "                m[2].errors[\"b\"] ** 2,\n",
    "            )\n",
    "        )\n",
    "\n",
    "    return [\n",
    "        (*np.mean(t, axis=0), *np.var(np.array(t)[:,:2], axis=0))\n",
    "        for t in (linear, soft_l1, arctan)\n",
    "    ]\n",
    "\n",
    "linear = []\n",
    "soft_l1 = []\n",
    "arctan = []\n",
    "\n",
    "for l, s, a in joblib.Parallel(-1)(compute(n) for n in n_pts):\n",
    "    linear.append(l)\n",
    "    soft_l1.append(s)\n",
    "    arctan.append(a)\n",
    "\n",
    "linear = np.transpose(linear)\n",
    "soft_l1 = np.transpose(soft_l1)\n",
    "arctan = np.transpose(arctan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0kAAAITCAYAAADSG1dlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABVA0lEQVR4nO3df7xcVXno/89DSAgQCJZgBQICFbGAASUCYq5FAYsoWFHiD/Qq1VKtXNRbr0VBFNT6436rV6oVkWrU0kgE0aAoAhWQqkiwIeGXSik0AVQSJCEmkkCe7x97TjLnnDknc3L2zN5z5vN+veZ1Zq+9Z806a/Y8M8+svdeOzESSJEmSVNim6gZIkiRJUp2YJEmSJElSE5MkSZIkSWpikiRJkiRJTUySJEmSJKmJSZIkSZIkNTFJkiT1jYj4XkS8aZT18yLiIx1uQ0TElyPidxHxs04+l6TOqkNMUWeYJPWpiLgvIo7twvMcHRHLO/08Ut116z3XeC7fdyPIzJdm5lcAIuLNEXFTBc2YAxwHzMzMwyPiQxHxL6M9YKT9JyKmRMRljfUZEUd3psmqG2NKPdQkpgxTp7b0KpMkCYiIbatug9Rv+vh993Tgvsz8fUn13QS8Afh1SfVtUR+/dqqxftsvI2JS1W0Yr8bIei3zkVo2StWIiG0i4qyI+M+IWBkRCyLij5rWfyMifh0RqyLixog4qGndCRFxZ0Q8FhEPRMR7ImJH4HvAHhGxpnHbo8XzTo2If2k856MRcUtE/HFj3b4RcUOj3msi4rMDv7i2+mWr+Ze1iDg8In7SqPOhxmOnNG2bEfGOiPgV8KtG2csjYnHjMT+OiFlN2/9d4397LCJ+ERHHlNT16mO+78b/vmu099GBD9qI+GJE/LZp/dci4l2N+9dHxFsj4k+BC4HnN/ro0aYqnxIR3208580R8ScjvHaj9eEeEbEwIh6JiHsi4q8a5W8BLm563puB9wOvaSzf1uq5RpKZ6zPz/2XmTcCTW9q+8f9/OCL+vfH//SAiZjStPyki7mj8P9c3+mlg3X2N12MJ8PuIeEbj9TwtIpZFcfjg2yLieRGxpFHHZ8fy/2j8jCm9G1PaeH3mRcTnI+KqiPg98KKI2CsivhkRDzf6/rMjtSUiXhYR/xERqxvv2Q811b1Poy/fFBH/HRErIuLsUdo5LyI+N9L/FRFHNfaBVY2/RzWtuz4iPhoR/w6sBfZrPPffRMSvGvV9OCL+pPH6rW7sx1NataVjMtNbH96A+4Bjh5S9E/gpMBPYDvgCML9p/V8COzXW/T9gcdO6h4D/0bj/FOC5jftHA8u30Ja/Bq4EdgAmAYcBOzfW/QT4VOM5Xwg8BvzLSHU3/1+Neo4EtgX2Ae4C3tW0bQLXAH8EbA88B/gtcESjHW9q1LcdcACwDNij8dh9gD+p+nX01ju3Vu+5RrnvuxLed8B/A4c17v8CuBf406Z1z2ncvx54a+P+m4GbhtQzD1gJHN74Hy4Bvr4VfXgj8E/AVOBQ4GHgxa2eF/jQQP+Odf8Zss1y4OgtbHM98J/AMxv9fz3w8ca6ZwK/pzgUcDLwXuAeYEpTGxYDezUeu0/j9byw8X++BPgD8C3gqcCejdf2z6p+/03E20j7BMaUno0pbbw+84BVwAsoBjp2BG4DPt24PxWYM0pbjgae3XjsLOA3wF809UUCX2z04yHA4wP/c4t2jvh/NV6L3wFvbKx7XWN516Y++2/goMb6yY3n/jawc6P8ceA6YD9gOnAn8KZuvsccSVKztwFnZ+byzHyc4oP71dEYvs7ML2XmY03rDomI6Y3HbgAOjIidM/N3mfnzMTzvBmBX4BmZ+WRm3pqZqyNib+B5wAcy8/HMvJEiALelUc9PM/OJzLyP4oPiz4Zs9rHMfCQz1wGnA1/IzJsb7fgKxZv0SIpfaLdr/I+TM/O+zPzPMfyP0kh835XzvrsB+LOIeFpj+bLG8r4UH7pjGaG5IjN/lplPUHzwHzrCdiP14V4UX2L+LjP/kJmLKUaP/ucY2tBJX87MXzb6fwGb/7/XAN/NzGsycwPw/1F8WTqq6bEXZOayxmMHfLjxf/6AIsman5m/zcwHgB9RfGlV9xhTejembOn1Afh2Zv57Zm6kSHT2AP5PZv6+8T4c8TykzLw+M5dm5sbMXALMZ3hfnpeZ6zLztsb/eMhW/F8vA36VmV9rvG7zgbuBE5seOy8z72is39Ao+2Rmrs7MO4DbgR9k5r2ZuYpiNLOrscQkSc2eDlzRGGJ+lOLXmieBP46ISRHx8SiG71dT/CoDMHCYxquAE4D7oxhSf/5ITxKbh+vXNILn14Crga9HxIMR8cmImEzxxv9dDj5u//52/5mIeGZEfKcxbL0a+Pum9g5YNuT//9uB/7/RB3tR/OJ0D/AuioD124j4erQ43EDaCr7vynnf3UDxK+kLKUZxrqf48P8z4EeNLxTtaj63Zy0wbYTtRuvDRzLzsaZt76cYWamDkf6/PWh6rRt9tozB7W5+7Qb8pun+uhbLI/WfOsOY0qMxpY3XBwb/r3sB9zeSlC2KiCMi4odRHJq3iiKhHtqX7ca/0bYdFEsahsbA2scSkyQ1Wwa8NDN3abpNbfwa+HrgFcCxFMOe+zQeEwCZeUtmvoLiEItvUfw6CcXw6SCZOa3p9t+ZuSEzz8vMAyl+sXw5xS+uD1Ecx7tj08P3brr/e4ph/aIhxQmMuzWt/zzFLxf7Z+bOFMf9x9DmDPn/Pzrk/9+h8QsImfmvmTmHIgAn8ImWvSiNje+7ct53NwD/g+JLzQ0Ukxm8gOILzQ0jPGZYP43FKH34IPBHEbFT0+Z7Aw90oh0lepCin4HihGqKL2HN7a5LWzUyY0qPxhS28Pq0eI5lwN7ResKKVm35V2AhsFdmTqc4VHZoX5ZhUCxpGBoDax9LTJL62+QoTrScGhFTKQ4H+WhEPB0gInaLiFc0tt2JYrh6JUUw+/uBSqKYgvbUiJjeGDJdDQz8wvIbYNchQ8WDRMSLIuLZjcC4mmLIfmNm3g8sAs5rPMccBg/V/hKYGsWJiJOBcyiG0Qfs1KhvTUQ8C3j7Fvrji8DbGr+0RETs2Kh7p4g4ICJeHBHbURxzv67pf5TaNeg91/hguxDfd+N+32Xmrxrr3wDckJmrG/3wKkb+QvMbYGZs5cnAo/ThMuDHwMcar/Ms4C3ASNN8/wbYJ7Y8w1Or/YeI2K4RwwGmNNZtzRefBcDLIuKYxmv7txT734+3oi51hzFluJ6NKYzy+ozgZxRJ6Mcb/+fUiHjBKG3ZiWKU+w8RcThFUtYJVwHPjIjXR8S2EfEa4EDgOx16vo4wSepvV1EEgIHbUyh+YfhBRDxGceLnEY1tv0oxVPoAxclzPx1S1xuB+6IYHn4bcCpAZt5NcczrvVEMe7ca1n4axbG+qykOC7iBYtgeijfwEcAjwAcb7aBR9yrgbyiSuwcofo1qniHnPY3HP0YRNC8drTMycxHwV8BnKU4wvIfixEcoAvbHgRUUw8tPBd43Wn1SC0Pfcx8CPoPvu7LedzcAKxtJysByACOdV/FvwB3AryNixWjtHMFoffg6il+BHwSuAD6YmdeOUM83Gn9XRsRo54C02n+gOKl8HcWhLFc37g/9FXeLMvMXFF8I/5Giz08ETszM9WOtS11jTBmix2PKll6fQTLzSYr36TMoJkJYTnFu4Uht+Rvg/MZ+cS6bRwpLlZkrKUYS/5Yi4Xsv8PLM3Jo+qUxk1n60S9okiukqn5GZb6i6LVK/8H0nqUzGFPUCR5IkSZIkqUllSVIUF7/6YRQXLbsjIt7ZYpuIiAuiuBDfkoh4bhVtlVRvxhNJZTGeSIIKD7eLiN2B3TPz51HMAHQrxQWt7mza5gTgf1FMR3kE8JnMPKJlhZL6lvFEUlmMJ5KgwpGkzHwoGxcpa1xL4i6GX0PiFcBXs/BTYJdG8JKkTYwnkspiPJEENTknKSL2obiK7s1DVu3J4ItNLac+F+OTVEPGE0llMZ5I/avVxae6KiKmAZcD72rMQb+19ZwOnA6w4447HvasZz2rpBZKGq9bb711RWbutuUtx6eMeGIskerNeCKpLKPFk0qTpMZFwy4HLsnMb7bY5AGKq30PmMkIVyzPzIuAiwBmz56dixYtKrm1krZWRNzfhecoJZ4YS6R6M55IKsto8aTK2e0C+Gfgrsz81AibLQT+Z2MWmSOBVZn5UNcaKaknGE8klcV4IgmqHUl6AcWVnZdGxOJG2fuBvQEy80KKK0mfQHG15LXAad1vpqQeYDyRVBbjiaTqkqTMvAmILWyTwDu60yJJvcp4IqksxhNJUIOJG6SJYMOGDSxfvpw//OEPVTelUlOnTmXmzJlMnjy56qZIPct4UjCeSONjLNlsa+KJSZJUguXLl7PTTjuxzz77UBzO3n8yk5UrV7J8+XL23Xffqpsj9SzjifFEKoOxpLC18aQW10mSet0f/vAHdt11174OQhHBrrvu6i9W0jgZT4wnUhmMJYWtjScmSVJJ+j0IgX0glcX3kn0glcH3UWFr+sEkSZIkSZKamCRJFXnNF37Ca77wk6qbIWkCMJ5IKoOxZDOTJGkCe+Mb30hEtH0D+Na3vkVEcPfdd1fcekl1YjzpPL+gqh/0Sixxdjupy2Z/5BpWrFm/aXmfs74LwIxpU1h0znGlPc9DDz3EAQccQHE5j/bNnz+f2bNnM3/+fM4777zS2iOpfMYTSWUwlgznSJLUZc1BqJ3yrXXppZdy6qmnblo++eSTOeecc3jhC1/I3nvvzbXXXjvsMWvWrOH666/n4osvZv78+YPWXXbZZRx55JEccsghzJkzh4cffrjU9koaO+OJpDIYS4YzSZImiJUrV3L11VdvWr7//vsHXQ9g6dKl7LLLLtx444185jOf4ZJLLhlWx7e//W2OPfZYDjnkEKZNm8att966ad2LXvQifvrTn3Lbbbdx3HHHsWDBgs7+Q5IqYzyRVIZejiUmSdIEcd999/Ge97yHRx99lLvuuotnP/vZm9atXbuWVatW8e53vxsorsK9yy67DKtj/vz5zJ07F4C5c+cO+sVm3rx5HH744RxyyCH80z/9E1OnTu3sPySpMsYTSWXo5VhikiRNEIcddhinnHIK3/jGN7j88ss55ZRTNq278847Oeyww5g0aRIAS5Ys4eCDDx70+EceeYSbb76Z448/HigC0aWXXkpm8tWvfpWf/exn/Nu//Ru33XYbBxxwAAcddFD3/jlJXWU8kVSGXo4lJknSBPKGN7yBr33ta6xdu5addtppU/nSpUs59NBDNy0vWbKEWbNmDXrsZZddxgknnMB2220HwH777cfuu+/Oj370I5YuXcpRRx3FtGnTuPzyy/nxj3886NcgSROP8URSGXo1lji7ndRlM6ZNaXki5IxpU8Zd93777ceTTz7J0UcfPah86dKlHHHEEZuWb7/99mG/1syfP5/bbruNffbZZ1PZypUrmT9/PmeccQYnn3wyl1xyCS95yUvYb7/92HHHHcfdXknjYzyRVAZjyXAx1in4esHs2bNz0aJFVTdDfeSuu+7iT//0T8f0mIFrYVz6188vtS033XQTRx11FNtsU81Acau+iIhbM3N2JQ0aB2OJqmA82azf4kmnXkf1J2PJYGONJ44kSRXp1IfgnDlzOlKvpPoynkgqg7Fks0rPSYqIL0XEbyPi9hHWHx0RqyJiceN2brfbKKn+jCWSymI8kQTVjyTNAz4LfHWUbX6UmS/vTnMk9ah5GEsklWMexhOp71U6kpSZNwKPVNkGSb3PWCKpLMYTSdAbU4A/PyJui4jvRcSIk59HxOkRsSgiFj388MPdbJ+k3mAskVQW44k0wdU9Sfo58PTMPAT4R+BbI22YmRdl5uzMnL3bbrt1q32SeoOxRFJZjCdSH6h1kpSZqzNzTeP+VcDkiJhRcbMk9RhjiaSyGE+k/lDrJCkinhYR0bh/OEV7V1bbKkm9xlgiqSzGE6k/VD0F+HzgJ8ABEbE8It4SEW+LiLc1Nnk1cHtE3AZcALw2J+LVb9Wfvvyy4tZBb3zjG4mItm8A3/rWt4gI7r777kF1XX/99bzxjW/saHu3lrFEfc94UhrjifqasWSTSqcAz8zXbWH9Zymm4ZQ0Rg899BAHHHAAY/3snj9/PrNnz2b+/Pmcd955m8pvu+02nvOc55TdzFIYS6TOMp4MWm88kbZSL8WSWh9uJ01YSxbA8lvg/pvg0wcXyyW79NJLOfXUUzctn3zyyZxzzjm88IUvZO+99+baa68d9pg1a9Zw/fXXc/HFFzN//vxB6xYvXswDDzzAEUccwX777cf1119fepslbQXjiaQyGEsGMUmSum3JArjyTHjy8WJ51bJieZzBaOXKlVx99dWblu+//3723XffTctLly5ll1124cYbb+Qzn/kMl1xyybA6vv3tb3PsscdyyCGHMG3aNG699dZN62677TZ22mknbr75Zi688EI+8IEPjKu9kkpgPJFUBmPJMCZJUrdddz5sWDe4bMO6onwc7rvvPt7znvfw6KOPctddd/HsZz9707q1a9eyatUq3v3udxdPt2EDu+yyy7A65s+fz9y5cwGYO3fupl9sNmzYwIoVK3j/+98PwKGHHsqKFSvG1V5JJTCeTDjnrvw/HT8nRBrGWDKMSZLUbauWj628TYcddhinnHIK3/jGN7j88ss55ZRTNq278847Oeyww5g0aRIAS5Ys4eCDDx70+EceeYSbb76Z448/HigC0aWXXkpmcvfdd/OMZzyDKVOmAPDzn/+cQw45ZFztlVQC44mkMhhLhjFJkrpt+syxlY/BG97wBr72ta+xdu1adtppp03lS5cu5dBDD920vGTJEmbNmjXosZdddhknnHAC2223HQD77bcfu+++Oz/60Y9YvHgx//Vf/8Xjjz/OmjVrOO+883jXu9417vZKGifjiaQyGEuGqXR2O6kvHXNucZxv87D25O2L8nHab7/9ePLJJzn66KMHlS9dupQjjjhi0/Ltt98+7Nea+fPnc9ttt7HPPvtsKlu5ciXz589nxx135OSTT+aoo45i3bp1fOADH+DII48cd3sljZPxRFIZjCXDmCRJ3TarOK6Wb59RnCA5fa8iCA2Uj9MnPvEJjjrqqEFln/rUpwYt33vvvcMe98Mf/rCU55fURcYTSWUwlgxjkiRVYdZcuPUrxf3Tvltq1XPmzCm1Pkk1ZzyRVAZjySAmSVJVSg5AkvqY8URSGYwlmzhxgyRJkiQ1MUmSJEmSpCYmSVJJMrPqJlTOPpDK4XvJPpDK4PuosDX9YJIklWDq1KmsXLmyr4NRZrJy5UqmTp1adVOknmY8MZ5IZTCWFLY2njhxg1SCmTNnsnz5ch5++OGqm1KpqVOnMnPm+C88J/Uz40nBeCKNj7Fks62JJyZJUgkmT57MvvvuW3UzJE0AxhNJZTCWjI+H20mSJElSk0qTpIj4UkT8NiJuH2F9RMQFEXFPRCyJiOd2u42SeoPxRFIZjCWSoPqRpHnA8aOsfymwf+N2OvD5LrRJUm+ah/FE0vjNw1gi9b1Kk6TMvBF4ZJRNXgF8NQs/BXaJiN270zpJvcR4IqkMxhJJUP1I0pbsCSxrWl7eKBsmIk6PiEURschZPCS10FY8MZZI2gK/m0h9oO5JUtsy86LMnJ2Zs3fbbbeqmyOpRxlLJJXFeCL1rronSQ8AezUtz2yUSdJYGU8klcFYIvWBuidJC4H/2ZhJ5khgVWY+VHWjJPUk44mkMhhLpD5Q6cVkI2I+cDQwIyKWAx8EJgNk5oXAVcAJwD3AWuC0aloqqe6MJ5LKYCyRBBUnSZn5ui2sT+AdXWqOpB5mPJFUBmOJJGjzcLuI2D4iDuh0YyRJkiSpaltMkiLiRGAx8P3G8qERsbDD7ZIkSZKkSrQzkvQh4HDgUYDMXAzs27EWSZIkSVKF2kmSNmTmqiFl2YnGSJIkSVLV2pm44Y6IeD0wKSL2B84EftzZZkmSJElSNdoZSfpfwEHA48C/AquAd3WwTZIkSZJUmS2OJGXmWuDsxk2SJEmSJrR2Zre7JiJ2aVp+SkRc3dFWSZIkSVJF2jncbkZmPjqwkJm/A57asRZJkiRJUoXaSZI2RsTeAwsR8XSc3U6SJEnSBNXO7HZnAzdFxA1AAP8DOL2jrZIkSZKkirQzccP3I+K5wJGNondl5orONkuSJEmSqtHOSBLAdsAjje0PjAgy88bONUuSJEmSqrHFJCkiPgG8BrgD2NgoTsAkSZIkSdKE085I0l8AB2Tm4x1uiyRJkiRVrp3Z7e4FJnfiySPi+Ij4RUTcExFntVj/5oh4OCIWN25v7UQ7JPU+44mkshhPpB715ZcVtxK0M5K0FlgcEdcBm0aTMvPM8TxxREwCPgccBywHbomIhZl555BNL83MM8bzXJImNuOJpLIYT6QxGEhITvtute3ogHaSpIWNW9kOB+7JzHsBIuLrwCuAoUFIkrbEeCKpLMYTDTeBkwG11s4U4F/p0HPvCSxrWl4OHNFiu1dFxAuBXwLvzsxlLbaR1N+MJ5LKYjyRtOVzkiJi/4i4LCLujIh7B27daBxwJbBPZs4CrgFGTNgi4vSIWBQRix5++OEuNU9SD2krnhhLJLXBeNJNJZ5nIrWrnYkbvgx8HngCeBHwVeBfSnjuB4C9mpZnNso2ycyVTbPqXQwcNlJlmXlRZs7OzNm77bZbCc2T1ENKiyfGEqnvGU9MSqS2zknaPjOvi4jIzPuBD0XErcC543zuW4D9I2JfiuDzWuD1zRtExO6Z+VBj8STgrnE+p6SJyXgiqSxdjyezP3INK9asB+BlvA9WA2d9lxnTprDonOPGU7WkrdROkvR4RGwD/CoizqAIGNPG+8SZ+USjvquBScCXMvOOiDgfWJSZC4EzI+IkilGsR4A3j/d5JU08xhNJZakingwkSO2WS+q8dpKkdwI7AGcCHwZeDLypjCfPzKuAq4aUndt0/33A+8p4LkkTm/FE6mE1mznMeCKpndntbmncXQOc1tnmSOoZNftSI0mSVJYRk6SI+H+Z+a6IuBLIoesz86SOtkySpPEymR/OPpGkLRptJOlrjb//XzcaIm3iB7gkSdLo/L7UUSMmSZl5a0RMAk7PzFO72CZJQxkItSV120fq1h6pxmZMm9JykoYZ06ZU0Jr6ueOhVQAcVHE71F9GPScpM5+MiKdHxJTMnFhTrNThA7wObRiqjm2S6s73jaRxGJjm+zVf+Annrvw/HLT79ErjiUmJ1N7sdvcC/x4RC4HfDxRm5qc61ipJkiRJtTaRE+p2kqT/bNy2AXbqbHMkSZIkbclETlDqoJ0pwM/rRkOk2qrBoVQGQvUa99kR1CCeSBo7Y1r/2WKSFBG7Ae+l2C+mDpRn5os72K6+UMc3XB3bVDX7RL3A/bQ31OF1qkMbJKkTyoxv27SxzSXA3cC+wHnAfcAtoz1AkiRJknpVO+ck7ZqZ/xwR78zMG4AbIsIkSR3jr5ySymI8kSRtjXaSpA2Nvw9FxMuAB4E/6lyTusMPTkllMZ5IkjSxtJMkfSQipgN/C/wjsDPw7o62SlLtmRhIKovxRFLdtHNO0s2ZuSozb8/MF2XmYZm5sOMtkyRJE9+SBey//m4OXL8UPn0wLFlQdYukwdxHe0PJr1M7SdK/R8QPIuItEfGUcT2bpInBDwz1AvfT+luyAK48kylsIABWLYMrz/S1Un24j/aGDrxO7Vwn6ZkRcTjwWuDsiLgT+Hpm/stWP2ufm/2Ra1ixZj3wvqLgrOJ6GTOmTWHROcf1d5saX2oms6H4UnPMuTBrbveev0lt+qRumgIRsDkQQWWvVT+r235am/bUbD+tTb/UzXXnw4Z1g8s2rCvKjSd9rxbvG/fR3tCB16mdc5LIzJ8BP4uIvwc+BXwFGHeSFBHHA58BJgEXZ+bHh6zfDvgqcBiwEnhNZt433uet+ot48YZvv7wbatGmmn2pqUOf1OIDYqiafmBUFk8qVof9tJ3n7Xp7araf1qFfahlPVi0fW3mX9Gs8qds+Uof3TR330dq8TjX6YbsTr9MWD7eLiJ0j4k0R8T3gx8BDwOFb/Yyb650EfA54KXAg8LqIOHDIZm8BfpeZzwA+DXxivM/rsGl9PfTN97X8UvPQN99XTYNqoBYfEEPV8AOjsngCHtJVUxtH2B9HKu8HK9as56RtbuKmKWdy73av56YpZ3LSNjdVG0+mzxxbeRdUGk8qVsvPnIo9xK5jKu+GOrxOZ5/3AdZe/o5B36fXXv4Ozj7vA11rwyAdiCXtnJN0G3AocH5mPjMz/y4zb93qZ9zscOCezLw3M9cDXwdeMWSbV1CMWgFcBhwTETGeJ63LF/FWH1RVq7pNf5wrRihf2dV21E3Vr8tQdfzAoKJ4UpcPibrtI3Voz4MbW++PI5X3g5O2uYmPT76YmdusYJuAmdus4OOTL650fzl79StZm1MGla3NKZy9+pUVtQioKJ7URR3ev3Vqz8fWz225j35sfX8favf2J/+VHWJwUrZDrOftT/5rNQ065lyYvP3gssnbF+VbqZ0kab/MfHdm/mSrn6W1PYFlTcvLG2Utt8nMJ4BVML5vYnX4Il7HD6o6tOnBnDFCuV9q6rSv1PQDo5J4UocPibrtI3VpzyefaL2ffvKJ/v1i895tF7TcX9+7bXWjn5esO5KzNryV5RtnsDGD5RtncNaGt3LJuiMraxMVxZM6qMv7t07tWbhxTst9dOHGOV1rQx3tEa2/T+8R1fywPXvhUzjz96cNep3O/P1pzF649XPObTFJyszc6tq7KCJOj4hFEbHo4YcfHnG7OnwRr+MHVR3a5Jea4erwugw10T8w2o0lUI8PibrtI3Vpz0TfT7dGHfbXVhZunMOc9Rew3+OXMGf9BRPqNRpLPKmDurx/69aeibyPbq06fJ9utmLN+pav03gOQWxnJKlTHgD2alqe2ShruU1EbAtMpzhBcpjMvCgzZ2fm7N12223EJ63DF/E9tmn9gTRSeTfUoU11+1IzY9qUlsP8M6ZN2fKDS+KXmraVFk/ajSVQjw+JOrx323neKtpTp/20DvGkDvtrj6gkntRB3T5z6hRP6qbqwxDr8H2606pMkm4B9o+IfSNiCsUU40MvUrsQeFPj/quBfxvvyFYdvohvM8JJZCOVd0Nd2lSnLzWLTvodF+z45UHD/Bfs+GUWnfS7rrXhN9H6S81vwi81Q1QST+rwIVGX9+6WnrfK+FYHdYgnn5/0+pb76+cnvb5rbegRlcSTOqjbZ47xpLVTt/9py8MQT93+p11rQx2+T3faiElSRPxjRFww0m28T9w4hvcM4GrgLmBBZt4REedHxEmNzf4Z2DUi7gH+N3DWeJ8XavBFvAMnl41bDdo00i+q3fyldZDRphDukt1P/ljL12X3kz/WtTb0gqriSS0+JGrw3q1je4wnw330gx9mh1d9jvVMJgGm78UOr/ocH/3gh7vWhl5Q5feTqtXuM6cG8aR2sQT46M5XtDwM8aM7X9G1NsyYNqXl9+kq+6Vso10naVHj7wsopsC8tLF8CnBnGU+emVcBVw0pO7fp/h8az1eaGdOmtDw+sasvamMO+fXf/Bsms4GYvle1c8vXpE0Dc/vf8ffFF8yD3l/xjH91mOq6Bq/LULV4D7VQVTxZuGYOC9fPGVbeNXXbR2rSHuPJCGbN5VffKX7nPOjd1c+qajypmZq8f+vUntrFEqhFPKllv5RsxCQpM78CEBFvB+Y0flkhIi4EftSd5pWvNi9qzT6ogHq2qUrTZxbX0WpV3k01e11q8x6qgdr0Rc32kdq1pw7qEk9qpjbvIW1Wt/dv3dpTB8aTYTrxg8toI0kDngLsDDzSWJ7WKJMmtmPOLS403HyITNWHRUrqTcYTSWUxngzTiR9c2kmSPg78R0T8EAjghcCHxv3MqqWDdp9edRPqY2A4/9tnwJOPQ9WHHUjqXTU4bEjSBGE86YotJkmZ+eWI+B5wRKPo7zLz151tllQTs+bCrY2Lqp/23WrbIvUQf3BpYdZcphhPJJXBwxA7rp2RJIBJwMON7Z8ZEc/MzBs716z+4JeI1uwXaex830iSVJ4tJkkR8QngNcAdwMZGcQI9nST5hUK9pI77ax3bVBX7or58bYazTyRpy9oZSfoL4IDMfLzDbZEkSX3IxE1S3bSTJN0LTAZMkiRJktRVJtGt2S+d1U6StBZYHBHX0ZQoZeaZHWuV+psnNPcGXydJmpD88i21lyQtbNwkSeotJvO9wddJNWfi2NpE7pd2pgD/Sjca0nUGZEllMZ5IklS5MpO2dma32x/4GHAgMHWgPDP3K60VkiRJklQT7Rxu92Xgg8CngRcBpwHbdLJRqpC/iEuSJKnPtZMkbZ+Z10VEZOb9wIci4lbg3A63beIzIZE0URnfJEk9rJ0k6fGI2Ab4VUScATwATOtssyQN4hdObYn7iNrlviKNne+bvtNOkvROYAfgTODDFIfcvamTjZJqxcAoqSzGE0nqCe3MbndL4+4aivORxi0i/gi4FNgHuA+Ym5m/a7Hdk8DSxuJ/Z+ZJZTy/pInDeCKpLMYTqceV+ENUVRMwnAVcl5n7A9c1lltZl5mHNm4GIEmtGE8klcV4IgmoLkl6BTBw/aWvAH9RUTsk9T7jiaSyGE8kAdUlSX+cmQ817v8a+OMRtpsaEYsi4qcR8RejVRgRpze2XfTwww+X2VZJ9VZqPDGWSH3NeCIJaO9isp8EPgKsA74PzALenZn/soXHXQs8rcWqs5sXMjMjIkeo5umZ+UBE7Af8W0Qszcz/bLVhZl4EXAQwe/bskeqT1IO6GU+MJdLEZjzRhFGHiWDq0IYOaWd2u5dk5nsj4pUUJzGeDNwIjJokZeaxI62LiN9ExO6Z+VBE7A78doQ6Hmj8vTcirgeeA7RMkiRNXMYTSWUxnrShbl9869Ye9YV2DrcbSKReBnwjM1eV8LwL2TyN+JuAbw/dICKeEhHbNe7PAF4A3FnCc0uaWIwnkspiPJEEtJckfSci7gYOA66LiN2AP4zzeT8OHBcRvwKObSwTEbMj4uLGNn8KLIqI24AfAh/PTIOQpKGMJ5LKYjyRBLR3naSzGuclrcrMJyPi9xSzv2y1zFwJHNOifBHw1sb9HwPPHs/zSJr4jCeSymI8kTSgnXOSAPYAjo2IqU1lX+1AeyRJkiSpUu3MbvdB4GjgQOAq4KXATZgkSZIkSZqA2jkn6dUUQ8+/zszTgEOA6R1tlSRJkiRVpJ0kaV1mbgSeiIidKabD3KuzzZIkSZKkarRzTtKiiNgF+CJwK7AG+EknGyVJkiRJVWlndru/ady9MCK+D+ycmUs62yxJkiRJqsaISVJEPCsz746I57ZY99zM/HlnmyZJkiRJ3TfaSNL/Bk4H/qHFugRe3JEWSZIkSVKFRkySMvP0xt8Xda85kiRJklStdq6TNBX4G2AOxQjSj4ALM/MPHW6bJEmSJHVdO7PbfRV4DPjHxvLrga8Bp3SqUZIkSZJUlXaSpIMz88Cm5R9GxJ2dapAkSZIkVamdi8n+PCKOHFiIiCOARZ1rkiRJkiRVZ7QpwJdSnIM0GfhxRPx3Y/npwN3daZ4kSZIkdddoh9u9vGutkCRJkqSaGPFwu8y8f7TbeJ40Ik6JiDsiYmNEzB5lu+Mj4hcRcU9EnDWe55Q0MRlPJJXFeCJpQDvnJHXC7cDJwI0jbRARk4DPAS8FDgReFxEHjrS9pL5lPJFUFuOJJKC92e1Kl5l3AUTEaJsdDtyTmfc2tv068ArAmfUkbWI8kVQW44mkAVWNJLVjT2BZ0/LyRpkkjZXxRFJZOhZPLn3+Mg568pdw/03w6YNhyYIyqpW0FTo2khQR1wJPa7Hq7Mz8dgee73TgdIC999677OolVaib8cRYIk1stY0nSxbAlWfCk48Xy6uWFcsAs+aW2SxJbehYkpSZx46zigeAvZqWZzbKRnq+i4CLAGbPnp3jfG5JNdLNeGIskSa22saT686HDesGl21YV5SbJEldV+fD7W4B9o+IfSNiCvBaYGHFbZLUm4wnksrSmXiyavnYyiV1VCVJUkS8MiKWA88HvhsRVzfK94iIqwAy8wngDOBq4C5gQWbeUUV7JdWX8URSWSqNJ9Nnjq1cUkdVNbvdFcAVLcofBE5oWr4KuKqLTZPUY4wnkspSaTw55tziHKTmQ+4mb1+US+q6Oh9uJ0mS1B9mzYUTL4BJ2xXL0/cqlj0fSapEJSNJkiRJGmLWXLj1K8X9075bbVukPudIkiRJkiQ1MUmSJEmSpCYmSZIkSZLUxCRJkiRJkpqYJEmSJElSE5MkSZIkSWpikiRJkiRJTUySJEmSJKmJSZIkSZIkNTFJkiRJkqQmJkmSJEmS1MQkSZIkSZKamCRJkiRJUhOTJEmSJElqUkmSFBGnRMQdEbExImaPst19EbE0IhZHxKJutlFSbzCeSCqL8UTSgG0ret7bgZOBL7Sx7Ysyc0WH2yOpdxlPJJXFeCIJqChJysy7ACKiiqeXNIEYTySVxXgiaUDdz0lK4AcRcWtEnF51YyT1NOOJpLIYT6QJrmMjSRFxLfC0FqvOzsxvt1nNnMx8ICKeClwTEXdn5o0jPN/pwOkAe++991a1WVI9dTOeGEukic14IqkdHUuSMvPYEup4oPH3txFxBXA40DJJysyLgIsAZs+eneN9bkn10c14YiyRJjbjiaR21PZwu4jYMSJ2GrgPvITihEpJGhPjiaSyGE+k/lDVFOCvjIjlwPOB70bE1Y3yPSLiqsZmfwzcFBG3AT8DvpuZ36+ivZLqy3giqSzGE0kDqprd7grgihblDwInNO7fCxzS5aZJ6jHGE0llMZ5IGlDbw+0kSZIkqQomSZIkSZLUxCRJkiRJkpqYJEmSJElSE5MkSZIkSWpikiRJkiRJTUySJEmSJKmJSZIkSZIkNTFJkiRJkqQmJkmSJEmS1MQkSZIkSZKamCRJkiRJUhOTJEmSJElqYpIkSZIkSU1MkiRJkiSpiUmSJEmSJDWpJEmKiP8bEXdHxJKIuCIidhlhu+Mj4hcRcU9EnNXlZkrqAcYTSWUxnkgaUNVI0jXAwZk5C/gl8L6hG0TEJOBzwEuBA4HXRcSBXW2lpF5gPJFUFuOJJKCiJCkzf5CZTzQWfwrMbLHZ4cA9mXlvZq4Hvg68olttlNQbjCeSymI8kTRg26obAPwlcGmL8j2BZU3Ly4EjRqokIk4HTm8sromIXwOrmjaZ3rQ8fci6GcCKsTV7i4Y+x3i3H239SOuGlo9luRf6ZEvbtFrXTtlE7pexlI/WD2Ptk6dvYX1Zxh1PWsSSX1BuX4zVWPePdh7j+2bs25TRL53+7KnDvtKqvOx9pZfjyWjfTQYv/2X4vmm9rtf7ZWvfN0PLjLFbXh5fPMnMjtyAa4HbW9xe0bTN2cAVQLR4/KuBi5uW3wh8dgzPf9FIyy3WLerA/39RmduPtn6kdaP1QRt9VPs+2Zp+aadsIvfLWMrbfQ91ok9atKUn4kkd9o92HuP7ppp+6fRnTx32lbHuG93YV1q0r7J4UnVf+L6pvl+29n0zWr/0ep+Uta+U3S8dG0nKzGNHWx8RbwZeDhyTjf9kiAeAvZqWZzbK2nXlKMtD13XCWJ9jS9uPtn6kdaP1QTvLZSu7T7a0Tat17ZRN5H4ZS3nV76FN+jyebE39ZceTfn/fjLSubu+bOuwrrcqr3lcGqTieVN0Xvm9aq+Pn8Jb6yX2lveWtFq3f/50VEccDnwL+LDMfHmGbbSlOmjyGIvjcArw+M+/oQHsWZebssuvtZfZJa/bLcFX3SZ3iSdV9UVf2S2v2y3BV94nxpP7sl+Hsk9bG2y9VzW73WWAn4JqIWBwRFwJExB4RcRVAFidOngFcDdwFLOhEgtRwUYfq7WX2SWv2y3BV90md4knVfVFX9ktr9stwVfeJ8aT+7Jfh7JPWxtUvlYwkSZIkSVJdVTWSJEmSJEm1ZJIkSZIkSU1MkiRJkiSpiUmSJEmSJDUxSZIkSZKkJiZJkiRJktTEJEmSJEmSmpgkSZIkSVITkyRJkiRJamKSJEmSJElNTJIkSZIkqYlJkiRJkiQ1MUmSJEmSpCYmSZIkSZLUxCRJkiRJkpqYJEmSJElSE5MkSZIkSWpikiRJkiRJTUySJEmSJKmJSZIkSZIkNTFJkiRJkqQmJkmSJEmS1MQkSZIkSZKamCRJkiRJUhOTJEmSJElqYpIkSZIkSU1MkiRJkiSpiUmSJEmSJDUxSZIkSZKkJiZJkiRJktSk0iQpIqZGxM8i4raIuCMizmuxzXYRcWlE3BMRN0fEPhU0VVKNGUsklcV4IgmqH0l6HHhxZh4CHAocHxFHDtnmLcDvMvMZwKeBT3S3iZJ6gLFEUlmMJ5KqTZKysKaxOLlxyyGbvQL4SuP+ZcAxERFdaqKkHmAskVQW44kkqH4kiYiYFBGLgd8C12TmzUM22RNYBpCZTwCrgF272khJtWcskVQW44mkbatuQGY+CRwaEbsAV0TEwZl5+1jriYjTgdMBdtxxx8Oe9axnldtQSVvt1ltvXZGZu3XyOYwlUn8wnkgqy2jxpPIkaUBmPhoRPwSOB5oD0QPAXsDyiNgWmA6sbPH4i4CLAGbPnp2LFi3qfKMltSUi7u/WcxlLpInNeCKpLKPFk6pnt9ut8SsNEbE9cBxw95DNFgJvatx/NfBvmTn02GBJfcxYIqksxhNJUP1I0u7AVyJiEkXCtiAzvxMR5wOLMnMh8M/A1yLiHuAR4LXVNVdSTRlLJJXFeCKp2iQpM5cAz2lRfm7T/T8Ap3SzXZJ6i7FEUlmMJ5KgBrPbSZIkSVKdmCRJkiRJUhOTJEmSJElqYpIkSZIkSU1MkiRJkiSpiUmSJEmSJDUxSZIkSZKkJiZJkiRJktTEJEmSJEmSmpSWJEXE0yPi2Mb97SNip7LqliRJkqRuKSVJioi/Ai4DvtAomgl8q4y6JUmSJKmbyhpJegfwAmA1QGb+CnhqSXVLkiRJUteUlSQ9npnrBxYiYlsgS6pbkiRJkrqmrCTphoh4P7B9RBwHfAO4sqS6JUmSJKlrykqSzgIeBpYCfw1cBZxTUt2SJEmS1DXbllTP9sCXMvOLABExqVG2tqT6JUmSJKkryhpJuo4iKRqwPXBtSXVLkiRJUteUlSRNzcw1AwuN+zuUVLckSZIkdU1ZSdLvI+K5AwsRcRiwrqS6JUmSJKlryjon6V3ANyLiQSCApwGvKaluSZIkSeqaUpKkzLwlIp4FHNAo+kVmbiijbkmSJEnqprJGkgCeB+zTqPO5EUFmfrXE+iVJkiSp40pJkiLia8CfAIuBJxvFCZgkSZIkSeopZY0kzQYOzMwsqT5JkiRJqkRZs9vdTjFZgyRJkiT1tLJGkmYAd0bEz4DHBwoz86SS6pckSZKkrigrSfpQSfVIkiRJUqXKmgL8hjLqkSRJkqSqlXJOUkQcGRG3RMSaiFgfEU9GxOoy6pYkSZKkbipr4obPAq8DfgVsD7wV+FxJdUuSJElS15SVJJGZ9wCTMvPJzPwycHxZdUuSJElSt5Q1ccPaiJgCLI6ITwIPUWICJkmSJEndUlYi80ZgEnAG8HtgL+BVJdUtSZIkSV1T1ux29zfurgPOK6NOSZIkSarCuEaSImJB4+/SiFgy9LaFx+4VET+MiDsj4o6IeGeLbY6OiFURsbhxO3c87ZU0MRlPJJXFeCIJxj+SNBA4Xr4Vj30C+NvM/HlE7ATcGhHXZOadQ7b7UWZuTf2S+ofxRFJZjCeSxpckZeZDETEJmJeZLxrrYykmeCAzH4uIu4A9gaFBSJJGZTyRVBbjiSQoYeKGzHwS2BgR07e2jojYB3gOcHOL1c+PiNsi4nsRcdDWPoek/mA8kVQW44nUv8qaAnwNsDQirqGY3Q6AzDxzSw+MiGnA5cC7MnP1kNU/B56emWsi4gTgW8D+I9RzOnA6wN577701/4OkHldGPDGWSALjidTvypoC/JvAB4AbgVubbqOKiMkUAeiSzPzm0PWZuToz1zTuXwVMjogZrerKzIsyc3Zmzt5tt922/j+R1JPKiifGEknGE0llTQH+lbE+JiIC+Gfgrsz81AjbPA34TWZmRBxOkdStHFdjJU04xhNJZTGeSIKSkqSI2B/4GHAgMHWgPDP3G+VhL6C4CO3SiFjcKHs/sHfjsRcCrwbeHhFPUFyD6bWZmWW0WdKEYjyRVBbjiaTSzkn6MvBB4NPAi4DT2MKhfJl5ExBb2OazwGdLaqOkCcp4IqksxhNJUN45Sdtn5nVAZOb9mfkh4GUl1S1JkiRJXVPWSNLjEbEN8KuIOAN4AJhWUt2SJEmS1DVljSS9E9gBOBM4DHgD8KaS6pYkSZKkrilrJOnJxlSYayjOR5IkSZKknlTWSNI/RMRdEfHhiDi4pDolSZIkqetKSZIy80UUs9o9DHwhIpZGxDll1C1JkiRJ3VTWSBKZ+evMvAB4G7AYOLesuiVJkiQ1LFkAnz4YPrRL8XfJgqpbNOGUdTHZPwVeA7yK4orTlwJ/W0bdkiRJkhqWLIArz4QN64rlVcuKZYBZc6tr1wRT1kjSl4DfAX+emUdn5ucz87cl1S1JkiQJ4LrzNydIAzasK8pVmlJGkjLz+WXUI0mSJGkUq5aPrVxbpbRzkiRJkiR12PSZYyvXVjFJkiRJknrFMefC5O0Hl03evihXaUpNkiJihzLrkyRJktRk1lw48QKYvhcQxd8TL3DShpKVNbvdUcDFwDRg74g4BPjrzPybMuqXJEmS1DBrrklRh5U1kvRp4M8ppv8mM28DXlhS3ZIkSZLUNaWMJAFk5rKIaC56sqy6JUmSJMHsj1zDijXrh5XPmDaFReccV0GLJqaykqRljUPuMiImA+8E7iqpbkmSJEnQMkEarVxbp6zD7d4GvAPYE3gAOLSxLEmSJEk9payLya4ATi2jLkmSJEmqUlmz230ZyKHlmfmXZdQvSZIkSd1S1jlJ32m6PxV4JfBgSXVLkiRJUteUdbjd5c3LETEfuKmMuiVJkiQVZkybMuLsdipPaVOAD7E/8NQO1S1JkiT1Jaf57o6yzkl6jMHnJP0a+Lsy6pYkSZKkbhp3khTFFWQPysz/LqE9kiRJklSpcV8nKTMT+G4JbZEkSZKkypV1MdmfR8TzSqpLkiRJkipT1sQNRwCnRsT9wO+BoBhkmlVS/ZIkSZLUFWUlSX9eUj2SJEmSVKmyDrf7SGbe33wDPlJS3ZIkSVI1liyATx8MH9ql+LtkQdUtUheUNZJ0UPNCREwCDiupbkmSJKn7liyAK8+EDeuK5VXLimWAWXOra5c6blwjSRHxvsY1kmZFxOrG7THgt8C3S2mhJEmSVIXrzt+cIA3YsK4o14Q2riQpMz+WmTsB/zczd27cdsrMXTPzfSW1UZIkSeq+VcvHVq4Jo5RzkkyIJEmSNOFMnzm2ck0YZU3cMGYRsVdE/DAi7oyIOyLinS22iYi4ICLuiYglEfHcKtoqqd6MJ5LKYjzRIMecC5O3H1w2efuiXBNaWRM3bI0ngL/NzJ9HxE7ArRFxTWbe2bTNS4H9G7cjgM83/kpSM+OJpLIYT7TZwOQM151fHGI3fWaRIDlpw4Q3riQpIm4FbgK+B1yfmX9o97GZ+RDwUOP+YxFxF7An0ByEXgF8NTMT+GlE7BIRuzceK0mA8URSeYwnGmbWXJOiPjTew+2OAK4AjgZuiIirIuKdEfHMsVQSEfsAzwFuHrJqT2BZ0/LyRpkktWQ8kVQW44nUv8Y1kpSZTwDXN25ExB7A8cBHIuJPgJsz829GqyMipgGXA+/KzNVb25aIOB04HWDvvffe2mok9bAy4omxRBIYT6R+V+o5SZn5IPAl4EsRsQ3w/NG2j4jJFAHoksz8ZotNHgD2alqe2Shr9dwXARcBzJ49O8feekm9rKx4YiyRZDyResvsj1zDijXrh5XPmDaFRecct1V1dmzihszcCPz7SOsjIoB/Bu7KzE+NsNlC4IyI+DrFoX2rPN5X0lDGE0llMZ6oWSe+fKt8rV6j0crbUeXsdi8A3ggsjYjFjbL3A3sDZOaFwFXACcA9wFrgtO43U1IPMJ5IKovxRJt04su3ekNlSVJm3gTEFrZJ4B3daZGkXmU8kVQW44kkKOlishHxzIi4LiJubyzPiohzyqhbkiRJkrqplCQJ+CLwPmADQGYuAV5bUt2SJEmS1DVlJUk7ZObPhpQ9UVLdkiSN3ZIF8OmD4UO7FH+XLKi6RZKkDpgxbcqYyttR1jlJKxrXRUqAiHg1jatVS5LUdUsWwJVnwoZ1xfKqZcUywKy51bVLUk+ZMW3KiLPbqT46MdNgWUnSOyiuA/CsiHgA+C/gDSXVLUnS2Fx3/uYEacCGdUW5SZKkNjnNd/8qJUnKzHuBYyNiR2CbzHysjHolSdoqq5aPrVyV8Bo0kuqqrNnt/j4idsnM32fmYxHxlIj4SBl1S5I0ZtNnjq1clfAaNJLqqqyJG16amY8OLGTm7ygusiZJUvcdcy5M3n5w2eTti/KqOJGEJPWMss5JmhQR22Xm4wARsT2wXUl1S5I0NgPnHV13fnGI3fSZRYJU1flITiQhST2lrCTpEuC6iPhyY/k04Csl1S1J0tjNmlufBMSJJCSpp5Q1ccMnImIJcEyj6MOZeXUZdUuS1POcSEKSekpZI0lk5veA75VVnyRJE8b0mcUhdq3K+5jXoJFUmiULSj3EupQkKSJOBj4BPBWIxi0zc+cy6pckqZedvfqVnJ0XskNsTgjW5hQ+uvqVfLTCdlXNab4llaID532WNbvdJ4GTMnN6Zu6cmTuZIEmSVLhk3ZGcteGtLN84g40ZLN84g7M2vJVL1h1ZddMkqfeNdt7nVirrcLvfZOZdJdUlSdKEs3DjHBaun1N1MyRp4unAeZ9lJUmLIuJS4FvA4wOFmfnNkuqXJKltsz9yzYjnuniIV82UfB6BpD7UgfM+y0qSdgbWAi9pKkvAJEmS1HWtEqTRylURrx8l9aTa/RB1zLmDYwmM+wLiZU0BfloZ9UiSpD7i9aOknlS7H6I6cAHxsma3mwq8BTgImDpQnpl/WUb9kiT1Mqe6HoHXj5JUlpIvIF7W4XZfA+4G/hw4HzgVcCIHSZJwqusRef0oSTVV1hTgz8jMDwC/z8yvAC8DjiipbkmSNBEdc25x3kCzcZ5HIEllKGskaUPj76MRcTDwa4oLy0qS1HUe3tYjOnAegSaO2k0OoL5SVpJ0UUQ8BfgAsBCYBvgzkCSpEn6B6iEln0egiWPFmvWctM1NvHfbBewRK3gwZ/DJJ+aycI3XG6taP/wQVdbsdhc37t4A7FdGnZIkSepfJ21zEx+ffDE7RPFlfGas4OOTL24cv/SyStvW7/rhh6hxJUkR8YbM/JeI+N+t1mfmp8ZTvyRJkvrTe7ddsClBGrBDrOe92y4APlZNo9Q3xjuStGPj707jbYgk9aQlCzyfQpI6YI9YMUL5yi63RP1oXElSZn4hIiYBqzPz0yW1SZJ6w5IFg6/wvWpZsQz9nSiZOEoqwYM5g5ktEqUHc1ecJF6dNu4pwDPzSeB1JbRFknrLdedvTpAGbFhXlPergcRx1TIgNyeOSxZU3TJJPebzk17P2hw8EcDanMLnJ72+ohapn5R1naR/j4jPRsT/iIjnDtxKqluS6mnV8rGV9wMTR0kl+egHP8wOr/ocTN8LCJi+Fzu86nN89IMfrrpp6gNlTQF+aONv86dgAi8uqX5Jqp/pMxsjJi3K+5WJo6QyOUV8fU3wQ6vLmgL8RWXUI0k95ZhzB5+TBDB5+6K8X5k4StLE1wfn5JZ1uB0R8bKIeG9EnDtwK6tuSaqlWXPhxAsGHQrCiRdMmA+IrXH26le2PIfg7NWvrKhFkqTS9cGh1aWMJEXEhcAOwIuAi4FXAz8ro25JqjUPBRnkknVH8tg2T/DebRewR6zkwdyVTz4xl4Ubj+SjVTdOklSOPji0uqxzko7KzFkRsSQzz4uIfwC+V1LdkqQesnDjHBaun1N1MyRJndIHh1aXdbjdwHjb2ojYA9gA7L6lB0XElyLitxFx+wjrj46IVRGxuHHzED5JwxhLJJXFeCK14Zhzi3Nwm02wc3LLGkn6TkTsAvxf4OcUM9t9sY3HzQM+C3x1lG1+lJkvH28DJU1o8zCWSCrHPIwn0ugGDjN3drvRZebAhPWXR8R3gKmZuaqNx90YEfuU0QZJ/auqWDL7I9ewYs36YeUzpk1h0TnHdbs5kkrgdxOpTRP8nNxSDreLiCUR8f6I+JPMfLydBGkMnh8Rt0XE9yLioBLrldRfSo8lrRKk0cr7wYxpU8ZULvUov5tIE1xZh9udCLwGWBARG4FLgQWZ+d/jrPfnwNMzc01EnAB8C9i/1YYRcTpwOsDee+89zqeVNMEYS7rEEbSam+AXf+wS44nUB0oZScrM+zPzk5l5GPB6YBbwXyXUuzoz1zTuXwVMjogZI2x7UWbOzszZu+2223ifWtIEYiyR2Hzxx1XLgNx88cclC6puWU8xnkj9ocyLyT49It4LfB14FvDeEup8WkRE4/7hFO1dOd56JfWXCR9LliyATx8MH9ql+OuXXrXSBxd/7AbjidQfyrqY7M3AZGABcEpm3tvm4+YDRwMzImI58MFGPWTmhRQXpX17RDxBMc34azMzy2izpImjr2PJwOjAwJffgdEB8DAqDdYHF38sg/HEeCJBeeck/c/M/MVYH5SZr9vC+s9STMMpSSOqKpbMmDZlxNntuma00QG/1KhZH1z8sQx9/d3EeCJtUtYU4GNOkCSp19VikgJHB9SuY84dPEoAE+7ijxon44m0SWnnJEmSKjDSKICjAxpq1lw48QKYvhcQxd8TL3CEQJs8xK5jKpcmMpMkSeplx5xbjAY0c3RAI5k1F959O3zo0eKvCZKafGz9XNbm4MOF1+YUPrbe/UT9Z1yH20XEyaOtz8xvjqd+SdIWDHzJ9do3ksZp4cY5sAHeu+0C9oiVPJi78skn5rJw4xwuqLpxUpeN95ykExt/nwocBfxbY/lFwI8BkyRJ6rRZc02KJJVi4cY5LFw/p+pmSJUbV5KUmacBRMQPgAMz86HG8u7AvHG3TpI0qtkfuWbEGfZqMbGEasN9RZLaV9Y5SXsNJEgNvwH2LqluSdIIWn3pHa1c/ct9RZLaV9Z1kq6LiKuB+Y3l1wDXllS3JEmSOqwW136TaqKs6ySdERGvBF7YKLooM68oo25JkiR1noddSpuVNZIE8HPgscy8NiJ2iIidMvOxEuuXJEmSpI4r5ZykiPgr4DLgC42iPYFvlVG3JEmSJHVTWRM3vAN4AbAaIDN/RTEtuCSpg0Y6V8BzCDSU+4okta+sw+0ez8z1EQFARGwLZEl1S5JG4DkEapf7iiS1r6yRpBsi4v3A9hFxHPAN4MqS6pYkSZKkrikrSToLeBhYCvw1cFVmnl1S3ZIkSZLUNWUdbve/MvMzwBcHCiLinY0ySZIkSeoZZY0kvalF2ZtLqluSJEmSumZcI0kR8Trg9cC+EbGwadVOwCPjqVuSJEmSqjDew+1+DDwEzAD+oan8MWDJOOuWJEmSpK4bV5KUmfcD9wPPL6c5kiRJklStUs5JiogjI+KWiFgTEesj4smIWF1G3ZIkSZLUTWVN3PBZ4HXAr4DtgbcCnyupbkmSJEnqmrKmACcz74mISZn5JPDliPgP4H1l1a8+sWQBXHc+rFoO02fCMefCrLlVt0pNZn/kGlasWT+sfMa0KSw657gKWiRJklSuspKktRExBVgcEZ+kmMyhrFEq9YslC+DKM2HDumJ51bJiGUyUaqRVgjRauSRJUq8pK5F5IzAJOAP4PbAX8KqS6la/uO78zQnSgA3rinJJkiSpS0oZSWrMcgewDjivjDrVh1YtH1u5JEmS1AFlzW738oj4j4h4JCJWR8Rjzm6nMZs+c2zlqsxJ29zETVPO5N7tXs9NU87kpG1uqrpJkiRJpSnrcLv/B7wJ2DUzd87MnTJz55LqVr845lyYvP3gssnbF+WqjZO2uYmPT76YmdusYJuAmdus4OOTLzZRkiRJE0ZZSdIy4PbMzJLqUz+aNRdOvACm7wVE8ffEC5y0oWbeN2UBO8TgSRp2iPW8b8qCilokSZJUrrJmt3svcFVE3AA8PlCYmZ8qqX71i1lz65UUOSX5MLuzckzlkiRJvaasJOmjwBpgKjClpDqlajkleWvTZxZ90apckiRpAigrSdojMw8uqS71qdpdpHS0Kcn7OUk65tzBySN47pgkSZpQyjon6aqIeElJdakqSxbApw+GD+1S/F3S3XNManeRUqckb81zxyRJ0gRX1kjS24H3RMTjwAYggHSGu95x9nkf4OyNF24+IX/VMtZe/g4+esVSPvrBD1fbuIo8xK7szooRyvtc3c4dkyRJKlEpI0mNKb+3ycztxzIFeER8KSJ+GxG3j7A+IuKCiLgnIpZExHPLaG/VIyZ19PYn/7XljGVvf/JfK2pR9T62fi5rc/ApdmtzCh9bb3JQR5XFE0kTirFEEowzSYqIZzX+PrfVrY0q5gHHj7L+pcD+jdvpwOfH015g88n4q5YBuflk/D5PlPaI4SMmRXn/zli2cOMcztrwVpZvnMHGDJZvnMFZG97Kwo1zqm6aWptHt+OJpIloHsYSqe+N93C7/00RIP6hxboEXjzagzPzxojYZ5RNXgF8tXH9pZ9GxC4RsXtmPrS1DfZk/NYezBnMbJEoPZi70s9zli3cOIeF602KekEl8UTShGMskQTjTJIy8/TG3Zdm5h+a10XE1PHU3bAnxYVqByxvlI0aiFauXMm8efMGlR100EE873nPY8OqX3MJpwx7zKGr7uRQYO3atSxYMHxUafbs2Rx88MGsWrWKK664Ytj65z//+RxwwAGsWLGC73znO8PWv/CFL2S//fbj17/+Nd///veHrT/mmGPYa6+9WLZsGdddd92w9ccffzxPe9rTuPfee7nxxhuHrX/5y1/OjBkz+MUvfsFPfvKTYetf+cpXMn36dG6//XYWLVo0bP1/PDGXj02+mF/GM1jMQQBsZBvujd15xrx5nHrqqUyePJlbbrmFO+64Y9jj3/zmNwPw4x//mF/+8peD1k2ePJlTTz0VgBtuuIH/+q//GrR+hx12YO7cucyYNoWn/+FedttmzaD1T267eVf6/ve/z69//etB63fddVdOPPFEAK688kpWrhw8+vW0pz2N448vfhT85je/yerVqwetnzlzJsceeywACxYsYO3atQAcP6Wo56GNO3PbE3sAcNyUXzKJjcyb9/Cmxz/zmc/kqKOOAhi230HTvrdhA5dccsmw9YceeiiHHnrohN73amLM8WTUWNLHr+dosWTu3LnssMMOLF68mMWLFw9b341YAnDttdeyfPngSVZ23nlnTj75ZKC7sWTAvvvuy5/92Z8BcMkll7Bhw4ZB640l7e17NVD+d5M+f02NJ8aToeoQT8qauOHHwNDD61qVdUxEnE4xqsWee+458oY77wmrW5Rv/5TONKxHfH/jEeQG+IvJt0LA40xmWT6VFbkzz+hSGxadc9yIgUjqlrZjiSRtgfFE6l1RjBZv5YMjnkbx68m/AK+nmNUOYGfgwsx8Vht17AN8p9V1liLiC8D1mTm/sfwL4OgtDWnPnj07W/0iAQy/QCgU13jp8ymMa3eNohqwT8oTEbdm5uwuPM8+lBhPRo0lkirRjXjS9e8mkioxWjwZ70jSnwNvBmZSnJc0kCQ9Brx/nHUDLATOiIivA0cAq8Z9zO9AInTd+cX1bqbPLC6C2ccJEuCX/hbskwmn/HgiqR8ZS6Q+MN5zkr4CfCUiXpWZl4/18RExHzgamBERy4EPApMbdV8IXAWcANwDrAVOG097N6nDNV6WLDBRk0pUWTyRNKEYSyRBeeckzYyInSlGkL5IcS7SWZn5g9EelJmv28L6BN5RUhvrY+ghfwPTkIOJkrSV+jaeSCqVsUQSlJck/WVmfiYi/hzYFXgj8DVg1CSpCrU4z8RpyDUWjjpKkiR1VVlJ0sC5SCdQXDvgjoiI0R5QlVYJ0mjlHbFq+djK1b8cdZQkSeq6bUqq59aI+AFFknR1ROwEbCyp7oln+giXZx2pXP1rtFFHSZIkdURZSdJbgLOA52XmWmAKnsg4orNXv5K1OWVQ2dqcwtmra3PRTdWFo46SJEldV1aSlMCBQOM4IHYEppZU94RzybojOWvDW1m+cQYbM1i+cQZnbXgrl6w7suqmqW4cdZQkSeq6ss5J+ieKw+teDJxPMcvd5cDzSqp/wlm4cQ4L18+puhmqu2PObX3x42POra5NkiRJE1xZI0lHZOY7gD8AZObvKA65q50Z01o3a6RyqVKz5sKJF8D0vYAo/p54gZM2SJIkdVBZI0kbImISxWF3RMRu1HTihq5N8y2VpQ4XP5YkSeojZY0kXQBcATw1Ij4K3AT8fUl1S5IkSVLXlDKSlJmXRMStwDEU10z6i8y8q4y6J6IZ06aMeEFbSZIkSdUq63A7MvNu4O6y6pvIPORPkiRJqq/SkiRJ5Zv9kWtGHHU02ZYkSeqMss5JktQBrRKk0colSZI0fiZJkiRJktTEJEmSJEmSmpgkSZIkSVITkyRJkiRJamKSJNXYSNfO8ppakiRJneMU4FKNOc23JElS9zmSJEmSJElNTJIkSZIkqUlkZtVtKF1EPAw8CqxqKp7etDx9yLoZwIqSmzH0Oca7/WjrR1o3tHwsy73QJ1vaptW6dsomcr+MpXy0fhhrnzw9M3fbwja104gl91NuX4zVWPePdh7j+2bs25TRL53+7KnDvtKqvOx9pZfjyaP4vqnb+2akNo13+/F+Dm+pn9xXtrw8vniSmRPyBlw00nKLdYs6/fzj3X609SOtG60P2uij2vfJ1vRLO2UTuV/GUt7ue6gTfVK3W5V9Mdb9o53H+L6ppl86/dlTh31lrPtGN/aVOt2q7gvfN9X3y9a+b0brl17vk7L2lbL7ZSIfbnflKMtD13Xj+ce7/WjrR1o3Wh+0s1y2svtkS9u0WtdO2UTul7GUV/0eqpMq+2Jr6i87nvT7+2akdXV739RhX2lVXvW+UidV94Xvm9bq+Dm8pX5yX2lveatNyMPtxioiFmXm7KrbUSf2SWv2y3D2yWb2RWv2S2v2y3D2yWb2RWv2y3D2SWvj7ZeJPJI0FhdV3YAask9as1+Gs082sy9as19as1+Gs082sy9as1+Gs09aG1e/OJIkSZIkSU0cSZIkSZKkJiZJkiRJktTEJEmSJEmSmpgkDRER+0XEP0fEZVW3pU4i4i8i4osRcWlEvKTq9tRBRPxpRFwYEZdFxNurbk+dRMSOEbEoIl5edVuqZDxpzXgynPFkZMaTgvFkOGNJa8aT1sYaS/oiSYqIL0XEbyPi9iHlx0fELyLinog4CyAz783Mt1TT0u4aY798KzP/Cngb8Joq2tsNY+yTuzLzbcBc4AVVtLdbxtIvDX8HLOhuK7vDeNKa8WQ440lrxpPNjCfDGUtaM54M1+lY0hdJEjAPOL65ICImAZ8DXgocCLwuIg7sftMqNY+x98s5jfUT1TzG0CcRcRLwXeCq7jaz6+bRZr9ExHHAncBvu93ILpmH8aSVeRhPhpqH8aSVeRhPBszDeDLUPIwlrczDeDLUPDoYS/oiScrMG4FHhhQfDtzT+GVmPfB14BVdb1yFxtIvUfgE8L3M/Hm329otY91XMnNhZr4UOLW7Le2uMfbL0cCRwOuBv4qICRVnjCetGU+GM560ZjzZzHgynLGkNePJcJ2OJduW29yesiewrGl5OXBEROwKfBR4TkS8LzM/VknrqtOyX4D/BRwLTI+IZ2TmhVU0riIj7StHAycD2zGxf6kZSct+ycwzACLizcCKzNxYQdu6zXjSmvFkOONJa8aTzYwnwxlLWjOeDFdaLOnnJKmlzFxJcWyrmmTmBcAFVbejTjLzeuD6iptRW5k5r+o2VM140prxZDjjyeiMJ8aTVowlrRlPRjaWWDKhhq3H6AFgr6blmY2yfme/DGeftGa/bGZftGa/DGeftGa/bGZfDGeftGa/DFdan/RzknQLsH9E7BsRU4DXAgsrblMd2C/D2Set2S+b2Ret2S/D2Set2S+b2RfD2Set2S/DldYnfZEkRcR84CfAARGxPCLekplPAGcAVwN3AQsy844q29lt9stw9klr9stm9kVr9stw9klr9stm9sVw9klr9stwne6TyMzyWitJkiRJPa4vRpIkSZIkqV0mSZIkSZLUxCRJkiRJkpqYJEmSJElSE5MkSZIkSWpikiRJkiRJTUySJEmSJKmJSZIkSZIkNTFJkiRJkqQmJkmSJEmS1MQkSZIkSZKamCRJkiRJUhOTJEmSJElqYpIkSZIkSU1MkiRJkiSpiUmSJEmSJDUxSZIkSZKkJiZJkiRJktTEJEmSJEmSmpgkSZIkSVITkyRJkiRJamKSJEmSJElNTJIkSZIkqYlJkiRJkiQ1MUmSJEmSpCYmSZIkSZLUxCRJkiRJkpqYJEmSJElSE5MkSZIkSWpikiRJkiRJTUySJEmSJKmJSZIkSZIkNTFJkiRJkqQmJkmSJEmS1MQkSZJUGxHx1xFxVUR8LiJWRMSDEXFc1e2SJPUXkyRJUp0cAhwJLASeCnwB+LtKWyRJ6jsmSZKkOpkFfDwzr87MjcCdI20YEftExEu61zRJUr8wSZIk1UJEBPBs4Mqm4oMZOVHaBzBJkiSVziRJklQX+wDbAr9oKnsOsDgijo6IH0TElRFxS0Q8G3g78JqIuD4i/igijoyImyPihxHxoe43X5I0UWxbdQMkSWqYBSxtHGY34DnAecBOwA7AnwPPAj4BfApYlpnvAYiIlwHnZeZVEeGPgJKkreaHiCSpLmYBiwcWImIG8DTg9kbRf2ThLmD3Fo//HHBCRFwCHN/htkqSJjBHkiRJtZCZHx6yvAKYDFCcrsShjfOWngk8BGwAJjU9ZFVmnhERU4Bbgau60W5J0sRjkiRJ6hWrKCZ1+GPgLcB9wMci4hvAXwFvjoiTKT7b5lXURknSBBCZWXUbJEkaVUQcDbx84PwjSZI6yXOSJEmSJKmJI0mSJEmS1MSRJEmSJElqYpIkSZIkSU1MkiRJkiSpiUmSJEmSJDUxSZIkSZKkJiZJkiRJktTEJEmSJEmSmpgkSZIkSVKT/x8J2O1wtVbMzwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1008x576 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2, 3, figsize=(14, 8), sharex=True, sharey=False)\n",
    "\n",
    "for k, (title, func) in enumerate((\n",
    "    (\"Least-squares\", linear),\n",
    "    (\"Least-squares with soft L1 norm\", soft_l1),\n",
    "    (\"Least-squares with arctan norm\", arctan),\n",
    ")):\n",
    "    ax[0, k].set_title(title)\n",
    "    for i, x in enumerate(\"ab\"):\n",
    "        ax[0, k].errorbar(\n",
    "            n_pts * 0.95 + 0.1 * i,\n",
    "            np.sqrt(n_pts) * func[0 + i],\n",
    "            np.sqrt(n_pts * func[4 + i]),\n",
    "            fmt=\"so\"[i],\n",
    "            label=f\"$\\sqrt{{n}}\\,\\Delta {x}$\",\n",
    "        )\n",
    "        ax[1, k].plot(\n",
    "            n_pts * 0.95 + 0.1 * i,\n",
    "            func[2 + i] / func[4 + i],\n",
    "            \"so\"[i],\n",
    "            label=f\"$\\sqrt{{n}}\\,\\Delta {x}$\",\n",
    "        )\n",
    "        ax[0, k].legend()\n",
    "plt.semilogx()\n",
    "for i in range(3):\n",
    "    ax[1, i].axhline(1, ls=\"--\", color=\"0.5\")\n",
    "    ax[0, i].set_ylim(-2, 2)\n",
    "    ax[1, i].set_ylim(0.7, 3)\n",
    "ax[0, 0].set_ylabel(\"bias and variance\")\n",
    "ax[1, 0].set_ylabel(\"estimated variance / true variance\")\n",
    "fig.supxlabel(r\"$n_\\mathrm{pts}$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normal least-squares fit has a smallest variance, which is equal to the minimum variance for this problem given by the Cramer-Rao bound. The robust fits use less information to achieve robustness, hence the variance is larger. The loss from the soft L1 norm in this case is nearly negligible, but for the arctan norm it is noticable.\n",
    "\n",
    "**Beware**: The variance estimate obtained from the fit is wrong for robust least-squares, since the robust least-squares is not even asymptotically a maximum-likelihood estimator. The estimate is significantly larger than the actual variance for the soft_l1 and arctan norms in this case."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}