{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Template fits: comparing two chi-square distributed test statistics\n", "\n", "The builtin binned cost functions in `iminuit.cost` provide a method `pulls` which returns the pull distribution for a given data set. The sum of pulls squared is asymptotically chi-square distributed, but as is explained in the documentation, it is better to use the minimum value of the cost function instead of the pulls to compute this test statistic. The reason is that the cost function has been designed so that the minimum is chi-square distributed and it approaches the asymptotic limit faster than a simple sum of pulls squared.\n", "\n", "We demonstrate this in this example for a Template fit. We generate random samples from a mixed model (gaussian peak over exponential background), which is fitted many times. The distribution of the cost function minimum and alternatively the sum of pulls squared is compared with a chi-square distribution by computing p-values based on the expected distribution. If the test statistic follows the chi-square distribution, the distribution of p-values must be uniform." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "from iminuit import cost, Minuit\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from scipy.stats import chi2\n", "from IPython.display import display" ] }, { "cell_type": "code", "execution_count": 17, "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 = 13.87 (χ²/ndof = 0.8) Nfcn = 101
EDM = 2.91e-08 (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", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0 1.01e3 0.05e3 0
1 x1 870 50 0
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x0 x1
x0 2.49e+03 -0.5e3 (-0.196)
x1 -0.5e3 (-0.196) 2.21e+03
\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-04-27T12:51:54.476322\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.1, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 13.87 (χ²/ndof = 0.8) │ Nfcn = 101 │\n", "│ EDM = 2.91e-08 (Goal: 0.0002) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ No Parameters at limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Below EDM threshold (goal x 10) │ Below call limit │\n", "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n", "│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n", "┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", "├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", "│ 0 │ x0 │ 1.01e3 │ 0.05e3 │ │ │ 0 │ │ │\n", "│ 1 │ x1 │ 870 │ 50 │ │ │ 0 │ │ │\n", "└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌────┬───────────────────┐\n", "│ │ x0 x1 │\n", "├────┼───────────────────┤\n", "│ x0 │ 2.49e+03 -0.5e3 │\n", "│ x1 │ -0.5e3 2.21e+03 │\n", "└────┴───────────────────┘" ] }, "metadata": {}, "output_type": "display_data" }, { "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 = 201.8 (χ²/ndof = 1.0) Nfcn = 100
EDM = 2.37e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0 1.02e3 0.05e3 0
1 x1 810 40 0
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x0 x1
x0 2.41e+03 -0.3e3 (-0.161)
x1 -0.3e3 (-0.161) 1.92e+03
\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-04-27T12:51:54.631153\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.1, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 201.8 (χ²/ndof = 1.0) │ Nfcn = 100 │\n", "│ EDM = 2.37e-05 (Goal: 0.0002) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ No Parameters at limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Below EDM threshold (goal x 10) │ Below call limit │\n", "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n", "│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n", "┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", "├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", "│ 0 │ x0 │ 1.02e3 │ 0.05e3 │ │ │ 0 │ │ │\n", "│ 1 │ x1 │ 810 │ 40 │ │ │ 0 │ │ │\n", "└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌────┬───────────────────┐\n", "│ │ x0 x1 │\n", "├────┼───────────────────┤\n", "│ x0 │ 2.41e+03 -0.3e3 │\n", "│ x1 │ -0.3e3 1.92e+03 │\n", "└────┴───────────────────┘" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xr = (0, 2) # xrange\n", "rng = np.random.default_rng(1)\n", "\n", "nmc = 1000\n", "trials = 1000\n", "\n", "data = {}\n", "data2 = {}\n", "\n", "first = True\n", "for trial in range(trials):\n", " for bins in (20, 200,):\n", " xdata = rng.normal(1, 0.1, size=1000)\n", " ydata = rng.exponential(size=len(xdata))\n", " xmix = np.append(xdata, ydata)\n", " xmix = xmix[(xr[0] < xmix) & (xmix < xr[1])]\n", "\n", " n, xe = np.histogram(xmix, bins=bins, range=xr)\n", "\n", " x = rng.normal(1, 0.1, size=nmc)\n", " y = rng.exponential(size=nmc)\n", " t = [\n", " np.histogram(x, bins=bins, range=xr)[0],\n", " np.histogram(y, bins=bins, range=xr)[0],\n", " ]\n", " c = cost.Template(n, xe, t)\n", " m = Minuit(c, 1, 1)\n", " m.migrad()\n", " assert m.valid\n", " assert m.accurate\n", " data.setdefault(bins, []).append(m.fmin.fval)\n", " data2.setdefault(bins, []).append(np.nansum(c.pulls(m.values) ** 2))\n", " # display one example fit\n", " if first:\n", " display(m)\n", " first = False\n", "\n", "for key in tuple(data):\n", " val = data[key]\n", " data[key] = np.array(val)\n", " val = data2[key]\n", " data2[key] = np.array(val)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGzCAYAAAAFROyYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA04ElEQVR4nO3deXyNZ/7/8fdJIosliViyzESTqrHVVlSDDiXfxjIqbWZoq0WVGE1M0dbS1laKmhZjH8ZYHsPojC+qpVqToqOCCjFaoZZYOpJohySWkYTcvz/6c749pJLDSXIdeT0fj/N4OPd9net87os471z3dd/HZlmWJQAAAIN4lHcBAAAANyOgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAFdSECRNks9n0/fffF9s2IiJC/fv3L/2iAOD/I6AAuOetXbtWvXv31v3336/KlSurfv36euWVV5SdnV1k+w0bNuihhx6Sr6+v6tSpo/Hjx+vatWtlWzRQwXmVdwEAzHfkyBF5eLjv7zPx8fEKCwvTc889pzp16ujgwYOaO3euNm3apH379snPz8/e9uOPP1ZsbKw6duyoOXPm6ODBg5o8ebLOnTunBQsWlONRABULAQVAsXx8fMq7hLuyZs0adezY0WFby5Yt1a9fP61cuVIDBw60b3/11VfVtGlTffrpp/Ly+uG/SH9/f02ZMkUvv/yyGjRoUJalAxWW+/5KBMAlvv/+e/Xq1Uv+/v6qUaOGXn75ZV29etWhzc1rUJYtWyabzaYvvvhCI0aMUK1atVSlShU9+eST+u677xxeu3fvXsXExKhmzZry8/NTZGSkBgwYUBaHZndzOJGkJ598UpKUlpZm33bo0CEdOnRI8fHx9nAiSS+99JIsy9KaNWtKvVYAP2AGBajgevXqpYiICE2dOlW7du3S7NmzdeHCBa1YsaLY1w4dOlTVq1fX+PHjdfLkSc2aNUuJiYl6//33JUnnzp3T448/rlq1amn06NEKDAzUyZMntXbt2mL7vnTp0i1BqSiVKlVSQEBA8Qd6k8zMTElSzZo17dv2798vSWrVqpVD27CwMP385z+37wdQ+ggoQAUXGRmpDz74QJKUkJAgf39/zZ8/336q43Zq1KihTz/9VDabTZJUWFio2bNnKycnRwEBAdq5c6cuXLigTz/91OFDf/LkycXWlZiYqOXLlxfbrkOHDtq2bVux7W72zjvvyNPTU7/+9a/t2zIyMiRJoaGht7QPDQ3V2bNnnX4fAHeGgAJUcAkJCQ7Phw4dqvnz52vTpk3FBpT4+Hh7OJGkRx99VDNnztSpU6fUtGlTBQYGSpI++ugjNWvWTJUqVSpxXSNHjtRzzz1XbLvq1auXuM8bVq1apSVLlmjkyJGqV6+efft///tfSUWvufH19VVubq7T7wXgzhBQgAruxx/QklS3bl15eHjo5MmTxb62Tp06Ds9vhIULFy5I+mF2Iy4uThMnTtTMmTPVsWNHxcbG6tlnny124W2jRo3UqFEjJ46kZP75z3/qxRdfVExMjN5++22HfTeu5snLy7vldVevXnW42gdA6SKgAHDw4xmR4nh6eha53bIse19r1qzRrl279OGHH+qTTz7RgAED9N5772nXrl2qWrXqT/adk5Njn9G4HW9vbwUFBZWo3gMHDuiJJ57Qgw8+qDVr1jgshJX+79RORkaGwsPDHfZlZGTo4YcfLtH7ALh7XMUDVHBHjx51eH7s2DEVFhYqIiLCZe/xyCOP6O2339bevXu1cuVKff3111q9evVtX/Pyyy8rNDS02MdTTz1VohqOHz+uLl26qHbt2tq0aVOR4ah58+aSfrjy6MfOnj2rb7/91r4fQOljBgWo4ObNm6fHH3/c/nzOnDmSpK5du9513xcuXFBgYKDDrMyND/miTqP8mCvXoGRmZurxxx+Xh4eHPvnkE9WqVavIdo0bN1aDBg20aNEiDR482D5DtGDBAtlsNocFtQBKFwEFqODS09P1xBNPqEuXLkpOTtZf/vIXPfvss2rWrNld9718+XLNnz9fTz75pOrWrauLFy9q8eLF8vf3V7du3W77WleuQenSpYtOnDihkSNHaseOHdqxY4d9X3BwsP7nf/7H/vz3v/+9nnjiCT3++ON6+umn9dVXX2nu3LkaOHCgGjZs6JJ6ABSPgAJUcO+//77GjRun0aNHy8vLS4mJifr973/vkr47dOigPXv2aPXq1crKylJAQIAefvhhrVy5UpGRkS55j5I4cOCAJGn69OlF1vjjgPKrX/1Ka9eu1cSJEzV06FDVqlVLr7/+usaNG1dm9QKQbNaN1WwAAACGYJEsAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBx3PI+KIWFhTp79qyqVavm1PeGAACA8mNZli5evKiwsDB5eNx+jsQtA8rZs2dv+SIvAADgHs6cOaOf//znt23jlgGlWrVqkn44QH9//3KuBgAAlERubq7Cw8Ptn+O345YB5cZpHX9/fwIKAABupiTLM1gkCwAAjENAAQAAxiGgAAAA47jlGhQAMJVlWbp27ZquX79e3qUAZc7T01NeXl4uuQUIAQUAXCQ/P18ZGRm6cuVKeZcClJvKlSsrNDRU3t7ed9UPAQUAXKCwsFDp6eny9PRUWFiYvL29uZEkKhTLspSfn6/vvvtO6enpqlevXrE3Y7sdAgoAuEB+fr4KCwsVHh6uypUrl3c5QLnw8/NTpUqVdOrUKeXn58vX1/eO+2KRLAC40N38xgjcC1z1M8BPEgAAMA4BBQAAGIc1KABQyiJGbyzT9zs5rXuZvp+rTJgwQQsWLNC5c+e0bt06xcbGllstNput3Guo6JhBAQCUmgkTJqh58+bFtktLS9PEiRP1xz/+URkZGeratWvpF6efrq8sa0DRmEEBAJS748ePS5J69uxpxOXZISEh5V1ChccMCgBUcIWFhZo+fboeeOAB+fj4qE6dOnr77bft+w8ePKhOnTrJz89PNWrUUHx8vC5dumTfv23bNj388MOqUqWKAgMD1a5dO506dUrLli3TxIkTdeDAAdlsNtlsNi1btuyW958wYYJ69Ogh6YcrQG4ElI4dO2rYsGEObWNjY9W/f3/784iICE2ZMkUDBgxQtWrVVKdOHS1atMjhNd9++62eeeYZBQUFqUqVKmrVqpV279592/psNpvWr19f4jHo37+/YmNj9e677yo0NFQ1atRQQkKCCgoKnPmrwI8wgwI4Y+vU8q7AeY+NKe8KYLgxY8Zo8eLFmjlzptq3b6+MjAwdPnxYknT58mXFxMQoKipKX375pc6dO6eBAwcqMTFRy5Yt07Vr1xQbG6tBgwbpr3/9q/Lz87Vnzx7ZbDb17t1bX331lTZv3qx//OMfkqSAgIBb3v/VV19VRESEXnjhBWVkZDhd/3vvvadJkybp9ddf15o1azRkyBB16NBB9evX16VLl9ShQwf97Gc/04YNGxQSEqJ9+/apsLCwxPUVNwY3bN26VaGhodq6dauOHTum3r17q3nz5ho0aJDTxwQCCgBUaBcvXtQf/vAHzZ07V/369ZMk1a1bV+3bt5ckrVq1SlevXtWKFStUpUoVSdLcuXPVo0cPvfPOO6pUqZJycnL0q1/9SnXr1pUkNWzY0N5/1apV5eXlddtTJlWrVlVgYKCkOzu10q1bN7300kuSpFGjRmnmzJnaunWr6tevr1WrVum7777Tl19+qaCgIEnSAw884FR9xY1BcHCwJKl69eqaO3euPD091aBBA3Xv3l1JSUkElDvk9Cmezz//XD169FBYWNgtU2AFBQUaNWqUmjRpoipVqigsLEx9+/bV2bNnHfo4f/68+vTpI39/fwUGBurFF190mCoDAJSNtLQ05eXlqXPnzj+5v1mzZvYPZklq166dCgsLdeTIEQUFBal///6KiYlRjx499Ic//OGOZkHuRtOmTe1/ttlsCgkJ0blz5yRJqampatGihT2c3InixuCGxo0by9PT0/48NDTUXgec53RAuXz5spo1a6Z58+bdsu/KlSvat2+fxo4dq3379mnt2rU6cuSInnjiCYd2ffr00ddff60tW7boo48+0ueff674+Pg7PwoAwB3x8/O76z6WLl2q5ORktW3bVu+//75+8YtfaNeuXXfdr4eHhyzLcthW1JqOSpUqOTy32WwqLCyU5JrjK6nb1QHnOR1QunbtqsmTJ+vJJ5+8ZV9AQIC2bNmiXr16qX79+nrkkUc0d+5cpaSk6PTp05J+SKKbN2/Wn/70J7Vp00bt27fXnDlztHr16ltmWgAApatevXry8/NTUlJSkfsbNmyoAwcO6PLly/ZtX3zxhTw8PFS/fn37thYtWmjMmDHauXOnHnzwQa1atUqS5O3trevXr99RbbVq1XKYjbl+/bq++uorp/po2rSpUlNTdf78+SL3l6S+ko4BXKvUr+LJycmRzWazn19MTk5WYGCgWrVqZW8THR0tDw8P7d69u8g+8vLylJub6/AAANw9X19fjRo1SiNHjtSKFSt0/Phx7dq1S0uWLJH0w4y3r6+v+vXrp6+++kpbt27V0KFD9fzzzys4OFjp6ekaM2aMkpOTderUKX366ac6evSofR1KRESE0tPTlZqaqu+//155eXklrq1Tp07auHGjNm7cqMOHD2vIkCHKzs526vieeeYZhYSEKDY2Vl988YVOnDih//3f/1VycnKJ6ytuDFA6SnWR7NWrVzVq1Cg988wz8vf3lyRlZmaqdu3ajkV4eSkoKEiZmZlF9jN16lRNnDixNEsFgFJj+p1dx44dKy8vL40bN05nz55VaGiofvvb30qSKleurE8++UQvv/yyWrdurcqVKysuLk4zZsyw7z98+LCWL1+u//znPwoNDVVCQoIGDx4sSYqLi9PatWv12GOPKTs7W0uXLnW4TPh2BgwYoAMHDqhv377y8vLS8OHD9dhjjzl1bN7e3vr000/1yiuvqFu3brp27ZoaNWpkX6ZQkvqKGwOUDpt18wk+Z158m1sBFxQUKC4uTt9++622bdtmDyhTpkzR8uXLHRYWSVLt2rU1ceJEDRky5Ja+8vLyHFJtbm6uwsPDlZOTY+8XKBNcZoyfcPXqVaWnpysyMvKuvmIecHe3+1nIzc1VQEBAiT6/S2UGpaCgQL169dKpU6f02WefORTx49XVN1y7dk3nz5//ycu8fHx85OPjUxqlAgAAA7k8oNwIJ0ePHtXWrVtVo0YNh/1RUVHKzs5WSkqKWrZsKUn67LPPVFhYqDZt2ri6HMClZiV9U94lOG2YczPiAGAEpwPKpUuXdOzYMfvzG4uLgoKCFBoaql//+tfat2+fPvroI12/ft2+riQoKEje3t5q2LChunTpokGDBmnhwoUqKChQYmKinn76aYWFhbnuyAAAgNtyOqDs3bvXYZHSiBEjJEn9+vXThAkTtGHDBkm65dsht27dqo4dO0qSVq5cqcTERHXu3FkeHh6Ki4vT7Nmz7/AQAADAvcbpgNKxY8dbbpzzYyVZcxsUFGS/Rh4AAOBmfJsxAAAwDgEFAAAYh4ACAACMQ0ABAADGKdVb3QMAVPZ3IL6H7x68aNEiTZo0Sf/+9781Y8YMDRs2rFTeZ8KECVq/fr1SU1MlSf3791d2drbWr19fKu/nbiIiIjRs2LBSG3+JgAIAcBO5ublKTEzUjBkzFBcXp4CAgPIuCaWIgAIAcAunT59WQUGBunfvrtDQ0PIux+3l5+fL29u7vMv4SaxBAYAKbs2aNWrSpIn8/PxUo0YNRUdH6/Lly5J+uPfVzdP4sbGxDt/4GxERocmTJ6tv376qWrWq7rvvPm3YsEHfffedevbsqapVq6pp06bau3fvbes4ffq0vb2/v7969eqlrKwsSdKyZcvUpEkTSdL9998vm82mkydP3tLHyZMnZbPZtHr1arVt21a+vr568MEHtX37dnubZcuWKTAw0OF169evl81mK+GI3X7MbnbhwgX16dNHtWrVkp+fn+rVq6elS5fa9+/Zs0ctWrSQr6+vWrVqpXXr1slms9lPL5Wk3uPHj6tnz54KDg5W1apV1bp1a/3jH/9weE1ERIQmTZqkvn37yt/fX/Hx8ZKkHTt26NFHH5Wfn5/Cw8P1u9/9zuFYzp07px49esjPz0+RkZFauXJlicfpbhBQAKACy8jI0DPPPKMBAwYoLS1N27Zt01NPPVWim27+2MyZM9WuXTvt379f3bt31/PPP6++ffvqueee0759+1S3bl317dv3J/stLCxUz549df78eW3fvl1btmzRiRMn1Lt3b0lS79697R+4e/bsUUZGhsLDw3+yntdee02vvPKK9u/fr6ioKPXo0UP/+c9/nDqmn+LsmI0dO1aHDh3Sxx9/rLS0NC1YsEA1a9aU9MPXx/zqV79So0aNlJKSogkTJujVV191uqZLly6pW7duSkpK0v79+9WlSxf16NFDp0+fdmj37rvvqlmzZtq/f7/Gjh2r48ePq0uXLoqLi9O//vUvvf/++9qxY4cSExPtr+nfv7/OnDmjrVu3as2aNZo/f/4tX/pbGjjFAwAVWEZGhq5du6annnpK9913nyTZZyqc0a1bNw0ePFiSNG7cOC1YsECtW7fWb37zG0nSqFGjFBUVpaysrCK/uT4pKUkHDx5Uenq6PXisWLFCjRs31pdffqnWrVvbv3y2Vq1aRfbxY4mJiYqLi5MkLViwQJs3b9aSJUs0cuRIp4/tZs6O2enTp9WiRQu1atVK0g8zGTesWrVKhYWFWrJkiXx9fdW4cWN9++23GjJkiFM1NWvWTM2aNbM/nzRpktatW6cNGzY4hI1OnTrplVdesT8fOHCg+vTpY58lq1evnmbPnq0OHTpowYIFOn36tD7++GPt2bNHrVu3liQtWbJEDRs2dKq+O8EMCgBUYM2aNVPnzp3VpEkT/eY3v9HixYt14cIFp/tp2rSp/c/BwcGSHD+0b2z7qd+809LSFB4e7jAr0qhRIwUGBiotLc3peqKioux/9vLyUqtWre6on6I4O2ZDhgzR6tWr1bx5c40cOVI7d+6070tLS1PTpk3l6+tbZO0ldenSJb366qtq2LChAgMDVbVqVaWlpd0yg3IjJN1w4MABLVu2TFWrVrU/YmJiVFhYqPT0dKWlpcnLy0stW7a0v6ZBgwa3nHIqDQQUAKjAPD09tWXLFn388cdq1KiR5syZo/r16ys9PV2S5OHhccupi4KCglv6qVSpkv3PN9ZGFLWtsLDQ5cfgrJIe008pbsxu1rVrV506dUrDhw/X2bNn1blzZ6dO45Sk3ldffVXr1q3TlClT9M9//lOpqalq0qSJ8vPzHdpVqVLF4fmlS5c0ePBgpaam2h8HDhzQ0aNHVbdu3RLXWBoIKABQwdlsNrVr104TJ07U/v375e3trXXr1kn64XRKRkaGve3169f11VdfubyGhg0b6syZMzpz5ox926FDh5Sdna1GjRo53d+uXbvsf7527ZpSUlLspyVq1aqlixcvOiwEvbEgtaRuN2ZFqVWrlvr166e//OUvmjVrlhYtWiTph+P+17/+patXrxZZe0nr/eKLL9S/f389+eSTatKkiUJCQopcRHyzhx56SIcOHdIDDzxwy8Pb21sNGjSwj98NR44cUXZ2drF93y0CCgBUYLt379aUKVO0d+9enT59WmvXrtV3331n/zDv1KmTNm7cqI0bN+rw4cMaMmRIqXw4RUdHq0mTJurTp4/27dunPXv2qG/fvurQocMtpyVKYt68eVq3bp0OHz6shIQEXbhwQQMGDJAktWnTRpUrV9brr7+u48ePa9WqVVq2bFmJ+y5uzG42btw4ffDBBzp27Ji+/vprffTRR/a2zz77rGw2mwYNGqRDhw5p06ZNevfddx1eX5J669Wrp7Vr19pnQJ599tkSzVaNGjVKO3fuVGJiolJTU3X06FF98MEH9nUr9evXV5cuXTR48GDt3r1bKSkpGjhwoPz8/Eo8XneKRbIAUNoMvrOrv7+/Pv/8c82aNUu5ubm677779N5776lr166SpAEDBujAgQPq27evvLy8NHz4cD322GMur8Nms+mDDz7Q0KFD9ctf/lIeHh7q0qWL5syZc0f9TZs2TdOmTVNqaqoeeOABbdiwwX7lTFBQkP7yl7/otdde0+LFi9W5c2dNmDDBftltcYobs5t5e3trzJgxOnnypPz8/PToo49q9erVkqSqVavqww8/1G9/+1u1aNFCjRo10jvvvGNf4FvSemfMmKEBAwaobdu2qlmzpkaNGqXc3Nxij6Vp06bavn273njjDT366KOyLEt169a1Xz0lSUuXLtXAgQPVoUMHBQcHa/LkyRo7dmyJxupu2CxnryUzQG5urgICApSTkyN/f//yLgcVyKw3XyjvEpw2bPLS4hvhrl29elXp6emKjIx0WPCIsnXy5ElFRkZq//79at68eXmXc0fc/Rhu97PgzOc3p3gAAIBxCCgAAMA4rEEBANwzIiIinL4LrmnuhWNwBQIKAFRUuRnFtzGNP18SWFFwigcAXIjffFHRuepngIACAC5w466pV65cKedKgPJ142fgx3cSvhOc4gEAF/D09FRgYKD9u2YqV65sv727sfJLfnt3Y/zojqswi2VZunLlis6dO6fAwEB5enreVX8EFABwkRvfsFsWX0XvEldzyrsC5/leLr4NylVgYGCx3zZdEgQUAHARm82m0NBQ1a5d26kvnys3u/9Y3hU4r+Hg8q4At1GpUqW7njm5gYACAC7m6enpsv+kS9OsLanlXYLThnXgLr0VBYtkAQCAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHbzMGAFfYOrW8KwDuKcygAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACM43RA+fzzz9WjRw+FhYXJZrNp/fr1Dvsty9K4ceMUGhoqPz8/RUdH6+jRow5tzp8/rz59+sjf31+BgYF68cUXdenSpbs6EAAAcO9w+kZtly9fVrNmzTRgwAA99dRTt+yfPn26Zs+ereXLlysyMlJjx45VTEyMDh06JF9fX0lSnz59lJGRoS1btqigoEAvvPCC4uPjtWrVqrs/IgAoB7OSvinvEoB7itMBpWvXruratWuR+yzL0qxZs/Tmm2+qZ8+ekqQVK1YoODhY69ev19NPP620tDRt3rxZX375pVq1aiVJmjNnjrp166Z3331XYWFhd3E4AADgXuDSNSjp6enKzMxUdHS0fVtAQIDatGmj5ORkSVJycrICAwPt4USSoqOj5eHhod27dxfZb15ennJzcx0eAADg3uXS7+LJzMyUJAUHBztsDw4Otu/LzMxU7dq1HYvw8lJQUJC9zc2mTp2qiRMnurLUew/fAwIAuIe4xVU8Y8aMUU5Ojv1x5syZ8i4JAACUIpcGlJCQEElSVlaWw/asrCz7vpCQEJ07d85h/7Vr13T+/Hl7m5v5+PjI39/f4QEAAO5dLj3FExkZqZCQECUlJal58+aSpNzcXO3evVtDhgyRJEVFRSk7O1spKSlq2bKlJOmzzz5TYWGh2rRp48pyKhR3vIJgWOdflHcJAFDqIkZvLO8S7sjJad3L9f2dDiiXLl3SsWPH7M/T09OVmpqqoKAg1alTR8OGDdPkyZNVr149+2XGYWFhio2NlSQ1bNhQXbp00aBBg7Rw4UIVFBQoMTFRTz/9NFfwAAAASXcQUPbu3avHHnvM/nzEiBGSpH79+mnZsmUaOXKkLl++rPj4eGVnZ6t9+/bavHmz/R4okrRy5UolJiaqc+fO8vDwUFxcnGbPnu2CwwEAAPcCpwNKx44dZVnWT+632Wx666239NZbb/1km6CgIG7KBpQRd5xeLu+pZQDlzy2u4gEAABWLSxfJAs5wx4W97miY15ryLuEOMIMCVHTMoAAAAOMQUAAAgHE4xQPAOO64sHcY/5sCLsUMCgAAMA4BBQAAGIdJySLMevOF8i4BAIAKjRkUAABgHAIKAAAwDqd4ABjHPW8uB8CVmEEBAADGIaAAAADjcIoHAIBS5L6nLMv3O7GYQQEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACM41XeBQAAUFIRozeWdwlOG8Yn7R1hBgUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA63jwEAuI1hXmvKuwSUEWZQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGcXlAuX79usaOHavIyEj5+fmpbt26mjRpkizLsrexLEvjxo1TaGio/Pz8FB0draNHj7q6FAAA4KZcHlDeeecdLViwQHPnzlVaWpreeecdTZ8+XXPmzLG3mT59umbPnq2FCxdq9+7dqlKlimJiYnT16lVXlwMAANyQy+8ku3PnTvXs2VPdu3eXJEVEROivf/2r9uzZI+mH2ZNZs2bpzTffVM+ePSVJK1asUHBwsNavX6+nn376lj7z8vKUl5dnf56bm+vqsgEAgEFcPoPStm1bJSUl6ZtvvpEkHThwQDt27FDXrl0lSenp6crMzFR0dLT9NQEBAWrTpo2Sk5OL7HPq1KkKCAiwP8LDw11dNgAAMIjLZ1BGjx6t3NxcNWjQQJ6enrp+/brefvtt9enTR5KUmZkpSQoODnZ4XXBwsH3fzcaMGaMRI0bYn+fm5hJSAAC4h7k8oPztb3/TypUrtWrVKjVu3FipqakaNmyYwsLC1K9fvzvq08fHRz4+Pi6uFAAAmMrlAeW1117T6NGj7WtJmjRpolOnTmnq1Knq16+fQkJCJElZWVkKDQ21vy4rK0vNmzd3dTkAAMANuXwNypUrV+Th4ditp6enCgsLJUmRkZEKCQlRUlKSfX9ubq52796tqKgoV5cDAADckMtnUHr06KG3335bderUUePGjbV//37NmDFDAwYMkCTZbDYNGzZMkydPVr169RQZGamxY8cqLCxMsbGxri4HAAC4IZcHlDlz5mjs2LF66aWXdO7cOYWFhWnw4MEaN26cvc3IkSN1+fJlxcfHKzs7W+3bt9fmzZvl6+vr6nIAAIAbslk/vsWrm8jNzVVAQIBycnLk7+/v8v5nvfmCy/sEAMCdDJu81OV9OvP5zXfxAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGKdUAsq///1vPffcc6pRo4b8/PzUpEkT7d27177fsiyNGzdOoaGh8vPzU3R0tI4ePVoapQAAADfk8oBy4cIFtWvXTpUqVdLHH3+sQ4cO6b333lP16tXtbaZPn67Zs2dr4cKF2r17t6pUqaKYmBhdvXrV1eUAAAA35OXqDt955x2Fh4dr6dKl9m2RkZH2P1uWpVmzZunNN99Uz549JUkrVqxQcHCw1q9fr6efftrVJQEAADfj8hmUDRs2qFWrVvrNb36j2rVrq0WLFlq8eLF9f3p6ujIzMxUdHW3fFhAQoDZt2ig5ObnIPvPy8pSbm+vwAAAA9y6XB5QTJ05owYIFqlevnj755BMNGTJEv/vd77R8+XJJUmZmpiQpODjY4XXBwcH2fTebOnWqAgIC7I/w8HBXlw0AAAzi8oBSWFiohx56SFOmTFGLFi0UHx+vQYMGaeHChXfc55gxY5STk2N/nDlzxoUVAwAA07g8oISGhqpRo0YO2xo2bKjTp09LkkJCQiRJWVlZDm2ysrLs+27m4+Mjf39/hwcAALh3uTygtGvXTkeOHHHY9s033+i+++6T9MOC2ZCQECUlJdn35+bmavfu3YqKinJ1OQAAwA25/Cqe4cOHq23btpoyZYp69eqlPXv2aNGiRVq0aJEkyWazadiwYZo8ebLq1aunyMhIjR07VmFhYYqNjXV1OQAAwA25PKC0bt1a69at05gxY/TWW28pMjJSs2bNUp8+fextRo4cqcuXLys+Pl7Z2dlq3769Nm/eLF9fX1eXAwAA3JDNsiyrvItwVm5urgICApSTk1Mq61FmvfmCy/sEAMCdDJu8tPhGTnLm85vv4gEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABin1APKtGnTZLPZNGzYMPu2q1evKiEhQTVq1FDVqlUVFxenrKys0i4FAAC4iVINKF9++aX++Mc/qmnTpg7bhw8frg8//FB///vftX37dp09e1ZPPfVUaZYCAADcSKkFlEuXLqlPnz5avHixqlevbt+ek5OjJUuWaMaMGerUqZNatmyppUuXaufOndq1a1dplQMAANxIqQWUhIQEde/eXdHR0Q7bU1JSVFBQ4LC9QYMGqlOnjpKTk4vsKy8vT7m5uQ4PAABw7/IqjU5Xr16tffv26csvv7xlX2Zmpry9vRUYGOiwPTg4WJmZmUX2N3XqVE2cOLE0SgUAAAZy+QzKmTNn9PLLL2vlypXy9fV1SZ9jxoxRTk6O/XHmzBmX9AsAAMzk8oCSkpKic+fO6aGHHpKXl5e8vLy0fft2zZ49W15eXgoODlZ+fr6ys7MdXpeVlaWQkJAi+/Tx8ZG/v7/DAwAA3Ltcfoqnc+fOOnjwoMO2F154QQ0aNNCoUaMUHh6uSpUqKSkpSXFxcZKkI0eO6PTp04qKinJ1OQAAwA25PKBUq1ZNDz74oMO2KlWqqEaNGvbtL774okaMGKGgoCD5+/tr6NChioqK0iOPPOLqcgAAgBsqlUWyxZk5c6Y8PDwUFxenvLw8xcTEaP78+eVRCgAAMFCZBJRt27Y5PPf19dW8efM0b968snh7AADgZvguHgAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHFcHlCmTp2q1q1bq1q1aqpdu7ZiY2N15MgRhzZXr15VQkKCatSooapVqyouLk5ZWVmuLgUAALgplweU7du3KyEhQbt27dKWLVtUUFCgxx9/XJcvX7a3GT58uD788EP9/e9/1/bt23X27Fk99dRTri4FAAC4KS9Xd7h582aH58uWLVPt2rWVkpKiX/7yl8rJydGSJUu0atUqderUSZK0dOlSNWzYULt27dIjjzzi6pIAAICbKfU1KDk5OZKkoKAgSVJKSooKCgoUHR1tb9OgQQPVqVNHycnJRfaRl5en3NxchwcAALh3lWpAKSws1LBhw9SuXTs9+OCDkqTMzEx5e3srMDDQoW1wcLAyMzOL7Gfq1KkKCAiwP8LDw0uzbAAAUM5KNaAkJCToq6++0urVq++qnzFjxignJ8f+OHPmjIsqBAAAJnL5GpQbEhMT9dFHH+nzzz/Xz3/+c/v2kJAQ5efnKzs722EWJSsrSyEhIUX25ePjIx8fn9IqFQAAGMblMyiWZSkxMVHr1q3TZ599psjISIf9LVu2VKVKlZSUlGTfduTIEZ0+fVpRUVGuLgcAALghl8+gJCQkaNWqVfrggw9UrVo1+7qSgIAA+fn5KSAgQC+++KJGjBihoKAg+fv7a+jQoYqKiuIKHgAAIKkUAsqCBQskSR07dnTYvnTpUvXv31+SNHPmTHl4eCguLk55eXmKiYnR/PnzXV0KAABwUy4PKJZlFdvG19dX8+bN07x581z99gAA4B7Ad/EAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYp1wDyrx58xQRESFfX1+1adNGe/bsKc9yAACAIcotoLz//vsaMWKExo8fr3379qlZs2aKiYnRuXPnyqskAABgiHILKDNmzNCgQYP0wgsvqFGjRlq4cKEqV66sP//5z+VVEgAAMIRXebxpfn6+UlJSNGbMGPs2Dw8PRUdHKzk5+Zb2eXl5ysvLsz/PycmRJOXm5pZKfVfz8kulXwAA3EVpfMbe6NOyrGLblktA+f7773X9+nUFBwc7bA8ODtbhw4dvaT916lRNnDjxlu3h4eGlViMAABXZmHdXlVrfFy9eVEBAwG3blEtAcdaYMWM0YsQI+/PCwkKdP39eNWrUkM1mc+l75ebmKjw8XGfOnJG/v79L+8b/YZzLBuNcNhjnssE4l53SGmvLsnTx4kWFhYUV27ZcAkrNmjXl6emprKwsh+1ZWVkKCQm5pb2Pj498fHwctgUGBpZmifL39+cHoAwwzmWDcS4bjHPZYJzLTmmMdXEzJzeUyyJZb29vtWzZUklJSfZthYWFSkpKUlRUVHmUBAAADFJup3hGjBihfv36qVWrVnr44Yc1a9YsXb58WS+88EJ5lQQAAAxRbgGld+/e+u677zRu3DhlZmaqefPm2rx58y0LZ8uaj4+Pxo8ff8spJbgW41w2GOeywTiXDca57Jgw1jarJNf6AAAAlCG+iwcAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEqZECZN2+eIiIi5OvrqzZt2mjPnj23bf/3v/9dDRo0kK+vr5o0aaJNmzaVUaXuzZlxXrx4sR599FFVr15d1atXV3R0dLF/L/iBs/+eb1i9erVsNptiY2NLt8B7hLPjnJ2drYSEBIWGhsrHx0e/+MUv+L+jBJwd51mzZql+/fry8/NTeHi4hg8frqtXr5ZRte7p888/V48ePRQWFiabzab169cX+5pt27bpoYceko+Pjx544AEtW7as1OuUVcGsXr3a8vb2tv785z9bX3/9tTVo0CArMDDQysrKKrL9F198YXl6elrTp0+3Dh06ZL355ptWpUqVrIMHD5Zx5e7F2XF+9tlnrXnz5ln79++30tLSrP79+1sBAQHWt99+W8aVuxdnx/mG9PR062c/+5n16KOPWj179iybYt2Ys+Ocl5dntWrVyurWrZu1Y8cOKz093dq2bZuVmppaxpW7F2fHeeXKlZaPj4+1cuVKKz093frkk0+s0NBQa/jw4WVcuXvZtGmT9cYbb1hr1661JFnr1q27bfsTJ05YlStXtkaMGGEdOnTImjNnjuXp6Wlt3ry5VOuscAHl4YcfthISEuzPr1+/boWFhVlTp04tsn2vXr2s7t27O2xr06aNNXjw4FKt0905O843u3btmlWtWjVr+fLlpVXiPeFOxvnatWtW27ZtrT/96U9Wv379CCgl4Ow4L1iwwLr//vut/Pz8sirxnuDsOCckJFidOnVy2DZixAirXbt2pVrnvaQkAWXkyJFW48aNHbb17t3biomJKcXKLKtCneLJz89XSkqKoqOj7ds8PDwUHR2t5OTkIl+TnJzs0F6SYmJifrI97mycb3blyhUVFBQoKCiotMp0e3c6zm+99ZZq166tF198sSzKdHt3Ms4bNmxQVFSUEhISFBwcrAcffFBTpkzR9evXy6pst3Mn49y2bVulpKTYTwOdOHFCmzZtUrdu3cqk5oqivD4Hy+1W9+Xh+++/1/Xr12+5nX5wcLAOHz5c5GsyMzOLbJ+ZmVlqdbq7Oxnnm40aNUphYWG3/FDg/9zJOO/YsUNLlixRampqGVR4b7iTcT5x4oQ+++wz9enTR5s2bdKxY8f00ksvqaCgQOPHjy+Lst3OnYzzs88+q++//17t27eXZVm6du2afvvb3+r1118vi5IrjJ/6HMzNzdV///tf+fn5lcr7VqgZFLiHadOmafXq1Vq3bp18fX3Lu5x7xsWLF/X8889r8eLFqlmzZnmXc08rLCxU7dq1tWjRIrVs2VK9e/fWG2+8oYULF5Z3afeUbdu2acqUKZo/f7727duntWvXauPGjZo0aVJ5lwYXqFAzKDVr1pSnp6eysrIctmdlZSkkJKTI14SEhDjVHnc2zje8++67mjZtmv7xj3+oadOmpVmm23N2nI8fP66TJ0+qR48e9m2FhYWSJC8vLx05ckR169Yt3aLd0J38ew4NDVWlSpXk6elp39awYUNlZmYqPz9f3t7epVqzO7qTcR47dqyef/55DRw4UJLUpEkTXb58WfHx8XrjjTfk4cHv4K7wU5+D/v7+pTZ7IlWwGRRvb2+1bNlSSUlJ9m2FhYVKSkpSVFRUka+JiopyaC9JW7Zs+cn2uLNxlqTp06dr0qRJ2rx5s1q1alUWpbo1Z8e5QYMGOnjwoFJTU+2PJ554Qo899phSU1MVHh5eluW7jTv599yuXTsdO3bMHgAl6ZtvvlFoaCjh5CfcyThfuXLllhByIxRafA+uy5Tb52CpLsE10OrVqy0fHx9r2bJl1qFDh6z4+HgrMDDQyszMtCzLsp5//nlr9OjR9vZffPGF5eXlZb377rtWWlqaNX78eC4zLgFnx3natGmWt7e3tWbNGisjI8P+uHjxYnkdgltwdpxvxlU8JePsOJ8+fdqqVq2alZiYaB05csT66KOPrNq1a1uTJ08ur0NwC86O8/jx461q1apZf/3rX60TJ05Yn376qVW3bl2rV69e5XUIbuHixYvW/v37rf3791uSrBkzZlj79++3Tp06ZVmWZY0ePdp6/vnn7e1vXGb82muvWWlpada8efO4zLi0zJkzx6pTp47l7e1tPfzww9auXbvs+zp06GD169fPof3f/vY36xe/+IXl7e1tNW7c2Nq4cWMZV+yenBnn++67z5J0y2P8+PFlX7ibcfbf848RUErO2XHeuXOn1aZNG8vHx8e6//77rbffftu6du1aGVftfpwZ54KCAmvChAlW3bp1LV9fXys8PNx66aWXrAsXLpR94W5k69atRf5/e2Ns+/XrZ3Xo0OGW1zRv3tzy9va27r//fmvp0qWlXqfNspgHAwAAZqlQa1AAAIB7IKAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHH+H2VC65XSQKkKAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGzCAYAAAAFROyYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAz60lEQVR4nO3de1xVZd7///cG5CRuEA+AE4aZqZSH0lKyxkOOqGRaODrlKKZmGXinlpqTechKx2nUydQmb0fsUY7deauVp3QctSnxhOKYoFOKYV856CSgNoLA+v3Rj323FZWNwL6Q1/Px2I+He61rrfVZl+J+c61rrW2zLMsSAACAQTzcXQAAAMCVCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKEAtMGPGDNlsNp09e/aGbSMiIjR8+PCqLwoAroOAAuCWsmbNGg0ePFh33HGH/P391bJlS7344ovKzc0ts/2nn36q++67T76+vmratKmmT5+uoqKiq9rl5uZq9OjRatSokerWravu3bvrwIEDVXw2QO3l5e4CAJjl2LFj8vCoub+7jB49Wk2aNNFvf/tbNW3aVIcPH9Y777yjjRs36sCBA/Lz83O03bRpkwYMGKBu3bpp4cKFOnz4sF5//XXl5ORoyZIljnYlJSWKiYnRoUOHNHHiRDVs2FCLFy9Wt27dlJycrBYtWrjjVIFbGgEFgBMfHx93l3BTVq9erW7dujkt69Chg+Li4vThhx9q1KhRjuUvvfSS2rZtqy1btsjL66f/Du12u95880298MILatWqlWOfu3bt0scff6yBAwdKkgYNGqS77rpL06dP18qVK6vn5IBapOb+mgTAZWfPntWgQYNkt9vVoEEDvfDCC7p06ZJTmyvnoCQmJspms+mrr77ShAkTHJc4Hn/8cZ05c8Zp2/379ys6OloNGzaUn5+fmjVrphEjRlTHqTlcGU4k6fHHH5ckpaWlOZalpqYqNTVVo0ePdoQTSXr++edlWZZWr17tWLZ69WqFhIToiSeecCxr1KiRBg0apE8++UQFBQVVcCZA7cYIClCLDBo0SBEREZo9e7Z2796tt99+W+fOndP7779/w23Hjh2r+vXra/r06Tp58qQWLFighIQEffTRR5KknJwc9erVS40aNdLLL7+soKAgnTx5UmvWrLnhvi9cuHBVUCpLnTp1FBgYeOMTvUJWVpYkqWHDho5lBw8elCR17NjRqW2TJk102223OdaXtr3vvvuuuvT1wAMP6L333tO//vUvtWnTxuW6AFwbAQWoRZo1a6ZPPvlEkhQfHy+73a7Fixc7LnVcT4MGDbRlyxbZbDZJP83LePvtt5WXl6fAwEDt2rVL586d05YtW5w+9F9//fUb1pWQkKAVK1bcsF3Xrl21Y8eOG7a70u9//3t5eno6Ls9IUmZmpiQpLCzsqvZhYWE6ffq0U9tf/vKXZbaTpNOnTxNQgEpGQAFqkfj4eKf3Y8eO1eLFi7Vx48YbBpTRo0c7wokkPfzww5o/f76+++47tW3bVkFBQZKk9evXq127dqpTp06565o0aZJ++9vf3rBd/fr1y73PUitXrtSyZcs0adIkp8ms//nPfySVPefG19dX+fn5Tm2v1e7n+wJQeQgoQC1y5d0mzZs3l4eHh06ePHnDbZs2ber0vjQsnDt3TtJPoxuxsbGaOXOm5s+fr27dumnAgAF66qmnbjjxNjIyUpGRkS6cSfn84x//0MiRIxUdHa033njDaV3p3TxlzR+5dOmS090+fn5+12z3830BqDwEFKAW+/mIyI14enqWudyyLMe+Vq9erd27d+uzzz7T559/rhEjRuiPf/yjdu/erYCAgGvuOy8vr1yjEN7e3goODi5XvYcOHdJjjz2me+65R6tXr3aaCCv93+WZzMxMhYeHO63LzMzUAw884NS29JLQle2kn+atAKhc3MUD1CLffPON0/tvv/1WJSUlioiIqLRjdO7cWW+88Yb279+vDz/8UEeOHNGqVauuu80LL7ygsLCwG75+fhfN9Rw/fly9e/dW48aNtXHjxjLDUfv27SX9dOfRz50+fVrff/+9Y31p2wMHDqikpMSp7Z49e+Tv76+77rqrXHUBKD9GUIBaZNGiRerVq5fj/cKFCyVJffr0uel9nzt3TkFBQU6jMqUf8je6Dbcy56BkZWWpV69e8vDw0Oeff65GjRqV2e7uu+9Wq1at9N577+nZZ591jBAtWbJENpvNaULtwIEDtXr1aq1Zs8ax/OzZs/r444/Vr1+/Gv/sGMBEBBSgFklPT9djjz2m3r17KykpSR988IGeeuoptWvX7qb3vWLFCi1evFiPP/64mjdvrvPnz2vp0qWy2+3q27fvdbetzDkovXv31okTJzRp0iR9+eWX+vLLLx3rQkJC9Ktf/crx/g9/+IMee+wx9erVS7/5zW/09ddf65133tGoUaPUunVrR7uBAweqc+fOevrpp5Wamup4kmxxcbFmzpxZKXUDuIIF4JY3ffp0S5KVmppqDRw40KpXr55Vv359KyEhwfrPf/7j1Pb222+34uLiHO+XL19uSbL27dvn1G779u2WJGv79u2WZVnWgQMHrCeffNJq2rSp5ePjYzVu3Nh69NFHrf3791f16TmRdM1X165dr2q/du1aq3379paPj4912223WVOnTrUKCwuvavfDDz9YI0eOtBo0aGD5+/tbXbt2vapPAFQem2X9/zPcAAAADMEkWQAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA49TIB7WVlJTo9OnTqlevnkvfJQIAANzHsiydP39eTZo0kYfH9cdIamRAOX369FVf7gUAAGqGU6dO6bbbbrtumxoZUOrVqyfppxO02+1urgYAAJRHfn6+wsPDHZ/j11MjA0rpZR273U5AAQCghinP9AwmyQIAAOMQUAAAgHEIKAAAwDg1cg4KAJjKsiwVFRWpuLjY3aUA1c7T01NeXl6V8ggQAgoAVJLCwkJlZmbqxx9/dHcpgNv4+/srLCxM3t7eN7UfAgoAVIKSkhKlp6fL09NTTZo0kbe3Nw+SRK1iWZYKCwt15swZpaenq0WLFjd8GNv1EFAAoBIUFhaqpKRE4eHh8vf3d3c5gFv4+fmpTp06+u6771RYWChfX98K74tJsgBQiW7mN0bgVlBZPwP8JAEAAOMQUAAAgHGYgwIAVSzi5Q3VeryTc2Kq9XiVZcaMGVqyZIlycnK0du1aDRgwwG212Gw2t9dQ2zGCAgCoMjNmzFD79u1v2C4tLU0zZ87Un//8Z2VmZqpPnz5VX5yuXV911oCyMYICAHC748ePS5L69+9vxO3ZoaGh7i6h1mMEBQBquZKSEs2dO1d33nmnfHx81LRpU73xxhuO9YcPH1aPHj3k5+enBg0aaPTo0bpw4YJj/Y4dO/TAAw+obt26CgoKUpcuXfTdd98pMTFRM2fO1KFDh2Sz2WSz2ZSYmHjV8WfMmKF+/fpJ+ukOkNKA0q1bN40bN86p7YABAzR8+HDH+4iICL355psaMWKE6tWrp6ZNm+q9995z2ub777/Xk08+qeDgYNWtW1cdO3bUnj17rlufzWbTunXryt0Hw4cP14ABA/TWW28pLCxMDRo0UHx8vC5fvuzKXwV+hhGUsmyf7e4KXNd9irsrAFBDTZkyRUuXLtX8+fP10EMPKTMzU0ePHpUkXbx4UdHR0YqKitK+ffuUk5OjUaNGKSEhQYmJiSoqKtKAAQP0zDPP6K9//asKCwu1d+9e2Ww2DR48WF9//bU2b96sv/3tb5KkwMDAq47/0ksvKSIiQk8//bQyMzNdrv+Pf/yjZs2apd/97ndavXq1xowZo65du6ply5a6cOGCunbtql/84hf69NNPFRoaqgMHDqikpKTc9d2oD0pt375dYWFh2r59u7799lsNHjxY7du31zPPPOPyOYGAAgC12vnz5/WnP/1J77zzjuLi4iRJzZs310MPPSRJWrlypS5duqT3339fdevWlSS988476tevn37/+9+rTp06ysvL06OPPqrmzZtLklq3bu3Yf0BAgLy8vK57ySQgIEBBQUGSKnZppW/fvnr++eclSZMnT9b8+fO1fft2tWzZUitXrtSZM2e0b98+BQcHS5LuvPNOl+q7UR+EhIRIkurXr6933nlHnp6eatWqlWJiYrRt2zYCSgVxiQcAarG0tDQVFBTokUceueb6du3aOT6YJalLly4qKSnRsWPHFBwcrOHDhys6Olr9+vXTn/70pwqNgtyMtm3bOv5ss9kUGhqqnJwcSVJKSoruvfdeRzipiBv1Qam7775bnp6ejvdhYWGOOuA6AgoA1GJ+fn43vY/ly5crKSlJDz74oD766CPddddd2r17903v18PDQ5ZlOS0ra05HnTp1nN7bbDaVlJRIqpzzK6/r1QHXEVAAoBZr0aKF/Pz8tG3btjLXt27dWocOHdLFixcdy7766it5eHioZcuWjmX33nuvpkyZol27dumee+7RypUrJUne3t4qLi6uUG2NGjVyGo0pLi7W119/7dI+2rZtq5SUFP3www9lri9PfeXtA1QuAgoA1GK+vr6aPHmyJk2apPfff1/Hjx/X7t27tWzZMknSkCFD5Ovrq7i4OH399dfavn27xo4dq6FDhyokJETp6emaMmWKkpKS9N1332nLli365ptvHPNQIiIilJ6erpSUFJ09e1YFBQXlrq1Hjx7asGGDNmzYoKNHj2rMmDHKzc116fyefPJJhYaGasCAAfrqq6904sQJ/e///q+SkpLKXd+N+gBVg0myAFDFTH+y66uvviovLy9NmzZNp0+fVlhYmJ577jlJkr+/vz7//HO98MILuv/+++Xv76/Y2FjNmzfPsf7o0aNasWKF/v3vfyssLEzx8fF69tlnJUmxsbFas2aNunfvrtzcXC1fvtzpNuHrGTFihA4dOqRhw4bJy8tL48ePV/fu3V06N29vb23ZskUvvvii+vbtq6KiIkVGRmrRokXlru9GfYCqYbOuvMBXA+Tn5yswMFB5eXmy2+2VfwBuMwbgokuXLik9PV3NmjW7qa+YB2q66/0suPL5zSUeAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcHnUPAFWtup9OfQs/Wfq9997TrFmz9P/+3//TvHnzNG7cuCo5zowZM7Ru3TqlpKRIkoYPH67c3FytW7euSo5X00RERGjcuHFV1v8SAQUAUEPk5+crISFB8+bNU2xsrAIDA91dEqoQAQUAUCNkZGTo8uXLiomJUVhYmLvLqfEKCwvl7e3t7jKuiTkoAFDLrV69Wm3atJGfn58aNGignj176uLFi5Kkbt26XTWMP2DAAKdv/I2IiNDrr7+uYcOGKSAgQLfffrs+/fRTnTlzRv3791dAQIDatm2r/fv3X7eOjIwMR3u73a5BgwYpOztbkpSYmKg2bdpIku644w7ZbDadPHnyqn2cPHlSNptNq1at0oMPPihfX1/dc8892rlzp6NNYmKigoKCnLZbt26dbDZbOXvs+n12pXPnzmnIkCFq1KiR/Pz81KJFCy1fvtyxfu/evbr33nvl6+urjh07au3atbLZbI7LS+Wp9/jx4+rfv79CQkIUEBCg+++/X3/729+ctomIiNCsWbM0bNgw2e12jR49WpL05Zdf6uGHH5afn5/Cw8P1X//1X07nkpOTo379+snPz0/NmjXThx9+WO5+uhkEFACoxTIzM/Xkk09qxIgRSktL044dO/TEE0/I1S+6nz9/vrp06aKDBw8qJiZGQ4cO1bBhw/Tb3/5WBw4cUPPmzTVs2LBr7rekpET9+/fXDz/8oJ07d2rr1q06ceKEBg8eLEkaPHiw4wN37969yszMVHh4+DXrmThxol588UUdPHhQUVFR6tevn/7973+7dE7X4mqfvfrqq0pNTdWmTZuUlpamJUuWqGHDhpKkCxcu6NFHH1VkZKSSk5M1Y8YMvfTSSy7XdOHCBfXt21fbtm3TwYMH1bt3b/Xr108ZGRlO7d566y21a9dOBw8e1Kuvvqrjx4+rd+/eio2N1T//+U999NFH+vLLL5WQkODYZvjw4Tp16pS2b9+u1atXa/HixcrJyXG5RldxiQcAarHMzEwVFRXpiSee0O233y5JjpEKV/Tt21fPPvusJGnatGlasmSJ7r//fv3617+WJE2ePFlRUVHKzs5WaGjoVdtv27ZNhw8fVnp6uiN4vP/++7r77ru1b98+3X///WrQoIEkqVGjRmXu4+cSEhIUGxsrSVqyZIk2b96sZcuWadKkSS6f25Vc7bOMjAzde++96tixo6SfRjJKrVy5UiUlJVq2bJl8fX1199136/vvv9eYMWNcqqldu3Zq166d4/2sWbO0du1affrpp05ho0ePHnrxxRcd70eNGqUhQ4Y4RslatGiht99+W127dtWSJUuUkZGhTZs2ae/evbr//vslScuWLVPr1q1dqq8iGEEBgFqsXbt2euSRR9SmTRv9+te/1tKlS3Xu3DmX99O2bVvHn0NCQiQ5f2iXLrvWb95paWkKDw93GhWJjIxUUFCQ0tLSXK4nKirK8WcvLy917NixQvspi6t9NmbMGK1atUrt27fXpEmTtGvXLse6tLQ0tW3bVr6+vmXWXl4XLlzQSy+9pNatWysoKEgBAQFKS0u7agSlNCSVOnTokBITExUQEOB4RUdHq6SkROnp6UpLS5OXl5c6dOjg2KZVq1ZXXXKqCgQUAKjFPD09tXXrVm3atEmRkZFauHChWrZsqfT0dEmSh4fHVZcuLl++fNV+6tSp4/hz6dyIspaVlJRU+jm4qrzndC036rMr9enTR999953Gjx+v06dP65FHHnHpMk556n3ppZe0du1avfnmm/rHP/6hlJQUtWnTRoWFhU7t6tat6/T+woULevbZZ5WSkuJ4HTp0SN98842aN29e7hqrAgEFAGo5m82mLl26aObMmTp48KC8vb21du1aST9dTsnMzHS0LS4u1tdff13pNbRu3VqnTp3SqVOnHMtSU1OVm5uryMhIl/e3e/dux5+LioqUnJzsuCzRqFEjnT9/3mkiaOmE1PK6Xp+VpVGjRoqLi9MHH3ygBQsW6L333pP003n/85//1KVLl8qsvbz1fvXVVxo+fLgef/xxtWnTRqGhoWVOIr7Sfffdp9TUVN15551Xvby9vdWqVStH/5U6duyYcnNzb7jvm0VAAYBabM+ePXrzzTe1f/9+ZWRkaM2aNTpz5ozjw7xHjx7asGGDNmzYoKNHj2rMmDFV8uHUs2dPtWnTRkOGDNGBAwe0d+9eDRs2TF27dr3qskR5LFq0SGvXrtXRo0cVHx+vc+fOacSIEZKkTp06yd/fX7/73e90/PhxrVy5UomJieXe94367ErTpk3TJ598om+//VZHjhzR+vXrHW2feuop2Ww2PfPMM0pNTdXGjRv11ltvOW1fnnpbtGihNWvWOEZAnnrqqXKNVk2ePFm7du1SQkKCUlJS9M033+iTTz5xzFtp2bKlevfurWeffVZ79uxRcnKyRo0aJT8/v3L3V0UxSRYAqprBT3a12+364osvtGDBAuXn5+v222/XH//4R/Xp00eSNGLECB06dEjDhg2Tl5eXxo8fr+7du1d6HTabTZ988onGjh2rX/7yl/Lw8FDv3r21cOHCCu1vzpw5mjNnjlJSUnTnnXfq008/ddw5ExwcrA8++EATJ07U0qVL9cgjj2jGjBmO225v5EZ9diVvb29NmTJFJ0+elJ+fnx5++GGtWrVKkhQQEKDPPvtMzz33nO69915FRkbq97//vWOCb3nrnTdvnkaMGKEHH3xQDRs21OTJk5Wfn3/Dc2nbtq127typV155RQ8//LAsy1Lz5s0dd09J0vLlyzVq1Ch17dpVISEhev311/Xqq6+Wq69uhs1y9V4yA+Tn5yswMFB5eXmy2+2Vf4Dqfix1ZTD4P0CgNrh06ZLS09PVrFkzpwmPqF4nT55Us2bNdPDgQbVv397d5VRITT+H6/0suPL5zSUeAABgHAIKAAAwDnNQAAC3jIiICJefgmuaW+EcKgMjKAAAwDgEFACoRPzmi9qusn4GCCgAUAlKn5r6448/urkSwL1KfwZ+/iThimAOCgBUAk9PTwUFBTm+a8bf39/xeHegNrAsSz/++KNycnIUFBQkT0/Pm9ofAQUAKknpN+xWx1fRA6YKCgq64bdNlwcBBQAqic1mU1hYmBo3buzSl88Bt4o6derc9MhJKQIKAFQyT0/PSvtPGqitmCQLAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADDOTQWUOXPmyGazady4cY5lly5dUnx8vBo0aKCAgADFxsYqOzvbabuMjAzFxMTI399fjRs31sSJE1VUVHQzpQAAgFtIhQPKvn379Oc//1lt27Z1Wj5+/Hh99tln+vjjj7Vz506dPn1aTzzxhGN9cXGxYmJiVFhYqF27dmnFihVKTEzUtGnTKn4WAADgllKhgHLhwgUNGTJES5cuVf369R3L8/LytGzZMs2bN089evRQhw4dtHz5cu3atUu7d++WJG3ZskWpqan64IMP1L59e/Xp00ezZs3SokWLVFhYWObxCgoKlJ+f7/QCAAC3rgoFlPj4eMXExKhnz55Oy5OTk3X58mWn5a1atVLTpk2VlJQkSUpKSlKbNm0UEhLiaBMdHa38/HwdOXKkzOPNnj1bgYGBjld4eHhFygYAADWEywFl1apVOnDggGbPnn3VuqysLHl7eysoKMhpeUhIiLKyshxtfh5OSteXrivLlClTlJeX53idOnXK1bIBAEAN4tKXBZ46dUovvPCCtm7dKl9f36qq6So+Pj7y8fGptuMBAAD3cmkEJTk5WTk5Obrvvvvk5eUlLy8v7dy5U2+//ba8vLwUEhKiwsJC5ebmOm2XnZ2t0NBQSVJoaOhVd/WUvi9tAwAAajeXAsojjzyiw4cPKyUlxfHq2LGjhgwZ4vhznTp1tG3bNsc2x44dU0ZGhqKioiRJUVFROnz4sHJychxttm7dKrvdrsjIyEo6LQAAUJO5dImnXr16uueee5yW1a1bVw0aNHAsHzlypCZMmKDg4GDZ7XaNHTtWUVFR6ty5sySpV69eioyM1NChQzV37lxlZWVp6tSpio+P5zIOAACQ5GJAKY/58+fLw8NDsbGxKigoUHR0tBYvXuxY7+npqfXr12vMmDGKiopS3bp1FRcXp9dee62ySwEAADWUzbIsy91FuCo/P1+BgYHKy8uT3W6v/ANsv/oOJeN1n+LuCgAAuC5XPr/5Lh4AAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOC4FlCVLlqht27ay2+2y2+2KiorSpk2bHOsvXbqk+Ph4NWjQQAEBAYqNjVV2drbTPjIyMhQTEyN/f381btxYEydOVFFRUeWcDQAAuCW4FFBuu+02zZkzR8nJydq/f7969Oih/v3768iRI5Kk8ePH67PPPtPHH3+snTt36vTp03riiScc2xcXFysmJkaFhYXatWuXVqxYocTERE2bNq1yzwoAANRoNsuyrJvZQXBwsP7whz9o4MCBatSokVauXKmBAwdKko4eParWrVsrKSlJnTt31qZNm/Too4/q9OnTCgkJkSS9++67mjx5ss6cOSNvb+9yHTM/P1+BgYHKy8uT3W6/mfLLtn125e+zqnWf4u4KAAC4Llc+vys8B6W4uFirVq3SxYsXFRUVpeTkZF2+fFk9e/Z0tGnVqpWaNm2qpKQkSVJSUpLatGnjCCeSFB0drfz8fMcoTFkKCgqUn5/v9AIAALculwPK4cOHFRAQIB8fHz333HNau3atIiMjlZWVJW9vbwUFBTm1DwkJUVZWliQpKyvLKZyUri9ddy2zZ89WYGCg4xUeHu5q2QAAoAZxOaC0bNlSKSkp2rNnj8aMGaO4uDilpqZWRW0OU6ZMUV5enuN16tSpKj0eAABwLy9XN/D29tadd94pSerQoYP27dunP/3pTxo8eLAKCwuVm5vrNIqSnZ2t0NBQSVJoaKj27t3rtL/Su3xK25TFx8dHPj4+rpYKAABqqJt+DkpJSYkKCgrUoUMH1alTR9u2bXOsO3bsmDIyMhQVFSVJioqK0uHDh5WTk+Nos3XrVtntdkVGRt5sKQAA4Bbh0gjKlClT1KdPHzVt2lTnz5/XypUrtWPHDn3++ecKDAzUyJEjNWHCBAUHB8tut2vs2LGKiopS586dJUm9evVSZGSkhg4dqrlz5yorK0tTp05VfHw8IyQAAMDBpYCSk5OjYcOGKTMzU4GBgWrbtq0+//xz/epXv5IkzZ8/Xx4eHoqNjVVBQYGio6O1ePFix/aenp5av369xowZo6ioKNWtW1dxcXF67bXXKvesAABAjXbTz0FxB56DUgaegwIAMFy1PAcFAACgqhBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxvNxdgIkWbPuXu0tw2bju7q4AAIDKwwgKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxXAoos2fP1v3336969eqpcePGGjBggI4dO+bU5tKlS4qPj1eDBg0UEBCg2NhYZWdnO7XJyMhQTEyM/P391bhxY02cOFFFRUU3fzYAAOCW4FJA2blzp+Lj47V7925t3bpVly9fVq9evXTx4kVHm/Hjx+uzzz7Txx9/rJ07d+r06dN64oknHOuLi4sVExOjwsJC7dq1SytWrFBiYqKmTZtWeWcFAABqNJtlWVZFNz5z5owaN26snTt36pe//KXy8vLUqFEjrVy5UgMHDpQkHT16VK1bt1ZSUpI6d+6sTZs26dFHH9Xp06cVEhIiSXr33Xc1efJknTlzRt7e3jc8bn5+vgIDA5WXlye73V7R8q9pwdSnK32fVW3c68vdXQIAANflyuf3Tc1BycvLkyQFBwdLkpKTk3X58mX17NnT0aZVq1Zq2rSpkpKSJElJSUlq06aNI5xIUnR0tPLz83XkyJEyj1NQUKD8/HynFwAAuHVVOKCUlJRo3Lhx6tKli+655x5JUlZWlry9vRUUFOTUNiQkRFlZWY42Pw8npetL15Vl9uzZCgwMdLzCw8MrWjYAAKgBKhxQ4uPj9fXXX2vVqlWVWU+ZpkyZory8PMfr1KlTVX5MAADgPhX6ssCEhAStX79eX3zxhW677TbH8tDQUBUWFio3N9dpFCU7O1uhoaGONnv37nXaX+ldPqVtruTj4yMfH5+KlAoAAGogl0ZQLMtSQkKC1q5dq7///e9q1qyZ0/oOHTqoTp062rZtm2PZsWPHlJGRoaioKElSVFSUDh8+rJycHEebrVu3ym63KzIy8mbOBQAA3CJcGkGJj4/XypUr9cknn6hevXqOOSOBgYHy8/NTYGCgRo4cqQkTJig4OFh2u11jx45VVFSUOnfuLEnq1auXIiMjNXToUM2dO1dZWVmaOnWq4uPjGSUBAACSXAwoS5YskSR169bNafny5cs1fPhwSdL8+fPl4eGh2NhYFRQUKDo6WosXL3a09fT01Pr16zVmzBhFRUWpbt26iouL02uvvXZzZwIAAG4ZLgWU8jwyxdfXV4sWLdKiRYuu2eb222/Xxo0bXTk0AACoRfguHgAAYBwCCgAAMA4BBQAAGIeAAgAAjFOhB7XBQNtnu7sC13Wf4u4KAACGYgQFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABjHy90FoHIs2PYvd5fgsnHd3V0BAMBUjKAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGMfL3QWg9op4eYO7S3DZyTkx7i4BAGoFRlAAAIBxCCgAAMA4XOIBXMBlKQCoHoygAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGcTmgfPHFF+rXr5+aNGkim82mdevWOa23LEvTpk1TWFiY/Pz81LNnT33zzTdObX744QcNGTJEdrtdQUFBGjlypC5cuHBTJwIAAG4dLgeUixcvql27dlq0aFGZ6+fOnau3335b7777rvbs2aO6desqOjpaly5dcrQZMmSIjhw5oq1bt2r9+vX64osvNHr06IqfBQAAuKW4/Kj7Pn36qE+fPmWusyxLCxYs0NSpU9W/f39J0vvvv6+QkBCtW7dOv/nNb5SWlqbNmzdr37596tixoyRp4cKF6tu3r9566y01adLkJk4HAADcCip1Dkp6erqysrLUs2dPx7LAwEB16tRJSUlJkqSkpCQFBQU5wokk9ezZUx4eHtqzZ0+Z+y0oKFB+fr7TCwAA3Loq9csCs7KyJEkhISFOy0NCQhzrsrKy1LhxY+civLwUHBzsaHOl2bNna+bMmZVZKgwwzmu1u0tw2YKige4uAQBqhRpxF8+UKVOUl5fneJ06dcrdJQEAgCpUqQElNDRUkpSdne20PDs727EuNDRUOTk5TuuLior0ww8/ONpcycfHR3a73ekFAABuXZUaUJo1a6bQ0FBt27bNsSw/P1979uxRVFSUJCkqKkq5ublKTk52tPn73/+ukpISderUqTLLAQAANZTLc1AuXLigb7/91vE+PT1dKSkpCg4OVtOmTTVu3Di9/vrratGihZo1a6ZXX31VTZo00YABAyRJrVu3Vu/evfXMM8/o3Xff1eXLl5WQkKDf/OY33MEDAAAkVSCg7N+/X927d3e8nzBhgiQpLi5OiYmJmjRpki5evKjRo0crNzdXDz30kDZv3ixfX1/HNh9++KESEhL0yCOPyMPDQ7GxsXr77bcr4XQAXCni5Q3uLsFlJ+fEuLsEAG5msyzLcncRrsrPz1dgYKDy8vKqZD7KgqlPV/o+cWvgLp7qQUABbk2ufH7XiLt4AABA7UJAAQAAxqnUB7UBQGVg3gwARlAAAIBxCCgAAMA4BBQAAGAc5qAAQCVg3gxQuRhBAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMw4PaABeM81rt7hJctqBooLtLAACXMYICAACMwwgKcItj1AdATcQICgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcbjNGABqqYiXN7i7BJednBPj7hJQTRhBAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDrcZAwBqDG6Nrj0IKACMM85rtbtLcNmCooHuLgG4pXCJBwAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAONxmDABAFaqJz26R3P/8FkZQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOXxYIAJVgnNdqd5fgsgVFA91dAnBNjKAAAADjEFAAAIBxuMQDALUUl6VgMkZQAACAcQgoAADAOFziAQDUGFyWqj0YQQEAAMYhoAAAAONwiQcAgCpUEy9L/STGrUdnBAUAABiHgAIAAIzj1oCyaNEiRUREyNfXV506ddLevXvdWQ4AADCE2wLKRx99pAkTJmj69Ok6cOCA2rVrp+joaOXk5LirJAAAYAi3BZR58+bpmWee0dNPP63IyEi9++678vf311/+8hd3lQQAAAzhlrt4CgsLlZycrClTpjiWeXh4qGfPnkpKSrqqfUFBgQoKChzv8/LyJEn5+flVUt+lgsIq2S8AADVFVXzGlu7TsqwbtnVLQDl79qyKi4sVEhLitDwkJERHjx69qv3s2bM1c+bMq5aHh4dXWY0AANRmU95aWWX7Pn/+vAIDA6/bpkY8B2XKlCmaMGGC431JSYl++OEHNWjQQDabrVKPlZ+fr/DwcJ06dUp2u71S943/Qz9XD/q5etDP1YN+rj5V1deWZen8+fNq0qTJDdu6JaA0bNhQnp6eys7OdlqenZ2t0NDQq9r7+PjIx8fHaVlQUFBVlii73c4PQDWgn6sH/Vw96OfqQT9Xn6ro6xuNnJRyyyRZb29vdejQQdu2bXMsKykp0bZt2xQVFeWOkgAAgEHcdolnwoQJiouLU8eOHfXAAw9owYIFunjxop5++ml3lQQAAAzhtoAyePBgnTlzRtOmTVNWVpbat2+vzZs3XzVxtrr5+Pho+vTpV11SQuWin6sH/Vw96OfqQT9XHxP62maV514fAACAasR38QAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAME6tDCiLFi1SRESEfH191alTJ+3du/e67T/++GO1atVKvr6+atOmjTZu3FhNldZsrvTz0qVL9fDDD6t+/fqqX7++evbsecO/F/zE1X/PpVatWiWbzaYBAwZUbYG3CFf7OTc3V/Hx8QoLC5OPj4/uuusu/u8oB1f7ecGCBWrZsqX8/PwUHh6u8ePH69KlS9VUbc30xRdfqF+/fmrSpIlsNpvWrVt3w2127Nih++67Tz4+PrrzzjuVmJhY5XXKqmVWrVpleXt7W3/5y1+sI0eOWM8884wVFBRkZWdnl9n+q6++sjw9Pa25c+daqamp1tSpU606depYhw8frubKaxZX+/mpp56yFi1aZB08eNBKS0uzhg8fbgUGBlrff/99NVdes7jaz6XS09OtX/ziF9bDDz9s9e/fv3qKrcFc7eeCggKrY8eOVt++fa0vv/zSSk9Pt3bs2GGlpKRUc+U1i6v9/OGHH1o+Pj7Whx9+aKWnp1uff/65FRYWZo0fP76aK69ZNm7caL3yyivWmjVrLEnW2rVrr9v+xIkTlr+/vzVhwgQrNTXVWrhwoeXp6Wlt3ry5SuusdQHlgQcesOLj4x3vi4uLrSZNmlizZ88us/2gQYOsmJgYp2WdOnWynn322Sqts6ZztZ+vVFRUZNWrV89asWJFVZV4S6hIPxcVFVkPPvig9d///d9WXFwcAaUcXO3nJUuWWHfccYdVWFhYXSXeElzt5/j4eKtHjx5OyyZMmGB16dKlSuu8lZQnoEyaNMm6++67nZYNHjzYio6OrsLKLKtWXeIpLCxUcnKyevbs6Vjm4eGhnj17KikpqcxtkpKSnNpLUnR09DXbo2L9fKUff/xRly9fVnBwcFWVWeNVtJ9fe+01NW7cWCNHjqyOMmu8ivTzp59+qqioKMXHxyskJET33HOP3nzzTRUXF1dX2TVORfr5wQcfVHJysuMy0IkTJ7Rx40b17du3WmquLdz1Oei2R927w9mzZ1VcXHzV4/RDQkJ09OjRMrfJysoqs31WVlaV1VnTVaSfrzR58mQ1adLkqh8K/J+K9POXX36pZcuWKSUlpRoqvDVUpJ9PnDihv//97xoyZIg2btyob7/9Vs8//7wuX76s6dOnV0fZNU5F+vmpp57S2bNn9dBDD8myLBUVFem5557T7373u+oouda41udgfn6+/vOf/8jPz69KjlurRlBQM8yZM0erVq3S2rVr5evr6+5ybhnnz5/X0KFDtXTpUjVs2NDd5dzSSkpK1LhxY7333nvq0KGDBg8erFdeeUXvvvuuu0u7pezYsUNvvvmmFi9erAMHDmjNmjXasGGDZs2a5e7SUAlq1QhKw4YN5enpqezsbKfl2dnZCg0NLXOb0NBQl9qjYv1c6q233tKcOXP0t7/9TW3btq3KMms8V/v5+PHjOnnypPr16+dYVlJSIkny8vLSsWPH1Lx586otugaqyL/nsLAw1alTR56eno5lrVu3VlZWlgoLC+Xt7V2lNddEFennV199VUOHDtWoUaMkSW3atNHFixc1evRovfLKK/Lw4HfwynCtz0G73V5loydSLRtB8fb2VocOHbRt2zbHspKSEm3btk1RUVFlbhMVFeXUXpK2bt16zfaoWD9L0ty5czVr1ixt3rxZHTt2rI5SazRX+7lVq1Y6fPiwUlJSHK/HHntM3bt3V0pKisLDw6uz/BqjIv+eu3Tpom+//dYRACXpX//6l8LCwggn11CRfv7xxx+vCiGlodDie3Arjds+B6t0Cq6BVq1aZfn4+FiJiYlWamqqNXr0aCsoKMjKysqyLMuyhg4dar388suO9l999ZXl5eVlvfXWW1ZaWpo1ffp0bjMuB1f7ec6cOZa3t7e1evVqKzMz0/E6f/68u06hRnC1n6/EXTzl42o/Z2RkWPXq1bMSEhKsY8eOWevXr7caN25svf766+46hRrB1X6ePn26Va9ePeuvf/2rdeLECWvLli1W8+bNrUGDBrnrFGqE8+fPWwcPHrQOHjxoSbLmzZtnHTx40Pruu+8sy7Ksl19+2Ro6dKijfeltxhMnTrTS0tKsRYsWcZtxVVm4cKHVtGlTy9vb23rggQes3bt3O9Z17drViouLc2r/P//zP9Zdd91leXt7W3fffbe1YcOGaq64ZnKln2+//XZL0lWv6dOnV3/hNYyr/55/joBSfq72865du6xOnTpZPj4+1h133GG98cYbVlFRUTVXXfO40s+XL1+2ZsyYYTVv3tzy9fW1wsPDreeff946d+5c9Rdeg2zfvr3M/29L+zYuLs7q2rXrVdu0b9/e8vb2tu644w5r+fLlVV6nzbIYBwMAAGapVXNQAABAzUBAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADj/H/aq+PQYxQUtQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for bins in data:\n", " plt.figure()\n", " plt.title(f\"bins = {bins}\")\n", " plt.hist(chi2(bins-2).cdf(data[bins]), bins=10, range=(0, 1), label=\"cost function\")\n", " plt.hist(chi2(bins-2).cdf(data2[bins]), bins=10, range=(0, 1), alpha=0.5, label=\"sum of pulls squared\")\n", " plt.legend();" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "When 20 bins are used, the number of counts per bin is large enough so that both test statistics are chi-square distributed. When 200 bins are used with samples of the same size, the density in some bins drops low enough so that we are not in the asymptotic limit and see deviations from the theoretical chi-square distribution. These deviations are larger for the sum of pulls squared." ] } ], "metadata": { "kernelspec": { "display_name": "py310", "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.9" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }