{ "cells": [ { "cell_type": "markdown", "id": "naked-recruitment", "metadata": {}, "source": [ "# Fit PDF with conditional variable\n", "\n", "In this example, we show an unusual fit where the total sample is not drawn form a single probability distribution, but each individual sample $x$ is drawn from a different distribution, whose parameters are determined by a conditional variable $y$.\n", "\n", "In our example, we are drawing samples $x$ from varying Gaussian distributions. The location of each Gaussian is a function of the conditional variable $y$, but all share the same width parameter $\\sigma$. We fit the shared parameter $\\sigma$, but also the parameters $a$ and $b$ which determine how the location of each gaussian depends on $y$, assuming a line function $\\mu = a + b y$.\n", "\n", "This tutorial reproduces a [corresponding one from RooFit](https://root.cern.ch/doc/master/rf303__conditional_8C.html)." ] }, { "cell_type": "code", "execution_count": 1, "id": "technological-economy", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "iminuit version 2.19.0\n" ] } ], "source": [ "import iminuit\n", "from iminuit.cost import UnbinnedNLL\n", "from iminuit import Minuit\n", "import numpy as np\n", "import numba as nb\n", "import boost_histogram as bh\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import norm\n", "from numba_stats import norm as norm_nb\n", "print(\"iminuit version\", iminuit.__version__)" ] }, { "cell_type": "code", "execution_count": 2, "id": "wicked-animal", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(1)\n", "\n", "# conditional variable: each sample is paired with a random y parameter\n", "y = rng.normal(0, 10, size=10000)\n", "y = y[np.abs(y) < 10] # truncate at 10\n", "\n", "# location of each gaussian is a function of y\n", "def mu(y, a, b):\n", " return a + b * y\n", "\n", "# draw samples from Gaussians whose locations depend on y\n", "truth = {\"a\": 0, \"b\": 0.5, \"sigma\": 1.0}\n", "x = rng.normal(mu(y, truth[\"a\"], truth[\"b\"]), truth[\"sigma\"])" ] }, { "cell_type": "markdown", "id": "removable-forward", "metadata": {}, "source": [ "The distribution in $x$ is more broad than the usual Gaussian because it is a convolution of many Gaussian distributions with varying means. We can visualise this by binning the data in $x$ and $y$." ] }, { "cell_type": "code", "execution_count": 3, "id": "subjective-sleep", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbSklEQVR4nO3deVxU9f4/8BeoLCMCAsKwgxtO7mhw0VIMEpevN29muXBzQagumkKaWubWTSp3zbTN7aLX0sxKTUMTl8ANJYUQA9FxAQxJEEFAOL8//HEuAzMswzAbr+fjcR4Pzuds78MIvP2c9/l8TARBEEBERERkpEx1HQARERFRc2KyQ0REREaNyQ4REREZNSY7REREZNSY7BAREZFRY7JDRERERo3JDhERERm11roOQB9UVlbizp07aNeuHUxMTHQdDhERETWAIAh48OABXFxcYGqquv+GyQ6AO3fuwN3dXddhEBERkRpu3rwJNzc3lduZ7ABo164dgCffLGtrax1HQ0RERA1RWFgId3d38e+4Kkx2APHRlbW1NZMdIiIiA1NfCQoLlImIiMioMdkhIiIio8Zkh4iIiIwakx0iIiIyakx2iIiIyKgx2SEiIiKjxmSHiIiIjBqTHSIiIjJqTHaIiIjIqDHZISIiIqPGZIeIiIiMGpMdIiKiZnT9+nWYmJggOTlZ16E0iZeXF9asWaPrMNTCZIeIiKgZubu7Izs7Gz169GjwMYsXL0afPn2aL6gWhskOERFRM2rVqhWkUilat26t9WuXlZVp/Zr6iMkOERkVuVyOCxcuKCxyuVzXYZGR2L59O+zt7VFaWqrQPnr0aPzzn/9UekzNx1jx8fEwMTHB0aNH0b9/f0gkEgwYMADp6ekAgK1bt2LJkiX47bffYGJiAhMTE2zduhUAcP/+fUybNg0dOnSAtbU1nnvuOfz222/itap6hL788kt4e3vDwsICn3/+OVxcXFBZWakQ1wsvvICpU6cCADIzM/HCCy/AyckJVlZWePrpp3HkyBFNfMv0ApMdIjIacrkcMpkM/fr1U1hkMhkTHtKIsWPHoqKiAj/88IPYdvfuXRw4cEBMHBrq3XffxcqVK3H+/Hm0bt1aPP6VV17BW2+9he7duyM7OxvZ2dl45ZVXxOvfvXsXP/30E5KSkuDr64ugoCDk5+eL583IyMC3336LvXv3Ijk5GWPHjsW9e/dw7NgxcZ/8/HwcOnQIEydOBAAUFRVhxIgROHr0KC5evIhhw4Zh1KhRRvNzw2SHiIxGXl4eiouLERsbi6SkJCQlJSE2NhbFxcXIy8vTdXhkBCwtLTFhwgRs2bJFbIuNjYWHhwcCAwMbda4PPvgAgwcPxlNPPYV58+YhISEBjx49gqWlJaysrNC6dWtIpVJIpVJYWlri1KlTOHv2LHbv3o3+/fujS5cuWLFiBWxtbbFnzx7xvGVlZdi+fTv69u2LXr16oX379hg+fDh27twp7rNnzx44ODhgyJAhAIDevXvjtddeQ48ePdClSxe8//776NSpk0JSZ8h0muycOHECo0aNgouLC0xMTLBv3z6F7VXddzWX5cuXi/t4eXnV2v7hhx9q+U6ISJ/IZDL4+vrC19cXMplM1+GQkQkPD8fPP/+M27dvA3jy2Gny5MkwMTFp1Hl69eolfu3s7AzgSS+RKr/99huKiopgb28PKysrccnKykJmZqa4n6enJzp06KBw7MSJE/Htt9+Kj9927NiBcePGwdT0SRpQVFSE2bNnQyaTwdbWFlZWVkhLSzOanh3tV0tV8/DhQ/Tu3RtTp07Fiy++WGt7dna2wvpPP/2EsLAwjBkzRqF96dKlCA8PF9fbtWvXPAETEVGL17dvX/Tu3Rvbt2/H0KFDkZqaigMHDjT6PG3atBG/rkqUatbVVFdUVARnZ2fEx8fX2mZrayt+3bZt21rbR40aBUEQcODAATz99NM4efIkVq9eLW6fPXs24uLisGLFCnTu3BmWlpZ46aWXjKbAWafJzvDhwzF8+HCV26VSqcL6999/jyFDhqBjx44K7e3atau1LxERUXOZNm0a1qxZg9u3byM4OBju7u4aPb+ZmRkqKioU2nx9fZGTk4PWrVvDy8urUeezsLDAiy++iB07diAjIwM+Pj7w9fUVt//666+YPHky/vGPfwB4klhdv369qbehNwymZic3NxcHDhxAWFhYrW0ffvgh7O3t0bdvXyxfvhyPHz+u81ylpaUoLCxUWIiIiBpqwoQJuHXrFr744otGFyY3hJeXF7KyspCcnIy8vDyUlpYiODgYAQEBGD16NH7++Wdcv34dCQkJePfdd3H+/Pl6zzlx4kQcOHAAmzdvFguTq3Tp0kUsaP7tt98wYcKEOnuZDI3BJDvbtm1Du3btaj3uevPNN7Fr1y4cO3YMr732GpYtW4a33367znPFxMTAxsZGXDSdkRMRkXGzsbHBmDFjYGVlhdGjR2v8/GPGjMGwYcMwZMgQdOjQAf/9739hYmKCgwcPYtCgQZgyZQq6du2KcePG4caNG3Bycqr3nM899xzs7OyQnp6OCRMmKGxbtWoV2rdvjwEDBmDUqFEICQlR6PkxdCaCIAi6DgJ48rzyu+++U/mPplu3bnj++eexfv36Os+zefNmvPbaaygqKoK5ubnSfUpLSxXGSCgsLIS7uzsKCgpgbW2t9j0QkW5duHAB/fr1E1/JVdVGpAlBQUHo3r071q1bp+tQWqzCwkLY2NjU+/dbpzU7DXXy5Emkp6fj66+/rndff39/PH78GNevX4ePj4/SfczNzVUmQkRERHX566+/EB8fj/j4eHz66ae6DocawCCSna+++gr9+vVD79696903OTkZpqamcHR01EJkRETU0vTt2xd//fUXPvroI5X/qSb9otNkp6ioCBkZGeJ6VTGWnZ0dPDw8ADzpotq9ezdWrlxZ6/jExEScOXMGQ4YMQbt27ZCYmIioqCiEhoaiffv2WrsPIiJqOYzpLaWWQqfJzvnz58XRGwEgOjoaADBp0iRxHpBdu3ZBEASMHz++1vHm5ubYtWsXFi9ejNLSUnh7eyMqKko8DxEREZFOk53AwEDUVx8dERGBiIgIpdt8fX1x+vTp5giNiIiIjITBvHpOREREpA4mO0RERGTUmOwQERGRUWOyQ0REREbNIMbZISIiaozLtwqafI6ebjaN2j8wMBDHjx8HAFy8eBF9+vRpcgzGrKysDF27dsWePXvQv3//Zr0We3aIqEnkcjkuXLigsMjlcl2H1Wg178MQ74F0Lzw8HNnZ2ejRo4fKfR49eoTJkyejZ8+eaN26tcppkuLj4+Hr6wtzc3N07txZHJKlLpcuXcKzzz4LCwsLuLu74+OPP1bzTmpLTEzEc889h7Zt28La2hqDBg1CSUlJncds2LABXl5esLCwgL+/P86ePStuMzMzw+zZszF37lyNxagKe3aISG1yuRwymQzFxcUK7RKJBGlpaeLgoPpO2X0Y2j2QfpBIJJBKpXXuU1FRAUtLS7z55pv49ttvle6TlZWFkSNH4vXXX8eOHTtw9OhRTJs2Dc7OzggJCVF6TGFhIYYOHYrg4GBs2rQJly9fxtSpU2Fra6tyCJeGSkxMxLBhwzB//nysX78erVu3xm+//QZTU9V9Jl9//TWio6OxadMm+Pv7Y82aNQgJCUF6ero4y8HEiRPx1ltvITU1Fd27d29SjHVhzw4RqS0vLw/FxcWIjY1FUlISkpKSEBsbi+LiYuTl5ek6vAareR+GeA9kONq2bYuNGzciPDxcZWK0adMmeHt7Y+XKlZDJZJg+fTpeeuklrF69WuV5d+zYgbKyMmzevBndu3fHuHHj8Oabb2LVqlVNjjkqKgpvvvkm5s2bh+7du8PHxwcvv/xynfNMrlq1CuHh4ZgyZQqeeuopbNq0CRKJBJs3bxb3ad++PQYOHIhdu3Y1Oca6MNkhoiaTyWTw9fWFr68vZDKZrsNRW9V9GPI9kHFITExEcHCwQltISAgSExPrPGbQoEEwMzNTOCY9PR1//fWX2rHcvXsXZ86cgaOjIwYMGAAnJycMHjwYp06dUnlMWVkZkpKSFO7B1NQUwcHBte7Bz88PJ0+eVDu+huBjLCLSC3K5vFZPioODAx8jUYuUk5MDJycnhTYnJycUFhaipKQElpaWSo/x9vaudUzVNnXnjLx27RoAYPHixVixYgX69OmD7du3IygoCCkpKejSpUutY/Ly8lBRUaH0Hq5cuaLQ5uLighs3bqgVW0OxZ4eIdK6qZqZfv34Ki0wmY6EwGbTu3bvDysoKVlZWGD58uK7DqdeyZcvEeK2srCCXy1FZWQkAeO211zBlyhT07dsXq1evho+Pj8IjKXVZWlrWqvvTNPbsEJHOVa+ZqXqElJaWhtDQUOTl5bF3hwzWwYMHUV5eDgBKe2NUkUqlyM3NVWjLzc2FtbW1yvOoOqZqW0O8/vrrePnll8V1FxcXVFRUAACeeuophX3r+s+Ig4MDWrVqpTSemrHk5+ejQ4cODYpPXUx2iEhvVNXMEBkLT09PtY4LCAjAwYMHFdri4uIQEBBQ5zHvvvsuysvL0aZNG/EYHx+fBj/CsrOzg52dnUKbl5cXXFxckJ6ertB+9epVlb1VZmZm6NevH44ePSq+Wl9ZWYmjR49i+vTpCvumpKSgb9++DYpPXXyMRUREpGW///47kpOTkZ+fj4KCAiQnJyM5OVnc/vrrr+PatWt4++23ceXKFXz66af45ptvEBUVJe7zySefICgoSFyfMGECzMzMEBYWhtTUVHz99ddYu3YtoqOjmxSriYkJ5syZg3Xr1mHPnj3IyMjAe++9hytXriAsLEzcLygoCJ988om4Hh0djS+++ALbtm1DWloa3njjDTx8+BBTpkxROP/JkycxdOjQJsVYH/bsEBGR0Wns6MfaNmLECIWi3KqeDUEQAADe3t44cOAAoqKisHbtWri5ueHLL79UGGMnLy8PmZmZ4rqNjQ1+/vlnREZGol+/fnBwcMDChQsVxtiJj4/HkCFDkJWVBS8vrwbHO2vWLDx69AhRUVHIz89H7969ERcXh06dOon7ZGZmKrxk8Morr+DPP//EwoULkZOTgz59+uDQoUMKRcuJiYkoKCjASy+91OBY1MFkh4iISMuuX79e7z6BgYG4ePGiyu2LFy/G4sWLFdp69epV52vcWVlZ6Ny5M1xdXRsaqmjevHmYN2+eyu3K7mn69Om1HltVt2bNGsyZM6dR9Uzq4GMsIiIiDfn0009hZWWFy5cv6zoUpQ4ePIhly5aJNT26VFZWhp49eyo8mmsu7NkhIiLSgB07dohzRenrG4S7d+/WdQgiMzMzLFiwQCvXYrJDRESkAeo8GiLt4GMsIiIiMmpMdoiIiMio8TEWEVEz4pxfRLrHZIeIqJlUzflVc94fiUSCtLQ0JjxEWsLHWEREzaT6nF9JSUlISkpCbGwsiouLa/X2EFHzYc8OEVEz45xfRLrFZIeIiIzPHdUjDzeYS+MmpwwMDMTx48cBABcvXkSfPn2aHoMB+9vf/oY5c+ZgzJgxug6Fj7GIiIg0JTw8HNnZ2ejRo0ed+wmCgBUrVqBr164wNzeHq6srPvjggzqPyc/Px8SJE2FtbQ1bW1uEhYWhqKioyTHHxMTg6aefRrt27eDo6IjRo0fXmuFcmd27d6Nbt26wsLBAz549a83SvmDBAsybNw+VlZVNjrGpmOwQERFpiEQigVQqRevWdT84mTlzJr788kusWLECV65cwQ8//AA/P786j5k4cSJSU1MRFxeH/fv348SJEwqTfKrr+PHjiIyMxOnTpxEXF4fy8nIMHToUDx8+VHlMQkICxo8fj7CwMFy8eBGjR4/G6NGjkZKSIu4zfPhwPHjwAD/99FOTY2wqPsYiIiLSorS0NGzcuBEpKSnw8fEB8GSW8/qOOXToEM6dO4f+/fsDANavX48RI0ZgxYoVcHFxUTueQ4cOKaxv3boVjo6OSEpKwqBBg5Qes3btWgwbNgxz5swBALz//vuIi4vDJ598gk2bNgEAWrVqhREjRmDXrl0YOXKk2vFpAnt2iEgn5HI5Lly4gAsXLiAtLa3Rx1Qtcrm8mSMl0qwff/wRHTt2xP79++Ht7Q0vLy9MmzYN+fn5Ko9JTEyEra2tmOgAQHBwMExNTXHmzBmNxldQUAAAsLOzqzOe4OBghbaQkBAkJiYqtPn5+dU5C7u2sGeHiLRO2fgzEokEDg4OjTqm6jiOWUOG5Nq1a7hx4wZ2796N7du3o6KiAlFRUXjppZfwyy+/KD0mJycHjo6OCm2tW7eGnZ0dcnJyNBZbZWUlZs2ahYEDB9ZZd5STkwMnJyeFNicnp1qxuLi44ObNm6isrISpqe76V9izQ0Rap2z8mfoSFo5ZQ4aoe/fusLKygpWVFYYPHw7gSUJRWlqK7du349lnn0VgYCC++uorHDt2rEGFweo6efKkGIuVlRV27NhRa5/IyEikpKRg165dGrmmpaWleL+6xJ4dItIZdcaf4Zg1ZEgOHjyI8vJyAE/+8AOAs7MzWrduja5du4r7yWQyAE96MKvqeKqTSqW4e/euQtvjx4+Rn58PqVTaoFj69++P5ORkcb1mz8z06dPFwmc3N7c6zyWVSpGbm6vQlpubWyuW/Px8tG3bVrx3XWGyQ9RC1JyjSZvzM9W8dkNrdDSp+jXVvXddfg/JMHl6etZqGzhwIB4/fozMzEx06tQJAHD16lWV+wNAQEAA7t+/j6SkJPTr1w8A8Msvv6CyshL+/v4NisXS0hKdO3eu1S4IAmbMmIHvvvsO8fHx9RZLV8Vz9OhRzJo1S2yLi4tDQECAwn4pKSno27dx4xU1ByY7RC2AqhoZbdS61FVrU1eNjqY4ODhAIpEgNDRU4dqNvXddfg/JuAQHB8PX1xdTp07FmjVrUFlZicjISDz//PNib8/Zs2fx6quv4ujRo3B1dYVMJsOwYcMQHh6OTZs2oby8HNOnT8e4ceOa9CYW8OTR1c6dO/H999+jXbt2Yt2NjY2N2CPz6quvwtXVFTExMQCevDo/ePBgrFy5EiNHjsSuXbtw/vx5fP755wrnPnnyJIYOHdqk+DSByQ5RC1C93kUmkyEtLQ2hoaHIy8tr9j/UNa9dRVu9Ih4eHkhLSxN7ZNS9d11+D0kNjRz9WJtMTU3x448/YsaMGRg0aBDatm2L4cOHY+XKleI+xcXFSE9PFx+BAcCOHTswffp0BAUFwdTUFGPGjMG6desUzm1iYoItW7Zg8uTJDY5n48aNAJ6MAF1d9fPI5XKFAuMBAwZg586dWLBgAd555x106dIF+/btUyhqvn37NhISEhAbG9vgWJoLkx2iFkSX9S66vLaHh4fGEhLWDJEmuLi44Ntvv1W5PTAwEIIgKLTZ2dlh586dKo/JyspC69atMXDgwEbFUvM6ysTHx9dqGzt2LMaOHavymHXr1mHy5Mn11v9oA5MdIjJ4VfU4uq4FAnRbC8UaIt379NNP8eWXXyIxMRE9e/bU6rUPHjyIiIgIdOnSRavXVcXR0RHR0dG6DgOAjpOdEydOYPny5UhKSkJ2dja+++47jB49Wtw+efJkbNu2TeGYkJAQhdEe8/PzMWPGDPz4449it97atWthZWWlrdsgIh1RVY+jq1qgquvrqhaKNUS6tWPHDpSUlACATj6DyMhIrV+zLm+99ZauQxDpNNl5+PAhevfujalTp+LFF19Uus+wYcOwZcsWcd3c3Fxh+8SJE5GdnS3O5zFlyhRERETU2dVHRMahZj0OoLtaIED9eiB1sIZI/7i6uuo6BFJBp8nO8OHDxUGWVDE3N1c5hkBzzhVCRIZBk/U4hnTtKqwhIqqf3o+gHB8fD0dHR/j4+OCNN97AvXv3xG3qzhVSWlqKwsJChYWINCstLa1R814RETUXvS5QHjZsGF588UV4e3sjMzMT77zzDoYPH47ExES0atVK7blCYmJisGTJkuYOn6hF0mUdDRGRMnqd7IwbN078umfPnujVqxc6deqE+Ph4BAUFqX3e+fPnK1SIFxYWwt3dvUmxEhmi+npd1Kl/0WUdDRGRMnqd7NTUsWNHODg4ICMjA0FBQWrPFWJubl6r0JmoJVH1JlFN6r7dow+1LEREVfS+Zqe6W7du4d69e3B2dgagOFdIlcbOFULUElX1vlTNHq5s4YziRI0TGBgIExMTmJiYKEy42RL87W9/q3OQRF3Tac9OUVERMjIyxPWsrCwkJyfDzs4OdnZ2WLJkCcaMGQOpVIrMzEy8/fbb6Ny5M0JCQgCgWecKITJ27H0hY5Z6L7XJ5+hu373Rx4SHh2Pp0qV11qjFx8dj9erVOHv2LAoLC9GlSxfMmTMHEydOrPPccrkcb7zxBo4dOwYrKytMmjQJMTExaN264X/KY2JisHfvXly5cgWWlpYYMGAAPvroI6UzrVe3e/duvPfee7h+/Tq6dOmCjz76CCNGjBC3L1iwAFFRUfjHP/6hMK2EvtBpROfPn0ffvn3FGVGjo6PRt29fLFy4EK1atcKlS5fw97//HV27dkVYWBj69euHkydPKjyC2rFjB7p164agoCCMGDECzzzzTK2JyIiItIlvorVcEokEUqm0zgQkISEBvXr1wrfffotLly5hypQpePXVV7F//36Vx1RUVGDkyJEoKytDQkICtm3bhq1bt2LhwoWNiu/48eOIjIzE6dOnxfHphg4diocPH9YZ7/jx4xEWFoaLFy9i9OjRGD16NFJSUsR9hg8fjgcPHuCnn35qVDzaotOeHWVzf1R3+PDhes9R31whRETawjfRqCHeeecdhfWZM2fi559/xt69e/F///d/So/5+eef8fvvv+PIkSNwcnJCnz598P7772Pu3LlYvHgxzMzMGnTt6jMQAMDWrVvh6OiIpKQkDBo0SOkxa9euxbBhwzBnzhwAwPvvv4+4uDh88skn2LRpEwCgVatWGDFiBHbt2oWRI0c2KBZt0r++JiIiA6WsForTN1BDFBQUwM7OTuX2qrm2nJycxLaQkBAUFhYiNVX9R3YFBQUAUO+1g4ODFdpCQkKQmJio0Obn54eTJ0+qHUtzMqi3sYiI9B1roaixvvnmG5w7dw6fffaZyn1ycnIUEh0A4npd48rVpbKyErNmzcLAgQPRo0ePRl+75nVdXFxw8+ZNVFZW6l3djn5FQ0REZES6d+8OKysrWFlZKZ0e6dixY5gyZQq++OILdO/e+ILopoiMjERKSgp27dqlkfNZWlqisrISpaWlGjmfJrFnh4iIqJkcPHgQ5eXlAJ4kA9UdP34co0aNwurVq/Hqq6/WeR6pVIqzZ88qtOXm5orbGmv69OnYv38/Tpw4ATc3t3qvXXWt6teued38/Hy0bdu21n3qA/bsEBERNRNPT0907twZnTt3VpgVPT4+HiNHjsRHH32EiIiIes8TEBCAy5cvKwykGxcXB2trazz11FMNjkcQBEyfPh3fffcdfvnlF3h7ezfo2kePHlVoi4uLQ0BAgEJbSkqK+Ha1vmGyQ0REpEXHjh3DyJEj8eabb2LMmDHIyclBTk4O8vPzxX2+++47dOvWTVwfOnQonnrqKfzzn//Eb7/9hsOHD2PBggWIjIxs1IwAkZGRiI2Nxc6dO9GuXTvx2iUlJeI+r776KubPny+uz5w5E4cOHcLKlStx5coVLF68GOfPn8f06dMVzn3y5EkMHTpUnW9Js+NjLCIiFaqPk8MxcwyLOgMCasu2bdtQXFyMmJgYxMTEiO2DBw9GfHw8gCdvSaWnp4vbWrVqhf379+ONN95AQEAA2rZti0mTJmHp0qXiPtevX4e3tzeOHTuGwMBApdfeuHEjANTavmXLFkyePBnAk8ELqxcYDxgwADt37sSCBQvwzjvvoEuXLti3b59CUfPt27eRkJCA2NhYdb4lzY7JDhFRDarmDuOYOaQJW7duxdatW+vcZ/LkyWLyUcXT0xMHDx5UeUxWVhZsbW3Ru3dvlfvUNbZdlaqEq7qxY8di7NixKo9Zt24dJk+eXG/9j64w2SEivVbVo6LNnhVlM7cDnL2d6vfpp5/iyy+/FMfF0aaDBw/inXfeQfv27bV6XQBwdHREdHS01q/bUEx2iEgv6Xo0Yo6XQ421Y8cOsfZFF/92li9frvVrVnnrrbd0du2GYLJDRHpJWe9KS+pZkcvlLfbeDVX1t61IvzDZISK91VJ7V+RyOWQyGYqLi8U2iUTCqSeI1MRXz4mI9ExeXh6Ki4sRGxuLpKQkxMbGori4uFYNERE1DHt2iIj0lEwmg6+vr67DIDJ4THaIqFGq15Jw7Bntaq5xf1gfRMaOyQ4RNZiqWhKOPdO8mnPcH9YHUUvAZIeIGqx6LYlMJgPAXgBtaM5xf2p+pmlpaQgNDUVeXh4/VzIaTHaIqNFYS6J9zf1mGj/TpgsMDMTx48cBABcvXkSfPn10G1AzmTdvHh4+fIj169frOpQG49tYRETUaHK5HBcuXFBY5HK5rsMSlaSkNnlRR3h4OLKzs8V5o3777TeMHz8e7u7usLS0hEwmw9q1a+s9T35+PiZOnAhra2vY2toiLCwMRUVFjY7n/v37iIyMhLOzM8zNzdG1a9c6p5wAgEuXLuHZZ5+FhYUF3N3d8fHHHytsnz17NrZt24Zr1641Oh5dYc8OERE1irI6H4C1PsCT74FUKhXXk5KS4OjoiNjYWLi7uyMhIQERERFo1apVrVnDq5s4cSKys7MRFxeH8vJyTJkyBREREdi5c2eDYykrK8Pzzz8PR0dH7NmzB66urrhx4wZsbW1VHlNYWIihQ4ciODgYmzZtwuXLlzF16lTY2toiIiICwJPHpyEhIdi4caNOR21uDCY7RETUKMpqt1jro9zUqVMV1jt27IjExETs3btXZbKTlpaGQ4cO4dy5c+jfvz8AYP369RgxYgRWrFgBFxeXBl178+bNyM/PR0JCAtq0aQMA8PLyqvOYHTt2oKysDJs3b4aZmRm6d++O5ORkrFq1Skx2AGDUqFF49913DSbZ4WMsIiJSS1Wdj6+vr5j0UP0KCgpgZ2encntiYiJsbW3FRAcAgoODYWpqijNnzjT4Oj/88AMCAgIQGRkJJycn9OjRA8uWLUNFRUWd1x40aBDMzMzEtpCQEKSnp+Ovv/4S2/z8/HDr1i1cv369wfHoEnt2iIh0oPo4OXyjreVISEjA119/jQMHDqjcJycnB46OjgptrVu3hp2dHXJychp8rWvXruGXX37BxIkTcfDgQWRkZOBf//oXysvLsWjRIpXX9vb2VmhzcnISt1XNqF7Vu3Tjxo16e4v0AZMdIiItUjWbe0uvdWkJUlJS8MILL2DRokUYOnRos1+vsrISjo6O+Pzzz9GqVSv069cPt2/fxvLly1UmOw1laWkJALXqtvQVkx0iIi2qOWZOVa3LyZMnFepfyLj8/vvvCAoKQkREBBYsWFDnvlKpFHfv3lVoe/z4MfLz8xWKn+vj7OyMNm3aoFWrVmKbTCZDTk4OysrKFB5VVb92bm6uQlvVevVr5+fnAwA6dOjQ4Hh0ickOEZGWVR8zpzlHRyb9kJqaiueeew6TJk3CBx98UO/+AQEBuH//PpKSktCvXz8AwC+//ILKykr4+/s3+LoDBw7Ezp07UVlZCVPTJyW6V69ehbOzs9JEp+ra7777LsrLy8Wi5ri4OPj4+IiPsIAnvVRt2rRB9+7dGxyPLrFAmYhIh6p6epKSkhQWPtYyDikpKRgyZAiGDh2K6Oho5OTkICcnB3/++ae4z9mzZ9GtWzfcvn0bwJPel2HDhiE8PBxnz57Fr7/+iunTp2PcuHENfhMLAN544w3k5+dj5syZuHr1Kg4cOIBly5YhMjJS3OeTTz5BUFCQuD5hwgSYmZkhLCwMqamp+Prrr7F27VpER0crnPvkyZN49tlnxcdZ+o49O0REOtbcoyO3RJY99KPHYc+ePfjzzz8RGxuL2NhYsd3T01N8k6m4uBjp6ekoLy8Xt+/YsQPTp09HUFAQTE1NMWbMGKxbt07h3CYmJtiyZQsmT56s9Nru7u44fPgwoqKi0KtXL7i6umLmzJmYO3euuE9eXh4yMzPFdRsbG/z888+IjIxEv3794ODggIULFyq8dg4Au3btwuLFi9X8rmgfkx0iIqJmsnjx4nqTgsDAQAiCoNBmZ2dX5wCCWVlZaN26NQYOHFjnuQMCAnD69OlGxderVy+cPHlS5TE//fQTTE1N8dJLL9V5bX3Cx1hEREQa8umnn8LKygqXL19u1uscPHgQERER6NKlS7NeR5mHDx9iy5YtaN3acPpLDCdSIiIiPbZjxw6UlJQAQLM/lqxed6NthtSjU4XJDhERkQa4urrqOgRSgckOEdWp+pgvHP+FiAwRkx0iUorjvxCRsWCyQ0RK1RzptwrncSIiQ8Nkh4hU4vgvRGQMmOwQGQG5XK7QA8PeFyKi/2GyQ2Tg5HI5ZDKZwuzDnEWbiOh/OKggkYHLy8tDcXExYmNjkZSUhNjYWBQXF9eqtSGi5hUYGAgTExOYmJggOTlZ1+Ho1Lhx47By5UpdhyHSabJz4sQJjBo1Ci4uLjAxMcG+ffvEbeXl5Zg7dy569uyJtm3bwsXFBa+++iru3LmjcA4vLy/xH1fV8uGHH2r5Toh0TyaTwdfXFzKZTNeh6FTqvdRaC7U8d28UNnlRR3h4OLKzs9GjRw+F9q1bt6JXr16wsLCAo6NjvYMCPnr0CJGRkbC3t4eVlRXGjBmD3NzcRsWSnZ2NCRMmoGvXrjA1NcWsWbOU7rd7925069YNFhYW6NmzJw4ePFjvuePj4+Hr6wtzc3N07twZW7duVdi+YMECfPDBBygoKGhUzM1Fp8nOw4cP0bt3b2zYsKHWtuLiYly4cAHvvfceLly4gL179yI9PR1///vfa+27dOlSZGdni8uMGTO0ET4RGbC0tDRcuHABFy5cMKrxg6rf14ULFyCXy3UdUosikUgglUoVplJYtWoV3n33XcybNw+pqak4cuQIQkJC6jxPVFQUfvzxR+zevRvHjx/HnTt38OKLLzYqltLSUnTo0AELFixA7969le6TkJCA8ePHIywsDBcvXsTo0aMxevRopKSkqDxvVlYWRo4ciSFDhiA5ORmzZs3CtGnTcPjwYXGfHj16oFOnTgqTn+qSTmt2hg8fjuHDhyvdZmNjg7i4OIW2Tz75BH5+fpDL5Qq1CO3atYNUKm3WWInIOBjr+EF13Rfrt3Tnr7/+woIFC/Djjz8iKChIbO/Vq5fKYwoKCvDVV19h586deO655wAAW7ZsgUwmw+nTp/G3v/2tQdf28vLC2rVrAQCbN29Wus/atWsxbNgwzJkzBwDw/vvvIy4uDp988gk2bdqk9JhNmzbB29tbfEwlk8lw6tQprF69WiGJGzVqFHbt2qXTqS2qGFTNTkFBAUxMTGBra6vQ/uGHH8Le3h59+/bF8uXL8fjx4zrPU1paisLCQoWFiFqGqvGDkpKSFBZDTwiU3Rfrt3QvLi4OlZWVuH37NmQyGdzc3PDyyy/j5s2bKo9JSkpCeXk5goODxbZu3brBw8MDiYmJGo0vMTFR4ToAEBISUud1GnqMn58fzp49i9LSUs0FrCaDeRvr0aNHmDt3LsaPHw9ra2ux/c0334Svry/s7OyQkJCA+fPnIzs7G6tWrVJ5rpiYGCxZskQbYRORHjLW8YOM9b4M2bVr11BZWYlly5Zh7dq1sLGxwYIFC/D888/j0qVLMDMzq3VMTk4OzMzMav3H3snJCTk5ORqNLycnB05OTo26jqpjCgsLUVJSAktLSwCAi4sLysrKkJOTA09PT43G3VgGkeyUl5fj5ZdfhiAI2Lhxo8K26Oho8etevXrBzMwMr732GmJiYmBubq70fPPnz1c4rrCwEO7u7s0TPBGRHqk5JpMx1Svpo8rKSpSXl2PdunUYOnQoAOC///0vpFIpjh07Vm/tjiGrSnqqD4uhK3qf7FQlOjdu3MAvv/yi0KujjL+/Px4/fozr16/Dx8dH6T7m5uYqEyEiImOlbEwmQHm9Us0kiANVqsfZ2RkA8NRTT4ltHTp0gIODg8ricalUirKyMty/f1+hdyc3N1fj9alSqbTWW171XUfVMdbW1mKCAwD5+fkAntyvrul1slOV6Pzxxx84duwY7O3t6z0mOTkZpqamcHR01EKERESGo/qYTNWHKKieyLDQWbMGDhwIAEhPT4ebmxuAJ0lAXl6eykc7/fr1Q5s2bXD06FGMGTNGPF4ulyMgIECj8QUEBODo0aMKr6XHxcXVeZ2AgIBar6crOyYlJQVubm56Ufiv02SnqKgIGRkZ4npWVhaSk5NhZ2cHZ2dnvPTSS7hw4QL279+PiooK8RminZ0dzMzMkJiYiDNnzmDIkCFo164dEhMTERUVhdDQULRv315Xt0VEpNeqxmRSRtkEsGlpaQgNDUVeXh6TnUbq2rUrXnjhBcycOROff/45rK2tMX/+fHTr1g1DhgwBANy+fRtBQUHYvn07/Pz8YGNjg7CwMERHR8POzg7W1taYMWMGAgICGvwmVpWqwQ2Liorw559/Ijk5GWZmZmJP08yZMzF48GCsXLkSI0eOxK5du3D+/Hl8/vnn4jnmz5+P27dvY/v27QCA119/HZ988gnefvttTJ06Fb/88gu++eYbHDhwQOHaJ0+eFB/d6ZpOk53z58+LHzbwv/qbSZMmYfHixfjhhx8AAH369FE47tixYwgMDIS5uTl27dqFxYsXo7S0FN7e3oiKilKoxyEiosYxhkJnR8+6Sx60afv27YiKisLIkSNhamqKwYMH49ChQ2jTpg2AJ08x0tPTFR4vrl69GqamphgzZgxKS0sREhKCTz/9VOG8Xl5emDx5MhYvXqzy2n379hW/TkpKws6dO+Hp6Ynr168DAAYMGICdO3diwYIFeOedd9ClSxfs27dPYVDE7OxshUdu3t7eOHDgAKKiorB27Vq4ubnhyy+/VKg/evToEfbt24dDhw6p9T3TNJ0mO4GBgRAEQeX2urYBgK+vL06fPq3psIiIiDTG2toaX331Fb766iul2728vGr9vbOwsMCGDRuUDroLPCn6zc3NRWBgYJ3Xru/vKACMHTsWY8eOVbm95ujIwJO/3xcvXlR5zJYtW+Dn59fonqjmYlDj7BAREemzTz/9FFZWVrh8+XKzXufYsWN47rnn6k12dKVNmzZYv369rsMQ6XWBMhGpr/rbNHy9mKj57dixAyUlJQDQ7I8BR44ciZEjRzbrNZpi2rRpug5BAZMdIiNjrNMhEOk7V1dXXYdAKjDZITIyyt6mAThOChG1XEx2iIyQMbxNQ0SkKUx2iIiMXFXNFmu3qKViskNEZKSU1W+xdotaIiY7RERGSln9Fmu3qCViskNEZMRYv0XEQQWJiIgMzuTJkzF69Ghdh2EwmOwQERFpQGBgoMLs4c11DDUeH2MRGRi5XF5rRmoiIlKNPTtEBkQul0Mmk6Ffv37iEhoaqts3bO5crL0QtTCTJ0/G8ePHsXbtWpiYmMDExATXr1/H8ePH4efnB3Nzczg7O2PevHl4/PhxncdUVFQgLCwM3t7esLS0hI+PD9auXavjOzRs7NkhMiB5eXkoLi5GbGwsZDKZ2M43bIj+p7i4GFeuXNHIubp16waJRFLvfmvXrsXVq1fRo0cPLF26FABQUVGBESNGYPLkydi+fTuuXLmC8PBwWFhYYPHixUqP6dChAyorK+Hm5obdu3fD3t4eCQkJiIiIgLOzM15++WWN3FdLw2SHyADJZDL4+vrqOgwivXTlyhX069dPI+dKSkpq0M+ajY0NzMzMIJFIIJVKAQDvvvsu3N3d8cknn8DExATdunXDnTt3MHfuXCxcuFDpMQDQqlUrLFmyRFz39vZGYmIivvnmGyY7amKyQ0RERqVbt25ISkrS2LnUlZaWhoCAAJiYmIhtAwcORFFREW7dulVnb+yGDRuwefNmyOVylJSUoKysDH369FE7lpaOyQ4RERkViURi0D2fu3btwuzZs7Fy5UoEBASgXbt2WL58Oc6cOaPr0AwWkx0iIiINMDMzQ0VFhbguk8nw7bffQhAEsXfn119/Rbt27eDm5qb0mKp9BgwYgH/9619iW2ZmphbuwHjxbSwiIlJLfnYR7t4oxN0bhcjPLtJ1ODrn5eWFM2fO4Pr168jLy8O//vUv3Lx5EzNmzMCVK1fw/fffY9GiRYiOjoapqanSYyorK9GlSxecP38ehw8fxtWrV/Hee+/h3LlzOr47w8Zkh4iISANmz56NVq1a4amnnkKHDh1QXl6OgwcP4uzZs+jduzdef/11hIWFYcGCBSqPkcvleO211/Diiy/ilVdegb+/P+7du6fQy0ONx8dYREREGtC1a1ckJiYqtHl5eeHs2bONOgYAtmzZgi1btii0xcTEiF9v3bq1acG2MEx2iPQYR0smImo6JjtEeqpqtOTi4mKFdp2OlkxEZICY7BBpQM0eGE2MaMzRkomINIPJDlETKeuBkUgkSEtL00hSwtGSiYiahskOURPV7IFJS0tDaGgo8vLyGp3sVO8hYn0OEZFmMNkh0pCm9sCo6iFifU7Tpd5LVVjvbt9dR5EQkS4w2SHSE8pqdFifQ0TUdEx2iPSM1mp07lys3ebSt/mvS0SkZUx2iHSEY+gQEWkHkx0iHeAYOkRE2sNkh0gHOIYOEZH2MNkh0iGOoUNE1Pw46zkREREZNfbsEBGBY/EQGTO1ena2bduGAwcOiOtvv/02bG1tMWDAANy4cUNjwRERERE1lVrJzrJly2BpaQkASExMxIYNG/Dxxx/DwcEBUVFRGg2QiKguqfdSFRZ9U5KSWmshIu1S6zHWzZs30blzZwDAvn37MGbMGERERGDgwIEIDAzUZHxERoPzXhER6YZaPTtWVla4d+8eAODnn3/G888/DwCwsLBASUlJg89z4sQJjBo1Ci4uLjAxMcG+ffsUtguCgIULF8LZ2RmWlpYIDg7GH3/8obBPfn4+Jk6cCGtra9ja2iIsLAxFRUXq3BZRs6kaV6dfv37o168fQkNDOaYOEZGWqJXsPP/885g2bRqmTZuGq1evYsSIEQCA1NRUeHl5Nfg8Dx8+RO/evbFhwwal2z/++GOsW7cOmzZtwpkzZ9C2bVuEhITg0aNH4j4TJ05Eamoq4uLisH//fpw4cQIRERHq3BZRs6k+rk5SUhKSkpKQlpbGMXWIiLRArcdYGzZswIIFC3Dz5k18++23sLe3BwAkJSVh/PjxDT7P8OHDMXz4cKXbBEHAmjVrsGDBArzwwgsAgO3bt8PJyQn79u3DuHHjkJaWhkOHDuHcuXPo378/AGD9+vUYMWIEVqxYARcXF3Vuj0gtf+Q+QJtbBfgj94HKffR+XJ2a82UpmytL2ZxaRER6TK1kp7CwEOvWrYOpqWLH0OLFi3Hz5k2NBJaVlYWcnBwEBweLbTY2NvD390diYiLGjRuHxMRE2NraiokOAAQHB8PU1BRnzpzBP/7xD6XnLi0tRWlpqcL9EBERkXFSK9nx9vZGdnY2HB0dFdrz8/Ph7e2NioqKJgeWk5MDAHByclJod3JyErfl5OTUiqF169aws7MT91EmJiYGS5YsaXKMRKRd+vi2FRHpP7VqdgRBUNpeVFQECwuLJgWkDfPnz0dBQYG4aKo3ioiIiPRPo3p2oqOjAQAmJiZYuHAhJBKJuK2iogJnzpxBnz59NBKYVCoFAOTm5sLZ2Vlsz83NFa8hlUpx9+5dheMeP36M/Px88XhlzM3NYW5urpE4iYiISL81Ktm5ePFJYaIgCLh8+TLMzMzEbWZmZujduzdmz56tkcC8vb0hlUpx9OhRMbkpLCzEmTNn8MYbbwAAAgICcP/+fSQlJaFfv34AgF9++QWVlZXw9/fXSBxERERk2BqV7Bw7dgwAMGXKFKxduxbW1tZNunhRUREyMjLE9aysLCQnJ8POzg4eHh6YNWsW/v3vf6NLly7w9vbGe++9BxcXF4wePRrAkzdbhg0bhvDwcGzatAnl5eWYPn06xo0bxzexiIiICICaBcpbtmzRyMXPnz+PIUOGiOtVj8kmTZqErVu34u2338bDhw8RERGB+/fv45lnnsGhQ4cU6oJ27NiB6dOnIygoCKamphgzZgzWrVunkfiIiOh/qkb+bugI4NVHDa/i4ODA8aVI69RKdh4+fIgPP/wQR48exd27d1FZWamw/dq1aw06T2BgoMpiZ+BJbdDSpUuxdOlSlfvY2dlh586dDQuciMiI1Jxny7JH88zU7uDgAIlEgtDQ0P9dy1ICu/b2Ko+pGjW8uLhYoV0ikXBATdI6tZKdadOm4fjx4/jnP/8JZ2dnmJiYaDouIiLSEx4eHkhLS1PspXlkDjdXd5XHVB81XCaTAXjSIxQaGoq8vDwmO6RVaiU7P/30Ew4cOICBAwdqOh4iMkbKRl1WNjoz6S0PDw+FBOXujYYNxqr3o4ZTi6DWODvt27eHnZ2dpmMhIiIi0ji1enbef/99LFy4ENu2bVMYa4eI9BjntDJoNetziKjh1Ep2Vq5ciczMTDg5OcHLywtt2rRR2H7hwgWNBEdERIalsW9sEWmDWslO1Tg3RESGyBDm2DK0nhy79va13tiSSCRwcHDQYVRET6iV7CxatEjTcRARkQFzc3Wv9cYWx9QhfaFWsgMA9+/fx549e5CZmYk5c+bAzs4OFy5cgJOTE1xdXTUZIxHpCut8qBFqvrFFpC/USnYuXbqE4OBg2NjY4Pr16wgPD4ednR327t0LuVyO7du3azpOIiIiIrWo9ep5dHQ0Jk+ejD/++ENh6oYRI0bgxIkTGguOiIiIqKnU6tk5d+4cPvvss1rtrq6uyMnJaXJQRMag+tsofDOFiEh31Ep2zM3NUVhYe/TMq1evokOHDk0OisiQtber/VYKwDdTiIh0Ra1k5+9//zuWLl2Kb775BsCTCTvlcjnmzp2LMWPGaDRAIkPjrOStFIBvphAR6Yragwq+9NJLcHR0RElJCQYPHoycnBwEBATggw8+0HSMRAaHb6U0nrKxb7rbN88s3kTUsqiV7NjY2CAuLg6nTp3CpUuXUFRUBF9fXwQHB2s6PiIiIqImUSvZuXnzJtzd3fHMM8/gmWee0XRMRGToOD4PEekRtV499/LywuDBg/HFF1/gr7/+0nRMRERERBqjVrJz/vx5+Pn5YenSpXB2dsbo0aOxZ88elJaWajo+IqNx+VaBwkItV0lKqsJCRM1LrWSnb9++WL58OeRyOX766Sd06NABERERcHJywtSpUzUdIxEREZHa1Ep2qpiYmGDIkCH44osvcOTIEXh7e2Pbtm2aio2IiIioyZqU7Ny6dQsff/wx+vTpAz8/P1hZWWHDhg2aio2IiIioydR6G+uzzz7Dzp07cerUKchkMkycOBHff/89PD09NR0fERERUZOolez8+9//xvjx47Fu3Tr07t1b0zERERERaYxayY5cLsepU6ewfPlyXLt2Dbt374arqyv+85//wNvbm2PvEJFRqjnKM0d4rtvdG4pzKOZnF+koEmrp1KrZ2bt3L0JCQmBpaYkLFy6Ir5wXFBRg2bJlGg2QiIiIqCnUSnb+/e9/Y9OmTfjiiy/Qpk0bsX3gwIG4cOGCxoIjImoJao67w7F3iDRLrWQnPT0dgwYNqtVuY2OD+/fvNzUmIiIiIo1RK9mRSqXIyMio1X7q1Cl07NixyUERERERaYpayU54eDhmzpyJM2fOwMTEBHfu3MGOHTswe/ZsvPHGG5qOkYiIiEhtar2NNW/ePFRWViIoKAjFxcUYNGgQzM3NMXv2bMyYMUPTMRIRUQMoq/Wx7NH4N8ZqvkUFAI6e1mrFRKQP1Ep2TExM8O6772LOnDnIyMhAUVERnnrqKVhZWWk6PiIiIqImUSvZqWJmZoannnpKU7EQERERaVyT5sYiIiIi0ndMdoiIiMioNekxFhERGT5lBclExoQ9O0RERGTU2LNDpCOXbxXUu09PNxstRKK/ak68SUSkDvbsEBERkVFjzw4RkRIN6VVStk93+8YP4kdEzUvve3a8vLxgYmJSa4mMjAQABAYG1tr2+uuv6zhqIiIi0hd637Nz7tw5VFRUiOspKSl4/vnnMXbsWLEtPDwcS5cuFdclEolWYySqqSH1OC3enYuK6+ZmuomDiIye3ic7HTp0UFj/8MMP0alTJwwePFhsk0gkkEqlDT5naWkpSktLxfXCQr52SUREZKz0/jFWdWVlZYiNjcXUqVNhYmIitu/YsQMODg7o0aMH5s+fj+Li4jrPExMTAxsbG3Fxd3dv7tCJqIUoSUlVWIhI9/S+Z6e6ffv24f79+5g8ebLYNmHCBHh6esLFxQWXLl3C3LlzkZ6ejr1796o8z/z58xEdHS2uFxYWMuEhIiIyUgaV7Hz11VcYPnw4XFxcxLaIiAjx6549e8LZ2RlBQUHIzMxEp06dlJ7H3Nwc5ubmzR4vERER6Z7BPMa6ceMGjhw5gmnTptW5n7+/PwAgIyNDG2ERERGRnjOYZGfLli1wdHTEyJEj69wvOTkZAODs7KyFqIiIiEjfGcRjrMrKSmzZsgWTJk1C69b/CzkzMxM7d+7EiBEjYG9vj0uXLiEqKgqDBg1Cr169dBgxGTO5XI68vDxxPS0tTYfREBFRfQwi2Tly5AjkcjmmTp2q0G5mZoYjR45gzZo1ePjwIdzd3TFmzBgsWLBAR5GSsZPL5ZDJZLXe+LOwlKC9nb2OoiKduZJZu81Gea0gEemOQSQ7Q4cOhSAItdrd3d1x/PhxHURELVVeXh6Ki4sRGxsLmUz2v/ZyMzi78o0+IiJ9ZBDJDpG+kclk8PX1Fdc5YjIRkf4ymAJlIiIiInWwZ4eIdCK1MEuxoYOPbgIhIqPHZIeIiLSq+huMDg4O8PDw0GE01BIw2SEyMMrqg3q62eggEqLGsWtvD4lEgtDQULFNIpEgLS2NCQ81KyY7RESkFW6u7khLSxPHqUpLS0NoaCjy8vKY7FCzYrJDRERa4+HhwcSGtI7JDlE9qo+YzNGSiTSv5s8V63hI05jsENVB2YjJEokEDg4OOoyKyDg4ODjUquEBWMdDmsdkh6gOykZM5v86iTTDw8NDoYYHYB0PNQ8mO0QNUH3E5Mu3CvRvxOQ7FxXXXfrqJg7SeyUpqbUb2+luqhPW8JA2cARlIiIiMmrs2SFqgD9yH6CNDnpz9K4HqaVTNss5NcndG4UK6/nZRTqKhIwZe3aIiIjIqLFnh4hIg64VKPb+dLTppKNItK9mLw2RvmDPDhERERk19uwQEVG92GtDhow9O0RERGTUmOwQERGRUWOyQ0REREaNNTtERuCPu4pjk3Rx0VEgTfFneu22Dj7aj8PIKB0xmaiFYbJDRER6p/pM6JyPjpqKyQ6RMao5VxaRgbBrb19rJnTOgk5NxWSHiIj0hpuru8JM6JwFnTSByQ6REapZwwMAXRytdBAJkaKGjNfDmdBJ0/g2FhERERk19uwQkf6q+YaWNt/O4gznOlO994ezoJMmsGeHiIiIjBp7doiIWrjSzGsK6+adOuooEqLmwZ4dIiIiMmrs2SEiakHy7pbrOgQirWOyQ1SNXC4Xx/cAFEdxJSIiw8Rkh+j/k8vlkMlkKC4uVmi3sJSgvZ29jqIiIqKmYrJD9P/l5eWhuLgYsbGxkMlk/2svN4Ozq7sOIyMioqZgskNUQ2s7N7Rx7CSuO+swFtIijqtDZLT4NhYREREZNSY7REREZNSY7BAREZFR0+tkZ/HixTAxMVFYunXrJm5/9OgRIiMjYW9vDysrK4wZMwa5ubk6jJhIf/1xt0hh0abUwqxaC2lH3t1yhYWoJdLrZAcAunfvjuzsbHE5deqUuC0qKgo//vgjdu/ejePHj+POnTt48cUXdRgtERER6Ru9fxurdevWkEqltdoLCgrw1VdfYefOnXjuuecAAFu2bIFMJsPp06fxt7/9TduhEtH/x54bItInet+z88cff8DFxQUdO3bExIkTIZfLAQBJSUkoLy9HcHCwuG+3bt3g4eGBxMTEOs9ZWlqKwsJChYWIiIiMk14nO/7+/ti6dSsOHTqEjRs3IisrC88++ywePHiAnJwcmJmZwdbWVuEYJycn5OTk1HnemJgY2NjYiIu7OweMIyIiMlZ6/Rhr+PDh4te9evWCv78/PD098c0338DS0lLt886fPx/R0dHiemFhIRMeIiIiI6XXPTs12draomvXrsjIyIBUKkVZWRnu37+vsE9ubq7SGp/qzM3NYW1trbAQERGRcTKoZKeoqAiZmZlwdnZGv3790KZNGxw9elTcnp6eDrlcjoCAAB1GSURERPpErx9jzZ49G6NGjYKnpyfu3LmDRYsWoVWrVhg/fjxsbGwQFhaG6Oho2NnZwdraGjNmzEBAQADfxCIyVn+m127r4KP9OBrhWkHtObc62nRSsicRNRe9TnZu3bqF8ePH4969e+jQoQOeeeYZnD59Gh06dAAArF69GqamphgzZgxKS0sREhKCTz/9VMdRExGRpqWlpSmsOzg4wMPDQ0fRkKHR62Rn165ddW63sLDAhg0bsGHDBi1FRERE2mTX3h4SiQShoaEK7RKJBGlpaUx4qEH0OtkhIqKWzc3VHWlpacjLyxPb0tLSEBoairy8PCY71CBMdoiISK95eHgwqaEmYbJDREQGqXodD2t4qC5MdogMjEXeJV2HQKRVd2/UmNLnkXmtOh7W8FBdmOwQEZFBqVnHwxoeqg+THSI9xl4cIuUsBFu42dsCAPJti3QbDOk9gxpBmYiIiKix2LNDRGSg8u6WK6w7OLbRUSRE+o3JDrVocrlc4bk/EREZHyY71GLJ5XL4dJPhUUmx2GZhKUF7O3sdRkVERJrGZIdarLy8PDwqKUbMus/h3bkrAKC9nT2cXd11HBkREWkSkx1q8bw7d8VTPfvoOgyiJqtZw0NET/BtLCIiIjJq7NkhoiZJLczSdQikYaWZ12q1mXfqqINIiDSDPTtERERk1JjsEBERkVFjskNERERGjckOERERGTUWKBMRkVGoOQq6g4MDZ0EnAEx2iKiaP+4qzh7dxdFKR5EQNZxde3tIJBKEhoYqtEskEqSlpTHhISY7RERk2Nxc3ZGWlibOcwc86eUJDQ1FXl4ekx1iskNERIbPw8ODSQ2pxAJlIiIiMmrs2SEiIoN390ahwnp+dpGKPaklYrJDRET1qjmFBKePIEPCx1hERERk1JjsEBERkVHjYywiIjJa1Qca5CCDLReTHSIiMjp27e1haak40CAHGWy5mOwQEZHRcXN1x6kjZ5H/1z0AwNWMq4icFc5BBlsoJjtEZNj+TFdc7+Cjmzia4FpBZq22jjaddBCJcXFzdYebq7uuwyA9wAJlIiIiMmrs2SHSIxZ5l7R2rZqTfhI1FcfiIX3Fnh0iIiIyauzZISLjd6VGTUw3/a+HybtbrrDu4NhGR5EQGT727BAREZFRY88OkY5osz6nRan5dhYAY/hVV7Onh4gazvB/AxA1kFwuR15enrhefWRVIiIyXkx2qEWQy+Xw6SbDo5JihXYLSwna29nrKCoiItIGva7ZiYmJwdNPP4127drB0dERo0ePRnq6Yhd1YGAgTExMFJbXX39dRxGTvsrLy8OjkmLErPscuw7Gi8v3x87AmYOOEREZNb3u2Tl+/DgiIyPx9NNP4/Hjx3jnnXcwdOhQ/P7772jbtq24X3h4OJYuXSquSyQSXYRLBsC7c1c81bOPrsMwaKmFWboOoelqvp1FREZNr5OdQ4cOKaxv3boVjo6OSEpKwqBBg8R2iUQCqVTa4POWlpaitLRUXC8sLGx6sKR3qtfosD6HiKjl0utkp6aCggIAgJ2dnUL7jh07EBsbC6lUilGjRuG9996rs3cnJiYGS5YsadZYSbeU1eiwPodIc2qOlkykzwwm2amsrMSsWbMwcOBA9OjRQ2yfMGECPD094eLigkuXLmHu3LlIT0/H3r17VZ5r/vz5iI6OFtcLCwvh7s66DWNSvUbHu3NXAEB7O3vW5xARtUAGk+xERkYiJSUFp06dUmiPiIgQv+7ZsyecnZ0RFBSEzMxMdOqkfJRUc3NzmJubN2u8pB9Yo0NERHr9NlaV6dOnY//+/Th27Bjc3Nzq3Nff3x8AkJGRoY3QiIiISM/pdc+OIAiYMWMGvvvuO8THx8Pb27veY5KTkwEAzs7OzRwdEemljNu129p7aj+OOlwr4NtgRNqk18lOZGQkdu7cie+//x7t2rVDTk4OAMDGxgaWlpbIzMzEzp07MWLECNjb2+PSpUuIiorCoEGD0KtXLx1HT0RERPpAr5OdjRs3AngycGB1W7ZsweTJk2FmZoYjR45gzZo1ePjwIdzd3TFmzBgsWLBAB9ESETXMg3uKFQTt7Ct1FAlRy6DXyY4gCHVud3d3x/Hjx7UUDekzzntFRA1R83eDg4MDPDw8dBQNaYteJztEDcF5r4ioPnbt7WFpKUFoaKhCu0QiQVpaGhMeI8dkhwyesjF1AI6rQ0T/4+bqjlNHziL/r3ti29WMq4icFY68vDwmO0aOyQ4ZDY6pQ0R1cXN1hxv/A9QiGcQ4O0RERETqYs8OEan0x92i2o0W2o/D2NV8OwsArK10EAiRkWLPDhERERk19uwQkUqZxfJabe4WljqIpA7KRkwmIqqGPTtERERk1JjsEBERkVHjYywySNVHTOZoyUREVBcmO2RwlI2YzNGSqU5/3VBc17NZ0JW5XXSr3n1crdy0EAmR4WOyQwZH2YjJHC2ZiIhUYbJDBkufR0y2yLtUq+2RQy8dREJERCxQJiIiIqPGnh0iLVHW20NERM2PyQ4REbVo1d/odHBw4AzoRojJDhE1yc38klpt7nZ6NsoykRJ27e1haSlBaGio2CaRSJCWlsaEx8gw2SEiohbJzdUdp46cRf5f9wAAVzOuInJWOPLy8pjsGBkmO0QtlLJ5rzpJ6v8Fr6wnh3Sj5lg8+jbuTmnmtVpt5p066iAS1dxc3eHGYSuMHpMd0nvVR0sGOGIyETWv+n7HsK7H8DDZIb2mbLRkgCMmNxdlvT1ELYWyGh5lWNdjeJjskN6pOe9VzdGSAY6YTESaV7OGp4qds5X4dVpaGkJDQ1nXY2CY7JBeUTXvla9fAJMbY5dxu3ZbZ1ftXd8A588yRDXreAythifftkiL0ZCmMNkhvcJ5r4iISNOY7JBe0ud5r6h+Nd/Y4rg7RKRLnBuLiIiIjBp7dog0gPNeaYmyuh4ionqwZ4eIiIiMGnt2iFoIfRtDp1Zdj47iaG4P7vH/lES6xp9CIiIiMmrs2SFqJNbnGIGaY+o0dJ8GjL3Dnhwi/cOfSiIiIjJq7NkhIr2QW1haq82Jb181mb7PjE6kDUx2iIiImqj6nH6qcLZ03WGyQ0TNT1kPjZ2d9uMgagZyuRzduslQUm1OP2U4W7ruMNkhrar5vx/+T4eIDFFaWprC1yUlxdiw5gt0/f9z+tV0NeMqImeFc7Z0HWGyQ1qjbEZz/k+n5TK7nq3rEBpPz2dGr1mfo+5xrOtRza69PSwtJQgNDVVot7SU4G9PB9Q5YzrpDpMd0pqaM5pnZVzF/Dcj+D8dIjIYbq7uOHXkLPL/uqfQbtfenomOHjOaZGfDhg1Yvnw5cnJy0Lt3b6xfvx5+fn66DouUMLQZzTmuTtMpe9NK3z24b1GrrZ3tI8V9lI2pU6hGj5W1c/3nUbaPGtTt/WkupZnX1DrOvFPHes9Tcx9NcXN1Z2JjYIwi2fn6668RHR2NTZs2wd/fH2vWrEFISAjS09Ph6Oio6/CMkrI3D1h/Q0RUt+q1Psp+Z/J3a/MwimRn1apVCA8Px5QpUwAAmzZtwoEDB7B582bMmzdPx9EZH2W1N4Bx1N/oWy+Ovs1n1VCaqsep2SPkZG2u0/PU7u0xwLqjBlDW+6NvdTzq9ghp4lrq9Bgpq/Wp+TtT1VtdxvC7VdcMPtkpKytDUlIS5s+fL7aZmpoiODgYiYmJSo8pLS1Faen/fvkVFBQAAAoLCzUeX05ODnJycjR+Xl1KT0/Ho5JiLPpoLTw7dQYA3MjMwJK5M3H48GH4+PioPA4AiosfouhBIYqLHwIAkpKSUFRUpHQfbXtc9FDr16zLw3peZdVX5cWP6t9JDUWtK9U67mFxWaPP87BEUOtaammt5HMuKal/Hy16YFKk0+uro0xDv0NKHyreuzrntbG2waHvf8Ff9/MBABnXMjB73psKvzPT09NRUlKMFR+uQ+eOnVXuZ4ikUimkUqnGz1v1d1sQ6vl5FQzc7du3BQBCQkKCQvucOXMEPz8/pccsWrRIAMCFCxcuXLhwMYLl5s2bdeYKBt+zo4758+cjOjpaXK+srER+fj7s7e1hYmKisesUFhbC3d0dN2/ehLW1tcbOq0+M/R55f4bP2O+R92f4jP0em/P+BEHAgwcP4OLiUud+Bp/sODg4oFWrVsjNzVVoz83NVdllZm5uDnNzxWf2tra2zRUirK2tjfIfcHXGfo+8P8Nn7PfI+zN8xn6PzXV/NjY29e5j8LOem5mZoV+/fjh69KjYVllZiaNHjyIgIECHkREREZE+MPieHQCIjo7GpEmT0L9/f/j5+WHNmjV4+PCh+HYWERERtVxGkey88sor+PPPP7Fw4ULk5OSgT58+OHToEJycnHQal7m5ORYtWlTrkZkxMfZ75P0ZPmO/R96f4TP2e9SH+zMRhPre1yIiIiIyXAZfs0NERERUFyY7REREZNSY7BAREZFRY7JDRERERo3JThN88MEHGDBgACQSicpBCeVyOUaOHAmJRAJHR0fMmTMHjx8/rvO8+fn5mDhxIqytrWFra4uwsDBx7ihdio+Ph4mJidLl3LlzKo8LDAystf/rr7+uxcgbx8vLq1a8H374YZ3HPHr0CJGRkbC3t4eVlRXGjBlTa6BLfXD9+nWEhYXB29sblpaW6NSpExYtWoSysrI6j9P3z3DDhg3w8vKChYUF/P39cfbs2Tr33717N7p16wYLCwv07NkTBw8e1FKkjRMTE4Onn34a7dq1g6OjI0aPHi3OH6fK1q1ba31WFhY1JzDVH4sXL64Vb7du3eo8xlA+P0D57xMTExNERkYq3V/fP78TJ05g1KhRcHFxgYmJCfbt26ewXRAELFy4EM7OzrC0tERwcDD++OOPes/b2J/hxmKy0wRlZWUYO3Ys3njjDaXbKyoqMHLkSJSVlSEhIQHbtm3D1q1bsXDhwjrPO3HiRKSmpiIuLg779+/HiRMnEBER0Ry30CgDBgxAdna2wjJt2jR4e3ujf//+dR4bHh6ucNzHH3+spajVs3TpUoV4Z8yYUef+UVFR+PHHH7F7924cP34cd+7cwYsvvqilaBvuypUrqKysxGeffYbU1FSsXr0amzZtwjvvvFPvsfr6GX799deIjo7GokWLcOHCBfTu3RshISG4e/eu0v0TEhIwfvx4hIWF4eLFixg9ejRGjx6NlJQULUdev+PHjyMyMhKnT59GXFwcysvLMXToUDx8WPeEtdbW1gqf1Y0bN7QUsXq6d++uEO+pU6dU7mtInx8AnDt3TuHe4uLiAABjx45VeYw+f34PHz5E7969sWHDBqXbP/74Y6xbtw6bNm3CmTNn0LZtW4SEhODRI9WTAzf2Z1gtGpmNs4XbsmWLYGNjU6v94MGDgqmpqZCTkyO2bdy4UbC2thZKS0uVnuv3338XAAjnzp0T23766SfBxMREuH37tsZjb4qysjKhQ4cOwtKlS+vcb/DgwcLMmTO1E5QGeHp6CqtXr27w/vfv3xfatGkj7N69W2xLS0sTAAiJiYnNEKFmffzxx4K3t3ed++jzZ+jn5ydERkaK6xUVFYKLi4sQExOjdP+XX35ZGDlypEKbv7+/8NprrzVrnJpw9+5dAYBw/Phxlfuo+n2krxYtWiT07t27wfsb8ucnCIIwc+ZMoVOnTkJlZaXS7Yb0+QEQvvvuO3G9srJSkEqlwvLly8W2+/fvC+bm5sJ///tfledp7M+wOtiz04wSExPRs2dPhcENQ0JCUFhYiNTUVJXH2NraKvSUBAcHw9TUFGfOnGn2mBvjhx9+wL179xo0UvWOHTvg4OCAHj16YP78+SguLtZChOr78MMPYW9vj759+2L58uV1PnpMSkpCeXk5goODxbZu3brBw8MDiYmJ2gi3SQoKCmBnZ1fvfvr4GZaVlSEpKUnhe29qaorg4GCV3/vExESF/YEnP5eG8lkBqPfzKioqgqenJ9zd3fHCCy+o/H2jL/744w+4uLigY8eOmDhxIuRyucp9DfnzKysrQ2xsLKZOnVrnpNOG9vlVycrKQk5OjsLnY2NjA39/f5Wfjzo/w+owihGU9VVOTk6tUZyr1nNyclQe4+joqNDWunVr2NnZqTxGV7766iuEhITAzc2tzv0mTJgAT09PuLi44NKlS5g7dy7S09Oxd+9eLUXaOG+++SZ8fX1hZ2eHhIQEzJ8/H9nZ2Vi1apXS/XNycmBmZlarbsvJyUnvPrOaMjIysH79eqxYsaLO/fT1M8zLy0NFRYXSn7MrV64oPUbVz6W+f1aVlZWYNWsWBg4ciB49eqjcz8fHB5s3b0avXr1QUFCAFStWYMCAAUhNTa33Z1UX/P39sXXrVvj4+CA7OxtLlizBs88+i5SUFLRr167W/ob6+QHAvn37cP/+fUyePFnlPob2+VVX9Rk05vNR52dYHUx2apg3bx4++uijOvdJS0urt4DOkKhzz7du3cLhw4fxzTff1Hv+6vVGPXv2hLOzM4KCgpCZmYlOnTqpH3gjNOYeo6OjxbZevXrBzMwMr732GmJiYvR2OHd1PsPbt29j2LBhGDt2LMLDw+s8Vh8+w5YuMjISKSkpddazAEBAQIDCJMgDBgyATCbDZ599hvfff7+5w2y04cOHi1/36tUL/v7+8PT0xDfffIOwsDAdRqZ5X331FYYPHw4XFxeV+xja52comOzU8NZbb9WZdQNAx44dG3QuqVRaq6K86g0dqVSq8piaRVmPHz9Gfn6+ymOaSp173rJlC+zt7fH3v/+90dfz9/cH8KRXQVt/KJvyufr7++Px48e4fv06fHx8am2XSqUoKyvD/fv3FXp3cnNzm+0zq6mx93fnzh0MGTIEAwYMwOeff97o6+niM1TGwcEBrVq1qvXmW13fe6lU2qj99cH06dPFlxUa+7/7Nm3aoG/fvsjIyGim6DTL1tYWXbt2VRmvIX5+AHDjxg0cOXKk0b2hhvT5VX0Gubm5cHZ2Fttzc3PRp08fpceo8zOsFo1V/7Rg9RUo5+bmim2fffaZYG1tLTx69EjpuaoKlM+fPy+2HT58WK8KlCsrKwVvb2/hrbfeUuv4U6dOCQCE3377TcORNY/Y2FjB1NRUyM/PV7q9qkB5z549YtuVK1f0tkD51q1bQpcuXYRx48YJjx8/Vusc+vQZ+vn5CdOnTxfXKyoqBFdX1zoLlP/v//5PoS0gIEAvC1wrKyuFyMhIwcXFRbh69apa53j8+LHg4+MjREVFaTi65vHgwQOhffv2wtq1a5VuN6TPr7pFixYJUqlUKC8vb9Rx+vz5QUWB8ooVK8S2goKCBhUoN+ZnWK1YNXamFujGjRvCxYsXhSVLlghWVlbCxYsXhYsXLwoPHjwQBOHJP9IePXoIQ4cOFZKTk4VDhw4JHTp0EObPny+e48yZM4KPj49w69YtsW3YsGFC3759hTNnzginTp0SunTpIowfP17r96fKkSNHBABCWlparW23bt0SfHx8hDNnzgiCIAgZGRnC0qVLhfPnzwtZWVnC999/L3Ts2FEYNGiQtsNukISEBGH16tVCcnKykJmZKcTGxgodOnQQXn31VXGfmvcoCILw+uuvCx4eHsIvv/winD9/XggICBACAgJ0cQt1unXrltC5c2chKChIuHXrlpCdnS0u1fcxpM9w165dgrm5ubB161bh999/FyIiIgRbW1vxLch//vOfwrx588T9f/31V6F169bCihUrhLS0NGHRokVCmzZthMuXL+vqFlR64403BBsbGyE+Pl7hsyouLhb3qXl/S5YsEQ4fPixkZmYKSUlJwrhx4wQLCwshNTVVF7dQr7feekuIj48XsrKyhF9//VUIDg4WHBwchLt37wqCYNifX5WKigrBw8NDmDt3bq1thvb5PXjwQPxbB0BYtWqVcPHiReHGjRuCIAjChx9+KNja2grff/+9cOnSJeGFF14QvL29hZKSEvEczz33nLB+/Xpxvb6fYU1gstMEkyZNEgDUWo4dOybuc/36dWH48OGCpaWl4ODgILz11lsKmf2xY8cEAEJWVpbYdu/ePWH8+PGClZWVYG1tLUyZMkVMoPTB+PHjhQEDBijdlpWVpfA9kMvlwqBBgwQ7OzvB3Nxc6Ny5szBnzhyhoKBAixE3XFJSkuDv7y/Y2NgIFhYWgkwmE5YtW6bQE1fzHgVBEEpKSoR//etfQvv27QWJRCL84x//UEgg9MWWLVuU/put3slriJ/h+vXrBQ8PD8HMzEzw8/MTTp8+LW4bPHiwMGnSJIX9v/nmG6Fr166CmZmZ0L17d+HAgQNajrhhVH1WW7ZsEfepeX+zZs0SvxdOTk7CiBEjhAsXLmg/+AZ65ZVXBGdnZ8HMzExwdXUVXnnlFSEjI0PcbsifX5XDhw8LAIT09PRa2wzt86v6m1VzqbqHyspK4b333hOcnJwEc3NzISgoqNZ9e3p6CosWLVJoq+tnWBNMBEEQNPdQjIiIiEi/cJwdIiIiMmpMdoiIiMioMdkhIiIio8Zkh4iIiIwakx0iIiIyakx2iIiIyKgx2SEiIiKjxmSHiIiIjBqTHSIiIjJqTHaIiIjIqDHZISIiIqPGZIeIjM6ff/4JqVSKZcuWiW0JCQkwMzPD0aNHdRgZEekCJwIlIqN08OBBjB49GgkJCfDx8UGfPn3wwgsvYNWqVboOjYi0jMkOERmtyMhIHDlyBP3798fly5dx7tw5mJub6zosItIyJjtEZLRKSkrQo0cP3Lx5E0lJSejZs6euQyIiHWDNDhEZrczMTNy5cweVlZW4fv26rsMhIh1hzw4RGaWysjL4+fmhT58+8PHxwZo1a3D58mU4OjrqOjQi0jImO0RklObMmYM9e/bgt99+g5WVFQYPHgwbGxvs379f16ERkZbxMRYRGZ34+HisWbMG//nPf2BtbQ1TU1P85z//wcmTJ7Fx40Zdh0dEWsaeHSIiIjJq7NkhIiIio8Zkh4iIiIwakx0iIiIyakx2iIiIyKgx2SEiIiKjxmSHiIiIjBqTHSIiIjJqTHaIiIjIqDHZISIiIqPGZIeIiIiMGpMdIiIiMmr/D9POppAufRacAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ax_x = bh.axis.Regular(100, -10, 10)\n", "ax_y = bh.axis.Regular(5, -10, 10)\n", "h = bh.Histogram(ax_x, ax_y)\n", "h.fill(x, y)\n", "for i, (a, b) in enumerate(ax_y):\n", " plt.stairs(h.values()[:,i], ax_x.edges, label=f\"[{a}, {b})\",\n", " fill=True, alpha=0.2)\n", "h1 = h[:, sum]\n", "plt.stairs(h1.values(), ax_x.edges, color=\"k\", label=\"total\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"events\")\n", "plt.legend(title=\"y interval\", frameon=False, handlelength=1.2);" ] }, { "cell_type": "markdown", "id": "copyrighted-plenty", "metadata": {}, "source": [ "## Fit with conditional variable\n", "\n", "The random distribution of $x$ depends on the value of $y$. We can exploit that information in the likelihood function to obtain a more accurate estimate of the parameters." ] }, { "cell_type": "code", "execution_count": 4, "id": "aware-fantasy", "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 = 1.93e+04 Nfcn = 130
EDM = 2.25e-06 (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", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 a -0.010 0.012
1 b 0.4993 0.0022
2 sigma 0.986 0.008 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", "
a b sigma
a 0.000142 1.31e-08 5.73e-08
b 1.31e-08 4.76e-06 1.02e-08
sigma 5.73e-08 1.02e-08 7.08e-05
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 1.93e+04 │ Nfcn = 130 │\n", "│ EDM = 2.25e-06 (Goal: 0.0002) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ No Parameters at limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Below EDM threshold (goal x 10) │ Below call limit │\n", "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n", "│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n", "┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", "├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", "│ 0 │ a │ -0.010 │ 0.012 │ │ │ │ │ │\n", "│ 1 │ b │ 0.4993 │ 0.0022 │ │ │ │ │ │\n", "│ 2 │ sigma │ 0.986 │ 0.008 │ │ │ 0 │ │ │\n", "└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌───────┬────────────────────────────┐\n", "│ │ a b sigma │\n", "├───────┼────────────────────────────┤\n", "│ a │ 0.000142 1.31e-08 5.73e-08 │\n", "│ b │ 1.31e-08 4.76e-06 1.02e-08 │\n", "│ sigma │ 5.73e-08 1.02e-08 7.08e-05 │\n", "└───────┴────────────────────────────┘" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def model(xy, a, b, sigma):\n", " x, y = xy\n", " mu = a + b * y\n", " # cannot use norm.pdf from numba_stats here, because it is not vectorized in mu\n", " return norm.pdf(x, mu, sigma)\n", "\n", "nll = UnbinnedNLL((x, y), model)\n", "\n", "m = Minuit(nll, 0.0, 0.0, 2.0)\n", "m.limits[\"sigma\"] = (0, None)\n", "m.migrad()" ] }, { "cell_type": "code", "execution_count": 5, "id": "aquatic-belgium", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAG0CAYAAADU2ObLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABf7klEQVR4nO3de1iUZf4/8PfDIKAwgIAwwKBYbZabmluJJwyTxbRcbaTdTNcsyw6e0LKyX2VaqWW7Ya7Z7n5brcx1kybJoqNHVDQPuWaWpYsJo+CR0yQgM8/vD5zHOTMzzHner+ua6+I5zdyPI8Nn7vtzf25BFEURREREREEqzNcNICIiIvIkBjtEREQU1BjsEBERUVBjsENERERBjcEOERERBTUGO0RERBTUGOwQERFRUGOwQ0REREGNwQ4REREFNQY7REREFNR8GuwsWrQIt9xyC+RyOZKTkzFmzBgcOXLE5JzGxkZMnToViYmJiImJwdixY1FdXW1yzokTJ3DHHXegU6dOSE5Oxpw5c9DS0uLNWyEiIiI/Fe7LF9+6dSumTp2KW265BS0tLXjmmWeQl5eHw4cPIzo6GgAwa9YsfPrpp1i3bh3i4uIwbdo0qFQq7NixAwCg0+lwxx13QKFQYOfOnTh16hQmTpyIDh06YOHChQ61Q6/X4+TJk5DL5RAEwWP3S0RERO4jiiLq6+uRlpaGsDA7/TeiHzl9+rQIQNy6dasoiqJYU1MjdujQQVy3bp10zg8//CACEMvKykRRFMWSkhIxLCxMrKqqks5ZsWKFGBsbKzY1NTn0uhUVFSIAPvjggw8++OAjAB8VFRV2/877tGfHXG1tLQAgISEBALBv3z5cunQJubm50jnXXXcdunbtirKyMvTv3x9lZWXo1asXUlJSpHOGDx+ORx99FN9//z369u1r8TpNTU1oamqStsXLC79XVFQgNjbWI/dGRERE7lVXV4eMjAzI5XK75/lNsKPX61FQUIBBgwbhhhtuAABUVVUhIiIC8fHxJuempKSgqqpKOsc40DEcNxyzZtGiRZg/f77F/tjYWAY7REREAaatFBS/mY01depUHDp0CGvXrvX4a82dOxe1tbXSo6KiwuOvSURERL7hFz0706ZNwyeffIJt27ZBqVRK+xUKBZqbm1FTU2PSu1NdXQ2FQiGd880335g8n2G2luEcc5GRkYiMjHTzXRAREZE/8mnPjiiKmDZtGj766CNs2rQJ3bt3Nzl+0003oUOHDti4caO078iRIzhx4gQGDBgAABgwYAC+++47nD59Wjrnq6++QmxsLHr27OmdGyEiIiK/5dOenalTp2LNmjUoLi6GXC6Xcmzi4uLQsWNHxMXFYfLkyZg9ezYSEhIQGxuL6dOnY8CAAejfvz8AIC8vDz179sSf//xnvPrqq6iqqsKzzz6LqVOnsveGiIiIIIiGqUi+eHEbCUUrV67EpEmTALQWFXz88cfx73//G01NTRg+fDjefPNNkyGqX375BY8++ii2bNmC6Oho3HfffVi8eDHCwx2L5erq6hAXF4fa2lomKBMREQUIR/9++zTY8RcMdoiIiAKPo3+//WY2FhEREZEnMNghIiKioMZgh4iIiIIagx0iIiIKagx2iIiIyC1ycnJQUFDg8PmrVq2yWBLKExjsEBERUVBjsENERERBjcEOERFRkMvJycH06dNRUFCAzp07IyUlBf/85z+h1Wpx//33Qy6X45prrsFnn30mXbN161b069cPkZGRSE1NxdNPP42WlhbpuFarxcSJExETE4PU1FT85S9/sXjdpqYmPPHEE0hPT0d0dDSysrKwZcsWb9yyCQY7REREIeCdd95BUlISvvnmG0yfPh2PPvoo7r77bgwcOBD79+9HXl4e/vznP+PXX3+FRqPByJEjccstt+C///0vVqxYgbfffhsvvfSS9Hxz5szB1q1bUVxcjC+//BJbtmzB/v37TV5z2rRpKCsrw9q1a3Hw4EHcfffduP322/Hzzz979d5ZQRmsoEwUTLRaLWJiYgAADQ0NiI6O9nGLKNiNWrYdZ+qbvP66XeSR2DB9sEPn5uTkQKfTobS0FACg0+kQFxcHlUqFd999FwBQVVWF1NRUlJWVYcOGDfjwww/xww8/SEs7vfnmm3jqqadQW1uLX3/9FYmJiVi9ejXuvvtuAMD58+ehVCoxZcoUFBYW4sSJE7jqqqtw4sQJpKWlSW3Jzc1Fv379sHDhQqxatQoFBQWoqalx6d/A0b/fPl0IlIiIKNCdqW9CVV2jr5vRpt69e0s/y2QyJCYmolevXtK+lJQUAMDp06fxww8/YMCAASZrWA4aNAgNDQ2orKzEhQsX0NzcjKysLOl4QkICevToIW1/99130Ol0uPbaa03a0dTUhMTERLffnz0MdoiIiNqhizwyIF63Q4cOJtuCIJjsMwQ2er2+/Y1Da8+qTCbDvn37IJPJTI4Zel+9hcEOERFROzg6lBRIrr/+enz44YcQRVEKgnbs2AG5XA6lUomEhAR06NABu3fvRteuXQEAFy5cwE8//YRbb70VANC3b1/odDqcPn0a2dnZPrsXgAnKREREZOaxxx5DRUUFpk+fjh9//BHFxcWYN28eZs+ejbCwMMTExGDy5MmYM2cONm3ahEOHDmHSpEkIC7sSVlx77bUYP348Jk6cCLVajfLycnzzzTdYtGgRPv30U6/eD3t2iIiIyER6ejpKSkowZ84c9OnTBwkJCZg8eTKeffZZ6ZwlS5agoaEBo0aNglwux+OPP47a2lqT51m5ciVeeuklPP7449BoNEhKSkL//v1x5513evV+OBsLnI1FFEw4G4sodDj695vDWERERBTUGOwQERFRUGOwQ0REREGNwQ4REREFNQY7REREFNQY7BAREVFQY7BDREREQY3BDhEREQU1BjtEREReptVqIQgCBEGAVqv1dXOCHoMdIiKiEJWTk4OCggJfN8PjGOwQERF5mU6nk37etm2byba/2rJlCwRBQE1Nja+b4jQGO0RERF6kVqvRs2dPaXvkyJHIzMyEWq32YauCG4MdImoXf8s9cPUbs7/dBwUntVqN/Px8aDQak/0ajQb5+fkeDXi0Wi0mTpyImJgYpKam4i9/+YvJ8ffeew8333wz5HI5FAoF7r33Xpw+fRoAcPz4cQwdOhQA0LlzZwiCgEmTJgEAPv/8cwwePBjx8fFITEzEnXfeiWPHjnnsPlzBYIeIgga/MZM/0+l0mDlzJkRRtDhm2FdQUOCxIa05c+Zg69atKC4uxpdffoktW7Zg//790vFLly7hxRdfxH//+1+sX78ex48flwKajIwMfPjhhwCAI0eO4NSpU1i6dCmA1iBq9uzZ2Lt3LzZu3IiwsDDcdddd0Ov1HrkPl4gk1tbWigDE2tpaXzeFKOA0NDSIAEQAYkNDg8/a8eGHH4qCIEhtMTwEQRAFQRA//PBDu9f7y31Q8Nq8ebPF/09rj82bN7v9tevr68WIiAjxgw8+kPadO3dO7Nixozhz5kyr1+zZs0cEINbX15u0/8KFC3Zf68yZMyIA8bvvvnNX821y9O83e3aIKOD5+hszkSNOnTrl1vOccezYMTQ3NyMrK0val5CQgB49ekjb+/btw6hRo9C1a1fI5XLceuutAIATJ07Yfe6ff/4Z48aNw1VXXYXY2FhkZmY6dJ03MdghIr/QnpyZ0tJSVFZW2jwuiiIqKipQWlra3mYSuSw1NdWt57mTVqvF8OHDERsbi/fffx979uzBRx99BABobm62e+2oUaNw/vx5/POf/8Tu3buxe/duh67zJgY7RBTwfPmNmchR2dnZUCqVEATB6nFBEJCRkYHs7Gy3v/bVV1+NDh06SIEIAFy4cAE//fQTAODHH3/EuXPnsHjxYmRnZ+O6666TkpMNIiIiAJhOAjh37hyOHDmCZ599FsOGDcP111+PCxcuuL397cVgh4gCnj9/YyYykMlkUlKvecBj2C4sLIRMJnP7a8fExGDy5MmYM2cONm3ahEOHDmHSpEkIC2sNA7p27YqIiAgsW7YM//vf//Dxxx/jxRdfNHmObt26QRAEfPLJJzhz5gwaGhrQuXNnJCYm4h//+AeOHj2KTZs2Yfbs2W5vf3v5NNjZtm0bRo0ahbS0NAiCgPXr15scN3Rpmz+WLFkinZOZmWlxfPHixV6+EyLyJV9+YyZyhkqlQlFREdLS0kz2K5VKFBUVQaVSeey1lyxZguzsbIwaNQq5ubkYPHgwbrrpJgBAly5dsGrVKqxbtw49e/bE4sWL8dprr5lcn56ejvnz5+Ppp59GSkoKpk2bhrCwMKxduxb79u3DDTfcgFmzZpn8jfYXgmgto89LPvvsM+zYsQM33XQTVCoVPvroI4wZM0Y6XlVVZXH+5MmTcfToUVx11VUAWoOdyZMn46GHHpLOk8vliI6OdrgddXV1iIuLQ21tLWJjY9t3U0QhRqvVIiYmBgDQ0NDg1O+eO5/HUL8EgEmisiEAausPibvug8gRhr87AFBSUoK8vDyP9OgEO0f/fod7sU0WRowYgREjRtg8rlAoTLaLi4sxdOhQKdAxMBRAIqLQZfjGPGPGDJOCbUqlEoWFhR79xkzkLOPAZsiQIQx0PCxgcnaqq6vx6aefYvLkyRbHFi9ejMTERPTt2xdLlixBS0uL3edqampCXV2dyYOIAp9KpcLhw4el7ZKSEpSXlzPQIb8THR0NURQhiiJ7Eb3Apz07znjnnXcgl8stPrRmzJiB3/3ud0hISMDOnTsxd+5cnDp1Cn/9619tPteiRYswf/58TzeZiHyA35iJyFzABDv/+te/MH78eERFRZnsN8767t27NyIiIvDwww9j0aJFiIyMtPpcc+fONbmurq4OGRkZnmk4ERER+VRABDulpaU4cuQI/vOf/7R5blZWFlpaWnD8+HGTypDGIiMjbQZCREREFFwCImfn7bffxk033YQ+ffq0ee6BAwcQFhaG5ORkL7SMiIiI/J1Pe3YaGhpw9OhRabu8vBwHDhxAQkICunbtCqB1iGndunUWS9EDQFlZGXbv3o2hQ4dCLpejrKwMs2bNwoQJE9C5c2ev3QcRkS2c0k7kez4Ndvbu3YuhQ4dK24Y8mvvuuw+rVq0CAKxduxaiKGLcuHEW10dGRmLt2rV44YUX0NTUhO7du2PWrFl+Wb2RiIiIfMOnRQX9BYsKErnOX4oKtvd5PNUDw54dIs9x9O93QOTsEBEREbmKwQ4RtYvxCsjbtm0z2Q4kwXIfRGSJwQ4RuUytVqNnz57S9siRI5GZmQm1Wu3DVjkvWO6DiKxjsENELjEsvGm8DhUAaDQa5OfnB0ygECz3QUS2MdghIqfpdDrMnDkT1uY3GPYVFBQ4NRTki2EkT9wHEfkfBjtE5LTS0lJUVlbaPC6KIioqKlBaWmrzHK1WC0EQIAgC1qxZ49AwkvE1Wq3WL+6DiPwfgx0ictqpU6fcet6ECRN8Mozk7vsgIv/EYIeInJaamurW83w1jOTu+yAi/8Rgh4iclp2dDaVSCUEQrB4XBAEZGRnIzs5u1+t4ehjJW/dBRL7FYIcoRLgz30Umk2Hp0qUAYBEoGLYLCwshk8nc8trtHUay9frevg8i8g0GO0TkEpVKhaKiIqSlpZnsVyqVKCoqgkqlcttreXIYyZv3QUS+4dOFQIkosKlUKuTm5iIuLg4AUFJSgry8PKknxFGCIFjN2xEEAUql0uPDSO66DyLyT+zZIaJ2MQ4IhgwZ4nKA4MgwkjXmQ0vR0dEQRRGiKDq16KYr9+HrYS1fvz5RoGCwQ0Q+t3r1ag4jEZHHMNghIp8bPXo0Dh8+LG2XlJSgvLycgQ4RuQWDHSLyC+4aDvM3XE2dyPcY7BCR2zGXpBVXUyfyDwx2iIg8gKupE/kPBjtEIcpa7wt7ZNyDq6kT+RcGO0REbsbV1In8C4MdIiI342rqRP6FwQ4RkZtxNXUi/8Jgh4jIjbRaLYYOHQrAsiq0AVdTJ/IuBjtE5BOu1J8JxJo1ri6DQUTuw2CHiLzOWv0Z421Hr8nMzERxcbHH2tleXAaDyD9w1XMi8qri4mJMmDDBYlr2yZMnbV5jqFljfo1Go8GECRM80k53GD16NO68806upk7kY+zZIQoRvhwCMn6tturPmJ/vSM0aZ16/PffuyvME6zIYRIGEwQ5RCHBl2MhdiouLTV7r7NmzbV6zY8cO6WdHatbY464lG3z5b0hE7cNhLKIgZ2sIyN6wkTtZG7JqS1VVlfRze2rR2Bv+ys/Pdzh3xtf/hkTUPuzZIQpijg4BeXJIy9lABwAUCoX0s6u1aNy1ZIM//BsSUfuwZ4coiLU1BGSwY8cOjBgxwqXXiI6OdimgsWfQoEHSz9nZ2VAqldBoNBavo5j4OmQxna9c98oWREVFAhDQ1NQE8Q8vId3Ga+gaLqDi3VkoLS1FTk6OxX1otVrExMQ43OYVK1bgiSeecPh8dzDPIWLyM5F1DHaIgpijQ0DGw0a+IgiCFGwY/8GWyWRYunQppn58ArLozibXhMuTTLZrmgE0N9k8bn5t+mOrMGNjPTrt3AgA6CKPxIbpg11q/9NPP42//e1veOWVV1y63llqtRozZsyQtkeOHAmlUomlS5dyWjuRGQY7REHM0SEg42EjX0lPTzfphRq1bDtO112E5nJeTKTiN3avb6k3TXyOjY1FXV2dxXnGAVC4PAl1l4C6S40AgKq6RvRf2Br4iKIeiomvo+rdWQ7fg7emwrsrF4koVAiiu/ufA1BdXR3i4uJQW1uL2NhYXzeHyG10Oh0yMzOtDgEZM/zfNx66aWhoQHR0tEuva/w8xj02tmS/+Amaw1qHnwyq6hptnt9SfxZhYWHQ6/UAWoekjIMSQRCQnt46gGV+78ZDXzKZDKkKBarqmmBLS/1ZyMJk0Ol1Vl/LnPH9NjQ0AIBL/6a23gvDe2preFIQBCiVSpSXl3NIi4Keo3+/2bNDFMQMQ0D5+fkWQYetYSN3MM4lEUXR4rUV971uMiR1ogEA7Acc6ZcrEZ/46fs2e1tEUURlZSXmz5+PF154weT1q96dJS3Z0NoDkotRy7bjTP2V1zcOtAw9QY5+WHr6+6MjU/ErKiqkXCQi8vFsrG3btmHUqFFIS0uDIAhYv369yfFJkyZBEASTx+23325yzvnz5zF+/HjExsYiPj4ekydPlr5NERGgUqlQVFRksWyBoefD3czr0QCW60NFxCYhXH7lYUwRGyU9UuQRaDr1MzRvTsKmgoHYVDDQqWGl3/zmN1bv3XzJhg3TB2PXM8OkR6/0OKTII9BSfxYt9WeRIo+AgMuBYXRnpD+2SnooJr7ucHvcwdE8rPZM2ScKNj7t2dFqtejTpw8eeOABm+PLt99+O1auXCltR0ZGmhwfP348Tp06ha+++gqXLl3C/fffjylTpmDNmjUebTtRIFGpVMjNzTVZtqB///5ISEgAcGUmT3tn91jLJTGfMZWU2JonoxcBUa+DTnsB6WlpEIQwiwTh1qGcPJfvOzU1FTk5ORb33tZ9bZg+2GQY6UhDA4Yt3YWqukYIYTK7ic+e5mgelqtT9omCkU+DnREjRrQ53TUyMtJm8uQPP/yAzz//HHv27MHNN98MAFi2bBlGjhyJ1157zeLbHFEoM/7jfuHCBfTq1UvaHjlyJBITE02CFFuze+zlklirRyOL6WwSHNQ0Xzmm016A5s1JOOJCflB6ejpOnjxpddjIkLeSnZ1tce+uLtnQRR4JUdRLCdOy6M4QwmRST4+xtvJ6XKXVajF06FAAtnOhzO/d0ed1R64Wkb/y+5ydLVu2IDk5GZ07d8Ztt92Gl156CYmJiQCAsrIyxMfHS4EOAOTm5iIsLAy7d+/GXXfdZfU5m5qa0NR0ZXze2owNomBmrarxuXPnLM5zZnaPIZfEvCfHkJtj6MVJSkxCVFQURFGPE6cuuHwPr776KiZMmGA1FwkACgsL3ZqLZN7bk/7YKoTLk3za0+OteycKdH5dQfn222/Hu+++i40bN+KVV17B1q1bMWLECKmrvaqqCsnJySbXhIeHIyEhwW7dkEWLFiEuLk56ZGRkePQ+iPyNo0m0zlQaNuSIGHpyDA8hrPWPrqEXZ2b309j1zDCn82/MjR492qF8HE/plhwPUXtByutpqT8L8fLsMENvz6BXtmBo4U6P5PWsXr3aZ/dOFGj8umfnnnvukX7u1asXevfujauvvhpbtmzBsGHDXH7euXPnYvbs2dJ2XV0dAx4iGxyd3WOeI2LoyTHQNbT+7M6aPtZykbxVRfjLJ3Kh0+lMXvvpnTrUXoLU29Na5LDZpKfLXUaPHo0777zTJ/dOFGj8Otgxd9VVVyEpKQlHjx7FsGHDoFAocPr0aZNzWlpacP78ebsfqJGRkRaJzkRkn7XZPX/8v704p225vCWi6/T3gKjWWheGnhxzxktBuIM78nHc8doXLlzA6YrTJlPqDXk93nh9b987wFwfChwBFexUVlbi3Llz0jfIAQMGoKamBvv27cNNN90EANi0aRP0ej2ysrJ82VSioGNtds/ZhmZU11/JOBY6We/B8GRNH3/x5JNPokqjMdlnyOsxDGsNLdyJMKE1e6A9S1MQkXN8Guw0NDTg6NGj0nZ5eTkOHDiAhIQEJCQkYP78+Rg7diwUCgWOHTuGJ598Etdccw2GDx8OALj++utx++2346GHHsJbb72FS5cuYdq0abjnnns4E4vIDkeqGhufa5jdY1jCwTD76ExDa6ATJgDJ8igAwMWLF1FTUyMNWwGWS0EEI41ZoGPMMKx12igwJCLv8Wmws3fvXmkaJQApj+a+++7DihUrcPDgQbzzzjuoqalBWloa8vLy8OKLL5oMQb3//vuYNm0ahg0bhrCwMIwdOxZvvPGG1++FKNA4EvCYz+45U9+E6vpmafaR/vLlyfIo7HrmSh6doYQ70JpLMmjQIGk7lBgHfACQkJCAhpYw6d+NiLzDp8FOTk6O3Q/bL774os3nSEhIYAFBIietXr0aTz75pElvhKHOzvnz56V9SqUShYWFFrN7rBUDNGaeSxKqzGebrd68GU/v1KGqrhGnarRQTn3HZkFFInKfgMrZISL3sDWTR6vVOjS7pz3FAIOVw0UOd25p3Xd5aKuaQ1tEHsdghyhEWZvJY7zvbz91wvxvt5hcc7re9krkwcba0hn2OFrk0LwSc0RsEvQiUKmphCAInNVE5AF+XVSQiHznbEMzquoaTR7ezjUxDzjaKmzoLuaLmY4cORKZmZkoLi62eY2jRQ43TB+MTQUDoXlzEjRvTkKXmAgAVwoR3la4E/0XbsRtHipGSBSK2LNDRHYZz7QC0O5lHhylVqsxY8YMaXvkyJEeW6nd/HXNFzMFWmdbTZgwwe617SlyaG1YyxPFCIlCEYMdIrLLfKZVe1cid4StgOPk5aEfT7G1mCnQWknaMCRlj7OF/pIu9+wYhrXS09JwpqHZpBfNkWE0IrKNwQ4RAYDNGjqeEh0dbTWocCTgUCqViIqKsjjeXobFTG1xtDaRMz54sHUhY0Ml4iMNDRi4eDNqL10JrAwr0L/yyituf32iUMBgh4gAwGYNHW9zJOBwZK0uV1hbEsPbiouLce78JZPKywAgAnj8y9N2ryUi65igTEQmRL0OLfVnkSKPgCI2yqKGjqc5GnB4IjCxtiSGtz355JPSz4Y8HsPDeN0tInIce3aIQoT5sJFWq7V6Xls1dKxNyW4rL8XWkJU1jgYczgQmjr5+dnY2lEolNBqN3Xo55eXlkMlkNv8N20Oj0UBhVnnZ2oKi27Ztc9uiqq68p0SBhD07RCFo1LLtuK1wJ9IfWyVNd7ZWQ0er1UIQBAiCAK1Wa3NKtlqtdlvbDAGHrWRgQRCQkZHRWqDPzWQyGZYuXSq9jvnrAlfq5XhS1buzpKnpmjcnQae1nP02cuRIk/fCVd54T4l8jcEOUQgyzs8xTHduK0enuLgY+fn5FgteajQa5Ofnu+2Po68DDpVK5VC9HF8w5PAYHhj9crtq8RhmvXn6PSXyNQY7RCHMPD8nRR5hsXilwZNPPmlzhhQAFBQUuK3on68DDpVKhcOHD0vbJSUlKC8v90qgk56ebrtXyyyHJ1yeJNXicfbfvq1Zb4B731MiX2KwQxQEzIebHGXIz9lUMBC7nhmGTQUDLRavNDD/9m/MeIaUuzgacBjycURRdOsyC47Uy/HEa7/66qsATHu1dA0X0FJ/1uQh6k2DkP379zv1Os7MerPFVxWuiZzFYIeI3MbdM6ScLdAXDKwtO2Gew2Mtj8fZf/v2znpjrg8FEs7GIgoBo5Ztx5n6JmnbUwt6+sPU7WBgvuyEPYY8nld/jMPShRvRRR6JDdMHt3lde2a92VtSIz8/3+e5TUTm2LNDFALO1De1e0FPu7kkZjOkPDW0FCwcGf4x7sVyJI/nQpOIqrpGk6DWHldnvTHXhwIRgx2iEBImAIrYqDaTka2xlktivO2NKdnu4mqOkztYG/5pawp5W3k8cR1EhF0+JIp6h+7N1Vlv7sj1IfI2BjtEIcSwqOeuZ4Zh9//7PU69U2DS+2KvR8ZaLgngH1OyA4W9qd4G1np6bOXxhG14Dm8Mk+O/L95psjK9o1yZ9ebLCtdErmKwQ0QO8+WU7EBnb/jHmCHRt7i42GS/p/7tnX1eT1S4JvI0JigTkVNCcYaUO7Q1/GNMo9FgwoQJFvsd+bc/09AsLR56W+FOCEJYm0nLzrynji6p4YkK10SuYrBDFGT++H97cU7bYrLPU7OvyHHODOuIomgzcbgtehHSyvXV9c1Wz9FqtYiJiQEANDQ0OPX8hlyf/Px8CIJgEvAEYv4WhQYOYxEFmbMNzSYzr1ydfUWusZX35OywjqMLpxp0kUdKieeGxOUw1+KlNvm6wjWRs9izQxSkwgRYJK12kUf6qDXU1vBPexmGqYx7bfq99KXN3p22mPf+mCesm9cCKikp4Wrp5LcY7BAFKcPMK/IP9oZ/AhXztyhQcBiLiMhLbA3/WONqzg4RWWLPDhHZZF7pNy8vz4etCQ7WloKwlejrLqfrG9F/4UZpOzGaH/0UWtizQ0RW2Vro0bz+CznPeLjn/ffft5rou3r1are9nl6EScL62QbX8niIAhXDe6IAN2rZdpyuuyjVVjnjhj9k9hZ6tFb/hVw3evRo3HnnnRaJvo2N7S8XkBQTAUG48p32dD1n5lFoYrBDFODO1Dehur5Zqq1i+GNWqamEIAhWZ9LY09ZCj6GSS2KYQu4Nnkr0/eDBm03e+/4LN6KqjjWXKPQw2CEKcBcvXgQgQNTroNO2LuwpC5M5tcinMUcWevQWbwYcRBS8GOwQBTC1Wo1z5+sRLk+CTnsBmjcnAbBMeHUGF3AkomDDYIcoQBmGm/CHlyyOGQc65itot4ULOPoHT/Zqma+fBQCKia+j6t1ZHnk9Il/jbCyiAOXowpI7duxw6nkNlX5t5eYEQ86O+ZR6ZwNCf2Z+b3V1dRAE4XJvnx7AlfWzwuVJqK5vRnV9M2QxnX3VZCKPY7BDFKAcHW6qqqpy6nkNlX4By8AmGAIdW1Pq1Wq1D1vlHtbuzXg7KSbCa+tnEfkTBjtEAcrR4SaFQuH0c9tb6NGd9V+8zTClXqPRmOzXaDTIz88P6IDH1r2dPHlS+vmDB2/GrmeGYVPBQGjenATNm5PQJSbC200l8jqfBjvbtm3DqFGjkJaWBkEQsH79eunYpUuX8NRTT6FXr16Ijo5GWloaJk6caPKLCwCZmZlSF63hsXjxYi/fCZF3jFq2Hf0XbkT/hRvx9M4WdJ3+HmTR9ocfBg0a5NJrqVQqHD58WNouKSlBeXk57r33Xqurevu7tqbUA0BBQYFfD2nZWlHdkXsznGeLMl0ZcO8pkaN8GuxotVr06dMHy5cvtzj266+/Yv/+/Xjuueewf/9+qNVqHDlyBH/4wx8szl2wYAFOnTolPaZPn+6N5hN53Zn6JqNKuE0QOnWGEGZZk8V4uKk9NVuCaaFHR6bUV1RUoLS01Iutcg9n87eMg56mJlZTpuDn09lYI0aMwIgRI6wei4uLw1dffWWy729/+xv69euHEydOoGvXrtJ+uVzuUlc9UaAy1NRJT0tDY2MTampqTOrqpKenO/THL5Q4muMUiFPvncnfUqvVmDFjhrTv7LmzCJcnXa7XdIW1ddECOdil0BZQOTu1tbUQBAHx8fEm+xcvXozExET07dsXS5YsQUtLi93naWpqQl1dncmDKJAYaupsKhiI/754J757+Q/StOGSkhJ8//33Pm6h/3E0xykQp9472uZjx45ZzesBgHPnz0k5S44mcdsaViPyNwET7DQ2NuKpp57CuHHjEBsbK+2fMWMG1q5di82bN+Phhx/GwoUL8eSTT9p9rkWLFiEuLk56ZGRkeLr5RB4VTMNNnuLIlPqMjAxkZ2d7uWXt19a9GaxcudJm7R5ZdGfM2vIr+jz3CWZsrIdu2BMmx4MhiZtCV0AEO5cuXcIf//hHiKKIFStWmBybPXs2cnJy0Lt3bzzyyCP4y1/+gmXLlqGpqcnm882dOxe1tbXSo6KiwtO3QEQ+5siU+sLCwoAMFB0tF2A+wcPkvDAZhE6dUXtJQLg8yaLuTqAkcRNZ4/fBjiHQ+eWXX/DVV1+Z9OpYk5WVhZaWFhw/ftzmOZGRkYiNjTV5EFHwszelvqioCCqVykctaz9b95aenm73Ol3DBanmTkv9WYh624FMICdxU2jz6+UiDIHOzz//jM2bNyMxMbHNaw4cOICwsDAkJyd7oYVEFGhUKhVyc3MRFxcHoDXHydvJt+ZLQWi1Wrc8r7V7GzRokLRtjfkSEemPrUK4PMnu6wRiEjeFNp8GOw0NDTh69Ki0XV5ejgMHDiAhIQGpqanIz8/H/v378cknn0Cn00mVYBMSEhAREYGysjLs3r0bQ4cOhVwuR1lZGWbNmoUJEyagc2eWPqfQxZk09gVzjpP5vRlLT0/HyZMnrebtCIKApCT7QY5BICZxU2jzabCzd+9eDB06VNqePXs2AOC+++7DCy+8gI8//hgAcOONN5pct3nzZuTk5CAyMhJr167FCy+8gKamJnTv3h2zZs2SnocokI1ath1n6k1zz07XN7Z5XXFxsUmS/siRI6FUKrF06dKAHqYJZp5c9NPYq6++igkTJlxeJ+vK6xnyepYvX44nSm3nOwqCAKVSGZBJ3BTafBrs5OTk2P0Fb+uX/3e/+x127drl7mYR+QVDAUFnTZgwweJ3xzCTJtDzUqh9Ro8ejaKiIsyYMcNk+rlSqURhYSFUKhVeOvgJai9ZXhvoSdwU2vw+QZko1IUJgCI2SnqkyCNMCgiaC+TlEMjzbC0DYgiCO3bsCACQmVXmDoYkbgpdfp2gTERAl5gI7P5/uQBa89wAICYmz+nnMZ5Jk5OT484mkh9rK3/LVs5SiiIFJy7/7IskbiJ3Ys8OUYjhTJrQUVxcbLUScnFxcZvXnm1oRvpjq5D+2Cq8sD8cg17ZglHLtnuyuUQew54dIjfQarWIiYkB0Nr74s+l8zmTJnTYyt+aMGFCm9fqRUhT0KvruVgoBTYGO0RBxnymjfF+zqQJLbbyt+wtK9FFHnn5PD00lysuR8QmQe/5yWJEHsNgh8iPGPcQ9XvpS5efx9bUYmdn0nhrSjR5l733dMP0wQAs/y+yd4cCGXN2iAKMecKp+eyq1atXB+VyCERErmKwQxRAHEk4HT16tN2pxUREoYbBDlEAmTBhgkkxOOBKwumHH34IURQRHR0d1MshkONs5ebYy9khCkbM2SHyIeO8iJGFW6TpvgBwpsEyR8JewmlBQQFGjx7NwIZM2MrfIgol7Nkh8hNnG5pRXd+McHkSwuXOzX4xLhhIZGArf2v16tU+ahGRbzDYIfIzol6HlvqzSJFHQBEbhaSYCIevZcFAMmYrf2v06NE+bBWR93EYi8jP6LQXoHlzEo5cLk6o1WoRU+DYtSwY6JhgnlJvfm9arVb6ub35W6frG9F/4UZpu4s8UpqqTuTPGOwQBRAWDCRf0otAVV2jr5tB5DQGO0QBxl0FA4kclRQTAUG4kvVwur6RFZUpoDDYIfJzxkUDn3nmGaxcuRInL5fxB1oTTgsLC1lHhxzm7DDeBw/ebLLeW/+FG9nDQwGFCcpEXqLVaiEIAgRBMMmjsEetVpsUEXz55ZdNjrNgIBFR2xjsELlBW0s4uKK4uBj5+fkWRQSNZ1yxYCB5iqH3x1CokiiQuRTsvPPOO/j000+l7SeffBLx8fEYOHAgfvnlF7c1jigQmPe+GJZwUKvV7XreJ5980mYRQQN3BFVERMHOpWBn4cKF6NixIwCgrKwMy5cvx6uvvoqkpCTMmjXLrQ0k8mdqtdpq74tGo0F+fn67Ah7z57Rmx44dLj8/EVGocCnYqaiowDXXXAMAWL9+PcaOHYspU6Zg0aJFrOBKIUOn02HmzJl2e18KCgoc7n252Oh8wmdVVZXT1xARhRqXgp2YmBicO3cOAPDll1/i97//PQAgKioKFy9edF/riPxYaWkpKisrbR53dgmH8+fPO90GhULh9DVERKHGpannv//97/Hggw+ib9+++OmnnzBy5EgAwPfff4/MzEx3to/Ibzm6NIO18/74f3txTtsCUdRLC3/Kojs73YZBgwZZ3R/MFYKJiJzlUs/O8uXLMWDAAJw5cwYffvghEhMTAQD79u3DuHHj3NpAIn/l6NIM1s4729CMqrpGk4U/hTDHZlUZr1rNmVhERG1zqWenrq4Ob7zxBsLCTGOlF154ARUVFW5pGJG/y87OhlKphEajcXkJB1Gvh05rOnyla7hg93XT09PtDp8ReYv5WlkA18si/+RSsNO9e3ecOnUKycnJJvvPnz+P7t27czoshQSZTIalS5ciPz/f6SUcWpORw6DTnofmzUkOv2ZJSQkGDRqEuLi49jafqN24VhYFCpeGsWzlAjQ0NCAqKqpdDSIKJCqVCkVFRUhLSzPZr1QqUVRUZLOycW1NrcOvYTxsxSKC5Cx3FQc0rgCeGB0ORWyUySNMaPs6RyuHE7mbUz07s2fPBtD64fv888+jU6dO0jGdTofdu3fjxhtvdGsDifydSqVCbm6u1NtSUlKCvLw8u0GJTq9z6JfPONAh8hfma2UBXC+L/JtTwc63334LoLVn57vvvkNERIR0LCIiAn369METTzzh3hYSBQDjwMadvS9KpRKLFy/G+PHj3fJ8REShyKlgZ/PmzQCA+++/H0uXLkVsbKxHGkVEwOLFi/HEE0+g0YVig0REdIVLCcorV650dzuIgp5x4r75TEZrHn30UebnEBG5gUvBjlarxeLFi7Fx40acPn0aer3e5Pj//vc/tzSOKFio1WrMmDFD2tbr9VZnBxjP6rIV6LBgIBGRc1wKdh588EFs3boVf/7zn5GamsokSiI7DIuFOhKgsIYOBbpKTSUEQUBDQ0O7Zn8RuZNLwc5nn32GTz/91GapeiJqZW+xUHOsoUNE5BkuBTudO3dGQkKCu9tCFHTaWizU2JAhQzzcGiKi0ORSUcEXX3wRzz//PH799dd2vfi2bdswatQopKWlQRAErF+/3uS4KIp4/vnnkZqaio4dOyI3Nxc///yzyTnnz5/H+PHjERsbi/j4eEyePBkNDQ3taheRuzi6WCgREXmOSz07f/nLX3Ds2DGkpKQgMzMTHTp0MDm+f/9+h55Hq9WiT58+eOCBB6xWmn311Vfxxhtv4J133kH37t3x3HPPYfjw4Th8+LBUqXn8+PE4deoUvvrqK1y6dAn3338/pkyZgjVr1rhya0RuZVgEVDHxdchirqxq7soK50RE5BqXgp0xY8a45cVHjBiBESNGWD0miiIKCwvx7LPPYvTo0QCAd999FykpKVi/fj3uuece/PDDD/j888+xZ88e3HzzzQCAZcuWYeTIkXjttdcsSvgTeZthsVAxpjPC5Um+bg4RUUhyKdiZN2+eu9thoby8HFVVVcjNzZX2xcXFISsrC2VlZbjnnntQVlaG+Ph4KdABgNzcXISFhWH37t246667rD53U1MTmpqapO26ujrP3QiFLK1Wi5iYGABA+uV9ol4HnfbKqubGPxMFA1l0Z6Q/tgq3Fe6EIIQhMdqlPzNEbuVSzg4A1NTU4P/+7/8wd+5cnD9/HkDr8JVGo3FLw6qqqgAAKSkpJvtTUlKkY1VVVRYrr4eHhyMhIUE6x5pFixYhLi5OemRkZLilzUS2GBL6ddoL0Lw5CZo3JyFsw3P4S15yG1cS+Qfjopjbtm0z2TYmhMkQLk9CdX0zquoacaah2aHriDzJpWDn4MGDuPbaa/HKK6/gtddeQ01NDYDWeiJz5851Z/s8Yu7cuaitrZUeFRUVvm4SBbmOl3PMDEpKSlBeXi4N0drjrlWriVylVqvRs2dPaXvkyJHIzMyEWq2W9nWRRyJFHoGW+rNoqT8rrYJeXVVt9zoib3Ap2Jk9ezYmTZqEn3/+WUoUBlr/I2/bts0tDVMoFACA6upqk/3V1dXSMYVCgdOnT5scb2lpwfnz56VzrImMjERsbKzJg6g9nA1I3LlYKJEnGYpimvfaazQa5OfnS4HLhumDsalgoNRzGRPeWllfp9fZvY7IG1wKdvbs2YOHH37YYn96errd4SNndO/eHQqFAhs3bpT21dXVYffu3RgwYAAAYMCAAaipqcG+ffukczZt2gS9Xo+srCy3tIOIKFTZK4pp2FdQUGB1aKq2ptbqc7Z1HZEnuBTsREZGWk3q/emnn9ClSxeHn6ehoQEHDhzAgQMHALQmJR84cAAnTpyAIAgoKCjASy+9hI8//hjfffcdJk6ciLS0NGk22PXXX4/bb78dDz30EL755hvs2LED06ZNwz333MOZWBSQOGRF/qStopiiKKKiogKlpaUWx8x7dBy9jsgTXAp2/vCHP2DBggW4dOkSgNbFC0+cOIGnnnoKY8eOdfh59u7di759+6Jv374AWofH+vbti+effx4A8OSTT2L69OmYMmUKbrnlFjQ0NODzzz83GTp7//33cd1112HYsGEYOXIkBg8ejH/84x+u3BYRERlxtCimq8UzWXSTvMXlooL5+flITk7GxYsXceutt6KqqgoDBgzAyy+/7PDz5OTk2F0zSBAELFiwAAsWLLB5TkJCAgsIEhF5gKEopqPnOTss5ejzE7WXS8FOXFwcvvrqK2zfvh0HDx5EQ0MDfve735nUxCEiosBmKIqp0WisfjEVBAFKpRLZ2dlQq9WYMWOGQ89rfB2RN7gU7FRUVCAjIwODBw/G4MGD3d0mopBhyNEh8kcymQxLly5Ffn4+BEEw+b8qCK1zywsLC1FcXIz8/HyH/i8bX8cZieQtLuXsZGZm4tZbb8U///lPXLjACrBERMFKpVKhqKjIYtKHUqlEUVERRo8ebXPGljWG66yth0jkKS4FO3v37kW/fv2wYMECpKamYsyYMSgqKjJZgoGIiIKDSqXC4cOHpW1DUUyVStXmjC1jxtcReZNLwU7fvn2xZMkSnDhxAp999hm6dOmCKVOmICUlBQ888IC720gUkIyTNZuamu2cSeT/jIecjItiOjOjisU0yVdcXhsLaB17HTp0KP75z3/i66+/Rvfu3fHOO++4q21EAcu8vP7Zc2d92Boiz+GMKgoE7VqOtrKyEmvWrMGaNWtw6NAhDBgwAMuXL3dX24gC0uAFG3C8uh4Y/bK02rksurNP20TkKW3N2CLyBy4FO3//+9+xZs0abN++Hddffz3Gjx+P4uJidOvWzd3tIwooOp0OFWdrES5PavM8omBgb8YWkb9wKdh56aWXMG7cOLzxxhvo06ePu9tEFLBKS0uh0+kQDkDU66DTms5W1DW0bu/fvx85OTnebyCRBxhmbM2YMcNiwdAO8kSkP7YKAHBb4U4IQhi6yCOxYTrLlpD3uBTsnDhxAtu3b8eSJUvwv//9D+vWrUN6ejree+89dO/enbV3KGQZJ2vqtBegeXNSm+cRBQOVSoXc3FzExcUBAJISk1DTDIgQpJ7O6nom6pNvuJSgrFarMXz4cHTs2BH79++XppzX1tZi4cKFbm0gUSBxtrw+UTAxnmmVnhgDRWwUUuQRaKk/i5b6swgTfNg4CmkuD2O99dZbmDhxItauXSvtHzRoEF566SW3NY4o0GRnZ0Omtl1zhGXyKVR88ODNiI6OhlarRUxMDACg30tfsneHfMKlnp0jR45gyJAhFvvj4uJQU1PT3jYRBSyZTIb4uHirx1gmn4jIN1wKdhQKBY4ePWqxf/v27bjqqqva3SiiQNaxY0cAgCzMNKBhmXwiIt9waRjroYcewsyZM/Gvf/0LgiDg5MmTKCsrwxNPPIHnnnvO3W0kCkgpihScuPxzSUkJ8vLy2KNDROQDLgU7Tz/9NPR6PYYNG4Zff/0VQ4YMQWRkJJ544glMnz7d3W0kCkjGuZgsk09E5DsuDWMJgoD/9//+H86fP49Dhw5h165dOHPmDF588UV3t4+IiPxAdHQ0RFGEKIqIjo526TpBaNcKRUQua9dyERERESbr/xARERH5G4bZREREFNQY7BAREVFQY7BDREREQY3BDhEREQW1diUoExERGRhmXhH5GwY7RO00atl2nKlvkrZP1zf6sDVERGSOwQ5RO52pb0JVHQMcIkedrm9E/4Ubpe0u8khsmD7Yhy2iYMdgh8hNwgQgWR4lbSdG89eLyBq9CH5BIK/ipzGRmyTLo7DrmWHStlarRUyB79pD5G+6yCNNtk/XN0LPFB/yAgY7RB7CZE0iU+ZDVf0XbmQPD3kFp54TERFRUGOwQ+QCrVYLQRAgCAJEUe/r5hARkR0MdoiIiCioMdghIiKioMZgh4iIiIIagx0iIiIKagx2iIiIKKj5fbCTmZkpzXoxfkydOhUAkJOTY3HskUce8XGriYjIGs5kJF/w+6KCe/bsgU6nk7YPHTqE3//+97j77rulfQ899BAWLFggbXfq1MmrbSQiIiL/5ffBTpcuXUy2Fy9ejKuvvhq33nqrtK9Tp05QKBTebhoREREFAL8fxjLW3NyM1atX44EHHoAgCNL+999/H0lJSbjhhhswd+5c/Prrr3afp6mpCXV1dSYPIluMu921Wq2vm0MUdERRz98x8ii/79kxtn79etTU1GDSpEnSvnvvvRfdunVDWloaDh48iKeeegpHjhyBWq22+TyLFi3C/PnzvdBiClbGQ6tNTc0+bAkREbUloIKdt99+GyNGjEBaWpq0b8qUKdLPvXr1QmpqKoYNG4Zjx47h6quvtvo8c+fOxezZs6Xturo6ZGRkeK7hFFTUajVmzJghbZ89dxbh8iRcvHjRh60iIiJbAmYY65dffsHXX3+NBx980O55WVlZAICjR4/aPCcyMhKxsbEmDyJHFBcXIz8/HxqNxuLYufPn7PYoEhGRbwRMz87KlSuRnJyMO+64w+55Bw4cAACkpqZ6oVUUauZurkHaoytN9smiO0s/FxQUYPTo0ZDJZN5uGlHAOtPQjPTHVgEAbivcCUEIQxd5JDZMH+zbhlHQCIhgR6/XY+XKlbjvvvsQHn6lyceOHcOaNWswcuRIJCYm4uDBg5g1axaGDBmC3r17+7DFFKx0EdEIlyfZPF5RUYHS0lLk5OR4r1FEAU4vQvq9qq5nDhy5X0AEO19//TVOnDiBBx54wGR/REQEvv76axQWFkKr1SIjIwNjx47Fs88+66OWUqgQ9TrotBdM9ukaWrdPnTrliyYRBZykmAgIQhhEUQ/NyZMAgIjYJOhFHzeMgk5ABDt5eXkQRcv//RkZGdi6dasPWkShTqe9AM2bk6we4xAqkW3GMxmnXfsr8vLy0NjYiJiYGABAv5e+ZO8OuV3AJCgT+QNZmO1cHEEQkJGRgezsbC+2iChwqNVq9OzZU9oeOXIkMjMzUVxc7MNWUSgIiJ4dIn8RFx+HukuW+w1FLgsLC5mcTGSFWq1Gfn6+RS+9RqPBhAkTfNQqChXs2SFqg3G3e5jQ+itj3sOjVCpRVFQElUrl1bYRBQKdToeZM2daTUcQRdFkP4t0kicw2CGyw7zb/ey5swCAuLg4aV9JSQnKy8sZ6BDZUFpaisrKSofONfyOsUgnuRODHSIbDN3u1goInr9wXvp5yJAhHLoissOVGYos0knuxGCHyAp73e5E5BxXZygWFBSYDCMTuYrBDpEVznS7E5F92dnZUCqVUiK/oyoqKvDll196qFUUShjsEFnBwoBE7iOTybB06VIAcDrgqaqq8kSTKMQw2CGygoUBidxLpVKhqKgIaWlpTl2nUCg81CIKJQx2iKxwtdudiGxTqVQ4fPiwtJ2UlNTm79igQYM83SwKAQx2iKxoT7c7EdlmPHPRkd8xznQkd2CwQ2SDvW73hIQEH7SIKLiMHj3apaEtImcx2CGyw6LbPTEJANAxKspXTSIKKrZ+x4jcicEOURuMu9EjIyN82BKi4MTfMfI0LgRKZMeoZdtxuu4i0h9bBQA409C6bo8ghLHgIJEHyaI7I/2xVbitcCcEIQxd5JHYMH2wr5tFAYrBDpEdZ+qbUF3fjHB5a9e6nvENkVcIYTKEy5NQXc+FQan9OIxFZEar1UIQBAiCAFHUAwBEvQ4t9WeRIo+AIjYKXeSRPm4lUXAKb/kVovYCWurPoqX+LER963IRXBiU2oM9O0QO0GkvQPPmJBxpaEB0dLSvm0MUtMrm32UyRJz+2CqEy5OkhUFVKpUPW0eBij07RETkN+zlwnFhUHIVgx0iIvKa6OhoiKIIURSd7iWtqKhAaWmph1pGwYzBDhERBQwu0kuuYLBDREQBg4v0kisY7BARkd+wt05WRkYGsrOzvdgaChYMdogcoExXupRjQERtM+TxfPjhhwBsBzyFhYVcGJRcwmCHiIj8gr3FdxMTEjntnFzGYIfIjPHU1qYmVm8l8iabi+927OirJlEQYLBDZEStVqNnz57S9tlzZwGweiuRN3FhUHI3BjtEl6nVauTn50Oj0VgcM1RvJSKiwMNghwitQ1czZ85k9VYioiDEYIcIQGlpKSorK+2ew+qtRESBicEOERyvysrqrUREgYfBDhEcr8rK6q1ERIEn3NcNIPIH2dnZUCqV0Gg0NvN2WL2VyHdO1zdCOfUdAEB6WhoEIQxd5JHYMH2wj1tGgYDBDhFap7p2e+hNiNXnTfdHd5Z+ZvVWIt/Ri0C4vLXmTnU961+RcxjsEF3WEt4J4XLrI7ut1Vvv9HKLiCgpJgKCEAZR1ENz8iQAICI2CXrbEyeJLPh1zs4LL7wAQRBMHtddd510vLGxEVOnTkViYiJiYmIwduxYVFdX+7DFFAzCBKCl/ixa6s8iPgJQxEaia3K8r5tFFDIMa2WJoohPZ96KXc8Mw6aCgdC8OQmaNyehSwwLDZJz/L5n57e//S2+/vpraTs8/EqTZ82ahU8//RTr1q1DXFwcpk2bBpVKhR07dviiqRQkusRE4JvFkwAARxoauPgnkZ+q1FRCEAQ08PeU2uD3wU54eDgUCoXF/traWrz99ttYs2YNbrvtNgDAypUrcf3112PXrl3o37+/t5tKREREfsivh7EA4Oeff0ZaWhquuuoqjB8/HidOnAAA7Nu3D5cuXUJubq507nXXXYeuXbuirKzM7nM2NTWhrq7O5EFERETBya+DnaysLKxatQqff/45VqxYgfLycmRnZ6O+vh5VVVWIiIhAfHy8yTUpKSmoqqqy+7yLFi1CXFyc9MjIyPDgXRARUXsZL9XS1MTZWOQcvw52RowYgbvvvhu9e/fG8OHDUVJSgpqaGnzwwQftet65c+eitrZWelRUVLipxURE5G5qtRo9e/aUts+eO+vD1lAg8vucHWPx8fG49tprcfToUfz+979Hc3MzampqTHp3qqurreb4GIuMjERkZKSHW0tERO2lVquRn59vd5Feorb4dc+OuYaGBhw7dgypqam46aab0KFDB2zcuFE6fuTIEZw4cQIDBgzwYSspkGi1WqmsgSjqfd0cIjKi0+kwc+bMNgMd4yEuImv8Oth54oknsHXrVhw/fhw7d+7EXXfdBZlMhnHjxiEuLg6TJ0/G7NmzsXnzZuzbtw/3338/BgwYwJlY1C6tBcxaa3xwOiuR75SWlqKysrLN81asWMGAh+zy62CnsrIS48aNQ48ePfDHP/4RiYmJ2LVrF7p06QIAeP3113HnnXdi7NixGDJkCBQKBdRqtY9bTURE7nDq1CmHznv66aeRmZnJz3+yya9zdtauXWv3eFRUFJYvX47ly5d7qUVEROQtqampDp+r0WiQn5+PoqIiqFQqD7aKApFf9+wQEVHoys7OhlKphCAIbZ5ryOspKCjgkBZZYLBDRER+SSaTYenSpQDgcMBTUVGB0tJSTzeNAgyDHQppLFRG5N9UKhWKioqQlpbm8DWO5vpQ6GCwQyHLVqGyixcv+qpJRGSFSqXC4cOHHT7fmVwfCg1+naBM5ClqtRpTPz4B2eiXkX55nyy6MwDg3PlzUKvVTHIk8iMymcxyX3RnpD+2StrWaS+gw6a/Ijs724sto0DAnh0KOYZCZbLozgiXJ0kPIezKhymTHIn8nxAmM/kdlkV3RkVFBcLDw6HVan3dPPIjDHYo5JgXKhP1OrTUn5UeuoYLTHIk8mPdkuMhai9Iv7OivvWLSUJCgo9bRv6Kw1gUcsyTF3XaC9C8OanN84jIP3z5RC50Oh3i4uIAAH2e+wQ1zUDHqCgft4z8FXt2KOQ4mrzIJEci/2WcwxMZGeHDllAgYM8OhRxDoTJbSwsKggClUskkRyKiIMGeHQo5xoXKzBkKlxUWFlqd/UFERIGHwQ6FJJVKhcSERIv9SqWSa+sQ+aHo6GiIoghRFBEdHe3r5lCAYbBDIatjx44m2yUlJSgvL2egQxSgjIemt23bxvIRJGGwQ3TZkCFDOHRFFMCqq6qln0eOHInMzEyo1Woftoj8BYMdIiIKOMbDWo2NTQAAnd60J0ej0SA/P58BDzHYISKiwKXT6VBTW2P1mCEYGjt2LOrq6rzbMPIrDHaIiChglZaWOpSbs2PHDi+0hvwVgx0iIgpYjlY6r6qq8nBLyJ8x2CEiooDlaKVzhULh4ZaQP2OwQyFPma5k7Q6iAJWdne3QLMpBgwZ5oTXkr7hcBIWEUcu240x9k8m+0/WNPmoNEbmLTCZDfFw8ai8BsujOSH9slclxXcMFVL07i2UlQhyDHQoJZ+qbUFXH4IYoGHXs2BG1lxohhMkQLk+yes62bduQl5fHoCdEcRiLQkqYAChio0weXeSRvm4WEbVDF3kkFLFRSJZHoKX+LFrqz0IwW+qXRQZDG3t2KKQky6Ow65lhvm4GEbnRhumDAQBarRYxMTEAgPTHVln08hiKDHL9u9DDnh0iIgoJotja21NQUMB1s0IMgx0KKaKohyAIEAQBWq3W180hIi8TRREVFRUoLS31dVPIixjsEBFRyHG0GCEFBwY7REQUchwtRkjBgcEOEREFhejoaLS0tECpVNo8RxAEZGRkIDs724stI19jsENEREFDJpNh6dKlVo8JggAAKCwsZL2dEMNgh4KWVquVkpFFUQ8AJpU3tm3bxhkZREFIpVIhMSHRYr9SqeS08xDFYIdCSnVVtfQzi4wRBa+OHTuabJeUlKC8vJyBTohisEMh4WJj61IROr1pT46hyBgDHqLgNmTIEA5dhTAGOxQSamtqre5nkTEiouDH5SIoaBkHLzq9zuZ/duMiYzk5OV5pGxF5h2El9NsKd0IQwtBFHiktL0Ghw697dhYtWoRbbrkFcrkcycnJGDNmDI4cOWJyTk5OjpSEang88sgjPmox+YvBCzag1//7GOmPrUL6Y6sgi+7c5jUsMkYUfAwroVfXN6OqrhFn6pt83STyAb/u2dm6dSumTp2KW265BS0tLXjmmWeQl5eHw4cPIzo6WjrvoYcewoIFC6TtTp06+aK55CfUajWOV9cjXJ7k1H9wFhkjCh5d5JEm26frG6EXbZxMQc+vg53PP//cZHvVqlVITk7Gvn37MGTIEGl/p06doFAovN088kM6nQ4zZ84E/vASAEDU66DTXrhyvOGCxTWCIECpVLLIGFEQMR+q6r9wI6rqGn3UGvI1vw52zNXWtiaZJiQkmOx///33sXr1aigUCowaNQrPPfec3d6dpqYmNDVd6cqsq6vzTIPJq7RaLWJiYgAA6Zf36bQXoHlzks1rWGSMKLQYFgMGgIaGBpNRAgpeARPs6PV6FBQUYNCgQbjhhhuk/ffeey+6deuGtLQ0HDx4EE899RSOHDlidyrxokWLMH/+fG80m/ycUqlEYWEha28QEQWxgAl2pk6dikOHDmH79u0m+6dMmSL93KtXL6SmpmLYsGE4duwYrr76aqvPNXfuXMyePVvarqurQ0ZGhmcaTn6rpKQEeXl57NEhIgpyfj0by2DatGn45JNPsHnzZrsLvAFAVlYWAODo0aM2z4mMjERsbKzJgwKfI3VyDN3XAIuMEYUiLhkTmvw62BFFEdOmTcNHH32ETZs2oXv37m1ec+DAAQCcWRNq1Go1evbsafcc40CHiEITl4wJTX4d7EydOhWrV6/GmjVrIJfLUVVVhaqqKly8eBEAcOzYMbz44ovYt28fjh8/jo8//hgTJ07EkCFD0Lt3bx+3nrxFrVYjPz8fGo3G7nlKpRKrV6/2UquIyJ8Y/m5wyZjQ5NfBzooVK1BbW4ucnBykpqZKj//85z8AgIiICHz99dfIy8vDddddh8cffxxjx47Fhg0bfNxy8hbDVHPDsg+2fPLJJygvL8fo0aO91DIi8hc6nQ41tTVWj3HJmNDg1wnKbf0By8jIwNatW73UGvJHpaWlqKysbPO8sLAwyGQyREdHt/n/ioiCS2lpKXQ6LhkTyvw62CFqi6NLPFRVVXm4JUTkr4w/JwxrZRnTNVxA1buzuGRMEGOwQwEtNTUViomvQxZjuvaV+VpYrLBNFLpaJ6z8AuDKWlm2z6NgxGCHApJxteSu096DYGehz4yMDOTl5XmraUTkZ7KzsyF7/zBa6k33y6I7QwhrLT+RkZHBJWOCGIMdCnhx8XGou2S5DhbQulzEG1wKgiikyWQyvD5Cgfz8fABX8kHTH1sl9fJwyZjgxmCHApLxrIkw4fKkwot1JutgZWRk4A0uBUFEAFQqFYqKijBjxgyLMhWJCYlQqe70UcvIG/x66jmRNeYFBM+eOwsAiIuLk/aVlJSgvLycgQ4RSVQqFQ4fPixtJyW29up07NjRV00iL2HPDgUUQwFBa9PHz184L/3MpSCIyBrjz4XIyAigudmHrSFvYc8OBQxHCwgSEREZY7BDAcPRAoJERM6o1FRCEARotVpfN4U8hMEOBQxHC36tWbMG0dHRHm4NEQUiQxV1URQhCPwTGCqYs0MBw9GCXywMRkTOMFRVHvTKFkRFRaKLPAobpg/2dbPIjRjsUMDIzs6GUqlEy7DHTSokG//MwmBE5KjWldAFqapyTTOA5iZcvNjo66aRmzHYoYBgXDE5PbqzzXLvLAxGRI5Qq9U4XXHC4ouTECbDufPnIAgCAKChoQEApM+fhoYGDpMHIAY7FBCMiwjGyuX4FaYVk2UyGTJTEqBSjfJRC4koUBhmdlaZTXgwrqhssG3bNvTv399kOy8vj1+qAgyzs8jvmRcRrKtvXeBGp70AzZuT8M87k/G/18dh+/MMdIiobc7M7Bw5ciSSkpJMtjMzM6FWqz3VPPIABjvkd+rq6iAIAgRBwLPPPov8/HyL8u7GWECQiJzh6MxOA71eb7Kt0WiQn5/PgCeAMNghv2Lei/Pyyy+3WUTQeIiLiKgt7Z2xafhMKigo4OdPgGCwQ37DsBSEvV4ca/bv3++hFhFRMDLM7DQkIbtCFEVUVFSgtLTUjS0jT2GwQ36hPUtBONslTUShTSaTYenSpQDQroAH4OdPoOBsLPILthIGFRNfhyyms8k+46miAIsIEpHzVCoVioqKMGPGDIveZEORQWO6hguoeneWxfPw8ycwMNghv2Dr25EsxnZNHYBFBInIdSqVCrm5uYiLiwMAhIW1DnYYigzaIwgClEolP38CBIexyC+09e1I1OvQUn/W5KHTXmARQSJql9jYWGmtrK5d4i0+Z0S9ZQKyYeiLnz+Bgz075BcMCYMajcZq3o6hpo5BRkYG3igshEql8mIriSiYbX9+FNRqtcnQlrVCg0qlEoX8/Ako7Nkhv+BMwmBJSQnKy8v5QUNEbqdSqXD48GFpOynRNNDh509gYrBDfsOQMJiWlmZxLCEhQfqZRQSJyJOMP18iIyNMjvHzJzBxGIv8ykpNMtIeXQmcPAkA6CBPggigU8dOLk1LJyJyVnR0tPR503/hRgBARGwSsl7+GsOW7gIAdJFHYsP0wT5rIzmHwQ75lTP1TThd3yyNkTO8ISJ/oBeBqrpGm8d1Oh1KS0tx6tQppKamSrO0zPexV8g3GOyQXwoTgGR5lLTdRR7pw9YQUagy/+w5Xd8Ivdm3MLVajZkzZ5rUCktMTAQAnDt3TtqnVCqxdOlS5vv4AIMd8ipr336sfdNJlkdh1zPDfNBCIqIrzIeq+i/ceLmHR8SWLVtQXFyMwsJCi+uMgxwDwwKiRUVFDHi8TBCZCIG6ujrExcWhtrYWsbGxvm5O0DL/9qOY+Doi4pIQHxePjh07ArjyrUkRy2CHiPyPIdgR9XrotOdNjtmqsmzMUIywvLycQ1pu4Ojfb/bskFcYFvk0jq1lMZ0hdOqM2ktA7SXbY+FERP7i4sWLAAQIYWFtVlm2xngB0ZycHLe3j6xjsEMe19Yin6JeBzTWIVWhANBaY4c5OkTkb3Q6HWqqTkAXEW2yXxbdGUKYc700XEDUuxjskMsczb8xLPJpvqinYUFPQ3XkdzZv5jcdIvJbpaWlOPF/Uy32G6osmy8gam9Yq7q6GjqdzuIz09HPVXIOiwqSS9RqNTIzMzF06FDce++9GDp0KDIzM6FWqy3ONXyDMSzqaXiYfxPiNx0i8mdtfUYZFhA1PIy/3JmbNWuWxWemM5+r5Bz27HiYI7UXBg4ciJ07d7b7HE8+t/E5P//8M1544QWLYanKykqMHTsWBQUF+LZLHhrRAQDQ0NAR6Y+tknpyRL0OOu2FK/9GDa0/Hz58GFu2bOE3GSLyS7YWLDZ8hhkYhrXMe3oM5xp6ezQaDcaOHYv58+fjwoULVmd1GZ/zm9/8xmOf6974++TLz/WgmY21fPlyLFmyBFVVVejTpw+WLVuGfv36OXStp2ZjOVp7QSaTQafTtfscTz63+TltsbZ4nkFL/VmTRT3NsRYFEfkjnU6HzMxMmwsWG7Tn888Rnvpc9+Q5nvpcd/Tvd1AEO//5z38wceJEvPXWW8jKykJhYSHWrVuHI0eOIDk5uc3rPRHsWJt9FMys5eMIYTKLXhyg7emZhoVAWYuCiPyN4bMdgNXPd8ue7QbU1Nba/Ex0ZLp6MPDU53pIBTtZWVm45ZZb8Le//Q0AoNfrkZGRgenTp+Ppp59u83p3BzuG6N+4RydQmActjnL3txjWoiAif2Wt1z4jIwOFhYUWf8j//e9/4957722zt8dZgRgkeeJzPWTq7DQ3N2Pfvn2YO3eutC8sLAy5ubkoKyuzek1TUxOampqk7draWgCt/2juYJh9BADyW8Yg9pYxbnleb3ClboQ5419c87FsRxlqUXz++efS+C8RkT/Izc3FwYMHsXPnTlRVVUGhUGDgwIGQyWQWf0cMf4DNPwuNP2td+dwNlydZ5AP5s7o961G/Z73bP9cN/95t9tuIAU6j0YgAxJ07d5rsnzNnjtivXz+r18ybN09E6xqTfPDBBx988MFHgD8qKirsxgoB37Pjirlz52L27NnStl6vx/nz55GYmCiNK7pDXV0dMjIyUFFREbTLUAT7PfL+Al+w3yPvL/AF+z168v5EUUR9fT3S0tLsnhfwwU5SUhJkMhmqq6tN9ldXV0OhUFi9JjIyEpGRphV64+PjPdVExMbGBuV/YGPBfo+8v8AX7PfI+wt8wX6Pnrq/uLi4Ns8J+KKCERERuOmmm7Bx40Zpn16vx8aNGzFgwAAftoyIiIj8QcD37ADA7Nmzcd999+Hmm29Gv379UFhYCK1Wi/vvv9/XTSMiIiIfC4pg509/+hPOnDmD559/HlVVVbjxxhvx+eefIyUlxaftioyMxLx58yyGzIJJsN8j7y/wBfs98v4CX7Dfoz/cX1DU2SEiIiKyJeBzdoiIiIjsYbBDREREQY3BDhEREQU1BjtEREQU1BjstMPLL7+MgQMHolOnTjaLEp44cQJ33HEHOnXqhOTkZMyZMwctLS12n/f8+fMYP348YmNjER8fj8mTJ6OhocEDd+CcLVu2QBAEq489e/bYvC4nJ8fi/EceecSLLXdOZmamRXsXL15s95rGxkZMnToViYmJiImJwdixYy0KXfqD48ePY/LkyejevTs6duyIq6++GvPmzUNzc7Pd6/z9PVy+fDkyMzMRFRWFrKwsfPPNN3bPX7duHa677jpERUWhV69eKCkp8VJLnbNo0SLccsstkMvlSE5OxpgxY3DkyBG716xatcrivYqKivJSi533wgsvWLT3uuuus3tNoLx/gPXPE0EQMHXqVKvn+/v7t23bNowaNQppaWkQBAHr1683OS6KIp5//nmkpqaiY8eOyM3Nxc8//9zm8zr7O+wsBjvt0NzcjLvvvhuPPvqo1eM6nQ533HEHmpubsXPnTrzzzjtYtWoVnn/+ebvPO378eHz//ff46quv8Mknn2Dbtm2YMmWKJ27BKQMHDsSpU6dMHg8++CC6d++Om2++2e61Dz30kMl1r776qpda7ZoFCxaYtHf69Ol2z581axY2bNiAdevWYevWrTh58qTF6sf+4Mcff4Rer8ff//53fP/993j99dfx1ltv4ZlnnmnzWn99D//zn/9g9uzZmDdvHvbv348+ffpg+PDhOH36tNXzd+7ciXHjxmHy5Mn49ttvMWbMGIwZMwaHDh3ycsvbtnXrVkydOhW7du3CV199hUuXLiEvLw9ardbudbGxsSbv1S+//OKlFrvmt7/9rUl7t2/fbvPcQHr/AGDPnj0m9/bVV18BAO6++26b1/jz+6fVatGnTx8sX77c6vFXX30Vb7zxBt566y3s3r0b0dHRGD58OBobG20+p7O/wy5xy2qcIW7lypViXFycxf6SkhIxLCxMrKqqkvatWLFCjI2NFZuamqw+1+HDh0UA4p49e6R9n332mSgIgqjRaNze9vZobm4Wu3TpIi5YsMDuebfeeqs4c+ZM7zTKDbp16ya+/vrrDp9fU1MjdujQQVy3bp2074cffhABiGVlZR5ooXu9+uqrYvfu3e2e48/vYb9+/cSpU6dK2zqdTkxLSxMXLVpk9fw//vGP4h133GGyLysrS3z44Yc92k53OH36tAhA3Lp1q81zbH0e+at58+aJffr0cfj8QH7/RFEUZ86cKV599dWiXq+3ejyQ3j8A4kcffSRt6/V6UaFQiEuWLJH21dTUiJGRkeK///1vm8/j7O+wK9iz40FlZWXo1auXSXHD4cOHo66uDt9//73Na+Lj4016SnJzcxEWFobdu3d7vM3O+Pjjj3Hu3DmHKlW///77SEpKwg033IC5c+fi119/9UILXbd48WIkJiaib9++WLJkid2hx3379uHSpUvIzc2V9l133XXo2rUrysrKvNHcdqmtrUVCQkKb5/nje9jc3Ix9+/aZ/NuHhYUhNzfX5r99WVmZyflA6+9loLxXANp8vxoaGtCtWzdkZGRg9OjRNj9v/MXPP/+MtLQ0XHXVVRg/fjxOnDhh89xAfv+am5uxevVqPPDAA3YXnQ6098+gvLwcVVVVJu9PXFwcsrKybL4/rvwOuyIoKij7q6qqKosqzobtqqoqm9ckJyeb7AsPD0dCQoLNa3zl7bffxvDhw6FUKu2ed++996Jbt25IS0vDwYMH8dRTT+HIkSNQq9VeaqlzZsyYgd/97ndISEjAzp07MXfuXJw6dQp//etfrZ5fVVWFiIgIi7ytlJQUv3vPzB09ehTLli3Da6+9Zvc8f30Pz549C51OZ/X37Mcff7R6ja3fS39/r/R6PQoKCjBo0CDccMMNNs/r0aMH/vWvf6F3796ora3Fa6+9hoEDB+L7779v83fVF7KysrBq1Sr06NEDp06dwvz585GdnY1Dhw5BLpdbnB+o7x8ArF+/HjU1NZg0aZLNcwLt/TNmeA+ceX9c+R12BYMdM08//TReeeUVu+f88MMPbSbQBRJX7rmyshJffPEFPvjggzaf3zjfqFevXkhNTcWwYcNw7NgxXH311a433AnO3OPs2bOlfb1790ZERAQefvhhLFq0yG/LubvyHmo0Gtx+++24++678dBDD9m91h/ew1A3depUHDp0yG4+CwAMGDDAZBHkgQMH4vrrr8ff//53vPjii55uptNGjBgh/dy7d29kZWWhW7du+OCDDzB58mQftsz93n77bYwYMQJpaWk2zwm09y9QMNgx8/jjj9uNugHgqquucui5FAqFRUa5YYaOQqGweY15UlZLSwvOnz9v85r2cuWeV65cicTERPzhD39w+vWysrIAtPYqeOsPZXve16ysLLS0tOD48ePo0aOHxXGFQoHm5mbU1NSY9O5UV1d77D0z5+z9nTx5EkOHDsXAgQPxj3/8w+nX88V7aE1SUhJkMpnFzDd7//YKhcKp8/3BtGnTpMkKzn6779ChA/r27YujR496qHXuFR8fj2uvvdZmewPx/QOAX375BV9//bXTvaGB9P4Z3oPq6mqkpqZK+6urq3HjjTdavcaV32GXuC37J4S1laBcXV0t7fv73/8uxsbGio2NjVafy5CgvHfvXmnfF1984VcJynq9Xuzevbv4+OOPu3T99u3bRQDif//7Xze3zDNWr14thoWFiefPn7d63JCgXFRUJO378ccf/TZBubKyUvzNb34j3nPPPWJLS4tLz+FP72G/fv3EadOmSds6nU5MT0+3m6B85513muwbMGCAXya46vV6cerUqWJaWpr4008/ufQcLS0tYo8ePcRZs2a5uXWeUV9fL3bu3FlcunSp1eOB9P4ZmzdvnqhQKMRLly45dZ0/v3+wkaD82muvSftqa2sdSlB25nfYpba67ZlC0C+//CJ+++234vz588WYmBjx22+/Fb/99luxvr5eFMXW/6Q33HCDmJeXJx44cED8/PPPxS5duohz586VnmP37t1ijx49xMrKSmnf7bffLvbt21fcvXu3uH37dvE3v/mNOG7cOK/fny1ff/21CED84YcfLI5VVlaKPXr0EHfv3i2KoigePXpUXLBggbh3716xvLxcLC4uFq+66ipxyJAh3m62Q3bu3Cm+/vrr4oEDB8Rjx46Jq1evFrt06SJOnDhROsf8HkVRFB955BGxa9eu4qZNm8S9e/eKAwYMEAcMGOCLW7CrsrJSvOaaa8Rhw4aJlZWV4qlTp6SH8TmB9B6uXbtWjIyMFFetWiUePnxYnDJlihgfHy/Ngvzzn/8sPv3009L5O3bsEMPDw8XXXntN/OGHH8R58+aJHTp0EL/77jtf3YJNjz76qBgXFydu2bLF5L369ddfpXPM72/+/PniF198IR47dkzct2+feM8994hRUVHi999/74tbaNPjjz8ubtmyRSwvLxd37Ngh5ubmiklJSeLp06dFUQzs989Ap9OJXbt2FZ966imLY4H2/tXX10t/6wCIf/3rX8Vvv/1W/OWXX0RRFMXFixeL8fHxYnFxsXjw4EFx9OjRYvfu3cWLFy9Kz3HbbbeJy5Ytk7bb+h12BwY77XDfffeJACwemzdvls45fvy4OGLECLFjx45iUlKS+Pjjj5tE9ps3bxYBiOXl5dK+c+fOiePGjRNjYmLE2NhY8f7775cCKH8wbtw4ceDAgVaPlZeXm/wbnDhxQhwyZIiYkJAgRkZGitdcc404Z84csba21ostdty+ffvErKwsMS4uToyKihKvv/56ceHChSY9ceb3KIqiePHiRfGxxx4TO3fuLHbq1Em86667TAIIf7Fy5Uqr/2eNO3kD8T1ctmyZ2LVrVzEiIkLs16+fuGvXLunYrbfeKt53330m53/wwQfitddeK0ZERIi//e1vxU8//dTLLXaMrfdq5cqV0jnm91dQUCD9W6SkpIgjR44U9+/f7/3GO+hPf/qTmJqaKkZERIjp6enin/70J/Ho0aPS8UB+/wy++OILEYB45MgRi2OB9v4Z/maZPwz3oNfrxeeee05MSUkRIyMjxWHDhlncd7du3cR58+aZ7LP3O+wOgiiKovsGxYiIiIj8C+vsEBERUVBjsENERERBjcEOERERBTUGO0RERBTUGOwQERFRUGOwQ0REREGNwQ4REREFNQY7REREFNQY7BAREVFQY7BDREREQY3BDhEREQU1BjtEFHTOnDkDhUKBhQsXSvt27tyJiIgIbNy40YctIyJf4EKgRBSUSkpKMGbMGOzcuRM9evTAjTfeiNGjR+Ovf/2rr5tGRF7GYIeIgtbUqVPx9ddf4+abb8Z3332HPXv2IDIy0tfNIiIvY7BDREHr4sWLuOGGG1BRUYF9+/ahV69evm4SEfkAc3aIKGgdO3YMJ0+ehF6vx/Hjx33dHCLyEfbsEFFQam5uRr9+/XDjjTeiR48eKCwsxHfffYfk5GRfN42IvIzBDhEFpTlz5qCoqAj//e9/ERMTg1tvvRVxcXH45JNPfN00IvIyDmMRUdDZsmULCgsL8d577yE2NhZhYWF47733UFpaihUrVvi6eUTkZezZISIioqDGnh0iIiIKagx2iIiIKKgx2CEiIqKgxmCHiIiIghqDHSIiIgpqDHaIiIgoqDHYISIioqDGYIeIiIiCGoMdIiIiCmoMdoiIiCioMdghIiKioPb/AXXeJYyavdOsAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# construct model representation for comparison with data histogram\n", "a, b, sigma = m.values\n", "\n", "# get expected content per bin from cdf, sum over the individual cdfs\n", "v = np.diff(np.sum(norm.cdf(ax_x.edges[:,np.newaxis],\n", " mu(y, a, b), sigma), axis=1))\n", "\n", "plt.stairs(v, ax_x.edges, label=\"model\", zorder=5, lw=2)\n", "plt.errorbar(ax_x.centers, h1.values(), h1.variances() ** 0.5,\n", " fmt=\"ok\", label=\"data\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"events\")\n", "plt.legend(frameon=False);" ] }, { "cell_type": "markdown", "id": "integrated-listening", "metadata": {}, "source": [ "## Fit without conditional variable\n", "\n", "We can also ignore the dependence of $x$ and $y$ and just fit the total $x$ distribution with a model built from the distribution of $y$ values. This also works in this case, but information is lost and therefore the parameter uncertainties become larger than in the previous case.\n", "\n", "On top of that, the calculation is much slower, because building the pdf is more expensive. We parallelise the computation with numba." ] }, { "cell_type": "code", "execution_count": 6, "id": "protecting-monte", "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 = -8.774e+04 Nfcn = 95
EDM = 1.33e-06 (Goal: 0.0002) time = 5.2 sec
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", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 a 0.002 0.029
1 b 0.500 0.005
2 sigma 0.98 0.04 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", "
a b sigma
a 0.000839 4.35e-06 (0.030) -3.16e-05 (-0.027)
b 4.35e-06 (0.030) 2.43e-05 -0.000141 (-0.718)
sigma -3.16e-05 (-0.027) -0.000141 (-0.718) 0.0016
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = -8.774e+04 │ Nfcn = 95 │\n", "│ EDM = 1.33e-06 (Goal: 0.0002) │ time = 5.2 sec │\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 │ a │ 0.002 │ 0.029 │ │ │ │ │ │\n", "│ 1 │ b │ 0.500 │ 0.005 │ │ │ │ │ │\n", "│ 2 │ sigma │ 0.98 │ 0.04 │ │ │ 0 │ │ │\n", "└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌───────┬───────────────────────────────┐\n", "│ │ a b sigma │\n", "├───────┼───────────────────────────────┤\n", "│ a │ 0.000839 4.35e-06 -3.16e-05 │\n", "│ b │ 4.35e-06 2.43e-05 -0.000141 │\n", "│ sigma │ -3.16e-05 -0.000141 0.0016 │\n", "└───────┴───────────────────────────────┘" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nb.config.THREADING_LAYER = 'workqueue'\n", "\n", "\n", "@nb.njit(parallel=True, fastmath=True)\n", "def model(x, a, b, sigma):\n", " mu = a + b * y\n", " total = np.zeros_like(x)\n", " for i in nb.prange(len(mu)):\n", " total += norm_nb.pdf(x, mu[i], sigma)\n", " return total\n", "\n", "\n", "nll = UnbinnedNLL(x, model)\n", "m2 = Minuit(nll, 0.0, 0.0, 2.0)\n", "m2.limits[\"sigma\"] = (0, None)\n", "m2.migrad()" ] }, { "cell_type": "code", "execution_count": 7, "id": "julian-border", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAADTCAYAAACBfoy1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA24ElEQVR4nO3deVhTZ9oG8DtBCJuAyhIQKu6IG64UqVulRe0iLlWprcuoHbValaqV1t0ZdbRarGs/p4pttW5ltFMdWmXEDdzFWqsWFVwBV1YVkLzfH5QzRggGJeQkuX/XlUvOyXNO3pP3zWOenE0hhBAgIiIiIiKSGaWxG0BERERERFQWFitERERERCRLLFaIiIiIiEiWWKwQEREREZEssVghIiIiIiJZYrFCRERERESyxGKFiIiIiIhkicUKERERERHJEosVIiIiIiKSJRYrRERkkmbNmgWFQoE7d+4YuylEZABDhw6Fr6+vsZtBRsZihYiIiIiIZKmasRtARERERPS0NWvWQKPRGLsZZGQsVoiIiIhIdqytrY3dBJIBHgZGsnHlyhWMGTMGjRs3hp2dHWrVqoV33nkHqampxm4aEcnYnTt30L9/fzg5OaFWrVoYP348Hj16ZOxmEdEz5OTkYMKECfD19YVKpYK7uztee+01nDx5EkDZ56zcvXsX77//PpycnODi4oIhQ4bg9OnTUCgUiI6OluKGDh0KR0dHXL16FW+++SYcHR1Ru3ZtrFixAgBw5swZvPrqq3BwcECdOnWwceNGrde5d+8eJk2ahObNm8PR0RFOTk7o0aMHTp8+bdD3hErjnhWSjWPHjiEhIQEDBw6Et7c3UlNTsWrVKnTp0gW///477O3tjd1EIpKh/v37w9fXF/Pnz8fhw4fx5Zdf4v79+/jmm2+M3TQiKseoUaOwbds2jB07Fv7+/rh79y4OHjyIc+fOoXXr1qXiNRoN3nrrLRw9ehSjR4+Gn58fduzYgSFDhpS5/qKiIvTo0QOdOnXCwoULsWHDBowdOxYODg747LPPMGjQIPTp0werV6/G4MGDERQUhLp16wIALl++jO3bt+Odd95B3bp1kZGRga+++gqdO3fG77//Di8vL4O+N/QEQSQTDx48KDUvMTFRABDffPONEVpERHI2c+ZMAUC8/fbbWvPHjBkjAIjTp08bqWVEpA9nZ2fx4Ycf6nx+yJAhok6dOtL0Dz/8IACIqKgoaV5RUZF49dVXBQCxbt06rWUBiHnz5knz7t+/L+zs7IRCoRCbNm2S5p8/f14AEDNnzpTmPXr0SBQVFWm1JyUlRahUKjFnzpzn2Fp6XjwMjGTDzs5O+ruwsBB3795FgwYN4OLiIu0SJiJ62ocffqg1PW7cOADArl27jNEcItKTi4sLjhw5gps3b+oVHxsbC2tra4wcOVKap1QqS+WAJ40YMULr9Ro3bgwHBwf0799fmt+4cWO4uLjg8uXL0jyVSgWlsvhrclFREe7evQtHR0c0btyY30mqGIsVko2HDx9ixowZ8PHxgUqlgqurK9zc3JCZmYmsrCxjN4+IZKphw4Za0/Xr14dSqeT5bkQyt3DhQvz222/w8fFB+/btMWvWLK2C4WlXrlyBp6dnqcPCGzRoUGa8ra0t3NzctOY5OzvD29sbCoWi1Pz79+9L0xqNBl988QUaNmyo9Z3k119/5XeSKsZihWRj3Lhx+Pvf/47+/ftjy5Yt+OWXX7B7927UqlWLly4kIr09/SWEiOSpf//+uHz5MpYtWwYvLy8sWrQITZs2xX/+859KWb+VlVWF5gshpL/nzZuHiIgIdOrUCd999x1+/vln7N69G02bNuV3kirGE+xJNrZt24YhQ4Zg8eLF0rxHjx4hMzPTeI0iItlLTk6WTooFgIsXL0Kj0fDO10QmwNPTE2PGjMGYMWNw69YttG7dGn//+9/Ro0ePUrF16tTB3r178eDBA629KxcvXqz0dm3btg1du3bF119/rTU/MzMTrq6ulf56pBv3rJBsWFlZaf2qAQDLli1DUVGRkVpERKag5FKkJZYtWwYAZX7ZISJ5KCoqKnU4lbu7O7y8vJCfn1/mMqGhoSgsLMSaNWukeRqNplQOqAxlfSfZunUrbty4UemvReXjnhWSjTfffBPffvstnJ2d4e/vj8TEROzZswe1atUydtOISMZSUlLw9ttvo3v37khMTMR3332Hd999Fy1btjR204hIh5ycHHh7e6Nfv35o2bIlHB0dsWfPHhw7dkzrCIsnhYWFoX379vj4449x8eJF+Pn54ccff8S9e/cAVO4hoG+++SbmzJmDYcOGoUOHDjhz5gw2bNiAevXqVdprkH5YrJBsLF26FFZWVtiwYQMePXqE4OBg7NmzB6GhocZuGhHJ2ObNmzFjxgxMnToV1apVw9ixY7Fo0SJjN4uIymFvb48xY8bgl19+QUxMDDQaDRo0aICVK1di9OjRZS5jZWWFnTt3Yvz48Vi/fj2USiV69+6NmTNnIjg4GLa2tpXWvk8//RR5eXnYuHEjNm/ejNatW2Pnzp2YOnVqpb0G6Uchnt7HRURERERkIrZv347evXvj4MGDCA4ONnZzqJKxWCEiIiIik/Dw4UOt+7IVFRXh9ddfx/Hjx5Genq71HJkHHgZGRERERCZh3LhxePjwIYKCgpCfn4+YmBgkJCRg3rx5LFTMFPesEBEREZFJ2LhxIxYvXoyLFy/i0aNHaNCgAUaPHo2xY8cau2lkICxWiIiIiIhIlnifFSIiIiIikiUWK0REREREJEsWeYK9RqPBzZs3Ub169Uq9gRCRpRFCICcnB15eXlAqzfe3D+YMosphDjmD+YCocuibDyyyWLl58yZ8fHyM3Qwis3Ht2jV4e3sbuxkGw5xBVLlMOWcwHxBVrmflA4ssVqpXrw6g+M1xcnIycmuITFd2djZ8fHykz5S5Ys4gqhzmkDOYD4gqh775wCKLlZLdtk5OTkw0RJXA3A+FYM4gqlymnDOYD4gq17PygWkeMEpERERERGaPxQoREREREcmSwYuVFStWwNfXF7a2tggMDMTRo0fLjd+6dSv8/Pxga2uL5s2bY9euXVrPDx06FAqFQuvRvXt3Q24CEREREREZgUGLlc2bNyMiIgIzZ87EyZMn0bJlS4SGhuLWrVtlxickJCA8PBzDhw/HqVOnEBYWhrCwMPz2229acd27d0daWpr0+P777w25GUREREREZAQGLVaWLFmCkSNHYtiwYfD398fq1athb2+PtWvXlhm/dOlSdO/eHZMnT0aTJk0wd+5ctG7dGsuXL9eKU6lUUKvV0qNGjRqG3AwiIiIiIjICgxUrBQUFOHHiBEJCQv73YkolQkJCkJiYWOYyiYmJWvEAEBoaWio+Pj4e7u7uaNy4MUaPHo27d++W25b8/HxkZ2drPYiIiIiISN4MVqzcuXMHRUVF8PDw0Jrv4eGB9PT0MpdJT09/Znz37t3xzTffIC4uDv/4xz+wb98+9OjRA0VFRTrbMn/+fDg7O0sP3syJiIiIiEj+TO4+KwMHDpT+bt68OVq0aIH69esjPj4e3bp1K3OZyMhIRERESNMlN6EhIiIiIiL5MtieFVdXV1hZWSEjI0NrfkZGBtRqdZnLqNXqCsUDQL169eDq6oqLFy/qjFGpVNLNm3gTJyIiIiIi02CwYsXGxgZt2rRBXFycNE+j0SAuLg5BQUFlLhMUFKQVDwC7d+/WGQ8A169fx927d+Hp6Vk5DSciIiIiIlkw6NXAIiIisGbNGqxfvx7nzp3D6NGjkZeXh2HDhgEABg8ejMjISCl+/PjxiI2NxeLFi3H+/HnMmjULx48fx9ixYwEAubm5mDx5Mg4fPozU1FTExcWhV69eaNCgAUJDQw25KUREREREVMUMes7KgAEDcPv2bcyYMQPp6ekICAhAbGysdBL91atXoVT+r17q0KEDNm7ciGnTpuHTTz9Fw4YNsX37djRr1gwAYGVlhV9//RXr169HZmYmvLy88Prrr2Pu3LlQqVSG3BQiIiIiIqpiCiGEMHYjqlp2djacnZ2RlZXF81eIXoClfJYsZTuJDM0cPkvmsA1EcqDvZ8mgh4ERERERERE9LxYrREREREQkSyxWiIiIiIhIllisEBERERGRLLFYISIiItnZv38/3nrrLXh5eUGhUGD79u3PXCY+Ph6tW7eGSqVCgwYNEB0drfX8rFmzoFAotB5+fn6G2QAiqhQsVoiIiEh28vLy0LJlS6xYsUKv+JSUFLzxxhvo2rUrkpKSMGHCBIwYMQI///yzVlzTpk2RlpYmPQ4ePGiI5hNRJTHofVaIiIiInkePHj3Qo0cPveNXr16NunXrYvHixQCAJk2a4ODBg/jiiy+0bhxdrVo1qNXqSm8vERkG96wQERGRyUtMTERISIjWvNDQUCQmJmrNS05OhpeXF+rVq4dBgwbh6tWr5a43Pz8f2dnZWg8iqjosVoiIiMjkpaenw8PDQ2ueh4cHsrOz8fDhQwBAYGAgoqOjERsbi1WrViElJQUdO3ZETk6OzvXOnz8fzs7O0sPHx8eg20FE2lisEBERkUXo0aMH3nnnHbRo0QKhoaHYtWsXMjMzsWXLFp3LREZGIisrS3pcu3atCltMRDxnhYiIiEyeWq1GRkaG1ryMjAw4OTnBzs6uzGVcXFzQqFEjXLx4Ued6VSoVVCpVpbaViPTHPStERERk8oKCghAXF6c1b/fu3QgKCtK5TG5uLi5dugRPT09DN4+InhOLFSKqMitWrICvry9sbW0RGBiIo0eP6oyNjo4udT8EW1tbrRghBGbMmAFPT0/Y2dkhJCQEycnJWjG+vr6l1rNgwQKDbB8RVZ7c3FwkJSUhKSkJQPGliZOSkqQT4iMjIzF48GApftSoUbh8+TKmTJmC8+fPY+XKldiyZQsmTpwoxUyaNAn79u1DamoqEhIS0Lt3b1hZWSE8PLxKt42I9MdihYiqxObNmxEREYGZM2fi5MmTaNmyJUJDQ3Hr1i2dyzg5OWndD+HKlStazy9cuBBffvklVq9ejSNHjsDBwQGhoaF49OiRVtycOXO01jNu3DiDbCMRVZ7jx4+jVatWaNWqFQAgIiICrVq1wowZMwAAaWlpWlfyqlu3Lnbu3Indu3ejZcuWWLx4Mf75z39qXbb4+vXrCA8PR+PGjdG/f3/UqlULhw8fhpubW9VuHBHpjeesEFWFvDzA0bH479xcwMHBuO0xgiVLlmDkyJEYNmwYgOJ7IuzcuRNr167F1KlTy1xGoVDovB+CEAJRUVGYNm0aevXqBQD45ptv4OHhge3bt2PgwIFSbPXq1XlfBTItzBno0qULhBA6n3/67vQly5w6dUrnMps2baqMphFRFeKeFSIyuIKCApw4cULrHghKpRIhISGl7oHwpNzcXNSpUwc+Pj7o1asXzp49Kz2XkpKC9PR0rXU6OzsjMDCw1DoXLFiAWrVqoVWrVli0aBEeP35cbnt5XwUiIiJ5YLFCRAZ3584dFBUVlXkPhPT09DKXady4MdauXYsdO3bgu+++g0ajQYcOHXD9+nUAkJZ71jo/+ugjbNq0CXv37sVf//pXzJs3D1OmTCm3vbyvAhERkTzwMDAikqWgoCCtq/h06NABTZo0wVdffYW5c+fqvZ6IiAjp7xYtWsDGxgZ//etfMX/+fJ2XI42MjNRaLjs7mwULERGREXDPChEZnKurK6ysrMq8B4K+55JYW1ujVatW0v0QSpar6DoDAwPx+PFjpKam6oxRqVRwcnLSehAREVHVY7FCRAZnY2ODNm3aaN0DQaPRIC4urtx7IDypqKgIZ86cke6HULduXajVaq11Zmdn48iRI+WuMykpCUqlEu7u7s+5NURERFRVeBgYEVWJiIgIDBkyBG3btkX79u0RFRWFvLw86epggwcPRu3atTF//nwAxZcbfvnll9GgQQNkZmZi0aJFuHLlCkaMGAGg+EphEyZMwN/+9jc0bNgQdevWxfTp0+Hl5YWwsDAAQGJiIo4cOYKuXbuievXqSExMxMSJE/Hee++hRo0aRnkfiIiISH8sVoioSgwYMAC3b9/GjBkzkJ6ejoCAAMTGxkonyF+9ehVK5f929t6/fx8jR45Eeno6atSogTZt2iAhIQH+/v5SzJQpU5CXl4cPPvgAmZmZeOWVVxAbGyvdPFKlUmHTpk2YNWsW8vPzUbduXUycOFHrfBQiIiKSL4Uo7yLmZio7OxvOzs7IysrisehUNcz0ngmW8lmylO0kGWHOkC1z2AYiOdD3s8RzVoiIiIiISJZYrBARERERkSyxWCEiIiIiIllisUJERERERLLEYoWIiIiIiGSJxQoREREREckSixUiIiIiIpIlFitERERERCRLLFaIiIiIiEiWWKwQEREREclVXh6gUBQ/8vKM3Zoqx2KFiIiIiIhkicUKERERERHJEosVIiIiIiKSJRYrREREREQkSyxWiIiIiIhIllisEBFZkLy8PCgUCigUCuRZ4FVliIjItLBYISIiIiIiWWKxQkREREREssRihagqFBX97+/9+7WniYiexpxBRASAxQqR4cXEAP7+/5vu2RPw9S2eT0T0NOYMIiIJixUiQ4qJAfr1A27c0J5/40bxfH75IKInMWcQEWlhsUJkKEVFwPjxgBClnyuZN2ECD+8gomLMGUREpRi8WFmxYgV8fX1ha2uLwMBAHD16tNz4rVu3ws/PD7a2tmjevDl27dql9bwQAjNmzICnpyfs7OwQEhKC5ORkQ24C0fM5cAC4fl3380IA164Vx1mIiuSD6Oho6RK7JQ9bW1utGH3ywb179zBo0CA4OTnBxcUFw4cPR25urkG2j+iFMGdo2b9/P9566y14eXlBoVBg+/btz1wmPj4erVu3hkqlQoMGDRAdHV0qpqLfS4jIuKoZcuWbN29GREQEVq9ejcDAQERFRSE0NBQXLlyAu7t7qfiEhASEh4dj/vz5ePPNN7Fx40aEhYXh5MmTaNasGQBg4cKF+PLLL7F+/XrUrVsX06dPR2hoKH7//fdSX2SepaCgAAUFBaXmK5VKVKtWTStOF4VCAWtr6+eKLSwshCjrFzQDxgKAjY3Nc8U+fvwYGo2mUmKtra2hUCgMGltUVISicn6BrEhstWrVoFQqKxablqYzRsufcRVpg0ajwePHj3XGWllZwcrKyuCx5Y33p1U0HwCAk5MTLly4IE2X9FcJffLBoEGDkJaWht27d6OwsBDDhg3DBx98gI0bN+rd9hLmkDPKej3mDOYMOeaMvLw8tGzZEn/5y1/Qp0+fZ8anpKTgjTfewKhRo7BhwwbExcVhxIgR8PT0RGhoKIDny0O6mEM+qGgswHzwPLHPlQ+ejC0oQMm7U1BQAPzZX8+7XlPLBwpR3ih7QYGBgWjXrh2WL18uNczHxwfjxo3D1KlTS8UPGDAAeXl5+Omnn6R5L7/8MgICArB69WoIIeDl5YWPP/4YkyZNAgBkZWXBw8MD0dHRGDhwYJntyM/PR35+vjSdnZ0NHx8fTJ06tcwCp2HDhnj33Xel6Xnz5qGwsLDMddepUwdDhw6VphctWoQHDx6UGevl5YWRI0dK01FRUcjKyioz1s3NDWPGjJGmV65cidu3b5cZa2trK72fubm52LhxI27evFlmrL29PSZPnixNR0dH48qVK2XGWltb49NPP5WmN27cWO5erJkzZ0p/b926Fb///rvO2MjISCkxbd++HadPn9YZO2nSJDg4OAAAdu7ciePHj+uMHT9+PFxcXAAAv/zyCxITE3XGjh49WvrPKT4+Hvv27dMZO2LECNSuXRsAcOjQIezZs0dn7JAhQ+Dr6wvExwNdu+qMk+zdC3TpgqSkJOzYsUNnWL9+/dC0aVMAwNmzZ7Ft2zadsb169UJAQAAA4I8//sD333+vM7ZHjx5o3749ACA1NRXr16/XGRsSEoLg4GAAwI0bN7B8+XIsWLAAWVlZcHJy0rkcUPF8EB0djQkTJiAzM7PM9emTD86dOwd/f38cO3YMbdu2BQDExsaiZ8+euH79Ory8vMpctznnjIKCAsybNw9Acb5wcHDAmjVrmDPAnCG3nPEkhUKBf/3rXwgLC9MZ88knn2Dnzp347bffpHkDBw5EZmYmYmNjAVQ8DwHmnQ8AwNnZGRMmTJCmmQ9cAMggHwA4evQo/vOf/8C6oACf/pm35336KQr/3Pbw8HA0atQIAMw6HxjsMLCCggKcOHECISEh/3sxpRIhISE6Oz8xMVErHgBCQ0Ol+JSUFKSnp2vFODs7IzAwsNwBNX/+fDg7O0sPHx+fF9k0Iv107Ih8d3fo+jVAAChUq4GOHauyVUbxPPkAKP4yXadOHfj4+KBXr144e/as9Jw++SAxMREuLi5SoQIUJ0ulUokjR47ofF1zzhlP/qq4f//+cn+JoyrGnPFCnvUd4nnzkDnnAyJTYLA9Kzdv3kTt2rWRkJCAoKAgaf6UKVOwb9++Mr8o2NjYYP369QgPD5fmrVy5ErNnz0ZGRgYSEhIQHByMmzdvwtPTU4rp378/FAoFNm/eXGZbdP0qcvv27TIrOVPbhfvgwQPUqFEDQPGXOxsbG+7ChUx24W7bBmX//gAAxRPvs/jzdcWWLVD261fh9cplF+69e/fg5ub2zF9FnicfJCYmIjk5GS1atEBWVhY+//xz7N+/H2fPnoW3t7de+WDevHlYv3691qFkAODu7o7Zs2dj9OjRZbbXXHPG9u3bMXHiRK1fTb29vbF48eJyf7Fmzqh4LHNG2bH65oyn6bNnpVGjRhg2bBgiIyOlebt27cIbb7yBBw8e4P79+xXOQ4D55oPniQWYD54n9oUPA8vLg03NmgCAgnv3gD/3FJn6YWD65gODnrMiFyqVCiqVqtR8GxsbrQ+HLvrEPE/sk8nhRWKf3r1cWet92pPJ1xRin/xQGC22Xz9g2zbgo4+0LkWq8PYGoqKgeOI47IqsV6lU6j3W5BD7PIKCgrS+UHTo0AFNmjTBV199hblz5xrsdQHzzBkxMTEYOHBgqS8WN27cwMCBA7Ft2za9zgtgzjBwLHOG7JhjPjBGrBw+4yaXD0pin/ieZ2NjA5QxTsw5HxjsMDBXV1dYWVkhIyNDa35GRgbUanWZy6jV6nLjS/6tyDqJjK5PH+DJY2937QJSUornW4jnyQdPs7a2RqtWrXDx4kUA+uUDtVqNW7duaT3/+PFj3Lt3z6JyRlFREcaPH1/mL6Al8yZMmMBDwuSCOeO56PoO4eTkBDs7u0rJQ0RU9QxWrNjY2KBNmzaIi4uT5mk0GsTFxWn9WvqkoKAgrXgA2L17txRft25dqNVqrZjs7GwcOXJE5zqJZOHJXzs6ddKetgDPkw+eVlRUhDNnzkiHfOmTD4KCgpCZmYkTJ05IMf/973+h0WgQGBhYGZtmEg4cOIDr5VwSVwiBa9eu4YCFXBLXJFh4zngez/oOURl5iIiMQBjQpk2bhEqlEtHR0eL3338XH3zwgXBxcRHp6elCCCHef/99MXXqVCn+0KFDolq1auLzzz8X586dEzNnzhTW1tbizJkzUsyCBQuEi4uL2LFjh/j1119Fr169RN26dcXDhw/1bldWVpYAILKysipvY40oNzdXoPjcS5Gbm2vs5lBZcnOFKL5LQvHfZqIin6WK5oPZs2eLn3/+WVy6dEmcOHFCDBw4UNja2oqzZ89KMfrkg+7du4tWrVqJI0eOiIMHD4qGDRuK8PBwg22nHG3cuFHKEeU9Nm7caOymUgnmDJGTkyNOnTolTp06JQCIJUuWiFOnTokrV64IIYSYOnWqeP/996X4y5cvC3t7ezF58mRx7tw5sWLFCmFlZSViY2OlmGflocreBrnj9wcTYeH5wKDFihBCLFu2TLz00kvCxsZGtG/fXhw+fFh6rnPnzmLIkCFa8Vu2bBGNGjUSNjY2omnTpmLnzp1az2s0GjF9+nTh4eEhVCqV6Natm7hw4UKF2mROiUYIJhuTYOGJpkRF8sGECROkWA8PD9GzZ09x8uRJrfXpkw/u3r0rwsPDhaOjo3BychLDhg0TOTk5Bt1Oudm7d69excrevXuN3VQqwZyhc9yW5IkhQ4aIzp07l1omICBA2NjYiHr16ol169aVWm95eaiyt0HuSrYFgNi1a5d4/PixsZtEZbHwfGDQ+6zIVXZ2NpydnSt8NRK5KtkeoPjKJ6+//rreJ1lRFcnLAxwdi//OzZWu5GHqzO2zpIupb2dRURF8fX1x48aNMs9bUSgU8Pb2RkpKCnOHXDBnyJY5bANQfNGNjz76CDeeuJCDt7c3li5dqtfFNqgKWXg+MNg5K1Q1YmJi4O/vL0337NkTvr6+iImJMWKriEhOrKyssHTpUgCQLrVZomQ6KiqKhQqRhYiJiUG/fv20ChWg+OqA/fr143cIkhUWKyaMyYaI9NWnTx9s27YNXl5eWvO9vb31vmwxEZk+Xh2QTA2LFRPFZENEFdWnTx/8/sQlcXft2oWUlBQWKkQWhFcHJFPDYsVEMdkQ0fN48lCvTp068dAvIguTlpZWqXFUBZ784Xn/fu1pC8BixUQx2RAREVFFldyrqrLiyMBiYoAnzk1Gz56Ar2/xfAvBYsVEMdkQERFRRXXs2BHe3t6lLrZRQqFQwMfHBx07dqzillEpMTFAv37AU+cm48aN4vkWUrCwWDFRTDZERERUUbw6oIkoKgLGjy++u8rTSuZNmGARh4SxWDFRTDZERET0PHh1QBNw4ABQzrnJEAK4dq04zsyxWDFhTDZERET0PHh1QJnT95xjCzg3uZqxG0Avpk+fPggJCeEd7ImIiKhCeHVAGdP3nGMLODeZe1bMAJMNERERVZSDgwOEEBBCwMHBwdjNoSd17Ah4ewM6zk2GQgH4+BTHmTkWK0REFoRfToiITICVFfDnucmlCpaS6aio4jgzx2KFiIiIiEhu+vQBtm0Dnjo3Gd7exfMt5PwinrNCRERERCRHffoAISHAn+cmY9cu4PXXLWKPSgnuWSEiIiIikqsnC5NOnSyqUAFYrBARERERkUyxWCEiIiIiIlniOStmoOTqPkRERERE5oR7VoiIiIiISJZYrBARERERkSyxWCEiIiIiIllisUJEVWbFihXw9fWFra0tAgMDcfToUb2W27RpExQKBcLCwrTmZ2RkYOjQofDy8oK9vT26d++O5ORkrZguXbpAoVBoPUaNGlVZm0REREQGxGKFqCo4OABCFD8cHIzdGqPYvHkzIiIiMHPmTJw8eRItW7ZEaGgobt26Ve5yqampmDRpEjp27Kg1XwiBsLAwXL58GTt27MCpU6dQp04dhISEIC8vTyt25MiRSEtLkx4LFy6s9O0jqlTMGUREAFisEFEVWbJkCUaOHIlhw4bB398fq1evhr29PdauXatzmaKiIgwaNAizZ89GvXr1tJ5LTk7G4cOHsWrVKrRr1w6NGzfGqlWr8PDhQ3z//fdasfb29lCr1dLDycnJINtIRERElYvFChEZXEFBAU6cOIGQkBBpnlKpREhICBITE3UuN2fOHLi7u2P48OGlnsvPzwcA2Nraaq1TpVLh4MGDWrEbNmyAq6srmjVrhsjISDx48KDc9ubn5yM7O1vrQURERFWP91khIoO7c+cOioqK4OHhoTXfw8MD58+fL3OZgwcP4uuvv0ZSUlKZz/v5+eGll15CZGQkvvrqKzg4OOCLL77A9evXkZaWJsW9++67qFOnDry8vPDrr7/ik08+wYULFxATE6OzvfPnz8fs2bMrvqFERERUqVisEJHs5OTk4P3338eaNWvg6upaZoy1tTViYmIwfPhw1KxZE1ZWVggJCUGPHj20bpL6wQcfSH83b94cnp6e6NatGy5duoT69euXue7IyEhERERI09nZ2fDx8amkrSMiIiJ9sVghIoNzdXWFlZUVMjIytOZnZGRArVaXir906RJSU1Px1ltvSfM0Gg0AoFq1arhw4QLq16+PNm3aICkpCVlZWSgoKICbmxsCAwPRtm1bnW0JDAwEAFy8eFFnsaJSqaBSqSq8nURERFS5eM4KERmcjY0N2rRpg7i4OGmeRqNBXFwcgoKCSsX7+fnhzJkzSEpKkh5vv/02unbtiqSkpFJ7OZydneHm5obk5GQcP34cvXr10tmWksPKPD09K2fjiIiIyGBYrBBRlYiIiMCaNWuwfv16nDt3DqNHj0ZeXh6GDRsGABg8eDAiIyMBFJ8036xZM62Hi4sLqlevjmbNmsHGxgYAsHXrVsTHx0uXL37ttdcQFhaG119/HUDxHpq5c+fixIkTSE1NxY8//ojBgwejU6dOaNGihXHeCCLSW0XuzVRYWIg5c+agfv36sLW1RcuWLREbG6sVM2vWrFL3XfLz8zP0ZhDRC+BhYERUJQYMGIDbt29jxowZSE9PR0BAAGJjY6WT7q9evQqlsmK/n6SlpSEiIgIZGRnw9PTE4MGDMX36dOl5Gxsb7NmzB1FRUcjLy4OPjw/69u2LadOmVeq2EVHlK7k30+rVqxEYGIioqCiEhobiwoULcHd3LxU/bdo0fPfdd1izZg38/Pzw888/o3fv3khISECrVq2kuKZNm2LPnj3SdLVq/CpEJGcK8eSZqBYiOzsbzs7OyMrK4v0WiF6ApXyWLGU7iQytIp+lwMBAtGvXDsuXLwdQfOioj48Pxo0bh6lTp5aK9/LywmeffYYPP/xQmte3b1/Y2dnhu+++A1C8Z2X79u06rzJY2dtAVCny8gBHx+K/c3PN5kax+n6WeBgYERERycrz3JspPz9f675LAGBnZ1fqvkvJycnw8vJCvXr1MGjQIFy9erXctvC+S0TGxWKFiIiIZKW8ezOlp6eXuUxoaCiWLFmC5ORkaDQa7N69GzExMVr3XQoMDER0dDRiY2OxatUqpKSkoGPHjsjJydHZlvnz58PZ2Vl68DLmRFWLxQoRERGZvKVLl6Jhw4bw8/ODjY0Nxo4di2HDhmmdC9ejRw+88847aNGiBUJDQ7Fr1y5kZmZiy5YtOtcbGRmJrKws6XHt2rWq2Bwi+hOLFSIiIpKVit6bCQDc3Nywfft25OXl4cqVKzh//jwcHR1Rr149na/j4uKCRo0a4eLFizpjVCoVnJyctB5EVHVYrBAREZGsVPTeTE+ytbVF7dq18fjxY/zwww/l3ncpNzcXly5d4n2XiGSMxQoRERHJTkXuzQQAR44cQUxMDC5fvowDBw6ge/fu0Gg0mDJlihQzadIk7Nu3D6mpqUhISEDv3r1hZWWF8PDwKt8+ItIPLy5OREREslPRezM9evQI06ZNw+XLl+Ho6IiePXvi22+/hYuLixRz/fp1hIeH4+7du3Bzc8Mrr7yCw4cPw83Nrao3j4j0xPus8NhToudmKZ8lS9lOIkMzh8+SOWwDmRgLv88K96wQEREREcmVgwNgefsWJDxnhYiIiIiIZInFChERERERyRKLFSIiIiIikiWDFSv37t3DoEGD4OTkBBcXFwwfPhy5ubnlLvPo0SN8+OGHqFWrFhwdHdG3b99SN4RSKBSlHps2bTLUZhARERERkZEYrFgZNGgQzp49i927d+Onn37C/v378cEHH5S7zMSJE/Hvf/8bW7duxb59+3Dz5k306dOnVNy6deuQlpYmPcLCwgy0FUREREREZCwGuRrYuXPnEBsbi2PHjqFt27YAgGXLlqFnz574/PPP4eXlVWqZrKwsfP3119i4cSNeffVVAMVFSZMmTXD48GG8/PLLUqyLiwvUarUhmk5ERERERDJhkD0riYmJcHFxkQoVAAgJCYFSqcSRI0fKXObEiRMoLCxESEiINM/Pzw8vvfQSEhMTtWI//PBDuLq6on379li7di2edauY/Px8ZGdnaz2IiIiIiEjeDLJnJT09He7u7tovVK0aatasifT0dJ3L2NjYaN1pFgA8PDy0lpkzZw5effVV2Nvb45dffsGYMWOQm5uLjz76SGd75s+fj9mzZz//BhERERERUZWr0J6VqVOnlnmC+5OP8+fPG6qtAIDp06cjODgYrVq1wieffIIpU6Zg0aJF5S4TGRmJrKws6XHt2jWDtpGIiIiIiF5chfasfPzxxxg6dGi5MfXq1YNarcatW7e05j9+/Bj37t3Tea6JWq1GQUEBMjMztfauZGRklHt+SmBgIObOnYv8/HyoVKoyY1Qqlc7niIiIiIhInipUrLi5ucHNze2ZcUFBQcjMzMSJEyfQpk0bAMB///tfaDQaBAYGlrlMmzZtYG1tjbi4OPTt2xcAcOHCBVy9ehVBQUE6XyspKQk1atRgMUJEREREZGYMcs5KkyZN0L17d4wcORKrV69GYWEhxo4di4EDB0pXArtx4wa6deuGb775Bu3bt4ezszOGDx+OiIgI1KxZE05OThg3bhyCgoKkK4H9+9//RkZGBl5++WXY2tpi9+7dmDdvHiZNmmSIzSAiIiIiIiMySLECABs2bMDYsWPRrVs3KJVK9O3bF19++aX0fGFhIS5cuIAHDx5I87744gspNj8/H6GhoVi5cqX0vLW1NVasWIGJEydCCIEGDRpgyZIlGDlypKE2g4iIiIiIjEQhnnXdXzOUnZ0NZ2dnZGVlwcnJydjNITJZlvJZspTtJDI0c/gsmcM2EMmBvp8lg93BnoiIiIiI6EUY7DAwOSvZmcSbQxK9mJLPkLnvoGXOIKoc5pAzmA+IKoe++cAii5WcnBwAgI+Pj5FbQmQecnJy4OzsbOxmGAxzBlHlMuWcwXxAVLmelQ8s8pwVjUaDmzdvonr16lAoFMZuTqXIzs6Gj48Prl27xmNoZcoc+0gIgZycHHh5eUGpNN+jSs0tZ5jjWDRH5thP5pAzmA/IGMyxn/TNBxa5Z0WpVMLb29vYzTAIJycnsxnE5src+shUfx2tCHPNGeY2Fs2VufWTqecM5gMyJnPrJ33ygWn+rEFERERERGaPxQoREREREckSixUzoVKpMHPmTKhUKmM3hXRgH5FccCyaBvYTVQWOM9Ngyf1kkSfYExERERGR/HHPChERERERyRKLFSIiIiIikiUWK0REREREJEssVkxcdHQ0XFxcnhmnUCiwfft2g7fHUvB9J1PEcWscfN9JjjgujYPve8WxWDFxAwYMwB9//CFNz5o1CwEBAcZrkIWQ2/vOpEb6kNu4tRRye9+ZLwiQ37i0FHJ7300hH1jkHezNiZ2dHezs7IzdDIvD951MEcetcfB9JzniuDQOvu8Vxz0rMvTTTz/BxcUFRUVFAICkpCQoFApMnTpVihkxYgTee+89rd2J0dHRmD17Nk6fPg2FQgGFQoHo6GhpmTt37qB3796wt7dHw4YN8eOPP1blZsmesd73ffv2oX379lCpVPD09MTUqVPx+PFj6XlfX19ERUVpLRMQEIBZs2ZJzwNA7969oVAopGmyDMwXxsF8QXLEfGAczAeGxWJFhjp27IicnBycOnUKQPFgdHV1RXx8vBSzb98+dOnSRWu5AQMG4OOPP0bTpk2RlpaGtLQ0DBgwQHp+9uzZ6N+/P3799Vf07NkTgwYNwr1796pik0yCMd73GzduoGfPnmjXrh1Onz6NVatW4euvv8bf/vY3vdt97NgxAMC6deuQlpYmTZNlYL4wDuYLkiPmA+NgPjAsFisy5OzsjICAAGmQx8fHY+LEiTh16hRyc3Nx48YNXLx4EZ07d9Zazs7ODo6OjqhWrRrUajXUarXWrsahQ4ciPDwcDRo0wLx585Cbm4ujR49W5abJmjHe95UrV8LHxwfLly+Hn58fwsLCMHv2bCxevBgajUavdru5uQEAXFxcoFarpWmyDMwXxsF8QXLEfGAczAeGxWJFpjp37oz4+HgIIXDgwAH06dMHTZo0wcGDB7Fv3z54eXmhYcOGFVpnixYtpL8dHBzg5OSEW7duVXbTTVpVv+/nzp1DUFAQFAqFFBMcHIzc3Fxcv369cjaKzB7zhXEwX5AcMR8YB/OB4fAEe5nq0qUL1q5di9OnT8Pa2hp+fn7o0qUL4uPjcf/+/VLVuT6sra21phUKhd7Vt6WQ4/uuVCohhNCaV1hYWOF2kPmS47i1BHJ835kvSI7j0hLI8X03l3zAPSsyVXL84xdffCEN8JJBHx8fX+q4xxI2NjbSCV5UcVX9vjdp0gSJiYlayeTQoUOoXr06vL29ARTvpk1LS5Oez87ORkpKitZ6rK2t2e8WjPnCOJgvSI6YD4yD+cBwWKzIVI0aNdCiRQts2LBBGuCdOnXCyZMn8ccff+is0H19fZGSkoKkpCTcuXMH+fn5Vdhq01fV7/uYMWNw7do1jBs3DufPn8eOHTswc+ZMREREQKks/ni++uqr+Pbbb3HgwAGcOXMGQ4YMgZWVVanXj4uLQ3p6Ou7fv//8bwCZJOYL42C+IDliPjAO5gPDYbEiY507d0ZRUZE06GvWrAl/f3+o1Wo0bty4zGX69u2L7t27o2vXrnBzc8P3339fhS02D1X5vteuXRu7du3C0aNH0bJlS4waNQrDhw/HtGnTpJjIyEh07twZb775Jt544w2EhYWhfv36WutZvHgxdu/eDR8fH7Rq1er5NpxMGvOFcTBfkBwxHxgH84FhKMTTB7MRERERERHJAPesEBERERGRLLFYISIiIiIiWWKxQkREREREssRihYiIiIiIZInFChERERERyRKLFRM0a9YsBAQESNNDhw5FWFhYucvEx8dDoVAgMzPToG172tNttSTsJ5ILjkX5Yx9RVeFYkz/2kbZqBl07VYmlS5dq3cG0S5cuCAgIQFRUlDSvQ4cOSEtLg7OzsxFaSAD7ieSDY1H+2EdUVTjW5M/S+4jFihnQZ2Da2NhArVZXQWtIF/YTyQXHovyxj6iqcKzJn6X3EQ8De0GHDh1Cly5dYG9vjxo1aiA0NBT3798HAOTn5+Ojjz6Cu7s7bG1t8corr+DYsWPSsiW77OLi4tC2bVvY29ujQ4cOuHDhgtZrLFiwAB4eHqhevTqGDx+OR48eaT3/5O7BoUOHYt++fVi6dCkUCgUUCgVSU1PL3D34ww8/oGnTplCpVPD19cXixYu11uvr64t58+bhL3/5C6pXr46XXnoJ//d//6cV88knn6BRo0awt7dHvXr1MH36dBQWFr7o21rp2E+m0U+WgGNR/mORfST/PjIXHGvyH2vsIxn0kaDndurUKaFSqcTo0aNFUlKS+O2338SyZcvE7du3hRBCfPTRR8LLy0vs2rVLnD17VgwZMkTUqFFD3L17VwghxN69ewUAERgYKOLj48XZs2dFx44dRYcOHaTX2Lx5s1CpVOKf//ynOH/+vPjss89E9erVRcuWLaWYIUOGiF69egkhhMjMzBRBQUFi5MiRIi0tTaSlpYnHjx9Lr3X//n0hhBDHjx8XSqVSzJkzR1y4cEGsW7dO2NnZiXXr1knrrVOnjqhZs6ZYsWKFSE5OFvPnzxdKpVKcP39eipk7d644dOiQSElJET/++KPw8PAQ//jHP6TnZ86cqdVWY2A/mUY/WQKORfmPRfaR/PvIXHCsyX+ssY/k0UcsVl5AeHi4CA4OLvO53NxcYW1tLTZs2CDNKygoEF5eXmLhwoVCiP8N4j179kgxO3fuFADEw4cPhRBCBAUFiTFjxmitOzAwUOcgFkKIzp07i/Hjx2st8/Qgfvfdd8Vrr72mFTN58mTh7+8vTdepU0e899570rRGoxHu7u5i1apVOt4RIRYtWiTatGkjTRs70QjBfiqLHPvJEnAslia3scg+Kk1ufWQuONZKk9tYYx+VZow+4mFgLyApKQndunUr87lLly6hsLAQwcHB0jxra2u0b98e586d04pt0aKF9LenpycA4NatWwCAc+fOITAwUCs+KCjohdt+7tw5rbYBQHBwMJKTk1FUVFRm2xQKBdRqtdQ2ANi8eTOCg4OhVqvh6OiIadOm4erVqy/cvsrEfjKNfrIEHIvyH4vsI/n3kbngWJP/WGMfyaOPWKy8ADs7u0pZj7W1tfS3QqEAAGg0mkpZ94t6sm1AcftK2paYmIhBgwahZ8+e+Omnn3Dq1Cl89tlnKCgoMEZTdWI/mUY/WQKORfmPRfaR/PvIXHCsyX+ssY/k0UcsVl5AixYtEBcXV+Zz9evXh42NDQ4dOiTNKywsxLFjx+Dv76/3azRp0gRHjhzRmnf48OFyl7GxsdGqmnWt98m2AcUnkTVq1AhWVlZ6tS0hIQF16tTBZ599hrZt26Jhw4a4cuWKXstWJfaTafSTJeBYlP9YZB/Jv4/MBcea/Mca+0gefcRLF7+AyMhING/eHGPGjMGoUaNgY2ODvXv34p133oGrqytGjx6NyZMno2bNmnjppZewcOFCPHjwAMOHD9f7NcaPH4+hQ4eibdu2CA4OxoYNG3D27FnUq1dP5zK+vr44cuQIUlNT4ejoiJo1a5aK+fjjj9GuXTvMnTsXAwYMQGJiIpYvX46VK1fq3baGDRvi6tWr2LRpE9q1a4edO3fiX//6l97LVxX2k2n0kyXgWJT/WGQfyb+PzAXHmvzHGvtIHn3EPSsvoFGjRvjll19w+vRptG/fHkFBQdixYweqVSuuARcsWIC+ffvi/fffR+vWrXHx4kX8/PPPqFGjht6vMWDAAEyfPh1TpkxBmzZtcOXKFYwePbrcZSZNmgQrKyv4+/vDzc2tzGMLW7dujS1btmDTpk1o1qwZZsyYgTlz5mDo0KF6t+3tt9/GxIkTMXbsWAQEBCAhIQHTp0/Xe/mqwn4yjX6yBByL8h+L7CP595G54FiT/1hjH8mjjxRCPHFLTCIiIiIiIpngnhUiIiIiIpIlFitERERERCRLLFaIiIiIiEiWWKwQEREREZEssVghIiIiIiJZYrFCRERERESyxGKFiIiIiIhkicUKERERERHJEosVIiIiIiKSJRYrREREREQkSyxWiIiIiIhIlv4f923xMZTfltYAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1, 3, figsize=(8, 2), constrained_layout=True)\n", "for par, axi in zip(m.parameters, ax):\n", " axi.set_title(par)\n", " t = truth[par]\n", " axi.axhline(t, ls=\"--\", color=\"0.5\")\n", " axi.errorbar([\"with\\n conditional\"], m.values[par],\n", " m.errors[par], fmt=\"ok\")\n", " axi.errorbar([\"without\\n conditional\"], m2.values[par],\n", " m2.errors[par], fmt=\"or\")\n", " axi.set_xlim(-0.5, 1.5)\n", " dt = 2 * m2.errors[par]\n", " axi.set_ylim(t - dt, t + dt)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.8.14 ('venv': venv)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.8" }, "vscode": { "interpreter": { "hash": "bdbf20ff2e92a3ae3002db8b02bd1dd1b287e934c884beb29a73dced9dbd0fa3" } } }, "nbformat": 4, "nbformat_minor": 5 }