Lecture1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Scientific Computing in Finance\n",
    "\n",
    "* MATH-GA 2048, Spring 2019\n",
    "* Courant Institute of Mathematical Sciences, New York University"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Instructors\n",
    "\n",
    "#### Yadong Li\n",
    "* MD, Head of Trading Central, Quantitative Analytics, Barclays\n",
    "* Ph.D. in Physics, M.S. in Computer Sciences, University of Wisconsin-Madison\n",
    "* Master of Financial Engineering, University of California, Berkeley \n",
    "\n",
    "#### Hongwei Cheng\n",
    "* Chief Risk Officer, Head of Research, Mill Hill Capital\n",
    "* Ph.D. in Applied Mathematics, New York University"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 1: Introduction\n",
    "\n",
    "## Topics\n",
    "\n",
    "* Introduction\n",
    "* Software development practices\n",
    "* Errors and floating point computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Objectives of this class\n",
    "\n",
    "* Teach fundamental principles of scientific computing\n",
    "* Teach the most common numerical techniques in quantitative Finance \n",
    "* Develop intuitions and essential skills using real world examples\n",
    "* Help students to start a career in quantitative Finance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## We expect you\n",
    "\n",
    "\n",
    "* Know basic calculus and linear algebra\n",
    "* Know basic stochastic calculus and derivative pricing theories\n",
    "* Have some prior programming experience\n",
    "* Are willing to learn through hard work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Main fields of quantitative finance\n",
    "\n",
    "* Derivative pricing and hedging\n",
    "* Risk management, regulatory capital and stress testing\n",
    "* Portfolio management\n",
    "* Quantitative strategy/Algo trading\n",
    "* Behavior models\n",
    "\n",
    "We will cover many topics in the first three fields."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Main topics in this class\n",
    "\n",
    "A rough schedule of the class:\n",
    "\n",
    "| Week | Main Topics | Practical Problems | Instructor |\n",
    "| :----: | :-----: | :-----: | :-----: |\n",
    "| 1 | Introduction and Error Analysis | Error and floating point numbers | Y. Li |\n",
    "| 2-3 | Linear Algebra | Portfolio optimization, PCA, least square | Y. Li |\n",
    "| 4-5 | Rootfinding and interpolation | Curve building | Y. Li |\n",
    "| 6 | Derivatives and integration | Hedging, Risk Transformation | Y. Li |\n",
    "| 7-8 | Monte Carlo simulation | Exotic pricing | Y. Li |\n",
    "| 9-10 | Optimization | Model calibration | H. Cheng |\n",
    "| 11-12 | ODE/PDE | Pricing | H. Cheng |\n",
    "| 13 | Entropy and Allocation | Advanced practical problems | Y. Li |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Text book\n",
    "\n",
    "D Bindel and J Goodman, 2009\n",
    "\n",
    "* Principles of scientific computing [online book](http://www.cs.nyu.edu/courses/spring09/G22.2112-001/book/book.pdf)\n",
    "\n",
    "  \n",
    "## References\n",
    "\n",
    "L Anderson and V. Piterbarg, 2010\n",
    "\n",
    "* Interest rate modeling, volume I: Foundations and Vanilla Models \n",
    "* Interest rate modeling, volume II: Vanilla Models\n",
    "  \n",
    "P. Glasserman 2010\n",
    "  \n",
    "* Monte Carlo method in Financial Engineering\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Lecture notes and homework\n",
    "\n",
    "Available online at: http://yadongli.github.io/nyumath2048\n",
    "\n",
    "\n",
    "## Grading\n",
    "* Homework 50%, can be done by groups of two students\n",
    "* Final 50% \n",
    "* Extra credits \n",
    " * extra credit homework problems\n",
    " * for pointing out errors in slides/homework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Software Development Practices\n",
    "\n",
    "Roberto Waltman: In the one and only true way. The object-oriented version of \"Spaghetti code\" is, of course, \"Lasagna code\". "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Modern hardware\n",
    "\n",
    "* Processor speed/throughput continues to grow at an impressive rate\n",
    " * Moore's law has ruled for a long time\n",
    " * recent trends: multi-core, GPU, TPU\n",
    "* Memory capacity continues to grow, but memory speed only grows at a (relatively) slow rate\n",
    "\n",
    "Memory hierarchy\n",
    "\n",
    "| Storage | Bandwidth | Latency | Capacity |\n",
    "| :-----: | --------: | --: | ---: | \n",
    "| L-1 Cache |  98GB/s | 4 cycles | 32K/core |\n",
    "| L-2 Cache |  50GB/s | 10 cycles | 256K/core |\n",
    "| L-3 Cache |  30GB/s | 40 - 75 cycles | 8MB shared |\n",
    "| RAM | 15GB/s | 60-100 ns | ~TB |\n",
    "| SDD | up to 4GB/s | ~0.1ms | $\\infty$ |\n",
    "| Network | up to 12GB/s | much longer | $\\infty$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## On high performance computing\n",
    "\n",
    "* Understand if the problem is computation or bandwidth bound\n",
    " * Most problems in practice are bandwidth bound\n",
    "* Computation bound problems:\n",
    " * caching\n",
    " * vectorization/parallelization/GPU\n",
    " * optimize your \"computational kernel\"\n",
    "* Bandwidth bound problems:\n",
    " * optimize cache and memory access\n",
    " * require highly specialized skills and low level programming (like in C/C++/Fortran)\n",
    "* Premature optimization for speed is the root of all evil\n",
    " * Simplicity, generality and scalability of the code are often sacrificed\n",
    " * Execution speed is often not the most critical factor in most circumstances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Before optimizing for speed\n",
    "\n",
    "* Correctness\n",
    "    - \"those who would give up correctness for a little temporary performance deserve neither correctness nor performance\"\n",
    "    - \"the greatest performance improvement of all is when system goes from not-working to working\"\n",
    "* Modifiability\n",
    "    - \"modification is undesirable but modifiability is paramount\"\n",
    "    - modular design, learn from hardware designers\n",
    "* Simplicity\n",
    "    - no obvious defects is not the same as obviously no defects\n",
    "* Other desirable features:\n",
    "    - Scalability, robustness, Easy support and maintenance etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In practice, what often trumps everything else is: \n",
    "\n",
    "## Time to market !!\n",
    "\n",
    "* First mover has significant advantage: Facebook, Google etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Conventional advice for programming\n",
    "\n",
    "* Choose good names for your variables and functions\n",
    "* Write comments and documents\n",
    "* Write good tests\n",
    "* Keep coding style consistent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "These are good advice, but quite superficial. \n",
    "\n",
    "* all of them have been compromised, even in very successful projects\n",
    "* they are usually not the critical differentiator between successes and failures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Best advice for the real world:\n",
    "\n",
    "by Paul Phillips \"We're doing it all wrong\":\n",
    "\n",
    "(quoting George Orwell's 1984)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### War is peace\n",
    " * Fight your war at the boundary so that you can have peace within\n",
    " * Challenge your clients to reduce the scope of your deliverable\n",
    " * Define a minimal set of interactions to the outside world\n",
    " \n",
    "#### Ignorance is strength\n",
    " * Design dumb modules that only do one thing\n",
    " * Keep it simple and stupid\n",
    " * Separation of concerns, modular\n",
    " \n",
    "#### Freedom is slavery\n",
    " * Must be disciplined in testing, code review, releases etc. \n",
    " * Limit the interface of your module (or external commitments)\n",
    " * Freedom $\\rightarrow$ possibilities $\\rightarrow$ complexity $\\rightarrow$ unmodifiable code $\\rightarrow$ slavery"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Less is more\n",
    "\n",
    "Sun Tzu (孙子): \"the supreme art of quant modeling is to solve problems without coding\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Choose the right programming tools\n",
    "\n",
    "* Programming paradigm: procedural, OOP, functional\n",
    "* Expressivity: fewer lines of code for the same functionality\n",
    "* Ecosystem: library, tools, community, industry support\n",
    "\n",
    "There is no single right answer\n",
    "\n",
    "* Different job calls for different tools\n",
    "* But most of the time, the choice is already made (by someone else)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Expressivity\n",
    "\n",
    "$$ $$\n",
    "\n",
    "<center>\n",
    "<img src=\"img/proglang.png\" width=500 height=500>\n",
    "</center>\n",
    "\n",
    "When you have the opportunity, choose the more expressive language for the job."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Eric Raymond, \"How to Become a Hacker\":\n",
    "\n",
    "\"Lisp is worth learning for the profound enlightenment experience you will have when you finally get it; that experience will make you a better programmer for the rest of your days, even if you never actually use Lisp itself a lot.\"\n",
    "\n",
    "- Clojure is a modern implementation of Lisp\n",
    "- Julia is a very promising new language, heavily influenced by Lisp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## We use Python for this class\n",
    "\n",
    "Python\n",
    "\n",
    "* a powerful dynamic and strongly typed scripting language \n",
    "* extremely popular in scientific computing\n",
    "* strong momentum in quantitative finance\n",
    "* highly productive, low learning curve\n",
    "\n",
    "Python ecosystem\n",
    "\n",
    "* numpy/scipy: high quality matrix and numerical library\n",
    "* matplotlib: powerful visualization library\n",
    "* Jupyter notebook: highly productive interactive working environment\n",
    "* Pandas: powerful data and time series analysis package\n",
    "* Machine learning frameworks: Tensorflow, PyTorch, Scikit-learn etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Python resources\n",
    "\n",
    "* Anaconda: a popular Python distribution with essential packages, best choice\n",
    "* Enthought Canopy: a commercial Python distribution, free for students\n",
    "* Google and Youtube\n",
    "\n",
    "Python version:\n",
    "\n",
    "* ** Python 3.6.x or higher ** : python is not backward compatible, do not use Python 2\n",
    "* With the latest Jupyter notebook, numpy/pandas/scipy etc.\n",
    "* either 32bit or 64bit version works\n",
    "\n",
    "Python setup instructions:\n",
    "\n",
    "* available at: http://yadongli.github.io/nyumath2048"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Errors and floating point computation\n",
    "\n",
    "\n",
    "James Gosling: \"95% of the folks out there are completely clueless about floating-point\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Absolute and relative error\n",
    "\n",
    "If the true value $a$ is approximated by $\\hat{a}$: \n",
    "\n",
    "* Absolute error: $\\hat{a} - a$\n",
    "* Relative error: $\\frac{\\hat{a} - a}{a}$\n",
    "\n",
    "Both measure are useful in practice.\n",
    "\n",
    "* Beware of small $|a|$ when comparing relative errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Sources of error\n",
    "\n",
    "Imprecision in the theory or model itself\n",
    "\n",
    "* George Box: \"essentially, all models are wrong, but some are useful\"\n",
    "\n",
    "Errors in numerical algorithm\n",
    "\n",
    "* Convergence error\n",
    " * Monte Carlo simulation\n",
    " * Numerical integration\n",
    "* Truncation and discretization\n",
    "* Termination of iterations\n",
    "\n",
    "Error due to machine representation\n",
    "\n",
    "* Floating point number\n",
    "* Rounding\n",
    " \n",
    "We cover the last category in this lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## IEEE 754 in a nutshell\n",
    "\n",
    "First published in 1985, the IEEE 754 standard includes:\n",
    "\n",
    "* Floating point number representations \n",
    "  * 16bit, 32bit, 64bit etc\n",
    "  * special representation for NaN, Inf etc\n",
    "* Rounding rules\n",
    "* Basic operations\n",
    "  * arithmetic: add, multiplication, sqrt etc\n",
    "  * conversion: btw formats, to/from string etc\n",
    "* Exception handling\n",
    "  * invalid operation, divide by zero, overflow, underflow, inexact\n",
    "  * the operation returns a value (could be NaN, Inf or 0), and set an error flag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Rounding under IEEE754\n",
    "\n",
    "Default rule: round to the nearest, ties to even digit\n",
    "* Example, round to 4 significant digit(bit): \n",
    " * decimal: $3.12450 \\rightarrow 3.1240$, $435150 \\rightarrow 435200$ \n",
    " * binary: $10.11100 \\rightarrow 11.00$, $1100100 \\rightarrow 1100000 $\n",
    "* Addition and multiplication are commutative with rounding, but not associative\n",
    "* Intermediate calculations could be done at higher precision to minimize numerical error, if supported by the hardware\n",
    "  * Intel CPU's FPU can use 80bit representation for intermediate calculations, even if input/outputs are in 64/32 bit floats\n",
    "  * Early version of Java (before v1.2) enforces intermediate calculation to be in 32/64bits for float/double for binary compatibility - a poor choice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Importance of IEEE 754  \n",
    "\n",
    "IEEE 754 is a tremendously successful standard:\n",
    "* Allows users to consistently reason about floating point computation\n",
    "  * before IEEE 754, floating point number calculation behaves differently across different platforms\n",
    "* IEEE 754 made compromise for practicality\n",
    "  * many of them do not make strict mathematical sense\n",
    "  * but they are necessary evils for the greater goods\n",
    "* Establishes a clear interface between hardware designers and software developers\n",
    "* Its main designer, William Kahan won the Turing award in 1989"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The IEEE 754 standard is widely adopted in modern hardware and software,\n",
    "\n",
    "* The vast majority of general purpose CPUs are IEEE 754 compliant\n",
    "* Software is more complicated  \n",
    "  * Compiled language like C++ depends on the hardware's implementation.\n",
    "  * Some software implemented their own special behaviours, e.g., the Java strictfp flag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## IEEE 754 floating point number format\n",
    "\n",
    "* IEEE floating point number standard\n",
    "\n",
    "$$\\underbrace{*}_{\\text{sign bit}}\\;\\underbrace{********}_{\\text{exponent}}\\overbrace{\\;1.\\;}^{\\text{implicit}}\\;\\underbrace{***********************}_{\\text{fraction}}\n",
    "$$\n",
    "\n",
    "$$(-1)^{\\text{sign}}(\\text{1 + fraction})2^{\\text{exponent} - p}$$\n",
    "\n",
    "* 32 bit vs. 64 bit format\n",
    "\n",
    "| Format | Sign | Fraction | Exponent | $p$ | (Exponent-p) range |\n",
    "| :----: | ---: | ----------: | -------: | --------: | :--: |\n",
    "| IEEE 754 - 32 bit | 1    |      23     |   8      |  127 | -126 to 127 |\n",
    "| IEEE 754 - 64 bit | 1    |      52     |   11     | 1023 | -1022 to 1023 |\n",
    "\n",
    "* an implicit leading bit of 1 and the binary point before fraction bits\n",
    "* 0 and max exponent reserved for special interpretations: 0, NaN, Inf etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Examples of floating numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "# need to do \"pip install bitstring\"\n",
    "%pylab inline\n",
    "\n",
    "import bitstring\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import fmt\n",
    "import math, sys\n",
    "\n",
    "pd.set_option(\"display.max_colwidth\", 80)\n",
    "\n",
    "def parts2f(sign, e, sig) :\n",
    "    return ((-1)**sign*sig*2.**e)\n",
    "\n",
    "def f2parts(v, flen) :\n",
    "    ps = {32: 9, 64 : 12}\n",
    "    ed = ps[flen]\n",
    "    eoff = 2**(ed - 2) - 1\n",
    "\n",
    "    f1 = bitstring.BitArray(float=v, length=flen)\n",
    "    signb = f1[:1]\n",
    "    sigb = f1[ed:]\n",
    "    eb = f1[1:ed]\n",
    "    sig = (1. + 1.*sigb.uint/(2**(flen-ed))) if eb.uint > 0 else 2.*sigb.uint/(2**(flen-ed))\n",
    "    e = eb.uint - eoff\n",
    "    \n",
    "    bins = np.array([signb.bin, eb.bin, sigb.bin])\n",
    "    vals = np.array([1-signb.uint*2, e, sig])\n",
    "    \n",
    "    fmt.displayDF(pd.DataFrame(np.array([bins, vals]).T, columns=['Bits', 'Decimal'], \n",
    "                               index=[\"Sign\", \"Exp\", \"Fraction\"]).T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "f2parts takes a floating number and its bit length (32 or 64) and shows the three parts. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sign</th>\n",
       "      <th>Exp</th>\n",
       "      <th>Fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bits</th>\n",
       "      <td>0</td>\n",
       "      <td>10000010</td>\n",
       "      <td>01010110011001100110011</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.337499976158142</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(10.7, 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sign</th>\n",
       "      <th>Exp</th>\n",
       "      <th>Fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bits</th>\n",
       "      <td>1</td>\n",
       "      <td>10000001101</td>\n",
       "      <td>0110111001010101000000000000000000000000000000000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>1.4309844970703125</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(-23445.25, 64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Range and precision\n",
    "\n",
    "* The range of floating point number are usually adequate in practice\n",
    " * the whole universe only have $10^{80}$ protons ...\n",
    "* The machine precision is the more commonly encountered limitation\n",
    " * machine precision $\\epsilon_m$ is the maximum $\\epsilon$ that makes $1 + \\epsilon = 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Min</th>\n",
       "      <th>Max</th>\n",
       "      <th>Machine Precision</th>\n",
       "      <th># of Significant Digits</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>32 bit</th>\n",
       "      <td>1.1755e-38</td>\n",
       "      <td>3.4028e+38</td>\n",
       "      <td>5.9605e-08</td>\n",
       "      <td>7.2247</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>64 bit</th>\n",
       "      <td>2.2251e-308</td>\n",
       "      <td>1.7977e+308</td>\n",
       "      <td>1.1102e-16</td>\n",
       "      <td>15.955</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prec32 = 2**(-24)\n",
    "prec64 = 2**(-53)\n",
    "\n",
    "f32 = [parts2f(0, -126, 1), parts2f(0, 254-127, (2.-prec32*2)), prec32, -np.log10(prec32)]\n",
    "f64 = [parts2f(0, -1022, 1), parts2f(0, 2046-1023, (2.-prec64*2)), prec64, -np.log10(prec64)]\n",
    "\n",
    "fmt.displayDF(pd.DataFrame(np.array([f32, f64]), index=['32 bit', '64 bit'], \n",
    "                           columns=[\"Min\", \"Max\", \"Machine Precision\", \"# of Significant Digits\"]), \"5g\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Examples of floating point representations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sign</th>\n",
       "      <th>Exp</th>\n",
       "      <th>Fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bits</th>\n",
       "      <td>1</td>\n",
       "      <td>00000000</td>\n",
       "      <td>00000000000000000000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>-127.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(-0., 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sign</th>\n",
       "      <th>Exp</th>\n",
       "      <th>Fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bits</th>\n",
       "      <td>0</td>\n",
       "      <td>11111111</td>\n",
       "      <td>10000000000000000000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>1.0</td>\n",
       "      <td>128.0</td>\n",
       "      <td>1.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(np.NaN, 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sign</th>\n",
       "      <th>Exp</th>\n",
       "      <th>Fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bits</th>\n",
       "      <td>1</td>\n",
       "      <td>11111111</td>\n",
       "      <td>00000000000000000000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>128.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(-np.Inf, 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sign</th>\n",
       "      <th>Exp</th>\n",
       "      <th>Fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bits</th>\n",
       "      <td>1</td>\n",
       "      <td>00000000</td>\n",
       "      <td>00000000000000001000111</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>-127.0</td>\n",
       "      <td>1.6927719116210938e-05</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(-1e-43, 32);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Floating point computing is NOT exact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b == c? False !\n",
      "b - c = -1.110223e-16\n"
     ]
    }
   ],
   "source": [
    "a = 1./3\n",
    "b = a + a + 1. - 1.\n",
    "c = 2*a\n",
    "print(\"b == c? {} !\".format(b == c))\n",
    "print(\"b - c = {:e}\".format(b - c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sign</th>\n",
       "      <th>Exp</th>\n",
       "      <th>Fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bits</th>\n",
       "      <td>0</td>\n",
       "      <td>01111111110</td>\n",
       "      <td>0101010101010101010101010101010101010101010101010100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>1.333333333333333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(b, 64);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><font size=2><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sign</th>\n",
       "      <th>Exp</th>\n",
       "      <th>Fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Bits</th>\n",
       "      <td>0</td>\n",
       "      <td>01111111110</td>\n",
       "      <td>0101010101010101010101010101010101010101010101010101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>1.3333333333333333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(c, 64);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Can't handle numbers too large or small"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inverse of 1.000000e-308 is 1.000000e+308\n"
     ]
    }
   ],
   "source": [
    "x = 1e-308\n",
    "print(\"inverse of {:e} is {:e}\".format(x, 1/x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inverse of 1.000000e-309 is inf\n"
     ]
    }
   ],
   "source": [
    "x = 1e-309\n",
    "print(\"inverse of {:e} is {:e}\".format(x, 1/x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Limited precision\n",
    "\n",
    "* Small number (in relative sense) is lost in addition and subtraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n"
     ]
    }
   ],
   "source": [
    "print(1. - 1e-16 == 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "print(1. + 1e-16 == 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integer arithmetic:\n",
      "a =  265252859812191058636308480000000\n",
      "a - 4294967295 == a ? False\n",
      "\n",
      "Floating point arithmetic:\n",
      "b =  2.6525285981219107e+32\n",
      "b - 4294967295 == b ? True\n"
     ]
    }
   ],
   "source": [
    "largeint = 2**32 - 1\n",
    "\n",
    "print(\"Integer arithmetic:\")\n",
    "a = math.factorial(30) # in Python's bignum\n",
    "print(\"a = \", a)\n",
    "print(\"a - {:d} == a ?\".format(largeint), (a - largeint) == a)\n",
    "\n",
    "print(\"\\nFloating point arithmetic:\")\n",
    "b = float(math.factorial(30))\n",
    "print(\"b = \", b)\n",
    "print(\"b - {:d} == b ?\".format(largeint), (b - largeint) == b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Subtle errors\n",
    "\n",
    "The limited precision can cause subtle errors that are puzzling even for experienced programmers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a =  [0. 1.]\n",
      "a/3 =  [0.         0.33333333]\n"
     ]
    }
   ],
   "source": [
    "a = np.arange(0, 3. - 1., 1.)\n",
    "print(\"a = \", a)\n",
    "print(\"a/3 = \", a/3.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b =  [0.         0.33333333 0.66666667]\n"
     ]
    }
   ],
   "source": [
    "b = np.arange(0, 1 - 1./3., 1./3.)\n",
    "print(\"b = \", b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Avoid equality test\n",
    "\n",
    "Equality test in floating point number is **always** a bad idea:\n",
    "* Mathematical equality often fails floating point equality test because of rounding errors \n",
    "* Often produce off-by-one number of loops when used to test end loop conditions\n",
    "\n",
    "It should become your instinct to:\n",
    "* always use error bound when comparing floating numbers, such as ```numpy.allclose()``` \n",
    "* try to use integer comparison as the end condition for loops\n",
    "  * integer representation and arithmetic in IEEE 754 is exact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Catastrophic cancellation\n",
    "\n",
    "Dramatic loss of precision can happen when very similar numbers are subtracted\n",
    "\n",
    "* Value the following function around 0: \n",
    "\n",
    "$$f(x) = \\frac{1}{x^2}(1-\\cos(x))$$\n",
    "\n",
    " * Mathematically: $\\lim_{x\\rightarrow 0}f(x) = \\frac{1}{2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The straight forward implementation failed spectacularly for small $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEJCAYAAABv6GdPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4nNWZ+P3vmaZebdmyJdly7102BFMcMM04EOqLIWSBBNgQloRlQ1kIyQbY3YQkEFh+m5BASGFxTHECxjSDKQZ33LvcJNmSLKtbbdp5/5h5xiNpJM2MpsjS/bkurosZPfPM8cwz99xzyn2U1hohhBD9iyneDRBCCBF5EtyFEKIfkuAuhBD9kAR3IYTohyS4CyFEPyTBXQgh+iEJ7kII0Q9JcBdCiH5IgrsQQvRDlng98eDBg3VhYWG8nl4IIc5ImzdvPqm1zunpuLgF98LCQjZt2hSvpxdCiDOSUupoMMdJt4wQQvRDEtyFEKIfkuAuhBD9UNz63IUQsedwOCgrK6O1tTXeTRE9SExMJD8/H6vVGtbjJbgLMYCUlZWRlpZGYWEhSql4N0d0QWtNdXU1ZWVljBo1KqxzSLeMEANIa2srgwYNksDexymlGDRoUK9+YUlwF2KAkcB+Zujt+yTBXUTVVyW17DxWH+9mDChltc18vLcy3s0QcSbBXUTVo8t38ov398W7GQPKX9Yd5c4/b8bhcse7KQEppbj//vt9t3/5y1/y05/+NKZt2LRpE/fee29Yj12wYEHABZiff/45U6ZMYebMmbS0tFBeXs7ixYu7PdeKFSt47LHHwmpHTyS4i6jRWnO0uolWhyveTRlQ2hxunG5NeV3fnBGTkJDAm2++ycmTJ+Py/E6nk6KiIp599tmInveVV17h4YcfZuvWrSQlJfHrX/+aO+64o9vHXHHFFbz99ts0NzdHtC0gwV1EUXWTnSa7C2cfzSD7KyNjP1rTFOeWBGaxWLjzzjt5+umnO/3t1ltv5fXXX/fdTk1NBeCTTz7hggsu4KqrrmL06NE89NBDvPLKK8ybN49p06Zx8OBBAKqqqrj22muZO3cuc+fO5YsvvgDgpz/9Kbfccgvz58/nlltu4ZNPPvFl1adOneK2225j2rRpTJ8+nTfeeAOA733vexQVFTFlyhR+8pOfdPtv+sMf/sCyZcv48Y9/zM033wzAG2+8wWWXXQbA008/ze233w7Ajh07mDp1Ks3NzSilWLBgAStWrAj79eyKTIUUUVNS48lGHC4d55YMLEZwN17/rvzH27vYfbwhos89eXg6P/nGlB6P+/73v8/06dN54IEHgj73tm3b2LNnD9nZ2YwePZrvfve7bNiwgd/85jc899xzPPPMM/zgBz/gvvvu49xzz6WkpIRLL72UPXv2ALB7927WrFlDUlISn3zyie+8jz/+OBkZGezYsQOA2tpaAJ588kmys7NxuVxcdNFFbN++nenTpwds23e/+13WrFnD4sWLue666zh8+DBZWVkkJCQA8IMf/IAFCxawfPlynnzySX73u9+RnJwMQFFREZ9//jk33HBD0K9FMCS4i6gp9QV3ydxjyen9Mu0puMdTeno63/72t3n22WdJSkoK6jFz585l2LBhAIwZM4ZLLrkEgGnTprF69WoAVq1axe7du32PaWho4NSpUwBceeWVAZ9r1apVLF261Hc7KysLgGXLlvHCCy/gdDopLy9n9+7dXQb3jsrLy8nJOV240WQy8fLLLzN9+nTuuusu5s+f7/vbkCFDOH78eFDnDYUEdxE1R6sluMeD3cjcq7sP7sFk2NH0wx/+kNmzZ3Pbbbf57rNYLLjdnva73W7sdrvvb0YWDJ5gadw2mUw4nU7fY9atW0diYmKn50tJSQm6bYcPH+aXv/wlGzduJCsri1tvvTWkOedJSUmdjj9w4ACpqamdAnlra2vQX3ChkD53ETVG5uh0S7dMLJ0JmTtAdnY2N9xwAy+++KLvvsLCQjZv3gzAW2+9hcPhCOmcl1xyCc8995zv9tatW3t8zMUXX8zzzz/vu11bW0tDQwMpKSlkZGRQWVnJu+++G1I7xo8fz5EjR3y36+vruffee/nss8+orq5uN66wf/9+pk6dGtL5gyHBXUSNkTk6nJK5x5LDL3PXum9/sd5///3tZs3ccccdfPrpp8yYMYO1a9eGlG0DPPvss2zatInp06czefJkfvvb3/b4mEcffZTa2lqmTp3KjBkzWL16NTNmzGDWrFlMnDiRm266qV03SjBSUlIYM2YMxcXFANx33318//vfZ/z48bz44os89NBDnDhxAoDVq1dzxRVXhHT+YKh4vflFRUVaNuvo387+z4+oaGglJy2BjY8sjHdzBoxvv7SBz/ZXAbDlxxeTlWLz/W3Pnj1MmjQpXk0bUJYvX87mzZt54oknujymsrKSm266iY8++ijg3wO9X0qpzVrrop6eX/rcRVS0OlxUNHj6HKXPPbb8fymV1DS3C+4idq6++mqqq6u7PaakpIRf/epXUXl+Ce4iKspqPV0yg1NttDokuMeS0+1mcKqNk6fsHK1pZkZBZrybNGB997vf7fbvc+fOjdpzS5+7iApjMG90Tqpv9oaIDbtLM3qwZ/FPaR8fVBXRI8FdRIUxmDomJ0VWqMaY0+UmLdFCTlpCj9MhRf8lwV1ExdGaZpJtZoamJ+LW4JLpkDHjcLmxmk2MyE7usyUIRPRJcBdRUVrTzIjsZKxmzyUmg6qx43RpLGbFyOxkSmta4t0cEScS3EVUHK32BHebBPeYs7vc2MwmCrKTOV7fQpuz71TlrK6uZubMmcycOZPc3Fzy8vJ8t/1Xo/ZGY2MjgwYN8pUdMCxevNhXFCyQVatW8c1vfjMibegLJLiLiNNaU+LN3C1mz24yTikeFjNG5j4iOxmt4Vht38neBw0axNatW9m6dSv//M//zH333ee7bbN5pmxqrX0lCMKRlpbGRRddxD/+8Q/ffbW1taxbty4qi4X6KgnuIuKqGttoc7oZOUi6ZeLB6HMfOchTdbCvlyEAKC4uZvLkydx8881MmTKF0tJSMjNPT+FcunSpb1phZWUl11xzDUVFRcybN49169Z1Ot+SJUvaFQN74403uOKKK0hMTGTdunV87WtfY9asWcyfP58DBw50evyjjz7KM88847s9ceJEysrKAPjTn/7EvHnzmDlzJnfffXevvoiiSea5i4g76g0mBdnJVBoLmWRANWb8B1Shm+D+7kNQsSOyT547DS7/77AeunfvXv785z9TVFTkKwQWyL333ssDDzzA2WefzZEjR1i8eDE7d+5sd8yiRYu46667qK2tJSsri6VLl/Jv//ZvAEyaNInPP/8ci8XCe++9x6OPPsrf/va3oNq4c+dOli9fzpdffumrS7906VJuuummsP7N0STBXUScMf1uRHYyNU2eflSpLxM7DpfGalbkpCWQaDWdMdMhx4wZQ1FRj6vqWbVqFfv2nd66sba2lpaWlnaVFRMSErjiiit48803Wbx4Mbt27WLhQk8JjLq6Or797W/7NvgIxapVq9i4caOvnS0tLRQUFIR8nliQ4C4i7mhNM0pBflYyO72bQTj76E/X/sjpdmMxm1BKeadDdhHcw8ywo8W/SJjJZGpX9My/fK7Wmg0bNvj66LuyZMkSnnrqKVpaWrj66quxWDzh7pFHHuHSSy/l7rvvpri42Ldbkj//0sP+z6+15vbbb+fxxx8P7x8ZQ9LnLiKutKaZ4RlJ2CwmbN4BVbtTumViQWvtzdw9H+0R2cln5CpVk8lEVlYWBw4cwO12s3z5ct/fFi5c2K5Eb1dlfS+66CJ27drFb3/7W5YsWeK7v76+nry8PABefvnlgI/1Lz28YcMGSktLfc+9bNkyXyXL6upqSkpKwv+HRlFQwV0pdZlSap9Sqlgp9VCAv49QSq1WSm1RSm1XSi2KfFPFmaKkppmCbM9PZIvJc4lJ5h4bRu18q8nzpVqQnUxJTd8v/RvIz3/+cy699FLOOecc8vPzffc///zzfPHFF76yvr///e8DPt5sNnPNNdfQ0NDAueee67v/wQcf5Ec/+hGzZ8/u8nW5/vrrqaysZOrUqbzwwguMHj0a8Oz69JOf/ISFCxcyffp0LrnkEiorKyP4r44grXW3/wFm4CAwGrAB24DJHY55Afie9/8nA0d6Ou+cOXO06J/mPP6hfuC1bVprrT/Zd0KPfHCF3nSkOs6tGhia2hx65IMr9P9+Uqy11vqPaw7pkQ+u0CcaWrXWWu/evTuezRMhCvR+AZt0D/FVax1U5j4PKNZaH9Ja24GlwFUdvyOAdO//ZwCR3xBQnBGa7U5OnmpjhHcanpFByibZsWG8zhbv6z7iDJoOKSIrmOCeB5T63S7z3ufvp8C3lFJlwErgXyLSOnHGMZa7F3in4Vkt3m4ZCe4xYRRps1lO97mDVIcciCI1oLoEeFlrnQ8sAv6ilOp0bqXUnUqpTUqpTVVVVRF6atGXHK32FKoa6Q0qFl/mLn3usXA6c/d8/PKzPO/DUb/pkPoM7H8fiHr7PgUT3I8B/hM58733+fsOsMzboLVAIjC444m01i9orYu01kU5OTnhtVj0acbPfyNjlBWqsWW8zlbvLKVEq5nc9ETf+5KYmEh1dbUE+D5Oa011dTWJiYlhnyOYee4bgXFKqVF4gvqNQMflWCXARcDLSqlJeIK7pOYDUElNM2mJFjKTrYB/cJdgEgung/vpvG3EoGRKvKV/8/PzKSsrQ345932JiYntZgmFqsfgrrV2KqXuAd7HM3PmJa31LqXUz/CM2r4F3A/8Xil1H57B1Vu1pAYDklEwTClP5mhkkDIVMjZ8UyH9g3t2Mp8f8ARzq9XKqFGj4tI2EVtBrVDVWq/EM1Dqf99jfv+/G5gf2aaJM1FJTTMThqb5bhtBxi7lB2LCeJ2NapzgCe6VDW20OlwkWs3xapqIMVmhKiLG5daU1bT4+tvhdHB3SuGwmDBeZ5tf5n4mVYcUkSPBXURMSU0zdpebMUNSffcZGaQMqMaG8Tr7Z+5jcjzvR/GJUwEfI/onCe4iYvZVNAIE7JaRAdXYCDSgOnZIKkqdfn/EwCDBXUTM/kpP8Bg39HTmbpXMPaaML1GrX+aeaDVTOCjF9/6IgUGCu4iYfZWNjMhOJtl2epze1+cuwT0mnAEyd4DxQ1PZJ8F9QJHgLiJmf0Uj4/26ZOD0ClW7dMvEhK/P3dT+oz1haBpHTjbR6ug7m2WL6JLgLiKizeni8MkmJuSmtrtfKYXVrCRzjxGjW8ZmUe3un5CbjlvLoOpAIsFdRMShqiacbt0pcwdPFil97rHRZebu/dKVfveBQ4K7iAgjaEzMTe/0N6tZyWyZGDGqbxrVOA0jB6VgM5uk330AkeAuImJfRSMWk2LU4JROf7OaJXOPFbsxoGpq3y1jNZsYnZPCfpkOOWBIcBcRsb+ykdE5Kb464v6sZpPUc4+RrmbLAEzITWN/pfS5DxQS3EVE7KvsPFPGYDErydxjxFfP3aw6/W380DSO1bXQ0OqIdbNEHEhwF73W1OaktKal3cpUfzazCYfUlokJh7vrzH1iruf9OSD97gOCBHfRawe80+vG53aTuUtVyJhwODuX/DUYv6z2VUjXzEAgwV302v4ANWX8Wc0mqeceI063G5MCs6lzt0xeZhIpNrNMhxwgJLiLXttb0Uii1dSu1K8/i9kkK1RjxO5yYwmQtQOYTIpxQ9OkgNgAIcFd9Np+72CqKUC2CGCTFaox43TpdrXcO5owNE0y9wFCgrvote5myoCsUI0lh8sdcKaMYUJuGtVNdqoa22LYKhEPEtxFr9R4A0VX/e1gTIWUbplYcLh0p9ID/iZ4B70le+//JLiLXjGCRFczZcA7FVIy95hwuNzYusncT8+YkeDe30lwF71iBPeeMndZoRobzm4GVAEGp9rITrFJ5j4ASHAXvbKvopH0RAtD0xO6PEZqy8SOw6Xb7cLUkVJKNu4YICS4i17ZV9HIxNx0lOo6oFjNJt/KSRFdDpc74AImfxOGprG/ohGt5ddUfybBXYRNa+2ZKdNhg46OrGblWzkpoiuo4J6bTpPdRVltS4xaJeJBgrsIW0VDK42tzm7728GziElWqMaG0627nQoJsnHHQCHBXYRt9/EGwJMJdsdmNmGX2jIxYXf2nLmP834ZG++f6J8kuIuwbS2tw2xSTM3rPrhbTAqnVIWMCae7+wFVgPREK2NyUthaWhejVol4kOAuwra1tI4JQ9NItlm6Pc5qkdkysRJMnzvArBFZbC2tk0HVfkyCuwiL263ZWlLHrBGZPR5rNXlWqEogib6eVqgaZhZkUt1kp7RGBlX7KwnuIiwHq07R2OZkZkEQwd2bSUrXTPQ5XG5slu67ZQDfl/KW0tpoN0nEiQR3EZYt3v7aWSOyejzWWDEpq1Sjz+lyB5W5TxiaRpLVzJYS6XfvryS4i7BsKakjPdHC6MEpPR5rDPDZpd896jwrVHv+WFvMJqblZ/i+pEX/I8FdhGVraR0zCjK7rOHuz9ctI8E96jwDqj2/J+DpmtlzvIE2pyvKrRLxIMG9n6psaOXIyaaonLupzcm+ioagumTgdHCXsr/RF+xsGYBZBZnYXW52RWm+e1VjG4eqZL/WeAnqKlBKXaaU2qeUKlZKPdTFMTcopXYrpXYppf4vss0Uofrh0q3c8+pXUTn39rJ63NoTHIJhrJiU6ZDR53T1vELVYHw5R6vf/an393L9b9fKL7Y46TG4K6XMwPPA5cBkYIlSanKHY8YBDwPztdZTgB9Goa0iSFWNbaw7XE1lQ3R22zEWvwQzUwbwbfsmwT367C53t9vs+RuansjwjMSoLWaqaGijusnO+sM1UTm/6F4wV8E8oFhrfUhrbQeWAld1OOYO4HmtdS2A1vpEZJspQvH+rgq0hvpmR1Tmlm8pqWXU4BSyUmxBHW9kkjIVMvqCqS3jb+aITLaURGc6ZH2zHYCVO8qjcn7RvWCCex5Q6ne7zHufv/HAeKXUF0qpdUqpyyLVQBG6d3d6Pkx2l5tme2QHy7TWbCmtCzprh9N97lJfJrrcbo3LHdxsGcOsgizKaluisqdqbbMD8CQbLvlij7lIDahagHHAAmAJ8HulVKdPv1LqTqXUJqXUpqqqqgg9tfBXfaqNdYdqyE1PBKCuxRHR8x+vb6WqsS2olakGq2TuMWHUzA8luM/0vo/R6Jqpa7aTm57IyVN2Nh6RrplYC+YqOAYU+N3O997nrwx4S2vt0FofBvbjCfbtaK1f0FoXaa2LcnJywm2z6MaHuytxuTU3zPW8ZbVN9oie3/gJP6sguJkyIFMhY8VYJBbsVEiAqcMzsJhUxLtmnC43Da1Orpw5nESriXelaybmggnuG4FxSqlRSikbcCPwVodj/o4na0cpNRhPN82hCLZTBGnlzgpGDkrmnDGDAKiPcOa+paSOBIuJicO6r+Huz1gxKYuYossYsA5mhaohyWZm0rD0iM+YaWh1AjA8I5EF44fw7s4K3PLLLaZ6vAq01k7gHuB9YA+wTGu9Syn1M6XUld7D3geqlVK7gdXAj7TW1dFqdH+zo6yeY3W9L+BU12zny+KTXD51GFnJnsHO2ubIZu5bS+uYlpcR0k9/o9aJlB+ILmMdgdUSWm/rzIJMtpfVRbRf3LjuMpNtXD4tlxONbXwVgV8HDa0Oviw+2evzDARBXQVa65Va6/Fa6zFa6ye99z2mtX7L+/9aa/2vWuvJWutpWuul0Wx0f9LQ6uDGF9by1Ht7e32uD3dX4nRrFk3LJSvZCkBdc+Qyd7vTzY5j9SH1t8PpTFKmQkaX8fpag1g17G/WiEya7C4OnIjczkzGdZeZbOXCiUOwWUys3FHR6/O+ur6Em/6wnp3H6nt9rv5OVqjG2bKNpTTZXRytae71ud7dWUFeZhLT8jLI8AX3yGXue8obsDvdzAyhvx1khWqsnO5zDz1zh8guZqrzy9zTEq2cPy6Hd3eW97prxvic/PGLI71tYr8nwb2XHvvHTl5aczisx7rcmj+tPQLQ682KG1odrDlwkkXTclFKkWAxk2wzRzRz9w2mhpi5W2WFakwYYxqhzHMHGDU4hYwka0QHVY3rzvgFuWhaLuX1rWwr690XiPE5eXvb8bCnb67ed4LbX97Y72vqSHDvhY1Havjz2qO8uaUsrMd/tKeS0poWpgxPp6qxjVZH+Bfbx3tOYHe5uWzqMN99mUlW31zjSPjyYDV5mUkMy0gM6XGn67lLcI8m4/UNdoWqQSnF3MJsvjxYHbFFb74+9yTP2M9Fk4ZiNSve3dm7rpmy2mYmDUvH7nLzf+tLwjrHyu3lfLz3BMs2lvZ88BlMgnsv/OqDfQDsrzwV1jS/P35xhLzMJG49pxCgV4OqK3eUk5ue2K7eS2ayjfqWyHTLOF1u1h6s5rxxg1EqtMzQV1vGKd0y0WS8vpYQgzvAeeMGU1bbwtHq3ncPgmeWlklBWqJnC8aMJCvzxw5m5Y7ysL9AtNYcq23hvHGDWTAhh7+uPxrWwri9FZ6xhec+Lu5VQtXXSXAP05fFJ1l3qIaZBZnYnW6OVIdWgXFPeQNrD1Vzy9dGUuitiR5u10xds51P9ldx+bTcdiV4M5Mjl7lvK6unsc3JueMGh/xYX20Zydyj6vQiptC+fAHf+/p5hGai1DbbyUiytrseF00bRlltS9izZqpOtdHmdJOflcRt80dR1djGOzuOh3QOp8vN/spGZhZkcqKxjb+uOxpWW84EEtzDoLXmVx/uJzc9kZ98w1NDbU95aDMNXv7iCIlWEzfOLSAvMwmAY2EG97e2HcfudHPdnPx292cl2yI2oLrmwEmUgvljQg/uRibpkPIDUWW8vqEOqAKMHpzC8IxE1hyIzMrxumaHbzquYdG0YSRZzby2KbxuTOPzkZeZxPnjBjMmJ4U/fnEkpF8CR6qbaXO6ueXskZwzZhC//fQgzXZnWO3p6/pNcP/z2iNc/OtPQxokeXvbcT7aUxnyc326v4rNR2u558KxTB6ejsWk2FsRfE3smiY7f996jGtm55OZbGNoeiIWk6KsNryfxMs2lTJ5WDpThme0uz8j2RqxAdU1xVVMHZ4RdLEwf1J+IDaM1zec4K6U4txxg/nyYHVE5rvXNTt8M7YMqQkWFk0bxort5WEFVOOXbX5WMkopbp0/iu1l9SH9EjA+pxOHpXH/JeM5ecrOn74MPXvfeayeP3x+KKQvlrv+somH39wR8nOFq08G9/L6Fu7729ag+6CLT5ziiXf2cODEKdYdCq6GxcGqU/zrsq08uXJPSG3TWvPrD/eTn5XEDUUFJFjMjMlJZW8ImfurG0poc7q5zdvXbjYphmcmhdUts/t4AzuPNXBDUX6nv2UlW6lr6X1lyFNtTraU1IXVJQN+hcNktkxUhTtbxnDuuBwaW51s7+WMFoC6FnunzB3ghqJ8TrU5eTeMOe/G5yMvy/NL99rZeaQnWngphGmRe8sbMZsUY4ekMmdkNgsm5PC7zw7S2BpaEvS/nx7kiXf28MHu4JLDyoZW3t9VyasbSvhgV3D/9ha7ix+9ti3suj99Lri73JofLt3K8i3H+GMQUwxdbs0Dr28j2WYmyWrmw909v3Baax77x04cLs2hqiZONLYG3b4Pd1eyvayeey8ch827EnDisDTfIE1PHC43f1l7lPPGDWbc0NNL+POzksLK3F/bXIrNbOKqmR0LdXpmKrjcmsa23v3sXHewGqdbc97Y3gV3WaEaXcbrG+psGcN8b8mKNQd63+9e2+QgM8na6f55o7IpHJTMa5tDn6lSVttMVrKV1ATPIG2yzcKSeSN4b2cFx4NMBPdWNDAmJ4UEixmAf714PHXNDl5acyTodmitWe9NIv/jrV1B/QpZ5e0hyMtM4pG/7wyqu3TljnJe21zG3X/dTH0Yv8D7XHD/f6uLWX+4huEZiSzfcqzH0fA/fnGYr0rq+Ok3pnD++MGs2n2ix0z17e3lfFFczbWzPdnuhiA3E3C7NU+vOkDhoGSumX06mE7MTedYXUtQdVyWf3WMioZWvnPuqHb3e4J7aJm73enm71uOcfHkoQG7SzK9P4vDuTD8rSk+SaLVxJzC0BYvGcwmhVIyzz3aHL3M3AelJjBleHpEBlXrWxxkBsjclVJcNyefdYdqOBriJISy2hbys5Lb3fftcwpRwAufBVfKak95IxNz0323p+dncvHkofzh80NBj08dOtnEyVNtXDM7j+P1rTz3cXGPj/lwdyUjspP53S1zqG2y87MVu3t8zLJNpQxOTeBEYxsPL98e8i/wPhXcNx+t5ZmPDnDljOE8cfVUqpvsfLy3630/jpxs4pcf7GPhpCFcNXM4l0zOpaKhlR3dLE1ubHXwxIrdTMvL4D+vmUqKzez7Fu7Jm1uOsae8gfsuHt9uuplRRGtfD9m70+Xmf1YXMz0/gwvGt6+KmZeZzInGtpDGDD7aU0lts4PrA3TJAL4PV2/ry3x+oIp5owb5sp1wWM0mWaEaZb7yA2Fm7uCZNbOlpJamXvzaszvdnGpz+pKLjq6dk49S8Prm0AZWj9W1+CYfGPIyk7hmdh6vbijp8Rd4Q6uDY3UtnYre3X/JeE7ZnUEFacAXL77/9bFcNyef3392iOJuSjecanPyZXE1l0weytS8DL63YAxvfnWMj/d23aVztLqJ9YdruG1+IfdfMoGVOyr4W4jz8vtMcG9odfCDpVsYlpHIE1dP5fxxOQxJS+D1Ln6+ud2aB97Yjs1s4smrp6GU4sKJQzCbFB/s6vpFe/rDA1SdauOJb04lwWJmTmE26w/3XOOsqc3JU+/vZUZBJt+YPrzd3yZ5M4GeBlXf2nackppm7vn62E5zxfO9/YjH64LvIlq2qZTc9ETOGxe4fHIk6suU17dwsKop7C4Zg9WkJHOPMl/hsBCqQnZ03tgcHC4d1GeiK8Yv2KwugvuwjCTOH5fDG5vLgh681VpTVtvs+5z4+/7Xx+J0a37fQ/ZuJF+T/DJ38Pzy/v+KCvjTl0eC2tB7/eFqBqcmMHpwCg9fPpGUBAs//vuuLjPrT/dVYXe5uXjyUADuuXAs44em8u9v7qShi77+1zeXYVLzPBH8AAAgAElEQVRwzew87jp/NPPHDuI/3t5N8YngNxzvE8Fda82jy3dSXt/Kb26cRXqiFYvZxLVz8lm9r4oTDZ0D3l/WHWXD4Rp+vHgyQ70bU2Sl2CgamcWHXQxy7D7ewMtfHuameSOY4V3sc9aobPZXnqKmh7rnv/v0IJUNbTy2eHK7ubsAQ9MTyEy2djsd0uXW/M/qYibmpvneZH/GRRtsv3tlQyuf7q/i2jl5mLsoFGVkTr3J3D/39r+GO5hqsFpMUs89yozX12oJr1sGoKgwiwSLyfe+h8Po3sgI0C1juL4on+P1rXwRZBdQdZOdVoc7YHAfOSiFq2YM56/rSqg+1XVJgr3lp2fKdHT/JRNItJr5z5XdF/Az+tvPGp2NUopBqQn86NIJrD1UzVvbAs+5/3B3BVnJVuaM9HRrJljMPHXdDE40tvLkis4TOlxuzeubyzh/fA7DMpIwmRS/vmEmiVYT9766pdv2+esTwf21zWW8te04P7xonO8FALh+Tj4ut+bNLe33BjlQ2ch/rtzDggk5neZ2Xzx5KPsqGynpsNLO7db8+B87yUq28cClE333nzUqG+i+3/14XQsvfH6Ib8wY3q59BqUUE3PTus3cV+4o51BVE/deNC7gCs/8bE9fYrD97m98VYZbw/VzCro8xuiW6U1N9zUHTjI4NYGJucHXbw/EYjJhl26ZqAqnnntHiVYz80Zl92pQta6HzB08n9PMZCuvBdk14z8NMpC7vz6WVqeLF7uZhLGnopGMJKtvlzJ/OWkJ3P31MazaU9ntF05pTQsVDa2c7Y0bAEvmjWBGfgZPvLOn02fN4XLz8d4TXDhxaLuu3BkFmdx1wRj+tqmU9zvMnllTfJLy+tZ2n+2h6Yk8dd0MdpcHP+U67sF957F6Hv37Ts4ZM4i7vz623d9G56RSNDKL1zaV+n7ytDpc/MurW0hNsPCL66Z3CpSXTM4F4IMOs2b+tPYIm4/W8u+LJrWbfzs9P5NEq6nbn6G/eG8vWsODl03o8piJuensq2gMWPXO7dY89/EBxg1J5bIpuQEfPzQtAXOQc9211ry2qcwz88C7ujWQjKTedcu43Zovik9y7thBIZcc6MhmVpK5R5mjl7NlDOeOHcyBE6eoqA++i9Cfr9xvUteZe4LFzDdn5vH+roqgBvyNz0VegMwdYOyQVK6YNow/fXmky4HRveUNTMxN6/Javn3+KPKzknh8xe4uu4vWeePEWaMH+e4zmxRPfHMaNU12nnyn/UDpxsM1NLQ6A/5av2/heKbmpfPgG9vbvdavbSolM9nKwslD2h2/cPJQ7ukQI7ujIlUoKFRFhRl63Y/PYcexerSmyw0gTjS2cuhkE1OGp5OWYOVIdRMVDa1MGJoWcB4twPayOsxmxZRhnkU9rQ4X24/VkZ5oZUJuGor2b+7u8gacbjfT8zpXO2xsc7DreAN5mUkUdJE1+LdzZkEmiR0GHmua2th/4hRjc1IZnJrQ5Tm2lNaSlmBh7JDus+T6Vgd7yhsYPTiFIWndF/HaeKSGnLQECgd1/SXQlSa7kx3H6hmTk0pON+0ORrD/NhG+4/UtlNQ0M7cwG3Mvvox7+75391kI9Dwjs5MZlhE4aBuMf1vRyKwuf5k0251sP1Yf8LOq0Ww8UktOWgKjuvksVDe1ceDEKUYNTmFogM9WcdUp6prtzBmZ1SmOlNQ0c7y+hYm5ab4vtiPVTVQ2tlI0MvB70uJwseNYPakJFiYNS8Pl1mw+WsvQ9MQuP7Pq9pWbtdZFXf4jvOKauRefOIXd6Wbc0NQuR/gHpSZgUoqqxjZqm+1UNLSSm57YZWAHT997Y6sTh9uNRnPw5CmUUozOSe30hgCkJ1potrs6VS3UaI5WN2M1mxie2f3Fl2zzzL1tbms/20WjKatrIdFiZlBq96s7Eywm2oJYol9R34rFpLr9ojBYzCrslaHGT8yMAPOVQ6WUQvL26DLytN79xoJkmxmLSVEf5liNcb1Zetg0JMVmIS3BQkVDK5rur9E2hxuzSXXb5ZRss5CdbKOiobXTZ7nN6catNcm27md8ZafYSEuwUFrTHLCKaWOrg/REa8A4kp+VRJLVzKGqJpze2FPTZCczydrll22S1UzhoGQaWh0cr2vl5Ck7Gk83Ua9prePyX/64KXrkgyv0n9ce0T25f9lWPeWx9/Scxz/Qlz79qW6xO7s9fltprR754Ar92qZS/cc1h/TIB1foZRtLujx+7cGTeuSDK/SHuyra3f/G5lI98sEV+m/dPNbQ1ObQhQ+t0M98uL/d/R/squjx+Q33L9uqz3pyVbfHHK46pQsfWqF/9f7eHs+ntdaLn/1c3/rS+qCO7eim36/VC3/1SViP7ejSpz/Vd/xpY0TOJQL79Qf79MgHV2i3293rc9376ld69s8+0E5X6Of6+bt79JiH3wmqHSu3H9cjH1yh391R3u1xt760Xi/6zWc9nm9HWZ0e+WDnz+H7O8v1yAdX6K+O1vR4ju2lnnP87O1d7e4vq23WIx9coV9ac6jLx24tqdWjHlqhH3htm955zHOepRuOdvt8brdbf++vm/SYh9/R5/zXR/ryZ7r/dwKbdBAxNm6Z+4lGzyKAb501osdjbygq4FSbk8ZWJ7+5cRaJ1u6/faflZZCbnshf1x3l5+/tCzjw6m9mQSY2S/t+95omO0+8s4fZIzK5bnbXjzUk2ywUDkppN6jqdLn573f3MHpwCt+c1XkFaUd5mUlUNrZ2O9f95S+PYDEpvnX2yB7PB54ZM3VhDKjWNNlZd6iGS6Z07isMh9VsktoyUeZ0u7GYVK/HR8Az4FndZGfTkeDWgPira3GQmWwNqh0XTx5KXmYSL33R/Wr0strOc9wDmZqXweVTc3nhs4PtNvPYW9GIUjB+aM/dgtPyM1gybwR//OJwu+381h/y9rePGtTVQ9sNlD6xYg9KwYUTu/8MKaX4r6unMyQtgWN1LQFLiYQjbsE90WrmyW9OC+oCmFuYxfVz8nnq+hlMCGLWhlKKhZOHsLW0DotZ8V/XdP88iVYzMwsyWe83Y+a/Vu6hocXBf10zvdPUx654Zsycng75t02lHKxq4sHLJwa1sCQ/KwmtobyLue71LQ6WbSrlGzOGMyTAiH8gmcm2sAZUP9hVgcutudxv84/esJhlnnu0OVw67NWpHX19whASLCZW7igP+bF1zfaAq1MDsZhN3Da/kA2Ha7rcF1VrHXB1alceuGwibU43v/lov+++vRUNjMxOJsVbuqAnD102keyUBP59+Q7f4Or6QzVkJFl7nDn2w4XjGDcklbWHqpk9IiuoLpaMZCvP3TSbS6cM5epZZ3hwHzU4haQe+r8MSimeun4GV84Y3vPBXldM8xz72OLJPQ7WAJw9Kpudx+ppbHWw9mA1r20u487zRwf1ZWKYmJvOkeommu1OTrU5efrDA8wtzOKSACPlgRgXb1cF05ZtLKXZ7uL2+aMC/j2QzCRrWGV/V+6sYER2MlOGp/d8cBCsJpME9yhzuNy9WsDkLyXBwoIJOby7syLkfU/rmgPXlenKDXMLSLGZu9yusrbZQYvDFXCOeyCjBqdw81kjeHVDqW/Rz94OZQd6kpFs5SffmMz2snr+vPYI4Fm8NLcwu8dkL8Fi5pfXz8BiUiyaFnxyNGdkFr+7pahTNc1wxS249zTY0ltfGzOILx66kOuLup4H7u+s0YNwa89Wco8s38GI7GTuvWhcSM85cVgaWnt2Zvr9Z4c4eaqNf180Keifyd0tZHK63Lz85RHOGpXN1LyMTn/vSlaylfoWR0gf0LpmO18Wn2TRtGER+YkPnoU1Un4guhwuN1ZL5D7Si6YN40RjW8iba9Q2B64r05X0RCvXFxXw9vbjARcsGp+HYIM7wL0XjSPJauYX7+2lxe7icHVTwMVL3Vk8fRgXjM/hl+/vY0tJLUeqmzl7dHbPD8TTPbPmwQv5p68F130aDXGf5x5NwfTRGWaPyMJqVjyyfCeHTjbx5NVTe+zb78hY1vzZ/ipe+OwQV0wbxqwRwRfbys1IxKQCL2T6YHclx+pauP3c4LN28KwSdGtobA2+VsiHuytxujWLpgWekx8Oi0lWqEab06UjmjRdOHEINouJlSGW561vtndZV6Yrt80vxOnW/CXAzkgdS/0GY1BqAt9bMIYPdlfyyvqjaE1ImTt4egye+OZUXFpzx583A933t3eUm5EY1paHkdKvg3sokmxmpudncvJUG9+cObzLei3dyc9KIsVm5rmPD+B0u3mgm0VPgVjNJoZlBK4O+eKaw4zITmbhpNAGOLPCKEGwckc5eZlJTAvhF0JPrGZZoRptdpe7V0XDOkpLtHL+uBze3Vke0i+/2mZHt6tTAxk5KIWFk4byyvqSTvuans7cg+tzN9w+fxS56Yn8/D1PSYFJIWbuAAXZydy3cDwnT7WRlmBhcoS6KWNBgrufCycOYXBqAo8unhzW400mxYTcNBwuzbfOHsnIMBYO5QWo6761tI7NR2u59ZzCLuvIdMXIoIKdMVPf4mBN8UkWTcuNWJcMeHZjksw9upwuHdb+qd1ZNC2X8vpWtga5gUerw0WLwxVSt4zhO+eOoqbJzvIO5UbKaltIS7SEvN4iyWbmXy8Zj8Plmd/e3SLE7tx+7ihm5Gdw/oSckD9/8STB3c/dC8aw5sGvB7U4qCszC7LISLLyLxeG1l9vyM9K6rSX6v98fID0REuXpX27E2rZ34/2VOJw6ZAGgoLhKfkrwT2aHBHO3MGz5N1qVrwb5KwZY+FbqN0y4KnzND0/g//95GC7a+VYCDNlOrp2dj6ThqUzPT8j6FlvHVnNJl7753N49sZZYT0+XiS4+1FKhdzP3tGPLp3Ah/edT3YYe42C56dnRUOrb5OS7WV1rNpzgjvOG01aYugfGGPWQrAbdqzcUcHwjERmFnQuxdAbnqmQ0i0TTZ6pkJH9SKcnWjlvXA4rd1QEtVlEMHVluqKU4gcXjaOkppnlX53O3j3TIIPvb/dnNimW3nE2v/3WnLAeb7BZTGdU1g4S3CMuyWYOeg56IPmZSbg1vkJCz6w6QGaylVvnF4Z1vlAy98ZWB58dqOKyqZGbJWOwSeYedQ6XG1uEu2UALp+ay7G6FraXdb0JjsG4zsLJ3MHTNTojP4PnVh/A4XL76riHMjmio4xka1jdRGc6Ce59jP90yK2ldXy8N/ysHUKrDPnx3hPYne6IzpIx9KbGjQiO0+2OyuyMiycPxWJSrNzZc9eML3MPM7grpfjhwvGU1rTwxuYy6podNNmDn+MuTpPg3scYfYtltS08s2o/WclW/umcwrDPZzYp0hMtQW/IOzQ9gdkhTN8MltVswhFEUTQRPocz8gOq4Pn1N3/sYN4Nomumzpe5h58pL5iQw8yCTJ77uJjD3n1Ww+1zH8gkuPcxxlz3FTvK+WRfFXeeP8a323u4slJsPc6WqW928Mm+Ki6fOizsgafuWM0mHAGq7InIcbgjP6BqWDQtl5Iaz6/J7gSzUUdPPNn7OI7VtfDMqgNAaAuYhIcE9z7GZjExND2Rz/ZXkZ1i49sRWOGWmWSltodumeVbymhzurstsNYbVhlQjbpozJYxLJo2jGSbmaUbut+kubbZjs1sIqmXExMuGJ/DrBGZfLa/CpDgHg4J7n2QcSHfdf7ooAsddScz2dZtbW6tNa9uKGV6fkZIpQ1CYTGZcLl1yHVKRPAivULVX1qilW9MH85b247T2MWmzuD5BRhsRcjuKKW4b+F4AFITQp/jLiS490njhqaRk5bALRGqS5GZ3H3mvqW0jn2VjSyZ13P55XDZvDVPpGsmeuwRri3T0ZKzRtDicHW5ETR4MvdwB1M7Om/cYOaNymbskNSIz94aCIK6EpRSlyml9imlipVSD3Vz3LVKKa2U6nELKNG1RxZNYuW95/l2d+qtrGRbtwOqr64vIdlm5hshVN0MlZFROqVrJmqcLo01inOxZ+RnMGlYOq9uKOnymLoQi4Z1RynFi/9UxEu3zo3I+QaaHoO7UsoMPA9cDkwGliilOq3PV0qlAT8A1ke6kQNNSoIlMttseWUkWWlodQZc/t/Q6uDt7ce5aubwXg/cdsfoC5a57tETzT538ATbm+YVsPNYAzu6mPMearnfnqQlWsNeEDjQBXMlzAOKtdaHtNZ2YClwVYDjHgd+DoS3ZbqIGmPmQkOAypD/2HqcVoc7ql0ygG+KngyqRk80Vqh2dNWsPBKtJl7dGDh7r2uxd7u/sYidYK6EPMB/iLzMe5+PUmo2UKC1fieCbRMR0tUqVa01/7e+hMnD0iNaATIQydyjL1orVP2lJ1pZPH04/9hyjKa29smC1tpby10GP/uCXn/NK6VMwK+B+4M49k6l1Cal1KaqqqrePrUIkq8yZIdB1e1l9ewpb2DJWSOiPmBlZJTS5x49Tld0Vqh2tGTeCJrsLt7uMLDa6nBjd7oH5FL/viiYK+EY4L+dUb73PkMaMBX4RCl1BDgbeCvQoKrW+gWtdZHWuignJ/R66SI8xoetvqV95v7qhhKSrGaumhm9gVSD0S1jl8w9ahwuHdU+d8PsEZlMGJrWaWC1rqV3dWVEZAVzJWwEximlRimlbMCNwFvGH7XW9VrrwVrrQq11IbAOuFJrvSkqLRYh823Y0XQ6cz/V5uStbcdZPH0Y6WHWrQmFEXScMhUyajwrVKM/ZVApxY3zCthWVs+u46cHVo3rqzerU0Xk9BjctdZO4B7gfWAPsExrvUsp9TOl1JXRbqDoPaP8qn8JglfXl9Bsd3Hz2bHZ49Eq3TJR5XJrtCYmmTvANbPySbaZefHz05taG5l7RhjlfkXkBTX3TWu9EljZ4b7Hujh2Qe+bJSIpLdGCSZ0u6tTmdPH7zw/xtdGDIl63vSsW6ZaJKmOg2hKDzB08ZXSXzBvBy18e4b6Lx1OQnewb08lKkcy9L5AVqgOAyaTISLL6PnxvbD7GicY2vv/1sTFrg00y96gygrsthhsyf/e8UZgUvPDZIaB3G3WIyJPgPkBkJduobbbjdLn53WcHmZ6fwfyxwe/k3lvGClWZChkdxvqBaNWWCWRYRhLXzs7nb5tKOdHY2uuNOkRkSXAfIDKSrdS3OHhnRzlHq5u5e8HYmNbrMGqeSHCPDmP1cTRrywRy1wVjcLrcvLTmCPUtDhKtpl5vVSkiQ4L7AJGVbKOmyc7/fnKQsUNSuWTy0Jg+v9VkBHfplokGYyzDeJ1jZdTgFC6fNoy/rjtKSXWzrE7tQyS4DxCZSVZ2lzewt6KRf75gTFQ25OiO1WIUDpPMPRqMsQzjdY6luxeM4VSbk/d3V0hp3j5EgvsAkZlsQ2vIy0yKyaKljizejFJmy0SHb7ZMjDN3gCnDM1gwIQetkcy9D5HgPkAYg1x3nj86ZnOh/clsmegyurvi8d4C3L3AM/NKBlP7jujVeBV9yvyxg9heNoQbigp6PjgKLGaZLRNNxusaixWqgcwblc3NZ42IyubqIjwS3AeIOSOz+cM/Zcft+X1VIWWbvagwyjrEK3MHePLqaXF7btGZdMuImPDVc3dK5h4Ndqd3nnucMnfR90hwFzEhhcOiy3hdY7lCVfRtciWImLDITkxRdbq2jHykhYdcCSImTi9iksw9Gk7PlpFuGeEhwV3EhMmkMJuUBPcoOT1bRj7SwkOuBBEzVrOSee5R4ozzPHfR98iVIGLGajLJCtUosftWqEq3jPCQ4C5ixmoxSeYeJcbraotxVUjRd8mVIGLGIn3uUeOQzF10IMFdxIzVbJKpkFHiiFM9d9F3yZUgYsZqlsw9WnxTIeNQFVL0TXIliJixmk2yQjVKnHEuHCb6HgnuImYsZpOvBoqILOMXkVn63IWXBHcRMzazksw9Shxujc1sium+uKJvk+AuYsZiNkmfe5Q4nG6pCCnakeAuYsYzFVK6ZaLB6dYyDVK0I8FdxIzNIpl7tNhdblnAJNqRq0HEjMUktWWixelyx2VzbNF3ydUgYsYqfe5R43BprBbplhGnSXAXMSPBPXocLrcsYBLtyNUgYsazQlW6ZaLB4XJLuV/RjlwNImYsZpNvJaWILKdLy1RI0Y4EdxEzVrMJu2TuUWGXzF10IFeDiBmrrFCNGqdLS10Z0Y4EdxEzVrMJh1OCezRIn7voKKirQSl1mVJqn1KqWCn1UIC//6tSardSartS6iOl1MjIN1Wc6SxmhcMt3TLR4HBrLBLchZ8erwallBl4HrgcmAwsUUpN7nDYFqBIaz0deB34RaQbKs58NpkKGTUOpxubdMsIP8F81c8DirXWh7TWdmApcJX/AVrr1VrrZu/NdUB+ZJsp+gOLyYTW4JLsPeKcblmhKtoL5mrIA0r9bpd57+vKd4B3e9Mo0T8ZKygle488zwpVCe7iNEskT6aU+hZQBFzQxd/vBO4EGDFiRCSfWpwBjBWUDpebRKs5zq3pXzwrVKVbRpwWzFf9MaDA73a+9752lFILgUeAK7XWbYFOpLV+QWtdpLUuysnJCae94gxmTNWTVaqRJ7NlREfBXA0bgXFKqVFKKRtwI/CW/wFKqVnA7/AE9hORb6boD4zZHLJKNfJkharoqMfgrrV2AvcA7wN7gGVa611KqZ8ppa70HvYUkAq8ppTaqpR6q4vTiQHM5g3udgnuEScrVEVHQfW5a61XAis73PeY3/8vjHC7RD9kZJZS0z3yZIWq6Ei+6kXMGJmllCCIPKdbMnfRnlwNImaMzNLulMw9krTWOFyyQlW0J1eDiBnJ3KPD6V0UJitUhT8J7iJmjMxSFjFFlvF6SuYu/MnVIGJG5rlHh/F6Sp+78CdXg4gZq2TuUWG8njJbRviT4C5ixtfnLpl7RDklcxcByNUgYsbirX0ii5giy9fnLrVlhB8J7iJmbBbJ3KPBCO42qQop/MjVIGLGyCylzz2yjAFVqecu/MnVIGJGBlSjQwZURSAS3EXMnA7u0i0TSaeDu3ycxWlyNYiYMTJLWaEaWcYKVQnuwp9cDSJmjBWUdqcE90hyOI0VqtItI06T4C5ixuarLSPdMpHkkMxdBCBXg4gZI7N0SOYeUcbrKQOqwp8EdxEzvqmQkrlHlDGGIZm78CdXg4gZpRRWs5KpkBFm95UfkMxdnCbBXcSU1WySDbIjzClTIUUAcjWImLKYlMxzjzCp5y4CkatBxJTNYpJumQhzSLeMCECCu4gpi0mCe6T5VqhKbRnhR64GEVNWi5KqkBHmq+cuVSGFH7kaRExZTSap5x5hdqnnLgKQ4C5iyjNbRjL3SJKdmEQgcjWImLLIPPeIc7jcmBSYJXMXfiS4i5iymE2yQjXCHG63TIMUncgVIWLKZlZSWybCHE7tK8omhEGuCBFTFpNJ6rlHmNPtlnK/ohMJ7iKmrBaTrxaKiAyHyy2DqaITuSJETFlNSmrLRJjDpbHKYKroQIK7iCmrWVaoRprD5ZYFTKITuSJETFnMskI10pwuLQuYRCcS3EVM2cyyQjXS7NLnLgII6opQSl2mlNqnlCpWSj0U4O8JSqm/ef++XilVGOmGiv5BMvfIc0pwFwH0eEUopczA88DlwGRgiVJqcofDvgPUaq3HAk8DP490Q0X/IH3ukedwaSn3KzoJ5ut+HlCstT6ktbYDS4GrOhxzFfAn7/+/DlyklJKrTXQiwT3yHC5ZoSo6swRxTB5Q6ne7DDirq2O01k6lVD0wCDgZiUaK/sNqVjS2Obn415/Guyn9RmltM3NGZsW7GaKPCSa4R4xS6k7gToARI0bE8qlFH7F4+nCO17eitfS7R8q4oal8Y/rweDdD9DHBBPdjQIHf7XzvfYGOKVNKWYAMoLrjibTWLwAvABQVFcmnewCaUZDJ8zfNjnczhOj3gumo2wiMU0qNUkrZgBuBtzoc8xbwT97/vw74WEtqJoQQcdNj5u7tQ78HeB8wAy9prXcppX4GbNJavwW8CPxFKVUM1OD5AhBCCBEnQfW5a61XAis73PeY3/+3AtdHtmlCCCHCJfOnhBCiH5LgLoQQ/ZAEdyGE6IckuAshRD8kwV0IIfohFa/p6EqpRmBfXJ48NIM5c8oonCltPVPaCWdOW6WdkddX2zpSa53T00ExLT/QwT6tdVEcnz8oSqlNZ0I74cxp65nSTjhz2irtjLwzqa2BSLeMEEL0QxLchRCiH4pncH8hjs8dijOlnXDmtPVMaSecOW2VdkbemdTWTuI2oCqEECJ6pFtGCCH6oT4R3JVS9yultFJqcLzbEohS6nGl1Hal1Fal1AdKqT67M4JS6iml1F5ve5crpTLj3aZAlFLXK6V2KaXcSqk+NyOhp03h+wql1EtKqRNKqZ3xbkt3lFIFSqnVSqnd3vf9B/FuUyBKqUSl1Aal1DZvO/8j3m0KV9yDu1KqALgEKIl3W7rxlNZ6utZ6JrACeKynB8TRh8BUrfV0YD/wcJzb05WdwDXAZ/FuSEdBbgrfV7wMXBbvRgTBCdyvtZ4MnA18v4++pm3AhVrrGcBM4DKl1NlxblNY4h7cgaeBB4A+2/mvtW7wu5lC327rB1prp/fmOjw7Z/U5Wus9Wuu+uogtmE3h+wSt9Wd49lDo07TW5Vrrr7z/3wjswbP3cp+iPU55b1q9//XZz3t34hrclVJXAce01tvi2Y5gKKWeVEqVAjfTtzN3f7cD78a7EWegQJvC97lAdKZSShUCs4D18W1JYEops1JqK3AC+FBr3Sfb2ZOor1BVSq0CcgP86RHg3/F0ycRdd+3UWv9Da/0I8IhS6mHgHuAnMW2gn57a6j3mETw/hV+JZdv8BdNOMbAopVKBN4AfdvhF3GdorV3ATO941XKl1FStdZ8e0wgk6sFda70w0P1KqWnAKGCbUgo83QdfKaXmaa0rot2ujrpqZwCv4NmVKm7Bvae2KqVuBRYDF8VzL9sQXtO+JphN4UWIlFJWPIH9Fa31m/FuT0+01nVKqdV4xu24H2kAAAJLSURBVDTOuOAet24ZrfUOrfUQrXWh1roQz0/f2fEI7D1RSo3zu3kVsDdebemJUuoyPGMYV2qtm+PdnjNUMJvCixAoTwb3IrBHa/3reLenK0qpHGOGmVIqCbiYPvx5705fGFA9E/y3UmqnUmo7nm6kPjmNy+t/gDTgQ+/Uzd/Gu0GBKKWuVkqVAV8D3lFKvR/vNhm8A9LGpvB7gGVa613xbVVgSqlXgbXABKVUmVLqO/FuUxfmA7cAF3qvy61KqUXxblQAw4DV3s/6Rjx97ivi3KawyApVIYTohyRzF0KIfkiCuxBC9EMS3IUQoh+S4C6EEP2QBHchhIiASBdxU0r9wlu8bI9S6lnvdNKgSXAXQojIeJkIFXFTSp2DZ/rodGAqMBe4IJRzSHAXQogICFTETSk1Rin1nlJqs1Lqc6XUxGBPByQCNiABTwGzylDaI8FdCCGi5wXgX7TWc4B/A/5fMA/SWq8FVgPl3v/e11rvCeWJo15bRgghBiJvkbRzgNf8ussTvH+7BvhZgIcd01pfqpQaC0zidMnuD5VS52mtPw/2+SW4CyFEdJiAOu8mP+14C6d1VzztamCdUVteKfUunlIdQQd36ZYRQogo8JY0PqyUuh48xdOUUjOCfHgJcIFSyuKtpnkBnjpHQZPgLoQQEdBFEbebge8opbYBuwh+R6/XgYPADmAbsE1r/XZI7ZHCYUII0f9I5i6EEP2QBHchhOiHJLgLIUQ/JMFdCCH6IQnuQgjRD0lwF0KIfkiCuxBC9EMS3IUQoh/6/wEqaX2UqGFFzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f(x) :\n",
    "    return (1.-np.cos(x))/x**2\n",
    "\n",
    "def g(x) :\n",
    "    return (np.sin(.5*x)**2*2)/x**2\n",
    "\n",
    "x = np.arange(-4e-8, 4e-8, 1e-9)\n",
    "df = pd.DataFrame(np.array([f(x), g(x)]).T, index=x, columns=['Numerical f(x)', 'True Value']);\n",
    "df.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## How did it happen?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cos(1.0999999999999999e-08) = 0.9999999999999999\n",
      "1-cos(1e-08) = 1.1102230246251565e-16\n",
      "true value of 1-cos(1e-08) = 6.0499999999999997e-17\n"
     ]
    }
   ],
   "source": [
    "x = 1.1e-8\n",
    "print(\"cos({:.16e}) = {:.16f}\".format(x, cos(x)))\n",
    "print(\"1-cos({:.0e}) = {:.16e}\".format(x, 1.-cos(x)))\n",
    "print(\"true value of 1-cos({:.0e}) = {:.16e}\".format(x, 2*sin(.5*x)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* Dramatic precision loss when subtracting numbers that are very similar \n",
    " * Most leading significant digits cancel\n",
    " * The result has much fewer number of significant digits\n",
    "* Catastrophic cancellation could easily happen in practice \n",
    " * Deltas are often computed by bump and re-value\n",
    " * Very small bump sizes should be avoided"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## A real example in Finance\n",
    "the CIR process is widely used in quant finance, eg, in the short rate and Heston model:\n",
    "\n",
    "$$ dr_t = \\kappa(\\mu - r) dt + \\sigma \\sqrt{r_t} dW_t\n",
    "$$\n",
    "\n",
    "the variance of $r_t$ at $t > 0$:\n",
    "\n",
    "$$\\text{Var}[r_t | r_0] = \\frac{r_0\\sigma^2}{\\kappa}(e^{-\\kappa t} - e^{-2\\kappa t}) + \\frac{\\mu\\sigma^2}{2\\kappa}(1-e^{-\\kappa t})^2$$\n",
    "\n",
    "how the catastrophic cancellation arises?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEbCAYAAAAh9sTfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8XOWV8PHfmabqIlmSqyy54U6zDaZjQmhLCwkJBMiSQHhTyIbNGxLYZMku6WVT3g0JNUASA6EkYMCEUExH2DJxwQ1k2bJky5ZsyXJRm/K8f9y5o5E0oxnJI0t35nw/H3/Q3LkaPVc8OnPm3HOfK8YYlFJKpRfXUA9AKaVU6mlwV0qpNKTBXSml0pAGd6WUSkMa3JVSKg1pcFdKqTSkwV0ppdKQBnd1VIjIKBFZKSKHRGTeUI9HqXSnwV0dLa3AvwBPDvVAlMoEGtzVUWGM8RtjGod6HEplCg3ualgQkR+LyC1DPY4jMRjHEC5lzU20TameNLirIScixcDngHuitt0sIpUi0iEiDw3Z4JIU6xhS5BfAnUlsU6obDe5qOLgeWG6MaYvatgv4AfCHIRlR/11P72NIhWXAEhEZl2CbUt1ocFdHjYgsB84D7hOR66OeuhB4PXpfY8xfjTFPA/tSPIYJIvKUiDSKyDYR+beo56aJSJOInBi1b6OInB1+vF1EbheRjSLSLCIPikh2vGPox5g+KSIV4a9LROQdEfkGgDGmHVgNnG/vH2ubUj1pcFdHjTHmImPMBGPMKcaYh6Kemg9sGeyfLyIu4FlgLTAR+Bhwi4icHx7fVuDbwJ9FJBd4EHjYGPNa1MtcgxVUpwHHAN9NwTEcB6wVkVnAq8DPjDG/jHp+U3gfEmxTKkKDuxoORgMHj8LPWQQUG2PuNMZ0GmOqgfuAq+wdjDH3AVXAe8B44Ds9XuO3xphaY0wT8EPg6vD2IzmG44AC4B3gpvAnlmgHw6+faJtSERrc1XDQDIwYyDeKyGsiYuL8e6vH7mXABBHZb/8D/gMY22O/+4B5wP8aYzp6PFcb9XUNMCHeMfRjbMcBM7AC9qwYhzkC2J/ENqUiNLiro6KPIGeAdVgljn4zxpxtjJE4/07vsXstsM0YMzrq3whjzEVR48wHfg08APyXiBT2eI3SqK8nY534JdYxJDM2ERkdfp2PAT8GviMinh4/czZWKSnRNqUiNLiro8IYI8BbgBv4H+CzdqADlgNnRe8vIp7wyUo34BaR7BhBr79WAgdF5NsikiMibhGZJyKLovb5DVBpjLkReB64u8drfFVEJoWD/neAv4S39zqGJB2L9YbThNUZ5AP+1X4y/DtYALzU1zaletLgro4KEckDQlgnKZ83xjwa9fQfgYtEJCdq23eBNuA24Nrw19/lCBhjgsDFwPHANmAvcD8wKjzGy4ALgC+Hv+UbwIkick3UyzwC/AOoBrZitWvGO4ZkHIeV9WOM6QR+Qvfs/RLgNWPMrqjvibVNqW5Eb5CtjgYROQkry20BzjPGNPR4/kdAgzHm10MxvmSIyHbgRmPMy3GeT/kxiMh7wA3GmA/62qZUTxrc1VEhIjcCQcAPHGuM+dYQD6nfEgV3pYYTLcuoo2U+VvnhUeB8ESka4vEoldY0c1dKqTSkmbtSSqUhDe5KKZWGjrRveMCKiopMeXn5UP14pZRypNWrV+81xhQn2m/Ignt5eTmVlZVD9eOVUsqRRKQmmf20LKOUUmlIg7tSSqUhDe5KKZWGNLgrpVQa0uCulFJpSIO7UkqlIQ3uKmVW1zRz14oqVtc0p2zfVO/X332T+d54r5fsvkcyHqXiGbI+d5VeVtc0c819FXQEQmR5XSy9cTELygri7nvVve8SDBl8nvj7rq5p5rP3VeAPhhLul+qfHe97P3tfBZ2BEG6X8IkTrDvs/e2fuwiGTGTb+FE51Le09doea197W8j0fzxK9UWDu0qJiup9dARCGKAzEKKiel/cIPX6hw34g9aCdf4+9rVfM9n9kvnZb33UmNTPjneMneGfEwgZnli9s9vz9jYRiF6Pr699o/V3PEr1RYO7SonFU8dEgpqIsHjqmLj7jsj2Rr72elxx943e7nX3vZ/bJQRCBo8r/n4As8aPBEAS/Ox4P8flEoIhQ3b4EwLANfdX4A+E8EZl3qtrmnttj7UvwJV3v0PI9H88SvVFg7tKifkTR5HlcdPmDzK9JL/P7HN/aycAXrfw5xtOjrvvlKK8yNc/u/LYuPstKCvgtOljeP3DvXxlybQ+f/bE0dZd8BaWF3DbhbP7lSUvKCvgnJklvL11L3+KGvfSGxdTUb2PxVPHRLYtKCuIuT3WtovmjefvG3b3+btQqr80uKuUWF3TTJs/yLTiPKobD9HuD5Ltdcfcd9V268ShP2goys+K+5qb6w9Evo7O9mM53BEEYHRO3/s1HbbeWOaMHzmgQDoyx0tBrq/b9y4oK4j5WrG2x9p2yvQxPLe+nvGj+3v7VaXi024ZlRKvbWnA6xa+ds4M/EHDurqWmPt1BIKsrd3PSVMKAdgUFcB72hj1XF1Ta58/v7bZer6p1d/nfs3hTw2dwYHdpMY+uZtK04rzAdjacCilr6symwZ3lRKvbWlkUXkhZx5jrUS6antTzP0+2HmAjkCIa06ejEv6Du6b6g9SlO8jy+NiRx/Bvd0fZM+BDgCaw5l5PHbm3hk+UdtfnYEQXrcM6HvjiQT3Rg3uKnU0uKsjtmt/G1v2HGTJzBIK83xML8mnMk5wt7efNr2IKUV5bKw/GPd1N9UfYPb4kUwqyKG2qS3ufnXNXc81tfYd3O3g3xkcYHAfhMy9KN/HyGyPBneVUhrc1RF7bUsjAGfPtLL2ReUFVNY0Ewr1Ln2s2t7M1KI8ivKzmD1+ZNzM3R8MUdVwiDnjRzK5MDdSdonFfs7ndiXO3MPB3z/AzN0fDOFzp/bPRkSYVpLP1obDKX1dldk0uKsjtmJLAxNH5zC9xCovLCwr5GB7gA8bumfloZChsqaJheXWCcXZ40eyc38bB9p718mrGw/TGQwxe/xISgtz+yzL1Iafmz1hZKTsEk/zYetnDTRz7wiE8KY4uINVmtHMXaWSBvdhaOW2ffzvqx854nL0jkCQd6r2cvbMYkSsWvSicutkqd0VY9vaeIj9rX4Whp+fE+453xyjNGNn9LPGj6C0IJeD7QFa4pwsrW1qJcvjYubY/MgJ03hSUXNPdVkGrODecLAj5hudUgOhwX2YWb29iavureB//vEh19xf0WeAHw5rklRub+ZwZ5AlM0si20oLcygZkdWr7m4Hezv4zxo/Aoh9UnXT7gP43C6mFedTWpgLEDd739HUSmlhLoV5WTQf9mNM/E6Yrm6ZgQf3rEEI7vannupGLc2o1HBccB8OAW0wPbJyB3ap2r6UPpbVNc1cfW8FP39xC9fc1/ebgL3/YPzeXtvSgM/t4tTpXVdWigiLygup7JG5V25voijfR/kYK1iPG5nN6Fxv7OBef5DpJfl43S5KC63+73h199qmNkoLcijM89IZDHG4Mxh3vEeaufuDg1WWsS7Ycno7ZLr/fTqJoy5iil64KdECUfb+Pa8GHM4CwRDvbWtCADv3jHc5ekX13kj22ZFgTZLVNc1cfV8FgQQLcA3Eii2NnDy1kFxf96m0qLyA59fXs3N/W+Sq0FU1TSwqL4yUb0SE2eNin1TdVH+AM2dYJ2jtzL02RuZujKG2qZWF5QUU5PoAqyMmP6v31DbGdGXuAy3LDEK3DFjH6HWLo+ruPf++Vm7bx3UPrEy40Js6OhwV3KMXbkq0QNTq7U1cfd97BEKJJ9pweRN4YnUddc1t3HreMbxVtZd3q5sYk+eLua/H1T3A9LUmybI1OyPBLNHiVP35XdQ1t1LVcIirFpX2es6uq1dub2Li8RPZ3dJObVMb1586pdt+s8eP5JGVNZGVEgH2Huqg8WAHs8Nlm5HZXkbnemOWZVra/BzsCDC5MJfC8O+q6XBn5A0h2qGOQNeiYUdQlkl1twxYa+eUjckbtsG957x4u2ovn39wFf5gCJcIpYVWu2owXBJL9PepBp+jgvviqWPI8rhoD4QIGZgT/uOP5devfBTJbPsKaMkuKzvYDncE+OVLH7KgrICvLJnOlYtKOf0nK3jw7W3892Xzuu1rjOG5dfWUjPAxb+JoXt3c0GfA2b6vq47b1+JU/f1d2C2QS2aV9Hpu1rgR5Gd5WLW9icuOn0hljVV/X1Te/fVmjR9Buz/E9n2HIxfz2CdYZ4dPuAKUFuRS29y7193uf59UkEuBHdzjnFS1O2Vg4DV3f9Ba9GswTCvOY+swrLlXbm/imvvfozNgBfJJBTnsaGqNfLoMWqvFcfFx41m+vj7yBpoOi6ANl8RvIBwV3BeUFbD0i4tZtmYXS9+r4ZGVtZw9syTyMd/28sY9vPnR3shjdx8rBb64YXdSy8qC/T96L4unFqX8f/R9b1bTeLCDu69dgIhQMiKby46fwOOVdfz7x49hdG5XBv/q5gbW72zhZ588lvPnjePkH73Mo6t2MH/S/F6vu7ulnber9lGY66OptZMffWJ+3LFXVO9N+ncBVr29tDCHqVELfNk8bhcnTB4dqbuv2tZErs8d6ZCx2Y831R+IBHe7TBMd3CcX5nZbjsBmZ/OTC3PJ9Vlr2cTrdbeDflG+b8BlmY5BytzB6ph5ZVPDoNX1k7VqexPPrt2F2yXsbG7jjY8aI/MiaAwul3Dlwkk8/c+dBEMGr8fF/1x5HAvKCvjcKeX8v1c+4vUPG8n2Dt9TerGC9spt+3h5UwMTRmeT6/Owats+/vq+tdZ+MmXg4SZhcBeRPwAXAw3GmHkxnhfgN8BFQCtwvTHm/VQP1GYvvFRamMMPnt/Eoytr+ezJkyPPVzUc5Ja/rGHexJHcet5MbvnLGjxuYda43ll+RyDIK5v2RB9Ln1ntZ+55l0DIkOWp4pEvpq7e33CwnXvfqOai+eO67XvDGVN4YnUdj6zcwVfOng5YWfuvX/6I0sIcPnHiRLxuFxfNH8+yNbv4zkWzyetRa3743e2EjOGB6xfyqbvf7fNj/6icrjcQTx9L7IL1RvDalkY+PmdsrzdX26LyQn718oe0tPpZtb2ZEyaPxtMjaE0vycftEjbVH+DiY62bV2yqP8DYkVmRMgvApMIcXtq4h1DICi42+yRraWFO5ER0vF53O+iPHZnN7pb2uMfWl8HqlgEruAdChh1NrZE3usFmz9O5E0bS0ubnydV13RKjsSOyOHXqGN6s2hu5wckvwoH8M4smx1wN8/9dfQKn/vgVHnhzG7/8zPFH5Tj6Ev23eOykUby4YTf//pc1BILWXJo/cRS7W9rYHV7CwhZ97suJZaZkMveHgN8Cf4zz/IXAjPC/k4Hfh/87qL5w2hRe29LI95/byMlTC5lWnE9Lq58bH64k2+vi3usWMmF0Dvd+biGfvuddfvb3zb3KGz94bhNbGw/z7Qtm8nhlHbtb2plUEHtlvscrawmEkqsnVlTv5boHViZ9t59fvfQR/mCIb50/q9v2WeNGcsaMIh5+Zzs3nj4Vn8fFii1dWbud3V190mT++v5Onl9Xz6ej6t+HOwIsrajh/LnjOGFyAadOG8Oza+v55nkzYwbk92uacQsEDdx05tQ+P8F87oFVBEKGlzftYXVNc8x9F5YXYAy89mEDm3cf4GvnzOi1T7bXzbTivG697hvDyw5EKy3IpTMYYs/BdsaP6vp/VNvUyuhcLyOyvRhj8LgkbnDfFxXc+7ooqi+DmVVPK+laQOxoBPfl6+r5t8f+GZnXALlRK3m6BD53ajlfXTI9ZrISbzXMUTlePr2olD+9W8O3LpjFuFHZg34s0D2Izxo3go8aDvHSxt3c/Xo1wZBBwscUvWZcMGSob2ljdI6PPQc6MOF9bjx9KufOKeFzD6yMlIFnjo1fBh6OEs5SY8wbQOyFQiyXAX80lgpgtIiMT9UA43G5hP/59HFkeV3c8tga2v1Bbn70fXbub+P31y5gQrhDY1F5If96SjkPv1vDym1dh/HMmp38qaKGm86cypfPns5Dn19E0Bh+8PymXj+rpdXPyxv3YIdDAzFLETYrWBtCpqu8Ec8z/9zJYyt38PE5YymP8Zo3nD6FPQc6eH79rl5Zu21hWQHTS/J5dNWObt/75Oo6DrQHuPGMqQBccuwEdjS1sn5n7xUbW1r9PL++ns+cNJniEVlU741f+62o3hepWYdCJu7xHV86Go9LuO/NakKmq7+9p+hlCDoDIbY2HuoV3Cfbve77ugflHU2tkedEhII8X9wLmaIz94GUZQJB6498MLplAKba7ZCDVHdfXdNsJTnPbuDyu97mK4+8HwnsAly7uIw/3nAS2V4XbrGO0/70tqCsgK8umZ505vqF06YQMoaH390+aMdy14oqVm7bx4d7DvLrlz/k0/e8y89f3MKnfv8Oc7/3Ipff9TZ3rdhKMHyMButE/1eXTMPnsY4x2+vid9cs4EdXzCcr6rjPnzeOk6aMYekXF/OF06fgc7t44K1tkddyglTU3CcCtVGP68Lb6lPw2n0aOzKbn1xxLF/682pO+ckrNB/28+Mr5vcKIt+6YCavbN7Dt55cywtfP5Od+1u5/a/rWVRewK3nzwSgbEweXz5rGr955SOuWlTKadOLAKsM8h9/W09Lm5+ffHI+2/Ye5oG3tvH3Dbu5cH7v97B1dfu7rYjYV3ljdU0z//74GgzwyqaGmBnwWccUM6Mkn/ve2MbIbC/r6lr46Sfnd8seRYSrFpXyg+c3sWX3QWaOG0EwZHjgrW2cOHl05DXPnzuO7zy9nmfX7uLYSaO7/Zyn1+ykIxDisydNJhAM8cIHuwkEQ73KKND9pGhfJ2hzfR7mThzF2tr9uF3C8ZNHx9xv1riRPLNmF/tbO6lvaccfNL3KaJF2yOa2bh8L65rbutXxC3N9cTP3ptZOPC5hTJ6PzmAIY0zcklIs9hvaYAX3kdleSkZkpbxjJhAMce8bW/nFPz6MlK7Kx+Ry3eIy6xNp+CTxJ06YGPcmI/1VWpjLBfPGsbSihpuXTO9VLuwPOyM/cfJocn0enltfzx/e3BbpzOnJAGfMKOLaxWX4gyG++cTayN2vvnXBLOumK7PGJnUjFfvTyexxI7j1yXX89tUqvn5u70+gfY275+/xraq9VG5v4owZxYNa5jmqJ1RF5CbgJoDJkycn2Ds5xSOycIvQfNiP2yUcE+OjU67Pw0+vOJbP3v8e//6XNaza3oTHJfzv1Sd2C5JfPnsaT6/ZyR3PfMALXz8Tn8fFU+/v5Pn19dx6/kw+s8gas8fl4rcrqvjCaVM4rrQrYHUEgtz6xDpKRmTzjY/P4NtPree8OWPj/g9cWlET+WMLBGOXekSEG8+YwrefWs/XHv0nxSN8XHHipF6vdcWJk/jZ37fw6Mod/Nelc3lp4x52NLVy24VdpZ5RuV7OnFHMc+vquf3C2ZHatTGGR1fuYP7EUcybOIolM0t4vLKO93d0rbserd1vBbnLj5/AdaeU9zlBF5UVsLZ2P0X5PrbsPhhzX7vlcVP9QepbrO6XnideJ4zORqT7VarBkKGuuZXz546LbCvI83brionWfLiTgjwfPo8LY6zv9/Rj+V5/wPqfNVgnVCE1a8zYQWXOBOsT0Z/frWFX1DkGl8CVC0v56pLpXH7CxKTLLf11w+lTWb5+N09U1nL9aVMSf0PU2BdPKaQwP4un3q/j969tjZsxC3DR/HF8bPZY/uNv6yNB/JZzj4kcw/hROUkdY1/H/akFk3hn6z5+88qHnDy1sFdC0zOQv1O1l+sfWoU/fDP1M2YUcagjwId7DtHS5keAu1/fOqgnaVMR3HcC0Y3Ok8LbejHG3AvcC7Bw4cKUfL6pqN6HsU97GBO3Fn7q9CI+Pmcsf9+wG7D+QHfub+tWD8z2uvmvS+fy+QdXcf9b1fzL/PF875kPOGlKIV86a1pkvy+dPY3HVu3gh8s38ZebFkeyv7terWLLnoM8eP0ilswqoWJbE8vX19NwsJ2SEd3rjm2dQd6saozUAfvKgO2yQ2tnEH8wxLq6ll7HWJjn47y5Y/nbP3dy24WzuP/NakoLc7oFPoBLjpvAK5sbeH9Hc6QXfU3tfjbvPsgPP2Gdkzh9RhEel/Dq5oaYwf3ZtbsYkeXhJ588Nu7dliLjyrdOijYc6OCa+ytiTubojpn6ljZ8Hle3W+wBZHncjB+Z3e2mHXsOWFm+fQWr/Xv4cE/s4Nh0uJPCXF8k8+6M88kkno6gdeXrYLVCAkwryWPZml39/lRhs9tZO6LKTqdNH8N1p5Txm1c+igS/6HLLYAWXBWUFnDh5NH94ezvXnVIeuY4heqzRJzqfXF3Hfz79QbdzANEEuPS4CVw4fxy3/GVN5Fi+cLp1fqhsTF7MTDkVxygifP/yeayp3c+X/7yaa04uY8nMYiaPyeP59bv4wXObCIQMLoGCXF/k/A5YN0Nftb2JORNGMbkwlw92tmAY/BuipyK4LwNuFpHHsE6kthhjBr0kY1s8dQw+j6vXpI1lzvgRvLTR6o4JhmL/YpfMLOH8uWP59Usf8dDb2zHArz5zfLeJmZ/l4ZZzj+G7T3/Ay5sa+PicsWzY1cLvXtvKFSdOjPR9f+2cGTyzZhf3vF7Nf148p9vPuWtFFY0HO/n+ZXM50B7o8yPw+zv2R87c2zXuWPtefdJknltXz09e2ExlTTN3XDyn1x/UuXPGkuVx8ezaXZHg/tjKWnK8bi49zupWGZHtZVF5Ia9taeiW+YP16eTvG3Zz3txxCQM7QHt4KYC+JnPxiCzG5PnYvPsA9S3tzBw7ImbQndRj6d/aqDZIW0GuL24rZHNrJwV53kjm3RkIkRv7GrGY7Dp91iBn7gfaA+w91EnxiPi3IIxlz4F27nx2QySwC/D508q545K5AJw0ZcxR79m+8YypfGXp+7y0cTcXzOsqY64K9877AyFErGSrvcd5kPPnjuWCueO5/W/rIn/fnzvV+qS4dET2oH3iiCc/y8NXl0zjm0+s47crqvjtiqpe+4SMVS4+b85Ynny/zmqqcLt4+Asnx7xx+mBeC5BMK+SjwNlAkYjUAd8DvADGmLuB5VhtkFVYrZCfH6zBxtKfGuGZx5RwzxvVCX+xV5w4iRc37KHhYAdet7C7pT1yCb3tqkWlPPj2Nn78wiZOn17ErU+soyDPxx1RQXxKUR6XHz+RP1fU8H/OnErJSCt739p4iHvesN4IrjulPOExLp46hixv4jewU6aOYXJhLg+9sx2fx8UxMdo/87M8nDOrhOfX7+aOS+bS2hlg2dpdXHLc+G73KV0yq5gfLd/cbfkAgNe3NHKwPcAlxyV3zvysmSXc+2bfv3MRCZ9Utcoy58S4KAqsIP7mR42Rx3aJprSgK7gXhk+o9myZBCtznzluRCTz7u+FTPbFOYNVc4euBcS2Nh5KKrivrmnm9S0N1Da38cIH9ZEyAMbqP/+XcHspDH7wi+X8ueMoLczh1y9/xLq6FtwuoWZfK//YuDvyZmkMzBo/ko/NKuG3K6oi5wBuOtO62fnkMblHPZDHs+dARyTREqwL+M6YXsRP/r45Mu7vXz6PBWUFfGphadI3Th8MCYO7MebqBM8b4KspG9EAJPs/OtlfbFXDoYSZssft4vYLZ3PjHys5+xcr2HOgg3uuW9DtYiOAr50znafX7OTu16u545I5GGO445kPyPa6uf3C2UkfXzLjdrmE02cU8ch7O/AHQtz48KqYZZBLjpvACx/s5r3qfWzbd5g2f5CrT+p+DuScWSX8aPlmXtvSwDUnl0W2P7uunsI8X+SEc6rGPmvcCB58ZzvBkGHWuJEx9yktyGXPgY7Izbdrm9sQIdIZBVbmHjJwoN3f6/9Fc6ufglxfJPPub8eMvf9gXmAUfcu9RFnd6u1NXHVfReRN58wZRfzg8vk0HuoYNldVul3CubPH8uDb29m822p3HZXj4eQphbyzdR+h8EVQ/3nxHBaUFXDq9KJhE8hj6Zlo2R1Ex5aOTnrcR+t4HHWFaiok84tNNlMuyPXiEuvd3CVQlN870yovyuMTJ0xk6Xs1fOmsqVRsa+Ltqn18/7K5/frYneyEsNei6asMsmRmCXk+N8+u28UHOw8wa9wIji/t3skyrTifSQU5rNjcFdxbOwO8vHEPV4Qvnkrl2GePHxk5adazDdJm19brmtuYXpJPbVMrE0bldMuko9eXiQ7uwZBhf2snhXlRNfcBBvfBzNzHjcwm1+dOeFemPQfa+eaT6yKB3SVw8tQxTB6Ty+QxucMmGILV925zCXzxjKncfM6MfvXODxfxkpXhOO6MC+7JSDbbrIjqmxeIWwv/2jnTrROdT63jvW1NTC3K47NR2XAqnT2zhPsSlEFyfG7OnTOWJ1fX4Q8abjitvNfJOxHhnFklPFFZF8mUX9nUQJs/yCXHTej1mkcqOqDHuwXe5Eg7ZGskuPe86MxeX6Znr/uBNj8hY2X29huTP9i/c/qD3QoJ1qevqcWxFxCzg2EgZPjDW9to6wzgcQkmXIIZrmu5nDGjmLtf3xqZk6dMsz71DceAmAynjFuDexzJZvjJnMwtG5PHWTOKeDW80Fbd/jbW1O4flAmS7BvTnPFWbznAn9/bwUXHToiZ4f/x3Rre29bEWccU8+zaXYwdmRX3YqQjcTDqDkQ3/bkyZjnJ7nW3O2Zqm1sjywLbCnPtzL17O6S9rkwqMndvP9onB2JacX6v9dB7dsEcMzafu689leZW/7ApwcRzNOvMqosG9yPQn0k7vWREJLgH4/S0p3JciV67I9B1Q4t4PfanTLNW4VyxuYHjS0fz2pZGrl1c1qsDJxUqo4JZ3K6a/Cx8Hhc7mlpp9wfZc6Cj19K+BXlWCaBnx4z9uDDP17VeSDD+TT1isTP3wVpbxjatOJ9la3fR1hkkJ7wYWvQCd4J13mRquD7vhGDplGw3nQzfZdscItnLss+fN46s8CXPw+Ej9GnTixOOJ9vr5tRpY1ixpYEXN+ymMxji0uNTX5KBruWc+xqPyyWUFljrhteFl/+N7nGHqJp7j7JMU1RwtzPvzkD/yjJ2ucjnTtwCeiSmFedjDFTvtUozm3cf4InPmYiEAAAX0ElEQVRK6yJwl0CW18Wp05I7oa0yl2buR8mCsgIe+eLw+Wia7HjOmVXCimc2cPfrWyktzOG4SaOGdDylhbnsaGqN2eMOkON1k+Vx9c7cw8G+IM8X+dTS31ZIe3+vZ5DLMiVda8y0dQb5wkOryPG5+dWnj2dXS9uwmD9q+NPgfhQNt4+myYzn7JklwAaqGw9zxQkTB3TVZCrHM7kwl9U1zV1L/RZ0D+4iQmFe7/Vl7Bp8Ya6PpkMDu9VepFtmkNdaLx9jBfdfvLiF+pY2Sgty+eMNJzGpoPfdpZSKR8syqk8NBzsiq2E+v75+yG98XFqQy8H2AB/sbCHL44rZTlqQ23tlyObWTrK9LnJ87kjm3d9b7R2NVkiADbsOIFgXaQVChv+8eLYGdtVvGtxVn6KX87VPvA4lu8b+dtU+SgtzY36SiJ25d0Y6aXwDvYgpeHQy9+jfsQvYGLXWvVLJ0uCu+mRf0DVcTgTb3TE797f1qrfbrDXdu7dC2itCAsP6IiYYfr9z5Uxac1d9Gm49ytGtj6Vx7ppVmOvtnbmHr06FqMx9gCdUBzu4D7ffuXImDe4qoeF0InhktpfRuV72t/p79bjbCvJ8tLT5u91spPlwZ+Tk60Azd/9RWFvGNpx+58qZtCyjHMcO0vGCu52h72/rKs00HY7K3Ae4KmRn0Fqe1jMIF3EplWoa3JXj5GdZFxEdag/EfL4gfOLU7nX3B0McaA9EtttlmXhr2MTTGQjhc7sGtR1UqVTR4K4cZXVNMyu3W+2Y3/nb+pitmdErQwLsD59cLQwvTeB2CSIDy9wHu1NGqVTRmaocpaJ6HyZ8Y2R/nNbMwh4rQ0ZfnQrWhU4+t2tA3TKDfTJVqVTRmaocxV6Js682wa7M3R/+b3hdmaj13X0eV/8zdw3uykG0W0Y5SjJtgqNzwytD2pn74e6ZOzCgzN0f1OCunEODu3KcRG2CWR43+VmeSMYevZa7zecZQFkmGDoqbZBKpYLOVJWWCvK6LmSyFwqzM3qwgvtA1pbRE6rKKXSmqrRUmOvrlrnnZ3nI8nStw+51D6DmHjRallGOoTNVpSVrfZmumnt0SQYGVnPvDAQ1c1eOoTNVpaXumbu/28lUsMoyHdoKqdKYzlSVlgryfJEumebDnRRG1dvBytz7W3P3a1lGOYjOVJWWCvN8HO4M0u4P0hS13K9tQN0ygVDk/qtKDXca3FVasteR2d/qp7m1s9sFTDDAi5iCIXyewb05tlKposFdpSV7HZldLW20dgZ7Ze5et+APmH69prZCKidJaqaKyAUiskVEqkTkthjPTxaRFSLyTxFZJyIXpX6oSiXPztyrGw8D9O6W8bgHmLlrWUY5Q8LgLiJu4C7gQmAOcLWIzOmx23eBx40xJwBXAb9L9UCV6g87mG9tPAR0BXvbgBcO08xdOUQyM/UkoMoYU22M6QQeAy7rsY8BRoa/HgXsSt0Qleo/uwyztcEK7r0zd9GFw1RaS2amTgRqox7XhbdF+y/gWhGpA5YDX4v1QiJyk4hUikhlY2PjAIarVHJG51g1dztzt2vwtoEuHKZryyinSNVMvRp4yBgzCbgI+JOI9HptY8y9xpiFxpiFxcXFKfrRSvXmcbsYleOlZl8rEKMs089WyFDIEAhpn7tyjmRm6k6gNOrxpPC2aDcAjwMYY94FsoGiVAxQqYEqzPMRCBlEYFRO98zd28+LmOwSjgZ35RTJzNRVwAwRmSIiPqwTpst67LMD+BiAiMzGCu5ad1FDqiB8VeqoHC+eHuUUn8dFIGQIhZJrh4wEdy3LKIdIOFONMQHgZuBFYBNWV8wGEblTRC4N7/Z/gS+KyFrgUeB6Y98LTakhYp9E7XkBE3Rl4MmeVLVLOJq5K6dI6mYdxpjlWCdKo7fdEfX1RuC01A5NqSNj19l7XsAEXRl4ZzBEtjfxVaeR4K6Zu3IInakqbdmZe8+TqRCVuSd5UtWuz2u3jHIKnakqbdkZe882SIjK3JMM7lqWUU6jM1WlrUjmHqMsY2fgyXbMdGhwVw6jM1WlLftE6pb6g6yuae723EDLMhrclVPoTFVpq+FgBwCvf9jINfdXdAvwdpBO9m5MekJVOY3OVJW2dh9oA6yFj/yBEBXV+yLP+fpZltGLmJTT6ExVaeusY0rI9rpwC3g9LhZPHRN5bsBlGc3clUMk1eeulBMtKCtg6Y2Lqajex+KpY1hQVhB5bqAXMWkrpHIKDe4qrS0oK+gW1G39Lctot4xyGp2pKiN5+9nn7g9aq2lkaXBXDqEzVWWkgXbLaFlGOYXOVJWR7AzczsgT6QwEAS3LKOfQmaoy0kDLMhrclVPoTFUZqasVMpjU/p2RhcNk0MakVCppcFcZydfPskyHXqGqHEZnqspIdgbenz53n9uFiGbuyhk0uKuMZGfgyXbL+IMhLckoR9HgrjKSiODrx02yOwMhPZmqHEVnq8pYXrf062YdGtyVk+hsVRnL53H1a+EwvYBJOYnOVpWx+hPcO4KauStn0dmqMpa3vzV3zdyVg+hsVRnL53HRkWRw92vmrhxGZ6vKWD538mUZzdyV0+hsVRnL59FWSJW+dLaqjNWfzF3LMsppkpqtInKBiGwRkSoRuS3OPp8WkY0iskFEHkntMJVKvX51ywS0FVI5S8Lb7ImIG7gL+DhQB6wSkWXGmI1R+8wAbgdOM8Y0i0jJYA1YqVTxul0c7ggktW+nZu7KYZKZrScBVcaYamNMJ/AYcFmPfb4I3GWMaQYwxjSkdphKpZ7P4+rX2jJZmrkrB0lmtk4EaqMe14W3RTsGOEZE3haRChG5IFUDVGqw+Dyufq0KqWUZ5SQJyzL9eJ0ZwNnAJOANEZlvjNkfvZOI3ATcBDB58uQU/WilBkYXDlPpLJnZuhMojXo8KbwtWh2wzBjjN8ZsAz7ECvbdGGPuNcYsNMYsLC4uHuiYlUqJ/nXLGA3uylGSma2rgBkiMkVEfMBVwLIe+zyNlbUjIkVYZZrqFI5TqZTrT7eMlmWU0yScrcaYAHAz8CKwCXjcGLNBRO4UkUvDu70I7BORjcAK4FZjzL7BGrRSqWCtLZP4NnvGGO2WUY6TVM3dGLMcWN5j2x1RXxvgG+F/SjlCspm7/QaQpcFdOYjOVpWx7G4ZKzeJz+6o0dvsKSfR4K4yli8crBOVZuzsXhcOU06is1VlLLuGnqjXPRLcPe5BH5NSqaLBXWUsOxNPVHf3a1lGOZAGd5Wx7Ew80YVMHZHMXf9clHPobFUZy87EE2Xu9vPaLaOcRGerylh2Jp5o8bCusoz+uSjn0NmqMpadiScqy9gnXLUso5xEZ6vKWN4kT6hqK6RyIp2tKmMl3Qppl2U0c1cOorNVZSw7E/dr5q7SkM5WlbHsTLwjyYuYtFtGOYnOVpWx+nsRk55QVU6is1VlrKS7ZQLaCqmcR2erylhJd8to5q4cSGeryliRbplkT6hqcFcOorNVZaz+tkJqt4xyEp2tKmPpRUwqnelsVRkrK8nM3R8M4XEJLpcu+aucQ4O7ylj9ydy1U0Y5jc5YlbHcLsHtkqRaIfVkqnIanbEqo/ncriRaIY0Gd+U4OmNVRvN5kgjugZCeTFWOozNWZTSv20Vn0PS5T2dQyzLKeXTGqoyWlVTmHtTMXTmOzliV0XweVxKtkAavR9sglbNocFcZzeuWpNZz18xdOU1SM1ZELhCRLSJSJSK39bHfJ0XEiMjC1A1RqcGTTOaurZDKiRLOWBFxA3cBFwJzgKtFZE6M/UYAXwfeS/UglRosybVC6kVMynmSmbEnAVXGmGpjTCfwGHBZjP2+D/wUaE/h+JQaVFa3TOLMXe/CpJwmmRk7EaiNelwX3hYhIicCpcaY51M4NqUGXVJ97toKqRzoiGesiLiAXwL/N4l9bxKRShGpbGxsPNIfrdQRS6YV0q9lGeVAyczYnUBp1ONJ4W22EcA84DUR2Q4sBpbFOqlqjLnXGLPQGLOwuLh44KNWKkW8bldya8tocFcOk8yMXQXMEJEpIuIDrgKW2U8aY1qMMUXGmHJjTDlQAVxqjKkclBErlULaLaPSVcIZa4wJADcDLwKbgMeNMRtE5E4RuXSwB6jUYEq2W0aDu3IaTzI7GWOWA8t7bLsjzr5nH/mwlDo6fB4ty6j0pDNWZTSv20VHH5m7MUYzd+VIOmNVRkvULRMMGYzR+6cq59EZqzJaorKMfbLVq5m7chidsSqjed0uQgYCcQK8ndVr5q6cRmesymh2LT1eO6S9XWvuyml0xqqMZmfk/kDsuzFp5q6cSmesymh2Lb0jGIz5fCS4a+auHEZnrMpoWeGMPF7HjJZllFPpjFUZzQ7a/jg3ybbLNbpwmHIanbEqo3kTZu5WuUYzd+U0OmNVRot0y8QJ7h16QlU5lM5YldEStULa5RqfR47amJRKBQ3uKqN53VbQjluWiWTu7qM2JqVSQYO7ymhZiS5i0lZI5VA6Y1VGszPyeJm7ve6MneEr5RQa3FVG84Zr6fEWD9PMXTmVzliV0XwJWiE79CIm5VA6Y1VGS9QK6ddWSOVQOmNVRotk7roqpEozOmNVRkuUueuqkMqpdMaqjJb4IqYQLgGPBnflMDpjVUbzRtZzj5+566Jhyol01qqM5nEJIvEz945ASOvtypF01qqMJiL43K4+L2LK0uCuHEhnrcp4Prerz+UHtCyjnEhnrcp4Pk/8zL0zqGUZ5Uw6a1XG6yu4+4MhbYNUjpTUrBWRC0Rki4hUichtMZ7/hohsFJF1IvKKiJSlfqhKDQ6v29Xn2jJallFOlHDWiogbuAu4EJgDXC0ic3rs9k9goTHmWOBJ4GepHqhSg8XniV9z124Z5VTJzNqTgCpjTLUxphN4DLgsegdjzApjTGv4YQUwKbXDVGrw9NUt06nBXTlUMrN2IlAb9bguvC2eG4AXYj0hIjeJSKWIVDY2NiY/SqUGkZW5m5jPac1dOVVKZ62IXAssBH4e63ljzL3GmIXGmIXFxcWp/NFKDZiVuQdjPqfdMsqpPEnssxMojXo8KbytGxE5F/gOcJYxpiM1w1Nq8Pk8Llo7AzGf6wxo5q6cKZlZuwqYISJTRMQHXAUsi95BRE4A7gEuNcY0pH6YSg0en8eFP25ZxuDVzF05UMJZa4wJADcDLwKbgMeNMRtE5E4RuTS828+BfOAJEVkjIsvivJxSw47XLX2fUNXMXTlQMmUZjDHLgeU9tt0R9fW5KR6XUkeNz+PWVkiVdnTWqoyXaOEwn1uO8oiUOnIa3FXG83mkz4XDNHNXTqSzVmW8Pi9i0lZI5VA6a1XGs7plegf3YMgQDBl8bvcQjEqpI6PBXWU8b5zM3Q74Xo/W3JXzaHBXGc/ncREIGUKh7r3uHeGAr62Qyol01qqMZ9fUe55UtTN3vc2eciKdtSrj2Zl5z+Bul2p0PXflRDprVcaLZO6B2MFdu2WUE+msVRkvkrkHYpdlNLgrJ9JZqzKeXXbp2Q7ZoWUZ5WA6a1XGi1uW0cxdOZjOWpXx7ODd0bMsE36cpZm7ciCdtSrj+eKUZTojFzHpn4lyHp21KuMl7JbRzF05kM5alfHiXcSkrZDKyXTWqowXr1smUpbRzF05kM5alfHi9bnbj3X5AeVEOmtVxovXLaOtkMrJdNaqjNfVLdN9VUi/XsSkHExnrcp4ehGTSkc6a1XG6wruwW7btRVSOZnOWpXxvG7rTks9yzKd4cf280o5iQZ3lfH66nP3uV2IaHBXzqPBXWU8u+zSq1smENJ6u3Isnbkq44kIPrer10VM/qAGd+VcSc1cEblARLaISJWI3Bbj+SwR+Uv4+fdEpDzVA1VqMHndEvMiJq23K6dKGNxFxA3cBVwIzAGuFpE5PXa7AWg2xkwHfgX8NNUDVWow+TyumK2Qmrkrp0pm5p4EVBljqo0xncBjwGU99rkMeDj89ZPAx0TPQimHWVu7n9U1zZHHew60c7gj0G2bUk7hSWKfiUBt1OM64OR4+xhjAiLSAowB9qZikEoNptU1zexv9dPc2sKnfv8OEwtyAKhrbgPgmvsrWHrjYhaUFQzlMJXql6P6mVNEbhKRShGpbGxsPJo/Wqm4Kqr3Rb42QK7PTa7PHdnmD4S67aOUEyQT3HcCpVGPJ4W3xdxHRDzAKKDXX4Mx5l5jzEJjzMLi4uKBjVipFFs8dQxZXhdugWyvix9fcSw/vuJYssPbvB4Xi6eOGephKtUvyZRlVgEzRGQKVhC/Cvhsj32WAf8KvAt8CnjVGGNQygEWlBWw9MbFVFTvY/HUMZHyS6xtSjlFwuAerqHfDLwIuIE/GGM2iMidQKUxZhnwAPAnEakCmrDeAJRyjAVlBb0CeKxtSjlFMpk7xpjlwPIe2+6I+roduDK1Q1NKKTVQ2sSrlFJpSIO7UkqlIQ3uSimVhjS4K6VUGtLgrpRSaUiGqh1dRBqBmgF+exHptbRBOh1POh0L6PEMZ+l0LJD88ZQZYxJeBTpkwf1IiEilMWbhUI8jVdLpeNLpWECPZzhLp2OB1B+PlmWUUioNaXBXSqk05NTgfu9QDyDF0ul40ulYQI9nOEunY4EUH48ja+5KKaX65tTMXSmlVB8cG9xF5Psisk5E1ojIP0RkwlCP6UiIyM9FZHP4mP4mIqOHekwDJSJXisgGEQmJiGO7GRLdGN5JROQPItIgIh8M9ViOlIiUisgKEdkYnmdfH+oxDZSIZIvIShFZGz6W/07Zazu1LCMiI40xB8Jf/xswxxjzpSEe1oCJyHlY6+AHROSnAMaYbw/xsAZERGYDIeAe4JvGmMohHlK/hW8M/yHwcaxbS64CrjbGbBzSgQ2QiJwJHAL+aIyZN9TjORIiMh4Yb4x5X0RGAKuBy534/yZ8r+k8Y8whEfECbwFfN8ZUHOlrOzZztwN7WB7WHdIcyxjzD2NMIPywAuuOV45kjNlkjNky1OM4QsncGN4xjDFvYN1rwfGMMfXGmPfDXx8ENmHdx9lxjOVQ+KE3/C8lscyxwR1ARH4oIrXANcAdifZ3kC8ALwz1IDJcrBvDOzKApDMRKQdOAN4b2pEMnIi4RWQN0AC8ZIxJybEM6+AuIi+LyAcx/l0GYIz5jjGmFFgK3Dy0o00s0fGE9/kOEMA6pmErmWNRajCJSD7wFHBLj0/yjmKMCRpjjsf6tH6SiKSkbJbUnZiGijHm3CR3XYp1p6jvDeJwjlii4xGR64GLgY8N93vQ9uP/jVMlc2N4NUTC9emngKXGmL8O9XhSwRizX0RWABcAR3zie1hn7n0RkRlRDy8DNg/VWFJBRC4AvgVcaoxpHerxqK4bw4uID+u+wMuGeEyKyEnIB4BNxphfDvV4joSIFNudcSKSg3UCPyWxzMndMk8BM7G6MmqALxljHJtZhW8ungXsC2+qcGr3j4h8AvhfoBjYD6wxxpw/tKPqPxG5CPg1XTeG/+EQD2nARORR4GyslQf3AN8zxjwwpIMaIBE5HXgTWI/19w/wH+F7PTuKiBwLPIw1x1zA48aYO1Py2k4N7koppeJzbFlGKaVUfBrclVIqDWlwV0qpNKTBXSml0pAGd6WUSlKqF2ATkb+LyH4Rea7H9odEZFt4YcQ1InJ8f19bg7tSSiXvIayLjFLl58B1cZ671RhzfPjfmv6+sAZ3pZRKUqwF2ERkWjgDXy0ib4rIrH683ivAwVSPEzS4K9WNiDwoIpeIyGgRWR6+IEupvtwLfM0YswD4JvC7FL3uD8P3d/iViGT195uH9doySg2B+ViXtj8D3G6MeWeIx6OGsfDiZacCT1irIgDWleaIyBVArKtNdyZxxfbtwG7Ah/Xm8e04rxWXBnelwkTEBcwGHgHu0sCukuAC9odXdewmvKDZgBY1M8bUh7/sEJEHsT4R9HtgSinLDGAXcD3wpfDKg0rFFV5qeJuIXAnWomYictyRvm74blP2ImmXM4BVIjW4K9VlPtbNEl7F+mP63BCPRw0z4QXY3gVmikidiNyAdbOgG0RkLbCBftyxS0TeBJ4APhZ+Pbtcs1RE1mMtjlYE/KDfY9WFw5SyhG9OvMcY8zsROQX4EzAr6vaHSjmGBnellEpDWpZRSqk0pMFdKaXSkAZ3pZRKQxrclVIqDWlwV0qpNKTBXSml0pAGd6WUSkMa3JVSKg39f2mNLy+usY+XAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "k = arange(-3e-8, 3e-8, 1e-9)*1e-7\n",
    "t = 1\n",
    "v = (1-exp(-k*t))/k\n",
    "plot(k, v, '.-')\n",
    "title(r'$\\frac{1}{\\kappa}(1-\\exp(-\\kappa))$')\n",
    "xlabel(r'$\\kappa$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAA0AAAASCAYAAACAa1QyAAAABHNCSVQICAgIfAhkiAAAAHZJREFUKJFjYKACCGFgYJjMwMBwmIGB4RMDA8N/BgaGJYQ0XYAq/MzAwHCdWE2ODAwMqgwMDIwMDAwOuDSxoPH3EzKVgYGBgYkYRaOaBlwTeuQGQDEDAwODBJS2ZGBgWABlv2FgYChBN6SBAZJ0cOEH5LiMzgAA6XoX52TB9a4AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$1$$"
      ],
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing(use_latex = True)\n",
    "\n",
    "k = sp.symbols(\"k\", positive = True)\n",
    "e = (1-sp.exp(-k))/k\n",
    "sp.limit(e, k, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd4XNWdv99zp6hXy92y5F4oBgsXEhx6QjZ2SLKBYEjCLi3JhueX7KYusCkkhOymLRtIwBDTYhswphpTYqrBliXLvRdZo+KiNlazpGnn98e5d4o0kiVbZTQ67/P4mbln7tw5tqX7mW8XUko0Go1GozlbjMHegEaj0WiGNlpINBqNRnNOaCHRaDQazTmhhUSj0Wg054QWEo1Go9GcE1pINBqNRnNOaCHRaDQazTmhhUSj0Wg054QWEo1Go9GcE1pINBqNRnNO2Ad7A/2JEGIJsCQtLe3O6dOnD/Z2NMOFAwfU44wZg7sPjeYcKSkpqZVSjjzTeWI49Nq65JJL5JYtWwZ7G5rhwhVXqMcPPhjMXWg054wQokRKecmZztOuLY1Go9GcE1pINBqNRnNOxLWQCCGWCCGWNTQ0DPZWNBqNJm6JayGRUr4upbwrIyNjsLei0Wg0cUtcC4lGo9Fo+p8hKSRCiFlCiEeFEC8KIb4z2PvRaDSaWKTE5eaR9w9T4nL36+cMuJAIIZYLIaqFELs7rF8nhDgghDgshPhpd9eQUu6TUn4buBH4dH/uV6PRaIYiJS43Sx8v5A/vHOCWJwr7VUwGwyJ5CrgufEEIYQMeAT4PzAaWCiFmCyEuEEKs7fBnlPmeLwJvAOsGdvsajUYT2xQdrePfn9+OxxcgIMHrC1BYWtdvnzfgle1Syo+EEPkdlucDh6WUpQBCiOeA66WUDwKLu7jOa8BrQog3gJX9t2ONRqMZOvzxnQP833uHg8eGAIfdYOHkEf32mbHSImU8UBF2XAks6OpkIcQVwFeABLqwSIQQdwF3AUycOLGv9qnRaDQxSdHROh54Yx87KkPlDgL49NQcvn/NdArysvrts2NFSHqFlPID4IMznLNMCHEcWOJ0OgsGYl8ajUYz0JS43Dz5yVHe2Hmcjg2vbIbodxGB2BGSKiA37HiCuXZOSClfB16/5JJL7jzXa2k0Gk0sUeJys6akgue3VOIPdO6ZaDcE919/fr+LCMSOkBQD04QQk1ACchNw87le1Or+O3Xq1HO9lEaj0cQMKzeX81+v7u4kIAJw2AQ3XJLLV+ZOGBARgUEQEiHEKuAKIEcIUQn8XEr5NyHE3cDbgA1YLqXcM9B702g0mljGskJWFVcQ3rh9sATEYjCytpZ2sb4Oncqr0Wg0USlxubn58ULafYGIdZuAm+ZPHBQBsRiSle09Rffa0mg08UBJWT0/XbOzk4jYDcGvvnQBD3z5gkETEYidGEm/0Bcxks2ldRSV1fOpKTmD+h+l0WiGJ49+cIT/fmt/MCPLEEpABsuNFY24FpJzzdpaXVzBj9bsRACPOA7zs8Xn4T7tYeHkETHxn6fRaOKXkrJ6fv/OATaV1gfXBqoupLfEtZCcK4eqmwGQgMcb4L5XdiElJDgMVtyxMKb+IzUaTXxQ4nKzotDFy9ur6DgJfaDqQnpLXMdIznWw1WXTctR1AMMQBKQSlf7uW6PRaIYnKzeXc8OjG3lpW2cRGci6kN4S10JyrsH2y6bmIICFk7O5//rzg+tCiGDfmoFq06zRaOKXEpebH7ywnXtf3kWgQ1qv0ya4ZcFEnv/Wpdy8IDbbPcW1a+tcg+2GIchIdjBtdBpL5+fyX6/uwh+AlAQbG4/UcuBEE794bQ9ef0C7uzQazVlR4nJz07JNeP2RJkgspPX2lLgWkr5okZKZ5ODUaS8tHj/+AEwfncrBk8384Z2DCAhmUnhMd1es/4drNJrY4aMDNfzwxR2dRMRyY8WqBdKRuBaSvmCh/RBzT+yj5fCXAcgfkcLBk6EgvEVAQpW7lRKXW4uJRqPplhKXm4ffO8j7B2qDa7GY1ttTtJB0x7aVPNjwIwIIxEsrmSt+ykUTZ/De/mp8piNTABOykqhwt7KqqJyXtlVqF5dGo+mS5R+X8qs39kUE0w1iM623p8R1sP1cs7ZwHwXAhoSAl4XGPi7JU4F3uyEwhEoF/sy0kYDO6NJoNF1T4nJz+1NF3L92X6c+WU6HMWRFBOLcIjnnGMm0a/Fv+CN26QMpqZepZCU7uHnBRGaMSaOwtC6YvbWyuBwp+38SmUajGVqUuNw8X1zOiyWVdOz2PpQC6t0R10JyzuTO571JP+KaIw8iRICfO56l7dRXYPRlFORlRfzHXzppBPtONPLErfOG9A+ERqPpO7pq9w5DL6DeHUNWSIQQKcCHwC+klGv763OyRTMSgYHEgQ/HiUKYcVmn8y6YkMGWcjcX52b211Y0Gs0QocTl5oXicl7YUhmRlDPY7d77i8GYR7IcWAxUSynPD1u/DngINY/kCSnlb89wqZ8AL/TbRk2axy7Ed/iv2PDhx0bC5M9EPW9CVhIeX4Da5nZKa5v5+FAdV84cFTc/KBqNpmeUuNwsfbwQTwy2e+8vBsMieQp4GHjGWhBC2IBHgGuBSqBYCPEaSlQe7PD+24A5wF4gsb836x8/j3u9t/F75zL+5riJf8udDxVFULYB8hdB7nwAZvsP8G+2V9m+0c93PrLjD0ge31DKyjt1BpdGM1woLK3lh6t3dhKReHJjRWMwBlt9JITI77A8HzgspSwFEEI8B1wvpXwQZb1EIIS4AkgBZgOtQoh1UspAx/P6gsxkBx8G5qiDhDQlIk99AfxeMGzwT3+A0bO5+P1vMMfuIbD5FebI/2Qr02n3BfjTPw7y79cO3WwMjUbTM3739n4eef9I8Hgo14X0lliJkYwHKsKOK4EFXZ0spbwXQAjxL0BtNBERQtwF3AUwceLZfwvISHJSRwZ+KRhra4B9a8HvUS8GfLDuB3Dx1xH+duwC/FKlCW/1Twfg48O1bD5ax43D4IdJoxmOFB2t41dr97KrqjG4Fqvt3vuLWBGSs0JK+VQ3ry0TQhwHljidzoKz/YzMZAcBDGrIZIxwQ8aFkScEAtDeHGyX4sVOYWAWs8ems/e4+sHy+iUrNpezZmv0YsUSlzuYSjwcfug0mnigxOVm+celrNt1go45WbHa7r2/iBUhqQJyw44nmGvnRF/02spIcgBwUmYxQtaDzaleEDaQfrDZIXce7H4RAXzf8x22yuk88dnp3L1yK+3eQPCHzNuhH1eJy83KzS5e3X6MgJQ47brxo0YT65S43LxYUsELWyq7TesdTr/HsVLZXgxME0JMEkI4gZuA1871oudc2Q44bAapCXaqZRaZ/jqoLwVbAixdpU645HZwJAfPzxGNjM9M4ppZo1lxx0KWLlD6KIgsVixxuVm6rJA1W6vwBSQBCW3eAH9854BuSa/RxCgrN5dz42ObWFVUESEiQ6Xde38xGOm/q4ArgBwhRCXwcynl34QQdwNvozK1lksp9wz03roiI8lBdXMmaZ5SqD8K2ZNg2mfBmQZIcLtAGDQ4RvIZ/07cud8ECBYtFpXW4/VL/vi1i4LfUgpL6/D4O+cHfHKkjpInCrVlotHEECUuN2tKKlhVXNGpvUk81oX0lsHI2lraxfo6YN0Ab6dHZCY7ONmURZLXDTUHIGcaCAFZ+UpYEjMgfQInRD6Xt2+knL3A3OD7L5iQyeYOLeYXTsoOPnfaBLPGprOjUllOHV1gGo1m8Chxubn58ULah1FdSG+JFddWv3CuExItMpMdnMT8Qak7BFmT1PPsfNXY8ZSLFlsak9wbSRA+vnHwe+wvXh98/6ScFI41tNHq8QfXslJUrOXa2aNZddel/GzJedgMAeh+XRpNrLClrJ4fv7ijk4jYDcGvvnQBD3z5gmEvIhDnQtIXMRKAzCQnJ2VY65NsU0iyJim3lruMFo8fAyUUDny4974XPH3yyBQAyupagmsfH1ZzCO77wqygC+zOReq6f7hhjv7h1GgGmb+8f4gbHt3EkRr1e2uI4R0H6Y64FpK+skg8vgDVMuzGnj3ZfJwE/nZoPokYNwcvjmB3z6xZVwVPn5SjhKS0JiQkB4rf5d8TX6e1dFNwbfGF4wA6dQjVaDQDR0lZPV97bCP/8/bBYMalVRey6q5LtRUShVhJ/+0XznVmOyj/6PsHqsmIEBLLIskPLo087yr2T7sRx+aHmVL7PjMnjg6+ZgnJ0Vo1WXFv4Zv8su6HCCSetWsoc/+C/KQ2ZuR+GqfNYHdVA0vmjDvrPWs0mt5T4nLz900uXtlRFRFQh+FXF9Jb4lpI+qKOpLC0joCU1JOGTxoYQmI0nVRWiRUrAcicyMy8T8GsOfD7qfD2vXDVfZA7n2SnnbEZiZTWKouk/aP/wy6Uz9UhvUzceB8QwGFz8lDq1Wwq/SwwK2Ifb+0+zpGaZhZOzgnuSxcwajR9w8pCF/e9ujuqN2A41oX0lrgWkr5g4eQROO0G5/sPYhBAADz7Zbj1NRg3Fwy7apWSmafe4D4KwoCjH8LTReq83Pl8OuEIUw6soXDVBC5s2YJEJX4ZSIRlQPs9fM7/Jle0vcu+oonMmn8NAK9sq+L7z283e/ccQgL+gC5g1GjOlRKXmxWFLl7aFln/rNN6e0dcC0lfuLYK8rJYccdCPB9sQBxVP2D4Par7b+58SBkJzSfhVAVkjFfrMiQMlG1g/4lGft1wD068iP1KQLzSwE4AM1ELUC1WDFSwft3rqzk9uoCCvCze2HUcULETr19GrZTXbVY0mt5R4nJz07JNeP2RZohO6+09OtjeAwrysrj0qi8h7ImqNYrNqVrIVxQpEZEBePZL6jh/UaiNimFA/iLce9/DiTdCNIyO427Gq3ZgAan6dW30zwzOfvebhYvWtyQLK024xOXma49t4ndvH+CWJwp1ZbxGcwbeP1DNt57d0klEdFrv2RHXQtKn5M6HW1+Hq+4Nuqso2xB6PdxK+earYE+AKVdD7nyyZl6OIJSN5ZMCHzak4VTCZE+E635Le1o+R+U4bvHcw04xI1hLUuFuBeCbE07wwbzNLEosxWETrLh9AQV5WRSW1uEzL25ZKRqNpjMlLjff/Ntm/vXJYmqbVRdvndZ77sS1a6vPyZ0fHGQFmNZHghIRy0oByLsUpl4LJ3YCMHNkIgg4OOrznBo5j8DpOrJmX8XMMekRA7ISc/LIdjaytWo63/70JArysmhq83K4ppm54iA/q/s1tlofj0snNwfuYeroawEiihftNl3MqNFE428bSvn1G/sifAEGw6vde38R10LSFzGSbsmdr6yTDtMSAZh8BexfC2/fB6fKwZ7EjDuWgzO58zUskrPJaDgGQHqy6jq8q6oBKWGhbR+G9AHgQM08Od7QSkaSg6mjUoOX+M4VU3Sbeo0mjBKXm4ffO8j7B2oj1gXgdBhaRPqAuBaSvkj/PSMdrRSLJPMHc9PDgISMiXByd/Rzg+/JRrTV47QZNLR6AdhecQqAg4kXgP95ACSCwsAs5jW0MXNMOuV1p4OXaPNGtnKw+gR5/QGd5aUZVpS43Dxf7OLFkqpOab06oN63xLWQDCpul/nE/AluKIenvxiKr0QjORvR6iYzyUbDaSUkjTte58eJe2hOyEWchoDhoMGfyFY5jRMNbQC46luYKw5ymWM/zWWXAjODlywsrQv2CdLNIDXDhZWby/mvV3bj71hZSPzPTx8MhqSQmDPbfwXsAZ6TUn4wqBuKxqRFKojuawci04G7FJKkbJABxiV5aWj1sr94PT+u/yUg8bXYQED1Bd9mzI4/8wv701BxI8yfSFtpISudD+AUPjwnX4KK2cHPuCRMNHQzSE28o6yQclZvqeyYF6nrQvqRwZhHshxYDFRLKc8PW78OeAg1j+QJKeVvu7mMBJqBRNR899jDyvLasRK2rVRFi+EB+Wgkq9byE5yt1J324t7zLgKJEOCQfppJpCbjAsYA37S9g2/3hzA/j5Tjm0gQXgRglz5aDn5AiikkYzISARidnsBfbinQv0CauMUaFtdxzo92Y/U/g2GRPAU8DDxjLQghbMAjwLUoYSgWQryGEpUHO7z/NmCDlPJDIcRo4I/ALQOw795jxU/m3Bw9IN+RJCUkY52nKW31kjJ9HqJMvSQEuEknreJdpFQpizbphbINbPZN5TpCM+OfrBrPpS43BXlZVJqpw3bD0L9Emrhl05Fafrh6RycR0W6sgWEwBlt9JITI77A8HzgspSwFEEI8B1wvpXwQZb10hRtI6I999ildBeQ7Yloko+wtNLR6yRiTD0CTkU5aoJHxsppA2ct4seHAj0RCQyWJDU5VpyLs3NJ+D9v2ZpJwSE1ZrHSrQHxdSzsAm0vreP9ANdfOHqOFRRMX/PbNfTz6YWnwWLUS0m6sgSRWYiTjgYqw40pgQVcnCyG+AnwOyERZN9HOuQu4C2DixCHybcTM9BppU0LSVKf6/7TmX01q6cuqMl76ec5/BVelHGVceylyy1P8QAoQYEgfu+RkZZmYgfU2r5qR0uYNsPFwLd9cXoQvIHnykzJW3rlQt1fRDFkKS2u5//W97D3eFFyz2r3rlN6BJVaEpFdIKV8CXjrDOcuEEMeBJU6ns2BgdnaOmBZJtmimud1HS52qKfFPvhp/2Zvg9xAw7LwiP8OMcdMYW/ooQgSwhV0i02ilJpCGwyxM/HuhK/jaBwdrghXwnrAK+KWPF+LxBUi0G6y4U6cHa2KbEpebJzaU8tbuE3TMydLt3geHWGmRUgXkhh1PMNfOib7qtTVgJGSAMMhEzS2xhCR51mfZeNly/ui7gQdG/JYjCbNpz/0MPlNCBNAu1XeCH3xGzUH58XUzKMjLosqMkQBMGZkS6vclVEV8YWkdXjM92OPX7VU0sUuJy81/vrSTGx/dxJtRRES3ex88YkVIioFpQohJQggncBPw2rletK9G7Q4YhgFJWaTJRgA8DSfwSBvpWSPxj5vHX/zX80HrJNKTHCROvpRven5Cu7QjBOwMqKmNX5qVitNmcOyUqjGpdJ8OjvodlZZI3ghVWS8lvLn7OFnJTgxTXWyG0OnBmphk5eZybnxsE6uKKiJqQwS6T1YsMBjpv6uAK4AcIUQl8HMp5d+EEHcDb6MytZZLKfcM9N5igqRsUvxKSGiupl5kMcYwSE9S/1XHT7UxY0wadS0e2nEGuwjPMY4AkOhvYdbYNF7aVsW1s0dzorGNL84ZR2lNC3UtHhpbfcG4yN82HCXBYTBtVCr7TzRx7ezR+tucJqYocbl5saSC54orIqYW6rqQ2GIwsraWdrG+Dlg3wNuJPZKzSfQpCyrZU0ujM5sxQFqi6r3l8QdIS7RzuLqJhca+4FAsu1CPRyqOs+dYBr6A5Bt/KyIg4cIJmbyy/RgnGlqpa/Ewa2w6QDAoX9esMroaW30D+3fVaLrBau9jdWaw0HUhsUesuLb6hSEXIwFIyibBo/pr5dDAaadyNaUlhjQ/PdHBwsk5bDXOw4sdnzTApoTGVXWMgPnVzcqpnzkmDafNYO9xZenMyc3AZrqzHHaDlnaV2XWoOpT90pESl5tH3j+sZ51oBoTisnp+tHpHJxHR80JikyGZtdVT+r37b3+QnI39+A4ARopTHE28EAhZJADpSXYK8rL40R3f5I1tuVxq28v4WQvh7//MtMwA8+2HmRvYQ71MJVs0IyoF2SkJ7DmmhGTh5BG0evws/6SM//nqhfy/VdvJSU3gZGO7ql9JckRsqcTl5sbHNhEISBIcuvGjpn/587uH+OM/DgaD6bouJPaJayEZkO6/fU1SFkarG4MA2TRyOHkkAClOG4ZQw7HSTVEpyMuiIO8rwFcgoKyKXI+Llc4XEP52kMp91f7+y8xNup91dSoQOTYjkStnjmL5J2XUNKnhPpdPH8marZUcrm6iIC87YktrSirxdxicpX+ZNX3NlrJ6fvvmfraEWb26LmRoENeurSGXtQWqA7CvlakJp7AJiZGm0nmFEEGrJNw6CWLYwJkGdQcx/B6sLF9DqBnwF/l3AzBXHCR391+Y5dsPwIcHawC4YoYSrEMnmztd+uDJkMtLN37U9DUlLjffe24bNzy2KUJEQNeFDBW0RRJrtCnRW2IvBj84MsYEX0pLtNPQ6g1mcHUiMR0SM8CwQ8AbXPZh53jWJcxtPsjzzl/h+MiP057IpxPuY3Op+i6xYFI2TpvgxZJKpo1OI6W6BPfe9yhNvZgtrmRSnDZaPH6W3zpP/1Jr+oyVhS7ue3V3p3khoOtChhJxLSRDjooiKPwrAN/x/x2A5KxxwZeVJdIadG11IiFddRie/UXY+yp1+deTU7qGqqsegsY5LKx6CLtQLjDh9/K5lEN8Uj+ZBLtBRf1pvAHJFpeb/378aZ6x/ZrpeLkIB2vEvezyzwAgPamLz9ZoekGJy83fN7l4eXtk3bFO6x2axLWQDLlge9kG1W4esKFu+GkjxwdfTjczt8IzuCJIzID2RrCPgoxcchbcAKVrmDJlBiMOOtkaCGseYNhxj5oP9TA+K4nCo/XBPP1L5F4c+DAEOKWPhcY+tgemA3Ckppnzxw+hLDhNTFHicrNmawUvFFcG2/VY6LTeoUtcx0iGXPpv/iKwJQAiGOOorzoSfNmKjXRpFSSmK9dY0wlIGwvJOWr9dC3ZKQmMMk4Fr8s1P8eRvxCA8ZlJLJw8ggS7+nEoDMwiYJ7px8ZmOQunTb1WWtPSV39bzTDDqgtZubmik4jotN6hTVwLyZAjdz7c+ho1025EStXGZOr629lfvB4IWSTdurbaGqHpOKSNhhRTSFpqyU5x8lXjQ7xWi8fkEUwdmcpccZBr61aQUl3CyjsXctXMkWyV09kRmALAztFf4qprl7DizoVMzE7mSE2zrinR9JoSl5t7X94VURei25vED3Ht2hqS5M7nsO8VshHYhMQhfbj3vgfzrqHVbAlfXt/C7HHpnd9rubZ87coiSQlZJDb3RxQYh5T7SkD1vk9IHpnEaucvoQU8a1fhWryKgrx8PjhQwxhDiYThTOa7VyrX4OSRKeyuamDpsk14/JJEXVOi6QGPbyjlN2/s03UhcUxcWyRDMv0XyJp9FR4c+KSBFztZs6+ixOXmH3tPAvC957ZHtwYS0+F0PXiaIW0MOFOVq6ylhlEVbyKE+iWWErzHdpFa+hY2IZVgoQRr4eQRZNnbmCBqAchLCKUDTxmZSqW7FY9f3RLavbpbsKZrSlxuvvm3zTwQJiJWXciquy7Vbqw4Iq6FZMjFSExmzrsG1+JVFE/+Dq7Fq5g57xoKS+uCrU98XbV7T0gH61c2bayaz5uSAy115IxVbgOfVN25UtPSyR6v3FcBSVCwCvKyWLEk9O81gpAITx6ZEuHblqjuwtrFpQmnxOXmR6u3c8OjG/noUG3Ea7ouJD7Rrq0YZea8a2DeNcHjhZNH4LQbeH2BrosCE8PcXWlm/UnyCDhdy/hcJRol+Xdx4an1ZCQ5yBitihDrU6dSe8V/q88EZhrmsMqRs6ClOnhJnz+yfbcEniuq4OVtVdrFpQFUu/f/emV3RKt3C10XEr8MSSERQhjAr4B0YIuU8ulB3lK/U5CXxYo7FnY/EjcxM/Q8bax6TMmBllo45YKUkSz41/+BFTdC8wloUIKRkzOGnDDRonqfcouNL4Aj7waXT532dPrI8LG++gYxfClxuXm+2MXqLVURA6d0XcjwYDDmkSwHFgPVUsrzw9avAx5CzSN5Qkr5224ucz1qimIdar77sED11urmFzEhzCJJVa1VSBkJdYfhVDlk5oXWTuyCU6blEWZ1AFBeqAL3AR+01EAgAIbBZdNG8tcPj+D1BbAZAm9AImVk25Ro89/1TPj4psTlZumywmC3aQtdFzJ8GAyL5CngYeAZa0EIYQMeAa5FCUOxEOI1lKg82OH9twEzgI1SyseEEC8C76IJubYcKZCQpp4nqxgJhh3GXqTWUnKUQDRYQlITukbph3BCdR9m9xolJm2nIDm7k1X0x38cZO+xBp4w26aUuNx87bFN+MO6BAPctGwTXp3lFZdsPFzLD1bv6CQilhtLp/QODwZjsNVHQoj8DsvzgcNSylIAIcRzwPVSygdR1ksE5mRFy8/i77/dDjESzSB52hgVaAdIGQHeFnC7YNYX1VrqKNWL66Rq5EirG/xeOLYNXv5W6HpmR2GaqyFZdQQOt4omZidz4EQTBXlZ7C9ez+73XuVCOZmtTMfjCyUEeP26c3A88ps39rFsQ2nwWKf1Dl9iJUYyHqgIO64EFnRz/kvAn4UQi4CPop0ghLgLuAtg4sRh8q3Icm1Z8REIVbdLP2Sa/w4pKshOWwMkZSkhOfQOrL5VCQqAMJQV4/eYrq+ZnT4uPdFOU5uX/cXryV97E9PwcaPTwS2ee9jBjKC7S5gpx4aeCR8XbDxcy/1r97L/RKgrtG73PryJFSHpFVLK08DtZzhnmRDiOLDE6XQWDMzOBhnLtdV2SjWAzJ0fKkqEyBiJxfgCOLweSj8IiQgGTL4C5iyFl+6MdH2FkZpgp90XoH7Pe0zDZxZQerk2+RBGKzSv/4jRF15DRpKDU6e9LL5wnL7JDGFKXG4e+/AI75j1TOHotN7hTazUkVQBYR0FmWCunRNDtY7krDm513zcA09/UYlJcriQdLBIAMbNNV/LNxcE2BPgiv+EKVerpeYwIakogg1/gIqiYPPIhGmXEzB/lPzYmD4pjxccv+Sy8r+St3Ypk1v3AHC8obUP/7KagaLE5eYnL+7ghkc3RhURndariRWLpBiYJoSYhBKQm4Cbz/WiQ67777ni+hj13SCgXFJlG2D2l0KvZ5paHWGRmEJiznxnypVKRHLnq2wtYQtldbkK4enFIANgc5K34AlAMHLmIrZvvpp5je9QU/AfpNbVYQgVF3GY3YNLE85jZ2UD/oCkxFVPcZmbhZNHBOeeZM2+KljHookdVm4u579e3R2ckGmh03o14QxG+u8q4Aogxwya/1xK+TchxN3A26hMreVSyj0DvbchT/4iZU34PWouSf6ikGsrdTQ4ktTz5BEESwoti6R8o3q85DYlIgCGod7fbArJjpWhgVl+D+MbSoBLaGzzErAnAjAuS7MDAAAgAElEQVRh/ATaMyaBy2zFgp3CwCw+f/4YVhVV8IvX9vBsoQtDwCW2Qzxr+xXT8OMpfZx36v7GoYTZOk04BihxuVm9pZzniyt1XYjmjAxG1tbSLtbXAesGeDvxhdk9mLINSkRy56u7ubCpwLkVN7HZVXqw3wPuMrAngcsUkhHTIq+ZMioUIzGsHxcBNidtEy6FLV6a2nwkeMxWKqfrmDLlcnhfFSuunvEntu7M4gcXjuNA8bukFb/KXDGLrXI6l8i9OPEhhLJcdny8lkd8NhLsBivv1GnCg0WJy83Sxwvx+HRdiKZnxIprq18YkqN2z5Xc+SGLAqCyWLmiGqtU3OTW19R6exMg4ZnrVZC++aQSnOxJkdezO+H4TiVCTSfU2sQFcO2vEGIG8DHN7T7SvCEh4bTqr2QIKG2yYTMEWXXbeNH5S0DSjpNbPPewVcwCQpbLJr869vp1mvBgUXS0jh+u3tlJRHRdiKY74lpIhl2MJBplG0L5t1bcJBy/BwxzRklWvnKNWVQUwbHtKnX46S+GXkufALnzSatVQ66a2rxM8pmpoKfrIrK8Wqp2MSb9Glp3vxIRN7k84QAXfu42xNvglQYbFj7Bjg0OkBKbThMeFP74zgH+773DwWNdF6LpKXEtJMPSIumINXUxPG4CYE8MrWVMhMZjkDM98r1lG5Q1A+BvB5+ZddWsMnesrK3mdh/JAbPdfJiQ+KVgUqACV2YSY7PToSJkfWzhPK5oKwPAIQIsWPhpFlYe5pMjdSydP/Gsblq6FcvZUXS0jgfW7mNHVajTs64LiRMqiiJd3f1EXAuJtkiIHjeByLVtz0JFIeR0iI9YwXtfW0hQsiYFhSTVFJKmNh+pMtwiqQVHMhX+HKYFKjmclcQ4lLurzZHBs5P+hw270rjdVRz8qFMnXMEKeMOqyu8Fm47UcvPjmwGC7Vn0DbB7Slxunvy4lDd2naBjr15dFxIHVBTBU4tVgowtQf3O95OYxEodSb8w7OpIuiJ3Piz6QeQPUfiaVYho2Du/79bXYZxZz2k4YPR50KSEJMFuw2k3aGhpJZ3T6hxLSJJzqEuezAxRwfjMRNXDC0iyC/IuugKA1uP7gh/VVFvJicY2oPt6k/3F69n09D3B8cMW7+2vRhLZjVjTNX8vdHHDoxtZG0VEdF1InHBgnfIkyEB0t3YfEtcWiaYHVBTBrhfV802PwIzPdxacS/8N1tyuvtkcfEs1cvS2QlUJ33e8TODYxQAEhB3jdL1ybaXk4GtNZaJRw0LXo9B0DJKyodXN5GwVa8luLaPOyGAEDbTWV3Gy0QnAiYa2qFvdX7yeKWtvxCCAp/Rx9rMqWHuSPyIleF6X81o0lLjcPL3xKK/tOB6xrtN64wzXJti+Sj0Xtki3dj8Q10KiXVs9oGyDCqaDEoiyDZ3N31MugnUnAdPFdegdWP0vfFsG8B17CYCW5PGktbig8RhNJHNx/Vsg4FNVT6q3tzcAkonJHoSAyeIYR9MKGNH0Ps11VbT7VAuXY6aQ7C9eH1Gs6N7zLnb8wXRha5Y9QHaKEqH8nBT+cMMcfSPsQInLzZqSCp7fUtmpuFCn9cYZFUWqcDjgUyJScKtqd6RjJGeHDrb3gK6C8R3PsYLzhk09HnkfpGqMYsMHQFv6JCUk9UdoTZpMUsfGzOYNLMHTwKxMHyNbGzkyeg4tjZtor1cdcaaOSuVITTN7Ct9h5ptfA2TQ+kgcdx6iLHI0sEVtczsAiXZD3ww7UOJyc/PjhbR3SOkFndYbl2z4oxIRi4wJ/SoiEOdCoukBXQXjuzondSy8+p1gyrAEfNKGTfjwZ02B4x+A34Nt1Ay8TaUgvdiQSAyEzewmfLqOL9hUoF0GJPVGNqJZuVouys3kcHUztcWrsQl147Osj9RJ8wCoYiT3iO8xvymPFpebgrwsapqUkFSbjxrFlrJ6frJmZ4SIaDdWHHPwHTj4ZujYsPerS8tCC4mmcxFjd+dY7VJqDwHqprTCfzW32d/GGDElePqIqfPYP+ubuPe+x9ix48lPalNtWl79LpU73+fOpr+AgIuP/IVKey5ZXjcAF0/M5MWSSho9KnMr3Po4dbgIAHtqDhvqJvPxOwd5xHGYFXcspMa0SOpbPLT7/CTYbX34DzQ0eWj9Qf60/lDwWNeFxDmuQnglbJ4QAi6+ud+tEYhzIdExkn4geYTyu5ZvCi5lCFWY6Bgdlj6cnMPMi64JxjCA4GjfdlcRdtPtZcePkAFGo4RkzgQ1dz7QeAwE1IsMar+wnJnzrmHz5hUAJARUVld4hlZNU2iefE1TOxOykvv27z2E2FxaxwNv7GOnrgsZPgTjImHzhGwJMOece9/2CJ3+q+kdhk11D/Z71A8qcJ5wAZA8dkbovPAOwxbmlMW0EeORGMHixPqUyYwSp8hMsjMpR2VfFaBSg+34g5lZKU1qGl+q4cEwS02sDK2a5nbs5uLJxuHp3ipxufmP57dz07LCCBEBXRcS9+x4PiQi1jyhfqwb6cgZhUQIYRNC/H4gNqMZIqSNVo+TrwBgqqjitEwgIXNc6JzwgVoWjmSwJzJqRDaezMk0JE3AtXgV/tFzSBIerkvaR0rRQ3wpsYQJopZqssmkmbbTqmp+VHs5AM5AK3csmgzA7756IQV5WdQ2tTN1VCoANU3R04fjmRWFLm58dBMvbavSdSHDjcPvws6wVF9rntAAiQj0wLUlpfQLIS4biM30FHPE7i2o/c+WUn5qkLc0vEgdA+yA3Pm0lm0myXuKelJItjkgMVNNaIwmJEKoWpLT9SR56kk673oy513DlpPK0vjV6V/Duz5+j7IsytMuYlTTe9QeKyNj5DhGUY9f2LF5WrhlwUSWfVRKfYsXKSU1ze1cOmUE+080DSuLpKuphTqgPkwo3wwrvqqKDg0HzP1Gv6f6RqOnMZJtQojXgNVAi7UopXyptx8ohFgOLAaqpZTnh61fBzyEmkfyhJTyt11dQ0q5AdgghPgSaiiWZiARpiFrS+B08niSGk7RYqSpteQRSkiSowiJ9bq7DFrrg52Gk7LGA2DHaz5KpIQLGj8CAQ3V5ZxurCMNaM6YQcapPeRlJnB1ShneD97koHsGd8hdTE+8jleMRE42Rlok8diDS9eFaKgogtfuDrUvkoEBSfWNRk+FJBGoA64KW5NAr4UEeAp4GHjGWhBC2IBHgGuBSqDYFC4b8GCH998mpTRTh7iZM8xu1/QxFUVw+B/q+fu/JpA5B4BWuzkv3uZUsZMTO6P/QCdnQ1WJep6Vr97rNqusJUG3jBBgMwslW+sqaDxxAIDTiaPJYA8HNq/jL76f4fT5oAh+aIfAtld5I+XnnGwcHyxmPJR4EffvSMUfkHHTg0vXhWiCfbT8pvUtjH6vXu+OHgmJlPJf++oDpZQfCSHyOyzPBw5LKUsBhBDPAddLKR9EWS+dEEJMBBqktLoFagaEiI7AXhxC3ezb7Wnqh7v2gHrdmn3SUUySs8Fjdgo2hcRbfVDN3xLglyAxQIIPO3Y8+MoKuaT2FSSQc/xDENC8710c5lAsKVVqqwh4WZS0n91VDqbs+Qk2/FyEk1fkPWxlOm3eAA+8sZd7vzB7SIvJXz84rOtChjslz4RExAquD3BcJJweZW0JISYIIV4WQlSbf9YIISb04T7GAxVhx5XmWnfcDjzZ1YtCiLuEEFuEEFtqamq6Ok3TW6xKeLN/T+vIiwDwOjMjm8J11SQuOawHlikkWRd8ljac+KSBBycl599H8eTvULZ4Fc0yiYyGfdgJIAQYqBto5uQCAmbmVxBDcCyjgCmnNmLHjyHAgZoZb7G1/BS3PFFIicsd9a/3jz0n+fN7h7p8fTBZt/M4V/7+fdbvqw6u2W2CmxdMZNVdl/LAly/QIjIcqCiC7StCxzbHoIoI9Ny19SSwErjBPP66uXZtf2yqJ0gpf36G15cJIY4DS5xOZ8EAbSv+6VAJ7z92DPY+QWqLi7LWy8g/U7uVpGzzMQsSVVr2zHnXsJ9Vwb5aC8JqT1xv5iAJFSf6sWHDx5SZczh1+DPYqop4wLOU79vXkJOaREFgDwG/C2GPnBlviGCHFjy+0ATG/UXrObXzTTLnfJ6jiefxnRVbERAsdBzMG7MV28lKcvDy9iqKyyLFTQA3XpLLA1++YHA2qBkcNj0CWBbpwBUddkdPhWSklDL82/9TQojv9+E+qoDcsOMJ5pomFgmrhK/btYwJwEzPHto/+QVln/6FqmLvqt2KZZGY1ojFzHkdihdNGh05TG3bgxCwJf0aRs/6NBM3/xI8LWSmpVKbPJbn2q5mSuAkdzav5ermxxBmUbsQsGHBY1ydeCFfTXZy/+t7aPMFEEJNYNxfvJ4pb9yIHT9tFc/wWNoDwPiIQsfBEpKVm8v5r1d3dwqkWwjU3JV/ntuXjgFNzLPrRdj3OiBCcZEBKjrsjp4KSZ0Q4uuAmazMUlTwva8oBqYJISahBOQmVCD9nNBNG/uf0ycPB2MUDunj+PEq8m/9TddvMIsSOwpJV7QljiapfRsAU7/xMFltFbAZ8JyGtgYSUrOgHny2BCSmr9a0PgwBn/3MZ/ismYo8Y0waP1y9g1aPj7kTMyl8Z21EN+Fxp7ZgeVQHsxV9cVk99768q1M9COh4yLDm6AZYcwcglYBc/PVBSfWNRk8r228DbgROAMeBrwJnFYAXQqwCNgEzhBCVQojbpZQ+4G7gbWAf8IKUcs/ZXL/DZy0RQixraGg488masyLr/FB8o2NH3qhYQtJcrXy9Z8CXOhaAo0YeWSPHgtOcO+JtgfZG/E6VdvwPzwX4pC0yZgLQeir4tCAvi3/9dD7jmnbx/hM/JevEJ1jDGH3Y2OSfhTWb8aGbLh6wm3SJy809L+3kxy/uMOeFlHUSEUOAU8dDhi8VRfDytwjmNQb8g5bqG40zWiRmau5XpJRf7IsPlFIu7WJ9HbCuLz5DM3B0jG/MjOKeiqDphHp0bew6syscMzPFnTiRSaCq4wE8LdDWQLVjLALYKqfzM++/cL/jKYQMEMCGEx+0RX6JmO7Zyyrnr3FU+qwJKwgBj/v+ia1yOoaZBeb1R6bW9kUtSrRrbC6tY+njhcH4zctbq8hMcSBQ1odhCO64bBJpSY64qoPR9IKKInjqCyr2CIOe6huNnla2LwX+NAD76VO0a2tg6Cq+EZXG4wSHZFmZXV0Iyf7i9RQcfwEEnN+yif3F65k560L1oqcF2hrJHJ1DwnEDry/Ay7ZrufLSK0g7WcjYMePI33QvtEUGqGXZxzjD0oYti2SU0YAtoFxa/oBkW/kpFl+oWr6s3FzOfa/sIiAhwW6w8s6eB+GDAXMzRtPuC2C3CW40XVNPfHyU8DCINyCpafLwtXm5TMxO1uKhgaJlIRGJgVTfaPQ0RvKJEOJh4HkiK9u39suu+gjd/TcGmXw5bPhD95ldJu6972GYXYINAmoi4pxL1Yse5doalTOKFXcs7PBN/4tQvR823dvJIsmadSWUPkJAKjnzS4FNSKbPmsN/jJ3Bwskj+Pvq1eRsW8v+0V+hZVQB976yK+gy8/QgCN/U5qOmuZ3fPrOF9/ZX4w9IbIbAZyqG1y9ZsbmcNVsrSbBH9y6/ur1q0LPGNDFA2Sewa3XoOAZSfaPRUyG5yHy8P2xNElnpHnNoiyQG6ckgLZOs2VfhLX0cpC8Uf7EnqRdP16opcInpFORldb7hmqnF4TESgJnnF8A6KMteRPXYK3E2HOXiyme5ePYsLr5oKmVv/ZnfNf0MgwDta1fzcO4fkDIUdBcCFk4eEeGmAoLPW71+7MeUeIX3v/JFyb5q8wZo8wb4l0/l4fFL9lQ1sLOyISayxjQxQDAuYhEbqb7R6EmMxAD+KqV8YQD206doiyRG6ckgLbqJvziSQ7GWhPTob05Sc01oixQSmtXNfdJVtzHp/K9ASx387llob4SyT8grvC/o7nJKL9k1RcDnMVCZ+1fNHAXA1x7bhD8gcdhUNMMXCOC0G1ycm8n3uvk7GQIMISKE5bniClbcsZB/njuBW54oxOsLDGrWmCYGiBoXGbj5Ir3ljFlbUsoA8OMB2Eufo+eRDH1mzruGS2/9TWQQ35lixloIWR4dsScq11kH15YlJKQqQSBBtZ6nvREOvRPM2gKQCN5qnsrFuZn84HMzSHbaGJ+ZRGFpHb6AVJaDX+LxBwhI5fZy1Z2O3IZNcN64kNgJ4MZ5uSyalhP8rHDrY8UdC/mPz87Qbq3hzv43OsdFBnC+SG/pqWtrvRDih3SOkdT3y640mu5wJEPTGYRECNXSvrWjRWK2F0k1Z6rYE9Q3vfYmGDkLMDO5gNf9l1Lsn4bzeCP3LZ7Nmq2V1DZ7uO78MaGt2JR1YRkYJxvbGJGagM0Q3LJAdeAFIiwNq4iwuKy+k/UR1U2nGV64NsGO59RzsxVRLMZFwumpkHzNfPxu2JoEJvftdvoW7dqKU5wpqhU9dO3aAuXeOpNFApCQpoQkXWVpiYu/gX/bSo6hbu5+v7IYclISqGluZ1ymitOkOG08c/sC/vD2ATaW1gXFJDPJwci0hIjWJZ0TAqKvaYY5wZG5PjWNdO6tMVN02B097f47qb830h/oYHuc4kwBr+lCSuxGSBIzosdI7ImRAmQJSauZKrzg2/gOvU9uYz02Eapy33Osgf0nmjjRoOadtHr9XJSbiccfYNqoVA5Vq67GR2tbSHTYSAv72GiWhrY+NJ3Y8AclIqC+qsdQ0WF3dBsjEUL8OOz5DR1e66YPhkbTj1hFidC9RdKVayt1VKiABJQYtTWqYVsAydkkjMjjyrGeiHhFTmoC45p2kbbl/5grDhKQUNfSzsmmNhw2EZwjH5CSxjYvGk2vOPg2HHwrdGzYY6rosDvOFGy/Kez5f3Z47bo+3otG0zOsNinQdYwEunZtWfERi4T0SIskKQsyJpDWfoLvXjk1aDWcHzjAE/KXzNz7ECucv2GuOEh1YzsnG9uZNjoNp93AZmZlpSc6+uAvqhk2uArh5W+HLcRuqm80zuTaEl08j3Ycc+gYSZxiWSTCFikqHYnq2qqG7A6hvYQ0OFUBp+tVnYojCdLHQ+Mx1dPIUO2Ep7fuwIkXg9Cck8PVzXh8AS4Yn8E3L82nsLSOWe+nk5bY0/CjZthTUQRPfyHk0orxVN9onMkikV08j3Ycc+j03zjFEo+EtEgXVUcSTYskvJNj88nIQLt1nfZG5QazmkpmTFC/2M2hosL28QuC354kgsLALHZXKYtndHoiBXlZfPfKqVpENL1jx6qQiAyBVN9onOknfo4QohFlfSSZzzGPE/t1ZxpNV1hC0p1by3pdBpTbKjEd/F44XRfFtWUJSb1yawFkmONxGiqD2VwZqckIAe3STquRwnY5FUeYkGg0vebwetjxvHo+RFJ9o9GtRSKltEkp06WUaVJKu/ncOh40J7AQYqIQ4hUhxHIhxE8Hax+aQcJybXWXsQWdq9tbzJHLUS2SJuXaCgqJOTCqoTJ42tiajQSk4Je+W8mUDfwp8XEcx7YAMDo9IfKa7U0qA6cHrfI1w5TyzbDiBjUSwXBAwa1DzhKx6Ok8kj7DvPlXCyF2d1i/TghxQAhxuAficAHwopTyNuDiftusJjYJurbOZJFYQmIG3IM1JFGC7QEfNB0LExI14CooJBVFpO17jlI5hkOBcUhgCR/yOPczVxxkVFqYRdJ0HI7vgHfvh6eXaDHRdKaiCF67W1nMoB6HSKpvNAbDmfsU8DDwjLVgzjx5BDUDvhIoFkK8BtiABzu8/zagEHhRCHEb8OwA7FkTSwRdW2ewSCzXl2sTHHoHDKc6jubaAmiogilXhd6bkK6EZMtTsO4HGAEf+cLgettGgGDQ/TMJB0hymvN9j3wAdUdC1/a1h1rlVxT1qFmlJs6pKIKnFgdn7cTifJHeMuBCIqX8SAiR32F5PnBYSlkKIIR4DrheSvkgsLjjNcx2LT83r/Ui8GTHczRxTNC1dQaLxHJtvX2PsjiEaYA3Hos8z6pFkf6QRQIq8L5rdai+BBBCBe6lsCOkDx82jiSbRnH5ZngpSu1r/iIzM2eJEhZ7Atz6uhaT4UrJMyERidH5Ir1lwF1bXTAeqAg7rsQanh2dt4D/J4R4FCiLdoIQ4i4hxBYhxJaampo+26gmBgi6ts5kkZhCEvACUgkFqJt9uLspIawGPcnM2qooglPlESICIDF4yb+I8iv/F4An/Z+jLmtOqFtrS3XoZEu4sqcoS8TXTsRAr4oiHUcZbpRvhu0rQscxOl+ktwzJPEUp5W7U3PjuzlkmhDgOLHE6nQUDszPNgNBb11ZH/N7IyYwRQmJaJGUbOr/PsPNEyr+xtW06h0cUMNKZw1h/vcrYKnvXFCzrOpnwxYfg1X+Dtd+DKdeoepSAT1UsJ41Q7o2AV9UMDNEgq6aXbPozaiABDLWiw+6IFYukCsgNO55grmk0nbFcW8e2d/9tvnp/6LmwqcwYK8Uy3B8dLkhWHUn+InWDt86/5Db2f/55flenJjTe/dw2KjIKWGTsYl7Fclwt1ncyoSyRzDwYMVUd73sd3vopjD5fnXL5T6C1Trk3ZCBkoWjim52rYf861M+ITfV8G0JFh90RKxZJMTBNCDEJJSA3Aef8L6ybNsYpVuffw+uh7OOuv82Xf0JwPjzA3K+r+pCOwe5orq0okxzfff8wUh4A1AyR0vYMZoomvtb0DL5NdvVR538FxuxQ13R9HLqu3xOKzWRM6Fxd31CpRDEOvp1qolD6kRk/k+qLycVfHxJdfXvKYKT/rgI2ATOEEJVCiNullD7gbuBtYB/wgpRyTx981hIhxLKGhoYzn6wZOtSXmk9k99/m8xepb32WVTHnZlj0g86/vOGxlvBge+78iPMXTh4R7KflsBskm4laNiGxY1YmX/6TkDDlLwKb+V3NsIfiLd5WdU1nCmCo2E3JU/D0F3W8JB4Jjsw1v9AE/EM61Tcag5G1tbSL9XXAugHejmYoMuPzUPgXFevoLm2yp/Phwy0Sy7UVBWuCoTVDJOWkg8AbagK1RBAw7BjZUyI//5+XwwvfgBn/BHtfVus+1YYevw/GXgjHt0e6uOLoBjPsqSiCp/5J/axCXKT6RiNWXFv9gnZtxSm581X6bE9qMnoyH96eoH65/Z5QplcXRMwQybsG965rSKv8iPasaaQkOEIWiMWsJZA2Fo5+EFrztUEgoGIkoy9QQoKIyxvMsGfzoyERiZNU32jESrC9X9CurTimg9vpnElIA2ca2J29elvWwq9jlx5S6veGgunhCKH2aLWoB/C2haySnClKTDJydeZWvFH2MexeEzqOk1TfaMS1kOjuv5oeY3Mqt0NvYxSTLicY0B99XvRzcheox7SxKhPM1xoSEnsSjL9YTXwMv8FUFMHa78Paf9dxk6FIMC5iET+pvtGIayHRFommR1QUQdMJaG/ofcA7ORtypqnn5tySTkwwbx72ROX68rapgDuAIxFypsPpWtU00trPk9fBlidhy3JVb6LFZOhgFadafdqEEVepvtGIayHRFommR4RnffW2pqOiKJRF9o+fRb/hW+0w3GXgaVEV85aQ2JOUkADUHlKPe19VmT1nuyfN4LJvrfo/A4bqfJHeEtdCotH0iI5pwr0JeJdtCA3OsirmO1Kx2WyXYp7XWKncW2BaJKZFU2cKiejwaymEDsIPFVyb1KAqMIsOE+I2LhJOXGdt6VG7mh7R0zThaOQvCmV8BUXo9SjnJKhzZEC1v/eGxUgy89R7aw+qtVMuSM5RGV9VJXBip7JSrL1qYpOKInh6sdkGxwZzb42rosPuiGuLRLu2ND3mbLPALBG66t6u3Rfh54yYCs7kMIskSd10sqdA7WGVFnx0A0z7LCz5X7jyHnXepkd0wWKs89HvQyNzJXFXdNgdcW2RaDQDQk9qVaxzDv1DZWxZFokjST0mj4CKQuUWaa2HyZer9eq95gWkLliMZQ68BYfeDh0b9mHljtRCotEMJPZEFXC3LBJ7orIyKgrVt9lXv6vWrbYt+YvMrsF+XbAYq7gK4eVvhy3Ed6pvNOLataXTfzUxhz1RiYg3zLVVtiE0ctUKyL94W6iJo5U2+o2Xh9XNaUhQUQRPfwHazILTYZDqG424FhIdI9HEHI5ENeDKG2aRWMF4ROi88JTfMReoRytNuCusQVlbntIDswaK7StCcZFhkuobDe3a0mgGEntSZIsUR1IoGL9jJWxbqW5M4W4sqyNxqxtSRih/fPWeyAwzqwguWL8glEgNw5vagHFoPexUTTuDqePDINU3GkNSSIQQs4FfAHXAu1LKFwd3RxpND7EndHZtQSgYP+fmzmnI1oyUVjdsW6GmLgojcrJi2Yaw5oAQDM7vWHl2ac2a7infDCtvUC5JwwFzvzFsUn2jMeBCIoRYDiwGqqWU54etXwc8BNiAJ6SUv+3mMp8H/iyl3CCEeA3QQqIZGjg6WCT2xMjXo2WABS2SejXMCzq3nQ8G5X2h9xk208LR43z7lIoilRRhxbVkYFil+kZjMGIkTwHXhS8IIWzAIyiBmA0sFULMFkJcIIRY2+HPKOBZ4CYhxO+AEQO8f43m7LEnmum/rWY1vTjze5LM1vatbkgdHVoPd3/lzodpn1NjiKd9DhBw0c2hIkjdZqVvqChSvc/CuxDobLqBFxIp5UdAfYfl+cBhKWWplNIDPAdcL6XcJaVc3OFPtfnnu8BPgdoB/itoNGePI0lZCJ7mztZIV4THSCxXmOHonMXlTFZCM2sxIGHK1aFGksOsrqHfKHk61DttGAfXOxIrWVvjgYqw40pzLSpCiHwhxDLgGeB3XZxzlxBiixBiS01NTZ9uVqM5a+wJ6rH1VEgUzkRiBiCUkLRUq7WAN3Qti/ZmSEhVs01AdSa2xOPqn/XsZmdlfumMr86Ub4btK0PHcTxfpLcMyWC7lLIMuOsM5ywTQhwHljidzoIB2ZhGc5HEKHMAACAASURBVCbspni0untukRg2JSatbmiuUVZH80koL4Txc0PntTep4VyWkDRUQptZQ5U25syfY7lt/O1qn/qbdiQbHwKsep/hV3TYHbFikVQBuWHHE8w1jSa+cJjiEe6m6glJWWpeSfNJVVeSOga2Ph1pOXiaTIvENOZPVYD7qHre1oOi3LINofRhHVOJZMdzKu0aYXb1HX5Fh90RK0JSDEwTQkwSQjiBm4DXzvWiuiBRE3NYFklbL1xboISk1Q0tNYBQjzX7Yfl18Pr3lKC0N6uRwY4kSBkJJ3eFRvz2REiszC9QjzqmojjygWqBIv3KnVVwq7bWOjDgQiKEWAVsAmYIISqFELdLKX3A3cDbwD7gBSnlnj74LN0iRRNbWHGN0+6QqPSEpCyV/ttSA94Wgq1UpB9KnlKdgVvrwZmq1jNy1cxwi+6ExIqLAEw3EyoXfEvfKEHFRV66k+C/d8A/7FN9ozHgMRIp5dIu1tcB6wZ4OxrNwGJZIe0NITdXT0jKgmNblctpzIVQtc2sRbGGanlUDUlCmjrOmKDOt+hKSCqKYPnn1HAueyJMvFStW8O6rHOGY1Gj1S0gYBZ66lTfLokV11a/oF1bmpgjPMDe02A7qAwsy001/hLlWrnkX0J1KDanuuGFWyQWaeO6FpJgw0izEr7RnDPuLlOPFUXw9BJ491fw5D/Bq3cPn4yuTx4KiYhO9e2WuBYS7drSxBzhcZHexkgsUkepm9ni/4VJl6sWKkvNtFTLIsk0hSR1jMrYam+Mft2Jnwo9tzlD4lZvBunLNqgmk0h1U9327PAYsHVoPexfGzrWqb7dEtdCoi0STcwRXvvRG4uko5BYjDDHSOfMUI8JlkUyQT1mT1Kpw11ZJFZH4ZSR6tu2hfuocm/lXdb5PfGe0VW+GV75TtiCTvU9E3EtJNoi0cQc9j6wSFJGRa63nQpZHB1dW97TKkDclZCcrlOPhkPdKNubQu9rPgk2OyAh71PqHIjvjC5rvohV+DlM54v0lrgWEm2RaGKO8AD72QiJsEWKSlK2inE0mmVXlmurxewcdHwnuD4x04ajcNo8zxKa9ibInKie1x+Fj/+k2qss+iF8fQ0g4Pyvxu+38/1vhHVR1nGRnhLXQqLRxBwRwfazEJLUUWAYnddPmR2GLCE5sQP16y2V0FiWhoWV8lteqI69LeoG2t6kssIA9r4K+15X2WDP3aKEb8wFoYB8vFG+OXK+iD1Bx0V6yJBskdJThBBLgCVTp04d7K1oNIpwIelt+i+oWEa09VPl6tFybeUvUjdCv0dldgV8KmhuT1BV2i9/W60LW+haLTWqPcro89U386LHQq8FW9YvgB2rlLvMCHvvUCc81VfYVNHhMJ4v0lvi2iLRri1NzBHuzjobi8TTEpkx1VFIrGC7NXXxqnthwbfVWlsjBALw3q8JWirh80saTEvDsl6seRvh9RO581Xn4rf+M74ytzb+X1iqL7rosJfEtZBoNDGHYQsFrXtjkdSa8y/qj0Sm31pC0mC5ttJD78mdD4t+AGPnqOO2BvjgN6FzEZFWhbXeWEWw0LFjnMDmVMtFy+InDfjweuXCs9At93tNXLu2usPr9VJZWUlbW9tgb0UzyCQmJjJhwgQcDsfAfKA9ETze3lkk5RtDz8MnIyabY3itGInl2gon0bTIXZ/AR2FTF0bOgPTxcOTdyGuMvQgOvq0+p+Mc8roj5ptl5D6GKq5CeOlbYQs61fds+P/tnX14VOW16H+L8BFCACVClS+JECEEyAcQUBuKIIgUBRIsWDiKUiyVettePwqC4tVr9T7ghYPUcoLQHBCxAoro4SC1gRIqmkRBQKAGMJLwoTZABCRAkvf8sWdnJiGTzGRmMpOwfs/DMzPv7P3uNS+TvWZ9vGs1akVSU4yksLCQ1q1b061bN8STLnVKo8QYQ1FREYWFhURHR9fPRZuFW5V6vcna6pZiKR775m7/Yg53dE88e8L6JV21Rwk4FYmtMGzKy6w4SdOWVh9527XVKcmyQKorixKdYsUQTFn15UIaUjkVO9XXdu9JE6slsab6ek2jdm3VFCMpKSkhKipKlchVjogQFRVVv5apbYl4o0hcYx6u6ahhTR3uLGNZI9V9n21FUuZSM0qaWDGB8/+CdjdZ47YiadHa6Rarqgy6JMPQWdbzUS9Wfj//H1YZlcwXGobba+9alxiRpvr6QqO2SGpDlYgCQfge2LERb3a2g3WDq+4m1/Iaa0Oia3zEFVuR2PWzbvut5aI6nGlZJR0T4NsvXBSJm3lskqfDtpesWIptgTSPtMbsgHWou72+/shSJGBZWFVdeIpXhLxFIiI3ichyEVnnMtZKRP5TRJaJyORgyucLYWFhJCQkVPzLz88Ptkhs27aNMWPGVDvetm1bEhMT6dmzJ0OGDOH99521iJYuXcrKlSsBOHjwIAkJCSQmJnL48GEWL15MbGwskyc32P8q/2K7n7yxSGrCDri3qCY+Ak7F8O0B69jhz0CXgZZ77fvjcM2N1s20uMpelJqu1yUZdq+x+qH87Tn47yetMvYASO1VcoPV0rcgG956wEr1vXDacgdqfxGfCahFIiIrgDHAt8aYPi7jo4B/B8KA14wxL7mbwxhzBJjmqkiAVGCdMeY9EfkLsDogHyDAtGzZkt27d3t9XmlpKU2b+ue/rqysjLAwz/YDpKSkVCiP3bt3M27cOFq2bMnw4cOZMWNGxXEbNmxgwoQJzJ07F4BXX32VDz/8kM6dO3t0HX9+vpDEdm15a5G4o6Uj4F5doB2geStnXKN9rOX+skuomDKIiHK08nUogtoUCUD7XnB0Z/XvXdMF0pbDpR9g64vQY3jlm7S9Z6Os1HLNJU6pnz0bBdmQMdpl5zpWPTFN9fWZQFskGcAo1wERCQP+CNwF9AbuE5HeItJXRN6v8q/DlVMCViteO4exLECyB4WSkhIefPBB+vbtS2JiIlu3bgUgIyODe+65h2HDhjF8+HBmzpzJxo1Wkb3x48fz0EMPAbBixQrmzJkDwLhx4+jfvz9xcXGkp6dXXCMyMpLHHnuM+Ph4du7cyebNm+nVqxdJSUm8/fbbHsmZkJDAM888w5IlSwB49tlnWbBgAZs2bWLRokX86U9/4vbbb2fGjBkcOXKEu+66i4ULF3L+/HkeeughkpOTSUxM5N1336328wHMnz+fgQMH0q9fP+bNmwdAfn4+sbGxTJ8+nbi4OEaOHMmFCxcAOHToEHfccQfx8fEkJSVx+PBht/MEFdu1VV8WiYjTvdXeUdzxGpcy87YiAesXuidyNW9VzXUct5MOcdbj6+Ph7y9dGS858ndHS99y6zH3z/UTU8ldUVmJeGI5KR4R0J99xpjtItKtynAycMhhaSAibwJjjTEvYlkvnlCIpUx2U4/uuU+/Ps3HR4oYfFMU/W+8tvYTauHChQskJCQAEB0dzTvvvMMf//hHRIS9e/dy8OBBRo4cyZdffgnAZ599xp49e2jXrh1vvvkmWVlZ3HPPPRw7dowTJ04AkJWVxaRJkwBLqbRr144LFy4wcOBA0tLSiIqK4vz58wwaNIiXX36ZkpISYmJiyMzMpEePHkycONFj+ZOSkpg/f36lsdGjRzNjxgwiIyN5/PHHAdi8eTNbt27luuuu46mnnmLYsGGsWLGCM2fOkJyczB133HHF59uyZQt5eXlkZ2djjOGee+5h+/btdO3alby8PNasWcOyZcv42c9+xvr165kyZQqTJ09m1qxZjB8/npKSEsrLy93OM2TIEN/+83yhaR1jJO6oUCQ1WBK2xdEh1nrdtqvzvVbXQXgb5xyexIx6j4Wc5dZu+SZN4JZfW3N8tsram1LR54Qr4yWuvVKAekklLsi2dvTbNGkGSf+mu9f9RDD8B51wWhNgKYVB7g4WkSjgBSBRRGY7FM7bwBIR+SnwnpvzHgYeBujatWt1h1Twf977gv3H3fRrcHC25DIHT56l3EATgV7Xt6Z1uPt9B707tmHe3XE1zlmda2vHjh08+uijAPTq1Ysbb7yxQpGMGDGCdu0sN0ZKSgqLFi1i//799O7dm9OnT3PixAl27tzJ4sWLAVi8eDHvvPMOAAUFBeTl5REVFUVYWBhpaWmAFc+Ijo4mJiYGgClTplSyXmrCuHbR85AtW7awceNGFixYAFgW2NGjR6/4fFu2bGHLli0kJiYCcO7cOfLy8ujatSvR0dEVCrh///7k5+dz9uxZjh07xvjx4wFrb0hN84SEIvG3RdK8FkUCTouk1XWWHKUllS0ST9xa4Mwiq5rqe3irVVHY9Vd+WDNoGWXFRLqlgK2neo+H/RsAE3jLYP+7ODdZCiRNgTELA3e9q4yQd0QbY4qAGVXGzgMP1nJeuoicAO5u3rx5f1/l+L6klHK7bbOxXtekSAJBq1ZOd0KnTp04c+YMmzdvZsiQIZw6dYq33nqLyMhIWrduzbZt2/jwww/ZuXMnERERDB06tCLFNTw83OO4SE3s2rWL2NhYr84xxrB+/Xp69uxZafyTTz6p9PmMMcyePZtf/vKXlY7Lz8+nRQvnXomwsLAK15a761U3T1BpVof035qozbXlymVHmrOIFRsoOlRFkdSSseVKdVlkEe2soH7HROfYLTNh8+8dtb7Coddo6zHtNbi2q9WJcMyiwFoGpfbndmRo6V4RvxIMRXIMcLVtOzvGgkZtlgNYbq3Jr33M5dJymjVtwr9PSvSLe6sqKSkprF69mmHDhvHll19y9OhRevbsyWeffXbFsYMHD2bRokVkZmZSVFTEhAkTmDBhAgDFxcVce+21REREcPDgQT7++ONqr9erVy/y8/M5fPgw3bt3Z82aNR7JuWfPHp5//nlee+01rz7fnXfeySuvvMIrr7yCiLBr164Ka6HqcU8//TSTJ08mMjKSY8eO1bjzvHXr1nTu3JkNGzYwbtw4Ll68SFlZmdt5OnRwF36rB0oc1u+JPdDVrTHuORG1BNsLsuHkHuv52qnODCVbaZz+2nuLxK0sUfDDKWcZe7DSju0ui2UX4din0KG3FWi/9X/BPxbD7lUQdVPglMmJzyEqBhLuaxibJRsYwUj/zQFiRCRaRJoDk4CNtZxTJ/xZtLH/jdey+heD+d8je7L6F4MDokQAHnnkEcrLy+nbty8TJ04kIyOj0i9wV1JSUigtLaVHjx4kJSVx6tQpUlIs98CoUaMoLS0lNjaWWbNmMXjw4GrnCA8PJz09nZ/+9KckJSXVeIPNysqqSP+dOXMmixcvrgiMe8rTTz/N5cuX6devH3FxcTz99NPVHjdy5Eh+/vOfc8stt9C3b18mTJjA2bNnqz3WZtWqVSxevJh+/fpx6623cvLkyTrNE1AKsuHQFuv5yrH+CTDXFiPJz3J6dexYREG2dXMF+MsUuHyh5jk8JSLKisWcO+kca97aGYhH4HyRVY4e4NQRyzrK3xG4gPv3x6EwB+InVb/JUvEZqYuf2+PJRdYAQ4HrgG+AecaY5SIyGliElf67whjzQoCub5dImZ6Xl1fpvQMHDnjtllEaL379Pgwdaj1u23ble1kvW9V3TbnlZhk2x7q5+cLRj2HFndBzNPz4d1feKAuyrZu0XV7Fjm24yhGdAke2QZ80mLCi7rLsfBU+mG2l/66fZo3F32dtADzztVWCpOwijF5gbWzMehn+9n+BcqAJDJ/r+3oAHN4Ge9dBWBj8cBoOvGvJ1HeC73NfRYjIp8aYAbUdF+isrfvcjG8CNgXy2ooSknRLcdxML/kvwHzW8ev/n5utYHfVzXXuAuOucnSIsxSJPywSsOIkYGVonfjcKnN/Q7zTCrKrDtt9U0ovAOXWjvuCbN+shr3rYf1DV46/+2ur+6NaJH4n5He2+4L2I1FCDnc1s3zh1GGH66jc6bqq7rqubp2qctzg6IroL0Xy3UHrsetg+HY/YKwsLRu7n4ktR+xYa3z3G767uD5aXP24u7VRfKZRKxIRuVtE0ouLi4MtiqI4cVcQsa7YVo6dkeSpleMqR12ytqojwhGv+fYANItwxkLAaudrU3bZeVPvkgwdHT1TXPeUuKOm8ioH3ndYPVX2wrg251L8Tsin//qCMeY94L0BAwZMD7YsihIw3LmuvMFWJIW5vrmWbIvk9FeWWyvK0cKh5bXQprOV9lt2+cqbercUa6zskuX2st+rWpa+INuqMFxeas3latUVZMPaB6jYlxIz0upxf308XCjSbK0A0qgVifZsV64a3FUG9hS7VW/eFvhqe93dbrYiMeUQ+SOr3hZYBRI3z4JR/6/6m3qXZLh/I7yeZqUB52fBN/th85MOxdPCqSzdVRg+ss1ZFr68zOqr4o/AvVIrjVqRqEWiKB5ypgDLHeRjuZLmkVb5kfLLljVw+ivne2WXLCXi7uZ+4y3QcxTsWw/ffOHomVLqPNe2TKSJI9usSWWrxi6GqW6seqdRx0hCnZMnTzJp0iS6d+9O//79GT16dEU5FG/JyMjg+PHjXp9nF1usbrxTp04kJCQQExNDamoq+/fvr3j/F7/4RcXrtWvXEhsby+233w7AfffdR79+/Vi4UEtQNBhu+onlKvI2zlIVEadV0qq91SyqaUvP523V3no05S5Np3Ce2yUZru1mjXXo7XR3Zb0M+duhWSurr4iWha9XGrVFEsquLWMM48eP54EHHuDNN61icp9//jnffPMNN998s9fzZWRk0KdPHzp27HjFe96Uinfld7/7XUXhxb/85S8MGzaMvXv30r59+0o72pcvX86yZcv48Y9/zMmTJ8nJyeHQoUMeX6fRl41vCPgjzmITEWVtSIzs4P28fdIg5zUXJeKwksYvdZ5b5njv5B7IuNslMG+g00D4yZN1l12pE43aIgnl9N+tW7fSrFmzSn084uPjK3ame1NCfd26deTm5jJ58mQSEhK4cOEC3bp14/e//z1JSUmsXbuWZcuWMXDgQOLj40lLS+OHH37wSt6JEycycuRI3njjDQCGDh1Kbm4uzz33HDt27GDatGk88cQTjBw5kmPHjpGQkEBWVhaHDx9m1KhR9O/fn5SUFA4etNJCp06dyowZMxg0aBBPPvlkjeXlU1NTGTVqFDExMTz5pPMmsXnzZpKSkoiPj6/YYe9uHsUD/JVNZpdssa0Lb+btkmzFUeysK3u/iWsv+pJiRzaYsawQDBVb949/Fvotfhsh+jPQG6pmkPjAvn376N+/+lqSdSmhvmTJEhYsWMCAAc5NqFFRURU1uoqKipg+3QoVzZ07l+XLl1dUGfaUpKSkCkVg88wzz5CZmVlx7ZkzZzJmzJiKqsbDhw9n6dKlxMTE8Mknn/DII4+QmZkJQGFhIR999BFhYWE1lpffvXs3u3btokWLFvTs2ZNHH32U8PBwpk+fzvbt24mOjubUKasp0wsvvFDtPK4FIZUAYyuSyDrWM7tYbLnIjLH+gZVO3PMuK4h+sdhZFuYKTGi3+G2kqCIB+O9ZcHJvzcdc/B6+2ecM8v2oT80599f3hbvcNn6sEW9LqLvDtbfIvn37mDt3LmfOnOHcuXPceeedXsvlbTmdc+fO8dFHH3HvvfdWjF28eLHi+b333lvhcqupvPzw4cOxrcrevXvz9ddfc/r0aYYMGUJ0dDRApfLz1c2j5XDqEdv15Fq40Ruq7v5v3sq5wfGio+Blh95QkGMdYysdjHWeBtnrnUatSPwaIykpdjbqMeXWax82b8XFxbFu3bpq3/NXCXXXX+FTp05lw4YNxMfHk5GRwbbq6kDVwq5duypZPLVRXl7ONddc47adcNWy8e7Ky1f9zKWlpbjD3TxKPVGQDXkfWM8/mG39oPLWOqgaV9n2krPkyoUz1uMN8ZWPAb95CxTvadSKxOP0X08sh6qF79Je8+kLO2zYMJ566inS09N5+OGHAas0e3Fxsdcl1MEqo15TVduzZ89yww03cPnyZVavXk2nTp28knf9+vVs2bKFl19+2eNz2rRpQ3R0NGvXruXee+/FGMOePXuIj4+/4lhPy8vbDB48mEceeYSvvvqqwrXVrl07r+dR/Ex+luV+AssyqaubyXVfTIdYyPmHNW+JQ5GEt71y74wqkKDRqIPtfsXPNZJEhHfeeYcPP/yQ7t27ExcXx+zZs7n++uvrVPrcDl7bwfaqPP/88wwaNIjbbruNXr16eSTjwoULK9J/X3/9dTIzM2nfvr1Xn3P16tUsX76c+Ph44uLi3Aa/PS0vb9O+fXvS09NJTU0lPj6+wo3n7TyKn+mW4p80Ylfa97QaU53OtzwBAOHX+D6v4jcCWkbeH4jITcAcoK0xZoK7sZoYMGCAyc3NrTSmZeQVV+qtjPzVgB+TUqz5cmD5HdBvIlx3M2Q+DzP+Adf38X1upUY8LSMfUItERFaIyLcisq/K+CgR+aeIHBKRWTXNYYw5YoyZVtuYoighgr+LUl46Zz3uecuKlwC0VIsklAh0jCQDWAKstAdEJAz4IzACKARyRGQjVpOrF6uc/5Ax5tsAy6goSihz3G4zbZwbFdW1FVIEurHVdhHpVmU4GThkjDkCICJvAmONMS8CYwIpj6IoDZBuKVTscG8SBuXlVkqwEjIEI9jeCShweV3oGKsWEYkSkaVAoojMdjdWzXkPi0iuiOR+9913fhRfUZR6pUsydBlkVRO+eZTV80Sk9vOUeiPk03+NMUXAjNrGqjkvXUROAHc3b968+i3kiqI0DK7vC98dsEqlhIdeyaOrnWBYJMeALi6vOzvGFEVRqqdtZyv1t7hQ4yMhSDAUSQ4QIyLRItIcmARsDMSFQrloo01hYSFjx44lJiaG7t2785vf/IZLly65Pf7MmTO8+uqrFa+PHz/OhAm1ZkB7hLuS8ooSdNp2th6/2a8ZWyFIoNN/1wA7gZ4iUigi04wxpcCvgQ+AA8BbxpgvAnT9kO7ZbowhNTWVcePGkZeXx5dffsm5c+eYM2eO23OqKpKOHTu6LbWiKI0GW5FcOquurRAk0Flb97kZ3wRsCuS1GwKZmZmEh4fz4IMPAlYdqYULFxIdHU10dDQffPABxcXFHDt2jClTpjBv3jxmzZrF4cOHSUhIYMSIERXVdvft20dGRgYbNmzg/Pnz5OXl8fjjj3Pp0iVWrVpFixYt2LRpE+3atWPZsmWkp6dz6dIlevTowapVq4iIiAjyaihKDdiKBNS1FYI06hIpoe7a+uKLL64oJd+mTRu6du1KaWkp2dnZrF+/nj179rB27Vpyc3N56aWX6N69O7t372b+/PlXzLlv3z7efvttcnJymDNnDhEREezatYtbbrmFlSut7Typqank5OTw+eefExsby/Lly+vl8ypKnYm83iq7AmqRhCAhn7XlCx5X//3tb8FNhdo6k5AAixb5NMWIESOIirLalqamprJjxw7GjRtX4zm33347rVu3pnXr1rRt25a7774bgL59+7Jnzx7APyXlFaVeCWsKbTpCcYHGSEIQtUiCSO/evfn0008rjX3//fccPXqUpk2bIlVy5au+rg7XkutNmjSpeN2kSZOK8utTp05lyZIl7N27l3nz5lFSUuLrR1GUwNPGsd1MLZKQQy0S8NlyqCvDhw9n1qxZrFy5kvvvv5+ysjIee+wxpk6dSkREBH/96185deoULVu2ZMOGDaxYsaLWcvGe4GtJeUUJCm07W1uZNUYScqhFEkTsUvJr164lJiaGm2++mfDwcP7whz8AkJycTFpaGv369SMtLY0BAwYQFRXFbbfdRp8+fXjiiSfqdN26lJRXlKAjjtvVuZPBlUO5gpAvI+8PGmIZ+YyMDHJzc1myZEmwRbkq0DLyIU5BNvx5NJRfttrpTn1fG1nVAyFRRj7YhPo+EkVRPCQ/C4yj82K5o/OiEjI0akUS6q6tmrAD4oqiYFUADmvh386Lit9o1MF2RVEaCXara392XlT8xlWtSIwxHqXUKo2bqyFO2CjokqwKJERp1K6tmmIk4eHhFBUV6U3kKscYQ1FREeHh4cEWRVEaLI3aIjHGvAe8N2DAgOlV3+vcuTOFhYVo0yslPDyczp07136goijV0qgVSU00a9aM6OjoYIuhKIrS4GnUri1FURQl8KgiURRFUXxCFYmiKIriE1dFiRQR+Q74Othy+IHrgH8FW4gQQtejMroeldH1qExd1uNGY0z72g66KhRJY0FEcj2pe3O1oOtRGV2Pyuh6VCaQ66GuLUVRFMUnVJEoiqIoPqGKpGGRHmwBQgxdj8roelRG16MyAVsPjZEoiqIoPqEWiaIoiuITqkgURVEUn1BFoiiKoviEKpIGjIjcJCLLRWRdTWNXGyLSW0TeEpE/iciEYMsTbEQkRUSWishrIvJRsOUJNiIyVESyHGsyNNjyBBsRiXWsxToR+VVd5lBFEiREZIWIfCsi+6qMjxKRf4rIIRGZVdMcxpgjxphptY01JPyxLsBdwCvGmF8B9wdM2HrAT9+TLGPMDOB94D8DKW+g8dP3wwDngHCgMFCy1gd++n4ccHw/fgbcVic5NGsrOIjIEKwv80pjTB/HWBjwJTAC6wueA9wHhAEvVpniIWPMt47z1hljKv3yrm6sIeCPdXE8zgN+AG41xtTpjyMU8PP35C1gmjHmbD2J73f89P34lzGmXER+BPx/Y8zk+pLf3/jr+yEi9wC/AlYZY97wVo6rth9JsDHGbBeRblWGk4FDxpgjACLyJjDWGPMiMKZ+JQwOflyXmY4/qLcDJWt94K/1EJGuQHFDViLg97+b00CLQMhZX/hrPYwxG4GNIvJfgNeKRF1boUUnoMDldaFjrFpEJEpElgKJIjLb3VgjwNt16SYi6cBKYH6AZQsGXq2Hg2nAnwMmUXDx9vuRKiL/AawClgRYtmDg7XoMFZHFjjXZVJcLqkXSgDHGFAEzahu72jDG5AMPB1uOUMIYMy/YMoQKxpi3aeCWqj8xxmwDtvkyh1okocUxoIvL686OsasdXZfK6HpUgWsxbwAAAyVJREFURtejMvW+HqpIQoscIEZEokWkOTAJ2BhkmUIBXZfK6HpURtejMvW+HqpIgoSIrAF2Aj1FpFBEphljSoFfAx8AB4C3jDFfBFPO+kbXpTK6HpXR9ahMqKyHpv8qiqIoPqEWiaIoiuITqkgURVEUn1BFoiiKoviEKhJFURTFJ1SRKIqiKD6hikRRFEXxCVUkilILIlImIrtF5HMR+UxEbq2Haz4kIntFZI+I7BORsY7x50TkjkBfX1G8QfeRKEotiMg5Y0yk4/mdwFPGmJ8E8Hqdgb8DScaYYhGJBNobY74K1DUVxRfUIlEU72iDVX7crpr6vv2GiCwRkamO5/ki8qLDkskVkSQR+UBEDovIDJfzt4vIfzmaEC0VkSZAB+AsVp8JjDHnbCUiIhkiMkFEBjjm3u2wXIzj/e4isllEPhWrC2Cvelwb5SpFq/8qSu20FJHdWB31bgCGeXjeUWNMgogsBDKwus+FA/uApY5jkoHewNfAZiAVeAf4BvhKRP4GvG2Mec91YmNMLpAAICLzHecCpAMzjDF5IjIIeNULeRWlTqgiUZTauWCMsW/atwArRaSPB+fZhfL2ApGOplJnReSiiFzjeC/bpQHRGuDHxph1IjIKGAgMBxaKSH9jzLNVLyAiE4EkYKTDBXYrsFZE7EMadOMmpWGgikRRvMAYs1NErgPaA6VUdg+HVzn8ouOx3OW5/dr+26sapDSO6xggG8gWkb9iNaV61vVAhzJ7FhhijClzuMXO2EpPUeoLjZEoihc4Yg5hQBGWO6q3iLRwWBjD6zBlsqPcdxNgIrBDRDqKSJLLMQmOa7nKcQ2wBrjfGPMdgDHmeyx32L2OY0RE4usgk6J4hVokilI7dowEQIAHjDFlQIGIvIUV8/gK2FWHuXOw2r32ALZixUe6AAtEpCNQAnzHlV0vxwI3AstsN5bDEpkM/ElE5gLNgDeBz+sgl6J4jKb/KkqQEJGhwOPGmDHBlkVRfEFdW4qiKIpPqEWiKIqi+IRaJIqiKIpPqCJRFEVRfEIViaIoiuITqkgURVEUn1BFoiiKoviEKhJFURTFJ/4HvOKbYrsD16sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def f(x):\n",
    "    return np.exp(x)\n",
    "\n",
    "def g(x):\n",
    "    return np.exp(x)\n",
    "\n",
    "def g2(x):\n",
    "    d = 1e-4\n",
    "    return (f(x+d) + f(x-d) - 2*f(x))/d/d\n",
    "\n",
    "x = 1.\n",
    "hs = 10.**np.arange(-12., -3., .05)\n",
    "\n",
    "fd = [(f(x+h)-f(x))/h for h in hs]\n",
    "fd2 = [.5*(f(x+h)-f(x-h))/h for h in hs]\n",
    "e1 = abs(fd - g(x))\n",
    "e2 = abs(fd2 - g(x))\n",
    "ho = sqrt(2.*f(x)*prec64/g2(x))\n",
    "\n",
    "loglog(hs, e1, '.-')\n",
    "loglog(hs, e2, '.-')\n",
    "axvline(x=ho, color='r')\n",
    "legend([\"Forward Difference\", \"Central Difference\", \"Optimal\"])\n",
    "xlabel(\"BumpSize\")\n",
    "ylabel(\"Error\");\n",
    "savefig(\"bumpsize.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Bump Size for Delta\n",
    "\n",
    "How to choose the $h$ so that the following finite difference approximation is the most accurate?\n",
    "\n",
    "$$ f'(x) \\approx \\frac{1}{h}(f(x+h) - f(x)) $$\n",
    "\n",
    "The best $h$ is when truncation error equals to the rounding error:\n",
    "\n",
    "$$ \\frac{1}{2} f''(x) h^2 = f(x) \\epsilon_m $$\n",
    "\n",
    "Thus the optimal $h^* = \\sqrt{\\frac{2 f(x)\\epsilon_m}{f''(x)}} $. \n",
    "\n",
    "It is however better off to simply use central difference:\n",
    "\n",
    "$$ f'(x) \\approx \\frac{1}{2h}\\left(f(x+h) - f(x-h)\\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## An Example of Numerical Delta\n",
    "\n",
    "For $f(x) = e^x$ at $x=1$:\n",
    "\n",
    "<center><img src = \"img/bumpsize.png\", height = 600, width = 600></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Condition number\n",
    "Consider the relative error of a function with multiple arguments: $f = f(x_1, ..., x_n)$:\n",
    "\n",
    "$$\n",
    "\\small\n",
    "df = \\sum_i \\frac{\\partial f}{\\partial x_i} dx_i \\iff\n",
    "\\frac{df}{f} = \\sum_i \\frac{x_i}{f}\\frac{\\partial f}{\\partial x_i} \\frac{dx_i}{x_i}\n",
    "$$\n",
    "\n",
    "Assuming input argument's relative errors are $\\epsilon_i$ (could be negative):\n",
    "\n",
    "$$\n",
    "\\small\n",
    "\\left| \\frac{\\Delta f}{f} \\right| = \\left| \\sum_i \\frac{x_i}{f}\\frac{\\partial f}{\\partial x_i} \\epsilon_i \\right|\n",
    "\\le \\sum_i \\left| \\frac{x_i}{f}\\frac{\\partial f}{\\partial x_i} \\right| \\left|\\epsilon_i\\right| \\le \\sum_i \\left| \\frac{x_i}{f}\\frac{\\partial f}{\\partial x_i} \\right| \\epsilon_m \\equiv k(f) \\epsilon_m\n",
    "$$\n",
    "\n",
    "where $\\epsilon_m$ is the maximum relative error of all inputs. \n",
    "\n",
    "* $k(f) = \\sum_i \\left| \\frac{x_i}{f}\\frac{\\partial f}{\\partial x_i} \\right|$ is defined as the condition number of a function, \n",
    "* it is the maximum growth factor of the relative error.\n",
    "* the calculation loses about $\\log_{10}(k(f))$ decimal digits of accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    " ## Well-posed and ill-posed problems\n",
    " \n",
    " Condition number is the systematic approach to detect potential numerical problems:\n",
    " \n",
    " * $f(x_1, x_2) = x_1 - x_2$ \n",
    " \n",
    "   $k(f) = \\frac{|x_1| + |x_2|}{|x_1 - x_2|}$, ill conditioned when $x_1 \\approx x_2$, catastrophic cancellation\n",
    "   \n",
    " * $f(x_1, x_2) = x_1 x_2$ \n",
    " \n",
    "   $k(f) = 2$, the multiplication is well conditioned\n",
    "  \n",
    "The problem is ill-posed if its condition number is large.\n",
    "\n",
    "* In practice, it is impossible to run condition number analysis over complicated calculations:\n",
    "* We have to rely upon monitoring and testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Well-posed problem can be unstable\n",
    "\n",
    "Consider the previous example of $f(x) = \\frac{1}{x^2}(1-\\cos(x))$:\n",
    "\n",
    "$$\n",
    "k(f) = \\left| 2 + \\frac{x \\sin{\\left (x \\right )}}{\\cos{\\left (x \\right )} - 1}  \\right|\n",
    "$$\n",
    "\n",
    " * $k(f) \\rightarrow \\infty$ near $x = 2\\pi$, it is ill-posed near $x = 2\\pi$ due to catastrophic cancellation.\n",
    " * $k(f) \\rightarrow 0$ near $x = 0$, it is well-posed near $x = 0$.\n",
    "\n",
    "The numerical algorithm can be unstable even if the problem itself is well-posed.\n",
    " * A better algorithm near $x=0$ is $f(x) = \\frac{2}{x^2}\\sin^2(\\frac{x}{2})$ \n",
    " \n",
    "But an ill-posed problem is always numerically unstable by definition.\n",
    "* Direct numerical calculation is unstable around $x = 2\\pi$\n",
    "* Does this help? $f(x) \\approx \\frac{1}{2x^2}(x-2\\pi)^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAE8AAAAYCAYAAAC7v6DJAAAABHNCSVQICAgIfAhkiAAAA71JREFUWIXt2FloXGUUB/Bfal1qlbqjgjBIrS244oJCfWhdAsUNfPPBBUWKqEhBpUWx+GBVqkFExJX2QVBpiVIVDUqrjVBrYsGg9cElFTVqtO5WbdPx4XxDJnfuvRMmaUpj/jAM92zfuec753znu0zhf4Ol+AC/YRDrcMpe9Wgfwpu4XgTsVHTiOxyxN53a06igilXjbPcQDOGyEpnV+AEzW1zjLOH7jWVCRyaBTnyGHfgV3bgB01pcnD0XvOOS3fkF/HOwG0vGuE4nBsRm5WJxcuRbPI8VeA6/JPoatLW4+P6YK152PPEStmC/An6X8H/GGNc5V8RgWZHAQpH+2Qw7Fl8l5avG6MR44hGx0ScW8OeIrHtqnNbbim1aqMBlIniPZeiX422R0v+Il3kHN2fkKhrLtp5WwQv4EX+jB5eW+NOR1pxbIvNAsn9hDq9LfjK0JX+qSb8e9yZ6e8maubgjKXbU0W5KtAGxu/fjGWwW40Q9KoqDt1409E3J/moRwCEsyPHlUXHCzmvicw92yT8oTk/2PzGy5B9OPj2Zo3NR4q1ssu4ITEefxqj3imw7JkfnqMxzRXHwqmJX69Ge6K9n6I+LGW+haCe1X7aRzxSB68t7oYRVaY3r0nOtul6UX5qzEn9zic0GrExKr2XovfgTh4/CRkVx8PrlN/xtoozrUS34Lc/IzUn0rhKfThATxZe4Jcm/gQNKdHaIrB8VbktGt2ocRJcYLtsOXImjC+xUFAfv5QKdblFareB8w1lUhhWGN+A9HNxE/huR0U1R242PRWnk4RrRq4aS7G7Rw87OyFWUHxh52JD4reCMpPtKE7laAlSVHz41bMfvzYRuTwb75Pe0LA7DIjwtAvmTkVlYMbHBOz7pdpfIXC02eyDJPtHE5rQk/3mZ0F3J2BaNjX80eFbjGFAxscFrEyf4YAF/Ef7FR2KTP8VOnFxic17yZ22NkD1V7hHzTa+Yj7INux4L5N84apn6V4nunkYV74rNn53hzRe3pa/FqT6Iu8VU8WCJzfPS//oaYXod81rcJ8puozgssug3nCmd+EP0vH4RyAvEfbIXb5U4MhFYK7K/XdzViV74qrizXyxKlghmD64Q77Axx94lIja5fXS54lGg9ttQJ79YBPALkWXbRanfiUMztismtmyJkeN7vJ+eZ4sx42ecliNfG4I35fBmiTGlaDKYlFgqAnLmGO3cqvzrzaTEQWLYXjcGGzPEnX1NllH0KWeyYBc+xIGip+1swcZJ4hr6kPi8NYUpTGHfxX/FDRnJBWA6zQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$2 \\sin^{2}{\\left (x \\right )}$$"
      ],
      "text/plain": [
       "     2   \n",
       "2⋅sin (x)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, k = sp.symbols('x, k', real = True)\n",
    "f = (1-sp.cos(x))/x**2\n",
    "k = sp.simplify(f.diff(x)*x/f)\n",
    "sp.limit(k, x, 0)\n",
    "sp.simplify(1-sp.cos(2*x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Ill-posed numerical algorithm\n",
    "\n",
    "Ill-posed problem are generally not suitable for numerical solutions.\n",
    "\n",
    "The Mandelbrot set is a famous example:\n",
    "\n",
    "* defined an iterative algorithm in complex number: $z_{i+1} = z_i^2 + c$\n",
    "* the color of $z_0$ is tied to the number of iterations for $|z_i|$ to go above a threshold\n",
    "* obviously unstable because of the iterative additions\n",
    "* showing exquisite details of numerical instability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><img src = \"img/mandelbrot.jpg\", height = 400, width = 400></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Rules of engagement for floats\n",
    "* Always allow errors when comparing numbers\n",
    "* Use double precision floating point number by default \n",
    " * unless you know what you are doing (really?)\n",
    "* Avoid adding and subtracting numbers of very different magnitudes\n",
    " * avoid using very small bump sizes when computing delta\n",
    "* Take care of the limits, especially around 0\n",
    " * many (careless) numerical implementations broke down around 0\n",
    "* Use well designed numerical libraries whenever possible\n",
    " * Even simple numerical algorithm can have spectacular pitfalls\n",
    "* Test and monitor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Portability\n",
    "\n",
    "* Floating point computation is usually not reproducible across different platforms\n",
    "* Same code will produce different results with different compiler and hardware\n",
    " * The difference can be significant in optimization and iterative methods\n",
    "* Design your unit test to be cross platform by giving enough error allowance\n",
    "* Test and monitoring are critical"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Assignments\n",
    "\n",
    "Required Reading:\n",
    "\n",
    "* Bindel and Goodman: Chapter 2\n",
    "\n",
    "Recommended Reading:\n",
    "\n",
    "* Review linear algebra (any elementary linear algebra book)\n",
    "* Andersen and Piterbarg: Chapter 1\n",
    "\n",
    "Homework\n",
    "\n",
    "* Setup your python environment following the instructions at: http://yadongli.github.io/nyumath2048/\n",
    "* Bindel and Goodman: Problem: 2.8\n",
    "* Complete [Homework Set 1](http://yadongli.github.io/nyumath2048/)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}