{ "cells": [ { "attachments": {}, "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)" ] }, { "attachments": {}, "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": 21, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:37.436843Z", "start_time": "2020-02-21T10:26:37.432080Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "JAX version 0.4.8\n", "numba version 0.57.0\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__}\")" ] }, { "attachments": {}, "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": 22, "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" ] }, { "attachments": {}, "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": 23, "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": 24, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:38.532308Z", "start_time": "2020-02-21T10:26:38.368563Z" } }, "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", "
Migrad
FCN = 496.2 Nfcn = 66
EDM = 1.84e-08 (Goal: 0.0001)
Valid Minimum Below EDM threshold (goal x 10)
No parameters at limit Below call limit
Hesse ok Covariance APPROXIMATE
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 66 │\n", "│ EDM = 1.84e-08 (Goal: 0.0001) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ Below EDM threshold (goal x 10) │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ No parameters at limit │ Below call limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Hesse ok │ Covariance APPROXIMATE │\n", "└──────────────────────────────────┴──────────────────────────────────────┘" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m1 = make_and_run_minuit(nll)\n", "m1.fmin" ] }, { "cell_type": "code", "execution_count": 25, "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", "
Migrad
FCN = 496.2 Nfcn = 26, Ngrad = 6
EDM = 1.84e-08 (Goal: 0.0001) time = 0.2 sec
Valid Minimum Below EDM threshold (goal x 10)
No parameters at limit Below call limit
Hesse ok Covariance APPROXIMATE
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 26, Ngrad = 6 │\n", "│ EDM = 1.84e-08 (Goal: 0.0001) │ time = 0.2 sec │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ Below EDM threshold (goal x 10) │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ No parameters at limit │ Below call limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Hesse ok │ Covariance APPROXIMATE │\n", "└──────────────────────────────────┴──────────────────────────────────────┘" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m2 = make_and_run_minuit(nll, grad=jax.grad(nll))\n", "m2.fmin" ] }, { "cell_type": "code", "execution_count": 26, "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", "
Migrad
FCN = 496.2 Nfcn = 66
EDM = 1.84e-08 (Goal: 0.0001)
Valid Minimum Below EDM threshold (goal x 10)
No parameters at limit Below call limit
Hesse ok Covariance APPROXIMATE
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 66 │\n", "│ EDM = 1.84e-08 (Goal: 0.0001) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ Below EDM threshold (goal x 10) │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ No parameters at limit │ Below call limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Hesse ok │ Covariance APPROXIMATE │\n", "└──────────────────────────────────┴──────────────────────────────────────┘" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m3 = make_and_run_minuit(jax.jit(nll))\n", "m3.fmin" ] }, { "cell_type": "code", "execution_count": 27, "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", "
Migrad
FCN = 496.2 Nfcn = 26, Ngrad = 6
EDM = 1.84e-08 (Goal: 0.0001) time = 0.2 sec
Valid Minimum Below EDM threshold (goal x 10)
No parameters at limit Below call limit
Hesse ok Covariance APPROXIMATE
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 26, Ngrad = 6 │\n", "│ EDM = 1.84e-08 (Goal: 0.0001) │ time = 0.2 sec │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ Below EDM threshold (goal x 10) │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ No parameters at limit │ Below call limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Hesse ok │ Covariance APPROXIMATE │\n", "└──────────────────────────────────┴──────────────────────────────────────┘" ] }, "execution_count": 27, "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": 28, "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", "
Migrad
FCN = 496.2 Nfcn = 82
EDM = 5.31e-05 (Goal: 0.0001) time = 1.0 sec
Valid Minimum Below EDM threshold (goal x 10)
No parameters at limit Below call limit
Hesse ok Covariance APPROXIMATE
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 496.2 │ Nfcn = 82 │\n", "│ EDM = 5.31e-05 (Goal: 0.0001) │ time = 1.0 sec │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ Below EDM threshold (goal x 10) │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ No parameters at limit │ Below call limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Hesse ok │ Covariance APPROXIMATE │\n", "└──────────────────────────────────┴──────────────────────────────────────┘" ] }, "execution_count": 28, "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": 29, "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": 30, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:26:45.142272Z", "start_time": "2020-02-21T10:26:45.033451Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqAAAAGwCAYAAAB2AtfDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHPUlEQVR4nO3de3zP9f//8ft7m53ZHLepMYqZw5RjkzKZ5pDQp/Jx3Cqn5IMcwiehnEPxoRSyEX1EoYMiZKphOUfmfBh9HHLaHMphe/7+6Of99c7MxvZ6Y7fr5fK+5PV8P1/P1+P13Nu7u9dpNmOMEQAAAGARF2cXAAAAgPyFAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoAizxhjlJaWJh41CwAArkUARZ45e/as/Pz8dPbsWWeXAgAA7iAEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApAigAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigAIAAMBSbs4uAPe+ykOWysXD29ll3JUOjG7q7BIAAMh1HAEFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARS4C/zwww9q1qyZSpYsKZvNpkWLFt2wb9euXWWz2TRhwoQsxxw1apRq1qypggULqkSJEmrRooV27tyZu4UDAJAJAihwFzh//ryqVq2q9957L8t+Cxcu1Nq1a1WyZMmbjrlq1Sq98sorWrt2rZYtW6bLly/rySef1Pnz53OrbAAAMkUABe4CjRs31vDhw9WyZcsb9vntt9/0r3/9S3PmzFGBAgVuOuaSJUsUGxurSpUqqWrVqoqPj1dKSoo2bNggSUpISJC7u7t+/PFH+zpvv/22SpQooWPHjt3+TgEA8i0CKHAPyMjIUPv27dWvXz9VqlTplsZITU2VJBUpUkSSFBkZqV69eql9+/ZKTU3Vpk2b9MYbb2j69OkKCAjItdoBAPkPARTZdrNrD+E8Y8aMkZubm3r06HFL62dkZKhXr1569NFHVblyZXv78OHDVbhwYXXu3Fnt2rVTTEyMnn766dwqGwCQT+X7APr3UHV1OT4+XjabLcvXgQMHnFY3cNWGDRs0ceJE+2f2Vrzyyivatm2b5s6d69Du7u6uOXPm6PPPP9eff/6pd999NzdKBgDkc/k+gN5Iq1atdOTIEfsrIiJCnTp1cmgLDg52dpk5kp6eroyMDGeXgVz2448/6vjx4ypVqpTc3Nzk5uamgwcPqk+fPgoJCbnp+t27d9fXX3+tlStX6v7777/u/dWrV0uSTp06pVOnTuV2+QCAfOiuCaCRkZHq0aOHXnvtNRUpUkSBgYEaOnSoQ5+UlBQ1b95cvr6+KlSokJ5//vlbvlnCy8tLgYGB9pe7u7u8vb0d2lxdXbM1VkhIiEaOHKkXX3xRBQsWVKlSpTR16lSHPlu3btUTTzwhLy8vFS1aVJ07d9a5c+eyHPfLL79UuXLl5Onpqfr162vmzJmy2Ww6c+aMJCk+Pl7+/v768ssvVbFiRXl4eCglJUXr1q1Tw4YNVaxYMfn5+alevXrauHGjw9i7d+/W448/Lk9PT1WsWFHLli3L/uTBUu3bt9cvv/yizZs3218lS5ZUv379tHTp0huuZ4xR9+7dtXDhQn3//fcqU6bMdX327t2rV199VdOmTVPt2rUVExPDP2IAALftrgmgkjRz5kz5+PgoKSlJb7/9tt566y17MMrIyFDz5s116tQprVq1SsuWLdO+ffvUqlUrJ1f9l/Hjx6tGjRratGmTunXrppdfftn+zMXz588rOjpahQsX1rp16zR//nwtX75c3bt3v+F4+/fv17PPPqsWLVpoy5Yt6tKli15//fXr+l24cEFjxozR9OnT9euvv6pEiRI6e/asYmJi9NNPP2nt2rUqV66cmjRporNnz0r6ay6feeYZubu7KykpSR988IH69+9/0328ePGi0tLSHF7IHefOnbOHS+mvn//mzZuVkpKiokWLqnLlyg6vAgUKKDAwUKGhofYxGjRooMmTJ9uXX3nlFc2ePVuffPKJChYsqKNHj+ro0aP6448/JP11xLxdu3aKjo7WCy+8oLi4OP3yyy8aP368pfsOALj3uDm7gJwIDw/XkCFDJEnlypXT5MmTtWLFCjVs2FArVqzQ1q1btX//fvup8VmzZqlSpUpat26datas6czS1aRJE3Xr1k2S1L9/f7377rtauXKlQkND9cknn+jPP//UrFmz5OPjI0maPHmymjVrpjFjxmR6x/GHH36o0NBQjR07VpIUGhqqbdu2acSIEQ79Ll++rPfff19Vq1a1tz3xxBMOfaZOnSp/f3+tWrVKTz31lJYvX64dO3Zo6dKl9udJjhw5Uo0bN85yH0eNGqU333wzhzOD7Fi/fr3q169vX+7du7ckKSYmRvHx8dkaY+/evTpx4oR9ecqUKZL+Ortwrbi4OMXGxmrEiBE6ePCgvv76a0lSUFCQpk6dqtatW+vJJ590+EwBAJATd10AvVZQUJCOHz8uSUpOTlZwcLDDdZkVK1aUv7+/kpOTnR5Ar63dZrMpMDDQofaqVavaw6ckPfroo8rIyNDOnTszDaA7d+68bp9q1ap1XT93d/fr5u3YsWMaNGiQEhISdPz4caWnp+vChQtKSUmx1xMcHOzwMPOIiIib7uPAgQPtwUiS0tLS7rrrZO9UkZGRMsZku39mN8j9ve1m4w0ePFiDBw92aHvmmWd08eLFbNcBAEBm7qoA+veHa9tstrvmejRn1e7l5XXdndExMTE6efKkJk6cqNKlS8vDw0MRERG6dOnSbW3Lw8NDHh4etzUGAAC4991V14BmJSwsTIcOHdKhQ4fsbdu3b9eZM2dUsWJFJ1Z2c2FhYdqyZYvDr0BMTEyUi4uLwzV81woNDdX69esd2tatW5et7SUmJqpHjx5q0qSJKlWqJA8PD4dTs1fn8siRI/a2tWvX5mSXAAAAbuieCaBRUVGqUqWK2rZtq40bN+rnn39Whw4dVK9ePdWoUcPZ5WWpbdu28vT0VExMjLZt26aVK1fqX//6l9q3b3/D3zjTpUsX7dixQ/3799euXbs0b948+7WAN3sWZLly5fTxxx8rOTlZSUlJatu2rby8vOzvR0VFqXz58oqJidGWLVv0448/ZnqDEwAAwK24ZwKozWbTF198ocKFC+vxxx9XVFSUypYtq08//fSG61w9Be7mdutXIiQkJNz2Q+m9vb21dOlSnTp1SjVr1tSzzz573R3Lf1emTBl99tlnWrBggcLDwzVlyhR7SLzZafCPPvpIp0+fVrVq1dS+fXv16NFDJUqUsL/v4uKihQsX6o8//lCtWrXUsWPH625uAgAAuFU2k5M7G+4xR48eVVBQkNatW3fLR0nj4uI0cuRIbd++/brrPK02YsQIffDBBw6XIThTWlqa/Pz8FNxrnlw8vJ1dzl3pwOimzi4BAIBcd1fdhJRbjDE6ePCgxo0bp4CAAIfffZ1T33zzjUaOHOmU8Pn++++rZs2aKlq0qBITEzV27Ngsnx0KAABwJ8iXATQ1NVWhoaEKCwvT3Llz5enpectjzZ8/Pxcry5ndu3dr+PDhOnXqlEqVKqU+ffpo4MCBTqsHAAAgO/L1KXjkLU7B3z5OwQMA7kX3zE1IAAAAuDsQQAEAAGApAigAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKZsxxji7CNyb0tLS5Ofnp9TUVBUqVMjZ5QAAgDsER0ABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApAigAAAAsRQAFAACApQigAAAAsBQBFAAAAJZyc3YBuPdVHrJULh7ezi4DAIB7xoHRTZ1dwm3hCCgAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAAHeh9PR0vfHGGypTpoy8vLz0wAMPaNiwYTLG2PsMHTpUFSpUkI+PjwoXLqyoqCglJSXddOz33ntPISEh8vT0VO3atfXzzz/nau0EUAAAgLvQmDFjNGXKFE2ePFnJyckaM2aM3n77bU2aNMnep3z58po8ebK2bt2qn376SSEhIXryySf1+++/33DcTz/9VL1799aQIUO0ceNGVa1aVdHR0Tp+/Hiu1U4ABQAAuAutXr1azZs3V9OmTRUSEqJnn31WTz75pMPRyjZt2igqKkply5ZVpUqV9M477ygtLU2//PLLDcd955131KlTJ73wwguqWLGiPvjgA3l7e2vGjBmSpISEBLm7u+vHH3+0r/P222+rRIkSOnbsWLZqJ4ACAADcherUqaMVK1Zo165dkqQtW7bop59+UuPGjTPtf+nSJU2dOlV+fn6qWrXqDfts2LBBUVFR9jYXFxdFRUVpzZo1kqTIyEj16tVL7du3V2pqqjZt2qQ33nhD06dPV0BAQLZqd8vJjiJ/i42N1ZkzZ7Ro0SJnlwIAQL43YMAApaWlqUKFCnJ1dVV6erpGjBihtm3bOvT7+uuv9c9//lMXLlxQUFCQli1bpmLFimU65okTJ5Senn5dkAwICNCOHTvsy8OHD9eyZcvUuXNnbdu2TTExMXr66aezXTtHQPOIzWZzCGpXl+Pj42Wz2bJ8HThwwGl1AwCAu8O8efM0Z84cffLJJ9q4caNmzpypcePGaebMmQ796tevr82bN2v16tVq1KiRnn/++du+ntPd3V1z5szR559/rj///FPvvvtujtYngFqsVatWOnLkiP0VERGhTp06ObQFBwfn2fYvX76cZ2MDAADr9OvXTwMGDNA///lPValSRe3bt9err76qUaNGOfTz8fHRgw8+qEceeUQfffSR3Nzc9NFHH2U6ZrFixeTq6nrdtZzHjh1TYGCgQ9vq1aslSadOndKpU6dyVHu+D6CRkZHq0aOHXnvtNRUpUkSBgYEaOnSoQ5+UlBQ1b95cvr6+KlSokJ5//vlsX2T7d15eXgoMDLS/3N3d5e3t7dDm6uqarbGOHDmipk2bysvLS2XKlNEnn3yikJAQTZgwwd7HZrNpypQpevrpp+Xj46MRI0YoPT1dL730kv2xDaGhoZo4caLD2Onp6erdu7f8/f1VtGhRvfbaaw6PdQAAAM514cIFubg4RjlXV1dlZGRkuV5GRoYuXryY6Xvu7u6qXr26VqxY4dB/xYoVioiIsLft3btXr776qqZNm6batWsrJibmptu9Vr4PoJI0c+ZM+fj4KCkpSW+//bbeeustLVu2TNJfk968eXOdOnVKq1at0rJly7Rv3z61atXKyVVLHTp00P/+9z8lJCTo888/19SpUzM9pD506FC1bNlSW7du1YsvvqiMjAzdf//9mj9/vrZv367Bgwfr3//+t+bNm2dfZ/z48YqPj9eMGTP0008/6dSpU1q4cGGW9Vy8eFFpaWkOLwAAkDeaNWumESNGaPHixTpw4IAWLlyod955Ry1btpQknT9/Xv/+97+1du1aHTx4UBs2bNCLL76o3377Tc8995x9nAYNGmjy5Mn25d69e2vatGmaOXOmkpOT9fLLL+v8+fN64YUXJP11kKpdu3aKjo7WCy+8oLi4OP3yyy8aP358tmvnJiRJ4eHhGjJkiCSpXLlymjx5slasWKGGDRtqxYoV2rp1q/bv328/NT5r1ixVqlRJ69atU82aNZ1S844dO7R8+XKtW7dONWrUkCRNnz5d5cqVu65vmzZt7B+aq9588037n8uUKaM1a9Zo3rx5ev755yVJEyZM0MCBA/XMM89Ikj744AMtXbo0y5pGjRrlMC4AAMg7kyZN0htvvKFu3brp+PHjKlmypLp06aLBgwdL+uto6I4dOzRz5kydOHFCRYsWVc2aNfXjjz+qUqVK9nH27t2rEydO2JdbtWql33//XYMHD9bRo0f10EMPacmSJfYbk0aMGKGDBw/q66+/liQFBQVp6tSpat26tZ588skb3mF/LQKo/gqg1woKCrIfSUxOTlZwcLDDdZkVK1aUv7+/kpOTnRZAd+7cKTc3N1WrVs3e9uCDD6pw4cLX9b0aUK/13nvvacaMGUpJSdEff/yhS5cu6aGHHpIkpaam6siRI6pdu7a9v5ubm2rUqJHlafiBAweqd+/e9uW0tLQ8vZ4VAID8rGDBgpowYYLDpXfX8vT01IIFC246TmY3P3fv3l3du3fPtP/gwYPtIfeqZ5555oan9TNDAJVUoEABh2WbzZaj6xjudD4+Pg7Lc+fOVd++fTV+/HhFRESoYMGCGjt2bLZ+NVdWPDw85OHhcVtjAACAex/XgN5EWFiYDh06pEOHDtnbtm/frjNnzqhixYpOqys0NFRXrlzRpk2b7G179uzR6dOnb7puYmKi6tSpo27duunhhx/Wgw8+qL1799rf9/PzU1BQkEMgvXLlijZs2JC7OwEAAPIlAuhNREVFqUqVKmrbtq02btyon3/+WR06dFC9evUyPbVtlQoVKigqKkqdO3fWzz//rE2bNqlz587y8vKSzWbLct1y5cpp/fr1Wrp0qXbt2qU33nhD69atc+jTs2dPjR49WosWLdKOHTvUrVs3nTlzJg/3CAAA5BcE0Juw2Wz64osvVLhwYT3++OP236f66aef3nCdq6fv3dxu/QqHhISEmz6UftasWQoICNDjjz+uli1bqlOnTipYsKA8PT2zHLtLly565pln1KpVK9WuXVsnT55Ut27dHPr06dNH7du3V0xMjP00/dW76gAAAG6HzfBwx1x39OhRBQUFOdyhnlNxcXEaOXKktm/fft01qjdy+PBhBQcHa/ny5WrQoMEtbTc3paWlyc/PT8G95snFw9vZ5QAAcM84MLqps0u4LdyElIuMMTp48KDGjRungIAAVa5c+ZbH+uabbzRy5Mgsw+f333+vc+fOqUqVKjpy5Ihee+01hYSE6PHHH7/l7QIAAOQ1AmguSk1NVWhoqMLCwjR37tybngrPyvz582/a5/Lly/r3v/+tffv2qWDBgqpTp47mzJmT7SOmAAAAzkAAzUX+/v45egbW7YqOjlZ0dLRl2wMAAMgN3IQEAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApmzHGOLsI3JvS0tLk5+en1NRUFSpUyNnlAACAOwRHQAEAAGApAigAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAlnJzdgG491UeslQuHt63vP6B0U1zsRoAAOBsHAEFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARR3hSlTpig8PFyFChVSoUKFFBERoW+//TbLdebPn68KFSrI09NTVapU0TfffGNRtQAAICsEUNwV7r//fo0ePVobNmzQ+vXr9cQTT6h58+b69ddfM+2/evVqtW7dWi+99JI2bdqkFi1aqEWLFtq2bZvFlQMAgL8jgOKu0KxZMzVp0kTlypVT+fLlNWLECPn6+mrt2rWZ9p84caIaNWqkfv36KSwsTMOGDVO1atU0efJkSdKOHTvk7e2tTz75xL7OvHnz5OXlpe3bt1uyTwAA5FcEUNx10tPTNXfuXJ0/f14RERGZ9lmzZo2ioqIc2qKjo7VmzRpJUoUKFTRu3Dh169ZNKSkpOnz4sLp27aoxY8aoYsWKeb4PAADkZ3d0AI2NjVWLFi2cXQYkxcfHy9/f36k1bN26Vb6+vvLw8FDXrl21cOHCG4bFo0ePKiAgwKEtICBAR48etS9369ZNdevWVbt27RQbG6uaNWvqX//6V57uAwAAuMMD6MSJExUfH59n4/89VF27HBkZKZvNdsNXZGRkntWFzIWGhmrz5s1KSkrSyy+/rJiYmNs+XT5jxgz98ssv2rhxo+Lj42Wz2XKpWgAAcCNuzi4gK35+fk7b9oIFC3Tp0iVJ0qFDh1SrVi0tX75clSpVkiS5u7s7rbZbdenSpbuy7qvc3d314IMPSpKqV6+udevWaeLEifrwww+v6xsYGKhjx445tB07dkyBgYEObVu2bNH58+fl4uKiI0eOKCgoKO92AAAASLrDj4Beewp+yZIlqlu3rvz9/VW0aFE99dRT2rt3r73vrFmz5Ovrq927d9vbunXrpgoVKujChQs53naRIkUUGBiowMBAFS9eXJJUtGhRe1uRIkWyNU5CQoJsNptWrFihGjVqyNvbW3Xq1NHOnTsd+k2ZMkUPPPCA3N3dFRoaqo8//jjLca9cuaIePXrY56N///6KiYlxuGQhMjJS3bt3V69evVSsWDFFR0dLkt555x1VqVJFPj4+Cg4OVrdu3XTu3DmH8ePj41WqVCl5e3urZcuWOnnyZLb210oZGRm6ePFipu9FRERoxYoVDm3Lli1zuGb01KlTio2N1euvv67Y2Fi1bdtWf/zxR57WDAAA7vAAeq3z58+rd+/eWr9+vVasWCEXFxe1bNlSGRkZkqQOHTqoSZMmatu2ra5cuaLFixdr+vTpmjNnjry9vZ1cvfT6669r/PjxWr9+vdzc3PTiiy/a31u4cKF69uypPn36aNu2berSpYteeOEFrVy58objjRkzRnPmzFFcXJwSExOVlpamRYsWXddv5syZcnd3V2Jioj744ANJkouLi/7zn//o119/1cyZM/X999/rtddes6+TlJSkl156Sd27d9fmzZtVv359DR8+/Kb7ePHiRaWlpTm8csvAgQP1ww8/6MCBA9q6dasGDhyohIQEtW3bVtJfP/+BAwfa+/fs2VNLlizR+PHjtWPHDg0dOlTr169X9+7d7X26du2q4OBgDRo0SO+8847S09PVt2/fXKsZAABk7o4+BX+tf/zjHw7LM2bMUPHixbV9+3ZVrlxZkvThhx8qPDxcPXr00IIFCzR06FBVr17dGeVeZ8SIEapXr54kacCAAWratKn+/PNPeXp6aty4cYqNjVW3bt0kSb1799batWs1btw41a9fP9PxJk2apIEDB6ply5aSpMmTJ2f6oPVy5crp7bffdmjr1auX/c8hISEaPny4unbtqvfff1/S/z3C6GooLV++vFavXq0lS5ZkuY+jRo3Sm2++mY3ZyLnjx4+rQ4cOOnLkiPz8/BQeHq6lS5eqYcOGkqSUlBS5uPzfv6fq1KmjTz75RIMGDdK///1vlStXTosWLbJ/VmbNmqVvvvlGmzZtkpubm9zc3DR79mzVrVtXTz31lBo3bpwn+wEAAO6iALp7924NHjxYSUlJOnHihP3IZ0pKij1UFC5cWB999JGio6NVp04dDRgwwJklOwgPD7f/+ep1hsePH1epUqWUnJyszp07O/R/9NFHNXHixEzHSk1N1bFjx1SrVi17m6urq6pXr26fl6syC+DLly/XqFGjtGPHDqWlpenKlSv6888/deHCBXl7eys5OdkebK+KiIi4aQAdOHCgevfubV9OS0tTcHBwlutk10cffZTl+wkJCde1Pffcc3ruuecy7d+hQwd16NDBoa1WrVr2634BAEDeuWtOwTdr1kynTp3StGnTlJSUpKSkJEm6LjD88MMPcnV11ZEjR3T+/HlnlJqpAgUK2P989U7rv4fFvODj4+OwfODAAT311FMKDw/X559/rg0bNui9996TdP1c5pSHh4f9V2VefQEAAPzdXRFAT548qZ07d2rQoEFq0KCBwsLCdPr06ev6rV69WmPGjNFXX30lX19fh+v97mRhYWFKTEx0aEtMTLzhMy79/PwUEBCgdevW2dvS09O1cePGm25rw4YNysjI0Pjx4/XII4+ofPny+t///nddPVcD/lU3+o1DAAAAOXVXnIIvXLiwihYtqqlTpyooKEgpKSnXnV4/e/as2rdvrx49eqhx48a6//77VbNmTTVr1kzPPvuskyrPnn79+un555/Xww8/rKioKH311VdasGCBli9ffsN1/vWvf2nUqFF68MEHVaFCBU2aNEmnT5++6XMsH3zwQV2+fFmTJk1Ss2bNHG5OuqpHjx569NFHNW7cODVv3lxLly696el3AACA7LorjoC6uLho7ty52rBhgypXrqxXX31VY8eOdejTs2dP+fj4aOTIkZKkKlWqaOTIkerSpYt+++23TMfNyMiQm9vtZfDY2Njbfih9ixYtNHHiRI0bN06VKlXShx9+qLi4uCzH7d+/v1q3bq0OHTooIiJCvr6+io6OlqenZ5bbqlq1qt555x2NGTNGlStX1pw5czRq1CiHPo888oimTZumiRMnqmrVqvruu+80aNCg29pHAACAq2zGGOPsIm6kdevWcnV11ezZs/Nk/NGjR2v27Nnatm3bLY9Rr1491a9fX0OHDs29wm5BRkaGwsLC9Pzzz2vYsGFOreWqtLQ0+fn5KbjXPLl43PqjsA6MbpqLVQEAAGe7I0/BX7lyRbt27dKaNWvUpUuXXB//woUL2rFjh+Li4m7rcTupqanau3evFi9enIvVZc/Bgwf13XffqV69erp48aImT56s/fv3q02bNpbXAgAAkBN35Cn4bdu2qUaNGqpUqZK6du2a6+NPnTpVUVFRqlq1qgYPHnzL4/j5+enw4cPy9fXNxeqyx8XFRfHx8apZs6YeffRRbd26VcuXL1dYWJjltQAAAOTEHX0KHnc3TsEDAIDM3JFHQAEAAHDvIoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApAigAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFI2Y4xxdhG4N6WlpcnPz0+pqakqVKiQs8sBAAB3CI6AAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApAigAAAAs5ebsAnDvqzxkqVw8vDN978DophZXAwAAnI0joAAAALAUARQAAACWIoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApAigAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAwul++OEHNWvWTCVLlpTNZtOiRYuyvW5iYqLc3Nz00EMP5Vl9AAAgdxFA4XTnz59X1apV9d577+VovTNnzqhDhw5q0KBBHlUGAADyAgEUTte4cWMNHz5cLVu2zNF6Xbt2VZs2bRQREeHQ/vvvvyswMFAjR460t61evVru7u5asWJFrtQMAABuHQEUd6W4uDjt27dPQ4YMue694sWLa8aMGRo6dKjWr1+vs2fPqn379urevTtHSwEAuAPkywAaGxurFi1aOLuMu05Or8/MK7t379aAAQM0e/Zsubm5ZdqnSZMm6tSpk9q2bauuXbvKx8dHo0aNsrhSAACQmXwZQCdOnKj4+Pg8Gz8+Pl7+/v6ZLkdGRspms93wFRkZmWd13QvS09PVpk0bvfnmmypfvnyWfceNG6crV65o/vz5mjNnjjw8PCyqEgAAZCXzw0f3OD8/P6dte8GCBbp06ZIk6dChQ6pVq5aWL1+uSpUqSZLc3d3zbNvp6emy2Wxycbl7/91x9uxZrV+/Xps2bVL37t0lSRkZGTLGyM3NTd99952eeOIJSdLevXv1v//9TxkZGTpw4ICqVKnizNIBAMD/d/cmkdtw7Sn4JUuWqG7duvL391fRokX11FNPae/evfa+s2bNkq+vr3bv3m1v69atmypUqKALFy7keNtFihRRYGCgAgMDVbx4cUlS0aJF7W1FihTJ9lhffvmlypUrJ09PT9WvX18zZ86UzWbTmTNnJP3fkdcvv/xSFStWlIeHh1JSUrRu3To1bNhQxYoVk5+fn+rVq6eNGzc6jL179249/vjj8vT0VMWKFbVs2bIc72teKFSokLZu3arNmzfbX127dlVoaKg2b96s2rVrS5IuXbqkdu3aqVWrVho2bJg6duyo48ePO7l6AAAg5dMjoNc6f/68evfurfDwcJ07d06DBw9Wy5YttXnzZrm4uKhDhw76+uuv1bZtW61evVpLly7V9OnTtWbNGnl7ezut7v379+vZZ59Vz5491bFjR23atEl9+/a9rt+FCxc0ZswYTZ8+XUWLFlWJEiW0b98+xcTEaNKkSTLGaPz48WrSpIl2796tggULKiMjQ88884wCAgKUlJSk1NRU9erV66Y1Xbx4URcvXrQvp6WlZWtfzp07pz179jjs2+bNm1WkSBGVKlVKAwcO1G+//aZZs2bJxcVFlStXdli/RIkS8vT0dGh//fXXlZqaqv/85z/y9fXVN998oxdffFFff/11tmoCAAB5J98H0H/84x8OyzNmzFDx4sW1fft2e6D58MMPFR4erh49emjBggUaOnSoqlev7oxy7T788EOFhoZq7NixkqTQ0FBt27ZNI0aMcOh3+fJlvf/++6pataq97eop6qumTp0qf39/rVq1Sk899ZSWL1+uHTt2aOnSpSpZsqQkaeTIkWrcuHGWNY0aNUpvvvlmjvdl/fr1ql+/vn25d+/ekqSYmBjFx8fryJEjSklJyfZ4CQkJmjBhglauXKlChQpJkj7++GNVrVpVU6ZM0csvv5zjGgEAQO7J9wF09+7dGjx4sJKSknTixAllZGRIklJSUuwBtHDhwvroo48UHR2tOnXqaMCAAc4sWZK0c+dO1axZ06GtVq1a1/Vzd3dXeHi4Q9uxY8c0aNAgJSQk6Pjx40pPT9eFCxfsIS85OVnBwcH28CnpumdtZmbgwIH28Cj9dQQ0ODj4putFRkbKGHPD9292w9jQoUM1dOhQh/EuX77s0CckJESpqak3rQUAAOS9fB9AmzVrptKlS2vatGkqWbKkMjIyVLlyZfuNQlf98MMPcnV11ZEjR3T+/HkVLFjQSRXnjJeXl2w2m0NbTEyMTp48qYkTJ6p06dLy8PBQRETEdfucUx4eHtxpDgAAbipf3oR01cmTJ7Vz504NGjRIDRo0UFhYmE6fPn1dv9WrV2vMmDH66quv5Ovra7/72plCQ0O1fv16h7Z169Zla93ExET16NFDTZo0UaVKleTh4aETJ07Y3w8LC9OhQ4d05MgRe9vatWtzp3AAAJDv5esAWrhwYRUtWlRTp07Vnj179P333zucQpZk/y06PXr0UOPGjTVnzhx9+umn+uyzz5xU9V+6dOmiHTt2qH///tq1a5fmzZtnP1X99yOef1euXDl9/PHHSk5OVlJSktq2bSsvLy/7+1FRUSpfvrxiYmK0ZcsW/fjjj3r99dfzcncAAEA+kq8DqIuLi+bOnasNGzaocuXKevXVV+039VzVs2dP+fj42H+veJUqVTRy5Eh16dJFv/32W6bjZmRk3PA39GRXbGxslg+lL1OmjD777DMtWLBA4eHhmjJlij0k3uw0+EcffaTTp0+rWrVq9nBdokQJ+/suLi5auHCh/vjjD9WqVUsdO3a87uYmAACAW2UzWd39cY9q3bq1XF1dNXv27DwZf/To0Zo9e7a2bdt2y2PUq1dP9evXd7i55mZGjBihDz74QIcOHbrl7eamtLQ0+fn5KbjXPLl4ZP7IqgOjm1pcFQAAcLZ8dRPSlStXtGvXLq1Zs0ZdunTJ9fEvXLigHTt2KC4u7qaPLMpKamqq9u7dq8WLF2fZ7/3331fNmjVVtGhRJSYmauzYsXfE9akAAABZyVen4Ldt26YaNWqoUqVK6tq1a66PP3XqVEVFRalq1aoaPHjwLY/j5+enw4cPy9fXN8t+u3fvVvPmzVWxYkUNGzZMffr0ydERUwAAAGfIl6fgYQ1OwQMAgMzkqyOgAAAAcD4CKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApAigAAAAsZTPGGGcXgXtTWlqa/Pz8lJqaqkKFCjm7HAAAcIfgCCgAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFJuzi4A977KQ5bKxcPbvnxgdFMnVgMAAJyNI6AAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApAigAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigMJpfvjhBzVr1kwlS5aUzWbTokWLbrpOQkKCqlWrJg8PDz344IOKj4/P8zoBAEDuIoDCac6fP6+qVavqvffey1b//fv3q2nTpqpfv742b96sXr16qWPHjlq6dGkeVwoAAHITARRO07hxYw0fPlwtW7bMVv8PPvhAZcqU0fjx4xUWFqbu3bvr2Wef1bvvvitJ+v333xUYGKiRI0fa11m9erXc3d21YsWKPNkHAACQcwRQ3DXWrFmjqKgoh7bo6GitWbNGklS8eHHNmDFDQ4cO1fr163X27Fm1b99e3bt3V4MGDZxRMgAAyMRdE0AjIyPVq1cvZ5eRrzn7Z3D06FEFBAQ4tAUEBCgtLU1//PGHJKlJkybq1KmT2rZtq65du8rHx0ejRo1yRrkAAOAG7poAmlf+HqquLh84cEA2my3LFzfA3JnGjRunK1euaP78+ZozZ448PDycXRIAALiGm7MLuFMFBwfryJEj9uVx48ZpyZIlWr58ub3Nz8/PGaXdlsuXL6tAgQLOLuOWBAYG6tixYw5tx44dU6FCheTl5WVv27t3r/73v/8pIyNDBw4cUJUqVawuFQAAZCFHR0AjIyPVo0cPvfbaaypSpIgCAwM1dOhQ+/tXjxpu3rzZ3nbmzBnZbDYlJCRI+usxOjabTUuXLtXDDz8sLy8vPfHEEzp+/Li+/fZbhYWFqVChQmrTpo0uXLjgsP0rV66oe/fu8vPzU7FixfTGG2/IGGN//+OPP1aNGjVUsGBBBQYGqk2bNjp+/HjOZ0WSq6urAgMD7S9fX1+5ubk5tF0berISGxurFi1aaNy4cQoKClLRokX1yiuv6PLly/Y+p0+fVocOHVS4cGF5e3urcePG2r17d5bj7tixQ3Xr1pWnp6cqVqyo5cuXOzzO6OrP49NPP1W9evXk6empOXPm6OTJk2rdurXuu+8+eXt7q0qVKvrvf//rMPb58+fVoUMH+fr6KigoSOPHj8/ZBOaBiIiI624mWrZsmSIiIuzLly5dUrt27dSqVSsNGzZMHTt2vOXPAAAAyBs5PgU/c+ZM+fj4KCkpSW+//bbeeustLVu2LMcbHjp0qCZPnqzVq1fr0KFDev755zVhwgR98sknWrx4sb777jtNmjTpum27ubnp559/1sSJE/XOO+9o+vTp9vcvX76sYcOGacuWLVq0aJEOHDig2NjYHNeWF1auXKm9e/dq5cqVmjlzpuLj4x1O4cfGxmr9+vX68ssvtWbNGhlj1KRJE4eQeq309HS1aNFC3t7eSkpK0tSpU/X6669n2nfAgAHq2bOnkpOTFR0drT///FPVq1fX4sWLtW3bNnXu3Fnt27fXzz//bF+nX79+WrVqlb744gt99913SkhI0MaNG7Pcx4sXLyotLc3hlZVz585p8+bN9n+w7N+/X5s3b1ZKSookaeDAgerQoYO9f9euXbVv3z699tpr2rFjh95//33NmzdPr776qr3P66+/rtTUVP3nP/9R//79Vb58eb344otZ1gEAACxmcqBevXqmbt26Dm01a9Y0/fv3N8YYs3//fiPJbNq0yf7+6dOnjSSzcuVKY4wxK1euNJLM8uXL7X1GjRplJJm9e/fa27p06WKio6Mdth0WFmYyMjLsbf379zdhYWE3rHfdunVGkjl79myW+9SzZ88bLl81ZMgQU7Vq1RuOk5WYmBhTunRpc+XKFXvbc889Z1q1amWMMWbXrl1GkklMTLS/f+LECePl5WXmzZuX6ZjffvutcXNzM0eOHLG3LVu2zEgyCxcuNMb8389jwoQJN62xadOmpk+fPsYYY86ePWvc3d0dtn3y5Enj5eWV6dxcNWTIECPpuldwr3mmdP+v7a+rrn4W/v6KiYmxz1u9evUctrFy5Urz0EMPGXd3d1O2bFkTFxfn8J6bm5v58ccf7W379+83hQoVMu+///5N5wAAAFgjx9eAhoeHOywHBQXd0inOa8cJCAiQt7e3ypYt69B27RE5SXrkkUdks9nsyxERERo/frzS09Pl6uqqDRs2aOjQodqyZYtOnz6tjIwMSVJKSooqVqyY4xpzU6VKleTq6mpfDgoK0tatWyVJycnJcnNzU+3ate3vFy1aVKGhoUpOTs50vJ07dyo4OFiBgYH2tlq1amXat0aNGg7L6enpGjlypObNm6fffvtNly5d0sWLF+Xt7S3pr2soL1265FBPkSJFFBoamuU+Dhw4UL1797Yvp6WlKTg4+Ib9IyMjHS6h+LvMbvKKjIzUpk2bbjje348Yh4SEKDU1Ncu6AQCAtXIcQP9+A4vNZrMHPReXv87oXxsqbnQK+dpxbDZbluNmx/nz5xUdHa3o6GjNmTNHxYsXV0pKiqKjo3Xp0qVsj5NXbnf/boePj4/D8tixYzVx4kRNmDBBVapUkY+Pj3r16nXb8+Th4cEd5wAA4KZy9TFMxYsXlySHu8evvSHpdiUlJTksr127VuXKlZOrq6t27NihkydPavTo0XrsscdUoUKFu+bmk7CwMF25csVh/06ePKmdO3fe8MhtaGioDh065HBX+Lp167K1vcTERDVv3lzt2rVT1apVVbZsWe3atcv+/gMPPKACBQo41HP69GmHPgAAALcqVwOol5eXHnnkEY0ePVrJyclatWqVBg0alGvjp6SkqHfv3tq5c6f++9//atKkSerZs6ckqVSpUnJ3d9ekSZO0b98+ffnllxo2bFiubTsvlStXTs2bN1enTp30008/acuWLWrXrp3uu+8+NW/ePNN1GjZsqAceeEAxMTH65ZdflJiYaJ/ray9TuNH2li1bptWrVys5OVldunRxCLK+vr566aWX1K9fP33//ffatm2bYmNj7Ue4AQAAbkeuJ4oZM2boypUrql69unr16qXhw4fn2tgdOnTQH3/8oVq1aumVV15Rz5491blzZ0l/HX2Nj4/X/PnzVbFiRY0ePVrjxo276ZgZGRlyc7u9x6HmxkPp4+LiVL16dT311FOKiIiQMUbffPPNDZ/Z6erqqkWLFuncuXOqWbOmOnbsaL8L3tPTM8ttDRo0SNWqVVN0dLQiIyMVGBioFi1aOPQZO3asHnvsMTVr1kxRUVGqW7euqlevflv7CAAAIEk2k9VdIPlAhQoV1LFjR/Xt2/eW1t+/f7/Kly+v7du3q1y5crlcXc4kJiaqbt262rNnjx544AGn1iL9dROSn5+fgnvNk4uHt739wOimTqwKAAA4W779TUhXH3y/c+dONWjQ4JbH+eabb9S5c2enhM+FCxfK19dX5cqV0549e9SzZ089+uijd0T4BAAAuJF8G0AbNWqk06dP6z//+Y8efvjhWx7nlVdeycWqcubs2bPq37+/UlJSVKxYMUVFRd0Rv7EIAAAgK/n+FDzyDqfgAQBAZritGQAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoACAADAUjZjjHF2Ebg3paWlyc/PT6mpqSpUqJCzywEAAHcIjoACAADAUgRQAAAAWIoACgAAAEsRQAEAAGApAigAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigAIAAMBSBFAAAABYigAKAAAASxFAAQAAYCkCKAAAACxFAAUAAIClCKAAAACwlJuzC8C9yxgjSUpLS3NyJQAAIKcKFiwom82WJ2MTQJFnTp48KUkKDg52ciUAACCnUlNTVahQoTwZmwCKPFOkSBFJUkpKivz8/JxczZ0lLS1NwcHBOnToUJ795b6bMT83xtxkjfm5MeYma8zP9QoWLJhnYxNAkWdcXP66xNjPz4+/zDdQqFAh5iYLzM+NMTdZY35ujLnJGvNjDW5CAgAAgKUIoAAAALAUARR5xsPDQ0OGDJGHh4ezS7njMDdZY35ujLnJGvNzY8xN1pgfa9nM1WflAAAAABbgCCgAAAAsRQAFAACApQigAAAAsBQBFAAAAJYigCLb3nvvPYWEhMjT01O1a9fWzz//nGX/+fPnq0KFCvL09FSVKlX0zTffOLxvjNHgwYMVFBQkLy8vRUVFaffu3Xm5C3kqt+cnNjZWNpvN4dWoUaO83IU8k5O5+fXXX/WPf/xDISEhstlsmjBhwm2PeafL7fkZOnTodZ+dChUq5OEe5J2czM20adP02GOPqXDhwipcuLCioqKu65+fv3eyMz/59XtnwYIFqlGjhvz9/eXj46OHHnpIH3/8sUOfe+2z43QGyIa5c+cad3d3M2PGDPPrr7+aTp06GX9/f3Ps2LFM+ycmJhpXV1fz9ttvm+3bt5tBgwaZAgUKmK1bt9r7jB492vj5+ZlFixaZLVu2mKefftqUKVPG/PHHH1btVq7Ji/mJiYkxjRo1MkeOHLG/Tp06ZdUu5Zqczs3PP/9s+vbta/773/+awMBA8+677972mHeyvJifIUOGmEqVKjl8dn7//fc83pPcl9O5adOmjXnvvffMpk2bTHJysomNjTV+fn7m8OHD9j75+XsnO/OTX793Vq5caRYsWGC2b99u9uzZYyZMmGBcXV3NkiVL7H3upc/OnYAAimypVauWeeWVV+zL6enppmTJkmbUqFGZ9n/++edN06ZNHdpq165tunTpYowxJiMjwwQGBpqxY8fa3z9z5ozx8PAw//3vf/NgD/JWbs+PMX/9j6B58+Z5Uq+Vcjo31ypdunSmAet2xrzT5MX8DBkyxFStWjUXq3SO2/05X7lyxRQsWNDMnDnTGMP3zt/9fX6M4XvnWg8//LAZNGiQMebe++zcCTgFj5u6dOmSNmzYoKioKHubi4uLoqKitGbNmkzXWbNmjUN/SYqOjrb3379/v44ePerQx8/PT7Vr177hmHeqvJifqxISElSiRAmFhobq5Zdf1smTJ3N/B/LQrcyNM8Z0lrzcl927d6tkyZIqW7as2rZtq5SUlNst11K5MTcXLlzQ5cuXVaRIEUl87/zd3+fnqvz+vWOM0YoVK7Rz5049/vjjku6tz86dggCKmzpx4oTS09MVEBDg0B4QEKCjR49mus7Ro0ez7H/1vzkZ806VF/MjSY0aNdKsWbO0YsUKjRkzRqtWrVLjxo2Vnp6e+zuRR25lbpwxprPk1b7Url1b8fHxWrJkiaZMmaL9+/frscce09mzZ2+3ZMvkxtz0799fJUuWtIeG/P6983d/nx8pf3/vpKamytfXV+7u7mratKkmTZqkhg0bSrq3Pjt3CjdnFwAgc//85z/tf65SpYrCw8P1wAMPKCEhQQ0aNHBiZbjTNW7c2P7n8PBw1a5dW6VLl9a8efP00ksvObEy64wePVpz585VQkKCPD09nV3OHedG85Ofv3cKFiyozZs369y5c1qxYoV69+6tsmXLKjIy0tml3ZM4AoqbKlasmFxdXXXs2DGH9mPHjikwMDDTdQIDA7Psf/W/ORnzTpUX85OZsmXLqlixYtqzZ8/tF22RW5kbZ4zpLFbti7+/v8qXL59vPjvjxo3T6NGj9d133yk8PNzent+/d6660fxkJj9977i4uOjBBx/UQw89pD59+ujZZ5/VqFGjJN1bn507BQEUN+Xu7q7q1atrxYoV9raMjAytWLFCERERma4TERHh0F+Sli1bZu9fpkwZBQYGOvRJS0tTUlLSDce8U+XF/GTm8OHDOnnypIKCgnKncAvcytw4Y0xnsWpfzp07p7179+aLz87bb7+tYcOGacmSJapRo4bDe/n9e0fKen4yk5+/dzIyMnTx4kVJ99Zn547h7LugcHeYO3eu8fDwMPHx8Wb79u2mc+fOxt/f3xw9etQYY0z79u3NgAED7P0TExONm5ubGTdunElOTjZDhgzJ9DFM/v7+5osvvjC//PKLad68+V37SIvcnp+zZ8+avn37mjVr1pj9+/eb5cuXm2rVqply5cqZP//80yn7eKtyOjcXL140mzZtMps2bTJBQUGmb9++ZtOmTWb37t3ZHvNukhfz06dPH5OQkGD2799vEhMTTVRUlClWrJg5fvy45ft3O3I6N6NHjzbu7u7ms88+c3iM0NmzZx365NfvnZvNT37+3hk5cqT57rvvzN69e8327dvNuHHjjJubm5k2bZq9z7302bkTEECRbZMmTTKlSpUy7u7uplatWmbt2rX29+rVq2diYmIc+s+bN8+UL1/euLu7m0qVKpnFixc7vJ+RkWHeeOMNExAQYDw8PEyDBg3Mzp07rdiVPJGb83PhwgXz5JNPmuLFi5sCBQqY0qVLm06dOt2VAcuYnM3N/v37jaTrXvXq1cv2mHeb3J6fVq1amaCgIOPu7m7uu+8+06pVK7Nnzx4L9yj35GRuSpcunencDBkyxN4nP3/v3Gx+8vP3zuuvv24efPBB4+npaQoXLmwiIiLM3LlzHca71z47zmYzxhhrj7kCAAAgP+MaUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKUIoAAAALAUARQAAACWIoACAADAUgRQAMjHEhISZLPZdObMmXy1bQDORQAFgHwiMjJSvXr1cmirU6eOjhw5Ij8/v3t221kpU6aMli9f7rTtA/mVm7MLAAA4j7u7uwIDA/PdtiXpl19+0enTp1WvXj2n1QDkVxwBBYBckJGRoVGjRqlMmTLy8vJS1apV9dlnn0mSjDGKiopSdHS0jDGSpFOnTun+++/X4MGD7WNMnz5dYWFh8vT0VIUKFfT+++87bOPw4cNq3bq1ihQpIh8fH9WoUUNJSUmSpNjYWLVo0cKhf69evRQZGWl/f9WqVZo4caJsNptsNpsOHDiQ6Wnwzz//XJUqVZKHh4dCQkI0fvx4h3FDQkI0cuRIvfjiiypYsKBKlSqlqVOn3nBusrvt+Ph4+fv76+uvv1ZoaKi8vb317LPP6sKFC5o5c6ZCQkJUuHBh9ejRQ+np6fbxL168qL59++q+++6Tj4+PateurYSEhJv+zL744gs1atRIBQoUuO49Y4yGDh2qUqVKycPDQyVLllSPHj1uOiaAbDIAgNs2fPhwU6FCBbNkyRKzd+9eExcXZzw8PExCQoIxxpjDhw+bwoULmwkTJhhjjHnuuedMrVq1zOXLl40xxsyePdsEBQWZzz//3Ozbt898/vnnpkiRIiY+Pt4YY8zZs2dN2bJlzWOPPWZ+/PFHs3v3bvPpp5+a1atXG2OMiYmJMc2bN3eoqWfPnqZevXrGGGPOnDljIiIiTKdOncyRI0fMkSNHzJUrV8zKlSuNJHP69GljjDHr1683Li4u5q233jI7d+40cXFxxsvLy8TFxdnHLV26tClSpIh57733zO7du82oUaOMi4uL2bFjR6Zzk91tx8XFmQIFCpiGDRuajRs3mlWrVpmiRYuaJ5980jz//PPm119/NV999ZVxd3c3c+fOtY/fsWNHU6dOHfPDDz+YPXv2mLFjxxoPDw+za9euLH9mNWrUMJ988kmm782fP98UKlTIfPPNN+bgwYMmKSnJTJ06NcvxAGQfARQAbtOff/5pvL297WHwqpdeesm0bt3avjxv3jzj6elpBgwYYHx8fBwC0gMPPHBdGBo2bJiJiIgwxhjz4YcfmoIFC5qTJ09mWsPNAqgxxtSrV8/07NnToc/fQ2CbNm1Mw4YNHfr069fPVKxY0b5cunRp065dO/tyRkaGKVGihJkyZUqmtWV323FxcUaS2bNnj71Ply5djLe3tzl79qy9LTo62nTp0sUYY8zBgweNq6ur+e233xzGbtCggRk4cOAN6zl8+LBxd3e3b/vvxo8fb8qXL28uXbp0wzEA3DquAQWA27Rnzx5duHBBDRs2dGi/dOmSHn74Yfvyc889p4ULF2r06NGaMmWKypUrJ0k6f/689u7dq5deekmdOnWy979y5Yr9Bp3Nmzfr4YcfVpEiRfJ0X5KTk9W8eXOHtkcffVQTJkxQenq6XF1dJUnh4eH29202mwIDA3X8+PHb3r63t7ceeOAB+3JAQIBCQkLk6+vr0HZ1W1u3blV6errKly/vMM7FixdVtGjRG27nyy+/VN26deXv75/p+88995wmTJigsmXLqlGjRmrSpImaNWsmNzf+twnkBv4mAcBtOnfunCRp8eLFuu+++xze8/DwsP/5woUL2rBhg1xdXbV79+7r1p82bZpq167tsP7VwOfl5ZVlDS4uLvbrS6+6fPlyDvck+/5+3aTNZlNGRkaejJvVts6dOydXV1f7vF7r2tD6d19++aWefvrpG74fHBysnTt3avny5Vq2bJm6deumsWPHatWqVZleMwogZwigAHCbKlasKA8PD6WkpGR5R3WfPn3k4uKib7/9Vk2aNFHTpk31xBNPKCAgQCVLltS+ffvUtm3bTNcNDw/X9OnTderUqUyPghYvXlzbtm1zaNu8ebNDWHJ3d3e4eSczYWFhSkxMdGhLTExU+fLlrwt4OZGdbd+Khx9+WOnp6Tp+/Lgee+yxbK1z7tw5rVy5UlOmTMmyn5eXl5o1a6ZmzZrplVdeUYUKFbR161ZVq1YtN0oH8jUCKADcpoIFC6pv37569dVXlZGRobp16yo1NVWJiYkqVKiQYmJitHjxYs2YMUNr1qxRtWrV1K9fP8XExOiXX35R4cKF9eabb6pHjx7y8/NTo0aNdPHiRa1fv16nT59W79691bp1a40cOVItWrTQqFGjFBQUpE2bNqlkyZKKiIjQE088obFjx2rWrFmKiIjQ7NmztW3bNodLAEJCQpSUlKQDBw7I19c30yDbp08f1axZU8OGDVOrVq20Zs0aTZ48+bo78nMqO9u+FeXLl1fbtm3VoUMHjR8/Xg8//LB+//13rVixQuHh4WratOl16yxZskTly5dXSEjIDceNj49Xenq6ateuLW9vb82ePVteXl4qXbp0rtQN5Hc8hgkAcsGwYcP0xhtvaNSoUQoLC1OjRo20ePFilSlTRr///rteeuklDR061H707M0331RAQIC6du0qSerYsaOmT5+uuLg4ValSRfXq1VN8fLzKlCkj6a8jiN99951KlCihJk2aqEqVKho9erT9qGR0dLTeeOMNvfbaa6pZs6bOnj2rDh06ONTYt29fubq6qmLFiipevLhSUlKu249q1app3rx5mjt3ripXrqzBgwfrrbfeUmxs7G3NT3a2favi4uLUoUMH9enTR6GhoWrRooXWrVunUqVKZdr/iy++yPL0uyT5+/tr2rRpevTRRxUeHq7ly5frq6++yvK6UgDZZzN/v2gIAIB71JUrVxQQEKBvv/1WtWrVcnY5QL7FEVAAQL5x6tQpvfrqq6pZs6azSwHyNY6AAgAAwFIcAQUAAIClCKAAAACwFAEUAAAAliKAAgAAwFIEUAAAAFiKAAoAAABLEUABAABgKQIoAAAALEUABQAAgKX+H3lJgqQp8RuHAAAAAElFTkSuQmCC", "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\");" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Conclusions:\n", "\n", "1. As expected, the best results with JAX are obtained by JIT compiling function and gradient and using the gradient in the minimization. However, the performance of the Numba JIT compiled function is comparable even without computing the gradient.\n", "\n", "1. JIT compiling the cost function with JAX but not using the gradient also gives good performance, but worse than using Numba for the same.\n", "\n", "3. Combining the JAX JIT with the JAX gradient calculation is very important. Using only the Python-computed gradient even 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." ] }, { "attachments": {}, "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": 31, "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", ")" ] }, { "attachments": {}, "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": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "30.4 ms ± 1.37 ms 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": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "89 ms ± 28.6 ms 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": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "187 µs ± 28.1 µ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)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "It is much faster... but only because the compilation cost is excluded here." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "496 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)" ] }, { "attachments": {}, "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." ] }, { "attachments": {}, "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": 36, "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": 37, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:25:43.646212Z", "start_time": "2020-02-21T10:25:43.512384Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0J0lEQVR4nO3dfXSU9Z3//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": 38, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:25:44.032210Z", "start_time": "2020-02-21T10:25:43.648365Z" } }, "outputs": [ { "data": { "text/plain": [ "Array(876.49545695, dtype=float64)" ] }, "execution_count": 38, "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": 39, "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", "
Migrad
FCN = 23.14 Nfcn = 91
EDM = 3.12e-05 (Goal: 0.0002)
Valid Minimum Below EDM threshold (goal x 10)
No parameters at limit Below call limit
Hesse ok Covariance accurate
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \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.039 (-0.530) -0.150 (-0.657)
x1 -0.039 (-0.530) 0.24 0.17 (0.230)
x2 -0.150 (-0.657) 0.17 (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 │ Below EDM threshold (goal x 10) │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ No parameters at limit │ Below call limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Hesse ok │ Covariance accurate │\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.039 -0.150 │\n", "│ x1 │ -0.039 0.24 0.17 │\n", "│ x2 │ -0.150 0.17 2.32 │\n", "└────┴──────────────────────┘" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = Minuit(cost, np.zeros(3))\n", "m.migrad()" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "ExecuteTime": { "end_time": "2020-02-21T10:25:44.566228Z", "start_time": "2020-02-21T10:25:44.065443Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAG3CAYAAAAU+jfPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYgElEQVR4nO3deZyNdf/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}$\");" ] }, { "attachments": {}, "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 }