{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Numerical Methods using Python (scipy)\n", "======================================\n", "\n", "Overview\n", "--------\n", "\n", "The core Python language (including the standard libraries) provide enough functionality to carry out computational research tasks. However, there are dedicated (third-party) Python libraries that provide extended functionality which\n", "\n", "- provide numerical tools for frequently occurring tasks\n", "\n", "- which are convenient to use\n", "\n", "- and are more efficient in terms of CPU time and memory requirements than using the code Python functionality alone.\n", "\n", "We list three such modules in particular:\n", "\n", "- The `numpy` module provides a data type specialised for “number crunching” of vectors and matrices (this is the `array` type provided by “`numpy`” as introduced in [14-numpy.ipynb](14-numpy.ipynb)), and linear algebra tools.\n", "\n", "- The `matplotlib` package (also knows as `pylab`) provides plotting and visualisation capabilities (see [15-visualising-data.ipynb](15-visualising-data.ipynb)) and the\n", "\n", "- `scipy` package (SCIentific PYthon) which provides a multitude of numerical algorithms and which is introduced in this chapter.\n", "\n", "Many of the numerical algorithms available through `scipy` and `numpy` are provided by established compiled libraries which are often written in Fortran or C. They will thus execute much faster than pure Python code (which is interpreted). As a rule of thumb, we expect compiled code to be two orders of magnitude faster than pure Python code.\n", "\n", "You can use the help function for each numerical method to find out more about the source of the implementation.\n", "\n", "SciPy\n", "-----\n", "\n", "`Scipy` is built on `numpy`. All functionality from `numpy` seems to be available in `scipy` as well. For example, instead of" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy\n", "x = numpy.arange(0, 10, 0.1)\n", "y = numpy.sin(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "we can therefor also use" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "import scipy as s\n", "x = s.arange(0, 10, 0.1)\n", "y = s.sin(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we need to import `scipy`:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import scipy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `scipy` package provides information about its own structure when we use the help command:\n", "\n", "```python\n", "help(scipy)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The output is very long, so we're showing just a part of it here:\n", "\n", " stats --- Statistical Functions [*]\n", " sparse --- Sparse matrix [*]\n", " lib --- Python wrappers to external libraries [*]\n", " linalg --- Linear algebra routines [*]\n", " signal --- Signal Processing Tools [*]\n", " misc --- Various utilities that don't have another home.\n", " interpolate --- Interpolation Tools [*]\n", " optimize --- Optimization Tools [*]\n", " cluster --- Vector Quantization / Kmeans [*]\n", " fftpack --- Discrete Fourier Transform algorithms [*]\n", " io --- Data input and output [*]\n", " integrate --- Integration routines [*]\n", " lib.lapack --- Wrappers to LAPACK library [*]\n", " special --- Special Functions [*]\n", " lib.blas --- Wrappers to BLAS library [*]\n", " [*] - using a package requires explicit import (see pkgload)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we are looking for an algorithm to integrate a function, we might explore the `integrate` package:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```python\n", "import scipy.integrate\n", "\n", "scipy.integrate?\n", "```" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "produces:\n", "\n", "```text\n", "=============================================\n", "Integration and ODEs (:mod:`scipy.integrate`)\n", "=============================================\n", "\n", ".. currentmodule:: scipy.integrate\n", "\n", "Integrating functions, given function object\n", "============================================\n", "\n", ".. autosummary::\n", " :toctree: generated/\n", "\n", " quad -- General purpose integration\n", " dblquad -- General purpose double integration\n", " tplquad -- General purpose triple integration\n", " nquad -- General purpose n-dimensional integration\n", " fixed_quad -- Integrate func(x) using Gaussian quadrature of order n\n", " quadrature -- Integrate with given tolerance using Gaussian quadrature\n", " romberg -- Integrate func using Romberg integration\n", " quad_explain -- Print information for use of quad\n", " newton_cotes -- Weights and error coefficient for Newton-Cotes integration\n", " IntegrationWarning -- Warning on issues during integration\n", "\n", "Integrating functions, given fixed samples\n", "==========================================\n", "\n", ".. autosummary::\n", " :toctree: generated/\n", "\n", " trapz -- Use trapezoidal rule to compute integral.\n", " cumtrapz -- Use trapezoidal rule to cumulatively compute integral.\n", " simps -- Use Simpson's rule to compute integral from samples.\n", " romb -- Use Romberg Integration to compute integral from\n", " -- (2**k + 1) evenly-spaced samples.\n", "\n", ".. seealso::\n", "\n", " :mod:`scipy.special` for orthogonal polynomials (special) for Gaussian\n", " quadrature roots and weights for other weighting factors and regions.\n", "\n", "Integrators of ODE systems\n", "==========================\n", "\n", ".. autosummary::\n", " :toctree: generated/\n", "\n", " odeint -- General integration of ordinary differential equations.\n", " ode -- Integrate ODE using VODE and ZVODE routines.\n", " complex_ode -- Convert a complex-valued ODE to real-valued and integrate.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following sections show examples which demonstrate how to employ the algorithms provided by `scipy`.\n", "\n", "Numerical integration\n", "---------------------\n", "\n", "Scientific Python provides a number of integration routines. A general purpose tool to solve integrals *I* of the kind\n", "\n", "$$I=\\int_a^b f(x) \\mathrm{d} x$$\n", "\n", "is provided by the `quad()` function of the `scipy.integrate` module.\n", "\n", "It takes as input arguments the function *f*(*x*) to be integrated (the “integrand”), and the lower and upper limits *a* and *b*. It returns two values (in a tuple): the first one is the computed results and the second one is an estimation of the numerical error of that result.\n", "\n", "Here is an example: which produces this output:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The numerical result is 17.864264 (+-1.55113e-11)\n" ] } ], "source": [ "from math import cos, exp, pi\n", "from scipy.integrate import quad\n", "\n", "# function we want to integrate\n", "def f(x):\n", " return exp(cos(-2 * x * pi)) + 3.2\n", "\n", "# call quad to integrate f from -2 to 2\n", "res, err = quad(f, -2, 2)\n", "\n", "print(\"The numerical result is {:f} (+-{:g})\"\n", " .format(res, err))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `quad()` takes optional parameters `epsabs` and `epsrel` to increase or decrease the accuracy of its computation. (Use `help(quad)` to learn more.) The default values are `epsabs=1.5e-8` and `epsrel=1.5e-8`. For the next exercise, the default values are sufficient.\n", "\n", "### Exercise: integrate a function\n", "\n", "1. Using scipy’s `quad` function, write a program that solves the following integral numerically: $I = \\int\n", "_0^1\\cos(2\\pi x) dx$.\n", "\n", "2. Find the analytical integral and compare it with the numerical solution.\n", "\n", "3. Why is it important to have an estimate of the accuracy (or the error) of the numerical integral?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise: plot before you integrate\n", "\n", "It is good practice to plot the integrand function to check whether it is “well behaved” before you attempt to integrate. Singularities (i.e. $x$ values where the $f(x)$ tends towards minus or plus infinity) or other irregular behaviour (such as $f(x)=\\sin(\\frac{1}{x}$) close to $x = 0$ are difficult to handle numerically.\n", "\n", "1. Write a function with name `plotquad` which takes the same arguments as the quad command (*i.e.* $f$, $a$ and $b$) and which \n", "- (i) creates a plot of the integrand $f(x)$ and \n", "- (ii) computes the integral numerically using the `quad` function. The return values should be as for the `quad` function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solving ordinary differential equations\n", "---------------------------------------\n", "\n", "To solve an ordinary differential equation of the type\n", "$$\\frac{\\mathrm{d}y}{\\mathrm{d}t}(t) = f(y,t)$$\n", "\n", "with a given $y(t_0)=y_0$, we can use `scipy`’s `odeint` function. Here is a (self explaining) example program (`useodeint.py`) to find \n", "\n", "$$y(t) \\quad \\mathrm{for}\\quad t\\in[0,2]$$\n", " given this differential equation:\n", "$$\\frac{\\mathrm{d}y}{\\mathrm{d}t}(t) = -2yt \\quad \\mathrm{with} \\quad y(0)=1.$$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEPCAYAAABY9lNGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHelJREFUeJzt3XeU1OW9x/H3l45UpQoICYIxiEgQ6cgCigurNAugYEdi\nLLmaEDXXe1k1iSX3aDQmekGMEEEsCIIoVRYVLCgiijQbCHoh1qOoEZbn/vEMui6zy+7O/PaZ8nmd\nM2en/Jz5OGeOX59uzjlERESKqxI6gIiIpCYVCBERiUsFQkRE4lKBEBGRuFQgREQkLhUIERGJK9IC\nYWZTzWynma0r5Zq7zGyLma01s85R5hERkbKLugXxD+CUkl40s8HAkc659sAE4N6I84iISBlFWiCc\nc88Dn5VyyTBgeuzal4AGZtYsykwiIlI2occgWgIfFHm8I/aciIgEFrpAiIhIiqoW+PN3AEcUedwq\n9twBzEybRomIVIBzziryz1VGC8Jit3jmAecCmFkP4HPn3M6S3sg5d9BbYaFjyxbHo486rrvOceKJ\njjp1HF26OCZOdCxf7vjuu4O/T6bfJk2aFDxDJt30feq7TNVbIiJtQZjZTCAHaGRm24BJQA3AOecm\nO+eeMrMhZvY2sBu4INHPrFIF2rXztzPO8M/9+9/wyiuweDFMnAhbtsCgQXDWWTBkCBxySKKfKiKS\neSItEM65s8twzeVRZgCoWRN69/a3G26AnTth3jy49164+GLIy4Pzz4eBA32BERGRLB2kbtYMxo+H\npUth0ybo2RN+9zs48kj4wx/gww9DJ4xeTk5O6AgZRd9n8ui7TB2WaB9VZTEzF3XWNWtgyhR4+GHf\n9XTVVXD88ZF+pIhIpMwMl8KD1GmjSxe45x545x3o3BlGjoR+/WDZMkiTOioikjRqQZRi716YNQtu\nusl3S+XnQ//+YBWqxSIilS+RFoQKRBnsLxQ33gjNm/uB7v79g0QRESkXFYhKsncvPPQQTJoExx4L\n//M/0L590EgiIqXSGEQlqVYNxo2DDRv8lNmePeE3v4HPPw+dTEQk+VQgKqBmTT8tdv16+Oor+NnP\n/OB2YWHoZCIiyaMupiRYtw4uu8x3QU2ZAh07hk4kIuKpiymwTp1gxQq/Grt/f7j+evj229CpREQS\nowKRJFWqwIQJ8PrrfoziuOPg2WdDpxIRqTh1MUVk7lzf7TRmDPzxj37cQkSksqmLKQUNH+5bE++9\nByecAG+8ETqRiEj5qEBEqHFjeOwxuPpqGDAAbr8d9u0LnUpEpGzUxVRJ3n0Xzj3XdzXNmOFXZIuI\nRE1dTGmgbVs/06lvX79DbEFB6EQiIqVTCyKAxYvhvPPgiivg2mt1SJGIREd7MaWh7dth9Gho0ACm\nT4dGjUInEpFMpC6mNNSqFSxfDh06+C6n114LnUhE5MfUgkgBjz0Gl14Kd98No0aFTiMimURdTBlg\n7Vq/duKcc/wBRRqXEJFkUIHIEP/6F5xxhh+XePBBqF8/dCIRSXcag8gQTZrAkiXQsqU/a+K990In\nEpFspgKRYmrU8GdLTJjgDyVavTp0IhHJVupiSmFPPAEXXwz33QfDhoVOIyLpKJEupmrJDiPJM2wY\ntGjh/27dCldeGTqRiGQTtSDSwPvvw5AhMGiQ3/BPM5xEpKw0iykLfPaZnwbbogVMm+bHKkREDkaz\nmLLAoYfCwoXw9dcwdCjs3h06kYhkOhWINFK7Nsye7VsRAwfCJ5+ETiQimUwFIs1UqwZTp0K/fn7r\n8A8+CJ1IRDKVZjGlITO49Va/sK5PH1i6FNq3D51KRDKNCkQa++1voWFD6N/fnzHRoUPoRCKSSVQg\n0tzFF0OtWnDSSfD003DccaETiUimUIHIAGPH+rOuBw2CBQuga9fQiUQkE6hAZIgzz/RrI4YMgblz\noVev0IlEJN1FPovJzHLNbKOZbTaza+K8Xt/M5pnZWjN7w8zOjzpTpho2zB9fOmwYrFgROo2IpLtI\nV1KbWRVgMzAQ+BBYDYx2zm0scs11QH3n3HVm1hjYBDRzzu0t9l5ZvZK6PJ55xp93/eCDvttJRLJX\nKq+k7gZscc5tdc7tAWYBxfcldUC92P16wCfFi4OUz4AB8Pjjfmxi6dLQaUQkXUVdIFoCRZdybY89\nV9TdQAcz+xB4Hfh1xJmyQp8+ftX1mDG+RSEiUl6pMEh9CvCac26AmR0JLDGzTs65r4pfmJ+f//39\nnJwccnJyKi1kOurbFx591A9gP/aYX30tIpmtoKCAgoKCpLxX1GMQPYB851xu7PG1gHPO3VrkmieB\nm51zK2OPlwHXOOdeKfZeGoOooGXL/JjEnDm+ZSEi2SOVxyBWA+3MrI2Z1QBGA/OKXbMVOAnAzJoB\nRwHvRpwrqwwcCDNmwMiR8MILodOISLqI/DwIM8sF7sQXo6nOuVvMbAK+JTHZzA4HHgAOj/0jNzvn\nHorzPmpBJOjpp+G882D+fOjePXQaEakMOjBIyuzJJ+HCC32xOP740GlEJGqp3MUkKebUU2HyZMjL\ng7feCp1GRFJZKsxikko2fDh8+aVfRPfss9C2behEIpKKVCCy1LhxvkicdBI89xy0LL46RUSyngpE\nFvvVr+CLL+Dkk31LonHj0IlEJJVokFq49lq/JceyZdCgQeg0IpJMmsUkCXEOLrsM3nwTFi6EQw4J\nnUhEkkUFQhK2bx+cey58+qk/T6JGjdCJRCQZVCAkKfbs+eHgoYcegqpVQycSkURpHYQkRfXqMGuW\nb0VMmOC7nkQke6lAyI/UquW7mN54A/7zP0OnEZGQVCDkAHXrwoIF/tChO+8MnUZEQtE6CImrcWNY\ntMhvD960qT94SESyiwqElKhNG7+p38CB0KiRzrcWyTbqYpJSdezojy4dOxZWrw6dRkQqkwqEHFSf\nPnDffTB0KGzeHDqNiFQWdTFJmQwdCh9/DKecAitXQosWoROJSNRUIKTMLrwQdu6E3Fy/uV/DhqET\niUiUtJJaysU5uOoqWLPGz3KqXTt0IhEpjbbakEq1bx+cc47fmuPhh7Ulh0gq01YbUqmqVIEHHoBP\nPoGrr9aWHCKZSgVCKqRmTZgzx58hcfvtodOISBQ0SC0V1rChX0jXq5c/snT06NCJRCSZVCAkIUcc\nAU895VdbN28OOTmhE4lIsqiLSRJ27LF+m/BRo2D9+tBpRCRZVCAkKQYMgDvugCFDYMeO0GlEJBnU\nxSRJc/bZ8MEHvkg8+yw0aBA6kYgkQusgJKmcgyuugI0b/diEzrYWCUsL5SSlFBbC6af7g4emT/fr\nJkQkDC2Uk5RStSrMnAnvvKNjS0XSmQqEROKQQ2D+fH+WxN//HjqNiFSEBqklMo0bw8KF/jyJFi1g\n+PDQiUSkPFQgJFJt28K8eTB4MDRrBj17hk4kImWlLiaJXNeuMG0ajBgBW7aETiMiZaUCIZViyBC4\n6Sbfkti1K3QaESkLFQipNOPHw5gxcNppsHt36DQicjCRFwgzyzWzjWa22cyuKeGaHDN7zczeNLPl\nUWeScG68EY4+2heKvXtDpxGR0kS6UM7MqgCbgYHAh8BqYLRzbmORaxoAq4BBzrkdZtbYOfdxnPfS\nQrkM8d13kJcH7dvD3/4GVqElPCJSFqm8UK4bsMU5t9U5tweYBQwrds3ZwGzn3A6AeMVBMkuNGn59\nxMqVcNttodOISEmiLhAtgQ+KPN4ee66oo4DDzGy5ma02s3ERZ5IUUL++36vp73+HGTNCpxGReFJh\nHUQ1oAswAKgDvGBmLzjn3g4bS6LWsiUsWOC3Cj/8cP9XRFJH1AViB9C6yONWseeK2g587Jz7FvjW\nzJ4FjgMOKBD5+fnf38/JySFHx5elvY4d4eGH/WFDy5b5w4dEpOIKCgooKChIyntFPUhdFdiEH6T+\nCHgZGOOc21DkmqOBvwK5QE3gJWCUc+6tYu+lQeoMNnMmXHstrFoFrVqFTiOSORIZpI60BeGcKzSz\ny4HF+PGOqc65DWY2wb/sJjvnNprZImAdUAhMLl4cJPPtP2woL0+HDYmkCp0HISnDObj8cti0SYcN\niSSLDgySjFFYCCNH+hbEtGlaIyGSqFReByFSLlWrwkMPwebN8F//FTqNSHZLhWmuIj+y/7Chnj2h\ndWu45JLQiUSykwqEpKQmTfxhQ337+vUSeXmhE4lkH3UxScpq1w7mzIHzz4dXXgmdRiT7qEBISuvR\nA+67D4YOhXffDZ1GJLuoi0lS3rBhsH27P2xo1Spo1Ch0IpHsoGmukjauuQaefx6WLoXatUOnEUkP\nWgchWWHfPjjnHH+exCOP+CmxIlI6rYOQrFClCjzwAHzyCVx9tV95LSLRUYGQtFKzJsyd67uZ7rgj\ndBqRzKZBakk7DRvC009Dr15+59ezzgqdSCQzqUBIWmrdGp58EgYN8ocN9e0bOpFI5il1kNrMagGn\nAn2BFsA3wJvAAufc+kpJ+EMWDVLLAZYsgbFjoaAAfv7z0GlEUk8ks5jM7AZ8cSgAXgV2AbXwZ0j3\nj93/jXNuXUU+uNxBVSCkBNOmQX6+XyNx+OGh04iklqgKRJ5zbkEpH9oUaO2cq5RNEFQgpDQ33eQH\nrwsKoF690GlEUkek6yDM7Ezn3KMHey5qKhBSGuf8rq/bt8O8eVC9euhEIqkh6gKxxjnX5WDPRU0F\nQg5m716/Z1Pz5jB1qg4bEoHoupgGA0OAs4CHi7xUH+jgnOtWkQ+sKBUIKYuvvoKBA2HAALj55tBp\nRMJLpECUNs31Q/zg9NDY3/2+BK6qyIeJRK1uXViwwE97bdLEr7gWkYopSxdTdefcnkrKU1oOtSCk\nzLZtgz594I9/hHHjQqcRCSeSvZjMbL6ZnVbCa23N7EYzu7AiHyoStdat/Yl0Eyf6FoWIlF9pYxDN\ngauBkcBnwL+A2sBPgLeBu51zT1ROTLUgpGJefBFOOw2eeMJvzSGSbaKexXQl8Bx+Ydw3wGbn3NcV\n+bBEqEBIRS1cCOedB888A8ccEzqNSOWKervvpsCj+IHp5vgiIZI2cnP9zq+5ubB1a+g0IumjTAcG\nmZkBg4ALgK7AI8BU59w70cb7UQa1ICQhd90Ff/ubP5WuSZPQaUQqR+QHBsX+y/x/sdte4FDgMTO7\nrSIfKhLClVfCmWfCkCHw5Zeh04ikvrKMQfwaOBf4GLgPmOuc22NmVYAtzrkjo4+pFoQkh3MwYQK8\n957fLrxmzdCJRKIV9SD1DcD9zrkDem/N7OfOuQ0V+eDyUoGQZCks9IcMVa0KDz2ks60ls0VaIFKF\nCoQk07ff+umvrVvDlCn+vGuRTBT5GIRIpqlVy28PvnGj345D/+8hciAVCMlader4VdYrVsCkSaHT\niKQenUktWa1hQ1i8GE480R80NHFi6EQiqUMFQrJekyawdKnfAbZePfjlL0MnEkkNKhAiQMuWvkj0\n6+e3DB87NnQikfBUIERi2raFRYv8gUN16sCIEaETiYQV+SC1meWa2UYz22xm15Ry3QlmtsfMRkad\nSaQkHTr4gesJE/zYhEg2i7RAxFZb3w2cAhwDjDGzo0u47hZgUZR5RMqiSxeYM8d3M61YETqNSDhR\ntyC64bfj2Bo7lW4WMCzOdVcAjwG7Is4jUia9e8OsWX7vpueeC51GJIyoC0RL4IMij7fHnvuembUA\nhjvn7gEqtNpPJAoDBsDMmXD66bByZeg0IpUvFRbK/QUoOjahIiEp46ST4J//9APWL74YOo1I5Yp6\nFtMOoHWRx61izxXVFZgVO3OiMTDYzPY45+YVf7P8/Pzv7+fk5JCTk5PsvCIHOOUUmDYNhg71O8B2\n6xY6kUjJCgoKKCgoSMp7RbpZn5lVBTYBA4GPgJeBMSXtAGtm/wDmO+cej/OaNuuToJ58Ei66yM9y\n6to1dBqRsknZzfqcc4XA5cBiYD0wyzm3wcwmmNkl8f6RKPOIJOLUU2HyZMjLgzVrQqcRiZ62+xYp\npzlz4NJL/aK6444LnUakdIm0ILSSWqScRozwhw6dcopfTNepU+hEItFQgRCpgDPOgH37YNAgeOop\nv7hOJNOoQIhU0FlnQfXqMHgwzJsH3buHTiSSXCoQIgkYMQJq1PDHlz7+OPTpEzqRSPKkwkI5kbSW\nlwczZvhisXx56DQiyaMCIZIEJ58Mjz4Ko0b52U0imUAFQiRJcnJg7lwYN84vqhNJdyoQIknUq5df\naX3RRfDII6HTiCRGg9QiSXbCCbBkiZ/d9MUXMH586EQiFaMCIRKBTp38YUODBsGnn8I1JZ6lKJK6\ntNWGSIR27PBF4tRT4ZZbwLSZvVSyRLbaUIEQidgnn8CQIb5Vce+9ULVq6ESSTVJ2N1cRgUaNYOlS\neO89GDMG/v3v0IlEykYFQqQS1KvnZzft3etXXX/5ZehEIgenAiFSSWrW9FNff/pT6NcPPvoodCKR\n0qlAiFSiatX8OMTIkX7NxMaNoROJlEzTXEUqmRlcfz20auVXX8+eDb17h04lciC1IEQCOf98mD7d\nb/L3+AGnsIuEpxaESECDBvnN/U47za+ZuOKK0IlEfqB1ECIp4P33/dYcQ4bAbbdprYQkjxbKiWSA\nTz+FM8+EQw6BmTP91FiRRGmhnEgGOOwwWLgQWrTwJ9Nt3Ro6kWQ7FQiRFFK9up8Ge8EF0LMnvPhi\n6ESSzVQgRFKMGfzHf8CUKTB0qO9uEglBYxAiKeyNN3yRGDcO8vOhiv6XTspJg9QiGWzXLr9WomlT\nmDYN6tcPnUjSiQapRTJY06awfDk0bw7du2t7Dqk8KhAiaaBGDbjnHpg4EU48EebODZ1IsoG6mETS\nzMsvwxlnwLnnwg03aFGdlE5jECJZZtcuOOssv6huxgw49NDQiSRVaQxCJMs0bQpLlsDPfgZdu8Kr\nr4ZOJJlIBUIkTVWvDnfcATffDLm5cPfdoEa2JJO6mEQywNtvw6hR8JOfwNSp0LBh6ESSKtTFJJLl\n2rWDVav8Pk5dusDq1aETSSZQgRDJEDVrwl//Cn/+M+TlwV/+oi4nSYy6mEQy0LvvwujR0KwZ3Hef\n/yvZKaW7mMws18w2mtlmM7smzutnm9nrsdvzZnZs1JlEMl3btvD889CxI3TuDPPnh04k6SjSFoSZ\nVQE2AwOBD4HVwGjn3MYi1/QANjjnvjCzXCDfOdcjznupBSFSAc895xfVnXwy3H471K0bOpFUplRu\nQXQDtjjntjrn9gCzgGFFL3DOveic+yL28EWgZcSZRLJK377w+uuwZ49vTbzwQuhEki6iLhAtgQ+K\nPN5O6QXgYuDpSBOJZKH69eEf/4Bbb4Xhw+G//9sXDJHSVAsdYD8z6w9cAPQp6Zr8/Pzv7+fk5JCT\nkxN5LpFMcvrp0KsXXHQRdOsG998Pv/hF6FSSTAUFBRQUFCTlvaIeg+iBH1PIjT2+FnDOuVuLXdcJ\nmA3kOufeKeG9NAYhkiTOwfTpfnfYiy/2LYpatUKnkiik8hjEaqCdmbUxsxrAaGBe0QvMrDW+OIwr\nqTiISHKZwXnnwbp1sGmTH5tYuTJ0Kkk1ka+DiM1MuhNfjKY6524xswn4lsRkM5sCjAS2Agbscc51\ni/M+akGIRGT2bLjiCr+N+J/+pJlOmUTbfYtIwj79FK66Clas8CuyTzstdCJJBhUIEUmaJUvg8sv9\nVuJ33gk//WnoRJKIVB6DEJE0c/LJfmyiRw9/1sRNN8G334ZOJSGoQIjIAWrWhN//3h9E9OqrcOyx\nsGhR6FRS2dTFJCIHtWABXHklHHcc3Hab315c0oO6mEQkUnl58OabcMIJvuvpqqv8oLZkNhUIESmT\n2rXhuutg/Xr45hs4+mh/5sR334VOJlFRgRCRcmnWDO69F5Yvh8WL4ZhjYM4cHU6UiTQGISIJWbIE\nfvtbqFMH/vAHGDAgdCIpSusgRCSowkJ45BGYNAlatvSFonfv0KkENEgtIoFVrQpjxsBbb/nDic45\nBwYPhldeCZ1MEqECISJJU60aXHABbN4MQ4f6syeGD1ehSFcqECKSdDVqwKWXwpYt0L8/jBjhV2g/\n84wGs9OJxiBEJHLffQczZvgT7Ro08NNlhw6FKvpf1MhpkFpE0kJhIcydCzffDF9/Db/7nR+7qFkz\ndLLMpQIhImnFOVi6FP78Z78x4CWX+C6pww8PnSzzaBaTiKQVMz8msXixH5f4+GPo0AHGjoWXXw6d\nTvZTC0JEUsJnn8H998Pdd0Pz5nDZZXD66X6LD6k4dTGJSMYoLIT58+F//xdWr/ZjFOPHQ6dOoZOl\nJxUIEclI27b5VsX99/tWxfjxMHo01KsXOln6UIEQkYxWWOjHK6ZM8ZsEnnqqX6190kl+cZ6UTAVC\nRLLGzp1+36cHH4StW2HUKF8sTjjBD37Lj6lAiEhW2rIFZs70i/Ccg7PP9gPbxx6rYrGfCoSIZDXn\n/ID2rFnw+ONQvTqMHOmLRba3LFQgRERinIPXXoPZs/1t926/F9Tw4dCnj98nKpuoQIiIlGDDBl8o\n5s+HTZv8gUZ5eX478hYtQqeLngqEiEgZ7NoFixbBU0/5v23a+EIxaBD06AG1aoVOmHwqECIi5bR3\nL7z0ki8WS5f6w466d4eBA30r4/jjM2MKrQqEiEiCvvgCnn0Wli3z+0Nt2+bHLPr08cendu2antt+\nqECIiCTZrl2wYgWsXAmrVsH69X76bK9e/ta7d3rsPqsCISISsa+/9lNp9xeMVav84UfdukGXLr5L\nqksXOOyw0El/TAVCRKSS7dvnZ0W98gqsWQOvvgpr10KjRr5Q7C8anTtDs2bh1mKoQIiIpIB9++Dt\nt32xWLPG39au9cXhmGOgY0f/d/+tcePoM6lAiIikKOf8/lFvvunHMdav/+F+7dq+UBx1FLRvD+3a\n+b9t2ybvGFYVCBGRNOMc7NjhC8WWLT++bdvmF/G1b/9DwWjd2q/baNMGmjYte5eVCoSISAbZswfe\nf/+HgvH++37n2v233bt9wShaNNq0gSOO8IWlZUuoW9e/V0oXCDPLBf6CP/96qnPu1jjX3AUMBnYD\n5zvn1sa5RgVCRARfILZt88Vi/9+tW2H7dt8q2bHDT8l96aUULhBmVgXYDAwEPgRWA6OdcxuLXDMY\nuNw5l2dm3YE7nXM94ryXCkQSFRQUkJOTEzpGxtD3mTz6LhPnnJ+WW6dOYgWiSrKDFdMN2OKc2+qc\n2wPMAoYVu2YYMB3AOfcS0MDMmkWcK+sVFBSEjpBR9H0mj77LxJn54pCoqAtES+CDIo+3x54r7Zod\nca4REZFKFnWBEBGRNBX1GEQPIN85lxt7fC3gig5Um9m9wHLn3MOxxxuBfs65ncXeSwMQIiIVUNEx\niKg3s10NtDOzNsBHwGhgTLFr5gGXAQ/HCsrnxYsDVPxfUEREKibSAuGcKzSzy4HF/DDNdYOZTfAv\nu8nOuafMbIiZvY2f5npBlJlERKRs0mahnIiIVK6UG6Q2s1wz22hmm83smhKuucvMtpjZWjPrXNkZ\n08nBvk8z62dmn5vZmtjt+hA504GZTTWznWa2rpRr9Nssg4N9l/pdlo+ZtTKzZ8xsvZm9YWZXlnBd\n+X6fzrmUueEL1ttAG6A6sBY4utg1g4EFsfvdgRdD507VWxm/z37AvNBZ0+EG9AE6A+tKeF2/zeR9\nl/pdlu/7bA50jt2vC2xKxn87U60FoYV1yVWW7xNAEwDKwDn3PPBZKZfot1lGZfguQb/LMnPO/Z+L\nbVHknPsK2MCB68nK/ftMtQKhhXXJVZbvE6BnrMm5wMw6VE60jKTfZnLpd1kBZvYTfOvspWIvlfv3\nGfU0V0l9rwKtnXNfx/bFmgscFTiTiH6XFWBmdYHHgF/HWhIJSbUWxA6gdZHHrWLPFb/miINcI95B\nv0/n3FfOua9j958GqptZip2qmzb020wS/S7Lz8yq4YvDP51zT8S5pNy/z1QrEN8vrDOzGviFdfOK\nXTMPOBe+X6kdd2GdAGX4Pov2QZpZN/zU508rN2ZaMUruG9dvs3xK/C71u6yQ+4G3nHN3lvB6uX+f\nKdXF5LSwLqnK8n0CZ5jZpcAe4BtgVLjEqc3MZgI5QCMz2wZMAmqg32a5Hey7RL/LcjGz3sA5wBtm\n9hrggN/jZzBW+PephXIiIhJXqnUxiYhIilCBEBGRuFQgREQkLhUIERGJSwVCRETiUoEQEZG4VCBE\nEmBmDWLz9UUyjgqESGIOBX4VOoRIFFQgRBJzM9A2dqjNraHDiCSTVlKLJMDM2gDznXOdQmcRSTa1\nIEREJC4VCBERiUsFQiQxXwL1QocQiYIKhEgCYmcUrDSzdRqklkyjQWoREYlLLQgREYlLBUJEROJS\ngRARkbhUIEREJC4VCBERiUsFQkRE4lKBEBGRuFQgREQkrv8Htv2R2Ngj0ZQAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "from scipy.integrate import odeint\n", "import numpy as N\n", "\n", "def f(y, t):\n", " \"\"\"this is the rhs of the ODE to integrate, i.e. dy/dt=f(y,t)\"\"\"\n", " return -2 * y * t\n", "\n", "y0 = 1 # initial value\n", "a = 0 # integration limits for t\n", "b = 2\n", "\n", "t = N.arange(a, b, 0.01) # values of t for\n", " # which we require\n", " # the solution y(t)\n", "y = odeint(f, y0, t) # actual computation of y(t)\n", "\n", "import pylab # plotting of results\n", "pylab.plot(t, y)\n", "pylab.xlabel('t'); pylab.ylabel('y(t)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `odeint` command takes a number of optional parameters to change the default error tolerance of the integration (and to trigger the production of extra debugging output). Use the help command to explore these:\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```python\n", "help(scipy.integrate.odeint)\n", "```\n", "\n", "will show:\n", "\n", "```\n", "Help on function odeint in module scipy.integrate.odepack:\n", "\n", "odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)\n", " Integrate a system of ordinary differential equations.\n", " \n", " Solve a system of ordinary differential equations using lsoda from the\n", " FORTRAN library odepack.\n", " \n", " Solves the initial value problem for stiff or non-stiff systems\n", " of first order ode-s::\n", " \n", " dy/dt = func(y, t0, ...)\n", " \n", " where y can be a vector.\n", " \n", " *Note*: The first two arguments of ``func(y, t0, ...)`` are in the\n", " opposite order of the arguments in the system definition function used\n", " by the `scipy.integrate.ode` class.\n", " \n", " Parameters\n", " ----------\n", " func : callable(y, t0, ...)\n", " Computes the derivative of y at t0.\n", " y0 : array\n", " Initial condition on y (can be a vector).\n", " t : array\n", " A sequence of time points for which to solve for y. The initial\n", " value point should be the first element of this sequence.\n", " args : tuple, optional\n", " Extra arguments to pass to function.\n", " Dfun : callable(y, t0, ...)\n", " Gradient (Jacobian) of `func`.\n", " col_deriv : bool, optional\n", " True if `Dfun` defines derivatives down columns (faster),\n", " otherwise `Dfun` should define derivatives across rows.\n", " full_output : bool, optional\n", " True if to return a dictionary of optional outputs as the second output\n", " printmessg : bool, optional\n", " Whether to print the convergence message\n", " \n", " Returns\n", " -------\n", " y : array, shape (len(t), len(y0))\n", " Array containing the value of y for each desired time in t,\n", " with the initial value `y0` in the first row.\n", " infodict : dict, only returned if full_output == True\n", " Dictionary containing additional output information\n", " \n", " ======= ============================================================\n", " key meaning\n", " ======= ============================================================\n", " 'hu' vector of step sizes successfully used for each time step.\n", " 'tcur' vector with the value of t reached for each time step.\n", " (will always be at least as large as the input times).\n", " 'tolsf' vector of tolerance scale factors, greater than 1.0,\n", " computed when a request for too much accuracy was detected.\n", " 'tsw' value of t at the time of the last method switch\n", " (given for each time step)\n", " 'nst' cumulative number of time steps\n", " 'nfe' cumulative number of function evaluations for each time step\n", " 'nje' cumulative number of jacobian evaluations for each time step\n", " 'nqu' a vector of method orders for each successful step.\n", " 'imxer' index of the component of largest magnitude in the\n", " weighted local error vector (e / ewt) on an error return, -1\n", " otherwise.\n", " 'lenrw' the length of the double work array required.\n", " 'leniw' the length of integer work array required.\n", " 'mused' a vector of method indicators for each successful time step:\n", " 1: adams (nonstiff), 2: bdf (stiff)\n", " ======= ============================================================\n", " \n", " Other Parameters\n", " ----------------\n", " ml, mu : int, optional\n", " If either of these are not None or non-negative, then the\n", " Jacobian is assumed to be banded. These give the number of\n", " lower and upper non-zero diagonals in this banded matrix.\n", " For the banded case, `Dfun` should return a matrix whose\n", " rows contain the non-zero bands (starting with the lowest diagonal).\n", " Thus, the return matrix `jac` from `Dfun` should have shape\n", " ``(ml + mu + 1, len(y0))`` when ``ml >=0`` or ``mu >=0``.\n", " The data in `jac` must be stored such that ``jac[i - j + mu, j]``\n", " holds the derivative of the `i`th equation with respect to the `j`th\n", " state variable. If `col_deriv` is True, the transpose of this\n", " `jac` must be returned.\n", " rtol, atol : float, optional\n", " The input parameters `rtol` and `atol` determine the error\n", " control performed by the solver. The solver will control the\n", " vector, e, of estimated local errors in y, according to an\n", " inequality of the form ``max-norm of (e / ewt) <= 1``,\n", " where ewt is a vector of positive error weights computed as\n", " ``ewt = rtol * abs(y) + atol``.\n", " rtol and atol can be either vectors the same length as y or scalars.\n", " Defaults to 1.49012e-8.\n", " tcrit : ndarray, optional\n", " Vector of critical points (e.g. singularities) where integration\n", " care should be taken.\n", " h0 : float, (0: solver-determined), optional\n", " The step size to be attempted on the first step.\n", " hmax : float, (0: solver-determined), optional\n", " The maximum absolute step size allowed.\n", " hmin : float, (0: solver-determined), optional\n", " The minimum absolute step size allowed.\n", " ixpr : bool, optional\n", " Whether to generate extra printing at method switches.\n", " mxstep : int, (0: solver-determined), optional\n", " Maximum number of (internally defined) steps allowed for each\n", " integration point in t.\n", " mxhnil : int, (0: solver-determined), optional\n", " Maximum number of messages printed.\n", " mxordn : int, (0: solver-determined), optional\n", " Maximum order to be allowed for the non-stiff (Adams) method.\n", " mxords : int, (0: solver-determined), optional\n", " Maximum order to be allowed for the stiff (BDF) method.\n", " \n", " See Also\n", " --------\n", " ode : a more object-oriented integrator based on VODE.\n", " quad : for finding the area under a curve.\n", " \n", " Examples\n", " --------\n", " The second order differential equation for the angle `theta` of a\n", " pendulum acted on by gravity with friction can be written::\n", " \n", " theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0\n", " \n", " where `b` and `c` are positive constants, and a prime (') denotes a\n", " derivative. To solve this equation with `odeint`, we must first convert\n", " it to a system of first order equations. By defining the angular\n", " velocity ``omega(t) = theta'(t)``, we obtain the system::\n", " \n", " theta'(t) = omega(t)\n", " omega'(t) = -b*omega(t) - c*sin(theta(t))\n", " \n", " Let `y` be the vector [`theta`, `omega`]. We implement this system\n", " in python as:\n", " \n", " >>> def pend(y, t, b, c):\n", " ... theta, omega = y\n", " ... dydt = [omega, -b*omega - c*np.sin(theta)]\n", " ... return dydt\n", " ...\n", " \n", " We assume the constants are `b` = 0.25 and `c` = 5.0:\n", " \n", " >>> b = 0.25\n", " >>> c = 5.0\n", " \n", " For initial conditions, we assume the pendulum is nearly vertical\n", " with `theta(0)` = `pi` - 0.1, and it initially at rest, so\n", " `omega(0)` = 0. Then the vector of initial conditions is\n", " \n", " >>> y0 = [np.pi - 0.1, 0.0]\n", " \n", " We generate a solution 101 evenly spaced samples in the interval\n", " 0 <= `t` <= 10. So our array of times is:\n", " \n", " >>> t = np.linspace(0, 10, 101)\n", " \n", " Call `odeint` to generate the solution. To pass the parameters\n", " `b` and `c` to `pend`, we give them to `odeint` using the `args`\n", " argument.\n", " \n", " >>> from scipy.integrate import odeint\n", " >>> sol = odeint(pend, y0, t, args=(b, c))\n", " \n", " The solution is an array with shape (101, 2). The first column\n", " is `theta(t)`, and the second is `omega(t)`. The following code\n", " plots both components.\n", " \n", " >>> import matplotlib.pyplot as plt\n", " >>> plt.plot(t, sol[:, 0], 'b', label='theta(t)')\n", " >>> plt.plot(t, sol[:, 1], 'g', label='omega(t)')\n", " >>> plt.legend(loc='best')\n", " >>> plt.xlabel('t')\n", " >>> plt.grid()\n", " >>> plt.show()\n", "\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise: using `odeint`\n", "\n", "1. Open a new file with name `testodeint.py` file in a text editor.\n", "\n", "2. Write a program that computes the solution *y*(*t*) of this ODE using the `odeint` algorithm:\n", " $$\\frac{\\mathrm{d}y}{\\mathrm{d}t} = -\\exp(-t)(10\\sin(10t)+\\cos(10t))$$\n", " from $t=0$ to $t = 10$. The initial value is $y(0)=1$.\n", "\n", "3. You should display the solution graphically at points $t=0$, $t=0.01$, $t=0.02$, ..., $t=9.99$, $t=10$.\n", "\n", "Hint: a part of the solution $y(t)$ is shown in the figure below.\n", "\n", "\"image\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Root finding\n", "------------\n", "\n", "If you try to find a $x$ such that\n", "$$f(x)=0$$\n", "then this is called *root finding*. Note that problems like $g(x)=h(x)$ fall in this category as you can rewrite them as $f(x)=g(x)−h(x)=0$.\n", "\n", "A number of root finding tools are available in `scipy`’s `optimize` module.\n", "\n", "### Root finding using the bisection method\n", "\n", "First we introduce the `bisect` algorithm which is (i) robust and (ii) slow but conceptually very simple.\n", "\n", "Suppose we need to compute the roots of *f*(*x*)=*x*3 − 2*x*2. This function has a (double) root at *x* = 0 (this is trivial to see) and another root which is located between *x* = 1.5 (where *f*(1.5)= − 1.125) and *x* = 3 (where *f*(3)=9). It is pretty straightforward to see that this other root is located at *x* = 2. Here is a program that determines this root numerically:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The root x is approximately x= 2.00000023842,\n", "the error is less than 1e-6.\n", "The exact error is -2.38419e-07.\n" ] } ], "source": [ "from scipy.optimize import bisect\n", "\n", "def f(x):\n", " \"\"\"returns f(x)=x^3-2x^2. Has roots at\n", " x=0 (double root) and x=2\"\"\"\n", " return x ** 3 - 2 * x ** 2\n", "\n", "# main program starts here\n", "x = bisect(f, 1.5, 3, xtol=1e-6)\n", "\n", "print(\"The root x is approximately x=%14.12g,\\n\"\n", " \"the error is less than 1e-6.\" % (x))\n", "print(\"The exact error is %g.\" % (2 - x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `bisect()` method takes three compulsory arguments: (i) the function *f*(*x*), (ii) a lower limit *a* (for which we have chosen 1.5 in our example) and (ii) an upper limit *b* (for which we have chosen 3). The optional parameter `xtol` determines the maximum error of the method.\n", "\n", "One of the requirements of the bisection method is that the interval \\[*a*, *b*\\] has to be chosen such that the function is either positive at *a* and negative at *b*, or that the function is negative at *a* and postive at *b*. In other words: *a* and *b* have to enclose a root.\n", "\n", "### Exercise: root finding using the bisect method\n", "\n", "1. Write a program with name `sqrttwo.py` to determine an approximation of $\\sqrt{2}$ by finding a root *x* of the function $f(x)=2 − x^2$ using the bisection algorithm. Choose a tolerance for the approximation of the root of 10−8.\n", "\n", "2. Document your choice of the initial bracket $[a, b]$ for the root: which values have you chosen for *a* and for *b* and why?\n", "\n", "3. Study the results:\n", "\n", " - Which value for the root *x* does the bisection algorithm return?\n", "\n", " - Compute the value of $\\\\sqrt{2}$ using `math.sqrt(2)` and compare this with the approximation of the root. How big is the absolute error of *x*? How does this compare with `xtol`?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Root finding using the `fsolve` funcion\n", "\n", "A (often) better (in the sense of “more efficient”) algorithm than the bisection algorithm is implemented in the general purpose `fsolve()` function for root finding of (multidimensional) functions. This algorithm needs only one starting point close to the suspected location of the root (but is not garanteed to converge).\n", "\n", "Here is an example:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The root x is approximately x= 2.000000000000006661\n", "The exact error is -6.66134e-15.\n" ] } ], "source": [ "from scipy.optimize import fsolve\n", "\n", "def f(x):\n", " return x ** 3 - 2 * x ** 2\n", "\n", "x = fsolve(f, 3) # one root is at x=2.0\n", "\n", "print(\"The root x is approximately x=%21.19g\" % x)\n", "print(\"The exact error is %g.\" % (2 - x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The return value[6] of `fsolve` is a numpy array of length *n* for a root finding problem with *n* variables. In the example above, we have *n* = 1.\n", "\n", "Interpolation\n", "-------------\n", "\n", "Given a set of *N* points $(x_i, y_i)$ with $i = 1, 2, …N$, we sometimes need a function $\\hat{f}(x)$ which returns $y_i = f(x_i)$ where $x == x_i$, and which in addition provides some interpolation of the data $(x_i, y_i)$ for all $x$.\n", "\n", "The function `y0 = scipy.interpolate.interp1d(x,y,kind=’nearest’)` does this interpolation based on splines of varying order. Note that the function `interp1d` returns *a function* `y0` which will then interpolate the x-y data for any given $x$ when called as $y0(x)$.\n", "\n", "The code below demonstrates this, and shows the different interpolation kinds." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEPCAYAAACneLThAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVNX7wPHPGTZZBRQVUBF3TVFcSNuEyvJrm21muaSW\ntpqZbVaGS5aV32+/suXbZmbZvmrfMi2j1DTNJTSXXAHFHVD2Zeb5/TGIqIDIAMPyvF+veTkz995z\nnxn14XLuOc8xIoJSSqm6yeLsAJRSSlUdTfJKKVWHaZJXSqk6TJO8UkrVYZrklVKqDtMkr5RSdZhr\nVZ/AGLMHOAbYgHwRiarqcyqllLKr8iSPPblHi0hqNZxLKaVUMdXRXWOq6TxKKaVOUx3JV4Alxpg1\nxpgx1XA+pZRShaqju+ZCEdlvjAnCnuy3iMjyajivUkrVe1We5EVkf+Gfh40xXwNRQFGSN8Zo8Ryl\nlKoAETFn26dKu2uMMV7GGJ/C597AFcCm0/cTEX2IEBsb6/QYaspDvwv9LvS7KPtRXlV9Jd8U+Lrw\nat0VmC8ii6v4nEoppQpVaZIXkd1A96o8h1JKqdLp0MYaJDo62tkh1Bj6XZyk38VJ+l2cO3MufTtV\nEoAx4uwYlFKqtjHGIOW48VodQyiVUjVUq1atSEhIcHYYqgxhYWHs2bOnwsfrlbxS9Vjh1aCzw1Bl\nKO3vqLxX8tonr5RSdZgmeaWUqsM0ySulVB2mSV4pVSuMGjWKp59+2tlhnNXAgQP54IMPnB1GEU3y\nSqk6JyYmhjlz5jjl3N9//z3Dhw8v177VEacOoVRKnWH37gQmT57Lvn02QkMtTJ8+kvDwsGpvQ1WC\nGlBkR5RSzlHS/79du/ZImzYTBTIERCBD2rSZKLt27Sl3u5XRxrp166RHjx7i5+cnt9xyiwwZMkQm\nT54sIiKpqaly9dVXS1BQkAQGBsrVV18t+/btExGRJ598UlxcXMTT01N8fX1l3LhxIiIyfvx4adGi\nhfj5+UmvXr1k2bJlpZ575MiRcvfdd0v//v3F19dXoqOjJSEhoWj7ihUrpHfv3uLv7y9RUVHy+++/\nF22Ljo6Wd999V0RE5s6dKxdddJE8/PDDEhAQIK1bt5ZFixaVGefpSsuRhe+fPceWZ6eqfGiSV8p5\nSvr/N3TolGLJWYqS9NChU8rdrqNt5OXlSVhYmLz88stSUFAgX3zxhbi5uRUl+aNHj8pXX30lOTk5\nkpGRIYMHD5ZBgwYVHV880Z4wf/58SU1NFavVKv/5z3+kWbNmkpubW+L5R44cKX5+frJ8+XLJy8uT\n8ePHy0UXXSQiIikpKRIQECDz588Xq9UqH3/8sQQEBEhKSsoZ5547d664u7vLu+++KzabTd544w0J\nCQkpM87TOZrktU9eKXWKfftsgPdp73qTnGyrtjZWrVpFQUEBDzzwAC4uLtx444307t27aHtgYCDX\nX389Hh4eeHt7M2nSJH777bcy27ztttvw9/fHYrEwYcIEcnNz2bZtW6n7X3XVVVx44YW4ubkxY8YM\nVq1axb59+/jf//5H+/btue2227BYLAwZMoSOHTuycOHCEtsJCwtj9OjRGGO4/fbb2b9/P4cOHSrX\n91AZNMkrpU4RGmoBMk97N5OQkPKnC0fbSE5OJjQ09JT3wsJO9udnZ2dz11130apVK/z9/enXrx9p\naWllzt6dNWsWnTt3JiAggICAAI4fP86RI0dK3b9FixZFz729vQkICCA5OZnk5ORTYjkR2759+0ps\np1mzZkXPPT09AcjIyCj1vJVNk7xS6hTTp4+kTZtYTibpTNq0iWX69JHV1kZwcPAZSTMxMbHo+axZ\ns9i+fTtr1qwhLS2t6Cr+RJI35tTZ/suXL+fFF1/kiy++IDU1ldTUVPz8/Mr8oZCUlFT0PCMjg9TU\nVEJCQggJCTmjlkxiYuIZP5TK4/Q4q4Im+UpiFSE1P5/EnByO5OVRYCv/r7ZK1STh4WEsWTKOoUNn\nERMTy9Chs1iyZNw5jYxxtI2+ffvi6urK7NmzKSgo4KuvvmL16tVF2zMyMvD09MTPz4+UlBSmTJly\nyvFNmzZl165dRa/T09Nxc3OjUaNG5OXlMW3aNNLT08uM4fvvv+f3338nLy+PyZMn06dPH0JDQxk4\ncCDbt2/nk08+wWq18umnn7Jlyxauueaacn8/pcVZJcrTcV+VD2rpjdccq1W+P3JE7tm2TXqsWSOe\nv/4qDX/7TUJXrJDAZcvENS5OOv7xhwzbvFk+OnBAjuXnOztkpc5Qk///rV27ViIjI8XPz0+GDBly\nyuia5ORkiY6OFh8fH+nQoYO89dZbYrFYxGq1iojIypUrpX379hIYGCjjx48Xm80mo0ePFj8/PwkJ\nCZEXX3xRwsPD5eeffy7x3CNHjpR77rlH+vfvLz4+PtKvXz/Zs+fkyKAVK1ZIz549xd/fX3r16nXK\n6JqYmJhTbrxefPHFp7RtsVhk586dJcZZktL+jijnjVetQnmODuTmMnvfPt7ev5/2np5c27gxlzRs\nSBdvb3xcT047yLfZ2JyVxR/Hj7PgyBGWHzvGzU2aMKF5czp7n35DSinn0CqUJRs1ahQtWrRg2rRp\nzg7F4SqUOhmqnDKtVmYlJfHK3r0MadKEZZGRdPDyKnV/N4uFbj4+dPPxYWxICIfz8vhvcjIxGzZw\nTaNGPBMeTjMPj2r8BEqp+kj75MtheVoa3dasYXNmJmt79uS19u3LTPAlCXJ3Z3KrVmyLiiLAzY2I\nP//kwwMH9CpKqRqoOm6IVhftrimDiDAzMZFX9u3jjXbtGBQUVGltr09PZ9iWLXT19ubdjh3xdnGp\ntLaVKi/trqn5HO2u0SR/mhP1NvYm20i+qT0u3cL5uVdPQqqgayXHauWe7dtZn57Ot127EtagQaWf\nQ6myaJKv+XRlqEq0e3cC/fvPZv7nD/FrzM1sz+5Ezt0LyE0+UCXna+DiwpwOHRjRrBkXrFvH35mn\nTx5RSinHVPmVvDFmAPB/2H+gvCsiz5+2vVKu5Hem7OTOhXditVkr1oAI1t930m1HAEuHPoB36hFa\nLP6ElUENkLBDdOxQsep5U6OnEhMec9b95h88yMM7d/J9165E+vpW6FxKnSu9kq/5avToGmOMBXgV\nuAxIBtYYY74Vka2Vfa7tKdvJzMtk1hWzyn2MJTsH/5XrCfxlFYG/rCLexZNrn5lOwJ/pXDlvDxdl\ne/Ne3t9keBrcr+zEsV5dON6rK1ltw8By9l+CXl/zOuv2rytXkh/atCleFgsD4uNZ0q0bET4+5f4c\nSilVmqoeQhkFbBeRBABjzCfAdUClJ/ncglyCfYO5JOySsnfctQu+/97+WL4cevSAgQNZM+FpopMP\nkf1uOw5+3ZrnuB0AQzqPRj/GzCu60XTFCnhnGqSmQt++cOGF9kfv3lDCaJslO5eQlZ9V7s9wfVAQ\n+SIMjI9nWWQk4YV1LpRSqqKqOsmHAknFXu/FnvgrXU5BDh4uJdwczcuDZcvgf/+zJ/a0NPjXv2DU\nKPjoI/D358vDh7n7n394uVULno9/nZ1MxV5BL5PWbaZy1+zHIDwM7rrL3ub+/fD777BiBTz6KGza\nBF26nEz6F14IzZrh6eZJem7ZU6dPN7hJEw7l5XFlfDy/R0bS2N3d4e9GKVV/1YjJUMXrTkRHRxMd\nHX3ObeRac2ngWjg6Zd++k1frS5dCp05w1VUwfz5ERhZ1tYgIzyck8Oq+fSyKiKCnry+XLxnH5Mmz\nSE62ERJiYfr0EuptBAfDjTfaHwBZWbBmjT3pz5kDY8ZAQACD2geyqZ0fNNsEnTuXq4sH4P7mzUnK\nzeWWzZv5MSIC13Iep5SqXtU5MzYuLo64uLhzPq5Kb7waY/oAU0RkQOHrx7HXW3i+2D6O33gtKODb\n95/AbdFiBm4HkpLgyith4ED7nyWMb0/Jz2fk1q0czMvji/POo0VlDl+02WDLFn77eCauK//ggkSB\nI0dO7eKJiiqxi+cEqwhXxcfT2dub/7RtW3mxKVWM3ngFq9WKSwXnqVRHknf0xmtVFx9zAXYAYYA7\nsAHodNo+JRbfOatDh0TmzRMZMkQkMFAOtQ+VRYN7iixfLnKWYmDLUlOl1cqVMmH7dsktLGhUFT74\n6wMZ+uVQ+4sDB0S++kpk4kSRPn1EvLxEevcWefBBkc8/F0lOPuP4o3l50mblSvm/jZtk6NApEh39\ntAwdOuWcllBTqiwV/v9XTVq1aiWzZs2SiIgI8ff3lyFDhhSt5rRw4ULp3r27+Pv7y4UXXijx8fFF\nx82cOVPatGkjvr6+ct5558nXX39dtG3u3Lly4YUXyoQJE6RRo0ZFRc/effdd6dSpkwQGBsqAAQNO\nWe7vwQcflCZNmoifn59ERETI33//LW+99Za4ubmJh4eH+Pr6yrXXXlsl30Fpf0fUlOX/gAHANmA7\n8HgJ28v3Sa1WkTVrRKZOFYmKEvHzE7n+epG33xbZu1deWP6CTPxxYplNZBQUyLh//pGQFSvk28OH\ny3deB3zx9xdy/SfXl7wxK0vkt99EnntO5OqrRQIDRcLDRYYNE3njDZH4eBGrVb7buk0sCxYJwUcq\nvFamUqWpDUn+/PPPlwMHDkhqaqp06tRJ3nzzTVm/fr00adJE1qxZIzabTebNmyetWrWSvLw8ERH5\n4osv5MCBAyIi8tlnn4m3t3fR67lz54qrq6u89tprYrVaJScnR7755htp166dbNu2TaxWq8yYMUMu\nuOACERH58ccfpVevXnL8+HEREdm6dWtRWyNHjiz6IVFVHE3yVd4nLyKLgA4VOjgtDRYvtvet//AD\nBAbau2Ceew4uugiK3ZTM3VmsT/7MGPji8GEe3bWLixs2ZGPv3gS6uVUopHPh5ebF3uN7+XnXzyXv\nEAoM7m1/2B7Aa1cS/mv/puFP3+D//DO4HU2jScOGXBZ9P39MzSfreVcKbG7stEUz5tmnmTRpRIXi\nauXfijaBbSr+wVS9YaZWTg0Xia14l9D48eNp2rQpANdccw3r169nw4YN3H333fTq1QuA4cOHFy3R\nd/HFF3PjiftlwM0338yzzz7L6tWri2q+h4aGcu+99wLg4eHBm2++yaRJk2jfvj0Ajz/+ODNmzCAp\nKQk3NzfS09PZvHkzUVFRdOhQsXTmLDXixuuwYVOZPn3kqTc433gDHnsMLr7Yntiffhpaty61jZyC\nnBKT/PK0NB7ftYsMq5V3O3Tg0oCAKvgEJevQuAP+Dfx5dvmz5T8otPBxbQf8j+fR4MdtXLb+TTJ6\njCPiiu2s3/Mnq4MastZ7F88u33vOMaVmpxLgGcDPI0r5waNUMY4k58pyIsEDeHl5kZycTEpKCu+/\n/z6zZ88G7Bdy+fn5JCcnAzBv3jxeeumlohWcMjMzT1nqr/jSfgAJCQmMHz+eiRMnFrVnjGHfvn3E\nxMRw//33c99995GYmMgNN9zArFmz8Kklc1lqRJKfP/9hVq2KPXXlmFtugdtvL/PmZHG5Bbk09GgI\n2G9a/nD0KDMTE9mfl8dTYWGMaNYMl2quLNc6oDWLhy92qI1hq6by+N8Pw3RX1r+9nA+n/8OeH87j\n71v6MGfEjHNu77eE33hy6ZMOxaSUMxljaNmyJU899RSTJk06Y3tiYiJjx47ll19+oW/fvgBERkae\ncvPy9CqTJ9q79dZbSzzn/fffz/3338+RI0e4+eabefHFF5k6dWqtqFZZI5I8eLNz51QmT57Fhx/G\nAvDvre/xzLJnyt1CZn4WT175Dk/s2sW8Awdo5u7Owy1acFNQUK0egjh9+khWrYpl586p5DwfyeAn\npzLnoYm8ueoQxPWHcxxu6mJcKl76QakaYsyYMQwaNIjLLruMqKgoMjMz+fXXX+nXrx+ZmZlYLBYa\nN26MzWbj/fffZ9OmTWW2d9dddzF58mS6detG586dOXbsGEuWLOGmm27izz//xGaz0aNHDzw9PWnQ\noAGWwpxSLcv3OaiGJHkAb5KTT66LelvkWG7oOoKGLi5n/LQUEdKtNrZl5/B3Vg5/ZmbyS1oGc60u\n3GCz8UNEBF1rya9SZ3NircwTY/cP7W3N25+/wrAjR2D4cPv4/xdeAD+/crXnanGlwFZQxVErVTlK\nu1Lu0aMH77zzDvfffz87duzA09OTiy66iH79+tGpUycmTpxInz59cHFxYcSIEVx00UVlnmfQoEFk\nZmYyZMgQEhMTadiwIf379+emm27i+PHjTJgwgd27d9OgQQOuvPJKHnnkEQDuuOMObr75ZgIDA4mO\njuarr76q9O/AUTWi1DAIkMnQoSev5F/du5endu8mV4Smbm54WCwYIFeEg3l5GKCDlxfdfHzo6ePD\nFYGBtPP0rBW/PjmiwGbj8r/+IiYggNiAAHjkEVi0CP77X/u9i7P4M/lP7vruLtaOXVsN0aqaTsfJ\n13x1op48ZNCmTWyJq7lnWq0czMsjXwSbCO4WC03d3E5ZT7W+2Z+bS8+1a3mvY0euDAy0z+odMwYu\nuAD+7/+gUaNSj12/fz0jvx3JX3f/VY0Rq5pKk3zNVyfqyQ8dOqvEBA/g7eJCa09POnh50cnbmzae\nnvU6wQMEe3jwaefODN+yha2ZmXDppRAfD40b22vofP65fUh9CVwtrtonr1Q9UiOu5J0dQ201Z/9+\nZiYmsqpHj5Pj/leuhDvugI4d4bXX7HV2itl8eDM3fnYjW+7b4oSIVU2jV/I1X524klcVMzo4mGsa\nNWLw33+Tbyu8ad23L6xfby+I1q0bzJ17ylW93nhVqn7RJF/LvdCmDZ4uLozetg3biWTu4QHPPGOf\nLfzKK/bSygkJgA6hVKq+0SRfy7kYw6edO7MnJ4cJO3ac+mtd9+7wxx/Qrx/07AmvvYYLBqtokleq\nvtA++ToiLT+ffhs2cGNQEE+3anXmDlu3wh13kEsBl1+yh2XPHaz2GFXNo33yNZ/2ySsA/N3c+DEi\ngg8PHmRGYdfMKTp2hN9+I+e6q/nmlcP2CVQF2jevVF2nSb4Oaebhwa/du/PRwYM8sWvXmT/9XVzI\nu+9urnjA395f36ePfeilUjVUeHg4S5cu5bnnnmPs2LHODqdWqt8DzuugYA8P4rp354r4eDKtVl5q\n2xZLsVnALhYXdvkLLFliX6rw8svhnnvgiSfsN2yVqoFKKkSmykev5OugIHd3lnbrxrqMDG76+28y\nrSdvtBZNhjLGPp5+wwb7o2dPWL3aiVErVfPYbLaz71TDaZKvowLc3PipWzd8XVzot349+3JzAfsQ\nylPGyYeEwDffwFNPwbXXwsMP2xcmV6oGmTp1KsOHDwfstd8tFgvz5s0jLCyMJk2a8OyzJ9dsEBFm\nzpxJ27ZtCQoKYsiQIaSmphZtHzx4MMHBwQQEBBAdHc3mzZuLto0aNYp7772Xq666Cl9f3wotnF3T\naJKvwzwsFuZ27MhNQUH0WruWRUeP4mJxOXMIpTEwZAhs3AjJyRARAb/+6pyglSrF6cUHV6xYwfbt\n2/npp5+YNm0a27ZtA+CVV15hwYIFLFu2jOTkZAICArjvvvuKjhs4cCA7d+7k0KFD9OjRg6FDh57S\n7scff8zkyZNJT08/a/XKWqE8awRW5YMavsZkXRGXmirNf/9dJm7/RyzTGpS984IFIqGhInffLXLs\nWPUEqJzirP//7POlHX9UUKtWreTnn3+WKVOmyPDhw0VEZM+ePWKxWCQ5Oblov6ioKPn0009FRKRT\np06ydOnSom3Jycni5uYmVqv1jPZTU1PFGFO0fuvIkSPl9ttvr3C8VaG0vyPKucarXsnXE/38/Vnf\nsyc7s3Ow9fgvvxb79fUM11wDmzbZh1h26WJfY1fVT5WV5qvA6csCZmRkAPbunOuvv57AwEACAwPp\n3Lkzbm5uHDx4EJvNxuOPP07btm3x9/cnPDwcY0yZSwPWdprk65HG7u583bUrZs+7DNuyhdu3bGFP\ndnbJO/v7w9tvw3vvwf33w4gRcPRo9QasVAW0bNmSH374gZSUFFJSUkhNTSUzM5Pg4GA++ugjFi5c\nyNKlS0lLS2PPnj3FexWA0hcqqa00yddDrimrWN+zO2ENGtBz7Vru/+ef0pP9ZZfZ++oDA+1X9V98\nUb3BKlUCKeO3g7vuuosnnniCxMREAA4fPsyCBQsASE9Px8PDg4CAADIzM5k0aVKdS+qn0yRfD7lY\nXPCywLTwcLZGReHl4kKvtWu5Oj6eBUeOkGM97cast7d9MZIvv4TJk+HGG0n8Yw3Dhk0lJiaWYcOm\nsnt3CbNslXJQaQn49PeLvx4/fjzXXXcdV1xxBQ0bNuSCCy5gdeHw4BEjRtCyZUtCQ0Pp0qULF1xw\nQdUFX0NUWe0aY0wsMAY4VPjWEyKyqIT9pKpiUCXzfc6X5IeS8fXwLXovy2rl00OHeO/AAf7KyODy\ngAAubtiQ8/386OTlhf+JevU5OaRNfBjrf+fwkO0l5nmNhkbHCI54m0eeu5b8hr4cKyggw2ol3Wol\nw2qloPivwoC7xYKviwt+Li74uroS4OpKsLs7oR4ehLi7E+zhgUctXny9NtHaNTVfjV3+rzDJp4vI\nf86ynyb5atb65dbsTttd+g5uDSHwfPDrDL6dwLM5YIO8NJB8EBtuxgcXFz9cBJqkHMb/2FG8M47i\nnnMUqzWdHJNDliWbdJdsjrkXcKwB2CyF/x4tbuDiBS6e4OoNrj4EN4ogrEkPknNzOZCXR7C7Ox28\nvOjg5UVHLy8ifXzo5uODl4tLtXxH9YUm+ZrP0SRf1WUN6nZnVy2184Gd57S/iJBaUMCR/HxybTZG\nj3mDtb/cRn66K97ZWbhwEB/2c3nnV5g89hI4cMD+OJgFB9Lszw8ftt/MbdbM/gj2h6ZNoVkz9noV\n8PQf03nv7h+gVSsK/P3Zk5PDtuxstmZlsS49nXf372dLVhbtPD3p7evL+X5+XBoQQOsGDYp+Vd+9\nO4HJk+eyb5+N0FAL06ePLHFJSaXqk6pO8vcbY4YDfwITReRYFZ9PlcO53mgyxtDI3Z1G7u4AdHTN\nZu0hG+BGJg3ZQUN2EErLyAjM+PElN2K12kfn7N9/8ofAgQOwdy9N9iYwel0yLLkRDhzANSuLtk2b\n0jY4mKtO/FBo1oyc4GA2hoSwJjeXuOPHeXr3btwtFi4NCKBrfgGzh89nz19PAd5AJqtWlbw4vFL1\niUPdNcaYJUDT4m8BAjwJrAKOiIgYY54BgkXkjhLakNjY2KLX0dHRREdHVzgmVfV2706gf//Z7Nw5\nlRMJtU2biidUEcHjGQ/u6XUPxhhc8wrwTc3CLzUL39RM/FKz8EvJLHrum5KFX2omPqlZbA5vxaK+\nUSzq0pXVnbsRkpiM/+Y9rD1ohePH6dBhLQMG9KnQ5/R282ZazDRcLHW3i0i7a2q+E39HcXFxp5RZ\nmDp1qnP75E85iTFhwEIRiShhm/bJ10InukaSk22EhDjeNfLppk/Zn7H/3A4SwSMjB6+jx1n87k/k\n7utCXrsATOeG/Hp+D7LSbPjvXcFDlzWngTn3QlNP/PwEex/aS6Bn4DkfW1tokq/5avKN12YicqDw\n+QSgt4jcVsJ+muSVw4YNm8r8+Q8D3rhQwGOW5+gYtYJZw28moWt7bggKYkLz5nT18Sl3m0EvBrH5\n3s0EeQdVXeBOpkm+5qvJK0O9YIyJN8ZsAPoBE6rwXKqemz59JG3axAKZWHHlWdtDfLE3iDUvzWTb\nnDm0FeHK+Hj+FR/P0tTUciU2F1NCMTelahld41XVGSV2ITVrYi+j/Omn5L79Nh92786spCT8XV15\nLjyc6ICAUtsL/U8oq+9cTahfaPV9iGqmV/I1X43trikvTfKqWsTFwciRMGAA1hdf5JOsLJ7evZu2\nnp4827o1PX19zzik5UstWTZqGWH+dXd0Tl1M8jExMQwfPpzRo0efsS0pKYnzzjuPY8eO1ZpyBjW5\nu0apmiM62r6ebV4eLpGRDN2xgy1RUQxq3JhrNm7kzq1bOZyXd8ohJdbeV7VaixYtOH78eK1J8JVB\nk7yqP/z87Ova/vvfcPPNuD/xBPc0bszWqCj8XF05b80aXtu3D2vhVVPRUolK1WKa5FX9c9118Ndf\n8M8/0Ls3fn//zX/atuXnbt347NAh+q5bx5bMTL3xWgPs3buXG2+8kSZNmhAUFMQDDzxwylKAcHI5\nwOLrse7YsYPzzz+fhg0bcv3115OWllbivqmpqYwePZrQ0FAaNWrEDTfcUL0fsBpoklf1U1CQvarm\nww9D//4wcyZdPT2J696d0c2acfH69aQGXUmeteDsbakqYbPZuPrqqwkPDychIYF9+/YxZMgQoOwq\nlAAffPABc+fO5cCBA7i4uDBu3LgS9x02bBjZ2dls2bKFQ4cOMWFC3RsEqDdelUpMhFGjICcH5s2D\nNm3YnZ1N16Uf0KZRJ77u1ovWnp7OjrJKnO3Gq6mkhaylArPYV61axXXXXcf+/fuxFKtKOnXqVHbu\n3Mm8efMA+9V569atyc/Px2KxEBMTQ9++fYsW996yZQvdu3cnJyeHxMTEon0PHjxIixYtSElJwc/P\nr1I+Z1Wo6QXKlKr5WraEJUvg1VehTx945hnCx46l3d43uazVq/RZt47X27XjpiZNnB1ptatIcq4s\nSUlJhIWFnZLgy6v4En5hYWHk5+efssQf2LuCAgMDa3SCrwzaXaMUgMUCDzwAv/0G77wDAwfS7HgB\nt/i78n3Xrjy2axf3/vPPmQuqqCrTokULEhMTT+lrB/D29iYrK6vo9f79Z5bDSEpKKnqekJCAu7s7\njRs3PqP9lJQUjh8/XsmR1yya5JUqrlMn+P136NOH+dM34//tj/Ty82Ndr14czc/nwvXrScrJcXaU\n9UJUVBTBwcE8/vjjZGVlkZuby++//0737t357bffSEpK4tixY8ycOfOMYz/88EO2bt1KVlYWsbGx\n3HzzzUV98Se6Ppo1a8a//vUv7r33XtLS0igoKGDZsmXV+hmrgyZ5pU7n5gaxsTzyYGdC//MO3Hor\nDY8f55POnbm1SRPOX7eOFce0anZVs1gsLFy4kO3bt9OyZUtatGjBZ599xuWXX87gwYOJiIigd+/e\nXHPNNacArkweAAAeOElEQVQcZ4xh+PDh3H777YSEhJCXl8fLL798yvYTPvjgA1xdXenYsSNNmzY9\nZb+6Qm+8KlWKS967hBl9n+LiN3+Azz+3d+MMGMCio0cZsXUrM8LDGRMS4uwwHVIXZ7zWNTrjVakq\n4mpxJd/DFV56CT74AO6+G+65hwENGrAsMpJ/JyXx0I4d2DRJqhpMk7xSpXCxuJyc8RoTY59AlZMD\n3bvTYcMGVvbowdr0dG7bvJlc27nXq1eqOmiSV6oUZ8x4bdgQ3nsPXnwRbryRgNhYfuzYkQIRBsTH\nc6xAJ06pmkeTvFKlcLG4UGArIXEPGmS/qt+8mQZ9+vCp1UoXb28uXr+e5Nzc6g9UqTJokleqFGUW\nKGvSBL7+Gh56CJf+/Xnl668Z0rgx/TZs0CGWqkbRJK9UKc5aoMwYuP12WLMGs2gRTwwfzr0eHvTb\nsIHd2dnVF6hSZdCyBkqV4pQbr2UJC4OffoJXXmFC//64v/oq0SL81K0b7by8qj5QB4SFhdWr2uq1\nUViYY4vWaJJXqhQuppQ++ZJYLPDgg3DFFdw3YgQel19OjNXKT5GRdPT2rtpAHbBnzx5nh6CqmCZ5\npUpRoZWhOneGlSu5c8YM3F5+mZixY+n1zXoyNmcSGlq47mx43V1OUNU82ievVCkqvDKUmxtMmcLl\nQ4fz0P+9wYrozqzbeRfz5z9M//6z2b07ofKDVaoUWtZAqVKM+34cczbMoYFrgwodn5GRjSWjAddn\nXcDS/iNp8+Nkfg84hLtHNj4+Fa9P/8H1HzCw3cAKH6/qhvKWNdAkr1QpCmwFHMupeCGy666byYoV\njwPQMWYVqaNyGfPIMla1tvHJwqcr1OaDPz7IJS0vYUzPMRWOS9UN1bJoiDHmJmAK0AnoLSLrim2b\nBIwGCoDxIrLYkXMpVd1cLa408mpU4eNbNfVhRXYDwJut319FA+suXn3BgwVPT6LRX0Ogb99zbtPb\nzVvXnVXnxNE++Y3A9cCvxd80xnQCBmNP/v8CXjc6TkvVM9Onj6RNm1ggE4CcH5vitnA7g1/5N1vv\nuw+efBLy8s6pTRdTzmGdShVy6EpeRLYBlJDArwM+EZECYI8xZjsQBfzhyPmUqk3Cw8NYsmQckyfP\nIjnZRkiIhelTR/KbpweXvfIKP7/zDh2jouwVLrt2LVebpZZaUKoUVTWEMhRYWez1vsL3lKpXwsPD\n+PDD2FPfK/zzsrFj+XnXLjpeeik8+ig89BC4uJTZ3lln4Sp1mrMmeWPMEqBp8bcAAZ4UkYWVEcSU\nKVOKnkdHRxPtxMWDlaoOtzdrBsBlwM8rVtBx7FhYsADefx9aty71uAoP61S1XlxcHHFxced83FmT\nvIj0r0A8+4AWxV43L3yvRMWTvFL1RVGi37WLn777jk5vvw3nnw/PPgt33mmvjXOaCk3QUnXC6RfA\nU6dOLddxlTkZqvi/yAXAEGOMuzEmHGgLrK7EcylVJ9zerBnPt25N9F9/8duoUfDrr/Df/8LVV8P+\n/Wfsrzde1blyKMkbYwYZY5KAPsB3xpgfAERkM/AZsBn4HrhXB8MrVbJhzZoxv1Mnbvr7bz5p1AhW\nrYJevSAy0r62bDF641WdK50MpVQNsTEjg6s2buTekBAea9kSs2YNjBhhT/izZ0NAAFPjpmIVK9Ni\npjk7XOVkupC3UrVMVx8fVvbowaeHDzNy61aye/aEdeugUSOIiIAlS/TGqzpnmuSVqkFCPTxYHhlJ\nnggXrV9PgsUCL79sX1v2jjsY8PJ3WLJ0QRJVfprklaphvF1c+KhTJ4Y2bcr5a9fyc2oqXH45xMfT\nIDOPB+6Za++3V6ocNMkrVQMZY3ioRQs+7tyZYVu2MHXPHgr8/PgxdijfjbzAvpj4U0+dc1kEVf9o\nkleqBosJCGBtz54sS0sjesMGjuHJ+gtbw4YNEB9vH1e/aZOzw1Q1mCZ5pWq4EA8PFnfrxqDGjflP\nbju2mBBo1gy+/RbGjYOYGJg1C6x6Q1adSYdQKlWLPPHHe7yRGUh0kza82q4doR4esHs3jBxp32Hu\nXAgPL6sJVUfoEEql6qBw1wKuy/gfXb296f7nn7yZnIytVSv45Re49lqIioJ33wW9cFKFNMkrVYu4\nWFwwks+08HB+6daNuQcOcNH69azOyICJE+3J/rXX7An/wAFnh6tqAE3yStUixSdDdfHxYUVkJGOD\ng7l+0yZGbNnCvnbt7MMru3e3P7780skRK2fTJK9ULXJ6PXmLMYwMDmZrVBQtPDyIWLOGR5OSODR5\nsv3G7KRJMHw4pKU5MWrlTJrklapFSitQ5uvqyozWrfmrVy+yrFY6rl7Nw40acXDNGvD3t5dF+Okn\nJ0SsnE2TvFK1yNlKDTdv0IBX27cnvlcvcmw2Om7cyO333sufc+bA6NH2IZdZWezencCwYVOJiYll\n2LCp7N6dUI2fQlWnqlr+TylVBVwtruVaNOREsp8aHs67+/dzk7c3wV98wejFi7k26nweTovgq31v\nAd5AJqtWxbJkyTjCw8Oq/DOo6qXj5JWqRRZsW8CMZTOYfMnkczrOKrAmx4Vfst1Ym1bApX9uwG/l\ncb7c05uCAjcgh37RX/HwxNsqFNclYZfg5+FXoWNVxZR3nLwmeaVqke1HtzNx8URsYqtwGyvWJmHz\nvJzwZq3Y1SocS8IOjiXuJDDnJ/qe3+ac2/vr4F9Mi57GqMhRFY5JnbvyJnntrlGqFmnXqB0Lbl3g\nUBvDvp/K/PkP8xdeDAl4n4jotcy+dDBp4VfTq3N7xoaEEOLhUe72xi4cS74t36GYVNXRG69K1TPT\np4+kTZtYIItPUkfy9teP8c0j97Dszdc4nJJClzVruOXvv1mXnl6u9nTd2ZpNk7xS9Ux4eBhLloxj\n6NBZxMTEcsHQLwna+Bs9L7mE1wYMYPf27fT18+PajRu5Oj6eP44fL7M9XXe2ZtM+eaXUSZs22SdP\ntWhBzptv8p7NxnOJiXTx9ubFNm04z9v7jEPG/zCe8IBwHuzzoBMCrr+0QJlS6tx16QJ//AERETSI\njOSeP/5gx/nnc2VgIDEbNnD3tm0cPG2hEl13tmbTJK+UOpW7OzzzDHzzDTz+OO4jRzLex4etUVF4\nubhw3urVvLx3L9bC38BdLC7lGruvnEOTvFKqZH36wPr14OcHEREE/vor/2nblhU9evD14cP0XbeO\nDenpuBjtk6/JHEryxpibjDGbjDFWY0yPYu+HGWOyjDHrCh+vOx6qUqraeXvDq6/CO+/AqFHwwAN0\nAH7p3p27Q0K4Ij6eX106kGut+Lh9VbUcvZLfCFwP/FrCth0i0qPwca+D51FKOdMVV9jXlD16FHr0\nwKxZw+jgYDb27k2q8eZNurMpI8PZUaoSOJTkRWSbiGwHSrrDe9a7vkqpWiQgAObPh2nT4JprIDaW\npsZwK1vpZUskesMG/i8pCZuOlqtRqrJPvlVhV80vxpiLqvA8SqnqNHgwbNgAf/4JffrQLOkoXW1J\nrOrRg08OHWLQpk2k5usM2JrirGUNjDFLgKbF3wIEeFJEFpZyWDLQUkRSC/vqvzHGdBaREn+fmzJl\nStHz6OhooqOjyxe9Uso5goPhu+/gnXcYNm4CPw+Jou2lDfgtMpJHdu6k19q1fHneeXT39XV2pHVG\nXFwccXFx53xcpUyGMsb8AkwUkXXnul0nQylVu73z+RNcOnUerYPawXvvQatWfHLwION27ODF1q0Z\nGRzs7BDrJGdMhio6mTGmsTHGUvi8NdAW2FWJ51JK1RAZLZow+/kbYOBA6N0b3nuPIU2a8Gv37sxM\nTGTCjh1FY+pV9XN0COUgY0wS0Af4zhjzQ+GmS4B4Y8w64DPgLhHRRSaVqoNcjAv5xgaPPAJLl8LL\nL8OgQXTOyGBljx7EZ2Rw/aZNZBToWHpncHR0zTci0kJEPEUkWET+Vfj+VyLSpXD4ZC8R+b5ywlVK\n1TQulmJVKLt2tZdFOO886N6dgO++Y1FEBE3c3Lh4wwb25uQ4N9h6SGe8KqUccsaShB4e8Oyz8OWX\n8OijuI0axdtNm3Jrkyb0Xb9ex9NXM03ySimHlFpP/oIL7EMtfX0x3brx6NatPN+6NZf99ddZyxer\nyqNJXinlEBeLCwVSSn/7ibIIc+bAmDHcNmUKc8LCuGbjRn5KSaneQOspTfJKKYeUa2Woyy+3l0XI\nyeGqfv340mrlti1b+PLw4eoJsh7TNV6VUg45o0++NA0b2sfRL1jAxUOG8OPddzMQsIowuEmTKo+z\nvtIreaWUQ04ZXVMe114L8fFEbtzIj1Om8MCWLXx+6FDVBVjP6ZW8UsohFaon37gxfPYZEZ9+yo+P\nPMKV06dzqGkwK2d9zb59NkJDLUyfPpLw8LAqibk+0SSvlHJIhVeGMgaGDKHbJZfw9UMTuWbwLXgl\n9SfptwuATFatimXJknGa6B2kC3krpRyyeOdibvj0Bpr7Na9wG/uTDxOT0Jbfb3yMQV+9xeqCrfwV\n5I6fXwrBwY0r1OYNnW7g2cuerXBMNV15a9dokldKOcQmNnak7MAmFV8dasTts1mzehyNww6S/WQ2\n/3npdXqs3MVPHRpxw8J/Y/P1Oaf2licu5/PNn/PjsB8rHFNNV94kr901SimHWIyF9o3aO9RG+4Am\nrDnSgiNHOsLj6dz1nC/dbWn8N+lR2vfoD4MGwZ132idYmbOvR5Scnky+VWvag46uUUrVANOnj6RN\nm1ggE/7xhUnt2PioK3/+73PYts1eC+eOO6BzZ/j3v+Es4+vdXdzJs+ZVS+w1nSZ5pZTThYeHsWTJ\nOIYOnUVMTCxDe7/Btx3Cee54GrPz8uDhh2HLFnj7bfukqnbt7CtULV4MtjO7idwsbprkC2mfvFKq\nxtqTnc2A+HgGNW7Ms61bYznRVZOWBh9/bE/6KSn2q/xRo6C5/ebv+v3rGb1gNOvvWu/E6KuWMxYN\nUUqpStXK05PlkZEsO3aMGzZt4tiJmvT+/nDPPbBuHXz1FezfDxERcNVV8M03uNuMXskX0it5pVSN\nl2ezMX7HDpampvJNly508vY+c6esLPj8c3jnHQr+2ca7Efnc9fpqe9dOHaRX8kqpOsPdYuGN9u15\nrGVLLtmwgTn793PGxaGXF9x+OyxbxsGFH4PVBhddBDExMH8+ZGc7J3gn0ySvlKo1RgcHs7RbN17Z\nu5frNm3iYF7JXTKmUyemXuUNSUlw333wwQf2/vpx4+w3busRTfJKqVqlq48Pq3v2pKu3NxFr1vB2\ncvIZC4UXja5xd4ebboJFi+z994GB9n77qCh46y1IT3fSp6g+2ievlKq11qenc//27eTabLzUti0X\n+/sDcCznGC3/ryXHHj925kFWq33o5Tvv2Bcev+EG+0SrPn3KNdGqptCyBkqpekFEmH/wIJP37KFN\ngwbEtmpFL28PAl8IJPvJs/TDHzgA8+bZE76bmz3ZDx9ur5IJ7N6dwOTJc2tkZUxN8kqpeiXfZuOD\ngwd5NiEBP1cXNqybxvHRX+PjWo7qLSKwbJl93P3ChTBgAPuvuoZLpqxjx65pgDeQSZs2NacypiZ5\npVS9ZBNhcUoK//rlDTyDLqS7h40LPG1097ARVI5873osneCFv+A++2NsGQ2YGXErH3e7lJyQBuBv\no1HT/US2C6WZq3Ceu43zPAT3cvby+Hn40bdFX8c+YKFqSfLGmBeAa4BcYCcwSkSOF26bBIwGCoDx\nIrK4lDY0ySulKt2Ir0ewN+sYhz3bcqhBO1I9WmARKw3z9+NVkIpnQRputixcpACLFGA1bliNB7ku\n3mS5BnAg2xObXzihhw/TddtWfvC0IdmZeHvtpVWH1mS6NiLVoyVZroE0y95MWPofeFlLuAdQzE+7\nfiL7yWzcXdwd/nzVleQvB5aKiM0YMxMQEZlkjOkMzAd6A82Bn4B2JWVzTfJKqeogIuzKyWFtejq7\nsrPZlZNDSn4+mTYbOTYb3hYLfq6uNHFzo52XF1+98iVL37gdS7oP33E12+jABGYwdOgsPvwwtqjd\nhJwc3k5O5o3kZIY2bcqM8HB8S+kiavBMA9IeT6OBawOHP0+1d9cYYwYBN4rIcGPM49gT/vOF234A\npojIHyUcp0leKVXj7N6dQP/+s9m5cyoNyedPevJak448sOr1Evvkj+Tl8eiuXSxNTeW9jh2JCQg4\nYx/PGZ4cffQoXm5eDsfnjBmvo4HvC5+HAknFtu0rfE8ppWqF4pUxe8S8xLsDr2BWwSrC00vukmns\n7s6cjh35b/v23LZlC6/s3XvGrFyLsZw5U7eKnfU2hDFmCdC0+FuAAE+KyMLCfZ4E8kXk44oEMWXK\nlKLn0dHRREdHV6QZpZSqVOHhYad0zfDRxfZx9WvWQAlX6gADGjViZWQk123axLasLGa3a1dUPdNi\nLBVeQSsuLo64uLhzPs7h7hpjzEhgDHCpiOQWvnd6d80iIFa7a5RStd6DD8L27fahlpbSO0OOFxQw\nMD6ejl5evNWhAxZjaDizIYkPJtKwQUOHw6iW7hpjzADgEeDaEwm+0AJgiDHG3RgTDrQFVjtyLqWU\nqhFefNFeDmH69DJ383N1ZVFEBDuys7lv+3ZExKEr+YpytE9+NuADLDHGrDPGvA4gIpuBz4DN2Pvp\n79XLdaVUneDmBp99Zp849b//lbmrj6srC7p2ZeWxY8xMTMRgEKo3FepkKKWUqojff4frr4cVK6Bt\n2zJ3Tc7Npe+6daT8/Ry7b5tDY6/GDp9e68krpVRVuuACiI2134jNzCxz1xAPD77p0oXMsDvZklW9\nde01ySulVEXdcw9ERsLYsfb6N2WI9PXFd9/H3LEjiYwTyxhWA03ySilVUcbAf/8LmzfD7Nln3d3r\naByR3g2YuHNnNQRnp0leKaUc4elpX0x8xgx7JcsyWIyFqc0b82NKCouOHq2W8DTJK6WUo8LD7XXp\nhwyB5ORSd7MYC94Ww5yOHRnzzz+k5udXeWia5JVSqjJceSXce699ucFS1p49MU7+0oAABjVuzEPV\n0G2jSV4ppSrLpEkQFAQPPVTi5uLj5J8ND2dJSgorjpVdnthRmuSVUqqyWCz2bpvFi+1/nr652IxX\nX1dX/t22Lff+8w8FtqqbBatJXimlKlPDhvYbsRMnwvr1p2w6vazB4KAgGru58XoZ/fiO0iSvlFKV\nrUsXeO01uPFGSEkpevv0JG+M4dV27ZiekEBKFd2E1SSvlFJVYfBge9mD224DqxUoKkVwym6dvL25\nKSiIZxMSqiQMTfJKKVVVnn8ecnOhcM2M0qpQPh0WxnsHDpCQk1PpIWiSV0qpquLqCp9+Cu+/D99+\nW2qSD/bw4L7QUJ7evbvSQ9Akr5RSValJE/j8cxgzhvBD+aXWk3+kRQt+TElh81mKnZ0rTfJKKVXV\nzj8fpk/n/95KxJSSxH1dXXmweXNmVHLfvNaTV0qp6iDCNxc2ol9QbwK+WWQvbnaa9IIC2vzxB8sj\nI2nv5VVmc1pPXimlahJjeP7WlrjvSYKXXipxF19XV8aFhlbqSBtN8kopVU3yPVzZ+fbz8MILEBdX\n4j7jQkP57uhRdmVXzuIimuSVUqqaGGPIDW0GH3xgHz+/d+8Z+/i7uXFPaCjPJyZWyjk1ySulVDUp\nGkLZvz888IC9YmVu7hn7PRAaymeHD3OwlGqW53ROh1tQSilVLqeMk3/sMQgJgfHjz9gvyN2dIU2a\n8Nq+fY6f0+EWlFJKlcspSd4YmDvX3jc/Z84Z+05o3pz3DxzA6uDoQ1eHjlZKKVVuxevJA+DnB19/\nDZdcAhER0KtX0ab2Xl5s6t0blxKGWp4Lh67kjTEvGGO2GGM2GGO+NMb4Fb4fZozJMsasK3y87lCU\nSilVB5RY1qBTJ3jjDXv//JEjp2zydXX8OtzR7prFwHki0h3YDkwqtm2HiPQofNzr4HmUUqrWK612\nDTfdBLfcArfeWlSxstLO6cjBIvKTSFHEq4DmxTY79juGUkrVMaUmeYAZM8Bmg6eeqtxzVmJbo4Ef\nir1uVdhV84sx5qJKPI9SStVKJdWTL+LqCp98Ah99ZF9ZqpKctcPHGLMEaFr8LUCAJ0VkYeE+TwL5\nIvJR4T7JQEsRSTXG9AC+McZ0FpGMks4xpbDWMkB0dDTR0dEV+ChKKVWzlXklD/ZFwL/4AgYOhM6d\noWPHok1xcXHElTJLtiwOFygzxowExgCXisiZo/rt+/wCTBSRdSVs0wJlSql64coPr+ShPg9xZdsr\ny97xnXfg3/+G1avB17fEXaqlQJkxZgDwCHBt8QRvjGlsjLEUPm8NtAV2OXIupZSq7c56JX/CnXfC\nxRfDqFHg4EWwo33yswEfYMlpQyUvAeKNMeuAz4C7RCTNwXMppVStdsY4+bLMnm2vQ+/gaBuHBmGK\nSLtS3v8KqLw7B0opVQeU+0oewMMDHnnE8XM63IJSSqlyOackX1nnrNazKaVUPaZJXiml6rAyx8lX\nEU3ySilVTfRKXiml6rDKSPK7dycwbNjU8p/TobMppZQqN0eT/O7dCfTvP5v58x8u9zFaT14ppapJ\nI89GDPlyCLd+eWu59g/wDODoo0eLXk+ePJedO6cC3uU+p8NlDRylZQ2UUvWFiJR/MlQhiznZ4RIT\nE0tc3ImumvKVNdAreaWUqibGGIwDVdhDQy1AJudyJa998kopVUtMnz6SNm1isSf68tEkr5RStUR4\neBhLloxj6NBZ5T5G++SVUqoWqpZSw0oppWo2TfJKKVWHaZJXSqk6TJO8UkrVYZrklVKqDtMkr5RS\ndZgmeaWUqsM0ySulVB2mSV4ppeowTfJKKVWHOZTkjTHTjDF/GWPWG2MWGWOaFds2yRiz3RizxRhz\nheOhKqWUOleOXsm/ICLdRCQS+B8QC2CM6QwMBjoB/wJeN8ZUvL5mPREXF+fsEGoM/S5O0u/iJP0u\nzp1DSV5EMoq99AZOrGt1LfCJiBSIyB5gOxDlyLnqA/0HfJJ+Fyfpd3GSfhfnzuFFQ4wxzwAjgDQg\npvDtUGBlsd32Fb6nlFKqGp31St4Ys8QYE1/ssbHwz2sAROQpEWkJzAfGVXXASimlyq/S6skbY1oA\n/xORCGPM44CIyPOF2xYBsSLyRwnHaTF5pZSqgCpf49UY01ZEdhS+HARsLXy+AJhvjHkJezdNW2B1\nRYNUSilVMY72yc80xrTHfsM1AbgbQEQ2G2M+AzYD+cC9uvyTUkpVP6cv/6eUUqrqOHXGqzFmgDFm\nqzHmH2PMY86MxZmMMe8aYw4aY+KdHYuzGWOaG2OWGmP+LrzJ/4CzY3IWY4yHMeaPwsmGG40xsc6O\nyZmMMRZjzDpjzAJnx+Jsxpg9xSailtgVXrSvs67kjTEW4B/gMiAZWAMMEZGtZR5YBxljLgIygHki\nEuHseJypcNZ0MxHZYIzxAdYC19XHfxcAxhgvEckyxrgAK4AHRKTM/9R1lTFmAtAT8BORa50djzMZ\nY3YBPUUk9Wz7OvNKPgrYLiIJIpIPfAJc58R4nEZElgNn/cuqD0TkgIhsKHyeAWyhHs+xEJGswqce\n2O+h1cv+VWNMc2Ag8I6zY6khDOXM385M8qFAUrHXe6nH/5nVmYwxrYDuwBlDb+uLwi6K9cABYImI\nrHF2TE7yEvAI9fSHXAkEWGKMWWOMGVPWjlqFUtVIhV01XwDjTyufUa+IiK2wNlRz4PzCulD1ijHm\nKuBg4W94pvBR310oIj2w/3ZzX2GXb4mcmeT3AS2LvW5e+J6q54wxrtgT/Aci8q2z46kJROQ48Asw\nwNmxOMGFwLWF/dAfAzHGmHlOjsmpRGR/4Z+Hga8pozaYM5P8GqCtMSbMGOMODME+iaq+0iuUk+YA\nm0XkZWcH4kzGmMbGmIaFzz2B/pyccFhviMgTItJSRFpjzxNLRWSEs+NyFmOMV+FvuhhjvIErgE2l\n7e+0JC8iVuB+YDHwN/aqlVucFY8zGWM+An4H2htjEo0xo5wdk7MYYy4EhgKXFg4PW2eMqY9XrwDB\nwC/GmA3Y70v8KCLfOzkm5XxNgeWF92pWAQtFZHFpO+tkKKWUqsP0xqtSStVhmuSVUqoO0ySvlFJ1\nmCZ5pZSqwzTJK6VUHaZJXiml6jBN8kopVYdpkldKqTpMk7xSpzHG9CpckMHdGONtjNlUHwuDqbpB\nZ7wqVQJjzDTAs/CRJCLPOzkkpSpEk7xSJTDGuGEvopcNXKAL0avaSrtrlCpZY8AH8AUaODkWpSpM\nr+SVKoEx5lvstcvDgRARGefkkJSqEFdnB6BUTWOMGQ7kicgnhQvOrzDGRItInJNDU+qc6ZW8UkrV\nYdonr5RSdZgmeaWUqsM0ySulVB2mSV4ppeowTfJKKVWHaZJXSqk6TJO8UkrVYZrklVKqDvt/J4VT\nRzYk92MAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import scipy.interpolate\n", "import pylab\n", "\n", "def create_data(n):\n", " \"\"\"Given an integer n, returns n data points\n", " x and values y as a numpy.array.\"\"\"\n", " xmax = 5.\n", " x = np.linspace(0, xmax, n)\n", " y = - x**2\n", " #make x-data somewhat irregular\n", " y += 1.5 * np.random.normal(size=len(x))\n", " return x, y\n", "\n", "#main program\n", "n = 10\n", "x, y = create_data(n)\n", "\n", "#use finer and regular mesh for plot\n", "xfine = np.linspace(0.1, 4.9, n * 100)\n", "#interpolate with piecewise constant function (p=0)\n", "y0 = scipy.interpolate.interp1d(x, y, kind='nearest')\n", "#interpolate with piecewise linear func (p=1)\n", "y1 = scipy.interpolate.interp1d(x, y, kind='linear')\n", "#interpolate with piecewise constant func (p=2)\n", "y2 = scipy.interpolate.interp1d(x, y, kind='quadratic')\n", "\n", "pylab.plot(x, y, 'o', label='data point')\n", "pylab.plot(xfine, y0(xfine), label='nearest')\n", "pylab.plot(xfine, y1(xfine), label='linear')\n", "pylab.plot(xfine, y2(xfine), label='cubic')\n", "pylab.legend()\n", "pylab.xlabel('x')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Curve fitting\n", "-------------\n", "\n", "We have already seen in [the numpy chapter](14-numpy.ipynb) that we can fit polynomial functions through a data set using the `numpy.polyfit` function. Here, we introduce a more generic curve fitting algorithm.\n", "\n", "Scipy provides a somewhat generic function (based on the Levenburg-Marquardt algorithm )through `scipy.optimize.curve_fit` to fit a given (Python) function to a given data set. The assumption is that we have been given a set of data with points $x_1, x_2, …x_N$ and with corresponding function values $y_i$ and a dependence of $y_i$ on $x_i$ such that $y_i=f(x_i,\\vec{p})$. We want to determine the parameter vector $\\vec{p}=(p_1, p_2, \\ldots,\n", "p_k)$ so that $r$, the sum of the residuals, is as small as possible:\n", "\n", "$$r = \\sum\\limits_{i=1}^N \\left(y_i - f(x_i, \\vec{p})\\right)^2$$\n", "\n", "Curve fitting is of particular use if the data is noisy: for a given $x_i$ and $y_i=f(x_i,\\vec{p})$ we have a (unknown) error term $\\epsilon_i$ so that $y_i=f(x_i,\\vec{p})+\\epsilon_i$.\n", "\n", "We use the following example to clarify this:\n", "$$f(x,\\vec{p}) = a \\exp(-b x) + c, \\quad\\mathrm{i.e.}\\quad \\vec{p}=\\mathtt{a,b,c}$$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimal parameters are a=2.524, b=1.04295, and c=0.378715\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEPCAYAAACneLThAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lGW6x/HvEzoBAgQMBCRgKAKiICxNhUTArrhnsWA4\nLFhgEVFWjsfFNQQO6+6q6NrFgiAiK2BBVERKiOIuIlWUBUESWuhSQgglkOf8MZMQQiaZTKbn97mu\nuTLlLXceyP2+81RjrUVERMJTRKADEBER31GSFxEJY0ryIiJhTEleRCSMKcmLiIQxJXkRkTBWapI3\nxlQzxqwwxqw1xvxojElxsd1Lxpgtxph1xpiO3g9VRETKqnJpG1hrTxljEq21OcaYSsC/jDFfWmu/\nz9/GGHMjEG+tbWWM6QZMBrr7LmwREXGHW9U11toc59NqOC4MRUdQ9QemO7ddAUQZY2K8FaSIiHjG\nrSRvjIkwxqwF9gKLrLUri2zSBNhZ6HWm8z0REQkgd+/k86y1nYCmQDdjTDvfhiUiIt5Qap18Ydba\nLGPMUuAG4D+FPsoELi70uqnzvfMYYzRRjoiIB6y1xpP93Old08AYE+V8XgPoB2wqstk8YLBzm+7A\nEWvtPheBBv0jJSUl4DEoTsUZqjEqTu8/ysOdO/nGwLvGmAgcF4VZ1tr5xpjhjpxt33S+vskY8wtw\nHBharqhERMQr3OlC+SNwZTHvv1Hk9UNejEtERLxAI16LkZCQEOgQ3KI4vSsU4gyFGEFxBhNT3vqe\nMp3MGOvP84mIhANjDNbDhtcy9a4REQFo3rw527dvD3QYYScuLo5t27Z59Zh+v5Pfn72fhpEN/XZO\nEfE+551loMMIO67KtTx38n6vk5+8arK/TykiUmH5Pcm/uvJVTp456e/TiohUSH5P8h0bdeSfP/7T\n36cVEamQ/J7kx/QYw/PfPa/6PBERP/B7kp+W8i9On85lcfpif59aRCqwoUOHMm7cuECH4Xd+T/Iz\n33+Mw/Ob81TqU/4+tYhIqVq0aEFqamqgw/CaAIx4jeTAkpms2LGaDfs3+P/0IiIVSGCmNThbn8Y7\nu/DCdy8E5PQi4jsZGdsZNGgCiYkpDBo0gYyMsg2aKu/++dauXUvnzp2Jiori7rvv5uTJc736nn76\naVq2bEmdOnW47LLLmDt3LgCDBw9mx44d3HrrrdSpU4dJkyaVuH1I8PN0mRashWz7u8GP2bp/r2v3\nZe+zIhJaHKnjQunp22x8/BgL2QV/6/HxY2x6+ja3jlve/fOdPn3axsXF2RdffNGeOXPGfvjhh7ZK\nlSo2OTnZWmvthx9+aPfu3WuttXb27Nk2MjKy4HXz5s1tamrqeccraXtvclWuzvc9y7ue7ujRyeC8\nf7Rh84bZ8UvHe6VwRMR/XCWjpKTxhRK0LUjUSUnu/Z2Xd/9833zzjW3SpMl57/Xs2bMgyRfVsWNH\nO2/ePGutI8kvWbKkxOMX3r6wnTt32o8++sjefffd1lrHxaZv375ux+2LJO/36pqkpEksWjSKFi3i\nGN19NK+teo0TuSf8HYaI+EBmZh4QWeTdSHbvzvPL/vl2795NkybnLzMdFxdX8Hz69Ol06tSJevXq\nUa9ePTZs2MDBgwddHs/d7Tdt2kTXrl3ZvXs3AMuXL6d58+YFz5cvX16m38Mb/J7kZ8xIoUULR2G3\nbdiWLrFdeP/H9/0dhoj4QJMmETjWDSrsOLGx7qWa8u6fr3HjxmRmnr8C6Y4dOwp+Dhs2jNdee43D\nhw9z+PBh2rdvXzB2xxhzwX4lbV9Y3759effdd0lKSgJgyZIlXHfddQD06NGDHj16lOn38IaAzyf/\naPdHeX65BkeJhIOJE4cQH5/CuUR9nPj4FCZOHOKX/fP16NGDypUr8/LLL3PmzBk+/vhjvv/+e8cR\njx8nIiKCBg0akJeXx9SpU/npp58K9o2JiSE9Pb3gdWnbF7VixQquvvpqAFJTU+nbty+rVq1i7Nix\nAclzAU/y17a4lqqVqrLglwWBDkVEyqlFizgWLRpFUtIkEhNTzque9cf++apUqcLHH3/M1KlTiY6O\nZs6cOfzud78DoG3btjz66KN0796dRo0asWHDhoKkDDB27FgmTpxI/fr1ef7550vdvqjbb7+dzz//\nnJdeeokzZ85Qr149YmNjycrKuuBbgj8ExaIhM9bP4O01b5M2JM1vsYiI5zTVcPFSU1NZvHgxf/3r\nX5kwYQKtW7dm4MCBbNu2jTlz5pCUlERsbKzL/cNiquHi3H3Z3ew4uoNl25cFOhQREY9FR0fTqlUr\nZsyYQZs2bRg4cCAABw4cIDIysuLeyQO8tfotPtr4EQsGnau2ycjYTnLyNDIz82jSJIKJE4eU+Wub\niHif7uR9wxd38kGT5E+fPU3Ll1ry4Z0f0rVJVzIyttOv38ts3ToBR5cqRwOMJ/VzIuJdSvK+EZZJ\nvvDdena7ldS9MpdF9y1i0KAJvP/+/3B+n9njJCVNYsaMFL/FLCIXUpL3DV8k+YAu5H3B3fq3B6n0\naHPmr/7Sa4MiREQqsoA2vCYnTytUHQOcacDZZX9mxD//x2uDIkREKrKAZsxi79ZXj2JPtW38fszV\nXhkUISJSkQW0uubc3XqhRH/a0C6rO+9te5dFiyaSnDyJ3bvziI2NYOJENbqKiJRFQBteXfWg+eiL\n39Pn00RW3L+C+PrxfotPRNyjhlffCEjvGmNMU2A6EAPkAW9Za18qsk1v4FMgf8KHj621fynmWC57\n15y7W3f0hR+3dBx7ju3hrdve8uT3EhEfUpL3jUAl+UZAI2vtOmNMLWA10N9au6nQNr2BMdba20o5\nlst+8kX9mvMrrV5uxbo/rKNZVDO39hER/1CS942ATGtgrd1rrV3nfJ4NbASaFLOpV8frRteM5v4r\n7+eZfz3jzcOKiFQoZepdY4xpDnQEVhTzcQ9jzDpjzBfGmHZeiI1HezzKzB9nsufYHm8cTkTEKzIy\nMgIdgtvcTvLOqpoPgUecd/SFrQaaWWs7Aq8AXlnltlGtRgy+YjB/+/Zv3jiciFQgmzdvplOnTkRF\nRfHyyy/ToUMHvvnmm3IfNyMjgxUrirvPLdmOHTuYNWtWuc9fVm71rjHGVAY+B7601r7oxvYZQGdr\n7aEi79uUlHNTEiQkJJCQkFDisfYf30/bV9uyethqmtdtXmqsIuJ7oVAnf//99xMVFcVzzz13wWct\nWrRgypQpXHvttWU+7uOPP87TTz/tUUxTp06lW7dutGtXfGVHfrmmpaWRlpZW8P6ECRN8O3eNMWY6\ncNBa+6iLz2Ostfucz7sCs621zYvZzu2G18LGLR3HjqM7mHb7tDLvKyLeFwpJvl+/fgwcOJB77733\ngs/KmuRfeOEFDhw4QExMDBERETz00EMexXTq1Cn++Mc/8tprrxX7eUAaXo0xVwFJwLXGmLXGmDXG\nmBuMMcONMcOcmw0wxvxkjFkLvADc5UkwrozpMYb5W+azYf8Gbx5WRMJUnz59WLp0KSNHjqROnTps\n2bKFFi1akJqayuDBg9mxYwe33norderUYdKkSSUeKysri1mzZtG/f3+OHTtGYmKix3FVq1aN06dP\nk51dtMbbh6y1fns4TueZZ//1rL39g9s93l9EvKc8f8v+kpCQYKdMmVLwunnz5nbJkiUFz1NTU906\nzsKFC+2IESOstdb279/f5uXllSuul156yX711VfFfuaqXJ3ve5R3AzqtQVmM/M1IXlzxIit2raBb\n026BDkdESmEmlL9XtU3xXZWQdaO6acWKFbzwwgs0adKETz75hBMnTlywutO8efOoVKkSy5Yto0OH\nDixYsIAnn3ySNm3aFHvM2NhYtmzZwnXXXeeV36M0IZPka1Spwbhe4xi7ZCxLBi8JyDJaIuI+XyZo\nf+nWrRs1atRg9OjRtGvXjldfffW8z3fs2EG7du1o2bIl48aN409/+hN169alWTPHAM7ly5cD0KNH\nj4J96taty+bNm/32O4TUvL1DOw0l81gmi9MXBzoUEQlhZblJ3LhxI23btgWgcuXz74ubNWtGy5Yt\n2b9/P3Xq1KFu3brcfPPN1KhRA3Ak98IJHuDEiRNERhZdK8N3QirJV46ozMTEiTyR+kTQt+yLSPCK\niYkhPT291O32799Pw4YNCy4KMTExHD9+bp2LTZs28cMPPzB//nx69eoFwOeffw7AqlWrGDt27AW5\n6tChQzRq1Mhbv0qpQirJAwxoN4CzeWd5/es3GDRoAomJKQwaNIGMjO2BDk1EgkjRu/XCr8eOHcvE\niROpX78+zz//vMtjrFixgp49exa87t2793kDoRYuXMgXX3yBtZaTJ08yd+5cYmJiAEfde1ZW1gVx\nrF+/nquuuqpcv1tZBHyNV0+8++10Hvh4NLkvZkBeFFrkW8S/QqGffHmsWbOGt956i/r163PnnXdy\nxRVXAHD48GEmTZrEU089Veoxtm3bxpw5c0hKSiI2Nrbg/fvvv5+333672H0C0k8+GC2cnE7u4fZw\nxUfOdyLZunUCycnTAhmWiISJiIgImjZtSnR0dEGCB6hXrx7R0dH8+uuvpR7jwIEDREZGnncnv3Ll\nSvr16+eTmF0JyTv5xMQU0n65Ae64C17eDGeqF7yfmjqh3McXkZKF+518SfLy8njrrbcYPnx4mfY7\ne/YskyZN4vHHH3e5je7knZo0iYBdl8PuztAtfyodLfItIr4XERFR5gQPjjv7hx9+2AcRlSwk7+QL\nlg08PBjuvxZeW0F8zOuqkxfxk4p8J+9LAVkZypu8leTh3LKB31T/iip1TrB41FwleBE/UZL3DSX5\nYmSdyqLNK234bOBndIntApy7AGRm5tGkybl1Y0XEO5TkfUNJ3oUpa6Ywdd1Ulg1dxrZtOxxVOVsn\nAJGoe6WI9ynJ+4YaXl0Y0nEIObk5zNowi+TkaYUSPKh7pYhUZCEzQVlJKkVU4sUbXiTp4ySa7xnE\nuQSfL5Ldu/MCEZpIWIqLi9MkgT4QF+f92oawSPIA18RdQ4+Le7Cx/XJIPc75iV7dK0W8adu2bYEO\nQdwUVpnvmb7PsLPxOuI6/BHIn0TIUSc/ceKQwAUmIhIgYdHwWlhyajI/7FxPnUVXsnt3HrGx6l0j\nIqGtwveuKez46eNc+uqlzBowi54X9yx9BxGRIFfhe9cUFlk1kr/1+RuPLHiEPKvGVhGp2MIuyQPc\n0+Eeqlaqyhur3jjv/YyM7ZqDXkQqlLCrrsn30/6fSHw3kfV/WE/j2o3PzXejQVIiEmJUXVOMyy66\njGFXDuORBY8AaJCUiFRIYZvkAZ7s9SRr9qzhi81fkJmZhwZJiUhFE9ZJvkaVGrx+8+uMnD+SmKZn\nONd3Pp8GSYlIeAv7DNcvvh9XN7uaOrcdJD4+BQ2SEpGKJGwbXgvbf3w/l712GVOvncY//7FSg6RE\nJKRoMJQbpqyZwhur32D5fcupFFEpIDGIiHhCvWvccG+nex119KteD3QoIiJ+U2qSN8Y0NcakGmM2\nGGN+NMYUuxKtMeYlY8wWY8w6Y0xH74daPsYYJt88mQlfTyAzKzPQ4YiI+IU7d/JngEette2BHsBI\nY8ylhTcwxtwIxFtrWwHDgclej9QL2jZsy4NdHmTUl6O0qo2IVAilJnlr7V5r7Trn82xgI9CkyGb9\ngenObVYAUcaYGC/H6hVjrxnLz7/+zAc/fRDoUEREfK5MdfLGmOZAR2BFkY+aADsLvc7kwgtBUKhe\nuTrTb5/O6K9Gs/vY7kCHIyLiU26vDGWMqQV8CDzivKP3yPjx4wueJyQkkJCQ4OmhPNY5tjMjuozg\nvnn3Mf+e+VrGTESCSlpaGmlpaV45lltdKI0xlYHPgS+ttS8W8/lkYKm1dpbz9Sagt7V2X5HtAtaF\nsqjcs7n0mNKD4Z2H80DnBwIdjoiIS/7oQvkO8J/iErzTPGCwM5juwJGiCT7YVKlUhem/nc4TqU+Q\ncTgj0OGIiPhEqXfyxpirgG+AHwHrfDwBxAHWWvumc7tXgBtwzBsw1Fq7pphjBc2dfL7n/v0c8zbP\nY+nvlxJhKsywAREJIRrxWg5n886S+G4ivS7qzbZ/ViYzM48mTTTlgYgEDyX5cvp6/TKu/ef15L29\nDA52RguKiEgw0bQG5fTWM6nkLfo7/PYPEJGLFhQRkXChJA+OBUVWj4IT0XDNX53vakEREQl9SvJA\nkyYRQA58+g50mQzNlqEFRUQkHKhOHs5f5LvlN3DrAzT/6jZSP39cdfIiEnBqePWCjIztJCdPY/fu\nPPZ1WEJMh2osvm+RulWKSMApyXtZ7tlcek3rxYC2AxjTc0zBBUDdK0UkEJTkfWD7ke10fbsrk3u9\nwWP3fOuoyiESda8UEX9TkveRuZvmMmjGUI4/9yOcbFrok+MkJU1ixoyUgMUmIhWH+sn7yO2X3k7d\nfW2g/8M4ZnPIp+6VIhIalORLcc3JfhC1Dbq+Wuhdda8UkdCgTFWKv/7f/TRbcSX0ngCN15BfJz9x\n4pAARyYiUjrVybshI2M7v3/6MVbWXcTNex7g2fEj1egqIn6jhlc/eWLJEyzftZyFgxZSpVKVQIcj\nIhWEGl79ZGLiRGpWqcmYhWMCHYqIiFuU5MugUkQlZv7XTL7a+hXvrH0n0OGIiJRK1TUe2HRwE72m\n9mLewHl0b9o90OGISJhTdY2fXdrgUt7p/w4DZg9g97HdgQ5HRMQlJXkP3dL6FkZ0GcF/zfovTp45\nGehwRESKpeqacrDWcueHd1Krai3eue0djPHo25SISIlUXRMgxhim9p/Kmj1reH7584EOR0TkAkry\n5VSrai0+G/gZL6x4gdkbZgc6HBGR8yjJe0GzqGZ8NvAzHpr/EMu2Lwt0OCIiBZTkvaRjo468/1/v\nM2DOADYe2BjocEREACV5r+oX349n+j7DTTNvYs+xPYEOR0REvWvKq7ilAWfsmM7cn+eS9vs0aler\nHegQRSTElad3TWVvB1ORZGRsp1+/l89bGvC771JYuPAhdhzdwa3TbyX2617syTRaG1ZEAkJ38uUw\naNAE3n//f3Ak+HyOpQFTJgyi49P9yNnXC+ZNBXK0NqyIeMSn/eSNMVOMMfuMMetdfN7bGHPEGLPG\n+XjSk0BCUWZmHucneMhfGnBCygxy3l0OMRug71igJlu3TiA5eZr/AxWRCsudhtepwPWlbPONtfZK\n5+MvXogrJDRpEgEcL/KuY2nAzMw8OB0DMxZAqy+g11/Q2rAi4m+lJnlr7bfA4VI2q5Dj+SdOHEJ8\nfArnEv25pQELLgAnouG9RXD5DOjxN60NKyJ+5a2M08MYs84Y84Uxpp2Xjhn0WrSIY9GiUSQlTSIx\nMYWkpEkFde7nXQCyG8H0eVTu+TStkyoFOmwRqUC80btmNdDMWptjjLkRmAu0drXx+PHjC54nJCSQ\nkJDghRACp0WLOGbMSCn2/UWLRpGcPIndu/OIjY3ggcGfMWjJPTRr3JQhHYf4P1gRCQlpaWmkpaV5\n5Vhu9a4xxsQBn1lrL3dj2wygs7X2UDGfhVXvGk9sOriJa9+9ln9c/w/uuuyuQIcjIiHAH/3kDS7q\n3Y0xMdbafc7nXXFcOC5I8OJwaYNLWTBoAde9dx3VK1en/6X9Ax2SiISxUpO8MWYmkABEG2N2AClA\nVcBaa98EBhhjRgC5wAlAt6eluDzmcj6/53Nuev8mACV6EfEZDYYKoFW7V3HzzJt58YYXufuyuwMd\njogEKU1rEKK6xHZh8X8v5voZ15OTm8O9ne4NdEgiEmaU5AOsQ0wHlv5+Kf3e60dObg4PdX0o0CGJ\nSBhRkg8CbRq04eshX9Nneh9ycnP436v+12/nLm4WTc2tIxI+VCcfRHZl7aLv9L7c1f4uxieMv2Bh\ncG8n5OJm0dQkaiLBpzx18kryQWZf9j76vdePfpf049nrniXCOAYl+yIhlzSLZnEDvEQkMHw6C6X4\nV86Bk7RZfiNvfzWbSx67gk2/bAYgOXlaoQQPEFnuWS1LmkVTRMKDknwQyb9b//C9cWS9vJntu1rQ\n6bm+/PDz+lITckbGdgYNmkBiYgqDBk0gI2N7qecraRZNEQkT1lq/PRynE1eSksZbyLZgHQ9zxnLD\ngzZqbEN7++9Hn/8Z1kK2TUoab9PTt9n4+DGFPs+28fFjbHr6thLP5+l+IuJfztzpUd5VnXyAFNeI\neu+975CWNuGCbVsOup7sy9dTeXY/dq16naJ18snJ0zyuW8+PI38SNfWuEQk+GgwVYlytDdu+vcFR\nfXJ+su5me9L/pvsZkTeCazv9AfvLJc6E7Gh0LU/duqtZNEUkPKjyNQBcNaIac8blIiR3tL+DuQPn\nsqH1IgY950jM+XfcqlsXEVeUBQLA1Z13VlYdl4uQAFzd7Gq+HvI1f//27zzy5SPkns0FSl6hSkQq\nNtXJB0B5+6cfOXmEez66h5NnTjL7jtk0qNlAdesiYUyDoUKMNwY2nc07y5OpT/LBhg/45K5P6Nio\no09jFpHAUZIPQd668569YTYj54/klRtf0UpTImFKSb6C+2HvD9w+63buan8XT137FJUitFi4SDhR\nkhcO5hzkzjl3UjmiMu/99j1iasUEOiQR8RLNXSM0qNmAhf+9kG5NunHlm1eyJH1JuY/pyVQJIhJc\ndCcfhhanL2bwJ4O5/8r7Gdd7HJUjyj7mTdMQiwQP3cnLefpe0pc1w9fw753/ps/0PmRmZZb5GL6Y\n9VJE/E9JPkw1qtWIrwZ9xXWXXEfnNzszf8v8Mu2vaYhFwoOSfBirFFGJP/f6M3PumMPwz4fzyJeP\nkJOb49a+mipBJDzoL7YCuCbuGn74ww8cyDlApzc68d2u70rdx9OpEtRYKxJc1PBawczZMIdRX47i\nvk73Ma73OKpVruZy27IO2FJjrYhvqJ+8lMne7L0M+2wY249uZ/rt07mi0RVl2t/VguJaM1bENzSf\nvJRJo1qN+PTuT3n3h3fp+15fRncbzeNXP+5WV0tXc+EvWjRKjbUiQUh18hWUMYYhHYewethqvt7+\nNV3e7ML3md+Xul9JXSvVWCsSfPTXV8E1i2rGV4O+4rGej9H/g/6Mmj+KrFNZLrcv6W5d89qLBJ9S\nk7wxZooxZp8xZn0J27xkjNlijFlnjNGctyHGGEPS5UlseHADJ8+cpN2r7fjoPx9RXPtJSXfrLVrE\nlbjoiYj4X6kNr8aYq4FsYLq19vJiPr8ReMhae7MxphvworW2u4tjqeE1BCzbvozhnw+nZf2WvHLT\nKzSLalbwmXrQiPifz3vXGGPigM9cJPnJwFJr7Szn641AgrV2XzHbKsmHiNNnT/Psv57lH9/9g9Hd\nRzOmxxhqVKkBeG8ufBFxT6CT/GfA36y1/3a+Xgz8r7V2TTHbKsmHmIzDGTy26DFW7V7FM/2e4Y52\nd2CMR//X/MpVN0+RUBRSXSjHjx9f8DwhIYGEhAR/hyBl0KJeCz6880O+3vY1jyx4hFe+f4UXbniB\nKxtfGejQXCqpm6cSvYSCtLQ00tLSvHIsX1TXbAJ6q7om/JzNO8s7a99hXNo4bmp5E0/1eYpGtRoF\nOqwLaFCWhBt/TDVsnI/izAMGOwPpDhwpLsFL6KsUUYkHOj/AppGbqF+jPu1fa8+4peM4evJooEM7\njwZliZzjThfKmcC/gdbGmB3GmKHGmOHGmGEA1tr5QIYx5hfgDeBBn0YsARdVPYpnr3uW1cNWs+Po\nDlq93Irn/v0cJ3JPlLifvyYv06AskXM0d00F4OtGyA37N/Dk0idZtXsV43qNY2inoRdMkeDPrpfq\n5inhRhOUiUv+THgrdq3gidQn2Hl0JxMSJnBn+zupFFEJ8H89ubp5SjhRkheXAtEIuTh9MSlpKRzM\nOcjYq8eS1CGJ6/r+hbS0CRdsm5iYQmrqhe+LyDla41VcCkQjZN9L+vLt0G+ZfPNk3lv/Hq1ebsXx\ndquh0qEiW6qeXMTX9BcW5gLVCGmMIbFFIksGL2Hm72YS2SmHSo/GQfdnoMpxNHmZiH+ouibMBVMj\n5LxVn/PQrMfYW3UHrbM78/awSXRv39WvMYiEItXJS4mCrRFy86+befG7F5n500xua3Mbf+z+Rzo2\n0uSlIq4oyUtIOnTiEG+ufpOXv3+ZSxtcyqPdH+XGVjcSYVSLKFKYkryEtNNnTzN7w2yeX/48Waey\nGN55OEM7DaVBzQaBDk0kKCjJS1iw1rIicwWvr3qdTzd9yi2tb2FElxH0vLhnSMx8KeJKeQckKslL\n2Dl04hDT1k1j8qrJVK9cnT90+QP3dLiHutXrBjo0kTLxRucHJXkJW3k2j6UZS5m8ejKLti7iplY3\nMaTjEPq06FMwmlYkmHljQKIGQ0nYijAR9LmkD3PumMPWh7dy1cVX8efUPxP3QhxPLHmCnw/+HOgQ\nRUoU6FlRleQlZETXjGZk15GsfGAlCwYtIPdsLgnvJtBzSk9e+f4V9mWXf4Zrf82UKRVHoGdFVXWN\nhLQzeWdYuHUhH/z0AZ9t/ozOjTsz8LKB/Lbtb6lfo36ZjhVMA8ckfKhOXqQcCvdaiGl6ll73NWXp\ngSUs3LqQXnG9uKv9XdzS+ha3Gmy1opT4SnkHJIbUGq8i3lLcHdKq5SksWjSJd26rz7yf5zFrwywe\n/OJBelzcg99e+lv6t+nPyYOni+3OFui6UwlfLVrEBexGQUleQlZy8rRCCR4gkq1bJ5Cc7LjzTro8\niaTLk8g+nc2CXxbwyaZPeHzR45zaFcmpjJGw6U5Iiy1Y5Ptc3en5d/KaKdO7fL2IjZxPSV5Clrt3\n3rWq1mJAuwEMaDeAgYOS+eC7LtD2SxjaG07XYuuWftz/t2ReHf8k332XckHd6cSJo/z0G/mWP5Or\nq3MV9+0r/yKrRO8j1lq/PRynE/GOpKTxFrIt2EKPbJuUNN7lPgkJ4wptm2dptMZyzV9snT82tXX+\nVsf2fbuv7TriFtvjhkdsUtJ4m56+zY+/ke+kp2+z8fFjCpVXto2PH+OT36+kc3nybybWOnOnR3lX\n30MlZE2cOIT4+BTOdU8rfY7687uzGdjbCZaN5tb995P+cDpDuw6lZd86/HLtTFZ2n8mk/zzDJxs/\n4cjJI0FmcZipAAAPGklEQVTTvdKTOFxXbU3zenwlnUvtHv6n6hoJWS1axLFo0SiSkycV6rVQ8tf+\niROHuKySia4ZzT0d7uGeDveQZ/P4cd+PLEpfxBur32DwJ4PJ3V2bU3sGwfbr4V+X8913T/u9msHT\n6g5/JteSzqV2jwDw9CuAJw9UXSNBIL/aIDFxnNtVMncPetLS/AtLn7GW+3pYnoi03Nvdtnv4avvl\nli9t1sksP0TuWRVVefbzdoz+rDYKJ5SjukZJXsQN59flW0uV45YWi23c0Gts76m9beRTkfY3b/7G\njv5ytJ3902y78+hO/8ThfCQmjitxv2Cpk8//vKwX2WCSH39Cgv/iL0+SV3WNiBsuqGbIrQkZ3bm6\nZx9mDEnh5JmTfJ/5Pf/e+W9m/DiDB+c/SI3KNeh5cU96XtyTHk17cEWjK6haqap34wDcqe7wpGrL\nU6WdK5B9xt0VVr2DPL06ePJAd/LigUDcORUXQ1nuhPPy8uzmg5vttLXT7PDPhtsOr3WwNZ+qaX/z\n5m/siM9H2HfWvGN/3PejPXP2jE/jkLILxt5BlONOXtMaSFALpvlkyjs0/fjp46zdu5aVmStZudvx\n2Ju9l46NOtKpUSfHo3En2jVsV+Idvy/W7NUApXNKmt4iMzOPtLQJF+yTmJhCauqF73uL5q6RsBXu\n88kcPnGYNXvWsHbvWtbuXcu6vevIOJxBmwZt6NSoE1fEXEGHmA50uKgDDSMb+iQGTy+k4XphSExM\ncZnIY2MjSvz/6KsyKU+SV3WNBDVPGxpD2fHTx+3HKz6xXUfcYmMf6GIbPt7M1n6qtr3o2Ytsn3f7\n2Ee+fMS+vfptu3zncnvkxJFyn8+TKohwrjbytHeQL8sEXze8GmNuAF7AMf/8FGvt00U+7w18CqQ7\n3/rYWvsXj646IoVUxH7V+3Yd4LF7vmXr1g/Iv7O+JH4c7829g6PVDvPT/p/4evvXvL7qdTYd3ESd\nanVo27AtbRs4Hw3b0ia6DbG1Y91aG9eTPvSlzRsUykoaS1FSo/KgQROCskxKTfLGmAjgFaAPsBtY\naYz51Fq7qcim31hrb/NBjFKBlfQHF0w8+Zruap/iEmj61v/jtb87ksWNrW4sOEaezWNX1i42HtjI\nxoMbWb9vPbM2zGLzr5vJPp1Ny/otaR3dmtbRrWlVvxWto1sTXz+ehjUbFlwAPLmQhvPIVU97BwVr\nmbhzJ98V2GKt3Q5gjPkA6A8UTfKe1ReJ4Drh+bPrn6c86VZX0j5lSRYRJoJmUc1oFtWM61tef95n\nR08eZcuhLWz5dQubf93MwvSFvLLyFbYe2sqZvDNcUu8SLql3CQ36N6TBr7dzcPPDcORSOBpNfPO/\nlnghDfdvWJ508wzWMim14dUY8zvgemvtMOfrQUBXa+3DhbbpDXwE7AIygcestf8p5li2tPNJxRNM\nPWg8UVrjcHEXsOTkaS73AXze2Hz4xGHSD6eTfjidrYe38sOO9Sz94TuyzBFOVcuiQWQ08dHxxNWN\no3lUc+LqxnFxnYtpFtWMi6Mu5tDuI1x33Ssh+2/mC778fxwMi4asBppZa3OMMTcCc4HWXjq2hLlQ\nr98t6c7b1R17w4aVXO4zZcq9Pq+iqlejHp1rdKZzbOdzb97j+HE27yx7svew7cg2Vm5ZxbS58zhw\n+iuIyqJWbAR7TuwBoNFDjWi8ey5k1aFB1TrcffP1/Cf3J47uPUxs7Vga1GxAhPH9XWyw9PIJ1m+d\n7iT5TKBZoddNne8VsNZmF3r+pTHmNWNMfWvtoaIHGz9+fMHzhIQEEhISyhiyhJtgrct0V0lf011d\nwM6eHexyn0Ani0oRlWhapym5v55lyPC5bN36GYUvNj8sfIj6sXXZeXQnO7N2sitrF5lZmaQf28q3\nK5eReSyTzKxMjp0+RqNajWhcq/H5P2s3pnGtxsTUiiEmMoaLIi+iRpUaHsUabCNQvTWaNy0tjbS0\ntPIHhHvVNZWAn3E0vO4BvgcGWms3Ftomxlq7z/m8KzDbWtu8mGOpukYuEOp94Uv6mn7vve8U2+e6\ne/fRHDhQOairO8r773LyzEn2Zu9lz7E9jp/Zewpe78new77j+9iXvY/9x/dTtVJVLoq8iJhajqTf\nsGZDxyPy3M8GNRvQsGZDomtGU7NKTa/EGCp8Wl1jrT1rjHkIWMi5LpQbjTHDHR/bN4EBxpgRQC5w\nArjLk2CkYgqVHjSulHTn7eouPz6+HjNnDgm6r/aFlfcbVvXK1WletznN6zYvcTtrLVmnsth33JHw\n92Xv40DOAQ4cP0DG4Qy+z/y+4PWBnAP8mvMrxhiia0RzpPEZ+P3XkBMNJ6LhRH04UY9Veav5eOPH\n1K9Rn3rV61GvRj3qVq9L7aq1C3oVBUs1j69pxKsEBV8M1Q8GodyoHKx3ydZacnJz+PXEr/zh0b/y\nZdotUOM41PwVqh+GGvu5pP03XN6tOYdOHOLwicMcOnGIo6eOciL3BFHVo6hVqRZ7t5/k9NG2cKo+\nnIwkqvp/GHzX1cTFNCOqehRR1aKoU60OtavVdvysWrvgdeUI/87tqGkNRIJYqF7AQuECVdYYc8/m\ncvTUUe4b+X/MW3gHVD8J1Y9CtSyodoAOv1lA35s7cfTkUY6eOsqx08fIOpXFsVOOn1mnsjh2+hhV\nK1WldtXa1K5Wm1pVa1G7qvOn83WtKrWIrBpJraq1iKzi/Ol8XbNKTSKrRFKzSk3H86rnnleJqFLs\nADYleRHxiVC4QHkSY0nz05Q20Vj+N4ns09kcO33M8fPUsfNeHz993PEz9/i517mOnzm5ORzPdf4s\n9Pr4aceylDWr1KRGlRoFib9mlZqsGrYq4F0oRSQMhcLc7/4euGSMIbJqJJFVI4khpkznLU7htoHG\nTSx/Sr6Ti5o0ICc3hxO5J8jJzaErXT0+vu7kRaTCCZaqKHfjUHWNiEgZBUNVlLuN28Ew4lVEJKQE\nQ1WUPwYChsdsQiIiIehc20Bh3p3UTEleRCRAJk4cQnx8CucSff5AwCFeO4fq5EUkJITrCFV32gbU\n8CoiYS1YesP4W/4F4P33xyvJi0j48nSKhVC++z//wlZLvWtEJHx50gsl2KYhLqsLp6n2jBpeRSTo\nedILxfViNNN8E6SXFX9hKzsleREJep70QgmfxWjKR9U1IhL0PFktK1gX1nbX+esseE4NryISlsKh\nR45614iIlCAY5qfxBvWTFxEJY+VJ8qFROSUiIh5Rw6tIGAnlwT/iG6quEQkT4dDQKMVTdY2IhPzg\nH/ENVdeIhIlQH/wDqm7yBSV5kTAR6oN/Qn2umWAVGv/6IlIqfyxA4UuqbvIN3cmLhAlPhv4Hk3Co\nbgpGSvIiYSQYFqf2VKhXNwUrlZ6IBIVQr24KVm71kzfG3AC8gOOiMMVa+3Qx27wE3IjjX2iItXZd\nMduon7yIuBQuc814m0/7yRtjIoBXgOuB9sBAY8ylRba5EYi31rYChgOTPQkmWKSlpQU6BLcoTu8K\nhThDIUbwPM786qbU1AnMmJHi8wQfKuVZHu5U13QFtlhrt1trc4EPgP5FtukPTAew1q4AoowxMV6N\n1I9C5R9ecXpXKMQZCjGC4gwm7iT5JsDOQq93Od8raZvMYrYRERE/U8OriEgYK7Xh1RjTHRhvrb3B\n+fpPgC3c+GqMmQwstdbOcr7eBPS21u4rciy1uoqIeMDThld3+smvBFoaY+KAPcDdwMAi28wDRgKz\nnBeFI0UTfHmCFBERz5Sa5K21Z40xDwELOdeFcqMxZrjjY/umtXa+MeYmY8wvOLpQDvVt2CIi4g6/\nzicvIiL+5ZOGV2PMDcaYTcaYzcaYx11s85IxZosxZp0xpqMv4ihNaXEaY3obY44YY9Y4H08GIMYp\nxph9xpj1JWwTDGVZYpxBUpZNjTGpxpgNxpgfjTEPu9guoOXpTpxBUp7VjDErjDFrnXEWO59CEJRn\nqXEGQ3k644hwnn+ei8/LXpbWWq8+cFw4fgHigCrAOuDSItvcCHzhfN4N+M7bcXgpzt7APH/HViSG\nq4GOwHoXnwe8LN2MMxjKshHQ0fm8FvBzkP7fdCfOgJenM46azp+VgO+ArsFWnm7GGSzl+UdgRnGx\neFqWvriTD5XBU+7ECRDQxmJr7bfA4RI2CYaydCdOCHxZ7rXO6TastdnARi4czxHw8nQzTghweQJY\na3OcT6vhaOMrWv8b8PJ0nru0OCHA5WmMaQrcBLztYhOPytIXST5UBk+5EydAD+dXoy+MMe38E1qZ\nBENZuitoytIY0xzHN48VRT4KqvIsIU4IgvJ0Vi+sBfYCi6y1K4tsEhTl6UacEPjy/AfwGMVfgMDD\nstRgqJKtBppZazvimL9nboDjCWVBU5bGmFrAh8AjzjvloFRKnEFRntbaPGttJ6Ap0C3QF29X3Igz\noOVpjLkZ2Of8Bmfw4rcKXyT5TKBZoddNne8V3ebiUrbxtVLjtNZm53/Ns9Z+CVQxxtT3X4huCYay\nLFWwlKUxpjKOxPmetfbTYjYJivIsLc5gKc9C8WQBS4EbinwUFOWZz1WcQVCeVwG3GWPSgX8CicaY\n6UW28agsfZHkCwZPGWOq4hg8VbSleB4wGApG1BY7eMrHSo2zcH2XMaYrji6nh/wbpuP0uL6yB0NZ\n5nMZZxCV5TvAf6y1L7r4PFjKs8Q4g6E8jTENjDFRzuc1gH7ApiKbBbw83Ykz0OVprX3CWtvMWnsJ\njlyUaq0dXGQzj8rS6ytD2RAZPOVOnMAAY8wIIBc4Adzl7ziNMTOBBCDaGLMDSAGqEkRl6U6cBEdZ\nXgUkAT8662ct8ASOHlZBU57uxEkQlCfQGHjXOKYjjwBmOcsvqP7W3YmT4CjPC3ijLDUYSkQkjKnh\nVUQkjCnJi4iEMSV5EZEwpiQvIhLGlORFRMKYkryISBhTkhcRCWNK8iIiYUxJXio0Y0wXY8wPxpiq\nxphIY8xPwTrJlognNOJVKjxjzP8BNZyPndbapwMckojXKMlLhWeMqYJjwroTQE+rPwoJI6quEYEG\nOJbZqw1UD3AsIl6lO3mp8Iwxn+KYw7sFEGutHRXgkES8xutTDYuEEmPMfwOnrbUfOKei/ZcxJsFa\nmxbg0ES8QnfyIiJhTHXyIiJhTEleRCSMKcmLiIQxJXkRkTCmJC8iEsaU5EVEwpiSvIhIGFOSFxEJ\nY/8Pi+zxhMGjLI0AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "from scipy.optimize import curve_fit\n", "\n", "\n", "def f(x, a, b, c):\n", " \"\"\"Fit function y=f(x,p) with parameters p=(a,b,c). \"\"\"\n", " return a * np.exp(- b * x) + c\n", "\n", "#create fake data\n", "x = np.linspace(0, 4, 50)\n", "y = f(x, a=2.5, b=1.3, c=0.5)\n", "#add noise\n", "yi = y + 0.2 * np.random.normal(size=len(x))\n", "\n", "#call curve fit function\n", "popt, pcov = curve_fit(f, x, yi)\n", "a, b, c = popt\n", "print(\"Optimal parameters are a=%g, b=%g, and c=%g\" % (a, b, c))\n", "\n", "#plotting\n", "import pylab\n", "yfitted = f(x, *popt) # equivalent to f(x, popt[0], popt[1], popt[2])\n", "pylab.plot(x, yi, 'o', label='data $y_i$')\n", "pylab.plot(x, yfitted, '-', label='fit $f(x_i)$')\n", "pylab.xlabel('x')\n", "pylab.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that in the source code above we define the fitting function $y = f(x)$ through Python code. We can thus fit (nearly) arbitrary functions using the `curve_fit` method.\n", "\n", "The `curve_fit` function returns a tuple `popt, pcov`. The first entry `popt` contains a tuple of the OPTimal Parameters (in the sense that these minimise equation (\\[eq:1\\]). The second entry contains the covariance matrix for all parameters. The diagonals provide the variance of the parameter estimations.\n", "\n", "For the curve fitting process to work, the Levenburg-Marquardt algorithm needs to start the fitting process with initial guesses for the final parameters. If these are not specified (as in the example above), the value “1.0“ is used for the initial guess.\n", "\n", "If the algorithm fails to fit a function to data (even though the function describes the data reasonably), we need to give the algorithm better estimates for the initial parameters. For the example shown above, we could give the estimates to the `curve_fit` function by changing the line\n", "\n", "```python\n", "popt, pcov = curve_fit(f, x, yi)\n", "```\n", "\n", "to\n", "\n", "```python\n", "popt, pcov = curve_fit(f, x, yi, p0=(2,1,0.6))\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "if our initial guesses would be *a* = 2, *b* = 1 and *c* = 0.6. Once we take the algorithm “roughly in the right area” in parameter space, the fitting usually works well.\n", "\n", "Fourier transforms\n", "------------------\n", "\n", "In the next example, we create a signal as a superposition of a 50 Hz and 70 Hz sine wave (with a slight phase shift between them). We then Fourier transform the signal and plot the absolute value of the (complex) discrete Fourier transform coefficients against frequency, and expect to see peaks at 50Hz and 70Hz." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEPCAYAAABoekJnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYHVWZ/z9v9r3TAQIhMSEsssmqIopKEAVkFFBUVFQW\nEUfZHBxHXOYHjI4z6oiKgtsgCoOigCKCyCIEBGUNYQ0hhDUhCZCts3XSy/v741R5q6ur7q26td3u\nfj/Pc5++S9065377nPOe9yzvEVXFMAzDMJplWNUZMAzDMAY2ZkgMwzCMTJghMQzDMDJhhsQwDMPI\nhBkSwzAMIxNmSAzDMIxMVGpIRGSGiNwmIo+LyKMicmbMdReKyCIRmS8i+5adT8MwDCOeERWn3w2c\nrarzRWQC8KCI3KyqT/oXiMi7gZ1UdRcReRPwY+DAivJrGIZhhKjUI1HV5ao633u+HlgATA9ddjRw\nmXfNvUCbiGxbakYNwzCMWFpmjkREdgD2Be4NfTQdeDHwein9jY1hGIZRES1hSLxhrauBszzPxDAM\nwxggVD1HgoiMwBmRy1X1DxGXLAVeE3g9w3sv6l4WOMwwDCMlqipZvt8KHsnPgSdU9fsxn18HfAJA\nRA4E1qjqiribqWppj7vvVq6/vrz00jzOPffcUtO7/Xbllluq/91V66Cq/PnPyh13VP/bW0GL665T\n/v736n97K2hx9dXKvHnV//bwIw8q9UhE5CDgeOBREXkIUODLwCxAVfWnqvonETlSRJ4GNgAnVZfj\nvtx0E1x1FRx5JEgmez7wueEGmDsX3vnOqnNSPX/4Azz+ONxxR9U5qZ5rroEVK+DGG6vOSfX85jfQ\n0+M0GWxUakhU9W5geILrTi8hO6np6IAFC+Bvf4ODDqo6N9XS0QEPPAAPPQT77Vd1bqqlowPuvBMW\nLoRdd606N9XS0eE6XM8/D7NmVZ2baunogL/8xRnWbQfZutNWGNoasKxdC7vvDj/7WdU56c+cOXNK\nTa9VtShbBzAtgqxdC7vtBpdcUnrSdalKi112gV/8ovSkC8cMSQY6OuCkk1wvowzWrYOnn052bdkV\npaMDTjyxPC3WrIFnnml8XRUNRtlarFzpevyNGApavPwyLFnS+LqqtCizvVi2zD3KwAxJBtauhR12\ncA18GVx7LRx/fDlppaVsLX7zG/jkJ8tJKy1la3HZZXDaaeWklZaytfjZz+Dss8tJKy1la/GDH8CX\nvlROWmZIMrB2LWy/PaxfDzktfqjLypVw//1ujLXVCGpRBitXwl13Oc+k1ahCi7/8BTZuLCe9NFSh\nxU03QVdXOemloQot/vQn6O0tPq3KDYmIXCIiK0TkkZjPDxaRNSIyz3t8tew8xtHRAVttBSNGwObN\nxae3apUzWK24Aqajw1WSDRvKMaqrVkF3N9x8c/FppaWjA6ZPL6/BWLUKOjvh9tvLSS8NVWjR0eE6\nGa1ET48z9NOmlavFK6+4zmfRVG5IgEuBwxtcc6eq7u89vl7vwtWr88tYI9auhUmTYOLEctzVVatg\n//3hj39sfO2LL7r8lcXatTBlCoweXU7P2Nfi+usbX/v88+VVXlWnxbRpTocyeoNpysWzz5bnuWzZ\n4jyDqVPLbTyTarF4sTPAZbBuHUyY4NqLsoa20miRlcoNiareBTRq/hPv0vinfyrPre3ogLY2V0DK\nqCirV8PHPga33NK413/OOfC+97meUNH09rrfP3Fi+VrcdFPjaz/3OTjuuHIa9U2bYORIGDMGxo1z\nHlrRpNHin/8ZTjihHK/Rrx/jxjmPvbu7+DTTaHHCCfDpT5erRVn1A9JpkZXKDUlC3uydRXKDiOxR\n78KnnoLzzis+Q11drnKMG1de4Vi1yu1L6O52DVY9Vq+GRx6Bb36z+HytXw/jx8Pw4eVqsc8+8Oqr\njQ3E6tVw771w4YXF52vtWtdgQLla7LefW7HUiNWr3RBYGUuTfS1EXPkow6im1eKGG+CKK4rPlz96\nMXq0q79ldHbTaJGVymNtJeBBYKaqbvTOJrkWeG3cxXvtdR6/+IXrFc6ZM6ewZX4dHa5giLgGo6yh\nrSlTag3UuHHx165dC2ec4byXL3+52Hz5lQTKHebbemsYO9YN1UyYUD9/Z57ptPjc54rNV5QW06YV\nm+aqVW4eorPTeaDD62zxDWpx6qnF5iusxfr1NSNbFKtWwcyZyQx4UIuPfazYfAWNql9/29uLTTNO\ni7lz5zJ37txc02p5Q6KBaMCqeqOIXCwiU1R1VdT1J598Hl/7WvFeiW9IoFZJiiZsSKZOjb92zRqY\nPdstGS4a322HcnvhQS3qGRJfizJc/Kq02Gor17HYuNGVxzh8Lf761+LzFdai6A6GqtNi++3d/Ex3\nt1sIE4evxYMPFpsviG4vijQkXV2uLEStEgt3sM8///zM6bXK0JYQMw8SPMRKRA4AJM6IgFunXcaS\n0KqGMNrbk1VKf816WVr4laQKLRqlN5i16OlxZcEffx/K5WLjRueNjRnTOD1/eHjGjMHZXqxZA5Mn\nu6G03l5nWIukco9ERH4FzAG2EpEXgHOBUXhBG4EPiMhngC5gE3BcvfvNmuVEVC02kGLZva3eXve7\nkjaea9bUtCiaoBZlDG1t2uT08Oen6qWn6vI3c+bg1GLNGtdYJ5mf6ux0emy3XTVaFN14+l4q1LSY\nPDn6Wt/ITZlSvhZltBe+FsGhNF+bIqjckKjqRxt8fhFwUdL7TZ3q3NlNm+rPIWQlavy3SDo6XIEY\nMSJ5b2v6dFdge3thWIG+Z9m9rdWr+1eSONavdz3UrbcenD1PX4sk6a1d6xrWyZOr0aLoxrPVtSiz\nvYgyqkUaklYZ2sqNMWPKKRxlNxjBgtGoIPrjsSNGuNUyRVfgsocw0mjhNxjjxjn3vmgXvwot/LH2\nJFq0tbmH77UXSauXi7a26oxqK2mRB4POkEA5hSM4eVamq+qnV69grFlTK7RlaVHmcE4zWog4LYre\npDkQtBg1yj2K3pjYCkNbcfhaTJzoliUXvd+q7PYijXeWB2ZImiTYw6iiktQriH4vHMrToqqeZ5Ih\njDKNatVaWLmopZdEi2HDXP3t6Cg2b63cXuSBGZImCfcwyhzCSOuRFB02puwlr2m1CDaeg1GLZj3V\nsrUowztr5TrSqu1FHlRuSBoFbfSuuVBEFnm72/dtdM8q5khayVWtoudZ5nBOWi3K9kjK1iJpg1F1\nuWjFhQcwONuLNB2MPKjckNAgaKO3m30nVd0F+DTw40Y3zFIwFi92gc4axQWqchVGo/TymiN59FF4\n85sbhyCpelK1XqUMeyTNanHffXDIIY0nqKvWooxycccdcOSRja8bClr8+c/wgQ80vq7s9iJoVIfE\nZHuCoI1HA5d5194LtAU3KUaRpWAsXOjOHb/yyvrXtfIQRrC31d7evBZPPgn33NN4d3yra5FHg7Fg\nAcyd2zhsfatrkYdRfeIJd5RBo93xVe2d8NMrQ4vHH4drroEHHqh/XSuXizyo3JAkYDrwYuD1Uu+9\nWLIUjCVLYOed4RvfqN8Tb2VXNa/G09fiP/+zfk+87FhbVWrx9bqHGAwtLf7zP+tf1wp7J+rlLU8t\nvvGN+te1cnuRB5VvSMyb8847jwcfhJdegrlz0wdtXLIEPvpRFxH0iSfgda+Lvs4/XwDK23jm96Aa\nFcQ1a1yMHXDfee655tJcssQdZ3vhhe4cix13jL5u3bpafKcqtGg0hDFrlnuetcE4/XQ4/3xYvtzt\nDo+iai0alYu8Gs9//Vd3pG1wEjlM1VrUO8c+Ty2+/GUXjn7LFresOoqq24vgqapDMmgjzgN5TeD1\nDO+9SM477zx++1u4+mpoJvDv0qVuXmDatPorOTZudJv9oJzeln/eByTrbe3hBdvPUkmWLoXXv941\nmkm1KKOSpNUi2GAsjS059Vm6FI44wkVOWL063pC0uhbB4Zxmj2xeutR1KrbaymkRZ0jCWhTdC8+i\nRZY68trXOg3WroVttul/TU+PMzJjx7rXVbQXixfXPhuSQRuB64BPAIjIgcAaVa1bBbL2MGbMcHML\n9RrPTZtqBSNJJfnWt9zcS7MEI9yWtSGxGS2SDOd8/etubLlZ0mqRR4NRlBb//u+waFFzeYLWLRd+\noMAxY9zrJI3nOefU9yIa0apadHY6HfzYf0nai7PPhmXLmssTpNMiDyo3JF7Qxr8BrxWRF0TkJBH5\ntIicCqCqfwKeFZGngZ8An210zzwKRqN7bNxYi+U1fnzjY1UvvxxOPrn5U+I2bKgVjKRhQaB8LZIU\n2ksvhVNOaf7EwrAWjTae5dlg1LtHb6877MxvPJNo8bOfuVMLmw1XUna5UE2mxaZN/RvPRlpcfLEb\nPsyiRdJRgjy06Olxw5zTptW/R7B+QGMtenvhBz9o/uwc1XRa5EHlhkRVP6qq26vqaFWdqaqXqupP\nvMi//jWnq+rOqrqPqs5rdM8yGs9gz9MPXV3v1MLly12l+slPmsvX+vXJh0zy6G319rp5pu23r38P\n1VqPC1yF2bQpPuSEqutpbd4Mv/xl+nxB+Vp0dtaGLerdo7PThe32A2Q2yltPD6xc6R5XXZU+X1C+\nFmvXuvI+cWJjQ+LXjyR527DBnaGxeLE7tTAtquVrsWKFG94bNSqdFo0a9lWrXD164AG47bb0+dqy\nxZXBkSPd6yHhkRRBswWjo8NV7ra2+q6qanRFiesZd3W5Cvj5z8Odd6bPF/R3VYsOhfHyy+67Y8bU\n12LzZleR/MZz2LD6x6quW+euOe20/LQouuf50kuu1zlsWH0twmVi7NjaAUtRvPKKu98nP+n2ZqSl\nt7dvz7OMcrFkiYsqDem0aOQ5rljh5p0+8YnmtOjsdOXQP8iqnhaq+XiqfqcT0mnR6P+0fLm77/HH\nN1dHgvUjSXp5MCgNSbPRTZcudf9AP8BfXOHavNlZ++CRpvUatBUrXG921ix48cXoaxoRLBx+Qx33\n+/KuJGncdqivhT8U8JrXNKdFT0/fIwLKWObZrBaNwtwvW5ZNC3/4yC+H9dLq7XWNiT8xXrYWo0bV\nP2ApqxZRjWecFp2drlPge9Fla+HnLa7+lqlFXgxKQzJ6dO1MkjSEC0bSHgbUd1f95aJZGs/gqo8R\nI1zFjPt9wfX7kya5fKWdj2i2twWNDUkWLfxK6XtA/vGyUb+vp6fvee55NRitokWaBmPDhr5Gp2wt\nfKMa56mWqUWwfkDxhiSshe/BxxnVMrXIi8oNiYgcISJPishTIvLFiM8PFpE1IjLPe3w1yX2bCcQW\nbjzT9MLHj2/cYGy/vfNO0k64+8MXwRMf4wpHT48bSvN7W8OGuWvThk9P09vK0nim9RrDlaTe/JT/\nf/J1GzfOabN5c7o0k5aLqMPUkmgxY0Y+DUa9+angUlyo/U/T6p+ljhRtVIO/r15aUVo0E7QxaWer\n2fbCDElCRGQY8ENcrK09gY+IyG4Rl96pqvt7jwZ7ix3N7CoOjv+m6WGAex3nIfgFY+RIN8SVdllf\nuJJAvAfkN2ZBo5NVi7SNZxItJkxwnuOqVenyFZwT8GmkhY9IseUiyqgm0WKbbVz+054PEtbCn5+K\nuk9Yi5Ejm/fai6wj06e752nPBwmuXvPTipufCmvhe0rNeO1FatFsZytcLobCqq0DgEWq+ryqdgFX\n4mJrhUl9+vqkSekbDH/8HtIP5/hDLHH39TewNdPLCPcwIL6XEdWYZdUi7XBOFVpE/b4itCiiXAwb\n5hqktJslqy4XRWgxapRbCbV8ebp8hbWoN5QW1mL4cJe3uGG3OMJapDEkSbSYNMnlLe2wW1iLRvNT\neVC1IQnH0VpCdBytN3sh5G8QkT2S3LiZw2peftntXIb0E8xJG89mhjHSNp7hvBWtRZpKsmxZMVrE\nNZ5la9FsuSjDqA4ULaqqI2mNaliLNENbZWrRaNFHHlRtSJLwIDBTVffFDYNdm+RLzfS2VqyAbb24\nwmldVX98OoqsDUbYbYd0jWczWrz8ck2LtENbRWqR1ZBk1aKVyoVpUSNq+DdOi6gyO2lSeqPaqu1F\nmnKRF1XH2loKzAy87hdHS1XXB57fKCIXi8gUVY0cXT/vvPMAeOYZuPfeORxzzJzEmQn2MIJnOQeX\n+UJ0pRw7NnnPM20YiDS9rahK0kzPc8WK5L2tqPHfVuiFF6FFo+GcLOWiaC2iVhqm0UI1uRZZ60ir\neyTd3c5wbLWVe5124UGVWpQetFFE1gFRUz0CqKrGhGtLzP3AziIyC1gGfBj4SCgP2/qxtUTkAEDi\njAjUDMkrr9RWVCQl2MMYNqwWiM0Px+yTdV7grrvS5SuqYMT1aPIYC1fta1T9+wV3sPvEeSRJtbjl\nluT5gmjvrEgturpcY+s3GG1trkz09taWINdLL80QxsMPJ88XpNciq0eyfr3rVPk9/2Z64VFa9Pb2\nrXt5ee1py0Uao/rqq65dCC6nzmO+aPNm9z/x25witCg9aKOqTlTVSRGPiTkYEVS1BzgduBl4HLhS\nVRcEY20BHxCRx0TkIeB7wHFJ7p22t7Vli6so/s5fiC8cVUwwh932uPTyGAtfu9atqAr+xjy06Olx\nFdA3UANBi1decUbENxojRrh7Ju3112vMNm+ulbeBoEXQGwHXWHV2OmMbJk25WL26tooP8uuFp9Ui\njVENdrSg/nLqNFr49/XLW9HlIi9SDW2JyFTgH31SVX0hawZU9c/ArqH3fhJ4fhFwUdr7pu1tvfyy\nW4YZ7GXGuatxk2dRS1n9Hax+SOe8ehhxywfz6HkGx8F9fC38VSrB9JI2nq++6iqcHwMorwajKi38\nHfM+abwzvwfuL9MeiFoEI0CEw6fH1ZGovIXPd2lWi2DDDvFa5DFHEjaqwU3Q4Xtv3FjzaH2K1mLm\nzL7v1VtunAeJJttF5CgRWQQ8C9wBPAfcWFy2spO2txXVYMS57nHrwqMajJUr+xaiadNcjz98bb1N\niml6W3nMC4QrCdTXIun4b1iLmTPdktfwbx+oWqSZLwprMXu2Ozws3KPNU4uscyRpy0VWLcI00iKr\nd5a1g5HGay9ai6TlIi+Srtr6GnAg8JSqzgYOBe4pLFc5kLa3FVdJogpGmuV8wfhG4DyeHXZwiwF8\nVOGtb40/DzxtzzPrvECaSpJmXiCsxejRrvf1QsCv7e2FN7whfh6plbVIs4ItrEVbm5t/evnl2nvd\n3bDXXvDgg9F5q9ojgfzqiO+xgysTGzb0zdvmzbDrru7U0ijyKBdFdjCabS9mznRL5oN7QDZsgJ12\n6tuGBEmjRV4kNSRdqroSGCYiw1T1duANxWUrO3l4JHFDW2nGwjs6+lYScIUgeGLZbbfBvffC3/8e\nnbesDUYzWoQrSR5ahBsM6K/FDTe4Sec4LbIO8xWpRRqjGlcunn669vrqq+HJJ+GemC5bK2uRto4E\nG08RdwJjsFxccYVrOO+7LzpvabSI81TzGvKMSq/Z9mLkSLcQI3hc9qWXus7XAw9E5y2NFnmR1JCs\nEZEJwJ3AFSLyfSDlPtByycsjSVMwkvQwoH/j+T//A4cfDvPnR+ctjdseN/6bVos0w3xJe1tR53sP\nNi3yKheq8O1v56tFVu+sSC3qdTB6e5OViyyT7Xl5JFkXpDQqFz09cMEF+ZWLvEhqSI4GNgH/AvwZ\nWAy8N48MNAra6F1zoYgs8na375vkvnn0tvxw9GHSrAtv5JG89JLrZX372/EFI23PM+tYeFotko7/\nNmowFi+Gp55yR/GmaTBaRYs080WNysWjj7rFG1/9an5aFOGR5FVH6jWe993nDOvZZ6drPNOWiyyr\ntqC2NDwqvbzKxR13uKXBn/lMPuUiLxIZElXdoKo9qtqtqr9U1Qu9oa5MJAnaKCLvBnZS1V2ATwM/\nTnLvPHpbcbF68uxh3H03HHQQ7LGHW2oaVSnT7iMpohdeT4u8PJK//hUOPhj23tsNY0T9vjx6nkVp\nkWYFW6NycdddcMghsM8+bl4ganK1lbXI0yO56y449FDYd1837Bm1xLYKj6SK9iKoRR7eWV4kXbX1\nfs8jWCsiHSKyTkRS7g+O5AAaB208GrgMQFXvBdpEJPQv7E9cb2vzZtezaW/vO1kV1cNIE26i2TmS\nv/3NGZLhw10D+sgj/e8R18PIul9g0yY44wynxUsv1d5Pq0VecyS+FqNGuYnVxx/vf4+ieuHr1rmz\n06dMcStnfNJokdaoJtFi4kQX0PGpp/rfoyjvbPVqd2rj1lv3/Z1F1pF6jaevxdSp7h4vRGw6SKNF\nmjmSV191pzZuu23f4weqbi9mznTfDy7O8GlZjwT4FnCUqrbluSGRZEEbw9csjbimH3G9rTvugFtv\nhQMP7HuMZdTkWb0GI81wTriSzJ7t1oZ3dzuP5C1vce/vs090LyNrJYnT4qab3ETufvv1XSWVtvFM\nqkW9BkO1Wi3++Ed47DHYffe+E/1FaZHEUw1qEbXrvSgtrr7ahfGZNavv5HaRWsQ1nkWUizTe2RVX\nuCHGbbaBhx5y74UjP/gUXS56elx9ffOb3YKEPMpFXiTdkLhCVRcUl4388EOkqML69XPo7Z3TZ5Ph\n4sXOiOyzjyugJ57oJvOeecZVnCBpCkY9V3Xrrfu+N3q0K4QLF7pe9xu89W977x291DPqDI56bnvS\nnufixa53s/32TosPfcjtUn7xRbcRKkgevfB16/qHrfGXvS5c6Hqa++zj3t97bzdPEKZILd7+ducN\n3X03vOc97vsrV/bfhFmUdzZtmrvvwoUuj7t623R9T/UjH+l7fVotknpnixfDnDnus7vvhne8w3kp\nXV39wwVNmODO5QiTdThn1iznJT/5pFu55G+w87U4OjRuERcWJI9yceih7u/dd7u246WXXH7Dek6Y\n4DyYMFmH+Xbc0e0leewx19n1N3/6WrzrXbVre3sb18nSY20FeEBEfoOLvPsPB09Vf5cx/YZBG73X\nr2lwzT/wDQm41Q0bNvT9xyxe7Cz8W94CF3n75Z97zhWM8O7TemPhaYYwdtyx//t77w0nneT++oVs\nxgy4/vr+12btbY0bVzvkZ0TgP754seuB778/nHWWe2/hQmdEwg2Uf058VHpplryGGwxwGnziE3DA\nAbX8zZjhlkSHKarnuXgxvO1t7rd/3Ts67bHHXGPu78T3KWq+SMTtG/n4x1359DtAM2Y4TzpMkVq8\n//3uXj/2ZiQfftj9n8LxxdIO5yQd5hs5El77WjjhBNfZ8SMAzJjRfy9JT098rz+rd7Z4Mbzzna4B\nv+Ya+PznnUe0b8SSnwkTogOyptUiXC4mTHCdjE99ymnhM2NG/0PyNm1yHbPw/ymoRemxtgJMAjYC\nh+FWa70XeE/m1ANBG0VkFC5o43Wha64DPgEgIgcCa/wgjo2I6mX4hmSvvdyu6lWr4gtG3HGYWfdO\nAPz61/DBD8KZZ9be23ZbN4kXJutku38eQbii+Fq8/vWu57d+ff1KUpQWv/sdvO99cNpptffSapG0\nFz56tOu1hY/b9bU48ECYN88Z3rzKRdKeJ8CNN8J73wuf/nTtvSgturtdHpOmV68XHp68Dna27rnH\nNdRF15GoDsZtt8G73+3ma3yitPD/z+HGMy/vbKedXAN+991Oqyrai7vvhsMOc4bVJ0qLqPrhp1d5\nrC1VPamIxFW1R0T8oI3DgEv8oI3uY/2pqv5JRI4Ukadxe1cS5yWql+EXjBEjXA/4b39L33jWW86n\n2veY27he+MSJ8IUv9H1v6tR0k2dJ905ATYv29tp7vhZjxrghpfvuy1eLMHFaTJ4MX/pS3/eitNiy\nxTVqfnC/YHpJe54iNS2C9/G1mDQJdt7ZGZMqtJgyBf793/u+F6WFP6wlobND02gxcqQbygt+puq0\n2Hlnl5epU93w6/z5LvpCmLy0iGo8t9kGwp3lKC3iGs803lnwuF3fIPX0uNGKHXd0dWTECDcEPn++\n6/iEidJC1QW2zDKPCG63v+8p++ShRV4kMiQicmHE22uBB1T1D1ky0Choo/f69GbuHe5lqLqCsNNO\n7vWhh8K117plt5/4RP/vp+mFDx/uKubmzX1Drcf1MKLwC0bQGPk951Gj+l6bZnUO9Neiu9vNheyw\ng3vta7FggVvVFqaZicSwUU2jRb3eVtLGs5EW/tzVxo1uDsA/f9vX4uGH4bjj+n+/SO8siigt1q3L\n3nhCTQv/s5UrXUPqz4X4WsyfD6dH1MI0jWeaOZI48tAirv6OHdt3KHzpUjfc7V8bLBdRo0FRWnR2\nurqb1FtKWy7ChiStFnmRdI5kDLAbcJX3+lhcAMd9ROQQVf1cEZnLStgjWbHCCeoX2lNOgV12cT2N\n73yn//fTNBhQKxxBQxLXw4hi3DhnjDo6apFl167tG9o+eG2aBiOsxQsvuF6O3yv/9KfdcJ9qbcI7\nSD23PZzeiBHxRjWpFu3t7votW2pGtJ4WSYcwoL8WzzzjDKpf2T/7WXjTm1z+9967//fjGs96mzOT\neqpRbLNN/w5GM1pElVlfCz/irO+Z+Zx+umtAOzpgzz37fz9qvmjzZvf/T9p4xnkkUUT1wvMuF35e\nwlqccYab5F+zxs3fhIkqF43aijBpysXUqf2Nalot8iLpHMnewCGq+gNV/QHwTpxheR9u3qQlCffC\nwwVj663dSpjOzugJ8TST7RDduKfpYUD/Xsbq1X2Ho3yCDVSSvDXSYvp0OOoo1wCEVylBtBZxjSdE\na5GmwRg2zP1/Xnml9l49LZrphfuEtdhpJ7cZcKutotOL0mLLFqdd+DTNESPcIxh0D9KVizFj3G8M\nblYtS4s993RzaLNnR/+foxrPuLT8vAXLrGo6Ldrb3f07O2vvpdHCnx8LH9AGjbV4wxucAdljj76L\nVnyiykWatqKryz2i8hbF1KmufvT21t5LWy7yIqlH0g5MwA1nAYwHpnhzHJvjv1Yt4Z5nuGAA/Nu/\nueWv4d4TuJ5wb2/fXrFqfC8jatwzTQ8Dar2MXXZxr9esie5hjBjhGq2urr7DXo3mSHyitPjKV1wl\nCQ8dQc0jCfaKu7qcblGVytciWKjTeCRQM6r+kFOcFs3OF/lEaXHuuW6/URRR3llcmfDzt2lTzftL\n23hCTQvHKiQnAAAgAElEQVRfzzgtonqe9RrPJFp8/evxQTTTaDFihCsvwTK7ebMrT+F5rzhEnIf2\nyiu1JepptPDzFlXGk2jxX//lVvNFkUaLKE/Vrx9ReYti1ChnvFavrq04TaNFniQ1JN8C5ovIXNwx\nu28HviEi44GY6lY9/lG5PlEFY4cdXDyjOPxehl/wu7rcPzq8JBSi/1l5eCRRBQNqhTFoSOoNYTTS\nYpddnGGNwh/n3bKlVunjelvQX4veXvc6vKy4HmHXPU6LkSPd/6Srq/Z/6e2NPhoYorXYdde+17zu\nde4RRZRRTaKFn/eNG2sHISXF18LPZ5wWUT1PX4eozlKUFm97W99r9tvPPaJI45FATQu/zKatH1Cb\nJ/ENSRot6hn8KC3Ck+oHHugeUaTRIlhms2rx8ss1Q5JGizxJGmvrEuAtuH0kvwfeqqr/68Xg+kL9\nb0cjIu0icrOILBSRm0SkLea650TkYRF5SETui7omDr/n4vP887XJ5aSEC0e9ghhuPOOWaNYjPAa8\nZk20qwr9C4ffeEalV7UWflC9qMYsjrBRTaNFZ6drrKPSy6rFyJHuERxeSaNFWi8VkmsxcqQrB8Hj\nb+M6F5Bdi/Hj+w+xFq1F0joS1XjWM3JZtUgzRwLVa5Endau1H0BRRPYHpuFClbwIbOe9l4VzgFtV\ndVfgNuBLMdf1AnNUdT9VPSBNAttt546u9Hnxxf47qxuRtvEM/rP8HkZSVxX6r0qp55GE06vXeOah\nRdh1T6NFmvkRn6hKUk+LYKWs12BUUS6CeWum55lUC5H+2hepxfDhrswF02umjqQhyqhGaTFqlFvC\nGwx4WaQW/m8LzlmUoUWwvahqaKtR/9BfCPqdwON/Ao8sHA380nv+S+CYmOuE5IsC+hBVMMKhPxoR\nbjDq9e7CY/V59TAaDW35xM0JQPVapJ0fgeRDW3564caslbQIG9VW1ELVhTsZCHUkiRYizWuxZYsL\ndxK18CSOYcP6/76y24s4LUaNcga13hG9WajbQKvqqd7THwFHq+ohwO24Sfd/zZj2VH+HuqouB6bG\nZQO4RUTuF5FPpUkgr0oSXIkRFdPHJ4+eZ1QPI244J6rnGVdog1p0d7vCt/326fKWRYtmPJI0Q1tR\nHkkSLTZscBqG46E1oopyUbQWr77qvM64hjaOsCEZDFq89JJLJ808FrRuexHlqeZJUpm+qqq/FZG3\nAu/AeSM/At5U70sicguwbfAtnGGImt6OOGUAgINUdZmIbIMzKAtUNeZE776xtnbbbQ7Ll88B3Ear\n0aPj/6lxhIdz4naOQvRwTh49jKilydC/t5XUbX/pJTceHLVgoB7hBiONFnl4JI28s2a08Icv0gw/\nQvpykXcvvCgt0na0oJo6EowA3Ipa+JHEy9DCj0oMybR48MHqgjb2eH//CfiZqt4gIl+v9wUAVX1X\n3GciskJEtlXVFSKyHRARHARUdZn39xUR+T3uDJNEhqSz08WmUW2+YEQ1nnG9hrJ7W80ObeXZYCTV\nIg+PJG6NPDQ/hFFWuch7XqDVtAj2wsuuI2m0KNqQVN1eJNGiyqCNS0XkJ8BxwJ9EZHSK78ZxHXCi\n9/wEoF+oFREZ550Vj7fU+DAgZhV3f8aMcf+s1aubG9aCdL3wMsd/Id3Q1uTJzrBu2lSNFkV7JGmG\nMPyNXD09+WkRF5oCiikXRSw8yLPxbFUt6k1+B7UYqO1FGi3yJKkx+BBwE3C4qq4BpgBNLfsN8E3g\nXSKyEDgU+G8AEZkmIn4w9W2Bu0TkIeAe4I+qenOaRKZNc4Ujr0pSr8HIo4fR3u4afP8+ebntIq6i\nrFgxcLTYdlvX4PsThHlpMWqUu8/KleU0nnlosf32bkjSX2abVou4xnPiRP/snuZWr0H55cLXwiev\ncrHNNq7j5sehGwh1JIsWeZI0+u9G4HeB18uAZfHfSHTPVbhQK+H3l+GFqFfVZ4GI+KvJ2W47F7M/\nT7e9XsEIFqJmhnNE3E7upUvdBsG8JhKhrxZp95BAei2yLv8dNcpttFqxwmlSz23PooV/sFgaytZi\n4kQ38euXh7y08DsYy5e7XnjcJsx6RBnVpAtEOjrSG69ttnHf85e7N/Lak2oxfLgrby+/7MrFIYek\nyxdEl4u4lV9RWoQ3CTdi++1dOe7tdR52Z2dyw3XRRdFnsDRD1uGplsevJGVMJIavrWcE6jFjRu3U\nuTTLPOuNx0J2LdL0wovQIk1vq2gtqiwXvb31hwqr0CKpUc1Di2HDXAO6dKn7ncOHx8enGuzlYvRo\nVydeftntym9ri184EtZi+fL+MeCaxQxJA9I0nm1tfUMs1Gv46uE3GKrp9pHUyxuUW0ny1GLpUheT\nqasrflgiqsEoUouqysXSpe5eEyb0DxDpM5S0qNfRgqGjxZIl+WuRhiFhSJYsgaefLn4VRltb3wit\nWQvG+vWupxW3TDfsGifpbb3wggubPtC08Htrcb2tqJAsjbR49lmnR15j4WVrEUczWjz9tBsi8QNk\npqEVykUczWjxxBNumGlq3M62OgwmLdJQmSERkQ+IyGMi0lMv3IqIHCEiT4rIUyLyxbTpbLcdXHKJ\nC3aXdvwR0vUwJk/Ot4fR6PvhHka9iT1wWlx0kQs6l2bHrs9g0+KCC9xRrs3ka7Bp8Y1vwIc/nDyE\neZDBpsV//AecfHK6uHA+g0mLNFTpkTyKO8/kjrgLRGQY8EPgcGBP4CN+/K+kbLedGxb50Y/SbzqD\n/pNn9cQP9zAauZpxTJ9ec1XT9jAaVZLubmdMTAv393vfS58nSDfZnlfPs0gtRo6E/2ky6FGalUoD\nQYvx451hbYbBVEfSUJkhUdWFqroIt9s9jgOARar6vKp2AVfiYnQlZs4cuP32/mHCk5JmXqDqHkYj\nV/Xww+G225pbsQUDq7fVSIujj4a//CV9mBifwVQuPvhBd/ZK2jAxPmkm2/OeF8hbi+OPh5tvbi5P\n0Hy58OdD/ZNR01CUFmlo9TmS6bhowz5LvPcSM2oUvOUtzWcg6+RZltU5ebuqY8fGn6WQhKyT7a2k\nxYQJ8MY3ps9P8PuDRYvJk+PPG0lCsx0Mf/VZmY1nIy223todOd0szZaLRqvP6jHoh7ZE5BYReSTw\neNT7+94i082T8PGbjSbPOjqSbRqrx9Spzk29//76E8FpGvY8SKPF5Mk1tz1Lb2v6dLfhav781tai\n3mR7UIstW9wjbWBEqDUYDz8MM2fGX1e1FkknmNetc3mNW31Wj8GgRbBcNNtWQG0F28MPV1dHUsa2\nTEe9WFsJWQoEi8kM771YgrG2wjFlmmHKFNeo+9QTf8QI16NYv941FBs2NOc6Dh/uJsMvuMAZkzgm\nTUpeaPMgjRYTJrjNUd3dbo5q1Ki+JzkmZexYd6+LL+4bqC9MK2sR7Hn6DUYzc1Tt7c4IXXZZ/HGv\n0NpaBD2SZucEwM1lvPoqXHMNLFwYf12VWvT21o9zFlUummH6dHjuOZfuxRfHXzdpklvWPHeuC9r4\n0kvw0582l2aYQg1JCuKq1f3AziIyC7eT/sPAR+rdKGhI8sAv+L29bhVHI3fQ73F1d7sC28zKD3C9\njDe+0Z2hXi+tcK+4yN7WlCmwalWy9ETc71+71hmUZisJOC323js+CjJUr0W9xnPMGOeVdXZmazBE\nnBaHHlp/mW6VWvgdh3rH2fpeexYtRoxwnvtxx7md7nFUqcWGDU6HuDYg6J1l0WL8eFfXPvOZ+l6/\nb7j8DvYFF8C558J3v5s9aGNlhkREjgF+AGwNXC8i81X13SIyDRdh+D2q2iMipwM344bhLlHVBWXm\nc/jwWoPY3t7YHfQNT3d3tsbz1FMbz2eEz5gu2m0fP95tCvTPAE+qxaZN2bQ47TR4Z79gOn0pWws/\nTIkfmqK7u3aWfRS+FlkaDIDPfQ6OOqr+NWVrEW48x4+P97iCXntWLb7wBbdkuR5VatEorQkTXN3o\n7s6uxVe+AiedVP+aoBaqtf9VHlRmSFT1WtwZ8OH3/xFry3v9Z6DJNVf54BeOyZMbFw6/l9HVla1g\nnHBC42vCva2i3XaRmuu+1VauMNYbrvK1yGpIPpXgOLOytRg50g1ZrFvnjMmECfWHq3wtsjYYn/1s\n42vK1mLiRPc/7upK1lD7cwNZtTjrrMbXBLXo7XXLX/NqPKNob3dthR8Is54Ww4bVGvesWnz+842v\nCWqxaZPr+KQ9uCuOVhnaaml8QzJ9uhO+3oFQfs8zqyFJQng1UNFuO9S0GD268Xn0vhYbNxavRbjn\nWZYWq1fXvNZ65OWRJCFYLvxAfvUCWGZFpOahJTFafv7K1sIP2NjscHMSRo92nasNG5JpUVW5yNsz\nM0OSAL+XkUR8/5+1ZUvxBcMfSkk61JQHvhYTJybXYsOGcgzJunWuJ+ifTd3MUso0+Fr4CwLqUWbj\nGTSqGzY4z6nIxhNqWtQ7WtanKi3KqB/QXHtRlha+R5J3R8sMSQL8nmcatz3rBHNSfHfVbyjqjdPn\nQXBoK6kW69cXr4U/9r5hgzMijYaa8sD3ztrayhvOSYI/9t7TU/ywlo9fLrZsaS0tgsM5ZWuRpr1Y\nvbr5DaFJCXskeWphhiQBfoORxIr7/6yyDInf4xo+vJzeVjNarFtXrhb+nEXR+Fok0b7MnuewYS4/\nHR3lDPFBTQvVdFrMmlVsvsaMcQZ18+bytUhTR9ascecPFUlwxVze3tlACNr4nIg8LCIPich9ZebR\nxy8YSV3VsnpbfnodHeW57aZFjbQNxlDQoszJ9iSItLYWZZaLkSPd/M3Gjfkb1ZYO2ujRC8xR1f1U\n9YDis9WftJWkrJ4n9O31l+W2t6oWvkdSphZphjCsXPTNm2lRTbnwjWqeWrR60Ea8zyuNCRacPEuy\nIqXMnqffeFYxkdhqWvgNhmlh5SKIaVGjKC1aPWgjgAK3iMj9IpJgR0H+NNvDaCYwX1pa2W0vWwt/\nDNi0sHIRZKBoMZA7W4VOtovILcC2wbdwhuErqvrHhLc5SFWXicg2OIOyQFXvirs471hbUBvCSDoW\n7i+DnDIlc9IN8XsYkyaVO5yTVItXX3UB5crQwq8kY8e23hBGW5s7V/uFF8otF83Ge0uLXy7GjWts\nHNra3GmMZZWLKoY8n3nGraD0z72Jo60N7r0XVqwoT4u//nUuN944l82bIa+IUq0etNHf6Y6qviIi\nv8edUZLIkORF2h7GfffBwQdnC9melFbvbd15JxxxBOy7b/F58z2Snp7ytPCNar14T+C0uOkmeN/7\nmj8bJw1+udi0qdxy0dbW+OjiyZPhuuvc2R/NHHOclqrqSJK9VpMnw1VXuegNRS//BafFDjvM4cAD\n5zBunAutcv752WNttcrQVuQ8iYiME5EJ3vPxwGFAnbinxeAXjHnzGjcC06bBPvvAr36VX/iBepQ9\nL5BGixkz3DkXl11W/IY4aG0tZs2C178efv7z4ve3QGtrMXu262T9+MemxU47wdveBt//fvH5guKM\naksHbcQNi/1eRNTL6xWqenPZeW1vh5UrXY/yu9+tf+2UKfXDnefNpEmwaFF5bnt7u3PDX3kFfvGL\n+tdOm+YqU1lMmuSGjvzIw0XT3u7Ccj/3XOOgkjvsAA88UHyefPzhnO5ut3m0aNrb4emnXXoHH1z/\n2t13h7//vfg8+QSHtupFTc6L9nYX5l+18ajEfvs5r70sihrma+mgjar6LFDCoEh9xoxx3sXOO7vG\nsZUo222fPNmNu7/97eWM6abB12LECNh228bXZ2XKFDdJ+u53FxsIsBna2mD5crc5s+hNf+C0WLkS\njj22uXNniqStzeWtTI/klVfg4x9v7tCuIimqvWiVoa2WZ8oUeM97Gl9XNmUvbRw+3BmTVtSi7CGM\nceNco9mKWlSx5BVMC6h1sIaSFmZIEvK618H73191Lvrj9zDKGtoCd8jUMceUk1Ya/Mn2srQQcUMT\n723Bg6PLLhcjRrgFFUceWXxaaSlbiwkT3IF0hx1WfFppKUoLi7WVkD//ueocROP3wp97rrxKfEej\nWAQV4WuxcWPjVVR5cc895aSTFl+LtWvL0+Khh8pJJy3BOlKGFiLw+OPFp9MMbW1upeHzz+erhRmS\nAc6kSa6CLFvmltkOZSZNgiefdENOhxxSdW6qZdIk17BPngxveUvVuamWSZPg7rvdno799qs6N9Uy\naRLceqs7tjrPZehVBm38logsEJH5InKNiEyKue4IEXlSRJ4SkS+Wnc9Wp63NTX4ff3zfw4vmzp1b\nWZ6qoq3Nuewnnlib8B2KOkBtCOOTn6xN+A51LT71qdpy46Guxamn5rv0uso5kpuBPVV1X2AR8KXw\nBSIyDPghcDiwJ/AREdmt1Fy2OBMmuD0a4eNoh2JFaWtzleOUU2rvDUUdwHkiw4bBySfX3hvKWowc\n2ff46qGsxdix8NGP5nvfKoM23qqqvd7Le4AZEZcdACxS1edVtQu4Eji6rDwOBIYPh0cegb32qjon\n1TN+PDz6qFumPdTZemtXLmZE1aohxmte44b5yporamV23dXtZ2pry/e+rbJq62Tgxoj3pwMvBl4v\n8d4zAuy5Z9U5aB1MixqmhUPEtPARcSvKcr+vquZ/V//mCYI2ishXgP1V9diI7x8LHK6qp3qvPwYc\noKpnxqRX3I8xDMMYpKhqphmTSoM2isiJwJHAO2IuWQrMDLye4b0Xl14JkXsMwzCMIFWu2joC+AJw\nlKpujrnsfmBnEZklIqOADwPXlZVHwzAMozFVzpH8AJiAO2NknohcDCAi00TkegBV7QFOx63wehy4\nUlUXVJVhwzAMoz+FzpEYhmEYg59WWbWVmCQbFEXkQhFZ5G12rDx6cFE00kJEdhWRv4lIp4icXUUe\nyyKBFh8VkYe9x10iMmgXTCfQ4ihPh4dE5D4ROaiKfJZB0g3NIvJGEekSkRaMqJcPCcrFwSKyxhsh\nmiciX018c1UdMA+c4XsamAWMBOYDu4WueTdwg/f8TcA9Vee7Qi22Bl4PfA04u+o8V6zFgUCb9/yI\nIV4uxgWe7wUsqDrfVWkRuO4vwPXA+6vOd4Xl4mDgumbuP9A8kiQbFI8GLgNQ1XuBNhEp4XSK0mmo\nhaq+qqoPAt1VZLBEkmhxj6qu9V7ew+Ddj5REi42BlxOAXgYnSTc0nwFcDbxcZuZKJqkWTa18HWiG\nJMkGxfA1SyOuGQzYZs0aabU4hegNsIOBRFqIyDEisgD4I25D8GCkoRYisj1wjKr+iCYb0QFC0jry\nZm9K4AYRSbx10aL/GkMKETkEOAl4a9V5qRL1TigVkbcCXwfq7vkaxHwPCM4XDGZj0ogHgZmqulFE\n3o07wfa1Sb440DySJBsUlwKvaXDNYCDVZs1BTiItRGRv4Ke4vUurS8pb2aTdxHsXsKOItNjBybmQ\nRIs3AFeKyLPAB4CLROSokvJXJg21UNX1/rCnqt4IjExaLgaaIUmyQfE64BMAInIgsEZVV5SbzVJI\nu1lzMPe0GmohIjOBa4CPq+riCvJYFkm02CnwfH9glKquKjebpdBQC1Xd0XvMxs2TfFZVB+Om5yTl\nYtvA8wNw20MSlYsBNbSlqj0i4m9QHAZcoqoLROTT7mP9qar+SUSOFJGngQ24YYxBRxItvILxADAR\n6BWRs4A9VHV9dTnPnyRaAP8OTAEuFhEBulT1gOpyXQwJtThWRD4BbAE2AR+qLsfFkVCLPl8pPZMl\nkVCLD4jIZ4AuXLk4Lun9bUOiYRiGkYmBNrRlGIZhtBhmSAzDMIxMVGJIROQSEVkhIo9EfPZ5EekN\nrhYQkS95IU8WiMhh5ebWMAzDqEdVHsmluHPY+yAiM3Dr2Z8PvLc7bjJwd1z4E3+y1DAMw2gBKjEk\n3tr1qHX838WdURLkaFz4+G5VfQ5YhNvubxiGYbQALTNH4m0CelFVHw19NFRCnhiGYQxIWmIfiYiM\nBb5MxjANdma7YRhGejTjMeWt4pHsBOwAPOyFKpgBzBORqaQP+WAPVc4999zK89AKD9PBtDAt6j/y\noEpDIt4DVX1MVbfTWqiCJcB+qvoybhv/cSIySkRmAzsD91WWa8MwDKMPVS3//RXwN+C1IvKCiITD\nmCg1I/ME8FvgCeBPuFg4NoRlGIbRIlQyR6KqH23w+Y6h1/8F/FehmRpkzJkzp+ostASmQw3TooZp\nkS+DKtaWiJizYhiGkQIRQQfJZLthGIYxQDFDYhiGYWTCDIlhGIaRCTMkhmEYRiZaJvqviHzLi+47\nX0SuEZFJgc8s+m8Knn66nHRUYfFgPrTWMIxEtFL035uBPVV1X1xgxi8BiMgeWPTfVOy9N2zcWHw6\njz4K73tf8ekYhtHatEz0X1W9VVV7vZf34EKhAByFRf9NTG8vbNrkHkWzYUM5BsswjNamqQ2JIjIG\neA/wNmB73EHxjwE3qOrjOeTrZODX3vPpwN8Dn1n03zps3tz3b9FplZGOYRitTWpDIiLn44zIXOBe\n4GVgDPBa4L89I/N5Ve13+mHC+38F6FLVXze82OhHZ2ffv0WnVUY6hmG0Ns14JPep6rkxn13gReyd\nGfN5XUTkROBI4B2Bt5cCrwm8rhv997zzzvvH8zlz5gy5UAjmkRiGUY+5c+cyd+7cXO9ZWYgUEdkB\n+KOq7uW9PgL4DvB2VV0ZuG4P4ArgTbghrVuAXaJioViIFHjuOZg9G+bNg/32Kzat3/wGPv5x2LKl\n2HQMwyiOPEKkNDtHMgP4CPBWQnMkwI2BSfO47/8KmANsJSIvAOfiDrYaBdziLcq6R1U/q6pPiIgf\n/bcLi/5bl7I9kq4uN8E/zHYkGcaQpZk5kktxnsH1wDfpO0dyBPAVETlHVe+Mu0dM9N9L61xv0X8T\nUvYcCTiDMnZs8ekZhtGaNOORfEdVH4t4/zHgdyIyiibnSIzsVGFIOjvNkBjGUCa1IYkxIsHPtwAl\n7a02wpQ9tFVWWoZhtC7NDG09ijvBsN9HgKrq3plzZTRNVR6JYRhDl2aGtt6Tey6M3DCPxDCMsmlm\naOv5IjJi5IN5JIZhlE3TizZF5EARuV9E1ovIFhHpEZGOPDNnpMc8EsMwyibL6v8f4vaSLALGAqcA\nFyX5YkwY+XYRuVlEForITSLSFvjMwsgnxDwSwzDKJtM2MlV9Ghiuqj2qeiluH0kSosLInwPcqqq7\nArdhYeSbwgyJYRhlk8WQbPT2jMz3DqX6l6T3iwojDxwN/NJ7/kvgGO+5hZFPgQ1tGYZRNlkMyce9\n758ObMAFVjw2w/2mquoKAFVdDkz13p8OvBi4zsLI16GzE0aMKM8jKSstwzBal6ZibUGf1VudwPn5\nZKdvEs18yaL/wuTJ5XkkZaVlGEY+FBH9t2lDIiIHAecBs4L3UdUdm7zlChHZVlVXiMh2uBhekCGM\n/FCksxPa2srzSMpKyzCMfAh3sM8/P7sfkGVo6xLgAlwE4DcGHkkR7+FzHXCi9/wE4A+B9z8sIqNE\nZDawM3Bf89ke3Gze7Br3sjySstIyDKN1adojAdaq6o3NfDEmjPx/A1eJyMnA87iVWlgY+XSU7ZFM\nmWIeiWEMdbIYkttF5NvA74B/9ElVdV6jL8aEkQd4Z8z1FkY+ITa0ZRhG2WQxJG/y/r4h8J7S95hc\no2T8CfCXXionLRvaMgwjy6qtQ/LMiJEPnZ2w3XbwzDPlpGUeiWEYWVZtnR3x9lrgQVWd33yWjCyU\nPdk+eTKsWVN8WoZhtC5ZVm29Afhn3ObA6cCncSFSfiYi/5ZD3owmsDkSwzDKJsscyQxgf1VdDyAi\n5wI3AG8HHgS+lT17RlrK8khUXRqTJtkciWEMdbJ4JFMJrNbCLc3dVlU3hd5PjIj8i4g8JiKPiMgV\n3t6R2KjARn/K8hK6u0EExo83j8QwhjpZDMkVwL0icq7njdwN/EpExuP2fKRCRLYHzsB5OXvjvKWP\nEBMV2IimLI9k82YYM8Y9zCMxjKFN04ZEVb8GnAqs8R7/rKr/oaobVPX4Jm87HBgvIiNwZ5wsJT4q\nsBFBZ6ebAC/aS+jsdEZk9GjzSAxjqJN6jkREJqlqh4hMAZ7xHv5nU1R1VTMZUdWXROQ7wAvARuBm\nVb3Vj7/lXbNcRKbWvdEQp6yhrc5OZ0TGjDFDYhhDnWYm238FvAc3oR4MVSLe66aCNorIZJz3MQu3\njPgqETme/lGALTxKHYIT4KpuHqOodHyPxIa2DGNok9qQqOp7vL+zc87LO4FnfI9GRH4PvIX4qMCR\nDPUw8p2dbgJcxE2IjxxZXDrmkRjGwKOIMPLSbPxDL4z8fFXdICIfA/YHvqeqLzR5vwNwEYXfiFv1\ndSlwPzATWKWq3xSRLwLtqnpOzD2GfDzHkSNh40Zob4fly2HChGLSmTcPTjkFLr8cPvQhePzxYtIx\nDKNYRARVzTR2kWXV1o9wx+3uA3weWAxc3uzNVPU+4GrgIeBh3FDZT4FvAu8SkYXAobgowUYEPT3u\nMWJE8ZPgvkdik+2GYWTZkNitqioiRwM/VNVLROSTWTKjqufT/7TFVcREBTb64s9biBS/LNeW/xqG\n4ZPFkKwTkS8BHwPeLiLDgIJG5I0k+Etyofi5C1v+axiGT5ahreNwcxmfVNXluJAp384lV0ZT+MNN\nUN7Qlk22G4bRzD4SUcdy3FG7AHiT7JcFr8kvm0YS/OEmKG9oy5b/GobRjEdyu4icISIzg296cbHe\nISK/xJ25bpRMFR7JyJG1SX7DMIYmzcyRHAGcDPxaRGbjwqOMwYU3uRm3BPih/LJoJKUKj0Sk5pWM\nG1dceoZhtC7NGJJuVb0YuFhERgJbA5tU1Y43qpgqPBKozZOYITGMoUkzQ1v3+U9UtUtVl+VpRESk\nTUSuEpEFIvK4iLzJQsknowqPpIy0DMNobZoxJAVFb/oH3wf+pKq7A/sAT2Kh5BNRxfJfsCXAhjHU\naWZoa5uY89oBUNUL4j5rhIhMAt6mqid69+oG1nqbHg/2LvslMBdnXIwAZQ9tjR/vntsSYMMY2jRj\nSGnmCPkAAA2BSURBVIYDEyjGM5kNvCoil+K8kQeAz+FOXrRQ8g0oe2hrq63cc1sCbBhDm2YMyTJV\n/Y/cc+IYgQv+eJqqPiAi38V5HolDyQ/l6L9VT7YbhtH6FBH9txlDUuQcyRLgRVV9wHt9Dc6QJA4l\nHzQkQ42qJtvNIzGMgUO4g33++eHwhulpZrL9v/wnItKeOQcBvOGrF0Xktd5bhwKPA9cBJ3rvnQD8\nIc90BwvmkRiGUQXNeCTnAL/1nv8FNxSVJ2cCV3h7VJ4BTsLNy/xWRE4Gngc+lHOag4KwR7JxY3lp\nmUdiGEOXrENbuQ9zqerDuMOtwlgo+QaEl/+uWlVOWrb81zCGNs0YkrEish9uWGyM9/wfBkVV5+WV\nOSMdNrRlGEYVNLVqi1rU3z4RgHGrqd6RNVNGc2zeDG3enn+bbDcMoyxSGxJVPaSIjBjZMY/EMIwq\naOqERBHZCvgosJv31gLgV6pa4Ki80Qhb/msYRhWkXv4rIrsDjwGvB54CFuEmxx8Tkd3qfdcoFvNI\nDMOogmY8kq8BZ6nqb4NvisixwH8Cx2bJkHf2+wPAElU9ytur8htgFvAc8CFVXZsljcGKeSSGYVRB\nMxsS9wobEQBVvQZ4XfYscRbwROC1Rf5NSFXRf80jMYyhTTOGZEOTnzVERGYARwL/G3j7aFzEX7y/\nx2RJYzBjQ1uGYVRBM0NbU2PCyAuwTcb8fBf4AhA8uMoi/ybEhrYMw6iCZjySnwETIx4T6OtJpEJE\n/glYoarzqb9jPjby71CnLI9E1RkO80gMw4Dm9pFkDxUZzUHAUSJyJDAWmCgilwPLk0b+haEdRr4s\nj6SrC0aMgGFeN8Q8EsMYOBQRRl5U03XwReSrwEWqujrm83cA41T1+qYzJXIw8Hlv1da3gJWq+k0R\n+SLQrqqRpyOKiKb9PYOJ3XeHa66BPfaAZctgv/1g+fL80+nogOnTYd069/raa+HSS+EPFpPZMAYc\nIoKqZoqb2MwcyaPA9SLSCcwDXgHGALsA+wK3At/IkqkQ/41F/k1EWR5JMB0wj8QwhjrNDG39AfiD\niOyCG46aBnQA/wecqqqbsmZKVe8A7vCer8Ii/yairCW5wXSKTsswjNanqRApAKq6CLer3d9EOCEP\nI2I0T3iyffNmNzEuOQf7D6bjp2WGxDCGLs2s2gJARH4lIpNEZDwuZMoTIvKF/LJmpCU45DRsmJsQ\n7+oqNh2wg60MY6jTtCEB9lDVDtwGwRuB2cDHc8mV0RRleQrhdGxoyzCGNlkMyUjvONxjgOtUtQvb\n41EZ3d1uCGtEYLCyKE/BJtsNwwiSxZD8BBdEcTxwp4jMwk26GxUQ9hLAPBLDMMohy2T7hcCFgbee\nFxE79Koiwl4CmEdiGEY5ZJls30pELhSReSLyoIh8n74xspq55wwRuU1EHheRR0XkTO/9dhG5WUQW\nishNIpIpncFIeEkuFOcp2PJfwzCCZBnauhK3GfFY4APe899kzE83cLaq7gm8GTjNOyzLQsk3oMqh\nLVv+axhDmyyGZJqqfk1Vn/UeXwe2zZIZVV3uBW1EVdfjjvCdgYWSb0iVQ1v+BH93d/5pGYbR+mQx\nJDeLyIdFZJj3+BBwU14ZE5EdcCFX7iEUSh6wUPIhqvRIRMwrMYyhTOrJdhFZh1vmK8DngMu9j4YD\n64F/zZopEZkAXI070ne9iISXFccuMx6q0X+r9EiCaU2YkH96hmHkRxHRf5uJtTXRfy4iU3DBGsfE\nfyMdIjICZ0Qu9+J6AaxIGko+aEiGElV6JGAT7oYxUAh3sM8/P/vJIFlWbZ2CC6z4Z+A87+//y5wj\n+DnwhKp+P/DedcCJ3vMTAAtYHqJqj8SWABvG0CXLHMlZwBuB51X1EGA/YG2WzIjIQcDxwDtE5CFv\nafERwDeBd4nIQuBQXGh5I0CVy3+LTMswjNan6Q2JQKeqdooIIjJaVZ8UkV2zZEZV78bNtURhoeTr\nUPbQ1jbblJOWYRitTxZDskREJgPXAreIyGrcwVNGBVQ9tGURgA1j6JIlRMr7vKfnicjtuF3tf84l\nV0Zqqp5sN4/EMIYuWTySf+CdaGhUiHkkhmFURZbJdqOFqNojscl2wxi6mCEZJFTtkdjyX8MYugwY\nQyIiR4jIkyLylIh8ser8tBq2/NcwjKoYEIZERIYBPwQOB/YEPuJFBTY8wsNNc+fOtcl2yD0UxEDG\ntKhhWuTLgDAkwAHAIlV93jvS90pcRGDDIzzcNHfuXJtsxxqMIKZFDdMiXwaKIZkOvBh4vcR7z/Co\nerK9VT0SwzCKJ5flv63Ee99bdQ6q4aGHIBzoeNw4uOuu/DV59lkYO7Z/Wj//Ofz1r/mmlZWFC+HB\nB6vORWtgWtQwLeB734OddsrnXqIaG5G9ZRCRA4HzVPUI7/U5gKrqN0PXtf6PMQzDaDFUVbJ8f6AY\nkuGAH7BxGXAf8BFVXVBpxgzDMIyBMbSlqj0icjpwM25e5xIzIoZhGK3BgPBIDMMwjNZloKzaqstQ\n36woIs+JyMPeGS73ee+1i8jNIrJQRG4Skbaq81kEInKJiKwQkUcC78X+dhH5kogsEpEFInJYNbku\nhhgtzhWRJd7ZPv75Pv5ng1mLGSJym4g8LiKPisiZ3vtDrmxEaHGG935+ZUNVB/QDZwyfBmYBI4H5\nwG5V56tkDZ4B2kPvfRP4N+/5F4H/rjqfBf32twL7Ao80+u3AHsBDuCHdHbxyI1X/hoK1OBc4O+La\n3Qe5FtsB+3rPJ+DmWHcbimWjjha5lY3B4JHYZkUQ+nuXRwO/9J7/Ejim1ByVhKreBawOvR33248C\nrlTVblV9DliEKz+DghgtwJWPMEczuLVYrqrzvefrgQXADIZg2YjRwt+Hl0vZGAyGxDYrguIOF7tf\nRE7x3ttWVVeAK0jA1MpyVz5TY357uKwsZWiUldNFZL6I/G9gKGfIaCEiO+A8tXuIrxdDQo+AFvd6\nb+VSNgaDITHgIFXdHzgSOE1E3oYzLkGG8qqKofzbLwZ2VNV9geXAdyrOT6mIyATgauAsrzc+ZOtF\nhBa5lY3BYEiWAjMDr2d47w0ZVHWZ9/cV3NHHBwArRGRbABHZDni5uhyWTtxvXwq8JnDdoC8rqvqK\negPfwM+oDVEMei1EZASu4bxcVf/gvT0ky0aUFnmWjcFgSO4HdhaRWSIyCvgwcF3FeSoNERnn9TQQ\nkfHAYcCjOA1O9C47AfhD5A0GB0Lfsd64334d8GERGSUis4GdcZtbBxN9tPAaS5/3A495z4eCFj8H\nnlDV7wfeG6plo58WuZaNqlcU5LQq4QjcSoRFwDlV56fk3z4bt1LtIZwBOcd7fwpwq6fLzcDkqvNa\n0O//FfASsBl4ATgJaI/77cCXcKtQFgCHVZ3/ErS4DHjEKyPX4uYIhoIWBwE9gboxz2snYuvFYNWj\njha5lQ3bkGgYhmFkYjAMbRmGYRgVYobEMAzDyIQZEsMwDCMTZkgMwzCMTJghMQzDMDJhhsQwDMPI\nhBkSY1AiImeKyBMicnnVecmLQNjv87zXJ4jID0LX3C4i+9e5x/+JyEoReX/B2TWGEAPihETDaILP\nAIeq6kvBN0VkuKr2VJSnPLhAVS8IvE61EUxVPyYiP885T8YQxzwSY9AhIj8CdgRuFJGzvJ78ZSJy\nF3CZiAwTkW+JyL1e5NNPBb77Q+8wn5tF5Aa/5y4iz4rIFO/560Xkdu/5OO9AqXtE5EERea/3/gki\nco2I3OgdovTNQBpHeNfOF5FbxPGUiGzlfS7eoUJbZdDgveIOOpsn7tC3xcGPm72vYURhHokx6FDV\nz4jI4cAcVV0tIufiDus5SFW3eIZjjaq+yYvPdreI3AzsD+yiqruLyDTgCeAS/7bhZLy/XwH+oqqf\n9MJw3ycit3qf7YML2d0FLBSRC3HhS34KvFVVXxCRyaqq3hDcx4DvA+8E5qvqygQ/98Mi8lbvuQA7\neRr8EfgjgIj8Brg9iXaG0QxmSIzBSr9Ajqq6xXt+GLCXiHzQez0J2AV4O/BrcBGVReS20P2iOAx4\nr4h8wXs9ilo06r+oC9eNiDyOO8VzCnCHqr7gpbPGu/ZSXLyj7wMne6+TcKWqnvmPTPbNMyLyb8BG\nVf1xwvsZRmrMkBhDhQ2B5wKcoaq3BC8QkX+q8/1uakPBY0L3OlZVF4XudSDO+/DppVbf+hklVV0i\n7rz1Q4A3Ah+tk5d6BCP/vhM4Fnhbk/cyjETYHIkxFLkJ+Kx3RgMisouIjAPuBI7z5lCmAYcEvvMs\n8Hrv+bGhewU9gn0bpH0P8DYRmeVd3x747BLg/4DfasZoqt79fwh8MOCJGUYhmEdiDFbqNcT/C+wA\nzBMRwR1udIyq/l5E3gE8jgvD/rfAd/4DuERE1gJzA+9/DfieiDyC65g9gzv/OzI/qvqqiJwK/D6Q\n9uHeNdfhzo34RfKfGZ0O7qyNKcC1XjpLVfU9Ge5rGLFYGHnDiEFELgX+qKq/Kym9NwDfUdWDYz4/\nF1ivqpmOyy37dxmDHxvaMox4SutlicgXgauAc+pcth74lL8hscl0/g+3qKCz2XsYRhjzSAzDMIxM\nmEdiGIZhZMIMiWEYhpEJMySGYRhGJsyQGIZhGJkwQ2IYhmFkwgyJYRiGkYn/DyaMeMlSwxcLAAAA\nAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import scipy\n", "import matplotlib.pyplot as plt\n", "pi = scipy.pi\n", "\n", "signal_length = 0.5 #[seconds]\n", "sample_rate=500 #sampling rate [Hz]\n", "dt = 1./sample_rate #time between two samples [s]\n", "\n", "df = 1/signal_length #frequency between points in\n", " #in frequency domain [Hz] \n", "t=scipy.arange(0,signal_length,dt) #the time vector\n", "n_t=len(t) #length of time vector\n", "\n", "#create signal\n", "y=scipy.sin(2*pi*50*t)+scipy.sin(2*pi*70*t+pi/4)\n", "\n", "#compute fourier transform\n", "f=scipy.fft(y)\n", "\n", "#work out meaningful frequencies in fourier transform\n", "freqs=df*scipy.arange(0,(n_t-1)/2.,dtype='d') #d=double precision float\n", "n_freq=len(freqs)\n", "\n", "#plot input data y against time\n", "plt.subplot(2,1,1)\n", "plt.plot(t,y,label='input data')\n", "plt.xlabel('time [s]')\n", "plt.ylabel('signal')\n", "\n", "#plot frequency spectrum \n", "plt.subplot(2,1,2)\n", "plt.plot(freqs,abs(f[0:n_freq]),\n", " label='abs(fourier transform)')\n", "plt.xlabel('frequency [Hz]')\n", "plt.ylabel('abs(DFT(signal))')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The lower plot shows the discrete Fourier transform computed from the data shown in the upper plot." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Optimisation\n", "------------\n", "\n", "Often we need to find the maximum or minimum of a particular function *f*(*x*) where *f* is a scalar function but *x* could be a vector. Typical applications are the minimisation of entities such as cost, risk and error, or the maximisation of productivity, efficiency and profit. Optimisation routines typically provide a method to minimise a given function: if we need to maximise *f*(*x*) we create a new function *g*(*x*) that reverses the sign of *f*, i.e. *g*(*x*)= − *f*(*x*) and we minimise *g*(*x*).\n", "\n", "Below, we provide an example showing (i) the definition of the test function and (ii) the call of the `scipy.optimize.fmin` function which takes as argument a function *f* to minimise and an initial value *x*0 from which to start the search for the minimum, and which returns the value of *x* for which *f*(*x*) is (locally) minimised. Typically, the search for the minimum is a local search, i.e. the algorithm follows the local gradient. We repeat the search for the minimum for two values (*x*0 = 1.0 and *x*0 = 2.0, respectively) to demonstrate that depending on the starting value we may find different minimar of the function *f*.\n", "\n", "The majority of the commands (after the two calls to `fmin`) in the file `fmin1.py` creates the plot of the function, the start points for the searches and the minima obtained:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: -2.023866\n", " Iterations: 16\n", " Function evaluations: 32\n", "Start search at x=1., minimum is [ 0.23964844]\n", "Optimization terminated successfully.\n", " Current function value: -1.000529\n", " Iterations: 16\n", " Function evaluations: 32\n", "Start search at x=2., minimum is [ 3.13847656]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEPCAYAAACneLThAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8TNf7B/DPCaKJJGhCEJGkoXRRqX4ttSZVe9Ha1VKq\nqguKFi3VSH2/LX76pYqW1k7Vt6pqi9aS2IuUUCWCLHZBqFqHzPP74ySZJCYzk8yde+/MPO/XK69m\nZu6c++S4febOc889RxARGGOMuSYPrQNgjDHmOJzkGWPMhXGSZ4wxF8ZJnjHGXBgnecYYc2Gc5Blj\nzIUpkuSFEG2EEElCiGQhxBgzrzcXQlwXQhzI/vlYif0yxhizrKS9DQghPADMBNACwHkA+4UQvxBR\nUoFNtxNRR3v3xxhjzHZKnMnXB3CCiNKJ6D6AHwB0MrOdUGBfjDHGikCJJB8E4Eyex2eznyvoeSFE\nohBivRDiSQX2yxhjzAq7yzU2+gNANSK6LYRoC2A1gMdV2jdjjLktJZL8OQDV8jyumv1cLiK6mef3\nWCHEbCHEo0SUWbAxIQRPpsMYY0VERGZL4kqUa/YDqC6ECBFCeALoCWBN3g2EEIF5fq8PQJhL8HmC\n1fQnOjpa8xj08sN9wX3BfaH/vrDE7jN5IsoSQgwB8Bvkh8Y8IjomhBgsX6a5ALoKId4GcB/AHQA9\n7N2vI6WlpWkdgm5wX5hwX5hwX5jovS8UqckT0UYANQs8NyfP77MAzFJiX4wxxmzHd7ya0b9/f61D\n0A3uCxPuCxPuCxO994WwVs9RmxCC9BYTY4zpmRAC5MALry4nPj5e6xB0g/vChPvChPvCRO99wUme\nMcZcGJdrGGPMyXG5hjHG3BQneTP0XmNTE/eFCfeFCfeFid77gpM8Y4y5MK7JM8aYk+OaPGOMuSlO\n8mbovcamJu4LE+4LE+4LE733BSd5xhhzYVyTZ4wxJ8c1ecYYc1FXrlh+nZO8GXqvsamJ+8KE+8KE\n+8JE675YsMDy67pM8gaD1hEwxpj+GY3AN99Y3kaXNfkffiD00PXaUYwxpr2NG4GxY4GDB52sJv/V\nV1pHwBhj+jd7NvDuu5a30WWST08HDh7Ubv9a19j0hPvChPvChPvCRKu+SE8Hdu0Ceva0vJ0uk/zb\nbwOzeEVYxhgr1Ny5QN++QJkylrfTZU0+I4Pw+OPAqVPAo49qHRFjjOnLvXtASAiwbRtQs6YTjpOv\nUAHo2BGYP1/rSBhjTH9WrQKeekomeGt0meQBWbL55hs5REhtXG804b4w4b4w4b4w0aIvvv4aeOcd\n27bVbZJv0ADw9QU2b9Y6EsYY048//5Sl7I4dbdtelzX5nJjmzgViY4Gff9Y4KMYY04l33gEqVgQm\nTDA9Z6kmr+skf/MmUK0acPgwULWqxoExxpjG/vlHXnD9808gKMj0vNNdeM3h4wO8+irw7bfq7pfr\njSbcFybcFybcFyZq9sXSpcALL+RP8NboOskD8gLsd98B9+9rHQljjGmHSN7h+vbbRXufrss1OZo1\nA957D+jSRaOgGGNMYzt3Am+8ARw7BogChRmnLdfkePttOWSIMcbcVc5ZfMEEb41TJPnOneWFhuRk\ndfbH9UYT7gsT7gsT7gsTNfri0iU50vC114r+XqdI8qVLA6+/bn3eZMYYc0Xz58tydblyRX+vU9Tk\nASAtDfjXv4AzZwAvL/XjYowxLWRlAY89Ju8XqlvX/DYOr8kLIdoIIZKEEMlCiDGFbDNDCHFCCJEo\nhIgo6j5CQ+VdsCtW2B0uY4w5jdhYoFKlwhO8NXYneSGEB4CZAFoDeApALyFErQLbtAUQTkQ1AAwG\nUKzCi1oXYLneaMJ9YcJ9YcJ9YeLovpg92/Z5asxR4ky+PoATRJRORPcB/ACgU4FtOgFYDABEtBdA\nWSFEYFF31LatvABx4IC9ITPGmP6lpAD79wPduxe/DSWSfBCAM3ken81+ztI258xsY1WJEsCbbzr+\nbD4yMtKxO3Ai3Bcm3Bcm3BcmjuyLOXPkiBp7rkOWVC4c5fTv3x+hoaEAgHLlyiEiIiK3I2vVisdn\nnwFTp0aibFnTV6Wc1/kxP+bH/NgVHt+9C8yZE5+9Sl7+13N+T0tLg1VEZNcPgIYANuZ5/CGAMQW2\n+QZAjzyPkwAEFtIeWdOjB9GMGVY3K7a4uDjHNe5kuC9MuC9MuC9MHNUXixcTtWpl27bZedNsjlai\nXLMfQHUhRIgQwhNATwBrCmyzBkA/ABBCNARwnYguFXeHORdgdTb6kzHGFFOUhUEsUWScvBCiDYAv\nIWv884hokhBiMOSny9zsbWYCaAPgFoABRGT28mlh4+TzIgKeflpedW7e3O7wmcoePAC2bwd27wYa\nNgSaNpU3vDHGpIMHgU6d5IXXkjYU1Z12PnlLvvoK2LUL+OEHFYJiijh+HJg6FVi9Wt730KQJ8Pvv\nwF9/yQ/rzz4DatfWOkrGtPfmm3ItjY8/tm17p5+gzJx+/YBff5VDKpWW9+KGu1OqL377TZ6xh4TI\nIWH79wPTpgF79gCpqUCbNkDLlvIMRq/4uDDhvjBRui/+/hv48Uc546QSnDbJly0LdO0KzJundSTM\nmlmz5IfyTz/JM5PsgVO5/P2Bd9+V5bc2bYCEBE3CZEwXFi8GWreWd7kqwWnLNYC8KeqVV2TdqkQJ\nBwfGimXUKGD9emDdOjn/hjW//CK/qq5dC9Sv7/j4GNMTIuCpp+RF16Jcb7RUrtHlOHlb1a0rP+1i\nY4GXXtI6GlbQ3Lkyue/ZY/vseZ06AUYj0LMncPQo8Mgjjo2RMT3Ztk3OF9+smXJtOm25Jocj5rPh\neqNJcftizx5Zmlm9uujTo77yClCnjqzZ6wkfFybcFyZK9kXOPDVFXRjEEqdP8j16AHv3yot3TB8u\nXAC6dZNzYNesWbw2pk4FvvhCtsWYO7hwAdi0CejbV9l2nbomn2PkSDnO+vPPHRQUs9n9+0BkpLxw\n9Mkn9rU1ZgyQkQEsWKBIaIzp2sSJwLlzxVscySXHyeeVnCyH550+zTfVaO3zz4H4eHmdxMPO74k3\nbgC1agFr1sgFYxhzVQ8eyFFn69fLUmVRueQ4+bwef1zeRLNqlTLtcb3RpCh9kZwsSyxz5tif4AHA\nz0+e3Qwfro8pLPi4MOG+MFGiL9aulfeQFCfBW+MSSR5Qb0ERZp7RKIc+mhsHb4/+/eUNb3v3Ktcm\nY3qj1Dw15rhEuQaQteDQUHkX7NNPKx8Xs+y77+SQyT17lL9nYfJk+S2Bb3xjrkiJcrPL1+RzREcD\nV68CM2cqHBSz6MIF4JlngC1b5H+VdumSHKVz+rQs4TDmSkaOBDw9gUmTit+Gy9fkcwwaBHz/PXDz\npn3tcL3RxJa+GDMGGDjQMQkeAAIDgRdflP+2WuLjwoT7wsSevrh1S05jMHiwcvEU5FJJvmpVeaeY\n1snAnSQkAJs3A+PGOXY/gwYB337r2H0wprZly4DGjYGwMMftw6XKNYCsyX/4oZzXRsm7xtjDiOT8\nGn37yiTsSEajnPtm1So5nQVjzo5IfvudNk1+U7WH25RrADld7T//8GgMNfz8M3D9OvD6647fl4eH\nLAnx2TxzFdu2AVlZQIsWjt2PyyV5Dw9Z37JnOCXXG00K6wuDARg9Wo6LV2sG0AEDgBUrZB1TC3xc\nmHBfmBS3L776ChgyxPEVB5dL8oBMBr/8IkfaMMeYNUvehNaypXr7rFoVaNRIuZveGNPK6dPyzvB+\n/Ry/L5eryefo2xeIiADef1+BoFg+N24A1asDcXFy7ms1LVgAbNggV85hzFl9+CFw9y4wfboy7bnN\nOPm8du8GXntNriuqxC32zOTTT4ETJ4AlS9Tfd0aG/AZx6RLPU8Sc0507cv3W3buBGjWUadOtLrzm\neP55oEwZubZoUXG90aRgX2RmAjNmyBvPtFCxIvDEE/Kildr4uDDhvjApal8sXy5XPVMqwVvjskle\nCGDYMJmQmHK++AJ4+WVZrtFKx45yQifGnA2RvOA6dKh6+3TZcg0ga14hIcD27cVfvIKZXL4sp/49\ncED2q1aOHJHLPaam8r0QzLns3CmHHCclKVtGdstyDSDXBx00SH5yMvtNnizXXtUywQPyYq8QwJ9/\nahsHY0U1Y4YcNqnmdUKXTvKAnL7z++/lTTu24nqjSU5fnD8vl/Nz9PQFthBCm5INHxcm3BcmtvbF\n2bNyCpD+/R0azkNcPslXqQK0aSMTFCu+zz+XB2eVKlpHInXowHV55ly++QZ49VX1Z1J16Zp8jr17\nZZnh5En17s50JadPy3sOkpLk6BY9MBhkLElJQKVKWkfDmGU51we3bZPXtZTmtjX5HA0ayOlq+cyv\neP79bzlVhF4SPCDn327dWq6JyZje/e9/8kTJEQneGrdI8oBcJ3TaNNu25XqjybJl8Vi1Chg1SutI\nHtahA7BunXr74+PChPvCxFpfEMncM2yYOvEU5DZJvmtXID0d2LdP60icy+LFcjTAo49qHcnDXnhB\nDo81GrWOhLHCbd0qy4tt22qzf7eoyeeYPl2uQbpihUOadzlJSXLtyZMngbJltY7GvJo15VdhR6xy\nz5gS2rYFunVz7JTcbjl3jTn//CNXYNm/37ErsbiKnj1l8vzoI60jKdybb8px8++9p3UkzikxUa7u\ndfQocOyYnFOlVy/54c6DFOx35IicqTUtzbFzLTnswqsQorwQ4jchxHEhxK9CCLPne0KINCHEISHE\nQSGEZgUTX1/gjTesz/zG9UZ5o1F8PBAREa91KBZFRqo3j40rHRf378sFpNu3B3bskCOU3n5brr41\nfLhM9jExcjtzXKkv7GWpL774Qk5hoOVkevbW5D8EsJmIagLYCqCwcz4jgEgiepaI6tu5T7sMHSpn\nT8zM1DIK/YuOlhdbvby0jsSy5s1lkue6vO3OnpUfjsnJ8sN80SK5AEzHjnJR9sREedPO3r1AVBRw\n7pzWETun8+fluhZvvaVtHHaVa4QQSQCaE9ElIUQlAPFE9NAgISFEKoB/EZHVZTwcWa7J0b+/rOXq\nuQyhpQMH5MiVEycAb2+to7GuRg3gp5/kepnMssREWSMeNkwmdEu31xuNwKRJclqQJUvsX4fU3Xz0\nkVzFTI1JEh1WkxdCZBLRo4U9zvN8CoDrALIAzCWiQlfqVCPJ//mnHGOdkiLnt2H5vfSSvEt4yBCt\nI7HNoEFA7draDVFzFrdvA889B4wdKxfVsVVcnLxT88svge7dHRefK7lxQ5a+9u2T/3U0S0m+pA1v\n3gQgMO9TAAjAx2Y2Lyw7NyaiC0KICgA2CSGOEdHOwvbZv39/hIaGAgDKlSuHiIgIREZGAjDVv+x9\n/NxzkViwAHjiiYdfT0xMxPDhwxXdn7M8nj07Hvv2AT/9JB9Pnz7dIf2v5OOKFYFt2yIxbJhj95e3\n9qqnv9/Wx2PGAJUrxyM4GABsf78QwK+/RqJlS+DAgeX4I3k+Lt++jKcffxrtG7VH5UqVdfH3afXY\nXL7YuzcSrVsDp0/H4/Rp5fef83taWhqsIqJi/wA4BiAw+/dKAI7Z8J5oACMtvE5q2L2bKCSEyGB4\n+LW4uDhVYtCjF18kmjPH9NgZ+uLMGSJ/f6KsLMfuxxn6ojAbNxIFBxNlZha/jR9XppDHk+GEsSC8\nBsJYUHj7cEpJTVEuUCdU8Li4fZuoUiWiP/9UL4bsvGk2p9p74XUNgP7Zv78G4JeCGwghvIUQPtm/\nlwHQCsARO/drt+efl1+jli9/+LWcT013s307cOqUXAg9hzP0RdWqQLlywF9/OXY/ztAX5ly9Cgwc\nKNfHLV+++O2s3jYexpdPAZ4AwgB4AqfqnML4/45XKlSnVPC4mDdPTqXy9NPaxFOQvUl+MoCWQojj\nAFoAmAQAQojKQoicG84DAewUQhwE8DuAtURUjEX5lDd2rJxdkUdmyFuvx48HPvkEKFVK62iKLjJS\nDvlkD3v/fXkzTosW9rVz7sY5meDz8gTO3zhvX8MuxGAA/u//ZG7RC7uSPBFlEtGLRFSTiFoR0fXs\n5y8Q0UvZv6cSUQTJ4ZO1iWiSEoEroUULOe3nqlX5n493w2yxZQtw8SLQp0/+552lL9QYL+8sfZFX\naqqc3ycmxv62gvyCAENOw9n/NQBV/HQy/7RG8h4Xy5bJhebrazpQPD+3mbvGHCHkJ+5nn8kzWXeV\ncxY/YQJQ0uqleH3i8fLmffGFvCtYiTnMJ46ciPBD4aZEbwAejQ/HxJET7W/cBWRlySGnejqLB9xs\nWgNzjEZ56/7nn8uhg+5owwZ5M8zhw+ouS6a00FDg1195Pd8cGRlyattjx+RU20pITUvF+P+Ox/kb\n51G2ZBXsXDMRP/4vDE56uUJRy5fLewp27VJ/7WGeu8aKVavknOl//OF+C0MTAf/6lzz76NJF62js\n07OnvNHntde0jkQfxo0Drl0DZs923D42bZI3Fx46BAQEOG4/evfggZxDadYsbW4ac/tFQ6x55RX5\n359/lv91xtprcf30k0z0OX1QkDP1RYMG8lZ8R3GmvrhxA5gzB/jgA8e0n9MXLVvKD9e33nLfkmd8\nfDyWLpXz/9h7cdsROMlDnr1PnChHlmRlaR2Neu7fl7deT5ni3GWaHA0bOjbJO5M5c4BWrdS52/I/\n/5HTUi9d6vh96dH9+8Cnn8ocosdKAJdrshEBjRrJW+N79VJ995qYOVOOvNi4UetIlHH3LuDvD1y+\n7Bxz7jjK/fvy+sSGDerNs5+YKD9UEhLkDJbuZM4cWfL99VftYuByjQ1yzuajo2V9zdXduCH/3smT\ntY5EOY88IuuiBw5oHYm2fvtNJnk1F1KJiJBTF7/2mnuNcLp7V17Pm6jjAUac5PNo0QKoUgUYNy5e\n61AcbsoUOQmZtUTgTHVoQNblf//dMW07S18sWyYnFHMkc30xapT8FmFtvQZXMmcOEBwcr6tx8QVx\nks9DCDlmfv58OWOfqzp3To640PPZR3G5e13+5k1ZptFitsgSJeSawJ9/LldEcnXXr8t8MXCg1pFY\nxjV5M7p3l1PXjnfRKTkGDAAqVnStUk2OkyflQhdnzmgdiTaWLQO+/x5Yv167GObNk+PF9+7VdkUk\nRxs9Wi4+9N13WkfC4+SLLCUFqFdPno1UrqxpKIrbu1cOl0xKUuYuSL0hAipUkOO2g4K0jkZ97drJ\nqSkcXa6xhAjo1EleH/n8c+3icKTUVHl/iV5yBF94LaLTp+Px+utySKUrMRrl8oeTJtme4J2lDp1D\nCMeNl9d7X2RkALt3y2X8HM1SXwgBfPstsHAhsLPQVSOc24cfyrVwK1fW/3HBSb4Q48YBa9bIW/1d\nxcKFcm6agpOQuRp3rcv/739yYW4fH60jkdMofPMN0K8f8M8/WkejrD175Ifp++9rHYltuFxjwcyZ\nciHe337T500ORXH9OvDEE3Jc/HPPaR2NY23aJIe1OXpWSr1p1Aj4+GNZstGLgQPl/zt6qFsrwWgE\nGjcG3n5bfoDpBZdrimnwYDn97o8/ah2J/SZMkItzu3qCB+T1lAMH3ON+hxwpKfKic8uWWkeS3/Tp\nwNat8luxK5g/X15zcKZvw5zkzcipsZUqJcfBjhghz4SdVUKCnCHvP/8p+nv1Xm80p1w5IDhY+ZWi\n9NwXq1YBnTurt+CLrX3h6yuHVQ4eLK8ZOLOMDDmR35w5+acB0fNxAXCSt6pRIzkF8bhxWkdSPAaD\n/Mr8xRdy1Im7cORNUXq0fr1+p8pu0kTOVPnGG849idn778s7etW8k1gJXJO3wbVrcjjYzz/L5OFM\n/v1veZFo/Xrnv65QFDNnAn/+Kc+6XN3ff8t1bi9eBMqU0Toa8wwGecLUvz8wZIjW0RTdli3yZOmv\nv/TZx1yTt1P58sDUqXKFnfv3tY7GdkePAl9+KROdOyV4QF57+OMPraNQx+bN8mKgHpNPDk9P4Icf\n5DKEiYlaR1M0d+/KqZS/+krffVwYTvJmmKux9eolx8Q6y80dWVnyzOPTT2V9urj0Xm8sTJ06ckUk\ng8H6trbSa19s2KD+iJri9EX16vKko0cPOf2Cs/j4Y3k8dehg/nW9Hhc5OMnbSAh5ZX3WLOeo9U6e\nLM+eBg/WOhJteHvLudRdfQ4Vo1GbJF9cr74qv3U4S8lm82Y5aOGbb7SOpPi4Jl9Eq1bJOSsOHpQj\nB/Ro1y450uKPP2St1l317y8TyqBBWkfiOAcOyG+Zx49rHYntbt2SUwK8/768GKtXV6/KKZTnzZNz\n5esZ1+QV1LkzEBkJvPee1pGYl5kpz5a++869EzzgHnV5ZzqLz1GmDLB6tRyOuHu31tGYRyS/BXfr\npv8Ebw0neTOs1dimTwd27JC3kesJkZxhsnPnwuuHRaX3eqMlSid5PfaFVkne3r6oWRNYsEAm0XPn\nlIlJSfPnA8nJcipha/R4XORVUusAbBUaGor09HStw8inRw/5o0futHCDNTkji0JCQpCWlqZpLEq6\nckUO6WvWTOtIiqd9e+Ddd+WsqNu3y5W99GDvXrn2cVycfmKyh9PU5LNrThpExFyFqx1Dy5bJKTdW\nr9Y6kuIjAnr2lB/Ey5bJhUe0dOaMnOBuzhz93lxmDtfkGXNBsbHOV48vSAg5O+rFi3LEjZafwbdu\nyWmahw93rgRvDSd5xmykp9orkbwLU6sJyZTsCy8vOYHZ/v1yTLoWsrLkrJJ16gAffFC09+rpuDDH\naWryjDGT48flfRChoVpHogw/P/nNpFkzOcHcqFHq7fvBAzknzd9/y6UTXe3ucK7JM7fhSsfQN9/I\nm/IWLtQ6EmWdPSvX6O3WTc6a6uiEe/++nDb4+nV5bcPLy7H7cxSuyTPmYuLi5P0arqZqVTl2futW\neXat5LQUBRkM8qLvrVtycSBnTfDWcJJnLmHfvn1Yv349tmzZ4rB96KX2SgTEx2ub5B3ZFxUqyCT/\n99/ywnJmpvL7OHNGfmMgAn76yb6hkno5LgpjV5IXQnQVQhwRQmQJIepa2K6NECJJCJEshBhjzz5d\nXWpqqtYhFEtWVhY2bdqk6j4PHTqEixcvAgCSk5PRvn177NixQ9UYtHDsmJybx1Xq8eZ4e8spROrU\nAWrXllNlK2X9erl6WMeOwMqVQOnSyrWtS0RU7B8ANQHUALAVQN1CtvEAcBJACIBSABIB1LLQJplT\n2POuJCUlhZYvX25xm/T0dPrhhx9Uish2y5Ytozt37tjVRmZmJk2aNIkWLFhACQkJD71+4sQJWrVq\nFcXExNAff/xBREQLFizIfT0pKYnWrFlTaPuucgzNmkU0YIDWUagnLo4oLEz+zdeuFb+djAyioUOJ\nqlYl2rFDsfB0IfvYNptT7TqTJ6LjRHQCgKXLI/UBnCCidCK6D+AHAJ3s2a+r+uabb9CzZ0+L21Sr\nVg23b9/G0aNH7d7f2rVrsWzZMnz66aeYPXu2XW1dunQJj9h5e+DChQsRFRWFPn364L///e9Dr69d\nuxZBQUEYMWIEpk6dCgDw9vbGlStXQEQ4evQoWuptkVMHcNV6fGEiI4HDh2VJJTxcjrw5c8b0+ohP\nRqD5a80R2T8y96f5a80x4pMRAOSiP+PGAbVqyaGSBw/K1archRpDKIMA5PknwVnIxM/yOHz4MIJt\nnPj91VdfxYgRI+xKzH///Te6d++O69evw9PTEwEBAWjfvj1CQkKK3Nbt27dRKntx0QsXLuDo0aPY\nsmULKlWqhCeffBIvvviiTe2kpKSga9euKFmyJK5du/bQ6yNGyP9pjx07hrCwMABAzZo1cfjwYSQl\nJeHkyZO4ffs2evfuXeS/wRbx8fGI1Di7Go2yHv/FF5qGoXpf+PgAs2fLGWC//FKWcV58UX4AVC7f\nGAln5+J2yO3c7b1SvVEzcxg6d5b91bmznLGzGIe3VXo4LiyxmuSFEJsABOZ9CgABGEdEax0VmLtZ\nu3YtXn75ZZu2LV26NAwGA27evAkfH59i7a9s2bJISEhA6eyCZFZWltXhhVeuXMG2bdsg8oxr8/f3\nR40aNVAmz5I5LVq0wI8//ojo6Ojc5G/p/c2bNwcAGI1GlLDhvvbVq1djXPaiuz4+PkhKSsI777xj\nw1/t/I4elWPKq1XTOhJthIYC06YB0dGynr5nD7B3bxfcuTEVeH1vbnZ6sL027jbpjM6dga+/BgID\nrbXsuqwmeSKy9/vvOQB5D8mq2c8Vqn///gjNvqpUrlw5RERE2BmC43399de4desWvLy8ULp0abzx\nxhtYvnw5MjMz4enpCQ8PDwwcOBCHDx/Gvn374O/vj8WLF+Pnn38GAOzfvx9jx461eX916tTB7t27\n0cqOeVCfeuopAMCOHTvQrFmz3D4vTEBAALp06fLQ83fu3MHN7KV+ypQpg4yMDFSsWBEGgwH37t2D\nn5+fxffnqFWrFjIyMuDv74+yZcua3Wbt2rUYMmQIzp07hxo1auDvv/+Gv7+/LX8uANNIiJwzr6I8\njoyMtOv99jwOCwnBwvHjsTf+CPypAtJT5yIkLEyzePTw+I03gOrV49G3L3A24wMMjn0Nt+k2Sl8s\njWVfjkKXDgLx8fE4dgwIDHRsPDnU+vtzfrdpwr3CivVF+QEQB+C5Ql4rAdOFV0/IC69PWGjL0oUF\nXdqxYwe1a9eOiIgOHDhAgwYNooSEBBo0aFDuNsOHD6dt27bRkCFD6MyZM0REtGjRotzXW7Vqla/N\nX375hdatW0djxoyhpUuXUp8+fSgpKSn39ZUrV9LMmTPtjv3777+nbt260YkTJ3Kf++uvv2j06NG0\nfv16iomJsamdL7/8koiIRo0aRStWrKBJkybR6tWrixTLlStXaMqUKTR37lzavXs3/fXXX/TZZ5/l\nvr5q1SqqV68etWzZkv79738TEdHixYvp+vXrNrWv52PIkrSUFHo/PJxuyhF/dBOg98PDKS0lRevQ\ndMNoNFKDrg0I0aAGXRuQ0WjUOiRVwcKFV3uT+8uQ9fY7AC4AiM1+vjKAdXm2awPgOIATAD600qal\nP8LCH6nMT3GMGDEiXzIiIhozZgxNmzYt9/HcuXPp7bffpq1bt1LFihWpa9eutHPnztzXW7Rokft7\nenp6btJQKbusAAAe6ElEQVStW7cuXbt2jdatW0e3b9/O3Wbz5s0P7ZOIaPLkyRQTE5PvZ8KECRQT\nE0NpaWlm479x4waFh4dTWloaZWRkUEhICGVkZBAR0dixY23qg2XLlpHBYLBpWyXlHV1jDQCy5//9\nuLi44r/ZDhN6985N8JQn0U/o3VuTeIi06wtLfvzlR/Jt5ksr16xUdb966AtLSd6uC69EtBrAQxOd\nEtEFAC/lebwRcrilw1gpJzuU0Wh8qJ599+5dGPLcrnf//n3cv38fjz/+OI4ePYoNGzbgzTffxNat\nWxEYGIiSJU3/FNWyC64ZGRnw8/NDuXLl0L59+3zt37lzJ18dPMfo0aNtinnDhg34z3/+g127dsHX\n1xeBgYFYuXIlvLy8EBISgsTERFy+fBlDbFyMs2vXroiPj1d1dMuhQ4ce6hdr0tOdb3y58dw5FPyX\nLgPAeP68FuHoVpcOXZBwMAGdX+qsdSi6wne8KqBTp07YunVr7uN169aha9euOHDgQO5ziYmJ6NKl\nC2bMmAE/Pz/07dsX7733Xu7NPIGBgbh16xYAICkpCYcOHcKGDRvQLHtFiHXr1uXbZ2ZmJipVqlTs\nmD08PBAVFQVAfps7e/YsateuDS8vL7Rr1w4tW7bEq6++ioyMjHwfVoXx9PRUffhinTp1UKFChSK9\nJ88/SZFpNYLCIygItwo8dwuAR5UqWoQDQLu+sEQIgUnRk/Jd2FeDHvsiL56gTCFfffUVrly5gvDw\ncNSsWRMNGjTAggULcPv2bWRlZaFEiRJ499138cknn6BKlSrw9fXFpUuXMHLkSADA/PnzERYWhqio\nKMyYMQM3b95E5cqVkZSUhOeffx5BQUGoV69e7v4++OADjBgxAkFBQcWO+euvv8aDBw+Qnp6OGjVq\nYPDgwbhz5w4+++wzNGzYEPfu3YOvr6/LjD0XQuDjjwkTJ2odSdGkp6biq5YtEXPqFMpAJvjo8HAM\n3bQJIdlDSZl7szRBmSIXXpX8gRNeeFVCZmamzfVvIqKBAwc6MBrXBIDaty/++7WsvaalpNCzPr1p\nZL0omtC7t+YXXfVQh9YLPfQFHFWTZ8opX748/P39cfXqVatDAvfv3+8yZ9dqs6dco6VHvMOQUmIp\nEn4HPLjIyoqADxcdGT58OFauXGlxm6ysLGzduhU99LqCuM7dvQtculS892pZe92zB3j+ef0keL3X\nodWk977QySHDAHkxdPDgwRa3uXz5MoYNG6ZSRK7n2Wfl3CXOZtcuoHFjraNgzoiTvJOpVKkSvFx1\ndQMV2JPktZw3fNcuoFEjzXb/EL3Poa4mvfcFJ3nmVp591vnq8nfvAocOAQ0aaB0Jc0Y8hJK5DSEE\njh4ldOgAnDypdTS227ULeO89ICFB60iYXvEar4xle/xx4OJFubScs9BbqYY5F07yzK2UKCGXk0tM\nLPp7taq97t6tv4uueq9Dq0nvfcFJnrmdunWdZ4QNkT6TPHMeXJNnbiPnGPruO2DHDmDRIq0jsi45\nWa6AdPq01pEwPeOaPGN5ONNYeR4fz+zFSZ65naefBk6cAO7cKdr7tKi96rVUo/c6tJr03hec5FV2\n5swZ+Pn52VR6Ksq2zHalSwM1awJHjmgdiXU8sobZi2vyrFh+/PFHTJ8+HYmJiWjQoEG++fT1Ku8x\nNGCAnAvmzTc1DsqCzEy5wElmJlCSpxJkFliqyTv1ofN/I0bg5oED+RYJICL41K2LUdOmqdaGO/L3\n98eIESOQlJTkFAm+IGe483XPHqB+fU7wzE6FzUGs1Q+KMJ987I8/0kZv73xrX8Z6e9PGlbav8ahE\nG6GhofR///d/9Mwzz5CPjw+98cYbdOnSJWrbti35+vpSy5YtcxebTktLIyEEZWVlERFRZGQkjR8/\nnho3bky+vr7UunVrunr1aqHbfvzxx9SoUSPy8fGhjh070tWrV6l3797k5+dH9evXp/T0dLPvzXn/\nvHnziIho4cKF1LhxYxoxYgSVK1eOwsPDaffu3bRw4UIKDg6mwMDAfAuNF+a7776jqKgom/tKS3mP\noe3bierXL9r71Z43/KOPiD75RNVd2kwPc6jrhR76Ahbmk3fqmnzrLl2wsXZt5BRxCMCvtWujVWfb\n13hUog0AWLVqFbZs2YLk5GSsWbMG7dq1w6RJk3DlyhVkZWVhxowZudsWXJ5s+fLlWLRoES5fvox7\n9+5h6tSphW67YsUKLFu2DOfPn8fJkyfRqFEjDBw4ENeuXUOtWrUQExNT6HsL2rdvHyIiIpCZmYle\nvXqhZ8+eSEhIwKlTp7BkyRIMGTIEt2/fLlI/OIuICFmTf/BA60gKx/V4pgSnTvJCCLT+4AP85u0N\nAPjV2xttRo0q0hqPSrQBAEOHDkVAQAAqV66Mpk2bokGDBnjmmWfg6emJV155BQctjNkbMGAAwsPD\nUbp0aXTv3h2JFm7HHDBgAEJDQ+Hr64u2bdsiPDwcUVFR8PDwQLdu3Szup6CwsDD069cPQgj06NED\nZ8+eRXR0NEqVKoWWLVvC09MTJ51pkpci8PUFgoKA48dtf4+a84bfvw/88QfQsKFquywSvc+hria9\n94VTJ3kg/5l4cc7AlWojMDAw93cvL6+HHt+8ebPQ9+ZdkNvb29vitvbsx1pbABAQEFDs9pyNnuvy\nBw8C4eFA2bJaR8KcndMn+Zwz8ZG+vsU6A1eqDb0pU6YMAOQrt1y8eFGrcHSpqNMbqDkeWu+lGr2P\nDVeT3vvC6ZM8IM/EK7/zTrHOwJVsw1ZUhKGgRdk2r4CAAAQFBWHp0qUwGo2YP38+Tp06pdi+jEYj\n7t27h/v37yMrKwv37t3DAz0XuM3Q852ver0Jijkfl0jyQgiMnjTJrjNwe9oo+B5rbeR9XcltC/r2\n228xZcoUBAQE4NixY2hsJWsU5e9YsmQJvLy88O6772Lnzp3w9vbGm3oedG5GTpK39bNNrdorkf6n\nM9B7HVpNeu8LvhmKuQ1zx1DVqsD27cBjj2kUlBmpqbJUc/484AKVQ6YCnqCMsUIUpS6vVu01p1Sj\n5wSv9zq0mvTeF5zkmVurW1cOVdQTvZdqmHPhcg1zG+aOofXrgenTgU2bNArKjDp1gLlzeeFuZjtL\n5RpO8sxtmDuGMjLkjJSZmfooj9y4AVSpIuPx9NQ6GuYsuCbPWCEqVgT8/ABbbuxVo/b6++/Ac8/p\nP8HrvQ6tJr33BSd55vbq1QP27dM6CknvN0Ex58NJnrm9evWA/futb6fGeGhnuQlK72PD1aT3vrAr\nyQshugohjgghsoQQdS1slyaEOCSEOCiE0Mk5E2OSrUne0R48APbulYuZMKYUe8/k/wTwCoBtVrYz\nAogkomeJqL6d+2RMUc89Bxw6ZH3aYUfXXv/8U86M6e/v0N0oQu91aDXpvS/sSvJEdJyITgCwNi5B\n2LsvZxQTE4N+/frZ1ca2bdsQHBxscZv4+Hi88MILKFeuHB7T062bTqJsWXnn619/aRsHj49njqBW\n4iUAm4QQ+4UQg5RsOD01FTF9+iA6KgoxffogPTVVkzYcISsrC0Rkdc6aMmXKYODAgfkWG2FFY0vJ\nxtG1V2epxwP6r0OrSe99YXX1SCHEJgCBeZ+CTNrjiGitjftpTEQXhBAVIJP9MSLaWdjG/fv3R2ho\nKACgXLlyiIiIMLtdemoqvmrZEjGnTqEMgFsAon//HUM3bUJIWJhNgSnRxuTJk/HVV1/hxo0bCAoK\nwuzZs2EwGPDZZ58BAH7++WdUr14dBw8exMKFCzFlyhScPXsWFStWxOjRo3Mn9tq2bRv69OmDoUOH\nYtq0aWjWrBnWrl0Lg8EAX19fCCGQnJycb/55AKhXrx7q1auHLVu22BSvO8v5ap3zP2bO43r1IrF/\nP1C9uvnX1Xi8axfQpk084uO12T8/dp7HOb+npaXBqsLWBSzKD4A4AHVt3DYawEgLr1tawzCfCb17\n0808a7MSQDcBmtC7t6XlEBVt4/jx4xQcHEwXL14kIqL09HRKSUmRbU+YQH379s23/YYNGyg1NZWI\niLZv307e3t508OBBIiKKj4+nkiVL0kcffUQGg4Hu3r1L8fHxFBwcbFMsmzdvprCwMJu2dUeFHVtE\nRLt3Ez37rOX3O3ItzzNniAICiIxGh+1CUXpY11Qv9NAXUGmNV/O31ArhLYTwyf69DIBWAI4osUPj\nuXMoU+C5MgCM58+r1kaJEiVgMBhw5MgRPHjwANWqVUOYhW8Abdu2zf2W0rRpU7Rq1Qo7duzI115M\nTAxKlSqF0qVL2/x3MPtERABJScDdu9rsf/duOT5eD3fdMtdi7xDKl4UQZwA0BLBOCBGb/XxlIcS6\n7M0CAewUQhwE8DuAtUT0mz37zeERFIRbBZ67BcCjShXV2ggPD8f06dMxYcIEBAYG4tVXX7W4AlNs\nbCyef/55+Pv7o3z58oiNjcWVK1dyX69QoQJKlSplc/xMGV5eQK1agIXldR1ae3W2i656r0OrSe99\nYe/omtVEFExEXkRUmYjaZj9/gYheyv49lYgiSA6frE1Ek5QIHAD6T5yI6PDw3CR9C0B0eDj6T5yo\nahs9e/bEjh07kJ6eDgAYM2YMgIcX3TAYDOjatStGjx6Ny5cv49q1a2jbtm2++VSKugAJU46W4+X5\nTlfmKE49rDEkLAxDN23C1N69ER0Vham9exfpgqkSbSQnJyMuLg4GgwGenp7w8vKCh4fs1sDAQKSl\npeUmcYPBAIPBgICAAHh4eCA2Nha//Wb5S01gYCCuXr2KGzduFLoNEeHevXswGAz5luVjRWMtyTtq\nPPStW8CxY8C//uWQ5h1C72PD1aT3vrA6ukbvQsLCEL10qWZt3Lt3Dx9++CGSkpJQqlQpNGrUCHPn\nzgUAdOvWDUuXLoW/vz8ee+wxJCQk4Msvv0S3bt1gMBjQoUMHdOrUyWL7NWvWRK9evfDYY4/BaDTi\n6NGjD42u2b59O6KionLP+r29vdG8eXNs3bq1WH+Tu2rQAJgyRf397tol57V/5BH1981cH081zNyG\ntWPIaJR3mx47BhT4HHWojz4CSpUCPv1UvX0y18JTDTNmAw8PefFzZ6F3cDhGXBwQFaXuPpn74CTP\nWB5NmwJ5RrTm44ja6z//AEeOON+kZHqvQ6tJ733BSZ6xPJo0UfdMfscOecGX6/HMUbgmz9yGLcfQ\nvXuyLn/+vFwxytFGjQJ8fYFPPnH8vpjr4po8YzYqXVpOPbxnjzr743o8czRO8owVUFhdXuna6/Xr\nwPHjQH0nXGFB73VoNem9LzjJM1aApYuvStq+HWjYUH57YMxRuCbP3Iatx9CNG0CVKsDVq45NwCNG\nABUqAGPHOm4fzD1wTV5Hzpw5Az8/P5uSTVG2Zcrx8wNq1gT++MOx++F6PFMDJ3mVBQcH48aNGzZN\nPFaUbdU2atQoPP744yhbtiyefPJJLFmyROuQFGWuZKNk7fXqVSAlxbnmq8lL73VoNem9L5x67poR\nn4zAgfQD+ZIgEaFuSF1M+3Saam24Ix8fH6xfvx41atTAvn370KZNG9SoUQMNGzbUOjRFNGkCLFwI\nZE8oqrht2+TdtTyrNHM0pz6Tb/yvxkgokYBtYdtyfxI8EtCkXhNV2wgLC8PUqVNRp04d+Pr6YtCg\nQcjIyEC7du3g5+eHVq1a4e+//wYApKenw8PDA0ajEQAQFRWFTz75BE2aNIGfnx/atGmDzMzMQrcd\nP348GjduDF9fX3Tq1AmZmZno06cPypYtiwYNGuD06dNm35vz/vnz5wMAFi1ahCZNmmDkyJEoX748\nqlevjj179mDRokWoVq0aKlWqhMWLFxf6N0dHR6NGjRoAgPr166Np06bYo9a4QxU0bSonDsvTfYrO\nG75xI9CypWLNqU7vc6irSe994dRJvkuHLqj9T2254iwAEFD7Zm10fqmzqm0AwKpVq7BlyxYkJydj\nzZo1aNeuHSZNmoQrV64gKysLM2bMyN22YPll+fLlWLRoES5fvox79+7lW5C74LYrVqzAsmXLcP78\neZw8eRKNGjXCwIEDce3aNdSqVQsxMTGFvregffv2ISIiApmZmejVqxd69uyJhIQEnDp1CkuWLMGQ\nIUNw+/Ztq3/7nTt3sH//fjz11FNWt3UWgYHy5+BB5ds2GoG1awErE5AypginTvJCCHzQ9wN4n/YG\nAHine2NUv1FFqmEr0QYADB06FAEBAahcuTKaNm2KBg0a4JlnnoGnpydeeeUVHLSQLQYMGIDw8HCU\nLl0a3bt3R6KF5YkGDBiA0NBQ+Pr6om3btggPD0dUVBQ8PDzQrVs3i/spKCwsDP369YMQAj169MDZ\ns2cRHR2NUqVKoWXLlvD09MTJkyettvPWW2/h2WefRatWrWzetzNo3x5Yt870WKnaa0IC8OijQHi4\nIs1pQu91aDXpvS+cOskD+c/Ei3MGrlQbgYGBub97eXk99PjmzZuFvjfv/PDe3t4Wt7VnP9baAoCA\ngIAitTdq1CgcPXoUK1assHm/zqJDh/xJXilr1gAdOyrfLmPmOH2SzzkT943zLdYZuFJt6E2ZMnJ5\n8rzlFktrzxZHdHQ0fv31V2zatAk+Pj6Ktq0HjRsDp07JeWwA5WqvrpDk9V6HVpPe+8Lpkzwgz8Tf\neeGdYp2BK9mGrYoy7r24Y+QDAgIQFBSEpUuXwmg0Yv78+Th16pRi+/r888+xfPlybN68GeXKlStW\njHpXqhTQujWwYYNybaamApcuOedUBsw5uUSSF0JgUvQku87A7WmjqItv531dyW0L+vbbbzFlyhQE\nBATg2LFjaNy4sc37sra/cePG4cyZM6hevTp8fX3h5+eHSZMUW6NdN156SV4kBZSpva5dK2v9JUrY\n3ZSm9F6HVpPe+4KnNWBuozjHUGYmEBoqz7737o23+6v5iy8CQ4YAL79sVzOai4+3vy9chR76wtK0\nBpzkmdso7jHUrBnw4YdAu3b27f/6daBaNeDCBSD7kgljiuC5axizg1KjbDZulB8YnOCZmjjJM2bF\nSy/JJB8XF29XO64wqiaH3uvQatJ7X3CSZ8yKWrXkSBsrg5MsunYNiI3lu1yZ+rgmz9yGPcfQiBGA\njw8wcWLx9j1tGrB/P/D998V7P2OW8IVXxmDfMXT0KPDCC0BaGvDII0V7r9Eovw0sWCBvsGJMaXzh\nlTE7PfkkUK1aPH74oejv3bwZ8PYGGjVSPi6t6L0OrSa99wUnecZs1LUr8OWXQFG/DMyaBbz7LuAC\ns2UwJ8TlGuY27D2GjEZ5Rj9nDtC8uW3vSU8H6tYFTp/moZPMcbhco5GYmBj069fPrja2bduG4OBg\ni9tMnToVtWvXhp+fH8LDw/PNR8+U4+EBvPceMH267e+ZMwfo25cTPNOOXUleCDFFCHFMCJEohPhJ\nCOFXyHZthBBJQohkIYSiC6qlpqWiz7A+iOofhT7D+iA1LVWTNhwhKysLRGTTnDVLlizB9evXERsb\ni5kzZ+J///ufChG6l/j4ePTrJ9d+TUmxvv2dO8C8ecA77zg+NrXpvQ6tJt33BREV+wfAiwA8sn+f\nBOBzM9t4ADgJIARAKQCJAGpZaJPMMfd8SmoKhbcPJ4wFYQIIY0Hh7cMpJTXFbBvmKNHGpEmTKCgo\niHx9falWrVq0detW2rhxI3l6epKnpyf5+PhQREQEEREtWLCAnnjiCfL19aXw8HCaM2dObjvx8fFU\ntWpVmjx5MlWqVIm6d+9OXl5eVKJECfLx8SFfX1+6cOGC1XiGDRtGw4YNszl+d1HYsWWruLg4IiIa\nPZrovfesb//WW0Q9e9q1S93K6Qumj77IPrbN59TCXijqD4CXASwx83xDALF5Hn8IYIyFdiz9Efn0\nHtrblJwnmJJ076G9be4ce9s4fvw4BQcH08WLF4mIKD09nVJS5AfEhAkTqG/fvvm237BhA6WmphIR\n0fbt28nb25sOHjxIRDLJlyxZkj766CMyGAx09+5dio+Pp+DgYJv/HiKiZ599Nt+HB5PsTfI5zpwh\nqliRKDa28G2WLSOqXp3o+nVFdsmYRZaSfEkFvxS8DsDcALMgAGfyPD4LQJHZtM/dOAf4F3jSEzh/\n47xqbZQoUQIGgwFHjhyBv78/qlWrZnH7tm3b5v7etGlTtGrVCjt27EBERERuezExMShVqpTNf0Ne\n0dHRICIMGDCgWO9n1lWtCvz8s5xJsm+LEfA5fyBfSe3WLcKqI3Xx055pKFtWw0AZA6wneSHEJgCB\neZ+CXPZ6HBGtzd5mHID7RKTI/Xz9+/dHaGgoAKBcuXK5CbCgIL8gwADAM8+TBqCKXxWb92VvG+Hh\n4Zg+fTomTJiAo0ePonXr1vjvf/+bb0m/vGJjY/Hpp58iOTkZRqMRd+7cwTPPPJP7eoUKFYqd4GfO\nnImlS5di586dxW7D1eXUT3Omhi3K47y118jISMybB7zeqxzev78XHxruye0BbEdpvDJ4GCIi7Nuf\nnh8X7BOt49HycWJiIoYPH67q/nN+T0tLg1WFneLb+gOgP4BdAEoX8npDABvzPFasXKOXmnyOf/75\nh3r16kX9+vUjIqKYmJh85Zp79+6Rt7c3rVq1irKysoiI6OWXX6bx48cTEZktzWzbts2mcs28efMo\nODiY0tLSihy3uyjs2LKVudrrd98ZqZlnAzLK4fNkBKhDhQaUlWW0a196p4c6tF7ooS9goVxj7+ia\nNgBGAehIRPcK2Ww/gOpCiBAhhCeAngDW2LPfHGGhYdg0cxN6/9MbUalR6P1Pb2yauQlhoWGqtZGc\nnIy4uDgYDAZ4enrCy8sLHh6yWwMDA5GWlpY7NttgMMBgMCAgIAAeHh6IjY3Fb7/9ZrH9wMBAXL16\nFTdu3Ch0m2XLlmHcuHHYtGkTQkJCbPzLWVGZWxhi4ECBNqM+wIbS3gCAjd7eePfrUfDwcO07n7Re\nJENPdN8XhWV/W34AnACQDuBA9s/s7OcrA1iXZ7s2AI5nb/+hlTYtfVLpzuHDh6l+/frk5+dH/v7+\n1KFDh9wRMFevXqUmTZpQ+fLl6bnnniMiolmzZlFgYCCVL1+e+vXrR7169bJ4Jk9ENHDgQPL396fy\n5cubHV0TFhZGnp6e5OvrmzsK5+2333bgX+2cHHUMGY1GGt5Ans0Pb9CAjEbXPotn+gMLZ/J8xytz\nG/YeQ5aWedu4ciV+ff11tFmwAK27dCn2PpyFHpa80ws99IWlO16VHF3DmNtq3aULDickoFXnzlqH\nwlg+fCbP3AYfQ8xV8dw1jDHmpjjJM2Yj3c9RoiLuCxO99wUnecYYc2Fck2dug48h5qpcYnRNSEiI\nTVPuMlYYvlGMuSOnKdfk3Dmqxk9cXJxq+9L7jyv1hU3zfFig99qrmrgvTPTeF06T5NWUmJiodQi6\nwX1hwn1hwn1hove+4CRvxvXr17UOQTe4L0y4L0y4L0z03hec5BljzIVxkjfD3tqtK+G+MOG+MOG+\nMNF7X+hyCKXWMTDGmLOhQoZQ6i7JM8YYUw6XaxhjzIVxkmeMMRfGSd4KIcT7QgijEOJRrWPRihBi\nihDimBAiUQjxkxDCT+uY1CSEaCOESBJCJAshxmgdj1aEEFWFEFuFEH8JIf4UQgzTOiatCSE8hBAH\nhBCKLGnqCJzkLRBCVAXQEnKJQ3f2G4CniCgCcgnHjzSORzVCCA8AMwG0BvAUgF5CiFraRqWZBwBG\nEtFTAJ4H8K4b90WO9wAc1ToISzjJWzYNcqFyt0ZEm4nImP3wdwBVtYxHZfUBnCCidCK6D+AHAJ00\njkkTRHSRiBKzf78J4BiAIG2j0k72SWA7AN9pHYslnOQLIYToCOAMEf2pdSw68zqAWK2DUFEQgDN5\nHp+FGye2HEKIUAARAPZqG4mmck4CdT1E0WlmoXQEIcQmAIF5n4L8B/sYwFjIUk3e11yWhb4YR0Rr\ns7cZB+A+EX2vQYhMJ4QQPgBWAngv+4ze7Qgh2gO4RESJQohI6Dg/uHWSJ6KW5p4XQjwNIBTAISHn\nN64K4A8hRH0iylAxRNUU1hc5hBD9Ib+avqBKQPpxDkC1PI+rZj/nloQQJSET/BIi+kXreDTUGEBH\nIUQ7AF4AfIUQi4mon8ZxPYRvhrKBECIVQF0iuqZ1LFoQQrQB8AWAZkR0Vet41CSEKAHgOIAWAC4A\n2AegFxEd0zQwjQghFgO4QkQjtY5FL4QQzQG8T0QdtY7FHK7J24ag469jKvgKgA+ATdnDxWZrHZBa\niCgLwBDIEUZ/AfjBjRN8YwC9AbwghDiYfSy00TouZhmfyTPGmAvjM3nGGHNhnOQZY8yFcZJnjDEX\nxkmeMcZcGCd5xhhzYZzkGWPMhXGSZ4wxF8ZJnjHGXBgnecYsEEL8SwhxSAjhKYQoI4Q4IoR4Uuu4\nGLMV3/HKmBVCiE8hJ6Hygpx+erLGITFmM07yjFkhhCgFYD+AOwAaEf9Pw5wIl2sYsy4AcoI2XwCP\naBwLY0XCZ/KMWSGE+AXAcgBhAKoQ0VCNQ2LMZm69aAhj1ggh+gIwENEP2Yt67xJCRBJRvMahMWYT\nPpNnjDEXxjV5xhhzYZzkGWPMhXGSZ4wxF8ZJnjHGXBgnecYYc2Gc5BljzIVxkmeMMRfGSZ4xxlzY\n/wMt18pkrenrOQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from scipy import arange, cos, exp\n", "from scipy.optimize import fmin\n", "import pylab\n", "\n", "def f(x):\n", " return cos(x) - 3 * exp( -(x - 0.2) ** 2)\n", "\n", "# find minima of f(x),\n", "# starting from 1.0 and 2.0 respectively\n", "minimum1 = fmin(f, 1.0)\n", "print(\"Start search at x=1., minimum is\", minimum1)\n", "minimum2 = fmin(f, 2.0)\n", "print(\"Start search at x=2., minimum is\", minimum2)\n", "\n", "# plot function\n", "x = arange(-10, 10, 0.1)\n", "y = f(x)\n", "pylab.plot(x, y, label='$\\cos(x)-3e^{-(x-0.2)^2}$')\n", "pylab.xlabel('x')\n", "pylab.grid()\n", "pylab.axis([-5, 5, -2.2, 0.5])\n", "\n", "# add minimum1 to plot\n", "pylab.plot(minimum1, f(minimum1), 'vr',\n", " label='minimum 1')\n", "# add start1 to plot\n", "pylab.plot(1.0, f(1.0), 'or', label='start 1')\n", "\n", "# add minimum2 to plot\n", "pylab.plot(minimum2,f(minimum2),'vg',\\\n", " label='minimum 2')\n", "# add start2 to plot\n", "pylab.plot(2.0,f(2.0),'og',label='start 2')\n", "\n", "pylab.legend(loc='lower left')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calling the `fmin` function will produce some diagnostic output, which you can also see above.\n", "\n", "##### Return value of `fmin`\n", "\n", "Note that the return value from the `fmin` function is a numpy `array` which – for the example above – contains only one number as we have only one parameter (here *x*) to vary. In general, `fmin` can be used to find the minimum in a higher-dimensional parameter space if there are several parameters. In that case, the numpy array would contain those parameters that minimise the objective function. The objective function $f(x)$ has to return a scalar even if there are more parameters, i.e. even if $x$ is a vector as in $f(\\mathbf{x})$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Other numerical methods\n", "-----------------------\n", "\n", "Scientific Python and Numpy provide access to a large number of other numerical algorithms including function interpolation, Fourier transforms, optimisation, special functions (such as Bessel functions), signal processing and filters, random number generation, and more. Start to explore `scipy`’s and `numpy`’s capabilities using the `help` function and the documentation provided on the web." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "scipy.io: Scipy-input output\n", "----------------------------\n", "\n", "Scipy provides routines to read and write Matlab `mat` files. Here is an example where we create a Matlab compatible file storing a (1x11) matrix, and then read this data into a numpy array from Python using the scipy Input-Output library:\n", "\n", "First we create a mat file in Octave (Octave is \\[mostly\\] compatible with Matlab):\n", "\n", "```octave\n", "octave:1> a=-1:0.5:4\n", "a =\n", "Columns 1 through 6:\n", " -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000 \n", "Columns 7 through 11:\n", " 2.0000 2.5000 3.0000 3.5000 4.0000\n", "octave:2> save -6 octave_a.mat a %save as version 6\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we load this array within python:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from scipy.io import loadmat\n", "mat_contents = loadmat('scipy/code/octave_a.mat')" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{'__globals__': [],\n", " '__header__': b'MATLAB 5.0 MAT-file Platform: posix, Created on: Mon Aug 8 12:21:36 2016',\n", " '__version__': '1.0',\n", " 'a': array([[-1. , -0.5, 0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. ]])}" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mat_contents" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-1. , -0.5, 0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. ]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mat_contents['a']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The function `loadmat` returns a dictionary: the key for each item in the dictionary is a string which is the name of that array when it was saved in Matlab. The key is the actual array.\n", "\n", "A Matlab matrix file can hold several arrays. Each of those is presented by one key-value pair in the dictionary.\n", "\n", "Let’s save two arrays from Python to demonstrate that:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import scipy.io\n", "import numpy as np\n", "\n", "# create two numpy arrays\n", "a = np.linspace(0, 50, 11)\n", "b = np.ones((4, 4))\n", "\n", "# save as mat-file\n", "# create dictionary for savemat\n", "tmp_d = {'a': a,\n", " 'b': b}\n", "scipy.io.savemat('data.mat', tmp_d)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This program creates the file `data.mat`, which we can subsequently read using Matlab or here Octave:\n", "\n", " HAL47:code fangohr$ octave\n", " GNU Octave, version 3.2.4\n", " Copyright (C) 2009 John W. Eaton and others.\n", " \n", "\n", " octave:1> whos\n", " Variables in the current scope:\n", "\n", " Attr Name Size Bytes Class\n", " ==== ==== ==== ===== ===== \n", " ans 1x11 92 cell\n", "\n", " Total is 11 elements using 92 bytes\n", "\n", " octave:2> load data.mat\n", " octave:3> whos\n", " Variables in the current scope:\n", "\n", " Attr Name Size Bytes Class\n", " ==== ==== ==== ===== ===== \n", " a 11x1 88 double\n", " ans 1x11 92 cell\n", " b 4x4 128 double\n", "\n", " Total is 38 elements using 308 bytes\n", "\n", " octave:4> a\n", " a =\n", "\n", " 0\n", " 5\n", " 10\n", " 15\n", " 20\n", " 25\n", " 30\n", " 35\n", " 40\n", " 45\n", " 50\n", "\n", " octave:5> b\n", " b =\n", "\n", " 1 1 1 1\n", " 1 1 1 1\n", " 1 1 1 1\n", " 1 1 1 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that there are other functions to read from and write to in formats as used by IDL, Netcdf and other formats in `scipy.io`.\n", "\n", "More → see [Scipy tutorial](http://docs.scipy.org/doc/scipy/reference/tutorial/io.html)." ] } ], "metadata": { "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.4.4" } }, "nbformat": 4, "nbformat_minor": 1 }