{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "
\n", " \n", " \"QuantEcon\"\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# General Purpose Packages" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Contents\n", "\n", "- [General Purpose Packages](#General-Purpose-Packages) \n", " - [Overview](#Overview) \n", " - [Numerical Integration](#Numerical-Integration) \n", " - [Interpolation](#Interpolation) \n", " - [Linear Algebra](#Linear-Algebra) \n", " - [General Tools](#General-Tools) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overview\n", "\n", "Julia has both a large number of useful, well written libraries and many incomplete poorly maintained proofs of concept.\n", "\n", "A major advantage of Julia libraries is that, because Julia itself is sufficiently fast, there is less need to mix in low level languages like C and Fortran.\n", "\n", "As a result, most Julia libraries are written exclusively in Julia.\n", "\n", "Not only does this make the libraries more portable, it makes them much easier to dive into, read, learn from and modify.\n", "\n", "In this lecture we introduce a few of the Julia libraries that we’ve found particularly useful for quantitative work in economics.\n", "\n", "Also see [data and statistical packages](data_statistical_packages.html) and [optimization, solver, and related packages](optimization_solver_packages.html) for more domain specific packages." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setup" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "hide-output": true }, "outputs": [], "source": [ "using InstantiateFromURL\n", "github_project(\"QuantEcon/quantecon-notebooks-julia\", version = \"0.5.0\")\n", "# github_project(\"QuantEcon/quantecon-notebooks-julia\", version = \"0.5.0\", instantiate = true) # uncomment to force package installation" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "hide-output": true }, "outputs": [], "source": [ "using LinearAlgebra, Statistics\n", "using QuantEcon, QuadGK, FastGaussQuadrature, Distributions, Expectations\n", "using Interpolations, Plots, LaTeXStrings, ProgressMeter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numerical Integration\n", "\n", "Many applications require directly calculating a numerical derivative and calculating expectations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adaptive Quadrature\n", "\n", "A high accuracy solution for calculating numerical integrals is [QuadGK](https://github.com/JuliaMath/QuadGK.jl)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "hide-output": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(value, tol) = quadgk(cos, -2π, 2π) = (-1.5474478810961125e-14, 5.7846097329025695e-24)\n" ] } ], "source": [ "using QuadGK\n", "@show value, tol = quadgk(cos, -2π, 2π);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is an adaptive Gauss-Kronrod integration technique that’s relatively accurate for smooth functions.\n", "\n", "However, its adaptive implementation makes it slow and not well suited to inner loops." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gaussian Quadrature\n", "\n", "Alternatively, many integrals can be done efficiently with (non-adaptive) [Gaussian quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature).\n", "\n", "For example, using [FastGaussQuadrature.jl](https://github.com/ajt60gaibb/FastGaussQuadrature.jl)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "hide-output": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "w ⋅ f.(x) = 0.6666666666666667\n" ] } ], "source": [ "using FastGaussQuadrature\n", "x, w = gausslegendre( 100_000 ); # i.e. find 100,000 nodes\n", "\n", "# integrates f(x) = x^2 from -1 to 1\n", "f(x) = x^2\n", "@show w ⋅ f.(x); # calculate integral" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The only problem with the `FastGaussQuadrature` package is that you will need to deal with affine transformations to the non-default domains yourself.\n", "\n", "Alternatively, `QuantEcon.jl` has routines for Gaussian quadrature that translate the domains." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "hide-output": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "w ⋅ cos.(x) = -3.0064051806277455e-15\n" ] } ], "source": [ "using QuantEcon\n", "\n", "x, w = qnwlege(65, -2π, 2π);\n", "@show w ⋅ cos.(x); # i.e. on [-2π, 2π] domain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Expectations\n", "\n", "If the calculations of the numerical integral is simply for calculating mathematical expectations of a particular distribution, then [Expectations.jl](https://github.com/QuantEcon/Expectations.jl) provides a convenient interface.\n", "\n", "Under the hood, it is finding the appropriate Gaussian quadrature scheme for the distribution using `FastGaussQuadrature`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "hide-output": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "E(f) = -6.991310601309959e-18\n" ] }, { "data": { "text/plain": [ "true" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Distributions, Expectations\n", "dist = Normal()\n", "E = expectation(dist)\n", "f(x) = x\n", "@show E(f) #i.e. identity\n", "\n", "# Or using as a linear operator\n", "f(x) = x^2\n", "x = nodes(E)\n", "w = weights(E)\n", "E * f.(x) == f.(x) ⋅ w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpolation\n", "\n", "In economics we often wish to interpolate discrete data (i.e., build continuous functions that join discrete sequences of points).\n", "\n", "The package we usually turn to for this purpose is [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl).\n", "\n", "There are a variety of options, but we will only demonstrate the convenient notations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Univariate with a Regular Grid\n", "\n", "Let’s start with the univariate case.\n", "\n", "We begin by creating some data points, using a sine function" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "hide-output": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd2BUVdow8Ofcmcmk98ykAgk99JZGAgEERFSKCCKguCIosurrsuquDdwV1311/d5dG64uirgaFEFQKaEnpAEh9A4J6TOZtEmdcs/3x2gMSYCUmbnn3vv8/oKbZO6ZmXPPc/ohlFJACCGE5IoTOgEIIYSQkDAQIoQQkjUMhAghhGQNAyFCCCFZw0CIEEJI1jAQIoQQkjUMhAghhGQNAyFCCCFZw0CIEEJI1jAQIoQQkjVWAuHJkyeFTgITrFar1WoVOhVsMZvNQieBLVarled5oVPBEEqpxWIROhVssVgsuH1m57ESCEeOHIlfGwA0NDQ0NDQInQqG8DxfWVkpdCrYYjQam5ubhU4FQ8xmc01NjdCpYEt1dTVWDjqPlUCIEEIICQIDIUIIIVnDQIgQQkjWMBAihBCSNQyECCGEZM3OgdBqtT722GOtrxiNxldeeWXOnDmvvvqq0Wi07+0QQgihHrJnIPz++++fffbZoqKi1hdTUlK0Wm1KSopGo9m8ebMdb4cQQgj1nD0DYVRU1OLFi9tcTE9PnzVrlouLy6xZs9LS0ux4O4QQQqjnlHZ8rZEjR7a/aDAYtFotAGi12tuvjH7xxRcJIW0uDhs2bNasWXZMJOPq6upu9aMqEzlTTc7Xcl5KelcwH+Qqi/0HeJ6vq6tzd3cXOiEMqaurM5lMjC+XrreQwgYobyRlTcRK4e5Q3t/FUTnWZDLV1dW5uLg46PXFqK6ujuM4lUoldEKE5+7urlAobv879gyEHaKU2sIbpfT2+0K5uLi0D4RKpbL9RQmzvdnWb9lKYf0lxf+7oGyw0GhfGu1DK5q5F3KV/bzo3WH06YFWL6WUIyL5ldAJYQjjn0mDBd6/qPzgIqdxBY0rDXGHZiu8mKtM0tCHI63TQ3kXe0/RY/wDEQR+Ji068yE4PBAGBATodLrw8PCKiorAwMDb/OZf/vIX/Nps+8x5enra/nuqkj6RZnVXwp57FEP9fvtwzDykl9ONl/mJqcqUyYpRAZL93Hieb2xsbPlAEABYLBa1Wu3m5iZ0QtriKWy6wr9yjI/XkmNzuEiv37Kl0QzfXec/uqT43/PwwzRFhIc9c6zJZKKUYiZprampydPTE1uEneTA5RO2fbTj4uJ2795NKd29e3dCQoLjbic9a3KtU3dalg/i9s9Uto6CAKDiYFII2TBB8dcx3N27LB+dxy2YkcDMPDy4z/rheT5liiJlsqJ1FAQALxU8NoA7fK9yYV8u7gdrpk7K3RhIdBwYCFevXg0AS5YsuXbt2sKFC/Pz8xctWuS424ldbm7u9OnTQ0NDQ0NDp0+fvvyrYz/doCfnqh4fyN2m8jw/iku/V/nvC/zSQ1YsWpBQLDwsOmA183D4XmW85natvT8O5z5JUszaY/nyCtbeECsII2c+EEJ4npdt12hGRsa90+76c3zUvf2DAeDHy2Vrj1xL+XnvPcmdakM3WWHaTsvEEPKXMXcYExYdnuf1er1tvhWyqa6uZqpr1EphyUFrtYluvUup7lwGPFtF79tjXTuGW9LPDnVxk8lkNBoDAgJ6/lKSUVFR4ePjg12jnYSBkAlTpkyZRsoXDQ1vufLVmaI9VLtv375OvoK+CWJ/sKwbyz3UV1K7BWEgbI+pQMhTWHrIWt5If5imdO1KNexMFZ38kyX1HuUI/54+9RgI28NA2CWSKjTFKzs7e3qUpvWV6VGa7Ozszr9CkCv8MFXxbJb1qJ6Jmg2SiX+e5a/U0m1TuxYFAWCoH/m/eMW8vdZqk2NShlCnYSBkAsdxlpvXlph5nuO69u0M8yefJinm7rUW12MsRM5woZquy7N+maxw69b084V9uZm9yJKDFh4zLBIUBkImJCUl/Xi5vPWVny6XJyUldfV17uvFLR/EPXnEar+kIdQxCw+PHrK+MUbR17v7fZv/G6OoNsG6PJw4g4Tk8HWEqDPWrVuXPGGCldKZ/bUA8NPl8o9Ole479E03XuqlEdzI7/kfCvhZvbGWgxzorZO8nxpWDO5RNlNxkDJZMfJ7y7xIMshXplMEkOAwEDJhwJARPn89lrZ/zca96QCQlJSU8/kbUVFR3XgpFQcfJyqWHLTeFcZ54NeLHOOEgf7rrPX4HGXPY1eoO/nTSMXvM62pM5jLr++9997GjRuFTkV3WK1WjuMkNv1w0qRJ//jHPxzxyszlPHl67bh1/NCor37/le2kKi8vr568WlIwSQomb56wrhsntdUUiAU8hccPW/8RZ7cNYlZFc/8+b33tP9stV7IIIYmJiTNmzLDLK/dQcXHxpEmTlixZInRCEBw8ePDQoUMOenEMhMI7XUk3XeFPP2DPic7/G6sYvsX8SH8Ou5uQ3aVc410UsMgeSwBtGutqVe/N2nY1955+Wgqw+pN/vTcq7vvvv2dh17TQ0NBRo0YJnQoEBQUFjguEOIwkvNdz+RdGKAJd7fmawW7w6ijFk0dwuxlkZxYe1uby68Yq7FjDevnll3tVXPjxobhV46J+Py7q54fj4erJtWvX2u8OCN0OBkKBnTDQHD19cpD9v4iV0ZyhCXYVYihE9vTZJb6XJ0wOtWdPw3ffffd8XD/u1wEtBSHPx/X99ttv7XgLhG4DA6HAXj/Ovzic694yrNtTEHhlFPdmHi6lQHbTaIG/nuDfHGvPsWdKqV6vj/C+qUskwttNr9fb8S4I3QYGQiEdr6C5BrpsoKO+hXmRXEUTHCzFRiGyj3+d4+M0ZFyQPZuDhJDIyMizemPri2f1xu7NmkaoGzAQCunV49ZXRjqkOWijIPDiCGwUIvuoMcG7p61vjLF/obFixYrXDp03NP6y2VpFg2nNoQsrVqyw+40Q6hDOGhVMpo6eq4JtUx1bF1ncj3vjBJ+to7G3PRwHoTv611n+nghusAPmIT///POVlZVT3v+/sRpPnkJWWf3Tz//h6aeftvuNEOoQBkLBrM21vjKKc3Fwm1zFwQvDuTfz+O3TcE0h6j4TDx+d5/fMcEgu4jhu3bp1q1atOnbsGCHEyzQ6YECIxBaDI5ZhIBTG+Wp6qpJun+aMrunfDeDezLPkVtDRgViyoG769ho/xA+G+DkwC4WGht5///0A0LuS3r3Lsiqa6+qJFgh1D44RCuODc/yygQ5vDtqoFfD8UO7vp3BfY9R9/zzL/36Ik4qL4f5kZADZhEfY38nUqVPv+DsbNmx48MEHnXAjUcNAKACjGb6+yq9wwNrBW1k+iEst5ksacPoo6rL6+vpMHa1shpkRzsuxLw5X/O8pHo9nur0FCxbc8Xd27Njx6aefdvsWL7zwQidvJGoYCAWw8TI/JZQLs9M+jZ3hqYJ5kdyGS1iuoM6qqqpatWqVv79/kK/PlOiwgUfeNpuanXb3iSHExwV2FWGOvZ1ly5bd8XeMRqOPj0+3b3HixIlO3kjUMBAK4KPz/Cpn9TK1eGow98kFHrdcQ51hNpunTJmiP7Bt/7yRl56+a/e9/Wq3vbd06VJnpuGJgdxnF7F39Ca7du2aP3/+3Llzt27dCr/2WE6dOvXAgQPLly+fO3fuli1bWv/+a6+9BgArVqxo3bd5q7+qqalZu3btAw88sHTp0rS0tPZ/Xltb+/bbb8+fP3/BggVvv/12bW3t7e8uIjhZxtn2FlMASAp29ryVkQFE4wp7iuiMCJwyg+5gy5YtpPT6/86PsW17Funr/snMkRO+2HrmzJmhQ4c6Jw0L+3IvHTWXNihC3J1zw05pssJD+520MNdNCV9Pumm+0Mcff/zee++pVKp//vOfc+bMabmu0+nWr1+fl5f3yiuvPPDAAy3X33jjjalTp65fv77DQb42f/XRRx+FhIS88sor586de/3112NjY9v8+UcffaRUKr/88ksA+OCDD9avX//HP/7xNncXEQyEzvbhef73QzhBYtGTg7n1F/gZETgVD93B8ePHk/sEcq0WMHi6KGPCfI8fP+60QOipgrl9uC8u8y+NYKjjSsXB0gFOenyV7RaQDB8+/LPPPps6depbb73V+vp9991HCBk1apTJZOr867f5q5ycnM8//1yhUAwbNsz2jza/n52d/dlnn6nVagB47LHHli9f3pO7MwUDoVPdqKOHS/kvk+154lLnPdSXeyHHXFTPhTtxeBKJkUqlqrW0bfc0WXil0qklxrKB3OKD1hdHCFNx7JCCwOzeggXmtWvXHjt2bM+ePT/99NPf//73luvu7p1tNTc0NNzqr3ieb1m7WVVVpVKp3Nzc2vx568WdVqu1q3dnFkNVLTn49CK/uL9gB8d7KGFhX+6zizhOiO4gOTl511VdU6tYWGJsytUZk5KSnJmMWA1xV+Jmub9ZvHhxSEjI4sWLL1++3KU/VKlUeXl5lNIdO3bc6ndGjx69efNmq9V6/vz5VatWWSwW2/WWf8TExGzYsMFkMjU3N2/YsCE2Nrbbb4Q1GAidhwJ8dYU+NkDIz3z5IO7TizhlBt3BtGnTRkyavmjb8fRCQ2Ft466r5Q9vPbbq+T/26tXLySlZNpD79AJOmfnFggULnnnmmdWrV7d0S3bS0qVL165du3z5cj8/v1v9zqpVq65fv/7ggw+uW7du9erVXl5eABATE/PII4/YfmHlypXNzc2LFi1asmSJ2Wx+6qmnevJemEIoZaJQJIS0bphLUnoZfeqI9fQDt2sPGo1GALBlQQcZv8Py4nDufuG6d7qE53m9Xq/VaoVOCEOqq6vVanX7biv74nn+oTc+Tfvhv+61hQMGDFi1atXMmTMdescOVTVDVIr56gKVv/qWv2MymYxGY0BAgN3vvnr16uDg4NWrV9v9lVFXbdu27fPPP9+2bZsjXhzHCJ1n0xV+SX/hw8/S/tymK/T+3kKnA7GN47hLw3/3zYonJoYIWT31U8O9vbgvL/PPDhX+2UFShXnLSZqtsCWff7iv8E3eeZFcajFvNAudDsS2PAOtMcMEQaOgzbKB3IZL2DuKHAgDoZP8VMgP9ycsTNf0U8OEEG5rPpYs6HY2XuYf6cfEWMWEEFJtgrNVTAziIEnCQOgkm67QJf1Y+bQXRpGvr2IgRLdk4eGba/wiNnIsAZgXSTZfk0KObW5uFu9iOwljIqNLXlUz7C/h5/Rh5dO+vzeXpaPljUKnA7FqVxHt600G+LDQIAQAWBDFpVwTcYuwqqrq5ZdfjoqK8vX08HJ3Gzx48FtvvVVfXy90utAvWCmapW3zdX5GBOfjInQ6fuWuhHt7cd9dl0IVGznCxsv8I2w0B23GBhELhRMGUcbCs2fPjhk1smbP119ODL/89F1XV039YGzAxZSP42Jjb9y4IWDCOn+40u1/UwKHNDGU1yXsy8v8YpaKFQBIVl7/2+t/XrBgwerVq48fPy50chBDGiywp5h/IJKhHEsAFkSRFBH2jlZXV8+6d+ZbY7QvjR/Q2+eXHVgGBHj+NXnw072Vs+6/X9o9pbZTnNjHUF6XqnwjvVxLp4Wx0ssEAJs2bXpp5qjZBd8mGy8qMrffk5z45ptvCp0oxIqfC/nYIHKbdXuCsPWOiq5J+O67784IUoyP8G//o3v7B/c3V3z22WfOT5XT2E5xYh+uI3S4Lfl0dm9OxUyVo6Ki4vdPr/zqnmHDtd62K/MGh9697i/33Xff8OHDhU0bYsF31+k8lpqDNsP9iYcSsnU0TsNQnfKOUlJSPo0PvdVPF0SH/SslpfNbtOzates///mPxWJZsmTJnDlzMjIyvvjii/LychcXl3nz5s2fP3/q1KkPPfTQ7t2777vvPkrp/v37q6qqFixY8PDDD0+dOnXZsmXffvttRETEiy++GBwcbHvNysrK999//9SpU25ubkOHDl22bFlAQEBtbe0HH3xw9OhRT0/P1sdc2HT40/aJaTnF6dFHH23zoy5/jg7GXHaXnq35DE2TAYB9+/aN8FO3REEACPF0vbef9jabECL5aLLCnmJ+FpMbD82P4sTVO2qxWAquX+vje8s9qfv6eVy4cKHzL/jxxx+//fbb//znPzMzMwHg888/nzJlypYtW9atW7dhwwbb70RGRr799tsbN2708PD497///dprr23atMn2o8bGxpSUlKFDh3700Uctr/mPf/xjwoQJX331le0Ypvfeew8APvzwQwD48ssv169f335f0w5/2j4xb7zxBgCsX7++w3QyBVuEjlXWCOeq6eRQhuqwNTU1Gve283a0Hi41NTWCpAcxZU8RP8KfaBy7fVs3LYgik3/i340Fhk6juC2O4win4OktE2zheZWqC2fRtDmG6eOPP7548WJqaurJkydbtsZOSkqyHRIye/ZshUIxevRos/mX7TOmT5+uUCjmzZv32GOPtbzmyZMns7OzW/5rO84+Jyfn008/9fDwAIBly5alpqa2TkaHP+0wMTa3+REjWKz3ScnWfH5mBOfC0sc8cODAE+U1/M2jLbllNQMHDhQqSYgdbPaL2gz0IVo3OFwmmoFCjuMGDR580VB3q184X1E3bNiwzr/g2rVrZ82alZ6e/qc//QkA/vKXv2zdutXHx+fxxx9v+R2VSmXbtNl2oGD7DZxtGzu3/NfT03Pjxo2pqampqanbt2+3tfZa/1WHr9D+3x0m5o4/YgSjOV4ytubzc/qwVX1NTEz0jRr81/RLZisPABRgw8kbF0wuDz74oNBJQwIz8fBTIXM5trUHo7jvRbUp0iOPPPLFqVuukdh46sajjz7a+VdrcwxTbm7uww8/HBcXd+zYMWh1QOCt7Nq1y2q1btmypXX0TUxM/Prrr5uamqqqql5//fWvv/4aAGJiYtavX9/Q0NDY2Nh+Ok+HP71VYiwWS1fT6XwYCB2o2gQ5ejo9nK0PWaFQbN++Xd9rZOx/Ds/fcjThP4d/bA5ITU319va+8x8jSdtbTKP9SKg7u4FwVm/yQ4GY5o6uXLnyoioo5Wxx+x+9f/Qa6TeySzNH2hzD9Lvf/e4Pf/jDE088UVtbO3bs2NZH9XbIarU++OCDeXl5K1eubLn42GOP8Ty/ZMmSZcuWabVa2ys/9dRTlNLFixcvX768/Ry6Dn/aYWJspzh1NZ3Oh8cwOdAXl/ntBXTLXYrO/4kTjmFqUVxcfOXKlf9XEjp9ZNSTg7uQSGfCY5jac9wxTI8ftg7zJ8+xfc7DwG8tKZMVIwN+KysYP4aptLR09uzZ/oaChUPDogO9rBRO62o2nip0GzQmJSXF19fXXqm9valTp7YZ6hMXPIZJrLbm03mR7Ib2sLCwsLAw/XX+P5f4JwcLnRokNAsPPxbyr41mvUy4rxfZfoO2DoSMCwkJOXLkyKZNm/77zTfn9pxTKBTDhw9/5t2XH3jgAYlV/cWL9UwvXvUWOFjKfz6xC1PCBHF3OLcszVprBm/WU4oc61AZ7e1JenuyXjTf35t7Psv62iimm61tKJXKpUuXLl26VMA0/PnPfxbw7owTU2YSl52FfJyG+DKzv+iteKogUUt2FYppAgJyhO0FbC14vZXxWlJQR4vqmRjTEZFJkyYJnQR2iSDfi9T3+XSuGIoVAJjVm9tWgMWK3P1USGdGsN4cBAAFgXsiuO2YY5H9iKOkFh0LD7uL+Pt6iePjvb83t6uIN2GbUMYuVFOTFYb7iyAQAsD9vcj2G5hfkd2Io6QWnQwdjfQiIbfcWYktWjcY7EsOlmIVW75+LKT39hJHFASA6eFclo7WmoVOB5IKDIQOsbOQnyGGXqYWs3tzPxRgFVu+fr7B3xMhmtLAUwXjtWR3EeZYZB+iyfri8nMhFVGxAgD39SI/3cAWoUzVmOB4BVs74t7R/b1xmBDZjZgKa7EoqqclDTQmSEzFyiBfQgicq8aSRY52F/GJwcRdVGup7utFdhbyFmwTyonZbF6zZs3s2bOXLFmycePG1jum9pCo8r5I/FxIp4dzCjHFQQCAGeFkVyGN9hVbulGP/VxIZ4qqAwMAQt1JpBfJ1NGkYMfm2MOHD+OydxacOXPmaE5OYNGZ+3oF1Okt7//pxw0bNuzcudPV1bXnL46B0P5+LqTzo8T35EwPJx+c458fJrICEfUQT2FXEb92jPiKgunhZE8xnxTswN0BX3/99TVr1pSUlDjuFg7S1NTk4uLCcdJ5nPfu3ZusUb856ZdNsOYMDHl467H333+/hxvg2eBeo3Zm4kGzyXx1gSpA3Z0/d+Zeo23UWyD0K3PxwypPlraYwb1G27PvXqNZOro83XpqrvgCYVoZfT7LenS20nF7jYpXRUWFj49Plw47ZFl5eXlkeNi5p6YoWx3tuPNK+dcNfvv37+/560unvsCIQ6V0iB/pXhQUlocSxgXhIgrZ+bmQF8U6+vbiNeSqkeoahU4Hcjyj0ejrqlLefMBxoLtLbW2tXV4fA6Gd/VwopmnobUwP53bhlHSZ2XGD3iuSnR/aUHIwKYRLLcYcK30RERFGUBbV3lTrOVpSHR0dbZfXF+UDwLKdhfQecdavAWBGBNlZiC1CGSltgBt1NFZUM5xbuzuc7CrCHCt9arX62WefXbXrVGldk+3K4RuGj/IKn3vuObu8vvgGBlh2zUhrzWI6IKaNoX7EQuFSDR3gI9a3gLoktZifHMopRVsfvjuCvHLcymMolIG1a9e6urpO/Nvbfd25WpPZVRvx3Q87Ro8ebZcXx0BoTz/doDPCOVHHkOnhZGchBkK5SC2mU8NE/F1HeJAANTlZSfsxf8wL6iGFQvHKK698EfU/r4VfHhvh06dPHztOrhRtVZBJqcV0WriIixUAuDuc7MZBF3mgAPtKeFEHQgCYHk5SS8X9FlAn5RtpHXWZN2FkZGSkfZcYYCC0GwsPaWX85FBxf6R3hXFHymiDReh0IMc7XUk9lCTSS9xRZHo4lyq+ZX6oO3YX0+mO6XITd6nNlGw9jfQiQXbY5UBI3ioYFUgOl+Goi/SJvV/UZmIwyTOA0SL6N4LuyHE5FgOh3ewtpneJv1gBgKlhOCVdFlKLRd8vCgBuSogNoml6B+4vg1hg4eFACX9XmENiFgZCu9lbwk91zJfkZFNCyb5ibBFKXJMVMsppcogUcuy0MLK/HOf9SVyOnvb2JFr77KfUlhQeAxYYzXDSQBO1oq9fA8DYQFJQR/VNQqcDOdKRcjrUj/iJcAuk9u4KoQcwEEqdQ3vyMRDax8FSPiaIuEniYVRyMCGEO1CCvaNSllrMi32Gc4vBvtDEQ74RuzGkbE8xPz3cUQELA6F97C2mDuq8FsSUULKvBIsVKdtTRKXRkw8ABCAx0Loft8mVrmoTnKmk4x3W5SaRJ0Fw+0okMlPGZkoo2YvDhNKlb4LrRpGdHX17SUEWHNiWsAMl/PhgonbYjCgMhHZQ0kDLGuho0e6s1l60H2m00uvY1yRRe4v5iSGcSkJP/0SNZV8JbrUmWftLqEOXaEvoURDO3mI6KZQT99ZqNyMAk0M57B2VKsks9WkR4U7dleR8NeZYaTpQSieFODDHYiC0A+kVK4DDhJK2v5TeFSq1HDs5lOzHHCtF5Y1Q0uDYwwwwENrB/lIp7NDRxpRQsh/7mqToRh1tttJBvlLLsRgIpepACT8xmFM4MsNiIOyp89XUhYMokW/Y2F4vT+LjQs5UYskiNftLJLKOvo3JodyhUt6KGVZyDpbSZEf2iwIGwp47UOLYzmsBTQ7B3lEJcvRwi1CC3SDYjeQZMMdKzf5SOtnBPfkYCHvqUBmdKMViBQCmhJF9uKxecg45vn4tlClh2DsqNSUNtKqZDvHDQMgwCnColJdqsTIphEsroxYMhRJytZZaKPSX6MHLk0Ox6iY1+0pocojD5+RjIOyR89XUU0V6eUqzWAl0hT6eJBf7miREqv2iNskhXGY5bbYKnQ5kPwdKHN4vCgB22xzTaDS+/fbbZ8+eHTp06AsvvODl5dXyo2eeeeb8+fO2f8+cOfO5556z100Fd7BEsr1MNhNDyKFSSW1BInMHJR0IfV1ggA/J0dOkYMm+R7k5UEpfGO7wBpvdbpCSkqLValNSUjQazebNm1uuU0qLiopSUlK2b9++ffv2lStX2uuOLJDwAKHNxBByqBT7mqTDCRPwhDU5lBzATUel4rrRSUt97BYI09PTZ82a5eLiMmvWrLS0tJbrBoPBarW+/PLL8+fPf+utt+rr6+11R8HZBggnSrruOSGYO1JOcUq6NFyqoRxAX28p59iJIVxaGVbdJGJ/CZ3klKU+dusaNRgMWq0WALRabWVlZcv1ysrK/v37P/nkkxqN5qOPPvrwww9ffvnlDl9h2bJlhLR9RMeMGfPwww/bK5H2daGWuCtcfay1NTV2e826ujoA4HlWnmQVQLCr65EC4wg/YZLE83xtba2rq6sgd2dTbW2tWq02mUxd/cOd15SJgVxNjXQqozYmk6murk6pVALAcHeSVe6qr6pxkff8h9raWgBQqVRCJ6RH9txQJwZZa2osPXkRT09PheIO23XbLRBSSm1hjFLauhwfMGDAO++8Y/v3smXLli1bdqtXGDZsWPuLvXv3Zva7zDRwSRpq3+TZXo2ptzwxGDIMqrEawQKhSqVi6gMRnEqlUiqV3fhMMgzKySG89D5MSmlLJvFXQX9vOGN0iQ2UdT+G6ldCJ6RH0vWKV4f3NMe2b1+1Z7dAGBAQoNPpwsPDKyoqAgMDW65funTJbDYPGTIEfv1ubvUKzz77bGdSzI7MSuu9vYi7uz0P+bZarQDg7u5ux9fsoSkR/KYr9CV3h52Acls8z7u5uTH1gQjOZDKp1Wo3N7eu/mGazvx2rNrdXUxPWWcolUqr1dqSSZJDrdlVykm9ZN0kbGhocHd3F3UgvFpLFcQ6VOuMZ99ueSUuLm737t2U0t27dyckJADAyZMnAaCpqWnNmjUFBQVms3nTpk3jx4+31x2FJYcBQpuJwVx6GY+7jordhWqq4kgfye0F2N6EEILDhBJwqIxOcNbELrsFwiVLlly7dm3hwoX5+fmLFi0CgNWrVzRPuuAAACAASURBVAPAsGHDFi1a9Oqrrz700ENGo/E2XaPiIu0VhK1p3EDjRk7hpqMid6jMGeuxWJAUzGXgDC/xSy9z3jIYu3WNenp6vvnmm62vpKamAgAhZPbs2bNnz7bXjRhxQOorCFubGEwOlzn2GBTkaIdLJXhYWIcC1BDuQfIMdEygLN6vVB0uo390/ApCG1l3o/fEoVKJryBszbasXuhUoB5JL5fRMnPMsWJX0kBrTc47LAwDYTcdLuMnyKlYOVyGo4QiVlBHzTztJ+kVhK1NCCaHyzDDitjhUpoY7OgdRn+DgbA7LtVQFUd6y2CA0CbUnfi6kHNVWLKI1eEymhQso4d9QgiXhjO8xCytnDqzpSGjZ8OO0sqc+iWxAPuaRC3NifMOWBDsBoGu5CxW3UTrcCkGQuall9NEORUrADAxhBzEQChazpyAx4iJIeQQ9o6Kk6EZCuvpCCfOzsNA2B0ybBEmanFtllhVNEFpAx3q4KNNWTMhmBzGqps4pZfxCVqicGKGxUDYZWWNUNlMBztrOhMjIr2IkiNXa7FkEZ/DZfx45xYrLEgOIYdwhpc4pTl9SBsDYZellfGJWoefmMygxGCSXo4Fi/iky2ymjE24B3FXkovVmGPF53AZTdI6tYSV3ePRc3Kbd9AiQUOOYCAUocNO3KqKKYlazLHiU2eG89V0rHMPA8dA2GWyDYRJwSQdZx+ITZ0ZLtXIdI+V8RgIRShDR8cGElfnbvKPgbBrakxwpZaOlmWxMsyflDRQfZPQ6UBdkaGjYwKJPA/nG6/FznzxSSvjnd/SkOXz0QOZOjpWrsWKgkC8hhwpx7mjYiJIscKIIX6koonqGoVOB+qKtDKa6PQhbVmW6D0g52IFAMZruQysYovK4VI6QX4zZWw4AvEako5VN/Ew83C8gsY6d4AQMBB2VboQtRV2JAaTNBwmFI9mK+QaaJxG1lU3HCYUkdwK2teL+Lg4+77yLdO7wVasxMu4WIkNIqcraYNF6HSgzjlWQQf6EE8Rn1LeU4nBOF9GTI4ItGkXBsIusBUrXjIuVtyUMMyfHNVjySIOGeV0vHPXY7FmXCA5g1U38TgiUI7FQNgF6WWy22K0vUQt9o6KRka5rDswAMBNCcP9SQ5W3UQio5xPwEDIuAwdTZR3/RoAxgfjxFHRyNTx42VfdRuvxfWv4nCllioEOt4OA2FnUYDMcl7m9WsASNRymTpqxYKFeZdrqFpBIjzknmOx6iYWR8oFa2lgIOysSzXUQ0XCZF+sBLpCiDs5XYmRkHVCDbewBqtuYiFgjsVA2FmZsh9uaYFbOIpCpg5zLMCvVbczeEgv8zLLqSADhICBsPOwft0iQUsydVissA7ndrXAYUL2VTVDQR0d4Y+BkG1Yv26RoCW4vwzjqpqhsJ4Ok9lhvLeCu2+zL1NHY4KIUqCIhIGwU6qaobCODheotsKaAT6k1kRLG4ROB7q1TB0dJ1yxwprxWHVj3pFyXsAuN3xQOiVDR2M0WKz8ggDEa0kGzsRj2JFyHpf6tOjvQxostLgeYyG7hN29Eov2Tsks5xOwX7SVeA2Hw4QsO1JOE7T4dP+CAMRpuCzMsawy85BroDFO32u7BT4qnZJRTuOxWGkFhwlZZuYht0LWe223F68lGRgIWXXCIMxe2y2wcL8zCw/HsFi5WUwQOVVJm6xCpwN1JM9A+wharDAoXkMyserGqgzhFk7YKAW8N/uuXbv27bff5l4p8rRGmmYsBo1G6BSxwl0Jg33J8QpcUsIiAXfoYFZL1c1VIXRSUDsZ5fT+3kLmWGwR3tLHH38cO2rE1ZSPhl3ZN+noe0MGD9qxY4fQiWJIPK4mZFWmjsZjILyZuxIG+ZLcCsyxLMrSCdzlhi3Cjl28ePGl55/74cFxff08bFfuL6x87NFHr16/7uPjI2zaGJGgIZuvUxgmdDpQOxnl9K1xWMdty7YRhLBdcKi9kgbazNN+3tgiZM+2bdtm9g1qiYIAMD7Cf6AH2bt3r4CpYkoCrqBgUkkDNfE0yguL+7biNdiHwaIj5TReI3AkwkDYsYqKinAvtzYXw7zdKioqBEkPg3p5EiVHrhmxZGFLRjmNE7pYYRPOl2GT4P2igIHwVqKios7oa9tcPKurjYqKEiQ9bMINOxiUUU5xzWuH+ngRQiAfq26MEXzKKGAgvJX58+fnVFm2Xyq1/ZcCvH/0Gh8YlpycLGi62IJ9TQzCmTK3EYc5ljHNVjhdSccEYiBkUkBAwE8///xBgeWerzNX/nxy8sb0dKLZsWOHSqUSOmkMiddgi5AtjBQrzMKqG2uOV9BBvsRD6FmbQt+fYTExMTm5eUFv5rwQXTJkYL+RI0cSguXLTUYFkiu11GgGL6wesOGEgQ5koFhhVryGfH0VZ3gxJIuNU33wibmd0zXKwaNjF8/GT6ljLhyMDCBH9XRyqPBZGYFtL0AGihVmjQkkF6ppnRk8serGhgwdfaCP8DkWu0ZvJ5OB6UyMiwsiuJcxO/DUzNtTK2C4PzmGy+qZwUiLEAPh7WRjILyTOA0GQoZgjr2jeC3mWFbcqKNWnvZhYM0rBsLbwY6mO4rXkkwdj+UKC2w7dPQVdIcO9mHVjR3snOrDRCLYVNJAG61YrNxBqDtxU5BrtViyCI+FHTrYF6chmTqcL8MEdnry8bG5pcxyGq/hmPiW2IZrsxiRiR0YnRDhQVw4ch2X1TMAA6EIZOloLBtfEuOwr4kRmGM7KRZzLAOarHCuipU1rxgIbwk3qu8kDIQsaLbCqUo6Lghz7J1hjmXB8Qoa7Ufc2FibhoGwY2YeTlXScWzUVhg3JpBcrKENFqHTIW+5BiZ26BAFXPPDAhb22m6BgbBjeQYa5UVw1W1nqBUwxI8cx7VZgsrGftFOGxNIzlbRJqvQ6ZC3HD2NYaYDAwNhxzLKcefiLsC+JsExsjBZFNyUMBhPqxcaU3O7MBB2LFvPULOdfdjXJLhslurX7MOqm7BKGmiTlUYxszgNA2HH2JnXKwq4NktYpQ1gNNH+PphjOwsDobAyy2kcS4vTMBB2QNcItVisdEUfL0KA3KjDkkUY2Xo+ToNno3QBBkJhZevZGtLGQNiBLB0fi8VKF8UE4bJ6wWTraCzuKdMVfb2JycrnXSu1WnHOjACYmjIKGAg7hMMt3YB7GQsIp4x2SX19/QsvvGBYHjB1WB8fb+/ly5cbDAahEyUjFh7yDGytecVlRx3I1tE/DFMInQqRiQ0ifzqKw4QCsFI4XoFVty6YN28euXIie9FYjYe6ptm8Ln3n1KnHsrOzVSpcL+UMpyppb0/izdKHjS3CtngKxypoDNavu2hcEDlVSZuxn8npzlTRcA/i6yJ0OkQiKyvrTMahf909XOOhBgAfteqtydG09NoPP/wgdNLkIlPH3OI0DIRtnaumWjcSoBY6HWLjroT+PuRkJfaOOlu2DuttXZCXlxcX5qdW/Fb0cYQkRQTk5eUJmCpZYWopvQ0GwrZwYXK3xQaRbBwmdDrW5h0wTq1W15vbdlzUmSxqNVZ+nYTBxWkYCNvCeQfdFqsh2XoMhM6GgbBLkpOTs0qqS4xNLVeMzZY913RTpkwRMFXyYWgGXSMd7MtWjsVA2FaWjsYy1mwXC1yb5Xw1Jiiqp0MYK1ZYFhkZ+cwfX1rw/dFtF0svGep2X9XN33L07vmLEhIShE6aLGTpaEwQYWgtPQBgIGzDaIb8OjrMn7FvSSQG+hBDE61ouvNvInvJ0dPRgUSJz3FXrFmz5pPN2/a79V1xvObZkpDXPvj0008/FTpRcpGj41kbIARcPtHGUT0dGUBUWKx0C0dgXBDJ1tOZEcxldKnCftHumTp16tSpU3kK/l+ap9zP0kR+qcvW01XRzC1OwyL/JgyO4opLnIZk46ajTpSl4zEQdhtHYEwgycGBbWehTE4ZBQyEbbD5JYlIrIbDiaNOYytWYoPwKe6+OA1OdXaeSzXUz4Vo3IRORzv4CN0kG+vXPROnITl6ymPB4hRXaqiHkoS4C50OMYsNIlnYh+EszPbkYyD8zXUjVRAS7sHi9yQWAWoIcCUXazASOgNrW/iLUZyGy9FTzK/OweziNAyEv8lmtbYiLris3mmwJ7/nNG7g7UIuY9XNKZg9zwAD4W+y8FR6e8Bl9U7DbEeTuMQGYY51hkYLXKqhowJYzLEYCH/DbLNdXHBZvXM0W+F8NR3NZLEiLrE4X8YpjlfQIX5EzdzSCQAMhC1MPJyupGMCsVjpqZEB5HINrTMLnQ6pyzXQgT7EDVcC9xhOHHWOLD27m3ZhIPzFSQPt5008sFjpMRcOhvmTXAOWLI6Vo2N0uEV0RgWQ89W00SJ0OqSO5RyLgfAXONxiR9jX5AQ4pG0vrgqI9iMnsOrmYCyXsRgIf4Ez0e0IZx84AU5ytqPYIBzYdqzSBmiw0ChvRnMsBsJf4Ex0O8IWoaPpm0i1ifb3wRxrHzjV2dGy9XyshjCbXzEQAgBUMnlElnhFeREzT4vrsWRxlKMGiAlit1gRHVz86mhZOhqnYTfcsJsyZ8rW0XHsHZElauOwd9SRjlYQ7Be1o34+pN5CSxuETod0Mb44DQMhgK3Zjv2idhUbhLtvO4pOp8ssNOJe23ZEAMYFkaN63HTUIXgKJwx0HMOL0+z5LBmNxldeeWXOnDmvvvqq0Wi843V24ACh3eGgiyN88sknISEhiaOG5y4OeW1eQk5OjtApko7YIA5zrIOcraLBbsRPLXQ6bs2egTAlJUWr1aakpGg0ms2bN9/xOiMoQI6OxjLcfy1G44JIbgW1YsFiP+++++47f/rDxilRuU8kn3tqyiN+DdOmTD5z5ozQ6ZIInOHlOOzPybdn6Z+enj5r1iwXF5dZs2alpaXd8TojLtdQLxeiZe+ILFHzdYEwD3K2CksW+7BYLG+99dYHM4YPDvQCAAIwZ1DIsqEh77zzjtBJk4iYIHJUj1U3h2B5BaGNPXdSMRgMWq0WALRabWVl5R2vtzF37tz2F+Pj45944gk7JrK9AzeUo3yUVVVM9NnW1dUBgMUihV0uRvu6HshviiA92myN5/mamhoXFxd7pUqkCgsLSaPRFgVbJPYKeP3EiaqqKqFSxQiTyVRXV8dxParWE4AgtXvOjZpB3lIYKayuruZ5XqVSCZ0QAIDMMvdFYU1VVcJ8sN7e3grFHXY4tWcgpJQSQmz/4Hn+jtfbmDt3Lmk3Ibx3795ubo5trOXVcgnB4Oi7dJItBDKSmB6K05DcSsWKnu2GyfO8q6urND6QnvDz86s3WS08Vbaa3FzbbHZ3d8cPR6FQWK3Wnn8OsUHkVJ3rKK0UWoVubm5ubm4sBEKjGW40cGOD1SqBRp/ah5X27BkIAwICdDpdeHh4RUVFYGDgHa+3sXjx4s6k2O6OGSyPDFS4ujLRcjebzQDg6uoqdELsIDGMfnLZ6urao8Ycz/NqtVoaH0hPREREDBgydPul0rmDQlsubj5XPP3hJ/HD4TjObDb3/HNICOFzK+kKVybPR+gi21PDQiDMrKIjA6xe7kznUnvG6Li4uN27d1NKd+/enZCQAAAnT57s8Do7mqxwvpqO9GciCkrMMD+SX0eNeAyFnXz44YdvZBW8m3XleGl12g3Dip/yityDn3/+eaHTJR04X8YRshnea7uFPQPhkiVLrl27tnDhwvz8/EWLFgHA6tWrO7zOjtwKOtgXz7JxCCUHIwPIUZySbicJCQl5Z84ejZrx5DnuP0bfyU/+MTs729vbW+h0SccIf3LViCeI2Rn7M2UAgFDKRDlFCOF53vldo++d4a/V0n8lsNIZYltn6eXldcffFIXV2dYAV/KnEd2vb/E8r9frbZOtEADM22e9R9O4sJ8ChwZbmEwmo9EYEBDQ85dK2G55a5xiYgjrBfcdVVRU+Pj4sNA1GvZfS+b9il6eTH+kcl88x/jGP2KHfU12l6WjY/ytQqdCsnAjCPu6UUd5ShmPgoCBEPeUcajYIJKtk8JkdEaUNFCTlfb2wJLaUXD3bfvK0YtjrxIRJNFxdI1Qi2fZOFIvT8IRUlCHJYt9ML6FvwTEavBgQnsSxQAhyDwQZur4GIaPyJKGGKxi2w/25DtapBexUjxBzG7Y31zNRtaBEPtFnQAHXewIc6wTxODu23Zi4SHPQMcwfOhEC1kHQuxocoI4nC9jJ1YKxyvoOAyEDoYzvOzlVCXt40m8hZ+4emfyDQM8FitOMS6InKykJpwx02Nnqmi4B/GV+66rDhcbhMOE9iGWflGQcyA8V021biSA4SOypMFDCVFe5HQlliw9lS2SeQdiF6shJwzUglW3HsvWUbEceC7fQIjDLU6DM/HsAnOsc3ipoJcnOYMniPVYjp7GiKTqJt9AiBPwnAbXZtkF5linwd7Rnqs2QXE9Heonjhwr30CYpaPxWKw4BU4c7TmjGQrq6DCRFCtihzm257J1dEwgUYgkw8o0ENaZ4ZqRDsNDJ5wi2peUN9LKZqHTIWZH9XRkAFHK9Hl1Npw42nMi6hcF2QbCYxV0hD9xkem7dzaOwNhAkoNV7B7IxJkyTjTUjxTV0yqsuvVAlo4XUZebTENBFg63OFeshmThpqM9IKIJeBKgIDAmkBytwKpbN1GAHL2YFqfJNBDiBDwniw0iaZd1169fZ+TYL9HJ0fNYdXMm7B3tiau11ENJQt1Fk2NlGghxSZYzHTp06KWZI/Me6zNx5BBNUNAHH3yA4bBLrhmpiiPhHphjnQcPTukJ0fXky/Fo9oI6SoH2Zv6ILGk4duzY7Jkz3p08aNpdkwDgcmXdk6+/1NjYuHr1aqGTJhpi2cJfSuI03BNpVgqAn3s3iG6pjxxbhLjFqDO9++67T4+KmBalsf23v7/nv+4e/ve//91qxdNlO0t0xYoEhLiDh4pcq8Wui+4QXZebHONBjh7nHTjP6dOnE8IDWl+JDvQyG6vLy8uFSpLoZGOOFQIuq++eJiucr6ajAsSUY+UYCDPLRVZbETUPD4/qZnPrK2YrX2+2uru7C5UkcWm2wtkqcZxlIzG4rL57jlfQIX7EVSF0OrpCdoHQxMOpSixWnGf69OlfnS5sfWXz+ZLRY8f5+voKlSRxyTXQgT7EXY6j+QLDE8S6R4yL02T3eJ000P4+xFMMR2RJw+rVq5N//PF3O048PDTcXaU4kF+xJb9m5+49QqdLNHJEWKxIw+gAcraKNllBXI0bweXo6X29RJZjZdcizMQtRp3L29s7Kytr5jN//p6GrS7wPD/sgdNnz40ZM0bodImG6GaiS4abEgb7klxcVt9FmeXiK2NlFwhxAp7zubi4PPvss9u3b//7N3sDF/9Nq9UKnSIxydbTOJwpI5A4PEGsi0oboMFCo7xFlmNlFwhxSZaA4jQksxyLlS7QNYLRRPv5YI4VBgbCrsrS8fFaIrr8Kq9AqGuEKhMdgMWKQPp4EZ7SwnosWTorU8fHaMRXrEgGBsKuytbTmCDxhRXxpbgnsnR8bBAWK0KK03DYKOw83PxBWH29iYmnRVh16zQxDhCC3AJhth4HCAWGa7O6RHQ7dEhPTBCHjcJOsvCQZxDToRMtZBYIdTRWhM12KYnHvqZOs1I4XkHH4ZpXQWHvaOedqqS9PImPi9Dp6DoZRQWewrEKbBEKbGwQOWmgzbjPaCecrqRhHsRPLXQ65A0DYedl6Gi8VpQFrIwC4ZkqGuJO/LFYEZSHEgb4kJOVWLLcGa55ZUGshpyspCY8kakTxDsnX0aBEM/4ZkQsVrE7BwcIWeChhL5e5KQBc+ydiTfHyigQirfZLjHY19RJGTqagDmWAZhjO0PXCBVNdJA4F6fJKBCKdF6v9GCx0hmGZtA10sG+mGOFhzm2M7J0fJyGcOLMsHIJhJXNUNpAh/iJ81uSlgE+xGimZY1Cp4NtmeU0NkisxYrEYCDsjEwxr3kVa7q7KktHxwURBRYrDCAA44JItg6nH9yOrX4tdCoQAMBAX1JtouVYdbst8Q4QgnwCYaaOx+EWdsRruEysYt9WRjlN0Mrl8WQcAYgJIllYdbs1K4VjFTRGtLMR5fKkZZTTeNE226UnXkMycKO1W7MtpRdvsSI9cRrcX+Z2TlXSCDGveZXFwby2YgWX0rMjVkNOGKiJBxesnHTkTBUupWdLvJasy8MWYQf27NmzY8eO7HyDW0i0YdJTAQEBQqeoO2RRDp2poqG4lJ4lXiro60XycG3WLWTgDGfGxGnI8QpqxlDYCs/zjz766FMPzQ08vXeJ+ezgrE+iBw08cuSI0OnqDlm0CHHhBIMStCSjHHv/OpZZTieE4CfDEG8VRHqRk5V0LG79+qtvv/0288fvdz0c76ZUAMCcQSGxF0uXLl166dIlIrYzfmTRIszEpfTsidcSnC9zK7i5GoMScGD7Ztu3b186IsIWBW1mDwypKys6e/asgKnqHlkEQvHugCdh8Xha/S1UNIG+CZfSMwerbm1UV1drPNoOOGk8XKqrqwVJT09IPxBWNEF5I43GYoUx/byJGY887UiWDpfSswhbhG0MHDgwt7Sm9ZU6k+VqVWP//v2FSlK3ST8QYrHCrFgNhyVLe7iUnk39fEijhRZj1e1XK1as+O9F3f58ve2/xmbL/+w5M/vB+VqtVtiEdYP0A2EmFiusSsC+po5klNN4XErPHgIQr+UyMMf+auDAgdt2/Pjmhfopm47M+y5n/Mb0PlPnrF+/Xuh0dYf0Z41mlNOXRiju/HvI6eI1ZHU2Tki/iYWH4xU4U4ZRtoHtByOFTgczkpOTPz9w+rFvL/x9ZFV0dHRgYKDQKeomiQdCW7GCS+nZNDaQnK2ijRZwk3g27IK8Strbi/i4CJ0O1JHxWvLHHKy63STHoJg8JnrCeHE3NiTeA3Oykvb2JL5YrDDJTQlD/MixCuxr+k16GR2PS31YNTaQnKmkDRah08GSI+VSODVT4oHwSDkWK0zDYcI2MjDHMsxWdTuOVbdWMnVSyLEyCITBov+SJAxXE7aRiafSsw2rbq3lG6mZp5Feos+xEg+EGeU0AQcIGRavJZl4us2vCuqohadR4i9WJAwPTmktU0fHS2KGsxTew60U1FEzT/t6Y7HCrggPolaQyzVYsgD80oEh5UdSAhKDser2G8mMPUn5qUsvo4lYrDBvvJakYxUbAGzzDrADg22h7kTNkSu1mGMBfjk+Wgo5VspxIlOHxYoIjNeSIxgIAQAgUyrFirTZDk4ROhXCqzPD5Vo6KkAKOVbKgTCtjCbiTBnmJQZjIAQAMJrhci0djaf8MA+rbjbZejoqgKjFvYDwF5INhLVmuGakIyVRW5G2YX6krIHqGoVOh9CydHR0AHGR7BMpHYnBJL0MA6FEVhDaSPaxyyynYwKxWBEBjkCcBicgwJFyXhrzDiRvuD8paaD6JqHTIbTMcl4yY0+SDRQZEvqSJG+8lsO+psxyPD5aHBQEYjUko1zWVTee2s55lUgEkcjbaC+9HKeMisZ42fc1WSlk62m8VIoVycOq26lKqnUjGjeh02En0nzwLDwc0+Op9KIRG0ROVdJGGW/heKaKhriTQFeh04E6J0n2Vbf0cpokoamI0gyEpypphCfxVwudDtQ57koY4keOyngLxyNlNBH7RcUjJoicrpJ11S1dWnPypRkI08uxWBEZmc/EO1xGJ4RgjhUNW9UtRy/fHJtWhi1C5klmvwP5GK+V9ewDrLqJTpKWpMm16na1lgKAlDbFlWYgTCujEyRUW5GDRC2XoaO8LAsWW7EigS38ZWV8MMmQ65ofiTUHQZKB8Eot5Qj0wWJFVDRuEOhKzlXLMRJivU2MErVcRjm1yjHDSm2mDEgyEKaV0Yk43CJCiVqZDhNKr34tB4GuEOJOTldijpUCaQZCiX1JMiHbTUelV7+WiaRgOR6comsEfRMd6iepHCvBQHgYO5rEKUlLDpXKrlgpa4SKJhrtizlWfBI0cqy6pZXxCRrCSSvDSi0QljTQWhMdhMWKCPX3IVYK143yKlnSy/jEYE5ixYpMJAXLceKoJDftktr7OVxKE4OxVBGrCSHksMxKFuzJF6++3gQArsmu6ibBpT5SC4Rp5dgvKmITgslhmfWOYk++qMktx9aZ4UINHRsktRwrtUB4uBSLFRGbECyvFmGtGa7W4qmZIjYxhBySU47N0NHRAcRVEofxtiapQGhohsJ6OgKLFdGK9iM1JlpcL5eS5UgZjQnCUzNFbGKwvGZ4pZfxUtpitIWkHsH0Mj5BSxQS/JrkggAkBXPyaRQeLuOTJDfvQFYG+ZImKy2ok0uOPVRKk0MkmGMl9ZbSyigWK2Inq97RNNxrW/wStXKpujVZ4YRBmts4K+31Qkaj8e233z579uzQoUNfeOEFLy+vlh8988wz58+ft/175syZzz33nL1u2sbhMvqPWAyE4jYxhKy/IIstHBstcKqSxkpu3oHcJIeQQ6V0ST+h0+F4mTo6zJ942C1oMMRuYSMlJUWr1aakpGg0ms2bN7dcp5QWFRWlpKRs3759+/btK1eutNcd26gzw/lqCU5nkpvh/qSskeoahU6H4x0ppyMCiLsUixVZmRAil2HCgyX8JIl2YNgtEKanp8+aNcvFxWXWrFlpaWkt1w0Gg9Vqffnll+fPn//WW2/V19fb645tZOjomEAJTmeSG45AopZLK5N+o/BgKZ8s0WJFVobIZobXwVI6QYoDhGDHrlGDwaDVagFAq9VWVla2XK+srOzfv/+TTz6p0Wg++uijDz/88OWXX+7wFZKTk9tfnDhxYie7UndeU8f60spKU3dSz4y6ujoAMJvNQidESGO8Xfbkk0k+zQDA83x1dbVKpRI6UfaXWuj+52hTZWWXjzmvqalRq9WNjTJoNXeOyWSqq6sjRLBaRay/289XGh+IYOjE+qqqKqvVascHp5knuRUe0S7GSrHtM+7j46NQ3KGF1KNA+NhjjxUVFQFAamoqpdSWESmlPP9bdX7AnIPIPAAAGk1JREFUgAHvvPOO7d/Lli1btmzZrV7txRdfbJ+Vg4ODPTw8OpOYIwaybjR4eIi7xLR9dJ18y1I1JQJWZXMeHkoA4Hm+sbFReh9IgwXO13ITw9VuSnVX/9ZsNqvVajc3N0ckTIxUKhWlVMBMMimM5FQrHhnEUIRoamry8PCwYyA8WkaG+UGQt7u9XtBpOO7OrdgeBcINGza0/DsgIECn04WHh1dUVAQGBrZcv3TpktlsHjJkCACoVKrbfDEzZszodp3OaIYLNebEMJVa5F2jJpMJANTqLheOUhIXAvl15npQ+6uB53kXFxfpfSAH9XR0oNXXozvvS/0ru6dKpAghJpNJwA9kSgT992WrWs3QeK/tqbFjIMyosE4OA7XYS9hbsFuHb1xc3O7duymlu3fvTkhIAICTJ08CQFNT05o1awoKCsxm86ZNm8aPH2+vO7Z2uIzGBOEAoUQoOYjTEGkPE+IAoZQM8yP6Rlom6b5qCQ8Qgh0D4ZIlS65du7Zw4cL8/PxFixYBwOrVqwFg2LBhixYtevXVVx966CGj0XibrtGeOFDCTwqV7JckQ5NCuAOSnol3oJROkm6xIjccgfHB3OFSyVbdbCsIx0txBaGN3drynp6eb775ZusrqampAEAImT179uzZs+11ow4dKKX/isdiRTomh5LfHZZssWI0w7kqGquRbLEiQ8kh5FAZnR8ldDocQ8IrCG2kEDyqTXC1lo7DFYQSMjqQlDRItq8pvYyOw558aUkOIftLJNuHIeEVhDZSCIQHSvh4DVFJ4a2gXygIJAVzB0uk2Sg8UMpLcsNGORsZQAxNtEiiqwkPltKJks6xUnhvB0opDhBKz+RQsl+iw4QHSqi069cyRACSQzhJNgolvMVoCynEj4OlWKxI0JRQsq9YgsVKjQku1tAYHCCUnClhZJ8UA+GRcokPEIIEAqGuEQrr6OhALFakJtqPNFhovlFqJUtaGY3T4BmEEjQllOyVYtVtbzF/V6jEC1jRP44HS/mkYA7PIJQeAjAplDtQJnQ67A0HCKWqnzdRcXChWmqxcG8xvStM4jlW9G/vQCmdLPXaimxNDiUHJNfXlFpMp2COlagpoVKbO2pohsu1NE7qPfniD4QldBIWKxI1JZTsL6VSKlfKGqG4Hg8Lk6wpYWSvtALh/hI+USv9Ofnifn9F9bSymQ7zw2JFmiK9iKuCXKmTzoK7vcV8cgj25EvW5FDuUClvlVAolEO/KIg9EO4uolPDOA6LFemaFALpFeI+UaS1PUV0WjjmV8kKdoMQd3LCIJ1IuLeYTg2Tfo4VdyBMLcZiReImhxLJBEIKsLeEl0OxImdSWvaTb6QNFhotgy43EQdCnsI+LFakblIIOVKhkEZf09kq6qYgUV6YY6XsrjCyVyo7Iu0ppneFyaLHTcSB8HgFDXEnoe5y+JrkK9gNwtzoUb0UIiH2i8pBcgiXraNNVqHTYQ+p8ugXBVEHwj2y+ZJkbpLGtKtIClXsvSXSX5iMvFQw1I9klIu+6sZTOFDCy2Spj5gDYRE/PVzE6UedNElj2V0k+mLFxMORMtwUVxamhZM9xaKvuuUaqNaNhHlgIGSY0Qx5Bpok6X1gkU2Mv+V8NTU0C52Onkkvo9F+xF8tdDqQ490dzu0sFH3VbW8xnSqbnnyxBsIDJXyshrhJeh9YZKMidGIIt1fkVey9xfxd2JMvD+OCSEkDLRb5kUx7i+XSLwp2PKHeOXie37FjR25ubmqFW0zSZIAYoVOEnGF6GNldRBeI+fjv1GL6bpx0dgZAt6EgMDWM21VEHx8o1kBSb4GjeiqfTXHF9D7Lysri4uL+8tTSuj3/HZn9ybfLk1esWGG1SmJ6FrqtGRFkVxEv3gq2bcPGeKlv2Iha3B1Odol5YHtvMR+jIV4SWcF7Z2JqEa5cuXKIqfyN+bG24uTZ2L5zN3/92Zgxy5cvFzhlyMEivYinipyupMP9RRlLUov4CcGc5DdsRC2mh3PPZZktvEIpzi99ZyGdIaepiKJ5q/X19T/9+OOLCf1bCkIPleLZ2L7ffPONkMlCzmLrHRU6Fd30YyG9t5coQzjqHq0bRHqRLJ1Yc+yuInpPhIxyrGgCocFg8FYRT5ebmrC9fNzKy8uFShJypunh3G5xria0UthdxMuqWEEAMCOc7BRnjj1dSRUEBvnKKMeKJhBqtdoGqqhoMLW+eKGirk+fPgKlCDnVpFByVE/rzEKno+uydDTcg4TLYz0WanF3OCfSYcKfC+XVHAQRBUK1Wr14yZI/7T/XZPlldkyxsemdzMvLli0TNmHIOTyUMDaIHCoTX8ny0w1+psyKFQQA8VpSYKSlDUKno+t+LuTviRBNaLALMU2Wee+99554oi75y62xoX6NFutRfeOf1rw5Z84codOFnOTucG5nIT8zQmSLEH4upB+OF1maUc8pCEwO5fYU84/2F1NQqTZBnoEmh8ir6iamQOju7v7VV1+dO3fuxIkT7u7uGxIStFqt0IlCznNfL3LXTv5fCSCiZ7SwnpY00FhcOCFLd0eQXUX00f5Cp6Mr9hTxScGy26tEfG83Ojo6Ojpa6FQgAQzyJR5KyDPQUQGiiSs/3aB3h+OR9DJ1dzh5IdsqrkUUO4voDJn1i4KIxggRAoD7e5EfCsQ0E+/HG/xMXDghV6HuJNKLpIvnJAqewq5Cfob8hrQxECIxub83t71ANMVKowXSyiiekSJns3pzIqq6Ha+g/mo5nh2NjygSkwQtKaqnN+rEEQsPlNJRgcTXReh0IOHM7kO2iafq9nMhlWFzEDAQInFREJjZi9txQxwly483+Ht74SMma0P9iAsHeQZx5Nit+fzs3nLMsXJ8z0jU7u9Ftouhr4nKcmEyau++XuQHMTQKrxlpeSNNkOUhrxgIkchMD+ey9bTadOffFNYxPVUrIFpO+1ShDs0WyTDhlut0dh+Ok2WGxUCIRMZdCUnBhP19R7fk8/MiZVmooJslaElpA71uZL1RuDWfn9NHphFBfOsIEbq/F7e9gPVzerfm068n4YYyCDgCow2Hn/zjgWGeTaNHj54/f75SyVzBW9YIF2tkt6FMC5nGfyRq9/fmdhbxZobbhHkGauJhVKBMixXUwmw2L1iw4MKaWSNOpbjn/PjhS78fM2ZMcXGx0Olqa8t1/t5enItcA4Jc3zcSM60bDPQhB0vZ7Wv6Pp+f14dgGETvv//+1bQ9exclPBPTd/noPt/MHRejqFm1apXQ6WprWwE/u7d8MywGQiRK8yK5zdfYbRJ+d53Oi8SHC8F33333+3FRKsVvmeH52H4///RTY2OjgKlqw9AMOTpZ7/wg33eORG1BFNmaz5uYDIXnq6nRDDG40TYC0Ov14d5ura94qZXuCqipqREqSe39eIOfGs65Mzdw6TwYCJEohXuQaD+yh8mDT7+9TudFYr8oAgCIioo6o69tfaXI2Mi7uAUGBgqVpPa2XKdz+8g6w2IgRGK1IIpLYbJ39Pvr/Fy5TkNHbTz55JNvH7lcUPPL+bzGZssLe88+/vjj7EwcrTPD4TJ+pvxOnGiNlS8Doa6aF8m9etzcaFEwdXbalVpa3kjHy3J7DtTe7Nmzi//2zuxXX4n2Uni4KHNKax9c8ujf/vY3odP1mx8K+KRg4iPvHXFZKkIQ6gqtG4wJJD8X8g+wNC1ly3U6R67bc6AOPf300wsXLjx27Fi1sT6zcMSryyNdXBjKH5uu8EsHMPQECULu7x+J2oIoLuUaW8OEX13lF/bFxwrdxN/ff9q0afMfmDN3dG+mcmxZI+To6f2y3xpe7u8fidoDfbg9xXytWeh0/Cq3ghrNkBjMUH0fMWVhX+7rqwwNbH99lZ/dm2NqcEEQGAiRiPmpIVFLfrzBSsnyxWV+aX/sFkW3NDmUlDTApRpWGoVfXuYX98cogIEQidyCvqz0jlp42HyNX9wP4yC6JY7AvEjyDRs59nw11TXBBOzAwECIxG52by6tjNcxsE3Hz4V8fx/S1xuLFXQ7C6NY6R3deJlf3I8oMMNiIERi56WCOb25Ly4LWbIUFxefOnVqw9mGR7GXCd1JjIaYeTheIXCjkKfw36t0cT/MsQAYCJEEPD6Q+/QiL0i5cuzYsbFjx44e2PehKYk75wSeXP9iU1OTEAlBokEAHunPfXZR4Ebh4TIaoIahftgeBMBAiCQgQUtUHKSVOTsUFhQUTJ0yZXGA6diy5D2LErIeTTj9/edPPfWUk5OBROeJgdw313ijoLOdv7zMY3OwBX4QSAqWDeT+fcHZVeyPP/54VqTvvMGhtkp1oLvL+zOGf/vfr8rLy52cEiQuIe6QHCLkSGGtGbYV4ILX3+AHgaRgSX/uxxt8VbNTb3ru3LnYUL/WV3zUqgH+7hcuXHBqOpAIrRjErXd61a3Fxsv8XWFciLtQ92cOBkIkBQFquDuC+8q5VWwvLy9Do6nNxYoGk5eXlzOTgcRoWjipNcFRvQBD2xTgg3P876Ox8P8NfhZIIp4YyH3i3Cr2zJkzvz5TZLL+dtP9+XrwDhgxYoQzk4HEiAA8PlCYRuHeYqrmcP+jm2AgRBIxKZQ0WpxaxV6wYEF08vRZKdnfni/Zn6//e8blPxy88ulnnykUCqelAYnX4wO57/P56rZ9Cg73wTn+aWwO3gw/DiQRBGDZQO7D886rYnMc983mb3VLPt3vPzrFHOw2ad6JM2enTZvmtAQgUQtyhWlh3FdXnNooLKij6WX8wzhf9Gay32wVScjyQVz/zebCei7Cw0ndPjtu8P2S5+y470Hn3A5JzIrB3DMZ1pXRztuf9uPz/CP9OQ8s+G+G9QIkHX5qWDqAe++086rY/3eGf2YIPkSom5JDiILAz4VO6s9vssJ/LvErsV+0HfxEkKT8z1Dui8t8pVPWUZyspFdr4YE++BChbiIAL4/i3jxhdc7tvrnKjwkk/XA73HbwGUaSEuZB7vbWzXv29QULFjz99NOpqamOu9c/z/BPRXNKfIZQDzzQh6sywf4ShzcKrRTeOsmvHoYzuTqADzGSlH379u1cOjQq+4vkuovB5w48uWDO7373O0rtX8rkG+kPBfzyQfgEoR7hCPxpBPdmnsMbhZuu8KHuMDkUm4MdwDFTJB0Wi+WxpUvfSYqcFqWxXZkfHTbz683b7rtvzpw59r3Xmlx+1RAuQG3fV0VytLAvtyaXP1JOx2sdFaWsFNbl8Z8kYnOwY1ifRdKRl5enqKtsiYIA4KFSLBoavmPHDvve6GIN3VXEP4+9TMgeVBy8MJz720kHNgo3XOJ7e8LEEGwOdgwDIZKOmpqaIPe2bbQgD3VNTY19b/Tno/zq4QpvlX1fFcnXYwO4EwY4YXDISKGJhzfz+DWjsd52SxgIkXQMGDDggqGuwXxTzTq3tHrQoEF2vEuOnubo6dOD8dlBdqNWwB+Hca8ec0ijcMMlfpAPJDis31UC8GFG0hERETFj1pznU0/XmSy2Kz9dLt96veqJJ56w411ePmZ9dRTnhsPryK5WRnPXjbCtwM6rYJutsC6PXzsGm4O3g08zkpTPPvts9erV8V98PsDPXVffXOQV9fm3u/r06WOv199XQm/Uwe8GYA0S2ZmKgw/GK5Yesk4L49ztVzD/Nc8aE0RigrA5eDvEETPLu4EQwvM8IXL/toxGIwDgOT4teJ7X6/VarbZLf1VZWXnu3Lng4OBdTX223YC999inXGmwwOitlr/FcLN7CxkIq6ur1Wq1m5ubgGlgislkMhqNAQEBQifEDpYctEZ4wLpxPWrA6fX67Ozsvn37NgX2v3sPPTlXFYyZ5bawYoskyN/fPzExsV+/fk8NUeqbYGu+fbqbXjpqHRtEhI2CSNre+f/t3X1QE+kdB/An2QRUTsPbGaypPezAaKWMingI41S5gTuLmHJtbX1tI1jPqWXsYB2r2HqDJ3Jjpx2mp+MQpKjUajs4eD2nuRzqGDzQsXdApZEoimPOl1ySKxACSdh9+kfu0lQgvJhkY/b7+YvZXYbfPPNjv8mz++y+ylQbOP1/Jvn9pLe3t7CwMPkbirK3Nq9atiQrfdHPpDeQgmPC1CiEM0ZEfp/BbLnKfmeWOPb51vxdekTPd9P2N/EvAwEkn0r2L2J2fMw2TmoaY/PmzZz+evNPl78UISGEvG94snfrqq0Zt+bMmePvSsMKPttCmMv+muiHiaINl4e457gI0OsihTq2ajkTgxX0EGDb54t7nOSP/57wNMb9+/ev/OPi73JS3ClICMlPTsh/JVqtVvu7xnCDIITwV57ODLLk7U8mf2/6zmb29dmiNxRCv4YNQcCIyF9fYw61sh99NrHPbgaD4VsvT58q+b/ri4sTojs7O/1aYBhCEEL4k4jJX7IlNQb6waTed3P8Nqd7Qo+8ihvQIUgSp4vOZks2Xhky9EygY2Uyman/2RevmPodMpnMr9WFIQQhCIJ8Kjn7GlN4daird2JZeMLAHWrlNKuYl/AcGQii5QmiQ+nMmg/ZL8b9TrG0tDTHtOjG+597tticQ2c6jEqlMiAlhhEEIQjFspmisjRmxQfsjc/Hm4W1d7gD/+Q+WsXMnY5JUQi2LcnivDmiHzQO9TjHdbxUKs369Yntl++/ffX2+3eenGh9sOpM8xs/3pyXlxfgSl94fg5ClmVVKpX3lr6+vtLS0oKCgv3797sXyQHwZes88XuZ4tWaoT93jX0nQu0dbt9NTvtdJkmGFAR+vLuUSYkRLT4/dNM8xqc3+xD5ZQt7KyH70w7917+35cOIRFPKayfq/37s2LHglPpC8+eC+vr6+kuXLnV2dnq/DVWtVg8MDGzbtu348ePTpk0rLCwcuQ4sqCeEYEH9MJNbUO/bv6xUqWXXfVP0m0VM5EgX/jq+oL+6wd7tJQ05zPzokOtJLKh/RjgtqB9RfTe3/RpbupD5xYKRv7poP6NvNbEZM0V/WMa8PIUQQsxms0wmk0oxoT8u/vxGOHfu3I0bNz6zsampSalURkREKJVKnU7nxz8HMDnfjhVdV0raLFRxxvXzj9lmEyWEcJQ8GSCfWuj2a2z2xaHXZ4s7vi8JwRQEAXrzFXHzGsnJO9zShqF3WrlPzJQSwlHS1UvPd3ObrrDbmtj3spi6lV+mIEyUP1cHL1y4cPhGi8Xi/jgvl8utVquPX09NTR2+MScnZ+/evf6qMPTZbDZCiMMx7uvj4Y7jOKvVyjB+vmNTRMif0ojRLvqbMfInl6WmQZGdFcVIufhIunLm0LWVjugI2uOrW/nU09MTGRk5ZQrOeV9yOp02my1EnhYZIDMIuZBJrpkljU8lP7ottTqJgxPFRdD5M7i0WPadFY5pDDWb/3e81Wp1uVz4RkgIiY6OlkjGSLrnDUKVSmU0Ggkh3tOh3iil7glPSinH+bowU11dPXxqNDY2VlD3/orFYoKpUS8cxzmdzgD1gExGFswiv00npgESG0kkYvcECUNISC+bp5RiatSb0+kUi8VCOFEoY4gyiRBCjHYyQ0pmSEWEiAmRDO9Yl8uFqVG38XyMft4grKmp8X1AXFycyWRSKBRmszk+Pt7Hkenp6bhG6G5ctK8Hx3FSqTTQAzL7hRpv6Vf4LiRUUEqFNiCJY4U+mmRCArh8oq2tjRCSkZGh0WgopRqNJjMzM3B/DgAAYBICGIS7du0ihGzatOnevXvr1q3r7u7esGFD4P4cAADAJOB9hKHlwoULhJA1a9bwXUiosNlsZWVlFRUVfBcSQqqrq1NTU9PT0/kuJFTcvXu3oaGhpKSE70JCyOHDh9evX4+XTowTniwTWlpaWlpaWviuIoTY7faqqiq+qwgtFy9e1Ov1fFcRQoxG47lz5/iuIrTU1dU9ffqU7ypeGAhCAAAQNAQhAAAIGoIQAAAELVRulklISMjNzeW7Cv51dHQQQhYsWMB3IaHC4XBotdrVq1fzXUgIuX79+qxZs3AfhIfZbL5169aKFSv4LiSENDY2pqWlRUdH810I/w4ePDjmP0uoBGF7e3trayvfVQAAQFjJz8+PiYnxfUyoBCEAAAAvcI0QAAAEDUEIAACChiAEAABB8+f7CGESiouLPU8JycvL27lzp/vnvr6+ioqKjo6OlJSU3bt3C+rFTDqdrra21mw2JyYmlpSUKBQKz67RhitcjdYGaA+0hwfOIX5AgT8cxxUUFFgsFrvdbrfbHQ6HZ1dVVVVlZaXD4aisrFSr1TwWGWSPHj3Kz8/X6/WDg4Nnz54tLi727PIxXOFqtDZAe6A93HAO8QtMjfLJYrGwLLtv3761a9eWl5f39/d7djU1NSmVyoiICKVSqdPpeCwyyB4/fpydnT1v3rzIyMjc3NyHDx96dvkYrnA1WhugPdAebjiH+AWCkE9WqzUpKamkpKSuri4qKuro0aOeXRaLRS6XE0LkcrnVauWvxmBbvHixe26HZdna2lrvVdI+hitcjdYGaA+0hxvOIX6Ba4TBplKpjEYjIUSr1SYnJx85csS9vaioqKioyHMYpdT9UipKKcdxvJQaNN5j4t5y8+ZNtVq9ZMkSlUrlOczHcIWr0dpAUO0xHNrDA+cQv0AQBltNTY3nZ4PB4HK53A9Uk0qlUqnUsysuLs5kMikUCrPZHB8fz0OhQeQ9JpTSqqoqvV5fWlrqfR8E8Tlc4Wq0NhBUe3hDezwD5xC/wNQonwYHBw8cOPDgwQOXy3X69OmsrCxCSFtbGyEkIyNDo9FQSjUaTWZmJt+VBk97e3tzc3NZWVlcXNzAwMDAwAD5akxGHK7wNrwN0B5oD284h/gFHrHGJ0ppQ0NDfX19f3//0qVLd+zYERUVlZOTo9VqbTZbeXl5V1dXUlLSnj17oqKi+C42SE6ePHnq1CnvLVqt1j0mIw4XX3UGx/A2QHugPbzhHOIXCEIAABA0TI0CAICgIQgBAEDQEIQAACBoCEIAABA0BCEAAAgaghAAAAQNQQgAAIKGIAQAAEFDEAIAgKAhCAEAQNAQhAAAIGgIQgAAEDQEIQAACBqCEAAABO2/92awop8u1nwAAAAASUVORK5CYII=" }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Interpolations\n", "using Plots\n", "gr(fmt=:png);\n", "\n", "x = -7:7 # x points, coase grid\n", "y = sin.(x) # corresponding y points\n", "\n", "xf = -7:0.1:7 # fine grid\n", "plot(xf, sin.(xf), label = \"sin function\")\n", "scatter!(x, y, label = \"sampled data\", markersize = 4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To implement linear and cubic [spline](https://en.wikipedia.org/wiki/Spline_%28mathematics%29) interpolation" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "hide-output": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "li(0.3) = 0.25244129544236954\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd3wT5/0H8Oe0TpLlLVneQ/KQLMmbPRIKlKYhYWSSQFoSMpq2tE1oJmQ2SWlo0l9KQpJmj6ZkkYSkhTgbMAYMtqZly3tbw0Memnf3+0PUccGAjSU9J+l5v/rqy9G4+1jI+uqeiVEUBRAEQRAkUjFgB0AQBEEQmFAhRBAEQSIaKoQIgiBIREOFEEEQBIloqBAiCIIgEQ0VQgRBECSioUKIIAiCRDRUCBEEQZCIhgohgiAIEtFQIUQQBEEiGl0KoVqthh2BFgiCIAgCdgp68Xg8sCPQC0EQJEnCTkEjFEV5vV7YKejF6/Wi5TOnjy6FsKSkBP2zAQDGx8fHx8dhp6ARkiQHBgZgp6CXkZERl8sFOwWNeDye4eFh2CnoZWhoCH05mD66FEIEQRAEgQIVQgRBECSioUKIIAiCRDRUCBEEQZCIhgohgiAIEtFQIUQQBEEimp8LIUEQmzdvnnzLyMjI9u3b161bt2PHjpGREf+eDkEQBEFmyZ+F8OOPP/7d737X1dU1+ca9e/eKxeK9e/cmJSW9//77fjwdgiAIgswey4/HkkgkqampO3bsmHzj4cOHH3vsMQ6Hs2bNmoceeuiWW24519O3b9+OYdgZNyqVytWrV/sxJM2NjY0BAM5+HUiKNDus3SO9XaO9Mbhgrrgsis2HETDYSJIcGxsbHR2FHYRGxsbGPB4PWoFogtvtHhsbw3EcdhAaGRsbYzKZbDYbdhD4+Hw+g3GBSz5/FsKSkpKzb7TZbGKxGAAgFovPv0SI1+s9uwBE2kJBvl928q/sJtx7TZ9+0vzvGE5MmiA5XZCqtlpfUL9eEJ+7KHXu8oylbIY//xHphvov2EFohP6viYtwf9nxXbu90+YcsDoGPKR3jrhkSdr8vDgJBs78G589+r8gwYdekwnTeREC/hlKUZSvvFEUdf4FEv/85z+fXQgjU3R0tO+HU32avx5/IT9B+q81/0jgxU88wOl1Vvec/KKp8tvuI48uuUfEF0JKGnAkSTqdzokXBAEAEASB4ziPx4MdZAokRX3Z+u2r6nfkwvyK5BIhPzGJLyQp8nDXsb/Vvewm3BsK16/N/7l/T+p2u8GkvxoEAOByuaKjo9EV4TQFvBAmJiaazeb09HSr1SoUhu3ntd+RFPXsiT3He079Ye4d81MrzriXy+JemrnoksyF/zLsu/3A3Q8uvKs8uRhKTgSZYLSZnj62m8fiPrrk3kJhweS78hOkNxfdYBps2Xn0OZ3FuG3er7ks1JKJ0EUAp0/4NpSYP3/+wYMHKYo6ePDgwoULA3e6UKfRaFavXp2ZmZmZmbl69eqHDz7VNdLzxurdZ1fBCRjANhSuf2jxH588+rf3DB8HMy2CnEFrqb/vu8duUFy9+6c7z6iCE/LiJc+v+guTwbjz4B+7R3qDnBBBziWAhXDbtm0AgE2bNrW0tGzYsKGtre3GG28M3OlCWnV19bylP/ki/rKhR4xDjxg1ZfIDmq/X8FbyWNwLPrckSfnyz57Z33Twy9ZvgxAVQc6mMet3/PDk9oV3L89acv5H4kzO/Qt+vyb/sl9/ea9psCU48RDk/DCa9KZiGEaSZMT2Ea5cufIr8Rpw6e0AgCTWD9msf6q/lS/pOFhZWTnNI7QPd/7uqweevGT7ub6MhyiSJC0Wi2+8FeIzNDREqz7COrPukUM7dyzaNqP2+e87qp4/9drLl/01Do+dZQC32z0yMpKYmDjL44QTq9UaGxuL+ginCa0sQwvV1dWgdA0AYJVzbwHzFZ37QZfi+urq6ukfISs2474Fv9tx6M+WcWvAYiLImVqG2h85tPORxffMtJf6ksyFq3KWPfzDTi+J5oEgkKFCSAsYhgGKYmBuluDTu5rGFg/1AJKY6fXx/NSKa2RXPvj9k04v2rUVCQYP4Xn8yK7bS39ZIlZdxNM3F93AZ/N3n3zF78EQZEZQIaSFRYsWgVMfS5h7ZSPuv8Y/9Ofe56817F60aNFMj3O9fF1WbPqLtW8EICOCnOll9duZMemXSZZf3NMZGLZ90V0n+9RfNE+3CwBBAgEVQlp44oknEk79NY36XNafc9ibdG2P/G7myWc3Xsx0q60Vt/3QWdVga/J7SASZ7GSf+tv2w3fPvXM2B4li85+45MGXat80oyZ9BB40WIYWKEDdsf/u1FMm61dtHzXbFi9e/OS9dwsOvMwrWhR7xc0zPdqBlm8+bvj8xZ/tYmAh/0UHDZY5Gx0Gy4y6x27+99a75945L7V89kd7Q/te+3DXw4v/eHFPD9xgmWefffatt97y+2GDgCAIBoMRZp+oy5Yte+aZZwJx5HBenSuEVLZ+R7GoWzy86H99/mRMgm+NDCI/3/rig5THFbfuDjCTN/QqybJ/N1fubzq4Ju+ygEVGItqzJ15cmrHAL1UQALCh8KpffP7rl/a/2lPTwWAwFi5cuHLlSr8ceZa6u7uXLVu2adMm2EEQ8N13333//fcBOjgqhPCNusderH3j0eyrOcJvGDEJE7czo+NFv/6L9eUdgx/8Pf6a306/FmIA+8PcX/3hqweXZiyM5852bDqCnKHe1qg269+98kV/HdA97rLu7/xHYdcpzRyKIMDLv/95SfYHH3zA58NfWT41NbW0tBR2CgS0t7cHrhCGfNNZGPi48fO5qeXpzW28ojNHxzD4AtGdTxG2voG3/0wR3ukfMyc282eS5XtqX/drUgQBAIAXa9/cXLQBZ3L8dcAdO3Z8XctzCytS15SBK3eAh2v+3QUee+wxfx0fQc4PFULInF7nxw1f3Chf79RV81QLzn4AxuEmbnmEdLlsrz1OedzTP/IvVNef6tM0DKBRM4g/VffUDDmHLnqk6JT27t0L1jza7L0lk/VhAtkDGCyw5uG9e/f68RQIch6oEEL2WdPBUrFKPDjC4EezRGlTPgZjcxJv3s7gRllf2k65HNM8Mo/F3VC4/m0d2gwZ8RuSol6ufevWkpv8OA6LoqixoYHLeX1/bv/XzyxjK4hnAQBAmG02m/11CgQ5P1QIYfIQnr31n2woXO/QVvGKzrciOcZkJWz8I1ucaXn+PnLMPs3jX5H703qbqXmozQ9ZEQSAg63f4Cx8UfpcvxyNHLOPHfvS9sqjNbdcsqHv0+NRimrstuGoNjY2DLq02dnZfjkLglwQKoQw/afl6/x4SX6C1KGp4qkuNH0ew+Ku/jWeX2L5+x8J+/m2OJ7AYXKukV35ru4DP2RFIp6bcL+uee/OsptnubkuMWAe/f4Ty+57+v50s7P+BL/80s/TL9m478Q7+NxKwbLyQSzX+0+wd9ttt93mr+QIcn5o1Cg0BEW8Z/h4+6K7PL1tgPCy06UXfg6Gxa7ezOBFWZ67W/irp1iJyRd8xtq8y6433NZh786MmbrdFUGmaX/Tl3nxOSqR/OKe7ultc2iqHNoqYsjKU86P/snVeF4JxuYAAH5fsrTHOvjcI8We3MUJpR7xnINXbFq1detWv8ZHkHNChRCar9t+SBYkKYQy+8F/clUz2Kkxevm1DG6U5e/bhHc8wU7OOv+DuSzuVQWr39V/cP+C388uLxLRSIr60PjZjkXbZvY0inK1GpzaKoe2ClCAp1oQt/4OPLsQMP6nLYrBYOzatWvr1q3Hjx9P6m80c4wVGxeG2WRwAIDL5Wpvb2exWBkZGWhfCFpBhRCa9wwf/6Z8CwDAoT0St+6OGT03atHlGC/K+sJ9iVse5WTmn//B6wtWb/j0tp7RvlTBha8gEWRKVd3H47lxhcILvNl8KK/H1Vjn0FY5ddWMmASeakHi5h3sNMn5n+XblZoYsjJ3/+r/Yj9dl385i8H0R3b4WlpaHnrooX37vxjniwFFxroHrr/26kceeSQ5Gf1J0gLqI4RDa6n3kkRZcpF3oJ8YsuE5ipkegV92afx1v7e+/JCrSXv+R0ax+Wvzf/5P/UcXGxZBwAfGz66RXXn+x5DO8fHa7wfefKp3x4aRr99nizNEv39W/MfnY3628YJVcAIzTpjLTkzjxH/VFqjZ00H2ww8/lC1a9i62YPzJVvC4FvxJP/xI/Uu27NK5C7XaC/zxBtT0l+85/yNpsgzQbKBCCMf+poNX5P4UA5hTc4SnnH9GS9E0cRXzEjc/aHvjTw5t1fkfeVXB6m87Do+4Ry8qLBLRXC5X02BL72jf0swp5rkCAMgx+/iJr6z/eLjvkY1jRw9wsuXi+18W/fZpwaXrp9ONfTaurHw9mfae4SMK0GIl5Nno7e1du+Gm4Ts+AZfcBvCo07fy48Cqu/s2vLp67frR0XD+q7znnntgR5gWVAghGHGPHuk69jPJcgCAQ1N19oIy04dLVcLb/zT0wd/Hj59vI5tYPGZB6pwDLd9c9ImQSDM8PPyHP/xBJBJxBbHX/emmhF4B6SUnP8A70D/6/SeWF+7ve3zzeN0hfsmS5EfeEd35lOCStcxJKwVeBFxWXtDcxWawT/aqZ/dLwLdz587BBbeBNOUU9+Ut7si57IUXXgh6qOCpra2FHWFaUB8hBAdbvl2QOicGjyZHhzy9bXh+yWyOxsnIE975Z+uLD1JeT9TCc+7cdGXeqqePPX+17IpZjn1HIoHX612xYkUNuwDcc4KdIOCyf/Pq3zRjVTe//fbbnq5mh7bKoa0iR4e5pwd/FmNMf36S4BKlt7d19cqbvmiurEiZ1V8HdPv27QO//vacd8+7bt++bdO/bDpw4MBrr73m9Xo3bdq0bt26qqqqN998s7+/n8PhXH311ddee+3KlSuvv/76gwcPXnHFFRRFffPNN4ODg9ddd90NN9ywcuXKLVu2fPDBBxkZGffee+9E9+TAwMDu3bs1Gg2Px1MqlVu2bElMTLTb7c8///yJEycEAsG6devOiDHlvWeHeeihhwAAt99++y9+8Ysz7prpyxho6IoQgs+bv1yd91MAgENXzZWVY6zZjh9jJ2eJfvP0yNcfjHzz4bkeU5SkYDKY6n7dLM+FRIKPP/64xoaBX74C4tNSWQcs1NKK6x+X9hk6d9xge/NJyuuOv3ZryqPvxl+7lSsr928VBABgbA5Holjiij7RWzvkGvbvwYPJ4/F09vSB+PRzPkKY09zcPP0Dvvjiizt37nzuueeOHj0KAHjjjTeWL1/+0UcfPfnkk6+/fnph4ZycnJ07d7711ltRUVH/+Mc/HnrooXfeecd3l8Ph2Lt3r1Kp3LNnz8Qxn3nmmaVLl7777rt79uxJSUl59tlnAQC+69S33377pZdeMplMZ8SY8t6zw/hWi33ppZemzEkr6Iow2LSWeoIkipMUAACHpoo/xz9rNrKEKaKtu6x77qccYzGX/2LKx1yZ+7NPTQdKxCq/nBEJYydOnADKnwEGk0uN52OfbGtgEKThP4LUY+nzr749GNP7uLIKT4NmYdbcytbvLzhIh7ZYLBaHxXARHsA6xwLlHueMdtgoKip69dVXV65c+dRTTwEAXnzxxYaGhsrKSrVa7fWeXpR/yZIlLBYLALB27Vomk1lWVubxeHx3rVq1islkXn311Zs3b544plqtPnbs2MR/xsbGAgCOHz/+yiuvREVFAQC2bNlSWfk/PS9T3jtlGJ/z3EUT6Iow2CaGyVAuh7tFx5XP8deRmbGJot/uchprhj56AUy13/IqybITvbWDziF/nREJVywWCxBuAEApsT/Vgd2c/refS5/9ewvlEvh/89spceUVTuPJ1dKVXzSfr/Ob5jAMKy0pAe0nz/mItpqKiorpH/DRRx9ds2bN4cOH77//fgDA448/vm/fvtjY2FtuuWXiMWw22zcFk8lk+jKcnYokf+zuFQgEb731VmVlZWVl5Weffea72pv8rCmPcPbPU4a54F00gQphUI26x450HfupZBkAwGE4zpEoGFx/7rjGiIoR/nqnp6d14J+7AEmccW8Um780Y8G/m7/y4xmRsHTJJZeA2k+BxyliHWI5Fd1sERjqYbZWL168ODgBWKI0jMmWkwKSIg3WxuCcNBC2bNkCvvzb1PdRJKh8bka1YePGjSkpKRs3bvQ1SJ46deqGG26YP39+TU0NAIAgzvyTP8OBAwcIgvjoo49Uqh+bhRYvXvzee+85nc7BwcGHH374vffeAwDMnTv3pZdeGh8fdzgcr7766hnHmfLec4Xxer0zzRl8qBAG1cHWb+enVsThsQAA5+zGi54Lg8sX3vEEOWa3vf4E5fWcce+avMv2Nx0kp7peRJAJq1atunJOPvvVNQ52v2m0DKi/AM9eds9vbsvKusBKRn7ElZU762sukyz/ovnLoJ3U7375y18uj7eDzx4/s5GG8IB//u7GedmXXXbZ9I923XXXbd26ddu20wux3nzzzXffffett95qt9srKir+8pe/nP/pBEFcc801dXV1d95558SNmzdvJkly06ZNW7ZsEYvFviP/6le/oihq48aNt912W1FR0RnHmfLeKcPMnTv3pptummlOCCh6AACQJAk7RcDdcWDbsZ6TFEWRHnf3fVcRI4NnPMBut9vt9tmfiPR6bG88aX7+PtI5fsZdt/7nD8e6T87+FMFBEERfXx/sFPQyODg4Pn7mP6vfeb3e37+87f5nf56fk7V8+fJ9+/YF+oxnGNcetbxw/4BjcPX7N4x7HOd5pMvlslqtgchw9913P/3007M8iN1uX7t2LcguBxt3g3u/A/d8A67bBVLlt9xyi9Pp9EvO6VixYkXQzhUI+/btW7NmTYAOjq4Ig6dntK93tL88uRgA4GqsZafmMARxAToXxmQlbLqXlSC27HmAHP+fGbuXS1cebD33eG4EAQAAwGQyh4T2lXh2Q0vbV199tXbt2iAH4OYVu9rq4xi8YrHym/ZDQT67H0VHR+/bt++rV3Zu4dYtPrJ9ydGHfyNqq/749VdeeQXHcdjpEADQqNFg+qbt0KWZi5gYEwDg0BwJRLvo/2Aw4q/73fD+Vy3P3SX81VPM2NPDHC7NXPRy3VtOr5PL4gY2ABLKmofaRpzDpTmrYQXAcB4nI8/VpLlcuvId/fuXS0N7Ha/ly5cvX+6fIeIX54EHHoB4dppDV4TB81X7D8uzlwIAAEk6dMe4yvkBPyWGxV65hT9nheW5bV5rr++2WDxGKZQf7jp2/qciEe5AyzeXjHJ4ucUQM3Bl5U5jzdzU0u6Rvr4xtGH9rCxbtgx2BPpChTBIWobaxz3jSpEMAOBq1TPjEi9uGcaLEL38WsGl6yzP3+vt7/Tdsjx76ddtPwTn7EgoIiji69bv5neOcLJlEGNw5XOc9TVMjLkkY/53HUcgJkHCGyqEQfJ12/fLs5b6ljdzaI7wioI0DN1HsOTK2NWbLbvvcXc0AACWZszXWuqHXfZgZkBCyPGeWjFTkJ1cMPtlj2aDnZpDuZ1eW++yzMUh3U2I0BwqhMFAAeqb9sM/yV7i+0+nrppXNIOdeP2CX74s/vrf2V55xNWi57K4c1PL0Fds5FwOtn6zjEjA82C2iwIAAIbhBWVO48kSsdIybusZ7YOcBzbfhkdhsO0R3aBCGAz11kYWg5UXLwEAeLqaAIN5wZ3lA4GrmJ+w8R7ba487jSdVPNlr37+9cePG++67T60O+TX+ET9yel3He06VtQ9woXYQ+nBl5a76GgbGuCRj4Tfth2HHoYXrrrsOdoRwgwphMHzV9sOK7Et8Pzs0R4J/OTgBzy8Vbnm499U/vXXrXRaG46P4xTuN3PJLf7Zz505YkRC6OdZzUhYv4Zv72JnT2o8+oLgFZa4mLUV4f5K1+FvUOgoA8K1Wg/gVKoQBR1LU951VP8k63SkYjIkT5zUSnbT+g+P3Ly2JIwpE5Ylg9YPE/Uce+NPTcDfLRujj+86qRew0To7C79tKXARGVAxLnOFuNShFhUMue/twJ+xE8E00kH777be33Xbb+vXrP/roIwDAwMDAY489dvXVV2/atGnnzp02mw0AUFVVdfvtt69du/baa699//33fU/8+uuvr7rqKri/Ba3Af6OHPa3FkMCNy4hJAwB4Ld2kY4yTWQAxz9dff10nUFyd98yTvfe9ndnbCdaChAyy4ur9+/dPXoEQiUwewnO859RN2Bz4HYT/5VuAG88tWpa56LuOI79QXQ8xDEV4hz/9R3DOhbHYsVee7+LPbDa/9NJLdXV127dvv+qqq5555pkVK1bce++9Ho/n448/fvbZZ//0pz+98cYbK1asuOqqq1pbW3/729/69gLU6/W7du0Kzm8RElAhDLjDndWL009PGXRoqniqBeCs1dyDaXBwEMQmd3LE94r+Ug5+lUNoW5kqEJcyNIR2pUBATZ86Jy6Ld9KI33g57CyncWXlgx/8PXb15mVZi5+u3g23EAIAWMLU4JzoglfkV1xxBYZhpaWlbrcbnGNDpSm3QLrpppvi4gK1rFUoQoUw4A53HXvykgd9Pzs0R2Iv/yXUOKCgoAC0/R+gyB626LKxJDnjy1amCrSeyL882GtoITT0fceRS8SlhP19Tlou7CyncTILiCErMWwrFBY4vM7W4Y6c2ExYYTAmS7B0Dayzn+GMvQwFAsHLL7+ckpICAHA4HCMjIwCAxx9/nM1mL1u27JZbbvnqq9M7z6AqeAbURxhYzUNtAICcuCwAADFs81q6OVIl3EhLliwpFXPBx9sB4RnxlrI5OvDdSynmU9dccw3cYAh0BEVUdZ+Y4+ThUhVg0ObDgcHA80tcDacwgF2auei7djTtZ2pTbqhE/y2Q6IA27/Uw9UPH0aUZC3w/O7RVXMU86AMQmEzm/v37V7OM4L5c7b8+dTHtC3rfq6ys9DWkIJGstk+bFp0c09aK02DixGRcWbnTeBIAsDhjXlX3cdhxaGrKDZVCYAskGkBNo4F1uOvY7+fc5vvZoTkiWHIl3Dw+aWlp+/fvb29vbzKZKht3L/y/bQqFAnYoBL7vO6uWZix01XxBkzfqBK58zvCnrwCSVAhl5jFr/5hZHJUEOxQElZWVZ/z/5Nv5fP62bdvOeMqaNWvWrDndlusbKePb3R6ZDF0RBlD/mNnmsCmEMgAA6Rj1dJi4BWWwQ/0oKytr+YoVC7iZh1q+h50FgY+kqCNdxxbHF5LjI1AWfDgPZkwCMzbR3dnIwBjz08qPdtfAToSEFVQIA+hw17GFaXMZGAMA4NRV4/nFGId2Ox/Nz5yvG+1weJ2wgyCQaS2GeG58YlcXnlcMd2DzlHyTKAAAi9LnHUGto4hfoUIYQIc6qxdnzPP97NBU8VQw59GfS3xuWb6TdaznJOwgCGRVXceXZMxzmdT0mUE4GVd2uhDOSSnVW4zjHgfsREj4QIUwUOyuEdNgS0VyCQCAcrtcJjVXMRd2qCmwM3LLrK7DaCRexKvqPrEgbY6rScPNK4GdZQocqdLb206Oj/BYXKVIfrz3FOxESPhAhTBQjnQfr0gu4TA5AACn8SQnq4DBj4YdagoYkzU/Ou9o90kvicZVR67e0f5Rz5iE4lNeDyspHXacKWBMFkeqdDXUAgAWpc+t6j4BOxESPlAhDJQjXccXpf+3XVR7hAtvoe0LSskpScO4arMOdhAEmqPdNfNSy90mLU7Ly0GfiUkUC9PmVHfXkBQJOxESJlAhDAgP4TnVp56fWg4AoAivU3+cp5gPO9Q54RJluZ15uOvYhR+KhKnqnpr5qeW07SD04cornMYaQFEivjA5KklnqYediHbQhoUXBxXCgNBYDDlxmTF4NADA3axlidKYcULYoc6Jky0v6Rqo7kJtTRHK6XXpLPVzUkpdTWo8rwh2nHNiCVMxNu7pawcALEyfc6QLjR2dGtqwcKZQIQyI6p6T81LLfT/D3YBwOjAOnhOf7fW6ukZ6YGdBIDjVry5IzMUHbIDJZCWmwI5zPlxZubO+BgCwKH3eoa5q2HFoKlw3LPR6vY8//vhVV121efPmd999l6Iofx0ZrSwTEMd6Tj648A8AAEBRDu1R0W/ovqwRLlGWevXHek6mFwRpZX2EPqq7T85PrXA2qemwJf35ceUVI9/vi/7J1bnxOR7C02nv9m1wFjRekvioYX9wzsVisK4qWD35lgMHDrz22mter3fTpk3r1q1buXLlli1bPvjgg4yMjHvvvTc5Odn3sJUrV1ZWVq5cufKBBx547733rFbrjTfeeNVVVw0MDOzevVuj0fB4PKVSuWXLlsTExOD8Ln5x7NjxL0ZTQeH1wGl/46G/v/76659//jmX64fJ2eiK0P96R/tH3KP5CVIAgLu9gcETsERB/Vu9CLhUqRrwHO+phR0EgaC65+SCtDkuk5rOI2V88LxiT0cj5XJgAJubWna8N9jvWApQVsdAcP434Bg44+wvvvjizp07n3vuuaNHj/pucTgce/fuVSqVe/bsOTutb8PCHTt2vPbaawCAZ555ZunSpe++++6ePXtSUlKeffbZQL9c/uVJyAIbnwdla8HCm8B933/dTbzwwgt+OTK6IvS/6p6aeSllGMAAAA4t3dtFfTg5Svm7PXtirC7CjTM5sOMgwdM81MbAsMzo1N5mXdza22HHuQCMw2Vn5LuaNFzFvDkppQdavj7jminQ2AzWr8tuDuYZJysqKnr11VdXrlz51FNP+W5ZtWoVk8m8+uqrN2/efPbjp7NhYSiJndRujzHAsjs//3zPXXfdNfsDoytC/zvec2pSB2EVVxUChZDBF0THJkn4yWgSRaSp7q5ZkDbH09eOcfl0HtI1gSuv8HUTViSXaMwGD+GBnSh4Hn300TVr1hw+fPiMhbMxDCPJKSaTnL1h4VtvvVVZWVlZWfnZZ5/563IqeDAMALDRs1XuqQEAgJgku93ulwOjQuhnHsKjNuvLU4oBAJ6+dorwcNLpssHp+eFSZSkVdwK1jkaYo901C9IqXI11dJ44MRlXVu6sPwEAEHCismMztZE0iWLjxo0pKSkbN240mUy+Ww4cOEAQxEcffaRSqS749Ck3LAwlHicbeCxR3f3MdAAAaJMp0GUAACAASURBVKqSy+V+OTAqhH5Wa9ZK43NiONFgYn1R+q1fPCWORFlkdR5DK1dFkhH3aMtQW0mS0tVE6xmEk7FTcyiC8Fp7AQDzUssiaq216667buvWrdu2bfPtNQgAIAjimmuuqauru/POOy/49Ck3LAwlvfXK4S8EHtYAIxnUf8M++PTvfvc7vxwY9RH62bH/nTgRt47unS4T8Nyi1H17RmO5vaP9KQIx7DhIMJzq0yhFcg7GcjXr4q7dCjvOdHELypzGGsHiK+aklO06tvtmxQ2wEwXJlVdeeeWV/7NV5JYtWyZPlriIDQtDSLY4Ia/hGXZmHPh7YW4M2P3hexUVFX45Mroi9LPq7pO+BWWIATMxaMZzCmEnmi5mTAKTwy+NyzsR9JF4CCw1fXUVKSXu7iZmbCIzOh52nOniysqd9ScBALLEXIvDZnMOwk6EBEN2dnbeYtmluXNMR/5jNBpXrVrlryOjQuhPPaN9Dq9DGp8NfONFVQsAgwk71AzgUmWZNzr4Q9IRWGp66+Ykl7hMmlBpF/XBC8rczRrK62FgjIrkkpq+OtiJ4HjggQdgRwi2BnJw+cJ1ubm5TKY/P1pRIfSno90181Mr/jtx4igvFMaLTsaRKhX9o7X9GrQTRSToHe13Eq7suEyXSY3Tfir9ZAy+gJWS7W7RAwDmppSejNRCuGzZMtgRgsrjdowzqOxM/892RYXQn2p66ypSSgAA5Jjd092C55fCTjQzuESJN9dnRKfpLAbYWZCAq+mrq0guxkjS3WrApUrYcWbG100IAJiTUnrKrCX9t9oWQluj40P5WKzvSsO/UCH0Gy9JqM268uRiAIBDexSXlWHsEJuZzhKlAYqqiM8/EalfsSNKTW9dRXKJu72BKUxhRMXAjjMzE7MJhfzEBG5cs70VdiIk4Mbco8o4SSCOjAqh3xisxvTo1Fg8BgDg0FaFXLuoD0eiULq4p/o0sIMggUVS5Kl+TXlysbNJzc2l744T58LJLCDsA8SwDQBQIS6ptaCFIMLfGOkuypwTiCOjQug3NX1q3+Ug5XK4m3XcwrmwE10MPEch6R1oG+4Y84zDzoIEUONAcyI3XshPdJnUeD7dlxidAoZxC0onWkfVNj3sQEjAuTBSIbskEEdGhdBvTvTWzkkpBQA462s4OYUMLv+CT6EhXKoiWwzyxHy01lp48/VnU16Pu72BIwmxDkIfrux066gisaBtpGPc44CdCAksDsXg44JAHBlNqPePUfdY61C7QiQDADg0h3nFi2Anukjs1BxyZKAsYempPs3CtJC8qEWm40Rf3YbCde62enZKVqh+aZNVDH3yMiBJDpMjicnWWgwTa1n4kVar/fDDD/1+WGSmTpw4EcXEA3RwVAj9o7ZfoxIV4kwO5fU4jSdj190BO9HFwjBOdqHKzf2b9diFH4yEJqfX1TjQVJykdFV+EFoTJyZjxsQz45Pc7UaQlqtKkJ/q0/i9EG7duvWuu+7617/+5d/DBoHb7Waz2ViIrO84HaYe3ZzSQHVmo0LoHzV9at/ECVdjHTs5O4QW6TgbLlVm9NqsnoEBx2ACL4R/EeRc1GZdfryUx+KOmuqiV90IO87F48rKnQ2nuGm5RYmFb5j8X64yMzND9HLQarXGxsay2WzYQfyEJNe+ve6xy3cG6PCoj9A/JmYQOjShsQHheeASpbdZV5RUWIu6CcPU6Q5Ct8vd3RJCqwCejSsr93UT5sXldI/0Drv8sykPQjddrSdJDEsV5gTo+KgQ+kH/mHnMMy6JywIU5dQf46oWwE40K+zMfI+5qzRRfqpPDTsLEhCn+rXlycWuFh0nTYpxuLDjXDyOROHt7yTH7EyMqRIV1vWjr27hSd1cVchKDNzxUSH0gxO9dRXJJRjAXK16RkwCKzHlws+hMYzJ4mTkqbx8NJswLNndI72jffkJua4mTUhOnJgEY7LwXJXbVAcAKEsuqu3Xwk6EBITObFAkFgTu+KgQ+kFNX51vJ16HpirU20V9cKkqpbvfTbh7R/thZ0H8TN2vU4kKWQxmCG3Gex64rMLTWAsAKEsuOonaMMKU0W0pkgSwpQ0VwtkiKaq2T+ubSu/UHuUVherEiclwqdLVoisRq071o4vCcFPbrysRKymXw2vu5GTJYMeZLa68wtVwElCUNC7H7hqxjttgJ0L8bGygp5tNyHMCOJsLFcLZah5qjeXGJPGFnu5mgGHslGzYifyAky13dzWViRToK3b4qevXlopVriY1O7MAY4X8qEJWYjKG88n+DgaGlYiV6Ktb+NHWf5sNojisQE0iBKgQzt6pPk2pWAXCqF0UAIBxuGxxZhEhONWnoQBa1z98DLvs/eOWvASJ06Thhn67qA+noIxs0QAAypKLUMd2+NF118mjMwN6ClQIZ6u2X1MmLgKnJ06EQ7uoDy5VxnV18Vjc9uEu2FkQv6nr16lEcibGdJnUeF5oj5SZgBeUEc0aAECpuAhdEYYfw0hHUXpZQE/ht0I4MjKyffv2devW7dixY2RkZPJdW7duXflff/vb3/x1RjogKEJrqS8RK73WHnJsOAx6XCZwJEpXi65YrESLjoaTOrO2RKwix0cIWx87Ixd2HP9gS1RkTwvpHM+MSSMpqme0D3YixG9It6OR5SjKXxrQs/itEO7du1csFu/duzcpKen999+fuJ2iqK6urr1793722WefffbZnXfe6a8z0kGDrSk5ShyLxzg0VVzVQhBGCxrhUqW7zVgiUqjR3KwwUtuvK01SuUxqjkSBMcNkYSmMgzPT81xNvotCFZr/Gk5ajEe4GEsUkxzQs/itEB4+fHjNmjUcDmfNmjWHDh2auN1msxEE8eCDD1577bVPPfXU2NiYv85IB6f6NGXJ/20XDc0NCM+FwY9mxgkVVHQduiIMF0OuYcu4NS9B4jKpw2DixGQMaZHLWAMAKBWr6vrRlkzhQ9NareAGfGa2374S2mw2sVgMABCLxQMDAxO3DwwM5OXl3XHHHUlJSXv27HnhhRcefPDBKY9w++23n71EbFlZ2fXXX++vkH53vPvUFTmrBrvaPf2dzqRs1/DwLA84OjoKACBJ0h/pZgtLz+c1NzMAw9jTmBIlhpKBJEm73c7lhvDqJ35nt9txHHe73TN9YlXvcXlc/oh9xNFYy13/m+FZv11pwu12O5KlxOcvYiuHpfysV/reCZtf7aLZ7XYAQBisNaodaMxLK5nNP6hAIGAymed/jN8KIUVRvjJGUdTkz/H8/Pxdu3b5ft6yZcuWLVvOdYSCgikWDkhLS6Ptv6WH9DYONZeIleDUIXZ+GYfLm/0xfb8sTX5lUqLw6KtV+fL6YVNmXDqcDCTJZrNp8oLQBJvNZrFYF/Ga6AcbSpIUTNcYOTKIp+cCRpiMlaMoiiXOBBTFGLZkCNMYGGZzDyZHJcHOBRP7v2AHmR2KaiSH1+Uumc0vMp0tOPxWCBMTE81mc3p6utVqFQqFE7c3NjZ6PB6FQgH++29zriPcddddobVpSJ1ZlxOXKYoVWhtqohf8jMf3w6ZuBEEAAPj+ONTs4YUV/V+8Xn7JL9WW+nXyy6FkIEmSx+PR5AWhCbfbjeM4jzfjL146W/1a2c8ZrY3c3GK+ICAbnELBYrEIgmDIyrE2PT8zt0isaBhploiyYeeCaXx8nM/nh3ohHOyst+FAkVPGxC5wSTdLfvtKOH/+/IMHD1IUdfDgwYULFwIA1Go1AMDpdD7yyCPt7e0ej+edd95ZtCh8JhjU9mnLxEWkY9TdZuTKKmDH8T9mbCID5ymZwjq0hGPoG3QO2xyD0rjs8Osg9OHKy53GGgBASZISjfAKD9qG7/MYsYGugsCPhXDTpk0tLS0bNmxoa2u78cYbAQDbtm0DAKhUqhtvvHHHjh3XX3/9yMjIeZpGQ86pfk2JWOXUH8fzijDcD+2iNMSRKET9FoIi+sbMsLMgs1LXry1OUjAwRrgWQjy/1N2ipzzuErESjfAKD9o+nTxOGoQT+a1pVCAQPPHEE5NvqaysBABgGLZ27dq1a9f660Q04fS6TAPNKpF87D9Ph9l40clwqdLVrC3OUqr7dcmSn8COg1w8tVlfnKQghiykc5ydnAU7jv8xeAJWSo67VZ+RX+L0uvrHLOIoEexQyKzUO3o2ZAejUyZMesuDT2sx5CVIcYrhaqzjKubBjhMouFTlbtahr9hhQGsxFCUpXI113LzicJrwOhlXVu40nsQAVpykUJvRJIrQ5hm2NeFelSQYvWmoEF6k2n5tmbjIaTzJzshjRMXAjhMoLFEa5fUoualoy9OQNuYZ7x3tz0uQOJs0Ydku6sOVV/g2rC9OUmpQIQxxjcZvExncGG4wPl1RIbxItf2a0mSVQxs+C22fCydHkWy2ObwOy7gVdhbkImnMBrkwn4kxXU3h2UHow8nII0aGiEFzCVoaMPRpO07J+WnBORcqhBfD4XW2DnUUxkud+mM8VQC3i6QDXKp0N+tUokLU1hS6NGZ9kUjhtfYCgmCJgvThAgGGcQtKnQ2ncuIyh1x2q2Pgwk9B6Epvb1OlBulLGyqEF0NrMRQk5lGtDazEFGZcmHfI4xKlu0VXIlah1tHQpbHoi5MKXaY6PD9Mdpw4l8ndhKh1NHRRHrcRGy3KC+xa2xNQIbwY6n5dSZLSoTkS9u2iAAB2moSwDygFmRqLAXYW5GK4CHfzYJssMd9lUuO5Ydsu6oPLyl2NdYAkitB4mVBmbqoZZzOzE3OCczpUCC9Gbb+uJEnh0FWH8cSJH2EYJ0uePjBicwwMuSJ9CcdQVG9tkMZnc1m4q1kbxh2EPszoeGZCkrvdiAaOhrS65sNyTiIGgjS8GRXCGXN6nS1DbblOBoPLY4kzYMcJBlyq9DbrCxMLdBYj7CzIjKnNBpWo0NPfgbHYrAQ4i6cHE1dW4TSezIuXWMatwy477DjIxTBYG+WJU6w+HSCoEM6YzmLMS5CSuhPcSLgcBACcnlavU4rkOks97CzIjGks+uIkpcukCft2UR/fJAoGxigUFmjMqD0/BFFUvcdSIg3eepyoEM5YnVlXkqSMhIkTE9gZ+R5zpypeij5WQg5BEfXWRoWoIFxXVjsbJ1vutfSQo8Mq9NUtNI33trTxQGFGWdDOiArhjNX165TsJMrl4KTnwc4SJBiLzcnIlY6QzUNtTq8LdhxkBkwDLclRSTFsgas5nKfST4YxWXhukbPhVJGoEI3wCkWGhh/SsCgeK3i7kKJCODMuwt002CLp6uMVLQrXdaqmhEuUWJtREpfZMNAEOwsyAxqzvlis8PS0MvjRzNhE2HGCxDeJQi4saBlqR1/dQo62R62IzQ7mGVEhnBm9xSiNz6a0x3lF4bOf1HRwJEpXs1YlKtSi1tGQorEYikQKl6mOmxfmMwgnw+UVLuNJnMGWxGUZbY2w4yAzUz/epUovD+YZUSGcGbVZXxwj9Q704xIF7CxBhecUujublAl5WtTpEjooQGnN9aqkQqcpUtpFfVgJYowX5eluQa2jIYcYGWxku4uCOFIGoEI4U3VmXYHdy1POB4yA7xVJKxjOYydlyNy43mokKQp2HGRaOoa7eGyuEI9zt+o4UiXsOEHFlVc4jTVKkRzN+QktnY1HAYuVEp0czJOiQjgDHsLTYDPlmFojZ7zoZLhUyetojePGtg63w86CTIvWUq8SFbq7mpixQmZ0POw4QcWVlTuNNUVJhTpLPUmRsOMg06VuO6bgpgb5pKgQzoDB1pgTk8bqbMXzgzeulz44UqWrWacSyVE3YajQWY1KkSxyJk5MhucWe7qaoylWIi+hZQh9dQsZ9UPNheJgdzyhQjgDdf26QioWl5VhbA7sLBDgEqW71aASylE3YajQW+qVIrnLVIdH0kgZH4zN4WTLXaa6oqRCLeomDBGUx22k7MW5i4N8XlQIZ0Bt1uX32SNtvOgERlQMMy5RTgnQEo4hYdhltzkGswWp7rZ6XKqCHQcCXFbuNJ5UitBXt5Bhb9P1cBkFSbIgnxcVwukiKMJoa5Q2t3HlFbCzQMORKEW9ZoIi+scssLMgF6C11CtEMqLTxBKlMfgC2HEg4MoqXPU1RWgrzdChMR2SMOPYTHaQz4sK4XQ1DjSLmTFxGXIGLxI/U3xwqcrVrFMKZWjlKvrTWeqVQnmkTZyYjJ2cSVGUeNxLUmT/mBl2HOTC9GaDMiE3+OdFhXC66vp18nFGxLaL+uBSlbtFpxTKdFY0JJ3utBaD6nQHYYQWQvDfJWbQbMLQQFH1rv6inPnBPzMqhNOl6dfndVu4ygWwg8DEjBNibI6MlahHc7PozUN4mgbbZLHZ7o5GXBJZMwgnm5hNqDWjNgy685g7m3hAlQmh7wkVwmkhKUrbr1NGZTJjImsy1tk4UmWWbbTd3uX0OmFnQc6pYaA5MyaN2dXKTsnGcB7sONDg+aXuVoMqPg815tNfU0OVgMmJ58YF/9SoEE5L63B7DMlIVi6BHQQ+XKIkW+tz47PrbSbYWZBz0loMKlGhq0kdUUuMno3B5bNTJRmDo31j5jHPOOw4yPlouk8qBZlQTo0K4bSo+3SyQTcvYnbiPQ9conQ3a5VCtNMbreks9UqRzNkYiVPpz4DLy73G2rx4icHaADsLcj4Ge7syrRTKqVEhnJa69qNyKpolTIEdBD5WUjrpcsr5qXo0XoauKEDprUZFnNTT08LJlsOOA5lvrTWlSIYWHaUzcszewHEVS+BcbKBCOC1am6lMEtHDZH6EYbhUmT9C6Cxo9W2a6rL3cJic2N5eTnouxsFhx4GMk55HjNrlvBT01Y3OLKaaYTaWHZ8F5ewsKGcNFa2trR9++GGztRXInNFZc2DHoQtcomC2tURzBJ32rqzYDNhxkDPprEalSB6ZS4xOAcO4srJcy6jB2kBSJAND3/7pSNN6VMYRw/rXQe+Jc3r55ZcLS+fec7DzBIOVPUzJLr3siy++gB2KFnCpytWiUwhletTpQku+qfQRPoNwMq6sgm3SJ/LiW4c7YGdBpqa3mRRBX1ltAiqEU2tsbPz1Xfc67/4GXP/XrBxqFCsavOmtTTf9Ynh4GHY0+NipEmLYqojNRm1N9KS3GBVx2Z7+Tk4WtE8WWuHKyl2NdQphARrhRU+U11NPDpYEdzPeyVAhnNonn3ziLbsKpMgAAAS7S8P6CZD/ZDBJ9fXXX8OORgMMBidLlj/OQB8rNDTqHusft6Rbh/FsOcYK9pqN9MQQxLKEKTIQi9ow6MnR0dDGB/JkaCs/oEI4NYvFAhIzAQAZRBPBIGrwxQAAkJBpsaDFpgEAAJeqUnvNlnGb3T0COwvyP+ptjfnxUo9Jh9pFJ+PKKqTWUfTVjZ6MpsMpzOgoNh9WAFQIpyaRSECnGgCgIL6Nd8aRvheqUy2RSCAnoweOROlp1ssS89DcLLrRW40KkcxlqsNzUSH8ES4vFzWa7K6RQSfq3aAdTa+2MDYHYgBUCKd27bXXJrT9AGo+TMD0Lo8UUBQ48LQsynXppZfCjkYLnMx8T3+HMiFXb0GFkF50FqMiOpOw9bMzIKziT1t4tpy09cjjJeiikHYoyuDsLc6GOSwfFcKpJSYm/vuLz+U/PO5itLcdUoOHixf2frF//342G3W6AAAAxmJz0nPzvTydFX2s0AgFKKPNJB12cyQKjIkmR03CYOK5xQVEFBrhRTdea08jnyzKnAcxAyqE5zRv3rwfjnw1ymPctXb9yc/fPXz4cG4u+or9I1yqlFpGjDYTQRGwsyCntQ13xuDR3BYT6iA8G1dekWcZRevL0E2XqRowmSkCMcQMqBCej850KI/kX79hY1lZGYZhsOPQC0eiZLU0JPGFrUPtsLMgp+ktRqVI7mpSc1EhPAtXVp7Z2GQabPEQHthZkB+pO2oU/HS4GVAhPB9110mlAM6SP/SH5xS6O02FCXloSDp96K1GuSCDsA+wU9GorjMx45N4vNh0PNE02AI7C/Ij/VCrMqUIbgZUCM9HN9pRkl4GOwVNYTiPnZRRwIhDhZA+DNaGvBECzy0CDPSnPQWuvEJGRulQNyFtkOOjDSwHrLW2J6C/lnNyepztmEMl+wnsIPTFkSpzh1xoBgVNjLrHzOPW5M4eNHHiXLiycollBL1j6WO4ua6Pi+UJ8+DGQIXwnDQN32d5WPw4Eewg9IVLlEkdXQOOwWGXHXYWBOitRllinhettX1uuFQl7bHozWioM11om6uknEQ2A/IIZ1QIz0ndfkzBTYadgtZwqdLbWi9LyDVYG2FnQYDealQIsking52MOranhrE56amFbo/TMm6FnQUBwLdTSmIB7BSoEJ6bZtBUlAy5C5fmGFExjJiEQl4yamuiA53FmOcA3LxigEY4nxtXVp5H8FHHNh1QhLfeayvKmQ87CCqE50BQRBNlL5FdCjsI3eFSpXQMoEnK0JEUZbSZsnpsqF30/LiyCqkVdRPSgqvT1CzAlGnw37GoEE6toaM20QMSUvJhB6E7XKrM6bUZbSaSImFniWitw+2JvHi8yYAK4fmxxBn5bo6uVw07CAIam6riGbw4PBZ2EFQIz6HW9IOCKURNTBeES4vwJmMCL75tuBN2lohWb22QR2cAgmCJ0mBnobvCzPLm4U4P6YUdJNLpu+vkMdmwUwCACuG5aKzGokR0OXhhzDghxubIBRkG1DoKld7akOvioMvB6YiTz0v2spvQtHrYDOM9RRm0mKiNCuEUKEDVe8wluUtgBwkNHIkiz4Oj0QdwGawNUvMwnlcCO0gI4OaX5g45dX162EEimtfa28AjVNkw19qegArhFDpsbRwvmSqtgB0kNOASpdQ6ikYfQDTmGTePW8VNTeiKcDownCfninXtx2EHiWjmphNjbEZWTAbsIACgQjil2oZv5GQ0xubADhIacKkqpbnNMm4bcY/CzhKh9FZjviCdxcJZCTCX8A8hqqy5+iHUNApTXftxOS+FQY9xGKgQTkHTq1FC3S45tLCS0oHLkReTWW9D0+rhMFgb8kk+Wllt+nIUlzq9rvp2I0mi0c5w6IdaVMkq2ClOQ4VwCrrx7uLsubBThA4MwyWKAhCDWkdhMVgbJbZxPB8VwmkZGxt76O//yBx2r7/9GkF0zG233Waz2WCHiiykY9TIGIW+1vYEVAjPZBsfGKPcufmLYQcJJbhUmWv3oPEyUFCAqrc1ZrV0oCvCabrmmmt2VjZQ3sy8dWscT7X8w0isXLnS60WzKYJnpEXbyWPIhPAXV/NBhfBMdaYfZC4OKyYBdpBQwpEos9t7662NFKBgZ4k4nfYePoOdgMcyYxNhZwkB1dXV/zluALe+3cmchzNaAD8ebHqhdpD56aefwo4WQfRNVVnsOC4Lhx3kNFQIz6TuPFEIe7vkkMNJkwpsNgGb12nvgZ0l4tRbGwpALBovOk11dXUgfwlgc4/hPx3GR9mUC2AMULiitrYWdrQIorfWKxIgb700GSqEZ9IOtxaloslYM8RgcLJlBRwR6iYMPr21QTrs5aIZhNOD4zhwjQIArMykWDe72HMEAACcIzhOl6uTsEcRXoPbUkKbDkKACuEZHF5nDzmmKFgKO0jo4UiVeeNMVAiDz2BtyOnq40iVsIOEhksvvZRpOgQGuwEAbI9YQh4HjmGg3r98+XLY0SKFu7vZJMBUdLreQIXwf2g6TkgcGD8FzZ2YMVyqyukbQDMogszpdXXau6QcITM6HnaW0JCTk/PgH34Ndq0EJ963jaVEETqwa+XNa1YsXEijC5Tw1tJYxWfiiTwavWNRIfwf6uYjco4IrbV9ETiZBeldfR32bqfXCTtLBDEOmLIZsVH5pbCDhJJHH3304NvPrzV/OF756Yhg/P2n73/11Vdhh4og6u46RQy99o5GhfB/aGwNRSI57BQhCWOx+anSbK6wcaAZdpYIYrA25I0DNHFipn7605/u27fvh4MnxjjMZfMVsONEFsN4V1F6OewU/wMVwh95SaLJO6iSorW2LxJHqswneGg2YTAZLEZJ7wAupcsKHaEFA1guM1bb8D3sIBHEa+tr5BHF2fB3pZ8MFcIfNZiNSS4qQVIEO0iowiVKiXUcdRMGk95skPGSGXwB7CChqjA+19Cvg50iglibTg6zsOxYWqy1PQEVwh/VNR2Sk9EYB42ivkicnMLszj50RRg0/WMW4PWkSejVyhRalDnzjY4+2CkiiKbtuJybzMDoVXrolQYubZ9WGZ8LO0UIY3D5qbFpXo/TMm6FnSUiGKwNuS42mko/G6qceU1cwm1FC0EEiXbQpEym3VQfVAhPowBlcPaitbZnCZeqCpjx6KIwOAxmg9Q2hkto97ESQmLxmBgG3mI8BDtIRCCd4w2MMfqstT0BFcLTOoe72V4iLZ9eXbghB5cqc+1kvRV1EwaDrlddwEvGcB7sIKFNHpWh6zgJO0VEGGvRtvGxwiTajcxHhfC0upYjMieHidbanh1coszusRjQeJnA85JEy2ivImsO7CAhrzCtuN7eDjtFRNA3V2WwY3ksLuwgZ0KF8DRN1ylFNL0GMoUihiA2jxlvspm8JAE7S5hrGmxJ9jJj89BImdkqyplvYrsI+wDsIOFPZzao6LTW9gRUCE/Tj7QXp6LlOfwgLkclZkS1DLXBDhLmdP166bCbk027VqaQk5sg7cfBcFMd7CDhjiQMbjPdZhD6oEIIAACDzuEhwpFbgDbj9QNcosxzs9Hq24FmaD9RwE1Gs31mj8Vg5uCJhpYq2EHCnKu72RSFqdLoeL2BCiEAAGg7a/LHMDwlG3aQcIBLVTn9w6ibMNAMQy3KVDRxwj8Kkwr1FvTVLbBam6q5TI6QT8fto1EhBAAAdWt1IZ6E1tr2C2ZCUq4bN/TrYQcJZyPuUTvhkMjQcoD+ocqaYwIjpGMUdpBwpu2qVURnw04xNVQIAQBAO9CgTEIL7/qNNL3I5hi0u0ZgBwlbhn6dZIziZslgBwkThSJ5UzTD3VoPO0g4TXnawAAAIABJREFU0452FmeUwU4xNVQIgYtwt3ntSuki2EHCB0+qzKX4RpsJdpCwVd9aLeOIMBYbdpAwIY4SMZjsrqbjsIOELWLQ3Mj1FucsgB1kaqgQAkO/PnOcislBV4R+w5EopYMugw11ugSE2WzWW42FqA3Dr2Sx2bpuNewUYcvSdHKYTbu1tif4sxCOjIxs37593bp1O3bsGBkZueDtNFHXfFgOotHyHH7EFmdIRkl9H+om9LOXX345JSWlYN4lJs/Aqy9+fPw4uoLxm8L0sga3mXK7YAcJT3Wt1TIe7dbanuDPWHv37hWLxXv37k1KSnr//fcveDtNaPr1CrTWtn9hmFIkr7c1UoCCHSV8/PWvf7390b/13f5F8s4vo0nG+8m3XLpilU6HthDyD6W4sCWW4+5AzRgBoR9sKhLTtw3Dn4Xw8OHDa9as4XA4a9asOXTo0AVvpwOSohqdNJ3jGdJEkhIegXXZ0aL+/uH1ep966ilw6zsgXVXk/SHaFeedd6Nj2dZdu3bBjhYmChLz2tiesSYt7CBhiHI5jIyRohzarbU9geXHY9lsNrFYDAAQi8UDAwMXvP0M69evP/vGBQsW3HrrrX4MeYZ2e2e0h8ST8gYHBwN3lukbHR0FAHi9XthBZosQZUpNVE1HrSCNP5vjkCQ5PDzM4XD8FSxEdXZ22txMkK4CAIgw/SiZCwAAhctr//1bmrx1IXK73aOjowzGbL/WiznxhuZqxeDP/ZIKrqGhIZIk2WxaDKdyNte18xhpeCqU92pMTAyTyTz/Y/xZCCmKwjDM9wNJkhe8/Qzr16/HzprJl5WVxeMFsPfO1KTOd7L44rTAnWJGfCUwoL9ykOTIc//jaR5qvix3xWwOQ5Ikl8sNhxdkduLj44FrFJBewGB52P1NnnUAADA+xOfz0YvDZDIJgpj96yAXFza1Ha7AOYBxgc9N+uPxeDwejyaF0NBzKp0ZHSeIg3L2s8vK2fxZCBMTE81mc3p6utVqFQqFF7z9DBs3bpxOYv/S9apV0VlcLl1WQ/d4PAAA+uSZDXls9usWwyx/F5IkcRwPjxdkNjIyMorleerj7ydXLBvke7RgEQAAHHlz1epV6MVhMBgej2f2r0NRWnF1XA3D2s3JzPdLMIh8fzU0KYR6i7E4o4DOb1R/9hHOnz//4MGDFEUdPHhw4cKFAAC1Wj3l7fShH+lQpdF0jmeok2WWdzgsLsINO0iYeOGFF6I+2lZheDTGiROG78Ge64q9zXfddRfsXOGjUFjQFIW5mlE3oV+RZD1d19qe4M9CuGnTppaWlg0bNrS1td14440AgG3btk15O01YHQNO0i3Jnwc7SHgSSIvTPczGgWbYQcLEwoULjdrasrRBXrd9hfaZXTcsOnbsWExMDOxc4SMrJmMI89ha0DYU/uTubWkUAFUGrfcL82fTqEAgeOKJJybfUllZOeXtNKHpPFkwSnFSJLCDhCdOtiz3325Dv14lQlsF+Ud6evpYPOOaS3+58r61qGvQ7xgYJkvIM7QasikKrTzsL6bGo9FMbiIvHnaQ86Hp9MbgULdWyzlJYNaDzZApYSx2Pkes7zoFO0j4IAbNTbhHnlEBO0jYUiQrmmM5nv4O2EHCh6b7lCI6C3aKC4joGqAbNClpPMczDChTiwyDrbBThI9Ow2GSxRLxRbCDhK1CYX5zPO5uRssU+I1hrLMondbtoiCSC6HD6+z02pW5aK3tAMrKm+sknDZHpE908xdd2zE5PxV2inAmTyxowsadLWi8jH8Qw7YGjoe2a21PiNxCqDcbssepqBwl7CDhDM9R5I6SejNadNQ/6gebFbTc4DtsxHNjo7iCtk4N7CBhot90YpzDyIqj6VrbEyK3EGqaqwooAVprO6AY3Kg8EKNrQ2tD+4HX0m3ikspM1EEYWAqRvAn3egf6YQcJB7WtR+V4CgboPvIocguhul9XlJAHO0X4UyTmG8xG2CnCwVhjbTuPKkhEC8QHVqFQ1iqOc6PZhP6gG2wuTlbBTnFhEVoISYpsdParstEMwoBT5S5scpsJioAdJOQ1Nh9N4cRHsWe1dityQYXCfBPudaHxMrNGuZ1GbKRYuhh2kAuL0ELYMtQW76FEeaiVKeDi8yoSXVTrYBvsICGOovS2hkIxmpEZcHnxki5ixN6Krghna7hV283DZEky2EEuLEILYV1LdYGDyYxDw9ADjiGIzfPg2pZq2EFCm6e/o4mPKdJKYAcJf2wmWxKf1UwOE/ZzbpWDTIe26YiUFc9m0mK90/OL0EKo7T6lEGTCThEp5DHZuq5a2ClCm8ukbolmFiaG/GLQIaFQWNCaJnK3otHOs6I161UJofGOjdBCqB/tKEpHw9CDRJFRZhzrgp0itFmbTg5i3qxYug9DDw+FiQXNMSzUTTgrFKV3m4sldJ9B6BOJhdAybnV53Vl5tF4NPZzkyy8ZIB2j7jHYQUIWRdX36QsSpAwsEv9gg69QWNBADrtaUCG8eM7elhY+UGXSfU0Zn0j8u9J01xaMUXiaFHaQSMFJTMlxMnStqJvwInm6W5piOcpQGIYeHlIEYgID5qFu0om+vV0kY8MhIYMfw4mGHWRaIrEQqlurZRxRGGxCHUJkvBQtKoQXy2mqa4nnFgoLYAeJIIXC/LbMVHerAXaQUKXuqVPF5MBOMV2RWAj1A02qJLSyWlDJxYV6S2NLSwtJkrCzhB6nSd2IjcrRSJkgKkwsaBUKUDfhRTOMddN8D8LJIq4QOr3OTu+wXBoaXbjh4bvvvvvbU28aiQHpghWiJPHu3bspioIdKnSQRFuPPpYbG8+NhR0lgsiF+Q0spxt1E14UYmTQiLvKckNgKr2PPzfmDQl6c332OCWQoCvCIKmpqVl15XrPL15ZQ76Z8dgnnf3gt3++wel0btu2DXa00ODuNLWIBIokNJU+qOSJ+c2OPkePh/K4MTYHdpwQ02Y8zGCyxQIx7CDTFXFXhOqWqgIyisGNgh0kUuzatcu96l6q+IpYd1whcRikysEtb+7cuZMg0KJr0+IyqZuF0fJE1EEYVHw2L1kg7klPdnc0wM4Seurajit5abBTzEDEFUJ9v74wHo0XDR6dTgcKLgEAjBFSIWYAAICMIquD7O9Hq/tPi6tJY2K7CoWogzDYFMKCtlQh6ia8CDp7S1FIDXKOrEJIUlS9s7c4G80gDJ6oqCgwNgAAMDHnE+weAADwuoFrlM9Hi0dfGEV4Rzrqez1D0viQGYAXNuSJ+SYeiboJZ4ryuI3YSGn+JbCDzEBkFcK24Y4YD5WE1toOolWrVoEfXgEA6JiLrDiR7WkHR96cX1EWFxcHO1oIcLcbW1OFufFSNiPiuvOhUwgLGt1Wd1s9IFEz/gxYmk8NcxgSYSjtFxZZhVDTdqxgnMlKCJku3DCwbdu2UqoN7F5HaL6McuBz9I/GHXj0+eefh50rNLga61qTE1C7KBRZsZk217AjUejuboadJZTUNR+Ws4ShtQpSKGWdPXXnSbkgHXaKyBITE3Ps2LH/u/mnV3a8we8ZLc4ab9Bry8rKYOcKDa4mtYnrRVPpoWBgmCwxry0r1Y26CWdCYzEqhSGw9dJkkVUIDSPtxaloI5tgY7PZW7du/fTTT7fc8IdBASFKTICdKDRQHre7s8no7FegQgiJQihriuOg8TIzQFEGr7VU+v/t3XdgG/XdP/DvaZ2GtyTLe9uSreGV4TiDvQpZLaWMhpI2bAql0JaShPIwSx9a+kChUPYqDS08P9KnQAjQgkecxLE1LVvy3pYlb1s6jbvfHwLXJduWfbq7z+uvIMu6d8TFH+s7Pt/1dOc4OxwqhOO+iemQL7ewmu4g3KVNL3PEIGhbdYb8XS1jmZl8Hl8pVdCdhaNKFOo2NE10WhG0gDgz04Pt/WKqhCG9tudxqBCaB41FsxSeAXsnaKOSKQUCvKelhu4gzEA4TT0ZyVolw0aZ2KREUdQ60YUksoCrj+4szGBp+yIPixXxGdaCgEOF0NR5qFigwPiw+o5OJYn51t5GulMwg89pdMbwYIKQRvF4XKI4YTQvF6YJz5Bp2KSLz6M7xVnjUCG0eRx6JfSpopkua1UrNRGaGqM7SLSjCG9gqNsR9MAEIb1KFOp2ZSxME54h29ygIXst3SnOGlcKoT/k7wqMlxTABCHNtMrijkQx0XqM7iDRjui0YlkFnZN9RUlM2o/FPjqlxin0Ex1muoMwQGBmvAMPlBdtojvIWeNKIWxzOzK8VFwuk7r+sJI6Kb+f759sPUp3kGhHOE292Zk58Zk406ZbWKZEoW6d6UMUFRpz0Z0l2tlaPlMgcSzOjMN4F+JKITR21mlIKU/KvP9DLCPii7LjMlr7mmAZ3qkRTpMzQahl2n4s9slPyBmZG/XnqQnotXY6pt5jOlkm3SkWgyuF0DJk0TJwCpeVdCpde6LY3+egO0j0Ir0zwdGBttA4LBmlHQ/jFSXm96QpiA4L3VminXW6uzS9nO4Ui8GJQkghyj43pM+GFqNRQatQd8hlPjusHT0pot0iyim2exw6WCkTBbRKTbsUg0+Ep0YG/W282XL1uXQHWQxOFML+qUFRMJRWCIdORAWtUtOGTftgvczJEe2mybwCClEqWTLdWQAqUahb/S5yZpKcmaA7S/TqaK3DMYEqgZE9LDlRCM09R9RenkCRSncQgBBCKbJknkA07O4m56bpzhKlCIexPUmqg90+0UGr0LS42wTZGqLTRneW6GXsrtfiKXSnWCROFEJLb2OxlEnHJbNesULdlZtBOIx0B4lG5OxUaMLdRk3CVvookSiOjxHJPDmZsJvwFCxjztLkErpTLBInCqF1slOfVkp3CvBvWoW6QyHztcI04QkQTpMoT2f1tOoU8IkwWmgVmvYEMRzSe1IUZSPHywrPpTvHIrG/EE77ZzykV12wju4g4N90ymI7b8Znb4RNFMcjnCZUoO2Z7C9KgnXO0aJEoXZQkwFXP+mboztLNBrstwV4KCdDS3eQRWJ/IbQNWXLnKHEWjDJFEXVSfu/siF8kDAz30J0l6vicxt7k+Nz4LMZ1LmYxnUJjcztEWUX+bjvdWaLLgQMH7rjjjt/+cW/6rNDj9tAdZ5HYXwiNHXUlvCTotR1VRHxRfkJOX2E+bKL4htDUGDk71YqmdLCDMJoUJOUOz7pCuWrovj2PJMkdO3ZcesOdz44VuBXS4R6sqFhbW1tLd67FYH8htLpbtXL4OBh1dEpNO0wTHodwGvECg83dBj1logof4xcm5nUqY2Fb/bx33333rc+b0ANH0fm3z8V7O7NvHf/272644QaSJOmOdtZYXghDVKg94NbnwwRh1NEqNG1oOtDroAgv3VmiCOE04YWlLR4HLBmNNjqlpk3o8w90UAE/3Vmiwv79+9F5tyChRBEanBYGW4RVaPVVHR5vSwvzTt5meSHs8HQqCCopn5Fdf9hNpyy2jjkEWUVEO/T1/zfCaR5Pz+BhPJVMSXcW8B+0So1trF2Ykg3dAcMmJydRfCpCqDL4aZIvNoiECCGUkDoxwby2AywvhMbOWnVAwpPF0R0EfJNckigTSj0FBdBiZl5o3EUS3lY0oYNx0eijUxS3uNuE+VqYJgxTq9WoowEhlE0dnQmpEULIO8UbbisqKqI52dljeSG0DppL4rLpTgFOTKco7lDEwHqZeT6nSVxU2uJxwGG8USgOj5VLEofTUmCaMOyWW26R1L8sNP+dwIcsgvPR3AR69YfXfXd7cjLz+gKyvBC2zPYbMivoTgFOTKtU2wNuKkAE3YN0Z4kKhNOIF5SaXXadkqkdOthNqyx2SEL+7lbEwPUgEVdUVPTx/31wleXRfhHZ+9Qv+bs1t63LeP755+nOtRhsLoRu75gv5M9Rr6c7CDgxnbLYOmoXayphdDSMcJpCeZqB6cFC2EoflXQKdctUNz9RGRjspDtLVNi0adMPb7o8V5j0xZ8eH+l2PPvss1KplO5Qi8HmQmjqPqyZ4wkV0GU0SuUn5IzOeQKFWhgdRQgFRwcQwhzUVGFSvpAH216jkVahsY624vk6GB39CkVZR21rcqs3bNggl8vpTrN4bC6Elu4jxZJUhGF0BwEnxsN4Gnlhh1zq77BQwQDdcWhGOE14UZll1K6HQyeiVU5C1iQxNZedB923w/zdLfYYXlnOGrqDLBWbC6F1ol2XoqM7BTgVnVLTMtUrSMn2d3H9gBvCaRIXGKxuO5y+FLUwhBUritrjRUSnFdrkIoSmzLXdYpIFzR9YWwiJkL83NF1SABOEUU2r0FjdME2IEEURHWZhvt7uhiWjUU2n0NhnB3i4NOjqpzsL/WyOmszYVKlQQneQpWJtIbQP2zK9ZEwW/HId1bQKTZunXaip8LUcpTsLnQIjvZhI3MPzyiVJ8Thse41eOmWx1d0qgmlChAKDXS3iYGk6G5bls7YQGjtqNVgiJoT+/VEtRiRTyZS9MkFoaiw0ydTW9UtHOE14QalltAUmCKNcsbzIOdbByyshOH82oddc51DKSpOZevTSQqwthLaRFm1SPt0pwOnplSVmT6tYXc7l0VHCacQLS62jrVAIo5xUKEmPTe1XxEEhnLXUO7EZAyv2vLKzEFKIshMjhlzotc0A+uRi66gd16wiOFsIKYrosOIFBuuoXQuFMOrplMUtQTcKBkPjLrqz0CY4NtLl9yhikuPwWLqzRAA7C2Hf5IAkGEorXE13EHB6BmWJyWUTa1b52pq42bAjMNDJi0kYF2HeoC8zLo3uOOA09Mpi66hdlKfl8odCr6nWmZ/JjnFRxNZCaO6oLSJwXkwC3UHA6alkyUKeYBjz8ROV/l4u9vX3OY3iwlKLq0WvLMYQbHuNdqXJWrOrBc/TcXk3oddc3xrLN0AhjGaW/ubimEy6U4AzpU8uMY+2iDWruHlOL9FuxgsNVrddC6fSM4FSqhDxhe60FM4WwtD0eMDV0+IdhE+EUc023WPIhDMIGUM/33SUg73WyJC/04YXlFpG7XpWrDvgAr2yxI6myOkxcmaS7iw08FkOjapLYkQyhSSJ7iyRwcJCOOOfdVPeoiLYSs8YemWJ2dUiytMGR/rI2Sm646wof5+Tn5Tsx0U9k/1qWOfMEHplscXdKsopITq52BHJa65ry1CUqdjTt4uFhdDaeyx3DolVOXQHAWcqLyFnzDc+GZzFC/S+tia646wowmnCC8taPc78hBwRH7a9MoM+ucQy2oLn6zi4Xob0zfq77S38mdJkKIRRzNxVpxWpoNc2g/AwTKvQWFxc3ERBOE3iQoPJ1WJIhnFRxsiNz57wTXkzczl4Wr3PdkRUoDe5W1kzQYhYWQitHie0LWacr6YJiyt9rce4086YCgWJbrsoT2cetbFmAR4X8DCsWFHUihMBVx9FeOmOs6K8lvpRdYlEIFbJlHRniRi2FUKSIp3BMW0hTBAyTHjhqECeiuGSwGAX3XFWiL+nVajKpMQSu9uhgyWjjKJXFls9DlFmIdHVQneWlUMF/ERbc1uSiE0ThIh9hbDd7UwiSEVOKd1BwNkplhd1TvT4goS4mEObKAiHES8sax/rSpElx4nY0KGDO76eJtRzanTU19YkzMg3jbeXsWiCELGvEJqcX6pRHCbC6Q4Czg7OF+Un5LR6HGLNKu40HSXaTXiBwQTjogxULC/qnOilctScOobCZ6mX6KvNLht8Ioxq1iGrNj6X7hRgMQzJJebRFrzA4O91cGHehQr4/X3teJ7WDCtlGCj8q1tHHN8/0EEFA3THWRFkyGs7PJqXy+fxVbJkutNEEtsKoX1uUJe1iu4UYDH0yhKLy46JcFG2xuc00R1n2fm7bMK0PISLLa4WdrTw5xpDcollrF2oyuJIa0CiwypIUll9Q2UqPd1ZIoxVhXDcOz5LEfmaDXQHAYthSC6xuVtDVEisqSQ4ME3oc5rERaW9k/1SoUQhldMdB5w1nbLYOmrH83V+boyOei31EsN644i1nF0ThIhlhdDUXlPkFwri4GcKI8WKYlJjVI6xDnHxKi70WgsfxmsehXFRptIri23uVgFHjqGgKK/lkERfbXRZS9k1QYjYVgh7G0sk6XSnAItnSNaaXDZhag4VCgZHB+iOs4wowhsY6hblFJtdNjZ16OCUeDwuWaroS4rxd9lZf4KYv9+JCUUjUgGGsLSYFLrjRBirCqFtolOfaqA7BVi8smSdacSKEGJ9A26i0yrKKsKEIugpw2ilKp15uoufqAwMdtKdZXl5zfUSw/rmEUs56z4OIjYVwkAo0ENOa9Wb6A4CFq9UpbWM2kmKZP0mivC46MjsaIAMZMTCYbxMFf7VjQtnE3rNdRLDepPLysoBDPYUwtZBUxqBxaYV0h0ELF4CHi+XJHZMdOPqcn+nlcWr0gmnSVxUah61wXpRRitN1ppHW1g/TRgcHaB8c6LMwqZhc0UKC0fd2FMITc5ajSAJem0zXWmyzjRi40liBKk5fpb+cCG9MwFXvzBLDTsImS5JkpiAxw2pEokOC4t75HpNtRLD+t7pQT5PwL4JQsSmQmgbtWvlRXSnAEtVqtIZXV9PE7J0dJRoN+N5WowvMLugpwzjGZK15tl+Hi5h8fIur7lOYqhuHjFXsG4HYVjECuH09PSePXu2b9++d+/e6enphV+68847L/ra73//+0hd8Rvsfpchv3qZXhysmNJkrdlloxDF4k0U4QnCcd+kxzuen5BDdxywJGUqndllE+Xp2NprLTTpCXqG8Xw9W8dFUQQL4b59+1Qq1b59+5KTk9999935xymK6u/v37dv3/79+/fv33/bbbdF6ooLDU70o1Aos2DNcrw4WEkKSVKsKKZ7sk+UWRSaHg9NuOlOFHmE04QXlhpHLIbkEh7GnlEZbipL1hlHrCL2HtLrtdSLtWspHs/kspXDJ8JTq62t3bp1q0gk2rp1a01NzfzjHo8nFArt3r37qquuevzxx2dnZyN1xYWMjn9pQlJMJF6OFwcrrDRZaxqxIgwTF5UTrDuwnpyZCE24RRkFzSMW9rWq4iClVCERil0qBVuPofCa6iSG9Z0TPTEimVKqoDvOshBE6oU8Ho9KpUIIqVSqsbGx+cfHxsYKCwtvueWW5OTkP/7xj88999zu3btP+AobN248/sHzzjvv7rvvPu3VjT3HCiXpC6/LUDMzMwihQIC1qyXPRGFM7tH+5k2KqlCmZspcP52UKxQK6Q4VMSFbA5ZZNDYxcWzIdE5y1SJu2snJSRzHvV729yU/Q36/f2ZmBqNvoVxxfOHhyZ5NfsLT5cDio6JUjI+Ph0Khpf/Doeam/b1tXmVObee/ShKKmPgzNj4+ns/nn/o5SyqEO3fu7O/vRwgdPHiQoqjwjUhRFLmgyUJRUdGTTz4Z/vOuXbt27dp1slfbu3fv8beySqWSyWSnTeLwDl6Us+VMnhnlwm8dC/4iS7E6o+L1tndlMhmpXzf2yVtSMc6mN2R6wCFRlxOC0JR/uiSlmHf2P74DgQCO4xKJZDniMZFQKKQoisabpCKt9NiI6cJcrWCkS5yWTVeMhXw+n0wmW3ohJFoPk4VlsoREe4vjvMwNTPyXyOOdfuBzSYXw1Vdfnf+zXC53uVwZGRlut1uh+PfvRA6HIxAIaLVahJBQKDzF/5iLL754cb/TeQPeIcyr116A44w/htDv9yOEWPAXWYoMPA3ni0YId5YiXSBXUe5+PDuH7lARM95hid+4pWm81ZBcIhEvZjAf/1rEszEUhmF+v5/GN2RVWtlr1nckhdsDfQ686hK6YiwkEolwHF96IZxpOSwr3yQUCa3u1l+su4utd13E5girqqoOHDhAUdSBAweqq6sRQiaTCSHk8/kefPDBnp6eQCDw1ltvrV+/PlJXnGfrqM/2C8RyFu5u4axylb55xIwQwtWVVBd7pl5Ckx5yblqYlmscsZar2LkAj4NSY1R8Ht+dqmLZMRSUnyDaLeLiNc6xToVUniiOpzvRcolYIdyxY0dnZ+c111zT3d193XXXIYTuvfdehJBer7/uuuv27t179dVXT09Pn2JodNGMXQ0lOFRBVilX6ZuGw4Wwguww0x0nYsIbJxCGNY+Yy1NgpQx7lCZrbdhUaGqMnJmkO0vE+OxHRTkanjSmacRcwerf2yK2WCYmJubRRx9d+MjBgwcRQhiGbdu2bdu2bZG60PHs4+2XqiqW7/XByqtMKX226RWSokS5JWh8hJyZ4MUk0B0qAsIbJzze8QliKjc+KiaTQESUJuuMLtva7GKiyybRs2RDs9dcH/67NI2YtxRcSnecZcT4PUwUolqD44YC6LXNKgqpPA6P7ZzoxvgCLEvjcxjpThQZhNOIFxqaR8xlybpFLJMBUasixdA0bMLzdazZREGFgr7WRrF+XYgK2UZby9h46MQ8xhfCnpE2WYhSZcMoE9tUqAxNwyaEEJarY0evtdCYiwoGhMmZxhEr7CBkmbSYFD5P4EpTsWZbPeE0CZIz+HFJLe62jNi0WFEM3YmWEeMLYXPbFxosHp3BAlnALBUphuYRC0IIyzcQrY0s6GjsczbjheEJQgtbezZyWWWKwcybCYz0UQQbtnh6LV+Nix4dMlamlNIdZ3kxvn7Yhq3axHy6U4DIq1AZTC5bkAxh8QpMLGPBwafhCUK3d2zGP5uTkEV3HBBhFSmGJpdFlFlAdNvpzrJkFOWzHJLo1yGEGoeMq1LL6A60vBhfCFt8Q/psaDHKQnF4bGpMStuYEyHEjgbc4SWjTcPmMpUOQzBByDarUsqaRyzCPC0Ljg/zd9t5sjiBMn02MNc12aNTFtOdaHkxuxBOeSfGkV+tgZUy7FT59eioWFPpa2V2IQyODiAeX6BIbR4xs3vdAWclihMUkqS+VAULjqHwWuolhmqEUPOIRavQ4HwR3YmWF7MLobn1n/lBMV/MvK4/4EyUqwxN4W31BaWB/g7Stywd21cG4TCKC0sRQk3D5soUlg80cVZ5isGDl349AAAfsElEQVTCn/X3tVNBZrcL9loOSQzrETfGRRHTC6Glt6lElk53CrBcSpO1rR6nnwxgQpEop5hwmuhOtHi+dhNeWNo7NRCiyKw4uGnZqUJlaHbbharMQJ+D7iyLFxjqRqGgMD0fIdQ4bFzFgd/bmF0IbVOdujSWL2fiMqlQkpeQ45joQAjhmkqilbFHMlGUv92C5xsah4yrOfD7NWeVq/RWt52XW0IweTeh11wX/jg4Mjs6TczkceDsaAYXQpIiO9CMXnMO3UHAMqpQGazjrQghsWaVz36U7jiLFBjuwXAJPyn52LBxVWo53XHAcokRybLiMjpTExm9m9BrrheHx0WHjatSy7jQ+YHBhdDRdUQe5CUqc+kOApZRRYrBOtaKEBKmZFEUFXT1051oMQinCS8qC1Eho8vK7p6NoDKl1Cr0+jtb0ILT6BgkODZCTo3hOcUoPEHI9h2EYQwuhOb2erVQTncKsLy0cnX/7OCMfxaF144ycxMF4TThBQa725kiU7G4hT9ACFWmlDZ5WvkJ8sBgF91ZFsNnrhNr1yIej6SopmFzJTdG8hlcCK3uVr1cTXcKsLyEfKE6viC8dlRcvIqRvdYoiuiw4AWljcPNMEHIenplccdEVyivhKGjo/PrRdvHO+Pw2GSp4rTfwgIMLoStgVF9XuRPNwTRpkyhPTx4DCGEF5X7u2xUwE93orPj72/nxyXx4xIbh0ysb1UFRHyROqnAmRLLxPUy5MxEYKATLyxFX08Q0p1ohTC1ELrH+mYRmVcIPWXYr1SuOzrUjBDiiaXCtDzG7VYOj4vOBbwdE12GZC3dccCyW51abuTP+jutjGuQ67U2iItXYUIRQohTK5yZWgjNrf9UIxnGi9h5iiBqpUlVGMJ6JvsQQnhxJeNGR4l2E15UZnRZiuVFrO/QARBCa9Mqjo7ZMaEo6B6kO8vZ8ZrrxIZqhJAv6Gv1OMuSudIanrGFcMCkjYX1olyxJq3iyFAzQkisqSSY1WuNDPk7W/B8/dEh42rYOMEN+Ym53oB3LL+AWaOjFOH1d7aIi1cjhJpGzBp5oVQooTvUCmFYISRJ8r333rv//vuNEx2xFCdmcQFCaE1q+ZGhJoSQKKOQnJ0OjbvoTnSm/L0OviKFJ4vjSKsqgBDCELY6tdyqEDNrvYzXdliUr+OJpQihhoFja9Mq6U60cphUCIeGhlavXn3lfb/7rUM4iIduvvcPu3btCoVCdOcCy64ipdQ22uoP+RGG4epyBo2OEk6TuKB0dM49SUzlJ8AYBlesSatoQhN+Rs1nzzfaRggdHoRCGK1uvfXWpsRq9It/GS6oSvTjEz9revmTxhdffJHuXGDZyYTS/MRck8uGwi1mmFMIfU4TXlh6eLCJIx06QNiq1DLzZIff7wtNjNKd5YxQwQDR1iTRViGEwvPxOfGZdIdaOYwphLOzs//48CO07b8QQtnUEaFfhURSdMX9+/btozsaWAmrU8uOzk8TOoxUKEh3otOjQkF/T6soT9sw2LgufRXdccDKiRPFZsdnduZlEp02urOcEaKtSZiWy4uJRwgdGmzk1MdBxKBC6PF4gpJEJI5FCPWhyhZsK0IIKXJdLsZMF4GlWJtaeWSwCSHEi4kXKFL9Pa10Jzo9f7ddqMoKifCmYfOa1Aq644AVVZVWaZaL/AxZL+O1HJLovxoXPTLYVJUOhTAqqVQqaXAaTQ4jhI6ILmkWnocQQv2WnJwceoOBlVGYlD9BTLrm3AghsWYVYWfA6Gi4xajJZctNyIrH4+iOA1bU6tSKJtLNjPUyJOm1Noh16xBC3qCv1eMs51hHXMYUQhzHr7/+evTW7cg/99VDnh70wa9uvPFGWnOBFcLDsIqU0sYhI2LObsLwVvqGwcaqNBgX5RyNvHA8OOeadZOzU3RnOQ2i08ZPUAjkKQihxiFjiUItEYjpDrWiGFMIEUK//e1vd+gTsL069NL16NnvxDxR/fsH7t22bRvducAKqUqrbBhsRAjhOcVBzyA5M0F3olOh/IS/vx3P0zYMHquCCULu4WHYqtSylpzk6J8m9Frq58dFubZeNIxJhVAqlb7xxhv2+k//fNeW/Q/d2Nlmu+uuu+gOBVZOVfqqY8Mmf8iPeHy8oNQX3ef0+rtbROl5Q/6JucBcQSJsnOCiNanl5jieP+pHR30LNk4cGWqqgkIY/dRq9TXXXLN582alUkl3FrCi4kSx+Qk5zSMWFD6JIrqPZPI5jHhhaV3/kXXpqzEEGye4aE1ahTnkno3u9TKB/g7EFwhTcxBCXRM9fIyfGZdOd6iVxrxCCLhsfcbauv4jCCFx8Wpf27FobmpMtJvwwjKYIOSyRHFCRnyGZbaHIrx0Zzkpr7l2/uNgw+Axrq0XDYNCCJhkQ8bauv7DFKL4CQqeLN7f3053ohOjCG9gqIdMz7G7HXD0EpdtyKxqTouL5t0+XnO9RP/VeXa1/Q3r09fSm4cWUAgBk6THpkqFUsdYB/rqnN4oHR0lOiyiLHWju6VYUcSdzsXgeBsy1h6REr5o7bUWdA+S3hlRthohNO6b6J0cKFXp6A5FAyiEgGHWZ6yp7z+Kwi1monU3IeE04YWlDYON69JW050F0Ck7PlMqkrb2ROmN6jXVSvTrEIYhhGr6GtamVQo5ebYdFELAMOsz1tT1H0YI4fn6wGAn6ZulO9EJ+JwmUYGhYaCRmzMuYKEN2esbvN3R2RTQa66XGL4aF63rP7I+g6NHnUMhBAyjVRS7vWMjsy5MKBLllhAOI92Jvon0zoTcQ04ZisVjM2LT6I4DaLYhe8OxREGgz0l3kG8KTXqCo/2ifD1CyBv0WUZbOLiDMAwKIWAYHoZVpVUeGmhECIk10dhihnCaRbnFXw4cPiezmu4sgH4lCvWUiNfjqKc7yDd5LYfE2rUYX4AQOjRw1JBcwtn5bCiEgHmqM9bUDRxB4SOZWo7QHeebiHYTXlBa29dwThYUQoB4GFZIKt9v/sfdd9/95ptvBgIBuhN9xWuuW7Be9PCGDC6uFw2DQgiYZ01quW20dS7gFagyMb4wMNJLd6L/QDhNPWlyHsbLS8imOwugmd/v/853vvPXtxrb44LPDKVf/8Sb5eXl/f39dOdCpHcm0OsQayoQQkEydHSwuTqdoxOECAohYCKxQGxILjkU7juqLieiqcUMOTMRmnAfIvrh4yBACD3zzDPvm0f6tn/SK8OKN16AfvIPW8Ylt99+O925kM96GC8qxURihFDziDkrPj1Jkkh3KNpAIQSMdG7Whn/21KKvdhNG0TQh4TTj+bov+w5BIQQIoffeew9d9jOSL0v0KtXUQYQQuvz+Dz/62OuludeM11Iv0a0L/7mmr2FDRhW9eegFhRAw0sbMqqZh02xgDi8s83fbqYCf7kRf8bWbhnOyiFCgKCmf7iyAfh6PB8mzEUKj5Bo570imfwRJ4oKimMnJSRpTUX6CcBjFuiqEEIWouoEjXJ4gRFAIAUPJhNIyla6u/whPLBVm5BPtZroTfYVwGBukvnMy10GjbYAQysvLQz1NCKEvhFe3xVJv9/5MM3wojh9UKBQ0pvK1HhNmFfGksQghu9sRK5RxsNH2QlAIAVOdn73xnz01KLx2NDpGR0OTHnJuumbcfk7WerqzgKhw6623ov/di0acISRxUWseSyt/Z/CxB278vkBAZwMXr6Vufh/9we4vzsveSGOYaACFEDDV+oy1Jpdtyj8t1qzy2Y/SHQchhAinaaxQPUlMlSjUdGcBUWHLli3P//qBhN+fh353qWt/zfDYh18I0q7G3QR93UepUNDXclSiX4cQIinqy97687M30BUmSkAhBEwlEYhXpZbV9R0WpudRvtng2AjdiRDhNB1V4psy1/EwGBcFX7n55pu72x2fPnP/KztvTitJvfbXv1Hc8EvPq4/6Wuj57c3fYREo0vjxcoRQ84hZLk3i+LgogkIIGO387I2f99YiDMPVlUQUjI4STuO/gv0X5myiOwiILvHx8RdccMGWyzefm73+s54avLBMcfND4395aq75i5UP4zX/e1z0s+4vL8yG2xUKIWCyqrRVLe62CWJSrKmk/cD6oGe4Q+AjUEinLKY3CYhaF+ac81n3lwghUWaR4tbHJv/fn2YPfbSiCSjKazkUPok3QAZr+w+fz/kJQgSFEDCaWICvTa2o7WsQayqJdjO9Df6JdlNtVsJleRfAelFwMqXJunHfZPdkH0JImJqj/PGT05/um/78bysWwN/bxpPIBMp0hNDhwWO5CdkKqXzFrh61oBACZjsve8PnPbU8WZxAmebvstOYZNZhrBdMXJx7Ho0ZQJTjYdj52Rs+76kJ/6dAkaq863dzRz+d/PsrKxNg4bjo591fXgAfBxFCUAgB061Nq2wf7xqZHaX9wPrDw03ZcZmpMSoaM4Dod0HOpk+7/z01yI9LUv74v4l288Rf/4Aoarmv7rUcChdCX9B3eKgJ+h+FQSEEzCbiiy7I2fRR56diTSUthZCiqO7ubuPnH/8rLnip+tKVDwCYRZ1UwMd4LW7H/CM8aazi1scCwz1jf34SkaHlu3RguIcK+oXp+Qihuv4jOkVxPB63fJdjECiEgPGuKLj4w47PBFnq0JgrNDW+kpduaGgoLy/Prdz05G8eM+HER8/tn5ubW8kAgIkuy7/w7+0fL3yEJ5YqbnmUnJ3yvPooFVyuc5q85nqJvhphGELo0+4vLoDlzV+DQggYLz8hJ1Ec3+gy44VlRNvKbaLo6uq66JLLTGvuRU90JF9SyicK//Bp580337xiAQBDXZZ3YU1fw4x/duGDmFAk/9GvkEDgfn43RSxLS26vuTY8Lur2jllHWzdmcrrR9kJQCAEbfCv/wg87Pl3hkyheeOGFmVXXoTXfwyjKFdPXgW1Gu17783sfDA8Pr1gGwESJ4vg1qeUHuv75jccxvkC+4z6BIm30j/eTczORvWhozBWacOO5WoTQhx0Hz8veIBGII3sJ5oJCCNjgwpxzGoeMvrwiX1vTCqw4CLPb7Vhhdbm37Z6x308KsU7eeiSJJ1OK29raViYAYK4tRZftd55oByGPl/i9u/DcktE//Cw0Hclxfq+lXqKtQjweSVH/aD94RcHFEXxxpoNCCNhAJpRuyFj7mcfEj0309zmX+3JUKOhra7o+lXdU9O5/DzzTHd8/GLyYCv9rmhqJi4MFCOA0ypJ1GMJMLtsJvoZh8VtvlJRtGn363tCYK1JX9Frqwvvojw41JYjj4ZiwhaAQApa4vOCiv7d/gmsqlq/FDOX3eU01Y2/+ZmjvNVMfvakq0n73k55Lcx4cjPFYsSsRQsj0j0wpaTAYlikAYJMrCi7Z7/z4ZF+Nu/iamI1bXM/cG3T1L/1a5MxkYKALLypHCP29/cBm+Dj4n+g8CgSACNIrSzAM68pUpdd8GXfJtRF8ZXJ2yms95DUf8ndYRDkaiWF9/NZd/LgkBUmuPdoRqrvKrTP4++uQszb26BuvvvsOn8+P4NUBW12ad/5rlnfGfZOJ4vgTPiFm01aeWDr67C8UNz0U3vOwaF5bA66pwISiMe+4ccR6/7q7l/Jq7AOFELDHFfkXfTjm3DnUTc5Nhw8dXYrQuMtnb/TaDvs7LKJ8vbRsY9L37+VJYuafwOPx3nz7ra1/+b786Jy65xV9pf72V4xpaWlLvC7giBiRbGNm1Uedn15b8p2TPUe65iKeLNb9wh75zj2iXO2ir+U11UlXX4AQ+kfHwfOyN0iFkkW/FCtBIQTscXnBxVd/8Ndv5xckOo2S0kX2jgp6hnzWw3PGmqCrT1yyWrb6QvkPfomJTry+7oveuiJV/lNPPbKE1IC7thZe9mDtb64u3s7DTjpLJdZWJe2QeF55JPHae8TFqxZxFYrw+jttSdffR1LU/7V/8vCmXy4hMjtBIQTsIRNKNxdc8mGP/Yf2xrMrhBTl727xWg55zXWIosT66vjNO/GcEsQ7zST6X1v3X6+/akmhAYdp5IXJUsXnPTUX5pxziqfhhaXyH+71vPJwwpW3S0rP+hBdX8tRUW4JTyw9MtgUL46DZTLHg0IIWOU76s03OD7a7OxNPIMnU8EA4TR6LYd81kO82ESJvlq+c8+ZT8ZYR1un/TNVaauXEhhw3PW67z1z7MXzszed+jBnUW6J4pZH3H96gCS8sjUXndUl5teLvtv6wfaiy5cUl6WgEAJWkUsSVysr/l9Pbff3toUSVFu2bLn88m/+yyd9cz77UZ+53td6TJiaIzZUK+/6nUCeerbX+mvrB99Rb4bD6MFSrEotkwolNX2HTtv/Wpier7zjN6N//CXlnY05Z9sZvj4VDPhaj8Vvv6XV4+ydGrgo59ylJmYj2D4BWOWTTz55+tanP1WgYIr+T3OGK268d8eOHRRFIYRC0+Oz9R+6X9g7/OCOucbPcXVFyu6XlXc+GXvutxdRBTsmus2jLZflX7AMfwnALTt033vD+i6FTt8IQqBMT77ztzN1/zf18Vtn+OKE0yhMyeHHJr5l++s1JdsFPFjSfALwiRCwRzAY/OGPfuTZ/kyArFNkuZD+AbTu+zW/33jo6QcL0ExwpEdcvFq29mL5Dfdj+FJXzb3Q/Pr3td+FJlVg6arTV79mfufQQGN1+umH2fkJyuQ7nxx9fjfpnU3YdhM63YCE11wvNlT3TPa1uNv2rr83QpHZBj4RAvZobm4e8IuR4fIW6vqWpKF7XG8c7Pn5+5fkjrRZ4i65NvXhvyTt+IWkbOPSq6DJZeud6t9ScElEYgNwnfbKNyz7zvDJvJgE5e1P+Hsd4395CpHkqZ5Kkj5rg0Rf/ab1r9/VbMX5oghkZSMohIA9pqamUFwyQsiFNMJgoiem8770O1b7rnxrBBNrKjF+ZMY/KET9senVG0t3CPnCiLwgABsz1/mCvqNDzWf4fJ4kRnnro6FJj+f1x6lQ8GRPI7pbeHFJLhwdHWreWnhZhMKyEBRCwB4ajQYbtCHfNELoQ96DbUm9Zlkq2XGouLg4glf5orc+RIXOz1nkPkUAjsfDsB/or/6T8Q2SOuUnvAUwkVh+438hhNzP7znZsU1ec73EsP7Ptve2FX0LNtGfAhRCwB7p6elXb/0WeuWHaG58jsoYCZ2TO/FIXPO+m266KVKXCJKhl0xv3VJ+A4ZgsSiIpPOyN8SIZP/r+PDMvwXjC+Q/+KVArhp99j5ybvr4J/gsh6aKNF/01X9Hc0XkkrIQFELAKi+++OId6zMFuzXosQ09ux9TxThf2/9qVlZWpF5/v/NjlUxZmVIaqRcEYN7dq295w7JvzHs2py/xeInf+wmepx197j5yZmLhV4IDnQjDnu358HuabXGipXYcZDcohIBVZDLZM8884xnoPvyXZ9prP7//ons+mas5k4XpZ6J/evB1y19+XHljRF4NgG/Iisu4ouCiZ5teObtvw7D4bTdJy89xPfOz0MQoQmhwcLCmpqbvnx8Yi3MHZoa+V7J9WeKyCBRCwEJxcXFr1qzJycm5JO88AU/wcefnS3/NEBV6tP6pH+ivzonPXPqrAXBC1+uvtrlbm4bNZ/uNsRdcFbNh88jvf3r3DddmaEq//cv/6aw98BvXoQvF64U82CZ3GlAIAZthCLuj8kd/Mr4xMrvUA07fsLwbI5JtV38rIsEAOCGcL/px5Y1PHX0+QJ50LejJxGzc8krb2HXS8eIHP825843DuoQhQn3bt2/r7u5ehqSsAoUQsJw6qeC6kit3f/GYL0gs+kVa3G372z++r+ouWCMDltv6jDW5CVnPNL54tt/Y2dm59+8N/5Vz99uDv/7RxEu1cqwjdvfsqmtffvnl5cjJJlAIAftdqdmcn5jzm4anFzdZOBuYe6T+dz9dfYtcciatvAFYqvuq7jK7bB84Pzqr73I4HChD/39J5/008zanwhiYuzBAxaPcNQ6HY5lysgYUQsAJ96y5bWBm+C8t/3u23zgbmLv381+tT1+zMXPdcgQD4HhSoeSxc/e8ZvlL84jlzL8rMTERTQxiiBxJ/LJDUFrPvx4hhCYGEhISlisoW0AhBJwg4ose3XT/39r+XtPXcObfNeOf/elne0sURbdV/nD5sgFwvLSYlF+tv/ehuicHZ4bP8FsqKyszcSJ/8gE+8taT9wcwIfJOopqXt20706MqOAsKIeAKhVT+2Kbd/9P4p9ct+85kjHTKP/3Tz/YakrU/rrwRpgbByitT6W/QX/2Lfz40MD10Js8XCAQ3/uH2JPJYy9tu6vBf0SdPYQ+v+fHVV1x2GTRXO40IF8JQKLRz586Fj0xPT+/Zs2f79u179+6dnj5B7wMAVoxaXvDCpb89OtR8/78emfHPnuKZTcPmH39yX0WK4fYK+CwIaLO18LKrirfe/snPa/oOnfqZJEW+0/K+NeR489rnHz43+8qJv9+dMfzFe68//fTTKxOV0bDwUW0R8f7773/++edtbW0HDx6cf/Cll17yer0333zzCy+8IJVKf/SjH504B4aRJIlx/ozT8O8KsbHQBuIrJEmOjo6qVKoIvmaQDD3X9ErDYOO1Jd/elFkdh//Hu9012ft802t90wM3lV1/btb6CF43UiYmJnAcl0igdeRX/H7/9PS0XC6nO8hyafO0/6r2iXMyq28qv56PneBAwfbxzt8c/kOMUPaztXekxqgQQm63Oz4+XiiEvvBnJJKF0Gg0+ny+vXv3LiyEN9xww0MPPZSVldXb2/vAAw+89tprJ84BhRAhBIXwOMtRCMOODDZ91PnZ4cFjWoWmPEU/6Ztye8dGZl3900Pf1353a9FlUbsNGQrhN7C+ECKEpvzTj9Y91TPVtylz3cbMdVqFhodhHu9431R/3cDRT7u/uLnsB5fknTc/hg+F8KxEshCGXXTRRQsL4ebNm//2t7/hOE4QxHe/+939+/efOAeGqdXq4wvhJZdcsmfPnsgmjGYzMzMIoZiYGLqDRAuSJD0ej1KpXKbXJ0L+Jre5bbw9AY9PxBOSJfLs2EypIKprzOTkJI7jYjGcCfwVv98/MzOTlJREd5Bl1z3dd3ik6YiryUOMh8iQABOkyVLy47OvzN8cJ/yP357HxsZiY2OhECKEEhISBILT/FK71N95d+7c2d/fjxBaWPwWoigqXN4oiiJPeYbkO++8c3whTEhIiI+PX2JIBuHxeAg+ES5AkqTf71/We+DSpAsuRRcs3+tHHEVR8IlwIb/fz+PxuPCDojQ+vjRDdxO63u314Hw8VnTS35gDgQB8Igzj808wmPwNSy2Er7766qmfIJfLXS5XRkaG2+1WKBSneGZZWRkMjYZvXLh955EkKRQK4Q1ZSPg1uoNEC4qiuPaGpApTTv0EuEnOyjJunzCZTAihqqqqAwcOUBR14MCB6urq5bscAAAAsAjLWAjvvfdehNCOHTs6Ozuvueaa7u7u6667bvkuBwAAACxC5BfLLA6sGg0LLybasmUL3UGixczMzMMPP/zEE0/QHSSKvPzyywaDYfXq1XQHiRbt7e0ffPDBPffcQ3eQKPLrX//62muvjeCR1OwGnWWiS0NDQ0PDWfQAY725ubkXXzzrNvzs9uGHH9rtdrpTRJH+/v53332X7hTR5e233x4ZGaE7BWNAIQQAAMBpUAgBAABwGhRCAAAAnBYti2VSUlIuvvhiulPQz2azIYS0Wi3dQaIFQRAHDx684oor6A4SRQ4fPpyamgrrIOa53W6r1XruuefSHSSKfPbZZ5WVlXASIULokUceOe0/lmgphGaz2Wg00p0CAAAAq2zevDkxMfHUz4mWQggAAADQAuYIAQAAcBoUQgAAAJwGhRAAAACnRenRo9xx5513zncJufzyy3/yk5+E/zw9Pf3EE0/YbDadTvfzn/+cUwcz1dTUvP766263Ozc395577snIyJj/0sneLrY62W0AtwfcHvPgZ0gEUIA+JElu377d4/HMzc3Nzc0RBDH/pRdffPHpp58mCOLpp59+6aWXaAy5wgYHBzdv3my3230+3759++688875L53i7WKrk90GcHvA7REGP0MiAoZG6eTxeEKh0O7du6+66qrHH398dnZ2/ku1tbVbt24ViURbt26tqamhMeQKGxoaOv/88zUaDY7jF198cV9f3/yXTvF2sdXJbgO4PeD2CIOfIREBhZBOY2NjhYWF99xzz9tvvy2TyZ577rn5L3k8HpVKhRBSqVRjY2P0ZVxpFRUV4bGdUCj0+uuvL9wlfYq3i61OdhvA7QG3Rxj8DIkImCNcaTt37uzv70cIHTx4sKio6Mknnww/vmvXrl27ds0/jaKo8KFUFEWRJElL1BWz8D0JP9LY2PjSSy+tWrVq586d8087xdvFVie7DTh1exwPbo958DMkIqAQrrRXX311/s8OhyMQCIQbqgmFQqFQOP8luVzucrkyMjLcbrdCoaAh6Apa+J5QFPXiiy/a7fY9e/YsXAeBTvl2sdXJbgNO3R4Lwe3xDfAzJCJgaJROPp/vwQcf7OnpCQQCb7311vr16xFCJpMJIVRVVXXgwAGKog4cOFBdXU130pVjNpsPHTr08MMPy+Vyr9fr9XrR1+/JCd8udjv+NoDbA26PheBnSERAizU6URT1wQcfvP/++7Ozs2vWrLnjjjtkMtlFF1108ODBmZmZxx9/vKOjo7Cw8L777pPJZHSHXSFvvPHGm2++ufCRgwcPht+TE75ddOVcGcffBnB7wO2xEPwMiQgohAAAADgNhkYBAABwGhRCAAAAnAaFEAAAAKdBIQQAAMBpUAgBAABwGhRCAAAAnAaFEAAAAKdBIQQAAMBpUAgBAABwGhRCAAAAnAaFEAAAAKdBIQQAAMBpUAgBAABwGhRCAAAAnPb/AYtuaUpaMrNiAAAAAElFTkSuQmCC" }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "li = LinearInterpolation(x, y)\n", "li_spline = CubicSplineInterpolation(x, y)\n", "\n", "@show li(0.3) # evaluate at a single point\n", "\n", "scatter(x, y, label = \"sampled data\", markersize = 4)\n", "plot!(xf, li.(xf), label = \"linear\")\n", "plot!(xf, li_spline.(xf), label = \"spline\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Univariate with Irregular Grid\n", "\n", "In the above, the `LinearInterpolation` function uses a specialized function\n", "for regular grids since `x` is a `Range` type.\n", "\n", "For an arbitrary, irregular grid" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "hide-output": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAGQCAIAAADZR5NjAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de1xUdf7H8e85M4OiCCgioGjmLSvvmRpeUstLKQlSltdSyW677vbLX+1m2nVrdc3dtjKzrS0va+jPa1mZpQmoeQ/U1DBBRVFCUEcZhjnnfH9/jEvmBUUGzlxezz/2gWcO53xs5wFv3+fM9yhSSgEAAADPUc0eAAAAwN8QsAAAADyMgAUAAOBhBCwAAAAPI2ABAAB4GAELAADAwwhYAAAAHkbAAgAA8DACFgAAgIcRsAAAADysugOWlDIzM7MyR9A0jcf74EIul8vsEeBFpJSappk9BbyIYRi6rps9BbyIruuGYVT1Wao7YJWWlt5+++2VOcKpU6f46YkLFRQUmD0CvIiu60VFRWZPAS/idDrtdrvZU8CLFBcXFxcXV/VZuEQIAADgYQQsAAAADyNgAQAAeBgBCwAAwMMIWAAAAB5GwAIAAAEhMzMzPj6+adOmzZo1S0pK2r9/f9Wdi4AFAAD8X3p6eu+4bnHFB9fd3/7rhDZtCn6M63L7Dz/8UEWns1bRcQEAALzH888///wdzR66Ndb9x0c73iClnDp16sqVK6vidDRYAADAz0kpt2zePLB51IUb72kRtWnTpio6Iw0WAADwMU6ns0aNGteypyFFZqFcn2dolzwdp0ofvEeDBQAAfIPT6fzLX/7SuHHj0Nq1YmJinnvuubNnz166m8sQG0/IaRnG4NVaxDzXyHX6vtOiZcduX/184sLdvvr5RFxcXBWNSoMFAAB8w4gRI05uT53ft3Xzurfm2h2vL5t7z8aN69evV1W1WBObf5GpeTL1uLH1F9kqTOkZrYy/Sf33nWpkTSGE2PiPN+IH9HMZclCLKEOKZfvz3s04ti5tURWNSsACAAA+YNu2bd+v+WLd6B61bBYhRGyd4H8ObHfvwk0PvvX5sRb3ZhbK9vWUXtHKpLaW7tFKqO3ib4+Li/tu4/dTpkx5c0m6xWLp1avXpq0rWrVqVUXTErAAAIAP2Lp1a1xsPXe6crOqyp031N+3b9vrD8V3iVSCrxZq2rZtu3z5crvdrihKSEhIlU5LwAIAAF4t3yFSjxuLD6kh+sV3qpfqRt9Y250xiimDlYOABQAAqoNhGJmZmYcPH27WrNmtt96qKOWlouMOkZpnrD8uvzsmjxXLntFql7ieH75dVFBcWr9WkHufcy79m+z8eXfeWS3jVwwBCwAAVLmMjIxx48adydl/Y3jtn06ebdLuto8++qhFixYX7nOsWK7Pk+vz5PrjMt8he0SrfWKU8a3U9hGKRRFC3FLy6OMjF3z0p+6tbq4fcvBU8YxNWXcMvK9Xr14m/Z3Ko8gqXQXiEk6nMzQ01Ol0XvcRCgoKwsLCbLZL7l5DoMrLy4uJiTF7CngLTdOKiooiIyPNHgTewuFwOJ3O8PBwswcJaGfOnLnpppv+cEu9UW0bCyEMKd/dmr2sQOzevTvfZfsuT6Yel9/lySKn7Bmt9o5R7oxR2tZV1Ms1XIsWLZo9e3ZWVlbTpk0ffvjhcePGqWrF1pziHiwAAOAPli5denNNzZ2uhBCqovy+S7MvF25q+uJqvc3AXtFqr2hl4q3qrXXLvWoohBBi2LBhw4YNq+qBK4+ABQAAqtaBAwfaNQi7aGP7qLDG4Qf/MtLmdTeoewIBCwAAVIkTDrHumLEuTy7Lrdv7bMlFrx6zlwyMre+X6UoQsAAAgAeddIr1eca6Y3JdnjxWLHtFq31ilPt+P2RUrxd+LjrXvG5t924/HD+9o7Dk0/79zZ226hCwAADAZUgpP//88+3btwcHB/fu3btr165X2vN0qUg9bqw9Jr/LkwfPyB7RSp+G6iet1A7nP/0nhGjx15n/SPrfZ4a1atCsbq29BWdXZBfN+deHfvx5FD5FCJ/HpwhxIT5FiIvwKcLr88svvwwZMqT44I99mtZ3asZnWcf7JT30wQcfWCznF1I/6xLpJ6T7CuC+U7JbA6VPQ7VPjNK5vmK9wqf6Dhw48Omnn+bk5LRo0WL06NGNGjWqvr/PBarnU4QELPg8AhYuRMDCRQhY12f48OHK7vRpd93qbqAcmp60eMsjz7/W/oHfrTtmfJcnMwrlbfWVPjFq34ZK1wZKUMWWSjATyzQAAAATlJaWLl+2bMsj3cvuQA+2Wp7u2uLRdxZ2a/ZEn4bKy7dZ4hpc/dl/gYz/NgAA4DdOFhYFCb1uzd9cLGoaHtzEOJ4eT3K4JvxnAgAAQghx0C6/PSq/PSbXHgk9qwYfs5c0rFOz7NV9BWdvbHqDieP5Ft+5ZAoAADztlxKx6KAxIV1vnqL1+ExLPS4Hxio77q/55PiHJ6/70aHp7t1OnHNO25iVnJxs7rQ+pMINlq7rycnJ//73vy/cOHHixL1797q/HjRo0B//+EfPTAcAADztnCbSjstvjxrfHJPZdtkrWr27kfKHW9Vb6/666uf06dMfO32699z/69aorlPXN58498yfpwwfPtzEsX1LxQLW0qVL165dm5ube+FGKWVubm5KSkpwcLAQouwDnAAAoKrl5uZmZWU1atSoRYsW5Tz2WDPE1gL5zVH57TFje4HsFKHc3Uh9N07tEnn5VRWCg4Pnzp27d++fd+zYUbNmzQ/i4vi8doVULGA1a9asYcOGU6ZMuXDjyZMndV2fPHlybm5ux44dn3766aCgoHIOIqV8++23y9khNjZ2wIABV3rV4XDYbDaWaUAZh8NRXFxs9hTwFpqm8ZbAhRwOR2lpafm/mMr39ttvL1y40IMjeYqu6ydPnnQ6HFZV0aVUrbb69etf9DctNUSJLhy6cOrCqoiaFhFsETdaxGlFLBFiiVmjm6RXr15//etfHQ6HEKKcMHpVNWvWvOq3VyxgdejQ4dKNhYWFLVu2fPzxxxs0aPDee+/NmjVr8uTJ5RxESrljx45ydnA4HH379r3Sqy6Xy+VyXfvM8Hu8JXAhTdN4S+BClX9LHD16tG/fvqNHj/bgVKh+33333bfffut+MyiKUpm3RI0aNa66jwc+RdiqVasZM2a4v05OTr7qHXCqql50C1eFuFwuFhrFhYqLi8PCLn5IOwKWpmmGYfCWQJmgoCCn01mZt0RQUFB0dHTHjh09OBWq36FDh6xWa1hYmKqq1bDQqAc+RfjTTz/t2bPH/TUX7wAAACoVsDIyMoQQJSUlL7300qFDh1wu1/z587t37+6h2QAAAHxSpQLWpEmThBBt27YdOXLklClTHnroIbvdziIZAAAgwF3PPVhr1qy58AtFURISEhISEjw5FwAAgM9iJXcAAAAPI2ABAAB4GAELAADAwwhYAAAAHkbAAgDAn/Xr16/sf1FtCFgAAPi/Bx980OwRAgsBCwAA/8cqldWMgAUAgP8ru1C4bt26CRMmDB06dMmSJUKIwsLCV1555f777x89evS0adNOnjwphNi4ceNjjz2WkJAwbNiwRYsWub/x22+/TUpKMvdv4UM88LBnAAACUIkuHlqrmz3Fr4KtYmEfy1V3y8/Pf//993/44YcXXnghKSlp5syZd99993PPPedyuZYuXfr3v//9tdde+/jjj+++++6kpKTs7Ozf//73w4YNE0Ls2bNnxowZVf/38BMELAAArodNFY+0Usye4ldW5ZqGiY+PVxSlY8eOpaWlQoiMjIzNmzeXvRoWFiaEmD179v79+9esWZORkaFpmvulMWPGhIeHV8Hg/omABQDA9bAoIuEG37vTplatWhf+MSQkZM6cOTExMUIIh8Nht9uFEK+++qrNZuvTp8/48eO/+eYb956kqwrxvXcGAADVo1gTKw8Z645JswepQj169Fi4cGFJSUlRUdGLL764cOFCIcSOHTtGjBjRrVu3bdu2CSF03YuuhPoKGiwAAH4j2y5XHZGrDhsbTsgukUp4DbMHqkpjx46dNWvW6NGjDcOIi4ubMGGCEGLcuHHPPPNM3bp1+/fv37lz5+nTp5s9pu9RpKzWYO50OkNDQ51O53UfoaCgICwszGazeXAq+LS8vDx3uQ0IITRNKyoqioyMNHsQeAuHw+F0Oq96eUszRPoJ+cUR4/PDstAp722sDmqs9ItVQ21i0qRJ0dHRkyZNqp6BUUWWL1/+8ccfL1++3G63K4oSEhJSpaejwQIABK58h/gy11h1WH5zzGgRqgxqrM7trXaKUFQvunkdPomABQAILFKIHQXnLwL+dFre3Ui9t7HydpwtKtjsyeBHCFgAgIBgd4k1R40vjsgvjhhhQcqgxsobt1t6Ris2Pu6FKkDAAgD4qvT09LS0NKfT2aVLl3vuuUe53EJQWWfE54dt6wq0zfnyjgbK4Cbqn9tbm4d66SVAp9OpKEpQUJDZg6CyCFgAAN9TWlr68MMPb/pyZXyr6CCL+qf3/j6tdYfly5fXrVtXCOHURerx83esOzRrv2jtd7eoy/uptb31l15RUdGMGTMWLlyYd+SwIWWzlq3GjBkzceLE2rVrmz0arpO3vtcAALiymTNn/py6eu3o7kEWVQjxhy7N/ufr3U/84Zm7J//riyNy7THj1rrKoMbq4rvUm2qVOJ3O8HDvvcFqz5498YPuvTfSMu/O2BvCWgkhfjp5dm7K7G4LFqz64osmTZqYPeBv9OvXb82aNZXf89qP46O48gwA8D0pKSnPdGvhTldCCFVRnuvectHixd8c1ZOaKgeG2TbEW5/voHaI8NJLgWVOnTo1ZPCgN26L+lP3VjeEnV9jvVVEyGu9b37qBuuQ++5zP9AmoDz77LNmj+ABBCwAgO85ceLEDWG/KaWiQ2raXMUf3F4ysoVav6ZZc1XYm2++eU+kpXvjepe+NLhldEtXwYcfflj9U5lr586dZo/gAVwiBAD4jGJNfJVrLMuRv9S+YW/B2djQXzPWz0Xn6oTXrVOnjonjXYeUlJR/3dHwSq8+eEujt1NSnnjiicqc4quvvvroo480TRs9enRiYuLGjRs/+eSTEydOBAUF3X///cOGDevXr99DDz20evXq+Ph4KeXatWuLiooefPDBESNG9OvXLzk5efHixY0bN37uueeio6PdxywsLHznnXcyMzODg4PbtGmTnJwcERFx5syZd999d+vWrSEhIYmJiReNcdlXLx1m6tSpQojHHnvs4YcfvuilyvxHqH40WAAAb1fkFPMOGEO/0WMWuGbvNeKilOnPPPp6+v5fis8/F+RsqTZl3d7k5GRz56woTdMOZR9sGl7rSjs0r1t73759lTzL7Nmzp02b9s9//nPTpk1CiI8//viuu+5asmTJ66+//u9//9u9z4033jht2rS5c+fWrl37gw8+mDp16vz5890vORyOlJSUNm3avPfee2XHnDlzZq9evRYsWPDee+/FxMT8/e9/F0LMmjVLCDFv3rz3338/KyvrojEu++qlw7zyyitCiPfff/+yc/oQGiwAgJfKKxYrDhlLc4zN+bJvQzWxqfJhT1td95MBbx535sSRfm/+rVt0iFVVvz966r7ho9y/m32IqqqKajGkuNLC8ZphVP7RcO3atfvwww/79ev3xhtvCCFmz569f//+NWvWZGRkaJrm3qdnz55Wq1UIkZCQYLFYOnXq5HK53C8NGDDAYrHcf//9Y8eOLTtmRkbG5s2by/4YFhYmhNiyZcu//vUv9ycfk5OTL7qH/bKvXnYYt3Je8gkELACAdzlol0tz5LIcY98peU9j9bHW6vJ+aq1Lfl+9+OKLycnJmzZtcrlcf+vcuWXLlmYMWymqqra++eb9J8/eGnn5K5t7C862bdu2kmd5+eWXt23b9vXXX69atWr69OmvvvqqzWbr06fP+PHjv/nmG/c+ZTHOYrEIIS5dUUxRFMMwyv4YEhIyZ84c93NgHQ6H3W6/6Lsue4RLv77sMFd9ySdwiRAA4BV2FcqXdxgdlmrdV2oHTsupHS15I23ze1uSbrxMunJr1KjR/fffP3z4cF9MV25jxoz5JPPwlV6dm3n44YcfruQpRo0aFRMTM2rUKPeFuR07dowYMaJbt27btm0TQui6Xv63f/XVV7quL1my5MKo16NHj4ULF5aUlBQVFb344osLFy4UQnTp0uX9998vLi52OByX3pt/2VevNIymaRWd09sQsAAAppFCfJ8vn9uit1ykDVmjn3HJd+IsR0fYZvewDIhVggLgd9STTz653xaZsufopS+9s/Wg0qJD5W/ufvDBBydOnDhp0qQJEyYIIcaNG/fMM888+uijZ86c6dy58/Tp08v/dl3XH3jggR9++OHJJ58s2zh27FjDMEaPHp2cnBwVFeU+8hNPPCGlHDVq1IQJE9q1a3fRcS776mWH6dKly5gxYyo6p7dRpJTVeT6n0xkaGup0Oq/7CAUFBWFhYZW/Jg2/kZeX566pASGEpmlFRUWRkZFmD4LyaIZYf1wuzTFWHJLhQSKxqTK0qdqxatascjgcTqczPDz8uo8wadKk6OjoSZMmeXCqC+Xl5SUkJNQ7eWh4m0a31K+jS7Er//TczCPBrW9LSUmpzOSV50/LgS5fvvzjjz9evny53W5XFCUkJKRKT8c9WACAauLQxNdHjWU58vPDRvNQZWhTde29aqswb18LtKrFxMRs2LBh/vz5//n00x+//tFisbRr127im5OTkpIu+3RF+AQCFgCgap0uFV8cMZbmyDVHjdvqK4lN1dc6W2NrEx1+ZbVaH3nkkUceecTsQS72/PPPmz2CryJgAQCqRL5DrDhkLDtkbDgu74xRE5sqs3vYImqYPRYqok+fPmaP4KsIWAAATzp0Vi7LkctyjMxCObCx+khLdVFfNYT7ZhFgCFgAAA/48dT5XHX4rLzvBvXZ9pa7Gyo1LGaPBZiEgAUAuE5SiO0Fcmm2seyQLNZEwg3KjK6WntGKhdurEPAIWACAitGlSDsul+UYyw/JWlaReIMyr7fltvp84A34FQELAHBNnLr45phclmOsPGQ0CVESm6pfDlRvCSdWAZdBwAIAlMfuEl8eMZbmyNW5RvsIJfEGdUpH6w0h5CqgPAQsAAhER44cEUI0btz4SjsUlIiVh41lOUZqnuwRrSQ2Vd+Os0XWrMYRvUlqaiprfvq63bt3V+fpCFgAEFiWLFnyP//zP67CfEURlvDIGTNmXPi0u9xz8osj8rPDRmqe7BWjPHCjOr+3GhZk4rzme/HFF1966aVjx46ZPQgqpV69enfeeWe1nY6ABQABZPny5U89Muq9e9rf3vBmIcT2vFOPj39YVdV2/ZKW5shlOcZBuxzcRH30JnXxXWpNFlkQQghRp06dN9980+wp4GMIWAAQQN54442Xe7W+veH55wffFhP+Wu+bH37u9dolQwbGqi90VAfGqjbV3BkBf0DAAoAAsnv37rjR3S7cEhdbz/V1+olRNu4wAjyIf6cAQKD4pUSoNWqfKnFduLGoxBVSuxbpCvAsAhYA+DmnLj47bAz7Vm+R4oq4rf+CXUcufPU/u48MGDDArNkAf8UlQgDwW5vy5bwsY9FBo0OEMqal+lEv26m7/tqzR4/CNbuHtIoWQnyWdXzTGTVtyTSzJwX8DQELAPzNobNyXpacd8BQFTGmpbpzqLVx7fPXAENiY3ft3v3WW2/NT0+XUvZ4KOGDP/yhTp065g4M+B8CFgD4CYcmPj9izM0yNp2QSTeqH/a0dI++zOKYISEhkydPNmE+IJAQsADAt+lSfHtMzssyPj9s9I5Rx7VSl9ytBnGHLWAqAhYA+Kq9p2TKQeOTLFnLIsa0VN/samsQbPZMAIQQBCwA8DlFTrE425ibZeScFUlNleX9LO3rscwC4F0IWADgG5y6+PywMfeATDtuxDdRX77N0idGUUlWgFciYAGAt3OvtrA422hfTxnTUl3Q2xZiM3smAOUiYAGAl8o9JxcckB/9ZChCPNRc2TLEemMdCivANxCwAMC72F1iSbbxSZaxu0g+2Eyd19vSJZJcBfgYAhYAeAVDim+Pybn/XW1h4q3qoCastgD4KgIWAJhs3yn5KastAP6FgAUA5igoEQt/NuZmGccdYlQL5YsBlpvDuRQI+AkCFgBUK6cuVh0x5mbJ9XlGfBP1jdstfRuy2gLgbwhYAFBNthfIuVlGykGjRagypqU6r7etDqstAH6KgAUAVevoOfl/2fKjnwynLh5qrmy8z9qM1RYAf0fAAoAq4V5tYW6WsatIPthMndPD0rUBuQoIFAQsAPCkC1dbuDNG/T2rLQABiYAFAJ7hXm1hbpYMZrUFIOARsACgUoqcYnG2MTfLyDkrkpoqS++2dIjgUiAQ6AhYAHA9Llpt4cVOlrtYbQHAfxGwAKBi9hTJeQeMj39itQUAV0TAAoBr4l5t4d8/GSWstgDgaghYAFAehyY+P2LMzTI2nZD3NlZndLXc1UghWAEoHwELAC7DkGLtMTk3y/jssNErRh3bSv2/u9QaFrPHAuAjCFgA8BsXrbYwg9UWAFQcAQsAhGC1BQAeRcACENBKDbE615iXJVfnGgNi1efaq/fEqlYWXgdQOQQsAAFqc76cd8BIOWi0ras83Er9sBerLQDwGAIWgMBy+Kycf0DOzTKEEKNbqtsTrE1CuBQIwMMIWAACwkWrLbwTx2oLAKoQAQuAPzOk2Hji/KXA2yOV0S3UlL5qLX7yAahi/JgB4J/2npJzs4z5B2R0sBjTUn1tmC2yptkzAQgYBCwAfqWgRHycXWPpRu1YsRjVQll9j+WWcK4EAqhuBCwAPmPdunX/+Mc/9u/f37BhwwceeGDChAkWy/m11S9YbUH2jrRO6chqCwDMxI8fAL5h5syZo4fc29eZPbtrZHK4fe5fXoiPjzcMY0+R/NNWvclC17QM4+5GSs4w5V+3n4tvQroCYCYaLAA+oKCg4JWpU5YndW5Rt7YQokXd2t0a1R30aXrzPy2RnRIfaq5svM/arI4ihNA0o8jsaQGAgAXAB2zatOnm8JrudOUWZFHjW0ZlF6YufugBEwcDgMuiQwfgA0pLS4NtF/+8CrZZ61lKTZkHAMpHwALgAzp26vTDiTOnna4LN64/VHDbbbeZNRIAlIOABcDb5djlhL2NrT1GTvj8h5xTxUKIs6XaK6n7jteMGDVqlNnTAcBlELAAeC8pxJx9RpcVWr9G6uHP3x/y++ceWJ3VZvbarvO3uNr3Wb9+fa1atcyeEQAuQ5FSVuf5nE5naGio0+m87iMUFBSEhYXZbDz1Hufl5eXFxMSYPQU8L9suk9P0Yk181Mty8wWLhZ46dSosLEy5woMENU0rKiqKjIysrjHh7RwOh9PpDA8PN3sQeAu73a4oSkhISJWehQYLgNdxF1ddV2j9G6np8dabf7sUe3h4+JXSFQB4CZZpAOBdsu1yfKpeoovUwdbWPOUGgG+iwQLgLcqKqwGxalo86QqAD6PBAuAVDtrl+FS91BBp8dabwohWAHwbDRYAkxlS/H230XWFlthUTRtMugLgD2iwAJhp/2k5LlW3qeL7+6zNQ4lWAPwEDRYAcxhSzNln9PhMu6+JuvZe0hUAv0KDBcAEP5+R41J1Q4iN8daWXBME4Heup8HSdX3s2LEXbrHb7S+88EJiYuKUKVPsdruHZgPgh3Qp3tptdFup3dtYXT+IdAXAP1U4YC1duvQPf/hDbm7uhRtTUlKioqJSUlIaNGiwaNEiz40HwK/8eErGrdSW5Bib7rM+115VCVcA/FSFA1azZs0ufbpqenr6kCFDgoKChgwZkpaW5qHZAPgPzRDTMozen2sjmqvfDbK24I4rAH6twvdgdejQ4dKNJ0+ejIqKEkJERUUVFhaWfwRd1xs3blzODnFxcW+//faVXi0sLHQ6nTyLEGUKCgosFovZU6A8++2WiTtrhdvk6p7FjYKNgl+q8Fyapp0+fbqan7IKb1ZSUuJ0OktLS80eBN7i7NmziqIUFxdf9xHq1atntV4lQXnmJncppfvRYFJKwzDK39lisXz33Xfl7BAcHBwREVHOuXjYMy5UWlpazhsG5tIMMXO3/Pse+UonJfkmRRE1qvyMmqaqKm8JlOFhz7hIUFBQJR/2rKpXvwDomYAVERGRn58fGxtbUFBQv379q+7fvHnz6z6X5b+u+wjwM7wfvNbuIjl2vR5RU2xLsDYJqaZrglJK3hK4EL81cBGLxaIoSlW/JSq7DlZGRoYQolu3bqtXr5ZSrl69Oi4uzhODAfBhLkO8stPou0p74hb1q4HVl64AwEtUNmBNmjRJCDF69OiDBw8OHz48Jydn5MiRnhgMgK/64aTsskLbnG/sTLSOa8VqxgACkVLNt4I6nc7Q0FCn03ndRygoKOAeLFwoLy8vJibG7CkghBCaId7cZczcrb96m2VCa3OilaZpRUVFkZGRppwdXoh7sHARu91eyXuwrgUruQPwjMxCOTZVjwoW2xOssbW5JgggoBGwAFSWyxAzzS6uAMCrELAAVEpmoXxkvR5TS+xIsDaiuAIAIQQBC8B1o7gCgCshYAG4HhmF8pH1eiOKKwC4HAIWgIqhuAKAqyJgAaiAzflyXKreLFTsTLQ2rEVxBQCXR8ACcE1KdPHSDv3D/cZfOlNcAcBVELAAXN33+XJsqt62rvLj/bbImmZPAwBej4AFoDwOTUzZrv/nZ+Ofd1juv5HiCgCuCQELwBVtypfjUvW2dZWMoRRXAFABBCwAl+HQxMs79blZxtt3WJIorgCggghYAC628YQcl6q3q6dkDrXVp7gCgIojYAH4lbu4mpcl345ThzaluAKA60TAAnDehhNyXKrevp6SmWSNqGH2NADgywhYAESxJl7Zqc/Lku/EqYkUVwBQaQQsINClH5fj0vQOFFcA4DkELCBwlRVXs7qrQ26guAIAjyFgAQEq7bgcl6p3jFB2JVnrUVwBgEcRsICA4y6u5h+Qs+LU+yiuAKAKELCAwPL1Uflomt41Utk11FqX4goAqgYBCwgUZ1zime/1NUflBz0t/RspZo8DAP6MqwNAQPgqV7ZdoiX1KJoAABeHSURBVFkUkZlkJV0BQFWjwQL83OlS8ewW/euj8sOelruJVgBQLWiwAH/2Va5st1QTQmQOtZKuAKDa0GAB/sldXK05Kj/qZbmrIdEKAKoVDRbgh748Itsu0YQQGUOtpCsAqH40WIBfOVUqntuirzkqP77T0pdoBQAmocEC/McXR2S7JefvuCJdAYCJaLAAf+Aurr45Kj/pbekTQ7QCAJPRYAE+b9UFd1yRrgDAG9BgAT6srLia19vSm2gFAF6DBgvwVZ8fPl9cZSZZSVcA4FVosADfU+QUf9qqf3tMzu9tuZNoBQDehwYL8DGfHTbaLj1/xxXpCgC8Ew0W4DPyHeKpjfqeIvl/d1m6NSBaAYD3osECfEPKQaP9UlfzULEj0Uq6AgAvR4MFeLt8h3hyo763SK7ob+0SSbQCAB9AgwV4tcXZRvulrhahYkci6QoAfAYNFuClTjjEkxv0/aflyv7W24lWAOBTaLAAb7Q42+iw1NUyTGxPIF0BgO+hwQK8y3GHeHKDnnVafjbA2rk+0QoAfBINFuBF3MVVqzCxLYF0BQA+jAYL8ArHHeKJdP1nu1w1wHob0QoAfBwNFmA+d3F1U7jYlkC6AgB/QIMFmCmvWDyxQc+2yy8GWDsRrQDAX9BgAaZZnG10WOZqHS62JpCuAMCv0GABJsgrFo+l64fOyq8GWjtGEK0AwN/QYAHVSgoxN8vosMx1S12xNYF0BQD+iQYLqD6HzspH0/R8h1g90NqBaAUA/osGC6gOUog5+4xOy7RO9ZWtCaQrAPBzNFhAlcu2y+Q0/Zwm0uKtt4QTrQDA/9FgAVVICjHrR6PLCm1grLqBdAUAAYMGC6gqOXY5Pk0v1kTqYOvNRCsACCQ0WIDnue+46rJC699ITY8nXQFAwKHBAjws2y7Hp+olukgdbG1NtAKAgESDBXiMu7jqukIbEKumxZOuACBw0WABnnHQLsen6qWGSIu33hRGtAKAgEaDBVRWWXE1MFZNHUy6AgDQYAGV8/MZOT5N1wyxId7aimgFABBC0GAB182QYs4+o9tK7Z5YNXUw6QoA8CsaLOB6/HxGjkvVDSE2xltbEq0AAL9FgwVUTFlxdW9jdf0g0hUA4DJosIAK+PGUHLter2ERm+6ztgglWgEALo8GC7gmmiGmZRg9P9OGNlW/G0S6AgCUhwYLuLo9RXJsqh4eJHYmWpuEEK0AAFdBgwWUx11c9VmlJd+krr6HdAUAuCY0WMAV7S6SY9frETXFtgSiFQCgAmiwgMtwF1d9V2mPtla/HEi6AgBUDA0WcLFdhXJcqh5RU2xPtDauTbQCAFQYAQv4lcsQr/9gvPujPr2L5ZFW9LsAgOtEwALOyyyUY1P1qGCxI9EaS3EFAKgEAhYgXIaYucuYuVt/9TbLhNYUVwCAyiJgIdC5i6voYLE9geIKAOAZBCwELoorAEAVIWAhQGUUyrHr9Ya1xI4EayOKKwCARxGwEHAorgAAVY2AhcCyOV+OS9WbhYqdidaGtSiuAABVgoCFQFGii5d26B/uN/7SmeIKAFC1CFgICN/ny3GpevNQkTGU4goAUOUIWPBz7uJqbpbxzzss999IcQUAqA4ELPizTflyXKretq6SMdQWWdPsaQAAAYOABf/k0MTLO/W5Wcbbd1iSKK4AANWLgAU/tPGEHJeqt6unZA611ae4AgBUOwIW/Iq7uJqXJd+OU4c2pbgCAJiDgAX/seGEHJeqt6+nZCZZI2qYPQ0AIIARsOAPyoqrd+LURIorAIDZCFjweVsKbf+brnWguAIAeA0CFnxYsSZe2al/vK/O7F5qwg0UVwAAb8HvJPiqtOOywzLt4Bmx7s5TpCsAgFehwYLvcRdX8w/IWXHqfTeoeXmG2RMBAPAbBCz4mNTjcnyq3jFCyRxqrccdVwAAr0TAgs844xKTNutf58o5PS39G/HAZgCA9+LOFfiG1bmy3RKtRBM7E62kKwCAl6PBgrc74xL/u1n/+qj8V0/L3UQrAIAvoMGCV1udK9su0YQQmUOtpCsAgK+gwYKXOl0qnt2if31UftTLcldDohUAwJdULGDZ7fZp06bt2bOnTZs2zz77bJ06dcpemjhx4t69e91fDxo06I9//KMnx0SA+fKIfHyDPjBWyRxqrWMzexoAACqoYgErJSUlKipq6tSp77///qJFi8aPH+/eLqXMzc1NSUkJDg4WQlgsFs9PisBwqlQ8t0Vfc1T+u5elL8UVAMA3VSxgpaenv/LKK0FBQUOGDJk6dWpZwDp58qSu65MnT87Nze3YsePTTz8dFBR0pYNIKefOnVvOWWJiYnr27HmlV51OZ0lJia7rFZocPuGro8rvN6sDGskt9xohNq2k5Jq+y/2WqOLR4DM0TeMtgQuVlJSUlpbylkCZkpISRVGs1uu/SyooKEhVr3IXe8WOfvLkyaioKCFEVFRUYWFh2fbCwsKWLVs+/vjjDRo0eO+992bNmjV58uQrHURKuWLFinLO0qZNm86dO1/pVYfDYbPZNE2r0OTwcqddyku7gtbnW9/t7OgZqQtNOK75/+GSkhKHw1GV08GXaJrmcDh4S6BMSUmJ0+ks55/9CDTugFWZq20229VvXqlYwJJSKori/sIwfn0+SatWrWbMmOH+Ojk5OTk5uZyDqKq6ZMmSCp33Qrquh4WFXcvfDb5i1RH5eLp+b2Nl1/2WEFuFfwiWlJTUrVu3KgaDL3L/64u3BMo4HA6n0xkeHm72IPAWVqtVUZSQkJCqPUuF9o6IiMjPz4+NjS0oKKhfv37Z9p9++snlct16661CCJvNRvrBNXLfcfXNUTmvt6V3DHdcAQD8RMXWwerWrdvq1aullKtXr46LixNCZGRkCCFKSkpeeumlQ4cOuVyu+fPnd+/evUqGhX/5/PB/17hKspKuAAD+pGIBa/To0QcPHhw+fHhOTs7IkSOFEJMmTRJCtG3bduTIkVOmTHnooYfsdnv5lwiBIqd4LF3/4/f6/N6W93tYarMcGwDAvyhSyuo8n9PpDA0NdTqd132EgoIC7sHyaZ8dNp7YYAxqrMzs5plolZeXFxMT44EDwS9omlZUVBQZGWn2IPAW3IOFi9jtdq+7BwuojHyH+N1GfcdJ+Z8+ll7RXBMEAPgtnkWIarI422i/1BVTS2QOtZKuAAD+jQYLVS7fIZ7cqO8tkiv6W7tEEq0AAP6PBgtVy11ctQgVOxJJVwCAQEGDhapywiGe3KDvPy1X9rfeTrQCAAQSGixUicXZRoelrpZhYnsC6QoAEHBosOBhxx3iyQ161mn52QBr5/pEKwBAIKLBgie5i6tWYWJbAukKABC4aLDgGccd4ol0/We7XDXAehvRCgAQ2Giw4AHu4uqmcLEtgXQFAAANFionr1g8sUHPtssvBlg7Ea0AABBC0GChMhZnGx2WuVqHi60JpCsAAH5Fg4XrkVcsHkvXD52VXw20dowgWgEA8Bs0WKgYKcQH+4z2S13dGijbE0hXAABcBg0WKuDQWflomp7vEF/fY+1AtAIA4AposHBNpBBz9hmdlmmd6itbE0hXAACUhwYLV5djl4+m63aXSI+33hxOtAIA4CposFAed3HVZYV2d0N1A+kKAIBrQ4OFK8qxy/FperEm1g8mWgEAUAE0WLiMsuKqfyOVy4IAAFQUDRYulm2X41P1El2kDra2JloBAFBxNFj4lbu46rpCGxCrpsWTrgAAuE40WDjvoF2OT9VLDZEWb70pjGgFAMD1o8HCr8XVwFg1dTDpCgCAyqLBCnQ/n5Hj03TNEBvira2IVgAAeAINVuAypJizz+i2UrsnVk0dTLoCAMBjaLAC1M9n5LhU3RBiY7y1JdEKAACPosEKOGXF1b2N1fWDSFcAAHgeDVZg+fGUHLter2ERm+6ztgglWgEAUCVosAKFZohpGUbPz7ShTdXvBpGuAACoQjRYAWFPkRybqocHiZ2J1iYhRCsAAKoWDZafcxdXfVZpyTepq+8hXQEAUB1osPzZ7iI5dr0eUVNsSyBaAQBQfWiw/JO7uOq7Snu0tfrlQNIVAADVigbLD+0qlONS9YiaYnuitXFtohUAANWNBsuvuIuru744X1yRrgAAMAUNlv/ILJRjU/WoYLEj0RpLtAIAwDwELH/gMsTMXcbM3fqrt1kmtKaVBADAZAQsn+curqKDxfYEiisAALwCAcuHUVwBAOCdCFi+KqNQjl2vN6wldiRYG1FcAQDgTQhYvofiCgAAL0fA8jE/nJRjU/XY2mJnorVhLYorAAC8EQHLZ5To4qUd+of7jb90prgCAMCrEbB8w/f5clyq3jxUZAyluAIAwNsRsLydu7iam2W8dYflgRsprgAA8AEELK+2KV+OS9Xb1lUyhtoia5o9DQAAuDYELC/l0MTLO/W5Wcbbd1iSKK4AAPApBCxvlH5cjkvTb6+v7EqyRdQwexoAAFBBBCzvUqyJ57fpiw/Kd7urCTdQXAEA4JMIWF5k4wk5NlVvX0/JGGqtzx1XAAD4LAKWV3DfcTUvS74TpyY2pbgCAMC3EbDMl35cjk/T29dTMpOs3HEFAIAfIGCZqVgTr+zU52VxxxUAAH6FgGWatONyfJreoZ6yK8laj+IKAAA/QsAygbu4mn9AvhunDqG4AgDA7xCwqlvacTkuVe8YoWQOpbgCAMA/EbCqzxmX+N/N+ueH5eweanwTiisAAPwWv+aryddHZbslmkMTu5OspCsAAPwbDVaVcxdXXx+VH/S09GukmD0OAACoclQpVWt1rmy7RBNCZA61kq4AAAgQNFhV5XSpeHaL/vVR+WFPy91EKwAAAgkNVpX4Kle2W3q+uCJdAQAQaGiwPMxdXK05Kj/qZbmrIdEKAIBARIPlSV8e+fWOK9IVAAABiwbLM06Viue26N8clZ/0tvSJIVoBABDQaLA84Isjst0STQiRMdRKugIAADRYlVJWXM3tbelNtAIAAEIIGqzK+Pzwf++4SrKSrgAAQBkarOtR5BR/2qp/e0zO7225k2gFAAB+iwarwj47bLRdev6OK9IVAAC4FA1WBeQ7xKTN+sZ8+Z8+ll7RRCsAAHB5NFjXanG20W6pK9gqModaSVcAAKAcNFhXl+8QT23UfyySK/tbu0QSrQAAwFXQYF3F4myj/VJX81CxI5F0BQAArgkN1hWdcIinNur7TsmV/a23E60AAMA1o8G6vMXZRoelrhahYnsC6QoAAFQMDdbFTjjEExv0rNMUVwAA4DrRYP2G+46rVmFiG8UVAAC4XjRY5x13iCfS9Z/t8vMB1s71iVYAAOD60WAJ8d87rm4KF9sSSFcAAKCyAr3ByisWT2zQD9rlqgHW24hWAADAEwK6wVqcbXRc5modLrYlkK4AAIDHBGiDlVcsHt+g59jlFwOsnYhWAADAowKxwVqcbXRY5ro5XGxNIF0BAADPC6wG69BZ+Wianu8QXw20dowgWgEAgCoRKA2WFGLOPqPTMq1TfWVrAukKAABUoYBosA6dlclpekGJ+PZeaweiFQAAqGJ+3mC5i6vbl2t3N1S3JZCuAABAdfDnBivHLpPT9HOaWD/YenM40QoAAFQT/2yw3MVVlxVav0ZqejzpCgAAVCs/bLCy7TI5TXdoInWwtTXRCgAAVDu/arDcxVXXFVr/RmpaPOkKAACYw38arIN2mZyqOw2RFm+9KYxoBQAATOMPDZa7uOq2QhsQq6YOJl0BAACT+ViDdfjw4c2bNzdt2rRDhw42m00IcdAux6fqpRRXAADAa1Q4YNnt9mnTpu3Zs6dNmzbPPvtsnTp1yt/uKUVFRU8++eQ3K5fdXD8k/5zTqBfz/pwPsqJ7Tt6mT2pr+d92qkq4AgAA3qHClwhTUlKioqJSUlIaNGiwaNGiq273lFGjRmkZqRvH9vxPYudvRnX/c+uQgfcOmpP284Z463PtSVcAAMCLVDhgpaenDxkyJCgoaMiQIWlpaVfd7hE///zzprVr/nrXrcFWi3tLv2YNRrSK6J/zcSsuCwIAAC9T4UuEJ0+ejIqKEkJERUUVFhZedfuldF2/5ZZbytmhW7du06dPv3DL9u3bb65fp4blN3GwQ3TYuj17CgoKKvpXgJ8pLCx035AHCCE0TTt9+rSi8E8vnFdSUuJ0OjVNM3sQeIuzZ88qilJSUnLdRwgPD7dar5KgKhywpJTun1xSSsMwrrr9Uqqqzp8/v5wd6tSpExYWduGW2NjY4+ecF+2WZy9p0KTBRXsiABUXF/M2QBlN0wzD4C2BMkFBQU6nk7cEyqiqqihKSEjIdR/BYrFcdZ8KB6yIiIj8/PzY2NiCgoL69etfdfulFEXp1KlThU7apUsXrU7ElwdO3NMiyr3ljFNbsDv3X1OTqC5gs9l4G6CMoii8JXAhd+bmLYEyNpvN/YOiSs9S4XuwunXrtnr1ainl6tWr4+LihBAZGRmX3e5BVqt13rx5U7bk/nntj0v2HZu1LXvAfzY+MP7xAQMGePZEAAAAlVfhgDV69OiDBw8OHz48Jydn5MiRQohJkyZddrtn9erV68d9+9uN/v33dW913XHfoq++nTlzpsfPAgAAUHmKlLI6z+d0OkNDQ53Oi2+ounZ//etfR4wY0aRJEw9OBZ82ceLEt956i5ua4Xbs2LEPP/xwypQpZg8Cb7Fly5bdu3ePGzfO7EHgLZYtWxYcHDxw4MAqPYvvPSpnwYIFJ06cMHsKeJF33323/M9VIKAUFBQsWLDA7CngRfbu3btq1Sqzp4AXSUtL27JlS1WfxfcCFgAAgJcjYAEAAHgYAQsAAMDDKrwOViUpihIZGTlmzJjrPkJhYeGrr74aHh7uwang02rVqvXII49wkzvczpw5c/r06cr8kIGfOXz4cF5eHm8JlNm1a5fNZjtw4MB1H+G111676oftqvtThEKInTt37tq1q5pPCgAA4BHx8fF169Ytfx8TAhYAAIB/4x4sAAAADyNgAQAAeBgBCwAAwMMIWAAAAB7mSwFL1/WxY8eaPQW8RVpaWnJyckJCwtNPP52bm2v2ODDfli1bxo8fn5CQMH78+G3btpk9DrxFdnb24MGDzZ4C5ps4cWK///rHP/5R1aer7nWwrtvSpUvXrl3L71G45eXl/e1vf5s+ffqNN964YsWKv/3tb2+99ZbZQ8FMuq6//vrrU6dO7dChQ3p6+owZMz799FOzh4L5zp49O336dKfTafYgMJmUMjc3NyUlJTg4WAhhsViq+ow+02A1a9Zs1KhRZk8Bb5GXl9e3b9/WrVvXqFGjf//+R44cMXsimEzX9T//+c8dO3YsKSmx2Wy1a9c2eyKYzzCM6dOnjxgxwuxBYL6TJ0/quj558uRhw4a98cYb586dq+oz+kyD1aFDB7NHgBfp1KlTp06dhBC6rn/yySe9e/c2eyKYLCgoqGvXrg6HY8iQIUKIauj/4f0WLlwYGxvbs2dPsweB+QoLC1u2bPn44483aNDgvffemzVr1uTJk6v0jD7TYAGX2rZt21NPPVW7du2nnnrK7FngFYKDg1euXDl27NhZs2aZPQtMtmPHjp07d44fP97sQeAVWrVqNWPGjBYtWoSGhiYnJ1fDbZoELPgkKeWcOXMWLFjwwgsvJCcnV8PVdHi5vLy8Dz74QAgRHBx8zz33HD582OyJYLIdO3ZkZGQMHDiwX79+Qoh+/frt3r3b7KFgmp9++mnPnj3ur202m81mq+ozErDgkzIzMzdt2vTqq69GREQ4HA6Hw2H2RDBZRETEqlWrdu3aJaX87rvvWrRoYfZEMFlycvKa/xJCrFmzpk2bNmYPBdOUlJS89NJLhw4dcrlc8+fP7969e1Wf0WfuwQIulJGRkZubm5iYWLbF/TMUASsoKOjll1+eNWvW8ePHGzduPGnSJLMnAuBF2rZtO3LkyClTppw7d65Lly6/+93vqvqMPOwZAADAw7hECAAA4GEELAAAAA8jYAEAAHgYAQsAAMDDCFgAAAAeRsACAADwMAIWAACAhxGwAAAAPIyABQAA4GEELAAAAA8jYAEAAHgYAQsAAMDDCFgAAAAe9v/RtxxRBnsq/wAAAABJRU5ErkJggg==" }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = log.(range(1, exp(4), length = 10)) .+ 1 # uneven grid\n", "y = log.(x) # corresponding y points\n", "\n", "interp = LinearInterpolation(x, y)\n", "\n", "xf = log.(range(1, exp(4), length = 100)) .+ 1 # finer grid\n", "\n", "plot(xf, interp.(xf), label = \"linear\")\n", "scatter!(x, y, label = \"sampled data\", markersize = 4, size = (800, 400))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this point, `Interpolations.jl` does not have support for cubic splines with irregular grids, but there are plenty of other packages that do (e.g. [Dierckx.jl](https://github.com/kbarbary/Dierckx.jl) and [GridInterpolations.jl](https://github.com/sisl/GridInterpolations.jl))." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Multivariate Interpolation\n", "\n", "Interpolating a regular multivariate function uses the same function" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "hide-output": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "interp_linear(3, 2) = 1.6094379124341003\n", "interp_linear(3.1, 2.1) = 1.6484736801441782\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "interp_cubic(3, 2) = 1.6094379124341\n", "interp_cubic(3.1, 2.1) = 1.6486586594237707\n" ] } ], "source": [ "f(x,y) = log(x+y)\n", "xs = 1:0.2:5\n", "ys = 2:0.1:5\n", "A = [f(x,y) for x in xs, y in ys]\n", "\n", "# linear interpolation\n", "interp_linear = LinearInterpolation((xs, ys), A)\n", "@show interp_linear(3, 2) # exactly log(3 + 2)\n", "@show interp_linear(3.1, 2.1) # approximately log(3.1 + 2.1)\n", "\n", "# cubic spline interpolation\n", "interp_cubic = CubicSplineInterpolation((xs, ys), A)\n", "@show interp_cubic(3, 2) # exactly log(3 + 2)\n", "@show interp_cubic(3.1, 2.1) # approximately log(3.1 + 2.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See [Interpolations.jl documentation](https://github.com/JuliaMath/Interpolations.jl#convenience-notation) for more details on options and settings." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear Algebra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Standard Library\n", "\n", "The standard library contains many useful routines for linear algebra, in\n", "addition to standard functions such as `det()`, `inv()`, `factorize()`, etc.\n", "\n", "Routines are available for\n", "\n", "- Cholesky factorization \n", "- LU decomposition \n", "- Singular value decomposition, \n", "- Schur factorization, etc. \n", "\n", "\n", "See [here](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) for further details." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## General Tools" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### LaTeXStrings.jl\n", "\n", "When you need to properly escape latex code (e.g. for equation labels), use [LaTeXStrings.jl](https://github.com/stevengj/LaTeXStrings.jl)." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "hide-output": false }, "outputs": [ { "data": { "text/latex": [ "an equation: $1 + \\alpha^2$" ], "text/plain": [ "L\"an equation: $1 + \\alpha^2$\"" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LaTeXStrings\n", "L\"an equation: $1 + \\alpha^2$\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### ProgressMeter.jl\n", "\n", "For long-running operations, you can use the [ProgressMeter.jl](https://github.com/timholy/ProgressMeter.jl) package.\n", "\n", "To use the package, you simply put a macro in front of `for` loops, etc.\n", "\n", "From the documentation" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "hide-output": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\r", "\u001b[32mComputing... 14%|█████▌ | ETA: 0:00:07\u001b[39m" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\r", "\u001b[32mComputing... 38%|██████████████▉ | ETA: 0:00:04\u001b[39m" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\r", "\u001b[32mComputing... 58%|██████████████████████▋ | ETA: 0:00:02\u001b[39m" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\r", "\u001b[32mComputing... 78%|██████████████████████████████▍ | ETA: 0:00:01\u001b[39m" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\r", "\u001b[32mComputing... 98%|██████████████████████████████████████▎| ETA: 0:00:00\u001b[39m" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\r", "\u001b[32mComputing...100%|███████████████████████████████████████| Time: 0:00:05\u001b[39m\n" ] } ], "source": [ "using ProgressMeter\n", "\n", "@showprogress 1 \"Computing...\" for i in 1:50\n", " sleep(0.1) # some computation....\n", "end" ] } ], "metadata": { "date": 1580349904.0568776, "download_nb": 1, "download_nb_path": "https://julia.quantecon.org/", "filename": "general_packages.rst", "filename_with_path": "more_julia/general_packages", "kernelspec": { "display_name": "Julia 1.3.0", "language": "julia", "name": "julia-1.3" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.3.0" }, "title": "General Purpose Packages" }, "nbformat": 4, "nbformat_minor": 2 }