Tutorial-Min_Max_and_Piecewise_Expressions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Removing `if` Statements from Expressions \n",
    "\n",
    "## Author: Patrick Nelson\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [Min_Max_and_Piecewise_Expressions.py](../edit/Min_Max_and_Piecewise_Expressions.py) Contains functions that can be used to compute the minimum or maximum of two values and to implement piecewise-defined expressions \n",
    "\n",
    "## Introduction:\n",
    "\n",
    "Conditional statements are a critical tool in programming, allowing us to control the flow through a program to avoid pitfalls, code piecewise-defined functions, and so forth. However, there are times when it is useful to work around them. It takes a processor time to evaluate the whether or not to execute the code block, so for some expressions, performance can be improved by rewriting the expression to use an absolute value function in a manner upon which we will expand in this tutorial. Even more relevant to NRPy+ are piecewise-defined functions. These inherently involve `if` statements, but NRPy+'s automatic code generation cannot handle these by itself, requiring hand-coding to be done. However, if it is possible to rewrite the expression in terms of absolute values, then NRPy+ can handle the entire thing itself. \n",
    "\n",
    "The absolute value is a function that simply returns the magnitude of its argument, a positive value. That is, \n",
    "\\begin{align}\n",
    "|x|&= \\left \\{ \\begin{array}{lll}x & \\mbox{if} & x \\geq 0 \\\\\n",
    "-x & \\mbox{if} & x \\leq 0 \\end{array} \\right. \\\\\n",
    "\\end{align}\n",
    "\n",
    "In C, this is implemented as `fabs()`, which merely has to make the first bit of a double-precision floating point number 0, and is thus quite fast. \n",
    "\n",
    "There are myriad uses for these tricks in practice. One example comes from GRMHD (and, by extension, the special cases of GRFFE and GRHD), in which it is necessary to limit the velocity of the plasma in order to keep the simulations stable. This is done by calculating the Lorentz factor $\\Gamma$ of the plasma and comparing to some predefined maximum $\\Gamma_\\max$. Then, if\n",
    "$$\n",
    "R = 1-\\frac{1}{\\Gamma^2} > 1-\\frac{1}{\\Gamma_{\\max}^2} = R_\\max,\n",
    "$$\n",
    "we rescale the velocities by $\\sqrt{R_\\max/R}$. In NRPy+, we instead always rescale by\n",
    "$$\n",
    "\\sqrt{\\frac{\\min(R,R_\\max)}{R+\\epsilon}},\n",
    "$$\n",
    "which has the same effect while allowing the entire process to be handled by NRPy+'s automatic code generation. ($\\epsilon$ is some small number chosen to avoid division by zero without affecting the results otherwise.) See [here](Tutorial-GRHD_Equations-Cartesian.ipynb#convertvtou) for more information on this. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#min_max): Minimum and Maximum\n",
    "    1. [Step 1.a](#confirm): Confirm that these work for real numbers\n",
    "1. [Step 2](#piecewise): Piecewise-defined functions\n",
    "1. [Step 3](#sympy): Rewrite functions to work with symbolic expressions\n",
    "1. [Step 4](#validation): Validation against `Min_Max_and_Piecewise_Expressions` NRPy+ module\n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='min_max'></a>\n",
    "\n",
    "# Step 1: Minimum and Maximum \\[Back to [top](#toc)\\]\n",
    "$$\\label{min_max}$$\n",
    "\n",
    "Our first job will be to rewrite minimum and maximum functions without if statements. For example, the typical implementation of `min(a,b)` will be something like this:\n",
    "```python\n",
    "def min(a,b):\n",
    "    if a<b:\n",
    "        return a\n",
    "    else:\n",
    "        return b\n",
    "```\n",
    "However, to take full advantage of NRPy+'s automated function generation capabilities, we want to write this without the `if` statements, replacing them with calls to `fabs()`. We will define these functions in the following way: \n",
    "$$\\boxed{\n",
    "\\min(a,b) = \\tfrac{1}{2} \\left( a+b - \\lvert a-b \\rvert \\right)\\\\\n",
    "\\max(a,b) = \\tfrac{1}{2} \\left( a+b + \\lvert a-b \\rvert \\right).}\n",
    "$$\n",
    "\n",
    "<a id='confirm'></a>\n",
    "\n",
    "## Step 1.a: Confirm that these work for real numbers \\[Back to [top](#toc)\\]\n",
    "$$\\label{confirm}$$\n",
    "\n",
    "For real numbers, these operate exactly as expected. In the case $a>b$,\n",
    "\\begin{align}\n",
    "\\min(a,b) &= \\tfrac{1}{2} \\left( a+b - (a-b) \\right) = b \\\\\n",
    "\\max(a,b) &= \\tfrac{1}{2} \\left( a+b + (a-b) \\right) = a, \\\\\n",
    "\\end{align}\n",
    "and in the case $a<b$, the reverse holds:\n",
    "\\begin{align}\n",
    "\\min(a,b) &= \\tfrac{1}{2} \\left( a+b - (b-a) \\right) = a \\\\\n",
    "\\max(a,b) &= \\tfrac{1}{2} \\left( a+b + (b-a) \\right) = b, \\\\\n",
    "\\end{align}\n",
    "\n",
    "In code, we will represent this as:\n",
    "```\n",
    "min_noif(a,b) = sp.Rational(1,2)*(a+b-nrpyAbs(a-b))\n",
    "max_noif(a,b) = sp.Rational(1,2)*(a+b+nrpyAbs(a-b))\n",
    "```\n",
    "For demonstration purposes, we will use `np.absolute()` and floating point numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The minimum of 5.0 and 10.0 is 5.0\n"
     ]
    }
   ],
   "source": [
    "import numpy as np               # NumPy: Python module specializing in numerical computations\n",
    "import matplotlib.pyplot as plt  # matplotlib: Python module specializing in plotting capabilities\n",
    "\n",
    "thismodule = \"Min_Max_and_Piecewise_Expressions\"\n",
    "# First, we'll write the functions. Note that we are not using sympy right now. For NRPy+ code generation,\n",
    "# use the expressions above.\n",
    "def min_noif(a,b):\n",
    "    return 0.5 * (a+b-np.absolute(a-b))\n",
    "\n",
    "def max_noif(a,b):\n",
    "    return 0.5 * (a+b+np.absolute(a-b))\n",
    "\n",
    "# Now, let's put these through their paces.\n",
    "a_number = 5.0\n",
    "another_number = 10.0\n",
    "print(\"The minimum of \"+str(a_number)+\" and \"+str(another_number)+\" is \"+str(min_noif(a_number,another_number)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feel free to test other cases above if you'd like. Note that we use a suffix, `_noif`, to avoid conflicts with other functions. When using this in NRPy+, make sure you use `sp.Rational()` and the `nrpyAbs()` function, which will always be interpreted as the C function `fabs()` (Sympy's `sp.Abs()` may get interpreted as $\\sqrt{zz^*}$, for instance).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='piecewise'></a>\n",
    "\n",
    "# Step 2: Piecewise-defined functions \\[Back to [top](#toc)\\]\n",
    "$$\\label{piecewise}$$\n",
    "\n",
    "Next, we'll define functions to represent branches of a piecewise-defined function. For example, consider the function \n",
    "\\begin{align}\n",
    "f(x) &= \\left \\{ \\begin{array}{lll} \\frac{1}{10}x^2+1 & \\mbox{if} & x \\leq 0 \\\\\n",
    "\\exp(\\frac{x}{5}) & \\mbox{if} & x > 0 \\end{array} \\right. , \\\\\n",
    "\\end{align}\n",
    "which is continuous, but not differentiable at $x=0$. \n",
    "\n",
    "To solve this problem, let's add the two parts together, multiplying each part by a function that is either one or zero depending on $x$. To define $x \\leq 0$, this can be done by multiplying by the minimum of $x$ and $0$. We also will need to normalize this. To avoid putting a zero in the denominator, however, we will add some small $\\epsilon$ to the denominator, i.e.,\n",
    "$$\n",
    "\\frac{\\min(x,0)}{x-\\epsilon}\n",
    "$$\n",
    "This $\\epsilon$ corresponds `TINYDOUBLE` in NRPy+; so, we will define the variable here with its default value, `1e-100`. Additionally, to get the correct behavior on the boundary, we shift the boundary by $\\epsilon$, giving us\n",
    "$$\n",
    "\\frac{\\min(x-\\epsilon,0)}{x-\\epsilon}\n",
    "$$\n",
    "\n",
    "The corresponding expression for $x > 0$ can be written as \n",
    "$$\n",
    "\\frac{\\max(x,0)}{x+\\epsilon},\n",
    "$$\n",
    "using a positive small number to once again avoid division by zero. \n",
    "When using these for numerical relativity codes, it is important to consider the relationship between $\\epsilon$, or `TINYDOUBLE`, and the gridpoints in the simulation. As long as $\\epsilon$ is positive and large enough to avoid catastrophic cancellation, these functional forms avoid division by zero, as proven [below](#proof).\n",
    "\n",
    "So, we'll code NumPy versions of these expressions below. Naturally, there are many circumstances in which one will want the boundary between two pieces of a function to be something other than 0; if we let that boundary be $x^*$, this can easily be done by passing $x-x^*$ to the maximum/minimum functions. For the sake of code readability, we will write the functions to pass $x$ and $x^*$ as separate arguments. Additionally, we code separate functions for $\\leq$ and $<$, and likewise for $\\geq$ and $>$. The \"or equal to\" versions add a small offset to the boundary to give the proper behavior on the desired boundary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deVxU1f/H8deZQRaXEv2akriUFq6ISyqmBiLhli1qZpqUu0luWZqa+5ZpaZa5pJmpaWKLX5dSUdRkUjGX1DL9llvhkpFLiixzfn/MwE8NFGRmLjCf5+MxD2Dmztw3d4bPHM49c47SWiOEEMJ9mIwOIIQQwrWk8AshhJuRwi+EEG5GCr8QQrgZKfxCCOFmPIwOkB3/+c9/dMWKFY2OIYQQ+cqePXv+1FqXuvX6fFH4K1asSHx8vNExhBAiX1FKncjseunqEUIINyOFXwgh3IwUfiGEcDP5oo9f5F0pKSmcPn2apKQko6MIA3h7e+Pv70+hQoWMjiJyQAq/yJXTp09TrFgxKlasiFLK6DjChbTWXLhwgdOnT/PAAw8YHUfkgHT1iFxJSkqiZMmSUvTdkFKKkiVLyn97+VCBLvwWC0yebPsqnEeKvvuS5955nFm/CmxXj8UCISGppKSY8PY2ERMDwcFGpxJCiDuzWCAsDJKTwdMTh9evAtvij42F5GSF1iaSkzWxsUYnEs5iNpsJCgqiRo0adOjQgatXrxIfH0///v1dmmPOnDksXrw414/z2muvUb16dV577TUHpLLZt28f69aty/h59erVTJkyxWGPLxzLVr8gLc321dH1q8C2+ENCwMsLrl9PwWSCkBAZdVBQ+fj4sG/fPgA6d+7MnDlzGDx4MPXq1XNpjj59+jjkcebNm8dff/2F2Wx2yOOBrfDHx8fTqlUrANq2bUvbtm0d9vjCsUJCbC399BZ/SIhjH7/AtviDg2HzZhP+/vMpXbozDRpYjY4kXKBJkyYcO3aM2NhY2rRpA8A///xDt27dqF+/PrVr1+brr78GIC0tjSFDhlCjRg0CAwOZNWsWAHv27OGxxx6jbt26REREkJCQwLlz56hbty4A+/fvRynFyZMnAahUqRJXr15lzJgxTJs2DYD33nuPatWqERgYyHPPPXfbHDdq27YtV65coW7duqxYsYIXX3yR6OjojNuLFi0KQGxsLCEhIbRv354qVarQuXNn0lfT2717N40aNaJWrVrUr1+fixcvMmrUKFasWEFQUBArVqxg0aJFREVFAXD8+HGaNWtGYGAgYWFhGb/Xiy++SP/+/WnUqBEPPvjgTTmEcwUH27p3xo93fDcPFOAWP0CjRoqpU315/vmVrF8fSevWrY2OVKANHDgwo+XtKEFBQcyYMSNb26amprJ+/XpatGhx0/UTJ06kWbNmLFy4kL///pv69evTvHlzFi9ezPHjx9m3bx8eHh789ddfpKSk8Morr/D1119TqlQpVqxYwYgRI1i4cCFJSUlcunSJ7du3U69ePbZv307jxo257777KFy48E37nDJlCr/99hteXl78/ffft81RpEiRjPutXr2aokWLZhzH9evXZ/n77t27l0OHDnH//ffz6KOPsmPHDurXr0/Hjh1ZsWIFjzzyCJcuXaJw4cKMGzeO+Ph43n//fQAWLVqU8TivvPIKkZGRREZGsnDhQvr3789XX30FQEJCAt999x0///wzbdu2pX379tl6LkTuBQc777xkgW3xp2vfvj3+/v68++67RkcRTnLt2jWCgoKoV68e5cuXp3v37jfdvmHDBqZMmUJQUBAhISEkJSVx8uRJNm3aRO/evfHwsLV/SpQowZEjRzh48CDh4eEEBQUxYcIETp8+DUCjRo3YsWMH27ZtY/jw4Wzbto3t27fTpEmTf2UKDAykc+fOLFmyJOPxs8pxt+rXr4+/vz8mk4mgoCCOHz/OkSNH8PPz45FHHgHgnnvuydh/ViwWC88//zwAL7zwAt99913GbU899RQmk4lq1apx9uzZu84q8pYC3eIHKFSoEFFRUQwbNowDBw4QGBhodKQCK7stc0e7sY8/M1prVq1aRUBAwB0fS2tN9erVsWQyhq5p06Zs376dEydO8OSTT/LWW2+hlMr0P8m1a9eybds2/vvf/zJx4kR+/PHHHOVI5+HhgdVq66a0Wq0kJydn3Obl5ZXxvdlsJjU1NduPm1037iO9K0nkfwW+xQ/Qs2dPChcubFhhEsaKiIhg1qxZGYVr7969AISHhzN37tyMgvnXX38REBDA+fPnMwp/SkoKhw4dAmznD5YsWcJDDz2EyWSiRIkSrFu3jsaNG9+0P6vVyqlTpwgNDeWtt97i4sWLXLlyJcsct1OxYkX27NkD2LqBUlJSbrt9QEAACQkJ7N69G4DLly+TmppKsWLFuHz5cqb3adSoEcuXLwdg6dKlmf4HIwoWtyj8JUqUIDIykqVLl3Lu3Dmj4wgXe/PNN0lJSSEwMJDq1avz5ptvAtCjRw/Kly9PYGAgtWrVYtmyZXh6ehIdHc3QoUOpVasWQUFBxMXFAbYirLWmadOmADRu3JjixYvj6+t70/7S0tLo0qULNWvWpHbt2vTv35/ixYtnmeN2evbsydatW6lVqxYWi+Wm8wGZ8fT0ZMWKFbzyyivUqlWL8PBwkpKSCA0N5fDhwxknd280a9YsPv74YwIDA/n000+ZOXNmto+tyJ9Ufvj3rV69ejq3C7EcOXKEKlWqMGbMGEaPHu2gZOKnn36iatWqRscQBpLXQN6llNqjtf7XuGa3aPGD7V/g1q1bM3v2bJlbRAjh1tym8INtuOG5c+cy+jOFEMIduVXhDwsLo2bNmrz77rsyQkEI4bbcqvArpRg4cCAHDhxgy5YtRscRQghDuFXhB3j++ecpVaqUDO0UQrgttyv83t7e9O3blzVr1nD06FGj4wghhMs5rfArpRYqpc4ppQ7ecF0JpdRGpdRR+1ff2z2Gs7z88ssUKlRIxiu7sd9++40GDRpQuXJlOnbseNMnYtNduHCB0NBQihYtmjGhGfz/J1jHjBlz08+3k5qaSuvWrfnPf/7DwYMHb7rttddeo0qVKgQGBvL0009nzO2THSNGjKBcuXIZk7cJkR3ObPEvAlrcct0wIEZr/RAQY//Z5UqXLs3zzz/Pxx9/TGJiohERhMGGDh3KoEGDOHbsGL6+vixYsOBf23h7ezN+/PiMGTfTbdiwgREjRnD16lU++uijbHUb9u3blypVqvDVV1/RsWPHjPl/wPYJ4oMHD3LgwAEefvhhJk+e/K/7jxkz5qaJ1dI98cQT7Nq1Kxu/sRD/z2mFX2u9DfjrlqufBD6xf/8J8JSz9n8nAwcOzPjDFa7lyCXlRo0adVPhHTFixB3/k9Nas3nz5oyZJiMjIzNmo7xRkSJFaNy4Md7e3jddHxERQUREBDNnzuTChQsMGjSIEydO8NBDD/Hnn39itVpp0qQJGzZsAGDs2LHce++9TJ8+ncaNG/PRRx/RqVMnLl68CMDjjz+eMZFaw4YNb3pTuJOGDRvi5+eX7e2FANdP0lZaa51g//4MUDqrDZVSvYBeAOXLl3d4kFq1ahEaGsqsWbMYNGjQHWcwFI7h6CXlunXrxjPPPMPAgQOxWq0sX76czZs3ExQUlOn2y5Yt47777qN48eIZz7m/vz+///57tve5ceNGYmNj6d+/PyVLlmTmzJkMGDCAoUOH0rdvX+rXr0+1atV4/PHHAf71SfHg4GC2b9+e6WMvXLiQjh07ZjuLEHfDsGqntdZKqSw7R7XW84B5YJuywRkZBg0aRNu2bVm1apX8sblIZkvK5abwV6xYkZIlS7J3717Onj1L7dq1qVChwm1n6/zzzz/vfodA8+bNCQ8PZ8yYMfTo0SOjj79Hjx6sXLmSOXPm3NW6BBMnTsTDw4POnTsD8OOPP/LCCy8AcObMGTw9PTP+u4mJiaFkyZK5+j2E+3J14T+rlPLTWicopfwAQ2dMa926NZUrV2bGjBlS+F3EGUvK9ejRg0WLFnHmzBm6devG5cuXs5xhctmyZVStWpW///6b1NRUPDw8OH36NGXLls32/pRSwP+f3E3/+erVqxndNFeuXKFYsWLZfsxFixaxZs0aYmJiMh6vZs2aGW8gY8aMoWLFirz44ovZfkwhsuLqwr8aiASm2L/+e+05FzKZTAwYMIBXXnmF77//noYNGxoZxy2kLykXG2sr+o5YYejpp59m1KhRpKSksGzZMsxm8x1b3KGhoURHR/Pcc8/xySef8OSTT+Y6x9ChQ+ncuTMVKlSgZ8+erFmzJlv3++abb5g6dSpbt27910peQjiF1topF+AzIAFIAU4D3YGS2EbzHAU2ASWy81h169bVznL58mVdvHhx/eyzzzptHwXZ4cOHjY6gtda6d+/eeujQodne/n//+59+5JFHdKVKlXT79u11UlKS1lrrr7/+Wr/55psZ21WoUEH7+vrqIkWK6LJly+pDhw5l+nixsbG6QYMGOjU1VWut9dNPP60XLlyYrSyVKlXS/v7+ulatWrpWrVq6d+/e/9pm9OjR+uOPP/7X9a+99pouW7asVkrpsmXL6tGjR2drn46UV14D4t+AeJ1JTXWbaZlv5/XXX+edd97h119/dcqJ5IIsL0zJa7VaqVOnDitXruShhx4yNIs7yguvAZE5t5+W+XbSP5yTvhC1yD8OHz5M5cqVCQsLk6IvRDZJ4cc2XLRdu3bMmzePK1euGB1H5EC1atX49ddfmT59utFRhMg3pPDbDRo0iIsXL2b66UghhChIpPDbNWzYkAYNGjBz5kysVqvRcYQQwmmk8N8gfe6WtWvXGh1FCCGcRgr/Ddq1a0e5cuV49913jY4ihBBOI4X/Bh4eHkRFRbFlyxb2799vdByRTwwcOJBt27bdcbspU6awdOlSFi1aRKlSpQgKCiIoKChjosDz58/TosWtE9oK4XhS+G/Rs2dPChcuLCt0iWy5cOEC33//PU2bNr3jtt9++23GxG0dO3Zk37597Nu3jx49egBQqlQp/Pz82LFjh1MzCyGF/xa+vr689NJLLFu2jDNnzhgdp0CynLIweftkLKdyPy/z3UzLDLBnzx4ee+wx6tatS0REBAkJCVy8eJGAgACOHDkCQKdOnZg/fz4ARYsWZdCgQVSvXp2wsDDOnz8PwKpVqzJa6be7/6VLl0hOTqZUqVK3zfXUU0+xdOnSHB4FIXIos4/z5rWLM6dsyMwvv/yioaFu1myjjotz6a7znZx+XD/uZJz2meCjzWPN2meCj447mbsD/Ntvv+natWtrrbVOS0vTDz74oD5+/HjG9Ae3Xg4dOqSTk5N1cHCwPnfunNZa6+XLl+uXXnpJa631hg0bdMOGDfVnn32mIyIiMvYD6CVLlmittR47dqzu16+f1lrrrl276tWrV2dsl9X9V61alTEVxMcff6zLlCmja9asqdu1a6dPnjyZsd3p06d1jRo1cnVMXE2mbMhaXJzWkyZpw+oIWUzZYHhRz87F1YU/Lk5rkylJQ4r28bFK8b+NnP7RT9o2SZvHmjVj0OaxZj1p26RcZ2jevLn+4Ycf9Pr163W7du3uuP2PP/6oixUrlvFmUKNGDR0eHp5xe8+ePXWJEiX0qVOnMq4zmUw6JSVFa22b56dWrVpaa63Dw8O1xWK56fEzu3/Pnj11nP2F9Oeff2bMDTRnzhwdGhqasV1ycrIuUaJETg+BoaTwZy4uTmsfH63NZttXI+pIVoVfVh/JRGwsgCeguH7dSmyscsgskgJCKobgafYkOS0ZT7MnIRVDcv2YOZ2WWWtN9erVsWSyBJjVauWnn36icOHCJCYm4u/vn+njpE+d7OPjQ1JS0h3vv2vXLj788EOAm+bR79GjB6+//nrGz0lJSfj4+OTwCIi8yNFrTziS9PFnIiQEvLwUkIrVep3g4OtGRyowgssFE9M1hvGh44npGkNwudz/JTz99NN888037N69m4iICIoVK5Zx4vTWS7Vq1QgICOD8+fMZhT8lJYVDhw4B8O6771K1alWWLVvGSy+9REpKCmAr6NHR0YDtzaNx48YAVK1alWPHjmVkyez+hw4dokqVKpjNZgASEhIytl+9evVNE5z98ssv1KhRI9fHRBgvfe0Js9lxa084irT4M5E+Z/zChSf46KMuHDrUhZCQfkbHKjCCywU7pOCn8/T0JDQ0lOLFi2cU1zttHx0dTf/+/bl48SKpqakMHDgQDw8PPvroI3bt2kWxYsVo2rQpEyZMYOzYsRQpUoRdu3YxYcIE7rvvPlasWAHYFvOZO3cuPXr04MiRI5nev1ixYjcN03zvvfdYvXo1Hh4elChR4qZpQrZs2ULr1q0ddmyEcZyx9oTDZNb/k9curu7jT2e1WnXjxo21v79/Rp+suFle6N9NS0vTtWrV0r/88ovT9lGkSJEsb3v00Ud1YmJilrc3b95c//HHH9naT5MmTfRff/2V43xGyguvAZE5sujjl66e21BKMWrUKE6fPi2Tt+VReWFa5unTp3Py5Mksb9+4cSN+fn53fJzz588zePBgfH19HRlPiH+RhVjuQGtNo0aN+OOPPzh69Cienp6G5MirZBEOIa+BvEsWYrlL6a3+kydP8umnnxodJ0/KD40H4Rzy3OdPUvizoUWLFtSrV4+JEydmjPIQNt7e3ly4cEEKgBvSWnPhwgW8vb2NjiJySEb1ZEN6q79t27YsW7aMyMhIoyPlGf7+/pw+fTpjCgPhXry9vbP8rIPIu6SPP5u01tSpU4d//vmHw4cP4+Eh75lCiLxN+vhzKb3Vf/To0Ywx3EIIkR9Jiz8HrFYrQUFBpKSkcPDgwWx9WEgIIYwiLX4HMJlMvPnmm/z8888ZH98XQoj8Rlr8OWS1WqlZsyYmk4n9+/djMsl7pxAib5IWv4OYTCZGjhzJwYMH+eqrr4yOI4QQOSYt/ruQlpZGtWrV8PHxYe/evRlT9AohRF4iLX4HMpvNjBw5kv379/Pf//7X6DhCCJEjhhR+pdQgpdQhpdRBpdRnSql899G/Tp06UalSJcaNGyefWhVC5CsuL/xKqbJAf6Ce1roGYAaec3WO3PLw8GDEiBHs2bOH9evXGx1HCCGyzaiuHg/ARynlARQG/jAoR6506dKFihUrSqtfCJGvuLzwa61/B6YBJ4EE4KLWesOt2ymleiml4pVS8Xl1HphChQoxfPhwdu7cycaNG42OI4QQ2eLyUT1KKV9gFdAR+BtYCURrrZdkdZ+8NqrnRsnJyVSuXJly5crx3XffyQgfIUSekZdG9TQHftNan9dapwBfAI0MyOEQnp6eDBs2jLi4OLZs2WJ0HCGEuCMjCv9JoKFSqrCyNY/DgJ8MyOEw3bp14/7772fcuHFGRxFCiDsyoo9/JxAN/AD8aM8wz9U5HMnb25uhQ4eydetWtm7danQcIYS4LfnkroNcu3aNBx54gBo1arBp0yaj4wghRJ7q4y+QfHx8eP3114mJiWHHjh1GxxFCiCxJ4Xeg3r17U6pUKcaPH290FCGEyJIUfgcqUqQIQ4YM4dtvv2Xnzp1GxxFCiExJ4Xewl19+mXvuiSAy8mcsFqPTCCFyy2KByZNx+d+z5ZSFidsmYjnl+B3LiuEO9uOPRbl27b8cOaIIDU1jyxYzwcFGpxJC3A2LBcLCIDkZPD0hJgaX/D1bTlkIXRTK9dTreBfyZnPkZoLLOW7H0uJ3sNhYsFo9AA+uX9ds2ZL3R00JITIXG2sr+mlptq+xsa7Zb8yvMVxPvQ4mSLGmEHvcsTuWwu9gISHg6akwmaxAMibTNqMjCSHuku3vGcxm29eQENfs98z3ZyANTJjwNHsSUtGxO5Zx/E5gscDmzWl89NELeHjs5tChQ3h6ehodSwhxFywWW0s/JMQ13TzHjx+nevXq1HmyDi37tiS0Yuhdd/NkNY5fCr8TrV+/nlatWjFjxgwGDBhgdBwhRB6ntaZNmzZs3bqVw4cPU758+Vw9nnyAywAtWrQgPDyccePGkZiYaHQcIUQet2rVKtatW8e4ceNyXfRvRwq/EymlePvtt0lMTGTixIlGxxFC5GEXL16kf//+BAUF0b9/f6fuSwq/k9WqVYuXXnqJWbNm8euvvxodRwiRR40YMYIzZ84wb948PDycO9JeCr8LjB8/Hg8PD4YNG2Z0FCFEHrRr1y5mz55NVFQUjzzyiNP3J4XfBe6//35ef/11Vq5cSVxcnNFxhBB5SGpqKr169cLPz48JEya4ZJ9S+F1kyJAh+Pn58eqrr8rC7EKIDDNnzmT//v3MmjWLe+65xyX7lMLvIkWKFGHChAl8//33rFy50ug4Qog84MSJE4waNYonnniCp59+2mX7lXH8LpSWlkadOnW4fPkyP/30E15eXkZHEkIYRGvNE088QWxsrEPG7GdGxvHnAWazmenTp/Pbb7/x/vvvGx1HCGGgL774grVr1zp9zH5mpMVvgFatWmGxWDh27BglS5Y0Oo4QwsUuXrxI1apVKV26NLt373ba8E1p8echb7/9NpcuXWLcuHFGRxFCGGDkyJGcOXOGuXPnOn3Mfmak8BugevXq9OzZk9mzZ/PLL78YHUcI4UK7d+/mgw8+oF+/ftSvX9+QDNLVY5CzZ89SuXJlwsPD+eKLL4yOI4RwgdTUVB555BHOnTvH4cOHuffee526P+nqyWNKly7NsGHD+PLLL9m2TebsF8IdvPfee+zbt4/33nvP6UX/dqTFb6CrV68SEBBAmTJl2LlzJyaTvA8LUVCdOHGCatWq0axZM1avXo1Syun7lBZ/HlS4cGEmTZpEfHw8n332mdFxhBBOorUmKioKgPfff98lRf92pPAbrHPnztSpU4fhw4dz7do1o+MIIZzgyy+/ZM2aNYwbN44KFSoYHUcKv9FMJhPTp0/n5MmTzJw50+g4QggHu3TpEq+88gq1atXKMyvxGVL4lVLFlVLRSqmflVI/KaVcsJJl3hUSEkLbtm2ZNGkS586dMzqOEMJBLKcstJjYgj/Mf7hknv3sMqrFPxP4RmtdBagF/GRQjjxj6tSpXLt2jbFjxxodRQjhAJZTFkIXhWLxtuDRzYM0vzSjI2VweeFXSt0LNAUWAGitk7XWf7s6R14TEBBAnz59mDt3Lj/95Pbvg0Lke98c+YbrqdfBBNqkiT0ea3SkDEa0+B8AzgMfK6X2KqU+UkoVuXUjpVQvpVS8Uir+/Pnzrk9pgFGjRlGkSBF69lzA5MlgsRidSIiCwWLBpX9TWmu2L94OaWDChKfZk5CKIa7ZeTa4fBy/Uqoe8D3wqNZ6p1JqJnBJa/1mVvcpqOP4M9Ov3xJmz34Gk8kbLy8TMTEQ7NZnQITIHYsFwsIgORk8PXHJ39R7773HgAEDeHnSy/g39iekYgjB5Vz/h5zVOH4jzjScBk5rrXfaf44GZDFauzJlOgIKq9VEcrImNlZJ4RciF2JjbUU/Lc32NTbWuYV/586dDBkyhCeeeIJZQ2flyQ9mujyR1voMcEopFWC/Kgw47OoceVXz5oXw8lJACkqlEBJidCIh8reQEFtL32y2fXXm39SFCxfo0KEDZcuW5ZNPPsmTRR+MG9XzCrBUKXUACAImGZQjzwkOhi1bzNSrtxqrtRleXj8YHUmIfC042Na9M368c7t5rFYrL7zwAmfPnmXlypX4+vo6Z0cOIHP15FGJiYlUq1YNPz8/du3alWfG/wohMjdx4kRGjhzJ7Nmz6du3r9FxAJmrJ9/x9fXl/fffZ+/evbzzzjtGxxFC3MaWLVsYNWoUnTp1ok+fPkbHuaM7tviVUq8AS7TWia6J9G/u2OJP98wzz7B+/XoOHDjAQw89ZHQcIcQtEhISqF27Nr6+vuzevZuiRYsaHSlDblr8pYHdSqnPlVItlNHTyrmZ999/Hy8vL3r16kV+6JYTwp2kpqby3HPPcfnyZaKjo/NU0b+dOxZ+rfVI4CFsn7R9ETiqlJqklKrk5GwCuP/++3n77beJjY1lwYIFRscRQtxg5MiRbNu2jblz51K9enWj42Rbtvr4ta2pecZ+SQV8gWil1FQnZhN2PXr0ICQkhCFDhvDHH38YHUcIAaxZs4a33nqLXr160aVLF6Pj5MgdC79SaoBSag8wFdgB1NRa9wXqAu2cnE8ASinmz5/P9evXMxZzEEIY5/jx43Tt2pXatWvny+nUs9PiLwE8o7WO0Fqv1FqnAGitrUAbp6YTGSpXrszYsWP58ssvWbVqldFxhHBb169fp0OHDlitVlauXIm3t7fRkXIsO338o7XWJ7K4TaaRdKHBgwdTu3ZtoqKiSEw0bJCVEG5t8ODBxMfHs2jRIipVyp+nOmUcfz7i4eHBggULOH/+PK+99prRcYRwO8uXL2f27Nm8+uqrPPXUU0bHuWtS+POZ2rVrM2TIEBYsWMDmzZuNjiOE2/j555/p0aMHjz76KJMnTzY6Tq7IlA350LVr1wgMDMRqtfLjjz9SuHBhoyMJUaD9888/NGjQgHPnzrF3717Kli1rdKRskSkbChAfHx/mz5/Pr7/+ypgxY4yOI0SBFncyjsZvNObQpUMsW7Ys3xT925HCn0+FhITQs2dPpk+fzp49e4yOI0SBZDllIeTjEPb57sOjuwdFAv61WGC+JIU/H5s6dSqlS5eme/fupKSkGB1HiALnw28+JMWaYls3V+WtdXNzQwp/Pla8eHE++OAD9u/fz7Rp04yOI0SBsnXrVlZMWYHSCrMy57l1c3NDTu4WAO3bt2fNmjUcOHCAhx9+2Og4QuR7e/bsITQ0lHLlyjFtxTT2Je4zbN3c3Mjq5K4U/gIgISGBatWqERgYyJYtW/Lscm9C5Ac///wzTZo0oUiRIuzYsSNfn8yVUT0FmJ+fH9OmTWPbtm3Mnz/f6DhC5FsnTpwgPDwcs9nMpk2b8nXRvx0p/AVEt27daNasGYMHr+SNNy5isRidSAjHslhg8mSc9to+e/Ys4eHhXLlyhQ0bNlC5cmXn7CgPkK6eAmTlytM8+2wJwAsfHxMxMcppC0sL4UoWC4SFQXIyeHo6ftH0v//+m9DQUI4cOcKmTZto1KiR4x7cQNLV4waOHfNHKW/ATFKSldhYoxMJ4RixsRzaUHwAABk/SURBVLain5Zm++rI1/bVq1d54oknOHToEF9++WWBKfq3I4W/AAkJAW9vhVJpaH2dYsXkg12iYAgJsbX0zWbb15AQxzxucnIy7du3Z8eOHSxdupSIiAjHPHAeJ4W/AAkOhpgYxahRqZQv340JE1qTkJBgdCwhcs322obx4x3XzZOWlkbXrl1Zv349c+fOpUOHDrl/0HxC+vgLqIMHD1K/fn3q16/Ppk2b8PDwMDqSEHmG1pq+ffsyd+5cpk6dWmCnOZc+fjdTo0YNPvzwQ7Zu3cro0aONjiNEnjJ8+HDmzp3LsGHDCmzRvx0p/AVYZGQk3bt3Z9KkSaxfv97oOELkCVOnTmXKlCn07t2bSZMmGR3HEFL4C7hZs2YRGBhIly5dOHnypNFxhDDU/PnzGTp0KB07duSDDz5AKWV0JEMYVviVUmal1F6l1BqjMrgDHx8foqOjSUlJ4dlnnyU5OdnoSEIY4vPPP6d37960bNmSxYsXYzabjY5kGCNb/AMAWazdBR566CEWLlzIzp07GTp0qNFxhHC5d1e+S6fZnajZqibR0dF4enoaHclQhhR+pZQ/0Br4yIj9u6P27dvTv39/ZsyYwapVq4yOI4TLTFkyhcH7B2N9zMrR4KPsv7Df6EiGM6rFPwN4HbAatH+39Pbbb1O/fn26devGsWPHjI4jhNPNmjWLN+a9AWbABMlpyQVmMZXccHnhV0q1Ac5prW/7sVKlVC+lVLxSKv78+fMuSleweXp68vnnn2M2m+nQoQPXrl0zOpIQTmG1WhkyZAj9+/eniX8TfDx9CtxiKrlhRIv/UaCtUuo4sBxoppRacutGWut5Wut6Wut6pUqVcnXGAqtChQp8+umn7Nu3jwEDBhgdRwiHS0pKolOnTkyfPp2oqCi2fLqFmK4xjA8dT0zXmHy3mIpTaK0NuwAhwJo7bVe3bl0tHGvYsGEa0IsXLzY6ihAOc+HCBd24cWMN6GnTpmmr1Wp0JEMB8TqTmirj+N3U+PHjeeyxx+jTpw+HDh0yOo4Qufbbb7/RqFEjdu3axfLly3n11Vfddpz+nRha+LXWsVrrNkZmcFceHh589tlnFC1alA4dOnDlyhWjIwlx1+Lj4wkODubcuXNs2rSJjh07Gh0pT5MWvxvz8/Pjs88+48iRI/Tu3Tu9+02IfGXdunU89thjeHt7s2PHDpo0aWJ0pDxPCr+ba9asGWPHjmXZsmXMmzfP6DhC5Mi8efNo27YtVapU4fvvv6dq1apGR8oXpPALhg8fTkREBFFRS4mK+l3W6xVO4cg1c7XWjBgxgt69e/P444+zdetWypQpk/sHdhMyH78AYP36v2nd2gutC+HtbWLzZpOs1yscxpFr5iYnJ9O9e3eWLFlCjx49+PDDD2W9iSzIfPzitvbtK25fr9eDpCQr69b9Y3QkUYA4as3cixcv0rJlS5YsWcKECROYN2+eFP27IIVfALY1TL28FCaTFUhm1ar+MtJHOIwj1sz9as9XBHQPYOv/trJ48WJGjBghwzXvkrxVCuD/1zSNjTVhMsUzYsQnPPXUCdauXYuXl5fR8UQ+9/+vL1vRz2k3z7x18+hj6YOurvGq5UXlkMrOiOk2pPCLDMHB6X+QTfHzW0hkZCTPP/88K1askH+nRa79/+sr+1JTU5k0aRJjNo1Bh2owQapOJfZ4rEy9kAvS1SMy1bVrV2bMmMEXX3xBr169ZIy/cLkTJ04QGhrK6NGjCX84XCZacyBpxoksDRgwgMTERMaOHYuvry/Tpk2TPlXhEp9//jm9evXCarXy6aef0qVLFyynLMQejyWkYoi09nNJCr+4rdGjR/PXX3/xzjvvULJkSYYPH250JFGAXblyhQEDBrBw4UIaNGjA0qVLqVSpEgDB5YKl4DuIFH5xW0opZsyYQWJiIiNGjMDX15e+ffsaHUsUQHv27KFTp04cO3aMESNGMHr0aAoVKmR0rAJJCr+4I5PJxMKFC7l48SL9+vWjePHidOrUyehYooCwWq1Mnz6dESNGULp0abZs2cJjjz1mdKwCTU7uimwpVKgQK1asoGnTpnTt2pV169YZHUkUAAkJCURERPD666/zxBNPsH//fin6LiCFX2Sbj48Pq1evJjAwkHbt2rF9+3ajI4l87L///S+BgYHs2LGDefPmER0dTYkSJYyO5Rak8Iscueeee/jmm2+oUKECbdq0Ye/evUZHEvnMtWvXiIqKom3btvj7+/PDDz/Qs2dPGTHmQlL4RY6VKlWKDRs2cO+99xIREcEvv/xidCSRD1hOWRgYPZAaLWrwwQcfMHjwYL7//nuqVKlidDS3I7Nzirt25MgRGjduTOHChdmxYwf+/v5GRxJ51NZftxL+aTgp1hRIg3dqvcOgDoOMjlXgyeycwuECAgL49ttvSUxMpHHj1xg58h+Zy7+Au5s59deuXcszg5+xFX0TmD3NJJVJcl5IcUdS+EWu1KlThylTtnLixAImTvSiWTMtxb+ASp9T/803bV/v9DwfPXqUNm3a0KZNG4r9WQwvDy+ZciGPkMIvcu3ixdqYTOlz+aexYsVZoyMJJ8junPpXrlzhjTfeoEaNGmzbto1p06bxy+Zf2PLiFsaHjiema4x8Atdg8gEukWu2ufxNJCdr0tJSmD+/MxERr9KyZUujowkHSp9TP30VrVvn1Ndas3z5cl577TV+//13IiMjmTx5Mn5+foBMuZCXSItf5Fr6XOvjxyu+/PISDz98gTZt2jBjxgyZ1bMA+f/n+d9LJ+7fv5+QkBCef/55ypQpQ1xcHIsWLcoo+iKP0Vrn+UvdunW1yD+uXLmin376aQ3onj176uvXrxsdSTjJhQsX9Msvv6xNJpMuWbKknjdvnk5NTTU6lrAD4nUmNVVa/MLhihQpQnR0NG+88Qbz588nIiKCCxcuGB1LOFBaWhpz587l4YcfZs6cObz88sscPXqUnj17YjabjY4n7kAKv3AKk8nEpEmTWLx4MXFxcTRs2JCff/7Z6FgilyynLPRZ0odqj1ejT58+VK9enb179zJr1ix8fX2NjieySQq/cKoXXniBzZs3c/HiRRo2bMjGjRuNjiTu0uIti2nyURPmHp3L0UZHGbtwLLGxsQQGBhodTeSQFH7hdI8++ii7du2iXLlytGzZktmzZxsdSeTAzp07adu2LZGjIkkjDUxgKmSiUOVCMr9OPuXywq+UKqeU2qKUOqyUOqSUGuDqDML1KlasSFxcHC1btqRfv35ERUWRmppqdCyRBa01sbGxhIeH07BhQ3bs2EGP5j1k3duCIrMzvs68AH5AHfv3xYBfgGq3u4+M6ik4UlNT9ZAhQzSgw8PDdWJiotGRxA2sVqtet26dbtSokQZ06dKl9dtvv60vXbqktdY67mScnrRtko47GWdwUpEdZDGqx/ChmsDXQPjttpHCX/AsWLBAFypUSAcEBOijR48aHcftpaWl6VWrVuk6depoQJcrV06///77+urVq0ZHE7mQVeE3tI9fKVURqA3szOS2XkqpeKVU/Pnz510dTThZt27d2LhxI+fPn6d+/fq8//6eHE/+Je7OjROtpaamsnTpUmrWrEm7du24dOkSCxYs4NixY/Tr1w8fHx+j4wonMGxaZqVUUWArMFFr/cXttpVpmQuu//3vf4SFjeTEiQUo5YW3t4mYGHXTp0KF41gsEPKChZSyWzCdbEKZlPf4/fdoatSowfDhw+nQoQMeHjKTS0GRp6ZlVkoVAlYBS+9U9EXBVqlSJSIjPwa80NrMtWtpfPDBIZnqwUnmrdtC8nNh6JBRpHWO4Pp95fnyyy/Zv38/nTp1kqLvJowY1aOABcBPWut3XL1/kfe0aOGNj48Zk8mKUiksXdqDli1byspeDqLtI3Q6d+7M4u3dwJwMpjQwJdPu1ZI89dRTmEwystudGPFsPwq8ADRTSu2zX1oZkEPkEemTf02YYGLr1kLMmPEcFouFmjVrMnz4cP755x+jI+ZL586dY+rUqQQEBBAaGsratWt5OqgunuZCKMx4FfIksmmo0TGFETI745vXLjKqx/0kJCTorl27ZowwiY6O1lar1ehYeV5aWpr+9ttvdbt27bSHh4cGdOPGjfUnn3yi//nnH621DMl0J2QxqkfW3BV52nfffUe/fv04cOAA4eHhzJo1i4CAAKNj5Tm///47H3/8MQsWLOD48eOULFmSyMhIevToQdWqVY2OJwyS1cldKfwiz0tNTeXDDz/kzTff5OrVqwwePJiRI0dStGhRo6MZxnLKwuZfN+P5hyfbP9vO2rVrsVqthIWF0bNnT5566im8vLyMjikMJoVf5Htnz55l2LBhLFq0CH9/f6ZPn06HDh3car4YrTVLYpfQbWs3UnUqpEGJNSXo3ao33bt3p1KlSkZHFHmIFH5RYMTFxREVFcXevXsJCwuje/ePOH68IiEhFMjx/1prdu/ezapVq/jiiy84VuYYNANMoDAxLmQsIx8baXRMkQdlVfhl0K7Idxo1asTu3buZO3cuQ4d+RUzMfUAaXl6wcaOmSZP8/7JOS0sjLi4uo9ifOnUKDw8PwsLCaFC7O0vTxoFORls9ue9qmNFxRT4jg3dFvmQ2m3n55ZcZMOALlPICzFy/bqVVq7cYNGgQ8fHx+e5DYCkpKWzcuJE+ffpQtmxZmjZtypw5c6hduzaffPIJ586d45tvvqH6PcMwfRoDW8Zj+jSGC/sK4L85wqmkq0fkaxYLhIVBcrLGbE4jOPhNLJZ3SE5OJiAggC5dutC5c2ceeOABo6PexHLKQuzxWILvD+by4cusWrWK1atXk5iYSJEiRWjdujXPPPMMrVq1olixYjffN+N3Bk/Pfy98LkQ66eMXBZbFArGxZPTxJyYmsmrVKpYsWcLWrVsB22IwnTt35tlnn6VkyZKG5v16z9c8u/ZZkq3JkAp8AsWvFKdt27Y888wzPP7443ecHO3W31mIzEjhF27pxIkTfPbZZ3z66accPnyYQoUK0bJlS7p06UKbNm3Yt8/H6QX077//JjY2lk2bNrFp0yaOlDqScXIWrXixfCRzu8zF09PTOQGE25LCL9ya1pr9+/ezdOlSli1bxh9//EHhwmFcv74WrQtRqJBm5cpEWrcuket5a5KSkrBYLBmFPj4+HqvVSpEiRXjssccoHFCf6MJvgSkZrJ7MfTSGXi2l2S4cTwq/EHZpaWnExsby+uuJ/PDDU9gGt6UAo/D0fIfy5ctTvnx5KlSokPE1/fty5cplfDBq3noLq/bE8kztpjzi55NR6Ldv305SUhJms5kGDRrQvHlzmjdvToMGDfD09GTyZBg5x4K1fCymkyFM6BPMG28YeEBEgSWFX4hb2E6SapKTwWy2EhX1JWbzLk6ePMmJEyc4efIkCQkJ/xodVKZMGbwfqsfxpjG2mS7TPOGTB+H0IapXr55R6Js2bco999yTxX7l5KxwPhnHL8QtbLOCKnsfv5ng4PZA+5u2uX79OqdPn77pzeDEiRN8deHv/5/eWCdT84lQvn1zI35+ftncr5ycFcaRwi/cm78FGseCfwjw7wrs5eVFpUqV/jUVQoP1FnrvWA/a1k8f9cTz2Sr66YKDpeAL40jhF27LcspC2OIwktOS8TR7EtM1huBy2avGtpOxMazaE0u7uiFyclbkK1L4hduKPR5LcloyaTqN5LRk2weqsln4wVb8peCL/EimbBD5nuWUhcnbJ2M5ZcnR/UIqhuBp9sSszHiaPQmpGOKcgELkMdLiF/labrprgssFE9M1htjjsYRUDMlRa1+I/EwKv8gT0ueuyWkBzm13TXC5YCn4wu1I4ReGy02rPb27Jv2+0l0jxJ1J4RcOcbctdshdq126a4TIOSn8Itdy02KH3LfapbtGiJyRwi8yGNnPLq12IVxHCn8Bc7fF2+h+dmm1C+E6UvjzICOKt/SzC+E+pPBnITcnK3N7XyOKt/SzC+E+CnThN6LlnNsTnUYVb2m1C+E+DCn8SqkWwEzADHyktZ7i6H0Y1XLO7YlOI4u3tNqFcA8uL/xKKTPwARAOnAZ2K6VWa60PO3I/RrWcHdFlIsVbCOFMRrT46wPHtNa/AiillgNPAg4t/Ea1nB3RZSLFWwjhTC5felEp1R5oobXuYf/5BaCB1jrqlu16Ab0AypcvX/fEiRM53lduTrIKIUR+l++WXtRazwPmgW3N3bt5DGk5CyHEvxkxH//vQLkbfva3XyeEEMIFjCj8u4GHlFIPKKU8geeA1QbkEEIIt+Tyrh6tdapSKgr4FttwzoVa60OuziGEEO7KkD5+rfU6YJ0R+xZCCHcna+4KIYSbkcIvhBBuxuXj+O+GUuo8kPOB/Db/Af50YBxHkVw5I7lyRnLlTEHNVUFrXerWK/NF4c8NpVR8Zh9gMJrkyhnJlTOSK2fcLZd09QghhJuRwi+EEG7GHQr/PKMDZEFy5YzkyhnJlTNulavA9/ELIYS4mTu0+IUQQtxACr8QQriZAlH4lVIdlFKHlFJWpVS9W257Qyl1TCl1RCkVkcX9H1BK7bRvt8I+eZyjM65QSu2zX44rpfZlsd1xpdSP9u3iHZ0jk/2NUUr9fkO2Vlls18J+DI8ppYa5INfbSqmflVIHlFJfKqWKZ7GdS47XnX5/pZSX/Tk+Zn8tVXRWlhv2WU4ptUUpddj++h+QyTYhSqmLNzy/o5ydy77f2z4vyuY9+/E6oJSq44JMATcch31KqUtKqYG3bOOS46WUWqiUOqeUOnjDdSWUUhuVUkftX32zuG+kfZujSqnIuwqgtc73F6AqEADEAvVuuL4asB/wAh4A/geYM7n/58Bz9u/nAH2dnHc6MCqL244D/3HhsRsDDLnDNmb7sXsQ8LQf02pOzvU44GH//i3gLaOOV3Z+f+BlYI79++eAFS547vyAOvbviwG/ZJIrBFjjqtdTdp8XoBWwHlBAQ2Cni/OZgTPYPuDk8uMFNAXqAAdvuG4qMMz+/bDMXvNACeBX+1df+/e+Od1/gWjxa61/0lofyeSmJ4HlWuvrWuvfgGPYln7MoJRSQDMg2n7VJ8BTzspq39+zwGfO2ocTZCyXqbVOBtKXy3QarfUGrXWq/cfvsa3bYJTs/P5PYnvtgO21FGZ/rp1Ga52gtf7B/v1l4CegrDP36UBPAou1zfdAcaWUnwv3Hwb8T2t9tzMC5IrWehvw1y1X3/gayqoORQAbtdZ/aa0TgY1Ai5zuv0AU/tsoC5y64efT/PsPoyTw9w1FJrNtHKkJcFZrfTSL2zWwQSm1x778pCtE2f/dXpjFv5fZOY7O1A1b6zAzrjhe2fn9M7axv5YuYnttuYS9a6k2sDOTm4OVUvuVUuuVUtVdFOlOz4vRr6nnyLrxZcTxAiittU6wf38GKJ3JNg45bnl26cVbKaU2AWUyuWmE1vprV+fJTDYzduL2rf3GWuvflVL3ARuVUj/bWwdOyQV8CIzH9oc6Hls3VLfc7M8RudKPl1JqBJAKLM3iYRx+vPIbpVRRYBUwUGt96Zabf8DWnXHFfv7mK+AhF8TKs8+L/RxeW+CNTG426njdRGutlVJOG2ufbwq/1rr5XdwtO8s8XsD2b6aHvaV210tB3imjUsoDeAaoe5vH+N3+9ZxS6kts3Qy5+oPJ7rFTSs0H1mRyk1OWy8zG8XoRaAOEaXsHZyaP4fDjlYns/P7p25y2P8/3YnttOZVSqhC2or9Ua/3Frbff+EagtV6nlJqtlPqP1tqpE5Jl43kxcgnWlsAPWuuzt95g1PGyO6uU8tNaJ9i7vc5lss3v2M5DpPPHdm4zRwp6V89q4Dn7iIsHsL1z77pxA3tB2QK0t18VCTjrP4jmwM9a69OZ3aiUKqKUKpb+PbYTnAcz29ZRbulXfTqL/bl8uUylVAvgdaCt1vpqFtu46nhl5/dfje21A7bX0uas3qwcxX4OYQHwk9b6nSy2KZN+rkEpVR/b37xT35Cy+bysBrraR/c0BC7e0M3hbFn+123E8brBja+hrOrQt8DjSilfe7fs4/brcsbZZ69dccFWsE4D14GzwLc33DYC24iMI0DLG65fB9xv//5BbG8Ix4CVgJeTci4C+txy3f3Auhty7LdfDmHr8nD2sfsU+BE4YH/h+d2ay/5zK2yjRv7nolzHsPVl7rNf5tyay5XHK7PfHxiH7Y0JwNv+2jlmfy096IJj1BhbF92BG45TK6BP+usMiLIfm/3YTpI3ckGuTJ+XW3Ip4AP78fyRG0bjOTlbEWyF/N4brnP58cL2xpMApNhrV3ds54RigKPAJqCEfdt6wEc33Leb/XV2DHjpbvYvUzYIIYSbKehdPUIIIW4hhV8IIdyMFH4hhHAzUviFEMLNSOEXQgg3I4VfCCHcjBR+IYRwM1L4hbgLSqlH7BPbeds/qXpIKVXD6FxCZId8gEuIu6SUmoDtE7s+wGmt9WSDIwmRLVL4hbhL9nl7dgNJ2D7an2ZwJCGyRbp6hLh7JYGi2Fa/8jY4ixDZJi1+Ie6SUmo1ttW4HsA2uV2UwZGEyJZ8Mx+/EHmJUqorkKK1XqaUMgNxSqlmWuvNRmcT4k6kxS+EEG5G+viFEMLNSOEXQgg3I4VfCCHcjBR+IYRwM1L4hRDCzUjhF0IINyOFXwgh3Mz/AV7jZsKPGPSeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "TINYDOUBLE = 1.0e-100\n",
    "def coord_leq_bound(x,xstar):\n",
    "    # Returns 1.0 if x <= xstar, 0.0 otherwise.\n",
    "    # Requires appropriately defined TINYDOUBLE\n",
    "    return min_noif(x-xstar-TINYDOUBLE,0.0)/(x-xstar-TINYDOUBLE)\n",
    "\n",
    "def coord_geq_bound(x,xstar):\n",
    "    # Returns 1.0 if x >= xstar, 0.0 otherwise.\n",
    "    # Requires appropriately defined TINYDOUBLE\n",
    "    return max_noif(x-xstar+TINYDOUBLE,0.0)/(x-xstar+TINYDOUBLE)\n",
    "\n",
    "def coord_less_bound(x,xstar):\n",
    "    # Returns 1.0 if x < xstar, 0.0 otherwise.\n",
    "    # Requires appropriately defined TINYDOUBLE\n",
    "    return min_noif(x-xstar,0.0)/(x-xstar-TINYDOUBLE)\n",
    "\n",
    "def coord_greater_bound(x,xstar):\n",
    "    # Returns 1.0 if x > xstar, 0.0 otherwise.\n",
    "    # Requires appropriately defined TINYDOUBLE\n",
    "    return max_noif(x-xstar,0.0)/(x-xstar+TINYDOUBLE)\n",
    "\n",
    "# Now, define our the equation and plot it.\n",
    "x_data = np.arange(start = -10.0, stop = 11.0, step = 1.0)\n",
    "y_data = coord_less_bound(x_data,0.0)*(0.1*x_data**2.0+1.0)\\\n",
    "        +coord_geq_bound(x_data,0.0)*np.exp(x_data/5.0)\n",
    "\n",
    "plt.figure()\n",
    "a = plt.plot(x_data,y_data,'k',label=\"Piecewise function\")\n",
    "b = plt.plot(x_data,0.1*x_data**2.0+1.0,'b.',label=\"y=0.1*x^2+1\")\n",
    "c = plt.plot(x_data,np.exp(x_data/5.0),'g.',label=\"y=exp(x/5)\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above shows the expected piecewise-defined function. It is important in applying these functions that each greater-than be paired with a less-than-or-equal-to, or vice versa. Otherwise, the way these are written, a point on the boundary will be set to zero or twice the expected value. \n",
    "\n",
    "These functions can be easily combined for more complicated piecewise-defined functions; if a piece of a function is defined as $f(x)$ on $x^*_- \\leq x < x^*_+$, for instance, simply multiply by both functions, e.g. \n",
    "```\n",
    "coord_geq_bound(x,x_star_minus)*coord_less_bound(x,x_star_plus)*f(x)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sympy'></a>\n",
    "\n",
    "# Step 3: Rewrite functions to work with symbolic expressions \\[Back to [top](#toc)\\]\n",
    "$$\\label{sympy}$$\n",
    "\n",
    "In order to use this with sympy expressions in NRPy+, we will need to rewrite the `min` and `max` functions with slightly different syntax. Critically, we will change `0.5` to `sp.Rational(1,2)` and calls to `np.absolute()` to `nrpyAbs()`. We will also need to import `outputC.py` here for access to `nrpyAbs()`. The other functions will not require redefinition, because they only call specific combinations of the `min` and `max` function. \n",
    "\n",
    "In practice, we want to use `nrpyAbs()` and *not* `sp.Abs()` with our symbolic expressions, which will force `outputC` to use the C function `fabs()`, and not try to multiply the argument by its complex conjugate and then take the square root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from outputC import nrpyAbs     # NRPy+: Core C code output module\n",
    "\n",
    "def min_noif(a,b):\n",
    "    # Returns the minimum of a and b\n",
    "    if a==sp.sympify(0):\n",
    "        return sp.Rational(1,2) * (b-nrpyAbs(b))\n",
    "    if b==sp.sympify(0):\n",
    "        return sp.Rational(1,2) * (a-nrpyAbs(a))\n",
    "    return sp.Rational(1,2) * (a+b-nrpyAbs(a-b))\n",
    "\n",
    "def max_noif(a,b):\n",
    "    # Returns the maximum of a and b\n",
    "    if a==sp.sympify(0):\n",
    "        return sp.Rational(1,2) * (b+nrpyAbs(b))\n",
    "    if b==sp.sympify(0):\n",
    "        return sp.Rational(1,2) * (a+nrpyAbs(a))\n",
    "    return sp.Rational(1,2) * (a+b+nrpyAbs(a-b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='validation'></a>\n",
    "\n",
    "# Step 4: Validation against `Min_Max_and_Piecewise_Expressions` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{validation}$$\n",
    "\n",
    "As a code validation check, we will verify agreement in the SymPy expressions for plane-wave initial data for the Scalar Wave equation between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [Min_Max_and_Piecewise_Expressions](../edit/Min_Max_and_Piecewise_Expressions.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ALL TESTS PASSED!\n"
     ]
    }
   ],
   "source": [
    "# Reset & redefine TINYDOUBLE for proper comparison\n",
    "%reset_selective -f TINYDOUBLE\n",
    "\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par   # NRPy+: parameter interface\n",
    "TINYDOUBLE = par.Cparameters(\"REAL\", thismodule, \"TINYDOUBLE\", 1e-100)\n",
    "\n",
    "import Min_Max_and_Piecewise_Expressions as noif\n",
    "all_passed=0\n",
    "\n",
    "def comp_func(expr1,expr2,basename,prefixname2=\"noif.\"):\n",
    "    passed = 0\n",
    "    if str(expr1-expr2)!=\"0\":\n",
    "        print(basename+\" - \"+prefixname2+basename+\" = \"+ str(expr1-expr2))\n",
    "        passed = 1\n",
    "    return passed\n",
    "\n",
    "a,b = sp.symbols(\"a b\")\n",
    "\n",
    "here = min_noif(a,b)\n",
    "there = noif.min_noif(a,b)\n",
    "all_passed += comp_func(here,there,\"min_noif\")\n",
    "\n",
    "here = max_noif(a,b)\n",
    "there = noif.max_noif(a,b)\n",
    "all_passed += comp_func(here,there,\"max_noif\")\n",
    "\n",
    "here = coord_leq_bound(a,b)\n",
    "there = noif.coord_leq_bound(a,b)\n",
    "all_passed += comp_func(here,there,\"coord_leq_bound\")\n",
    "\n",
    "here = coord_geq_bound(a,b)\n",
    "there = noif.coord_geq_bound(a,b)\n",
    "all_passed += comp_func(here,there,\"coord_geq_bound\")\n",
    "\n",
    "here = coord_less_bound(a,b)\n",
    "there = noif.coord_less_bound(a,b)\n",
    "all_passed += comp_func(here,there,\"coord_less_bound\")\n",
    "\n",
    "here = coord_greater_bound(a,b)\n",
    "there = noif.coord_greater_bound(a,b)\n",
    "all_passed += comp_func(here,there,\"coord_greater_bound\")\n",
    "\n",
    "import sys\n",
    "if all_passed==0:\n",
    "    print(\"ALL TESTS PASSED!\")\n",
    "else:\n",
    "    print(\"ERROR: AT LEAST ONE TEST DID NOT PASS\")\n",
    "    sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-Min_Max_and_Piecewise_Expressions.pdf](Tutorial-Min_Max_and_Piecewise_Expressions.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Min_Max_and_Piecewise_Expressions.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-Min_Max_and_Piecewise_Expressions.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Min_Max_and_Piecewise_Expressions\")"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}