diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index a87ddb48..f0dac19d 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -44,7 +44,7 @@ repos: hooks: - id: pydocstyle files: src/iminuit/[^_].*\.py - additional_dependencies: [toml] # update when https://github.com/PyCQA/pydocstyle/pull/608 goes in! + additional_dependencies: [tomli] # Python linter (Flake8) - repo: https://github.com/pycqa/flake8 diff --git a/doc/notebooks/template_model_mix.ipynb b/doc/notebooks/template_model_mix.ipynb new file mode 100644 index 00000000..628dce49 --- /dev/null +++ b/doc/notebooks/template_model_mix.ipynb @@ -0,0 +1,788 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fitting a mixture of templates and parametric models\n", + "\n", + "The class `iminuit.cost.Template` supports fitting a mixture of templates and parametric models. This is useful if some components have a simple shape like a Gaussian peak, while other components are complicated and need to be estimated from simulation or a control measurement.\n", + "\n", + "In this notebook, we demonstrate this usage. Our data consists of a Gaussian peak and exponential background. We fit the Gaussian peak with a parametric model and use a template to describe the exponential background." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "from iminuit import Minuit\n", + "from iminuit.cost import Template, ExtendedBinnedNLL\n", + "from numba_stats import norm, truncexpon\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We generate the toy data." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAp7klEQVR4nO3deXRUZZ7G8aeyB0gqbNm6gwSGVUGQJUagsSFtRpaGkW6lzSAIDT2QMEKOsjSb7MthE0QYFgE9YLrpaRhGaGwI2iiEgDEoCkSUINiaoGOTsHQSktz5w2O1RdKShKrUW8n3c06dQ269de/vvgTr8X3fe6/NsixLAAAABvHxdAEAAAC3I6AAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIzj5+kCaqK8vFxffPGFQkJCZLPZPF0OAACoAsuydO3aNUVHR8vH54fHSLwyoHzxxReKiYnxdBkAAKAGLl++rB//+Mc/2MYrA0pISIikb08wNDTUw9UAAICqKCwsVExMjON7/Id4ZUD5blonNDSUgAIAgJepyvIMFskCAADjEFAAAIBxCCgAAMA4XrkGBQBgnrKyMt26dcvTZcCDfH195efn55JbgBBQAAB37fr16/r8889lWZanS4GHNWjQQFFRUQoICLir/RBQAAB3paysTJ9//rkaNGig5s2bcwPNesqyLJWUlOirr75Sbm6u2rRpc8ebsf0QAgoA4K7cunVLlmWpefPmCg4O9nQ58KDg4GD5+/vrs88+U0lJiYKCgmq8r2pHmyNHjmjw4MGKjo6WzWbTnj17nN63LEuzZ89WVFSUgoODlZCQoPPnzzu1+eabb5SUlKTQ0FCFhYVpzJgxun79eo1PAgDgeYycQNJdjZo47ae6H7hx44buv/9+rVu3rtL3ly1bpjVr1mjDhg3KzMxUw4YNlZiYqKKiIkebpKQkffTRRzp48KBef/11HTlyROPGjav5WQAAgDql2gHl0Ucf1YIFC/Rv//ZvFd6zLEurV6/WzJkzNWTIEHXu3FmvvPKKvvjiC8dIy9mzZ3XgwAFt3rxZcXFx6t27t9auXau0tDR98cUXd31CAABUxcMPP6xJkya5bf+jRo3S0KFD3bZ/T7h48aJsNptOnTrl9mO5dA1Kbm6u8vLylJCQ4Nhmt9sVFxenjIwMDR8+XBkZGQoLC1P37t0dbRISEuTj46PMzMxKg09xcbGKi4sdPxcWFrqybACAOzxvr+XjFdTu8eBWLr1RW15eniQpIiLCaXtERITjvby8PIWHhzu97+fnpyZNmjja3G7x4sWy2+2OF08yBgDURyUlJZ4uodZ4xZ1kp0+froKCAsfr8uXLni4JAFAHlJaWKiUlRXa7Xc2aNdOsWbMc93J59dVX1b17d4WEhCgyMlJPPvmkrly54vT5jz76SIMGDVJoaKhCQkLUp08fffrpp5Ue6+TJk2revLmWLl3q2LZgwQKFh4crJCREv/71rzVt2jR16dLF8f5300QLFy5UdHS02rVrJ0k6ffq0+vXrp+DgYDVt2lTjxo1zutiksumroUOHatSoUY6fW7ZsqUWLFmn06NEKCQlRixYttHHjRqfPnDhxQl27dlVQUJC6d++u7OzsKvft3XJpQImMjJQk5efnO23Pz893vBcZGVnhL7i0tFTffPONo83tAgMDHU8u5gnGAABX2b59u/z8/HTixAm98MILWrlypTZv3izp28un58+fr/fff1979uzRxYsXnb7g//rXv+onP/mJAgMDdfjwYWVlZWn06NEqLS2tcJzDhw/rZz/7mRYuXKipU6dKknbs2KGFCxdq6dKlysrKUosWLbR+/foKn01PT1dOTo7jwpIbN24oMTFRjRs31smTJ7Vr1y4dOnRIKSkp1T7/FStWOILHhAkTNH78eOXk5Ej69uZ7gwYNUseOHZWVlaXnn39ezz77bLWPUVMuXYMSGxuryMhIpaenOxJgYWGhMjMzNX78eElSfHy8rl69qqysLHXr1k3St39x5eXliouLc2U5ALxUy2n7Kmy7uGSgBypBXRcTE6NVq1bJZrOpXbt2On36tFatWqWxY8dq9OjRjnatWrXSmjVr1KNHD12/fl2NGjXSunXrZLfblZaWJn9/f0lS27ZtKxxj9+7deuqpp7R582Y98cQTju1r167VmDFj9PTTT0uSZs+erT//+c8VbrvRsGFDbd682XFn1k2bNqmoqEivvPKKGjZsKEl68cUXNXjwYC1durTCMosfMmDAAE2YMEGSNHXqVK1atUpvvvmm2rVrp507d6q8vFxbtmxRUFCQ7r33Xn3++eeO73N3q/YIyvXr13Xq1CnHCt7c3FydOnVKly5dks1m06RJk7RgwQLt3btXp0+f1lNPPaXo6GjHSuYOHTroX//1XzV27FidOHFCR48eVUpKioYPH67o6GhXnhsAAD/owQcfdLp/S3x8vM6fP6+ysjJlZWVp8ODBatGihUJCQtS3b19J0qVLlyRJp06dUp8+fRzhpDKZmZn65S9/qVdffdUpnEhSTk6Oevbs6bTt9p8lqVOnTk63jT979qzuv/9+RziRpF69eqm8vNwx+lFVnTt3dvzZZrM5zXKcPXtWnTt3drrZWnx8fLX2fzeqHVDeffddde3aVV27dpUkpaamqmvXrpo9e7YkacqUKZo4caLGjRvnSJoHDhxwOsEdO3aoffv26t+/vwYMGKDevXtXmPcCAMBTioqKlJiYqNDQUO3YsUMnT57U7t27Jf1joWpV7prbunVrtW/fXi+//HKNH6T4/SBSVT4+PhWei1TZ8W8PVzabTeXl5dU+njtUO6A8/PDDsiyrwmvbtm2Svj25efPmKS8vT0VFRTp06FCFIa8mTZpo586dunbtmgoKCvTyyy+rUaNGLjkhAACqKjMz0+nn48ePq02bNjp37pz+7//+T0uWLFGfPn3Uvn37CusnO3furLfffvsHg0ezZs10+PBhffLJJ3r88ced2rZr104nT550an/7z5Xp0KGD3n//fd24ccOx7ejRo/Lx8XEsom3evLm+/PJLx/tlZWX68MMP77jv24/zwQcfON1o9fjx49Xax93wiqt4AABwh0uXLik1NVU5OTl67bXXtHbtWj3zzDNq0aKFAgICtHbtWl24cEF79+7V/PnznT6bkpKiwsJCDR8+XO+++67Onz+vV199tcI0S3h4uA4fPqxz587pV7/6lWMR7cSJE7VlyxZt375d58+f14IFC/TBBx/c8ZEBSUlJCgoK0siRI/Xhhx/qzTff1MSJEzVixAjH+pN+/fpp37592rdvn86dO6fx48fr6tWr1eqbJ598UjabTWPHjtWZM2e0f/9+LV++vFr7uBsEFABAvfXUU0/p73//u3r27Knk5GQ988wzGjdunJo3b65t27Zp165d6tixo5YsWVLhy7lp06Y6fPiwrl+/rr59+6pbt27atGlTpWtSIiMjdfjwYZ0+fVpJSUkqKytTUlKSpk+frmeffVYPPPCAcnNzNWrUqDs+YK9BgwZ644039M0336hHjx76xS9+of79++vFF190tBk9erRGjhypp556Sn379lWrVq3005/+tFp906hRI/3v//6vTp8+ra5du2rGjBlOl0i7m826fZLKCxQWFsput6ugoIBLjoE6iKt4vEtRUZFyc3MVGxt7V0+vhfSzn/1MkZGRevXVVz1dSo390O9Ddb6/XXqZMQAAqJqbN29qw4YNSkxMlK+vr1577TUdOnRIBw8e9HRpRiCgAADgATabTfv379fChQtVVFSkdu3a6b//+7+dnmdXnxFQAADwgODgYB06dMjTZRiLRbIAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAcJtRo0Zp6NChtX7c559/Xl26dKn145qI+6AAANyiskcWuJMrH4fwwgsvyAufBFOnEFAAALiN3W73dAn1HlM8AIB66w9/+IM6deqk4OBgNW3aVAkJCbpx40aFKZ5r164pKSlJDRs2VFRUlFatWqWHH35YkyZNcrRp2bKlFi1apNGjRyskJEQtWrTQxo0bnY43depUtW3bVg0aNFCrVq00a9Ys3bp1q5bO1rsQUAAA9dKXX36pX/3qVxo9erTOnj2rt956S4899lilUzupqak6evSo9u7dq4MHD+rtt9/We++9V6HdihUr1L17d2VnZ2vChAkaP368cnJyHO+HhIRo27ZtOnPmjF544QVt2rRJq1atcut5eiumeAAA9dKXX36p0tJSPfbYY7rnnnskSZ06darQ7tq1a9q+fbt27typ/v37S5K2bt2q6OjoCm0HDBigCRMmSPp2tGTVqlV688031a5dO0nSzJkzHW1btmypZ599VmlpaZoyZYrLz8/bEVAAAPXS/fffr/79+6tTp05KTEzUI488ol/84hdq3LixU7sLFy7o1q1b6tmzp2Ob3W53hI7v69y5s+PPNptNkZGRunLlimPb7373O61Zs0affvqprl+/rtLSUoWGhrrh7LwfUzwAgHrJ19dXBw8e1J/+9Cd17NhRa9euVbt27ZSbm1vjffr7+zv9bLPZVF5eLknKyMhQUlKSBgwYoNdff13Z2dmaMWOGSkpK7uo86ioCCgCg3rLZbOrVq5fmzp2r7OxsBQQEaPfu3U5tWrVqJX9/f508edKxraCgQB9//HG1jnXs2DHdc889mjFjhrp37642bdros88+c8l51EVM8QAA6qXMzEylp6frkUceUXh4uDIzM/XVV1+pQ4cO+uCDDxztQkJCNHLkSD333HNq0qSJwsPDNWfOHPn4+Mhms1X5eG3atNGlS5eUlpamHj16aN++fRXCEP6BERQAQL0UGhqqI0eOaMCAAWrbtq1mzpypFStW6NFHH63QduXKlYqPj9egQYOUkJCgXr16qUOHDgoKCqry8X7+859r8uTJSklJUZcuXXTs2DHNmjXLladUp9gsL7xVXmFhoex2uwoKClhcBNRBld2B1JV3CYVrFRUVKTc3V7GxsdX6wvZmN27c0I9+9COtWLFCY8aM8XQ5Rvmh34fqfH8zxQMAwB1kZ2fr3Llz6tmzpwoKCjRv3jxJ0pAhQzxcWd1FQAEAoAqWL1+unJwcBQQEqFu3bnr77bfVrFkzT5dVZxFQAAC4g65duyorK8vTZdQrLJIFAADGIaAAAADjEFAAAC7hhReFwg1c9XtAQAEA3BVfX19J4pbtkCTdvHlTUsXb/lcXi2QBAHfFz89PDRo00FdffSV/f3/5+PD/vvWRZVm6efOmrly5orCwMEdwrSkCCgDgrthsNkVFRSk3N5dny0BhYWGKjIy86/0QUAAAdy0gIEBt2rRhmqee8/f3v+uRk+8QUAAALuHj41NvbnUP9yOgAKhVPGcHQFWwkgkAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADCOywNKWVmZZs2apdjYWAUHB6t169aaP3++LMtytLEsS7Nnz1ZUVJSCg4OVkJCg8+fPu7oUAADgpVweUJYuXar169frxRdf1NmzZ7V06VItW7ZMa9eudbRZtmyZ1qxZow0bNigzM1MNGzZUYmKiioqKXF0OAADwQn6u3uGxY8c0ZMgQDRw4UJLUsmVLvfbaazpx4oSkb0dPVq9erZkzZ2rIkCGSpFdeeUURERHas2ePhg8f7uqSAACAl3H5CMpDDz2k9PR0ffzxx5Kk999/X++8844effRRSVJubq7y8vKUkJDg+IzdbldcXJwyMjIq3WdxcbEKCwudXgAAoO5y+QjKtGnTVFhYqPbt28vX11dlZWVauHChkpKSJEl5eXmSpIiICKfPRUREON673eLFizV37lxXlwoAAAzl8hGU3//+99qxY4d27typ9957T9u3b9fy5cu1ffv2Gu9z+vTpKigocLwuX77swooBAIBpXD6C8txzz2natGmOtSSdOnXSZ599psWLF2vkyJGKjIyUJOXn5ysqKsrxufz8fHXp0qXSfQYGBiowMNDVpQIAAEO5PKDcvHlTPj7OAzO+vr4qLy+XJMXGxioyMlLp6emOQFJYWKjMzEyNHz/e1eUA8AItp+2rdpuLSwa6qxwABnB5QBk8eLAWLlyoFi1a6N5771V2drZWrlyp0aNHS5JsNpsmTZqkBQsWqE2bNoqNjdWsWbMUHR2toUOHurocAADghVweUNauXatZs2ZpwoQJunLliqKjo/Wb3/xGs2fPdrSZMmWKbty4oXHjxunq1avq3bu3Dhw4oKCgIFeXAwAAvJDN+v4tXr1EYWGh7Ha7CgoKFBoa6ulyAFRDVaZzqoIpHsD7VOf7m2fxAAAA4xBQAACAcQgoAADAOAQUAABgHJdfxQMA3+eqRbEA6hdGUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjMNlxgDqrMouceYZPoB3YAQFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHF4mjEAr8STioG6jREUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjOPn6QIAeIeW0/ZV2HZxycA7tqlNnj4+ANdhBAUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcdwSUP7617/q3//939W0aVMFBwerU6dOevfddx3vW5al2bNnKyoqSsHBwUpISND58+fdUQoAAPBCLg8of/vb39SrVy/5+/vrT3/6k86cOaMVK1aocePGjjbLli3TmjVrtGHDBmVmZqphw4ZKTExUUVGRq8sBAABeyM/VO1y6dKliYmK0detWx7bY2FjHny3L0urVqzVz5kwNGTJEkvTKK68oIiJCe/bs0fDhw11dEgAA8DIuH0HZu3evunfvrl/+8pcKDw9X165dtWnTJsf7ubm5ysvLU0JCgmOb3W5XXFycMjIyKt1ncXGxCgsLnV4AAKDucnlAuXDhgtavX682bdrojTfe0Pjx4/Wf//mf2r59uyQpLy9PkhQREeH0uYiICMd7t1u8eLHsdrvjFRMT4+qyAQCAQVweUMrLy/XAAw9o0aJF6tq1q8aNG6exY8dqw4YNNd7n9OnTVVBQ4HhdvnzZhRUDAADTuDygREVFqWPHjk7bOnTooEuXLkmSIiMjJUn5+flObfLz8x3v3S4wMFChoaFOLwAAUHe5PKD06tVLOTk5Tts+/vhj3XPPPZK+XTAbGRmp9PR0x/uFhYXKzMxUfHy8q8sBAABeyOVX8UyePFkPPfSQFi1apMcff1wnTpzQxo0btXHjRkmSzWbTpEmTtGDBArVp00axsbGaNWuWoqOjNXToUFeXAwAAvJDLA0qPHj20e/duTZ8+XfPmzVNsbKxWr16tpKQkR5spU6boxo0bGjdunK5evarevXvrwIEDCgoKcnU5AADAC7k8oEjSoEGDNGjQoH/6vs1m07x58zRv3jx3HB4AAHg5nsUDAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBw/TxcAwHu1nLbP0yUAqKMYQQEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwjp+nCwBgppbT9nm6BAD1GCMoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGcXtAWbJkiWw2myZNmuTYVlRUpOTkZDVt2lSNGjXSsGHDlJ+f7+5SAACAl3DrfVBOnjyp//qv/1Lnzp2dtk+ePFn79u3Trl27ZLfblZKSoscee0xHjx51ZzkAUOH+LheXDPRQJQB+iNtGUK5fv66kpCRt2rRJjRs3dmwvKCjQli1btHLlSvXr10/dunXT1q1bdezYMR0/ftxd5QAAAC/itoCSnJysgQMHKiEhwWl7VlaWbt265bS9ffv2atGihTIyMirdV3FxsQoLC51eAACg7nLLFE9aWpree+89nTx5ssJ7eXl5CggIUFhYmNP2iIgI5eXlVbq/xYsXa+7cue4oFQAAGMjlIyiXL1/WM888ox07digoKMgl+5w+fboKCgocr8uXL7tkvwAAwEwuDyhZWVm6cuWKHnjgAfn5+cnPz09/+ctftGbNGvn5+SkiIkIlJSW6evWq0+fy8/MVGRlZ6T4DAwMVGhrq9AIAAHWXy6d4+vfvr9OnTztte/rpp9W+fXtNnTpVMTEx8vf3V3p6uoYNGyZJysnJ0aVLlxQfH+/qcgAAgBdyeUAJCQnRfffd57StYcOGatq0qWP7mDFjlJqaqiZNmig0NFQTJ05UfHy8HnzwQVeXAwAAvJBb74Pyz6xatUo+Pj4aNmyYiouLlZiYqJdeeskTpQAAAAPZLMuyPF1EdRUWFsput6ugoID1KICb3H5Ds7qKG7UBtac63988iwcAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYByP3KgNgFnqyz1PAHgPRlAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADCOn6cLAABPajltX4VtF5cM9EAlAL6PERQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxuFpxkA9VNkTfAHAJIygAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh8uMAeAOKrss++KSgR6oBKg/GEEBAADGIaAAAADjuDygLF68WD169FBISIjCw8M1dOhQ5eTkOLUpKipScnKymjZtqkaNGmnYsGHKz893dSkAAMBLuXwNyl/+8hclJyerR48eKi0t1W9/+1s98sgjOnPmjBo2bChJmjx5svbt26ddu3bJbrcrJSVFjz32mI4ePerqcoB6j9vaA/BGLg8oBw4ccPp527ZtCg8PV1ZWln7yk5+ooKBAW7Zs0c6dO9WvXz9J0tatW9WhQwcdP35cDz74oKtLAgAAXsbta1AKCgokSU2aNJEkZWVl6datW0pISHC0ad++vVq0aKGMjIxK91FcXKzCwkKnFwAAqLvceplxeXm5Jk2apF69eum+++6TJOXl5SkgIEBhYWFObSMiIpSXl1fpfhYvXqy5c+e6s1QAcGBaDPA8t46gJCcn68MPP1RaWtpd7Wf69OkqKChwvC5fvuyiCgEAgIncNoKSkpKi119/XUeOHNGPf/xjx/bIyEiVlJTo6tWrTqMo+fn5ioyMrHRfgYGBCgwMdFepAADAMC4fQbEsSykpKdq9e7cOHz6s2NhYp/e7desmf39/paenO7bl5OTo0qVLio+Pd3U5AADAC7l8BCU5OVk7d+7U//zP/ygkJMSxrsRutys4OFh2u11jxoxRamqqmjRpotDQUE2cOFHx8fFcwQMAACS5IaCsX79ekvTwww87bd+6datGjRolSVq1apV8fHw0bNgwFRcXKzExUS+99JKrSwEAAF7K5QHFsqw7tgkKCtK6deu0bt06Vx8eAADUATzNGKhDuDy29tze1zzdGHAtHhYIAACMQ0ABAADGIaAAAADjsAbFVZ63V7KtoPbrAACgDmAEBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOFxmDBiKW6l7l8oeM8DfGVBzjKAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHy4wrc/uTiV31VGJXPfGYJycDAOo4RlAAAIBxCCgAAMA4BBQAAGAc1qB4mrvWu6DO4VbqAOoTRlAAAIBxCCgAAMA4TPG4U2WXA3sbLo32eheDnqywrWXRTg9UAgBVxwgKAAAwDgEFAAAYh4ACAACMwxoUb+CqtSzeeEmzN9bsIrevHfH0uhHWssChHv+7RO1hBAUAABiHgAIAAIzDFE9dVRcucfZSt9/xtbKpkZoMiVc6xTKNKRYAdRMjKAAAwDgEFAAAYBwCCgAAMA5rUOqKmqw5ceft52tzDUxVjuWu86pkv5WuOfFCdeU8PKnCeqRKnj5dlTa1yvRLiHlsRr3BCAoAADAOAQUAABiHKR54F3dOHdXi0La7hvVrOi1Tm9M5nr4jrSfv0Hv737tUST1VuHS8xr8vtTl9Y9pUkWn14I4YQQEAAMYhoAAAAOMQUAAAgHFYg1IV9em28a66XNmdn3PFsTw8/1xxHQbz4d9X0zUxt68n8cZLpatWsxfeDqAqx68r60Jq89YHdRgjKAAAwDgEFAAAYBwCCgAAMA5rUEzj6Tlh/INpt+v3At5wH5aqqMm9Ujx+Dqav56gj9zC647Fr+1im/T27ECMoAADAOAQUAABgHKZ4UH11ZDrCJTzcF56eUvH4tEYt8cbzrPy2+jXYkaenMEzjDTVWhelTgmIEBQAAGIiAAgAAjENAAQAAxmENCuqnujKPDCfeuFbEXWq1L0z79+Sqx2+4c12GN/SZh9elMIICAACMQ0ABAADGYYoHQL3mldNCpk0PeKO60odVmZrx0nNlBAUAABiHgAIAAIzj0YCybt06tWzZUkFBQYqLi9OJEyc8WQ4AADCEx9ag/O53v1Nqaqo2bNiguLg4rV69WomJicrJyVF4eLinygIAwLt56ZqT23lsBGXlypUaO3asnn76aXXs2FEbNmxQgwYN9PLLL3uqJAAAYAiPjKCUlJQoKytL06dPd2zz8fFRQkKCMjIyKrQvLi5WcXGx4+eCgm9XKBcWFrqnwGLLPfsFAHiPyr5j6tP3gxu+Y7/73rasO/ejRwLK119/rbKyMkVERDhtj4iI0Llz5yq0X7x4sebOnVthe0xMjNtqBADUc0vqxlRJjbnx/K9duya7/Yf37xX3QZk+fbpSU1MdP5eXl+ubb75R06ZNZbPZXHqswsJCxcTE6PLlywoNDXXpvvEP9HPtoJ9rB/1cO+jn2uOuvrYsS9euXVN0dPQd23okoDRr1ky+vr7Kz8932p6fn6/IyMgK7QMDAxUYGOi0LSwszJ0lKjQ0lH8AtYB+rh30c+2gn2sH/Vx73NHXdxo5+Y5HFskGBASoW7duSk9Pd2wrLy9Xenq64uPjPVESAAAwiMemeFJTUzVy5Eh1795dPXv21OrVq3Xjxg09/fTTnioJAAAYwmMB5YknntBXX32l2bNnKy8vT126dNGBAwcqLJytbYGBgZozZ06FKSW4Fv1cO+jn2kE/1w76ufaY0Nc2qyrX+gAAANQinsUDAACMQ0ABAADGIaAAAADjEFAAAIBx6mVAWbdunVq2bKmgoCDFxcXpxIkTP9h+165dat++vYKCgtSpUyft37+/lir1btXp502bNqlPnz5q3LixGjdurISEhDv+veBb1f19/k5aWppsNpuGDh3q3gLriOr289WrV5WcnKyoqCgFBgaqbdu2/LejCqrbz6tXr1a7du0UHBysmJgYTZ48WUVFRbVUrXc6cuSIBg8erOjoaNlsNu3Zs+eOn3nrrbf0wAMPKDAwUP/yL/+ibdu2ub1OWfVMWlqaFRAQYL388svWRx99ZI0dO9YKCwuz8vPzK21/9OhRy9fX11q2bJl15swZa+bMmZa/v791+vTpWq7cu1S3n5988klr3bp1VnZ2tnX27Flr1KhRlt1utz7//PNarty7VLefv5Obm2v96Ec/svr06WMNGTKkdor1YtXt5+LiYqt79+7WgAEDrHfeecfKzc213nrrLevUqVO1XLl3qW4/79ixwwoMDLR27Nhh5ebmWm+88YYVFRVlTZ48uZYr9y779++3ZsyYYf3xj3+0JFm7d+/+wfYXLlywGjRoYKWmplpnzpyx1q5da/n6+loHDhxwa531LqD07NnTSk5OdvxcVlZmRUdHW4sXL660/eOPP24NHDjQaVtcXJz1m9/8xq11ervq9vPtSktLrZCQEGv79u3uKrFOqEk/l5aWWg899JC1efNma+TIkQSUKqhuP69fv95q1aqVVVJSUlsl1gnV7efk5GSrX79+TttSU1OtXr16ubXOuqQqAWXKlCnWvffe67TtiSeesBITE91YmWXVqymekpISZWVlKSEhwbHNx8dHCQkJysjIqPQzGRkZTu0lKTEx8Z+2R836+XY3b97UrVu31KRJE3eV6fVq2s/z5s1TeHi4xowZUxtler2a9PPevXsVHx+v5ORkRURE6L777tOiRYtUVlZWW2V7nZr080MPPaSsrCzHNNCFCxe0f/9+DRgwoFZqri889T3oFU8zdpWvv/5aZWVlFe5WGxERoXPnzlX6mby8vErb5+Xlua1Ob1eTfr7d1KlTFR0dXeEfBf6hJv38zjvvaMuWLTp16lQtVFg31KSfL1y4oMOHDyspKUn79+/XJ598ogkTJujWrVuaM2dObZTtdWrSz08++aS+/vpr9e7dW5ZlqbS0VP/xH/+h3/72t7VRcr3xz74HCwsL9fe//13BwcFuOW69GkGBd1iyZInS0tK0e/duBQUFebqcOuPatWsaMWKENm3apGbNmnm6nDqtvLxc4eHh2rhxo7p166YnnnhCM2bM0IYNGzxdWp3y1ltvadGiRXrppZf03nvv6Y9//KP27dun+fPne7o0uEC9GkFp1qyZfH19lZ+f77Q9Pz9fkZGRlX4mMjKyWu1Rs37+zvLly7VkyRIdOnRInTt3dmeZXq+6/fzpp5/q4sWLGjx4sGNbeXm5JMnPz085OTlq3bq1e4v2QjX5fY6KipK/v798fX0d2zp06KC8vDyVlJQoICDArTV7o5r086xZszRixAj9+te/liR16tRJN27c0Lhx4zRjxgz5+PD/4K7wz74HQ0ND3TZ6ItWzEZSAgAB169ZN6enpjm3l5eVKT09XfHx8pZ+Jj493ai9JBw8e/KftUbN+lqRly5Zp/vz5OnDggLp3714bpXq16vZz+/btdfr0aZ06dcrx+vnPf66f/vSnOnXqlGJiYmqzfK9Rk9/nXr166ZNPPnEEQEn6+OOPFRUVRTj5J2rSzzdv3qwQQr4LhRaPmXMZj30PunUJroHS0tKswMBAa9u2bdaZM2escePGWWFhYVZeXp5lWZY1YsQIa9q0aY72R48etfz8/Kzly5dbZ8+etebMmcNlxlVQ3X5esmSJFRAQYP3hD3+wvvzyS8fr2rVrnjoFr1Ddfr4dV/FUTXX7+dKlS1ZISIiVkpJi5eTkWK+//roVHh5uLViwwFOn4BWq289z5syxQkJCrNdee826cOGC9ec//9lq3bq19fjjj3vqFLzCtWvXrOzsbCs7O9uSZK1cudLKzs62PvvsM8uyLGvatGnWiBEjHO2/u8z4ueees86ePWutW7eOy4zdZe3atVaLFi2sgIAAq2fPntbx48cd7/Xt29caOXKkU/vf//73Vtu2ba2AgADr3nvvtfbt21fLFXun6vTzPffcY0mq8JozZ07tF+5lqvv7/H0ElKqrbj8fO3bMiouLswIDA61WrVpZCxcutEpLS2u5au9TnX6+deuW9fzzz1utW7e2goKCrJiYGGvChAnW3/72t9ov3Iu8+eablf739ru+HTlypNW3b98Kn+nSpYsVEBBgtWrVytq6davb67RZFuNgAADALPVqDQoAAPAOBBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGOf/AbM+raqHP0BfAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "rng = np.random.default_rng(1)\n", + "\n", + "s = rng.normal(0.5, 0.05, size=1000)\n", + "b = rng.exponential(1, size=1000)\n", + "b = b[b < 1]\n", + "\n", + "ns, xe = np.histogram(s, bins=100, range=(0, 1))\n", + "nb, _ = np.histogram(b, bins=xe)\n", + "n = ns + nb\n", + "\n", + "plt.stairs(nb, xe, color=\"C1\", fill=True, label=\"background\")\n", + "plt.stairs(n, xe, baseline=nb, color=\"C0\", fill=True, label=\"signal\")\n", + "plt.legend();" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Fitting a parametric component and a template\n", + "\n", + "We now model the peaking component parametrically with a Gaussian. A template fit is an extended binned fit, so we need to provide a scaled cumulative density function like for `iminuit.cost.ExtendedBinnedNLL`. To obtain a background template, we generate more samples from the exponential distribution and make a histogram." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Migrad
FCN = 87.01 (chi2/ndof = 0.9) Nfcn = 194
EDM = 4.93e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0_n 990 40 0
1 x0_mu 0.4951 0.0020 0 1
2 x0_sigma 0.0484 0.0018 0
3 x1 630 40 0
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
x0_n x0_mu x0_sigma x1
x0_n 1.43e+03 0.000324 (0.004) 0.0185 (0.274) -441 (-0.284)
x0_mu 0.000324 (0.004) 3.87e-06 -5.26e-08 (-0.015) -0.000347 (-0.004)
x0_sigma 0.0185 (0.274) -5.26e-08 (-0.015) 3.21e-06 -0.0185 (-0.251)
x1 -441 (-0.284) -0.000347 (-0.004) -0.0185 (-0.251) 1.69e+03
" + ], + "text/plain": [ + "┌─────────────────────────────────────────────────────────────────────────┐\n", + "│ Migrad │\n", + "├──────────────────────────────────┬──────────────────────────────────────┤\n", + "│ FCN = 87.01 (chi2/ndof = 0.9) │ Nfcn = 194 │\n", + "│ EDM = 4.93e-05 (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 │ x0_n │ 990 │ 40 │ │ │ 0 │ │ │\n", + "│ 1 │ x0_mu │ 0.4951 │ 0.0020 │ │ │ 0 │ 1 │ │\n", + "│ 2 │ x0_sigma │ 0.0484 │ 0.0018 │ │ │ 0 │ │ │\n", + "│ 3 │ x1 │ 630 │ 40 │ │ │ 0 │ │ │\n", + "└───┴──────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", + "┌──────────┬─────────────────────────────────────────┐\n", + "│ │ x0_n x0_mu x0_sigma x1 │\n", + "├──────────┼─────────────────────────────────────────┤\n", + "│ x0_n │ 1.43e+03 0.000324 0.0185 -441 │\n", + "│ x0_mu │ 0.000324 3.87e-06 -5.26e-08 -0.000347 │\n", + "│ x0_sigma │ 0.0185 -5.26e-08 3.21e-06 -0.0185 │\n", + "│ x1 │ -441 -0.000347 -0.0185 1.69e+03 │\n", + "└──────────┴─────────────────────────────────────────┘" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# signal model: scaled cdf of a normal distribution\n", + "def signal(xe, n, mu, sigma):\n", + " return n * norm.cdf(xe, mu, sigma)\n", + "\n", + "# background template: histogram of MC simulation\n", + "rng = np.random.default_rng(2)\n", + "b2 = rng.exponential(1, size=1000)\n", + "b2 = b2[b2 < 1]\n", + "template = np.histogram(b2, bins=xe)[0]\n", + "\n", + "# fit\n", + "c = Template(n, xe, (signal, template))\n", + "m = Minuit(c, x0_n=500, x0_mu=0.5, x0_sigma=0.05, x1=100)\n", + "m.limits[\"x0_n\", \"x1\", \"x0_sigma\"] = (0, None)\n", + "m.limits[\"x0_mu\"] = (0, 1)\n", + "m.migrad()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABQI0lEQVR4nO3de3xT9f0/8NdpoLfQGwV6SxFUQLk4b4O1UJzan3ytKHxLFScw5pigolDwBjqVOcdNxCIT8MIUHeq0i6isYxMUqAPBIW4KDmH0q2164dqWpi2U9Pz+qOc0l5PkJD1JTprX8/HI49FzcnLyzmmSvvv5fN6fjyCKoggiIiIiHYkKdQBEREREzpigEBERke4wQSEiIiLdYYJCREREusMEhYiIiHSHCQoRERHpDhMUIiIi0h0mKERERKQ7PUIdgD/a29tRXV2NhIQECIIQ6nCIiIhIBVEUcebMGWRmZiIqynMbSVgmKNXV1cjOzg51GEREROSHyspKmEwmj8eEZYKSkJAAoOMFJiYmhjgaIiIiUqOxsRHZ2dny33FPwjJBkbp1EhMTmaAQERGFGTXDMzhIloiIiHSHCQoRERHpDhMUIiIi0h0mKERERKQ7TFCIiIhId5igEBERke4wQSEiIiLdYYJCREREusMEhYiIiHSHCQoRERHpDhMUIiIi0h0mKERERKQ7TFCIiIhId5igEBERke4wQSEi3bFarRAEAYIgwGq1hjocIgoBJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdnxOUnTt34uabb0ZmZiYEQcCmTZsc7hdFEU888QQyMjIQFxeH/Px8HD582OGYU6dOYcqUKUhMTERycjJmzJiBpqamLr0QIiIi6j58TlCsVit+9KMf4YUXXlC8f/ny5Xj++eexbt067NmzB0ajEePGjUNra6t8zJQpU3DgwAF89NFH2Lx5M3bu3ImZM2f6/yqIiIioWxFEURT9frAg4L333sPEiRMBdLSeZGZm4oEHHsCDDz4IAGhoaEBaWhpee+013H777fjmm28wdOhQfP7557j66qsBAFu2bEFBQQGqqqqQmZnp9XkbGxuRlJSEhoYGJCYm+hs+EemU9BkHgLKyMtxwww0wGAwhjoqIusqXv9+ajkGpqKhAbW0t8vPz5X1JSUkYNWoUdu/eDQDYvXs3kpOT5eQEAPLz8xEVFYU9e/Yonvfs2bNobGx0uBFR92Q2mzF06FB5u6CgAAMGDIDZbA5hVEQUbJomKLW1tQCAtLQ0h/1paWnyfbW1tejXr5/D/T169EDv3r3lY5wtWbIESUlJ8i07O1vLsIlIJ8xmM4qKimCxWBz2WywWFBUVMUkhiiBhUcWzcOFCNDQ0yLfKyspQh0REGrPZbJg7dy6Uep2lfcXFxbDZbMEOjYhCQNMEJT09HQBQV1fnsL+urk6+Lz09HceOHXO4//z58zh16pR8jLOYmBgkJiY63IioeykvL0dVVZXb+0VRRGVlJcrLy4MYFRGFiqYJysCBA5Geno5t27bJ+xobG7Fnzx7k5OQAAHJyclBfX499+/bJx3z88cdob2/HqFGjtAyHiMJITU2NpscRUXjr4esDmpqacOTIEXm7oqICX375JXr37o3+/fujuLgYTz/9NAYNGoSBAwfi8ccfR2Zmplzpc+mll+J//ud/cNddd2HdunVoa2vDfffdh9tvv11VBQ8RdU8ZGRmaHkdE4c3nMuPt27fj2muvddk/ffp0vPbaaxBFEU8++SReeukl1NfXY8yYMVizZg0GDx4sH3vq1Cncd999+PDDDxEVFYVJkybh+eefR69evVTFwDJjou7HZrNhwIABsFgsiuNQBEGAyWRCRUUFS46JwpQvf7+7NA9KqDBBIeqepCoeAA5JiiAIAIDS0lIUFhaGJDYi6rqQzYNCRNQVhYWFKC0tdenuNZlMTE6IIgxbUIhIdziTLFH3xBYUItItq9UKQRAgCAKsVqviMfbJyNixY5mcEEUgJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdJihERESkOz1CHQARkbPq+hb55wPVDYiPP48UYzSykuNCGBURBRMTFCLSFUt9C65/doe8XbR2N6KiYxHX04CtD1zDJIUoQrCLh4h05bT1nOL+ljab2/uIqPthgkJERES6wwSFiIiIdIcJChEREekOB8kSUchZ6lvk8SVHjjWFOBoi0gMmKEQUVDabTf55586dGDYyD2Of2eHhEUQUidjFQ0RBYzabMeSSS+XtgoICXDl0MJoP7QphVESkR2xBIaKgMJvNKCoqgiiKDvtPHqsBNi1G34mPIn5ILgAgKjoWFzyyORRhEpFOsAWFiALOZrNh7ty5LsmJvVPbXoLYbnN7PxFFFiYoRBRw5eXlqKqq8niM7cwJnK06EKSIiEjvmKAQUcDV1NSoOq7urUfRfq5V8b72c60YYUqGIAiwWq1ahkdEOsQEhYgCLiMjI9QhEFGYYYJCRAGXl5cHk8kEQRBCHQoRhQkmKEQUcAaDAatWrQp1GEQURpigEFFQFBYWYuWLG2Doleqw35DQB6njHwxRVESkV5wHhYiCJr/gFmQciEHVqskAgL5FixA38AqI59twMsSxEZG+sAWFiIJKiDLIP8dmD3fYJiKSMEEhooCyWq0QBAGCIKC5WZvy4APVDbDUt2hyLiLSJyYoRKQ7cT09t6oUrd2N/Gd3MEkh6sY4BoWIdKelzYaSyZfj4n69AABHjjVhzhufuRxz2noOWclxoQiRiAKMCQoRhZxzi0lcTwN+PLA3kw+iCMYEhYhCrqWtc5HA0ntykNU3hckJUYRjgkJEujIsMwlGI5MTokjHQbJEFDQ19coLAXobFEtEkYctKEQUNKebzzlsl96Tg/h4I2LQhsHLtX8+q9WKXr06Bto2NTXBaDRq/yREFBBMUIgoZDq6c4ywWrWZH4WIug928RAREZHuMEEhIiIi3WGCQkRERLrDBIWIAqrabjr6o8ebQhgJEYUTJihEFDCW+hZc/+wOefu1Xd+FMBoiCidMUIgoYE5bz3k/iIhIgeYJis1mw+OPP46BAwciLi4OF110EX77299CFEX5GFEU8cQTTyAjIwNxcXHIz8/H4cOHtQ6FiIiIwpTmCcqyZcuwdu1a/P73v8c333yDZcuWYfny5Vi9erV8zPLly/H8889j3bp12LNnD4xGI8aNG4fWVuVZJomIiCiyaJ6g7Nq1CxMmTMBNN92EAQMGoKioCDfccAP27t0LoKP1pKSkBL/+9a8xYcIEXHbZZXj99ddRXV2NTZs2aR0OEYUBm61zscCdO3c6bEvE9s59rZVfO2wTUfejeYKSm5uLbdu24dtvvwUA/Otf/8Knn36KG2+8EQBQUVGB2tpa5Ofny49JSkrCqFGjsHv3bsVznj17Fo2NjQ43IuoezGYzhg4dKm8XFBRgwIABMJvN8r6tZR+gZv298vbx0kWwrJuBrWUfBDVWIgoezae6X7BgARobG3HJJZfAYDDAZrPhd7/7HaZMmQIAqK2tBQCkpaU5PC4tLU2+z9mSJUvwm9/8RutQiSjE3n//fUydOtVhjBoAWCwWFBUVobS0FAAwf9Z0l2NsZ05g3qzp6J9qRGFhYdBiJqLg0LwF5Z133sHGjRvx5ptv4osvvsCGDRuwYsUKbNiwwe9zLly4EA0NDfKtsrJSw4iJKFQefvhhl8QDgLxv7ty5mDt3ruIxHQcCxcXFil1CRBTeNG9Beeihh7BgwQLcfvvtAIARI0bgu+++w5IlSzB9+nSkp6cDAOrq6pCRkSE/rq6uDpdffrniOWNiYhATE6N1qEQUYhaLxe19oiiiqqrKyxlEVFZW4u9//7vcjUxE3YPmLSjNzc2IinI8rcFgQHt7OwBg4MCBSE9Px7Zt2+T7GxsbsWfPHuTk5GgdDhFFAHfdw0QUvjRvQbn55pvxu9/9Dv3798ewYcOwf/9+rFy5Er/85S8BAIIgoLi4GE8//TQGDRqEgQMH4vHHH0dmZiYmTpyodThEpDNR0bH4qqoeJw7vx7XXXqvJOaWWWSLqPjRPUFavXo3HH38c9957L44dO4bMzEzMmjULTzzxhHzMww8/DKvVipkzZ6K+vh5jxozBli1bEBsbq3U4RBRizuXBcQOvAADk5eXBZDLBYrEojjERBAFZWVkA4PYYyejRozWOmohCTRA9fep1qrGxEUlJSWhoaEBiYmKowyEiN0pefgMPzp8HW9NJeZ8hoQ9WPLsSxXdNg9lsRlFREQA4JCCCIACAXMWjdIy9Pd9aEB9vRIoxGlnJcfJ+q9WKXr16AQCamppgNBo1fHVE5Ctf/n5zLR4iCgiz2Yz5s6Y7JCdAR3nw/FnTYTabUVhYiNLSUmRmZjocYzKZUFpaisLCQrfHGHqlyj8Xrd2N8as/Rf6zO2CxWz2ZiMIXExQi0pzNZvNcHozO8uDCwkIcPHhQ3l9WVoaKigqHuU2cj1nzxrvImLHG5ZwtbTYuUEjUTTBBISLNlZeXeywRFsWO8uDy8nIAHZV+krFjxzpsS+z3XTUqF0KU6zFE1H1oPkiWiCKXpb4Fp63nsPfAf1UdX1NTE+CIiChcMUEhIk1Y6lsweunHAIDW74+peoz9ZI2B4LwI4Q033KDYOkNE+sMuHiLShP3YjxjTMBgS+rg9VhAEZGdnIy8vL2DxqFmEkIj0iwkKEWlOiDKg9/Uzle/7oYS4pKQkYK0ZUvmy81T60iKETFKI9I8JChEFRPyQXPSd+KhDOTDgWEIcCJ4qiKR9XGCQSP+YoBBRwMQPyXUoB1YqIdbaF3t2+VRBRET6xASFiALKvhzYXQmxlo4fq1N1HCuIiPSNCQoRdSt9+6WpOi7QFURE1DVMUIioW7lyVC5MJpM8GNdZMCqIiKjrmKAQUbdiMBiwatUqAHBJUoJRQURE2mCCQkTdjppFCIlI3ziTLBF1S4WFhcjPz0dSUhKAjgoiziRLFD7YgkJEqlitVgiCAEEQYLVaFY9pP9eK75aNx3fLxqP9XGuQI+x4/hGmZDlGNYsQEpE+MUEhIiIi3WGCQkRhJyU+GrF2jSGtlV9DbOfMsETdCcegEFHY+XzH33Fm40Py9vHSRUhMTUPPnGkhjIqItMQEhYjCztSpU13W2mk8eQzYvCJEERGR1tjFQ0RhR2khQEBpHxGFKyYoRKQJS32L4v64nqycISLfsYuHiDRx2nrOYbv0nhzExxsRgzYMXh6ioIgobDFBIaKAGJaZBKPR6HbOFCIiT9jFQ0Rhx91CgETUfTBBIaKQMxqNEEURoijCaDR6PObPf/4zAKUkhUkLUXfCBIWIwoq7hQBT09KROv7BEEVFRFpjgkJEYaewsBAHDx6Ut8vKyvD8e/9A/KCfhDAqItISExQiCkvOCwFGcSFAom6FCQoRERHpDhMUIiIi0h3Og0JEASVV3xAR+YItKERERKQ7bEEhIr9Z6lvkKe6PHGsKcTRE1J0wQSEiv1jqWzB66cehDoOIuil28RCRX5wXByQi0hITFCJSxWazyT/v3LnTYVsitns/JpDsn//19/6KypPsdiIKV0xQiMgrs9mMoUOHytsFBQUY95MRaD60S97XfGgXatbf63DMgAEDYDabgxLjgV0fOTz/vdNuxYjhw/HS6296XOOHiPRJEMOw/q+xsRFJSUloaGhAYmJiqMMh6tbMZjOKiopcSoUFQYAoiug78VEAwPFNi10eKy3oV1paisLCQk3jslqt6NWrFwBg48aNmDp1qmI5syAIAXl+IvKdL3+/maAQkVs2mw0DBgxAVVWV22OieqVCAGBrOql4vyAIMJlMqKiocJievqvsE5SsrCxYLJagPj8R+c6Xv9/s4iEit8rLyz0mJwDQ3nTSbXICAKIoorKyEuXl5VqHJ3OXnATr+YlIe0xQiMitmpoaXZ4rHJ+fiHzDBIWI3MrIyNDlucLx+YnIN0xQiMitvLw8mEwmebCrkqheqTD0SnV7vyAIyM7ORl5eXiBCBNAxBsVdjMF4fiLSHhMUInLLYDBg1apVAOCSAEjbqfmz0Dt/luLjpWNKSkoCOkB1+fLlHu8P9PMTkfaYoBCRR4WFhSgtLUVmZqbD/rSMTPSd+Cjih+Qifkgu+k581KUlxWQyBaXEd8KECVj54gaX5zck9MHKFzewxJgoDDFBISKvCgsLcfDgQXm7rKwMW3b/G/FDcuV98UNykTFjjcMxFRUVQUsO8gtucXj+vkWLkHX3euQX3BKU5ycibTFBISJV7LtIxo4dq9hlIkR5PyaQ7J8/Nnu4wzYRhRcmKERERKQ7PUIdABGRP4xGo+PU9vUNoQuGiDTHFhQiIiLSHSYoROQXS32Ly764nhzzQUTaYBcPEfnltPUcAOBXYwZi4hVZAIAYtGGw5ylJiIhUCUgLisViwdSpU5Gamoq4uDiMGDEC//znP+X7RVHEE088gYyMDMTFxSE/Px+HDx8ORChEFEDt51rx+M3DMMKUjIHJPZCZHBeyWFKM0YotOEotPUSkf5onKKdPn8bo0aPRs2dP/PWvf8XBgwfx7LPPIiUlRT5m+fLleP7557Fu3Trs2bMHRqMR48aNQ2trq9bhEFGEyEqOw+Y5Y+TtX+ReAKCzpYeIwovmXTzLli1DdnY2Xn31VXnfwIED5Z9FUURJSQl+/etfY8KECQCA119/HWlpadi0aRNuv/12rUMioghh34JzYd9eAOpCFwwRdYnmLSgffPABrr76atx6663o168frrjiCrz88svy/RUVFaitrUV+fr68LykpCaNGjcLu3bsVz3n27Fk0NjY63IiIiKj70jxBOXr0KNauXYtBgwbhb3/7G+655x7MmTMHGzZsAADU1tYCANLS0hwel5aWJt/nbMmSJUhKSpJv2dnZWodNREREOqJ5gtLe3o4rr7wSixcvxhVXXIGZM2firrvuwrp16/w+58KFC9HQ0CDfKisrNYyYiIiI9EbzBCUjIwNDhw512HfppZfi+++/BwCkp6cDAOrqHPuG6+rq5PucxcTEIDEx0eFGRERE3ZfmCcro0aNx6NAhh33ffvstLrigY0T9wIEDkZ6ejm3btsn3NzY2Ys+ePcjJydE6HCIiIgpDmlfxzJs3D7m5uVi8eDFuu+027N27Fy+99BJeeuklAIAgCCguLsbTTz+NQYMGYeDAgXj88ceRmZmJiRMnah0OERERhSHNE5Qf//jHeO+997Bw4UI89dRTGDhwIEpKSjBlyhT5mIcffhhWqxUzZ85EfX09xowZgy1btiA2NlbrcIiIiCgMBWSq+/Hjx2P8+PFu7xcEAU899RSeeuqpQDw9ERERhTkuFkhERES6wwSFiIiIdIcJChEREekOExQiUsVoNEIURYiiCKPR6PcxgWT//LFx8UF/fiLSDhMUIiIi0h0mKERERKQ7TFCIiIhId5igEBERke4wQSEiv4ntNvnnnTt3wmazeTiaiEg9JihE5Je9n/wVNevvlbcLCgowYMAAmM3mEEbl6sixJnxtaZBvlvqWUIdERCoIoiiKoQ7CV42NjUhKSkJDQwMSExNDHQ5RxDGbzZg0qQiA49eHIAgAgNLSUhQWFoYgsk5/O1CLWW/sc9kf19OArQ9cg6zkuBBERRTZfPn7zRYUIvKJzWbD3Llz4ZycAID0/05xcXHIu3vcJSAtbTactp4LcjRE5CsmKETkk/LyclRVVbm9XxRFVFZWory8PIhREVF3wwSFiHxSU1Oj6XFEREqYoBCRTzIyMjQ9johICRMUIvJJXl4eTCYTAEHxfkEQkJ2djby8vOAGRkTdChMUIvKJwWDAqlWrFO+TqnhKSkpgMBiCGRYRdTNMUIhIFUt9izyXyOBR1+P2hc/B0CvV4RiTyaSLEmMiCn89Qh0AEemfpb4Fo5d+7LT3YmTMWIOqVZMBAGVlZbjhhhvYckJEmmALChF55W7eECGqMxkZO3YskxMi0gwTFCIiItIdJihERESkO0xQiIiISHeYoBAREZHuMEEhom4pxRiNuJ7Kg3Yt9S1BjoaIfMUyYyLqlrKS47D1gWscKpA27bfglU8ruJoxURhggkJE3VZWchyykuPk7a8tDSGMhoh8wS4eIlJktVohCAIEQUBzszXU4RBRhGGCQkRERLrDBIWIvKqpb1Xc724QKhFRV3EMChF5dbq5Y1Dpr8YMxMQrsuT9MWjD4OWhioqIujMmKESk2sX9emF4VpK8bbVybAoRBQa7eIiIiEh3mKAQERGR7jBBIaKI0n6uFT8bdQEEQWAXFZGOMUEhIiIi3WGCQkRERLrDKh4i8pvRaIQoiqEOg4i6IbagEBERke4wQSEiIiLdYYJCREREusMEhYgi1oHqBljqW0IdBhEpYIJCRIqq7f5wHz3eFMJItJNijHbYLlq7G/nP7mCSQqRDTFCIyIWlvgXXPfOxvL3uT3+B2G5z+QMfbrKS41z2tbTZcNp6LgTREJEnTFCIyMW775aiZv298vbx0kWwrJuBb3Z9FMKoiCiSMEEhIgdmsxnzZ02Hremkw37bmROYP2s6zGZziCIjokjCBIWIZDabDXPnzvU4+VpxcTFsNlsQoyKiSMQEhYhk5eXlqKqqcnu/KIqorKxEeXl5EKMiokjEBIWIZDU1NZoeR0TkLyYoRCTLyMjQ9DgiIn8xQSEiWV5eHkwmEwRBULxfEARkZ2cjLy8vyJFpR2zvHD/TWvm1wzYR6QcTFCKSGQwGrFq1yuMxJSUlMBgMQYpIW1vLPlAsn95a9kEIoyIiJUxQiAhWqxWCIEAQBIwbNw4rX9wAQ69Uh2MMCX2w8sUNKCwsDFGUXcPyaaLwwgSFiFzkF9yCjBlr5O2+RYuQdfd65BfcEsKo/MfyaaLwE/AEZenSpRAEAcXFxfK+1tZWzJ49G6mpqejVqxcmTZqEurq6QIdCRD4Qojq7cWKzhztshxuWTxOFn4AmKJ9//jlefPFFXHbZZQ77582bhw8//BDvvvsuduzYgerq6rBtNiYi/WP5NFH4CViC0tTUhClTpuDll19GSkqKvL+hoQHr16/HypUrcd111+Gqq67Cq6++il27duGzzz4LVDhEFMFYPk0UfgKWoMyePRs33XQT8vPzHfbv27cPbW1tDvsvueQS9O/fH7t371Y819mzZ9HY2OhwI6LAsdS3uOyL62kI29WMvZVPA+FfPk3U3QQkQXn77bfxxRdfYMmSJS731dbWIjo6GsnJyQ7709LSUFtbq3i+JUuWICkpSb5lZ2cHImwi+sFp6zmH7dJ7crD1gWuQlRwXooi6xr582l2SEs7l00TdkeYJSmVlJebOnYuNGzciNjZWk3MuXLgQDQ0N8q2yslKT8xKROsMyk8I2OZEUFhaitLQUmZmZDvsNCX0wb+lajoMj0hnNE5R9+/bh2LFjuPLKK9GjRw/06NEDO3bswPPPP48ePXogLS0N586dQ319vcPj6urqkJ6ernjOmJgYJCYmOtyIiHxVWFiIgwcPytuPPLcBWXevx8hrbwxhVESkpIfWJ7z++uvx1VdfOey78847cckll+CRRx5BdnY2evbsiW3btmHSpEkAgEOHDuH7779HTk6O1uEQETmw78a55IqREGr/G8JoiMgdzROUhIQEDB8+3GGf0WhEamqqvH/GjBmYP38+evfujcTERNx///3IycnBT37yE63DISIiojCkeYKixnPPPYeoqChMmjQJZ8+exbhx47BmzRrvDyQiIqKIEJQEZfv27Q7bsbGxeOGFF/DCCy8E4+mJiIgozHAtHiIiItIdJihE5LBI3s6dO9HORfOIKMSYoBBFOLPZjKFDh8rbBQUFuH/iaDQf5tITRBQ6IRkkS0T6YDabUVRUBFEUHfafOlYLbF4RoqiIiNiCQhSxbDYb5s6d65KcdFDaR0QUPGxBIYpQ5eXlqKqq8nrcJ598AqPRGISIQufIsSZ8bWmQt1OM0WE/tT9RuGOCQhShampqND0uHKXEd6zO/MqnFXjl0wp5f1xPQ1gvjkjUHbCLhyhCZWRkaHpcOMpIVl7QtKXN5rKiMxEFFxMUogiVl5cHk8kEQRDcHpOemYW8vLwgRkVE1IEJClGEMhgMWLVqlcdjHlm0xGFxPSKiYGGCQhTBCgsLsfLFDTD0SnXYb0jog74TH0V+wS0hioyIIh0HyRJFuPyCW5BxIAZVqyYDAPoWLULcwCsgRHXPlhOj0SiXVttX7hCRvrAFhYgckpHY7OHdNjkhovDBBIWIiIh0hwkKERER6Q4TFCKKaO3nWvHdsvH4btl4tJ9rDXU4RPQDJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3eFEbUQRyFLfIi+Gd+RYU4ijCZ0UYzTieirP+WKpb8HwrKQgR0REEiYoRBHGUt+C0Us/9npcXE8DUozRQYgodLKS47B5zhgMXt6xXXpPDv5+qB6vfFrB1YyJQowJClGE8faHt/SeHMTHG5FijEZWclyQogqdTLvXOCwzCUdPnw9hNEQkYYJCRA6GZSbBaDSGOgwiinAcJEtERES6wwSFiBAVHYuvquohiiJbT4hIF5igEBERke4wQSGiiGaz2eSfd+7ciXa7bSIKHSYoRBSxzGYzhg4dKm8XFBTg/omj0XxoVwijIiKACQpRROIKvh3JSVFRESwWi8P+U8dqcXzTYuz95K8hioyIACYoRBSBbDYb5s6dC1EUFe7t2Pf6c79x6P4houBigkJEEae8vBxVVVUejzlZV4PXzX/F15YGWOpbghQZEUk4URsRRZyamhpVxz38xk4Y9wmI62nA1geuiYiZdYn0gi0oRBRxMjIyVB1n6JUCAGhps3FtHqIgY4JCRBEnLy8PJpMJgiC4PcaQ0AcxpmFBjIqI7DFBIaKIYzAYsGrVKgBwSVKk7d7Xz4QQZQh6bETUgQkKEUWkwsJClJaWIjMz02F/WkYm+k58FPFDcuV97edaMcKUDEEQYLVagx0qUURigkJEEauwsBAHDx6Ut8vKyrBl978dkhMiCg0mKEQU0QyGzm6csWPHOmwTUegwQSGKMEpzesT1NCDFGB2CaIiIlHEeFKJuzmq1olevXgCApqYmnLaeg9jeOUPqwitF3HTjGM7xQUS6whYUogiz95O/omb9vfL2vdNuxejLL4XZbA5hVEREjpigEEWQ999/H88tuAe2ppMO+y0WC4qKipikEJFuMEEhiiAPP/wwpMXw7EmL5hUXF3OBPHC1ZyI9YIJCFEEsFovb+0RRRGVlJcrLy4MYUfg5UN3ABQSJgoCDZInIgdqF9CJV0drdiIqO5QKCRAHGFhQicqB2Ib3uKsUYjbiejnOh2Fc9tVZ+DbHdxgUEiQKMCQpRBMnKygKgvECeIAjIzs5GXl5ecIPSmazkOGyeM0beLupT41D1dLx0ESzrZqD50K5QhEcUMZigEEWQ5cuXK+6XFsgrKSnhTKoAMu26bZ779RyXqifbmRM4vmkxtpZ9EOzQiCIGExSibs6+KiclJQXFi1+AoVeqwzEmkwmlpaUoLCwMdnghZzQaIYoiRFGE0Wh0PUB0rXqSLFu0kFVPRAHCBIWom7FarRAEAYIg4M0338TQoUPl+woKCvB6yW+RdM10eV9ZWRkqKioiMjnpqtpqC3r06MEVjokCQBBFD/8e6FRjYyOSkpLQ0NCAxMTEUIdDpCv2U9sLggDXj7gA+7lQmpqalFsOIpj9NVSD15BIHV/+frMFhagbU/7/I+z+J9G9A9WcF4VIa5onKEuWLMGPf/xjJCQkoF+/fpg4cSIOHTrkcExraytmz56N1NRU9OrVC5MmTUJdXZ3WoRAR+cXXcSVFa3cj/9kdTFKINKR5grJjxw7Mnj0bn332GT766CO0tbXhhhtucOijnTdvHj788EO8++672LFjB6qrq8O+/9u+35/90UThy2w2O4zbUYvzohBpS/OZZLds2eKw/dprr6Ffv37Yt28fxo4di4aGBqxfvx5vvvkmrrvuOgDAq6++iksvvRSfffYZfvKTn2gdEhGRKmazGUVFRW66xjoZEvog+Zpf4OTmFUGKjCjyBHyq+4aGBgBA7969AQD79u1DW1sb8vPz5WMuueQS9O/fH7t372aCQqQlQfBYJkudbDYb5s6d6zU56TPpCcRfeBXE82046fFIIuqKgCYo7e3tKC4uxujRozF8+HAAQG1tLaKjo5GcnOxwbFpaGmpraxXPc/bsWZw9e1bebmxsDFjMRBSZysvLUVVV5fU4QYiCEGWAiLYgREUUuQJaxTN79mx8/fXXePvtt7t0niVLliApKUm+ZWdnaxQhUfe29PmXXCZlMyT0Qer4B0MUkX6pXSTRZq0PbCBEBCCACcp9992HzZs345NPPoHJZJL3p6en49y5c6ivr3c4vq6uDunp6YrnWrhwIRoaGuRbZWVloMImCnv2FSiJyclIv/P38nbfokXIuns94gexK9WZ2kUSDcZkAMoLCBKRdjRPUERRxH333Yf33nsPH3/8MQYOHOhw/1VXXYWePXti27Zt8r5Dhw7h+++/R05OjuI5Y2JikJiY6HAjIlfOFSj3TrsVta/eJ2/HZg+HEMW1dpTk5eXBZDLJ6xK5E5N1KZoP7VJcQJBr8xBpR/MEZfbs2fjjH/+IN998EwkJCaitrUVtbS1aWjrmB0hKSsKMGTMwf/58fPLJJ9i3bx/uvPNO5OTkcIAskR1fS9elChSLxeKw33mhOwCIio7FV1X17tefiUAGgwGrVq0CAI9JSst/P8fxTYsVFxCcN/PnnGqASCOaD5Jdu3YtAOCnP/2pw/5XX30Vv/jFLwAAzz33HKKiojBp0iScPXsW48aNw5o1a7QOJajsm9V37tyJG264QXerwlrqW7zO05BijEaW3UquFB7UVqCwG8KzwsJClJaWYs6cOQ6JXkZmFmqqO7YN/9zo9TxcQJCo67gWjwbMZrPLF5rJZMKqVat0MwGdpb4Fo5d+7PW4uJ4GbH3gGiYpOmC/Hoy3tV62b9+Oa6+91us5+xYtQvxFVwMANt8/BsOzkrQJtpuRvmOAjsUUR48eLW+rUVZWhhtvvDFQ4RGFLa7FE0TumtUtFguKiopgNptDFJkjtTNccjbM8MQKFG3Zt36OHTvW59ZQd1MmEJF6TFC6wFOzurSvuLiYzb0UcL5WoFBguatIJCL1Aj6TbDhyHqvhblyGt4mdRFFEZWUlysvLXcbkkCulMTIcE6NOXl4e0jIyUVdT7fG4mKxLgxRR92I0GiGKImw2GwYMGACLxeJxvM/o0aODGB1R98QExYnSWA134zLUNqurPS6SuRsjwzEx6hgMBiz4zVLMm/lzj8dJJcZxPQ1IMUYHI7RuRar0KSoq8njcf+qaEH9GZIJN1AXs4nGiNP7C3bgMtc3qao+LZO7GvXBMjHr5Bbeg78RHXWeOtdsuvScHm+8fw6SvCwoLC7HyxQ0er3PR2t0Yv/pT5D+7A5b6lmCHSNQtMEFR0H6uFd8tG4/vlo1H+7lWt8d5m9hJEARkZ2cjLy8PgO/zWrhz2HJCPs/ew9X42tLAL8FuyLl0Xc1YpvghuciY0Vmy37doEbLu+YM858nIQZkYnpXE5KSL8gtucbnO9jP2SjPLMsEm8h8TlC7wNLGTtF1SUqLpfCiW+hZc/+wOeZv/qXVPzjPCFhQUYMCAAaqqwuxniuXMsYFjf13bW5scZuyVZpZtPrQrFKERdQtMULpImtgpMzPTYb/JZEJpaanm86CwK6T786V03VLfgq8tDfja0oAjx5qCHSr94OTmFYozyx7ftJjT3xP5iYNkNVBYWIj8/HyHiZ3UziTry2Rc1P15K10XBAGTJk0CAHxbdRz/b/WeYIcYEaSqHXdSjNGI66muZWrZooW4/5d36G5maSK9Y4Kika5O7ETu2bcMKFVFdKcp/NWUrktON7PFLFSykuOwec4YDF7u/djaagunGiDyAxMU0r3iP30p/+xcdtzdpvBnSXr4yPThvcTfK5HvOAZFgf2CatJofC34U5WhJFDxhQPnsTbdbQp/lqR3T/y9EvmOCYqTrWUfoGb9vfK2NBrfn4Fu9mXFb775pmJVxvvvv69JfIGsFtCqPFqr83RnakrXSYc8/F7SM7PkqQaISD0mKHbMZjPmz5quOBp//qzpXVr4b+rUqYpVGVOnTtUkPlYLdA/2pevUPTyyaAnHpBH5gQnKDzxVT0i6svCfpwUFtYpv2aKFXJiwG3A7U2lCHyx5/iWvj4+KjsUFj2zGBY9sRlR0bKDCjHhSpY8oinjOze+r78RHkV9wS4giJApvTFB+4MvCf1qyTzg8jUvxFh/QWS1A4U9pptKsu9djbP44ed++PbsiavyRnrn7fcUPycWRY03yXDWcTJFIPVbx/EAPC/8VFBQgIzMLDy9a4vBfVwzacO2116o6x94D/0WfQVe47D9yrAnt51pR+VzHImfZ80rd/nftXNab3NOXV6EttTF3R84zwrYc3oOJ194t77t32q0wJPRB7+tnIn5ILoCOaqWWNsekhQsDBoe7GXydq9DWTbsKqXa/j3ApgbfH+ZsoGJig/KCrC/95m9hJrZpqC+bN/Dn6TnxU/qMTI7apfvzKfxzDmspPuxSD8xfqB3df3aXzUdc1H/4MJzevcNkvjT+S3i8tbTaUTL4cF/frJR8Tjn8Au6uWNhum/2Gvw75wKYEnCjZ28fzA14X/1OjKeJBT216Sm+8d/iP2UC1gSOiDGNMwv59TSUubDSfOdDZLd6U8Wqsy60hUv/1Vj/fbv18u7tcLw7OS5Bv/8OlbuJTAEwUbE5QfqKme8GXhP+fF3nxlO3MCZ6sO+PSY3tfP1HxhuOZDuzDx2lHyti+L1tnryuJ34c5+vRx/xyI4V2653O/H+4WISM+YoNjxVD2x8sUNHhf+s/8jVPLyG4qLvfnK1nTaZd+8p593Wy0gdQkBHWM3vls2Ht8tG4/2c61+PX/zoV04vmkxjtU6jrtRWrTOE18Wv1PDftChHhfIs5/v5bDlBEYv/RgFz27FCFMyRpiSUfDs1oCsPt1Wf0zT8xEBygl2NQf7UhBwDIqT/IJbkHEgBlWrJgPoGI0fN/AK5Bdc4/Yx9tOti+02WNbN12Q8iqFXinxOyR/3n0T6nb+HZfXPHOJT03LiPAOtp8eJ7Tac2qZc0iotWldcXIwJEyZ4bFVSs/hdcXExPvz0S9Ux24+R0Tt36+VIzfpadr8YjMmanYsIcL+URLSt85+enTt3ql4clcgXTFAUuBuN7459//HZqgOwnTnh9Tn69OmDkydPuk1kpPEkzYd24dTWF+X9x0sXObSgqIkPgPJ5nCpA7Hl7HfZl154WQVNbvv3Fnl0AHMfX+BpzdyPNZ9KR9M7w+r6Kybo0SJGRsxRjNIxGIy54ZHOoQ9GU0tiY5kO7UGX3uSwoKIDJZMKqVas8tjIT+YoJisaUumWUTJ48GWvWrIEgCIpJSu/rZ6Ll8B4c37RY4TlcxyM4l+Pa81YBIrEv41X7OtyVXduXIapx/FgdgPTOmH/oXnIXs32XllIpshblyUqrJMegDYNNfQFoW17p/Fz2XVdClAG9r5+peD3saT3+iNTLSo7D1geucfkdznnjM6/vw/ZzrRhhSgagj5Jd+8/unm8du2XdfS6l7trS0tKgJCnhVuas9F3C6jrvmKBoTOqW8eamm27Cddddhzlz5jiMzZBaCOIGjYJl3Qyv53E3UZf9/tNuumoAdFQF/ZAgLS+6DAs++LYjDpWvQ6tF0Pr2SwMqO3721L0kObXtJcQNGhWwP8rumrbVlHzbVyd1TKYWrdhV5e257MUPyUXfiY/i1NYXHRLUKGNvtFtPKZ6XgisrOa5b/8HRqts30rj7fLO83DsOktVYjGkYDAl9vB43evRoFBYW4uDBg/K+NW+8K88+qbar6KzlG5d9zYc/c1hQsL2l0f0J7FpvLuzb+V+It9fhS9l1VlaW1/LtK0d1dtmoee2BrlpxV/bpPAmaM+dqpXun3Yqq1VNR/WJnsikt8PjGW+/ga0sDPq84pSqm+CG5DrOVJubchii763q8dBGq183A7o/LVJ2PyBe+dPtSJ0/fJSwv94wJipN+CTGK+9VWXEjN8d5I/2HY/6dx1ahcuUVAbReLzVrvsu/k5hVey1K9UfM61JZdL1++vOOcTkmKtF1SUoLaM50fVNWvXeVxweKuWqm99QzaWx2rjWxnTmD5QzNx3X3LfRr0a99i1Lj7HZx3Xjiy6SRmTZ8SEeXbFFxd7fYl8hUTFCf9EmOx46Gfytu/yL0AAPB5xSnV81hIzfEu5cB22weqPZfrqe1ikSo3tFiTxXltF7evQ0XZdeXJzj/IDe3RWLH2VfRNS3c4JiMzS+6zlv6T+NWYgVg+bayqeE98+Azaz7W6dJ+I7TbFffaUSicPW07I5cHNzVbF57Q/j/1kc2oWc1RyfNNi2FqVn8sf0vN3ZWFL0o7z+9DWanUp/3d+TzU2NsrvQ6u1471hX7ou7fOHmvO4dlF2bPva7atVzN5iVDvpoz/xKH1PdHVOo2BwjluPMarBMSgKBmakyl/0fztQi9f/WYdXPq3AK59WyMd46z+MH5KLmAt+5FKuLP0HfNv6/QDcj2mQulg8NakaEvogbuAVLtUu/lJa28Xd6/BUdv3yG2/hnnvvdzlv0tifA39ZKZ8nYfDVGHXddQ6PvbhfL9x61Y14wmSCxWLx+ge/+fBnDrOsHi9dhKjYBACiwz5DQh9sHbQSw++a5vf4EufrbF+90Lt3b6+LObpz1vIN4i/SbjkBtRVWFFhbyz5Azfp58rZzBR6g/J7KysoKWozOzGYz5syZI2/bfyfEDRrl8TtJEASYTCafZtvWIsZAVRGpGRsG6G8siVLceotRLbageOHuF6qm/1BNubK7MQ1quljsK3262qUjkapkmg/tcohF4q2s2Ww2Y9b0KS7x2M6cwKkfkhPpPK025f5ZNbP6SpS6s9x1qcyfNR1ms9mv8SVS9YLzc0nVC++//76qeJUoddNpgU3toWM2mzF/1nTXz4HddvPhzxTfU9XV1UGJ0Zm7LkrpO6Hl8B6330n23bWBHCCr9aSPnqgdH6K3sSRKsegtRrWYoOiYpy6Wh595CR+tegBRe14LyHPbr+2ilr/dHEo8zeqbOv5Bv8/rT9eHt+oFANi4caPfMQVqgjWtKqzIN2o/B/Xb/6C43/5xweqmUxOzVDmn9J1kMpkCXmLsbdJHgF2b3Q0TFBWUpo2X5i7wNF5BjajoWLz3RRW+qqpHdZPrB8+5cqNv0SJk3b0e0352G04f/TfqagLz35btzAm0VOxXdaylvgV7D1ejR48ePnVz2F/D1pZmAEBrS7PcT5z70+sVX3vcRT/27cX8QOr62LVTudnWPiGz73tXU71w/PhxpKSmuq1W8hxXu+pkMCo6Fv0fel+zCivSnrfJCSW2Ju/VW2vXroXNZtNsoU2l81itVlWfXalyzvk7qaysDBUVFS7JifNzfX+yyWX8xt7D1fLnfe/hao9jPNRO+uiuiiiYi5Xaj3eR/j5osfxIpOEYFB3wVsXhrosl0E34aroepP7OQHzgUuKjER8TLW/HZg9Hy+E9XR5vc+L4MQAXOOxzHgtg3/cu2rzPfQIA5weOgXjS966eE39+SnGG3LieBsVuJ08TtwWrqZ3c0/JzuWDBAixZusxhEXN/x1y4G7uxbNky1eeQKnnsv5PGjh3r8l5Tei6l97j990bR2t2KE9lJ4yfUXlel44I1boW0xRaUAJKmKr/gkc1+zWTqjS9N+P50jUhdD55ehz/9mmk/W6zqemQmx2HznDHydlGfGk3G2/Tp289h2934Eqnvve2UulaqeDfN31GxCYiK9TyrrtLYn5Y2G0omX47N949ByeTLHZ/LTfdfMJrayTOtu9Ya6k+j/rRjia+WC3ZOnTpVdSxqKnm8jWWxf4+rIY2fUHtdnY8L5rgV0hYTFBW8lbE6l+cGS15eHkwmk9duhT6TnuiYAG7QT+R9niZPk/jS9aD2OGmNIefH/Gf/XojtNrQ7NcOmJXS2oGxc7Xmqd7Xa2ztfl5pZa8/8a4tLIuBMel3Ozd+Pr34Dpvv/iMxZ61XF5jz25+J+vTA8KwkX93NNcNQ2tVNwqf1cdoUoihABzL5/Dv71/SmPpaTexm6oHTNm/9l1R+1YFvvPn0RpOgB7eXl5SMvIdHu/UtdmqMatKJVq+/JaqQO7eLxQKhV0LmNVKs/1ly9ryEjVLkVFRR7PGdf/MghRBojo7KpYvnw5pk6d6nYtIMB91wPQuVbMkWNNPpU5975+JoQog8tjls2bjqjYBMx5sbOp2Lnk8litNk3n902fLL+uqLheXmetbVfRYiO9LsCx+fu1o3GIijaonpLfduYEvn9mgsvvPsUYrdjlY98FptTUTsFn/7n09PnqMlFEbbUF4xa8iNj+lyFGbMO3y/8XgOP6NGrHxHhj/x53R81zSe/x1PEPukwR4O37Zvr8RVj+kHIlkTQGpbW1VfVr16okv7nZCkFIBtAxYP7hhx+W77t32q1epz4A1K/XE27rEHUFExQPpFJB5y+Y9tYzLscqLWKnFamLRYlU7fLg/HkOXRSGXqkuXRZR0bH4qqoew7OSAACxsbEuawE5c/e6pHEz7hYPc2b/xePuMe2tZ2B1GspSXd0xiG7u3LkoKSnx+jwAkJraMY/NqVPuByFKryvh6gmqzglA/kJ1uM5BWl1ZaTE64IfFC5cH9KnJD4WFhSgtLfX6+dKCNC7EXZl8V8fEKL3Hnb9L/HkuTwuYuvu+ATIV16RS+r7zJR4txw1NnTpV9d+M+bOmo3+qEaOuu5Hr9ShgF48b/pbM+lOeay/W7h8Utc2At95ahAvv7myN6Fu0yKHpXxLX04AUY+d/3M5rAXmi9LrUdI8Adl1MQ3JVP0Z+Dh/LeJcuXYq6ujpUVFR4PxiA9cAnqmOJiu2F9Dt/L29LVUWBSE6UfvdZyXEYnpXkcMuM0C+ucODL56srpHEh7mY57sqYGPvPrhpajb/x9D2qVNlo/7n057V3NW77Lh1f/2YUFxfjRKNy95zS/CXBrEYKNSYobvjbLNqVReyaD+3C0XWz5O3jpYtg3XC310FlWclx+Etx58yumxffhT/P7hxcWnpPDjbfP0YxE1fbJaD0utQuaCgIUXLTsNrH2LMv4/XmnnvugcFgUP262lsaERWXqOrY46WLUPvqffK2t0nrJHE9fe92kRYU3Fr2gc+PJf2wfx8GYkyKNC6k+dAuhwVCCwoKkJXdHyUvv4GUCy9DWkamX89v/9n1RJpavSvPZc/b96h9TO2tTQ6fy4KCAmT3v0DVa9eiJL/50C5MvHaUX4+Vupg+/Ns2t8ccOdZZnv3yG285LEZq/1qlY6TuME/nCZfp79nF40ZXmvz8WcTOXbfHqeO1EJ2aPJ1bQgA4/Cc9LDMJRqNR875v6XXZj5NR9Ti7cuWuLPA36baf4ZW1v/d4jD9jMIzDrsWZf6orD1ZqRlYaN2TfJdfSZsPz036Cix+ox9ayDzBv1nSHVaTdPpddE7C7ga+B+D1T4DiPSenqGBX72aSd1dVUY97Mn/t9bqDjsxvX0wCr1eryHrfUt2B4VpLL1Orto34BUUW3rzd1bz2K7HmlAOBxXJ5SV1FNtQXzZv4cfSc+CniJx9eSfPvPe+r4BxWf31erPvwcxqHKy4d46063f62eWrqcp7NwN2ZJT9iC4kZXmvyWTxuLzfePcbg5l4na8zZTqSAIMOzdgPfvdd8SEgxKTcmqHmc3U6raBceUTLltEl56/U30S3f83WRkdm3tEnflwd74UokgVeMU3zUNzynMkOsJZ8fsHv74xz8iM9OxCiUtI9Oh/D91/IOqytQNCX3Qd+KjiBs0SnWXaWLObTD06u1TzAZjMlrabFj6v53VOz9NOgGxvbPrwbkLwl0JvD9aK79G+/k2h21fvn88zX6rZtFTJfbPf9qH7mpPvH0vquka93V4QfPZzt+bXruKmKC44U+poNRc+PPCG13GCiiViUrUzFRaW23B6aP/xvCspIAlJ55eq7umZDVisi7t/PmHRRB9jUtqhr1r2s9w+NB/5PvKysrwn2+89POreF3O/dpqnLV843I9pK4ZT91y+QW3qH4ub7NjUviYMGGCw5iUsrIybNn9b4fy//hBP3EZX+Fcpm4/9smXLtPG3e84NNwl9/aeQMRkXYrmQ7sw77bORT3fWHQ3LOtmYO8nf3X7OH8+T0qOly6C5fdTHLfXzUDz4c9UPd7d7LfSNcwvuMWneJw/7+0tjT49Xoma8m01v2dfhhcodQkOGDBAd3PCMEFxw5cF6+z5M4On2m4PT91OUlO/KIo+NdVJj/vzn/8MwH2S4u/ChH0nPgpDbGc8ahZBtKc0M6r99R07diwSExNdXrv96/KUYrorD1aj9f/2e5zgzV2SkmJ0nCFXDS78F56cP5fO712DweAyEaLSzNHuZpP2tcu03dpZ2XbnA4s8Htt34qM4+92/cHzTYpcSf9uZE3huwT0oefkNt2MefP08uSW2uzz3yc0rVFfgKc1+K11D53EZ9jfn1+VuQseuUlO+rfb3rOY4bwuf6ilJYYLigbsF69w1ufrTXAio7/YI5OJvUlmkcxN0ZpYJz730us8LE0pN0Ep9ou6agKNiE5CU7HgtujozqrdFB49vWuz32hjeKoDcNblmOc2QqwYX/iMlXekyffdYP7ddH+q6j0Q8+MB8zH1rH4DgrzWjtgLvxIfPuI2n+E9fYvzqT+VbwbNbMcKUjBGmZMx5o7OVxtfqQ2dRsQmIT3AsyZau+/FNix3WeLO/htL2iQ+fUfU8XekqkhLpSZMmobGx6y1DWuAgWS/yC25BxoEYVK2aDKCjaTBu4BVoP9fqsi+/QHmQkzdSt4e7JjxBEGAymQK++FthYSHy8/ORlNTxQSorK8MNN9wAg8GA7du3q16YULoenv4riB+Si5gLfuRyDd/+5RXIvbS/y/Pb83VgqKffoaS18mvEDbwCFzyyGWK7DZZ1Mzw2qUbFJXpt3pWaXGP7X+Zy36CsPhBFETabDQMGDIDFYlF8TcH63VNwuLx36xtcjlGa98jdXEjevjs8aTv2X7efQyHKgNbv/626WyG2/2UuY7F8+Tz5Q6rAU9PNYh+PJ0qvQYgy+FV9KOkz6QnEX3iVw3fbmjfexe/2nodl9c8cnsv5+aMzhsjb3l6rVl1FQMcilQ8++GDIJ35kC4oXzs3x3ppcPZ3HXbmpp26PYC/+ptQEDfjWxRCbPVxVF4bSNXT3/F3l/Fy2is/djh1x+H246fIyDrtW1fN6a3K170p07l7jwn/kja9dpvY8dX3Y36/mPJ7GYnUlRm/Ufg7VjA1Teg3VPzymK9WH0kze9p/hxvp6h9Lo46WLULV6KqpfnOEYs934G2+JmH1Xkbu/NWpfx4IFC3QxJoUJihfOzfGl9+R4rMjxdJ6tD1zjtrJH74u/+drFYL/QndLNn2uolebDn8FS+rTHsSPS76NfWrrDMXLX1SB18x6oaYJ3172ml9896Zu/VTPe3ptqu4/aTlu8jsVyF2NXK318qcDzNDbM3biM82dO4Pj7SzB5sG+dDd7iWTBnpstztbeeQXur03gep/E3is+l0J3ubqFRX7oE9TAmhV08KkjN8ZKvLQ0ep593Jys5zqECx3l9lfghuUi++EocXtFRY++uiyOQ3HWfSFVN7roigI4PStbd6yFEGRDX04AfD+ztseJI6RrGxwdmXg/puaTmZk+k0sT4Ibl4a8ndcrPsG+9swtNf9kCrraMp2FvTupomV4mn7jWKLEprLjmztVodu2UGjUKGXVeNL90Bzp9D6fnVdB9F9UpF05d/8xir9HlKHZ7ntjvJn24g6TUIUQaH8yb3TkX9KfcDWaV4hCiDw3V0RwDw/p/e8Pr9J3HXxR0fb8T58+cxYMAATdZGAly7iuyfW5raAHB8v6jtFpOmuCguLsaECRNC8l3EBCWElNZXiUEbBv8w74+eFn/zujChIGDFsyvlcThKi1zpgS/lerH9L3O4/v9bkI9r/1+U/PvaOmhlx1pNgOLEa2pG59sLVPcW6Ze7RSA3/HIkUn+YjPGk9RzufmOffIzzQpvSwnPJ1/yi87zXz/Q4gZin96b03/fF/XrJ73FAeQr3hB/9Dxr+4XkZCunzJPS/DM/cdgV+9kNx5ObFdyE+vrPCz9vnyYHd982m/Ra89PE38l2ekhP7eNpbmlQtciqKIqqqqvCb3/wGixYt8jq5nqcuf60WbpR0dBUtkLftF1wEOlr+nRe89aU0WqvFFP3FBCXEnFtVAN/XcggWtwsTJvTBimdXoviHVTn1zNdyPecWHSMg/76G3zUN/VONLgvC9Ujsi5Tr7vI4868SzgobeZT+SVFK7qVjtpZ9gHnLl7j8AZdKb+2b+qN6RLsuqqdycUvpv29373HpPKKtzcNZ7OL74fM0/IJ+bt/j7p4LQpRDV4fz943Uov3Wnu8g/vcfuOOOO7zG03x4j+rZoyWDBg1SXADSlwVDtZwuoLi4GAvnzHS5nlJX1taCS/BtqlFxwVtfhWqaAyYo5JNbby3Cmm/j5W6ovkWL0Hvw1bj11uu8PDJ07P9LVdsHq/Y4pa6ZYSPz0Nja+R+xXluTSB+U/klROiY9IRo3PvWox9YF++4LTxU6EqXWG6VFRe3f42veeBdL/x0tV/qo4e/n6Y13NuE3n53FkZWdr8HT943asXK+LBJqf+6f/vSnbq+F2nNoZePGjR4Tj98smIcoiJr80xOqaQ6YoJBPpIUJpW6ozYvvQlbfFF3/Abb/L9Vmy8G47b/Hsdoaj2Np1I4dAVy7ZvS4pgWFPzXdA86l7d6qDe27cyRKCbX9e/yqUbkQvt4PQF2Zc1c+T/9bkI9RuS0YvLJj29P3zZFjTbj0ssuQ2i8DJ4+5/49f7RgMiXOpv7troYaasXxq4unTpw+OHz/u8ThvXV1OZwWgHE96ZlbIpjlgFY8flEqG1TbjdwfSoGFRFDFyUKZfyUmwr2FWchyGZyXhR/17Y83vVwPwPGuuNNDXl64ZX2fxJfKF2mZ2+25M51lqlcjdOT/clD7P9u9x+3EjakqIu/p58vZ9I53zlU8rMGHNbkTl3un5/CpLk+3Zl/q7uxZq+DtDuT0RwJQpU7wep1bHjLzuk6VHFi0J2Xg4tqD4QW2/MbkXymsolfU69yVnZpnw0JOL5fU5+DslPVHbzL582lj8ONd1puIjx5pcVrTVglRCbPvHH3DqeK28Pz0zC48sWhLwz5PzOb3Fk5Scgl/epm78iSGhD5Y/swK9e/fGW2+9hYyMDOTl5XXpD7b0/TP7vvtRazf5ZUrvVAgQcepU53IEUVFRaG93HH8z57GnccuoISgpKfE7Bnvxg0Yh1jTM7XglX9cr0lJIE5QXXngBzzzzDGpra/GjH/0Iq1evxsiRI0MZkmpq+o3Js1Bew8LCQkyYMAHl5eWoqanR5IuHKJC8dQ9IXRE/L7xR8X2sVDHkT6ul0nlSh+fhb68sxNGv/qmLz1P8kFy8v+oBnD76b5d4bDYbHsvMQk21xe3jhdgE9J3wCHq0NWPFU485HGsymbBq1SoUFha6rcJypjSuR+n7B4DDvtzcXOzatQs1NTX4pj4KGyriMPLay5F3VRZMJhOqqizw1PrhjX2pdtygUR1Vjk2nYeiVIu8PJUEMUdnAn/70J/z85z/HunXrMGrUKJSUlODdd9/FoUOH0K9fP4+PbWxsRFJSEhoaGpCYmBikiImIQstsNsul/vZf3VJ3pbeJ/Sz1LZq0Wmp1Hi18bWnA+NWfuuzffP8YeR4QZ9J1dC5plq7jyhc3AIBiBYzztXa+Fkq0uD5v7/0eC8xfYWnhCNw+sj/MZjMmTSqCPwmKVCrtbr00e56uoz98+fsdsjEoK1euxF133YU777wTQ4cOxbp16xAfH48//OEPoQqJiEjXpO6BrKwsh/1qZx2WxmJ5Gm+ihlbnCRXpOprcXMf7f3kHnn3qUcWWKmlfcXExbDaby7VQugXi+hQWFmLe0rUwJPTx+bFpGZkuyYnS9PihHlsZki6ec+fOYd++fVi4cKG8LyoqCvn5+di9e7fL8WfPnsXZs2fl7YaGjgW29LLiIhFRsOTn5+Pf//43du3ahdraWqSnpyM3NxcGgyEivxObzjSi/Wyz4v7GRuWB8IDn67hlyxaPFVPSBGZbtmwJWoVLc9MZtJ9txof/PILmpjMAgLqEQci4czXuuKAFFyecx8KFC3HypPvqnZSUFLz22mtIvegy/OyVzx2um/UssLRwBC7s2znwNzk+GglRbWhsVDffjRrSe1RV540YAhaLRQQg7tq1y2H/Qw89JI4cOdLl+CeffFJERzsWb7zxxhtvvPEW5rfKykqvuUJYVPEsXLgQ8+fPl7fb29tx6tQppKamui0V9VdjYyOys7NRWVnJ8S0BxOscHLzOwcHrHBy8zsETqGstiiLOnDnjsjiqkpAkKH369IHBYEBdXZ3D/rq6OqSnp7scHxMTg5iYGId9ycnJgQwRiYmJ/AAEAa9zcPA6Bwevc3DwOgdPIK61NBOvNyEZJBsdHY2rrroK27Ztk/e1t7dj27ZtyMnJCUVIREREpCMh6+KZP38+pk+fjquvvhojR45ESUkJrFYr7rzT8yyARERE1P2FLEGZPHkyjh8/jieeeAK1tbW4/PLLsWXLFqSlpYUqJAAd3UlPPvmkS5cSaYvXOTh4nYOD1zk4eJ2DRw/XOmQTtRERERG5w8UCiYiISHeYoBAREZHuMEEhIiIi3WGCQkRERLoTkQnKCy+8gAEDBiA2NhajRo3C3r17PR7/7rvv4pJLLkFsbCxGjBiBsrKyIEUa3ny5zi+//DLy8vKQkpKClJQU5Ofne/29UAdf38+St99+G4IgYOLEiYENsJvw9TrX19dj9uzZyMjIQExMDAYPHszvDhV8vc4lJSUYMmQI4uLikJ2djXnz5qG1tTVI0YannTt34uabb0ZmZiYEQcCmTZu8Pmb79u248sorERMTg4svvhivvfZawOMMyVo8ofT222+L0dHR4h/+8AfxwIED4l133SUmJyeLdXV1isf/4x//EA0Gg7h8+XLx4MGD4q9//WuxZ8+e4ldffRXkyMOLr9f5jjvuEF944QVx//794jfffCP+4he/EJOSksSqqqogRx5efL3OkoqKCjErK0vMy8sTJ0yYEJxgw5iv1/ns2bPi1VdfLRYUFIiffvqpWFFRIW7fvl388ssvgxx5ePH1Om/cuFGMiYkRN27cKFZUVIh/+9vfxIyMDHHevHlBjjy8lJWViY899phoNptFAOJ7773n8fijR4+K8fHx4vz588WDBw+Kq1evFg0Gg7hly5aAxhlxCcrIkSPF2bNny9s2m03MzMwUlyxZonj8bbfdJt50000O+0aNGiXOmjUroHGGO1+vs7Pz58+LCQkJ4oYNGwIVYrfgz3U+f/68mJubK77yyivi9OnTmaCo4Ot1Xrt2rXjhhReK586dC1aI3YKv13n27Nnidddd57Bv/vz54ujRowMaZ3eiJkF5+OGHxWHDhjnsmzx5sjhu3LgARiaKEdXFc+7cOezbtw/5+fnyvqioKOTn52P37t2Kj9m9e7fD8QAwbtw4t8eTf9fZWXNzM9ra2tC7d+9AhRn2/L3OTz31FPr164cZM2YEI8yw5891/uCDD5CTk4PZs2cjLS0Nw4cPx+LFi2Gz2YIVdtjx5zrn5uZi3759cjfQ0aNHUVZWhoKCgqDEHClC9XcwLFYz1sqJEydgs9lcZqtNS0vDf/7zH8XH1NbWKh5fW1sbsDjDnT/X2dkjjzyCzMxMlw8FdfLnOn/66adYv349vvzyyyBE2D34c52PHj2Kjz/+GFOmTEFZWRmOHDmCe++9F21tbXjyySeDEXbY8ec633HHHThx4gTGjBkDURRx/vx53H333Xj00UeDEXLEcPd3sLGxES0tLYiLiwvI80ZUCwqFh6VLl+Ltt9/Ge++9h9jY2FCH022cOXMG06ZNw8svv4w+ffqEOpxurb29Hf369cNLL72Eq666CpMnT8Zjjz2GdevWhTq0bmX79u1YvHgx1qxZgy+++AJmsxl/+ctf8Nvf/jbUoZEGIqoFpU+fPjAYDKirq3PYX1dXh/T0dMXHpKen+3Q8+XedJStWrMDSpUuxdetWXHbZZYEMM+z5ep3/+9//4v/+7/9w8803y/va29sBAD169MChQ4dw0UUXBTboMOTP+zkjIwM9e/aEwWCQ91166aWora3FuXPnEB0dHdCYw5E/1/nxxx/HtGnT8Ktf/QoAMGLECFitVsycOROPPfYYoqL4P7gW3P0dTExMDFjrCRBhLSjR0dG46qqrsG3bNnlfe3s7tm3bhpycHMXH5OTkOBwPAB999JHb48m/6wwAy5cvx29/+1ts2bIFV199dTBCDWu+XudLLrkEX331Fb788kv5dsstt+Daa6/Fl19+iezs7GCGHzb8eT+PHj0aR44ckRNAAPj222+RkZHB5MQNf65zc3OzSxIiJYUil5nTTMj+DgZ0CK4Ovf3222JMTIz42muviQcPHhRnzpwpJicni7W1taIoiuK0adPEBQsWyMf/4x//EHv06CGuWLFC/Oabb8Qnn3ySZcYq+Hqdly5dKkZHR4ulpaViTU2NfDtz5kyoXkJY8PU6O2MVjzq+Xufvv/9eTEhIEO+77z7x0KFD4ubNm8V+/fqJTz/9dKheQljw9To/+eSTYkJCgvjWW2+JR48eFf/+97+LF110kXjbbbeF6iWEhTNnzoj79+8X9+/fLwIQV65cKe7fv1/87rvvRFEUxQULFojTpk2Tj5fKjB966CHxm2++EV944QWWGQfK6tWrxf79+4vR0dHiyJEjxc8++0y+75prrhGnT5/ucPw777wjDh48WIyOjhaHDRsm/uUvfwlyxOHJl+t8wQUXiABcbk8++WTwAw8zvr6f7TFBUc/X67xr1y5x1KhRYkxMjHjhhReKv/vd78Tz588HOerw48t1bmtrExctWiRedNFFYmxsrJidnS3ee++94unTp4MfeBj55JNPFL9vpWs7ffp08ZprrnF5zOWXXy5GR0eLF154ofjqq68GPE5BFNkORkRERPoSUWNQiIiIKDwwQSEiIiLdYYJCREREusMEhYiIiHSHCQoRERHpDhMUIiIi0h0mKERERKQ7TFCIiIhId5igEBERke4wQSEiIiLdYYJCREREusMEhYiIiHTn/wN1ymumgQVVkgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "m.visualize()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The fit succeeded and the statistical uncertainty in the template is propagated into the result. You can play with this demo and see what happens if you increase the statistic of the template.\n", + "\n", + "Note: the parameters of a parametric components are prefixed with `xn_` where `n` is the index of the component. This is to avoid name clashes between the parameter names of individual components and for clarity which parameter belongs to which component." + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Extreme case: Fitting two parametric components\n", + "\n", + "Although this is not recommended, you can also use the `Template` class with two parametric components and no template. If you are in that situation, however, it is simpler and more efficient to use `iminuit.cost.ExtendedBinnedNLL`. The following snipped is therefore just a proof that `iminuit.cost.Template` handles this limiting case as well." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Migrad
FCN = 92.87 (chi2/ndof = 1.0) Nfcn = 190
EDM = 1.67e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0_n 990 40 0
1 x0_mu 0.4964 0.0019 0 1
2 x0_sigma 0.0487 0.0016 0
3 x1_n 629 29 0
4 x1_mu 0.97 0.14 0
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
x0_n x0_mu x0_sigma x1_n x1_mu
x0_n 1.3e+03 0.000424 (0.006) 0.0123 (0.209) -264 (-0.252) -0.369 (-0.076)
x0_mu 0.000424 (0.006) 3.44e-06 -9.8e-08 (-0.032) 0.000293 (0.005) -1.64e-05 (-0.065)
x0_sigma 0.0123 (0.209) -9.8e-08 (-0.032) 2.64e-06 -0.0129 (-0.271) -1.64e-05 (-0.075)
x1_n -264 (-0.252) 0.000293 (0.005) -0.0129 (-0.271) 851 0.337 (0.085)
x1_mu -0.369 (-0.076) -1.64e-05 (-0.065) -1.64e-05 (-0.075) 0.337 (0.085) 0.0183
" + ], + "text/plain": [ + "┌─────────────────────────────────────────────────────────────────────────┐\n", + "│ Migrad │\n", + "├──────────────────────────────────┬──────────────────────────────────────┤\n", + "│ FCN = 92.87 (chi2/ndof = 1.0) │ Nfcn = 190 │\n", + "│ EDM = 1.67e-05 (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 │ x0_n │ 990 │ 40 │ │ │ 0 │ │ │\n", + "│ 1 │ x0_mu │ 0.4964 │ 0.0019 │ │ │ 0 │ 1 │ │\n", + "│ 2 │ x0_sigma │ 0.0487 │ 0.0016 │ │ │ 0 │ │ │\n", + "│ 3 │ x1_n │ 629 │ 29 │ │ │ 0 │ │ │\n", + "│ 4 │ x1_mu │ 0.97 │ 0.14 │ │ │ 0 │ │ │\n", + "└───┴──────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", + "┌──────────┬───────────────────────────────────────────────────┐\n", + "│ │ x0_n x0_mu x0_sigma x1_n x1_mu │\n", + "├──────────┼───────────────────────────────────────────────────┤\n", + "│ x0_n │ 1.3e+03 0.000424 0.0123 -264 -0.369 │\n", + "│ x0_mu │ 0.000424 3.44e-06 -9.8e-08 0.000293 -1.64e-05 │\n", + "│ x0_sigma │ 0.0123 -9.8e-08 2.64e-06 -0.0129 -1.64e-05 │\n", + "│ x1_n │ -264 0.000293 -0.0129 851 0.337 │\n", + "│ x1_mu │ -0.369 -1.64e-05 -1.64e-05 0.337 0.0183 │\n", + "└──────────┴───────────────────────────────────────────────────┘" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# signal model: scaled cdf of a normal distribution\n", + "def signal(xe, n, mu, sigma):\n", + " return n * norm.cdf(xe, mu, sigma)\n", + "\n", + "# background model: scaled cdf a an exponential distribution\n", + "def background(xe, n, mu):\n", + " return n * truncexpon.cdf(xe, xe[0], xe[-1], 0, mu)\n", + "\n", + "# fit\n", + "c = Template(n, xe, (signal, background))\n", + "m = Minuit(c, x0_n=500, x0_mu=0.5, x0_sigma=0.05, x1_n=100, x1_mu=1)\n", + "m.limits[\"x0_n\", \"x1_n\", \"x0_sigma\", \"x1_mu\"] = (0, None)\n", + "m.limits[\"x0_mu\"] = (0, 1)\n", + "m.migrad()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABAjklEQVR4nO3de3RU5aH+8WdnkAQTEgyQGxMBOVqE2uMFoRFDpc2SU6onnBCv4KEu6xV7AE9FqIottYIsj1wsytF6PYBUQpBKo7bFIKlSUNTz8wBi1aghkKClJCQKSLJ/f+CMM8lMZs/Mnpk9yfez1l4re887e97Zc3vyvvvdr2GapikAAAAHSUl0BQAAADoioAAAAMchoAAAAMchoAAAAMchoAAAAMchoAAAAMchoAAAAMchoAAAAMfplegKRKK9vV379u1T3759ZRhGoqsDAAAsME1Thw8fVkFBgVJSum4jScqAsm/fPhUWFia6GgAAIAJ1dXVyu91dlknKgNK3b19JJ55gZmZmgmsDAACsaG5uVmFhofd3vCtJGVA83TqZmZkEFAAAkoyV0zM4SRYAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQWA47S2tsowDBmGodbW1kRXB0ACEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjhB1QtmzZoksvvVQFBQUyDEPPP/+83+2maWrevHnKz89Xnz59VFJSor/97W9+ZQ4ePKgpU6YoMzNT/fr103XXXaeWlpaonggAAOg+wg4ora2t+ud//mctX7484O2LFi3SsmXLtGLFCm3btk3p6emaMGGCjhw54i0zZcoU7dy5U3/605+0ceNGbdmyRTfccEPkzwIAAHQrhmmaZsR3NgytX79ekyZNknSi9aSgoED/+Z//qZ/97GeSpKamJuXm5uqpp57SlVdeqd27d2vEiBF64403NGrUKEnSSy+9pIkTJ2rv3r0qKCgI+bjNzc3KyspSU1OTMjMzI60+AIfyfMYlqaqqShdffLFcLleCawUgWuH8ftt6Dkptba0aGhpUUlLi3ZaVlaUxY8Zo69atkqStW7eqX79+3nAiSSUlJUpJSdG2bdsC7vfo0aNqbm72WwB0T5WVlRoxYoR3feLEiRoyZIgqKysTWCsA8WZrQGloaJAk5ebm+m3Pzc313tbQ0KCcnBy/23v16qXs7GxvmY4WLFigrKws71JYWGhntQE4RGVlpcrLy1VfX++3vb6+XuXl5YQUoAdJilE8c+fOVVNTk3epq6tLdJUA2KytrU0zZsxQoF5nz7aZM2eqra0t3lUDkAC2BpS8vDxJUmNjo9/2xsZG7215eXk6cOCA3+3Hjx/XwYMHvWU6Sk1NVWZmpt8CoHupqanR3r17g95umqbq6upUU1MTx1oBSBRbA8rQoUOVl5enTZs2ebc1Nzdr27ZtKioqkiQVFRXp0KFD2rFjh7fMK6+8ovb2do0ZM8bO6gBIIvv377e1HIDk1ivcO7S0tOiDDz7wrtfW1uqdd95Rdna2Tj31VM2cOVP33nuvTj/9dA0dOlR33323CgoKvCN9zjzzTP3Lv/yLrr/+eq1YsUJfffWVbr31Vl155ZWWRvAA6J7y8/NtLQcguYU9zHjz5s0aP358p+3Tpk3TU089JdM0dc899+jRRx/VoUOHdOGFF+rhhx/WGWec4S178OBB3XrrrXrhhReUkpKiyZMna9myZcrIyLBUB4YZA91PW1ubhgwZovr6+oDnoRiGIbfbrdraWoYcA0kqnN/vqK6DkigEFKB78ozikeQXUgzDkCRVVFSorKwsIXUDEL2EXQcFAKJRVlamioqKTt29brebcAL0MAQUAHHV2toqwzBkGIZaW1s73V5WVqZdu3Z516uqqlRbW0s4AXoYAgoAx/E9x2TcuHGccwL0QAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAHHV1tbm/XvLli1+6wDgQUABEDeVlZUaMWKEd33ixIkaMmSIKisrE1grAE5EQAEQF5WVlSovL1d9fb3f9vr6epWXl/uFlPT0dJmmKdM0lZ6eHu+qAnAAAgqAmGtra9OMGTNkmman2zzbZs6cSXcPAC8CCoCYq6mp0d69e4Pebpqm6urqVFNTE8daAXAyAgqAmNu/f7+lcuPHj1dra2vA21pbW2UYhgzDCFoGQPdBQAEQc/n5+YmuAoAkQ0ABEHPFxcVyu90yDCNwAZ/ta9+s0xsfH4xTzQA4FQEFQMy5XC4tXbpUkjqHFMOQfE6enbdhpy5bsZWQAvRwBBQAcVFWVqaKigoVFBT4be+fk6f+l/zMuz7ujAGSpA8PtMS1fgCchYACIG7Kysq0a9cu73pVVZWWrX9NJ5/+Xe+2kjNzE1E1AA7TK9EVANCzuFwu79/jxo3TCzv/HrDcyzsb/NYLMoKcvwKgWyKgAIip1tZWZWRkSJJaWkJ325w28MSVY6v3fKbqPZ95t7cfOxKbCgJwJAIKgLja8ck3J7+ufbNOr37U7Hf7eYOztfamIr9zUF7e2aBN79bFrY4AEo+AAiBudnxyUFN/u927Pm/DTqX0TutU7vwh2Tp/SLbfNgIK0LMQUADEzUef+V8Bdn7pSKX1OVkFGYa+tzhBlQLgSAQUAAlz2ahCpaenc+l6AJ0wzBgAADgOAQVAt8UEg0DyIqAAAADHIaAAAADHIaAAAADHIaAAAADHIaAAAADHIaAAAADHIaAASDpr36zTGx8fDF0QQNKyPaC0tbXp7rvv1tChQ9WnTx8NGzZMv/rVr2SapreMaZqaN2+e8vPz1adPH5WUlOhvf/ub3VUB0E0My8nwW5+3YacuW7GVkAJ0Y7YHlPvvv1+PPPKIfvOb32j37t26//77tWjRIj300EPeMosWLdKyZcu0YsUKbdu2Tenp6ZowYYKOHGE6dQCdnT8kWyt/Mtq7Pu6MAZLkN+MxgO7F9oDy+uuvq7S0VD/60Y80ZMgQlZeX6+KLL9b27SdmMDVNU0uWLNFdd92l0tJSfec739Ezzzyjffv26fnnn7e7OgCSQFtbm/fvLVu2+K17nO3O8v7tPvKJzPbOZQB0H7YHlAsuuECbNm3S+++/L0n63//9X/3lL3/RD3/4Q0lSbW2tGhoaVFJS4r1PVlaWxowZo61btwbc59GjR9Xc3Oy3AEgOOz75phvmz7sbO91eWVmpESNGeNcnTpyoIUOGqLKyMmiZ+2dNU/2K67S9+sUY1RpAotk+m/GcOXPU3Nys4cOHy+Vyqa2tTb/+9a81ZcoUSVJDQ4MkKTc31+9+ubm53ts6WrBggX75y1/aXVUAMfbGxwc19bfbvetb3v/c7/YNGzZo6tSpfueoSVJ9fb3Ky8tVUVEhSSovL+9Upu3w51o852ZdePpAlZWVxegZAEgU21tQnnvuOa1atUqrV6/WW2+9paeffloPPPCAnn766Yj3OXfuXDU1NXmXuro6G2sMIFY6niMyv3Sk37kks2fP7hQ8JHm3zZgxQzNmzAhYxmPmzJkBu4QAJDfbW1Buv/12zZkzR1deeaUk6ayzztInn3yiBQsWaNq0acrLy5MkNTY2Kj8/33u/xsZGnX322QH3mZqaqtTUVLurCiDOLhtV6LdeX18ftKxpmtq7d2+IPZqqq6vTH//4R283MoDuwfYWlC+++EIpKf67dblcam9vlyQNHTpUeXl52rRpk/f25uZmbdu2TUVFRXZXB0APEKx7GEDysj2gXHrppfr1r3+tP/zhD/r444+1fv16Pfjgg/q3f/s3SZJhGJo5c6buvfde/f73v9e7776rf//3f1dBQYEmTZpkd3UAJJjvaJstW7YoLS1NpmmqurratsfwtMwC6D5s7+J56KGHdPfdd+uWW27RgQMHVFBQoBtvvFHz5s3zlpk9e7ZaW1t1ww036NChQ7rwwgv10ksvKS0tze7qAEig7dUvav/jd3rXJ06cKLfbraVLl6q0tFRut1v19fUBzzExDEODBg2SpKBlPMaOHWt/5QEklGF29al3qObmZmVlZampqUmZmZmJrg6AACorKzV5crkk/68YwzAkyW+EjiS/AGKljK+Wlhalp6d32t7a2qqMjIwuywCIn3B+v5mLB4Dt2traNGPGDHUMJ9I3IWPmzJkqLS1VRUWFCgoK/Mq43W5VVFSorKxMZWVlAcu4MvrHrP4AEo+AAsB2NTU1XY7AMc0To29qampUVlamXbt2eW+rqqpSbW2t37VNOpa5Y/HTyr/u4dhUHoAjEFAA2G7//v1hlXO5XN5t48aN81v38N02/JzRMlI6lwHQfRBQANjO9xpHdpSLlJU5fgA4EwEFgO2Ki4vldrslGQFvNwxDhYWFKi4ujlkdrMzxA8C5CCgAbOdyubR06dKAt3lG6CxZsiRgV44dKisrVV5e3ulKtZ45fggpgPMRUADERFlZmWYtfKTTaBvfETqx4BlB1NUcP8zfAzgfAQVAzIwe/0O/0TaBRujYLZwRRACci4ACIKZ8R9sEG6Fjp3BHEAFwJgIKgG7FKSOIAETH9rl4ACDe1r5Zp7Q+J2tYToZ3BFFXc/y43e6YjiACED0CCoCkc9pA/zl15m3YqZTeJyYbXXtTkZYuXary8nIZhhFwjp9YjiACYA+6eAAknfMGZ2vlT0Z71+eXjtT4bw2UJH14oCXo/D2xHkEEwD4EFABJ6bzB2d6/LxtVqAkj8/xutzLHDwDnIqAAsKS1tVWGYcgwDLW2tkZcJpbajx3RVWMGex/fyhw/AJyJgAIAAByHgALANjs+Oej9e+2bdXp5Z0PMHqvjRIDtXBkW6FYYxQPAFm98fFBTf7vdu+47ssZuGzZs0OzZs73rEydOVHZOvozRU2LyeADij4ACwBYfHmjxW59fOlJpfU5WQYah7y2297GmTp3a6RonBw80SBsfsPeBACQMAQVATFw2qlDp6ekxOVk20AXYpEDbACQrzkEBAACOQ0ABAACOQ0ABAACOQ0ABkHDp6ekyTVOmaSo9PT1omePHj8vtdnvn1AHQfRFQACQNl8ulpUuXSlKAkEJoAboTAgqApBJsIsD+uXnqf8nPElQrAHYjoABIOoEmAly2/jWdfPp3E1grAHYioABISh0nAkxhIkCgWyGgAAAAx+FKsgBiyjNCBwDCQQsKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAEva2tq8f2/ZssVv3cNsD10mlhL9+ADsQ0ABEFJlZaVGjBjhXZ84caKGDBmiyspK77bt1S9q/+O3dFkmlgI9/siRI7Vu3bou5/gB4EyGmYTj/5qbm5WVlaWmpiZlZmYmujpAt1ZZWany8vJOQ4U9c+FUVFRIkiZPLpcUvExZWZmt9WptbVVGRoYkadWqVZoyZWpcHx9A+ML5/SagAAiqra1NQ4YM0d69ewPebhiGBg0aJEldlnG73aqtrfW7+mu0fAPKoEGDVF9fH9fHBxC+cH6/6eIBEFRNTU3Q4CFJpmlq7969IcvU1dWppqYmFlWUpKDhJF6PD8B+BBQAQe3fv9+R+0rGxwcQHgIKgKDy8/Mdua9kfHwA4SGgAAiquLhYbrfbe7JpR57zO9xut6TgZQoLC1VcXByzeg4aNChoHaXYPz4A+xFQAATlcrm0dOlSSeoUADzrS5cu9ZbpyFNmyZIlMT1BddGiRQHr6BHrxwdgPwIKgC6VlZWpoqJCBQUFftvdbrd3+G5ZWZlmLXxEroz+QcvEUmlpacA6uvoO0KyFjzDEGEhCDDMGYInncydJVVVVuvjii/VWXZM+PNAiSXp5Z4P+/L8fa+/SK/zKxKrlwneYcUtLi9LT0/3qeMfip/Xsvn66v/xsXTn61JjUAUB4wvn97hWnOgFIcr5BY9y4cXqrrkmXrdjqV8ZI8S8T724V38cbfs5oGQ0fxvXxAdiHgAIgIp6Wk/HfGqgJI/MkSQUZhr63OJG1AtBdEFAARGXCyDxvF0pra2vcHjc9Pb3T5fcBdB+cJAsAAByHgAIAABwnJgGlvr5eU6dOVf/+/dWnTx+dddZZevPNN723m6apefPmKT8/X3369FFJSYn+9re/xaIqAGKo/dgRXTVmsAzDiGv3DoDuz/aA8o9//ENjx47VSSedpBdffFG7du3Sf/3Xf+mUU07xllm0aJGWLVumFStWaNu2bUpPT9eECRN05MgRu6sDAACSkO0nyd5///0qLCzUk08+6d02dOhQ79+maWrJkiW66667VFpaKkl65plnlJubq+eff15XXnml3VUCAABJxvYWlN///vcaNWqULrvsMuXk5Oicc87RY4895r29trZWDQ0NKikp8W7LysrSmDFjtHXr1kC71NGjR9Xc3Oy3AACA7sv2gPLRRx/pkUce0emnn66XX35ZN998s/7jP/5DTz/9tCSpoaFBkpSbm+t3v9zcXO9tHS1YsEBZWVnepbCw0O5qAwAAB7E9oLS3t+vcc8/Vfffdp3POOUc33HCDrr/+eq1YsSLifc6dO1dNTU3epa6uzsYaAwAAp7E9oOTn52vEiBF+284880x9+umnkqS8vBNXnGxsbPQr09jY6L2to9TUVGVmZvotAACg+7I9oIwdO1Z79uzx2/b+++9r8ODBkk6cMJuXl6dNmzZ5b29ubta2bdtUVFRkd3UAAEASsn0Uz6xZs3TBBRfovvvu0+WXX67t27fr0Ucf1aOPPipJMgxDM2fO1L333qvTTz9dQ4cO1d13362CggJNmjTJ7uoAAIAkZHtAOf/887V+/XrNnTtX8+fP19ChQ7VkyRJNmTLFW2b27NlqbW3VDTfcoEOHDunCCy/USy+9pLS0NLurAwAAklBMJgu85JJLdMkllwS93TAMzZ8/X/Pnz4/FwwMAgCTHXDwAAMBxCCgAAMBxCCgAAMBxCCgALElPT5dpmjJNU+np6RGXiSXfx0/rc3LcHx+AfQgoAADAcQgoAADAcQgoAADAcWJyHRQAPYPZ3ub9e8uWLbr44ovlcrkSWKPOXt7pP0v6sJwMnT8kO0G1AWAVAQVARLZXv6j9j9/pXZ84caLcbreWLl2qsrKyBNbshGE5GZKk6j2fqXrPZ363rb2piJACOJxhmqaZ6EqEq7m5WVlZWWpqamJmYyABKisrNXlyuST/rw/DMCRJFRUVjggpb3x8UB8eaPGuv7yzQdV7PtPCsrN05ehTE1gzoGcK5/ebFhQAYWlra9OMGTPUMZxIkmma3glBS0tLE97dc/6Q7E4tJR1bUwA4EyfJAghLTU2N9u7dG/R20zRVV1enmpqaONYKQHdDQAEQlv3799taDgACIaAACEt+fr6t5QAgEM5BAWCJ54TT9tTBys7J18EDgVtIDMOQ2+1WcXFxnGsIoDshoAAI6Y2PD+qyFVu9664LrpWev69TOc8oniVLliT8BFkAyY2AAiAkz1Dd8d8aqAkj8ySdpe3fPVUrF/9SnzV+05Lidru1ZMkSRwwxBpDcCCgALJswMs97/ZArR9+oX9x8lbKysiRJVVVVjrySLIDkxEmyACLmG0bGjRtHOAFgGwIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKgIBaW1tlGIYMw9CRL79IdHUA9DAEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAA9SvuxI7pqzGAZhqHW1tZEVwdAEAQUAADgOAQUAADgOL0SXQEAySs9PV2maSa6GgC6IVpQAACA4xBQAACA4xBQAPRYa9+s0xsfH0x0NQAEQEAB0GMMy8nwW5+3YacuW7GVkAI4EAEFQI9x/pBsrfzJaO/6uDMGSJI+PNCSqCoBCIKAAiCgtrY279/vvb1dZntbF6WTx3mDs71/l5yZm8CaAOgKAQVAJ5WVlRoxYoR3/f5Z01S/4jptr34xgbUC0JMQUAD4qaysVHl5uerr6/22tx3+XIvn3KzKysoE1QxAT0JAAeDV1tamGTNmdHnxtZkzZ/p1/wBALBBQAHjV1NRo7969XZQwVVdXp5qamrjVCUDPREAB4LV//35bywFApAgoALzy8/NtLedE3XV0EtDdEFAAeBUXF8vtdksyAt5uGIYKCwtVXFwc34rZhNFJQPIgoABQa2urDMNQr169dP/99wcsYxgnQsuSJUvkcrniWT1bMDoJSC4EFAB+SktLNWvhI3Jl9Pfb7na7VVFRobKysgTVLHKMTgKST8wDysKFC2UYhmbOnOndduTIEU2fPl39+/dXRkaGJk+erMbGxlhXBYBFo8f/UPnXPexdr6qqUm1tbVKGE4nRSUAyimlAeeONN/Tf//3f+s53vuO3fdasWXrhhRe0du1avfrqq9q3b1/SfvEB3ZWR8k03zrhx45KyW8eD0UlA8olZQGlpadGUKVP02GOP6ZRTTvFub2pq0uOPP64HH3xQ3//+93XeeefpySef1Ouvv66//vWvsaoOgB6sJ4xOArqbmAWU6dOn60c/+pFKSkr8tu/YsUNfffWV3/bhw4fr1FNP1datWwPu6+jRo2pubvZbAMAqz+gkz4m+nSX36CSgO4pJQFmzZo3eeustLViwoNNtDQ0N6t27t/r16+e3PTc3Vw0NDQH3t2DBAmVlZXmXwsLCWFQbQDflcrm0dOlSSQoaUpJ1dBLQXdkeUOrq6jRjxgytWrVKaWlptuxz7ty5ampq8i51dXW27BdAz1FWVqaKigoVFBT4bXf1HaBZCx/hPDjAYWwPKDt27NCBAwd07rnnqlevXurVq5deffVVLVu2TL169VJubq6OHTumQ4cO+d2vsbFReXl5AfeZmpqqzMxMvwUAwlVWVqZdu3Z51+9Y/LQG3fS4Ro//YQJrBSCQXnbv8Ac/+IHeffddv23XXnuthg8frjvuuEOFhYU66aSTtGnTJk2ePFmStGfPHn366acqKiqyuzoA4Me3G2f4OaNlNHyYwNoACMb2gNK3b199+9vf9tuWnp6u/v37e7dfd911uu2225Sdna3MzEz99Kc/VVFRkb773e/aXR0AAJCEbA8oVixevFgpKSmaPHmyjh49qgkTJujhhx8OfUcAANAjxCWgbN682W89LS1Ny5cv1/Lly+Px8AAAIMkwFw8AvzlotmzZonbmpAGQYAQUoIerrKzUiBEjvOsTJ07UTyeN1Rd/48rOABInIeegAHCGyspKlZeXd5rl9+CBBmnjAwmqFQDQggL0WG1tbZoxY0ancHJCoG0AED8EFKCHqqmp0d69e0OWq66uVnp6ehxqBADfoIsH6KH2799va7lk9vJO/3nAhuVk6Pwh2QmqDQCJgAL0WPn5+baWS0anDTzRMlS95zNV7/nM77a1NxURUoAEIqAAPVRxcbHcbrfq6+uDnIci9c/NV3FxcZxrFj/nDc7W2puK9OGBFu+2l3c2qHrPZ/rwQAsBBUggAgrQQ7lcLi1dulTl5eUyDKNDSDEkmfr3Wff4zV3THZ0/JLtTEOnYmgIg/jhJFujBysrKVFFRoQE5/jOJn5ydo4GTfs4svwAShhYUoIcrPPcipV65RFp6hSRpYPkv1GfoOTJSXBqWk5HYysVAenp60C4tAM5BQAF6uA8PtMhI+aYbZ9GtVyitz8mMZAGQUAQUAH4uG1XIdU8AJBznoAAAAMchoADo0VpbW2UYhgzDUGtra6KrA+BrBBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQASumdpme3fSLTNJmHB4AjEFAA9GhtbW3ev7ds2aJ2n3UAiUNAAdBjVVZWasSIEd71iRMn6qeTxuqLPa8nsFYAJAIK0CP5TpB35MsvEl2dhKisrFR5ebnq6+v9th880KDPnr9P26tfTFDNAEgEFAA9UFtbm2bMmCHTNAPcemLbM4t/6df9AyC+CCgAepyamhrt3bu3yzJ/b9yvmpqaONUIQEe9El0BAIi3/fv3Wyq3ruZdNZx8moblZOj8IdkxrhUAXwQUAD1Ofn6+pXJrd7fohdZ3T/x9UxEhBYgjungA9DjFxcVyu90yDCNwAcNQ/9x8PfjTKzX+WwMlSR8eaIljDQEQUIAeaMcnB71//3l3YwJrkhgul0tLly6VpE4hxTAMGZIeffg3urpoqCaMzEtADQEQUIAe5o2PD2rqb7d717e8/7kkaVhORqKqlBBlZWWqqKhQQUGB33a3262KigqVlZV5t7UfO6KrxgyWYRhqbW2Nd1WBHolzUIAepmNXxfzSkRo5OKdHnl9RVlamkpISZWVlSZKqqqp08cUXy+VyJbhmAAgoQA932ajCHj3/jm8YGTduHOEEcAi6eIAeyGz3n3+GC5IBcBoCCtDN+V7WvrW1VdurX9T+x2/x3j5x4kQNGTJElZWVCawlAPgjoAA9yIYNG7R4zs1qa/m73/b6+nqVl5cTUgA4BgEF6EFmz54tz1wzvjxz0sycOZPuHgCOQEABepCOM/f6Mk1TdXV1zD8jZnsGnICAAsCP1XlqACCWCCgA/FidpwYAYomAAvQggwYNkhR4/hnDMFRYWKji4uL4VioJMCwbiD8CCtCDLFq0KOB2z3w0S5Ys4UJlHby55Y8MywYSgIAC9CClpaWatfARuTL6+20PNP9MT5Geni7TNGWaZsAr6i6fN5Nh2UACEFCAbs63O2LLli0aNe5i5V/3sHdbVVWVamtre2Q4sYZh2UAiEFCAbsZ3iOzq1as1YsQI720TJ07UTyeN1ZcfvuHdxvwzkfEMy+7VqxczHAMxwGSBQDc2depU73/7HgcPNEgbH0hQjQDAGlpQgG6sYzj5emvc65Fs6LYBEs/2gLJgwQKdf/756tu3r3JycjRp0iTt2bPHr8yRI0c0ffp09e/fXxkZGZo8ebIaGxvtrkpcdZyQDUByqqys9OsWA5AYtgeUV199VdOnT9df//pX/elPf9JXX32liy++2O9He9asWXrhhRe0du1avfrqq9q3bx8n6AFIuMrKSpWXl3c5JQCA+DDMwG3Atvnss8+Uk5OjV199VePGjVNTU5MGDhyo1atXq7y8XJL03nvv6cwzz9TWrVv13e9+N+Q+m5ublZWVpaamJmVmZsay+pa1trYqIyNDktTS0hJwuCIQD77vRSt4v57Q1tamIUOGaO/evSHL9s/Nl86for9/fS4PxxCwJpzf75ifg9LU1CRJys7OliTt2LFDX331lUpKSrxlhg8frlNPPVVbt24NuI+jR4+qubnZbwFgReCrxsoIsr0Hq6mpsRROFi5cqGXrX9PJp4f+ZwpA5GIaUNrb2zVz5kyNHTtW3/72tyVJDQ0N6t27t/r16+dXNjc3Vw0NDQH3s2DBAmVlZXmXwsLCWFYb6NYMwwgWW3o0q5Mk5uTkKIVh2UDMxTSgTJ8+Xf/3f/+nNWvWRLWfuXPnqqmpybvU1dXZVEOg+/EdgTLp2lvlysj2u93tdmvlypXxrpbjWZ0kMS8vTxLz8wCxFrOAcuutt2rjxo2qrq6W2+32bs/Ly9OxY8d06NAhv/KNjY3eD35HqampyszM9FsAdNZxBMrzTz4k37PMPFeNLS0tTUDtnK24uFhut9s7L1EwY8eO1fbqF5mfB4gx2wOKaZq69dZbtX79er3yyisaOnSo3+3nnXeeTjrpJG3atMm7bc+ePfr0009VVFRkd3WApBXu0PVgI1DaWw96//ZcNTbU/DM9kcvl0tKlSyWpU0jxXd+4caMWz7k54Pw8kydP5lIDgE1sDyjTp0/XypUrtXr1avXt21cNDQ1qaGjQl19+KUnKysrSddddp9tuu03V1dXasWOHrr32WhUVFVkaweNUHec7obkX8dTW1qYZM2YEuTCbfzkEV1ZWpoqKChUUFPhtHzRokPfv2bNnq6v5eSSOM2AH2wPKI488oqamJl100UXKz8/3Lr/73e+8ZRYvXqxLLrlEkydP1rhx45SXl5fUTaMdm9Vp7kW8WR2B8tprr8WhNsmtrKxMu3bt8q5XVVVp586d3nUr10jhOAPRs30uHiuXVUlLS9Py5cu1fPlyux8+7jzN6h2ft2c69p46hT3iy+oIlGAj5eDPd/LEcePGhX1/jjMQPebiiUJXzepMx454CncECmKL4wxEj4AShVDN6p7p2GtqauJYK/RE4YxAQfg8JxUfP37861GJHGcg1ggoUbDarG61HBCprkagdCyHyPke5458jzvHGYgeASUAq8M7rTare8rZNeMxMyf3DOGODAs2AsWV0T8m9eupysrKNGvhI52Oq+9xZyQfED0CShRCNasbhqHCwkIVFxfHuWZIdpGODCs89yLNf+Zl73rJrMUadPMTenbbJ1zzxEajx/9Q+dc97F2fdO1Pdez4N4GEkXxA9AgoUbByYaclS5bQ3IuwBLvgmmdkWLAfvTc+PqjLVmzVLza+5922R4UyUlwalmN9dmOENiwnQ0bKN5/r5598SJ81+o/cCfV6AegaASVKwZrV3W43Q4wRtmhGhn14oEWSNO6MAd5t80tHau1NRTp/SHan8ojc+UOytfIno7ssw0g+IDoEFBsEurBTbW2tpXDC+STwZXVkWK9evYK+X0rOzPX+fdmoQsJJBKxMBXDe4NDHlZF8QOQIKDbpeGEnunUQCUZ8dU+8rkD4CCgBxGpeHbv2y7w/3ZfVkWFILryuQPgIKB3YOa+Ob/fN6tWrA+53w4YNCaufVQyPjh8rI8PgRIzkA+xGQPER6egJK6ZOnRpwv1OnTnVE/eAMVkaGITkwkg+IDgHla7GeVyfYfn23d9Vdw7w/PUdXI8OeeeYZ73qw90tan5NDnuCJ6PmeSBvowm2M5AOiQ0D5mhPm1emqu8YJ9UP8BBoZ9uCDD2rOnDnebVwMzDk6XrgtnJF8AAIjoHzNKfPqBOquaW1t1fjx4y3dv6v6JeM5IMlYZ7v4dgv84x//0OWXX96pe29vfb0mTy7Xbff/t17e2dBxF4gj3wu3dfeRfD35c4n4IaB8Ldx5dTqyct0EK6LtrmG0QPc0e/bsgN17Mk1Jppb9+i69svtEQOGqsQC6AwLK12Ixr04056sE665JxLw/DI9OvI4tJx21Hf5c04Z+yVVjHWLtm3Vas/1TvfHxwURXBUhaBJSv2T2vTsfhwJEI1l1jpX52NcHaNaw5EcOjEyXQsY9Hk/jQtC8IJwnSsdVq3oadmlP5ri5bsZWQAkSIgOLDrnl1gg0HDleg7pqVK1fGbd4fu4Y1Mzw6PvLy8hJdhR6r49w880tHavy3Bkr6Zo4kAOEhoHQQzbw6UtfDgX0NGDDAcneNb1fIKaeconfffTei+oXTxWLXsOZo90O30AmDBg0KeR2UsWPHxqk2CMR3bp7LRhVqwsjuGxj5XCIeCCgBRDOvTqjhwB433nijpNDdNYG6Rs4666yw6xduF4tdw5qj2U9P6hYKxPfE62XLlknq+mJt3XnUSDKw60R5p+vpn0vEDwHFZlaHIQ8bNixkd1KwrpF9+/Z12l9X5zhs2LChyy6WQPeLdti1pz6RDo8Op1sonud8JGp4ZVlZmRY+/KT6Dcj1297x4mBwlvZjR3TVmMFdvl+cNmS3q/o4pbvWaccMsUFAsZnVYb55eXlddidZ6RqRgo8U8t1uZT+RPg+7hjX77ifZr5obqPk72ibxNz4+qIc/HqD0q5d6tw0s/4Xyrv1NVPsFrEr2zyWSDwElgJ37mrx/e4YLhlo8Z+qHGq7s4TlfIFh3ktWuotdee63Ttg0bNvg1wX7++edB7x8spNg57Lqr8ycC7SeZr5obqPk7NzdXQ4cO9dsWbpO450TL7w3/pgXl8u+comPP3RbVfgGrkvlzieTUK9EVcKLyR7Z6/563YadSeqdZut9PLhyqf8rJ0GW33qXFc27WiRlOfQPAN+uVb+9TWp9DOvLlF95b175Zp5GDc3T+kGzLXSwNDZ2vHjp16tSQJ+mG4hl27ekC8t1fuMOuFy1apKlTp1rej1Ou6hsuT/N3x2P/97//vVNZT5N4uKOvSs7M1aqv/374npmdHivS/SL2fD/fyShZP5dIXrSgBFBxc5H37/mlI7Ww7Kwul59ceOK/49/+pVZzKt9V5T/cGjhprlwZ/l9Evuue6yTM27DTb9tlK7bq3o27tPuQtZdmT5NLa7Z/qu0ffebdFmk46dhFEM2w644jj5577jnL+7HabXT11VertbXVti4VK/3awfZrdfSWh+dkysmTJ6u5uTlgmTc+PuhtoQt0GXua2p1tWE6GzPZvXoPZv/mdypb8udN7rON7qrm5OWbnVEXzHg+32zeW54nE6vMNhzGTUFNTkynJbGpqSnRVvLbX/t18dtsnfsvjm/7vxHXIJfOOxU+bq17/qFOZJzfv9pa5e+0b5uA7NpqD79honnr7BtPVd4D3tkCLq+8A89TbN5gDJ/3cdGX077Ks1SUnr8Bct26d33PzHG9JZlVVlXn8+PEuj8W6devMQYMG+e3X7Xab//M//2NpP8ePHzfdbrdpGEbAOvpuX7VqVafH6t+/v5mdnd3p8Ts+r45aWlq85VtaWiw/r3Xr1pnV1dURH/OqqqqA7yfPe8F3eXXnp5b3W11d3eXzRWytW7fOHJib7/+Z9fmctrS0BHxP+a573oeh3ptWRfMet/K5LCws9H6u7apzOHWM5rkjPsL5/SagxJCVD0THMr5BZ9bCFaZkfL34fhmc2DZr4QqfMtGHE9/9z1q4ImCIenLz7k4ha3vt373PZ926dQG/wAzD8Nse6gvCs5+O++q4LdiXZbDH7+pLrKvXK9TzmjlzZsTH+4knnuhUl2e3fWIOvmOj+eMntvkdZ986hlpWr17d5TFG7AR7v/guq1atCvqeSkRACfUeX7duXcjPpe/nKxaBwEodI3nuiB8Cis0CvbEjCR+RCPTfQmFhod9/NPaGE5//Smb8zhx8x0azcFbFN489qyLgf/Z3r33DUvDxDTodg0/HH+BArSOFhYV+LTGRLMHeN8FaikIdZ8MwzIEDB0Zcn4ULF5pbPzjgF/p+/MQ2c/AdG81nt33SqZ5WW2toQUkMq5/LAQO6biGVZG7cuNE8fvx42K2YwQTaj5XQ69s60tV3UqjH6iic70grn0PfFhwrzz1WIv3N6AkIKDZLZEAxzeAfrGi6Fawsdyx+OmQLyq9e2NkpxFhZCmdVBAw+vkHnyc27O3WTzbzvYfOUgXlRPa9AXSp2dN8MHDjQcotOx8XVd4A5cNLPOwU/39Ypj3Cb2hFfdn8uI+2y7CjYe3zVqlWW6+IJvaF+7K12w4TzHRlNMI+0WyhSBJTgCCg2S3RACbaf1atXW/5iCdQsG2rxdD2Eeh7ba//uF2KsLB1bUDzn33QMLb7r/S/5mS1f+Dfe/YBfyAraTRZm983MmTMjOs6+i2/XWqBw4hFOUzviK5zPZaRLuK+z1a7XUIun2zDariKPcL4jrR7Xjl2b0XQLRYqAEhwBJUaPJ33z30KimwtNM7z/1AoLC8P6T0k60fVg9Xn5Ho9QS7D/wrZ+cKBTi4nverQtJ54la9w089TbN1g8GdkwM0/JtrTfux9eY85auKJTPVPS+popaRmh9xFmy4fVpnbEV6xbNn1/YK28X+zsCvZ8doN9J4XbDRPO92gkLSjRdgtFKtG/GU5GQLFRoB8Bu5pcAwkneYdq6vcsnn5s3/0MGjTI0n9OVp5XoGMU6svA6nG1st9IluycfHPWwhXm3Q+vsWV/nhFVg+/YaLpn/M67fWD5L8xTb99g3rnmr2Htz+p/WHzpOY/Vz6VdS6jQYEdgsjpCJ5zHCnSOWVffN1aPq+/vQrzO1wp17pydvxnJ3hJDQLGJlTPxfT/AdjQXhvvm66qpv6v9BLtfuM/L6jEKNBIgnOMazSgZ+/bb9YiqYKOewh19E84XT7J/WXVXVj9fdiyhul2i7XIK57st3G7nSL9vrH7fRdotFC7fYx/r34xk/8wTUGwQSbOoHc2FkfxHHOwEsFBv4khaPiI9RtGMPAp3lEyg/1bs2O+dd95pFhQUBHxevkK1gllZrL72yf5l1Z1Z/XxFu4Q6cTXaFpRwug1j0Vpj5bj6rkfy3KNtQQmnizuc5xrqsZKx1ZSAYoNoPmiRvtmjOdO845vWd72rH65Izh0J9xj5nssSzXG1Mkom0LBMO/brWXwDSjhDJyP5D4uLTyW/SH+4wv1xi+YCa1Y/u1bY2b3V1feo73EN9I9DpBeXi4QdQdTqb0a8RyPFAgHFBtE0i0bSXBjtmeaR/lBFcuGvcFsDfC9EFs1xDTZKJtSFrSLdb7AvtK4eK9IL0AV7LEbkJLeOwTTUhQdl9YfdQpdpoIsaxmokny87u7daWloiCvyRXlwu3Nc02MX2wl2s/GbYPRopUf/chPP7zVw8QVidd8KO+ybLNOae5xVuPfLy8jrtIxKlpaUB5wYaNGhQxPvsar+B+L5GvvPweISaFyTY/EZdPZYTXntEb+XKlQHno1q5cqV3fdXKlRqY6/8ZOblvP6WkZfhtc2X01xVzFutI/jm6/uZbg353+G7/+c9/Hvbnz/PZjdV7PJQtW7bo2LFjXT52qO/NYJ9vK3OKBeL7+OHMv9WVUK9LLH4jIpnPKO5iFJJiKp7noISTjiNtLrSjn9SOFhQrzaCRNGf6vk52HFcr3VnhPq+O+7WyVFVVRdzkGu5jcVXY5NTxcxnqSq7Byrzy/2q92675xQrviLHcq+4L633Ub0Cu9+++WaeYoabJaGpqitt7PNiSkpLS6bFjcXE5K+w+r8jqb4bd59IksquILh6bxKs5Pl5nmnfFSjOoXaOaoj2udjQ3h9qvlcXTNRTJeyHcx2Jene4h0gs8dtzmmbPr1vnLIv5x7PrCh6Hm+uo8eq2riVDtXMIdgWfl4nJWhPP9F+13ZEd2/kYk4sJ1vggoNrJ6vY5oLpDllLlV7Jz3J9TxiOa4RjIUO9ScPh1bYqwsXY0ACvWfUbgBhRaU7sGugOIRzUnnj2/6v4AXFvRMuxDOjOqBrv/ju27/Yv0CilLniz4GmvT02W2f+LVUhTMfV6gl1DWegl3M7eDBg7Z+T1h9Hp7BBrFAQInR43X1RormxXTS3Cp2DFW0ejziebXFjvtdu3ZtVNPKWx2ebPULwwmvPZJPNKNmurqi87PbPrF8EcNgV1D2XAxx1esfmdk5+ab9s67LTOmTablssPmuPMvAST83XRn9Az6HaC7oePt/PWGuev0jv4A06dqfmqcMzPUrl5HZz0zP7Of//Dp0bwVb7O4q8v0+tBsBxWbxmFfBKXOr2HGxp2guMharM8utnHkfzpn/4TYvd8Uprz2SU6SjZuy6wJuVrk4rF1iLZPnhldeZgS+gGGgJ3i0VvCvrxNJ3VGnEdfRMhBruhKrWl9DdbZ4lvC7B2Hz/EFBsFq8fUifMrWJHU7LTA0pXJ7mFOiHY83rE46Q15tVBMMePHzerq6vN1atXm9XV1QFbBK2+N6P9zFvt6gx1QclIlurq6rBOXO3Y0nD8+HHzz3/+c4gLO4bXnSTJzM75pjXJ050Uq3NyUjL6m1ljp5gDLr3dzL3qPm+Xm+9y6u0bzNyr7jMzi64Ie/92t+ASUJKYU68S2F26IsINFsFej1gcD6e+9nCWYD/0vudUDRgwIOr3Ziy6OiP9PIV6DpGMigt3RI7VCzp29dm1eyLJQF1Fnm4p39ahaCdatfMcOAJKEnPylUG7Q1dEuGfDd/V6xPICUE577eEMVi/K5unGjPa9aXdXZySfp2DPNRYj8KzcJ9ixsPLZjXZepEDHoqvjY9fIIztHERJQEDPJ3hVB1wySVagRGB1bFex6bya6q9Plcll+DnZ0S3X1HKI9pna2oIR6DgMGDLA0J5mdr58VBBTEVDJ3RdA1g2QVSRiw672ZyK5O36G2oZ6Dnd1SwZ5DNMfUjrmKInkO0TwW56CEiYCCaHSHrir0PE64oGMgTvs82dUtFavnEM1cRZE8h1BLqK6rRI7iYS4e9DjB5gqJdG4OIB6szqMTzXxXkXDa5ylUfUpLSy3vy+1263e/+52ys7P17LPPavPmzVHPWROsfv3791d2drbfNpfL1ak+4T6HUOyeq8hWtkajMP3mN78xBw8ebKamppqjR482t23bZul+tKDADnTNIJk4fSSd0z5P0Y4Y2rhxo7l27dpO5/3YdQEzKxeqDNa9ZVdXkV1dV5E8b0d38axZs8bs3bu3+cQTT5g7d+40r7/+erNfv35mY2NjyPsSUAD0RE7rTklW0cw95pRjbUdXUSKeQ1IElNGjR5vTp0/3rre1tZkFBQXmggULQt6XgAKgp2LkmD2imXss0a1VXT0HK0si3y/h/H4bpmmaFnqCbHXs2DGdfPLJqqio0KRJk7zbp02bpkOHDmnDhg1+5Y8ePaqjR49615ubm1VYWKimpiZlZmbGq9oA4AhtbW2qqanR/v37lZ+fr+Li4k7nKyC0YMdx8+bNGj9+fMj7V1dX66KLLop9RbvgeQ719fWaNWuWPv/8cwX7Wc/OztZzzz2niy66KGHvl+bmZmVlZVn6/e4Vpzr5+fzzz9XW1qbc3Fy/7bm5uXrvvfc6lV+wYIF++ctfxqt6AOBoLpcr4T+M3UGw47h//35L97daLpZ8n0OfPn1UXl4uwzD8QophGJKkxx57TD/4wQ8SUc2IJMUonrlz56qpqcm71NXVJbpKAIBuyqkjpkLxjBAaNGiQ33ZHjMiJQEJaUAYMGCCXy6XGxka/7Y2NjcrLy+tUPjU1VampqfGqHgCgBysuLpbb7VZ9fX3A7hLDMOR2u1VcXJyA2nWtrKxMpaWl3aILMCEtKL1799Z5552nTZs2ebe1t7dr06ZNKioqSkSVAACQdKLbZOnSpZK+6R7x8KwvWbLEsT/6nm6fq666KqHnm0QrYV08t912mx577DE9/fTT2r17t26++Wa1trbq2muvTVSVAACQ1P26S5JRQrp4JOmKK67QZ599pnnz5qmhoUFnn322XnrppU4nzgIAkAjdqbskGSVkmHG0whmmBAAAnCGc3++kGMUDAAB6FgIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwnIRd6j4anovfNjc3J7gmAADAKs/vtpWL2CdlQDl8+LAkqbCwMME1AQAA4Tp8+LCysrK6LJOUc/G0t7dr37596tu3b6epsKPV3NyswsJC1dXVMc9PDHGc44PjHB8c5/jgOMdPrI61aZo6fPiwCgoKlJLS9VkmSdmCkpKSIrfbHdPHyMzM5AMQBxzn+OA4xwfHOT44zvETi2MdquXEg5NkAQCA4xBQAACA4xBQOkhNTdU999yj1NTURFelW+M4xwfHOT44zvHBcY4fJxzrpDxJFgAAdG+0oAAAAMchoAAAAMchoAAAAMchoAAAAMfpkQFl+fLlGjJkiNLS0jRmzBht3769y/Jr167V8OHDlZaWprPOOktVVVVxqmlyC+c4P/bYYyouLtYpp5yiU045RSUlJSFfF5wQ7vvZY82aNTIMQ5MmTYptBbuJcI/zoUOHNH36dOXn5ys1NVVnnHEG3x0WhHuclyxZom9961vq06ePCgsLNWvWLB05ciROtU1OW7Zs0aWXXqqCggIZhqHnn38+5H02b96sc889V6mpqfqnf/onPfXUUzGvp8weZs2aNWbv3r3NJ554wty5c6d5/fXXm/369TMbGxsDln/ttddMl8tlLlq0yNy1a5d51113mSeddJL57rvvxrnmySXc43z11Veby5cvN99++21z9+7d5o9//GMzKyvL3Lt3b5xrnlzCPc4etbW15qBBg8zi4mKztLQ0PpVNYuEe56NHj5qjRo0yJ06caP7lL38xa2trzc2bN5vvvPNOnGueXMI9zqtWrTJTU1PNVatWmbW1tebLL79s5ufnm7NmzYpzzZNLVVWVeeedd5qVlZWmJHP9+vVdlv/oo4/Mk08+2bztttvMXbt2mQ899JDpcrnMl156Kab17HEBZfTo0eb06dO9621tbWZBQYG5YMGCgOUvv/xy80c/+pHftjFjxpg33nhjTOuZ7MI9zh0dP37c7Nu3r/n000/HqordQiTH+fjx4+YFF1xg/va3vzWnTZtGQLEg3OP8yCOPmKeddpp57NixeFWxWwj3OE+fPt38/ve/77fttttuM8eOHRvTenYnVgLK7NmzzZEjR/ptu+KKK8wJEybEsGam2aO6eI4dO6YdO3aopKTEuy0lJUUlJSXaunVrwPts3brVr7wkTZgwIWh5RHacO/riiy/01VdfKTs7O1bVTHqRHuf58+crJydH1113XTyqmfQiOc6///3vVVRUpOnTpys3N1ff/va3dd9996mtrS1e1U46kRznCy64QDt27PB2A3300UeqqqrSxIkT41LnniJRv4NJOVlgpD7//HO1tbUpNzfXb3tubq7ee++9gPdpaGgIWL6hoSFm9Ux2kRznju644w4VFBR0+lDgG5Ec57/85S96/PHH9c4778Shht1DJMf5o48+0iuvvKIpU6aoqqpKH3zwgW655RZ99dVXuueee+JR7aQTyXG++uqr9fnnn+vCCy+UaZo6fvy4brrpJv385z+PR5V7jGC/g83Nzfryyy/Vp0+fmDxuj2pBQXJYuHCh1qxZo/Xr1ystLS3R1ek2Dh8+rGuuuUaPPfaYBgwYkOjqdGvt7e3KycnRo48+qvPOO09XXHGF7rzzTq1YsSLRVetWNm/erPvuu08PP/yw3nrrLVVWVuoPf/iDfvWrXyW6arBBj2pBGTBggFwulxobG/22NzY2Ki8vL+B98vLywiqPyI6zxwMPPKCFCxfqz3/+s77zne/EsppJL9zj/OGHH+rjjz/WpZde6t3W3t4uSerVq5f27NmjYcOGxbbSSSiS93N+fr5OOukkuVwu77YzzzxTDQ0NOnbsmHr37h3TOiejSI7z3XffrWuuuUY/+clPJElnnXWWWltbdcMNN+jOO+9USgr/g9sh2O9gZmZmzFpPpB7WgtK7d2+dd9552rRpk3dbe3u7Nm3apKKiooD3KSoq8isvSX/605+Clkdkx1mSFi1apF/96ld66aWXNGrUqHhUNamFe5yHDx+ud999V++88453+dd//VeNHz9e77zzjgoLC+NZ/aQRyft57Nix+uCDD7wBUJLef/995efnE06CiOQ4f/HFF51CiCcUmkwzZ5uE/Q7G9BRcB1qzZo2ZmppqPvXUU+auXbvMG264wezXr5/Z0NBgmqZpXnPNNeacOXO85V977TWzV69e5gMPPGDu3r3bvOeeexhmbEG4x3nhwoVm7969zYqKCnP//v3e5fDhw4l6Ckkh3OPcEaN4rAn3OH/66adm3759zVtvvdXcs2ePuXHjRjMnJ8e89957E/UUkkK4x/mee+4x+/btaz777LPmRx99ZP7xj380hw0bZl5++eWJegpJ4fDhw+bbb79tvv3226Yk88EHHzTffvtt85NPPjFN0zTnzJljXnPNNd7ynmHGt99+u7l7925z+fLlDDOOlYceesg89dRTzd69e5ujR482//rXv3pv+973vmdOmzbNr/xzzz1nnnHGGWbv3r3NkSNHmn/4wx/iXOPkFM5xHjx4sCmp03LPPffEv+JJJtz3sy8CinXhHufXX3/dHDNmjJmammqedtpp5q9//Wvz+PHjca518gnnOH/11VfmL37xC3PYsGFmWlqaWVhYaN5yyy3mP/7xj/hXPIlUV1cH/L71HNtp06aZ3/ve9zrd5+yzzzZ79+5tnnbaaeaTTz4Z83oapkk7GAAAcJYedQ4KAABIDgQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOP8fyoeT2U0GM4sAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "m.visualize()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The result is identical when we fit with `iminuit.cost.ExtendedBinnedNLL`." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Migrad
FCN = 92.87 (chi2/ndof = 1.0) Nfcn = 190
EDM = 1.67e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0_n 990 40 0
1 x0_mu 0.4964 0.0019 0 1
2 x0_sigma 0.0487 0.0016 0
3 x1_n 629 29 0
4 x1_mu 0.97 0.14 0
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
x0_n x0_mu x0_sigma x1_n x1_mu
x0_n 1.3e+03 0.000424 (0.006) 0.0123 (0.209) -264 (-0.252) -0.369 (-0.076)
x0_mu 0.000424 (0.006) 3.44e-06 -9.8e-08 (-0.032) 0.000293 (0.005) -1.64e-05 (-0.065)
x0_sigma 0.0123 (0.209) -9.8e-08 (-0.032) 2.64e-06 -0.0129 (-0.271) -1.64e-05 (-0.075)
x1_n -264 (-0.252) 0.000293 (0.005) -0.0129 (-0.271) 851 0.337 (0.085)
x1_mu -0.369 (-0.076) -1.64e-05 (-0.065) -1.64e-05 (-0.075) 0.337 (0.085) 0.0183
" + ], + "text/plain": [ + "┌─────────────────────────────────────────────────────────────────────────┐\n", + "│ Migrad │\n", + "├──────────────────────────────────┬──────────────────────────────────────┤\n", + "│ FCN = 92.87 (chi2/ndof = 1.0) │ Nfcn = 190 │\n", + "│ EDM = 1.67e-05 (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 │ x0_n │ 990 │ 40 │ │ │ 0 │ │ │\n", + "│ 1 │ x0_mu │ 0.4964 │ 0.0019 │ │ │ 0 │ 1 │ │\n", + "│ 2 │ x0_sigma │ 0.0487 │ 0.0016 │ │ │ 0 │ │ │\n", + "│ 3 │ x1_n │ 629 │ 29 │ │ │ 0 │ │ │\n", + "│ 4 │ x1_mu │ 0.97 │ 0.14 │ │ │ 0 │ │ │\n", + "└───┴──────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", + "┌──────────┬───────────────────────────────────────────────────┐\n", + "│ │ x0_n x0_mu x0_sigma x1_n x1_mu │\n", + "├──────────┼───────────────────────────────────────────────────┤\n", + "│ x0_n │ 1.3e+03 0.000424 0.0123 -264 -0.369 │\n", + "│ x0_mu │ 0.000424 3.44e-06 -9.8e-08 0.000293 -1.64e-05 │\n", + "│ x0_sigma │ 0.0123 -9.8e-08 2.64e-06 -0.0129 -1.64e-05 │\n", + "│ x1_n │ -264 0.000293 -0.0129 851 0.337 │\n", + "│ x1_mu │ -0.369 -1.64e-05 -1.64e-05 0.337 0.0183 │\n", + "└──────────┴───────────────────────────────────────────────────┘" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def total(xe, x0_n, x0_mu, x0_sigma, x1_n, x1_mu):\n", + " return signal(xe, x0_n, x0_mu, x0_sigma) + background(xe, x1_n, x1_mu)\n", + "\n", + "c = ExtendedBinnedNLL(n, xe, total)\n", + "m = Minuit(c, x0_n=500, x0_mu=0.5, x0_sigma=0.05, x1_n=100, x1_mu=1)\n", + "m.limits[\"x0_n\", \"x1_n\", \"x0_sigma\", \"x1_mu\"] = (0, None)\n", + "m.limits[\"x0_mu\"] = (0, 1)\n", + "m.migrad()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/SklEQVR4nO3de3xU9Z3/8fdkMAnkxiWQCwmEuhQBr9WKgNBa8oAu4sIvZNUKFl1W7IpdAusFUJS1FpS6FlSU4lK1D6VU0xEqsrQWq6SC2KW69VZqK9Ukksg1IYMJODm/P3DGTDIzOTNzZuZM5vV8PObxyDnznTOf880M+fC9OgzDMAQAAGAjaYkOAAAAoDMSFAAAYDskKAAAwHZIUAAAgO2QoAAAANshQQEAALZDggIAAGyHBAUAANhOr0QHEIn29nZ98sknysnJkcPhSHQ4AADABMMwdPz4cRUXFystLXQbSVImKJ988olKS0sTHQYAAIhAbW2tSkpKQpZJygQlJydH0ukbzM3NTXA0AADAjObmZpWWlvr+joeSlAmKt1snNzeXBAUAgCRjZngGg2QBAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKAACwHRIUAABgOyQoAADAdkhQAACA7ZCgAAAA2yFBAQAAtkOCAgAAbIcEBQAA2A4JCgAAsB0SFAC243a75XA45HA45Ha7Ex0OgAQgQQEAALZDggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKAACwHRIUAABgOyQoAADAdkhQAACA7ZCgAAAA2yFBAQAAtkOCAgAAbIcEBQAA2A4JCgAAsB0SFAAAYDskKAAAwHZIUAAAgO2QoAAAANsJO0HZuXOnrrjiChUXF8vhcGjz5s1+zxuGobvuuktFRUXq3bu3ysvL9cEHH/iVOXLkiGbNmqXc3Fz17dtXc+fOVUtLS1Q3AgAAeo6wExS3263zzjtPa9euDfj8qlWr9NBDD2ndunXas2ePsrKyNGXKFLW2tvrKzJo1S++++65eeuklbd26VTt37tS8efMivwsAANCjOAzDMCJ+scOh559/XjNmzJB0uvWkuLhY//Ef/6FbbrlFktTU1KSCggI9+eSTuvrqq/X+++9r1KhR+sMf/qCLLrpIkrR9+3ZNnTpVdXV1Ki4u7vZ9m5ublZeXp6amJuXm5kYaPgCb8n7HJWnbtm2aPHmynE5ngqMCEK1w/n5bOgZl//79amhoUHl5ue9cXl6exowZo927d0uSdu/erb59+/qSE0kqLy9XWlqa9uzZE/C6bW1tam5u9nsA6JlcLpdGjRrlO546darKysrkcrkSGBWAeLM0QWloaJAkFRQU+J0vKCjwPdfQ0KBBgwb5Pd+rVy/179/fV6azlStXKi8vz/coLS21MmwANuFyuVRZWan6+nq/8/X19aqsrCRJAVJIUsziWbJkiZqamnyP2traRIcEwGIej0cLFixQoF5n77mqqip5PJ54hwYgASxNUAoLCyVJjY2NfucbGxt9zxUWFurTTz/1e/7zzz/XkSNHfGU6y8jIUG5urt8DQM9SU1Ojurq6oM8bhqHa2lrV1NTEMSoAiWJpgjJs2DAVFhZqx44dvnPNzc3as2ePxo4dK0kaO3asjh07pr179/rKvPzyy2pvb9eYMWOsDAdAEjlw4ICl5QAkt17hvqClpUV//etffcf79+/XW2+9pf79+2vIkCGqqqrSvffeq+HDh2vYsGFatmyZiouLfTN9Ro4cqW9/+9u64YYbtG7dOp06dUo333yzrr76alMzeAD0TEVFRZaWA5Dcwp5m/Morr+iyyy7rcn7OnDl68sknZRiG7r77bq1fv17Hjh3TpZdeqkcffVRf/epXfWWPHDmim2++WS+88ILS0tI0c+ZMPfTQQ8rOzjYVA9OMgZ7H4/GorKxM9fX1AcehOBwOlZSUaP/+/Uw5BpJUOH+/o1oHJVFIUICeyTuLR5JfkuJwOCRJ1dXVqqioSEhsAKKXsHVQACAaFRUVqq6u7tLdW1JSQnICpBhaUADYDivJAj0TLSgAbMvtdsvhcMjhcMjtdgcs0zEZmThxIskJkIJIUAAAgO2QoAAAANshQQEAALZDggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKAACwHRIUAABgOyQoAADAdkhQAACA7ZCgAAAA2yFBAQAAtkOCAgAAbIcEBQAA2A4JCgAAsB0SFAAAYDskKAAAwHZIUADElcfj8f28c+dOv2MA8CJBARA3LpdLo0aN8h1PnTpVZWVlcrlcCYwKgB2RoACIC5fLpcrKStXX1/udr6+vV2VlpV+SkpWVJcMwZBiGsrKy4h0qABsgQQEQcx6PRwsWLJBhGF2e856rqqqiuweADwkKgJirqalRXV1d0OcNw1Btba1qamriGBUAOyNBARBzBw4cMFXusssuk9vtDvic2+2Ww+GQw+EIWgZAz0GCAiDmioqKEh0CgCRDggIg5iZMmKCSkhI5HI6Azwc7DyB1kaAAiDmn06k1a9ZI6pqMOBwOv8GzI5dtj2tsAOyJBAVAXFRUVKi6ulrFxcV+50tKSjRg2i0JigqAXZGgAIibiooKvffee77jbdu2af/+/eoz/JIERgXAjkhQAMSV0+n0/Txx4kS/YwDw6pXoAAD0bG63W9nZ2ZKklpYWU68pW/xil3PvLvumlWEBsDlaUAAAgO2QoAAAANuhiwdAXHWcRjxy2XalpWcmMBoAdkULCgAAsB0SFAAAYDskKAAAwHZIUAD0WOyADCQvEhQAAGA7JCgAAMB2SFAAAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYjuUJisfj0bJlyzRs2DD17t1bZ555pn7wgx/IMAxfGcMwdNddd6moqEi9e/dWeXm5PvjgA6tDAdCDdN7DB0DPZnmCcv/99+uxxx7TI488ovfff1/333+/Vq1apYcffthXZtWqVXrooYe0bt067dmzR1lZWZoyZYpaW1utDgcAACQhyxOUXbt2afr06br88stVVlamyspKTZ48WW+88Yak060nq1ev1p133qnp06fr3HPP1c9+9jN98skn2rx5s9XhAEgCRrvH93Nr7Tt+x8HKeDxdywDoOSxPUMaNG6cdO3boL3/5iyTp//7v//T73/9e//iP/yhJ2r9/vxoaGlReXu57TV5ensaMGaPdu3cHvGZbW5uam5v9HgB6hhP7dunAhpt8xwerl6t+3Vyd2LcrZJmysjK5XK64xgogfixPUBYvXqyrr75aZ511ls444wxdcMEFqqqq0qxZsyRJDQ0NkqSCggK/1xUUFPie62zlypXKy8vzPUpLS60OG0CMhBo7cuKD13Vw8wp5Wg77nfccP6SDm1foxL5dOrFvV8Ay9fX1qqysJEkBeijLE5Rnn31WzzzzjDZu3Kg//vGPeuqpp/TAAw/oqaeeiviaS5YsUVNTk+9RW1trYcQAEuXYK0+EfP7wb3+iI7/9ScDnvAPvq6qq6O4BeqBeVl/w1ltv9bWiSNI555yjjz76SCtXrtScOXNUWFgoSWpsbFRRUZHvdY2NjTr//PMDXjMjI0MZGRlWhwogwTq3inTW3s3zhmGotrZWv/nNb3zdyAB6BstbUE6cOKG0NP/LOp1Otbe3S5KGDRumwsJC7dixw/d8c3Oz9uzZo7Fjx1odDoAUEKx7GEDysrwF5YorrtAPf/hDDRkyRKNHj9abb76pBx98UP/yL/8iSXI4HKqqqtK9996r4cOHa9iwYVq2bJmKi4s1Y8YMq8MBYDNp6ZkaevtWtX78JzX+fKkl1/S2zALoOSxPUB5++GEtW7ZMN910kz799FMVFxfrxhtv1F133eUrc9ttt8ntdmvevHk6duyYLr30Um3fvl2ZmZlWhwMgwTpPD+497AI50pzKKBktZ06+PMcPBX1tWvYAOdR9V9D48eOtCheATTiMjku8Jonm5mbl5eWpqalJubm5iQ4HQBAul0tXzpnnl2A4c/LVf9I89RkxzjdDJ5iBM063sAQq43A4fANlW1palJWV1aWM2+1WdnZ2yDIA4iecv9/sxQMgJlwulyorK0NOIe4zYpwGzlgqZ/YAvzLOnHwNnLFUfUaMC1pm8ODBMb8HAIlDggLAch6PRwsWLFCoBtojO9bLaPeoz4hxKpr7qO/8wMrlGvy9DeozYpzvXKAy7777bmyCB2ALJCgALFdTU6O6urqQZTzHD6mt7nSS4Uhz+s5nlp7td+zVuYzT2bUMgJ6DBAWA5Q4cOGCqnKflaIwjAZCsSFAAWK7jIoyhOLP7xTSOjivM7ty5kxVngSRCggLAchMmTFBJSYkcDkfQMs6cfGWUjI5ZDC6XS6NGjfIdT506lQ0GgSRCggLAck6nU2vWrAlZpv+keQHHmljBO4Oovr7e7zwbDALJgwQFQExUVFSouro65BTiWAg1g4gNBoHkQYICIGYqKiq6nUJste5mEHk3GKypqYlZDACiR4ICIKbMTCG2ktkZRGbLAUgMEhQAPYrZGURmywFIDBIUAD1KdzOIHA6HSktLNWHChDhHBiAcJCgAktLIZdv9fi5b/KIk/xlEnZMU7/Hq1atZiRawORIUAD2OdwZRcXGx3/mSkhJVV1eroqIiQZEBMKtXogMAgFioqKhQeXm58vLyJEnbtm3T5MmTaTkBkgQtKABMcbvdcjgccjgccrvdEZeJpc7v3zEZmThxIskJkERIUAAAgO2QoABISkb7lyvBtta+43cMIPkxBgWAZTrPrElLz4zJ+5z44HUde+UJ3/HB6uVy5uRry/mh9/8BkDxIUAAkncNbH+hyznP8kGbPnp2AaADEAl08AADAdkhQAPQYgXYwBpCcSFAAAIDtkKAAAADbIUEB0GME2yAQQPJhFg+AhEtLz9TQ27eaKnNi3y4d3LwiTpEBSBRaUAAklT4jxmngjKVyZg/wO+/MydfTTz+doKgAWI0WFABJp8+IccoYep7q1lwlSRpYuVy9h12g6dMnJTgyAFahBQVAUnKkfbnxX2bp2X7HAJIfCQoAALAdEhQAAGA7jEEBEFNmZugAQGe0oAAAANshQQEAALZDggIAAGyHBAUAANgOCQoAUzwej+/nnTt3+h17Ge1fnmutfcfvOB7MxAggOZCgAOiWy+XSqFGjfMdTp05VWVmZXC6XX5kDG27yHR+sXq76dXN1Yt+uuMR4Yt+uLjGOHj1av/zlL2UYhrKysuISBwBrOAzDMBIdRLiam5uVl5enpqYm5ebmJjocoEdzuVyqrKxU538qvDsHV1dXS1LAMl4DZyxVnxHjLI2r/WSran9cKUkaMO0WHd76QJcyHWOsqKiw9P0BhC+cv98kKACC8ng8KisrU11dXcDnHQ6HBg8eLElBy0inN/Ib/L0Nli5H3zFBcWYPkKflcNAYS0pKtH//fjmdLIcPJFI4f7/p4gEQVE1NTcjEwzAM1dXVhSwjSZ7jh9RW967V4X15/SDJiXQ6xtraWtXU1MTs/QFYjwQFQFAHDhyw7FqelqOWXSsSVt4LgNgjQQEQVFFRkWXXcmb3s+xakbDyXgDEHgkKgKAmTJigkpIS32DTzrzjO0KVkU6PQckoGR2rMOXMHhD0OYfDodLSUk2YMCFm7w/AeiQoAIJyOp1as2aNJHVJQLzHa9as8ZUJpv+keZYOkO2s7zevD3jeG+Pq1asZIAskGRIUACFVVFSourpaxcXFfudLSkp803e9ZTq3ZDhz8mMyxbizPsMv0cAZS7u8f8cYASSXXokOAID9VVRUqLy8XHl5eZKkbdu2afLkyX6tEhUVFSra+bnq1lwlSRpYuVy9h10Q05aTjvqMGKeMoef5vf/+TXfScgIkKVpQAJjS8Q/9xIkTA/7h75iMZJaeHbfkJNj7k5wAyYsWFAARK1v8YqJDANBDkaAASEpp6ZkaevvWRIcBIEbo4gEAALZDggIAAGyHBAUAANhOTBKU+vp6zZ49WwMGDFDv3r11zjnn6H//9399zxuGobvuuktFRUXq3bu3ysvL9cEHH8QiFAAx1H6yVR/dP00f3T9N7SdbEx0OgB7E8gTl6NGjGj9+vM444wz9z//8j9577z3913/9l/r1+3IfjlWrVumhhx7SunXrtGfPHmVlZWnKlClqbeUfOAAAEINZPPfff79KS0v1xBNP+M4NGzbM97NhGFq9erXuvPNOTZ8+XZL0s5/9TAUFBdq8ebOuvvpqq0MCAABJxvIWlF/96le66KKL9M///M8aNGiQLrjgAj3++OO+5/fv36+GhgaVl5f7zuXl5WnMmDHavXt3wGu2tbWpubnZ7wEAAHouyxOUDz/8UI899piGDx+uX//61/q3f/s3/fu//7ueeuopSVJDQ4MkqaCgwO91BQUFvuc6W7lypfLy8nyP0tJSq8MGAAA2YnmC0t7erq997WtasWKFLrjgAs2bN0833HCD1q1bF/E1lyxZoqamJt+jtrbWwogBAIDdWJ6gFBUVadSoUX7nRo4cqY8//liSVFhYKElqbGz0K9PY2Oh7rrOMjAzl5ub6PQAAQM9leYIyfvx47du3z+/cX/7yFw0dOlTS6QGzhYWF2rFjh+/55uZm7dmzR2PHjrU6HAAAkIQsn8WzcOFCjRs3TitWrNCVV16pN954Q+vXr9f69eslSQ6HQ1VVVbr33ns1fPhwDRs2TMuWLVNxcbFmzJhhdTgAACAJWZ6gfP3rX9fzzz+vJUuW6J577tGwYcO0evVqzZo1y1fmtttuk9vt1rx583Ts2DFdeuml2r59uzIzM60OBwAAJKGY7GY8bdo0TZs2LejzDodD99xzj+65555YvD0AAEhy7MUDAABshwQFAADYDgkKAACwnZiMQQHQ82RlZckwjJBl0tIzNfT2rXGKyH7vD8A6tKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKgIgZ7R7fz6217/gdA0A0HEZ3w/JtqLm5WXl5eWpqamJnYyBBXC6XrpwzT56Ww75zzpx89Z80T31GjEtgZKH9/b7LEx0CkLLC+ftNCwqAsLlcLlVWVvolJ5LkOX5IBzev0Il9uxIUGYCeggQFQFg8Ho8WLFgQck2UIzvW090DICokKADCUlNTo7q6upBlPMcPqa3u3ThFBKAnIkEBEJYDBw6YKudpORrjSAD0ZCQoAMJSVFRkqpwzu1+MIwHQk5GgAAjLhAkTVFJSIofDEbSMMydfGSWj4xgVgJ6GzQIBmFK2+EXfz20XfVdG3YqgZftPmidHmjMeYQHooWhBARC2PiPGaeCMpXJmD/A778zJ18AZS229DgqA5EALCoCI9BkxThlDz1PdmqskSQMrl6v3sAtoOQFgCVpQAESsYzKSWXo2yQkAy5CgAAAA2yFBAQAAtkOCAgAAbIcEBQAA2A4JCgAAsB0SFAAAYDskKAAAwHZIUAAE5Ha75XA45HA45Ha7Ex0OgBRDggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtkKAASClMnwaSAwkKAACwHRIUAABgO70SHQCA5JWWnqmht29NdBgAeiBaUAAAgO2QoAAAANshQQEAALZDggIAAGyHBAVAShm5bHvAnwHYCwkKgIA8Ho/v5507d8po94QoDQDWIkEB0IXL5dKoUaN8x1OnTlX9urk6sW9XAqMCkEpIUAD4cblcqqysVH19vd95z/FDOrh5BUkKgLggQQHg4/F4tGDBAhmGEbTMkR3r6e4BEHMkKAB8ampqVFdXF7KM5/ghtdW9G6eIAKQqEhQAPgcOHDBVztNyNMaRAEh1JCgAfIqKikyVc2b3i3EkAFIdCQoAnwkTJqikpEQOhyNoGWdOvjJKRscxKmt1HD/TWvuO33RqAPZBggLAx+l0as2aNSHL9J80T440Z5wistaJfbt0YMNNvuOD1ctVVlYml8uVwKgABEKCAkBut1sOh0MOh0NTpkxRdXW1nNkD/Mo4c/I1cMZS9RkxLkFRRufEvl06uHmFPC2H/c7X19ersrKSJAWwGRIUAF1UVFSoaO6jvuOBlcs1+HsbkjY5Mdo9OrJjfeDnvphSXVVVRXcPYCMxT1Duu+8+ORwOVVVV+c61trZq/vz5GjBggLKzszVz5kw1NjbGOhQAYejYjZNZenbSdutIUlvdu/IcPxT0ecMwVFtbq5qamjhGBSCUmCYof/jDH/STn/xE5557rt/5hQsX6oUXXtBzzz2nV199VZ988okqKipiGQqAFGZ2WrTZadYAYi9mCUpLS4tmzZqlxx9/XP36fTklsampSRs2bNCDDz6ob33rW7rwwgv1xBNPaNeuXXr99ddjFQ6AFGZ2WrTZadYAYi9mCcr8+fN1+eWXq7y83O/83r17derUKb/zZ511loYMGaLdu3cHvFZbW5uam5v9HgBgVkbJaDlz8oM+73A4VFpaqgkTJsQxKgChxCRB2bRpk/74xz9q5cqVXZ5raGhQenq6+vbt63e+oKBADQ0NAa+3cuVK5eXl+R6lpaWxCBtAD+VIc6r/pHmBn/tizZfVq1fL6UzecTZAT2N5glJbW6sFCxbomWeeUWZmpiXXXLJkiZqamnyP2tpaS64LIHX0GTFOA2cs7TJ9uqSkRNXV1YyDA2yml9UX3Lt3rz799FN97Wtf853zeDzauXOnHnnkEf3617/WyZMndezYMb9WlMbGRhUWFga8ZkZGhjIyMqwOFUCK6TNinDKGnqe6NVdJOj19ev+mO2k5AWzI8gRl0qRJevvtt/3OXX/99TrrrLN0++23q7S0VGeccYZ27NihmTNnSpL27dunjz/+WGPHjrU6HADw03n6NMkJYE+WJyg5OTk6++yz/c5lZWVpwIABvvNz587VokWL1L9/f+Xm5ur73/++xo4dq0suucTqcAAAQBKyPEEx48c//rHS0tI0c+ZMtbW1acqUKXr00Ue7fyEAAEgJcUlQXnnlFb/jzMxMrV27VmvXro3H2wMAgCTDXjwAAMB2SFAA+G2St3PnTjbNA5BwJChAinO5XBo1apTveOrUqSorK9OJD9h6AkDiJGSQLAB7cLlcqqyslGEYfufr6+tl1D2QoKgAgBYUIGV5PB4tWLCgS3IiKeA5AIgnEhQgRdXU1Kiurq7bcgXfWaG0dGu2rQAAs0hQgBR14MABU+U8LUdjHAkAdMUYFCBFFRUVmSrnzO4X40gSq2zxi13O/f2+yxMQCYCOaEEBUtSECRNUUlIih8MRtIwzJ18ZJaPjGBUAnEaCAqQop9OpNWvWSFKXJMV73H/SPL/N9QAgXkhQgBRWUVGh6upqFRcX+50vKSnRwBlL1WfEuARFBiDVMQYFSHEVFRVasPNzac1VkqSBlcvlGHaB+vTQlpO09EwNvX1rosMA0A1aUAD4deNklp5Ntw6AhCNBAQAAtkOCAgAAbIcEBUBKaz/Zqo/un6aP7p+m9pOtiQ4HwBdIUAAAgO2QoAAAANshQQEAALZDggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtsJsxAHb4BWA7tKAAAADbIUEBkNKMdo/v59bad/yOASQOCQqAlHVi3y4d2HCT7/hg9XLVr5srl8uVwKgASCQoQEpyu91yOBxyOBxyu92JDichTuzbpYObV8jTctjvvOf4IVVWVpKkAAlGggIg5RjtHh3ZsT5kmaqqKnk8dPcAiUKCAiDltNW9K8/xQ0GfNwxDtbW1qqmpiWNUADoiQQGQcjwtR02VO3DgQIwjARAMCQqAlOPM7meq3MIX/q6yxS/GOBoAgZCgAEg5GSWj5czJD1nGmZOvjJLRcYoIQGckKABSjiPNqf6T5oUs03/SPDnSnHGKCEBnJChAChq5bHvAn1NJnxHjNHDGUjmzB/idd+bka+CMpeozYpzvHNOygfhjLx4AKavPiHHKGHqe6tZcJUkaWLlcvYddQMsJYAO0oABIaR2TkczSs0lOAJsgQQEAALZDggL0cIHGT7BBHgC7YwwKkGJcLleXDfKcOfnqP2me38BQAEgkWlCAFLJlyxZVVlYG3CDv4OYVOrFvV4IiAwB/JChACrnttttkGEbQ54/sWE93j6T2k6366P5p+uj+aWo/2ZrocICURIICpJD6+vqQz3uOH1Jb3btxigYAgiNBAeDH7EZ6ABBLJCgA/JjdSC+VeDxfdnvt3LnT7xhAbJCgAClk8ODBcjgcQZ9ng7yuTnzwukaNGuU7njp1qsrKyuRyuRIYFdDzkaAAKWTVqlUhn2eDvK4Ob32gy9id+vp6VVZWkqQAMUSCAvRwHbsj+vXrp2effdbUBnmpIi09U0Nv36qht29VWnqmqdd4Z0JVVVXR3QPECAkK0MN0XDl248aNXbonFi5cqLxvzPGdG1i5XIO/tyElk5NoGIah2tpa9erVix2OgRhgJVmgB5s9e3aXdU/q6+tl1D3oO2aDPAB2RAsK0IMFWpQt1EJtAGAXlicoK1eu1Ne//nXl5ORo0KBBmjFjhvbt2+dXprW1VfPnz9eAAQOUnZ2tmTNnqrGx0epQACAiZlbTDTUbCkD0LE9QXn31Vc2fP1+vv/66XnrpJZ06dUqTJ0/266NduHChXnjhBT333HN69dVX9cknn6iiosLqUOIq0I6xAJLPiX27/DZTBJAYlo9B2b59u9/xk08+qUGDBmnv3r2aOHGimpqatGHDBm3cuFHf+ta3JElPPPGERo4cqddff12XXHKJ1SEBgCkn9u3Swc0rui3nzMnXz9at0axZs+IQFZCaYj4GpampSZLUv39/SdLevXt16tQplZeX+8qcddZZGjJkiHbv3h3rcICUQjeEeUa7R0d2rO+2XP7MuzT4exs0ffr0OEQFpK6YzuJpb29XVVWVxo8fr7PPPluS1NDQoPT0dPXt29evbEFBgRoaGgJep62tTW1tbb7j5ubmmMUMIDW11b0rz/FD3ZZzONKY9QTEQUxbUObPn6933nlHmzZtiuo6K1euVF5enu9RWlpqUYRAz/b0008HXJRtwLRbEhSRfZndJNHjPhbbQABIimGCcvPNN2vr1q363e9+p5KSEt/5wsJCnTx5UseOHfMr39jYqMLCwoDXWrJkiZqamnyP2traWIUNJL3OK8cWXv+I79i3KNtwxnp1ZnaTRGdWX0lsIAjEmuUJimEYuvnmm/X888/r5Zdf1rBhw/yev/DCC3XGGWdox44dvnP79u3Txx9/rLFjxwa8ZkZGhnJzc/0eALpyuVxdVo5teOJm3zGLsgWXUTJazpz87ssNHqkT+3axgSAQY5aPQZk/f742btyoLVu2KCcnxzeuJC8vT71791ZeXp7mzp2rRYsWqX///srNzdX3v/99jR07lhk8QAdut1vZ2dmSpJaWFmVlZYUs73K5VFlZ2WUhNk/L4S5lvfvP4EuONKf6T5rX7Syez/72Bx3e+kCX8/X19Zo5c6Ykc78vAKFZ3oLy2GOPqampSd/85jdVVFTke/ziF7/wlfnxj3+sadOmaebMmZo4caIKCwuT/n8eNPcikTwejxYsWNDtKrFmFiBLZX1GjNPAGUu7jtvpcHzslScCvrZj3fP9B6IXky6eQI/rrrvOVyYzM1Nr167VkSNH5Ha75XK5go4/SQaBmtVp7kU81dTUqK6urttybfXvxyGa5NZnxDgVzX3UdzywcrnfcaAWqc5ee+21mMQGpBL24omSt1m9vr7e73x9fb0qKytJUhAXBw4cMFWOGSjmdBynE8m4nWBLJgAwjwQlCqGa1b3nqqqqaO5FzBUVFZkq552BgthK5hZhwC5IUKLQXbO6YRiqra1VTU1NHKNCKpowYYJKSkq6XTk2Y/DIOEXUs3gHFQ+5dYupmT7jx4+PQ1RAz0aCEgWzzepmywGRcjqdWrNmjaTQy9szxTg63pk+Utd67njsdFLPQLRIUKJgtlndbDkgGhUVFaqurlZxcbHf+c4zUhAd70yfzvU8ePDgBEUE9EwkKAG43W45HA45HA653e6g5bprVnc4HCotLdWECRPCuq5V8SG5RTJ1vaKiQu+9957veNu2bRr8bz/V0Nu3aujtW5WWnhmTWFNNnxHjutTzn/70J98xSw0A0SNBiUKoZnXv8erVq2nuRdginbpetvhFnb38Jd/xjS+doFsnRjp+r48ePapzzjnHd8xSA0D0SFCiFKxZvaSkRNXV1aqoqEhQZEhWTF1PDiOXbff9PGvWLH5fgMVIUCwQqFl9//79ppITumvQkZmp6zNnzuTzEmPeWTvRdIux1AAQHRIUi3Rs7p04cSLdOoiImanrSB4sNQBEjgQFsBGmpPdM/F6B8JGgBBCrjf+sui4bE/ZcTEnvmfi9AuEjQenEyo3/Oo4v2bhxY8DrbtmyJWHxmcX06PgxM3UdyaPzUgMAzCNB6SCWsydmz54d8LqzZ8+2RXywBzNT12FPLDUAWIsE5Qux3vgv1HXtEB/sI9TU9aeffrrb11sxAwXd61jPgVaWZakBIDokKF9I1MZ/HROOUONJ2JgwtQSbuj5t2jTfOcYf2UeglWXNLjUAIDASlC/YYeO/YONJ3G63LrvsMlPXCBVfMo4BScaYrdJ56vqWLVviPv4I5qXSUgOp/L1E/JCgfCHajf+ysrJkGIYMw1BWVlbEcUQ7noTZAj3Tli1bAo4/qqur08yZMzXo/92hssUvJig6ALAeCcoXwt34z4xImt+7G09iZXxmMT068W677baQY5aO7Fgvo536BNBzkKB8weqN/zpPBw5Hd+NJ4jlbwKppzYmYHt2TdG456cxz/JDa6t6NUzQIpOPePCOXbadFC4gSCUoHVm38F2w6cLgCjSd5+umnTcVnRR+xVdOaU216dKC6j0ef/aljn8bkugCQCCQonUSz8Z8UejpwuLzjSTp2hfTr109vv/12RPGF08Vi1bTmaK9Dt5B5zqy+iQ4BKYLvJeKBBCWAaEbjdzcd2Cs/P9/UeJJAXSPnnHNO2PGF28Vi1bTmaK6T6t1C3oHXn3/+ecjxUV4Zg0fGKTIEkirrz6T69xLxQ4JiMbPTkK+66ipJoceTBJu58cknn3S5XqguhGDX8XaxBHpdtNOuvfFEOj06nG6heHapJGJ6ZajxUR050nrutNZkZebzYrcpu6HisUt3rd3qDLFBgmIxs9N8L7/88pDjXaZPn95t14gUfKZQx/NmrtNZtNOuw9XxOsm+am6g5u9om8QrKiqUP32J0rL6+53veNxa+w4zeRAzyf69RPIhQbFYd9OVvcaPHx9yvIvZrqLXXnuty7nOC3odOnQo6OuDJSlWTrsePHhwWNdJ5lVzAzV/FxQUaNiwYX7nImkS7zNinIrmPuo7zh17pTpW68Hq5apfN1cn9u2K/AaAIJL5e4nkRIISQKDpgp0fwZjd7M07biTYeBezXSwNDQ1dzgXamDBcVk67XrVqVVjXscOqvpEI1vx9+PBhHTlyxO9cpE3iHbtxmnc/K0+L/3U9xw/p4OYVJCmwXLJ+L5G8SFAiFChp8T4WvZERsDl+8ODBpq9vtuuksLBQkn+3QqQziDp3PUQz7brzzKNnn33W9HXM3vs111wjt9ttWZeKmX7tYNcNd/aWd9XhmTNnqrm52dRrwsHCbfZw1h1f/mdm2HUPqLm5uctnrPNnKlAZq8ZcRPMZD7fbN5bjRGL1/Ya9kKAEYMVo/M7N8QMrlyvtO2t91x39g1dUtvjFoK01c7Y1y5mTH/T63q6RyZMnR7UoXEeBuh4imXYdqJtj4cKFuu+++0xdx0z3kleg/Wms6lIxc1/e65rtkgskUDddtFi4LfFO7NulAxtu8h0frF7e5Xsa6DNlxXc5UqE+47FYbdvqGNGzkKDEUMfm+MzSs8OaZeFIc6r/pHlBnzcMQ60XXquiyrs0c+bMqLt0vAJ1PYQz7TrUKP/vfve7pq5jtptMCtydZWWXild3sxe2bNkS0XWlwN10UtdWunB5Wo5GHBOic2LfLh3cvEKelsN+5zvOwAtnll48mPmMW7nadixiJEnpWUhQbKzPiHEaOGOpnNkD/M47c/I1cMZS9R4+Rkd2rLf0Pb1dD1ded6OG3varsF5rZpS/WaG6l55++umwrxvNLAMz9/XMM8+Edc2OvN10VnNm94vJdRGa0e4J+r3s+Bm69dZbI56lZzWzM3SmT59uyWrbsYyRWUQ9BwmKCe0nW/XR/dP00f3T1H6yNei5SHTXnRSoq2jw9zaoz4hxaqt7V57jwWfoRMNz/JA+2/9myG6ojo8hi36pXr16dTvKv6Pu1i+ZMmVKwO6ladOmRXRP3lkGv/nNbwLfc5B+bTOzFw4ePKiBAwd2O3srkPb2dtP/qKalZ2rIrVtCdv9Jp5PYjJLRYceC6Jn9XpppKXnssccsmabuFeg6brfb1HfXO0PHbLevmZjDGRcS7SyieK5+m6jtLnoaEpQkEKyrKNZN+B73sZhde+Sy7QETn87nzl7+ku/4xpdOdBlzEolAXSqh+rXNzkqYNWuWpNCLqQUybdo0ZfQr1KD/d4epLp3uuv8kqf+keSzcliBWfi8XL15s2ZiqYJ/xcLonvd+F7rp9YzFOJJpZRIxbSU4kKDEU66Wvw2nCd+bka8C0W8K7/hd7u1h9HwXfWRHRdU588Lol4206d6l016/9wQcfmLpusObvAQMGqH///kFedVq404O76/7rM2KcqevAelZ3rVkxpirUZ3z27NmmYzEzkydW40QiXTyScSvJiwTFhI7TNb2rdQY6F28ZJaO7beqXpPyZd53uFhp+ie9c5z9sgRhGu+n7MluuY9dDuPV67JUnTL1Hd7674XUNve1XKlv8oobe9itded2NQfu1DcPQf/7ooZD11XH2QkVFhRyVD/qeG1i5XFn/8lP1ufYxU7GFMz04VPcfEsfs9zIaVm7YaWYcl9kZOuGOEwmn2yWSWUSJGrcSi9WkUxEJSjcCTRWse3i2PvnJXL9zVq3gGc7YFjNN/ZLUe8i5XZr7+37z+m5fd+iX95i6r851FIq368FsvXYs03lGRKQ63peZ8QLtLYdDvrd3RtWZd5zupgrUJWe2u8Vz/JA+/tF00+Oaopkphtgw+72MVucxF8HGOEQzBV4Kb4aO2XEivXr10saNG8Pqduludp/32q2tX3534rX6bce6D3RfVi59kEpjWUhQQgg2VbC99bjaW1v8zsVyBc9QXSxBm/oD/I+/43WyR38z4Os66+6+gtVRZx27HsKq1y/K5Fw0PeT1O0rLzFFaZnbIMr77+mCP6esOmHYLXSowJdj3Mha6G5sR7cqugWboeHfaNgxDWVlZEb1XoCUCuut2CTa7L9gimIlY/TZeSx+kAhKUIEJNFQwl2hU8I+k6CtTU3/HY7OtCCXRfZuvI18U0YlzE9ep+93emyuVNnKOS7z+t4hs3WHpdSUrLzFbh9Y/4jmPZpcLGf8kvnO9XNLxjLqJdATaQrVu3drswY6BYzIi02yXQLKI//elPvuNoVr+NVCQreUfaxZRKXUUkKEFEOoU3mhU8A3V7mO06irSpP5yuh873ZbaOHI403/tEWq/tnzUrrXdut+VyLpgaVpeK2etKp38fDU/c7DuOZZcKG//1DLHscus45iKaFWBDSUtLC2vhtWjey8tMt0vHmI4ePapzzjnHdxzv1W+jWck73C6mVJuNRIISRDRTBSN5bbBuj0i7jmIxg8h7X95xMo0/X2rudR2mK0dTr1mjL+u2TCR/EMxc1ytQV1bncUOB6t57buAMc3Ummfvdx3qmGOyp47iQYCvSersQOq5xEm7i4J2Sb3bcQ6hxIuG67LLL5Ha7u33vUF1FsVj9tvN4k0B1Hy4zXUxWz0ZKhrEsJChBRDNVMNzXmun2sMPmb977CjcO73TljteIRJ/hY0yPt7Hiut3x1kM43XKRjE2ww+8e0Qs2hqnj9P9AZQKNqfKOC5k+fbrpGTpLly4Nuysj0Gak3XUrBBsnEomdO3fq5MmTId873qvfdnz/cDYIDaW730ssZiMlQ1cRCUoQkU4VjGQFTzPdHone/M17X+HM2PHKGDzyy5+jrNdIx9uEe10z2urfj6hbLtz3SvTvHtboM/ySwNPCO0z/D1Sm85iqgZXL5bj6ES16I0ODr11leobOD3/4Q7/j/Pz8bls5xo8fH1G3QudxIpGaOnWq8vPz/Y7NLi4Xyeq33elcF4cORbeSt9kuJqtnIyVLVxEJShCRThWMZAVPs90eocpF2tRvtuuh/6R5+uyDPaZm7HQ0cMZSOTO/HOVvRb12Hm/jzMyKuEsl2HXNaP37mxF3y8XqMwJ76fy9DDRWzEwZq1aT7ri8vnHJdZICd304HA798pe/1G9/+9uIuxWs2jSwvb29y3vPnj1bVVVVpl5vdvXb7gTrYolUOF1MVs5GSqaF60hQQgjWHB+oyTWa6aZmuz1iufmb1RsThqqPeNVrqPfyNq0f3Lwi4v2UupsBZGXXDBv/IZBoPhe9z/y68qcvUVqW/yrHadkDlD99iRa+3qvbRQy9m4p69+PyJjdDFv3Sb9sKK4W7Qec111xjeoxFsHEZobpYzAi0mvTgwYNlGIZmzpwZdL8e7/E111xj6n2i7SryxtPc3GzyzmKrV6IDsLs+I8YpY+h5qltzlaTTzau9h12g9pOtXc5FOmLf2+0RqpsnHpu/BbtXR5pTrR//yfTsGzP1EU29ev/HGe19dUxKWmvfUe9hF2jo7VtltHtUv25uyPtN652r9s9Cf4m9XTOZQ87t+vov7sHMe7HxX89h5rMbqEyw15n5tyOYU5/+LervfMfPeOexWOF8n8Ll3aDTzPdQkoZd94AObLqz2wSj87iMyZMny+l0RrXY3datW/Xtb39bbrdbeXl5kk53MV1yySW+pMX7Xp3f/5JLvuz+y8/P1+HDhwPeg8PhUElJSdRdRV6PPfaYbrnlFstawSJFgmJCOE2ukV6//6R5Orh5RdAy8dr8zYqmZLP1Eet6DfVen32wR0d++xPfuYPVy+XMyVf/SfPUZ8S4bn8fWaMv0/H/7b4fvLt6s9PvHsnHzOcnGO9nM9rvvKflqE7s2xXV9ylSZr+HB6uXK6PfI754Aul8D1OnTvXdg+E5FXGM/7ajVWk12/3+Q3Ttoy8r84Yb/N5rwIABfsnH1KlTlZb2ZSdHsPEusegqWrx4sR555BGtWbMm7HE6ViJBsQlvV8SR3/7Eb0xDxy95IvWkLoYTH7yuw1sf6HLeO3bE26UU6veR1jvb1D+MZurN7r972Fuwz093uvtsmv3Onzpar6bfb+xy3tT3KXtAVFtY9Bk+Rpklo03de+d4OvIu8xDsNXnjZ4UVV3f3Fejfn8OHAyxh0Gn8TSBp2QPUf9I8LXojQ4veCL4LuiS1fvz3bq/n5R2TEslMJ6s4DCvmSMVZc3Oz8vLy1NTUpNxcc4tshSPUVvex5ml1W9Z1ZCWzXRGDv7fBFvEGE+59BPt9xKI+7Pq7h/0E+qx07B7truvDzGfTVFdn9gA5FHqfLKu+T2avG869d3xdKGbu0yvUd9fq7q6BlcuVXjRC9Q9/J+R7h1M/HXm7jvbv329Zd084f78ZJGszdt38zczsm2Toigh3Snew30cs6sOuv3vYS7Cp7Z/97Q++c/0s+Gya+YznnPdtUy0XVnyfOgs2A8/s2LBwlkxobzms7POnmCob6rsb6UraQeNqbfFb3TrQMged79NsciJZt5lipGhBCSCRLSh25+2nTdauCPd7r+rQCz/qtlz+Fbcqa9Q3ui2X7PWB5BKsK8KrY/eFVZ/NUNcxPKdi/n2SI00yvuzqCHUPZr/fORdNN9VF21H+FbfK4Twjqjo1G58Z3d2Dd4kFK8b+bNy4Ud/5zneivo4U3t9vxqAgLKFG/ScDq6d0J3t9IHmYXXG69/AxcqQ5LftsdjfTx4xovk9mujDCfZ9wNgnteO3MIedGVadWjuXr7h4ObX9YDlnT/hDtZoqRoosHYUvmrggzK9mGO603mesDySOSFaet+mwGu048vk9pvc4I+N6BmIknnDEYXh3vIZo6jXQl7c7M3IPRelztrS1RvY8VmylG9f508XRFF0/PFk4zOWAXVndPWsVu36fu4omke8fKe+guPjMiuYdIruVwOCyfxcMgWSCE7lbNJTmBHdlhxelA7PZ96jae4WNMX8uZk6/86YuV1jtb7vdeVevHf4p6dehwVtKWw/9PdCT30G08wTZhzclP6BRjKcFjUNauXasf/ehHamho0HnnnaeHH35YF198cSJDQoroM2Kceg8fc7rZvOWonNn9lFEymu4Z2JZdVpwOxG7fp1DxGO2ebuvRkZmjgdNvV3tri46+/N9+Za0YBB8sPkl+59KLz9LJT/4c0T2Y4f28ONKcAeOpqPinqK4frYR18fziF7/Qd7/7Xa1bt05jxozR6tWr9dxzz2nfvn0aNGhQyNfSxQMgFdmtOyVZmalHKfQMmETXtRVdRd3dw9/vuzyq6weSFF08Dz74oG644QZdf/31GjVqlNatW6c+ffropz/9aaJCAgBb83UPdBpoSfdkeLqrRzObo1q5GWgkgt2DGcnyeUlIF8/Jkye1d+9eLVmyxHcuLS1N5eXl2r17d5fybW1tamtr8x03NTVJUsx2XGxvOxGT6wJAtDLLzlfR9Q+r7ZP35Wk5Jmd2X2UUj5Qjzcm/XWEIVY+f7d9rasbUZ/v3KrPk7DhF3FXHe/j8+BEde/UpGa3Hg78gI1v5l1cp84tune4+L7H4G+u9ppnOm4QkKIcOHZLH41FBQYHf+YKCAv35z3/uUn7lypX6z//8zy7nS0tLYxYjAAChHKy+J9EhhKetRYdc95ounrc6dqEcP37ct7tzMEmxUNuSJUu0aNEi33F7e7uOHDmiAQMG+HZytEpzc7NKS0tVW1sbk/EtOI16jg/qOT6o5/ignuMnVnVtGIaOHz+u4uLibssmJEHJz8+X0+lUY2Oj3/nGxkYVFhZ2KZ+RkaGMjAy/c3379o1liMrNzeULEAfUc3xQz/FBPccH9Rw/sajr7lpOvBIySDY9PV0XXnihduzY4TvX3t6uHTt2aOzYsYkICQAA2EjCungWLVqkOXPm6KKLLtLFF1+s1atXy+126/rrr09USAAAwCYSlqBcddVVOnjwoO666y41NDTo/PPP1/bt27sMnI23jIwM3X333V26lGAt6jk+qOf4oJ7jg3qOHzvUdVLuxQMAAHo29uIBAAC2Q4ICAABshwQFAADYDgkKAACwnZRMUNauXauysjJlZmZqzJgxeuONN0KWf+6553TWWWcpMzNT55xzjrZt2xanSJNbOPX8+OOPa8KECerXr5/69eun8vLybn8vOC3cz7PXpk2b5HA4NGPGjNgG2EOEW8/Hjh3T/PnzVVRUpIyMDH31q1/l3w4Twq3n1atXa8SIEerdu7dKS0u1cOFCtba2xina5LRz505dccUVKi4ulsPh0ObNm7t9zSuvvKKvfe1rysjI0D/8wz/oySefjHmcMlLMpk2bjPT0dOOnP/2p8e677xo33HCD0bdvX6OxsTFg+ddee81wOp3GqlWrjPfee8+48847jTPOOMN4++234xx5cgm3nq+55hpj7dq1xptvvmm8//77xnXXXWfk5eUZdXV1cY48uYRbz1779+83Bg8ebEyYMMGYPn16fIJNYuHWc1tbm3HRRRcZU6dONX7/+98b+/fvN1555RXjrbfeinPkySXcen7mmWeMjIwM45lnnjH2799v/PrXvzaKioqMhQsXxjny5LJt2zbjjjvuMFwulyHJeP7550OW//DDD40+ffoYixYtMt577z3j4YcfNpxOp7F9+/aYxplyCcrFF19szJ8/33fs8XiM4uJiY+XKlQHLX3nllcbll1/ud27MmDHGjTfeGNM4k1249dzZ559/buTk5BhPPfVUrELsESKp588//9wYN26c8d///d/GnDlzSFBMCLeeH3vsMeMrX/mKcfLkyXiF2COEW8/z5883vvWtb/mdW7RokTF+/PiYxtmTmElQbrvtNmP06NF+56666ipjypQpMYzMMFKqi+fkyZPau3evysvLfefS0tJUXl6u3bt3B3zN7t27/cpL0pQpU4KWR2T13NmJEyd06tQp9e/fP1ZhJr1I6/mee+7RoEGDNHfu3HiEmfQiqedf/epXGjt2rObPn6+CggKdffbZWrFihTweT7zCTjqR1PO4ceO0d+9eXzfQhx9+qG3btmnq1KlxiTlVJOrvYFLsZmyVQ4cOyePxdFmttqCgQH/+858DvqahoSFg+YaGhpjFmewiqefObr/9dhUXF3f5UuBLkdTz73//e23YsEFvvfVWHCLsGSKp5w8//FAvv/yyZs2apW3btumvf/2rbrrpJp06dUp33313PMJOOpHU8zXXXKNDhw7p0ksvlWEY+vzzz/W9731PS5cujUfIKSPY38Hm5mZ99tln6t27d0zeN6VaUJAc7rvvPm3atEnPP/+8MjMzEx1Oj3H8+HFde+21evzxx5Wfn5/ocHq09vZ2DRo0SOvXr9eFF16oq666SnfccYfWrVuX6NB6lFdeeUUrVqzQo48+qj/+8Y9yuVx68cUX9YMf/CDRocECKdWCkp+fL6fTqcbGRr/zjY2NKiwsDPiawsLCsMojsnr2euCBB3Tffffpt7/9rc4999xYhpn0wq3nv/3tb/r73/+uK664wneuvb1dktSrVy/t27dPZ555ZmyDTkKRfJ6Liop0xhlnyOl0+s6NHDlSDQ0NOnnypNLT02MaczKKpJ6XLVuma6+9Vv/6r/8qSTrnnHPkdrs1b9483XHHHUpL4//gVgj2dzA3NzdmrSdSirWgpKen68ILL9SOHTt859rb27Vjxw6NHTs24GvGjh3rV16SXnrppaDlEVk9S9KqVav0gx/8QNu3b9dFF10Uj1CTWrj1fNZZZ+ntt9/WW2+95Xv80z/9ky677DK99dZbKi0tjWf4SSOSz/P48eP117/+1ZcAStJf/vIXFRUVkZwEEUk9nzhxoksS4k0KDbaZs0zC/g7GdAiuDW3atMnIyMgwnnzySeO9994z5s2bZ/Tt29doaGgwDMMwrr32WmPx4sW+8q+99prRq1cv44EHHjDef/994+6772aasQnh1vN9991npKenG9XV1caBAwd8j+PHjyfqFpJCuPXcGbN4zAm3nj/++GMjJyfHuPnmm419+/YZW7duNQYNGmTce++9ibqFpBBuPd99991GTk6O8fOf/9z48MMPjd/85jfGmWeeaVx55ZWJuoWkcPz4cePNN9803nzzTUOS8eCDDxpvvvmm8dFHHxmGYRiLFy82rr32Wl957zTjW2+91Xj//feNtWvXMs04Vh5++GFjyJAhRnp6unHxxRcbr7/+uu+5b3zjG8acOXP8yj/77LPGV7/6VSM9Pd0YPXq08eKLL8Y54uQUTj0PHTrUkNTlcffdd8c/8CQT7ue5IxIU88Kt5127dhljxowxMjIyjK985SvGD3/4Q+Pzzz+Pc9TJJ5x6PnXqlLF8+XLjzDPPNDIzM43S0lLjpptuMo4ePRr/wJPI7373u4D/3nrrds6cOcY3vvGNLq85//zzjfT0dOMrX/mK8cQTT8Q8Todh0A4GAADsJaXGoAAAgORAggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKAACwnf8PnGX4Ai7lLKIAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "m.visualize()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "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.9.16" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "bdbf20ff2e92a3ae3002db8b02bd1dd1b287e934c884beb29a73dced9dbd0fa3" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/doc/tutorials.rst b/doc/tutorials.rst index 1c1a2bcb..5587a726 100644 --- a/doc/tutorials.rst +++ b/doc/tutorials.rst @@ -17,6 +17,7 @@ Important for most users are only the first two entries. notebooks/interactive notebooks/simultaneous_fits notebooks/template_fits + notebooks/template_model_mix notebooks/conditional_variable notebooks/scipy_and_constraints notebooks/external_minimizer diff --git a/src/iminuit/cost.py b/src/iminuit/cost.py index 93b19aa6..95c75913 100644 --- a/src/iminuit/cost.py +++ b/src/iminuit/cost.py @@ -48,7 +48,6 @@ the histogram, the sum of weights and the sum of squared weights is needed then, see class documentation for details. """ - from .util import ( describe, make_func_code, @@ -56,19 +55,49 @@ PerformanceWarning, _smart_sampling, ) +from .typing import Model, LossFunction import numpy as np -from collections.abc import Sequence +from numpy.typing import ArrayLike, NDArray +from collections.abc import Sequence as ABCSequence import abc -import typing as _tp +from typing import ( + List, + Tuple, + Union, + Sequence, + Collection, + Dict, + Any, + Iterable, + Optional, + overload, +) import warnings +__all__ = [ + "CHISQUARE", + "NEGATIVE_LOG_LIKELIHOOD", + "BohmZechTransform", + "chi2", + "multinominal_chi2", + "poisson_chi2", + "template_chi2_jsc", + "template_chi2_da", + "template_nll_asy", + "Cost", + "CostSum", + "Constant", + "BinnedNLL", + "UnbinnedNLL", + "ExtendedBinnedNLL", + "ExtendedUnbinnedNLL", + "Template", + "LeastSquares", +] + CHISQUARE = 1.0 NEGATIVE_LOG_LIKELIHOOD = 0.5 -# correct ArrayLike from numpy.typing generates horrible looking signatures -# in python's help(), so we use this as a workaround -_ArrayLike = _tp.Sequence - def _safe_log(x): # guard against x = 0 @@ -107,7 +136,10 @@ class BohmZechTransform: :meta private: """ - def __init__(self, val: _ArrayLike, var: _ArrayLike): + _scale: np.ndarray + _obs: np.ndarray + + def __init__(self, val: ArrayLike, var: ArrayLike): """ Initialize transformer with data value and variance. @@ -124,11 +156,17 @@ def __init__(self, val: _ArrayLike, var: _ArrayLike): np.divide(val, var, out=self._scale, where=var > 0) self._obs = val * self._scale + @overload + def __call__(self, val: ArrayLike) -> Tuple[NDArray, NDArray]: + ... + + @overload def __call__( - self, val: _ArrayLike, var: _ArrayLike = None - ) -> _tp.Union[ - _tp.Tuple[np.ndarray, np.ndarray], _tp.Tuple[np.ndarray, np.ndarray, np.ndarray] - ]: + self, val: ArrayLike, var: ArrayLike + ) -> Tuple[NDArray, NDArray, NDArray]: + ... + + def __call__(self, val, var=None): """ Return precomputed scaled data and scaled prediction. @@ -144,12 +182,14 @@ def __call__( (obs, pred) or (obs, pred, pred_var) """ s = self._scale + assert isinstance(val, np.ndarray) if var is None: return self._obs, val * s + assert isinstance(var, np.ndarray) return self._obs, val * s, var * s**2 -def chi2(y: _ArrayLike[float], ye: _ArrayLike[float], ym: _ArrayLike[float]) -> float: +def chi2(y: ArrayLike, ye: ArrayLike, ym: ArrayLike) -> float: """ Compute (potentially) chi2-distributed cost. @@ -174,7 +214,7 @@ def chi2(y: _ArrayLike[float], ye: _ArrayLike[float], ym: _ArrayLike[float]) -> return np.sum(_z_squared(y, ye, ym)) -def multinominal_chi2(n: _ArrayLike[float], mu: _ArrayLike[float]) -> float: +def multinominal_chi2(n: ArrayLike, mu: ArrayLike) -> float: """ Compute asymptotically chi2-distributed cost for binomially-distributed data. @@ -202,7 +242,7 @@ def multinominal_chi2(n: _ArrayLike[float], mu: _ArrayLike[float]) -> float: return 2 * np.sum(n * (_safe_log(n) - _safe_log(mu))) -def poisson_chi2(n: _ArrayLike[float], mu: _ArrayLike[float]) -> float: +def poisson_chi2(n: ArrayLike, mu: ArrayLike) -> float: """ Compute asymptotically chi2-distributed cost for Poisson-distributed data. @@ -229,9 +269,7 @@ def poisson_chi2(n: _ArrayLike[float], mu: _ArrayLike[float]) -> float: return 2 * np.sum(mu - n + n * (_safe_log(n) - _safe_log(mu))) -def template_chi2_jsc( - n: _ArrayLike[float], mu: _ArrayLike[float], mu_var: _ArrayLike[float] -) -> float: +def template_chi2_jsc(n: ArrayLike, mu: ArrayLike, mu_var: ArrayLike) -> float: """ Compute asymptotically chi2-distributed cost for a template fit. @@ -266,14 +304,10 @@ def template_chi2_jsc( p = 0.5 - 0.5 * mu * beta_var beta = p + np.sqrt(p**2 + n * beta_var) - return poisson_chi2(n, mu * beta) + np.sum( # type:ignore - (beta - 1) ** 2 / beta_var - ) + return poisson_chi2(n, mu * beta) + np.sum((beta - 1) ** 2 / beta_var) -def template_chi2_da( - n: _ArrayLike[float], mu: _ArrayLike[float], mu_var: _ArrayLike[float] -) -> float: +def template_chi2_da(n: ArrayLike, mu: ArrayLike, mu_var: ArrayLike) -> float: """ Compute asymptotically chi2-distributed cost for a template fit. @@ -297,12 +331,10 @@ def template_chi2_da( n, mu, mu_var = np.atleast_1d(n, mu, mu_var) k = mu**2 / mu_var beta = (n + k) / (mu + k) - return poisson_chi2(n, mu * beta) + poisson_chi2(k, k * beta) # type:ignore + return poisson_chi2(n, mu * beta) + poisson_chi2(k, k * beta) -def template_nll_asy( - n: _ArrayLike[float], mu: _ArrayLike[float], mu_var: _ArrayLike[float] -) -> float: +def template_nll_asy(n: ArrayLike, mu: ArrayLike, mu_var: ArrayLike) -> float: """ Compute marginalized negative log-likelikihood for a template fit. @@ -348,19 +380,19 @@ def template_nll_asy( # precision. Fall back to plain numpy for float128 which is not currently supported # by numba. try: - from numba import njit as _njit - from numba.extending import overload as _overload + from numba import njit as jit + from numba.extending import overload as nb_overload - @_overload(_safe_log, inline="always") + @nb_overload(_safe_log, inline="always") def _ol_safe_log(x): return _safe_log # pragma: no cover - @_overload(_z_squared, inline="always") + @nb_overload(_z_squared, inline="always") def _ol_z_squared(y, ye, ym): return _z_squared # pragma: no cover _unbinned_nll_np = _unbinned_nll - _unbinned_nll_nb = _njit( + _unbinned_nll_nb = jit( nogil=True, cache=True, error_model="numpy", @@ -373,15 +405,15 @@ def _unbinned_nll(x): return _unbinned_nll_np(x) _multinominal_chi2_np = multinominal_chi2 - _multinominal_chi2_nb = _njit( + _multinominal_chi2_nb = jit( nogil=True, cache=True, error_model="numpy", )(_multinominal_chi2_np) - def multinominal_chi2(n: _ArrayLike[float], mu: _ArrayLike[float]) -> float: # noqa - n, mu = np.atleast_1d(n, mu) # type:ignore - if mu.dtype in (np.float32, np.float64): # type:ignore + def multinominal_chi2(n: ArrayLike, mu: ArrayLike) -> float: # noqa + n, mu = np.atleast_1d(n, mu) + if mu.dtype in (np.float32, np.float64): return _multinominal_chi2_nb(n, mu) # fallback to numpy for float128 return _multinominal_chi2_np(n, mu) @@ -389,15 +421,15 @@ def multinominal_chi2(n: _ArrayLike[float], mu: _ArrayLike[float]) -> float: # multinominal_chi2.__doc__ = _multinominal_chi2_np.__doc__ _poisson_chi2_np = poisson_chi2 - _poisson_chi2_nb = _njit( + _poisson_chi2_nb = jit( nogil=True, cache=True, error_model="numpy", )(_poisson_chi2_np) - def poisson_chi2(n: _ArrayLike[float], mu: _ArrayLike[float]) -> float: # noqa - n, mu = np.atleast_1d(n, mu) # type:ignore - if mu.dtype in (np.float32, np.float64): # type:ignore + def poisson_chi2(n: ArrayLike, mu: ArrayLike) -> float: # noqa + n, mu = np.atleast_1d(n, mu) + if mu.dtype in (np.float32, np.float64): return _poisson_chi2_nb(n, mu) # fallback to numpy for float128 return _poisson_chi2_np(n, mu) @@ -405,17 +437,15 @@ def poisson_chi2(n: _ArrayLike[float], mu: _ArrayLike[float]) -> float: # noqa poisson_chi2.__doc__ = _poisson_chi2_np.__doc__ _chi2_np = chi2 - _chi2_nb = _njit( + _chi2_nb = jit( nogil=True, cache=True, error_model="numpy", )(_chi2_np) - def chi2( - y: _ArrayLike[float], ye: _ArrayLike[float], ym: _ArrayLike[float] - ) -> float: # noqa - y, ye, ym = np.atleast_1d(y, ye, ym) # type:ignore - if ym.dtype in (np.float32, np.float64): # type:ignore + def chi2(y: ArrayLike, ye: ArrayLike, ym: ArrayLike) -> float: # noqa + y, ye, ym = np.atleast_1d(y, ye, ym) + if ym.dtype in (np.float32, np.float64): return _chi2_nb(y, ye, ym) # fallback to numpy for float128 return _chi2_np(y, ye, ym) @@ -423,7 +453,7 @@ def chi2( chi2.__doc__ = _chi2_np.__doc__ _soft_l1_loss_np = _soft_l1_loss - _soft_l1_loss_nb = _njit( + _soft_l1_loss_nb = jit( nogil=True, cache=True, error_model="numpy", @@ -435,12 +465,12 @@ def _soft_l1_loss(z_sqr): # fallback to numpy for float128 return _soft_l1_loss_np(z_sqr) - @_overload(_soft_l1_loss, inline="always") + @nb_overload(_soft_l1_loss, inline="always") def _ol_soft_l1_loss(z_sqr): return _soft_l1_loss_np # pragma: no cover _soft_l1_cost_np = _soft_l1_cost - _soft_l1_cost_nb = _njit( + _soft_l1_cost_nb = jit( nogil=True, cache=True, error_model="numpy", @@ -513,7 +543,7 @@ def verbose(self): def verbose(self, value: int): self._verbose = value - def __init__(self, args: _tp.Tuple[str, ...], verbose: int): + def __init__(self, args: Sequence[str], verbose: int): """For internal use.""" self._func_code = make_func_code(args) self._verbose = verbose @@ -538,7 +568,7 @@ def __radd__(self, lhs): """ return CostSum(lhs, self) - def __call__(self, *args): + def __call__(self, *args: float) -> float: """ Evaluate the cost function. @@ -558,6 +588,10 @@ def __call__(self, *args): print(args, "->", r) return r + @abc.abstractmethod + def _call(self, args: Sequence[float]) -> float: + ... + class Constant(Cost): """ @@ -577,11 +611,11 @@ def __init__(self, value: float): def _ndata(self): return 0 - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: return self.value -class CostSum(Cost, Sequence): +class CostSum(Cost, ABCSequence): """Sum of cost functions. Users do not need to create objects of this class themselves. They should just add @@ -608,7 +642,7 @@ class CostSum(Cost, Sequence): __slots__ = "_items", "_maps" - def __init__(self, *items): + def __init__(self, *items: Union[Cost, float]): """Initialize with cost functions. Parameters @@ -616,7 +650,7 @@ def __init__(self, *items): *items : Cost Cost functions. May also be other CostSum functions. """ - self._items = [] + self._items: List[Cost] = [] for item in items: if isinstance(item, CostSum): self._items += item._items @@ -628,12 +662,12 @@ def __init__(self, *items): args, self._maps = merge_signatures(self._items) super().__init__(args, max(c.verbose for c in self._items)) - def _split(self, args): + def _split(self, args: Sequence[float]): for component, cmap in zip(self._items, self._maps): component_args = tuple(args[i] for i in cmap) yield component, component_args - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: r = 0.0 for comp, cargs in self._split(args): r += comp._call(cargs) / comp.errordef @@ -651,7 +685,7 @@ def __getitem__(self, key): return self._items.__getitem__(key) def visualize( - self, args: _ArrayLike, component_kwargs: _tp.Dict[int, _tp.Dict] = None + self, args: Sequence[float], component_kwargs: Dict[int, Dict[str, Any]] = None ): """ Visualize data and model agreement (requires matplotlib). @@ -674,8 +708,6 @@ def visualize( """ from matplotlib import pyplot as plt - args = np.atleast_1d(args) - n = sum(hasattr(comp, "visualize") for comp in self) w, h = plt.rcParams["figure.figsize"] @@ -703,7 +735,9 @@ class MaskedCost(Cost): __slots__ = "_data", "_mask", "_masked" - def __init__(self, args, data, verbose): + _mask: Optional[NDArray] + + def __init__(self, args: Sequence[str], data: NDArray, verbose: int): """For internal use.""" self._data = data self._mask = None @@ -720,7 +754,7 @@ def mask(self): return self._mask @mask.setter - def mask(self, mask): + def mask(self, mask: Optional[ArrayLike]): self._mask = None if mask is None else np.asarray(mask) self._update_cache() @@ -730,7 +764,7 @@ def data(self): return self._data @data.setter - def data(self, value): + def data(self, value: ArrayLike): self._data[...] = value self._update_cache() @@ -747,7 +781,7 @@ class UnbinnedCost(MaskedCost): __slots__ = "_model", "_log" - def __init__(self, data, model: _tp.Callable, verbose: int, log: bool): + def __init__(self, data, model: Model, verbose: int, log: bool): """For internal use.""" self._model = model self._log = log @@ -767,7 +801,7 @@ def _ndata(self): # unbinned likelihoods have infinite degrees of freedom return np.inf - def visualize(self, args: _ArrayLike, model_points: int = 0, nbins: int = 50): + def visualize(self, args: Sequence[float], model_points: int = 0, nbins: int = 50): """ Visualize data and model agreement (requires matplotlib). @@ -830,8 +864,8 @@ def scaled_pdf(self): def __init__( self, - data: _ArrayLike, - pdf: _tp.Callable, + data: ArrayLike, + pdf: Model, verbose: int = 0, log: bool = False, ): @@ -861,7 +895,7 @@ def __init__( """ super().__init__(data, pdf, verbose, log) - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: data = self._masked x = self._model(data, *args) x = _normalize_model_output(x) @@ -905,8 +939,8 @@ def scaled_pdf(self): def __init__( self, - data: _ArrayLike, - scaled_pdf: _tp.Callable, + data: ArrayLike, + scaled_pdf: Model, verbose: int = 0, log: bool = False, ): @@ -941,7 +975,7 @@ def __init__( """ super().__init__(data, scaled_pdf, verbose, log) - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: data = self._masked ns, x = self._model(data, *args) x = _normalize_model_output( @@ -961,6 +995,10 @@ class BinnedCost(MaskedCost): __slots__ = "_xe", "_ndim", "_bztrafo" + _xe: Union[NDArray, Tuple[NDArray, ...]] + _ndim: int + _bztrafo: Optional[BohmZechTransform] + n = MaskedCost.data @property @@ -968,14 +1006,24 @@ def xe(self): """Access bin edges.""" return self._xe - def __init__(self, args, n, xe, verbose, *updater): + def __init__( + self, + args: Sequence[str], + n: ArrayLike, + xe: Union[ArrayLike, Sequence[ArrayLike]], + verbose: int, + *updater, + ): """For internal use.""" - if not isinstance(xe, _tp.Iterable): + if not isinstance(xe, Iterable): raise ValueError("xe must be iterable") shape = _shape_from_xe(xe) self._ndim = len(shape) - self._xe = _norm(xe) if self._ndim == 1 else tuple(_norm(xei) for xei in xe) + if self._ndim == 1: + self._xe = _norm(xe) # type:ignore + else: + self._xe = tuple(_norm(xei) for xei in xe) n = _norm(n) is_weighted = n.ndim > self._ndim @@ -983,6 +1031,7 @@ def __init__(self, args, n, xe, verbose, *updater): if n.ndim != (self._ndim + int(is_weighted)): raise ValueError("n must either have same dimension as xe or one extra") + xei: NDArray for i, xei in enumerate([self._xe] if self._ndim == 1 else self._xe): if len(xei) != n.shape[i] + 1: raise ValueError( @@ -1003,7 +1052,7 @@ def _ndata(self): return np.prod(self._masked.shape[: self._ndim]) @abc.abstractmethod - def _pred(self, args: _tp.Sequence[float]) -> np.ndarray: + def _pred(self, args: Sequence[float]) -> Union[NDArray, Tuple[NDArray, NDArray]]: ... # pragma: no cover def _update_cache(self): @@ -1012,7 +1061,9 @@ def _update_cache(self): ma = _replace_none(self._mask, ...) self._bztrafo = BohmZechTransform(self._data[ma, 0], self._data[ma, 1]) - def prediction(self, args: _tp.Sequence[float]) -> np.ndarray: + def prediction( + self, args: Sequence[float] + ) -> Union[NDArray, Tuple[NDArray, NDArray]]: """ Return the bin expectation for the fitted model. @@ -1023,12 +1074,12 @@ def prediction(self, args: _tp.Sequence[float]) -> np.ndarray: Returns ------- - np.ndarray + NDArray Model prediction for each bin. """ return self._pred(args) - def visualize(self, args: _tp.Sequence[float]): + def visualize(self, args: Sequence[float]): """ Visualize data and model agreement (requires matplotlib). @@ -1062,7 +1113,7 @@ class BinnedCostWithModel(BinnedCost): :meta private: """ - __slots__ = "_xe_shape", "_model", "_model_arg" + __slots__ = "_xe_shape", "_model", "_model_xe" def __init__(self, n, xe, model, verbose): """For internal use.""" @@ -1072,15 +1123,15 @@ def __init__(self, n, xe, model, verbose): if self._ndim == 1: self._xe_shape = None - self._model_arg = _norm(self.xe) + self._model_xe = _norm(self.xe) else: self._xe_shape = tuple(len(xei) for xei in self.xe) - self._model_arg = np.row_stack( + self._model_xe = np.row_stack( [x.flatten() for x in np.meshgrid(*self.xe, indexing="ij")] ) - def _pred(self, args: _tp.Sequence[float]): - d = self._model(self._model_arg, *args) + def _pred(self, args: Sequence[float]) -> NDArray: + d = self._model(self._model_xe, *args) d = _normalize_model_output(d) if self._xe_shape is not None: d = d.reshape(self._xe_shape) @@ -1096,20 +1147,23 @@ class Template(BinnedCost): """ Binned cost function for a template fit with uncertainties on the template. - This cost function is for a mixture model. Samples originate from two or more - components and we are interested in estimating the yield that originates from each - component. In high-energy physics, one component is often a peaking signal over a - smooth background component. Templates are shape estimates for these components which - are obtained from Monte-Carlo simulation. Even if the Monte-Carlo simulation is exact, - the templates introduce some uncertainty since the Monte-Carlo simulation produces - only a finite sample of events that contribute to each template. This cost function - takes that additional uncertainty into account. - - There are several ways to approach this problem. Barlow and Beeston [1]_ found an - exact likelihood for this problem, with one nuisance parameter per component per bin. - Solving this likelihood is somewhat challenging though. The Barlow-Beeston likelihood - also does not handle the additional uncertainty in weighted templates unless the - weights per bin are all equal. + This cost function is for a mixture of components. Use this if the sample originate + from two or more components and you are interested in estimating the yield that + originates from one or more components. In high-energy physics, one component is often + a peaking signal over a smooth background component. A component can be described by a + parametric model or a template. + + A parametric model is accepted in form of a scaled cumulative density function, while + a template is a non-parametric shape estimate obtained by histogramming a Monte-Carlo + simulation. Even if the Monte-Carlo simulation is asymptotically correct, estimating + the shape from a finite simulation sample introduces some uncertainty. This cost + function takes that additional uncertainty into account. + + There are several ways to fit templates and take the sampling uncertainty into + account. Barlow and Beeston [1]_ found an exact likelihood for this problem, with one + nuisance parameter per component per bin. Solving this likelihood is somewhat + challenging though. The Barlow-Beeston likelihood also does not handle the additional + uncertainty in weighted templates unless the weights per bin are all equal. Other works [2]_ [3]_ [4]_ describe likelihoods that use only one nuisance parameter per bin, which is an approximation. Some marginalize over the nuisance parameters with @@ -1127,7 +1181,8 @@ class Template(BinnedCost): All methods implemented here have been generalized to work with both weighted data and weighted templates, under the assumption that the weights are independent of the data. This is not the case for sWeights, and the uncertaintes for results obtained with - sWeights will only be approximately correct [6]_. + sWeights will only be approximately correct [6]_. The methods have been further + generalized to allow fitting a mixture of parametric models and templates. .. [1] Barlow and Beeston, Comput.Phys.Commun. 77 (1993) 219-228 .. [2] Conway, PHYSTAT 2011 proceeding, https://doi.org/10.48550/arXiv.1103.0354 @@ -1137,14 +1192,14 @@ class Template(BinnedCost): .. [6] Langenbruch, Eur.Phys.J.C 82 (2022) 5, 393 """ - __slots__ = "_bbl_data", "_impl" + __slots__ = "_model_data", "_model_xe", "_xe_shape", "_impl" def __init__( self, - n: _ArrayLike, - xe: _ArrayLike, - templates: _tp.Sequence[_tp.Sequence], - name: _tp.Collection[str] = None, + n: ArrayLike, + xe: Union[ArrayLike, Sequence[ArrayLike]], + model_or_template: Collection[Union[Model, ArrayLike]], + name: Optional[Sequence[str]] = None, verbose: int = 0, method: str = "da", ): @@ -1161,14 +1216,17 @@ def __init__( Bin edge locations, must be len(n) + 1, where n is the number of bins. If the histogram has more than one axis, xe must be a collection of the bin edge locations along each axis. - templates : collection of array-like - Collection of arrays, which contain the histogram counts of each template. The - template histograms must use the same axes as the data histogram. If the - counts are represented by an array with dimension D+1, where D is the number - of histogram axes, then the last dimension must have two elements and is - interpreted as pairs of sum of weights and sum of weights squared. + model_or_template : collection of array-like or callable + Collection of models or arrays. An array represent the histogram counts of a + template. The template histograms must use the same axes as the data + histogram. If the counts are represented by an array with dimension D+1, where + D is the number of histogram axes, then the last dimension must have two + elements and is interpreted as pairs of sum of weights and sum of weights + squared. Callables must return the model cdf evaluated as xe. name : collection of str, optional - Optional name for the yield of each template. Must have length K. + Optional name for the yield of each template and the parameter of each model + (in order). Must have the same length as there are templates and model + parameters in templates_or_model. verbose : int, optional Verbosity level. 0: is no output (default). 1: print current args and negative log-likelihood value. @@ -1179,40 +1237,55 @@ def __init__( this parameter explicitly in code that has to be stable. For all methods except the "asy" method, the minimum value is chi-square distributed. """ - M = len(templates) + M = len(model_or_template) if M < 1: - raise ValueError("at least one template is required") - if name is None: - name = [f"x{i}" for i in range(M)] - else: - if len(name) != M: - raise ValueError("number of names must match number of templates") + raise ValueError("at least one template or model is required") shape = _shape_from_xe(xe) ndim = len(shape) - temp = [] - temp_var = [] - for ti in templates: - t = _norm(ti) - if t.ndim > ndim: - # template is weighted - if t.ndim != ndim + 1 or t.shape[:-1] != shape: - raise ValueError("shapes of n and templates do not match") - temp.append(t[..., 0]) - temp_var.append(t[..., 1]) + + npar = 0 + args = [] + self._model_data: List[ + Union[ + Tuple[NDArray, NDArray], + Tuple[Model, float], + ] + ] = [] + for i, t in enumerate(model_or_template): + if isinstance(t, Collection): + tt = _norm(t) + if tt.ndim > ndim: + # template is weighted + if tt.ndim != ndim + 1 or tt.shape[:-1] != shape: + raise ValueError("shapes of n and templates do not match") + t1 = tt[..., 0] + t2 = tt[..., 1] + else: + if tt.ndim != ndim or tt.shape != shape: + raise ValueError("shapes of n and templates do not match") + t1 = tt + t2 = tt.copy() + # normalize to unity + f = 1 / np.sum(t1) + t1 *= f + t2 *= f**2 + self._model_data.append((t1, t2)) + args.append(f"x{i}") else: - if t.ndim != ndim or t.shape != shape: - raise ValueError("shapes of n and templates do not match") - temp.append(t) - temp_var.append(t) - - nt = [] - nt_var = [] - for t, tv in zip(temp, temp_var): - f = 1 / np.sum(t) - nt.append(t * f) - nt_var.append(tv * f**2) - self._bbl_data = (nt, nt_var) + assert isinstance(t, Model) + par = describe(t)[1:] + npar = len(par) + self._model_data.append((t, npar)) + args += [f"x{i}_{x}" for x in par] + + if name is None: + name = args + else: + if len(args) != len(name): + raise ValueError( + "number of names must match number of templates and model parameters" + ) known_methods = { "jsc": template_chi2_jsc, @@ -1236,16 +1309,42 @@ def __init__( super().__init__(name, n, xe, verbose) - def _pred(self, args: _tp.Sequence[float]): - ntemp, ntemp_var = self._bbl_data - mu = 0 - mu_var = 0 - for a, nt, vnt in zip(args, ntemp, ntemp_var): - mu += a * nt - mu_var += a**2 * vnt + if self._ndim == 1: + self._xe_shape = None + self._model_xe = _norm(self.xe) + else: + self._xe_shape = tuple(len(xei) for xei in self.xe) + self._model_xe = np.row_stack( + [x.flatten() for x in np.meshgrid(*self.xe, indexing="ij")] + ) + + def _pred(self, args: Sequence[float]) -> Tuple[NDArray, NDArray]: + mu: NDArray = 0 # type:ignore + mu_var: NDArray = 0 # type:ignore + i = 0 + for t1, t2 in self._model_data: + if isinstance(t1, np.ndarray) and isinstance(t2, np.ndarray): + a = args[i] + mu += a * t1 + mu_var += a**2 * t2 + i += 1 + else: + assert isinstance(t2, int) + d = t1(self._model_xe, *args[i : i + t2]) + d = _normalize_model_output(d) + if self._xe_shape is not None: + d = d.reshape(self._xe_shape) + for i in range(self._ndim): + d = np.diff(d, axis=i) + # differences can come out negative due to round-off error in subtraction, + # we set negative values to zero + d[d < 0] = 0 + mu += d + mu_var += np.ones_like(mu) * 1e-300 + i += t2 return mu, mu_var - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: mu, mu_var = self._pred(args) ma = self.mask @@ -1261,10 +1360,10 @@ def _call(self, args): ma = mu > 0 return self._impl(n[ma], mu[ma], mu_var[ma]) - def _errordef(self): + def _errordef(self) -> float: return NEGATIVE_LOG_LIKELIHOOD if self._impl is template_nll_asy else CHISQUARE - def prediction(self, args: _tp.Sequence[float]): + def prediction(self, args: Sequence[float]) -> Tuple[NDArray, NDArray]: """ Return the fitted template and its standard deviation. @@ -1282,14 +1381,14 @@ def prediction(self, args: _tp.Sequence[float]): Returns ------- - y, yerr : np.ndarray, np.array + y, yerr : NDArray, np.array Template prediction and its standard deviation, based on the statistical uncertainty of the template only. """ mu, mu_var = self._pred(args) return mu, np.sqrt(mu_var) - def visualize(self, args: _ArrayLike): + def visualize(self, args: Sequence[float]): """ Visualize data and model agreement (requires matplotlib). @@ -1316,8 +1415,10 @@ def visualize(self, args: _ArrayLike): cx = cx[self.mask] plt.errorbar(cx, n, ne, fmt="ok") - mu, mu_err = self.prediction(args) - plt.stairs(mu + mu_err, xe, baseline=mu - mu_err, fill=True, color="C0") + mu, mu_err = self.prediction(args) # type: ignore + # need fill and line so that bins with mu_err=0 show up + for fill in (False, True): + plt.stairs(mu + mu_err, xe, baseline=mu - mu_err, fill=fill, color="C0") class BinnedNLL(BinnedCostWithModel): @@ -1341,7 +1442,11 @@ def cdf(self): return self._model def __init__( - self, n: _ArrayLike, xe: _ArrayLike, cdf: _tp.Callable, verbose: int = 0 + self, + n: ArrayLike, + xe: Union[ArrayLike, Sequence[ArrayLike]], + cdf: Model, + verbose: int = 0, ): """ Initialize cost function with data and model. @@ -1369,7 +1474,7 @@ def __init__( """ super().__init__(n, xe, cdf, verbose) - def _pred(self, args): + def _pred(self, args: Sequence[float]) -> NDArray: p = super()._pred(args) ma = self.mask if ma is not None: @@ -1377,7 +1482,7 @@ def _pred(self, args): scale = np.sum(self._masked[..., 0] if self._bztrafo else self._masked) return p * scale - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: mu = self._pred(args) ma = self.mask if ma is not None: @@ -1416,9 +1521,9 @@ def scaled_cdf(self): def __init__( self, - n: _ArrayLike, - xe: _ArrayLike, - scaled_cdf: _tp.Callable, + n: ArrayLike, + xe: Union[ArrayLike, Sequence[ArrayLike]], + scaled_cdf: Model, verbose: int = 0, ): """ @@ -1445,7 +1550,7 @@ def __init__( """ super().__init__(n, xe, scaled_cdf, verbose) - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: mu = self._pred(args) ma = self.mask if ma is not None: @@ -1454,6 +1559,7 @@ def _call(self, args): n, mu = self._bztrafo(mu) else: n = self._masked + # assert isinstance(n, np.ndarray) return poisson_chi2(n, mu) @@ -1467,6 +1573,8 @@ class LeastSquares(MaskedCost): __slots__ = "_loss", "_cost", "_model", "_ndim" + _loss: Union[str, LossFunction] + @property def x(self): """Get explanatory variables.""" @@ -1513,7 +1621,7 @@ def loss(self): return self._loss @loss.setter - def loss(self, loss: _tp.Union[str, _tp.Callable[[np.ndarray], np.ndarray]]): + def loss(self, loss: Union[str, LossFunction]): self._loss = loss if isinstance(loss, str): if loss == "linear": @@ -1522,19 +1630,20 @@ def loss(self, loss: _tp.Union[str, _tp.Callable[[np.ndarray], np.ndarray]]): self._cost = _soft_l1_cost else: raise ValueError("unknown loss type: " + loss) - else: - assert hasattr(loss, "__call__") + elif isinstance(loss, LossFunction): self._cost = lambda y, ye, ym: np.sum( loss(_z_squared(y, ye, ym)) # type:ignore ) + else: + raise ValueError("loss must be str or LossFunction") def __init__( self, - x: _ArrayLike, - y: _ArrayLike, - yerror: _ArrayLike, - model: _tp.Callable, - loss: _tp.Union[str, _tp.Callable] = "linear", + x: ArrayLike, + y: ArrayLike, + yerror: ArrayLike, + model: Model, + loss: Union[str, LossFunction] = "linear", verbose: int = 0, ): """ @@ -1590,7 +1699,7 @@ def __init__( super().__init__(describe(self._model)[1:], data, verbose) - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: x = self._masked.T[0] if self._ndim == 1 else self._masked.T[: self._ndim] y, yerror = self._masked.T[self._ndim :] ym = self._model(x, *args) @@ -1600,7 +1709,7 @@ def _call(self, args): def _ndata(self): return len(self._masked) - def visualize(self, args: _ArrayLike, model_points: int = 0): + def visualize(self, args: ArrayLike, model_points: int = 0): """ Visualize data and model agreement (requires matplotlib). @@ -1664,9 +1773,9 @@ class NormalConstraint(Cost): def __init__( self, - args: _tp.Union[str, _tp.Iterable[str]], - value: _ArrayLike, - error: _ArrayLike, + args: Union[str, Iterable[str]], + value: ArrayLike, + error: ArrayLike, ): """ Initialize the normal constraint with expected value(s) and error(s). @@ -1712,7 +1821,7 @@ def value(self): def value(self, value): self._value[:] = value - def _call(self, args): + def _call(self, args: Sequence[float]) -> float: delta = self._value - args if self._covinv.ndim < 2: return np.sum(delta**2 * self._covinv) @@ -1721,7 +1830,7 @@ def _call(self, args): def _ndata(self): return len(self._value) - def visualize(self, args: _ArrayLike): + def visualize(self, args: ArrayLike): """ Visualize data and model agreement (requires matplotlib). @@ -1760,7 +1869,7 @@ def visualize(self, args: _ArrayLike): plt.ylim(-n + 0.5, 0.5) -def _norm(value: _ArrayLike[float]) -> np.ndarray: +def _norm(value: ArrayLike) -> NDArray: value = np.atleast_1d(value) dtype = value.dtype if dtype.kind != "f": @@ -1785,7 +1894,7 @@ def _normalize_model_output(x, msg="Model should return numpy array"): def _shape_from_xe(xe): - if isinstance(xe[0], _tp.Iterable): + if isinstance(xe[0], Iterable): return tuple(len(xei) - 1 for xei in xe) return (len(xe) - 1,) @@ -1797,7 +1906,7 @@ def _shape_from_xe(xe): } -def __getattr__(name: str) -> _tp.Any: +def __getattr__(name: str) -> Any: if name in _deprecated_content: new_name, obj = _deprecated_content[name] warnings.warn( diff --git a/src/iminuit/typing.py b/src/iminuit/typing.py index 9b9dbedf..e503f3ba 100644 --- a/src/iminuit/typing.py +++ b/src/iminuit/typing.py @@ -1,12 +1,41 @@ -""" -Types for iminuit. +"""Types for iminuit. These are used by mypy and similar tools. """ -import typing as _tp +from typing import ( + Protocol, + Optional, + Collection, + List, + Union, + runtime_checkable, +) +from numpy.typing import NDArray +import numpy as np + +# LossFunction = Callable[[np.ndarray], np.ndarray] -UserBound = _tp.Optional[_tp.Collection[_tp.Optional[float]]] +# Used by Minuit +UserBound = Optional[Collection[Optional[float]]] # Key for ValueView, ErrorView, etc. -Key = _tp.Union[int, str, slice, _tp.List[_tp.Union[int, str]]] +Key = Union[int, str, slice, List[Union[int, str]]] + + +@runtime_checkable +class Model(Protocol): + """Type for user-defined model.""" + + def __call__(self, x: np.ndarray, *args: float) -> np.ndarray: + """Evaluate model at locations x and return results as an array.""" + ... # pragma: no cover + + +@runtime_checkable +class LossFunction(Protocol): + """Type for user-defined loss function for LeastSquares clas.""" + + def __call__(self, z: NDArray) -> NDArray: + """Evaluate loss function on values.""" + ... # pragma: no cover diff --git a/tests/test_cost.py b/tests/test_cost.py index cec6eef1..2efb3845 100644 --- a/tests/test_cost.py +++ b/tests/test_cost.py @@ -16,9 +16,9 @@ _soft_l1_loss, PerformanceWarning, ) +from iminuit.util import describe from typing import Sequence import pickle -from sys import version_info as pyver try: # pytest.importorskip does not work for scipy.stats; @@ -48,6 +48,14 @@ def norm_pdf(x, mu, sigma): return np.exp(norm_logpdf(x, mu, sigma)) +def norm_cdf(x, mu, sigma): + from math import erf + + z = (x - mu) / (np.sqrt(2) * sigma) + + return (1 + np.vectorize(erf)(z)) * 0.5 + + def mvnorm(mux, muy, sx, sy, rho): C = np.empty((2, 2)) C[0, 0] = sx**2 @@ -87,11 +95,11 @@ def pdf(x, mu, sigma): def cdf(x, mu, sigma): - return norm(mu, sigma).cdf(x) + return norm_cdf(x, mu, sigma) def scaled_cdf(x, n, mu, sigma): - return n * norm(mu, sigma).cdf(x) + return n * norm_cdf(x, mu, sigma) def line(x, a, b): @@ -110,6 +118,12 @@ def test_norm_pdf(): assert_allclose(norm_pdf(x, 3, 2), norm.pdf(x, 3, 2)) +@pytest.mark.skipif(not scipy_available, reason="scipy.stats is needed") +def test_norm_cdf(): + x = np.linspace(-3, 3) + assert_allclose(norm_cdf(x, 3, 2), norm.cdf(x, 3, 2)) + + def test_Constant(): c = Constant(2.5) assert c.value == 2.5 @@ -346,7 +360,6 @@ def test_ExtendedUnbinnedNLL_pickle(): assert_equal(c.data, c2.data) -@pytest.mark.skipif(not scipy_available, reason="scipy.stats is needed") @pytest.mark.parametrize("verbose", (0, 1)) def test_BinnedNLL(binned, verbose): mle, nx, xe = binned @@ -528,7 +541,6 @@ def test_BinnedNLL_pickle(): assert_equal(c.data, c2.data) -@pytest.mark.skipif(not scipy_available, reason="scipy.stats is needed") @pytest.mark.parametrize("verbose", (0, 1)) def test_ExtendedBinnedNLL(binned, verbose): mle, nx, xe = binned @@ -1221,7 +1233,36 @@ def test_Template_pickle(): assert_equal(c.data, c2.data) -@pytest.mark.skipif(pyver < (3, 7), reason="module getattr requires Python-3.7+") +def test_Template_with_model_template_mix(): + n = np.array([3, 2, 3]) + xe = np.array([0, 1, 2, 3]) + t = np.array([0.1, 0, 1]) + + c = Template(n, xe, (t, scaled_cdf)) + assert describe(c) == ["x0", "x1_n", "x1_mu", "x1_sigma"] + + m = Minuit(c, 1, 1, 0.5, 1) + m.limits["x0", "x1_n"] = (0, None) + m.limits["x1_sigma"] = (0.1, None) + m.limits["x1_mu"] = (xe[0], xe[-1]) + m.migrad() + assert m.valid + + +def test_Template_with_models(): + n = np.array([3, 2, 3]) + xe = np.array([0, 1, 2, 3]) + + c = Template(n, xe, (lambda x, n, mu: n * expon_cdf(x, mu), scaled_cdf)) + assert describe(c) == ["x0_n", "x0_mu", "x1_n", "x1_mu", "x1_sigma"] + + m = Minuit(c, 1, 0.5, 1, 0.5, 1) + m.limits["x0_n", "x1_n", "x1_sigma"] = (0, None) + m.limits["x0_mu", "x1_mu"] = (xe[0], xe[-1]) + m.migrad() + assert m.valid + + def test_deprecated(): from iminuit import cost