{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Error Propagation using Random Gaussian Numbers\n",
    "\n",
    "Example calculation of propagating uncertainties, both when adding and multiplying number, and also in the general case. The propagation can be done both analytically (using the error propagation formula) and also using simulation.\n",
    "\n",
    "The example is based on FIRST doing the error propagation analytically, and then verifying it by running a so-called Monte-Carlo (MC) program, which uses random numbers for propagating errors.\n",
    "\n",
    "## References:\n",
    "- Barlow: page 48-61\n",
    "- Bevington: page 36-48\n",
    "\n",
    "## Author(s), contact(s), and dates:\n",
    "- Author: Troels C. Petersen (NBI)\n",
    "- Email:  petersen@nbi.dk\n",
    "- Date:   12th of November 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "DO THE FOLLOWING ANALYTICAL EXERCISE FIRST!!!\n",
    "\n",
    "1. A class of students estimate by eye, that the length of the table in Auditorium A is $L = (3.5\\pm 0.4)$m, and that the width is $W = (0.8\\pm 0.2)$m.\n",
    "\n",
    "   Assuming that there is no correlation between these two measurements, calculate ANALYTICALLY what the Perimeter (P), area (A), and diagonal (D) length is including (propagated) uncertainties. Repeat the calculation, given that the correlation between length and width is $\\rho(L,W) = 0.5$ - not an unreasonable number, given that they are estimated by the same (uncertain) scale.\n",
    "   \n",
    "NOTE: This is a complete standard problem, that you will be asked to solve again and again in the course. For this reason, make sure that you understand how to do it, and become good at doing it reasonably fast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining the parameters:\n",
    "mu1   =  3.5\n",
    "sig1  =  0.4\n",
    "mu2   =  0.8\n",
    "sig2  =  0.2\n",
    "rho12 =  0.5           # Correlation parameter!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not (-1.0 <= rho12 <= 1.0): \n",
    "    raise ValueError(f\"Correlation factor not in interval [-1,1], as it is {rho12:6.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note on analytic solutions:\n",
    "\n",
    "Python includes symbolic algebra in the package *SymPy*, which seems to be both simple and powerful (and in Python). In addition, printing with Latex can also be included (see below), which (in combination) is very nice.\n",
    "\n",
    "Below is a SymPy and Latex example with the hope that it will wet your appetite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.display import Latex\n",
    "\n",
    "def lprint(*args,**kwargs):\n",
    "    \"\"\"Pretty print arguments as LaTeX using IPython display system \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    args : tuple \n",
    "        What to print (in LaTeX math mode)\n",
    "    kwargs : dict \n",
    "        optional keywords to pass to `display` \n",
    "    \"\"\"\n",
    "    display(Latex('$$'+' '.join(args)+'$$'),**kwargs)\n",
    "    \n",
    "def myDiff(formula):\n",
    "    return sqrt((formula.diff(L) * dL)**2 + (formula.diff(W) * dW)**2)\n",
    "\n",
    "def myDiffWithCorr(formula, name = \"\", printNow = False):\n",
    "    dd = sqrt((formula.diff(L) * dL)**2 + (formula.diff(W) * dW)**2 + 2*(formula.diff(L)*formula.diff(W)*(sigCorr**2)))\n",
    "    if(printNow):\n",
    "        lprint(latex(Eq(symbols('sigma_'+name), dd)))\n",
    "    fd = lambdify((L,dL,W,dW,sigCorr),dd)\n",
    "    return dd, fd\n",
    "    \n",
    "def diff_and_print(formula, name = \"\"):\n",
    "    # Calculate uncertainty and print original relation/formula and the uncertainty\n",
    "    dd = myDiff(formula)\n",
    "    lprint(latex(Eq(symbols(name),formula)))\n",
    "    lprint(latex(Eq(symbols('sigma_'+name), dd)))\n",
    "    \n",
    "def lambdifyFormula(formula, *args, name = \"\"):\n",
    "    # Turn expression into numerical functions \n",
    "    f = lambdify((L,W),formula)\n",
    "    d = myDiff(formula)\n",
    "    fd = lambdify((L,dL,W,dW),d)\n",
    "    return f, fd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$P = 2 L + 2 W$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\sigma_{P} = \\sqrt{4 \\sigma_{L}^{2} + 4 \\sigma_{W}^{2}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$A = L W$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\sigma_{A} = \\sqrt{L^{2} \\sigma_{W}^{2} + W^{2} \\sigma_{L}^{2}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$D = \\sqrt{L^{2} + W^{2}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\sigma_{D} = \\sqrt{\\frac{L^{2} \\sigma_{L}^{2}}{L^{2} + W^{2}} + \\frac{W^{2} \\sigma_{W}^{2}}{L^{2} + W^{2}}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$C = (8.6 \\pm 0.9)\\,\\mathrm{m}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$A = (2.8 \\pm 0.8)\\,\\mathrm{m}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$D = (3.6 \\pm 0.4)\\,\\mathrm{m}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\sigma_{P} = \\sqrt{4 \\sigma_{L}^{2} + 8 \\sigma_{LW}^{2} + 4 \\sigma_{W}^{2}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\sigma_{A} = \\sqrt{L^{2} \\sigma_{W}^{2} + 2 L W \\sigma_{LW}^{2} + W^{2} \\sigma_{L}^{2}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\sigma_{D} = \\sqrt{\\frac{L^{2} \\sigma_{L}^{2}}{L^{2} + W^{2}} + \\frac{2 L W \\sigma_{LW}^{2}}{L^{2} + W^{2}} + \\frac{W^{2} \\sigma_{W}^{2}}{L^{2} + W^{2}}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$P = (8.6 \\pm 1.1)\\,\\mathrm{m}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$A = (2.8 \\pm 0.9)\\,\\mathrm{m}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$D = (3.6 \\pm 0.4)\\,\\mathrm{m}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Import SymPy: \n",
    "from sympy import * \n",
    "    \n",
    "# Define variables:\n",
    "L,W,P,A,D = symbols(\"L, W, P, A, D\")\n",
    "dL,dW,dP,dA,dD = symbols(\"sigma_L, sigma_W, sigma_P, sigma_A, sigma_D\")\n",
    "\n",
    "# Define relations:\n",
    "# Perimeter:\n",
    "P = 2*L + 2*W\n",
    "# Area:\n",
    "A = L*W\n",
    "# Diagonal\n",
    "D = sqrt(L**2 + W**2)\n",
    "\n",
    "# Try writing a simple function to not repeat yourself! (See cell above)\n",
    "diff_and_print(P,\"P\")\n",
    "diff_and_print(A,\"A\")\n",
    "diff_and_print(D,\"D\")\n",
    "\n",
    "dP = myDiff(P)\n",
    "dA = myDiff(A)\n",
    "dD = myDiff(D)\n",
    "\n",
    "# Turn expressions into numerical functions \n",
    "fP, fdP = lambdifyFormula(P,\"P\")\n",
    "fA, fdA = lambdifyFormula(A,\"A\")\n",
    "fD, fdD = lambdifyFormula(D,\"D\")\n",
    "\n",
    "# Define values and their errors\n",
    "vL, vdL = mu1,sig1\n",
    "vW, vdW = mu2,sig2\n",
    "\n",
    "# Numerically evaluate expressions and print \n",
    "vP = fP(vL,vW)\n",
    "vdP = fdP(vL,vdL,vW,vdW)\n",
    "vA = fA(vL,vW)\n",
    "vdA = fdA(vL,vdL,vW,vdW)\n",
    "vD = fD(vL,vW)\n",
    "vdD = fdD(vL,vdL,vW,vdW)\n",
    "\n",
    "lprint(fr'C = ({vP:.1f} \\pm {vdP:.1f})\\,\\mathrm{{m}}')\n",
    "lprint(fr'A = ({vA:.1f} \\pm {vdA:.1f})\\,\\mathrm{{m}}')\n",
    "lprint(fr'D = ({vD:.1f} \\pm {vdD:.1f})\\,\\mathrm{{m}}')\n",
    "\n",
    "#Adding correlations (and also derivation, printing and lambdifying)\n",
    "sigCorr = symbols(\"sigma_LW\")\n",
    "rho = symbols(\"rho_LW\")\n",
    "\n",
    "dP, fdP = myDiffWithCorr(P, \"P\", True)\n",
    "dA, fdA = myDiffWithCorr(A, \"A\", True)\n",
    "dD, fdD = myDiffWithCorr(D, \"D\", True)\n",
    "\n",
    "sCorr = sqrt(rho*dL*dW)\n",
    "fSC = lambdify((rho,dL,dW),sCorr)\n",
    "\n",
    "vSigmaCorr = fSC(rho12,vdL,vdW)\n",
    "\n",
    "# Numerically evaluate expressions and print \n",
    "vdP = fdP(vL,vdL,vW,vdW,vSigmaCorr)\n",
    "vdA = fdA(vL,vdL,vW,vdW,vSigmaCorr)\n",
    "vdD = fdD(vL,vdL,vW,vdW,vSigmaCorr)\n",
    "\n",
    "lprint(fr'P = ({vP:.1f} \\pm {vdP:.1f})\\,\\mathrm{{m}}')\n",
    "lprint(fr'A = ({vA:.1f} \\pm {vdA:.1f})\\,\\mathrm{{m}}')\n",
    "lprint(fr'D = ({vD:.1f} \\pm {vdD:.1f})\\,\\mathrm{{m}}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # Modules to see files and folders in directories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('../../../External_Functions')\n",
    "from ExternalFunctions import Chi2Regression\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error propagation - Simulation\n",
    "\n",
    "Now we want to try to see, if we can solve the above error propagation problem using simulation. The method is relatively straight forward: You simply take \"realistic\" values of the input parameters x (here Length (x1) and Width (x2)), calculate the resulting value y (here Perimeter, Area, and Diagonal), and do this many times. The resulting distribution of y should be centered around the value y(x1,x2), and the standard deviation should reflect the uncertainty in y from the uncertainties in the input variables.\n",
    "\n",
    "This is a much more clumsy of calculating the uncertainty, but comes with the advantage, that if the resulting uncertainty is not Gaussian, then one can actually see this (i.e. it avoids the assumptions used in the usual error propagation formula)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we set the parameters of the program:\n",
    "N_exp = 10000           # Number of \"experiments\" (i.e. drawing from random distributions)\n",
    "save_plots = False\n",
    "r = np.random\n",
    "r.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define parameters for two random numbers (Gaussianly distributed):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1   =  3.5\n",
    "sig1  =  0.4\n",
    "mu2   =  0.8\n",
    "sig2  =  0.2\n",
    "rho12 =  0.5           # Correlation parameter!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not (-1.0 <= rho12 <= 1.0): \n",
    "    raise ValueError(f\"Correlation factor not in interval [-1,1], as it is {rho12:6.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we calculate numbers that allows the transform from uncorrelated variables `u` and `v` to correlated random numbers `x1` and `x2` below (see Barlow page 42-44 for method):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the parameters needed for the transformation:\n",
    "theta = 0.5 * np.arctan( 2.0 * rho12 * sig1 * sig2 / ( np.square(sig1) - np.square(sig2) ) )\n",
    "sigu = np.sqrt( np.abs( (((sig1*np.cos(theta))**2) - (sig2*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )\n",
    "sigv = np.sqrt( np.abs( (((sig2*np.cos(theta))**2) - (sig1*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )\n",
    "\n",
    "# Produce random numbers with the (possible) correlation:\n",
    "u = r.normal(0.0, sigu, N_exp)\n",
    "v = r.normal(0.0, sigv, N_exp)\n",
    "x1_all = mu1 + np.cos(theta)*u - np.sin(theta)*v\n",
    "x2_all = mu2 + np.sin(theta)*u + np.cos(theta)*v\n",
    "x12_all = np.array([x1_all, x2_all])\n",
    "\n",
    "#y_all = 2*x1_all + 2*x2_all\n",
    "#y_all = x1_all * x2_all\n",
    "y_all = np.sqrt(x1_all**2 + x2_all**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the above is nothing more than a matrix multiplication written out! Also note that the absolute value is taken before the square root to avoid `np.sqrt(x)` with `x<0`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "Plot both input distribution and resulting 2D-histogram on screen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4kAAAIzCAYAAACtCFzyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3yV5f3/8dcnIQSQJCyFyBBUFBBQMFH6cxCBlIQhijgQyyiFKqKtgIItBUTE2GrFylBQv3EwZYlgkWVwoQhKK8MByFAQR4GwhIzr98c5HBMySE7COcnh/ezjPMh939d9XZ/7DtR8ci1zziEiIiIiIiICEBbsAERERERERKTsUJIoIiIiIiIiPkoSRURERERExEdJooiIiIiIiPgoSRQREREREREfJYkiIiIiIiLioyRRRM5aZrbJzBKCHUcwmdnNZrbbzA6bWat8rjszuzgIcSWY2belVNe5ZvalmVUqQtl/mtndpdFuUZhZmpn9oYR1PGdmfyvkeqHfQzPbYWYdShJDjrpmmtlNRSjX0sw+LI02RUTkV2Z2qZltyPFJN7M/m1kNM1tuZl97/6xeWD1KEkUkJOX3g6+Z9TWz908eO+cuc86lnaaeht4fsiucoVCD7UlgsHOuqnPus2AFcYaT0RHA/znnfilC2X8AfzWzimcollLnnLvbOfdoUcqaWaqZjTsTcZhZS+By4I3TlXXO/Rc4YGZdz0QsIiJnK+fcl865K5xzVwBXAkeBBXj+W7jSOdcYWOk9LpCSRBGRICoDyecFwKYgx3DGmFkk0Ad4rSjlnXN7gS+AG89kXCHqj8B055wrYvnp3ntEROTMaA9sc87tBLoBL3vPvwwUOupDSaKInLVy9jaa2VVmts47LGOfmf3TW+xd758HvEMyf2NmYWY20sx2mtkPZvaKmcXkqLe399rPZva3U9oZY2Zzzew1M0sH+nrbXmNmB8xsr5lNzNmT5e1lG+QdInLIzB41s4u896Sb2ZyCer4KitXMIs3sMBAO/MfMthXhfUWa2ZNmtsv7jp4zs8reawlm9q2ZDfW2s9fM+uW4t6aZvemN9xMzG3eyV9fMTr7j/3jf8e057iuovk5mttn7Pr4zs2EFhH01cMA59633vhreOLt6j6ua2VYz653jnjSgcyHv4XUz+97MDprZu2Z2WY5rqWY2ycyWeGP72MwuynE90cy+8N47EbAC2qhkZsfMrJb3eKSZZZpZtPd4nJlNyNHmuBz3Puh9X3vM7Pc5zg8EegEPed/zmzmavMLM/uuNa7Z5h+aaWS0zW+z9u/k/M3vPzAr62SEZWO29L9JbvkWO9s/zPtO5Od5ze/Mk8iIiUvruAGZ6v67t/UXoyV+InlfYjcH+DbaISFnxDPCMc+5VM6sKNPeevx74BqjmnMsE8P7g3Re4AfgBeAWYCPzOzJoBk4EkYC0wHqh7SlvdgFuB3kAk0Ax4AFgH1AP+DQwCJuS4JwnPsJH6wKfA/8PzA//PwBqgJ7/+hjCnvvnF6pz7HVDVzBxwuXNuaxHe0RPAhcAVQAYwAxgFPOy9XgeI8T5vIjDXzBY65/YDk4Aj3jINgbeBnQDOuetPjcM8c0ULq+9F4Dbn3HvmmVfRqICYWwBfnjxwzv3P+/17xTzDIx8DNjjnXslxzxbglkLew7+B3wMnvO9kuvednNQTz/frUzzfk8eAO7wJ3zzvvW8Ag4G7gVdPbcA594uZfQK09d5zvfd9XeNt/3rg6VPvM7MkYBie3x5/A0zLUedUM/t/wLfOuZGn3HqbN+ZfgA/w/J15DhgKfAucTOzaAHl6Cs3sHDzfgy+9bR03s1nAXcDwHO9lhXPuR2+Z78wsA7gU+O+pdYqIlBcdbzjH/fy/rIC0tf6/xzfh+f/qk6Y656aeWs77y+Mb+fW/0cWiJFFEQtlCM8vMcVwRzw/u+ckALjazWs65n4CPCqm3F/BP59x2ADN7GNjo7enqAbzpnDvZSzYKuP+U+9c45xZ6vz4GrM9xbYeZPY8nOciZJD7hnEsHNpnZRmBZjvb/DbQi/ySxwFhPJr1FYWYGDABaOuf+5z03Hk+iePI/QBnAWG+9b3l7Ki/1Jju3AM2dc0eBzWb2MpBwmmbzrQ/P9yYDaGZm//EmjfsLqKMacCjnCefcMjN7Hc+cjJp4EsmcDnnvy5dz7qWTX5vZGGC/mcU45w56T893zq31Xp8OnOyV7gRsds7N9V6bgCcJK8hqoK2ZvQG0BB73Hr8DxAPv5XPPbXjmX27MEV/PQto46V/OuT3ee97k16Q3A4gFLvAm8Pm1Cb++r5zv+mU8if3Dzrls4HfA30+5r9B3LSJSHvz8vyzWvt0gIG2Fx379i3MurghFk4FPnXP7vMf7zCzWObfXzGLx/OK4QBpuKiKh7CbnXLWTHzy9cwXpD1wCfOEdDtmlkLLn4+0F89qJ55dutb3Xdp+84E2Kfj7l/t05D8zsEu+Qvu/NMwR1PFDrlHv25fj6WD7HVf2ItTjOBaoA671DDw8AS/m1hwng51MSz6PeuM71tpnzuXO9gwIUVB94ks5OwE4zW21mvymgjv1AVD7np+LpLf4/59yp358o4EB+lZlZuJmlmNk27/dqh/dSzu/X9wXEfOrfDUfh72E1nkS6NfA5sBzPLw/aAFu9v8w4Va42yP29L0xBMf8D2AosM7PtZlbQQgcn35fvXTvnPsbTe9zWzJoAFwOLTrmvwHctIiIl0pNfh5qC5/9/+3i/7sNpFhlTkigiAjjnvnbO9cQzRv8JPD0g55DP0DpgD54FX05qAGTiSdz24hkyCoB55uzVPLW5U46n4FkspbFzLhr4CwXMVfNDYbEWx094ktHLciTeMc65gpLTnH70tlkvx7n6xWw/F+fcJ865bni+XwuBOQUU/S+e5N/HzMKB5/EMvb3H8q6q2hT4TwH13YlnuHAHPENhG56stghh7yXHc3t7Zwt7Dx/i6Tm9GVjtnNuM5/vXGe/cv9O14S2fU1EXlfEUdu6Qc26oc+5CoCswxMza51PuCLCNU941nt7Eu/D0Is7NucKsmZ2Pp3f/S0REyjEHZAfof0VhZlXwTNOYn+N0CpBoZl97r6UUVoeSRBERwMzuMrNzvcPiTvZsZOFJcLLxzMU7aSbwgJk18s5fHA/M9vZ6zQW6mtn/884HeITTJxBRQDpw2Nvjck+pPVjhsRaZ971MA542s/MAzKyumXUswr1ZeP5DNcbMqnifsfcpxfaR+x0XyMwqmlkv7xDPDDzvrqDJIGuBamaWc17oX7x//h7PFiCveBPHk9rimfeXnyjgOJ7e4Sp43mdRLQEuM7Pu5lnV9n488y7z5e2FXg/cy69J4Yd4VgQtKEmcg2cxpGbeHxJGn3K9yO8ZwMy6mNnF3oT25Hsu6F2/hefd5fQqniT3LjxJeU4JwCrn3PGixiMiIqfnnDvqnKuZYxoEzrmfnXPtnXONvX/+r7A6lCSKiHgk4ZnvdxjPIjZ3OOd+8f6g/hjwgXeYZRvgJTw//L6LZ3GQX4D7AJxzm7xfz8LTq3MIz7j/wn4QHoanh+oQnkRsdik+V4Gx+mE4nqGHH3mHWq7A09NVFIPx9Lx9741nJrnfyRjgZe87vq0I9f0Oz/zNdDyLv9yVXyHn3Akg9eR1M7sSGAL09iavT+D5JfAI7/VYPAsJLcyvPjyJzk7gO2Azhc9dPTWWn/AsWJSCJ8lsjGeRmMKsBiLwJLsnj6P4ddXdU9v4N565rKvwfK9WnVLkRTxzOQ+YWUHPmFNjPN/nw3gWSJpcyN6iU4Fe3oTyZDzf4pkH7Mg7n7EXnsVxRETKOUeWyw7IJ1Cs6NsZiYhIcXl77w7gGUr6TbDjKSvM7AmgjnOuz2kLl7ytc/EkKK2cc8dOU/YpPHtKTT7TcYUiM5sBzMmxMBNm9hKwJ+eKqubZGmOqc66guaQiIuXGlZdHug+XnrqQ+ZlR6fxv1hdx4ZoS0eqmIiKlzDx78K3EM8z0STyLjuwIZkzB5h1iWhHPu4jHs1DQHwLRtnfLhSZFLFvYaqNyGs65O3Mem1lDoDue1XdzlvscUIIoIiHBMycxtDregjbc1DwbBa81s/+Y2SYzeySfMmZm/zLPRsf/NbPWwYhVRKSYuuFZMGYPnuF6dzgN24jCMy/xCJ55c09xmpXVpHwzs0eBjcA/1IsuIlK+BG24qXfOwjnOucNmFgG8D/zJOfdRjjKd8Myd6QRcjWej66uDErCIiIiIiMgpWl8e6d5bWuA6ZKWq6vm7AjLcNGg9ic7jsPcwwvs5NWPtBrziLfsRnhXqYgMZp4iIiIiIyNkkqHMSvUuOr8ezwe4k78a7OdUl96bA33rP7Q1MhCIiIiIiIgVzOLJCbFZJUJNE7/LjV5hZNWCBmTV3zm3MUSS/vcXy/Q6Y2UBgIMA555xzZZMmRVqjQEREREREgmz9+vU/OefODXYc4lEmVjd1zh0wszQ8+5TlTBK/BernOK6HZyGI/OqYimePJuLi4ty6devOTLAiIiIiIlKqzGxnsGMoCa1uWkrM7FxvDyJmVhnoAHxxSrFFQG/vKqdtgIPOOQ01FREREREROUOC2ZMYC7zsnZcYhmfz3cVmdjeAc+454C08K5tuBY4C/YIVrIiIiIiIyKkckBViPYlBSxKdc//llM11veefy/G1A+4NZFwiIiIiIiJnszIxJ1FERERERKS80pxEERERERERCVnqSRQREREREfGTg5DbJ1E9iSIiIiIiIuKjnkQREREREZESyA52AKVMPYkiIiIiIiLio55EERERERERPzlcyO2TqJ5EERERERER8VFPooiIiIiIiL8cZIVWR6J6EkVERERERORXShLzMWbMGMwsz6datWp+1TdhwgTS0tJKHNeOHTsws1Kp62wzZ84cEhMTqV27NtWqVSMhIYEPP/ww2GGJiIiISDnn8KxuGohPoGi4aQFiYmJYunRprnMVKvj3uiZMmEDfvn1JSEgoUUyxsbGsWbOGZs2alaies9HEiRO58MILuffee6lSpQrPPfccN9xwA+vWraNFixbBDk9EREREpMxQkliAChUq0KZNm2CHkUtkZGSZi6m8WLBgATVr1vQdJyQkUL9+fSZPnsyUKVOCGJmIiIiIlG9GFhbsIEqVhpv66WTP4NSpU6lfvz41a9ZkyJAhZGf/2hF8cpjqzp07eeSRR3zHffv2zVVXamoqZsaXX35JUlISVapUITY2ljfeeAOADRs25Br2WtBw0xUrVtCmTRsqV67M+eefz/jx4/OUmTZtGk2aNKFy5crUrVuXO+64g+PHj5faeymrciaIABUrVqRRo0bs2LEjOAGJiIiIiJRR6kksRGZmZq7jsLAwwsJ+zas3btzIwoULeeGFF0hLSyMlJYUbbriBrl27ArBmzRoAbr75Zjp37swf/vAHAM4999x827vrrrvo2LEjw4YN46uvvuLYsWMAXHLJJaxZs4a9e/fSvXv3fO9NS0sjOTmZHj168Mgjj/DFF1/w8MMPc+655zJgwAAA3n33XQYOHMjw4cNJSkpi3759vP766xw/fpzIyMgSvKny5+jRo2zZssX3bkRERERE/OGA7BBb3VRJYgF+/vlnIiIicp3r2LFjrnmKx44dY+bMmcTExNCxY0dmzJjBypUrfUniyaGhkZGR1KtX77RDRTt16sQjjzwCQIcOHXznq1SpQps2bQrt9Xr44Ydp3bo1M2bMwMzo2LEj+/btIyUlxZcIrV27lkqVKpGSkuK77/bbby/C2wg9KSkpnDhxgsGDBwc7FBERERGRMkVJYgFiYmJYsWJFnnM5NWnSJNe5Bg0asG/fPr/bvO222/y678iRI6xdu5bx48eTlZXlOx8fH8/jjz/OoUOHiIqKonnz5vzyyy8MHDiQXr16cdVVV1G5cmW/4y2v0tLSePzxx3nmmWdo2LBhsMMRERERkXJOcxLPEhUqVCAuLi7Xp3HjxrnKREVF5ToODw8nIyPD7zbr1q3r13379+8nOzubESNGEBER4fucHJr63XffAZCUlMRzzz3Hxx9/TLt27ahRowaDBg3KlViGum3bttGjRw/69evHoEGDgh2OiIiIiEiZo57EMsTfLTaqV6+OmTF27FiSkpLyXM/ZW/bHP/6RP/7xj/z8888888wzPProoyQkJPjdi1meHDhwgC5dutCqVSsmTZoU7HBEREREJAQ4Qq8nUUliAERFRXHkyJEzVv8555xDfHw833zzDXFxcUW6p2bNmjzyyCM88cQTvp7GUJaRkcEtt9xChQoVmDt3bp75piIiIiIi4qEksQCZmZl89NFHec7HxcUVu8evRYsWzJ07l06dOhEbG0tMTAyxsbFFvn/z5s2kp6ezd+9e33GlSpWAXxfHGT9+PMnJyURGRvoWzvnss89Yt24d8+fPB+DJJ59k+/btdOjQgerVqzNjxgyysrJISEgo1vOUR/fccw8ffPABr776Klu2bPGdj46OplmzZkGMTERERETKu2ynnsSzwsGDB/nNb36T5/yPP/5IrVq1ilXXuHHjGDBgAN26dePQoUP06dOH1NTUIt8/aNAgVq9e7Tu+9957fV8751lvt3379ixdupQxY8aQmppKxYoVad68OX369PGVbdWqFUuXLmXOnDkcO3aMSy+9lNdff51WrVoV63nKoxUrVnD8+PE8w2rbtm1b4L6TIiIiIiJnIzuZZISSuLg4t27dumCHISIiIiIiRWBm651zRZs3VcZc1rKim7G4dkDauuKCbwPyntSTKCIiIiIi4qdQXLhGW2CIiIiIiIiIj3oSRURERERE/OQwskKs7y20nkZERERERERKRD2JIiIiIiIiJRBqW2CoJ/EUqampmBmtW7fOdf7aa6/FzJg4cWKQIitfXn75Za6++mqqV69OdHQ01157LatWrSpRnVu2bCEiIoK4uNwLOmVnZ/O3v/2N888/n0qVKhEfH59ry5CTpk2bRrNmzahSpQqXXHIJTz/9dIniEREREREJRUoSC7B9+3Z2794NePZG3Lx5c5AjKl/27dtHUlISqampzJ8/n8aNG5OUlMT69ev9rnPo0KFERUXlOf/Pf/6Tf/zjHzz44IMsWLCAmJgYOnXqxPbt231lXn/9dQYOHEj37t1ZvHgxvXv3ZujQoUybNs3veERERERETq5uGohPoChJLEDHjh158803AVi8eDGJiYlBjqh8eeihh3jkkUfo1q0bHTp04KWXXqJu3bpMnz7dr/refvttdu3aRZcuXfJce/rppxk4cCAPPPAAycnJzJs3j7CwMKZMmeIrM3fuXK6++mrGjRtHu3btGDlyJO3atWPBggV+P6OIiIiISChSkliAG2+8kUWLFgGwaNEibrzxxjxlVqxYQZs2bahcuTLnn38+48ePz3X9888/59Zbb6Vu3bpERkZy6aWXMmnSpFxl+vbtS0JCAlOnTqV+/frUrFmTIUOGkJ2dXeRYMzMziYuLo0ePHrnO33777bRs2ZITJ04Uua4zxcyIiooiPT292PdmZWUxdOhQHnvsMcLCcv+VPXToEHv27OHqq6/2nYuJiaFVq1akpaX5zmVmZhIdHZ3r3qioKJxzxY5HRERERORXRpYLC8gnUJQkFiAxMZGPP/6YH3/8kVWrVpGcnJzrelpaGsnJyTRq1IiFCxcyfPhwxo0bl2v44rZt22jQoAH/+te/WLZsGYMGDWLIkCHMnDkzV10bN25k4cKFvPDCCwwcOJCnn36aJUuWFDnWChUq8Morr7B48WJmzZoFwPz581mwYAGvvPIKFStWLMGbKJnMzEz+97//8dRTT7F582buvPPOYtfx/PPPEx0dTbdu3fJcO5kAn/qMFStW5JtvvvEd//73v+fdd99lyZIlHDp0iGXLlrFs2TLuvvvuYscjIiIiIhLKtLppAapUqcJ1113Hgw8+SKtWrahRo0au6w8//DCtW7dmxowZmBkdO3Zk3759pKSkMGDAAABuuukmbrrpJgCcc1xzzTW88847zJo1i549e/rqOnbsGDNnziQmJoaOHTsyY8YMVq5cSdeuXYscb7NmzRg3bhyDBw+mZcuWDBo0iNGjR3PFFVeUwtvwz+HDh31zCCtVqsScOXNo165dseo4ePAgo0ePZt68efler1mzJlFRUbnmjGZkZLBx48ZcvZadO3fmxRdfpHv37pw4cYLw8HAmTZqUb+IpIiIiIlJUDsgOsb630HqaUnbjjTfy8ssv5xlqeuTIEdauXUv37t3JysoiMzOTzMxM4uPj2b59O4cOHQLg6NGjjBgxgkaNGhEREUFERARvvPEG+/bty1VfkyZNiImJ8R03aNAgT5miGDJkCE2bNiU+Pp4LLriAESNG+PHUpadKlSp88sknLF++nJ49e9K3b1/Wrl1brDrGjh3LVVddxfXXX19gmd/97ndMnjyZDz74gJ9++okHH3yQI0eO5Bqaunr1au69915GjRpFWloajz32GA888ADz58/3+/lEREREREKRehIL0bVrVzp27OjrDTxp//79ZGdnM2LEiHwTse+++44mTZrw0EMP8eqrrzJu3Dhat25NZGQko0aN4ocffshV/tQVO8PDw8nIyCh2vGFhYfTs2ZP333+fO+64g/Dw8GLXUZrCwsJ821V06NCBnTt3MnLkSJYtW1ak+3fu3MmkSZNYuXIlhw8fBjzDV7Ozszl8+DCVK1cmPDycsWPHsnHjRq699loA4uLiGDBggG/oLcCwYcO4+eab+etf/wpA27Zt2bVrFyNGjKB79+6l+dgiIiIicpYJ5MqjgaAksRC1a9dm6dKlec5Xr14dM2Ps2LEkJSXlud6wYUPAMy9w2LBh3Hfffb5rZ3IRmT179jBy5Ejatm3LmDFjuOWWW2jQoMEZa6+4WrdunWc+ZmG++eYbjh8/7kv+coqKimL58uV06NCBmjVrsnr1anbu3Mkvv/zCJZdcQq9evWjevLmv/KZNm7j99ttz1dG8eXOmTJlCVlZW0BNqEREREZGyQkmiH8455xzi4+P55ptv8mzsntOxY8eoXLmy73jfvn188MEHXHbZZWckrv79+9O0aVNWrlxJYmIiffv2ZeXKlZgF/jcbzrk87a5fv546derkKXvw4EH27t1LrVq1qFWrlu98q1ateO+993KVHT9+PNu3b+eFF16gRYsWua5dcMEFAOzdu5c333yTf/zjH75rsbGxfP7557nKb9q0iTp16ihBFBERERG/OWcBXXk0EJQk+mn8+PEkJycTGRnpW2Dms88+Y926db55bu3bt2fChAk0bNiQsLAwHn30UWrWrHlG4pkyZQppaWls2LCB8PBwXnrpJVq0aMEzzzzDn//85zPSZmHi4+O59dZbadmyJQCzZ8/mnXfeITU1NU/ZBQsW0K9fP0aPHs2YMWN852NiYvL0Ip533nn88MMPuc6///77bNiwgebNm/P9998zduxYGjZsSL9+/Xxl+vfvz6hRo7jwwgu57rrr+OSTT5g2bRrDhw8v3QcXERERESnnlCT6qX379ixdupQxY8aQmppKxYoVad68OX369PGVmThxIgMHDqRfv37ExMQwdOhQNm3axIYNG0o1lm3btvHggw/y6KOPcumllwKeIa9///vfGTJkCElJSTRp0qRU2zydtm3bMn36dB577DEAmjZtyrx5887I/L+KFSsybdo0vvrqK8455xy6dOnC3//+dyIjI31lhg8fTkREBC+++CJPPPEE9erVY/To0Tz00EOlHo+IiIiInF2yQ2xOooXiZuJxcXFu3bp1wQ5DRERERESKwMzWO+cKnsdVhl3SorJ7dlGjgLSVdOGWgLwn9SSKiIiIiIj4yQFZIbazYGg9jYiIiIiIiJSIehJFRERERET8Fnqrm4bW04iIiIiIiEiJKEk8RWpqKmZG69atc52/9tprMTMmTpwYpMjKn3fffZfWrVtTuXJlWrZsyfLly/2qZ9q0aTRr1owqVapwySWX8PTTTxe7rTlz5pCYmEjt2rWpVq0aCQkJfPjhh37FIyIiIiJykgOyCQvIJ1CUJBZg+/bt7N69G4Aff/yRzZs3Bzmi8mX37t106dKFJk2a8NZbbxEfH0+3bt34+uuvi1XP66+/zsCBA+nevTuLFy+md+/eDB06lGnTphWrrYkTJ1K3bl2ef/555syZQ40aNbjhhhv4/PPPS+2ZRURERERCgZLEAnTs2JE333wTgMWLF5OYmBjkiMqXKVOmULVqVVJTU7nhhhuYOnUqsbGxxe6JnTt3LldffTXjxo2jXbt2jBw5knbt2rFgwYJitbVgwQJSU1O56aab+O1vf8usWbOoVq0akydPLrVnFhEREZGzU5azgHwCRUliAW688UYWLVoEwKJFi7jxxhvzlFmxYgVt2rShcuXKnH/++YwfPz7X9c8//5xbb72VunXrEhkZyaWXXsqkSZNylenbty8JCQlMnTqV+vXrU7NmTYYMGUJ2dnaRY01PT6dKlSq5etfA0xtqZixdurTIdZWWVatWkZiYSMWKFQEIDw8nKSmJVatWFauezMxMoqOjc52Liooi5/6eRWmrZs2aueqoWLEijRo1YseOHcWKR0REREQk1ClJLEBiYiIff/wxP/74I6tWrSI5OTnX9bS0NJKTk2nUqBELFy5k+PDhjBs3Lleitm3bNho0aMC//vUvli1bxqBBgxgyZAgzZ87MVdfGjRtZuHAhL7zwAgMHDuTpp59myZIlRY41Ojqam2++mddeey3X+RkzZlCnTp2g9IJu3bqViy++GPAksQAXXXQRW7duLVY9v//973n33XdZsmQJhw4dYtmyZSxbtoy77767RG0dPXqULVu2cNlllxUrHhERERGRnBxGFmEB+QSKtsAoQJUqVbjuuut48MEHadWqFTVq1Mh1/eGHH6Z169bMmDEDM6Njx47s27ePlJQUBgwYAMBNN93ETTfdBIBzjmuuuYZ33nmHWbNm0bNnT19dx44dY+bMmcTExNCxY0dmzJjBypUr6dq1a5Hj7dOnD0lJSezatYsGDRoAMHPmTO68807Cw8NL+jqK7eDBg0RHRzN79mzuuOMO5syZQ3R0NL/88gsnTpzw9fqdTufOnXnxxRfp3r07J06cIDw8nEmTJhrPYVcAACAASURBVNGtW7cStZWSksKJEycYPHhwqT2ziIiIiEgoUE9iIW688UZefvnlPENNjxw5wtq1a+nevTtZWVlkZmaSmZlJfHw827dv59ChQ4Cnt2rEiBE0atSIiIgIIiIieOONN9i3b1+u+po0aUJMTIzvuEGDBnnKnE6HDh04//zzmT59OgAbNmxg8+bN9O7d259HLzVVq1YlKiqKqlWr+nX/6tWruffeexk1ahRpaWk89thjPPDAA8yfP9/vttLS0nj88cd56qmnaNiwoV9xiYiIiIiclO3CAvIJFCWJhejatSsdO3b09QaetH//frKzsxkxYoQv+YuIiKB79+4AfPfddwA89NBDTJkyhSFDhrB69Wo++eQTkpOTyczMzFVfVFRUruPw8HAyMjKKFWtYWBh33XWXL0mcMWMGLVq04PLLLy9WPaUlJiaG9PR0OnfuTHp6OsnJyaSnp1OpUqUi9yICDBs2jJtvvpm//vWvtG3bluHDh9OvXz9GjBjhV1vbtm2jR48e9OvXj0GDBpXa84qIiIiIhAoNNy1E7dq18130pXr16pgZY8eOJSkpKc/1k71T8+fPZ9iwYdx3332+aydOnDhj8fbp04cnnniCTz/9lFmzZnH//fefsbZO5+KLL2bbtm25zm3fvt03d7CoNm3axO23357rXPPmzZkyZQpZWVmEh4cXua0DBw7QpUsXWrVqlWcBIRERERERfzgI6HzBQAhakmhm9YFXgDpANjDVOffMKWUSgDeAb7yn5jvnxgYyzvycc845xMfH88033xAXF1dguWPHjlG5cmXf8b59+/jggw/O2GIpTZs2JT4+nsGDB7Nnzx569ep1Rtopinbt2vHyyy+TkZFBREQE2dnZLF26lM6dO+cpe/DgQfbu3UutWrWoVatWrmuxsbF59jLctGkTderU8c21LEpbGRkZ3HLLLVSoUIG5c+cSERFxBp5aRERERKT8C2ZPYiYw1Dn3qZlFAevNbLlz7tRd699zznUJQnyFGj9+PMnJyURGRvoWmPnss89Yt26db75c+/btmTBhAg0bNiQsLIxHH300z1YMpa1Pnz4MHjyYxMREYmNjz2hbhbnnnnt49tln6du3L3/4wx+YPn06e/bsyXehmAULFtCvXz9Gjx7NmDFjcl3r378/o0aN4sILL+S6667jk08+Ydq0aQwfPrxYbd1zzz188MEHvPrqq2zZssV3Pjo6mmbNmpX+CxARERGRs4IjsHsYFsbMqgEvAM3xdHL+HvgSmA00BHYAtznn9hdWT9D6RZ1ze51zn3q/PgRsAeoGK57iat++PUuXLmXjxo3ccsst9OzZk7feeivXVhkTJ06kdevW9OvXj/vvv5/evXvnOzy1NHXp4smng9mLCFC/fn2WLFnC5s2bSU5O5uOPP2bhwoU0bty4WPUMHz6cxx9/nJkzZ9KlSxdefPFFRo8ezahRo4rV1ooVKzh+/Di33XYbv/nNb3wfzUsUERERkRDyDLDUOdcEuBxPjjUCWOmcawys9B4XynJuSh4sZtYQeBdo7pxLz3E+AZgHfAvsAYY55zadrr64uDi3bt26MxJrWTd16lQeeOABvv/++zwL4oiIiIiIlEVmtt45V/A8rjKsUYuqbsz8lgFpq+8lawp8T2YWDfwHuNDlSPLM7EsgwTm318xigTTn3KWFtRP0hWvMrCqeRPDPORNEr0+BC5xzh82sE7AQyLcryswGAgMB3z6BZ5MdO3bwxRdfMHbsWO68804liCIiIiIiZ5cLgR+B/zOzy4H1wJ+A2s65veAZzWlm552uoqAuw2NmEXgSxOnOuTwb3znn0p1zh71fvwVEmFmtU8t5r091zsU55+LOPffcMxp3WTRmzBi6du1K06ZNeeKJJ4IdjoiIiIiIlL5aZrYux2dgjmsVgNbAFOdcK+AIRRhamp9grm5qwIvAFufcPwsoUwfY55xzZnYVnqT25wCGWW6kpqaSmpoa7DBERERERM4qzkFW4Da6/6mQYbnfAt865z72Hs/FkyTuM7PYHMNNfzhdI8EcbnoN8DvgczPb4D33F6ABgHPuOaAHcI+ZZQLHgDtcWZhEKSIiIiIiUoY45743s91mdqlz7kugPbDZ++kDpHj/fON0dQVzddP3nXPmnGvpnLvC+3nLOfecN0HEOTfROXeZc+5y51wb59yHgYovNTWVyy+/nCpVqlCvXj26d+/O+vXr85SbMGECaWlpRa43ISGBvn37ll6gZdi7775L69atqVy5Mi1btmT58uUlqm/Lli1ERETkuzfl6dqaM2cOiYmJ1K5dm2rVqpGQkMCHHwbsr5OIiIiIhCwjO0CfIrgPmG5m/wWuAMbjSQ4TzexrINF7XKigzkksq6ZOncqAAQPo3r07ixcvZvz48WRlZfHRRx/lKVvcJPFssXv3brp06UKTJk146623iI+Pp1u3bnz99dd+1zl06NB8F+QpSlsTJ06kbt26PP/888yZM4caNWpwww038Pnnn/sdj4iIiIhIWeKc2+Bdp6Wlc+4m59x+59zPzrn2zrnG3j//d7p6gr66aVk0YcIE7r77bkaPHu0717t3b06cOBHEqMqXKVOmULVqVVJTU6lYsSLXX389aWlpTJw4kWeeeabY9b399tvs2rWLLl26sHnz5mK3tWDBAmrWrOm7JyEhgfr16zN58mSmTJlSsocVERERkbOWI6BzEgMitJ6mlOzcuZPatWvnOV+xYkXf12aGmbFz504eeeQR33HOoaTOOUaOHEmtWrWoXr16rg3giyozM5O4uDh69OiR6/ztt99Oy5Yty2ziumrVKhITE33vLDw8nKSkJFatWlXsurKyshg6dCiPPfYYYWF5/8oWpa2cCSJ4vpeNGjVix44dxY5HRERERCSUKUnMR7NmzZg8eTKLFi0iIyMj3zJr1qxhzZo11KlTh/79+/uO//a3v/nKTJs2jZSUFIYNG8Zrr73GypUrWbt2bbFiqVChAq+88gqLFy9m1qxZAMyfP58FCxbwyiuv5Epcy5KtW7dy8cUXA5Ce7tn+8qKLLmLr1q3Fruv5558nOjqabt26lVpbR48eZcuWLVx22WXFjkdEREREJKcswgLyCRQlifl49tlnycrKolu3blSvXp0ePXqwevXqXGXatGlDmzZtiIyMpF69er7jiy66yFdmwoQJ9OrVixEjRtC5c2dmz57N8ePHix1Ps2bNGDduHIMHD2bz5s0MGjSI0aNHc8UVV5T4Wc+UgwcPEh0dzezZs4mJieH1118nOjqaX375pVi9nwcPHmT06NGkpBQ8v9aftlJSUjhx4gSDBw8u9rOJiIiIiIQyJYn5aNOmDV999RUvvPACHTt25O2336Zdu3a8+uqrRa4jIyODL7/8krZt2/rO1atXj0suucSvmIYMGULTpk2Jj4/nggsuYMQIv/bFDLiqVasSFRVF1apV/bp/7NixXHXVVVx//fWl1lZaWhqPP/44Tz31FA0bNvQrLhERERERAIeR7QLzCRQliQWIiYmhf//+zJs3j927d3PllVfmWsjmdH7++Weys7OpXr16rvOnHhdVWFgYPXv25OjRo9xxxx2Eh4f7VU+gxMTEkJ6eTufOnUlPTyc5OZn09HQqVapU5CGyO3fuZNKkSfzlL3/h8OHDHD58mMzMTLKzszl8+DBZWVnFbmvbtm306NGDfv36MWjQoFJ/bhERERGR8k5JYhFUq1aN3r17s2PHDl9icjo1a9YkLCyM/fv35zp/6nFR7dmzh5EjR9K2bVvGjBnDrl27/KonUC6++GK2bduW69z27dt9cweL4ptvvuH48eNce+21REVFERUVxfTp0/nss8+IiorinXfeKVZbBw4coEuXLrRq1YpJkyb5+WQiIiIiIrlpTuJZ4Icffshzbvv27cTGxubpwYuKiuLIkSN5ykdERNC0adNccxm/++47vvrqK79i6t+/P02bNmXlypVceeWV9O3bF+ecX3UFQrt27Vi+fLlv4Z/s7GyWLl1Ku3bt8pQ9ePAgX3zxBT/99FOu861ateK9997L9UlOTubSSy/lvffeIz4+vshtZWRkcMstt1ChQgXmzp1LRETEmXp0EREREZFyTfsk5qN9+/Zcc801dOrUiapVq/L+++8zceJE/vKXv+Qp26JFC+bOnUunTp2IjY0lJiaG2NhYAO677z7uvfdemjZtSsuWLRk/fjyRkZHFjmfKlCmkpaWxYcMGwsPDeemll2jRogXPPPMMf/7zn0v8vGfCPffcw7PPPkvfvn35wx/+wPTp09mzZ0++C8UsWLCAfv36MXr0aMaMGeM7HxMTw7XXXpur7HnnnccPP/yQ63xR2rrnnnv44IMPePXVV9myZYvvfHR0NM2aNSvFJxcRERGRs4kDskNsn0Qlifm47777mD59OvPmzePo0aM0atSIlJQU7r///jxlx40bx4ABA+jWrRuHDh2iT58+pKamAjBw4EB27drFk08+SWZmJvfffz8VKhTvlW/bto0HH3yQRx99lEsvvRSAhg0b8ve//50hQ4aQlJREkyZNSvzMpa1+/fosWbKEP/3pTyQnJ9O4cWMWLlxI48aNg9LWihUrOH78OLfddluue9u2bUtaWlqpxyQiIiIiUl5ZWR6y6K+4uDi3bt26YIchIiIiIiJFYGbrnXNxwY7DH/Wbx7g/vd4mIG092GxZQN5TaPWLioiIiIiISIlouKmIiIiIiIifQnFOYmg9jYiIiIiIiJSIehJFRERERERKIAsLdgilSj2JIiIiIiIi4qMksQCpqalcfvnlVKlShXr16tG9e3fWr1+fp9yECROKtYVCQkICffv2Lb1Ay7B3332X1q1bU7lyZVq2bMny5cuLXcecOXNITEykdu3aVKtWjYSEBD788EO/2iqNeEREREREcnLOyHZhAfkEipLEfEydOpUBAwbQvXt3Fi9ezPjx48nKyuKjjz7KU7a4SeLZYvfu3XTp0oUmTZrw1ltvER8fT7du3fj666+LVc/EiROpW7cuzz//PHPmzKFGjRrccMMNfP7558Vqq7TiEREREREJdZqTmI8JEyZw9913M3r0aN+53r17c+LEiSBGVb5MmTKFqlWrkpqaSsWKFbn++utJS0tj4sSJPPPMM0WuZ8GCBdSsWdN3nJCQQP369Zk8eTJTpkwpclulFY+IiIiIyKmytLpp6Nu5cye1a9fOc75ixYq+r80MM2Pnzp088sgjvuOcQ0mdc4wcOZJatWpRvXp1Ro0aVexY0tPTqVKlCtOmTct1fvv27ZgZS5cuLXadgbBq1SoSExN97yw8PJykpCRWrVpVrHpyJojg+R40atSIHTt2FKut0opHRERERCTUKUnMR7NmzZg8eTKLFi0iIyMj3zJr1qxhzZo11KlTh/79+/uO//a3v/nKTJs2jZSUFIYNG8Zrr73GypUrWbt2bbFiiY6O5uabb+a1117LdX7GjBnUqVOHxMTE4j9gAGzdupWLL74Y8CS6ABdddBFbt24tUb1Hjx5ly5YtXHbZZcVq60zFIyIiIiJnNwdkYwH5BIqSxHw8++yzZGVl0a1bN6pXr06PHj1YvXp1rjJt2rShTZs2REZGUq9ePd/xRRdd5CszYcIEevXqxYgRI+jcuTOzZ8/m+PHjxY6nT58+vPfee+zatct3bubMmdx5552Eh4f7/6Bn0MGDB4mOjmb27NnExMTw+uuvEx0dzS+//FKiYbspKSmcOHGCwYMHF6utMxWPiIiIiEioUZKYjzZt2vDVV1/xwgsv0LFjR95++23atWvHq6++WuQ6MjIy+PLLL2nbtq3vXL169bjkkkuKHU+HDh04//zzmT59OgAbNmxg8+bN9O7du9h1BVrVqlWJioqiatWqJa4rLS2Nxx9/nKeeeoqGDRv61VZpxiMiIiIiAkaWCwvIJ1CUJBYgJiaG/v37M2/ePHbv3s2VV16ZayGb0/n555/Jzs6mevXquc6felwUYWFh3HXXXb4kccaMGbRo0YLLL7+82HUFSkxMDOnp6XTu3Jn09HSSk5NJT0+nUqVKueZ2FtW2bdvo0aMH/fr1Y9CgQcVuq7TjEREREREJVUoSi6BatWr07t2bHTt2kJWVVaR7atasSVhYGPv37891/tTjourTpw+bNm3i008/ZdasWWW+F/Hiiy9m27Ztuc5t377dNy+wOA4cOECXLl1o1aoVkyZN8qut0oxHREREROQkB2Q7C8gnUJQk5uOHH37Ic2779u3ExsbmmQMYFRXFkSNH8pSPiIigadOmueYyfvfdd3z11Vd+xdS0aVPi4+MZPHgwe/bsoVevXn7VEyjt2rVj+fLlvoV/srOzWbp0Ke3atctT9uDBg3zxxRf89NNPea5lZGRwyy23UKFCBebOnUtERIRfbRUnHhERERGRs5mSxHy0b9+eu+++m0WLFrFq1SrGjh3LxIkTGTBgQJ6yLVq0YO7cubzzzjt88cUX7N2713ftvvvuY/r06aSkpPDWW29x++23ExkZ6Xdcffr0Yc2aNbRr147Y2Fi/6wmEe+65h/T0dPr27cs777zDwIED2bNnT64FZ05asGABTZs2ZeLEifnW88EHHzBq1Ci2bNnCRx99xEcffcTmzZuL1VZx4hERERERKY4swgLyCZQKAWupHDmZ3M2bN4+jR4/SqFEjUlJSuP/++/OUHTduHAMGDKBbt24cOnSIPn36kJqaCsDAgQPZtWsXTz75JJmZmdx///1UqOD/K+/SpQuDBw8u872IAPXr12fJkiX86U9/Ijk5mcaNG7Nw4UIaN25crHpWrFjB8ePHue2223Kdb9u2LWlpaUVuq7TiEREREREJdeacC3YMpS4uLs6tW7cu2GGUuqlTp/LAAw/w/fffExUVFexwRERERERKhZmtd87FBTsOf9S5rIbrPaN9QNr6xxVzA/Ke1JNYDuzYsYMvvviCsWPHcueddypBFBERERGRM0ZzEsuBMWPG0LVrV5o2bcoTTzwR7HBERERERCSHbMIC8gkUJYnlQGpqKhkZGSxfvpwaNWoEOxwREREREQlhGm4qIiIiIiLiJ+cgK4B7GAaCehJFRERERETER0miiIiIiIiI+Gi4qYiIiIiISAlka7ipiIiIiIiIhCr1JIqIiIiIiPjJYWS70Op7C62nERERERERkRJRT6KIiIiIiEgJZKE5iSIiIiIiIhKi1JMoIiIiIiLiJ4dWNxUREREREZEQpp5EERERERERv2l1UxEREREREQlh6kkUEREREREpgWytbioiIiIiIiKhSj2JIiIiIiIifnIOsrS6qYiIiIiIiIQq9SSKiIiIiIiUgFY3FRERERERkZClnkQRERERERE/OYxszUkUERERERGRUBW0JNHM6pvZO2a2xcw2mdmf8iljZvYvM9tqZv81s9bBiFVERERERKQg2VhAPoESzOGmmcBQ59ynZhYFrDez5c65zTnKJAONvZ+rgSneP0VEREREROQMCFqS6JzbC+z1fn3IzLYAdYGcSWI34BXnnAM+MrNqZhbrvVdERERERCSoHGhO4plgZg2BVsDHp1yqC+zOcfyt91x+dQw0s3Vmtu7HH388E2GKiIiIiIiEvKCvbmpmVYF5wJ+dc+mnXs7nFpdfPc65qcBUgLi4uHzLiIiIiIiIlDbtk1iKzCwCT4I43Tk3P58i3wL1cxzXA/YEIjYREREREZGzUdB6Es3MgBeBLc65fxZQbBEw2Mxm4Vmw5qDmI4qIiIiISJnhQm+fxGAON70G+B3wuZlt8J77C9AAwDn3HPAW0AnYChwF+gUhThERERERkbNGMFc3fZ/85xzmLOOAewMTkYiIiIiISPE4COgehoEQWjMsRUREREREpESCvrqpiIiIiIhIeRZqcxLVkygiIiIiIiI+6kkUERERERHxk0M9iSIiIiIiIhLC1JMoIiIiIiJSAqHWk6gkUUREREREJESY2Q7gEJAFZDrn4sysBjAbaAjsAG5zzu0vqA4NNxUREREREfGTw8h2gfkUww3OuSucc3He4xHASudcY2Cl97hAShJFRERERERCWzfgZe/XLwM3FVZYw01FRERERERKIJuAzUmsZWbrchxPdc5NPaWMA5aZmQOe916v7ZzbC+Cc22tm5xXWiJJEERERERGR8uGnHENIC3KNc26PNxFcbmZfFLcRDTcVEREREREJEc65Pd4/fwAWAFcB+8wsFsD75w+F1aEkUURERERExF+OMrNwjZmdY2ZRJ78GfgtsBBYBfbzF+gBvFFaPhpuKiIiIiIiEhtrAAjMDT643wzm31Mw+AeaYWX9gF3BrYZUoSRQREREREfGTg+JuT3HGOOe2A5fnc/5noH1R69FwUxEREREREfFRT6KIiIiIiEgJlJWexNKinkQRERERERHxUU+iiIiIiIiInxxFW3m0PFFPooiIiIiIiPioJ1FERERERKQEnHoSRUREREREJFSpJ1FERERERKQEslFPooiIiIiIiIQo9SSKiIiIiIj4yTntkygiIiIiIiIhTD2JIiIiIiIiJaDVTUVERERERCRkqSdRRERERETEb6Y5iSIiIiIiIhK61JMoIiIiIiJSApqTKCIiIiIiIiFLPYkiIiIiIiJ+cmifRBEREREREQlh6kkUERERERHxlwPngh1E6VJPooiIiIiIiPioJ1FERERERKQEstGcRBEREREREQlR6kkUERERERHxk0P7JIqIiIiIiEgIU0+iiIiIiIiI30z7JIqIiIiIiEjoUk+iiIiIiIhICWifRBEREREREQlZ6kkUEREREREpAa1uKiIiIiIiIiFLPYkiIiIiIiJ+ck49iSIiIiIiIhLC1JMoIiIiIiJSAtonUUREREREREKWkkQRERERERHx0XBTERERERGREnAu2BGULvUkioiIiIiIiE9Qk0Qze8nMfjCzjQVcTzCzg2a2wfsZFegYRURERERECuOcBeQTKMEebpoKTAReKaTMe865LoEJR0RERERE5OwW1CTROfeumTUMZgwiIiIiIiL+cgS2ly8QysOcxN+Y2X/M7N9mdlmwgxEREREREQllwR5uejqfAhc45w6bWSdgIdA4v4JmNhAYCNCgQYPARSgiIiJlSmLYrcW+Z3n262cgksA63XOHwjOKlFUhtrhp2e5JdM6lO+cOe79+C4gws1oFlJ3qnItzzsWde+65AY1TREREREQkVJTpnkQzqwPsc845M7sKT1L7c5DDEhERERER8XCE3JzEoCaJZjYTSABqmdm3wGggAsA59xzQA7jHzDKBY8AdzoXaVpUiIiIiIiJlR7BXN+15musT8WyRISIiIiIiUjaFWDdWmZ6TKCIiIiIiIoFVpuckioiISOnyZ+VPKHsrY/r7HMFsv7TeYbCf/UzTKq1SHmlOooiIiEgZd7YmEmfrc4tI6VKSKCIiIiIiUgKhtrSm5iSKiIiIiIiIj3oSRURERERE/OQIvTmJ6kkUERERERERH/UkioiIlFOFrQJpFSLOePsdq/bJ93z20aP5nj9dTMtOzMhzzt+VPH9b8c4851xmRrHiKqh8aQrUSp6FtaPFbkRKyAEh1pOoJFFEROQsYhUi8k3GyptAJMFlVagndaH+fCLlgZJEERERERGREtDqpiIiIiIiIhKy1JMoIiIiIiJSEupJFBERERERkVClnkQRESkXArUKZFnk7wqfpVVfWJUqpdp+QQpaLTVU+PN9LO49/vw7KKiN4talFVTl7GUht0+ikkQREZEQ5DIz8v3BvDQTzrBzipc8FrQ1BhSciBZ0vrC6ClLcrS782bKjIKWd6BdXQUlasOMSkbJJSaKIiIiIiEhJaE6iiIiIiIiIhCr1JIqIiIiIiPjLEXJzEtWTKCIiIiIiIj7qSRQRERERESmJEJuTqCRRRCSItGT8mVXaKzcW53vy24p3Frv+glbZ9Jc/MUjJBPudF9Z+cVZjFZHyy8zCgXXAd865LmZWA5gNNAR2ALc55/YXVoeSRBERKRfO5qQ5vFq1fM9nHThQ4D0Fbd9Q0Pnw6jEF1pV9pHjbTVjFivlfKGwLjAK208jafzD/8n7s3VjQthmn2+qiNFiFiHyTtKQaAwq8p7Dvb37/HvxJUEvz39XZ/G9UBMrUnMQ/AVuAaO/xCGClcy7FzEZ4j4cXVoHmJIqIiIiIiIQAM6sHdAZeyHG6G/Cy9+uXgZtOV496EkVEREREREqi7MxJnAA8BETlOFfbObcXwDm318zOO10l6kkUEREREREpH2qZ2bocn4EnL5hZF+AH59z6kjainkQREREREZHy4SfnXFwB164BbjSzTkAlINrMXgP2mVmstxcxFvjhdI0oSRQRkaAo7squp1uUIxArN+YXcyAWPZGyqbBFZ0qjfGH8WaQmv7+/BS2KBLD0f9OK3YbIWasMDDd1zj0MPAxgZgnAMOfcXWb2D6APkOL9843T1aUkUUQkiLQa4JlV0Pv1J+Es7nYaVrH4yWNBW2BkHz6S7/nCfsAvaGXMCufH5t/GgfxXEQUIa1A3/wv/v717D7Lsqu9D//3NS5rRaPQaIQn0ANsCQ/wQWNHlUUlhi4eEcXBijDFlDC7HcjAkJhcqF4gd1712GZKyU+SWMVi2MfIDg0qEa4wFEighBIONJCIeQhYWBNCgt9BI85JmenrdP6ZpjzRr90y/zuk+/flUnZo+v73PXuucfXr6rPPb67fun7OC+jG3nSTZfHw3vH7L5m68Pbhr8FBt//5ufLDi68C5Gqq4miTtwFQ/PnAOl7JC7VyGvixZ6uVggFXlbUmurKqfT/LNJEf9D8EgEQAAYKFakrailsBIa+0TST4x8/P9SS6ez+MVrgEAAGCWTCIAAMAitBUwJ3EpySQCAAAwSyYRgIk1VKzjaAVBFPmA4d+DcRbcmm9VZBiZCcskGiQCMBbj/EA3VIFyIdZtPaEbr4GqnEnS9u7rxjecfNK89q+Nw3/GN5x3zuC2MYch/QAAIABJREFUngPPeurgto0PPNyNr8sp/WOdcWI3vunu4YqkGaoWum3g9d0wx0eYgaqrQ1Vap/7u7/v7D5yPJIP9navibLeNLXNUUN0//D6d75IvQ1+MTO/esyTVfIHJYpAIAACwGCusuulimZMIAADALJlEAACARagJm5MokwgAAMAsmUQAAICFalHdFACWwnyrJx5t2QpYS5by90clU+CxDBIBWBXa1IF5L5sx9MF4riUwhh6z/pQ5lkOYr3PO7Mfvf7AbPvCD39WNrzswPdjEvscd140f90D/uU9tXj94rI39FSWyf2CpiyF7vqe/ZEaSbL5jTze++0n9Nrbc2V+WI5njw82+gaU8nvGPuvF24OBgG3Wg/zoefKB/DofeP3MtY9I29ZfZOLhz5+BjBo818Puz2tYdXIl9gqRUNwUAAGByySQCAAAsxoTNSZRJBAAAYJZMIgAAwGJMWCbRIBGARTtadcSlKjYx1M6G07cvyfFhks23iunQ/orHwOQzSARgRVm3ZUs3Pr1377yPNVTFdKiNJKlNm/obtmzuhg+ecXI3vv+U4wfb2Hzbfd341BPP6MYf3t7v054zh2eNHBx4GqdM9b/uvv+pc3wkeGq/Mue22/vVVTfu7sf3nTZcQfXA1n4V021f2dWNTx8/vKTDru9/XDd+wjd3Dz6mpx7qV1xNkpzWr9S6fmO/X9M7+1VP59L299+/cy1nce3+93bjlrmAZTZhmURzEgEAAJglkwgAALBQLdZJBAAAYHLJJAIAACxCTdicRINEAI7ZQotfKJoBo+F3DVgKBokALNpc1RaHKowm/VL6L9j0iu6+G847Z979GnrMI9/Tr36ZJJvu6Ve03Pf4fvXNTQ883I3vOnv4NbnzWWd149u/OFAVdHt/dkgbLhaak/73wW78vu/r/+mv/u5Jkunj+vG7L+rPwTnh9n4bBweOkyQn7ujH95y7tRufq1Lq8Tv7r+Mjp/er2h5370Dl3M3DFWqzr3/eBx8zUN10rqq9c1XhvWb3FUfEFjpA7P0eGmzCPE1YJtGcRAAAAGaNdZBYVe+uqnuq6ksD26uq/t+quq2qvlBVzxh1HwEAANaScWcS35Pkkjm2X5rk/JnbZUneOYI+AQAArFlHHSRW1euq6pTlaLy19skk355jl5ck+eN2yN8kObmq+hM5AAAAxqDaaG6jciyFa85Mcn1VfS7Ju5Nc01obVRefkOT2w+7vmIndOaL2AVacuQpK9ApQLORYcxWiAcZvNRWWWcr/s4DROOogsbX2K1X1q0lekOTnkvxOVV2Z5A9ba19d5v71yqZ1B6hVdVkOXZKac889dzn7BLBmDVVbbPuHK5jOu42tJ/Tb2NaPJ8n05v6gtqb6VS7nsve8bd34Q+f1/2Sum+q3/fAc1+Bs6he6zD0/1L/AZ+qE/vOY3jJckvSh7+5X/2zr5yhjOuCEb/SPtfnOfnXTqYFT9fDpw+ejbeg/93WP9Nveeufw99VTm/v92nJn/3267va7+306MDXYRp3WP8FTX/3f3fj6k0/uH2eO3525Kp8uxHwGZAZvME+t///OanVMcxJnMod3zdymkpyS5Kqq+k/L2LfkUObw8PrlZye5Y6CPl7fWLmytXXj66acvc7cAAAAm07HMSfw3VXVjkv+U5K+TfH9r7TVJfijJTyxz/z6U5Gdnqpw+M8mDrTWXmgIAACtDG+FtRI5lTuL2JP+itfaNw4OttemqevFiGq+qP0/y3CTbq2pHkl9LsnHm+O9KcnWSFyW5LcneHLrcFQAAgGVyLHMS/8Mc225ZTOOttZ8+yvaW5LWLaQMAAGBZjTDLNwrjXicRAACAFeRYLjcFYJUYKjW/UisVXvqk/3PcXYA1Y7X9/wCMj0EiwCoz9IFuIeumDa2HOLQMRZJ89Nu/f0Ts0rP/zeD+be++ftsnndh/wDlndsMHTzx+sI31ux7uxvc/rv88amr4uqA9Z/T/NO46b37XEq1/ZLgc+iOn9o/15Iu+0Y3fdne/avfUXZsH23jt867txt9x43O78S1fHH599549vyU4Tvhq/311ypeHL2Da9SN7uvFT/7K/7MqGfcPnY8ud/ffDvU/vH+usu/vvxXpo12AbQ4aWupje3X9+tWl4TdI2Nf+lZZby/4elYhDKWjDKhe5HweWmAAAAzJJJBAAAWAyZRAAAACaVTCIAAMBiTFgm0SARYATmKhox36IO4yxAkSQv3PqqI2LrTj5pDD0BAJaDQSIAR2j798/vAdu2Dm87bWAAeaBfGXOoIun6ff39k2TXk/sVJfecOf9ZFY+c0o+f+NQHuvF9N57ajf/wj35usI0HDwxXJe35tz9wXTf+racMdDbJld/4oW785y74TDd+x1OHB/p//a3v6sZ37ehXBR2q3loHhyu+nvjf+ud98739Cp8Htq4fPNaBE/sVQ8+67t5ufPqkftXTdXNUN21bjuvH7+z/7gxVDD64c+dgG8nSVQZVYRSWTzXVTQEAAJhgMokAAACL0YavlFiNZBIBAACYJZMIAACwGBM2J9EgEWANGKqIWhv6BT6AlWsh1ZJfsOkV3Xib6hcGmutYwOQzSAQYgaX8sLWUS2bUpv4gcXrv3m583ZZ+Fcjpb35ruAPf+6T+Yzb32x6qYnrv0/ttJ8n+fpHN7D27f6xNDwxXxmxP3tONf+/J3+7G/9f2flXQHz3l84NtXHXfP+7GX7r9+m78b3Z/Tzd+/4F+xcwkeeqpd3fjZ2/qP4/tG4cree47c1M3/j8eeHI3vn7g9Z0a7m7awCnZuKf/UWXb1/YNHmv9roeHG+q1vWH+s2/qrvv6Gzb1X6vp3f33lS9qYDKobgoAAMDEkkkEAABYDJlEAAAAJpVMIgAAwEK1yZuTaJAIrBkLqQg4X0MVBJPk2v3vXZI2hp7H+pNPXpLjAyvHXP9vASwXg0SACbfhvHMGt03fe3//Md/dr0iahwYqYG4bKC+a5JFTju/Gpzb3y1k+cnI/ftyDw1/TPnJKdeMnfKN/rKFKmkkytCDAto39ipn//Nk3dONDFUyT5I/O/Z/d+BvuekY3/sHP/lA3/rp/8vHBNn7vS/+kG9+xp1+N9Zv3nTp4rDNOeagbr4f6HyMOHtc/V9u/MNhEjtvZr0Q7ZN3Dw0s3tI39EzxYUXdHv1JpOzA12MZQtdJ1W/slXBey1MRSfrG1VF9SAQMmLJNoTiIAAACzZBIBAAAWQyYRAACASSWTCAAAsAiTVt1UJhEAAIBZMokAAxZSer429KsnJv3lMWrT8P7X7L5i3u0Dq884l7kYxbI9wOpjkAisGUu1FuJc2tSBbjtzfRCbr3VbtvTb3r+//4Bt/ZL8SbJuqMT/vv5yDzntlP7u5w2v0bj+4f7SBnvO6P8JOn7ndL+N7cMXv2zor0aQ/f3VHvLImcNLG5w9sNzDQwf6S3l87q6zu/EDU8PrbPzq8Q9041/fc1o3fuKZu7vxL+7qt50k+x/a1I1/M/2lLg48cNzgsXY8cHo3vumh/jl5/Kf6yz1semDgfZXk4MDyFJt29F+rTA2fw2zun6uh9/XBe/pLYMz1Jc7QtrZ/f/cLHusdAquJQSIAAMBimJMIAADApJJJBAAAWKimuikAAAATTCYRYIV64dZXjbsLAMCxmLBMokEiwIChaqgLqVI431LyCxkgrjuhX/X04C1fHXxM/cCT+8e6/e5ufP8ZJ3bjx9+xa7CN+5/Rr4g6ZKiK6cNzHGbfWf2KqOtOf6QbP/nEfYPHuvuBbcMNdaxf32/7pC3DbewYeDJfuP0J3fjBR/qVUv926rzBNjac0K/+OXXX5m58y53DFxdt3NuPb7mn/9yHKtrOZeO3+40cPG1rv40d/YqkSdIevLcbn949UAZ3geazTM0oqisvhGUuYHJU1fFJPpnkuBwa513VWvu1qjo1yfuTPDHJ15O8rLU2UDr6EJebAgAALEYb0W1ujyT5kdbaDya5IMklVfXMJG9Kcl1r7fwk183cn5NBIgAAwCrXDvnOwrobZ24tyUuSfOfShyuS/PjRjmWQCAAAsECVQ9VNR3E7al+q1lfVTUnuSfKx1trfJjmjtXZnksz8+7ijHccgEQAAYHXYXlU3HHa77PCNrbWDrbULkpyd5KKq+r6FNKJwDcASW0hhm8dat6VfhAYAWIFGV930vtbahUfbqbW2s6o+keSSJHdX1VmttTur6qwcyjLOySARYAWa3jtQTjLJ+pNP7sYPPvBgN75u6wmDx6oH++20M7d345vu7lcxbRv71TeTZNvX+xVG9561qRufOq668QPb+vEkmd7Sr6a59fp+Jc89z+xX5UySi7/r1m78ozf+wOBjejadO1w47lOf/Efd+NYd/ec4PfDyHjyu/xomyZb+yz7o9C8MP2B6Q//Co3VT/ddx09f71UXbg8NVcGtj/yPJ+o0b+33a2X+/z6VNHRjc1qs+ahkaYDWpqtOTHJgZIG5O8rwk/zHJh5K8KsnbZv79i6MdyyARAABg9TsryRVVtT6HphVe2Vr7cFV9JsmVVfXzSb6Z5KiXPBkkAgAALNQxFpVZ9m609oUkT+/E709y8XyOpXANAAAAs2QSAQAAFmMFZBKXkkwiAAAAs2QSgTVvaMmK2tCvqggwX3MtjdOrrAqsMhOWSTRIBBiBuQac1+5/7xGxS079hcH9p3fvmVfbc+2//qQTu/F6aKobb9uGl9MYsmFXf2mFzRv6yz3s/J7juvGzPt3vU5IcPK6/RsS3n9rff90tw8/j2m8dMec/SbLx4X5/Dx7f/2Rw39/1lxFJksd/ur90RFvfb2PPmf0LfzY9NNhENu3ptzE98LoPLXORJJu/uKO/4cDAkhLbBt5XW/pLkhxq5Phu+OA3BtoegbZ/eMkMgElmkAgAALAIK6G66VIyJxEAAIBZMokAAACLIZMIAADApJJJBADWvF710VFVOFb5FFa5lonLJBokAgxoUweW7APaCza9Yn5t798/vG2qX3Fxw+n9aprTe/YONzRQUTL7Hh5+TEft7VcwTZLd3/+4bvzb39uvSLr9S/0qpge2Dl/8su+0/rbTvnSw36cn9NtOkjNu6D9m32n9x2x4pP/JYKiK6FztH/dg/1gn3zb/KpvH37GrG6+HBqrdDr0X5tAODFScvf+BbnjOSruP679/h97vcxka3NWGjd1qwkODtKX8P2DoOHMNEAHGxSARAABgEVQ3XUJVdUlV3VpVt1XVmzrbn1tVD1bVTTO3/zCOfgIAAKwVY8skVtX6JO9I8vwkO5JcX1Ufaq19+TG7/s/W2otH3kEAAIBjIZO4ZC5Kcltr7Wuttf1J3pfkJWPsDwAAwJo3zjmJT0hy+2H3dyT5Pzr7PauqPp/kjiRvbK3dPIrOAaNztMINq6m639BzmatK4iWn/sJydQcAGIFJm5M4zkFir+zbY1/ezyU5r7W2u6pelOT/S3J+92BVlyW5LEnOPffcpewnMOFGMQgdqpI4NECc3jtckXTdli3d+NS9981r/ySDVSizZXM3XAf6lT8feeKpg02c+Nlv9uN/3a9a2c7sV7kcajtJtpxxYje+fl+/jU27Ng0ea8j2q7/WjT/yA+d148d94RvDB9s48KXBtq39+MB5GnqtkiS33zXQdv9Pf9t2wuCh2s4Hu/Gh9+lClo6YHmhj3ZYtuWb3FfM+3nyM84uo1fQlGLB2jPNy0x1Jzjns/tk5lC2c1Vp7qLW2e+bnq5NsrKruX8TW2uWttQtbaxeefvrpy9VnAACAR2sjuo3IOAeJ1yc5v6qeVFWbkrw8yYcO36Gqzqyqmvn5ohzq7/0j7ykAAMAaMbbLTVtrU1X1uiTXJFmf5N2ttZur6l/NbH9XkpcmeU1VTSXZl+TlrbUJu+IXAABYtUac5RuFcc5J/M4lpFc/Jvauw37+nSS/M+p+AQAArFVjHSQCLLWVXCn1BZtecURs3dbhYiEAwMpX6VfkXM0MEoGxWwvV/YaqPU7v3pNr97/3iPhcg922v1+xc8Pp/UqX03uGK6UOOjDVj/eLnmbTPXuGjzVUsXNA27i+G6+7+tVbk2SwVulDu/ttPPGMwWPVZ/srLU1v6p/D4771UP9AQxVMM1zJs/bu68cHqs1O3/LVwTaGvoAYej+0OY6VpPs+HfLCra/qtz1H1d421X9fAzB6BokAAACLMWFzEsdZ3RQAAIAVRiYRAABgEUomEQAAgEklkwgAALAYE5ZJNEgEGLPe0hiw0vQq7g5V7a2BSrAArA4GicBEGfdyGkPtDy1pMfgheyCeDC8VsJClLqbu7S8rseHxZ/Xbvuvebvzgzp2DbWz43vP7j7nt69344NqRp50y2MZQv4aWjlh301cGj5WBAc66k0/qxoeex1zmPYg60D/ncy0bcXDnzu77cSFfSsx3eYq2/8C8l3ZJxv/7C7BgE5ZJNCcRAACAWTKJAAAAC9VUNwUAAGCCySQCAAAsxoRlEg0SAQYcrcjGEMU3WCkuPeM14+4CAKuQQSLACAxVK52rauR8B5tDVSvnamPdli3d+MF7+lVPh2w475zBbYOVRwcqfLb9+/sHuv+BwTbqpBO78YPfumvwMfM1dced3fj6k09e0PE++u3fP+Z9h87t0PlLknUn9LetP+WkfOTudx5z28nwFybzrXoKMKnMSQQAAGBiGSQCAAAwy+WmAAAAi+FyUwAAACaVTCLAKjNURGSoOA4r01Axmmv3v3fEPTlkodV8AZi8wjUGiQADlnIpi6FKnrVpY67ZfcWStNGmDiz78htDA5upb9w++Jihweu6rSd04wd37uwfaO/e4Y4NPGao8uhgGxmuGDpUyXPoWEs5aB9qe65zvtSDPku7AKwdBokAAAAL1WJOIgAAAJNLJhEAAGAxZBIBAACYVDKJwEQ5WrGOlVh8Y6jPS9XXhbwmL9z6qiVpm/lTZRRgdamobgrAAgxVMF3IgGCcA92hKpsLeczBnTu7z2Wu12TouV9y6i8MPuaj3/79o/Rwce0P7X+012qpnvtSGcX7aiV+SQPAkQwSAQAAFmPCMonmJAIAADBLJhEAAGARqk1WKlEmEQAAgFkyiQAAAAvVMnFzEg0SAVaooUqXtWHj4GOu3f/e5erOsrDcAwCsPNUm7PrZJLnwwgvbDTfcMO5uACzYXIOncQ4SF7JEw2pb5gOA0auqG1trF467HwtxwvZz2tP+2b8dSVs3/NEbRvI6mZMIAADALJebAgAALMaEXZwpkwgAAMAsmUQAAIBFqAnLJBokAivaQgqlwHLwXgRgrTBIBFiBVuqgY6n7tVKfJwDMy4RlEs1JBAAAYJZMIgAAwEK1yZuTKJMIAADALJlEAACAxZBJBAAAYFLJJAIrmuqXrBTeiwD0VMxJBAAAYIWpqnOq6r9X1S1VdXNV/fJM/NSq+lhV/f3Mv6cc7VgGiQAAAIvR2mhuc5tK8obW2lOTPDPJa6vqaUnelOS61tr5Sa6buT8ng0QAAIBVrrV2Z2vtczM/70pyS5InJHlJkitmdrsiyY8f7VjmJAIAACzCSpuTWFVPTPL0JH+b5IzW2p3JoYFkVT3uaI83SAQAAFgdtlfVDYfdv7y1dvnhO1TV1iQfSPL61tpDVTXvRgwSATjC89f95OA2VT4B4DAto1wn8b7W2oVDG6tqYw4NEP+stfZfZ8J3V9VZM1nEs5Lcc7RGDBIBWFYGlQCw/OpQyvAPk9zSWvvPh236UJJXJXnbzL9/cbRjGSQCAACsfs9J8sokX6yqm2Zib8mhweGVVfXzSb6ZZPhyoRljHSRW1SVJ/kuS9Un+oLX2tsdsr5ntL0qyN8mrv1OxBwAAYCWo6XH3IGmtfSrJ0ATEi+dzrLEtgVFV65O8I8mlSZ6W5Kdn1vE43KVJzp+5XZbknSPtJAAAwBozzkziRUlua619LUmq6n05tIbHlw/b5yVJ/ri11pL8TVWd/J1Jl3Md+NZbb81zn/vcR8Ve9rKX5Zd+6Zeyd+/evOhFLzriMa9+9avz6le/Ovfdd19e+tKXHrH9Na95TX7qp34qt99+e175ylcesf0Nb3hDfuzHfiy33nprfvEXf/GI7b/yK7+S5z3vebnpppvy+te//ojtv/mbv5lnP/vZ+fSnP523vOUtR2x/+9vfngsuuCAf//jH8xu/8RtHbP+93/u9POUpT8lf/uVf5rd/+7eP2P4nf/InOeecc/L+978/73znkWPtq666Ktu3b8973vOevOc97zli+9VXX50tW7bkd3/3d3PllVcesf0Tn/hEkuS3fuu38uEPf/hR2zZv3pyPfOQjSZJf//Vfz3XXXfeo7aeddlo+8IEPJEne/OY35zOf+cyjtp999tn50z/90yTJ61//+tx0002P2v7kJz85l19+qKjTZZddlq985SuP2n7BBRfk7W9/e5LkZ37mZ7Jjx45HbX/Ws56Vt771rUmSn/iJn8j999//qO0XX3xxfvVXfzVJcumll2bfvn2P2v7iF784b3zjG5PkiPdd4r3nvbc633ufbzfn7Hx3zqxz8nDbmy/ls7PbvtOW994nknjv+X/Pe+9w3nvee8nC3nur3gpbAmOxxpZJzKGFHW8/7P6Omdh890mSVNVlVXVDVd1w4MCBJe0oAADAWlGHknRjaLjqJ5O8sLX2L2fuvzLJRa21f33YPn+V5K0z19emqq5L8u9aazfOdewLL7yw3XDDDXPtAgAArBBVdeNcSzusZFtPPaf94MVHZo6Xw6eveuNIXqdxZhJ3JDnnsPtnJ7ljAfsAAACwRMY5SLw+yflV9aSq2pTk5Tm0hsfhPpTkZ+uQZyZ58GjzEQEAAEamJWltNLcRGVvhmtbaVFW9Lsk1ObQExrtbazdX1b+a2f6uJFfn0PIXt+XQEhg/N67+AgAArAVjXSextXZ1Dg0ED4+967CfW5LXjrpfAAAAx6pUNwUAAGBSjTWTCAAAsOrJJAIAADCpZBIBAAAWqGJOIgAAABNMJhEAAGChRryG4SjIJAIAADBLJhEAAGARzEkEAABgYskkAgAALIZMIgAAAJNKJhEAAGARzEkEAABgYskkAgAALFRLMj1ZqUSZRAAAAGbJJAIAACzGZCUSZRIBAAD4BzKJAAAAi6C6KQAAABNLJhEAAGAx2mSlEmUSAQAAmCWTCAAAsAjmJAIAADCxZBIBAAAWqsU6iQAAAEwumUQAAIAFqiSluikAAACTSiYRAABgMabH3YGlJZMIAADALINEAAAAZrncFAAAYBEUrgEAAGBiySQCAAAsVJu5TRCZRAAAAGbJJAIAACxYS8xJBAAAYFLJJAIAACxCTVYiUSYRAACAfyCTCAAAsBjmJAIAADCpZBIBAAAWqiU1Pe5OLC2ZRAAAAGbJJAIAACyGOYkAAABMKplEAACAxZisRKJMIgAAAP9AJhEAAGARypxEAAAAJpVMIgAAwGLIJAIAADCpZBIBAAAWqiWZHncnlpZMIgAAALNkEgEAABao0iauuulYBolVdWqS9yd5YpKvJ3lZa+2Bzn5fT7IrycEkU621C0fXSwAAgLVnXJebvinJda2185NcN3N/yA+31i4wQAQAAFak1kZzG5FxDRJfkuSKmZ+vSPLjY+oHAAAAhxnXIPGM1tqdSTLz7+MG9mtJrq2qG6vqspH1DgAA4FhNWCZx2eYkVtXHk5zZ2fTv53GY57TW7qiqxyX5WFX9XWvtkwPtXZbksiQ599xz591fAAAAlnGQ2Fp73tC2qrq7qs5qrd1ZVWcluWfgGHfM/HtPVX0wyUVJuoPE1trlSS5PkgsvvHCyygsBAAArk3USl8yHkrxq5udXJfmLx+5QVSdU1Ynf+TnJC5J8aWQ9BAAAWEWq6t1VdU9Vfemw2KlV9bGq+vuZf0852nHGNUh8W5LnV9XfJ3n+zP1U1eOr6uqZfc5I8qmq+nySzyb5q9baR8fSWwAAgAHV2khux+A9SS55TGw+K0skGdM6ia21+5Nc3InfkeRFMz9/LckPjrhrAAAAq1Jr7ZNV9cTHhF+S5LkzP1+R5BNJ/q+5jjOWQSIAAMDEGGHl0QV41MoSM0VB52SQCAAAsDpsr6obDrt/+UwBzyVlkAgAALBgI13D8L7W2oXzfMwxrSxxuHEVrgEAAGD5HXVliceSSQQAAFiolhUzJ7Gq/jyHitRsr6odSX4th1aSuLKqfj7JN5P85NGOY5AIAAAwAVprPz2w6YiVJebiclMAAABmySQCAAAsxvS4O7C0ZBIBAACYJZMIAACwCLVCCtcsFZlEAAAAZskkAgAALIZMIgAAAJNKJhEAAGChWpJpmUQAAAAmlEwiAADAgjVzEgEAAJhcMokAAACLIZMIAADApJJJBAAAWAyZRAAAACaVTCIAAMBCWScRAACASSaTCAAAsGAtadPj7sSSkkkEAABglkwiAADAYqhuCgAAwKSSSQQAAFgo1U0BAACYZDKJAAAAi2FOIgAAAJNKJhEAAGAxZBIBAACYVDKJAAAAC9ZkEgEAAJhcMokAAAAL1ZJMT4+7F0tKJhEAAIBZMokAAACLYU4iAAAAk0omEQAAYDFkEgEAAJhUMokAAAAL1pJpmUQAAAAmlEwiAADAQrWkNeskAgAAMKEMEgEAAJjlclMAAIDFULgGAACASSWTCAAAsBhNJhEAAIAJJZMIAACwUK0l05bAAAAAYELJJAIAACyGOYkAAABMqrEMEqvqJ6vq5qqarqoL59jvkqq6tapuq6o3jbKPAAAAx6JNT4/kNirjyiR+Kcm/SPLJoR2qan2SdyS5NMnTkvx0VT1tNN0DAABYm8YyJ7G1dkuSVNVcu12U5LbW2tdm9n1fkpck+fKydxAAAOCYNHMSR+gJSW4/7P6OmRgAAADLZNkyiVX18SRndjb9+9baXxzLITqxwSF6VV2W5LKZu49U1ZeOoQ0mx/Yk9427E4yc8772OOdrj3O+Njnva89Txt2BBWtJpicrk7hsg8Sydna1AAAFS0lEQVTW2vMWeYgdSc457P7ZSe6Yo73Lk1yeJFV1Q2ttsCAOk8c5X5uc97XHOV97nPO1yXlfe6rqhnH3gX+wktdJvD7J+VX1pCTfSvLyJK8Yb5cAAAAeo42u8ugojGsJjH9eVTuSPCvJX1XVNTPxx1fV1UnSWptK8rok1yS5JcmVrbWbx9FfAACAtWJc1U0/mOSDnfgdSV502P2rk1y9gCYuX3jvWKWc87XJeV97nPO1xzlfm5z3tWfVnvOWpE3YnMRqE1auFQAAYFS2rTutPXPDC0fS1scO/PmNo5ivu5LnJAIAAKxsrZmTuJJV1SVVdWtV3VZVbxp3f1h+VfXuqrrHkidrR1WdU1X/vapuqaqbq+qXx90nll9VHV9Vn62qz8+c9/973H1iNKpqfVX9r6r68Lj7wmhU1der6otVdZOKl2tDVZ1cVVdV1d/N/H1/1rj7tNZNTCaxqtYneUeS5+fQ8hnXV9WHWmtfHm/PWGbvSfI7Sf54zP1gdKaSvKG19rmqOjHJjVX1Mb/rE++RJD/SWttdVRuTfKqqPtJa+5txd4xl98s5VMBu27g7wkj9cGvNOolrx39J8tHW2kuralOSLePu0HxN2pzEScokXpTkttba11pr+5O8L8lLxtwnlllr7ZNJvj3ufjA6rbU7W2ufm/l5Vw59eHzCeHvFcmuH7J65u3HmNll/kTlCVZ2d5EeT/MG4+wIsj6raluSfJvnDJGmt7W+t7Rxvr5ikQeITktx+2P0d8cERJlpVPTHJ05P87Xh7wijMXHZ4U5J7knystea8T763J/l3SSZrsg9H05JcW1U3VtVl4+4My+67ktyb5I9mLi3/g6o6Ydydmrc2PZrbiEzSILE6Md8yw4Sqqq1JPpDk9a21h8bdH5Zfa+1ga+2CJGcnuaiqvm/cfWL5VNWLk9zTWrtx3H1h5J7TWntGkkuTvLaq/um4O8Sy2pDkGUne2Vp7epI9SdQWGbOJmZOYQ5nDcw67f3aSO8bUF2AZzcxJ+0CSP2ut/ddx94fRaq3trKpPJLkkiaJVk+s5Sf5ZVb0oyfFJtlXVn7bWfmbM/WKZzaybndbaPVX1wRyaUvTJ8faKZbQjyY7Drg65KqtskLgrD1zz8XbV9hE1N5K5upM0SLw+yflV9aQk30ry8iSvGG+XgKVWVZVD8xZuaa3953H3h9GoqtOTHJgZIG5O8rwk/3HM3WIZtdbenOTNSVJVz03yRgPEyTdzmeG61tqumZ9fkOT/GXO3WEattbuq6vaqekpr7dYkFydZVcXoWmuXjLsPS21iBomttamqel2Sa5KsT/Lu1trNY+4Wy6yq/jzJc5Nsr6odSX6ttfaH4+0Vy+w5SV6Z5Isz89OS5C2ttavH2CeW31lJrpipZL0uyZWtNUsiwOQ5I8kHD30fmA1J3tta++h4u8QI/OskfzZT2fRrSX5uzP1Z86o10/YAAAA4ZJIK1wAAALBIBokAAADMMkgEAABglkEiAAAAswwSAQAAmGWQCAAAwCyDRAAAAGYZJAKwolXVP66qL1TV8VV1QlXdXFXfN+5+AcCkqtbauPsAAHOqqt9IcnySzUl2tNbeOuYuAcDEMkgEYMWrqk1Jrk/ycJJnt9YOjrlLADCxXG4KwGpwapKtSU7MoYwiALBMZBIBWPGq6kNJ3pfkSUnOaq29bsxdAoCJtWHcHQCAuVTVzyaZaq29t6rWJ/l0Vf1Ia+2/jbtvADCJZBIBAACYZU4iAAAAswwSAQAAmGWQCAAAwCyDRAAAAGYZJAIAADDLIBEAAIBZBokAAADMMkgEAABg1v8P/ovNRdSE6DEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(14, 8))\n",
    "counts, xedges, yedges, im = ax.hist2d(x1_all, x2_all, bins=[120, 80], range= [[0.0, 6.0], [-1.0, 3.0]], cmin=1)\n",
    "ax.plot([0.0, 6.0], [0.0, 0.0], \"--k\")   # NOTE: This draws a line from [x1, x2], [y1, y2] with dashed line (\"--\") and in black (\"k\")\n",
    "fig.colorbar(im) # ticks=[-1, 0, 1]\n",
    "\n",
    "ax.set(title='Histogram of lengths (x) and widths (y)',\n",
    "       xlabel='x', \n",
    "       ylabel='y',\n",
    "       aspect='equal', # NOTE: This forces the x- and y-axis to have the SAME scale!!!\n",
    "      )\n",
    "\n",
    "d = {'Entries': len(x12_all),\n",
    "     'Mean x': x1_all.mean(),\n",
    "     'Mean y': x2_all.mean(),\n",
    "     'Std x': x1_all.std(ddof=1),\n",
    "     'Std y': x2_all.std(ddof=1),\n",
    "    }\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "add_text_to_ax(0.02, 0.97, text, ax, fontsize=15);\n",
    "\n",
    "fig.tight_layout()\n",
    "fig\n",
    "\n",
    "if save_plots :\n",
    "    fig.savefig(\"Dist_2Dgauss.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we bin `y_all` and fit it with a Gaussian distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian(x, N, mu, sigma):\n",
    "    return N * 1.0 / (sigma*np.sqrt(2*np.pi)) * np.exp(-0.5* (x-mu)**2/sigma**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAFlCAYAAAApo6aBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYE0lEQVR4nO3dXaxd51kn8P9DUopbqNxOncichHGQQjUJEi1jBUaWUMcBGsZV0wuCXAlkMRmFiwwUzYzoCReDuLB0pJEQSDNFilIYI0qD+agSEQSEmIohgganlClJmqlpTGInxKaDBQUrKOGZC69Oj+1jn2P7Pd77+Px+krXXeva79n7Okj/+Xu/6qO4OAABX7mtm3QAAwLVCsAIAGESwAgAYRLACABhEsAIAGESwAgAY5PpZN5Ak73znO3vHjh2zbgMAYFVPP/3033T3tpXeWzVYVdW7kvzqstI3J/mvSX5pqu9IcjTJD3T3307bPJDk3iRvJPmx7v7di33Hjh07cvjw4VV/EACAWauqv7rQe6tOBXb389397u5+d5J/neQfk3wyyWKSJ7r71iRPTOupqtuS7E1ye5K7kny0qq674p8CAGDOXeo5Vncm+cvu/qskdyc5MNUPJPngtHx3koe7+7XufiHJkSR3jGgWAGCeXWqw2pvkE9Pyjd39SpJMrzdM9YUkLy3b5thUO0tV3VdVh6vq8MmTJy+xDQCA+bPmYFVVX5vkA0l+bbWhK9TOeyBhdz/Y3Tu7e+e2bSue/wUAsKFcyhGr70vyme5+dVp/taq2J8n0emKqH0ty87Ltbkry8pU2CgAw7y4lWH0oX50GTJJHk+yblvcleWRZfW9Vvbmqbklya5KnrrRRAIB5t6b7WFXVW5J8T5IfWVZeSnKwqu5N8mKSe5Kku5+pqoNJnk3yepL7u/uNoV0DAMyhNQWr7v7HJP/inNqXcuYqwZXG70+y/4q7AwDYQDzSBgBgEMEKAGAQwQoAYBDBCgBgEMEKAGCQNV0VCLAZ7Vo6lOOnTp9VW9i6JU8u7p5RR8C8E6wALuD4qdM5urTnrNqOxcdm1A2wEZgKBAAYxBEr4JplKg+42gQr4JplKg+42kwFAgAM4ogVsKksbN1y3lEr04PAKIIVsKmsFKBMDwKjmAoEABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABjEfayATW+lm4Z+pQ5wKQQrYNNz13VgFMEKuCbsWjqU46dOn1VzxAm42gQr4Jpw/NTpHF3aM+s2gE3OyesAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAINcv5ZBVbU1yUNJvjVJJ/n3SZ5P8qtJdiQ5muQHuvtvp/EPJLk3yRtJfqy7f3d048DmtWvpUI6fOn1WbWHrlhl1A/BVawpWSX4uye909/dX1dcmeUuSn0zyRHcvVdViksUkH6mq25LsTXJ7km9M8vtV9S3d/cY69A9sQsdPnc7RpT2zbgPgPKtOBVbV25J8V5KPJUl3/1N3n0pyd5ID07ADST44Ld+d5OHufq27X0hyJMkdoxsHAJg3aznH6puTnEzyi1X1Z1X1UFW9NcmN3f1KkkyvN0zjF5K8tGz7Y1PtLFV1X1UdrqrDJ0+evKIfAgBgHqwlWF2f5NuT/Hx3vyfJP+TMtN+F1Aq1Pq/Q/WB37+zundu2bVtTswAA82wtwepYkmPd/elp/ddzJmi9WlXbk2R6PbFs/M3Ltr8pyctj2gUAmF+rBqvu/uskL1XVu6bSnUmeTfJokn1TbV+SR6blR5Psrao3V9UtSW5N8tTQrgEA5tBarwr80SQfn64I/GKSH86ZUHawqu5N8mKSe5Kku5+pqoM5E75eT3K/KwIBgM1gTcGquz+bZOcKb915gfH7k+y/gr4AADYcd14HABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhkrTcIBSDJwtYt2bH42Hm1Jxd3z6gjYJ4IVgCXYKUAdW7QAjYvU4EAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAg1w/6wYALmbX0qEcP3X6rNrC1i0z6mZlC1u3ZMfiY+fVnlzcPaOOgFkRrIC5dvzU6Rxd2jPrNi5qpQB1btACNgdTgQAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIOsKVhV1dGq+lxVfbaqDk+1d1TV41X1hen17cvGP1BVR6rq+ap633o1DwAwTy7liNW/7e53d/fOaX0xyRPdfWuSJ6b1VNVtSfYmuT3JXUk+WlXXDewZAGAuXclU4N1JDkzLB5J8cFn94e5+rbtfSHIkyR1X8D0AABvCWoNVJ/m9qnq6qu6bajd29ytJMr3eMNUXkry0bNtjU+0sVXVfVR2uqsMnT568vO4BAObIWh9ps6u7X66qG5I8XlWfv8jYWqHW5xW6H0zyYJLs3LnzvPeBzWcjPBcQ4GLWFKy6++Xp9URVfTJnpvZerart3f1KVW1PcmIafizJzcs2vynJywN7Bq5RG+G5gAAXs+pUYFW9taq+4SvLSb43yV8keTTJvmnYviSPTMuPJtlbVW+uqluS3JrkqdGNAwDMm7UcsboxySer6ivjf6W7f6eq/jTJwaq6N8mLSe5Jku5+pqoOJnk2yetJ7u/uN9alewCAObJqsOruLyb5thXqX0py5wW22Z9k/xV3BwCwgbjzOgDAIIIVAMAgghUAwCCCFQDAIIIVAMAgghUAwCCCFQDAIIIVAMAgghUAwCCCFQDAIIIVAMAgghUAwCCCFQDAIIIVAMAgghUAwCCCFQDAIIIVAMAgghUAwCCCFQDAIIIVAMAgghUAwCCCFQDAIIIVAMAg18+6AYBr0cLWLdmx+Nh5tScXd8+oI+BqEKwA1sFKAercoAVce0wFAgAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAM4j5WwEzsWjqU46dOn1Vb2LplRt0AjCFYATNx/NTpHF3aM+s2AIYyFQgAMIhgBQAwiGAFADDImoNVVV1XVX9WVb81rb+jqh6vqi9Mr29fNvaBqjpSVc9X1fvWo3EAgHlzKUesPpzkuWXri0me6O5bkzwxraeqbkuyN8ntSe5K8tGqum5MuwAA82tNwaqqbkqyJ8lDy8p3JzkwLR9I8sFl9Ye7+7XufiHJkSR3jGkXAGB+rfWI1c8m+Ykk/7ysdmN3v5Ik0+sNU30hyUvLxh2bamepqvuq6nBVHT558uQlNw4AMG9WDVZV9f4kJ7r76TV+Zq1Q6/MK3Q92987u3rlt27Y1fjQAwPxayw1CdyX5QFX9uyRfl+RtVfXLSV6tqu3d/UpVbU9yYhp/LMnNy7a/KcnLI5sG2IgWtm7JjsXHzqs9ubh7Rh0Bo60arLr7gSQPJElVvTfJf+nuH6yq/5ZkX5Kl6fWRaZNHk/xKVf1Mkm9McmuSp8a3DrCxrBSgzg1awMZ2JY+0WUpysKruTfJiknuSpLufqaqDSZ5N8nqS+7v7jSvuFABgzl1SsOruTyX51LT8pSR3XmDc/iT7r7A3AIANxZ3XAQAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABrl+1g0A175dS4dy/NTps2oLW7fMqBuA9SNYAevu+KnTObq0Z9ZtAKw7U4EAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAg3hWIMAMLWzdkh2Lj51Xe3Jx94w6Aq6EYAUwQysFqHODFrBxmAoEABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhk1WBVVV9XVU9V1Z9X1TNV9dNT/R1V9XhVfWF6ffuybR6oqiNV9XxVvW89fwAAgHmxliNWryXZ3d3fluTdSe6qqu9Mspjkie6+NckT03qq6rYke5PcnuSuJB+tquvWo3kAgHmyarDqM748rb5p+tVJ7k5yYKofSPLBafnuJA9392vd/UKSI0nuGNo1AMAcWtM5VlV1XVV9NsmJJI9396eT3NjdryTJ9HrDNHwhyUvLNj821c79zPuq6nBVHT558uSV/AwAAHNhTcGqu9/o7ncnuSnJHVX1rRcZXit9xAqf+WB37+zundu2bVtbtwAAc+ySrgrs7lNJPpUz5069WlXbk2R6PTENO5bk5mWb3ZTk5SvuFABgzq3lqsBtVbV1Wt6S5LuTfD7Jo0n2TcP2JXlkWn40yd6qenNV3ZLk1iRPjW4cAGDeXL+GMduTHJiu7PuaJAe7+7eq6o+THKyqe5O8mOSeJOnuZ6rqYJJnk7ye5P7ufmN92gcAmB+rBqvu/t9J3rNC/UtJ7rzANvuT7L/i7gAANhB3XgcAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGGQtzwoEWLNdS4dy/NTps2oLW7fMqJuNaWHrluxYfOy82pOLu2fUEbBWghUw1PFTp3N0ac+s29jQVgpQ5wYtYD6ZCgQAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABjEQ5iBy7Zr6VCOnzp9Vm1h65YZdQMwe4IVcNmOnzqdo0t7Zt0GwNwwFQgAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMMiqwaqqbq6qP6iq56rqmar68FR/R1U9XlVfmF7fvmybB6rqSFU9X1XvW88fAABgXqzliNXrSf5zd/+rJN+Z5P6qui3JYpInuvvWJE9M65ne25vk9iR3JfloVV23Hs0DAMyTVR/C3N2vJHllWv77qnouyUKSu5O8dxp2IMmnknxkqj/c3a8leaGqjiS5I8kfj24eYLNY2LolOxYfO6/25OLuGXUErGTVYLVcVe1I8p4kn05y4xS60t2vVNUN07CFJH+ybLNjUw2Ay7RSgDo3aAGzt+aT16vq65P8RpIf7+6/u9jQFWq9wufdV1WHq+rwyZMn19oGAMDcWlOwqqo35Uyo+nh3/+ZUfrWqtk/vb09yYqofS3Lzss1vSvLyuZ/Z3Q92987u3rlt27bL7R8AYG6s5arASvKxJM91988se+vRJPum5X1JHllW31tVb66qW5LcmuSpcS0DAMyntZxjtSvJDyX5XFV9dqr9ZJKlJAer6t4kLya5J0m6+5mqOpjk2Zy5ovD+7n5jeOcAAHNmLVcF/lFWPm8qSe68wDb7k+y/gr4AADYcd14HABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGOT6WTcAwOVZ2LolOxYfO6/25OLuGXUECFYAG9RKAercoAVcXYIVsCa7lg7l+KnTZ9UWtm6ZUTcA80mwAtbk+KnTObq0Z9ZtAMw1J68DAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAziPlbAedwMFODyCFbAedwMFODymAoEABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhk1WBVVb9QVSeq6i+W1d5RVY9X1Rem17cve++BqjpSVc9X1fvWq3EAgHmzlmcF/s8k/z3JLy2rLSZ5oruXqmpxWv9IVd2WZG+S25N8Y5Lfr6pv6e43xrYNjOKBywDjrBqsuvsPq2rHOeW7k7x3Wj6Q5FNJPjLVH+7u15K8UFVHktyR5I/HtAuM5oHLAOOs5YjVSm7s7leSpLtfqaobpvpCkj9ZNu7YVAPgKljYuiU7Fh87r/bk4u4ZdQSby+UGqwupFWq94sCq+5LclyTf9E3fNLgNgM1ppQB1btAC1s/lXhX4alVtT5Lp9cRUP5bk5mXjbkry8kof0N0PdvfO7t65bdu2y2wDAGB+XG6wejTJvml5X5JHltX3VtWbq+qWJLcmeerKWgQA2BhWnQqsqk/kzInq76yqY0l+KslSkoNVdW+SF5PckyTd/UxVHUzybJLXk9zvikAAYLNYy1WBH7rAW3deYPz+JPuvpCkAgI1o9MnrAMwZVwrC1SNYAVzjXCkIV49nBQIADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAzizuuwSexaOpTjp06fV1/YumUG3QBcmwQr2CSOnzqdo0t7Zt0GwDXNVCAAwCCCFQDAIKYCATahha1bsmPxsfNqTy7unlFHcG0QrAA2oZUC1LlBC7h0pgIBAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAZxHysAkrhpKIwgWME1aNfSoRw/dfqs2sLWLTPqho3CTUPhyglWcA06fup0ji7tmXUbAJuOc6wAAAZxxAo2ONN+rKeVzrv6St25V3A+wQo2ONN+rKcLhSfnXsHKTAUCAAwiWAEADCJYAQAM4hwrmFMrnZS+EieqA8wPwQrmlJPSATYewQqAS3ah2zCsNM5tGdhMBCsALtlaw5LbMrDZrFuwqqq7kvxckuuSPNTdS+v1XQDMJw92ZrNZl2BVVdcl+R9JvifJsSR/WlWPdvez6/F9sJE4KZ3NxIOd2WzW64jVHUmOdPcXk6SqHk5ydxLBik3PSelsdo5icS1br2C1kOSlZevHknzHOn0XrKtLOcLkHwZY3Up/TnYtHbrsk+Ev9LxMfx6Zheru8R9adU+S93X3f5jWfyjJHd39o8vG3Jfkvmn1XUmeH97I+d6Z5G+uwvdsVPbP6uyji7N/VmcfXZz9c3H2z+quxj76l929baU31uuI1bEkNy9bvynJy8sHdPeDSR5cp+9fUVUd7u6dV/M7NxL7Z3X20cXZP6uzjy7O/rk4+2d1s95H6/VImz9NcmtV3VJVX5tkb5JH1+m7AADmwrocseru16vqPyb53Zy53cIvdPcz6/FdAADzYt3uY9Xdv53kt9fr8y/TVZ163IDsn9XZRxdn/6zOPro4++fi7J/VzXQfrcvJ6wAAm9F6nWMFALDpbIpgVVV3VdXzVXWkqhZn3c+8qapfqKoTVfUXs+5lHlXVzVX1B1X1XFU9U1UfnnVP86aqvq6qnqqqP5/20U/Puqd5VFXXVdWfVdVvzbqXeVRVR6vqc1X12ao6POt+5k1Vba2qX6+qz09/H/2bWfc0T6rqXdPvna/8+ruq+vGr3se1PhU4PV7n/2TZ43WSfMjjdb6qqr4ryZeT/FJ3f+us+5k3VbU9yfbu/kxVfUOSp5N80O+hr6qqSvLW7v5yVb0pyR8l+XB3/8mMW5srVfWfkuxM8rbufv+s+5k3VXU0yc7udp+mFVTVgST/q7sfmq64f0t3n5p1X/No+rf/eJLv6O6/uprfvRmOWP3/x+t09z8l+crjdZh09x8m+b+z7mNedfcr3f2ZafnvkzyXM08XYNJnfHlafdP069r+X9slqqqbkuxJ8tCse2Hjqaq3JfmuJB9Lku7+J6Hqou5M8pdXO1QlmyNYrfR4Hf8oclmqakeS9yT59Gw7mT/TNNdnk5xI8nh320dn+9kkP5Hkn2fdyBzrJL9XVU9PT+fgq745yckkvzhNJz9UVW+ddVNzbG+ST8ziizdDsKoVav4nzSWrqq9P8htJfry7/27W/cyb7n6ju9+dM09auKOqTCtPqur9SU5099Oz7mXO7erub0/yfUnun05T4Izrk3x7kp/v7vck+YckzhlewTRN+oEkvzaL798MwWrVx+vAaqbzhn4jyce7+zdn3c88m6YnPpXkrhm3Mk92JfnAdA7Rw0l2V9Uvz7al+dPdL0+vJ5J8MmdO5eCMY0mOLTsS/Os5E7Q43/cl+Ux3vzqLL98Mwcrjdbgi04nZH0vyXHf/zKz7mUdVta2qtk7LW5J8d5LPz7ar+dHdD3T3Td29I2f+DjrU3T8447bmSlW9dbo4JNMU1/cmcaXypLv/OslLVfWuqXRnEhfQrOxDmdE0YLKOd16fFx6vs7qq+kSS9yZ5Z1UdS/JT3f2x2XY1V3Yl+aEkn5vOIUqSn5yeLsAZ25McmK7E+ZokB7vbLQW4FDcm+eSZ/8fk+iS/0t2/M9uW5s6PJvn4dJDgi0l+eMb9zJ2qekvO3AXgR2bWw7V+uwUAgKtlM0wFAgBcFYIVAMAgghUAwCCCFQDAIIIVAMAgghUAwCCCFQDAIIIVAMAg/w8qNN7wZbQ55QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xmin, xmax = 0.0, 7.0\n",
    "\n",
    "fig2, ax2 = plt.subplots(figsize=(10, 6));\n",
    "counts, bin_edges, _ = ax2.hist(y_all, 100, range=(xmin, xmax), histtype='step');\n",
    "bin_centers = (bin_edges[1:] + bin_edges[:-1])/2\n",
    "s_counts = np.sqrt(counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the distribution of \"whatever you put into it\" (initially x1-2*x2), which shows what output you get and what uncertainty to expect (given by the width - think about this!). We can thus get the result by simply recording the mean and width (SD):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Mean = 3.593,    SD = 0.415\n"
     ]
    }
   ],
   "source": [
    "print(f\"  Mean = {y_all.mean():5.3f},    SD = {y_all.std(ddof=1):5.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, we are in principle not even sure, if this distribution is Gaussian, so in order to check this, we fit it. Note that in the second line (defining the Minuit object to minimised) we give \"reasonable\" starting values, which is key to making fitting work (more on that later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit value: N = 698.28075 +/- 6.99326\n",
      "Fit value: mu = 3.59303 +/- 0.00416\n",
      "Fit value: sigma = 0.41390 +/- 0.00298\n"
     ]
    }
   ],
   "source": [
    "Chi2_object = Chi2Regression(gaussian, bin_centers[counts>0], counts[counts>0], s_counts[counts>0])\n",
    "minuit = Minuit(Chi2_object, pedantic=False, N=5000, mu=3.5, sigma=0.4, limit_sigma=(0, 10000), print_level=0) #   \n",
    "minuit.migrad()                          # Perform the actual fit\n",
    "\n",
    "for name in minuit.parameters:\n",
    "    print(\"Fit value: {0} = {1:.5f} +/- {2:.5f}\".format(name, minuit.values[name], minuit.errors[name]))\n",
    "\n",
    "chi2 = minuit.fval\n",
    "N_var = 3                                # Number of variables (N, mu, sigma)\n",
    "N_dof = len(counts[counts>0]) - N_var    # Number of degrees of freedom\n",
    "\n",
    "from scipy import stats\n",
    "chi2_prob =  stats.chi2.sf(chi2, N_dof) # The chi2 probability given N_DOF degrees of freedom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot the fit on top of the histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "xaxis = np.linspace(xmin, xmax, 1000)\n",
    "yaxis = gaussian(xaxis, *minuit.args)\n",
    "ax2.plot(xaxis, yaxis);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXhN1/rA8e/OdDJVBkkkQiREKYkYgppCy71FWkJpSwdC65pLQ+nwQ6lZW1XVFhVjNZWrpbTpoFzRagkiUUMNSRFBDEFE5vX7I8mpI0ESSY4c7+d59nPOWWvttd59QvvaWXstTSmFEEIIIYQQIp+ZsQMQQgghhBDifiIJshBCCCGEEDeRBFkIIYQQQoibSIIshBBCCCHETSRBFkIIIYQQ4iYWxg4AwMXFRXl7exs7DCGEEEII8QDZs2fPBaWU663l90WC7O3tTUxMjLHDEEIIIYQQDxBN0/4urlymWAghhBBCCHGTuybImqY10DQt9qbjqqZpYzRNc9Y07SdN044WvDrddM4bmqYd0zTtiKZpT1TsJQghhBBCCFF+7pogK6WOKKWaKqWaAi2AdOBrYCKwRSlVH9hS8BlN0xoBzwGNga7AIk3TzCsofiGEEEIIIcpVaadYdAaOK6X+BnoCKwrKVwAhBe97Al8qpTKVUgnAMaBVeQQrhBBCCCFERSttgvwcsLbgfQ2lVDJAwatbQbkncOqmc04XlBnQNG2IpmkxmqbFpKSklDIMIYQQQgghKkaJE2RN06yAHsC6uzUtpkwVKVBqsVIqUCkV6OpaZHUNIYQQQgghjKI0d5C7AXuVUucKPp/TNM0DoOD1fEH5aaD2TefVAs7ca6BCCCGEEEJUhtIkyP34Z3oFwEZgQMH7AcCGm8qf0zRNp2maD1Af2HWvgQohhBBCCFEZSrRRiKZptsC/gP/cVDwL+ErTtMHASaAvgFLqT03TvgIOAjnACKVUbrlGLYQQQgghRAUpUYKslEoHqt9SdpH8VS2Kaz8dmH7P0QkhhBBCCFHJZCc9IYQQQgghbiIJcjnr1KkTI0eONHYYQgghhBCijKp8gjxw4EA0TStyLF26tFT9LF++HHt7+3uOZ/369cycOfOe+xEPhujoaHr06EGNGjXQNI1t27YVafP111/ToEEDrK2tCQwMZPfu3Qb1SUlJBAcHY2trS82aNZk9e3aRPt555x3c3d2xs7MjJCSE8+fPG9Tv2LGDZs2aYW1tTaNGjfj+++/L9TqFEEKIqqTKJ8gAXbp0ITk52eB4/vnnjRKLs7MzDz30kFHGFlXP1atX8fPzY86cOcXWHzhwgGeeeYYXXniBvXv30qRJE7p3786VK1f0bfr06cPVq1f57bffeO+995g0aRJr1qzR13/22WfMmTOHRYsWER0dTVJSksHfj5SUFIKDg2nbti179+4lJCSEXr16ceLEiYq7cCGEEOJ+ppQy+tGiRQtVVgMGDFDBwcG3rd+6dasCVEREhGrQoIF66KGH1IABA1R2drZSSqnw8HBF/kYmBkfHjh0N+gkPD1d2dnbq119/VU2bNlU6nU75+vqqCxcuKKWUeuKJJ/Tnjhgxokgcubm5asqUKapWrVrKzs5OdejQQcXGxhq0WbdunWrcuLGytrZW7u7uauDAgWX+XkTVkpKSogC1detWg/JXX31VNWvWTP85IyND2dvbqyVLliillNq3b58C1P79+/VtXn75ZdWuXTv954CAADV27Fj95x07dihAHT16VCml1AcffKBcXFxUTk6Ovk39+vXVW2+9Va7XKIQQQtxvgBhVTG5qEneQS2Lp0qVERkayevVqVq9eTWRkJADPPvssycnJzJ8/H1tbW/0d6PXr1xfpIycnh7CwMGbMmMGBAwd444039HVr164lOTmZNm3aFDv+1KlTWb16NeHh4cTGxtK+fXueeOIJrl27BsC5c+fo378//fv359ChQ3z77bc0aNCgAr4JUZXs3r2bDh066D/rdDpatWpFTEyMvt7BwYEmTZro2wQFBbF3716UUmRmZhIfH2/QR5s2bbC0tDToo23btpibm+vbdOjQQV8vhBBCPGhMIkGOiorC3t7e4IiPjzdoM2nSJPz8/OjRowetWrXSz+O0sbHB3d0dBwcHNE3D3d0dd3d3nJ2di4yTmZnJ9OnT6datG76+vgwaNIjq1fNXv3NycsLd3R0rK6si52VkZDBnzhwWLFhAly5d8PX1Zfr06Sil2Lx5MwCnT58mOzubHj164O3tTWBgIBMnTizvr0pUMSkpKbi4uPDNN9/g5ORETEwMLi4upKSk6OurV69Oeno6Pj4+vP7667i4uHDjxg3S0tK4ePEieXl5uLi4MG7cOOrWrUtGRgZOTk4Gfbi4uLBr1y4cHR3ZsGGDwRhCCCHEg6ZE6yDf74KCgli8eLFBmZeXl8FnX19f/XtnZ2cuXbpU6nE0TaNdu3alPu/o0aPcuHGDPn36oGmavvzGjRv6eZ5NmjShbdu2dOjQgX//+9+0bduWfv364ebmVurxhOlQSqFpGg899BBeXl7Y2tqS/xshw3pzc3Nq166Nm5tbkXrI/7Pr6uqKl5cXFhYWxfZhZ2dHnTp1cHBwMKgXD452s34hKfVGkXJPRxt+nfi4ESISQgjjMIkE2dbW1iABLo6FheGlliUBsLW1RafTlfq8Qps2baJ27doGZYV3qi0tLYmOjua3335j27ZtLFq0iNmzZ/Pnn3/i5ORU5jFF1ebm5kZKSgqdO3dm//79AFy8eJH69evr6y9cuIBOp2P79u0ArFq1ChsbG+zt7bGyssLMzIyUlBQmTJjAhAkTyMvLIzU1FVdXV4MxGjdurB9j1apV+npRBeXmwO6lsH8tXEoAOxfw7wNtRoJ1tduelpR6g8RZwUXKvSdurshohRDivmMSUyzKg5WVFTk5ORXSd/369bG2tiY5ORlfX1+D4+apHGZmZrRv3563336b6OhokpOT2bVrV4XEJKqGwMBAoqOj9Z8zMzPZtWsXgYGB+vorV64QFxenb7N9+3aaN2+OpmnodDr8/PwM+ti5cyfZ2dkGffz222/k5v6zI3x0dLS+XlQxaSmw7N8QNQHMLCDgWXD0gv/Ngc+C4OwBY0cohBD3PZO4g5yZmcnZs2cNygrnIpeUr68vmZmZbNiwgX/9619YWFgUO5+4OFlZWfopG1lZWaSnp+vjcXd3x9ramvHjxxMWFoalpSUtWrTg5MmTfPnll4wZM4aGDRuye/dufvjhB7p160b16tVZvXo1lpaWNGzYsMTXIKqetLQ0jh07RmpqKgDHjh3D0dERLy8vnJ2deeWVV1i0aBHvvvsuvXv3Zt68eeh0Ovr27QtAs2bNePTRRxkxYgQLFy7k0KFDrFy5kmXLlunHGDZsGGFhYbRv3x5vb2/GjBmjnwsP8PzzzzNlyhRGjx7NyJEjWbVqFSdPniQ0NLTyvxBRajdPi3DkGhFW0/DSzjMuezSxFx/j11cKpkac/B3WDYSVPWHQD+By59+6CSHEA624pS0q+7jXZd4oZpm2wiWqCpd5S0lJ0Z8THBysBgwYUKSvsWPHKldX1zsu81acwjGKOwrl5OSoqVOnKm9vb2Vpaam8vLxUaGioPq4jR46orl27KhcXF2VjY6MCAgLU119/XebvRVQNt/uzEx4erm/z3//+V9WvX19ZWVmpFi1aqD/++MOgj1OnTqlu3brplwecOXNmkXEmT56s3NzclI2NjerZs6c6e/asQf327dtVQECAsrKyUg0bNlSbN2+ukOsV5a/OhE35b3JzlVrdV6mpLkod32ZYVyjlL6Vm11VqQXOlMtNu31cJy4UQoqrjNsu8aeo+eBgnMDBQyZJSQghRet4TN+fPG/790/xpFd3nQatXDOtulrAdVvSAwFB48oPi+7rdGEIIYWI0TdujlCoyp1DmIAshRFV37Rz88i74doGWL9+5rU8QtBkBMcvgtNyYEEKI4pjEHGQhhHigbXkHcjOh2xy4aSlJT0ebYleg8HVoy8/26yBqIgz+yeAcIYQQkiALIUSV5qMl5y/n9uhwqF7PoO52axd7T9wMz/wfbBwJhzfDI09WRqhCCFFlyBQLIYSowkZafAPmOmj3aulODOgHTt4Q/R7cB8+iCCHE/UQSZCGEqKqunKan2a8QOAjsS7nrprkFtBsDZ/bCiW0VEp4QQlRVVTpBXr58OZqm0atXL31ZbGwsmqbh7e1tvMCEKKGhQ4dSr149rK2t8fT0ZPjw4Vy9erVUfWiaVuQ4cOCfzSCSk5P125bb29vz1FNPcerUKX39jh07aNeuHdWrV8fe3p7WrVsTFRVVbtcoKtDuz9FQ8OjQsp3ftD/YucKuxeUblxBCVHFVOkEGsLGxISYmRp9UREREFNnOWYj7VUBAACtWrODw4cOsW7eO6Ohohg8fXup+lixZQnJysv64eYOZF198kYSEBKKioti5cyfp6emEhITot1u3trZm5MiRbN++nfj4eEJCQujZsycHDx4st+sUFSA7A/au4Oe8Fvk75ZWFhQ6avwR/RUHqyfKNTwghqrAqnyCbmZkRHBzMhg0bAIiMjKRPnz4Gbfbv30/nzp2xtbWlTp06TJo0yWBb6ffff58mTZpgZ2dH9erVeemll/Q74wFMmTIFPz8/PvzwQ2rUqIGbmxsffvhh5VygMGnDhg3T73DXtm1bhgwZwi+//FLqfhwdHXF3d9cfFhb5z9+mp6fzyy+/MGnSJJo3b46/vz/z589n79697NmzB8jfarpfv340btwYHx8f3njjDezt7dmxY0e5XqsoZ4c3QfpFVub+6976aVGwY+Ke5fcckhBCmIoqnyADPPvss0RERLB79248PDyoWbOmvu7ixYs8/vjjBAYGsn//flatWsWaNWt4//339W0uX77MtGnT2L9/P99//z0HDhxg2LBhBmOcOHGCo0ePEh0drd+69+ZfUwtxr86cOcNXX31FixYtSn3umDFjcHV1pU2bNnzzzTf68uzsbJRS6HQ6fZm1tTUA+/btK9JPbm4ua9asITU1lWbNmpXhKkSliYuAarX4La/xvfXjWDt//eT9EWjklU9sQghRxZlEgtyxY0fi4+P55JNPePbZZw3qFi5ciK+vL7Nnz6Z+/foEBQURFhbGsmXL9G2mTZtGz5498fX1pVWrVowaNYrvv//eoB8zMzPee+89Hn74Yd58803y8vKIjY2tlOsTpm3RokXY2tri6emJg4MDX3zxRanOnzFjBpGRkXz//fe0b9+eXr168eOPPwLg4OBAQEAACxYs4OrVq1y/fp1p06ZhYWFBSkqKQT+1atVCp9MxbNgw1q9fT8uWLcvtGkU5SzsPx7ZAk2dQ5fGfcf9n4OppWmlH7r0vIYQwASaRIJuZmdGzZ09WrlxZZHpFXFwce/fuxd7eXn+EhYVx4sQJfZtffvmFLl264Onpib29PcOHDyctLc2gn8LkAUCn02Fra2swDUOIsnr++eeJjY3l22+/5cSJE7zzzjulOv+NN97g0UcfJTAwkLlz59KjRw+DKUArV64kISEBR0dHnJyccHFxwdXVFTMzw7/+0dHR7N69m1GjRjFs2DCOHj1aLtcnKkB8JKhcCHiufPpr2B0s7ehpLtNqhBACTGijkNGjR9OyZUtq1KhRpO7JJ59k7ty5xZ73999/ExwczKBBg5g1axaOjo5s3LiRsLAwg3aFczpvpmTtUFEOHBwccHBw4OGHH8bBwYGOHTsybtw43N3dy9RfixYtiIiI0H9u0qQJcXFxXLp0CTMzM3Q6HR999JHBVCQAHx8fAJo1a8bvv//OBx98wKJFi8p+YaLi7F8LHk3BtQFw7N77s7KDhsEEx22CnMz8h/eEEOIBZjIJsq+vL76+vkXK/f39iYiIoG7dukXumAHExMSQlZXFggULMDc3B/LnggphDObm5iiluHHjhkF5WloaFy5cwMXFBXt7+zv2cfjw4WKXOXR2dgby7yjn5uYSFBR0xziuX79e+gsQFe/icTgbB0/MLHMXxW1B3dHMhxVW6XD0J9lZTwjxwDOZBPl2Ro4cyUcffcSgQYMYO3YslpaW/Pbbb5w8eZKpU6dSv3598vLyWLRoEU8++STbt29n1apVxg5bPACOHz/OsmXLCA4OxtPTk4SEBMLCwmjatKn+bm6hyMhIQkNDCQ8PZ+DAgfryjRs3kpSURPv27dHpdGzYsIGIiAg2bdqkb7Nz506uXr1Kw4YN2bdvH+PHj2fgwIH6JPr999+nZs2aBAQEYG5uzvr16/n5559Zv359ZXwNorQOF/xs7yGJLXYL6twnYN4SOLRREmQhxAPP5BNkFxcXtmzZwoQJE2jXrh3m5ub4+fkxcuRIIP/Xzx9++CGzZ89mwoQJPPbYY/zf//0fo0aNMnLkwtTZ2NgQGxvLsmXLuHTpEi4uLnTp0oUZM2aUuA9LS0sWLlzI+PHjAXjkkUdYt24d3bp107e5ceMGw4cP59SpU7i5uREaGsrUqVP19VZWVkyfPp3ExESUUtSvX58VK1YQEhJSfhcrys/hzeARUPa1j2/H3AIefgKOfA+52WBuWb79CyFEFaLdD/NoAwMDVUxMjLHDEEKI+1a7Wb+QlZrMH7oRfJDzNB/l9gbyp0sUe0e4LA5uhK9ehAHfgs8/U3C8J24mcVZw+YwhhBD3EU3T9iilAm8tN/k7yEIIYQqSUm+Q2CcLNinCRr9GWI17XP+4OPUeB3Nd/l1kn9vPURdCCFNnEsu8CSHEA+HwZnDyBrdGFdO/zh7qdswf5z747aIQQhiLJMhCCFEFWJMJCduhQXfQtIobqEF3SP0bzh+suDGEEOI+JwmyEEJUAa3NDkNuJvh2rtiBHu6a/3rs54odRwgh7mOSIBdj+fLld11rVgghKlMHs7j8+cFebSt2oGoe4PoIHP+lYscRQoj7WJVPkG/cuMHEiRPx9vbG2tqaunXrMnjw4CLtOnXqpF/aTYj7RXZ2NiNHjsTZ2Zlq1aoRGhpa5g06rly5gpeXV5F/3B08eJA+ffrg5eWFpmksX77coH7btm1ommZwuLi4lPWSRAUJMouDOm3ByrbiB6v3GPy9E7Jv3L2tEEKYoCqfIL/66qt88803fP755xw6dIhPP/2UnJwcY4clRIlMmjSJdevWERkZSVRUFNu2bWPMmDFl6mvEiBFFNhgBuHbtGt7e3sybNw8bG5vbnv/nn3+SnJxMcnIyBw/K/NP7ypXTPGyWVPHTKwrVezx/OsfJnZUznhBC3GeqfIL83//+l9dff53OnTvj4+PDv//9b1asWKGv79SpE5qm8b///Y+PP/5Yf4fs5rtox44dIygoCGtra1q2bMlff/1lhCsRD5rc3FwWL17MW2+9xeOPP07btm159913WbVqVanvIkdERHDu3DlCQ0OL1LVu3Zp58+bxzDPPFLvdeiE3Nzfc3d1xd3fHzc2t1NcjKtDxrfmv9cppveO7qdMWzK1kmoUQ4oFV5RPkatWq8eOPP5KRkVFs/fr160lOTqZNmzaEhobq75A9++yz+jYvvPACNjY27NmzhwkTJrBw4cLKCl88wE6cOMGlS5fo0KGDviwoKIjMzEzi4+NL3E9SUhLjx49nyZIl9xRPy5YtqVmzJsHBwRw4cOCe+hLl7PgvnFOOFbe8262s7KB2azi+rXLGE0KI+0yVT5AXLlzIjz/+iKurK8HBwXz88cdcuXJFX+/s7Iy7uztWVlbY2trq75AV/qo5Pj6eP/74g/nz59O4cWP69OlDv379jHU54gGSkpIC5G+H3rdvX1q3bq2f+1tYdzdKKQYOHMjrr7+Ot7d3meLw8PBg6dKlfP3116xdu5a8vDyCgoI4d+5cmfoT5UwpSIzmt7zGFbu8263qPQbn4iHtfOWNKYQQ94kqnyAHBwdz8uRJVqxYQb169Zg5cyZNmjQpcYJx7NgxzM3NadCggb7M39+/osIVQq9wm3dN0/Dw8MDLy6vUfXz00UdkZGQwYsSIMsfRoEEDBg8eTNOmTenYsSPr16/HysqKlStXlrlPUY4uHoPrKfyR90jljuvTMf/1718rd1whhLgPVPkEGcDe3p7evXuzYMECDh48SEZGBuHh4cYOS4g7Kpznm5KSwoIFC1i3bh0XLlwAwNXVtUR9/PLLL/z+++/Y2NhgbW3NK6+8wvXr17G2tubrr78uU1w2NjY0atSIxMTEMp0vylniDgB25TWs3HE9AsDSDhIlQRZCPHhMIkG+WbVq1fDw8ODatWsG5VZWVsWubuHr60tubi5HjhzRl8n8S1EZ6tati5OTE9HR0fqy7du3o9PpivwWIy0tjcTERNLS0gzKP/74Y+Lj44mNjSU2NpapU6dia2tLbGwsXbp0KVNcOTk5HDt2rMxTNkQ5+/s3sHPjhPKo3HHNLaF2K7mDLIR4IFkYO4B71aNHD3r06EHr1q3R6XREREQQHx9f5EE7X19ftm3bxsmTJ3Fzc8PS0hJzc3P8/f1p1aoVY8aM4f333+fIkSOsXbvWSFcjHiTm5uYMGTKE6dOn4+fnh7W1NW+//TYvvvgidnZ2Bm0jIyMJDQ0lPDycgQMH6ss9PT0N2nl4eKBpGg0b/nO3MSsrS79sW15eHidPniQ2NlY/H/+9997D29ubJk2akJ6ezrx587h27RovvPBCxV28KBml8hPUOm3hYiXOPy7k3Q5+eRdHrt29rRBCmJAqfwe5TZs2LFy4kLZt2xIYGMjGjRuJiIigffv2Bu3Gjx+Pi4sLjzzyCDY2NqxatUpft2bNGtLT02nevDkzZsyQDUVEpZk6dSp9+/bl6aefpmvXrnTs2JH58+eX6xhnzpyhWbNmNGvWjBs3bjB58mSaNWvGp59+CuRvVvLaa6/h7+9P586duXDhAtu2bcPDo5LvWIqiUv+Gq0ng3f7ubStCnfxxW5kdNs74QghhJFrhg0LGFBgYqGJiYowdhhBC3F/2rYENw2HYb3h/kEjirODKHT8nE2Z58XnGYwx+98vKHVsIISqBpml7lFKBt5ZX+SkWQghhStrN+oWk1PwtnudafEkXc3uaf3ACT0e7u5xZASx0UKslrU4cqvyxhRDCiEqUIGua5ggsBfwABQwCjgARgDeQCDyjlLpc0P4NYDCQC4xWSv1Q3oELIYQpSkq98c+d4g/fBLeOJPR7yngB1WlLo4Q5kHkNdA8ZLw4hhKhEJZ2D/CEQpZRqCAQAh4CJwBalVH1gS8FnNE1rBDwHNAa6Aos0TTMv78CFEMKkXT0DlxPzH9AzptqtMNcUJO0xbhxCCFGJ7poga5pWDQgCPgdQSmUppVKBnsCKgmYrgJCC9z2BL5VSmUqpBOAY0Kq8AxdCCJN2alf+q1cb48bh2SL/9fRu48YhhBCVqCR3kOsCKUC4pmn7NE1bqmmaHVBDKZUMUPDqVtDeEzh10/mnC8oMaJo2RNO0GE3TYkq6650QQjwwTu8Gcx24G3lnTxsnjuZ5wilJkIUQD46SJMgWQHPgE6VUM+A6BdMpbqO4xTqLLJWhlFqslApUSgWWdNewsti2bRuapul3KLuTKVOm4OfnV2GxCCFEiZ2Oyd/NzsLK2JGwN69+fsJ+H6x6JIQQlaEkCfJp4LRS6o+Cz5HkJ8znNE3zACh4PX9T+9o3nV8LOFM+4Rbv2LFj9O7dG0dHRx566CHatWvH3r17S93PuHHj+N///lekfMmSJbRs2ZJq1arh6upK3759ZRtecVdxcXH06tULDw8P7OzsaNGiBevXrzdo89dffxEcHIyTkxOOjo688MILXL58uVTjaJpW5Lh5N8jk5GT69euHm5sb9vb2PPXUU5w6dcqgj88//1y/RriPjw8zZswo+4WLe5eTBcmxUKulsSMBYJ/yhRuX4OJxY4cihBCV4q4JslLqLHBK07QGBUWdgYPARmBAQdkAYEPB+43Ac5qm6TRN8wHqA7vKNeqbJCcn07ZtWy5fvsymTZvYtWsXzz33HMnJyaXuy97enurVqxcp37FjBy+//DK//vorP//8M5cuXaJr167Fbl0tRKHY2Fi8vb356quviI+P57nnnqNv3776f4Tl5OTQo0cPzMzM+PXXX/nhhx+Ij49n0KBBpR5ryZIlJCcn64+bd9J78cUXSUhIICoqip07d5Kenk5ISAiFa6D/73//Y8iQIbz22mscOnSIefPmMW3aNJYvX14u34Mog3MHICcDahVZmtMo9ubVz38j85CFEA8KpdRdD6ApEAPEAd8ATkB18levOFrw6nxT+7eA4+QvBdftbv23aNFCldWoUaNUtWrVVFpaWrH1W7duVYCKiIhQDRo0UA899JAaMGCAys7O1reZOXOmIn8aiGrcuPFdx9yzZ48CVGxsbJnjFg8mf39/NWbMGKWUUgcPHlSA+vPPP/X1GzduVJqmqXPnzpW4T0CtW7eu2Lrr168rTdPU5s2b9WVxcXEKULt371ZKKTV37lzl7e1tcF7Hjh3VyJEjSxyDKD91JmxS6vfPlJpcTanLJ40djlJKKe8JG5WaUUupb8cYOxQhhChXQIwqJjct0TJvSqlYlT9fuIlSKkQpdVkpdVEp1VkpVb/g9dJN7acrpeoppRoopb6/9zT+9qKioujWrRt2dndeRH/p0qVERkayevVqVq9eTWRkpL5u1KhRJCcnExYWVqIxz507B4Cjo2PZAxcPHKUUKSkpODk5AZCVlQWATqfTt7G2tkYpxf79+0vV95gxY3B1daVNmzZ88803+vLs7GyUUkXGANi3bx8A7du358yZM/o72/Hx8cTHx9OtW7cyXKUoF0kxYO8ODrWMHQkACrP81SzkQT0hxAOipOsg37dOnjxJrVp3/5/IpEmT8PPzo0ePHrRq1Yrdu//5D72dnR3u7u7Y29vftZ/c3FxmzJhBnz59qFOnzj3FLh4sixcv5saNG7z88ssANGjQgBo1ajBv3jwyMjK4dOkS7733HpqmUZqVXWbMmEFkZCTff/897du3p1evXvz4448AODg4EBAQwIIFC7h69SrXr19n2rRpWFhY6Md49NFH+fLLL3nqqaewtLSkefPmTJ8+ne7du5f/lyBK5vTu/OkVWnHPPBtJrZZw/s/8DUOEEMLEVfkEGfIfUrobX1s5EYsAACAASURBVF9f/XtnZ2cuXbp0h9a3N2rUKC5evMiSJUvKdL54MG3fvp1x48bxxRdfULNmTSD/Tu7atWuJiorCzs4OT09POnXqlP+rHbOS/9V84403ePTRRwkMDGTu3Ln06NGDDz/8UF+/cuVKEhIScHR0xMnJCRcXF1xdXfVjHD58mNGjRzNz5kz27t3LsmXLePPNN9m4cWP5fgmiRJy4CpdO3DcP6OnVbgUqD87sM3YkQghR4ap8guzl5cXp06fv2s7CwnBXbVWG5YomTJhAVFQUP/30k0yvECW2a9cuevTowaefflrkruxjjz1GQkICycnJpKSk8PTTTwPok+iyaNGihcEqK02aNCEuLo4LFy5w/vx5pk+fTkpKin6MWbNm0bp1a0aMGIG/vz8vvvgioaGhzJ49u8wxiLJralawUsT9liDLhiFCiAdIlU+Qn3jiCb7//nvS09MrdJwpU6awdu1afvnlFzw9i+x7IkSxYmNj6datG/PmzeP555+/bbvCJdjWrl2Lvb09gYGGqxekpaWRmJhIWlraXcc8fPgw3t7eRcqdnZ1xdHRk3bp15ObmEhQUBEBqamqRO9YWFhbcuHGjBFcoylszs6OgmUPNpsYOxZCtMzj5wJlYY0cihBAVzuLuTe5vEydO1M+fnDZtGs7OzmzZsgVvb2+Cg4NL1MfZs2eB/CQkJydH/9nZ2RkrKytmzZrFhx9+yLfffou1tbW+3sHBARsbm4q5MFHlHThwgC5dujBixAiefPJJ/Z8bKysrnJ2dgfyHTG1tbalTpw5bt25l9uzZTJgwAVtbW4O+IiMjCQ0NJTw8nIEDB+rLN27cSFJSEu3bt0en07FhwwYiIiLYtGmTvs3OnTu5evUqDRs2ZN++fYwfP56BAwfqk+iuXbsyevRoVqxYQVBQEPHx8SxdupSRI0dW7BckihWgnQC3R8Dqzg8eG0XNZvkbmAghhImr8gmyp6cnO3bsYMKECXTr1o2cnBz8/f1ZtGhRifvw8PAo9vPWrVvp1KkTn376KampqXTo0MGg3a3JihA3i4yM5OLFi0ybNo1p06bpyzt27Mi2bdsAuHjxIq+88grnzp2jdu3aTJkyhXHjxpV4DEtLSxYuXMj48eMBeOSRR1i3bp3BChQ3btxg+PDhnDp1Cjc3N0JDQ5k6daq+/j//+Q/Xrl1j+vTpDB06FHd3d0aMGMH//d//3eM3IEpNKfzMEqBmT2NHUryaTeHP9XD9ItgVXTNeCCFMhVaWubjlLTAwUMXEyF0JIcQDLvUUzPeD7vOg1SvGjkbPe+JmEmcFQ8J2WPEUvPBf8O1i7LCEEOKeaZq2RylVZFemKj8HWQghTEZywfzems2MG8fteATkv8pKFkIIEycJshBC3C/OxJKjzKBGY2NHUjxrB3CuJw/qCSFMniTIQghxvzizj6OqFljexw//1mwmCbIQwuRV+Yf0hBDCJCgFybHE5/nxiLFjuYWnow3eEzcDMNhcx/9ZniZw4hfoHN35deLjRo5OCCHKnyTIQghxP7hyGtIvEq98eMbYsdzCIAlOdITla4gZ5IL3MlkrWwhhmqr8FAtN03BxcSEjI0Nf5ufnx5QpUwDo1KkTmqZhaWmJl5cXr7zyCmfOnDHoY+DAgWiaVuSIjIyszEsRosKNHj262D/bY8eOpXnz5lhaWtKpUyfjBPegK3hA70Cej5EDuQuPJoD2zwOFQghhgqp8ggxw+fJlIiIiblsfGhpKQkICy5cv58iRI7Rp04bLly8btOnSpQvJyckGx1NPPVXRoQtRaX766Sfi4uKKrcvKyuLll18u8eY6ogKciQXNnIOqjrEjuTPdQ+BSX1ayEEKYNJNIkIODg++4MYitrS21atXi8ccfZ+PGjVy4cIFPPvnEoI1Op8Pd3d3g0Ol0FR26EJXi0qVLDB8+nCVLlhRb//HHHzN8+HBq1apVyZEJveRYcHuETKyMHcndeTSVBFkIYdJMIkF+9tlnOXHiBHv27LlrW0dHR1q3bs2WLVsqITIh7g/Dhg3jP//5D/Xr1zd2KKI4SuXfQfZoauxISsajCVxLxpmrxo5ECCEqhEkkyDqdjkGDBvHxxx+XqL27u3uRechRUVHY29vrj3r16lVEqEJUulWrVpGYmMjYsWONHYq4natnIP3CPxtx3O/c/QF4xOxvIwcihBAVw2RWsRg6dCj+/v7Mmzfvrm01TePWLbaDgoJYvHix/rOFhcl8NeIBdurUKcLCwti6dSvm5ubGDkfczrkD+a/u/sAlo4ZSIjUKEmTtpJEDEUKIimEyWaCPjw8dO3YkPDz8rm2Tk5Px9PQ0KLO1tcXX17eiwhPCKPbs2cOFCxdo0aKFQXn//v354osvWL9+vZEiEwbOxue/1mgMRBs1lBKxqw4P1aRRqtxBFkKYJpNJkAFGjBjBq6++ipXV7R9yuXTpEr///jtvvfVWJUYmhHF06dKFgwcPGpQ98sgjzJs3j969exspKlHEuQPgWAesqxk7kpJz9+ORK4eMHYUQQlQIk0qQu3btSl5eHocPHzYoT09P5/Tp0/z1119MmjQJV1dXhg8fbqQohag89vb2NGzYsEh5zZo1DVasOHbsGGlpaVy4cIG0tDRiY2OxsrKiUaNGlRnug+tsvH5eb5Xh7o/vX1sgOwMsrY0djRBClCuTeEivkJmZGUOHDiUvL8+gPDw8HG9vb1566SUaNmzIzp07cXJyMlKUQtx/Xn75ZZo1a0ZERAR79uyhWbNmdO/e3dhhPRiyrsPF41UyQbbUciHl8N3bCiFEFVPl7yDf+rDd+PHjGT9+vP7ztm3b7trH8uXLyzkqIe5ft/6dgZL9PREV5PwhQEENP2NHUjoFD+px7gDUrCLL0wkhRAmZ1B1kIYSocgof0HOvYgmysw/Xle6f+IUQwoRIgiyEEMZ07gDoquU/pFeVmJlzWHlJgiyEMEmSIAshhDGdjc+fXqFpxo6k1A7lecHZA/k7AQohhAmRBFkIIYwlLw/O/Vn1plcUOKi8IfMKpMqGIUII0yIJcglMmTIFP7+q+T8wIcR9LDURstKq3gN6BQ7leeW/KdwJUAghTESVT5AHDhyIpmlomoZOp6NBgwbMnDmT3NxcY4cmxF1lZ2czcuRInJ2dqVatGqGhoVy/fr1UfURFRREYGIidnR2enp5MmTKlyEoV77zzDu7u7tjZ2RESEsL58+f1dTf/Hbr5CA4OLpdrFHdwtnCL6aqZIB9WtQFN5iELIUxOlU+QIX+3sOTkZI4cOcKIESN46623mDdvnrHDEuKuJk2axLp164iMjCQqKopt27YxZsyYEp9/7NgxQkJC6N69O3FxcSxatIj58+ezYMECfZvPPvuMOXPmsGjRIqKjo0lKSuL555/X13/44YckJyfrj7///hsHBweefvrpcr1WUYyz8aCZgVvV3JDlBtZQ3VcSZCGEyTGJBFmn0+Hu7o63tzejR4+mc+fOfPPNN/r6xMRENE1j06ZNdO/eHVtbW9zd3dmyZQsAJ0+e5Mknn8TOzg4nJydCQ0NJS0srMs7rr79OtWrV8PT0ZOnSpZV2fcI05ebmsnjxYt566y0ef/xx2rZty7vvvsuqVatKfBf5p59+wsHBgalTp1KvXj169uzJ4MGDWbJkib7NJ598wn/+8x969+5N8+bNmT9/Pj///DPHjh0DwMHBAXd3d/3x+++/k5ubyzPPPFMh1y1ucu4AVK8PljbGjqTs3P3gbJyxoxBCiHJlEgnyrWxsbMjOzi5SPm7cOEJCQoiLi2P58uXY2dkB0L9/f1JTU/ntt9/YuHEj27dvN9hsBODw4cMkJyeza9cuJkyYwJAhQzhwQObdibI7ceIEly5dokOHDvqyoKAgMjMziY8v2R25rKwsrKysDMqsra05dOgQGRkZ+r5uHqNNmzZYWloSExNTbJ9Lliyhb9++2Nvbl+GqRKmcPVBlp1fo1Wic/5Be5jVjRyKEEOXGpBLkvLw8vvvuO6KioujSpUuR+meeeYYhQ4bg6+tL165defTRRzlw4AC//vorH3/8MQEBAXTo0IF33nmHZcuWkZWVpT/XzMyMjz76iIYNGzJ69GgCAwNZtmxZZV6eMDEpKSkAuLi40LdvX1q3bo2Li4tB3d0EBQWRlJTE6tWrycvL49ChQ6xZs4a8vDwuXbrExYsXycvLw8XFhXHjxlG3bl0yMjJwcnIqdoyEhAS2bNnC4MGDy+9CRfFupMKVk1X2AT09t8b5r+dly2khhOkwiQQ5KioKe3t7rK2t6dWrFwMGDGDKlClF2t18F63QsWPHMDMzM1ilIiAggKysLE6e/GfpIk9PTxwdHfWfGzduzPHjx8v3QsQDpfBBOk3T8PDwwMvLq9R9NGvWjPnz5zNs2DCsrKzo2LEj/fv3B/L/UXfzGK6urnh5eWFhYVHsdtMAS5cu5eGHH6Zdu3ZlvCpRYucP5b9W+QT5kfzX8weNG4cQQpQjk0iQg4KCiI2N5fjx46Snp7NkyRKsra2LtHNycipSppRCu80C/bcrv9t5QpSEm5sbkH+3eMGCBaxbt44LFy4A4OrqWuJ+Ro8ezZUrVzh58iRJSUn4+PhgZWVF9erVcXFxwczMjJSUFCZMmMC2bduwsLAgNTW1yBg5OTmEh4cTGhpafhcpbq8woSxMMKsqxzpgaScJshDCpJhEgmxra4uvry+1a9fG3Ny8VOf6+vqSm5trMJ94//79WFpaUrt2bX1ZUlISqamp+s8HDx7E19f33oMXD6y6devi5OREdHS0vmz79u3odDr8/f0N2qalpZGYmFjsw6OQf7e4Zs2aWFpasnbtWtq1a4elpSU6nQ4/Pz+DMXbu3El2djaBgYEGfWzevJnz58/z0ksvleNVits6fwisHgKHWsaO5N6YmYFbQ0mQhRAmxcLYARibv78/bdq0YeTIkSxcuJCrV68yefJkQkNDDR5+ysvL49VXX+XNN9/kxx9/JCYmhuXLlxsvcFHlmZubM2TIEKZPn46fnx/W1ta8/fbbvPjii/oHSAtFRkYSGhpKeHg4AwcONKhbsWIFTZs2xc7OjqVLl7J9+3Z++uknff2wYcMICwujffv2eHt7M2bMGLp06VLkH3hLliyhW7dueHh4VNg1i3/sjfkNlevO0298Z1Du6Vh1VrTwdLTBe+JmZls8xOPme2k5cbO+/NeJjxs5OiGEKLsHPkEG+OKLLxg+fDht2rTBysqKnj17MnfuXIM2DRs2pEaNGgQGBlKtWjU+//xzGjWqmmuXivvH1KlTSUtL4+mnnyY3N5fevXszf/78UvWxf/9+xo4dS3p6On5+fmzYsIHHHntMXz906FDOnj3LsGHDuHbtGv/+97/57LPPDPo4ffo0UVFRrFu3rlyuS9yFUnjn/Y1z4NMk9qi6G7Lok+Cdf8MP20h8uxXYu+JdkCgLIURVpd3uYZ3KFBgYqG635JQQQpica+fgvYeh2xxo/R9jR3Pvjm+FVSHw0kao2xHviZtJnFV1E38hxIND07Q9SqnAW8tNYg6yEEJUKabygF6hGoVLvR0ybhxCCFFOJEEWQojKpk+QTWSalp0r2FaH838aOxIhhCgXkiALIURlO3+QFFUN7FyMHUn50LT8ZF/uIAshTIQkyEIIUdnOH+KvvNp3b1eVFCbI98FzLUIIca8kQRZCiMqUlwfnD/OXquLrH9/K7RHISoPUk3dvK4QQ9zlJkIUQojJdOQnZ1zmiTPAOMsg0CyGESajyCbKmabi4uJCRkaEv8/PzY8qUKcYL6j52+vRpnn32WZycnLC3t6dNmzZcvnxZX79r1y46duxItWrVcHV1ZeTIkWRmZurrs7KyCAsLo1atWtja2tK0aVM2bNhQ4vFzcnIICwujQYMG2NraUqdOHcaNG0d6erpBu5kzZ1K3bl1sbGxo2LAhS5cuNahfunQpnTp1wt7evtK2/P7000/x8fEhJyenXPstXFNbp9NRq1YtFi5cWKrzo6KiCAwMxM7ODk9PT6ZMmcLNyzcmJyfTr18/3NzcsLe356mnnuLUqVP6+h07dtCuXTuqV6+Ovb09rVu3JioqqtyuT9yiIIH8K88E7yCD7KgnhDAJJUqQNU1L1DQtXtO0WE3TYgrKnDVN+0nTtKMFr043tX9D07RjmqYd0TTtiYoKvtDly5eJiIio6GGqvIyMDDp37kxKSgqbN28mLi6ON954AwuL/P1iUlNTCQ4Opn79+uzZs4eIiAg2btzIxIkT9X3MmjWLVatWsWLFCv78809CQkLo06cPR48eLVEMmZmZxMfHM2PGDPbv38+yZcv46quvGDVqlL7NypUrmTp1KvPmzePQoUOMHTuWIUOGsHXrVn2ba9eu0bVrV0aPHl1O386d5ebmMm/ePMLCwvTfV3lYtmwZo0aNYsSIERw4cIBvv/2WJk2alPj8Y8eOERISQvfu3YmLi2PRokXMnz+fBQsW6Nu8+OKLJCQkEBUVxc6dO0lPTyckJESfRFtbWzNy5Ei2b99OfHw8ISEh9OzZk4MHJdGpEAUJpMlNsbCuBg61JUEWQpgGpdRdDyARcLmlbA4wseD9RGB2wftGwH5AB/gAxwHzO/XfokULVVaAeuqpp1SrVq30ZY0bN1aTJ08uc5+mKjw8XDk6OqqrV68WW//dd98pTdNUWlqavmzBggXKwcFBZWdnK6WUCg4OVgMHDtTX5+TkKHNzcxUZGVnmuN577z3l6Oio/zxixAjVqVMngzZ16tRR8+bNK3LuunXrVP4f47IDVEJCwh3bfPnll8rFxUWlp6ff01i38vb2VlOnTi3z+YsWLVJubm4GZa+99ppq3LixUkqp69evK03T1ObNm/X1cXFxClC7d+++bb/Ozs7qs88+K3Nc4g4iByv1vp+qM2GTsSMpf6v7KLWorWlemxDCJAExqpjc9F6mWPQEVhS8XwGE3FT+pVIqUymVABwDWt3DOHf17LPPcuLECfbs2VORw1R527Zto3379kyePBl3d3f8/f0Nfp2flZWFpmlYWVnpy6ytrbly5QoJCQkABAUFER0dzcmTJ1FK8eWXX/LQQw/Rtm3bMsd17tw5nJz0v4AgKCiI/fv38+ef+Wuqbt26lQsXLtClS5cyj3GvZs+ezejRo7GxsSm3PhMTE0lMTNRvYe7h4UHPnj1JTEwscR9ZWVkGPy/I/5kdOnSIjIwMsrOzUUqh0+kM6gH27dtXpL/c3FzWrFlDamoqzZo1K9uFiTs7d9B0Ngi5lVsjuPAXFpTvNCQhhKhsJU2QFfCjpml7NE0bUlBWQymVDFDw6lZQ7gmcuunc0wVlFUan0zFo0CA+/vjjihymyktOTmb79u2cPXuWzZs3M2bMGF577TUiIyMBaNWqFTqdjrlz55KTk8Pp06f57LPPAEhJSQFg/Pjx9OvXD29vb6ysrBg7dizfffcdHh4eZYrpzJkzLF26lLCwMH3ZM888w/Tp02nevDmWlpY89dRTrFmzhoCAgHv8Bsrmhx9+4OjRo4wcObJc+01OTgZgxowZTJgwgY0bN5KWlkbPnj3Jy8srUR9BQUEkJSWxevVq8vLyOHToEGvWrCEvL49Lly7h4OBAQEAACxYs4OrVq1y/fp1p06ZhYWGh/5kWqlWrFjqdjmHDhrF+/XpatmxZrtcrgNxsuPCXaSfIuVl4a2eNHYkQQtyTkibI7ZRSzYFuwAhN04Lu0La4J6aKLIypadoQTdNiNE2LufV/1GUxdOhQvvrqKy5dunTPfZmqvLw88vLyWLp0KS1atGDw4MH06tWL1atXA+Dh4cHy5cuZP38+1tbWNGrUiJ49ewJgZpb/RyUyMpI1a9bw9ddfExMTw9ChQwkJCSnVXc9CaWlphISE0LVrV0aMGKEvj46OZvr06Xz++efs3buXmTNnMmDAAHbv3n3vX0IBe3t7/QHQuHFj/efo6GiDtrNnz2bIkCEGd7kLzZgxw6CvkydLvsRVYRI8ePBg+vbtS8uWLVm4cCFxcXHExcWVqI9mzZoxf/58hg0bhpWVFR07dqR///7APz+zlStXkpCQgKOjI05OTri4uODq6qqvLxQdHc3u3bsZNWoUw4YNK/G8cnFn7Wb9gvfEzXhP3EyXt5dBXjZjtmbh6Vh+v424b9TIX8migXbayIEIIcQ9Km7exZ0OYAowDjgCeBSUeQBHCt6/AbxxU/sfgDZ36vNe5yCvW7dOKaVU9+7d1bx582QO8m0899xzqmHDhgZlEydOVIGBgUXanjlzRt24cUP9+OOPClB///23Uip/LvBHH31k0LZ58+bqzTffLFUs6enpqlOnTio4OFhlZWUZ1HXs2FGFhYUZlPXu3Vv179+/SD9lnYN89OhR/QGobdu26T/fPM94165dysrKSp06darYfi5evGjQV+Fc7ZI4fPiwAtSXX36pL8vIyFCA2rSpdHM4c3NzVVJSksrKylKLFy9WVlZWRb7XixcvqsuXL6v09HRlYWGhVqxYcdv+Hn/8cTVs2LBSxSCKZzAfN/6/Sk2uptSZ/cYLqCJl3VBqipOa/9YAY0cihBAlwm3mIN/1cXxN0+wAM6XUtYL3/wamAhuBAcCsgtfCtb42Al9omvY+UBOoD+y6pyy+hEaMGMGrr75aZE6myNesWTM2b95MZmamfk7qyZMn8fQsOgOmcMrE2rVr8fHxwcvLC8hf6eLWO48WFhbcuHHDoKxwioa9vT0uLobb6WZmZtKrVy8sLCyIjIzE0tLSoL6kY9wLX19fg8916tTB29u7SLuZM2fy/PPPU6tW8SsOODs74+zsXKYY6tWrR7Vq1Thx4oS+rPAO9K0/k7S0NC5cuICLi4v+rvfNzMzMqFmzJpD/M2vXrl2R77UwzpUrV5Kbm0tQ0O1/EWRubs7169fLdF3iDs4fAs0MXB42diQVw9IaqtejwTm5gyyEqNpKsl5VDeDrgrVmLYAvlFJRmqbtBr7SNG0wcBLoC6CU+lPTtK+Ag0AOMEIplVsh0d+ia9eu5OXlcfjw4coYrsrp378/77zzDqNGjWLcuHHs37+f9evXs3LlSn2byMhIvLy8cHNzY/369axYsYLw8HB9fdeuXZkzZw4NGjTAx8eHDRs2sHv3bmbOnGkw1unTp/Hx8WHAgAEsX75cX56dnU2fPn04d+4c69evJzU1VV/n6uqKubk5Xbt2ZfHixbRt25aAgAB+/fVXvvnmGz799FN927Nnz3L27Fn91I7Y2FgAGjVqVG7/QDpy5AgbN27kwIED5dLfrSwsLBg4cCDz58+nefPm1KlTh/Hjx+Pv74+/v79B28jISEJDQwkPD2fgwIEGdStWrKBp06bY2dmxdOlStm/fzk8//aSv37lzJ1evXqVhw4bs27eP8ePHM3DgQP0/CN5//31q1qxJQEAA5ubmrF+/np9//pn169dXyHU/0M4fBOd6+YmkqXJtgO95eWBaCFHFFXdbubKP8ppioZRSc+bMUYBMsbiNLVu2qKZNmyqdTqfq1aun5s+fb1A/d+5c5ebmpiwtLVWjRo1UeHi4Qf3ly5fV0KFDVc2aNZWNjY3y9/dXa9asKTJOQkKCAtSAAQOKLS/uKFxqLSMjQ02cOFHVqVNHWVtbq4cffrhInJMnT75jH+Vh0KBBKiQkpNz6K86NGzfUsGHDlLOzs3JwcFDdunVTx48fL9IuPDxcAUV+HkopNXbsWOXk5KR0Op1q0aJFkekZW7ZsUXXr1lWWlpbK09NTTZgwQWVmZurrP/roI+Xn56fs7e2VnZ2datq0qVq5cmW5X+uDymCKxYLmSn35vPGCqQxbpqnsSY5KZWfeva0QQhgZt5lioSlV5Pm5ShcYGKhiYmKMHYYQeklJSdStW5ft27fTunVrY4cjqjDviZtJnBUMOZkw3R06hMHjbxs7rIoTHwn/HQzDfzfd1TqEECZD07Q9SqnAW8ur/FbTQlSE06dPM3fuXEmORfm5eBxUHrg0MHYkFcu14PpSZKqbEKLqKr89c4UwIa1bt5bkWJSvwoTR1cQT5Oq+5CoN85Qjxo5ECCHKTO4gCyFEZbjwF6CBS31jR1KxLG04pdzkDrIQokqTBFkIISpDyhFwqgOWJrhByC2OKs/86xVCiCpKEmQhhKgMKUdMf/5xgWPKEy4chdwcY4cihBBlIgmyKJVPPvmEJk2aYG9vT/Xq1enRowdHjpTuTlGnTp3QNM3gWLhwob4+KyuLN998kzp16mBtbU2LFi3Yvn27vv7ixYt07dqVmjVrYm1tzcMPP8wHH3xQbtd4O59++ik+Pj7k5JTf//Szs7MZOXIkzs7OVKtWjdDQ0FJv0JGUlERwcDC2trbUrFmT2bNnF2nzzjvv4O7ujp2dHSEhIZw/f77Yvnbv3o2lpSVPPvlkma5H3EZuDlw8Bq4mukHILY7meUJeNlxOMHYoQghRJpIgi1Jxd3dn9uzZxMbGEh0djbW1NU888QS5uaXbCyY0NJTk5GT9MWjQIH3dnDlzWLJkCZ999hkHDhygY8eOdO/enTNnzgCgaRpPPvkkmzZt4siRI8yePZvJkyezdOnScr3Wm+Xm5jJv3jzCwsKwsCi/Z1snTZrEunXriIyMJCoqim3btjFmzJhS9dGnTx+uXr3Kb7/9xnvvvcekSZNYs2aNvv6zzz5jzpw5LFq0iOjoaJKSknj++eeL9JOens7gwYNp3rz5PV+XuEXq35CbCa4NjR1JpTiqCnaelHnIQoiqqrjFkSv7KOtGIVu3blUeHh6qR48eyt3dXX300UeqZcuWys3NTf3888/6NoBKSUnRnxccHFxkAwtRNnFxcQpQR44cKfE5HTt2VCNGjLhtfatWrdT48eP1n3Nzc5WLi4uaO3fubc/p3bu36tevX4ljuBkl2GDkyy+/VC4uLio9Pb1MYxQnJydHOTs7qw8//FBftnr1aqXT6VRaWlqJ+ti3b58C1P79+/VlL7/8smrXrp3+c0BAgBo7P1acywAAIABJREFUdqz+844dOxSgjh49atDX0KFD1bvvvqsGDBiggoODy3pZ4hZ1JmxS6tBmpSZXU+rkLmOHUykemRCZf73/m2PsUIQQ4o64zUYhVf4O8tmzZwkLC6Nfv36MHj2aOXPm8NJLL/Hee+8ZOzSTd/36dZYsWUKNGjWoVatWqc6NiIjAxcUFf39/Zs6caTBtISsrC51Op/9sZmaGlZUV+/btK7avPXv2EB0dTYsWLcp2ISUwe/ZsRo8ejY1N+T1gdeLECS5dukSHDh30ZUFBQWRmZhIfH1+iPnbv3o2DgwNNmjQx6GPv3r0opfR93TxGmzZtsLS05ObNeb777jv++OMPJkyYUA5XJoq4UDAN6QGZYpGONTh4yYN6Qogqq8onyK6urgQFBfGvf/0LV1dXOnXqROfOnUlMTDR2aCYrPj4e+/9v787Do6wO9o9/TxKyESCETEgIS4gJUFbRKFhUKqLFF31BW6tdKC59UcRarVZjS6X2VxWt1n2timLdl1YrVitSBKkL4AayCAKShC0BAwQSSDLn98dMYgIBsszkzHJ/rivXzDx55sk9A4Q7T85zTkoKnTp14q233uLdd98lOTm52c+fNGkSzz//PPPmzeOyyy7jlltu4YYbbqj//OjRo3nmmWdYt24dtbW1PPDAA2zbto3S0tJGx/nxj39MQkICxx13HFOnTuXXv/51wF5jQ2+99RZr1qzh8ssvD+hx615Peno65557LiNGjCA9Pb3R55pzjG7durF371769u3LtddeS3p6OpWVlVRUVLB9+3a8Xi/p6elcc8015ObmUlVVRdeuXeu/RllZGZdccgmPPfZYQIePSAOlq6FTFiR2cZ2k/Xj6a4iFiIStsC/IdWf0kpKS6u8nJiZSWVnpMlZE69+/f/0Y5P79+3PRRRexf//+Zj//4osvZsyYMQwdOpSpU6dy4403cs8992D9y57/4Q9/YNCgQeTl5ZGQkMDLL7/MuHHjiIlp/Nf1zjvv5JNPPuGvf/0rDzzwAH//+9+bnSElJaX+A2DQoEH1jxcuXNho31tvvZUpU6bQtWvXg45z8803NzrWxo0bm52h7vUaY8jKyqJ3797Nfm7DYxhjiI2NpVevXmRkZNQf98Cv4fF46N27N3FxcY32mTJlCpMnT2b48OEt/vrSTKWrIT06zh7X8/T3zWThbdn1CSIioSBiTxc1LAYH8nq97R0nosTHx5OXl0deXh4vvvgiaWlp/OMf/+BHP/pRq4537LHHsmfPHsrKyvB4PKSmpvLaa6+xd+9edu3aRWZmJqNGjaJ//8ZTZGVmZpKZmcnAgQPZtGkT06dP55xzzmnW1/z000/r7+fn5/PGG2+QnZ0NUH8LviEMixYtYvbs2U0e59JLL230unv06NHs152RkQH4zgLfc889ABQVFQG+34w09xhlZWUkJCTUz/Tx1FNPkZSUREpKCvHx8cTExFBaWsp1113Hddddh9frpby8vP5rzJs3jzlz5nD77bcDvpk1wPeDZklJCd26dWv2a5KmWN8iIUcffGFkRPMMgJoqKN8IaX1dpxERaZGILch1UlNTAaioqKj/9XVRUVF9OZG2qTure+DUZDU1NRQXF5OSklL/vh/KqlWr6Nix40H7JScnk5yczNdff80HH3zAJZdccshjxMbGtmh6tLy8vEaP+/TpQ05OzkH73XLLLfz0pz895BjrtLQ00tLSmv11G8rNzaVr164sXLiw/uztggULSEhIYMiQIY32raiooKysjPT09Pqz3gAFBQXs3LmTzz//vH4c8oIFCzjmmGMwxpCQkMDgwYNZuHAhZ599NgDvv/8+1dXVFBQUAL4fAhrOQnL99ddTXl7Ogw8+WP/vR1ovix2wvyJqxh/Xq5uxo3S1CrKIhJ2IL8j5+fl06tSJF154gWuvvZZ//vOfrFy5MqgXdEWyK6+8kjPOOIP8/Hx27drFbbfdRmxsLKeeemqj/YqLi+nbty+TJ0/miSeeqN++du1aZs+ezZlnnonH42Hx4sX8/ve/5/LLL68/279t2zbefPNNRo0aRWlpKVdccQX9+vXjvPPOA+C1117j66+/5sQTTyQ1NZXFixdzxx138POf/zygr3X16tW89tprLF++PKDHrRMbG8uUKVO46aabGDx4MImJiUyfPp1JkybRsWPHRvu+9NJLXHjhhcyaNYsLLrigfvvw4cMZOXIk06ZN47777mPlypXMnj2bxx9/vH6fqVOncvXVV3PiiSeSk5PDlVdeydixY+t/SMjPb7z0cZcuXaiurmbAgOiYkizY8mJKfHeiZIq3enU/EJSugv7j3GYREWmhiC/IycnJPPTQQ/zmN7/hzjvvZMKECYwdO9Z1rLBVUVHBpZdeyubNm+nYsSPHHHMM//73v5s9fjYhIYG5c+dy7733snfv3vrCds011zTa7+677+aSSy4hMTGRcePGcccdd9TPbNGpUyeee+45ZsyYwZ49e+jVqxeXXXYZ06dPb9Vrajget6HbbruNs846K6hF8Y9//CMVFRX84Ac/oLa2lnPOOYe77rqrRcd48cUXmTJlCiNHjiQ1NZUbb7yx0TzHl156KVu2bGHq1Kns3r2b008/nYcffjjQL0UOId/4C3KUrKJXL7GL78JEzWQhImHIHKoctKeCggLbcMopEddKSkrIzc1lwYIFjBgxwnUcCWPPTD+bn3T6BK5dD01cExGJcgrnsGHmeJg9Aap2wZT/uI4kItIkY8xSa23BgdvDfhYLkWAoLi7mz3/+s8qxtFleTIlveEWUlONGPAN8Z5BD4ESMiEhLRPwQC5HWGDFihMqxBESeKYH0ka5juOHpD9V7YGcxpPZynUZEpNl0BllEJFj2lJFmKqLvAr06DWeyEBEJIyrIIiLBUreSXLRN8VanviBrRT0RCS9RVZBzcnLqF0MQEQm6ujOn0TaDRZ3kNOjoUUEWkbATVQV58eLFXHbZZa5jhL2HHnqIPn36kJSUxCmnnMLatWtbfawJEyZgjKHhLCZVVVVcdNFFDBo0iNjY2Ebz/h5o586d9O7du9HiGcHy0EMP0bdvX2pqagJ2zOrqai6//HLS0tLo3LkzF154YYsWPAHfjBvjx48nOTmZHj16cOuttx60z4033khmZiYdO3Zk4sSJbNu2rf5z7733HqNGjaJbt26kpKQwYsQI3nzzzTa/NgFKV1NhE6FL0wvNRIW6C/VERMJIVBVkj8dDcnKy6xhh7c033+Tyyy/nd7/7HR999BHJycmcddZZjVZia65HH32U8vLyg7bX1tYSFxfHVVddxfHHH3/YY0ybNo2+fYO/SldtbS233347V199NXFxgbu29YYbbuDFF1/kpZde4s0332T+/PlceeWVLTrGD3/4Q3bt2sV///tf7rjjDm644Qaefvrp+s8//PDD3HbbbTzwwAMsXLiQkpKSRvMkJyYmcvnll7NgwQKWLVvGxIkTmTBhAitWrAjY64xaZav5yvaIzhks6nj6ayYLEQk/1lrnH8cee6xtixtuuMH27NnTJiQk2KOOOsrec889jT7fv39/C1jA/vnPfz7o+evWrbMnnXSSTUhIsAUFBXb69Om2Y8eO1lpr169fbwH7i1/8wqanp9uZM2fa0047zaalpdnZs2dba62tqamxF198se3bt6+Nj4+3vXr1sjNmzLC1tbVtel2haMKECfbss8+uf1xcXGwB+/bbb7foOGvXrrVHHXWUXbp0qQXs4sWLm9xv/PjxdvLkyU1+7rnnnrNjx461s2bNqv/zag3Arl+//rD7PPfcczY9Pd3u3bu31V/nQDU1NTYtLc3efffd9dv+9re/2YSEBFtRUdGsY3zyyScWsJ999ln9tl/84hd21KhR9Y+HDRtmr7rqqvrH7733ngXsmjVrDnnctLQ0+/DDD7fk5UhTbh9gX5o+3nWKdtfnute/ffDhI9bO6Gztzk3uAomIHAKwxDbRTcP+DPIrr7zCbbfdxoMPPsiqVat45JFH6NSpU6N93nvvPTZv3kzPnk3/mnPSpEkYY1i8eDHTp0/nvvvuO2ifMWPGcP3111NYWMjkyZOZMWMGt9xyC+A7uxgbG8usWbNYtWoVDz74IHfddVdErla2ePFiTjrppPrH2dnZ5Obm0pKFXmpra5k0aRIzZ84kLS2tVTlKSkr4zW9+w1//+tdWPb+lbr31Vq644gqSkpICdsx169axY8eORu/nySefzL59+1i2bFmzjrF48WK6dOnC0KFDGx3j448/xlpbf6yGX+OEE06gQ4cOTf6Z1dbW8vTTT1NeXs7w4cPb8OqEqp2wexNrvdmuk7jl8Y+/1jhkEQkjYV+Q169fT2pqKuPGjSMnJ4cxY8YcNGY1PT2dzMxMYmNjD3r+F198waJFi7j33nsZMmQIEyZM4Pzzzz9ov4kTJ3LGGWcAcM4553D66aezYcMGAOLj43n44YcZPXo0ffv2Zfz48UyYMIF//etfAX+9rpWWlpKens4999xDRkYGRUVFpKenU1pa2uxj3HzzzfTo0YMf/vCHrcpgreWCCy7g2muvJScnp1XHaIm33nqLNWvWcPnllwf0uHXvWXp6Oueeey4jRowgPT290eeac4xu3bqxd+9e+vbty7XXXkt6ejqVlZVUVFSwfft2vF4v6enpXHPNNeTm5lJVVUXXrl0P+ho9e/YkISGBqVOn8sorr3DccccF9PVGnbI1AKyxUV6Q01WQRST8hH1BPvvsswHIz8/n//7v/3jqqafYv39/s5+/Zs0aYmJiGDRoUP22hvfrJCUl1Z89TEpKIjExkcrKyvrPP/TQQxQUFODxeEhJSeHZZ5+loqKitS8rpBljSEtLo0+fPiQkJGBbMLbw448/5v777+f+++9v9de/9957qaqqYtq0aa0+RkpKSv0H+P7M6x4vXLiw0b633norU6ZMoWvXrgcd5+abb250rI0bNzY7Q937ZowhKyuL3r17t/h1WGsxxhAbG0uvXr3IyMho9OfR8Gt4PB569+5NXFxck39mCxcuZPHixfzyl79k6tSprFmzpsV5pAF/IVxrezgO4lhKBiSm6kI9EQkrYb+SXm5uLl999RXvvPMO8+fP54orruCZZ55p9tnbuoLRFs8//zy/+tWvuOOOOxg9ejRJSUn87ne/Y+vWrW06bijyeDyUlpZy1VVX8bOf/QyA7du34/F4mvX8BQsWUFpaSp8+fYBvC9yoUaO47LLLuPPOO494jHnz5vHBBx/U/8BSW1tLTU0NiYmJPPvss/U/NB3Op59+Wn8/Pz+fN954g+xs35m+ulvwDWFYtGgRs2fPbvI4l156KT/60Y/qH/fo0fwylJGRAfjOAt9zzz0AFBUVATT7/czIyKCsrIyEhAQWLFgAwFNPPUVSUhIpKSnEx8cTExNDaWkp1113Hddddx1er5fy8vKDvkbdxY7Dhw/ngw8+4M477+SBBx5o9uuRA5Suhth4imyG6yRuGeMbZlH2peskIiLNFvYFGaifSeGss85ixIgRnHfeeVRVVZGYmHjE5/br14/a2lqWL1/OsGHDAFi+fHmLvv57773HCSec0OhX8OvWraNjx44teyFhoKCggIULF3LVVVcBvrHA69ato6CgoNF+NTU1FBcXk5KSUj9sAGDy5MmMGzeu/nFJSQljx47lmWee4cQTT2xWhvvvv5+ZM2fWP/773//On/70J5YuXdqo3B5OXl5eo8d9+vRpcrjGLbfcwk9/+tNDjl9PS0tr9Tjq3NxcunbtysKFC+vH+y5YsICEhASGDBnSaN+KigrKyspIT09vNKVdQUEBO3fu5PPPP68fh7xgwQKOOeYYjDEkJCQwePBgFi5cWP+Dw/vvv091dfVBf2YNxcbGtni6OTlA2ZfQLZ/aPQcP7Yo6nv6w6g3XKUREmi3sC/Ls2bPZv38/3/3ud4mJieH5558nPz+/vhxXVlayc+dOwHemcffu3WzZsoXY2Fg8Hg+DBg1i1KhRXHHFFdx3332sW7eOl19+uUUZ+vXrx+zZs5k7dy69e/fm8ccfZ9WqVRx77LEBf72uTZ06lTPPPJO//vWvjBw5ksLCQgYMGMApp5zSaL/i4mL69u3L5MmTeeKJJ+q3d+3atdFQhbo/pz59+tC9e/f67StWrGD//v3s2rWLmJgYPv30U1JSUsjLyzuoBGdlZWGMYcCAwC7nu3r1al577bUW/8DUXLGxsUyZMoWbbrqJwYMHk5iYyPTp05k0adJBP1y99NJLXHjhhcyaNavRGPvhw4czcuRIpk2bxn333cfKlSuZPXs2jz/+eP0+U6dO5eqrr+bEE08kJyeHK6+8krFjx9b/kPCXv/yFHj16MGzYMGJjY3nllVeYO3cur7zySlBed9QoXQU9hkPzR91EjOzUJHIK59Q/vji2lt93KON/bvkHb1w/0WEyEZFmampqi/b+aMs0b6+++qodMWKE7dSpk+3cubMdO3asXbZsWf3nZ82aVT/FW8OPPn361O9TN81bfHy8Pe644+zvf/97m5aWZq39dpq3w93ft2+fvfjii21qaqpNTU21v/zlL+0ll1xiR48e3erXFcoeeOAB26tXL5uQkGBHjx5tv/zyy4P2qXt/DjVF24H7HTjNW58+fQ76MzvU+9nWad4O5aKLLrITJ04M+HEb2rdvn502bZpNTU21nTp1spMnT25yire6v8ezZs066HNFRUX2jDPOsImJiTYzM9PecsstB+0zY8YMm5GRYZOSkuyECRPsli1b6j9377332sGDB9uUlBTbsWNHe/TRR9dPYSittH+vtTO6WPufWxpPeRatvnzb2hmd7Q8Lb3edRESkEQ4xzZuxITB5e0FBgW3JNGHBNmPGDF5//XWWLl3qOoo4UlJSQm5uLgsWLGDEiBGu40i42fw5PHwS/HAWOX9LYMPM8a4TuVW+Ee4awm+rL+bmm/7iOo2ISD1jzFJr7UFjDsN+iEUgvPrqq1RXVzN8+HC++uorHn74YQoLC13HEoeKi4v585//rHIsrVN3QZqnP7DBZZLQ0LkndOhIXk2J6yQiIs2iggxUVVXxu9/9jqKiIrKysrjkkksCPuethJcRI0aoHEvrla4GEwPd8lBBBmJiID2fvGIVZBEJDyrIwHnnncd5553nOoaIRIrSVdC1L8QluE4SOjwDyNv0b9cpRESaJewXChERCTllX4InsLOqhD1PP3qYHVC1y3USEZEjUkEWEQmk2mrYvhY8/VwnCS11PzCUaYVGEQl9KsgiIoG0Yz14ayC9v+skoaXu/fAvwS0iEspUkEVEAqlste/Wo4LcSNcc9tm4b98fEZEQpov0REQCYNTMeZSUVzIt9jV+0wEG3ruevWwmOzXJdbTQEBvHepvFgFIVZBEJfSrIIiIBUFJe6VsQ5OV/wMZerPjDD1xHCjlrbQ8VZBEJCxpiISISSKWrIF0X6DVlrc2GbzZAdaXrKCIih6WCLCISKF6vb5YGjT9u0lpvNmB9s3yIiIQwFWQRkUDZWQQ1lSrIh7DGZvvuaJiFiIQ4FWQRkUCpK36a4q1J622WbwluFWQRCXEqyCIigaIp3g5rPx18S3BrqjcRCXHNLsjGmFhjzCfGmNf9j9OMMW8bY9b4b7s22Pd6Y8xaY8xqY8z3gxFcRCTklK6Gjh5ITnOdJHR5+usMsoiEvJacQf4VsLLB40LgHWttPvCO/zHGmIHA+cAgYBzwgDEmNjBxRURCWOlqDa84Ek9/2P6Vb0luEZEQ1ayCbIzpCYwHHm2weQLwpP/+k8DEBtufs9bus9auB9YCxwcmrohIqLK+oQMaXnF46f3BW+1bkltEJEQ19wzyXcC1gLfBtu7W2s0A/tsM//ZsoKjBfsX+bY0YY6YYY5YYY5aUlpa2OLiISCjxUA5VO1WQj6Tu/dE4ZBEJYUcsyMaYM4Ft1tqlzTymaWKbPWiDtY9YawustQUej6eZhxYRCU15MZt8d1SQD69uEZXSVW5ziIgcRnOWmh4F/K8x5n+ARKCzMeZvwFZjTJa1drMxJgvY5t+/GOjV4Pk9gU2BDC0iEmryTInvjsYgH15CCnTpBaVfuk4iInJIRzyDbK293lrb01qbg+/iu3nW2p8BrwGT/btNBl71338NON8Yk2CM6QvkAx8FPLmISAjJNyWQ0Bk6ZbqOErKyU5PIKZzDuzu6svyzj8gpnENO4RxGzZznOpqISCPNOYN8KDOBF4wxFwMbgXMBrLVfGGNeAFYANcA0a21tm5OKiISwPFPiG15hmhplJgCLCsf47ry5CJY8zoY/nAExMeQUznEbTETkAC0qyNba+cB8//3twKmH2O8m4KY2ZhMRCRv5MSWQfpbrGOHB08+3JPfOjdA1x3UaEZGDaCU9EZG22rsDj9npK35yZHXjtDUOWURClAqyiEhblfmLnmeA2xzhom6mD81kISIhSgVZRKSt6pZO1hRvzZOc5luSW3Mhi0iIUkEWEWmr0tVU2njo0tt1kvDhGfDtDxYiIiFGBVlEpK3KVvOV7QEx+pbabOn9fGOQ7UHrSImIOKfv5iIibVW6mrW2h+sU4cUzAPbthIqtrpOIiBxEBVlEpC32VcDOItZ4e7pOEl48WnJaREKXCrKISFv4Z7DQGeQWqpvxQ1O9iUgIUkEWEWmL+oKc7ThImEnpDglddAZZREKSCrKISFuUroKYOL623V0nCS/G+KbFK9MZZBEJPSrIIiJtUfoldMujhjjXScKPp5/OIItISFJBFhFpi9JVWiCktTwDYE8pqex2nUREpBEVZBGR1qqugm/WQ7oKcqv437c8U+I4iIhIYyrIIiKtteMrsF6dQW4t//uWF7PJcRARkcZUkEVEWqtu/KwKcut06QUdksnXGWQRCTEqyCIirVX6JZgY6JbnOkl4ivG9dxpiISKhRgVZRKS1SldBah/okOQ6SfjyDCAvRgVZREKLCrKISGuVffntinDSOp5+ZJvtsE8zWYhI6FBBFhFpjdoaKFvjm8tXWq/uBwwtGCIiIUQFWUSkNb7ZAN5qnUFuq7op8kpVkEUkdGjpJxGRFhg1cx4l5ZWcHrOYR+JhwvOlfPbcHLJTNQ65VdL6st/GEl+22nUSEZF6KsgiIi1QUl7JhpnjYeEqeAde/cNFkNDJdazwFduBDTaTfqUqyCISOjTEQkSkNUpXQ+dsleMAWGOzfe+niEiIUEEWEWmN0tVaICRA1tps35Ld1VWuo4iIACrIIiIt5/X6Zl1IV0EOhLXebN+S3Tu+ch1FRARQQRYRabmdRVC9V2eQA2StzfbdqVu6W0TEMRVkEZGWqpuzVwU5INbZLN+S3ZrqTURChAqyiEhL1Z3p1BzIAbGPeN+S3TqDLCIhQgVZRKSlSldDcjokp7lOEjk8A7SanoiEDBVkEZGWKl2ts8eB5unnW7q7tsZ1EhERFWQRkZaxUKYp3gIuvb9v6e5vNrhOIiKigiwi0hIeyqFqpwpyoNWdkdc4ZBEJASrIIiItkBezyXdHBTmw0vN9t2VaUU9E3FNBFhFpgXxT7LujRUICK7Gzb+luLTktIiFABVlEpAX6mWJI7AKdMl1HiTzp/VSQRSQkqCCLiLRAfkwJeL4DxriOEnnqpnrzel0nEZEoF+c6gIhI2LCW/qYIMk5wnSSiZKcmkVM4h5/E7ufmDnsZ9dvZlOAhOzWJRYVjXMcTkSikgiwi0lwVW0k1eyDjO66TRJT6Evx1V5j1GIsu7gn5Y8kpnOM2mIhELQ2xEBFprm0rfLcqyMFRd+GjpnoTEcdUkEVEmmvbSt+tRwU5KDp28y3hraneRMQxFWQRkebatpIy2xlSPK6TRC5Pf81kISLOqSCLiDTXtpWs8fZ0nSKy1RVka10nEZEopoIsItIc1kLpKlZbFeSgSu8PVeVQsc11EhGJYirIIiLNsbMI9lfwpe3lOklkyxjgu627IFJExAEVZBGR5tjmm1nhS2+24yARLmOQ77bugkgREQdUkEVEmsN/RvNLDbEIrhSPbyaLbV+4TiIiUeyIBdkYk2iM+cgY85kx5gtjzI3+7WnGmLeNMWv8t10bPOd6Y8xaY8xqY8z3g/kCRETaRekq6JTFLlJcJ4l83QfqDLKIONWcM8j7gDHW2mHA0cA4Y8xIoBB4x1qbD7zjf4wxZiBwPjAIGAc8YIyJDUZ4EZF2s22FFghpLxkDYdsqDF7XSUQkSh2xIFufCv/DDv4PC0wAnvRvfxKY6L8/AXjOWrvPWrseWAscH9DUIiLtyVsLpV9qgZD2kvEdqN5DT1PqOomIRKlmjUE2xsQaYz4FtgFvW2s/BLpbazcD+G8z/LtnA0UNnl7s3yYiEp6+2QA1lTqD3F78F+r1N8WOg4hItGpWQbbW1lprjwZ6AscbYwYfZnfT1CEO2smYKcaYJcaYJaWlOksgIiGs1DeDhQpyO/FP9dbfFB1hRxGR4GjRLBbW2nJgPr6xxVuNMVkA/tu6Wd2LgYYThfYENjVxrEestQXW2gKPR8u2ikgIq5uT19PfbY5okdAJUnvTP0YFWUTcaM4sFh5jTKr/fhIwFlgFvAZM9u82GXjVf/814HxjTIIxpi+QD3wU6OAiIu1m2yro0ttX3KR9ZAykn4ZYiIgjcc3YJwt40j8TRQzwgrX2dWPM+8ALxpiLgY3AuQDW2i+MMS8AK4AaYJq1tjY48UVE2sG2lRpe0d4yBnLU6rehZj/ExbtOIyJR5ogF2Vr7OTC8ie3bgVMP8ZybgJvanE5ExLXaati+BvLHuk4SXTIG0sHUwva1vnmRRUTakVbSExE5nB3roHa/b25eaT91Z+zrxn+LiLQjFWQRkcOpW9HNM8BtjmiT3o9qG6uCLCJOqCCLiBzOtpWA0QwW7S0unvU2E7aqIItI+1NBFhE5nNKVkJYLHZJcJ4k6q20vnUEWESdUkEVEDmfrCs1g4chqby8o/xr2VbiOIiJRRgVZRORQ9u+FHV9B98Mo8xPrAAAdKklEQVQtHirBstr615yqW8lQRKSdqCCLiBxK6UqwXshUQXahviBrmIWItDMVZBGRQ9n6he+2+yC3OaJUkfVAh2RdqCci7U4FWUTkULZ+AR06QmqO6yRRyRLjm15PZ5BFpJ2pIIuIHMqW5b5V3GL0rdKZjIEqyCLS7vRdX0SkKdbC1uUaXuFa94GwpxQqSl0nEZEoooIsItKUXZugqlwzWLhWN8Ve6Uq3OUQkqqggi4g0pf4CPRVkpzL8Z/Dr/jxERNqBCrKISFO2LvPddh/oNke0S8mAjh7feHARkXYS5zqAiEgoevs/8/iOTefEP7zXaHt2qpacbi/ZqUnkXP8Gsztk0vXj9zjrgzmNPreocIzDdCISyVSQRUSa0KdmPT0HHMeGn4x3HSVq1Rfgf38AHz7EhhmnQ2wHAHIK5xzmmSIibaMhFiIiB6quItds1gp6oSJzCNTuh7I1rpOISJRQQRYROVDpKuKMV1O8hYrMIb7bLcvc5hCRqKGCLCJyIM1gEVq65UNswrcXToqIBJkKsojIgbZ+QaWNh7Rc10kEIDbONx+yziCLSDtRQRYROdDW5ay2PSEm1nUSqZM52DfVm7Wuk4hIFFBBFhFpyL/E9Cpvb9dJpKHMobC3DHZvcZ1ERKKACrKISEMVW2HvdlZZFeSQUjcefKsWDBGR4FNBFhFpyF/AVJBDTN2Ue1s+d5tDRKKCCrKISEP+C8FWaohFaEnsAqm9teS0iLQLFWQRkYY2fw5derOTFNdJ5ECZQzWThYi0CxVkEZGGtnwOWUNdp5CmdB8M29fC/j2uk4hIhFNBFhGps283bP/Kd6ZSQk/mEMDCtpWuk4hIhFNBFhGps2U5YCFrmOsk0pT6C/U0zEJEgksFWUSkTt0MCRpiEZpS+0BCZxVkEQk6FWQRkTqbP4eOHuiU5TqJNMUY3zhkzYUsIkEW5zqAiIhLo2bOo6S8EoA34hdSarOYfP0bZKcmOU4mTcocAp/8DYPXdRIRiWAqyCIS1UrKK9kwczzU7IObfw7fPYcNY8e7jiWHkjkEqveQY7a6TiIiEUxDLEREwDczgrdGF+iFuh5HAzDErHccREQimQqyiAjA5s98t5riLbR5BkBsAoNjVJBFJHhUkEVEwDeDRXwn6NrXdRI5nNgOkDlYZ5BFJKhUkEVEwDeDRdZQiNG3xZCXdTSDYtaDVxfqiUhw6H8CERFvrW/qMA2vCA9Zw+hsKuEbnUUWkeBQQRYR2b4WqvfqAr1w4b9Qj02fuM0hIhFLBVlEZLNW0Asrnu+wz8bB5k9dJxGRCKWCLCKy+VOIS4T0fq6TSHPExbPK9oZNKsgiEhwqyCIiJR/7xh/HdnCdRJppubev78y/ta6jiEgEUkEWkagWg9c3B3L2Ma6jSAsss31h307Ysc51FBGJQCrIIhLV8kwJVO+BHsNdR5EWWOb1z1etC/VEJAhUkEUkqg2L+cp3p4fOIIeTL20viI3XhXoiEhQqyCIS1YaadZDQGbrluY4iLVBNHGQM1IV6IhIUKsgiEtWGxqzzzX+sFfTCT4+jdaGeiATFEf9HMMb0Msb8xxiz0hjzhTHmV/7tacaYt40xa/y3XRs853pjzFpjzGpjzPeD+QJERFqtZh/fMV/rAr0wlJ2aROEHcbBvJ9/77WPkFM4hp3AOo2bOcx1NRCJAXDP2qQGuttZ+bIzpBCw1xrwNXAC8Y62daYwpBAqB64wxA4HzgUFAD2CuMaaftbY2OC9BRKSVti4n3tRq/HEYWlQ4Bjanw8OPMv/HnWHoeAByCuc4TiYikeCIZ5CttZuttR/77+8GVgLZwATgSf9uTwIT/fcnAM9Za/dZa9cDa4HjAx1cRKTNSj723eoMcnjKGAgdkqF4ieskIhJhWjTozhiTAwwHPgS6W2s3g69EAxn+3bKBogZPK/ZvO/BYU4wxS4wxS0pLS1ueXESkrTZ9QpntDF16uU4irREbB1lHQ4kKsogEVrMLsjEmBXgZuNJau+twuzax7aArKKy1j1hrC6y1BR6Pp7kxREQCp+RjPvfmgmnq25aEhZ4FsGUZ1OxznUREIkizCrIxpgO+cvy0tfYV/+atxpgs/+ezgG3+7cVAw9MxPYFNgYkrIhIg+yqgbDWf21zXSaQtehZA7X5fSRYRCZDmzGJhgMeAldbavzT41GvAZP/9ycCrDbafb4xJMMb0BfKBjwIXWUQkADZ/BtbLZ96jXCeRtsgu8N1qHLKIBFBzZrEYBUwClhlj6mZk/y0wE3jBGHMxsBE4F8Ba+4Ux5gVgBb4ZMKZpBgsRCTmbfBfofe7VGeSw1iUbOmVpHLKIBNQRC7K19j2aHlcMcOohnnMTcFMbcomIBFfxYujSm+1VXVwnkbbKPlZnkEUkoLR0lIhEH2uh6CPopRkoI0LPAvhmPezZ7jqJiEQIFWQRiT47i2H3ZhXkSFE3DrlkqdscIhIxVJBFJPoU+68bVkGODD2Gg4nROGQRCRgVZBGJPkUfQVwSdB/sOokEQkIKeL6jccgiEjAqyCISfYo+8l3YFdvBdRIJlJ7H+odYHLQulYhIi6kgi0h0qa6ELZ9Dr+NcJ5FAyi6AqnJyzWbXSUQkAqggi0h02fQJeGugp8YfR5ReIwA4NuZLx0FEJBKoIItIdCnSBXoRKb0fJHWlwKggi0jbqSCLSHQp+gjScqFjuuskEkgxMdBrJAUxq10nEZEIoIIsItHDWt8Ub/5fx0uE6T2So2I2Q0Wp6yQiEuZUkEUkenyzHvaUQk9doBeRep/guy36wG0OEQl7KsgiEj2KFvtudQY5MvU4mn22A2xUQRaRtlFBFpHoUfQhxKdAxndcJ5FgiEvgU3sUbHzfdRIRCXMqyCISPb7+r+/scUys6yQSJEu8/WDzZ7B/j+soIhLGVJBFJDrs2Q6lK6HPd10nkSBa7O3vm+e6ZKnrKCISxlSQRSQ61P3aPedEtzkkqD725gMGNn7oOoqIhDEVZBGJDl8vgrhE6DHcdRIJol2kQMZAjUMWkTaJcx1ARKRdfL3IN71bXILrJBJE2alJPLUpi4mxixhW+E+8/vNA2alJLCoc4zidiIQLFWQRiXxVO2HLMjj5N66TSJAtKhwDn2+HV+ay7le9IWsYADmFcxwnE5FwoiEWIhL5Nn4I1gt9RrlOIu2h7kLMDe+5zSEiYUtnkEUk4j31/DOcb2MZ8sh2qmh8JjE7NclRKgmaLtmQlgvrF8IJ01ynEZEwpIIsIhFvYPVyOvQ+llW/OMd1FGkvfU+G5a9AbQ3E6r86EWkZDbEQkci2fw9DzTrNfxxtck6Cfbtgy2euk4hIGFJBFpHI9vX7dDC1vjOKEj3q/rzXL3SbQ0TCkgqyiES29fPZZ+Og9wmuk0h7SskAzwBYv8B1EhEJQyrIIhLZ1s3nY28/iE92nUTaW85JsPEDqNnvOomIhBkVZBGJXHu2w5ZlLPIOcp1EXOh7MlTvgU0fu04iImFGBVlEItf6dwFY5B3sOIg4kXMiYDQOWURaTAVZRCLX+nchoTOf21zXScSF5DTIHFz/g5KISHOpIItI5Fr3LuScSC2xrpOIK7mnwMYPSKbKdRIRCSMqyCISmb75Gr5ZD31Hu04iLuWdCt5qToj5wnUSEQkjWl5IRCLGqJnzKCmvBOBHsf/htg4w9lWj5aSjWe8ToEMyJ9d87jqJiIQRFWQRiRgl5ZVsmDne9+CF56Eoi7kzpoAxboOJO3EJkHMiJ61WQRaR5tMQCxGJPLU18NV8yBurcixw1KnkxmyBbza4TiIiYUIFWUQiT/FHsG8n5J/mOomEgrxTfbdr33GbQ0TChgqyiESeNf+GmDjI/Z7rJBIKuuVRbNPhq3muk4hImFBBFpHIs2Yu9BoJiV1cJ5FQYAwLaof6pv2rrXadRkTCgAqyiESWXZtg6zLIH+s6iYSQd71DYf9u2PiB6ygiEgZUkEUksqyd67vNP91tDgkpC71DITYevnzTdRQRCQMqyCISWda8DZ16QMZA10kkhOwlEfqeDKvfAGtdxxGREKeCLCIRowM1sG6+b3iFpneTA/UbBzvWQdka10lEJMSpIItIxDgh5gvYtwv6j3cdRUJR/zN8t6vfcJtDREKeCrKIRIzTY5ZAh46a3k2a1qUnZA7ROGQROSIVZBGJDF4vp8Uu9Q2v6JDoOo2Eqv7/A0Ufwp7trpOISAhTQRaRyFCylO6mHAac6TqJhLL+Z4D1+haTERE5BBVkEYkMq/5JtY3V9G5yeFlH+2Y5WfW66yQiEsLijrSDMeZx4Exgm7V2sH9bGvA8kANsAH5krf3G/7nrgYuBWuAKa+1bQUkuIlFr1Mx5lJRXNthimRf/PGVxgzk+KdVZLgkDxsB3zoKlT8C+3ZDQyXUiEQlBzTmD/AQw7oBthcA71tp84B3/Y4wxA4HzgUH+5zxgjIkNWFoREaCkvJINM8d/+3F1PrkxWzh+3CTX0SQcDJoItfvgS52/EZGmHfEMsrV2gTEm54DNE4Dv+e8/CcwHrvNvf85auw9Yb4xZCxwPvB+YuCIiTVjxKmBggKZ3k6ZlpyaRUzgHAIOXDxJS+eSFh/h/c9JYVDjGcToRCTVHLMiH0N1auxnAWrvZGJPh354NNFzovti/7SDGmCnAFIDevXu3MoaIRD1rYdlL0Oe70LmH6zQSog4qwW8sZNzHT3J1+Q43gUQkpAX6Ir2mlq5qck1Pa+0j1toCa22Bx+MJcAwRiRpbv4Cy1TD4HNdJJJwMmgg1VZwa84nrJCISglpbkLcaY7IA/Lfb/NuLgV4N9usJbGp9PBGRI1j+EphYGDjRdRIJJ71GQkom/xP7oeskIhKCWluQXwMm++9PBl5tsP18Y0yCMaYvkA981LaIIiKHYC0sf9m3cl7HdNdpJJzExMDA/+V7MZ9C1S7XaUQkxByxIBtjnsV3kV1/Y0yxMeZiYCZwmjFmDXCa/zHW2i+AF4AVwJvANGttbbDCi0iUK14C5RthyA9dJ5FwNORcEk01rHzNdRIRCTHNmcXix4f41KmH2P8m4Ka2hBIRaZblL0FsgmavkNbpeRzrvJnkfvYcDP+Z6zQiEkJaO4uFiIhbtdWw/BXIPw0Su7hOI+HIGObGj2HKhmcYVfgEJXx7wXh2apKmfxOJYirIIhKe1vwb9mzTmT9pkynTCuHuZ1h0RimMvqB+e92cySISnQI9zZuISPv45G+Q0h3yTnOdRMJZ1z7Q50T47FnfRZ8iIqggi0gY8lDuWyZ42PkQq1+ESRsNOx92fOW76FNEBBVkEQlDZ8cuBFsLR2t4hQTAwAkQlwSfzHadRERChAqyiIQXa/lR7LvQawR4+rlOI5EgsTMM+YFvyfKqna7TiEgIUEEWkfBS9CF5MZtg+CTXSSSSFFwM1Xvhs+ddJxGREKCCLCLh5aNH2GWTYNDZrpNIJMk+BnoMhyWP6WI9EVFBFpEwsmszrHiVF2u/BwkprtNIpCm4GEpXwcb3XScREcdUkEUkfCydBd5aZtdqajcJgsHnQEIXWPyY6yQi4pgKsoiEh5r9sGQW5J/O1zbTdRqJRPEd4eifwIpX6c4O12lExCEVZBEJDyv+4Vs5b8QU10kkko28FGwtF8a95TqJiDikgiwioc9aeP9+6JYHuWNcp5FI1jUHBk7gJ7FzoWqX6zQi4ogKsoiEvq/egc2fwqhfQYy+bUmQffcKOptKWPqE6yQi4ojWaBWRkDVq5jxKyit5Pv6P9DJpjH6hC9UvzCE7Ncl1NIlk2cfwfu1ATvjgQRhxKcTFu04kIu1MBVlEQlZJeSUbpqbBrFUwbiZrRk5wHUmixMO14zlh959h2QswXEuai0Qb/a5SRELbwjsguRscM9l1Eoki871HQ9YwePc2qK12HUdE2pkKsoiErKHmK1j7NoycCvHJruNIVDFwynQo/xo++ZvrMCLSzjTEQkRCk7VcH/csJKfD8Ze4TiNRJjs1iZzH9/JKfB6Z//wjp7yUyj7iyU5NYlGhZlIRiXQqyCISmta+wwmxK2D0bZDY2XUaiTL1JXhdCsyewOqzt8GIS8gpnOM2mIi0Cw2xEJHQ4/XC3D+w0euBYy90nUaiWd/RkHOSbyxy1U7XaUSknaggi0joWfYibF3G7TXnaYotccsYOP3/wd7tvpIsIlFBBVlEQkvVLnj7Bsg6mn96R7pOIwI9hvumevvwIXLNJtdpRKQdqCCLSGiZPxMqtsL4v2D1LUpCxak3QIdkfh/3lOskItIO9L+PiISOLcvgw4eg4ELoeazrNCLfSsmA0ddySuxnsPJ112lEJMhUkEUkNHi9MOdqSOrqO1snEmqOv4SV3t4w59ewd4frNCISRCrIIhIaPngAij6E0//kK8kioSYunmuqL4E9ZfDWb12nEZEgUkEWEfe2roB3boQBZ8Kw812nETmkL2xfOPEq+OxZ+PIt13FEJEhUkEXErZp98Mr/QWIXOOtu37RaIqFs9LWQMRBevRx2b3WdRkSCQAVZRNx6ewZsXQ7/ex90THedRuSwslOTyJk+l9OLJlNZUc77t00kt/CfjJo5z3U0EQkgLTUtIu589jx8+CCMuBT6j3OdRuSI6pegBvi0Myf8YyrrTvmMnH8PdxdKRAJOBVlE3Nj8GfzzCuhzIid/eiob351z0C7ZqUkOgok009E/gQ2LYMHtnB5zJTDedSIRCRAVZBFpf7s2w3M/heRucO4sNv5pMRtmqlxIGBp/O5Su4u7i+6H4LM3fLRIhNAZZRNpX5Tfwt3N8t+c/7VuAQSRcdUiCHz9Hqe0Cz54HO9a7TiQiAaAzyCLSbsbc8i9uq7yBIWYdF1Rfx/v3bAI2aSiFhLcUD9clzuDBiuvYdddpnL//92zCd8FpdmpS43HLIhIWVJBFpH1U7eKWyj9QELMGzn2CZwdNdJ1IJGCe/e0kKBlI6uyz+W/aHXDBHEjtRU7hwWPrRST0aYiFiATfnu3w5FkcY9bADx4FlWOJRNnHws//DpXl8MT/QNka14lEpJVUkEUkuMrWwuPfh20rmVL9axjyQ9eJRIIn+1j4+T9g/154dCwjzErXiUSkFVSQRSR41rwNfx0DlTtg0t/5j1dzxUoUyD4GfjEXUjJ4Kv5mWDILrHWdSkRaQAVZRAKvZj+880d4+lzo2humzIecUa5TibSftL5w8b/5wDsQXr8SXroIqna6TiUizaSL9ESk1UbNnEdJeWWjbf3NRu5NfJh+dj0MnwRn3AbxyY4SijiU1JXJ1dex/vQvYd5NULIExt8J+WNdJxORI1BBFpFWKymv/HaBj8pymD8TPnqEb2xHfrH/aua+fyy8/5/6/TWdm0QbSwycdDXknAT/uAye/gEM/gGc/ifo3MN1PBE5BBVkEWmb/XthyePw3p2wdzscewFdx/yeRzt2c51MJHT0Oh6mLoL37oKFt8OqOTDiEhh1JSSnuU4nIgdQQRaR1qnYxmWx/4C7fwV7SqHvyXDa/4MeR7tOJhIyslOTDpgLeTA9zW1cFfcSZ793D3vfe5gXa0fzRO33+dpm1j9Hi4uIuGVsCFxZW1BQYJcsWeI6hogcSW0NbFgAHz8FK/8J3mrIPQVGXwd9TnCdTiS8bP0C/nsvLHsJvDWQfxoMPY/vPG1YOfMHrtOJRAVjzFJrbcFB21WQReSwKsvh6//C6jd8vxau3AGJXeDonzLm3Vzm3fJ/rhOKhLfdW2Dxo/DJ07B7E3tJ5O3aY5hfO4yF3qGU0QXQmWWRYDhUQdYQCxEB6mak2Etvs41BZgPDYr7ihJgVDDbriTWW3TaJud5j+Fft8bxbNYx98+N10Z1IIHTKhDHT4Xu/hY3/JXnZi0xY+ToT9v7X9/nuQ6DXcVzzfgcozYZu+RCjWVpFgiloZ5CNMeOAu4FY4FFr7cxD7aszyCLtqLYa9pTB7s3wzXrY4fv48OMljEjaBPt2+faL6QA9j4O+J/nGF/c8DuIS3GYXiRZeL2z5HNbOhfULYNMn3/7bjEuC9DxI7wfp/aHbUb4ZMTpl+W7171Sk2dp1iIUxJhb4EjgNKAYWAz+21q5oan8V5DB2qL8/h/x71dL9W/OcCPoatTVQu9831re27mO/79bb4H7tfti/B/bt9n3sr4B9Fb7/UPdX+GaXqCiFPdt89w+U0p3Fu7pw3PEnQdZQyBwKGQOhQ+JhXo+ItBuvl5/OfIoee5bT3xRxlNlEntlEtikjxhzwPSS5G6R0h8RUSEptfJvYGeISoUOSr0jHJfn+ndfdxiZATJzvDHVMHJhYiIn13/dvq38c69sGYEzjW5Ew0d5DLI4H1lpr1/m/+HPABKDJguzM/SOhfGMTnwhQWWrNc0L1a0jY2UsiySmpkJDi+w+z21G+C+lSukNHj+/Xul1zfB/xHTm3cA4bzhzvOraINCUmhqd/O/ng7fv38rPbX8BUbCbT7CCTHWTu+ob03TvpYr6hC8V0Nnvowh5STFX758Y0KM0HlugjfK7+sUp3xMsaBhf9y3WKRoJVkLOBogaPi4ERDXcwxkwBpvgfVhhjVgcpy5GkA2WOvnak03sbPM14b3cB21p0UHNrq/NEGv3dDR69t8Gj9za49P4GzeZ0Ljau3ts+TW0MVkFu6se9RqcnrbWPAI8E6es3mzFmSVOn1qXt9N4Gj97b4NL7Gzx6b4NH721w6f0NnlB8b4N1GWwx0KvB457ApiB9LRERERGRgAlWQV4M5Btj+hpj4oHzgdeC9LVERERERAImKEMsrLU1xpjLgbfwTfP2uLX2i2B8rQBwPswjgum9DR69t8Gl9zd49N4Gj97b4NL7Gzwh996GxEp6IiIiIiKhQkvxiIiIiIg0oIIsIiIiItJA1BZkY8w4Y8xqY8xaY0yh6zyRxBjzuDFmmzFmuesskcYY08sY8x9jzEpjzBfGmF+5zhQpjDGJxpiPjDGf+d/bG11nijTGmFhjzCfGmNddZ4k0xpgNxphlxphPjTFamjaAjDGpxpiXjDGr/N97T3CdKVIYY/r7/87WfewyxlzpOhdE6Rjkli6FLS1jjDkZqABmW2sHu84TSYwxWUCWtfZjY0wnYCkwUX93284YY4CO1toKY0wH4D3gV9baDxxHixjGmF8DBUBna+2ZrvNEEmPMBqDAWquFLALMGPMksNBa+6h/Zq5ka22561yRxt/NSoAR1tqvXeeJ1jPI9UthW2v3A3VLYUsAWGsXADtc54hE1trN1tqP/fd3AyvxrVwpbWR9KvwPO/g/ou8MQpAYY3oC44FHXWcRaS5jTGfgZOAxAGvtfpXjoDkV+CoUyjFEb0FuailslQwJK8aYHGA48KHbJJHDPwTgU3xrdL9trdV7Gzh3AdcCXtdBIpQF/m2MWWqMmeI6TATJBUqBWf7hQY8aYzq6DhWhzgeedR2iTrQW5CMuhS0SyowxKcDLwJXW2l2u80QKa22ttfZofKt/Hm+M0RChADDGnAlss9YudZ0lgo2y1h4DnAFM8w91k7aLA44BHrTWDgf2ALpuKcD8Q1f+F3jRdZY60VqQtRS2hC3/+NiXgaetta+4zhOJ/L9CnQ+McxwlUowC/tc/TvY5YIwx5m9uI0UWa+0m/+024O/4hhJK2xUDxQ1+m/QSvsIsgXUG8LG1dqvrIHWitSBrKWwJS/4LyR4DVlpr/+I6TyQxxniMMan++0nAWGCV21SRwVp7vbW2p7U2B9/323nW2p85jhUxjDEd/Rft4v/1/+mAZhEKAGvtFqDIGNPfv+lUQBdFB96PCaHhFRCkpaZDXZgthR12jDHPAt8D0o0xxcAMa+1jblNFjFHAJGCZf6wswG+ttW84zBQpsoAn/VdSxwAvWGs1HZmEg+7A330/PxMHPGOtfdNtpIjyS+Bp/wm1dcCFjvNEFGNMMr5ZxS5xnaWhqJzmTURERETkUKJ1iIWIiIiISJNUkEVEREREGlBBFhERERFpQAVZRERERKQBFWQRERERkQZUkEVEREREGlBBFhERERFp4P8DQgafgXEMJS8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = {'Entries': len(y_all),\n",
    "     'Mean': y_all.mean(),\n",
    "     'Std': y_all.std(ddof=1),\n",
    "     'Chi2': chi2,\n",
    "     'NDF': N_dof,\n",
    "     'Prob': chi2_prob, \n",
    "     'N': [minuit.values['N'], minuit.errors['N']], \n",
    "     'mu': [minuit.values['mu'], minuit.errors['mu']], \n",
    "     'sigma': [minuit.values['sigma'], minuit.errors['sigma']],\n",
    "    }\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "add_text_to_ax(0.02, 0.97, text, ax2, fontsize=14);\n",
    "\n",
    "fig2.tight_layout()\n",
    "fig2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "if save_plots:\n",
    "    fig2.savefig(\"Dist_ErrorProp.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------------------------------------------------------------------------- \n",
    "\n",
    "# Questions:\n",
    "\n",
    "0. First solve the problem of obtaining the Perimeter, Area & Diagonal with uncertainty ANALYTICALLY.\n",
    "\n",
    "Example answer 0: See cell 4. No need to differentiate by hand, SymPy can do it easily for you.\n",
    "\n",
    "---\n",
    "\n",
    "1. Now look at the program, and assure yourself that you understand what is going on. Put in the correct expression for y in terms of x1=L and x2=W in order to calculate the circumference, area, and diagonal length, and run the program. Does the output correspond well with the results you expected from your analytical calculations to begin with?\n",
    "\n",
    "Example answer 1:\n",
    "Yes, both the mean value and the uncertainty are in a very good agreement to the results from analytical calculations. In the area calculation, there is a slight discrepancy, simply because multiplying two Gaussians does not give (exactly) a Gaussian... only when adding! \n",
    "Also, one little thing to note is, that the width can actually come out negatively in rare cases (i.e. beyond 4 sigma low). When producing a simulation of cases, one should of course be aware of such things!\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "2. Imagine that you wanted to know the central value and uncertainty of y1 and y2, given the\n",
    "   same above PDFs for `x1`=$L$ and `x2`=$W$:\n",
    "   \n",
    "     `y1 = log(square(x1*tan(x2))+sqrt((x1-x2)/(cos(x2)+1.0+x1)))`\n",
    "     \n",
    "     `y2 = 1.1+sin(20*x1)`\n",
    "\n",
    "   Get the central value of y, and see if you can quickly differentiate this with\n",
    "   respect to `x1` and `x2`, and thus predict what uncertainty to expect for y using\n",
    "   the error propagation formula. It is (for once) OK to give up on the first expression :-)\n",
    "   Next, try to estimate the central value and uncertainty using random numbers\n",
    "   like above - do you trust this result more? And are the distributions Gaussian?\n",
    "   \n",
    "   \n",
    "Example answer 2:\n",
    "The first expression \"y1\" looks horrible, and it is! This is an obvious case where simulation is probably a better way forward, and at least a good cross check.\n",
    "The second expression \"y2\" is simple, but made to make the error propagation formula break down! The function varies wildly, and as a result, the error propagation formula assumption is broken (badly). The resulting error distribution is also very far from Gaussian, as can be seen below.\n",
    "\n",
    "For code solution example, see cell 19 below.\n",
    "\n",
    "\n",
    "### Advanced questions:\n",
    "\n",
    "3. Try to generate `x1` and `x2` with e.g. a quadratic correlation, and see that despite\n",
    "   not having any linear correlation, the result on perimeter, area, and diagonal\n",
    "   length is still affected.\n",
    "\n",
    "---------------------------------------------------------------------------------- "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$y_{1} = \\frac{\\sqrt{x_{1} - x_{2}}}{x_{1} + \\cos{\\left(x_{2} \\right)} + 1} + \\log{\\left(\\sqrt{x_{1} + \\tan{\\left(x_{2} \\right)}} \\right)}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$y_{2} = \\sin{\\left(20 x_{1} \\right)} + 1.1$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\sigma_{y 1} = \\sqrt{\\sigma_{x 1}^{2} \\left(- \\frac{\\sqrt{x_{1} - x_{2}}}{\\left(x_{1} + \\cos{\\left(x_{2} \\right)} + 1\\right)^{2}} + \\frac{1}{2 \\left(x_{1} + \\tan{\\left(x_{2} \\right)}\\right)} + \\frac{1}{2 \\sqrt{x_{1} - x_{2}} \\left(x_{1} + \\cos{\\left(x_{2} \\right)} + 1\\right)}\\right)^{2} + \\sigma_{x 2}^{2} \\left(\\frac{\\sqrt{x_{1} - x_{2}} \\sin{\\left(x_{2} \\right)}}{\\left(x_{1} + \\cos{\\left(x_{2} \\right)} + 1\\right)^{2}} + \\frac{\\frac{\\tan^{2}{\\left(x_{2} \\right)}}{2} + \\frac{1}{2}}{x_{1} + \\tan{\\left(x_{2} \\right)}} - \\frac{1}{2 \\sqrt{x_{1} - x_{2}} \\left(x_{1} + \\cos{\\left(x_{2} \\right)} + 1\\right)}\\right)^{2}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\sigma_{y 2} = 20 \\sqrt{\\sigma_{x 1}^{2} \\cos^{2}{\\left(20 x_{1} \\right)}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$y1 = (1.1 \\pm 0.1)\\,\\mathrm{m}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$y2 = (1.9 \\pm 5.1)\\,\\mathrm{m}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Mean = 1.101,    RMS = 0.705\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAFlCAYAAAApo6aBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAY30lEQVR4nO3db6ye510f8O+PpE0NBTlZnGCcBAfNMCWVaDsrKouEWMPW0CKcFwsyEsxDkbxN2YBp0nB4McQLS540oTGNbrKAzWjQYP50sVr+BWcVWgUNTinQJM1qSEiceLHpCKUjCkv47cW5M47tc3Ke43Mdn+ccfz7S0XM/13Pdz/M7l26f5+vr/lfdHQAA1u7LNroAAICtQrACABhEsAIAGESwAgAYRLACABhEsAIAGOTajS4gSW688cbevXv3RpcBALCiJ5544k+6e8dSr81FsNq9e3dOnTq10WUAAKyoqv54uddm2hVYVf+iqp6sqs9W1Ueq6h1VdUNVPVpVn58er1/U/6GqOl1Vz1TVB0b8EgAA827FYFVVu5J8X5K93f2uJNck2Z/kUJKT3b0nycnpearqjun1O5Pcm+TDVXXN+pQPADA/Zj14/dok26rq2iRfnuSlJPuSHJteP5bkvml5X5KHu/u17n42yekkd40rGQBgPq0YrLr7xST/NsnzSc4m+bPu/vUkN3f32anP2SQ3TavsSvLCorc4M7VdoKoOVtWpqjp1/vz5tf0WAABzYJZdgddnYRbq9iRfk+Qrquq732qVJdouudNzdx/t7r3dvXfHjiUPrAcA2FRm2RX4rUme7e7z3f1/k/xSkr+T5OWq2pkk0+O5qf+ZJLcuWv+WLOw6BADY0mYJVs8neV9VfXlVVZJ7kjyd5ESSA1OfA0kemZZPJNlfVddV1e1J9iR5fGzZAADzZ8XrWHX3p6rqF5J8OsnrSX43ydEk70xyvKoeyEL4un/q/2RVHU/y1NT/we5+Y53qBwCYG9V9yeFPV9zevXvbBUIBgM2gqp7o7r1LveZegQAAgwhWAACDCFYAAIMIVgAAg6x4ViAAXK67jzyWF1959YK2Xdu35ZOH3r9BFcH6EqwAWDcvvvJqnjvyoQvadh/6+AZVA+vPrkAAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEFWDFZV9Q1V9ZlFP1+sqh+oqhuq6tGq+vz0eP2idR6qqtNV9UxVfWB9fwUAgPmwYrDq7me6+93d/e4kfzvJXyT5aJJDSU52954kJ6fnqao7kuxPcmeSe5N8uKquWaf6AQDmxmp3Bd6T5A+7+4+T7EtybGo/luS+aXlfkoe7+7XufjbJ6SR3jSgWAGCerTZY7U/ykWn55u4+myTT401T+64kLyxa58zUdoGqOlhVp6rq1Pnz51dZBgDA/Jk5WFXV25N8R5KfX6nrEm19SUP30e7e2917d+zYMWsZAABzazUzVt+W5NPd/fL0/OWq2pkk0+O5qf1MklsXrXdLkpfWWigAwLxbTbD6rvz1bsAkOZHkwLR8IMkji9r3V9V1VXV7kj1JHl9roQAA8+7aWTpV1Zcn+XtJ/vGi5iNJjlfVA0meT3J/knT3k1V1PMlTSV5P8mB3vzG0agCAOTRTsOruv0jyNy5q+0IWzhJcqv/hJIfXXB0AwCbiyusAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAg1y70QVsFncfeSwvvvLqBW27tm/LJw+9f4MqAgDmjWA1oxdfeTXPHfnQBW27D318g6oBAOaRXYEAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDzBSsqmp7Vf1CVX2uqp6uqm+qqhuq6tGq+vz0eP2i/g9V1emqeqaqPrB+5QMAzI9ZZ6x+LMmvdvffSvKNSZ5OcijJye7ek+Tk9DxVdUeS/UnuTHJvkg9X1TWjCwcAmDcrBquq+qok35zkJ5Oku/+yu19Jsi/JsanbsST3Tcv7kjzc3a9197NJTie5a3ThAADzZpYZq69Lcj7Jf66q362qn6iqr0hyc3efTZLp8aap/64kLyxa/8zUdoGqOlhVp6rq1Pnz59f0SwAAzINZgtW1Sd6b5D9293uS/J9Mu/2WUUu09SUN3Ue7e293792xY8dMxQIAzLNZgtWZJGe6+1PT81/IQtB6uap2Jsn0eG5R/1sXrX9LkpfGlAsAML9WDFbd/b+SvFBV3zA13ZPkqSQnkhyY2g4keWRaPpFkf1VdV1W3J9mT5PGhVQMAzKFrZ+z3z5P8TFW9PckfJfneLISy41X1QJLnk9yfJN39ZFUdz0L4ej3Jg939xvDKAQDmzEzBqrs/k2TvEi/ds0z/w0kOr6EuAIBNx5XXAQAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGEawAAAYRrAAABhGsAAAGmSlYVdVzVfUHVfWZqjo1td1QVY9W1eenx+sX9X+oqk5X1TNV9YH1Kh4AYJ6sZsbq73b3u7t77/T8UJKT3b0nycnpearqjiT7k9yZ5N4kH66qawbWDAAwl9ayK3BfkmPT8rEk9y1qf7i7X+vuZ5OcTnLXGj4HAGBTmDVYdZJfr6onqurg1HZzd59Nkunxpql9V5IXFq17Zmq7QFUdrKpTVXXq/Pnzl1c9AMAcuXbGfnd390tVdVOSR6vqc2/Rt5Zo60sauo8mOZoke/fuveR1AIDNZqYZq+5+aXo8l+SjWdi193JV7UyS6fHc1P1MklsXrX5LkpdGFQwAMK9WDFZV9RVV9ZVvLif5+0k+m+REkgNTtwNJHpmWTyTZX1XXVdXtSfYkeXx04QAA82aWXYE3J/loVb3Z/2e7+1er6neSHK+qB5I8n+T+JOnuJ6vqeJKnkrye5MHufmNdqgcAmCMrBqvu/qMk37hE+xeS3LPMOoeTHF5zdQAAm4grrwMADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADDJzsKqqa6rqd6vqY9PzG6rq0ar6/PR4/aK+D1XV6ap6pqo+sB6FAwDMm9XMWH1/kqcXPT+U5GR370lycnqeqrojyf4kdya5N8mHq+qaMeUCAMyvmYJVVd2S5ENJfmJR874kx6blY0nuW9T+cHe/1t3PJjmd5K4x5QIAzK9ZZ6z+XZJ/leSvFrXd3N1nk2R6vGlq35XkhUX9zkxtF6iqg1V1qqpOnT9/ftWFAwDMmxWDVVV9e5Jz3f3EjO9ZS7T1JQ3dR7t7b3fv3bFjx4xvDQAwv66doc/dSb6jqj6Y5B1Jvqqq/muSl6tqZ3efraqdSc5N/c8kuXXR+rckeWlk0QAA82jFGavufqi7b+nu3Vk4KP2x7v7uJCeSHJi6HUjyyLR8Isn+qrquqm5PsifJ48MrBwCYM7PMWC3nSJLjVfVAkueT3J8k3f1kVR1P8lSS15M82N1vrLlSAIA5t6pg1d2fSPKJafkLSe5Zpt/hJIfXWBsAwKbiyusAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAINcu9EFbGa7tm/L7kMfv6Ttk4fev0EVAQAbSbBag6UC1MVBCwC4etgVCAAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwyIrBqqreUVWPV9XvVdWTVfUjU/sNVfVoVX1+erx+0ToPVdXpqnqmqj6wnr8AAMC8mGXG6rUk7+/ub0zy7iT3VtX7khxKcrK79yQ5OT1PVd2RZH+SO5Pcm+TDVXXNehQPADBPVgxWveBL09O3TT+dZF+SY1P7sST3Tcv7kjzc3a9197NJTie5a2jVAABz6NpZOk0zTk8k+ZtJfry7P1VVN3f32STp7rNVddPUfVeS3160+pmpDYAt7O4jj+XFV169oG3X9m0bVA1sjJmCVXe/keTdVbU9yUer6l1v0b2WeotLOlUdTHIwSW677bZZygBgjr34yqt57siHNroM2FCrOiuwu19J8oksHDv1clXtTJLp8dzU7UySWxetdkuSl5Z4r6Pdvbe79+7YseMySgcAmC+znBW4Y5qpSlVtS/KtST6X5ESSA1O3A0kemZZPJNlfVddV1e1J9iR5fHThAADzZpZdgTuTHJuOs/qyJMe7+2NV9VtJjlfVA0meT3J/knT3k1V1PMlTSV5P8uC0KxEAYEtbMVh19+8nec8S7V9Ics8y6xxOcnjN1QEAbCKuvA4AMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMMhM9wpkdru2b8vuQx+/pO2Th96/QRUBAFeKYDXYUgHq4qAFAGxNdgUCAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAziOlYArNrdRx7Li6+8ekHbru3bNqgamB+CFQCr9uIrr+a5Ix/a6DJg7tgVCAAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIgrr18Bu7Zvy+5DH7+k7ZOH3r9BFQEA60GwugKWClAXBy0AYPOzKxAAYBDBCgBgkBV3BVbVrUl+OslXJ/mrJEe7+8eq6oYkP5dkd5Lnknxnd//ptM5DSR5I8kaS7+vuX1uX6gFYV3cfeSwvvvLqJe27tm/bgGpg/s1yjNXrSf5ld3+6qr4yyRNV9WiSf5TkZHcfqapDSQ4l+cGquiPJ/iR3JvmaJL9RVV/f3W+sz68AwHp58ZVX89yRD210GbBprLgrsLvPdvenp+U/T/J0kl1J9iU5NnU7luS+aXlfkoe7+7XufjbJ6SR3jS4cAGDerOoYq6raneQ9ST6V5ObuPpsshK8kN03ddiV5YdFqZ6Y2AIAtbeZgVVXvTPKLSX6gu7/4Vl2XaOsl3u9gVZ2qqlPnz5+ftQwAgLk1U7CqqrdlIVT9THf/0tT8clXtnF7fmeTc1H4mya2LVr8lyUsXv2d3H+3uvd29d8eOHZdbPwDA3FgxWFVVJfnJJE93948ueulEkgPT8oEkjyxq319V11XV7Un2JHl8XMkAAPNplrMC707yPUn+oKo+M7X9UJIjSY5X1QNJnk9yf5J095NVdTzJU1k4o/BBZwQCAFeDFYNVd/+PLH3cVJLcs8w6h5McXkNdAFxhS12zyvWqYHXcKxCAJK5ZBSMIVhtk1/Ztl9yIedf2bUvesBkA2BwEqw2yVIC6OGgBrBe7/WB9CFZzbrk/fma24Ooy+m+B3X6wPgSrObfUHz8zW3D18bcANodV3dIGAIDlCVYAAIMIVgAAgzjGahNyqQYAmE+C1SbkUg3Aari0Alw5gtUWsdQs1pvtZrLg6ubSCnDlCFZbxHLhyUwWbF0OC4D5I1gBbFIOC4D5I1jNkeX+9wkAbA6C1RwxfQ8Am5vrWAEADGLGCmDOrOXyCA4pgI0lWAHMmbVcHsEhBbCx7AoEABjEjBVJlt/14H+/ADA7wYokS+96cD0cAFgdweoq5L5hMD/8e4StRbC6CrlvGMwP/x5haxGstjinXgPAlSNYbXEOPgeAK0ewYlnLzXZdHNacUQgACwQrlrVUMLr7yGNLhi1nFAKAYMUqmYWCy+cMQNj6BCuuGLsM184YXlmjx9sZgLD1CVasi+WOz7LLcG1cyHXtVhOWZh1vgRd4k2DFuvCFsmAtX7h2G62P9QinAi/wphWDVVX9VJJvT3Kuu981td2Q5OeS7E7yXJLv7O4/nV57KMkDSd5I8n3d/WvrUjlsAmv5wrXbaHlmiIB5NcuM1X9J8h+S/PSitkNJTnb3kao6ND3/waq6I8n+JHcm+Zokv1FVX9/db4wtm61iqV2Gy/XzpcmbRs8QLbcdLjVD6KK7wFtZMVh1929W1e6Lmvcl+ZZp+ViSTyT5wan94e5+LcmzVXU6yV1JfmtMuWw1s4al5S7zME9hayvturvaZoRW83tt1TEAxrjcY6xu7u6zSdLdZ6vqpql9V5LfXtTvzNQGa7LUl9lqZihmDQpr7Xe5u+7mLZTN0zFD8zY2szKzBVen0Qev1xJtvWTHqoNJDibJbbfdNrgMuNBSQWEtFzsdffyT46mWt1nHxswWXJ0uN1i9XFU7p9mqnUnOTe1nkty6qN8tSV5a6g26+2iSo0myd+/eJcMXrKdZv/hGzzxs9ZkMZ0ICV7PLDVYnkhxIcmR6fGRR+89W1Y9m4eD1PUkeX2uRMKulvpiTtX05j555uBIzGasJN7OGmdWcaOBMSOBqNcvlFj6ShQPVb6yqM0l+OAuB6nhVPZDk+ST3J0l3P1lVx5M8leT1JA86I5D1MutFSK9GqzlGatYwMzoQjp6dmvWm4QDraZazAr9rmZfuWab/4SSH11IUzMIX5uY2enZqrSc4AIzgyutwlVjNtZquxGdv5c8Frl6CFVwlNnKGb6M+26wmcKUJVrDJmZUBmB+CFWxyZmUA5seXbXQBAABbhWAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMIhgBQAwiGAFADCIYAUAMMi6Bauqureqnqmq01V1aL0+BwBgXqxLsKqqa5L8eJJvS3JHku+qqjvW47MAAObFes1Y3ZXkdHf/UXf/ZZKHk+xbp88CAJgL6xWsdiV5YdHzM1MbAMCWde06vW8t0dYXdKg6mOTg9PRLVfXMOtWy2I1J/uRyV65/M7CSrWFN48mSjOlYxnO8IWPq7+n/Zxsd70qM6dcu98J6BaszSW5d9PyWJC8t7tDdR5McXafPX1JVneruvVfyM7cy4zmeMR3LeI5nTMcynuNt9Jiu167A30myp6pur6q3J9mf5MQ6fRYAwFxYlxmr7n69qv5Zkl9Lck2Sn+ruJ9fjswAA5sV67QpMd/9ykl9er/e/TFd01+NVwHiOZ0zHMp7jGdOxjOd4Gzqm1d0r9wIAYEVuaQMAMMiWC1Yr3UqnFvz76fXfr6r3bkSdm8kMY/otVfVnVfWZ6edfb0Sdm0VV/VRVnauqzy7zum10FWYYT9vnKlXVrVX136vq6ap6sqq+f4k+ttMZzTiettNVqKp3VNXjVfV705j+yBJ9NmYb7e4t85OFA+X/MMnXJXl7kt9LcsdFfT6Y5FeycK2t9yX51EbXPc8/M47ptyT52EbXull+knxzkvcm+ewyr9tGx46n7XP1Y7ozyXun5a9M8j/9LV338bSdrm5MK8k7p+W3JflUkvdd1GdDttGtNmM1y6109iX56V7w20m2V9XOK13oJuL2RIN1928m+d9v0cU2ugozjCer1N1nu/vT0/KfJ3k6l949w3Y6oxnHk1WYtrsvTU/fNv1cfND4hmyjWy1YzXIrHbfbWZ1Zx+ubpinZX6mqO69MaVuWbXQ82+dlqqrdSd6ThRmBxWynl+EtxjOxna5KVV1TVZ9Jci7Jo909F9voul1uYYOseCudGfvw12YZr08n+dru/lJVfTDJf0uyZ90r27pso2PZPi9TVb0zyS8m+YHu/uLFLy+xiu30LawwnrbTVeruN5K8u6q2J/loVb2ruxcfa7kh2+hWm7Fa8VY6M/bhr81ye6Ivvjkl2wvXL3tbVd145UrccmyjA9k+L09VvS0LIeBnuvuXluhiO12FlcbTdnr5uvuVJJ9Icu9FL23INrrVgtUst9I5keQfTmcLvC/Jn3X32Std6Cay4phW1VdXVU3Ld2Vhu/rCFa9067CNDmT7XL1pvH4yydPd/aPLdLOdzmiW8bSdrk5V7ZhmqlJV25J8a5LPXdRtQ7bRLbUrsJe5lU5V/ZPp9f+UhavBfzDJ6SR/keR7N6rezWDGMf0HSf5pVb2e5NUk+3s6JYNLVdVHsnAG0I1VdSbJD2fhwEvb6GWYYTxtn6t3d5LvSfIH0zEsSfJDSW5LbKeXYZbxtJ2uzs4kx6rqmiyE0OPd/bF5+L535XUAgEG22q5AAIANI1gBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAzy/wCL4fSlc01xngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y1,y2,x1,x2 = symbols('y_1,y_2,x_1,x_2')\n",
    "dy1,dy2,dx1,dx2 = symbols('sigma_y_1,sigma_y_2,sigma_x_1,sigma_x_2')\n",
    "\n",
    "y1 = log(sqrt(x1+tan(x2)))+sqrt((x1-x2))/(cos(x2)+1+x1)\n",
    "lprint(latex(Eq(symbols('y_1'),y1)))\n",
    "\n",
    "y2 = 1.1 + sin(20*x1)\n",
    "lprint(latex(Eq(symbols('y_2'),y2)))\n",
    "\n",
    "dy1 = sqrt((y1.diff(x1) * dx1)**2 + (y1.diff(x2) * dx2)**2)\n",
    "lprint(latex(Eq(symbols('sigma_y_1'),dy1)))\n",
    "\n",
    "dy2 = sqrt((y2.diff(x1) * dx1)**2 + (y2.diff(x2) * dx2)**2)\n",
    "lprint(latex(Eq(symbols('sigma_y_2'),dy2)))\n",
    "\n",
    "fy1 = lambdify((x1,x2),y1)\n",
    "fy2 = lambdify((x1,x2),y2)\n",
    "\n",
    "fdy1 = lambdify((x1,dx1,x2,dx2),dy1)\n",
    "fdy2 = lambdify((x1,dx1,x2,dx2),dy2)\n",
    "\n",
    "vx1,vdx1 = mu1,sig1\n",
    "vx2,vdx2 = mu2,sig2\n",
    "\n",
    "vy1 = fy1(vx1,vx2)\n",
    "vdy1 = fdy1(vx1,vdx1,vx2,vdx2)\n",
    "vy2 = fy2(vx1,vx2)\n",
    "vdy2 = fdy2(vx1,vdx1,vx2,vdx2)\n",
    "\n",
    "lprint(fr'y1 = ({vy1:.1f} \\pm {vdy1:.1f})\\,\\mathrm{{m}}')\n",
    "lprint(fr'y2 = ({vy2:.1f} \\pm {vdy2:.1f})\\,\\mathrm{{m}}') \n",
    "\n",
    "u = r.normal(0.0, sigu, N_exp)\n",
    "v = r.normal(0.0, sigv, N_exp)\n",
    "x1_all = mu1 + np.cos(theta)*u - np.sin(theta)*v\n",
    "x2_all = mu2 + np.sin(theta)*u + np.cos(theta)*v\n",
    "x12_all = np.array([x1_all, x2_all])\n",
    "\n",
    "y_all = 1.1 + np.sin(20*x1_all)\n",
    "\n",
    "xmin, xmax = 0.0, 3.0\n",
    "\n",
    "fig3, ax3 = plt.subplots(figsize=(10, 6));\n",
    "counts, bin_edges, _ = ax3.hist(y_all, 100, range=(xmin, xmax), histtype='step');\n",
    "bin_centers = (bin_edges[1:] + bin_edges[:-1])/2\n",
    "s_counts = np.sqrt(counts)\n",
    "\n",
    "print(f\"  Mean = {y_all.mean():5.3f},    RMS = {y_all.std(ddof=1):5.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}