{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Scientific Computing in Finance\n",
    "\n",
    "* MATH-GA 2048, Spring 2020\n",
    "* Courant Institute of Mathematical Sciences, New York University"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Instructors\n",
    "\n",
    "#### Wujiang Lou\n",
    "* Senior trader, Global Fixed Income, HSBC\n",
    "* Ph.D. in Aeronautics and Astronautics, Nanjing University\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 | H. Cheng |\n",
    "| 2-3 | Linear Algebra | Portfolio optimization, PCA, least square | H. Cheng |\n",
    "| 4-5 | Rootfinding and interpolation | Curve building | H. Cheng |\n",
    "| 6 | Derivatives and integration | Hedging, Risk Transformation | H. Cheng |\n",
    "| 7-8 | Monte Carlo simulation | Exotic pricing | H. Cheng |\n",
    "| 9-10 | Optimization | Model calibration | W. Lou |\n",
    "| 11-13 | ODE/PDE | Pricing | W. Lou |"
   ]
  },
  {
   "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 only ( version 3.7.x is not supported currently) ** : 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",
    "* George Box: \"essentially, all models are wrong, but some are useful\"\n",
    "\n",
    "Truncation Errors: introduced when exact mathematical formulae are represented by approximations.\n",
    "* Convergence error\n",
    " * Monte Carlo simulation\n",
    " * Numerical integration\n",
    "* Truncation and discretization\n",
    "* Termination of iterations\n",
    "\n",
    "Roundoff Errors: occur because computers have a limited ability to represent numbers\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": 6,
   "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": 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>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": 8,
   "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": 9,
   "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": 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>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": 11,
   "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": 12,
   "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": 31,
   "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>00000000000000000000001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decimal</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>-127.0</td>\n",
       "      <td>2.384185791015625e-07</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table></font></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f2parts(-1e-45, 32);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Floating point computing is NOT exact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 15,
   "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": 16,
   "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": 17,
   "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": 18,
   "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": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n"
     ]
    }
   ],
   "source": [
    "print(1. - 1e-16 == 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "print(1. + 1e-16 == 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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": 22,
   "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": 23,
   "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": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEJCAYAAABv6GdPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8XNWZ8PHfmaJuVUsuki3LvTdkA6FjL5jgmBACi5MsC4SSEDYkL3U3vMDCsptsspAl4V1CICG7ycaBgINxqDamGtzAvcq2sGTZalax+pTz/jG6Y5Upd2auNCrP9/Px5+MZ3Zk5lu8888w5z3mu0lojhBBiaLHFewBCCCGsJ8FdCCGGIAnuQggxBElwF0KIIUiCuxBCDEES3IUQYgiS4C6EEEOQBHchhBiCJLgLIcQQ5IjXC48cOVJPmDAhXi8vhBCD0rZt22q01rnhjotbcJ8wYQJbt26N18sLIcSgpJT6wsxxMi0jhBBDkAR3IYQYgiS4CyHEEBS3OXchRP9zuVyUl5fT1tYW76GIMJKSkigoKMDpdEb1eAnuQgwj5eXljBgxggkTJqCUivdwRBBaa2praykvL6eoqCiq55BpGSGGkba2NnJyciSwD3BKKXJycmL6hiXBXYhhRgL74BDr/5MEd9GvdpU38PmxungPY0hZt7eSivrWeA9DDDAS3EW/enTtHv719X3xHsaQcucfP+PXHx6J9zBMU0px9913+2//7Gc/45FHHunXMWzdupXvf//7UT324osvDrgB88MPP2TWrFnMnz+f1tZWTpw4wfLly0M+19q1a3n44YejGkc4EtxFvyo71Uq72xvvYQwZWmvaXF7KTg2ezD0xMZFXXnmFmpqauLy+2+2muLiYp556ytLn/cMf/sA999zD9u3bSU5O5oknnuDWW28N+Zgrr7ySNWvW0NLSYulYQIK76Ecdbi+Vp9tweXS8hzJkeLy+3+XxQTQt43A4uO2223jyySd7/ezGG2/kz3/+s/92WloaAO+99x4XXXQR1113HVOnTuWBBx7gD3/4A4sXL2bOnDkcPnwYgOrqaq655hoWLVrEokWL+PjjjwF45JFHuO2227jsssu44YYbeO+99/xZdVNTEzfddBNz5sxh7ty5vPzyywB897vfpbi4mFmzZoXNrp977jlefPFFHn30Ub75zW8C8PLLL7Ns2TIAnnjiCW6++WYAdu3axezZs2lpaUEpxcUXX8zatWuj/n0GI6WQot+cbGhDa3B7JHO3itsI7nWRZ37//Noe9lY0WjqemWPTefgrs8Ie973vfY+5c+dy3333mX7uHTt2sG/fPrKzs5k4cSK33HILmzdv5j//8z/5xS9+wc9//nPuuusufvjDH3L++edz7NgxLr/8cvbt800Dbtu2jY8++ojk5GTee+89//M+9thjZGRksGvXLgDq6nxrQo8//jjZ2dl4PB6WLFnCzp07mTt3bsCx3XLLLXz00UcsX76cr3/96xw9epSsrCwSExMB+MEPfsDFF1/M6tWrefzxx/nVr35FSkoKAMXFxXz44Ydcd911pn8XZkhwF/2mvDMAGQFJxM7V+UHZ2Oamsc1FelJ0G176W3p6OjfccANPPfUUycnJph6zaNEixowZA8CkSZO47LLLAJgzZw4bNmwAYN26dezdu9f/mMbGRk6fPg3AihUrAr7WunXrWLVqlf92VlYWAC+++CLPPvssbrebEydOsHfv3qDBvacTJ06Qm3umcaPNZuOFF15g7ty53H777Zx33nn+n+Xl5VFRUWHqeSMhwV30m/LOqQOXZO6WcXeZ4jpe10r6GPPB3UyG3Zd+8IMfsHDhQm666Sb/fQ6HA6/Xd35oreno6PD/zMiCwRcsjds2mw232w2A1+vlk08+CRjEU1NTA45Da92r7PDo0aP87Gc/Y8uWLWRlZXHjjTdGVHOenJzc6/hDhw6RlpbWK5C3tbWZ/oCLhMy5i35zvM4X3N0y524Zl/fMB6Xx+x0ssrOzue6663j++ef9902YMIFt27YB8Oqrr+JyuSJ6zssuu4xf/vKX/tvbt2+P+DF1dXU0NjaSmppKRkYGlZWVvPHGGxGNY+rUqZSWlvpvNzQ0cNddd/HBBx9QW1vbbV3h4MGDzJ49O6LnN0OCu+g3xqKf2yuZu1W6Ze6DaFHVcPfdd3ermrn11lt5//33Wbx4MZs2bQqabQfz1FNPsXXrVubOncvMmTN55plnwj7mwQcfpK6ujtmzZzNv3jw2bNjAvHnzWLBgAbNmzeLmm2/uNo1iRmpqKpMmTaKkpASAH/7wh9xxxx1MnTqV559/ngceeICqqioANmzYwJVXXhnR85uhtI5PFlVcXKzlYh3Dy/XPfsKnR06Rkexkx8OXxXs4Q8Kx2hYu/Klvvvm2CyfyT1+eEfL4ffv2MWNG6GOENVavXs22bdv4l3/5l6DHVFZW8o1vfIP169cH/Hmg/y+l1DatdXG415c5d9Fv/Jm7zLlbpuu0THkUFTOi71x99dXU1taGPObYsWP8x3/8R5+8vgR30S88Xs2Jet8Ck0uqZSzTc0FVDCy33HJLyJ8vWrSoz15b5txFv6g63Ybbq8lJTZDM3UJG5VFOasKgnHMXfUeCu+gX5Z1ZZWFOCl4NXsneLWHsGSjMSaGmqYM2lyfOIxIDhQR30S+MKYMJI33VDy6pmLGE8S3I+L1K9i4MEtxFvzCCTmG2LwhJrbs1jD49E3J8v9dymXcXnSS4i35RXtdKTmoCI5J8a/gS3K1h7BkozPH1KRnoi6q1tbXMnz+f+fPnM3r0aPLz8/23u+5GjcXp06fJycmhqamp2/3Lly/nlVdeCfq4devW8dWvftWSMQwEUi0j+kV5XQv5Wck47b5t3jItYw3jQ7IgKxmHTXG8fmCXQ+bk5Ph3jT7yyCOkpaVxzz33dDtGa43WGpstutxzxIgRXHrppbz66qv+Do11dXVs2rSp287QoU4yd9Evjte3kp+ZjMPuO+Ukc7eGUS2T6LAzOiNpwGfuwZSUlDB79my+853vsHDhQsrKysjMzPT/fNWqVf6ywsrKSr72ta9RXFzM4sWL+fTTT3s938qVK7s1A3v55Ze58sorSUpK4tNPP+Xcc89lwYIFnHfeeRw6dKjX4x988EF+/vOf+29Pnz6d8vJyAH73u9+xePFi5s+fzx133OHvhTPQSOYu+pzWmor6Vi6dlofD1pm5SzmkJYxqGYddkZ+ZHNmc+xsPwMld1g5o9By44sdRPXTv3r389re/5ZlnnvE3Agvk+9//Pvfddx/nnHMOpaWlLF++nN27d3c75sorr+T222+nrq6OrKwsVq1axb333gvAjBkz+Oijj7Db7bz55ps8+OCD/OlPfzI1xt27d7N69Wo2btzo70u/atUqvvGNb0T1b+5LEtxFn6tt7qDN5aUgKxmnkblLKaQljA9Jh81GflYynxwOvSNyIJs0aZKpTT3r1q3jwIED/tt1dXW0trZ266yYmJjIlVdeySuvvMLy5cvZs2cPS5YsAaC+vp4bbrjBf4GPSKxbt44tW7ZQXOzb/d/a2sq4ceMifp7+IMFd9Dkjm8zPSqHd7avDlo1M1jCmt5x2RUFmMpWNbbg8Xv+HaEhRZth9pWuTMJvNRte+V13b52qt2bx5MwkJCSGfb+XKlfzsZz+jtbWVr33tazgcvnD3ox/9iMsvv5w77riDkpIS/9WSuuraerjr62utufnmm3nsscei+0f2I5lzF33OmAfOz0zG0blIJpfas4ZRLeOw2yjI8m0QO9lgvu/4QGWz2cjKyuLQoUN4vV5Wr17t/9nSpUt5+umn/beDtfVdunQpe/bs4ZlnnmHlypX++xsaGsjPzwfghRdeCPjYrq2HN2/eTFlZmf85X3zxRX8ny9raWo4dOxb9P7QPmQruSqllSqkDSqkSpdQDAX4+Xim1QSn1uVJqp1Lqy9YPVQxWRgVH12oZaftrDeND0mlT5Gf5piWGSq37T37yE5YtW8aSJUsoKCjw3//000/z8ccf+9v6/vrXvw74eLvdztVXX01jY2O3lr33338/9957b8g2vtdeey2VlZUsWLCA559/nokTJwK+qz49/PDDLF26lLlz53LZZZdRWVlp0b/YYkbZUbA/gB04DEwEEoAdwMwexzwLfLfz7zOB0nDPe9ZZZ2kxPDz0l1169sNvaq21fu9AlS68f63eWnoqzqMaGn770RFdeP9aXdvUro9WN+nC+9fqF7ccC3r83r17+3F0IlaB/r+ArTpMfNVam8rcFwMlWusjWusOYBVwVc/PCCC98+8ZgPUXBBSDVnmdrwwSfBkmyJy7VbpWy4zJTAKkBYHwMbOgmg+UdbldDpzd45hHgLeVUv8ApAJLLRmdGBKO17dS0Dll4JBqGUudmZaxkeiwMyo9cdDWugtrmcncVYD7er4zVwIvaK0LgC8D/6OU6vXcSqnblFJblVJbq6urIx+tGJSO17VSkOXbHu+wS527lYxvQMbvNT8zOWzmruN09TURmVj/n8wE93KgayFnAb2nXb4NvNg5oE+AJGBkzyfSWj+rtS7WWhfn5uZGN2IxqDS0ujjd7u4yLSM7VK1kXPjE2ByWn5USMrgnJSVRW1srAX6A01pTW1tLUlJS1M9hZlpmCzBFKVUEHAeuB3puxzoGLAFeUErNwBfcJTUX/ku/5funZaRaxkpujxeHTaHUmcz9zd0n8Ho1NlvvL90FBQWUl5cj35wHvqSkpG5VQpEKG9y11m6l1J3AW/gqZ36jtd6jlHoU36rtGuBu4NdKqR/im7K5UUtqIDhT427Mufsbh0nmbgm3V/s/MMH3e3Z5NFWn2xmd0TvrczqdFBUV9ecQRZyY2qGqtX4deL3HfQ91+fteIHjRqBi2jCkCY1rG2MQkmbs1XB6vf6oLznxDOl7fEjC4i+FDdqiKPnW8rpUkp43sVN9WcYdk7pZye3pk7plDayOTiJ4Ed9GnjtY0MyEn1T8n7JSWv5Zye73+8lKAcdkpKOX7vYvhTYK76FOHq5uYlJfmv21Udci0jDVcHu3fGAaQ5LQzLiuFw9US3Ic7Ce6iz7S7PRw71cKk3C7B3S6Nw6zk9nTP3AEm5aZSUtUU5BFiuJDgLvpMaU0LXu0LNgZ/4zDZxGQJV49qGYBJuWkcqW7CK7uAhzUJ7qLPHK72ZY/dMnebtB+wkrtHtQzApLw02t1e6TEzzElwF33mcOfUwMQAmbu0H7BGz2oZOPNhany4iuFJgrvoM4erm8jPTCYl4cx2CqUUdpuSahmL+KZlur+NJ+cZwV0WVYczCe6izxyubu5WKWNw2BQuqZaxhG9apnvmnp2aQFaKUzL3YU6Cu+gTWmtfGWSXKRmD026TzN0igaZlwDc1IxUzw5sEd9EnTja20dLh6baYanDYlVTLWMTlDXwxbKNiRgxfEtxFnzhc5ZvvDRjcbTZ/q1oRG7dH+zeGdTUpL5Wapg7qWzriMCoxEEhwF32ipOo04AsyPTklc7eMK8AmJuhaMSOLqsOVBHfRJw5XNzMiyUFuWmKvn/mmZSRzt4Lbq/3lpV1JOaSQ4C76hG8xNc3fMKwrp0zLWMZ3sY7eb+OCrGQS7DYJ7sOYBHfRJ4zgHogsqFrHFaRaxmG3MWFkin8jmRh+JLgLy51uc1HZ2O7fTNOTw2aTxmEWcXt7tx8wTM5Lkzn3YUyCu7DckWqjUqb3Yip0LqjKJiZLBKtzB9+8+7FTLbS7Pf08KjEQSHAXlvM3DAuWucsmJsu4PIHr3MEX3D1ezbHaln4elRgIJLgLy5VUNeGwKcZnpwT8ucOmpHGYRdzewHXuIBUzw50Ed2G5w9VNFOakBM0onXabtPy1iG9aJvDv2ejGKfPuw5MEd2G5w9XNQStlQKplrORrPxA4c09NdDAmI0l6zAxTEtyFpVweL1/UBu4GaZBqGWt4vBqtCVjnbpiUmybTMsOUBHdhqbJTLbg8OmTmLtUy1jDWLYJVy4CvYulwVRNay4fpcCPBXVjKmN8NVuMOUi1jFWPdIti0DPj+H5o7PFQ2tvfXsMQAIcFdWOpgpa9h2MQgNe4ATrlYhyWMdYuQ0zKdH7IHOv9fxPAhwV1YavfxBgpzUkhPcgY9RhqHWcNYtwiVuc8akwH4/l/E8CLBXVhqd0UDs8dmhDzGYZcFVSsY6xbBSiEBMlKcjMtOZk+FBPfhRoK7sEx9Swdlp1qZnR86uDttsqBqBePbT7BNTIY5+Rnsksx92JHgLiyzp6IR8AWTUOw2WVC1gplqGYDZ+RmUnWqlocXVH8MSA4QEd2EZIzucNTY95HFOu7QfsIJRLRNqQRXwT5PtlqmZYUWCu7DMruMNFGQlk5WaEPI4h11J+wELGB+QoRZUAf80mSyqDi8S3IVl9hwPv5gKvkzTt7tSAnwszsy5h34bZ6cmkJ+ZLPPuw4wE92GoodXF/pONlj5nY5uL0toW5hSED+5GpikVM7E5Uy0TOnMHmJ2f7l8TsdK2L+rwyLewAclUcFdKLVNKHVBKlSilHghyzHVKqb1KqT1Kqf+1dpjCSj9+Yx/XP/uppc+557gvcISrlIEzpXtSMRObM3Xu4d/Gc/IzOFrTTGObdYuqZadauOa/NrJqyzHLnlNYJ+xZoZSyA08DVwAzgZVKqZk9jpkC/CNwntZ6FvCDPhirsIDXq3lnbxX1LS463NYFV2M+d3aYxVQ4U7onmXtszJZCAszq/NA1PoStUHXa19Lgnb2Vlj2nsI6ZzH0xUKK1PqK17gBWAVf1OOZW4GmtdR2A1rrK2mEKq+w83kBNk+9N2dBqXRa363gDYzOSyElLDHuskWlK29/YuExsYjIYayFWbmZq7Dx/Nh6upaXDbdnzCmuYCe75QFmX2+Wd93U1FZiqlPpYKfWpUmqZVQMU1np335ksy8rgvruiwZ8dhmPMEUvFTGzcJtoPGHJHJDI6PcnSRVXj/Olwe/m4pNay5xXWMBPcA505Pd+VDmAKcDGwEnhOKZXZ64mUuk0ptVUptbW6ujrSsQoLrNtX5Q8GVgX3pnY3R2uaw25eMjg7qzuk1j02ZhqHdTU7P8PSckjj/HHaFev3ydTMQGPmrCgHxnW5XQBUBDjmVa21S2t9FDiAL9h3o7V+VmtdrLUuzs3NjXbMIkonGlrZe6KRS6fnAVi2uLbneANah9+ZavBn7jLnHhOXiZa/Xc3Jz+BITTNN7dZMoRjTMpdMy+Pd/VV45ZvYgGImuG8BpiilipRSCcD1wJoex/wFuARAKTUS3zTNESsHKmL37n7fUsjVCwqAM2/OWO3uLLGblR9+MRWkWsYq/szdxJw7+MohtYZ9J6xZVG1odZGSYGfZ7NFUnW7vk1JLEb2wZ4XW2g3cCbwF7ANe1FrvUUo9qpRa0XnYW0CtUmovsAG4V2stk3AW2neikS9qY7vQ8fp9VYzLTuaswizAummZ3ccbGJWeSN6IJFPHO6VaxhKRVMvAmW9Wu8qtmZppaHWRnuTk4ml5KAXrYpya2Vlez4mGVkvGJkzWuWutX9daT9VaT9JaP95530Na6zWdf9da6/+jtZ6ptZ6jtV7Vl4MebjrcXr713CYeW7s36udo7fDwcUkNS6aPIiPZ12vdqkZSu03uTDX4M3cJ7jExqmXM1LkD5KUnkTsi0bIeMw2tLjKSnWSnJrBwfJb/m2G0bvvvbdzz0g5LxiZkh+qg8M7eSmqbfe10o7XxcA3tbi9LZuSR4LCR7LRbkrm3dLg5XN1kavOSwZhzl6sxxcafuZuccwdf9m7VoqoR3AGWzMhj1/EGKhvbonqudreHk41tfFxSG/M3VOEjwb2flNe1cN0zn1Be1xLxY40dgBUxfGVdt6+K1AQ7i4uyAchIdlqyoLq3ohFvBIupcKZaRjL32Pgbh5mslgFfxUxJVZMldemNbW7SjeA+fRRA1Nn7yYYzHwp/2lIW4sjAPjxUzU2/3Uy72xPV6w9FEtz7yX+9d5jNpafYGGE9cNmpFj48VEN2agKn29xRVTporXl3fyUXTMkl0WEHfMHdisx9S2kdAPPG9ap8DepMtYxk7rHwt/yNIHNfMD4Tr4bPj9XH/PqNXTL3qaPSyM9MZv2+6IJ7Rb0vuGenJvDStvKIy2Tf2nOSDQeqWf3Z8ahefyiS4N4Pqk638dK2cgAOVUV2oeI/bSnDpuDWCyYCcDKK7H1PRSOVje0smZHnvy892WFJcN98tJZJuankjgi/M9XgbxwmpXMxcZu8WEdXxYVZ2BRsOnoq5tdvaHWRnuwAQCnF0hl5fFxSQ5sr8uz5ZKPvvP7ORROpPt0e8TeAQ5VNAPzqgyPSyKyTBPd+8MLHpbg8XnJHJHKoqsn049weLy9tK+OiqbksHO/LjI0MJxJv7TmJTcEl088Ed1/mHttXc49Xs7W0jrMn5kT0OIdN2g9Ywd84LIJpmRFJTmaNzWDTkdiK2dweL03tbn/mDrBkxihaXR4+OBj5BkXjvF65eDyj0hNZtTmyZmQlVU2MTk/iaE0zb+85GfHrD0US3PvY6TYX//PpF1wxezTnTszxZxhmbDhQTWVjO9cvHs/YzGSAiEvFtNas2VHBuZNyGNml70t6sjPmOve9FY2cbndzduc8vlkOaflrCbfXi02BzWQppOHsomw+L6uPKsM2NLb5EoOuwf3cSTlkpyawZkfPPY7hnWhoJTPFyYgkJ9eeNY73D1ZTUW/uXD/V3EFtcwc3nTeBwpwUnnn/sFwrgGEe3PefbORL/7bedN3vxpIaXtoa2WLPHzcf43Sbm+9cNIkpeWkcr2+l2eS8+arNx8gdkcil0/MYlZ6EUpFn7ruON/BFbQsr5o3tdn+GBcF901Ff9nd2UWSZu1M2MVnC7dGmNzB1dfbEHDrcXnbGUO9unDtdg7vTbuPLc0azfl9VxAu2J+rbGJPhS2D+dtE4vBpeNPleK+n8Njxt9Ahuu3AiO8ob+CSCbyZaa36x/hBHa8xV6by0tYzLnnzf0vbJfWHIBfenN5Tw/EdHwx6nteb//mU3FQ1trN0ZPtNoanfz/VXbeejVPaanE9rdHp778CjnTc5hbkEmU0alAXCkOvxJdKKhlQ0Hqrj2rAKcdhsJDhsj0xIjztzXbK/AaVcsmzWm2/0ZyU5Ot7tjmp/cdPQUhTkpjM4wt3nJYGy6kWqZ2Lg82r8hLBKLJmShFDFNzTQECO4AX5k7llaXJ+I2wBUNbYztPI/GZadwwZSRvLilzNT5aaxjTRk1gmsWFjAyLZH/eu+w6dc+VNXEf7xzkAde3mkq43/ls+McrGziyXcOhj1Wa82Df9nFm7tPmB6PVYZUcN9YUsNP3zrAk+8cDPuVc/Xnx9lSWkdaooP3TcwRPvvBEWqa2ml1edh/0tyi6F8+P07V6Xa+c9EkACbn+YK7mUXVl7aW49W+LMYwNiOJEw3mM3evV7N25wkumppLRkr3N2F6ku92tNm716vZUnoq4ikZOJO5S+Ow2Li93qgy98yUBKaNGhHToqoR3NN7BPdFE7IZk5HEaxFOzZxoaGVM5pkk4fpF46loaOODQ+Hfm4cqm0hJsDM2I4kkp52bzpvAh4dqTNfzf/aFr+Jr09FTYat9mtrdbP3iFGmJDn63sTRsK4f9J0/z+0+Pce9LO/t99+2QCe4tHW7uf2UnIxIdNLW7ee9A8P+kxjYX//r6fuaPy+R7l0xm/8nT3epse6psbOPXHxyhuHPb/ufH6sKOx+PV/OqDI8wam875k0cCUJiTisOmwi6qdri9/GHTF5w/eSSFOan++8dkJEcU3LeUnuJkYxtf6TElA2cyrmgrZg5Unqa+xcXiCKdkQFr+WsXl0aabhvV0zsQctn1RF/UHbLDM3WZTLJ87hvcPVlPf0mHquVo7PNS3uPzTMgB/M3MUI9MS+d3G0rCPL6lqYnJeGkr5fhffOqeQtEQHz7xvLnv/7FgdmSlOJuam8m9v7Av5zfyTw7W4PJqffn0uGclOHnp1d8hs/7UdFdhtCpfXyz+9sqtf1wKGTHD/6VsHKDvVyq9uOIuRaaEXdZ585yC1ze08dtVsLpnu604ZaoX/ibcP4vZ6eeK6+YxMS+QzEzXCa3dWcKS6me9ePMl/0jntNopGpoZdVP3rrgoqG9v59gVF3e4fk5nEifpW0yfImh0VJDvt/M3MUb1+Fmtw39yZ9UWTuUu1jDXcHq/pdr89nV2UTavLE3V/92DBHWDFvHxcHs2bu81VrRib88Z2ydwTHDb+/txC3jtQzaHK0N90D1Wd9n8rNsb0rXMK+euuE2EfC/DZsXoWjs/igWXTOVzdzKoQm6jeP1hFSoKdJTNGcf+y6WwprWP154Fr67XWvLazgvMmj+Tey6ez4UA1f9nef3X4QyK4by09xQsbS7nh3EK+NGkkV84Zw/p9VZwOsOCxt6KR320s5RuLxzOnIINpo0YwKj0x6NTM/pONvLStjL8/dwLjc1I4qzCTz8Jk7m6PlyffOcj00SP48uzuc91TRqVREmJaRmvNcx8eZXJeGhdN6d4WeUxGEs0dHn+lQiguj5fXd51gyYw8UhIcvX5uTNNEuyi06Wgt+ZnJjMtOifixcoFsa7i9OqIa964WdX4obzoS3dSMcd4ECu6z89OZkJPCaybWssC3mAowOj252/3fPKeQRIeN33wcfA2tsc1FZWM7U/JGdLv/9gsnkprg4Ml1oefFG1pclFQ1cVZhFn8zcxSLi7L5+bqDAWOH1pr3DlTzpUkjSXDYuK54HPPGZfKvr+8P+D76vKyeslOtrJg3lhu/NIGzCrN4ZM1eqk5H16IhUoM+uLe5PNz3552MzUjmvmXTAVgxfyztbm+vLnVer+bhNbvJSHZy7+XTAN/mi4um5vLhoeqAmeS/vb6ftEQHd146GYCF47P4orbFf6m6QF7+rJzS2hbuvmxarzK1yXkjOHaqJeiawKdHTrGnopFvn1/U67HG11Yzc3cfl9RQ1+LqVSVjiCVz11qz+egpfyuDSEnLX2u4PF7TTcN6GpmWyOS8NH/FU6QaWl0kOGwkOe29fqaUYsW8sXxyuNZUIAuUuYNvt+o1ZxXw8mfHqQ3yfjMqZaZ0ydwBslITuPn8Il7fdTLk3PtnZb5EbcH4TJRS/OgX+s56AAAgAElEQVTLM6hp6uBX7/fuWH60ppnyulYumuZLumw2xb9cNZva5naeeLv3h8hrOypIcNi4bNYo7DbFT66ZS6vLw8Ov7gk6HisN+uD+xDsHOVLTzI+vmUNaoi9DXTAui/zMZNZs7545/HHLMbaU1vHAFdPJTEnw33/R1Dwa29zsKO8+3fLhoWreP1jN95dM8R+/sHPe3ViE6and7eGp9SXMG5fJ0i47Qg2T89LwaoKWXT3/0RGyUxO4ekHPKxmeOflPmCiHXLOjghFJDv+J2JOxoBpNcD9c3UxNU0dUUzIgF8i2itujTbf7DeTsomy2ltZFVTHV2NnuN5ivzBuLV8Nfd4avEjHWuwJVXd18XhEdbi+//zTwpqaSzinOyT2CO8AtFxSRkewMWdXy+Rd12BTMK/BtEpw3LpMV88by3EdHeiVRxrf7rt+o5xRk8K2zC/nvT0rZXnYmfng6ixkumZbr/z1Nzkvjh0un8sbuk7y+q++rZwZ1cN94uIZff3iElYvHc0GXX7jNplg+bwwfHqqhrtm3qFNR38q/vb6fL03K4bricd2e5/zJI7EpeO/AmamZDreXR9bsYXx2Cn93bqH//jn5GThsKui8+5+2lHG8vpV7Lpvqn2vvaoq/Yqb3vPuR6ibW7aviW+cUBsyIjMw9XAOxNpeHt/dUsmzWaH8vmZ5iydz99e0R7kw1OKXlryWirZYxnD0xh6Z2N3ujuMiGryNk7+k+w5RRI5g+eoSpDU0nGloZmZYQ8FydnJfGpdPz+J9PSwN+2z1UdZoEhy3g9GB6kpPbL5rI+v1VQadSPztWz/TR6aQmnvm33Hv5NLSGf/nrvm7Hvn+wmokjUxmf0/217ls2jVHpSdz35x3+xmWbjtRSfbqdFfO6J2m3XlDENQsLGB/FdGakVLx2chUXF+utW7dG/sA3HoCTu3B7fZswbEoxpyADe49A2tzhZtfxBopGppI3IpEDJ0/T2OZmbkEGSQFOot0V3S8Vd7y+lbK6FqaNGkFWlywfYFdFAzYFs8Z074To0ZrtZfUkOW3MHJOOCnD5Wa/WbC495Zuvzur+H3y0pomq0+0sGJ9FQoA3rUaz6Wjgx3ZV09xOSVUT00ePIDM5IeAxGt/UyuiMJAqzUwMeE8yhKt/vcuH4zID/RjM+NebsQ/w7RGj7Tjbi9uiIOnJ21eHx8tmxOgqzU7pVqpix90QjXq1D9vE33kPzCjJJDpCsGML9OxpaXew72cjEkam9Lgiz/2QjHR4vc/MDN67zaM32Y3UkJziYOab7lcI0vvYZI9MSKRrZ/T1QXtdCeX0rM0ank5HsxKs1W784Rd6IJCbk9H6/1LV0cKDytP+cPlLTRE1TB2cVZvWKTQCMngNX/Djo7yQUpdQ2rXVxuOMGZeau0Rypacbl8TI5Ly3gLy8lwU6S005tUzs1TR3Ut7oYl5UcMLADZKY4ae5w4/J4aXd7OF7fSlZKQq/ADjAi0UFzuxtvj+uEVza24fJ4GZeVEjTo2ZQiyWGntaN7FuL2eqluamdkWmLAwA6gUCTYbbS7Q89VVzW2k2C3BVzs6vpcDpst4q/kGu1r9ZrkiDqw+14fZId4bLSGQHHDrAS7jSSHzdQCfU8eb/gpIaOZXLh59w63lwRH8FCUnuwgxWnnREMbusd7rtXlCfnBYVeKsZnJNLa5aOix6NnS4cGjNWlJvb+BjM1MJtFh42hNM16taWxz4dX02i9iyEpJYGRaIhX1rTS1u6lt7iA7xRk4sPcXrXVc/px11lk6Wi9tLdOF96/VT284FPK4J985oCc8sFbPefhNffXTH2m3xxv02O3H6nTh/Wv16s/K9e3/vVVPe/B1XXaqOeCxa7Yf14X3r9U7yur899U3d+j5//yW/tZzn4Yd/y2/26KX/Md73e57at1BXXj/Wr3vREPIx1799Ed65bOfBP35keomXXj/Wv2L9QfDjuPSn23Qd/x+W9jjujpcdVoX3r9W//cnpRE9rqcZ//cN/dhre2J6juHu2mc26muf2RjTc9z30g495+E3tcvtiehxF/zkXX3XHz8Le9xt/71FL3j0bd3mcgc9ZvZDb+qHX90d8nmM9/yG/ZX++5rbXbrw/rX6qXWhz/XWDrc++/F1+qpffqS93jMx4PeflurC+9fq0pqmgI/bsL9SF96/Vv/y3UP60df26Ck/el23tAf/d5xqatdnPfa2nv/Pb+nC+9fq9ftOhhxXtICt2kSMHXSZ+xe1zTz86m4WF2Vz+4WTQh67Yt5YtIY2l5d///pc7CEyjTn5GWSnJvDU+kO8ueck/3DpFAqCTBkEWlR96t1D1Le6eOCK6WH/DZPz0ijt/OYBvl1vz398lEum5TJ9dOiLTI/JDL2RadWWY9htqte6QiDpUfR0N1qxXjw18EKtWQ6bkk1MMXJ7vFFvYjJcNC2Xxja3qb0bXfna/Qb/Zmj4xtmFnGruCNqO4HSbi9PtbsaEaWGxYt5Y8jOTeWr9If8+j8NVvqKEQIupXSU57dx92VS2l9V3WwP47It6clITgs5/Xzwtj2WzRvOLdw/x+q4TnF2UTXJC8G8JWakJPHrVbOpafH3uz58c23skVoMuuK/deQKbTfHk384PGawBJuamcf2icTy8YiaTe9TB9mSzKS6YMpIjNc1MzE31908PZGxGEqPTk/xviCPVTfxuYyl/WzyOWSauJTolLw23V/svJ/a7jaXUt7j4wdKpYR87NiOJiiAbmTrcXv68tZylM/LISw/f7yWaC3as31fF1FFpUdW3d+W026T9QIzcXh31JibDBVNG4rQr1kdwcWuv1zdNEWraz//8k0eSn5nM/24KXO1iJCpjMkPP+Sc4bNxxySQ+O1bPh4dqACipNnrKhA7uANcsLGB2fjo/fmO/f0r082N1LCzMClj4YHjoKzNRKE40tHGRiYTmy3PGcPN5Rdx5yeSQU039YdAF9+9dMpm3f3gh+WFOBsOPr5nLN88uDH8gvi3PSsGjK2aH/I9RSrGwMJNtnZn7v76+rzM7mGbqdYwNF4cqm2hqd/PrD49wybRcU1czGp2RTLvbS12Ai1u/vfcktc0drFw83tQ4Ig3uDa0utpSeYsmM3jteI+WwK6mWiVEs7QcMI5KcnF2Uw/oILo5xut2N1oE3MPVksylWLh7HxsO1Act/jba+4TJ3gGvPGsfYjCR+vu4gWmsOVTbhsKluLTpCjeOh5bM40dDGrz44zKnmDo7UNLNwfFbIx43NTOYHS6f0uh5CKA99ZSa3Xhg8Oewvgy64AxGv7Jt15ZwxfHz/pZw/ZWTYYxeOz+J4fSsvbytn3b4q7rx0sumrEU3K852Mh6qa/Fn7XSaydsDfOS9Qr+s/bj5GfmYyF04x93Uw0uuofnCwGrdXs8TkSR6Kw2aTC2THKJb2A10tmZFHSVWT6QtTNwZpGhbMdcXjsNuU/1rAXfkzdxPB3Ze9T+azY/V8VFLDoaomikammt7ItbgomyvnjOGZ9w/768yNi+CEctuFE/no/kuZlBv+G8JAMiiDe19RnSvrZizo/MT/p9W7GJ+dwk3nTTD9OikJDvIzk9lRVs9znVn7fJPXIB3jv2hH93n30ppmPi6pZeXicaYv3mD0dPeanPt+d38V2akJ/n97LJySuccslvYDXRkXt15n8vqnofrKBJKXnsTSGXn8eWs5HT0qvU40tKEUjDIxjQhwbXFBZ/Z+yN8wLBIPXDEdr4bH1u7FYVPMLQj/voskLgwkEtyjNDs/3V+W+E9fnh50s1AwU0alsX5/FXURZO1wJnPvuXvuj50LqdeaWEg1pCc58WpoMnFhBbfHy4YDVVw8LTfsWocZDrtN2g/EKJb2A12Nz0lhSl4a7+43N+/uz9xD7FDtaeXi8dQ2d/D23u7NxE7Ut5I3ItH0vyPRYeeOSyaz7Ys6jtY092o7EM647BRuvaCIdreXmWPTQy6QDnYS3KOU6LBzzqQcLpgykstnjY748cZJGUnWDr6eIE676pa5t7s9/oVUsxkQdNmlGmD+vqfPy+qpb3H5s7xYOWxK2g/EKNb2A10tmTGKTUdOmZqmizRzB7hwSm7AhdUTDW0RT7Ma2TvA5FGhCyUC+e7FkyOavhysJLjH4Dd/X8xvblwUcrU9mNn5GdgUEWXt4FsYGpXua/1reHV7BbXNHXzrHHMLx4b0CFoQrNtXidOuuHBq+PUIM5z2yDdQie580zLWvIWXzMjD7dWmLm7tD+5BNvQEYrMpvnnOeDYeru3W7qCiobVXw7BwEh127rx0CgCzxoYuHQ4kLdHBu/dcxN2XRfbeG2wkuMfAYbdF/bV4+dyxvH/vJRFl7YaxGclUdGbuXq/m1x8cYcaYMxcFMcvIvMxka+v3VXF2UQ4jIvgqHorDrqQUMkZub+x17oaF47PISnHyrol591DtfkP55uJCUhLs/PpDX8dFrXW3a6dGYuXicXxw7yVRL3ImOuxRJWWDiQT3OLHbVNS14mMyk/xz7u8drOJQVRO3XVgU8cnqD+5hMvcvapspqWriUguqZAxOm00WVGPkm5ax5i1stykumZbHhgNVYb9RNbS6sNsUqRHOV2ekOLl+0Xhe21FBRX0rDa0uWl0eU5UyPSmlejXwEt1JcB+ExmQkc7KhDa9X88z7RxibkcTyuYH7toeS3tnVL9y0jHFdyaUW1LcbHHYlC6oxclmwQ7WrS2fkUdfiCnsxmoZWl6+3UBSZ77cvKEIDv/noKBWdrasHYyXKYCDBfRAak5GEy6NZv7+KzUdPcfP5RVFND5lt+/vu/iqm5KVZmik57DZZUI2RVaWQhgun5uKwqbAXiW5odUc8JWPIz0zmK3PH8MfNxzhQ6Zt7D9THXcROgvsgZHyNffyvexmR5OB6kztSe0pLdGC3qZDBva65g01Ha7k0wIVHYuG0SeYeC611Z2dG697C6UlOzp6YzVt7Toa8Tq+vl3v0ay+3XTiJ5g4P/7nuEOBbQxLWk+A+CBlfY0trW/xXeo+GUor0JAeNrcHr3P+66wQujw56ub5oSfuB2BjfeqyclgH4ytyxHK1pZmd58EvTNZpsGhbMzLHpXDBlJKW1LThsyvTObhEZCe6DkJG5J9ht3PSlCTE9V7j+Mqs/P860USN6XeggVg5pHBYT41uPVaWQhivmjCHBYWP158eDHtMYY+YO+Du6jkpPsmRTnOhNgvsglJ2aQFaKk2vOKjDV/TGUUG1/j9W2sO2LOr66IN/ysjGntPyNiZG5W7WJyZCR7GTpjDxe21ER9MPXbLvfUM6bnMPcggwm5kZ2FTBhnqngrpRappQ6oJQqUUo9EOK4ryultFIq7CWgRPSUUrxx14U8smJmzM8VKnNf/flxlIKr5ls7JQOd7QdkWiZq7s7Aa0X7gZ6uXlBAbXMHH3W21u1Kax3znDv4zuH/uflsfrlyYUzPI4ILe2YopezA08AVwExgpVKqV1RRSo0Avg9ssnqQorfRGUkR97MJJL2zeVhPWmv+sv045xTl9EmpmlM2McXE+NZjZbWM4aKpuWSlOHklwNRMS4cHt1fHHNzBV/ceyS5XERkzH/uLgRKt9RGtdQewCrgqwHGPAf8OhL5gohhQgmXu28vqOVrTzNUL8wM8KnYOm02mZWJgfDA6LayWMSQ4bCyfO5a395zkdI/dy9H0lRHxYebMyAfKutwu77zPTym1ABintV5r4dhEPzB6uvcsffvL58dJdNhYNjvypmhmSPuB2BhTWn2RuQN8dUE+7W4vb+7u3sUx2tYDov+ZCe6Bzh5/JFBK2YAngbvDPpFStymltiqltlZXh29QJPpeepITl0fT6vL473N5vLy28wRLZ46KqK1rJJwy5x6TvqqWMSwcn0lhTgp/2d59asboINpX54Wwjpkzoxzo2iS8AKjocnsEMBt4TylVCpwDrAm0qKq1flZrXay1Ls7NHdrtNgeLQLtUPzhYzanmDr62oG+mZMC4QLZk7tHy17n3URmhUoqvzs9n4+FaTnZpLy3TMoOHmeC+BZiilCpSSiUA1wNrjB9qrRu01iO11hO01hOAT4EVWuutfTJiYalAwf2Vz4+TnZrAhSYuCBwto/1AqJ2QIrgz0zJ9V8381QX5aA2vdsneJbgPHmHPDK21G7gTeAvYB7yotd6jlHpUKbWirwco+lbPC3ZUnW7j7T0nWTFvbJ+U2RmMjFN6ukfH5Z+W6bsNQEUjUykuzOKPm4/5L8UowX3wMPXu1Vq/rrWeqrWepLV+vPO+h7TWawIce7Fk7YPHmZ7uvhYEf9xUhsujueHcyC78ESkj45SKmei4/dMyfbsP8YYvTaC0toX3DvqaiTW2uVEKRiRF1/JC9B/ZoTrMdW372+H28vtNX3DJtFwm9vGV3o2eKFIxEx1jE1NfZu4AV8wezaj0RH77cSngaz2QlugwfRF2ET8S3Ie5rnPur+86QfXpdm48r6jPX9fYNi8VM9FxefumcVhPTruNvzunkA8P1VBSddqS3amif0hwH+aMy+Y1tLr47cZSJuamckGEl+uLhjEt45KKmaj4M/c+npYBWLl4PAkOGy9sLJXgPohIcB/m7DbFiCQH7x+sZkdZPTd+aUK/fOU2Mk7J3KPj6uNNTF3lpCVy1byxvLztOGWnWiS4DxIS3AUZyU52lNUzItHBNQsL+uU1jYxTgnt0jD0CfVnR1NWN502g1eXhUFWTBPdBQoK78O82vG7ROFKjvPBHpIyMU6ZlouPuo5a/wcwam8HiomxAdqcOFhLcBRnJTpSiz8sfuzIyTsnco+Pqw5a/wRgXhpFOjoODFKsKrpgzmnnjMinM6b8LJxgZp5RCRqcvW/4G8zczR/HlOaO5YErfL7iL2ElwF9xw7oR+f02nbGKKSX9Wyxgcdhv/75tn9dvridjItIyIC4e/WkYy92j01QWyxdAhwV3EhZFxumTOPSp93fJXDH5yZoi48Ne5S7VMVPrqAtli6JDgLuLCIdUyMfE3DpPMXQQhZ4aIC6mWiY3b60Up3w5jIQKR4C7iQqplYuPy6D5v9ysGNzk7RFw4pOVvTNweb7/WuIvBR4K7iAun9JaJidurZTFVhCTBXcSFQ6plYuLyeGUxVYQkZ4eIizPTMpK5R8Pt0TItI0KS4C7i4sy0jGTu0XB5vf3aekAMPnJ2iLg4My0jmXs03B4trQdESBLcRVwY88UyLRMdt9crrQdESHJ2iLg4c4FsmZaJhssj1TIiNAnuIi6MnZUumZaJiluqZUQYcnaIuFBK4bQrydyj5PZKtYwITYK7iBuHzSYLqlFyebzSfkCEJGeHiBuHXUn7gShJnbsIR4K7iBuHTUn7gSi5vFo6QoqQJLiLuHHYbdJ+IEqyoCrCkbNDxI3TpqTOPUpuKYUUYUhwF3HjsNukWiZKLq9k7iI0OTtE3DjsSurcoyQLqiIcCe4ibpw2ydyj5fZI4zARmpwdIm4cdqmWiZbLK43DRGimgrtSaplS6oBSqkQp9UCAn/8fpdRepdROpdR6pVSh9UMVQ43DbpNpmSjJZfZEOGGDu1LKDjwNXAHMBFYqpWb2OOxzoFhrPRf4M/DvVg9UDD1Om7QfiJavWka+eIvgzJwdi4ESrfURrXUHsAq4qusBWusNWuuWzpufAgXWDlMMRTItEz1ftYxk7iI4M8E9Hyjrcru8875gvg28EcugxPDgtNtwySamqPiqZSRzF8E5TBwTKD0ImG4ppb4FFAMXBfn5bcBtAOPHjzc5RDFUSfuB6GitcXs1TtnEJEIw89FfDozrcrsAqOh5kFJqKfAjYIXWuj3QE2mtn9VaF2uti3Nzc6MZrxhCHHabNA6LgtFJUzJ3EYqZs2MLMEUpVaSUSgCuB9Z0PUAptQD4Fb7AXmX9MMVQ5LQrafkbBePbjlTLiFDCBnettRu4E3gL2Ae8qLXeo5R6VCm1ovOwnwJpwEtKqe1KqTVBnk4IP4dsYoqKsU4h/dxFKGbm3NFavw683uO+h7r8fanF4xLDgK+fu2TukZLMXZghH/0ibpw2afkbDePbjsy5i1Dk7BBxI3Xu0TF29Uq1jAhFgruIG6dUy0RFMndhhpwdIm4cNqmWiYaxTiE7VEUoEtxF3Pgu1iHBPVLGOoX0lhGhyNkh4sZpV9J+IApSLSPMkOAu4sZhs6E1eGRqJiLGOoVMy4hQJLiLuDEyT1lUjYy//YBMy4gQ5OwQcWNknrKoGhmXv1pGMncRnAR3ETdG5iktCCLj9lfLyNtXBCdnh4gbp39aRjL3SJyplpHMXQQnwV3EjbEJR1oQRMYlmbswQc4OETdG5im17pGRUkhhhgR3ETdG5inVMpGRTUzCDDk7RNw4pFomKtJ+QJghwV3EjZF5SuYeGWkcJsyQs0PEjb/OXebcIyItf4UZEtxF3Ei1THQkcxdmyNkh4sbIPKXOPTJSLSPMkOAu4safuUtwj4hcIFuYIWeHiBt/4zCZlomIZO7CDAnuIm6cNsnco+Gfc5cFVRGCBHcRN/46dymFjIjLq3HYFEpJcBfBSXAXceNvHCabmCLi9nhlSkaEJcFdxI20/I2Oy6NlMVWEJWeIiBuHbGKKitsrmbsIT4K7iBt/4zCplomI26NlA5MIS84QETfS8jc6vmkZydxFaBLcRdw4pOVvVHzTMvLWFaHJGSLiRi6QHR3ftIxk7iI0Ce4ibqRaJjouj1eqZURYcoaIuJELZEfH7ZXMXYQnwV3EjVIKu01Jy98IuTwy5y7CkzNExJXDpqRaJkJuqZYRJpgK7kqpZUqpA0qpEqXUAwF+nqiU+lPnzzcppSZYPVAxNDntNpmWiZBsYhJmhA3uSik78DRwBTATWKmUmtnjsG8DdVrrycCTwE+sHqgYmhx2mZaJlMuj/RvAhAjGzBmyGCjRWh/RWncAq4CrehxzFfC7zr//GViipGWdMMFhk8w9Um6vV9r9irAcJo7JB8q63C4Hzg52jNbarZRqAHKAGisGKYYup12xdmcFW0tPxXsog8YXp1oYk5Ec72GIAc5McA+UIvRMtcwcg1LqNuA2gPHjx5t4aTHU3X7hRDZLYI/IlFFpfP2sgngPQwxwZoJ7OTCuy+0CoCLIMeVKKQeQAfR6x2qtnwWeBSguLpbv4oIbzyvixvOK4j0MIYYcM3PuW4ApSqkipVQCcD2wpscxa4C/7/z714F3tdYSvIUQIk7CZu6dc+h3Am8BduA3Wus9SqlHga1a6zXA88D/KKVK8GXs1/floIUQQoRmZloGrfXrwOs97nuoy9/bgGutHZoQQohoSbGsEEIMQRLchRBiCJLgLoQQQ5AEdyGEGIIkuAshxBCk4lWOrpSqBr6I8uEjGbitDQbq2AbquGDgjm2gjgsG7tgG6rhg6IytUGudG+6guAX3WCiltmqti+M9jkAG6tgG6rhg4I5toI4LBu7YBuq4YPiNTaZlhBBiCJLgLoQQQ9BgDe7PxnsAIQzUsQ3UccHAHdtAHRcM3LEN1HHBMBvboJxzF0IIEdpgzdyFEEKEMOiDu1LqHqWUVkqNjPdYDEqpx5RSO5VS25VSbyulxsZ7TABKqZ8qpfZ3jm21Uioz3mMyKKWuVUrtUUp5lVJxr2gId1H4eFFK/UYpVaWU2h3vsXSllBqnlNqglNrX+f94V7zHZFBKJSmlNiuldnSO7Z/jPaaulFJ2pdTnSqm1Vj7voA7uSqlxwN8Ax+I9lh5+qrWeq7WeD6wFHgr3gH7yDjBbaz0XOAj8Y5zH09Vu4GvAB/EeiMmLwsfLC8CyeA8iADdwt9Z6BnAO8L0B9DtrBy7VWs8D5gPLlFLnxHlMXd0F7LP6SQd1cAeeBO4jwCX94klr3djlZioDZHxa67e11u7Om5/iu6rWgKC13qe1PhDvcXQyc1H4uNBaf0CAq5zFm9b6hNb6s86/n8YXrPLjOyof7dPUedPZ+WdAvCeVUgXAlcBzVj/3oA3uSqkVwHGt9Y54jyUQpdTjSqky4JsMnMy9q5uBN+I9iAEq0EXhB0SgGgyUUhOABcCm+I7kjM6pj+1AFfCO1nqgjO3n+BJUr9VPbOpiHfGilFoHjA7wox8B/wRc1r8jOiPU2LTWr2qtfwT8SCn1j8CdwMMDYVydx/wI39foP/THmCIZ2wBh6oLvojelVBrwMvCDHt9g40pr7QHmd64zrVZKzdZax3XdQim1HKjSWm9TSl1s9fMP6OCutV4a6H6l1BygCNihlALf9MJnSqnFWuuT8RxbAP8L/JV+Cu7hxqWU+ntgObCkv69zG8HvLN7MXBRe9KCUcuIL7H/QWr8S7/EEorWuV0q9h2/dIt6L0ucBK5RSXwaSgHSl1O+11t+y4skH5bSM1nqX1jpPaz1Baz0B35txYX8F9nCUUlO63FwB7I/XWLpSSi0D7gdWaK1b4j2eAczMReFFF8qXZT0P7NNaPxHv8XSllMo1KsOUUsnAUgbAe1Jr/Y9a64LOGHY98K5VgR0GaXAfBH6slNqtlNqJb+pooJSF/RIYAbzTWab5TLwHZFBKXa2UKgfOBf6qlHorXmPpXHQ2Lgq/D3hRa70nXuPpSin1R+ATYJpSqlwp9e14j6nTecDfAZd2nlvbOzPSgWAMsKHz/bgF35y7pWWHA5HsUBVCiCFIMnchhBiCJLgLIcQQJMFdCCGGIAnuQggxBElwF0IIC1jd1E0p9e+djc72KaWe6iw3NU2CuxBCWOMFLGrqppT6Er7y0rnAbGARcFEkzyHBXQghLBCoqZtSapJS6k2l1Dal1IdKqelmnw7frtUEIBFfs7PKSMYjwV0IIfrOs8A/aK3PAu4B/p+ZB2mtPwE2ACc6/7yltY6oLfCA7i0jhBCDVWcTtS8BL3WZLk/s/NnXgEcDPOy41vpypdRkYAZn2nK/o5S6sPPbgSkS3IUQom/YgPrOi/Z009lYLVRztauBT40+9EqpN/BdBMV0cJdpGSGE6AOdLY+PKqWuBV9zNaXUPJMPPyIVqRsAAABySURBVAZcpJRydHbbvIgIr9YkwV0IISwQpKnbN4FvK6V2AHswf0WvPwOHgV3ADmCH1vq1iMYjjcOEEGLokcxdCCGGIAnuQggxBElwF0KIIUiCuxBCDEES3IUQYgiS4C6EEEOQBHchhBiCJLgLIcQQ9P8BfidF5fJpGD0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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": 25,
   "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 Cox-Ingersoll-Ross (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": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEbCAYAAAAh9sTfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8XOWV8PHfmabqIlmSq2y54U5zwXRMSCgbStgUCJAlgfBmE7LJ5k2BJEt209tukndDQg2kGAgloRoIxRCasGXighvYsmXJli3ZkuWiNuV5/7hzRyNpRjMjz1i6M+f7+fiD5s7V6Lni0Zkz5577XDHGoJRSKru4hnoASiml0k+Du1JKZSEN7koplYU0uCulVBbS4K6UUllIg7tSSmUhDe5KKZWFNLir40JERonIKhE5IiLzh3o8SmU7De7qeGkH/gl4dKgHolQu0OCujgtjjN8Y0zzU41AqV2hwV8OCiPxIRL481OM4Fpk4hnApa16ibUr1pcFdDTkRKQc+BdwZte1mEakRkS4RuX/IBpekWMeQJj8HvpvENqV60eCuhoPrgRXGmI6obXuA7wO/G5IRpe56+h9DOjwJLBOR8Qm2KdWLBnd13IjICuBDwN0icn3UUxcDr0bva4z5izHmceBAmscwQUQeE5FmEdkhIv8W9dx0EWkRkVOj9t0vIueFH+8UkVtFZJOItIrIfSKSH+8YUhjTR0WkOvx1hYi8JSJfATDGdAJrsH5vxNumVF8a3NVxY4y5xBgzwRhzujHm/qinFgBbM/3zRcQFPAWsAyYCHwC+LCIXhse3HfgGsFxECoH7gPuNMa9Evcw1wIXAdOAE4NtpOIYTgXUiMht4GfiJMeZ/op7fDJzU53tibVMqQoO7Gg5GA4ePw89ZDJQbY75rjOk2xtQCdwNX2TsYY+4G3gfeBsYD3+rzGr82xtQbY1qAHwBXh7cfyzGcBJQAbwI3hT+xRDscfv1E25SK0OCuhoNWYMRgvlFEXhERE+ff6312nwJMEJGD9j/gm8DYPvvdDcwH/tcY09Xnufqor+uACfGOIYWxnQTMxArYs2Mc5gjgYBLblIrQ4K6OiwGCnAHWY5U4UmaMOc8YI3H+ndVn93pghzFmdNS/EcaYS6LGWQz8ErgX+E8RKe3zGpVRX0/GOvFLrGNIZmwiMir8Oh8AfgR8S0Q8fX7mHKxSUqJtSkVocFfHhTFGgGrADfw3cJ0d6IAVwLnR+4uIJ3yy0g24RSQ/RtBL1SrgkIh8Q0QKRMQtIvNFZHHUPr8C1hhjbgSeAe7o8xpfEJFJ4aD/TeDP4e39jiFJJ2G94bRgdQb5gH+xnxSRPGAh8MJA25TqS4O7Oi7CJyiDWAHsOWPMn6Ke/gNwiYgURG37NtAB3AJcG/762xwDY0wQuBQ4GdgB7AfuAUaFx3g5cBHwufC3fAU4VUSuiXqZB4C/AbXhf98f4BiScRJW1o8xphv4Mb2z98uAV4wxe6K+J9Y2pXoRvUG2Oh5EZAnwMFad+CJjzN4+z/8QaDLG/HIoxpcMEdkJ3GiMeTHO82k/BhF5G7jBGPPuQNuU6kuDuzouRORGIAB0AkuMMV8Z4iGlLFFwV2o40bKMOl4WYJUfHgbOF5GKIR6PUllNM3ellMpCmrkrpVQW0uCulFJZ6Fj7hgetrKzMVFVVDdWPV0opR1qzZs1+Y0x5ov2GLLhXVVVRU1MzVD9eKaUcSUTqktlPyzJKKZWFNLgrpVQW0uCulFJZSIO7UkplIQ3uSimVhTS4K6VUFtLgrtJmTV0rt6/cxpq61rTtm+79Ut03me+N93rJ7nss41EqniHrc1fZZU1dK9fcXU1XIESe18XyG5eycEpJ3H2vuustgiGDzxN/3zV1rXzy7mr8wVDC/dL9s+N97yfvrqY7EMLtEj5yinWHvb/+Yw/BkIlsGz+qgMa2jn7bY+1rbwuZ1Mej1EA0uKu0qK49QFcghAG6AyGqaw/EDVKvvteEP2gtWOcfYF/7NZPdL5mf/fr7zUn97HjH2B3+OYGQ4ZE1u3s9b28Tgej1+AbaN1qq41FqIBrcVVosnTYmEtREhKXTxsTdd0S+N/K11+OKu2/0dq974P3cLiEQMnhc8fcDmD1+JACS4GfH+zkulxAMGfLDnxAArrmnGn8ghDcq815T19pve6x9AT52x5uETOrjUWogGtxVWiyYOIo8j5sOf5AZFcUDZp8H27sB8LqFP91wWtx9p5YVRb7+6cdOjLvfwiklnDljDK++t5/PL5s+4M+eONq6C96iqhJuuXhOSlnywiklnD+rgje27+ePUeNefuNSqmsPsHTamMi2hVNKYm6Pte2S+eN5buPeAX8XSqVKg7tKizV1rXT4g0wvL6K2+Qid/iD5XnfMfVfvtE4c+oOGsuK8uK+5pfFQ5OvobD+Wo11BAEYXDLxfy1HrjWXu+JGDCqQjC7yUFPp6fe/CKSUxXyvW9ljbTp8xhqc3NDJ+dKq3X1UqPu2WUWnxytYmvG7hi+fPxB80rG9oi7lfVyDIuvqDLJlaCsDmqADe16ao5xpa2gf8+fWt1vMt7f4B92sNf2roDg7uJjX2yd10ml5eDMD2piNpfV2V2zS4q7R4ZWszi6tKOecEayXS1TtbYu737u5DdAVCXHPaZFwycHDf3HiYsmIfeR4XuwYI7p3+IPsOdQHQGs7M47Ez9+7widpUdQdCeN0yqO+NJxLcmzW4q/TR4K6O2Z6DHWzdd5hlsyooLfIxo6KYmjjB3d5+5owyppYVsanxcNzX3dx4iDnjRzKppID6lo64+zW09jzX0j5wcLeDf3dwkME9A5l7WbGPkfkeDe4qrTS4q2P2ytZmAM6bZWXti6tKqKlrJRTqX/pYvbOVaWVFlBXnMWf8yLiZuz8YYlvTEeaOH8nk0sJI2SUW+zmf25U4cw8Hf/8gM3d/MITPnd4/GxFhekUx25uOpvV1VW7T4K6O2cqtTUwcXcCMCqu8sGhKKYc7A7zX1DsrD4UMNXUtLKqyTijOGT+S3Qc7ONTZv05e23yU7mCIOeNHUllaOGBZpj783JwJIyNll3haj1o/a7CZe1cghDfNwR2s0oxm7iqdNLgPQ6t2HOB/X37fEZejdwWCvLltP+fNKkfEqkUvrrJOltpdMbbtzUc42O5nUfj5ueGe8y0xSjN2Rj97/AgqSwo53BmgLc7J0vqWdvI8LmaNLY6cMI0nHTX3dJdlwAruTYe7Yr7RKTUYGtyHmTU7W7jqrmr++2/vcc091QMG+OGwJknNzlaOdgdZNqsisq2ytICKEXn96u52sLeD/+zxI4DYJ1U37z2Ez+1ienkxlaWFAHGz910t7VSWFlJalEfrUT/GxO+E6emWGXxwz8tAcLc/9dQ2a2lGpYfjgvtwCGiZ9MCqXdilavtS+ljW1LVy9V3V/Oz5rVxz98BvAvb+mfi9vbK1CZ/bxRkzeq6sFBEWV5VS0ydzr9nZQlmxj6oxVrAeNzKf0YXe2MG98TAzKorxul1Ullr93/Hq7vUtHVSWFFBa5KU7GOJodzDueI81c/cHM1WWsS7Ycno7ZLb/fTqJoy5iil64KdECUfb+fa8GHM4CwRBv72hBADv3jHc5enXt/kj22ZVgTZI1da1cfXc1gQQLcA3Gyq3NnDatlEJf76m0uKqEZzY0svtgR+Sq0NV1LSyuKo2Ub0SEOeNin1Td3HiIc2ZaJ2jtzL0+RuZujKG+pZ1FVSWUFPoAqyOmOK//1DbG9GTugy3LZKBbBqxj9LrFUXX3vn9fq3Yc4Lp7VyVc6E0dH44K7tELNyVaIGrNzhauvvttAqHEE224vAk8sqaBhtYOvvahE3h9237eqm1hTJEv5r4eV+8AM9CaJE+u3R0JZokWp0rld9HQ2s62piNctbiy33N2Xb1mZwsTT57I3rZO6ls6uP6Mqb32mzN+JA+sqouslAiw/0gXzYe7mBMu24zM9zK60BuzLNPW4edwV4DJpYWUhn9XLUe7I28I0Y50BXoWDTuGsky6u2XAWjtnypiiYRvc+86LN7bt59P3rcYfDOESobLUalcNhktiif4+VeY5KrgvnTaGPI+LzkCIkIG54T/+WH750vuRzHaggJbssrKZdrQrwP+88B4Lp5Tw+WUz+NjiSs768Urue2MH/3X5/F77GmN4en0jFSN8zJ84mpe3NA0YcHYe6KnjDrQ4Vaq/C7sFctnsin7PzR43guI8D6t3tnD5yROpqbPq74urer/e7PEj6PSH2HngaORiHvsE65zwCVeAypJC6lv797rb/e+TSgopsYN7nJOqdqcMDL7m7g9ai35lwvTyIrYPw5p7zc4WrrnnbboDViCfVFLArpb2yKfLoLVaHB8+aTwrNjRG3kCzYRG04ZL4DYajgvvCKSUs/+xSnly7h+Vv1/HAqnrOm1UR+Zhve3HTPl57f3/ksXuAlQKf37g3qWVlwf4fvZ+l08rS/j/67tdqaT7cxR3XLkREqBiRz+UnT+Dhmgb+/YMnMLqwJ4N/eUsTG3a38dN/PpEL54/jtB++yIOrd7Fg0oJ+r7u3rZM3th2gtNBHS3s3P/zIgrhjr67dn/TvAqx6e2VpAdOiFviyedwuTpk8OlJ3X72jhUKfO9IhY7Mfb248FAnudpkmOrhPLi3stRyBzc7mJ5cWUuiz1rKJ1+tuB/2yYt+gyzJdGcrcweqYeWlzU8bq+slavbOFp9btwe0Sdrd28Pf3myPzImgMLpfwsUWTePwfuwmGDF6Pi//+2EksnFLCp06v4v+99D6vvtdMvnf4ntKLFbRX7TjAi5ubmDA6n0Kfh9U7DvCXd6y19pMpAw83CYO7iPwO+DDQZIyZH+N5AX4FXAK0A9cbY95J90Bt9sJLlaUFfP+ZzTy4qp5PnjY58vy2psN8+c9rmT9xJF/70Cy+/Oe1eNzC7HH9s/yuQJCXNu+LPpYBs9pP3PkWgZAhz7ONBz6bvnp/0+FO7vp7LZcsGNdr3xvOnsojaxp4YNUuPn/eDMDK2n/54vtUlhbwkVMn4nW7uGTBeJ5cu4dvXTKHoj615t+/tZOQMdx7/SI+esdbA37sH1XQ8wbiGWCJXbDeCF7Z2swH547t9+ZqW1xVyi9efI+2dj+rd7ZyyuTRePoErRkVxbhdwubGQ3z4ROvmFZsbDzF2ZF6kzAIwqbSAFzbtIxSygovNPslaWVoQOREdr9fdDvpjR+azt60z7rENJFPdMmAF90DIsKulPfJGl2n2PJ03YSRtHX4eXdPQKzEaOyKPM6aN4bVt+yM3OPl5OJB/YvHkmKth/r+rT+GMH73Eva/t4H8+cfJxOY6BRP8tnjhpFM9v3Mu//3ktgaA1lxZMHMXetg72hpewsEWf+3JimSmZzP1+4NfAH+I8fzEwM/zvNOC34f9m1GfOnMorW5v53tObOG1aKdPLi2lr93Pj72vI97q467pFTBhdwF2fWsTH73yLnz63pV954/tPb2Z781G+cdEsHq5pYG9bJ5NKYq/M93BNPYFQcvXE6tr9XHfvqqTv9vOLF97HHwzx9Qtn99o+e9xIzp5Zxu/f3MmNZ03D53GxcmtP1m5nd1cvmcxf3tnNM+sb+XhU/ftoV4Dl1XVcOG8cp0wu4YzpY3hqXSNf/dCsmAH5nbpW3AJBAzedM23ATzCfunc1gZDhxc37WFPXGnPfRVUlGAOvvNfElr2H+OL5M/vtk+91M728qFev+6bwsgPRKksK6Q6G2He4k/Gjev4f1be0M7rQy4h8L8YYPC6JG9wPRAX3gS6KGkgms+rpFT0LiB2P4L5ifSP/9tA/IvMaoDBqJU+XwKfOqOILy2bETFbirYY5qsDLxxdX8se36vj6RbMZNyo/48cCvYP47HEjeL/pCC9s2ssdr9YSDBkkfEzRa8YFQ4bGtg5GF/jYd6gLE97nxrOmccHcCj5176pIGXjW2Phl4OEo4Sw1xvwdiL1QiOVy4A/GUg2MFpHx6RpgPC6X8N8fP4k8r4svP7SWTn+Qmx98h90HO/jttQuZEO7QWFxVyr+cXsXv36pj1Y6ew3hi7W7+WF3HTedM41/Pm8H9n15M0Bi+/8zmfj+rrd3Pi5v2YYdDAzFLETYrWBtCpqe8Ec8T/9jNQ6t28cG5Y6mK8Zo3nDWVfYe6eGbDnn5Zu23RlBJmVBTz4Opdvb730TUNHOoMcOPZ0wC49MQJ7GppZ8Pu/is2trX7eWZDI59YMpnyEXnU7o9f+62uPRCpWYdCJu7xnVw5Go9LuPu1WkKmp7+9r+hlCLoDIbY3H+kX3Cfbve4HegflXS3tkedEhJIiX9wLmaIz98GUZQJB6488E90yANPsdsgM1d3X1LVaSc5TG7ni9jf4/APvRAK7ANcuncIfblhCvteFW6zjtD+9LZxSwheWzUg6c/3MmVMJGcPv39qZsWO5feU2Vu04wHv7DvPLF9/j43e+xc+e38pHf/sm877zPFfc/ga3r9xOMHyMButE/xeWTcfnsY4x3+viN9cs5IdXLiAv6rgvnD+OJVPHsPyzS/nMWVPxuV3c+/qOyGs5QTpq7hOB+qjHDeFtjWl47QGNHZnPj688kc/9aQ2n//glWo/6+dGVC/oFka9fNIuXtuzj64+u49kvncPug+3c+pcNLK4q4WsXzgJgypgi/vXc6fzqpfe5anElZ84oA6wyyDf/uoG2Dj8//ucF7Nh/lHtf38FzG/dy8YL+72HrGw72WhFxoPLGmrpW/v3htRjgpc1NMTPgc08oZ2ZFMXf/fQcj872sb2jjJ/+8oFf2KCJctbiS7z+zma17DzNr3AiCIcO9r+/g1MmjI6954bxxfOvxDTy1bg8nThrd6+c8vnY3XYEQn1wymUAwxLPv7iUQDPUro0Dvk6IDnaAt9HmYN3EU6+oP4nYJJ08eHXO/2eNG8sTaPRxs76axrRN/0PQro0XaIVs7en0sbGjt6FXHLy30xc3cW9q78biEMUU+uoMhjDFxS0qx2G9omQruI/O9VIzIS3vHTCAY4q6/b+fnf3svUrqqGlPIdUunWJ9IwyeJP3LKxLg3GUlVZWkhF80fx/LqOm5eNqNfuTAVdkZ+6uTRFPo8PL2hkd+9tiPSmdOXAc6eWca1S6fgD4b46iPrIne/+vpFs62brswem9SNVOxPJ3PGjeBrj67n1y9v40sX9P8EOtC4+/4eX9+2n5qdLZw9szyjZZ50BPdYfx0xf+sichNwE8DkyZNj7ZKy8hF5uEVoPerH7RJOiPHRqdDn4SdXnsgn73mbf//zWlbvbMHjEv736lN7Bcl/PW86j6/dzW1PvMuzXzoHn8fFY+/s5pkNjXztwll8YrE1Zo/Lxa9XbuMzZ07lpMqegNUVCPK1R9ZTMSKfr3xwJt94bAMfmjs27v/A5dV1kT+2QDB2qUdEuPHsqXzjsQ188cF/UD7Cx5WnTur3WleeOomfPreVB1ft4j8vm8cLm/axq6WdWy7uKfWMKvRyzsxynl7fyK0Xz4nUro0xPLhqFwsmjmL+xFEsm1XBwzUNvLOrZ931aJ1+K8hdcfIErju9asAJunhKCevqD1JW7GPr3sMx97VbHjc3Hqaxzep+6XvidcLofER6X6UaDBkaWtu5cN64yLaSIm+vrphorUe7KSny4fO4MMb6fk8Ky/f6A9b/rEydUIX0rDFjB5W5E6xPRH96q449UecYXAIfW1TJF5bN4IpTJiZdbknVDWdNY8WGvTxSU8/1Z05N/A1RY186tZTS4jwee6eB376yPW7GLMAlC8bxgTlj+eZfN0SC+JcvOCFyDONHFSR1jAMd90cXTuLN7Qf41Uvvcdq00n4JTd9A/ua2/Vx//2r84Zupnz2zjCNdAd7bd4S2Dj8C3PHq9oyepE1HcG8AohudJwF7Yu1ojLkLuAtg0aJFafl8U117AGO/lxgTtxZ+xowyPjh3LM9t3AtYf6C7D3b0qgfme93852Xz+PR9q7nn9Vr+acF4vvPEuyyZWsrnzp0e2e9z503nodW7+MGKzfz5pqWR7O/2l7exdd9h7rt+MctmV1C9o4UVGxppOtxJxYjedceO7iCvbWuO1AEHyoDtskN7dxB/MMT6hrZ+x1ha5OND88by13/s5paLZ3PPa7VUlhb0CnwAl540gZe2NPHOrtZIL/ra+oNs2XuYH3zEOidx1swyPC7h5S1NMYP7U+v2MCLPw4//+cS4d1uKjKvYOinadKiLa+6pjjmZoztmGts68HlcvW6xB5DncTN+ZH6vm3bsO2Rl+fYVrPbv4b19sYNjy9FuSgt9kcy7O84nk3i6gtaVr5lqhQSYXlHEk2v3pPypwma3s3ZFlZ3OnDGG606fwq9eej8S/KLLLZkKLgunlHDq5NH87o2dXHd6VeQ6huixRp/ofHRNA//x+Lu9zgFEE+CykyZw8YJxfPnPayPH8pmzrPNDU8YUxcyU03GMIsL3rpjP2vqD/Ouf1nDNaVNYNqucyWOKeGbDHr7/9GYCIYNLoKTQFzm/A9bN0FfvbGHuhFFMLi3k3d1tGDJ/Q/R0BPcngZtF5CGsE6ltxpiMl2RsS6eNwedx9Zu0scwdP4IXNlndMcFQ7F/sslkVXDhvLL984X3uf2MnBvjFJ07uNTGL8zx8+YIT+Pbj7/Li5iY+OHcsG/e08ZtXtnPlqRMjfd9fPH8mT6zdw52v1vIfH57b6+fcvnIbzYe7+d7l8zjUGRjwI/A7uw5GztzbNe5Y+169ZDJPr2/kx89uoaaulds+PLffH9QFc8eS53Hx1Lo9keD+0Kp6CrxuLjvJ6lYZke9lcVUpr2xt6pX5g/Xp5LmNe/nQvHEJAztAZ3gpgIEmc/mIPMYU+diy9xCNbZ3MGjsiZtCd1Gfp3/qoNkhbSaEvbitka3s3JUXeSObdHQhRGPsasZjsOn1ehjP3Q50B9h/ppnxE/FsQxrLvUCfffWpjJLAL8Okzq7jt0nkALJk65rj3bN949jQ+v/wdXti0l4vm95QxV4d75/2BECJWstXZ5zzIhfPGctG88dz61/WRv+9PnWF9Ulw+Ij9jnzjiKc7z8IVl0/nqI+v59cpt/Hrltn77hIxVLv7Q3LE8+k6D1VThdvH7z5wW88bpmbwWIJlWyAeB84AyEWkAvgN4AYwxdwArsNogt2G1Qn46U4ONJZUa4TknVHDn32sT/mKvPHUSz2/cR9PhLrxuYW9bZ+QSettViyu5740d/OjZzZw1o4yvPbKekiIft0UF8allRVxx8kT+VF3H/zlnGhUjrex9e/MR7vy79UZw3elVCY9x6bQx5HkTv4GdPm0Mk0sLuf/Nnfg8Lk6I0f5ZnOfh/NkVPLNhL7ddOo/27gBPrtvDpSeN73Wf0mWzy/nhii29lg8AeHVrM4c7A1x6UnLnzM+dVcFdrw38OxeR8ElVqyxzfoyLosAK4q+93xx5bJdoKkt6gntp+IRq35ZJsDL3WeNGRDLvVC9ksi/OyVTNHXoWENvefCSp4L6mrpVXtzZR39rBs+82RsoAGKv//J/C7aWQ+eAXy4XzxlFZWsAvX3yf9Q1tuF1C3YF2/rZpb+TN0hiYPX4kH5hdwa9XboucA7jpHOtm55PHFB73QB7PvkNdkURLsC7gO3tGGT9+bktk3N+7Yj4Lp5Tw0UWVSd84PRMSBndjzNUJnjfAF9I2okFI9n90sr/YbU1HEmbKHreLWy+ew41/qOG8n69k36Eu7rxuYa+LjQC+eP4MHl+7mztereW2S+dijOG2J94l3+vm1ovnJH18yYzb5RLOmlnGA2/vwh8IcePvV8csg1x60gSefXcvb9ceYMeBo3T4g1y9pPc5kPNnV/DDFVt4ZWsT15w2JbL9qfWNlBb5Iiec0zX22eNGcN+bOwmGDLPHjYy5T2VJIfsOdUVuvl3f2oEIkc4osDL3kIFDnf5+/y9a2/2UFPoimXeqHTP2/pm8wCj6lnuJsro1O1u46u7qyJvOOTPL+P4VC2g+0jVsrqp0u4QL5ozlvjd2smWv1e46qsDDaVNLeXP7AULhi6D+48NzWTilhDNmlA2bQB5L30TL7iA6sXJ00uM+XsfjqCtU0yGZX2yymXJJoReXWO/mLoGy4v6ZVlVZER85ZSLL367jc+dOo3pHC29sO8D3Lp+X0sfuZCeEvRbNQGWQZbMqKPK5eWr9Ht7dfYjZ40ZwcmXvTpbp5cVMKilg5Zae4N7eHeDFTfu4MnzxVDrHPmf8yMhJs75tkDa7tt7Q2sGMimLqW9qZMKqgVyYdvb5MdHAPhgwH27spLYqquQ8yuGcycx83Mp9CnzvhXZn2Herkq4+ujwR2l8Bp08YweUwhk8cUDptgCFbfu80l8Nmzp3Hz+TNT6p0fLuIlK8Nx3DkX3JORbLZZHdU3LxC3Fv7F82dYJzofW8/bO1qYVlbEJ6Oy4XQ6b1YFdycogxT43FwwdyyPrmnAHzTccGZVv5N3IsL5syt4pKYhkim/tLmJDn+QS0+a0O81j1V0QI93C7zJkXbI9khw73vRmb2+TN9e90MdfkLGyuztNyZ/MLVz+pluhQTr09e08tgLiNnBMBAy/O71HXR0B/C4BBMuwQzXtVzOnlnOHa9uj8zJ06dbn/qGY0BMhlPGrcE9jmQz/GRO5k4ZU8S5M8t4ObzQVsPBDtbWH8zIBEn2jWnueKu3HOBPb+/ikhMnxMzw//BWHW/vaOHcE8p5at0exo7Mi3sx0rE4HHUHopv+VBOznGT3utsdM/Wt7ZFlgW2lhXbm3rsd0l5XJh2ZuzeF9snBmF5e3G899L5dMCeMLeaOa8+gtd0/bEow8RzPOrPqocH9GKQyaWdUjIgE92CcnvZ0jivRa3cFem5oEa/H/vTp1iqcK7c0cXLlaF7Z2sy1S6f068BJh5qoYBa3q6Y4D5/Hxa6Wdjr9QfYd6uq3tG9JkVUC6NsxYz8uLfL1rBcSjH9Tj1jszD1Ta8vYppcX8+S6PXR0BykIL4YWvcCdYJ03mRauzzshWDol280mw3fZNodI9rLsC+ePIy98yfNw+Ah95ozyhOPJ97o5Y/oYVm5t4vmNe+kOhrjs5PSXZKBnOeeBxuNyCZUl1rrhDeHlf6N73CGq5t6nLNMSFdztzLs7kFpZxi4X+dyJW0CPxfTyYoxiua6TAAAYA0lEQVSB2v1WaWbL3kM8UmNdBO4SyPO6OGN6cie0Ve7SzP04WTilhAc+O3w+miY7nvNnV7DyiY3c8ep2KksLOGnSqCEdT2VpIbta2mP2uAMUeN3keVz9M/dwsC8p8kU+taTaCmnv7/VkuCxT0bPGTEd3kM/cv5oCn5tffPxk9rR1DIv5o4Y/De7H0XD7aJrMeM6bVQFspLb5KFeeMnFQV02mczyTSwtZU9fas9RvSe/gLiKUFvVfX8auwZcW+mg5Mrhb7UW6ZTK81nrVGCu4//z5rTS2dVBZUsgfbljCpJL+d5dSKh4ty6gBNR3uiiwe9MyGxiG/8XFlSSGHOwO8u7uNPI8rZjtpSWH/lSFb27vJ97oo8LkjmXeqt9o7Hq2QABv3HEKwLtIKhAz/8eE5GthVyjS4qwFFL+drn3gdSnaN/Y1tB6gsLYz5SSJ25t4d6aTxDfYipuDxydyjf8cuYFPUWvdKJUuDuxqQfUHXcDkRbHfH7D7Y0a/ebrPWdO/dCmmvCAkM64uYYPj9zpUzac1dDWi49ShHtz5WxrlrVmmht3/mHr46FaIy90GeUM10cB9uv3PlTBrcVULD6UTwyHwvowu9HGz39+txt5UU+Wjr8Pe62Ujr0e7IydfBZu7+47C2jG04/c6VM2lZRjmOHaTjBXc7Qz/Y0VOaaTkalbkPclXI7qC1PK0nAxdxKZVuGtyV4xTnWRcRHekMxHy+JHzi1O519wdDHOoMRLbbZZl4a9jE0x0I4XO7MtoOqlS6aHBXjrKmrpVVO612zG/9dUPM1szolSEBDoZPrpaGlyZwuwSRwWXume6UUSpddKYqR6muPYAJ3xjZH6c1s7TPypDRV6eCdaGTz+0aVLdMpk+mKpUuOlOVo9grcQ7UJtiTufvD/w2vKxO1vrvP40o9c9fgrhxEu2WUoyTTJji6MLwypJ25H+2duQODytz9QQ3uyjk0uCvHSdQmmOdxU5zniWTs0Wu523yeQZRlgqHj0gapVDroTFVZqaSo50Ime6EwO6MHK7gPZm0ZPaGqnEJnqspKpYW+Xpl7cZ6HPE/POuxe9yBq7kGjZRnlGDpTVVay1pfpqblHl2RgcDX37kBQM3flGDpTVVbqnbn7e51MBass06WtkCqL6UxVWamkyBfpkmk92k1pVL0drMw91Zq7X8syykF0pqqsVFrk42h3kE5/kJao5X5tg+qWCYQi919VarjT4K6ykr2OzMF2P63t3b0uYIJBXsQUDOHzZPbm2EqliwZ3lZXsdWT2tHXQ3h3sl7l73YI/YFJ6TW2FVE6S1EwVkYtEZKuIbBORW2I8P1lEVorIP0RkvYhckv6hKpU8O3OvbT4K0L9bxuMeZOauZRnlDAmDu4i4gduBi4G5wNUiMrfPbt8GHjbGnAJcBfwm3QNVKhV2MN/efAToCfa2QS8cppm7cohkZuoSYJsxptYY0w08BFzeZx8DjAx/PQrYk74hKpU6uwyzvckK7v0zd9GFw1RWS2amTgTqox43hLdF+0/gWhFpAFYAX4z1QiJyk4jUiEhNc3PzIIarVHJGF1g1dztzt2vwtsEuHKZryyinSGamxioy9j0TdTVwvzFmEnAJ8EcR6ffaxpi7jDGLjDGLysvLUx+tUknyuF2MKvBSd6AdiFGWSbEVMhQyBELa566cI5mZ2gBURj2eRP+yyw3AwwDGmLeAfKAsHQNUarBKi3wEQgYRGFXQO3P3pngRk13C0eCunCKZmboamCkiU0XEh3XC9Mk+++wCPgAgInOwgrvWXdSQKglflTqqwIunTznF53ERCBlCoeTaISPBXcsyyiESzlRjTAC4GXge2IzVFbNRRL4rIpeFd/u/wGdFZB3wIHC9se+FptQQsU+i9r2ACXoy8GRPqtolHM3clVMkdbMOY8wKrBOl0dtui/p6E3Bmeoem1LGx6+x9L2CCngy8Oxgi35v4qtNIcNfMXTmEzlSVtezMve/JVIjK3JM8qWrX57VbRjmFzlSVteyMvW8bJERl7kkGdy3LKKfRmaqyViRzj1GWsTPwZDtmujS4K4fRmaqyln0idWvjYdbUtfZ6brBlGQ3uyil0pqqs1XS4C4BX32vmmnuqewV4O0gnezcmPaGqnEZnqspaew91ANbl1P5AiOraA5HnfCmWZfQiJuU0OlNV1jr3hAryvS7cAl6Pi6XTxkSeG3RZRjN35RBJ9bkr5UQLp5Sw/MalVNceYOm0MSycUhJ5brAXMWkrpHIKDe4qqy2cUtIrqNtSLctot4xyGp2pKid5U+xz9wet1TTyNLgrh9CZqnLSYLtltCyjnEJnqspJdgZuZ+SJdAeCgJZllHPoTFU5abBlGQ3uyil0pqqc1NMKGUxq/+7IwmGxbkym1PCjwV3lJF+KZZkuvUJVOYzOVJWT7Aw8lT53n9uFiGbuyhk0uKucZGfgyXbL+IMhLckoR9HgrnKSiOBL4SbZ3YGQnkxVjqKzVeUsr1tSulmHBnflJDpbVc7yeVwpLRymFzApJ9HZqnJWKsG9K6iZu3IWna0qZ3lTrblr5q4cRGerylk+j4uuJIO7XzN35TA6W1XO8rmTL8to5q6cRmerylk+j7ZCquyls1XlrFQydy3LKKdJaraKyEUislVEtonILXH2+biIbBKRjSLyQHqHqVT6pdQtE9BWSOUsCW+zJyJu4Hbgg0ADsFpEnjTGbIraZyZwK3CmMaZVRCoyNWCl0sXrdnG0K5DUvt2auSuHSWa2LgG2GWNqjTHdwEPA5X32+SxwuzGmFcAY05TeYSqVfj6PK6W1ZfI0c1cOksxsnQjURz1uCG+LdgJwgoi8ISLVInJRugaoVKb4PK6UVoXUsoxykoRlGSDWUnh9F8H2ADOB84BJwGsiMt8Yc7DXC4ncBNwEMHny5JQHq1Q66cJhKpslM1sbgMqox5OAPTH2ecIY4zfG7AC2YgX7XowxdxljFhljFpWXlw92zEqlRWrdMkaDu3KUZGbramCmiEwVER9wFfBkn30eB5YBiEgZVpmmNp0DVSrdUumW0bKMcpqEs9UYEwBuBp4HNgMPG2M2ish3ReSy8G7PAwdEZBOwEviaMeZApgatVDpYa8skvs2eMUa7ZZTjJFNzxxizAljRZ9ttUV8b4Cvhf0o5QrKZu/0GkKfBXTmIzlaVs+xuGSs3ic/uqNHb7Ckn0eCucpYvHKwTlWbs7F4XDlNOorNV5Sy7hp6o1z0S3D3ujI9JqXTR4K5ylp2JJ6q7+7UsoxxIg7vKWXYmnuhCpq5I5q5/Lso5dLaqnGVn4okyd/t57ZZRTqKzVeUsOxNPtHhYT1lG/1yUc+hsVTnLzsQTlWXsE65allFOorNV5SxvkidUtRVSOZHOVpWzkm6FtMsymrkrB9HZqnKWnYn7NXNXWUhnq8pZdibeleRFTNoto5xEZ6vKWalexKQnVJWT6GxVOSvpbpmAtkIq59HZqnJW0t0ymrkrB9LZqnJWpFsm2ROqGtyVg+hsVTkr1VZI7ZZRTqKzVeUsvYhJZTOdrSpn5SWZufuDITwuweXSJX+Vc2hwVzkrlcxdO2WU0+iMVTnL7RLcLkmqFVJPpiqn0RmrcprP7UqiFdJocFeOozNW5TSfJ4ngHgjpyVTlODpjVU7zul10B82A+3QHtSyjnEdnrMppeUll7kHN3JXj6IxVOc3ncSXRCmnwerQNUjmLBneV07xuSWo9d83cldMkNWNF5CIR2Soi20TklgH2+6iIGBFZlL4hKpU5yWTu2gqpnCjhjBURN3A7cDEwF7haRObG2G8E8G/A2+kepFKZklwrpF7EpJwnmRm7BNhmjKk1xnQDDwGXx9jve8BPgc40jk+pjLK6ZRJn7noXJuU0yczYiUB91OOG8LYIETkFqDTGPJ3GsSmVcUn1uWsrpHKgZGZsrDaBSGOwiLiAXwD/N+ELidwkIjUiUtPc3Jz8KJXKkGRaIf1allEOlMyMbQAqox5PAvZEPR4BzAdeEZGdwFLgyVgnVY0xdxljFhljFpWXlw9+1EqlidftSm5tGQ3uymGSmbGrgZkiMlVEfMBVwJP2k8aYNmNMmTGmyhhTBVQDlxljajIyYqXSSLtlVLZKOGONMQHgZuB5YDPwsDFmo4h8V0Quy/QAlcqkZLtlNLgrp/Eks5MxZgWwos+22+Lse96xD0up48Pn0bKMyk46Y1VO87pddA2QuRtjNHNXjqQzVuW0RN0ywZDBGL1/qnIenbEqpyUqy9gnW72auSuH0RmrcprX7SJkIBAnwNtZvWbuyml0xqqcZtfS47VD2tu15q6cRmesyml2Ru4PxL4bk2buyql0xqqcZtfSu4LBmM9Hgrtm7sphdMaqnJYXzsjjdcxoWUY5lc5YldPsoO2Pc5Nsu1yjC4cpp9EZq3KaN2HmbpVrNHNXTqMzVuW0SLdMnODepSdUlUPpjFU5LVErpF2u8Xli3dZAqeFLg7vKaV63FbTjlmUimbv7uI1JqXTQ4K5yWl6ii5i0FVI5lM5YldPsjDxe5m6vO2Nn+Eo5hQZ3ldO84Vp6vMXDNHNXTqUzVuU0X4JWyC69iEk5lM5YldMStUL6tRVSOZTOWJXTIpm7rgqpsozOWJXTEmXuuiqkciqdsSqnJb6IKYRLwKPBXTmMzliV07yR9dzjZ+66aJhyIp21Kqd5XIJI/My9KxDSertyJJ21KqeJCD63a8CLmPI0uCsH0lmrcp7P7Rpw+QEtyygn0lmrcp7PEz9z7w5qWUY5k85alfMGCu7+YEjbIJUjJTVrReQiEdkqIttE5JYYz39FRDaJyHoReUlEpqR/qEplhtftGnBtGS3LKCdKOGtFxA3cDlwMzAWuFpG5fXb7B7DIGHMi8Cjw03QPVKlM8Xni19y1W0Y5VTKzdgmwzRhTa4zpBh4CLo/ewRiz0hjTHn5YDUxK7zCVypyBumW6Nbgrh0pm1k4E6qMeN4S3xXMD8GysJ0TkJhGpEZGa5ubm5EepVAZZmbuJ+ZzW3JVTJTNrY92lIOZfgohcCywCfhbreWPMXcaYRcaYReXl5cmPUqkMsjL3YMzntFtGOZUniX0agMqox5OAPX13EpELgG8B5xpjutIzPKUyz+dx0d4diPlcd0Azd+VMycza1cBMEZkqIj7gKuDJ6B1E5BTgTuAyY0xT+oepVOb4PC78ccsyBq9m7sqBEs5aY0wAuBl4HtgMPGyM2Sgi3xWRy8K7/QwoBh4RkbUi8mScl1Nq2PG6ZeATqpq5KwdKpiyDMWYFsKLPttuivr4gzeNS6rjxedzaCqmyjs5alfMSLRzmc8fqKVBqeNPgrnKezyMDLhymmbtyIp21KucNeBGTtkIqh9JZq3Ke1S3TP7gHQ4ZgyOBzu4dgVEodGw3uKud542TudsD3erTmrpxHg7vKeT6Pi0DIEAr17nXvCgd8bYVUTqSzVuU8u6be96SqnbnrbfaUE+msVTnPzsz7Bne7VKPruSsn0lmrcl4kcw/EDu7aLaOcSGetynmRzD0QuyyjwV05kc5alfPsskvfdsguLcsoB9NZq3Je3LKMZu7KwXTWqpxnB++uvmWZ8OM8zdyVA+msVTnPF6cs0x25iEn/TJTz6KxVOS9ht4xm7sqBdNaqnBfvIiZthVROprNW5bx43TKRsoxm7sqBdNaqnBevz91+rMsPKCfSWatyXrxuGW2FVE6ms1blvJ5umd6rQvr1IiblYDprVc7Ti5hUNtJZq3JeT3AP9tqurZDKyXTWqpzndVt3WupblukOP7afV8pJNLirnDdQn7vP7UJEg7tyHg3uKufZZZd+3TKBkNbblWPpzFU5T0TwuV39LmLyBzW4K+dKauaKyEUislVEtonILTGezxORP4eff1tEqtI9UKUyyeuWmBcxab1dOVXC4C4ibuB24GJgLnC1iMzts9sNQKsxZgbwC+An6R6oUpnk87hitkJq5q6cKpmZuwTYZoypNcZ0Aw8Bl/fZ53Lg9+GvHwU+IHoWSjnMuvqDrKlrjTzed6iTo12BXtuUcgpPEvtMBOqjHjcAp8XbxxgTEJE2YAywPx2DVCqT1tS1crDdT2t7Gx/97ZtMLCkAoKG1A4Br7qlm+Y1LWTilZCiHqVRKksncY2XgZhD7ICI3iUiNiNQ0NzcnMz6lMq669kDkawMU+twU+tyRbf5AqNc+SjlBMsG9AaiMejwJ2BNvHxHxAKOAlr4vZIy5yxizyBizqLy8fHAjVirNlk4bQ57XhVsg3+viR1eeyI+uPJH88Davx8XSaWOGephKpSSZssxqYKaITAV2A1cBn+yzz5PAvwBvAR8FXjbG9MvclRqOFk4pYfmNS6muPcDSaWMi5ZdY25RyioTBPVxDvxl4HnADvzPGbBSR7wI1xpgngXuBP4rINqyM/apMDlqpdFs4paRfAI+1TSmnSCZzxxizAljRZ9ttUV93Ah9L79CUUkoNljbxKqVUFtLgrpRSWUiDu1JKZSEN7koplYU0uCulVBaSoWpHF5FmoG6Q315Gdi1tkE3Hk03HAno8w1k2HQskfzxTjDEJrwIdsuB+LESkxhizaKjHkS7ZdDzZdCygxzOcZdOxQPqPR8sySimVhTS4K6VUFnJqcL9rqAeQZtl0PNl0LKDHM5xl07FAmo/HkTV3pZRSA3Nq5q6UUmoAjg3uIvI9EVkvImtF5G8iMmGox3QsRORnIrIlfEx/FZHRQz2mwRKRj4nIRhEJiYhjuxkS3RjeSUTkdyLSJCLvDvVYjpWIVIrIShHZHJ5nXxrqMQ2WiOSLyCoRWRc+lv9K22s7tSwjIiONMYfCX/8bMNcY87khHtagiciHsNbBD4jITwCMMd8Y4mENiojMAULAncBXjTE1QzyklIVvDP8e8EGsm9GsBq42xmwa0oENkoicAxwB/mCMmT/U4zkWIjIeGG+MeUdERgBrgCuc+P8mfK/pImPMERHxAq8DXzLGVB/razs2c7cDe1gRMW7r5yTGmL8ZYwLhh9VYd7xyJGPMZmPM1qEexzFK5sbwjmGM+Tsx7o7mRMaYRmPMO+GvDwObse7j7DjGciT80Bv+l5ZY5tjgDiAiPxCReuAa4LZE+zvIZ4Bnh3oQOS7WjeEdGUCymYhUAacAbw/tSAZPRNwishZoAl4wxqTlWIZ1cBeRF0Xk3Rj/LgcwxnzLGFMJLAduHtrRJpboeML7fAsIYB3TsJXMsThcUjd9V0NHRIqBx4Av9/kk7yjGmKAx5mSsT+tLRCQtZbOk7sQ0VIwxFyS56wPAM8B3MjicY5boeETkX4APAx8Y7vegTeH/jVMlc2N4NUTC9enHgOXGmL8M9XjSwRhzUEReAS4CjvnE97DO3AciIjOjHl4GbBmqsaSDiFwEfAO4zBjTPtTjUT03hhcRH9Z9gZ8c4jEpIich7wU2G2P+Z6jHcyxEpNzujBORAuAC0hTLnNwt8xgwC6srow74nDFm99COavDCNxfPAw6EN1U7tftHRD4C/C9QDhwE1hpjLhzaUaVORC4BfknPjeF/MMRDGjQReRA4D2vlwX3Ad4wx9w7poAZJRM4CXgM2YP39A3wzfK9nRxGRE4HfY80xF/CwMea7aXltpwZ3pZRS8Tm2LKOUUio+De5KKZWFNLgrpVQW0uCulFJZSIO7UkolKd0LsInIcyJyUESe7rP9fhHZEV4Yca2InJzqa2twV0qp5N2PdZFRuvwMuC7Oc18zxpwc/rc21RfW4K6UUkmKtQCbiEwPZ+BrROQ1EZmdwuu9BBxO9zhBg7tSvYjIfSLyYREZLSLPhi/IUmogdwFfNMYsBL4K/CZNr/uD8P0dfiEieal+87BeW0apIbAA69L2J4BbjTFvDvF41DAWXrzsDOARa1UEwLrSHBG5Eoh1tenuJK7YvhXYC/iw3jy+Eee14tLgrlSYiLiAOVgL0d2ugV0lwQUcDK/q2Et4QbNBLWpmjGkMf9klIvdhfSJIeWBKKctMrJUfrwc+F155UKm4wksN7xCRj4G1qJmInHSsrxu+25S9SNoVDGKVSA3uSvVYgHWzhJex/pg+NcTjUcNMeAG2t4BZItIgIjdg3SzoBhFZB2wkhTt2ichrwCPAB8KvZ5drlovIBqzF0cqA76c8Vl04TClL+ObE+4wxvxGR04E/ArOjbn+olGNocFdKqSykZRmllMpCGtyVUioLaXBXSqkspMFdKaWykAZ3pZTKQhrclVIqC2lwV0qpLKTBXSmlstD/B/4YNzJrClqsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAA0AAAASCAYAAACAa1QyAAAABHNCSVQICAgIfAhkiAAAAHZJREFUKJFjYKACCGFgYJjMwMBwmIGB4RMDA8N/BgaGJYQ0XYAq/MzAwHCdWE2ODAwMqgwMDIwMDAwOuDSxoPH3EzKVgYGBgYkYRaOaBlwTeuQGQDEDAwODBJS2ZGBgWABlv2FgYChBN6SBAZJ0cOEH5LiMzgAA6XoX52TB9a4AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$1$$"
      ],
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 27,
     "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": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd0XNW97z/7TFGvlrtlyb1QDBYuJHHoCbnBIckNBEMCuQRIY73k3tQLpJEQ8m7a4wYSMOBQYhswphpTYqrBliXLvRdZo+KiNlazpGn7/bHPmSKNZAm10Wh/1vKaOXvOnNm2pfOdXxdSSjQajUaj+bgYQ70BjUaj0QxvtJBoNBqNpk9oIdFoNBpNn9BCotFoNJo+oYVEo9FoNH1CC4lGo9Fo+oQWEo1Go9H0CS0kGo1Go+kTWkg0Go1G0ye0kGg0Go2mT9iHegMDiRBiKbA0LS3t9pkzZw71djQjhYMH1eOsWUO7D42mj5SUlNRKKUef7TwxEnptXXTRRXLr1q1DvQ3NSOHSS9Xje+8N5S40mj4jhCiRUl50tvO0a0uj0Wg0fUILiUaj0Wj6RFwLiRBiqRBieUNDw1BvRaPRaOKWuBYSKeWrUso7MjIyhnorGo1GE7fEtZBoNBqNZuAZlkIihJgjhHhYCPG8EOI7Q70fjUajiUVKXG4eevcIJS73gH7OoNeRCCFWANcA1VLKc8PWrwYeAGzAY1LK33d1DSnlfuDbQggDeHSAt6zRaDTDjhKXm2WPFuLzB3DaDVbetpiCvKwB+ayhsEieAK4OXxBC2ICHgM8Bc4FlQoi5QojzhBDrOvwZY77nC8CHwNuDu32NRqOJbYqO1fGfz+7A4wsQkOD1BSgsrRuwzxt0i0RK+YEQIr/D8kLgiJSyFEAI8QxwrZTyfpT1Eu06rwCvCCFeA1YN3I41Go1m+PDntw7yv+8cCR4bAhx2g8VTRw3YZ8ZKi5SJQEXYcSWwqKuThRCXAl8GEoD1XZxzB3AHwOTJk/trnxqNRhOTFB2r477X9rOzMlTuIIBPTs/hB1fOHDC3FsSOkIgoa132bpFSvge8190FpZTLhRAngKVOp7OgT7vTaDSaGKXE5eYfHx3jtV0nOt00bYYYcBGB2BGSSiA37HgScLyvF5VSvgq8etFFF93e12tpNBpNLFHicrO2pIJnt1biD3T+3m03BPdee+6AiwjEjpAUAzOEEFOAKuAG4Ma+XtTq/jt9+vS+Xkqj0WhihlVbyvn5y3s6CYgAHDbBdRfl8uX5kwZFRGBo0n9XA5cCOUKISuCXUsrHhRB3Am+i0n9XSCn3DvbeNBqNJpaxrJDVxRWEN24fKgGxGIqsrWVdrK+ni8C5RqPRjHRKXG5ufLSQdl8gYt0m4IaFk4dEQCyGZWV7T9G9tjQaTTxQUlbPz9bu6iQidkPwmy+ex31fOm/IRARiJ0YyIPRHjGRLaR1FZfV8YlrOkP5HaTSakcnD7x3l/75xIJiRZQglIEPlxopGXAtJX7O21hRX8OO1uxDAQ44j/OKac3Cf8bB46qiY+M/TaDTxS0lZPX986yCbS+uDa4NVF9Jb4lpI+srh6mZAFbR4vAHueWk3UkKCY2D71mg0mpFLicvNykIXL+6oouMk9MGqC+ktcR0j6etgq0/NyFHXAQxDEJBKVAa6b41GoxmZrNpSznUPb+KF7Z1FZDDrQnpLXAtJX4Ptn5qegwAWT83m3muDjYoRQgT71gxWm2aNRhO/lLjc/PC5Hdz94m4CHdJ6nTbBTYsm8+y3LubGRbHZ7imuXVt9DbYbhiAj2cGMsWksW5jLz1/ejT8AKQk2Nh2t5eDJJn71yl68/oB2d2k0mo9FicvNDcs34/VHmiCxkNbbU+JaSPqjRUpmkoPTZ7y0ePz4AzBzbCqHTjXzp7cOIQg1BPOY7q5Y/w/XaDSxwwcHa/jR8zs7iYjlxopVC6QjcS0k/cFi+2Hmn9xPy5EvAZA/KoVDp0JBeIuAhCp3KyUutxYTjUbTLSUuNw++c4h3D9YG12IxrbenaCHpju2ruL/hxwQQiBdWMV/8jAsmz+KdA9X4TEemACZlJVHhbmV1UTkvbK/ULi6NRtMlKz4s5Tev7Y8IphvEZlpvT4nrYHtfs7ZwHwPAhoSAl8XGfi7KU4F3uyEwhEoF/vSM0YDO6NJoNF1T4nLzzSeKuHfd/k59spwOY9iKCMS5RdLnGMmMq/Bv/DN26QMpqZepZCU7uHHRZGaNS6OwtC6YvbWquBwpB34SmUajGV6UuNw8W1zO8yWVdOz2PpwC6t0R10LSZ3IX8s6UH3Pl0fsRIsAvHU/TdvrLMPZTFORlRfzHXzxlFPtPNvLYLQuG9Q+ERqPpP7pq9w7DL6DeHcNWSIQQKcAHqDb06wbqc7JFMxKBgcSBD8fJQpj1qU7nnTcpg63lbi7MzRyorWg0mmFCicvNc8XlPLe1MiIpZ6jbvQ8UQzGPZAVwDVAtpTw3bP1q4AHUPJLHpJS/P8ulfgo8N2AbNWkevxjfkb9jw4cfGwlTPx31vElZSXh8AWqb2ymtbebDw3VcNntM3PygaDSanlHicrPs0UI8MdjufaAYCovkCeBB4ClrQQhhAx4CrkKN3S0WQryCEpX7O7z/VuB8YB+QONCb9U9cwN3eW/mjczmPO27gu7kLoaIIyjZC/hLIXQjAXP9Bvmt7mR2b/HznAzv+gOTRjaWsul1ncGk0I4XC0lp+tGZXJxGJJzdWNIZisNUHQoj8DssLgSNSylIAIcQzwLVSyvtR1ksEQojLgBRgLtAqhFgvpQx0PK8/yEx28H5gnjpISFMi8sTnwe8Fwwb/9icYO5cL3/068+weAlteYp78b7Yxk3ZfgL/86xD/edXwzcbQaDQ94w9vHuChd48Gj4dzXUhviZUYyUSgIuy4EljU1clSyrsBhBDfAGqjiYgQ4g7gDoDJkz/+t4CMJCd1ZOCXgvG2Bti/Dvwe9WLAB+t/CBd+DeFvxy7AL1Wa8Db/TAA+PFLLlmN1XD8Cfpg0mpFI0bE6frNuH7urGoNrsdrufaCIFSERUdY6pzl0PEHKJ7p5bbkQ4gSw1Ol0FnzcjWUmOwhgUEMm44QbMs6PPCEQgPbmYLsUL3YKA3OYOz6dfSfUD5bXL1m5pZy126IXK5a43MFU4pHwQ6fRxAMlLjcrPixl/e6TnW5WsdrufaCIFSGpBHLDjicBx/t60f7otZWR5ADglMxilKwHm1O9IGwg/WCzQ+4C2PM8AviB5ztskzN57DMzuXPVNtq9geAPmbdDP64Sl5tVW1y8vOM4ASlx2nXjR40m1ilxuXm+pILntlZ2m9Y7kn6PY6WyvRiYIYSYIoRwAjcAr/T1on2ubAccNoPUBDvVMotMfx3Ul4ItAZatVidc9E1wJAfPzxGNTMxM4so5Y1l522KWLVL6KIgsVixxuVm2vJC126rwBSQBCW3eAH9+66BuSa/RxCirtpRz/SObWV1UESEiw6Xd+0AxFOm/q4FLgRwhRCWqDuRxIcSdwJuoTK0VUsq9g723rshIclDdnEmapxTqj0H2FJjxGXCmARLcLhAGDY7RfNq/C3fuzQDBosWi0nq8fsmfv3pB8FtKYWkdHn/n/ICPjtZR8lihtkw0mhiixOVmbUkFq4srOrU3ice6kN4yFFlby7pYXw+sH+Tt9IjMZAenmrJI8rqh5iDkzAAhICtfCUtiBqRP4qTI55L2TZSzD5gffP95kzLZ0qHF/OIp2cHnTptgzvh0dlYqy6mjC0yj0QwdJS43Nz5aSPsIqgvpLbHi2hoQ+joh0SIz2cEpzB+UusOQNUU9z85XjR1Pu2ixpTHFvYkE4ePrh77PgeINwfdPyUnheEMbrR5/cC0rRcVarpo7ltV3XMwvlp6DzVA5B7pfl0YTG2wtq+cnz+/sJCJ2Q/CbL57HfV86b8SLCMS5kPRHjAQgM8nJKRnW+iTbFJKsKcqt5S6jxePHQAmFAx/ufe8ET586OgWAsrqW4NqHR9Qcgns+PyfoArt9ibrun66bp384NZoh5m/vHua6hzdztEb93hpiZMdBuiOuhaS/LBKPL0C1DLuxZ081H6eAvx2aTyEmzMOLI9jdM2vO5cHTp+QoISmtCQnJweK3+c/EV2kt3Rxcu+b8CQCdOoRqNJrBo6Ssnq8+son/efNQMOPSqgtZfcfF2gqJQqyk/w4IfZ3ZDso/+u7BajIihMSySPKDS6PPuZwDM67HseVBptW+y+zJY4OvWUJyrFZNVtxX+Dq/rvsRAoln3VrK3L8iP6mNWbmfxGkz2FPVwNJ5Ez72njUaTe8pcbn552YXL+2sigiow8irC+ktcS0k/VFHUlhaR0BK6knDJw0MITGaTimrxIqVAGROZnbeJ2DOPPjjdHjzbrj8HshdSLLTzviMREprlUXS/sH/YhfK5+qQXiZvugcI4LA5eSD1CjaXfgaYE7GPN/ac4GhNM4un5gT3pQsYNZr+YVWhi3te3hPVGzAS60J6S1wLSX+weOoonHaDc/2HMAioEvynvwS3vAIT5oNhV61SMvPUG9zHQBhw7H14skidl7uQTyYcZdrBtRSunsT5LVuRqMQvA4mwDGi/h8/6X+fStrfZXzSZOQuvBOCl7VX84NkdZu+ew0jAH9AFjBpNXylxuVlZ6OKF7VUR6zqtt3fEtZD0h2urIC+LlbctxvPeRsQxs5eL36O6/+YuhJTR0HwKTldAxkS1LkPCQNlGDpxs5LcNd+HEizigBMQrDewEMMKaw0hU0MqBj/WvruHM2AIK8rJ4bfcJQMVOvH4ZtVJet1nRaHpHicvNDcs34/VHmiE6rbf36GB7DyjIy+Liy7+IsCeq1ig2p2ohX1GkREQG4OkvquP8JaE2KoYB+Utw73sHJ94I0TA6jruZqNqBBaTq17XJPzs4+91vFi5a35IsrDThEpebrz6ymT+8eZCbHivUlfEazVl492A133p6aycR0Wm9H4+4FpJ+JXch3PIqXH530F1F2cbQ6+FWys0vgz0Bpl0BuQvJmn0JglA2lk8KfNiQhlMJkz0Rrv497Wn5HJMTuMlzF7vErGAtSYW7FYCbJ53kvQVbWJJYisMmWPnNRRTkZVFYWofPvLhlpWg0ms6UuNzc/PgW/uMfxdQ2qy7eOq2378S1a6vfyV0YHGQFmNZHghIRy0oByLsYpl8FJ3cBMHt0Igg4NOZznB69gMCZOrLmXs7scekRA7ISc/LIdjayrWom3/7kFArysmhq83Kkppn54hC/qPsttlofj0onNwbuYvrYqwAiihftNl3MqNFE4/GNpfz2tf0RvgCDkdXufaCIayHpjxhJt+QuVNZJh2mJAEy9FA6sgzfvgdPlYE9i1m0rwJnc+RoWydlkNKimx+nJquvw7qoGpITFtv0Y0geAAzXz5ERDKxlJDqaPSQ1e4juXTtNt6jWaMEpcbh585xDvHqyNWBeA02FoEekH4lpI+iP996x0tFIskswfzM0PAhIyJsOpPdHPDb4nG9FWj9Nm0NDqBWBHxWkADiWeB/5nAZAICgNzWNDQxuxx6ZTXnQleos0b2crB6hPk9Qd0lpdmRFHicvNssYvnS6o6pfXqgHr/EtdCMqS4XeYT8ye4oRye/EIovhKN5GxEq5vMJBsNZ5SQNO58lZ8k7qU5IRdxBgKGgwZ/ItvkDE42tAHgqm9hvjjEpxwHaC67GJgdvGRhaV2wT5BuBqkZKazaUs7PX9qDv2NlIfE/P30oGJZCIoS4FPgNsBd4Rkr53pBuKBpTlqgguq8diEwH7lJIkrJBBpiQ5KWh1cuB4g38pP7XgMTXYgMB1ed9m3E7/8qv7E9CxfWwcDJtpYWsct6HU/jwnHoBKuYGP+OiMNHQzSA18Y6yQspZs7WyY16krgsZQIZiHskK4BqgWkp5btj61cADqHkkj0kpf9/NZSTQDCSipivGHlaW185VsH2VKloMD8hHI1m1lp/kbKXujBf33rcRSIQAh/TTTCI1GecxDrjZ9ha+Pe/DwjxSTmwmQXgRgF36aDn0HimmkIzLSARgbHoCf7upQP8CaeIWa1hcxzk/2o018AyFRfIE8CDwlLUghLABDwFXoYShWAjxCkpU7u/w/luBjVLK94UQY4E/AzcNwr57jxU/mXdj9IB8R5KUkIx3nqG01UvKzAWIMvWSEOAmnbSKt5FSpSzapBfKNrLFN52rCc2M/0fVRC52uSnIy6LSTB22G4b+JdLELZuP1vKjNTs7iYh2Yw0OQzHY6gMhRH6H5YXAESllKYAQ4hngWinl/SjrpSvcQMJA7LNf6Sog3xHTIhljb6Gh1UvGuHwAmox00gKNTJTVBMpexIsNB34kEhoqSWxwqjoVYeem9rvYvi+ThMNqymKlWwXi61raAdhSWse7B6u5au44LSyauOD3r+/n4fdLg8eqlZB2Yw0msRIjmQhUhB1XAou6OlkI8WXgs0AmyrqJds4dwB0AkycPk28jZqbXaJsSkqY61f+nNf8KUktfVJXx0s8z/ku5POUYE9pLkVuf4IdSgABD+tgtpyrLxAyst3nVjJQ2b4BNR2q5eUURvoDkHx+Vser2xbq9imbYUlhay72v7mPfiabgmtXuXaf0Di6xIiQiylqXUzmklC8AL3R3QSnlciHECWCp0+ks6OP+BgfTIskWzTS3+2ipUzUl/qlX4C97HfweAoadl+SnmTVhBuNLH0aIALawS2QardQE0nCYhYn/LHQFX3vvUE2wAt4TVgG/7NFCPL4AiXaDlbfr9GBNbFPicvPYxlLe2HOy001Ct3sfGmKlRUolkBt2PAk43teL9levrUEjIQOEQSZqboklJMlzPsOmT63gz77ruG/U7zmaMJf23E/jMyVEAO1SfSf44afVHJSfXD2LgrwsqswYCcC00Smhfl9CVcQXltbhNdODPX7dXkUTu5S43Pz3C7u4/uHNvB5FRHS796EjVoSkGJghhJgihHACNwCv9PWi/TVqd9AwDEjKIk02AuBpOIlH2kjPGo1/wgL+5r+W91qnkJ7kIHHqxdzs+Snt0o4QsCugpjZ+cU4qTpvB8dOqxqTSfSY46ndMWiJ5o1RlvZTw+p4TZCU7MUx1sRlCpwdrYpJVW8q5/pHNrC6qiKgNEeg+WbHAUKT/rgYuBXKEEJXAL6WUjwsh7gTeRGVqrZBS7h3svcUESdmk+JWQ0FxNvchinGGQnqT+q06cbmPWuDTqWjy04wx2EZ5nHAUg0d/CnPFpvLC9iqvmjuVkYxtfmDeB0poW6lo8NLb6gnGRxzceI8FhMGNMKgdONnHV3LH625wmpihxuXm+pIJniisiphbqupDYYiiytpZ1sb4eWD/I24k9krNJ9CkLKtlTS6Mzm3FAWqLqveXxB0hLtHOkuonFxv7gUCy7UI9HK06w93gGvoDk648XEZBw/qRMXtpxnJMNrdS1eJgzPh0gGJSva1YZXY2tvsH9u2o03WC197E6M1joupDYI1ZcWwPCsIuRACRlk+BR/bVyaOCMU7ma0hJDmp+e6GDx1By2GefgxY5PGmBTQuOqOk7A/Opm5dTPHpeG02aw74SydOblZmAz3VkOu0FLu8rsOlwdyn7pSInLzUPvHtGzTjSDQnFZPT9es7OTiOh5IbFJrGRtDQgD3v13IEjOxn5iJwCjxWmOJZ4PhCwSgPQkOwV5Wfz4tpt5bXsuF9v2MXHOYvjnvzMjM8BC+xHmB/ZSL1PJFs2ISkF2SgJ7jyshWTx1FK0ePys+KuN/vnI+/2f1DnJSEzjV2K7qV5IcEVsqcbm5/pHNBAKSBIdu/KgZWP769mH+/K9DwWC6rguJfeJaSAal+29/k5SF0erGIEA2jRxJHg1AitOGIdRwrHRTVArysijI+zLwZQgoqyLX42KV8zmEvx2kcl+1v/si85PuZX2dCkSOz0jkstljWPFRGTVNarjPJTNHs3ZbJUeqmyjIy47Y0tqSSvwdBmfpX2ZNf7O1rJ7fv36ArWFWr64LGR7EtWtr2GVtgeoA7GtlesJpbEJipKl0XiFE0CoJt06CGDZwpkHdIQy/J1iYYwg1A/4C/x4A5otD5O75G3N8BwB4/1ANAJfOUoJ1+FRzp0sfOhVyeenGj5r+psTl5vvPbOe6RzZHiAjoupDhgrZIYo02JXpL7cXgB0fGuOBLaYl2Glq9wQyuTiSmQ2IGGHYIeIPLPuycyLqI+c2HeNb5Gxwf+HHaE/lkwj1sKVXfJRZNycZpEzxfUsmMsWmkVJfg3vcOpakXstWVTIrTRovHz4pbFuhfak2/sarQxT0v7+k0LwR0XchwIq6FZNhRUQSFfwfgO/5/ApCcNSH4srJEWoOurU4kpKsOw3O/APtepi7/WnJK11J1+QPQOI/FVQ9gF8oFJvxePptymI/qp5JgN6ioP4M3INnqcvN/H32Sp2y/ZSZeLsDBWnE3u/2zAEhP6uKzNZpeUOJy88/NLl7cURWxrtN6hydxLSTDLthetlG1mwdsqBt+2uiJwZfTzcyt8AyuCBIzoL0R7GMgI5ecRddB6VqmTZvFqENOtgXCmgcYdtxjFkI9TMxKovBYfTBP/yK5Dwc+DAFO6WOxsZ8dgZkAHK1p5tyJwygLThNTlLjcrN1WwXPFlcF2PRY6rXf4EtcxkmGX/pu/BGwJgAjGOOqrjgZftmIjXVoFienKNdZ0EtLGQ3KOWj9TS3ZKAmOM06GmZlf+Ekf+YgAmZiaxeOooEuzqx6EwMIeAeaYfG1vkHJw29VppTUt//W01IwyrLmTVlopOIqLTeoc3cS0kw47chXDLK9TMuB4pVRuT6Ru+yYHiDUDIIunWtdXWCE0nIG0spJhC0lJLdoqTrxjv47VaPCaPYvroVOaLQ1xVt5KU6hJW3b6Yy2ePZpucyc7ANAB2jf0il1+1lJW3L2ZydjJHa5p1TYmm15S43Nz94u6IuhDd3iR+iGvX1rAkdyFHfC+RjcAmJA7pw73vHVhwJa1mS/jy+hbmTkjv/F7LteVrVxZJSsgisbk/oMA4rNxXAqr3f0Ty6CTWOH8NLeBZtxrXNaspyMvnvYM1jDOUSBjOZL53mXINTh2dwp6qBpYt34zHL0nUNSWaHvDoxlJ+99p+XRcSx8S1RTIs03+BrLmX48GBTxp4sZM193JKXG7+te8UAN9/Zkd0ayAxHc7Ug6cZ0saBM1W5ylpqGFPxOkKoX2IpwXt8N6mlb2ATUgkWSrAWTx1Flr2NSaIWgLyEUDrwtNGpVLpb8fjVLaHdq7sFa7qmxOXm5se3cF+YiFh1IavvuFi7seKIuBaSYRcjMZm94Epc16ymeOp3cF2zmtkLrqSwtC7Y+sTXVbv3hHSCY1zSxqv5vCk50FJHznjlNvBJ1Z0rNS2d7InKfRWQBAWrIC+LlUtD/16jCInw1NEpEb5tieourF1cmnBKXG5+vGYH1z28iQ8O10a8putC4hPt2opRZi+4EhZcGTxePHUUTruB1xfouigwMczdlWbWnySPgjO1TMxVolGSfwfnn95ARpKDjLGqCLE+dTq1l/5f9ZnAbMMcVjl6DrRUBy/p80e275bAM0UVvLi9Sru4NIBq9/7zl/ZEtHq30HUh8cuwFBIhhAH8BkgHtkopnxziLQ04BXlZrLxtcfcjcRMzQ8/TxqvHlBxoqYXTLkgZzaL/+B9YeT00n4QGJRg5OePICRMtqvcrt9jEAjj6dnD59BlPp48MH+urbxAjlxKXm2eLXazZWhUxcErXhYwMhmIeyQrgGqBaSnlu2PrVwAOoeSSPSSl/381lrkXNea9HTVccEajeWt38IiaEWSSpqrUKKaOh7gicLofMvNDayd1w2rQ8wqwOAMoLVeA+4IOWGggEwDD41IzR/P39o3h9AWyGwBuQSBnZNiXa/Hc9Ez6+KXG5Wba8MNht2kLXhYwchsIieQJ4EHjKWhBC2ICHgKtQwlAshHgFJSr3d3j/rcAsYLOU8hEhxPPA22hCri1HCiSkqefJKkaCYYfxF6i1lBwlEA2WkNSErlH6PpxU3YfZs1aJSdtpSM7uZBX9+V+H2He8gcfMtiklLjdffWQz/rAuwQA3LN+MV2d5xSWbjtTywzU7O4mI5cbSKb0jg6EYbPWBECK/w/JC4IiUshRACPEMcK2U8n6U9RKBOVnR8rP4B263w4xEM0ieNk4F2gFSRoG3BdwumPMFtZY6RvXiOqUaOdLqBr8Xjm+HF78Vup7ZUZjmakhWHYHDraLJ2ckcPNlEQV4WB4o3sOedlzlfTmUbM/H4QgkBXr/uHByP/O61/SzfWBo81mm9I5dYiZFMBCrCjiuBRd2c/wLwVyHEEuCDaCcIIe4A7gCYPHmEfCuyXFtWfARC1e3SD5nmv0OKCrLT1gBJWUpIDr8Fa25RggIgDGXF+D2m62t2p49LT7TT1OblQPEG8tfdwAx8XO90cJPnLnYyK+juEmbKsaFnwscFm47Ucu+6fRw4GeoKrdu9j2xiRUhElLUo/UDNF6Q8A3yzuwtKKZcLIU4AS51OZ0Ef9zc8sFxbbadVA8jchaGiRIiMkVhMLIAjG6D0vZCIYMDUS2HeMnjh9kjXVxipCXbafQHq977DDHxmAaWXq5IPY7RC84YPGHv+lWQkOTh9xss150/QN5lhTInLzSPvH+Uts54pHJ3WO7KJlTqSSiCsoyCTgON9vehwrSP52JzaZz7uhSe/oMQkOVxIOlgkABPmm6/lmwsC7Alw6X/DtCvUUnOYkFQUwcY/QUVRsHlkwoxLCJg/Sn5szJySx3OOX/Op8r+Tt24ZU1v3AnCiobUf/7KawaLE5eanz+/kuoc3RRURndariRWLpBiYIYSYAlQBNwA39vWiw677b19xfYj6bhBQLqmyjTD3i6HXM02tjrBITCExZ74z7TIlIrkLVbaWsIWyulyF8OQ1IANgc5K36DFAMHr2EnZsuYIFjW9RU/BfpNbVYQhlUDrM7sGlCeewq7IBf0BS4qqnuMzN4qmjgnNPsuZeHqxj0cQOq7aU8/OX9wQnZFrotF5NOEOR/rsauBTIMYPmv5RSPi6EuBN4E5WptUJKuXew9zbsyV+irAm/R80lyV8Scm2ljgVHknqePIpgSaFlkZRvUo8X3apEBMDJh6vPAAAgAElEQVQw1PubTSHZuSo0MMvvYWJDCXARjW1eAvZEACZNnER7xhRwma1YsFMYmMPnzh3H6qIKfvXKXp4udGEIuMh2mKdtv2EGfjylj/JW3eMcTpir04RjgBKXmzVby3m2uFLXhWjOylBkbS3rYn09sH6QtxNfmN2DKduoRCR3obqbC5sKnFtxE5tdpQf7PeAuA3sSuEwhGTUj8popY0IxEsP6cRFgc9I26WLY6qWpzUeCx2ylcqaOadMugXdVkGvNrL+wbVcWPzx/AgeL3yat+GXmizlskzO5SO7DiQ8hlOWy88N1POSzkWA3WHW7ThMeKkpcbpY9WojHp+tCND0jVlxbA8KwHLXbV3IXhiwKgMpi5YpqrFJxk1teUevtTYCEp65VQfrmU0pwsqdEXs/uhBO7lAg1nVRrkxfBVb9BiFnAhzS3+0jzhoSEM6q/kiGgtMmGzRBk1W3neeevAUk7Tm7y3MU2MQcIWS6b/erY69dpwkNF0bE6frRmVycR0XUhmu6IayEZcTGSaJRtDOXfWnGTcPweMMwZJVn5yjVmUVEEx3eo1OEnvxB6LX0S5C4krVYNuWpq8zLFZ6aCnqmLyPJqqdrNuPQrad3zUkTc5JKEg5z/2VsRb4JXGmxc/Bg7NzpASmw6TXhI+PNbB/nfd44Ej3VdiKanxLWQjEiLpCPW1MXwuAmAPTG0ljEZGo9DzszI95ZtVNYMgL8dfGbWVbPK3LGytprbfSQHzHbzYULil4IpgQpcmUmMz06HipD1sZVzuLStDACHCLBo8SdZXHmEj47WsWzh5I9109KtWD4eRcfquG/dfnZWhTo967qQOKGiKNLVPUDEtZBoi4TocROIXNv+NFQUQk6H+IgVvPe1hQQla0pQSFJNIWlq85Eqwy2SWnAkU+HPYUagkiNZSUxAubvaHBk8PeV/2Lg7jW+6ioMfdfqkK1gBb4hoZUXds/loLTc+ugUg2J5F3wC7p8Tl5h8flvLa7pOdirZ0XUgcUFEET1yjEmRsCep3foDEJFbqSAaEEVdH0hW5C2HJDyN/iMLXrEJEw975fbe8ChPMek7DAWPPgSYlJAl2G067QUNLK+mcUedYQpKcQ13yVGaJCiZmJqoeXkCSXZB3waUAtJ7YH/yoptpKTja2Ad3Xmxwo3sDmJ+8Kjh+2eOdANZLIbsSarvlnoYvrHt7EuigioutC4oSD65UnQQaiu7X7kbi2SDQ9oKIIdj+vnm9+CGZ9rrPgXPxdWPtN9c3m0BuqkaO3FapK+IHjRQLHLwQgIOwYZ+qVayslB19rKpONGha7Hoam45CUDa1upmarWEt2axl1RgajaKC1vopTjU4ATja0Rd3qgeINTFt3PQYBPKWPcoDVwdqT/FEpwfO6nNeiocTl5slNx3hl54mIdZ3WG2e4NsOO1eq5sEW6tQeAuBYS7drqAWUbVTAdlECUbexs/p52Eaw7CZgursNvwZpv8G0ZwHf8BQBakieS1uKCxuM0kcyF9W+AgE9U/UO9vb0BkExO9iAETBXHOZZWwKimd2muq6Ldp1q4HDeF5EDxhohiRffet7HjD6YLW7PsAbJTlAjl56Twp+vm6RthB0pcbtaWVPDs1spOxYU6rTfOqChShcMBnxKRgltUuyMdI/l46GB7D+gqGN/xHCs4b9jU49F3QarGKDZ8ALSlT1FCUn+U1qSpJHVszGzewBI8DczJ9DG6tZGjY+fR0riZ9voqAKaPSeVoTTN7C99i9utfBWTQ+kiccA6iLHI0sEVtczsAiXZD3ww7UOJyc+OjhbR3SOkFndYbl2z8sxIRi4xJAyoiEOdCoukBXQXjuzondTy8/J1gyrAEfNKGTfjwZ02DE++B34NtzCy8TaUgvdiQSAyEzewmfKaOz9tUoF0GJPVGNqJZuVouyM3kSHUztcVrsAl147Osj9QpCwCoYjR3ie+zsCmPFpebgrwsapqUkFSbjxrF1rJ6frp2V4SIaDdWHHPoLTj0eujYsA+oS8tCC4mmcxFjd+dY7VJqDwPqprTSfwW32t/EGDUtePqo6Qs4MOdm3PveYfz4ieQntak2LS9/j8pd73J7099AwIVH/0alPZcsrxuACydn8nxJJY0elbkVbn2cPlIEgD01h411U/nwrUM85DjCytsWU2NaJPUtHtp9fhLstn78BxqePLDhEH/ZcDh4rOtC4hxXIbwUNk8IARfeOODWCMS5kOgYyQCQPEr5Xcs3B5cyhCpMdIwNSx9OzmH2BVcGYxhAcLRvu6sIu+n2suNHyABjUUIyb5KaOx9oPA4C6kUGtZ9fwewFV7Jly0oAEgIqqys8Q6umKTRPvqapnUlZyf379x5GbCmt477X9rNL14WMHIJxkbB5QrYEmNfn3rc9Qqf/anqHYVPdg/0e9YMKnCNcACSPnxU6L7zDsIU5ZTFt1EQkRrA4sT5lKmPEaTKT7EzJUdlXBajUYDv+YGZWSpOaxpdqeDDMUhMrQ6umuR27uXiqcWS6t0pcbv7r2R3csLwwQkRA14XEPTufDYmINU9oAOtGOnJWIRFC2IQQfxiMzWiGCWlj1ePUSwGYLqo4IxNIyJwQOid8oJaFIxnsiYwZlY0ncyoNSZNwXbMa/9h5JAkPVyftJ6XoAb6YWMIkUUs12WTSTNsZVTU/pr0cAGeglduWTAXgD185n4K8LGqb2pk+JhWAmqbo6cPxzMpCF9c/vJkXtlfpupCRxpG3YVdYqq81T2iQRAR64NqSUvqFEAVCCCGl7HJq4WBijti9CbX/uVLKTwzxlkYWqeOAnZC7kNayLSR5T1NPCsk2ByRmqgmN0YRECFVLcqaeJE89SedcS+aCK9l6SlkavznzW3jbxx/NgZnlaRcwpukdao+XkTF6AmOoxy/s2Dwt3LRoMss/KKW+xYuUkprmdi6eNooDJ5tGlEXS1dRCHVAfIZRvgZVfUUWHhgPmf33AU32j0dMYyXbgZSHEGqDFWpRSvtDbDxRCrACuAaqllOeGrV8NPICaR/KYlPL3XV1DSrkR2CiE+CJqKJZmMBGmIWtL4EzyRJIaTtNipKm15FFKSJKjCIn1ursMWuuDnYaTsiYCYMdrPkqkhPMaPwABDdXlnGmsIw1ozphFxum95GUmcEVKGd73XueQexa3yd3MTLyal4xETjVGWiTx2INL14VoqCiCV+4MtS+SgUFJ9Y1GT4UkG6gDLg9bk0CvhQR4AngQeMpaEELYgIeAq1Bjd4uFEK+gROX+Du+/VUpppg5xI3Dbx9iD5uNSUQRH/qWev/tbApnzAGi1m/PibU4VOzm5K/oPdHI2VJWo51n56r1us8paEnTLCAE2s1Cyta6CxpMHATiTOJYM9nJwy3r+5vsFTp8PiuBHdghsf5nXUn7JqcaJwWLGw4kXcO/OVPwBGTc9uHRdiCbYR8tvWt/CGPDq9e7okZBIKf+jvz5QSvmBECK/w/JC4IiUshRACPEMcK2U8n6U9dIJIcRkoEFK2dhfe9P0gIiOwF4cQt3s2+1p6oe79qB63Zp90lFMkrPBY3YKNoXEW31Izd8S4JcgMUCCDzt2PPjKCrmo9iUkkHPifRDQvP9tHOZQLClVaqsIeFmSdIA9VQ6m7f0pNvxcgJOX5F1sYyZt3gD3vbaPuz8/d1iLyd/fO6LrQkY6JU+FRMQKrg9yXCScHmVtCSEmCSFeFEJUCyFOCSHWCiEm9eM+JgIVYceV5lp3fBP4R1cvCiHuEEJsFUJsramp6eo0TW+xKuHN/j2toy8AwOvMjGwK11WTuOSwHlimkGSd9xnacOKTBh6clJx7D8VTv0PZNatplklkNOzHTgAhwEDdQDOnFhAwM7+CGILjGQVMO70JO34MAQ7UzHiLbeWnuemxQkpc7qh/vX/tPcVf3znc5etDyfpdJ7jsj++yYX91cM1uE9y4aDKr77iY+750nhaRkUBFEexYGTq2OYZURKDnrq1/AKuA68zjr5lrV/XTPqL1De82sC+l/OVZXl8uhDgBLHU6nQV92ZwmjA6V8P7jx2HfY6S2uChr/RT5Z2u3kpRtPmZBokrLnr3gSg6wOthXa1FY7Ynr9RwkoeJEPzZs+Jg2ex6nj3waW1UR93mW8QP7WnJSkygI7CXgdyHskTPjDRHs0ILHF5rAeKBoA6d3vU7mvM9xLPEcvrNyGwKChY5DeWO2YjtZSQ5e3FFFcVmkuAng+otyue9L5w3NBjVDw+aHAMsiHbyiw+7oqZCMllKGf/t/Qgjxg37cRyWQG3Y8CTjej9fX9CdhlfB1u5czCZjt2Uv7R7+i7JO/UlXsXbVbsSwS0xqxmL2gQ/GiSaMjh+ltexECtqZfydg5n2Tyll+Dp4XMtFRqk8fzTNsVTAuc4vbmdVzR/AjCLGoXAjYueoQrEs/nK8lO7n11L22+AEKoCYwHijcw7bXrseOnreIpHkm7D5gYUeg4VEKyaks5P395T6dAuoVAzV359/n96RjQxDy7n4f9rwIiFBcZpKLD7uipkNQKIb4GmMnKLEMF3/uLYmCGEGIKUAXcgAqk9wndtHHgOXPqSDBG4ZA+TpyoIv+W33X9BrMosaOQdEVb4liS2rcDMP3rD5LVVgFbAM8ZaGsgITUL6sFnS0Bi+mpN68MQ8JlPf5rPmKnIs8al8aM1O2n1+Jg/OZPCt9ZFdBOecHorlkd1KFvRF5fVc/eLu6Oa5DoeMoI5thHW3gZIJSAXfm1IUn2j0dPK9luB64GTwAngK+ZarxFCrAY2A7OEEJVCiG9KKX3AncCbwH7gOSnl3o9z/Q6ftVQIsbyhoeHsJ2s+FlnnhuIbHTvyRsUSkuZq5es9C77U8QAcM/LIGj0enObcEW8LtDfid6q04395zsMnbXSqdGo9HXxakJfFf3wynwlNu3n3sZ+RdfIjrGGMPmxs9s8J+lgfuOHCQbtJl7jc3PXCLn7y/E5zXkhZJxExBDh1PGTkUlEEL36LoMc/4B+yVN9onNUiMVNz/11K+YX++EAp5bIu1tcD6/vjMzSDR8f4xuwo7qkImk6qR9emrjO7wjEzU9yJk5kCqjoewNMCbQ1UO8YjgG1yJr/wfoN7HU8gZIAANpz4oC3yS8RMzz5WO3+Lo9JnTVhBCHjU929skzMxzCwwrz8ytbY/alGiXWNLaR3LHi0Mxm9e3FZFZooDgbI+DENw26emkJbkiKs6GE0vqCiCJz6vYo8w5Km+0ehpZfu1wF8GYT/9inZtDQ5dxTei0niC4JAsK7OrCyE5ULyBghPPgYBzWzZzoHgDs+ecr170tEBbI5ljc0g4YeD1BXjRdhWXXXwpaacKGT9uAvmb74a2yAC1LPsQZ1jasGWRjDEasAWUS8sfkGwvP80156uWL6u2lHPPS7sJSEiwG6y6vedB+GDA3IzRtPsC2G2C603X1GMfHiM8DOINSGqaPHx1QS6Ts5O1eGigaHlIRGIg1TcaPY2RfCSEeBB4lsjK9m0Dsqt+Qnf/jUGmXgIb/9R9ZpeJe987GGaXYIOAmog472L1oke5tsbkjGHlbYs7fNP/AlQfgM13d7JIsuZcBqUPEZBKzvxSYBOSmXPm8V/jZ7F46ij+uWYNOdvXcWDsl2kZU8DdL+0Ousw8PQjCN7X5qGlu5/dPbeWdA9X4AxKbIfCZiuH1S1ZuKWfttkoS7NG9yy/vqBryrDFNDFD2EexeEzqOgVTfaPRUSKxeVveGrUkiK91jDm2RxCA9GaRlkjX3crylj4L0heIv9iT14plaNQUuMZ2CvKzON1wztTg8RgIw+9wCWA9l2UuoHn8ZzoZjXFj5NBfOncOFF0yn7I2/8oemX2AQoH3dGh7M/RNShoLuQsDiqaMi3FRA8Hmr14/9uBKv8P5XvijZV23eAG3eAN/4RB4ev2RvVQO7KhtiImtMEwME4yIWsZHqG42exEgM4O9SyucGYT/9irZIYpSeDNKim/iLIzkUa0lIj/7mJDXXhLZIIaFZ3dynXH4rU879MrTUwR+ehvZGKPuIvMJ7gu4up/SSXVMEfA4Dlbl/+ewxAHz1kc34AxKHTUUzfIEATrvBhbmZfL+bv5MhwBAiQlieKa5g5W2L+ff5k7jpsUK8vsCQZo1pYoCocZHBmy/SW86atSWlDKAyqoYdeh7J8Gf2giu5+JbfRQbxnSlmrIWQ5dERe6JynXVwbVlCQqoSBBJU63naG+HwWxGVsRLBG83TuTA3kx9+dhbJThsTM5MoLK3DF5DKcvBLPP4AAancXq66M5HbsAnOmRASOwFcvyCXJTNygp8Vbn2svG0x//WZWdqtNdI58FrnuMggzhfpLT11bf1LCPEjOsdI6gdkVxpNdziSoeksQiKEamnf2tEiMduLpJozVewJ6pteexOMngOYmVzAq/6LKfbPwHmikXuumcvabZXUNnu4+txxoa3YlHVhGRinGtsYlZqAzRDctEh14AUiLA2riLC4rL6T9RHVTacZWbg2w85n1HOzFVEsxkXC6amQWDUj3wtbk8DU/t1O/6JdW3GKM0W1ooeuXVug3Ftns0gAEtKUkKSrLC1x4dfxb1/FcdTN3e9XFkNOSgI1ze1MyFRxmhSnjae+uYg/vXmQTaV1QTHJTHIwOi0honVJ54SA6GuaEU5wZK5PTSOdf0vMFB12R0+7/04Z6I0MBDrYHqc4U8BrupASuxGSxIzoMRJ7YqQAWULSaqYKL/o2vsPvkttYj02Eqtz3Hm/gwMkmTjaoeSetXj8X5Gbi8QeYMSaVw9Wqq/Gx2hYSHTbSwj42mqWhrQ9NJzb+SYkIqK/qMVR02B3dxkiEED8Je35dh9e66YOh0QwgVlEidG+RdOXaSh0TKiABJUZtjWrYFkByNgmj8rhsvCciXpGTmsCEpt2kbf1f5otDBCTUtbRzqqkNh00E58gHpKSxzYtG0ysOvQmH3ggdG/aYKjrsjrMF228Ie/7fHV67up/3otH0DKtNCnQdI4GuXVtWfMQiIT3SIknKgoxJpLWf5HuXTQ9aDecGDvKY/DWz9z3ASufvmC8OUd3YzqnGdmaMTcNpN7CZWVnpiY5++ItqRgyuQnjx22ELsZvqG42zubZEF8+jHcccOkYSp1gWibBFikpHorq2qiG7Q2gvIQ1OV8CZelWn4kiC9InQeFz1NDJUO+GZrTtx4sUgNOfkSHUzHl+A8yZmcPPF+RSW1jHn3XTSEnsaftSMeCqK4MnPh1xaMZ7qG42zWSSyi+fRjmMOnf4bp1jikZAW6aLqSKJpkYR3cmw+FRlot67T3qjcYFZTyYxJ6he7OVRU2D5xUfDbk0RQGJjDnipl8YxNT6QgL4vvXTZdi4imd+xcHRKRYZDqG42z/cTPE0I0oqyPJPM55nHigO5Mo+kKS0i6c2tZr8uAclslpoPfC2fqori2LCGpV24tgAxzPE5DZTCbKyM1GSGgXdppNVLYIafjCBMSjabXHNkAO59Vz4dJqm80urVIpJQ2KWW6lDJNSmk3n1vHQ+YEFkJMFkK8IoRYIYT42VDtQzNEWK6t7jK2oHN1e4s5cjmqRdKkXFtBITEHRjVUBk8bX7OJgBT82ncLmbKBvyQ+iuP4VgDGpidEXrO9SWXg9KBVvmaEUr4FVl6nRiIYDii4ZdhZIhY9nUfSb5g3/2ohxJ4O61cLIQ4KIY70QBxmAq9JKW8F5g7YZjWxSdC1dTaLxBISM+AerCGJEmwP+KDpeJiQqAFXQSGpKCJt/zOUynEcDkxAAkt5n0e5l/niEGPSwiySphNwYie8fS88uVSLiaYzFUXwyp3KYgb1OExSfaMxFM7cJ4AHgaesBXPmyUOoGfCVQLEQ4hXABtzf4f23AtuBu4UQXwWeHoQ9a2KJoGvrLBaJ5fpybYbDb4HhVMfRXFsADVUw7fLQexPSlZBsfQLW/xAj4CNfGFxr2wQQDLp/OuEgSU5zvu/R96DuaOjavvZQq/yKoh41q9TEORVF8MQ1wVk7sThfpLcMupBIKT8QQuR3WF4IHJFSlgIIIZ4BrpVS3g9c0/EaZruWX5rXeh74R8dzNHFM0LV1FovEcm29eZeyOIRpgDcejzzPqkWR/pBFAirwvntNqL4EEEIF7qWwI6QPHzaOJl+oXizfAi9EqX3NX2Jm5ixVwmJPgFte1WIyUil5KiQiMTpfpLcMumurCyYCFWHHlVjDs6PzBvB/hBAPA2XRThBC3CGE2CqE2FpTU9NvG9XEAEHX1tksElNIAl5AKqEAdbMPdzclhNWgJ5lZWxVFcLo8QkQAJAYv+JdQftn/A+Af/s9SlzUv1K21pTp0siVc2dOUJeJrJ2KgV0WRjqOMNMq3wI6VoeMYnS/SW2IlTzFaDmeX6cVSyj2oufFdIqVcLoQ4ASx1Op0FfdyfJpborWurI35v5GTGCCExLZKyjZ3fZ9h5LOW7bGubyZFRBYx25jDeX68ytsreNgXLuk4mfOEBePm7sO77MO1KVY8S8KmK5aRRyr0R8KqaASvIqt1f8c3mv6IGEsBwKzrsjlgRkkogN+x4EnC8i3M1Ix3LtXV8h7rxdvWLWH0g9FzYlIUgA5390eGCZNWR5C9RN3i/RwnAhV/jwNjP84cXPYDkzme289K4Apa0b+JMxQpcGdPIUx+kPiczD0ZNV8f7X4XDG2DsuXBiB1zyU2itC7k3LAsFzBkUPtP9NTwzeDRdsGsNHFhP8GfE5hxWRYfdEStCUgzMEEJMAapQrVn6/C+smzbGKVbn3yMboOzDrm+45R8RnA8PMP9rqj6k47f9aK6tKJMc3373CFIeBNQMkdL2DGaLJr7a9BS+zXb1Ued+GcbtVNd0fRi6rt8Tis1kTOpcXd9QqQrTrBkUZ5lnrxlmlH5gxs+kEpALvzYsuvr2lEEXEiHEauBSIEcIUYkKmj8uhLgTeBOVqbVCSrm3Hz5Lt0iJR+pLzSey+xtu/hLV6deaDz+vCzdCeKwlPNjeYZLj4qmjcNqN4AyRZDNRyyYkSLMy+ZKfwoPfCn2+za5caYY9FG/xtqrrOlPA06piNyVPqHMsbI5hncWjCSM4Mtf8QhPwD+tU32gMRdbWsi7W1wPrB3k7muHIrM9B4d/UDbq7tMmezocPt0gs11YUrAmG1gyRlFMOAq+pCdQSQcCwY2RPi/z8f18Bz30dZv0b7HtRrftUG3r8Phh/vnJ3yYC6wVgs+WFc3WhGLBVF8MS/qZ9ViItU32jEimtrQNCurTgld6FKn+1JULon8+HtCeqX2+8JZXp1QcQMkbwrce++krTKD2jPmkFKgkNZIOHMWQpp4+HYe6E1XxsEAipGMvY8JSQIZZH4TTFp0ZmGccGWh0MiEiepvtGIlfTfAUEIsVQIsbyhoeHsJ2uGF7kL+/dbe0IaONPA7uzV27IWfw279JBSv08F0zsihNqj1aIewNsWskpypikxyciFJf+l1tLGqyJKzfCm7EPYszZ0HCepvtGIayHR3X81PcbmVG6H3tZ0TLmEYEB/7DnRz8ldpB7TxqtMMF9rSEjsSTDxQjXx0YqRnH89nNoNL30X1v2nrjMZjgTjIhbxk+objbgWEm2RaHpERRE0nYT2BnjyC727cSdnQ84M9dycW9KJSebNw56oXF/eNhVwB3AkQs5MOFMLx7er9i3pZi3ujpWwdYWqN9FiMnywilOtPm3CUP/3cZLqG424FhJtkWh6RHjxYXhNR0+oKAplkf3rF9Fv+Fa9iLsMPC2qYt4SEnuSEhKA0vchKz/SDfZx9qQZWvavC6VxD9P5Ir0lroVEo+kRVpqwNQ+iNxk1ZRtDg7OsivmOVGwx26WY5zVWKvcWmBaJadF4mpSQTLtc7cNCiLjL8olbXJtVPRConyd7QtzGRcKJ66wtXUei6RE9TROORv6SUMZXUIRejXKOWSUvA6r9vTcsRpKZF7pGVr76/G+8BjtXQWUJnNwF+14O7VUTm1QUwZPXmG1wbDD/lrgqOuyOuLZItGtL02M+bhaYJUKX3921+yL8nFHTwZkcZpEkqZuOVX+SlR96zzX/Dy67Sx1vfqj38RvN4PLBH0MjcyVxV3TYHXFtkWg0g0JPalWscw7/S2VsWRaJI0k9Jo9Sj1bsxKJ6n/nkLFX8mqHl4Btw+M3QsWEfUe7IuLZINJqYw55o1pG0ho4riqCiUB2/eVek1ZG/JJQNFocV0XGBqxBe/HbYQnyn+kYjroVEp/9qYg57ohIRb5hrq2xjaORqx4B97sJQ2ujXXxxRN6dhQUURPPl5aDMz7UZAqm804tq1pVukaGIOR6IacOUNs0jCg/HRrI5x56lHK024K6xZJkmjVJt6PdNk4NmxMhQXieMWKGcjroVEo4k57EmRLVIcSWfPGrM6Ere6IWWU8sdX74081yqCC9YvCCVScV6/MKQc3gC7VNPOYOr4CBQRGKZCIoSYC/wKqAPellI+P7Q70mh6iD2hs2sLug/YWzNSWt2wfaWauiiMyMmKZRvDmgNCMDi/c5WeuDgQlG+BVdcpl6ThgPlfHzGpvtEY9BiJEGKFEKJaCLGnw/rVQoiDQogjQoifneUynwP+KqX8DnDzgG1Wo+lvHB0sEnvi2d8TtEjq1TAvUDew8Ir38KC8hWGD7avgnd/q1OH+pKIIXv5eKK4lAyMq1TcaQxFsfwK4OnxBCGEDHkIJxFxgmRBirhDiPCHEug5/xgBPAzcIIf4AjBrk/Ws0Hx97opn+22pW04uzvyfJbG3f6la9uCzC4ym5C2HGZ9UY4hmfBQRccGOoCFK3WekfKopU77O6w+o4TueL9JZBFxIp5QdAfYflhcARKWWplNIDPANcK6XcLaW8psOfavPP94CfAbWD/FfQaD4+jiQIeMHT3DNrBCJjJJYrzHB0zuJyJiuhmXMNIGHaFSErZYTVNQwYJU+GeqeNkD5aPSFW0n8nAhVhx5XmWlSEEPlCiOXAU8AfujjnDiHEViHE1poaPSRIEyPYE9Rj6+mQKJyNxHwGUGcAACAASURBVAxAKCFpqVZrAW/oWhbtzZCQqmabgOpMbInHFb/o2c2uogg2/km7waJRvgV2rAodx/F8kd4SK8H2aPa97OpkKWUZcEd3F5RSLhdCnACWOp3Ogr5tT6PpJ+ymeLS6e26RGDYlJq1uaK5RVkfzKSgvhInzQ+e1N6nhXJaQNFRCm1lDlTbu7J9juW387Wqf+pt2JJseAMy4yAgsOuyOWLFIKoHcsONJwPEh2otGM3A4TPEId1P1hKQsOFOvBGTceZA6DrY9GWk5eJpMi8Q05k9XgPuYet7Wg6Lcso2h9GEdU4lk5zMq7RphdvUdeUWH3RErQlIMzBBCTBFCOIEbgFf6elHdtFETc1gWSVsvXFughKTVbc5yF+qx5gCsuBpe/b4SlPZmNTLYkQQpo9WURWu2SU+EJDzzy7DpmIrF0fdUCxTpV+6sglu0tdaBoUj/XQ1sBmYJISqFEN+UUvqAO4E3gf3Ac1LKvf3wWbpFiia2sOIaZ9whUekJSVkq/belBrwtBD2/0g8lT6j03tZ6cKaq9YxcNTPcojshseIiADPNhMpF39I3SlBxkRduJ/jvHfCP+FTfaAx6jERKuayL9fXA+kHejkYzuFhWSHtDyM3VE5Ky4Pg25XIadz5UbTdrUayhWh7VqiMhTR1nTFLnW3QlJBVFsOKzajiXPREmX6zWpYw8ZyQWNVrdAgJmoadO9e2SWHFtDQjataWJOcID7D0NtoPKwLLcVBMvUq6Vi74RqkOxOdUNL9wisUib0LWQBBtGmpXwjeaccXeZeqwogieXwtu/gX/8G7x858jJ6ProgZCI6FTfbolrIdGuLU3MER4X6W2MxCJ1TGjw1ZRLVAuVZWZaqmWRZJpCkjpOZWy1N0a/7uRPhJ7bnCFxqzeD9GUbVZNJpLqpbn96ZFTJH94AB9aFjnWqb7fEtZBoi0QTc4TXfvTGIukoJBajzDHSObPUY4JlkUxSj9lTVOpwVxaJ1VE4ZbT6tm3hPqbcW3mf6vyeeM/oKt8CL30nbEGn+p6NuBYSbZFoYg57P1gkKWMi19tOhyyOjq4t7xkVIO5KSM7UqUfDoW6U7U2h9zWfApsdkJD3CXUOxHdGlzVfxCr8HKHzRXpLXAuJtkg0MUd4gP3jCImwRYpKUraKcTRWqWPLtdVidg46sQtcH5lpw1E4Y55nCU17E2ROVs/rj8GHf1HtVZb8CL62FhBw7lfi99v5gdfCuijruEhPiWsh0Whijohg+8cQktQxYBid10+bHYYsITm5E/XrLZXQWJaGhZXyW26O+PW2qBtoe5PKCgPY9zLsf1Vlgz1zkxK+ceeFAvLxRvmWyPki9gQdF+khsdIiZUAQQiwFlk6fPn2ot6LRKMKFpLfpv6BiGdHWT5erR8u1lb9E3Qj9HpXZFfCpoLk9QVVpv/httS7CWs+31Kj2KGPPVd/Mix4JvWbFRXIXwc7Vyl3WsW39cCY81VfYVNHhCJ4v0lvi2iLRri1NzBHuzvo4FomnJTJjqqOQWMF2a+ri5XfDom+rtbZGCATUfBLLUgmOiUX15oKQ9WLN2wivn8hdqDoXv/Hf8ZW5tel/w1J90UWHvSSuhUSjiTkMWyho3RuLpNacf1F/NDL91hKSBsu1lR56T+5CWPJDGD9PHbc1wHu/C52LiLQqrPXGKkI9UzvECWxOtVy0PH7SgI9sUC48C91yv9fEtWurO7xeL5WVlbS1tQ31VjRDTGJiIpMmTcLhcAzOB9oTwePtnUVSvin0POhmWqgKFSEUI/n/7Z15eFRVtuh/iwiEEECJ0MocMTIEyMCgKKEZJNI0CAQQFFtxoDtX2teDbTcqiOKz9X3gg4vI5QGhc0FEBQTRh4h2oAFFQ1oiINDEQDAREDsMgiZAkn3/OHVSlZBKUqmqVCWs3/fxVdWuc/ZZtamcVWvYa9muLVdCHRb5sU9gu0vXhVZdoHlbyP572TlujIXDH1rXKd+HPD/bcbIpK0dd5dhn8M5vXAY01bcm1GtFUlmMJC8vj2bNmtGpUyekOl3qlHqJMYb8/Hzy8vKIjIysnYs2DLUq9XqStdUpwVI89s3d/sUc6uieeP6E9Uu6fI8ScCoSW2HYlBRbcZJrmlh95G3XVtt4ywKpqCxKZIIVQzDFFZcLqUvlVOxUX9u9Jw0gpLGm+taAeu3aqixGUlhYSEREhCqRqxwRISIionYtU9sS8USRuMY8XNNRQ65xuLOMZY1U9H22FUmxS80oaWDFBH78N7S8yRq3FUnjZk63WHll0L4fDJpuPR/+Utn3cz6xyqikvVg33F771rjEiDTV1xvqtUVSFapEFAjA98COjXiysx2sG1xFN7km11obEl3jI67YisSun3XH7y0XVXaaZZW0iYVTX7koEjfz2PSbCttetmIptgXSKNwaswPWwe72OvappUjAsrDKu/AUjwh6i0REbhKRFBFZ6zLWVET+W0SWisjkQMrnDSEhIcTGxpb+y8nJCbRIbNu2jZEjR1Y43qJFC+Li4ujSpQsDBw7k/fedtYgWL17MihUrADh06BCxsbHExcWRnZ3NggUL6NatG5Mn19n/Kt9iu588sUgqww64N64gPgJOxXDqoHXs0GehfV/LvfbDcbi2o3UzPVduL0pl12vfDzJXW/1Q/j4bPvizVcYeAKm6Sm6gWvrmpsPbD1qpvgVnLHeg9hfxGr9aJCKyHBgJnDLG9HAZHw78JxACLDPGvOxuDmPMEeARV0UCJAFrjTHvichbwCq/fAA/06RJEzIzMz0+r6ioiGuu8c1/XXFxMSEh1dsPkJCQUKo8MjMzGTNmDE2aNGHo0KEkJyeXHrdhwwZGjx7N888/D8CiRYv44IMPqh2D8OXnC0ps15anFok7mjgC7hUF2gEaNXXGNVp1s9xfdgkVUwxhEY5Wvg5FUJUiAWjVFb7ZVfF717aHcSlw6SfY+hLcPLTsTdres1FcZLnm4u6vnT0buemQOsJl5zpWPTFN9fUaf1skqcBw1wERCQFeA34BdAfuFZHuItJTRN4v96/1lVMCViteO4ex2E+yB4TCwkIeeughevbsSVxcHFu3bgUgNTWVCRMmMGrUKBITE3nsscfYuNEqsjd27FgefvhhAFJSUpgxYwYAY8aMoXfv3kRHR7NkyZLSa4SHh/Pss89y6623smvXLjZv3kzXrl0ZMGAA77zzTrXkjI2N5dlnn2XhwoUAPPfcc8ydO5dNmzYxf/58li1bxuDBg0lOTubIkSPcfffdzJs3jx9//JGHH36Yvn37EhcXx7vvvlvh5wOYM2cOffv2pVevXsyaNQuAnJwcunXrxtSpU4mOjiYxMZGCggIAvv76a+68805iYmKIj48nOzvb7TwBxXZt1ZZFIuJ0b7VyFHe81qXMvK1IwPqFXh25GjWt4DqO20nraOvx9bHwj5evjJcc+YejpW+J9Zjxt9qJqWQsL6tEqmM5KdXCrz/7jDHbRaRTueF+wNcOSwMReRMYbYx5Cct6qQ55WMokk1p0z/3z2Bk+O5LPbTdF0LvjdVWfUAUFBQXExsYCEBkZyfr163nttdcA2LdvH4cOHSIxMZHDhw8DsGvXLvbu3UvLli1588032bFjB3fffTfffvstJ06cAGDnzp1MmjQJgOXLl9OyZUsKCgro27cv48aNIyIigh9//JEePXowe/ZsCgsLiYqKIi0tjZtvvpmJEydWW/74+HjmzJlTZmzEiBEkJycTHh7On/70JwA2b97M1q1buf7663n66acZMmQIy5cv5+zZs/Tr148777zzis+3ZcsWsrKySE9PxxjD3Xffzfbt2+nQoQNZWVmsXr2apUuXcs8997Bu3Truv/9+Jk+ezPTp0xk7diyFhYWUlJS4nWfgwIFe/M95yTU1jJG4o1SRVGJJ2BZH627W6xYdnO81vR5CmzvnqE7MqPto2J1i7ZZv0AD6/9aa44uV1t6U0j4nXBkvce2VAtRKKnFuurWj36ZBQ4j/le5e9xGB8B+0xWlNgKUUbnV3sIhEAC8CcSLylEPhvAMsFJFfAu+5Oe/XwK8BOnToUNEhpTz/3lccOO6mX4OD84WXOXTyPCUGGgh0vaEZzULd7zvo3qY5s0ZFVzpnRa6tnTt38vjjjwPQtWtXOnbsWKpIhg0bRsuWlhsjISGB+fPnc+DAAbp3786ZM2c4ceIEu3btYsGCBQAsWLCA9evXA5Cbm0tWVhYRERGEhIQwbtw4wIpnREZGEhUVBcD9999fxnqpDOPaRa+abNmyhY0bNzJ37lzAssC++eabKz7fli1b2LJlC3FxcQBcuHCBrKwsOnToQGRkZKkC7t27Nzk5OZw/f55vv/2WsWPHAtbekMrmCQpF4muLpFEVigScFknT6y05igrLWiTVcWuBM4usfKpv9larorDrr/yQhtAkwoqJdEoAW091HwsHNgDG/5bBgXdxbrIUiL8fRs7z3/WuMgKhSCr6ueP2jmSMyQeSy439CDxU2UWMMUtE5AQwqlGjRr1rIqgrPxQWUWK3bTbW68oUSU2p7ObctKnTndC2bVvOnDnD5s2bGThwIKdPn+btt98mPDycZs2asW3bNj7++GN27dpFWFgYgwYNKk1xDQ0NLRMXqWnW0p49e+jWrZtH5xhjWLduHV26dCkz/vnnn5f5fMYYnnrqKX7zm9+UOS4nJ4fGjZ17JUJCQigoKHC7bu7mCSgNa5D+WxlVubZcuexIcxaxYgP5X5dTJFVkbLlSURZZWEsrqN8mzjnWfxps/ouj1lcodB1hPY5bBtd1sDoRjpzvX8ugyP7cjgwt3SviUwKhSPIAV9u2HXA8AHKUUpXlAJZba/Kyz7hcVELDaxrwn5PifOLeKs/AgQNZtWoVQ4YM4fDhw3zzzTd06dKFL7744opj+/fvz/z580lLSyM/P5/x48czfvx4AM6dO8d1111HWFgYhw4d4rPPPqvwel27duXo0aNkZ2fTuXNnVq9eXS059+7dywsvvMCyZcs8+nx33XUXr776Kq+++ioiwp49e0qthfLHzZw5k8mTJxMeHs63335b6c7z5s2b065dOzZs2MCYMWO4ePEixcXFbudp3dpd+K0WKHRYvyf2Qge3xnj1Casi2J6bDif3Ws/XTHFmKNlK48wxzy0St7JEwE+nnWXswUo7trssFl+Eb/8Jrbtbgfbb/xd8sgAyV0LETf5TJie+hIgoiL23bmyWrGMEIv13NxAlIpEi0giYBGys4pwa4cuijb07XseqR2/jj4ldWPXobX5RIgCPPfYYxcXF9OzZk4kTJ5KamlrmF7grCQkJFBUVcfPNNxMfH8/p06dJSLDcA8OHD6eoqIhevXoxc+ZMbrvttgrnCA0NZcmSJfzyl79kwIABdOzY0a1sO3bsKE3/nTZtGgsWLGDo0KEefb6ZM2dy+fJlevXqRY8ePZg5c2aFxyUmJnLffffRv39/evbsyfjx4zl//nyFx9qsXLmSBQsW0KtXL26//XZOnjxZo3n8Sm46fL3Fer5itG8CzFXFSHJ2OG1+OxaRm27dXAHeuh8uF1Q+R3UJi7BiMRdOOscaNXMG4hH4Md8qRw9w+ohlHeXs9F/A/YfjkLcbYiZVvMlS8RqpiZ+72pOLrAYGAdcD3wGzjDEpIjICmI+V/rvcGPOin65vl0iZmpWVVea9gwcPeuyWUeovPv0+DBpkPW7bduV7O16xqu+aEsvNMuQZ6+bmDd98Bsvvgi4jYMAfrrxR5qZbN2m7vIod23CVIzIBjmyDHuNg/PKay7JrEXz4lJX+u+4RayzmXmsD4NljVgmS4oswYq61sXHHK/D3/w2UAA1g6Azv1wMgexvsWwshIfDTGTj4riVTz/Hez30VISL/NMb0qeo4f2dt3etmfBOwyZ/XVpSgpFOC42Z6yXcB5vOOX///2mwFu8tvrnMXGHeVo3W0pUh8YZGAFScBK0PrxJdWmfsbY5xWkF112O6bUlQAlFg77nPTvbMa9q2DdQ9fOf7ub63uj2qR+Jyg39nuDdqPRAk63NXM8obT2Q7XUYnTdVXRdV3dOuXluNHRFdFXiuT7Q9Zjh9vg1AHAWFlaNnY/E1uObqOt8cw3vHdxfbqg4nF3a6N4Tb1WJCIySkSWnDt3LtCiKIoTdwURa4pt5dgZSdW1clzlqEnWVkWEOeI1pw5CwzBnLASsdr42xZedN/X2/aCNo2eK654Sd1RWXuXg+w6rp1wmomtzLsXn1OM6FJZFArzXp0+fqYGWRVH8hjvXlSfYiiQvwzvXkm2RnDlqubUiHC0cmlwHzdtZab/Fl6+8qXdKsMaKL1luL/u98mXpc9OtCsMlRdZcrlZdbjqseZDSfSlRiVaP+xtioCBfs7X8SL1WJNqzXblqcFcZuLrYrXqztsDR7TV3u9mKxJRA+M+seltgFUjcPB2G/5+Kb+rt+8EDG+H1cVYacM4O+O4AbP6zQ/E0dipLdxWGj2xzloUvKbb6qvgicK9USb1WJGqRKEo1OZuL5Q7yslxJo3Cr/EjJZcsaOHPU+V7xJUuJuLu5d+wPXYbD/nXw3VeOnilFznNty0QaOLLNGpS1auximOrGqnXqdYwk2Dl58iSTJk2ic+fOdO/enREjRpSWQ/GU1NRUjh/3fF+nXWyxovG2bdsSGxtLVFQUSUlJHDhwoPT9Rx99tPT1mjVr6NatG4MHDwbg3nvvpVevXsybpyUo6gw3/dxyFXkaZymPiNMqadrKahZ1TZPqz9u0lfVoSlyaTuE8t30/uK6TNda6u9PdteMVyNkODZtafUW0LHytUq8tkmB2bRljGDt2LA8++CBvvmkVk8vMzOS7777jlltu8Xi+1NRUevToQZs2ba54z5NS8a784Q9/KC28+NZbbzFkyBD27dtHq1atyuxoT0lJYdGiRQwePJiTJ0/y6aefcuzYsWpfp96Xja8L+CLOYhMWYW1IDG/t+bw9xsHuZS5KxGEljV3sPLfY8d7JvZA6yiUwb6BtX/j5n2suu1Ij6rVFEszpv1u3bqVhw4Zl+njExsaW7kz3pIT62rVrycjIYPLkycTGxlJQUECnTp2YPXs2AwYMYM2aNSxdupS+ffsSExPDuHHj+OmnnzySd+LEiSQmJvLGG28AMGjQIDIyMpg9ezY7d+4kOTmZJ598ksTERE6dOkVsbCw7duwgOzub4cOH07t3bxISEjh0yEoLnTJlCn/84x8ZPHgwf/nLXyotL5+UlMTw4cOJioriz3923iQ2b95MfHw8MTExpTvs3c2jVANfZZPZJVts68KTedv3s+IodtaVvd/EtRd94TlHNpixrBAMpVv3j38R/C1+6yH6M9ATymeQeMH+/fvp3bviWpI1KaG+cOFC5s6dS58+zk2ooaGh7Ny5E4D8/HymTrVCRTNmzCAlJaW0ynB1iY+PL1UENs8++yxpaWml1542bRojR44srWo8dOhQFi9eTFRUFJ9//jmPPfYYaWlpABw+fJiPP/6YkJCQSsvLZ2ZmsmfPHho3bkyXLl14/PHHCQ0NZerUqWzfvp3IyEhOn7aaMr344osVzuNaEFLxM7YiCa9hPbOL5ywXmTHWP7DSibv8wgqiXzznLAtzBSa4W/zWU1SRAHwwHU7uq/yYiz/Ad/udQb6f9ag85/6GnvALt40fK8XTEurucO0tsn//fmbMmMHZs2e5cOECd911l8dyeVpO58KFC3z66adMmDChdOzixYulzydMmFDqcqusvPzQoUOxrcru3btz7Ngxzpw5w8CBA0u7LrqWn69oHi2HU4vYrifXwo2eUH73f6Omzg2OFx0FL1t3h9zd1jG20sFY52mQvdap14rEpzGSwnPORj2mxHrtxeat6Oho1q5dW+F7npZQd4frr/ApU6awYcMGYmJiSE1NZVtFdaCqYM+ePWUsnqooKSnh2muvddtOuHzZeHfl5ct/5qKiIowxFZa/dzePUkvkpkPWh9bzD5+yflB5ah2Uj6tse9lZcqXgrPV4Y0zZY8Bn3gLFc+q1Iql2+m91LIfyhe/GLfPqCztkyBCefvppli5dWupy2r17Nz/99JPHJdQBmjVrVmlV2/Pnz3PjjTdy+fJlVq1aRdu2bT2Sd926dWzZsoVXXnml2uc0b96cyMhI1qxZw4QJEzDGsHfvXmJiYq44trrl5W369+/PtGnTOHr0aKlrq2XLlh7Po/iYnB2W+wksy6SmbibXfTGtu8HuT6x5Cx2KJLTFlXtnVIEEjHodbPcpPq6RJCKsX7+ejz76iM6dOxMdHc1zzz1HmzZtalT6fMqUKSQnJ5cG28vzwgsvcOuttzJs2DC6du1aLRnnzZtXmv77+uuvk5aWRqtWrTz6nKtWrSIlJYWYmBiio6PdBr+rW17eplWrVixZsoSkpCRiYmJK3XiezqP4mE4JvkkjdqVVF6sx1ZkcyxMAEHqt9/MqPsOvZeR9gYjcBDwDtDDGjHc3Vhl9+vQxGRkZZca0jLziSq2Vkb8a8GFSijXfbki5E3pNhOtvgbQXIPkTuKGH93MrlVLdMvJ+tUhEZLmInBKR/eXGh4vIv0TkaxGZXtkcxpgjxphHqhpTFCVI8HVRyksXrMe9b1vxEoAmapEEE/6OkaQCC4EV9oCIhACvAcOw2u7uFpGNWE2uXip3/sPGmFN+llFRlGDmuN1m2jg3KqprK6jwd2Or7SLSqdxwP+BrY8wRABF5ExhtjHkJGOlPeRRFqYN0SqB0h3uDECgpsVKClaAhEMH2tkCuy+s8x1iFiEiEiCwG4kTkKXdjFZz3axHJEJGM77//3ofiK4pSq7TvB+1vtaoJ3zLc6nlSQeq3EjgCkf5b0TfAbcTfGJMPJFc1VsF5S0TkBDCqUaNGFW8hVxSlbnBDT/j+oFUqJTT4Sh5d7QTCIskD2ru8bgd4XrZWUZSrhxbtrNTfc3kaHwlCAqFIdgNRIhIpIo2AScBGf1womIs22uTl5TF69GiioqLo3Lkzv/vd77h06ZLb48+ePcuiRYtKXx8/fpzx46vMgK4W7krKK0rAadHOevzugGZsBSH+Tv9dDewCuohInog8YowpAn4LfAgcBN42xnzlp+sHdc92YwxJSUmMGTOGrKwsDh8+zIULF3jmmWfcnlNekbRp08ZtqRVFqTfYiuTSeXVtBSH+ztq61834JmCTP69dF0hLSyM0NJSHHnoIsOpIzZs3j8jISCIjI/nwww+5ePEiR48e5b777mPWrFlMnz6d7OxsYmNjGTZsWGm13f3795OamsqGDRsoLi5m//79PPHEE1y6dImVK1fSuHFjNm3aRMuWLVm6dClLlizh0qVL3HzzzaxcuZKwsLAAr4aiVIKtSEBdW0FIvS6REuyura+++uqKUvLNmzenQ4cOFBUVkZ6ezqpVq8jMzGTNmjVkZGTw8ssv07lzZzIzM5kzZ84Vc+7fv5833niD9PR0nnnmGcLCwtizZw/9+/dnxQprO09SUhK7d+/myy+/pFu3bqSkpNTK51WUGhN+g1V2BdQiCULqddHGalf//f3vwU2F2hoTGwvz51d6SGUVbEWEYcOGERFhtS1NSkpi586djBkzptI5Bw8eTLNmzWjWrBktWrRg1KhRAPTs2ZO9e/cCvikpryi1Ssg10LwNnMvVGEkQohZJAImOjqZ8DbAffviB3NxcQkJCrlAyFSmd8riWXG/QoEHp6wYNGlBUZO0KnjJlCgsXLmTfvn3MmjWLwsJCbz+Kovif5o7tZmqRBB1qkUCVloO/GDp0KNOnT2fFihU88MADFBcX88QTTzBlyhTCwsL46KOPOH36NE2aNGHDhg0sX768ynLx1cHbkvKKEhBatLO2MmuMJOhQiySA2KXk16xZQ1RUFLfccguhoaH89a9/BWDAgAH86le/IjY2lnHjxtGnTx8iIiK444476NGjB08++WSNrluTkvKKEnDEcbu6cDKwcihXEPRl5H1BXSwjn5qaSkZGBgsXLgy0KFcFWkY+yMlNh7+NgJLLVjvdKe9rI6taICjKyAeaYN9HoihKNcnZAcbRebHE0XlRCRrqtSIJdtdWZdgBcUVRsCoAhzT2bedFxWfU62C7oij1BLvVtS87Lyo+46pWJO72cShXF1dDnLBe0L6fKpAgpV67tiqLkYSGhpKfn683kascYwz5+fmEhoYGWhRFqbPUa4vEGPMe8F6fPn2mln+vXbt25OXloU2vlNDQUNq1a1f1gYqiVEi9ViSV0bBhQyIjIwMthqIoSp2nXru2FEVRFP+jikRRFEXxClUkiqIoildcFSVSROR74Fig5fAB1wP/DrQQQYSuR1l0Pcqi61GWmqxHR2NMq6oOuioUSX1BRDKqU/fmakHXoyy6HmXR9SiLP9dDXVuKoiiKV6giURRFUbxCFUndYkmgBQgydD3KoutRFl2PsvhtPTRGoiiKoniFWiSKoiiKV6giURRFUbxCFYmiKIriFapI6jAicpOIpIjI2srGrjZEpLuIvC0i/yUi4wMtT6ARkQQRWSwiy0Tk00DLE2hEZJCI7HCsyaBAyxNoRKSbYy3Wish/1GQOVSQBQkSWi8gpEdlfbny4iPxLRL4WkemVzWGMOWKMeaSqsbqEL9YF+AXwqjHmP4AH/CZsLeCj78kOY0wy8D7w3/6U19/46PthgAtAKJDnL1lrAx99Pw46vh/3ADXasKhZWwFCRAZifZlXGGN6OMZCgMPAMKwv+G7gXiAEeKncFA8bY045zltrjCnzy7uisbqAL9bF8TgL+Am43RhzRy2I7hd8/D15G3jUGPNDLYnvc3z0/fi3MaZERH4G/F9jzOTakt/X+Or7ISJ3A9OBhcaYNzyV46rtRxJojDHbRaRTueF+wNfGmCMAIvImMNoY8xIwsnYlDAw+XJdpjj+od/wla23gq/UQkQ7AubqsRMDnfzdngMb+kLO28NV6GGM2AhtF5P8DHisSdW0FF22BXJfXeY6xChGRCBFZDMSJyFPusCmj8gAAA7FJREFUxuoBnq5LJxFZAqwA5vhZtkDg0Xo4eAT4m98kCiyefj+SROT/ASuBhX6WLRB4uh6DRGSBY0021eSCapEEF1LBmFvfozEmH0iuaqwe4Om65AC/9ps0gcej9QAwxszykyzBgKffj3eo45ZqFXi6HtuAbd5cUC2S4CIPaO/yuh1wPECyBBO6LmXR9SiLrkdZan09VJEEF7uBKBGJFJFGwCRgY4BlCgZ0Xcqi61EWXY+y1Pp6qCIJECKyGtgFdBGRPBF5xBhTBPwW+BA4CLxtjPkqkHLWNrouZdH1KIuuR1mCZT00/VdRFEXxCrVIFEVRFK9QRaIoiqJ4hSoSRVEUxStUkSiKoiheoYpEURRF8QpVJIqiKIpXqCJRlCoQkWIRyRSRL0XkCxG5vRau+bCI7BORvSKyX0RGO8Zni8id/r6+oniC7iNRlCoQkQvGmHDH87uAp40xP/fj9doB/wDijTHnRCQcaGWMOeqvayqKN6hFoiie0Ryr/LhdNfV9+w0RWSgiUxzPc0TkryKyS0QyRCReRD4UkWwRSXY5f7uIrBeRA2J1qWsAtAbOY/WZwBhzwVYiIpIqIuNFpI/DSsp0WC7G8X5nEdksIv8Uqwtg11pcG+UqRav/KkrVNBGRTKyOejcCQ6p5Xq4xpr+IzANSgTscc3wFLHYc0w/oDhwDNgNJwHrgO+CoiPwdeMcY857rxMaYDCAWQETmOM4FWAIkG2OyRORWYJEH8ipKjVBFoihVU2CMsW/a/YEVItKjGufZhfL2AeHGmPPAeREpFJFrHe+luzQgWg0MMMasFZHhQF9gKDBPRHobY54rfwERuQeIBxIdLrDbgTUipZXE63TjJqVuoIpEUTzAGLNLRK4HWgFFlHUPh5Y7/KLjscTluf3a/tsrH6Q0jusYIB1IF5GPsJpSPed6oIhEA88DA40xxQ632Flb6SlKbaExEkXxAEfMIQTIx3JHdReRxiLSAst68JR+jnLfDYCJwE4RaSMi8S7HxDqu5SpHC+BN4AFjzPcAjja6R0VkguMYEZGYGsikKB6hFomiVI0dIwGr+9yDxphiIFdE3gb2AlnAnhrMvQt4GegJbMeKj7QH5opIG6AQ+J4ru16OAToCS203lsMSmQz8l4jMABpiKZsvayCXolQbTf9VlAAhIoOAPxljRgZaFkXxBnVtKYqiKF6hFomiKIriFWqRKIqiKF6hikRRFEXxClUkiqIoileoIlEURVG8QhWJoiiK4hWqSBRFURSv+B/2FbckOT01+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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": 29,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAE8AAAAYCAYAAAC7v6DJAAAABHNCSVQICAgIfAhkiAAAA71JREFUWIXt2FuIlVUUB/Df2NUs7GZU9PA9mFR0JYuifMjKCelGPfXShSKkGyFUKEjSQ0YpQkFBNywIFBQLK2oo1JzAzMmHoXypGKMay7IbZaXj9LD2YY7f+b7vTOfMjEzNHz4Oe9322muvvfbahwn8b7AAH+NX7MI6nH1QPRpHeBd3iICdg7XYieMPplOjjQyDWDHCdo/GAK6rkHkV32NKi3NcKHy/s0roBNwldvNz7MEv6E6Kk1qcnNEL3inJ7mUl/JnYj/ltzrMW/WKzCjEvOfItXsMSvIyfE301Olqc/DCcIRY7kliFbTikhN8l/J/c5jwXixgsLBOYLdI/n2En46ukfHObTowknhLZML2EP0Nk3fMjNN927FC+UaVYKIL3TI5+Pd4Xi/hLZO1G3JOTyzQe23pahpX4AX9iK66t8GcZvsNZFTJPJPtXFvC6Eu+mHL0j+TOY9OvxaKJ3VsxZiIeS4vI62t2J1i9293G8iC2inahHpjx460VB35zsvyICOIArCnx5WvPAERuwT/FFcV6y/5kDM2lZ8qkoW69KvKVN5j0Ah6JXY9R7RLadVKBzYm6cKQ/eoNjVenQm+ts5+rOix5stykntyxfyKSJwvUULSliR5rg9jWuna5Xiy3Fq4m+psNmApUnprRy9B7/juGHYyJQHr09xHdkhjnE9Bku+xTm5GYneVeHTaaKj6MN9Sf4dHF6hs0f0lcPCA8nodo2N6HxDt/Ny3IhpJXYy5cF7vUSnWxytVnCpoSyqwhJDG/Ahjmoi/43I6Ka4Nxn9VByNItwqatVAkt0vatjMnFym+sIowobEbwXnJ903msjVEmBQtFLNsBu/NRN6MBnsVVzT8jgWc/GCCOSPOb3M2Abv1KTbXSFzi9js/iT7XBObk5L8F1VCjyRj2zQW/uHgJY09YWZsg9chbvBdJfy5+FskxzRRlvaqzr4zkz9raoT8rbJI9Dc9oj/KF+x6XCNu4jxqGfdHhe5oYxAfiM3PN9GXi9fS15gjArxIrCXf29XjkvS7vkaoX/xteEwcu03issijz1CmrBT9WHeid2AWLhLBf6/CkbHAGpH9neKtTvR3b4o3+9XiyBLB3IobxBo2FdibI2JTWEcXK28Fat+GOvl54sH8pciy3eKoP4xjcrYzY3tsiZZjJz5K4+lp/BPOLZCvNcGbC3hTRZtS1hn8J7FABOSCNu3cn+zMatujcYQjRbO9rg0bk0UvuzrP+Nf/EIwz7MMnOELUtL0t2DhdPEOfFH9vTWACExi/+AcITBnNL+PAYQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$2 \\sin^{2}{\\left (x \\right )}$$"
      ],
      "text/plain": [
       "     2   \n",
       "2⋅sin (x)"
      ]
     },
     "execution_count": 29,
     "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}