01 - metric objects for the perceptron in TensorFlow.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget --no-cache -O init.py -q https://raw.githubusercontent.com/rramosp/metricobjects/master/init.py\n",
    "import init; init.init(force_download=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'2.1.0'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "try:\n",
    "    import einsteinpy\n",
    "except ModuleNotFoundError:\n",
    "    !pip install einsteinpy\n",
    "import sys\n",
    "if 'google.colab' in sys.modules:\n",
    "    !pip install --upgrade sympy\n",
    "    print (\"setting tensorflow version in colab\")\n",
    "    %tensorflow_version 2.x\n",
    "    import sympy as sy\n",
    "    sy.init_printing(use_latex='matplotlib')\n",
    "\n",
    "else:\n",
    "    import sympy as sy\n",
    "    sy.init_printing(use_latex=True)\n",
    "\n",
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "\n",
    "import tensorflow as tf\n",
    "tf.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from einsteinpy.symbolic import MetricTensor, ChristoffelSymbols, RiemannCurvatureTensor, RicciTensor\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# The transformation induced by a single layer perceptron\n",
    "\n",
    "- $dim\\_u$: the number of input neurons\n",
    "- $dim\\_x$: the number of ouput neurons\n",
    "\n",
    "\n",
    "- if $dim\\_u<dim\\_x$ we have an **embedding** of a flat manifold in a larger, also flat, ambient space\n",
    "- if $dim\\_u>dim\\_x$ we have a **projection** of a flat manifold in a smaller, also flat, projection space\n",
    "\n",
    "\n",
    "- **input**: $\\mathbf{u} = u^\\alpha \\in \\mathbb{R}^{dim\\_u}$ (indices as greek letters)\n",
    "- **output**: $\\mathbf{x} = x^i \\in \\mathbb{R}^{dim\\_x}$ (indices as latin letters)\n",
    "- **weights**: $\\mathbf{W} = W^i_\\alpha \\in \\mathbb{R}^{dim\\_x\\;\\times\\; dim\\_u}$\n",
    "- **biases**: $\\mathbf{b} = b^i \\in \\mathbb{R}^{dim\\_x}$\n",
    "- **activation**: $\\text{a}: \\mathbb{R} \\rightarrow \\mathbb{R}$ (or vectorized component-wise)\n",
    "\n",
    "\n",
    "## Metric objects\n",
    "\n",
    "$$\\begin{align*}\n",
    "&&&\\text{general expression}&&\\text{perceptron}\\\\\n",
    "\\text{The transformation function}\\;\\;\\;: \\;\\;\\;&x^i&= \\;\\;\\;&x^i(u^\\alpha, u^\\beta,...)&s^i= \\;\\;\\;& W^i_\\gamma u^\\gamma + b^i\\\\\n",
    "&&&&x^i= \\;\\;\\;& \\text{a}(s^i)\\\\\n",
    "\\text{The Jacobian}\\;\\;\\;: \\;\\;\\;&\\mathbf{J} \\;= \\; J^i_\\alpha &= \\;\\;\\;&\\frac{\\partial x^i}{\\partial u^\\alpha}&= \\;\\;\\; &\\dot{\\text{a}}(s^i) W^i_\\alpha \\\\\n",
    "\\text{The Hessian}\\;\\;\\;: \\;\\;\\;&\\mathbf{H} \\;= \\; H^i_{\\alpha\\beta} &= \\;\\;\\;&\\frac{\\partial x^i}{\\partial u^\\alpha\\partial u^\\beta} &= \\;\\;\\;&\\ddot{\\text{a}}(s^i) W^i_\\alpha W^i_\\beta \\\\\n",
    "\\text{The Metric Tensor}\\;\\;\\;: \\;\\;\\;&g_{\\alpha\\beta} &= \\;\\;\\;& \\mathbf{J}^T\\mathbf{J} = \\sum_{k=0}^2 J^k_\\alpha J^k_\\beta&= \\;\\;\\; &\\sum_{k=0}^2 \\dot{\\text{a}}(s^k)^2 W^k_\\alpha W^k_\\beta \\\\\n",
    "\\text{The Derivatives of the Metric Tensor}\\;\\;\\;: \\;\\;\\;&g_{\\alpha\\beta},\\gamma &= \\;\\;\\;&\\frac{\\partial g_{\\alpha\\beta}}{\\partial u^\\gamma}  &= \\;\\;\\;  &\\sum_{k=0}^2 2\\dot{\\text{a}}(s^k)\\ddot{\\text{a}}(s^k) W^k_\\alpha W^k_\\beta W^k_\\gamma\\\\\n",
    "\\text{The Christoffel Symbols}\\;\\;\\;: \\;\\;\\;&\\Gamma^\\beta_{\\mu\\nu} &= \\;\\;\\;& \\frac{1}{2}g^{\\beta\\alpha}\\Big(g_{\\alpha\\mu,\\nu} + g_{\\alpha\\nu,\\mu} - g_{\\mu\\nu,\\alpha}  \\Big)\\\\\n",
    "\\text{The Derivatives of the Christoffel Symbols}\\;\\;\\;: \\;\\;\\;&\\Gamma^\\beta_{\\mu\\nu,\\gamma}  &= \\;\\;\\;&\\frac{\\partial \\Gamma^\\beta_{\\mu\\nu} }{\\partial u^\\gamma}\\\\\n",
    "\\text{The Riemman Curvature Tensor}^{(*)}\\;\\;\\;: \\;\\;\\;&R^i_{jkm} &= \\;\\;\\;&\\Gamma^i_{jm,k} - \\Gamma^i_{jk,m} + \\Gamma^i_{rk}\\Gamma^r_{jm} - \\Gamma^i_{rm}\\Gamma^r_{jk}\\\\\n",
    "\\text{The Ricci Tensor}^{(*)}\\;\\;\\;: \\;\\;\\;&R_{jm} &= \\;\\;\\;& R^i_{jim}\\\\\n",
    "\\text{The Ricci Scalar}^{(*)}\\;\\;\\;: \\;\\;\\;&R &= \\;\\;\\;& g^{jm}R_{jm} \\\\\n",
    "\\end{align*}$$\n",
    "\n",
    "$$_\\text{(*) keeping standard notation with latin indices but should use greek indices as it refers to the output manifold}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://digitalcommons.latech.edu/cgi/viewcontent.cgi?article=1008&context=mathematics-senior-capstone-papers\n",
    "\n",
    "https://math.stackexchange.com/questions/1809324/visual-explanation-of-the-indices-of-the-riemann-curvature-tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( W, \\  u, \\  \\left[\\begin{matrix}W_{0, 0} u_{0, 0} + W_{0, 1} u_{1, 0}\\\\W_{1, 0} u_{0, 0} + W_{1, 1} u_{1, 0}\\\\W_{2, 0} u_{0, 0} + W_{2, 1} u_{1, 0}\\end{matrix}\\right], \\  \\left[\\begin{matrix}\\frac{1}{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1}\\\\\\frac{1}{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1}\\\\\\frac{1}{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1}\\end{matrix}\\right]\\right)$"
      ],
      "text/plain": [
       "⎛                           ⎡           1           ⎤⎞\n",
       "⎜                           ⎢───────────────────────⎥⎟\n",
       "⎜                           ⎢ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎥⎟\n",
       "⎜                           ⎢ℯ                   + 1⎥⎟\n",
       "⎜                           ⎢                       ⎥⎟\n",
       "⎜      ⎡W₀₀⋅u₀₀ + W₀₁⋅u₁₀⎤  ⎢           1           ⎥⎟\n",
       "⎜      ⎢                 ⎥  ⎢───────────────────────⎥⎟\n",
       "⎜W, u, ⎢W₁₀⋅u₀₀ + W₁₁⋅u₁₀⎥, ⎢ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎥⎟\n",
       "⎜      ⎢                 ⎥  ⎢ℯ                   + 1⎥⎟\n",
       "⎜      ⎣W₂₀⋅u₀₀ + W₂₁⋅u₁₀⎦  ⎢                       ⎥⎟\n",
       "⎜                           ⎢           1           ⎥⎟\n",
       "⎜                           ⎢───────────────────────⎥⎟\n",
       "⎜                           ⎢ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎥⎟\n",
       "⎝                           ⎣ℯ                   + 1⎦⎠"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dim_u = 2\n",
    "dim_x = 3\n",
    "\n",
    "W = sy.MatrixSymbol(\"W\", dim_x,dim_u)\n",
    "u = sy.MatrixSymbol(\"u\", dim_u,1)\n",
    "linear  = lambda s: s\n",
    "sqr     = lambda s: s**2\n",
    "sigmoid = lambda s: 1/(1+sy.exp(-s))\n",
    "a = sigmoid\n",
    "s = sy.Matrix([i for i in W*u])\n",
    "x = sy.Matrix([a(i) for i in s])\n",
    "W, u, s, x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  the jacobian of the transformation\n",
    "\n",
    "$$\\mathbf{J} = J^i_\\alpha = \\frac{\\partial x^i}{\\partial u^\\alpha}$$\n",
    "\n",
    "in this case\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\frac{\\partial x^i}{\\partial u^\\alpha} =&  \\frac{\\partial \\text{a}(s^i)}{\\partial u^\\alpha} \\\\\n",
    " =& \\frac{\\partial \\text{a}(s^i)}{\\partial s^i} \\frac{\\partial s^i}{\\partial u^\\alpha}\\\\\n",
    " =& \\dot{\\text{a}}(s^i) W^i_\\alpha\\\\\n",
    "\\end{align*}$$\n",
    "\n",
    "denoting $\\frac{\\partial \\text{a}(s^i)}{\\partial s^i}=\\dot{\\text{a}}(s^i)$\n",
    "\n",
    "\n",
    "$$\\mathbf{J} = J^i_\\alpha = \\dot{\\text{a}}(s^i) W^i_\\alpha$$\n",
    "\n",
    "and observing that $ \\frac{\\partial s^i}{\\partial u^\\alpha} =  \\frac{\\partial}{\\partial u^\\alpha}(W^i_\\gamma u^\\gamma + b^i)=W^i_\\alpha$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} W_{0, 0}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{2}} & \\frac{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} W_{0, 1}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{2}}\\\\\\frac{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} W_{1, 0}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{2}} & \\frac{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} W_{1, 1}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{2}}\\\\\\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} W_{2, 0}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}} & \\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} W_{2, 1}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡  -W₀₀⋅u₀₀ - W₀₁⋅u₁₀          -W₀₀⋅u₀₀ - W₀₁⋅u₁₀      ⎤\n",
       "⎢ ℯ                  ⋅W₀₀     ℯ                  ⋅W₀₁  ⎥\n",
       "⎢──────────────────────────  ──────────────────────────⎥\n",
       "⎢                         2                           2⎥\n",
       "⎢⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞   ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞ ⎥\n",
       "⎢⎝ℯ                   + 1⎠   ⎝ℯ                   + 1⎠ ⎥\n",
       "⎢                                                      ⎥\n",
       "⎢  -W₁₀⋅u₀₀ - W₁₁⋅u₁₀          -W₁₀⋅u₀₀ - W₁₁⋅u₁₀      ⎥\n",
       "⎢ ℯ                  ⋅W₁₀     ℯ                  ⋅W₁₁  ⎥\n",
       "⎢──────────────────────────  ──────────────────────────⎥\n",
       "⎢                         2                           2⎥\n",
       "⎢⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞   ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞ ⎥\n",
       "⎢⎝ℯ                   + 1⎠   ⎝ℯ                   + 1⎠ ⎥\n",
       "⎢                                                      ⎥\n",
       "⎢  -W₂₀⋅u₀₀ - W₂₁⋅u₁₀          -W₂₀⋅u₀₀ - W₂₁⋅u₁₀      ⎥\n",
       "⎢ ℯ                  ⋅W₂₀     ℯ                  ⋅W₂₁  ⎥\n",
       "⎢──────────────────────────  ──────────────────────────⎥\n",
       "⎢                         2                           2⎥\n",
       "⎢⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞   ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞ ⎥\n",
       "⎣⎝ℯ                   + 1⎠   ⎝ℯ                   + 1⎠ ⎦"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J = sy.Matrix([[j.diff(i) for  i in u] for j in x])\n",
    "J"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "just as exercise, we can do this by hand if we use the sigmoid. Observe its differential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEQAAAAlCAYAAAD7u09NAAAABHNCSVQICAgIfAhkiAAAAztJREFUaIHt2E+IHEUUx/FPjBJXA2L0pHEXk5MSMZLgRVdFSSAKogdRwcNCxINgFPRictkkBxU8aEBYFCHoTS+SQ9CD//APKuouruhFRdBNFIPB/2CM8fB6SFHObDozXT27y36hYbq6qt/rN69eVf1WTE5OWgSswwYcKG3ojNIGGmIbLm/D0GIIyPXYi+2YxrkljZ1Z8uUN8TY+wwS+LW1soQRkRndftuIQRrUQDD2cKMHFeAy3YCVex/34sXq+cZ6xa0VQWqGNGnIpPsUcrsUNuBBTNceP4XARz7rQRkCm8DwexZdieuzFTTXHfy6W3VktrDSlp8yYqAPj2JG0r8SfNd/xCzY17FdPSgfkSvyq+wf9Xdh2X5QOyDGxb/gBvxe21Qila8gHOIoXcRXWYwueacF2X5R26qjYdp+HN0VBfRLf49/CtvuijX3Ix7ixBTuNsCDTdpgsByRjOSAZnRpyYqheLCA6AVkxVC8WEEttyqzDrfM8vwRv4QuhsdyRd1hqATmV1PgPHqr6bMVTMgUuD8j5QqNY35yPrdFNanwZDyd9DovNIXGcOII16UvygOzEQXxd3V8nlO45UXgnmvK+AB2pcYs4JvyBPdgldso5m8Sp+7u0MQ3IObhXaBcdVgs94kH81YeT+zHZx7huzFS+5NdFSZ9capzFN7gne9cavID7ciPp1v1mkQXvJW0Hq4v4uJIMIjPSW2o8gLvFgRJW4RU8jvfzzmmGjOMTw9mTDCoz0ltq/AhXY0RsL/bjDXEC/x9pQMa0KOZmDCoz0ltqPISzxNS6BnfitsrGDK5IX5JOmREn07NfdlZXh1Ui4x5J2rbhneS+CZmR3lJjp/aN4F2n2GqkATkilt1BmMJLyf0TYhrsS9rmsjGlZcbOsvpTnc5pQKYNvqz+XF0dfqvuv5pnTGmZcYP4E2plf5o+r+EyXJC0rRbVfWPVd7T6PdqEpxWlZcZx8W21SA3Oiop8V9K2WWTOtJiDu6vfewZ28yQlZcazcTueqzsglxB342lRC46Lg9AgJ+GJmv1KyYzb8aHIwlrkKfmqSNW1DTo1TI7hgdMZ0E1k3telbbHy7OkOWGrH/4FZDkjGf7Ahq+kg1KaPAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{e^{- t}}{\\left(1 + e^{- t}\\right)^{2}}$"
      ],
      "text/plain": [
       "    -t    \n",
       "   ℯ      \n",
       "──────────\n",
       "         2\n",
       "⎛     -t⎞ \n",
       "⎝1 + ℯ  ⎠ "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = sy.symbols('t')\n",
    "sigmoid(t).diff(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "substituting $t$ with $s^i$ to obtain $\\dot{\\text{a}}(s^i)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ4AAAAmCAYAAAA1HCJjAAAABHNCSVQICAgIfAhkiAAABdpJREFUeJztm2uIVVUUx39Ths2MZtq7kVERgsjKUPowZkSUYJFTBpZQYKMDEWRBfWmiGPVDBRJlBEOZiIJQkZQftBJ7P6bnSGP4oZdlvshymEmnGh99WPt49+zZ55597nnNvewfHObuc9ZdZ69z191r7z3/W9fZ2YnHkzdnFN0BB+qB/Vp7AXBCa7cAH+bao+TUYkyxqIbEGwTGau1lwCGgQWu/lHenElKLMcWiGhIP4CRQBzQBk4CvgPHquAl4o0K/3wIT1esW4NVk3YxFLcbkTLUk3t9AI3AfsE61xwOLgc3APxX4HANMAI6o9pVAb+KeulOLMTkzpugOKHZi78s8ZC7UjzzQO4E5wLXIh7QUaFO284GFwDlImdph+GoEXgSGgPfVPX/Qrs9Q7zHtNkX4Ma9HxRKQNKapwJNIeb5bnbssJCabbZyYUievEa8J2AD8CfQhZeQi7fpM5CGZR/BB9QO3A58Dx5DRoQU4DnyvbLYB7cD9wCJLHxYiZacdaFX+d2nXZyOjg2kX5cckKpaApDHtoZSgAWEx2WzjxJQ6eSTeNGTesQ+4DrgBOB/oiuGjH3iI0oR7QLVfttg+FuJ7MrBXvT6BzKv6VHsOcAXws8Uuyk+lpBGTSVhMUaQVkzN5lNou4BWgQzu3CpnHuNKPPNCdqj0AXAi8ZtitArYDPRYfvyMPeBfyhdsGvIWUp1+B3cApi12Un0pJIyaTsJiiSCsmZ+oy3kCeggzzg8gqLuBMpLycl+K92oDlQDcywnYjJX6but4IvAD8C3xM+DzGZneV5svVTxqUi+lc4GlkBdwFrC7jx2ZbVExA9om3ANgIzLJc+w/4LcN7L0Mm1r+MMl+jpR+FxpR1qR1Cvk0Hkclznqwdpb6SUDMxZV3Pu5E9pY3ANcB04GZk6V4te4ieDMj6wz+C7EVNoLR3thqZzJ4s8z5PjZPHqvZr4MYc7uOpIny58xSCTzxPIfjE8xRCMMdz2d32eFIjGPHqMjoagANau5WSDq0O+X/iRxne38dUzNGMKKh3IyKFRaZN1qvackrbY1Sn0rYWY0qb48DDyPbZxcA3wFbgaGBgzvEmIg9xeoqdyEppC8WpbWsxpiS8DjyitQ9QEj8cBA4jz+k0ZuJ1IJn5U4qdykJpC8WqbWsxpjCuB7YgsrZTwBKLzUrgcaTvJrMQUche/aReahuQMnFbzI6loR6uhDC1bRLSVA9XQhYx2ViPqIY6HWzHIXKpDeqw0Yvo/u5B/h0aMEm9p918g/6Qb0Ey+lPDpgl4CrgVydwdwANISQZR3JbDRWkbV7a+Cbva9rkQ2yhfAVGxpBnTVEbK0cNistm6xpSUreoASdgwtiAjfpB4Y4E3ETnWZ6axXmrnIpNAfWtlGvmoh+PK1iFcbZtUuu5CGjHtYeToGBaTzVYnd+m6hS+Rkb8emf+uB95DBCIj0Ee8KYz8XUCe6mEoL1sPFLiBNDtMbWuzjfIVlzRislGpgjiNmJKyHzgLuBS4BLgL+A6pDAD3os1X9cSrp1Q+QRJxHjISLtfOB+phVxYb7WfVYRJHtg7yweil8YkytlG+4pJGTDbCYooiKqYOhg8eY5GEflQ7Nx9RH1fKoPpbD3wS0o/T6Il3mNIyHuBq5Jsdph5OkzZkUXMBUt5N2fpmRJrdiswlymHa6hLvuL6SEBVTIEefjSSAi3Q9sH2XeDF1MXw0fgaZPq3Rzu1zCyuUYLvkDxdjPfF6GL5UzlM9vE4dAYEsO+Ao7qtF0zaJryRExdSHzP9cMG3jxvSXOgIGVPtHx/u7MANJ3kNRhjB8OHwHuJzSD3CKVA+vJb3fAqTpKwnVGtM4pPzPRD73ZvW62bCbi+SQE3oC9SIrk2C57tXDHpDy3qOOemCFer1SszkbuAP7b4KtmJulK4DnkTnBCbx6uFZZEsP2A2R7pBxLgS+QKumEWTLfRkrp5Bgd83iGgAfjvMH276E1lnMeTzliq3G8AtlTCD7xPIXgE89TCP8DOhQUDQcPJ9EAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}}$"
      ],
      "text/plain": [
       "    -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    \n",
       "   ℯ                      \n",
       "──────────────────────────\n",
       "                         2\n",
       "⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞ \n",
       "⎝ℯ                   + 1⎠ "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigmoid(t).diff(t).subs({t: s[dim_x-1]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "to we can build the Jacobian by hand for the sigmoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} W_{0, 0}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{2}} & \\frac{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} W_{0, 1}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{2}}\\\\\\frac{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} W_{1, 0}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{2}} & \\frac{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} W_{1, 1}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{2}}\\\\\\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} W_{2, 0}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}} & \\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} W_{2, 1}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡  -W₀₀⋅u₀₀ - W₀₁⋅u₁₀          -W₀₀⋅u₀₀ - W₀₁⋅u₁₀      ⎤\n",
       "⎢ ℯ                  ⋅W₀₀     ℯ                  ⋅W₀₁  ⎥\n",
       "⎢──────────────────────────  ──────────────────────────⎥\n",
       "⎢                         2                           2⎥\n",
       "⎢⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞   ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞ ⎥\n",
       "⎢⎝ℯ                   + 1⎠   ⎝ℯ                   + 1⎠ ⎥\n",
       "⎢                                                      ⎥\n",
       "⎢  -W₁₀⋅u₀₀ - W₁₁⋅u₁₀          -W₁₀⋅u₀₀ - W₁₁⋅u₁₀      ⎥\n",
       "⎢ ℯ                  ⋅W₁₀     ℯ                  ⋅W₁₁  ⎥\n",
       "⎢──────────────────────────  ──────────────────────────⎥\n",
       "⎢                         2                           2⎥\n",
       "⎢⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞   ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞ ⎥\n",
       "⎢⎝ℯ                   + 1⎠   ⎝ℯ                   + 1⎠ ⎥\n",
       "⎢                                                      ⎥\n",
       "⎢  -W₂₀⋅u₀₀ - W₂₁⋅u₁₀          -W₂₀⋅u₀₀ - W₂₁⋅u₁₀      ⎥\n",
       "⎢ ℯ                  ⋅W₂₀     ℯ                  ⋅W₂₁  ⎥\n",
       "⎢──────────────────────────  ──────────────────────────⎥\n",
       "⎢                         2                           2⎥\n",
       "⎢⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞   ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞ ⎥\n",
       "⎣⎝ℯ                   + 1⎠   ⎝ℯ                   + 1⎠ ⎦"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J_hand = sy.Matrix([[sigmoid(t).diff(t).subs({t: s[i]})*W[i,alpha] for i in range(dim_x)] for alpha in range(dim_u)]).T\n",
    "J_hand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J==J_hand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "get numerical evaluations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sy_eval(expr, W, u, nW, nu):\n",
    "    return np.r_[expr.subs({W[i,alpha]: nW[i, alpha] for i,alpha in itertools.product(range(0,dim_x), range(0,dim_u))})\\\n",
    "                     .subs({u[alpha,0]: nu[alpha,0] for alpha in range(0,dim_u)})].astype(float)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.19254481],\n",
       "       [0.2680553 ],\n",
       "       [0.86977336]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sy_eval(x, W, u, nW=np.random.normal(size=W.shape), nu=np.random.normal(size=u.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.04185382, -0.13996193],\n",
       "       [ 0.15571656,  0.05539929],\n",
       "       [-0.06997518, -0.18963061]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sy_eval(J, W, u, nW=np.random.normal(size=W.shape), nu=np.random.normal(size=u.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### now in TensorFlow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Variable 'W:0' shape=(3, 2) dtype=float64, numpy=\n",
       " array([[ 0.7329587 ,  0.12032583],\n",
       "        [ 1.06797699, -0.4415987 ],\n",
       "        [-2.04648562,  1.4915178 ]])>,\n",
       " <tf.Variable 'Variable:0' shape=(2, 1) dtype=float64, numpy=\n",
       " array([[-0.81454065],\n",
       "        [ 0.17518983]])>)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tW = tf.Variable(np.random.normal(size=(dim_x,dim_u)), name=\"W\")\n",
    "tu = tf.Variable(np.random.normal(size=(dim_u,1)))\n",
    "tb  = tf.Variable(np.zeros((dim_x,1)))\n",
    "tW, tu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transformation(tW, tu):\n",
    "    \n",
    "    tsigmoid = lambda x: 1/(1+tf.exp(-x))\n",
    "    tlinear = lambda x: x\n",
    "    ta = tsigmoid\n",
    "    ts = tf.matmul(tW,tu) + tb\n",
    "    tx = ta(ts)\n",
    "    return tx, ts, ta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "tx, ts, ta = transformation(tW, tu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "check transformation function evals ok in sympy and tensorflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(3, 1), dtype=float64, numpy=\n",
       "array([[0.35986623],\n",
       "       [0.27943332],\n",
       "       [0.87305498]])>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.35986623],\n",
       "       [0.27943332],\n",
       "       [0.87305498]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sx = sy_eval(x, W, u, nW=tW.numpy(), nu=tu.numpy())\n",
    "sx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(tx, sx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "compute the Jacobian and check it with sympy's version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(3, 2), dtype=float64, numpy=\n",
       "array([[ 0.16884622,  0.02771856],\n",
       "       [ 0.21503753, -0.08891605],\n",
       "       [-0.22681196,  0.16530489]])>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_jacobian(tW, tu, transformation):\n",
    "    with tf.GradientTape(persistent=True) as t:\n",
    "        tx, ts, ta = transformation(tW, tu)\n",
    "    tJ = t.jacobian(tx, tu)[:,0,:,0]\n",
    "    return tJ\n",
    "tJ = get_jacobian(tW, tu, transformation)\n",
    "tJ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.16884622,  0.02771856],\n",
       "       [ 0.21503753, -0.08891605],\n",
       "       [-0.22681196,  0.16530489]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sJ = sy_eval(J, W, u, nW=tW.numpy(), nu=tu.numpy())\n",
    "sJ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(sJ, tJ.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## the metric tensor of the embeded surface or the projection produced by the transformation\n",
    "\n",
    "$$g_{\\alpha\\beta} = \\sum_{k=0}^2 J^k_\\alpha J^k_\\beta = \\mathbf{J}^T\\mathbf{J} $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "observe that we must keep the sumation symbol since Einstein's convention does not kick in automatically over the $k$ index, due to the way the matrix multiplication with the transpose must be expressed:\n",
    "$$\n",
    "\\Bigg[\\begin{matrix}\n",
    "J^i_\\mu & J^j_\\mu & J^l_\\mu\\\\\n",
    "J^i_\\nu & J^j_\\nu & J^l_\\nu\\\\\n",
    "\\end{matrix}\\Bigg]\n",
    "\\times\n",
    "\\Bigg[\\begin{matrix}\n",
    "J^i_\\mu & J^i_\\nu \\\\\n",
    "J^j_\\mu & J^j_\\nu \\\\\n",
    "J^l_\\mu & J^l_\\nu \\\\\n",
    "\\end{matrix}\\Bigg]\n",
    "=\n",
    "\\Bigg[\\begin{matrix}\n",
    "J^i_\\mu J^i_\\mu+J^j_\\mu J^j_\\mu+J^l_\\mu J^l_\\mu & J^i_\\mu J^i_\\nu+J^j_\\mu J^j_\\nu+J^l_\\mu J^l_\\nu \\\\\n",
    "J^i_\\nu J^i_\\mu+J^j_\\nu J^j_\\mu+J^l_\\nu J^l_\\mu & J^i_\\nu J^i_\\nu+J^j_\\nu J^j_\\nu+J^l_\\nu J^l_\\nu \\\\\n",
    "\\end{matrix}\\Bigg]\n",
    "=\n",
    "\\sum_{k=0}^2 J^k_\\alpha J^k_\\beta\n",
    "$$\n",
    "\n",
    "for $\\alpha, \\beta \\in \\{\\mu, \\nu\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "in this case\n",
    "\n",
    "\n",
    "$$g_{\\alpha\\beta} = \\sum_{k=0}^2 \\dot{\\text{a}}(s^k)^2 W^k_\\alpha W^k_\\beta$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAG0AAAAZCAYAAAA7S6CBAAAABHNCSVQICAgIfAhkiAAAAvRJREFUaIHt2k+IVlUYx/HPjKahQ6QO+GehFQUtLP+NG0EZkNnMMhS3764yImiXCL4T6FJwpwhqtBKEIFHRjaOgC0mDEKtN/yDLGjQIVLScFue+zZ07V+bce3x75433C5dz33vOc875nec5555779vXbDb16C76c+fDmMwd33SiQz3+ZdB0f0y2MuaXFL6EcUz8Fz3r8VTuYyw7b2BNK6PMaeNotrtHPWblvik/DMs5rb+kcI85ThWnfSSsq2+V5K3O8j57Fp3qAF2lrYrTNmXpFyV5m7P0Rlp3OkZXaavqtAn8VJI3lKVfJveoM3SVtlinLcVLuP6U/JawORONFeg6bbFO25ilZcsHIVJ/w+3C9d34Hg+FQdlasX+p9jHU0bYNn+Nn4X7XqNFubW2xTmut+WXR+AqWmBmJu3AIB7ABV3FOuLHHkGofSx1tA7iJD/CgRptJ2mKdtiFLy5aI0Swtrvkf4gSO4mu8j1/wbmSbqfax1NF2FntwCk9qtJmkLdZpr+MxfixcX4i3s/O86AVCBF8olL+ALRHtpdpXoaq2VJK1xTrtEZ7Da7lri3EMa7Pf+WgcxDzcKdRzBysi2qtqf0L9e0tVbamkjk3pa6wyzgvPK5eFh8wBbMdXwrRehO8i62oHreD7q4btXNc2g1in7ccL2ClE8y18jJO4KwiezJWfwN9YXqhnOX6NaK+q/Rv4E2ci6i5SVVsqqWMTvTw+FHZKq4TIG8IR/JHVMVwo/0jYjY0Uro8IO6XZqGL/It7EYdyLqLtIVW2ppI5N6Uzblx3fCjfpuhzEp7iGK3hHGJjDz9h+q7CROJjQ16oM4NXsvF/Yqq8XZmbZW5UiMdoG8XuZcd5pP5j6fkP697STWIa9WCk814yavktr4Dheztqvag+n8XxiX6syhIu532PZ8YmpzVBDmrb897Rp9HX47wZj2IF16m0i5jJt09bp72mjeM//z2G0UVvs7rFdbJ69SNfSNm2dnmk9atBzWhfSc1oX8g9qzN0GF9L/TAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left[ u_{0, 0}, \\  u_{1, 0}\\right]$"
      ],
      "text/plain": [
       "[u₀₀, u₁₀]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 0}^{2}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 0}^{2}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 0}^{2}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{4}} & \\frac{e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 0} W_{2, 1}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 0} W_{1, 1}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 0} W_{0, 1}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{4}}\\\\\\frac{e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 0} W_{2, 1}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 0} W_{1, 1}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 0} W_{0, 1}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{4}} & \\frac{e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 1}^{2}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 1}^{2}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 1}^{2}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{4}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡     -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀    2    -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁⋅u₁₀    2    -2⋅W₀₀⋅u₀₀\n",
       "⎢    ℯ                      ⋅W₂₀    ℯ                      ⋅W₁₀    ℯ          \n",
       "⎢    ──────────────────────────── + ──────────────────────────── + ───────────\n",
       "⎢                              4                              4               \n",
       "⎢     ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞      ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞      ⎛ -W₀₀⋅u₀₀\n",
       "⎢     ⎝ℯ                   + 1⎠      ⎝ℯ                   + 1⎠      ⎝ℯ        \n",
       "⎢                                                                             \n",
       "⎢ -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀            -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁⋅u₁₀            -2⋅W₀₀⋅u\n",
       "⎢ℯ                      ⋅W₂₀⋅W₂₁   ℯ                      ⋅W₁₀⋅W₁₁   ℯ        \n",
       "⎢─────────────────────────────── + ─────────────────────────────── + ─────────\n",
       "⎢                            4                                 4              \n",
       "⎢   ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞         ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞         ⎛ -W₀₀\n",
       "⎣   ⎝ℯ                   + 1⎠         ⎝ℯ                   + 1⎠         ⎝ℯ    \n",
       "\n",
       " - 2⋅W₀₁⋅u₁₀    2        -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀            -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁⋅\n",
       "            ⋅W₀₀        ℯ                      ⋅W₂₀⋅W₂₁   ℯ                   \n",
       "─────────────────       ─────────────────────────────── + ────────────────────\n",
       "               4                                    4                         \n",
       " - W₀₁⋅u₁₀    ⎞            ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞         ⎛ -W₁₀⋅u₀₀ - W₁₁⋅\n",
       "           + 1⎠            ⎝ℯ                   + 1⎠         ⎝ℯ               \n",
       "                                                                              \n",
       "₀₀ - 2⋅W₀₁⋅u₁₀               -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀    2    -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁\n",
       "              ⋅W₀₀⋅W₀₁      ℯ                      ⋅W₂₁    ℯ                  \n",
       "──────────────────────      ──────────────────────────── + ───────────────────\n",
       "                   4                                  4                       \n",
       "⋅u₀₀ - W₀₁⋅u₁₀    ⎞          ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞      ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u\n",
       "               + 1⎠          ⎝ℯ                   + 1⎠      ⎝ℯ                \n",
       "\n",
       "u₁₀            -2⋅W₀₀⋅u₀₀ - 2⋅W₀₁⋅u₁₀        ⎤\n",
       "   ⋅W₁₀⋅W₁₁   ℯ                      ⋅W₀₀⋅W₀₁⎥\n",
       "─────────── + ───────────────────────────────⎥\n",
       "        4                                 4  ⎥\n",
       "u₁₀    ⎞         ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞   ⎥\n",
       "    + 1⎠         ⎝ℯ                   + 1⎠   ⎥\n",
       "                                             ⎥\n",
       "⋅u₁₀    2    -2⋅W₀₀⋅u₀₀ - 2⋅W₀₁⋅u₁₀    2     ⎥\n",
       "    ⋅W₁₁    ℯ                      ⋅W₀₁      ⎥\n",
       "───────── + ────────────────────────────     ⎥\n",
       "       4                              4      ⎥\n",
       "₁₀    ⎞      ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞       ⎥\n",
       "   + 1⎠      ⎝ℯ                   + 1⎠       ⎦"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mt = MetricTensor(sy.Array(J.T*J), list(u))\n",
    "mt.tensor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAArCAYAAAAUsNCoAAAABHNCSVQICAgIfAhkiAAADAFJREFUeJztnW2MFdUZgJ8VGirrFkQL1hVYo5BSJYVCsAEl0KKRxgDFaqGxKSCkJm3VWv9AE7qVP9pgQmuIRKylpTWiYqtt2KCl2NpWIhVIWWtara6VRW2qwO5WBKH0x3sme3Z2Pu7Mnbkz99z3SW52Z+bcM+e98+zl5ZwzZ5ra29tRFEVRFEWJ4qwM6xoLPAv8DfgrcIPZfzZw2Cq3ADhtbc8Efp9hO2pBI8VqExa36zTS9W6kWG3UbfevdyPFapOZ20MzahDAKeB24ABwAfAisAP4LzDMKrcSeAcYDrxvth/IsB21oJFitYmK22Ua6Xo3Uqw26rb717uRYrXJzO0sexjeMg0CeBv4DzDKbP8PaAJazb69QIt5zQO2pzznPuBc8/tMYFvKepJSRKxQXLweUXG7jLotqNvuoW4L6nYFZJkw2EwDhgBvmu0+oBlYDjxktluApcATwAcpzjEUGAEcMduTgYPpm5yaWsQK2ce7Bejw7fsMcAa437d/KXAcGG3t88fdKKjb6rarqNvqdiRJhyQOhLznGvrHgEYBPwNWWcd7kA/temAWMAO5GDcDK0yZ+cBi4GNI988u3zmagY3Ah8Bu05ZXrOOXW+/xl304oh77WKUxeuQVaxuwFukSW2L2TQyJN6isTVi8R02dNncCJ4CRvv13AFuBf0fEXe/k6XYbya5RPbsdF2vQcXU7X/J0O845l9yOa0Mbjrvt72FoNRW+axq2HRhjHZ+CBOx/eRdkGPAr4G7gz9b7eoBFwPPImFAf0jVzCnjJlOkwgdwC3BjQ1sVIV84qYKE5b6d1fDr9mZu/bFQ9fuJi9Mgz1i76JfUIizeorE1YvEeQPxCPccgfy2bf/tlIVnqv2Q6Lu8zEeQ35ut1FsmtUz27HxRp0XN1OT9FuxznnkttxbejCcbfthOFiZKylG7gSmAOcD2yqsK4mpLvkd0hWY9MD3Eb/xJFes705oJ7VIee8iP5ulNNItnTUbM8CLgNeCykbVU8a8o41iKh4owiL9ygDBbsdmQizl4GZ6neA3wB/JzruslKt15Dd9Q7DRbeToG6nowxuxznnkttp2uCU23YXzibgx8Aaa986ZPymEmYBX0Zu21hk9n0VyaZ6kEC9iRe9yLjKo7461gHPAPsD6j+EfIidSKLTATyJdM+8AbyMjOMElY2qJw15xxpEVLxRhMVrZ6ojkJnA1wIft/ZPAK4D5prtqLjLSrVeQzbXOwrX3E6Kup2OMrgd55xLbqdpg1NuN5mFm8YjXSTHkdmiHkOQ7pnzklacghXArcAeJGveg3S3eRM8moH7kLGa5xg8fmTjL9tp1ZWknryIi3Uk0nU0D/lSWB9Rl7/s01T2uS1ARB6KjHUtQsSaC/zc1LER6UK7oopYi6QMXkP6axRE2d2Oi1XdzoayuB3lY9DxenY7Llbn3fYShgVIV8W0gDIngX9lfeIKWIlMDnm9ZHXlQRGxzkYWIxmD3Jf7LWSMa5rZPw7pEltOf5Y9FvFkNDK2tw54LOIc1yFjaGcB9wAPJo6mOsroNajbedeVxm2AXyLd+ruAL8WcQ90ORt3Ot640bid1NbS8lzDMB36NZD198bEpDjAZ6aJaC9wETEL+p3IpMqv3e8AypHvLG0P7BCKqvQDIRIIXABmKrCw2Fzhmys5EJmfVCvW6MUnjNkiy0AJ8jeiEQd1WiiKp20ldjSzvjY3sQcZGtgJTgUuAq5GujbzWalCKxZuIcweSTXrdmj3m57eBDQz8Qk2yAMgMZHZxN/KF1oHc4lRL1OvGJI3bIMvn9lZQv7qtFEVSt5O6GlneE+sIkrGOoP9e2fXIxAt7fExxB28hkZPIbVkex8zPM8iEqjDiFgC5EJHOoxsZX6sl6nVjUq3bcajbSlEkdTupq5Hl7bsk/gJ8rqImKy7Qh9xu4+dEyH6belrcRr1uPKpxu55QtxuPQt3WrislKZUuAHKYgZlsK4MXUUnCMiR7nlNFHYqSBeq2Ui8kdTWyvCYMShKiFgDZxUDRXkBWOWsFzkG6T3fm30RFyRx1W6lX4lxN5HaWj7dW3CdsAZBOZJbue1bZU8hqY7uRxPQH1HYWuaIk5bfAp5F74A8BNyCTC9VtpV6JcrWJhG5rwqAk4Y8E90pNQtawP+7b/5R5KUo9MC9gn7qt1Dthrn6ShG7rkISSBS8jt/koimuo24qrJHZbEwaljHQhk8Ds10/Msd0Bx7bUvIWKko4u1G2lTtEhCaWMbGDws92nII95/SnypWtzAEWpD9RtpW7RhEEpIxsC9i1DvlS3ICvyKUo9om4rdYuXMFTyuE1FyZJaLKCjXitFoG4rTuLNYWjK6TUcef6At70QWbbU254F/CHH82tMxb6akafm3RtwrBaoBxpTXi91u/hr0Ogx5fUKdTvvSY/HkZUBPVYC7yAXz9t+IOc2ZI2LMeXFd5H72F3ERQ9cjCkv1O36wsWY8iLUbX/CcC7yIV6S4cm9LK4VeQbBXuQRsi3Ifc/bq6h7H9JmkEdwbquiriS4GFM1PIYs9mEzAbnPt6P2zRlEHl6Dmx64GFM1qNvueOBiTNWQ2G1/wrAG2AH8M8NG9SFdHMuBh8x2C7AUeAL4IGW9Q5EntXlP75oMHKyqpZXjYkxhzEYW8ehGxk2XBZS5C8lKR1j71gOrM2zHFuSP/dkU783Da3DTAxdjCkPdDsdFD1yMKYxc3LYThuFIt0zSx74eQJYG9r8uNMd7TIOuBx5BnjffAtwMbDZl5pvftwGfDzhHM/JHtRn4itk3EXjFKnO5OW9Q2bi6Ko3FI4uY2hBpH7H2hcUUVLbSmILYArRXUA5kPfFO4DYGrwjmcRB4DbjJbC8E/mFeRZPWa8jG7TaSX7uyux0XU9BxdTt7inY7jY9ldzuuHW00sNt2wvAFJBP5k69MK/Io43eBo0i3zRjr+BTkA/K/vCdc9SDPHXgeeB/J6mYia1a/ZMp0II9KvgW4MaCdi5F/eFeZoKD/onhMRz6AoLJxdVUai0cWMXUBK3z7wmIKKltpTNWyA/lfzONIl14YTyGZOsBngSVIu9ebdq3NuF2VEuY11MbtLpJfu7K7HRdT0HF1O3uKdjuNj2V3O64dXTSw23bCcBXwIgNv17kYGZvpBq5EHr96PrApQcN7kCzHm1DSa7Y3B5RdHVL3RcCb5vfT5uco5A8BZIbrZUi2FFQ2rq6kZBFTEGExxZFFTNXyAjADOBuJeSySYd+JfC53FdSuIK+h9m6HUc9uJ0Hdzp6i3U7jY9ndTnOehnHbXrhpPIMzsk1Id9caa986ZLynUnqQD9NbsawXGA086iu3DngG2B9QxyHkw+2kP8npAJ5EuoPeQNbFPhNSNq6upGQRUxBhMcWRRUzVchj4CNINmPV4ajUEeQ21dTuKenU7Kep29hTtdhofy+52mvM0jNtN7e3t3u87ka6Ir5vt8Wb7OAO7NIYg3TnnZdBYjxXArcitHPvMz1b6Z2o2A/cBJ4DngIcj6vKX7ayirmqIi2kkcDcyO3cT0gUUhr/s0ySLaQ0DvzyGIUKftPbNN++Nog/4JsHr209Axr4mM7B7rmj8XkNt3a722tmUxe24mNTt2lC021n6WBa349rR0G7bCcMvkEzDG29fAGwFpgW87ySysENerAR2Aa+XrK6ytCNpXaPMy+MepLvyR9a+bsInx3hEiXcFkhRdgNzmVRb8XoO6nTXqdjGo2/mjblvYQxL7GXjrxYdIBvS2OWEtebCkdVVDkTG9Z14evWb71cxaJBN/uinXFyoM9hrU7axRt4tB3c4fddvCHjfZCUyiv8tqD3Jf6VZgKrIwyNXARvSx2I3EOcgs5CnIdR9nfh/nK3cV4lDZ8HsN6rYiqNuKq+Titi3QQWTG5BKzfQQZHxmBPKf9ADJec4jo2zQUt5iO/E9mPzKT9vvmd3v27EeBL1L9TPo88HsN6rYiqNuKq+Titj2HAeBa4IfApyjuNg+l/vgGch/xNUU3JAT1WkmLuq24SmK3h8yZM8fefhXpdXgLOJZlyxSnmQrcjywSU0bUayUt6rbiKond9vcwKIqiKIqiDEInwSiKoiiKEosmDIqiKIqixKIJg6IoiqIosfwf7b79h8h+BBUAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 0}^{2}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 0}^{2}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{4}} + \\frac{e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 0}^{2}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{4}}$"
      ],
      "text/plain": [
       " -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀    2    -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁⋅u₁₀    2    -2⋅W₀₀⋅u₀₀ - 2⋅\n",
       "ℯ                      ⋅W₂₀    ℯ                      ⋅W₁₀    ℯ               \n",
       "──────────────────────────── + ──────────────────────────── + ────────────────\n",
       "                          4                              4                    \n",
       " ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞      ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞      ⎛ -W₀₀⋅u₀₀ - W₀\n",
       " ⎝ℯ                   + 1⎠      ⎝ℯ                   + 1⎠      ⎝ℯ             \n",
       "\n",
       "W₀₁⋅u₁₀    2\n",
       "       ⋅W₀₀ \n",
       "────────────\n",
       "          4 \n",
       "₁⋅u₁₀    ⎞  \n",
       "      + 1⎠  "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mt[0,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metric tensor from sympy's Jacbians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.12619385, -0.05193324],\n",
       "       [-0.05193324,  0.03600009]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sG = np.r_[[sy_eval(mt[alpha, beta], W, u, nW=tW.numpy(), nu=tu.numpy()) for alpha, beta in itertools.product(range(dim_u), range(dim_u))]].reshape(dim_u,dim_u)\n",
    "sG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metric tensor from tensorflow Jacobians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_metric_tensor(tW, tu, transformation):\n",
    "    with tf.GradientTape(persistent=True) as t:\n",
    "        tx, ts, ta = transformation(tW, tu)\n",
    "    tJ = t.jacobian(tx, tu)[:,0,:,0]\n",
    "    tG = tf.matmul(tf.transpose(tJ), tJ)\n",
    "\n",
    "    return tG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(2, 2), dtype=float64, numpy=\n",
       "array([[ 0.12619385, -0.05193324],\n",
       "       [-0.05193324,  0.03600009]])>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tG = get_metric_tensor(tW, tu, transformation)\n",
    "tG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metric tensor by directly applying $g_{\\alpha\\beta} = \\sum_{k=0}^2 \\dot{\\text{a}}(s^k)^2 W^k_\\alpha W^k_\\beta$ in tensorflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_metric_tensor_PERCEPTRON(tW, tu, transformation):\n",
    "    with tf.GradientTape(persistent=True) as t:\n",
    "        tx, ts, ta = transformation(tW, tu)\n",
    "\n",
    "    mG = np.r_[[[tf.reduce_sum(t.gradient(tx, ts)**2 * tW[:,alpha:alpha+1] * tW[:,beta:beta+1]).numpy() for alpha in range(dim_u)] for beta in range(dim_u)]]\n",
    "\n",
    "    return mG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.12619385, -0.05193324],\n",
       "       [-0.05193324,  0.03600009]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mG = get_metric_tensor_PERCEPTRON(tW, tu, transformation)\n",
    "mG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(True, True)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(sG, tG.numpy()), np.allclose(tG.numpy(), mG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### the inverse metric tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(2, 2), dtype=float64, numpy=\n",
       "array([[19.50242854, 28.13393853],\n",
       "       [28.13393853, 68.36334741]])>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tG_ = tf.linalg.pinv(tG)\n",
    "tG_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "check they are inverses from each other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(tG.numpy().dot(tG_.numpy()), np.eye(tG_.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## the Hessian\n",
    "\n",
    "$$H^i_{\\alpha \\beta} = \\frac{\\partial x^i}{\\partial u^\\alpha \\partial u^\\beta} = \\frac{\\partial}{\\partial u^\\alpha}J^i_\\beta$$\n",
    "\n",
    "In this case:\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\frac{\\partial x^i}{\\partial u^\\alpha \\partial u^\\beta} =& \\frac{\\partial }{\\partial u^\\alpha}\\Big[\\dot{\\text{a}}(s^i) W^i_\\beta\\Big] \\\\\n",
    " =& \\frac{\\partial }{\\partial u^\\alpha}\\Big[\\dot{\\text{a}}(s^i)\\Big] W^i_\\beta + \\dot{\\text{a}}(s^i)\\frac{\\partial }{\\partial u^\\alpha}\\Big[W^i_\\beta\\Big]\\\\\n",
    " =&\\frac{\\partial\\dot{\\text{a}}(s^i) }{\\partial s^i}\\frac{\\partial s^i}{\\partial u_\\alpha} W^i_\\beta\\\\\n",
    " =&\\ddot{\\text{a}}(s^i) W^i_\\alpha W^i_\\beta \n",
    "\\end{align*}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\left[ - \\frac{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} W_{0, 0}^{2}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 0}^{2}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{3}}, \\  - \\frac{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} W_{0, 0} W_{0, 1}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 0} W_{0, 1}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{3}}\\right] & \\left[ - \\frac{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} W_{0, 0} W_{0, 1}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 0} W_{0, 1}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{3}}, \\  - \\frac{e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} W_{0, 1}^{2}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{0, 0} u_{0, 0} - 2 W_{0, 1} u_{1, 0}} W_{0, 1}^{2}}{\\left(e^{- W_{0, 0} u_{0, 0} - W_{0, 1} u_{1, 0}} + 1\\right)^{3}}\\right]\\\\\\left[ - \\frac{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} W_{1, 0}^{2}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 0}^{2}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{3}}, \\  - \\frac{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} W_{1, 0} W_{1, 1}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 0} W_{1, 1}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{3}}\\right] & \\left[ - \\frac{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} W_{1, 0} W_{1, 1}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 0} W_{1, 1}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{3}}, \\  - \\frac{e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} W_{1, 1}^{2}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{1, 0} u_{0, 0} - 2 W_{1, 1} u_{1, 0}} W_{1, 1}^{2}}{\\left(e^{- W_{1, 0} u_{0, 0} - W_{1, 1} u_{1, 0}} + 1\\right)^{3}}\\right]\\\\\\left[ - \\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} W_{2, 0}^{2}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 0}^{2}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{3}}, \\  - \\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} W_{2, 0} W_{2, 1}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 0} W_{2, 1}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{3}}\\right] & \\left[ - \\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} W_{2, 0} W_{2, 1}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 0} W_{2, 1}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{3}}, \\  - \\frac{e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} W_{2, 1}^{2}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{2}} + \\frac{2 e^{- 2 W_{2, 0} u_{0, 0} - 2 W_{2, 1} u_{1, 0}} W_{2, 1}^{2}}{\\left(e^{- W_{2, 0} u_{0, 0} - W_{2, 1} u_{1, 0}} + 1\\right)^{3}}\\right]\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡⎡    -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    2       -2⋅W₀₀⋅u₀₀ - 2⋅W₀₁⋅u₁₀    2     -W₀₀⋅u₀₀ -\n",
       "⎢⎢   ℯ                  ⋅W₀₀     2⋅ℯ                      ⋅W₀₀     ℯ          \n",
       "⎢⎢- ────────────────────────── + ──────────────────────────────, - ───────────\n",
       "⎢⎢                           2                              3                 \n",
       "⎢⎢  ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞      ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞        ⎛ -W₀₀⋅u₀₀\n",
       "⎢⎣  ⎝ℯ                   + 1⎠      ⎝ℯ                   + 1⎠        ⎝ℯ        \n",
       "⎢                                                                             \n",
       "⎢⎡    -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    2       -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁⋅u₁₀    2     -W₁₀⋅u₀₀ -\n",
       "⎢⎢   ℯ                  ⋅W₁₀     2⋅ℯ                      ⋅W₁₀     ℯ          \n",
       "⎢⎢- ────────────────────────── + ──────────────────────────────, - ───────────\n",
       "⎢⎢                           2                              3                 \n",
       "⎢⎢  ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞      ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞        ⎛ -W₁₀⋅u₀₀\n",
       "⎢⎣  ⎝ℯ                   + 1⎠      ⎝ℯ                   + 1⎠        ⎝ℯ        \n",
       "⎢                                                                             \n",
       "⎢⎡    -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    2       -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀    2     -W₂₀⋅u₀₀ -\n",
       "⎢⎢   ℯ                  ⋅W₂₀     2⋅ℯ                      ⋅W₂₀     ℯ          \n",
       "⎢⎢- ────────────────────────── + ──────────────────────────────, - ───────────\n",
       "⎢⎢                           2                              3                 \n",
       "⎢⎢  ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞      ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞        ⎛ -W₂₀⋅u₀₀\n",
       "⎣⎣  ⎝ℯ                   + 1⎠      ⎝ℯ                   + 1⎠        ⎝ℯ        \n",
       "\n",
       " W₀₁⋅u₁₀              -2⋅W₀₀⋅u₀₀ - 2⋅W₀₁⋅u₁₀        ⎤  ⎡   -W₀₀⋅u₀₀ - W₀₁⋅u₁₀ \n",
       "        ⋅W₀₀⋅W₀₁   2⋅ℯ                      ⋅W₀₀⋅W₀₁⎥  ⎢  ℯ                  ⋅\n",
       "──────────────── + ─────────────────────────────────⎥  ⎢- ────────────────────\n",
       "               2                                3   ⎥  ⎢                      \n",
       " - W₀₁⋅u₁₀    ⎞        ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞    ⎥  ⎢   ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁\n",
       "           + 1⎠        ⎝ℯ                   + 1⎠    ⎦  ⎣   ⎝ℯ                 \n",
       "                                                                              \n",
       " W₁₁⋅u₁₀              -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁⋅u₁₀        ⎤  ⎡   -W₁₀⋅u₀₀ - W₁₁⋅u₁₀ \n",
       "        ⋅W₁₀⋅W₁₁   2⋅ℯ                      ⋅W₁₀⋅W₁₁⎥  ⎢  ℯ                  ⋅\n",
       "──────────────── + ─────────────────────────────────⎥  ⎢- ────────────────────\n",
       "               2                                3   ⎥  ⎢                      \n",
       " - W₁₁⋅u₁₀    ⎞        ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞    ⎥  ⎢   ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁\n",
       "           + 1⎠        ⎝ℯ                   + 1⎠    ⎦  ⎣   ⎝ℯ                 \n",
       "                                                                              \n",
       " W₂₁⋅u₁₀              -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀        ⎤  ⎡   -W₂₀⋅u₀₀ - W₂₁⋅u₁₀ \n",
       "        ⋅W₂₀⋅W₂₁   2⋅ℯ                      ⋅W₂₀⋅W₂₁⎥  ⎢  ℯ                  ⋅\n",
       "──────────────── + ─────────────────────────────────⎥  ⎢- ────────────────────\n",
       "               2                                3   ⎥  ⎢                      \n",
       " - W₂₁⋅u₁₀    ⎞        ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞    ⎥  ⎢   ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁\n",
       "           + 1⎠        ⎝ℯ                   + 1⎠    ⎦  ⎣   ⎝ℯ                 \n",
       "\n",
       "             -2⋅W₀₀⋅u₀₀ - 2⋅W₀₁⋅u₁₀              -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    2      \n",
       "W₀₀⋅W₀₁   2⋅ℯ                      ⋅W₀₀⋅W₀₁     ℯ                  ⋅W₀₁     2⋅\n",
       "─────── + ─────────────────────────────────, - ────────────────────────── + ──\n",
       "      2                                3                                2     \n",
       "₀    ⎞        ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞        ⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞      \n",
       "  + 1⎠        ⎝ℯ                   + 1⎠        ⎝ℯ                   + 1⎠      \n",
       "                                                                              \n",
       "             -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁⋅u₁₀              -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    2      \n",
       "W₁₀⋅W₁₁   2⋅ℯ                      ⋅W₁₀⋅W₁₁     ℯ                  ⋅W₁₁     2⋅\n",
       "─────── + ─────────────────────────────────, - ────────────────────────── + ──\n",
       "      2                                3                                2     \n",
       "₀    ⎞        ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞        ⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞      \n",
       "  + 1⎠        ⎝ℯ                   + 1⎠        ⎝ℯ                   + 1⎠      \n",
       "                                                                              \n",
       "             -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀              -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    2      \n",
       "W₂₀⋅W₂₁   2⋅ℯ                      ⋅W₂₀⋅W₂₁     ℯ                  ⋅W₂₁     2⋅\n",
       "─────── + ─────────────────────────────────, - ────────────────────────── + ──\n",
       "      2                                3                                2     \n",
       "₀    ⎞        ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞        ⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞      \n",
       "  + 1⎠        ⎝ℯ                   + 1⎠        ⎝ℯ                   + 1⎠      \n",
       "\n",
       " -2⋅W₀₀⋅u₀₀ - 2⋅W₀₁⋅u₁₀    2⎤⎤\n",
       "ℯ                      ⋅W₀₁ ⎥⎥\n",
       "────────────────────────────⎥⎥\n",
       "                         3  ⎥⎥\n",
       "⎛ -W₀₀⋅u₀₀ - W₀₁⋅u₁₀    ⎞   ⎥⎥\n",
       "⎝ℯ                   + 1⎠   ⎦⎥\n",
       "                             ⎥\n",
       " -2⋅W₁₀⋅u₀₀ - 2⋅W₁₁⋅u₁₀    2⎤⎥\n",
       "ℯ                      ⋅W₁₁ ⎥⎥\n",
       "────────────────────────────⎥⎥\n",
       "                         3  ⎥⎥\n",
       "⎛ -W₁₀⋅u₀₀ - W₁₁⋅u₁₀    ⎞   ⎥⎥\n",
       "⎝ℯ                   + 1⎠   ⎦⎥\n",
       "                             ⎥\n",
       " -2⋅W₂₀⋅u₀₀ - 2⋅W₂₁⋅u₁₀    2⎤⎥\n",
       "ℯ                      ⋅W₂₁ ⎥⎥\n",
       "────────────────────────────⎥⎥\n",
       "                         3  ⎥⎥\n",
       "⎛ -W₂₀⋅u₀₀ - W₂₁⋅u₁₀    ⎞   ⎥⎥\n",
       "⎝ℯ                   + 1⎠   ⎦⎦"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = sy.Matrix([[[J[i ,alpha].diff(u[beta,0]) for alpha in range(dim_u)] for beta in range(dim_u)] for i in range(dim_x)])\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.03468516,  0.00569407],\n",
       "        [ 0.00569407,  0.00093476]],\n",
       "\n",
       "       [[ 0.10130854, -0.04189015],\n",
       "        [-0.04189015,  0.01732119]],\n",
       "\n",
       "       [[-0.34631993,  0.25240458],\n",
       "        [ 0.25240458, -0.18395728]]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sH= np.r_[[[[sy_eval(H[i,alpha][beta], W, u, nW=tW.numpy(), nu=tu.numpy())[0] for i in range(dim_x)] for alpha in range(dim_u)] for beta in range(dim_u)]].astype(float).T\n",
    "sH"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now in tensorflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "@tf.function\n",
    "def get_hessian(tW, tu, transformation):\n",
    "    with tf.GradientTape(persistent=True) as t1:\n",
    "        with tf.GradientTape(persistent=True) as t2:\n",
    "            tx, ts, ta = transformation(tW, tu)\n",
    "        _tJ = t2.jacobian(tx, tu)[:,0,:,0]\n",
    "    tH = t1.jacobian(_tJ, tu)[:,:,:,0]\n",
    "    return tH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(3, 2, 2), dtype=float64, numpy=\n",
       "array([[[ 0.03468516,  0.00569407],\n",
       "        [ 0.00569407,  0.00093476]],\n",
       "\n",
       "       [[ 0.10130854, -0.04189015],\n",
       "        [-0.04189015,  0.01732119]],\n",
       "\n",
       "       [[-0.34631993,  0.25240458],\n",
       "        [ 0.25240458, -0.18395728]]])>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tH = get_hessian(tW, tu, transformation)\n",
    "tH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose (tH.numpy(), sH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### the christoffel symbols\n",
    "\n",
    "\n",
    "$$\\Gamma^\\beta_{\\mu\\nu} = \\frac{1}{2}g^{\\beta\\alpha}\\Bigg(\\frac{\\partial g_{\\alpha\\mu}}{\\partial u^\\nu} + \\frac{\\partial g_{\\alpha\\nu}}{\\partial u^\\mu} - \\frac{\\partial g_{\\mu\\nu}}{\\partial u^\\alpha}  \\Bigg)$$\n",
    "\n",
    "\n",
    "\n",
    "$$\\begin{align}\n",
    "\\text{The Metric Tensor}\\;\\;\\;: \\;\\;\\;&g_{\\alpha\\beta} &= \\;\\;\\;& \\sum_{k=0}^2 \\dot{\\text{a}}(s^k)^2 W^k_\\alpha W^k_\\beta\n",
    "\\end{align}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{\\partial g_{\\alpha\\beta}}{\\partial u^\\gamma} =  \\sum_{k=0}^2 2\\dot{\\text{a}}(s^k)\\ddot{\\text{a}}(s^k) W^k_\\alpha W^k_\\beta W^k_\\gamma$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(3, 1), dtype=float64, numpy=\n",
       "array([[-0.5759448 ],\n",
       "       [-0.94727428],\n",
       "       [ 1.92824448]])>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# del(get_christoffel_symbols_v1)\n",
    "# del(get_christoffel_symbols_v2)\n",
    "#            tx, ts, ta = transformation(tW, tu)\n",
    "\n",
    "i,j,k = tf.Variable(-1), tf.Variable(-1), tf.Variable(-1)\n",
    "\n",
    "@tf.function\n",
    "def get_christoffel_symbol(tg_, tdg_du, i, j ,k):\n",
    "    return 0.5*tf.reduce_sum([tg_[i,a]*(tdg_du[k,a,j] + tdg_du[j,a,k] - tdg_du[k,j,a]) for a in range(tg_.shape[0])])    \n",
    "\n",
    "@tf.function\n",
    "def get_christoffel_symbols(tW, tu, transformation):   \n",
    "    # using TF's full symbolic engine \n",
    "    with tf.GradientTape(persistent=True) as gtape2:\n",
    "        with tf.GradientTape(persistent=True) as gtape1:\n",
    "            tx, _, _ = transformation(tW, tu)\n",
    "\n",
    "        tJ = gtape1.jacobian(tx, tu)[:,0,:,0]\n",
    "\n",
    "        tg = tf.matmul(tf.transpose(tJ), tJ)\n",
    "        tg_ = tf.linalg.inv(tg)   \n",
    "        \n",
    "    # let TF compute the derivatives of the metric tensor wrt source variables\n",
    "    tdg_du = gtape2.jacobian(tg, tu)[:,:,:,0]\n",
    "\n",
    "    r = []\n",
    "    for vi,vj,vk in itertools.product(*[list(range(u.shape[0]))]*3):\n",
    "        i.assign(vi)\n",
    "        j.assign(vj)\n",
    "        k.assign(vk)\n",
    "        r.append(get_christoffel_symbol(tg_, tdg_du, i, j, k))\n",
    "    \n",
    "    return tf.reshape(r, [tu.shape[0]]*3)\n",
    "\n",
    "@tf.function\n",
    "def get_christoffel_symbols_PERCEPTRON(tW, tu, transformation):   \n",
    "    # using the perceptron expressions for the derivative of the metric tensor\n",
    "    # which just require the first and second derivatives of the activation function\n",
    "    with tf.GradientTape(persistent=True) as gtape2:\n",
    "\n",
    "        with tf.GradientTape(persistent=True) as gtape1:\n",
    "            tx, ts, ta = transformation(tW, tu)\n",
    "\n",
    "        tJ = gtape1.jacobian(tx, tu)[:,0,:,0]\n",
    "        tg = tf.matmul(tf.transpose(tJ), tJ)\n",
    "        tg_ = tf.linalg.inv(tg)   \n",
    "        tdx_ds = gtape1.gradient(tx, ts)\n",
    "        \n",
    "    \n",
    "    tdx_dsds = gtape2.jacobian(tdx_ds, ts)[:,0,:,0]\n",
    "    tdx_dsds = tf.linalg.diag_part(tdx_dsds)\n",
    "\n",
    "    # compute manually the derivatives of the metric tensor wrt source variables\n",
    "    tdg_du = tf.reshape(\n",
    "        tf.stack([tf.reduce_sum((tW[:,alpha:alpha+1] * tW[:,beta:beta+1] * tW[:,gamma:gamma+1])[:,0]*2*tdx_ds[:,0]*tdx_dsds) \\\n",
    "                   for alpha, beta, gamma in itertools.product(*[list(range(dim_u))]*3)]),\n",
    "        (dim_u, dim_u, dim_u))    \n",
    "\n",
    "    r = []\n",
    "    for vi,vj,vk in itertools.product(*[list(range(u.shape[0]))]*3):\n",
    "        i.assign(vi)\n",
    "        j.assign(vj)\n",
    "        k.assign(vk)\n",
    "        r.append(get_christoffel_symbol(tg_, tdg_du, i, j, k))\n",
    "    \n",
    "    return tf.reshape(r, [tu.shape[0]]*3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 0.23398133  0.0096741 ]\n",
      "  [ 0.0096741  -0.00869439]]\n",
      "\n",
      " [[-1.47620478  1.28079288]\n",
      "  [ 1.28079288 -0.89929744]]]\n",
      "867 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -n 1 -r 1\n",
    "tcf2 = get_christoffel_symbols(tW, tu, transformation)\n",
    "print(tcf2.numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 0.23398133  0.0096741 ]\n",
      "  [ 0.0096741  -0.00869439]]\n",
      "\n",
      " [[-1.47620478  1.28079288]\n",
      "  [ 1.28079288 -0.89929744]]]\n",
      "824 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -n 1 -r 1\n",
    "tcf1 = get_christoffel_symbols_PERCEPTRON(tW, tu, transformation)\n",
    "print(tcf1.numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(get_christoffel_symbols_PERCEPTRON(tW, tu, transformation), get_christoffel_symbols(tW, tu, transformation))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transformation(tW, tu):\n",
    "    tlinear = lambda x: x\n",
    "    ta = lambda x: 1/(1+tf.exp(-x))\n",
    "    ts = tf.matmul(tW,tu) + tb\n",
    "    tx = ta(ts)\n",
    "    return tx, ts, ta\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Variable 'Variable:0' shape=(3, 2) dtype=float64, numpy=\n",
       " array([[-1.15159902, -2.99051546],\n",
       "        [-0.38944001,  0.34125733],\n",
       "        [-1.39846121, -0.18306849]])>,\n",
       " <tf.Variable 'Variable:0' shape=(2, 1) dtype=float64, numpy=\n",
       " array([[-0.61250047],\n",
       "        [ 0.75652911]])>)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tW = tf.Variable(np.random.normal(size=(dim_x,dim_u)))\n",
    "tu = tf.Variable(np.random.normal(size=(dim_u,1)))\n",
    "tb  = tf.Variable(np.zeros((dim_x,1)))\n",
    "tW, tu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### check with `sympy` on a  2D manifold embedded in 3D ambient space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transformation_sphere(tW, tu):\n",
    "    tx = tf.stack([tf.sin(tu[0])*tf.cos(tu[1]), \n",
    "                   tf.sin(tu[0])*tf.sin(tu[1]), \n",
    "                   tf.cos(tu[0])])\n",
    "    return tx, None, None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "tu = tf.Variable(np.array([[-0.64853741], [ 0.45577495]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "r, t, p = sy.symbols('r theta phi')\n",
    "\n",
    "x, y, z = sy.symbols('x y z')\n",
    "# define the metric \n",
    "metric = [[0 for i in range(3)] for i in range(3)]\n",
    "metric[0][0] = 1\n",
    "metric[1][1] = 1\n",
    "metric[2][2] = 1\n",
    "# creating metric object of 3D ambient space\n",
    "mh = MetricTensor(metric, (x, y, z))\n",
    "\n",
    "# ------------------\n",
    "# the transformation\n",
    "# ------------------\n",
    "xf = sy.sin(t)*sy.cos(p)\n",
    "yf = sy.sin(t)*sy.sin(p)\n",
    "zf = sy.cos(t)\n",
    "\n",
    "esyms = [t, p]\n",
    "hfuncs = [xf, yf, zf]\n",
    "\n",
    "me = [[0,0],[0,0]]\n",
    "for ki,kj in itertools.product(*([range(len(esyms))]*len(esyms))):\n",
    "    me[ki][kj] = (sum([hfuncs[m].diff(esyms[ki])*hfuncs[n].diff(esyms[kj])*mh[m,n] \\\n",
    "                for m,n in itertools.product(*([range(len(hfuncs))]*len(esyms)))])).simplify()\n",
    "me = MetricTensor(me, (t,p))\n",
    "cf = ChristoffelSymbols.from_metric(me)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.        ],\n",
       "       [0.        , 0.36484187]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sme = me.tensor().subs({t: tu.numpy()[0,0], p: tu.numpy()[1,0]})\n",
    "sme = np.r_[[sme[i,j] for i,j in itertools.product(range(dim_u), range(dim_u))]].reshape(dim_u, dim_u).astype(float)\n",
    "sme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.        ],\n",
       "       [0.        , 0.36484187]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tme = get_metric_tensor(tW, tu, transformation_sphere).numpy()\n",
    "tme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.        ,  0.        ],\n",
       "        [ 0.        ,  0.48138579]],\n",
       "\n",
       "       [[ 0.        , -1.31943681],\n",
       "        [-1.31943681,  0.        ]]])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scf = cf.tensor().subs({t: tu.numpy()[0,0], p: tu.numpy()[1,0]})\n",
    "scf = np.r_[[scf[i,j,k] for i,j,k in itertools.product(range(dim_u), range(dim_u), range(dim_u))]].reshape(dim_u, dim_u, dim_u).astype(float)\n",
    "scf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[-5.55111512e-17,  5.55111512e-17],\n",
       "        [ 5.55111512e-17,  4.81385791e-01]],\n",
       "\n",
       "       [[-7.60756315e-17, -1.31943681e+00],\n",
       "        [-1.31943681e+00,  0.00000000e+00]]])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mcf = get_christoffel_symbols(tW, tu, transformation_sphere).numpy()\n",
    "mcf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "check on many random points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (200 of 200) |######################| Elapsed Time: 0:00:38 Time:  0:00:38\n"
     ]
    }
   ],
   "source": [
    "from progressbar import progressbar as pbar\n",
    "from time import time\n",
    "ttimes = []\n",
    "stimes = []\n",
    "for _ in pbar(range(200)):\n",
    "    tu = tf.Variable(np.random.normal(size=(dim_u,1)))\n",
    "    s_start = time()\n",
    "    sme = me.tensor().subs({t: tu.numpy()[0,0], p: tu.numpy()[1,0]})\n",
    "    sme = np.r_[[sme[i,j] for i,j in itertools.product(range(dim_u), range(dim_u))]].reshape(dim_u, dim_u).astype(float)\n",
    "    scf = cf.tensor().subs({t: tu.numpy()[0,0], p: tu.numpy()[1,0]})\n",
    "    scf = np.r_[[scf[i,j,k] for i,j,k in itertools.product(range(dim_u), range(dim_u), range(dim_u))]].reshape(dim_u, dim_u, dim_u).astype(float)\n",
    "    s_elapsed = time()-s_start\n",
    "    \n",
    "    t_start = time()\n",
    "    tme = get_metric_tensor(tW, tu, transformation_sphere).numpy()\n",
    "    tcf = get_christoffel_symbols(tW, tu, transformation_sphere).numpy()\n",
    "    t_elapsed = time() - t_start\n",
    "\n",
    "    stimes.append(s_elapsed)\n",
    "    ttimes.append(t_elapsed)\n",
    "    \n",
    "    assert np.allclose(sme, tme)\n",
    "    assert np.allclose(scf, tcf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3hc1bW33z1do15suciWXHFvcgEM2IBpcagBAgECJFxSILkkX0j5wgc3hYQbUiC55FKCISEhkBBCCMUNW7jg3rtly0WSLVu9a+r+/jhzjmakkTQjjQqj/T6PH2vmtDXn7LPOOr+99tpCSolCoVAo4hdTfxugUCgUit5FOXqFQqGIc5SjVygUijhHOXqFQqGIc5SjVygUijjH0t8GtCUrK0vm5eV1e/vGxkYSExNjZ1CMUHZFx0C1Cwaubcqu6BiodkH3bNuxY0eFlHJI2IVSygH1Lz8/X/aEtWvX9mj73kLZFR0D1S4pB65tyq7oGKh2Sdk924DtsgO/qqQbhUKhiHOUo1coFIo4Rzl6hUKhiHOUo1coFIo4Rzl6hUKhiHOUo1coFIo4Rzl6hUKhiHOUox8gHCmrZ+uJqv42Q9EHtHh8/W1Cv1Ne72LHKdXe+wrl6AcIv/2okMfe2dffZih6md3FNUz/rxWcrW3ub1P6lVc2nuC+V7b1txmDBuXoBwgtHh9NbhXpxTsl1U14fJLzda7+NqVfaXR5VXvvQ5SjHyC4fX5aPP7+NkPRy3h8/pD/Bysev8QX+KfofZSjHyB4fH5cXhXhxDser+bY3IPd0XvVA68vUY5+gOD1SVwqoo97dAfv9Q3uSFa92fQtytEPEDw+P26fH796lY1rlIPT8AQedG7v4D4PfYVy9AMEt0+90g8GlKPXaD0PgzuweXZ1IS+tK+r140Tk6IUQ1wohjgghjgkhvh9m+beFEAeFEHuFEB8JIXKDlvmEELsD/96NpfHxhN7wB1OO9ZrD57jn5S1opbQHB0YkO8gdnHrgaaw+dI61R873+nG6nGFKCGEGngOuAkqAbUKId6WUB4NW2wXMlVI2CSG+BvwC+HxgWbOUclaM7Y479AbvGkSvsrtO17C+sAKX14/Dau5vc/oEw8ENouscDo96gwW09tAXD7tIIvr5wDEpZZGU0g28AdwYvIKUcq2UsinwcTOQE1sz4x/9xh9MHbL6TT6YbnYVyWq41XkAtPPQF293kcwZOxIoDvpcAizoZP0vAx8GfXYIIbYDXuApKeU7bTcQQjwIPAiQnZ1NQUFBBGaFp6GhoUfb9xZd2dXYog2gWf/JZkYm913XSX+er6KT2m8u+HgDKXYRsmygXkfomW3HT7gBOHj4CAVNsdVmB+o5C2dXVbU2Mnjzlm2UpfbP29xAOF/1DU14WkQ7O2JtW0wnBxdC3A3MBRYFfZ0rpSwVQowF1ggh9kkpjwdvJ6V8EXgRYO7cuXLx4sXdtqGgoICebN9bdGWXWLcS8DBjdj7Tc1IHjF29yarqfXDqNHMXXMiItIQBY1dX9MS2groDcPIkeWPHs/iSMQPGrt4knF2/2rcBamqZPms2+bkZA8auvsb0yWrsDiuLFy8K+T7WtkUSOpYCo4I+5wS+C0EIsQT4IXCDlNIY3y2lLA38XwQUALN7YG/cYkg3g2jQlJ5aN5hS7JR0o6H/frd3sHdKywGj0W8DJgghxgghbMAdQEj2jBBiNvACmpM/H/R9uhDCHvg7C1gIBHfiKgLonVODqQyCR2n0gxal0Wt4vP4+6Zjv0tFLKb3Aw8AK4BDwNynlASHEj4UQNwRWexpIAv7eJo1yMrBdCLEHWIum0ceFo99bUsOyDSdisi8ppdHwuxPRbyisYN6Tq6lpcsfEnr7C6IwdVBG9Sq+E1pHBg93RuwMDJXubiDR6KeUHwAdtvns86O8lHWz3CTC9JwYOVN7YVsxbO0r4Ugx0Vm/QaNjuRPTv7T1Deb2Lw2X1XDg2s93ylzecoMXj46HLx/fIzljjNuSqwXOzq0hWQ73ZtAZ4fRHoqJGx3aSywYXb649J9b3guifRRvRSStYXVgBwsqIx7Dof7jvLe3vPdt/AXsI1GDV6vZjXIPrN4WiV7Qbvm43PL5Gyb0YHK0ffTSoaNJmkye3t8b6CX92ijW5PVzVRWqOlqp2oDO/oG1xe6ls83TewlxjMGr13kNc0cqsHnuHgB0pnrCIMFQ1aYlFzDCZPCL7Q0ZZA2HisEoAku6XDiL6+xUt9S88fSLFmcGbdqBGhoM4DBFUy9cteL2YY0zz6wURFveboG2Ps6KON6Dceq2BYioNpI1M4WdEUdp36Fg8NLi9SSoQQYdfpDwZjZ6yh0Q+i3xwOr19p9MHt3u3z4zD13sAxFdF3g2a3z3DwsZBuPEG5xNGWQNh5upoFYzPIy0zkZGVju8hASkmDy4tfxuahFEuMiN43sOzqTVQnpNYmVZni0DbQ2+1BOfpuoMs2EBvpJvj1tSXKzti6Zg9ZSXbyshJxef2U1bWELG9y+9B9/0DT6QendKPK8wb/9sF9HoIi+l6+B5Sj7wbBjj4WExzrr7EQXUQvpaTJ4yPBamZMViLQPvOmwdX6xjHQdPrBGNV5lTbdp5HsQCb0PPTuA085+m5Q2dA6MCkWjj5Yuokmond5/UgJCbZWR9828ybYudc1D6yI3jWI8+i9g9jBeUMi+sF7HoLbvZJuBiAh0o0nxumVUUT0+kPGaTMzLMWB3WJqF9EHyzUDLaJ3Bx5qg8nRK+kmtL0Ppre5tgS3gd5+w1OOvhsEO/pGV6yzbiLfn94R7LSZMZkEozOcnK4KzbwJlm7qBppGPwizbvS3t/6Wbk5WNLK3pKZfjh2iTQ/iiF5p9AOcigY3FpOWphj7PPrIL7iec59g07Jkh6U6KKtzhazTECzdDLCIfjDmUg+UrJtfrjzCI2/s7pdjx1KjL693sfCpNRwpq++pWX2OR0k3A5uKBpdRPz0mnbEBh2cSWkS/4kAZv1p5pMvtDOkmMA1fdoqDc7WhWTfBcs1Ayrrx+aVRPmIwRfQDpdZNXYuXM7XN/TJfb0jWTQ/LFJ+oaKS0ppnC858+R+9S6ZUDm4oGF9kpdmwWE00x1OiT7BZcXj//3nOGP28+1eV2wRo9wLAUB+UNrpD6O/V9kHVT0eDi/le2Utng6nrlACGDRQaRo2+dM7Z/NfoWt48Wj5+65r5/y4tlRK8HL5/G8t7BEX1v91MpR98NKhvcZCbacdrMMZVukh1WXB4fVY1u6lu8XUZb+rETbHpEb8fnlyEOV5duUhOsvRbR7yutZe2RcvZEofkOXkcfqG/i79/f3ByQ/c7Vt3SxZuyJpUav90E1R1k6ZCDQl+MJlKPvBhUNLrKSbTit5tikV7aJ6Ksa3Xj9ssunfFM7R+8ACBk0Vd/iwWkzk+a09lr01hTokA5OO+0K9yDskAuWq/pbutE78stq+9nR9/Ahr/c7uT6Vjj7ozUZF9AMLj89PdZM2GtVpt8S0BEKyw0KLx0dlo+Ywu8qSMbJurK2dsQDngjpkG1xekuwWUhy9F9HrdlRHMfHJYEyxC72x+1m6CUgd5+r6w9HHLo++VbqJztHvPF3Nr7e3RHX8kuqmmPZpuFVn7MClpklrWBmJNpy2ziP6SBuFodE7LLR4/FQHHH1DF5p6syeCiN7lJdlhIdlh6TWNXj8H+gMqEtx9qE8OFAbSiFD94Xy+PvJ+lVgRyxGhepuOVqPfdqKKvRU+yiP8/ScqGrnkv9ey4VhF1DZ2RF++1SpHHyX6DZJkt5DQiXRz03Mb+e1HxyLapzdIuqlsdBm1yrtyzG07Y7OS7JhNIiTzpr7FS5LD2quOvjFwTqqikW68g0+68QRlV/X3b9aDhP6UbhKs5p5r9C3d0+j1eyfSe6I4MD5lf2ldVMfpDJVHP4DRB0g5bZYOO2P9fsn+0loKjp5vtywcugNIdlhDIpzgwU7hMDpjA+mVZpNgSJI9JKJvaPGQHJBuemvAlK7RVze5kVKy7mh5lzNvhTbyT5++2h3035xos/RrRO/3y36VbtwB2SrRbu436UZ/MEQqZ+qy5ImKhqiO0xmq1s0ARo/oE+1mnLbwGn1diwevX3LwTF1EDVmPapIdodMDdNUImz0+HFYTJlNrjfnsFHvIzdtgSDfWXo/oKxvd7Cut5YvLtrLuaHmn27gGYdaN/juddnO/lkAIrqfUPxp94DzE4IHXXelGv28jvSf0RIMTHUzu0x2URj+AaXR3HdHrJRJcXj9Hz3U9kCM40guma+nGi7PNNtkpjpCbt75F64xNdlhocHljMsdtOzsCEX1Vo9u4Eaq60OsNp2fr+ev7pwVdknPaLPj6YFahjghus+fq+l6j16u1Om3mHndKG44+yrdCXbqJ9C23NaKPoaMPrnWjpJuBRaMrOKI3h53MoyJIq95XUtvlPj0+P2aTMLR2nUg0el220RmW6gjRXRtavCQFOmOhazmoOzR5Wh29Pn9tV8fRnXui3TJoIvrWSFa7Zv2VS687uewUe7sBdn2Bx5BuLD1+yOuOuiXKNOfmKDV6PXCpaHBTG6MqsH1Z80c5+igxHL3NQoLN0mlED7AnIkcvsZoFdmvo5YhEo09o83DITnFQ1+Kl2e3D75c0uL0kO6ykJFiB3ilV3ORqfQ3Wq2d2ZbueN5w8iBy9/jv1N7eeyjf3vbKVt3eWRL2drmfnZia2G2DXF7iDHng9vfZ6O+vtiD74DbWjuZmjxeP1G4Gakm4GGMGZLrrs0La2uK7nTRqWzL7SrkeLenx+rGYTDot20R1WEwlWc5cafZPb1+4tQE+xPFfXQqPbi5QEOmM159IbOn1jUD/FvkBWQlfHCU4pHSyO3ojo7YGbuwe/2+PzU3CknLd3lka9rd6Gx2Rqcxj0tXwT/GbTXxp9dyL6zEQbAEUx6pB1+/w4rCYsJqGkm4FGo9EZazGcbFObHv+KBhcmAYsmDuFIWX2XGQEenx+b2WRE9BlOG0kRpEM2h5Fu9MZY3eQ2oh29MxZ6p7BZcIqp3ifR4Or8OHrDTorB6/unBT2CN6SbHvxu/c1sx6nqqPejZ5zkBSaraTv9ZG+jF/HrafaRPh8yRJ91o9eoans/1DZ7+N1Hhe32V93kZuaoNEwCTpTHKKIPBHhWs0lF9AONJpcPs0lgt5iMjtC28k1Fg5uMRBsXDEvG45OcCejWHeHxSixmgT0Q0Wck2bS89y7kjyaPt11EnxSkxesPiiSHhTSn5uijGdQUKY0uL1lJdgBD7+1qsFewox9sA6b0dtOTB1xNwNE3e3zsK+1aHgxGd/Rjh2iO/mxt5+0z1gT3z/REvmpy+4z2Fqs8+n/tLuVXq47yxtbTId9XNbrJTnEwKsNJUYykG7dXBhy9UOmVA41Gt+ZchWjtPG07aKqiwUVWkp0UI4ruQq8OPNn1iD7daSPZ3nVEr0k3oVk3SfaAo28JcvR2C+OGJGE2CQ6djd2Aj2A7ctITQr5r6GJCluCKnYNFunG36Yz19uDm1kdoA2w9URXVtnpgMirdic1iMgYD9RXB0k1PHnbB90c0M7NBx9LNliLtXL60/oRhp98vqW7ykJFoJS8zMWaZNx6fH7vFhM3S+5lnytFHSZPLZ3Sm6R2hjW0i78oGF5mBqBwi06ttZhN2i3Y5MhNtJDusNHQhs7SE6YzVHX29y2u8liY7rDisZiYMTWJ/F9FfWW1L1I63ye1jVIYz5LuIpRuHJt30R130vsbjDY3oYyHdmE2CLUWVUW2rO7lEu5mc9ARKqvs2ovf4/FgCb8Vub/evvd7GUhOsHUo39S2esJOStEb0re1USsnmokpGpiVQWtPMe3vPBNbR0pIzEu1kp9hDki16gturBXg2s9LoBxwNbq/RmaZHZm1fGysa3GQl2SPWxY3O2IDenpFoJymSiN7TvjPWSKMMiuhTE7Tvpo5IZf+Z1oj+o0Pn+NexVimnpsnN5b8s4E+bTnZ63HZ2uL2MSHMgAuO2hqU4Ik6vTLJbkLI1xzye8RjatHbNeibdaNdtfl4G209WR5Ui2RRUI2lUupPi6r6N6L0+Taq0mjX3091rr1euHJJs79DR//jfB7n9hU3tvg8X0R8730Blo5uHrxjPhKFJvLTuBFJKKhs1x56RaCXdaaO6yROTwMTj82O1CKyWAaLRCyGuFUIcEUIcE0J8P8zybwshDgoh9gohPhJC5AYtu1cIURj4d28sje8PmlxeI6LvSLqpbHCRmWiPOKL3+iRWizAi+oxEqzHAqVNbwkT0ifZWjV5PHdMlpGkjUyivd3E+0Pn2xrZi3ivyGAN31h45T7PH12XUH4zb68fjkyTbLaQ7tY7gC4YlR67RB87RYJBvWrNuep5eqUs3l00cQr3L22U/UDB6znmCtX8ierfeCRlo7911cvp9NSTJTrPH1875Nrq8vL/vLLXNnpD25fX5jYds8L25OSCBXTwuk3svzuPg2Tp2nq4xBkulO22kOq24vf6YTHTiHkidsUIIM/AccB0wBbhTCDGlzWq7gLlSyhnAW8AvAttmAE8AC4D5wBNCiPTYmd/3NAalNLZ2xrY2lma3j0a3j6xkW2vuehcRvTtcRN9F1o3PL3F7/UaJYh1tPybN0TfrWTe6o08FYP8ZzZGfrmzC44ezAce/+qBWm+d4FFkFrROUW0h3Wkm2WxiRltB1RO/1YxKtdXoGg6M3OiFjkHWjO/pJw5IBOBtFcbKmIEc/KsNJTZOnT6eZ1LPM9Ii+u6NjdZuHJNvxy/YPzuX7y4zfGiyvBmfJBd+bm4sqGZbiYHSGk5tnjyTZbuG1TSeNdOnMRLsRzERTkrsj9PNgM5sGhHQzHzgmpSySUrqBN4Abg1eQUq6VUurvf5uBnMDf1wCrpJRVUspqYBVwbWxM7x+a3F5DBw8X0ev6XVZAfoH287a2rY/j8fmxmkzkZSbylcvGsmTKUE2jd3k7HCav7yPB1v4SJtm1ujZ1LR6sZoEj0Mk7eXgKQmgV+KSUnA50wp2qaMTl9fFxoD5NUXlDxK+mjUF6b2aSnZHpCaREkBqqy1W2QFQ3GFIs20X0Pbi5a5s9pDgsjMrQOsGjyZxp9viwmU1YzCZGpWt9K8VVfRfVewLZJjazpvV199rXB0k30H7Q1Nu7StqtC62yTYqNkHtsT3EN+XnpCCFItFv4XH4OH+wro/C8ljefnmglPZC9FgtH7/b6sVm0Nxt3L2fdWLpehZFAcdDnErQIvSO+DHzYybYj224ghHgQeBAgOzubgoKCCMwKT0NDQ4+274qK6iYSfY0UFBRQ49Ia6O79h8io00oSH6/RGtGZE0dY33gchxkOHjvBhFEeCgoK+NmWZlLtgodmOVr3WdWMzQzr133MRU44uOMc50q0SGP5RwU4rYK21LRoxy45WUSBvzhkmdnvpuh0KYlWgcMs+fjjj41l2U5BwZ7jjPYWG30LKz7Zxc7dggaXlxlZZvZW+HhnxVrSHV3HAaUNmh0njx3liiECrx+Kzhbj8vpZvWYtFlN72wGKTrkw4aeo8CgAH2/4hKHO1uP19nXsCd217eAp7ZoeP3IQgO27duMuieQWbM+REy3YhZ9je7cDsGHnARYNdUdk17ETLqwmPwUFBZyr1drA8vVbOZ/dPVu6ou35KjnjwuvxUXS8EIB1GzaSmRB9d+HuE9r5rD+vOfQ1H68nza7tp8kj+eRYEzlJgpIGydqNm8hN0QKzskatzaZYJXVuwYo1BVhNUFrdzOwMr2HrRJMm8bywVmujB3Zs4WSdtu3Hm7ZTnhkqm0ZLVU0zPrug2StxNxJyjmLd/mN6ZYUQdwNzgUXRbCelfBF4EWDu3Lly8eLF3bahoKCAnmzfFfKT1YwZNZTFi2do8sTaFRx3J7NwSB7XTB2G9+A52Lydyy+ay8xRaaRv+ojUrCySkqqZPX8hhStWkpvhDLHx1/s3kJFoY/Hi+cZ3Zc7TvHlkH7PmXciItIR2dpysaISCAmZOnczi/JyQZUP3bcCZZCPJYSWruTbkWAvKdrGlqIpRk2bBWq2TypY5kgqfxGE9zX8unc2X/7id7PEzuHh8VpfnY09xDWzYyLzZ07liUjYAr2w8wduFB8lfsJD0wACutqys3kdiZRkzpk2B/buZM3ce44cmG8t7+zr2hO7aVriuCA4dYkH+LNi5hSlTp7N4Sna3bHj1xFaGm91cu+QSUjasICFzBElJFRHZ9UHFHpJrtHVnNLr58aZVpI4cx+JLxnTLlq5oe77+cXYXKe5apk8ZD/v3kD9vgTF4Kxp2uo8gjh5j7ozJ/O3oXubMvZDRmdobyrHzDciPPmbhpBze3F7MBdNmceHYTAAOnKmF9RvIcFooafQxc+6FeHx+5MoCLp09mcVzRxnHeL9sMxuPVeKwmrhmyeUcLqvjqa3ryZ0whcUzhvfovDh2r2NYppMGlxeXx8/ixRcby2Ld/iN5jJYCo4I+5wS+C0EIsQT4IXCDlNIVzbafJppcrbnriTYzl18whF2na/jKaztYdfAcJYEMhqzA66Q24YcWeWw+UYmUUFzd3K4Wta5X6iR10ZHbdtKRkG3tFqMztm3p41mj0iira2FLoOPJZoYTFU1sLqpkXl4GU0akAHA8wlzhxiCNPvj40Hm9G7c3NKV0MAya0ouYJcYgvbKmyUNqoA9oRFoCZ2oi1+ibPX6jEz/dacVpMxvtti/weFs7IaH756GuxUuSrXWEerB0o9fvyc3SHH9DGOkmza69bda3eA3pSpeydO65MA/QRqsDMdXo3T4/NosZq9nU69JlJBH9NmCCEGIMmpO+A/hC8ApCiNnAC8C1Usrg2TZWAD8L6oC9GvhBj63uJ6SUNLq9JAbSK4UQvHL/fNxePzc9t5HvvrWHZo+PCUOTyDYcfWsd+E3HtXxnn19SWt1sRDF6p0wwegdqR/nozR5dow/j6B0Wiqua8PmlkXGjM3u0dine3X0GIWBSupm9JTWcr3dx/cwRDEtx4LSZOX4+snoeTa72D5xIso201LIgjX4wOHpvaAmEntzcdc0eY5Da8FSHlnWT28VGAZrdXqMTXAihpVj2pUavpxWae9Y/U9+izbWg14gKTrHUR4DnBer51AfdR3qQlG44eo/RX6W/EegsmTyUEakOMgMjv/WHa02MNHprIM203ztjpZRe4GE0p30I+JuU8oAQ4sdCiBsCqz0NJAF/F0LsFkK8G9i2CvgJ2sNiG/DjwHefSlo8fvySdqNRbRYTz9wxiya3jzFZSbz+HxdiCTTi4Cn8PjleYZQiOFnZGjFrHZOhWrYeFdd14Cyb3aGDb4JJDkT09S1eUhJCl08ZnoLNYuJwWT3DUxyMSDIZ84ZeODYDIQRjhyRGPMw7fERvDVkWDj2it5kHT9aNXo7aYVQs7EF6ZbPHaEvD0xI4W9tMZbOfH//7YLsie21p9oSm5Y7KSOjbiN4vsZhM2Cxam+/ueahv8RiDASG0sJnu6HMz20f0TWEi+tNVTVjNgmEprX1nABazif+5aw7/77NaoqHDaibBag4ZmdxdjKybPsijj0ijl1J+AHzQ5rvHg/5e0sm2y4Bl3TVwINEYNLtUWyZmJ7PmO4vJTLQZDQ+0HPaTFY3UuuDouUa+fMkYXt5wQtPYL9DW8Xj9xoOhdbvWgU/haE1rDB/RN7i8eHx+ku2hEb3NYmLaiBR2nq5hVIaTYU4tcndYTUwfmQbAuCFJbD9Z3eX50OxozboJPn5ntkNrxoFe9mGwZN1YgwYKdffm9vslNU1u0hI0GWFEqoPqJg/LT0pWnTrBrfk5hgQXjia3zwgkQHv4rzl8nrO1zQxPbd8fFGs83jbpld08D1WNbtKcViPzLHjgoiHdGBF9kHQTeBtOc2iOvq7FQ3FVEznpTsxhkgfmjA7NCE93WqmOiaPXs49MqtbNQKIpaL7YcIxMSwhx8tAa0R8LZON8ZvpwEm1mTla2RlDubmj0zUGjG9uSZLfQ0OLVUvAS2tuqyze5mU6yA5kuc3MzDBll/JAkztQ2RzRJSWtfQXuNXr+5vD5/u5GLmj5pMiSr3oroy+tdXUa4fYU7kEZrMeuRbPfsanB78UuMiF7vrP/kjHa+uxo81ez2hbTTW/NH4Zfw5rbiTraKHW2lm+6mmZ6qaiI302kUAwxuY1WNblITrCTazFjNIuQ+ai/daBF92zIeHZHmtMVEuvHo6ZWqBMLAQo/ok8JE9B2ha/TlTdoTe/yQJHIzEznVRrqxtZFuMhPtmASUdZAfHTzopS2JdgvewATQbTV6gNmjtch9dIaTYYnacS8al2ksnzYyFSnhQAQjZPVJR8Jp9HpE/8uVR7npuY0h2xnSTS9q9B/uO8vFT33EHzedivm+u4PeL9EayXYviqttaq3xAhhReGMgyDzTRU59S5vSGaMznVw6IYs3txX3yUMxFmMoGl1eyutd5GYmBkk3wRG9m8wkG0III/DRadsZWxfQ6EdnRPY2k55oNaqH9gTXQBoZq2ilKYwe3RXJgaJdZY1+kuwWUhIsjMlKDInovYELHozNYmJEWkLIesHor6ZJjjAafdB3+ujcYOaPySDZYWFObjrpDhOvP7CALy1sTa3TR9BGUv620e0LeQ2H4Kwb7WbYU1zDkXOhdfmNiL6XBkytOXyOh17ficcne6ViZ3fQBgoJ4y2muze3rg+nBTJARqSF6spdZeCEm4LyrgWjOVvbYgya603cvoBG34O3Ob3zNDfTaQwIDK5gWdHgIisxOCGifWdssk1gMQnO1DRT2+xhdKQRfYKty6wbKSVX/fpjXtl4osPleoBns/R+1o1y9FHQ6GqvR3eFrrWfrvczMi0BIQS5mU6Kq5qM6Mnjk0bdj2Dy2kT+wawrrGDy8JSwEXuw/to2vRJgaLKDvU9czcXjtDz5i8dnhUhAQ5LtjEh1sDeCaRCbgoq86WhlnFsj+tNVTUhJSE2V1s5Y7Xe/vbOUuT9d3a4SaHd5a0cJQ5MdzMxJjdnUbz2ldaKJgHTTzbcYvaCZHtEPS9UcfYpN61jtUrrx+NpJjFdOzibZbmH1ofMdbBU7vD4/tmDpphtvNvp9kZeZaDy0mttINxmBMRx6urFOUyA4MZsEyQ4L6wsrACJ39E5rl52xxVXNFJ5v6P4ohtEAACAASURBVPDB6fNLpMTQ6JV0M4DoXkSv3YwlDX4jHS4vMxGvX3KmpgUppVHrpi25mc6wEX1tk4cdp6q5YtKQsMcMdvThHgSgpdV1xvSc1E4jerfXz6qD52hoaS3yFrzvJJs2cYrL6zOkhOC650ZnbOAB9/HRcioaXByMUfS9p7iW/Lx0Jg1L6fCtKNb4/bLTukZ6OWqzSSBELCJ67draLWZGZSQwc4iFkWkROPowU1BazSYuHJfJhmO9H9G3e+B14zzo13R0pjO8dNOoSTegvfXWhUg3XiOw+cqicUa7zEmPzNGnBzT6jsqTQOvb8IEz4duz/nCzWZR0E1Pe33uW+U+u7lHxJn0yjbaOrTP0iNrtg5EBRz9uqJYJcKiszijRag3T25+XmUhts6ddx8+6wnJ8fskVk4aGPWZSF9JNJMzISeNERWOHM96/+skJ/uNP21lxoKzjzJ8WL6XVzehlc9r1SwRJNzqR9At0RWWDi9KaZmbmpJKXlUhFg6tPinY9t/YYC3++hnMdTM2nOzghArnTnUSym45Xdnjz6/pwWtC1/ftXLuauybbA4KmOHb3H58frl2H7di6dkEVxVXOHb5GxQs826Uke/anKRjITbaSESa/0+SXVTW4j9z3FYWmXXqm32a8uGsc/v76QR6+5gCnDO85UCibNacUvCTsDXEugiqbu6MvrXZyvb98e9Ahe76vwS6IqNR0tg8bRf7D/LOfrXWwu6n4avxHRRyPdBN2MIwPZEdNGpmKzmNh2osq4mcNJN3oOcNuIdO3h86Q7rcwaFb4QaHBKZbism0iYHtDpgx3vrtPV3Pg/GzhZ0cgfP9E6OBvdPqNIVzD66/KpoCj+VJuIPrhDTqejCEgnkrlBdclpRk4aeYFzeKqXo3q3188fN52i3uXl2Y8Kw67jCZSjBu3B3pEjP3a+njtf2syfN4fvRD4bcOTBbWtYqgOHRTAyLYGyupYOO1U7y9a6JFDyYn1hRUhRu4Ij56Mqg9wV7kCaqa0HZYpPVTYZg5vMJoHVLIyRsdVNbqRsnT+5nXTTZhzBzFFpPHT5eEwd1GVqi9430i4AO1rOnJ+s4oV1RewvrTVkyYNh2rQ76L7vaZppJMSNoy+vd/HQX3ayv6K9I5BSGlOEbSiM/tW0ttnDnS9uZuMxTcvrTkQPrRG93WJm9qg0tp6sMkZLhpNu9JGzbSPhtUfOs2jikLA5vxAa0Sd3IN10he7od5fUGN+9t/cse0pqufn3Gymtaeahy8cBrWV329rQ4PJyOuBgs5LsodJNm/RKIbRjBjv6v249za9WHjE6U/eV1DL9v1awr4u+gz0lNQihPVD1c3gygig12gmmg/lw/1kqGlxMH5nKm9uKKSpvP7LYEyTRWS2mDp2x/qBaeeBcu2Uur4+/7yjh0glZ7XR20DJw/BJjEFwwPr80UoTDOfoxWYmMTEvg5Q0nmPGjlbzw8XEOl9Vx/6vb+O/lhzv59dGhJx8YndLd0KdPVTYZo15BG8ikZ9MYZYUD0k3bzthw0lU0tFawbN3n+sJyHvjjdprcPl7ZeIL9Z2q5KlDHSG/TLq+PTccrjY5YAHuQhNWbZUDixtEnOywsP1DGkar2N2tRRSMVDS4sJsH6gLOOhpUHythUVMmKA+cwCYxe/sjsanW0wRrggjEZ7C+tNTrW2qZXQmvnUHA0uvrgOaqbPNwwa0SHxwzV6LsX0acn2pg6IoUP9p01vtt6ooqhyXaqmzzkZjr59lUX8K0lE/nsjPa2GBF9ZRNOm5nZo9NCfocr0BlrMZswCZg6IoVLJmRReL4et9ePX0qefP8Qv1tzjM/+bgOHztax+tA5PD7JB/vPtjteMHtLahk/JIkku8V4K+oqov/D+iImP76cL7y0mac+PMyvVh6htsnD8fIGvvDS5pDp6P5R6OZrf94Rsv1rm06Rl+nk5fvmYjObeHlD+2wL/S0G6FS60SPArSerjFRKnX/vOUt5vYsHLh0bdls9AydcBP7Nv+5i6W/XA+HTcoUQXDZxCCcqGrGaTfx61VEe++d+pIQ1h8/HrMPQkG66mXGl9/vkBpUrcFjNuAIRfeuMUK0afYPLa7ylNLm97eZxiAa9b0SP6N1eP4+9s59RGQk8fesMztW5qGnycNG4TEZlJBjX842txdz50mZWHCgLepPv2ZtNpMSNo9fnRD1V1/5k6dH85+eNoqi8MerX0BUHzhlP3USbpcuOzGBCIvqgKpTzx2Til9pkBxA+ondYzQxPdYREo3/dVszwVAeLJobX54OPaRLRvX205fa5o9hfWseBM7XUtXg4cKaWO+aP5vm75/DbO2ZjNgn+c8kEvrBgdLttMxJtnK5s4lh5A6MznORlOgPZN9rNpk+MDNoD8Lppw5k6IgWPT3L0XD1ljZIGl5dvLZmIX0o+3K89bEGTrvx+yfL9ZWwpqgyp7y+lZG9JDTNytLECTpuF7BQ7JyoaaXR5efL9g9z03MaQvoe/by/mp+8fYm5uOiXVzSzbcILn1h7jP/60nYf+spNPjlfyy5VHAPjLllP8+7iHD/eXGR2v+0tr2X6qmrsvzGVosoMrJg1lxYFz7TTX4JpGtk464A6erSPZbsHnlxQcbc2C8fklf1hfxMTsJC6bEL6yqN7GStu08bLaFj7cf9bQ9ztKKPjetRfwj69dxAffvBSLSbD9VDVzc9Opb/GyuUiLRv9nTSHzn1zNu3vO8OfNp7j59xs5X9dCk9vL7wuOUdXYeeqhu01nbGejTIurmth2sipEJimu0vp9QiN6k6HR6xF9VlJrYUGPTxoRc3OYmdmiISOQtqkHD3/deppTlU08tnQKN88eaZRRmD4ylanDU40Eg7VHtGv5k/cOGZMChbzZKEcfGdNGpnKyrv2UYltOVJKVZOeLF+UBsKGwNapv8fjYdrKqwx70RpeXdYXl3LUgl1mj0khLjE4KSbJZEAKsJshKai3ZOyc3DYtJsOFYx44eNJ1+S1EVS3+7nq//ZQfrC8u5be6oDmUbALvFhMUkSHZYI9Ydw3HjrBHYzCb+vr2EHaeq8Uu4cEwG104bzsxRaZ1ue8ucHCob3aw7Ws7oDCejM5y4vH5DUgiObj/6P4v42qJxTB2hyUUHz9RxIlAn/brpw8gfnc4H+86y63Q1qQlWDpfV8/TKI3z1zzv4/Iubuem5jcZNsru4hooGN3NyW+3Ly0xk1+lqrnt2PS+tP8Hu4hqWBSLulQfK+P7b+7h0QhZ/fmAB6757OUefvI5n7pjN1pNVHC6r57KJQ1h18By/WXWUJ/51gKwE7ZzuD0gsL284QaLNzO3ztEKtV0/NpqLBxa7ToWUktEg2oNGbw2v0UkoOnKnjM9OHk5VkZ+XBVvnmd2sKOVxWz0OXj+8w2BgecPRtc+n/uasUv4R/fO1ifnDdJBaOzwy3OWlOG/m5GQxLdfD49VOYmZPK8/fkk2A1886uUh55cze/XHkUifaG8Ng7+9l1uoZ/7z3LWztK+MXyI9z/6rZ2k+uEngctf9xuMXPJ+CyWbTgRdvrKd3aVcukv1nLb85v4xl93Gd//8ZOTAWmutfM0wWqmxeOjpsltjDHRNfrkNhMANfVQusnNcDIjJ5X/LTjOsfP1/PajQi4cm8HiC4ZgMZu4b2Ee6U4rFwxLZuqIFE5UNFJc1cSm45XMHp1GaU0zvy/Q5q8IqeLZzZm2IqF3ZhroJ6aNSOGtHZo+mZFo40xNM1uKqig4Us4lE7KYmJ3E8FQHb+0o4ba5Oaw+dJ7/evcApTXNXDFpKF9bPA6/XzJ/TIZxIxUcKcft9XPdtGE8fMV4I1qIFJNJSzVMtPhCbk6nzcL0nFQKDmtP+XCdsaA5qc1FVZhNgjWHzyOA2+fmhF1XRwhBksMSIuF0hzSnjaunZvPPXaWG9DV7dGQzQV42IYtJw5I5XFZPbqaT0YHo669bT1N4rgG/bJXA9Iaem+EkyW5h5+lqymv9OG1mxg1J4srJ2YZG/I0rxvPT9w/xvwXHuWhsJtdOG8YT7x7g79tL+MKC0fzPmmOkOa3cOKt1fpu8zES2nKgiNcHKmw9eyKufnGTZhhOkJlh5avlhpo1M5fm7842h9AA3zByB36+lvl49JZuFT63h2Y8KWTAmg1tHNfPoumb2lNQybmgS7+09w90X5hqprFdMGorNbGLFgTLm5mUY+wzR6DuI6M/UtlDb7GFajtZh//rW06w8UEazx8ezHxVyy+yR3DCzc9kuxWHhnV2lnKtr4VtXTSTFYeGtHcXMzU1n1qg0ZnXxkNb5/LzRfH6e9ra2aOIQ3t5ViknAo9dcwH9cOpbnPz5OZpKN1zadYsX+MiSSzEQb+0pqmP/kR5gEpDqtZFvdDJ9UzwXDko38cb220zN3zOL6323ga3/ZwepvLzKugZSS59Ye44LsZPLz0nl9y2lOVjRytraF1zaf4ksLx4TMX+Cwmjlwpo55T67GYdXGceidprp8Wt/iYUiyPexcy9FgMgkeWzqF21/YxGee3YDDauK/bphq3N9fuWwsX7woF4fVzNIZw/nN6qM89PpOXF4/37xyAq9uPMnyA2UAxgxTAH/fUczYIYncFNR2Y0V8OfpAB+Inxyv41cqjxgCdCUOT+Opl4xBC8NVF43ji3QM8s7qQ59YeY/zQJL5xxXhe+LiINQGn+/zd+Vw7bRgAKw6UkZloY25eBmaTMF4HoyHZYSHD0v5p/d1rJnH3y1uA8OmVAA9cOpbJw1P4/LxR1LV4KKttiSjfV7vhu9cRG8zXF49n47EK3tt7ljmj0yK+QYQQfGXRWL715h5GZyaSG+hveGZ1IUOS7Xx+7ihunhP6wDKZBFdPyea9vWdJtfqYPjIds0mwZPJQ/nv5YSwmwR3zR/PHTScpq23hpzdPY2xWIv/aXcpvPyokzWnlo8Pn+fZVE0MecvPHZLCusJyX753HlBEppDltLD9Qxo/fO8iCMRk8f3e+Mal6MDfNbr3hHvvsFA6eqeOHSyezeeN6RmUksLekhkaXF69fct/Feca6yQ4rF4/P5IN9ZdwyJ4dkh4WPDp3nZGWj8dZiMZtwByK4n39wiH/sLMHrlyyZrHXgTRmews2zR7K3tJYHX9P6A6aPTOUnN03rUjq8dMIQPjlewWubT1HR4OLKyUM5Xt7IU7eE1/Uj4YsX53K2tpkffGayMYHHN6+cAGiJEM9+VIiU8J2rJzJpWAoFR89jFoKaZg8r95/humfX8ds7Zxu/T3/gZSXZ+dnN07n/1W0s319mPKDXFVZQeL6BX902k0sC5Rl++1EhnxyvJDfTyaPXXBBin8Ni5nRVLTaLCZfXz5Aku/HW23Z+hGZPzyJ60NrUjbNG8PHRcl770gImDWt9uxBCGNLY2CFJ3DRrJG/vKsVuMXHR2ExqmtzGQKrgQYP/W3CceXkZ3Dy780CuO8SVo588PAUBPPn+YSoaXPzX9VOYNTqdmTmpxs1x5/zRvLS+iGc/KiQ308mbX7mI1AQrt8zJ4WRFI4+9s58/bTrJtdOG4fNL1hWWc+Wk7E6lkq744sV51J1p3zl30bhMHls6mR/9+6AxyrEt44cmMX5oEqBFLUOTHWHXa4tebqGnTBmRwnvfvJQn/nWAz0wfFtW2188YQXWjh+tnDCfFYeWbV4xn7JAkls4Y3qFUdfdFuby9q5QGFyydo0We44cmkZfpJCtJm4f3RzdMpcntY9wQ7bw8es0k7nxpM1//y06SHRbuDXK6AJ/Lz+GWOSONNnDBsGSevnUmSXYL10zNjqjP5c75of0QM3LS2H6yis1FlVw1Oduokqhz14JcHnxtO9c9u974buH4TB5ZojlHp83MkXN1/GF9ES+sK+KKSUMpr3fx1o4ShIDJw5Nx2iz88f55PPHuAeaPyeDzc0e1q3IajufumgPA/6wp5Jcrj/Le3rPMH5MR8uCKlovHZfGvhy8Ju+y6acN5ZrWWUvrZGSPIy0pkSdDMWe+trOHlYzb+79v7yPmy9sAPLsu9aOIQ8jKd/GnTKT47YwTrC8t5esURhiTbuX7mCGwWE0smD+XtXaUk2S0su++idgGHXgX11vwcHrhkTEg6ZXDtpdomD3XNng7vt2j49e2zcHl9XQ6gfPiK8byzu5QFYzNxWM0smZyNPfBA0lKMtXPh9Uv+M9A+Yk1cOfpEu4VhiYKzDS6unDSU+4Lqt+jYLCYeWzqZx/91gN/fNce44GOyEhmTlcgXFozm6RVHOHa+nroWLzVNHhZfEH4EaqR8ddE4CgrCVwa87+I8LhmfZTitWHHvxXk9jlp0RqYl8Id750a9ncVs4ktB09N9++oLOllbY/aoNKYMT+Hg2TpmBjpUhRD84d55RuSjT1moc9G4TJY/cilna1oYnekMexO3dea35vcsapqVk8b7e7Xsny+HmYLvqinZbPnBlXx8VBvcNm5oEnNz0w07vrVkIl95bTs/ff8Q+bnpvHhPPg0uL7c9vwmL2WQ4jzSnjWfvmN0tG7+2eDw7TlVjEoLffWF22HTMWDAxWwtGEu2WsFMCJtkEv759Fp95dr1R3C643IDJJLj7wlx++v4hrn1mHYXnG0iyW/jJTVONjJQHLxvL3pJafnnbzLAlmPUsovsvzmNsm3tJTzeua/Hy9x3FeP2Sz0wfTvnRsh79brNJRDRKfuyQJH5752zGBM5NssPKFZOG8uH+spCy1ReNzTTelmJNXDl6gNwUE2cbfTx0xfgO17l22nCumTosbCT3+XmjeHZ1IX/85BTpiTZMQhsx2FsIIZiQndz1ilHSNgL9tCCE4IFLx/Ddt/aQn9vaH6C/1XTEpGEpIa/Pvc2MHE2CmTYyhfljMsKuMzTFwW1B848Gc8mELP721Yt44eMivnvtBVjMJtKcNt55aGGP8vmDMZu0GdB6GyEEr9w3zyi/HI4xWYn88raZrDtazhcvzjUkLJ3b8kfxzOpCaps9/ObzM/nM9OEhfSb5uRl88v0rOnz7unrqMMYMSQx7L+kSZl2Lh9c2nyI/N52pI1IpONqdX9s92qYg35qfw/IDZWQl2Ulz2kiyW/g/V0/stePHnaO/Js/KkjkT200W0JaOGkxWkp3P5efw5y2nyEy0a5k2zvATXCt6h1vm5OCoLDSKdQ1EZuSkMWlYMt9aMjGqdNtgpo5I5bd3hkbriXZL2P6CgU4ktdyXzhjO0g4m1E51Wln17ctIcVg7/P2dnefO3tB0jf4vW7Q0yG9f1XsONVKunJzNjseuMnL99z5xdY8y5Lri09eiumBMqpnFHQwmiZTHPzuF3cU1HDpbxxcvinAiTkVMcVp7r9HHggSbmeWPXNbfZsQVvTW7VbLDQqLNzJ7iGublpXPdtPAPm75Gd/JArzp5iENHHwsSbGZevCefJ949wM096MBSKBT9j8VsYvkjl2G3miJOZog3lKPvgFEZTpbdN6+/zVAoFDEg0mkC45W4GhmrUCgUivYoR69QKBRxjnL0CoVCEecoR69QKBRxjnL0CoVCEecoR69QKBRxjnL0CoVCEecoR69QKBRxjnL0CoVCEedENDJWCHEt8CxgBv4gpXyqzfLLgGeAGcAdUsq3gpb5gH2Bj6ellDfEwnCFYrDj8XgoKSmhpUWbNjA1NZVDhw71s1XtUXZFT2e2ORwOcnJysFojr6nfpaMXQpiB54CrgBJgmxDiXSnlwaDVTgP3Ad8Js4tmKeWsiC1SKBQRUVJSQnJyMnl5eQghqK+vJzk59iWve4qyK3o6sk1KSWVlJSUlJYwZ034ehI6IRLqZDxyTUhZJKd3AG8CNbQ5+Ukq5F+i9acwVCkUILS0tZGZmdrtMsuLThxCCzMxM4y0uUiKRbkYCwdMjlQALojiGQwixHfACT0kp32m7ghDiQeBBgOzsbAoKCqLYfSgNDQ092r63UHZFx0C1CwaObampqTQ0NBiffT4f9fX1/WhReJRd0dOVbS0tLdG1QSllp/+AW9F0ef3zPcD/dLDuq8Ctbb4bGfh/LHASGNfZ8fLz82VPWLt2bY+27y2UXdExUO2ScuDYdvDgwZDPdXV1fXr86upq+dxzzxmfv/Od78gpU6bI73znO/1qV6QMVLuk7Nq2ttdeSimB7bIDvxpJRF8KBM+HlhP4LtIHSWng/yIhRAEwGzge6fYKhWJgUlNTw+9//3u+/vWvA/Diiy9SVVWF2dw7c9Mquk8kGv02YIIQYowQwgbcAbwbyc6FEOlCCHvg7yxgIXCw860UCsWnge9///scP36cWbNmcdVVV9HQ0EB+fj5vvvlmf5umaEOXEb2U0iuEeBhYgZZeuUxKeUAI8WO0V4V3hRDzgH8C6cD1QogfSSmnApOBF4QQfrSHylMyNFtHoVDEgB/9+wD7iqtjGk1PGZHCE9dP7XD5U089xf79+9m9ezcASUlJxt+KgUVEefRSyg+AD9p893jQ39vQJJ22230CTO+hjQqFQqHoAWoqQYUiDnji+qkDOi9c0b+oEggKhaJbJCcnD9j0REUoKqJXKBTdIjMzk4ULFzJt2jSuu+66/jZH0QnK0SsUim7z+uuvG38//fTT/WiJojOUdKNQKBRxjnL0CoVCEecoR69QKBRxjnL0CoVCEecoR69QKBRxjnL0CoVCEecoR69QKOKG5uZmFi1ahM/n69b2ixcvZvv27TG2SuPVV1/l4Ycfjnj9O+64g8LCwpgcWzl6hUIRNyxbtoxbbrklLkolf+1rX+MXv/hFTPalHL1CoegWjY2NLF26lJkzZzJt2jTefPNN1qxZw0033WSss2rVKr7whS8AWnXLRx99lKlTp7JkyRK2bt3K4sWLGTt2LO++q1U+f/XVV7nxxhtZvHgxEyZM4Ec/+hEAjz/+OM8884yx3x/+8Ic8++yz7Wz6y1/+wo03ajOdSil59NFHmTZtGtOnTzfKJxcUFLB48WLuueceJk2axF133aVPkmSwbNkyHnnkEePzSy+9xLe+9a2QdXw+H/fdd5+x/9/85jdA6FtBRUUFeXl5xjbFxcXtflu48whw6aWXsnr1arxeb0TXozPUyFiFIh748PsklO4Ccwxv6WHT4bqnOly8fPlyRowYwfvvvw9AbW0tKSkpfP3rX6e8vJwhQ4bwyiuvcM899wCaQ7viiit4+umnufnmm3nsscdYtWoVBw8e5N577+WGG24AYOvWrezfvx+n08m8efNYunQpX/rSl7jlllt45JFH8Pv9vPHGG2zdujXEHrfbTVFRkeFY3377bXbv3s2ePXuoqKhg3rx5XHbZZQDs2rWLLVu2MHHiRBYuXMjGjRu55JJLjH3dfvvtPPnkkzz99NNYrVZeeeUVXnjhhZDj7d69m9LSUvbv3w9oE7F0RbjfdurUqXbnEcBkMjF+/Hj27NlDfn5+l/vuDBXRKxSKbjF9+nRWrVrF9773PdavX09qaipCCO655x7+/Oc/U1NTw6ZNm7jqqqsAsNlsXHvttca2ixYtwmq1Mn36dE6ePGns96qrriIzM5OEhARuueUWNmzYQF5eHpmZmezatYuVK1cye/ZsMjMzQ+ypqKggLS3N+LxhwwbuvPNOzGYz2dnZLFq0iG3btgEwf/58Ro4ciclkYtasWSHHB+3t44orruC9997j8OHDeDwepk8Prbg+duxYioqK+MY3vsHy5ctJSUnp8pyF+23hzqPO0KFDOXPmTNcXowtURK9QxAPXPUVzH5cpnjhxIjt37uSDDz7gscce48orr+Txxx/n/vvv5/rrr8fhcHDbbbdhsWhuxmq1IoQAtGjVbrcbfwfLE/o6bT8/8MADvPrqq5SVlfGlL32pnT0JCQm0tLREZLt+bACz2RxWHnnggQf42c9+xqRJk7j//vvbLU9PT2fPnj2sWLGC559/nr/97W8sW7YMi8WC3+8HaGdPuN8W7jzqMlFLSwsJCQkR/abOUBG9QqHoFmfOnMHpdHL33Xfz6KOPsnPnTgBGjBjBiBEj+OlPfxrWQXbFqlWrqKqqorm5mXfeeYeFCxcCcPPNN7N8+XK2bdvGNddc02679PR0fD6f4VwvvfRS3nzzTXw+H+Xl5axbt4758+dHbMeCBQsoLi7m9ddf584772y3vKKiAr/fz+c+9zl++tOfGr8/Ly+PHTt2APDWW291+ds6Oo8AR48eZdq0aRHb3BEqolcoFN1i3759PProo5hMJqxWK//7v/9rLLvrrrsoLy9n8uTJUdesnz9/Pp/73OcoKSnh7rvvZu7cuYAm/Vx++eWkpaV1mFVz9dVXs2HDBpYsWcLNN9/Mpk2bmDlzJkIIfvGLXzBs2DAOHz4csS233347u3fvJj09vd2y0tJS7r//fiN6//nPfw7Ad77zHW6//XZefPFFli5d2uVvW7FiRdjzeO7cORISEhg2bFjE9naIlHJA/cvPz5c9Ye3atT3avrdQdkXHQLVLyoFj28GDB0M+19XV9ZMl7XnooYfkH/7wBylldHa98sor8qGHHgq7zOfzyZkzZ8qjR492uP2OHTvk3XffHdGxIrFr6dKlcvXq1RHtL5bU1dXJX//618Y5bEvbay+llGhzeIf1q0q6USgUMSU/P5+9e/dy9913x2yfBw8eZPz48Vx55ZVMmDChw/XmzJnD5Zdf3u0BUzo1NTVMnDiRhIQErrzyyh7tq7ukpaVx7733xmRfSrpRKBQxRdenu8N9993Hfffd1+77KVOmUFRUFNE+wnXURktaWhpHjx7t8X56Qnf6NzpCRfQKhUIR5yhHr1B8ipFtRnQq4p/uXHPl6BWKTykOh4PKykrl7AcRUkoqKytxOBxRbac0eoXiU0pOTg4lJSWUl5cD2uCaaB1AX6Dsip7ObHM4HOTk5ES1P+XoFYpPKVarlTFjxhifCwoKmD17dj9aFB5lV/TE2jYl3SgUCkWcoxy9QqFQxDnK0SsUCkWcE5GjF0JcK4Q4IoQ4JoT4fpjllwkhdgohvEKIW9ssu1cIURj4F5thovruTgAAFHJJREFUXgqFQqGImC4dvRDCDDwHXAdMAe4UQkxps9pp4D7g9TbbZgBPAAuA+cATQoj21YEUCoVC0WtEEtHPB45JKYuklG7gDeDG4BWklCellHsBf5ttrwFWSSmrpJTVwCrg2hjYrVAoFIoIiSS9ciRQHPS5BC1Cj4Rw245su5IQ4kHgQYDs7GwKCgoi3H17GhoaerR9b6Hsio6BahcMXNuUXdExUO2C2Ns2IPLopZQvAi8CzJ07Vy5evLjb+9In/h1oKLuiY6DaBQPXNmVXdAxUuyD2tkUi3ZQCo4I+5wS+i4SebKtQKBSKGBCJo98GTBBCjBFC2IA7gHcj3P8K4GohRHqgE/bqwHcKhUKh6CO6dPRSSi/wMJqDPgT8TUp5QAjxYyHEDQBCiHlCiBLgNuAFIcSBwLZVwE/QHhbbgB8HvlMoFApFHxGRRi+l/AD4oM13jwf9vQ1Nlgm37TJgWQ9sVCgUCkUPUCNjFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs5Rjl6hUCjiHOXoFQqFIs6JyNELIa4VQhwRQhwTQnw/zHK7EOLNwPItQoi8wPd5QohmIcTuwL/nY2u+QqFQKLrC0tUKQggz8BxwFVACbBNCvCulPBi02peBainleCHEHcB/A58PLDsupZwVY7sVCoVCESGRRPTzgWNSyiIppRt4A7ixzTo3An8M/P0WcKUQQsTOTIVCoVB0l0gc/UigOOhzSeC7sOtIKb1ALZAZWDZGCLFLCPGxEOLSHtqrUCgUiigRUsrOVxDiVuBaKeUDgc/3AAuklA8HrbM/sE5J4PNxYAFQDyRJKSuFEPnAO8BUKWVdm2M8CDwIkJ2dnf/GG290+wc1NDSQlJTU7e17C2VXdAxUu2Dg2qbsio6Bahd0z7bLL798h5RybtiFUspO/wEXASuCPv8A+EGbdVYAFwX+tgAVBB4ibdYrAOZ2drz8/HzZE9auXduj7XsLZVd0DFS7pBy4tim7omOg2iVl92wDtssO/Gok0s02YIIQYowQwgbcAbzbZp13gXsDf98KrJFSSiHEkEBnLkKIscAEoCiCYyoUCoUiRnSZdSOl9AohHkaL2s3AMinlASHEj9GeIO8CLwOvCSGOAVVoDwOAy4AfCyE8gB/4qpSyqjd+iEKhUCjC06WjB5BSfgB80Oa7x4P+bgFuC7PdP4B/9NBGhUKhUPQANTJWoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6BUKhSLOUY5eoVAo4hzl6Hubxsq+OY6UHS8r2went/SNHQMJvx8++C4Ub+tvSxSKfkU5+p4gJVQeB58n/PJPfge/mgilO3rXjgP/hN9MhXMH2i/zeeGvd8Kya2DDM50/EGJBxTF44TIo+rh3jxMJRWtg6wuw6v/1tyV9h5Tw8S9g/z/625K+weuCna+Bqz62+z3wDhSuiu0++5H4cvS97cSCOfRveHYm/G4OvPUlLXoMprFCu+H8XvjoJ9Hv/9wB+McDWoPrjJZaLWqtK4W3vgz734Z3vg4N5dryIx9AbTEMnwGrn9D+xeI8ueph31vQUtf6nZSw/Htwdg+8dT/Ulvb8OD1h6x+0/09vgpJeftj2F2d2w4ffh+Ya7fO2P8DaJ+Hd/9TaYGec3Ai/vwg2P993947PCxt+owUn7sae7UtKePcb8O7D2j0QDXVnGHquQAvSPC3aedQpXAV/vw/+9kWoOQ3rf6Wd4+C2/inD0t8GxAxPM7x8FSOTFoDnQjCZoXgL1BRD7sWaI6w6AWMXQ+pIbRu/T3NK5w/CyLmQNUFznCc+hqYqcGZCc5XmsBrOgSMVhlwAzizNuWddAHO+CDv/BKsfh0v/D9SfgzO7tAeBuxHy74cdr5BtmwUlSbD1JWiqgCGTYOgUSMoGAQyfDYmZWmMq+DlseQGQsO/vsPNKGH0hODMgcSiMuQzqzkDJNihcCY3lcOUT8NGPNAcLUFEI9/5b20/qaHhgDXz4KGx8VmvcCx+B8wcYWfI+bDsGZjskDdWOU3EMzuzUzunwmZB3KZhMWqMv2we2JFj+fe28OVJh0vWQNhqkD46thnn/AXv+Ci8u0vY3806YeJ22D9B+47n94G4CnwuaqyFjLCQPh4PvQN0Zxp05D9kNkDEGpF/7V3lM+10TroaRc1qvfUsdnN0NWRMheZj2XfUpOLoc5j8Ie96Ajb+Bz70MFjt43Zrzd2ZA9jRt/TM7Ncc3djEMmw7eFu0h6/fA1JvBnqyt5/NoDuLDD8HigPFLtOMeelf7f9h07W+Edo1HzoGS7XBqA0y5GbLGB2yuBWEGqxPOHwg8FKXmvCqOwKlPtLY2fgmMWgAn1mkP7Ck3aW+IFUdh9MXw1zug8Tyc+Jix9olQ+j6MuhBKtsK6X8J1T2kPuaMfwrgrYPRF0FID25fB2p9r52P597R2lpAOwhT4J1r/J/C3PQnGLNauefUJrZ1njIN5X9bO3/lD2rVIzdHazbAZYDJhddfCnjdhZD588qx2vwDYUyD/Pm1dkxlMFu2cSB80nIfELG3/dWfA2wxmG5itWntOGgprfwZ734Ts6bDndbjgWrhgqfZb68tgyo3aeWo4r7UZd6PW7hrLYcX/ZUpTJby4WrsWtafhkm9D7kL4x5e1c19zGpZdq/kOgMPvw+U/gMk3aO3i1Cda2829WLsvT2+GbS/B0Mmw4KuQOARcddpxHanw/9s78+CqqyuOf04SEpZgQlgCypKAYhWlbCNY9kFapHXX1k7r0tJxrDJ1nQ7WqbV2nFY7bcfOOFJaqGuFlkrNjDpAFZSioAEji2wBw2YMlE0gSJZ3+se5L3mJeSEv5P3er2/uZ+ZNfu++3/L9nXvuuefe9/K7B7fD4Z1QNAHOObfBn/i0zPzx5AE7Zvw9HRMTYxANMgtuA2PGjNHS0tLEDzy2H/51pzUIyXAZSpx7y8y2TFsjLX/eHMmwCvjimAUAgJ4XwKxl1jhKZsOHL375uAn3weQ58PRlcHS3lWV3hx5F5oD1p2MvYoEu6hijb4epP4P1z1nDOLonvr7x98L0X0LZy/Y+K9s6os751qinP2bOE4nAa/fBumfbdt9RcvIsIJw80FjWOQ+m/wrKl8Pe960jBOv8frzaOtkP5psDH6+0Rp3dzTqPL44Rt24AOudTX1NNZqQm/j5delhgyMiyhhups/Lcvhacqg9bY7xnA5TOtyyyU1frpE8dgRo31O/ay/Y/ebDx3NndzTdqXcaZmWP3m5EJtdWmv1M3qK+xazRBmt5bVudGn0GsgddWmwYwTbXVX76/giFW55Fa589xfDUnz+p+2c+J1J4i4/xpcP2fYPkv4MMX7P5i6y0zx86pERg6A659xvxrS4nrUNVdSxu3o39PHmjUHbVddSujhs750LWAyJE9ZGhdY/nEB6xDLf2rdextbYexRG0yZhbM+DXMm2odZma21UtzmtQD0PsitudPZuinr0D3QvPbTYvtsx7FcMsSs8nyR2D4zTDmB/DaA9ZRtEa33m4kdYa42jnfEs2amCmnrC7WadzyCitXrmTKlCltsUQDIrJOVce0+FlbAr2IzACeAjKBv6jqb5p9ngM8D4wGDgHfUdUK99lDwCygHviJqi5t7VrtDvSOsiVPMaL7UWuUfYdbUN39rmV6PYfAzrfMWaNBomAI9L2kMVhl5VjWk9cfqg9ZVt+tD2RmmcNXbbKh3aU3WhYLFkD3vQ8VqyxDHzAOcntbMAKoPsyGNxYwfMh5cOGV0CXfhrBHKmzEUHfaAuLR3eaoI79v2U8sNdU2XXKkAnatsPspnmSOFc02Y9m42PRk58KUOU33ObDFhs59L2X1XmX85WOt0zlSYbbqUWRZfKcudq2K1RYceg217LL6kNk2OjICu4dj+8xeXfIby+vrrDHvec8aWqeulo31G2HOnpFpQbRqs2WsX/km9Cji7beWM3lwVwskkkFDkOxRZPPPB7daQ9F6CzgDx9kI41C5XbdLAQyebBlxpB62LzV7nDpiHc6QaXbuvWvt/P2+avuWv2nniNRZ5pbV2fTXnLDziLChdhDDr7vfbLZxsQXkYddC5QbLxoddZ/dWWQafrIKe51u2+dEiyxwzsyF/kGk/XgXnjrBAI2Kv3EK719MnTPPud80feg6xUUafi8z+ZS9asB70NTh9glWr3mHiFTMbfI7//MHut3AYDLsediyFg9usDi66yqbzEiFSb/cUqYe8AeaD+9dbQBx4udVBdq7V4961prvmJHuO1TFw+l2w512z68QH3UgBC4rVh13iVe86bDEfOV5pPnlOfxtN1NfYaOy/283Ol95k9Qb2o4ctJebbxZNsJLj1NbN9bh/4+FVLpAaMtTZeeAkrV69hysTxFgcA1s61EcPIWy1ZikRM84Bxje1/1wrLwDMy4dyR1v72vGdts1sfuOQGOLzLEqCaavO1nFybWsvr3xiDjleZ33UtMI3Fk21E7+joQI+qtvrCgvtOYDCQDXwEXNxsn7uAuW77ZmCR277Y7Z8DFLvzZLZ2vdGjR+vZsGLFirM6Pll4XYkRVl2q4dXmdSVGWHWptk8bUKpx4mpbvoy9DChX1V2qWgMsBK5pts81wHNuezEwTUTElS9U1dOq+glQ7s7n8Xg8noA449SNiNwIzFDVH7n3twBjVXV2zD6b3D773PudwFjgUWCNqr7oyucDb6jq4mbXuAO4A6CwsHD0woUL231DJ06cIDc3t93HJwuvKzHCqgvCq83rSoyw6oL2aZs6dWrcqZtQ/OpGVecB88Dm6BOdm4qlPXNbQeB1JUZYdUF4tXldiRFWXdDx2toydbMfGBDzvr8ra3EfEckC8rAvZdtyrMfj8XiSSFsC/QfABSJSLCLZ2JetJc32KQFuc9s3Am+5LwdKgJtFJEdEioELgPc7RrrH4/F42sIZp25UtU5EZgNLsV/gLFDVzSLyGPYtbwkwH3hBRMqBw1hngNvv78DHQB1wt6rWJ+lePB6Px9MCbZqjV9XXgdeblT0Ss/0FcFOcYx8HHj8LjR6Px+M5C9LrWTcej8fj+RKhewSCiBwEdp/FKXoBZ3iaU0rwuhIjrLogvNq8rsQIqy5on7ZBqtq7pQ9CF+jPFhEpjfdb0lTidSVGWHVBeLV5XYkRVl3Q8dr81I3H4/GkOT7QezweT5qTjoF+XqoFxMHrSoyw6oLwavO6EiOsuqCDtaXdHL3H4/F4mpKOGb3H4/F4YvCB3uPxeNKctAn0IjJDRLaJSLmIzEmhjgEiskJEPhaRzSJyjyt/VET2i0iZe81Mkb4KEdnoNJS6sgIRWS4iO9zfHgFrujDGLmUi8rmI3JsKm4nIAhE54B69HS1r0T5i/NH53AYRGRX/zEnR9VsR2equvURE8l15kYicirHb3GTpakVb3LoTkYeczbaJyDcC1rUoRlOFiJS58sBs1kqMSJ6fxVuR5P/pRRtWwQpQSz9glNvuDmzHVtp6FHgwBLaqAHo1K3sSmOO25wBPpLguPwMGpcJmwCRgFLDpTPYBZgJvYAvFjgPWBqzr60CW234iRldR7H4pslmLdUc7Vp3rSF3NPv8d8EjQNmslRiTNz9Ilo2/LKliBoKqVqrrebR8HtgDntX5UyoldIew54NoUapkG7FTVs/nv6Hajqu9gD+aLJZ59rgGeV2MNkC8i/YLSparLVBtW3V6DPQY8cOLYLB6BrTrXmi4REeDbwMvJuHZrtBIjkuZn6RLozwP2xrzfRwiCq4gUASOBta5otht6LQh6eiQGBZaJyDqxlb0AClW10m1/BhSmRhpgTz6NbXxhsFk8+4TJ736IZX1RikXkQxF5W0QmpkhTS3UXFptNBKpUdUdMWeA2axYjkuZn6RLoQ4eI5AL/BO5V1c+BZ4AhwAigEhs2poIJqjoKuBK4W0QmxX6oNlZMyW9uxdY7uBr4hysKi80aSKV94iEiD2OPAX/JFVUCA1V1JHA/8DcROSdgWaGru2Z8l6YJReA2ayFGNNDRfpYugT5UK1mJSCesAl9S1VcAVLVKVetVNQL8mRQtkq6q+93fA8ASp6MqOhR0fw+kQhvW+axX1SqnMRQ2I759Uu53InI78C3gey444KZFDrntddg8+NAgdbVSd2GwWRZwPbAoWha0zVqKESTRz9Il0LdlFaxAcHN/84Etqvr7mPLYObXrgE3Njw1AWzcR6R7dxr7M20TTFcJuA14NWpujSZYVBps54tmnBLjV/SpiHHAsZuiddERkBvBT4GpVrY4p7y0imW57MLay266gdLnrxqu7MKw6dwWwVVX3RQuCtFm8GEEy/SyIb5mDeGHfTG/HeuKHU6hjAjbk2gCUuddM4AVgoysvAfqlQNtg7BcPHwGbo3YCegJvAjuAfwMFKdDWDVtnOC+mLHCbYR1NJVCLzYXOimcf7FcQTzuf2wiMCVhXOTZ3G/WzuW7fG1z9lgHrgatSYLO4dQc87Gy2DbgySF2u/Fngzmb7BmazVmJE0vzMPwLB4/F40px0mbrxeDweTxx8oPd4PJ40xwd6j8fjSXN8oPd4PJ40xwd6j8fjSXN8oPd4PJ40xwd6j8fjSXP+B8MH58e+CbkuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ttimes, label=\"tf\")\n",
    "plt.plot(stimes, label=\"sympy (only subs)\")\n",
    "plt.legend(); plt.grid();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "p37",
   "language": "python",
   "name": "p37"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}