HYUNSEONG-KIM/Gaussian Hypergeometric function for negative real values in python.ipynb

## Gaussian Hypergeometric function for negative real values in python.ipynb
{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ${}_2F_1$, Hypergeometric function extension for negative inputs\n",
    "\n",
    "Gauss Hypergeometric function ${}_2F_1(a, b, c; z)$ is defined on $|z|<1$ region in real value, but for complex domain with an analytic continuation, we can extend them to larger domain of complex plane.\n",
    "\n",
    "There are many implementation routines on various system and language of ${}_2F_1$ function.\n",
    "\n",
    "* scipy: python, opensource\n",
    "* Maxima: C, opensource\n",
    "* SageMath: using Maxima, python-sage modified, opensource\n",
    "* mpmath: python, opensource\n",
    "* GSL: C, opensource\n",
    "* NAG: C, C++, commerical\n",
    "* Numerical Recipes: C, C++, Fortran, licensed\n",
    "* Matlab: matlab, commerical\n",
    "* Mathematica: Wolfram Language, commerical\n",
    "* Maple: Maple, commerical\n",
    "\n",
    "**Supporting analytic continuation**\n",
    "\n",
    "* mpmath\n",
    "* Maxima\n",
    "* NAG\n",
    "* Matlab\n",
    "* Mathematica\n",
    "* Maple\n",
    "\n",
    "But scipy does not support analytic continuation feature for ${}_2F_1$ function and does GSL. \n",
    "Simple way is using a mapmath pr Maxima library on project but without additional library simple extended version can be used.\n",
    "\n",
    "This notebook presents two features.\n",
    "\n",
    "1. Extending scipy-routine for negative input with analytic continuation.\n",
    "2. Wolframe Engine usage in Python for special function and compare with 1. result."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extending Scipy implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import hyp2f1, gamma\n",
    "from math import isclose, fabs\n",
    "from sys import float_info\n",
    "\n",
    "EPS = float_info.epsilon*1000"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `scipy.special.hyp2f1` implementation defined on $|z|<1$. However, with analytic continuation we can expand it to $|z|>1$ region of complex plane.\n",
    "\n",
    "See details in \n",
    "\n",
    "> Pearson, J.W., 2009. Computation of hypergeometric functions (Doctoral dissertation, University of Oxford).\n",
    "\n",
    "The above dissertation result was applied to [NAG C library special function $2F1$](https://www.nag.com/numeric/cl/nagdoc_latest/html/s/s22bec.html) implementation. \n",
    "\n",
    "The opensource implementations are `scipy.special.hyp2f1` from scipy and `gsl_sf_hyperg_2F1` from GNU Scientific library.\n",
    "However, both of them only support $|z|<1$ region.\n",
    "\n",
    "In this section, we will not accurately see a detatils and error analysis of the analytical continuation but for practical implmentation to use in application.\n",
    "Below code used `scipy`'s one but its structure also can be applied to GSL implementation too.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hyper2F1(a, b, c, z):\n",
    "        if isclose(z, 1.0, abs_tol=EPS):\n",
    "            if c>a+b:\n",
    "                return gamma(c) * gamma(c-(a+b))/ (gamma(c-a) * gamma(c-b))\n",
    "        # coefficient integer test\n",
    "        ab_abs = fabs(a-b)\n",
    "        ab_int = int(ab_abs)\n",
    "        ab_test = isclose(ab_abs, ab_int, abs_tol =  EPS)\n",
    "\n",
    "        cab_abs = fabs(c-a-b)\n",
    "        cab_int = int(cab_abs)\n",
    "        cab_test = isclose(cab_abs, cab_int, abs_tol =  EPS)\n",
    "        \n",
    "        if not ab_test and not cab_test:\n",
    "            a1 = a2 = 0\n",
    "            b1 = b2 = 0\n",
    "            c1 = c2 = 0\n",
    "            w1 = w2 = 0\n",
    "            co1 = co2 = 0\n",
    "\n",
    "            if z< -1: #case 1\n",
    "                w = 1/(1-z)\n",
    "\n",
    "                co1 = (w**a) * (gamma(c)*gamma(b-a)/(gamma(b) * gamma(c-a)))\n",
    "                co2 = (w**b) *(gamma(c) * gamma(a-b)/(gamma(a) * gamma(c-b)))\n",
    "\n",
    "                a1 = a\n",
    "                b1 = c-b\n",
    "                c1 = a-b+1\n",
    "                a2 = b\n",
    "                b2 = c-a\n",
    "                c2 = b-a+1\n",
    "\n",
    "                w1 = w**a\n",
    "                w2 = w**b\n",
    "\n",
    "            elif z<0: #case 2\n",
    "                w = z/(z-1)\n",
    "\n",
    "                co1 = 1\n",
    "\n",
    "                a1 = a\n",
    "                b1 = c-b\n",
    "                c1 = c\n",
    "\n",
    "                w1 = (1-z)**(-a)\n",
    "\n",
    "            elif z<=1/2: #case 3\n",
    "                w  = z\n",
    "\n",
    "                co1 = 1\n",
    "\n",
    "                a1 = a\n",
    "                b1 = b\n",
    "                c1 = c\n",
    "\n",
    "                w1 = 1\n",
    "\n",
    "            elif z<=1: #case 4\n",
    "                w = 1-z \n",
    "\n",
    "                co1 = (gamma(c) * gamma(c-a-b))/(gamma(c-a) * gamma(c-b))\n",
    "                co2 = (gamma(c) * gamma(a+b-c))/(gamma(a) * gamma(b))\n",
    "\n",
    "                a1 = a\n",
    "                b1 = b\n",
    "                c1 = a+b-c+1\n",
    "                a2 = c-a\n",
    "                b2 = c-b\n",
    "                c2 = c-a-b+1\n",
    "\n",
    "                w1 = 1\n",
    "                w2 = w**(c-a-b)\n",
    "\n",
    "            elif z<=2: #case 5\n",
    "                w =1 - 1/z\n",
    "\n",
    "                co1 = (gamma(c) * gamma(c-a-b))/(gamma(c-a) * gamma(c-b))\n",
    "                co2 = (gamma(c) * gamma(a+b-c))/(gamma(a) * gamma(b))\n",
    "                \n",
    "                a1=a\n",
    "                b1 = a-c+1\n",
    "                c1 = a+b-c+1\n",
    "                a2=c-a\n",
    "                b2=1-a\n",
    "                c2=c-a-b+1\n",
    "\n",
    "                w1 = z**(-a)\n",
    "                w2 = z**(a-c)*(1-z)**(c-a-b)\n",
    "            \n",
    "            else: #case 6\n",
    "                w =1/z\n",
    "\n",
    "                co1 = (gamma(c) * gamma(b-a))/(gamma(b) * gamma(c-a))\n",
    "                co2 = (gamma(c) * gamma(a-b))/(gamma(a) * gamma(c-b))\n",
    "                \n",
    "                a1 =a\n",
    "                b1 = a-c+1\n",
    "                c1 = a-b+1\n",
    "                a2 = b-c+1\n",
    "                b2= b\n",
    "                c2 = b-a+1\n",
    "\n",
    "                w1 = (-z)**(-a)\n",
    "                w2 = (-z)**(-b)\n",
    "                \n",
    "            result = w1*co1*hyp2f1(a1, b1, c1, w) + 0 if co2 ==0 else w2*co2*hyp2f1(a2,b2,c2,w)\n",
    "        else:\n",
    "            result =  (1-z)**(-a) * hyp2f1(a, c-b, c ,z/(z-1))\n",
    "            # Error can be occured\n",
    "            #raise NotImplementedError(\"Not implemented\")\n",
    "        return result "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.010886764704915206"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hyper2F1(2, 4, 6, -15) # extended scipy function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Error\n",
    "a1, b1, c1, z1 =(10, 5, -300.5, 0.5)\n",
    "a2, b2, c2, z2 =(2.25, 3.75, -0.5, -1)\n",
    "cor1, cor2 = -3.852027081523919 * 1E32, -0.631220676949703"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6.518494497351803e+156"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hyper2F1(a1, b1, c1, z1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hyper2F1(a2, b2, c2, z2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Wolfram Engine from python\n",
    "\n",
    "This uses [Wolfram client library for python](https://reference.wolfram.com/language/WolframClientForPython/) and [Wolfram Engine](https://www.wolfram.com/engine/).\n",
    "Client libray is open sourced and Wolfram engine is opend freely for student, individual researcher and prototyper.\n",
    "If you have institution license of Mathematica, you already installed Wolfram Engine.\n",
    "\n",
    "You can use special function implemenations in Wolframe Language on python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "from wolframclient.evaluation import WolframLanguageSession\n",
    "from wolframclient.language import wl, wlexpr\n",
    "session = WolframLanguageSession()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "hyper2F1_wol = session.function(wlexpr('Hypergeometric2F1'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3.8520270815239185e+32"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hyper2F1_wol(a1, b1, c1 , z1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.5339784518435058e+108 \n",
      " -2242928664381.3975\n"
     ]
    }
   ],
   "source": [
    "print(hyp2f1(a1+1.2, b1+1.5, c1+2 , 0.4), \"\\n\", hyper2F1_wol(a1+1.2, b1+1.5, c1+2 , 0.4))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Test cases from Pearson, real valued\n",
    "\n",
    "import math\n",
    "\n",
    "test_cases = [\n",
    "    (0.1, 0.2, 0.3, 0.5),\n",
    "    (-0.1, 0.2, 0.3, 0.5),\n",
    "    (1E-8, 1E-8, 1E-8, 1E-6),\n",
    "    (2+1E-9, 3, 5, -0.75),\n",
    "    (-2, -3, -5 + 1E-9, 0.5),\n",
    "    (-1, -1.5, -2-1E-15, 0.5),\n",
    "    (500, -500, 500, 0.75),\n",
    "    (500, 500, 500, -0.6),\n",
    "    (-1000, -2000, -4000.1, -0.5),\n",
    "    (-100, -200, -300+1E-9, 0.5*math.sqrt(2)),\n",
    "    (300, 10, 5, 0.5),\n",
    "    (5, -300, 10, 0.5),\n",
    "    (10, 5, -300.5, 0.5),\n",
    "    (2.25, 3.75, -0.5, -1),\n",
    "]\n",
    "\n",
    "correct = [\n",
    "    1.046432811217352,\n",
    "    0.956434210968214,\n",
    "    1.000000000000010,\n",
    "    0.492238858852651,\n",
    "    0.474999999913750,\n",
    "    0.625000000000000,\n",
    "    9.332636185032189 * 1E-302,\n",
    "    8.709809816217217 * 1E-103,\n",
    "    5.233580403196932 * 1E94,\n",
    "    2.653635302903707 * 1E-31,\n",
    "    3.912238919961547 * 1E98,\n",
    "    1.661006238211309 * 1E-7,\n",
    "    -3.852027081523919 * 1E32,\n",
    "    -0.631220676949703\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------------------------------------------------\n",
      "Correct: 1.046432811217352\n",
      "Scipy ext\n",
      "result:1.046432811217352\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "Wolfram\n",
      "result:1.046432811217352\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "-------------------------------------------------------\n",
      "Correct: 0.956434210968214\n",
      "Scipy ext\n",
      "result:0.9564342109682141\n",
      "error_abs:1.1102230246251565e-16\n",
      "error_rel:1.160793928001864e-16\n",
      "Wolfram\n",
      "result:0.9564342109682142\n",
      "error_abs:2.220446049250313e-16\n",
      "error_rel:2.321587856003728e-16\n",
      "-------------------------------------------------------\n",
      "Correct: 1.00000000000001\n",
      "Scipy ext\n",
      "result:1.00000000000001\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "Wolfram\n",
      "result:1.00000000000001\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "-------------------------------------------------------\n",
      "Correct: 0.492238858852651\n",
      "Scipy ext\n",
      "result:0.492238858852651\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "Wolfram\n",
      "result:0.492238858852651\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "-------------------------------------------------------\n",
      "Correct: 0.47499999991375\n",
      "Scipy ext\n",
      "result:0.47499999991375\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "Wolfram\n",
      "result:0.47499999991374997\n",
      "error_abs:5.551115123125783e-17\n",
      "error_rel:1.168655815607105e-16\n",
      "-------------------------------------------------------\n",
      "Correct: 0.625\n",
      "Scipy ext\n",
      "result:0.6250000000000002\n",
      "error_abs:2.220446049250313e-16\n",
      "error_rel:3.552713678800501e-16\n",
      "Wolfram\n",
      "result:0.6250000000000001\n",
      "error_abs:1.1102230246251565e-16\n",
      "error_rel:1.7763568394002506e-16\n",
      "-------------------------------------------------------\n",
      "Correct: 9.332636185032189e-302\n",
      "Scipy ext\n",
      "result:0.0\n",
      "error_abs:9.332636185032189e-302\n",
      "error_rel:1.0\n",
      "Wolfram\n",
      "result:9.332636185032189e-302\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "-------------------------------------------------------\n",
      "Correct: 8.709809816217216e-103\n",
      "Scipy ext\n",
      "result:8.709809816216975e-103\n",
      "error_abs:2.409061978389464e-116\n",
      "error_rel:2.7659180042070667e-14\n",
      "Wolfram\n",
      "result:8.709809816216975e-103\n",
      "error_abs:2.409061978389464e-116\n",
      "error_rel:2.7659180042070667e-14\n",
      "-------------------------------------------------------\n",
      "Correct: 5.233580403196932e+94\n",
      "Scipy ext\n",
      "result:5.23358040319729e+94\n",
      "error_abs:3.5867757562151065e+81\n",
      "error_rel:6.853388082132311e-14\n",
      "Wolfram\n",
      "result:5.2335804031969554e+94\n",
      "error_abs:2.3714219875802357e+80\n",
      "error_rel:4.531165674137065e-15\n",
      "-------------------------------------------------------\n",
      "Correct: 2.653635302903707e-31\n",
      "Scipy ext\n",
      "result:2.6536353029036913e-31\n",
      "error_abs:1.5764607723654192e-45\n",
      "error_rel:5.9407589680480845e-15\n",
      "Wolfram\n",
      "result:2.6536352681765626e-31\n",
      "error_abs:3.4727144509184485e-39\n",
      "error_rel:1.308663043153848e-08\n",
      "-------------------------------------------------------\n",
      "Correct: 3.912238919961547e+98\n",
      "Scipy ext\n",
      "result:3.912238919961548e+98\n",
      "error_abs:1.2141680576410807e+83\n",
      "error_rel:3.103512036154005e-16\n",
      "Wolfram\n",
      "result:3.912238919961547e+98\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "-------------------------------------------------------\n",
      "Correct: 1.661006238211309e-07\n",
      "Scipy ext\n",
      "result:1.661006238211307e-07\n",
      "error_abs:2.117582368135751e-22\n",
      "error_rel:1.2748792385127438e-15\n",
      "Wolfram\n",
      "result:1.661006238211309e-07\n",
      "error_abs:0.0\n",
      "error_rel:0.0\n",
      "-------------------------------------------------------\n",
      "Correct: -3.852027081523919e+32\n",
      "Scipy ext\n",
      "result:-6.518494497351803e+156\n",
      "error_abs:6.518494497351803e+156\n",
      "error_rel:1.6922244728282102e+124\n",
      "Wolfram\n",
      "result:-3.8520270815239185e+32\n",
      "error_abs:7.205759403792794e+16\n",
      "error_rel:1.8706408992696095e-16\n",
      "-------------------------------------------------------\n",
      "Correct: -0.631220676949703\n",
      "Scipy ext\n",
      "result:-0.6312206769497023\n",
      "error_abs:6.661338147750939e-16\n",
      "error_rel:1.0553105104131005e-15\n",
      "Wolfram\n",
      "result:-0.6312206769497027\n",
      "error_abs:2.220446049250313e-16\n",
      "error_rel:3.517701701377002e-16\n"
     ]
    }
   ],
   "source": [
    "result = []\n",
    "for param, cor in zip(test_cases, correct):\n",
    "    abs_cor = math.fabs(cor)\n",
    "    element = []\n",
    "    a = param[0]\n",
    "    b = param[1]\n",
    "    c = param[2]\n",
    "    z = param[3]\n",
    "\n",
    "    print(\"-------------------------------------------------------\")\n",
    "    print(f\"Correct: {cor}\")\n",
    "    print(\"Scipy ext\")\n",
    "    try:\n",
    "        r_scipy = hyper2F1(a, b, c, z)\n",
    "\n",
    "        r_error = math.fabs(cor-r_scipy)\n",
    "        \n",
    "\n",
    "        \n",
    "        print(f\"result:{r_scipy}\")\n",
    "        print(f\"error_abs:{r_error}\")\n",
    "        print(f\"error_rel:{r_error/abs_cor}\")\n",
    "\n",
    "        element.append((r_scipy, r_error, r_error/abs_cor))\n",
    "\n",
    "    except Exception as e:\n",
    "        print(\"Error\")\n",
    "    print(\"Wolfram\")\n",
    "    try:\n",
    "        r_wol = hyper2F1_wol(a, b, c, z)\n",
    "\n",
    "        r_error_wol = math.fabs(cor-r_wol)\n",
    "\n",
    "        \n",
    "        print(f\"result:{r_wol}\")\n",
    "        print(f\"error_abs:{r_error_wol}\")\n",
    "        print(f\"error_rel:{r_error_wol/abs_cor}\")\n",
    "\n",
    "        element.append((r_wol, r_error_wol, r_error_wol/abs_cor))\n",
    "    except Exception as e:\n",
    "        print(\"Error\")\n",
    "    \n",
    "    result.append(element)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "scipy_rel_error = [x[0][2] for x in result]\n",
    "wol_rel_error = [x[1][2] for x in result]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\qwqwh\\AppData\\Local\\Temp\\ipykernel_18372\\575945362.py:8: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.\n",
      "  axs[0].set_ylim( scipy_rel_error[12]-y_EPS, scipy_rel_error[12]+y_EPS)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGsCAYAAAD+L/ysAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4W0lEQVR4nO3df3xT9d3//2da0tBqW1YQSrCUCsVOqFgUGCgCk6LoEOYQlclPN3dN5IdMB0xFEKUDBwLjAsU5OudV5OOlFLc5aBVsleFEsII/BhQLitCr8lWblkob2vP9wzWjJv0J6XnTPu63W28277xzzuu8TJon75wkDsuyLAEAABgsxO4CAAAA6kNgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGa3GBJTc3V6NGjZLb7ZbD4VBmZmajbn/q1ClNnjxZycnJatOmjcaMGeM35+WXX1ZqaqouuugiRUVFaeDAgdq6dWut23zhhRfkcDgCbgsAANSvxQWWkydPqk+fPlq9enWTbl9ZWanw8HDNmDFDw4cPDzgnNzdXqampevXVV7V7924NGzZMo0aN0nvvvec398iRI7r//vs1ePDgJtUDAAAkR0v+8kOHw6FNmzbVWNmoqKjQQw89pP/5n//R119/rd69e2vJkiUaOnSo3+0nT56sr7/+ukGrNL169dJtt92m+fPn+8YqKys1ZMgQTZkyRW+++WaDtwUAAGpqcSss9ZkyZYp27NihF154QXv37tWtt96qG264QQcPHmzyNquqqlRSUqKYmJga448++qguuugi3XXXXWdbNgAArVobuwtoTocOHdKGDRt09OhRud1uSdL999+vLVu2aP369Vq8eHGTtrts2TKdPHlS48aN843t2LFDzz77rPLy8s5F6QAAtGqtKrDs2bNHlmWpZ8+eNcbLy8vVvn37Jm1zw4YNWrBggTZv3qyOHTtKkkpKSnTnnXfqmWeeUYcOHc66bgAAWrtWFViqqqoUGhqq3bt3KzQ0tMZ1F154YaO3t3HjRt1111168cUXa5yge+jQIR0+fFijRo2qsW9JatOmjfbv36/u3bs38SgAAGh9WlVgSUlJUWVlpYqKis76XTsbNmzQ1KlTtWHDBt100001rktKStK+fftqjD300EMqKSnRypUrFRcXd1b7BgCgtWlxgaW0tFT5+fm+ywUFBcrLy1NMTIx69uypn/70p5o4caKWLVumlJQUnThxQtu2bVNycrJuvPFGSdJHH32kiooKffnllyopKfGdh3LFFVdI+jasTJw4UStXrtQPfvADFRYWSpLCw8MVHR2ttm3bqnfv3jXqateunST5jQMAgAawWpjt27dbkvx+Jk2aZFmWZVVUVFjz58+3unXrZjmdTis2Ntb68Y9/bO3du9e3jfj4+IDbqDZkyJA69xHIpEmTrNGjR5/VsZ06dcp65JFHrFOnTp3VdloSehIYffFHT/zRk8Doiz8TetKiP4elpfF4PIqOjlZxcbGioqLsLscI9CQw+uKPnvijJ4HRF38m9KTVfQ4LAAA4/xBYAACA8VrMSbdVVVU6duyYIiMj5XA47C4nKDweT43/gp7Uhr74oyf+6Elg9MVfMHtiWZZKSkrkdrsVElL7OkqLOYfl6NGjvF0YAIDz1GeffaaLL7641utbzApLZGSkpG8PuKWeJOX1epWVlaURI0bI6XTaXY4R6Elg9MUfPfFHTwKjL/6C2ROPx6O4uDjf83htWkxgqX4ZKCoqqkUHloiICEVFRfEg+jd6Ehh98UdP/NGTwOiLv+boSX2nc3DSLQAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACM1+jAkpubq1GjRsntdsvhcCgzM7PO+ZMnT5bD4fD76dWrl29Oenp6wDmnTp1q9AEBAICWp9GB5eTJk+rTp49Wr17doPkrV67U8ePHfT+fffaZYmJidOutt9aYFxUVVWPe8ePH1bZt28aWBwAAWqBGf1vzyJEjNXLkyAbPj46OVnR0tO9yZmamvvrqK02ZMqXGPIfDodjY2MaWAwAAWoFGB5az9eyzz2r48OGKj4+vMV5aWqr4+HhVVlbqiiuu0KJFi5SSklLrdsrLy1VeXu677PF4JH37Fdherzc4xdus+rha6vE1BT0JjL74oyf+6Elg9MVfMHvS0G06LMuymroTh8OhTZs2acyYMQ2af/z4ccXFxSkjI0Pjxo3zjb/99tvKz89XcnKyPB6PVq5cqVdffVXvv/++EhMTA25rwYIFWrhwod94RkaGIiIimnQ8AACgeZWVlWn8+PEqLi5WVFRUrfOaNbCkpaVp2bJlOnbsmMLCwmqdV1VVpb59++raa6/VqlWrAs4JtMISFxenEydO1HnA5zOv16vs7GylpqbK6XTaXY4R6Elg9MUfPfFHTwKjL/6C2ROPx6MOHTrUG1ia7SUhy7L0xz/+URMmTKgzrEhSSEiI+vXrp4MHD9Y6x+VyyeVy+Y07nc4WfwdrDcfYWPQkMPrij574oyeB0Rd/wehJQ7fXbJ/DkpOTo/z8fN111131zrUsS3l5eercuXMzVAYAAEzX6BWW0tJS5efn+y4XFBQoLy9PMTEx6tq1q+bNm6fPP/9czz33XI3bPfvssxowYIB69+7tt82FCxfqBz/4gRITE+XxeLRq1Srl5eXpv//7v5twSAAAoKVpdGB59913NWzYMN/l2bNnS5ImTZqk9PR0HT9+XJ9++mmN2xQXF+ull17SypUrA27z66+/1t13363CwkJFR0crJSVFubm56t+/f2PLAwAALVCjA8vQoUNV13m66enpfmPR0dEqKyur9TZPPvmknnzyycaWAgAAWgm+SwgAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGC8RgeW3NxcjRo1Sm63Ww6HQ5mZmXXOnzx5shwOh99Pr169asx76aWXdNlll8nlcumyyy7Tpk2bGlsaAABooRodWE6ePKk+ffpo9erVDZq/cuVKHT9+3Pfz2WefKSYmRrfeeqtvzs6dO3XbbbdpwoQJev/99zVhwgSNGzdO//znPxtbHgAAaIHaNPYGI0eO1MiRIxs8Pzo6WtHR0b7LmZmZ+uqrrzRlyhTf2IoVK5Samqp58+ZJkubNm6ecnBytWLFCGzZsaGyJAACghWl0YDlbzz77rIYPH674+Hjf2M6dO3XffffVmHf99ddrxYoVtW6nvLxc5eXlvssej0eS5PV65fV6z23Rhqg+rpZ6fE1BTwKjL/7oiT96Ehh98RfMnjR0m80aWI4fP66///3vysjIqDFeWFioTp061Rjr1KmTCgsLa91WWlqaFi5c6DeelZWliIiIc1OwobKzs+0uwTj0JDD64o+e+KMngdnZlypLOuRxyOOVopxS9yhLIQ7byvEJRk/KysoaNK9ZA0t6erratWunMWPG+F3ncNT8P2FZlt/YmebNm6fZs2f7Lns8HsXFxWnEiBGKioo6ZzWbxOv1Kjs7W6mpqXI6nXaXYwR6Ehh98UdP/NGTwOzuy9YP/09P/PVDHS097Ru7+MI2mvujXrq+V6c6bhk8wexJ9Ssk9Wm2wGJZlv74xz9qwoQJCgsLq3FdbGys32pKUVGR36rLmVwul1wul9+40+ls8Q+81nCMjUVPAqMv/uiJP3oSmB192fLBcU1/IU/X5e/Sqp0bdekXR7T/onitGXSbppd6tfbOK3VD787NWtOZgtGThm6v2T6HJScnR/n5+brrrrv8rhs4cKDfMlNWVpYGDRrUXOUBAGCryipLj23ep+vyd2ndS4vU99h+XeA9pb7H9mvd/y7SdYd26fHN+1RZZdldqi0aHVhKS0uVl5envLw8SVJBQYHy8vL06aefSvr2pZqJEyf63e7ZZ5/VgAED1Lt3b7/rZs6cqaysLC1ZskT/+te/tGTJEr322muaNWtWY8sDAOC89E7Blzpa4tU9OzcqRDVDSYgs/fIf/0+flXj1TsGXNlVor0YHlnfffVcpKSlKSUmRJM2ePVspKSmaP3++pG9PrK0OL9WKi4v10ksvBVxdkaRBgwbphRde0Pr163X55ZcrPT1dGzdu1IABAxpbHgAA56WiklOSpEu/OBLw+ktPHKkxr7Vp9DksQ4cOlWXVvhyVnp7uNxYdHV3vWcBjx47V2LFjG1sOAAAtQsfItpKk/RfFq++x/X7X7+8QX2Nea8N3CQEAYID+CTG6ONKpNYNuU5Vqvku2Sg6tHTROcZFO9U+IsalCexFYAAAwQGiIQw+NTtbr3fvp52Mf1m53kkrDwrXbnaS7xz6s17v304OjkxVqwgey2KDZP+kWAAAEdkPvzlp755ValOnUT7r3943HRTq1dnSyrW9pthuBBQAAg9zQu7Ou7tFByQuyJEnpU/ppcOJFrXZlpRovCQEAYJgzw0n/hJhWH1YkAgsAADgPEFgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8RodWHJzczVq1Ci53W45HA5lZmbWe5vy8nI9+OCDio+Pl8vlUvfu3fXHP/7Rd316erocDoffz6lTpxpbHgAAaIHaNPYGJ0+eVJ8+fTRlyhT95Cc/adBtxo0bp//7v//Ts88+qx49eqioqEinT5+uMScqKkr79++vMda2bdvGlgcAAFqgRgeWkSNHauTIkQ2ev2XLFuXk5OiTTz5RTEyMJKlbt25+8xwOh2JjYxtbDgAAaAUaHVga65VXXtFVV12lpUuX6s9//rMuuOAC3XzzzVq0aJHCw8N980pLSxUfH6/KykpdccUVWrRokVJSUmrdbnl5ucrLy32XPR6PJMnr9crr9QbvgGxUfVwt9fiagp4ERl/80RN/9CQwE/ri9Z4+43evvA7Ltlqqazjzv8HYdn2CHlg++eQTvfXWW2rbtq02bdqkEydO6J577tGXX37pO48lKSlJ6enpSk5Olsfj0cqVK3X11Vfr/fffV2JiYsDtpqWlaeHChX7jWVlZioiICOox2S07O9vuEoxDTwKjL/7oiT96EpidfSmvlKqforduzZIr1LZSaghGT8rKyho0z2FZVpNjm8Ph0KZNmzRmzJha54wYMUJvvvmmCgsLFR0dLUl6+eWXNXbsWJ08ebLGKku1qqoq9e3bV9dee61WrVoVcLuBVlji4uJ04sQJRUVFNfWQjOb1epWdna3U1FQ5nU67yzECPQmMvvijJ/7oSWAm9KWs4rT6LNomSXr/4R8qIizo6wt1CmZPPB6POnTooOLi4jqfv4Pegc6dO6tLly6+sCJJ3//+92VZlo4ePRpwBSUkJET9+vXTwYMHa92uy+WSy+XyG3c6nS3+gdcajrGx6Elg9MUfPfFHTwKzsy9Oy/GdOuwNLNWC0ZOGbi/on8Ny9dVX69ixYyotLfWNHThwQCEhIbr44osD3sayLOXl5alz587BLg8AAJwHGh1YSktLlZeXp7y8PElSQUGB8vLy9Omnn0qS5s2bp4kTJ/rmjx8/Xu3bt9eUKVP00UcfKTc3Vw888ICmTp3qezlo4cKF2rp1qz755BPl5eXprrvuUl5env7rv/7rHBwiAAA43zV6jendd9/VsGHDfJdnz54tSZo0aZLS09N1/PhxX3iRpAsvvFDZ2dmaPn26rrrqKrVv317jxo3TY4895pvz9ddf6+677/ad55KSkqLc3Fz179//bI4NAAC0EI0OLEOHDlVd5+mmp6f7jSUlJdV5ZvGTTz6pJ598srGlAACAVoLvEgIAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4LS6w7NixQ6NGjZLb7ZbD4VBmZmZQ95ebm1vn/rxer+bMmaPk5GRdcMEFcrvdmjhxoo4dO3ZW+505c6auvPJKuVwuXXHFFQHn7Nu3T0OGDFF4eLi6dOmiRx99VJZlndV+63Lq1ClNnjxZycnJatOmjcaMGRO0fQEAWpcWF1jKysrUp08frV69uln2d/LkyTr3V1ZWpj179ujhhx/Wnj179PLLL+vAgQO6+eab69yuw+HQ4cOHa73esixNnTpVt912W8DrPR6PUlNT5Xa7tWvXLv3+97/X7373Oy1fvrzBx9ZYlZWVCg8P14wZMzR8+PCg7QcA0ApZLURxcbElySouLvaNSbI2bdpUY155ebn1wAMPWG6324qIiLD69+9vbd++/ZzUEGh/gbzzzjuWJOvIkSN1bqugoKDebT3yyCNWnz59/MbXrFljRUdHW6dOnfKNpaWlWW6326qqqvKNvfLKK1bfvn0tl8tlJSQkWAsWLLC8Xm+9+63PpEmTrNGjR5/1dgAALVug5+9AWswKS2RkpIqLixUZGVnnvClTpmjHjh164YUXtHfvXt1666264YYbdPDgwWaqVCouLpbD4VC7du2Cto+dO3dqyJAhcrlcvrHrr79ex44d863cbN26VXfeeadmzJihjz76SE8//bTS09P1+OOPB60uAADO1NDn7xYTWBwOh6KiouRwOGqdc+jQIW3YsEEvvviiBg8erO7du+v+++/XNddco/Xr1zdLnadOndLcuXM1fvx4RUVFBW0/hYWF6tSpU42x6suFhYWSpMcff1xz587VpEmTdMkllyg1NVWLFi3S008/HbS6AAA4U0Oev6UWFFgaYs+ePbIsSz179tSFF17o+8nJydGhQ4ckSYcPH5bD4ajz5957723S/r1er26//XZVVVVpzZo1Na4bOXJkjZokqVevXn5jjfHd//nWv0+4rR7fvXu3Hn300Rr7+PnPf67jx4+rrKxMkjR06NA6e9GUugAAaKw2dhfQnKqqqhQaGqrdu3crNDS0xnXVT7xdunTRxx9/XOd2vve97zV6316vV+PGjVNBQYG2bdvmt7ryhz/8Qd98843vcmJiol599VV16dKl0fuSpNjYWN9KSrWioiJJ/1lpqaqq0sKFC3XLLbf43b5t27aSpOeee84XXgIJCWlVmRcAYJNWFVhSUlJUWVmpoqIiDR48OOAcp9OppKSkc7rf6rBy8OBBbd++Xe3bt/ebEyiYxMfHq1u3bk3a58CBA/Wb3/xGFRUVCgsLkyRlZWXJ7Xb7ttm3b1/t379fPXr0qHU7Xbt2bdL+AQA4l1pcYCktLVV+fr7vckFBgfLy8hQTE6OePXvqpz/9qSZOnKhly5YpJSVFJ06c0LZt25ScnKwbb7zxnO6va9euOn36tMaOHas9e/bor3/9qyorK30rHzExMb4w0Vj5+fkqLS1VYWGhvvnmG+Xl5UmSLrvsMoWFhWn8+PFauHChJk+erN/85jc6ePCgFi9erPnz5/teEpo/f75+9KMfKS4uTrfeeqtCQkK0d+9e7du3T4899liT6vroo49UUVGhL7/8UiUlJb66avusGAAAGqRZ3rPUjLZv325J8vuZNGmSZVmWVVFRYc2fP9/q1q2b5XQ6rdjYWOvHP/6xtXfv3qDsr6CgIOD1kup8O7XqeVvzkCFDAm7zzNvs3bvXGjx4sOVyuazY2FhrwYIFNd7SbFmWtWXLFmvQoEFWeHi4FRUVZfXv399at25dk3phWZYVHx8fsC4AAM6Gw7KC+NGnAAAA5wBnTAIAAOMRWAAAgPFaTGCxLEsejyeoX+4HAADOrYY+f7eYwFJSUqLo6GiVlJTYXUrQeL1ebd68WV6v1+5SjEFPAqMv/uiJP3oSmAl9Kas4rW5z/6Zuc/+msorTttVRLZg9aejzd4sJLAAAoOUisAAAAOMRWAAAgPEILAAAwHgEFgAAYLwW911CAACg+VRWVtb57qHQ0FC1adPG9z12TUVgAQAATVJaWqqjR4/W+xkqERER6ty5c5O/8FcisAAAgCaorKzU0aNHFRERoYsuuijgCoplWaqoqNAXX3yhgoICJSYmKiSkaWejEFgAAECjeb1eWZaliy66SOHh4bXOCw8Pl9Pp1JEjR1RRUaG2bds2aX+cdAsAAJqsIeemNHVVpcY2znoLAAAAQUZgAQAAxiOwAAAA4xFYAACA8QgsAACgyer7DJaGzqkPgQUAADRaaGioJKmioqLeuWVlZZIkp9PZ5P3xOSwAAKDR2rRpo4iICH3xxRdyOp0B37psWZbKyspUVFSkdu3a+UJOk/Z3NsUCAIDWyeFwqHPnziooKNCRI0fqnNuuXTvFxsae1f4ILAAAoEnCwsKUmJhY58tCTqfzrFZWqhFYAABAk4WEhDT54/YbtZ+g7wEAAOAsEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDy+/BBAUFRWWXqn4EsVlZxSx8i26p8Qo9AQh91lAThPNcsKS25urkaNGiW32y2Hw6HMzMwa11uWpQULFsjtdis8PFxDhw7Vhx9+2BylAQiCLR8c15C0bN3xzNua+UKe7njmbQ1Jy9aWD47bXRqA81SzBJaTJ0+qT58+Wr16dcDrly5dquXLl2v16tXatWuXYmNjlZqaqpKSkuYoD8A5tOWD4/rl87uV9N4OvfznX+nD5WP18p9/paS8Hfrl87sJLQCapFleEho5cqRGjhwZ8DrLsrRixQo9+OCDuuWWWyRJf/rTn9SpUydlZGToF7/4RXOUCOAcqKyy9Njmfbouf5fWvbRIIbIkSX2P7de6/12ku8c+rMc3hyn1slheHgLQKLafw1JQUKDCwkKNGDHCN+ZyuTRkyBD94x//qDWwlJeXq7y83HfZ4/FIkrxer7xeb3CLtkn1cbXU42sKehKYXX35Z8GXOlri1aqdG31hpVqILP3yH/9PP+neXzvzizQgIaZZa+O+4o+eBGZCX7ze02f87pXXYdUxO/iC2ZOGbtP2wFJYWChJ6tSpU43xTp066ciRI7XeLi0tTQsXLvQbz8rKUkRExLkt0jDZ2dl2l2AcehJYc/dl9wmHpFBd+kXgx+6lJ74dz3rzn/r/PrbnDzD3FX/0JDA7+1JeKVU/RW/dmiVXqG2l1BCMnpSVlTVonu2BpZrDUXN52LIsv7EzzZs3T7Nnz/Zd9ng8iouL04gRIxQVFRW0Ou3k9XqVnZ2t1NRUOZ1Ou8sxAj0JzK6+tC/4Us8dfFf7L4pX32P7/a7f3yFekjRi8ABbVli4r9RETwIzoS9lFaf163e2SZKuv36EIsLsfboOZk+qXyGpj+2BJTY2VtK3Ky2dO3f2jRcVFfmtupzJ5XLJ5XL5jTudzhb/wGsNx9hY9CSw5u7LwB4ddXGkU2sG3aZ1/7uoxstCVXJo7aBxiot0amCPjradw8J9xR89CczOvjit/zw+vq3D9qdrScHpSUO3Z/sHxyUkJCg2NrbGMlNFRYVycnI0aNAgGysD0FihIQ49NDpZr3fvp5+PfVi73UkqDQvXbneS7h77sF7v3k8Pjk7mhFsAjdYska20tFT5+fm+ywUFBcrLy1NMTIy6du2qWbNmafHixUpMTFRiYqIWL16siIgIjR8/vjnKA3AO3dC7s9beeaUWZTr1k+79feNxkU6tHZ2sG3p3ruPWABBYswSWd999V8OGDfNdrj73ZNKkSUpPT9evf/1rffPNN7rnnnv01VdfacCAAcrKylJkZGRzlAfgHLuhd2dd3aODkhdkSZLSp/TT4MSLWFkB0GTNEliGDh0qy6r9HQEOh0MLFizQggULmqMcAM3gzHDCx/IDOFu2n8MCAABQHwILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8YwJLSUmJZs2apfj4eIWHh2vQoEHatWuX3WUBAAADGBNYfvaznyk7O1t//vOftW/fPo0YMULDhw/X559/bndpAADAZkYElm+++UYvvfSSli5dqmuvvVY9evTQggULlJCQoLVr19pdHgAAsFkbuwuQpNOnT6uyslJt27atMR4eHq633nor4G3Ky8tVXl7uu+zxeCRJXq9XXq83eMXaqPq4WurxNQU9CcyEvni9p8/43Suvw7Ktluoazvwv6EltTOhLa3r8NHSbDsuy7O3Cvw0aNEhhYWHKyMhQp06dtGHDBk2cOFGJiYnav3+/3/wFCxZo4cKFfuMZGRmKiIhojpIB1KG8Uvr1O9/+m2hp/9NyhdpcEHAeaU2Pn7KyMo0fP17FxcWKioqqdZ4xgeXQoUOaOnWqcnNzFRoaqr59+6pnz57as2ePPvroI7/5gVZY4uLidOLEiToP+Hzm9XqVnZ2t1NRUOZ1Ou8sxAj0JzIS+lFWcVp9F2yRJ7z/8Q0WE2buga0JPTENPAjOhL63p8ePxeNShQ4d6A4sRLwlJUvfu3ZWTk6OTJ0/K4/Goc+fOuu2225SQkBBwvsvlksvl8ht3Op0t/oHXGo6xsehJYHb2xWk5vlOHGX9uuK/4oyeB8fjxF4yeNHR7Rpx0e6YLLrhAnTt31ldffaWtW7dq9OjRdpcEAABsZkZkk7R161ZZlqVLL71U+fn5euCBB3TppZdqypQpdpcGAABsZswKS3FxsaZNm6akpCRNnDhR11xzjbKyslimBAAA5qywjBs3TuPGjbO7DAAAYCBjVlgAAABqQ2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYr8UFlh07dmjUqFFyu91yOBzKzMwM6v5Onz6thx56SAkJCQoPD9cll1yiRx99VFVVVUHdLwAArUkbuws418rKytSnTx9NmTJFP/nJT4K+vyVLluipp57Sn/70J/Xq1UvvvvuupkyZoujoaM2cOTPo+wcAoDVocSssqampeuyxx3TLLbcEvL6iokK//vWv1aVLF11wwQUaMGCA3njjjSbvb+fOnRo9erRuuukmdevWTWPHjtWIESP07rvvNnmbtamqqtJ7773H6s0Z6ElgJvQlIqyNDv/2Jh3+7U2KCLP/30Ym9MQ09CQwE/rC4ycAq4WoqqqyiouLraqqKt+YJGvTpk015o0fP94aNGiQlZuba+Xn51tPPPGE5XK5rAMHDjRpv2lpaVZ8fLy1f/9+y7IsKy8vz+rYsaOVkZHR5GOpTXFxsSXJKi4uPufbPl/Rk8Doiz964o+eBEZf/AWzJ4GevwOxP7adIw6HQ1FRUXXOOXTokDZs2KCjR4/K7XZLku6//35t2bJF69ev1+LFixu93zlz5qi4uFhJSUkKDQ1VZWWlHn/8cd1xxx1NOg4AAFqThjx/Sy3wHJa67NmzR5ZlqWfPnjXGy8vL1b59e0nS4cOHlZCQUOd2pk2bptWrV0uSNm7cqOeff14ZGRnq1auX8vLyNGvWLLndbk2aNCk4BwIAQCvTqgJLVVWVQkNDtXv3boWGhta47sILL5QkdenSRR9//HGd2/ne977n+/2BBx7Q3Llzdfvtt0uSkpOTdeTIEaWlpRFYAAA4R1pVYElJSVFlZaWKioo0ePDggHOcTqeSkpIavM2ysjKFhNQ8dzk0NDQoJya5XC498sgjcrlc53zb5yt6Ehh98UdP/NGTwOiLPxN64rAsy7Jt70FQWlqq/Px8Sd8GlOXLl2vYsGGKiYlR165ddeedd2rHjh1atmyZUlJSdOLECW3btk3Jycm68cYbG72/yZMn67XXXtPTTz+tXr166b333tPdd9+tqVOnasmSJef68AAAaJVaXGB54403NGzYML/xSZMmKT09XV6vV4899piee+45ff7552rfvr0GDhyohQsXKjk5udH7Kykp0cMPP6xNmzapqKhIbrdbd9xxh+bPn6+wsLBzcUgAALR6LS6wAACAlqfFfXAcAABoeQgsAADAeC0msFiWJY/HI17hAgDg/NHQ5+8WE1hKSkoUHR2tkpISu0sJGq/Xq82bN8vr9dpdijHoSWAm9KWs4rS6zf2bus39m8oqTttWRzUTemIaehKYCX1pTY+fhj5/t5jAAgAAWi4jAktaWpr69eunyMhIdezYUWPGjNH+/fvtLgsAABjCiMCSk5OjadOm6e2331Z2drZOnz6tESNG6OTJk3aXBgAADGDER/Nv2bKlxuX169erY8eO2r17t6699lqbqgIAAKYwIrB8V3FxsSQpJiam1jnl5eUqLy/3XfZ4PJK+PTGopZ5AVn1cLfX4moKeBGZCX7ze02f87pXXYe87+EzoiWnoSWAm9KU1PX4auk3jPunWsiyNHj1aX331ld58881a5y1YsEALFy70G8/IyFBEREQwSwTQAOWV0q/f+fbfREv7n5YrtJ4bAPBpTY+fsrIyjR8/XsXFxYqKiqp1nnGBZdq0afrb3/6mt956SxdffHGt8wKtsMTFxenEiRN1HvD5zOv1Kjs7W6mpqXI6nXaXYwR6EpgJfSmrOK0+i7ZJkt5/+IeKCLN3QdeEnpiGngRmQl9a0+PH4/GoQ4cO9QYWo14Smj59ul555RXl5ubWGVakb7/qOtDXXDudzhb/wGsNx9hY9CQwO/vitBzfqcOMPzfcV/zRk8B4/PgLRk8auj0jOmBZlqZPn65NmzbpjTfeUEJCgt0lAQAAgxgRWKZNm6aMjAxt3rxZkZGRKiwslCRFR0crPDzc5uoAAIDdjPgclrVr16q4uFhDhw5V586dfT8bN260uzQAAGAAI1ZYDDvvFwAAGMaIFRYAAIC6EFgAAIDxCCwAAMB4BBYAAGA8AgsAtFKrtx/SzJ2hWr39kN2lAPUisABAK7Tq9YNaue2QJIdWbjukVa8ftLskoE4EFgBoZVa9flDLsw/UGFuefYDQAqMRWACgFQkUVqoRWmAyAgsAtBJ1hZVqhBaYisACAK1AQ8JKNUILTERgAYBW4MkGhpWmzgeCjcACAK3Afak9gzofCDYCCwC0AjOuS9TsBoaQ2ak9NeO6xCBXBDQOgQUAWomGhBbCCkxFYAGAVqSu0EJYgckILADQysy4LlHTf9ijxhhhBaYjsABAK/TLod19v98zJIGwAuMRWACglfvFtQl2lwDUi8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeASWeqx6/aAS5v5Nq14/aHcpWr39kGbuDNXq7YfsLsWYvtCTwEzqiyStfcP+OkzqiUn3FVOY1BOT7isSj59qBJY6rHjtgJZnH5AlaXn2Aa147YCttazcdkiSQyu3HbK9FhP6Qk9qr8WEvpz5R/b32/JtfSIypSfVtZhwX6mssny/7zr8VY3Lzc2UnlTXYsJ9hcePP4dlWfbdS79jzZo1euKJJ3T8+HH16tVLK1as0ODBgxt0W4/Ho+joaBUXFysqKuqsa5mesUd/2Xvcb3zU5Z31+/F9z3r71NIy6qCWwFa9flDLs/3/qM1O7akZ1yU2Wx2SOT0xqZYtHxzXo5l7daz0tG/s4kinHhqdrBt6d262OiRzemJSLa3t8dPQ529jVlg2btyoWbNm6cEHH9R7772nwYMHa+TIkfr000+bvZba/gdJ0l/2Htf0jD3UYmMtptRBLYHV9sdW+vZfz835L0VTemJSLVs+OK5fPr9bl77/tm8sY8M8JeXt0C+f360tHwSuMRhM6YlJtfD4qZ0xKywDBgxQ3759tXbtWt/Y97//fY0ZM0ZpaWn13v5crbCseO2AVrz27zuEZclVWRFw3rRr4nXP4G5N3k9DrHnzsP77rSP1zmtNtZhSB7UE9tRbR/T7Nw/XO2/64G76r2vig1aHZE5PTKqlssrS6JU56vHBO3psy2oNuPfPkqRdq+9UW2+FZt/8Kx3q1V+ZM4coNMQRtDokc3piUi3nw+OnPDRMcvznvjFreKJmDe95Vvtq6PO3EYGloqJCERERevHFF/XjH//YNz5z5kzl5eUpJyfH7zbl5eUqLy/3XfZ4PIqLi9OJEyeaHFhWbz/079fpvuU6Xa7Mvz7YpG0BANDSjPnR4ypv46oxNvOH3XXvsO5N3qbH41GHDh3qDSxtmryHc+jEiROqrKxUp06daox36tRJhYWFAW+TlpamhQsX+o1nZWUpIiKiSXWs3BkqKbj/qgAAoCVZuS1fl3yzv8m3Lysra9A8IwJLNYejZliwLMtvrNq8efM0e/Zs3+XqFZYRI0Y0eYXlk/CaKyzloWEa86PHa51/z6A43T0orkn7qs+6f3ymNf/4rMHzW0MtptRBLWbXQS212/3p17onI0/pLy7Q5YX+50K8H9tTU259RGvGX6Eru7YLSg2SWT0xpRZT6qivlvLQML+xmT/soRvPcoWlIYwILB06dFBoaKjfakpRUZHfqks1l8sll8vlN+50OuV0OptUx30jkuQICfnPOSwOh9/SV7VZwxM18yxft6vLzFtiZUVF/6eWOrSWWkypg1pqr8MR3a7WEwbPFOx3O5jSE9NqGdixk6KyP9WaQbdp3f8uUoj+c0ZAlRxaO2icoi+K0cCrLg3qOSwm9cSUWs7nx8/ZnsPS0OdsI94lFBYWpiuvvFLZ2dk1xrOzszVo0KBmrWXW8J4adXndb+sbdXnns/4fRC3ndx3UEtiM6xI1O7XufTTXWzNN6YlJtYSGOPTQ6GS93r2f7h77sHa7k1QaFq7d7iTdPfZhvd69nx4cnRz0E24lc3piUi08fupmxEm30rdva54wYYKeeuopDRw4UOvWrdMzzzyjDz/8UPHx9Z8NzeewtK5aTKmDWgJrbZ8jcb7VsuWD43ps8z4dLfH6xuIinXqQz2ExopbW9vg5r94lVG3NmjVaunSpjh8/rt69e+vJJ5/Utdde26DbnuvAIn3nLc46N0tf1NLy6qCWwL77R9eOP7bVTOmJSbVUVlnamV+krDf/qRGDB2hgj47NsrISiCk9MamW1vT4afDzt9VCFBcXW5Ks4uLic7rdla8dsLrN+au18rUD53S7TbF868dW/Jy/WMu3fmx3Kcb0hZ4EZkpf6ElgpvSloqLCyszMtCoqKmytw7LM6YllmXNfaS09aejzt1ErLGcjGCsspvF6vXr11Vd14403NvnE4paGngRGX/zRE3/0JDD64i+YPTnvPpofAACgNgQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPFsDyyHDx/WXXfdpYSEBIWHh6t79+565JFHVFFRYXdpAADAEG3sLuBf//qXqqqq9PTTT6tHjx764IMP9POf/1wnT57U7373O7vLAwAABrA9sNxwww264YYbfJcvueQS7d+/X2vXriWwAAAASQYElkCKi4sVExNT55zy8nKVl5f7Lns8HkmS1+uV1+sNan12qT6ulnp8TUFPAqMv/uiJP3oSGH3xF8yeNHSbDsuyrHO+97Nw6NAh9e3bV8uWLdPPfvazWuctWLBACxcu9BvPyMhQREREMEsEAADnSFlZmcaPH6/i4mJFRUXVOi9ogaW2QHGmXbt26aqrrvJdPnbsmIYMGaIhQ4boD3/4Q523DbTCEhcXpxMnTtR5wOczr9er7Oxspaamyul02l2OEehJYPTFHz3xR08Coy/+gtkTj8ejDh061BtYgvaS0L333qvbb7+9zjndunXz/X7s2DENGzZMAwcO1Lp16+rdvsvlksvl8ht3Op0t/g7WGo6xsehJYPTFHz3xR08Coy/+gtGThm4vaIGlQ4cO6tChQ4Pmfv755xo2bJiuvPJKrV+/XiEhtr/bGgAAGMT2k26PHTumoUOHqmvXrvrd736nL774wnddbGxsg7dT/cpW9cm3LZHX61VZWZk8Hg+p/9/oSWD0xR898UdPAqMv/oLZk+rn7frOULE9sGRlZSk/P1/5+fm6+OKLa1zXmNNrSkpKJElxcXHntD4AABB8JSUlio6OrvV6494l1FRVVVU6duyYIiMj5XA47C4nKKpPLP7ss89a7InFjUVPAqMv/uiJP3oSGH3xF8yeWJalkpISud3uOk8JsX2F5VwJCQnxW6FpqaKiongQfQc9CYy++KMn/uhJYPTFX7B6UtfKSjXObgUAAMYjsAAAAOMRWM4jLpdLjzzySMDPn2mt6Elg9MUfPfFHTwKjL/5M6EmLOekWAAC0XKywAAAA4xFYAACA8QgsAADAeAQWAABgPALLeSAtLU39+vVTZGSkOnbsqDFjxmj//v12l2WUtLQ0ORwOzZo1y+5SbPX555/rzjvvVPv27RUREaErrrhCu3fvtrssW50+fVoPPfSQEhISFB4erksuuUSPPvqoqqqq7C6t2eTm5mrUqFFyu91yOBzKzMyscb1lWVqwYIHcbrfCw8M1dOhQffjhh/YU20zq6onX69WcOXOUnJysCy64QG63WxMnTtSxY8fsK7iZ1HdfOdMvfvELORwOrVixollqI7CcB3JycjRt2jS9/fbbys7O1unTpzVixAidPHnS7tKMsGvXLq1bt06XX3653aXY6quvvtLVV18tp9Opv//97/roo4+0bNkytWvXzu7SbLVkyRI99dRTWr16tT7++GMtXbpUTzzxhH7/+9/bXVqzOXnypPr06aPVq1cHvH7p0qVavny5Vq9erV27dik2Nlapqam+72hrierqSVlZmfbs2aOHH35Ye/bs0csvv6wDBw7o5ptvtqHS5lXffaVaZmam/vnPf8rtdjdTZZIsnHeKioosSVZOTo7dpdiupKTESkxMtLKzs60hQ4ZYM2fOtLsk28yZM8e65ppr7C7DODfddJM1derUGmO33HKLdeedd9pUkb0kWZs2bfJdrqqqsmJjY63f/va3vrFTp05Z0dHR1lNPPWVDhc3vuz0J5J133rEkWUeOHGmeogxQW1+OHj1qdenSxfrggw+s+Ph468knn2yWelhhOQ8VFxdLkmJiYmyuxH7Tpk3TTTfdpOHDh9tdiu1eeeUVXXXVVbr11lvVsWNHpaSk6JlnnrG7LNtdc801ev3113XgwAFJ0vvvv6+33npLN954o82VmaGgoECFhYUaMWKEb8zlcmnIkCH6xz/+YWNlZikuLpbD4Wj1K5ZVVVWaMGGCHnjgAfXq1atZ991ivvywtbAsS7Nnz9Y111yj3r17212OrV544QXt2bNHu3btsrsUI3zyySdau3atZs+erd/85jd65513NGPGDLlcLk2cONHu8mwzZ84cFRcXKykpSaGhoaqsrNTjjz+uO+64w+7SjFBYWChJ6tSpU43xTp066ciRI3aUZJxTp05p7ty5Gj9+fKv/MsQlS5aoTZs2mjFjRrPvm8Bynrn33nu1d+9evfXWW3aXYqvPPvtMM2fOVFZWltq2bWt3OUaoqqrSVVddpcWLF0uSUlJS9OGHH2rt2rWtOrBs3LhRzz//vDIyMtSrVy/l5eVp1qxZcrvdmjRpkt3lGcPhcNS4bFmW31hr5PV6dfvtt6uqqkpr1qyxuxxb7d69WytXrtSePXtsuW/wktB5ZPr06XrllVe0fft2XXzxxXaXY6vdu3erqKhIV155pdq0aaM2bdooJydHq1atUps2bVRZWWl3ic2uc+fOuuyyy2qMff/739enn35qU0VmeOCBBzR37lzdfvvtSk5O1oQJE3TfffcpLS3N7tKMEBsbK+k/Ky3VioqK/FZdWhuv16tx48apoKBA2dnZrX515c0331RRUZG6du3q+7t75MgR/epXv1K3bt2Cvn9WWM4DlmVp+vTp2rRpk9544w0lJCTYXZLtrrvuOu3bt6/G2JQpU5SUlKQ5c+YoNDTUpsrsc/XVV/u93f3AgQOKj4+3qSIzlJWVKSSk5r/NQkNDW9XbmuuSkJCg2NhYZWdnKyUlRZJUUVGhnJwcLVmyxObq7FMdVg4ePKjt27erffv2dpdkuwkTJvidL3j99ddrwoQJmjJlStD3T2A5D0ybNk0ZGRnavHmzIiMjff8Sio6OVnh4uM3V2SMyMtLvHJ4LLrhA7du3b7Xn9tx3330aNGiQFi9erHHjxumdd97RunXrtG7dOrtLs9WoUaP0+OOPq2vXrurVq5fee+89LV++XFOnTrW7tGZTWlqq/Px83+WCggLl5eUpJiZGXbt21axZs7R48WIlJiYqMTFRixcvVkREhMaPH29j1cFVV0/cbrfGjh2rPXv26K9//asqKyt9f3djYmIUFhZmV9lBV9995bvBzel0KjY2Vpdeemnwi2uW9yLhrEgK+LN+/Xq7SzNKa39bs2VZ1l/+8herd+/elsvlspKSkqx169bZXZLtPB6PNXPmTKtr165W27ZtrUsuucR68MEHrfLycrtLazbbt28P+Ddk0qRJlmV9+9bmRx55xIqNjbVcLpd17bXXWvv27bO36CCrqycFBQW1/t3dvn273aUHVX33le9qzrc1OyzLsoIfiwAAAJqOk24BAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMN7/DxGGJxTxVoizAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(3, 1, sharex=True)\n",
    "fig.subplots_adjust(hspace=0.1)\n",
    "\n",
    "y_EPS = 1E5*EPS\n",
    "\n",
    "y_lower = max((fabs(min(scipy_rel_error)), fabs(min(wol_rel_error))))\n",
    "\n",
    "axs[0].set_ylim( scipy_rel_error[12]-y_EPS, scipy_rel_error[12]+y_EPS)\n",
    "axs[1].set_ylim( scipy_rel_error[6]-8*EPS, scipy_rel_error[6]+4*EPS)\n",
    "axs[2].set_ylim( -y_EPS, y_lower+y_EPS)\n",
    "\n",
    "axs[1].legend([\"Scipy ext\", \"Wolframe\"])\n",
    "\n",
    "axs[0].spines.bottom.set_visible(False)\n",
    "axs[1].spines.top.set_visible(False)\n",
    "axs[1].spines.bottom.set_visible(False)\n",
    "axs[2].spines.top.set_visible(False)\n",
    "\n",
    "axs[0].xaxis.tick_top()\n",
    "axs[2].xaxis.tick_bottom()\n",
    "\n",
    "markerline1, stemlines1, baseline1 = axs[0].stem(np.arange(1, len(scipy_rel_error)+1), scipy_rel_error)\n",
    "markerline2, stemlines2, baseline2 = axs[0].stem(np.arange(1, len(scipy_rel_error)+1), wol_rel_error,  markerfmt='D')\n",
    "markerline3, stemlines3, baseline3 = axs[1].stem(np.arange(1, len(scipy_rel_error)+1), scipy_rel_error)\n",
    "markerline4, stemlines4, baseline4 = axs[1].stem(np.arange(1, len(scipy_rel_error)+1), wol_rel_error,  markerfmt='D')\n",
    "markerline5, stemlines5, baseline5 = axs[2].stem(np.arange(1, len(scipy_rel_error)+1), scipy_rel_error)\n",
    "markerline6, stemlines6, baseline6 = axs[2].stem(np.arange(1, len(scipy_rel_error)+1), wol_rel_error,  markerfmt='D')\n",
    "\n",
    "markerline1.set_markerfacecolor('red')\n",
    "markerline3.set_markerfacecolor('red')\n",
    "markerline5.set_markerfacecolor('red')\n",
    "\n",
    "axs[0].grid()\n",
    "axs[1].grid()\n",
    "axs[2].grid()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "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.8.16"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "38df5b977d726286e1d8f7be28e0a7ab528a11cffc8b179ff55f1f026c4c21da"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}