Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Uncertainty Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p align=\"center\">\n",
    "  <img src=\"./figures/numerical_analysis.png\" alt=\"Image description\" />\n",
    "</p>\n",
    "\n",
    "**Numerical analysis** is the study of algorithms that use numerical approximation for problems involving mathematical analysis (Trefethen, 1992). Numerical analysis primarily involves studying numerical uncertainty, validating a program’s results, and the optimisation of operations. \n",
    "\n",
    "We focus on **numerical uncertainty** in order to get insight on the detection and quantification of numerical errors and their origin. The two methods of investigating numerical uncertainty are static and dynamic analysis. Static analysis analyses a program without executing it, thus it depends greatly on the algorithmic analysis of the program. While effective, static analysis cannot be scaled up to large code bases, which limits its scope. Dynamic analysis analyses a program during or after execution and while efficient, it is not generally as thorough as static analysis and is dependent on the program being analysed and its input data. Nonetheless, it is scalable and can be applied efficiently on large code bases.\n",
    "\n",
    "**Stochastic arithmetic**, a form of dynamic analysis, analyses a program’s execution to measure numerical error by executing the program multiple times with random numerical perturbations.\n",
    "Numerical uncertainty is determined by quantifying the errors that can occur in floating point arithmetic, these errors typically being rounding and truncation errors. \n",
    "More specifically, one type of rounding error, absorption errors, are caused by mathematical operations between large and small numbers where the shifting of decimal points in the mantissa causes a loss of precision. In turn, another rounding error, catastrophic cancellation, is the loss of leading significant digits caused by subtracting two approximately equal operands.\n",
    "Two stochastic arithmetic methods that compute numerical uncertainty are DSA (Discrete Stochastic Arithmetic) (Vignes, 2004) and MCA (Monte Carlo Arithmetic) (Parker, 1997). Both techniques leverage randomness to model loss of accuracy inherent to floating point arithmetic due to its finite precision. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**MCA** is a variant of stochastic arithmetic that leverages randomness in fundamental floating-point operations, introducing a probabilistic element. We use it over DSA due to its scalability to large code bases. It simulates computations for a given virtual precision using the following perturbation: \n",
    "\n",
    "$$ inexact(x)=x+2^{e_{x}-t}\\xi $$\n",
    "\n",
    "where \n",
    "$$\n",
    "e_{x} \\textnormal{ is the exponent in the floating point representation of } x \\\\\n",
    "        t \\textnormal{ is the virtual precision and } \\\\\n",
    "        \\xi \\textnormal{ is a random uniform variable of } (-\\frac{1}{2},\\frac{1}{2})\n",
    "$$\n",
    "\n",
    "\n",
    "MCA has three modes that allow perturbation to be introduced in the function input (Precision Bounding---PB), the function output (Random Rounding---RR) or both (full $mca\\_mode$).  Full $mca\\_mode$ combines RR and PB in the following equation: \n",
    "\n",
    "$$mca\\_mode(x \\circ y)=inexact_{RR}(inexact_{PB}(x) \\circ inexact_{PB}(y))$$\n",
    "\n",
    "Random Rounding tracks rounding errors, while Precision Bounding tracks catastrophic cancellations, and full $mca\\_mode$ is a combination of both approaches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Verificarlo** (Denis et al., 2015) is a clang-based compiler that replaces floating point operations by a generic call to one of several configurable floating point models available. The source code is available on [Github](https://github.com/verificarlo/verificarlo), with libraries compiled with Verificarlo available on Dockerhub.\n",
    "\n",
    "\n",
    "<p align=\"center\">\n",
    "  <img src=\"./figures/verificarlo.png\" alt=\"Image description\" />\n",
    "</p>\n",
    "\n",
    "Verificarlo requires recompilation of a program in order to instrument the floating point operations, which is not always feasible on large, complex code bases such as Tensorflow or PyTorch.\n",
    "This is especially difficult when the clang compiler used by Verificarlo is incompatible with the one used by the program or if the program calls upon many third-party tools which would also require recompilation. \n",
    "\\\n",
    "Nonetheless, MCA can be leveraged in programs through the use of the **Fuzzy ecosystem** (Kiar et al., 2020), a collection of Python software packages compiled with Verificarlo.\n",
    "Fuzzy enables the instrumentation of different libraries by Verificarlo in order to quantify the numerical uncertainty present in code using these libraries.\n",
    "The [Fuzzy ecosystem](https://github.com/verificarlo/fuzzy) with instrumented packages available such as the Python interpreter, the Numpy and SciPy packages and the libmath library.\n",
    "\n",
    "\n",
    "<p align=\"center\">\n",
    "  <img src=\"./figures/tool_differences.png\" alt=\"Image description\" />\n",
    "</p>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Verrou** (Févotte et al., 2016) is a tool which also uses Monte Carlo Arithmetic to monitor the accuracy of floating point operations without needing to instrument the source code or recompile it. \n",
    "The Verrou tool is located at [Github](https://github.com/edf-hpc/verrou) and its documentation [here](https://edf-hpc.github.io/verrou/vr-manual.html#vr-cr.start-instrumentation).\n",
    "\n",
    "Verrou is based upon Valgrind (Nethercote et al., 2007), a framework allowing for plug-in tools which use dynamic binary instrumentation (DBI), i.e. the process of modifying instructions of a binary program while it executes.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A comparison of the most popular software implementations of MCA (Verificarlo and Verrou; discussed above) and DSA (CADNA; not discussed in this tutorial, consult (Jezequel et al., 2013)) is shown in the table below:\n",
    "\n",
    "|                          | **Verificarlo** | **Verrou** | **CADNA** |\n",
    "|:-------------------------------:|:---------------:|:----------:|:---------:|\n",
    "| **Theoretical Approach**        |      MCA        |     MCA    |    DSA    |\n",
    "| **Practical Approach**          |    Compiler     |     DBI    | Manual Implementation |\n",
    "| **Languages**                   |  C, C++, Fortran| All        | C, C++, Fortran |\n",
    "| **Perturbation Type**           | RR, PB, full $mca\\_mode$| RR, PB, full $mca\\_mode$ | RR |\n",
    "| **Parallelization**             | Yes             | No         | Conditionally upon Synchronisation |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practical Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import significantdigits as sd\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verificarlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before beginning, please pull and open up an interactive container of the `verificarlo/fuzzy:nilearn` Docker image. \\\n",
    "When doing so, bind the local `data` directory in this repository to a directory called `/data` in the Docker image, to properly save the files. \\\n",
    "For more practical examples, please consult this [tutorial](https://github.com/verificarlo/fuzzy/tree/master/tutorial) on the Fuzzy repository and this [tutorial](https://github.com/verificarlo/verificarlo_tutorial) on the Verificarlo repository."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quadratic equation below adapted from the Parker MCA paper can be a suprising example of instability.\n",
    "$$ 7169x^2 - 8686x + 2631 = 0 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Roots of the quadratic equation are:\n",
      "[0.60536166] [0.60624387]\n"
     ]
    }
   ],
   "source": [
    "#Sample implementation\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "# Define the quadratic equation as a function\n",
    "def quadratic_eq(x):\n",
    "    return 7169 * x**2 - 8686 * x + 2631\n",
    "\n",
    "# Use fsolve to find roots, providing initial guesses\n",
    "root1 = fsolve(quadratic_eq, 0)  # First guess near zero\n",
    "root2 = fsolve(quadratic_eq, 1)  # Another guess near one\n",
    "\n",
    "print(\"Roots of the quadratic equation are:\")\n",
    "print(root1, root2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can run the above code located in [```quadratic_eq.py```](https://github.com/InesGP/fuzzy-tutorial/blob/main/scripts/quadratic_eq.py) file in the Fuzzy Docker image and then execute the code cell below to compare your instability to that from the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Execute this cell with your MCA instability\n",
    "\"\"\"\n",
    "\n",
    "# Specify the file path to the .npy file\n",
    "file_path = './data/mca_quadratic_results.npy'\n",
    "# Load the .npy file\n",
    "data = np.load(file_path)\n",
    "print(\"Roots of the quadratic equation are:\")\n",
    "print(\"r1\\t\\t\\t\\tr2\")\n",
    "for r1, r2 in data:\n",
    "    print(f\"{r1:<10.16f} {r2:<10.16f}\")\n",
    "print()\n",
    "# Inspect the loaded data\n",
    "print(\"Computed average: \", np.mean(data,  dtype=np.float64, axis=0))\n",
    "print(\"Computed standard deviation: \", np.std(data, dtype=np.float64, axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Roots of the quadratic equation are:\n",
      "r1\t\t\t\tr2\n",
      "0.6053616574595290 0.6062438663216425\n",
      "0.6053616574596200 0.6062438663217239\n",
      "0.6053616574596471 0.6062438663217581\n",
      "0.6053616574596571 0.6062438663217253\n",
      "0.6053616574595766 0.6062438663216633\n",
      "0.6053616574594421 0.6062438663216593\n",
      "0.6053616574595408 0.6062438663217650\n",
      "0.6053616574596802 0.6062438663216190\n",
      "0.6053616574596514 0.6062438663218352\n",
      "0.6053616574596260 0.6062438663218589\n",
      "\n",
      "Computed average:  [0.60536166 0.60624387]\n",
      "Computed standard deviation:  [7.03905826e-14 7.67795601e-14]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Compare this cell with previously generated instability to your own just above\n",
    "\"\"\"\n",
    "\n",
    "# Specify the file path to the .npy file\n",
    "file_path = './data/previous_mca_quadratic_results.npy'\n",
    "# Load the .npy file\n",
    "data = np.load(file_path)\n",
    "print(\"Roots of the quadratic equation are:\")\n",
    "print(\"r1\\t\\t\\t\\tr2\")\n",
    "for r1, r2 in data:\n",
    "    print(f\"{r1:<10.16f} {r2:<10.16f}\")\n",
    "print()\n",
    "# Inspect the loaded data\n",
    "print(\"Computed average: \", np.mean(data,  dtype=np.float64, axis=0))\n",
    "print(\"Computed standard deviation: \", np.std(data, dtype=np.float64, axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another adapted example from the MCA paper demonstrates the impact a machine's floating point precision can have. \\\n",
    "Here, we have an iterative formula where the convergence point changes depending on the numerical precision it is executed in. \n",
    "\n",
    "$$\n",
    "x_0 = 1.5100050721319 \\\\\n",
    "x_{k+1} = (3x_k^4 - 20x_k^3 + 35_k^2 - 24) / (4x_k^3 - 30x_k^2 + 70x_k - 50)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The iterative formula converges to: 0.37228132326901436\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3.0000000000000426"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Sample implementation \n",
    "\n",
    "# Define the iterative formula\n",
    "def iterative_formula(x_k):\n",
    "    numerator = 3 * x_k**4 - 20 * x_k**3 + 35 * x_k**2 - 24\n",
    "    denominator = 4 * x_k**3 - 30 * x_k**2 + 70 * x_k - 50\n",
    "    return numerator / denominator\n",
    "\n",
    "# Initialize variables\n",
    "x_k = 1.5100050721319  # Initial guess\n",
    "iterations = 30\n",
    "\n",
    "# Perform iterations\n",
    "results = [x_k]  # To store all iteration results\n",
    "for _ in range(iterations):\n",
    "    x_k = iterative_formula(x_k)\n",
    "    results.append(x_k)\n",
    "\n",
    "results[-1]  # The value after 30 iterations\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can run the above code located in [`convergence_eq.py`](https://github.com/InesGP/fuzzy-tutorial/blob/main/scripts/convergence_eq.py) file in the Fuzzy Docker image and lower the virtual precision to single precision in order to emphasize the instability that can occur. Double precision will still have instability, but it will be less noticeable.\n",
    "To modify the virtual precision of the script, in the Docker image, run this command:\n",
    "```\n",
    "for i in {1..10}; do VFC_BACKENDS=\"libinterflop_mca.so -m rr --precision-binary32=24 --precision-binary64=53\" python convergence_eq.py; done\n",
    "```\n",
    "After, you can execute the code cell below and compare your instability to that from the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Execute this cell with your MCA instability\n",
    "\"\"\"\n",
    "\n",
    "f = open('./data/mca_convergence_results.txt', 'r')\n",
    "data = []\n",
    "print(\"Convergence values after 30 iterations: \")\n",
    "for value in f:\n",
    "    i = np.float64(value.strip())\n",
    "    data.append(i)\n",
    "    print(i)\n",
    "print()\n",
    "# Inspect the loaded data\n",
    "print(\"Computed average: \", np.mean(data,  dtype=np.float64, axis=0))\n",
    "print(\"Computed standard deviation: \", np.std(data, dtype=np.float64, axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Convergence values after 30 iterations: \n",
      "3.0\n",
      "3.000000000000071\n",
      "2.999999999999943\n",
      "3.0000000000000284\n",
      "1.0\n",
      "2.999999999999929\n",
      "3.0000000000001137\n",
      "2.9999999999999147\n",
      "2.9999999999999716\n",
      "3.000000000000071\n",
      "\n",
      "Computed average:  2.8000000000000043\n",
      "Computed standard deviation:  0.6000000000000014\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Compare this cell with previously generated instability to your own just above\n",
    "\"\"\"\n",
    "\n",
    "f = open('./data/previous_mca_convergence_results.txt', 'r')\n",
    "data = []\n",
    "print(\"Convergence values after 30 iterations: \")\n",
    "for value in f:\n",
    "    i = np.float64(value.strip())\n",
    "    data.append(i)\n",
    "    print(i)\n",
    "print()\n",
    "# Inspect the loaded data\n",
    "print(\"Computed average: \", np.mean(data,  dtype=np.float64, axis=0))\n",
    "print(\"Computed standard deviation: \", np.std(data, dtype=np.float64, axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compute the numerical stability of a [Hilbert matrix](https://en.m.wikipedia.org/wiki/Hilbert_matrix). The Hilbert matrices are considered ill-conditioned matrices and are famously difficult to use in numerical computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hilbert Matrix:\n",
      "[[1.         0.5        0.33333333 0.25       0.2       ]\n",
      " [0.5        0.33333333 0.25       0.2        0.16666667]\n",
      " [0.33333333 0.25       0.2        0.16666667 0.14285714]\n",
      " [0.25       0.2        0.16666667 0.14285714 0.125     ]\n",
      " [0.2        0.16666667 0.14285714 0.125      0.11111111]]\n",
      "\n",
      "Inverse of Hilbert Matrix:\n",
      "[[ 2.500e+01 -3.000e+02  1.050e+03 -1.400e+03  6.300e+02]\n",
      " [-3.000e+02  4.800e+03 -1.890e+04  2.688e+04 -1.260e+04]\n",
      " [ 1.050e+03 -1.890e+04  7.938e+04 -1.176e+05  5.670e+04]\n",
      " [-1.400e+03  2.688e+04 -1.176e+05  1.792e+05 -8.820e+04]\n",
      " [ 6.300e+02 -1.260e+04  5.670e+04 -8.820e+04  4.410e+04]]\n"
     ]
    }
   ],
   "source": [
    "from scipy.linalg import inv\n",
    "\n",
    "# Function to generate a Hilbert matrix of size n x n\n",
    "def hilbert_matrix(n):\n",
    "    return np.fromfunction(lambda i, j: 1.0 / (i + j + 1), (n, n))\n",
    "\n",
    "# Define the size of the matrix\n",
    "n = 5  # Example size (can change as needed)\n",
    "\n",
    "# Generate the Hilbert matrix\n",
    "H = hilbert_matrix(n)\n",
    "\n",
    "# Calculate the inverse using scipy.linalg.inv\n",
    "H_inv = inv(H)\n",
    "\n",
    "# Display the Hilbert matrix and its inverse\n",
    "print(\"Hilbert Matrix:\")\n",
    "print(H)\n",
    "\n",
    "print(\"\\nInverse of Hilbert Matrix:\")\n",
    "print(H_inv)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can run the above code located in [`hilbert_eq.py`](https://github.com/InesGP/fuzzy-tutorial/blob/main/scripts/hilbert_eq.py) file in the Fuzzy Docker image, and then execute the code cell below to observe the instability in this equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Execute this cell with your MCA instability\n",
    "\"\"\"\n",
    "\n",
    "# Specify the file path to the .npy file\n",
    "file_path = './data/mca_hilbert_results.npy'\n",
    "# Load the .npy file\n",
    "data = np.load(file_path)\n",
    "\n",
    "# Inspect the loaded data\n",
    "print(\"Computed average: \", np.mean(data,  dtype=np.float64, axis=0))\n",
    "print(\"Computed standard deviation: \", np.std(data, dtype=np.float64, axis=0))\n",
    "\n",
    "\n",
    "sigdig=sd.significant_digits(data , reference=np.mean(data, axis=0, dtype=np.float64), basis=10)\n",
    "gc = sns.heatmap(sigdig, cmap='RdYlGn', vmax=14, )\n",
    "gc.set_title(\"Significant Digit Heatmap for Hilbert Matrix\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed average of values: \n",
      " [[ 2.500e+01 -3.000e+02  1.050e+03 -1.400e+03  6.300e+02]\n",
      " [-3.000e+02  4.800e+03 -1.890e+04  2.688e+04 -1.260e+04]\n",
      " [ 1.050e+03 -1.890e+04  7.938e+04 -1.176e+05  5.670e+04]\n",
      " [-1.400e+03  2.688e+04 -1.176e+05  1.792e+05 -8.820e+04]\n",
      " [ 6.300e+02 -1.260e+04  5.670e+04 -8.820e+04  4.410e+04]]\n",
      "Computed standard deviation of values: \n",
      " [[2.02083719e-11 3.80836269e-10 1.64921341e-09 2.49758091e-09\n",
      "  1.22432262e-09]\n",
      " [3.86034906e-10 7.27126768e-09 3.14947527e-08 4.77026912e-08\n",
      "  2.33850690e-08]\n",
      " [1.69050749e-09 3.18448282e-08 1.37936768e-07 2.08922775e-07\n",
      "  1.02421404e-07]\n",
      " [2.57602153e-09 4.86036840e-08 2.10565624e-07 3.18932901e-07\n",
      "  1.56356269e-07]\n",
      " [1.27087175e-09 2.39665129e-08 1.03827982e-07 1.57264195e-07\n",
      "  7.70972591e-08]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAGzCAYAAABD8k8yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8MElEQVR4nO3de1yUZf7/8feAMigK5hFPaOoqqRsUHsLUPCWZ6Wp5SNsCNXe3zK1IK9pNsdqw3DUtzdRM7WC65iHX3TBPZT48S5QdNI+YJih5RgVk7t8f/Zivc9/AMDA0YK/n43E/dK77nms+cw/DfPhc13WPzTAMQwAAAEXw83UAAACg/CNhAAAAbpEwAAAAt0gYAACAWyQMAADALRIGAADgFgkDAABwi4QBAAC4RcIAAADc+s0nDE2bNlVcXFyZPsZnn30mm82mzz77zKX9vffeU3h4uCpXrqwaNWpIkrp166Zu3bqVaTzlRWmea1xcnJo2berVeFAyU6ZMUbNmzeTv76/IyEhfh+PiyJEjstlsWrBggbMtMTFRNpvN5TibzabHHnvsV46uYuM9+Ntz3SYMe/bs0aBBg9SkSRMFBgaqYcOGuvPOO/XGG2/4OjRJ0t69exUXF6fmzZtr7ty5mjNnjq9DcrFlyxYlJibq7NmzxTo+Li5ONpvNuVWrVk3NmjXToEGDtGzZMjkcjjKN99KlS0pMTLQkZYXJT+I++uijAvfHxcWpWrVqXozQytNzXB59+umnevrpp3X77bdr/vz5evnll8v08dy9LhXhg/+7775TYmKijhw5Uqzj8xMcPz8//fjjj5b958+fV5UqVUr83D197+C3q5KvAygLW7ZsUffu3RUWFqbRo0crNDRUP/74o7Zt26bp06dr7NixzmP37dsnP7+yzZu6du2qy5cvKyAgwNn22WefyeFwaPr06WrRooWz/dNPPy3TWIpry5YtmjRpkuLi4pzVD3fsdrvefvttSdLly5eVlpam//znPxo0aJC6deumjz/+WMHBwc7jS/Nc586d65KEXLp0SZMmTZKkClOhKck5Lm82bNggPz8/zZs3z+Xnu7xo0qSJLl++rMqVK/s6FKfvvvtOkyZNUrdu3Tz6C91ut+vDDz/U008/7dK+fPnyUsVT0veO+T2I6991mTD84x//UEhIiHbu3Gn5RXzy5EmX23a7vczj8fPzU2BgYIFxmOMrj790i6tSpUr64x//6NL20ksvafLkyUpISNDo0aO1ZMkS577SPNfy9AHwW3by5ElVqVLFaz+3hmHoypUrqlKlilf6s9lslveer1y5cqVU5+nuu+8uMGFYtGiR+vbtq2XLlpU2xGLJyspSUFAQ78HfoOtySOLgwYNq06ZNgX+11a1b1+V2QXMYvv76a91xxx2qUqWKGjVqpJdeeknz58+XzWZzKSM2bdpU99xzjzZv3qwOHTooMDBQzZo107vvvuvSn3kOQ9OmTTVx4kRJUp06dWSz2ZSYmCip4HH9K1euKDExUS1btlRgYKDq16+ve++9VwcPHnQe889//lOdOnVSrVq1VKVKFUVFRRVYbs8vW65cuVJt27aV3W5XmzZtlJyc7DwmMTFR48ePlyTdeOONzmGG4pZQzZ599ln17t1bS5cu1Q8//OBsL+i5pqWlqX///goKClLdunX15JNPas2aNZY5INeOnx45ckR16tSRJE2aNMkZb/459aZPPvlEXbp0UVBQkKpXr66+ffvq22+/dTnm66+/VlxcnJo1a6bAwECFhoZq5MiR+vnnn53HuDvH+a/T0qVL1bp1a1WpUkXR0dHas2ePJGn27Nlq0aKFAgMD1a1bN8tr88UXX2jw4MEKCwuT3W5X48aN9eSTT+ry5csux+WX+A8dOqSYmBgFBQWpQYMGeuGFF+Tui2xtNpvmz5+vrKwsZ/z5cwWuXr2qF198Uc2bN5fdblfTpk313HPPKTs726WP/PfQmjVr1K5dO1WpUkWzZ88u1mtRHAXNYSjKBx98oFatWikwMFBRUVHatGmT5Zjjx49r5MiRqlevnvP9884777gck/+eX7x4sf7+97+rYcOGqlq1ql5//XUNHjxYktS9e3fneSvOcMDw4cOVmpqqvXv3OtvS09O1YcMGDR8+3HJ8Tk6OJkyYoKioKIWEhCgoKEhdunTRxo0bXc5PUe+d/J+PgwcP6u6771b16tX1wAMPOPddWyGZOHGi/Pz8tH79epc4/vSnPykgIEBfffWV2+eI8u26rDA0adJEW7du1TfffKO2bdt6dN/jx48738gJCQkKCgrS22+/XWgl4sCBAxo0aJBGjRql2NhYvfPOO4qLi1NUVJTatGlT4H2mTZumd999VytWrNCsWbNUrVo13XzzzQUem5eXp3vuuUfr16/X/fffr8cff1wXLlzQ2rVr9c0336h58+aSpOnTp6t///564IEHlJOTo8WLF2vw4MFavXq1+vbt69Ln5s2btXz5cj366KOqXr26Xn/9dd133306evSoatWqpXvvvVc//PCDPvzwQ7322muqXbu2JDl/sZTEgw8+qE8//VRr165Vy5YtCzwmKytLPXr00IkTJ/T4448rNDRUixYtcvkFV5A6depo1qxZeuSRRzRw4EDde++9klToOb3WhQsXlJmZaWk3f7BJv0xSjY2NVUxMjF555RVdunRJs2bNUufOnfXll186f3muXbtWhw4d0ogRIxQaGqpvv/1Wc+bM0bfffqtt27bJZrMV6xx/8cUXWrVqlcaMGSNJSkpK0j333KOnn35ab775ph599FGdOXNGr776qkaOHKkNGzY477t06VJdunRJjzzyiGrVqqUdO3bojTfe0LFjx7R06VKX55WXl6e77rpLt912m1599VUlJydr4sSJunr1ql544YVCz917772nOXPmaMeOHc6hqE6dOkmSHn74YS1cuFCDBg3SU089pe3btyspKUnff/+9VqxY4dLPvn37NGzYMP35z3/W6NGj1apVq0IfM19Br1lpff7551qyZIn++te/ym63680339Rdd92lHTt2OH+PZGRk6LbbbnMmdHXq1NEnn3yiUaNG6fz583riiSdc+nzxxRcVEBCgcePGKTs7W71799Zf//pXvf7663ruued00003SZLz36J07dpVjRo10qJFi5yvy5IlS1StWjXLe1z6ZW7D22+/rWHDhmn06NG6cOGC5s2bp5iYGO3YsUORkZHFeu9cvXpVMTEx6ty5s/75z3+qatWqBcb397//Xf/5z380atQo7dmzR9WrV9eaNWs0d+5cvfjii4qIiHD/IqB8M65Dn376qeHv72/4+/sb0dHRxtNPP22sWbPGyMnJsRzbpEkTIzY21nl77Nixhs1mM7788ktn288//2zUrFnTkGQcPnzY5b6SjE2bNjnbTp48adjtduOpp55ytm3cuNGQZGzcuNHZNnHiREOScerUKZd47rjjDuOOO+5w3n7nnXcMScbUqVMtsTscDuf/L1265LIvJyfHaNu2rdGjRw+XdklGQECAceDAAWfbV199ZUgy3njjDWfblClTLM+3KLGxsUZQUFCh+7/88ktDkvHkk08628zP9V//+pchyVi5cqWz7fLly0Z4eLjl/MXGxhpNmjRx3j516pQhyZg4cWKx4s1/TYrarn0+Fy5cMGrUqGGMHj3apZ/09HQjJCTEpd38WhiGYXz44YeWn5WizrEkw263u+ybPXu2IckIDQ01zp8/72xPSEiw9FNQDElJSYbNZjPS0tKcbbGxsYYkY+zYsc42h8Nh9O3b1wgICLD8fJoV9LqnpqYakoyHH37YpX3cuHGGJGPDhg3Otvz3UHJycpGPY463qG3MmDHO4w8fPmxIMubPn+9sy3/vXSv/vrt27XK2paWlGYGBgcbAgQOdbaNGjTLq169vZGZmutz//vvvN0JCQpznPf/nq1mzZpbXYunSpZaf56Jc+7ti3LhxRosWLZz72rdvb4wYMcL5HK597levXjWys7Nd+jpz5oxRr149Y+TIkc62ot47+ef72WefLXDfte9BwzCMPXv2GAEBAcbDDz9snDlzxmjYsKHRrl07Izc3t1jPFeXbdTkkceedd2rr1q3q37+/vvrqK7366quKiYlRw4YNtWrVqiLvm5ycrOjoaJflYTVr1nSW4cxat26tLl26OG/XqVNHrVq10qFDh7zyXJYtW6batWu7TNTMd+3SsGvHfM+cOaNz586pS5cuSklJsdyvV69ezsqE9MtfE8HBwV6LuSD5M9svXLhQ6DHJyclq2LCh+vfv72wLDAzU6NGjyyyuCRMmaO3atZatd+/eLsetXbtWZ8+e1bBhw5SZmenc/P391bFjR5cqyLWvxZUrV5SZmanbbrtNkgp8PQrTs2dPl5Jvx44dJUn33Xefqlevbmm/9vW7NoasrCxlZmaqU6dOMgxDX375peWxrp1dn//Xc05OjtatW1fsePP973//kyTFx8e7tD/11FOSpP/+978u7TfeeKNiYmKK3X9gYGCBr9natWs9jvVa0dHRioqKct4OCwvTH/7wB61Zs0Z5eXkyDEPLli1Tv379ZBiGy89BTEyMzp07Z3l9Y2NjvTYfQ/plWOLAgQPauXOn89+ChiMkyd/f3zlnwuFw6PTp07p69aratWvn0c+hJD3yyCPFOq5t27aaNGmS3n77bcXExCgzM1MLFy5UpUrXZTH7N+e6fRXbt2+v5cuXKycnR1999ZVWrFih1157TYMGDVJqaqpat25d4P3S0tIUHR1tab92JcO1wsLCLG033HCDzpw5U7on8P8dPHhQrVq1cvuGW716tV566SWlpqa6lNPN682lso+5IBcvXpQklw86s7S0NDVv3twSc2Hn3ht+//vfq1evXpb2999/3+X2/v37JUk9evQosJ9rV3+cPn1akyZN0uLFiy2TbM+dO1fs2MyvU0hIiCSpcePGBbZf+/odPXpUEyZM0KpVqyyvqzkGPz8/NWvWzKUtf9ioJPNW0tLS5OfnZ3ndQkNDVaNGDaWlpbm033jjjR717+/vX+BrVlq/+93vLG0tW7bUpUuXdOrUKfn5+ens2bOaM2dOocugza+3p8/NnVtuuUXh4eFatGiRatSoodDQ0EJ/JiVp4cKF+te//qW9e/cqNze3RHFVqlRJjRo1Kvbx48eP1+LFi7Vjxw69/PLLhf6uRcVz3SYM+QICAtS+fXu1b99eLVu21IgRI7R06VLnpMPS8vf3L7DdcDNhzJu++OIL9e/fX127dtWbb76p+vXrq3Llypo/f74WLVpkOd4XMX/zzTeSyvbDvyzlLx977733FBoaatl/bUI3ZMgQbdmyRePHj1dkZKSqVasmh8Ohu+66y6NlaIW9Tu5ev7y8PN155506ffq0nnnmGYWHhysoKEjHjx9XXFzcr7YUrqBktSDe/Au8LOWftz/+8Y+KjY0t8BjzvJmyeG7Dhw/XrFmzVL16dQ0dOrTQZeHvv/++4uLiNGDAAI0fP15169aVv7+/kpKSXCZMu2O32z1aen7o0CFngp0/SRfXh+s+YbhWu3btJEknTpwo9JgmTZrowIEDlvaC2n4NzZs31/bt25Wbm1voMqZly5YpMDBQa9ascZmcOX/+/BI/bnF/2RfXe++9J5vNpjvvvLPQY5o0aaLvvvtOhmG4PH5xzr234zXLH8KpW7dukX/dnjlzRuvXr9ekSZM0YcIEZ3v+L9BrlVXMe/bs0Q8//KCFCxfqoYcecrYXVrJ3OBw6dOiQy2TU/NUsJbmSX5MmTeRwOLR//36XyXwZGRk6e/asmjRp4nGfv4aCXqMffvhBVatWdU5GrV69uvLy8kpV4Sjt6z58+HBNmDBBJ06c0HvvvVfocR999JGaNWum5cuXuzym+Y8lb/4cOhwOxcXFKTg4WE888YRefvllDRo0yDmZEhXbdTmHYePGjQX+tZw/tlrULOyYmBht3bpVqampzrbTp0/rgw8+8HqcxXHfffcpMzNTM2bMsOzLf47+/v6y2WzKy8tz7jty5IhWrlxZ4scNCgqSJK9chXDy5Mn69NNPNXTo0ALLvvliYmJ0/Phxl3kmV65c0dy5c90+Rv7M7bK6amJMTIyCg4P18ssvu5R28506dUrS//31b/75mzZtmuU+3jzH1yooBsMwNH369ELvc+3Pl2EYmjFjhipXrqyePXt6/Ph33323JOtznjp1qiQVOKO/PNi6davL2P6PP/6ojz/+WL1795a/v7/8/f113333admyZc6K2bXyfwbcKe3r3rx5c02bNk1JSUnq0KFDoccV9HOwfft2bd261eU4b753pk6dqi1btmjOnDl68cUX1alTJz3yyCNlsqoFv77rssIwduxYXbp0SQMHDlR4eLhycnK0ZcsWLVmyRE2bNtWIESMKve/TTz+t999/X3feeafGjh3rXFYZFham06dPl/lfsmYPPfSQ3n33XcXHx2vHjh3q0qWLsrKytG7dOj366KP6wx/+oL59+2rq1Km66667NHz4cJ08eVIzZ85UixYt9PXXX5focfMnf/3tb3/T/fffr8qVK6tfv37OX3YFuXr1qnPs/8qVK0pLS9OqVav09ddfq3v37m4vf/3nP/9ZM2bM0LBhw/T444+rfv36+uCDD5wX3inq3FepUkWtW7fWkiVL1LJlS9WsWVNt27b1eFltYYKDgzVr1iw9+OCDuvXWW3X//ferTp06Onr0qP773//q9ttv14wZMxQcHKyuXbvq1VdfVW5urho2bKhPP/1Uhw8ftvRZknNcHOHh4WrevLnGjRun48ePKzg4WMuWLSt0jkpgYKCSk5MVGxurjh076pNPPtF///tfPffccyVaShsREaHY2FjNmTNHZ8+e1R133KEdO3Zo4cKFGjBggLp3716q51dW2rZtq5iYGJdllZKcV0GUfkl+N27cqI4dO2r06NFq3bq1Tp8+rZSUFK1bt06nT592+ziRkZHy9/fXK6+8onPnzslut6tHjx6Wa8QU5fHHH3d7zD333KPly5dr4MCB6tu3rw4fPqy33npLrVu3ds4pkrz33vn+++/1/PPPKy4uTv369ZMkLViwQJGRkXr00Uf173//26P+UA75YmlGWfvkk0+MkSNHGuHh4Ua1atWMgIAAo0WLFsbYsWONjIwMl2PNyyoN45clgF26dDHsdrvRqFEjIykpyXj99dcNSUZ6errLffv27Wt5fPNywdIsqzSMX5bI/e1vfzNuvPFGo3LlykZoaKgxaNAg4+DBg85j5s2bZ/zud78z7Ha7ER4ebsyfP7/Q5WPXLr0q6jy8+OKLRsOGDQ0/Pz+3SyzNy92qVq1qNG3a1LjvvvuMjz76yMjLy3N7ngzDMA4dOmT07dvXqFKlilGnTh3jqaeeMpYtW2ZIMrZt2+byeOYlXVu2bDGioqKMgIAAt0ss81+TpUuXFvp8ClomunHjRiMmJsYICQkxAgMDjebNmxtxcXEuy/GOHTtmDBw40KhRo4YREhJiDB482Pjpp58KjKmwc1zQ65S/RHDKlClun8t3331n9OrVy6hWrZpRu3ZtY/To0c7ls9cuMcx/ngcPHjR69+5tVK1a1ahXr54xceLEAl+z4p6n3NxcY9KkSc6f2caNGxsJCQnGlStXXI4r7D3k6ePlM583T5ZVjhkzxnj//fed76NbbrmlwKWPGRkZxpgxY4zGjRs73489e/Y05syZ4zzG3c/X3LlzjWbNmhn+/v5ul1gW9rvC3XN3OBzGyy+/bDRp0sT5fFavXu3Re6eo831tP1evXjXat29vNGrUyDh79qzLcdOnTzckGUuWLCkyfpR/NsP4FWfnVWBPPPGEZs+erYsXLxY66QxlY9q0aXryySd17NgxNWzY0NfhXFfi4uL00UcfufzFCQAFuS7nMJSW+fK5P//8s9577z117tyZZKGMmc/9lStXNHv2bP3ud78jWQAAH7ou5zCUVnR0tLp166abbrpJGRkZmjdvns6fP6/nn3/e16Fd9+69916FhYUpMjJS586d0/vvv6+9e/f6bNIpAOAXJAwFuPvuu/XRRx9pzpw5stlsuvXWWzVv3jx17drV16Fd92JiYvT222/rgw8+UF5enlq3bq3Fixdr6NChvg4NAH7TmMMAAEA5sWnTJk2ZMkW7d+/WiRMntGLFCg0YMKDAY//yl79o9uzZeu211yxffGY2c+ZMTZkyRenp6YqIiNAbb7xR5LLcgjCHAQCAciIrK0sRERGaOXNmkcetWLFC27ZtU4MGDdz2uWTJEsXHx2vixIlKSUlRRESEYmJiLJcyd4eEAQCAcqJPnz566aWXNHDgwEKPOX78uMaOHasPPvig0CsAX2vq1KkaPXq0RowYodatW+utt95S1apV9c4773gUGwkDAABlKDs7W+fPn3fZrv2SQE84HA49+OCDGj9+vNq0aeP2+JycHO3evdvlcuZ+fn7q1auX5aqf7pSbSY8/3nWrr0OosH7af8nXIVRYZ84whac0Ll12fwwKlpPj6wgqtvvz9pVp/7ZHbvNaXxPr3eVyxVDpl+/0SExM9LivV155RZUqVdJf//rXYh2fmZmpvLw81atXz6W9Xr162rt3r0ePXW4SBgAAygubn/e+BiAhIUHx8fEubdd+UWBx7d69W9OnT1dKSsqv/jUFEkMSAACUKbvdruDgYJetJAnDF198oZMnTyosLEyVKlVSpUqVlJaWpqeeeqrQb5atXbu2/P39lZGR4dKekZGh0NBQjx6fhAEAABObn81rm7c8+OCD+vrrr5WamurcGjRooPHjx2vNmjUF3icgIEBRUVFav369s83hcGj9+vWKjo726PEZkgAAwMSbH/SeuHjxog4cOOC8ffjwYaWmpqpmzZoKCwtTrVq1XI6vXLmyQkND1apVK2dbz549NXDgQD322GOSpPj4eMXGxqpdu3bq0KGDpk2bpqysrCK/ubkgJAwAAJj4KmHYtWuXy1fA5899iI2N1YIFC4rVx8GDB5WZmem8PXToUJ06dUoTJkxQenq6IiMjlZycbJkI6U65udIjqyRKjlUSJccqidJhlUTJsUqidMp6lUTlJzt7ra/c1zZ7rS9fosIAAICJL1YhlHckDAAAmPhqSKI8Y5UEAABwiwoDAAAmVBisSBgAADAhYbBiSAIAALhFhQEAABMqDFYkDAAAmJAwWDEkAQAA3KLCAACACRUGKxIGAABMSBisSBgAADDh0tBWzGEAAABuUWEAAMCEIQkrEgYAAExIGKwYkgAAAG5RYQAAwIQKgxUJAwAAJiQMVgxJAAAAt6gwAABgQoXBioQBAAATEgYrhiQAAIBbVBgAADChwmBFwgAAgAkJgxUJAwAAJnz5lBVzGAAAgFtUGAAAMGFIwoqEAQAAExIGK48ThszMTL3zzjvaunWr0tPTJUmhoaHq1KmT4uLiVKdOHa8HCQAAfMujhGHnzp2KiYlR1apV1atXL7Vs2VKSlJGRoddff12TJ0/WmjVr1K5duyL7yc7OVnZ2tmubwyG7H1MqAAC+R4XByqOEYezYsRo8eLDeeustywxSwzD0l7/8RWPHjtXWrVuL7CcpKUmTJk1yaXuieajiW9T3JBwAAMoEf79a2QzDMIp7cJUqVfTll18qPDy8wP179+7VLbfcosuXLxfZT0EVhpODulJhKKGf9l/ydQgV1pkzxf7xRwEuFf1WRxFycnwdQcV2f96+Mu2/0ZyBXuvr2J9WeK0vX/KowhAaGqodO3YUmjDs2LFD9erVc9uP3W6X3W53aTtHsgAAKCf8uQ6DhUcJw7hx4/SnP/1Ju3fvVs+ePZ3JQUZGhtavX6+5c+fqn//8Z5kECgDAr8WfOQwWHiUMY8aMUe3atfXaa6/pzTffVF5eniTJ399fUVFRWrBggYYMGVImgQIAAN/xeFnl0KFDNXToUOXm5iozM1OSVLt2bVWuXNnrwQEA4AsMSViV+MJNlStXVv36rGoAAFx//JlWZ8GVHgEAMKHCYEUOBQAA3KLCAACACRUGKxIGAABMWFZpxZAEAABwiwoDAAAm/hQYLEgYAAAwYUjCiiEJAADgFhUGAABMWCVhRcIAAIAJQxJWDEkAAAC3qDAAAGDCKgkrEgYAAEwYkrAiYQAAwIRJj1bMYQAAAG5RYQAAwIQKgxUVBgAATPz9vLd5YtOmTerXr58aNGggm82mlStXuuxPTExUeHi4goKCdMMNN6hXr17avn17kX0mJibKZrO5bOHh4R6eERIGAADKjaysLEVERGjmzJkF7m/ZsqVmzJihPXv2aPPmzWratKl69+6tU6dOFdlvmzZtdOLECee2efNmj2NjSAIAABNfDUn06dNHffr0KXT/8OHDXW5PnTpV8+bN09dff62ePXsWer9KlSopNDS0VLGRMAAAYOLNZZXZ2dnKzs52abPb7bLb7aXqNycnR3PmzFFISIgiIiKKPHb//v1q0KCBAgMDFR0draSkJIWFhXn0eAxJAABQhpKSkhQSEuKyJSUllbi/1atXq1q1agoMDNRrr72mtWvXqnbt2oUe37FjRy1YsEDJycmaNWuWDh8+rC5duujChQsePS4VBgAATLw5JJGQkKD4+HiXttJUF7p3767U1FRlZmZq7ty5GjJkiLZv3666desWePy1Qxw333yzOnbsqCZNmujf//63Ro0aVezHpcIAAICJN1dJ2O12BQcHu2ylSRiCgoLUokUL3XbbbZo3b54qVaqkefPmFfv+NWrUUMuWLXXgwAGPHpeEAQCACszhcFjmSBTl4sWLOnjwoOrXr+/R45AwAABg4m+zeW3zxMWLF5WamqrU1FRJ0uHDh5WamqqjR48qKytLzz33nLZt26a0tDTt3r1bI0eO1PHjxzV48GBnHz179tSMGTOct8eNG6fPP/9cR44c0ZYtWzRw4ED5+/tr2LBhHsXGHAYAAEx89eVTu3btUvfu3Z238+c+xMbG6q233tLevXu1cOFCZWZmqlatWmrfvr2++OILtWnTxnmfgwcPKjMz03n72LFjGjZsmH7++WfVqVNHnTt31rZt21SnTh2PYiNhAADAxFfXYejWrZsMwyh0//Lly932ceTIEZfbixcvLm1YkhiSAAAAxUCFAQAAE0+/A+K3gIQBAAATvq3SqtwkDMf2XfJ1CBVWVlbh410o2tWrvo4Av1WVys1vX6B4+JEFAMDEnwKDBQkDAAAmfgxJWDCtAwAAuEWFAQAAE4YkrEgYAAAw8dGFHss1hiQAAIBbVBgAADBhSMKKhAEAABM/xiQsSBgAADChwmDFHAYAAOAWFQYAAEwYkbAiYQAAwIQhCSuGJAAAgFtUGAAAMOG7JKxIGAAAMGFIwoohCQAA4BYVBgAATFglYUXCAACAiT9zGCwYkgAAAG5RYQAAwIQhCSsSBgAATFglYUXCAACAiR8D9hacEgAA4BYVBgAATFglYUXCAACACZMerRiSAAAAblFhAADAhFUSViQMAACYMCRhxZAEAABwiwoDAAAmrJKwImEAAMCEIQkrhiQAAIBbVBgAADBhlYQVCQMAACZ+zGGwIGEAAMCECoOV1+cw/Pjjjxo5cmSRx2RnZ+v8+fMuW47h8HYoAADAS7yeMJw+fVoLFy4s8pikpCSFhIS4bO+ePe3tUAAAKBE/m81r2/XC4yGJVatWFbn/0KFDbvtISEhQfHy8S9uXN7fzNBQAAMrE9fRB7y0eJwwDBgyQzWaTYRiFHmNzc6LtdrvsdrtLW4CNFZ4AAJRXHn9K169fX8uXL5fD4ShwS0lJKYs4AQD41TAkYeVxwhAVFaXdu3cXut9d9QEAgPLOz+bnte164fGQxPjx45WVlVXo/hYtWmjjxo2lCgoAAJQvHicMXbp0KXJ/UFCQ7rjjjhIHBACAr11PQwnewoWbAAAwIWGwun4GVwAAQJmhwgAAgAkVBisSBgAATPwowFuQMAAAYEKFwYoUCgAAuEWFAQAAEyoMVlQYAAAw8dWVHjdt2qR+/fqpQYMGstlsWrlypcv+xMREhYeHKygoSDfccIN69eql7du3u+135syZatq0qQIDA9WxY0ft2LHDo7gkEgYAAMqNrKwsRUREaObMmQXub9mypWbMmKE9e/Zo8+bNatq0qXr37q1Tp04V2ueSJUsUHx+viRMnKiUlRREREYqJidHJkyc9is1mlJMvfth6Y7ivQ6iwsrLKxUtYIV254usIKracXF9HUHE5HL6OoGIblL2vTPvfcOw5r/XVo9HLJbqfzWbTihUrNGDAgEKPOX/+vEJCQrRu3Tr17NmzwGM6duyo9u3ba8aMGZIkh8Ohxo0ba+zYsXr22WeLHQ9zGAAAMPHmHIbs7GxlZ2e7tNntdtnt9lL1m5OTozlz5igkJEQRERGFHrN7924lJCQ42/z8/NSrVy9t3brVo8djSAIAgDKUlJSkkJAQly0pKanE/a1evVrVqlVTYGCgXnvtNa1du1a1a9cu8NjMzEzl5eWpXr16Lu316tVTenq6R49LhQEAABNvVhgSEhIUHx/v0laa6kL37t2VmpqqzMxMzZ07V0OGDNH27dtVt27d0oZaJCoMAACYeHOVhN1uV3BwsMtWmoQhKChILVq00G233aZ58+apUqVKmjdvXoHH1q5dW/7+/srIyHBpz8jIUGhoqGfnpMQRAwAAn3M4HJY5EvkCAgIUFRWl9evXuxy/fv16RUdHe/Q4DEkAAGDiJ99cuOnixYs6cOCA8/bhw4eVmpqqmjVrqlatWvrHP/6h/v37q379+srMzNTMmTN1/PhxDR482Hmfnj17auDAgXrsscckSfHx8YqNjVW7du3UoUMHTZs2TVlZWRoxYoRHsZEwAABg4qsrPe7atUvdu3d33s6f+xAbG6u33npLe/fu1cKFC5WZmalatWqpffv2+uKLL9SmTRvnfQ4ePKjMzEzn7aFDh+rUqVOaMGGC0tPTFRkZqeTkZMtESHe4DsN1gOswlBzXYSgdrsNQclyHoXTK+joMOzJe8FpfHepN8FpfvsQcBgAA4BZDEgAAmPDlU1YkDAAAmJAwWDEkAQAA3KLCAACAiadfS/1bQMIAAIAJQxJWpFAAAMAtKgwAAJj46kqP5RkJAwAAJgxJWJEwAABgwqRHK84IAABwiwoDAAAmDElYlZuE4cwZvkCppPgSm5K7etXXEeC3yo/6brlmY0jCgjMCAADcKjcVBgAAygs//p62IGEAAMCEIQkrzggAAHCLCgMAACZch8GKhAEAABMbBXgLzggAAHCLCgMAACYMSViRMAAAYMKQhBUJAwAAJlQYrDgjAADALSoMAACYcOEmKxIGAABMuDS0FWcEAAC4RYUBAAAThiSsSBgAADBhlYQVZwQAALhFhQEAABOb/H0dQrlDwgAAgAlDElacEQAA4BYVBgAATPguCSsSBgAATBiSsCJhAADAhOswWHFGAACAW1QYAAAw4bskrEgYAAAwYUjCijMCAADcosIAAIAJqySsSBgAADDhOgxWnBEAAOAWFQYAAEwYkrAiYQAAwIQhCSvOCAAAcIsKAwAAJgxJWPkkYcjOzlZ2drZLW67hUGVeIABAOcCFm6w8PiOXL1/W5s2b9d1331n2XblyRe+++67bPpKSkhQSEuKy/Tv7tKehAABQJmyG97brhc0wjGI/nR9++EG9e/fW0aNHZbPZ1LlzZy1evFj169eXJGVkZKhBgwbKy8srsp+CKgwbG0dRYSghh8PXEVRcV6/6OoKKzXEd/TJExXLv5X1l+wDGRu/1Zevuvb58yKNP6GeeeUZt27bVyZMntW/fPlWvXl233367jh496tGD2u12BQcHu2wkCwCAcsNweG+7Tng0h2HLli1at26dateurdq1a+s///mPHn30UXXp0kUbN25UUFBQWcUJAMCv5zr6oPcWj/6sv3z5sipV+r8cw2azadasWerXr5/uuOMO/fDDD14PEAAA+J5HFYbw8HDt2rVLN910k0v7jBkzJEn9+/f3XmQAAPgKFQYLjyoMAwcO1IcffljgvhkzZmjYsGHyYA4lAADlk4/mMGzatEn9+vVTgwYNZLPZtHLlSue+3NxcPfPMM/r973+voKAgNWjQQA899JB++umnIvtMTEyUzWZz2cLDwz0+JR4lDAkJCfrf//5X6P4333xTDqbsAwBQIllZWYqIiNDMmTMt+y5duqSUlBQ9//zzSklJ0fLly7Vv375iVffbtGmjEydOOLfNmzd7HBtXegQAwMxHf/z26dNHffr0KXBfSEiI1q5d69I2Y8YMdejQQUePHlVYWFih/VaqVEmhoaGlio21jAAAmHlxSCI7O1vnz5932czXIiqpc+fOyWazqUaNGkUet3//fjVo0EDNmjXTAw884PHlECQSBgAAylRBVzdOSkoqdb9XrlzRM888o2HDhik4OLjQ4zp27KgFCxYoOTlZs2bN0uHDh9WlSxdduHDBo8fz6EqPZel/NVr5OoQKi2kjJceVHkuHKz3CV8r8So+XP/ZaV9l+d1kqCna7XXa7vcj72Ww2rVixQgMGDLDsy83N1X333adjx47ps88+KzJhMDt79qyaNGmiqVOnatSoUcW+H3MYAAAw8+KyyuIkB57Izc3VkCFDlJaWpg0bNniULEhSjRo11LJlSx04cMCj+zEkAQCAmcPhvc2L8pOF/fv3a926dapVq5bHfVy8eFEHDx50fg9UcZEwAABQTly8eFGpqalKTU2VJB0+fFipqak6evSocnNzNWjQIO3atUsffPCB8vLylJ6ervT0dOXk5Dj76Nmzp/OCipI0btw4ff755zpy5Ii2bNmigQMHyt/fX8OGDfMoNoYkAAAw89GVHnft2qXu3f/v2y3j4+MlSbGxsUpMTNSqVaskSZGRkS7327hxo7p16yZJOnjwoDIzM537jh07pmHDhunnn39WnTp11LlzZ23btk116tTxKDYmPV4HmPRYckx6LB0mPcJXynzS4/kl3usreKj3+vIhhiQAAIBbDEkAAGDGl09ZkDAAAGBiGHle68vmtZ58iyEJAADgFhUGAADMmE1uQcIAAIAZcxgsGJIAAABuUWEAAMCMCoMFCQMAAGYkDBYkDAAAmJEwWDCHAQAAuEWFAQAAM5ZVWpAwAABgxpCEBUMSAADALSoMAACYUWGwIGEAAMCMhMGCIQkAAOAWFQYAAMxYJWFBwgAAgBlDEhYMSQAAALfKTYXh0mVfR1Bx+ZH2ARWOn83XEaBIVBgsyk3CAABAucEcBgsSBgAAzByGryModyhmAwAAt6gwAABgxpCEBQkDAABmJAwWDEkAAAC3qDAAAGDGpEcLEgYAAMwYkrBgSAIAALhFhQEAADMqDBYkDAAAmDGHwYIhCQAA4BYVBgAAzBiSsCBhAADAjCEJCxIGAADMqDBYMIcBAAC4RYUBAAAzKgwWJAwAAJgYhvfmMNi81pNvMSQBAADcosIAAIAZQxIWJAwAAJiRMFgwJAEAANyiwgAAgBkXbrIgYQAAwIwhCQuGJAAAgFtUGAAAMKPCYEHCAACAGXMYLEgYAAAwo8JgwRwGAADgFhUGAADMqDBYkDAAAGDGHAYLhiQAAIBbJAwAAJg5HN7bPLBp0yb169dPDRo0kM1m08qVK537cnNz9cwzz+j3v/+9goKC1KBBAz300EP66aef3PY7c+ZMNW3aVIGBgerYsaN27Njh6RkhYQAAwMJHCUNWVpYiIiI0c+ZMy75Lly4pJSVFzz//vFJSUrR8+XLt27dP/fv3L7LPJUuWKD4+XhMnTlRKSooiIiIUExOjkydPehSbzTCMcjFQ85G9la9DqLD8SPuACsfP5usIKrYBl/aVaf+OVQ97rS+//m+X6H42m00rVqzQgAEDCj1m586d6tChg9LS0hQWFlbgMR07dlT79u01Y8YMSZLD4VDjxo01duxYPfvss8WOx+NJj99//722bdum6OhohYeHa+/evZo+fbqys7P1xz/+UT169HDbR3Z2trKzs13acg2HKtv45AMAlANenPRY0Gee3W6X3W4vdd/nzp2TzWZTjRo1Ctyfk5Oj3bt3KyEhwdnm5+enXr16aevWrR49lkef0MnJyYqMjNS4ceN0yy23KDk5WV27dtWBAweUlpam3r17a8OGDW77SUpKUkhIiMu2Iu+0R4EDAFBmvDgkUdBnXlJSUqlDvHLlip555hkNGzZMwcHBBR6TmZmpvLw81atXz6W9Xr16Sk9P9+jxPEoYXnjhBY0fP14///yz5s+fr+HDh2v06NFau3at1q9fr/Hjx2vy5Mlu+0lISNC5c+dctoH+NT0KHACAiqCgz7xr/+IvidzcXA0ZMkSGYWjWrFleirRoHg1JfPvtt3r33XclSUOGDNGDDz6oQYMGOfc/8MADmj9/vtt+CirFMBwBACgvjDzvDUl4a/ghX36ykJaWpg0bNhRaXZCk2rVry9/fXxkZGS7tGRkZCg0N9ehxPf6Uttl+manj5+enwMBAhYSEOPdVr15d586d87RLAADKF4fhvc2L8pOF/fv3a926dapVq1aRxwcEBCgqKkrr16//v6fmcGj9+vWKjo726LE9qjA0bdpU+/fvV/PmzSVJW7dudZmVefToUdWvX9+jAAAAKHe8WGHwxMWLF3XgwAHn7cOHDys1NVU1a9ZU/fr1NWjQIKWkpGj16tXKy8tzzkOoWbOmAgICJEk9e/bUwIED9dhjj0mS4uPjFRsbq3bt2qlDhw6aNm2asrKyNGLECI9i8yhheOSRR5SXl+e83bZtW5f9n3zySbFWSQAAAKtdu3ape/fuztvx8fGSpNjYWCUmJmrVqlWSpMjISJf7bdy4Ud26dZMkHTx4UJmZmc59Q4cO1alTpzRhwgSlp6crMjJSycnJlomQ7nAdhusA12EAKh6uw1A6ZX0dhqvvP+i1vir98T2v9eVLfPkUAABmPhqSKM/42xQAALhFhQEAALM8z74D4reAhAEAABPDy8shrwcMSQAAALeoMAAAYMakRwsSBgAAzBiSsGBIAgAAuEWFAQAAE29++dT1goQBAAAzB8sqzUgYAAAwo8JgwRwGAADgFhUGAABMuHCTFQkDAABmDElYMCQBAADcosIAAIAZFQYLEgYAAEyYw2DFkAQAAHCLCgMAAGZ5XLjJjIQBAAAThiSsGJIAAABuUWEAAMCMVRIWJAwAAJgxJGFBwgAAgAlfb23FHAYAAOAWFQYAAMwYkrAgYQAAwIzrMFgwJAEAANyiwgAAgAkXbrIiYQAAwIxVEhblJmHIyfF1BBVXpXLzKlY8fgzKlQrnD/jt4KMGAAAThiSsSBgAADDhwk1WFBQBAIBbVBgAADBhSMKKhAEAABMHQxIWJAwAAJhQYbBiDgMAAHCLCgMAACaGg++SMCNhAADAhGWVVgxJAAAAt6gwAABgwqRHKxIGAABMGJKwYkgCAAC4RYUBAAAThiSsSBgAADBxkDBYMCQBAADcosIAAIAJkx6tSBgAADBhDoMVCQMAACYkDFbMYQAAAG5RYQAAwIQ5DFYkDAAAmPBtlVYMSQAAALdIGAAAMDHyDK9tnti0aZP69eunBg0ayGazaeXKlS77ly9frt69e6tWrVqy2WxKTU112+eCBQtks9lctsDAQI/ikkgYAACwMByG1zZPZGVlKSIiQjNnzix0f+fOnfXKK6941G9wcLBOnDjh3NLS0jy6v8QcBgAAyo0+ffqoT58+he5/8MEHJUlHjhzxqF+bzabQ0NDShEbCAACAmTe/SyI7O1vZ2dkubXa7XXa73WuP4c7FixfVpEkTORwO3XrrrXr55ZfVpk0bj/pgSAIAABNvzmFISkpSSEiIy5aUlPSrPZdWrVrpnXfe0ccff6z3339fDodDnTp10rFjxzzqhwoDAABlKCEhQfHx8S5tv2Z1ITo6WtHR0c7bnTp10k033aTZs2frxRdfLHY/JAwAAJh489LQv/bwgzuVK1fWLbfcogMHDnh0P4YkAAAw8dWyyl9DXl6e9uzZo/r163t0PyoMAACY+OrLpy5evOjyl//hw4eVmpqqmjVrKiwsTKdPn9bRo0f1008/SZL27dsnSQoNDXWugnjooYfUsGFD5zyJF154QbfddptatGihs2fPasqUKUpLS9PDDz/sUWwkDAAAlBO7du1S9+7dnbfz5z7ExsZqwYIFWrVqlUaMGOHcf//990uSJk6cqMTEREnS0aNH5ef3fwMIZ86c0ejRo5Wenq4bbrhBUVFR2rJli1q3bu1RbDbDMEqdRhmGIZvNVqo+FtlalTaM36xKpH0l5segXKlw/krOr3S/Mn/zBlzaV6b9/9Desw/TorTc+Z3X+vIlr7zd7Xa7vv/+e290BQCAz13PcxhKyqO/Tc3LQvLl5eVp8uTJqlWrliRp6tSpRfZT0EUscuVQZeZgAgBQLnmUMEybNk0RERGqUaOGS7thGPr+++8VFBRUrKGJpKQkTZo0yaXtXtXUfartSTgAAJQJb17p8Xrh0RyGyZMna86cOXr77bfVo0cPZ3vlypX11VdfFXsCRUEVhhUhUVQYSog5DCXHGHzpcP5KjjkMpVPWcxi+vTnca321+Xqv1/ryJY/e7s8++6yWLFmiRx55ROPGjVNubm6JHtRutys4ONhlI1kAAKD88vhTun379tq9e7dOnTqldu3a6Ztvvin1CgkAAMoTh8N72/WiRMXsatWqaeHChVq8eLF69eqlvLw8b8cFAIDPXE8f9N5SqtHv+++/X507d9bu3bvVpEkTb8UEAADKmVJPl2vUqJEaNWrkjVgAACgXWCRhxfx6AABMGJKwImEAAMCEhMGKtYwAAMAtKgwAAJhQYbAiYQAAwISEwYohCQAA4BYVBgAATKgwWJEwAABgQsJgxZAEAABwiwoDAAAmVBisSBgAADAhYbBiSAIAALhFhQEAABPD4NunzEgYAAAwYUjCioQBAAATEgYr5jAAAAC3qDAAAGBChcGKhAEAABMSBiuGJAAAgFtUGAAAMKHCYEXCAACACQmDFUMSAADALSoMAACYUGGwImEAAMDEwZWhLRiSAAAAblFhAADAhCEJKxIGAABMSBisSBgAADAhYbBiDgMAAHCLCgMAACZUGKxshmGweKQI2dnZSkpKUkJCgux2u6/DqXA4fyXHuSs5zl3pcP5QEBIGN86fP6+QkBCdO3dOwcHBvg6nwuH8lRznruQ4d6XD+UNBmMMAAADcImEAAABukTAAAAC3SBjcsNvtmjhxIhN/SojzV3Kcu5Lj3JUO5w8FYdIjAABwiwoDAABwi4QBAAC4RcIAAADcImEAAABukTAAAAC3SBjcmDlzppo2barAwEB17NhRO3bs8HVIFcKmTZvUr18/NWjQQDabTStXrvR1SBVGUlKS2rdvr+rVq6tu3boaMGCA9u3b5+uwKoRZs2bp5ptvVnBwsIKDgxUdHa1PPvnE12FVSJMnT5bNZtMTTzzh61BQTpAwFGHJkiWKj4/XxIkTlZKSooiICMXExOjkyZO+Dq3cy8rKUkREhGbOnOnrUCqczz//XGPGjNG2bdu0du1a5ebmqnfv3srKyvJ1aOVeo0aNNHnyZO3evVu7du1Sjx499Ic//EHffvutr0OrUHbu3KnZs2fr5ptv9nUoKEe4DkMROnbsqPbt22vGjBmSJIfDocaNG2vs2LF69tlnfRxdxWGz2bRixQoNGDDA16FUSKdOnVLdunX1+eefq2vXrr4Op8KpWbOmpkyZolGjRvk6lArh4sWLuvXWW/Xmm2/qpZdeUmRkpKZNm+brsFAOUGEoRE5Ojnbv3q1evXo52/z8/NSrVy9t3brVh5Hht+bcuXOSfvngQ/Hl5eVp8eLFysrKUnR0tK/DqTDGjBmjvn37uvzuAySpkq8DKK8yMzOVl5enevXqubTXq1dPe/fu9VFU+K1xOBx64okndPvtt6tt27a+DqdC2LNnj6Kjo3XlyhVVq1ZNK1asUOvWrX0dVoWwePFipaSkaOfOnb4OBeUQCQNQjo0ZM0bffPONNm/e7OtQKoxWrVopNTVV586d00cffaTY2Fh9/vnnJA1u/Pjjj3r88ce1du1aBQYG+joclEMkDIWoXbu2/P39lZGR4dKekZGh0NBQH0WF35LHHntMq1ev1qZNm9SoUSNfh1NhBAQEqEWLFpKkqKgo7dy5U9OnT9fs2bN9HFn5tnv3bp08eVK33nqrsy0vL0+bNm3SjBkzlJ2dLX9/fx9GCF9jDkMhAgICFBUVpfXr1zvbHA6H1q9fz3goypRhGHrssce0YsUKbdiwQTfeeKOvQ6rQHA6HsrOzfR1GudezZ0/t2bNHqampzq1du3Z64IEHlJqaSrIAKgxFiY+PV2xsrNq1a6cOHTpo2rRpysrK0ogRI3wdWrl38eJFHThwwHn78OHDSk1NVc2aNRUWFubDyMq/MWPGaNGiRfr4449VvXp1paenS5JCQkJUpUoVH0dXviUkJKhPnz4KCwvThQsXtGjRIn322Wdas2aNr0Mr96pXr26ZJxMUFKRatWoxfwaSSBiKNHToUJ06dUoTJkxQenq6IiMjlZycbJkICatdu3ape/fuztvx8fGSpNjYWC1YsMBHUVUMs2bNkiR169bNpX3+/PmKi4v79QOqQE6ePKmHHnpIJ06cUEhIiG6++WatWbNGd955p69DAyo8rsMAAADcYg4DAABwi4QBAAC4RcIAAADcImEAAABukTAAAAC3SBgAAIBbJAwAAMAtEgYAAOAWCQMAAHCLhAEAALhFwgAAANz6f36S6CgeGZTqAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Compare this cell with previously generated instability to your own just above\n",
    "\"\"\"\n",
    "\n",
    "# Specify the file path to the .npy file\n",
    "file_path = './data/previous_mca_hilbert_results.npy'\n",
    "# Load the .npy file\n",
    "data = np.load(file_path)\n",
    "\n",
    "# Inspect the loaded data\n",
    "print(\"Computed average of values: \\n\", np.mean(data,  dtype=np.float64, axis=0))\n",
    "print(\"Computed standard deviation of values: \\n\", np.std(data, dtype=np.float64, axis=0))\n",
    "\n",
    "\n",
    "sigdig=sd.significant_digits(data , reference=np.mean(data, axis=0, dtype=np.float64), basis=10)\n",
    "gc = sns.heatmap(sigdig, cmap='RdYlGn', vmax=14, )\n",
    "gc.set_title(\"Significant Digit Heatmap for Hilbert Matrix\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Significant digits** (Sohier et al., 2021) are generally defined as the number of digits in common between multiple MCA iterations. They are a metric used to quantify numerical uncertainty, and are implemented at this [Github](https://github.com/verificarlo/significantdigits/tree/main/significantdigits)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verrou"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For Verrou, run the Docker image containing Verrou and execute the above scripts again. Results should have a similar variability as Verificarlo and Verrou are both implementing the same type of perturbation, through different implementations.\n",
    "\n",
    "**Quadratic Example**\n",
    "```\n",
    "valgrind --tool=verrou --rounding-mode=random python quadratic_eq.py\n",
    "```\n",
    "**Convergence Example**\n",
    "* As before, we run this example completely in single precision (aka float) in order to emphasize the instability that can occur.\n",
    "```\n",
    "valgrind --tool=verrou --rounding-mode=random --float=yes python convergence_eq.py\n",
    "```\n",
    "**Hibert Example**\n",
    "```\n",
    "valgrind --tool=verrou --rounding-mode=random python hilbert_eq.py\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "* L. N. Trefethen, “The definition of numerical analysis,” Cornell University, Tech. Rep., 1992.\n",
    "* J. Vignes, “Discrete stochastic arithmetic for validating results of numerical software,” Numerical Algorithms, vol. 37, pp. 377–390, 2004.\n",
    "* D. S. Parker, Monte Carlo Arithmetic: Exploiting Randomness in Floating-Point Arithmetic. University of California (Los Angeles). Computer Science Department, 1997.\n",
    "* C. Denis, P. D. O. Castro, and E. Petit, “Verificarlo: Checking Floating Point Accuracy through Monte Carlo Arithmetic,” in 2016 IEEE 23nd Symposium on Computer Arithmetic (ARITH). Los Alamitos, CA, USA: IEEE Computer Society, jul 2016, pp. 55–62. [Online]. Available: https://doi.ieeecomputersociety.org/10.1109/ARITH.2016.31\n",
    "* G. Kiar, P. de Oliveira Castro, P. Rioux, E. Petit, S. T. Brown, A. C. Evans, and T. Glatard, “Comparing Perturbation Models for Evaluating Stability of Neuroimaging Pipelines,” The International Journal of High Performance Computing Applications, vol. 34, no. 5, pp. 491–501, 2020.\n",
    "* F. Févotte and B. Lathuilière, “Verrou: Assessing Floating-Point Accuracy Without Recompiling,” 2016.\n",
    "* N. Nethercote and J. Seward, “Valgrind: A Framework for Heavyweight Dynamic Binary Instrumentation,” ACM Sigplan notices, vol. 42, no. 6, pp. 89–100, 2007.\n",
    "* F. Jezequel, J.-L. Lamotte, and O. Chubach, “Parallelization of discrete stochastic arithmetic on multicore architectures,” in 2013 10th International Conference on Information Technology: New Generations. IEEE, 2013, pp. 160–166.\n",
    "* D. Sohier, P. D. O. Castro, F. Févotte, B. Lathuilière, E. Petit, and O. Jamond, “Confidence Intervals for Stochastic Arithmetic,” ACM Transactions on Mathematical Software (TOMS), vol. 47, no. 2, pp. 1–33, 2021."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tutorial Expansion Ideas\n",
    "\n",
    "* Significant digit viewer as done here (https://github.com/yohanchatelain/significantdigit-viewer)\n",
    "* Neuroscience example\n",
    "* Lead to applications and work done?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tutorial2",
   "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.11.0rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}