{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adaptive PDE discretizations on Cartesian grids\n",
    "## Volume : Algorithmic tools\n",
    "## Part : Generalized acuteness\n",
    "## Chapter : Riemannian norms and the Voronoi vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook illustrates a construction of acute stencils applicable to two and three dimensional Riemannian metrics. It allows to compute Riemannian distance maps, by solving the corresponding eikonal equations, efficiently in a single pass over the discretization set, using the Fast-Marching method, a Dijkstra-like algorithm. The present notebook is only illustrative, as this algorithms are more efficiently implemented in C++, due to its serial nature.\n",
    "\n",
    "The stencil construction presented in this notebook is based on the following mathematical tools:\n",
    "* Voronoi vectors. Those vectors $e\\in Z^d \\setminus \\{0\\}$ whose Voronoi cell intersects the Voronoi cell of the origin. \n",
    "* The notion of acute angle associated with a scalar product.\n",
    "\n",
    "Here the scalar product and norm on $R^d$ are associated with a $d\\times d$ symmetric positive definite matrix $M$, and defined as follows:\n",
    "$$\n",
    "    <u,v>_M := <u,Mv> \\qquad \\|u\\|_M := \\sqrt{<u,M u>}\n",
    "$$\n",
    "\n",
    "**Comparison with the Stern-Brocot based stencils.**\n",
    "Another stencil construction is presented in the notebook [SternBrocot](SternBrocot.ipynb). This one distinguishes itself on the following points:\n",
    "* It only applies to *Riemannian* metric. In constrast, the other construction extends to *Finslerian* metrics.\n",
    "* It applies in dimension $3$. In constrast, the other construction is limited to dimension two.\n",
    "\n",
    "**Academic publication.** The contents of this notebooks are related with the following publication:\n",
    "* Mirebeau, J.-M. (2014). Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction. SIAM Journal on Numerical Analysis, 52(4), 1573–1599."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[**Summary**](Summary.ipynb) of volume Algorithmic tools, this series of notebooks.\n",
    "\n",
    "[**Main summary**](../Summary.ipynb) of the Adaptive Grid Discretizations \n",
    "\tbook of notebooks, including the other volumes.\n",
    "\n",
    "# Table of contents\n",
    "  * [1. Stencil and acuteness](#1.-Stencil-and-acuteness)\n",
    "    * [1.1 Two dimensional stencils](#1.1-Two-dimensional-stencils)\n",
    "    * [1.2 Three dimensional stencils](#1.2-Three-dimensional-stencils)\n",
    "  * [2. Two dimensional construction](#2.-Two-dimensional-construction)\n",
    "  * [3. Three dimensional construction](#3.-Three-dimensional-construction)\n",
    "\n",
    "\n",
    "\n",
    "**Acknowledgement.** Some of the experiments presented in these notebooks are part of \n",
    "ongoing research with Ludovic Métivier and Da Chen.\n",
    "\n",
    "Copyright Jean-Marie Mirebeau, Centre Borelli, ENS Paris-Saclay, CNRS, University Paris-Saclay"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Importing the required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.143349Z",
     "iopub.status.busy": "2024-04-30T08:46:15.143024Z",
     "iopub.status.idle": "2024-04-30T08:46:15.151497Z",
     "shell.execute_reply": "2024-04-30T08:46:15.151007Z"
    }
   },
   "outputs": [],
   "source": [
    "import sys; sys.path.insert(0,\"..\") # Allow import of agd from parent directory (useless if conda package installed)\n",
    "#from Miscellaneous import TocTools; print(TocTools.displayTOC('VoronoiVectors','Algo'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.154090Z",
     "iopub.status.busy": "2024-04-30T08:46:15.153890Z",
     "iopub.status.idle": "2024-04-30T08:46:15.259434Z",
     "shell.execute_reply": "2024-04-30T08:46:15.259042Z"
    }
   },
   "outputs": [],
   "source": [
    "from agd.Metrics import Riemann\n",
    "from agd import Selling\n",
    "from agd import LinearParallel as lp\n",
    "from agd import FiniteDifferences as fd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.261221Z",
     "iopub.status.busy": "2024-04-30T08:46:15.261078Z",
     "iopub.status.idle": "2024-04-30T08:46:15.413529Z",
     "shell.execute_reply": "2024-04-30T08:46:15.413212Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Stencil and acuteness\n",
    "\n",
    "In the context of semi-Lagrangian discretizations of eikonal equations, a stencil at a discretization point $x$ is defined as a polygonal surface enclosing $x$, whose vertices are also discretization points.\n",
    "\n",
    "We specialize to cartesian grids discretizations, and assume w.l.o.g. that the discretization point considered is the origin. Up to rescaling, the polygonal surface vertices should therefore lie on $Z^d$.\n",
    "With these conventions the stencil is said *acute*, w.r.t. the scalar product associated to a matrix $M$, iff for any vertices $u,v$ of a common facet of the polygonal surface one has \n",
    "$$\n",
    "    <u,M v> \\, \\geq 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.415234Z",
     "iopub.status.busy": "2024-04-30T08:46:15.415126Z",
     "iopub.status.idle": "2024-04-30T08:46:15.416979Z",
     "shell.execute_reply": "2024-04-30T08:46:15.416758Z"
    }
   },
   "outputs": [],
   "source": [
    "direction2 = [np.cos(np.pi/6),np.sin(np.pi/6)] # Arbitrary non-zero vector\n",
    "direction3 = [1.5,4.2,10.7] # Arbitrary non-zero vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.418394Z",
     "iopub.status.busy": "2024-04-30T08:46:15.418300Z",
     "iopub.status.idle": "2024-04-30T08:46:15.420321Z",
     "shell.execute_reply": "2024-04-30T08:46:15.420085Z"
    }
   },
   "outputs": [],
   "source": [
    "aX0 = np.linspace(-1,1); aX1=aX0\n",
    "X = np.array(np.meshgrid(aX0,aX1,indexing='ij')) # Coordinate system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Two dimensional stencils\n",
    "\n",
    "We represent two dimensional stencils as a list of vectors. The $4$-element diamond stencil, and $8$-element square stencil are most classical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.421677Z",
     "iopub.status.busy": "2024-04-30T08:46:15.421586Z",
     "iopub.status.idle": "2024-04-30T08:46:15.423604Z",
     "shell.execute_reply": "2024-04-30T08:46:15.423365Z"
    }
   },
   "outputs": [],
   "source": [
    "e0,e1=np.eye(2)\n",
    "diamond2 = [e0,e1,-e0,-e1,e0]\n",
    "box2 = [e0,e0+e1,e1,e1-e0,-e0,-e0-e1,-e1,-e1+e0,e0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.424923Z",
     "iopub.status.busy": "2024-04-30T08:46:15.424846Z",
     "iopub.status.idle": "2024-04-30T08:46:15.426900Z",
     "shell.execute_reply": "2024-04-30T08:46:15.426633Z"
    }
   },
   "outputs": [],
   "source": [
    "def PlotStencil2(stencil):\n",
    "    plt.plot(*np.array(stencil).T)\n",
    "    plt.scatter(*np.array(stencil).T)\n",
    "    plt.scatter(0,0,color='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.428278Z",
     "iopub.status.busy": "2024-04-30T08:46:15.428198Z",
     "iopub.status.idle": "2024-04-30T08:46:15.568092Z",
     "shell.execute_reply": "2024-04-30T08:46:15.567789Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/QAAAF0CAYAAACNCnKTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABP6ElEQVR4nO3deXgUZb728bsJWQBJs0SSIEjAg+wqywCJrIJhGZCRo4BABF9gxBkHI3JE3FgcjXBcRwdcBmFUBlADoyOLBFk9hE0SFARkDluQBAShO4okkNT7ByetTXeWDl3p7fu5rrou+8lT1U8VZf36rqquthiGYQgAAAAAAASUar4eAAAAAAAA8ByBHgAAAACAAESgBwAAAAAgABHoAQAAAAAIQAR6AAAAAAACEIEeAAAAAIAARKAHAAAAACAAEegBAAAAAAhABHoAAAAAAAIQgR6opIULF8pisTimqKgoxcXFqXfv3kpLS9OpU6dc5pkxY4YsFosPRlu1evXqpV69el31ck6cOKEZM2YoOzv7qpflLWPHjlVCQoJTW0JCgsaOHeuT8QAAfG/btm268847df311ysyMlKxsbFKTEzUI4884uuh+YW5c+dq4cKFvh6Gw4YNG2SxWLRhwwZHW6h8RkPwIdADV2nBggXKzMxURkaG/vrXv+qWW27R7Nmz1apVK61du9ap7/jx45WZmemjkQaeEydOaObMmX4V6J966iktX77c18MAAPiJFStWKCkpSXa7XXPmzNGaNWv06quv6tZbb9XSpUt9PTy/4G+BvkOHDsrMzFSHDh18PRTgqlX39QCAQNe2bVt16tTJ8fo///M/9fDDD6tbt24aOnSoDh48qNjYWElSo0aN1KhRI18NFV5www03+HoIAAA/MmfOHDVt2lSfffaZqlf/5aP1iBEjNGfOHB+OrHznz59XzZo1fT2MKhcdHa2uXbv6ehiAV3CFHjDB9ddfrxdffFH5+fl68803He3ubudaunSpkpOTFR8frxo1aqhVq1Z67LHH9NNPPzn1Gzt2rK655hrt379f/fr1U61atRQfH6/nn39ekrR161Z169ZNtWrV0o033qi///3vLuPas2ePhgwZorp16yoqKkq33HKLS7+S29AWL16sJ554Qg0bNlR0dLT69u2rAwcOOPU1DENz5sxRkyZNFBUVpQ4dOmjVqlUV3k4ffvihunTpIqvVqpo1a6pZs2b6f//v/znG8Zvf/EaSdN999zm+2jBjxgzH/Dt37tQdd9yhevXqKSoqSu3bt9cHH3zg9B4lX41Yv369HnjgAcXExKh+/foaOnSoTpw44TKmf/zjH0pMTNQ111yja665Rrfccovmz5/v9O9w5S33AIDQdebMGcXExDiF+RLVqjl/1L548aIeffRRxcXFqWbNmurWrZu2b9/u8tWt0m7/LqlpR44ccbR5+jni66+/VnJysmrXrq0+ffpIkgoLC/XnP/9ZLVu2VGRkpK699lrdd999+v7778td/0OHDmnEiBFq2LCh4+sGffr0cdxdl5CQoL1792rjxo2OWv7rOmq32zVlyhQ1bdpUERERuu6665SamuoyfovFogcffFDvvfeeWrVqpZo1a+rmm2/Wp59+6jKm/fv365577lFsbKwiIyN1/fXX695771VBQYEk97fcA4GKK/SASQYOHKiwsDBt2rSpzH4HDx7UwIEDlZqaqlq1amn//v2aPXu2tm/frnXr1jn1vXjxooYOHaqJEyfqv/7rv/SPf/xD06ZNk91uV3p6uqZOnapGjRrptdde09ixY9W2bVt17NhRknTgwAElJSWpQYMG+stf/qL69evr/fff19ixY3Xy5Ek9+uijTu/1+OOP69Zbb9Xf/vY32e12TZ06VYMHD9a+ffsUFhYmSZo5c6ZmzpypcePG6a677lJOTo4mTJigoqIitWjRosz1zszM1PDhwzV8+HDNmDFDUVFROnr0qGOdO3TooAULFui+++7Tk08+qd/+9reS5LjDYf369erfv7+6dOmiN954Q1arVUuWLNHw4cN1/vx5l++0jx8/Xr/97W/1j3/8Qzk5Ofqv//ovjR492mkbP/3003rmmWc0dOhQPfLII7JardqzZ4+OHj1a5roAAEJXYmKi/va3v2nSpEkaNWqUOnTooPDwcLd9J0yYoHfffVdTpkzR7bffrj179mjo0KHKz8+v9Pt78jmisLBQd9xxh+6//3499thjunTpkoqLizVkyBBt3rxZjz76qJKSknT06FFNnz5dvXr10s6dO1WjRo1S33/gwIEqKirSnDlzdP311+v06dPasmWLzp07J0lavny57rrrLlmtVs2dO1eSFBkZKenyHQI9e/bU8ePH9fjjj+umm27S3r179fTTT+vrr7/W2rVrnU5srFixQjt27NCsWbN0zTXXaM6cObrzzjt14MABNWvWTJK0e/dudevWTTExMZo1a5aaN2+u3NxcffLJJyosLHS8NxA0DACVsmDBAkOSsWPHjlL7xMbGGq1atXK8nj59ulHW/3bFxcXGxYsXjY0bNxqSjN27dzv+NmbMGEOSkZ6e7mi7ePGice211xqSjF27djnaz5w5Y4SFhRmTJ092tI0YMcKIjIw0jh075vSeAwYMMGrWrGmcO3fOMAzDWL9+vSHJGDhwoFO/Dz74wJBkZGZmGoZhGGfPnjWioqKMO++806nf//zP/xiSjJ49e5a6noZhGC+88IIhyfG+7uzYscOQZCxYsMDlby1btjTat29vXLx40al90KBBRnx8vFFUVGQYxi//Tn/4wx+c+s2ZM8eQZOTm5hqGYRiHDh0ywsLCjFGjRpU57jFjxhhNmjRxamvSpIkxZsyYMucDAASn06dPG926dTMkGZKM8PBwIykpyUhLSzPy8/Md/fbt22dIMh5++GGn+RctWmRIcqojpX1eKKlphw8fdjuWinyOeOedd5zmWbx4scvnC8P4pQbPnTu3zHWXZLzyyiul9jEMw2jTpo3bzwVpaWlGtWrVXD5LffTRR4YkY+XKlY42SUZsbKxht9sdbXl5eUa1atWMtLQ0R9ttt91m1KlTxzh16lSp4yn5rLN+/XpHW3mf0QB/xS33gIkMwyi3z6FDhzRy5EjFxcUpLCxM4eHh6tmzpyRp3759Tn0tFosGDhzoeF29enX9x3/8h+Lj49W+fXtHe7169dSgQQOnK8vr1q1Tnz591LhxY6dljh07VufPn3d5WN8dd9zh9Pqmm26SJMcyMzMzdeHCBY0aNcqpX1JSkpo0aVLuepfcTj9s2DB98MEH+u6778qdp8S///1v7d+/3/Hely5dckwDBw5Ubm6uy9cDylufjIwMFRUV6Y9//GOFxwEAQP369bV582bt2LFDzz//vIYMGaJvv/1W06ZNU7t27XT69GlJl+8sk+RSN4cNG+b2dv2K8uRzhHT5WT+/9umnn6pOnToaPHiwUz295ZZbFBcXV+Zt6fXq1dMNN9yg//7v/9ZLL72krKwsFRcXV3jsn376qdq2batbbrnF6b379evn9pb43r17q3bt2o7XsbGxTp93zp8/r40bN2rYsGG69tprKzwOIJAR6AGT/PTTTzpz5owaNmxYap8ff/xR3bt317Zt2/TnP/9ZGzZs0I4dO7Rs2TJJ0s8//+zUv2bNmoqKinJqi4iIUL169VyWHRERoQsXLjhenzlzRvHx8S79SsZ35swZp/b69es7vS65Ra1kTCX94+LiXJbpru1KPXr00D//+U9dunRJ9957rxo1aqS2bdtq8eLF5c578uRJSdKUKVMUHh7uNP3hD3+QJMcHqIquT8n3BHloIQCgMjp16qSpU6fqww8/1IkTJ/Twww/ryJEjjgfjlVY3q1ev7lKjKqoynyOio6Od2k6ePKlz584pIiLCpabm5eW51NNfs1gs+vzzz9WvXz/NmTNHHTp00LXXXqtJkyZV6GsEJ0+e1FdffeXyvrVr15ZhGOXWculyPS9Zz7Nnz6qoqIhajpDCd+gBk6xYsUJFRUVl/h77unXrdOLECW3YsMFxNl2S43tn3lS/fn3l5ua6tJc8GC4mJsbj5UlSXl6ey9/y8vIq9OC4IUOGaMiQISooKNDWrVuVlpamkSNHKiEhQYmJiaXOVzLWadOmaejQoW77lPcd/iuVnMk/fvy4y10MAAB4Ijw8XNOnT9fLL7+sPXv2SHKum9ddd52j76VLl1xOqpecvC8oKHD6zveVAdfTzxHuHrRX8rDY1atXu53n11fE3WnSpInj4bHffvutPvjgA82YMUOFhYV64403ypw3JiZGNWrU0DvvvFPq3z1Rr149hYWF6fjx4x7NBwQyAj1ggmPHjmnKlCmyWq26//77S+1XUlivfEDLr5+M7y19+vTR8uXLdeLECae7Bt59913VrFnT459v6dq1q6KiorRo0SKn2/e2bNmio0ePevQk+MjISPXs2VN16tTRZ599pqysLCUmJrpcRS/RokULNW/eXLt379Zzzz3n0bhLk5ycrLCwMM2bN6/MkwkAAPxabm6u2zvgSm53L6m5JSf4Fy1a5HhgrSR98MEHunTpktO8JTX0q6++cnxFTZL+9a9/OfXzxueIQYMGacmSJSoqKlKXLl0qPJ87N954o5588kmlp6dr165djvZfX0W/8r2fe+451a9fX02bNr2q95akGjVqqGfPnvrwww/17LPPenxCAAhEBHrgKu3Zs8fxna9Tp05p8+bNWrBggcLCwrR8+fIyv8OVlJSkunXrauLEiZo+fbrCw8O1aNEi7d692+vjnD59uj799FP17t1bTz/9tOrVq6dFixZpxYoVmjNnjqxWq0fLq1u3rqZMmaI///nPGj9+vO6++27l5ORoxowZFbrl/umnn9bx48fVp08fNWrUSOfOndOrr77q9N2/G264QTVq1NCiRYvUqlUrXXPNNWrYsKEaNmyoN998UwMGDFC/fv00duxYXXfddfrhhx+0b98+7dq1Sx9++KFH65OQkKDHH39czzzzjH7++Wfdc889slqt+uabb3T69GnNnDnTo+UBAEJDv3791KhRIw0ePFgtW7ZUcXGxsrOz9eKLL+qaa67RQw89JElq1aqVRo8erVdeeUXh4eHq27ev9uzZoxdeeMHlNviBAweqXr16GjdunGbNmqXq1atr4cKFysnJcernjc8RI0aM0KJFizRw4EA99NBD6ty5s8LDw3X8+HGtX79eQ4YM0Z133ul23q+++koPPvig7r77bjVv3lwRERFat26dvvrqKz322GOOfu3atdOSJUu0dOlSNWvWTFFRUWrXrp1SU1OVnp6uHj166OGHH9ZNN92k4uJiHTt2TGvWrNEjjzzi8UmGl156Sd26dVOXLl302GOP6T/+4z908uRJffLJJ3rzzTfLveMACDQEeuAq3XfffZIuf2e9Tp06atWqlaZOnarx48eX+0CW+vXra8WKFXrkkUc0evRo1apVS0OGDNHSpUvVoUMHr46zRYsW2rJlix5//HH98Y9/1M8//6xWrVppwYIFLj/xVlGzZs1SrVq1NHfuXL333ntq2bKl3njjDb3wwgvlztulSxft3LlTU6dO1ffff686deqoU6dOWrdundq0aSPp8nf93nnnHc2cOVPJycm6ePGipk+frhkzZqh3797avn27nn32WaWmpurs2bOqX7++WrdurWHDhlV6fZo3b67XXntNo0aNUvXq1dW8eXNNmjSpUssDAAS/J598Uh9//LFefvll5ebmqqCgQPHx8erbt6+mTZumVq1aOfrOnz9fsbGxWrhwof7yl7/olltuUXp6ukaMGOG0zOjoaK1evVqpqakaPXq06tSpo/Hjx2vAgAEaP368o583PkeEhYXpk08+0auvvqr33ntPaWlpql69uho1aqSePXuqXbt2pc4bFxenG264QXPnzlVOTo4sFouaNWumF198UX/6058c/WbOnKnc3FxNmDBB+fn5atKkiY4cOaJatWpp8+bNev755/XWW2/p8OHDqlGjhq6//nr17dvXo7v9Stx8883avn27pk+frmnTpik/P19xcXG67bbbFBER4fHyAH9nMSryGG4AAAAApkhISFCvXr20cOFCXw8FQIDhKfcAAAAAAAQgAj0AAAAAAAGIW+4BAAAAAAhAXKEHAAAAACAAEegBAAAAAAhABHoAAAAAAAJQSP4OfXFxsU6cOKHatWvLYrH4ejgAAMgwDOXn56thw4aqVo3z7d5AvQcA+BMzan1IBvoTJ06ocePGvh4GAAAucnJy1KhRI18PIyhQ7wEA/sibtT4kA33t2rUlXd6Q0dHRPh4NAACS3W5X48aNHTUKV496DwDwJ2bU+pAM9CW33UVHR1PgAQB+hVvDvYd6DwDwR96s9XxJDwAAAACAAESgBwAAAAAgABHoAQAAAAAIQAR6AAAAAAACEIEeAAAAAIAARKAHAAAAACAAheTP1gHwjqJiQ9sP/6BT+RfUoHaUOjetp7Bq/OQWAPgzjt3wBPsLPMH+UvVMvUK/adMmDR48WA0bNpTFYtE///nPcufZuHGjOnbsqKioKDVr1kxvvPGGS5/09HS1bt1akZGRat26tZYvX27C6AGUZfWeXHWbvU73vL1VDy3J1j1vb1W32eu0ek+ur4cGoApR6wMLx254gv0FnmB/8Q1TA/1PP/2km2++Wa+//nqF+h8+fFgDBw5U9+7dlZWVpccff1yTJk1Senq6o09mZqaGDx+ulJQU7d69WykpKRo2bJi2bdtm1moAuMLqPbl64P1dyrVdcGrPs13QA+/v4sANhBBqfeDg2A1PsL/AE+wvvmMxDMOokjeyWLR8+XL97ne/K7XP1KlT9cknn2jfvn2OtokTJ2r37t3KzMyUJA0fPlx2u12rVq1y9Onfv7/q1q2rxYsXV2gsdrtdVqtVNptN0dHRlVshIEQVFRvqNnudywG7hEVSnDVKX0y9jVusAA8EQ23yp1ovBcc29Zbyjt2SFBsdqbWTe3LshoqKDfV9aaNO2gtK7cP+ghLl7S98NvyFGXXJr75Dn5mZqeTkZKe2fv36af78+bp48aLCw8OVmZmphx9+2KXPK6+8UupyCwoKVFDwyw5mt9u9Om4glGw//EOZHwgNSbm2C9p++Acl3lC/6gYGICCYVesl6n1Zyjt2S9JJe4HazVhTRSNCoGN/QUXx2dBcfvWU+7y8PMXGxjq1xcbG6tKlSzp9+nSZffLy8kpdblpamqxWq2Nq3Lix9wcPhIhT+WV/IPS0H4DQYlatl6j3ZeGYDMDXOA6Zw6+u0EuXb9f7tZJvBPy63V2fK9t+bdq0aZo8ebLjtd1up8gDldSgdpRX+wEIPWbUeol6X5aKHpMXjP2NujSrZ/Jo4O+2HfpB9y3cUW4/9hdIFd9f+GxoDr8K9HFxcS5n30+dOqXq1aurfv36Zfa58kz+r0VGRioyMtL7AwZC0HV1aijMYlFRGY/fiLde/pkSALiSWbVeot6XpXPTeoq3RinPdkHujt4l33HtceO1If8dV0g9bryW/QUVVtH9hc+G5vCrW+4TExOVkZHh1LZmzRp16tRJ4eHhZfZJSkqqsnECoer42fMaNX9rmWFeknpS4AGUglrvG2HVLJo+uLXbv5UcracPbs2xG5Kc95cr9wj2F1yJ44tvmRrof/zxR2VnZys7O1vS5Z+qyc7O1rFjxyRdvjXu3nvvdfSfOHGijh49qsmTJ2vfvn165513NH/+fE2ZMsXR56GHHtKaNWs0e/Zs7d+/X7Nnz9batWuVmppq5qoAIe/42fO65+2tyvnhZyXUr6m0O9sq3up861StiDBJ0pIdOXov84gPRgmgqlHrA0f/tvGaN7qDYqOd72KIs0Zp3ugO6t823kcjgz8q2V/irqj17C9wh+OLDxkmWr9+vaHLDzZ0msaMGWMYhmGMGTPG6Nmzp9M8GzZsMNq3b29EREQYCQkJxrx581yW++GHHxotWrQwwsPDjZYtWxrp6ekejctmsxmSDJvNVtlVA0JKzg8/Gd1mf240mfqp0XPOOuPEufOGYRjGpaJiY8u/Txv/zDpubPn3aePipSLjuRXfGE2mfmo0mfqp8e6Ww74dOBBAArU2+WutN4zA3aZms/9c6DhOr9t30rhUVOzrIcGPXVnr2V9QFo4vZTOjLlXZ79D7E36XFqi4K6/ML/59V8Vba5Ta3zAMPb9qv97cdEiS9MyQNkpJTKii0QKBi9rkfWxT984XXlLrpz+TJH0zq59qRvjVI5UABDCOL2Uzoy751XfoAfgXT8O8dPnJ1I8NaKn7ezSTJD318V5uvwcAAABMQKAH4FZlwnwJQj0AAABgPgI9ABdXE+ZLEOoBAAAAcxHoATjxRpgvQagHAAAAzEOgB+DgzTBfglAPAAAAmINAD0CSOWG+BKEeAAAA8D4CPQBTw3wJQj0AAADgXQR6IMRVRZgvQagHAAAAvIdAD4SwqgzzJQj1AAAAgHcQ6IEQ5YswX4JQDwAAAFw9Aj0QgnwZ5ksQ6gEAAICrQ6AHQow/hPkShHoAAACg8gj0QAjxpzBfglAPAAAAVA6BHggR/hjmSxDqAQAAAM8R6IEQ4M9hvgShHgAAAPAMgR4IcoEQ5ksQ6gEAAICKI9ADQSyQwnwJQj0AAABQMQR6IEgFYpgvQagHAAAAykegB4JQIIf5EoR6AAAAoGwEeiDIBEOYL0GoBwAAAEpHoAeCSDCF+RKEegAAAMA9Aj0QJIIxzJcg1AMAAACuCPRAEAjmMF+CUA8AAAA4I9ADAS4UwnwJQj0AAADwC9MD/dy5c9W0aVNFRUWpY8eO2rx5c6l9x44dK4vF4jK1adPG0WfhwoVu+1y4cMHsVQH8TiiF+RKEesA/Ue8BAKh6pgb6pUuXKjU1VU888YSysrLUvXt3DRgwQMeOHXPb/9VXX1Vubq5jysnJUb169XT33Xc79YuOjnbql5ubq6ioKDNXBfA7oRjmSxDqAf9CvQcAwDdMDfQvvfSSxo0bp/Hjx6tVq1Z65ZVX1LhxY82bN89tf6vVqri4OMe0c+dOnT17Vvfdd59TP4vF4tQvLi7OzNUA/E4oh/kShHrAf1DvAQDwDdMCfWFhob788kslJyc7tScnJ2vLli0VWsb8+fPVt29fNWnSxKn9xx9/VJMmTdSoUSMNGjRIWVlZZS6noKBAdrvdaQICFWH+F4R6wPeo9wAA+I5pgf706dMqKipSbGysU3tsbKzy8vLKnT83N1erVq3S+PHjndpbtmyphQsX6pNPPtHixYsVFRWlW2+9VQcPHix1WWlpabJarY6pcePGlVspwMcI864I9YBvUe8BAPAd0x+KZ7FYnF4bhuHS5s7ChQtVp04d/e53v3Nq79q1q0aPHq2bb75Z3bt31wcffKAbb7xRr732WqnLmjZtmmw2m2PKycmp1LoAvkSYLx2hHvA96j0AAFWvulkLjomJUVhYmMvZ+VOnTrmcxb+SYRh65513lJKSooiIiDL7VqtWTb/5zW/KPGMfGRmpyMjIig8e8DOE+fKVhHpJenPTIT318V5JUkpigg9HBQQ/6j0AAL5j2hX6iIgIdezYURkZGU7tGRkZSkpKKnPejRs36t///rfGjRtX7vsYhqHs7GzFx8df1XgBf0WYrziu1ANVj3oPAIDvmHaFXpImT56slJQUderUSYmJiXrrrbd07NgxTZw4UdLlW+O+++47vfvuu07zzZ8/X126dFHbtm1dljlz5kx17dpVzZs3l91u11/+8hdlZ2frr3/9q5mrAvgEYd5zXKkHqh71HgAA3zA10A8fPlxnzpzRrFmzlJubq7Zt22rlypWOp9jm5ua6/EatzWZTenq6Xn31VbfLPHfunH7/+98rLy9PVqtV7du316ZNm9S5c2czVwWocoT5yiPUA1WLeg8AgG9YDMMwfD2Iqma322W1WmWz2RQdHe3r4QAuCPPeYRiGnl+1X29uOiRJemZIG0I9/Ba1yfvYpu6dL7yk1k9/Jkn6ZlY/1Yww9foOgBDC8aVsZtQl059yD8AzhHnv4Tv1AAAACGYEesCPEOa9j1APAACAYEWgB/wEYd48hHoAAAAEIwI94AcI8+Yj1AMAACDYEOgBHyPMVx1CPQAAAIIJgR7wIcJ81SPUAwAAIFgQ6AEfIcz7DqEeAAAAwYBAD/gAYd73CPUAAAAIdAR6oIoR5v0HoR4AAACBjEAPVCHCvP8h1AMAACBQEeiBKkKY91+EegAAAAQiAj1QBQjz/o9QDwAAgEBDoAdMRpgPHIR6AAAABBICPWAiwnzgIdQDAAAgUBDoAZMQ5gMXoR4AAACBgEAPmIAwH/gI9QAAAPB3BHrAywjzwYNQDwAAAH9GoAe8iDAffAj1AAAA8FcEesBLCPPBi1APAAAAf0SgB7yAMB/8CPUAAADwNwR64CoR5kMHoR4AAAD+hEAPXAXCfOgh1AMAAMBfEOiBSiLMhy5CPQAAAPwBgR6oBMI8CPUAAADwNdMD/dy5c9W0aVNFRUWpY8eO2rx5c6l9N2zYIIvF4jLt37/fqV96erpat26tyMhItW7dWsuXLzd7NQAHwjxKEOqBX1DvAQCoeqYG+qVLlyo1NVVPPPGEsrKy1L17dw0YMEDHjh0rc74DBw4oNzfXMTVv3tzxt8zMTA0fPlwpKSnavXu3UlJSNGzYMG3bts3MVQEkEebhilAPUO8BAPAVi2EYhlkL79Klizp06KB58+Y52lq1aqXf/e53SktLc+m/YcMG9e7dW2fPnlWdOnXcLnP48OGy2+1atWqVo61///6qW7euFi9eXKFx2e12Wa1W2Ww2RUdHe7ZSCFmEeZTFMAw9v2q/3tx0SJL0zJA2SklM8O2gEFACuTZR7wPL+cJLav30Z5Kkb2b1U82I6j4eEYBgwfGlbGbUJdOu0BcWFurLL79UcnKyU3tycrK2bNlS5rzt27dXfHy8+vTpo/Xr1zv9LTMz02WZ/fr1K3OZBQUFstvtThPgCcI8ysOVeoQq6j0AAL5jWqA/ffq0ioqKFBsb69QeGxurvLw8t/PEx8frrbfeUnp6upYtW6YWLVqoT58+2rRpk6NPXl6eR8uUpLS0NFmtVsfUuHHjq1gzhBrCPCqKUI9QRL0HAMB3TL8HwmKxOL02DMOlrUSLFi3UokULx+vExETl5OTohRdeUI8ePSq1TEmaNm2aJk+e7Hhtt9sp8qgQwjw8VRLqJenNTYf01Md7JYnb7xH0qPcAAFQ9067Qx8TEKCwszOVM+qlTp1zOuJela9euOnjwoON1XFycx8uMjIxUdHS00wSUhzCPyuJKPUIJ9R4AAN8xLdBHRESoY8eOysjIcGrPyMhQUlJShZeTlZWl+Ph4x+vExESXZa5Zs8ajZQLlIczjahHqESqo9wAA+I6pt9xPnjxZKSkp6tSpkxITE/XWW2/p2LFjmjhxoqTLt8Z99913evfddyVJr7zyihISEtSmTRsVFhbq/fffV3p6utLT0x3LfOihh9SjRw/Nnj1bQ4YM0ccff6y1a9fqiy++MHNVEEII8/AWbr9HqKDeAwDgG6YG+uHDh+vMmTOaNWuWcnNz1bZtW61cuVJNmjSRJOXm5jr9Rm1hYaGmTJmi7777TjVq1FCbNm20YsUKDRw40NEnKSlJS5Ys0ZNPPqmnnnpKN9xwg5YuXaouXbqYuSoIEYR5eBuhHqGAeg8AgG+Y+jv0/orfpYU7hHmYid+pR3moTd7HNnWP34kGYBaOL2ULqN+hBwIJYR5m4zv1AAAA8DYCPUIeYR5VhVAPAAAAbyLQI6QR5lHVCPUAAADwFgI9QhZhHr5CqAcAAIA3EOgRkgjz8DVCPQAAAK4WgR4hhzAPf0GoBwAAwNUg0COkEObhbwj1AAAAqCwCPUIGYR7+ilAPAACAyiDQIyQQ5uHvCPUAAADwFIEeQY8wj0BBqAcAAIAnCPQIaoR5BBpCPQAAACqKQI+gRZhHoCLUAwAAoCII9AhKhHkEOkI9AAAAykOgR9AhzCNYEOoBAABQFgI9ggphHsGGUA8AAIDSEOgRNAjzCFaEegAAALhDoEdQIMwj2BHqAQAAcCUCPQIeYR6hglAPAACAXyPQI6AR5hFqCPUAAAAoQaBHwCLMI1QR6gEAACAR6BGgCPMIdYR6AAAAEOgRcAjzwGWEegAAgNBGoEdAIcwDzgj1AAAAoYtAj4BBmAfcI9QDAACEJtMD/dy5c9W0aVNFRUWpY8eO2rx5c6l9ly1bpttvv13XXnutoqOjlZiYqM8++8ypz8KFC2WxWFymCxcumL0q8CHCPFA2Qj18jXoPAEDVMzXQL126VKmpqXriiSeUlZWl7t27a8CAATp27Jjb/ps2bdLtt9+ulStX6ssvv1Tv3r01ePBgZWVlOfWLjo5Wbm6u0xQVFWXmqsCHCPNAxRDq4SvUewAAfKO6mQt/6aWXNG7cOI0fP16S9Morr+izzz7TvHnzlJaW5tL/lVdecXr93HPP6eOPP9a//vUvtW/f3tFusVgUFxdn5tDhJwjzgGdKQr0kvbnpkJ76eK8kKSUxwYejQrCj3gMA4BumXaEvLCzUl19+qeTkZKf25ORkbdmypULLKC4uVn5+vurVq+fU/uOPP6pJkyZq1KiRBg0a5HJG/0oFBQWy2+1OE/wfYR6oHK7UoypR7wEA8B3TAv3p06dVVFSk2NhYp/bY2Fjl5eVVaBkvvviifvrpJw0bNszR1rJlSy1cuFCffPKJFi9erKioKN166606ePBgqctJS0uT1Wp1TI0bN67cSqHKEOaBq0OoR1Wh3gMA4DumPxTPYrE4vTYMw6XNncWLF2vGjBlaunSpGjRo4Gjv2rWrRo8erZtvvlndu3fXBx98oBtvvFGvvfZaqcuaNm2abDabY8rJyan8CsF0hHnAOwj1qErUewAAqp5p36GPiYlRWFiYy9n5U6dOuZzFv9LSpUs1btw4ffjhh+rbt2+ZfatVq6bf/OY3ZZ6xj4yMVGRkZMUHD58hzAPexXfqYTbqPQAAvmPaFfqIiAh17NhRGRkZTu0ZGRlKSkoqdb7Fixdr7Nix+sc//qHf/va35b6PYRjKzs5WfHz8VY8ZvkWYB8zBlXqYiXoPAIDvmPqU+8mTJyslJUWdOnVSYmKi3nrrLR07dkwTJ06UdPnWuO+++07vvvuupMvF/d5779Wrr76qrl27Os7216hRQ1arVZI0c+ZMde3aVc2bN5fdbtdf/vIXZWdn669//auZqwKTEeYBc3GlHmai3gMA4BumBvrhw4frzJkzmjVrlnJzc9W2bVutXLlSTZo0kSTl5uY6/Ubtm2++qUuXLumPf/yj/vjHPzrax4wZo4ULF0qSzp07p9///vfKy8uT1WpV+/bttWnTJnXu3NnMVYGJCPNA1SDUwyzUewAAfMNiGIbh60FUNbvdLqvVKpvNpujoaF8PJ6QR5oGqZxiGnl+1X29uOiRJemZIG0K9H6A2eR/b1L3zhZfU+unPJEnfzOqnmhGmXt8BEEI4vpTNjLpk+lPugdIQ5gHf4Dv1AAAAwYFAD58gzAO+RagHAAAIfAR6VDnCPOAfCPUAAACBjUCPKkWYB/wLoR4AACBwEehRZQjzgH8i1AMAAAQmAj2qBGEe8G+EegAAgMBDoIfpCPNAYCDUAwAABBYCPUxFmAcCC6EeAAAgcBDoYRrCPBCYCPUAAACBgUAPUxDmgcBGqAcAAPB/BHp4HWEeCA6EegAAAP9GoIdXEeaB4EKoBwAA8F8EengNYR4IToR6AAAA/0Sgh1cQ5oHgRqgHAADwPwR6XDXCPBAaCPUAAAD+hUCPq0KYB0ILoR4AAMB/EOhRaYR5IDQR6gEAAPwDgR6VQpgHQhuhHgAAwPcI9PAYYR6ARKgHAADwNQI9PEKYB/BrhHoAAADfIdCjwgjzANwh1AMAAPgGgR4VQpgHUBZCPQAAQNUj0KNchHkAFUGoBwAAqFrVfT0A+JeiYkPbD/+gU/kX1KB2lK6rU0Oj5hPm4V5RUZE2b96s3NxcxcfHq3v37goLC/P1sOBDJaFekt7cdEhPfbxXkpSSmKDCS8V6L/OIjv5wXk3q1VRKYoIiqnNeGahqRcWG47+3HfpBPW68VmHVLD4cEfwZtR6eKLxU7Pjvv//PEY3r3oxabzLTt+7cuXPVtGlTRUVFqWPHjtq8eXOZ/Tdu3KiOHTsqKipKzZo10xtvvOHSJz09Xa1bt1ZkZKRat26t5cuXmzX8kLJ6T666zV6ne97eqoeWZOuet7eq9wsbCPNwa9myZUpISFDv3r01cuRI9e7dWwkJCVq2bJmvhwYfc3elfvTftqrlU6v0zIp9ejfzqJ5ZsU8tn1qltJXf+Hi08BbqfWBYvSdXfV/a6Hh938Id6jZ7nVbvyfXhqOCvqPXwRNrKb9R+Vobj9ezPDlDrq4CpgX7p0qVKTU3VE088oaysLHXv3l0DBgzQsWPH3PY/fPiwBg4cqO7duysrK0uPP/64Jk2apPT0dEefzMxMDR8+XCkpKdq9e7dSUlI0bNgwbdu2zcxVCXqr9+Tqgfd3Kdd2wam9yLh8Fv/+Hs0I83BYtmyZ7rrrLh0/ftyp/bvvvtNdd91FoYdLqP/i32f0q4uCkqRiQ3pz02EKfRCg3geGklp/0l7g1J5nu6AH3t9FqIcTaj08kbbyG7256bCuKPXU+ipgMQzjyu3uNV26dFGHDh00b948R1urVq30u9/9TmlpaS79p06dqk8++UT79u1ztE2cOFG7d+9WZmamJGn48OGy2+1atWqVo0///v1Vt25dLV68uELjstvtslqtstlsio6OruzqBY2iYkPdZq9zCfO/Fm+N0hdTb+OWPKioqEgJCQkuBb6ExWJRo0aNdPjwYW7JgwouFqnFU6vL7FPNIu1/ZkDI35IXyLWJeu//yqv1Fklx1Hr8H2o9PFF4qVgtn1rlcuL+16j1l5lRl0zbooWFhfryyy+VnJzs1J6cnKwtW7a4nSczM9Olf79+/bRz505dvHixzD6lLVOSCgoKZLfbnSb8YvvhH8oM85KUa7ug7Yd/qKIRwZ9t3ry51AIvSYZhKCcnp9zbbREa3t96tNw+xYZ4eF4Ao94HhvJqvSFqPX5BrYcn3ss8UmaYl6j1ZjIt0J8+fVpFRUWKjY11ao+NjVVeXp7befLy8tz2v3Tpkk6fPl1mn9KWKUlpaWmyWq2OqXHjxpVZpaB1Kr/sMO9pPwS33NyK3ZJZ0X4Ibkd/OO/VfvA/1PvAQK2HJ6j18AS13rdMv+fBYnG+bcswDJe28vpf2e7pMqdNmyabzeaYcnJyKjz+UNCgdpRX+yG4xcfHe7UfgluTejW92g/+i3rv36j18AS1Hp6g1vuWaYE+JiZGYWFhLmfST5065XLGvURcXJzb/tWrV1f9+vXL7FPaMiUpMjJS0dHRThN+0blpPcVbo1TWN+birVHq3LRelY0J/qt79+5q1KhRqR+qLRaLGjdurO7du1fxyOCPbmtZ+rG5RDXL5Z+1Q2Ci3geG8mq9RdR6/IJaD0+kJCaovEdvUOvNY1qgj4iIUMeOHZWRkeHUnpGRoaSkJLfzJCYmuvRfs2aNOnXqpPDw8DL7lLZMlC+smkXTB7eWpFILfU9+oxb/JywsTK+++qok16tnJa9feeUVHpIDHT97XvcuKP+J5BO6Nw35h+QEMup9YPh1rb9SyZF8+uDW1HpIotbDMxHVq2lC96Zl9qHWm8fUrTp58mT97W9/0zvvvKN9+/bp4Ycf1rFjxzRx4kRJl2+Nu/feex39J06cqKNHj2ry5Mnat2+f3nnnHc2fP19Tpkxx9HnooYe0Zs0azZ49W/v379fs2bO1du1apaammrkqQa9/23jNG91BcVbnW+1qRVw+UC/ZkcODLOAwdOhQffTRR7ruuuuc2hs1aqSPPvpIQ4cO9dHI4C+Onz2ve97eqpwfflZC/Zoa1aWxy9n7ahbp/h5NNW2g+5CBwEG9DwwltT42OtKpPc4apXmjO6h/W26fxi+o9fDEtIGtdX+Ppi4XB6n15jP1Z+skae7cuZozZ45yc3PVtm1bvfzyy+rRo4ckaezYsTpy5Ig2bNjg6L9x40Y9/PDD2rt3rxo2bKipU6c6PhCU+Oijj/Tkk0/q0KFDuuGGG/Tss896dFDhZ2xKV1RsaPvhH3Qq/4Ia1I7SbxLqas5nB/TWpkOSpGeGtOF2GTgUFRVp8+bNys3NVXx8vLp3787ZeriE+cW/76p4aw0VXirWe5lHdPSH82pSr6ZSEhM4W/8rgV6bqPeBI//CRbWbsUaStGDsb9SDu/BQBmo9PHHufKFumXX57qqp/VpoXPdm1PpfMaMumR7o/REF3jOGYSht1X5CPYBylRbmUT5qk/exTd07X3hJrZ/+TJL0zax+qhlR3ccjAhAsOL6ULaB+hx7Bw2KxaNqAlvp9j2aSpKc+3svt9wBcEOYBAACqFoEeFUKoB1AWwjwAAEDVI9Cjwgj1ANwhzAMAAPgGgR4eIdQD+DXCPAAAgO8Q6OExQj0AiTAPAADgawR6VAqhHghthHkAAADfI9Cj0gj1QGgizAMAAPgHAj2uCqEeCC2EeQAAAP9BoMdVI9QDoYEwDwAA4F8I9PAKQj0Q3AjzAAAA/odAD68h1APBiTAPAADgnwj08CpCPRBcCPMAAAD+i0APryPUA8GBMA8AAODfCPQwBaEeCGyEeQAAAP9HoIdpCPVAYCLMAwAABAYCPUxFqAcCC2EeAAAgcBDoYTpCPRAYCPMAAACBhUCPKkGoB/wbYR4AACDwEOhRZQj1gH8izAMAAAQmAj2qFKEe8C+EeQAAgMBFoEeVI9QD/oEwDwAAENgI9PAJQj3gW4R5AACAwEegh88Q6gHfIMwDAAAEBwI9fIpQD1QtwjwAAEDwMC3Qnz17VikpKbJarbJarUpJSdG5c+dK7X/x4kVNnTpV7dq1U61atdSwYUPde++9OnHihFO/Xr16yWKxOE0jRowwazVQBQj1QNUgzMMM1HsAAHzHtEA/cuRIZWdna/Xq1Vq9erWys7OVkpJSav/z589r165deuqpp7Rr1y4tW7ZM3377re644w6XvhMmTFBubq5jevPNN81aDVQRQj1gLsI8zEK9BwDAd6qbsdB9+/Zp9erV2rp1q7p06SJJevvtt5WYmKgDBw6oRYsWLvNYrVZlZGQ4tb322mvq3Lmzjh07puuvv97RXrNmTcXFxZkxdPhQSaiXpLc2HdJTH++VJKUkJvhwVEDgI8zDLNR7AAB8y5Qr9JmZmbJarY7iLkldu3aV1WrVli1bKrwcm80mi8WiOnXqOLUvWrRIMTExatOmjaZMmaL8/Pwyl1NQUCC73e40wT9xpR7wLsI8zES9BwDAt0y5Qp+Xl6cGDRq4tDdo0EB5eXkVWsaFCxf02GOPaeTIkYqOjna0jxo1Sk2bNlVcXJz27NmjadOmaffu3S5n+38tLS1NM2fO9HxF4BNcqQe8gzAPs1HvAQDwLY+u0M+YMcPlATVXTjt37pR0OZRdyTAMt+1XunjxokaMGKHi4mLNnTvX6W8TJkxQ37591bZtW40YMUIfffSR1q5dq127dpW6vGnTpslmszmmnJwcT1YbPsCVeuDqEOZxNaj3AAAEBo+u0D/44IPlPmE2ISFBX331lU6ePOnyt++//16xsbFlzn/x4kUNGzZMhw8f1rp165zO1rvToUMHhYeH6+DBg+rQoYPbPpGRkYqMjCxzOfA/XKkHKocwj6tFvQcAIDB4FOhjYmIUExNTbr/ExETZbDZt375dnTt3liRt27ZNNptNSUlJpc5XUtwPHjyo9evXq379+uW+1969e3Xx4kXFx8dXfEUQMAj1gGcI8/AG6j0AAIHBlIfitWrVSv3799eECRO0detWbd26VRMmTNCgQYOcnnjbsmVLLV++XJJ06dIl3XXXXdq5c6cWLVqkoqIi5eXlKS8vT4WFhZKk//3f/9WsWbO0c+dOHTlyRCtXrtTdd9+t9u3b69ZbbzVjVeAHuP0eqBjCPKoa9R4AAN8y7XfoFy1apHbt2ik5OVnJycm66aab9N577zn1OXDggGw2myTp+PHj+uSTT3T8+HHdcsstio+Pd0wlT8qNiIjQ559/rn79+qlFixaaNGmSkpOTtXbtWoWFhZm1KvADhHqgbIR5+Ar1HgAA3zHlKfeSVK9ePb3//vtl9jEMw/HfCQkJTq/dady4sTZu3OiV8SHwcPs94B5hHr5EvQcAwHdMu0IPmIEr9YAzwjwAAEDoItAj4BDqgcsI8wAAAKGNQI+ARKhHqCPMAwAAgECPgEWoR6gizAMAAEAi0CPAEeoRagjzAAAAKEGgR8Aj1CNUEOYBAADwawR6BAVCPYIdYR4AAABXItAjaBDqEawI8wAAAHCHQI+gQqhHsCHMAwAAoDQEegQdQj2CBWEeAAAAZSHQIygR6hHoCPMAAAAoD4EeQYtQj0BFmAcAAEBFEOgR1Aj1CDSEeQAAAFQUgR5Bj1CPQEGYBwAAgCcI9AgJhHr4O8I8AAAAPEWgR8gg1MNfEeYBAABQGQR6hBRCPfwNYR4AAACVRaBHyCHUw18Q5gEAAHA1CPQISYR6+BphHgAAAFeLQI+QRaiHrxDmAQAA4A0EeoQ0Qj2qGmEeAAAA3kKgR8gj1KOqEOYBAADgTQR6QIR6mI8wDwAAAG8j0AP/h1APsxDmAQAAYAbTAv3Zs2eVkpIiq9Uqq9WqlJQUnTt3rsx5xo4dK4vF4jR17drVqU9BQYH+9Kc/KSYmRrVq1dIdd9yh48ePm7UaCDGEengbYR7BjnoPAIDvmBboR44cqezsbK1evVqrV69Wdna2UlJSyp2vf//+ys3NdUwrV650+ntqaqqWL1+uJUuW6IsvvtCPP/6oQYMGqaioyKxVQYgh1MNbCPMIBdR7AAB8p7oZC923b59Wr16trVu3qkuXLpKkt99+W4mJiTpw4IBatGhR6ryRkZGKi4tz+zebzab58+frvffeU9++fSVJ77//vho3bqy1a9eqX79+3l8ZhKSSUC9Jb206pKc+3itJSklM8OGoEEgI8wgF1HsAAHzLlCv0mZmZslqtjuIuSV27dpXVatWWLVvKnHfDhg1q0KCBbrzxRk2YMEGnTp1y/O3LL7/UxYsXlZyc7Ghr2LCh2rZtW+ZyCwoKZLfbnSagPFypR2UR5hEqqPcAAPiWKYE+Ly9PDRo0cGlv0KCB8vLySp1vwIABWrRokdatW6cXX3xRO3bs0G233aaCggLHciMiIlS3bl2n+WJjY8tcblpamuO7fVarVY0bN67kmiHUEOrhKcI8Qgn1HgAA3/Io0M+YMcPlITZXTjt37pR0OQhdyTAMt+0lhg8frt/+9rdq27atBg8erFWrVunbb7/VihUryhxXecudNm2abDabY8rJyangGgOEelQcYR7BgnoPAEBg8Og79A8++KBGjBhRZp+EhAR99dVXOnnypMvfvv/+e8XGxlb4/eLj49WkSRMdPHhQkhQXF6fCwkKdPXvW6az9qVOnlJSUVOpyIiMjFRkZWeH3Ba7Ed+pRHsI8ggn1HgCAwOBRoI+JiVFMTEy5/RITE2Wz2bR9+3Z17txZkrRt2zbZbLYyC/GVzpw5o5ycHMXHx0uSOnbsqPDwcGVkZGjYsGGSpNzcXO3Zs0dz5szxZFUAjxHqURrCPIIN9R4AgMBgynfoW7Vqpf79+2vChAnaunWrtm7dqgkTJmjQoEFOT7xt2bKlli9fLkn68ccfNWXKFGVmZurIkSPasGGDBg8erJiYGN15552SJKvVqnHjxumRRx7R559/rqysLI0ePVrt2rVzPAUXMBO33+NKhHmEMuo9AAC+ZcrP1knSokWLNGnSJMcTau+44w69/vrrTn0OHDggm80mSQoLC9PXX3+td999V+fOnVN8fLx69+6tpUuXqnbt2o55Xn75ZVWvXl3Dhg3Tzz//rD59+mjhwoUKCwsza1UAJ1ypRwnCPEC9BwDAlyyGYRi+HkRVs9vtslqtstlsio6O9vVwEKAMw1Daqv16a9MhSdIzQ9oQ6kMIYR7eRm3yPrape+cLL6n1059Jkr6Z1U81I0y7vgMgxHB8KZsZdcmUW+6BUMDt96GLMA8AAAB/QKAHrgKhPvQQ5gEAAOAvCPTAVSLUhw7CPAAAAPwJgR7wAkJ98CPMAwAAwN8Q6AEvIdQHL8I8AAAA/BGBHvAiQn3wIcwDAADAXxHoAS8j1AcPwjwAAAD8GYEeMAGhPvAR5gEAAODvCPSASQj1gYswDwAAgEBAoAdMRKgPPIR5AAAABAoCPWAyQn3gIMwDAAAgkBDogSpAqPd/hHkAAAAEGgI9UEUI9f6LMA8AAIBARKAHqhCh3v8Q5gEAABCoCPRAFSPU+w/CPAAAAAIZgR7wAUK97xHmAQAAEOgI9ICPEOp9hzAPAACAYECgB3yIUF/1CPMAAAAIFgR6wMcI9VWHMA8AAIBgQqAH/ACh3nyEeQAAAAQbAj3gJwj15iHMAwAAIBgR6AE/Qqj3PsI8AAAAghWBHvAzhHrvIcwDAAAgmBHoAT9EqL96hHkAAAAEO9MC/dmzZ5WSkiKr1Sqr1aqUlBSdO3euzHksFovb6b//+78dfXr16uXy9xEjRpi1GoDPEOorjzAPVB3qPQAAvlPdrAWPHDlSx48f1+rVqyVJv//975WSkqJ//etfpc6Tm5vr9HrVqlUaN26c/vM//9OpfcKECZo1a5bjdY0afFBHcCoJ9ZL01qZDeurjvZKklMQEH47KvxHmgapFvQcAwHdMCfT79u3T6tWrtXXrVnXp0kWS9PbbbysxMVEHDhxQixYt3M4XFxfn9Prjjz9W79691axZM6f2mjVruvQFghWhvuII80DVot4DAOBbptxyn5mZKavV6ijuktS1a1dZrVZt2bKlQss4efKkVqxYoXHjxrn8bdGiRYqJiVGbNm00ZcoU5efnl7msgoIC2e12pwkIJNx+Xz7CPFD1qPcAAPiWKVfo8/Ly1KBBA5f2Bg0aKC8vr0LL+Pvf/67atWtr6NChTu2jRo1S06ZNFRcXpz179mjatGnavXu3MjIySl1WWlqaZs6c6dlKAH6GK/WlI8wDvkG9BwDAtzy6Qj9jxoxSH2RTMu3cuVPS5fBxJcMw3La7884772jUqFGKiopyap8wYYL69u2rtm3basSIEfroo4+0du1a7dq1q9RlTZs2TTabzTHl5OR4sNaA/+BKvSvCPOB91HsAAAKDR1foH3zwwXKfMJuQkKCvvvpKJ0+edPnb999/r9jY2HLfZ/PmzTpw4ICWLl1abt8OHTooPDxcBw8eVIcOHdz2iYyMVGRkZLnLAgIBV+p/QZgHzEG9BwAgMHgU6GNiYhQTE1Nuv8TERNlsNm3fvl2dO3eWJG3btk02m01JSUnlzj9//nx17NhRN998c7l99+7dq4sXLyo+Pr78FQCCBKGeMA+YiXoPAEBgMOWheK1atVL//v01YcIEbd26VVu3btWECRM0aNAgpyfetmzZUsuXL3ea126368MPP9T48eNdlvu///u/mjVrlnbu3KkjR45o5cqVuvvuu9W+fXvdeuutZqwK4LdC+fZ7wjzgH6j3AAD4limBXrr8ZNp27dopOTlZycnJuummm/Tee+859Tlw4IBsNptT25IlS2QYhu655x6XZUZEROjzzz9Xv3791KJFC02aNEnJyclau3atwsLCzFoVwG+FYqgnzAP+hXoPAIDvWAzDMHw9iKpmt9tltVpls9kUHR3t6+EAV80wDKWt2q+3Nh2SJD0zpE1Q3n5PmEcwozZ5H9vUvfOFl9T66c8kSd/M6qeaEab86BGAEMTxpWxm1CXTrtADqDqhcKWeMA8AAAA4I9ADQSKYQz1hHgAAAHBFoAeCSDCGesI8AAAA4B6BHggywRTqCfMAAABA6Qj0QBAKhlBPmAcAAADKRqAHglQgh3rCPAAAAFA+Aj0QxAIx1BPmAQAAgIoh0ANBLpBCPWEeAAAAqDgCPRACAiHUE+YBAAAAzxDogRDhz6GeMA8AAAB4jkAPhBB/DPWEeQAAAKByCPRAiPGnUE+YBwAAACqPQA+EIH8I9YR5AAAA4OoQ6IEQ5ctQT5gHAAAArh6BHghhvgj1hHkAAADAOwj0QIirylBPmAcAAAC8h0APoEpCPWEeAAAA8C4CPQBJ5oZ6wjwAAADgfQR6AA5mhHrCPAAAAGAOAj0AJ94M9YR5AAAAwDwEegAuvBHqCfMAAACAuQj0ANy6mlBPmAcAAADMR6AHUKrKhHrCPAAAAFA1qvt6AAD8W0mol6S3Nh3SUx/vlSSlJCaoqNjQ9sM/6FT+BTWoHaXr6tTQqPmEeQDwZ0XFhuO/tx36QT1uvFZh1Sw+HBH82ZW1vnPTeuwvKBXHl6pn2hX6Z599VklJSapZs6bq1KlToXkMw9CMGTPUsGFD1ahRQ7169dLevXud+hQUFOhPf/qTYmJiVKtWLd1xxx06fvy4CWsAoIS7K/WPpX+lbrPX6Z63t+qhJdm65+2t6v3CBsI8EGKo94Fl9Z5c9X1po+P1fQt3qNvsdVq9J9eHo4K/Wr0n16XWs7+gNBxffMO0QF9YWKi7775bDzzwQIXnmTNnjl566SW9/vrr2rFjh+Li4nT77bcrPz/f0Sc1NVXLly/XkiVL9MUXX+jHH3/UoEGDVFRUZMZqAPg/V4b6JTtylGu74NSnyLh8Vvb+Hs0I80CIoN4HjtV7cvXA+7t00l7g1J5nu6AH3t/Fh244Kdlfrqz17C9wh+OL71gMwzDK71Z5CxcuVGpqqs6dO1dmP8Mw1LBhQ6Wmpmrq1KmSLp+dj42N1ezZs3X//ffLZrPp2muv1Xvvvafhw4dLkk6cOKHGjRtr5cqV6tevX4XGZLfbZbVaZbPZFB0dfVXrB4SaS0XFunnmGv1UWPqH6nhrlL6Yehu3WAEeCPTaRL33b0XFhrrNXucSzn4tNjpSayf35NgNFRUb6vvSRpdw9mvsLyhR3v5ikRTHZ0NJ5tQlv/kO/eHDh5WXl6fk5GRHW2RkpHr27KktW7bo/vvv15dffqmLFy869WnYsKHatm2rLVu2lFrgCwoKVFDwyw5mt9vNWxEgyO04crbMMC9JubYL2n74ByXeUL+KRgUgUFDvfWP74R/KDPOSdNJeoHYz1lTRiBDo2F9QUYb4bGgmv3nKfV5eniQpNjbWqT02Ntbxt7y8PEVERKhu3bql9nEnLS1NVqvVMTVu3NjLowdCx6n8sj8QetoPQGih3vsGx2QAvsZxyBweXaGfMWOGZs6cWWafHTt2qFOnTpUekMXifBuGYRgubVcqr8+0adM0efJkx2u73U6RByqpQe0or/YD4H+o98GnosfkBWN/oy7N6pk8Gvi7bYd+0H0Ld5Tbj/0FUsX3Fz4bmsOjQP/ggw9qxIgRZfZJSEio1EDi4uIkXT4rHx8f72g/deqU4yx+XFycCgsLdfbsWaez9qdOnVJSUlKpy46MjFRkZGSlxgXAWeem9RRvjVKe7YLcPYCj5HtSnZtS4IFARb0PPhU9dvMTU5CkHjdey/6CCqvo/sJnQ3N4dMt9TEyMWrZsWeYUFVW5My9NmzZVXFycMjIyHG2FhYXauHGjo3h37NhR4eHhTn1yc3O1Z8+eMgs8AO8Jq2bR9MGtJV0+QP9ayevpg1tT4IEARr0PPhy74Qn2F3iC/cW3TPsO/bFjx5Sdna1jx46pqKhI2dnZys7O1o8//ujo07JlSy1fvlzS5VvvUlNT9dxzz2n58uXas2ePxo4dq5o1a2rkyJGSJKvVqnHjxumRRx7R559/rqysLI0ePVrt2rVT3759zVoVAFfo3zZe80Z3UJzV+QN9nDVK80Z3UP+28aXMCSDYUO8DB8dueIL9BZ5gf/Ed055y//TTT+vvf/+743X79u0lSevXr1evXr0kSQcOHJDNZnP0efTRR/Xzzz/rD3/4g86ePasuXbpozZo1ql27tqPPyy+/rOrVq2vYsGH6+eef1adPHy1cuFBhYWFmrQoAN/q3jdftreO0/fAPOpV/QQ1qX76VirOvQGih3gcWjt3wBPsLPMH+4hum/w69P+J3aQEA/oba5H1sUwCAPzGjLvnNz9YBAAAAAICKI9ADAAAAABCACPQAAAAAAAQgAj0AAAAAAAGIQA8AAAAAQAAi0AMAAAAAEIBM+x16f1byS312u93HIwEA4LKSmhSCvyZrGuo9AMCfmFHrQzLQ5+fnS5IaN27s45EAAOAsPz9fVqvV18MICtR7AIA/8mattxgheCmguLhYJ06ckGEYuv7665WTk6Po6GhfD8tv2O12NW7cmO1yBbaLe2wX99gu7rFd3CvZLt98841atGihatX4Rpw3FBcX68CBA2rdujX73BX4f9E9tot7bBf32C7usV3cM6vWh+QV+mrVqqlRo0aOWx6io6PZ2dxgu7jHdnGP7eIe28U9tot71113HWHei6pVq6brrrtOEvtcadgu7rFd3GO7uMd2cY/t4p63az2fGgAAAAAACEAEegAAAAAAAlBIB/rIyEhNnz5dkZGRvh6KX2G7uMd2cY/t4h7bxT22i3tsF/Owbd1ju7jHdnGP7eIe28U9tot7Zm2XkHwoHgAAAAAAgS6kr9ADAAAAABCoCPQAAAAAAAQgAj0AAAAAAAGIQA8AAAAAQAAKqUD/7LPPKikpSTVr1lSdOnUqNM/YsWNlsVicpq5du5o70CpWme1iGIZmzJihhg0bqkaNGurVq5f27t1r7kCr2NmzZ5WSkiKr1Sqr1aqUlBSdO3euzHmCdX+ZO3eumjZtqqioKHXs2FGbN28us//GjRvVsWNHRUVFqVmzZnrjjTeqaKRVy5PtsmHDBpd9w2KxaP/+/VU4YnNt2rRJgwcPVsOGDWWxWPTPf/6z3HlCYV/xdLuEwr5iNuq9e9R796j3l1Hr3aPWO6PWu+fLWh9Sgb6wsFB33323HnjgAY/m69+/v3Jzcx3TypUrTRqhb1Rmu8yZM0cvvfSSXn/9de3YsUNxcXG6/fbblZ+fb+JIq9bIkSOVnZ2t1atXa/Xq1crOzlZKSkq58wXb/rJ06VKlpqbqiSeeUFZWlrp3764BAwbo2LFjbvsfPnxYAwcOVPfu3ZWVlaXHH39ckyZNUnp6ehWP3FyebpcSBw4ccNo/mjdvXkUjNt9PP/2km2++Wa+//nqF+ofKvuLpdikRzPuK2aj37lHv3aPeU+tLQ613Ra13z6e13ghBCxYsMKxWa4X6jhkzxhgyZIip4/EXFd0uxcXFRlxcnPH888872i5cuGBYrVbjjTfeMHGEVeebb74xJBlbt251tGVmZhqSjP3795c6XzDuL507dzYmTpzo1NayZUvjsccec9v/0UcfNVq2bOnUdv/99xtdu3Y1bYy+4Ol2Wb9+vSHJOHv2bBWMzvckGcuXLy+zT6jsK79Wke0SavuKmaj37lHvf0G9v4xa7x61vmzUevequtaH1BX6ytqwYYMaNGigG2+8URMmTNCpU6d8PSSfOnz4sPLy8pScnOxoi4yMVM+ePbVlyxYfjsx7MjMzZbVa1aVLF0db165dZbVay13HYNpfCgsL9eWXXzr9W0tScnJyqdshMzPTpX+/fv20c+dOXbx40bSxVqXKbJcS7du3V3x8vPr06aP169ebOUy/Fwr7ytVgX6l6wXT89gbqfWjUe2q9e9R67wiFfeVqeGNfIdCXY8CAAVq0aJHWrVunF198UTt27NBtt92mgoICXw/NZ/Ly8iRJsbGxTu2xsbGOvwW6vLw8NWjQwKW9QYMGZa5jsO0vp0+fVlFRkUf/1nl5eW77X7p0SadPnzZtrFWpMtslPj5eb731ltLT07Vs2TK1aNFCffr00aZNm6piyH4pFPaVymBf8Y1gO357A/U+NOo9td49ar13hMK+Uhne3FeqmzC+KjVjxgzNnDmzzD47duxQp06dKrX84cOHO/67bdu26tSpk5o0aaIVK1Zo6NChlVpmVTB7u0iSxWJxem0Yhkubv6nodpFc108qfx0DdX8pj6f/1u76u2sPdJ5slxYtWqhFixaO14mJicrJydELL7ygHj16mDpOfxYq+4on2Ffco967R713j3rvOWq9e9T6qxcq+4onvLmvBHygf/DBBzVixIgy+yQkJHjt/eLj49WkSRMdPHjQa8s0g5nbJS4uTtLlM27x8fGO9lOnTrmcgfM3Fd0uX331lU6ePOnyt++//96jdQyU/aU0MTExCgsLczkTXda/dVxcnNv+1atXV/369U0ba1WqzHZxp2vXrnr//fe9PbyAEQr7ireE+r4iUe9LQ713j3pfcdR696j13hEK+4q3VHZfCfhAHxMTo5iYmCp7vzNnzignJ8epsPkjM7dL06ZNFRcXp4yMDLVv317S5e8Zbdy4UbNnzzblPb2lotslMTFRNptN27dvV+fOnSVJ27Ztk81mU1JSUoXfL1D2l9JERESoY8eOysjI0J133uloz8jI0JAhQ9zOk5iYqH/9619ObWvWrFGnTp0UHh5u6nirSmW2iztZWVkBu294QyjsK94S6vuKRL0vDfXePep9xVHr3aPWe0co7CveUul95aofqxdAjh49amRlZRkzZ840rrnmGiMrK8vIysoy8vPzHX1atGhhLFu2zDAMw8jPzzceeeQRY8uWLcbhw4eN9evXG4mJicZ1111n2O12X62G13m6XQzDMJ5//nnDarUay5YtM77++mvjnnvuMeLj44Nqu/Tv39+46aabjMzMTCMzM9No166dMWjQIKc+obC/LFmyxAgPDzfmz59vfPPNN0ZqaqpRq1Yt48iRI4ZhGMZjjz1mpKSkOPofOnTIqFmzpvHwww8b33zzjTF//nwjPDzc+Oijj3y1CqbwdLu8/PLLxvLly41vv/3W2LNnj/HYY48Zkoz09HRfrYLX5efnO44fkoyXXnrJyMrKMo4ePWoYRujuK55ul1DYV8xGvXePeu8e9Z5aXxpqvStqvXu+rPUhFejHjBljSHKZ1q9f7+gjyViwYIFhGIZx/vx5Izk52bj22muN8PBw4/rrrzfGjBljHDt2zDcrYBJPt4thXP4pm+nTpxtxcXFGZGSk0aNHD+Prr7+u+sGb6MyZM8aoUaOM2rVrG7Vr1zZGjRrl8tMSobK//PWvfzWaNGliREREGB06dDA2btzo+NuYMWOMnj17OvXfsGGD0b59eyMiIsJISEgw5s2bV8UjrhqebJfZs2cbN9xwgxEVFWXUrVvX6Natm7FixQofjNo8JT/BcuU0ZswYwzBCd1/xdLuEwr5iNuq9e9R796j3l1Hr3aPWO6PWu+fLWm8xjP97KgEAAAAAAAgY/GwdAAAAAAABiEAPAAAAAEAAItADAAAAABCACPQAAAAAAAQgAj0AAAAAAAGIQA8AAAAAQAAi0AMAAAAAEIAI9AAAAAAABCACPQAAAAAAAYhADwAAAABAACLQAwAAAAAQgAj0AAAAAAAEoP8PyPvlRVHBeUcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1200x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[12,4])\n",
    "plt.subplot(1,2,1)\n",
    "plt.title(\"Diamond stencil\"); plt.axis('equal')\n",
    "PlotStencil2(diamond2)\n",
    "plt.subplot(1,2,2)\n",
    "plt.title(\"Square stencil\"); plt.axis('equal')\n",
    "PlotStencil2(box2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute th minimum scalar product between consecutive elements of a stencil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.569709Z",
     "iopub.status.busy": "2024-04-30T08:46:15.569585Z",
     "iopub.status.idle": "2024-04-30T08:46:15.571729Z",
     "shell.execute_reply": "2024-04-30T08:46:15.571497Z"
    }
   },
   "outputs": [],
   "source": [
    "def MinStencilScal2(stencil,m):\n",
    "    s = np.array(stencil[:-1]).T\n",
    "    m=np.expand_dims(m,axis=2)\n",
    "    return lp.dot_VAV(s,m,np.roll(s,axis=1,shift=1)).min()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both the diamond and the box stencil are *acute*, a.k.a. the minimal scalar product is non-negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.573100Z",
     "iopub.status.busy": "2024-04-30T08:46:15.573001Z",
     "iopub.status.idle": "2024-04-30T08:46:15.575270Z",
     "shell.execute_reply": "2024-04-30T08:46:15.575031Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1.0)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MinStencilScal2(diamond2,np.eye(2)), MinStencilScal2(box2,np.eye(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, as soon as some non-diagonal anisotropy is introduced, the diamond stencil looses acuteness.\n",
    "The matrix corresponding to a needle-like Riemannian norm has two eigenspaces: the line spanned by the given direction, and the orthogonal space. The corresponding two eigenvalues are prescribed as parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.576687Z",
     "iopub.status.busy": "2024-04-30T08:46:15.576586Z",
     "iopub.status.idle": "2024-04-30T08:46:15.578234Z",
     "shell.execute_reply": "2024-04-30T08:46:15.578015Z"
    }
   },
   "outputs": [],
   "source": [
    "riemann2_2 = Riemann.needle(direction2,1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.579632Z",
     "iopub.status.busy": "2024-04-30T08:46:15.579552Z",
     "iopub.status.idle": "2024-04-30T08:46:15.581766Z",
     "shell.execute_reply": "2024-04-30T08:46:15.581533Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.299038105676658, 0.45096189432334155)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MinStencilScal2(diamond2,riemann2_2.m), MinStencilScal2(box2,riemann2_2.m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The box stencil looses acuteness when anisotropy becomes more pronounced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.583088Z",
     "iopub.status.busy": "2024-04-30T08:46:15.582997Z",
     "iopub.status.idle": "2024-04-30T08:46:15.584543Z",
     "shell.execute_reply": "2024-04-30T08:46:15.584323Z"
    }
   },
   "outputs": [],
   "source": [
    "riemann2_3 = Riemann.needle(direction2,1,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.585822Z",
     "iopub.status.busy": "2024-04-30T08:46:15.585741Z",
     "iopub.status.idle": "2024-04-30T08:46:15.588023Z",
     "shell.execute_reply": "2024-04-30T08:46:15.587798Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.4641016151377544, -0.4641016151377553)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MinStencilScal2(diamond2,riemann2_3.m), MinStencilScal2(box2,riemann2_3.m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.589374Z",
     "iopub.status.busy": "2024-04-30T08:46:15.589293Z",
     "iopub.status.idle": "2024-04-30T08:46:15.727244Z",
     "shell.execute_reply": "2024-04-30T08:46:15.726960Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+oAAAF0CAYAAAC0bcPIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHdElEQVR4nO3dd3hUZd7/8c8QkgktQwmkrJHiYkCKQiIQeGiCoQh2AWEjWFBsiMBPioWyu7JgdxEsi7JYEBVYG6BRaT4EBAwoiogKgpqAtAmgBAj37w8eZh0mbWBO5szM+3Vdc13OyX1O7nMy8vl+zzk5cRhjjAAAAAAAgC1UCvYEAAAAAADAf9GoAwAAAABgIzTqAAAAAADYCI06AAAAAAA2QqMOAAAAAICN0KgDAAAAAGAjNOoAAAAAANgIjToAAAAAADZCow4AAAAAgI3QqAOSZs+eLYfD4fWqW7euunTpovfeey/Y0yvWokWLNHHixGBPw4vD4fCa07Jly+RwOLRs2bJS1yvvODt77bXX9OSTTwZ7GgAQtsjqwAjXrH744Yf1n//8J9jTKNUvv/yiiRMnasOGDcGeCkIAjTrwBy+99JJycnK0atUqPf/884qKilLfvn317rvvBntqPhYtWqRJkyYFexpecnJydMstt/i9XuvWrZWTk6PWrVtbMKuKQaMOABWDrD474ZrVodKoT5o0iUYd5VI52BMA7KR58+ZKT0/3vO/Zs6dq1aqluXPnqm/fvkGcWWho167dGa0XFxd3xusCACILWX12yGogNHBFHShFbGysYmJiFB0d7bV83759uuOOO/SnP/1JMTExatSoke6//34VFhZKko4cOaJWrVrpz3/+s9xut2e9/Px8JSYmqkuXLioqKirx+/72228aPXq0GjZsqNjYWNWuXVvp6emaO3euJGnIkCF65plnJMnrFsDt27dLkowxmjFjhi666CJVqVJFtWrV0rXXXqsffvjB6/t06dJFzZs319q1a9WxY0dVrVpVjRo10j/+8Q+dOHHCa+yBAwc0atQoNWrUSE6nU/Xq1VPv3r31zTffeMacfjtdeRV3O90PP/ygAQMGKDk5WU6nUwkJCerWrVuZZ6HXrVunAQMGqEGDBqpSpYoaNGig66+/Xj/++KPP2J9//lm33nqrUlJSFBMTo+TkZF177bXatWuXpP/eZnnquJY03y5duuj999/Xjz/+6PXzOOXo0aP629/+piZNmsjpdKpu3bq68cYb9euvv/p9rAAA3sjq/wqVrB4yZIiqV6+u7777Tr1791b16tWVkpKiUaNGeX4+p5T1czy1T4cPH9a///1vz3Hu0qVLqXOYNGmS2rZtq9q1aysuLk6tW7fWrFmzZIzxGfvaa68pIyND1atXV/Xq1XXRRRdp1qxZnq83aNBAQ4YM8VmvS5cunnksW7ZMF198sSTpxhtv9Mzzjz+LdevW6fLLL1ft2rUVGxurVq1a6Y033ih1PxC+uKIO/EFRUZGOHz8uY4x27dqlRx55RIcPH9bAgQM9Y44cOaKuXbvq+++/16RJk9SyZUutXLlSU6ZM0YYNG/T+++8rNjZWb7zxhtLS0nTTTTdp/vz5OnHihAYNGiRjjObOnauoqKgS5zFy5Ei9/PLL+tvf/qZWrVrp8OHD2rRpk/bu3StJevDBB3X48GG99dZbysnJ8ayXlJQkSbrttts0e/ZsDR8+XFOnTtW+ffs0efJktW/fXhs3blRCQoJnnfz8fA0aNEijRo3ShAkTtHDhQo0bN07Jycm64YYbJEkHDx7U//zP/2j79u0aM2aM2rZtq0OHDmnFihXKy8tTkyZNAvpzkKTevXurqKhI06ZN07nnnqs9e/Zo1apVOnDgQKnrbd++XampqRowYIBq166tvLw8zZw5UxdffLG+/vprxcfHSzrZpF988cU6duyYxo8fr5YtW2rv3r364IMPtH//fq9jVJYZM2bo1ltv1ffff6+FCxd6fe3EiRO64oortHLlSt13331q3769fvzxR02YMEFdunTRunXrVKVKFb+PDwBEKrI69LNako4dO6bLL79cN998s0aNGqUVK1bor3/9q1wulx566CFJ5fs5Sidv57/kkkvUtWtXPfjgg5JO3gFQmu3bt+u2227TueeeK0lavXq17r77bv3888+e7y9JDz30kP7617/q6quv1qhRo+RyubRp06ZiLwCUpnXr1nrppZd044036oEHHtBll10mSTrnnHMkSUuXLlXPnj3Vtm1bPfvss3K5XHr99dfVv39//fbbb8WeCECYMwDMSy+9ZCT5vJxOp5kxY4bX2GeffdZIMm+88YbX8qlTpxpJ5sMPP/QsmzdvnpFknnzySfPQQw+ZSpUqeX29JM2bNzdXXnllqWPuvPNOU9z/wjk5OUaSeeyxx7yW79y501SpUsXcd999nmWdO3c2ksyaNWu8xl5wwQWmR48enveTJ082kkx2dnapc5JkJkyY4Hm/dOlSI8ksXbq01PVOH7dnzx7PcTtbx48fN4cOHTLVqlUzTz31lGf5TTfdZKKjo83XX39d4rqnPhfbtm0rdb7GGHPZZZeZ+vXr+2xj7ty5RpKZP3++1/K1a9caST6fLwBA8cjq8MnqwYMHF/vz6d27t0lNTfW89+fnWK1aNTN48GC/52KMMUVFRebYsWNm8uTJpk6dOubEiRPGGGN++OEHExUVZQYNGlTq+vXr1y/2e3fu3Nl07tzZ8/5U9r/00ks+Y5s0aWJatWpljh075rW8T58+JikpyRQVFfm9Xwht3PoO/MGcOXO0du1arV27VosXL9bgwYN15513avr06Z4xn3zyiapVq6Zrr73Wa91TZzo//vhjz7J+/frp9ttv1//7f/9Pf/vb3zR+/HhdeumlZc6jTZs2Wrx4scaOHatly5bp999/L/c+vPfee3I4HPrLX/6i48ePe16JiYm68MILfZ7WmpiYqDZt2ngta9mypdeZ4sWLF+v8889X9+7dyz2Ps1G7dm2dd955euSRR/T4448rNzfX5/a+khw6dEhjxozRn//8Z1WuXFmVK1dW9erVdfjwYW3evNkzbvHixeratauaNm1q1W5IOvnzqFmzpvr27ev187jooouUmJho26fnAoBdkdUnhXJWSydvVz/9mQKn75M/P0d/ffLJJ+revbtcLpeioqIUHR2thx56SHv37tXu3bslSdnZ2SoqKtKdd955xt+nPL777jt98803GjRokCR5fSZ69+6tvLw8bdmyxdI5wH5o1IE/aNq0qdLT05Wenq6ePXvqueeeU2Zmpu677z7PbVx79+5VYmKi1+8gS1K9evVUuXJlzy1vp9x00006duyYKleurOHDh5drHk8//bTGjBmj//znP+ratatq166tK6+8Ulu3bi1z3V27dskYo4SEBEVHR3u9Vq9erT179niNr1Onjs82nE6nV8Hx66+/em7NqggOh0Mff/yxevTooWnTpql169aqW7euhg8froMHD5a67sCBAzV9+nTdcsst+uCDD/TZZ59p7dq1qlu3blD2adeuXTpw4IDn9yf/+MrPz/f5eQAASkdWnxTKWS1JVatWVWxsrNcyp9OpI0eOeN77+3Msr88++0yZmZmSpBdeeEH/+7//q7Vr1+r++++XJM9xPfUsGauP66ln44wePdrn83DHHXdIEvVCBOJ31IEytGzZUh988IG+/fZbtWnTRnXq1NGaNWtkjPEKjt27d+v48eOe34GWpMOHDysrK0vnn3++du3apVtuuUVvv/12md+zWrVqmjRpkiZNmqRdu3Z5ztj37dvX64EwxYmPj5fD4dDKlSvldDp9vl7csrLUrVtXP/30k9/rnY369et7HtTy7bff6o033tDEiRN19OhRPfvss8Wu43a79d5772nChAkaO3asZ3lhYaH27dvnNbY8+3SqgDj9wTb+hGV8fLzq1KmjJUuWFPv1GjVqlHtbAIDikdWhk9X+8Ofn6I/XX39d0dHReu+997xOFpz+593q1q0rSfrpp5+UkpJS4vZiY2N9agXpZL1QnjmeGjNu3DhdffXVxY5JTU0tczsIL1xRB8pw6smlp/6x7tatmw4dOuTzj/mcOXM8Xz9l2LBh2rFjhxYsWKBZs2bpnXfe0RNPPOHX909ISNCQIUN0/fXXa8uWLfrtt98k/TfET7/Vrk+fPjLG6Oeff/Zccfjjq0WLFn59f0nq1auXvv32W33yySd+rxsI559/vh544AG1aNFCn3/+eYnjHA6HjDE+Bc6//vUvnyf39urVS0uXLi31VrIGDRpIkr744guv5e+8847P2NOvbJzSp08f7d27V0VFRcX+PAheADh7ZHXoZLU//Pk5lpTDxXE4HKpcubLXwwJ///13vfzyy17jMjMzFRUVpZkzZ5a6vQYNGvjUCt9++61PjVHS5yE1NVWNGzfWxo0bi/08pKenc2I/AnFFHfiDTZs26fjx45JO3m61YMECZWdn66qrrlLDhg0lSTfccIOeeeYZDR48WNu3b1eLFi306aef6uGHH1bv3r09vxv2r3/9S6+88opeeuklNWvWTM2aNdNdd92lMWPGqEOHDj6/a/ZHbdu2VZ8+fdSyZUvVqlVLmzdv1ssvv6yMjAxVrVpVkjwhPnXqVPXq1UtRUVFq2bKlOnTooFtvvVU33nij1q1bp06dOqlatWrKy8vTp59+qhYtWuj222/367iMGDFC8+bN0xVXXKGxY8eqTZs2+v3337V8+XL16dNHXbt29ftYl+aLL77QXXfdpeuuu06NGzdWTEyMPvnkE33xxRdeV8pPFxcXp06dOumRRx5RfHy8GjRooOXLl2vWrFmqWbOm19jJkydr8eLF6tSpk8aPH68WLVrowIEDWrJkiUaOHKkmTZro4osvVmpqqkaPHq3jx4+rVq1aWrhwoT799FOf792iRQstWLBAM2fOVFpamipVqqT09HQNGDBAr776qnr37q177rlHbdq0UXR0tH766SctXbpUV1xxha666qqAHj8ACGdkdfFCJav9Ud6fo3TyWC9btkzvvvuukpKSVKNGjRJPhl922WV6/PHHNXDgQN16663au3evHn30UZ8T/Q0aNND48eP117/+Vb///ruuv/56uVwuff3119qzZ48mTZokScrKytJf/vIX3XHHHbrmmmv0448/atq0aZ4TR6ecd955qlKlil599VU1bdpU1atXV3JyspKTk/Xcc8+pV69e6tGjh4YMGaI//elP2rdvnzZv3qzPP/9cb775ZkCOKUJI8J5jB9hHcU+Sdblc5qKLLjKPP/64OXLkiNf4vXv3mmHDhpmkpCRTuXJlU79+fTNu3DjPuC+++MJUqVLF5wmgR44cMWlpaaZBgwZm//79Jc5n7NixJj093dSqVcs4nU7TqFEjc++995o9e/Z4xhQWFppbbrnF1K1b1zgcDp+nk7/44oumbdu2plq1aqZKlSrmvPPOMzfccINZt26dZ0znzp1Ns2bNfL7/4MGDfZ5gvn//fnPPPfeYc88910RHR5t69eqZyy67zHzzzTeeMQrQk2R37dplhgwZYpo0aWKqVatmqlevblq2bGmeeOIJc/z48VK39dNPP5lrrrnG1KpVy9SoUcP07NnTbNq0qdgnsu7cudPcdNNNJjEx0URHR5vk5GTTr18/s2vXLs+Yb7/91mRmZpq4uDhTt25dc/fdd5v333/fZ7/27dtnrr32WlOzZk3Pz+OUY8eOmUcffdRceOGFJjY21lSvXt00adLE3HbbbWbr1q2l7g8A4CSy2lsoZ/XgwYNNtWrVfJZPmDDB5yn5Zf0cT9mwYYPp0KGDqVq1qpHk9bT14rz44osmNTXV87ObMmWKmTVrVrF/7WXOnDnm4osv9mR4q1atvJ7cfuLECTNt2jTTqFEjExsba9LT080nn3zi89R3Y07+NZgmTZqY6Ohon5/Fxo0bTb9+/Uy9evVMdHS0SUxMNJdccol59tlnS90XhCeHMcZU5IkBAAAAAABQMn5HHQAAAAAAG6FRBwAAAADARmjUAQAAAACwEUsb9RUrVqhv375KTk6Ww+Hw+dMKxVm+fLnS0tIUGxurRo0aFfs3GOfPn68LLrhATqdTF1xwgRYuXGjB7AEAQFnIegAAAs/SRv3w4cO68MILNX369HKN37Ztm3r37q2OHTsqNzdX48eP1/DhwzV//nzPmJycHPXv319ZWVnauHGjsrKy1K9fP61Zs8aq3QAAACUg6wEACLwKe+q7w+HQwoULdeWVV5Y4ZsyYMXrnnXe0efNmz7Jhw4Zp48aNysnJkST1799fBQUFWrx4sWdMz549VatWLc2dO9ey+QMAgNKR9QAABEblYE/gj3JycpSZmem1rEePHpo1a5aOHTum6Oho5eTk6N577/UZ8+STT5a43cLCQhUWFnrenzhxQvv27VOdOnXkcDgCug8AAJwJY4wOHjyo5ORkVaoUvo+QsSrrJfIeAGBv/mS9rRr1/Px8JSQkeC1LSEjQ8ePHtWfPHiUlJZU4Jj8/v8TtTpkyRZMmTbJkzgAABNLOnTt1zjnnBHsalrEq6yXyHgAQGsqT9bZq1CX5nPE+dWf+H5cXN6a0M+Xjxo3TyJEjPe/dbrfOPfdc7dy5U3FxcYGYNgAAZ6WgoEApKSmqUaNGsKdiOSuyXiLvAQD25k/W26pRT0xM9Dlbvnv3blWuXFl16tQpdczpZ97/yOl0yul0+iyPi4sjuAEAthLut2hblfUSeQ8ACA3lyXpb/RJcRkaGsrOzvZZ9+OGHSk9PV3R0dKlj2rdvX2HzBAAAZ4asBwCgbJZeUT906JC+++47z/tt27Zpw4YNql27ts4991yNGzdOP//8s+bMmSPp5FNfp0+frpEjR2ro0KHKycnRrFmzvJ7wes8996hTp06aOnWqrrjiCr399tv66KOP9Omnn1q5KwAAoBhkPQAAFjAWWrp0qZHk8xo8eLAxxpjBgwebzp07e62zbNky06pVKxMTE2MaNGhgZs6c6bPdN99806Smppro6GjTpEkTM3/+fL/m5Xa7jSTjdrvPdNcAAAioUM0mu2a9MaF7TAEA4cmfXKqwv6NuJwUFBXK5XHK73fzOGgDAFsimwOOYAgDsxJ9cstXvqAMAAAAAEOlo1AEAAAAAsBEadQAAAAAAbIRGHQAAAAAAG6FRBwAAAADARmjUAQAAAACwERp1AAAAAABshEYdAAAAAAAboVEHAAAAAMBGaNQBAAAAALARGnUAAAAAAGyERh0AAAAAABuhUQcAAAAAwEZo1AEAAAAAsBEadQAAAAAAbIRGHQAAAAAAG6FRBwAAAADARmjUAQAAAACwERp1AAAAAABshEYdAAAAAAAboVEHAAAAAMBGaNQBAAAAALARGnUAAAAAAGyERh0AAAAAABuhUQcAAAAAwEZo1AEAAAAAsBEadQAAAAAAbMTyRn3GjBlq2LChYmNjlZaWppUrV5Y4dsiQIXI4HD6vZs2aecbMnj272DFHjhyxelcAAEAJyHsAAALH0kZ93rx5GjFihO6//37l5uaqY8eO6tWrl3bs2FHs+Keeekp5eXme186dO1W7dm1dd911XuPi4uK8xuXl5Sk2NtbKXQEAACUg7wEACCxLG/XHH39cN998s2655RY1bdpUTz75pFJSUjRz5sxix7tcLiUmJnpe69at0/79+3XjjTd6jXM4HF7jEhMTrdwNAABQCvIeAIDAsqxRP3r0qNavX6/MzEyv5ZmZmVq1alW5tjFr1ix1795d9evX91p+6NAh1a9fX+ecc4769Omj3NzcgM0bAACUH3kPAEDgVbZqw3v27FFRUZESEhK8lickJCg/P7/M9fPy8rR48WK99tprXsubNGmi2bNnq0WLFiooKNBTTz2lDh06aOPGjWrcuHGx2yosLFRhYaHnfUFBwRnsEQAAOB15DwBA4Fn+MDmHw+H13hjjs6w4s2fPVs2aNXXllVd6LW/Xrp3+8pe/6MILL1THjh31xhtv6Pzzz9c///nPErc1ZcoUuVwuzyslJeWM9gUAABSPvAcAIHAsa9Tj4+MVFRXlczZ99+7dPmfdT2eM0YsvvqisrCzFxMSUOrZSpUq6+OKLtXXr1hLHjBs3Tm632/PauXNn+XcEAACUiLwHACDwLGvUY2JilJaWpuzsbK/l2dnZat++fanrLl++XN99951uvvnmMr+PMUYbNmxQUlJSiWOcTqfi4uK8XgAA4OyR9wAABJ5lv6MuSSNHjlRWVpbS09OVkZGh559/Xjt27NCwYcMknTzz/fPPP2vOnDle682aNUtt27ZV8+bNfbY5adIktWvXTo0bN1ZBQYGefvppbdiwQc8884yVuwIAAEpA3gMAEFiWNur9+/fX3r17NXnyZOXl5al58+ZatGiR56mueXl5Pn9j1e12a/78+XrqqaeK3eaBAwd06623Kj8/Xy6XS61atdKKFSvUpk0bK3cFAACUgLwHACCwHMYYE+xJVLSCggK5XC653W5uiwMA2ALZFHgcUwCAnfiTS5Y/9R0AAAAAAJQfjToAAAAAADZCow4AAAAAgI3QqAMAAAAAYCM06gAAAAAA2AiNOgAAAAAANkKjDgAAAACAjVQO9gQAIFQ1evoxy7/HD8NHWf49AAAAYC806gBQjIpowsujrHnQyAMAAIQfGnUAEc0uDfmZKm7+NO8AAAChjUYdQMQI9aa8vE7fTxp3AACA0EKjDiAsRUpTXh5/PBY07QAAAPZHow4gLNCYlw9NOwAAgP3RqAMISTTmZ+/UMaRhBwAAsBcadQAhgcbcOjTsAAAA9kKjDsC2aM4rFg07AACAPdCoA7ANGnN7oGEHAAAILhp1AEFFc25fNOwAAADBQaMOoMLRnIeWRk8/RrMOAABQgWjUAVQImvPQRrMOAABQcWjUAViKBj18cCs8AAAVI9D1E9kdemjUAQQczXl44+o6AABnp6JrpZK+H3luXzTqAAKC5jyy0KwDAFA2u9dHf5wfuW4vNOoAzordAwjWoVkHAOC/Qr0momm3Fxp1AH4L9SBC4NCsAwAiVTjXQzyXJvho1AGUWzgHEs4czToAIBJEYh1ExgcPjTqAMkViMME/BDkAINxQ/5zE1fXgoFEHUCzCCf6iWQcAhDJqn9KR8xWrktXfYMaMGWrYsKFiY2OVlpamlStXljh22bJlcjgcPq9vvvnGa9z8+fN1wQUXyOl06oILLtDChQut3g0gYjR6+jGCCoDfyHsAoehU3UPtUz4cp4pjaaM+b948jRgxQvfff79yc3PVsWNH9erVSzt27Ch1vS1btigvL8/zaty4sedrOTk56t+/v7KysrRx40ZlZWWpX79+WrNmjZW7AoQ9QgqBwGcoMpH3AEIJzfnZ4bhVDIcxxli18bZt26p169aaOXOmZ1nTpk115ZVXasqUKT7jly1bpq5du2r//v2qWbNmsdvs37+/CgoKtHjxYs+ynj17qlatWpo7d2655lVQUCCXyyW32624uDj/dgoIM/xjCytwa5z/QjmbyHsAdke9E3hkvf/8ySXLrqgfPXpU69evV2ZmptfyzMxMrVq1qtR1W7VqpaSkJHXr1k1Lly71+lpOTo7PNnv06FHqNgsLC1VQUOD1AiIdZ5IBBAJ5D8CuuHKOUGZZo75nzx4VFRUpISHBa3lCQoLy8/OLXScpKUnPP/+85s+frwULFig1NVXdunXTihUrPGPy8/P92qYkTZkyRS6Xy/NKSUk5iz0DQhuBhYrAZyxykPcA7ITmvOJwjK1l+VPfHQ6H13tjjM+yU1JTU5Wamup5n5GRoZ07d+rRRx9Vp06dzmibkjRu3DiNHDnS876goIDwRsThH1MAViLvAQQLNQ7CkWWNenx8vKKionzOfO/evdvnDHlp2rVrp1deecXzPjEx0e9tOp1OOZ3Ocn9PIJwQXggW/oxLZCDvAQQD9Y09kPXWsezW95iYGKWlpSk7O9treXZ2ttq3b1/u7eTm5iopKcnzPiMjw2ebH374oV/bBCIBt30BqAjkPYCKRH2DSGHpre8jR45UVlaW0tPTlZGRoeeff147duzQsGHDJJ28Re3nn3/WnDlzJElPPvmkGjRooGbNmuno0aN65ZVXNH/+fM2fP9+zzXvuuUedOnXS1KlTdcUVV+jtt9/WRx99pE8//dTKXQFCBuEFO+FMe2Qg7wFYidrG3sh6a1jaqPfv31979+7V5MmTlZeXp+bNm2vRokWqX7++JCkvL8/rb6wePXpUo0eP1s8//6wqVaqoWbNmev/999W7d2/PmPbt2+v111/XAw88oAcffFDnnXee5s2bp7Zt21q5K4DtEWIAgoW8B2AFahtEMkv/jrpd8XdVEU4IMYQCzrSXjWwKPI4pEHqoa0ITOV8+/uSS5U99B2ANggwAAIQL6hrAG406EIIIMwAAEOqoZ4CS0agDIYRAAwAAoY56BigbjToQAgg0hDqeCAsAoJ4Byo9GHbAxAg0AAIQ66hnAfzTqgE0RagAAIFRRxwBnh0YdsBmCDQAAhCrqGCAwaNQBmyDYAABAqKKOAQKrUrAnAIBwQ2Tgcw4A4afR04/x7ztgAa6oA0FEsAEAgFBEDQNYi0YdCBICDgAAhBrqF5yOP79qDRp1oIIRcAAAINRQvwAVi0YdqCAEHAAACDXUL0Bw0KgDFYCQAwAAoYTaBeXBbe/WoVEHLEbQAQCAUEHdAtgDjTpgEYIOAACECuoWwF5o1AELEHYAACAUULPgTHHbu7Vo1IEAIuwAAEAooGYB7I1GHQgQAg8AANgd9QoQGmjUgbNE4AEAALujXkEgcdu79WjUgbNA6AEAADujVgFCU6VgTwAIVQQfAACwM2oVWIGr6RWDK+qAnwg94MwQ7ABQMahVgNBHow74geADAAB2RZ0Cq3HSveLQqAPlRPgBAAA7okYBwg+NOlAGwg8AANgVdQoqClfTKxaNOlAKwg8AANgRNQoQ3ix/6vuMGTPUsGFDxcbGKi0tTStXrixx7IIFC3TppZeqbt26iouLU0ZGhj744AOvMbNnz5bD4fB5HTlyxOpdQYQhAAGg/Mh7oGI0evoxahRUOK6mVzxLG/V58+ZpxIgRuv/++5Wbm6uOHTuqV69e2rFjR7HjV6xYoUsvvVSLFi3S+vXr1bVrV/Xt21e5uble4+Li4pSXl+f1io2NtXJXEEEIQCDwCPjwRt4D1qM+QbCQ4cHhMMYYqzbetm1btW7dWjNnzvQsa9q0qa688kpNmTKlXNto1qyZ+vfvr4ceekjSyTPsI0aM0IEDB854XgUFBXK5XHK73YqLizvj7SD8EICANQj5soVyNpH3gLWoTxBMZHjg+JNLll1RP3r0qNavX6/MzEyv5ZmZmVq1alW5tnHixAkdPHhQtWvX9lp+6NAh1a9fX+ecc4769OnjcwYeOBOEIAD4j7wHrMNVdAQbTXrwWPYwuT179qioqEgJCQleyxMSEpSfn1+ubTz22GM6fPiw+vXr51nWpEkTzZ49Wy1atFBBQYGeeuopdejQQRs3blTjxo2L3U5hYaEKCws97wsKCs5gjxDOCEEAODPkPRB41CUALH/qu8Ph8HpvjPFZVpy5c+dq4sSJevvtt1WvXj3P8nbt2qldu3ae9x06dFDr1q31z3/+U08//XSx25oyZYomTZp0hnuAcEcYAtbibHxkIO+BwKAugV2Q38FlWaMeHx+vqKgon7Ppu3fv9jnrfrp58+bp5ptv1ptvvqnu3buXOrZSpUq6+OKLtXXr1hLHjBs3TiNHjvS8LygoUEpKSjn2AuGMIASAs0feA4FBXQI7oUkPPst+Rz0mJkZpaWnKzs72Wp6dna327duXuN7cuXM1ZMgQvfbaa7rsssvK/D7GGG3YsEFJSUkljnE6nYqLi/N6IbIRhgAQGOQ9cPaoS2AnNOn2YOmt7yNHjlRWVpbS09OVkZGh559/Xjt27NCwYcMknTzz/fPPP2vOnDmSTob2DTfcoKeeekrt2rXznJ2vUqWKXC6XJGnSpElq166dGjdurIKCAj399NPasGGDnnnmGSt3BWGEMASAwCLvgTNDTQKgJJY26v3799fevXs1efJk5eXlqXnz5lq0aJHq168vScrLy/P6G6vPPfecjh8/rjvvvFN33nmnZ/ngwYM1e/ZsSdKBAwd06623Kj8/Xy6XS61atdKKFSvUpk0bK3cFYYJABCoWZ+UjA3kP+Id6BHZFbtuHpX9H3a74u6qRiVAEKh6BX35kU+BxTGFH1COwKzLbev7kkuVPfQfsgFAEAADBRC0CwB+WPUwOsAuCEQgOzswDwEnUIrA7Mtt+uKKOsEUoAgCAYKIWQSigSbcnrqgjLBGMQHAR+gAiHbUIQgF5bV9cUUfYIRgBAECwUIcgVNCk2xtX1BFWCEcAABAs1CEAAoUr6ggbhCNgD5yhBxBpqEEQashq++OKOsICAQkAAIKBGgShhiY9NNCoI+QRkIB9EP4AIgk1CEINOR06uPUdIY2ABAAAFY36A6GIJj20cEUdIYuQBOyFAgBAJKD+QCgio0MPV9QRkghJAABQkag9EKpo0kMTV9QRcghKwH4oAgCEM2oPhCryOXTRqCOkEJQAAKAiUXsgVNGkhzZufUfIICgBe6IQABCOqDsABBNX1BESCEsAAFBRqDsQ6jiJHvpo1GF7hCVgXxQCAMINdQdCHdkcHrj1HQAAABGPBh3hgCY9fHBFHbZGaAL2RTEAIFxQbyAckMvhhUYdtkVoAvZFMQAgXFBvIByQy+GHW99hS4QmAACwErUGwgVNenjiijpsh+AE7I2CAECoo9ZAuCCTwxeNOmyF4AQAAFai1kC4oEkPbzTqsA2CE7A/igIAoYxaA+GCPA5//I46AKBcKAoAhCoadIQT8jgy0KjDFghQAABgBWoMhAsa9MjCre8IOgIUsD+KAwChiBoD4YIcjjxcUUdQEaChzZw4oSPf/6CigoOKiquh2PMayVGJ83/hhuIAQCiixggMsj74yOHIZPn/ZTNmzFDDhg0VGxurtLQ0rVy5stTxy5cvV1pammJjY9WoUSM9++yzPmPmz5+vCy64QE6nUxdccIEWLlxo1fRhIQI0tB3e+KV2Tvq78qc/q1/nvKr86c9q56S/6/DGL4M9NQBBQN7DLho9/Rg1RoCQ9cFHkx65LG3U582bpxEjRuj+++9Xbm6uOnbsqF69emnHjh3Fjt+2bZt69+6tjh07Kjc3V+PHj9fw4cM1f/58z5icnBz1799fWVlZ2rhxo7KystSvXz+tWbPGyl0B8AeHN36p3S/+W0UH3F7Liw64tfvFfxPgYYQCAeVB3sMuaNADh6wPPjI4sjmMMcaqjbdt21atW7fWzJkzPcuaNm2qK6+8UlOmTPEZP2bMGL3zzjvavHmzZ9mwYcO0ceNG5eTkSJL69++vgoICLV682DOmZ8+eqlWrlubOnVuueRUUFMjlcsntdisuLu5Mdw9ngSANXebECe2c9Hef4P6jqJo1lTJhPLfGhTgKhIoVytlE3sMOqC0Ch6wPPjI4PPmTS5b9n3X06FGtX79emZmZXsszMzO1atWqYtfJycnxGd+jRw+tW7dOx44dK3VMSduUpMLCQhUUFHi9EDwEaWg78v0PpQa3JBUdOKAj3/9QQTMCEEzkPeyA2iKwyPrgokmHZGGjvmfPHhUVFSkhIcFreUJCgvLz84tdJz8/v9jxx48f1549e0odU9I2JWnKlClyuVyeV0pKypnsEgKAIA19RQUHAzoO9kSRgPIi7xFs1BaBR9YHD/mLUyy/V8XhcHi9N8b4LCtr/OnL/d3muHHj5Ha7Pa+dO3eWe/4AvEXF1QjoONgPRQLOBHmPYKBJtwZZX/F+GD6K/IUXy/48W3x8vKKionzOfO/evdvnDPkpiYmJxY6vXLmy6tSpU+qYkrYpSU6nU06n80x2AwFEmIaH2PMaKaqmq8zfW4s9r1EFzgqBQpEAf5H3CAZqCmuR9RWL7EVxLLuiHhMTo7S0NGVnZ3stz87OVvv27YtdJyMjw2f8hx9+qPT0dEVHR5c6pqRtwh4I1PDhqFRJda6+stQxda6+gofLABGCvEdFo6awHllfcWjSURJL/+8aOXKk/vWvf+nFF1/U5s2bde+992rHjh0aNmyYpJO3qN1www2e8cOGDdOPP/6okSNHavPmzXrxxRc1a9YsjR492jPmnnvu0YcffqipU6fqm2++0dSpU/XRRx9pxIgRVu4KzgKBGn6qXdhC9W4arKiaLq/lUTVrqt5Ng1XtwhZBmhnOBsUCzhR5j4pCTVFxyHrrkbsojWW3vksn/7TK3r17NXnyZOXl5al58+ZatGiR6tevL0nKy8vz+hurDRs21KJFi3TvvffqmWeeUXJysp5++mldc801njHt27fX66+/rgceeEAPPvigzjvvPM2bN09t27a1clcAnKbahS1UtUWzk0+GLTioqLgaij2vEWfXQxTFAs4GeY+KQJNe8ch665C7KIulf0fdrvi7qhWHUAXsj2LBHsimwOOYhg/qCYQTcjdy+ZNLll5RR2QjVAH7o1gAYHfUEwgXZC78wX0rAAAAsCWadIQLmnT4i0YdliBYAfujaABgZ9QSCBfkLc4EjToARCCKBgB2RpOOcEHe4kzxO+oIOMIVsDeKBgB2Rh2BcEDW4mxxRR0BRbgC9kbhAMDOqCMQDshaBAKNOgBECAoHAHZGk45wQNYiULj1HQFDwAL2ReEAwM6oIRDqyFkEGlfUASDMUTwAsDOadIQ6chZWoFFHQBCyAADAX9QPCHU06bAKt77jrBGygH1RQACwK+oHhDLyFVbjijoAhCmKCAB2RZOOUEa+oiLQqOOsELSAPVFEALArageEMvIVFYVb3wEgzFBEALArmnSEKrIVFY0r6jhjhC1gPxQSAOyKugGhimxFMNCoA0CYoJAAYFc06QhVZCuChVvfcUYIXMBeKCQA2BU1A0IRuYpg44o6AIQ4igkAAAKHXIUdcEUdAEIYxQQAO+NqOkIJmQo74Yo6/EboAvZAQQHAzqgXEErIVNgNV9QBIARRUACwM5p0hAryFHbFFXX4heAFgo+iAoCdUSsgVJCnsDOuqANACKGoAGBnNOkIBWQpQgGNOsqN8AWCh6ICAICzR54iVNCoA4DNUVQACAWc0IedkaUINfyOOgDYGIUFgFBAkw47I0sRiriijnIhgIGKR2EBIBRQI8CuyFGEMsuuqO/fv19ZWVlyuVxyuVzKysrSgQMHShx/7NgxjRkzRi1atFC1atWUnJysG264Qb/88ovXuC5dusjhcHi9BgwYYNVuAEBQUFwgVJD3AOyIHEWos6xRHzhwoDZs2KAlS5ZoyZIl2rBhg7Kyskoc/9tvv+nzzz/Xgw8+qM8//1wLFizQt99+q8svv9xn7NChQ5WXl+d5Pffcc1btBgBUOIoLhBLyPrJxNR1288PwUeQowoIlt75v3rxZS5Ys0erVq9W2bVtJ0gsvvKCMjAxt2bJFqampPuu4XC5lZ2d7LfvnP/+pNm3aaMeOHTr33HM9y6tWrarExEQrpo5iEMJAxaG4QCgh7yMb9QHshgxFOLHkinpOTo5cLpcntCWpXbt2crlcWrVqVbm343a75XA4VLNmTa/lr776quLj49WsWTONHj1aBw8eDNTUASAouAKAUETeA7ADMhThyJIr6vn5+apXr57P8nr16ik/P79c2zhy5IjGjh2rgQMHKi4uzrN80KBBatiwoRITE7Vp0yaNGzdOGzdu9Dk7/0eFhYUqLCz0vC8oKPBjbwDAWhQXCFXkfeTiajrsgPxEOPPrivrEiRN9Huxy+mvdunWSJIfD4bO+MabY5ac7duyYBgwYoBMnTmjGjBleXxs6dKi6d++u5s2ba8CAAXrrrbf00Ucf6fPPPy9xe1OmTPE85MblciklJcWf3QYAy1BkwI7Ie5SGJh12QH4i3Pl1Rf2uu+4q84mrDRo00BdffKFdu3b5fO3XX39VQkJCqesfO3ZM/fr107Zt2/TJJ594nV0vTuvWrRUdHa2tW7eqdevWxY4ZN26cRo4c6XlfUFBAeJcTYQxYhyIDdkXeoyTUBQg2shORwq9GPT4+XvHx8WWOy8jIkNvt1meffaY2bdpIktasWSO326327duXuN6p0N66dauWLl2qOnXqlPm9vvrqKx07dkxJSUkljnE6nXI6nWVuCwAqCoUG7Iy8B2A35CYijSUPk2vatKl69uypoUOHavXq1Vq9erWGDh2qPn36eD0BtkmTJlq4cKEk6fjx47r22mu1bt06vfrqqyoqKlJ+fr7y8/N19OhRSdL333+vyZMna926ddq+fbsWLVqk6667Tq1atVKHDh2s2BUACCgeeINwQt5HFq6mI1jITUQiSx4mJ518Uuvw4cOVmZkpSbr88ss1ffp0rzFbtmyR2+2WJP3000965513JEkXXXSR17ilS5eqS5cuiomJ0ccff6ynnnpKhw4dUkpKii677DJNmDBBUVFRVu0KAAQEhQbCEXkPwCrkJiKZwxhjgj2JilZQUCCXyyW3213m78RFOs6eA4FBsYGykE2BxzENHOoBVCQyE+HKn1yy5NZ3hAdCGQgMCg4AAMqHzAROsuzWdwCIdBQbAMIBJ+5REchMwBuNOgBYgIIDQDigSYfVyEugeDTqABBgFB0AAJSOrARKR6MOAAFC0QEgnHA1HVYhL4Gy0agDQABQdAAAUDqyEig/GnUAOAsUHQAAlI6sBPxHow4AZ4jCA0C44rZ3BAI5CZw5GnUUi4AGSkfxAQBA8chI4OzRqAOAHyg+AIQ7TtbjTJGRQODQqANAOVGAAADgi3wEAo9GHcX6YfgozqgD/4cCBAAAX+QjYB0adQAoBUUIgEjCSXqUB9kIWI9GHQCKQRECAIA3shGoODTqAHAaChEAAP6LXAQqHo06APwfChEAAP6LXASCh0YdAEQxAgD8fjpOIROB4KNRBxDRKEYAADiJTATsg0YdQESiGAEA4CQyEbAfGnUAEYeCBAAA8hCwMxp1ABGDggQAAPIQCAU06ijRD8NH8WAZhAUKEgBApCMLgdBCow4grFGYAAAiGTkIhCYadQBhicIEABDJyEEgtNGoAwgrFCYAgEhGDgLhgUYdpeL31BEqKEwAAJGKDATCD406gJBGcQIAiFRkIBC+aNQBhCwKFABApCH7gMhQyaoN79+/X1lZWXK5XHK5XMrKytKBAwdKXWfIkCFyOBxer3bt2nmNKSws1N133634+HhVq1ZNl19+uX766SerdgMiEGA/PwwfxecSsAnyHqgYZB8QWSxr1AcOHKgNGzZoyZIlWrJkiTZs2KCsrKwy1+vZs6fy8vI8r0WLFnl9fcSIEVq4cKFef/11ffrppzp06JD69OmjoqIiq3YFgE1QpAD2Q96HD/59tZ9TucfPBog8ltz6vnnzZi1ZskSrV69W27ZtJUkvvPCCMjIytGXLFqWmppa4rtPpVGJiYrFfc7vdmjVrll5++WV1795dkvTKK68oJSVFH330kXr06BH4nYEkHiqH4KJAAeyJvAcCj8wDIFl0RT0nJ0cul8sT2pLUrl07uVwurVq1qtR1ly1bpnr16un888/X0KFDtXv3bs/X1q9fr2PHjikzM9OzLDk5Wc2bNy91u4WFhSooKPB6AbA/riIA9kbehx/+zQ0eMg/AH1lyRT0/P1/16tXzWV6vXj3l5+eXuF6vXr103XXXqX79+tq2bZsefPBBXXLJJVq/fr2cTqfy8/MVExOjWrVqea2XkJBQ6nanTJmiSZMmnfkOQRJX1VFxKFSA0EDeA2eHvANQEr+uqE+cONHn4S+nv9atWydJcjgcPusbY4pdfkr//v112WWXqXnz5urbt68WL16sb7/9Vu+//36p8ypru+PGjZPb7fa8du7cWc49BlCRuJoA2AN5H9n4d9ha/N45gPLw64r6XXfdpQEDBpQ6pkGDBvriiy+0a9cun6/9+uuvSkhIKPf3S0pKUv369bV161ZJUmJioo4ePar9+/d7nWXfvXu32rdvX+J2nE6nnE5nub8vSsZVdViBYgWwF/IeCCxyDoC//GrU4+PjFR8fX+a4jIwMud1uffbZZ2rTpo0kac2aNXK73aUG7On27t2rnTt3KikpSZKUlpam6OhoZWdnq1+/fpKkvLw8bdq0SdOmTfNnV3AWaNYRCBQtgH2R9yDrzx45B+BsOIwxxooN9+rVS7/88ouee+45SdKtt96q+vXr69133/WMadKkiaZMmaKrrrpKhw4d0sSJE3XNNdcoKSlJ27dv1/jx47Vjxw5t3rxZNWrUkCTdfvvteu+99zR79mzVrl1bo0eP1t69e7V+/XpFRUWVa24FBQVyuVxyu92Ki4sL/M5HAMIbZ4rCBSheqGYTeR/eyPvyI98AlMWfXLLkYXKS9Oqrr2r48OGeJ7Zefvnlmj59uteYLVu2yO12S5KioqL05Zdfas6cOTpw4ICSkpLUtWtXzZs3zxPakvTEE0+ocuXK6tevn37//Xd169ZNs2fPLndoIzA40w5/UcAA4Ym8D2/kfenINgBWseyKup1xhj1wCG+UhSIGKB+yKfA4poFF5pNpAM6OLa6oIzJwph3FoZABgPATqZlPpgEIBhp1nLVIDW74opgBgPB26t/5cM59sgyAHdCoIyBo1iMXBQ0ARJ5wadjJMAB2RaOOgKFZjywUNwCA07PAznUAuQUglNCoI6Bo1sMbRQ4AoDTBbtzJKQDhgkYdAUezHn4ofAAAZ6I8+VFWzUAGAYhENOqwRLj87lokozACAFQE8gYAfNGow1JcXQ8tFEsAAABA8NGow3JcXbc3mnMAAADAXmjUUWG4um4fNOcAAACAfdGoo0JxdT04aMwBAACA0EGjjqCgYbcezTkAAAAQmmjUEVR/bCZp2s8OjTkAAAAQHmjUYRtcZfcPjTkAAAAQnmjUYTtcZfdFUw4AAABEDhp12NrpDWokNO405QAAAEBko1FHSAm3xp2mHAAAAMDpaNQR0oprdO3UvNOIAwAAAPAXjTrCjr/NsT+NPY03AAAAAKvRqCPi0XwDAAAAsJNKwZ4AAAAAAAD4Lxp1AAAAAABshEYdAAAAAAAboVEHAAAAAMBGaNQBAAAAALARGnUAAAAAAGyERh0AAAAAABuxrFHfv3+/srKy5HK55HK5lJWVpQMHDpS6jsPhKPb1yCOPeMZ06dLF5+sDBgywajcAAEApyHsAAAKvslUbHjhwoH766SctWbJEknTrrbcqKytL7777bonr5OXleb1fvHixbr75Zl1zzTVey4cOHarJkyd73lepUiWAMwcAAOVF3gMAEHiWNOqbN2/WkiVLtHr1arVt21aS9MILLygjI0NbtmxRampqseslJiZ6vX/77bfVtWtXNWrUyGt51apVfcYCAICKRd4DAGANS259z8nJkcvl8oS2JLVr104ul0urVq0q1zZ27dql999/XzfffLPP11599VXFx8erWbNmGj16tA4ePBiwuQMAgPIh7wEAsIYlV9Tz8/NVr149n+X16tVTfn5+ubbx73//WzVq1NDVV1/ttXzQoEFq2LChEhMTtWnTJo0bN04bN25UdnZ2idsqLCxUYWGh531BQUE59wQAAJSEvAcAwBp+XVGfOHFiiQ+AOfVat26dpJMPijmdMabY5cV58cUXNWjQIMXGxnotHzp0qLp3767mzZtrwIABeuutt/TRRx/p888/L3FbU6ZM8TzkxuVyKSUlxY+9BgAgspD3AAAEl19X1O+6664yn7jaoEEDffHFF9q1a5fP13799VclJCSU+X1WrlypLVu2aN68eWWObd26taKjo7V161a1bt262DHjxo3TyJEjPe8LCgoIbwAASkDeAwAQXH416vHx8YqPjy9zXEZGhtxutz777DO1adNGkrRmzRq53W61b9++zPVnzZqltLQ0XXjhhWWO/eqrr3Ts2DElJSWVOMbpdMrpdJa5LQAAQN4DABBsljxMrmnTpurZs6eGDh2q1atXa/Xq1Ro6dKj69Onj9QTYJk2aaOHChV7rFhQU6M0339Qtt9zis93vv/9ekydP1rp167R9+3YtWrRI1113nVq1aqUOHTpYsSsAAKAE5D0AANawpFGXTj6ptUWLFsrMzFRmZqZatmypl19+2WvMli1b5Ha7vZa9/vrrMsbo+uuv99lmTEyMPv74Y/Xo0UOpqakaPny4MjMz9dFHHykqKsqqXQEAACUg7wEACDyHMcYEexIVraCgQC6XS263W3FxccGeDgAAZJMFOKYAADvxJ5csu6IOAAAAAAD8R6MOAAAAAICN0KgDAAAAAGAjNOoAAAAAANgIjToAAAAAADZCow4AAAAAgI3QqAMAAAAAYCM06gAAAAAA2AiNOgAAAAAANkKjDgAAAACAjdCoAwAAAABgIzTqAAAAAADYCI06AAAAAAA2QqMOAAAAAICN0KgDAAAAAGAjNOoAAAAAANgIjToAAAAAADZCow4AAAAAgI3QqAMAAAAAYCM06gAAAAAA2AiNOgAAAAAANkKjDgAAAACAjdCoAwAAAABgIzTqAAAAAADYCI06AAAAAAA2QqMOAAAAAICN0KgDAAAAAGAjNOoAAAAAANiIZY363//+d7Vv315Vq1ZVzZo1y7WOMUYTJ05UcnKyqlSpoi5duuirr77yGlNYWKi7775b8fHxqlatmi6//HL99NNPFuwBAAAoC3kPAEDgWdaoHz16VNddd51uv/32cq8zbdo0Pf7445o+fbrWrl2rxMREXXrppTp48KBnzIgRI7Rw4UK9/vrr+vTTT3Xo0CH16dNHRUVFVuwGAAAoBXkPAEDgOYwxxspvMHv2bI0YMUIHDhwodZwxRsnJyRoxYoTGjBkj6eTZ9ISEBE2dOlW33Xab3G636tatq5dffln9+/eXJP3yyy9KSUnRokWL1KNHj3LNqaCgQC6XS263W3FxcWe1fwAABEKoZxN5DwBA6fzJpcoVNKcybdu2Tfn5+crMzPQsczqd6ty5s1atWqXbbrtN69ev17Fjx7zGJCcnq3nz5lq1alWJwV1YWKjCwkLPe7fbLenkgQIAwA5OZZLF58+DjrwHAEQqf7LeNo16fn6+JCkhIcFreUJCgn788UfPmJiYGNWqVctnzKn1izNlyhRNmjTJZ3lKSsrZThsAgIA6ePCgXC5XsKdhGfIeABDpypP1fjXqEydOLDYA/2jt2rVKT0/3Z7NeHA6H13tjjM+y05U1Zty4cRo5cqTn/YkTJ7Rv3z7VqVOnzG2HsoKCAqWkpGjnzp3c8ldOHDP/ccz8xzHzXyQcM2OMDh48qOTk5GBPhbwPIZHw/0agccz8xzHzH8fMf5FwzPzJer8a9bvuuksDBgwodUyDBg382aRHYmKipJNn0ZOSkjzLd+/e7TnrnpiYqKNHj2r//v1eZ9l3796t9u3bl7htp9Mpp9Pptay8T6YNB3FxcWH7YbcKx8x/HDP/ccz8F+7HzC5X0sn70BPu/29YgWPmP46Z/zhm/gv3Y1berPerUY+Pj1d8fPwZTagsDRs2VGJiorKzs9WqVStJJ58ku3z5ck2dOlWSlJaWpujoaGVnZ6tfv36SpLy8PG3atEnTpk2zZF4AAEQa8h4AgOCy7HfUd+zYoX379mnHjh0qKirShg0bJEl//vOfVb16dUlSkyZNNGXKFF111VVyOBwaMWKEHn74YTVu3FiNGzfWww8/rKpVq2rgwIGSTp59uPnmmzVq1CjVqVNHtWvX1ujRo9WiRQt1797dql0BAAAlIO8BAAg8yxr1hx56SP/+978970+dNV+6dKm6dOkiSdqyZYvniaySdN999+n333/XHXfcof3796tt27b68MMPVaNGDc+YJ554QpUrV1a/fv30+++/q1u3bpo9e7aioqKs2pWQ5XQ6NWHCBJ/bAFEyjpn/OGb+45j5j2NmX+R9cPH/hv84Zv7jmPmPY+Y/jpk3y/+OOgAAAAAAKL9KwZ4AAAAAAAD4Lxp1AAAAAABshEYdAAAAAAAboVEHAAAAAMBGaNTDyN///ne1b99eVatWVc2aNcu1jjFGEydOVHJysqpUqaIuXbroq6++snaiNrJ//35lZWXJ5XLJ5XIpKytLBw4cKHWdIUOGyOFweL3atWtXMRMOkhkzZqhhw4aKjY1VWlqaVq5cWer45cuXKy0tTbGxsWrUqJGeffbZCpqpffhzzJYtW+bzmXI4HPrmm28qcMbBs2LFCvXt21fJyclyOBz6z3/+U+Y6fMYQych7/5H3ZSPr/UfWlx9Z7z8a9TBy9OhRXXfddbr99tvLvc60adP0+OOPa/r06Vq7dq0SExN16aWX6uDBgxbO1D4GDhyoDRs2aMmSJVqyZIk2bNigrKysMtfr2bOn8vLyPK9FixZVwGyDY968eRoxYoTuv/9+5ebmqmPHjurVq5d27NhR7Pht27apd+/e6tixo3JzczV+/HgNHz5c8+fPr+CZB4+/x+yULVu2eH2uGjduXEEzDq7Dhw/rwgsv1PTp08s1ns8YIh157z/yvnRkvf/Iev+Q9WfAIOy89NJLxuVylTnuxIkTJjEx0fzjH//wLDty5IhxuVzm2WeftXCG9vD1118bSWb16tWeZTk5OUaS+eabb0pcb/DgweaKK66ogBnaQ5s2bcywYcO8ljVp0sSMHTu22PH33XefadKkidey2267zbRr186yOdqNv8ds6dKlRpLZv39/BczO3iSZhQsXljqGzxhwEnlfPuR92ch6/5H1Z46sLx+uqEewbdu2KT8/X5mZmZ5lTqdTnTt31qpVq4I4s4qRk5Mjl8ultm3bepa1a9dOLperzP1ftmyZ6tWrp/PPP19Dhw7V7t27rZ5uUBw9elTr16/3+oxIUmZmZonHKCcnx2d8jx49tG7dOh07dsyyudrFmRyzU1q1aqWkpCR169ZNS5cutXKaIS3SP2OAv8h78r40ZL3/yHrrRfpnTOLW94iWn58vSUpISPBanpCQ4PlaOMvPz1e9evV8lterV6/U/e/Vq5deffVVffLJJ3rssce0du1aXXLJJSosLLRyukGxZ88eFRUV+fUZyc/PL3b88ePHtWfPHsvmahdncsySkpL0/PPPa/78+VqwYIFSU1PVrVs3rVixoiKmHHIi/TMG+Iu8J+9LQ9b7j6y3XqR/xiSpcrAngNJNnDhRkyZNKnXM2rVrlZ6efsbfw+FweL03xvgsCyXlPWaS775LZe9///79Pf/dvHlzpaenq379+nr//fd19dVXn+Gs7c3fz0hx44tbHs78OWapqalKTU31vM/IyNDOnTv16KOPqlOnTpbOM1TxGUO4Ie/9R94HFlnvP7LeWpH+GaNRt7m77rpLAwYMKHVMgwYNzmjbiYmJkk6esUpKSvIs3717t88ZrFBS3mP2xRdfaNeuXT5f+/XXX/3a/6SkJNWvX19bt271e652Fx8fr6ioKJ+zw6V9RhITE4sdX7lyZdWpU8eyudrFmRyz4rRr106vvPJKoKcXFiL9M4bwRN77j7wPDLLef2S99SL9MybRqNtefHy84uPjLdl2w4YNlZiYqOzsbLVq1UrSyd+5Wb58uaZOnWrJ96wI5T1mGRkZcrvd+uyzz9SmTRtJ0po1a+R2u9W+fftyf7+9e/dq586dXsVPuIiJiVFaWpqys7N11VVXeZZnZ2friiuuKHadjIwMvfvuu17LPvzwQ6Wnpys6OtrS+drBmRyz4uTm5oblZyoQIv0zhvBE3vuPvA8Mst5/ZL31Iv0zJomnvoeTH3/80eTm5ppJkyaZ6tWrm9zcXJObm2sOHjzoGZOammoWLFjgef+Pf/zDuFwus2DBAvPll1+a66+/3iQlJZmCgoJg7EKF69mzp2nZsqXJyckxOTk5pkWLFqZPnz5eY/54zA4ePGhGjRplVq1aZbZt22aWLl1qMjIyzJ/+9KewPWavv/66iY6ONrNmzTJff/21GTFihKlWrZrZvn27McaYsWPHmqysLM/4H374wVStWtXce++95uuvvzazZs0y0dHR5q233grWLlQ4f4/ZE088YRYuXGi+/fZbs2nTJjN27FgjycyfPz9Yu1ChDh486Pn3SpJ5/PHHTW5urvnxxx+NMXzGgNOR9/4j70tH1vuPrPcPWe8/GvUwMnjwYCPJ57V06VLPGEnmpZde8rw/ceKEmTBhgklMTDROp9N06tTJfPnllxU/+SDZu3evGTRokKlRo4apUaOGGTRokM+fzfjjMfvtt99MZmamqVu3romOjjbnnnuuGTx4sNmxY0fFT74CPfPMM6Z+/fomJibGtG7d2ixfvtzztcGDB5vOnTt7jV+2bJlp1aqViYmJMQ0aNDAzZ86s4BkHnz/HbOrUqea8884zsbGxplatWuZ//ud/zPvvvx+EWQfHqT9Zc/pr8ODBxhg+Y8DpyHv/kfdlI+v9R9aXH1nvP4cx//db+QAAAAAAIOj482wAAAAAANgIjToAAAAAADZCow4AAAAAgI3QqAMAAAAAYCM06gAAAAAA2AiNOgAAAAAANkKjDgAAAACAjdCoAwAAAABgIzTqAAAAAADYCI06AAAAAAA2QqMOAAAAAICN0KgDAAAAAGAj/x9jcN03sA+54AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1200x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[12,4])\n",
    "plt.subplot(1,2,1)\n",
    "plt.title(\"Box stencil is acute\"); plt.axis('equal')\n",
    "plt.contourf(*X,riemann2_2.norm(X),levels=[0.,1.]); plt.scatter(0,0,color='black'); \n",
    "plt.subplot(1,2,2)\n",
    "plt.title(\"Box stencil is not acute\"); plt.axis('equal')\n",
    "plt.contourf(*X,riemann2_3.norm(X),levels=[0.,1.]); plt.scatter(0,0,color='black'); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Three dimensional stencils\n",
    "\n",
    "We describe a three dimensional stencil as a set of triangles. The diamond and box stencils admit obvious generalizations, with 8 and 48 triangules respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.728700Z",
     "iopub.status.busy": "2024-04-30T08:46:15.728617Z",
     "iopub.status.idle": "2024-04-30T08:46:15.731462Z",
     "shell.execute_reply": "2024-04-30T08:46:15.731225Z"
    }
   },
   "outputs": [],
   "source": [
    "e0,e1,e2 = np.eye(3)\n",
    "diamond3 = [[e0*s0,e1*s1,e2*s2] for s0,s1,s2 in itertools.product([-1,1],repeat=3)]\n",
    "box3 = [ [u*su, u*su+v*sv, u*su+v*sv+w*sw] \n",
    "        for u,v,w in itertools.permutations([e0,e1,e2]) \n",
    "        for su,sv,sw in itertools.product([-1,1],repeat=3)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.732777Z",
     "iopub.status.busy": "2024-04-30T08:46:15.732684Z",
     "iopub.status.idle": "2024-04-30T08:46:15.734696Z",
     "shell.execute_reply": "2024-04-30T08:46:15.734468Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 48)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(diamond3),len(box3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.736019Z",
     "iopub.status.busy": "2024-04-30T08:46:15.735925Z",
     "iopub.status.idle": "2024-04-30T08:46:15.737847Z",
     "shell.execute_reply": "2024-04-30T08:46:15.737609Z"
    }
   },
   "outputs": [],
   "source": [
    "def MinStencilScal3(stencil,m):\n",
    "    def MinFacet(u,v,w):\n",
    "        return min(lp.dot_VAV(u,m,v),lp.dot_VAV(v,m,w),lp.dot_VAV(w,m,u))\n",
    "    return min(MinFacet(*s) for s in stencil)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the diamond and box stencils are acute for the isotropic metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.739149Z",
     "iopub.status.busy": "2024-04-30T08:46:15.739047Z",
     "iopub.status.idle": "2024-04-30T08:46:15.742103Z",
     "shell.execute_reply": "2024-04-30T08:46:15.741869Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1.0)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MinStencilScal3(diamond3,np.eye(3)),MinStencilScal3(box3,np.eye(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Acuteness of the diamond stencil is immediately lost for non-diagonal anisotropies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.744006Z",
     "iopub.status.busy": "2024-04-30T08:46:15.743923Z",
     "iopub.status.idle": "2024-04-30T08:46:15.745708Z",
     "shell.execute_reply": "2024-04-30T08:46:15.745446Z"
    }
   },
   "outputs": [],
   "source": [
    "riemann3_2 = Riemann.needle(direction3,1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.747014Z",
     "iopub.status.busy": "2024-04-30T08:46:15.746938Z",
     "iopub.status.idle": "2024-04-30T08:46:15.749940Z",
     "shell.execute_reply": "2024-04-30T08:46:15.749720Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.0032742967703527, 0.08245274594433738)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MinStencilScal3(diamond3,riemann3_2.m),MinStencilScal3(box3,riemann3_2.m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Acuteness for the box stencil only holds for very moderate anisotropies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.751335Z",
     "iopub.status.busy": "2024-04-30T08:46:15.751255Z",
     "iopub.status.idle": "2024-04-30T08:46:15.752981Z",
     "shell.execute_reply": "2024-04-30T08:46:15.752752Z"
    }
   },
   "outputs": [],
   "source": [
    "riemann3_3 = Riemann.needle(direction3,1,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.754344Z",
     "iopub.status.busy": "2024-04-30T08:46:15.754264Z",
     "iopub.status.idle": "2024-04-30T08:46:15.757254Z",
     "shell.execute_reply": "2024-04-30T08:46:15.757020Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-2.6753981247209406, -1.446792677481767)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MinStencilScal3(diamond3,riemann3_3.m),MinStencilScal3(box3,riemann3_3.m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Two dimensional construction\n",
    "\n",
    "We propose a stencil construction based on the preliminary computation of an $M$-obtuse superbase, see notebook [Selling](TensorSelling.ipynb).\n",
    "It can be shown that, for generic $M$, the stencil consists of the Voronoi vectors of $Z^2$ for the distance $\\|\\cdot\\|_M$, see the main introduction. \n",
    "\n",
    "The stencil is guaranteed to be acute, and to hold no more than $6$ distinct vertices, for any positive definite matrix $M$. This is independent of the condition number of $M$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.758558Z",
     "iopub.status.busy": "2024-04-30T08:46:15.758480Z",
     "iopub.status.idle": "2024-04-30T08:46:15.760252Z",
     "shell.execute_reply": "2024-04-30T08:46:15.760026Z"
    }
   },
   "outputs": [],
   "source": [
    "def MakeStencil2(m):\n",
    "    e0,e1,e2 = Selling.ObtuseSuperbase(m).T\n",
    "    return [e0,-e2,e1,-e0,e2,-e1,e0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.761538Z",
     "iopub.status.busy": "2024-04-30T08:46:15.761460Z",
     "iopub.status.idle": "2024-04-30T08:46:15.763205Z",
     "shell.execute_reply": "2024-04-30T08:46:15.762988Z"
    }
   },
   "outputs": [],
   "source": [
    "riemann2_5 = Riemann.needle(direction2,1,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.764444Z",
     "iopub.status.busy": "2024-04-30T08:46:15.764371Z",
     "iopub.status.idle": "2024-04-30T08:46:15.827074Z",
     "shell.execute_reply": "2024-04-30T08:46:15.826833Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGxCAYAAABvIsx7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAqklEQVR4nO3deXRTdeL//1foytaw1G4KZRlkF6EVWpRFkSIICi6AYAU/DK6IWBwVGRTwqxV0VBAUdRgBZQAVUBFEqgLqoSBgUUFEVLSILUuFpqh0gffvD3/NGNLdpk1uno9zco55931v3jem7ZObpTZjjBEAAICF1KntBQAAAFQ3AgcAAFgOgQMAACyHwAEAAJZD4AAAAMshcAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCB9Vi+vTpstlsOnbsWIlf79Spk/r27VulfY8dO1YtWrRwGXvsscf05ptvVmj7H374QTabTU8++WSVbr+sfS5atMg5tmjRItlsNv3www/lbr9ixQp17NhRdevWlc1m065du6ptbai6X375RSNHjlRERIRsNpuGDh1a6ty+ffvKZrM5L6GhoerQoYP+3//7fyooKHCZW9LjxV+U9P1bUyrzcwLWQ+DA602bNk2rV692GfPlH1xHjx5VcnKyWrdurfXr1ys9PV3nn39+bS8Lkh555BGtXr1aTz/9tNLT0zV79uwy57dq1Urp6elKT0/X66+/rjZt2mjatGmaMGGCy7zo6Gilp6fryiuv9OTyvVJJ3781xZd/TuCvC6ztBQDlad26dW0voVp98803Kiws1I033qg+ffpUyz5/++031atXr1r2VRsKCwtls9kUGFi7P5J2796t1q1ba/To0RWaX7duXSUkJDivDxw4UB06dNDixYs1d+5chYaGSpJCQkJc5vkTq33/wndwBge1YtOmTbLZbFq2bJmmTp2qmJgYhYWF6fLLL9e+fftc5p59ittms+nXX3/V4sWLnU8PVOTprzNnzujRRx9V8+bNFRoaqvj4eH3wwQcuc7799lvdfPPNatOmjerVq6dzzz1XQ4YM0Zdfflkdh62xY8fqkksukSSNGDHCbe1vv/22EhMTVa9ePTVs2FD9+/dXenq6yz6Knw787LPPdN1116lx48Zl/hI5evSo7rjjDnXo0EENGjRQRESELrvsMn388ccVWnOLFi00ePBgrV+/Xt26dVPdunXVrl07/ec//3Gbu3v3bl199dVq3LixQkNDdeGFF2rx4sUuc4r/37/yyiuaPHmyzj33XIWEhOjbb7/V2LFj1aBBA3399dcaMGCA6tevr+joaD3++OOSpK1bt+qSSy5R/fr1df7557vtuzS//PKL7rjjDp177rkKDg5Wq1atNHXqVOXn50v631NI77//vvbu3et8XG3atKlC+y8WGBioCy+8UAUFBTpx4oRzvLSnqPbv369Ro0YpIiJCISEhat++vebPn1/i/fXf//5X999/v6Kjo9WgQQMNGTJEhw8fVl5enm655RaFh4crPDxcN998s06ePOmyj/nz56t3796KiIhQ/fr11blzZ82ePVuFhYUu8/r27atOnTpp+/bt6tWrl+rVq6dWrVrp8ccf15kzZ9zWVJXvX0+spySl/Zz44YcfFBgYqNTUVLdtPvroI9lsNr3++uuS/ve9lpGRoWuuuUZhYWGy2+268cYbdfToUbftV6xYocTERNWvX18NGjTQgAEDlJGRUeY64UEGqAYPP/ywkWSOHj1a4tc7duxo+vTp47y+ceNGI8m0aNHCjB492qxdu9YsW7bMNG/e3LRp08YUFRU5544ZM8bExsY6r6enp5u6deuaQYMGmfT0dJOenm727NlT6toOHDhgJJlmzZqZSy65xKxcudK8/vrr5qKLLjJBQUFmy5YtzrmbN282kydPNm+88YbZvHmzWb16tRk6dKipW7eu+frrr932+fLLLzvHXn75ZSPJHDhwoNS1fPvtt2b+/PlGknnsscdc1r506VIjySQlJZk333zTrFixwsTFxZng4GDz8ccfu93XsbGx5v777zdpaWnmzTffLPU2v/76a3P77beb5cuXm02bNpl33nnHjBs3ztSpU8ds3Lix1O2KxcbGmvPOO8906NDBLFmyxLz33nvm+uuvN5LM5s2bXW6nYcOGpnXr1mbJkiVm7dq15oYbbjCSzKxZs5zziv/fn3vuuea6664zb7/9tnnnnXdMTk6OGTNmjAkODjbt27c3c+bMMWlpaebmm282ksyUKVPM+eefbxYuXGjee+89M3jwYCPJ7Nixo8z1//777+aCCy4w9evXN08++aTZsGGDmTZtmgkMDDSDBg0yxhhz6tQpk56ebrp27WpatWrlfFzl5uaWut8+ffqYjh07uo3Hx8ebRo0auTyGS3q87Nmzx9jtdtO5c2ezZMkSs2HDBjN58mRTp04dM336dLf7KzY21owdO9asX7/eLFiwwDRo0MBceumlpn///ubee+81GzZsMLNmzTIBAQHmrrvuclnTPffcY55//nmzfv168+GHH5qnn37ahIeHm5tvvtntmJo2bWratGljFixYYNLS0swdd9xhJJnFixe7rakq37+eWE9Jyvo5MWzYMNO8eXOXdRpjzPXXX29iYmJMYWGhMcb1e+0f//iHee+998xTTz1l6tevb7p27WoKCgqc2z766KPGZrOZ//u//zPvvPOOWbVqlUlMTDT169cv8+cTPIfAQbWoauAU/4Ip9tprrxlJJj093TlW0g/I+vXrmzFjxlRobcW/XGJiYszvv//uHHc4HKZJkybm8ssvL3XboqIiU1BQYNq0aWPuuecet31WNnCM+d+xv/76686x06dPm5iYGNO5c2dz+vRp53heXp6JiIgwPXv2dI4V39cPPfRQRQ6/xGMqLCw0/fr1M8OGDSt3fmxsrAkNDTU//vijc+z33383TZo0MbfeeqtzbOTIkSYkJMRkZma6bD9w4EBTr149c+LECWPM/46/d+/ebrc1ZswYI8msXLnSOVZYWGjOOeccI8l89tlnzvGcnBwTEBBgUlJSylz/ggULjCTz2muvuYzPmjXLSDIbNmxwjpUWLSUpnltYWGgKCwtNVlaWeeihh4wks2DBApe5JT1eBgwYYM477zy3iJowYYIJDQ01v/zyizHmf/fXkCFDXOZNmjTJSDITJ050GR86dKhp0qRJqes+ffq0KSwsNEuWLDEBAQHO2yk+Jklm27ZtLtt06NDBDBgwwHn9r37/Vvd6SlPaz4ni9a9evdo5dujQIRMYGGhmzJjhHCv+Xvvz974x//vHyKuvvmqMMSYzM9MEBga6hWVeXp6Jiooyw4cPL3etqH48RYVaddVVV7lcv+CCCyRJP/74Y7Xf1jXXXON8TYQkNWzYUEOGDNFHH32k06dPS5KKior02GOPqUOHDgoODlZgYKCCg4O1f/9+7d27t9rXVGzfvn36+eeflZycrDp1/vdt2aBBA1177bXaunWrfvvtN5dtrr322grvf8GCBerWrZtCQ0MVGBiooKAgffDBBxU+pgsvvFDNmzd3Xg8NDdX555/v8v/pww8/VL9+/dSsWTOXbceOHavffvvN7am20tZvs9k0aNAg5/XAwED97W9/U3R0tLp27eocb9KkiSIiIsp9rHz44YeqX7++rrvuOrd1SXJ7mrIy9uzZo6CgIAUFBSk6OlozZ87UlClTdOutt5a53alTp/TBBx9o2LBhqlevnoqKipyXQYMG6dSpU9q6davLNoMHD3a53r59e0lye+Fy+/bt9csvv7g8TZWRkaGrrrpKTZs2VUBAgIKCgnTTTTfp9OnT+uabb1y2j4qKUvfu3V3GLrjgghLv56p+/3pqPRXVt29fdenSxeXpwAULFshms+mWW25xm3/2a7KGDx+uwMBAbdy4UZL03nvvqaioSDfddJPL/8vQ0FD16dOn0k91onoQOKgWxS8OLQ6FsxUVFSkoKMhtvGnTpi7XQ0JCJEm///57Na/wjx+UJY0VFBQ4fxmkpKRo2rRpGjp0qNasWaNt27Zp+/bt6tKli0fWVCwnJ0fSH++2OVtMTIzOnDmj48ePu4yXNLckTz31lG6//Xb16NFDK1eu1NatW7V9+3ZdccUVFT6ms/8/SX/8v/rz9jk5OaWuv/jrFVl/vXr1XEJUkoKDg9WkSRO3ucHBwTp16lSZa8/JyVFUVJRsNpvLeEREhAIDA93WVRmtW7fW9u3b9emnn+r1119Xly5dlJqaquXLl5e7pqKiIj377LPOQCq+FMfd2R+5cPbxBwcHlzlefL9kZmaqV69eOnTokObMmaOPP/5Y27dvd/5yP/sxUJH/16XNrcj3ryfXUxkTJ07UBx98oH379qmwsFAvvfSSrrvuulJ/TvxZYGCgmjZt6nzsHD58WJJ00UUXuf3/XLFiRakfnwHP4l1UqBaRkZGSpEOHDjn/u5gxRllZWYqPj6+NpTllZ2eXOBYcHKwGDRpIkl599VXddNNNeuyxx1zmHTt2TI0aNfLY2op/iGdlZbl97eeff1adOnXUuHFjl/Gzf2GX5tVXX1Xfvn31/PPPu4zn5eVVcbUla9q0aanrl6Tw8HCX8YquvzrWtW3bNhljXG7zyJEjKioqcltXZRS/WF3645fbpZdeqo4dO2rSpEkaPHiw83F1tsaNGysgIEDJycm68847S5zTsmXLKq/rz9588039+uuvWrVqlWJjY53jtfXZS96ynlGjRun+++/X/PnzlZCQoOzs7FL/X2RnZ+vcc891Xi8qKlJOTo7z+7b4MfTGG2+4HBNqF2dwUC0uu+wy2Ww2rVixwu1r69evl8Ph0OWXX15tt1eVf8GtWrXK5V/7eXl5WrNmjXr16qWAgABJf/zSLf5XaLG1a9fq0KFDf33RZWjbtq3OPfdc/fe//5Uxxjn+66+/auXKlc53VlVFScf0xRdfuD1l9Ff169dPH374oTNoii1ZskT16tWrtbdJ9+vXTydPnnT7PJQlS5Y4v15dmjZtqscff1yHDx/Ws88+W+q8evXq6dJLL1VGRoYuuOACxcfHu11KOnNRFcVR9+fHgDFGL730UrXs35vXU9bPidDQUN1yyy1avHixnnrqKV144YW6+OKLS5y7dOlSl+uvvfaaioqKnO+AHDBggAIDA/Xdd9+V+P+ytv9x5684g4Nq0bp1a02YMEFPPPGETpw4oUGDBqlu3bravn27Hn/8ccXHx2vUqFHVdnudO3fWpk2btGbNGkVHR6thw4Zq27ZtmdsEBASof//+SklJ0ZkzZzRr1iw5HA7NmDHDOWfw4MFatGiR2rVrpwsuuEA7d+7UE088ofPOO6/a1l6SOnXqaPbs2Ro9erQGDx6sW2+9Vfn5+c77s/ht0lUxePBgPfLII3r44YfVp08f7du3TzNnzlTLli1VVFRUbcfw8MMP65133tGll16qhx56SE2aNNHSpUu1du1azZ49W3a7vdpuqzJuuukmzZ8/X2PGjNEPP/ygzp0765NPPtFjjz2mQYMGVWt4F9/eU089pSeffFJ33nmnwsLCSpw3Z84cXXLJJerVq5duv/12tWjRQnl5efr222+1Zs0affjhh9Wynv79+ys4OFg33HCD7rvvPp06dUrPP/+821OeNaUm11Pez4k77rhDs2fP1s6dO/Xvf/+71P2sWrVKgYGB6t+/v/bs2aNp06apS5cuGj58uKQ/Pkph5syZmjp1qr7//ntdccUVaty4sQ4fPqxPP/1U9evXd/k5g5rBGRxUmzlz5ui5557TZ599plGjRmnIkCFavHix7rzzTm3cuNH52oDquq02bdpo5MiRuuiii8p9UackTZgwQf3799fEiRM1atQoFRUVae3atS7/apszZ45uvPFGpaamasiQIXr77be1atWqGvmwslGjRunNN99UTk6ORowYoZtvvllhYWHauHGj87NzqmLq1KmaPHmyFi5cqCuvvFL//ve/tWDBgr+0z5K0bdtWW7ZsUdu2bXXnnXdq6NCh2r17t15++WX94x//qNbbqozQ0FBt3LhRo0eP1hNPPKGBAwdq0aJFuvfee7Vq1apqv706dero8ccf1y+//KJnnnmm1HkdOnTQZ599pk6dOumf//ynkpKSNG7cOL3xxhvVelapXbt2WrlypY4fP65rrrlGd911ly688ELNnTu32m7DW9dT3s+Jc889V5dccomaNGlS5j/AVq1apa+//lrXXHONHnroIQ0ZMkQbNmxw+Zk2ZcoUvfHGG/rmm280ZswYDRgwQPfdd59+/PFH9e7du9qPDeWzmT+fDwcAwE8cOXJEsbGxuuuuu0r8sxzTp0/XjBkzdPTo0b/0Wi3UDp6iAgD4lZ9++knff/+9nnjiCdWpU0d33313bS8JHsBTVAAAv/Lvf/9bffv21Z49e7R06VKXd0jBOniKCgAAWI5Hz+B89NFHGjJkiGJiYmSz2Sr0Z+s3b96suLg4hYaGqlWrVlqwYIHbnJUrV6pDhw4KCQlRhw4dtHr1ag+sHgAA+CqPBs6vv/6qLl26aN68eRWaf+DAAQ0aNEi9evVSRkaGHnzwQU2cOFErV650zklPT9eIESOUnJyszz//XMnJyRo+fLi2bdvmqcMAAAA+psaeorLZbFq9erWGDh1a6pz7779fb7/9tsvfx7ntttv0+eefOz+UbMSIEXI4HHr33Xedc4o/c2DZsmUeWz8AAPAdXvUuqvT0dCUlJbmMDRgwQAsXLlRhYaGCgoKUnp6ue+65x21OWZ83kZ+fr/z8fOf1M2fO6JdfflHTpk1r7OPiAQDAX2OMUV5enmJiYlz+MHFJvCpwsrOz3f6OUWRkpIqKinTs2DFFR0eXOqekvzNULDU1lU+RBADAIg4ePFjuJ8x7VeBI7n+Ar/gZtD+PlzSnrDMxU6ZMUUpKivN6bm6umjdvroMHD5b6MeoAAMC7OBwONWvWTA0bNix3rlcFTlRUlNuZmCNHjjj/NH1Zc84+q/NnISEhbn9sUJLCwsIIHAAAfExFXl7iVR/0l5iYqLS0NJexDRs2KD4+XkFBQWXO6dmzZ42tEwAAeDePnsE5efKkvv32W+f1AwcOaNeuXWrSpImaN2+uKVOm6NChQ1qyZImkP94xNW/ePKWkpGj8+PFKT0/XwoULXd4ddffdd6t3796aNWuWrr76ar311lt6//339cknn3jyUAAAgA/x6BmcHTt2qGvXrurataskKSUlRV27dtVDDz0kScrKylJmZqZzfsuWLbVu3Tpt2rRJF154oR555BHNnTtX1157rXNOz549tXz5cr388su64IILtGjRIq1YsUI9evTw5KEAAAAf4pd/qsHhcMhutys3N5fX4AAA4CMq8/vbq16DAwAAUB0IHAAAYDkEDgAAsBwCBwAAWA6BAwAALIfAAQAAlkPgAAAAyyFwAACA5RA4AADAcggcAABgOQQOAACwHAIHAABYDoEDAAAsh8ABAACWQ+AAAADLIXAAAIDlEDgAAMByCBwAAGA5BA4AALAcAgcAAFgOgQMAACyHwAEAAJZD4AAAAMshcAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCBwAAWA6BAwAALIfAAQAAlkPgAAAAyyFwAACA5RA4AADAcggcAABgOQQOAACwnBoJnOeee04tW7ZUaGio4uLi9PHHH5c6d+zYsbLZbG6Xjh07OucsWrSoxDmnTp2qicMBAABezuOBs2LFCk2aNElTp05VRkaGevXqpYEDByozM7PE+XPmzFFWVpbzcvDgQTVp0kTXX3+9y7ywsDCXeVlZWQoNDfX04QAAAB/g8cB56qmnNG7cOP39739X+/bt9cwzz6hZs2Z6/vnnS5xvt9sVFRXlvOzYsUPHjx/XzTff7DLPZrO5zIuKivL0oQAAAB/h0cApKCjQzp07lZSU5DKelJSkLVu2VGgfCxcu1OWXX67Y2FiX8ZMnTyo2NlbnnXeeBg8erIyMjFL3kZ+fL4fD4XIBAADW5dHAOXbsmE6fPq3IyEiX8cjISGVnZ5e7fVZWlt599139/e9/dxlv166dFi1apLffflvLli1TaGioLr74Yu3fv7/E/aSmpsputzsvzZo1q/pBAQAAr1cjLzK22Wwu140xbmMlWbRokRo1aqShQ4e6jCckJOjGG29Uly5d1KtXL7322ms6//zz9eyzz5a4nylTpig3N9d5OXjwYJWPBQAAeL9AT+48PDxcAQEBbmdrjhw54nZW52zGGP3nP/9RcnKygoODy5xbp04dXXTRRaWewQkJCVFISEjlFg8AAHyWR8/gBAcHKy4uTmlpaS7jaWlp6tmzZ5nbbt68Wd9++63GjRtX7u0YY7Rr1y5FR0f/pfUCAABr8OgZHElKSUlRcnKy4uPjlZiYqBdffFGZmZm67bbbJP3x9NGhQ4e0ZMkSl+0WLlyoHj16qFOnTm77nDFjhhISEtSmTRs5HA7NnTtXu3bt0vz58z19OAAAwAd4PHBGjBihnJwczZw5U1lZWerUqZPWrVvnfFdUVlaW22fi5ObmauXKlZozZ06J+zxx4oRuueUWZWdny263q2vXrvroo4/UvXt3Tx8OAADwATZjjKntRdQ0h8Mhu92u3NxchYWF1fZyAABABVTm9zd/iwoAAFgOgQMAACyHwAEAAJZD4AAAAMshcAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCBwAAWA6BAwAALIfAAQAAlkPgAAAAyyFwAACA5RA4AADAcggcAABgOQQOAACwHAIHAABYDoEDAAAsh8ABAACWQ+AAAADLIXAAAIDlEDgAAMByCBwAAGA5BA4AALAcAgcAAFgOgQMAACyHwAEAAJZD4AAAAMshcAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCBwAAWA6BAwAALCewthcAAICntZr7r2rd3/cTJ1fr/lD9aiRwnnvuOT3xxBPKyspSx44d9cwzz6hXr14lzt20aZMuvfRSt/G9e/eqXbt2zusrV67UtGnT9N1336l169Z69NFHNWzYMI8dAwDA+1V3yFT0dgge7+PxwFmxYoUmTZqk5557ThdffLFeeOEFDRw4UF999ZWaN29e6nb79u1TWFiY8/o555zj/O/09HSNGDFCjzzyiIYNG6bVq1dr+PDh+uSTT9SjRw+PHg8AwDvUVMxURPFaCB3vYTPGGE/eQI8ePdStWzc9//zzzrH27dtr6NChSk1NdZtffAbn+PHjatSoUYn7HDFihBwOh959913n2BVXXKHGjRtr2bJl5a7J4XDIbrcrNzfXJaIAAN7Lm4KmNASOZ1Xm97dHX2RcUFCgnTt3KikpyWU8KSlJW7ZsKXPbrl27Kjo6Wv369dPGjRtdvpaenu62zwEDBpS6z/z8fDkcDpcLAMD7tZr7L+fFF/jKOv2BR5+iOnbsmE6fPq3IyEiX8cjISGVnZ5e4TXR0tF588UXFxcUpPz9fr7zyivr166dNmzapd+/ekqTs7OxK7TM1NVUzZsyohiMCAHgSgYDqUiMvMrbZbC7XjTFuY8Xatm2rtm3bOq8nJibq4MGDevLJJ52BU9l9TpkyRSkpKc7rDodDzZo1q/RxAACqH1EDT/Bo4ISHhysgIMDtzMqRI0fczsCUJSEhQa+++qrzelRUVKX2GRISopCQkEqsHADgSVaOmlZz/8VrcbyAR1+DExwcrLi4OKWlpbmMp6WlqWfPnhXeT0ZGhqKjo53XExMT3fa5YcOGSu0TAFCzfO31NPBtHn+KKiUlRcnJyYqPj1diYqJefPFFZWZm6rbbbpP0x9NHhw4d0pIlSyRJzzzzjFq0aKGOHTuqoKBAr776qlauXKmVK1c693n33Xerd+/emjVrlq6++mq99dZbev/99/XJJ594+nAAAJVAzKC2eDxwRowYoZycHM2cOVNZWVnq1KmT1q1bp9jYWElSVlaWMjMznfMLCgp077336tChQ6pbt646duyotWvXatCgQc45PXv21PLly/XPf/5T06ZNU+vWrbVixQo+AwcAvABRw9NU3sDjn4PjjfgcHACoXkSNOwKn+lXm9zd/iwoAUCVETemIm9pH4AAAKoWwgS8gcAAA5SJq4GsIHABAiYiaquHpKe9A4AAAXBA2sAICBwBA1MByCBwA8GOETfXi6SnvQeAAgJ8hauAPCBwA8BOEjWdx9sa7EDgAYHGEDfwRgQMAFkTU1CzO3ngfAgcALISwqXnEjXcicADAAggbwBWBAwA+iqipfZy98V4EDgD4GMIGKB+BAwA+grDxLpy98W4EDgB4OcLG+xA33o/AAQAvRdh4J+LGNxA4AOBFiBrvRtz4DgIHALwAYQNULwIHAGoRYeM7OHvjWwgcAKgFhI1vIW58D4EDADWIsPE9xI1vInAAoAYQNr6JuPFdBA4AeBBh47uIG99G4ACABxA2vo248X0EDgBUI8LG9xE31kDgAEA1IGwA70LgAMBfQNhYC2dvrIPAAYAqIGysh7ixFgIHACqBsLEm4sZ6CBwAqADCxrqIG2sicACgDISNtRE31kXgAEAJCBvrI26sjcABgD8hbPwDcWN9dWriRp577jm1bNlSoaGhiouL08cff1zq3FWrVql///4655xzFBYWpsTERL333nsucxYtWiSbzeZ2OXXqlKcPBYBFtZr7L+LGTxA3/sHjgbNixQpNmjRJU6dOVUZGhnr16qWBAwcqMzOzxPkfffSR+vfvr3Xr1mnnzp269NJLNWTIEGVkZLjMCwsLU1ZWlsslNDTU04cDwGIIG/9C3PgPmzHGePIGevTooW7duun55593jrVv315Dhw5VampqhfbRsWNHjRgxQg899JCkP87gTJo0SSdOnKjSmhwOh+x2u3JzcxUWFlalfQDwbUSNfyFsrKEyv789eganoKBAO3fuVFJSkst4UlKStmzZUqF9nDlzRnl5eWrSpInL+MmTJxUbG6vzzjtPgwcPdjvD82f5+flyOBwuFwD+iTM2/oe48U8eDZxjx47p9OnTioyMdBmPjIxUdnZ2hfbxr3/9S7/++quGDx/uHGvXrp0WLVqkt99+W8uWLVNoaKguvvhi7d+/v8R9pKamym63Oy/NmjWr+kEB8FmEjf8hbvxXjbyLymazuVw3xriNlWTZsmWaPn263nrrLUVERDjHExISlJCQ4Lx+8cUXq1u3bnr22Wc1d+5ct/1MmTJFKSkpzusOh4PIAfwIYeOfiBv/5tHACQ8PV0BAgNvZmiNHjrid1TnbihUrNG7cOL3++uu6/PLLy5xbp04dXXTRRaWewQkJCVFISEjlFg/A5xE2/ou4gUefogoODlZcXJzS0tJcxtPS0tSzZ89St1u2bJnGjh2r//73v7ryyivLvR1jjHbt2qXo6Oi/vGYAvo/X2fg34gZSDTxFlZKSouTkZMXHxysxMVEvvviiMjMzddttt0n64+mjQ4cOacmSJZL+iJubbrpJc+bMUUJCgvPsT926dWW32yVJM2bMUEJCgtq0aSOHw6G5c+dq165dmj9/vqcPB4AXI2pA3KCYxwNnxIgRysnJ0cyZM5WVlaVOnTpp3bp1io2NlSRlZWW5fCbOCy+8oKKiIt1555268847neNjxozRokWLJEknTpzQLbfcouzsbNntdnXt2lUfffSRunfv7unDAeCFCBtIxA1cefxzcLwRn4MDWAdxA8LGf1Tm9zd/iwqATyJsIBE3KB2BA8CnEDYoRtygLAQOAJ9A2ODPiBuUh8AB4NUIG5yNuEFFePyviQNAVRE3OBtxg4riDA4Ar0PY4GyEDSqLwAHgNQgblIS4QVUQOABqHWGD0hA3qCpegwOgVhE3KA1xg7+CMzgAagVhg9IQNqgOBA6AGkXYoCzEDaoLT1EBqDHEDcpC3KA6cQYHgMcRNigPcYPqRuAA8BjCBuUhbOApPEUFwCOIG5SHuIEncQYHQLUibFARxA08jcABUC0IG1QEYYOawlNUAP4y4gYVQdygJnEGB0CVETaoKOIGNY3AAVBphA0qirBBbeEpKgCVQtygoogb1CbO4ACoEMIGlUHcoLYROADKRNigMggbeAueogJQKuIGlUHcwJtwBgeAG8IGlUXcwNsQOABcEDeoDMIG3orAASCJsEHlETfwZrwGBwBxg0ojbuDtOIMD+DHCBpVF2MBXEDiAHyJsUBXEDXwJgQP4GeIGlUXYwBcROICfIGxQFcQNfBUvMgb8AHGDqiBu4Ms4gwNYGGGDqiBsYAUEDmBR/hg35swZnfrue5125CkgrKFCW7eSrQ4nqiuDuIFVEDiAxfhj2EjSr59/qZxVb+r0iVznWEAju5peM1T1u3SuxZX5BsIGVlMj/7R57rnn1LJlS4WGhiouLk4ff/xxmfM3b96suLg4hYaGqlWrVlqwYIHbnJUrV6pDhw4KCQlRhw4dtHr1ak8tH/AZ/hw3R/6z2CVuJOn0iVwd+c9i/fr5l7W0Mt9A3MCKPB44K1as0KRJkzR16lRlZGSoV69eGjhwoDIzM0ucf+DAAQ0aNEi9evVSRkaGHnzwQU2cOFErV650zklPT9eIESOUnJyszz//XMnJyRo+fLi2bdvm6cMBvFKruf/y27gxZ84oZ9WbZc7JWfWWzJkzNbMgH/L9xMnEDSzLZowxnryBHj16qFu3bnr++eedY+3bt9fQoUOVmprqNv/+++/X22+/rb179zrHbrvtNn3++edKT0+XJI0YMUIOh0Pvvvuuc84VV1yhxo0ba9myZW77zM/PV35+vvO6w+FQs2bNlJubq7CwsGo5TqC2+GvYFPt9/7fKnud+lvdsURNuU902f6uBFfkGwga+yOFwyG63V+j3t0fP4BQUFGjnzp1KSkpyGU9KStKWLVtK3CY9Pd1t/oABA7Rjxw4VFhaWOae0faampsputzsvzZo1q+ohAV7Dn8/a/NlpR161zrM6ztrAX3g0cI4dO6bTp08rMjLSZTwyMlLZ2dklbpOdnV3i/KKiIh07dqzMOaXtc8qUKcrNzXVeDh48WNVDArwCYfM/AWENq3WelRE28Cc18i4qm83mct0Y4zZW3vyzxyuzz5CQEIWEhFRqzYA3ImzchbZupYBGdrcXGP9ZQKNGCm3dqgZX5V0IG/gjj57BCQ8PV0BAgNuZlSNHjridgSkWFRVV4vzAwEA1bdq0zDml7ROwAuKmZLY6ddT0mqFlzml6zdV++3k4xA38lUe/44ODgxUXF6e0tDSX8bS0NPXs2bPEbRITE93mb9iwQfHx8QoKCipzTmn7BHwZr7UpX/0unRXxf2MU0MjuMh7QqJEi/m+MX34ODq+1gb/z+FNUKSkpSk5OVnx8vBITE/Xiiy8qMzNTt912m6Q/Xh9z6NAhLVmyRNIf75iaN2+eUlJSNH78eKWnp2vhwoUu7466++671bt3b82aNUtXX3213nrrLb3//vv65JNPPH04QI0ibCqufpfOqte5I59kLM7aAFINBM6IESOUk5OjmTNnKisrS506ddK6desUGxsrScrKynL5TJyWLVtq3bp1uueeezR//nzFxMRo7ty5uvbaa51zevbsqeXLl+uf//ynpk2bptatW2vFihXq0aOHpw8HqDHETeXZ6tTx67eCEzbA/3j8c3C8UWXeRw/UNMIGVUHcwB9U5vc3f4sK8CLEDSqLsAFKRuAAXoCwQVUQN0DpCByglhE3qCzCBigfgQPUEsIGVUHcABVD4AC1gLhBZRE2QOX43wdEALWMuEFlETdA5XEGB6ghhA0qi7ABqo7AAWoAcYPKIGyAv46nqAAPI25QGcQNUD04gwN4CGGDyiBsgOrFGRzAA4gbVAZxA1Q/zuAA1YiwQWUQNoDnEDhANSFuUFGEDeB5PEUFVAPiBhVF3AA1gzM4wF9E3KAiCBugZhE4QBURNqgIwgaoHTxFBVQBcYOKIG6A2sMZHKCSiBuUh7ABah+BA1QQYYPyEDaA9+ApKqACiBuUh7gBvAtncIByEDcoC2EDeCcCBygDcYPSEDaAdyNwgBIQNigLcQN4PwIHOAtxg9IQNoDvIHCAPyFuUBLCBvA9BA4gwgYlI2wA38XbxOH3iBuUhLgBfBtncODXiBucjbABrIHAgd8ibvBnhA1gLQQO/BJxg2KEDWBNvAYHfoe4QTHiBrAuzuDAbxA2KEbYANZH4MAvEDeQCBvAnxA4sDziBoQN4H8IHFgacePfCBvAf3n0RcbHjx9XcnKy7Ha77Ha7kpOTdeLEiVLnFxYW6v7771fnzp1Vv359xcTE6KabbtLPP//sMq9v376y2Wwul5EjR3ryUOCDiBv/RtwA/s2jZ3BGjRqln376SevXr5ck3XLLLUpOTtaaNWtKnP/bb7/ps88+07Rp09SlSxcdP35ckyZN0lVXXaUdO3a4zB0/frxmzpzpvF63bl3PHQh8DnHjvwgbAJIHA2fv3r1av369tm7dqh49ekiSXnrpJSUmJmrfvn1q27at2zZ2u11paWkuY88++6y6d++uzMxMNW/e3Dler149RUVFeWr58GHEjX8ibAD8mccCJz09XXa73Rk3kpSQkCC73a4tW7aUGDglyc3Nlc1mU6NGjVzGly5dqldffVWRkZEaOHCgHn74YTVs2LDEfeTn5ys/P9953eFwVP6A4PUIG/9E2AAoiccCJzs7WxEREW7jERERys7OrtA+Tp06pQceeECjRo1SWFiYc3z06NFq2bKloqKitHv3bk2ZMkWff/6529mfYqmpqZoxY0bVDgQ+gbjxP4QNgLJUOnCmT59ebixs375dkmSz2dy+ZowpcfxshYWFGjlypM6cOaPnnnvO5Wvjx493/nenTp3Upk0bxcfH67PPPlO3bt3c9jVlyhSlpKQ4rzscDjVr1qzcNcA3EDf+hbABUBGVDpwJEyaU+46lFi1a6IsvvtDhw4fdvnb06FFFRkaWuX1hYaGGDx+uAwcO6MMPP3Q5e1OSbt26KSgoSPv37y8xcEJCQhQSElLmPuCbiBv/QdgAqIxKB054eLjCw8PLnZeYmKjc3Fx9+umn6t69uyRp27Ztys3NVc+ePUvdrjhu9u/fr40bN6pp06bl3taePXtUWFio6Ojoih8IfB5x4z+IGwCVZTPGGE/tfODAgfr555/1wgsvSPrjbeKxsbEubxNv166dUlNTNWzYMBUVFenaa6/VZ599pnfeecflTE+TJk0UHBys7777TkuXLtWgQYMUHh6ur776SpMnT1bdunW1fft2BQQElLsuh8Mhu92u3Nzccs8OwTsRN/6BsAHwZ5X5/e3Rz8FZunSpJk6cqKSkJEnSVVddpXnz5rnM2bdvn3JzcyVJP/30k95++21J0oUXXugyb+PGjerbt6+Cg4P1wQcfaM6cOTp58qSaNWumK6+8Ug8//HCF4ga+j7ixPsIGwF/l0TM43oozOL6LuLE2wgZAWbzmDA5QnYgb6yJsAFQ3Agc+gbixJsIGgKcQOPB6xI31EDYAPI3AgVcjbqyFsAFQUwgceC3ixjoIGwA1jcCBVyJurIGwAVBbCBx4HeLG9xE2AGobgQOvQtz4NsIGgLcgcOA1iBvfRdgA8DYEDoAqI2wAeCsCB16Bsze+hbAB4O0IHNQ64sY3EDUAfAmBg1pF3Hg/wgaALyJwUGuIG+9G2ADwZQQOagVx470IGwBWQOAAkETYALAWAgc1jrM33oWwAWBFBA5qFHHjHYgaAFZH4KDGEDe1j7AB4C8IHMAPEDYA/A2BgxrB2ZvaQdgA8FcEDmAxRA0AEDioAZy9qRmEDQD8D4EDjyJuPI+wAQB3BA7gg4gaACgbgQOP4exN9SNsAKBiCBzAyxE1AFB5BA7gpQgbAKg6AgcewdNTVUPUAED1IHAAL0DYAED1InCAWkLUAIDnEDiodjw9VTqiBgBqBoED1ADCBgBqFoEDeAhRAwC1h8ABqhFRAwDeoY4nd378+HElJyfLbrfLbrcrOTlZJ06cKHObsWPHymazuVwSEhJc5uTn5+uuu+5SeHi46tevr6uuuko//fSTB48EleFvv+S/nzjZeQEAeAePnsEZNWqUfvrpJ61fv16SdMsttyg5OVlr1qwpc7srrrhCL7/8svN6cHCwy9cnTZqkNWvWaPny5WratKkmT56swYMHa+fOnQoICKj+AwHOQswAgHfzWODs3btX69ev19atW9WjRw9J0ksvvaTExETt27dPbdu2LXXbkJAQRUVFlfi13NxcLVy4UK+88oouv/xySdKrr76qZs2a6f3339eAAQOq/2Dg9wgaAPAtHnuKKj09XXa73Rk3kpSQkCC73a4tW7aUue2mTZsUERGh888/X+PHj9eRI0ecX9u5c6cKCwuVlJTkHIuJiVGnTp1K3W9+fr4cDofLBZ5lhSDgqScA8F0eO4OTnZ2tiIgIt/GIiAhlZ2eXut3AgQN1/fXXKzY2VgcOHNC0adN02WWXaefOnQoJCVF2draCg4PVuHFjl+0iIyNL3W9qaqpmzJjx1w4IfoGYAQBrqHTgTJ8+vdxY2L59uyTJZrO5fc0YU+J4sREjRjj/u1OnToqPj1dsbKzWrl2ra665ptTtytrvlClTlJKS4rzucDjUrFmzMo8B/oGgAQBrqnTgTJgwQSNHjixzTosWLfTFF1/o8OHDbl87evSoIiMjK3x70dHRio2N1f79+yVJUVFRKigo0PHjx13O4hw5ckQ9e/YscR8hISEKCQmp8G2ienw/cbLXfaoxQQMA/qHSgRMeHq7w8PBy5yUmJio3N1effvqpunfvLknatm2bcnNzSw2RkuTk5OjgwYOKjo6WJMXFxSkoKEhpaWkaPny4JCkrK0u7d+/W7NmzK3s48LDioKiN0CFmAMB/2YwxxlM7HzhwoH7++We98MILkv54m3hsbKzL28TbtWun1NRUDRs2TCdPntT06dN17bXXKjo6Wj/88IMefPBBZWZmau/evWrYsKEk6fbbb9c777yjRYsWqUmTJrr33nuVk5NT4beJOxwO2e125ebmKiwszDMHjzJVZ/AQMgDgHyrz+9ujn4OzdOlSTZw40fmOp6uuukrz5s1zmbNv3z7l5uZKkgICAvTll19qyZIlOnHihKKjo3XppZdqxYoVzriRpKefflqBgYEaPny4fv/9d/Xr10+LFi3iM3B8CFECAPAkj57B8VacwQEAwPdU5ve3R/9UAwAAQG0gcAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCBwAAWA6BAwAALIfAAQAAlkPgAAAAyyFwAACA5RA4AADAcggcAABgOQQOAACwHAIHAABYDoEDAAAsh8ABAACWQ+AAAADLIXAAAIDlEDgAAMByCBwAAGA5BA4AALAcAgcAAFgOgQMAACyHwAEAAJZD4AAAAMshcAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCBwAAWA6BAwAALIfAAQAAlkPgAAAAyyFwAACA5RA4AADAcjwaOMePH1dycrLsdrvsdruSk5N14sSJMrex2WwlXp544gnnnL59+7p9feTIkZ48FAAA4EMCPbnzUaNG6aefftL69eslSbfccouSk5O1Zs2aUrfJyspyuf7uu+9q3Lhxuvbaa13Gx48fr5kzZzqv161btxpXDgAAfJnHAmfv3r1av369tm7dqh49ekiSXnrpJSUmJmrfvn1q27ZtidtFRUW5XH/rrbd06aWXqlWrVi7j9erVc5sLAAAgefApqvT0dNntdmfcSFJCQoLsdru2bNlSoX0cPnxYa9eu1bhx49y+tnTpUoWHh6tjx4669957lZeXV+p+8vPz5XA4XC4AAMC6PHYGJzs7WxEREW7jERERys7OrtA+Fi9erIYNG+qaa65xGR89erRatmypqKgo7d69W1OmTNHnn3+utLS0EveTmpqqGTNmVP4gAACAT6r0GZzp06eX+kLg4suOHTsk/fGC4bMZY0ocL8l//vMfjR49WqGhoS7j48eP1+WXX65OnTpp5MiReuONN/T+++/rs88+K3E/U6ZMUW5urvNy8ODBSh41AADwJZU+gzNhwoRy37HUokULffHFFzp8+LDb144eParIyMhyb+fjjz/Wvn37tGLFinLnduvWTUFBQdq/f7+6devm9vWQkBCFhISUux8AAGANlQ6c8PBwhYeHlzsvMTFRubm5+vTTT9W9e3dJ0rZt25Sbm6uePXuWu/3ChQsVFxenLl26lDt3z549KiwsVHR0dPkHAAAALM9jLzJu3769rrjiCo0fP15bt27V1q1bNX78eA0ePNjlHVTt2rXT6tWrXbZ1OBx6/fXX9fe//91tv999951mzpypHTt26IcfftC6det0/fXXq2vXrrr44os9dTgAAMCHePSD/pYuXarOnTsrKSlJSUlJuuCCC/TKK6+4zNm3b59yc3NdxpYvXy5jjG644Qa3fQYHB+uDDz7QgAED1LZtW02cOFFJSUl6//33FRAQ4MnDAQAAPsJmjDG1vYia5nA4ZLfblZubq7CwsNpeDgAAqIDK/P7mb1EBAADLIXAAAIDlEDgAAMByCBwAAGA5BA4AALAcAgcAAFgOgQMAACyHwAEAAJZD4AAAAMshcAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCBwAAWA6BAwAALIfAAQAAlkPgAAAAyyFwAACA5RA4AADAcggcAABgOQQOAACwHAIHAABYDoEDAAAsh8ABAACWQ+AAAADLIXAAAIDlEDgAAMByCBwAAGA5BA4AALAcAgcAAFgOgQMAACyHwAEAAJZD4AAAAMshcAAAgOUQOAAAwHI8GjiPPvqoevbsqXr16qlRo0YV2sYYo+nTpysmJkZ169ZV3759tWfPHpc5+fn5uuuuuxQeHq769evrqquu0k8//eSBIwAAAL7Io4FTUFCg66+/XrfffnuFt5k9e7aeeuopzZs3T9u3b1dUVJT69++vvLw855xJkyZp9erVWr58uT755BOdPHlSgwcP1unTpz1xGAAAwMfYjDHG0zeyaNEiTZo0SSdOnChznjFGMTExmjRpku6//35Jf5ytiYyM1KxZs3TrrbcqNzdX55xzjl555RWNGDFCkvTzzz+rWbNmWrdunQYMGFDuehwOh+x2u3JzcxUWFvaXjw8AAHheZX5/B9bQmirkwIEDys7OVlJSknMsJCREffr00ZYtW3Trrbdq586dKiwsdJkTExOjTp06acuWLSUGTn5+vvLz853Xc3NzJf1xRwEAAN9Q/Hu7IudmvCpwsrOzJUmRkZEu45GRkfrxxx+dc4KDg9W4cWO3OcXbny01NVUzZsxwG2/WrFl1LBsAANSgvLw82e32MudUOnCmT59eYiz82fbt2xUfH1/ZXTvZbDaX68YYt7GzlTVnypQpSklJcV4/c+aMfvnlFzVt2rTc/VqRw+FQs2bNdPDgQZ6iqwHc3zWL+7vmcZ/XLH++v40xysvLU0xMTLlzKx04EyZM0MiRI8uc06JFi8ruVpIUFRUl6Y+zNNHR0c7xI0eOOM/qREVFqaCgQMePH3c5i3PkyBH17NmzxP2GhIQoJCTEZayi7+qysrCwML/75qhN3N81i/u75nGf1yx/vb/LO3NTrNKBEx4ervDw8EovqCJatmypqKgopaWlqWvXrpL+eCfW5s2bNWvWLElSXFycgoKClJaWpuHDh0uSsrKytHv3bs2ePdsj6wIAAL7Fo6/ByczM1C+//KLMzEydPn1au3btkiT97W9/U4MGDSRJ7dq1U2pqqoYNGyabzaZJkybpscceU5s2bdSmTRs99thjqlevnkaNGiXpj3IbN26cJk+erKZNm6pJkya699571blzZ11++eWePBwAAOAjPBo4Dz30kBYvXuy8XnxWZuPGjerbt68kad++fc53NUnSfffdp99//1133HGHjh8/rh49emjDhg1q2LChc87TTz+twMBADR8+XL///rv69eunRYsWKSAgwJOHYxkhISF6+OGH3Z62g2dwf9cs7u+ax31es7i/K6ZGPgcHAACgJvG3qAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCBwAAWA6B4wceffRR9ezZU/Xq1avwJzgbYzR9+nTFxMSobt266tu3r/bs2ePZhVrI8ePHlZycLLvdLrvdruTkZJ04caLMbcaOHSubzeZySUhIqJkF+5jnnntOLVu2VGhoqOLi4vTxxx+XOX/z5s2Ki4tTaGioWrVqpQULFtTQSq2jMvf5pk2b3B7LNptNX3/9dQ2u2Hd99NFHGjJkiGJiYmSz2fTmm2+Wuw2PcXcEjh8oKCjQ9ddfr9tvv73C28yePVtPPfWU5s2bp+3btysqKkr9+/dXXl6eB1dqHaNGjdKuXbu0fv16rV+/Xrt27VJycnK5211xxRXKyspyXtatW1cDq/UtK1as0KRJkzR16lRlZGSoV69eGjhwoDIzM0ucf+DAAQ0aNEi9evVSRkaGHnzwQU2cOFErV66s4ZX7rsre58X27dvn8nhu06ZNDa3Yt/3666/q0qWL5s2bV6H5PMZLYeA3Xn75ZWO328udd+bMGRMVFWUef/xx59ipU6eM3W43CxYs8OAKreGrr74ykszWrVudY+np6UaS+frrr0vdbsyYMebqq6+ugRX6tu7du5vbbrvNZaxdu3bmgQceKHH+fffdZ9q1a+cyduutt5qEhASPrdFqKnufb9y40Ugyx48fr4HVWZsks3r16jLn8BgvGWdw4ObAgQPKzs5WUlKScywkJER9+vTRli1banFlviE9PV12u109evRwjiUkJMhut5d7/23atEkRERE6//zzNX78eB05csTTy/UpBQUF2rlzp8tjU5KSkpJKvW/T09Pd5g8YMEA7duxQYWGhx9ZqFVW5z4t17dpV0dHR6tevnzZu3OjJZfo1HuMlI3DgJjs7W5Kcf8G9WGRkpPNrKF12drYiIiLcxiMiIsq8/wYOHKilS5fqww8/1L/+9S9t375dl112mfLz8z25XJ9y7NgxnT59ulKPzezs7BLnFxUV6dixYx5bq1VU5T6Pjo7Wiy++qJUrV2rVqlVq27at+vXrp48++qgmlux3eIyXzKN/iwqeM336dM2YMaPMOdu3b1d8fHyVb8Nms7lcN8a4jfmTit7nkvt9J5V//40YMcL53506dVJ8fLxiY2O1du1aXXPNNVVctTVV9rFZ0vySxlG6ytznbdu2Vdu2bZ3XExMTdfDgQT355JPq3bu3R9fpr3iMuyNwfNSECRM0cuTIMue0aNGiSvuOioqS9Me/CqKjo53jR44ccftXgj+p6H3+xRdf6PDhw25fO3r0aKXuv+joaMXGxmr//v2VXqtVhYeHKyAgwO3MQVmPzaioqBLnBwYGqmnTph5bq1VU5T4vSUJCgl599dXqXh7EY7w0BI6PCg8PV3h4uEf23bJlS0VFRSktLc35F+ALCgq0efNmzZo1yyO36Qsqep8nJiYqNzdXn376qbp37y5J2rZtm3Jzc9WzZ88K315OTo4OHjzoEpn+Ljg4WHFxcUpLS9OwYcOc42lpabr66qtL3CYxMVFr1qxxGduwYYPi4+MVFBTk0fVaQVXu85JkZGTwWPYQHuOlqM1XOKNm/PjjjyYjI8PMmDHDNGjQwGRkZJiMjAyTl5fnnNO2bVuzatUq5/XHH3/c2O12s2rVKvPll1+aG264wURHRxuHw1Ebh+BzrrjiCnPBBReY9PR0k56ebjp37mwGDx7sMufP93leXp6ZPHmy2bJlizlw4IDZuHGjSUxMNOeeey73+VmWL19ugoKCzMKFC81XX31lJk2aZOrXr29++OEHY4wxDzzwgElOTnbO//777029evXMPffcY7766iuzcOFCExQUZN54443aOgSfU9n7/OmnnzarV68233zzjdm9e7d54IEHjCSzcuXK2joEn5KXl+f8OS3JPPXUUyYjI8P8+OOPxhge4xVF4PiBMWPGGElul40bNzrnSDIvv/yy8/qZM2fMww8/bKKiokxISIjp3bu3+fLLL2t+8T4qJyfHjB492jRs2NA0bNjQjB492u0ts3++z3/77TeTlJRkzjnnHBMUFGSaN29uxowZYzIzM2t+8T5g/vz5JjY21gQHB5tu3bqZzZs3O782ZswY06dPH5f5mzZtMl27djXBwcGmRYsW5vnnn6/hFfu+ytzns2bNMq1btzahoaGmcePG5pJLLjFr166thVX7puK32Z99GTNmjDGGx3hF2Yz5/1+JBAAAYBG8TRwAAFgOgQMAACyHwAEAAJZD4AAAAMshcAAAgOUQOAAAwHIIHAAAYDkEDgAAsBwCBwAAWA6BAwAALIfAAQAAlvP/AY9jloe3b5FMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Unit ball for a norm of Riemannian type\"); plt.axis('equal')\n",
    "plt.contourf(*X,riemann2_5.norm(X),levels=[0.,1.]); plt.scatter(0,0,color='black'); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.828547Z",
     "iopub.status.busy": "2024-04-30T08:46:15.828469Z",
     "iopub.status.idle": "2024-04-30T08:46:15.895842Z",
     "shell.execute_reply": "2024-04-30T08:46:15.895568Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGxCAYAAACqUFbqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9lklEQVR4nO3deXhTdaL/8U8obVq0BKG0DVsLI6sgsggURTZF2dz3EcFRFBxUQK5adIQ6apVRxFER8CKIuDAKOKMgwh0o4KUKaHtlEVwGBKEVQWgBpYX2+/uDXyOhaZuUbCd9v54nz2NOvif5np6meXOSHG3GGCMAAACLqBXqCQAAAPiCeAEAAJZCvAAAAEshXgAAgKUQLwAAwFKIFwAAYCnECwAAsBTiBQAAWArxAgAALIV4wRn7/PPPdc0116hZs2ay2+1KSkpSWlqaHnzwQbdx06dP19y5c0MzSQ+ysrJks9mUlZXlWjZ58mTZbDav1n/ppZd07rnnKiYmRjabTYcOHQrMROGTnTt3avDgwapfv75sNpvGjh1b4djU1FTZbDbX5ayzzlLnzp318ssv6/STj3v6fakp+vTpoz59+gT9cX/99VdNnjy5Rv7MUbnaoZ4ArG3JkiW68sor1adPH02ZMkVOp1N5eXnauHGj3n33XT3//POusdOnT1dCQoJGjBgRugmfonPnzsrOzla7du18Xjc3N1f333+/7rrrLg0fPly1a9dWfHx8AGYJX40bN06ff/65Xn/9dSUnJ8vpdFY6/qKLLtJzzz0nSdq7d6+mTp2q++67T4WFhZo4caJr3Jn8vljd9OnTQ/K4v/76qzIyMiQpJPGE8EW84IxMmTJFzZs31yeffKLatX//dbr55ps1ZcqUEM6sanXr1lWPHj2qte6WLVskSSNHjlS3bt38Mp9ff/1VderU8ct9hcJvv/2m2NhYr49cBcrmzZvVrVs3XX311V6Nr1evntvvwaWXXqpmzZpp5syZbvFyJr8vVlcTgw3hjbeNcEYOHDighIQEt3ApU6vW779eqamp2rJli1avXu06RJ+amuq6vbCwUBMmTFDz5s0VExOjxo0ba+zYsTp69KjbfdpsNo0ZM0Zvvvmm2rZtqzp16qhjx4766KOPyj3+tm3bdMsttygpKUl2u13NmjXT7bffrqKiIknVfxugT58+uu222yRJ3bt3l81mczua9Prrr6tjx46KjY1V/fr1dc011+jrr792u48RI0bo7LPP1qZNmzRgwADFx8erf//+FT7md999pzvuuEMtW7ZUnTp11LhxYw0dOlSbNm3yas6+/Nw+/fRT9e/fX/Hx8apTp4569uypJUuWuI2ZO3eubDabli9frj/96U9q2LCh6tSpo6KiIvXp00ft27dXdna2evbsqbi4OKWmpmrOnDmSTh6t69y5s+rUqaMOHTpo2bJlXm3Drl27dNtttykxMVF2u11t27bV888/r9LSUkm/78/vvvtOH3/8sev3bOfOnV7df5m6deuqVatW+umnn9yWV/T7snHjRl155ZWqX7++YmNj1alTJ/3jH//w+PNauXKlRo4cqQYNGqhu3bq6/fbbdfToUeXn5+vGG29UvXr15HQ6NWHCBB0/ftztPjIyMtS9e3fVr19fdevWVefOnTV79uxyb2+lpqZqyJAhWrZsmTp37qy4uDi1adNGr7/+usc5rVq1SqNHj1ZCQoIaNGiga6+9Vnv37nUb6+ltI3/P53Q7d+5Uw4YNXY9Vtj9HjBihtWvXymaz6Z133im33rx582Sz2bRhwwZJvz/XtmzZov79++uss85Sw4YNNWbMGP36669u6xpjNH36dF1wwQWKi4vTOeeco+uvv17/+c9/Kp0rQsAAZ+Cuu+4yksx9991nPvvsM1NcXOxx3JdffmlatGhhOnXqZLKzs012drb58ssvjTHGHD161FxwwQUmISHBTJ061fzP//yPefHFF43D4TD9+vUzpaWlrvuRZFJTU023bt3MP/7xD7N06VLTp08fU7t2bfP999+7xuXm5pqzzz7bpKammhkzZph///vfZv78+ebGG280hYWFxhhjVq1aZSSZVatWudabNGmSqeppsWXLFvPYY48ZSWbOnDkmOzvbfPfdd8YYY55++mkjydxyyy1myZIlZt68eaZFixbG4XCYb775xnUfw4cPN9HR0SY1NdVkZmaaf//73+aTTz6p8DFXr15tHnzwQfP++++b1atXm8WLF5urr77axMXFmW3btlU6X19+bllZWSY6Otp06dLFLFiwwHzwwQdmwIABxmazmXfffdc1bs6cOUaSady4sbn77rvNxx9/bN5//31z4sQJ07t3b9OgQQPTunVrM3v2bPPJJ5+YIUOGGEkmIyPDdOjQwbzzzjtm6dKlpkePHsZut5s9e/ZUOv99+/aZxo0bm4YNG5oZM2aYZcuWmTFjxhhJZvTo0cYYYwoKCkx2drZJTk42F110kev37NixYxXeb0pKihk8eLDbsuPHj5vk5GTToUMHt+Wefl9WrlxpYmJiTK9evcyCBQvMsmXLzIgRI1y/G6f/vJo3b24efPBBs3z5cvPss8+aqKgoc8stt5jOnTubJ5980qxYscI8/PDDRpJ5/vnn3R5/xIgRZvbs2WbFihVmxYoV5q9//auJi4szGRkZ5bapSZMmpl27dmbevHnmk08+MTfccIORZFavXl1uTi1atDD33Xef+eSTT8x///d/m3POOcf07dvX7T579+5tevfuHdD5nO7YsWNm2bJlRpK58847Xfuz7LnWqVMnc9FFF5Vb78ILLzQXXnih6/rw4cNNTEyMadasmXnqqafM8uXLzeTJk03t2rXNkCFD3NYdOXKkiY6ONg8++KBZtmyZefvtt02bNm1MUlKSyc/Pr3CuCD7iBWdk//795uKLLzaSjCQTHR1tevbsaTIzM83hw4fdxp533nnl/gAaY0xmZqapVauW2bBhg9vy999/30gyS5cudS2TZJKSklwBYowx+fn5platWiYzM9O1rF+/fqZevXpm3759Fc69uvFizO9/+E+d88GDB01cXJwZNGiQ29hdu3YZu91ubr31Vtey4cOHG0nm9ddfr/KxPDlx4oQpLi42LVu2NOPGjatyvLc/tx49epjExES3fXfixAnTvn1706RJE1dIlm3/7bffXu6xevfubSSZjRs3upYdOHDAREVFmbi4OLdQyc3NNZLM3//+90rn/8gjjxhJ5vPPP3dbPnr0aGOz2cz27dtdyzwFSUVSUlLMoEGDzPHjx83x48fNDz/84HoB++ijj9zGevp9adOmjenUqZM5fvy429ghQ4YYp9NpSkpKjDG//7zuu+8+t3FXX321kWSmTp3qtvyCCy4wnTt3rnDeJSUl5vjx4+aJJ54wDRo0cAv8lJQUExsba3744QfXst9++83Ur1/f3HPPPa5lZXO699573e57ypQpRpLJy8tzLfMUL/6ejyc///yzkWQmTZpU7ray+efk5LiWrV+/3kgyb7zxhmtZ2XPtxRdfdFv/qaeeMpLMp59+aowxJjs722M07t6928TFxZmHHnqo0rkiuHjbCGekQYMGWrt2rTZs2KBnnnlGV111lb755hulp6erQ4cO2r9/f5X38dFHH6l9+/a64IILdOLECdfl8ssv93iYvm/fvm4fjk1KSlJiYqJ++OEHSSc/O7J69WrdeOONrsPOwZCdna3ffvut3AeSmzZtqn79+unf//53uXWuu+46r+77xIkTevrpp9WuXTvFxMSodu3aiomJ0bffflvuLamKVPVzO3r0qD7//HNdf/31Ovvss13joqKiNGzYMP3444/avn27V/N3Op3q0qWL63r9+vWVmJioCy64QI0aNXItb9u2rSS55lCRlStXql27duU+XzRixAgZY7Ry5cpK16/M0qVLFR0drejoaKWkpOi1117TSy+9pMGDB1e63nfffadt27bpj3/8oyS5/e4OGjRIeXl55X5eQ4YMcbtetv2nP1bbtm3L/UxWrlypSy+9VA6HQ1FRUYqOjtbjjz+uAwcOaN++fW5jL7jgAjVr1sx1PTY2Vq1atfL4c77yyivdrp9//vmSvNsngZiPt2655RYlJibqlVdecS176aWX1LBhQ910003lxpftpzK33nqrJGnVqlWSTv4dstlsuu2229z2ZXJysjp27Mg3nsIM8QK/6Nq1qx5++GG999572rt3r8aNG6edO3d69aHdn376SV999ZXrBaTsEh8fL2NMuQBq0KBBufuw2+367bffJEkHDx5USUmJmjRp4p+N89KBAwckyeO3Wxo1auS6vUydOnVUt25dr+57/Pjx+stf/qKrr75aH374oT7//HNt2LBBHTt2dG13Vbz5uRljKpy/pHLbUNE3eerXr19uWUxMTLnlMTExkqRjx45VOvcDBw74NC9fXHzxxdqwYYM+++wzvfnmm0pNTdWYMWP06aefVrpe2WdiJkyYUO53995775Wkcr+7FW2/p+Wn/kzWr1+vAQMGSJJee+01/e///q82bNigRx99VJLK/Q5Uta8rG2u32z3e56kCOR9v2e123XPPPXr77bd16NAh/fzzz/rHP/6hu+66y7UNZWrXrl1uDsnJyZJ+/9356aefZIxRUlJSuf352WefefUPMQQP3zaC30VHR2vSpEl64YUXtHnz5irHJyQkKC4ursIP8CUkJPj0+PXr11dUVJR+/PFHn9Y7U2V/HPPy8srdtnfv3nLb4cu3cubPn6/bb79dTz/9tNvy/fv3q169er5P1oNzzjlHtWrVqnD+Uvl9EaxvFjVo0MCnefnC4XCoa9eukk5+ALt79+7q2LGj7r33XuXm5rp98PxUZY+Znp6ua6+91uOY1q1bV3tep3r33XcVHR2tjz76SLGxsa7lH3zwgV/u36rzGT16tJ555hm9/vrrOnbsmE6cOKFRo0aVG3fixAkdOHDALWDy8/Ml/f68TUhIkM1m09q1a8vFjySPyxA6HHnBGfH0giLJ9VbGqW8RVPQvrSFDhuj7779XgwYN1LVr13KXU7+V5I24uDj17t1b7733XlD/tZSWlqa4uDjNnz/fbfmPP/6olStXVvptoqrYbLZyfzyXLFmiPXv2VPs+T3fWWWepe/fuWrRokdt+Ki0t1fz589WkSRO1atXKb4/ni/79+2vr1q368ssv3ZaXfbOkb9++fnusli1b6qGHHtKmTZu0YMGCCse1bt1aLVu21P/93/95/L3t2rWr3879Y7PZVLt2bUVFRbmW/fbbb3rzzTf9cv/hOp+qjgI5nU7dcMMNmj59umbMmKGhQ4e6vT11qrfeesvt+ttvvy3p9/PHDBkyRMYY7dmzx+O+7NChg5+2Cv7AkReckcsvv1xNmjTR0KFD1aZNG5WWlio3N1fPP/+8zj77bD3wwAOusR06dNC7776rBQsWqEWLFoqNjVWHDh00duxYLVy4UJdcconGjRun888/X6Wlpdq1a5eWL1+uBx98UN27d/dpXlOnTtXFF1+s7t2765FHHtG5556rn376Sf/61780c+bMgJxQrl69evrLX/6iiRMn6vbbb9ctt9yiAwcOKCMjQ7GxsZo0aVK173vIkCGaO3eu2rRpo/PPP19ffPGF/va3v/n9rbHMzExddtll6tu3ryZMmKCYmBhNnz5dmzdv1jvvvBOyc7iMGzdO8+bN0+DBg/XEE08oJSVFS5Ys0fTp0zV69Gi/R9WECRM0Y8YMZWRk6MYbb3R7kT7VzJkzNXDgQF1++eUaMWKEGjdurF9++UVff/21vvzyS7333nt+mc/gwYM1depU3Xrrrbr77rt14MABPffccyE7GhCs+cTHxyslJUX//Oc/1b9/f9WvX18JCQlu/6B54IEHXH8fyr6Of7qYmBg9//zzOnLkiC688EKtW7dOTz75pAYOHKiLL75Y0smTFd5999264447tHHjRl1yySU666yzlJeXp08//VQdOnTQ6NGj/bp9qD7iBWfkscce0z//+U+98MILysvLU1FRkZxOpy699FKlp6e7PpAonTxXQ15enkaOHKnDhw8rJSVFO3fu1FlnnaW1a9fqmWee0axZs7Rjxw7FxcWpWbNmuvTSS30+8iJJHTt21Pr16zVp0iSlp6fr8OHDSk5OVr9+/VyfMwiE9PR0JSYm6u9//7sWLFiguLg49enTR08//bRatmxZ7ft98cUXFR0drczMTB05ckSdO3fWokWL9Nhjj/lx9lLv3r21cuVKTZo0SSNGjFBpaak6duyof/3rX+U+bBpMDRs21Lp165Senq709HQVFhaqRYsWmjJlisaPH+/3xzv77LP1+OOP689//rPeeust3X777R7H9e3bV+vXr9dTTz2lsWPH6uDBg2rQoIHatWunG2+80W/z6devn15//XU9++yzGjp0qBo3bqyRI0cqMTFRd955p98eJxznM3v2bP3Xf/2XrrzyShUVFWn48OFu/5uRbt26KTU1VXFxcRUe3Sx7i+v+++/Xk08+qbi4OI0cOVJ/+9vf3MbNnDlTPXr00MyZMzV9+nSVlpaqUaNGuuiii/x2Mkr4h82Y084oBACARXz11Vfq2LGjXnnlFdcHpU81YsQIvf/++zpy5EgIZodA4cgLAMByvv/+e/3www+aOHGinE5n2Pw/0xAcfGAXAGA5f/3rX3XZZZfpyJEjeu+99yz9/wWD73jbCAAAWApHXgAAgKUQLwAAwFKIFwAAYCkR922j0tJS7d27V/Hx8SE7oRYAAPCNMUaHDx9Wo0aNKvzfcpSJuHjZu3evmjZtGuppAACAati9e3eVZw+PuHgpO+377t27vf4/9gIAgNAqLCxU06ZNvfrft0RcvJS9VVS3bl3iBQAAi/HmIx98YBcAAFgK8QIAACyFeAEAAJZCvAAAAEshXgAAgKUQLwAAwFKIFwAAYCnECwAAsBTiBQAAWEpA42XNmjUaOnSoGjVqJJvNpg8++KDS8VlZWbLZbOUu27ZtC+Q0AQCAhQT0fw9w9OhRdezYUXfccYeuu+46r9fbvn2726n9GzZsGIjpAQAACwpovAwcOFADBw70eb3ExETVq1fPq7FFRUUqKipyXS8sLPT58QAAgHWE5WdeOnXqJKfTqf79+2vVqlWVjs3MzJTD4XBdmjZtGqRZAgCAUAireHE6nZo1a5YWLlyoRYsWqXXr1urfv7/WrFlT4Trp6ekqKChwXXbv3h3EGQMAgGAL6NtGvmrdurVat27tup6Wlqbdu3frueee0yWXXOJxHbvdLrvdHqwpAgCAEAurIy+e9OjRQ99++22opwEAAMJE2MdLTk6OnE5nqKcBAADCREDfNjpy5Ii+++471/UdO3YoNzdX9evXV7NmzZSenq49e/Zo3rx5kqRp06YpNTVV5513noqLizV//nwtXLhQCxcuDOQ0AQCAhQQ0XjZu3Ki+ffu6ro8fP16SNHz4cM2dO1d5eXnatWuX6/bi4mJNmDBBe/bsUVxcnM477zwtWbJEgwYNCuQ0AQCAhdiMMSbUk/CnwsJCORwOFRQUuJ3oDgAAhC9fXr/D/jMvAAAApyJeAACApRAvAADAUogXAABgKcQLAACwFOIFAABYCvECAAAshXgBAACWQrwAAABLIV4AAIClEC8AAMBSiBcAAGApxAsAALAU4gUAAFgK8QIAACyFeAEAAJZCvAAAAEshXgAAgKUQLwAAwFKIFwAAYCnECwAAsBTiBQAAWArxAgAALIV4AQAAlkK8AAAASyFeAACApRAvAADAUogXAABgKcQLAACwFOIFAABYCvECAAAshXgBAACWQrwAAABLIV4AAIClEC8AAMBSiBcAAGApxAsAALAU4gUAAFgK8QIAACyFeAEAAJZCvAAAAEshXgAAgKUQLwAAwFKIFwAAYCm1Qz0BAEDNUHyiVG9m79QPv/yqlPp1NCwtVTG1+Te0lZSUGq3f8Yv2HT6mxPhYdWteX1G1bEGfR0DjZc2aNfrb3/6mL774Qnl5eVq8eLGuvvrqStdZvXq1xo8fry1btqhRo0Z66KGHNGrUqEBOEwAQYJlLt+q1tTtUan5f9tTSrzWyV3OlD2oXuonBa8s25ynjw63KKzjmWuZ0xGrS0Ha6or0zqHMJaPIePXpUHTt21Msvv+zV+B07dmjQoEHq1auXcnJyNHHiRN1///1auHBhIKcJAAigzKVbNXONe7hIUqmRZq7ZocylW0MzMXht2eY8jZ7/pVu4SFJ+wTGNnv+llm3OC+p8AnrkZeDAgRo4cKDX42fMmKFmzZpp2rRpkqS2bdtq48aNeu6553TdddcFaJYAgEApPlGq19buqHTMrDU71DmlvqKjgv/2A6pWWiqlL9ok4+E2I8kmKePDrbqsXXLQ3kIKq8+8ZGdna8CAAW7LLr/8cs2ePVvHjx9XdHR0uXWKiopUVFTkul5YWBjweQIAvPNm9s5yR1xOZyTd8+YXQZkP/M9Iyis4pvU7flHaHxoE5THDKl7y8/OVlJTktiwpKUknTpzQ/v375XSWf08tMzNTGRkZwZoiAMAHX+d79w/KBmfFqPE5cQGeDarj4K/F2v3Lb1WO23f4WJVj/CWs4kWSbDb3Q07GGI/Ly6Snp2v8+PGu64WFhWratGngJggAqNJPhcf0atb3Wpyz16vx9/b5g+7s1SLAs0J1ZH9/QLe89lmV4xLjY4Mwm5PCKl6Sk5OVn5/vtmzfvn2qXbu2GjTwfCjKbrfLbrcHY3oAgCqURcvb63ep+ESpV+vUsknD0lIDOzFU28GjxZXebpOU7Dj5telgCat4SUtL04cffui2bPny5eratavHz7sAAMKDp2jpmnKOxl3WSqu3/6RZa3dWuO7IXs0530uYWrY5T/e/m1Ph7WXviUwa2i6o53sJaLwcOXJE3333nev6jh07lJubq/r166tZs2ZKT0/Xnj17NG/ePEnSqFGj9PLLL2v8+PEaOXKksrOzNXv2bL3zzjuBnCYAoJoqi5aef2ggm82mi85NkM1mK3eel1o2cZ6XMLZsc57GvJ2jE6VGV13QSJe3S9Jfl3zt9nXp5BCd58Vmyj5UEgBZWVnq27dvueXDhw/X3LlzNWLECO3cuVNZWVmu21avXq1x48a5TlL38MMP+3SSusLCQjkcDhUUFKhu3br+2AwAwGm8iZbTcYZd6zg9XKbeeIGiatkCeoZdX16/AxovoUC8AEDgVCdaYC0VhUug+fL6HVafeQEAhCeipWYIVbj4ingBAFSIaKk5rBIuEvECAPCAaKlZrBQuEvECADgF0VLzWC1cJOIFACCipaayYrhIxAsA1GhES81l1XCRiBcAqJGIlprNyuEiES8AUKMQLbB6uEjECwDUCEQLpMgIF4l4AYCIRrSgTKSEi0S8AEBEIlpwqkgKF4l4AYCIQrTgdJEWLhLxAgARgWiBJ5EYLhLxAgCWRrSgIpEaLhLxAgCWRLSgMpEcLhLxAgCWQrSgKpEeLhLxAgCWQLTAGzUhXCTiBQDCGtECb9WUcJGIFwAIS0QLfFGTwkUiXgAgrBAt8FVNCxeJeAGAsEC0oDpqYrhIxAsAhBTRguqqqeEiES8AEBJEC85ETQ4XiXgBgKAiWnCmanq4SMQLAAQF0QJ/IFxOIl4AIICIFvgL4fI74gUAAoBogT8RLu6IFwDwI6IF/ka4lEe8AIAfEC0IBMLFM+IFAM4A0YJAIVwqRrwAQDUQLQgkwqVyxAsA+IBoQaARLlUjXgDAC0QLgoFw8Q7xAgCVIFoQLISL94gXAPCAaEEwES6+IV4A4BREC4KNcPEd8QIAIloQGoRL9RAvAGo0ogWhQrhUH/ECoEYiWhBKhMuZIV4A1ChEC0KNcDlzxAuAGoFoQTggXPyDeAEQ0YgWhAvCxX+IFwARiWhBOCFc/It4ARBRiBaEG8LF/4gXABGBaEE4IlwCo1YwHmT69Olq3ry5YmNj1aVLF61du7bCsVlZWbLZbOUu27ZtC8ZUAVjMT4XHNPlfW9RryirNXbdTxSdK1TXlHL11V3e9NypNF52bQLggJAiXwAn4kZcFCxZo7Nixmj59ui666CLNnDlTAwcO1NatW9WsWbMK19u+fbvq1q3rut6wYcNATxWAhXCkBeGMcAksmzHGBPIBunfvrs6dO+vVV191LWvbtq2uvvpqZWZmlhuflZWlvn376uDBg6pXr57Pj1dYWCiHw6GCggK3+AEQGYgWhDvCpXp8ef0O6JGX4uJiffHFF3rkkUfclg8YMEDr1q2rdN1OnTrp2LFjateunR577DH17dvX47iioiIVFRW5rhcWFp75xAGEHaIFVkC4BEdA42X//v0qKSlRUlKS2/KkpCTl5+d7XMfpdGrWrFnq0qWLioqK9Oabb6p///7KysrSJZdcUm58ZmamMjIyAjJ/AKFHtMAqCJfgCcq3jU7/42KMqfAPTuvWrdW6dWvX9bS0NO3evVvPPfecx3hJT0/X+PHjXdcLCwvVtGlTP80cQKgQLbASwiW4AhovCQkJioqKKneUZd++feWOxlSmR48emj9/vsfb7Ha77Hb7Gc0TQPggWmA1hEvwBTReYmJi1KVLF61YsULXXHONa/mKFSt01VVXeX0/OTk5cjqdgZgigDBBtMCKCJfQCPjbRuPHj9ewYcPUtWtXpaWladasWdq1a5dGjRol6eTbPnv27NG8efMkSdOmTVNqaqrOO+88FRcXa/78+Vq4cKEWLlwY6KkCCAGiBVZFuIROwOPlpptu0oEDB/TEE08oLy9P7du319KlS5WSkiJJysvL065du1zji4uLNWHCBO3Zs0dxcXE677zztGTJEg0aNCjQUwUQREQLrIxwCa2An+cl2DjPCxDeiBZYHeESGGFznhcAKEO0IBIQLuGBeAEQUEQLIgXhEj6IFwABQbQgkhAu4YV4AeBXRAsiDeESfogXAH5BtCASES7hiXgBcEaIFkQqwiV8ES8AqoVoQSQjXMIb8QLAJ0QLIh3hEv6IFwBeIVpQExAu1kC8AKgU0YKagnCxDuIFgEdEC2oSwsVaiBcAbogW1DSEi/UQLwAkES2omQgXayJegBqOaEFNRbhYF/EC1FBEC2oywsXaiBeghiFaUNMRLtZHvAA1BNECEC6RgngBIhzRApxEuEQO4gWIUEQL8DvCJbIQL0CEIVoAd4RL5CFegAhBtADlES6RiXgBLI5oATwjXCIX8QJYFNECVIxwiWzEC2AxRAtQOcIl8hEvgEUQLUDVCJeagXgBwhzRAniHcKk5iBcgTBEtgPcIl5qFeAHCDNEC+IZwqXmIFyBMEC2A7wiXmol4AUKMaAGqh3CpuYgXIESIFqD6CJeajXgBgoxoAc4M4QLiBQgSogU4c4QLJOIFCDiiBfAPwgVliBcgQIgWwH8IF5yKeAH8jGgB/ItwwemIF8BPiBbA/wgXeEK8AGeIaAECg3BBRYgXoJqIFiBwCBdUhngBfES0AIFFuKAqxAvgJaIFCDzCBd4gXoAqEC1AcBAu8BbxAlSAaAGCh3CBL4gX4DRECxBchAt8RbygxigpNVq/4xftO3xMifGx6ta8vtsfSKIlvJWUlGjt2rXKy8uT0+lUr169FBUVFeppwQeenoMrtuYTLvAZ8YIaYdnmPGV8uFV5Bcdcy5yOWE0a2k6dmp1DtIS5RYsW6YEHHtCPP/7oWtakSRO9+OKLuvbaa0M4M3jL03OwXp1oFf52XKVGhAt8UisYDzJ9+nQ1b95csbGx6tKli9auXVvp+NWrV6tLly6KjY1VixYtNGPGjGBMExFq2eY8jZ7/pdsfTUnKKzimUfO/1EXPrNTcdTtVfKJUXVPO0Vt3ddd7o9J00bkJhEsYWLRoka6//nq3cJGkPXv26Prrr9eiRYtCNDN4q6Ln4KFfT4bLhannEC7wScDjZcGCBRo7dqweffRR5eTkqFevXho4cKB27drlcfyOHTs0aNAg9erVSzk5OZo4caLuv/9+LVy4MNBTRQQqKTXK+HCrTCVjTpQadSFawlJJSYkeeOABGVN+D5YtGzt2rEpKSoI9NXjJm+fg7l9+Ddp8EBkCHi9Tp07VnXfeqbvuuktt27bVtGnT1LRpU7366qsex8+YMUPNmjXTtGnT1LZtW911113605/+pOeee87j+KKiIhUWFrpdgDLrd/xS7l97njx4WSuiJQytXbu23BGXUxljtHv37iqP5iJ0vHkO5hcWaf2OX4I0I0SCgMZLcXGxvvjiCw0YMMBt+YABA7Ru3TqP62RnZ5cbf/nll2vjxo06fvx4ufGZmZlyOByuS9OmTf23AbC8fYerDhdJ+vlIUYBngurIy8vz6zgEn7fPQW/HAVKA42X//v0qKSlRUlKS2/KkpCTl5+d7XCc/P9/j+BMnTmj//v3lxqenp6ugoMB12b17t/82AJaXGB/r13EILqfT6ddxCK6SUqPv9x3xaizPQfgiKN82Ov1QvDGm0sPznsZ7Wi5JdrtddrvdD7NEJOrYxCF77Voq+v/fIjqdTVKy4+RXNhF+evXqpSZNmmjPnj0eP/dis9nUpEkT9erVKwSzQ0VKSo2WbsrT3//9rb6tIl54DqI6AnrkJSEhQVFRUeWOsuzbt6/c0ZUyycnJHsfXrl1bDRo0CNhcEXmOHS/RqLe+rDRcJGnS0HZ8yyFMRUVF6cUXX5RU/h8vZdenTZvG+V7CREmp0Yf/t1dXTFuj+97J0bf7jqhubG0NPf/kkbHTn2U8B1FdAY2XmJgYdenSRStWrHBbvmLFCvXs2dPjOmlpaeXGL1++XF27dlV0dHTA5orIcux4ie5+8wut+eZnxUVH6cHLWsrpcD8sneyI1au3ddYV7XnLIZxde+21ev/999W4cWO35U2aNNH777/PeV7CQEXRMu7SVlr7cD+9dGtnzbits5J5DsJPbMbTsVg/WrBggYYNG6YZM2YoLS1Ns2bN0muvvaYtW7YoJSVF6enp2rNnj+bNmyfp5Fel27dvr3vuuUcjR45Udna2Ro0apXfeeUfXXXddlY9XWFgoh8OhgoIC1a1bN5CbhjB1erjMveNCdW/RoMoz7CK8cYbd8OPp7aG6sbV158UtNOKiVDniosuN5zmIivjy+h3wz7zcdNNNOnDggJ544gnl5eWpffv2Wrp0qVJSUiSd/JbAqed8ad68uZYuXapx48bplVdeUaNGjfT3v//dq3ABKgoXSYqqZVPaH3jr0aqioqLUp0+fUE8D8j1ayvAchL8E/MhLsHHkpeaqLFwAnLnqRgvgjbA68gIEA+ECBA7RgnBDvMDyCBcgMIgWhCviBZZGuAD+R7Qg3BEvsCzCBfAvogVWQbzAkggXwH+IFlgN8QLLIVwA/yBaYFXECyyFcAHOHNECqyNeYBmEC3BmiBZECuIFlkC4ANVHtCDSEC8Ie4QLUD1ECyIV8YKwRrgAviNaEOmIF4QtwgXwDdGCmoJ4QVgiXADvES2oaYgXhB3CBfAO0YKainhBWCFcgKoRLajpiBeEDcIFqBzRApxEvCAsEC5AxYgWwB3xgpAjXADPiBbAM+IFIUW4AOURLUDliBeEDOECuCNaAO8QLwgJwgX4HdEC+IZ4QdARLsBJRAtQPcQLgopwAYgW4EwRLwgawgU1HdEC+AfxgqAgXFCTES2AfxEvCDjCBTUV0QIEBvGCgCJcUBMRLUBgES8IGMIFNQ3RAgQH8YKAIFxQkxAtQHARL/A7wgU1BdEChAbxAr8iXFATEC1AaBEv8BvCBZGOaAHCA/ECvyBcEMmIFiC8EC84Y4QLIhXRAoQn4gVnhHBBJCJagPBGvKDaCBdEGqIFsAbiBdVCuCCSEC2AtRAv8BnhgkhBtADWRLzAJ4QLIgHRAlgb8QKvES6wOqIFiAzEC7xCuMDKiBYgshAvqBLhAqsiWoDIRLygUoQLrIhoASIb8YIKES6wGqIFqBmIF3hEuMBKiBagZqkVyDs/ePCghg0bJofDIYfDoWHDhunQoUOVrjNixAjZbDa3S48ePQI5TZyGcIFVlJQaffh/e3XFtDW6750cfbvviOrG1ta4S1tp7cP99MClLQkXIAIF9MjLrbfeqh9//FHLli2TJN19990aNmyYPvzww0rXu+KKKzRnzhzX9ZiYmEBOE6cgXGAFHGkBaraAxcvXX3+tZcuW6bPPPlP37t0lSa+99prS0tK0fft2tW7dusJ17Xa7kpOTAzU1VIBwQbgjWgBIAYyX7OxsORwOV7hIUo8ePeRwOLRu3bpK4yUrK0uJiYmqV6+eevfuraeeekqJiYkexxYVFamoqMh1vbCw0H8bUYMQLghnRAuAUwUsXvLz8z0GR2JiovLz8ytcb+DAgbrhhhuUkpKiHTt26C9/+Yv69eunL774Qna7vdz4zMxMZWRk+HXuNQ3hgnBFtADwxOd4mTx5cpWxsGHDBkmSzWYrd5sxxuPyMjfddJPrv9u3b6+uXbsqJSVFS5Ys0bXXXltufHp6usaPH++6XlhYqKZNm1a5HTiJcEE4IloAVMbneBkzZoxuvvnmSsekpqbqq6++0k8//VTutp9//llJSUleP57T6VRKSoq+/fZbj7fb7XaPR2RQNcIF4YZoAeANn+MlISFBCQkJVY5LS0tTQUGB1q9fr27dukmSPv/8cxUUFKhnz55eP96BAwe0e/duOZ1OX6eKShAuCCdECwBf2IwxJlB3PnDgQO3du1czZ86UdPKr0ikpKW5flW7Tpo0yMzN1zTXX6MiRI5o8ebKuu+46OZ1O7dy5UxMnTtSuXbv09ddfKz4+vsrHLCwslMPhUEFBgerWrRuoTbM0wgXhgmgBUMaX1++Anuflrbfe0v33368BAwZIkq688kq9/PLLbmO2b9+ugoICSVJUVJQ2bdqkefPm6dChQ3I6nerbt68WLFjgVbigaoQLwgHRAuBMBPTISyhw5KVihAtCjWgBUJGwOfKC8EG4IJSIFgD+RLzUAIQLQoVoARAIxEuEI1wQCkQLgEAiXiIY4YJgI1oABAPxEqEIFwQT0QIgmIiXCES4IFiIFgChQLxEGMIFwUC0AAgl4iWCEC4INKIFQDggXiIE4YJAIloAhBPiJQIQLggUogVAOCJeLI5wQSAQLQDCGfFiYYQL/I1oAWAFxItFES7wJ6IFgJUQLxZEuMBfiBYAVkS8WAzhAn8gWgBYGfFiIYQLzhTRAiASEC8WQbjgTBAtACIJ8WIBhAuqi2gBEImIlzBHuKA6iBYAkYx4CWOEC3xFtACoCYiXMEW4wBdEC4CahHgJQ4QLvEW0AKiJiJcwQ7jAG0QLgJqMeAkjhAuqQrQAAPESNggXVIZoAYDfES9hgHBBRYgWACiPeAkxwgWeEC0AUDHiJYQIF5yOaAGAqhEvIUK44FRECwB4j3gJAcIFZYgWAPAd8RJkhAskogUAzgTxEkSEC4gWADhzxEuQEC41G9ECAP5DvAQB4VJzES0A4H/ES4ARLjUT0QIAgUO8BBDhUvMQLQAQeMRLgBAuNQvRAgDBQ7wEAOFScxAtABB8xIufES41A9ECAKFDvPgR4RL5iBYACD3ixU8Il8hGtABA+CBe/IBwiVxECwCEH+LlDBEukYloAYDwRbycAcIl8hAtABD+iBcvlZQard/xi/YdPqbE+Fh1bOLQqLe+JFws5PR92K15fUXVsrluI1oAwBoCGi9PPfWUlixZotzcXMXExOjQoUNVrmOMUUZGhmbNmqWDBw+qe/fueuWVV3TeeecFcqqVWrY5TxkfblVewTHXMnvtWio6UUq4WISnfeh0xOovg9uqxIhoAQALCWi8FBcX64YbblBaWppmz57t1TpTpkzR1KlTNXfuXLVq1UpPPvmkLrvsMm3fvl3x8fGBnK5HyzbnafT8L2VOW150olSSdG+fFoRLmKtoH+YVHNO9b+e4rhMtAGANAY2XjIwMSdLcuXO9Gm+M0bRp0/Too4/q2muvlSS98cYbSkpK0ttvv6177rknUFP1qKTUKOPDreVe9E719vrdurdvS9fbDwgv3uxDm6QH+rfUHRc3J1oAwALC6jMvO3bsUH5+vgYMGOBaZrfb1bt3b61bt85jvBQVFamoqMh1vbCw0G/zWb/jF7e3GTzJKzimS6asUlxMlN8eF/7zW3FJlfvQSOreogHhAgAWEVbxkp+fL0lKSkpyW56UlKQffvjB4zqZmZmuIzz+tu9w5S96ZfYc+i0gj4/g8XZfAwBCz+d4mTx5cpWxsGHDBnXt2rXak7LZ3N+CMcaUW1YmPT1d48ePd10vLCxU06ZNq/3Yp0qMj/Vq3F+GtNN5jer65THhX1v2FuqvH22tcpy3+xoAEHo+x8uYMWN08803VzomNTW1WpNJTk6WdPIIjNPpdC3ft29fuaMxZex2u+x2e7UeryrdmteX0xGr/IJjHj8zYZOU7IjViJ6pfOYlTF2YWl//vfY/Ve7Dbs3rB3tqAIBq8jleEhISlJCQEIi5qHnz5kpOTtaKFSvUqVMnSSe/sbR69Wo9++yzAXnMykTVsmnS0HYaPf9L2SS3F7+yVJk0tB3hEsbYhwAQeWoF8s537dql3Nxc7dq1SyUlJcrNzVVubq6OHDniGtOmTRstXrxY0sm3i8aOHaunn35aixcv1ubNmzVixAjVqVNHt956ayCnWqEr2jv16m2dlexwf1sh2RGrV2/rrCvaOytYE+GCfQgAkSWgH9h9/PHH9cYbb7iulx1NWbVqlfr06SNJ2r59uwoKClxjHnroIf3222+69957XSepW758eUjO8VLmivZOXdYuucKzsyL8sQ8BIHLYjDGVnQLDcgoLC+VwOFRQUKC6dfkQLQAAVuDL63dA3zYCAADwN+IFAABYCvECAAAshXgBAACWQrwAAABLIV4AAIClEC8AAMBSiBcAAGApxAsAALAU4gUAAFgK8QIAACyFeAEAAJZCvAAAAEshXgAAgKUQLwAAwFKIFwAAYCnECwAAsBTiBQAAWArxAgAALIV4AQAAlkK8AAAASyFeAACApRAvAADAUogXAABgKcQLAACwFOIFAABYCvECAAAshXgBAACWQrwAAABLIV4AAIClEC8AAMBSiBcAAGApxAsAALAU4gUAAFgK8QIAACyFeAEAAJZCvAAAAEshXgAAgKUQLwAAwFKIFwAAYCnECwAAsBTiBQAAWArxAgAALIV4AQAAlkK8AAAASyFeAACApQQ0Xp566in17NlTderUUb169bxaZ8SIEbLZbG6XHj16BHKaAADAQgIaL8XFxbrhhhs0evRon9a74oorlJeX57osXbo0QDMEAABWUzuQd56RkSFJmjt3rk/r2e12JScnB2BGAADA6sLyMy9ZWVlKTExUq1atNHLkSO3bt6/CsUVFRSosLHS7AACAyBV28TJw4EC99dZbWrlypZ5//nlt2LBB/fr1U1FRkcfxmZmZcjgcrkvTpk2DPGMAABBMPsfL5MmTy32g9vTLxo0bqz2hm266SYMHD1b79u01dOhQffzxx/rmm2+0ZMkSj+PT09NVUFDguuzevbvajw0AAMKfz595GTNmjG6++eZKx6SmplZ3PuU4nU6lpKTo22+/9Xi73W6X3W732+MBAIDw5nO8JCQkKCEhIRBz8ejAgQPavXu3nE5n0B4TAACEr4B+5mXXrl3Kzc3Vrl27VFJSotzcXOXm5urIkSOuMW3atNHixYslSUeOHNGECROUnZ2tnTt3KisrS0OHDlVCQoKuueaaQE4VAABYREC/Kv3444/rjTfecF3v1KmTJGnVqlXq06ePJGn79u0qKCiQJEVFRWnTpk2aN2+eDh06JKfTqb59+2rBggWKj48P5FQBAIBF2IwxJtST8KfCwkI5HA4VFBSobt26oZ4OAADwgi+v32H3VWkAAIDKEC8AAMBSiBcAAGApxAsAALAU4gUAAFgK8QIAACyFeAEAAJYS0JPUhULZaWsKCwtDPBMAAOCtstdtb04/F3HxcvjwYUlS06ZNQzwTAADgq8OHD8vhcFQ6JuLOsFtaWqq9e/cqPj5eNpvNr/ddWFiopk2bavfu3RF59t5I3z4p8reR7bO+SN/GSN8+KfK3MVDbZ4zR4cOH1ahRI9WqVfmnWiLuyEutWrXUpEmTgD5G3bp1I/IXskykb58U+dvI9llfpG9jpG+fFPnbGIjtq+qISxk+sAsAACyFeAEAAJZCvPjAbrdr0qRJstvtoZ5KQET69kmRv41sn/VF+jZG+vZJkb+N4bB9EfeBXQAAENk48gIAACyFeAEAAJZCvAAAAEshXgAAgKUQLwAAwFKIl0rs3LlTd955p5o3b664uDj94Q9/0KRJk1RcXFzpesYYTZ48WY0aNVJcXJz69OmjLVu2BGnWvnnqqafUs2dP1alTR/Xq1fNqnREjRshms7ldevToEdiJVlN1ts9K+0+SDh48qGHDhsnhcMjhcGjYsGE6dOhQpeuE8z6cPn26mjdvrtjYWHXp0kVr166tdPzq1avVpUsXxcbGqkWLFpoxY0aQZlp9vmxjVlZWuX1ls9m0bdu2IM7Ye2vWrNHQoUPVqFEj2Ww2ffDBB1WuY6V96Ov2WW3/ZWZm6sILL1R8fLwSExN19dVXa/v27VWuF+x9SLxUYtu2bSotLdXMmTO1ZcsWvfDCC5oxY4YmTpxY6XpTpkzR1KlT9fLLL2vDhg1KTk7WZZdd5vqfRoaT4uJi3XDDDRo9erRP611xxRXKy8tzXZYuXRqgGZ6Z6myflfafJN16663Kzc3VsmXLtGzZMuXm5mrYsGFVrheO+3DBggUaO3asHn30UeXk5KhXr14aOHCgdu3a5XH8jh07NGjQIPXq1Us5OTmaOHGi7r//fi1cuDDIM/eer9tYZvv27W77q2XLlkGasW+OHj2qjh076uWXX/ZqvNX2oa/bV8Yq+2/16tX685//rM8++0wrVqzQiRMnNGDAAB09erTCdUKyDw18MmXKFNO8efMKby8tLTXJycnmmWeecS07duyYcTgcZsaMGcGYYrXMmTPHOBwOr8YOHz7cXHXVVQGdj795u31W239bt241ksxnn33mWpadnW0kmW3btlW4Xrjuw27duplRo0a5LWvTpo155JFHPI5/6KGHTJs2bdyW3XPPPaZHjx4Bm+OZ8nUbV61aZSSZgwcPBmF2/iXJLF68uNIxVtyHZbzZPivvP2OM2bdvn5FkVq9eXeGYUOxDjrz4qKCgQPXr16/w9h07dig/P18DBgxwLbPb7erdu7fWrVsXjCkGRVZWlhITE9WqVSuNHDlS+/btC/WU/MJq+y87O1sOh0Pdu3d3LevRo4ccDkeV8w23fVhcXKwvvvjC7WcvSQMGDKhwW7Kzs8uNv/zyy7Vx40YdP348YHOtrupsY5lOnTrJ6XSqf//+WrVqVSCnGVRW24fVZdX9V1BQIEmVvu6FYh8SLz74/vvv9dJLL2nUqFEVjsnPz5ckJSUluS1PSkpy3WZ1AwcO1FtvvaWVK1fq+eef14YNG9SvXz8VFRWFempnzGr7Lz8/X4mJieWWJyYmVjrfcNyH+/fvV0lJiU8/+/z8fI/jT5w4of379wdsrtVVnW10Op2aNWuWFi5cqEWLFql169bq37+/1qxZE4wpB5zV9qGvrLz/jDEaP368Lr74YrVv377CcaHYhzUyXiZPnuzxA1SnXjZu3Oi2zt69e3XFFVfohhtu0F133VXlY9hsNrfrxphyywKlOtvni5tuukmDBw9W+/btNXToUH388cf65ptvtGTJEj9uRcUCvX1SaPef5Ns2eppXVfMN9T6sjK8/e0/jPS0PJ75sY+vWrTVy5Eh17txZaWlpmj59ugYPHqznnnsuGFMNCivuQ29Zef+NGTNGX331ld55550qxwZ7H9YOyL2GuTFjxujmm2+udExqaqrrv/fu3au+ffsqLS1Ns2bNqnS95ORkSSdL1Ol0upbv27evXJkGiq/bd6acTqdSUlL07bff+u0+KxPI7QuH/Sd5v41fffWVfvrpp3K3/fzzzz7NN9j70JOEhARFRUWVOwJR2c8+OTnZ4/jatWurQYMGAZtrdVVnGz3p0aOH5s+f7+/phYTV9qE/WGH/3XffffrXv/6lNWvWqEmTJpWODcU+rJHxkpCQoISEBK/G7tmzR3379lWXLl00Z84c1apV+cGq5s2bKzk5WStWrFCnTp0knXyfe/Xq1Xr22WfPeO7e8GX7/OHAgQPavXu324t9IAVy+8Jh/0neb2NaWpoKCgq0fv16devWTZL0+eefq6CgQD179vT68YK9Dz2JiYlRly5dtGLFCl1zzTWu5StWrNBVV13lcZ20tDR9+OGHbsuWL1+url27Kjo6OqDzrY7qbKMnOTk5Id1X/mS1fegP4bz/jDG67777tHjxYmVlZal58+ZVrhOSfRiwjwJHgD179phzzz3X9OvXz/z4448mLy/PdTlV69atzaJFi1zXn3nmGeNwOMyiRYvMpk2bzC233GKcTqcpLCwM9iZU6YcffjA5OTkmIyPDnH322SYnJ8fk5OSYw4cPu8acun2HDx82Dz74oFm3bp3ZsWOHWbVqlUlLSzONGzeOiO0zxlr7zxhjrrjiCnP++eeb7Oxsk52dbTp06GCGDBniNsYq+/Ddd9810dHRZvbs2Wbr1q1m7Nix5qyzzjI7d+40xhjzyCOPmGHDhrnG/+c//zF16tQx48aNM1u3bjWzZ8820dHR5v333w/VJlTJ12184YUXzOLFi80333xjNm/ebB555BEjySxcuDBUm1Cpw4cPu55nkszUqVNNTk6O+eGHH4wx1t+Hvm6f1fbf6NGjjcPhMFlZWW6veb/++qtrTDjsQ+KlEnPmzDGSPF5OJcnMmTPHdb20tNRMmjTJJCcnG7vdbi655BKzadOmIM/eO8OHD/e4fatWrXKNOXX7fv31VzNgwADTsGFDEx0dbZo1a2aGDx9udu3aFZoNqIKv22eMtfafMcYcOHDA/PGPfzTx8fEmPj7e/PGPfyz3tUwr7cNXXnnFpKSkmJiYGNO5c2e3r2gOHz7c9O7d2218VlaW6dSpk4mJiTGpqanm1VdfDfKMfefLNj777LPmD3/4g4mNjTXnnHOOufjii82SJUtCMGvvlH01+PTL8OHDjTHW34e+bp/V9l9Fr3mn/o0Mh31o+/+TBQAAsIQa+W0jAABgXcQLAACwFOIFAABYCvECAAAshXgBAACWQrwAAABLIV4AAIClEC8AAMBSiBcAAGApxAsAALAU4gUAAFjK/wORiATEP84MfQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Stencil for a norm of Riemannian type\"); plt.axis('equal')\n",
    "PlotStencil2(MakeStencil2(riemann2_5.m)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.897242Z",
     "iopub.status.busy": "2024-04-30T08:46:15.897160Z",
     "iopub.status.idle": "2024-04-30T08:46:15.899509Z",
     "shell.execute_reply": "2024-04-30T08:46:15.899267Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8230854637602008"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MinStencilScal2(MakeStencil2(riemann2_5.m),riemann2_5.m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Three dimensional construction\n",
    "\n",
    "We propose a stencil construction based on the preliminary computation of an $M$-obtuse superbase, see notebook [Selling](TensorSelling.ipynb).\n",
    "It can be shown that, for generic $M$, the stencil consists of the Voronoi vectors of $Z^3$ for the distance $\\|\\cdot\\|_M$, see the main introduction. \n",
    "\n",
    "The stencil is guaranteed to be acute, and to hold precisely than $14$ distinct vertices and $24$ facets, for any positive definite matrix $M$. This is independent of the condition number of $M$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.900942Z",
     "iopub.status.busy": "2024-04-30T08:46:15.900862Z",
     "iopub.status.idle": "2024-04-30T08:46:15.902877Z",
     "shell.execute_reply": "2024-04-30T08:46:15.902639Z"
    }
   },
   "outputs": [],
   "source": [
    "def MakeStencil3(m):\n",
    "    e0,e1,e2,e3 = Selling.ObtuseSuperbase(m).T\n",
    "    return [(u,u+v,u+v+w) for u,v,w in itertools.permutations([e0,e1,e2,e3],3)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.904196Z",
     "iopub.status.busy": "2024-04-30T08:46:15.904119Z",
     "iopub.status.idle": "2024-04-30T08:46:15.905792Z",
     "shell.execute_reply": "2024-04-30T08:46:15.905573Z"
    }
   },
   "outputs": [],
   "source": [
    "riemann3_5 = Riemann.needle(direction3,1,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.907077Z",
     "iopub.status.busy": "2024-04-30T08:46:15.906999Z",
     "iopub.status.idle": "2024-04-30T08:46:15.909765Z",
     "shell.execute_reply": "2024-04-30T08:46:15.909512Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.078434290817091"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MinStencilScal3(MakeStencil3(riemann3_5.m),riemann3_5.m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-30T08:46:15.911069Z",
     "iopub.status.busy": "2024-04-30T08:46:15.910991Z",
     "iopub.status.idle": "2024-04-30T08:46:15.915757Z",
     "shell.execute_reply": "2024-04-30T08:46:15.915514Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(array([ 0., -1., -2.]), array([-1., -1., -3.]), array([ 0.,  0., -1.])),\n",
       " (array([ 0., -1., -2.]), array([-1., -1., -3.]), array([-1., -1., -2.])),\n",
       " (array([ 0., -1., -2.]), array([1., 0., 0.]), array([ 0.,  0., -1.])),\n",
       " (array([ 0., -1., -2.]), array([1., 0., 0.]), array([1., 0., 1.])),\n",
       " (array([ 0., -1., -2.]), array([ 0., -1., -1.]), array([-1., -1., -2.])),\n",
       " (array([ 0., -1., -2.]), array([ 0., -1., -1.]), array([1., 0., 1.])),\n",
       " (array([-1.,  0., -1.]), array([-1., -1., -3.]), array([ 0.,  0., -1.])),\n",
       " (array([-1.,  0., -1.]), array([-1., -1., -3.]), array([-1., -1., -2.])),\n",
       " (array([-1.,  0., -1.]), array([0., 1., 1.]), array([ 0.,  0., -1.])),\n",
       " (array([-1.,  0., -1.]), array([0., 1., 1.]), array([0., 1., 2.])),\n",
       " (array([-1.,  0., -1.]), array([-1.,  0.,  0.]), array([-1., -1., -2.])),\n",
       " (array([-1.,  0., -1.]), array([-1.,  0.,  0.]), array([0., 1., 2.])),\n",
       " (array([1., 1., 2.]), array([1., 0., 0.]), array([ 0.,  0., -1.])),\n",
       " (array([1., 1., 2.]), array([1., 0., 0.]), array([1., 0., 1.])),\n",
       " (array([1., 1., 2.]), array([0., 1., 1.]), array([ 0.,  0., -1.])),\n",
       " (array([1., 1., 2.]), array([0., 1., 1.]), array([0., 1., 2.])),\n",
       " (array([1., 1., 2.]), array([1., 1., 3.]), array([1., 0., 1.])),\n",
       " (array([1., 1., 2.]), array([1., 1., 3.]), array([0., 1., 2.])),\n",
       " (array([0., 0., 1.]), array([ 0., -1., -1.]), array([-1., -1., -2.])),\n",
       " (array([0., 0., 1.]), array([ 0., -1., -1.]), array([1., 0., 1.])),\n",
       " (array([0., 0., 1.]), array([-1.,  0.,  0.]), array([-1., -1., -2.])),\n",
       " (array([0., 0., 1.]), array([-1.,  0.,  0.]), array([0., 1., 2.])),\n",
       " (array([0., 0., 1.]), array([1., 1., 3.]), array([1., 0., 1.])),\n",
       " (array([0., 0., 1.]), array([1., 1., 3.]), array([0., 1., 2.]))]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MakeStencil3(riemann3_5.m)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}