lab02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ASTR 2100, Winter 2020\n",
    "\n",
    "## Lab 2: approximating functions with polynomials and piecewise splines\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Jupyter notebooks with class notes and assignments: https://github.com/a-kravtsov/a211w20\n",
    "\n",
    "#### Slack workspace for this class: https://a211w20.slack.com\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Instructor: Andrey Kravtsov (kravtsov@uchicago.edu)\n",
    "\n",
    "#### TAs: Dimitrios Tanoglidis (dtanoglidis@uchicaago.edu), Georgios Zakharegkas (gzakharegkas@uchicago.edu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import packages needed by the codes below. Run this cell first before using these codes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# use jupyter \"magic\" command to tell it to embed plot into the notebook \n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from codes.plot_utils import plot_pretty\n",
    "plot_pretty(fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: Using splines for interpolation and estimates of  derivatives and integrals\n",
    "\n",
    "\n",
    "**Background.** I introduced splines as a way to construct piecewise polynomial approximation that is continuous and smooth and also has smooth first derivative (see the <a href=\"https://github.com/a-kravtsov/a211w20/blob/master/03_approximation.ipynb\"><tt>03_approximation.ipynb</tt></a> notebook). Splines not only provide a general flexible and accurate way to approximate a function represented by some tabulated values at the training points, but can also automatically provides approximations for the function first three derivatives and its integral. Moreover, splines can be used as a tool for numerical solution (finding roots) of equations and inverting functions. In this sense splines are one of the \"Swiss army knife\"-type numerical methods that are very useful to have in one's toolbox. \n",
    "\n",
    "In this exercise you will practice using splines for these tasks. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What you are learning in this exercise\n",
    "\n",
    "* To use polynomials and splines to approximate a function using SciPy routines\n",
    "* Use splines to obtain accurate and smooth estimate of the function derivative\n",
    "* Use spline to obtain accurate estimate of the function integral for integration intervals that fall within the interval, in which the approximating spline was constructed.  \n",
    "* Use splines to invert a complicated function (which is equivalent to finding a root of a function). Inversion is when we instead of estimating $y=f(x)$, find $x$ that satisfies such equation given a $y$ value. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data used in this exercise\n",
    "\n",
    "For this exercise we will use a function $P(k)$ that is tabulated in a data file. \n",
    "\n",
    "For the purposes of this exercise it is not needed to know what this function is. This is something called \"power spectrum\" - the concept we will come back to when we will discuss Fourier transforms (also see below task 2). This particular power spectrum is often used in astrophysical calculations related to galaxies and the spatial structures in galaxy distribution. The tabulated values in the file were computed using a sophisticated code <a href=\"https://camb.info/\">CAMB</a> that solves a complicated system of differential equations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in function that will be used for this exercise\n",
    "k, Pk = np.loadtxt('data/test_matterpower_logintk1000.dat', usecols=(0, 1),  unpack=True)\n",
    "\n",
    "# their log10 values that will be used \n",
    "lk, lPk = np.log10(k), np.log10(Pk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAH4CAYAAAC/lWwoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAa6QAAGukBho3xbQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3db4xc13nf8d/DWe5SokguKcv/Zy0tlbp2ndompdrtrtHY5DpNgKJoQkpv/MJFbdIpGiRFXDJ+012/iU27LYq0QE06MNrCBSqTSoICbtBwaXcdrIFWIh0ZbpPY5kqhVEeVbHJJiRR3ydmnL+bOcmY4O5yZM3fOPbPfD7CYnTv33DlarGZ/POe555i7CwAAIFVbYncAAAAgBGEGAAAkjTADAACSRpgBAABJI8wAAICkEWYAAEDSCDMAACBphBkAAJA0wgwAAEgaYQYAACSNMAMAAJJGmAEAAEkjzAAAgKQRZgAAQNIIMwAAIGkjsTuAdJiZZ9+uRu0IAGBYjEqSu1vIRQgz6Nro6Oho7D4AANK3utqffxsTZtCN1dHR0dGVlZXY/QAADIGxsTGt9iHRUDMDAACSRpgBAABJI8wAAICkEWYAAEDSCDMAACBphBkAAJA0wgwAAEgaYQYAACSNRfPQtcXFxYbn5XJZExMTkXoDANjszN3vfRYgycxWSqXSaKVSaTg+Ozurubm5OJ0CACSrtgKwu4+FXIeRGXSlVCppYWGh4Vi5XI7UGwAAGJlBF8xshb2ZAAD90q+RGQqAAQBA0ggzAAAgaYQZDNzylev6+tcW9fzFV2N3BQAwBKiZQcf6VTNz/Dee0ve/96Ik6Y++89l+dA0AkCBqZpCsWpABAKAfCDMAACBphBkAAJA0wgwAAEgaYQYAACSNMAMAAJJGmAEAAEkjzAAAgKQRZgAAQNIIMwAAIGmEGQAAkDTCDAAASBphBgAAJI0wAwAAkkaYAQAASSPMAACApBFmAABA0ggziMrdY3cBAJC4kdgdQHoWFxcbnpfLZU1MTPR0LXfJrB+9AgBsVoQZdKVSqWh6errh2OzsrObm5nq6nq+5tIU0AwDoHWEGXSmVSlpYWGg4Vi6Xe76ei2kmAEAYwgy6NjU11bdrUTIDAAhFATCiogAYABCKMIOoCDMAgFCEGURFlgEAhCLMIC7SDAAgEGEGUa2RZQAAgQgziIqaGQBAKMIM4iLLAAACEWYQ1RojMwCAQIQZxEWYAQAEIswgKrIMACAUYQZRMc0EAAhFmEFchBkAQCDCDKIiywAAQhFmEBVhBgAQijCDqFg0DwAQijCDqAgzAIBQhBlERZYBAIQizCAqRmYAAKEIM4iKMAMACEWYQVxkGQBAIMIMomIFYABAKMIMoiLLAABCjcTuANKzuLjY8LxcLmtiYqK3i5FmAACBCDPoSqVS0fT0dMOx2dlZzc3N9XQ9ppkAAKEIM+hKqVTSwsJCw7Fyudz7BckyAIBAhBl0bWpqqm/X4tZsAEAoCoARFdNMAIBQhBnERZYBAAQizCAqppkAAKEIM4iKLAMACEWYQVSMzAAAQhFmEBVZBgAQijCDqBiZAQCEIswgKsIMACAUYQZRkWUAAKEIM4iKkRkAQCjCDKIiywAAQhFmEBUjMwCAUIQZxEWWAQAEIswgKjaaBACEIswgKqaZAAChCDOIiywDAAhEmEFUjMwAAEIRZhAVYQYAEIowg6jIMgCAUIQZREWYAQCEIswgKqaZAAChCDOIijADAAg1ErsDSM/i4mLD83K5rImJiZ6uRZYBAIQizKArlUpF09PTDcdmZ2c1NzfX0/UYmQEAhCLMoCulUkkLCwsNx8rlcs/XI8wAAEIRZtC1qampvl2LMAMACEUBMKIiywAAQhFmEBVhBgAQijCDuEgzAIBAhBlEtUaYAQAEIswgKrIMACAUYQZxkWYAAIEIM4iKaSYAQCjCDOIiywAAAhFmEBWL5gEAQhFmENUaWQYAEIgwg7gYmQEABCLMYODM7ny/xtAMACAQYQYDZ/VpBgCAQJtq12wzOyLpqKTx7NAFSU+5+5l7tNsn6XOSLmeHJiWdHJZ2MVEADAAItWnCjJmdlLTs7vuz5+OSTkg6bWbzkg67+3KLdkey8/a7+1Jd23NmNuPuRzd4vyTaRcE0EwCgjzZFmDGzY6oGmeO1Y1lwOZpNeRyR9FVJh5va7ZN0UtJMLSDU2prZAUlXzOy8u59KsV0sJhMLzAAA+mWz1MwcrQ8yTWrHD5nZwabXTqgaguabG2VhaD47p1kq7aJjmgkAEGrow4yZHZK0Yc1I9sf+QvZ0punlg5KebXP5C5LGs/dIsV0U9fW/vhavHwCA4TD0YUbVIthjZna6zTm1qZl9tQN1f/iX7j593cXscT0EpdIuqro040w3AQACbYYwU3MoK4htpXa8vgD48ezxojZWCxCPJdgumi31YYYsAwAItBkKgE+p+gf/mVZ3K2Vqf+TrRzcmO7h27Rbo+pCUSrtoGqeZSDMAgDBDH2ayAHN4o9ezO4Fqf+S/UPdSq9GajexJsF089WGGaSYAQKChDzMd+Fz2+KWmkZtu/vDXj3ik0q4nq6urGhsba3vOyspK29frVwBmZAYAEGoz1czcJRuVOSRpvs2t2+iz+s0MqJkBAITa7CMzpyVdcPdWd/lcbnGsE6m068no6Og9R17uxbbUFwCTZgAAYTbtyIyZnVW1zuTABqd0UoPS6txU2kVj3M0EAOijTRlmsn2a9kg60OYOp5p2dSa1epVWoyOptIuKkRkAQKhNN82Ubcj4WG3DyTZqt2nvbXNO/e7bqbUrBLIMACDUphqZyfZeOrpRkDGz+v2LzmaP7e4Wqq3x8kyC7QqBkRkAQKhNE2ayO5eOtwky46qbqsk2bVxW3RYHLdRW313f+ymVdkVBmAEAhNoUYSYLKifUtHiemY1nX5OqrjdzvqnpKUmT2eut1G7rbt4XKZV2UdQHGLIMACDUpggzqoaUg5KumJnXviRdyb4uSjqmps0as7VnllQNQg3M7JiqIyJHm19LpV0RMDIDAAg19GEmq4PpZP8iSXq2xbH9ksbN7LSZTWZfJ1QNBwfajHak0i4qwgwAIJTxx6QzWc3N+oaUWa3K0LTr8Noro6Ojo6GL5v3qL/2ublxflSR98shH9OQnPtSP7gEAEjM2NqbV1dVVd2+/T849bLpbs3vl7hfUw63NqbQbpPr8TJgGAIQa+mkmFBtZBgAQijCDqBiZAQCEIswgKsIMACAUYQZxkWUAAIEIMxi8utGYNUZmAACBCDOIimkmAEAowgziIssAAAIRZhAV00wAgFCEGQycb/gEAIDuEWYQFTUzAIBQhBlEtbZGmAEAhCHMAACApBFmEBUjMwCAUIQZDJw3bpsdryMAgKEwErsDSM/i4mLD83K5rImJiZ6uRZQBAIQizKArlUpF09PTDcdmZ2c1NzfX0/WYZgIAhCLMoCulUkkLCwsNx8rlcu8XJMsAAAIRZtC1qampsAvUl8yQZgAAgSgARlTONBMAIBBhBlFxMxMAIBRhBlExzQQACEWYQVS+FrsHAIDUEWYwcI1r5jEyAwAIQ5hBVIQZAEAowgyiIssAAEIRZhAVIzMAgFCEGQxc/R1MhBkAQCjCDKIiywAAQhFmEBUjMwCAUIQZREWWAQCEIswgKkZmAAChCDMYuMZF8+L1AwAwHAgziIqRGQBAKMIMoiLMAABCEWYQFWEGABCKMIPBc2/1LQAAPSHMICpGZgAAoQgziMrXYvcAAJA6wgyiqt+nCQCAXhBmEBdZBgAQaCR2B5CexcXFhuflclkTExMdt68vk1lbI80AAMIQZtCVSqWi6enphmOzs7Oam5vr6XpMMwEAQhFm0JVSqaSFhYWGY+VyufcLkmUAAIEIM+ja1NRU367FNBMAIBQFwBi4xqklwgwAIEzUkRkz2ylpT/b0srtfi9kfDB5r5gEAQg0kzGSh5QlJM5L2qRpgxlucJ0lLki5I+mNJpwk4w41pJgBAqFzDjJn9iqSjkg5KMt0JKpclLUv6WfY4LunB7HGPpP2SDks6ZWbzkr7o7t/Os68AACBNuYQZM3tY0hlVR2FOqToqM+/uV7u4xi5VR3IOSjpnZuclfcrdn+t7hxENIzMAgFB9LwA2s3+u6ujLU5J2u/tn3P3pboKMJLn7VXc/4+6fkbRb0mlJ3zaz3+p3nzFgDfW/hBkAQJi+hhkz+4qkx9x9j7t/udsAs5Es2HzJ3fdI+pCZ/ft+XBfxEWUAAKH6FmayEZnz7v5kv67Zirs/IWnJzD6b5/tgMJhmAgCE6mfNzBl3f76P19uQu3/ZzB4ZxHshZ2QZAECgvo3MDCrIxHo/9I/X1cmwNxMAINRAVwA2s51m9lkz+8Ag3xfF5UwzAQACDXo7g9+T9CVJz2YL6TUws181s6fM7GMD7hci4WYmAECoQYeZy5KOS/rTVi9mt3A/KWl3FmoeHmDfEAHTTACAUAPfaDK7ZfuxdtsU1IWaEwPsGiJgmgkAEGrQYea0mX1qoxdbTT1h+NRPLa0xzwQACDTQMOPu5yRZtmfTOjPbZWY/knTFzCpm9iMz+9kg+4Y4fC12DwAAqctjO4OH273u7l+V9GBTke9pSXslnZP0PVU3pTwt6Ui/+4diWVsjzQAAwuSx0eSSmbmks6oGknPu/kL9Ce7+VTP7tJlJ0nlJj0kab1dHg+HECsAAgFC57Jqt6sjKx1Xd9VpmtqxquJlXdffsF7JA80VVd9b+AkFmcyLMAABC5RFmXNUpo92SDqoaag5IekLSYakh3JyT9BVJ+3PoBxJAmAEAhMojzDyfTSu9oGr9y5clycw+qNbhxlS9y2le2egNozSbB7dmAwBC9T3MuPujGxz/nlqHm5nsca+ygl8zW1I12Jxy95YL7GE4VCgABgAEyqtm5p42CDczkp6U9EFVw82Tkh6M1Ue0tri42PC8XC5rYmKip2sxzQQACBUtzDSrCzdfkiQzOxC3R2ilUqloenq64djs7Kzm5uY6au9Ni+QxzQQACFWYMNMsW2APBVMqlbSwsNBwrFwu93w9RmYAAKEKG2ZQXFNTU327FmEGABCqbysAm9kH+nWtDt/v4UG+H/LBCsAAgFD93M7gcTP7bB+vt6HsfQ4N4r3QX837SjIyAwAI1bcwk+259KiZ/Xcz29Gv69Yzs51m9seS9rr7v8zjPTBYFAADAEL1daNJd/+MpKcl/aWZ/Y6Z7ezHdbMQ80VJz0v6hrv/Wj+ui/gqhBkAQKC+75rt7qdU3TjyMUlXzOwpM/uH3QYbM3vYzH7FzJ6SdEXVtWcec/ff63efEU/zrdoAAHQrl7uZ3H1J0sfN7KCqq/o+LcmzPZmelbQkaVnSz+qaPShpUtK4sl20Vd3qYF7Sx7lVezitVQgzAIAwud6a7e7zkubNbJeqq/s+oWpgmdmgybKqQee07uzTdDXPPmKwmkdiuJsJABBqIOvMZIHkTPa1Lgs5eyRdJrRsTmtMMwEAAkVdNC8LMISYTYxpJgBAqL4XAAPdYq0ZAEAIwgwGq0VuYa0ZAECIXKeZzOxjkvZJWnL338/zvZCuytqaSuRqAECPcgkzZvZBVW+pHq87dkXSp939D/J4T6SLtWYAACHy+ufwOUm7VV0npva1R9IZM/vHOb0nEkURMAAgRN/DTLbtgCQddPctqoaa/ZK+LOmapFNm9q5+vy/SRQEwACBEHiMzByV9yt2/JVVvv3b377n7cUmPSPqWpJM5vC8S4C0qgFlrBgAQIo8w88hGxb7uvuzuM6rurv2xHN4bCVqrsAowAKB3eYQZ6+Cc35Z0OIf3RoKYZgIAhMgjzFy+1wnufkbVPZoAVRiZAQAEyCPM8M9sbKhVeUzlNmEGANC7mCuVdTIdhU3gNmEGABAgjzCzp8PzlnN4byTodqUSuwsAgITlsQLwbjP7qaSvSnpG0ry7X8vhfRDJ4uJiw/NyuayJiYmer8c0EwAgRF57M+2RdKz2xMyWVN3e4KzuhBtqaxJUqVQ0PT3dcGx2dlZzc3M9X5NpJgBAiLzCzHFJe1VdQG8y+36vpCPSerjZk21tcM7dX2h1ETP7kbv/XE59RA9KpZIWFhYajpXL5Y7bt9qHiZEZAECIPMLMsrt/ufbEzHapGmpm1BhuJOlUds6yGkduXsher52HApmamurr9W7fpmYGANC7PMLMN+qfuPtVSU9nXxuFm92qLqJ3KDtnWdKzYipqU7jNOjMAgAB9DzPu/pl7vN5puJkRYWZTYJoJABAi5jozktY3onza3T/j7o/qzijNuchdQx5axFMKgAEAIaKHmWZ14WYmdl8wGBVqZgAAAQoXZppcjd0B5I+RGQBAiKKHmeOxO4D8sdEkACBEocOMu381dh+QP27NBgCEKHSYwfDxFhXATDMBAEIQZhAdt2YDAEIQZhAdIzMAgBCEGURXqVAzAwDoHWEGA9Vin0lGZgAAQQgziI6aGQBACMIMomNkBgAQIo9dswvPzE5LOunu8x2cu0/S5yRdzg5NZm3PDEO7ImCdGQBAiE0TZsxsUtVduY+r+of+qQ7aHJF0QtJ+d1/Kjo1LOmdmM+5+NOV2RcEKwACAEEMfZszsoKSTkpYknc4eJztoty9rN1MLCJLk7stmdkDSFTM77+6nUmwXTYsK4Nu3GJkBAPRu6Gtm3H3e3fe6+0z2B325w6YnJC23mopy92VJ89k5qbYrDGpmAAAhhj7MBDgo6dk2r1+QNG5mhxJtVxi3mWYCAAQgzLRQ94d/qc1pF7PHmdTaFQ23ZgMAQhBmWns8e7zY5pxagHgswXbRtF40j5oZAEDvCDOt3bNAWHdugR5PsF2hMDIDAAhBmGmt9oe/k2LhPQm2KxQKgAEAIYb+1uwedfOHv37EI5V2PVtdXdXY2Fjbc1ZWVrq6JhtNAgBCMDKD6BiZAQCEYGSmtcv3PiXpdj0bHR3teuSlXov6X8IMACAIIzOtdbqwXvO5qbQrFAqAAQAhCDPttaszqdWrtBodSaVdIdymZgYAEIAw01ptbZa9bc6pBYgLCbYrFEZmAAAhCDOtnc0e290tVFvj5ZkE20XjdavmjYxUf/2omQEAhCDMtJBt2rgsaV+b02qr755JrV1RlLIww8gMACAEYWZjpyRNmtlGq+wekjTv7s37IqXSLrqRkZIktjMAAITZjGGm9kf/8XYnuftxVWtSTjS/ZmbHVB0ROZpquyIolbJpJnbNBgAEGPowY2aTZnbezC6a2RXdmZI5ZmZXstfObtB8v6RxMzudXWfSzE6oGg4OtBntSKVdVCNMMwEA+sC81TbGaGBm+3Rn1+mlrFZlaNp1ysxWRkdHR0MWzbv++ooO/fK/lSS9+a079crL1/TAjm06/c1/2q9uAgASMTY2ptXV1VV3b79Pzj2wAnAH3P2Ceri1OZV2sYyUaiMz1MwAAHo39NNMKK7a3Uy3mGYCAAQgzCCa2t1MjMwAAEIQZjBQDYvmbd2SHZMq3NEEAOgRYQbRbM1GZiTuaAIA9I4wg2hqNTMSa80AAHpHmEE0Iw0jM9TNAAB6Q5jBQNWvalRbAVhis0kAQO8IM4hm69Y7IzOEGQBArwgziKa+ZoZpJgBArwgziGaEAmAAQB8QZhBNfc0Mt2YDAHpFmMFg1S2aR80MAKAfCDOIpmGdGWpmAAA9YtdsdG1xcbHheblc1sTERNfXGSmxAjAAIBxhBl2pVCqanp5uODY7O6u5ubmur8UKwACAfiDMoCulUkkLCwsNx8rlcsft60pmGmtmbjHNBADoDWEGXZuamurLdRprZhiZAQD0hgJgRDPCrdkAgD4gzCCaxpoZppkAAL0hzCCaxl2zGZkBAPSGMIOB8roK4BFqZgAAfUCYQTQjWxmZAQCEI8wgmhFqZgAAfUCYQTT1NTNMMwEAekWYwUDVrZknM9OWLSaJaSYAQO8IM4jG7M7t2Ww0CQDoFWEGUdUWzmNkBgDQK8IMoipldTNsNAkA6BVhBlHV7mhiZAYA0CvCDAarfttsUTMDAAhHmEE8Zus1M9yaDQDoFWEGUdXWmmFkBgDQK8IMoipRMwMACESYwUA1lcysFwAzzQQA6BVhBtGY6guACTMAgN4QZhDVSKlaM1Nho0kAQI8IM4iKkRkAQCjCDKLaujUbmSHMAAB6RJjBQDUXAJdYZwYAEGgkdgeQnsXFxYbn5XJZExMTXV/HrH47A2pmAAC9IcygK5VKRdPT0w3HZmdnNTc319P11mtm2GgSANAjwgy6UiqVtLCw0HCsXC73fL07KwATZgAAvSHMoGtTU1O9N24qmmHXbABAKAqAEY2Z1RUAUzMDAOgNYQZR1aaZGJkBAPSKMIOoKAAGAIQizCAqbs0GAIQizGCgmtbMW6+ZuXWLkRkAQG8IM4jH6mpmmGYCAPSIMIOoRka4mwkAEIYwg6hKrDMDAAhEmMFA+QaL5rECMACgV4QZRGMylUp3amaagw4AAJ0gzCCq2siMRBEwAKA3hBlEte2+revf33zjVsSeAABSRZhBVPdvH1v//vr1lYg9AQCkijCDgWoui3nggbow8zphBgDQPcIMojGTtm8nzAAAwhBmENX2HYQZAEAYwgyi2k7NDAAgEGEGg9VUNFMfZl67dnPQvQEADAHCDOIxU2lkix7YsU2SdO3qG5E7BABIEWEG0e3adZ8k6SphBgDQA8IMots1Xg0zjMwAAHpBmEF0O3cRZgAAvRuJ3QGkZ3FxseF5uVzWxMRER21bbSW5a/x+SdLV5RuhXQMAbEKEGXSlUqloenq64djs7Kzm5ua6vpZZ9ZGRGQBACMIMulIqlbSwsNBwrFwuB11z567a3Uw35e6yWsoBAKADhBl0bWpqqq/X2/5ANcxUKmtauXlL2+4b7ev1AQDDjQJgDJQ37zSpps0mr68OsjsAgCFAmEE0pup00v1sNgkACECYQXTbHyDMAAB6R5hBdNu336mRYbNJAEC3CDOIrlYALDEyAwDoHmEGg9Vi1bztD9SNzBBmAABdIswgmtpyMqOjIxoZqf4qMs0EAOgWYQbRmdl6ETBhBgDQLcIMCmF7dnv29dcIMwCA7hBmMFCtFs2TxMgMAKBnhBnEU7cHU23hPMIMAKBbhBkUQm1LA+5mAgB0izCDQrifMAMA6BFhBoVQKwC+wTQTAKBLhBkMVOvy3zvTTK8zMgMA6BJhBtHU1f+uFwC/cWNVa2sbRR4AAO5GmEEh1LY0cK8GGgAAOkWYQSHs2HHf+vdXl29E7AkAIDWEGQzUBmvm6S1v27n+/ct/dXVAvQEADAPCDKKpr5l5y9t2rX9PmAEAdIMwg0LYvn1M47vvlyS9dOly5N4AAFJCmEFhPDz5JknS0o9fjdwTAEBKCDMojMm9D0mSnr/46oYbUgIA0GwkdgeQnsXFxYbn5XJZExMTnTVuE1IeefTNkqRrV9/Qz376ut700I6e+wgA2DwIM+hKpVLR9PR0w7HZ2VnNzc31cDVreDb56EPr3y/9+FXCDACgI4QZdKVUKmlhYaHhWLlc7su1y+96UKXSFlUqa3ph6VX9rb892ZfrAgCGG2EGXZuamsrlulu3lvT2d4zrxUuX9dKLV3J5DwDA8KEAGAN1r7red07skSS99Jfcng0A6AxhBtGY3X3sHeXdkqQXL13mjiYAQEcIMyiU2sjM66/d1NWrb0TuDQAgBYQZFEo5CzMSKwEDADpDmEGhlN91J8xc/OErEXsCAEgFYQYD5WpfB7Nj5316JFsJ+H9+9+IgugQASBxhBtFYqwpgSR+e3itJeu57l7S8fGOQXQIAJIgwg8L5uwf+uiRpreL6k2/9ReTeAACKjjCDwnnXw2/S3p+r7tP07bN/Frk3AICiI8xgsDpcOqY2OvNn//snevX/XcuxQwCA1BFmEM0GJTOSpI/8wl9b//5P/scPB9AbAECqCDMopLe+fVzvfs9bJUl/8I3zWlm5FblHAICiIsygsJ78xIclST999TV9/WvfjdwbAEBREWaGjJntM7PTZnYy+zprZodi96sXH57eq32PPyxJevq/PKMffP+luB0CABQSYWaImNkRSeckHXf3o+5+VNJhSZ8zs5Nxe1fVzeaRZqbfPP6L2v7AmNylf/U7f6Q3bqzm2DsAQIoIM0PCzPZJOinpsLsv1Y67+7KkA5KOZGGnQNpUAGceevMO/ZPfPCBJevknV/Uvjv2+Xn/tZt4dAwAkhDAzPE5IWnb3+eYXskAzn52TnI/OvEcHfvG9kqQffP8lHf+Np/STl65E7hUAoCgIM8PjoKRn27x+QdJ4ivUzZqZ/9tt/T7/09/+mJGnpx6/q1z75H/Wf/8N3df36SuTeAQBiI8wMgbqAstTmtNqujTM5d6etLkpmGpRKW/Trn53RJ498RCMjW7S6eltf/9p39cnDp/Tv/vVZ/fn/+SutrfV4cQBA0kZidwB98Xj22G6b6VrQeSznvnSs3aJ5rc83PfmJD+lDf2dSX/ndb+u5C5f0+usr+uYfPqdv/uFz2rnrPn1g/4Te+7536JFHH9Ijex/Sjh3b8uk8AKAwCDPDYbKDcy5nj+N5dmQQHp58SF/8N0/oB8+9pP/2X5/Td7/zI62s3Na1q2/oO9/6C32nbnPK8d33681v2amH3rJDb37LTu3es107d96nHbu2rT9u3z6mbdu2amzbVm3dWor4XwYA6AVhZjjUAspyB+fuybMjg/S+979T73v/O3X9+oou/K8X9KfnL+m5C5f0f+uKg5ev3NDylRv64Z+/3NE1S6Ut2nbf1mq4GRvR2LatGtm6RSOlkkojWzQyskUjIyWNjGzJnpfWj9VeL5W2yMy0ZYvJzGQm2RbTFrOGxzvnaP3c9TZN58uqI1PSnREt0/o3jcdbDHnd1bbpnObX64fNmg+ZNb2vWretf4uNzumnHC6Z43VT+u9Ppa/5/ADy6GseP9NY3vM33qZt943G7gZhZkh0E1CCRmZWV1c1NjbW9pyVlcEW5W7fPqaPfPTd+shH3y1JunFjVS9cfFXPL72ql39yVa+8fE2vvHJNr7x8TVeX31ClsrbhtSqVNV1/fUXXX6ewGADu5eR/+keaePjB2N0gzGCwdj+4Xb/+W9Ua5PK78vkf4KvM6z0AAAiNSURBVP77R/Xen3+H3vvz77jrNXfXjRureu3qG7p27aZeu/aGblxf1c2bt7Ry85Zu3rylm2/c0srK7ezxlm7fXlPl9ppu317T7UpFt2+tqXK7otuVtTuv3ao9r2htzeVrrjWvPrp79ZhX37/+teqj7joXANA5wsxwuHzvU/pjdHQ0aORlx45t+uV/8P4+9qg7Zqbt28e0ffuY3vr2aN24pzsB6M5j9fj6CQ3Pva7dnYs0Hlt/pbntBq83ntPUtun9Nnq9oU9N79dPecW/HLqay0Vz6ackz+Mnm8clc/oB5HHZXH6mEb317btid0ESYWZYdFIr08u5iMTMVCoNz7w6AOSJdWaGS7t6mFpdzcBGcQAAGATCzHCorSGzt805taBzIee+AAAwUISZ4XA2e2x3V1NtLZpncu4LAAADRZgZAtnmksuS9rU5rbZK8Jn8ewQAwOAQZobHKUmTZrbRasCHJM27e7v9mwAASA5hZki4+3FVa2dONL9mZsdUHbk5Ouh+AQCQN8LMcNkvadzMTpvZZPZ1QtUQc6AoozJjY2P3XEUYSBm/4xh2RfsdZ52ZIeLuy5JmzGyfpIPZ4bPZqA0AAEOJMDOE3P2CuAUbALBJMM0EAACSRpgBJF26dElzc3O6dOlS7K6so0+dKWKfiqiIPyf61Jki9qloCDOApBdffFGf//zn9eKLL8buyjr61Jki9qmIivhzok+dKWKfioYwAwAAkkaYAQAASTN3j90HJMLMXJJGR0eDrrO6uqp+XKef3F23bt3S1q1bZWaxuyOJPnWqiH3id7wz9KkzRexTv37Ha9dx96D/MMIMOlYLM5JWo3YkPyVJldidaEKfOlPEPhVREX9O9KkzRexTP4xKhBkAALDJUTMDAACSRpgBAABJI8wAAICkEWYAAEDSCDMAACBphBkAAJA0wgwAAEgaYQYAACSNMAMAAJJGmAEAAEkjzGAgzGyfmZ02s5PZ11kzOxS7X0A/mNl5M5vc4LWDZnZs0H0CQmSf1wc7PDf65/vIIN8Mm5OZHZF0QtJ+d1/Kjo1LOmdmM+5+NGoHgXD7JF00s1OSzkpaljQu6UlJByXtj9g3oCNZID8o6bikSUlPddCmEJ/vhBnkysz2STopaab2iy5J7r5sZgckXTGz8+5+Klongf45kn3VLEk6UP+7DxRNNgJzUtXf19PZY8uRxqZ2hfl8Z9ds5MrMzkp6zN139/I6kAIzuyLpWVX/VStJ85LOuvuX4vUK6I2ZnZZ0SNJhdz/T5rzCfL5TM4O8HVT1Q34jFySNUz+DxD3r7jPubtnXDEEGm0BhPt8JM8hN3S9wuyH2i9njTM7dAQD0SdE+3wkzyNPj2ePFNufU/kd4LOe+AAD6p1Cf7xQAI0/3LCCTdDl7HM+zI8AgmNkJ3fm9H5d0muJ2DKlCfb4TZpCn2i/wcgfn7smzI0DOJrMg8wV3X5buuj31cNzuAX1XqM93ppmQp25+gRmZQcrm3f14LchI1dtTJX1a0qEs6ADDpFCf74QZAAi00cJg7n5B1X+5HstGagDkgDCDPF2+9ynA0JvPHp+I2gugvwr1+U6YQZ46mUvt5VwgJbUPfbY0wDAp1Oc7YQaD0G54vTbvWqiUD+Sgk7s/gNQU4vOdMIM81dYY2NvmnNr/CBdy7gvQd9mO2FfuUeDLnXoYRoX6fCfMIE9ns8d2H+a1f60+k3NfgDzMqPqBfaTNOQR2DKNCfb4TZpAbd59Xda50X5vTaqtIbriZGVBgz6ha4HugzTm11U+fyr87wGAU7fOdMIO8nVJ1QbGN6gUOqbpGR7v9PYCimpe0nN2CfRczO6jqyMypjc4BElaYz3fCDHLl7sdVnVu9q6bAzI6pmuxbrtEBFF22MN7ZVjUz2Qf8aVV//48Pum9AgFo4ebzdSUX6fDd3H8T7YBPLFgs7reovdu1D/aiqqf0w/2JF6rIdhD+n6lTSkqp/DD4n6RsbLagHFEVd8B5XtQam/g6lZVV/py+7+127Xxfl850wg4Exs326Uz+wlM25AkMjm1bap+oHO9On2DRif74TZgAAQNKomQEAAEkjzAAAgKQRZgAAQNIIMwAAIGmEGQAAkDTCDAAASBphBgAAJI0wAwAAkkaYAQAASSPMAACApBFmAABA0ggzAAAgaYQZAACQNMIMAABIGmEGAAAkjTADAH1gZhfNzOu+JmP3CdgszN1j9wEAhkIWYC5Kkrtb5O4AmwYjMwDQP+PZ43zUXgCbDGEGAPrnYPZ4OmovgE2GMAMA/TOTPTIyAwwQNTMA0Cdm5pKW3X137L4AmwkjMwDQB2a2L/uWURlgwEZidwAAhkStXuZs/UEzOyjpcPZ0j6RPu/vyIDsGDDtGZgCgP+6qlzGzI5Im3f2oqrdsH5J0JELfgKFGzQwA9EFzvUwWZC67+5ns+XlJ+yTtd/cL8XoKDB+mmQAgUHO9TBZklty9vn7mgCQxxQT0H9NMABBuvV7GzE7o7iAjd18myAD5IMwAQLhavcxJST+T9GzEvgCbDmEGAMIdlLQkaX/2/HkzO2tm423aAOgTCoABIEBWL3Ne0qnsriVlIeaKqtNNe2P2D9gMGJkBgDB3rS+T1cbMS5qsH50xs0MD7huwKRBmACDMk9lj88q/e6S77l56UgD6jmkmAAiQrS9z13SSmV2UpNrxbFRmufkuJwDhGJkBgB5lWxVIrfdjmlc2OpN5kiAD5INF8wAgzLKqt2Q3Oy5pj5mdrnsOIAdMMwEAgKQxzQQAAJJGmAEAAEkjzAAAgKQRZgAAQNIIMwAAIGmEGQAAkDTCDAAASBphBgAAJI0wAwAAkkaYAQAASSPMAACApBFmAABA0ggzAAAgaYQZAACQNMIMAABIGmEGAAAkjTADAACSRpgBAABJI8wAAICkEWYAAEDSCDMAACBp/x9s/Lr68+D5zgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 525x525 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAH6CAYAAAC59R7OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAa6QAAGukBho3xbQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXjU1dn/8ffJJBP2HdkMSwDZlF0RWWRJQHEXxKVa28e19WltXUBQBBdQUFrb2kXa8vRnNxEQtaIgEIgaRHZQ9lVAZCesIcvk/P6YSQoIkSQzc2b5vK7re8XMTGY+thfx5nzPuW9jrUVERESktBJcBxAREZHopCJCREREykRFhIiIiJSJiggREREpExURIiIiUiYqIkRERKRMVESIiIhImaiIEBERkTJRESEiIiJloiJCREREykRFhIiIiJSJiggREREpExURIiIiUiYqIkRERKRMVESIiIhImSS6DiAlM8bYwD/mOQ0iIiKxxAtgrTXleRMVEVHC6/V6XWcQEZHYkJcXnL+XqoiIfHler9ebm5vrOoeIiMSI5ORk8oJQSWhPhIiIiJSJiggREREpExURIiIiUiYqIkRERKRMtLEySmRlZZ3xfUpKCo0bN3aURkREBIy19vtfJc4YY3I9Ho/X5/Od8fjo0aMZM2aMm1AiIhLVik5nWGuTy/M+WomIAh6Ph8zMzDMeS0lJcZRGRETETysREc4Yk6s+ESIiEkzBWonQxkoREREpExURIiIiUiYqIkRERKRMVESIiIhImaiIEBERkTJRESEiIiJloiJCREREykTNpkQkph09ksO6NbvZ+fUhdu86zP59RzmVU8CpU/l4PIaq1SpQrXolUprUomWr+lzSuh5Vq1V0HVskKqiIEJGY4iso5MtVu/j8s02sXLaDHdsPlurnjYHLOqZwdf/W9OpziQoKkRKoY2WEU8dKkQuzZdM+Zv1nNZ/M38DRIznfeT7BY6hXvzr1G1SnUiUvyRWS8PkKOXokh+zDJ9n59UEKCgrP+JnkConccEsnBt95OTVqVArXv4pIyAWrY6WKiAhXVERkZGSc8bimeIr4Vx0+zdzIjClL2bh+zxnPVayYROfLm3JZxxTatW9Ek2Z1SErynPe98nIL2LZ1P8sXb2fBvPVnrGAkV0hk6A+6MfQHV5CYeP73EIkWKiLihKZ4inxXXl4Bsz/4kulvLWHvnqPFjyd5PfTofQl90lrTqUsTvMllu2NrrWXzxr289eYXLPx0U/HjzVtexOMjr6VZ87rl/ncQcUlFRJzQSoTIf/kKCpk7aw3/+n8L2bf3WPHj9RtW56bBnek/sG3Q9zBs2bSPP/0mg69W7wIgMTGBhx/tx3U3dQzq54iEk4qIOKE9ESJ+y5ds543fzT/jNkPzlhdx+93duKp3Szye0J1YLyy0vD99OX+b9Cm5uQUA3Di4Ew8+0hdPok7KS/RREREnVERIvNuzO5s//XY+XyzcUvxY46a1+eF9Pbiqd0uMMWHLsmP7QcaMmMG332QD0OWKpjzzwo1UqOgNWwaRYFARESdUREi8KijwMePtZfzz/xYW/+2/Rs1K/OiBXqRd2y6kKw8lOXokh7HPvs/qFTsBaN8phedevkWFhEQVFRFxQkWExKNtW/bz6tgP2bp5PwAeTwI339aZO+/tTuXK5fqdFxQFBT4mjvuIBXPXA3Bph4t5fvytVKykQkKig4qIOKEiQuKJz1fIO1OW8uZfsyjI959IatOuIT97Ij3iTkT4CgqZ+NJHzJ+zDoAOnVJ4/pXBeL3q4SeRT0VEnFARIfHi293ZTBz3EWtWfwNAcnIiP364Nzfc0omEhPDteygNn6+QX788i3mz1wLQu18rhj97fcTmFSkSrCJCJbOIOGWtZfbML5n0u/nk5OQDcEnr+jz5zCAublzLcbqSeTwJ/HL4NRw/lssXC7fwScYGatepwoP/29d1NJGw0EpEhNNKhMSykyfz+O0rH5M5z7+3IMFjuOve7txx95VRdXTy1Kl8Rvzybdav+RaAh3/ej5uGdHacSuT8dDsjTqiIkFj19bYDjB31Pjt3HALg4sa1ePKZQVzSur7jZGVzJPskjz/yb77ZeZiEBMOYl2/h8itTXccSOScVEXFCRYTEooyP1/LbVz8m95T/6GaftDb8/In0qD/dsHvXYR596J8cP3aKipW8/PqPd9GkWR3XsUS+Q0VEnFARIbEkL7eAN343nw/fXwVAYpKHh3/el0E3dghr06hQWr1iByMfm4bPV0j9BtX59Rs/0ARQiTgqIuKEZmdIrDh04DjPP/MeG9b69w1cVL8aTz9/Y9TevijJrA9W85sJHwPQrn0jxv3qNh39lIiiIiJOaIqnxILNG/by3Mh3ObDfPzTriu6pPPH0tUEflhVJJr0+nxlvLwMgfdCl/HL4wJhZbZHopyOeccTj8ZCZmXnGYykpKY7SiJTOp/M3MHHcR8Wtq2+/pxs/vK9nzPdSuO8nV7NrxyGWLNrGnA+/4uKUmgz9QTfXsUSCSisREU57IiRaFRZa/vW3hfzzb58DkOT18Mvh19A3vY3jZOFz4kQuT/z032zfdgCAJ56+lv4D2zlOJRK8lYjoOYgdpYwx440xD7rOIRJOebkFjH/ug+IColbtyrzy2zviqoAAqFw5mTHjb6FW7coA/Prl2Sz9YpvjVCLBoyIihIwxnYFhQA3XWUTC5djRHEY+PpVP5m8AoGXrevxm0t20atvAcTI36tWvzouvDqFylWR8vkL/BNCVO13HEgkKFRGhNdV1AJFw2rM7m8d++u/i+RdX9WrJhN/eQZ26VR0nc6tZ87o8O+5mkrweTuXk8+yT01m2WCsSEv1URISIMWYYsNx1DpFw2bRhD4/99F/sCnSgvGlwZ0Y+fwMVKiQ5ThYZ2ndMYfS4m/F6E8nNLWDMiHdZEGj3LRKtVESEgDEmFWgOzHGdRSQcFn++lWE/n8LhQycBeOCRPjz08754PPoVc7ouVzTjhVcHU7FiEgX5PsY/9wF//WMmPl+h62giZaI/4aExHhjuOoRIOMydtYbnRs7gVE4+iUkeRoy5nltv76qeCOfRvmMKL//mdmrXrQLAtH8vYdST0zl2NMdxMpHSUxERZIGTGFOstdmus4iE2oy3lzJx3EcU+ixVqiQzbuIQevdr7TpWxLukdX1+O+ke2rVvBMCKpV/z8wf+wbYt+x0nEykdFRFBZIypAaRba6e5ziISStZa/t+fP2PS6wuAwBHO1+/gso5qgnahatWuzEu/Hsr1N3cEYM+3R/jlw/8k4+O1jpOJXDgVEcE1HnjAdQiRUPL5Cnl94lze+vsiABo0qsHE399J09S6jpNFn6QkD488lsYvhg8kyeshN7eAV178kMl/+gQ1ApRooCIiSIwxacAy3caQWJaf72PC8zOLp3CmtqjLq6/fSf2GaoVSHgOvu4xXX7+Ti+r5j8JO/ddiXhs/WxsuJeKpiAieh6y1k0Lxxnl5eSQnJ5d4iYRabm4+z498t7iJVLv2jRj/m9uLuzFK+VzSuj6/mXQ3LS6pB8DHH37FpNfnO04lUjIVEUFgjNFpDIlpp3LyGPPUjOKWzVd0T+XFV4dQpWoFx8liS42alXn5N0Npe5l/w+X701fw3jS1m5HIpSKinAKtrbHWbg3VZ3i9XnJzc0u8RELl5Mk8Rj35DiuX7QCgx9UteebFm9REKkQqV05m9LibaZRSE4A//34BmzfsdZxK5NxURJTfCGutViEkJh0/doqnH5/KV6t3AdAnrQ0jRt9AUpLHcbLYVq16RZ57+VYqVEzC5yvk1XEfkpdX4DqWyHeoiCiHwCpEmjFmy7ku/Kc1AEac9vgQh5FFLtjRIzmM+OVU1q/5FoD0QZfyxNPX4knUr41waJRSk/t/ejUAX287yNR/LXacSOS79NugHKy1y621Na21zc91AW8HXvrSaY+rh4REvMOHTjD80Sls3uhfRr/u5g78YthAtbEOs0E3dqBDJ3/vjan/WszBA8cdJxI5k34jiMgZDh44zvCfT2H71gMA3DK0C4/8Mo2EBLWxDjdjDA/+rC/GQO6pAt78y2euI4mcQUWEiBTbv/coT/7sLXYGJnHefnc3Hnikj+ZgOJTa4iLSr70U8M8p+Xa3WtFI5FAREVqpga/NnaYQuQD79x5l2KNT+PYb/3+k7rmvB/c+0FMFRAS460fdSfAYCgst099a4jqOSDEVEUFmjEkzxiwzxhwG0gIPPxjYVLms6EioSCTZv+8Ywx6dwp7dRwD48UO9uOve7iogIkS9+tXpm9YG8DehOnzohONEIn4qIoLMWjvXWtslsOHSnHY1DzyuzjESUfbvO8bwswqIoT/o5jiVnG3InZcDkJ/nUwMqiRgqIkTi2IH9x3jqtFsYP3pQBUSkappal25X+e+M/mfGCk4cV5M5cU9FhEic8hcQb7M7UEDc+0BPbr9bBUQkG/qDKwA4eSKPme+tdJxGREWESFw6eOA4T/3ibb7ZdRiAH97XgzvuudJxKvk+bS9rxKXtLwbgvWnL1cVSnEt0HUAuTFZW1hnfp6Sk0LhxY0dpJJodOnCc4Y9O4Zud/gLinvt6cOe93R2nkgs15K7L+Wr1Lg4dPMGCuesYMOgy15EkjmklIgr4fD569ux5xjV58mTXsSQKHTp4guG/eLu4gLj7x1dxlwqIqHL5lamkNK4FwPS3lmKtdZxI4plWIqKAx+MhMzPzjMdSUlIcpZFodST7JCMfe5tdgUZSP/hRd37w46scp5LSSkgwDL7jcl6bMJsd2w+yZNE2ruie+v0/KBICRlVsZDPG5Hq9Xq/GfUt5HD92iqd+8TZbNu0D4I57ruSH9/dQH4golZdXwI+G/pnDh07QvlMK439zu+tIEmWSk5PJy8vLs9Yml+d9dDtDJMadPJnHqGHTiwuIW2/vqgIiynm9idw0xN+3bvWKnWxcv8dxIolXKiJEYtipU/mMeeqd4nHe193cgft/erUKiBgw6KYOVKiYBMA//2+h4zQSr1REiMSovNwCXnj6Pb5cuQuA9Gvb8dNfpKmAiBFVq1bgxls7AbD4860sX7LdbSCJSyoiRGJQQYGPcaP/U/wflt79WvHosIEa5x1jht7djRo1KwHwxu/mk5ub7ziRxBsVESIxxldQyIQXZvLFwi0AdO/ZgiefGYTHoz/usaZy5WTufaAnADu2H2TS6wvcBpK4o98qIjGksNDy6/Gz+HT+RgC6dmvKU2OuJzHR4ziZhMrA6y4rnqnx4XurmPH2UseJJJ6oiBCJEdZafv+rucybvRaADp1SeObFm/B61Q4mlhljeHzkNdRvUB2ASa8v4L3pmvIp4aEiQiQGWGuZ9PoCPnx/FQBt2jVk9Eu3kJyc5DiZhEPVahV5+bWhXFSvKgB/+k0GH7yrAV0SeioiRGLAm3/J4t2pywBo2aoeL7wymIqVvI5TSTjVa1Cdl167nTp1/YXE7381l9kffOk4lcQ6FREiUe7fby7irb8vAqBpszq8+OoQKlcpVxM6iVING9Xg5deGUrtOFQBe/9UcNm/Y6ziVxDK1vY5wRW2vMzIyznhcUzwF4J0pS/nz7xcA0CilJq/87g5q1qrsNpQ4t2P7QX7+4N/JPVVASpNa/OH/7tXmWjmD2l7HEU3xlHP58L1VxQVE/QbVeenXQ1VACACNm9bmwUf6ArDz60PM0m0NCRFt244CmuIpZ5s/Zx2v/2oOAHXqVuXl14ZS96KqjlNJJLnmhvbMfH8VWzft419/W0j6te200VaCTrczIpymeMrZFmVt4YVn3qXQZ6lRsxKvvH4HF6fUch1LItDSL7Yx6snpADw6bADXXN/ecSKJFLqdIRKHVi3fwbjR71Pos1SpkszYiUNUQMh5dbmiKU2b1QHg3anL0F8aJdhURIhEifVrv+W5ETPIz/ORXCGR5ycMJrXFRa5jSQQzxnDzbV0A+HrbQdZ8+Y3jRBJrVESIRIHtW/cz6snp5OTkk5jkYfS4W2hzaUPXsSQKXN2/VXHPEPWNkGBTESES4XbvOszIx6Zx/NgpEjyGp0ZfT6euTVzHkihRoaKXPv1bA/Dpgg2cOKH9VRI8KiJEItj+fccY8dhUDh86AcAvh19Dj94tHaeSaDPw+ssAyD1VQOa89Y7TSCxRESESobKzT/L0Y1PZt+coAD95tB9p17RznEqi0SWt69M01b/BUrc0JJhURIhEoBPHcxn1xHR27jgEwA/v78mNgzs7TiXRyhjDwOv8qxEb1+9h25b9jhNJrFARIRJhTp3KZ8xTM9i80T/zYPAdXbnjnm6OU0m06zegLYlJ/tbXs2dqNUKCQ0WESATJz/cxdtR7fLV6F+DvOnjfT67GGOM4mUS7atUrclWvFgBkfLyWvLwCx4kkFqiIEIkQPl8hr7w4k6VfbAegd79W/O9jaSogJGiKbmkcO3qKzz/d7DiNxALNzogSWVlZZ3yvKZ6xxVrL716dw6fzNwJwRfdUnnxmEB6P6nwJno5dmnBR/Wrs23OUD95dSe9+rVSkSrnoN1QU0BTP2Gat5c+/X1B8n/qyjhcz8vkbNLpZgi4hwTDoxg4AfLVqF8uXbHcbSKKeBnBFuKIBXBkZGWc8rpWI2PHPvy3kH5MXAtCydT1e+vVQKlcu10wckfPKOZnHfXf9hcOHTnJR/Wq8/tcfUrVqBdexJMyCNYBLRUSE0xTP2Pbu1GW88bv5ADRuWptXfncH1apXdJxKYt2CeesZ/9wHAFzVqyXPvHijbmvEGU3xFIlycz76qriAqN+gOuN+dZsKCAmLPv1bc02gi+XCTzcx4+1ljhNJtFIRIeLA559u5rUJswGoXacK4359G7XrVHGcSuLJQz/vVzwm/K9/ymT1ih2OE0k0UhEhEmarV+zgpef+Q6HPUrVaBcZOHEKDhjVcx5I4U6FCEs+8eBOVqyRT6LNMHPcRp07lu44lUUZFhEgYbdqwhzEj3iU/z0eFikk8P2EwTQJ/GxQJt0YpNXn0yQEA7Nt7jPemLXecSKKNigiRMNm14xDPPDGdnJN5JCYmMOrFm2jdtoHrWBLneva5hMs6XgzA+9OXU1Dgc5xIoomKCJEw2L/3KCMfm8rRIzkYA8NGXUfny5u6jiWCMYZbb78cgEMHT/BZ5ibHiSSaqIgQCbEj2Sd5+vFp7N93DID/fTydXn1bOU4l8l+XX9mM+g2rA/C+bmlIKaiIEAmhkyfzGDXsvyO9f/Rgr+KOgSKRwuNJ4IZbOgGwbs1uvt52wHEiiRYqIkRCJC+3gOdHvsum9f6R3rfe3pWhP7jCcSqRc+s/sC2Jif7/JMz56CvHaSRaqIgQCQFfQSHjn5/JquX+s/fpgy7l/p9qpLdEruo1KtHtquaAf1S4NljKhdAUzyihKZ7Rw1rL7ybOYeGn/g1q3Xu14NEnBqiAkIiXdu2lZH2yicOHTrL0i+1c2aO560gS4bQSEQU0xTO6TP7TJ8UTOTt0SuGpZ6/Hk6g/ahL5unZrSs1alQCYq1sacgG0EhEFPB4PmZmZZzyWkpLiKI2UZOq/FjPt30sAaNmqHqPG3Yw3WX/MJDokJnrom96Wd6Ys5YuFWziSfZLqNSq5jiURTL/dokSPHj1cR5DvMeuD1Uz+0ycAXNy4Fs+/MlgjvSXqpF/bjnemLKWgoJCMj9dxy9AuriNJBNMaa5AYYx40xsw57VpmjHnQdS4Jj88WbOR3r84BoE7dqoydOIQa+hucRKGmqXVp1aY+ANPfWkJuruZpyPmpiAgCY8wcAGttetEFPAC8YYzZYozRdKUYtmLp14x/YSaFhZZq1SsyduIQLqpXzXUskTK760fdATh44DjvT1/hOI1EMhUR5WSMGQ9MtdZOOv1xa+1yYAKQCkx1kU1Cb8Pab3n+6XcpyPdRsWISL7wymMZNa7uOJVIul1+ZSrv2jQD495uLOLD/mONEEqlURJTfg8B4Y0zaOZ5bEviaptWI2LNj+0FGDZvOqZx8EpM8PDvuFi5pXd91LJFyM8bwwCN9SEgw5JzM44+vZbiOJBFKRUT51Qhct53juezT/rlWeOJIOOzdc4SRj0/l2NFTJCQYnhp9PR27qG+HxI5WbRoUt8Je+OkmFn222XEiiUQqIsovHf9ti+HneK549cFauzVsiSSksg+f4OnHpnFw/3EAfv7kAHr0buk4lUjw/fCBntSpWxWAN/+ahbXWcSKJNCoiyslaO9daO9xam32Op28PfJ10juckCp04kcszT0znm12HAbjvJ1cz8LrLHKcSCY1Klbzcee+VAGzbsp8vFurvQnImFREhEtgDMQTYyrlXKSTK5Obm89yIGWzZtA+A2+66giF3Xu44lUhopV3Tjtp1qwDwzpQl3/NqiTcqIkIgUEDMA5YDXc6zSiFRxFdQyMtjPuDLlbsAuOaG9vz4oV6OU4mEntebyI2BvRFfrtzF19s1Jlz+S0VEkBhjahhj3gj0jDiMfwWifzAKiLy8PJKTk0u8JHQKCy2vTZjNoqwtAPS4uiX/+1iaBmpJ3Bhw3aXFY8I/em+14zQSSVREBIm1Ntta+1Cg2ZQB5gCHA30kJEpZa/nLHxYwd9YaADp1bcKwUdfh8eiPjsSPGjUr0+PqSwCYO3sNp3LyHCeSSKHfhCESaD41HBhW1NGyrLxeL7m5uSVeEhpT/vEFM95eBkCrNvUZ9eJNeL0aOSPxZ9BNHQA4cTyXzHkbHKeRSKEiIrSKTmWkGWOGOE0ipTbzvZX8vz9/BkBKk1o8P2EwFSt5HacSceOyDhcXd2Od+f5Kx2kkUqiICKHAfoiiM1EjXGaR0smct57f/2ouABfVq8rYibdRrXpFx6lE3DHGFK9GbFq/l43r9zhOJJFARUQ5GGPSAhM7S5rWWVREdA5HJim/ZYu38erYD7EWqteoyNhf3Ubdi6q6jiXiXP8BbUmu4L+dN/NdrUaIiojyGg+kAW+U8Jqidtfq0hIF1n21mxeeeY+CgkIqVvLy4itDuDhFHctFAKpUrUCf/m0A/2rdsaM5jhOJayoiymdp4Ou0El6TGvi6PMRZpJy2b93Ps8PfIfdUAV5vImNevoUWreq5jiUSUa6/uSMAubkFTJ+y9HteLbFORUT5zAEmWWvPNXwLY0xn/js/46WwpZJS27M7m6cfn8bxY6dI8BhGjLme9h1TXMcSiTgtWtXjiu7+vxu9N2052YdPOE4kLqmIKAdr7TQAY8yw87ykqEfEcGutViIi1KGDJxj5+DQOHfT/Mvzl8Gu4smcLx6lEItc99/UA4FROPlP/udhxGnFJRUQ5WWsfAjDGLDPGDDPGDDHGPGiMWQZ0BdKttRPcppTzOX7sFKOenMa33/gbiz74v31Ju6ad41Qika3FJfXo2cfffOqDd1dx9Ij2RsQrFRFBYK2dYK3tgn/fQ9EuvAestTWttXMdRpMSnDqVz5gRM9i6eT8Ad/7wSm4Z2sVxKpHocMc9/umeeXkFzPnoK8dpxBUVEUEUGAs+KXDp9kUEKyjw8dLo/7Bm9TeAvxtf0RKtiHy/5i0vou2lDQGY+d4qCgut40TigooIiTuFhZaJ42ax+HP/qdve/Vrx01/010AtkVK6LnBS49tvslm+ZLvbMOKEhgBEiaysrDO+T0lJoXHjxo7SRC9rLX/6bQYL5q4DoMsVTXni6UEaqCVSBj37XMIbv5vP0SM5fPj+Krp2a+Y6koSZfnNGAZ/PR8+ePc+4Jk+e7DpWVPrX3z7nP++sAKBNu4Y888KNJCV5HKcSiU5ebyIDBl0KwBcLt3DwwHHHiSTctBIRBTweD5mZmWc8lpKiHgal9f705fzj/xYC0LRZHZ4bfwsVKmqglkh5XHtDe6b9ewmFPsvsmV9y173dXUeSMFIRESV69NCmv/KYP2cdf/xNBgD1G1TnxYlDqFpNA7VEyqvhxTXp1LUJK5Z+zUfvr2boD64gMVGre/FCtzMk5i3+fCsTx30EQM1alRg7cQi161RxnEokdlx/i3+D5YH9x8jK3OQ4jYSTigiJaV+t3sXYUe/j8xVSuUoyL746hIYX13QdSySmdLuqOfUbVAdgxtRljtNIOKmIkJi1dfM+xjw1g7y8ApKT/QO1Ultc5DqWSMzxeBK4cXBnADas/ZZ1X+12nEjCJeRFhDGmmjGmaeCqFurPEwHYveswzzwxjRPHc/F4Ehj5/I1c2v5i17FEYtaA6y6lYiX/RuUZUzXdM14ErYgIFAv3G2OmGGM2GWMOGmN8wGFgS+A6bIzxBZ6fYoy5T4WFBNvBA8cZ+dhUDh86CcDjI68tnjooIqFRuXIy11x/GQBZmZvYu+eI40QSDsba8rUqNcbcCjwEpAEG2Ip/hsQhIBs4GPhaA6gd+FoL6AykAhaYC7xsrZ1frjAxyBiT6/V6vbm5ua6jRIVjR3MY9rMpbN92AICfPNqveJlVREJr77dH+J87/0JhoWXwHV25/6d9XEeS80hOTiYvLy/PWptcnvcp8xFPY0xTYBr+YmASMBSYa6294PLTGFMdSMdfgMwLTL6831q7qqy5JH6dyslj9PAZxQXE3T++SgWESBjVa1Cd7r1akJW5iVkffMldP7qKSpXUiyWWlel2hjHmSfyrDVOAmtbah62100tTQABYa49Ya6dZax8GagJTgfnGmMfLkkviV36+jxdHvc+6Nf4NXTcO7sRdP1LTG5Fwu/k2/yTcE8dzyZi9xnEaCbVSFxHGmD8BXa21tay1r5S2cDifQEExwVpbC+hmjPljMN5XYp/PV8irYz9k2eLtAPRNb8NDP+ungVoiDrS7rBHNW/pPQb3/zgrKe8tcIlupiojACsQya+3tIcoDgLV2KLDVGPNEKD9Hop+1lj+8No9PMjYAcEX3VB4bcQ0JCSogRFwwxnDj4E4A7Pz6EKuW73ScSEKptHsipllrt4UkyVmsta8YYzQSLkBTPM/t73/N4sP3/Fto2rVvxIjnblDLXRHHru7fmr/8IZNjR0/x/jvL6dhFv6tiValWIsJVQLj6vEilKZ7nNuPtZfz7zUUApLaoy5iXbqFChSTHqUQkOTmJgdf5j3t+kbVFxz1jmAZwRQFN8fyuubPWMOl1/4ngBo1q8MIrQ6hStYLjVCJS5PqbOzL9rSUUFjEGv1oAACAASURBVFo+fG8VP36ot+tIEgLl7hPxnTf0N496EP9xz5VBffM4pD4R37Xos828MOo9Cn2WWrUrM/H3d1K/YQ3XsUTkLM+NmMGirC1Uq16Rv09/CK9Xf2+NFMHqExGKttd/ASYAS8/VjdIYMzjQrbJfCD5bYtzqlTsZN+Y/FPosVapWYOzEISogRCJU0XTPo0dyNN0zRoWiiDgEDAfOuQoR6CdxO1AzUEw0DUEGiUGbN+5lzFMzyM/zkVwhkefH30rT1LquY4nIeXTq2pQGjfxF/gfvamE6FoVkAFegf0RXa+3REl5TVEyMD0UGiS27dh7imSemkXMyj8TEBEa9eBNtLm3oOpaIlCAhwTDoxg4ArP3yG7Zt2e84kQRbKIqIqcaY+8/3pAZuSWnt33eMpx+bypHsHIyBJ54eRJcrdPpXJBoMGHQpSV7/sWutRsSeoBcR1tp5gAkM5ipmjKlujNnEmZM8Dwb78yW2HD2Sw9OPT2Xf3mMAPPJYGlf3b+04lYhcqGrVK9K7bysAMj5ey4kT2iQeS8o6O6NpSc9ba/8M1D5r8+RUoDkwD1iBf+LnVPwnOUS+4+TJPJ4dNp2dXx8C4If39+S6mzo6TiUipXXdzf4/t6dy8pk/Z53jNBJMZV2J2BpYTZhljLnvXEVFoJBobozpF5jW2RWoYa0dENgv0SIwuEtdSOQ78vIKeOHpd9mwbg8Atwztwh33dHOcSkTKonXbBqQG5mnMfHel5mnEkPLczjDAAPxjwLcYYw4aY94yxtxfVFQECokBwAPASyVttBQp4vMVMuGFmaxctgOAtGvacf9P+2iglkiUMsZwXWCD5fatB1i/5lvHiSRYylpEWCAV6AI8BWTgH+U9FHiD04oKYAv+Exhzyh9XYp21ltcnzik+U35lj+b8YthADdQSiXJ90ttQsaK/Lf3M97TBMlaUtYjYZq3dbq1dETjOmW6tTeDcRcWfCOx/MMb80Rhzq05oyPn83xufMuuDLwG4rOPFjBhzA57EkJxEFpEwqlTJS98BbQH4ZP4Gjh7JcZxIgiHoba/PeHNjOgFpQHrgK/hXMQC24l+dmKT22OdX1PY6IyPjjMdjcYrntH8v5q9//ASAFpfU4+XfDKVy5XJ1ZBWRCLJ18z4e+Z83AXjgkT7centXx4niV7DaXoe0iPjOh/mLinTgdqBT4OHD1traYQsRZYwxuR6Px+vz+c54fPTo0YwZM8ZNqBCYPfNLXhs/G4BGKTV59fU7qFGzsuNUIhJsj/3kX6xbs5tGF9fkz//8H+11ciRYRURYp6FYa1fgP945AcAY0z+cnx+tYn2KZ9Ynm/jtKx8DULtuFcZOHKICQiRGDbqpA+vW7OabXYdZtXwnHbvE1opqvAnrSoSUXqxP8Vy5bAejhk2nIN9H1WoVeOX1O2jStI7rWCISIrm5+dx96xscP3aKXn0vYeRzN7qOFJcieYqnyAXZsO5bnh85g4J8HxUqJvHChMEqIERiXHJyEunXtgNg4SebOXTwhONEUh6lKiKMMWFtF6gJn7Fr59cHeXbYO+Tk5JOY5OHZsTfTqm0D17FEJAyKhnL5fIXMnvml4zRSHqVdibjcGPNESJKcJfA5Q8LxWRJe+/YeZeRj0zh6JIeEBMPwUdfRqWsT17FEJEwublyLDp39eyFm/WcVPl+h40RSVqUqIgIdKFsYY2YbY6qGIpAxppox5mOgubX21VB8hriTnX2Spx+fxoH9/oFaP3sinZ59LnGcSkTC7bqb/KsR+/YeY9nibY7TSFmVek+EtfZhYDrwtTFmXLAaRwWKh5eBbcDb1tqfBON9JXKcOJHLs09OZ9cO/0Ct/3m4N9dc395xKhFx4cqeLahZqxIAM99b5TiNlFWZNlZaayfhH6jVFf9o7ynGmFtKW1AYY5oGOlhOAQ7j7x3R1Vr7l7LkksiVl1vACyPfZdOGvQAMufNybrvrCsepRMSVpCQPA667DIAln29l7x7NYoxGZe4TYa3dCgwwxqThH+c9HbDGmGxgKf6OlNnAwdN+rDb+mRs1CEz1xN8Sey4wwFo7r6x5JHL5Cgp5+bkPWLViJwADr7uM/3m4t+NUIuLatTe05+1/fIG1MOs/X3LvAz1dR5JSKnezKWvtXGBuYNx3Ov55GamBfz6XbPwFxlT8ba/nahx47LLW8ptXZvP5Z5sB6NG7JT97PF1d6kSEevWrc/mVqSz+fCuzZ67mBz/uTmKix3UsKYWgdawMFALTAlexQHFRCzikYiG+WGv5yx8ymfPRGgA6dG7MsFHXaaCWiBQbdFMHFn++lcOHTvL5Z5vp1aeV60hSCiH/bW6tPWKt3aYCIv68/c/FvDNlKQAtW9fj2XE3400Oa6d1EYlwXbs1o+5F/sN+M9/VBstoo9/oUSIrK+uM7yN9iueH76/ib5M+BSClcS1emDCYSpW8jlOJSKTxeBK49ob2vPnXLFYt38HuXYdpeHFN17HkApV7JcIY088Y84Qx5tZgBJLv8vl89OzZ84xr8uTJrmOd16fzN/D6xDkA1L2oKmMnDqF6jUqOU4lIpBpw3WUkJPj3SX384VeO00hplHklIjDWey7+ExZFjx0GHrDWzghCNgmIpimey5dsZ8ILM7EWqlWvyLhf3UbdekFpJSIiMap2nSpcfmUqXyzcwscffsU9/9NDe6eiRHluZ8zjtAIioBYwzRjzoLX2r+V4bzlLjx49XEf4XuvXfssLz7xHQUEhFSt5efHVwVzcuJbrWCISBa65/jK+WLiFw4dOsPjzrXTv1cJ1JLkAZSr1Ap0lAdKstQlATaAL8ApwFJhkjNEwhDjy9bYDjHpyOqdy8knyehjz0s20bFXfdSwRiRKXX5lKrdqVAZg1c7XjNHKhyrpelAbcb63NgOITGCustcOBZkAG8EaQMkYFY8yDxphlxpgtgWuqMSYuBojt/fYITz8+jePHTpGQYBgx+gbad4rcTZ8iEnk8iQmkD7oUgKWLthXP15HIVtYiopm19p1zPWGtzbbWpuMf1NWv7NGihzHmDfwDw7pYa5vjX5U5BEw1xswxxpx92ydmHD50gpGPTeXggeMA/GL4QC1DikiZDBjkb4NdWGiL+8tIZCtrEXEh7QafAm4r4/tHDWPMMCA7sAoDFBdSDwGT8K/a/NlVvlA6cTyXZ56czu5vsgF44JE+pF97qeNUIhKtGjaqUTwi/OOZX1JYaB0nku9T1iLi0Pe9wFo7DX/761j30OkFxFmKHh8SmDESM3Jz8xnz1Ay2btoHwO33dOPW27s6TiUi0e6a6/2rEXu+PcKq5Tscp5HvU9YiQuUhENjzMO18z1trs4HlgW/PN0sk6hQU+Bg3+j98tXoXAINu7MC992twjoiU31W9WlK1WgUAZn2gDZaRLtQHcWN9ylIqMMwYM7WE12wNfO0chjwhV1ho+fXLs1m80P+v1btvK376y/4aqCUiQeFNTqTfgLYALPx0M0eyTzpOJCUpaxFxoYf/s8v4/tFmSAmbJ4sej/r/Lay1TPrdfDI+XgtA58ub8sQzg/B41BRGRIJn4HX+WxoF+b7i3zcSmcr627+mMeaAMeYlY8ytxph4bUk4Cf/tjOGBWxfnUrRRYOt5no8a/35zEe9N99+dad2uAc+8cCNJSRrbKyLB1ax5XVq1bQDArA++xFrdQY9U5fkrZC1gGDAVOGyM2WSM+eNZRUVM/z8fOIVxm7V2wrmeN8Z05r8rES+FL1nwfTBjBX//q38IWJNmtXl+/K1U1EAtEQmRog2WO7YfZP2abx2nkfMpTxExHP/RxW349z40Bx7itKICSDPG3GeMaXq+Nwm8LlaNCHydUMJKxffKy8sjOTm5xCuUFsxdxx9emwdAvfrVePHVIVStVjGknyki8e3qfq2pWDEJ0AbLSFbWIiLbWvuKtfZha20L/G2vb8O/vH96UVEz8NgWY8xBY8wUY8z9ZxUVzcucPoIFViGGAHNLOAIa8ZYs2sqrYz/CWqhRsxJjJ95GnbpVXccSkRhXsZKXq/u3BiAzYz0nTuQ6TiTnUtYi4u3Tvwm0vZ7+PUVF0WNv8N+iYjaxe8tjKrA80L2zXLxeL7m5uSVeobD2y28YO+p9fL5CKlX28uKrQ2iUUjMknyUicraB17cHIPdUAZnz1jtOI+dSpiLCWvvw9zx/oUVFzPROOJ0xZg7+0xj9XWcpq21b9jN6+Dvk5hbg9Sby3Mu30rzlRa5jiUgcadWmPk2b1QFg9gdfOk4j5xKWs3klFBXzwvH54RSYo1EL6F+efRAu7f4m2z9Q63guCR7DyOdv4NIOF7uOJSJxxhjDwMAGy43r97B18z7HieRsTg74n1ZUxNRKhDHmQaBrYBBXVBYQhw4c5+nHp3L40AkAHh9xLd2uisltKyISBfoPbEuS13+UfPZMrUZEmkjoEnTEdYBgCMzGeMha2+U8z48Pc6RSO3bsFE8/MY09u/3/lzz8837FneNERFyoWq0iPXq3BCBj9lpyc/MdJ5LTRUIREbUnF4oETmIML6GAqMF/+0VEpFM5eYwZ/g7btx4A4K4fdeemITHRqVtEotzA6/wbLI8fzyUrM5a7AkQf50WEtTaqx2QHCoTxnDX23BhTI3Cl4u8XscxFvguRn+9j7LPvs/ar3QDccEtH7v7xVY5TiYj4te+UQv2G1QH1jIg0zouIGLAMSMPfYMsWXcDhwLUFf2fPiGx7XVho+dVLH7H0i+0A9ElrzcOPaqCWiESOhARTPE/jy5W72P1NVG45i0kqIsohsM8h9QJfvjSUWcrCWssfX5vHgrn+89dduzXj8ZHXkpCgAkJEIkvaNe2KfzfNnfWV4zRSREVEOVhrh1trzQVeEVc6/+P/FvLBuysBaHtZI55+4UYSEzVQS0QiT526VenUtQkA82atobAwVvsURhcVEXHqvWnL+dffPgegaWodxrx8CxUqJDlOJSJyfunXXgrAvr3HWLV8h+M0Aioi4tK82Wv4028zAKjfsDpjXx1C1aoVHKcSESlZ954tqFLFP3Bwzke6pREJEl0HkAuTlZV1xvcpKSk0bty41O+zKGsLv3p5FgA1a1Vm3K9uo1adKkHJKCISSt7kRPqkteGDd1eSlbmJE7/MpXKV0E4xlpJpJSIK+Hw+evbsecY1efLkUr/Plyt38tLo/1Dos1SpkszYiUNo0DCi21eIiJwhfZD/lkZeXgGZGRrK5ZqzlQhjTD/8DZhqBb5mA1uttRmuMkUqj8dDZmbmGY+lpKSU6j22bNrHmBEzyMsrIDk5kefG30qz5nWDGVNEJORatqpH02Z12L7tAHM++opBN3ZwHSmuhbWIMMY0xT/Js2i65dlnCW2gP8Ec/C2kvw5buAjXo0ePcv38m3/5jJMn8vB4EnjmxZtoe1mjICUTEQkfYwxp17bjL3/IZP2ab9mx/SCNm9Z2HStuhe12hjGmGbAcf9Olh4EB+EeBn34NDTx3DFhujFGJGSTDR19Pp65NeOKZQXTt1sx1HBGRMus3oC0JnqKeEWscp4lvxtrwnLU1xkwBnrLWbrvA16cCL1lrbw9tsshmjMn1er3e3Nzccr+XtVadKEUkJjw3YgaLsrZQq3Zl3pz6EJ5EbfErjeTkZPLy8vKsteXamRrO/9XNhRYQANbarfjbRkuQqIAQkVhRtMHy0METLFuy3W2YOBbOIuJgGX5GLclEROQ7ruieSvUaFQH1jHAprCsRxphbSvHi+/nuxksRERESEz30TW8L+PvfHD2S4zhRfArn6YzhwDxjzAT8myuzgUOBr3Dmcc+ugee6hDGfiIhEkQGDLuXdqcsoyPexYO46bhzc2XWkuBO2IsJaewToaox5EBgC3Hael24FXrbWvhKubCIiEn2aNa9Li0vqsXnjXuZ8tEZFhANhbzZlrZ2Ev1dE0bHPIocChYaIiMgFSb+2HZs37mXzxr1s27JfTfTCzOmZGGvtttMuFRAiIlIqfdPbkJjkAeDjD7XBMtwi+mBtYHOliIjIOVWtVpHuPZsDkPHxWvLzfY4TxZdIn+KZDvzFdYhIEKwpniIisSbtmkv5dP5Gjh7JYcmirVzVq6XrSHEj3LMzqpXyR7RLhv9O8Tzd6NGjGTNmjJtAIiIRpMvlTalVuzKHDp5gzodfqYgIo7AVEcaYpUCncH1eLAnGFE8RkVjlSUyg/8B2TP3XYhYv2srhQyeoWauy61hxIWyzMwCMMXPwt7JeciEvB4Zba+N6PFswZ2eIiMSqXTsO8cDdkwF44JE+3Hp7V8eJIluwZmeEe0/Eg8CwC+0BYYxJC3EeERGJARc3rkWbdg1Zt2Y3cz78iluGdtG8oDAI6+mMwACu5qX4ka2hyiIiIrGlaCjX9m0H2LRhr+M08cHFEc/zdar8Dmvtw6EMIiIisaN331YkJ/sX2OeoZ0RYhL2IUFMpEREJhcpVkrmqt/9kxoK568jLLXCcKPZFdLMpERGR0ii6pXH8eC6LsjY7ThP7gl5EGGPuN8a8VIaf62iMubUMvSREREQA6NCpMRfV9/9nZM5HuqURaqFYiXgYOHihLzbGNDPGbAKWAdOAw8aYDiHIJSIiMS4hwZB2TTsAli/5mgP7jzlOFNtCUUTMtda+WorXTwVSgQFATeAnQEYIcomISBwoKiIKCy3z56xznCa2hWRPhDGmWuDWxGxjzBJjzFvnWl0wxlTH39p6grV2nrX2SGBU+MsaviUiImXRoGEN2rVvBMDcWWsIZ1PFeBOKIuIlYBv+FYZ0oAswFFhujHn8rNcOBSww5azHJwV+VkREpNSKViN2bD/I5o3qGREqoehYORR/y+qHgaVANlADf1HwtDFmnrV2ZeC1RY2nzmgqZa09YtRq7Aya4ikicuF69WnFH1/LIC+vgLmz1tCyVX3XkWJSKFYihgDNrLV/ttausNZuC3ydAKQBt5/9A9bao+d4n0MhyBaViqZ4nn5NnjzZdSwRkYhVuUoy3Xu1AGDB3PXk5/scJ4pNoViJ2Ha+hlLW2uXGmAdOe6hGCe+jm1gBmuIpIlJ6ade0I3Peeo4eyWHpom3FRYUETyiKiJrf8/yFTuWsVd4gsaRHjx6uI4iIRJVOXZpQs1ZlDh86wdxZa1REhEAobmccNsb8wRhT9fQHAyc2ZgNvnfZwraLnznptM/z7KkRERMrEk5hAvwFtAFj8+RaOHslxnCj2hKKIGA7cAWQbYzYFjnhuAg7j31yZDWCMGYx//8Q2vjuU62XgTyHIJiIicaT/QP8pjYKCQjLnrXecJvYEvYgI7IdIxd8wqjn+I57Ngen4b3V0McZ8DDwEDLfWtgBGGGMeN8b0M8b8EehsrVXDKRERKZdmzevSvOVFgL9nhARXKPZEYK3NJtDnwRjTCdh62mbLVwLX6YYCc/FvtNyCv3uliIhIufUf2I4tm/axcf0edmw/SOOmF7o1T75PyKd4Bo53ljj+21q73FpbC2hurW1prd0W6lwiIhIf+qa3JsHj32Y3b7ZWI4IpLKPAAxM6+xljOpb0OhUPIiISbDVqVqZrt2YAZHy8jsJCdRAIlpAWEcaYJ4wxPvwTOucAy4wxPmPMH0L5uSIiIqdLC2ywPLD/GKtX7HScJnaEcgDXZmAC/qOaR067DPCwMWajMaZJKD5fRETkdN2uak6lyl4AMuasdZwmdoRqJSIDWA50sdYmWGtrnXYl4N9I+TUwLUSfLyIiUsybnEivPq0A+GzBRk6dynecKDYEvYgIjPCeYq0daq1dca7XWGunWWvT8d/eeCLYGURERM7Wb2BbAHJO5rEoa4vjNLEhFEc8b7PWDryQF1prHzbGnD0GXM5BUzxFRMrn0vYXc1G9quzbe4z5H6+lT//WriNFvVDczijxOOc5qL3199AUTxGR8ktIMPRN969GLF28jezDJxwnin6hWIko7dkZnbX5HpriKSISHP0GtGXKP76g0GfJnLeBm4Z0dh0pqoWiiLjglQVjTPXSvD6eaYqniEj5NW5am5at6rFpw14yPl6rIqKcQnE7Y0op9jm8zZlTPaOeMWaqMSbNdQ4RETm3olsaG9fvYdeOQ47TRLdQDOCazn8neN4X6FZZLXA1DXz/kjHmYOD17wQ7Q7gZY1KNMQ8aY7bgn0xaw3UmERE5tz5pp7fBVs+I8ghJnwhr7UPASuDP+LtVHg5cWwLfDw98HRqKzw8XY0xaoHB4I/DQVpd5RETk+9WsVZnOXZsCMH/OWrXBLoeQtb221t6GfxpnBv/tVGmAbcBD1toB3zeYK9JZa+daa5tba9OttZOAbNeZRETk+/Ub4L+lsXfPUdZ++Y3jNNErJKPAi1hr5+If8Y0xpnq0Fw0iIhIbuvdqQcWKSeTk5DPv47Vc2uFi15GiUlimeAKcr4AwxjQNVwYRERGAChWS6HH1JQB8On8DebkFjhNFp7AVESWY4zqAiIjEn6I22CeO57L4c21pKwunRUSgT0SqywwiIhKf2ndMoXadKgBkfKxTGmVR6iLCGNPMGOMLxgXogK6IiDjh8STQN70NAEsWbeXokRzHiaJPqTdWWmu3GWOKukyWd/2nNlCtnO8hIiJSJv0GtGXav5dQUFDIJ/M3cP3NHV1HiirlOZ3R2Vq7srwBippOyfnl5eWRnJxc4mtyc3PDlEZEJHY0a16XZs3rsm3LfjJmr1URUUpl3RORHYwCIuBwkN5HRESk1Ip6Rqxbs5vdu/SfpNIoUxFhra0VrADW2hbBeq9Y5fV6yc3NLfESEZGy6ZvehqKb9PPnrHMbJspEwhFPERERZ2rXqULHLk0AmPfxWqxVG+wLpSJCRETiXtEtjW+/yWb9mm8dp4keKiJERCTu9ejdkuQK/rMG6hlx4VREiIhI3KtYyctVvVoCkJmxnvx8n+NE0UFFhIiICP+9pXHs6CmWLtrmOE10UBERXEUtvC93mkJEREqtU5cm1KxVCYCMObqlcSFURJSDMSbVGLPMGLPFGHMY6Bx4apgx5nDgOQ0YExGJAp7EBPqk+dtgf7FwC8ePnXKcKPKpiCgHa+1Wa20Xa21za21Na6057aoZeC7ddU4REbkw/dL9tzTy83x8umCj4zSRT0WEiIhIQPNLLqJx09qATmlcCBURIiIiAcYY+gc2WH61ahd79xxxnCiyqYiIEllZWWdcO3bscB1JRCQm9QmMBwfInLfeYZLIpyIiCvh8Pnr27HnGNXnyZNexRERi0kX1qtGufSMAFsxVEVGS8owClzDxeDxkZmae8VhKSoqjNCIisa9PWhvWrP6GbVv28/W2AzRpVsd1pIhkNGgkshljcr1er1eTOkVEwic7+yQ/uOWPFPosd9xzJfc+0NN1pKBKTk4mLy8vz1qbXJ730e0MERGRs9SoUYnOXZsCsGDeOk32PA8VESIiIufQJ601AHt2H2HDuj2O00QmFREiIiLn0L1XS7xe/9bBBXPWOU4TmVREiIiInEOlSl669fCPRPpk/np8vkLHiSKPiggREZHzKJqlcfjQSVav2Ok4TeRRESEiInIeXbs1o3IV/wGGBXN1S+NsKiJERETOw+tNpEfvlgBkfbKJvLwCx4kii4oIERGREhTd0jhxPJelX2xznCayqIgQEREpQftOKdSsVQnQLY2zqYgQEREpgceTQO++/p4RX2Rt5eTJPMeJIodmZ0SJrKysM75PSUmhcePGjtKIiMSXPulteG/6cvLyCvj80030H9jOdaSIoJWIKKApniIibrVqU5/6DasDmux5Oq1ERAFN8RQRccsYQ5/+bXjr74tYvnQ72dknqVGjkutYzqmIiBI9evRwHUFEJK71SWvNW39fRKHP8tmCjVx/c0fXkZzT7QwREZEL0KRZHZo1rwvolEYRFREiIiIXqE9//ymNNau/Yd/eo47TuKciQkRE5AJdHSgiADLnaYOliggREZELVK9Bddpe2hCABSoiVESIiIiURlEb7K2b9rFj+0HHadxSESEiIlIKvfpeQoLHALBgXnxvsFQRISIiUgo1alamU5cmgL/xlLXWcSJ3VESIiIiUUtEGy2+/yWbj+j2O07ijIkJERKSUrurdkiSvB4jvNtgqIkREREqpcuVkruieCsAnGevx+QodJ3JDba+jhKZ4iohElr5pbcjK3MShgyf4cuUuOnaJv9/JWomIApriKSISeS6/MpVKlb0AZMbpKQ2tREQBTfEUEYk83uREevRuyZyP1vDZgo385Bf98Xrj6z+r8fVvG8U0xVNEJPL0SWvDnI/WcPx4LssXb+fKni1cRwor3c4QEREpow6dGlOjZiUA5sfhZE8VESIiImXkSUygd79WAHyRtYWck3mOE4WXiggREZFy6NPfP0sjN7eAzz/b7DhNeKmIEBERKYfW7RpQr341ABbE2S0NFREiIiLlYIwpnuy5fMnXHMk+6ThR+KiIEBERKac+af5ZGj5fIZ9lbnScJnxURIiIiJRT09S6NG1WB4ivWRoqIkRERIKgT7r/lsZXq3axf+9Rx2nCQ0WEiIhIEFwdOOoJkJmxwWGS8FERISIiEgT1G9agTbuGQPyc0lDb6yihKZ4iIpHv6rTWrFuzmy2b9rHz64OkNKntOlJIaSUiCmiKp4hIdOjdtxUJCQaABfNif4OlViKigKZ4iohEh5q1KtOhc2NWLP2azLnrufvHV2GMcR0rZFRERAlN8RQRiQ5909uwYunXfLPrMJs27OWS1vVdRwoZ3c4IEmNMZ2PMVGPMG4FrjjFmiOtcIiISXlf1akmS1wNAZoxvsFQREQTGmAeBecBwa+1D1tqHgNuAEcaYN9ymExGRcKpcJZkrrkwF/Ec9fb5Cx4lCR0VEORljOgNvALdZa7cWPW6tzQb6Aw8GigwREYkTVwfaYB88cJw1q3c5ThM6KiLKbzyQba2de/YTgUJibuA1IiISJ67onkrFSl4A5sdwG2wVEeWXBiwt4fnlQA3tjxARiR/JyUlc1asFAJ8t2Eh+vs9xizjGRwAAEfRJREFUotBQEVEOpxUGW0t42ZbA1/QQxxERkQjSNzBL4/ixUyxbvN1tmBBREVE+lwe+binhNUUFRtcQZxERkQjSsXMTqteoCEDmvNg8paEionxSL+A1hwJfa4QyiIiIRBZPYgK9+/qHcn3+2WZO5eQ5ThR8KiLKp6gwyL6A19YKZRAREYk8V6f5b2nknipgUVZJi9bRSUVE+ZSmMNBKhIhInGnTriEX1a8GwPwYbDylIiIK5OXlkZycXOIlIiKRJyHBcHU/f8+IZV9s5+iRHMeJgktFRPkc+v6XiIhIPCs6peHzFfJZ5kbHaYJLRUT5XMheiLK89gxer5fc3NwSLxERiUxNU+vQuGltADJjbDy4iojgKGm/Q9G+Ca1aiIjEIWMMfQIbLL9cuZP9+445ThQ8KiLKp6gHRPMSXlNUYCwPcRYREYlQV/f374uwFj7JiJ3VCBUR5TMn8LWkUxpFvSSWhDiLiIhEqIaNatCqbQMgtm5pqIgoh8DQrWygcwkvK+pqOS30iUREJFL1CaxGbNqwl107Y+MOt4qI8psEpBpjzte9cggw9/Qx4SIiEn96921FQoIBIDNGJnuqiCgna+1w/HsjvjPu2xgzDP9KxUPhziUiIpGlVp0qtO+UAvgbT1lrHScqPxURwdEF/7jvqcaY1MA1Hn/x0F+rECIiAhSf0vhm52G2bNznOE35qYgIAmtttrU2HXgJSAtcc6y1za21OpUhIiIA9OjdksQkDxAbbbBVRASRtXa5tXZS4JrrOo+IiESWKlUrcPmVzQDIzFhPYWF039JQESEiIhJGRac0Du4/zlerdzlOUz4qIkRERMLoiquaU6FiEhD9jadURESJrKysM64dO3a4jiQiImVQoUISV/bwNzr+bMFGfAWFjhOVnYqIKODz+ejZs+cZ1+TJk13HEhGRMurd139L40h2DqtX7nScpuwSXQeQ7+fxeMjMzDzjsZSUFEdpRESkvLp0a0qlyl5Onsjjk4z1dOraxHWkMjGx0Owilhljcr1er1fjvkVEYsurYz9k3uz/39799cZ1lXsc/z1xYlNKYWLaUkqTqo76B4pOhRPERWgCiS3eQEwvjoTERWPxBhzlqulV5bwDh3ukE0dH595pz2nALSWxQUICCrXTOiqU0+IMLU3rcdOHi712PJ6MZ8Z7/qw9s78fyZp4Mp55PNN0frPWs9b6g+778hf0i//5mfbuHerZY4+MjKhSqVTcfaSd+2E6AwCACI6dSKY0PvrwU/322juRq8mGEAEAQATfOfKovnTfFyRJV155M3I12RAiAACIYN++IR099rgk6fVfvaVK5bPIFe0eIQIAgEiOnXhSkvTxvza0/Ju34xaTASECAIBInvnOQX2ldI8k6dU+nNIgRAAAEMnQ3j06evwJSdIbi29pY2MzckW7Q4gAACCi42FK45NPNnX19euRq9kdQgQAABE9/R+PaP/ovZKkK//bX2dpECIAAIhoaGiPnv1hMqXxm9dW9cmtSuSKWkeIAAAgsnTjqY2Nz/TG66uRq2kdZ2f0icXFxW3fHzhwQAcPHoxUDQCgk7759MO6/4H79MH7H+nKK3/SD04+FbukljAS0Qc4xRMABtuePaZjYUrj2hvX9fHH/XFeEiMRfYBTPAFg8B078ZT+++KSNiu39etfvaWTP3o6dklNESL6xNGjR2OXAADooie++ZAe+vpX9N7f/qkrr7zZFyGC6QwAAHLAzPTsD5M9I5avvq2PPvo0ckXNESIAAMiJ9CyNzz77XK9d+UvkapojRAAAkBOHHn9Q33hkv6T+OB6cEAEAQE6Y2Z3RiN8tv6Ny+VbkihojRAAAkCPHwx4Rn992vfZqvqc0CBEAAOTIo4/dr0cf+6ok6dVX8n2WBiECAICcSbfB/v3vbmj9g39FrmZnhAgAAHIm7Ytwl3756p8jV7MzQgQAADnzyIFRjT3+oKR8r9IgRAAAkEPHw8ZTf/j9u3r/7x9GrqY+QkSfWFxc3Pa1trYWuyQAQBc9G6Y0JOmX/5fPKQ1CRB/gFE8AKJ6vP1zSE089JEm6ktNVGhzA1Qc4xRMAiunYiSf15z+9pzf/+J7e+2tZDz1cil3SNubusWtAA2a2MTw8PLyx0R9nywMAOuf9v3+on0xdkCT9dPpZ/fg/v9eR+x0ZGVGlUqm4+0g798N0BgAAOfXA176sb337YUn5XKVBiAAAIMfSjadW/vL/evfGzcjVbEeIAAAgx77/gydklvw5bw2WhAgAAHLsq/d/Sd9+5hFJ+ZvSIEQAAJBz6ZTG29c/0DtvfxC5mi2ECAAAcu77xx/Xnj3JnEaeRiMIEQAA5Fxp/716ZvygpCRE5GV7BkIEAAB94NiJJ7Vnj+mBB+/TrVuV2OVIYrOp3GOzKQCAJH1yq6KNjU2V9t/b9n11arMptr0GAKAP3PPFYd3zxeHYZWxDiOgTi4uL274/cOCADh48GKkaAACYzsg9M9sYGhoavn379rbrX3jhBZ07dy5OUQCAvsZ0RoFwiicAII8Yicg5GisBAJ3GKZ4AACAqQgQAAMiEENFhZjZvZhOx68CWkZERjYy0NWKHPsFrXRy81vlAiOgAMxszs9NmtiLplKRS7JoAAOg2QkQbzGwiBIe5cNVqzHoAAOgllni2wd0vSzqUfm9mkxHLAQCgpxiJKJC1tTWdO3dOa2trhXrsWIr6fBfxtZbi/d681r3Hc76FEFEgN27c0IsvvqgbN24U6rFjKerzXcTXWor3e/Na9x7P+RZCBAAAyIQQAQAAMmHb6w4ys3klSzyn3P1Sh+6zpRdoeLj58bDurs3NTe3bt09m1nZtuxHzsSuViqTWnqNOKurzXcTXWor3e/NaF+e17uRjp8+fu7f1C7A6Y0BU0v8imhva3Ny83fxmXRHzsdX6U9RRRX2+i/haS/F+b17r3uv357wj6YsQkXPtpkQAALqFnggAAJDJwIeIsB21t/l1M/bvAQBA3gz8dIa7XzCza23eTbkjxQAAMEAGPkRIkrsvx64BAIBBM/DTGQAAoDsIEQAAIBNCBAAAyIQQ0Vlj4fK7UasAAKAH2Pa6DWY2JmleUknSaLhMlSWtSlp398kI5QEA0FWECABAXzCzcUlnJa2Hq8YkzXXqrCLsHiECAJB7ZnZa0qykw+6+Gq4rSXpZ0jV3n45ZX1ERIgAAuRZGIJYkTbr75Zq/K0m6KWna3S/EqK/ICBFADTOblbTC/5CAfDCzBUlH3H1/lr9H97A6A6gSPvHMaHuTLAZAOEdnycxWwte8mZ2KXRdaMiGp0fEFy5JKvJ69R4gAtpuPXQA6z8zmJB1y98PufkjSYSXNefNmthCGxJFDVcFgtcHNVsIlK+F6jBABBGY2o+QTDQZIeF3L7n4mvc7dy6ER74KST7k/j1Ufmkr33VlpcJs0YBzpci2oQYgAdGfPj0OSFmLXgo6brg4QNdLrT5nZRK8Kwq6MNb/JnSWfjCj1GCECSMxq6w0FAyIMhe+4h4C7l7U1+sRQeD6lwaDcwm1Hu1kI7kaIQOGF9ef/Fd5QMFjGJM2YWaNel3QofLwH9WD3dhMMGInoMUIECi001E2y493AO9WgeXI3n3QBVNkbuwAgsllJz8cuAl1zQUlj3tUGI01pM16j7n/Es978JoiFEIHCCo10S0xjDK7w2k7t9PdhX5B0JOKlnhSF3drNv0/+LfcYIQJFNu3uO77BoBDOhsvzhMnca9TvkPZNMGrRY/REoJDC1tasxiiwMApxStLlBktAEV86zXSowW3SgME+Lz1GiEBuhW2Kvc2vm3Xud1yS0pMAUVjzkpbdnaWd+Zbu3dJolUa6l8TVLteCGkxnILfc/YKZNdovvxX1hqjPMo2RL2GZ7Vybd1Nu9QCmcGBTWdLJNh8TXebul82srMZLcNNdLVll1WOc4olCCaMQL2vnudNRJUOj5arbnGEJaPelI0RtKLcyuhTO0Tgi6SR9EP0hTD/OKDn/5K7X2MxcybQUo0o9RogAqoQ3mNNKgsP52PWgs8KIx7S7H45dC3bHzFaUTD9N1Vw/o6RB9jBTlL1HTwSAQghLencMEOHTLvLrsJLjvufNbCx8zUqaVjKqRICIgJ4IAAMvTJWcaRAgSmLL5FwLU0+T4bVMD0tbYGVNXIQIYLu0y7vRcjL0kRAQZlWz6VTVNtijSj7NLvW4NGTg7stiKWdu0BOBwgvD3LNKAkT1p9FVJQ2Wz4f/caEPhbn0Vo6TnnT3y92uBxgkhAgAA6uqq78V+1mtAewOIQIAAGTC6gwAAJAJIQIAAGRCiAAAAJkQIgAAQCaECABtM7NWllACGDCECACZmNnN9Mh1SSux64nJzFZqjqAnVKEQCBEAMgnHbnOQlSR3P6SqXU45xwFFQYgAkFnYyZM3zES62ym7XqIwCBEA0BnpoVDzUasAeogQAQCdMRkuGYlAYbDtNYC2pAdcubvFriWm0GBaDr0iQCEwEgEAbTKz8fBHRiFQKHtjFwBg8JnZKUnPVV1VkjTv7hca/ExJyRHto5LWlRznPS/poqSz4WYTkuYa3U+PpP0QC9VXhmPmp8K3o0qOleekUAwMQgSArjKzOUmj7j5Vc/2MmS1JOln7xhoCxHVJF9Ofq7puVtJj7l42sxlVLa2M6K5+CDM7LUnuPh3qPC3pqqTzvS8P6A56IgC0pVFPhJnNSxoP+yjU+9m6f29mC0o+3e+vDhjhjXlO0pS7X+rgr9GW2n6IUOd6WmMIS+OSDodlscBAoCcCQFeEofxTSkYOdvKSpLHwSb3ahKTVOkP/18LldGeqbF9tP0QIEKs1IeekkkBEgMBAIUQA6JYz4fLaTjeoelO9EzSq3pTr9Q6k1422XV3n3OmHMLNZJQFiW4Olu5fphcAgIkQA6JYj4bLZjpZlaVt4SG9fqnPbUs1t8iDth5iT9A81CE3AoCFEAOi40ARZLwQ0Mioln9olXVIyzVF7H+mn/pdarGOpUY1mNhe+5sNXloOzJpSEmvQcketmtlCndmDgECIAdFzN0H2rUw/rVX9+Plze2UI6jFSclXSmWW+BmY2F5szxBm/mL0tacPfpsALkjKSl3QSJ6n4Id1929/OSHlMSLHYMMMCgYIkngG5ZVbK3Q7NP5CVpW3+ElCyHnJRUCis41pVMe5xsFCBCYJgPj71jD0JofixVNz+6+6qZXVQyLTG508/WuGt/iLD09LKkCTMrpYHKzE7laUUJ0AmECADdMqvkDfk5SXXf+MMmVJJUu1nUZPhULyVTGy0Jb9iT4b5nlKwOqWdK9fsqlpQEmFalG2jV7lRZPTVTfVtCBAYK0xkAuiLsInlZjd+Up5WsZrhryWZVwOiGCdUfqVgNjz1R5+/qGVf9paglVYWU8LvMZagTyDVGIgC0a1RKphLqvJlOSZoP/QlTNRtHzYWfrTd1sBB+7pK2eiXSn11RspNlpiWTVT0P6w1u1rQpsipo1Dsv47KkH1d9/1ztjp3AICBEAMgk7FQ5qq033JtmtqpkKmJV2ppeCJ/Ef262bVPLq/VGIMLPnTezae08HTFnZstKgknW5Z6NQkirzaBl1R9hOCNpNPRzpN8DA4cQASCTnbay3uG2l9RiP0BY8TCvpKfiYs3oRUlJs+ZzkmaUjFh04+yMpiMRYUOpusd+h5oZecDAoycCQN68rHAyZ+2URdj5cdndzyiZBhmrWmbZqnQao1FQYHdJoAWECAC5EfoVStphNUeNdGfIXW0QVRVM6k1Z5HFHTCC3CBEAciP0N6yqtamAtHGxXmNjM6uqPxKRLs3Mcp9A4RAiAOTNpJKNmupuQx22q55R0jNxOOMqjTltne1R+9gECKBFNFYCyJUwGnEoLKGcDSs60j6GdApiwd3rNjW2+BjnzWzazE6H/SzSqZQJbZ2BAaAJc/fYNQBAx4T9J8aUjDSkmz4tKwkeF6puV1LVEeRKAspLzc7lALCFEAEAADKhJwIAAGRCiAAAAJkQIgAAQCaECAAAkAkhAgAAZEKIAAAAmRAiAABAJoQIAACQCSECAABkQogAAACZECIAAEAmhAgAAJAJIQIAAGRCiAAAAJkQIgAAQCaECAAAkAkhAgAAZEKIAAAAmRAiAABAJoQIAACQCSECAABkQogAAACZECIAAEAm/wbBnabkJr2oAQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 525x525 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from codes.plot_utils import plot_line_points\n",
    "plot_pretty(fontsize=12)\n",
    "\n",
    "plot_line_points(k, Pk, xlabel=r'$k$', ylabel=r'$P(k)$', figsize=3)\n",
    "plot_line_points(lk, lPk, xlabel=r'$\\log_{10} k$', ylabel=r'$\\log_{10} P(k)$', figsize=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots above show that $P(k)$ is changing very rapidly for small values of $k$ and slowly for large values. Functions that change rapidly are more difficult to approximate accurately because more training points are needed in the areas of rapid change. So it is often useful to examine how $\\log_{10}$ of the function behaves as a function of $\\log_{10}$ of the parameter (i.e., examine function in logarithmic space). This is what is done in the second plot above. We can see that in this case, $\\log_{10} P(k)$ changes relatively slowly and so approximating $\\log_{10} P(k)$ as a function of $\\log_{10} k$ is preferrable. \n",
    "\n",
    "Overall, $\\log_{10} P(k)$ looks like two power law functions of different exponent smoothly joined in the middle, but notice little \"wiggles\" around $\\log_{10} k\\sim -1$. These are real features of this function and make it not entirely trivial to approximate with, say, a polynomial of low order. \n",
    "\n",
    "Next, let's see how <tt>k</tt> and <tt>lk</tt> are spaced. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.04304196 0.04303239 0.04305963 0.04303769 0.04302909 0.04304367\n",
      " 0.04303872 0.04302599 0.04304582 0.04303598]\n",
      "0.0003967464299742818 0.04306143598772039\n"
     ]
    }
   ],
   "source": [
    "# np.diff computes difference of 2 subsequent elements of lk\n",
    "dlk = np.diff(lk) \n",
    "# if we look at the first 10 elements, it may seem that l (and lk) are evenly spaced\n",
    "print(dlk[:10])\n",
    "# however, this is not the case\n",
    "print(dlk.min(), dlk.max())\n",
    "# spacing of data points is not uniform overall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polynomial interpolation\n",
    "\n",
    "In class on Friday I showed that straitforward solution of the linear system of equations for the polynomial coefficients by computing the inverse of the Vandermonde matrix as in the code below (reproduced from the <tt>03_approximation.ipynb</tt> notebook) is numerically unstable, and fails in interpolation of $d_L$ using as few as 20 training points. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Om0tr are values of Om0 at ntr training points\n",
    "# V is the Vandermonde matrix (see 03 notebook for the definition)\n",
    "ntr = np.size(Om0tr)\n",
    "V = np.vander(Om0tr, N=ntr, increasing=True) # set up the Vandermonde matrix\n",
    "\n",
    "# dLtr are tabulated values of d_L at Om0tr points\n",
    "a = np.dot(np.linalg.inv(V), dLtr) \n",
    "\n",
    "# np.poly1d initializes a 1d polynomial class with input coefficients\n",
    "# [::-1] reverses order of the coefficients, as required by this routine\n",
    "pfit = np.poly1d(a[::-1]) \n",
    "# now compute polynomial approximation at test Om0 values in vector Om0test\n",
    "dLpoly = pfit(Om0test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if we instead using <tt>numpy.linalg.solve</tt> routine to solve the same system of equations, namely if we replace the following line in the <tt>03_approximation.ipynb</tt>\n",
    "\n",
    "    a = np.dot(np.linalg.inv(V), dLtr) # a = V^{-1} f\n",
    "\n",
    "with \n",
    "\n",
    "    a = np.linalg.solve(V, dLtr)\n",
    "    \n",
    "this fixes the problem at least until the number of training of points is much larger than 20. \n",
    "This is because <tt>numpy.linalg.solve</tt> uses a more numerically stable algorithm with LU decomposition of the Vandermonde matrix to solve the system of equations (the routine uses fortran-based LAPACK library, see <a href=\"https://software.intel.com/en-us/node/468876\">here</a>). \n",
    "\n",
    "Likewise, if we replace that line with \n",
    "\n",
    "    a = np.linalg.lstsq(V, dLtr, rcond=None)[0] \n",
    "   \n",
    "this solves the problem as well because this routine uses numerically robust Singular Value Decomposition (SVD) method to solve the system of normal equations. Note that we can control the order of approximating polynomial in this case, by limiting the size of the Vandermonde matrix to the required order <tt>poly_order</tt> as shown here: \n",
    "    V = np.vander(xtr, N=poly_order, increasing=True)\n",
    "\n",
    "Let's see if this will also work for approximating function $P(k)$ that was read above. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Task 1a. (10 points)** Based on the code in the 03 notebook reproduced above, write a function that constructs interpolating (using <tt>np.linalg.solve</tt>) or approximating (using <tt>np.linalg.lstsq</tt>) polynomial for input vector of test points and corresponding tabulated function values with the following interface:\n",
    "\n",
    "    def poly_approx(xtr, ftr, xt, method='interpolation', poly_order=1):\n",
    "        \"\"\"\n",
    "        function constructing an approximating polynomial for the input \n",
    "        vector of training values of the approximated function ftr, \n",
    "        tabulated at the vector of training points xtr and compute the values\n",
    "        of constructed polynomial at the input vector of test values xt\n",
    "    \n",
    "        Parameters:\n",
    "        -----------\n",
    "        xtr, ftr: array like\n",
    "                    vectors of training points and tabulated function values\n",
    "        xt: array like\n",
    "                    vector of test points\n",
    "        method: string\n",
    "                defining the method to be used in construction of the polynomial\n",
    "                'interpolation' - construct interpolating polynomial using np.linalg.solve\n",
    "                'regression' - construct approximating polynomial np.linalg.lstsq\n",
    "        poly_order: integer\n",
    "                    order of the approximating polynomial. This parameter is ignored\n",
    "                    for method='interpolation', but is used for 'regression'\n",
    "                    must be 0<=poly_order<=np.size(xtr)\n",
    "        Returns:\n",
    "        --------\n",
    "        ft: array like\n",
    "            values of the polynomial at xt \n",
    "        a:  array like\n",
    "            polynomial coefficients\n",
    "        \"\"\"\n",
    "\n",
    "Test this function using $\\sin(x)$ as a  test function as shown below and demonstrate that your routine works (i.e., produces sensible results). Try changing the number of training points used <tt>ntr</tt> from 20 to smaller values. Identify and report the smallest <tt>ntr</tt> for which you get *interpolating* polynomial that still looks qualitatively as a sine (for example, it has a clear minimum and maximum in the $[0,2\\pi]$ interval). Try to come with your interpretation of the meaning of the minimum value that you find and briefly describe it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, xmax, ntr = 0., 2.*np.pi, 20\n",
    "xtr = np.linspace(xmin, xmax, ntr)\n",
    "ftr_int = np.sin(xtr)\n",
    "ftr_fit = np.sin(xtr) + np.random.normal(scale=0.05, size=np.size(xtr))\n",
    "\n",
    "nt = 1000 # should be different from nt\n",
    "xt = np.linspace(xmin, xmax, nt)\n",
    "\n",
    "poly_int, a_int = poly_approx(xtr, ftr_int, xt, method='interpolation')\n",
    "\n",
    "poly_order=17\n",
    "poly_fit, a_fit = poly_approx(xtr, ftr_fit, xt, method='regression', poly_order=poly_order)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from codes.plot_utils import plot_line_points\n",
    "\n",
    "plot_line_points(xt, poly_int, xp=xtr, yp=ftr_int, points=True,\n",
    "                 xlabel=r'$x$', ylabel=r'$f(x)$', figsize=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_line_points(xt, poly_fit, xp=xtr, yp=ftr_fit, points=True,\n",
    "                 xlabel=r'$x$', ylabel=r'$f(x)$', figsize=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1b. (10 points)** Now apply the routine you implemented to interpolate $P(k)$ read in below. If you don't manage to get the working implementation in 1a, you can use <a href=\"https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html\"><tt>np.polyfit</tt></a> routine for this exercise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nkt = 1000\n",
    "lkt = np.linspace(lk.min(), lk.max(), nkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_int, a_int = poly_approx(lk, lPk, lkt, method='interpolation', poly_order=np.size(lkt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_line_points(lkt, poly_int, xlabel=r'$x$', ylabel=r'$f(x)$', figsize=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens? Try doing this for $\\log_{10} k$ and $\\log_{10} P(k)$ instead of $P(k)$. Then try this using least squares regression option. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_int, a_int = poly_approx(lk, lPk, lkt, method='regression', poly_order=np.size(lk))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_line_points(lkt, poly_int, xlabel=r'$x$', ylabel=r'$f(x)$', figsize=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do you get a sensible result? Think about what the routine is doing and try to come up with an explanation for the outcome of this exercise. (7 points)\n",
    "\n",
    "Try regression with smaller poly_order values and find polynomial order for which at least somewhat sensible approximation can be obtained. What is that order? Think about a possible explanation for why this works for small order (3 points). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SciPy routines for spline interpolation and approximation\n",
    "\n",
    "Spline interpolation is implemented in the SciPy library <a href=\"http://scipy.github.io/devdocs/interpolate.html\"><tt>scipy.interpolate</tt></a> in routines <a href=\"https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.interpolate.interp1d.html\"><tt>interp1d</tt></a> and <a href=\"https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html\"><tt>UnivariateSpline</tt></a>. Scipy   has more spline routines for more specialized tasks (see the list <a href=\"http://scipy.github.io/devdocs/interpolate.html\">here</a> and examples of their use <a href=\"https://docs.scipy.org/doc/scipy-1.0.0/reference/tutorial/interpolate.html\">here</a>), but the routines listed here are most frequently used for 1-dimensional function interpolation/approximation. \n",
    "\n",
    "<a href=\"https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.interpolate.interp1d.html\"><tt>interp1d</tt></a> is a routine written in python, although parts of it are written in <a href=\"https://cython.org/\">Cython</a> and use a couple of <a href=\"\">Fitpack library routines</a> and can perform only spline interpolation (i.e. approximation matching the function values at the training values to machine precision). \n",
    "\n",
    "<a href=\"https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html\"><tt>UnivariateSpline</tt></a> is a class-based python wrapper around routines <a href=\"https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splev.html\"><tt>splev</tt></a>, <a href=\"https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splrep.html#scipy.interpolate.splrep\"><tt>splrep</tt></a>, and <a href=\"https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splprep.html#scipy.interpolate.splprep\"><tt>splprep</tt></a> that use old, efficient fortran-based Fitpack library of interpolation routines. This function allows construction of the spline that is not interpolating data, but approximating it using the least squares regression discussed in class on Friday. The effective weights of tabulated function data points are set by the input parameter <tt>s</tt>.   \n",
    "\n",
    "For <tt>s=0</tt> <tt>UnivariateSpline</tt> is equivalent to routine\n",
    "SciPy <a href=\"http://scipy.github.io/devdocs/generated/scipy.interpolate.InterpolatedUnivariateSpline.html#scipy.interpolate.InterpolatedUnivariateSpline\"><tt>InterpolatedUnivariateSpline</tt></a> that constructs spline that matches the function values at the training points to machine precision.\n",
    "\n",
    "For <tt>s>0</tt> <tt>UnivariateSpline</tt> constructs spline that satisfies condition for the $N+1$ training points $x_i$:\n",
    "\n",
    "$$\\sum\\limits_{i=0}^N w_i [f_i-{\\rm spline}(x_i)]^2 \\leq s$$\n",
    "\n",
    "by solving the system of normal equations for the spline coefficients. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are my routines that compute coefficients of the interpolating \n",
    "natural cubic spline and the spline itself, given the vector of these coefficients.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cubic_spline_coefficients(xi, fi):\n",
    "    \"\"\"\n",
    "    compute coefficients of the interpolating natural cubic spline\n",
    "    see Appendix in the note for the derivation and details of the algorithm\n",
    "    \n",
    "    Parameters: \n",
    "        xi, fi: numpy float vectors\n",
    "                tabulated points and function values\n",
    "            \n",
    "    Returns:\n",
    "        a, b, c, d: numpy float vectors\n",
    "                cubic spline coefficients \n",
    "                dx = x - xi[i]\n",
    "                fx = a[i] + dx*(b[i] + c[i]*dx + d[i]*dx*dx)\n",
    "\n",
    "    \"\"\"\n",
    "    delx = np.diff(xi); delf = np.diff(fi)\n",
    "    # form matrices to solve for spline coefficients\n",
    "    vx = np.zeros_like(xi)\n",
    "    # form rhs vector using python's array slicing \n",
    "    vx[1:-1:] = 3.*(delf[1::]/delx[1::] - delf[:-1:]/delx[:-1:])\n",
    "    # construct 3 diagonals\n",
    "    nx = np.size(xi)\n",
    "    diags = np.zeros((3, nx))\n",
    "    diags[1,0] = 1.; diags[1,-1] = 1.\n",
    "    diags[1,1:-1:] = 2. * (delx[1::] + delx[:-1:])\n",
    "    diags[0,1:] = delx[:]\n",
    "    diags[2,1:-1] = delx[1:]\n",
    "    # solve for coefficients c using Thomas algorithm for tri-diagonal matrices\n",
    "    # see https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)\n",
    "    ac, bc, cc, dc = map(np.array, (diags[0,:], diags[1,:], diags[2,:], vx)) # copy arrays\n",
    "    for k in range(1, nx):\n",
    "        mk = ac[k] / bc[k-1]\n",
    "        bc[k] = bc[k] - mk * cc[k-1] \n",
    "        dc[k] = dc[k] - mk * dc[k-1]\n",
    "    \n",
    "    c = np.zeros_like(bc)\n",
    "    c[-1] = dc[-1] / bc[-1]\n",
    "\n",
    "    for k in range(nx-2, -1, -1):\n",
    "        c[k] = (dc[k]-cc[k]*c[k+1])/bc[k]\n",
    "\n",
    "    # now get the rest of the coefficients\n",
    "    b = delf[::]/delx[::] - (c[1::] + 2.*c[:-1:])*delx[::]/3.\n",
    "    d = (c[1::] - c[:-1:])/(3.*delx[::]) \n",
    "    a = fi\n",
    "    return a, b, c, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cubic_spline(xt, xi, a, b, c, d):\n",
    "    \"\"\"\n",
    "    piecewise linear approximation of f(x) given input of tabulated values of xi and fi\n",
    "    note that xi are expected in ascending order\n",
    "    \n",
    "    Returns:\n",
    "        vector of spline values at test points x \n",
    "\n",
    "    \"\"\"\n",
    "    n = np.size(xi) - 1\n",
    "    fxt = np.empty_like(xt)\n",
    "    for j, x in enumerate(xt):\n",
    "        for i in range(n):\n",
    "            if (x >= xi[i]) and (x <= xi[i+1]):\n",
    "                # reusing computations is always a good idea, but here we also can return dfdx\n",
    "                dx = x - xi[i]\n",
    "                fxt[j] = a[i] + dx*(b[i] + c[i]*dx + d[i]*dx*dx)\n",
    "        \n",
    "    return fxt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the test that the spline interpolation produces good results for interpolation of <tt>Plk</tt> as a function of <tt>lk</tt>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b, c, d = cubic_spline_coefficients(lk, lPk)\n",
    "\n",
    "spline_int = cubic_spline(lkt, lk, a, b, c, d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_line_points(lkt, spline_int, \n",
    "                 xlabel=r'$\\log_{10} k$', ylabel=r'$\\log_{10} P(k)$', figsize=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that once we have the spline coefficients for a given interval $a_i$, $b_i$, $c_i$, and $d_i$, we can compute derivatives of the interpolating spline as well as its integral. For interpolating spline \n",
    "\n",
    "$$s_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3,$$\n",
    "\n",
    "we have for derivatives\n",
    "\n",
    "\\begin{eqnarray}\n",
    "f^\\prime(x\\in[x_i,x_{i+1}])\\approx s^\\prime_i(x)&=&b_i + 2c_i(x-x_i)+3d_i(x-x_i)^2,\\\\\n",
    "f^{\\prime\\prime}(x\\in[x_i,x_{i+1}])\\approx s^{\\prime\\prime}_i(x)&=&2c_i+6d_i(x-x_i),\\\\\n",
    "f^{\\prime\\prime\\prime}(x\\in[x_i,x_{i+1}])\\approx s^{\\prime\\prime\\prime}_i(x)&=&6d_i\n",
    "\\end{eqnarray}\n",
    "\n",
    "It is also straightforward to use interpolating spline to estimate the integral of the approximated function. For example, the integral in the interval $[x_i,x_{i+1}]$ (using change of variables $x\\rightarrow\\xi=x-x_i$:\n",
    "\n",
    "$$\\int\\limits_{x_i}^{x_{i+1}}f(x)dx\\approx \\int\\limits_{x_i}^{x_{i+1}}s_i(x)dx = \\int\\limits_{0}^{x_{i+1}-x_i}s_i(\\xi)d\\xi=a_i(x_{i+1}-x_i)+\\frac{1}{2}b_i(x_{i+1}-x_i)^2+\\frac{1}{3}c_i(x_{i+1}-x_i)^3+\\frac{1}{4}d_i(x_{i+1}-x_i)^4$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1c. (10 points)** Add calculation of the first, second, and third derivatives based on the spline interpolation to the routine <tt>cubic_spline</tt> above (3 points). Write a routine to estimate integral of $f(x)$ in the interval $[a,b]$ using interpolating cubic spline for $f(x)$ and its coefficient vectors $a$, $b$, $c$, $d$. (4 points).\n",
    "\n",
    "Plot the first three derivatives of $\\log_{10} P$ for some function (for example, the sine) and comment on their behavior (smoothness, continuity). Test the accuracy of the integral for a function for which you know the exact answer (e.g., $\\sin x$ or $e^x$.) (3 points). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1d. (5 points)** Suppose we want to compute the *logarithmic slope* of $P(k)$:\n",
    "\n",
    "$$\\frac{d\\log_{10}P}{d\\log_{10} k}=\\frac{d\\ln P}{d\\ln k}=\\frac{k}{P}\\,\\frac{dP}{dk}.$$\n",
    "\n",
    "Compute this derivative using your routine implemented in 1c using <tt>lk</tt> and <tt>lPk</tt> vectors and plot this derivative as a function of <tt>lk</tt>. Compare your plot to the left panel of Figure 4 in <a href=\"https://arxiv.org/pdf/1407.4730\">this paper</a>. It should look roughly similar, but with some differences. Describe the differences. Briefly describe what causes them. \n",
    "\n",
    "Note that if you don't complete 1c, you can construct spline interpolation of <tt>lk</tt> and <tt>lPk</tt> using SciPy <tt>UnivariateSpline</tt> with <tt>s=0</tt> and use its <tt>derivative()</tt> method to compute the derivative (see Task 1e below). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1e (5 points)** Use SciPy <a href=\"https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html\"><tt>UnivariateSpline</tt></a> routine to approximate <tt>lk</tt> and <tt>lPk</tt> and this class'  <tt>derivative()</tt> method to construct approximation for the derivative, as illustrated below, to compute the derivative $$d\\log_{10}P/d\\log_{10} k$$. Compute and plot the derivative for <tt>s=0.</tt>, <tt>s=0.0001</tt>, <tt>s=0.001</tt>, and <tt>s=0.01</tt> and comment on how the derivative changes for different values of the smoothing parameter <tt>s</tt>. Try to explain the behavior you observe. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example of use of derivative method of UnivariateSpline \n",
    "from scipy.interpolate import UnivariateSpline\n",
    "\n",
    "ntr = 4\n",
    "xtr = np.linspace(0.,2.*np.pi, ntr)\n",
    "spl = UnivariateSpline(xtr, np.sin(xtr), s=0.) # spline approximating sine\n",
    "dspl = spl.derivative() # spline approximating dsin(x)/dx=cos(x)\n",
    "\n",
    "nt = 100\n",
    "xt = np.linspace(0.,2.*np.pi, nt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_line_points(xt, dspl(xt), xlabel=r'$x$', ylabel=r'$f^\\prime(x)$', \n",
    "                 figsize=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spl = UnivariateSpline(lk, lPk, s=0.001)\n",
    "dspl = spl.derivative()\n",
    "\n",
    "plot_line_points(lkt, dspl(lkt), xlabel=r'$\\log_{10} k$', ylabel=r'$\\log_{10}P(k)$', \n",
    "                 figsize=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1f. (5 points)** Using spline to invert function. Using spline we can invert a function or, for $y=f(x)$ compute $x$ for a given $y$. For this we need to generate vectors of $x$ and $y$ and construct interpolating spline for $y$ as the argument and $x$ as a function value. The spline can then be used for input $y$ values to estimate corresponding $x$. This, of course, assumes that function is monotonic and that a unique value of $x$ can be obtained for a given $y$. \n",
    "\n",
    "Construct interpolating spline for $y=\\sin(x)$ on the interval $[0,\\pi/2]$ and use it to compute $x$ for $y=1$, $y=0$, $y=0.5$. Compare the values you get to the output of <tt>np.arcsin</tt> for these values of $y$. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}