{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic Least Squares Fit  Example in Python using iminuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitting example using iminuit (https://iminuit.readthedocs.io). iminuit can be installed with `pip install iminuit`. [Minuit](https://en.wikipedia.org/wiki/MINUIT) is a robust numerical minimization program written by CERN physicist Fred James in 1970s. It is widely used in particle physics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from iminuit import Minuit\n",
    "from pprint import pprint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "xd = np.array([1., 2., 3., 4., 5.])\n",
    "yd = np.array([1.7, 2.3, 3.5, 3.3, 4.3])\n",
    "yd_err = np.array([0.5, 0.3, 0.4, 0.4, 0.6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define fit function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x, theta0, theta1):\n",
    "    return theta0 + theta1*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define $\\chi^2$ function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Minuit finds the minimum of a multi-variate function. We need to define a $\\chi^2$ function which is then minimized by iminuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi2(theta0, theta1):\n",
    "    fy = f(xd, theta0, theta1)\n",
    "    diffs = (yd - fy) / yd_err\n",
    "    return np.sum(diffs**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize minuit and perform the fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Minuit(chi2, theta0=1, theta1=0.)\n",
    "m.errordef = Minuit.LEAST_SQUARES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <th colspan=\"5\" style=\"text-align:center\" title=\"Minimizer\"> Migrad </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:left\" title=\"Minimum value of function\"> FCN = 2.296 </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center\" title=\"No. of function evaluations in last call and total number\"> Nfcn = 32 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:left\" title=\"Estimated distance to minimum and goal\"> EDM = 9.35e-23 (Goal: 0.0002) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center\" title=\"No. of gradient evaluations in last call and total number\">  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Valid Minimum </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center;background-color:#92CCA6;color:black\"> No Parameters at limit </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Below EDM threshold (goal x 10) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Below call limit </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\"> Covariance </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\"> Hesse ok </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Is covariance matrix accurate?\"> Accurate </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Is covariance matrix positive definite?\"> Pos. def. </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Was positive definiteness enforced by Minuit?\"> Not forced </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th title=\"Variable name\"> Name </th>\n",
       "        <th title=\"Value of parameter\"> Value </th>\n",
       "        <th title=\"Hesse error\"> Hesse Error </th>\n",
       "        <th title=\"Minos lower error\"> Minos Error- </th>\n",
       "        <th title=\"Minos upper error\"> Minos Error+ </th>\n",
       "        <th title=\"Lower limit of the parameter\"> Limit- </th>\n",
       "        <th title=\"Upper limit of the parameter\"> Limit+ </th>\n",
       "        <th title=\"Is the parameter fixed in the fit\"> Fixed </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 0 </th>\n",
       "        <td> theta0 </td>\n",
       "        <td> 1.2 </td>\n",
       "        <td> 0.5 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 1 </th>\n",
       "        <td> theta1 </td>\n",
       "        <td> 0.61 </td>\n",
       "        <td> 0.15 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th> theta0 </th>\n",
       "        <th> theta1 </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> theta0 </th>\n",
       "        <td> 0.211 </td>\n",
       "        <td style=\"background-color:rgb(131,131,250);color:black\"> -0.0646 <strong>(-0.919)</strong> </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> theta1 </th>\n",
       "        <td style=\"background-color:rgb(131,131,250);color:black\"> -0.0646 <strong>(-0.919)</strong> </td>\n",
       "        <td> 0.0234 </td>\n",
       "    </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "┌─────────────────────────────────────────────────────────────────────────┐\n",
       "│                                Migrad                                   │\n",
       "├──────────────────────────────────┬──────────────────────────────────────┤\n",
       "│ FCN = 2.296                      │              Nfcn = 32               │\n",
       "│ EDM = 9.35e-23 (Goal: 0.0002)    │                                      │\n",
       "├──────────────────────────────────┼──────────────────────────────────────┤\n",
       "│          Valid Minimum           │        No Parameters at limit        │\n",
       "├──────────────────────────────────┼──────────────────────────────────────┤\n",
       "│ Below EDM threshold (goal x 10)  │           Below call limit           │\n",
       "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n",
       "│  Covariance   │     Hesse ok     │ Accurate  │  Pos. def.  │ Not forced │\n",
       "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n",
       "┌───┬────────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n",
       "│   │ Name   │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │\n",
       "├───┼────────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n",
       "│ 0 │ theta0 │    1.2    │    0.5    │            │            │         │         │       │\n",
       "│ 1 │ theta1 │   0.61    │   0.15    │            │            │         │         │       │\n",
       "└───┴────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n",
       "┌────────┬─────────────────┐\n",
       "│        │  theta0  theta1 │\n",
       "├────────┼─────────────────┤\n",
       "│ theta0 │   0.211 -0.0646 │\n",
       "│ theta1 │ -0.0646  0.0234 │\n",
       "└────────┴─────────────────┘"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.migrad()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta0 = 1.162066 +/- 0.459550\n",
      "theta1 = 0.613945 +/- 0.153005\n"
     ]
    }
   ],
   "source": [
    "for p in m.parameters:\n",
    "    print(f\"{p} = {m.values[p]:.6f}\" \\\n",
    "          f\" +/- {m.errors[p]:.6f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th> theta0 </th>\n",
       "        <th> theta1 </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> theta0 </th>\n",
       "        <td> 0.211 </td>\n",
       "        <td style=\"background-color:rgb(131,131,250);color:black\"> -0.0646 <strong>(-0.919)</strong> </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> theta1 </th>\n",
       "        <td style=\"background-color:rgb(131,131,250);color:black\"> -0.0646 <strong>(-0.919)</strong> </td>\n",
       "        <td> 0.0234 </td>\n",
       "    </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "┌────────┬─────────────────┐\n",
       "│        │  theta0  theta1 │\n",
       "├────────┼─────────────────┤\n",
       "│ theta0 │   0.211 -0.0646 │\n",
       "│ theta1 │ -0.0646  0.0234 │\n",
       "└────────┴─────────────────┘"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# covariance matrix\n",
    "m.covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th> theta0 </th>\n",
       "        <th> theta1 </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> theta0 </th>\n",
       "        <td> 1 </td>\n",
       "        <td style=\"background-color:rgb(131,131,250);color:black\"> -0.919 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> theta1 </th>\n",
       "        <td style=\"background-color:rgb(131,131,250);color:black\"> -0.919 </td>\n",
       "        <td> 1 </td>\n",
       "    </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "┌────────┬───────────────┐\n",
       "│        │ theta0 theta1 │\n",
       "├────────┼───────────────┤\n",
       "│ theta0 │      1 -0.919 │\n",
       "│ theta1 │ -0.919      1 │\n",
       "└────────┴───────────────┘"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# correlation matrix\n",
    "m.covariance.correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.21118631 -0.06460345]\n",
      " [-0.06460345  0.02341046]]\n"
     ]
    }
   ],
   "source": [
    "print(repr(m.covariance))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot data along with fit function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "xf = np.linspace(0., 6., 1000)\n",
    "theta0 = m.values[\"theta0\"]\n",
    "theta1 = m.values[\"theta1\"]\n",
    "yf = f(xf, theta0, theta1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbyElEQVR4nO3de5TVZb3H8fcXzRQxO+W1C5jHFPOek2mRYGparUzNWyJlkVOp1clTKVErrSBPmWaeOgWdUpC8ZC3vpnnyUioqpKcQNU0FK1EU0KPIcHvOH8+Me4ABGZzZz2/v/X6tNYuZ3+yZ+e4lzoff7/M8vx0pJSRJqpoBpQeQJKknBpQkqZIMKElSJRlQkqRKMqAkSZVkQEmSKsmAkiRVUtGAioitI+KCiJgbEYsiYmZEDC85kySpGtYv9YMj4rXAbcAfgQ8Cc4FtgadKzSRJqo4odSeJiBgPDE8pvbvIAJKkSit5ie9Q4M6IuCQinoqIeyPi5IiIgjNJkiqi5BnUos53zwEuBXYHzgNOSyn9Zw+PbwfaATbeeOM9hw4dWqdJJUl9Zfr06U+nlDZfm8eWDKjFwLSU0ru6HRsPHJZS2nFNX9vW1pamTZvW3yNKkvpYRExPKbWtzWNLXuJ7Api50rH7gcEFZpEkVUzJgLoN2GGlY9sDswrMIkmqmJIBdQ6wd0SMjYjtIuJI4PPAjwrOJEmqiGIBlVK6m7yS7yhgBjAO+Drw41IzSZKqo9hGXYCU0jXANSVnkCRVk/fikyRVkgElSaokA0qSVEkGlCSpkgwoSVIlGVCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSQaUJKmSDChJUiUZUJKkSjKgJEmVZEBJkirJgJIkVZIBJUmqJANKklRJBpQkqZIMKElSJRlQkqRKMqAkSZVkQEmSKsmAkiRVkgElSaqkYgEVEadHRFrpbU6peSRJ1bJ+4Z//IDCi28fLCs0hSaqY0gG1NKXkWZMkaRWlO6htI+KfEfFoRFwcEduu7oER0R4R0yJi2ty5c+s5oySpgJIBdSdwPHAwcAKwFXB7RLy+pwenlCaklNpSSm2bb755/aaUJBVR7BJfSum67h9HxFTgEeDjwNlFhpIkVUbpS3wvSSk9D9wHvLX0LJKk8ioTUBGxITAUeKL0LJLU30aMyG9avZL7oM6KiOER8ZaIeCdwGbAxcEGpmSRJ/eC55+AXv4D99+/Vl5VcZv4m4CJgM2AuMBXYO6U0q+BMkqS+sHQp/O53MHkyXH45vPgibLddr75FyUUSx5T62ZKkfpAS3HNPDqWLLoInn4TXvQ4+8QkYNQre+U4YsPYX7kpv1JUkNbrHH4cpU3IwzZwJG2wAH/oQHHccfOAD+eN1YEBJknrvuefg17/OoXTzzfns6d3vhp/8BI46Cv7lX17xjzCgJElrZ+lSuOGGWq+0aFHulU4/PZ8tbbvamwGtEwNKkrR6KcGf/lTrlZ56KvdKn/xkrVeK6JcfbUBJklY1e3atV7r//lqvNGoUvP/969wr9YYBJUnKnnsOLrssh9Itt+Szp2HD4Kc/hSOP7JNeqTcMKElqZUuW1HqlK67IvdJb3wpnnJF7pbe8pdhoBpQktZquXmnSpNwrzZ0Lr389jB6dL+HttVe/9Uq9YUBJUqvoqVc65JAcSgcfXJdeqTcMKElqZt17pZtvzseGDYMJE3Kv9NrXlpxujQwoSWo2q+uVvvnN4r1SbxhQktQMUoLp02v7lSraK/WGASVJjWzWrFqv9MAD8OpX1/YrVbBX6g0DSpIazbPPrrhfCeA974FTTql8r9QbBpQkNYIlS+D663MoXXll7pW23x6+9S0YObJheqXeMKAkqc6mTIGpU6GjA7bZBsaNyxmzipRg2rQcShdfnHulzTaDT30qX8J7xzsarlfqDQNKkupoyhRob8/hBLlCam/P778UUj31St33K73qVUVmrzcDSpLqaOxYWLhwxWMLF8LYMcsZuegXK/ZK++4L//7vcMQRTdMr9YYBJUl1NHv2ao4/Tr5019UrHXdcvv7XwgwoSaqjwW9OzJq9am80eJMFcOOdTd8r9caA0gNIUkuYNQvGjWPc4i8xkBdW+NTAgYlx//W6htxM258MKEnqL88+Cz/7GQwfni/Xfe1rjNx+GhNG38lGr14GwJAhMGFC9LyKr8V5iU+S+tKSJfDb39b2K3V0wA47wLe/nZfpbbMNI4GJD+eHd92/VasyoCTplUoJ7r67tl/p6afzfqX29rw0vK3NS3frwICSpHX12GO1/UoPPpj3K334wzmUDjqoZfYr9ZfKBFREjAHGAz9KKZ1ceh5J6tGCBbX74N16az42fDh8+ct5v9KmmxYdr5lUIqAiYm+gHfhz6VkkaRVr0Sup7xUPqIjYFJgCfBL4RuFxJCnrqVfafHN7pToqHlDABOCylNJNEWFASSrrscfgwgtzMP31r/ZKBRUNqIg4AdgOOG4tHttOvgzI4MGD+3kySS1lwQL41a9yKP3hD/nY8OHwla/YKxVULKAiYgfyoohhKaUlL/f4lNIE8tkWbW1tqZ/Hk9TsFi+u9UpXXZV7paFDa699MWRI6QlbXskzqH2AzYD7onYddz1g34j4DLBxSqmj1HCSmlBKcNddtV7pmWdyr/TpT+dLeHvuaa9UISUD6nJg2krHfgE8RD6zWlzvgaTSRozIf3p3gT726KO5V7rwwtwrbbhhrVd63/vslSqqWECllBYAC7ofi4gXgHkppRklZpLURObPr/VKf/xjPjZiBJx6KnzkI/ZKDaAKq/gkqW9075WuvDJ/vOOOMH587pVcYNVQKhVQKaURpWeQ1GB66pW22AI++9l8Ce/tb7dXalCVCihJWmtdvdLkyfDQQ/ZKTciAktQ4VtcrnXaavVITMqAkVdvixXDddbX9SvZKLcOAklRUj0vrU4I778yhdMklTdkruZXg5RlQkqrjkUdq+5W6eqVDD82hdOCB9kotxoCSVNSgJfPZb+6lMGwy3HZbPjMaMQLGjMm90mteU3pEFWJASaq/rl5p0iR+c8fVbJAWw6veBt/5Tu6V3vzm0hOqAgwoSfWREkydmi/fdeuVrnjDidyw5SgmTtuj4Xsl9S0DSqqIKVPy7++OjvwCrV031W54Xb3S5Mnw8MO1XuljH4MDD+RHB3T+GjKbtBIDSqqAKVPyC7V2dN6/f9as/DE0aEjNnw+XXppDqXuv9NWv2itprRlQUgWMHQsLF654bOHCfLxhAmrxYrj22hxKV1+dP36bvZLWnQElVcDs2b07XhldvVLXfqV58/J+pRNPzEvD97BX0rozoKQKGDw4X9br6Xgl/e1vtf1KDz8MG2204n6l9f3VolfOv0VSBYwblzun7pf5Bg7Mxytj3rxar3T77fnMaL/98nXIww+3V1KfM6CkCujqmUaPzgslhgypyCq+jo5ar3TNNblX2mknOPNMOPZYeyX1KwNKqoiRI2HixPx+0fu0pQR33JFD6dJL85nTllvCSSflS3i7795nvVLTLq1XnzCgJGVdvdLkyfn9jTaCww7LoXTAAX3eKzXd0nr1OQNKamU99UrvfS98/eu5V9pkk3770U2xtF79yoCSWs2aeqWRI+FNb6rLGA27tF51Y0BJraB7r3TJJflOD/3UK62thltar7ozoKRm9re/5VC68MK69Eq90RBL61WUASU1m3nz8lnS5Mn5rKmOvVJvVHZpvSrDgJKaQUdH7pO6eqUlS2DnneE//iPvV6pTr9RblVlar0oyoKRGlVJeede1X2n+fNhqK/jc5/IlvN128z54amgGlNRoHn641is98kgubrp6pf339z54ahrF/iZHxEnAp4FtOg/dB3w7pXRNqZmk0lZ7meuZZ2r7lbp6pf33h298I4dTRXolqS/1KqAi4q/AfwMXpJTmvMKf/XfgVOAhYADwceDyiNgzpfTnV/i9pca3ul7pu9/NvdIb31h6Qqlf9fYMagnwHeBbEXEt8DPg2pTS8t7+4JTSFSsdGhsRnwX2AQwotSZ7JeklvQqolNJOEbE3MBo4CvgQMCcizgd+nlL627oMERHrAUcCg4Db1+V7SA3toYdqr69kryQB69BBpZSmAlMj4gvA0eSwGgOcFhG3kM+qfp1S6ni57xURuwB3ABsCzwOHpZT+sprHtgPtAIPdaq5m8Mwztf1KU6eu2CsdfjgMGlR6QqmoSCm98m8SsT3wDeCjQAIWAJOBs1NKq72zVkRsAAwGNgWOAE4ARqSUZqzp57W1taVp06a94rmluuvogKuvzqF07bW5V9pll3ym1KK90ogR+U/3QbWGiJieUmpbm8e+ousGnZfmDiGfRR1MDqebgA7gZOCEiDi2h74JgJTSYuDhzg+nR8Q7gC92fj+pOaQEt91W65UWLICtt4bPf77WK0laxToFVEQMJYfIKGAL4CngLGBiVw8VEdsBlwLfBXoMqB4MAF69LjNJlfPQQ7X9So8+mnulww+v9UrrrVd6wkrwzEmr09tl5qOBTwJ7dx66EZgAXJFSWtr9sSmlhyPih+ROqqfvdSZwDfA4sAlwLDAC+GBvZpIqZeVeacCAHEZnnJEXPdgrSWutt2dQE4E5wJnks6XHXubxM8ldVE+2Ai7s/PNZ8tLy96eUru/lTFJZPfVKu+4K3/te7pXe8IbSE0oNqbcBdThwVUpp2do8OKV0F3DXaj53fC9/tlQdq+uVvvCFfAlv111LTyg1vN7ug7q8n+aQGsPqeqWPfSy/pIW9ktRn3P0nvZynn671Snfeaa8k1YkBJfVk0aIVe6WlS+2VpDozoKQuy5fXeqVf/arWK/3bv9krSQUYUNJf/1rrlR57DDbeuLZfyV5JKsaAUmvq6pUmTYK77sq90gEHwLe+BYceaq8kVYABpdbRU6+0225w1lnw0Y/aK0kVY0CpsvrkJqLde6VLL4Vnn81B9MUv5kt4u+zSB5NK6g8GlJpTT73SRz6SQ2m//eyVpAZgQKl5PP00XHxxDqaVe6XDDsshJalhGFBqbIsWwVVX5VC67roVe6Vjj83LxCU1JANKjWf5cvjjH2v7leyVpKZkQKlxPPhgrVeaNcteSWpyBpQqbdPFc+G8zl7p7rtzr3TggTBuXN6vZK8kNS0DStWzaBFceSXj/zKZveb/Fu5YCrvvDt//ft6vZK8ktQQDStXQvVe69FKmPPdBTuO/+DtvZMjWSxj3pQ0YObL0kJLqyYBSWSv3SoMGMWWP79F+1wks7Mid0qwnNqC9PT/ckJJax4DSA6gFzZ0L550He+0FQ4fCd76T/7zwQpgzh7GzP/NSOHVZuBDGji00r6QiPINSfbz4Ym2/0m9/m/crraZXmj2752+xuuOSmpMBpf6zfDn84Q+1/UrPPQdvfCOcckpeGr7zzj1+2eDB+WpfT8cltQ4DSn3vgQdyKE2Z8lKv9NJ+pREjXna/0rhx0N6eL+t1GTgwH5fUOgwo9Y2nnqrdB2/atLxf6X3vg/Hj836lgQPX+lt1LYQYPRo6OmDIkBxOLpCQWosB1QD65GUn+sOLL8KVV9Z6pWXLYI894Oyzc6+01Vbr/K1HjoSJE/P7lXvekurCgFLvdPVKkybBZZflXulNb4IvfSlfwttpp9ITSmoSBpTWzv3313ql2bNzr3TEETmUhg/3PniS+lyxgIqIMcDhwA5ABzAVGJNSmlFqJq1kdb3SmWfChz/cq15Jknqr5BnUCODHwN1AAN8EboyIt6WU5hWcq7X1Y68kSb1RLKBSSgd1/zgiRgHPAu8GrioyVKtavhxuvTWHkr2SpIqoUge1CfnWS/NLD9Iy1tQrjRiRL+lJUiFVCqhzgXuBO3r6ZES0A+0Ag72lwLp76im46KIcTNOn58UNFe2VXF4utbZKBFREnA0MA4allJb19JiU0gRgAkBbW1uq43iN78UX4Yorcihdf33uld7+djjnnNwrbbll6QklaRXFAyoizgGOAfZLKT1Sep6msXw53HJLrVf6v/+DN78ZvvzlfAnvbW8rPaEkrVHRgIqIc4GjyeH0QMlZmsbMmbVe6fHHYZNNVtyvZK8kqUGU3Af1I2AUcCgwPyK61i8/n1J6vtRcDamrV5o0Cf70p9wrHXQQfPe7cMghleqVJGltlTyDOrHzz/9Z6fgZwOn1HaUB9dQr7bkn/OAHcMwx9kqSGl7JfVBR6mc3LHslSS2k+CIJvbwhL8zkfU9Ohm3slSS1Dn+zVdWTT8IPfsCUbcZyy7SBjHp8HNvMvZspJ90Gc+bAz38O++1nOElqWp5BVcnChbVe6YYbmLLsKNrjv1nIRgDMWrQl7b/YEvbxxfskNT//+V3a8uXw+9/DJz6Rb8R67LEwYwZ85SuM3fp8FqaNVnj4woUwdmyhWSWpjjyDKuW++2r7lf7+99wrHXlk7pX23RcGDGD2mT1/6ezZ9R1VkkowoOppzpzaffDuuSfvVzr4YDjrrLxfaaMVz5YGD4ZZs1b9Nt6KUFIr8BJff1u4MIfSBz6QX8LilFNyMJ17Lvzzn3D11XD00auEE8C4cavusR04MB+XpGbnGVR/WL4834p78mT49a/zfqXBg+HUU+G442DHHdfq23QthBg9Gjo6YMiQHE4ukJDUCgyovrRyr/Sa18BRR+VQ6uyVemvkSJg4Mb/vy09IaiUG1CvVy15JkrR2DKh1sXAhXH75S/uVWL4c2tpyr3TMMbDFFqUnlKSGZ0CtrWXLVuyVnn8+90qnnZaXhg8dWnpCSWoqBtTLmTGj1iv94x+5Vzr66BxK73mPtxqSpH5iQPVkzhz45S9zMN17L6y/fu6Vzj4bPvQheyVJqgMDqssLL9R6pd/9LvdK73gH/PCHuVfafPPSE0pSS2ntgOrqlSZNgt/8JvdKQ4bAmDF5abi9kiQV05oBtXKvtOmm+Sxp1CgYNsxeSZIqoHUC6oknavuVuvdK55yTe6UNNyw9oSSpm+YOKHslSWpYzRdQy5bBTTflULJXkqSG1TwB9Ze/5FD65S/tlSSpCTR2QD3xRG2/0v/+b+6V3v9+eyVJagKNGVDz5sFBB8GNN+Zeaa+94Lzz8h0emrBX8i7mklpRYwbUo49CSvDVr+ZeaYcdSk8kSepjjRlQO+wAM2faK0lSE2vM3/CDBhlOktTkiv6Wj4h9I+LKiPhHRKSIOL7kPJKk6ih9GjIImAF8AXix8CySpAop2kGllK4FrgWIiPNLziJJqpbSZ1BrLSLaI2JaREybO3du6XEkSf2sYQIqpTQhpdSWUmrbvAn3OkmSVtQwASVJai0GlCSpkgwoSVIlFV3FFxGDgO06PxwADI6I3YF5KaXZxQaTJBVX+gyqDbin820j4IzO979ZcihJUnml90HdDETJGSRJ1VT6DEqSpB4ZUJKkSjKgJEmVZEBJkirJgJIkVZIBJUmqJANKklRJBpQkqZIMKElSJRlQkqRKMqAkSZVkQEmSKsmAkiRVkgElSaokA0qSVEkGlCSpkgwoSVIlGVCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSQaUJKmSDChJUiUVD6iIODEiHo2IRRExPSLeU3omSVJ5RQMqIo4GzgXGA3sAtwPXRcTgknNJksorfQZ1CnB+SmliSun+lNLngCeAzxaeS5JUWLGAiogNgD2BG1b61A3Au+o/kSSpStYv+LM3A9YDnlzp+JPAASs/OCLagfbODzsiYkb/jlc5mwFPlx6igFZ83j7n1tGKz3uHtX1gyYDqlZTSBGACQERMSym1FR6prlrxOUNrPm+fc+toxecdEdPW9rElO6ingWXAlisd3xKYU/9xJElVUiygUkqLgenAgSt96kDyaj5JUgsrfYnvbGByRNwF3AZ8BngD8JOX+boJ/T1YBbXic4bWfN4+59bRis97rZ9zpJT6c5CXHyDiROArwNbADOCLKaVbiw4lSSqueEBJktST0ht1JUnqkQElSaqkhguoVru5bETsGxFXRsQ/IiJFxPGlZ+pPETEmIu6OiOciYm5EXBURO5eeq79FxEkR8efO5/1cRNwRER8sPVc9df63TxHxn6Vn6S8RcXrnc+z+1hLbaiJi64i4oPP/60URMTMihq/paxoqoFr05rKDyItHvgC8WHiWehgB/Jh8u6v3AkuBGyPidSWHqoO/A6cCbwfagN8Dl0fErkWnqpOI2Jt8p5g/l56lDh4kLwrretul7Dj9LyJeS16pHcAHgR2BzwFPrfHrGmmRRETcCfw5pXRCt2MPAZellMaUm6w+IuJ54OSU0vmlZ6mXiBgEPAscmlK6qvQ89RQR84AxKaWflp6lP0XEpsCfgE8B3wBmpJROLjtV/4iI04EjUkpNf1Wgu4gYDwxPKb27N1/XMGdQ3ly2ZW1C/ns6v/Qg9RIR60XEMeSz51bYtD6B/I/Mm0oPUifbRsQ/O6uKiyNi29ID1cGhwJ0RcUlEPBUR90bEyRERa/qihgko1nxz2a3qP47q5FzgXuCOwnP0u4jYpfMsuYO8Wf2wlNJfCo/VryLiBGA74GulZ6mTO4HjgYOBE8i/u26PiNeXHKoOtgVOBB4BDiL/f30mcNKavqj0nSSk1YqIs4FhwLCU0rLS89TBg8DuwKbAEcAFETEipdSUd+6PiB3IffKwlNKS0vPUQ0rpuu4fR8RU8i/tj5PvrNOsBgDTulUx90TEW8kBtdpFMY10BuXNZVtIRJwDfBR4b0rpkdLz1ENKaXFK6eGU0vTO/5HvBb5YeKz+tA/5ysh9EbE0IpYCw4ETOz9+ddnx+l9K6XngPuCtpWfpZ08AM1c6dj+wxgVuDRNQ3ly2dUTEudTC6YHS8xQ0AGjmX9KXk1ew7d7tbRpwcef7i4tMVUcRsSEwlPwLvJndxqqvA7U9MGtNX9Rol/jW9eayDatzFdt2nR8OAAZHxO7AvJTS7GKD9ZOI+BEwilyqzo+Irn7x+c5/bTaliDgTuAZ4nLww5Fjykvum3QuVUloALOh+LCJeIP/dbtbLmmcBVwGzgS2ArwMbAxeUnKsOziF3bWOBS8jbhD4PfHVNX9RQy8yh9W4uGxEjgJ5WN12QUjq+rsPUQUSs7i/kGSml0+s5Sz1FxPnAfuTS/FnyfqDvpZSuLzlXvUXEzTT3MvOLgX3JlzbnAlOBr6eUVr781XQ6N56PJ59JzSZ3T+elNYRQwwWUJKk1NEwHJUlqLQaUJKmSDChJUiUZUJKkSjKgJEmVZEBJkirJgJIkVZIBJUmqJANKklRJBpQkqZIMKKmgiFg/Im6LiBciYuhKn2uPiBQR3yw1n1SS9+KTCouIIeTXfpoFvDOl1BEROwF3k19iZkSLvGCjtALPoKTCUkqzgNHAbsD3I2Ij8ksSLAJGGk5qVZ5BSRURET8GPkt+Ac53AR9JKf2m7FRSOQaUVBGdr646A/hXYGJKqb3wSFJRXuKTqmM3YHDn+ztHRKO94rXUpwwoqQIi4jXARcDTwFhgH+CMokNJhfkvNKkaJgBDgANTSr+PiD2A0yLixpTSTYVnk4qwg5IKi4jRwM+A8SmlsZ3HXkteev4qYNeU0jPFBpQKMaCkgjo3504nh9HwlNLSbp/bB7gVuC6ldEiZCaVyDChJUiW5SEKSVEkGlCSpkgwoSVIlGVCSpEoyoCRJlWRASZIqyYCSJFWSASVJqiQDSpJUSf8PKLIe0zS7+qYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel(\"x\", fontsize=18)\n",
    "plt.ylabel(\"y\", fontsize=18)\n",
    "plt.errorbar(xd, yd, yerr=yd_err, fmt=\"bo\")\n",
    "plt.xlim(0., 6.)\n",
    "plt.ylim(0., 6.)\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.plot(xf, yf, color=\"red\")\n",
    "plt.tight_layout()\n",
    "# plt.savefig(\"basic_chi2_fit_plot1.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the $\\chi^2$ per degree of freedom and p-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2/ndf = 2.3 / 3, p-value = 0.513\n"
     ]
    }
   ],
   "source": [
    "from scipy import stats\n",
    "\n",
    "chi2 = m.fval\n",
    "\n",
    "n_data_points = xd.size\n",
    "n_fit_parameters = 2\n",
    "n_dof = n_data_points - n_fit_parameters\n",
    "pvalue = 1 - stats.chi2.cdf(chi2, n_dof)\n",
    "\n",
    "print(f\"chi2/ndf = {chi2:.1f} / {n_dof}, p-value = {pvalue:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error ellipse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABhtElEQVR4nO3ddXzV1RvA8c/ZveuNBaOb0d2tdAoC0iiKIiAoiIGK/VOwUSxAuhEEVEC6WxjdXaNZ9248vz++YzAXDFgB5/168Vp8z/3ec7Zxn3vqOUpE0DRN0zQAh+yugKZpmpZz6KCgaZqmJdJBQdM0TUukg4KmaZqWSAcFTdM0LZE5uytwr/z8/KR48eLZXQ1NS+L48eMAlC1bNptromkp2717900RyXO3cg9dUChevDgBAQHZXQ1NS6JJkyYAbNiwIVvroWmpUUqdT085PXykaZqmJXroegqalhN9+OGH2V0FTcsQOihoWgZo0aJFdldB0zKEHj7StAywb98+9u3bl93V0LQHpnsKmpYBhg0bBuiJZu3hp3sKmqZpWiIdFDRN07REOihomqZpiXRQ0DRN0xLpiWZNywBffPFFdldB0zKEDgqalgEaNGiQ3VXQtAyhh480LQNs27aNbdu2ZXc1NO2BZVpPQSk1BWgPXBeRSilcV8CPQDsgGugrInsyqz6alpnef/99QO9T0B5+mdlTmAa0SeN6W6B0wr8BwLhMrIumaZqWDpkWFERkExCcRpGOwAwx7AC8lVIFMqs+mqZp2t1l55xCIeDiHV8HJnxP0zRNyyYPxUSzUmqAUipAKRVw48aN7K6OpmnaIys7l6ReAorc8XXhhO8lIyITgAkAtWrVksyvmqbdmzFjxmR3FTQtQ2RnUFgMvKaU+h2oC4SJyJVsrI+m3bdq1apldxU0LUNk5pLUuUATwE8pFQh8AjgCiMh4YBnGctRTGEtSX8ysumhaZluzZg2gD9vRHn6ZFhREpNddrgvwamY9v6ZlpZEjRwI6KGgPv4diolnTNE3LGjooaJqmaYl0UNA0TdMS6aCgaZqmJdKpszUtA/z222/ZXQVNyxA6KGhaBihbtmx2V0HTMoQePtK0DLBkyRKWLFmS3dXQtAemewqalgFGjx4NQIcOHbK5Jpr2YHRPQdM0TUukg4KmaZqWSAcFTdM0LZEOCpqmaVoiPdGsaRlg5syZ2V0FTcsQOijkcGI5gkSOAwkDsQBWEGvCRwsoBzCVAHMZlLkUmMuAuThKOWZ31R8rRYoUuXshTXsI6KCQQ4nEIZG/QNQkULnAXBKUI+AKDmbAEZTZCAzW0xC3DsGW8GhHxFwCzOVQTrXAqTaYSqKUysYWPdrmzZsHQI8ePbK5Jpr2YHRQyIEkPgAJ+wBsZ8G1C8rzPZSDV9qPkTiwngHrScR6EqwnIX47ErvYKOCQG3GsjXKqDU51wFwapfSUUkYZN24coIOC9vDTQSEHEREkYhREzwSHgiifKSjnRul6rFLO4FgeHMtzqz8gImA7D/G7kPidxse4FcZFBz/EuSXKpRU41UUp/aegaZoOCjmL9SBEzwDXZ1CeH6Ec3B/odkopMBc35hjcugEgtksQ/y8StxFi/0Zi5oLyRlxaoFzagFM9lHLKgMZomvYw0kEhJ7GHA6Bcuz1wQEiNMhUygo7rM4jEQtxmJHYlxC5HYhaAyoW4tES5tE8IEKZMqYemaTmTDgo5iUQZH5VbljydUi7g0hLl0tKYk4jbhsQuh9gVSMxCcMiDuLQzAoRjFT1RrWmPAR0UchKJNj6qzOklpEUpZ3BpinJpmtCD2IjELIHoOUj0dDAVRVw6oFw7oczFsrx+Od2CBQuyuwqaliF0UMhJJDbhk+x9R270IFqjXFoj9nCIXY3ELoGocUjUr4hTQ5RbL3BupieoE/j5+WV3FTQtQ+g1iTmJYw3jY9zqTHsKu9iJtcVjtdvuXhhQDrlQbl1w8J2GyrMR5TEUrKeR0NeQG02RyJ8R29VMq+/DYtq0aUybNi27q6FpD0yJSHbX4Z7UqlVLAgICsrsamcYe1BPswSi/lfc9hm+xW5l+ZiOrrx4k1mYhzmYh3m4lzmbBIkYwMCsTxdz98PfMh79Hfkp55qOkRz4KuHrjcJf9CyJWiNuARM+F+M2Ayeg1uD0LTvUfy7mHJk2aALBhw4ZsrYempUYptVtEat2tnO775zDKrScS9g7E7wDn+vf8+MOhFxl5aBGnI69RO7c/fs65cHEw42xyxCnho7ODI+GWaE5HXuNgyAVWXTmQ+Hg3kxMVvApTJ3cp6viVomyugpj+EySUMoNLC5RLC8R6AYmZB9F/IHGrwVQS3HobK5wcPB7456FpWtbSPYUcRiQWuf4kONfHwfvHe3rs5FPrmHhqLX7OnrxbsSNP5C2frsdFWmM5E3Gd05FXORVxlX0h5zkZcQUAT7MLtXL7Uzu3P/X8SlPYLXcq9Y6D2GVI9GywHDBWULl0Qrk/Z+RkesTpnoKW0+mewkNKKRfErRtETUIsh1GOFdP1OKvdxoRTa2mQpwyfV+mBh6NLup/Tw+xCFZ+iVPEpmvi94LhIAoJPszPoNDtvnmL9tcMA+Hvko1WBqrQqUIVCbr531NsZXDujXDsjlgNI1GyIWYDEzEGc6qPc+oJzY51aQ9NyOB0UciDlPhCJWYiEjwTfOekaow+zRCMIDfzK3FNASI2vs0fCi39VRIQL0UHsuHGCNVcPMu7kKsadXEUlryK0LlCV5gUq4+fsebv+jlVQ3lUQ+7sQPR+Jno2EDgSTP7j3BddORhDRNC3H0cNHOZREz0PCP0J5fY9ybX/X8qcjrtJr60+MqtqTlgWqZGrdrsSEsOrKAVZdOcDJiCs4oKiV2592BavTLH9FXExJ02SIxBs7pqOmgPUoOPgak9Juz6IcfFN5lodLdLSxx8TNLWs2HmravUrv8JEOCjmUiA0J6pKwEmkFyiHtF5vdQWcYtGsSv9buR+3c/llUSzgTeY1VVw6w8vJ+LsUE425ypnmBynQoVJMq3kWT9HJEBOJ3INFTIG4jkDDk5P4yylw09SfRNO2B6TmFh5xSJsj1IRLcG4kcjcr1UZrlo2xxAMlWCmW2kh75eKV0SwaWasHekHMsvbSb1VcOsDgwgGLufnQoVIt2Bavh55LLCBDO9VHO9RHrKSRqGsQsRGLmIy7tUR6vPLST0mPHjgVg8ODB2VwTTXswuqeQw9nDR0H0dJTXDyjXp1ItFxYfTfsNX9GuUA1GVOyUdRVMQbQ1jrVXD7E4MID9oecxKQca561A16J1qemb9LAfsV1LCA5zQWLAubURHBwrZF8D7oNefaTldLqn8IhQnsMRy0Ek/ANwLJvqO2kvJzdaFqjCisv7GFK2DR7mB59svl9uZmc6FK5Jh8I1OR95g78Cd7Hk0m7WXTtECfc8dClaj3YFq+Ph6IIy5UPlehfx6I9EzYDoGUjQCsS5Ccp9MMqpWra1Q9MeR3p9YA6nlBPKewwoVyTkNcQemWrZLkXqEmOLZ8XlfVlWv7sp5pGH18u1Y2mT9/i4Uhdczc58d3QJT234ii8P/8WpCCNFhnLwxcFzGCrPBpTHGxC/Dwnujj34RST+8ekZalp208NHDwmJ24GE9AXnVijvH1NcpioivLB9LHF2C3MbDr1ruor0sFhtzFyzmy2HznItJBIfT1eK5/OhRP7clCzgS4n8vhTO443ZlP7nOhIWyMIL/7Lqyn7i7Fbq5i5FnxJPUju3f2K7xB4FMb8jUZPAHmScDuf+Gsq57gO3KTPo4SMtp9Orjx5BEjUZifga5TEE5TEkxTIrL+/nowPzeKfC03QtWu+Bnu/wuav8b+YqTl0OonKJAhT28yIkMpqzV4O5FnK7x2I2OVA8nw81yxShQYVi1CpTBFdnx7vePyw+mr8Cd/H7+W0ExUVQNldBniv+BM3zV8LsYBzuIxID0fOQqIlgvwGOdYz257DgoIOCltPpoPAIEhEk7F2I/cvoLbi0TbHMawFTOBIWyPxGb5DHJdc9P09MvIXxS7Yze+0e/LzcGdGrGY2rJF3mGhkTx7lrIZy9EsTZq8EcD7zB3pOXiLVYcTSbqFGqEA0qFqdBhWKULJA7zQ148XYryy/vY9bZTZyPukkBVx96F2/I04Vq4Wp2SmhXrLERLmoC2K8bwcHzdZRT7Xtun6Y9jnRQeESJxCPBz4PlCCr3bJRj5WRlLkYF0XvrjzTMU5avqj97T/fff/oyH09fycUboXR5ojJDOz+Bp2v6dh/HWazsPXWJ7UfOs+3wOU5fCQIgv48nzaqXomXNMlQuXgAHh5QDhF3sbL5+jJlnN3Eg9AJejm70KFafbkXr4+XkltD+W8HhN6Pn4FQf5fE6yqnGPbVT0x43Oig8wsQWZGxsw4rKvQBlyp+szLTTGxh7chXfVn+OxvnSt7zTbhfavD8RJ7OJT/q0onbZIg9Uz6vBEWw/co5NB8+w7ch5LFYb+X08aVGjNC1rlqFS8fyp9iD2h5xn+pmNbLlxDFeTE52L1KF38YbkdfECbgWHuQk9hyBwegLlMRTlVPWB6ny/vvvuOwDefvvtbHl+TbubHBEUlFJtgB8BEzBJRL76z/WiwHTAO6HMeyKyLK176qBgEMsxJLgXmAqjfGegHHySXLfabbywfSw3YsOYVv9VCrr5pHKn226ERtJ6xESGdmpE39YZOywTERPHpgNnWL37BNuPGgGiUO5ctK5djk4NKlI4j3eKjzsVcZUZZzax+uoBHFC0L1SD50s2TkzGJ/ZoiJmDRE4ACTXOdfB4HeWYvgyxGUXPKWg5XbYHBaWUCTgBtAQCgV1ALxE5ckeZCcBeERmnlKoALBOR4mndVweF2yRuKxIyEMylUb7TUQ5J5w/OR97gxR3jyO/izaR6A3Ezpz0MJCJ0+3wGHq7OTBveM9PqHREdy4b9p1kRcJx/j15AEOqVL0aXJ6rwZOWSKa5kuhQdzKyzm1lyaTc2sdOqQBVeLNmE4h55jbrbIyF6JhI1GSQcXJ4ygoO5eKa14046KGg5XXqDQmbuU6gDnBKRMyISD/wOdPxPGQFuvZJ5AZczsT6PHOXcEOXzC1hPICH9ku1hKOaRhy+q9eJM5DU+PjAfu9jTvp9SPF2/IgfOXOHs1eBMq7enmwsd6lfk1yHP8M+ofgxoV4/Tl4N4+7cltPtgEuOWbONKcHiSxxRy8+Xdih3568m36VmsAeuvHabHlh/5YN9cTkdcRTl4oDwGofKsBfdXIG4dcrMt9rAP9XGhmnYPMrOn0BVoIyIvJ3zdB6grIq/dUaYAsArwAdyBFiKyO4V7DQAGABQtWrTm+fPnM6XODyuJXY2EDgXHaiifycmS5/1+bivfH/uHviWbMLhMqzTvFRQeRdsRk+jdvDrDnnkyM6udhNVmZ+uhsyzYfIBtR86hUDSqVILezapTu2yRZHMPIfGRzD23lfnntxNti6dJvgr0829G2VwFARDbTSRqHET/DigjI6vHwEzLyqp7ClpOlxOGj9ITFN5MqMNopVR9YDJQSST1t7R6+ChlErMMCXsTnOqgfH5DKdfb10T44vCf/B0YwKeVu9GuUPU07/XW+MXsP3OFJZ+/lK79BhntclAYi7Yc4s8tBwmJjKFM4Tw817wGrWuVxdFsSlI2LD6aeee38fv5bURaY2mUpxwDSjWnnFchAMQaiET9AjF/gXJFufcH9xeT/HwyQtu2xvLg5cuXZ+h9NS2j5ISgUB/4VERaJ3w9AkBEvryjzGGMwHEx4eszQD0RuZ7afXVQSJ3E/G2c7+xULyEw3M5/ZLFbGRowlX0h5/mh5gvU8yud6n32nrpEv9Hz6dW0GsO7N82KqqcozmJl+c5jzF67h9NXgsjj5U7PptXo0qgKudyT5naKsMQw/8J25p7bSrglhmb5KjGwdAtK3JpzsJ5CIn6AuNXgkM9IpeHa0chGq2mPgZwQFMwYE83NgUsYE829ReTwHWWWA/NEZJpSqjywFigkaVRKB4W0ScyfSNh7xvp9n/FJAkOkJZYBOydwOTqY8XX6J76bTsnX89Yzb8M+Jr7ZjZqlC2dF1VMlImw7cp7Za3ez4+gFXJ0d6dSgIs82r0nB3Ekn1yMtscw5t4U557YQa7PQrlB1XvZvnrj6SuIDkIivjHOkzeVRnu+inBtkR7M0LUtle1BIqEQ7YAzGctMpIjJKKfUZECAiixNWHE0EPDAmnd8RkVVp3VMHhbuTmEVI2AhwaoDyGZckMNyIDaffjvHE2y1MrPsKRdxzp3iP6Nh4eoyciXJQzPugT7YMI6XkROANZq3ZzYpdxxGEVrXK8kLLWpQpnCdJuZD4SKaf2cSCCzuwi9C5SG1e9G+Kn7OncdhP7DIkcjTYAo2zoz3eQTmm3nu6m88//xyAjz5K+9wLTcsuOSIoZAYdFNJHohci4e8n9BjGJRlDPxd5nf7/TsDd7MykugPxSyUVRsCJiwz4YUG2DyOl5FpIBHPW7WXh5gNEx1moX6EYL7SslWxS+lpsGFNOrWPxpd04OpjoWawhz5d8Eg+zi3FMaPRMJHIsSBS4djeWsZpSDpRp0RPNWk6ng4J2u8fgWNOYY3DwTLx2OPQig3dNJr+LF+Pq9MfX2SPFe9waRho9sANNq+W8U9HCo2JZsPkAc9fvJSg8mkrF89OvbR2erJz0MJ+LUUH8dnI1q64ewNvRjZdLNadzkdo4OpgRewgS+QtEzwXlgvIYDG7Po5RTGs+clA4KWk6ng4IGgMT8g4QNB3NZlO/kJEsy9wSf4fWA6RR282VcnZfxdnJP9vg4i5V+o+dz/loIM97tRYn8Gbek88j24wSs3E/eon4U8M9HQf/85C7og4PDvW+fibNYWbL9MNNXBXApKJxyRfLSv11dmlT1TxIcjoZd4qfjy9kdfIYibrkZXKYVzfJVQimFWM8gEV9D3HowFUF5vgvOLdNM5neLDgpaTqeDgpZI4jYgIUMSUmJMTZIraWfQKd7aPYOi7nkYW7tfYuK5O10NjuDZL2fj5e7CjHd74ZHOBHmpiYmMYcr7c/n71xX89+/PycWRAiXzUcA/HyUrF6N2m2qUr1cGkzl9q4QsNhvLdx5j0vJ/CbwRRpnCeejfri5Nq5ZKTMQnImy7eYKfjy/nTOR1KnkVYUjZNlT3LWFcj9uCRHwJ1pNGNtZc79/1eFAdFLScTgcFLQmJ34WEDAAHH5TPVJS5WOK1HTdP8vaemZT0yMsvtfuRyzH5Gv7dJwJ55ccFNKxYgu9feTrVTKd3s2vlPsYM/I0bF4Po+GobXvisBxHBkVw+fZXLp69x+dRVrpwxPj9/JBC7zY6Htzs1W1Whdpvq1G5TDd/8d8/jZLXZWbHLCA4XrodSqpAf/dvWpXn10ol1t4mdfy7t4beTa7gRF07jvBUYUrYNRd39ELFCzB9IxBgjp5JrN5THG6nON3Tp0gWAhQsX3tfPRdMymw4KWjJiOYQEvwTKEeUzMcm73603jvPOnln4e+bnp1p9UxxK+n3DPr6Zt55nm9fgzS5PpmtY5U7r5mzmy+d+oki5Qrw58RUqNSyXZvnI0Ch2rz7AruV72bVyH8FXQgAoXaMEjbs3pM1LTfHyS/u8CKvNzqrdx5m07F/OXQuhbOE8DOnciPrliyXWP9YWz+xzW5hxZhPxdivditajn38zvJzcEHt4wnzDrIT5hiHg9hxK5YzVWJqWXjooaCkS6ykkuB9IOMp7DMq5ceK1rTeO897e2RRw9eGnWi+S39U76WNF+Hb+Bn7fsI/+7eoyqMO9re/f+Md2Rvb4nq9Xf0yN5snPgUiz3iKc3n+OXcv3sWNpAEe2n8DR2ZHG3evz9OA2lKtTKs0gZbPbWRlwnHGLt3EpKJzaZYswtFMjKha/PZR2My6CCSfXsDgwAA9HF172b0bXovUwO5gQ62kk/AuI3wwmf1SuD1HODe+pDZqWnXRQ0FIltqvGUJL1JCrXpyi3HonX9gaf5a09M3EzO/NzrRcTdwTfYrcLn81azeLthxnauRF9W6U/xXZMZAzd8r1Mm5ea8drP/R6oDecOX2TJuJWsmbmJ6IgYStcoQYdBrWnaqxEubqnPeVisNhZuPsDEZf8SEhlDyxqlebVjQ4rmvT0kdTLiCj8eW87OoFMUdfNjaLm2PJEnoVcTtwGJGAW2C+DcAuX5HspclBEjRgDw5ZdfpvS0mpbtdFDQ0iT2SCT0deOdr/srxnh5wjvtE+FXeH33NCx2K2Nq9qWSd9LDdmx2Ox9OXcHKgOO826MpPZpUS/fzftrlW47uOMnci+Pva5XRf0VHxLB29maWjFvJ2YMX8PB2p22/ZnQa0pa8RfOk+rio2HhmrtnNzDW7ibdY6dywMgOeqoeflzFsJiJsvXGcH48v43zUTerkLsUb5drh75nf2N8QNdVIuCdWcO9Psw4LAQc90azlWDooaHclYkHCP4WYP8DlaZTXF4lr8y9FBzMkYAo34yL4utqz1M9TJsljLTYb7078hw37T/Nxn5Z0alApXc+5dvZmvurzEz9s/vyucwr31hbh8NZj/P3rCjYt2AFA4+716fpmB8rU9E/1cUHhUUxavpOFmw7g6GiiT4uavNCyVuIObqvdxsIL/zLh1BqirHE8U7QOA0q1xNvJzehxRXwLsUto1uU6mIqyfuP+e55r0bSsoIOCli4iAlHjkcgfwLEWyufXxFPcbsZFMCxgGqcjr/FhpWd4qlDSc5DjLVbeHL+EbUfOMaJnM7o1vvtRmFHh0fQu+gr+VYvz3fpPM6S38F/XL9zgz5+Ws2zSGqLDY6jevDLd3n6aWq2qpvqCfeF6KL/8vYU1e07i5+XO4Kcb0KFeBUwJ9QuNj2biqTUsvPAvHo4uDCzVgs5F6hjzDfG7aNq0DUgM65b2QuX6KMsO99G09NJBQbsnErPUSKRnKojymZD4ohZpjeWdPbMICD7Dq2Va83yJpKuO4ixW3pm4lM0Hz/J65yd4odVd/+ZYOW093700lldGv0CXN9pnVpOICovinwlrWPTjPwRdDqFk1WL0eKcTjbvXx2RKed/D/tOX+X7hJg6evULpQn682eVJ6pa/vXz3VMRVRh9dyu7gM5TyyM9b5dtTM3dJmjRpDLbrrFvgChIH7v1Q7q8kO9tC07KLDgraPZP43UjIIACUzy8opzoAxNutfHZgAauuHqBb0Xq8Wb49JnX7Hb7FauPDaStYvfsE/drUYfDTDdIcQhERPu74NXvWHGD83m8pUjb1bK0ZwRJvYd2cLcz/9m8uHL1EwVL56fluJ1r0eRJHp+RLS0WE1XtO8POfW7gUFE6jSiV4o8uTibu5RYT11w7z47FlXIkNpWX+yhz/+i9cTI7MnP4DEvENxP4NDoWMVUouzTO1fZqWHjooaPdFrBeMlUm2i6hcn6HcjE1ZdrHz8/EVzD63hcZ5K/B51e64mG7nBrLZ7XwxZy1/bj1EzybVeLtbkzQ3uAVdCaF/pTcoVKYgYzZ/nu4dyw/Cbrez9a9dzP1iISf3nCVP4dx0e/tp2r7cPMUVS/EWK79v2Mek5TuJiYun25NVGfBUPbw9jM19sTYLs85uYvqZjSil6OffjF7FG+LkYDY2C4b/D6wnwLkZyvMDlLlIsufQtKyig4J238QebhzvGb8N3F9GebyVeBjNvHPb+P7YP1T0Ksx3NfokSaQnIoxZtJmZa3bTrk45Pnm+FY6pDNMAbJi3lVG9xtDhlVYM+fXlLJugFRECVu1n7heLOLj5KN55vej2VgeefrVNisEhJCKacUu2s2jLQTxcnXilQwO6PlEFs8noLV2ODmHM8X/YcO0IRdxy83b5DtTPUwYRC0TPQCJ/BrEbifbcX7qnRHuallF0UNAeiLEyaRTEzAHnpiiv0SgHIwCsv3aYj/fPJ7ezBz/UfCHJXgYRYcqKXfy6eCsNKxbnm/7t0zyLYeI7M5n/3WI6D23HoB/6ZvnKnYObjzJr5AL2rD6Ab35ven/QhXb9m6c4rHTy0k2++2MDu45fxL9Abt7u3oS65YoCMGzYMK7EhGJ7tgoXom/SJF8F3izXnvyu3sYqpfBRELfS2Pjm9b/EoTlNyyo6KGgZQqJnI+EjweyP8h6XOARyOPQib+2ZgcVu4+vqz1Ird9Jln39uOcioOWupWDw/Pw7umDjkkuz+Iox7Yxp//rSMrm92YMC3fbJlSeehLUeZ8sFcDm4+Sv7ieXj+0x40e7ZRsglpEWH9/tP8sGAjl4LCaVqtFG91fZLeXToCsGrdGuac3cLk0+tRCvr5N6N38YZGiu64DUj4Z8bBPi6dUbneTZK1VtMykw4KWoaRuK3GRjdMCRPQxi7my9EhvLFnOhejgni/Umfa/2fJ6vp9pxgxeRmF/Lz4ZUhnCvimnKdIRPhlyGQWj11Jz3c78dIXvbMlMNwaVpr6wRxO7jlL0fKF6Pt5Lxp1rpOsPnEWK7PW7GHyin8REYLWTyO/ryebNm4E4EpMCKOPLmXT9aOU8MjLuxWepoZvSURijEN9oiaD8kB5vgeunfXeBi3T6aCgZSixnkVCXkmYgP4Y5dYTgAhLDO/tm8OuoNO8WLIJA0u3wOGOlUm7TwTyxri/cXNx4qdXOyU7NvMWu93OT4Mn8c+E1fR+/xn6ft4z214oRYQti/5l6ke/c/HYJcrU8ueV0S9Q+YnyycpeDY5gzKJN/Prx6zg5mli4eBmNq9w+4Gfz9aN8d3QpV2JCaFewOq+Xa4uPkwdiOYGEfwyWPeBUF5XrfyhzyaxuqvYY0UFBy3DGBPQbRmoM157GckvlhNVu4+sjf/N3YACN8pTjsyrd8XC8fS70yUs3GfLLn0TGxPHVy0/RqFKJFO9vt9sZM3ACyyevpfPQdrzy/QuZsrktvWw2G2tnbWbaR79zIzCIxt3r0//rPuQrljyw1axbn4vXQyncZgCNKpVgePcmFMnjDRhZWKec3sCss5txMzvxWpk2PF24JgoS0nN/CxKbMBH9sp6I1jKFDgpaphCxGbufoyaAYw2U908oU15EhIUX/2X00aUUdvPlu+p9KOZx+8Xzemgkw8b+zYnAG7zdvQk9U8mXZLfbGffGNP76eTmt+zbljQkDs2S5alpio+P449vFzPvmL0SEbm89TY/3OuHqfjvwDRgwALsITXoOYvySbVhtdvq2rs2LrWvj7GgG4Gzkdb46/BcA4+v0T+xNiO2GkWQvdhmYS6FyjUQcqyXpcWnag9JBQctUErMMCR8ByhPl/SvKyUhxsSf4DCP2zTU2vFXpwRN5b+c3io6N54Opy9l44Aw9m1Tjza6NE5d1Jrm3CDM+nc+szxfQsFNt3p8zDCeX7H/3fP3iTSa9N4v1c7eSu6APL3/1HM16N0rWm7kRGsn3CzexMuA4RfJ4817PZtSvYOyKFhEirLEpHmQkseuNXFT2q0y82R2LuQ6vleuQBS3THgc6KGiZTizHkNDBYLtmjIm7dQXgakwow/fO4kT4FQaWbsGLJZskviu22e38uGgzs9buoVGlEnzZrx3uqbzg//nTMsYOm0qVxhX47K93cPdKfvBPdji09Rhjh03l5O4zlK9Xmld/6kfZWsmT7v179Dxf/r6OC9dDaVmzDG93bUweb48U7nib2CORyDFcDV3Ab1drccFSik+qvkTx/6Qw17R7pYOCliXEHoKEvgnxW8G1e8I8gwuxtnhGHfqTlVf20zhvBT6p3DXJPMOCTfv5et56iuXzYfTApymWL+UjNtfN2cw3fX+lWIXCjPpnBH6FUj4OM6vZ7XZWz9jI5BGzCb0eTmTZK5SpXYpp06cmKRdnsTJ9VQBTVuzE0WzitY4N6fpklcREe6mR+H1I+AeMu+CMr6s/PcuO0MtXtQeS3qCgBy21B6IcfFA+E8F9IMTMR4J6INbzuJic+KxKd4aVa8eWG8d4YfuvnIy4kvi4rk9W5ZchzxAcEUOfr+awYf/pFO/frPcTjFw6gitnrjG0/gecPXQhq5qWJgcHB1r3bcrUYz/SeWg7jh09zj9zV7BiyjrsdntiOWdHMwOeqsf8j56ncon8fD1vPS9+O48TgTdSvff2Gyf46OgxrrlNZUdUJdzsh5AbbY2khQ/Zmzjt4aN7ClqGkdj1SNg7gBXl9RXKpTUA+0LO8f6+uURYYnm3Ysck+xmuBIczfMJSjpy/xktt6jCoQ/0U30Wf2neWD576krjoOD5dNJxqTdN3fkNWqVe7HoEnr1A+vA4V6pdhyC8vU6p60lVWIsKKXcf57o8NhEfH8nzLWvRvVw8XJzN2sSdOLFvtNkbsm0u0LY6aviV5poAfuWL+B5YDRh6lXP9DmfJlRzO1h1im9hSUUsvv53Hao025NEX5/QVmfyR0CPbwUYjEU82nODMbvEYl7yJ8dnABXxz6kzibBYACvrmY/FZ3OjesxJQVOxnyy1+ERsYku3epaiX4efso/Ar5MqLNSNbO3pzFrUubi7sLpaqVYPjUV7l86iqv1n6XX4dOISo8OrGMUoq2dcqx6NO+PFW3AlNX7qLHyJnsPHaBa7FhWO02AMwOJvqXaoa3oxsv+TfF260yynceyvNdiNuK3GyHRC/UvQYtU6TaU1BK1UjxAihgqYgUyLRapUH3FHI+kXgjfXT0DHCsivIegzIVwmq3MeHUGqad2Ui5XAX5slpvCrndHif/c8tBvpq3Hr9c7nw74CkqFMuf7N4RIZH8r8t37N9wmOc+6kqfT7pl616GW5o0aQLAhg0biAyNYsoHc1k6fhW+Bbx59ceXaPRM3WSb8f49doFRc9YQeCOUqs18cCgUx+fVelDQzYffz20l3BJD/1LNkzxOrOeQsPfBEgBOT6C8RqJM2fJfUXvIPPBEs1LKBmzECAL/VU9EUk5mk8l0UHh4SOwK4wUMB5TXlyiXloCxy/fTA38gwEeVu9A0X8XExxw+d5XhE5cSFB7N8G6N6fJElWQvppZ4Cz++MpGV09bToGNt3pn+Gu65svcwm2HDhgEwZsyYxO8d33WKHwb+xul952jQsTav/dyPPIWTTpTHxluZ8M92Zq7ZjVtZK57l7VT0K4RN7HQvVp86uUsllo20xuJmcjL+Q0bPQiJHAyaU5whw7apTZWhpyoigcAjoLCInU7h2UUSyJTm8DgoPF7FeQEKHgfUQuD2P8nwHpZy4HB3C+/vnciQskB7F6jOkbFucHIxNXqGRMXw4dQXbjpyjTa2yfPBsi2TLVkWEv39Zwbg3p1G4TAH+9+c7FC5TMBtamDab1caiMf8w/ZN5mMwmXhzViw6DWiVLtHf0wjU+m7ma49evUbNWfj5p35bCPrd7USLCG7unE2+38mGlLhR08zF+tmHvg2VnQq9hFMqUvHelaZAxQaErcFBEjqdwrZOI/PXAtbwPOig8fIzhpO8gehqYKxrDSeZiWOxWfj6+kt/Pb6VcrkJ8Wa1X4nCS3S5MXbmLcUu2USSPN98MaE/pQn7J7r1v/SE+7/49NquN9+cMo07b6lncuvS5cvYaPw6ayO5V+ylfrzRvTHiFEpWKJiljsdmYvjKACct24OHqzHs9m9KyRhmUUkYQDAzgx2PLEIQhZdvyTJE6gED0HCTyW8CMyvWBkYFV9xq0/9D7FLQcR2LXGudAY0Xl+hzlapzPvOHaET4/uAA7woeVnqF5/sqJjwk4cZH3Jy8jMiaed3o2pWP9isle8K6eu86nz3zLmf3neWlUL3q82ynLXxSfe+45AGbNmpVqGRFh3ZwtjHtjKlFh0fQa8Qy93u+c7OyGU5du8unMVRw5f41m1UoxolczcucyNu5djQll5KFF7Aw6Re3c/nxY6RkKuPog1vPGz9ayO2GF0ucoU8rJB7XHU4YGBaXUU0BFIHH3kYh89kA1vE86KDzcxHbZ2Oxm2QOuXVCeH6Ec3LgcHcIH++dyOCyQLkXq8nq5driYjBfLoPAoPpiynJ3HL9KuTjlG9GqebDgpNjqO0S+PY8PvW6n7VA2GT30VL7+UU3Vnhjsnmu8m7GY4496YxtrZmylRuSjDp75K6RpJM6RabXZmrglg/NIduDs78l6vZrSqWRYwgsufF3fy0/HlKBRvlH+KDoVqAnbjpLeI70G5Gof5uLTN4JZqD6sMCwpKqfGAG9AUmAR0BXaKSL+MqOi90kHh4SdiRSJ/gahxYCqO8v4B5VgBi93KuJOrmXV2M6U88jOqWs/EU91sdjtTVuzkt6U7KJzHi69ffoqyRfL+577GPMOE4TPI5efJiFmvU7VJxZSqkOHuJSjcsn1JAGNemUDo9TB6vdeZ3h92wek/p9SduRLEJ9NXcvj8NVrWLMN7PZvhk3Bg0aXoYD4/uJA9IWdplKcc71fqjJ+zJ2I9bewXsRwEl/ZGqnMH7wxqqfawysigcEBEqtzx0QNYLiJPZFRl74UOCo8OiduOhA0HewjK821wewGlHNh+4wSfHvyDaGs8b5VvT8fCtRKHg3afCOT9KcsIjYplWOcn6Nm0WrKholN7zzKq1w9cOnmVZz/swnMfdc30TKv3ExTAWGI77o1prJ6xkeKVijB86quUqZk0j5LVZmfaql1M+GcH3u4ufPRcS56obPQs7GJn3vntjD2xEheTEyMqdqJZ/kqIWCFqghF8HXyN1V/O2fJfVsshMjIo/CsidZVSO4BngCDgsIiUSvOBmUQHhUeL2IORsA8gbm3CCpqvUSY/bsaG8+nBBewMOkWzfJV4v1LnxMyiIZEx/G/GKjYdPMMTlUvw6fOtE9893xITGcPPQyazevpGKj9ZnhGzXk+2HDQj3W9QuGXH0t38MPA3o9cwojPPftgl2VzDicAbfDhtBacu3aRzo0q82aVx4jDa2cjrfHLgD46FX6JtweoML98BD0cXxHLYCLzWU+D2bMLqr2xZTa5ls4wMCh8BPwPNgV8BASaJyEcZUdF7pYPCo0dEIGYuEv4lOHgkvKttgl3szD63hbEnVuHn7Mn/qnSnhm+JxMfM27CPHxZtxsvdhc9faE3d8sWS3Xv1zI38NHgijs6OvDVpEA071cmUNowYMQKAL7/88r7vERESydhhUzm15yy/7voqxXTh8RYr45duZ/ryfymY14fPX2xDNf9CgJEeY8rp9Uw9swE/Z08+qdyVWrn9EYk15hmip4GpBMr7O5Rj5WT31h5tGRkUnEUk7tbnGJPNsbe+l9V0UHh0ieUEEvYWWI+DWx+U53CUcuFIWCAf7Z9HYHQwL5RszIBSzTE7GMNBxy9eZ8SU5Zy7GkyfFjV59ekGOCUcanNL4MkrfNF7DCd3n6Ftv+YM+uEFXD1y7rvlqPDoNDfjRUfE8F6nrzlpsnGjVG5ebF2bAe3r4Ziw9+Fw6EU+OfAHF6ODeLZ4I14p0xInB3PCcN27YL+J8hgC7v1Rypzq82iPlowMCntEpMbdvpdVPD09pWbNmkm+1717dwYPHkx0dDTt2rVL9pi+ffvSt29fbt68SdeuXZNdHzRoED169ODixYv06dMn2fW33nqLDh06cPz4cQYOHJjs+ocffkiLFi3Yt29f4s7WO33xxRc0aNCAbdu28f777ye7PmbMGKpVq8aaNWsYOXJksuu//fYbZcuWZcmSJYwePTrZ9ZkzZ1KkSBHmzZvHuHHjkl1fsGABfn5+TJs2jWnTpiW7vmzZMtzc3Bg7dizz589Pdv3WkMh3333H0qVLk1xzdXVl+XIjFdbnn3/O2rVrk1zPnTs3CxcuBIx309u3b09yvXDhwonLOIcNG8a+fXvAFgi2q6BcKV2+GRMn/UG0NY7GPdtz/MQJ3MzOlHDPg7PJkWrVqvHlN98yZuEmvv3kHcyWaEoW8MUlYeilfv36fPnll1jiLdSr3JCzJ87h5OJI0fKFcMvlRvPmzfnoI6PT27ZtW2JikuZdat++PW+//TZwe4joTtnytxevUJFmTIe88OlQktMFC5DfHEPw7mW4OBkv8nYRAqODcO5Vgxp1a9MppgS/jvwOsIL1HNiDQXnww4/TqV6rrf7bGzaMffv2JblepkwZJkyYABgn6504cSLJ9WrVqiXuYH/uuecIDAxMcv3W3x5Aly5dCAoKSnI9q//2Nm7c+GAJ8ZRS+ZVSNQFXpVR1pVSNhH9NMFYj3ZVSqo1S6rhS6pRS6r1UynRXSh1RSh1WSs1Jz321R5kDmIqCuSxghdg1SNQkXE2OVPctQUmPvMTZLBwNv8TNuAhEBFcnR0b0ak7N0oWxWG0cuXCd66GRSe7q6ORI8UpF8K9WHBE4tfccV89ex26zp1yNe9S7d28OHz6cIfe6KycBiwOSy0r/wS34dkB7rodFcvTCNW6GRQHgoBRF3f0YUrYNN+PCGXlwITfiwgEzmEuB2R8kBgkbhsT8pZPraYnS2tH8AtAXqAXcOV4TDkwXkUVp3lgpE3ACaAkEAruAXiJy5I4ypYH5QDMRCVFK5RWR62ndVw8fPT6MSeiPIG41ONVNmIQuyPXYMD47uJCdQad4Ik85Pqj0DL7OxolmN8Oi+N/MVWw9fI76FYrxaZ9WyU47iwqP5tfXp7B6+kZK1yzJezOHUrRcoQeq64NONN+L4wGnGd1vLMN+G0iFemWw2WwERcTw8bQV7Dx+kabVSvHxcy3xSjhD+mZcBJ8fXMj2mydomKcsH1Xqgq+zB2K7hIQON5LruTxlpOR2yLq9HVrWysjhoy4isvA+KlAf+FREWid8PQJARL68o8w3wAkRmZTe++qg8HgxJqEXIREjAYXK9Qm4PI0gzDu/nV9PrMTd7Mz7FTvTOF+FxMcs2HSAHxZuwsnRxIhezWldq2yye29euIMfBv5GXHQc/b54lk5D2953xtWsCgrhQRGM6j2GGs2r0OOdjohI4pJcu12YvW4PP/+1hdy53BjZty01yxQGjJ/J/Avb+fn4CjwdXfm0clfq+pVGxJawdPUncMhvTEI71UyrCtpDKiPPU9iqlJp86wwFpVQFpVR6Nq4VAi7e8XVgwvfuVAYoo5TaqpTaoZRqk9KNlFIDlFIBSqmAGzdSP7FKe/QopVBuXVC5/wZzGSRsOBL6OkrC6FW8IdPrDyaPcy6G753FZwcXEGmNRSlFt8ZVmfvBcxTN68OIycsYMXkZYVGxSe79RJd6TDz4PdWbV2bcm9N4u+mnXDp1JZWaZK9DW48BxrnVzq5O9HinY7IyDg6KPi1qMv2dnjg7mhk4ZgHjlmzDarOjlKJHsQZMqz+YXI6uDAmYyk/Hl2MVQXkMQvn+DsoBCX4WifzZ2OegPZbS01NYDkwFPhCRqspYrrBXRNJc05aQUK+NiLyc8HUfoK6IvHZHmaWABegOFAY2AZVFJDS1++qewuPLeFc7GYn8ERy8UblGolyaYrFbmXRqHdPPbCSfixcfV+5KzdzG5q7EjV9Ld+Dt6crHz7WkUaXkJ6KtnrGRscOmYo238tKo3vfca8jMnkJMZAyfdP6WsBvh2G12Jh78HgCbzZYs2+otUbHxfDNvPUt2HKGaf0G+eKkd+X09AYi1xTPm2DIWXdxJ+VyFGFWtJ4XdciP2SCT8U4hdDI41jV6D6cGG1bScIyN7Cn4iMh+wA4jxFsKWjsddAu5Mr1044Xt3CgQWi4hFRM5izEGUTse9tceQUiaUxwBU7gXg4IuEDsQe9j5mYhlUphUT6g7E7GBi0K5JfH90KbE2C2aTAy+3rcuM93rh5ebC0F//4rOZq4mIibvjvopWLzRh0qHvqdasEuPenMZbTT4h8GT6ew3NmzenefPmmdFsXD1c+Wb1x1RrWokrZ65x5cw1gFQDAoC7ixP/e6E1I19sw4nAG/T6YhabDp4BwMXkxHsVO/F19WcJjA6iz9ZfWHVlP8rBAwfv71Be34L1GHKzIxK7MlPapOVc6ekpbAC6AKtFpIZSqh7wtYg0vsvjzBgv8s0xgsEuoLeIHL6jTBuMyecXlFJ+wF6gmogEpXRP0D0FzSASb4yDR00ChwIJG97qEWON55cTK/jjwg6KuvnxceWuVPExUlTHW6z89s8Opq8KIK+3Bx8+14IGFYr/577CmpmbGDtsKvGx8fT9rCfPDHsq09NkpNe//+wmMjSaJj0aJKvT2YPnKV6paLK0H+evhfDepH84HniDPi1q8lqnhol7Gq7EhPDR/nkcCL3A04Vr8Xb59riYnJKeg+HaG5VrBMY2Je1hlZETzTUwdjRXAg4BeYCuInIgHZVoB4wBTMAUERmllPoMCBCRxcr46x0NtMHofYwSkd/TuqcOCtqdJH6vsSHLdg7cnkN5vI1ycGPnzVOMPLSI67Fh9C7RiIGlWuCckHX14NkrfDpjFWevBtOpQSXe6Poknq5JX/BuXg7mp8ET2b44gDK1/Hlz4iv4Vy2e9Q1MwZ2Ty7ecPxrIoOrDeaJrPd6Y8AoubknbE2ex8v2CTfyxaT9VShbgq35PJQ4n3XlMqr9HPr6o1osSHnkTzsH4AaIng7lcwjkYSbO5ag+PjE6dbQbKYhzNeVxELA9exfujg4L2XyIxCWkcpoOpmNFrcKpFpDWWn48t58/AXRR3z8PHlbtSydsY0YyzWJmQ0Gvw83Lnkz6tqF+h2H/uK2ycv41fX59KeFAE3d7qQJ9PuuHsmvwdc9u2RorqW5upspqIMPfLP5n20e/4VyvOp4uGk69Y8vMUVgYcZ+TsNZhNDnzet02S+ZVbiQhjbPG8W6EjTxUy9qdK3AYk9B0gHpXrU5RrpyxqlZaRMjooNACKA4l74kVkxoNU8H7poKClRuJ3GgfN2C4ZGVc930ApV3bcPMmoQ4u4ERtO7xKNGFCqReJZDUfOX+XTGat4uV3dxPMK/is8OIKJw2eyfMYGChXPy+vjBlCjRZUkZbJyn0Ja/v1nN188+yOOTmY+mv9WiqnDz18L4Z2JSzl56SYvtq7NoA4NMJuM6cUbseF8tH8ee0LO0qFQTYZX6GAMJ9muIqFvgWUXuD6D8vwY5ZC952Jr9yYjh49mAv7APm5PMIuIDH3QSt4PHRS0tIg9yjiaMnpOQq/hC5RTbSItsfx0fDl/Be6imLsfH1bqQlUfo2dgsdkwOzikeVpbRHQsn4//h/2rD2BfeYSWzzdmwLd98M7jBeScoABw8fglPun8LZdOXmHwmBfp+Gryld6x8Va+nb+eP7ceolaZInzVrx2+CfmWrHYbk06vY+rpDZT0yMtX1XpTzCMPIjYk8leI+hVMJVHeP6Icy2R187T7lJFB4ShQQXLIPngdFLT0kLgdSPj7Rh4lt+dQHm+hHNzZefMUow4t4mpsGN2K1mNwmVa4me8+gRoRE8eZy0Es2nyAHQEnsf8egCcOvPzVc7R5qSnNmjUDckZQAIgKi+KrPj+zY+luOg1pyyvfv5DiaqUl2w/zxdy1eLm78k3/9lQpWSDx2vYbJ/jkwHwsdhsfVr59TKqRWO8tsEcap7u5ds6ydmn3LyODwh/AUBHJEbt6dFDQ0svoNXwP0TPBVNjY1+DcgGhrHONOrmL++R3kd/FiRKXO1PNL/0roL+eupaCzM/t/XMXBzUcpX680/8avw9XDJccEBTD2MUx4eyaLfvyHhp3rMGLW0BTnQ45fvM7bvy3hWmgk73RvQpcnqiT2mq7FhDJi31wOhV2kV7GGDCnbBrODCbFdTzhWdSe4djVOd1Muye6t5RwPvE9BKbVEKbUY8AOOKKVWKqUW3/qXkZXVtMygHNxxyPURync2YEaijYyYbmZn3irfgQl1B+BscmRowFT+d2ABYfHRKd7HYr29LSc0MoYzV4IpUbogozf8j3emvcaVM9cJ3RODjyUvkaFRWdG0dDGZTAz6oS+DfujLtr928U6Lzwi7GZ6sXNkieZk14lnqlivKF3PX8enMVcTGGzua87l681vd/nQvWp+557cyaOckrseGoUx5Ub7TwP0ViFmABHVHrOeytoFapkgrIV5jjNVGXwPv3HkJY59C3cyvXnK6p6DdD5EYkFiUg0+S78fZLEw5vZ4ZZzfh5ejKm+Xb0TJ/8iM+LTYba3afZMmOwxTL58uAdnXx8TTG4CNDo5j20e8sGbeSXH65GPBNH1r0eTLNOYqstnnhDr7q8xN5ivjxxbL3KeifP1kZu1347Z/tTFz2L+WK5OW7ge0pmNsr8fqqK/sZdehPXEyOjKrak1q5jWNDb69Oshgrv1xSzFajZbPMPk/hgIhUSe0xmUkHBS0znAi/wqhDizgbeZ7WeQswrMg5LKoHfwU4cvF6KFdDIvBwceLJKv48Vbd8ivc4tfcsP706kaM7TlKhQVle++klStfIOev6D287zkdPf4XJbOLL5R9QqnqJFMttOniGD6euwGxy4PtXOiSe7AbGsZ/v7Z3N+aibDC3Xll7FGqKUQmyXkdDXwbIf3F825nBUztjwpxkeOCgopQYBg4GSwOk7LnkCW0XkuYyo6L3SQUHLLNa4PZy78Qv7Qq/TMFcghy5X470ZVXm5bV0GdWiQ5mNvrT5at24dq6ZtYPKI2YTdjKD1i015aVQvfPJ5Z34D0uHi8Uu813okkSFRfPb3uykuWQVj2errY//iSnAEn/ZpRds65RKvRVnj+OzgAtZfO0zbgtUZUbETLibHhM1uXxgrv5waorx/QDl4Z1HLtLvJiKDgBfgAXwJ3HpATISLBGVLL+6CDgpbRxHIEiV0FCMpckiCrG5dv/sgHZ0qhrpbE9YIXTSqVTjEwnLp0k2L5fGjZwsh7dGuiOSosilmfL+Svn5fh5OJE7w+60Pn1djg5O2Zhy1J2IzCI91p/zpUz1/lg7rBUz60OjYzh7QlL2HPyEgOfqseAp+rdTtMtdqae2cBvJ9dQLlchvq3+LPlcvQGQ6D+MxHqm/CjvsSjHlPd/aFkrQzev5SQ6KGgZSSQOCXsfrGdx8FuUkCn0M3CswvqIqow+uoSgiGiKB5akf8OGtKh6+wUuOjaeDh9NwdvDlXP/jMfTzTnZ6qPAk1f47e3p7Fiym4L++Rj43QvUf7pWts83hAdF8EH7LzkRcJq3pwymZZ+UU5lZrDZGzl7Dkh1HaFunHJ881zLJGdibrx/l4/3zcTaZ+bJab6r7GkNSEr8PCX0NJALl9RXKpW2WtEtLnQ4KmpZOYg9Fwj4G20VwLAc44eD1PwAiLDGMO7mKhRd2ktvRkzcrPkXzfJUSX9Q3HTzDt/PWs2HKl/h4uhHw7zby/uekN4CAVfsZ/+Y0zh8JpHrzygz4pk+qY/pZJToihk86fc2+9YcZOrY/HV5plWI5EWHyip2MXbyNav4F+WFQx8RT3QDORV7n7T2zuBQTzPDyHXimqLEGxVi2OgQse8F9IMrjDZS6v0OMtAeng4KmpYOILXFC1B7xE0RPBZe2OHh9keTa4dCLfHn4L05EXKGeX2mGl3+aIu65AWN3cJWadbgaHEHlbq/zcts6PNusRpJ31AA2q42lv61m+ifziAiOpHH3+rzwvx4UKZt9ZxbEx8bzWbfR/PvPHt74bSDt+rdItezKgON8Mn0lBf28+OW1zhTMffvozkhLLB8dmMfWG8d5tngjhpRtg4NyMOYZwj+DmPng3Azl9R3KIXnQ1DKfDgqalk4idrBdRiJ/NrKAOpZFOTe5fd0eAcoDm9hZeOFfxp9cjUVsPF/iSV4o2RhnkyNjx44lNDKGoFylWb//NEXzejO8WxMaVkreG4gKi2LB90tZ+MNS4qLjaPl8E/p80i3FBHZZIT7OwqfPfEvAin28OWkQbV5smmrZ3ScDeWPcYtycHfl1SGf8C/olXrOJnR+O/sP8C9tpkq8Cn1XpbuRNEoHo2UjEKDD7o7zHocxFUn0OLXPooKBp90AkDqJnodz7/ef7diS4JyhPVK5PUOai3IwNZ8yxZay6eoDCbr68Vb4DDfPcnmvYduQc387bwPnrITSqVIK3ujamWD6f/z4loTfC+P2rv1g8diVit9N+YCt6vd8Z3/zJy2a2+Nh4Pur4NXvXHGT4tFdTnWMAOBF4g9d++ZN4i5UfB3eiqn/BJNfnndvG98f+obxXIb6r0Qc/ZyNFt8RtNZatYkL5/IJyqp2ZTdL+QwcFTbtPd55XIGKD6FlI5BgQK8pjMLj3Qykndt48xbdHF3M+6iYNPP0ZUrYN/n7GUJDFauP3DfuY+M8OYuOt9GpWnZfb1U12bgMYq4Fmf76A5VPW4eTsSIdBrej29tNZvow1LiaODzt8xYENh3l35lCa9WqUatlLN8N49edFXA+J5JsB7ZMdcbrp+lE+3P87Po7ufF/zBfw98wEg1rNIyCtgC0Tl+h/KrWumtkm7TQcFTctAYruKRHwJscvBVMLoNTg3wGK3MufcVl7v/AICjF44medKPJmYmjsoPIqxi7fx17ZDeLu7MujpBnRuWAlTCuc/B568wqzP/2D9nC2Yncy07dec7sOfJm/RrBtWio2O44OnvuDw1uP8tH0UZWr6p1o2ODyaIb/+ycnAm3zWtzVtapdLcv1o2CXe3DODOJuFb6s/l3huttjDjFPd4reCWz+U59t6o1sW0EFB0zKBxG02Jk5t58HlKZTneyhTPho82YhL0cHkG9mJgq4+vFm+PU/kKZfY4zh64Rrf/bGRvacuUaqQH291eZK65Yul+ByBJ68w76s/WT1zEwAtnnuSnu91onCZgimWz2gRIZEMqPIWbrlcGbf7G5xcnFItGxkTxxvjFrPnVCAjejWn6xNJEx1cjQnl9d3TCIwK4uPKXWldsCoAIlZjjiF6Njg3T5iAds/Udj3udFDQtEwiEgdRE5HI8aAcUR5DafrUNEDx3cLJfHt0CWcjr1PfrwxvlGtHcY+8CY8T1u09xZhFm7gUFE7DisUZ9swTSSZr73T9wg3++G4JyyatwRJn5clu9eg+vGOa794zyq4Ve3m/3Rd0f/tp+n/TJ82ysfFW3pm4lC2HzjK0cyP6tko6VxBuiWH4npnsDTnH62Xb8myJJxKvSdQMYxe0uRzKZzzKlDwnk5YxdFDQtEwm1gtIxOcQt5FmXYLAXJwNG3djtdv448IOJpxaQ6zNQo9iDXjZvxkejsba/niLld837GPS8p1Ex8bTuVElXmlfn9y5Un6nHHI9jEU/LGXx2JVER8RQvXllug/vSM2WVTJ1E9wPA8azfPI6ftj8ORUbpL0r2WK18dG0FazafYL+7erySvv6SeoWZ7Pw6cE/WHv1EL2KNeT1cm1xSNizYCTUewOUB8rnN5RjhUxr0+NMBwVNywIiAnHraNq8K0gc65YNQnm+gzLlIzguknEnV7E4cDc+Tm68UroVHQrXxJTwYhgSGcPEf3awYNMBnJ3MvNCqFs82r4GrU8qpMKLCovhnwhoW/fgPQZdDKFm1GF3f6ECTng1wTOUxDyI6IoaBVd/C5Ghm8pEfUjyk5042u51Rs9fy17ZDvNCqFkM7NUoSGOxi54dj/zDv/HZa5a/CJ1W64uhg7OUQy1EkZCBIOMp7TJIlwVrG0EFB07LQ1KkTkbj19O10wBhSch8M7i+glBNHwy7x/dGl7A89TxnPArxZvj01fG+v1jl/LYSf/tzM+v2nyePlzivt69OhfsXEc5P/Kz7Owro5W1j4/RLOHb6IbwEfOr7ahrYvN8cnr1eKj7lf63/fyhe9x/Dd+k+p2jjl5Hl3stuFr+et549N+xn4VD0Gtq+f5LqIMOPsJn49sZK6uUvzdfXeiSffie2aERisx1C5PkW59czQtjzudFDQtGxgDCl9CXFrE1YpfYByfhIRYfXVA/x8fAXXYsNolq8iQ8q2pZCbb+Jj9566xI9/bubAmSsUz+/Lax0b0rSqf6pDRCJCwKr9LPh+CXtWH8DRycyT3erTYVBrKtQvkyFDSzGRMXTN24+2/Zrz2s/97v4AjMDwv1mrWLL9CK93foIXWiV/HVoSuJtRhxZR3qsQP9R8AW8nY+hM7FHGXob4TeA+AOXxpk6NkUF0UNC0LHTz5k0A/PyMSWOJ24CEfwG2c0Z6B88RKHMxYm3xzDq7mRlnN2Gz2+lRrAEv+TdNnG8QETbsP83Pf2/l3NVgKpcowNDOjahZunCaz3/h2CWWjFvJqukbiA6PoWTVYjw9qDXNejfC1cP1gdr26TPfcGznKeZcGI9DCktpU2Kz2/lgynJW7T7Buz2a0qNJtWRlNl47wgf7f6egqw8/13rxdpZVsSLh/4OYecYKL6+vUSr1FVBa+uigoGlZ6NZ5CndmSRWJh+gZSOSvIPHg/iLKfRDKwZ0bseGMO7mKfy7txcvRlQGlW9CpcG3MDsa4vdVmZ+mOI4xfup3roZE0qlSC1zo2pEzhtPcsxETFsn7OFv4eu4Iz+8/jlsuVpj0aUrd9Tao1q4Sr+72fo7x29ma+6vMTP24dSYX66U+DbbHZeGfCUjYeOMMnfVrRsUHy4ac9wWd5a88MPMwu/FLrJYp5GO0TEYiagESOBsfaKJ9xKIdcyR6vpZ8OCpqWhVIKCreI7ToS8R3E/gUOeVGeb4JLJ5Ry4FjYJcYcW8aekLMUd8/DkLJtaZSnbOLQT0y8hXnr9zF15S4iYuJoUaM0L7ete9fgICIc3XGCxeNWsu2vXcRExuLoZKbyk+Wp3aY6tdtWp2i5QncdYrp8+irf9x/P/g2HeWfaa7R8PvX0FymJs1h5Y9zf7Dp+kV+HPEOdckWTlTkRfpmhAVNRyoHxdfpTzP32El2JWYqEvWvkTPKZgjKlvHxXuzsdFDQtC6UVFG6R+H1IxEiwHADHyijPD1BONRARNl0/ys/HV3Ah+iY1fUswpGxbKnjdHjIKj4pl5prd/L5hH1Gx8TStVooB7epStkjeu9YtPs7CoS3H2LV8LzuX7+HC0UsA5C+eh3J1S1OgZD4K+uengL/xMXdBH0SEP39cxrSPfsfkaGLgt8/T9uXm9zVPERkTR99vfycoPJrZI3onOff5ljOR1xi0cxJmZWJ8nf6JGWjhVs6kweCQH+U7FWXKmk18jxodFDQtC6UnKEBCRtbYv5GI0WC/nrAr+m2UqRBWu41FF3cy6dRaQi3RtCpQhUGlWyWZjA6PimX2uj3MXb+PyJg4mlT1p3+7upQvmi/ddb12/gY7l+8lYOU+zh26wNVzN7Db7InXnVwc8fB2J/hqKPU61OT1sf3xK5Q7jTve3YXrITz31VwK5c7FlOE9Ulx2ezriKq/snISLyZHxdfonabfE70ZCBhh7GXynoczZexbFw0gHBU3LQukNCreIPRqJmghRk4xvuPdDufdHObgTaY1l5plNzDm3FbvY6Vq0Hi/5N8XLyS3x8RHRscxZt5c56/YSERPHE5VL8FLrOskylqaHzWrj+oWbXD59lcunrnL59DVuXgqiQcc6NOnRIMM2yG0+eIZh4/6mTe1yjOzbJsX7ngi/zOBdk3E3OTO+bn8KuN7OGCuWI0jIS4BC+UzWm9zukQ4KmpaF5s2bB0CPHj3u6XFiu2z0GmKXGPMNHm+AayeUMnE9NowJp9ayNHA37mZnni/ZmB7F6uNiur0SJyImjnnr9zF73R7ComKpUboQL7SsRcOKJXBwyN4jP1MycdkOxi3ZzoiezejWuGqKZY6FXeLVXZPxdHRlUt2B+LncnmAW6xkkuC9IlDHH4JTyPbTkdFDQtIeIxO819jdY9hl5gDzfRTk3BIxhlV9OrGTrjePkdvakn39TOhaulbgbGCAmzsKiLQeZvXYPV0Mi8C+Qmxda1aJ1rbI4mnNOBlK7XXjtlz/Zf+Yy8z/sQyG/lDfbHQ69yOBdkynk5stvdfrj6Xh7Wa3YLiHBz4M9BOUzVQeGdNJBQdOy0MWLFwEoUuT+TxQTEYhdZizDtAWC05NGygzHMgDsDT7L2JOr2B9yngKuPgwo1Zw2Baslps0AYxnoyoDjzFgVwKnLQeTz8aB3sxp0blgJjxTOcsgOV4LD6f75TMoXzcv417um2qP59+ZJ3tg9g8reRfip1os4m27PQ4jtChLcB+zBOjCkkw4KmpaF7nVOIS3G/oaZSOQ4kEhw7YryGIoy5UVE2HHzJONOruJY+GVKeOTllVItaJKvYpIxehFh2+FzTFsVwO6Tgbi7ONG5YSV6NatOAd/sX+//55aDfD57De/1bEb3VIaRAFZd2c+H++fRLF8lRlXrmSQAGoHhObCHGquSHKukeh9NBwVNy1IZGRRuEXuIERiiZyfkU+oHbi+hHNwREdZfO8z4k6s5F3WDcrkK0b9UMxrdcYbDLUfOX2XWmj2s3nMCgBY1yvBc8xpULJ59aapFhFd/vvswEsCcc1sYc2wZ3YrW4+3yHZIGP9vlhMAQZqxKcqycFdV/KOmgoGlZKDOCwi1iPW9sfotbCQ6+KPdB4NYTpZyx2m2suLyPSafXcTkmhDKeBXjJvylN8lVITE19y5XgcH5fv49FWw4SFRtPNf+C9GhSjWbVSmXLvMOV4HA6fzKNzo0q826PpmmW/fHYMmaf28Lw8h3oVuw/SfYSA0MEyne6XpWUCh0UNC0LZWZQuEXi9yOR30P8dnAogPJ4DVw7o5TZCA5X9jP19HouRgdR0iMvL/k3pXn+ykmGXMDYTLZ4+2F+X7+PwJth+Hq68XT9CjzTqDKF83hnWv1TMmLyMrYfOcfKrwbg7GhOtZxN7LyzZxbbbp7gx1p9qZO7VJLrYg1Egp8FiUX5zkI5ls7sqj90dFDQtCyUFUHhFonbZgQHywEjE6vH6+DSBqUcsImd1VcOMPX0es5G3aCYux99SzahdYGqiXmVbrHbhe1Hz7No8wE2HTyDzS7UK1+MLk9U5skqJXG8y/kJGWHH0fMM/mkRX/ZrR+taaedVirLG8fKO8dyIC2dKvUEUdU+a8kKs543AgN0IDOaSmVjzh48OCpqWhZYsWQJAhw4dsuT5jMN91iCRY8B6EswVjODg3ASlFHaxs/7aYSafXs+piKsUcPXh+RJP0r5QjSSreG65HhrJX1sP8efWg1wLicTPy522tcvRulZZyhfNm2knvNntQvuPJlM8nw9jh3a5a/lL0cH03T4WHyd3ptQblJhd9haxnjaGkjChfOegzMlzLT2udFDQtMeAiA1ilyCRP4PtIjhWNYKDU0OUUogIW24cY+rpDRwKu0huZ096FWvIM0Xr4GFOnjHVarOz9fBZ/txyiG2Hz2G12ymSx5uWNcvQskZpyhTOk+EB4sc/NzNjdQCbv38VN5e7p8jeE3yGV3dNoXHe8nxV/dlk18Vy3Fiu6uCJyv0HysE3hbs8fnRQ0LQsdPz4cQDKlk1/aumMJGKBmEXGaiX7ZXCskRAc6iUGh4DgM0w/s5GdQafwNLvQtWg9ehRrgK+zR4r3DIuKZd3ek6zec5Jdxy9gswtF83rTokYZnqxckorF82FK5/kKqTl7NZj3Jv1DSEQ0K78akO6AM/3MRn49sZKvq/Wmaf5KyX8e8fuNoSTHqsZyVX0eQ84ICkqpNsCPgAmYJCJfpVKuC7AAqC0iab7i66Cg5URZOaeQFpF4iFmQEByugWMdlOdQlFOdxDJHwgKZfmYjG64dwcnBRPtCNeldvFGSzKT/FRIZw/p9p1i9+wS7jl/ELkIuN2fqlS9G/QrFaVChGHm8Uw4uKbHYbMxYtZsJy3bg6mTm0+db06Sqf7ofb7Xb6Lt9LCHxUcxrNCzZMBKAxCxBwt4y9nnkGpVpQ2APi2wPCkopE3ACaAkEAruAXiJy5D/lPIF/ACfgNR0UtIdRTgkKt4jEQfQ8JGqCkY3VqS7KY0iS4HAu8jqzz21h2aW92MRO0/wV6VPiySQpu1MSGhnDv8cusO3IObYfOc/NsCgAShfyo175YpQu5EeJ/L4Uz++LewrDQUcvXOOzmas5HniDFjVK826PpuTO5X7PbTwSFshL28fxTNE6vFOhY4pl7BE/QNQ44+Q79xfv+TkeJTkhKNQHPhWR1glfjwAQkS//U24MsBoYDrytg4L2MMppQeEWkViInn87ODjWMZayOtVNfOd8Mzaceee3s/Div0RaY6npW4LnSjxJfb/SyfY6JL+/cPLSTbYdOce2w+fZf+YyFqst8Xp+H09KFPCleD5fvD1cOH8thJUBx/HxcOXdns1oXv3Blo6OPrqE+ed3MLneK1TyTp5iRMSOhA6BuLUon99Qzvd2SNCjJCcEha5AGxF5OeHrPkBdEXntjjI1gA9EpItSagOpBAWl1ABgAEDRokVrnj9/PlPqrGn3K6cGhVuSB4eaCcHhdmrsSGssf1/cxdzz27geG0ZRNz+6FK3LU4VqkMsxfec8W2w2Lt0I4+zVYM5cCebs1SDOXAnm3LVgYuOt+Hi40qSqP693foJc93E06H9FWePovPFbqvkW55vqz6Xcdns0EtwTbJdRuRc9tiuScnxQUEo5AOuAviJyLq2gcCfdU9ByopweFG4xhpX+SAgOV42JWPdB4Nw0MThY7FbWXj3EHxd2cDD0As4OjrQuWJWuRepSzqvQfT2v3S5YbLY0N6jdr5+Pr2DOuS383Xg4eV1STpch1gtI0DNgKoTKPQ+lHjwgPWxyQlBIc/hIKeUFnAYiEx6SHwgGnk4rMOigoOVEa9asAaBFixbZXJP0MSakFxnBwRZopOv2GATOrTCmAw3Hwy+z8MK/rLiyj1ibhUpeRehStC7N81fGJYX9DtnhUnQwz2waTT//pgwonfrPX2LXI6EDwbULKtcXj93Ec04ICmaMiebmwCWMiebeInI4lfIb0D0FTctSIhaIXYpEjgfbWTCVRHkMBJf2KHX7RT/CEsOyS3tZcHEH56Nu4mF2oXWBqnQoXJPyuQpl+wvs6wHTOBVxlb8bD0+2c/tO9ogfIepXYzWSW7csrGH2S29QeLBFxmkQESvwGrASOArMF5HDSqnPlFJPZ9bzalp22LdvH/v27cvuatwzpRxRrp1RfstQXmNAOSJh7yI3WyNRMxF7NACejq70KN6A+Y3eYGztfjyRtxxLL+2m7/ax9N76E7PPbiE4LjLtJ8tETxeuyY24cA6GXkiz3K15FAn/DLEGZlHtHi5685qmZYCHZU7hboz0GeuQqN+MU+CUl5GR1e05lClfkrIRlhhWXTnA0ku7ORwWiEk50DBPWVoXqMoTecslOTY0s22+fpS39sxkWv3Bd11SK7YryI3W4NIUB+8fs6iG2S+9PYWMn/XRNO2hpZQCl+Yol+ZI/B4kagpETUCipiAu7VHuL6IcywFG76FL0bp0KVqXM5HXWHppDysu72PT9aO4mZxonK8CrQpUpW7uUmkO6WSE0HijR+Pl6Hb3NpoKgMcAJPInJK43yrluptbtYaODgqZpKVJONVBONYyVO9HTIGYhEvsn4tQQ5fY8ODdGJexjKOmRj6Fl2/JqmdbsDT7Lyiv7WX/tMMsv78PL0Y1m+SvRIn9lqvkUS3K2dEYJtRgb6Hyc0rkJzv1liF6IRIwEpz8xpkA10EFB07S7UOaiqFwfIx5DIfp3JHqWsYrHVBTcehureRyMpaAm5UCt3P7Uyu3POxWeZsfNk6y8vJ/ll/fy58WduJqcqOFbgtq5/amTuxT+HvkeeJLaardxJCwQJwczrukcslLKBXK9i4QOhZgF4NbzgerwKNFBQdO0dFEO3uDxCrj3g9hVSPRMJOIriBiDuHZAuT2b5NQzRwczT+QtzxN5yxNtjWNX0Gn+DTrJzpun2XrDSCDo6+RB7dz+VPMpTinP/JT0yItnOjfKgbFkduShRRwPv0yXInXvLcA4twbHykjUNHDtke0rqHIKPdGsaRlg27ZtADRo0CCba5K1xHIEiZ4DMYuBWCM7q1uvhEN/nFN93NWYUHYFnWZn0Cl2BZ0iOD4q8Vo+Fy/8PfLh75kff498uJudibdbibNZiLVbiLNZibdbuRobyuLAALwc3Xi3wtMpZku9a/2jFyHh7xmH8tyRF+pRlO37FDKLDgqalvOIPczYDBc9B2znQXmDW1eUa8+7ppUQEa7GhnI64hqnI68lfjwXeR2L2FJ9nEk50K5gdV4v1y7daTiSP3cMcr0RODfGwfv7+7rHw0IHBU3LQo9rT+G/ROwQvx2JngtxawEbOD2BcuuZkEoj/SPWVruNwOggYm0WnE2OODs44mQy45Lw0VGZMmTIxx4+EqLnovJsQplSTx/+sNNBQdOy0KOyTyEjie0qxPyBRM83znZwyAsu7VCuHcBcKceM4YvlABLUFeX1Pcq1fXZXJ9Nk+45mTdMeb8qU3zjDIc96lPev4FgFomcjQV2Qm62xR/yEWM9kdzWBhHQeKmfkcspuevWRpmmZSikzuLREubQ05h5iVyGxSyDqVyTqF8RcEeX6FDi3RpmTn4mQ6SQ6oaL3ftDPo0gHBU3Tsoxy8AK3bii3bojtGsQuM47NjPgGIr5BzJVQLq2N1UvmYllTqcSgcPfd0I8DHRQ0TcsWypQP3F9Eub+IWC9A3CokdiUSORoiRyPmskaAcG4C5vJJUnpnFBEbErcloUK6pwA6KGhahhgzZkx2V+GhpsxFwfwyyv1lxHYZYlcisauQyJ8h8idQHohTTZRjbXCqA44Vk6T2vh9iOYGEvw+WA+DcAsylMqg1DzcdFDQtA1SrVi27q/DIUKaCt3sQthsQ/y8SvwvidyJxGxMKuSKO1Y3gYC4D5tJg9k9zw9wtIvEQ9ZtxhoTyQHl9Dy5P5ZjVUNlNBwVNywAP28lrDwtlygOu7ROXiootCCwBSPxOiN8FUdMQLAmlHRBTsYQAURIwARYQK2BN+BhvpAS3ngKXDqhcH6AcfLOlbTmVDgqalgFGjhwJ6KCQ2ZQpN5haG3MNJJwcZzsP1pOI5QRYTxr/4tYAdsARlBnjpc5sLDt18EF5/4ZyaZqNLcm5dFDQNO2hpZSjMRdgLoVyaZv4fRE7oPSQ0H3QQUHTtEfOrXMetHunf3KapmlaIh0UNE3TtER6+EjTMsBvv/2W3VXQtAyhg4KmZYCyZctmdxU0LUPo4SNNywBLlixhyZIl2V0NTXtguqegaRlg9OjRAHTo0CGba6JpD0b3FDRN07REOihomqZpiXRQ0DRN0xLpoKBpmqYl0hPNmpYBZs6cmd1V0LQMoYOCpmWAIkWy4WxhTcsEevhI0zLAvHnzmDdvXnZXQ9MemO4paFoGGDduHAA9evTI5ppo2oPRPQVN0zQtkQ4KmqZpWiIdFDRN07REOihomqZpifREs6ZlgAULFmR3FTQtQ+igoGkZwM/PL7uroGkZIlOHj5RSbZRSx5VSp5RS76Vw/U2l1BGl1AGl1FqlVLHMrI+mZZZp06Yxbdq07K6Gpj2wTAsKSikT8CvQFqgA9FJKVfhPsb1ALRGpAiwAvsms+mhaZtJBQXtUZGZPoQ5wSkTOiEg88DvQ8c4CIrJeRKITvtwBFM7E+miapml3kZlBoRBw8Y6vAxO+l5p+wPKULiilBiilApRSATdu3MjAKmqapml3yhFLUpVSzwG1gG9Tui4iE0SklojUypMnT9ZWTtM07TGSmauPLgF3po4snPC9JJRSLYAPgMYiEpeJ9dE0TdPuIjODwi6gtFKqBEYw6An0vrOAUqo68BvQRkSuZ2JdNC1TLVu2LLuroGkZItOCgohYlVKvASsBEzBFRA4rpT4DAkRkMcZwkQfwh1IK4IKIPJ1ZddK0zOLm5pbdVdC0DJGpm9dEZBmw7D/f+/iOz1tk5vNrWlYZO3YsAIMHD87mmmjag8kRE82a9rCbP38+8+fPz+5qaNoD00FB0zRNS6SDgqZpmpZIBwVN0zQtkQ4KmqZpWiIlItldh3uilLoBnM+Cp/IDbmbB8+Qkus2Ph8etzY9beyHlNhcTkbumhHjogkJWUUoFiEit7K5HVtJtfjw8bm1+3NoLD9ZmPXykaZqmJdJBQdM0TUukg0LqJmR3BbKBbvPj4XFr8+PWXniANus5BU3TNC2R7ilomqZpiXRQ0DRN0xI99kFBKdVGKXVcKXVKKfVeCtedlVLzEq7/q5Qqng3VzFDpaHNfpdQNpdS+hH8vZ0c9M4pSaopS6rpS6lAq15VS6qeEn8cBpVSNrK5jRktHm5sopcLu+B1/nFK5h4VSqohSar1S6ohS6rBS6vUUyjxSv+d0tvnef88i8tj+wzjn4TRQEnAC9gMV/lNmMDA+4fOewLzsrncWtLkv8Et21zUD2/wkUAM4lMr1dhjngyugHvBvdtc5C9rcBFia3fXMwPYWAGokfO4JnEjh7/qR+j2ns833/Ht+3HsKdYBTInJGROKB34GO/ynTEZie8PkCoLlKOBHoIZWeNj9SRGQTEJxGkY7ADDHsALyVUgWypnaZIx1tfqSIyBUR2ZPweQRwFCj0n2KP1O85nW2+Z497UCgEXLzj60CS/1ATy4iIFQgDcmdJ7TJHetoM0CWhi71AKVUkheuPkvT+TB419ZVS+5VSy5VSFbO7MhklYYi3OvDvfy49sr/nNNoM9/h7ftyDgpayJUBxEakCrOZ2T0l7dOzByIVTFfgZ+Ct7q5MxlFIewEJgmIiEZ3d9ssJd2nzPv+fHPShcAu58F1w44XspllFKmQEvIChLapc57tpmEQkSkbiELycBNbOobtklPX8HjxQRCReRyITPlwGOSim/bK7WA1FKOWK8OM4WkUUpFHnkfs93a/P9/J4f96CwCyitlCqhlHLCmEhe/J8yi4EXEj7vCqyThBmch9Rd2/yfcdanMcYqH2WLgecTVqfUA8JE5Ep2VyozKaXy35obU0rVwXgteGjf7CS0ZTJwVES+T6XYI/V7Tk+b7+f3bM7oij5MRMSqlHoNWImxKmeKiBxWSn0GBIjIYowf+kyl1CmMibue2VfjB5fONg9VSj0NWDHa3DfbKpwBlFJzMVZh+CmlAoFPAEcAERkPLMNYmXIKiAZezJ6aZpx0tLkrMEgpZQVigJ4P+ZudhkAf4KBSal/C994HisIj+3tOT5vv+fes01xomqZpiR734SNN0zTtDjooaJqmaYl0UNA0TdMS6aCgaZqmJdJBQdM0TUukg4KmJVBKeSulBid83kQptfQeH99XKVUwHeV8lVKrlVInEz763G+dNS2j6aCgabd5Y2TFvV99gbsGBeA9YK2IlAbWJnytaTmC3qegaQmUUrcyxh4HLEAUcBOoBOwGnhMRUUrVBL4HPBKu98XYSDQNI21CDFAfGA50AFyBbcDAhMcfB5qIyJWE3eMbRKRsVrVT09Kig4KmJUjINLlURCoppZoAfwMVgcvAVowX+X+BjUBHEbmhlOoBtBaRl5RSG4C3RSQg4X6+IhKc8PlMYL6ILFFKhYqId8L3FRBy62tNy26PdZoLTbuLnSISCJCQRqA4EIrRc1idkFLGBKSWP6epUuodwA3wBQ5jZKBNlNBz0O/MtBxDBwVNS13cHZ/bMP6/KOCwiNRP64FKKRdgLFBLRC4qpT4FXBIuX1NKFbhj+Oh6xldd0+6PnmjWtNsiMI41TMtxII9Sqj4YqYvvOLjkzsffCgA3E/Ldd73jHndm3n0BY5hK03IE3VPQtAQiEqSU2qqMw+5jgGsplIlXSnUFflJKeWH8HxqDMTQ0DRivlLo10TwROARcxUhZfstXwHylVD/gPNA90xqlafdITzRrmqZpifTwkaZpmpZIBwVN0zQtkQ4KmqZpWiIdFDRN07REOihomqZpiXRQ0DRN0xLpoKBpmqYl+j8Qjq32fCu7JwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m.draw_contour('theta0','theta1', bound=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_error_ellipse(cov, center, nstd, **kwargs):\n",
    "    \"\"\"\n",
    "    return error ellipse patch representing the covariance matrix\n",
    "    cov: 2x2 covariance matrix\n",
    "    centre: mean values (mu1, mu2)\n",
    "    nstd: number of standard deviations\n",
    "    \"\"\"\n",
    "\n",
    "    # eigenvalues and eigenvectors of the covariance matrix\n",
    "    # (in ascending order of the eigenvalue)\n",
    "    eigvals, eigvecs = np.linalg.eigh(cov)\n",
    "   \n",
    "    # x, y component of the eigenvector (eigvecs[:,1]) with the larger eigenvalue\n",
    "    x, y = eigvecs[:,1][0], eigvecs[:,1][1]\n",
    "    \n",
    "    # angle of the ellipse    \n",
    "    theta = np.arctan2(y, x)\n",
    "    \n",
    "    # width = 2 times radius in x, height = 2 times radius in y \n",
    "    height, width = 2 * nstd * np.sqrt(eigvals)\n",
    "        \n",
    "    return Ellipse(xy=center, width=width, height=height,\n",
    "                   angle=np.degrees(theta), **kwargs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArlElEQVR4nO3deXhdV33u8e9Ps2QNtgZLdmx5jBM7dkhiQ0lCQmjiQoG2NwOEMtW0kBtCoS1DS2gupb29cFumpNymNPC0aZJyGwhDSUsgEAgJIXViB3CcOJ7t2LKtwbIsyZql1T/WOdHxiSRrOGfvdeT38zz7kbWHs39n+0iv1t5rr23OOUREREKTF3cBIiIiY1FAiYhIkBRQIiISJAWUiIgESQElIiJBUkCJiEiQFFAiIhKkYALKzK40s++aWZOZOTPbdIb1rzKzfzezo2bWY2bbzOz3IypXRESyLJiAAsqB7cAfAb2TWP8y4FngBmAt8A/AXWb29qxVKCIikbEQR5Iws27gD51zd09xu68D+c6567NSmIiIRCakFlQmVAIn4i5CRERmriDuAjLFzN4MXA1cPs7ym4CbAObMmbP+/PPPj7A6EREZz9atW9ucc3Xp82dFQJnZ5cDXgA85554aax3n3F3AXQAbNmxwW7ZsibBCEREZj5kdHGt+zp/iM7PXAA8Bn3TO/UPc9YiISGbkdECZ2ZX4cPqUc+72mMsREZEMCuYUn5mVAysT3+YBjWZ2EdDunHvRzD4DvMo5d3Vi/auA/wTuBL5mZg2JbYedc61R1i4iIpkXUgtqA/CLxFQK/GXi33+VWL4AWJGy/iagDPgocDRlejqackVEJJuCaUE55x4FbILlm8b4ftNY64qISO4LqQUlIiLyEgWUiIgESQElIiJBUkCJiEiQFFAiIhIkBZSIiARJASUiIkFSQImISJAUUCIiEiQFlIiIBEkBJSIiQVJAiYhIkBRQIiISJAWUiIgESQElIiJBUkCJiEiQFFAiIhIkBZSIiARJASUiIkFSQImISJAUUCIiEiQFlIiIBEkBJSIiQVJAiYhIkBRQIiISpGACysyuNLPvmlmTmTkz2zSJbdaZ2U/NrDex3SfNzCIoV0REsiyYgALKge3AHwG9Z1rZzCqBHwLNwCsT230M+HAWaxQRkYgUxF1AknPue8D3AMzs7kls8g6gDPg951wvsN3Mzgc+bGZfcM65rBUrIiJZF1ILaqouBR5PhFPSD4CFwNJYKhIRkYzJ5YBqwJ/eS9Wcsuw0ZnaTmW0xsy2tra1ZL05ERGYmlwNqSpxzdznnNjjnNtTV1cVdjoiInEEuB9QxoD5tXn3KMhERyWG5HFBPAleYWUnKvI3AEeBALBWJiEjGBBNQZlZuZheZ2UX4uhoT3zcmln/GzB5J2eRrQA9wt5mtNbPrgI8D6sEnIjILBBNQwAbgF4mpFPjLxL//KrF8AbAiubJz7iS+xbQQ2AL8PfB54AvRlSwiItkS0n1QjwLjjgLhnNs0xrxngSuzV5WIiMQlpBaUiIjISxRQIiISJAWUiIgESQElIiJBUkCJiEiQFFAiIhIkBZSIiAQpmPugZGwjI9DcDE1N0NUFfX3Q2+u/DgyAGeTljU75+X4qKYHS0pdPZWVQUQHFxXG/MxGRiSmgAtXZCf/1X3D4MAwNZf71i4uhqgoqK/3X1KmoKPP7ExGZKgVUgPbuhcceg8HB7O2jvx9aWvyUrqTk5cE1b56fbNyxPkREMksBFZjBQfjpT7PTapqsvj4/pYdXQQFUV0Nt7ehUXe1PLYqIZJoCKjAnTsQbThMZGnp5qysvz7eskoFVVwc1Nf46mIjITCigAlNdDYWF2T29l0kjI3D8uJ927vTz8vL8+5g/3091dT7ERESmQgEVmIICWL589Jd9LhoZgbY2Pz3/vJ9XVORbWMnQmj/f9ygUERmPAipAr3gF7NuXO62oyRgYgCNH/JQ0Zw4sWAALF/qvVVXx1Sci4VFABWjuXHjd6+Dhh+OuJLtOnYI9e/wEo4GVDC0FlsjZTQEVqKVLYf162Lo17kqikx5YZWWjrasFC3xwi8jZQwEVsPXr/Wm+bdviriQePT0vD6zUU4IKLJHZTQEVuFe/GsrL4cknwbm4q4lXT4+/iXnvXv99MrCSoaXAEpldFFA5YO1a/8v4Jz+B4eG4qwlHemCVlo6GlQJLJPcpoHLE8uU+pB55xF+rkZfr7fW9H/ft89+XlvqgWrzYT6Wl8dYnIlOjgMohDQ1w/fV+KKSDB+OuJny9vae3sGprR8Oqvl7jCoqETgGVY0pK4PWvh+3bYfNmnfKbiuTNw7/4hb9xeNGi0cDSTcMi4VFA5ai1a/31lkcegY6OuKvJPQMDp58OrKk5vXWlAXBF4qeAymE1NXDddbBlCzz7rHr5zURyPMFf/tK3rs45ZzSw5syJuzqRs1NQfyea2S1mtt/M+sxsq5ldcYb1325mvzSzHjM7Zmb3mVlDVPWGoKDAd0W/9lp/jUVmbmAA9u/3z+T613+FBx7wp1OPHPHjDIpINMwF8me3md0I3AfcAvws8fU9wBrn3ItjrH858BjwUeA7QD1wJ3DCOXf1RPvasGGD27JlS0brD8HIiG9Jbd0a7iM7cl1h4emtq/LyuCsSyX1mttU5t+Fl8wMKqM3ANufc+1Lm7QYecM7dOsb6HwU+6JxbkjLvPcCXnHMT/tqYrQGV1NkJjz8OTU1xVzL7zZs3GlYLFujalch0jBdQQVyDMrMiYD3wubRFDwOXjbPZE8Cnzey3gP8AaoC3Ad/LVp25orIS3vQm371682bo7o67otnrxAk/bdvmW1eLF/txFBsb/bUsEZm+IAIKqAXygea0+c3ANWNt4Jx70szeBvwrUIp/Lz8Efm+s9c3sJuAmgMbGxsxUHbgVK/wvy23b/MX/2fT4jhANDo72DMzL8zcJL13qJ3VjF5m6nD0hYWZrgC8B/xvf+noD0AD841jrO+fucs5tcM5tqKuri67QmOXnw8UXw403wvnn6+bUqIyMwOHD8LOfwX33wXe+A7/6FZw8GXdlIrkjlBZUGzCM7+iQqh44Ns42twJPOec+m/h+m5mdAh43s0845w5np9TcVFYGV14JF1zgB55NfXCgZF9Li582b/bXrZItq7PobyWRKQsioJxzA2a2FdgIfCNl0Ubgm+NsVoYPtVTJ73O2ZZhtNTXw5jfDgQP+l6X+oo9e8rrVL37hewEmw6qhQZ0sRFIFEVAJXwDuNbOn8B0gbgYWAl8GMLN7AJxz706s/yDwFTN7P/ADYAFwO/DMWN3S5XTJC/nPPee7pQ8MxF3R2am72w9btX07FBfDkiX+/2bRIn+Pm8jZLJgfAefc/WZWA9yGD5vtwBudc8lhURvT1r/bzCqAPwQ+D5wEfgz8WXRV57a8PFi3Ds4914fUjh25eSNqYe9JyjqbKeztxOXlgeXhklNePsMFRQyWVDJQUslIYXHc5Y6rvx927fJTQYEPqaVLfWgVh1u2SNYEcx9UlGb7fVDT1dXlTzvt2pUbQVXa2czCF35MYV/XpLcZKShioKSSwVIfWAOllQyUVtFfVh1seOXl+dN/y5b5wNLQSzLbBH+jbpQUUBPr7IRnnoHdu8Me32/JL79DaWdLxl5vsKScvvLalKmO4aLwHiJVVzd63WrevLirEZm5oG/UlbBUVsJVV8Ell/hTf3v2hBlUfeW1GQ2owr5uCvu6qWg78NK8wZJy+irm05uY+sprcfnx/ti0tvrp6af9U4OXL/f3vCmsZLZRC0rO6ORJH1R794YVVPkDPSz95Xco7ItuqAyXl0dfeS09VQvpqVpAT1UDLr8wsv1PZN48H1QrVkBVVdzViEyeTvGlUEBNT0eHD6p9+8IJquJT7Sz51XfJG4qnG+JoYC1IhFYYgVVTMxpWFRVxVyMyMQVUCgXUzHR0+KGT9uwJozNFSWcL57zwo0hbUuMJMbDq6nxQLV+u0dclTAqoFAqozOjq8kG1c2f8QZU32M+C3Y9R0bY/3kLSmdFbMZ/u6sV0VzfSXx7vQ7vq60fDSuMDSigUUCkUUJl16pQfZ+6FF+J/DtXcI89Tv+9JbCR9kJEwDBeV0j1vMafmLaJ73uJYu7YvWDAaViUlsZUhooBKpYDKjt5eP3L688/HO3J68al2Fr7wCMWnTsRXxGSY0VtRR3d1I6fmLaavIp6B+fLy/E3BydHvC+O/hCZnGQVUCgVUdvX3+yf7Pvec/3csRkaoPrKd2hefia0DxVQNF5b4llV1I6fmLWK4MPpmTUGBHwJr5Ur/bKv8/MhLkLOQAiqFAioaQ0N+VIrt233HijjkD/RSd+Bp5jbvDKfr4WSY0VtZT1fNUrpqlzFYEn1XvKIiP3rFypX+2VZ6VItkiwIqhQIqei++6E//xfWYj+Lu49Tv/TllJ4/GU8AM9ZfX+LCqWUp/eU3k+y8t9acAV66E+fMj373McgqoFAqo+Bw/7k//xdVFvaJ1H/P3b57S+H2hGSypSLSsltJb2RB506aiwg8wvGqVH3VEZKYUUCkUUPHr6fGdKZ5/Hvr6ot23jQwzr2k7NYd/Rf5gxDvPsOHCErprltBZu4yeuefg8qK9aNTQ4INqxQp1rpDpU0ClUECFY2jID0r77LPRX6ey4UGqjzxH9eFtOR9UACP5hZyqXkxn7TK6q5dEOmZgQYHvAbhqFZxzjq5XydQooFIooMJ06JC/TtXUFO1+bXiQ6qbtVDdtI38wrm6HmeXyC+iubqSzbgXd1Y2RtqzKy/21qvPO05iAMjkKqBQKqLC1t49epxqO8H7bl1pUTc+SP9Ab3Y6zbCS/0J8GrFtB97zFkT5Xfv58H1QrVvhegSJjUUClUEDlht5efy9V1NepbHiIucdeoObwryjoPxXdjiMwUlBEV+0yTtatpGdudH3H8/NHTwEuWqRTgHI6BVQKBVRuGR4evU51IsrBIUZGqGrZRXXTs+GPSjENw0WldNatoLNuBb2V9ZHtt6xstBegnmEloIA6jQIqdx0+7K9THT4c7X7nnDhM9eFtzDkR8Y4jMlhSQef8lZycfy4DZXMj229dnQ+qlSuhOL5hCSVmCqgUCqjcd+IE7NjhR6oYiHAko6JTJ6huepaqlt3BDkg7U72V8zk5fxWd81cyUhDNhaO8PFiyxIfV4mgvk0kAFFApFFCzx9CQf9Lvjh3Qkrmnv59R/mAf8448x7yjz8+qDhWpXF4+3TVL6Khfxal5iyO7cFRaOtoLsLo6kl1KzGIJKDMrAl5wzi3P2k6mQQE1O7W1+aDasye60dRtZJjKlj2J61Tt0ew0BkNFZXTOX0lH/XkMzInuwlFNjQ+qlSv1SJDZLK6AKgZ6nXNBNdgVULPb4KDvVLFjhx9aKSplJ5qoPrKd8vYXc2tg2inqK6/lZMN5dNatiGzE9fx8P3Dt6tX+OVYyu2QtoMxs3wSL84DFzrmgBu1XQJ09Wlp8N/V9+6J7mGJBXzdzj73A3GMvUDDQE81OY+Dy8uiuXsLJ+lWR3l81b54PqlWrdG/VbJHNgOoCPgG8OMbiIuDfFFASt4EB36Fix44Iu6qPjFDe/iLzjj4/a3v/JQ0VldHRcD4dDeczVFIeyT4LCvzTgNes0QjruS6bAfU48CXn3NfHWDalU3xmdgvwMWAB8Bzwx865xydYvwi4DXgXsBBoBj7nnPu7ifajgDq7HTvmW1X790c3UkVhb6dvVTXvnLWdKgAwo7u6kRMLVkfasaKmxreqzj1Xg9bmomwG1FuAdufcI2MsywPe5Zz7l0m8zo3AfcAtwM8SX98DrHHOjdU6w8y+BSwC/hzYDdQDpc65RyfalwJKwI9OsWuXD6vOzoh2OjJCxfH9zDu6g7KOmB6OFZHBknI6GlbT0XAew0VlkeyzsNB3qFizxoeW5Ibgu5mb2WZgm3PufSnzdgMPOOduHWP93wC+AaxwzrVNZV8KKEnX1OSD6uDB6J5TVdTTwdyjO6hq2TVrBqkdS/Ja1YkFa+iZd05k+50/37eqVqzwpwMlXDMOKDMrB9YCVUAXsM85dyxDxRUBPcDvOue+kTL/74G1zrnXjrHNncAq4Cng3UAv8BDwCedc9xjr3wTcBNDY2Lj+4MGDmShdZpmeHti501+r6n7ZpyhLRkaoaD9I1bGdzOk4jMXxJMeIDJRW0bFgNSfrV0XWA7CoyHeoWL1aQyuFatoBlbiOdDuwCd/pIdUx4BHgfuB7bprNMTNbCDQBr3XOPZYy/5PAO5xz542xzfeBqxL7/ytgLvAlfCvshon2pxaUnIlz/vEfO3b4x9VHdaIhf6CXqpbdVDXvmtX3Vbm8fDrrVtB+zrpIH2Hf0OBP/y1b5ruuSxjGC6jJNHw/B/xP4AfAo8AA0IDvzFACvBN4B/CCmX1orGtRWZIHOODtzrmTAGb2h8APzKzeOdccUR0yC5lBY6Ofenr8zb87d2a/B+BwUSntiy6kfdGFFHe3Mbd5F5Ute2bFAxVT2cgwVc27qGreRc/chRxfdCGnqhuzvt9jx/xUUuJbVWvW6LH1IZtMQN0I/JNz7r3JGWZWgw+otwD7El/fjw+HjznnvjjFOtqAYXwnh1T1+FbaWI4CTclwStiR+NqI79EnMmNlZXDhhX5qbfUdK/buzf4jQPrLa2kur6V52aspb3+RqpZdlLe/OOtOAZZ1HKGs4wgDZXNpX7iWk/Wrsv404L4+P+jwtm3+CcCrV/vHgWgMwLBM5lNQCjw53kLn3AHgs2Z2O/Bp4HNmtt0598PJFuGcGzCzrcBGfMeHpI3AN8fZ7AngLWZWnnLNaVXiqy4wSVbU1fnp0kt9h4pdu/ypwKxmRl4e3bVL6a5dSv5gH5Ute6hq3kVJ95T6BgWvqKeDhj0/o+7g03Q0rKb9nLWR9P5ravJTWZkfVumCC/y/JX6TuQb1E3xL5Z0p82qAVuAa59yP09b/FlDtnLtqSoX4bub34ruXPwHcDPwBcIFz7qCZ3QPgnHt3Yv1yfIvpv4BP4a9B/SOwwzn3lon2pWtQkkm9vf4U4K5d0Q6tVHyqncqWPVS27qWwryu6HUfE5eWlXKeqjWy/eXm+59+6dVAb3W7PajO5BvV/gYcSLzCZU3ffA6Z6ig/n3P2J4LsNf6PuduCNzrlka6gxbf1uM7sG3zHiaeAE8B3g41Pdt8hMlJb6X2br1vmA2rXLjwWY9VOAc6ppXfYqWpe9ipLOFqpafVjNlhuBbWSEqubdVDXvpqdqAccXX8Sp6sVZ3+/IiP//273bj/t34YX+WqSeAhy9SXUzN7OPAH8LbAX+Dn8j7T7GbkHdC/yGcy66R3ROkVpQkm0jI/7U386dvhdgZJeNnKOs4whVrXuoaNtP3lCED8uKQF9FHW2LL6a7dmmk+62qgrVr/SlA3VOVeZm4D+oa4A5gNTAE5APfxrdeuoBK4PXAlcAdzrkPZ6b0zFNASZT6+kZPAbZFedloZIQ5HYepbN1LxfGDsyqs+udU09Z4MV21yyNt2hQX+w4VF1wAc+ZEtttZLyMjSZiZAW/A9+y7irTTbvibbb8KfMw5F9ETeaZOASVxaW/3QbVnj+++HhUbGWZO+yEfVu0HseGIhnbPsoGyubQuWR95UOXl+YFqL7xQ16kyIStDHSWuGS0HyoFO4HnnXPAnwBVQErfkjcC7dvnegFENWgtgw0OUnzhE+fEDs6Zl1T9nHm2N6+mqi/7ZqA0NPqiWLNF1qukKfiy+KCmgJCT9/f6+ql27on1sPQAjI5SdPErF8f1UtB3I+edX9c+ppmXZr0XSmSJdZeXodSqNqD41CqgUCigJVUfH6I3AXTH0HC/pbKHi+AEqjh+gqKcj+gIypLu6kZblr2agbG7k+y4qgvPP92FVHs2jsXKeAiqFAkpyQXOzD6q9e/29VlEr6umg4vgByo8foLQz6qbdzLm8PE4suIC2JesZKYj+0bt5eX7Mv3Xr9EDFM1FApVBASS5xzo90sHevf8jiQAyXjAr6T70UVmUnj+bUcEvDhSW0LtlAx4LVsV0kqq/316mWLtV1qrEooFIooCRXDQ/7zhV79vj7q4Zi6Ixnw4PM6TjCnBOHKG8/lDOjWPTPqebYysvprVoQWw2VlXDRRX6gWo37N0oBlUIBJbPB4CAcOODDqqkpwpuB0xT1dLwUVr51FWGXxGk4vvgVtC55ZawJUV4Or3iFv1alx34ooE6jgJLZpq8P9u3zpwGPHYvu+VXpbHiIspNHKG8/RPmJQxT2dsZTyBn0Vs7nyPlXM1hSEWsdyZHy16w5u0eoUEClUEDJbNbTMxpWzTE/dKaw9yTl7YeYc+IwZZ3HgrrnaqSgiKPnXkFX3Yq4S6GkxHemuOAC3wvwbKOASqGAkrNFd7cPq337YrjHKp1zlHS3UXbyqH8GVCCBdWLhGppXvibuMgAfTmvX+rAqLo67mugooFIooORs1NU12rKKdEzA8ThH8anjzEk8sDDOwGprvIS2pS/7/RibwkIfVBdeeHYElQIqhQJKznadnaP3WLW3x11NQmpgnTxKaVdLpI8OOXL+r9M5f2Vk+5uMZIvqwgtn96k/BVQKBZTIqI6O0XusggmrhIK+bkq7Wl6aSrrbsjbQ7UDZXPZteGtWXnumiopGnzk2G4NqJg8sFJFZbO5cWL/eT52dvuv6gQO+g0Xcf78OlZTTVVI+OgiscxT1nqSku+2lqfjUcfIH+2e+s7jf7AQGBmDrVti+3bem1q49O8b7UwtKRMbU2+tHWj9wwN9nFeWI61OVN9hPUV8nhX2dFPV2+n8nvhb0n5pw28GScnqqFtK69JUMFefGQ55KSkaDajZ0T9cpvhQKKJGpGRz0I1js3++/xjHc0nTZ8BB5wwPkjQyDG8FGRjDnp6GispwJpbGUlcGGDX4E9VweQkmn+ERk2goL/QP6li/3I1Y0NY2eCoxjINupcPkFDOcXEHADcNp6euCxx+C55+DSS2Hhwrgryiy1oERk2pzz91ft3+/DqjPMgSPOGkuXwqtf7cf8yyVqQYlIxpn5kbrr6/0vxvb20ZZVEPdanWUOHPCDCK9dC5dckvs9/hRQIpIx1dV+uuQSP4pFMqyOHg26k9ysMjIC27b5B1+uXw+rV+fuyOkKKBHJivJy/5f82rV+MNtkj8DDh8PuEThb9PXBE0+MXp9avDjuiqZOASUiWVdS4nuanXeef4bVoUM+rA4d8r9IJXs6OuChh2DRIrj8cqiqiruiyQuq4Wdmt5jZfjPrM7OtZnbFJLd7jZkNmdn2bNcoIjNTUOAfhf6618G73gXXXuu7Ss+fn9tdpUN3+DA88AA880x8zw6bqmBaUGZ2I3AHcAvws8TXh8xsjXPuxQm2mwfcAzwCnBNFrSKSGWZQV+enSy7xranDh33L6vDh8Luw55rhYdiyxQ8afOWV/o+CkAXTzdzMNgPbnHPvS5m3G3jAOXfrBNt9C/gVYMANzrm1Z9qXupmL5Ia2Nt8r7dAh3509kF9Xs4KZf/7UK18Z/7BJQXczN7MiYD3wubRFDwOXTbDdLUA98NfA/8pagSISi9paP11yCfT3+xuEX3zRt656euKuLrc558f2O3AArrgizE4UQQQUUAvkA+nP/2wGrhlrAzNbB/wF8Grn3LCd4eS1md0E3ATQ2Ng403pFJGLFxaOjWYBvXR06NNq6ypXrKqHp7vadKNau9feyhdQlPZSAmhIzKwbuBz7qnNs/mW2cc3cBd4E/xZfF8kQkAsnW1cUX+7EBk9euDh1S62o6tm+HY8fg6qvD6ekXSkC1AcP403Wp6oFjY6y/AFgN/LOZ/XNiXh5gZjYEvNE593C2ihWRsBQVnd66On58NKyam9W6mqy2NvjWt3x39FWr4q4mkIByzg2Y2VZgI/CNlEUbgW+OsUkTsC5t3i2J9a8FDmShTBHJETU1frroIt+6amqCI0f8dOJE3NWFbXAQHn3Uj/5xxRXxnvILIqASvgDca2ZPAU8ANwMLgS8DmNk9AM65dzvnBoHT7nkysxag3zmne6FE5CVFRf6+q2XL/Pd9fT6ojh5VYE1k5044dQo2boyvl18wAeWcu9/MaoDb8KfwtuNP1R1MrKKeDSIyYyUlp58O7OsbDaujR8N77H2cDh+GBx+EN7zBP3sqasHcBxUl3QclIuNJBlYytBRYUFEBb3pT9h7jEfR9UCIioSgpefkpwWPHRq9hnY2B1dUFDz8Mv/M70Z7uU0CJiEygpMQ/CHDpUv99f//ppwSPH4+zuui0t/vOExs3RrdPBZSIyBQUF48dWEeP+i7tx4/P3seJ7N8Pzz4L69L7UGeJAkpEZAbSA2tkxIdUa6sf4aKlxT/yYrZ47jkFlIhITsrLGx2hfc0aP29gwAdWamjl6mgXnZ3+mlxDQ/b3pYASEcmyoiI45xw/JfX0+JEbUqfu7vhqnIqowlUBJSISg7IyaGz0U1J//+mBdfw4nDwZ1mNGCgpOD9qs7iua3YiIyJkUF7+8pTU05E+rdXb6sDp5cvTfp05FW19VlR9Mtrg4mv0poEREAlZQANXVfkqXDK/04Orq8qfhMjFIbkEBLFgAixbB6tX++6gooEREctRE4QX+lGFv78unoSEfXulTYSGUlvp7v0pLYc4c39kjPz/a95WkgBIRmaWKi/00d27clUxPQM9OFBERGaWAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAEhGRIAUVUGZ2i5ntN7M+M9tqZldMsO51ZvawmbWaWZeZbTaz346yXhERyZ5gAsrMbgTuAD4NXAz8HHjIzBrH2eS1wI+BNyXW/x7w7YlCTUREcoc55+KuAQAz2wxsc869L2XebuAB59ytk3yNp4DHnXMfmWi9DRs2uC1btsyoXhERyQwz2+qc25A+P4gWlJkVAeuBh9MWPQxcNoWXqgBOZKouERGJTxABBdQC+UBz2vxmoGEyL2BmHwAWAfeOs/wmM9tiZltaW1tnUquIiEQglICaETO7Hvgs8Hbn3MGx1nHO3eWc2+Cc21BXVxdtgSIiMmWhBFQbMAzUp82vB45NtKGZ3YBvNb3bOfdgdsoTEZGoBRFQzrkBYCuwMW3RRnxvvjGZ2Vvx4bTJOfdA9ioUEZGoFcRdQIovAPcmeuI9AdwMLAS+DGBm9wA4596d+P5t+HD6KPCYmSWvVQ0459ojrl1ERDIsmIByzt1vZjXAbcACYDvwxpRrSun3Q92Mr//2xJT0U+CqbNYqIiLZF0xAATjn7gTuHGfZVRN9LyIis0sQ16BERETSKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgKaBERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgqSAEhGRICmgREQkSAooEREJkgJKRESCpIASEZEgBRVQZnaLme03sz4z22pmV5xh/dcm1uszs31mdnNUtYqISHYFE1BmdiNwB/Bp4GLg58BDZtY4zvrLgO8l1rsY+AzwJTO7PpqKRUQkm4IJKODDwN3Oua8453Y45z4IHAXeP876NwNHnHMfTKz/FeBfgI9GVK+IiGRREAFlZkXAeuDhtEUPA5eNs9mlY6z/A2CDmRVmtkIREYlaQdwFJNQC+UBz2vxm4JpxtmkAfjTG+gWJ1zuausDMbgJuSnzbb2bbZ1JwBGqBtriLmEDo9UH4NYZeH4RfY+j1Qfg1hlDfkrFmhhJQWeecuwu4C8DMtjjnNsRc0oRCrzH0+iD8GkOvD8KvMfT6IPwaQ64viFN8+PQeBurT5tcDx8bZ5tg46w8R/18DIiIyQ0EElHNuANgKbExbtBHfS28sT46z/hbn3GBmKxQRkagFEVAJXwA2mdl7zWy1md0BLAS+DGBm95jZPSnrfxk4x8xuT6z/XmAT8LlJ7OuuDNeeDaHXGHp9EH6NodcH4dcYen0Qfo3B1mfOubhreImZ3QL8KbAA2A78iXPuscSyRwGcc1elrP9a4IvABcAR4G+cc1+OtmoREcmGoAJKREQkKaRTfCIiIi9RQImISJBmRUBlY5DZqb5mpuozs+vM7GEzazWzLjPbbGa/nbbOJjNzY0wlEdV41Tj7Pz9tvevN7Hkz6098vTai+u4ep75TU30Pk6ztSjP7rpk1JV5j0yS2WWdmPzWz3sR2nzQzS1snk8dvSjUmjs+/m9lRM+sxs21m9vtjrBPLMTSzpePs+w1p62VsQOlp1PipcWp0ZjZ/Ku9jkvXdamZPm1ln4vfHg2a2dhLbRfpZnBLnXE5PwI3AIPA+YDXwJaAbaBxn/WXAqcR6qxPbDQLXT/c1M1zfHcDHgVcBK4G/wN8jdkXKOpsS76EhdYrwGF4FOGBNWg35Ketcir8n7c8Tr/nnie9/LYL6qtKPDbAX+OepvIcp1PdG/CDHNwA9wKYzrF+Jv4/v68DaxHZdwEeycfymWeMngL8GLgeW48fEHALeHsgxXJrY9+vT9l00lZ/1LNdYPsbn8FHgJ1N5H1Oo7wfAexKfqXXAtxOfs+qQPotTek/Z3kHW3wBsBr6SNm838Jlx1v8bYHfavK8CT073NTNZ3ziv8RTw+ZTvNwHdMR7D5C+m2gle837gh2nzfgT8/6iPIf6XrAMum8p7mOax7J7EL673A51Aacq824AmRjsuZez4TafGcbb7OvDNQI5h8hf7hgnWOePPepTHEFiM/2MzNeTP+D5mUGN5Yn+/Fepn8UxTTp/isywMMjvN18xkfWOpAE6kzSs1s4NmdtjM/sPMLp5KbRmqcUviFNAjZva6tGXjHec4juH7gOecc2Pd9D3Re8iWS4HHnXO9KfN+gL/vb2nKOjM+fhlWycs/hxDPMUz6lpm1mNkTZnZD2rLQBpT+A/zx++YYyyZ6H9NVgb+MM9b/WVLQn8WcDigmHmS2YZxtGsZZPznI7HReM5P1ncbMPgAsAu5Nmb0T+H3gd4DfBfqAJ8zs3CnWN90ak49BuR64LlHPI3b6daHxjnOkx9DMqoC3Al9JWzSZ95At4x2b5LKJ1pnq8csIM3szcDWn39QZ5zHsxj9a5634U2+PAPeb2TtT1jnTz3pkzCwf/zN7r3OuP2XRZN7HdN0B/BI/6s54gv4snjWDxeYi8w9f/Cxwo3PuYHK+c+5JUj50ZvZz/Afxg8CHsl2Xc24n/pdR0pNmthT4GPB4tvc/Re/E/yGWGvC59h5iZWaXA18DPuSceyo5P85j6JxrAz6fMmuLmdXib/S/L5v7nqY34E/xnfaHUrbeh5l9AXgN8Brn3PB0Xyduud6CysYgs9N5zUzWB0CimX8v8G7n3IMTrZv4AG4BptOCytT73Zy2//GOc2THMOF9+Osm7ZNYN/09ZMt4xya5bKJ1pnr8ZsTMXgM8BHzSOfcPk9gkqmM4mX2HNKD0TcDPnXPPT2LdGR1DM/si/szKrzvn9p1h9aA/izkdUC4Lg8xO8zUzWR9m9lZ8OG1yzj1wpv0kuoReSNozsLJZ4xguStv/eMc5kmMIYGavAl7By0/vjecipnEMp+FJ4Ao7/baAjfjhug6krDPj4zcTZnYlPpw+5Zy7fZKbXUQ0x3Ay+w5iQGkzWwi8iQg+h+bHME2G0wuT2CTsz2K2e2Fke8J3QR4A3ovvAnkH/rzuksTye4B7UtZPdj29PbH+exPbp3czH/c1s1zf2/BdYf+I07udVqes8xf4bqnL8R/mf0ps86qIjuEfA/8D/1feBcBn8D2RrktZ5zL8X6ofB84Hbk3UON1u5pOuL2W7rwK7xnnNM76HKdRXnvh/uAjf/fiTiX83JpZ/BngkZf0q/F+f/4bv2nsdvidVatfejB2/adZ4Ff7n5LNpn8O6QI7h7wFvT3wezsNfxxnAj9856Z/1bNaYst1twEmgbIxlZ3wfU6jv7xOfo19P+z8rT1kn9s/ilN5TtncQxQTcgk/7fvxf21emLHsUeDRt/dcCzyTW3w/cPJXXzGZ9ie/dGFPqOl8EDiZerwXfo+bSqI4h/vz4bqAXaMdfb3jjGK95A/BC4gduB9P4xTWD/+MKfIj96TivN6n3MMnarhrn/+zuxPK7gQNp26wDHsN3cDmK/6PDsnj8plRj4vux1k9dJ7ZjiP/F/jw+gDrxp7jfOcbrnvFnPcv/z5bY753jvOak3sck6xurNodvATNBjZF+FqcyabBYEREJUk5fgxIRkdlLASUiIkFSQImISJAUUCIiEiQFlIiIBEkBJSIiQVJAiYhIkBRQIiISJAWUSA4ws9eb2aNm1p14nPf/Sxs/TWTWUUCJBM7MPgJ8Hz8MzZ8ADwIfwI9JKDJraagjkYCZ2TX4p5n+qXPucynzvw+8Dj94a2dc9Ylkk1pQIoEyszx8K+kXnP5QO/AD5BbhR6AWmZX0RF2RcL0eWIN/Llj6qY6BxNeqaEsSiY4CSiRcN+KfJvx44jHgqZJPOO2KtiSR6OgalEigzOwg0HiG1c5xzh2Joh6RqCmgRAKUaDG1At8G7hxjla8D/c65BZEWJhIhneITCdPyxNennXM/Sl1gZsuAecDXUuYV4DtSvAvf+embwAecc33RlCuSeerFJxKm8sTXsa4x3ZD4en/KvE/gu52vA87Fd67426xVJxIBBZRImJL3NlWmzjSzIuD9wE7gP1MWvRf4tHOuyTnXCnwK2GRm+RHUKpIVCiiRMD0P9OC7mqf6P8BS4EPOuWEAM5sLLAZ+mbLeM0BFYl2RnKRrUCIBcs71mNlXgQ+Z2X3AT4HfBK4FPuacezhl9YrE146UeR1py0RyjgJKJFwfAxzwDnwwbQV+0zn3/bT1ktepqoBjiX/PTVsmknPUzVxkFjCzF/Etq/sT3/8G8AAwL3kqUCTX6BqUyOzwVeBWM1toZnX4ThJ3K5wkl+kUn8js8GmgFngO/4fnA8CfxVqRyAzpFJ+IiARJp/hERCRICigREQmSAkpERIKkgBIRkSApoEREJEgKKBERCZICSkREgvTflbFzEA9tYaUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.patches import Ellipse\n",
    "fig, ax = plt.subplots()\n",
    "cov = np.array(m.covariance)\n",
    "el_1sigma = get_error_ellipse(cov, (theta0, theta1), 1, fc='red', alpha=0.4)\n",
    "el_2sigma = get_error_ellipse(cov, (theta0, theta1), 2, fc='blue', alpha=0.4)\n",
    "ax.add_artist(el_2sigma)\n",
    "ax.add_artist(el_1sigma)\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "ax.set_xlim(0., 2.2)\n",
    "ax.set_ylim(0., 1.2)\n",
    "plt.xlabel(r\"$\\theta_0$\", fontsize=18)\n",
    "plt.ylabel(r\"$\\theta_1$\", fontsize=18)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"basic_chi2_fit_plot2.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error band around fit function (propagation of fit parameter uncertainties)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "# error propagation of the uncertainties of the fit parameters (see lecture slides)\n",
    "Var_yf = cov[1][1] * xf**2 + 2. * cov[1][0] * xf + cov[0][0]\n",
    "sigma_yf = np.sqrt(Var_yf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqUElEQVR4nO3de5TcZZ3n8fc396TT6c6NXAidCIGEO0IjFwUCigrZcdTV8RIRFGkVUAZndkfMeLzshPHszOpxHTy7YWYH1s0ZnTM74yynE7koiAMESZThjoYk3eTSCbl30l1dl372j2//qErTSaq6q+r3q67P65ycpCvVVU+H0J88z/f7PI+FEBAREUmaMXEPQEREZCgKKBERSSQFlIiIJJICSkREEkkBJSIiiaSAEhGRRFJAiYhIIsUaUGY2z8zuN7M3zCxlZi+Z2VVxjklERJJhXFxvbGbNwBPAvwHLgTeAU4HdcY1JRESSw+I6ScLM7gauCiG8M5YBiIhIosW5xPdB4Gkz+4mZ7TazZ83sdjOzGMckIiIJEecMKjXwy+8B/whcAPwA+GoI4W+GeH4b0AbQ0NBw0dKlS6s0UhERKZeNGzfuCSHMLua5cQZUGtgQQri84LG7gQ+FEM483ue2traGDRs2VHqIIiJSZma2MYTQWsxz41zi2wm8NOixl4GWGMYiIiIJE2dAPQEsGfTYGUBHDGMREZGEiTOgvgdcamYrzWyxmX0U+DJwT4xjEhGRhIgtoEIIz+CdfH8EvACsAr4O/DCuMYmISHLEtlEXIITQDrTHOQYREUkmncUnIiKJpIASEZFEUkCJiEgiKaBERCSRFFAiIpJICigREUkkBZSIiCSSAkpERBJJASUiIomkgBIRkURSQImISCIpoEREJJEUUCIikkgKKBERSSQFlIiIJJICSkREqiMExpaQO7FeWCgiInUgl4M9e2DTJhpharGfpoASEZHKSKdh1y7YtAlSKWhsLOnTFVAiIlJePT3w+uuwdSuEAE1NMG1ayS+jgBIRkfI4eNBDaccOGDcOpk+HsWOH/XIKKBERGb7+fti3z5fx9u+HiRNh9mwwG/FLK6BERKR0mYzXl37/e68vNTTASSeV9S0UUCIiUrzeXti+HTZv9tnTtGnDqi8VQwElIiIndugQdHR4OJlBc7PXmSpIASUiIkMLwetLr70Ge/fChAkwcyaMqc4ZDwooERE5WjYLb7zh9aUjR2DKlLLXl4oRW0CZ2TeBbwx6eFcIYW4MwxERkb6+fH0pk/HaUgzBFIl7BvUqsKzg41xM4xARqV+HD0Nnp/8YM8Y31la4vlSMuEeQDSF0xTwGEZH6E4LvW9q82Zfzxo+van2pGHEH1KlmtgPoA54GvhZC2DzUE82sDWgDaGlpqd4IRURGk1zOA2nTJujuhsmTY13GO544A+pp4CbgFeAk4M+BJ83s7BDC3sFPDiGsBlYDtLa2hiqOU0Sk9vX1wc6d3pGXTsdeXypGbAEVQlhX+LGZrQc2AzcC341lUCIio83hw35wa2enf9zU5HuYakDcS3xvCiEcNrMXgdPjHouISE0LAQ4c8PrS7t3e8DBjRqLqS8VITECZ2SRgKfBo3GMREam0Zcv858ceK+OLRvWl117zkx8mTy7bwa1xiHMf1F8DDwCdeA3q60ADcH9cYxIRqUnpdL6+1NfnFwMmqb50+DD84hewbt2Jn1sgzhnUAuAfgFnAG8B64NIQQkeMYxIRqR1HjuTrS9HFgE1NcY/KZbPw9NOwdq1PE/v64JRTSnqJOJskPh7Xe4uI1KyovrR1K3R1eX2puXlEFwOWTQjw6qvQ3g4PPeTn9zU1wQc+ANdfD+ecAxdfXPTLJaYGJSIix5HLwZ49vn/p4MFk1Ze6uuBnP/PZ0ubNvun3iis8lN75Tv94GBRQIiJJFtWXNm/2iwEbG2HOnLhHla8rrV0LGzf67On88+Guu+Daa8tyR5QCSkQkiY4cgW3bfCkPfKmsQhcDFi2bhfXrPZR++ct8XamtDa67DhYsKOvbKaBERJIiBF++27IlX1+aPj3e+lII8MorHkoPPuj3Qw2uK1VomVEBJSISt8L60qFDMGlS/PWlri5vC1+71gOzTHWlUiigRETiEvqh4/Wj60tx7l86fBh+/nMPpd/8xmdPF1wAX/savOc9VV9iVECJiFTbkSPQY9CX8uWzOOtLUV2pvR0ef9zrSi0t8PnPe13p5JPjGRcKKBGR6hi8f6nvEq8xzZoVz1heecVD6cEH/V6oqK60fDmcfXYi2tcVUCIilTT4fLyovjSu8jWctxiqrnTllV5XuvzyqtSVSqGAEhGphOj+pc2b4z0fr7CutHGjP3bBBbBypdeVGhurP6YiKaBERMppqPuXqn0+3rHqSl/4Qux1pVIooERERioEr+Ns2RLf/UshwMsv5/crJbSuVAoFlIjIcOVyHkivvQbd3fGcj7dzZ76utHUrTJiQ36+UwLpSKRRQIiKlSqVgxw6vL2Uy3iJezfrS4cPwyCP5/UoAb387rFiR+LpSKRRQIiLF6u722tLrr/vyXVOTL+dVQzYLTz3loVTDdaVSKKBERI6nv9/rOa+95vcbTZgAM2eOqL60Zt101j/fQF/GWPQfzmHVbdtZcd3+tz4xBHjpJQ+lhx7ycTQ3wx/+oS/h1WBdqRQKKBGRoWQysGuXB1NPD0yZUpZlvDXrptO2aiF9GQ+4jq6JtK1aCJAPqaHqSoX7lao1a4tZfXyVIiLF6u31ay62bPEmiKamstaXVt5zMj2po08n70mNZeXfzGdF398fXVe68EL41Kfg3e+u7bpSX58f75TNEiAU+2kKKBER8Gsutm715oexYytWX+rcNeHYj//FX+TrStdfD/Pnl/39qyKX81lnKuXLlA0NsGgRzJrFIegu9mUUUCJSv3I5rytt2uTn5FXhmouWOWk6uia+9fEpe+GH99VuXSmV8llSLucBP2cOzJ3rHY6TJ7/5tKKnTyigRKQeDb5GferUyl+jvnMnrF3Lqsw42vhLemh487emTMqx6q5uv/yvVmSzHkh9ff5xczOccYZfsNjYWJZNygooEakfQx1DVMlrLobYr7TiwgvhirXc0v4hejNjWTg3fewuviQJwetzPT3e2Thxos+QZs/2P8cJQy9djoQCSkRGtxD8mvItW/xU8UofQ5TNwpNP5vcrpdOwcCF88Yu+X2n+fFYA93b2AvDY6t9VZhzlkMl4yGYy/uc1Y4bXkqZP97pShZciFVAiMjpls34M0aZN/k22oaFypz2EAC++mN+vdOCAL3l96EPe7HDWWbVRV+rv9xlSr4cnkybBKaf4LGnatKq3tyugRGR0GdwmPm1a5epLO3bk9yt1dPgy11VXeShddllt7Fca3Nwwa5bXkpqafO9XjBLzp2dmdwF3A/eEEG6PezwiUkNC8Dbxjo6Kt4nT3Z2/X6lwv9INN/g5eFOnlv89yymXO7q5YepUWLzYl+8aG/3PLiESEVBmdinQBjwX91hEpIYMvq124sTKtIkXUVdKtFTKlzn7+z20C1vAJ02Ke3THFHtAmVkTsAb4LPCNmIcjIrUglfLryzdt8gJ+JW6rHaquNH16bdSVhmoBX7LEZ0lTp1b3nqoRiD2ggNXAP4UQHjUzBZSIHNuhQ94m/vrrHg7NzeVfxtuxw0Np7VpvR6+FulLUAt7b67PKqAX8pJN8llSBFvBqiPVP2sxuARYDnyriuW34MiAtLS0VHpmIJEZ/v5/28Nprfpp3GU4Tf4vu7vx+pd/+1h+78EL49KeTW1fKZHyWlE57WM+cWdUW8GqILaDMbAneFPGuEELmRM8PIazGZ1u0traWclqGiNSivj5fxotOeyh3m3gmk68r/epX/o1+0SK49VavK82bV773KocQPJCiFvCJE2HBAu+6mzatpm/OPZY4Z1CXAbOAFy2f9GOBK83sC0BDCKEvrsGJSEy6u/PLeFDe0x6iulJ7u9eVDh70GceHP+xLeGeemayZRzrtoZTJeHfdzJlw+un5FvAkjbUC4gyonwIbBj3298Dv8ZlVutoDEonbsmX+82OPxTmKGPT3+2kPmzf7ct748eU97WH79vx+pc5On31ceSUsXw6XXpqculLhRtkQPIQWLvRZUmNjcsZZJbF9tSGEA8CBwsfM7AiwL4TwQhxjEpEqG3xoa5kuBQS8oSKqKz37rD920UVw441+v1JS6kp9fd4CnsCNsnGrrzgWkWQYvIw3bVp5lvEK60qPP+4fv+1tcNttXleaO3fk7zFS0V1JmQwQ/OPTTvPlu4RtlI1bogIqhLAs7jGISIVE3XhbtpR3GW+outKMGfCRj3hdaenS+Gs1g48TOukkn8GNGwtXXBHv2BIsUQElIqNQpbrxklxXimZJvb0ejoXHCU2b5qFcm1uTqkoBJSKVUbipFnxT7UiX8ZJcV4pmSf39HkBz5/qRQk1NiT5OKMkUUCJSPrkc7Nnjm2oPHizPptpMBp54Ir9fKSl1pWiWlEr5MuO0ad4CHh26WiPHCSWZAkpERq63N9+Nl8n4TKbIZbxlbWcAgy7uCwFeeCF/Dl5S6kqDa0lz5w770NW620owDAooERmeEPwA1c7O8l1xsW2b15XWrcvXlaJz8OKoKxXWkkCzpCpTQIlIaTKZ/BUXhw+P+IqLqdn9XL3/X+Dmv4d//3d/nYsugptugmuuqX5daahZkmpJsVBAiUhxDh/2zrmOjvxNtcPtxovqSu3t/PPzTzAhpGH8qXD77fD+91e3rtTf74EU1ZIaG73jLtqXpFlSbBRQInJsQ+1dam4e3mbSEOD55335rqCu9K+zP8dDMz7OvT+aXL260uDTG046yQ+H1SwpURRQIgmxZg2sX+/fOxctglWrYMWKmAYTXQi4ebMPaCRHEEV1pbVrveU8qistXw6XXMI9t57lz7PfHf91RmLwGXcNDTq9oQYooEQSYM0aaGvLX4Da0eEfQxVDqrDpYedOX9qaNs1nFaU6dAgefthDqbCu9JnPVK+ulMn4LCk6CTw64665GSZPrvz7y4gpoEQSYOVK/wd+oZ4ef7ziAZVOw+7dPls6csSXuGbNKn25raCuxL/9m398ahXrSiEcPUuaOBFaWvL3JdXZSeCjgf6LiSRAZ2dpj5fFwYO+/LZtW36jaanLeFFdae1anzEV7ldavhyWLKlsXSmb9VlS4a2yp53ms6Q6uC9ptFNAiSRAS4sv6w31eFllMn7Sw+bNvgw3YcLwDmzdts1Dad26fF1p2TLfr3TJJZWdraRSHkr9/T7+efM8WJuaRuWtsvVMASWSAKtWec2pcJlvyhR/vCy6u71FvLPTv7E3NpY+Wzp4MF9Xeu45n520tsJnPwtXX125ulLUBl64WXbJkvxmWc2SRi0FlEgCRHWmm2/2RomFC8vQxZfN5mdLBw/6rKbUFvF0On8OXjXrSum0z5Ky2Xwb+Jlnqg28ziigRBJixQq4917/9YjOaevu9qOHog21JZyLB3hd6bnnPJQeecTDbeZM+OhHfQmvjHWlNeums/75BvoyxqLlZ7PqxldZsWyHd9ktWpRvcFAbeF1SQImMBuWYLUV1pbVr/dcTJ/rS3fXXwzveUfa60pr2JtruXkhfxutfHbsm0faDc2HJElZ8VrMkUUCJ1LZDh7y29PrrPlsqtbY0VF3p4ovhc5/zcGpoKO9402mf4eVyrLznLHr6jg7Qnt4xrPz2JFZ8trxvK7VJASVSa6LDWrds8YAq9fih49WVrrvOD0YtlxC8ueHIkfwJDqeeCrNm0fnG0LOkirbWS01RQInUghB8trN9e37fUimzpcK60sMPe7BVqK5ELueBFB2LMX26h1Jz81Ezsqq11kvNUkCJJFkq5ac8bNniPeil7lvats1Pdli3rrJ1pcHHCs2dmz98dcKEIT+l4q31UvMUUCJJEwLs2etrXbt3+8ymlNlStepKhRtmJ03yqc/s2UV33VWktV5GFQWUSFJ0d0PvWEj1wTPPeKt1sWfipdNeT4rqStmsH/nzpS/5fqVy1JVCyG+YjY5GWrrUZ3RTpw5ribBsrfUyKimgROKUTnt7eNTwkLoExo4rbrYUgp8UHu1XiupKH/uYL+GdccbI60pRPSmV8teaPduvPJ8+XSeCS8UpoESqrb8f9u/3mlBXlwdNtJl2XBFnyb3+en6/0vbtvrwW1ZUuvnjkdaXCA1gL60nNzTrrTqoqtoAys9uAzwOLBh56EfiLEEJ7XGMSqajubg+kzk7/5j9p0lsaHh5bfYxL+w4c8FnS4LrSLbeUp65UsD+JCRNg/nxfFmxq0ikOEpuSAsrMfgf8HXB/CKFrhO+9Dfgz4PfAGOBG4KdmdlEI4bkRvrZIMqRSvoS3dasHwLhxxd9NdKy60pe/7HWl4d5wWzi2qMlh8mR/7ehoIR3AKglQ6gwqA/wl8F/MbC3wt8DaEEJ/qW8cQvjXQQ+tNLMvApcBCiipXdmsL+F1dvqGWrPiz8OrdF2pt9dDKdpHtXSpv35Dg0JJEqekgAohnG1mlwI3A38E/AHQZWb3Af8rhPDacAZhZmOBjwJTgSeH8xoiserv9/buHTu8LpTL+aaeYrvwOjt9r1K560rRLbM9Pf7r5mY45xxfWpwyZXivKVIlJf+tDyGsB9ab2R3Ax/Cwugv4qpn9Ep9V/d8QQt+JXsvMzgWeAiYBh4EPhRCeP8Zz24A2gBZtNZckCMFnI4V1pQkTij926MCB/H6l558/uq50zTXDD5CoHTzaATtrFixe7KGkqyqkhlgIYeQvYnYG8A3gE0AADgA/Ar4bQjjmyVpmNgFoAZqAjwC3AMtCCC8c7/1aW1vDhg0bRjxukWHp6fGlu44OD4KxY72ZoNi60q9+5aH0xBO+HLh4sc+URlJXKtyjBN4OfvLJ3g4+ceLwXrNKli3zn7UPqj6Y2cYQQmsxzx1RF9/A0twH8FnU+/FwehToA24HbjGzTw5RbwIghJAGNg18uNHMLgbuHHg9keRIpWDvwOkOBw54KJVaV2pv97pSd7fPaj7+8XxdaTgKZ0pjxngoLV3qoXSM44VEasmwAsrMluIhcgNwErAb+Gvg3qgOZWaLgX8E/iswZEANYQyQ7H/uSf1Ip2HfPt93tGdPac0O4GG2dq3XlgrrSsuX+1LecNq3hwqlM8/0ZcUaDSXNnORYSm0zvxn4LHDpwEOPAKuBfw0hZAufG0LYZGb/Ha9JDfVa3wHagdeBRuCTwDJgeSljEimrTCa/iXb37vwVEbNnF9fsMLiuNGaMh1Fbm4fTcOpKQy3f1XgoiRSj1BnUvUAX8B18trT1BM9/Ca9FDWUu8H8Gfj6It5ZfF0J4sMQxiYxM1Ba+fTvs2uUdeaV04A1VVzr9dLjjDq8rzZ5d+pgKu+9Ay3dSl0oNqA8DD4QQcsU8OYTwa+DXx/i9m0p8b5HyyWZ9trNjh3fh5XK+WbXYqyyOVVf6xCe8rnT66cMbV+E+pZkz/XVmzEh8o4NIJZS6D+qnFRqHSOUNFUqTJvmspNj7lSpRVyo80aG5Gc4918NJLeFS53RYrIxu5QilAwfgoYc8mF54oTx1pUzGT4nI5bzxYulSX8bT5lmRNymgZPTJZPKhtGvX8EKpr+/oulIuN/K6UjbrS4GZjI/ntNO8I7CxsfTXEqkDCigZHdJpD6Xt2737LrrltZRQ6u8/+hy8qK70yU8Ov67U3+/Ld6mUb+RdsCB/FbrOvhM5LgWU1K6+Pt+ntGOHn+wQgi+RFdvoEOnoyNeVduzwZonCc/CGU1fq6fFgMvNAWrCg+COQRARQQEmt6e31UNq2zVvDobSr0SNRXam9HV580QPtHe+AL3zBz94ZTi2or89DKZfzkLzgAm92UFu4yLAooCTZok2qe/d6KB065GFSyj6lyFB1pTPOgD/+Y3jf+4ZXV8rlfEyZjI9pyRI1O4iUiQJKkqe/H7q7WfbeCdDXx2N/+ZSHUinHDBW+VlRXevhhn+HMnu11peXL/aDWUhUeNxTVlebP10V/ImWmgJJkyOX8PqVdu7wOlE5D6hKv2QznhO+h6krXXON1pdbWke1XCsFD7qyzfClPdSWRilBASXyizruurvwepQkTfKY0bhyMG1/a6x04AA8+6ME0uK509dUeUqUqXMJraPBQmj1bm2hFqkABJdXV05NvB9+712cjkyeX1g5eKKortbfDk08eXVd6//u9TjUchw/7WMeOhVNO8buVGhu1hCdSRQooqawQfD9R1HnX3e1BNJzOu0h/Pzz7bH6/UjnqSuCzpIMH/fVnzfITw2fMGP516yIyIvo/T8ovm/VlscJ60tixvkQ23BtjAbZuzdeVdu4sT10pCtBUypftzjgD5sxRF55IAiigpDxSKZ99bN/ul/sNricNU1NmD/z4xx5ML73ks69LLoFbb/X9SsOpK4EvDR465L+eN8+X8Zqbh7fMKCIVoYCS4QnBl9aie5QOHPDHp0wZfj0p0tcHjz/O3Zv+nHccegSez5anrtTf76GUTvts7uyzfUanqyxEEkkBJcWLWsHfeMOX7lKp4e9PGqywrvTww6w58gG+yr1s42QWzuph1Q27WHHd/uG9dirly3hmvmfp5JN1Fp5IDVBAyfFFS3fReXe5HIwf7x1t06aN/PUH15WmTGHNGd+i7aU76El7m3nHngbaVi0EKD6kcjkPpXTax3reeT7z0rFDIjVDASVHK+y6277dfw0jawUfbP/+/H6lIepKKz968ZvhFOlJjWXlPSefOKBSqfxxSFF7uE54EKlJCijJX563a5fPYtLp/NLdcM6nG0oqlT8H7wT7lTp3DT3LOdbjR9WWGhvh/PN93ONL3OgrIomigKpXQ22YLUPX3VH6++G3v83vVzpyxGtVK1Z4a/gx9iu1zEnT0fXWxoWWOemjH4hqS2PG+ExpwQLNlkRGEQVUvYiO7Nmzx+tJPT3+jXw4p4KfyNatfrLDz372Zl3pzf1KF110wv1Kq27bTtuqhfSk8s+bMinHqtu2v3mQLH193ol3zjkeeqotiYw6CqjRLGpw2LnTb5nN5Xx2NHWq/yinffv8fqXCutKll8Jtt/l+pRLOrovqTDd/exF9GWPh3DSrPt/BitbNsBefLZ1yijrxREY5BVQNWLbMf37ssRM8MZpdRHuTurtHftbd8aRS8PjjHkpPPeUBuGQJ3Hmn36803P1KeEjd+y+zIZfz6zYmTYJTz4S5c7VvSaROKKBqXTqdv6aiq8sbHqJjhcrV4FAoqiu1t8PPf+51pTlz4FOf8iW8004b+Xtks14fyyyC8RP8RPJKBKyIJJoCqtYMdYKDmddgpk2r3N1EW7bk9yt1dXld6d3v9lC68MLyvG9Pj39t48d7A0VzM4wZCzNH/tIiUntiCygzuwv4MLAE6APWA3eFEF6Ia0yJFfohm4OXfu/1pEzGQ2mkh6+eyLHqSrffXnJd6Zj6+30GmMn4LOmii2DmTA88TZhE6lqcM6hlwA+BZwADvg08YmZnhRD2xTiu+EVXiu/f7x13+08HgodTOdvAh1LButJRosNax4zx9vCoRVxEZEBsARVCeF/hx2Z2A3AQeCfwQCyDilO0WXb3bg+ivj7/5j1lSn7DaXNzZd67vx9+8xsPpUrVlSC/PNnb61+XWsRF5DiSVINqxBd1hnkiaI0pnCV1db11s2xTU+XHcLy60kUXla8pIZv1ZbxczgPpvPN8OU8t4iJyHEkKqO8DzwJPDfWbZtYGtAG0tLRUb1TldLxZUrk3yx7Lvn35c/BeftlrPZdcUt66UqS311vdx42Dt73N9y+VcBHgCdvqRWRUS0RAmdl3gXcB7woh5IZ6TghhNbAaoLW1NVRxeMM31CwJfMmuWrMk8LrSL3/pobR+vc9kli6Fr3zF60ozy9gmF4KHcF+f15QuuMDb3XVtuoiUKPbvGmb2PeDjwNUhhM1xj2fEkjBLgnxdqb0dfvGLfF3phht8Ce/UU8v7ftEyXn+/TnoQkbKINaDM7PvAx/BweiXOsQxb4Sxp505fQqt2LanQ5s35utKuXd6Kfs01sHy571cq92bXaO/ShAm+d2n+/PIuE4pI3YpzH9Q9wA3AB4H9ZjZ34LcOhxAOxzWuoqTT+VlSV1d8s6RIVFdqb4dXXvG60qWXwpe/DFddVf7AKLzeoqnp6L1LIiJlEucM6taBn38+6PFvAd+s7lBOoPD0hp07/Weofi2p0FB1pTPPhD/5E3jve8tbV4pkMr6MF4Iv4S1YEM/XLiJ1Ic59UMkuTkRn3EWzpOgSv4aGeGZJUP26UiRaxps40TftzpunA1tFpOJib5JIjP5+/ya8b5/Pkg4c8McnTqz86Q0nsLD3Fd6778fwB2uqU1eCo5fxmpu1jCciVVffAZVK+Tfhri6fKWUy+avOK3nGXTH27oUHH2TNPxi/3Hk7P+KvWDnxP7Hqo+tZccfsyjUiFC7jLVgALS06gkhEYlFfAZXL+cbRaJZ06JAv1U2cWNmTwIuVSvnu1LVr4emnWZP7I9rs7+hhMgAdfXNpe+AP4LyONy/1K5vCbrwlS/zeJXXjiUiMRn9A9fb6jCCaJeVy+fuS4p4lgS+lbdzooVRYV/r0p1n5wA/o2TP5qKf3pMay8p6TyxNQ0abaVMqX8S680OtrcQe1iAijMaCyWZ8l7dnjs6QjR3yWVKlbZYfrtdc8lH72s3xdqfB+pTFj6Lxv6GOBOneN8HDVwrPxTj4ZFi70GaQ21YpIgtR+QIXgy1MHD3og7d3r33jHjUtGLanQnj35c/BefdVnKpddBnfcAVde+ZYltZY5aTq63tot1zInPbz3j2pu48b5CeXz53twi4gkUG0GVAgeRG+8kT9OCJI3S4K31JXI5eCss+BP/9T3K82YccxPXXXbdtpWLaQnlV9ymzIpx6rbthf//iH4jDKVgsZGnY0nIjWjNr9LHTkCv/61F/QbGpLXZTZUXWnuXPj0p30J721vK+plojrTzd9eRF/GWDg3zarbthdXf8rlfFaZzfp7L1rkdSYt44lIjajNgArBZx5JmwUMVVd6z3vguuuGvV9pxXX7ufdfZgPw2OrfnfgTCm+qXbjQT3wo4YoLEZGkSNh3+BpUYl2pYg4f9lrc5Ml+U+2cOfmbeEVEapACajgK60rr1/uSXpF1pbIqPO1h5kw4+2x/7yTV4EREhkkBVaxc7ui6Uk+P13ZuvNGPHFq0qHpjyWb9KKbo0NaWFm+AEBEZRRRQJ7JpU76utHu315WuvdabHd7+9urOVqIr1MeP16GtIjLqKaCGsmePB9LatfC733ld6fLL4c474YorqnsEUAiQy/pynpmH4uzZOu1BREY9BVSkt/fo/UqFdaX3vc/3V1VTYZv4uCUeipdfrjZxEakb9R1QUV2pvR0efdTrSvPmwU03+RJeNetKkXTa60tjx+bbxKcOtIkrm0SkjtRnQA2uK02d6t1311/vJy3E0QV35Ij/mDTJ28TnzlWbuIjUtfoJqGPVlb7yFa8rxdFsEIIv46XTvoR41llqExcRGTC6AyppdaVIdJp4f//Rp4mLiMibRl9A5XKwYYOHUlLqSpFUytvEx47VaeIiIicwegJq0yZvdnjwweTUlSLRMUQNDXDeeX4FSNLOERQRSZja/i6ZxLpSpL/fl/EyGd+3dM45Xl9Sm7iISFEshBD3GErWOm9e2LBwITzzjAfB2Wf7TOm9742vrhTRMUQiIsdkZhtDCK3FPLc2Z1BdXT4T+cxn/CqLOOtKkei22vHj4fTTvflBxxCJiAxbbQbUKafAT37iFxbGKQSvL/X25m+rPekkHUMkIlIGtRlQkyfH2/QwuL503nm+tKj6kohI2cQaUGZ2JfCnwEXAfOAzIYT74hzTcWUyHkyQry9NnRrvmERERqm4Z1BTgReA/z3wI5mi+tKECbrmQkSkSmINqBDCWmAtgJndF+dY3mJwfUnXXIiIVFXcM6iimVkb0AbQMmdO5d6ov9/bxLNZb3g4/3xoblZ9SUSkymomoEIIq4HVAK1LlpR/81ZhfamlZeCaC9WXRETiUjMBVTGqL4mIJFJ9BlRUX+rp8VPEVV8SEUmc+gqowv1Lqi+JiCRa3PugpgKLBz4cA7SY2QXAvhBCZ9neKDofD7R/SUSkRsQ9g2oFHi34+FsDP+4Hbhrxq6dSPmNSfUlEpObEvQ/qMaD862tRfSnav6Tz8UREak7cM6jyGXw+3rnn6nw8EZEaVvsBlc16MPX35/cv6f4lEZGaV7sB1dcH+/b5/UuLF+v+JRGRUaZ2AyqX8/uXZs+GcbX7ZYiIyNBq8zv71KmwbJnqSyIio1iMt/6NgJnCSURklKvNgBIRkVFPASUiIomkgBIRkURSQImISCIpoEREJJEUUCIikkgKKBERSSQFlIiIJJICSkREEkkBJSIiiaSAEhGRRFJAiYhIIimgREQkkRRQIiKSSAooERFJJAWUiIgkkgJKREQSSQElIiKJpIASEZFEij2gzOxWM9tiZikz22hmV8Q9JhERiV+sAWVmHwO+D9wNvB14ElhnZi1xjktEROIX9wzqK8B9IYR7QwgvhxC+BOwEvhjzuEREJGaxBZSZTQAuAh4a9FsPAZdXf0QiIpIk42J871nAWGDXoMd3Ae8Z/GQzawPaBj7sM7MXKju8xJkF7Il7EDGox69bX3P9qMeve0mxT4wzoEoSQlgNrAYwsw0hhNaYh1RV9fg1Q31+3fqa60c9ft1mtqHY58ZZg9oD5IA5gx6fA3RVfzgiIpIksQVUCCENbASuHfRb1+LdfCIiUsfiXuL7LvAjM/s18ATwBWA+8D9O8HmrKz2wBKrHrxnq8+vW11w/6vHrLvprthBCJQdy4gGY3Qr8Z2Ae8AJwZwjh8VgHJSIisYs9oERERIYS90ZdERGRISmgREQkkWouoOrtcFkzu9LM/p+ZbTezYGY3xT2mSjKzu8zsGTM7ZGZvmNkDZnZO3OOqNDO7zcyeG/i6D5nZU2a2PO5xVdPAf/tgZn8T91gqxcy+OfA1Fv6oi201ZjbPzO4f+P86ZWYvmdlVx/ucmgqoOj1cdirePHIH0BvzWKphGfBD/Lira4As8IiZzYhzUFWwDfgz4EKgFfgF8FMzOy/WUVWJmV2KnxTzXNxjqYJX8aaw6Me58Q6n8sysGe/UNmA5cCbwJWD3cT+vlpokzOxp4LkQwi0Fj/0e+KcQwl3xjaw6zOwwcHsI4b64x1ItZjYVOAh8MITwQNzjqSYz2wfcFUL4n3GPpZLMrAn4DfA54BvACyGE2+MdVWWY2TeBj4QQRv2qQCEzuxu4KoTwzlI+r2ZmUDpctm414n9P98c9kGoxs7Fm9nF89lwPm9ZX4//IfDTugVTJqWa2Y6BU8WMzOzXuAVXBB4GnzewnZrbbzJ41s9vNzI73STUTUBz/cNm51R+OVMn3gWeBp2IeR8WZ2bkDs+Q+fLP6h0IIz8c8rIoys1uAxcCfxz2WKnkauAl4P3AL/r3rSTObGeegquBU4FZgM/A+/P/r7wC3He+T4j5JQuSYzOy7wLuAd4UQcnGPpwpeBS4AmoCPAPeb2bIQwqg8ud/MluD15HeFEDJxj6caQgjrCj82s/X4N+0b8ZN1RqsxwIaCUsxvzex0PKCO2RRTSzMoHS5bR8zse8AngGtCCJvjHk81hBDSIYRNIYSNA/8jPwvcGfOwKukyfGXkRTPLmlkWuAq4deDjifEOr/JCCIeBF4HT4x5Lhe0EXhr02MvAcRvcaiagdLhs/TCz75MPp1fiHk+MxgCj+Zv0T/EOtgsKfmwAfjzw63Qso6oiM5sELMW/gY9mT/DWe6DOADqO90m1tsQ33MNla9ZAF9vigQ/HAC1mdgGwL4TQGdvAKsTM7gFuwIuq+80sqi8eHvjX5qhkZt8B2oHX8caQT+It96N2L1QI4QBwoPAxMzuC/90ercuafw08AHQCJwFfBxqA++McVxV8D6+1rQR+gm8T+jLwteN9Uk21mUP9HS5rZsuAobqb7g8h3FTVwVSBmR3rL+S3QgjfrOZYqsnM7gOuxovmB/H9QH8VQngwznFVm5k9xuhuM/8xcCW+tPkGsB74eghh8PLXqDOw8fxufCbVideefhCOE0I1F1AiIlIfaqYGJSIi9UUBJSIiiaSAEhGRRFJAiYhIIimgREQkkRRQIiKSSAooERFJJAWUiIgkkgJKREQSSQElIiKJpIASiZGZjTOzJ8zsiJktHfR7bWYWzOzbcY1PJE46i08kZma2EL/7qQO4JITQZ2ZnA8/gV8wsq5MLG0WOohmUSMxCCB3AzcD5wH8zs8n4lQQpYIXCSeqVZlAiCWFmPwS+iF/AeTnwH0MI/xzvqETio4ASSYiB21VfAE4D7g0htMU8JJFYaYlPJDnOB1oGfn2OmdXajdciZaWAEkkAM5sG/AOwB1gJXAZ8K9ZBicRM/0ITSYbVwELg2hDCL8zs7cBXzeyREMKjMY9NJBaqQYnEzMxuBv4WuDuEsHLgsWa89Xw8cF4IYW9sAxSJiQJKJEYDm3M34mF0VQghW/B7lwGPA+tCCB+IZ4Qi8VFAiYhIIqlJQkREEkkBJSIiiaSAEhGRRFJAiYhIIimgREQkkRRQIiKSSAooERFJJAWUiIgkkgJKREQS6f8DXOXOrmhu4mkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plt.xlabel(\"x\", fontsize=18)\n",
    "plt.ylabel(\"y\", fontsize=18)\n",
    "plt.errorbar(xd, yd, yerr=yd_err, fmt=\"bo\")\n",
    "ax.fill_between(xf, yf - sigma_yf, yf + sigma_yf, alpha=0.2, color='red')\n",
    "plt.xlim(0., 6.)\n",
    "plt.ylim(0., 6.)\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.plot(xf, yf, color=\"red\")\n",
    "plt.tight_layout()\n",
    "# plt.savefig(\"basic_chi2_fit_plot3.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}