{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Automatic differentiation with JAX\n", "\n", "Here we look into automatic differentiation, which can speed up fits with very many parameters.\n", "\n", "iminuit's minimization algorithm MIGRAD uses a mix of gradient descent and Newton's method to find the minimum. Both require a first derivative, which MIGRAD usually computes numerically from finite differences. This requires many function evaluations and the gradient may not be accurate. As an alternative, iminuit also allows the user to compute the gradient and pass it to MIGRAD.\n", "\n", "Although computing derivatives is often straight-forward, it is usually too much hassle to do manually. Automatic differentiation (AD) is an interesting alternative, it allows one to compute exact derivatives efficiently for pure Python/numpy functions. We demonstrate automatic differentiation with the JAX module, which can not only compute derivatives, but also accelerates the computation of Python code (including the gradient code) with a just-in-time compiler.\n", "\n", "[Recommended read: Gentle introduction to AD](https://www.kaggle.com/borisettinger/gentle-introduction-to-automatic-differentiation)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fit of a gaussian model to a histogram\n", "\n", "We fit a gaussian to a histogram using a maximum-likelihood approach based on Poisson statistics. This example is used to investigate how automatic differentiation can accelerate a typical fit in a counting experiment.\n", "\n", "To compare fits with and without passing an analytic gradient fairly, we use `Minuit.strategy = 0`, which prevents Minuit from automatically computing the Hesse matrix after the fit." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:37.436843Z", "start_time": "2020-02-21T10:26:37.432080Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/hdembinski/Extern/iminuit/venv/lib/python3.10/site-packages/jax/_src/lib/__init__.py:33: UserWarning: JAX on Mac ARM machines is experimental and minimally tested. Please see https://github.com/google/jax/issues/5501 in the event of problems.\n", " warnings.warn(\"JAX on Mac ARM machines is experimental and minimally tested. \"\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "JAX version 0.3.2\n", "numba version 0.56.4\n" ] } ], "source": [ "# !pip install jax jaxlib matplotlib numpy iminuit numba-stats\n", "\n", "import jax\n", "from jax import numpy as jnp # replacement for normal numpy\n", "from jax.scipy.special import erf # replacement for scipy.special.erf\n", "from iminuit import Minuit\n", "import numba as nb\n", "import numpy as np # original numpy still needed, since jax does not cover full API\n", "\n", "jax.config.update(\"jax_enable_x64\", True) # enable float64 precision, default is float32\n", "\n", "print(f\"JAX version {jax.__version__}\")\n", "print(f\"numba version {nb.__version__}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We generate some toy data and write the negative log-likelihood (nll) for a fit to binned data, assuming Poisson-distributed counts.\n", "\n", "**Note:** We write all statistical functions in pure Python code, to demonstrate Jax's ability to automatically differentiate and JIT compile this code. In practice, one should import JIT-able statistical distributions from jax.scipy.stats. The library versions can be expected to have fewer bugs and to be faster and more accurate than hand-written code." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:37.594856Z", "start_time": "2020-02-21T10:26:37.585943Z" } }, "outputs": [], "source": [ "# generate some toy data\n", "rng = np.random.default_rng(seed=1)\n", "n, xe = np.histogram(rng.normal(size=10000), bins=1000)\n", "\n", "\n", "def cdf(x, mu, sigma):\n", " # cdf of a normal distribution, needed to compute the expected counts per bin\n", " # better alternative for real code: from jax.scipy.stats.norm import cdf\n", " z = (x - mu) / sigma\n", " return 0.5 * (1 + erf(z / np.sqrt(2)))\n", "\n", "\n", "def nll(par): # negative log-likelihood with constants stripped\n", " amp = par[0]\n", " mu, sigma = par[1:]\n", " p = cdf(xe, mu, sigma)\n", " mu = amp * jnp.diff(p)\n", " result = jnp.sum(mu - n + n * jnp.log(n / (mu + 1e-100) + 1e-100))\n", " return result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check results from all combinations of using JIT and gradient and then compare the execution times." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:37.890967Z", "start_time": "2020-02-21T10:26:37.886224Z" } }, "outputs": [], "source": [ "start_values = (1.5 * np.sum(n), 1.0, 2.0)\n", "limits = ((0, None), (None, None), (0, None))\n", "\n", "\n", "def make_and_run_minuit(fcn, grad=None):\n", " m = Minuit(fcn, start_values, grad=grad, name=(\"amp\", \"mu\", \"sigma\"))\n", " m.errordef = Minuit.LIKELIHOOD\n", " m.limits = limits\n", " m.strategy = 0 # do not explicitly compute hessian after minimisation\n", " m.migrad()\n", " return m" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:38.532308Z", "start_time": "2020-02-21T10:26:38.368563Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n" ] }, { "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 = 496.2 Nfcn = 66
EDM = 1.84e-08 (Goal: 0.0001) time = 0.2 sec
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok APPROXIMATE Pos. def. Not forced
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 66 │\n", "│ EDM = 1.84e-08 (Goal: 0.0001) │ time = 0.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 │APPROXIMATE│ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m1 = make_and_run_minuit(nll)\n", "m1.fmin" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:39.371830Z", "start_time": "2020-02-21T10:26:38.797460Z" } }, "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 = 496.2 Nfcn = 26, Ngrad = 6
EDM = 1.84e-08 (Goal: 0.0001) time = 0.5 sec
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok APPROXIMATE Pos. def. Not forced
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 26, Ngrad = 6 │\n", "│ EDM = 1.84e-08 (Goal: 0.0001) │ time = 0.5 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 │APPROXIMATE│ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m2 = make_and_run_minuit(nll, grad=jax.grad(nll))\n", "m2.fmin" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:39.510553Z", "start_time": "2020-02-21T10:26:39.373728Z" } }, "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 = 496.2 Nfcn = 26, Ngrad = 6
EDM = 1.88e-08 (Goal: 0.0001)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok APPROXIMATE Pos. def. Not forced
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 26, Ngrad = 6 │\n", "│ EDM = 1.88e-08 (Goal: 0.0001) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ No Parameters at limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Below EDM threshold (goal x 10) │ Below call limit │\n", "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n", "│ Covariance │ Hesse ok │APPROXIMATE│ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m3 = make_and_run_minuit(jax.jit(nll), grad=jax.grad(nll))\n", "m3.fmin" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:40.573574Z", "start_time": "2020-02-21T10:26:40.229476Z" } }, "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 = 496.2 Nfcn = 26, Ngrad = 6
EDM = 1.88e-08 (Goal: 0.0001) time = 0.1 sec
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok APPROXIMATE Pos. def. Not forced
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 26, Ngrad = 6 │\n", "│ EDM = 1.88e-08 (Goal: 0.0001) │ time = 0.1 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 │APPROXIMATE│ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m4 = make_and_run_minuit(jax.jit(nll), grad=jax.jit(jax.grad(nll)))\n", "m4.fmin" ] }, { "cell_type": "code", "execution_count": 8, "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 = 496.2 Nfcn = 82
EDM = 5.31e-05 (Goal: 0.0001) time = 0.9 sec
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok APPROXIMATE Pos. def. Not forced
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 82 │\n", "│ EDM = 5.31e-05 (Goal: 0.0001) │ time = 0.9 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 │APPROXIMATE│ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from numba_stats import norm # numba jit-able version of norm\n", "\n", "@nb.njit\n", "def nb_nll(par):\n", " amp = par[0]\n", " mu, sigma = par[1:]\n", " p = norm.cdf(xe, mu, sigma)\n", " mu = amp * np.diff(p)\n", " result = np.sum(mu - n + n * np.log(n / (mu + 1e-323) + 1e-323))\n", " return result\n", "\n", "m5 = make_and_run_minuit(nb_nll)\n", "m5.fmin" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:45.031931Z", "start_time": "2020-02-21T10:26:40.674388Z" } }, "outputs": [], "source": [ "from timeit import timeit\n", "\n", "times = {\n", " \"no JIT, no grad\": \"m1\",\n", " \"no JIT, grad\": \"m2\",\n", " \"jax JIT, no grad\": \"m3\",\n", " \"jax JIT, grad\": \"m4\",\n", " \"numba JIT, no grad\": \"m5\",\n", "}\n", "for k, v in times.items():\n", " t = timeit(\n", " f\"{v}.values = start_values; {v}.migrad()\",\n", " f\"from __main__ import {v}, start_values\",\n", " number=1,\n", " )\n", " times[k] = t" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:45.142272Z", "start_time": "2020-02-21T10:26:45.033451Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqAAAAGwCAYAAAB2AtfDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGcUlEQVR4nO3deXhN597/8c9OIpFBEmJItDGWiLHmprSiQgxV9NfyKJK0xoMTaqynWpQmFKecDlpFQlGlh1bLoaJCG8RQUSlSFKEnOKbEVCRZvz9c9mNXRJCsHfJ+Xde+aq11r3t97zuhn6wpFsMwDAEAAAAmcbB3AQAAAChaCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoCoxhGMrIyBCvmgUAALcigKLAXLhwQV5eXrpw4YK9SwEAAIUIARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqQigAAAAMBUBFAAAAKYigAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwlZO9C8Cjr/a4tXJwcbN3GTDRkckd7F0CAKAQ4wwoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqQigAAAAMBUBFAAAAKYigAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCoCKAAAAExFAAUAAICpCKAACkSlSpVksVhu+wwaNCjH9suXL1ejRo3k7e0td3d3Pfnkk/r8889NrhoAYAYnexcA4NG0fft2ZWVlWZeTk5PVunVrvfzyyzm2L1WqlN58803VqFFDzs7O+u677/Tqq6+qbNmyCg0NNatsAIAJOAMKoECUKVNGvr6+1s93332nqlWrqkWLFjm2Dw4OVpcuXRQYGKiqVatqyJAhqlu3rn766SdJ0v79++Xm5qbFixdb91m6dKlcXV21d+9eU8YEAMgfBFAABe7atWtauHChXnvtNVkslru2NwxD69evV0pKip599llJUo0aNTRt2jQNHDhQqampOn78uAYMGKApU6aoZs2aBT0EAEA+shiGYdi7CDwcLBaLVqxYoc6dO+epfUZGhry8vOQ/dKkcXNwKtjgUKkcmd7BZXrp0qV555RWlpqaqfPnyd9wvPT1djz32mK5evSpHR0d9/PHHeu2112zaPP/888rIyJCzs7McHR21Zs2aPIVaAEDhUeTvAf1rqLq5fP78eb366qu57nv48GFVqlSp4IsEHnJz585Vu3btcg2fklSiRAklJSXp4sWLWr9+vYYNG6YqVaooODjY2mbevHmqXr26HBwc9OuvvxI+AeAhVOQD6J1069ZNbdu2tS6/+OKLql27tt555x3rujJlytijtPuWlZUli8UiBwfuvIB5jh49qri4OC1fvvyubR0cHPTEE09Ikp588knt27dP0dHRNgF09+7dunTpkhwcHJSWliY/P7+CKh0AUEAemiQSHBysyMhIjRo1SqVKlZKvr6/Gjx9v0yY1NVWdOnWSh4eHPD091bVrV508efK+jufq6mrzAIWzs7Pc3Nxs1jk6Ouapr0qVKikqKkqvvfaaSpQooQoVKmj27Nk2bfbs2aPnnntOrq6u8vHxUb9+/XTx4sVc+125cqWqVaum4sWLq2XLlpo/f74sFovOnz8vSYqNjZW3t7dWrlypmjVrysXFRampqdq+fbtat26t0qVLy8vLSy1atNDPP/9s0/eBAwf07LPPqnjx4qpZs6bWrVuX98kDbhETE6OyZcuqQ4cOd2/8F9nZ2bp69ap1+ezZs4qIiNCbb76piIgI9ejRQ1euXMnPcgEAJnhoAqgkzZ8/X+7u7kpMTNR7772nd955xxqMsrOz1alTJ509e1YbN27UunXr9Pvvv6tbt252rvqG6dOnq1GjRtq1a5cGDhyov/3tb0pJSZEkXbp0SaGhoSpZsqS2b9+uZcuWKS4uToMHD75jf4cPH9ZLL72kzp07a/fu3erfv7/efPPN29pdvnxZU6ZM0Zw5c/Trr7+qbNmyunDhgsLDw/XTTz9p69atqlatmtq3b68LFy5IujGXL774opydnZWYmKhPPvlEo0ePvusYr169qoyMDJsPirbs7GzFxMQoPDxcTk62F1zCwsI0ZswY63J0dLT17+2+ffs0ffp0ff755+rZs6e1zYABA+Tv76+xY8fqH//4h7KysjRixAjTxgMAyB8P1SX4unXraty4cZKkatWq6cMPP9T69evVunVrrV+/Xnv27NHhw4fl7+8vSVqwYIFq1aql7du3q3HjxvYsXe3bt9fAgQMlSaNHj9b777+vDRs2KCAgQIsXL9aff/6pBQsWyN3dXZL04YcfqmPHjpoyZYrKlSt3W3+ffvqpAgICNHXqVElSQECAkpOT9e6779q0u379uj7++GPVq1fPuu65556zaTN79mx5e3tr48aNev755xUXF6f9+/dr7dq11nv2oqKi1K5du1zHGB0drQkTJtzjzOBRFhcXp9TU1NseJJJuXLG49XaQS5cuaeDAgTp+/LhcXV1Vo0YNLVy40PpD5IIFC7R69Wrt2rVLTk5OcnJy0sKFC9W8eXM9//zzd/3+BAAUHg9dAL2Vn5+fTp06JUnat2+f/P39reFTkmrWrClvb2/t27fP7gH01totFot8fX1taq9Xr541fEpSs2bNlJ2drZSUlBwDaEpKym1jatKkyW3tnJ2db5u3kydPauzYsYqPj9epU6eUlZWly5cvKzU11VqPv7+/zQMjQUFBdx3jmDFjNGzYMOtyRkaGzdcDRU+bNm10pxdtxMfH2yxPmjRJkyZNumNfYWFhCgsLs1nXpEkTXbt27YHrBACY66EKoMWKFbNZtlgsys7OtlM198Zetbu6ut72lHB4eLjOnDmjmTNnqmLFinJxcVFQUNAD/4/cxcVFLi4uD9QHAAB49D1U94DmJjAwUMeOHdOxY8es6/bu3avz588X+pdUBwYGWp/svSkhIUEODg4KCAjIcZ+AgADt2LHDZt327dvzdLyEhARFRkaqffv2qlWrllxcXHT69Gmbeo4dO6a0tDTruq1bt97LkAAAAO7okQmgISEhqlOnjnr06KGff/5Z27ZtU1hYmFq0aKFGjRrZu7xc9ejRQ8WLF1d4eLiSk5O1YcMG/f3vf1evXr1yvPwuSf3799f+/fs1evRo/fbbb1q6dKliY2Ml6a7vRaxWrZo+//xz7du3T4mJierRo4dcXV2t20NCQlS9enWFh4dr9+7d+vHHH3N8wAkAAOB+PDIB1GKx6JtvvlHJkiX17LPPKiQkRFWqVNGXX355x31uXgL/69O59yI+Pl4Wi0VHjhy57z7c3Ny0du1anT17Vo0bN9ZLL72kVq1a6cMPP7zjPpUrV9ZXX32l5cuXq27dupo1a5Y1JN7tMvjcuXN17tw5NWjQQL169VJkZKTKli1r3e7g4KAVK1boypUratKkifr06XPbw00AAAD3q0j/Ks4TJ07Iz89P27dvv++zpDExMYqKitLevXtvu8/TbO+++64++eQTm9sQ7IlfxVl0/fVXcQIAcKuH6iGk/GIYho4ePapp06apXLlyql279n33tXr1akVFRdklfH788cdq3LixfHx8lJCQoKlTp+b67lAAAIDCoEgG0PT0dAUEBCgwMFBLlixR8eLF77uvZcuW5WNl9+bAgQOaNGmSzp49qwoVKmj48OE2L/YGAAAojIr0JXgULC7BF11cggcA5OaReQgJAAAADwcCKAAAAExFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMZTEMw7B3EXg0ZWRkyMvLS+np6fL09LR3OQAAoJDgDCgAAABMRQAFAACAqQigAAAAMBUBFAAAAKYigAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFRO9i4Aj77a49bKwcXN3mUAAPDIODK5g71LeCCcAQUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqQigAAAAMBUBFAAA4CEUHR2txo0bq0SJEipbtqw6d+6slJQU6/YjR47IYrHk+Fm2bFmejjFgwABZLBbNmDEjX2sngAIAADyENm7cqEGDBmnr1q1at26drl+/rjZt2ujSpUuSJH9/f6Wlpdl8JkyYIA8PD7Vr1+6u/a9YsUJbt25V+fLl8712AigAAMBDaM2aNYqIiFCtWrVUr149xcbGKjU1VTt37pQkOTo6ytfX1+azYsUKde3aVR4eHrn2/ccff+jvf/+7Fi1apGLFitlsW7BggTw8PHTgwAHruoEDB6pGjRq6fPlynmongAIAADwC0tPTJUmlSpXKcfvOnTuVlJSk3r1759pPdna2evXqpZEjR6pWrVq3bQ8LC1P79u3Vo0cPZWZmatWqVZozZ44WLVokNze3PNVKAEWeRUREqHPnzvYuAwAA/EV2draGDh2qZs2aqXbt2jm2mTt3rgIDA/X000/n2teUKVPk5OSkyMjIO7b59NNPlZaWpsjISPXu3Vvjx49Xw4YN81wvAbSAWCwWff3117ctx8bG3vGG4JufI0eO2K1uAADw8Bk0aJCSk5O1ZMmSHLdfuXJFixcvvuvZz507d2rmzJnWvHInJUuW1Ny5czVr1ixVrVpVb7zxxj3VSwA1Wbdu3WxuBg4KClLfvn1t1vn7+xfY8a9fv15gfQMAAPMNHjxY3333nTZs2KDHH388xzZfffWVLl++rLCwsFz7+vHHH3Xq1ClVqFBBTk5OcnJy0tGjRzV8+HBVqlTJpu2mTZvk6OiotLQ064NPeVXkA2hwcLAiIyM1atQolSpVSr6+vho/frxNm9TUVHXq1EkeHh7y9PRU165ddfLkyfs6nqurq83NwM7OznJzc7NZ5+jomKe+0tLS1KFDB7m6uqpy5cpavHixKlWqZPOqBIvFolmzZumFF16Qu7u73n33XWVlZal3796qXLmyXF1dFRAQoJkzZ9r0nZWVpWHDhsnb21s+Pj4aNWqUDMO4rzEDAID8ZxiGBg8erBUrVuiHH35Q5cqV79h27ty5euGFF1SmTJlc++zVq5d++eUXJSUlWT/ly5fXyJEjtXbtWmu7zZs3a8qUKfr222/l4eGhwYMH31PtTvfU+hE1f/58DRs2TImJidqyZYsiIiLUrFkztW7dWtnZ2dbwuXHjRmVmZmrQoEHq1q2b4uPj7Vp3WFiYTp8+rfj4eBUrVkzDhg3TqVOnbms3fvx4TZ48WTNmzJCTk5Oys7P1+OOPa9myZfLx8dHmzZvVr18/+fn5qWvXrpKk6dOnKzY2VvPmzVNgYKCmT5+uFStW6LnnnrtjPVevXtXVq1etyxkZGfk/aAAAIOnGZffFixfrm2++UYkSJXTixAlJkpeXl1xdXa3tDh48qE2bNmn16tU59lOjRg1FR0erS5cu8vHxkY+Pj832YsWKydfXVwEBAZKkCxcuqFevXoqMjFS7du30+OOPq3HjxurYsaNeeumlPNVOAJVUt25djRs3TpJUrVo1ffjhh1q/fr1at26t9evXa8+ePTp8+LD10viCBQtUq1Ytbd++XY0bN7ZLzfv371dcXJy2b9+uRo0aSZLmzJmjatWq3db2lVde0auvvmqzbsKECdY/V65cWVu2bNHSpUutAXTGjBkaM2aMXnzxRUnSJ598YvOTT06io6Nt+gUAAAVn1qxZkm5czb1VTEyMIiIirMvz5s3T448/rjZt2uTYT0pKivUJ+rwYMmSI3N3dFRUVJUmqU6eOoqKi1L9/fwUFBemxxx67ax8EUN0IoLfy8/Oznknct2+f/P39be7LrFmzpry9vbVv3z67BdCUlBQ5OTmpQYMG1nVPPPGESpYseVvbmwH1Vh999JHmzZun1NRUXblyRdeuXdOTTz4p6cZrHNLS0tS0aVNreycnJzVq1CjXy/BjxozRsGHDrMsZGRkFej8rAABFWV5vjYuKirKGxfvp568PR8+bN++2NsOGDbPJAHdDAJVue8GqxWJRdna2narJf+7u7jbLS5Ys0YgRIzR9+nQFBQWpRIkSmjp1qhITEx/oOC4uLnJxcXmgPgAAwKOvyD+EdDeBgYE6duyYjh07Zl23d+9enT9/XjVr1rRbXQEBAcrMzNSuXbus6w4ePKhz587ddd+EhAQ9/fTTGjhwoOrXr68nnnhChw4dsm738vKSn5+fTSDNzMy0/mYFAACAB0EAvYuQkBDVqVNHPXr00M8//6xt27YpLCxMLVq0yPHStllq1KihkJAQ9evXT9u2bdOuXbvUr18/ubq65vreLunGfa47duzQ2rVr9dtvv+mtt97S9u3bbdoMGTJEkydP1tdff639+/dr4MCBOn/+fAGOCAAAFBUE0LuwWCz65ptvVLJkST377LMKCQlRlSpV9OWXX95xn5uX752c7v8Oh/j4+Lu+lH7BggUqV66cnn32WXXp0kV9+/ZViRIlVLx48Vz77t+/v1588UV169ZNTZs21ZkzZzRw4ECbNsOHD1evXr0UHh5uvUzfpUuX+x4PAADATRaDlzvmuxMnTsjPz8/mCfV7FRMTo6ioKO3du/e2e1Tv5Pjx4/L391dcXJxatWp1X8fNTxkZGfLy8pL/0KVycMnb74YFAAB3d2RyB3uX8EB4CCkfGYaho0ePatq0aSpXrtwdfxdrXqxevVpRUVG5hs8ffvhBFy9eVJ06dZSWlqZRo0apUqVKevbZZ+/7uAAAAAWNAJqP0tPTFRAQoMDAQC1ZsuSul8Jzs2zZsru2uX79uv73f/9Xv//+u0qUKKGnn35aixYtyvMZUwAAAHsggOYjb29vm98EVNBCQ0MVGhpq2vEAAADyAw8hAQAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqQigAAAAMBUBFAAAAKYigAIAAMBUBFAAAACYymIYhmHvIvBoysjIkJeXl9LT0+Xp6WnvcgAAQCHBGVAAAACYigAKAAAAUxFAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKmc7F0AHn21x62Vg4ubvcsAACDPjkzuYO8SHmmcAQUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqQigAAAAMBUBFAAAIAfR0dFq3LixSpQoobJly6pz585KSUmxaRMcHCyLxWLzGTBgQK79Ll++XG3atJGPj48sFouSkpIKcBSFEwEUAAAgBxs3btSgQYO0detWrVu3TtevX1ebNm106dIlm3Z9+/ZVWlqa9fPee+/l2u+lS5fUvHlzTZkypSDLL9Sc7F0AAABAYbRmzRqb5djYWJUtW1Y7d+7Us88+a13v5uYmX1/fPPfbq1cvSdKRI0dy3B4fH682bdpo/fr1euaZZyRJ7733nqZNm6Y9e/aoXLly9ziSwoczoAAAAHmQnp4uSSpVqpTN+kWLFql06dKqXbu2xowZo8uXLz/QcYKDgzV06FD16tVL6enp2rVrl9566y3NmTPnkQifUiEPoBEREercubO9y4Bu/NTn7e1t7zIAALCL7OxsDR06VM2aNVPt2rWt61955RUtXLhQGzZs0JgxY/T555+rZ8+eD3y8SZMmqWTJkurXr5969uyp8PBwvfDCCw/cb2FRqC/Bz5w5U4ZhFFj/sbGxGjp0qM6fP3/bcnBwsDZu3HjHfVu0aKH4+PgCqw0AABQegwYNUnJysn766Seb9f369bP+uU6dOvLz81OrVq106NAhVa1a9b6P5+zsrEWLFqlu3bqqWLGi3n///fvuqzAq1AHUy8vLbsdevny5rl27Jkk6duyYmjRpori4ONWqVUvSjW+Mh821a9ceyroBALCnwYMH67vvvtOmTZv0+OOP59q2adOmkqSDBw8+UACVpM2bN0uSzp49q7Nnz8rd3f2B+itMHppL8GvWrFHz5s3l7e0tHx8fPf/88zp06JC17YIFC+Th4aEDBw5Y1w0cOFA1atS4r3sxSpUqJV9fX/n6+qpMmTKSJB8fH+u6v97/cSfx8fGyWCxav369GjVqJDc3Nz399NO3vcZh1qxZqlq1qpydnRUQEKDPP/88134zMzMVGRlpnY/Ro0crPDzc5paF4OBgDR48WEOHDlXp0qUVGhoqSfrHP/6hOnXqyN3dXf7+/ho4cKAuXrxo039sbKwqVKggNzc3denSRWfOnMnTeAEAeFQYhqHBgwdrxYoV+uGHH1S5cuW77nPzlUp+fn4PdOxDhw7p9ddf12effaamTZsqPDxc2dnZD9RnYVKoA+itLl26pGHDhmnHjh1av369HBwc1KVLF+sXIywsTO3bt1ePHj2UmZmpVatWac6cOVq0aJHc3NzsXL305ptvavr06dqxY4ecnJz02muvWbetWLFCQ4YM0fDhw5WcnKz+/fvr1Vdf1YYNG+7Y35QpU7Ro0SLFxMQoISFBGRkZ+vrrr29rN3/+fDk7OyshIUGffPKJJMnBwUH//Oc/9euvv2r+/Pn64YcfNGrUKOs+iYmJ6t27twYPHqykpCS1bNlSkyZNuusYr169qoyMDJsPAAAPq0GDBmnhwoVavHixSpQooRMnTujEiRO6cuWKpBshceLEidq5c6eOHDmilStXKiwsTM8++6zq1q1r7adGjRpasWKFdfns2bNKSkrS3r17JUkpKSlKSkrSiRMnJElZWVnq2bOnQkND9eqrryomJka//PKLpk+fbuLoC5bFKMibLB9QRESEzp8/n2OwOn36tMqUKaM9e/ZYbwY+d+6c6tatq44dO2r58uWKjIzU//7v/96x/9zuAb3VkSNHVLlyZe3atUtPPvnkPY0hPj5eLVu2VFxcnFq1aiVJWr16tTp06KArV66oePHiatasmWrVqqXZs2db9+vatasuXbqkVatW5divr6+vRowYoREjRki68c1apUoV1a9f3zpfwcHBysjI0M8//5xrjV999ZUGDBig06dPS7pxQ3V6errNsf/nf/5Ha9asuW1ubjV+/HhNmDDhtvX+Q5fKwcX+PwQAAJBXRyZ3kMViyXFbTEyMIiIidOzYMfXs2VPJycm6dOmS/P391aVLF40dO1aenp7W9haLxbqPdCNvvPrqq7f1O27cOI0fP17vvPOOPvnkE+3Zs0c+Pj6Sbtwa2L17d23btk316tXL/wGbrFDfA3qrAwcO6O2331ZiYqJOnz5tPfOZmppqDaAlS5bU3LlzFRoaqqefflpvvPGGPUu2cetPQjdPy586dUoVKlTQvn37bG5ilqRmzZpp5syZOfaVnp6ukydPqkmTJtZ1jo6Oatiw4W2n5xs2bHjb/nFxcYqOjtb+/fuVkZGhzMxM/fnnn7p8+bLc3Ny0b98+denSxWafoKCg296H9ldjxozRsGHDrMsZGRny9/fPdR8AAAqru52j8/f3z/WB5Tv1ExERYQ2jOXn77bf19ttv26x78cUXdfXq1bse62Hx0FyC79ixo86ePavPPvtMiYmJSkxMlCTrg0I3bdq0SY6OjkpLS7vtNxXYU7Fixax/vvkTlRn3cvz1huUjR47o+eefV926dfWvf/1LO3fu1EcffSTp9rm8Vy4uLvL09LT5AAAA/NVDEUDPnDmjlJQUjR07Vq1atVJgYKDOnTt3W7vNmzdrypQp+vbbb+Xh4aHBgwfbodp7FxgYqISEBJt1CQkJqlmzZo7tvby8VK5cOW3fvt26Lisr666X2iVp586dys7O1vTp0/XUU0+pevXq+s9//nNbPTcD/k1bt27N63AAAABy9VBcgi9ZsqR8fHw0e/Zs+fn5KTU19bbL6xcuXFCvXr0UGRmpdu3a6fHHH1fjxo3VsWNHvfTSS3aqPG9Gjhyprl27qn79+goJCdG3336r5cuXKy4u7o77/P3vf1d0dLSeeOIJ1ahRQx988IHOnTt3x/tVbnriiSd0/fp1ffDBB+rYsaPNw0k3RUZGqlmzZpo2bZo6deqktWvX3vXyOwAAQF49FGdAHRwctGTJEu3cuVO1a9fW66+/rqlTp9q0GTJkiNzd3RUVFSXpxstgo6Ki1L9/f/3xxx859pudnS0npwfL4BEREQoODn6gPjp37qyZM2dq2rRpqlWrlj799FPFxMTk2u/o0aPVvXt3hYWFKSgoSB4eHgoNDVXx4sVzPVa9evX0j3/8Q1OmTFHt2rW1aNEiRUdH27R56qmn9Nlnn2nmzJmqV6+evv/+e40dO/aBxggAAHBToX4Kvnv37nJ0dNTChQsLpP/Jkydr4cKFSk5Ovu8+WrRooZYtW2r8+PH5V9h9yM7OVmBgoLp27aqJEyfatZabMjIy5OXlxVPwAICHzpHJHexdwiOtUF6Cz8zM1G+//aYtW7aof//++d7/5cuXtX//fsXExKhdu3b33U96eroOHTp0x1clFaSjR4/q+++/V4sWLXT16lV9+OGHOnz4sF555RXTawEAALgXhfISfHJysho1aqRatWppwIAB+d7/7NmzFRISonr16t32moN74eXlpePHj8vDwyMfq8sbBwcHxcbGqnHjxmrWrJn27NmjuLg4BQYGml4LAADAvSjUl+DxcOMSPADgYcUl+IJVKM+AAgAA4NFFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqSyGYRj2LgKPpoyMDHl5eSk9PV2enp72LgcAABQSnAEFAACAqQigAAAAMBUBFAAAAKYigAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJjKyd4F4NFXe9xaObi45UtfRyZ3yJd+AACA/XAGFAAAAKYigAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQPHQ2bdqkjh07qnz58rJYLPr6669zbb98+XK1bt1aZcqUkaenp4KCgrR27VpzigUAALchgOKhc+nSJdWrV08fffRRntpv2rRJrVu31urVq7Vz5061bNlSHTt21K5duwq4UgAAkBMCKB467dq106RJk9SlS5c8tZ8xY4ZGjRqlxo0bq1q1aoqKilK1atX07bffSpL++9//ytfXV1FRUdZ9Nm/eLGdnZ61fv75AxgAAQFHmZO8CALNlZ2frwoULKlWqlCSpTJkymjdvnjp37qw2bdooICBAvXr10uDBg9WqVSs7VwsAwKOnSJ4BjYiIUOfOne1dxkMnL/dbPgymTZumixcvqmvXrtZ17du3V9++fdWjRw8NGDBA7u7uio6OtmOVAAA8uopkAJ05c6ZiY2MLrP/Y2Fh5e3vnuBwcHCyLxXLHT3BwcIHVBWnx4sWaMGGCli5dqrJly9psmzZtmjIzM7Vs2TItWrRILi4udqoSAIBHW5G8BO/l5WW3Yy9fvlzXrl2TJB07dkxNmjRRXFycatWqJUlydnYusGNnZWXJYrHIwaFI/tyhJUuWqE+fPlq2bJlCQkJu237o0CH95z//UXZ2to4cOaI6derYoUoAAB59RTKJ3HoJfs2aNWrevLm8vb3l4+Oj559/XocOHbK2XbBggTw8PHTgwAHruoEDB6pGjRq6fPnyPR+7VKlS8vX1la+vr8qUKSNJ8vHxsa67eV9iXqxcuVLVqlVT8eLF1bJlS82fP18Wi0Xnz5+X9H9nXleuXKmaNWvKxcVFqamp2r59u1q3bq3SpUvLy8tLLVq00M8//2zT94EDB/Tss8+qePHiqlmzptatW3fPYy1MvvjiC7366qv64osv1KFDh9u2X7t2TT179lS3bt00ceJE9enTR6dOnbJDpQAAPPqKZAC91aVLlzRs2DDt2LFD69evl4ODg7p06aLs7GxJUlhYmNq3b68ePXooMzNTq1at0pw5c7Ro0SK5ubnZre7Dhw/rpZdeUufOnbV79271799fb7755m3tLl++rClTpmjOnDn69ddfVbZsWV24cEHh4eH66aeftHXrVlWrVk3t27fXhQsXJN14SOfFF1+Us7OzEhMT9cknn2j06NF3renq1avKyMiw+RSEixcvKikpSUlJSZJuzEVSUpJSU1MlSWPGjFFYWJi1/eLFixUWFqbp06eradOmOnHihE6cOKH09HRrmzfffFPp6en65z//qdGjR6t69ep67bXXCqR+AACKuiJ5Cf5W/+///T+b5Xnz5qlMmTLau3evateuLUn69NNPVbduXUVGRmr58uUaP368GjZsaI9yrT799FMFBARo6tSpkqSAgAAlJyfr3XfftWl3/fp1ffzxx6pXr5513XPPPWfTZvbs2fL29tbGjRv1/PPPKy4uTvv379fatWtVvnx5SVJUVJTatWuXa03R0dGaMGFCfgwvVzt27FDLli2ty8OGDZMkhYeHKzY2VmlpadYwKt0YX2ZmpgYNGqRBgwZZ199sHx8frxkzZmjDhg3y9PSUJH3++eeqV6+eZs2apb/97W8FPiYAAIqSIh9ADxw4oLfffluJiYk6ffq09cxnamqqNYCWLFlSc+fOVWhoqJ5++mm98cYb9ixZkpSSkqLGjRvbrGvSpMlt7ZydnVW3bl2bdSdPntTYsWMVHx+vU6dOKSsrS5cvX7aGtn379snf398aPiUpKCjorjWNGTPGGgYlKSMjQ/7+/vc0rrwIDg6WYRh33P7XB8zi4+Pv2t/169dt1lWqVMnmDCkAAMg/RT6AduzYURUrVtRnn32m8uXLKzs7W7Vr17Y+KHTTpk2b5OjoqLS0NF26dEklSpSwU8X3xtXVVRaLxWZdeHi4zpw5o5kzZ6pixYpycXFRUFDQbWO+Vy4uLjw5DgAA7qpI3wN65swZpaSkaOzYsWrVqpUCAwN17ty529pt3rxZU6ZM0bfffisPDw8NHjzYDtXaCggI0I4dO2zWbd++PU/7JiQkKDIyUu3bt1etWrXk4uKi06dPW7cHBgbq2LFjSktLs67bunVr/hQOAACKvCIdQEuWLCkfHx/Nnj1bBw8e1A8//GBzCVmSLly4oF69eikyMlLt2rXTokWL9OWXX+qrr76yU9U39O/fX/v379fo0aP122+/aenSpdZLz3894/lX1apV0+eff659+/YpMTFRPXr0kKurq3V7SEiIqlevrvDwcO3evVs//vhjjg84AQAA3I8iHUAdHBy0ZMkS7dy5U7Vr19brr79ufajnpiFDhsjd3d36e8Lr1KmjqKgo9e/fX3/88UeO/WZnZ8vJ6cHuboiIiMj1pfSVK1fWV199peXLl6tu3bqaNWuWNSTe7TL43Llzde7cOTVo0MAarm99KbuDg4NWrFihK1euqEmTJurTp89tDzcBAADcL4uR29Mcj6ju3bvL0dFRCxcuLJD+J0+erIULFyo5Ofm++2jRooVatmyp8ePH53mfd999V5988omOHTt238fNTxkZGfLy8pL/0KVycMmfV1YdmXz7OzwBAMDDpUg9hJSZmanffvtNW7ZsUf/+/fO9/8uXL2v//v2KiYm56yuLcpOenq5Dhw5p1apVubb7+OOP1bhxY/n4+CghIUFTp04tFPenAgAA5KZIXYJPTk5Wo0aNVKtWLQ0YMCDf+589e7ZCQkJUr149vf322/fdj5eXl44fPy4PD49c2x04cECdOnVSzZo1NXHiRA0fPvyezpgCAADYQ5G8BA9zcAkeAADkpEidAQUAAID9EUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqQigAAAAMBUBFAAAAKYigAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCqLYRiGvYvAoykjI0NeXl5KT0+Xp6envcsBAACFBGdAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmcrJ3AXj01R63Vg4ubve9/5HJHfKxGgAAYG+cAQUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqQigAAAAMBUBFA+FTZs2qWPHjipfvrwsFou+/vrru+4THx+vBg0ayMXFRU888YRiY2MLvE4AAHB3BFA8FC5duqR69erpo48+ylP7w4cPq0OHDmrZsqWSkpI0dOhQ9enTR2vXri3gSgEAwN0QQPFQaNeunSZNmqQuXbrkqf0nn3yiypUra/r06QoMDNTgwYP10ksv6f3335ck/fe//5Wvr6+ioqKs+2zevFnOzs5av359gYwBAADcQADFI2nLli0KCQmxWRcaGqotW7ZIksqUKaN58+Zp/Pjx2rFjhy5cuKBevXpp8ODBatWqlT1KBgCgyHhoAmhwcLCGDh1q7zKKtIfpa3DixAmVK1fOZl25cuWUkZGhK1euSJLat2+vvn37qkePHhowYIDc3d0VHR1tj3IBAChSHpoAWlD+GqpuLh85ckQWiyXXDw+1PPymTZumzMxMLVu2TIsWLZKLi4u9SwIA4JHnZO8CCit/f3+lpaVZl6dNm6Y1a9YoLi7Ous7Ly8sepT2Q69evq1ixYvYuo8D5+vrq5MmTNutOnjwpT09Pubq6WtcdOnRI//nPf5Sdna0jR46oTp06ZpcKAECRc09nQIODgxUZGalRo0apVKlS8vX11fjx463bb541TEpKsq47f/68LBaL4uPjJd14NY7FYtHatWtVv359ubq66rnnntOpU6f073//W4GBgfL09NQrr7yiy5cv2xw/MzNTgwcPlpeXl0qXLq233npLhmFYt3/++edq1KiRSpQoIV9fX73yyis6derUvc+KJEdHR/n6+lo/Hh4ecnJysll3a5DJTUREhDp37qxp06bJz89PPj4+GjRokK5fv25tc+7cOYWFhalkyZJyc3NTu3btdODAgVz73b9/v5o3b67ixYurZs2aiouLs3lF0c2vx5dffqkWLVqoePHiWrRokc6cOaPu3bvrsccek5ubm+rUqaMvvvjCpu9Lly4pLCxMHh4e8vPz0/Tp0+9tAu0sKCjotoeJ1q1bp6CgIOvytWvX1LNnT3Xr1k0TJ05Unz597vv7BQAA5N09X4KfP3++3N3dlZiYqPfee0/vvPOO1q1bd88HHj9+vD788ENt3rxZx44dU9euXTVjxgwtXrxYq1at0vfff68PPvjgtmM7OTlp27Ztmjlzpv7xj39ozpw51u3Xr1/XxIkTtXv3bn399dc6cuSIIiIi7rm2grBhwwYdOnRIGzZs0Pz58xUbG2tzCT8iIkI7duzQypUrtWXLFhmGofbt29uE1FtlZWWpc+fOcnNzU2JiombPnq0333wzx7ZvvPGGhgwZon379ik0NFR//vmnGjZsqFWrVik5OVn9+vVTr169tG3bNus+I0eO1MaNG/XNN9/o+++/V3x8vH7++edcx3j16lVlZGTYfPLLxYsXlZSUZP3h5vDhw0pKSlJqaqokacyYMQoLC7O2HzBggH7//XeNGjVK+/fv18cff6ylS5fq9ddft7Z58803lZ6ern/+858aPXq0qlevrtdeey3fagYAAHdg3IMWLVoYzZs3t1nXuHFjY/To0YZhGMbhw4cNScauXbus28+dO2dIMjZs2GAYhmFs2LDBkGTExcVZ20RHRxuSjEOHDlnX9e/f3wgNDbU5dmBgoJGdnW1dN3r0aCMwMPCO9W7fvt2QZFy4cCHXMQ0ZMuSOyzeNGzfOqFev3h37yU14eLhRsWJFIzMz07ru5ZdfNrp162YYhmH89ttvhiQjISHBuv306dOGq6ursXTp0hz7/Pe//204OTkZaWlp1nXr1q0zJBkrVqwwDOP/vh4zZsy4a40dOnQwhg8fbhiGYVy4cMFwdna2OfaZM2cMV1fXHOfmpnHjxhmSbvv4D11qVBz93X1/DOP/vm/++gkPD7fOcYsWLWzq2bBhg/Hkk08azs7ORpUqVYyYmBibbU5OTsaPP/5oXXf48GHD09PT+Pjjj+86XwAA4P7d8z2gdevWtVn28/O7r8uWt/ZTrlw5ubm5qUqVKjbrbj0jJ0lPPfWULBaLdTkoKEjTp09XVlaWHB0dtXPnTo0fP167d+/WuXPnlJ2dLUlKTU1VzZo177nG/FSrVi05Ojpal/38/LRnzx5J0r59++Tk5KSmTZtat/v4+CggIED79u3Lsb+UlBT5+/vL19fXuq5JkyY5tm3UqJHNclZWlqKiorR06VL98ccfunbtmq5evSo3NzdJN+6LvHbtmk09pUqVUkBAQK5jHDNmjIYNG2ZdzsjIkL+/f6775FVwcLDN7RZ/ldMDYcHBwdq1a9cd+/vr2eVKlSopPT39geoEAAB3d88B9K8PsFgsFmvQc3C4cUX/1qBwp0vIt/ZjsVhy7TcvLl26pNDQUIWGhmrRokUqU6aMUlNTFRoaqmvXruW5n4LyoON7EO7u7jbLU6dO1cyZMzVjxgzVqVNH7u7uGjp06APPk4uLC0+RAwCAu8rX1zCVKVNGkmyeHr/1gaQHlZiYaLO8detWVatWTY6Ojtq/f7/OnDmjyZMn65lnnlGNGjUemgdKAgMDlZmZaTO+M2fOKCUl5Y5nbgMCAnTs2DGbJ723b9+ep+MlJCSoU6dO6tmzp+rVq6cqVarot99+s26vWrWqihUrZlPPuXPnbNoAAADcr3wNoK6urnrqqac0efJk7du3Txs3btTYsWPzrf/U1FQNGzZMKSkp+uKLL/TBBx9oyJAhkqQKFSrI2dlZH3zwgX7//XetXLlSEydOzLdjF6Rq1aqpU6dO6tu3r3766Sft3r1bPXv21GOPPaZOnTrluE/r1q1VtWpVhYeH65dfflFCQoJ1rm+9TeFOx1u3bp02b96sffv2qX///jZB1sPDQ71799bIkSP1ww8/KDk5WREREdYz3AAAAA8i3xPFvHnzlJmZqYYNG2ro0KGaNGlSvvUdFhamK1euqEmTJho0aJCGDBmifv36Sbpx9jU2NlbLli1TzZo1NXnyZE2bNu2ufWZnZ8vJ6cFeh5ofL6WPiYlRw4YN9fzzzysoKEiGYWj16tV3fGeno6Ojvv76a128eFGNGzdWnz59rE/BFy9ePNdjjR07Vg0aNFBoaKiCg4Pl6+urzp0727SZOnWqnnnmGXXs2FEhISFq3ry5GjZs+EBjBAAAkCSLkduTHUVAjRo11KdPH40YMeK+9j98+LCqV6+uvXv3qlq1avlc3b1JSEhQ8+bNdfDgQVWtWtWutUg3HkLy8vKS/9ClcnBxu+9+jkzukI9VAQAAeyuyvwnp5ovvU1JS1KpVq/vuZ/Xq1erXr59dwueKFSvk4eGhatWq6eDBgxoyZIiaNWtWKMInAADAnRTZANq2bVudO3dO//znP1W/fv377mfQoEH5WNW9uXDhgkaPHq3U1FSVLl1aISEhD91vLAIAAEVPkb8Ej4LDJXgAAJATHmsGAACAqQigAAAAMBUBFAAAAKYigAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFQEUAAAAJiKAAoAAABTEUABAABgKgIoAAAATEUABQAAgKkIoAAAADCVxTAMw95F4NGUkZEhLy8vpaeny9PT097lAACAQoIzoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAAABMRQAFAACAqQigAAAAMBUBFAAAAKYigAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCoCKAAAAEzlZO8C8OgyDEOSlJGRYedKAADAvSpRooQsFkuB9E0ARYE5c+aMJMnf39/OlQAAgHuVnp4uT0/PAumbAIoCU6pUKUlSamqqvLy87FzNwy0jI0P+/v46duxYgf1jUJQwn/mHucw/zGX+Yj4fXIkSJQqsbwIoCoyDw41bjL28vPjLn088PT2Zy3zEfOYf5jL/MJf5i/ksnHgICQAAAKYigAIAAMBUBFAUGBcXF40bN04uLi72LuWhx1zmL+Yz/zCX+Ye5zF/MZ+FmMW6+KwcAAAAwAWdAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwFQEUefbRRx+pUqVKKl68uJo2bapt27bl2n7ZsmWqUaOGihcvrjp16mj16tU22w3D0Ntvvy0/Pz+5uroqJCREBw4cKMghFCr5PZ/Lly9XmzZt5OPjI4vFoqSkpAKsvnDJz7m8fv26Ro8erTp16sjd3V3ly5dXWFiY/vOf/xT0MAqN/P7eHD9+vGrUqCF3d3eVLFlSISEhSkxMLMghFBr5PZe3GjBggCwWi2bMmJHPVRdO+T2XERERslgsNp+2bdsW5BBwKwPIgyVLlhjOzs7GvHnzjF9//dXo27ev4e3tbZw8eTLH9gkJCYajo6Px3nvvGXv37jXGjh1rFCtWzNizZ4+1zeTJkw0vLy/j66+/Nnbv3m288MILRuXKlY0rV66YNSy7KYj5XLBggTFhwgTjs88+MyQZu3btMmk09pXfc3n+/HkjJCTE+PLLL439+/cbW7ZsMZo0aWI0bNjQzGHZTUF8by5atMhYt26dcejQISM5Odno3bu34enpaZw6dcqsYdlFQczlTcuXLzfq1atnlC9f3nj//fcLeCT2VxBzGR4ebrRt29ZIS0uzfs6ePWvWkIo8AijypEmTJsagQYOsy1lZWUb58uWN6OjoHNt37drV6NChg826pk2bGv379zcMwzCys7MNX19fY+rUqdbt58+fN1xcXIwvvviiAEZQuOT3fN7q8OHDRSqAFuRc3rRt2zZDknH06NH8KboQM2M+09PTDUlGXFxc/hRdSBXUXB4/ftx47LHHjOTkZKNixYpFIoAWxFyGh4cbnTp1KpB6cXdcgsddXbt2TTt37lRISIh1nYODg0JCQrRly5Yc99myZYtNe0kKDQ21tj98+LBOnDhh08bLy0tNmza9Y5+PioKYz6LKrLlMT0+XxWKRt7d3vtRdWJkxn9euXdPs2bPl5eWlevXq5V/xhUxBzWV2drZ69eqlkSNHqlatWgVTfCFTkN+X8fHxKlu2rAICAvS3v/1NZ86cyf8BIEcEUNzV6dOnlZWVpXLlytmsL1eunE6cOJHjPidOnMi1/c3/3kufj4qCmM+iyoy5/PPPPzV69Gh1795dnp6e+VN4IVWQ8/ndd9/Jw8NDxYsX1/vvv69169apdOnS+TuAQqSg5nLKlClycnJSZGRk/hddSBXUXLZt21YLFizQ+vXrNWXKFG3cuFHt2rVTVlZW/g8Ct3GydwEAUFhdv35dXbt2lWEYmjVrlr3Leai1bNlSSUlJOn36tD777DN17dpViYmJKlu2rL1Le2js3LlTM2fO1M8//yyLxWLvch56//M//2P9c506dVS3bl1VrVpV8fHxatWqlR0rKxo4A4q7Kl26tBwdHXXy5Emb9SdPnpSvr2+O+/j6+uba/uZ/76XPR0VBzGdRVZBzeTN8Hj16VOvWrXvkz35KBTuf7u7ueuKJJ/TUU09p7ty5cnJy0ty5c/N3AIVIQczljz/+qFOnTqlChQpycnKSk5OTjh49quHDh6tSpUoFMo7CwKx/M6tUqaLSpUvr4MGDD1407ooAirtydnZWw4YNtX79euu67OxsrV+/XkFBQTnuExQUZNNektatW2dtX7lyZfn6+tq0ycjIUGJi4h37fFQUxHwWVQU1lzfD54EDBxQXFycfH5+CGUAhY+b3ZnZ2tq5evfrgRRdSBTGXvXr10i+//KKkpCTrp3z58ho5cqTWrl1bcIOxM7O+L48fP64zZ87Iz88vfwpH7uz9FBQeDkuWLDFcXFyM2NhYY+/evUa/fv0Mb29v48SJE4ZhGEavXr2MN954w9o+ISHBcHJyMqZNm2bs27fPGDduXI6vYfL29ja++eYb45dffjE6depUpF7DlN/zeebMGWPXrl3GqlWrDEnGkiVLjF27dhlpaWmmj89M+T2X165dM1544QXj8ccfN5KSkmxe0XL16lW7jNFM+T2fFy9eNMaMGWNs2bLFOHLkiLFjxw7j1VdfNVxcXIzk5GS7jNEsBfH3/K+KylPw+T2XFy5cMEaMGGFs2bLFOHz4sBEXF2c0aNDAqFatmvHnn3/aZYxFDQEUefbBBx8YFSpUMJydnY0mTZoYW7dutW5r0aKFER4ebtN+6dKlRvXq1Q1nZ2ejVq1axqpVq2y2Z2dnG2+99ZZRrlw5w8XFxWjVqpWRkpJixlAKhfyez5iYGEPSbZ9x48aZMBr7ys+5vPkaq5w+GzZsMGlE9pWf83nlyhWjS5cuRvny5Q1nZ2fDz8/PeOGFF4xt27aZNRy7yu+/539VVAKoYeTvXF6+fNlo06aNUaZMGaNYsWJGxYoVjb59+1oDLQqexTAMwz7nXgEAAFAUcQ8oAAAATEUABQAAgKkIoAAAADAVARQAAACmIoACAADAVARQAAAAmIoACgAAAFMRQAEAAGAqAigAFGHx8fGyWCw6f/58kTo2APsigAJAEREcHKyhQ4farHv66aeVlpYmLy+vR/bYualcubLi4uLsdnygqHKydwEAAPtxdnaWr69vkTu2JP3yyy86d+6cWrRoYbcagKKKM6AAkA+ys7MVHR2typUry9XVVfXq1dNXX30lSTIMQyEhIQoNDZVhGJKks2fP6vHHH9fbb79t7WPOnDkKDAxU8eLFVaNGDX388cc2xzh+/Li6d++uUqVKyd3dXY0aNVJiYqIkKSIiQp07d7ZpP3ToUAUHB1u3b9y4UTNnzpTFYpHFYtGRI0dyvAz+r3/9S7Vq1ZKLi4sqVaqk6dOn2/RbqVIlRUVF6bXXXlOJEiVUoUIFzZ49+45zk9djx8bGytvbW999950CAgLk5uaml156SZcvX9b8+fNVqVIllSxZUpGRkcrKyrL2f/XqVY0YMUKPPfaY3N3d1bRpU8XHx9/1a/bNN9+obdu2Klas2G3bDMPQ+PHjVaFCBbm4uKh8+fKKjIy8a58A8sgAADywSZMmGTVq1DDWrFljHDp0yIiJiTFcXFyM+Ph4wzAM4/jx40bJkiWNGTNmGIZhGC+//LLRpEkT4/r164ZhGMbChQsNPz8/41//+pfx+++/G//617+MUqVKGbGxsYZhGMaFCxeMKlWqGM8884zx448/GgcOHDC+/PJLY/PmzYZhGEZ4eLjRqVMnm5qGDBlitGjRwjAMwzh//rwRFBRk9O3b10hLSzPS0tKMzMxMY8OGDYYk49y5c4ZhGMaOHTsMBwcH45133jFSUlKMmJgYw9XV1YiJibH2W7FiRaNUqVLGRx99ZBw4cMCIjo42HBwcjP379+c4N3k9dkxMjFGsWDGjdevWxs8//2xs3LjR8PHxMdq0aWN07drV+PXXX41vv/3WcHZ2NpYsWWLtv0+fPsbTTz9tbNq0yTh48KAxdepUw8XFxfjtt99y/Zo1atTIWLx4cY7bli1bZnh6ehqrV682jh49aiQmJhqzZ8/OtT8AeUcABYAH9Oeffxpubm7WMHhT7969je7du1uXly5dahQvXtx44403DHd3d5uAVLVq1dvC0MSJE42goCDDMAzj008/NUqUKGGcOXMmxxruFkANwzBatGhhDBkyxKbNX0PgK6+8YrRu3dqmzciRI42aNWtalytWrGj07NnTupydnW2ULVvWmDVrVo615fXYMTExhiTj4MGD1jb9+/c33NzcjAsXLljXhYaGGv379zcMwzCOHj1qODo6Gn/88YdN361atTLGjBlzx3qOHz9uODs7W4/9V9OnTzeqV69uXLt27Y59ALh/3AMKAA/o4MGDunz5slq3bm2z/tq1a6pfv751+eWXX9aKFSs0efJkzZo1S9WqVZMkXbp0SYcOHVLv3r3Vt29fa/vMzEzrAzpJSUmqX7++SpUqVaBj2bdvnzp16mSzrlmzZpoxY4aysrLk6OgoSapbt651u8Vika+vr06dOvXAx3dzc1PVqlWty+XKlVOlSpXk4eFhs+7msfbs2aOsrCxVr17dpp+rV6/Kx8fnjsdZuXKlmjdvLm9v7xy3v/zyy5oxY4aqVKmitm3bqn379urYsaOcnPjfJpAf+JsEAA/o4sWLkqRVq1bpscces9nm4uJi/fPly5e1c+dOOTo66sCBA7ft/9lnn6lp06Y2+98MfK6urrnW4ODgYL2/9Kbr16/f40jy7q/3TVosFmVnZxdIv7kd6+LFi3J0dLTO661uDa1/tXLlSr3wwgt33O7v76+UlBTFxcVp3bp1GjhwoKZOnaqNGzfmeM8ogHtDAAWAB1SzZk25uLgoNTU11yeqhw8fLgcHB/373/9W+/bt1aFDBz333HMqV66cypcvr99//109evTIcd+6detqzpw5Onv2bI5nQcuUKaPk5GSbdUlJSTZhydnZ2ebhnZwEBgYqISHBZl1CQoKqV69+W8C7F3k59v2oX7++srKydOrUKT3zzDN52ufixYvasGGDZs2alWs7V1dXdezYUR07dtSgQYNUo0YN7dmzRw0aNMiP0oEijQAKAA+oRIkSGjFihF5//XVlZ2erefPmSk9PV0JCgjw9PRUeHq5Vq1Zp3rx52rJlixo0aKCRI0cqPDxcv/zyi0qWLKkJEyYoMjJSXl5eatu2ra5evaodO3bo3LlzGjZsmLp3766oqCh17txZ0dHR8vPz065du1S+fHkFBQXpueee09SpU7VgwQIFBQVp4cKFSk5OtrkFoFKlSkpMTNSRI0fk4eGRY5AdPny4GjdurIkTJ6pbt27asmWLPvzww9ueyL9XeTn2/ahevbp69OihsLAwTZ8+XfXr19d///tfrV+/XnXr1lWHDh1u22fNmjWqXr26KlWqdMd+Y2NjlZWVpaZNm8rNzU0LFy6Uq6urKlasmC91A0Udr2ECgHwwceJEvfXWW4qOjlZgYKDatm2rVatWqXLlyvrvf/+r3r17a/z48dazZxMmTFC5cuU0YMAASVKfPn00Z84cxcTEqE6dOmrRooViY2NVuXJlSTfOIH7//fcqW7as2rdvrzp16mjy5MnWs5KhoaF66623NGrUKDVu3FgXLlxQWFiYTY0jRoyQo6OjatasqTJlyig1NfW2cTRo0EBLly7VkiVLVLt2bb399tt65513FBER8UDzk5dj36+YmBiFhYVp+PDhCggIUOfOnbV9+3ZVqFAhx/bffPNNrpffJcnb21ufffaZmjVrprp16youLk7ffvttrveVAsg7i/HXm4YAAHhEZWZmqly5cvr3v/+tJk2a2LscoMjiDCgAoMg4e/asXn/9dTVu3NjepQBFGmdAAQAAYCrOgAIAAMBUBFAAAACYigAKAAAAUxFAAQAAYCoCKAAAAExFAAUAAICpCKAAAAAwFQEUAAAApiKAAgAAwFT/H83kOhc7/cJwAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot as plt\n", "\n", "x = np.fromiter(times.values(), dtype=float)\n", "xmin = np.min(x)\n", "\n", "y = -np.arange(len(times))\n", "plt.barh(y, x)\n", "for yi, k, v in zip(y, times, x):\n", " plt.text(v, yi, f\"{v/xmin:.1f}x\")\n", "plt.yticks(y, times.keys())\n", "for loc in (\"top\", \"right\"):\n", " plt.gca().spines[loc].set_visible(False)\n", "plt.xlabel(\"execution time / s\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Conclusions:\n", "\n", "1. As expected, the best results are obtained by JIT compiling the function and the gradient.\n", "\n", "2. JIT compiling the cost function with Jax but not using the gradient gives a negligible performance improvement. Numba is able to do much better.\n", "\n", "3. JIT compiling the gradient is very important. Using the Python-computed gradient even drastically reduces performance in this example.\n", "\n", "In general, the gain from using a gradient is larger for functions with hundreds of parameters, as is common in machine learning. Human-made models often have less than 10 parameters, and then the gain is not so dramatic. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computing covariance matrices with JAX\n", "\n", "Automatic differentiation gives us another way to compute uncertainties of fitted parameters. MINUIT compute the uncertainties with the HESSE algorithm by default, which computes the matrix of second derivates approximately using finite differences and inverts this.\n", "\n", "Let's compare the output of HESSE with the exact (within floating point precision) computation using automatic differentiation." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:27:38.715871Z", "start_time": "2020-02-21T10:27:37.907690Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sigma[amp] : HESSE = 100.0, JAX = 100.0\n", "sigma[mu] : HESSE = 0.0100, JAX = 0.0100\n", "sigma[sigma]: HESSE = 0.0071, JAX = 0.0071\n" ] } ], "source": [ "m4.hesse()\n", "cov_hesse = m4.covariance\n", "\n", "\n", "def jax_covariance(par):\n", " return jnp.linalg.inv(jax.hessian(nll)(par))\n", "\n", "\n", "par = np.array(m4.values)\n", "cov_jax = jax_covariance(par)\n", "\n", "print(\n", " f\"sigma[amp] : HESSE = {cov_hesse[0, 0] ** 0.5:6.1f}, JAX = {cov_jax[0, 0] ** 0.5:6.1f}\"\n", ")\n", "print(\n", " f\"sigma[mu] : HESSE = {cov_hesse[1, 1] ** 0.5:6.4f}, JAX = {cov_jax[1, 1] ** 0.5:6.4f}\"\n", ")\n", "print(\n", " f\"sigma[sigma]: HESSE = {cov_hesse[2, 2] ** 0.5:6.4f}, JAX = {cov_jax[2, 2] ** 0.5:6.4f}\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Success, HESSE and JAX give the same answer within the relevant precision.\n", "\n", "**Note:** If you compute the covariance matrix in this way from a least-squares cost function instead of a negative log-likelihood, you must multiply it by 2.\n", "\n", "Let us compare the performance of HESSE with Jax." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6.59 ms ± 595 µs per loop (mean ± std. dev. of 3 runs, 1 loop each)\n" ] } ], "source": [ "%%timeit -n 1 -r 3\n", "m = Minuit(nll, par)\n", "m.errordef = Minuit.LIKELIHOOD\n", "m.hesse()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "14.8 ms ± 523 µs per loop (mean ± std. dev. of 3 runs, 1 loop each)\n" ] } ], "source": [ "%%timeit -n 1 -r 3\n", "jax_covariance(par)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The computation with Jax is slower, but it is also more accurate (although the added precision is not relevant).\n", "\n", "Minuit's HESSE algorithm still makes sense today. It has the advantage that it can process any function, while Jax cannot. Jax cannot differentiate a function that calls into C/C++ code or Cython code, for example.\n", "\n", "Final note: If we JIT compile `jax_covariance`, it greatly outperforms Minuit's HESSE algorithm, but that only makes sense if you need to compute the hessian at different parameter values, so that the extra time spend to compile is balanced by the time saved over many invocations. This is not what happens here, the Hessian in only needed at the best fit point." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "104 µs ± 12.8 µs per loop (mean ± std. dev. of 3 runs, 1 loop each)\n" ] } ], "source": [ "%%timeit -n 1 -r 3 jit_jax_covariance = jax.jit(jax_covariance); jit_jax_covariance(par)\n", "jit_jax_covariance(par)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is much faster... but only because the compilation cost is excluded here." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "285 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n" ] } ], "source": [ "%%timeit -n 1 -r 1\n", "# if we include the JIT compilation cost, the performance drops dramatically\n", "@jax.jit\n", "def jax_covariance(par):\n", " return jnp.linalg.inv(jax.hessian(nll)(par))\n", "\n", "\n", "jax_covariance(par)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With compilation cost included, it is much slower.\n", "\n", "Conclusion: Using the JIT compiler makes a lot of sense if the covariance matrix has to be computed repeatedly for the same cost function but different parameters, but this is not the case when we use it to compute parameter errors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fit data points with uncertainties in x and y\n", "\n", "Let's say we have some data points $(x_i \\pm \\sigma_{x,i}, y_i \\pm \\sigma_{y,i})$ and we have a model $y=f(x)$ that we want to adapt to this data. If $\\sigma_{x,i}$ was zero, we could use the usual least-squares method, minimizing the sum of squared residuals $r^2_i = (y_i - f(x_i))^2 / \\sigma^2_{y,i}$. Here, we don't know where to evaluate $f(x)$, since the exact $x$-location is only known up to $\\sigma_{x,i}$.\n", "\n", "We can approximately extend the standard least-squares method to handle this case. We use that the uncertainty along the $x$-axis can be converted into an additional uncertainty along the $y$-axis with error propagation,\n", "\n", "$$\n", "f(x_i \\pm \\sigma_{x,i}) \\simeq f(x_i) \\pm f'(x_i)\\,\\sigma_{x,i}.\n", "$$\n", "\n", "Using this, we obtain modified squared residuals\n", "\n", "$$\n", "r^2_i = \\frac{(y_i - f(x_i))^2}{\\sigma^2_{y,i} + (f'(x_i) \\,\\sigma_{x,i})^2}.\n", "$$\n", "\n", "We demonstrate this with a fit of a polynomial." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:25:43.510168Z", "start_time": "2020-02-21T10:25:43.371319Z" } }, "outputs": [], "source": [ "# polynomial model\n", "def f(x, par):\n", " return jnp.polyval(par, x)\n", "\n", "\n", "# true polynomial f(x) = x^2 + 2 x + 3\n", "par_true = np.array((1, 2, 3))\n", "\n", "\n", "# grad computes derivative with respect to the first argument\n", "f_prime = jax.jit(jax.grad(f))\n", "\n", "\n", "# checking first derivative f'(x) = 2 x + 2\n", "assert f_prime(0.0, par_true) == 2\n", "assert f_prime(1.0, par_true) == 4\n", "assert f_prime(2.0, par_true) == 6\n", "# ok!\n", "\n", "# generate toy data\n", "n = 30\n", "data_x = np.linspace(-4, 7, n)\n", "data_y = f(data_x, par_true)\n", "\n", "rng = np.random.default_rng(seed=1)\n", "sigma_x = 0.5\n", "sigma_y = 5\n", "data_x += rng.normal(0, sigma_x, n)\n", "data_y += rng.normal(0, sigma_y, n)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:25:43.646212Z", "start_time": "2020-02-21T10:25:43.512384Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0J0lEQVR4nO3dfXSU9Z3//9ckkAkIMxhIMmG5iwUb02CR+xiwFqOBL/LVhdXaA1uCHD3NL1Ihu6ukqwYoEnAtsLoY1GUTu8jBulsVSg1irOnxW5BIyh5DBHULhZVMsKXMIGwmkMzvD5oxAwnJDLnmumbyfJwz5+S6rplr3lwkmVc+d5fN7/f7BQAAYEFxZhcAAADQGYIKAACwLIIKAACwLIIKAACwLIIKAACwLIIKAACwLIIKAACwLIIKAACwrD5mF3CtWltbdfLkSQ0cOFA2m83scgAAQDf4/X6dPXtWQ4cOVVxc5+0mUR9UTp48qeHDh5tdBgAACMOJEyc0bNiwTo9HfVAZOHCgpEv/UIfDYXI1AACgO7xer4YPHx74HO9M1AeVtu4eh8NBUAEAIMp0NWyDwbQAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyDA0qo0aNks1mu+JRWFgoSWpqalJhYaEGDx6sAQMGaN68eWpsbDSyJAAAEEUMDSo1NTVqaGgIPPbs2SNJuu+++yRJy5Yt086dO/X666+rurpaJ0+e1Ny5c40sCQAARBGb3+/3R+rNli5dql/+8pf67LPP5PV6lZycrG3btulv/uZvJEmHDx/WTTfdpL1792rq1KndOqfX65XT6ZTH42FlWgAAokR3P78jNkalublZW7du1YMPPiibzaYDBw7owoULys3NDTwnIyNDI0aM0N69ezs9j8/nk9frDXoAAIDYFLGg8uabb+rMmTPKz8+XJLndbiUkJGjQoEFBz0tNTZXb7e70PKWlpXI6nYEHd04GACB2RSyobNmyRbNmzdLQoUOv6TzFxcXyeDyBx4kTJ3qoQgAAYDURuXvyH/7wB7377rv6xS9+EdjncrnU3NysM2fOBLWqNDY2yuVydXouu90uu91uZLkAAESt880XlfnUbklS/ao89U+IyEe9YSLSolJeXq6UlBTNnj07sG/ChAnq27evqqqqAvuOHDmi48ePKzs7OxJlAQAAizM8ZrW2tqq8vFwLFy5Unz5fv53T6dTixYtVVFSkpKQkORwOLVmyRNnZ2d2e8QMAAGKb4UHl3Xff1fHjx/Xggw9ecWzDhg2Ki4vTvHnz5PP5lJeXpxdeeMHokgAAQJSI6DoqRmAdFQAAvhYtY1Qst44KAABAqAgqAADAsggqAADAsggqAADAsggqAADAsggqAADAsggqAADAsggqAADAsggqAADAsggqAADAsggqAADAsggqAADEkJbWr2/ht//o6aDtaERQAQAgRlTWNSh3fXVgO7+8RtPWvafKugYTq7o2BBUAAGJAZV2DCrbWqtHrC9rv9jSpYGtt1IYVa977GQCAXuh888WwXtfS6lfJjkPqqJPHL8kmacWOeuWMHqL4OFtI5+6fYG5UIKgAAGARmU/tNuS8fklub5PGrngn5NceWzu75wsKAV0/AADAsmhRAQDAIupX5YX1uv1HTyu/vKbL51UsmqTJ6UlhvYdZCCoAAFhEuONBpo9JVpozUW5PU4fjVGySXM5ETR+THPIYFbPR9QMAQJSLj7OpZE6mpEuhpL227ZI5mVEXUiSCCgAAMWFmVprKFoxXisMetN/lTFTZgvGamZVmUmXXhq4fAABixMysNOWMHhKY3VOxaFJUdve0R4sKAAAxpH0omZyeFNUhRSKoAAAACyOoAAAAyyKoAAAAyyKoAAAAyyKoAAAAyyKoAAAAyyKoAAAAyyKoAAAAyzI8qHzxxRdasGCBBg8erH79+mns2LH66KOPAsf9fr+eeuoppaWlqV+/fsrNzdVnn31mdFkAACAKGBpU/vznPysnJ0d9+/bV22+/rfr6ev30pz/V9ddfH3jOM888o+eee06bN2/Whx9+qOuuu055eXlqamoysjQAABAFDL3Xz7p16zR8+HCVl5cH9qWnpwe+9vv92rhxo5544gndc889kqSf/exnSk1N1ZtvvqkHHnjAyPIAAIDFGdqismPHDk2cOFH33XefUlJSdMstt+jll18OHD969Kjcbrdyc3MD+5xOp6ZMmaK9e/d2eE6fzyev1xv0AAAAscnQoPL73/9eZWVlGjNmjHbv3q2CggL96Ec/0iuvvCJJcrvdkqTU1NSg16WmpgaOXa60tFROpzPwGD58uJH/BAAAYCJDg0pra6vGjx+vNWvW6JZbbtHDDz+shx56SJs3bw77nMXFxfJ4PIHHiRMnerBiAABgJYaOUUlLS1NmZmbQvptuukn/+Z//KUlyuVySpMbGRqWlpQWe09jYqHHjxnV4TrvdLrvdbkzBAABEuf4JfXRs7Wyzy+gxhrao5OTk6MiRI0H7Pv30U40cOVLSpYG1LpdLVVVVgeNer1cffvihsrOzjSwNAICQnW++qFHLd2nU8l0633zR7HJ6BUNbVJYtW6Zbb71Va9as0f3336/9+/frpZde0ksvvSRJstlsWrp0qVavXq0xY8YoPT1dTz75pIYOHap7773XyNIAAEAUMDSoTJo0SW+88YaKi4u1atUqpaena+PGjZo/f37gOY899pjOnTunhx9+WGfOnNG0adNUWVmpxMREI0sDACDgfPNFZT61W5JUvypP/RMM/XiMunrMZPi//O6779bdd9/d6XGbzaZVq1Zp1apVRpcCAACiDPf6AQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQCgm1pa/YGv9x89HbQNYxBUAADohsq6BuWurw5s55fXaNq691RZ12BiVbGPoAIAQBcq6xpUsLVWjV5f0H63p0kFW2sJKwbqY3YBAABEwvnmi50eO9t0IfD1bz79Ujmjhyg+zibpUndPyY5D6qiTxy/JJmnFjvqg14Sif8KVH8WXdzFNH5Mc1rljgc3v90d1B5vX65XT6ZTH45HD4TC7HACARY1avsvsEjp0bO3soO3KugaV7DgU1HqT5kxUyZxMzcxKi3R5hunu5zddPwAAWARdTFei6wcAENWu1qXT3kdP3BG03dLq15zn/59OnfV18gop1WFXyZxM/X+v/q7L829eMF4TR12v/21ukST1S4jvVl1t9ZvRxRQN6PoBAEQ1q3bpWM3lXUxmo+sHAABEvehsBwIA4C/qV+WF9br9R08rv7ymy+dVLJqkyelJ2lPfqKd3fRLUVeRyJKr4/2TozszUwL62rpxQu1pCrae3IKgAAKJauGMvpo9JVpozUW5PU4fjQmySXM7EwNTge8b9lWZkpGjsinckXQoMHU0bjlQ9vQVdPwCAXik+zqaSOZmSLoWA9tq2S+ZkBoWC9l9PTk/q0cAQTj29AUEFANBrzcxKU9mC8Upx2IP2u5yJKlswPuLrllitHiug6wcA0KvNzEpTzughXXbp9NZ6zEaLCgCg1zOySyccVqvHTAQVAABgWQQVAABgWQQVAABgWQQVAABgWYYGlRUrVshmswU9MjIyAsebmppUWFiowYMHa8CAAZo3b54aGxuNLAkAAEQRw1tUvvWtb6mhoSHw+OCDDwLHli1bpp07d+r1119XdXW1Tp48qblz5xpdEgAAiBKGr6PSp08fuVyuK/Z7PB5t2bJF27Zt04wZMyRJ5eXluummm7Rv3z5NnTrV6NIAAIDFGd6i8tlnn2no0KG64YYbNH/+fB0/flySdODAAV24cEG5ubmB52ZkZGjEiBHau3ev0WUBAIAoYGiLypQpU1RRUaFvfvObamho0MqVKzV9+nTV1dXJ7XYrISFBgwYNCnpNamqq3G53p+f0+Xzy+b6+c6XX6zWqfAAAYDJDg8qsWbMCX998882aMmWKRo4cqZ///Ofq169fWOcsLS3VypUre6pEAABgYRGdnjxo0CDdeOON+vzzz+VyudTc3KwzZ84EPaexsbHDMS1tiouL5fF4Ao8TJ04YXDUAADBLRIPKV199pf/+7/9WWlqaJkyYoL59+6qqqipw/MiRIzp+/Liys7M7PYfdbpfD4Qh6AAAQCf0T+ujY2tk6tna2+idwX99IMPQq//3f/73mzJmjkSNH6uTJkyopKVF8fLy+//3vy+l0avHixSoqKlJSUpIcDoeWLFmi7OxsZvwAAABJBgeV//mf/9H3v/99/elPf1JycrKmTZumffv2KTk5WZK0YcMGxcXFad68efL5fMrLy9MLL7xgZEkAAFheW8sNJJvf7/ebXcS18Hq9cjqd8ng8dAMBwF+cb76ozKd2S5LqV+XRTQHL6e7nN/f6AQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQAYAY1NLqD3y9/+jpoG0gmhBUACDGVNY1KHd9dWA7v7xG09a9p8q6BhOrAsJDUAGAGFJZ16CCrbVq9PqC9rs9TSrYWktYQdTpY3YBAIBg55svhvW6lla/SnYcUkedPH5JNkkrdtQrZ/QQxcfZQjp3/wQ+LmAOvvMAwGIyn9ptyHn9ktzeJo1d8U7Irz22dnbPFwR0A10/AIAecb75okYt36VRy3eF3SoEXI4WFQCwmPpVeWG9bv/R08ovr+nyeRWLJmlyelJY7wFEGkEFACwm3PEg08ckK82ZKLenqcNxKjZJLmeipo9JDnmMCmAWun4AIEbEx9lUMidT0qVQ0l7bdsmcTEIKogpBBQBiyMysNJUtGK8Uhz1ov8uZqLIF4zUzK82kyoDw0PUDADFmZlaackYPCczuqVg0ie4eRC1aVAAgBrUPJZPTkwgpiFoEFQAAYFkEFQAAYFkEFQAAYFkEFQAAYFkRCypr166VzWbT0qVLA/uamppUWFiowYMHa8CAAZo3b54aGxsjVRIAALC4iASVmpoavfjii7r55puD9i9btkw7d+7U66+/rurqap08eVJz586NREkAACAKGB5UvvrqK82fP18vv/yyrr/++sB+j8ejLVu2aP369ZoxY4YmTJig8vJy/fa3v9W+ffuMLgsAAEQBw4NKYWGhZs+erdzc3KD9Bw4c0IULF4L2Z2RkaMSIEdq7d2+n5/P5fPJ6vUEPAAAQmwxdmXb79u2qra1VTc2Vd/N0u91KSEjQoEGDgvanpqbK7XZ3es7S0lKtXLmyp0sFAAAWZFiLyokTJ/Too4/q1VdfVWJiYo+dt7i4WB6PJ/A4ceJEj50bABC+ltav79m8/+jpoG0gXIYFlQMHDujUqVMaP368+vTpoz59+qi6ulrPPfec+vTpo9TUVDU3N+vMmTNBr2tsbJTL5er0vHa7XQ6HI+gBADBXZV2DctdXB7bzy2s0bd17qqxrMLEqxALDgsodd9yhjz/+WAcPHgw8Jk6cqPnz5we+7tu3r6qqqgKvOXLkiI4fP67s7GyjygIA9LDKugYVbK1Vo9cXtN/taVLB1tqQw8r55osatXyXRi3fpfPNF3uyVEQhw8aoDBw4UFlZWUH7rrvuOg0ePDiwf/HixSoqKlJSUpIcDoeWLFmi7OxsTZ061aiyAACdCCcUtLT6VbLjkDrq5PFLsklasaNeOaOHdPvGiIQTtGfoYNqubNiwQXFxcZo3b558Pp/y8vL0wgsvmFkSAPRamU/t7vFz+iW5vU0au+KdHj83egeb3++P6tFOXq9XTqdTHo+H8SoAcA1GLd9ldglXqF+Vp/4Jpv5NDYN09/Ob/30AgKRLoSBU+4+eVn75lUtQXK5i0SRNTk/q1jnPN1/UxNVVXT8RvQJBBQAgSWG1XEwfk6w0Z6LcnqYOx6nYJLmciZo+JrnbY1SA9rh7MgAgbPFxNpXMyZR0KZS017ZdMieTkIKwEVQAANdkZlaayhaMV4rDHrTf5UxU2YLxmpmVZlJliAV0/QAArtnMrDTljB4SmN1TsWgS3T3oEbSoAAB6RPtQMjk9ybCQwoJwvQtBBQAAWBZBBQAAWBZBBQAAWBZBBQAAWBZBBQAAWBZBBQAAWBZBBQBgKS2tXy/Gv//o6aBt9D4EFQCAZVTWNSh3fXVgO7+8RtPWvafKugYTq4KZCCoAAEuorGtQwdZaNXp9QfvdniYVbK0lrPRSLKEPAOhxoa4Y29LqV8mOQx3egdmvSzc4XLGjXjmjh8h3sSWk9wnnrtCwDpvf74/qzj+v1yun0ymPxyOHw2F2OQDQa51vvqjMp3abXcYVjq2dbXYJ6EB3P7/p+gEAAJZFexgAoMd99MQdIXW57D96WvnlNV0+r2LRJGX9lUMTV1eF9T6IPvzvAgB6XP+EPiEFiOljkpXmTJTb09ThOBWbJJczUdPHJAeNUQn1fRB96PoBAPSI/gl9dGztbB1bOzvk8BAfZ1PJnExJl0JJe23bJXMyFR93+VHEOoIKAMASZmalqWzBeKU47EH7Xc5ElS0Yr5lZaZJYEK63YdYPAMBSzjZd0NgV70i6NCZl+pjkQEtKZV2DSnYcClprJc2ZqJI5mYEgg+jArB8AQFRq370zOT0pKKSwIFzvwwgkAEBEdbVIW/vjbV+HsiBcOONYGJBrXXT9AAAiatTyXWaXcAUWhYs8un4AAEDUo60LABBR9avyrnr8fPPFKxZ0C2VBuMnpST1SJ6yBoAIAiKhQxoO0LegWyoJwrLUSW+j6AQBYHgvC9V4EFQBAVOjugnCILYYGlbKyMt18881yOBxyOBzKzs7W22+/HTje1NSkwsJCDR48WAMGDNC8efPU2NhoZEkAgCg2MytN7xZ9J7BdsWiSPnh8BiElhhkaVIYNG6a1a9fqwIED+uijjzRjxgzdc889OnTokCRp2bJl2rlzp15//XVVV1fr5MmTmjt3rpElAQCiXGcLwiE2GTqYds6cOUHbTz/9tMrKyrRv3z4NGzZMW7Zs0bZt2zRjxgxJUnl5uW666Sbt27dPU6dONbI0AAAQBSI2RqWlpUXbt2/XuXPnlJ2drQMHDujChQvKzc0NPCcjI0MjRozQ3r17Oz2Pz+eT1+sNegAAgNhkeFD5+OOPNWDAANntdv3whz/UG2+8oczMTLndbiUkJGjQoEFBz09NTZXb7e70fKWlpXI6nYHH8OHDDf4XAAAAsxgeVL75zW/q4MGD+vDDD1VQUKCFCxeqvr4+7PMVFxfL4/EEHidOnOjBagEgupxvvqhRy3dp1PJdXd5DB4hGhi/4lpCQoNGjR0uSJkyYoJqaGv3zP/+zvve976m5uVlnzpwJalVpbGyUy+Xq9Hx2u112u73T4wAAIHZEfB2V1tZW+Xw+TZgwQX379lVVVVXg2JEjR3T8+HFlZ2dHuiwAiAhaQLrWP6GPjq2drWNrZ3NXYxjbolJcXKxZs2ZpxIgROnv2rLZt26b3339fu3fvltPp1OLFi1VUVKSkpCQ5HA4tWbJE2dnZzPgBAACSDA4qp06d0g9+8AM1NDTI6XTq5ptv1u7du3XnnXdKkjZs2KC4uDjNmzdPPp9PeXl5euGFF4wsCQAARBFDg8qWLVuuejwxMVGbNm3Spk2bjCwDAABEKe71AwAALIugAgAALIugAgAALIugAgAxjinRiGZMUAcARJW2dVbQO9CiAgAALIugAgAALIugYjD6hgEACB9BBQAAWBZBBQAAWBZBBQCiWEurP/D1/qOng7aBWEBQAYAoVVnXoNz11YHt/PIaTVv3nirrGkysCuhZBBUAiKCeagGprGtQwdZaNXp9QfvdniYVbK0lrCBmsOAbAERIZV2DSnYcCmznl9co1WHXj//PTbozM7Xb52lp9atkxyF1FHH8kmySVuyoV87oIYqPswXNOOzO7MP+CXw0wDr4bgSACGhrAbk8XDR6fXp0+8EefS+/JLe3SWNXvHPFsYmrq7p8Pau+wkoIKgDQzvnmi8p8arckqX5VXlDrQrhrIV2tBQTA1RFUAKCb2gJMtKhYNEmT05N0vvlioCXloyfuoGsHUYXvVgCIMTZJLmeipo9JVnycLehY/4Q+BBVEFb5bAaCb6lflhfW6/UdPK7+8psvntbWAdMee+kYt3X7wiu6ktlhSMifzipACRCOCCgB0U7gtEdPHJCvNmSi3p6nDcSpXawHpzD3j/kr2PnEq2XEoaIqyy5mokjmZmpmVFth3+ZToUN4HMBvrqACAweLjbCqZkynp6xaPNtfSAjIzK03vFn0nsF2xaJI+eHxGUEhhUThEO4IKAETAzKw0lS0YrxSHPWi/y5mosgXjg8JFKNqHm8npSUHbLAqHWEDXDwBEyMysNOWMHhJY36Ri0aQuu2G6mhLd2WJuoS4KFwoG4yKS+G4zGH3DANq7WgtIR0KZEt2dxdzaXG1RuK6wIBwiia4fA9E3DADAtaFFxSCdLZfd1jd8LX3SAHqPrqZEd7aYmxFTogEzEFS6EM6S2Ub2DUv0DwO9SSg/7+0XczNiSjRgBj7xumDEktnX0jcs0T8MoGttU6ILttbKJgWFFRaFQzRhjAoAxCijpkQDkUSLShfCWTKbvmEAVhHOlGjASgxtUSktLdWkSZM0cOBApaSk6N5779WRI0eCntPU1KTCwkINHjxYAwYM0Lx589TY2GhkWSFp6/MN5dHWN9zZrwGbpLS/9A2Hc34ACEWoU6IBKzE0qFRXV6uwsFD79u3Tnj17dOHCBd111106d+5c4DnLli3Tzp079frrr6u6ulonT57U3LlzjSzLcEYtlw3AeJevfdR+G0DkGfrneWVlZdB2RUWFUlJSdODAAd12223yeDzasmWLtm3bphkzZkiSysvLddNNN2nfvn2aOnWqkeUZqq1vuDs3DGvvfPPFwADe+lV5tKAAEVRZ16CSHYcC2/nlNUrr4mcWgLEiOpjW4/FIkpKSLo3LOHDggC5cuKDc3NzAczIyMjRixAjt3bs3kqUZojs3DANgDdwXB7CmiP253traqqVLlyonJ0dZWVmSJLfbrYSEBA0aNCjouampqXK73R2ex+fzyef7+heJ1+s1rOaeQN8wEDnhrHsksfYRYGUR++kpLCxUXV2dPvjgg2s6T2lpqVauXNlDVQGIJUaseySx9hFgpoh0/TzyyCP65S9/qV//+tcaNmxYYL/L5VJzc7POnDkT9PzGxka5XK4Oz1VcXCyPxxN4nDhxwsjSAcDS+if00bG1s3Vs7WxabhCTDP2u9vv9WrJkid544w29//77Sk9PDzo+YcIE9e3bV1VVVZo3b54k6ciRIzp+/Liys7M7PKfdbpfdbu/wGIDeLZx1jyTWPgKszNCgUlhYqG3btumtt97SwIEDA+NOnE6n+vXrJ6fTqcWLF6uoqEhJSUlyOBxasmSJsrOzo3rGDwBzhNuiEMn74rS1gADoHkO7fsrKyuTxeHT77bcrLS0t8HjttdcCz9mwYYPuvvtuzZs3T7fddptcLpd+8YtfGFkWAARh7SPAugzv+ulKYmKiNm3apE2bNhlZCgBcVbhrHwEwFiOvAOAvuC8OYD3cPRkA2mHtI8BaCCoAAMCyCCoAAMCyGKMCADGOKdGIZrSoAAAAyyKoAAAAyyKoWExL69drz+w/ejpoGwCA3oYxKgYLpW+4sq5BJTsOBbbzy2uUxmJTQNjON18M3FG5flUeN+0DohAtKhZRWdeggq21QStiSpLb06SCrbWqrGswqTIgdp1vvqhRy3dp1PJdOt980exyAHSAPy96ULi/6Fpa/SrZcajDm6H5deleIyt21Ctn9JCwFp/ir0gAQLTiE6wHtTUx9zS/JLe3KbCsd6iYlggAiFZ0/QAAAMuiRaUH1a/KC+t1+4+eVn55TZfPq1g0SZPTk8J6DwAAohFBpQeFOxZk+phkpTkT5fY0dThOxaZLt5rnLq4AgN6Grh8LiI+zqWROpqRLoaS9tu2SOZmEFABAr0NQsYiZWWkqWzBeKQ570H6XM1FlC8azjgrQDT0x3bht7aNja2czYw6wAH4KLWRmVppyRg8JzO6pWDSJ7h4AQK9Gi4rFtA8lk9OTCCkAgF6NoAIAACyLoAIAACyLoAIAACyLoAIAACyLoAIAACyLoAIgZrW0fr3W8/6jp4O2AUQHggoiqicW5AK6Y099o3LXVwe288trNG3de6qsawjsI8gA1seCbwBi0tLtB6+4d5bb06SCrbUqWzBeklSy41DgWH55jdKciSqZk8lK0ICFEFQAWFI4LW7tX9NR24hfl+6ftfw/P5bnfy90GmQ2PjBOd2amdvgeLKsPRBY/cQCuyfnmi8p8arckqX5VXo99kLeds6f5JZ353wudHpOkR7cf7PT19avyrnp+ggzQs/iJgiUY9WEH9LSuAtSxtbMjVAnQO/BpAMCSumq5MKrFBYC1GBpUfvOb3+if/umfdODAATU0NOiNN97QvffeGzju9/tVUlKil19+WWfOnFFOTo7Kyso0ZswYI8sCEAW6alXrKMicb76oiaureuT9KxZN0uT0pB45F4DwGTo9+dy5c/r2t7+tTZs2dXj8mWee0XPPPafNmzfrww8/1HXXXae8vDw1NTUZWRaAGNA/oc8VD3uf+KDnXH7v8bbtQf37XnGs/XPSnImaPia5w/fo6gGgZxn6UzVr1izNmjWrw2N+v18bN27UE088oXvuuUeS9LOf/Uypqal688039cADDxhZmmX1T+hDHzcQhsq6hqDpxpJks0n+dlN7XH+ZfixJBVtrZVPw7KC28FIyJ1PxcZ1FGQCRZFr8P3r0qNxut3JzcwP7nE6npkyZor1793YaVHw+n3w+X2Db6/UaXisAa6usa1DB1torphu3X7+tYtEkTR+THAggZQvGq2THITV6v/594mIdFcByTAsqbrdbkpSaGrxWQWpqauBYR0pLS7Vy5UpDa4t2zKBBNAp3peKWVr9KdhzqcN2U9m5KGyjfxZbA9m03JmvHIzmasuY9SdLmBeOVM3qI4uNsQbXw8wOYK+p+AouLi1VUVBTY9nq9Gj58uIkVAegJRs/iaQsknfnh1toO99MVC5jLtHv9uFwuSVJjY2PQ/sbGxsCxjtjtdjkcjqAHAACITaa1qKSnp8vlcqmqqkrjxo2TdKl15MMPP1RBQYFZZQEwSVfrpnRm/9HTyi+v6fJ5mxeM1203Jgftaz+d+aMn7qCbB7AgQ38qv/rqK33++eeB7aNHj+rgwYNKSkrSiBEjtHTpUq1evVpjxoxRenq6nnzySQ0dOjRorRUAnYul8Ujh1j59TLLSnIlye5quOk4lZ/SQq74H04sBazL0p/Kjjz7Sd7/73cB229iShQsXqqKiQo899pjOnTunhx9+WGfOnNG0adNUWVmpxMREI8sCEEPi42wqmZPZ6XRjf7vnAYg+ho5Ruf322+X3+694VFRUSJJsNptWrVolt9utpqYmvfvuu7rxxhuNLAkma2k3X3T/0dNB20C4ZmalqWzBeKU47EH7Ux380QNEO9MG06JnnG++qFHLd2nU8l1hT++MlMq6BuWurw5s55fXaNq691RZ12BiVYgVM7PS9G7RdwLbFYsmaU/RbSZWBKAnEFQQEW0LcrVfXEuS3J4mFWyt1Z76xk5eCXRf++6dyelJdPcAMYCRY+g2Ixbk8uvSOII1v/ok7PeJ1ADIWBq42pMu785rv/orAFwrftOi24xakMsvBbW0hHr3WysvyBXr4eby++vkl9cojWXoAfQgun4AhKWr7jzGHgHoCbH15x0MZfSCXO21dRxsfGCc7sxMvepzET4ju/NW7KgP3DsnVLHW8gQgfPw2QLd19uHR1YfdhJHXK9Vh1ymvr8sbx7Vp+7Ar/dVhzchIueqHXWfvz4dd14zsznN7mzR2xTthvd7K3XkAIovf5LhmfNgBAIxCUAF6MaO78yoWTdLk9KSw3qMn9E/oQ2AFohxBJQZFerpoKB92Z5suaMqa9yRJy3LHaMO7n3X5GrM/7GKZUffXsUlyORMtP1WZIANYH7N+YowZq7+23cytO4+BiX0Dr3twWrrSnInq7GPMJintLx92obwHN5frWR2tftx2fx1JV/z/tW2XzMm0dEgBEB0IKjFkT31jVE0X5cMuunV2fx2XM1FlC8azjgqAHsGfnRYT6nTR9s9/etcnUTddtO3DrmTHoaCA5WLRsKgwMytNOaOHBAY8VyyaZPnuHgDRhaBiMdcyg+bUWV+nx6w8g4YPu/BZYfl67q8DwEh0/cA0mU/tDox74MMudKGMR4qmu2wDQHu0qFhMqNNFzzdfDOneOMygiQ1ty9df3tXXNh6pt44RYRYPEHsIKhYT6liQ9k3/1/fvqzPnL0T1dNHeJJLL17d/r67elxlTAKyE30hR7PI71/75/IUOn8cMGmsya0XfrlrgaJEAYCWMUYlSnd25tiNMF+0Zlw9cbb8NADAGLSomMqLp/3KbF4wPdAF09/1o+r/S5a1X+eU1SuvGFOqrzcqJ5PL17ccyffTEHfwfA4ga/LYykVFN/+39cGttyK+h6T9YuANXuwo3Zi1fz8q9AKIJXT/oNc43Xwz5cbbpwlUHrkqXBq6ebboQ9Lq3Dn5x1VWC3zr4RdgtaqzoC6A34c8qExnd9C/RzN+eES1YoS6k1xZuHt1+UFJ43wP9E/qwoi+AXoNPMBMZ2fTftp9mfmsLJzy1dc2xoi+A3oBPsCjU1vRfsLU2KJRIV3YF4GvhtF5EYhzRtWBFXwCxjqASpa7W9L98Vkaga8Fq2lYOPd98MeIhIJyWpY9X3KXc9dU65fV12nqV6kjUnqLbAiEhnFk5AICOEVSiWGdN/76LLSZXFjsGJvbVyv/7rau2Xq34v5kamNg3sP9aZ+VYjRVufAig92LWT5Sj6d94ba1XKQ570P7OFtKLpVk53bnxYVsr2bG1sxkPBaDHEVSAbpiZlaZ3i74T2K5YNEkfPD6j09k1oYYbK+ps9eO2KdYd3aUZAHoaf/7AEqLhrrehtl6FMiun/Zid+lV5PdYy0bZWSyg3JZTCu/FhKGh5AdBd/LaAKXrLuAezu+Y6GrDc1U0JuyPU9WMuZ/VQCsA6LNH1s2nTJo0aNUqJiYmaMmWK9u/fb3ZJMFB3xj0AACBZoEXltddeU1FRkTZv3qwpU6Zo48aNysvL05EjR5SSkmJ2eehh4d43B+H55wfGaen2g1dc77Z2nY0PjNOdmalXvI4p1gCswvSgsn79ej300ENatGiRJGnz5s3atWuX/u3f/k3Lly83uTp0xIi7PjPuoXNXu95XG3vS0urXml99ctXrXfqrw5qRkXLF9Z4w8nqlOuxXXT8mmqZYA4hepv5mb25u1oEDB1RcXBzYFxcXp9zcXO3du7fD1/h8Pvl8X89C8Hq9hteJYEYt1Ma4h45193qHOvbkWq63X9EzxRpAdDM1qPzxj39US0uLUlODm55TU1N1+PDhDl9TWlqqlStXRqK8qBUNM2gQ/eiiAxAJUddWXlxcrKKiosC21+vV8OHDTayo9zH6rs+Mewh2tet9vvlioCXl8jtl99T1Ptt0QVPWvCdJ2rxgfNhdcwAQDlODypAhQxQfH6/Gxsag/Y2NjXK5XB2+xm63y263d3gMkWHkXZ8Z93Cl7l7vy++UbcT1vu3G5JgdCwTAmkydnpyQkKAJEyaoqurrvvXW1lZVVVUpOzvbxMpghFhaWj4acL0BxALT11EpKirSyy+/rFdeeUWffPKJCgoKdO7cucAsIMSWWFha3kq6us8O1xtAtDO9Dfd73/uevvzySz311FNyu90aN26cKisrrxhgi9gRytLykWDU8vVWYbXrDQChsMRv5EceeUSPPPKI2WUggsxeWr634XoDiFaWCCoIH1ORAQCxzPQxKgAAAJ2hRQXoJlqvACDyCCqAgbobblpav17pZP/R0wx2BYC/oOsHMFllXYNy11cHtvPLazRt3XuqrGswsSoAsAaCCmCiyroGFWytVaPXF7Tf7WlSwdZawgqAXo+uH+AanW++GNbrWlr9KtlxqMPl7f26tHrsih31Yd9bJ9bWgwHQO/GbDLhGbYvF9TS/JLe3KbBQW6gY+AsgFtD1AwAALIsWFeAa1a/KC+t1+4+eVn55TZfPq1g0SZPTk8J6DwCIdgQV4BqFOxZk+phkpTkT5fY0dThOxaZLNw80e6oy68cAMBNdP4BJ4uNsKpmTKelSKGmvbbtkTibrqQDo1QgqgIlmZqWpbMF4pTjsQftdzkSVLRivmVlpJlUGANZA1w9gsplZacoZPSQwu6di0STTu3sAwCoIKjCFlcY9WGH5+vbvNzk9qcff30rXGwBCQdcPejWWrwcAayOooNdi+XoAsD66fhDVWL4eAGIbv00R1Vi+HgBiG10/AADAsmhRQVRj+XoAiG0EFUS1WF++HgB6O7p+0CuxfD0ARAeCCnotlq8HAOuj6we9GsvXA4C10aKCXs/o5esBAOEjqAAAAMsiqAAAAMsiqAAAAMtiMC1gAf0T+rDsPgB0gBYVAABgWYYFlaefflq33nqr+vfvr0GDBnX4nOPHj2v27Nnq37+/UlJS9A//8A+6eDG8u+ECAIDYY1jXT3Nzs+677z5lZ2dry5YtVxxvaWnR7Nmz5XK59Nvf/lYNDQ36wQ9+oL59+2rNmjVGlQUAAKKIYS0qK1eu1LJlyzR27NgOj7/zzjuqr6/X1q1bNW7cOM2aNUs/+clPtGnTJjU3NxtVFgAAiCKmjVHZu3evxo4dq9TU1MC+vLw8eb1eHTp0qNPX+Xw+eb3eoAcAAIhNpgUVt9sdFFIkBbbdbnenrystLZXT6Qw8hg8fbmidAADAPCEFleXLl8tms131cfjwYaNqlSQVFxfL4/EEHidOnDD0/QAAgHlCGkz7d3/3d8rPz7/qc2644YZuncvlcmn//v1B+xobGwPHOmO322W32zs9DgAAYkdIQSU5OVnJyck98sbZ2dl6+umnderUKaWkpEiS9uzZI4fDoczMzB55DwAAEN0Mm558/PhxnT59WsePH1dLS4sOHjwoSRo9erQGDBigu+66S5mZmfrbv/1bPfPMM3K73XriiSdUWFhIiwkiilVhAcC6bH6/32/EifPz8/XKK69csf/Xv/61br/9dknSH/7wBxUUFOj999/Xddddp4ULF2rt2rXq06f7+cnr9crpdMrj8cjhcPRU+QAAwEDd/fw2LKhECkEFAIDo093Pb+71AwAALIugAgAALIugAgAALIugAgAALIugAgAALIugAgAALIugAgAALIugAgAALIugAgAALIugAgAALIugAgAALIugAgAALIugAgAALKuP2QVcq7abP3u9XpMrAQAA3dX2ud32Od6ZqA8qZ8+elSQNHz7c5EoAAECozp49K6fT2elxm7+rKGNxra2tOnnypAYOHCibzWZqLV6vV8OHD9eJEyfkcDhMrcXquFah4XqFhuvVfVyr0HC9uq+ra+X3+3X27FkNHTpUcXGdj0SJ+haVuLg4DRs2zOwygjgcDr6Bu4lrFRquV2i4Xt3HtQoN16v7rnatrtaS0obBtAAAwLIIKgAAwLIIKj3IbrerpKREdrvd7FIsj2sVGq5XaLhe3ce1Cg3Xq/t66lpF/WBaAAAQu2hRAQAAlkVQAQAAlkVQAQAAlkVQAQAAlkVQMZjP59O4ceNks9l08OBBs8uxpGPHjmnx4sVKT09Xv3799I1vfEMlJSVqbm42uzRL2LRpk0aNGqXExERNmTJF+/fvN7skSyotLdWkSZM0cOBApaSk6N5779WRI0fMLisqrF27VjabTUuXLjW7FMv64osvtGDBAg0ePFj9+vXT2LFj9dFHH5ldliW1tLToySefDPqd/pOf/KTLe/p0hqBisMcee0xDhw41uwxLO3z4sFpbW/Xiiy/q0KFD2rBhgzZv3qwf//jHZpdmutdee01FRUUqKSlRbW2tvv3tbysvL0+nTp0yuzTLqa6uVmFhofbt26c9e/bowoULuuuuu3Tu3DmzS7O0mpoavfjii7r55pvNLsWy/vznPysnJ0d9+/bV22+/rfr6ev30pz/V9ddfb3ZplrRu3TqVlZXpX/7lX/TJJ59o3bp1euaZZ/T888+Hd0I/DPOrX/3Kn5GR4T906JBfkv93v/ud2SVFjWeeecafnp5udhmmmzx5sr+wsDCw3dLS4h86dKi/tLTUxKqiw6lTp/yS/NXV1WaXYllnz571jxkzxr9nzx7/d77zHf+jjz5qdkmW9Pjjj/unTZtmdhlRY/bs2f4HH3wwaN/cuXP98+fPD+t8tKgYpLGxUQ899JD+/d//Xf379ze7nKjj8XiUlJRkdhmmam5u1oEDB5SbmxvYFxcXp9zcXO3du9fEyqKDx+ORpF7/fXQ1hYWFmj17dtD3GK60Y8cOTZw4Uffdd59SUlJ0yy236OWXXza7LMu69dZbVVVVpU8//VSS9F//9V/64IMPNGvWrLDOF/U3JbQiv9+v/Px8/fCHP9TEiRN17Ngxs0uKKp9//rmef/55Pfvss2aXYqo//vGPamlpUWpqatD+1NRUHT582KSqokNra6uWLl2qnJwcZWVlmV2OJW3fvl21tbWqqakxuxTL+/3vf6+ysjIVFRXpxz/+sWpqavSjH/1ICQkJWrhwodnlWc7y5cvl9XqVkZGh+Ph4tbS06Omnn9b8+fPDOh8tKiFYvny5bDbbVR+HDx/W888/r7Nnz6q4uNjskk3V3evV3hdffKGZM2fqvvvu00MPPWRS5Yh2hYWFqqur0/bt280uxZJOnDihRx99VK+++qoSExPNLsfyWltbNX78eK1Zs0a33HKLHn74YT300EPavHmz2aVZ0s9//nO9+uqr2rZtm2pra/XKK6/o2Wef1SuvvBLW+VhCPwRffvml/vSnP131OTfccIPuv/9+7dy5UzabLbC/paVF8fHxmj9/ftj/WdGmu9crISFBknTy5Endfvvtmjp1qioqKhQX17tzdHNzs/r376//+I//0L333hvYv3DhQp05c0ZvvfWWecVZ2COPPKK33npLv/nNb5Senm52OZb05ptv6q//+q8VHx8f2NfS0iKbzaa4uDj5fL6gY73dyJEjdeedd+pf//VfA/vKysq0evVqffHFFyZWZk3Dhw/X8uXLVVhYGNi3evVqbd26NazWYLp+QpCcnKzk5OQun/fcc89p9erVge2TJ08qLy9Pr732mqZMmWJkiZbS3eslXWpJ+e53v6sJEyaovLy814cUSUpISNCECRNUVVUVCCqtra2qqqrSI488Ym5xFuT3+7VkyRK98cYbev/99wkpV3HHHXfo448/Dtq3aNEiZWRk6PHHHyekXCYnJ+eKqe6ffvqpRo4caVJF1nb+/PkrfofHx8ertbU1rPMRVAwwYsSIoO0BAwZIkr7xjW9o2LBhZpRkaV988YVuv/12jRw5Us8++6y+/PLLwDGXy2ViZeYrKirSwoULNXHiRE2ePFkbN27UuXPntGjRIrNLs5zCwkJt27ZNb731lgYOHCi32y1Jcjqd6tevn8nVWcvAgQOvGLtz3XXXafDgwYzp6cCyZct06623as2aNbr//vu1f/9+vfTSS3rppZfMLs2S5syZo6efflojRozQt771Lf3ud7/T+vXr9eCDD4Z3wmudhoSuHT16lOnJV1FeXu6X1OEDfv/zzz/vHzFihD8hIcE/efJk/759+8wuyZI6+x4qLy83u7SowPTkq9u5c6c/KyvLb7fb/RkZGf6XXnrJ7JIsy+v1+h999FH/iBEj/ImJif4bbrjB/4//+I9+n88X1vkYowIAACyLgQAAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCy/n+djjKcH6FahwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.errorbar(data_x, data_y, sigma_y, sigma_x, fmt=\"o\");" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:25:44.032210Z", "start_time": "2020-02-21T10:25:43.648365Z" } }, "outputs": [ { "data": { "text/plain": [ "DeviceArray(876.49545695, dtype=float64)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# define the cost function\n", "@jax.jit\n", "def cost(par):\n", " result = 0.0\n", " for xi, yi in zip(data_x, data_y):\n", " y_var = sigma_y ** 2 + (f_prime(xi, par) * sigma_x) ** 2\n", " result += (yi - f(xi, par)) ** 2 / y_var\n", " return result\n", "\n", "cost.errordef = Minuit.LEAST_SQUARES\n", "\n", "# test the jit-ed function\n", "cost(np.zeros(3))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:25:44.059729Z", "start_time": "2020-02-21T10:25:44.034029Z" } }, "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 = 23.14 Nfcn = 91
EDM = 3.12e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0 1.25 0.15
1 x1 1.5 0.5
2 x2 1.6 1.5
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x0 x1 x2
x0 0.0223 -0.0388 (-0.530) -0.15 (-0.657)
x1 -0.0388 (-0.530) 0.24 0.172 (0.230)
x2 -0.15 (-0.657) 0.172 (0.230) 2.32
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 23.14 │ Nfcn = 91 │\n", "│ EDM = 3.12e-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.25 │ 0.15 │ │ │ │ │ │\n", "│ 1 │ x1 │ 1.5 │ 0.5 │ │ │ │ │ │\n", "│ 2 │ x2 │ 1.6 │ 1.5 │ │ │ │ │ │\n", "└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌────┬─────────────────────────┐\n", "│ │ x0 x1 x2 │\n", "├────┼─────────────────────────┤\n", "│ x0 │ 0.0223 -0.0388 -0.15 │\n", "│ x1 │ -0.0388 0.24 0.172 │\n", "│ x2 │ -0.15 0.172 2.32 │\n", "└────┴─────────────────────────┘" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = Minuit(cost, np.zeros(3))\n", "m.migrad()" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:25:44.566228Z", "start_time": "2020-02-21T10:25:44.065443Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAG3CAYAAAAU+jfPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYgElEQVR4nO3deZyNdf/H8deZfQYzDGZRmLGUpOz7GopiuivRphvtfijcd3e0jShUWiVLdVOE9pCaFkJ1WyZS2Uu2MIMwgzHruX5/nOZwzHrOnDnXOTPv5+NxHq79fM41Y67P+a4WwzAMRERERDzEz+wAREREpHJR8iEiIiIepeRDREREPErJh4iIiHiUkg8RERHxKCUfIiIi4lFKPkRERMSjlHyIiIiIRyn5EBEREY9S8iEiIiIepeRDBMjKyuKuu+6iXr16hIeH06FDB9auXWt2WCIiFZKSDxEgNzeXuLg4vv/+e06ePMno0aNJSEjg9OnTZocmIlLhWDSxnEjh6tSpw7Jly2jdurXZoYiIVCgq+RApxG+//cbx48dp1KhRgX1Wq5XatWvz3HPPmRBZ+UhOTmbkyJFcfvnlVKlShXr16jFo0CB27dpV4NitW7cycOBAGjRoQFhYGLVq1aJbt24sW7asVO91+vRpEhMT6du3L5GRkVgsFubNm1fiec888wwWi4VmzZoVecz5PxtnPtPQoUOxWCxFvg4ePFiqz+aqrKwsHnnkEerUqUNoaCjt27fn66+/LtW5v/32G7feeisXX3wxYWFhNGnShIkTJ5KRkVHo8Zs2beL6668nMjKSsLAwmjVrxquvvurOjyNSogCzAxDxNmfPnmXw4MGMHz+eiIiIAvs3bNjAsWPH6NevnwnRlY9nn32WH374gYEDB3LllVeSkpLCa6+9RqtWrVi3bp3DA3/fvn2cOnWKIUOGUKdOHTIyMvjoo4+4/vrrmT17Nvfdd1+x73Xs2DEmTpxIvXr1aN68OatWrSoxvj///JPJkydTpUqVYo87/2eTmJhY6s90//3307t3b4drGYbBAw88QFxcHBdddFGJMZbF0KFD+fDDDxk9ejSNGzdm3rx5XHfddXz77bd06dKlyPMOHDhAu3btiIiIYOTIkURGRrJ27VoSExPZuHEjS5YscTj+q6++IiEhgZYtW/LEE09QtWpVdu/ezZ9//lmun0+kAENE7LKzs41+/foZt99+u2G1Wgs95oknnjDq16/v2cDK2Q8//GBkZWU5bNu1a5cRHBxs3HHHHSWen5ubazRv3ty49NJLSzw2MzPTOHz4sGEYhpGcnGwAxty5c4s955ZbbjF69uxpdO/e3bj88suLPO78n01ZP9N3331nAMYzzzxT4rFlsX79egMwnn/+efu2s2fPGg0bNjQ6duxY7LnPPPOMARhbtmxx2P7Pf/7TAIzjx4/bt6WlpRnR0dHGjTfeaOTl5bn3Q4g4SdUuUqG9+eabhISE0LlzZ/bt22ffbhgGV111FbVq1eLIkSOArcj+zjvvxGKx8Pbbb2OxWAq95vLlyx1KPa666iq6devGpk2buPbaa6lWrRoXXXQRr7zySvl+ODfq1KkTQUFBDtsaN27M5Zdfzvbt20s839/fn7p163Ly5MkSjw0ODiYmJqbUsa1Zs4YPP/yQl19+ucRjz//ZlPUzLVy4EIvFwu23317qWF3x4Ycf4u/v71BiFBISwt13383atWs5cOBAkeemp6cDEB0d7bA9NjYWPz8/h8+/cOFCUlNTeeaZZ/Dz8+PMmTNYrVY3fxqR0lHyIRVa27Ztefjhh1m3bh3Tpk2zb58xYwarVq1i+vTpREVFAbai98OHD/PBBx8QEFB4jWRKSgo//fQT1113nX3br7/+ysmTJ0lISKB169ZMmzaN2NhYxowZw6+//lq+HxDIycnh2LFjpXo587AxDIPU1FRq1apV6P4zZ85w7Ngxdu/ezUsvvcQXX3xBr1693PWxAMjLy2PUqFHcc889XHHFFcUeW9jP5kIlfaZ8OTk5vP/++3Tq1Im4uLgij3HHff/pp5+45JJLCA8Pd9jerl07ADZv3lzkuT169ADg7rvvZvPmzRw4cID33nuPmTNn8uCDDzpUU33zzTeEh4dz8OBBLr30UqpWrUp4eDjDhw8nMzOz2Psh4nYml7yIeMQ111xjL8LevXu3UaVKFeOGG26w79+7d68BGCEhIUaVKlXsrzVr1jhc56233jJCQ0ONjIwMwzAM49ChQwZg1K5d2zhw4ID9uG3bthmA8fbbb5f7Z/v2228NoFSvPXv2lPq68+fPNwDjrbfeKnT//fffb7+un5+fcfPNNzsU85dGSdUur732mhEREWEcOXLEMAyj2GqXC382rnymfMuWLTMA4/XXXy/yGHfd98svv9zo2bNnge1bt241AGPWrFnFxjpp0iQjNDTU4f0ee+yxAsddeeWVRlhYmBEWFmaMGjXK+Oijj4xRo0YZgHHrrbcW+x4i7qYGp1IptGjRgpkzZ2K1WrnrrrsIDg5m5syZ9v3169fHKEWv888//5yrrrqK0NBQAHvJRmJiIhdffLH9uMDAQIACxf7n69GjB/fccw+DBw8u9j23b9/OoEGD2Lt3L/PmzWPAgAEO+5s3b17qnhGlre7YsWMHI0aMoGPHjgwZMqTQY0aPHs3NN9/MoUOHeP/998nLyyM7O7tU1y+Nv/76iyeffJInnniC2rVrl3j8hT+bC5XmM+VbuHAhgYGBDBo0qMhj3HXfz549S3BwcIHtISEh9v3FiYuLo1u3bgwYMICaNWuyfPlyJk+eTExMDCNHjrQfd/r0aTIyMnjggQfsvVtuuukmsrOzmT17NhMnTqRx48al+jwiZWZ29iPiCe+8844BGGPHjjUAY/78+U5fIzs72wgPDzdmzJhh3zZt2jQDMP7880+HY/O/OW/atKnI63Xv3r1UcQwbNswYP3680/G66vDhw0aDBg2MunXrGgcPHiz1eVdffbXRtm3bIhvqFqa4ko8HHnjAaNSokUOj0aJKPgr72ZzPmc906tQpIywszOjfv3+pP0dZlKXkY9GiRUZoaKhDqZthGMbQoUONsLAw49ixYw7vAxirV692OHb16tUeK6UTyaeSD6kU8rtVvvjii/Tv37/E0obCfP/996Snpzu0Kfjll1+IiYkp0BXz559/JiAggKZNm5YtcGD//v307NmzyP3Z2dkcP368VNeqXbs2/v7+Re5PS0vj2muv5eTJk3z33XfUqVOn1HHefPPN3H///ezatYtLL7201OcV5rfffmPOnDm8/PLLHDp0yL49MzOTnJwc9u7dS3h4OJGRkUDhPxtXP9Onn35KRkYGd9xxR7HHueu+x8bGFjqOyOHDhwGKjff111+nZcuWDqVuANdffz3z5s3jp59+snchrlOnDlu3bi3QODW/zdOJEydK9VlE3EENTqVSyH8YVq9endmzZ7t0jeXLl9O0aVOHBoi//vorzZs3L3DsL7/8wiWXXOJQnJ6cnMyVV15JeHg4DzzwgEMjxK1bt9K1a1eqV69O69at+eGHHwC49tpr+fbbb7nnnnuoWrUqf/31V4H3+t///kdsbGypXsX1nMjMzCQhIYFdu3bx2WefOZ045VcPpKWlOXVeYQ4ePIjVauXBBx8kPj7e/lq/fj27du0iPj6eiRMn2o8v7Gfj6md69913qVq1Ktdff32xx7nrvrdo0YJdu3bZe67kW79+vX1/UVJTU8nLyyuwPScnB7BNG5Avf6TeCxOd/OSuNFVbIu6ikg+pFN544w3A9o3QmW/z5/v888/p37+/fT0vL4/t27dz9dVXFzj2559/pmXLlvb17OxsbrrpJh599FHuueceZs2axZtvvsl9991HdnY2CQkJjB49mpUrV/Lxxx+TkJDA7t27+eKLL0psG+KOtgd5eXnccsstrF27liVLltCxY8cir3HkyBH7t+V8OTk5vPPOO4SGhjo84DMyMti/fz+1atUqsYfJ+Zo1a8Ynn3xSYPvjjz/OqVOneOWVV2jYsKF9+4U/G2c/U76jR4/yzTffcNtttxEWFlbsse5q83HzzTczbdo05syZw7///W/ANuLp3Llzad++PXXr1gUKv5eXXHIJX331Fbt27eKSSy6xX3PRokX4+flx5ZVX2rcNGjSIqVOn8tZbbzmUpL355psEBATYe86IeITZ9T4i5e333383wsLCDMBo166dS9f4448/DMBYtWqVfdv27dsNwHj33Xcdjs3IyDD8/f2NyZMn27etWrXKiIuLs69brVbj4osvNubPn2+sWbOmwKBlHTp0MBYuXGgYRunbhpTFQw89ZABGQkKCMX/+/AKv891www1Gz549jQkTJhhvvPGGMWnSJKNJkyYGYLzwwgsOx+b3CElMTHTYPn36dGPSpEnG8OHDDcC46aabjEmTJhmTJk0yTp48WWSchbX5KOxn4+xnOj8uwEhKSirplrnVwIEDjYCAAOPhhx82Zs+ebXTq1MkICAhwaJ9R2L1cvXq14e/vb0RFRRkTJ040ZsyYYVx77bUGYNxzzz0F3ueuu+4yAGPQoEHGjBkzjIEDBxqAR9sUiRiGYSj5kArNarUa3bt3N2rUqGEMGzbMqFq1qlMNIvPld/nMycmxb3v//fcLHV1yw4YNBmB89tln9m2LFi0yunTp4nBchw4djPnz5xuLFy8usO+WW24xpk2bZhiGZ5KP7t27F9tV9HyLFi0yevfubURHRxsBAQFGjRo1jN69extLliwpcN2iko/69eu71C21sOSjsJ+Ns58pX4cOHYyoqCgjNze3mLvlfmfPnjX+/e9/GzExMUZwcLDRtm3bAglQUfdy/fr1xrXXXmvExMQYgYGBxiWXXGI888wzBe6HYdga5k6YMMGoX7++ERgYaDRq1Mh46aWXyvGTiRROs9pKhTZjxgxGjhzJO++8Q2BgILfddhu7d++mQYMGTl3nuuuuo2rVqrz//vsuxbF69WqGDh3Knj177Nvq1q3LlClTqF+/PnfeeSd79+617+vUqROjRo3itttuK3WX3MqqrD8bEfE8NTiVCmvv3r2MGzeOhIQE7rzzTvsImZs2bXL6Wj169GDMmDEux9KxY0dycnKYM2cOOTk5zJgxw96boX379gC89tpr5Obm8sEHH7B9+3b69u3r8vtVJmX92YiI5yn5kArJMAzuvvtuAgMD7b1b8oeUfuyxx5gzZw5nzpwp9fX+85//lKrBYlGCgoL46KOPmD59OjVr1uSXX36hU6dO9n1Lly5l0aJF1KxZkylTprB06VJq1Kjh8vtVJmX92YiI56naRSqk2bNn88ADD/DOO+9w55132rfPmzePJ554gqNHj3Lq1Cn7SKQiIuI5Sj5ERETEo1TtIiIiIh6l5ENEREQ8SsmHiIiIeJTXDa9utVo5dOgQ1apVw2KxmB2OiIiIlIJhGJw6dYo6derg51d82YbXJR+HDh2yz2UgIiIivuXAgQMFZlq+kNclH9WqVQNswYeHh5scjYiIiJRGeno6devWtT/Hi+N1yUd+VUt4eLiSDxERER9TmiYTanAqIiIiHqXkQ0RERDxKyYeIiIh4lNe1+SgNwzDIzc0lLy/P7FB8jr+/PwEBAerGLCIipvG55CM7O5vDhw+TkZFhdig+KywsjNjYWIKCgswORUREKiGfSj6sVit79uzB39+fOnXqEBQUpG/wTjAMg+zsbI4ePcqePXto3LhxiQPBiIiIuJtPJR/Z2dlYrVbq1q1LWFiY2eH4pNDQUAIDA9m3bx/Z2dmEhISYHZKIiFQyPvm1V9/Wy0b3T0REzKSnkIiIiHiUkg8RERHxKCUfHmIYBvfddx+RkZFYLBaqV6/O6NGjzQ5LRETE43yqwakvS0pKYt68eaxatYoGDRrg5+dHaGiofX9cXByjR49WQiIiIhWekg8P2b17N7GxsXTq1MnsUERExJcYBnx0D8R3hRZ3gH+g2RGVme8nH4YBOSYNOBYYBqUYZ2To0KG8/fbbgG22v/r16xMXF0eLFi14+eWX6dGjB/v27WPMmDGMGTMGsFXTiIiIsHslbPkQti+DxtdAeB2zIyoz308+cjJgskk/iEcPQVCVEg975ZVXaNiwIXPmzCE5ORl/f38GDhxo3//xxx/TvHlz7rvvPu69997yjFhERHyJYcC3z9iW29xVIRIPqAjJhw+IiIigWrVq+Pv7ExMTU2B/ZGQk/v7+VKtWrdD9IiJSSe36Eg5utJW0dxljdjRu4/vJR2CYrQTCrPcWEREpD1bruVKPdvdCtWhz43Ej308+LJZSVX2IiIj4lB3LIOUXCKoKnR4yOxq30jgfXiIoKIi8vDyzwxAREW9gzYNvp9iWOwyHKjXNjcfNlHx4ibi4ONasWcPBgwc5duyY2eGIiIiZtn4CR7dDSAR0HGl2NG6n5MNLTJw4kb1799KwYUNq165tdjgiImKWvFxY9XepR8dREFrd1HDKg8XwsgEl0tPTiYiIIC0tjfDwcId9mZmZ7Nmzh/j4eE0FXwa6jyIiXmzzQvh0OIRGwuhfILia2RGVSnHP7wup5ENERMRb5OXAqqm25c4P+Uzi4SwlHyIiIt7ipwVwch9UibJ1r62glHyIiIh4g9wsWDPNttx1bIUeRkLJh4iIiDfY+Dak/wnV6kDrYWZHU66UfIiIiJgt5yx894Jtudu/ILBidwZQ8iEiImK25LfgdApE1IOW/zQ7mnKn5ENERMRMWafh+xdty93/AwFB5sbjAZU2+cjIziVu3HLixi0nIzvX7HBERKSyWj8LMv6CyAbQ/Dazo/GISpt8iIiImC7jOPzwqm25x3jw9/35Xkuj0iYfedZzA7tu2HPcYd1TevTowejRoz3+viIi4iV+eAWy0iDqcmh2s9nReEylTD6Sthym94ur7etD5ybT5dmVJG05bGJUxVu1ahUWi4WTJ0+aHYqIiLjDqRRYP9u23OsJ8Ks8j+TK80n/lrTlMMMXbCI1Pcthe0paJsMXbPLqBERERCqQ1c9B7lmo2x4u6Wt2NB5VKZKPjOxcMrJzOZWZQ+LSrRRWwZK/bcLSbZzKzHF7I9QzZ87wz3/+k6pVqxIbG8sLL7zgsH/+/Pm0adOGatWqERMTw+23386RI0cA2Lt3L1dddRUANWrUwGKxMHToUACSkpLo0qUL1atXp2bNmvTv35/du3e7NXYREXGz43/Aprdty70SwWIxNx4PqxTJR9Mnv6Tpk19yxYSvCpR4nM8AUtIzuWLCVzR98ku3xvDwww+zevVqlixZwldffcWqVavYtGmTfX9OTg6TJk3i559/5tNPP2Xv3r32BKNu3bp89NFHAOzcuZPDhw/zyiuvALakZuzYsfz444+sWLECPz8/brzxRqxWq1vjFxERN/p2ClhzoVFviOtsdjQeVzma1Zrs9OnTvPXWWyxYsIBevXoB8Pbbb3PxxRfbj7nrrrvsyw0aNODVV1+lbdu2nD59mqpVqxIZGQlAVFQU1atXtx87YMAAh/f673//S+3atdm2bRvNmjUrx08lIiIuSdkCv35gW+71pLmxmKRSlHxsm9iHbRP7MG9Y21IdP29YW7ZN7OO299+9ezfZ2dm0b9/evi0yMpJLL73Uvr5x40YSEhKoV68e1apVo3v37gDs37+/2Gv/9ttv3HbbbTRo0IDw8HDi4uJKdZ6IiJhk5STAgMtvhNjmJR5eEcelqhTJR1hQAGFBAXRtXJvYiBCKqlmzALERIXRtXJuwIM8VCp05c4Y+ffoQHh7Ou+++S3JyMp988gkA2dnZxZ6bkJDA8ePHeeONN1i/fj3r168v1XkiImKC/etgVxJY/OGqx82OxjSVIvnI5+9nITGhKUCBBCR/PTGhKf5+7m3407BhQwIDA+2JAcCJEyfYtWsXADt27OCvv/5i6tSpdO3alSZNmtgbm+YLCrINt5uXl2ff9tdff7Fz504ef/xxevXqxWWXXcaJEyfcGruIiLiJYcCKibbllndArUbmxmOiSpV8APRtFsvMwa2ICg922B4TEcLMwa3o2yzW7e9ZtWpV7r77bh5++GFWrlzJli1bGDp0KH5/9+muV68eQUFBTJ8+nT/++IOlS5cyadIkh2vUr18fi8XCZ599xtGjRzl9+jQ1atSgZs2azJkzh99//52VK1cyduxYt8cvIiJu8PsK2PcD+AdD93FmR2OqSpd8gC0B+WZsd/v6vGFt+f6RnuWSeOR7/vnn6dq1KwkJCfTu3ZsuXbrQunVrAGrXrs28efP44IMPaNq0KVOnTmXatGkO51900UU89dRTjBs3jujoaEaOHImfnx+LFy9m48aNNGvWjDFjxvD888+X22cQEREXWa2w4inbcrt7IeIic+MxmcUwDM+PK16M9PR0IiIiSEtLIzw83GFfZmYme/bsIT4+npCQkDK9T0Z2rr077baJfTzaxsNs7ryPIiJSCls+hg+HQVA1eOhnqFKz1Kf6yvOquOf3hbzzE3hAWFAAe6f2MzsMERGp6PJy4dtnbMudRjqVeFRUlbLaRURExGM2vwt//Q5hNaHjCLOj8QpKPkRERMpLdgasmmpb7vovCK5mbjxeQsmHiIhIeVk/C04dgoh60OZus6PxGj6ZfHhZG1mfo/snIuIBZ/6C71+yLfd8HALVwD+fTyUfgYGBAGRkZJgciW/Lv3/591NERMrBd9MgKx1iroArBpodjVfxqd4u/v7+VK9e3T76Z1hYGJZKNg1xWRiGQUZGBkeOHKF69er4+/ubHZKISMV0Yi9seMO23Psp8POp7/rlzqeSD4CYmBiAAsOPS+lVr17dfh9FRKQcrHwarDnQ4Cpo1MvsaLyOzyUfFouF2NhYoqKiyMnJMTscnxMYGKgSDxGR8nToJ/j1A9vy1U+V+XJ51nPt9DbsOU7XxrXdPgeZp/lc8pHP399fD1EREfEuhgFfJ9qWrxgEsc3LdLmkLYdJXLrVvj50bjKxESEkJjQt1ylBypsqoURERNxl9wrYsxr8g2w9XMogacthhi/YRGp6lsP2lLRMhi/YRNKWw2W6vpl8tuRDRETEq1it8PUE23Lbe6FGfTKyc126VJ7VIHHpVgobGMEALMCEpdvo3KiWS1UwZs8Po+RDRETEHX59H1J/heAI6PZvAPuEcO5mACnpmVwx4SuXzjd7bjNVu4iIiJRVTqathwtA1zEQFmluPF5OJR8iIiJltWEOpB2A8Iug/QP2zdsm9nHtcnuOM3RuconHzRvWlnbxvpfoKPkQEREpi7Mn4LsXbMtXPQqBofZdrrat6Nq4NrERIaSkZRba7sMCxESE+Gy3W1W7iIiIlMV3L0LmSYhqCs1vc8sl/f0sJCY0BWyJxvny1xMTmvpk4gFKPkRERFx3cj+sn21b7v0U+Llv/Km+zWKZObgVUeHBDttjIkKYObiVT4/zoWoXERERV33zFORlQVxXaHy12y/ft1ksnRvVsvdqmTesrc9WtZxPJR8iIiKuOJAMWz4ELNBnMpTTRKfnJxrt4iN9PvEAJR8iIiLOMwz48lHbcos7IPZKc+PxMUo+REREnLXtU/hzAwSGlXkY9cpIyYeIiIgzcjLPTR7X+SEI992Gn2ZR8iEiIuKMDbPh5D6oFgudRpkdjU9S8iEiIlJaZ47Bmmm25V5PQlAVc+PxUUo+RERESmvVVMhKh5gr4cpbzY7GZyn5EBERKY2ju+DH/9qW+zwDfnqEusrpO3fw4EEGDx5MzZo1CQ0N5YorruDHH3+07zcMgyeffJLY2FhCQ0Pp3bs3v/32m1uDFhER8bivnwAjDy7tB/HdzI7GpzmVfJw4cYLOnTsTGBjIF198wbZt23jhhReoUaOG/ZjnnnuOV199lVmzZrF+/XqqVKlCnz59yMzMdHvwIiIiHvHHKtiVBH4BcPVEs6PxeU4Nr/7ss89St25d5s6da98WHx9vXzYMg5dffpnHH3+cf/zjHwC88847REdH8+mnn3LrrQXrx7KyssjKyrKvp6enO/0hREREyo01D778eyyPtvdArUbmxlMBOFXysXTpUtq0acPAgQOJioqiZcuWvPHGG/b9e/bsISUlhd69e9u3RURE0L59e9auXVvoNadMmUJERIT9VbduXRc/ioiISDnYvBBSf4WQCOj+iNnRVAhOJR9//PEHM2fOpHHjxnz55ZcMHz6cBx98kLfffhuAlJQUAKKjox3Oi46Otu+70Pjx40lLS7O/Dhw44MrnEBERcb+s07Bykm25238gLNLceCoIp6pdrFYrbdq0YfLkyQC0bNmSLVu2MGvWLIYMGeJSAMHBwQQHB5d8oIiIiKd9/xKcToUa8dDuXrOjqTCcKvmIjY2ladOmDtsuu+wy9u/fD0BMTAwAqampDsekpqba94mIiHiLjOxc4sYtJ27ccjKycx13ntgL/5tuW75mEgSY80U5LCiAvVP7sXdqP8KCnCoz8FpOJR+dO3dm586dDtt27dpF/fr1AVvj05iYGFasWGHfn56ezvr16+nYsaMbwhUREfGQrx6HvCyI7w5N+psdTYXiVAo1ZswYOnXqxOTJkxk0aBAbNmxgzpw5zJkzBwCLxcLo0aN5+umnady4MfHx8TzxxBPUqVOHG264oTziFxERcb8/VsP2ZWDxh75TwWIp8yUzsnNp+uSXAGyb2KfClGK4wqlP3rZtWz755BPGjx/PxIkTiY+P5+WXX+aOO+6wH/Of//yHM2fOcN9993Hy5Em6dOlCUlISISEhbg9eRESkKC4/7PNyIWm8bbnt3RDdtPjjxWlOp139+/enf/+ii58sFgsTJ05k4kQNwiIiIj5o0zw4shVCa0CP8WZHUyFpYHoREZF8Gcdh5dO25aseU9facqLkQ0REJN+qqXD2BEQ1hdbDzI6mwlLyISIiAliO7oDkN20rfaeCf+VtEFrelHyIiIhgEPT1o7ZZay9LgAbdzQ6oQlPyISIild7Vfhvx37sa/IPh6klmh1PhKfkQEZFKLYgcHg9YYFvpNBIi44s/QcpMyYeIiFRqd/t/QX2/I1irxkCXsWaHUyko+RARkUrLcuowIwM+ASCnZyIEVzU5ospByYeIiFRagd9OpIoli43WxuRdPtDscCoNJR8iIlI57VuL5dcP+F9eU0Zn/x8b9p4gz2qYHVWloE7MIiJS+eTlkvThGzyV9SqHqQnA0LnJxEaEkJjQlL7NYk0OsGJTyYeIiFQ6SUveZfjRmziM4/DpKWmZDF+wiaQth93+nueXqmzYc7xSl7Ko5ENERHxSRnZusftPZebYl9fsOkrnRrXw97OQl57KhGQ/bI9+i8M5xt9bJizdZj/eWYXNnpu05TCJS7fa1yt7KYvFMAyvSr3S09OJiIggLS2N8PBws8MREREvFTduuUvnPeC/lFl517s5mnP2Tu3nsJ605TDDF2ziwodtflozc3CrCpGAOPP8VsmHiIhUIgY55fzoO79EJs9qkLh0a4HEwxZJ2UpZCith8RUq+RAREa9TUpVKUcfkWQ0Spv/AkVNZhZ5jwaAG6RwnosTrzxrcijZxNQBo8/SKEo/3tAtLWMymkg8REfFpTZ/8slyua2DhOBFYsGKU0OfigQWbyiUGUfIhIiKVUEmJR2G2Tezj9Dkb9hxn6NzkEo+bN6wt7eIjSzyuolDyISIiXseVBz0497BPO5vDM8u3O1TRxISHMP66JlzdNLrAOa60sejauDaxESGkpGUW2u7DAsREhNC1cW2Xetb4KiUfIiLidVxtTOnsw75nkyiumPAVYEtI3J0E+PtZSExoyvAFm7CAQ0z575KY0LRSJR6gQcZERKQCyX/Yw4UjeNgam4Ljw/78h367+MhySQL6Notl5uBWRIUHO2yPiQipMN1snaXkQ0REKpQiH/bhwaY97Ps2i+Wbsd3t6/OGteX7R3pWysQDVO0iIiIVUN9msXSOyWPDK4M5TQiR7W6lU78hplZveKKUxVeo5ENERCqk0JVP0sv/J+ItKbS++rZK/bD3Nko+RESk4tm9koBtH5FnWHg0527w8zc7IjmPkg8REalYcjJh+b8AeCfvGrYYDUwOSC6k5ENERCqW716A439grRrDC7kDzY5GCqHkQ0REKo6ju+D7lwDIuXoKpwkzOSApjJIPERGpGAwDlo8Faw40voa8JglmRyRFUPIhIiIVw8+LYe93EBAK100Di3q3eCslHyIi4vsyjsNXj9mWezwCNeqbG48US8mHiIj4vm8SIeMvqH0ZdBxpdjRSAiUfIiLi2/athU3v2JYTXgb/QFPDkZIp+RAREd+Vmw2fjbEtt/on1OtgbjxSKprbRUREfNe6GXB0O4TVhN5POX16WFAAe6f2K4fApDgq+RAREd90Yi+seta2fM0zEBZpajhSehbDMAyzgzhfeno6ERERpKWlER4ebnY4IiLijQwD5t8If3wLcV1hyDJ1rTWZM89vlXyIiPiAjOxc4sYtJ27ccjKyc80Ox3w/L7IlHgEhkPCKEg8fo+RDRER8y+kjkDTettxjHNRsaG484jQlHyIi4lu+eAQyT0LMldBxlNnRiAuUfIiIiO/Y+QVs/Rgs/nD9dPBXp01fpORDRER8Q2Y6LP+XbbnTSKjTwtRwxHVKPkRExDeseArSD0KNeOg+zuxopAyUfIiIiPfbvw6S37QtJ7wCQWHmxiNlouRDRES8W04mLP27YWnLO6FBd3PjkTJT8iEiIt7tuxfg2C6oGg3XTDI7GnEDJR8iIuK9UrfC9y/alq97HkJrmBuPuIWSDxER8U7WPFt1izUXmvSHy643OyJxEyUfIiLindbPgoMbITgcrpumIdQrECUfIiLiff7aDSv+bt9x9UQIjzU3HnErJR8iIuJdrHnw6f9B7llo0ANaDzU7InEzJR8iIj4gz2rYlzfsOe6wXuGsnw0H1kFQVdsQ6qpuqXCUfIiIeLmkLYfp/eJq+/rQucl0eXYlSVsOmxhVOflrN6yYaFu+ZhJUr2duPFIulHyIiHixpC2HGb5gE6npWQ7bU9IyGb5gU8VKQApUtwwzOyIpJ5oOUESknGVk57p0Xp7VIHHpVgqrYDEACzBh6TY6N6qFv5/zVRNhQV72CFg/S9UtlYSX/eaJiFQ8TZ/8slyuawAp6ZlcMeErl87fO7WfewMqC1W3VCqqdhERkSJlZOcSN245ceOWu1yCUyJ7dUumqlsqCZV8iIiUs20T+7h03oY9xxk6N7nE4+YNa0u7+EiX3sMrqLql0lHyISJSzlxtW9G1cW1iI0JIScsstN2HBYiJCKFr49outfnwCqpuqZRU7SIi4qX8/SwkJjQFbInG+fLXExOa+m7ioeqWSkvJh4iIF+vbLJaZg1sRFR7ssD0mIoSZg1vRt5kPDzuu6pZKS9UuIiJerm+zWDo3qmXv1TJvWFvfrmoBOLpT1S2VmEo+RER8wPmJRrv4SN9OPPJy4OP7bNUtDXuquqUSUvIhIiKetfo5OLwZQqrDP15XdUslpORDREQ850AyfDfNttz/RQj34TYr4jIlHyIi4hnZZ+CT+8CwwhUDodkAsyMSk5Qp+Zg6dSoWi4XRo0fbt2VmZjJixAhq1qxJ1apVGTBgAKmpqWWNU0REfN1XT8DxP6BaHbjuebOjERO5nHwkJycze/ZsrrzySoftY8aMYdmyZXzwwQesXr2aQ4cOcdNNN5U5UBER8WG/fQ0/vmVbvuF1CK1hbjxiKpeSj9OnT3PHHXfwxhtvUKPGuV+gtLQ03nrrLV588UV69uxJ69atmTt3Lv/73/9Yt26d24IWEREfknEcloywLbd/ABpeZW48YjqXko8RI0bQr18/evfu7bB948aN5OTkOGxv0qQJ9erVY+3atYVeKysri/T0dIeXiIhUEIYBn42G06lQ6xLoPcHsiMQLOJ18LF68mE2bNjFlypQC+1JSUggKCqJ69eoO26Ojo0lJSSn0elOmTCEiIsL+qlu3rrMhiYhIOcmznptVZsOe4w7rpfLL+7BtCfgFwI2zITDUzRGKL3Iq+Thw4AAPPfQQ7777LiEhIW4JYPz48aSlpdlfBw4ccMt1RUSkbJK2HKb3i6vt60PnJtPl2ZUkbTlcuguk/QmfP2xb7v4IXNSqHKIUX+RU8rFx40aOHDlCq1atCAgIICAggNWrV/Pqq68SEBBAdHQ02dnZnDx50uG81NRUYmJiCr1mcHAw4eHhDi8RETFX0pbDDF+widT0LIftKWmZDF+wqeQExGqFT4dDVhpc1Aa6jC3HaMXXODW3S69evfj1118dtg0bNowmTZrwyCOPULduXQIDA1mxYgUDBtj6b+/cuZP9+/fTsWNH90UtIiIlysjOdem8PKtB4tKtFFbBYmCbUXfC0m10blSryGHeA9a9RtCeNRgBoZzo8yqtHvsSgB8f70Wtqu4pORff5VTyUa1aNZo1a+awrUqVKtSsWdO+/e6772bs2LFERkYSHh7OqFGj6NixIx06dHBf1CIiUqKmT35ZLtc1gJT0TPtEdxdqZvmDj4MmggUePXs7i17fa9/X5ukV7J3ar1ziEt/h9lltX3rpJfz8/BgwYABZWVn06dOH119/3d1vIyJSqYQFBfjEQzuMTF4NfI0gSx5f5LVlUV5Ps0MSL2QxDMPJpsvlKz09nYiICNLS0tT+Q0SkDFytdtmw5zhD5yaXeNy8YW1pFx/psC3os1EE/LIQa7U6ZN6zBkJrkJGdS5unVwCqdqnInHl+u73kQ0REvENYkGt/4rs2rk1sRAgpaZmFtvuwADERIXRtXNuxzceWj+CXhYAFvwFvEBZR220xScWiieVERMSBv5+FxISmgC3ROF/+emJCU8fE48Q+WDbGttz1XxDXpdzjFN+l5ENERAro2yyWmYNbERUe7LA9JiKEmYNb0bdZ7LmNebnw8X22brUXt4Ue4zwcrfgalX+JiEih+jaLpXOjWvZeLfOGtS1Y1QKw5nk4sA6CqsGAN8E/0IRoxZeo5ENERIp0fqLRLj6yYOKx73+w5jnbcv+XoEacS++TkZ1L3LjlxI1b7nJDWfEdSj5ERMQ1Z0/AR/eCYYXmt8GVA82OSHyEkg8REXGeYcCy0ZD+J9SIh+ueNzsi8SFKPkRExHmb3oFtn9pmqx3wFgRXMzsi8SFKPkRExDkpW+CL/9iWr3oMLm5tbjzic5R8iIhI6WWdgg+GQG4mNLoaOo82OyLxQUo+RESkdAwDPhsDf/0O1erAjbPBT48RcZ5+a0REpFT8N8+HXz8Aiz/c/F+oUrNU5+VZzw3SvmHPcYd1qZyUfIiISIkus+wj6Ku/Ry7t9QTU71iq85K2HKb3i6vt60PnJtPl2ZUkbTlcHmGKj1DyISIixarCWV4LfBVLXhY0vgY6PVSq85K2HGb4gk2kpmc5bE9Jy2T4gk1KQCoxDa8uIiJFysjKYXLgWzT0O0xetTpk9XsNcq2Atdjz8qwGiUu3FjorroFtgroJS7fRuVEt/P0sDqOalmaEU82O69sshmF4VeVbeno6ERERpKWlER4ebnY4IiKV2vjHxjIl8C1yDT9uyX6CjcalZocEwN6p/cwOQS7gzPNb1S4iIlK4lF+ZEPAOAM/l3uI1iYf4PpVbiYhIQZnp8P4Qgi05rMhryRt5/fjx8V6lru7YsOc4Q+cml3jcvGFtaRcfSUZ2Lm2eXgHg1PuIb9JPV0REHBkGLHsQju+G8Ivp9cDH7AmLdOoSXRvXJjYihJS0zELbfViAmIgQujaujb+fxaH77ZaD6fbtUjGp2kVERBytex22fmKbt+Xm/4KTiQeAv5+FxISmgC3ROF/+emJCU/z9LOqOWwkp+RARkXP2fg9fPWFb7jMZ6rV3+VJ9m8Uyc3ArosKDHbbHRIQwc3Ar+jaLVXfcSkq9XURExCb9EMzuBmeOwhWD4KY5YCm56qOkrrGnMnNoP3klALMGt7J3r82zGvR+cXWBxCOfBYgOD+Hrsd1cqoJRuxHPcub5reRDREQgNxvmXQd/JkN0M7j7awgKK9WpceOWl3NwrlF3XM9SV1sREXHOl+NtiUdIBNwyv9SJh4grVCYlIlLZbV4EyW/alm96AyIbOHX6tol9it1fVDdaZ7vjSsWh5ENEpDI7/DN8Ntq23P0RuKT4RKIwzrStCAsKsB/vbHdcqThU7SIiUlllHIf37oTcTGh0NXQf59G3d6Y7rlQsSj5ERCojax58fC+c3AfV69t6tvh5/pFQmu64UvGo2kVEpDJaNRV+/wYCQuCWBS4NJOYufZvF0rlRLa6Y8BVga+OhqpaKTSUfIiKVzfZlsOY523LCKxB7pbnxgEOi0S4+UolHBafkQ0SkMkndCh/fb1tudz80v9XceKRSUvIhIlJZnPkLFt0GOWcgvhv0ecbsiKSSUvIhIlIZ5OXAB0NsDUxrxMHAt8E/0OyopJJS8iEi4kUysnOJG7ecuHHLS5wzxSlfPgZ7v4OgqnDrIlMbmIoo+RARqeg2vQMbZtuWb5oD0U3NjUcqPXW1FREpo4zsXJo++SVgG2rcq2ZT3b8OPhtrW77qMWji+cnWwoICNMmbOFDJh4hIRZX2J7w3GKw50PQf0O1hsyMSAZR8iIhUTFmnYOEtcOYoRF8BN8wEi8bOEO+g5ENEpKKx5sFH90DqFqgSBbctgqAqZkclYqfkQ0SkovnqCdiVZBs6/bbFUL2u2RGJOFDyISJSkfz4X1g3w7Z8w0y4uLW58YgUQsmHiIiPKXIskN3fwvJ/25avehya3WROgCIlUPIhIlIRHN0F7w8BIw+uvAW6/dvsiESK5EWd0UVExCVn/oKFgyArDep2gOun+1zPFo0FUrmo5ENExJflnIXFt8GJPVC9Ptz6LgQEmx2VSLGUfLig3OZeEBFxggUrQUuHw4H1EBIBt78PVWqZHZZIiVTtIiLiox4LeJeAnV+AfxDcuhCimpgdkkipqORDRMQHDfP/gnsCvrCt3DAT4rqYG5CIE5R8iIh4kTyrYV/esOe4w3o+/x3LeCJgAQDZVz0JV9zssfhE3EHJh4iIl0jacpjeL662rw+dm0yXZ1eStOXwuYP2rydo6QP4WQzm5/Ymt8ODJkQqUjZKPkREyqg0pRUlSdpymOELNpGanuWwPSUtk+ELNtkSkGO/w6JbseRm8k1eSybkDvG5LrUioAanIiJlkrTlMIlLt9rXh85NJjo8mEevu4yrm0aX6hp5VoPEpVspLGUxAAswYcmv9Ap7mMCzx8mJbsGofaPIw79UPe7CgvSnXryLfiNFRFyUX1pxYdKQmp7FQ4s3u+19DCDlVA4/ZoZSxxLFgH33c5YQANo8vaLE8zV4l3gbJR8iUqFlZOfS9MkvAdg2sU+BUgBXx+oprrSivOw26vBI7n0cI8KD7yrifko+RKRSy09MfEHdfg+T1LotGdm59hKPHx/vpWoV8Tn6jRUR8XIWrMRU8adL+w74+zk2MA0LClDyIT5Hv7EiUqltm9jHpfM27DnO0LnJJR43b1hb2sVHlnjc19tSGb14c4FqHAtWwELijS3ticeFvWu6Nq5dICkR8WZKPkSkUnO11KBr49rERoSQkpZZaLsPCxATEVLqxOAfLS4iOMCPxKVbHbrbxoQaJA5oQ99msUDhvWtiI0JITGhqP0bE22mcDxERF/j7WUhMaArYEo3z5a8nJjR1qkSib7NYvr0mlUWBk3glcDrz2/zB908kOCQeJY4FIuIDlHyIiLiob7NYZg5uRVS44xT2MREhzBzcqtCSiIzs3CJfWb98QujykXT0385xI5zL+txDVm4eGdm5nMrMKXYsEIAJS7dxKjOn2Pco6iXiSap2cYHqW0UkX99msXRuVIsrJnwF2Np4FPc3oajeNd38fubNwGlYLFbez+3OxNw7eeqZlaWOwwBS0jPtcThLY4GIJ6nkw0mlmntBRCqV8xONdvGRTn8ZaWfZzuzAlwiy5PFZXnvG5d6LoT/PUoGp5MMJRY1mmF/fWlQxq4hIvgt71/gd/ongd+/Dkp1NXqNr6JDwFtYp3wGOY3i4u3eNiJkqXfJRHqMZ2udeWLqNzo1quVQFo376IpWDw//11G2weCBkn4a4rvjf8g5hRqDDsfnHu7t3jYiZKt0Tr7xGM1R9q4g45djvMP8GOHsCLmoDty2CwFAo4gtSfu+a4Qs2YQGHBMTV3jUiZlGlooiIp/21G97uD6dTIboZ3PEBBFcr8TRXeteIeKNKV/LhLaMZikgldfwPeDsBTh2G2pfBP5dAWOn/Zjjbu0bEGzlV8jFlyhTatm1LtWrViIqK4oYbbmDnzp0Ox2RmZjJixAhq1qxJ1apVGTBgAKmpqW4Nuizy61CdfeXXtxb139sCxP5d3+rK9UWkfFzYNf78dY87sRfmJUD6Qah1KQxZClVqOX2ZsvauETGbU8nH6tWrGTFiBOvWrePrr78mJyeHa665hjNnztiPGTNmDMuWLeODDz5g9erVHDp0iJtuusntgXtaWUYzzMjOJW7ccuLGLddgPiIe5FVd40/u/zvx+BNqNoYhy6BqlOfjEPECTn3lTkpKclifN28eUVFRbNy4kW7dupGWlsZbb73FwoUL6dmzJwBz587lsssuY926dXTo0MF9kZsgv761wNwLmldBxOt4Vdf4kwdgXn9I2w+RDW2JR7Voz7y3iBcqU3l/WloaAJGRtvrKjRs3kpOTQ+/eve3HNGnShHr16rF27dpCk4+srCyyss49yNPT08sSUrlTfauIZ7lSWuhVXePTDtraeJzcBzXiYehnEK4vKlK5uZx8WK1WRo8eTefOnWnWrBkAKSkpBAUFUb16dYdjo6OjSUlJKfQ6U6ZM4amnnnI1DFOovlXEc8qje7zHusaf2Hcu8ahe/+/Eo45L7ylSkbicfIwYMYItW7bw/ffflymA8ePHM3bsWPt6eno6devWLdM1RURM99duePt6WxuPGvG2qpaIi0s8LSwoQOP+SIXnUvIxcuRIPvvsM9asWcPFF5/7zxQTE0N2djYnT550KP1ITU0lJiam0GsFBwcTHBxc6D4REVe6x5veNf7oLnjnelt32pqNbb1aVOIhYudU8mEYBqNGjeKTTz5h1apVxMfHO+xv3bo1gYGBrFixggEDBgCwc+dO9u/fT8eOHd0XtYhUGq50Rff0UOQOpRWp22yJx5mjtnE8hixVrxaRCzj1v3rEiBEsXLiQJUuWUK1aNXs7joiICEJDQ4mIiODuu+9m7NixREZGEh4ezqhRo+jYsaPP93QREd9h2lDkh3+Gd26As8ch5gq4cwlUqene9xCpAJwa52PmzJmkpaXRo0cPYmNj7a/33nvPfsxLL71E//79GTBgAN26dSMmJoaPP/7Y7YGLiBTH40OR/7nR1rj07HGo0wr+uVSJh0gRnK52KUlISAgzZsxgxowZLgclIuIOHusa/8dqWHy7bXbauu1tc7WERLj3PUQqEI3rLSIVWrl3jd/+GXw4DPKyIb473LoQgqu69z1EKhglHyIirvrpXVg6EgwrNOkPN/8XAtR7T6QkSj5ERFyx9nX4crxtucVgSHgF/D3zJ1VjgYivU/IhIuIMw4Bvn4E1z9vWO46Ea54Gi0Y6Fiktp3q7iGu8akpvEXGdNQ8+//e5xKPn40o8RFygko9ylrTlMIlLt9rXh85NJlaz4Ir4npyz8NE9sOMzwAL9pkHbe8yOSsQnKflwQWnrW71qSm+RCiAjO9c+0dy2iX1cGv3UtTc+DgtvgT83gH8w3DQbLr/RM+8tUgEp+SiBK9N5g5dN6S1SSZRLcnJiLywYAH/9bhu747bFUL9T2a8rUonpCVaC8pjOGzw4pbeIuO7QT/DuIDhzBCLqwh0fQlQTs6MS8XlKPkRECvPbN/D+PyHnDERfYRu1NFzVpCLuoOSjBK5M5w1eMKW3iLgu+S34/GEw8qBBDxg0H0LCzY5KpMJQ8lECV+uMPT2lt4i4gTUPvnwM1s+0rTe/DRJehYAgc+MSqWA0zkc5yZ/SG85N4Z2vXKf0FhHXZKbDolvPJR49n4AbZirxECkHSj5cVYoZfj0+pbdIBZSRnUvcuOXEjVvucu+zEp3cD//tA799BQGhMPBt6PZvDR4mUk5U7eKK3Cxbn/82d0HT64s91GNTeotIoUocl+dAMiy+Dc4charRcNsiuKi15wIUqYRU8uGK5Dfhj2/h/Tth5TNgtRZ7eLlP6S0irvl5MczrZ0s8Yq6Ae1cq8RDxACUfrmh3P3QYYVte8xwsvt1WXywiviEvB74YB5/cD3lZcOl1MCwJIi42OzKRSkHJhyv8A6DvZLhhlm2o5V1fwJu94NhvZkcmIiU5cwzm33iuYWm3/8At70JwVXPjEqlElHyURYvb4K4kCL8Iju2CN3rCrvIZEVVE3ODQZpjTA/Z+B0FV4ZYF0PMx8NOfQhFP0v+4srqoFdy3Cup1hKx0W0PUNc+X2A5ERDzsl/dtPVrSDkBkQ7hnBVyWYHZUIpWSkg93qBoF/1wKbe4GDFj5tK31fMZxsyMTqVDyrOe6uG/Yc9xhvUi52fDFI/DxvZCbCY2vsTUs1RwtIqZR8uEuAUHQ/0W4fvrf7UCSYHZ3OLjJ7MhEKoyE6T/Yl4fOTabLsytJ2nLYvq1AcnJ8H8ztC+tn2TZ2/ZdtVtrQ6p4KWUQKoXE+3K3VPyG2hW1CqhN74L99COj9NHARBcc6rRjKZRpzkUIcOZXlsJ6SlsnwBZuYObgVAIlLt9r3DZ2bTKzlBIkBfvStUh1unAWXXuvJcEWkCHpKlIfYK+H+1fDp/8GOzwj68j+8GtiRcTn3mh2ZiGlcHZ30VGZOkfsMbCn9uI9+Je1sToF5lFKMCIbnjOblrhdzdXwzKCQGJcsinqf/deUlJMLWkn7tDIxvErmetTS17MNypBFc3Nzs6ESKVR6lWfnXczcDOHm28ATFwA8weGjZQVh2sNBjSpq5WsmJiPvpf1V5slig00iyolty4p3BNPI7hDH3arhmErS7T/NGiHhE8f/PSkqKih2aXURcouTDA6x129MvazLPBc6mNz/BF/+B3SvhHzOgSi2zwysXagciFyqphKG8SkZExPvoieAhxwnnnpx/s+vaAwStSLT1hpnZ2dYIruFVZocnUu5KSkCLSk5OZebQfvLKYs7Mb/lRvHnD2tIuPrLE40Sk/Cn58CgLuW3uJahBV/jwLji2E+bfAJ0ehJ5P2LrrilRSRSUnxY3lYcGKgYXqltOkGVULNDi1HQMxESGaTVrEi2icDw/In9J779R+tj+wMc1so6K2uct2wP9etQ3NnrLF1DhFvE3SlsP0fnF1kftjOM6sequYOqAFULD8I389MaGpEg8RL6LkwyxBYdD/JVuPmNBISP3VNufEdy9AnmtdEkUqkqQthxm+YBOp6VkX7LECBmMCP+H7G7LoO/x5+ra5lJmDWxEVHuxwZExECDMHt6Jvs1iPxS0iJVO1i9kuS4C67WHZQ7Dzc1gxEXZ8bmsLUqsxoMab4rtcHdsjz2qQuHRrodUo/N19dmHQzdzVohf+OXkAdLukNktHdra3D5k1uBWdG9XC38/iEIf+/4iYT/8LvUHVKLh1Ify82DYHxcEfYVYX6D0B2t1vdnQiLiu/HiwWUs/kccWEr4o84oEFhU9toK6zIuZTtYu3sFigxW3wf/+DBlfZJsBKGgfzrsNybKfZ0YmIiLiNSj68TcTFcOcnsHEufPk47F9LyFs9eND/embmXW92dCJOKWlsj0IZBj+u+IB/roko8dALu89mZOfS5ukVAPz4eC9VsYh4Kf3P9EYWi60nTKPesPxfWH77irGBH9Lffy1+f9aGBp3MjlDKSUVr3+N0/Md+h89G03nP98TyKilEYhQxhkdMePHdZ8OCAnz+/olUVKp28WbV68Ht75P1jzkcM8K5xO8gwe9cB8v/BZnpZkdnV2Aa82LGZRApVG42rH4eZnaCvd/hHxhCYptswFLk8GHjr2ui7rMiPkrJh7ezWMi7fAC9s57ng9xuWDAg+U2Y0Q5+/RAMcx/0F47DMHRuMl2eXcnX21JNjEp8yt4fYHZX+PZpyMuChr1gxDr63nxPod1n813dNNrDgYqIuyj58DIZ2bnEjVtO3LjlDt0DT1KNh3MfIPP2TyCyAZw6DB/dDfP6QepWU2ItahyGlLRMRi/ebEpM4kPSDtpG+p13HRzdAWG14KY3YfBHUCMOgL7NYvlmbHf7KbMGtzIpWBFxJ1WI+hhrXDcYvhbWToc1L8C+H2BWV9ssuT3GQWh1p65XHuMwXLjNlfdQXb25LqxKc+vQ5LlZsHYGrJkGOWcAC7QeCr2ehLCCc6+c/75t4mq4JwYRMZX+wvuiwBDo9jBceSt8+ShsXwrrZ8KWD21jgzS/HfxKV6jliZlE83sfOMNTYzG40sCzojUKvVDSlsMkLj1XmjZ0bjKxESEkJjQt+0ihu76ydSE/vtu2Xrc9XPsc1GlRtuuKiE9RtYsvq14Xbplv65pbszGcOQpLRsCcbvDHKrOjEx9UXFXa8AWbSNpy2LULp26FBTfDwoG2xKNqNNw4G+76UomHSCVUsb6yVVYNe8Lw/8H6WbDmeUj5Fd75BzS6Gq6eCNFNizzVpXEYsBXFD52bXOrj8wvOX761hRoKlrPyqkqzABOWbrMPWV4allOHCFwzFf9fFmExrOAXAB2GQ7f/QEi4S3GKiO9T8lFRBARB5wehxR2w5jlbj5jfv4bdK6DlYLjqMagWU+C0wqoMSvPwal2/BtHhwRxJzypi/g1H+Q+vKZ/voGeTqBIfXoXFUNGqN8pLeVWlGUBKemaxQ5rnq0oGDwQs427/LwiwZP8d2D+gVyLUbFgu8YmI79Bf84qmSk249llbA9RvJtjag2x6x9Ytt/390OnBQhv1nc8bHl6F0Zwc3i+YbG7zX8mogE+oaTkFQLL1Eibn3MEng0abG5yIeA0lHz7ApZ4HNRva2oPsXw9fPQZ/JsP3L8GGN6Hj/0GH/3O6Z4z4hvKuSrtwSHMAcrMI2DyfgLUv43fK1i7EGtmInJ6JXN74Wt61uKenTFhQgJJQkQpAyYeXK3PPg3rt4e6vYecX8O1kSP0VVj9rax/ScRR0eACCqzmc4uzD61Rmjn0a89Io9OElbuNq9VTXxrWJjQghJS2z0Ko0CxATccGQ5rnZ8NN8+O4FSD9o2xZ+EXT7N34t7yTYP9ClWFyl5ETENyj58GJfb0tl9OLNBR4E+T0PZg5uVboExGKBJtfBJX1hxzL4dgoc3W4bUXLd67aSkLb3QKhtDIWytK0orh1IoQ8v8biiugr7+1lITGjK8AWbsOA4Xkv+TysxoantZ5ebBZsX2pKOtAO2ndVioeu/oNU/IaDwUUlFREDJR7lztufB+cc/s3y7W3seANCoHzToS9hvy2DVFPjrd1j5NHz/sm2gp44jILyOUzGf79HrLmP04s0lP7zEK/VtFsvMwa1IXLrVobttTH5pW6MwW/Xduplw+u8h9KvGQNex0GqIbQwaEZESKPkoZ2VpvHnkVFaR+8reePNmaHoDbP3E9jA5shXWvgbrZ0PzW6HzQ1CrsdPXvbppdPEPr7IOUiXlrm+zWDo3qmX/3Zo3rC1do3Px3/A6LJ0L2baGpIRfBB1HQpthEBhqYsQi4muUfFRm/gFw5UC44mb47Wv44WXbcO0/zYefFsCl10H7+yC+u63qppQKfXipqqVUynVYcyfkv2cTy366bPsc//ffh7y/u8zWbmJLTpvdbOviLSLiJCUf5czZxpsZ2blODUfulsabFgtcco3ttX+9LQnZ+TnsXG571W5i67p75S0QXLXYSzV98ku2Tezj8MBsFx+pxKMUnGlcXK5DvOfl4L9tCe8FPU97vx3w89/b63WyJR2Nryn18P0iIoVR8lHOnH0onP/Nt0ZYICczcjzbeLNee6i3CI7uhA1zYPMi24yjy8fCN09ByztsjVM1UJRb5Q9rXubGxWVxKhU2zoONcwk+dZj2fpBr+MFlCQR0GmH73TCBerCIVDxKPrzIhd98T2TkFHqcRxpv1r4U+r1gm2l080JbInL8D1vvmHWvQ/0utpFTm14PqGcDeHZY8/Pfq6T3LTYBtubCzm9g87u27thW2++cUSWKV9O6sDC3J9/edBsBGl1WRNxIf1G8RFHffAvj0cabIRG2uTja3Q+7V8KG2bb2Ifu+t70+f5igpjfS0tKQn4xGnEuNfIM721iYNTJsSdV0hZUaNLQcZKD/akKnj4Ezqed21O0A7e7lbKPreOmpb8sStohIkZR8uFF5fPO90KzBrezfgEv7fm5pD+DnB417215pf8LPi2yNUk/sJWDzO3wSDL9ZL2JpXkcsfzXEqNmo7O9ZzlwdwM1bGoU67VQqbF9K8OZFrAjeaNt2BgirZWvP0+J2iGlm2+7i77KISGlYDMMozTPPY9LT04mIiCAtLY3wcN+a9TJu3HKzQyhUudWXW62w7wdyN75Dzq+fEpo/gRhgjb6S5/5symfWjnz11GCvmxSuqJKm/BSiqDYW+QnL+d2I8xOWbpfUdikWV4Y1P79h8o+P9yr6/p4+Qtjuz2Hrp7aeTIYVsLXl+Nbakm6DHiL4smsL9Fop1watIlIhOfP81l8UcZ2fH8R3JfuijrT/8Rr6+CfT328d3QO24Jf6C+MCf2Eci8mb9w5c1h8uvdbWc8ZN83x4eur4kkacffnWFlzdNNrpB7VLw5qfJywowPE90/6EXUmwbQns/d6ecABwURuym/yDTstrcYwItl3aBwL0Z0BEPEt/ddyovCf0ghK+5ZroFGF8mNedD/O6s+1fbfDf+Rk/fvYWHfy24X9oIxzaCCuegur1bMO8X9IX4rqUaRhub5l9Nz9heGjxZsC134NSD2teGGse7P/RlnDs+so2YNz56rSCy2+0TWlfoz652bkcW14+905EpDS87ynmw8pzQq/87QW+5XqjsJrktRzCHZ/EUJuT/PCPMwT98TX8sRpO7rf1nNkwBwKrQHxXiO9me0VdXiHGj3AlKdo7tV/pR4Y1DCwn9jDI/1s6+W0l9JURcPbEuf0WP7i4ra2k6fIboUZcGT6NiIj7eflTrHIozYRevuoo1cltfQtBHe+D7DO2BGRXEuz6Ek6n/L2cZDs4NNJWGhLfDep3tnX39fMv8tqulDCUV2mJOxQ5MqwFOLHXVoWy5zvY+z2h6X/yXP6EsWex9UpqdDVc0gca9YYwzRosIt5LyYeXKG5Cr3HXNrEX6XujsKAAtk3sU/KDPaiKbXbdJtfZGqum/AJ71the+/4HZ4/D9qW2F0BQVajTEi5qDRe3gYvaQPi5EgBXSoC2TexDntWg94uri519Nzo8hK/HdrNXdbjSKNQV/n4WqnOK5n5/0OnPX/DfuBkOboQzRx2OM/wCSc5twHrrZdw79B5C4jvahssvBZ/trSMiFYaSDy9S1DffrNw8kyMrB35+UKeF7dX5QcjLgYObYO/fycifGyH7NOz9zvbKVy0WoppC1GXn/q19qS2xKYX8hOWp6y8vtqRpwvVNqRYSaN9e1kahhcrNgmO/2UaQPboDjmwnJGULm0P22vaf97HxC7AlX3FdIK4LZ2NaM2iS7YC765U+8XC1e7GIiDsp+fAylXZOFP/Av4d2bw/dHrY1ojy6A/78EQ7+aEtMjmyDU4dtr93nD6xlgRr1bW0basRB9fq29epxtn/DahboYVPi1PEXPIhLUzVWoFGoYdjaYpzcX/D11++2EWMNx8Qyv8XLbmss9a/sQsDFbWwlPzFXOE5X70JPH2eGcNeQ5iJSnpR8SLko88PLzx+iL7e9Wg+xbcs6Dalb4eh2OLLdlowc2W6rkjix1/Yq9FoBUKU2VKkFVaLsy31DqtO9ayj/98UJ0oyqjO7ViM5x/vgH7YP9+89LWCxgsdA3wsrMq4OZ8EMmKRnnkoyYkBwSG+2m72/LYPMxyDgGZ/6y/ZuTUfznDI6AqCa2LshRl5FZozHt5h4nnSpsu75PocOaZ2TnOjW8Orjevbi0vL4RtIh4Ff3FELfwSDuC4KrnSkfOd/ooHNsJJ/bByX3n/bvXVkpizT1XYnKBUGBu/vha3//9KkZf4GrDwobAJhyhOlGcpJ2xA//fix6rz6gShaV6PVs34/xXZDzUvgyqxTiUylizc0mn+LYzF7atcWYW5CJjxLnuxRdSKYmIOEPJh5SZ6e0Iqta2veK6FNyXmwVnjtlKR85/nT4CWenkZp7h21/3EEYWHeoG459zFvKybFUmABjnli0WCAyDgBA4kElVMml7+SX4B7e09S4Jq2krXQmrRWZQDfq8sY0UI5KfHk1QyYCIyHnK7S/ijBkzeP7550lJSaF58+ZMnz6ddu3aldfbiUm8Yir44gQEQ8RFtlchsrNzuXfT38OIDy3dMOJZ2bnclj/0+E2FV41Ys3PZZxwvQ+CFS9pyuECbEzjX7iR/lNULeaq3johIaZRL8vHee+8xduxYZs2aRfv27Xn55Zfp06cPO3fuJCoqqjzeUsrA08OUl1ZFLS0o7n4X15ajNPd7yuc76NkkqsD9bl2/BtHhwcV2L3a6t46IiIvK5a/7iy++yL333suwYcMAmDVrFsuXL+e///0v48aNczg2KyuLrKxzvQ3S09PLIyQphrcMU36hitqOoLT329m2HGW53wYlDOEuIuJGbh/LOjs7m40bN9K7d+9zb+LnR+/evVm7dm2B46dMmUJERIT9VbduXXeH5PPye47sndqvwpYGiPk0zoeIeIrbn2THjh0jLy+P6GjHeufo6Gh27NhR4Pjx48czduxY+3p6eroSEA8r7wnx1I7AUXH3OyM7117iceEkgu6636cyc2g/eSUAswa3crlaTETEVaZ/jQ4ODiY42PWZTaXsynNCPLUjKKi09/vCSQTL4353u6S2StNExOPcXu1Sq1Yt/P39SU1NddiemppKTEyMu99OTJQ/6icUnACvVFPBi1N0v0WkonB78hEUFETr1q1ZseJcYzmr1cqKFSvo2LGju99OTJY/THlUuGPpVUxEiPndbH1QSe17dL9FpCIol/LWsWPHMmTIENq0aUO7du14+eWXOXPmjL33i1QsRU4Fb+I38IzsXHuvkm0TSzd+h6/wxvstIuKMcvmLfMstt3D06FGefPJJUlJSaNGiBUlJSQUaoUrFUWknxDOJ7reI+LJy+zo4cuRIRo4cWV6Xr7A0m6iIiFR0bm/zISIiIlKcilMRLuICV0qaSnOOR2b5FRHxUSr5EHGzpC2H6f3iavv60LnJdHl2JUlbDpsYlYiI91DyIeJG+bP8pqZnOWzPn+VXCYiIiKpdRArQLL8iIuVLf81ELqBZfkVEypeqXURERMSjVPIhcgHN8isiUr6UfIhcoKLP8quB7ETEbKp2EXETzTorIlI6Sj5E3EizzoqIlEzVLlIhmTnCqGadFREpnpIPcQtvakeQtOUwiUu32teHzk0mNiKExISmHit5KO9ZZ73pfouIOEvVLlKhaIRRERHvp5IP8TreOMKoRhcVEXEf/UUVr+ONI4yqikNExH1U7SIiIiIepZIP8ToaYVREpGJT8iFep6KPMCoiUtmp2kUqDI0wKiLiG5R8SIWiEUZFRLyfql2kwtEIoyIi3k0lH1IhlfcIoyIi4jolHyIiIuJRSj5ERETEo9TmQ6QcaOI3EZGiqeRDREREPErJh4iIiHiUkg8RERHxKCUfIiIi4lFKPkRERMSjlHyIiIiIRyn5EBEREY9S8iEiIiIepeRDREREPEojnEqFpBFGRUS8l0o+RERExKOUfIiIiIhHKfkQERERj1LyISIiIh6l5ENEREQ8SsmHiIiIeJSSDxEREfEoJR8iIiLiUUo+RERExKOUfIiIiIhHKfkQERERj1LyISIiIh6l5ENEREQ8SsmHiIiIeJSSDxEREfGoALMDuJBhGACkp6ebHImIiIiUVv5zO/85XhyvSz5OnToFQN26dU2ORERERJx16tQpIiIiij3GYpQmRfEgq9XKoUOHqFatGhaLxdRY0tPTqVu3LgcOHCA8PNzUWHyB7lfp6V45R/er9HSvnKP75Zzi7pdhGJw6dYo6derg51d8qw6vK/nw8/Pj4osvNjsMB+Hh4fqldILuV+npXjlH96v0dK+co/vlnKLuV0klHvnU4FREREQ8SsmHiIiIeJSSj2IEBweTmJhIcHCw2aH4BN2v0tO9co7uV+npXjlH98s57rpfXtfgVERERCo2lXyIiIiIRyn5EBEREY9S8iEiIiIepeRDREREPErJh4iIiHiUkg8XZGVl0aJFCywWC5s3bzY7HK+zd+9e7r77buLj4wkNDaVhw4YkJiaSnZ1tdmheY8aMGcTFxRESEkL79u3ZsGGD2SF5pSlTptC2bVuqVatGVFQUN9xwAzt37jQ7LJ8wdepULBYLo0ePNjsUr3Xw4EEGDx5MzZo1CQ0N5YorruDHH380Oyyvk5eXxxNPPOHwN33SpEmlmkCuKEo+XPCf//yHOnXqmB2G19qxYwdWq5XZs2ezdetWXnrpJWbNmsWjjz5qdmhe4b333mPs2LEkJiayadMmmjdvTp8+fThy5IjZoXmd1atXM2LECNatW8fXX39NTk4O11xzDWfOnDE7NK+WnJzM7NmzufLKK80OxWudOHGCzp07ExgYyBdffMG2bdt44YUXqFGjhtmheZ1nn32WmTNn8tprr7F9+3aeffZZnnvuOaZPn+76RQ1xyueff240adLE2Lp1qwEYP/30k9kh+YTnnnvOiI+PNzsMr9CuXTtjxIgR9vW8vDyjTp06xpQpU0yMyjccOXLEAIzVq1ebHYrXOnXqlNG4cWPj66+/Nrp372489NBDZofklR555BGjS5cuZofhE/r162fcddddDttuuukm44477nD5mir5cEJqair33nsv8+fPJywszOxwfEpaWhqRkZFmh2G67OxsNm7cSO/eve3b/Pz86N27N2vXrjUxMt+QlpYGoN+lYowYMYJ+/fo5/I5JQUuXLqVNmzYMHDiQqKgoWrZsyRtvvGF2WF6pU6dOrFixgl27dgHw888/8/3333Pttde6fE2vm9XWWxmGwdChQ3nggQdo06YNe/fuNTskn/H7778zffp0pk2bZnYopjt27Bh5eXlER0c7bI+OjmbHjh0mReUbrFYro0ePpnPnzjRr1szscLzS4sWL2bRpE8nJyWaH4vX++OMPZs6cydixY3n00UdJTk7mwQcfJCgoiCFDhpgdnlcZN24c6enpNGnSBH9/f/Ly8njmmWe44447XL5mpS/5GDduHBaLpdjXjh07mD59OqdOnWL8+PFmh2ya0t6r8x08eJC+ffsycOBA7r33XpMil4pgxIgRbNmyhcWLF5sdilc6cOAADz30EO+++y4hISFmh+P1rFYrrVq1YvLkybRs2ZL77ruPe++9l1mzZpkdmtd5//33effdd1m4cCGbNm3i7bffZtq0abz99tsuX7PSz+1y9OhR/vrrr2KPadCgAYMGDWLZsmVYLBb79ry8PPz9/bnjjjvK9EPwFaW9V0FBQQAcOnSIHj160KFDB+bNm4efX6XPdcnOziYsLIwPP/yQG264wb59yJAhnDx5kiVLlpgXnBcbOXIkS5YsYc2aNcTHx5sdjlf69NNPufHGG/H397dvy8vLw2Kx4OfnR1ZWlsO+yq5+/fpcffXVvPnmm/ZtM2fO5Omnn+bgwYMmRuZ96taty7hx4xgxYoR929NPP82CBQtcLrGt9NUutWvXpnbt2iUe9+qrr/L000/b1w8dOkSfPn147733aN++fXmG6DVKe6/AVuJx1VVX0bp1a+bOnavE429BQUG0bt2aFStW2JMPq9XKihUrGDlypLnBeSHDMBg1ahSffPIJq1atUuJRjF69evHrr786bBs2bBhNmjThkUceUeJxgc6dOxfotr1r1y7q169vUkTeKyMjo8DfcH9/f6xWq8vXrPTJR2nVq1fPYb1q1aoANGzYkIsvvtiMkLzWwYMH6dGjB/Xr12fatGkcPXrUvi8mJsbEyLzD2LFjGTJkCG3atKFdu3a8/PLLnDlzhmHDhpkdmtcZMWIECxcuZMmSJVSrVo2UlBQAIiIiCA0NNTk671KtWrUCbWGqVKlCzZo11UamEGPGjKFTp05MnjyZQYMGsWHDBubMmcOcOXPMDs3rJCQk8Mwzz1CvXj0uv/xyfvrpJ1588UXuuusu1y9alu43ldmePXvU1bYIc+fONYBCX2Izffp0o169ekZQUJDRrl07Y926dWaH5JWK+j2aO3eu2aH5BHW1Ld6yZcuMZs2aGcHBwUaTJk2MOXPmmB2SV0pPTzceeugho169ekZISIjRoEED47HHHjOysrJcvmalb/MhIiIinqWKeBEREfEoJR8iIiLiUUo+RERExKOUfIiIiIhHKfkQERERj1LyISIiIh6l5ENEREQ8SsmHiIiIeJSSDxEREfEoJR8iIiLiUUo+RERExKP+H5rG5EiIL44EAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.errorbar(data_x, data_y, sigma_y, sigma_x, fmt=\"o\", label=\"data\")\n", "x = np.linspace(data_x[0], data_x[-1], 200)\n", "par = np.array(m.values)\n", "plt.plot(x, f(x, par), label=\"fit\")\n", "plt.legend()\n", "\n", "# check fit quality\n", "chi2 = m.fval\n", "ndof = len(data_y) - 3\n", "plt.title(f\"$\\\\chi^2 / n_\\\\mathrm{{dof}} = {chi2:.2f} / {ndof} = {chi2/ndof:.2f}$\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We obtained a good fit." ] } ], "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": "python3", "version": "3.10.8 (main, Oct 13 2022, 09:48:40) [Clang 14.0.0 (clang-1400.0.29.102)]" }, "vscode": { "interpreter": { "hash": "bdbf20ff2e92a3ae3002db8b02bd1dd1b287e934c884beb29a73dced9dbd0fa3" } } }, "nbformat": 4, "nbformat_minor": 2 }