{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "view-in-github"
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "XGh1ScX5EZVp",
"tags": []
},
"source": [
"# Interpolation Methods"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "b9K_MnFJEZVr"
},
"source": [
"Due to the discrete, and sometimes sparse, nature of experiments and observations, data taking procedures will always produce discrete data as well. Even, as we have seen before, information only can be discretely presented into a computer due to the binary representation. However, when we are dealing with physical models, continuous and smooth properties are of course preferred. Interpolation techniques allow then to recover a continuous field from sparse datasets. Throughout this section we shall cover some of these interpolation methods.\n",
"\n",
"---\n",
"\n",
"\n",
"## Bibliography: \n",
"[1a] Gonzalo Galiano Casas, Esperanza García Gonzalo [Numerical Computation](https://www.unioviedo.es/compnum/expositive/handbook/metnum.pdf) - [GD](https://drive.google.com/file/d/1gwHYLx5aBf3oXRhmGeG-n1toTtAUhvfX/view?usp=sharing) - [Web page with notebooks](https://www.unioviedo.es/compnum/index.php/english-menu) \n",
"[1g] Mo Mu, [MATH 3311: Introduction to Numerical Methods](https://www.math.ust.hk/~mamu/courses/230/course.htm) [GD](https://drive.google.com/drive/folders/1vv38SS7Zpw8mOHG7Kq_3lZr6o9H4xOYP?usp=sharing) - [Demostration Error in Lagrange Polynomials](https://www.math.ust.hk/~mamu/courses/231/Slides/CH03_1B.pdf) \n",
"[1h] Zhiliang Xu, [ACMS 40390: Fall 2016 - Numerical Analysis](https://www3.nd.edu/~zxu2/ACMS40390-F16.html) [GD](https://drive.google.com/drive/folders/1Z9Fws1v34PuPdGAp_mBCWmTvZM_Qd6Xk?usp=sharing)\n",
"\n",
"\n",
"https://github.com/restrepo/Calculus/blob/master/Differential_Calculus.ipynb\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ooE-OhsOEZVs"
},
"source": [
"- - -\n",
"- [NumPy polynomials](#NumPy-Polynomials)\n",
"- [Linear Interpolation](#Linear-Interpolation)\n",
" - [Steps](#Steps-LI)\n",
" - [Example 1](#Example-1)\n",
"- [Lagrange Polynomial](#Lagrange-Polynomial)\n",
" - [Derivation](#Derivation)\n",
" - [Steps](#Steps-LP)\n",
" - [Activity](#Activity-LP)\n",
"- [Divided Differences](#Divided-Differences)\n",
" - [Example 2](#Example-2)\n",
"- [Hermite Interpolation](#Hermite-Interpolation)\n",
" - [Derivation in terms of divided differences](#Derivation-in-terms-of-divided-differences)\n",
" - [Example 3](#Example-3)\n",
"- - -"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"colab_type": "code",
"id": "8BwVVom3EZVu",
"outputId": "5d379375-42ad-47ce-d4f1-ea45a27b490b"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"%pylab inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "jLrHQj5MEZV0"
},
"outputs": [],
"source": [
"from IPython.display import display, Markdown, Latex, Image \n",
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "qA6vPaJBEZV2"
},
"outputs": [],
"source": [
"from scipy import interpolate"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "rIsrgFTCEZV4"
},
"outputs": [],
"source": [
"import numpy as np\n",
"# JSAnimation import available at https://github.com/jakevdp/JSAnimation\n",
"#from JSAnimation import IPython_display\n",
"from matplotlib import animation\n",
"from IPython.core.display import Image "
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "izA9D35WdJjg"
},
"source": [
"Pretty print inside colaboratory"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "vU_ALGnRhKna"
},
"outputs": [],
"source": [
"import IPython\n",
"\n",
"def setup_typeset():\n",
" \"\"\"MathJax initialization for the current cell.\n",
" \n",
" This installs and configures MathJax for the current output.\n",
" \"\"\"\n",
" IPython.display.display(IPython.display.HTML('''\n",
" \n",
" \n",
" '''))\n",
" \n",
"def Polynomial_to_LaTeX(p):\n",
" \"\"\" Small function to print nicely the polynomial p as we write it in maths, in LaTeX code.\"\"\"\n",
" coefs = p.coef[::-1] # List of coefficient, sorted by increasing degrees\n",
" res = \"\" # The resulting string\n",
" for i, a in enumerate(coefs):\n",
" if int(a) == a: # Remove the trailing .0\n",
" a = int(a)\n",
" if i == 0: # First coefficient, no need for X\n",
" if a > 0:\n",
" res += \"{a} + \".format(a=a)\n",
" elif a < 0: # Negative a is printed like (a)\n",
" res += \"({a}) + \".format(a=a)\n",
" # a = 0 is not displayed \n",
" elif i == 1: # Second coefficient, only X and not X**i\n",
" if a == 1: # a = 1 does not need to be displayed\n",
" res += \"x + \"\n",
" elif a > 0:\n",
" res += \"{a} \\;x + \".format(a=a)\n",
" elif a < 0:\n",
" res += \"({a}) \\;x + \".format(a=a)\n",
" else:\n",
" if a == 1:\n",
" # A special care needs to be addressed to put the exponent in {..} in LaTeX\n",
" res += \"x^{i} + \".format(i=\"{%d}\" % i)\n",
" elif a > 0:\n",
" res += \"{a} \\;x^{i} + \".format(a=a, i=\"{%d}\" % i)\n",
" elif a < 0:\n",
" res += \"({a}) \\;x^{i} + \".format(a=a, i=\"{%d}\" % i)\n",
" return \"$\" + res[:-3] + \"$\" if res else \"\" "
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "udTzmy-CEZV7"
},
"source": [
"- - -"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "Z-txzMoJEZV7",
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"## NumPy Polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "bA0x7_nCEZV9"
},
"source": [
"In Numpy there is an implementation of Polynomials. The object is initialized giving the polynomial coefficients: \n",
"\n",
"More information about this\n",
"* [LaTeX print](https://perso.crans.org/besson/publis/notebooks/Demonstration%20of%20numpy.polynomial.Polynomial%20and%20nice%20display%20with%20LaTeX%20and%20MathJax%20(python3).html): "
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"colab_type": "code",
"id": "OEDBQReUEZV9",
"outputId": "dc526b35-cc2c-414b-b2e9-7248e6b93397"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"1 x + 2 x - 3\n"
]
}
],
"source": [
"p = np.poly1d([1, 2, -3])\n",
"print(p)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\thola\n",
"mundo\\\n"
]
}
],
"source": [
"print('\\thola\\nmundo\\\\')"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'$(-3) + 2 \\\\;x + x^{2}$'"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Polynomial_to_LaTeX(p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Copy an Paster in a MarkDown cell: $(-3) + 2 \\;x + x^{2}$"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$(-3) + 2 \\;x + x^{2}$"
],
"text/plain": [
""
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Latex(Polynomial_to_LaTeX(p))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The numpy polynomial is automatically a function of its variable $x$"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.29"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p(1.3)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "HuqkJBhWEZV_"
},
"source": [
"By default, the assigned the attribute `variable` is `x`:"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"colab_type": "code",
"id": "jT6FBw47EZWA",
"outputId": "88b68268-2d6d-412c-9fbb-b5b051e0cf8d"
},
"outputs": [
{
"data": {
"text/plain": [
"'x'"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p.variable"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "x-7Tef0-EZWC"
},
"source": [
"which can be assigned at initialization "
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"colab_type": "code",
"id": "A8ZMzlaTEZWD",
"outputId": "8df22d7e-3db3-4d15-b61f-566993847465"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"1 t + 2 t - 3\n"
]
}
],
"source": [
"q = np.poly1d([1, 2, -3],variable='t')\n",
"print(q)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Change of variable"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"1 x + 2 x - 3\n"
]
}
],
"source": [
"p = np.poly1d([1, 2, -3])\n",
"print(p)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 1, 2, -3])"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p.coef"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"1 t + 2 t - 3\n"
]
}
],
"source": [
"pp=np.poly1d(p.coef,variable='t')\n",
"print(pp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Polynomial can be added but not multiplied, simplified or expanded"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p1(x)= \n",
"1 x + 1\n",
"********************\n",
"p2(x)= \n",
"-1 x + 1\n",
"********************\n",
"p1(x)+p2(x)= \n",
"2\n"
]
}
],
"source": [
"p1=np.poly1d([1,1])\n",
"print('p1(x)={}'.format(p1))\n",
"print('*'*20)\n",
"p2=np.poly1d([-1,1])\n",
"print('p2(x)={}'.format(p2))\n",
"print('*'*20)\n",
"print('p1(x)+p2(x)={}'.format(p1+p2))"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"poly1d([2])"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p1+p2"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
"2\n"
]
}
],
"source": [
"print(p1+p2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "hDgtm5VCEZWF"
},
"source": [
"The object have in particular methods for \n",
"__Integration__:"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"1 x + 2 x - 3\n"
]
}
],
"source": [
"print(p)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"colab_type": "code",
"id": "nMJq0CrOEZWG",
"outputId": "836fdb43-e142-4740-fdbf-94221970b2ad"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 3 2\n",
"0.3333 x + 1 x - 3 x\n"
]
}
],
"source": [
"p = np.poly1d([1, 2, -3])\n",
"print( p.integ() )"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.33333333, 1. , -3. , 0. ])"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p.integ().coef"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "JVolgx0wEZWI"
},
"source": [
"__Derivatives__"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"colab_type": "code",
"id": "fVS5U1bBEZWJ",
"outputId": "3e77e527-c41f-45dc-b7a2-0f575d102095"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
"2 x + 2\n"
]
}
],
"source": [
"print( p.deriv() )"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "R8f3r0yOEZWN"
},
"source": [
"__roots__:"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 70
},
"colab_type": "code",
"id": "Md4Q2gmNEZWO",
"outputId": "8a03f49c-1d06-4170-8c07-58dfc99e947a"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" \n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-3. 1.]\n"
]
}
],
"source": [
"setup_typeset() #active colab pretty print\n",
"print(p.roots)"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$p(-3.0)$=1.7763568394002505e-15"
],
"text/plain": [
""
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Latex( '$p({})$={}'.format(round(p.roots[0],1),\n",
" p((p.roots[0] ) ) ) )"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$p(-3.0)$=1.7763568394002505e-15"
],
"text/plain": [
""
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Latex( f'$p({round(p.roots[0],0)})$={p(p.roots[0])}' )\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "1NLa6wh1EZWQ"
},
"source": [
" It is possible to define polynomial by given the list of roots and "
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"colab_type": "code",
"id": "E67K42bQEZWR",
"outputId": "81efc332-dc52-471c-c562-81ff1a4fb00f"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 4 2\n",
"1 x - 6.221e+04 x + 9.698e+07\n"
]
}
],
"source": [
"p=np.poly1d([-246.2,-40,40,246.2],r=True)\n",
"print(p)"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.4901161193847656e-08"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p(-40)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "GSuuNnvMEZWU"
},
"source": [
"For further details check the official help:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "n5GOuKoeEZWV"
},
"outputs": [],
"source": [
"np.poly1d?"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "3cs7-2egEZWW"
},
"source": [
"__Activity__: Movement with uniform acceleration\n",
"1. Define a polynomial for the movement with uniform acceleration:\n",
"\\begin{align}\n",
"x(t)=x_0+v_0 (t-t_0)+\\tfrac{1}{2} a (t-t_0)^2 \\,,\n",
"\\end{align}\n",
"2. Use the previous formula expressed as polynomial of degree 2, to solve the following problem with `np.poly1d`: \n",
" * A car departs from rest with a constant acceleration of $6~\\text{m}\\cdot\\text{s}^{-2}$ and travels through a flat and straight road. 10 seconds later a second pass for the same starting point and in the same direction with an initial speed of $10~\\text{m}\\cdot\\text{s}^{-1}$ and a constant acelleration of $10~\\text{m}\\cdot\\text{s}^{-2}$. Find the time and distance at which the two cars meet.\n",
" \n",
"_Hint_. \n",
"\\begin{align}\n",
"x(t)=x_0-v_0t_0+\\frac{1}{2}at_0^2 +(v_0-at_0)t+\\tfrac{1}{2} a t^2 \n",
"\\end{align}\n"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [],
"source": [
"x0=0\n",
"t0=0\n",
"v0=0\n",
"a=6 #m/s^2\n",
"x_1=np.poly1d([0.5*a,v0-a*t0,x0-v0*t0+0.5*a*t0**2],variable='t')\n",
"x0=0\n",
"t0=10 #s\n",
"v0=10 #m/s\n",
"a=10 #m/s^2\n",
"x_2=np.poly1d([0.5*a,v0-a*t0,x0-v0*t0+0.5*a*t0**2],variable='t')"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"3 t\n"
]
}
],
"source": [
"print(x_1)"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"5 t - 90 t + 400\n"
]
}
],
"source": [
"print(x_2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "smN1y6pyEZWX"
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sea el tiempo de encuentro definido $t_e$\n",
"$$x_1(t_e)=x_2(t_e)\\,,$$\n",
"que se interpreta como\n",
"$$x_1(t_e)-x_2(t_e)=0\\,,$$\n",
"Puedo definir un nuevo polinomio\n",
"$$\n",
"x(t)=x_1(t)-x_2(t)\n",
"$$\n",
"Entonces $t_e$, es la raíz del polinomio $x(t)$\n"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"-2 x + 90 x - 400\n"
]
}
],
"source": [
"x=x_1-x_2\n",
"print(x)"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([40., 5.])"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x.roots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La raíz 5 no es física porque el tiempo incial para el segundo carro es 10 "
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_1[40]"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_2[40]"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" \n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$x_1(t)=$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"3 t\n"
]
},
{
"data": {
"text/latex": [
"$x_2(t)=$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"5 t - 90 t + 400\n"
]
},
{
"data": {
"text/latex": [
"meeting time $t_{\\rm end}=$ 40 s; meeting point $x_{\\rm end}=$ 4800 m"
],
"text/plain": [
""
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEKCAYAAADTgGjXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyCklEQVR4nO3dd3xUZdr/8c+VQhJCCCWQQDq919BBIyiCumIXFMXKWlbXn6uurrtu0+fRdVdZy+7KqqtYwK7YQFCCAiK911BDIHRSSE+u3x8ZlsADJJlkciYz1/v1mhdzysx8uQlz5Zz7nPsWVcUYY4xxR4DTAYwxxjRcVkSMMca4zYqIMcYYt1kRMcYY4zYrIsYYY9xmRcQYY4zbgpwOUJ+ioqI0KSnJrdceP36c8PDwug3UgFl7nMra4yRri1P5QnssX778kKq2OtM2vyoiSUlJLFu2zK3XpqWlkZqaWreBGjBrj1NZe5xkbXEqX2gPEdl1tm12OssYY4zbrIgYY4xxmxURY4wxbrMiYowxxm1eV0REJFBEVorIF67lZBH5SUTSReQ9EWnkWh/iWk53bU9yNLgxxvghrysiwC+BjZWWnwGeV9UOwFHgdtf624GjrvXPu/YzxhhTj7yqiIhIHHAp8KprWYCRwIeuXd4ErnA9H+daxrV9lGv/One8qJQFmSWUlJV74u2NMcajvlmfRcaRfI+8t7fdJzIFeASIcC23BI6paqlreQ8Q63oeC2QAqGqpiGS79j9U+Q1FZDIwGSA6Opq0tLQah1p1oJRX1xbTOOg7+kV7W5M5Iy8vz6229FXWHidZW5zK6fYoLFUemJdP/+gg7uwVUufv7zXfiCJyGXBAVZeLSGpdva+qTgWmAqSkpKg7N/0MLyvnjfWzWF8QyYOpA+oqWoPmCzdQ1SVrj5OsLU7ldHu8vyyDwrI1PHD5AAYktajz9/em01nDgMtFZCcwg4rTWH8HmonIiWIXB2S6nmcC8QCu7ZHAYU8ECwoMYHhsEPM2HyAru9ATH2GMMR7x/tIM2rUKJyWxuUfe32uKiKo+pqpxqpoEjAe+U9UbgXnANa7dJgGfuZ7PdC3j2v6denCu3xGxQZQrfLRij6c+whhj6lT6gVyW7TrK9SnxeKjL2HuKyDn8GnhQRNKp6PN4zbX+NaCla/2DwKOeDBEdHsDgdi14b2kG5eU2L70xxvu9tzSDoADhqn5xHvsMrywiqpqmqpe5nm9X1YGq2kFVr1XVItf6QtdyB9f27Z7ONX5AAruP5LN4h0fOmhljTJ0pLi3n4xWZjOramlYRdd+hfoJXFhFvNaZHDBGhQby3NMPpKMYYc07fbtzP4ePFjB+Q4NHPsSJSA6HBgVzZN5av12WRnV/idBxjjDmrGUsziGkaynmdzjgNSJ2xIlJD16XEVxwmrrQOdmOMd8o8VsD3Ww9ybUocgQGe6VA/wYpIDfWIjaRnbCQzlmTgwYvBjDHGbSdOuV+XEu/xz7Ii4oYJAxPYvD+XlRnHnI5ijDGnKC0r54NlGYzo2Ir4Fo09/nlWRNxweZ+2NG4UyPSfdjsdxRhjTjF/y0H2ZRdyw0DPH4WAFRG3NAkJYlyftnyxZh85hdbBbozxHtOXZBDVJIRRXaPr5fOsiLhp/IAECkrK+GzVXqejGGMMAFnZhXy3aT/XpsQRHFg/X+9WRNzUKy6Sbm2aMv2n3dbBbozxCh8sy6BcYfyA+jmVBVZE3CYiTBiUwIZ9OazNzHY6jjHGz5WVKzOWZjC8QxSJLcPr7XOtiNTCuD5tCQsO5F3rYDfGOOz7rQfJPFbA+HrqUD/BikgtNA0N5vLebfls1V7rYDfGOOqdxbuJatKI0d1i6vVzrYjU0o2DXR3sKzOr3tkYYzxgX3YB323az3Up8TQKqt+vdSsitdQrrhk9YyN5xzrYjTEOmbEkA6XiRuj6ZkWkDtw4KIFNWbms2H3U6SjGGD9TWlbOjKW7Oa+e7lA/nRWROvCz3m2JCAnincXWwW6MqV/fbjrA/pwibhxU/0chYEWkToSHBHFlv1i+WLuPo8eLnY5jjPEj7/y0m5imoYzs0tqRz7ciUkduGJRAcWk5Hy63IeKNMfVj9+F8vt9ykOsHxBNUT3eon86KSB3pEtOUlMTmvPPTLpuD3RhTL975aReBAeJIh/oJVkTq0E1DEtl5OJ8F6YecjmKM8XGFJWW8tyyD0d2iiYkMdSyHFZE6NKZHDC3DG/HW4l1ORzHG+Lgv1+zjWH4JNw1OdDSHFZE6FBIUyPUD4vl2434yjxU4HccY48PeWryLdq3CGdK+paM5rIjUsRsGJaBgE1YZYzxm7Z5sVmUc46bBiYh4dg71qlgRqWNxzRszqktrZizdTXFpudNxjDE+6O3FuwgLDuSqfnFOR7Ei4gkTBydyKK+Yr9ftczqKMcbHZOeX8NnqTK7o25bIsGCn41gR8YTzOrYisWVj3vrROtiNMXXrg+UZFJaUM9HhDvUTrIh4QECAcNPgRJbtOso6m7DKGFNHysuVaT/uYkBSc7q3jXQ6DmBFxGOuTYknLDjQjkaMMXVm/paD7D6Sz81DkpyO8l9WRDwkMiyYK/vF8umqTBtPyxhTJ95YtJPWESGM6VG/E0+dixURD7p5SCJFpeW8vyzD6SjGmAZux6HjzN9ykBsHJRLs0DhZZ+I1SUQkXkTmicgGEVkvIr90rW8hInNEZKvrz+au9SIiL4hIuoisEZF+zv4N/q8uMU0ZlNyCtxbvoszG0zLG1MK0H3cSHChMGFS/c6hXxWuKCFAK/EpVuwGDgXtFpBvwKPCtqnYEvnUtA4wFOroek4F/1n/kqt0yNIk9Rwv4btMBp6MYYxqo40WlfLhsD5f0bEPrCOfGyToTrykiqrpPVVe4nucCG4FYYBzwpmu3N4ErXM/HAdO0wmKgmYi0qd/UVbuoWzRtIkN5Y9EOp6MYYxqoj1fsIbeo1Ks61E8Qb5wXXESSgO+BHsBuVW3mWi/AUVVtJiJfAE+r6gLXtm+BX6vqstPeazIVRypER0f3nzFjhluZ8vLyaNKkiVuv/WJbMR9uLeHJYWHERXhN3a6V2rSHL7L2OMna4lS1bY9yVX6zoICwQOGJIaGODHNywQUXLFfVlDNtC6rvMFURkSbAR8ADqppTucFUVUWkRlVPVacCUwFSUlI0NTXVrVxpaWm4+9reA4r5/H+/ZX1JKyam9nTrPbxNbdrDF1l7nGRtcaratsf8LQfJOr6E56/vxQV9nR/m5HRe9WuxiARTUUDeUdWPXav3nzhN5frzROdCJlC5hynOtc7rNA9vxJV9Y/lk5R6O5dvlvsaY6vvPwh20igjh0p5tnY5yRl5TRFynql4DNqrqc5U2zQQmuZ5PAj6rtP5m11Vag4FsVfXawapuGZZEYUk505fY5b7GmOrZdjCPtM0HmTgokUZBXvN1fQpvSjUMuAkYKSKrXI9LgKeBi0RkK3ChaxngK2A7kA78G7jHgczV1iWmKUPateStH3dSWmaj+xpjqvbmop00CgzghkHOTX9bFa/pE3F1kJ+tx2jUGfZX4F6Phqpjtw5LYvJby5m9fj+X9vK6C8mMMV4ku6CED5fv4bLebWgVEeJ0nLPypiMRnzeqazQJLRrz+kK73NcYc24fLMsgv7iM24YlOx3lnKyI1KPAAOGWoUks33WUVRnHnI5jjPFSpWXl/GfhTgYmt6BHrHeM1ns2VkTq2XUD4okICeK1BXY0Yow5s9nr95N5rIDbh3v3UQhYEal3TUKCGD8wnq/W7iPzWIHTcYwxXui1BdtJbNmYC7tGOx2lSlZEHDBpaBKqyrRFO52OYozxMit2H2XF7mPcOjSJwID6vzu9pqyIOCCueWPG9mzDu0t2c7yo1Ok4xhgv8tqCHUSEBnFtineN1ns2VkQccvvwZHILS/nA5hoxxrjsOZrPrHVZ3DAwgfAQr7kD45ysiDikX0Jz+iU047WFO2yuEWMMAG8s3AnAzUOTHM1RE1ZEHHTniHZkHClg9vosp6MYYxyWXVDC9CW7uaxXG2KbhTkdp9qsiDhodPcYEls25pXvt+ONQ/IbY+rPjCW7OV5cxp0j2jkdpUasiDgoMEC4Y3gyqzOOsWzXUafjGGMcUlxacXPhsA4tvf7mwtNZEXHYNf3jad44mKnfb3c6ijHGIZ+v3ktWTmGDOwoBKyKOC2sUyE2DE5m7cT/bDuY5HccYU89UlX//sJ3O0RGc36mV03FqzIqIF7hpSBLBgQG8+oMNhWKMv/lh6yE2ZeVyx4hkR6a+rS0rIl6gVUQIV/eL46MVeziQW+h0HGNMPfrX/G20jghhXJ9Yp6O4xYqIl7hzRDIlZeX/vU7cGOP71uw5xqJth7l9eLLXzlxYlYaZ2ge1a9WEMd1jeGvxLvJsKBRj/MIr87cTERrk1TMXVqXKIiIiLarxaFYPWX3eXee3J7ewlOk/7XY6ijHGw3YeOs7X6/YxcXAiEaHBTsdxW3UGZ9nrepyrxycQaLil1Ev0jm/GkHYteW3BDiYNTWqwh7fGmKpN/WE7QYEB3DosyekotVKdb6mNqtpOVZPP9gAOezqov7grtT1ZOYV8uirT6SjGGA85kFvIh8v3cHW/OFpHhDodp1aqU0SG1NE+phrO6xhFtzZNeWX+NsptYEZjfNIbC3dSUlbO5PMa3s2Fp6uyiKhqldecVmcfUz0iwt2p7dl28DjfbLCBGY3xNdkFJbz14y7G9oghOSrc6Ti1Vu2T7iKSIiKfiMgKEVkjImtFZI0nw/mrS3q2IallY16et80GZjTGx7y9eBe5RaXck9rB6Sh1oiY9t+8A/wGuBn4GXOb609SxwICKo5G1mdksSD/kdBxjTB0pKC7j9QU7SO3cqsENtHg2NSkiB1V1pqruUNVdJx4eS+bnruwbR5vIUF6el+50FGNMHXlv6W4OHy/m3gt84ygEalZEfi8ir4rIBBG56sTDY8n8XKOgAO4c0Y7F24+wfNcRp+MYY2qpuLScqd9vZ2BSCwYktXA6Tp2pSRG5FegDjKHiNNaJU1rGQ8YPrBgm/uV525yOYoyppU9XZbI3u5C7L2jvdJQ6VZOZ4AeoamePJTH/R+NGQdw2LJm/zdnC+r3ZdG/rG+dQjfE3ZeXKP9O20a1NU1Ib4HDv51KTI5FFItLNY0nMGU0alkREaBAvfWd9I8Y0VF+s2cuOQ8e5f1SHBjnc+7nUpIgMBlaJyGa7xLf+NA0N5tahSXy9LovNWblOxzHG1FC5Ki99l07n6AhGd4txOk6dq0kRGQN0BEbjRZf4isgYV2FLF5FHnc7jCbcNTya8USAv2ZVaxjQ4y/eXsfVAHveO7EBAgG8dhUANikjly3q95RJfEQkEXgbGAt2ACb54yq1Z40bcNCSJL9bstSl0jWlAysuVmdtKaNcqnEt7tnE6jkdUZyj4FXWxj4cMBNJVdbuqFgMzgHEOZfGoO0YkExIUYPeNGNOAzN24n4zccu5N7UCgDx6FQPWuzupaRd+HAE5dNhQLZFRa3gMMciiLR0U1CeHGQYm8sWgn94/sSJIPjLljjC9TVV78Lp1WYcK4Pm2djuMx1SkiXaqxT1ltg3iKiEwGJgNER0eTlpbm1vvk5eW5/dq60jOonACUx99dwJ29QhzN4g3t4U2sPU6ytqiw8kApazOLuLGjsuCH752O4zFVFhGn+z2qkAnEV1qOc637L1WdCkwFSElJ0dTUVLc+KC0tDXdfW5fWlGzgzR938tQNAxw9GvGW9vAW1h4nWVtUHIX89aUFJLQI5IJkfLo9GvrUeUuBjiKSLCKNgPHATIczedRdqe0IChBetPtGjPFaczceYF1mDveN7ECQj/aFnNCgi4iqlgK/AGYDG4H3VXW9s6k8q3VEKBMHJ/LJyj3sOHTc6TjGmNOoKlPmbiGxZWOu7BvrdByPq3EREZFw16W1XkFVv1LVTqraXlWfcjpPffj5+e1oFBTAi99tdTqKMeY0czbsZ/3eHH5xQQeCAhv07+nVUp1LfANE5AYR+VJEDgCbgH0iskFEnhUR3xnTuIFoHRHKxEGJfLoyk+1234gxXqO8XJkyd6vfHIVA9Y5E5gHtgceAGFWNV9XWwHBgMfCMiEz0YEZzBj8/vz0hQYFMmWtHI8Z4i1nrs9iwL4f7R3b0i6MQqF4RuVBV/wzkqGr5iZWqekRVP1LVq4H3PJbQnFGriBBuGZbE52v22phaxniBsnLl+TlbaN8qnCv85CgEqlFEVLXE9fTj07eJyODT9jH16OfntaNJoyCen7PF6SjG+L3PV+9l64E8Hryos8/enX4m1ekTuU5EngYiRKSriFR+zVTPRTNVada4EbcNT2bW+izW7sl2Oo4xfqukrJwpc7fQJSaCsT18b6Tec6nO6ayFwAagOfAckC4iK0TkC6DAk+FM1W4fkUxkWDDPzdnsdBRj/NbHK/aw83A+vxrd2SdH6j2X6tyxnglME5FtqroQQERaAklUXKllHNQ0NJifn9+Ov8zazPJdR+if6DtzNxvTEBSVlvHCt+n0jovkwq6tnY5T76pzOksAThQQ1/PDqrpcVY9X3sc445ahSUQ1CeGZWZtRVafjGONX3lm8m8xjBTx8cRefm7WwOqp1ia+I3CciCZVXikgjERkpIm8CkzwTz1RH40ZB3D+qA0t2HGH+loNOxzHGb+QVlfLyvHSGtm/J8I5RTsdxRHWKyBgqRumdLiJ7XTcZbge2AhOAKar6hgczmmoYPyCB+BZh/GXWZsrL7WjEmPrw2g87OHy8mEfGVGewc99UnUt8C1X1H6o6DEgERgH9VDVRVe9U1ZUeT2mq1CgogAcv6sSGfTl8uXaf03GM8XlHjhfz7x+2c3H3aPrEN3M6jmNqdEulqpao6j5VPeahPKYWLu8dS5eYCP72zWZKysqrfoExxm3/mJdOfnEpD43u7HQUR9X6vnwR+XVdBDG1FxggPDS6MzsP5/Pe0oyqX2CMcUvmsQKmLd7FVf3i6Bgd4XQcR1VnZsNTiMj7lReBPsAzdRXI1M6orq0ZkNScKXO3cmXfWMJDavxPbIypwt++qbgv6/9d1MnhJM5z50gkR1Wvcz2uBebWdSjjPhHhsUu6ciiviH//sN3pOMb4nA17c/hkZSa3Dk0itlmY03Ec504ROX3OjsfrIoipO/0SmjO2RwxTv9/Owdwip+MY41OenrWJpqHB3JNqs2BADYqIiPxdRERVd1Rer6pH6j6Wqa2HL+5McWk5L3xrQ8UbU1cWbD3E91sO8osLOhDZONjpOF6hJkciucBMEQkHEJGLRWRhFa8xDmnXqgkTBibw7pLdNnGVMXWgvFx5etZGYpuFcdOQRKfjeI1qFxFV/S0wHUhzFY8HgUc9FczU3v2jOhIWHMjTX9sQZ8bU1qerMlmXmcNDF3ciNNhrZgh3XE1OZ40C7gSOA1HA/ar6g6eCmdprFRHC3ant+WbDfhZvP+x0HGMarILiMp6dvZlecZGM6+0/E05VR01OZz0O/E5VU4FrgPdEZKRHUpk6c/vwZNpGhvLklxtsOBRj3PTqD9vZl13Iby/t5ndDvVelJqezRqrqAtfztcBY4ElPBTN1IzQ4kEfGdGFdZsVlicaYmjmQW8g/529jTPcYBibbVAunc/uOdVXdR8U4WsbLXd67Lb3iInl29mYKisucjmNMg/L8nC2UlJXz6Fj/HWTxXGo17Imq2syGDUBAgPDbS7uRlVPIK99vczqOMQ3Ghr05vLc0g5sGJ5EUFe50HK9U4yIiIj/zRBDjWQOTW3Bpzzb8a/429h6z2m9MVVSVP32xnsiwYH45qqPTcbxWXdyxbhqIR8d2QRWemWWX/BpTlVnrsli8/Qi/Gt3Zbiw8B3eKiF2a0EDFt2jM5PPa8dmqvSzbaQMNGHM2hSVlPPXVRrrERDB+QLzTcbyaO0XErhNtwO5ObU9M01D++Lld8mvM2by2YAd7jhbwxGXdCAqs9YwZPs1ax880bhTEo2O7sDYzmw9X7HE6jjFeJyu7kJfnpTOmewxDO/jnvOk1YUXED43r05b+ic155utNZBeUOB3HGK/y1FcbKS1XfnNJV6ejNAjuFJH9dR1CRJ4VkU0iskZEPhGRZpW2PSYi6SKyWUQurrR+jGtduojYGF41ICL88fLuHM0v5vk5W5yOY4zX+HHbYT5fvZe7z29PQsvGTsdpEGpcRFT1Ig/kmAP0UNVewBbgMQAR6QaMB7oDY4B/iEigiAQCL1Nx13w3YIJrX1NNPWIjuXFQItN+3MnGfTlOxzHGcSVl5fxh5nrimodxd2p7p+M0GF5xOktVv1HVUtfiYiDO9XwcMENVi1zzmKQDA12PdFXdrqrFwAzXvqYGfjW6E5Fhwfz+s/WoWie78W/TftzF5v25PHFZNxultwa8ooic5jbga9fzWCCj0rY9rnVnW29qoFnjRvx6TBeW7DzCZ6v2Oh3HGMccyC1kypwtnN+pFRd1i3Y6ToMSVNMXuCalKlTVGg3CJCJzgZgzbHpcVT9z7fM4UAq8U9Nc5/jcycBkgOjoaNLS0tx6n7y8PLdf681aq9IuMoAnPllN8KEthAdX7zYgX20Pd1l7nNQQ2+JfqwspKC5jbHQe8+fPr9P3bojtURNVFhERCaCiX+JGYABQBISIyCHgS+AVVU2v6n1U9cIqPucW4DJglJ48t5IJVL7TJ861jnOsP/1zpwJTAVJSUjQ1NbWqqGeUlpaGu6/1dq07ZXP5Swv4Kb8VfxrXo1qv8eX2cIe1x0kNrS0Wph9i8ayfuH9UR8Zf1KnO37+htUdNVed01jygPRWd3TGqGq+qrYHhVPRfPCMiE2sTQkTGAI8Al6tqfqVNM4HxIhIiIslAR2AJsBToKCLJItKIiiI3szYZ/FmP2EhuHpLEW4t3sTrjmNNxjKk3RaVl/O7TdSS2bMw91pnuluoUkQtV9c9AjqqWn1ipqkdU9SNVvRp4r5Y5XgIigDkiskpE/uX6jPXA+8AGYBZwr6qWuTrhfwHMBjYC77v2NW761ehOtGoSwuOfrqXM7mQ3fuKV+dvZfug4fxrXwzrT3VRlEVHVE3ejfXz6NhEZfNo+blHVDq4jnD6ux12Vtj2lqu1VtbOqfl1p/Veq2sm1zQaFrKWI0GB+d1k31mXmMO3HnU7HMcbjdh0+zkvz0rm0VxvO79TK6TgNVpVFRESuE5GngQgR6erqIzlhqueimfp2Wa82nNepFX+dvdmGizc+TVX5zSdrCQkM4InL7Baz2qjO6ayFVJxOag48B6SLyAoR+QKwbxofIiI8dUUPylR5wu4dMT7s4xWZLEw/zCNjuxDdNNTpOA1alVdnqWomME1EtqnqQgARaQkkATYxhY+Jb9GYBy/qxP98tYlZ67IY27ON05GMqVOH84p48ssN9E9szo0DE5yO0+BV53SWAJwoIK7nh1V1uaoer7yP8Q23DUume9umPDFzvQ3QaHzOn7/YQF5RKU9f1ZOAAPvqqq1qXeIrIveJyCklW0QaichIEXkTmOSZeMYJQYEBPH1VLw7nFfH013awaXzH/C0H+XTVXu5O7UDH6Ain4/iE6hSRMUAZMF1E9orIBhHZAWwFJgBTVPUND2Y0DugZF8kdI9oxfcluFqUfcjqOMbWWW1jCYx+toX2rcLsnpA5V5xLfQlX9h6oOAxKBUUBfVU1U1TtVdaXHUxpHPHhRJ5Kjwvn1x2s4XlRa9QuM8WLPzNrEvpxC/nJNb7snpA7VaABG1/0g3YFnRaQ3/HdsKuODQoMDeebqXmQcKeDZ2ZudjmOM237cdpi3F+/m9mHJ9E9s7nQcn+LOKL63Aw8DN4nISKBPnSYyXmVgcgsmDUnkzR93snTnEafjGFNj+cWl/PqjNSS2bMyvRnd2Oo7PcaeI5KrqMVV9CBhNxaCMxoc9MqYLbSPDePiD1eQX22kt07D8ZdZmdh/J55mrexHWyE5j1TV3isiXJ56o6qPAtLqLY7xReEgQz17Ti52H8/nLLDutZRqORdsO8cainUwaksjgdi2djuOTql1EROTvIiIn5v44QVVfrPtYxtsM7RDFLUOTeGPRTrtayzQIuYUlPPzBGpKjwnl0bFen4/ismhyJ5AIzRaQxgIhcLCILq3iN8SG/HtOF5KhwHv5wDQWlNiSK8W5PfrGRfdkF/PXa3nYay4OqXURU9bfAdGC+q3g8CDzqqWDG+4Q1CuSv1/ZmX3YB724sdjqOMWf13ab9vLcsg5+f396uxvKwmpzOGgXcCRwHooD7VfUHTwUz3ql/YnPuOr89P2SWMmtdltNxjPk/DuUV8ciHa+gSE8EDF3Z0Oo7Pq8nprMeB36lqKnAN8J7rEl/jZx64sBOJTQN47OM1HMgpdDqOMf+lqjz60RpyCkuZMr4PIUF2GsvTanI6a6SqLnA9XwuMBZ70VDDjvRoFBfDzXiEUlJTx0IdrKLeZEI2XeHfJbuZuPMCvx3ShS0xTp+P4BXcu8QVAVfdRMQSK8UNtmwTw+KXd+H7LQd60mRCNF9h2MI8/f7GBER2juHVoktNx/IbbRQRAVW1SKj82cVACI7u05n+/3sSGvTlOxzF+rKi0jF/OWElocMXFHzbEe/2pVREx/k1E+Ms1vYgMC+a+6SvsbnbjmGe+3sy6zBz+cnUvm6mwnlkRMbUS1SSEKdf3Yfuh4/xh5nqn4xg/9O3G/by+cAe3DE1idPcYp+P4HSsiptaGdYjintT2vL9sD5+tynQ6jvEjWdmFPPTBarq2acqjY7s4HccvWRExdeKBCzvRP7E5j3+yjh2Hjjsdx/iB0rJy7p+xkqLScl66oa/NEeIQKyKmTgQHBvDChL4EBQp3v72cwpIypyMZH/e3OVtYsuMIT17Rg/atmjgdx29ZETF1JrZZGM9f34dNWbk88dk6p+MYH/btxv38M20bEwbGc1W/OKfj+DUrIqZOXdC5Nb+4oAPvL9vDB8synI5jfFDGkXwefH813do05fc/6+50HL9nRcTUuf93USeGtGvJ7z5bx/q92U7HMT6ksKSMe99dQXm58s+J/awfxAtYETF1LjBAeGFCX5qFNeKut5dzLN9G/DW1p6r87tN1rNmTzd+u601iy3CnIxmsiBgPaRURwj8n9mN/dhH3TV9JmY2vZWrp7Z9288HyPdw/soPdD+JFrIgYj+mb0Jw/juvOD1sP8bdvbFpd477lu47wp8/Xc0HnVjxwYSen45hKvKqIiMivRERFJMq1LCLygoiki8gaEelXad9JIrLV9ZjkXGpzLhMGJjBhYAL/SNvG56v3Oh3HNEB7jxVw19sriG0WxpTxfW1cLC8T5HSAE0QkHhgN7K60eizQ0fUYBPwTGCQiLYDfAymAAstFZKaqHq3f1KY6/nB5N7buz+WhD1aT2LIxveKaOR3JNBD5xaXcOW0ZBcVlvHPHICLDgp2OZE7jTUcizwOPUFEUThgHTNMKi4FmItIGuBiYo6pHXIVjDjCm3hObagkJCuRfN/UnqkkId05bxn6byMpUg6ry8Adr2LAvhxcm9KFTdITTkcwZeEUREZFxQKaqrj5tUyxQ+WaDPa51Z1tvvFRUkxBenZRCbmEpk6ctszvaTZX+/u1Wvly7j8fGdmFkl2in45izqLfTWSIyFzjTJRWPA7+h4lSWJz53MjAZIDo6mrS0NLfeJy8vz+3X+iJ32+OO7kG8uDKbiS/N4Z4+IQSIb5zftp+Pk+qiLRbtLWXqmiKGtQ2iY9lu0tIa7o2rvv6zUW9FRFUvPNN6EekJJAOrpeILJQ5YISIDgUwgvtLuca51mUDqaevTzvK5U4GpACkpKZqamnqm3aqUlpaGu6/1Re62RyrQtO12nvxyIz8VxPDYJV3rOpoj7OfjpNq2xeLth3ljzhIGt2vBf24bRKMgrzhh4jZf/9lw/F9HVdeqamtVTVLVJCpOTfVT1SxgJnCz6yqtwUC2a1re2cBoEWkuIs2pOIqZ7dTfwdTM7cOTuXlIIq98v523F+9yOo7xIukH8vj5W8uJbxHGKxNTGnwB8Qdec3XWWXwFXAKkA/nArQCqekRE/gwsde33J1U94kxEU1MiwhOXdSPzaAFPfLaONpGhjOpq57z93YHcQm59YwnBgcIbtw4ksrFdidUQeF2Zdx2RHHI9V1W9V1Xbq2pPVV1Wab/XVbWD6/Ef5xIbdwS5ho7vERvJve+uYPku+x3An+UUljDp9aUcyi3m1UkDiG/R2OlIppq8rogY/xEeEsR/bhlAm8gwbntjGVv35zodyTigsKSMydMq/v3/dVN/+sQ3czqSqQErIsZRLZuEMO22gYQEBXDz60vIPFbgdCRTj8rKlQffX8Xi7Uf467W9Ob9TK6cjmRqyImIcF9+iMW/eNpC8olImvvoTB3LtZkR/UF6uPPrRGr5am8VvL+3KFX3tVq+GyIqI8Qpd2zTljVsHsD+nkJteXcLR4zZ8vC9TVf70xQY+WL6HX47qyB0j2jkdybjJiojxGv0TW/DqzSnsOHycm19fQk5hidORjAeoKs/O3swbi3Zy54hkHriwo9ORTC1YETFeZWiHKP41sR8b9+Uw6fUl5Foh8SmqyvNztvCPtG3cMCiB31zSFfGRUQv8lRUR43VGdonmpRv6sXZPth2R+BBV5bk5W3jhu3SuT4nnyXE9rID4ACsixiuN6RHDyze6CslrVkgauhMF5MXv0hk/IJ7/vaqnzQviI6yIGK91cfeKQrIuM5uJr/5kne0NlKryv19v+m8B+Z8rrYD4Eisixqtd3D2GqTf3Z3NWLtdP/ZEDNhdJg1JWrvzmk3VM/X47Nw9JtALig6yIGK83sks0b9w6kMyjBVz7yo9kHMl3OpKphpKych58fxXTl+zmntT2/PHy7lZAfJAVEdMgDGnfkrfvGMSx/BKu/uci1u/NdjqSOYfjRaXc/uYyPlu1l4cv7swjY7pYJ7qPsiJiGoy+Cc354K4hBAYI17+ymIXph5yOZM4gu0gZP7Xi3+fpq3py7wUdnI5kPMiKiGlQOkVH8PE9Q4ltFsYt/1nCpysznY5kKtl+MI8nFxew9UAuU2/qz/iBCU5HMh5mRcQ0OG0iw3j/riH0T2zOA++t4m/fbKa8XJ2O5fcWpR/iyn8sorBUmX7nYJsjxk9YETENUmRYMNNuG8R1KXG8+F06v5i+goLiMqdj+a13f9rNza8vIbppCE8MCaNvQnOnI5l6YkXENFiNggJ45upe/OaSLny9LotrX1lkV27Vs+LScp74bB2/+WQtwztG8dHdQ2nV2L5W/In9a5sGTUSYfF57Xr05hV2H8vnZSwv4YetBp2P5hf05hUz492Km/biLO0ck8+rNKUSE2pS2/saKiPEJo7pGM/O+4bSOCGHS60t4eV669ZN40OLth7nsxQVs3JfDixP68vil3QgKtK8Tf2T/6sZnJEeF88k9w7i0V1uenb2ZSf9ZwsHcIqdj+ZSycmXK3C3c8O/FNAkJ4pN7hvGz3m2djmUcZEXE+JTwkCBeGN+H/7myJ0t2HGHs33/g+y12eqsuZGUXcsO/FzNl7lbG9Ynl8/uG0zkmwulYxmFWRIzPERFuGJTAzF8Mp0V4MDe/voTff7bOrt5yk6ry6cpMRj8/n7WZ2fzt2t48f30fmoQEOR3NeAH7KTA+q3NMBDN/MZy/zNrM6wt38P3WQ/z12t70T7TLT6vrcF4Rv/10HV+vy6J/YnP+em1vkqPCnY5lvIgdiRifFhocyBM/68b0OwdTXFrONf9axB9mrrcZE6ugqny4fA8XPjefbzce4NdjuvD+z4dYATH/hx2JGL8wpH1LZj0wgr/O3sybP+5k1ros/nB5dy7uHm0DA55mx6HjPP7JWhZtO0z/xOb871U96RRtfR/mzKyIGL8RERrMH8f14Iq+sTz28Vruens5IzpG8cRl3ehoX5LkFpbw0nfpvL5wB6FBgTx5RQ9uGJhgw7ebc7IiYvxO34TmfH7fcN76cRdT5m5hzN9/YOKgBO4b1ZGoJiFOx6t3pWXlfLB8D3/7ZjOH8oq5tn8cD1/cmdZNQ52OZhoAKyLGLwUHBnDb8GSu6BvLc3M289biXXy4fA+3D0/mjvPa0dQP7rwuL1e+XLuP5+ZsYceh4/RPbM5rkwbQO76Z09FMA2JFxPi1FuGNePKKntwyNJnn52zhhe/SmbZ4F7cMTeKWoUk0a9zI6Yh1rrSsnK/WZfGPeelsysqlU3QTpt7Un4u6Wf+QqTkrIsYAHVo34eUb+3F3ZjZT5m5hytytTP1+OzcMTGDS0CTiWzR2OmKt5ReX8vGKTP79w3Z2Hc6nfatwnr++N5f3jiXQ+j2Mm7ymiIjIfcC9QBnwpao+4lr/GHC7a/39qjrbtX4M8HcgEHhVVZ92JLjxKT1iI3l10gA2ZeXwyvzt/GfRTl5buINRXVpz05AkRnSIanAdzdsP5vH24t18sDyD3MJSesVF8q+J/RjdLabB/V2M9/GKIiIiFwDjgN6qWiQirV3ruwHjge5AW2CuiHRyvexl4CJgD7BURGaq6ob6T298UZeYpjx/fR8evrgz7/60m+lLdjN34xJim4VxRd+2XNk3jg6tmzgd86yy80v4Yu1ePl6RyfJdRwkKEMb2bMOkIYn0T2xup61MnfGKIgLcDTytqkUAqnrAtX4cMMO1foeIpAMDXdvSVXU7gIjMcO1rRcTUqbbNwnjo4s7cN6oDs9Zl8fGKTP6Zto2X522jS0wEo7vHcHH3aFSdHzH4QG4hczccYPb6LBZtO0RJmdKhdRMeGdOZa/rF2dVWxiO8pYh0AkaIyFNAIfCQqi4FYoHFlfbb41oHkHHa+kH1EdT4p5CgQMb1iWVcn1gO5BQyc/Vevlm/nxe/28oL326lWYhwwYFVDOsQxYCk5iS0aOzx3/aP5RezcvcxFqYfYkH6ITZl5QKQ0KIxtwxN4vLesfSIbWpHHcaj6q2IiMhcIOYMmx535WgBDAYGAO+LSLs6+tzJwGSA6Oho0tLS3HqfvLw8t1/ri/y9PToAHbpATnJjVh0sZXVWEXPXZ/LJykwAwoMhOTKQhIgA2oQLbZoE0LpxABHB1PhLvaBUOVKg7D1ezr7j5WTmlrMjp5wD+RVHP0EB0LFZAFd3DKZP6yDimoDIAQ6nH2B+el3/zavm7z8bp/P19qi3IqKqF55tm4jcDXysFecElohIORAFZALxlXaNc63jHOtP/9ypwFSAlJQUTU1NdSt/Wloa7r7WF1l7nHQ5Fe1x3nnns+VALit3H2N1xjFWZRxjzu48SspOnupqFBhA66YhtGwSQpOQQJqEBBESFMiJulJapuQWlXK8qJRj+cXszykir6j0lM+LbRZGv+RIesc3o3d8JP0SmhMaHFiPf+Nzs5+NU/l6e3jL6axPgQuAea6O80bAIWAm8K6IPEdFx3pHYAkgQEcRSaaieIwHbnAgtzH/FRAgdIlpSpeYpkwYmABU3JORcbSA9AN57DmaT1ZOIfuzCzmSX8LxolIO5eZTVHpyiPrAAKFJaDBNQgLpFB3BiI6tiIkMpU1kKO1bNSE5KpxwG4LdeBFv+Wl8HXhdRNYBxcAk11HJehF5n4oO81LgXlUtAxCRXwCzqbjE93VVXe9MdGPOLigwgOSocBv91vgsrygiqloMTDzLtqeAp86w/ivgKw9HM8YYcw42n4gxxhi3WRExxhjjNisixhhj3GZFxBhjjNusiBhjjHGbFRFjjDFusyJijDHGbeINo4/WFxE5COxy8+VRVNxFbypYe5zK2uMka4tT+UJ7JKpqqzNt8KsiUhsiskxVU5zO4S2sPU5l7XGStcWpfL097HSWMcYYt1kRMcYY4zYrItU31ekAXsba41TWHidZW5zKp9vD+kSMMca4zY5EjDHGuM2KiDHGGLdZEakGERkjIptFJF1EHnU6T30TkddF5IBr0rAT61qIyBwR2er6s7mTGeuLiMSLyDwR2SAi60Xkl671/toeoSKyRERWu9rjj671ySLyk+v/zHsi0sjprPVFRAJFZKWIfOFa9um2sCJSBREJBF4GxgLdgAki0s3ZVPXuDWDMaeseBb5V1Y7At65lf1AK/EpVuwGDgXtdPw/+2h5FwEhV7Q30AcaIyGDgGeB5Ve0AHAVudy5ivfslsLHSsk+3hRWRqg0E0lV1u2sGxhnAOIcz1StV/R44ctrqccCbrudvAlfUZyanqOo+VV3hep5LxZdFLP7bHqqqea7FYNdDgZHAh671ftMeIhIHXAq86loWfLwtrIhULRbIqLS8x7XO30Wr6j7X8ywg2skwThCRJKAv8BN+3B6u0zergAPAHGAbcExVS127+NP/mSnAI0C5a7klPt4WVkRMrWnFdeJ+da24iDQBPgIeUNWcytv8rT1UtUxV+wBxVBy5d3E2kTNE5DLggKoudzpLfQpyOkADkAnEV1qOc63zd/tFpI2q7hORNlT8FuoXRCSYigLyjqp+7Frtt+1xgqoeE5F5wBCgmYgEuX4D95f/M8OAy0XkEiAUaAr8HR9vCzsSqdpSoKPrCotGwHhgpsOZvMFMYJLr+STgMwez1BvXOe7XgI2q+lylTf7aHq1EpJnreRhwERX9RPOAa1y7+UV7qOpjqhqnqklUfE98p6o34uNtYXesV4PrN4spQCDwuqo+5Wyi+iUi04FUKoa03g/8HvgUeB9IoGJ4/etU9fTOd58jIsOBH4C1nDzv/Rsq+kX8sT16UdFZHEjFL6Xvq+qfRKQdFRehtABWAhNVtci5pPVLRFKBh1T1Ml9vCysixhhj3Gans4wxxrjNiogxxhi3WRExxhjjNisixhhj3GZFxBhjjNusiBhjjHGbFRFj6oGIxInI9WdYnyQiBa6xp8722jARWSUixSIS5dGgxtSQFRFj6scooN9Ztm1zjT11Rqpa4Nq+1wO5jKkVKyLGeJjrLvfngGtcRxTtzrFvuIh86Zrkad2Zjl6M8SY2AKMxHqaqC0RkKRXDYKyrYvcxwF5VvRRARCI9HtCYWrAjEWPqR2dgUzX2WwtcJCLPiMgIVc32cC5jasWKiDEe5uoMz640MdFZqeoWKvpO1gJPisgTns5nTG3Y6SxjPC+JanaKi0hb4Iiqvi0ix4A7PJjLmFqzImKM520CokRkHTBZVRedY9+ewLMiUg6UAHfXR0Bj3GVFxBgPU9U8KqaNrc6+s4HZnk1kTN2xPhFjnFUGRFbnZkMgmJMTYRnjFWxSKmOMMW6zIxFjjDFusyJijDHGbVZEjDHGuM2KiDHGGLdZETHGGOM2KyLGGGPcZkXEGGOM26yIGGOMcdv/B4TupUORFUeiAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"setup_typeset() #active colab pretty print\n",
"def x(x0,t0,v0,a):\n",
" return np.poly1d( [0.5*a,v0-a*t0,x0-v0*t0+0.5*a*t0**2],\n",
" variable='t' )\n",
"x1=x(0,0,0,6)\n",
"x2=x(0,10,10,10)\n",
"display(Latex('$x_1(t)=$'))\n",
"print(x1) \n",
"display(Latex('$x_2(t)=$'))\n",
"print(x2)\n",
"t=np.linspace(0,45,100)\n",
"plt.plot(t,x2(t)-x1(t))\n",
"plt.xlabel('$t$ [s]')\n",
"plt.ylabel('$x_2(t)-x_1(t)$ [m]')\n",
"plt.grid()\n",
"#plt.plot(t,(x2-x1)(t))\n",
"#plt.grid()\n",
"#Physica solution is the one after 10s\n",
"Latex(r'meeting time $t_{\\rm end}=$ %g s; meeting point $x_{\\rm end}=$ %g m' \n",
" %( (x2-x1).r[0] , x2( (x2-x1).r[0] ) )) \n",
"#plt.plot(t,x2(t),'ro')\n",
"#plt.plot(t,x2(t),'k-')\n",
"#plt.xlim(10,15)\n",
"#plt.ylim(0,200) "
]
},
{
"cell_type": "code",
"execution_count": 77,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
"10 x - 90\n"
]
}
],
"source": [
"v=x_2.deriv()\n",
"print(v)"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
"10\n"
]
}
],
"source": [
"print(v.deriv())"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "vWJZfq_YEZWe",
"tags": []
},
"source": [
"## Linear Interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "wGj1U7PmEZWe"
},
"source": [
"When we have a set of discrete points of the form $(x_i, y_i)$ for $1\\leq i \\leq N$, the most natural way to obtain (approximate) any intermediate value is assuming points connected by lines. Let's assume a set of points $(x_i, y_i)$ such that $y_i = f(x_i)$ for an unknown function $f(x)$, if we want to approximate the value $f(x)$ for $x_i\\leq x \\leq x_{i+1}$, we construct an equation of a line passing through $(x_i,y_i)$ and $(x_{i+1},y_{i+1})$.\n",
"\n",
"The linear equation is\n",
"\n",
"$$y=mx+b$$\n",
"\n",
"where\n",
"\n",
"$$m=\\frac{y_{i+1}-y_i}{x_{i+1}-x_i} $$\n",
"\n",
"and $b$ is obtained by evaluating with either $(x_i,y_i)$ or $(x_{i+1},y_{i+1})$\n",
"\n",
"$$y=\\frac{y_{i+1}-y_i}{x_{i+1}-x_i}x+b$$\n",
"\n",
"$$b=y_i-\\frac{y_{i+1}-y_i}{x_{i+1}-x_i}x_i$$\n",
"\n",
"\n",
"\\begin{align}\n",
"%$$ \n",
"f(x)\\approx y = &\\frac{y_{i+1}-y_i}{x_{i+1}-x_i}(x-x_i) + y_i \\\\\n",
"=&\\frac{y_{i+1}-y_i}{x_{i+1}-x_i}x+\\left[y_i-\\frac{y_{i+1}-y_i}{x_{i+1}-x_i}x_i\\right] \\\\\n",
"%$$\n",
"\\end{align}\n",
"\n",
"and this can be applied for any $x$ such that $x_0\\leq x \\leq x_N$ and where it has been assumed an ordered set $\\left\\{x_i\\right\\}_i$."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "WsnL9PoDEZWi",
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"### Steps LI"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "w4MSYnO6EZWi"
},
"source": [
"Once defined the mathematical basis behind linear interpolation, we proceed to establish the algorithmic steps for an implementation.\n",
"\n",
"1. Establish the dataset you want to interpolate, i.e. you must provide a set of the form $(x_i,y_i)$.\n",
"2. Give the value $x$ where you want to approximate the value $f(x)$.\n",
"3. Find the interval $[x_i, x_{i+1}]$ in which $x$ is embedded.\n",
"4. Use the above expression in order to find $y=f(x)$."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "0Gmn6NAdjz1W"
},
"source": [
"To make the linear interpolation we will use"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "bzCHHksPEZWn"
},
"outputs": [],
"source": [
"import scipy as sp\n",
"sp.interpolate.interp1d?"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "epbBukAVjz1Z"
},
"source": [
"The option `kind` specifies the kind of interpolation: \n",
"* `'linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic'`,\n",
" where `'zero', 'slinear', 'quadratic' and 'cubic'` refer to a spline\n",
" interpolation of zeroth, first, second or third order) \n",
"* or as an\n",
" integer specifying the order of the spline interpolator to use.\n",
"\n",
"Default is `'linear'` correponding to integer 1."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "OZXQtU4WEZWj",
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"### Example 1"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "Jb7eChUUEZWk"
},
"source": [
"Sample the function $f(x) = \\sin(x)$ between $0$ and $10$ using $N=10$ points (9 intervals). Plot both, the interpolation and the original function."
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "5G1t56pnjz1a"
},
"outputs": [],
"source": [
"import scipy as sp\n",
"from scipy import interpolate "
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "A0mxSeRjk1T6"
},
"source": [
"`sp` reemplaza completamente a numpy con `np`"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "sB0lNqKIEZWn"
},
"source": [
"Interpolation with 9 equal intervals"
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 85
},
"colab_type": "code",
"id": "BaCvQYBoEZWp",
"outputId": "00b9b521-818a-4f8b-c652-df4ec980a4c9"
},
"outputs": [],
"source": [
"n_points = 10\n",
"x=np.linspace(0, 2*np.pi, n_points)\n",
"f=interpolate.interp1d( x,np.sin(x),kind='linear' )"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "IFSu7iE6jz1g"
},
"source": [
"Plotting the results adding the real function with enough points"
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 283
},
"colab_type": "code",
"id": "3Ggn5WacEZWq",
"outputId": "c33ed8a9-bd0c-4418-85e9-45e32cb2ab0d"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuUAAAGKCAYAAACrcD/sAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAACA/ElEQVR4nOzde3yO9R/H8dd3M5uZw7CR05wPcxpWFEJUQiWHcgw5VNKBhOikKKfQQSSJQnT4VYRSWEUqcj4f55jzaUeb7fv7Y7M2Rg7brnvb+/l43I/tOt3X+97X4ePyub6XsdYiIiIiIiLOcXM6gIiIiIhIdqeiXERERETEYSrKRUREREQcpqJcRERERMRhKspFRERERBymolxERERExGEqykVEREREHKaiXERERETEYSrKRSTbM8ZsNsY0cjoHuFaW/2KMCTXGNL3BYzP8cxpjKhpj1hljwowxz1znsdeV1xjzljHmuWvY7y9jTJXrySIiWZPREz1FJLswxoQCPa21PzudJS04/Xmu9fxO50yW42PgnLW2Xzqfxw9YB5Sz1kb9x74PA49Ya9ukZyYRcX26Ui4i4gBjTI7sfH6HBACbM+A83YCF/1WQJ5oHNDbGFEnfSCLi6lSUi0i2d2kbRuLyAGPMBmPMWWPMXGOMV+K2osaYr40xx40xey9tgzDGDDbG7E5skdhijHnokvcdZIzZAESkVhgnz/IfOT4DSgLzjTHhxpiB15DtsvMnrnsxMetpY8wnyc5R2RgTYow5k9i+8cBVfoapfu7UcqbyOa96nqv9HFLJkep7GWOWAo2B9xNzVEjl2EHGmEOJn2G7MabJ9Y5LovuAX5IdO9oY822y5THGmCXGmJzW2mjgb+DeK/1sRSR7UFEuIpK6h4FmQGmgOtDNGOMGzAfWA8WAJsBzxpjkBdVuoAGQDxgGzDTG3JJsewegBZDfWnvhRnIAWGu7APuB+621PsDYa8h2pfN3IqEoLAtUAF4yxngkvt9iwB94GphljKl4hZypfu5Lc1prRyc/6DrOk+rP4Vrfy1p7F/Ab0Dcxx45Ljq0I9AVutdbmSfx5hF7hs/5XnmrA9mTLo0i4Gl7TGPNE4nGtrbUxidu3AjWuci4RyQZUlIuIpO5da+1ha+0pEgq9IOBWwM9a+7q1NsZauwf4CGh/8SBr7ZeJx8Vba+cCO4HbLnnfA9fY2nClHKn5z2xXOf/7ietOASNIKNzrAj7AyMT3Wwp8n7jtMtfwua/kWs9zLT+H68p8iTjAEwg0xnhYa0Ottbuvsv/V8uQHwi4uWGtPAuOBGcCLQHNr7dlk+4clHiMi2ZiKchGR1B1J9n0kCcVeAFA0sTXijDHmDDAEKHxxR2PMoyZhho+L26sChZK914E0yJGa/8x2lfMnX7cPKJr4OmCtjb9kW7HUTn4Nn/tKrvU81/JzuK7MyVlrdwHPAa8Bx4wxc4wxRa9yyNXynAbyXLL/WhKuoL9orb10DPIAZ/4ro4hkbSrKRUSu3QFgr7U2f7JXHmttcwBjTAAJV6f7AgWttfmBTYBJ9h5pOeVV8ve6arb/OH+JZN+XBA4nvkoktuwk33bo0oOv4XNf7TNf83muwU29l7V2trW2Pgn/wLEktJ3ciA0ktAEBYIypBkwi4Ur5Y6nsX5mEtiMRycZUlItIduNhjPFK9rqeWUj+AsISbwjMZYxxN8ZUNcbcmrg9NwnF3HEAY0x3Eq4Yp5ejQJlrzHY1TxljihtjCgBDgbnAnyRcAR5ojPEwCXN03w/MSeX4//rcyXNe6nrO819u+L1MwhzmdxljPIFoIAqI/4/DrmQh0DDxfYuR0N7yBNAHqGaSzXeeeINobeCnGzyXiGQRKspFJLtZSELBdfH12rUeaK2NA1qS0D+8FzgBTCXh5kastVuAt4GVJBSi1YAVaZb8cm+RcFPmGaDf1bL9h9kk3By5h4QbNocn3oR4PwkziZwAPgAetdZuu/Tga/jcSTmNMQMuOfaaz/NfbvK9PIGRiccdIeFG0RevN0OiT4Hmxph8JPx6G2etnWetjQTGkNC3f9H9QIi19vANnktEsgg9PEhEJBszLvJgn6zGGPMmcMxaO+E/9vsT6GGt3ZQhwUTEZWXHh0eIiIikK2vtkGvcr056ZxGRzEHtKyIiIiIiDlP7ioiIiIiIw3SlXERERETEYSrKRUREREQcphs9gUKFCtlSpUpl+HkjIiLInTt3hp9X/qUxcJ7GwHkaA+dpDJynMXBedhiDv//++4S11i+1bSrKgVKlSrF69eoMP29ISAiNGjXK8PPKvzQGztMYOE9j4DyNgfM0Bs7LDmNgjNl3pW1qXxERERERcZiKchERERERh6koFxERERFxmHrKRURERNJIbGwsBw8eJDo62ukomU6+fPnYunWr0zHShJeXF8WLF8fDw+Oaj1FRLiIiIpJGDh48SJ48eShVqhTGGKfjZCphYWHkyZPH6Rg3zVrLyZMnOXjwIKVLl77m49S+IiIiIpJGoqOjKViwoArybMwYQ8GCBa/7f0tUlIuIiIikIRXkciO/BlSUi4iIiGQh7u7uBAUFUaVKFWrUqMHbb79NfHz8VY8JDQ1l9uzZGZRQUqOiXERERMQps2ZBqVLg5pbwddasm37LXLlysW7dOjZv3sxPP/3EokWLGDZs2FWPUVHuPJcsyo0x04wxx4wxm66w3Rhj3jXG7DLGbDDG1Eq2rasxZmfiq2vGpRYRERG5DrNmQe/esG8fWJvwtXfvNCnML/L392fKlCm8//77WGsJDQ2lQYMG1KpVi1q1avH7778DMHjwYH777TeCgoIYP378FfeT9OOqs69MB94HPr3C9vuA8omvOsAkoI4xpgDwKhAMWOBvY8w8a+3pdE8sIiIicj2GDoXIyJTrIiMT1nfqlGanKVOmDHFxcRw7dgx/f39++uknvLy82LlzJx06dGD16tWMHDmSsWPH8v333yfGiEx1P0k/LlmUW2t/NcaUusouDwKfWmst8IcxJr8x5hagEfCTtfYUgDHmJ6AZ8Hk6RxZxTLy1HImJ4dyFC8RYS0x8PDHWcj7xq6cxFPDwoKCHBwVy5CCXu7vTkUVEBGD//utbnwZiY2Pp27cv69atw93dnR07dtzUfpJ2XLIovwbFgAPJlg8mrrvS+ssYY3oDvQEKFy5MSEhIugS9mvDwcEfOK/9yxTHw//lnykydiuexY5z392dPz54cadqUg8B24DBwBDia+DoOxF7H+3sCeYGCQCmgdLKvhYCMnjPAFccgu9EYOE9j4Ly0GoN8+fIRFhZ2TfvmLl4ctwMHLlsfX7w4Edf4HleSPMPevXtxc3MjV65cvPXWW/j6+rJ8+XLi4+Px8/MjLCyMyMhILly4kHTcyJEjU90vPcXFxaX7OTJSdHT0df2ayqxF+U2z1k4BpgAEBwfbRo0aZXiGkJAQnDiv/MvVxsDOnMknY8bQ09OTvMHB+AQGcip/ftbFx3PaLeEWEAPckjMnAV5e3OnlRYCnJwFeXuTPkYOcbm54GkNONzdyJn49Hx/PydhYTl64wKnY2KTv90dHszYigh9i/y3p87m7U8PHh7t8fWnq68ttefLg4Za+t5642hhkRxoD52kMnJdWY7B169ZrfwDOW28l9JAnb2Hx9sbtrbdu+iE6F48/fvw4AwYM4OmnnyZv3rxER0cTEBBAvnz5+OSTT4iLiyNPnjwULlyYqKiopOOutF96yioPD7rIy8uLmjVrXvP+mbUoPwSUSLZcPHHdIRJaWJKvD8mwVCI3Yd3x43T74QfWjx0LVaok3IkfHw979+L7ww80ueUWWlWowCP16uGXP3+anfdETAybIyPZFBHB5ogIVoWFMSw0lNdCQ/Fxd6dR/vw09fXlbl9fAnPnTrPziohkexf7xocOTWhZKVkSRoy46X7yqKgogoKCiI2NJUeOHHTp0oX+/fsD0KdPH9q0acOnn35Ks2bNyJ3453r16tVxd3enRo0adOvW7Yr7SfrJrEX5PKCvMWYOCTd6nrXW/mOM+RF40xjjm7jfPcCLToUUuRprLWvDw/nmxAnmHDjArvh46NkTdu6ETz+F9eth+3aIiuI0sCTx9Zy7O8HBwTRp0oQmTZpwxx134OXldcM5CuXMScOcOWmYrNA/FRtLyJkz/Hz6ND+fPs33J08CUC13bjoXLkxHf3+K38Q5RUQkUadOaXpTJyS0gVxJ+fLl2bBhQ9LyqFGjAPDw8GDp0qUp9k1tP0k/LlmUG2M+J+GKdyFjzEESZlTxALDWTgYWAs2BXUAk0D1x2yljzBvAqsS3ev3iTZ8iruLchQt8dvQokw8fZlNEBMZa2LABfvsNli+Ho0fpQcJ/8ywB/gAuJDs+Li6OP//8kz///JM333wTT09P6tWrx1133UWTJk0IDg4mR46b+61dwMOD1n5+tPbzA2BfdDQLTp5k5tGjDNqzh8F79tA4f346Fy5MGz8/8t7k+URERLI7l/yb1Frb4T+2W+CpK2ybBkxLj1wiN2NNWBiTDx9m9tGjRMTHE+TtTfAvv7B63Dg4dw4AHy8vpuTMSYeYGABeA8Jz5WJ5374ssZalS5eydu1aEn4LJDh//jxLly5l6dKlvPTSS+TNm5eGDRsmFelVq1a96Uc+B3h50adYMfoUK8buqChmHj3KzKNHeWz7dvrs3MmjhQvTv0QJKnp739R5REREsiuXLMpFsoo4a/nq+HHePnCAVWFh5HJzo6O/P43DwnitY0d27dqVtG+NGjX44osvqLBqVYr+Qp8RI2jWqRPNEvc7efIkISEhLF26lCVLlrB9+/YU5zx37hzz589n/vz5QMKDIxo3bpzU7lKmTJmb+kxlc+Xi1VKleCUggL/Cwvj4n3+YceQIH/3zDw8ULMiAEiWoly/fTf9DQEREJDtRUS6SDuKt5cvjx3k9NJQtkZFUzJWLd8uVo3Phwsz9+GN6PPcc58+fT9r/8ccfZ/z48eTKlQsqVLhqf2HBggVp06YNbdq0AeDQoUNJBfqSJUs4ePBgiv2PHTvG3LlzmTt3LgClSpVKuop+1113UaRIkRv6jMYY6uTNS528eRleujQTDx1i4qFDfHfyJHXy5GFAiRI85OeHu4pzERGR/6SiXCQNxSdeGR+WWIwHenszJzCQtn5+RISF0btLl6TiGMDHx4cpU6bQocNVO7auqlixYnTp0oUuXbpgrWXXrl0sWbIkqaXlZOJNmheFhoYybdo0pk1L6PIKDAxMuoresGFD8t/AzC7+OXMyrHRpBpUsyfQjRxh34ADttmyhWu7cjCpThmYFCujKuYiIyFWoKBdJA9ZavjtxgqF797IlMpLKyYpxd2NYu3YtDz/8cOrtKhUqpFkOYwzly5enfPnyPPHEE8THx7Nhw4akIv2XX34hIiIixTFbtmxhy5YtvPfee7i5uVG7du2kK+n16tXD+zr6xL3d3elTrBiPFy3KV8ePM3TPHppv3Mhd+fMzumxZameh+WdFRETSUvo+FUQkG9gWEcG9Gzbw0ObNxANzAgPZeOutPOLvjxswefJkbr/99hQF+eOPP87KlSvTtCBPjZubG0FBQTz//PMsWLCA06dPs3z5coYNG0bDhg3JmTNniv3j4+NZtWoVo0aN4p577sHX15dGjRrxxhtv8PvvvxMbe23PDnU3hkf8/dly2228W64cGyIiCP77bzpt2cLeqKj0+KgiIpLIx8cHgMOHD9O2bdsMPfe8efMYOXLkVfcJDQ1l9uzZGZKnUaNGrF69+qr7TJgwgchkD3Bq3rw5Z86cSedkl1NRLnKDzl24wIBdu6i2ejV/nTvHO+XKsTE4mEf8/XE3hnPnztG+fXuefPLJpP5xHx8fZs+ezeTJkxP6xzOYh4cH9erV45VXXiEkJITTp0+zePFiBg0aRHBw8GUtJjExMfzyyy+88sor1KtXjwIFCtCiRQvGjRvHunXriI+Pv+r5crq58XTx4uyuU4ehJUvyzYkTVPrrLwbv3k3kVebRFRGRm1e0aFG++uqrdD3HhQsXUiw/8MADDB48+KrH3EhRful50tKlRfnChQtvqJXzZqkoF7lO1lo+O3KEin/9xdsHD9K1cGF21KnDM8WLkyPxkfRr166ldu3afPHFF0nHVa9enb///vum+sfTmre3N3fffTcjR45k1apVnDx5km+++Ya+ffsSGBh42f7h4eEsXLiQ559/npo1a+Lv70+7du2YPHkyO3fuTDFVY3J5c+RgeJky7KxThw7+/ow6cIDqq1bx8yk9RkBEJL2EhoZStWpVAKZPn07r1q1p1qwZ5cuXZ+DAgUn7LV68mNtvv51atWrRrl07wsPDAXj99de59dZbqVq1Kr179076M75Ro0Y899xzBAcH884776Q45/Tp0+nbty8A3bp145lnnuGOO+6gTJkySf9AGDx4ML/99htBQUGMHz+euLg4XnjhBRo2bEj16tX58MMPAQgJCaFBgwY88MADBAYGEhoaSqVKlejUqROVK1embdu2ScX0kiVLqFmzJtWqVeOxxx5LMZnCRU8++STBwcFUqVKFV199FYB3332Xw4cP07hxYxo3bgwkTIhw4sQJAMaNG0fVqlWpWrUqEyZMSPq5Vq5cmV69elGlShXuueceotLif4Gttdn+Vbt2beuEZcuWOXJe+df1jsGuyEh755o1lmXL7K2rV9s/z55NsT0+Pt5+8MEH1tPT0wJJr969e9vIyMg0TJ4xDh8+bGfNmmUfe+wxGxAQkOIzpfYqUaKE7datm/3000/toUOHrvi+S0+dsuX/+MOybJm9d9kyeyImJgM/lVxKfxY5T2PgvLQagy1btiR9/19/Zt7M62py585trbV27969tkqVKtZaaz/55BNbunRpe+bMGRsVFWVLlixp9+/fb48fP24bNGhgw8PDrbXWjhw50g4bNsxaa+3JkyeT3rNz58523rx51lprGzZsaJ988slUz/3JJ5/Yp556ylprbdeuXW3btm1tXFyc3bx5sy1btqy1NuFn3aJFi6RjPvzwQ/vGG2/Yc+fO2ejoaFu7dm27Z88eu2zZMuvt7W337NmT9HkAu3z5cmuttd27d7djxoyxUVFRtnjx4nb79u3WWmu7dOlix48fn5R11apVKT7PhQsXbMOGDe369euttdYGBATY48ePJ+W5uLx69WpbtWpVGx4ebsPCwmxgYKBds2aN3bt3r3V3d7dr16611lrbrl07+9lnn132s0j+a+EiYLW9Qj2qK+Ui18Bay6RDh6i+ahXrwsOZUqECf9SqxW158ybtc7FdpU+fPpe1q3z44YeOtKvcrFtuuYWOHTvy8ccfs3fvXnbt2sWUKVN45JFH8Et82mdyBw4cYPr06Tz66KMUK1aMypUr89RTT/G///2PU8muijf29WV9cDBDSpbkZ6DSX38x6+jRK15pFxGRm9ekSRPy5cuHl5cXgYGB7Nu3jz/++IMtW7ZQr149goKCmDFjBvv27QNg2bJl1KlTh2rVqrF06VI2b96c9F6PPPLINZ2zVatWuLm5ERgYyNGjR1PdZ/HixXz66afUq1ePOnXqcPLkSXbu3AnAbbfdRunSpZP2LVGiBPXq1QOgc+fOLF++nO3bt1O6dOmk+7S6du3Kr7/+etl5vvjiC2rVqkXNmjXZvHkzW7ZsuWr25cuX89BDD5E7d258fHxo3bo1v/32GwClS5cmKCgIgNq1axMaGnpNP4+r0ewrIv/hQHQ0PbZv56fTp7nb15ePK1akhJdXin0yanYVJxljKFu2LGXLlqVXr15Ya9m0aVPS/Oi//PILYWFhKY7Ztm0b27Zt44MPPsAYQ82aNZNmdmnQoAEjypSh7P79TPHyovPWrXx+9CjTKlXC/5IbUEVE5OZ5enomfe/u7s6FCxew1nL33Xfz+eefp9g3OjqaPn36sHr1akqUKMFrr71GdHR00vbcuXNf9zmvdOHFWst7773HHXfcQZ5ks3SFhIRcdp5L73261ul29+7dy9ixY1m1ahW+vr5069Ytxee5Xpf+LNOifUVXykWuwFrLjCNHqLpqFSvOnuWD8uX5sXr1FAW5tZZJkyZRt25dR2ZXcZIxhmrVqvHcc88xf/58Tp06xcqVKxk+fDh33XVXij+wIOFntWbNGsaOHct9992Hr68vDRo04Nfp0xkZEcHYUqX4+fRpqq9axWL1motIFnClNoW0eKWVunXrsmLFiqS/wyIiItixY0dSwVqoUCHCw8PT9IbRPHnypLiIc++99zJp0qSkGb527Nhx2fS9F+3fv5+VK1cCMHv2bOrXr0/FihUJDQ1N+gyfffYZDRs2THHcuXPnyJ07N/ny5ePo0aMsWrToinkuatCgAd9++y2RkZFERETwzTff0KBBg5v78FeholwkFadjY2m9eTPdtm2juo8PG269lSeLFUvxL/Lk7SoxMTGA87OrOClHjhzUrVuXoUOHsmTJEk6fPs3PP//MkCFDqFOnDm5uKf+4iY2NZfny5cyYMYPGDRvySpUq1J46Fc6e5d4NG+i/cycxs2ZBqVLg5pbwddYsRz6biEhW5efnx/Tp0+nQoQPVq1fn9ttvZ9u2beTPn59evXpRtWpV7r33Xm699dY0O2f16tVxd3enRo0ajB8/np49exIYGEiDBg2oWrUqjz/++BVnW6lYsSITJ06kcuXKnD59mieffBIvLy8++eQT2rVrR7Vq1XBzc+OJJ55IcVyNGjWoWbMmlSpVomPHjkktMAC9e/emWbNmSTd6XlSrVi26devGbbfdRp06dejZsyc1a9ZMs5/DpYx6OCE4ONj+1xyW6SEkJIRGjRpl+HnlX6mNwepz52i3ZQuHzp/nzdKl6VeixGWPis8O7Spp7ezZs/z6669J7S6bNm1KfcecOaFPH3jwQfLv2MHTb7xBt4MHKQPg7Q1TpkCnThkZPcvTn0XO0xg4L63GYOvWrVSuXPnmA2VDYWFhKdpXLhUaGkrLli2v/PeHi0nt14Ix5m9rbXBq++tKuUgiay2TDx2i3tq1xFnLbzVrMqBkyRQFubXW0YcBZWb58uXj/vvvZ8KECWzcuJEjR47w+eef06JFC8qUKfPvjjExMGECvPQSZwoX5o0pUyh37728CxAZCUOHOvQJRERE0o9u9BQBwi9c4IkdO5h17BjNChRgZuXKFPTwSLHPuXPn6NWrV4q5x318fJgyZYpLzT2eWRQuXJj27dtTpEgRGjVqRGhoKEuWLGHp0qUsWbKEoytWwPbtMGQIdvBgnq1UCf/336f9/v1ORxcREQeUKlUq01wlvxEqyiXb2xoRQdvNm9kaGckbpUoxJCAAN7WrZLhSpUrRo0cPevTogbWWrcWKseSff/h4wADW9+oF7dvTsXRp8kycSAunw4qIiKQxta9ItvY7cNuaNRyPjWVx9eq8VKpUioL84uwqalfJWMYYAseM4Wlvb5bGx1Pxww9h+HBspUo8MGIE3ySbK1dERCQrUFEu2ZK1lnEHDvASUMnbm7XBwTQtUCDFPld6GNDnn3+eLWdXyXCdOsGUKRQICGAh4LdsGTz9NPFxcbQ5dIjJyf6RJCIiktmpKJdsJzY+nid27OD53bupD/wSFESxS+bUXrt2LbVr107RP16jRg3+/vtv2rdvn8GJs7FOnSA0lDLWMv/338l18CA8+SR2yxaePHiQ57ZtI04zSImISBagolyylTOxsTTfuJEp//zD4JIleQ3wdndP2n6lhwE98cQT/PHHH2pXcVCdOnWYPXs25uxZGDAA/vc/3jlyhHabNxMVF+d0PBGRLKNUqVKcOHHisvVffvkllStXvmw+75tx5swZPvjgg6Tlw4cP07Zt2zR7/8xERblkG7ujorh97Vp+OXOGaRUr8laZMil+A1zpYUCff/45kyZNwivZkzzFGa1atWL8+PEQFwfvvQfvv883x49z9/r1nEp8EpyIiCSw1hIfH59m7/fxxx/z0UcfsWzZsjR7z0uL8qJFi6bp00MzExXlki38ee4cddes4VhMDD/VqEH3W25JsX3t2rXUqlXrsnaVNWvWqF3FxTz77LM8++yzCQtffw3DhvHnmTPUW7uWfYmPhRYRya5CQ0OpWLEijz76KFWrVuXAgQOMGTOGW2+9lerVq/Pqq68m7duqVStq165NlSpVmDJlylXf9/XXX2f58uX06NGDF154genTp9O3b9+k7S1btiQkJARIuKA1dOhQatSoQd26dTl69CgAR48e5aGHHqJGjRrUqFGD33//ncGDB7N7926CgoJ46aWXCA0NpWrVqgBER0fTvXt3qlWrRs2aNZP+MTB9+nRat25Ns2bNKF++PAMHDkzLH6FjNCWiZHk/nzpFq02bKJwzJz9Ur055b++kbRfbVZ577rmkq+OQ0K4yfvx4XR13UW+//Tb79u3j22+/hV9+Ie755zkwbhy3r1nDwmrVCLrKE+FERDLKczt3si48PE3fM8jHhwnly191n507dzJjxgzq1q3L4sWL2blzJ3/99RfWWh544AF+/fVX7rzzTqZNm0aBAgWIiori1ltvpU2bNhQsWDDV93zllVdYunQpY8eOJTg4mOnTp1/x/BEREdStW5cRI0YwcOBAPvroI1566SWeeeYZGjZsyDfffENcXBzh4eGMHDmSTZs2sW7dOsLCwjh58mTS+0ycOBFjDBs3bmTbtm3cc8897NixA4B169axdu1aPD09qVixIk8//TQlSpS4/h+oC9GVcsnSvj5+nBYbN1I2Vy5W1KyZoiA/d+4cr7/+utpVMiF3d3dmzZrFbbfdBoBdt464Pn2Ii43lznXr+OnUKYcTiog4JyAggLp16wKwePFiFi9eTM2aNalVqxbbtm1j586dALz77rtJV7MPHDiQtP5m5cyZk5YtWwJQu3ZtQkNDAVi6dClPPvkkkPDneL58+a76PsuXL6dz584AVKpUiYCAgKSivEmTJuTLlw8vLy8CAwPZt29fmmR3kq6US5Y19fBhHt+xg7p58/J9tWr4JntC59q1a2nXrh27d+9OWqeHAWUu3t7ezJ8/n7p167J3716it23D58knKfrRRzTfuJHZlSvTzt/f6Zgiko391xXt9JI7d+6k7621vPjiizz++OMp9gkJCeHnn39m5cqVeHt706hRI6KvowUwR44cKfrVkx/r4eGBSXzmh7u7OxcuXLjRj3JFnslmTUuvc2Q0XSmXLGn0/v302rGDewoUYHGNGkkFubWWDz74gLp166YoyB9//HHNrpIJ+fv7s2jRInx9fQE4sXUr9umnqZ0rF+23bGHmkSMOJxQRcda9997LtGnTCE9sozl06BDHjh3j7Nmz+Pr64u3tzbZt2/jjjz+u631LlSrFunXriI+P58CBA/z111//eUyTJk2YNGkSAHFxcZw9e5Y8efIQFhaW6v4NGjRg1qxZAOzYsYP9+/dTsWLF68qZmagolyzFWsug3bsZtGcP7f39+a5qVXInTnl49uxZHnnkEZ566qmkdpVcuXIlPQxI7SqZU8WKFfnuu+/ImTMnADvWrsXjpZe4M29eHt22jY8OH3Y4oYiIc+655x46duzI7bffTrVq1Wjbti1hYWE0a9aMCxcuULlyZQYPHpzU7nKt6tWrR+nSpQkMDOSZZ56hVq1a/3nMO++8w7Jly6hWrRq1a9dmy5YtFCxYkHr16lG1alVeeumlFPv36dOH+Ph4qlWrxiOPPML06dNTXCHPaozVgzcIDg62q1evzvDzhoSE0KhRoww/b1ZlreXpnTuZePgwTxYtynvly+Oe+N9na9as4eGHH76sXeX555+nS5cuTkUW0u73wZw5c+jQoUPScvuuXTnbvz+LTp3ivXLl6Fu8+E2fI6vSn0XO0xg4L63GYOvWrVSuXPnmA2VDYWFh5MlCN+qn9mvBGPO3tTY4tf1d8kq5MaaZMWa7MWaXMWZwKtvHG2PWJb52GGPOJNsWl2zbvAwNLo5JXpAPKFGCiYkF+cV2ldtvvz1FQX7xYUCZ/U5t+Vf79u156623kpbnzJhB0Ndf06pQIZ7etYsx+/c7mE5EROTqXO5GT2OMOzARuBs4CKwyxsyz1m65uI+1tl+y/Z8GaiZ7iyhrbVAGxRUXcGlBPrpMGYwxnD17ll69evHll18m7evj48NHH32kucezqEGDBrF3796k+Xbfev11PgwIwKtuXQbu2UNUfDwvBwQk3YAkIiLiKlzxSvltwC5r7R5rbQwwB3jwKvt3AD7PkGTicqy1PLNrFxMPH+b54sWTCvI1a9ZQu3btFAV5jRo1+Pvvv1WQZ2HGGCZOnEizZs2S1j31+ON0PXSIbkWK8GpoKMMSp+YSERFxJa5YlBcDDiRbPpi47jLGmACgNLA02WovY8xqY8wfxphW6ZZSHHexIH//0CGeL16cMWXLAly1XUWzq2R9OXLk4IsvviAoKAiACxcu8HCbNjx7/jyPFSnCsH37GKVWFhERcTEud6OnMaYt0Mxa2zNxuQtQx1rbN5V9BwHFrbVPJ1tXzFp7yBhThoRivYm1dncqx/YGegMULly49pw5c9LnA11FeHg4Pj4+GX7erMAC7wHfAO2AJ4HIiAjGjh2b9JhfSJjL+vnnn+euu+5K9X00Bs5LrzE4ceIEffr04fjx4wAUKlSI9z/4gCl+fiwF+gJt0vysmZN+HzhPY+C8tBqDfPnyUa5cuTRIlP3ExcXhnjhjWlawa9cuzp49m2Jd48aNr3ijJ9Zal3oBtwM/Jlt+EXjxCvuuBe64yntNB9r+1zlr165tnbBs2TJHzpvZxcfH2/47d1qWLbP9d+608fHxds2aNbZs2bKWhHrdAjYoKMju2LHjqu+lMXBeeo7Bhg0bbN68eZN+TVSvXt2eOHPGPrRxo2XZMvvhoUPpdu7MRL8PnKcxcF5ajcGWLVuu+5jY2DN2w4ZWNjb2TJpkyKzOnTvndIQ0ldqvBWC1vUI96ortK6uA8saY0saYnEB74LJZVIwxlQBfYGWydb7GGM/E7wsB9YAtlx4rmduIffsYd/AgTxcrxpgyZZg0adJlDwN64oknWLlyJeUdepqauIZq1arx9ddfkyNHwj3tGzZsoOPDD/NZ+fI0L1CAJ3bs4DM9YEhEHHbixDxOnvyWEyfmOx3lujVq1Ij0nFZ63rx5jBw58qr7hIaGMnv27HTLkFFcrii31l4g4X+WfwS2Al9YazcbY143xjyQbNf2wJzEf3VcVBlYbYxZDywDRtpks7ZI5vf+wYO8HBrKo4ULM8zPj/bt26d4GFCePHmYO3cukyZN0sOABICmTZvy0UcfJS0vXryY5/r25asqVbgrf366bdvGl8eOOZhQRLK7I0empfgq/3rggQcYPPiy2bFTUFGejqy1C621Fay1Za21IxLXvWKtnZdsn9estYMvOe53a201a22NxK8fZ3R2ST+zjh7l6V27eLBgQZ6KjOTW4OAUs6sEBQUlPSRIJLlu3brxyiuvJC1PnTqV8aNG8V21atyRLx8dt27l+xMnHEwoItnJunVNCQkxSa+zZ38H4OzZFSnWr1vX9LrfOyIighYtWlCjRg2qVq3K3LlzAXj99de59dZbqVq1Kr17977Y5kujRo3o168fwcHBVK5cmVWrVtG6dWvKly+f9ITN0NBQKlWqRKdOnahcuTJt27YlMjLysnMvXryY22+/nVq1atGuXTvCw8Mv26dRo0Y8++yzBAUFUbVqVf766y8ATp06RYcOHahevTp169Zlw4YNAEyfPp2+fRNuK+zWrRvPPPMMd9xxB2XKlOGrr74CYPDgwfz2228EBQUxfvz46/6ZuQqXLMpFLjX/xAm6bt1K4/z5abx8OQ3uuCNFu0qfPn1YuXKlbq6RK3rttddSPL116NChfDd3LguqVSPIx4eHt2xh5SU35IiIpIeAgKG4uXknLSfMAP3vVwA3N28CAl667Nj/8sMPP1C0aFHWr1/Ppk2bkqaI7du3L6tWrWLTpk1ERUXx/fffJx2TM2dOVq9ezRNPPMGDDz7IxIkT2bRpE9OnT+fkyZMAbN++nT59+rB161by5s3LBx98kOK8J06cYPjw4fz888+sWbOG4OBgxo0bl2rGyMhI1q1bxwcffMBjjz0GwKuvvkr16tXZsGEDb775Jo8++miqx/7zzz8sX76c77//PukK+siRI2nQoAHr1q2jX79+qR6XGagoF5cXcvo07TZvpoa3N/nGjOG5Pn0ua1eZOHGi2lXkqowxTJ06lcaNGyet6969O2tXrGB+5ZKM4GUe3rCSrRERDqYUkezA17cx1ap9n6IwT87NzZtq1Rbg69vout+7WrVq/PTTTwwaNIjffvuNfPnyAbBs2TLq1KlDtWrVWLp0KZs3b0465oEHHkg6tkqVKtxyyy14enpSpkwZDhxImKW6RIkS1KtXD4DOnTuzfPnyFOf9448/2LJlC/Xq1SMoKIgZM2awb9++VDN26NABgDvvvJNz585x5swZli9fnvQckbvuuouTJ09y7ty5y45t1aoVbm5uBAYGcvTo0ev++bgyFeXi0v4OC+OBTZsoZgynevXi22Q9YzVr1uTvv/9Wu4pcs5w5c/K///2PwMBAAGJiYmjVqhWHt02jZvyv1GUl927YwMHoaIeTikhW5+vbmMDAubi5pbyg5ObmRWDg3BsqyAEqVKjAmjVrqFatGi+99BKvv/460dHR9OnTh6+++oqNGzfSq1cvopP9Oefp6Zl4brek7y8uX7hwAeCyJyFfumyt5e6772bdunWsW7eOLVu28PHHqXcR/9d7XU3yfClvK8z8VJSLy9oTFUXzDRvwiIriQPv2hCb2l0FCu8rvv/+u2VXkuuXPn5+FCxdSpEgRAM6cOcOff74KwLO5Qjhz4QLNNmzgdGyskzFFJBu4cOEMkANww80tFwllWY7E9Tfm8OHDeHt707lzZ1544QXWrFmTVIAXKlSI8PDwpF7s67F//35WrkyY8G727NnUr18/xfa6deuyYsUKdu3aBST0tu/YsSPV97rY5758+XLy5ctHvnz5aNCgAV988QUAISEhFCpUiLx5815Ttjx58hAWFnbdn8nV5HA6gEhqTsbGcu+6dZwJCyOmd2/45x8g4Tfe1KlTdXVcbsrp0z34/PN/p0KMiTkPQFzEH8yzd0IkrF8BefPfRa2gJU7FFJEs7siRj4mPj8THpwZlyoxiz55BhIev58iRaRQp0vmG3nPjxo288MILuLm54eHhwaRJk8ifPz+9evWiatWqFClShFtvvfW637dixYpMnDiRxx57jMDAQJ588skU2/38/Jg+fTodOnTg/PmEP1OHDx+e6pO0vby8qFmzJrGxsUybljDjzGuvvcajjz5K9erV8fb2ZsaMGdecrXr16ri7u1OjRg26deuWafvKXe6Jnk4IDg626TnH5pWEhITQqFGjDD+vq4uKi+P2339nQ3Q0tn9/2LQJSGhX+eKLL9L0Zk6NgfOcGIPTp5excWNL4uMvnz3gomg8+T7vB0wI6kYOt6z9n4r6feA8jYHz0moMtm7dSuXKla9p340bW5E//50UL/4cxrhhbRwHDkzg7NnfqFbt25vOklZCQ0Np2bIlmxL/Pr4ZjRo1YuzYsQQHX/5Qy7CwMPLkyXPT53AVqf1aMMZc8YmeWftvGsl04uLjqff996yPjcWOGJFUkF9sV9HsKpIWruUmq3+KfsbEc2Xou3NnlutbFBHXUK3at5Qo0R9jEsoxY9wpWfJ5lyrIJeOoKBfnzZoFpUpx1hgqP/cca/Plg8mT4ZdfNLuKpJsr3WR1/jzs3NmB7hXaMbhkST785x/GHTzoUEoREeeVKlUqTa6SQ8L/SKR2lVzUUy5OmzULevfm78hI7n3oIU62bg3/+x98+WW6tKuIJJf8JqvYWIObWxxxcTB37sf4+jZnxEMPsTsqihd276aslxet/PwcTiwiIlmVrpSLs4YO5YfISOrWq8fJvn3ht99g4kT6+PioXUXSXfKbrKpV+5ajR73x8oJmzaBTp0789eefzKhUidvy5KHT1q38nQXu7heR9KeWN7mRXwMqysVRR/bto0P58lx46SXYtg2fESOYGx/PxIgItatIunN3z0fZsmOoXXs1RYq05O679/DVVwWJjITo6Gjuv/9+DoeG8l21avh5eHD/xo0c0BzmInIVXl5enDx5UoV5Nmat5eTJk9ddx6h9RRxjraXDLbdwZsQIOHOGwi+9xPLz5ykHULKk0/EkG7j0Zio/v8L07/8HdevWBU5y4sQJ7rvvPlauXMmC6tW5Y80aWm7cyPKaNcmTQ398isjlihcvzsGDBzl+/LjTUTKd6OjoLHNBzsvLi+LFi1/XMfpbRRwzbuJEQl5+GXx84OmnmXX6dEJB7u0NI0Y4HU+yqXLlyjFv3jzuuusuzp8/z86dO3nwwQf5+eef+aJKFVps2ED7LVv4rmrVLD9VoohcPw8PD0qXLu10jEwpJCSEmjVrOh3DMfobRRyxceNGBh4/DpUrw4gRDDh+nCbGQEAATJkCnTo5HVGysTvuuIOZM2cmPfp5xYoVdOvWjbvz5+f98uVZeOoUz+/e7XBKERHJSlSUS4aLjo7mnhkziG/cGKZOJSgiguHHjkF8PISGqiAXl9C2bVvGjBmTtDx37lyGDBnCE8WK0a94cd49dIiPE580KyIicrPUviIZ7uF33+VIy5bw0094fv01s/7+G09PT6djiVymf//+7N27l4kTJwIwatQoSpcuzehevdgUEcGTO3ZQ2dubO/LlczipiIhkdrpSLhlq0k8/Mb9aNdiyBcaMYdzbbxMYGOh0LJFUGWN45513uP/++5PW9enTh8U//MCcwEBKenrSetMmDmpGFhERuUkqyiXDbDtyhL5nz0JYGLz0Ei3uuYcnn3zS6VgiV+Xu7s7nn3+e9AS6+Ph4Hn74YUI3buS7atWIiI+n9ebNRMfFOZxUREQyMxXlkiFi4+Np8NNPxOfJAy+9hL+HBx9//HHSjXQirix37tzMnz+fgIAAACIiImjRogV5Tp5kZuXKrAoLo/eOHZqXWEREbpiKcskQzb/7jhMlSsDbb8POnUybNo3ChQs7HUvkmhUpUoRFixaRP39+AI4cOULz5s1pmCMHw0qV4rOjR5lw8KCzIUVEJNNSUS7p7u116/jZ1xf+9z/46SeeeuopWrRo4XQsketWuXJlvvnmGzw8PADYvHkzbdq0YeAtt9C6UCEG7N7NT6dOOZxSREQyIxXlkq5Wnz7NC8eOwfr18MEHVK5cOcU0cyKZTaNGjfjkk0+SlpcuXcoTjz/O9EqVCMydm0e2bGFvVJSDCUVEJDNSUS7p5lRsLHevXIk9exaGDcPDzY3Zs2eTK1cup6OJ3JROnToxfPjwpOUZM2YwbsQIvq1alXhraasbP0VE5DqpKJd0EWctzVas4IyHB7z6Kpw+zVtvvUVQUJDT0UTSxJAhQ+jRo0fS8muvvcbyL77gs8qVWRMeTt+dOx1MJyIimY2KckkXA7duZRXAu+/C1q00adKEfv36OR1LJM0YY5g0aRL33HNP0rqePXvivX49Q0qW5OMjR/TETxERuWYqyiXNfXf8OOOOHYPvv4fvv6dAgQLMmDEDNzf9cpOsxcPDgy+//JLq1asDcOHCBVq3bs3DERE09fXlqR07WBMW5nBKERHJDFQlSZraGxVFx40bYfv2hKvkwJQpUyhWrJjDyUTSR968eVmwYEHSr/Fz585xf4sWjPP1xT9nTtps3syp2FiHU4qIiKtTUS5pJjoujgfWriUyMhJeew1iY+nRowdt2rRxOppIuipevDgLFiwgT548ABw4cICurVoxo3RpDp0/T5etW4nXg4VEROQqXLIoN8Y0M8ZsN8bsMsYMTmV7N2PMcWPMusRXz2Tbuhpjdia+umZs8uztuZ072RQTA2+9BUeOUK5cOSZMmOB0LJEMUaNGDb788kvc3d0BWLt2LWN79GBcmTIsPHWKEfv2OZxQRERcmcsV5cYYd2AicB8QCHQwxgSmsutca21Q4mtq4rEFgFeBOsBtwKvGGN8Mip6tzT56lA+PHIE5c+D333F3d2fWrFn4+Pg4HU0kw9x77718+OGHScsLFy5k05tv0tnfn1dDQ/lZDxYSEZErcLminIRiepe1do+1NgaYAzx4jcfeC/xkrT1lrT0N/AQ0S6eckmhrRAQ9t26FjRth6lQAhg0bxm233eZwMpGM16NHD4YOHZq0/OHkyVRcuJDK3t502rqVf86fdzCdiIi4KlcsyosBB5ItH0xcd6k2xpgNxpivjDElrvNYSSMRcXG03riR8+fOweuvQ1wc9evXZ/Dgy7qORLKNN954g44dOyYtv/zCC3Tfv5+wuDg6bd1KnPrLRUTkEsa62F8Oxpi2QDNrbc/E5S5AHWtt32T7FATCrbXnjTGPA49Ya+8yxgwAvKy1wxP3exmIstaOTeU8vYHeAIULF649Z86cdP9slwoPD8/U7R0WeAv4yVoYMADWrCF37txMnTqVIkWKOB3vmmT2McgKsuoYxMTEMHDgQNavXw8kTJ/48McfM6tECboC3RxNl1JWHYPMRGPgPI2B87LDGDRu3Phva21wattyZHSYa3AIKJFsuXjiuiTW2pPJFqcCo5Md2+iSY0NSO4m1dgowBSA4ONg2atQotd3SVUhICE6cN618/M8//LR9O0yfDmvWAAnTH7Zv397ZYNchs49BVpCVx+DWW2/ljjvuYNu2bcTGxrLoued4cOFCPo2K4tEaNbjL1zVuecnKY5BZaAycpzFwXnYfA1dsX1kFlDfGlDbG5ATaA/OS72CMuSXZ4gPA1sTvfwTuMcb4Jt7geU/iOkljWyIi6LtjBznWr4eZMwHo2LFjiv+yF8nufH19WbhwIf7+/gCcOnWK9d26UT5nTjpu2cIR9ZeLiEgilyvKrbUXgL4kFNNbgS+stZuNMa8bYx5I3O0ZY8xmY8x64BkS/yfYWnsKeIOEwn4V8HriOklDUXFxPLJ5M3FhYVwYNgzi4wkICGDixIlORxNxOaVLl+b777/H29sbgNBt2/AcNYpzcXF0Vn+5iIgkcrmiHMBau9BaW8FaW9ZaOyJx3SvW2nmJ379ora1ira1hrW1srd2W7Nhp1tpyia9PnPoMWVn/3bvZFBlJ7Ouvw+nTuLm58dlnn5E/f36no4m4pFtvvZXPP/8cYwwAG+fNo9LixSw5c4Y3NX+5iIjgokW5uK4vjx1j8uHDuM2dC6tWAfDiiy/SoEEDh5OJuLYHHniAd955J2l57YgRVDpwgNdCQ/nlzBnngomIiEtQUS7XbG9UFL22b8drzx7iP/oISLgC+OqrrzqcTCRzePrpp+nXr1/S8rbevSl4/jydt27lZGysg8lERMRpKsrlmsTGx9NhyxaioqKIHjoU4uLInTs3s2bNwsPDw+l4IpnG2LFjad26dcJCdDTHn3mGI9HR9Ny+HVebolZERDKOinK5Ji/v3cufYWHEvPkmHDkCwLvvvkv58uUdTiaSubi5uTFz5kzq1q2bsGLHDpg6lW9PnODDw4edDSciIo5RUS7/afGpU4w6cACvxYvhl18AaN26Nd27d3c4mUjmlCtXLr777jvKlCkDwIXPPyfn+vU8t2sXmyMiHE4nIiJOUFEuV3U8Joau27bhc+IE0W+/DUDRokWZMmVK0kwSInL9/P39WbRoEQUKFABriRk2jAtnz9Ju40ai4uKcjiciIhlMRblckbWWHtu3cyI6mvBBgyAmBoBPP/2UggULOpxOJPOrUKEC3333HZ6ennD6NHHDh7M1Opr+O3Y4HU1ERDKYinK5osmHDzP/5Enshx/Cnj0APP/88zRp0sThZCJZR/369ZkxY0bCwqpV8OWXTD56lG+PH3c2mIiIZCgV5ZKqLRER9N+1C5+tW4n78ksAatSowYgRIxxOJpL1PPLII4wcOTJh4aOPYMcOOq5dy6Hz550NJiIiGUZFuVzmfHw8HbdswURHEz50KFiLl5cXs2fPTvhvdhFJcwMHDuTxxx+H2FgYPpyouDjuXrKEeE2TKCKSLagol8sM2bOH9RERRL32Gpw+DSTMrRwYGOhsMJEszBjD+++/T/PmzeHAAfjgA7Z6e/Pkzz87HU1ERDKAinJJYfGpU4w7eBDvxYvhjz8AaN68OX369HE4mUjWlyNHDubOnUvNmjXh++9hxQqmWMtXa9Y4HU1ERNKZinJJcnH6wzwnTxKZOP2hv78/06ZN0/SHIhnEx8eH77//nhIlSsCYMRAeTodNm9i1f7/T0UREJB2pKBcgYfrDXonTH4YNHJg0/eEnn3xC4cKFHU4nkr0ULVqUhQsXktdaGD2aCyVLcvvHH3Pu3Dmno4mISDpRUS4AfHLkCN+dPImZNi1p+sOnnnoqob9VRDJc1apV+d///keOv/+Gb7/lROPGNB4wgNjYWKejiYhIOlBRLuyJiuLZnTvJs2sXsbNnA1C5cmXGjBnjcDKR7K1JkyZMnToVJk+G/ftZc/fdPFayJNYYKFUKZs1yOqKIiKQRFeXZXJy1PLp1K7HnzxOWOP2hh4cHs2fPJleuXE7HE8n2unbtyquDB8OIEeDry8zHH2cEwL590Lu3CnMRkSxCRXk2N2b/flacO8f5MWPg2DEA3nzzTYKCgpwNJiJJXn31VboeOgTTpkGjRrx8zz0sBoiMhKFDnY4nIiJpQEV5NrY2LIxXQkPJ9eef8NNPQMJ/l/fv39/hZCKSnDGGKRER3DV3LqxfD08/TTd/f84CaFYWEZEsQUV5NhUdF0fnrVtxDw8n6s03AShQoAAzZszAzU2/LERcTc6AAD6Pj8d31Chwc+OfQYN43hgoWdLpaCIikgZUfWVTQ/buZUtkJNHDhkHiNGsfffQRxYoVcziZiKRqxAj8vb2Z/M8/MHEi1KrFxw89xOL27Z1OJiIiaUBFeTa05PRpxh88iMeCBbBqFQA9evSgdevWDicTkSvq1AmmTOHhgADaLFwIK1dC7950XbaMs2fPOp1ORERukorybObshQt037YNr+PHiX3vPQDKlSvHhAkTnA0mIv+tUycIDeWDo0cpMG0aREdzpHt3+r/wgtPJRETkJqkoz2b67drFwehool95Bc6fJ0eOHMyePRsfHx+no4nINfL392fSiBEwbhxUqsS08+f58ccfnY4lIiI3QUV5NjL/xAk+OXIEZs+GbdsAeO2117j11lsdTiYi1+vhhx+mrb9/wsxJXbrw6IgRamMREcnEVJRnEydiYui5bRs59+/HTp8OQP369Rk8eLCzwUTkhk2cOJECM2fC6dMce+wxnhs0yOlIIiJyg1SUZxNP7dzJifPniRk2DC5cIG/evMycORN3d3eno4nIDfL392fSmDEwahSUKsV0d3e1sYiIZFIqyrOBOUeP8sXx48RPmwZ79gAwefJkAgICHE4mIjfr4Ycfpm3p0vDNN9CmDV3GjFEbi4hIJuSSRbkxppkxZrsxZpcx5rL+CmNMf2PMFmPMBmPMEmNMQLJtccaYdYmveRmb3PX8c/48T2zfjvv27TBnDgCdOnWiQ4cODicTkbQyceJECn79NfzzD8cfe4xn1ZYmIpLpuFxRboxxByYC9wGBQAdjTOAlu60Fgq211YGvgNHJtkVZa4MSXw9kSGgXZa2l5/bthJ0/T9zw4RAfT0BAABMnTnQ6moikIX9/fz4YNy6hjaVIEWZ4eqqNRUQkk3G5ohy4Ddhlrd1jrY0B5gAPJt/BWrvMWhuZuPgHUDyDM2YK044cYeGpU8RPngwHD+Lm5sbMmTPJly+f09FEJI09/PDDtK1YEb7+Glq1ovO4cWpjERHJRFyxKC8GHEi2fDBx3ZX0ABYlW/Yyxqw2xvxhjGmVDvkyhf3R0Ty7fTtm3bqEXlNgyJAh1K9f39lgIpJuJk6cSMFvv4X9+znRvTvPqI1FRCTTMNZapzOkYIxpCzSz1vZMXO4C1LHW9k1l385AX6ChtfZ84rpi1tpDxpgywFKgibV2dyrH9gZ6AxQuXLj2nMR+64wUHh6eLg/tscDzcXGsi4nBdu8OR49SqVIl3nvvPXLkyJHm58vM0msM5NppDNLWsmXLeP3LL+G99+CHHxjl58dtt9121WM0Bs7TGDhPY+C87DAGjRs3/ttaG5zaNles0A4BJZItF09cl4IxpikwlGQFOYC19lDi1z3GmBCgJnBZUW6tnQJMAQgODraNGjVKu09wjUJCQkiP8354+DBrd+yASZPg6FFy587N/PnzKVeuXJqfK7NLrzGQa6cxSFuNGjViy5YtfDVnDnTqxOixY9n9+ONXbVvTGDhPY+A8jYHzsvsYuGL7yiqgvDGmtDEmJ9AeSDGLijGmJvAh8IC19liy9b7GGM/E7wsB9YAtGZbcBYRGRdFv+3ZYvRrmzwfg3XffVUEuko1MnDiRgt9/D3v3crJrV54eMsTpSCIi8h9crii31l4goSXlR2Ar8IW1drMx5nVjzMXZVMYAPsCXl0x9WBlYbYxZDywDRlprs01RHm8tXTZuJDoqCsaOBaBNmzZ0797d4WQikpH8/f354J13YORIKFCAz/Ll02wsIiIuzhXbV7DWLgQWXrLulWTfN73Ccb8D1dI3neuadOgQyyMj4YMP4OhRihUrxpQpUzDGOB1NRDLYww8/zJdffslXM2dC1650GjeO3XXravYlEREX5XJXyuXG7ImKov/27bBqFSxYAMCMGTMoUKCAw8lExCkTJ06k4A8/wO7dnOzalb5qYxERcVkqyrOAeGt5ZM0aYqKjYcwYAAYMGECTJk0cTiYiTkrRxpI/PzPVxiIi4rJUlGcB7+zbx+rY2IS2lePHCQoKYvjw4U7HEhEX8PDDD9M2KAhmzYJ77qHTxIl6qJCIiAtSUZ7J7Y6KYuCuXfDnn7BwIV5eXsyePRtPT0+no4mIi5g4cSIFFy1KaGN59FH6vvii05FEROQSKsozsXhrab1yJReio+HttwF4++23qVy5ssPJRMSVJLWxjBoFvr5qYxERcUEqyjOxsTt2sMHNLeEhQceP06JFC5588kmnY4mIC3r44YdpW6NGQhvLvffS6f331cYiIuJCVJRnUnsiIxmyb1/CbCsLF+Lv78+0adM0/aGIXFGKNpauXdXGIiLiQlSUZ0LWWlqGhBAXG5v0kKBPPvkEf39/h5OJiCvz9/dn0rvvpmhj+eGHH5yOJSIiqCjPlIavW8dWb2/48EM4doy+ffvSvHlzp2OJSCbQrl27FG0snd9/n/DwcKdjiYhkeyrKM5nd4eG8dvQorFkD8+cTGBjI6NGjnY4lIpnIBx98kNDGsncvJ7t25d2PP3Y6kohItqeiPBOx1nL3Dz8Qby2MHk3OnDmZNWsWuXLlcjqaiGQifn5+/7axFCjAT2XLqo1FRMRhKsozkcG//sreQoVgyhQ4epS33nqLoKAgp2OJSCbUrl072lWvDl98AS1b0uWddzQbi4iIg1SUZxKbjx1jTEQErFsH331H06ZNee6555yOJSKZ2Pvvv0/B77+H/fs50bUrTw8a5HQkEZFsS0V5JhAfH889P/yAdXODMWMo4OvL9OnTcXPT8InIjfP392fShAkwejT4+/OZp6faWEREHKKqLhN46ttvOVyyJEydCocPM3XqVIoVK+Z0LBHJAtq1a0fDQoXg66/hoYd49O231cYiIuIAFeUu7q9du5js4QGbNsE339CjRw8eeughp2OJSBby7LPPUvC77+DQIY5368azAwc6HUlEJNtRUe7CLly4QPMffgBPTxg9mnJlyjBhwgSnY4lIFuPr68uk8eNhzBgoVowZOXKojUVEJIOpKHdhnT/6iJNVq8L06eT45x9mz56Nj4+P07FEJAtq164d7SpUgO++gzZteHT0aLWxiIhkIBXlLuqHlSuZW6QIbN8Oc+fy2muvceuttzodS0SysIkTJ1Lgq6/g2DGOd+3Ksy+84HQkEZFsQ0W5CwoLC6PdL79AnjwwejT177iDwYMHOx1LRLI4Pz8/Jo8bB2+/DQEBzLBWbSwiIhlERbkLemjMGMLr1oWZM8l74gQzZ87E3d3d6Vgikg20a9eOdqVLw6JF0L49XYcPVxuLiEgGUFHuYmZ8/TVLqlWDPXtg1iwmTZpEQECA07FEJBt5//33KTB3Lpw+zbHu3XlObSwiIulORbkLOXjwIL3XroUCBWDUKDo98ggdO3Z0OpaIZDP+/v5MHjMGJkyAsmWZHh2tNhYRkXSmotxFxMfH0/K114hp2hTmziXg/HkmTpzodCwRyabatWtHu6JFYckS6NKFbsOGqY1FRCQdqSh3EW+OH8/6pk1h/37MZ58xc+ZM8uXL53QsEcnGJk6cSIHZsyEigqOPPkq/AQOcjiQikmWpKHcBa9eu5dUDB8DfH8aMYeiAAdSvX9/pWCKSzfn5+TF55Eh4912oXJlPzp1TG4uISDpRUe6wyMhIWr38MvGtWsG333KbtzevvPKK07FERICENpa2hQrB8uXw2GN0e/lltbGIiKQDFeUOe27QIPY//DD88w/es2czc+ZMPDw8nI4lIpLkg4kTKfDZZxATw9HOnenXv7/TkUREshyXLMqNMc2MMduNMbuMMZc9NccY42mMmZu4/U9jTKlk215MXL/dGHNvhga/TitXruSj8+ehZEkYO5b3Ro+mfPnyTscSEUnBz8+PD0eMgA8+gBo1+OTUKRYtWuR0LBGRLOU/i3JjzC/GmFoZESbxfO7AROA+IBDoYIwJvGS3HsBpa205YDwwKvHYQKA9UAVoBnyQ+H6uZdYsjpQowejp0+GRR2DBAlqXKkX37t2dTiYikqq2bdvSNk8eWLUKevWi++DBnDlzxulYIiLXZ9YsdhUrRpwxUKoUzJrldKIk13Kl/ADwpzFmhjGmWHoHAm4Ddllr91hrY4A5wIOX7PMgMCPx+6+AJsYYk7h+jrX2vLV2L7Ar8f1cx6xZ2F696HrkCGcGDoRTpygyeTJT7r6bhI8gIuKaPpg4kQLTp4MxHO3cmf7PP+90JBGRazdrFsd69eL2qCga58jBgX37oHdvlynM/7Mot9Z2BuoBZYEdxpjXjTG50zFTMRL+IXDRwcR1qe5jrb0AnAUKXuOxzho6lIlRUSzu2BHKloUJE/gsPJyCI0c6nUxE5Kr8/Pz4cNgwmDIFbr2VTw4fVhuLiGQadsgQHouO5sRrr/Hb2LHcA8RHRsLQoU5HAyDHtexkrf0LqG+MeQQYCfQ0xrwMTLPW2vQMmF6MMb2B3gCFCxcmJCQkQ87bcP9+SgL5ly3jTFwcz//+O00Bu38/v2RQBvlXeHh4ho29pE5j4LzrGYNChQrR8MwZflm/Hvr0ofMzzzBj/Hh8fHzSN2QWp98HztMYOC+9x2Dr/v0saNUKqleHkSOZQMLVaZepway11/UCPIHBwBlgHdD0et/jP97/duDHZMsvAi9ess+PwO2J3+cATgDm0n2T73e1V+3atW2GCQiwFuwhsH3BRoO1kLBeMtyyZcucjpDtaQycd71jcOzYMetbvbrlhx8sr79uu3Xvnj7BshH9PnCexsB56TkGW7ZssZ5FilgWLrSMGmWfuVh/ZXANBqy2V6hHr3n2FWNMTmPMbUBPoBJwDqgO/GiMmW+MKX0z/zhIZhVQ3hhT2hiTk4QbN+ddss88oGvi922BpYkfdB7QPnF2ltJAeeCvNMqVNkaMAG9vigLvkfAvHLy9E9aLiGQCfn5+THn5ZZg+HRo0YPrevSxcuNDpWCIiqYqJiaFjp06cf/55sJYK48YlzBACLlWDXcvsK5ONMatJKML/AEYARUm40fIB4A7AAhuMMffdbCCb0CPel4Sr3FuBL6y1mxN72R9I3O1joKAxZhfQn4Qr91hrNwNfAFuAH4CnrLVxN5spTXXqlNCPGRCANQYCAhKWO3VyOpmIyDVr27Ytba2FbdvgmWfo0b+/ZmMREZf08ssvs65IEQgOJsfHH/OVmxteLliDXUtPeR3gTxKmKfwD2JZ4VTq5B4wxo4F3Sbg6fVOstQuBhZeseyXZ99FAuyscO4KEfzi4rk6doFMnfgkJoVGjRk6nERG5IR+89x5LmjXj9MiRHGnXjv79+zNt2jSnY4mIJFm6dCmjp02DTz6B9esZfeedVPv6a6djpepaZl+paa19wlr7ibV2ayoF+UVfAWXSNp6IiLgqPz8/PhoyBGbOhCZN+GTbNrWxiIjLOHXqFF0efRSeew48PKj/xx88+/TTTse6orR8oud6Eh74IyIi2USbNm1oExMDu3dDv370ePZZtbGIiOOstTz++OMcrlAB6tXDe+5cvhg3Djc3l3yYPZCGRblNeGDP4rR6PxERyRwmvfsu+T/6CAoU4EirVvTv39/pSCKSzc2YMYOvfvoJnnkGtm7lswcf5JZbbnE61lW57j8XREQkU/Dz82PqCy/A3LnQogWfbNigNhYRcczu3bt5+umn4emnIXdu2u7ZQ+sHL304vOtRUS4iIjetTZs2tA4Ph/37YcAAevTtqzYWEclwsbGxdOrUifAaNaBJEwouWsT0115zOtY1UVEuIiJpYvI77yS0sfj7c6RlS7WxiEiGGz58OH9u3pxwc+fu3Szo2JHcuXM7HeuaqCgXEZE04efnx9TnnoNvvoHWrflk9Wq1sYhIhlmxYgXDhw+HPn3A15e+ERHUCQ52OtY1U1EuIiJppk2bNrQ5fRoOH4YXXqBnnz5qYxGRdHfu3Dk6d+5MfO3acN99lFyxgglPPeV0rOuiolxERNLUpPHjyTd1KpQowT/NmqmNRUTSXd++fQk9dgyefx63Awf4uWtX3N3dnY51XVSUi4hImvLz8+Pjp56C77+Hdu345I8/1MYiIulmzpw5fPbZZ9C7N/j5MSx3bsoHBDgd67qpKBcRkTTXpk0bHjp2DE6ehIED6fnkk2pjEZE0t3//fp544gmoUQMefJDALVt4qXVrp2PdEBXlIiKSLqaMG0feqVOhVCn+adpUbSwikqbi4uLo0qULZ8+fhxdeIMfRoyzp3NnpWDdMRbmIiKSLQoUKMe3xx+GHH6BjRz757Te1sYhImhk9ejS//vorPPYYFCvGO8WKUSR/fqdj3TAV5SIikm7atGnDQ4cPw5kzMGgQPZ94Qm0sInLTVq9ezSuvvAJVqkCbNgQfOkSf+vWdjnVTVJSLiEi6mjJuHPmmTYNy5fincWP69evndCQRycQiIiLo2LEjF9zcYOBAcp49y09t2zod66apKBcRkXRVqFAhpvXoAUuWQOfOTP/lF7WxiMgN69evHzt37oTu3aFkSaZWrEh+T0+nY900FeUiIpLuWrduTasDByA8PKGN5fHH1cYiItft22+/5aOPPoLKlaFdOxqFhdElMNDpWGlCRbmIiGSIj8aOJe/06VCxIv80aKA2FhG5LocPH6Znz57g4QEDB+IdGcm3zZo5HSvNqCgXEZEMUahQIT7p2hV++QW6dWP6smUsWLDA6VgikgnEx8fTrVs3Tp48CV27QqlSzKhenXweHk5HSzMqykVEJMO0bt2aVvv2QVQUDBpEL83GIiLX4N133+Wnn36CihWhfXvus5a2mfCpnVejolxERDLUR6NHJ7SxVK7MP/Xrq41FRK5qw4YNDBo0KKFtZdAgfGJjmZ3Jpz9MjYpyERHJUIUKFWJaly7w66/QvTvTly5VG4uIpCoqKoqOHTsSExMDnTtD6dLMrFmT/FmobeUiFeUiIpLh2rRuTavQ0IQ2lsGD1cYiIqkaPHgwmzdvhvLloVMnHvTy4sHChZ2OlS5UlIuIiCNStLHUq6c2FhFJYdGiRbz77rsJbSuDB5PPWj6pXdvpWOlGRbmIiDgiRRvLY4+pjUVEkhw7dozu3bsnLHTpAmXKMDMoCN8s2LZykYpyERFxTIo2lkGD6KWHColke9ZaevbsydGjR6FCBejYkfb589OyUCGno6UrFeUiIuKopDaWwEDNxiIifPjhh8yfPz+pbaWgmxuTqlRxOla6U1EuIiKOuqyNRQ8VEsm2tm3bRv/+/RMWunaF0qX5rHr1LDnbyqVcqig3xhQwxvxkjNmZ+NU3lX2CjDErjTGbjTEbjDGPJNs23Riz1xizLvEVlKEfQEREbshlbSyajUUk24mNjaVjx45ERUVBpUrQoQNd/fy4r2BBp6NlCJcqyoHBwBJrbXlgSeLypSKBR621VYBmwARjTP5k21+w1gYlvtald2AREUkbH40eTd5p0xJmY2nYUG0sItnMtGnTWLt2LXh4YAYPpnCOHLxTsaLTsTKMqxXlDwIzEr+fAbS6dAdr7Q5r7c7E7w8DxwC/jAooIiLpo1ChQnzStSssWwZduzL911/VxiKSTSxbtoy5c+cmLHTvjg0IYEbVquTLkcPZYBnI1YrywtbafxK/PwJcdXZ4Y8xtQE5gd7LVIxLbWsYbYzzTKaeIiKSD1q1b89D+/RAWlvBQoSef5PTp007HEpF0dOrUKbp06YK1FipXhkceoWeRItxboIDT0TKUsdZm7AmN+RkoksqmocAMa23+ZPuettZe1leeuO0WIAToaq39I9m6IyQU6lOA3dba169wfG+gN0DhwoVrz5kz50Y/0g0LDw/Hx8cnw88r/9IYOE9j4DxXG4OzZ8/S6YMPiHjxRZgxg2ZHjjBo0CCnY6UrVxuD7Ehj4AxrLcOGDeOXX34BT0/cPv6YAkWKMN3dndxOh0sHjRs3/ttaG5zatgz/PwFrbdMrbTPGHDXG3GKt/SexwD52hf3yAguAoRcL8sT3vniV/bwx5hNgwFVyTCGhcCc4ONg2atTouj/LzQoJCcGJ88q/NAbO0xg4zxXH4NO4ONosXgydO/NDnz70jYigRYsWTsdKN644BtmNxsAZ06dPTyjIAXr2JL5YMebWqMFdvqlek83SXK19ZR7QNfH7rsB3l+5gjMkJfAN8aq396pJttyR+NST0o29Kz7AiIpI+WrduTetDh+DUqYQ2lj591MYiksXs3r2bp59+OmEhKAjatuXpYsWyZUEOrleUjwTuNsbsBJomLmOMCTbGTE3c52HgTqBbKlMfzjLGbAQ2AoWA4RmaXkRE0syHb79NvqlToXRp/rnnHs3GIpKFxMbG0qlTJ8LDw8HbmxxDh1I0Pp6RZco4Hc0xLnVLq7X2JNAklfWrgZ6J388EZl7h+LvSNaCIiGSYQoUKMe3xx2mzYAG0b8+MZ56h3YIFWbqNRSS7GD58OH/++ScA5qmniC9UiCGAt7u7s8Ec5GpXykVERJK0bt2aNkePwrFj8OKL9Hz6abWxiGRyv//+O8OHJzYz1K2Lbd6cgSVLUsXZWI5TUS4iIi5t8rhx5P/wQyhalCMPPKA2FpFM7Ny5c3Tq1In4+HjIm5ecQ4ZQzdub10qVcjqa41SUi4iISytUqBDTnn0WvvwSWrVixpYteqiQSCbVt29fQkNDAfB4/nlsnjx8Wrkynm4qSfUTEBERl/fQQw/RLiwMQkPhhRfo2a+f2lhEMpk5c+bw2WefJSw0bEjsnXfyaunSBOXJ42wwF6GiXEREMoUPJkzAd/Jk8PXlyCOPqI1FJBPZv38/TzzxRMJCwYLkHDyY2/LkYVCJEs4GcyEqykVEJFMoVKgQHw8aBDNmQJMmzNi/nwXGQKlSMGuW0/FE5Ari4uLo0qULZ8+eBSDXa6+Rw9ubzypXJofaVpLoJyEiIpnGQw89xMN79sDWrfDss/QoWJDT+/ZB794qzEVc1JgxY/j1118BMA89RFTVqrxdtiwVvL0dTuZaVJSLiEimMnH/fgq89RZ4enJ0wAD6ADYyEoYOdTqaiFxi9erVvPzyywkLJUrg/tRT3FegAI8XLepsMBekolxERDKVQocOMfXAAfjwQ6hblzkPPJDw+Ob9+52OJiLJRERE0KlTJy5cuADu7uR+803yeXryccWKGGOcjudyVJSLiEjmUrIkDwE9vvsO/voLnnySV0qU4LMCBZxOJiLJ9OvXjx07dgDg0aMHEcWLM6ViRW7x9HQ4mWtSUS4iIpnLiBHg7c0ka2kwejScPw9Dh/JYeDhLly51Op2IAN9++y0fffRRwkLlysS1b0/XwoVp7efnbDAXpqJcREQyl06dYMoUPAICmH/yJCUnTICKFbnQuTOtW7dm8+bNTicUydYOHz5Mz549Exa8vPB5801KeHnxTvnyzgZzcSrKRUQk8+nUCUJDyWctyz/9FO+QEOjYkbMlS9K8eXP++ecfpxOKZEvx8fF069aNkydPApB7wAAi8udnRuXK5MuRw+F0rk1FuYiIZGolSpRg8YMPYo4cgSFD2H/yJC1btiQ8PNzpaCLZzrvvvstPP/2UsFC/PhFNmjCgRAka5s/vaK7MQEW5iIhkevVq1uTtAgXAzw+eeYY1a9bQvn37hFkfRCRDbNiwgUGDBiUsFCxIrpdfppaPD8NLl3Y2WCaholxERLKEfs2a0fLMGbjnHmjcmAULFvDMM89grXU6mkiWFxUVRceOHYmJiQFjyPPmm+DlxezAQHLqqZ3XRD8lERHJMv730EMUPXMG+vUDPz8mTZrE22+/7XQskSxv8ODBSTdZ5+jQgbAKFZhQrhwV9dTOa6aiXEREsgwPNzeW3X03Oby8YMgQcHPjhRde4Msvv3Q6mkiW9cMPP/Duu+8mLJQti+3Zk4cKFaLXLbc4GyyTUVEuIiJZSoXcuZlUqRIEBUHHjgB06dKFFStWOBtMJAs6duwY3bp1S1jw9MRn1CgKe3rykZ7aed1UlIuISJbTo3hx2uTPD926QWAg58+f58EHH2Tnzp1ORxPJMqy19OzZk6NHjwKQq39/wgsWZEalShT08HA4XeajolxERLIcYwwfV61K8Zw5cXv1Vcidm5MnT3Lfffdx/Phxp+OJZAkffvgh8+fPT1i44w6i7rmHASVK0LRAAWeDZVIqykVEJEvKlyMHX1avjvH3x33AAAB2797Ngw8+SFRUlMPpRDK3bdu20b9//4SFQoXweuUVamr6w5uiolxERLKsuvny8Xrp0sQ1agTNmgGwcuVKunTpQnx8vLPhRDKpmJgYOnbsmPCPWzc3vN98E/dcufg8MBBPTX94w/STExGRLG1QyZI0yp+fnAMGQPHiAHz99dcMHDjQ4WQimdPLL7/M2rVrAXDv1o3I8uX5oEIFTX94k1SUi4hIluZuDDMrV8bH0xP/996DxBvQ3n77bSZOnOhwOpHMZdmyZYwZMyZhoUYN4jt3pkvhwjxapIizwbIAFeUiIpLlFfP0ZFqlShzLn5+yo0YlrX/mmWf+vVFNRK7q1KlTdOnSJeEpuXnz4vn665Tz9uaD8uWdjpYlqCgXEZFs4cFChXi6WDF216xJ+e7dAYiPj6d9+/asXr3a4XQirs1ayxNPPMGhQ4cA8Hj5ZWy+fMwNDMQnRw6H02UNKspFRCTbGFO2LLV9fDjWrRslbr0VgMjISFq2bEloaKiz4URc2IwZM/59Mm6bNsQGBzO2bFlq5snjbLAsxKWKcmNMAWPMT8aYnYlffa+wX5wxZl3ia16y9aWNMX8aY3YZY+YaY3JmXHoREXF1nm5uzK1SBWsMBcaPx9fPD4CjR4/SvHlzTp8+7XBCEdeze/dunn766YSFChVwe/JJHixYkL7FijkbLItxqaIcGAwssdaWB5YkLqcmyloblPh6INn6UcB4a2054DTQI33jiohIZlM2Vy6mVqzI+thY7v3yS3LmTLh+s3XrVlq3bs358+cdTijiOmJjY+nUqRPh4eGQOzceb7xB0Vy5mFapEsYYp+NlKa5WlD8IzEj8fgbQ6loPNAm/Mu4CvrqR40VEJPto5+9Pn6JFmWMtz86Zk7Q+JCSEnj17JtzIJiIMHz6cP//8EwAzcCDx/v7MqVKFAomzGEnacbWivLC19p/E748Aha+wn5cxZrUx5g9jTKvEdQWBM9baC4nLBwH9v4qIiKTq7bJlqenjw8d+fgx6++2k9TNnzuSVV15xMJmIa/j9998ZPnx4wkKbNtg772RkmTLUy5fP2WBZlMnoqwHGmJ+B1CazHArMsNbmT7bvaWvtZX3lxphi1tpDxpgywFKgCXAW+COxdQVjTAlgkbW26hVy9AZ6AxQuXLj2nGRXSjJKeHg4Pj4+GX5e+ZfGwHkaA+dl5zE4RMJfBKWspdSECSycl3SbEi+88ALNmzfPkBzZeQxchcYgpYiICHr27MmRI0egcmV47z3ucHNjuDGkV9NKdhiDxo0b/22tDU5tW4bPYWOtbXqlbcaYo8aYW6y1/xhjbgGOXeE9DiV+3WOMCQFqAl8D+Y0xORKvlhcn4c/bK+WYAkwBCA4Oto0aNbrBT3TjQkJCcOK88i+NgfM0Bs7L7mOQ49gx2m/ZQrOxY7GxsSxatAiAcePG0bRpU+655550z5Ddx8AVaAxSevTRRxMK8rx5McOGUdzTkwV16pA/HdtWsvsYuFr7yjyga+L3XYHvLt3BGONrjPFM/L4QUA/YYhMu+S8D2l7teBERkeQeSewvH3foEJ0mTyYoKAiAuLg42rZty/r1650NKJLB5syZw2effQbGwIsv4l6oEN/UqJGuBbm4XlE+ErjbGLMTaJq4jDEm2BgzNXGfysBqY8x6EorwkdbaLYnbBgH9jTG7SOgx/zhD04uISKY0rlw5bsuThyf37eOdb76hePHiAISFhdGiRQsOHjzocEKRjLF//36eeOKJhIWOHaFuXd6rUIHamo883blUUW6tPWmtbWKtLW+tbWqtPZW4frW1tmfi979ba6tZa2skfv042fF7rLW3WWvLWWvbWWs1r5WIiPwnTzc3vqpShZzG0OfYMb5asIC8efMCcOjQIVq2bMm5c+ccTimSvuLi4ujSpQtnz56FoCB47DHa+fryeNGiTkfLFlyqKBcREXFKCS8vPg8MZEtkJO/kyMFXX39NjsTHh69fv56HH36Y2NhYh1OKpJ8xY8bw66+/QqFC8PLLlHJzY1qVKpqPPIOoKBcREUl0d4ECvFG6NJ8fO8a2SpWYMmVK0rYff/yRPn36aA5zyZJWr17Nyy+/DB4e8NpreOTLx4Jbb8UnR4bPCZJtqSgXERFJ5sWSJbm/YEH6795NxdatEwqVRFOnTmXkyJEOphNJexEREXTq1IkLFy5A375QpQqfValCYO7cTkfLVlSUi4iIJONmDJ9WqkRJT0/abd7MU0OH0qVLl6TtQ4YMYfbs2Q4mFElb/fv3Z8eOHdC8OTzwAI/nzcsjRVJ7pIykJxXlIiIil8jv4cHXVapw6sIF2m/dyqQpU2jcuHHS9u7duyf03opkct9++21Cm1alSvDss1SJimJizZpOx8qWVJSLiIikIihPHqZUqEDImTO8eOAA//vf/wgMDAQgJiaGVq1asW3bNodTity4w4cP07NnT/D1hddfJ3d0NCFNmuCuGzsdoaJcRETkCroUKUL/4sV579Ah/hcVxcKFCymS+N/6p0+f5r777uPo0aMOpxS5fvHx8XTv3p2TZ87Aq69i8uZlYXAwhXLmdDpatqWiXERE5CpGlSnD3b6+PLljB4fz5+f7778nd+INcKGhodx///1EREQ4nFLk+rz77rssXrwYnngCatTgRU9P7tR85I5SUS4iInIVOdzcmBMYSHFPT1pv3kyRqlWZM2cObm4Jf4WuWrWKjh07EhcX53BSkWuzYcMGBg0aBPfcA23bUjs0lBGNGjkdK9tTUS4iIvIfCnh4MK9aNcLj4nho0yaa3ncf7733XtL2efPm0a9fP81hLi4vKiqKTp06EVO+PDz/PHl27uSX9u2djiWoKBcREbkmVXLn5rNKlVgVFsYTO3bw5JNPMmDAgKTt7733HhMmTHAuoMg1GDx4MJuOH4c33sAcP85P9euT28vL6ViCinIREZFr1srPj2GlSjHj6FEmHDzIqFGjaNeuXdL2559/nq+//trBhCJX9sMPP/DulCkwYgR4ePBqTAx1Kld2OpYkUlEuIiJyHV4KCKB1oUIM2L2bH06f5tNPP+WOO+4AwFpL586dWblypcMpRVI6fvw4Xbt3h6FDoVQp6v74I6906+Z0LElGRbmIiMh1cDOGTytXJsjHh0e2bGH7hQt89913lC9fHoDo6GgeeOABdu3a5XBSEWDWLGxAAD38/TnWsiXUr0+ezz5j3iuvYDQfuUtRUS4iInKdcru7M79aNfK5u9Ny40Zi8+Rh4cKFFCpUCIATJ07QvHlzTpw44XBSydZmzYLevflw/37mN20KnTrB/Pl8Ubgwfn5+TqeTS6goFxERuQFFPT35vlo1TsfGcv/GjdxSujTz5s3DK/GmuZ07d9KqVSuio6MdTirZVfyQIYyOjOSpwEB44QVYu5a+77xDs7lznY4mqVBRLiIicoOC8uRhbpUqrA0Pp/PWrdSpW5eZM2cmtQWsWLGCrl27Eh8f73BSyW5OnDjB/fv3M6hYMeLffBOOHaPSa68xJi4O9u93Op6kQkW5iIjITWhRsCDjy5Xj2xMnGLxnD23atGHs2LFJ27/44gtefPFFBxNKdrNixQpq1qzJwnz5YNQosJagQYP46dw5vABKlnQ6oqRCRbmIiMhNerpYMZ4qWpQxBw7w0eHD9OvXj6eeeipp++jRo5k8ebKDCSU7iI+PZ/To0TRs2JCDx4/Dm29CoUJ0GjKEvw4fpjiAt3fClIjiclSUi4iI3CRjDBPKleO+AgV4cscOFp06xTvvvMP999+ftM9TTz3FggULHEwpWdmJEye4//77GTRoEHHWJkx9WKkSQzZuZGZkJB7GQEAATJmScMOnuBwV5SIiImkgh5sbcwMDCfLxod3mzayOiODzzz8nODgYSLiK+cgjj/D33387nFSymhUrVhAUFMTChQsTVjz1FDRowGt+fowYMABCQyE+PuGrCnKXpaJcREQkjeTJkYOF1atzS86ctNiwgYPGMH/+fEqVKgVAREQELVu2ZN++fc4GlSwhPj6eUaNG0bBhQw4dOpSwsl07aN2a54oW5dVq1ZwNKNdFRbmIiEga8s+Zkx9r1MDdGO5dv554X18WLlxI/vz5AThy5AgtWrTgzJkzjuaUzO3EiRO0bNmSwYMHExcXB0Duli2hTx/a+vnxduLDrCTzUFEuIiKSxsrmysWi6tU5eeECzTZs4JZy5fj222/x8PAAYPPmzbRp04aYmBiHk0pmdLFdZdGiRUnrKj/6KDHPP0/9fPn4rFIl3PS0zkxHRbmIiEg6qJUnD99UqcK2yEge3LSJOvXr88knnyRtX7p0Kb169cJa62BKyUxSbVcBOr35Jvt69KCKjw/zq1bFy93dwZRyo3I4HUBERCSralqgAJ9WqkSHrVvptHUrX3TsyN69e3n55ZcB+PTTTyldujSNGjVyNqi4vBMnTvDoo4+muDpeoEABXps5k5d9fCiWMyc/VK9O/sT/jZHMR1fKRURE0lH7woWZUK4c/ztxgp7bt/PikCH06NEjafuwYcP44YcfHEwori61dpXbb7+db/76ixF585InRw5+qlGDwjlzOphSbpaKchERkXT2bPHiDCtViulHjvD0rl188MEH3HPPPUnbx44dy88//+xgQnFFyR8GlLxdZeDAgcxevJhHjx0jDvipenUCvLycCyppwqWKcmNMAWPMT8aYnYlffVPZp7ExZl2yV7QxplXitunGmL3JtgVl9GcQERFJzcsBAQwqUYJJhw8zeP9+vvjiC6pXrw5AXFwcbdq0YePGjQ6nFFeR4mFAibOrFChQgO+//54X3niD5lu2cOrCBX6oXp1KuXM7nFbSgksV5cBgYIm1tjywJHE5BWvtMmttkLU2CLgLiAQWJ9vlhYvbrbXrMiCziIjIfzLG8FaZMjxTrBjjDx5kzMmTLFiwgGLFigFw7tw5WrRoweHDhx1OKk677GFAwB133MG6deu44557uHfDBvZGRzO/WjVq58njYFJJS65WlD8IzEj8fgbQ6j/2bwssstZGpmcoERGRtGCMYUK5cvS65RZG7N/PjLg4FixYgLe3NwAHDhygRYsWhIWFOZxUnHCl2VUGDhxISEgIPkWK0HT9ejZFRPC/KlVomDj3vWQNrlaUF7bW/pP4/RGg8H/s3x74/JJ1I4wxG4wx440xnmmeUERE5CYYY5hUoQKdCxfmpb17WVKgAK+++iruidPYrVu3jvbt23PhwgWHk0pGSu1hQBfbVUaNGkU4JBXk31atyn0FCzobWNKcyej5UY0xPwNFUtk0FJhhrc2fbN/T1trL+soTt90CbACKWmtjk607AuQEpgC7rbWvX+H43kBvgMKFC9eeM2fODX+mGxUeHo6Pj0+Gn1f+pTFwnsbAeRoDZ8QBbwC/AL3Onyffzz8zduzYpO33338//fr1w+ghMBnCyd8HGzdu5PXXX+fEiRNJ66pUqcIrr7yCv78/YcAAYC/wOlDXkZTpLzv8WdS4ceO/rbXBqW601rrMC9gO3JL4/S3A9qvs+yww5SrbGwHfX8t5a9eubZ2wbNkyR84r/9IYOE9j4DyNgXPOx8XZhzZutCxbZt8KDbVDhgyxQNJr5MiRTkfMNpz4fRAXF2ffeust6+7unmLcBw4caGNiYqy11p6KibG1V62yOUNC7IITJzI8Y0bKDn8WAavtFepRV2tfmQd0Tfy+K/DdVfbtwCWtK4lXyjEJlxVaAZvSPqKIiEjayOnmxtzAQO4CXty7lxw9e9KhY8ek7YMHD8aJ/8mV9HfixAlatGjBiy++mKJdZcGCBYwaNQoPDw9Ox8Zyz4YNbIyI4H9Vq9JcLStZmqsV5SOBu40xO4GmicsYY4KNMVMv7mSMKQWUIOF//ZKbZYzZCGwECgHDMyK0iIjIjfJwc2MI0K1IEV7ft4+ir77KnQ0bJm3v2rUrv/32m3MBJc0tX76coKCgFA+Nuji7SvPmzQE4HhPD3evXsyE8nK+rVKGFCvIsz6WKcmvtSWttE2tteWttU2vtqcT1q621PZPtF2qtLWatjb/k+LustdWstVWttZ2tteEZ/RlERESulzvwccWKPFm0KG8fPkyl99+nYuXKAMTExNCqVSu2b9/ubEi5afHx8YwcOZJGjRqlmF1l0KBBhISEUKJECQAOREfTYO1aNkdG8k3VqrQsVMipyJKBXKooFxERya7cjGFi+fL0K16cKSdOUPPTT/ErnDAJ2alTp2jevDnHjh1zOKXcqKu1q4wcORIPDw8AtkdGUm/tWv6JiWFx9epqWclGVJSLiIi4CGMMb5cty9CSJZkTHk7Ql1+SK29eAPbs2cMDDzxAZKQezZHZXEu7CsDfYWHUX7uW8/HxhAQF0UDzkGcrKspFRERciDGG4WXKMKpMGX6Ki6Ps119jEp/a+Oeff9K5c+ekK63i2q7UrnLxYUAX21UAQk6fpvG6deR2c2N5zZrU1JM6sx0V5SIiIi5oYMmSfFapEts9PCjyxReQ2Ff8zTff8MILLzicTv7L8ePHU21XufgwoIvtKgDfnThBsw0bKOHpyYpatSif+IRXyV5UlIuIiLiozkWKsLBaNcJz58bn00+hVCkAxo8fz3vvvedsOLmi3377jZo1a6ZoV7n99ttZt24dLVq0SLHvewcP0nrTJmr4+PBrzZoU89TDyLMrFeUiIiIurGmBAvwaFESefPnwmDwZatQA4Nlnn+W77672OA/JaPHx8bz11ls0btz4snaVX375JUW7yoX4ePru2MEzu3bRsmBBltSoQcFkV88l+1FRLiIi4uKC8uRhZa1alM2XDzN2LNx1F9ZaOnTowF9//eV0POHfdpUhQ4b8Z7vKuQsXuH/TJiYePszzxYvzv6pV8cmRw6no4iJUlIuIiGQCAV5erKhdmzp588LLL0OvXkSdP8/999/P3r17nY6XraXWrnJxdpVL21VCo6K4Y80afj59mikVKjC2XDncjcnoyOKCVJSLiIhkEgU8PPglOJhHvL2hY0d4802ORUZy3333cerUKafjZTtXa1e5dHYVgD/OnqXOmjUcPH+eH6pXp1fRohkdWVyYinIREZFMJKebG3Nuu40Bbm5QqxZMmsT28+d5qGhRzhuTcDPorFlOx8zyrtSusmDBgsvaVay1fPzPPzRatw4fd3f+qFWLJr6+TkUXF6WiXEREJBMac+edvB4eDt7eMHEivwYH0w04vW8f9O6twjwd/fbbb9f0MCCAyLg4um/bRs/t22mQPz9/1qpFpdy5MzqyZAIqykVERDKplx96iCEDB8K+fTB8OHO6d6egmxu3RkYyqE8fFi9erCeApqHk7SqHDx9OWj9o0KBU21V2REZSd80aPj16lFcDAvihenUK5cyZ0bElk9CtviIiIpnY8N27Ofbss0zt1w8efRQbFMTqESNYfewYo++9Fw8PD26//XaaNGlCkyZNuO2221K0Vsi1OX78OI8++miKq+MFCxbk008/vezqOMCXx47RY/t2chrDourVubdAgYyMK5mQrpSLiIhkYiYggEmxsYwfPZoyI0ZA2bIwdSrceScAsbGx/Prrr7z66qvUr18fX19fmjdvzttvv83atWuJj493+BO4viu1q6xdu/aygvx8fDzP7dzJw1u2UCV3btYGB6sgl2uiolxERCQzGzGCHN7ePAfs/vln1vTuTflDh2DYMPIPHw5eXil2j4iIYNGiRQwYMIBatWrh7+9Pu3btmDx5Mjt27MBa68jHcEXX266yPjyc2/7+m3cOHeK54sX5JSiIEpf8/EWuRO0rIiIimVmnTglfhw6F/fup6eHB5jx5eLVkSUYC5ZYs4bFDh9j9448sWbKE0NDQFIefPHmSr776iq+++gqA4sWL06RJE+666y6aNGlCsWLFMvbzuIjraVe5EB/PqAMHGBYaSkEPD+ZXrUrLQoUyOrJkcirKRUREMrtOnf4tzgEP4E2gqa8vXbZu5VV/fwYNHcrWDz/kn/37WbJkCUuWLGHp0qUcO3YsxVsdPHiQGTNmMGPGDAAqVqyYVKA3atSIggULZuAHc8Zvv/1G+/btU1wdr1evHnPmzKF48eIp9t0aEUHXbdtYFRZGe39/3i9fnoLq2ZcboPYVERGRLOouX1823Hor7f39Gb5vH0GrV3PA15eePXvy+eefc+TIETZu3MiECRO4//77yZs372XvsX37diZNmkTbtm3x8/OjVq1aDBgwgEWLFhEeHu7Ap0o/V2pXGTx4MMuWLUtRkMdby7gDB6i5ejV7oqL4IjCQzwMDVZDLDdOVchERkSysoIcHn1auTJfChXl8xw4arltHr1tuYXSZMuT38KBq1apUrVqVZ599lgsXLvD3338nXUlfsWIF58+fT3ovay1r165l7dq1vP322+TIkYM6deokzexSp04dPD09Hfy0N+7MmTM0b96cH3/8MWldwYIF+eyzz7jvvvtS7PvnuXP03bmT1WFhPFCwIFMqVqSwpjqUm6SiXEREJBu4u0ABNt16K6+FhjLuwAHmnzzJu+XK0dbPD2MMQFKRXadOHYYMGUJ0dDS///57UqvLqlWrkp5eCXDhwgVWrFjBihUreP3118mVKxcNGjRI6kmvWbMm7u7uTn3ka/bbb7/Rq1cvTpw4kbQutXaVozExvLhnD58cOULRnDmZVbkyHfz9k35+IjdDRbmIiEg24e3uzuiyZeng70/P7dt5eMsWGuTLx5iyZamTSuuKl5cXd911F3fddRcAZ8+e5ddff2Xp0qUsWbKEjRs3ptg/KiqKxYsXs3jxYgB8fX1p1KhR0pX0ihUrulQBGx8fz6hRo3jppZdSTA05ePBg3njjDXLkSCiTYuPjmXjoEK+GhhIVH8+gEiUYGhBAnhwqoyTt6FeTiIhINlMzTx7+rFWLaUeO8GpoKHXXrKGtnx9vli5NeW/vKx6XL18+7r//fu6//34Ajh07xrJly5KupO/evTvF/qdPn+abb77hm2++AaBo0aJJRX6TJk0oWbJk+n3I/3D8+HG6dOly1XYVay0/njrFgN272RwZSbMCBZhQrhwVr/IzErlRKspFRESyoRxubvQuWpSO/v6MO3iQMQcO8O2JEzxRtCgvBwTgfw090v7+/jzyyCM88sgjAOzbty+pQF+yZAlHjhxJsf/hw4eZOXMmM2fOBKBcuXJJBXrjxo3x8/NL+w+ail9//ZUOHTqkuJmzatWqLFq0iOLFi2OtZfHp07wWGsof585R2suL76pW5f6CBV3qSr9kLSrKRUREsjGfHDl4pVQpHi9alNdDQ5l06BDTjxzhiaJFebZYMYpfx8NvAgICeOyxx3jsscew1rJ161aWLl3K0qVLWbZsGWfOnEmx/65du9i1axdTpkwBoHr16kmtLnfeeSd58uRJy4961XaVpk2bUqxYMX48dSqpGA/w9GRKhQp0LVKEnG6asE7Sl4pyERERoXDOnEysUIFnihdnWGgo4w8cYMLBg3Tw9+f5EiWo4eNzXe9njCEwMJDAwED69u1LXFwca9euTbqK/ttvvxEVFZXimA0bNrBhwwbGjx+Pu7s7t912W1KRXrduXbxu4umYx48fp3Pnzkn97vBvu8o9zZox6pdfeGntWv44d46SKsbFASrKRUREJElFb29mBwbyZunSvHPoEB8dPsxnR49yt68vA0qU4G5f3xtq4XB3dyc4OJjg4GAGDhzI+fPn+eOPP5KK9D///JMLFy4k7R8XF8fKlStZuXIlw4cPx8vLi/r16ye1u9SuXfuaZ3ZJrV2lXr16vDNzJouAx//4gwNAyfPn+bBCBbqpGBcHqCgXERGRy5TKlYvx5crxSkAAHx4+zLuHDnHvhg2U9vKic+HCdC5cmAo3ccOjp6cnDRs2pGHDhgwbNoywsDCWL1+eNEf6unXrUuwfHR3Nzz//zM8//wwk3HTasGHDpCvpgYGBl/1jIT4+npEjR/Lyyy+naFd5ZNQoYu+7j7r79nHBWpr6+tLz/HlerFMHDxXj4hAV5SIiInJFvh4eDA4IoF+JEnxx7BifHT3KiH37eGPfPm7Lk4fOhQvziL//Nd0YejV58uThvvvuS5r55MSJE4SEhCTdOLpjx44U+589e5Z58+Yxb948AAoXLvzvzC5hYfiMG0fngwdZDODmBlWr4nX33eRv2ZK5bm4UOHuWZ4sV4/GiRSnv7U1ISIgKcnGUinIRERH5T55ubnQpUoQuRYpw+Px5Pj92jJlHj/LMrl3027WL+vny0dTXl6a+vgTnyUOOmyxwCxUqRNu2bWnbti0ABw4cSLppdMmSJRw6dCjF/kePHuXzzz/n888/B8DDw4PYunWhfn244w7w9SUeqFmgAB38/Wnn54dXJniwkWQfLlWUG2PaAa8BlYHbrLWrr7BfM+AdwB2Yaq0dmbi+NDAHKAj8DXSx1sZkQHQREZFso6inJ8+XKMHzJUqwKTyc2ceO8cOpU7wcGsrLoaHkc3enUf78NPX1pV6+fFTy9ibXTRbAJUqUoGvXrnTt2hVrLTt27Egq0JctW8YpNzcIDEx6xVaoAJ6eEBFB5T/+4KVt22j5+efk1QN/xEW52q/MTUBr4MMr7WCMcQcmAncDB4FVxph51totwChgvLV2jjFmMtADmJT+sUVERLKnqj4+vOnjw5tlynA8JoZlZ87w8+nT/Hz6NN+dPAmAG1A2Vy6q5M5N1dy5qeLtTblcuSjo4UFBDw/yuLtf082jZ2Jj2Xf+PKHR0ezLnZt9d9/Nhfr18erfH2ISrsG5xcTAjh3Ez5tHnlWrmLluHQ/ExoIxoIJcXJhL/eq01m4F/us35m3ALmvtnsR95wAPGmO2AncBHRP3m0HCVXcV5SIiIhnAL2dOHvb352F/fwD2REXxd1gYmyIi2BwRwebISOafOEHcJcflMIYCOXJQ0MMDH3d3Yq0lJj6emGRfw+PiCI9LeaSXmxulvLxonD8/dfLmpW7evNSoWRP27GErUBrIe3FnB58eKnItjLXW6QyXMcaEAANSa18xxrQFmllreyYudwHqkFCA/2GtLZe4vgSwyFpb9Qrn6A30BihcuHDtOXPmpMMnubrw8HB8rnPeV0lbGgPnaQycpzFwXnYagxhgP3AUOAeEAWeTfR9FwhVDj2RfPQBPwA8onOyVH7j0Mp7/zz9TcexY3M+fT1oX5+nJ9gEDONa06RVzZacxcFXZYQwaN278t7U2OLVtGX6l3BjzM1AklU1DrbXfZVQOa+0UYApAcHCwbdSoUUadOklISAhOnFf+pTFwnsbAeRoD52kM0lCjRlC5MgwdCvv3Q8mSuI8YQWCnTgRe5TCNgfOy+xhkeFFurb3yP1OvzSGgRLLl4onrTgL5jTE5rLUXkq0XERGR7KRTp4SXSCaSGSfkXAWUN8aUNsbkBNoD82xCH84yoG3ifl2BDLvyLiIiIiJyo1yqKDfGPGSMOQjcDiwwxvyYuL6oMWYhQOJV8L7Aj8BW4Atr7ebEtxgE9DfG7CJhWsSPM/oziIiIiIhcL1ebfeUb4JtU1h8GmidbXggsTGW/PSTMziIiIiIikmm41JVyEREREZHsSEW5iIiIiIjDVJSLiIiIiDhMRbmIiIiIiMNUlIuIiIiIOExFuYiIiIiIw1SUi4iIiIg4TEW5iIiIiIjDVJSLiIiIiDjMWGudzuA4Y8xxYJ8Dpy4EnHDgvPIvjYHzNAbO0xg4T2PgPI2B87LDGARYa/1S26Ci3EHGmNXW2mCnc2RnGgPnaQycpzFwnsbAeRoD52X3MVD7ioiIiIiIw1SUi4iIiIg4TEW5s6Y4HUA0Bi5AY+A8jYHzNAbO0xg4L1uPgXrKRUREREQcpivlIiIiIiIOU1HuEGNMM2PMdmPMLmPMYKfzZDfGmGnGmGPGmE1OZ8mujDEljDHLjDFbjDGbjTHPOp0puzHGeBlj/jLGrE8cg2FOZ8qOjDHuxpi1xpjvnc6SHRljQo0xG40x64wxq53Okx0ZY/IbY74yxmwzxmw1xtzudCYnqH3FAcYYd2DH/9u7v5e95ziO4883Qza/DrCWW9mBHHBgWhZjRPNz4QzFgdROEFGK1P4DCeVkIzKWzEqZHyvym9akxCSN2r3oVlpMSubl4P4e3EfcK11v9/19Purqur7fo1ddXV2vPtf787mA9cA0sBu4NclXrcFGpKrWAYeA55Kc151njKpqBbAiyWdVdSKwB7jJz8HkVFUBy5IcqqpjgA+Ae5N80hxtVKrqfmA1cFKSDd15xqaqvgdWJ1ns52P/b1XVs8D7STZX1bHA0iQHm2NNnCvlPS4Evk2yL8kfwDbgxuZMo5LkPeDn7hxjluSHJJ8Nr38F9gJn9KYal8w6NFweMzxcqZmgqpoCrgc2d2eROlTVycA6YAtAkj/GWMjBUt7lDGD/nOtpLCMasao6C1gFfNocZXSG0YnPgRlgVxLfg8l6DHgQ+Ks5x5gFeKuq9lTVxu4wI7QS+Al4Zhjj2lxVy7pDdbCUS2pVVScA24H7kvzSnWdskhxOcj4wBVxYVY5zTUhVbQBmkuzpzjJylyS5ALgWuGsYb9TkLAEuAJ5Ksgr4DRjlXjtLeY8DwJlzrqeGe9KoDHPM24GtSV7pzjNmw8/F7wDXNEcZk7XADcNM8zbgiqp6vjfS+CQ5MDzPADuYHTHV5EwD03N+pXuZ2ZI+OpbyHruBs6tq5bCh4Rbg1eZM0kQNmwy3AHuTPNqdZ4yq6rSqOmV4fTyzm8+/bg01IkkeSjKV5CxmvwfeTnJbc6xRqaplw0ZzhpGJqwBP5ZqgJD8C+6vqnOHWlcAoN/wv6Q4wRkn+rKq7gTeBo4Gnk3zZHGtUqupF4HLg1KqaBjYl2dKbanTWArcDXwwzzQAPJ9nZF2l0VgDPDidCHQW8lMRj+TQmy4Eds2sELAFeSPJGb6RRugfYOixU7gPuaM7TwiMRJUmSpGaOr0iSJEnNLOWSJElSM0u5JEmS1MxSLkmSJDWzlEuSJEnNLOWSJElSM0u5JEmS1MxSLkmSJDWzlEuS5qWqLquqVNV1c+6trKqZqnq8M5skLXT+o6ckad6q6m3guCRrq+pk4CPgO+DGJId700nSwmUplyTNW1VdCrwHXA08ACwHLklyqDWYJC1wlnJJ0hGpql3AxcBBYE2S6d5EkrTwOVMuSTpS3wJLgU0Wckn6b7hSLkmat6raCDwB7AV+T3JRcyRJWhQs5ZKkeamq9cBO4E7gG+Bj4Lokr7cGk6RFwFIuSfpXVXUu8CHwZJJHhnu7gJOSrGkNJ0mLgKVckvSPqup04FNgN3Bzhi+OqloHvAtsSPJaY0RJWvAs5ZIkSVIzT1+RJEmSmlnKJUmSpGaWckmSJKmZpVySJElqZimXJEmSmlnKJUmSpGaWckmSJKmZpVySJElqZimXJEmSmv0NvGro+6fsIlwAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(df.X,df.Y,'ro')\n",
"plt.plot(df.X,f(df.X),'b-')\n",
"plt.plot(5,f(5),'y*')\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "BP0sX_6fEZXA",
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"### Polynomial object in `numpy`\n",
"In `numpy` it is possible to define polynomials friom either its coefficients o its roots with `np.poly1d`"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "yKOPE5AcEZXC"
},
"source": [
"Define a two degree polynomial from its roots:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "S_V2zHPnjz19"
},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "l-sKCHWojz2A"
},
"source": [
"We try to make a fit by using an inverted parabola passing trough the tree points and using the roots as a guess.\n",
"In fact, we can try with a polynomial of degree two with roots at 1 and 22.\n",
"\n",
"With $k$ we can flip the curve and reduce the maximum without change the roots"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "uV6xujqPEZXC"
},
"outputs": [],
"source": [
"k=-0.1\n",
"P=k*np.poly1d([0,23],r=True)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"colab_type": "code",
"id": "9ztBor9-EZXE",
"outputId": "44ab36ca-8f74-4113-aec0-e8d6072beab0"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2\n",
"-0.1 x + 2.3 x - 0\n"
]
}
],
"source": [
"print(P)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"colab_type": "code",
"id": "oGBYy5vAm6LU",
"outputId": "b2e42386-2d86-4092-bcfe-c5088642f126"
},
"outputs": [
{
"data": {
"text/plain": [
"array([23., 0.])"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P.roots"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 286
},
"colab_type": "code",
"id": "PUM8A4EmEZXG",
"outputId": "413d70fc-0b17-497a-a84f-2b96300bc306"
},
"outputs": [
{
"data": {
"text/plain": [
"(-15.0, 20.0)"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlOUlEQVR4nO3deZgU1dn+8e/DIrLI5jJBVmNQQ1DxHcRdATeieIGJaxCVaMAoSgQVFY3EfSMY4xJBXIjgaNzFLWgG0fdn1AFBEBRRwaCgEVwYN5R5fn+c5nXAAWborj693J/rqqunq5e6LZunq0+dOsfcHRERKS71YgcQEZHsU/EXESlCKv4iIkVIxV9EpAip+IuIFCEVfxGRIpR28Tez9mZWbmbzzOwNMxuWWt/azKaa2dup21bpxxURkUywdPv5m1kboI27zzSzLYAZQH/gZGCFu19tZucDrdx9ZJp5RUQkA9I+8nf3pe4+M/X3SmA+0BboB9ydetrdhC8EERHJAWkf+a/1ZmadgOlAV+B9d2+ZWm/Ap2vur/OawcBggMaNG5e2b98+Y3k2pKqqinr18ueUh/ImS3mTpbzJWrBgwSfuvnWdXuTuGVmAZoQmn1+l7n+2zuOfbuw9SktLPVvKy8uztq1MUN5kKW+ylDdZQIXXsWZn5KvNzBoCDwKT3P2h1OqPUucD1pwX+DgT2xIRkfRlorePAROA+e7+52oPPQaclPr7JODRdLclIiKZ0SAD77EPMBCYY2azUusuBK4G7jezU4DFwDEZ2JaIiGRA2sXf3V8EbD0PH5ju+4uISOblz+lsERHJGBV/EZEipOIvIlKEVPxFRIqQir+ISBFS8RcRKUIq/iIiRUjFX0SkCKn4i4gUIRV/EZEipOIvIlKEVPxFRIqQir+ISBFS8RcRKUIq/iIiRUjFX0SkCKn4i4gUIRV/EZEipOIvIlKEVPxFRIqQir+ISBFS8RcRKUIZKf5mdoeZfWxmc6utG21mH5jZrNRyWCa2JSIi6cvUkf9dQJ8a1o91926p5ckMbUtERNKUkeLv7tOBFZl4LxERSZ65e2beyKwTMMXdu6bujwZOBr4AKoAR7v5pDa8bDAwGKCkpKS0rK8tIno2prKykWbNmWdlWJihvspQ3WcqbrF69es1w9+51epG7Z2QBOgFzq90vAeoTfl1cAdyxsfcoLS31bCkvL8/atjJBeZOlvMlS3mQBFV7Hmp1Ybx93/8jdV7t7FTAe6JHUtkREpG4SK/5m1qba3SOBuet7roiIZFeDTLyJmd0L9AS2MrMlwCVATzPrBjiwCBiSiW2JiEj6MlL83f34GlZPyMR7i4hI5ukKXxGRIqTiLyJShFT8RUSKkIq/iEgRysgJX5FC5g5ffw0rVvywfPMNNGny46VxY2jWDMxipxbZMBV/EaCqChYuhNdfD8ucObBgASxfHor9t9/W/r2aNIEddqh5adUquf8GkbpQ8Zei9NlnMHUqPPsszJ4Ns2fvxzffhMfq1YPOnWGnnWDPPaF16x8vjRqFXwNffbX28uWXsHRp+OKYORMefBBWr/5hu9tvD717h6VXLygpifKfL6LiL8XBHebOhSefDMv//m8oyi1awP/8Dxx++FIOO6wdu+wCXbqEo/dMWLUK3nsvfBm8+WbY7v33w/jx4fEuXX74MjjkEGjaNDPbFdkYFX8pWO5QUQF33gmPPw5LloT13brByJFw2GGwxx7QoAFMm7aQnj3bZTzDZpvBjjuG5Ygj4Nxzw5fOa6/Bv/4VljvugJtuCucKjjoKTjwRDjgg/AIRSYqKvxSczz+HyZNh3DiYNSscxffpA6NHh9u2bePmq18funcPy3nnhV8HL74IkybBP/4Bd90FHTrACSfAwIGh+Ukk03RsIQXBHV5+GU45BbbdFk4/Pay/9dbQBv/gg+Gx2IW/JpttFpp9JkyAZcvg3nvhF7+Aq6+Gn/88nHe4//61zx2IpEvFX/KaOzz6KJSWhiJ5330wYAC8+mo44XraadC8eeyUtdekCRx3XDgvsWQJjBkTTk4fe2z4IpgwIfxSEEmXir/kJXeYMiU0nfTvDytXhqP8Dz8MzT3du+d/X/s2bWD4cHjjDXjgAdhiCzj11NBj6IEH2vLll7ETSj5T8Ze84g5PPx1O1B5xBHz6aTihO39+/h3l11b9+vDrX4eT108/DT/9Kdx8c2c6dYIrrwxdTEXqSsVf8sbzz8O++8IvfwkffRS6S771Fpx8cuixU+jM4NBDw3648caZ9OgBo0aF5qAHHghfjCK1peIvOe/jj0Ovl549YfHi0Lzz9tuhCaRhw9jp4th55y944gmYNg1atoSjj4YDDwzXMojUhoq/5KyqqtB+v+OO4UTuqFGh6J92WughI+F6gBkz4JZbwpXK3brBWWeF5jCRDVHxl5w0ezbssw8MGRIK2uzZcPnlYeA0WVuDBvD734eriIcMgZtvDsNTTJigpiBZPxV/ySmVlXDOOaHr5jvvwMSJ4SrYn/88drLct+WWofDPnBmuEzj1VOjbN1w7ILIuFX/JGRUVsNtuoW/7KaeEsXAGDsz/LpvZtuuuUF4Of/1r+OLs2hUeeih2Ksk1Kv4SXVVVKPh77x2GTn7+ebjttjB6pmyaevVg6NAwhlCnTqGr6KBB8MUXsZNJrlDxl6iWLQtdN885J/TbnzUL9t8/dqrCsdNO8NJLcPHFoQltl11g+vTYqSQXZKT4m9kdZvaxmc2ttq61mU01s7dTt5rGQtbyzDOhiWL6dPjb30JfdR3tZ17DhnDppWHwuAYNQpfZUaM0VlCxy9SR/11An3XWnQ885+6dgedS90VYtSoMbdynD2yzTWjrHzJEbftJ22uv8Mvqt78NVwb37asuocUsI8Xf3acDK9ZZ3Q+4O/X33UD/TGxL8tsnn8BBB8H114eRN195JfRMkexo1gxuvz380nruOdh99zB2kBQf8wx1BDazTsAUd++auv+Zu7dM/W3Ap2vur/O6wcBggJKSktKysrKM5NmYyspKmjVrlpVtZUIh5F20qAkXXrgzn3zSiJEj3+TAAz+OlO7HCmH/1tWcOc255JKufPNNPS644E322++TDKX7sWLcv9nUq1evGe7evU4vcveMLEAnYG61+5+t8/inG3uP0tJSz5by8vKsbSsT8j3vk0+6N2/u/pOfuP/733EybUi+799NtWSJe48e7uD+xz+6r16dkbf9kWLdv9kCVHgda3aSvX0+MrM2AKnb3DnMk6xxhxtuCO3LP/1paObZY4/YqWSNtm1D19pBg8JJ4f79w0xoUviSLP6PASel/j4JeDTBbUkOWrUqnMg9+2zo1y/0NmnfPnYqWdfmm4ehIG66CZ56KnS1Xbo0dipJWqa6et4LvATsaGZLzOwU4GrgYDN7GzgodV+KRGVlAw49NAy7fOGFoRtn06axU8n6mMEZZ8ATT4RhNfbdN9xK4crIKOjufvx6HjowE+8v+eXjj+EPf+jG+++HC4sGDoydSGrrkEPCkBCHHRYG1nvqqTDkhhQeXeErGfWf/4RmgyVLGvP44yr8+ahHj9BEt9lm4YKwadNiJ5IkqPhLxixcCPvtF9qLr732dQ49NHYi2VQ77QT/7/+FE8J9+sDDD8dOJJmm4i8ZMXduKPyVlaHZYJdd1GUk37VrBy+8EJp9jjoqXBwmhUPFX9JWURFmlDIL4/SUlsZOJJmy5Zbw7LNh7uDf/Q7Gjo2dSDJFxV+CSZPC2L/16oXbSZNq9bLp06F3b2jePLQTd+mSaEqJoGlTePTRcPQ/fHiYMEbyX0Z6+0iemzQJBg+Gr74K9xcvDvcBBgxY78umTw/twR07hqPDtm2zkFWiaNgQJk+G774L8wQ0bPjDR0Tyk478JYzvu6bwr/HVV2H9erz6arhqt1OncIWoCn/ha9gQ7rsPDj88XLx3552xE0k6VPwF3n+/TuvnzAltwFttBVOnhmGZpTg0ahQu2DvkkDDV5j33xE4km0rFX6BDh1qvX7AADj4YGjcOQwLriL/4bL45PPII9OoFJ50Ufg1I/lHxF7jiCmjSZO11TZqE9dW8/34Yi3/16tDGv912WcwoOaVxY3jssTAMxIABmiA+H6n4S/jXO25cOHNrFm7HjVvrZO+yZXDggWEC8H/+E37+84h5JSc0bQpTpoRRWo89NkzLKflDxV+CAQNg0SKoqgq31Qr/8uWhqWfpUo31ImvbYgt48skwG9tRR8Frr8VOJLWl4i8b9OWXYZCvt98Ofb332it2Isk1LVqEL4BWrcJnZfHi2ImkNlT8Zb1Wrw4/AF59NZzUO1BjtMp6bLtt+FX4zTfwy1/CinVn9Jaco+Iv63XeeeFo/4YbwmQsIhvyi1+EXkDvvBNmBPvmm9iJZENU/KVGt94Kf/4znHkmnHVW7DSSLw44IMzh8MILcOKJ4RSS5CYN7yA/8tRT4RL+vn01kJfU3bHHwpIlcM45YdrOMWNiJ5KaqPjLWmbPhmOOgV13hXvvhfr1YyeSfDR8eLgu5M9/Dl8A3brFTiTrUrOP/J8PPwxH+y1awOOPQ7NmsRNJvjILhf9XvwpfBP/+d+vYkWQdKv4ChC6dRxwBn30WJvHWsA2Srvr14e9/D0f9l1/ehQULYieS6lT8haoqOOEEmDUrdOncddfYiaRQNGkSpoBs0KCKfv3CFeKSG1T8hauuCl30xowJF+mIZFLHjjB69DzefhsGDlQPoFyRePE3s0VmNsfMZplZRdLbk7p55hm4+GL4zW9g2LDYaaRQdev2GWPHhsHgLr00dhqB7PX26eXun2RpW1JLixaFot+1axjHzSx2IilkQ4fCzJnwpz+F8wD9+8dOVNzU7FOkvv469MRYvToMx9u0aexEUujMwsWDu+8emn/mzYudqLhlo/g78E8zm2FmmvUzB7jD6aeHERjvuQd+9rPYiaRYbL75Dwcb/frBp5/GTlS8zN2T3YBZW3f/wMy2AaYCZ7r79GqPDwYGA5SUlJSWlZUlmmeNyspKmuVRR/ZM5n3ssW0ZO3YHBg5cxG9/uygj77muYt6/2ZDveefMac7w4d0oLf2UK6+cQ70ca4PIt/3bq1evGe7evU4vcvesLcBo4Jz1PV5aWurZUl5enrVtZUKm8r70knvDhu59+rh//31G3rJGxbp/s6UQ8t5yizu4X3NN9vNsTL7tX6DC61iPE/2+NbOmZrbFmr+BQ4C5SW5T1u+jj8KEG+3awaRJGrpB4jrttPB5vPBCeOml2GmKT9I/tkqAF81sNvAK8IS7P53wNqUGay7kWr4cHnwQWutqe4nMDMaPhw4d4Ljj1P6fbYl29XT3dwFdL5oDxowJk67fdpumYZTc0bIllJXBPvvAb38bTgary3F25NhpFknCq6+Gn9a//jX87nex04isrUcPuOaacJX5zTfHTlM8VPwL3MqVcPzx0KZN+ImtoyrJRWefDYcfDiNGhAvBJHkq/gVu6FB4771wgrdVq9hpRGpmBnfdBVtvHSaDWbkydqLCp+JfwCZPDlPqXXQR7Ldf7DQiG7bVVmECoXffDT2BEr4Eqeip+BeoNf+A9tknDNwmkvMmTWK/gZ34U9XFTJ4Mdw5R/88kqfgXoO++CwO21asXmnsaaLJOyXWTJsHgwbB4MRdwJb15jrPG78y7Yx+NnaxgqfgXoNGj4eWXw0idHTvGTiNSC6NGwVdfAVCfKu7iZOqzmpMuaMPq1ZGzFSgV/wIzbVqYnOWUU8JE7CJ54f3317rbniXcyFm8+G0Pxo6NlKnAqfgXkMpKGDQItt8e/vKX2GlE6qBDhx+tOpGJ9G/8DKNGwRtvRMhU4FT8C8h558HixaHLnMbnl7xyxRVhwt9qrEkTbrt+JS1ahPH/V62KlK1AqfgXiGefDRNlnH126OEjklcGDPjhJJVZuB03jm1OP4rbbgtzT1x+eeyQhUXFvwB88UVo499hB/0DkTw2YECYW7SqKtwOGADAkUeGI/8rrwxDlUhmqPgXgHPPhSVLQnNP48ax04hk3o03hiFKBg4MU5BK+lT889w//xl+LY8YAXvtFTuNSDJatoQ774S33oILLoidpjCo+Oexzz+HU0+FnXaCSy+NnUYkWQcdBGecEXqylZfHTpP/VPzz2IgR8MEHcPfdYWJskUJ3zTXws5+FocnV/JMeFf889fTTMGFC6N7Zo0fsNCLZ0bRpmJDonXfgT3+KnSa/qfjnoTXNPV26hKEcRIpJ795h1q/rr4dZs2KnyV8q/nnoggtg6dLQu6dRo9hpRLLvuutgyy1D84/G/tk0Kv555qWX4G9/g7POgt13j51GJI7WrUP3z4qKcCt1p+KfR777DoYMgbZt1btH5JhjwtSPF10UrgmTulHxzyNjx8KcOXDTTbDFFrHTiMRlBrfcEuat+P3vNfNXXSVe/M2sj5m9ZWYLzez8pLdXqJYt25zRo6F/f+jXL3YakdzQoUMYE+7pp8MUkFJ7iRZ/M6sP3Az8EugCHG9mXZLcZiFyhxtu6Ez9+mrfFFnXGWfAHnvAH/4Ay5fHTpM/kj7y7wEsdPd33X0VUAbouLWOHngAXn55Sy67DNq3j51GJLfUrw/jx8Onn4YLH6V2ki7+bYH/VLu/JLVOaunzz0PPns6dVzJ0aOw0Irlp553DBY933w3PPRc7TX4wT/AsiZkdBfRx91NT9wcCe7j70GrPGQwMBigpKSktKytLLE91lZWVNGvWLCvbSscNN3Tm8ce3ZcyYF+nWLX86NOfL/l1DeZOVjbyrVtVj0KDdadCgittvr6Bhw02vbfm2f3v16jXD3bvX6UXuntgC7AU8U+3+BcAF63t+aWmpZ0t5eXnWtrWpXnrJ3cx92LD8yFud8iZLeWs2ZYo7uF93XXrvk2/7F6jwOtbnpJt9XgU6m9l2ZrYZcBzwWMLbLAjffx/69G+7LVx2Wew0Ivnh8MOhb98w7s+HH8ZOk9sSLf7u/j0wFHgGmA/c7+6airkWbr0VXn899O5Rn36R2rvhhnBB5Lnnxk6S2xLv5+/uT7r7Du6+vbtfkfT2CsF//wt//CMcfHCYwk5Eam/77cPJ38mTYfr02Glyl67wzUEXXwwrV4YjGLPYaUTyz/nnhznghw4NTajyYyr+Oea118K0jGeeGYZsFpG6a9Lkh+FQbrkldprcpOKfQ9xDn/6ttoJLLomdRiS/9e8PhxwSfkl/9FHsNLlHxT+HlJXBiy/ClVeGCatFZNOZhQ4TX38dmoFkbSr+OeLLL0PvhNJSGDQodhqRwrDjjjB8eJj46KWXYqfJLSr+OeKqq8Jk7DfeGMYqEZHMuOiiMAfG0KGa9as6Ff8c8O67YT7SE06AvfeOnUaksDRrFv59zZwJd94ZO03uUPHPASNGQIMGcM01sZOIFKZjjw0HVhdfDJWVsdPkBhX/yKZOhUceCT9Nt902dhqRwmQGY8bAsmVw7bWx0+QGFf+Ivv8+TECx/fZw9tmx04gUtj33DL8Arr8eliyJnSY+Ff+IJkyAefPCh7FRo9hpRArfVVeFk76jRsVOEp+KfyQrV4YLufbdV3PyimTLdtuFX9sTJ4YTwMVMxT+SMWPCVYfXX6/xe0Sy6cILw1X0I0aEq+qLlYp/BEuXwnXXwdFHh4mnRSR7WrSA0aNh2jR4/PHYaeJR8Y/gkkvCeONXXRU7iUhxGjw4XP177rnh32IxUvHPsnnzwone008PvXxEJPsaNgy/vhcsgNtui50mDhX/LBs5MlxxeNFFsZOIFLe+faFXr9AE9NlnsdNkn4p/Fk2bBlOm/HDCSUTiWXPh14oVcEURzjGo4p8lVVWhfbF9+zBmv4jEt9tucNJJYUDFxYtjp8kuFf8sue8+qKiAyy+Hxo1jpxGRNS69NPwKGD06dpLsUvHPgm+/DU09u+4aRu4UkdzRvn3ogDFxIsyfHztN9qj4Z8HNN8OiRaF3QT3tcZGcc8EFYd7fiy+OnSR7VIoS9sUX4WTSIYfAwQfHTiMiNdl663DF74MPhubZYpBY8Tez0Wb2gZnNSi2HJbWtXDZ2bOhNcOWVsZOIyIYMHw5bbhmaaItB0kf+Y929W2p5MuFt5Zzly0NXsiOPDHPzikjuat48FP6pU+G111rGjpM4Nfsk6LrrwqxBl10WO4mI1Mbpp0O7djB+/E8LftA384T+C81sNHAy8AVQAYxw909reN5gYDBASUlJaVlZWSJ51lVZWUmzZs0Se/8VKzbjN7/Zg/32+4RRo9LvQpB03kxT3mQpb3KmTGnDmDE7ctllc9h33+Wx49RKr169Zrh79zq9yN03eQGeBebWsPQDSoD6hF8XVwB3bOz9SktLPVvKy8sTff8zz3SvX9/97bcz835J58005U2W8ibnu+/c27X70rt2df/++9hpageo8DrW77Safdz9IHfvWsPyqLt/5O6r3b0KGA/0SGdb+eT998NgUYMGwc9+FjuNiNRFgwYwaNB7zJ0L994bO01ykuzt06ba3SMJvwiKwpo2/mLqMyxSSHr2/C/dusEf/wirVsVOk4wkT/hea2ZzzOx1oBdQFFOUL1wId94JQ4ZAhw6x04jIpqhXL3TPfu89uP322GmS0SCpN3b3gUm9dy4bPRo226x4+gqLFKo+fWC//cIv+UGDCm9MLnX1zKC5c2HyZDjzTPjJT2KnEZF0mIVB35Ytg/HjY6fJPBX/DLrkEthiCzjvvNhJRCQTevaEAw6Aq6+Gb76JnSazVPwzZMYMeOihHy4RF5HCcMklsHRp4R39q/hnyCWXQOvWcHZRnNYWKR49e8L++xfe0b+KfwbMnAlPPBGO+ps3j51GRDLJLBzcffhhYfX8UfHPgMsug5YtYejQ2ElEJAm9eoWeP1ddVThH/yr+aXr9dXjkERg2DFq0iJ1GRJJQ/eh/woTYaTJDxT9Nl18eevgMGxY7iYgkqXdv2HffcPT/7bex06RPxT8N8+bBAw+Efv2tWsVOIyJJWnP0/8EHhXH0r+KfhiuuCPN+qoePSHE48EDYZ5/COPpX8d9ECxZAWRmccQZstVXsNCKSDWuO/pcsgTvuiJ0mPSr+m+jKK6FRozDps4gUj4MOgr33DjUgn4/+Vfw3wbvvwj33wGmnwTbbxE4jItlUKEf/Kv6b4KqrwoQP554bO4mIxHDwwbDnnuGq3+++i51m06j419HixXDXXfC730GbNht9uogUILMwbPv77+fvbF8q/nV09dVhooeRI2MnEZGYDj8cunYNNaGqKnaaulPxr4M1bXyDBkG7drHTiEhM9erB+efD/Pnw2GOx09Sdin8djBkTvuHPPz92EhHJBcceC9ttF84DusdOUzcq/rW0YkUYz/v446FTp9hpRCQXrOn48corMG1a7DR1o+JfSzffDF9+qVm6RGRtgwZBSUk4+s8nKv618NVXcOONP5zgERFZY/PNwxAvU6eGGf3yhYp/LdxxB3zyiXr4iEjNfv/7MKR7Ph39p1X8zexoM3vDzKrMrPs6j11gZgvN7C0zOzS9mPF8/3040bv33mE41w2aNCmcEKhXL9xOmpSFhCISW/PmYZyvhx6CN9+MnaZ20j3ynwv8CphefaWZdQGOA34B9AFuMbP6aW4rivvvh0WLwlG/2QaeOGkSDB4crgJzD7eDB+sLQKRIDBsWxvu69trYSWonreLv7vPd/a0aHuoHlLn7t+7+HrAQ6JHOtmJwD/8ju3SBvn038uRRo8LJgeq++iqsF5GCt802cOqp8Pe/w3/+EzvNxplnoHOqmU0DznH3itT9m4B/u/s9qfsTgKfc/YEaXjsYGAxQUlJSWlZWlnae2qisrKRZs2YbfM4rr7Rm5MhdGDnyTfr0WbbB5x7QuzdWw750M57/17/Sygq1y5tLlDdZypusTc27bFkjTjhhD/r3/5ChQxcmkKxmvXr1muHu3Tf+zGrcfYML8CyheWfdpV+150wDule7fxNwQrX7E4CjNrat0tJSz5by8vKNPueAA9zbtXP/9ttavGHHju7hx8LaS8eO6QVNqU3eXKK8yVLeZKWT98QT3Zs0cf/vfzOXZ2OACt9IfV132Wizj7sf5O5da1ge3cDLPgDaV7vfLrUub7z8Mjz/fOjCtdlmtXjBmmm9qmvSJKwXkaIxcmRo8f3rX2Mn2bCkuno+BhxnZo3MbDugM/BKQttKxDXXQMuWYfTOWhkwAMaNg44dw5nhjh3D/QEDkowpIjmmSxc44ohwYei6pwFzSbpdPY80syXAXsATZvYMgLu/AdwPzAOeBs5w99Xphs2Wt96CRx4JXbe22KIOLxwwIHQNqqoKtyr8IkVpxAhYvhwmToydZP3S7e3zsLu3c/dG7l7i7odWe+wKd9/e3Xd096fSj5o9110XumyddVbsJCKSj/bfH0pLYezY3B3uWVf4ruPDD0NXrUGDNEWjiGwas3D0v2ABPPFE7DQ1U/Ffx003hWnZNDG7iKTjqKOgffswQkAuUvGv5ssv4bbboH9/2H772GlEJJ81bBiu+n3++dwc8E3Fv5qJE8O4/cOHx04iIoXg1FNDp5FcPPpX8U+pqoIbboDu3WGffWKnEZFC0KJF6C5+//1hsvdcouKf8uST4eTM8OEbGcBNRKQOhg0LtzfeGDfHulT8U8aODZOyH3VU7CQiUkg6dICjjw7TwH7xRew0P1DxB2bPhn/9C4YODSdpREQyacSIUPhvvz12kh+o+BOO+ps0CcPvi4hkWvfu4cKvv/wlTBCVC4q++C9bBvfeGy7qatUqdhoRKVQjRoSTvg/8aGD7OIq++N9yS7ioa81JGRGRJPTtCzvsELp9ZmAalbQVdfH/+mu49dYwAl/nzrHTiEghq1cvDBFfUQEvvBA7TZEX/3vugU8+Cf9DRESSduKJ0Lp1bnT7LNri7x5O9O62GxxwQOw0IlIMmjQJF309/HD8i76Ktvi/+mpr5s8PR/26qEtEsuX008PtrbfGzVG0xf8f/2hHmzZw7LGxk4hIMenQIQweOW5cOO8YS1EW/3nzoKKiNUOH1nJ+XhGRDDrrrDCI5OTJ8TIUZfG/+WZo2LCq9vPziohk0P77w847h0neY3X7LLri//nncPfd0Lv3x2y9dew0IlKMzMLR/+zZ8bp9Fl3xnzgxTNpy5JEfxI4iIkXsN7+J2+2zqIp/VVWYpnHPPWHHHVfGjiMiRaxJkzDZyyOPxOn2WVTF/9lnw5j9Q4fGTiIiErp9usfp9plW8Tezo83sDTOrMrPu1dZ3MrOvzWxWavlb+lHTd9NNsM02GrNfRHJDx47Qr1+cbp/pHvnPBX4FTK/hsXfcvVtqOS3N7aTtvfdgypQwbHOjRrHTiIgEa7p93ntvdrebVvF39/nu/lamwiTp1lvDwEpDhsROIiLygwMOCN0+b7wxu90+GyT43tuZ2WvAF8BF7l5jhyYzGwysmUal0swS/TJp3/7//twK+CTJbWWY8iZLeZOlvLVQb9MPx3es6ws2WvzN7FngJzU8NMrdH13Py5YCHdx9uZmVAo+Y2S/c/UczWLr7OGBcXUJngplVuHv3jT8zNyhvspQ3WcqbLDOrqOtrNlr83f2gur6pu38LfJv6e4aZvQPsANQ5oIiIZF4iXT3NbGszq5/6+6dAZ+DdJLYlIiJ1l25XzyPNbAmwF/CEmT2Temh/4HUzmwU8AJzm7ivSSpp5WW9qSpPyJkt5k6W8yapzXvNcmExSRESyqqiu8BURkUDFX0SkCBVt8Tez0Wb2QbUhKA6LnakmZtbHzN4ys4Vmdn7sPBtjZovMbE5qn+Zk7y4zu8PMPjazudXWtTazqWb2duq2VcyM1a0nb05+fs2svZmVm9m81NAvw1Lrc3L/biBvTu5fADPb3MxeMbPZqcx/Sq3fzsxeTtWK+8xsg1NVFW2bv5mNBird/frYWdYn1WNqAXAwsAR4FTje3edFDbYBZrYI6O7uOXtBj5ntD1QCE929a2rdtcAKd7869SXbyt1Hxsy5xnryjiYHP79m1gZo4+4zzWwLYAbQHziZHNy/G8h7DDm4fwHMzICm7l5pZg2BF4FhwHDgIXcvS42nNtvd1ztkXNEe+eeJHsBCd3/X3VcBZUC/yJnynrtPB9btfdYPuDv1992EApAT1pM3J7n7Unefmfp7JTAfaEuO7t8N5M1ZHlSm7jZMLQ70JvSuhFrs42Iv/kPN7PXUz+qc+Bm6jrbAf6rdX0KOfzAJH8J/mtmM1NAd+aLE3Zem/l4GlMQMU0s5/fk1s07AbsDL5MH+XScv5PD+NbP6qa70HwNTgXeAz9z9+9RTNlorCrr4m9mzZja3hqUfcCuwPdCNMBzFmJhZC8i+7v4/wC+BM1JNFnnFQ1torreH5vTn18yaAQ8Cf1h3WJdc3L815M3p/evuq929G9CO0EKwU13fI8mB3aKr7dAUZjYemJJwnE3xAdC+2v12qXU5y90/SN1+bGYPEz6YNQ35nWs+MrM27r401Q78cexAG+LuH635O9c+v6l26AeBSe7+UGp1zu7fmvLm8v6tzt0/M7NywoW2Lc2sQerof6O1oqCP/Dck9QFc40jC3AS55lWgc+os/mbAccBjkTOtl5k1TZ00w8yaAoeQm/u1Jo8BJ6X+PglY36CFOSFXP7+pk5ETgPnu/udqD+Xk/l1f3lzdv/B/w+e0TP3dmNAhZD5QDqyZqmqj+7iYe/v8nfCTzoFFwJBqbZI5I9XF7AagPnCHu18RN9H6WRjH6eHU3QbA5FzMa2b3Aj0Jw/Z+BFwCPALcD3QAFgPH5MqQJOvJ25Mc/Pya2b7AC8AcoCq1+kJCO3rO7d8N5D2eHNy/AGa2C+GEbn3CAfz97n5p6t9fGdAaeA04ITXIZs3vU6zFX0SkmBVts4+ISDFT8RcRKUIq/iIiRUjFX0SkCKn4i4gUIRV/EZEipOIvIlKE/j9yrs6T0FBGbwAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(df.X,df.Y,'ro')\n",
"x=np.linspace(-8,30)\n",
"plt.plot(x,P( x),'b-')\n",
"plt.grid()\n",
"plt.xlim(-8,30)\n",
"plt.ylim(-15,20)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "tlOZC1tQEZXM"
},
"source": [
"__Activity__: Scipy interpolation.\n",
"\n",
"Define an interplation function which passes throgh the three points by using `interp1D` of Scipy with several linear and quadratic curves between the points."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "QCe6jFV_jz2M"
},
"source": [
""
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "uFPhevxyEZXY"
},
"outputs": [],
"source": [
"#np.poly1d?"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "oKqUP0SaEZXZ",
"tags": []
},
"source": [
"## Lagrange Polynomial"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "P09CyJ2SEZXa"
},
"source": [
"Algebraic polynomials are very special functions as they have properties like differentiability (unlike linear interpolation) and continuity that make them useful for approximations like interpolation. A Polynomial is defined as a function given by the general expression:\n",
"\n",
"$$P_n(x) = a_nx^n + a_{n-1}x^{n-1} + \\cdots + a_1 x + a_0$$\n",
"\n",
"where $n$ is the polynomial degree.\n",
"\n",
"Another important property of polynomials is given by the [Weierstrass Approximation Theorem](http://en.wikipedia.org/wiki/Stone%E2%80%93Weierstrass_theorem), which states given a continuous function $f$ defined on a interval $[a,b]$, for all $\\epsilon >0$, there exits a polynomial $P(x)$ such that\n",
"\n",
"$$|f(x) - P(x)|<\\epsilon\\ \\ \\ \\ \\ \\mbox{for all }\\ x\\ \\mbox{ in }\\ [a,b].$$\n",
"\n",
"This theorem guarantees the existence of such a polynomial, however it is necessary to propose a scheme to build it."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "jN_ZqgxNEZXs",
"tags": []
},
"source": [
"### Derivation"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "d2qULdqzEZXt"
},
"source": [
"Let's suppose a well-behaved yet unknown function $f$ and two points $(x_0,y_0)$ and $(x_1,y_1)$ for which $f(x_0) = y_0$ and $f(x_1) = y_1$. With this information we can build a first-degree polynomial that passes through both points by using the last equation in sec. [Linear Interpolation](interpolation.ipynb#Linear-Interpolation), we have\n",
"\n",
"$$P_1(x) = \\left[ \\frac{y_{1}-y_0}{x_{1}-x_0} \\right]x + \\left[ y_0 - \\frac{y_{1}-y_0}{x_{1}-x_0}x_0 \\right]$$\n",
"\n",
"We can readily rewrite this expression like:\n",
"\\begin{align}\n",
"P_1(x) =& \\frac{y_{1}}{x_{1}-x_0} x- \\frac{y_0}{x_{1}-x_0} x + y_0 -\\frac{y_{1}}{x_{1}-x_0}x_0 + \\frac{y_0}{x_{1}-x_0}x_0 \\nonumber\\\\\n",
"=& \\left[1 - \\frac{x}{x_{1}-x_0} + \\frac{x_0}{x_{1}-x_0}\\right]y_0+\n",
"\\left[\\frac{x}{x_{1}-x_0} -\\frac{x_0}{x_{1}-x_0}\\right]y_1 \\nonumber\\\\\n",
" =& \\left[\\frac{x_1-x_0}{x_{1}-x_0} - \\frac{x}{x_{1}-x_0} + \\frac{x_0}{x_{1}-x_0}\\right]y_0+\n",
"\\left[\\frac{x}{x_{1}-x_0} -\\frac{x_0}{x_{1}-x_0}\\right]y_1 \\nonumber\\\\\n",
"=& \\left[\\frac{x_1-x}{x_{1}-x_0}\\right]y_0+\n",
"\\left[\\frac{x-x_0}{x_{1}-x_0}\\right]y_1 \\,.\n",
"\\end{align}\n",
"In this way\n",
"$$P_1(x) = L_0(x)f(x_0) + L_1(x)f(x_1)$$\n",
"\n",
"where we define the functions $L_0(x)$ and $L_1(x)$ as:\n",
"\n",
"$$L_0(x) = \\frac{x-x_1}{x_0-x_1} \\mbox{ and } L_1(x) = \\frac{x-x_0}{x_1-x_0}$$\n",
"\n",
"Note that\n",
"\n",
"$$L_0(x_0) = 1,\\ \\ \\ L_0(x_1) = 0,\\ \\ \\ L_1(x_0) = 0,\\ \\ \\ L_1(x_1) = 1$$\n",
"\n",
"implying:\n",
"\n",
"$$P_1(x_0) = f(x_0) = y_0$$\n",
"\n",
"$$P_1(x_1) = f(x_1) = y_1$$\n",
"\n",
"Although all this procedure may seem unnecessary for a polynomial of degree 1, a generalization to polynomials of larger degrees is now possible."
]
},
{
"cell_type": "markdown",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"### General case"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "t5benGvIEZXt",
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"Let's assume again a well-behaved and unknown function $f$ sampled by using a set of $n+1$ data $(x_m,y_m)$ ($0\\leq m \\leq n$).\n",
"We call the set of $[x_0,x_1,\\ldots,x_n]$ as the _node_ points of the _interpolation polynomial in the Lagrange form_, $P_n(x)$, where:\n",
"$$f(x)\\approx P_n(x)\\,,$$\n",
"\n",
"$$P_n(x) = \\sum_{i=0}^n f(x_i)L_{n,i}(x) = \\sum_{i=0}^n y_iL_{n,i}(x)$$\n",
"\n",
"Such that\n",
"$$f(x_i)= P_n(x_i)\\,,$$\n",
"\n",
"We need to find the _Lagrange polynomials_, $L_{n,i}(x)$, such that \n",
"$$L_{n,i}(x_i) = 1\\,,\\qquad\\text{and}\\,,\\qquad L_{n,i}(x_j) = 0\\quad\\text{for $i\\neq j$}$$ \n",
"A function that satisfies this criterion is\n",
"\n",
"$$L_{n,i}(x) = \\prod_{\\begin{smallmatrix}m=0\\\\ m\\neq i\\end{smallmatrix}}^n \\frac{x-x_m}{x_i-x_m} =\\frac{(x-x_0)}{(x_i-x_0)}\\frac{(x-x_1)}{(x_i-x_1)}\\cdots \\frac{(x-x_{i-1})}{(x_i-x_{i-1})}\\underbrace{\\frac{}{}}_{m\\ne i}\n",
"\\frac{(x-x_{i+1})}{(x_i-x_{i+1})} \\cdots \\frac{(x-x_{n-1})}{(x_i-x_{n-1})}\\frac{(x-x_n)}{(x_i-x_n)} $$\n",
"Please note that in the expansion the term $(x-x_i)$ does not appears in both the numerator and the denominator as stablished in the productory condition $m\\neq i$.\n",
"\n",
"Moreower\n",
"$$L_{n,i}(x_i) = \\prod_{\\begin{smallmatrix}m=0\\\\ m\\neq i\\end{smallmatrix}}^n \\frac{x_i-x_m}{x_i-x_m} =1$$\n",
"and, for $j\\ne i$\n",
"$$L_{n,i}(x_j) = \\prod_{\\begin{smallmatrix}m=0\\\\ m\\neq i\\end{smallmatrix}}^n \\frac{x_j-x_m}{x_i-x_m} =\\frac{(x_j-x_0)}{(x_i-x_0)}\\cdots \\frac{(\\boldsymbol{x_j}-\\boldsymbol{x_j})}{(x_i-x_j)}\\cdots\\frac{(x_j-x_n)}{(x_i-x_n)}=0.$$\n",
"\n",
"\n",
"Then, the polynomial of $n$th-degree $P_n(x)$ will satisfy the definitory property for a interpolating polynomial, i.e. $P_n(x_i) = y_i$ for any $i$ and it is called the _interpolation Polynomial in the Lagrange form_.\n",
"\n",
"Check [this implementation in sympy](./LagrangePoly.ipynb) [[View in Colaboratory](https://colab.research.google.com/github/restrepo/ComputationalMethods/blob/master/material/LagrangePoly.ipynb)] where both the interpolating polynomial and the Lagrange polynomials are defined."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Further details at:**\n",
"[Wikipedia](https://en.wikipedia.org/wiki/Lagrange_polynomial)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "HAld-e8YEZXu",
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"### Example:\n",
"Obtain the Lagrange Polynomials for a Interpolation polynomial of degree 1."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "9WWq5hcOEZXu"
},
"source": [
"$i=0$, $n=1$\n",
"$$ L_{1,0}=\\prod_{\\begin{smallmatrix}m=0\\\\ m\\neq 0\\end{smallmatrix}}^1 \\frac{x-x_m}{x_i-x_m}=\\prod_{\\begin{smallmatrix}m=1\\end{smallmatrix}}^1 \\frac{x-x_m}{x_0-x_m}=\\frac{x-x_1}{x_0-x_1}$$\n",
"$i=1$, $n=1$\n",
"$$ L_{1,1}=\\prod_{\\begin{smallmatrix}m=0\\\\ m\\neq 1\\end{smallmatrix}}^1 \\frac{x-x_m}{x_i-x_m}=\\prod_{\\begin{smallmatrix}m=0\\end{smallmatrix}}^0 \\frac{x-x_m}{x_1-x_m}=\\frac{x-x_0}{x_1-x_0}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"Until now we can only guarantee that $P_n(x_i)=f(x_i)$. To calculate the function from any $x$ in the interpolation interval we have the following Theorem (See [here](https://www3.nd.edu/~zxu2/acms40390F12/Lec-3.1.pdf))\n",
"### Theorem\n",
"Suppose $X_0,\\ldots,x_n$ are distinct numbers in the interval $[a,b]$ and $f\\in[a,b]$.Then for each $x$ in $[a,b]$, a number $\\xi(x)$ between $x_0,\\ldots,x_n$, and hence in $[a,b]$, exists with\n",
"$$\n",
"f(x)=P_n(x)+E_n(x)\\,,\\qquad \\text{such that } E_n(x_i)=0\\,,\n",
"$$\n",
"where the formula for the error bound is given by\n",
"$$\n",
"E_n(x) = {f^{n+1}(\\xi(x)) \\over (n+1)!} \\cdot \\prod_{i=0}^{n}\\left(x-x_{i}\\right)\\,,\n",
"$$\n",
"$f^{(n+1)}$ is the $n+1$ derivative of $f$\n",
"\n",
"For a demostration see [1d] → https://www.math.ust.hk/~mamu/courses/231/Slides/CH03_1B.pdf\n",
"\n",
"The specific calculation of the bounded error is to find the $\\xi$ and $x$ such that\n",
"$$\n",
"\\left|f(x)-P(x)\\right| \\leq \\max_{\\xi\\in[a,b]}\\left|\\frac{f^{(n+1)}(\\xi)}{(n+1) !}\\right| \\cdot \\max_{x\\in[a,b]}\\left|\\prod_{i=0}^{n}\\left(x-x_{i}\\right)\\right|\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "me-qxh8sEZXv"
},
"source": [
"### Exercise-interpolation\n",
"Obtain the Lagrange Polynomials for a Interpolation polynomial of degree 2."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "QrN1hmlnEZXv"
},
"source": [
"### Implementation in Scipy"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "IRQNq9xrEZXw"
},
"outputs": [],
"source": [
"from scipy import interpolate"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "SP7TER3gEZXy"
},
"outputs": [],
"source": [
"#interpolate.lagrange?"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 141
},
"colab_type": "code",
"id": "-p_85Kiz5Xi-",
"outputId": "f6fe740e-5ba0-4407-d20f-50672c7926ef"
},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" X Y\n",
"0 3.0 8.0\n",
"2 10.0 6.5\n",
"1 21.3 3.0"
]
},
"execution_count": 102,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df"
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.0"
]
},
"execution_count": 103,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df.X[0]"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "FEj64uDmEZX8"
},
"source": [
"### Steps LP"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "t6ks3grEEZX8"
},
"source": [
"Once defined the formal procedure for constructing a Lagrange Polynomial, we proceed to describe the explicit algorithm:\n",
"\n",
"1. Give the working dataset $(x_i, y_i)$ and stablish how many points you have.\n",
"2. Define the functions $L_{n,i}(x)$ in a general way.\n",
"3. Add each of those terms as shown in last expression.\n",
"4. Evaluate your result wherever you want."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "FbuCBnS4EZX9"
},
"source": [
"### Activity ###\n",
"\n",
"Along with the professor, write an implementation of the previous algorithm during classtime."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "-1r8xBk8EZX9"
},
"source": [
"### Activity LP"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ZkuEZPFTEZX-"
},
"source": [
"
\n",
" \n",
"
\n",
"\n",
"One of the very first evidences of the existence of dark matter was the flat rotation curves of spiral galaxies. If we assume the total budget of mass of a galaxy is entirely made of luminous matter, the orbital circular velocity of stars around the galaxy plane should decay according to a keplerian potential. However this is not the case and the circular velocity barely decreases at larger radius, thus indicating the presence of a new non-visible matter component (dark matter). When it is necessary to determine how massive is the dark matter halo embedding a galaxy, an integration of the circular velocity is required. Nevertheless, due to the finite array of a CCD camera, only a discrete set of velocities can be measured and interpolation techniques are required.\n",
"\n",
"\n",
"In this activity we will take a discrete dataset of the circular velocity as a function of the radius for the galaxy [NGC 7331](http://es.wikipedia.org/wiki/NGC_7331) and perform both, a linear and a Lagrange interpolation. You can download the dataset from this [link](https://raw.githubusercontent.com/sbustamante/ComputationalMethods/master/data/NGC7331.dat).\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "6BmEFS28EZX_"
},
"source": [
"[Video](https://upload.wikimedia.org/wikipedia/commons/transcoded/3/33/Galaxy_rotation_under_the_influence_of_dark_matter.ogv/Galaxy_rotation_under_the_influence_of_dark_matter.ogv.360p.webm)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "fseOFy-DEZX_"
},
"source": [
"\n",
"**TRIVIA** \n",
"To which of two curves the real data approach better?\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "XaM0WuPQEZX_"
},
"source": [
"import os\n",
"os.remove('trivia_results.txt')"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "pcbVK4pkEZX_",
"outputId": "c24c3075-c96b-4809-b64a-f3ce6b049926"
},
"outputs": [
{
"name": "stdin",
"output_type": "stream",
"text": [
"A: to the curve \"velocity goes to zero when distance goes to infinity\"\n",
"B: to the curve \"velocity goes to high constant when distance goes to infinity\"\n",
" A\n"
]
}
],
"source": [
"f=open('trivia_results.txt','a')\n",
"AB=input(r'''A: to the curve \"velocity goes to zero when distance goes to infinity\"\n",
"B: to the curve \"velocity goes to high constant when distance goes to infinity\"\n",
"''')\n",
"f.write( '{}\\n'.format(AB) )\n",
"f.close()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "xQ1a6G2ZEZYB",
"outputId": "be642608-3aa3-4645-c2d7-a93e5e8bd029"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"C\n",
"A\n",
"\n"
]
}
],
"source": [
"fr=open('trivia_results.txt')\n",
"print( fr.read())\n",
"fr.close()"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "WsMD7hIjEZYD"
},
"source": [
"os.remove('trivia_results.txt')"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "IAYkvn-_EZYD",
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"### Lets us check! "
]
},
{
"cell_type": "code",
"execution_count": 133,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "JiCl9H_TEZYE",
"outputId": "329bce3d-5d4f-40fd-f726-dc0755207e29"
},
"outputs": [
{
"data": {
"text/html": [
"