tango.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# True Asymptotic Natural Gradient Optimization\n",
    "# paper: https://arxiv.org/abs/1712.08449\n",
    "# this notebook recreates Figure 1.\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare and check some gaussian distribution functions.\n",
    "\n",
    "def pdf(x, mu, sigma):\n",
    "    return np.exp(-(x - mu)**2 / (2 * sigma**2)) / np.sqrt(2 * np.pi * sigma**2)\n",
    "\n",
    "def nll(x, mu, sigma):  # negative log-likelihood\n",
    "    return (x - mu)**2 / (2 * sigma**2) + np.log(2 * np.pi * sigma**2) / 2\n",
    "\n",
    "def d_nll_x(x, mu, sigma):  # gradient of nll wrt. x\n",
    "    return (x - mu) / sigma**2\n",
    "\n",
    "def d_nll_mu(x, mu, sigma):  # gradient of nll wrt. mu\n",
    "    return (mu - x) / sigma**2\n",
    "\n",
    "def d_nll_sigma(x, mu, sigma):  # gradient of nll wrt. sigma\n",
    "    return (1 - (x / sigma - mu / sigma)**2) / sigma\n",
    "\n",
    "def _check(f, g):\n",
    "    args = np.random.uniform(size=(3, 1000)) * 0.999 + 0.001\n",
    "    return np.allclose(f(*args), g(*args))\n",
    "\n",
    "def _grad(f, args, ind, eps=1e-7):\n",
    "    a = np.array(args)\n",
    "    b = a.copy()\n",
    "    a[ind] -= eps\n",
    "    b[ind] += eps\n",
    "    return f(*b) / (2 * eps) - f(*a) / (2 * eps)\n",
    "\n",
    "assert _check(lambda *args: np.exp(-nll(*args)), pdf), \"nll or pdf is incorrect\"\n",
    "assert _check(lambda *args: _grad(nll, args, 0), d_nll_x), \"d_nll_x is incorrect\"\n",
    "assert _check(lambda *args: _grad(nll, args, 1), d_nll_mu), \"d_nll_mu is incorrect\"\n",
    "assert _check(lambda *args: _grad(nll, args, 2), d_nll_sigma), \"d_nll_sigma is incorrect\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare transformations to and from modified parameter space.\n",
    "# we transform the space of (µ, σ) into (µ, ln σ) to make the loss convex.\n",
    "\n",
    "def to_modified(typical):\n",
    "    params = typical.copy()\n",
    "    params[1] = np.log(params[1])\n",
    "    return params\n",
    "\n",
    "def to_typical(params):\n",
    "    typical = params.copy()\n",
    "    typical[1] = np.exp(typical[1])\n",
    "    return typical\n",
    "\n",
    "def grad(x, params):\n",
    "    mu, sigma = to_typical(params)\n",
    "    return np.array([\n",
    "        d_nll_mu(x, mu, sigma),\n",
    "        d_nll_sigma(x, mu, sigma) * sigma,  # chain rule\n",
    "    ], float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# configure hyperparameters and the problem to solve.\n",
    "\n",
    "sgd_rate = 1e-3\n",
    "gamma = 1e-2  # γ\n",
    "delta_time = 1e-4  # δt\n",
    "\n",
    "max_iters = 100_000\n",
    "\n",
    "init_typical = np.array([0, 1], float)\n",
    "target_typical = np.array([10, 1], float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged in 82826 iterations\n"
     ]
    }
   ],
   "source": [
    "# perform TANGO\n",
    "params = to_modified(init_typical)\n",
    "xy_tango = []\n",
    "\n",
    "vk = np.zeros_like(params)\n",
    "\n",
    "for i in range(max_iters):\n",
    "    # Definition 1\n",
    "    yk = np.random.normal(*target_typical)\n",
    "    yk_tilde = np.random.normal(*to_typical(params))\n",
    "    \n",
    "    # Equation (2)\n",
    "    gk = grad(yk, params)\n",
    "    gk_tilde = grad(yk_tilde, params)\n",
    "    \n",
    "    # Equation (3)\n",
    "    m1dt = 1 - delta_time  # should be previous delta_time if non-constant\n",
    "    rescale = vk @ gk_tilde\n",
    "    vk = m1dt * vk + gamma * gk - gamma * m1dt * rescale * gk_tilde\n",
    "    \n",
    "    # Equation (4)\n",
    "    params -= delta_time * vk  # should be current delta_time if non-constant\n",
    "    \n",
    "    current_typical = to_typical(params)\n",
    "    xy_tango.append(current_typical)\n",
    "    \n",
    "    # early stopping\n",
    "    if np.allclose(current_typical, target_typical, atol=0.01):\n",
    "        print(f\"converged in {i} iterations\")\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged in 57231 iterations\n"
     ]
    }
   ],
   "source": [
    "# compare to plain stochastic gradient descent.\n",
    "params = to_modified(init_typical)\n",
    "xy_sgd = []\n",
    "\n",
    "for i in range(max_iters):\n",
    "    yk = np.random.normal(*target_typical)\n",
    "    \n",
    "    gk = grad(yk, params)\n",
    "    \n",
    "    params -= sgd_rate * gk\n",
    "    \n",
    "    current_typical = to_typical(params)\n",
    "    xy_sgd.append(current_typical)\n",
    "    \n",
    "    # early stopping\n",
    "    if np.allclose(current_typical, target_typical, atol=0.01):\n",
    "        print(f\"converged in {i} iterations\")\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged in 72003 iterations\n"
     ]
    }
   ],
   "source": [
    "# compare to averaged stochastic gradient descent.\n",
    "params = to_modified(init_typical)\n",
    "xy_asgd = []\n",
    "\n",
    "vk = np.zeros_like(params)\n",
    "\n",
    "for i in range(max_iters):\n",
    "    yk = np.random.normal(*target_typical)\n",
    "    \n",
    "    gk = grad(yk, vk)\n",
    "    \n",
    "    vk -= gamma * gk\n",
    "    \n",
    "    params = (1 - delta_time) * params + delta_time * vk\n",
    "    \n",
    "    current_typical = to_typical(params)\n",
    "    xy_asgd.append(current_typical)\n",
    "    \n",
    "    # early stopping\n",
    "    if np.allclose(current_typical, target_typical, atol=0.01):\n",
    "        print(f\"converged in {i} iterations\")\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAFlCAYAAAB/dUv3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAOwwAADsMBx2+oZAAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3wU1drA8d+mQAgQSqjSCZ2QLElooRfpBAvSe1X0AmKBFwQE5Qp4FTsXkCJFEJUuRVAUlSst0lto0gRCgADp5bx/nGRgSSgJSWaTfb757Cczs7MzzwaSefbMOc+xKKUUQgghhBBJnMwOQAghhBD2RZIDIYQQQtiQ5EAIIYQQNiQ5EEIIIYQNSQ6EEEIIYUOSAyGEEELYkORACCGEEDYkORBCCCGEDUkOhBBCCGFDkgMhhBBC2DAtOVi5ciUtW7akQIECWCwW4uPjbZ4/ceIEzZs3J0+ePJQvX5758+ebFKkQQgjhWExLDiIjI2nRogVjx45N8VxcXBwdOnSgSJEi7N69mwkTJjBs2DB++uknEyIVQgghHIvF7ImXfvnlF5o3b05cXBwuLi4ArF27lq5duxIaGkr+/PkB6Nu3L7du3WL16tVmhiuEEELkeC5mB5CaXbt2UadOHSMxAGjZsmWqrQzJ4uLibG5NJCYmcufOHfLnz4/FYsnUeIUQQgh7oJQiOjqaggUL4uSU/psDdpkcXL16lWLFitlsK1q0KKGhoQ98zdSpU5k8eXJmhyaEEELYvbCwMAoXLpzu19tlcpCeOx3jx49nzJgxxnpkZCRFihQhLCyMPHnyZGR4QgghhF2KiorC09MTNze3JzqOXSYHxYsX59ixYzbbQkNDKVq06ANf4+rqiqura4rtefLkkeRACCGEQ3nS2+l2Weegbt267Nmzhzt37hjbfv75Z+rVq2diVEIIIYRjMK3l4Pr165w7d46TJ08CsH//fpydnalUqRJt27alVKlSDBw4kEmTJrFz506WLVvGxo0bzQpXCGEnlFLExsaaHYYQpsqVK1emdrY3LTlYu3YtAwYMMNYDAgIA2LZtG82aNeOHH35g2LBh+Pv7U7x4cWbNmkXLli3NClcIYQfi4uI4c+YMCQkJZocihKmcnZ2pUKFCqrfTM4LpdQ4yS1RUFO7u7kRGRkqfAyFyAKUU58+fJy4ujqeeeuqJhmkJkZ0lJiZy6dIlXF1dKVOmjE0LQkZd++yyQ6IQQtwvISGBiIgISpcuLQm/cHjFihXjwoULJCYm4uzsnOHHl9RbCJEtJN9KyKxmVCGyk+Tfg/vnJcookhwIIbIVqXgqROb/HkhyIGyc4xw72IEFC554mh2OEHbLarVitVqpUaMGLi4uxnqvXr2MfX788UcsFgtLliyxee3bb7+NxWLht99+M7Z99tln9O/f31i/cuUKAwcOpGLFivj6+uLj48NLL71EWFiYsc/69eupU6cOVatWpWLFigwdOpTw8PDMe9PCYUhyIAwWLJSjHA1pCMB1rmNJ+jrIQZOjE8K+7Nu3j3379rFhwwYKFixorC9dutTYZ/78+TRr1ox58+aleH358uVtqrreKzIykiZNmlC+fHlCQkLYv38/e/bsoWLFily8eBGATZs2MWzYMObMmcPx48c5ceIEuXPnpkOHDumqMivEvSQ5EABEEvnQ533wIYaYLIpGiOzv+vXrbNy4ka+//ppDhw5x6tQpm+efe+45oqOjWbVqVYrXfv311xQqVIiJEycanc1y5crFG2+8gY+PDwDvvvsu48ePp3bt2gC4uLjwwQcfcPbsWX7++edMfncip5PRCgLAuIUwilFUpjKd6UxpStvs44YbCvlEIuxEdDTcd8HNUF5e8AT16ZcuXUrbtm0pWbIkPXr0YMGCBbz77rvG8xaLhWnTpjFy5EiCgoJsXhscHPzIirDBwcF8+umnNtty5cqFv78/wcHBUhdGPBFJDgSJJBJNNAAzmWlsT04E/uZvylMe0Lcekl+TvCyEKU6dAm/vzDv+oUNQs2a6Xz5//nymTZsGwKBBg+jYsSNTpkyxqc/QunVrSpUqxfz581O8/t4OZ9988w3vvfce4eHhTJw40Sggl1qnNLmlIDKCJAeCC1wAYCxjU32+HOXIT35uc9vY5oSTtCIIc3l56Qt4Zh4/nYKDgzl06BBDhw41LuBXr15l8+bNtGvXzmbf6dOn07lzZ0aMGGFsq127tk3C0K1bN7p160b//v2JiIgAwM/Pjx07dmC1Wo39YmNjCQ4OZuTIkemOXQiQ5EAAdagDQBBBD9znJjeZylQmMjGrwhLi4dzcnuiTfWaaN28er732mtFyAPDpp58yf/78FMmBv78/jRo1YtasWTRt2hSAnj17MmPGDKZOncrYsWONfgeRkXf7Bo0bN47BgwcTGBiI1WolPj6e1157jbJly9KiRYsseJciJ5MOiYJ+9AOgPvUfuI8TTkxgAgpFZzoDcIc7D9xfCEcVHR3NsmXLbIY0gv70v2HDBq5du5biNVOnTjVGIQDkzZuX7du3ExISQqVKlbBarTRo0IAiRYrwzDPPANC+fXtmzZrFoEGDqFq1KpUrVyYqKooNGzZILQjxxGRuBZFmCoVTUl4ptxZEVomJieH06dNUrFiR3Llzmx2OEKZ60O9DRl37pOVApNm9HRETSeQSl1jPehMjEkIIkZEkORB8wicc5WiaXpM89NEZZ0pRik504n/8LzPCE0IIkcUkOXBw8cQzilEcIm29vt/jvRTbAgnkG77BGWf2sjejQhRCCJHFJDlwcFe4gkJRkpJpet0QhhBFFO64U5rStEQXXOlOdxJJJIAA6Y8ghBDZlAxldHB/8icA+ciX5te64UYEEcb6/UWRnHAinnicyfi5xoUQQmQeaTlwcOc5D0B1qj/xsRSKcYyjO92NbS6SfwohRLYjf7kd3Ku8CkBuMmZo2FSmAtCVrjzHc4BOGqTUshBCZB/ScuDghjIUd9wz/LjP8qxRjvkt3srw4wthD1auXIm/vz9Wq5Xq1avTsmVLEhMTAfjoo4+4evXqEx2/f//+fPbZZxkRKpAypv/+97/MnDnzIa/IOmfPnqVIkSLGutVqJSoq6omOmdZ/A4vFwp079lHc7ebNm8yYMcO080ty4OBmM9um30BGGsYwAP7Nv3mGZzjAgUw5jxBmuHz5Mi+++CIrV65k3759HD16lPfff9+oTpgRyUFGuz+mF198kVdffTVTzhUfH/9Er9+3b98TF7Czx3+DxyXJgcixylOeBSwAYA1r8MWX13nd5KhETnKVqxx+wNc//GOz7xGOPHDfOOKM/cII4yqPvqD8888/uLi44OnpaWzz8/PDYrEwZcoULl26RJcuXbBarezbt487d+4wcOBAvL298fb2ZvLkycbrLl68yPPPP0+tWrXw8fFhwoQJd+M+coRWrVpRpUoVunTpQlycjvWnn36iQYMG1K5dG29vbxYsWGC85ssvv6RGjRpYrVZq1arFzp07U43p7bff5vXX7/5OTps2DW9vb3x9falfv77NXA73xlOvXj28vb3p3bs39evXZ/16XQStf//+jBgxgnbt2uHr6wtA7969CQgIwMfHh44dO9pcrD///HMqVapE48aNmTdvns157v0UHxISQocOHahTpw6+vr588cUXNvtNnz6devXqUaFCBb766iuAVN/v/VauXEm1atUIDAy0mU77YeeMioqiW7du1KhRA19fX1q3bm28ZsGCBfj6+uLr60tAQABnz54FYPPmzTRq1Ah/f3/q1avH9u3bAfjll1+wWq0MHz4cX19fvL29CQ4OBnTidvPmTaxWKwEBASliz3Qqh4qMjFSAioyMNDsUhxav4hWpfC1QC8wOTWQz0dHR6siRIyo6OtrYNklNSvX/Fwo1Qo2web2bcnvgvhfVRWO/j9RHapKa9Mh4EhIS1HPPPacKFSqknnnmGTVjxgx14cIF4/ly5cqpgwcPGutvvvmm6tWrl0pISFB37txRVqtVrVixQimlVLNmzdSMGTOMfa9evaqUUqpfv36qQYMGKjIyUsXHx6vAwEC1fPlypZRS169fV/Hx8UoppcLCwlS5cuXUpUuXlFJKeXh4qIsX9XuKjY1Vt2/fTjWmSZMmqddee00ppdTChQtV/fr1VXh4eIrj38vPz08tXrxYKaXUnj17lJOTk1q3bp0Rb+3atY3zKaVUaGiosfzee++pl19+WSml1P79+1XJkiXV5cuXlVJKvfTSS8rT09PYF1C3b99W8fHxKiAgQB09elQppVRERISqVauW2rt3r7HfRx99pJRS6siRIypfvnwqLi4u1fd7rytXrqjChQurY8eOKaWUmj59+mOdc+XKlerpp582jhMWFqaUUmrbtm3Ky8vL+DeIiIhQERER6tSpU6pBgwbGzzUkJEQ99dRTKjY2Vm3btk25uLio3bt3K6WUmjVrlmrbtq1SSqkzZ87Y/Dzul9rvg1IZd+2TDokOrhe96EMf2tI2U47vjLNR7+BDPuQ1XgNI8alOiPQYznBe4IVUnytMYZv1vex9YO2NohQ1lnvTmwQSHnluJycnvv/+e44dO8avv/7Kpk2bmDp1Knv27KFSpUop9t+6dSsff/wxTk5O5M2bl759+7J161batWvHjh072LJly914it6N57nnnjOa1+vWrUtISAgAYWFhDBo0iBMnTuDi4sK1a9c4fPgwJUuWpEWLFvTt25dOnTrRrl07qlSp8sj3s379el566SU8PDwAKFSoUIp9bt26xaFDh+jZsyegZ5T08fGx2adr167ky3d3aPTSpUtZvHgxMTExREVFUaJECUB/au7QoQPFixcHYOjQoaxYsSLFOY8fP87hw4fp3v3uKKjbt29z5MgR/Pz8AIxJrqpXr46LiwuXL1+mdOnSD32/f/75J35+flStWtU4/5gxYx55zsDAQI4dO8bw4cNp2rQp7du3B+CHH36gb9++lCypa8a4u+u+XJs2beLkyZM0adLE5vznz+uRYlWrVjVaBho0aMB//vOfh8adVSQ5cGAKxdd8jT/+mZYc3Gs0o+lNb9rQhnGMYxe7WMWqTD+vyLmKJX09jhrUeKz9kkuDP65q1apRrVo1hg0bRtu2bVm7di2jR49OsZ9KZY67x5k90c3NzVh2dnY27uW/+OKLdOrUie+//x6LxYKfnx/R0dGAbi7fu3cvv/zyC+3bt+fdd9+1udClJrX4UtvHYrE8NO57E4Pff/+dzz77jB07dlC0aFHWrl3LlClTHvt8yfsVKVIk1dsCyR70M3rUcdN7ziNHjvDzzz+zdetW3nzzTfbt2/fA4ymlaNu2LYsWLUrx3Llz59IVe1aQPgcOLArdE7gQKT8hZJZiFONjPgZgNav5lV+z7NxCZKSLFy/yxx9/GOs3btzgzJkzeHl5AeDh4UF4eLjx/NNPP82cOXNQShEREcGSJUto1aoV+fLlIzAw0GbUQGho6CPPf+PGDcqVK4fFYmH79u3s378f0B0BT506RUBAAK+//jpdunRh165dqcZ0r06dOjFr1ixu3boF6A5xCQm2LSgFChSgRo0aLFu2DIC//vqLgwcPPjRGDw8PChcuTGxsLLNnzzaea968ORs2bDD6INzf5yBZ1apVcXd3t7m4hoSEcP369Yf+fB71fhs0aMBff/3FiRMnAN1P43HOeeHCBSwWC0FBQfznP/9BKcX58+fp1KkTixYt4vLlywBERkYSGRlJ69at2bRpE4cO3S1Rn/zv8ajYIyMjTUsWJDlwYDe4AUBlKmfpeZvQhEEMAqAZzaQGgsiW4uPjmTJlClWqVMFqtdK4cWP69etH586dARgxYgQDBgwwOsNNmDABZ2dnatWqRb169QgKCqJLly4ALFmyhB07dlCzZk18fX0fa/jitGnTeOONN6hfvz4LFy6kXr16ACQkJDBgwAC8vb2xWq3s3bvXaMm4P6Z7Jd+GaNCgAVarlfbt2xMTE5PivIsWLWLmzJn4+/sza9YsfH19KVCgQKoxtmvXjkqVKlGtWjXatGmD1Wo1nvPx8WHcuHEEBgbSuHFjSpUqleoxXFxcWLduHStWrMDHx4eaNWsyZMiQxxrm+LD3W6xYMebMmUOnTp0IDAzExeVuQ/rDznnw4EECAwPx8fHBz8+PPn364OPjQ9OmTRk3bhytW7fG19eXpk2bEhoaSuXKlVmyZAmDBw/G19eX6tWr8/HHHz8y9sKFC9OrVy9q1aplSodEi3rctp1sJqPmtM7J9rMfK1YOc/ixm1wz0r1JwTKW2VRWFOJ+D5q/XmStiIgI3N3dsVgsHDlyhGbNmnH8+PFU+yiIzPOg34eMuvZJy4EDO85xAPKT35TzJ5DARCYC0IMeUixJiGzgjz/+oHbt2vj4+NC9e3fmzp0riUEOJC0HDmwVq3iO54ghhlzkMi2Oe1sQZCZH8SDSciDEXdJyIDJNQxqyiU2mJgZgmxAMZKCJkQghhABJDhxaMYrRhjZmhwHAQXSP5wUsoCpVTY5GCCEcmyQHwi54400PegBwghPEYx9jfYUQwhFJcuDA9rOf3ew2OwzD13zNSU4C4IorFiy8wismRyWEEI5HkgMHNpOZTGKS2WHY8MLLZv1zPic06QtgLnON+gxC2IPw8HDc3d0ZPHiw2aGk28KFC42aC/c7ePAgLVq0MCYGqlOnjk1Bny1bttCkSRMqVqxIQEAAdevWZc6cOcbz5cuXp1q1avj6+lK5cmU6d+7Mjh07Mv09iScjyYEDCyecAqRevMRMk5lss55cIteChaEMTVEzXwgzLVu2DH9/f77//nsiIjJ2+nN7KKXbs2dPRo0axf79+zl06BArV66kWDFdsvrHH3+kX79+TJs2jdOnT7Nnzx6+/fZbTp48aXOM7777jv379xMSEsLAgQNp3749O3fuNOPtiMckyYEDu8UtPPAwO4wUJjCBWGIfOqzxTd7MwoiEeLB58+YxduxYGjZsaDNxUJUqVdi7d6+xvmDBAl54QU8SdfnyZbp27UrdunXx8fFh4sSJxn7ly5dn6tSpNG/enH79+nH58mWaN2+Ov78/NWvWZMSIEUYd/9jYWIYOHUqVKlVo1KgRr7zyik0LwOLFi6lXrx5+fn40bdrU+MQfGxvLsGHDqFKlCi1atHhoOd9z587ZTGJUpkwZIzmYMmUKEydOJDAw0Hi+XLlyzJgx44HH69y5M8OHD7ebCYZE6iQ5cGD22nJgwYIrrgBEE81qVqfY533e5yY3szo0YWeuXoXDh+H4cdvtx4/r7Ull+wG4fVtvO3wY4uLubv/7b73t77/Tfv6DBw9y6dIl2rZty8CBA23mB+jfvz8LFiww1hcsWMCAAQMA6NevH6+88gq7du0iODiYXbt2sWrV3UnIzp07x88//8zSpUspWLAg69atY+/evRw4cIDTp0/z/fffAzB79mzOnTvHkSNH2Lp1K3v27DGO8ccff7B8+XK2b99OcHAw7777rjF74ezZszlz5gyHDx9m/fr1D00OJk6cSJMmTWjZsiXjx4/nr7/+Mp4LDg42yjanRZ06dTh8+HCaXyeyjiQHDiyccLtsObhXbnLTmc6oe76SfcAHJkYm7MEXX4C3NzRvbru9eXO9/Ysv7m775Re9zdsb7p3X6JVX9LZX0tH3dd68efTt2xdnZ2c6depESEgIx5MylX79+rFixQpiY2M5deoUp06dok2bNkRERPDzzz8zYsQIrFYrAQEBnDx5kmPHjhnHHTBggDHzYWJiImPGjMHX15fatWuzZ88eY56Abdu20adPH1xcXHBzc6NHjx7GMdasWcP+/fupV68eVquVf/3rX4SGhhIbG8u2bdvo168frq6uuLu707t37we+x9dee41Tp04xePBgbty4QePGjfnmm2+M5++dobFXr15YrVZKlChBZGTkA4+ZQ2vv5SgyZbMDu8Utu2w5eBSFwoKFd3mXPORhHOPMDkmYZPhweOEFcLnvL9m2bRAfD0WL3t3WrBkk96O7d/tnn8G0aXDPTMOPJTY2lqVLl5IrVy5jlsLo6Gjmz5/P9OnTKVWqFH5+fqxdu5b9+/fTu3dvnJ2dSUxMxGKxsHv3blxdXVM99r3THn/44YeEhYWxc+dO3NzcGD16tDE1c/IUyqlRSjFw4EBjiuT7n0uL4sWL06NHD3r06EHZsmVZunQp3bp1o3bt2uzcudOYUGnp0qWAThgSExMfeLzdu3fj7e2dphhE1pKWAwdmxUp5ypsdRrpsZzsA4xlPGGEmRyPMUqwY1KwJVe+rm1W1qt6edGscgPz59baaNeHea3K5cnpbuXJpO/eaNWvw8vLi4sWLnD17lrNnz/L777+zaNEioyPhwIEDmT9/PosWLaJ///5JceSncePGTJs2zTjWxYsXuXDhQqrnuXHjBiVKlMDNzY0rV67w7bffGs81b96cJUuWEB8fT3R0tM0n+uQphM+fPw/oFojdu/XQ5ZYtW7J48WLi4+OJiori66+/fuD7XLVqFXFJ92Hi4+M5cOCAMS31hAkTmDJlCn/++aex/8NaDJJ/brNmzTJmihT2SVoOHNhmNpsdQro1prGxPJ7x/Jf/mhiNcETz58837uEnq1WrFkWLFmXDhg0EBQXRuXNnXnrpJSpXrkz16tWN/ZYuXcro0aOpVasWAHnz5mX27Nk2Hf+SjRgxghdeeAGr1UqpUqVo1aqV8dyLL77I/v37qVmzJmXKlMHf39+4ODdp0oR///vfdO7cmYSEBGJjY+nYsSN16tRh6NChHDhwgBo1alC2bFmaNm3KmTNnUn2fK1euZOzYseTOnZuEhATq1q3L5Ml6RFHbtm2ZN28eb7zxBpcuXaJo0aLkypWLzz77DHd3d+MYXbp0IXfu3ERERFCjRg02bNhA/fr10/mTF1lBJl4S2dZ1ruOJp7HegAb8m3/TjGbmBSUyjUy8lLrbt2+TP39+YmJiCAoK4oUXXsjWNRfE48nsiZek5UBkW/fXO/gf/6M5umdaKKEUoYgZYQmRpVq1akVsbCxRUVG0atXKuH0hxJOQ5MBBhRJKEEHMYQ61qGV2OOl2jGNUo1qK7U1owhGOmBCREFlLigmJzCAdEh3ULW7xJ38SS6zZoTyRqlQ1hjg2pCED0OPIj3KUM6R+D1UIIcTDSXLgoKKIAiAPOac/xu/8znzm8zmfA1CRinzHdyZHJYQQ2Y9dJwc3b95k0KBBlChRgnz58hEYGMj27dvNDitHyInJQbLhDDeWX+AFYogxMRohhMh+7Do5GD16NLt372b16tXs37+funXr0rFjR27ckFn5nlROTg4ArnK3bq4bbiZGIoQQ2Y9dJwc7d+5k4MCB1K9fHy8vL9555x1u375tlCcV6ReJHgudU5ODohS1KbX8sEmchBBC2LLr5KBBgwasWbOGa9eukZCQwPz583nqqadSLbsZFxdHVFSUzUM8WHLLgTvuj9gze0ue/tkJJxaz2ORoRE5itVqxWq3UqFEDFxcXY/3+wkj24Pr165k2C+LJkycpUaJEphz7YUqXLm3MRzFgwAB27NjxRMdbuXKlzcRVjs6uk4NPP/2UIkWKULRoUXLnzs17773HDz/8YFN3PNnUqVNxd3c3Hp6enqkcUSRrQAPWsc6Y/TCnmsjdqXD70pfRSMlWkTH27dvHvn372LBhAwULFjTWk+cXuFdyOWWzPElykFWxJyQkpPu1CxYssJk2Oj0kObBl18nBxx9/TEhICFu2bGH37t306NGDoKAgwsJS1tIfP348kZGRxiO1fcRdJShBRzqaHUaWuM1t2tMegJnMxIKFRB48KYwQT2rr1q34+/vzyiuv0KBBA9auXUvv3r3573/vlvkeNWoU7777LqAncXrjjTeoW7cuVquVnj17cuvWrRTHTf6UPn78ePz8/KhSpQpbt241nu/evTt16tTBx8eHoKAgQpOmn3zxxRcJCwvDarUaUyw3atSITZs2Ga995plnWLJkCQC9e/dm1KhRtG3bFj8/v4ce+2Fu3rzJs88+S7Vq1WjVqhW9e/dm7NixAHz55Zd06NCB3r17ExAQQHBwMDNmzDB+BvXq1bO5WP/666/UqlWLevXqMWrUKJvJo+59L7du3WLw4MHUrVsXHx8fXn75ZSPBadSoEWPHjqVx48Z4eXkxYsQIANauXcuGDRuYOnUqVqvVZqpth6XsVGRkpHJ1dVW//vqrzfbKlSurTz755LFeD6jIyMjMClFkM5VUJUUqXyJ7iI6OVkeOHFHR0dF3N15RSp1NWr6klDqftHw+aV0lPX8lafm0Uupa0nKIUupG0vIxpdStpOXDSqk7jx/XmTNnlKenp822LVu2KCcnJ7Vjxw5jW69evdSsWbOM9ZEjR6p33nlHKaXU5MmT1XvvvWc8N3HiRDV69OgU5woJCVGAWr9+vVJKqXXr1ilvb2/j+dDQUGP5nXfeUSNHjjReV7x4cZtjNWzYUG3cuNFY79y5s1q8eLERa0BAgLpz5+4PIi3HTjZixAg1bNgwpZRS165dU2XLllVjxoxRSik1d+5clS9fPnXy5Elj/6tXrxrLv/32m/L19VVKKRUVFaVKlCihfvvtN6WUUkuXLlWAOnr0aIr3MmDAAPX1118rpZRKTExU/fr1M64ZDRs2VN26dVPx8fEqIiJClSlTRu3Zs8d4z/f++9i7VH8fVMZd++y2QmJcXBxxcXE4OzvbbHdycnroVKDi8ZzlLKGEUoc6ZoeSZU5wgggiyE9+m+0xxJAbqdWfLX0B7AXWAdOAG8AiYBxQCPgYeAXwB94G+gAvACOBzsAYoC/QHJgNdEradwvQ6MlCq169Og0aNHisfVevXk1kZCTLly8HdEtC1funmkxSoEABOnToAOh+WSEhIcZzixYtYunSpUY55dQmcnpc3bp1I2/evE907G3btjF37lwAPD09CQoKsnm+adOmxgyPAHv27OG9997j+vXruLi4cOjQIRISEjhy5AgFCxakUSP9j9KzZ0+GDRuW6jnXrFnDnj17mD59OqDnGvDw8DCe7969O87Ozri7u+Pr60tISAj+/v6P+VNxHHabHHh4eNCwYUNGjx7NJ598gqenJwsXLuTMmTO0bt3a7PCyva/4iuUs5yhHzQ4ly1iwkI98KE7zTW0AACAASURBVBSnOY0X+o+SG24kkogFi8kRijQbDiT3PR4LJN+2/jeQ/LniMzAG5SwGkq8Ta8CYfmMb8FTS8l4gjdM3p+b+vlEuLi4299Wjo6ONZaUUc+bMoUmTJo88rpvb3aG5zs7ORpP5L7/8wpw5c/j9998pUqQIK1euZMaMGQ88zsPiuT/+tB773vdlsTz49+rec0RHR9O1a1d+++03rFYr169fx9PTk7i4OJtbCPceOzWJiYmsX7+esmXLpvr8g35+wpZd9zn45ptvqFixIkFBQVitVjZv3syqVatspj4V6RNFVI4dxvg4KlIRhaIQhQA9mkFkQ8W4eyEvCSR/mC2dtE7S88WSliuAMZFnJaBg0nJVMBqUagB3PzBnGC8vL2MehGvXrrFx40bjuaCgID744ANjlNWdO3c4fPhwmo5/48YNPDw8KFy4MLGxscyZM8d4zsPDg4iICJtk4N54Tp069dDe/g879sM0b96chQsXArpT5Lp16x64b2RkJPHx8UaLxKeffmo8V7NmTcLDw40Yly9fTkRERKrHCQoKYtq0acZ7DQsL4+TJk4+M1cPDg/Dw8Md6X47Arv8ilipVimXLlnHlyhXu3LnD7t27ad++vdlh5QiOnhwkO8tZAOojc8uLzPXSSy9x7tw5fHx8GDJkCPXr3/0/N378eGrUqGF0+AsMDOTgwYNpOn7Hjh0pV64c1apVo02bNkZHQoBixYrxwgsv4O3tbXRI/L//+z9++OEH/P39mTRpEnXr1k3XsR9m8uTJXLx4kZo1a9KvXz8aNmxIgQIFUt23cOHCTJw4kYCAAJo0aYK7+91h1m5ubixbtoxhw4ZRr149Dh48SKlSpVI9zqeffkpiYiJWq5VatWrRunVrzp0798hY+/bty1dffUXt2rWlQyJgUQ9qm8nmMmpO65xqCEM4wxm2svXROzuARBKJIoq8mfGRUWSIB81fL+xXXFwciYmJ5M6dm/DwcAIDA/n8889p1qyZ2aFlew/6fcioa5/d9jkQmSuKqBxfACktnLnb8TWGGHKRy8RohMgZrl27RseOHUlISCA6Opo+ffpIYpBNSHLgoOS2woPlJjev8Rr/IXMqygnhKEqWLMnevXvNDkOkg133ORCZpxSlqEQls8OwG/cXRfqAD2T0ghDCYUnLgYP6hE/MDsGuWLCgUGxjG//iXxxG9xSPIEL6IdiZHNpNSog0yezfA0kOhLhHc5pziEOc4Qy1qGXURRDmc3V1xcnJiWvXrlGkSJGHjp8XIidTSnHt2jWcnJxwdc2c+XEkORAiFRWoQAR6HLUFC6GEUsSomCPM4OTkRJkyZTh//jy3b982OxwhTJX8++DklDm9AyQ5cFD96U972tOVrmaHYrf+4R9KJlXSKUpRTnCCylQ2OSrH5u7uTuXKlYmLizM7FCFMldySllkkOXBQW9mKDz5mh2HXSlCC4xynKrrGfRWq4IorN7gh/RBM5OTkJHUOhMhkMlrBQUUSKUMZH0MVqtj0OYgjjnzkw4KF4Qw3MTIhhMg8khw4KCmClDaJJNKNbjbbZjGLAQwwKSIhhMg8khw4IIUimmhpOUgDCxaWs5xEEtnEJmP7QhZiwcKzPGtidEIIkbEkOXBA0eipWSU5SDsLFtrQBoUihhhj+2pWS9EkIUSOIcmBA4pCTwsrycGTyUUualDDZtuf/GlSNEIIkXEkOXBA7riznOXUopbZoWR7hzmMQjGGMQCc57zJEQkhxJOTKZuFyCDJtxXiibeZ5VEIIbJKRl37pOVAiAzmIuVDhBDZnCQHDug2twkmmDikylxGundmRwsWRjFKbjMIIbIlSQ4cUDDB+OPPda6bHUqOkjyzY7KP+ZiylCWUUBOjEkKItJPkwAFFEgnIaIXMcoMbNuvFKMYd7pgUjRBCpJ0kBw5IhjJmroIURKEII8zYlp/8/MVfJkYlhBCPT5IDBxRFFC644ErmzAMutMIU5h/+Mdb98DMxGiGEeHySHDigKKKk1SCLlKCEUZFSCCGyCxlz5YAkOchaubk7vbAFC1/yJbnIRX7y8wzPmBiZEEKkTpIDB+SGG5WoZHYYDuU2t8lPfgAGM9jY/iVfMohBZoUlhBCpkgqJQmSh1CZniiVW+n8IITKEVEgUIhuKI46TnLSph5CLXCZGJIQQKUlyIEQWcsEFL7wAiCDC2G5J+koeZiqEEGaS5MABfc7n/Jt/mx2Gw3PHnR3sSLFNCCHMJsmBA9rBDnay0+wwBNCABoxiFC/wgrHNgoW61DUxKiGEo5PkwAHJUEb7MpOZrGAF7/COsW03u1nOchOjEkI4MkkOHJAkB/bpLd5CoZjABAB60IMEEkyOSgjhiCQ5cECSHNi3KUzBJakESW1qmxyNEMIRSXLggCKJlOTAzt3iFgAHOZhqbQQhhMhMkhw4oCiipFe8nctDHptbChYsvM/7JkYkhHAkkhw4oIlMpDOdzQ5DPIITTnzER8b6m7xpUxtBCCEyi5RPFiIbuP/WQiKJcrtBCJGClE8WwoEkkkgpShnrM5lpYjRCiJxOkgMHFEKINE9nMxYsXOACl7kMwGu8RiKJJkclhMipJDlwMIkkUoUqbGCD2aGIdChOcT7mYwCcceY4x02OSAiRE0ly4GCiiQbADTeTIxHpNYIRxnI1qmHBIsWShBAZSpIDB5M865/UOcje7r+l4IKLdFAUQmQYSQ4cTHLLgSQH2ZsFC38mfd2rEIVMikgIkZO4mB2AyFrJyUFucpsciXhS9agHQBxxnOEMVajCTW6aHJUQIieQ5MDBxBADSJ+DnMQFFypTmZ/4iXzkI4wwFIoiFDE7NCFENiXJgYNJIIHCFJbbCjlQC1rY9Dv4gz8IJNDEiIQQ2ZX0OXAwtahFGGF44WV2KCITeOJpLDekITWoISMZhBBpZvfJQXBwMC1btsTd3Z1ChQrRtWtXs0MSwm5d4xoKhQceABzlKC64oMiRVdKFEJnErpODo0eP0qJFCxo1asTu3bvZsWMH3bt3NzssIexeOOGsZa2x3o52WLBIVUUhxGOx64mXnn/+eTw8PFiwYEGaXysTL6VuL3uZxSxmMxtnnM0OR2SyG9ygMIWN9bnMZTCDTYxICJGZcvzESwkJCWzatIkKFSrQrFkzihcvztNPP82BAwdS3T8uLo6oqCibh0gphBDmMx8n+/2nFxmoEIU4xjFjvTjFCSfcxIiEENmB3V4hQkNDiYyM5P3336dHjx5s3LiRMmXK0LJlS8LDU/5xmzp1Ku7u7sbD09MzlaOKaKLJTW6ppudAqlKVeOIJI4wggihIQbNDEkLYObtNDhIT9b3RLl26MGzYMPz8/Jg9ezYWi4W1a9em2H/8+PFERkYaj7CwsKwOOVuIIUZqHDggZ5xtbi9sYIN0UhRCPJDd1jkoUqQIzs7OVK1a1djm6upKxYoVOX/+fIr9XV1dcXV1zcoQs6VooiU5cGAXuEBpStOBDoCeo0FakYQQ97PbloNcuXJRu3ZtTp48aWyLj4/n7NmzlC1b1sTIsrcYYqR0sgMrRSl+5Edj3QknY1ZHaUkQQiSz25YDgFdffZVBgwbRvHlz6tSpwyeffAJAUFCQyZFlX9JyIJ7maRJJtOmU6pL0p0ASBCEE2Hly0LNnT0JDQ/m///s/bty4QUBAAFu3bsXDw8Ps0LKtVrSiDGXMDkOYzIIFhWIVqzjFKd7gDUAnB3KbQQhh13UOnoTUORDi8W1kI+1pz0hG8hEfmR2OECKdMuraJ8mBEALAaDHIRz5uc1tuQQmRDeX4Ikgic4QTzi1umR2GsEOu6NE+d7iDBQt5yMNlLpsclRDCDJIcOJiXeZk+9DE7DGGHoolmPvNttpWkJLe5bVJEQgizSHLgYJIrJApxPyecGMAAFIoIIviWbwHwwIONbJRJm4RwIOkeraCU4vLly8TFxdlslxoE9i2aaJtKeUKkxh13utDFWG9Pe0CGOgrhKNLcchAWFkb37t1xc3OjdOnSVKhQweYh7JuUTxZpkUgii1hkrO9jn4nRCCGySpqTg3/9619cunSJX3/9lTx58vDDDz+wcOFCqlatynfffZcZMYoMJLcVRFpYsNCHPmxjGwC1qS11EIRwAGm+rbB161Y2bdqEn58fTk5OVKxYkbZt21K0aFEmTZrEs88+mxlxigwiw9NEejSjGYtZbHRmlWJJQuRsaW45iIuLo1ChQgAULVqUixcvAlCpUiUOHjyYsdGJDKdQkhyIdOlNb8YwBtCdFzexyeSIhBCZJc0tB7Vr12bv3r1UqFCBwMBAJk+ezO3bt1m8eLHNDIrCPu1hj9khiGxsIhOZznQA2tEOgGtcwxNPM8MSQmSwNLccTJ06lQIFCgAwbdo0XFxc6NWrFyEhIcydOzfDAxRC2A933Ekk0WbESxGK8DqvmxiVECKjSflkIUS63dvvYC978cPPxGiEEKaXT46IiODs2bOcPn3a5iHs21jGsotdZochcogLXKAylQHwx58QQkyOSAiREdLc5+DAgQMMGjSI4OBgQBdDslgsxveEhIQMD1JknI/4iBrUoC51zQ5F5AClKMVxjuOU9DmjClUAXR9BRjMIkX2lueWgb9++lC1blh07dnDq1CnOnDnD6dOnje/CfimUFEESGc6CBYUiiCBj2//xfyZGJIR4Umnuc5A/f3727duHl5dXZsWUIaTPQUrJicFqVtOZzmaHI3KgeOKN2R2l1LIQWc+0Pgft2rVj1y65Z50dRRMNIC0HItO43HOn0oKFUYwyMRohRHqlueUgPDycXr16UbFiRWrWrImrq6vN8wMHDszQANNLWg5SuspVilOcX/iFpjQ1OxyRQykUXenKd3xnrAshskZGXfvS3CFx9erVbNmyhVy5cuHp6YnFcrfTkcVisZvkQKSU3HIgcyuIhzoLFAHyAt8ArYB8wLtgNASMAP4LXAOGApuBA8BosGyz8O2v37L+/fV0Wt+Jrhu70m5FOwYsGABrgZ3AVOAH4BIwBNgOxAMtgBOAG1AWSEQmlhfCBGn+tRszZgxvv/024eHhnD17ljNnzhgP6ZBo3/KRj7d4izKUMTsUYZaDwLGk5Q/QF3WAlsDypOUGwM+ABX3hPob+S7EWuJm0HIG+mOcGKqIv4h5AvaRjeELLei0BOF/wPCGVQ2hNa/1xJDk3vQScSVpeD6xJWn4L+DRpuT0wI2n5X2BUbF6GTmIA7qTtRyCEeLQ031bw9PRk586dVKpUKbNiyhByW0E4HIW+UOYH/kBfPHsBE4AbwGdJ60WAj4GxgBXoDiwFfAFvdAJRDn2xz4BP7vcOafTCi8lMphe9HvyCG0nvpTA6SSkB1EAnDW2BRkBldALRGnAHfgLqAAOA/yS99lDSe0pz+6gQ2VdGXfvSnByMGzcOi8XC1KlT033SrCDJgcjRooA8wHHga2AysBDdGnAQmIduqv8K+A2IBNoAt9AX0yy+YEYSSV7yGusZ0g9BAQnAj+iEIRadHCwBTgEB6BaOY8AUYCUQClwHqj356YWwR6b1OTh37hzr169n7dq1eHt7p+iQuGjRonQHIzJXQtKXK65SoCa7iQLmoi9+F9Gf8K8D0ehP1xPRn6LLJ+0/KOkB0Pie43hkQaypcMed13md//AfQLcm3OY2+ciX/oNa0H/B2t+zbV3Sdx9034U8Sft5JX3/Fp04BQPfoVtXXkcnTfmQ/g1CJEnzr4KrqyvPPvssAQEBuLm54ezsbPMQ9utnfiY3ubnBDbNDEQ/zD/oTrgKao1sAXNH33k+hL3TLk7b5olsGnIGngGZZH+7jep/3bVoM8pM/807mAiTf+bRCUk4Cg4HVScuR6NsmAC8DI5OW1wCHMy80IbKDNLccLFiwIDPiEFkghhhA6hzYpeSOgW2ArkA7YBzQCSiE/k09D0aDT5esDjDjhBNOAfTMrhOYwDu8k3Unz40eBQHQ957tbwDJjaAzgf5ATaAH8CpQF4gBGegjHIU0ojkQGcpoJ5KnH/keeC1p+Td0pzqABcArScujgVpJyznkTpAHHkYLwru8a3I0SXyA6knLv6ATh3j0bRhP9JDNgsBpdItDMEj5BpGTpbnloEKFCja1DZJZLBZy586Nl5cXffr0oWvXrhkSoMg40UTjggvOyO2fLHceKIPuHzAIfZEpCkkfoOEd7l787XsgUIazYKE3vVnMYrNDucsp6TE7af02us9HeWAjetTHNSAE3bdBqpGLHCbNycFLL73E9OnTadu2LQEBAQDs2bOHzZs3M3LkSM6dO0e/fv0IDw9nyJAhGR6wSD+ZdCkLKfSnSy/0J82y6AtJbfTQwnigSdIDsl+rwI0bcPQoXL0KkZFw/jyEhcHt23cfERGQmAhXrsA//0ClSuDrC+7uEBPDnStPk2/VFgCWJH2pWt7g5wf+/lCjhn5NmTJgdn+m/EDvpOV2wD70X8+t6NESnYEV6FsTz5oRoBAZK83Jwe+//8706dMZPHiwzfZ58+axZs0a1q5di7+/Px999JEkB3YmmmhJDjLbHnR1v5rA8+iqgr2B39EJQi4gOxURvXIFVqyAAwfg4kW4dAnOndPJQTJnZyheHEqWhPz5IV8+/b1YMbBYoGFDKFFCJxOHD0NcHLi6krdUKdSUyXzR+iQv109qNWjfHoKDYd26u+ewWKBAgbvnK10avL31IzxcJyYuLuDhoWMoXFgnFQ0bQt67wyczjBN3R4X8CxietHwA3RoE8CZ6FEWzjD+9EFkhzXUO8ubNy759+6hcubLN9pCQEKxWKxEREZw+fRpvb28iIyMzNNi0kDoHKX2Y9HWBC2aHkrOcQFf6a4NOBIoBH6KHHD5F9msViIyE1athyRL48UfIkwfq1oVSpfSjdGmoWBGqV9cX49xP3oflJCdZzGKmMIUf+ZGnVSu4dg1OnoS//4abN3WSoJReP3AADh2CQoWgfHndQnH9uk5mwsJ0YuHqCuXKQa5cULs2BAZCnTqQPKNs/vx6n4ym0CWlh6BbijoBX6ArSQqRyUyrc1CmTBlmzZrFhx9+aLN91qxZlCmjy/KGhobi6emZ7qBE5hjJSIYxzOwwcoYwdCnfXsCvwCJ0cvBfMGr9lDIntHS5cwd+/x2WLYOVKyEqCtq0gcWLoXNnfSsgE1WiErnIBUBrWrPZsplmRZuRq2gDaNAg7Qe8dAl++gnOnIHoaNi9G8aO1bc7kuXOrVsf/PygUSPd0nD2rH6vNWvqloj0sKD7J4Aekloy6bEfXeVxFVK1Udi9NLccbN26leeff57ChQtjtVqxWCz89ddf3Lhxg++//56WLVvy1Vdfce7cOSZMmJBZcT+StByIDBePnl+gNXABPQfBaXTxHGey19if69fh55/1BfR//4ODB/Wn77p1oXdv6NZN3xbIYvcX50okMeMKdiUkwOnTOmEA3Q8iOFg/du7UtzvuVb8+DBumk6NChZ78/MHoKo3voqtZXgTGP/lhhbiXaeWTAW7evMnSpUs5efIkSimqVKlCz549KViwYLoDyWiSHIgMcwr9abAMUAH4El3jP7uMe1dKX/zWrNEXxuPHYf9+3Uzv768/NTdI+oReurTZ0fIMz7DGmIVJm8MchpCJfZhu3dL9IcqX17dVDhyAr7/Wt1fi46FyZejVC4YO1bdSntQSdHXLwejqls3QM1IK8YRMTQ6yA0kOUlrPeq5znb421V9EqhLQrQPl0B0LSwGfoCc2eoKKv1kiMhL27YNjx/T3tWv1ffpKlXRzuZeXbkJv3jxjPhFngtRG1iSQgFNWN8+EhsKff8Ivv8BXX+m+D506wZAh8PTTT95nQaH7JgxEF1qaiC66VPRhLxLiwbI0OZg/fz69evUid+7czJ8//6H7DhxoH12xJTlIaQhDOMtZtrDF7FDsV3JrwBx0ueIQ9Hj25CqF9uzIEZg1S1/Ebt/WHQmrVdN9B7p2BatVtxZkI7HEkoc8JCbVOR7EIL7kS3OCiY6G776DuXNh+3bdZ+Hpp+G556Bt2ydvUTiNTkS3o2emPIge8ZC9/smEybI0OahQoQJ79uzB09OTChUqPPhgFgunT59OdzAZSZKDlHrTm1vcYi1rzQ7FPl1CVyPci57y9yJ3q+bZm5gY3cnu999hzx6dGBw9ChUqwEsvwfPP6yZyp+zUESJ1YYRRkYrc4haQwf0Q0uvkSdi2Td+q2bIFYmPBxwfatdPfa9bUD5d0ZpSz0PNn/AqcRLdcyZ8x8Rjs5rZCXFwc+/fvp1y5chQtaj9tYZIcpNSFLjjhxApWmB2KfRmFnrnwefQUx53RpXLtyc2bsGOHTgZ++00nBjExur5Agwb6nnjbttCsWY5ICFLTkIbsYAeVqcwJTpgdzl137uiWhI0bYfNm3a8jPl4P+5w4EQYNSl8Rp+RpuRsCTwNvA3HcnQNCiFSYlhy8+OKL+Pv7M2TIEOLi4mjYsCF79uzBzc2NVatW0aZNm3QHk5EkOUipAx0oQhG+4iuzQzHfYXTBon7A5+i+BR1NjchWQoK+4KxaBb/+qkcTKAVVqugOhI0b6+9eXtnuVkF6JZJolP6ezGQmMtHkiB4gLk53bpw7F+bM0UMl587VLQrpcZm7CYE3sA2oljGhipwno659af6IsWbNGvz9/QFYvXo1oaGhXLlyhSlTpvDWW2+lOxCR+aRCInAI3QnsLHrUQSJ6ul57SAzi4/XQwhdfhKeeghYt9HDDli31ve4rV/RIg3nzoH9/3cHQQRIDACecjAmbJjEJS9JXBBEmR3YfV1fdv+Pzz+Gvv/S22rXhhRdg/XpdYTItn8lKoCd/yo1uPaiKvt0wLoPjFuIeaU4Obt68adw+2LBhA927d6do0aJ07dqVo0ePZniAIuNEEeW4yYECrgN1gD/QHb22Y35tgsRE3UIwZIi+RdCqle4dP2KE7kNw6BB8+KHu9GZC3QF79Cmf2qwXprBJkTwGb29dR2LxYl1XoVMnXbXRywtGjoS9e3VS+DjyAcPQHRRDwZg/7Td0B0YhMlCa/zSWL1+enTt3EhkZyYYNG2jbti0AYWFhuGdyFTXxZLrQhWaOWOx9LDAV3clwL/oergVze4ErBd98oy8STZvqPgRjxsCJE3r44fjxeqSBSOEVXkGhiCEG0CMaLFgYbkxyYGecnKBnT91f5MoV3ZGxRw89xDQgQHciffttXZjqcXVBz+SZiJ7lM7mPcdwDXyFEmqS5z8GyZcsYOHAgbm5u1KpVi19//RWLxcKMGTP46aef2Lx5c2bFmibS58DBhaGr0Q1BlzmOBp4xNaK7LlyA4cP15EIDB8KoUVCrltlRZUvlKc/f/G2s16QmhzhkYkRpEB+v+5IsX65vFeXKBe+9B336pK1TaTi6ZPd+4Dl0fxp7r8UhMo2poxUuX77MpUuX8PX1xTmpF+6uXbvw8PCgmp182pHkwEHFolsEzqArzgWjJ0KyB4mJumPam29C0aJ6uXlzs6PK1m5xizGMYSMbjSQhuV9CthIWBuPGwZdf6tsOydUYk+areSxX0IlwP2A1esRNs8wIVtgzuxnKaK8kOXBQLYAOwGuYP+zr1Ck9idGNG+DmppuR//oLXn0VpkzJ9MmMHM1WtlKf+rzP+wxnOMUpbnZIaXf4MMyfr0s3h4XpPgodO0JQEKRlMruBgB/wCrqvjR13yxAZS5KDR5DkIKViFOM93mMQg8wOJWOdRd9zHYHuxV0WPQeCGS5cgM8+04nA0aNQuLCunBcdre8vjxmje66LTDGVqbyFHjX1Nm8ziUkmR5ROMTGwdKnul/LLL3rbc8/pxLJu3cc7hgKOossyH0XPDSJyPEkOHkGSg5TccGM2s+lHP7NDyRjR6JaBg0BfYAdZe681MhL+7//0pzxnZ6hRQxcqKlJED1vr0EF3NnzS+vvisYUTTsF7KljFEYeL3de9foTbt2HFCpg9W88g+cUX+pbD44gDfkS3pm1I2tY+c8IU9kGSg0eQ5MCWQuGEE9/wDV3panY4T04BjYDewEtJ61k5+uCPP2DAALh6FcaOhagouHhRT/Pbu7e+jSBMc5KTVKaysf4JnzCQgeQlr4lRPaHERJg8Wd+SatkSKlbUVTE7dXq8BHQU8BTwJtlnRlGRZg6XHDzzzDOsWbOGLVu20KpVq0fuL8mBrSiicMedNawhiCCzw0m/m8BmoBuwBaiCrm6YFaKi9Ce4b7+FH36A1q11B7K0dBoTWeZB8y9kyw6L9/r+e10U6++/dU2McuXgnXd0J8ZHFcVSwAmgKbpCqPkzdIsMZlqFRDMsWLCAqKgos8PI1qKJBsj+RZCOoKe0DUfXm8+KxCAhARYs0KWLhw7VQ9C+/x42bZLEwI4pFDe5yf/4n832ujzmPXt79fzzsGyZvoV16pRuRejTR8+8+fvv+v/rg1jQkzhNQycGu+CekaBCGOw+Ofj777+ZNGkS8+bNMzuUbC0KnVzlya5Tu70KLAIC0VPbFsiCcyoFGzboUriDB+sZ986c0UnBc885VOni7KoABahPfRTKmPZ5N7u5wAWTI8sgFSro1qsff9SVGBs3hsBAOHDgwa/JB/RPWp4ALMj8MEX2Y9fJQWJiIv369WPy5MmULv3w9q+4uDiioqJsHuKubNty8E/S99pA/qTlrHgLR45Amza6U6GXly5jPGeOnvNAZEsWLGxhCwBlKMMd7pgcUQZ6+mk9Z0NwsE5aAwL0PA6Psho9R8NZoAfY2zQVwjx2nRzMnDmTfPnyMWDAgEfuO3XqVNzd3Y2HZ1rGBDuAcpTjH/7BF1+zQ3l8J4GK6GFYfYFns+CcR4/qceU+PhAaqmdEXL0aqlfPgpOLzNaKu/2V8pP/gf0SsiUnJz1M9o8/9G2GLl30/92HyQPkQlcULQS4IwmCAOw4OTh69CgffPABc+bMeaz9x48fT2RkpPEICwvL5AizF2ecKUEJcpHL7FAe7W/0ZDKV0BXfsqro5vbt0KCBniBn8WLYswearCcq+QAAIABJREFUNMmik4usspSlNutFKGJSJJnE2Vm3cvXtC88+Cy+//Oh5G/yBL9CjGPyBFZkfprBvdpsc7Ny5k8uXL1O2bFlcXFxwcdFjldu0aUOvXr1S7O/q6kqePHlsHiKb+hZ4C92zuilZM0RxzRo9+qBlS/3Jq0cP/UdW5Dg96YlCGeWWwwhjHjmsT5Ozs66LMG+e7jxbvboeZfOowWluwL+BTuiZH89meqTCTtntUMabN29y4YJtp6FatWoxd+5c2rZt+8g+CDKU0dZZzrKa1bzIi/bZ7yABnRD0Rc9XH41u4swKX34Jw4bBoEEwa5YkBQ4kmmijk+5pTlPBtNKamejGDT0F+JIlusPirFlQs+ajX/cKcI67Mz6KbCHHD2UsWLAg3t7eNg/QU0Y/KjEQKR3iEK/yKgk8ZJiTWRT6f+IN4HjSclYkBomJupzxkCF60pvZsyUxcDD3Vk+sSEUTI8lEhQrp22T/+x/cuaNnAB0+XA/JfZgZwJfoycxmgz3+6RCZx26TA5Gx7Ha0win0rYMw4L9k3bTKERF6vPjMmbBwoS4iI0MTHY4LLjZFkUYzGgsWYogxMapMUr8+7N4NX32l/89366YT5AdxR89o+icwE7iVJVEKO5Gtio7b6R2QbCGaaFxxxRk7+mQci54Mpj5ZW/r44kU9y93ff8NPP+mmVuHQoonmEpeM1gNvvAkhxOSoMoGzsx7JUK4ctGoF/frpFrOHzRDaBDiAnsdkHHpK6KpZEq0wkbQcOIgoouyrANJCoDn6f+AMICtGnsbEwPLlerhXRIQuPSuJgQByk9tmWONJTmLBwklOmhhVJmrSRHfCXbdOj9AJeUQilAs9xHEv8JDGBpFzSHLgIKKJto9bClHozoadgOFkftuVUrBrlx7O9dRT0LOnHpWwaxdUqpTJJxfZSXnKo1BEE83apF54lamcs4ol3atdO100ycUF/P31fA0Pkw89r0l1YDo6wRc5liQHDiJ54iXTDUHPDucJpByR+mSU4v/bu/v4muv+geOvs5nduqlt7kKjGxKXm+SeRJGbIiKRm+W2UopUlxD9IpTUlVy5CVFKyZVqEXI3DYlhXW5yH1eWmdjYzLbz/f3xtp0ddr9zznc75/302MPO93zPOe9v2dn7fG7eb5YskX733btLQlCnDjRtChs2wNixMpXw6adQtqyDX1y5C198eZiHM29nFEtyq4JJGWrWlK27/fpJm/H+/aXOR16S0aZNbk6TAw9xF3fxOI+bF8A5ZFfC68DzTnj+5GQYNEi+atWShVbR0VJWdscOKYf86qvaKEnlm4FBfeozjWmZx+Ywx8SInMTPT7Y3fvUVREZCo0a592YAmAQ8ACwHRjg/ROV6xbbOQVFpnYNixArUR0YNnJEYnDwpjZCOHoVly6BzZye8iPJk61hHBzqwk500prHZ4TjPxYuyWPfQIakYeueduZ8fAfwFPOWK4FR+uH2dA+UmziH/yhYD/Z3w/FFR0mTm6lUpd6yJgXKCB3kQA4PlLHfP6YUM5crJIsVq1eDee+HLL3OvqtgFSQyigM7gjjtAPZUmB8p5EoE6wPdIvfabHPz869bJdqxmzSRJ0AWGysne4R1AOjy+xEu8wAsmR+QEZcvCpk1SQvzxx+GFF/Iuu1wWaA0loXWLyh+dVvAQ4YRTilLMZ75rXvAiUA7YAjQFfB343BkLD4cNk0VUixfLimulnMzAwCubz1RZCym5lWXLZJHis8/Ce+9J58fcHEGmDj/D8R8GVL7otIIqkDjiuMpV17zYKaAGsBspoOLIxCA2VhKCQYNg5Eip9qaJgXIRCxasWBnAAF7iJbvjf/KniZE5Sd++Uhvko4+kYFJ6HjWUvZEuqroZqMTT5MBDJJHkmq2MqUjVwyXAPxz4vFeuwNSpMnWwYwesXQszZ2ovBOVyFix8wie8zduc4lTm8Vu4hV3sMjEyJ+nVC77/XtYfjByZe4JQA3gX6ZMyCLjgkgiVE2hy4CEuc9n5ycFR5FPDIaArjilwZBjwww/SRW7KFGmUdPCgbFFUymRVqYqBkblI8XVeNzkiJ+nQQaYYFi6U4knnz+d+fgJS7Exz9xJLkwMP4ZKRgxrAWCCsCM+RkCBTBeHhsiCqTh3o0kUquB08CBMmQGCgY+JVykGsWEkllQAC3Hc3Q8+esHUrHDggtRAiInI+tybwBdKPYQJc6/umShBNDjyEU5OD00hBlNNIQZTCrjE4cADq1pWFhidPSv+DTp2kB8KXX2oBI1WslaIUX/EVAHvZa3I0TnLvvbJluEkT6NpVdjKkpuZ8/klgHdJ1VZUoupLLQySRRCBO+sQdjKwvKFOE59i0ST6Z1Kolbz4VKjgmNqVc6C/+oiIVaUAD993BULGiJOtLl8KIEfD779LEycfnxnNrAduQ6qhfAz1wbQdWVWg6cuAhoojiSZ507JMmAaOR1svvUritS+np8MYb0L49tG0L69drYqBKrArY/u36409HOpoYjZP17y89SzZuhHHjcj7PAvyCNFr7n4tiU0WmyYGHqEENbuZmxz7pJWA/kI8+LdnasUNax775JsyaJV3hcusrr1QJsJnNgHRCXctaalLT5IicqGlT6cvwzjuwcmXO5zVDFixXRaYaVLGnyYEqnNXISMEaZIdCQfz0E/TuLZUNvb1lTcHzz4NFxxtVydeGNhgYHOc4AMc5zipWmRyVEw0aBEOGyCLiI0dyPi8I+A5ogBRJU8WaJgceIJlk5jGPM4X+iH+dRGTh4XeFeOyMGVLy+Ngx+M9/YPNmWfmslJsJIyyzMNJGNmLB4rpCZK72wQeyYHjUqNxLLXcCNiPVU60uik0ViiYHHuA85xnOcI5xrOhPFo98AohBFhcVxPz5Uqdg1ixZdNi9u44WKLdWmcpc5CL+SBlbX4eWCy1G/Pxg9mxYswYGDMi5UFIpZPHyV0ArIM11IaqC0eTAAySRBFD03QpW4H5k8WFBy6MuWgTDh8PEibL9SSkPUZayvMVbmbf/4dDSocVI27ZSSXHFCplisOYyNNAAGIXulyvGNDnwAJe5DFD0OgdewOdAvwI8Jj0dXnwRnnoKXn0VJk0qWgxKlVBRRAEQQwwWLExlqskROUGnTvDNN/D557JIMSd3AI8DG4AvXRSbKhBNDjxAxshBkZKDfwHvAHcDlfL7wknQo4esZv70U+mNoNMIykM1pzlppBFJJACv8RpX3LF0YMeOsgNp/HjYvTv3c1cjuxhUsaPJgQfIGDko0rTCbUgp1PxKTJSyx5GRsjuhX0GGG5RyT95404pWTGISAA/ipj1CXnoJmjeXXUlnz+Z83gzgn0iCcMlFsal80eTAAxRp5CARWAt0QeYI8+PiRfn08N//SuXDli0L/rpKubGMBk1b2Yrl2p/TnDY5Kgfy9obly2VasWtXuJTDb34LsijxQeATF8an8qTJgQeoQAV60IPSlC74gyOR1qsJ+Tw/NhbatYMTJ2Sb4j/cdPGVUkWURBIrsRUOqkY1trDFfcouV6oEP/4IR4/KYuSctjiWAn5CKiiqYsNiGLltSi25kpOTCQgIICkpCX9/f7PDKZlSkamEy5CvGYnt26WTore3bGm6/XbnxqeUG/iDP2hEI+KzdCdayUr88KMTnUyMzEHWrIHOnWWnUm4LkvcjCcJKcHQxV0/iqN99OnKgspcA3IV0VMtPYvDee9CihSQE27ZpYqBUPlWnOuc4x5u8mXmsBz3oTGde4RUTI3OQhx6Cjz6CyZPl75xUAloCfq4KTOVGRw5U9gxki1EH8m6otHSpFD6ZPh3GjtUdCUoV0QAGsJSlAHzO5/Shj8kROcBrr0kBtJ9+ksWKOTkDHAbauCow9+Ko332aHHiAmcxkAxuIICJ/D4hEkoP8/HDu3y893keMgJkzixClUiqrlaykJz0zb5f4tQipqfDooxARIaXTu3fP/rzngb+A5a4Mzn3otILKtzOcIY64/D9gO/BhPs67fBl69YJ69WDatMKGp5TKRg96kE46TWgCgAUL85lvclRF4OMjOxgefVS2OP76a/bnTQe+gJKeC5V0mhx4gCSS8r+NMQUYi1RCzMtzz8Gff8oPvE9BiiAopfLDCy+2sz3z9jCGmRiNAwQGSnnl5s1ltDG7Hgz+wFngXiDaxfGpTJoceIDLXM5fAaSTQBjyA5nXv4zFi6VfwuLFcOutRYxQKZUTCxasWHmapznJSXaz2zFN1Mzi5QVz5sCuXVJWPTuhSHnlMBfGpexocuAB8j1yUA2YBdTN47yYGHjmGamC1q2bAyJUSuXGgoU5zOEMZ7iHe7iN27Bg4SQnzQ6tcO6+W8orf/EFHD9+4/1eyAimFWnxrFxOFyR6gC50IZRQFrM455N+Qeb4mubxZAkJsgCxQgXYsEGnE5RyoTTS8MeftOt6HaeQUrgiZ2ZKTYU6daBJE/jss+zPeRlpD7/alYGVbLogUeXbZS7nPXLwE2TZZp09w4ChQ+HvvyXj18RAKZcqRSlSScXAII00vK69hfviyzd8Y3J0BeTjIw2aPv9caqNkZzLkd5OVciwdOfAAhzlMaUpzK3msDUgHvHO5f/ZseP55WLsWHnjAkSEqpQrJG2+sWAF4hVeYRgnaOWQY0ub58GHYt08WLF7vGDASWAoEuzi+EkhHDlS+3cEdOScGaUBbpLlSbonBsWOyxmDCBE0MlCpG0kkniigApjM9cwFjiWCxwMKFEB8Pr7+e/TkhwJ1QUi7JXWhyoKAPUCOPc158EWrUkCpnSqlipTnNWcOazNveeJecoklVqsCMGVI98dChG+8vC7wH+EKW9hPKyTQ58ABf83X2W5+uAIeAEcAduTzBmjXw7bfwr39B6RK26EkpD9GRjqSTTilKMYEJHOaw2SHl3+DBULs2jBmTffdGAxnhnOXiuDyYrjnwAEEE8T7vM5jB9ndEAL2QWublcnjw1atSAbFOHSl5qpQq9uKJJ4QQAK5wBV98TY4oHzZulHbvK1ZAz5433v87UBNp8axypGsOVL5YsXKZy5Sl7I13dgGOknNiAPDBB3DyJLz7rpMiVEo5WnCWlXt+JaXN4f33w8CBsug5Ppv5gzuRAm1fuTowz6TJgZu7xCUAylDG/o5twMdIm9ScnDkjbVZfflnWGyilSoxznMv8/t/828RICuCdd6SCYt++YM1mBeIPwCZXB+WZNDlwcwkkANkkB6eQHQq5dVceORKCg+EVN+gpr5SHCSY48+d/37U/xV5ICHz1FaxbB8uW3Xj/RPLXFE4VmSYHbi6RRCCb5KA3ubdEXbECVq6E+fOz33uslCr2ylCGBBLwx5/61MeS66eBYqJZMwgPl51RiYn291mQaold0a6NTqbJgZvLSA7s1hyMAcbn8qDz5+HZZ6UpitY0UKpEK0MZxjEu83aJSBDefBOSk6Xz6/VuATohRduU02hy4Ob88acDHShPedvBbuTeQ+Hll2Xe7513nB2eUsoFQgjhClcyb0cRxQ/8YGJEeahcGT78ED75BHbssL/vH8CzkGVJhXIC3croaf5C2qHmlBbu2wcNGsDSpdCvnwsDU0q5StbRAyvW4jmaYBjQtq0sTNyyRaopZvgZuA9pM3+LOeEVV26/lXHq1Kk0atSIoKAgKleuTHh4OHFxcWaHVfL1AZ7J5f5Jk6B+fXjiCRcFpJRytbWszfy+He1MjCQXFgtMmQJbt8oCxawaA3vRxMCJiu3IQefOnenbty+NGzcmISGB5557jsDAQDZs2JCvx+vIQQ4uAAlA9Wzu27MHGjaEb76Bbt1cHJhSytWGMIQZzGAXu+hAB9JIwzvXJism6NZNyirv2QN+WWo2XAKigA5mBVY8Oep3X7FNDq63bds2WrRowYULFyhXLreqPUKTAzGPeXx37Q+7kIJHt+dw8qOPwokTsHu3/RCei1itkJAAZcvKkgew1UIpXx68i9l7llLuwIrVLiE4yEFqUcvEiK5z/LiMZr70EkycaDv+LTAUOALXb8byZG4/rXC9c+fO4efnR2AO2+pSU1NJTk62+1JwkpP8wR9yYz4wIYcTo6NlxGDSJKcnBosWwejR9tWYIyLkl/9NN0ntpQzh4bL1uWtX27ELF6BHDynDfvCgU0NVyu154cVv/JZ5uza1SaYYvX/WqCG7Fj78UHYwZOgM/A9NDJykRCQHKSkpvPHGGwwcOJBSpbIvrD1lyhQCAgIyv4KDtfE3yFbGzBoHHwGfZHOSYcC4cbIQ8ZFHHPbaZ8/KusbnnrPvpTJ/PkRGQtYlJI0awfLl8OOPcPPNtuPTp8t0Y9Zurqmpkkhs3iyJQoYJE+DOOyWhyOr06eyLrSmlxN3cjYHBLGZRjWoc4hC/87vZYdmMHAmXLknztwylgIOAbqpyimI/rZCenk6fPn04ceIEGzduJCgoKNvzUlNTSUtLy7ydnJxMcHCwx08rDGIQf/EXqyNWw0WgbzYnrVgBvXrBpk1w332Ffq30dJkWuOkm29P27QtNmsAPP8h0AUii4IzBichI2LBBEofx1+o4XLkCAQHSTHLRIts6yyNH4M8/JZmolFsJaaU8UG9689W1JgaBBGaWYTfVxInw/vvw++9QsaIci0DKwH+JNmS6xiPWHFitVgYOHMiePXvYvHkzN2f9SJkHXXMgetKTUpRi+dzlkmVf3/I0MVFapT74ICxeXOjXmTBBRgQ6dIAlS+TYlSvyiT0goNBPW2RpabJN+vBhaNkS7rjWmnryZJlBCQ6Gc1n2S8+aJcdatYKaNU0JWSnTbWEL92H7oHCCE/zGb3Shi3lBJSZCrVrw0EOwcKF5cRRzbp8cGIbB4MGDiYyMJDIykkoF/HinyYHoQAeqU50FLMj+hBdflEIjhw5BaGi+n/fsWamqnLEE5P335We3Z0+46y4HBO5kVquMHMTGQuPGcswwZGblwAGYMQNeeEGOb9ggeVNYmCRBPj5yPCUFfEtAJ1ylCiOJJHaxi0QSeYEXOMxhfuVX7uEe84JatAiGDoWjR+HWW+XYf5ARhBze4jyN2y9IHDFiBN999x2fffYZALGxscTGxpKerjUzCyKBBO5fer9UFLtedLTM4c2Yke/E4OpVGDIEqle3H2gYNUqG8ktCYgCyG6JqVVtiADLVsXcvJCXBiBG246mpsj4iIgKyLnlp3FjWR7z5pu3Y+fNSr+V//9N1DqpkCyCA1rSmM535lm8BaExjBjHIvKD69ZP3qnnzbMduBuqivRYcrNiOHFhymJQ+fvw4YWFheT5eRw7EPvYRvCGYWw7fAsOz3JGaCi1ayGR8ZKRt72A+hIdD06YwYIC5UwZmW7tWpivq1JFW9ADffmsrEXHypCRRAJ9+KiMrjRvDvfeaE69SRfEzP9OKVgAsYQn96W9OIGPGwKpV8sNnwpbr4s7tpxWKSpODPEyaJFsBdu/O9eP+iROyU2D0aNlqrHJntcpWzKNHZY1DRm2GAQPk/axbN9uajBMn4Omn4fbbZQt3xiip1VqgXE0plznIQe7iLhaxyLwRhMhIaNNGFiZmLCLqDDwCjMjtgZ7B7acVlAO1AVZmub17t4yFT5+ea2Jw5YqMEPz8swyXq7x5ecEtt8h7V9aiTUuWyLbL+fNtx9LTZYR01y5ZOJlh8GDpOzN4sP25a9dKXYeUFOdfh1LZqU1trFjpS18sWBifa3tXJ2naFIKCJNvOMBho7fpQ3JmOHLixVFKJMCJov7Q9ZRqXgTrX7hgwQBos7d6d7UfUrFsNd+6UEYPSpV0Xt6eLjJT/7uXK2RKE06ehWjX5/j//ge7dbd9v2SIbToYPz/75lHK0WGKpTGUAjnOcMMJcG8BLL8HHH8ue5IyaNgYUx/5RrqYjBypPccTxqOVR9gzYY0sMLl6UAgRDh96QGCQkyC+jDz6wHbv3Xk0MXK11a5nGyTpycMstUko6Oloa1WWIj4eoKFi2zHbMMODuu2UdRESE7XhiovzvV6qoKlGJ7WwHoAlNSMHFw1kTJ8onmA8/lNuzgYauDcHdaXLgxv7mbzpHdKbprU3JbOX+xRcyqd33xmpIgwfLSF317JoyKVNZLLIzokED6TORYcgQqeOwebPtWHo69O8v07FZF4wuXiyPrVzZfidFRARs3CjbOpXKr6Y0ZRaziCMOP/xIJx3DVVsGypaFYcPgo4/kH3Mn4EPXvLSn0GkFN7aVrTwW+xj/3f5fgrtfG3pr0kRWwGX9qHnN8ePg768VA93VmTMy8nD2LAwaZDtevTqcOiVbUf/v/+TYL7/IDtfq1WVBakavs0uXpLaFLhJXAAYGvehFTWqyhjXEEMNWttKSls5/8e3boXlz2bVw++1wFZlW8HH+SxdnOq2g8vQ3fxMfHI9/92v/QPbtk8nsa+PV6enS+yAjPaxRQxMDd1a5MnTubJ8YgOyaiI2VWhUZMv5NREbaTyt17SoJ5MiRtmN//y3Tv+vWydSF8hwWLKxgBTOYQSXkzaMVraQLrLP94x8yNbpzJ5wF/IGfnf+ynkJHDtzYEpYQ2iWUh255CMs8CzzzjLyDHzrElate9OsnPQ9+/VXmqJXKy/btUkyzcmUplQ2yrrV1ayketXs3NLw29zttmkx5NGsGr7wixwwDfvtNHh8crCMQ7mY722lOcwCuchUfZ3+Mb99ehrJWfSst6esAHlx7BXTkQOXDec4z+V+TsYy1yEe6pUul9J+XF1arbFVcu1YTA5V/zZrBwIG2xACko+alS9Kjom5d2/EKFWTU4exZ27ELF+QDX2ioLH/J8N13UvNh2jT714uLs9/mqYq3ZjRjLnMBSQ4A565DGDZMFs2cPgW1AB25chhNDtxYsDWYxgGN4Xakf0JqauaYckCA/Ey11r3BygEsFhkJ8MnyQfGpp6QN98yZtmNlysjIwbp1tqqSIDtl9u+Xlt1Z3X23JBjTp9uOxcRIj4s5cyA52XZcK6sXD8MYxiUu8Tu/Y8GCF15Yrv2xYuUCF/J+kvzq3l1W2S5ZAk8CzzvuqT2dTiu4s1NAdWB9Cgy8jfSHuzOnzmyeftq+R4BSxVVkpDTIuusuGXEAWLNGtrmfOSP1HzJ+vO+7TxKH8HBbQnLxIqxfL9Un69YFPz9zrsMTneEMVahid+wQh6hFrczbVqxYilqcYMwYyUJ3nYKyFll74MG0fHIeNDkArMBpYMUcGDeacUPjmDmvDDt32t5olXIX69bJwvWwMFl4CbYF7QB79thKgGesh2jaFF59VY4ZhpxTpYpMe2gJa8fJSAKuchU/7DO0Ik87ZJRT3nsMkmpAs6I9XUmnaw5U3mKBuPMwbRLpw55m/6kyzJmjiYFyTw8+KGtuMxIDkDUSyclShj9rpfCQEGm3nbUs+Pnzsn6iUiX46ivb8VWrZGr7rbfsX+/sWV0PkV8ZUwu++GJgkHhtccAoRtGMZpziVOGfvHFjmc9ackxKxeu6A4fQkQM3tmTyEh75qDXlvVrBf/+LUa68rg5XKgdpadIw68wZ6bRZoYIcX7YMFiyQXhnr1tnODw2VCpXTp8PYsXJs3z5JLCpXlumNjLcebaaVvZWspCc97W4/yqMFf6LmzaFeQ/j3HPDO+3R3ptMKefD05MDAIDA9gLkDvXmyx1IsPQrxA6eUytHWrbb1EPXqybHVqyVR+PNP+cpY49CmjSy4fOopKS4Fsrtj2TIZqXjgAamA6YlGMYp/8a/M24Vah/Dqq7LuYMExCLPAbQ4OsgTR5CAPnp4cJJLIw5sepsf4O1iaMo8NGyyUKWN2VEp5prVrZWqjZk3btEd0NDz0kGzX3LPHNt03ZowkGS1ayIgFyMjDBx9AxYrQqhVUrWrOdThLHHHEEMMVrtCZznk/4Hpbt8rWq9uS4J/+0qXRQ2lykAdPTw6OcIQrda/wzalKbGsVYteARylVfKSlyZRDxrTDd99J0b+KFeHZZ+XYhQvSeTMuTupD9Oolx997T6Y1wsJg2zbbc86dK9MgTZrYkg6rVbacFvepxayFlPI9inD1qiwimTxZmjJ5MEf97tMNbW7qrPEXLffcx/Sh3zB9TFezw1FK5eD6bcUPPyxfWZUvLyWu09Ptm2bdd5/8nXJdU8SPP4YDB+Dll23JQUQE9O4tayUykg+QqY34eNnJ0aaNHDMM+TJjnUQwwZnfe11bM7+b3QQSyJ3cmf2DSpeWaok/nJEqiY+5IFA3p8mBm7p44A/6Rj/OqD6l8a2b9/lKqeLP21u+MjRsaCtXndUvv8jfWROJBg2kieHZs/adPX/8UaYxOnSwJQdHjshaitBQWLnSth30xx9l7USNGlJ/KENKinxwd4Q7uAMr1szEoB71aESjzPutWDnJSRJJpB71bA9s2RLeTYEv0OTAATQ5cFNX98TzxsRJ+K718ELjSnmwrJ/8q1WT0tfX++QT+TvrBHNoqBQd/OsvmbLIsG0bzJsnz5WRHFit0t7A3x9mz7a9xs8/yzbQ0FBZR5ERS3q6fJ/b9IYFi139g050Yg1r5Jqy7MCfy1yGMYx44lnX9RJ93ngXZj0LVMvjv4zKi645cFOf//M/rFz6J1/+8QwWr2I+yaiUKlEMw/bLPT0dvv9e1kO0aCHbQAFWrJDaEJcuSbOuDLffLp1Ax4611Y44cADefVfqT4wZI3+DjHL4+0vZbQPDLjHITvm/y7Nw1Kt0nT4an8qe2btZiyCpXF35ujaV/npYEwOllMNl/dTv7Q3dusGQIbbEAOCxx2DXLvvEAODDD2W0oqetvAGXL0vCsGaNfWGpoUOhbFl45BHbaELiJYPOT56nZucDsKe+7eQ99WHO05SJ6UD9v7McV4WiIwduaktwNKlefrSPuyvvk5VSqhg6dEjWPwQF2RZf/vEHDB4sdSEenLGOobdeaxE6ZD58PES+NyzO7QZZjOlWxjx4enLA999DuXLadlEp5dYMDCxYpAV92hq6efXAKHVBdjB4IE0O8uDxyYFSSimPo2sOlFJKKeUUmhyaC4W2AAALXElEQVQopZRSyo4mB0oppZSyo8mBUkoppexocqCUUkopO5ocKKWUUsqOJgdKKaWUsqPJgVJKKaXsaHKglFJKKTuaHCillFLKjiYHSimllLKjyYFSSiml7GhyoJRSSik7mhwopZRSyo4mB0oppZSyo8mBUkoppexocqCUUkopO5ocKKWUUsqOJgdKKaWUsqPJgVJKKaXsaHKglFJKKTuaHCillFLKjiYHSimllLJT7JODadOmUaVKFQICAnjkkUeIjY01OySllFLKrRXr5GDRokW8+eabzJ49m6ioKBISEnj88cfNDksppZRyaxbDMAyzg8hJo0aN6NSpE1OmTAHg2LFj3HbbbURHR9OgQYNcH5ucnExAQABJSUn4+/u7IlyllFLKVI763VfKgTE5VEpKCnv37uXtt9/OPFazZk3CwsLYsWPHDclBamoqaWlpmbeTkpIA+Q+llFJKeYKM33lF/dxfbJOD+Ph4rFYrFSpUsDseGhrK2bNnbzh/ypQpTJ48+YbjwcHBTotRKaWUKo4SExMJCAgo9OOLbXJQ0Kzntdde45VXXsm8ffnyZUJDQzl37lyR/gOVVMnJyQQHBxMfH++R0yqefP2efO2g16/X79nXn5SUREhICEFBQUV6nmKbHISEhODl5XXDKEFcXNwNowkAPj4++Pj43HA8ICDAI/+BZPD399fr99Dr9+RrB71+vX7Pvn4vr6LtNyi2uxV8fX2pX78+GzduzDx2/PhxTpw4QdOmTU2MTCmllHJvxXbkAGDkyJGMGjWKe+65h5o1a/Liiy/SunXrPHcqKKWUUqrwvCdNmjTJ7CBy0rBhQywWCxMnTmT27NnUrl2bpUuX5nsuxcvLi7Zt2+Lt7e3kSIsnvX7PvX5PvnbQ69fr1+sv6vUX6zoHSimllHK9YrvmQCmllFLm0ORAKaWUUnY0OVBKKaWUHY9IDubPn0+LFi0oV64coaGh9OzZk2PHjpkdllN5YjfLqVOn0qhRI4KCgqhcuTLh4eHExcWZHZZpunfvjsViYf369WaH4lK7d++mffv2BAQEcNNNN9G7d2+zQ3KZCxcuMHjwYCpVqkRQUBAtWrRgy5YtZoflFCtXrqR9+/aUK1cOi8ViVz4f4Pfff+f+++/H39+fsLAwFi5caFKkjpfbte/Zs4fevXtTpUoVAgMDadiwIStWrCjwa3hEcrB582YGDhxIZGQkP/30E1euXKFTp06kpqaaHZpTeGo3y61btzJ69Gh+/fVXVq1axf79+z3iurOzaNEij+wrcuDAAdq1a0erVq3YuXMnUVFR9OnTx+ywXGb06NHs3LmTb775hr1799KkSRO6du3K33//bXZoDpeUlES7du149dVXb7gvNTWVLl26EBISws6dO5kwYQLDhw/np59+MiFSx8vt2qOjo6latSrLly8nJiaG8PBw+vTpw6ZNmwr2IoYH+vPPPw3A2Lt3r9mhOEXDhg2NcePGZd4+evSoARjR0dEmRuV6UVFRBmBcuHDB7FBc6sSJE0a1atWMU6dOGYCxbt06s0NymR49ehiDBg0yOwzT1KlTx5g1a1bm7YSEBAMwtm3bZmJUzrVx40YDMFJTUzOPrVq1yvD19TUSEhIyj/Xv39/o1q2bGSE6TXbXnp0OHToYL774YoGe2yNGDq537tw5AG6++WaTI3G8jG6W7dq1yzyWtZulJzl37hx+fn4EBgaaHYrLWK1WBg4cyOTJk6latarZ4bhUeno6a9asoUaNGrRt25aKFSvy4IMPsm/fPrNDc5nmzZuzatUqzp07R3p6OgsXLqRKlSrUrVvX7NBc6pdffuHee++lTJkymcfat2/vce+BGc6dO1fg33celxwYhsH48ePp2LGjW755FrSbpbtKSUnhjTfeYODAgZQqVawLgTrUrFmzCAoKIjw83OxQXC4uLo6kpCTefvttnnjiCVavXk21atVo3749Fy9eNDs8l/jggw8ICQkhNDQUX19f3nrrLSIiIorchKekOXv2bLbvgZ64Bunrr7/mwIED9OvXr0CPK9HJwYgRI7BYLDl+tW3b9obHjBkzhpiYGBYtWuT6gF3A0JpWpKen8+STTwLwzjvvmByN6xw4cICZM2cyb948s0MxhdVqBeCxxx5j+PDhNGrUiLlz52KxWPj2229Njs413n//fQ4fPsy6devYuXMnTzzxBI888gjx8fFmh+ZS+j4ooqKiCA8PZ8GCBdSoUaNAjy3RH6mmTZvG+PHjc7zf19fX7va4ceP48ssviYyMpHLlys4OzxQF7WbpbqxWK4MGDeLgwYNs3rzZoz4x7dixg9jYWKpXr253vGPHjvTp04fPPvvMpMhcIyQkBG9vb2rVqpV5zMfHh5o1a3Lq1CkTI3ON5ORkJk6cyPr162nTpg0gJegjIiJYtmwZzz33nMkRuk7FihU5ePCg3bG4uDhCQ0NNisj1du7cSefOnXn77bfp27dvgR9fopOD8uXLU758+XydO3nyZBYsWMDmzZsLnEGVJFm7WbZv3x7wnG6WhmEwZMgQtm/fTmRkpFuuKclN9+7dady4sd2xevXqMXfuXB566CGTonKd0qVL07BhQ44cOZJ5LC0tjRMnTtyQMLmj1NRUUlNTb6in7+XllTmq4imaNGnCzJkzuXTpUuYHhA0bNrj9e2CG6OhoOnbsyPjx4xk+fHjhnqSwqyRLkrfeesvw9/c3Vq9ebZw5cybzKyUlxezQnOLjjz82goKCjJUrVxp79uwx7r//fqN169Zmh+V0w4YNM0JCQowdO3bY/X9OS0szOzTT4GG7FT777DPDz8/P+PTTT41Dhw4Zzz77rFGxYkXj4sWLZofmEi1btjSaNGlibN++3Th8+LDx2muvGaVLlzb2799vdmgOFx8fb0RHRxvz5883AOPXX381oqOjjcTERCMlJcW47bbbjF69ehm//fab8fHHHxs+Pj7G+vXrzQ7bIXK79piYGCM4ONh45pln7N4HC7pryyOSg1tvvdUAbvjauHGj2aE5zdSpU41KlSoZfn5+RteuXY0zZ86YHZLTZff/GDCOHz9udmim8bTkwDAM47333jOqVatmBAUFGW3btjViYmLMDsllTp8+bfTp08eoUKGCERgYaDRu3NiIiIgwOyynWLRoUa7v6wcPHjTuu+8+w9fX16hevbqxYMECcwN2oNyu/fXXX8/2voEDBxboNbQro1JKKaXslOjdCkoppZRyPE0OlFJKKWVHkwOllFJK2dHkQCmllFJ2NDlQSimllB1NDpRSSillR5MDpZRSStnR5EAppZRSdjQ5UEoppZQdTQ6UUkopZUeTA6WUUkrZ0eRAKVUgbdu2ZezYsQwdOpQyZcoQFhbGDz/8wOnTp3nggQcIDAykZcuWnDx5MvP88ePH2z1HWFgYCxYsMCN8pVQ+aHKglCqwefPmUbduXXbv3k2XLl3o378/Q4YM4YUXXmDXrl0YhsHo0aPNDlMpVUiaHCilCqxly5aMGjWKO+64gwkTJnD+/HkeeOABunbtSu3atXn++efZtGmT2WEqpQpJkwOlVIHVq1cv8/uKFSsCcPfdd9sdO3/+POnp6S6PTSlVdJocKKUKzMfHJ/N7i8WS4zHDMPDy8sIwDLvHp6amuiBKpVRhaXKglHKq0NBQYmNjM2/HxcXZ3VZKFT+aHCilnKpNmzasXLmSDRs2EBMTw5AhQ/D19TU7LKVULkqZHYBSyr0NGTKE3bt306NHD8qWLcuMGTPYu3ev2WEppXJhMa6fDFRKKaWUR9NpBaWUUkrZ0eRAKaWUUnY0OVBKKaWUHU0OlFJKKWVHkwOllFJK2dHkQCmllFJ2NDlQSimllB1NDpRSSillR5MDpZRSStnR5EAppZRSdjQ5UEoppZSd/wf65zUeOdqLzwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x384 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the results.\n",
    "x = np.linspace(0, 10, 101)\n",
    "y_truth = 4.1 * np.sqrt(1 - (x / 5 - 1)**2) + 1  # FIXME: scale is a rough estimate\n",
    "\n",
    "fig = plt.figure(figsize=(6, 4), dpi=96)\n",
    "fig.patch.set_facecolor(\"white\")\n",
    "\n",
    "ax = fig.gca()\n",
    "ax.set_xlabel(\"mu\")\n",
    "ax.set_ylabel(\"sigma\")\n",
    "ax.plot(*np.array(xy_tango).T, label=\"TANGO\",\n",
    "        color=\"red\", linewidth=1, linestyle=\"-\")\n",
    "ax.plot(*np.array(xy_sgd).T, label=\"Stochastic gradient descent\",\n",
    "        color=\"lime\", linewidth=1, linestyle=(0, (5.3, 2.6)))\n",
    "ax.plot(*np.array(xy_asgd).T, label=\"Averaged SGD\",\n",
    "        color=\"blue\", linewidth=1, linestyle=(0, (1.8, 2.4)))\n",
    "ax.plot(x, y_truth, label=\"True natural gradient\",\n",
    "        color=\"magenta\", linewidth=1, linestyle=(0, (0.9, 1.3)))\n",
    "ax.legend(loc=\"upper right\")\n",
    "ax.set_xlim(-2, 12)\n",
    "ax.set_ylim(0, 10)\n",
    "\n",
    "plt.show(fig)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}