{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing Function Inner Products, Norms & Metrics\n",
    "\n",
    "**Randall Romero Aguilar, PhD**\n",
    "\n",
    "This demo is based on the original Matlab demo accompanying the  <a href=\"https://mitpress.mit.edu/books/applied-computational-economics-and-finance\">Computational Economics and Finance</a> 2001 textbook by Mario Miranda and Paul Fackler.\n",
    "\n",
    "Original (Matlab) CompEcon file: **demmath02.m**\n",
    "\n",
    "Running this file requires the Python version of CompEcon. This can be installed with pip by running\n",
    "\n",
    "    !pip install compecon --upgrade\n",
    "\n",
    "<i>Last updated: 2022-Oct-23</i>\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('seaborn-dark')\n",
    "import scipy.integrate as integrate\n",
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Class function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the class **function**. An object of class **function** operates just as a lambda function, but it supports several function operations: sum, substraction, multiplication, division, power, absolute value, integral,  norm, and angle.\n",
    "\n",
    "This example illustrates how it is possible to overwrite the methods of the function class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class function:\n",
    "    def __init__(self, func):\n",
    "        self.f = func     \n",
    "    \n",
    "    def __call__(self, *args):\n",
    "        return self.f(*args)\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        return function(lambda *args: self.f(*args) + other.f(*args))\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        return function(lambda *args: self.f(*args) - other.f(*args))\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        return function(lambda *args: self.f(*args) * other.f(*args))\n",
    "   \n",
    "    def __pow__(self, n):\n",
    "        return function(lambda *args: self.f(*args) ** n)\n",
    "   \n",
    "    def __truediv__(self, other):\n",
    "        return function(lambda *args: self.f(*args) / other.f(*args))\n",
    "   \n",
    "    def integral(self, l, h):\n",
    "        return integrate.quad(self.f, l, h)[0]\n",
    "    \n",
    "    def abs(self):\n",
    "        return function(lambda *args: np.abs(self.f(*args)))\n",
    "    \n",
    "    def norm(self, l, h, p=2):\n",
    "        return  ((self.abs()) ** p).integral(l, h) ** (1/p)\n",
    "    \n",
    "    def angle(self, other, l, h):\n",
    "        fg = (self * other).integral(l, u)\n",
    "        ff = (self**2).integral(l, u)\n",
    "        gg = (other**2).integral(l, u)\n",
    "        return np.arccos(fg*np.sqrt(ff*gg))*180/np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute inner product and angle\n",
    "\n",
    "Define the functions $f(x) = 2x^2-1$ and $g(x)= 4x^3-3x$, both over the domain $[-1,1]$. Compute their inner product and angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "∫(f*g)(x)dx = 0.00\n",
      "∫(f^2)(x)dx = 0.93\n",
      "∫(g^2)(x)dx = 0.97\n",
      "Angle in degrees = 90°\n"
     ]
    }
   ],
   "source": [
    "l, u = -1, 1\n",
    "f = function(lambda x: 2 * x**2 - 1)\n",
    "g = function(lambda x: 4 * x**3 - 3*x)\n",
    "\n",
    "fg = (f*g).integral(l, u)\n",
    "ff = (f**2).integral(l, u)\n",
    "gg = (g**2).integral(l, u)\n",
    "angle = f.angle(g, l, u)\n",
    "\n",
    "print(f'∫(f*g)(x)dx = {fg:.2f}')\n",
    "print(f'∫(f^2)(x)dx = {ff:.2f}')\n",
    "print(f'∫(g^2)(x)dx = {gg:.2f}')\n",
    "print(f'Angle in degrees = {angle:.0f}°')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Function Norm\n",
    "\n",
    "Now compute the norm of function $f(x) = x^2 - 1$ over the domain $[0, 2]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "∥f∥₁ = 2.000\n",
      "∥f∥₂ = 1.751\n"
     ]
    }
   ],
   "source": [
    "l, u = 0, 2\n",
    "f = function(lambda x: x ** 2 - 1)\n",
    "\n",
    "print(f'∥f∥₁ = {f.norm(l, u, 1):.3f}')\n",
    "print(f'∥f∥₂ = {f.norm(l, u, 2):.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute function metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "∥f-g∥₁ = 0.883\n",
      "∥f-g∥₂ = 1.000\n"
     ]
    }
   ],
   "source": [
    "l, u = 0, 1\n",
    "\n",
    "f = function(lambda x: 5 + 5*x**2)\n",
    "g = function(lambda x: 4 + 10*x - 5*x**2)\n",
    "\n",
    "print(f'∥f-g∥₁ = {(f-g).norm(l, u, 1):.3f}')\n",
    "print(f'∥f-g∥₂ = {(f-g).norm(l, u, 2):.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Illustrate Function metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x23c63c32190>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEECAYAAAAlEzNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6JUlEQVR4nO3deXxU5bkH8N+s2WaSyb7vCwECZGEVAigCAtYF1AAl0Ootbb2K9WLrSmpRFCu1rXC9rdhSRZQI0kpQQHZkJyEJBMgCSSb7PpNkZjIzmZlz/wiJRiCZJHPmzPJ8Px8/mjnJnMcDOc+c933e5+UxDMOAEEKI0+JzHQAhhBBuUSIghBAnR4mAEEKcHCUCQghxcpQICCHEyVEiIIQQJ8daIigsLERmZuZtr1++fBnLly/HsmXLsGbNGuh0OrZCIIQQYgYhG2+6detW7N27F25ubv1eZxgG69atw/vvv4/IyEjs2rULtbW1iImJYSMMQgghZmDliSAiIgKbN2++7fWKigrIZDJ8/PHHWLFiBZRKJSUBQgjhGCtPBPPnz0dNTc1trysUCuTn52PdunWIjIzEr371KyQlJWHatGm3fW9zcycboRFCiEPz95cO+WesOlksk8kQGRmJuLg4iEQipKeno6ioyJohEEII+RGrJoLw8HCo1WrI5XIAQG5uLuLj460ZAiGEkB9hZWjox3JycqDRaJCRkYENGzZg7dq1YBgGKSkpmD17tjVCIIQQchc8W+0+SnMEhBAydDY/R0AIIcT2UCIghBAnR4mAEEKcHCUCQghxcpQICCHEyVEiIIQQJ0eJgBBCnJxVFpQRYq+6jSa0qPVoVeuhM5hgMDEwmBiI+DxIXYWQiIWQuYkgdaVfJWK/6G8vIQCMJgZlzSpca1ShrEmFsmY1qhRdUHR1m/XzXq5ChHu7IVzmhsRACcaHeGJUgAQiAT10E9tHK4uJ02pV63HiZisuyBXIrVKiXWsAAHiIBYj390CUjzsCpC7w9xDD10MMVxEfQj4fQj4PeqMJKp0RKp0BbRo9atu1qFJ0Qd6mQZNKDwAQC3gYH+KJ9FhfzIz1RZjMbaBwCLGI4awspkRAnIpKZ8DR0hYcLG5CbrUSJgYIkIgxJdIbkyO9MS5EihBPV/B4vGGfo1mlw5W6DhTWdeC8XIGbLRoAQKyfOxaNCcSisYHwcRdb6n+JkH4oERByF5WtGnxRUId9VxvQ1W1CmMwV8xIDMHeUP2J93Ud04x9MjbIL35W34VBxM67Ud0DA52FmrC8yUkKQGubF6rmJ86FEQMiPXKnrwNazcpytVEAk4GFeYgAemxCMsUFSTm7AFa0afHWlAd9ca4SiqxtJwVKsmhSOmXG+4FNCIBZAiYCQW641dOLDM3KcrmiDt5sIGakheHR8sM0MyWi7jfj6WiO2X6xBbbsWsX7ueDY9BvdEe9MTAhkRSgTE6TV16vD+yXIcLG6Gl6sQKyaG4YmUULiLBVyHdkcGE4MjJc34+5lKVCu1mBghw3Mzo5EYOPRfZkIASgTEiekNJnx+qRb/OCeH0cRgxaRwZE4Mg8TFPiqku40m7Cmsx9azcrRrDVgyIRjPpEfbTfzEdlAiIE7pakMn/nCgBBWtGsyM9cXzs2PstlRTpTPgwzNyZOfXwsddjN/OicO9cb40XETMRomAOBW9wYR/nJPj4wvV8PUQ45W5CZge48N1WBZxraETG74tRWmzGnMS/PDy/fHwchNxHRaxA5QIiNOoaNXg1a+vo6xZjQfHBuJ/Zsc6XJsHg4nBjtwa/O10JXzcRfjDgkRMjJBxHRaxcTa1VWVhYSEyMzPvenzdunXYtGkTW6cnDuyba41YteMSmlV6/OmRsfj9A6McLgkAgJDPw6rJ4di2PBmuIgGe3nUZW76rgMFkk5/diB1jJRFs3boVr732GnQ63R2P79y5E6WlpWycmjgwbbcRG74txe/3lyAxQIIdmamYGevLdVisSwyU4tPMVDw8LggfX6jGmi+vQKkxrwcSIeZgJRFERERg8+bNdzyWn5+PwsJCZGRksHFq4qBaVDr88ovL+M+VBvxscjg+eGICAqQuXIdlNW4iAV6dl4Cs+QkorG3Hyh2XUNxIw6fEMlhJBPPnz4dQePujelNTE7Zs2YKsrCw2TkscVEmjCqt25KOiVY13HxqD/06PhpDvnFU0P0kKwtalyTAxwH/tLMShkmauQyIOwKoDqwcOHIBCocDq1avR3NwMrVaLmJgYLF682JphEDtyvKwF674phpebCFuXJmNUgITrkDg3JkiKT1ak4HdfXcMr+66joUOLFRPDqMSUDJtVE8HKlSuxcuVKAMCePXtQXl5OSYDc1Z7L9dh4qAxjg6V49+Gx8POwjfYQtsDHXYz/fXw8Xt9fjPdPVqChQ4f/uTcWAid9UiIjY5VEkJOTA41GQ/MCxCwMw2Db+Wr83+lKTI/2wcafjIaryDZbRHDJRcjHhgdHI/BEBXbk1aBZrceGRYm0GQ4ZMlpHQGyKiWHw5+Pl2HmpFgtGByBrfgKEdGMb1Gd5Nfjz8XLcE+2Nd34yhhKnE7OpdQSEDJWJYfDO4RvYeakWS1ND8fqCUZQEzLQ8LQyvzo3H2QoFnv93ETR6I9chETtCv2XEJvQmgT2X6/GzyeH4n9kx1J9/iB4ZH4w/LByF/Jp2PPvlFaj1Bq5DInaCEgHh3I+TwNMzoqgCZpgWjA7EWz8Zg6v1HXj+31eh7aYnAzI4SgSEUwzDYNPRm5QELOi+eD+sX5iIwtp2vPDVVegMJq5DIjaOEgHh1N/PyLGroA4rJoZRErCgeYkBeG1eAs7LlXgp5xq6jZQMyN1RIiCc6dlIpgoPJwVhzcxoSgIW9pOkILw4Jw6nytuw/mApTLZZIEhsgOO1bCR24ZtrjXjv2E3MjvPFS3PjKQmw5LHkEHTqDPjgVCUCJGI8OzOG65CIDaJEQKzuvFyB9QdLMTFChjcXjXbavkHW8rPJ4Wjq1OGTizXwl7hgaWoo1yERG0OJgFjVzRY1Xtx7DVE+bnj3oTFwEdLoJNt4PB5euC8OLWo93jt2E/4SMeYk+HMdFrEh9FtIrKZVrcfz/y6Cq0iAvzyaRBuzW5GAz8MbCxMxLsQTv99fgqsNtHKffI8SAbEKbbcRa/9zFQpNN/786FgEebpyHZLTcRUJsOnhMfB1F+GF/1xFU+edN44izocSAWEdwzB489tSXGvoxJuLEjE6cOi9UIhleLuL8adHk6DRG/HCV7TgjPSgREBYtyOvFgeLm/HrGVGYFefHdThOL87PA28sSkRxowp/OFAKG+07SayIEgFh1flKBTafLMecBD/8bHI41+GQW2bG+uKZ9GgcLm3GJxdruA6HcIwSAWFNjbILr3x9HdG+7siaP4rWCtiYzElhuD/BHx+cqsDFKgXX4RAOUSIgrNAZTHhx7zUAwKaHx8JdTP3xbQ2Px8Nr8+MR4e2GV/cVo5Emj50WJQLCiveO3URpsxp/WDAKYTI3rsMhd+EhFuKPD42FzmDCy9STyGlRIiAWd/B6E/ZcrsfKSWGYEePLdThkENG+7sh6IAFX6jvx/skKrsMhHKBEQCyqsk2Dtw6VYUKIJ349PYrrcIiZ5iT4IyMlBDsv1eLkzVauwyFWRomAWIy224iXc65DfGtTddpm0r6smRmDBH8PrD9QQovNnAxrv6mFhYXIzMy87fV9+/bh8ccfx9KlS5GVlQWTicYkHcWW7ypwo0WN1xeMQqDUhetwyBD1JnC90YR13xTDaKL1Bc6ClUSwdetWvPbaa9Dp+n+q0Gq1+Mtf/oJPPvkEO3fuhEqlwrFjx9gIgVjZmYo2ZOfXISMlBNOjfbgOhwxTlI87fjcnDpdq2rHtfBXX4RArYSURREREYPPmzbe9LhaLsXPnTri59VSRGAwGuLjQJ0d7p9Do8YcDJYj1c6d+9w5g0ZhAPDA6AB+dleNqfQfX4RArYCURzJ8/H0Lh7Z0l+Xw+/Px6Wgxs374dGo0G06dPZyMEYiUMw+CNg6VQ6Qx4c+FoaivtAHg8Hl6cEwc/iQuy9pdQPyInYPXfWpPJhHfeeQenT5/G5s2babWpnfv3lQZ8V96GZ2bGIM7fg+twiIVIXIT4/QMJqFJ0YTOVlDo8qyeCrKws6HQ6fPDBB31DRMQ+1bVr8dfj5ZgYIUNGSgjX4RALmxThjWWpofiioA7nK6kFhSOzSiLIyclBdnY2rl69it27d6O0tBSrVq1CZmYmDh06ZI0QiIWZGAZvfFsKAMianwA+Pdk5pKdnRCHaxx3rD5agQ9vNdTiEJTzGRnvQNjfTDkq2bFdBHf545AZemRuPR8cHcx0OYVFxYyd+tiMfC8cEIuuBUVyHQwbh7z/0/T5oZo8MWY2yC++fKMfUSG88Mi6I63AIyxIDpVg5ORw5VxtxtrKN63AICygRkCHp3W1MwOfh1XnxNNnvJJ6aGokoHze89W0Z1HoD1+EQC6NEQIYk52oj8qrbsWZWDO077ERchHysmz8KjZ06bKEqIodDiYCYTaHR4/0T5ZgQ4klDQk5ofIgnMlJDsbuwHnnVSq7DIRZEiYCY7c/Hy6HWG/HKvHiqEnJST8+IQoiXK94+VAa9gfqEOQpKBMQs5+UK7L/ehJWTwxHjSwvHnJWbSIDfzYmDXNGFT3Npr2NHQYmADErbbcTGw2UIl7niySkRXIdDODY92gf3xfvhn+erUNvexXU4xAIoEZBB/fN8FWqUWrx0fzz1EiIAgOdnx4DPAzYdvQkbXYpEhoB+q8mAbrSo8cnFGiwaE4DJkd5ch0NsRJCnK1bfE4VT5W04foN2NLN3lAjIXTEMg3cOl0EiFuA3s2K5DofYmKUpIYjz88Cmozeg0VOHUntGiYDc1cHiZhTUduCZ9GjI3EVch0NsjFDAx0v3x6FJpcfWs3KuwyEjQImA3JFGb8T7J8sxOlCCh2jNALmLCaFeeDgpCJ/n1eBGs5rrcMgwUSIgd/SvC1VoVumx9t5YWjNABvTMzGhIXIT449EbNHFspygRkNvUKHtqxBeMDsCEUC+uwyE2TuYmwq9nRCG/ph3HaOLYLlEiILf5y/FyCPk8PDszmutQiJ14eFwwYv3c8dcT5bTi2A5RIiD9nKtsw4mbrXhySgT8JS5ch0PshJDPw/OzYlHXrsXOS7Vch0OGiBIB6WMwmvDesXKEyVyxPC2M63CInZkS5Y30GB/883wVWtV6rsMhQ0CJgPT5oqAOFW0aPD87FmJaQUyG4blZMdAaTPjb6UquQyFDQL/tBADQ3tWNj85WYeqtT3WEDEekjzueSA7BV1caUNKk4jocYibWEkFhYSEyMzNve/3o0aNYsmQJMjIy8MUXX7B1ejJE285XQ6Uz4LmZMbTrGBmR/5oWAU9XIf58nPoQ2QtWEsHWrVvx2muvQafT9Xu9u7sbb7/9Nv75z39i+/btyM7ORnNz8x3fgyoPrKeuXYsvCmrx4NhAxPlTi2kyMp6uIqy+Jwp51e04QeWkVsMwDIrqO4b1s6wkgoiICGzevPm212/evImIiAh4eXlBLBYjLS0Nubm5d3wP6nVuPX87XQk+j4dfTo/iOhTiIBZPCEa0rzs2f1cBg5E+1FnD2UoFfv5ZwbB+lpVEMH/+fAiFwtteV6lUkEqlfV97eHhApbrzOOKFKgUboZEfKWlUYf/1JixNDUWglMpFiWUI+Tz894xoVCm6sPdqI9fhOIULciXEguEN61p1slgikUCt/r4fiVqt7pcYfuhKXQd0NDzEus3flcPLVYhVk8K5DoU4mJmxPhgf4omtZ+TQdlN3UrZdqlEiKdhzWD9r1UQQGxsLuVwOpVIJvV6P3NxcpKSk3PF79cbhj3cR85yrbMN5uRJPTo2A1PX2JzhCRoLH4+GZ9Gi0qPXIzq/jOhyH1qk1oKRJhbTw4bWEsUoiyMnJQXZ2NkQiEV566SU89dRTWLp0KZYsWYLAwMA7B8YD8qqV1gjPKZkYBu+frECIpwsemxDCdTjEQaWEeWF6tA8+vlCNDm031+E4rILadpgYIC1cNqyf5zE2Wt+14M8n4CoS4MOMCVyH4pC+udaI3+8vwRsLE/HA6ACuwyEOrKxZhZ9+cgmZk8KpfxVL/nqiHNn5tTj63/cgPEQ25J+32QVlaeEyFNV30NgiC/S3Vn4mBkgwL9Gf63CIg4v3l2D+6ABk59eiqVM3+A+QIcur7pkfcBUJhvXzNpsIJobL0G1kUFTfyXUoDuerogbUd+jwdHoU7TVArOKX90TCaGLwj3NVXIficFS6W/MDYcNvGW+ziWBCqCf4PCCX5gksSmcwYdv5KkwI8cRU2oyeWEmYzA2Lxwfjqyv1kLdpuA7HofTOD6QOc6IYsOFEIHERIjFQikuUCCxqz+V6NKv0+NX0KGolQazqyakREAv5+Ntp2t/Yki5Vt0Mk4GHcMEtHARtOBAAwMdwLV+o7aZ7AQrTdRvzrfBXSwr0wMULGdTjEyfh6iLEsNRRHSptxo4X2N7aUvJp2JAVJhz0/AAADFo9v2bJl0Dd45plnhn3ywaSGy/DJxRpcruvAZBrGGLFdBXVo03TjnZ9EcR0KcVLL0sKQnV+Hf5ytwts/Gc11OHZPpTOguLETP5sSMaL3GTARTJ48eURvPlLJoZ4Q8HoyHiWCkVHrDfjkYg2mRnojeQSTSoSMhMxNhIyUEGw7X40bLRGI86MmhyNRWNvRs35ghL/TAw4NTZ48GZMnT4ZMJkNpaSkuXLiAGzduwN/fv+8YmzzEQowOkiKvSsnqeZzBF/l1UHZ141fTI7kOhTi5ZWlhcBcL8I+zNFcwUnnVSgj5PIwPGf78ADBIIrh58yZ++ctfYsuWLVAqlQgICIBarcZ7772HX//61ygrKxvRyc2RGibD1YZOdNE8wbCpdAZ8mluDGTE+GDuCCSVCLKH3qeBwaQvNFYxQbrUSScEjmx8ABhka+vrrr7Fp06Y7NoZrb2/Hv/71Lzz33HMjCmAwaeFe+ORiNS7XdmBKFA0PDcfnebXo0Brwy3voaYDYhu/nCuR4+ydjuA7HLvX2F3pyhPMDwCBPBGvWrIFUKu3XMRQA6uvr4eXlxXoSAIDkUK9b8wRK1s/liDq1Bnx2qQaz43yRGHjnTq+EWBs9FYzcpZqe9QOWqAA0q3x02bJlKCkpAQAcPHgQTz755IhPbC53sQBjgqTIq2632jkdya6COqh0RvzXVHoaILZlWVoYPGiuYNhyq5VwEfJHtH6gl1m9h9977z28+uqr8PX1hVAoxI4dO0Z84qFIDZfh09waaPRGuItHNhbmTLq6jfgsrwbTo30wKlDCdTiE9NP7VPDP89V4qkVNFURDlFetxPgQT4iFI18OZtY79DYo1ev14PF4EAisezOeGO4Fo4nB5Tp6KhiKf1+uR7vWgJ9PoU1niG3qfSr4iJ4KhkSh0aOsWY1JFloYalYi+M1vfoPXXnsNH330ERYsWIDly5db5OTmGh/iBQGfR8NDQ6A3mPBpbg3Swr0wIZTWDRDbJHMT4fHkEBwtbUEl9SAyW++9cOIw9x/4MbMSQXZ2NsaNGwcAeOCBB7B161aLnNxc7mIBxgZJaaOaIdh3tQHNKj1+boGKAkLYtCwtFGIhH9svVnMdit3IrVbCXSTAaAsN+Q6YCH73u9/h+PHjcHNz6/d6UFAQDh8+jBdeeMEiQZgjLdwL1xo6odHTeoLBGEwMPr5Yg7FBUkymnkLExvm4i/FQUhC+udaERtqvwCy5VUqkhHlBKLBMu7gB3+XNN99EWVkZHnnkESxfvhxr1qzBqlWr8PDDD6O8vBwbNmywSBDmSAuTwcj0tFwlA/u2uAl17Vr8fEoEdRgldmHFxDAwDIPP8mq4DsXmNXXqIFd0WbRx5IBVQ2KxGL/4xS/wi1/8ApWVlVAoFPD19UVEhPWHG8aHekJ4a57gnmgfq5/fXpgYBv86X404Pw+kx9J1IvYhxMsV8xID8O/L9fj5lAjI3ERch2SzevdomTiC/Qd+zKzy0R93IRWJRAgKCsLChQshEt3+B2YymfD666+jpKQEYrEYb775JiIjv69j37t3L7Zt2wY+n48lS5aYNfnsJuqZJ7hEC8sGdLysBRVtGmxYlEi7jxG7smpyOPZfb8Ku/Dr8glbB31VetRKerkLE+1uuJNysAaaSkhJUVlbCz88PtbW1OHv2LE6dOoVXXnnljt9/+PBh6PV6ZGdnY+3atdi4cWO/43/84x+xbds2fP7559i2bRva280b7kkL98L1hk6odAazvt/ZMAyDbeerEeHthjkJtBcxsS+xfh6YGeuL7PxamgscQG6VEqlhPZWUlmJWIujo6MCmTZuwdOlSvP322+Dz+Xj33XdRU3Pn8by8vDykp6cDAJKTk1FUVNTv+KhRo9DZ2Qm9Xg+GYcwex04Lp3mCgVysUqK4SYXMiWEW/UtCiLWsmhyOdq0B/7lSz3UoNqm2vQt1HTqLlY32MisRdHZ2oq2tDQCgUCjQ2dmJ7u5uaLXaO36/SqWCRPL9Y4tAIIDB8P2n+Pj4eCxZsgSLFi3C7Nmz4elp3hLp8SGeEAt4uEhtqe9oe24NfNxFWDAmkOtQCBmW8SGeSA3zwo7cGnQbTVyHY3Pyqm6tH7BwNaBZieDZZ5/FE088gUceeQQZGRl49tlnsW3bNjz22GN3/H6JRNKvUZ3JZIJQ2DMdUVxcjOPHj+PIkSM4evQo2trasH//frOCdRUJMD7UixLBHZQ1q3CuUoGlqaFwscCSc0K48rMp4WhS6bH/WhPXodici9VK+LiLEOPrbtH3NWuy+N5778WsWbPQ3NyMgIAA8Hg8zJw5867fn5qaimPHjmHhwoUoKChAQkJC3zGpVApXV1e4uLhAIBDAx8cHHR0dZgc8OUKGD05Vok2jh4+72Oyfc3Q78mrhKuRj8fhgrkMhZESmRnpjVIAEn1ysxoNJgVT0cAvDMMitUmJiuMziZeFmf3Tk8/n47W9/a1YAc+fOhVgs7ptTePnll5GTk4Ps7GyEhoYiIyMDy5cvx7Jly9DZ2YlHH33U7IB7e2vk0lNBn6ZOHQ5eb8LD44LgRWV3xM7xeDxkTgyDXNGFMxVtXIdjM+RtXWhR65HGwiJRs54IevU2nxsMn8/H+vXr+70WGxvb99/Lli3DsmXLhnLqPomBUniIBbhYpcS8xIBhvYejyc6vhYlhsDQ1lOtQCLGIOQl+eP+kGDvyajEjxpfrcGxC7/qBSRaeKAaG8EQA9PQZ4pqQz0NauIzmCW5R6Qz4srAe98X7IUzmNvgPEGIHhAI+lqaGIrdKiZImFdfh2ITcaiUCpS4Ik7la/L0HTQQqlQrbt2/HmjVrcODAAfzmN79Bdnb2bbuWWdOkCBlq27Wobe/iLAZbsbeoAWq9ET+dGMZ1KIRY1CPjguEm4lPbCfR0DOiZH/BipW3MgIngyy+/xPPPP98zZpeZiTfeeANPPvkk9Ho9nnvuOezatcviAZljcqQMAHBRruTk/LbCYDTh87xapIR6Iok2pScORuoqxENJQThY3IwmJ29Gd7NFjXatweJlo70GnCPw8/O7Y8vp8ePHIzMzEydOnGAlqMFE+7jDz0OMi1VKPOLEVTKHS1vQ0KnDb+fEcR0KIaxYmhqKXQV12FVQh/9Oj+Y6HM70DoVbeiFZrwGfCGbNmgWgp/b/h44ePdrvuLXxeDxMjJAht1pp9gS2o2EYBp/m1iDS2w0zYqi5HHFMYTI3zI7zw57L9ejqdt62E7lVSoTLXBHkafn5AcDMyeJXXnkFu3btgl6vxxtvvIHt27ezEsxQTIqQoU3TjZstzrmrUX5tO0qaVFieFkp11sShLU8LRYfWgH1XG7kOhRMGE4NLNe1IY+lpADAzEXz22Wf47rvvcO+998Lf3x/btm1jLSBz9W64cqFKwW0gHMm+VAdPVyEWUjsJ4uDGh3giKViKz/NqYDQ53whASZMKar3RYvsT34lZiSAnJwcVFRVYtWoVDhw4gLy8PNYCMleQpyvCZa5OWUba0KHF8RsteGRcEFxFAq7DIYRVPB4Py9PCUK3U4lR5K9fhWN0Fec+HXc6fCE6dOoXPPvsMq1evxubNm/Huu++yFtBQTIrwRn5NOwxO9ilhV0FPZ8bHkkM4joQQ67g33g/Bni7YkVfLdShWd6FKiXh/D/h6sNdSZ8BE0NjYMyb317/+FVKpFAAQHh6OTz/9tN9xrkyKkEGtN+JaQyencViTttuIr67UY1acH4JZmjgixNYI+Tw8kRKK/Jp2lDrRAjNttxGFte2sDgsBgySCjz76CO+88w6uXbuG7u5uAD3VKsXFxdiwYQM+/PBDVoMbTG8p1UUnmic4cL0J7VoDMlLoaYA4l4eSAuEi5GNXQR3XoVhNQW07uo0MpkR6s3qeAdcRvPrqqygsLMQ//vEPXLx4ESaTCa6urkhNTcXy5cuRnJzManCDkbmLMCpAgotVSjw11fG3tmMYBtn5dYj390BqmOX2KyXEHni6irBgdAD2X2/CM+nRTtFg8bxcCZGAhxSWf98HTASffvopVqxYgczMTPzpT39iNZDhmhQhQ3Z+LbTdRoefOL1U044bLWqsm5fAyjJzQmzd48kh+M+VBuRcbcQKJ2irckGuwPgQT7ixfG8bcGgoOzsbx48fx+9//3ucOnWq3z+2YlKEDN1Gxim2r9x5qRZerkLMS6T9iIlzSgiQICXUE7sL6hy+lLRNo0dps5r1YSFgkCeC5557DocPH0Zrayu+/vrrfsdmzJjBamDmSgnzgpDfs33l1CjHXWFb167FyZutWDkp3OGffAgZyOMpoXhl33WcrWxz6BbVvb3UJrM8UQwM8kRw5MgRGI3GvoniH3r55Zfx8ssvsxaYudxEAowLljr8eoLdBXXgAVgywXl7KxECAPfG+cJfIkZ2vmNPGl+oUkDqIkRioJT1cw34RPDMM8/0+7etmhThja1n5Wjv6nbICSRttxFfFTXg3ng/1nqNEGIvhIKeLVn/fkYOeZsGkT6W3b/XFjAMg/NyJSZGyCDgsz8fOOATQWho6KD/2IJJETIwAPJqHHOe4NviZnRoDXicSkYJAQA8Mj4YQj4PuwvruQ6FFVWKLjR26jDlVst9tg1phzJbNTZYCjcRHxfljrme4MvL9YjxdUdKKJWMEgIAfh5izEnwQ05RAzR6x+tKev7W/IA1JooBlhKByWRCVlYWMjIykJmZCblc3u/45cuX+zavX7NmDXS6kW06IRLwkRLmhQsOOE9wvbET1xo6sWRCMJWMEvIDGSmhUOuN2H/d8bqSXqxSIMTTBaFe1hkKZiURHD58GHq9HtnZ2Vi7di02btzYd4xhGKxbtw5vv/02Pv/8c6Snp6O2duT9Q6ZEeqNK0YWGDu2I38uWfFlQD1chn7qMEvIjScFSjA6UIDu/zqH2JTGYGFysUmJypLfVPvyxkgjy8vKQnp4OAEhOTkZRUVHfsYqKCshkMnz88cdYsWIFlEolYmJiRnzO3keoc5WOMzzUqTXgYHET5o8OgMRlwHl9QpwOj8fDYxNCUNGqQUFtB9fhWMy1hk6o9UZMttKwEMBSIlCpVJBIJH1fCwQCGAwGAIBCoUB+fj6WL1+Obdu24dy5czh79uyIzxnj644AiRjnHWie4JtrjdAaTFQySshdzE30h4dYgD2XHWfS+IJcAR7AeqO5H2IlEUgkEqjV6r6vTSYThMKeT7QymQyRkZGIi4uDSCRCenp6vyeG4eLxeJgS6Y0LVUqHWHHIMAy+LKzHmCApRluhjpgQe+QmEmDhmEAcLW2Gsuv29U726IJcgcRACWRWLIVnJRGkpqbi5MmTAICCggIkJCT0HQsPD4dare6bQM7NzUV8fLxFzjs1yhsdWgOKG+2/LXV+bTsq2jT0NEDIIBaPD4beyOCba/Y/aazRG3G5vhOTIqw3LAQMsqBsuObOnYvTp09j6dKlYBgGb731FnJycqDRaJCRkYENGzZg7dq1YBgGKSkpmD17tkXOOznCGzwA5+QKjA32tMh7cuXLgnpIXYSYN4r6ChEykDh/D4wLlmJPYT2WpYbadXXdpZqeEQ1rrR/oxUoi4PP5WL9+fb/XYmNj+/572rRp2L17t8XPK3MXITFQgvOVCrtuS92q1uNoWQseSw6hvkKEmOHR8cFYf7CU9U3e2XZBroSLkI8JVl4z5BALyn5oSqQ3Ltd3QqUzcB3KsOUUNcBgYrB4PA0LEWKOuaP8IXER4N92Pml8Xq5AcqgnXITWvTU7XCKYGuUNo4lBXrV9tpswMQz+fbkeaeFeiPZ1vB4qhLDBVSTAojGBOFrWAqXGPieNm1U6lLdqrLaa+IccLhH0bOLAt9sy0vNyBeo6dPQ0QMgQPTI+GN1GBjlXG7gOZVh6OyhPtvJEMeCAiUAk4CMtXIZzlW1chzIse680wMtViNlxflyHQohdifPzwPgQT/znSoNdrjQ+L1fA202E+AAPq5/b4RIBAEyN9Ea1UosaZRfXoQyJQqPH8RutWDgmEGIrjxES4ggWjw9GlaLL7oaGGYbBBbkSkyJk4HNQ9eSQd5spUT2PVhfsbHho//UmGEwMHhoXxHUohNilOQl+8HQV2t1K45utGrSo9Zhs5bLRXg6ZCCK93RAkdcG5W61c7QHDMPjPlQYkBUsR52f9R0NCHIHrrZXGx+xs0ri3RxpX2+06ZCLg8XiYEuWNi1UKGOyk3cSV+k5UtGrwcBI9DRAyEg8nBcFgYrC/uInrUMx2tqINMb7uCJS6cHJ+h0wEQM88gUpnxLUG+2g3sfdKA9xEfMxNpJXEhIxEnL8HRgdKsNdOJo213UYU1LZjapT1q4V6OWwi6Jl0Ac7bQVtqtd6Ab0uaMHeUPzzE1G6akJF6KCkIN1rUKG5ScR3KoPJq2qE3MphGicDyvNxEGBMkxTk7mDA+XNKMrm4THh5HawcIsYT5iQFwEfKx94rtryk4W9EGFyEfKWEyzmJw2EQA9LSbuFrfgU6tbbeb+OpKA6J93DEumNpNE2IJUlch7o33w8HiZmi7bXtP43OVCqSGeVm9rcQPOXQimBrpDSMDXKxWch3KXd1sUeNKfSceHhdk110TCbE1DyUFolNnwIkbrVyHcld17VrIFV2czg8ADp4IkoKl8BALbHqeYG9RA4R8HhaOCeA6FEIcSlq4DCGeLthbZLvDQ70dEKZxVDbay6ETgVDAx8Rb7SZssXpAbzDh66uNmBXnC293MdfhEOJQ+DweHkwKwsUqJeratVyHc0dnKxUIkrogyseN0zgcOhEAPauM6zp0qFHa3l+EU+WtaNca8BCtHSCEFQ+ODQQAfH3V9nYvMxhNuFilxNQob86HhR0+EfSWZNli9dC+q43wl4g5aTtLiDMI9nTF5EgZcq42wGRjowJX6juh1hs5LRvt5fCJIEzmhlAvV5ytsK1upG0aPc5UKrBgdAAEfJokJoQtDyUFob5D19fm2Vacq2yDgAer7098Jw6fCICep4LcaiX0BhPXofQ5cL0JRhODRbceXQkh7JgV5wepixA5NjZpfLZSgaRgT0hduV9EykoiMJlMyMrKQkZGBjIzMyGXy+/4fevWrcOmTZvYCKGfe6J90NVtQkGt7bSm/fpqI0YHShDjSw3mCGGTi5CPB0YH4FhZi82sKVJo9ChuVHFeNtqLlURw+PBh6PV6ZGdnY+3atdi4ceNt37Nz506UlpaycfrbTIyQQSTg4UyFbcwTlDWrUNqs7pvIIoSw68GxgdAbGRwubeY6FADAebkSDGAT8wMAS4kgLy8P6enpAIDk5GQUFRX1O56fn4/CwkJkZGSwcfrbuIkESA3zwhkb2bVs39VGCPk8zEuktQOEWMPoQAmifNzwzTXbqB46V9kGL1chEgNto5sAK4lApVJBIpH0fS0QCGAw9DySNTU1YcuWLcjKymLj1Hd1T7QPKlo1aOjgtozUYDThwPUmzIjxgcxNxGkshDgLHo+HhWMCUVDbwfnOhSaGwdlKBaZEettMoQgriUAikUCtVvd9bTKZIBT2TIgcOHAACoUCq1evxocffoh9+/Zhz549bITRzz23Vu6d4bh66GylAm2abhoWIsTKFowOAA89OwFyqaxJjTZNt83MDwAsJYLU1FScPHkSAFBQUICEhIS+YytXrsSePXuwfft2rF69Gg8++CAWL17MRhj9RPq4IcTThfN5gq+vNULmJsI90dwuKSfE2QR5uiItQoZvrjVy2mmgd4h6mg3dA1hJBHPnzoVYLMbSpUvx9ttv4+WXX0ZOTg6ys7PZOJ1ZeDwepkX74GKVEt1GbspI27u6cfJmK+Yn+kMkcIrKXUJsysLRAahRanGlnrsNq06Xt2F0oAR+HrbTVoaVAlY+n4/169f3ey02Nva277PGk8APTYvywZeF9SiobedkEcehkmZ0GxkaFiKEI/cl+OGdIzfwzbVGjA/xtPr527u6caW+Az+fEmH1cw/EqT6WTuK4jPTra42I9XPHqADJ4N9MCLE4D3HPPgWHSpo5WWB6Xq6AiQGm29CwEOBkicBdLEByqBcnE8aVbRoU1XfiwbG07wAhXFo0JgAdWgNOcXAfOF3RUzY6Jsg2ykZ7OVUiAHrKSMs5KCM9cL0JfB7wAG1OTwinJkV4w89DjG+s3JHUxDA4W6HA1CjbKRvt5YSJoGdu4IwVN6thGAYHrjdhYrgMfhIXq52XEHI7AZ+HB0YH4FRFG5Sabqud93qjCoqubkyPsa1hIcAJE0G0jzuCpC5W7UZaVN+J2nYtHhhNK4kJsQWLxgTCaGLwbYn1Wk6cKW8DD8C0SEoEnOPxeLgn2gcX5NYrIz1wvQliAQ/3xvtZ5XyEkIHF+Xsg3t/Dqi0nTle0ISlYCpm77XUUcLpEAPQMD2m6jSis7WD9XAajCYdKmpEe6wuJC/ftZgkhPRaNCcTVhk7I2zSsn0uh0eNaQ6fNLiR1ykQwMUIGIZ9nleqh81VKKLq68QA1mCPEpsxL9AcPwLfF7A8Pna1UgAFscn4AcNJE4CEWIiXMyyrlYweuN0HqIrTZTwKEOCt/iQvSwr1woLiJ9ZYTZyra4OMustk1RE6ZCABgRkxPN9LadvY6EXZ1G3HiRgvmJPhBLHTaS02IzZqfGIAqRReKm1SsncNo6uk2Oi3aB3wbXUPktHenGTG+AHr6frDlxI1WdHWbqFqIEBt1X4IfhHweDl5nb3ioqL4DHVqDza0m/iGnTQQR3m6I8HbDKRYTwYHrTQiUuiAlzIu1cxBChs/TtacT8KGSJphYGh46U9GzSf2USBkr728JTpsIgJ7hodxqJTR6o8XfW6HR41xlG+Yn+tvs4yAhBJif6I8mlR75NezsaX6mQoHxIZ7wdLW9stFeTp8Iuo0MLlZZfpXxoZIWGBnQsBAhNm5mrC/cRHwcLLb8hjUtKh2Km1Q2Xyzi1IkgOdQLHmIBK8NDB643IdbPHfH+tlklQAjp4SoSYFacH46Utlh8kenpW5WJtlo22supE4FIwMfUKG+cKm+zaPlYbXsXrtR30NoBQuzEA4k9HUnPWbgH2Xc32xAkdUGcn4dF39fSnDoRAD3DQy1qPUosWD7Wu0BlHiUCQuzClEgZvFyFFh0e0hlMOC9XYEaMj823nnf6RHBPtA94gEWHhw6VNGNcsBQhXq4We09CCHuEAj7uH+V/q+TbMsUjudVKaA0mpMf6WuT92OT0icDHXYyxwVKLJYLKNg3KmtW4fxTtO0CIPZmX6A+twYSTN1ot8n7f3WyFm4iPtHCZRd6PTawkApPJhKysLGRkZCAzMxNyubzf8X379uHxxx/H0qVLkZWVBZOJm83ke02P9sHVhk60qvUjfq9DJc3gAbg/gRIBIfYkOdQLARIxDlhgeIhhGJwqb8OUSG+42EFXAVYiPHz4MPR6PbKzs7F27Vps3Lix75hWq8Vf/vIXfPLJJ9i5cydUKhWOHTvGRhhmS7+1ytgSTegOlTQjOcwLAVLagIYQe8Ln8TA/MQBnKxVQdo1sw5qyZjUaO3V99xZbx0oiyMvLQ3p6OgAgOTkZRUVFfcfEYjF27twJNzc3AIDBYICLC7c3zYQAD/hLxCMeHrrRokZFqwZzaViIELs0L9EfRhODEzdaRvQ+vfeSe2y8bLQXK4lApVJBIvm+fl4gEMBgMPSckM+Hn1/PBi3bt2+HRqPB9OnT2QjDbDweD9OjfXBerhhRHfGh4p59ie+jDWgIsUujAiQIk7nicMnIEsF35a0YGySFn4fYQpGxi5VEIJFIoFar+742mUwQCoX9vn7nnXdw+vRpbN682SZKq2bE+EKtNw57mTnDMDhU0oy0cBl87eQPnxDSH4/Hw/0J/rhYpRj2fsataj2u1ndihp08DQAsJYLU1FScPHkSAFBQUICEhIR+x7OysqDT6fDBBx/0DRFxbXKkDGIBb9jDQyVNKlQrtTQsRIidu3+UP4wMcGyYw0OnK9rAAHZRNtqLlb0T586di9OnT2Pp0qVgGAZvvfUWcnJyoNFokJSUhN27d2PixIlYtWoVAGDlypWYO3cuG6GYzU0kwMQIGb4rb8Xzs2OG/JRyqKQZAj7tS0yIvUvw90CEtxsOlzTj0fHBQ/75U+VtCJCIkeBv26uJf4iVRMDn87F+/fp+r8XGxvb9d3FxMRunHbGZsb7YePgGKto0iPE1/w+xd1hoSqQMMjfb7TBICBlcz/CQH/51oRoKjR7e7uYP9eoNJpyvVGDBmACbGPI2l+0XuFpRb6nXiSEuKCmq70R9h46GhQhxEPeP8oeJAY6VDW14KK9GCU230a7mBwBKBP0ESF0wOlCC724OLREcKmmGSMDD7DgaFiLEEcT5eSDS2w2HSoeWCE7dbIOLkI+JdrCa+IcoEfzIzFhfFNV3osXMVcYmhsGR0mZMi/KBxIWVkTZCiJXxeDzMGeWPS9VKszsOMAyD78pbMSXSG64iAcsRWhYlgh+ZGesLBsApM58KCms70KTSYx4NCxHiUOYmDG146GarBvUdOrsbFgIoEdwm3t8DwZ4uOGlmIjhS2gyxgIcZsfb3h08IubtYP3dE+bjhcKl5G9v3NqtLp0Rg/3g8HmbG+uJClXLQdrQmhsGxshZMi/KBh5iGhQhxJL2Lyy5Vt5s1VHz8RgvGBUvhJ7G/PmOUCO5gZqwvdAYTLsgH3q3oan0nmlR63JdAk8SEOKL7R/mDAXB0kEnjxk4drjeqMNOOFpH9ECWCO0gN84LERTDo8NCR0hYI+Ty76TBICBmaWD8PRPu6Dzo81HuvmGWnlYOUCO5AKODjnigffHezDUbTnfcyZhgGx8qaMSXSG1JXGhYixFHNTfBHQU07mlW6u37PyRutiPB2Q5SPbbTMGSpKBHcxK84Xiq5uFNV33PF4cZMKdR06GhYixMHNGeUHBsCxsjuPEKh0BuRWKzEr1teuVhP/ECWCu5gW5QMBn3fX4aEjpS0Q8GC3Y4KEEPPE+HogyscNx8ruPDx0pqINBhODWXH2ey+gRHAXUlch0sK87pgIGIbB0dJmTIyg3kKEOIP74v1wqaYdCs3t1UPHb7TCx12EpGBPDiKzDEoEA5gZ64vKti7I2zT9Xr/Roka1Uov7aF9iQpzCffE9i8t+3IdMbzDhTEUb0mN8IeDb57AQQIlgQDNvPer9+KngSGkL+Dxgth0/ChJCzJcQ4IEQL1cc/dEq47waJdR6o10PCwGUCAYU7OmKeH+P2xLB0bIWpIR5wWcI7WkJIfaLx+Phvng/XKxSolNr6Hv9xI1WuAr5mBQh4y44C6BEMIhZsb4orO1A262xwfLWng3qaV9iQpzLffF+MJh6GssBPZ0FvrvZimnRPnbXZO7HKBEMYnZ8T+lYb2vq3hWGtBMZIc5lbLAUARJx3z3geqMKTSo9ZjlA5SAlgkEk3GpCd/zWJNHRshZMCPGEvx32EyGEDB+f17MV7Tm5Ahq9ESdv9JSQT7fDJnM/RolgEDxez4YzF+QKXKxSoKxZTYvICHFS9yX4QWcwYf/1Rhy/0YrkMC+HKCGnRGCGWXG+0BsZPL3rCjzEAsyhslFCnNKEEC9E+bhh4+EbKG/VOMyCUlaa5JhMJrz++usoKSmBWCzGm2++icjIyL7jR48exf/+7/9CKBRiyZIleOKJJ9gIw2KSQ72waGwggqUueCw5BL4eVC1EiDMS8Hn4ZEUq9l1txPlKBR4YHcB1SBbBYxjmzl3VRuDbb7/F0aNHsXHjRhQUFODvf/87/u///g8A0N3djYULF2L37t1wc3PDsmXL8Le//Q3+/v0/ZTc3d1o6LEIIcXj+/tIh/wwrTwR5eXlIT08HACQnJ6OoqKjv2M2bNxEREQEvLy8AQFpaGnJzc7FgwYJ+7zGc/xlCCCFDx8ocgUqlgkQi6ftaIBDAYDD0HZNKv7/Je3h4QKVSsREGIYQQM7CSCCQSCdRqdd/XJpMJQqHwjsfUanW/xEAIIcS6WEkEqampOHnyJACgoKAACQkJfcdiY2Mhl8uhVCqh1+uRm5uLlJQUNsIghBBiBlYmi3urhkpLS8EwDN566y1cu3YNGo0GGRkZfVVDDMNgyZIl+OlPf2rpEAghhJiJlURgLkcrMx2Jwa7Fvn378PHHH0MgECAhIQGvv/46+HzHXAYy2LXotW7dOnh5eeGFF17gIErrGOxaXL58GRs3bgTDMPD398e7774LFxfHXPU+2LXYu3cvtm3bBj6fjyVLlmD58uUcRmsdhYWF2LRpE7Zv397v9SHfOxkOHTx4kHnxxRcZhmGY/Px85le/+lXfMb1ez9x///2MUqlkdDods3jxYqapqYmrUFk30LXo6upi5syZw2g0GoZhGOb5559nDh8+zEmc1jDQtej1+eefM0888QTz7rvvWjs8qxroWphMJuahhx5iKisrGYZhmC+++IK5efMmJ3Faw2B/L6ZPn84oFApGp9P13Tsc2Ycffsg8+OCDzOOPP97v9eHcOzn9SGlumalYLO4rM3VUA10LsViMnTt3ws2tZ2Nsg8HgsJ/6gIGvBQDk5+ejsLAQGRkZXIRnVQNdi4qKCshkMnz88cdYsWIFlEolYmJiuAqVdYP9vRg1ahQ6Ozuh1+vBMIzd7h9sroiICGzevPm214dz7+Q0EVCZ6fcGuhZ8Ph9+fj39jbZv3w6NRoPp06dzEqc1DHQtmpqasGXLFmRlZXEVnlUNdC0UCgXy8/OxfPlybNu2DefOncPZs2e5CpV1A10LAIiPj8eSJUuwaNEizJ49G56e9rt1pDnmz5/fV435Q8O5d3KaCKjM9HsDXYver9955x2cPn0amzdvduhPOwNdiwMHDkChUGD16tX48MMPsW/fPuzZs4erUFk30LWQyWSIjIxEXFwcRCIR0tPTb/uU7EgGuhbFxcU4fvw4jhw5gqNHj6KtrQ379+/nKlRODefeyWkioDLT7w10LQAgKysLOp0OH3zwQd8QkaMa6FqsXLkSe/bswfbt27F69Wo8+OCDWLx4MVehsm6gaxEeHg61Wg25XA4AyM3NRXx8PCdxWsNA10IqlcLV1RUuLi4QCATw8fFBR0cHV6Fyajj3TlZaTJhr7ty5OH36NJYuXdpXZpqTk9NXZvrSSy/hqaee6iszDQwM5DJcVg10LZKSkrB7925MnDgRq1atAtBzQ5w7dy7HUbNjsL8XzmSwa7FhwwasXbsWDMMgJSUFs2fP5jpk1gx2LTIyMrB8+XKIRCJERETg0Ucf5TpkqxrJvZPT8lFCCCHcc8xCdEIIIWajREAIIU6OEgEhhDg5SgSEEOLkKBEQQoiTo0RACCFOjhIBIYQ4OUoEhJhpx44dWLt2LQDgxRdfxI4dOziOiBDLoAVlhAzB008/DU9PT+j1erz33ntch0OIRVAiIGQICgoKkJGRgT179mDs2LFch0OIRVAiIMRMer0eK1aswJIlS7B7927s2LEDYrGY67AIGTGaIyDETJs2bcLs2bORkZGBmTNn4k9/+hPXIRFiEfREQAghTo6eCAghxMlRIiCEECdHiYAQQpwcJQJCCHFylAgIIcTJUSIghBAnR4mAEEKc3P8DjpYWivjGnKwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(l,u,200)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set(xlabel='x', \n",
    "       ylabel='|f(x)-g(x)|', \n",
    "       xlim=[0, 1], \n",
    "       ylim=[0, 1.6])\n",
    "plt.plot(x, (f-g).abs()(x))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Demonstrate Pythagorean Theorem\n",
    "\n",
    "Again, define the functions $f(x) = 2x^2-1$ and $g(x)= 4x^3-3x$, both over the domain $[-1,1]$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "l,u = -1, 1\n",
    "f = function(lambda x: 2 * x**2 - 1)\n",
    "g = function(lambda x: 4 * x**3 - 3*x)\n",
    "\n",
    "ifsq = (f**2).integral(l,u)\n",
    "igsq = (g**2).integral(l,u)\n",
    "ifplusgsq = ((f+g)**2).integral(l,u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "∫f²(x)dx = 0.9333\n",
      "∫g²(x)dx = 0.9714\n",
      "∫(f+g)²(x)dx = 1.9048\n"
     ]
    }
   ],
   "source": [
    "print(f'∫f²(x)dx = {ifsq:.4f}')\n",
    "print(f'∫g²(x)dx = {igsq:.4f}')\n",
    "print(f'∫(f+g)²(x)dx = {ifplusgsq:.4f}')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3.9.7 ('base')",
   "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}