physgrad-comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Example comparing Different Optimizers\n",
    "\n",
    "The previous section has made many comments about the advantages and disadvantages of different optimization methods. Below we'll show with a practical example how much differences these properties actually make.\n",
    "[[run in colab]](https://colab.research.google.com/github/tum-pbs/pbdl-book/blob/main/physgrad-comparison.ipynb)\n",
    "\n",
    "\n",
    "## Problem formulation\n",
    "\n",
    "We'll consider a very simple setup to clearly illustrate what's happening: we have a two-dimensional input space $\\mathbf{x}$, a mock \"physical model\" likewise with two dimensions $\\mathbf{y}$, and a scalar loss $L$, i.e. \n",
    "$\\mathbf{x} \\in \\mathbb{R}^2$, \n",
    "$\\mathbf{y}:  \\mathbb{R}^2 \\rightarrow  \\mathbb{R}^2 $, and \n",
    "$L:  \\mathbb{R}^2 \\rightarrow  \\mathbb{R} $.\n",
    "The components of a vector like $\\mathbf{x}$ are denoted with $x_i$, and to be in sync with python arrays the indices start at 0.\n",
    "\n",
    "Specifically, we'll use the following  $\\mathbf{y}$ and $L$:\n",
    "\n",
    "$\\quad \\mathbf{y}(\\mathbf{x}) = \\mathbf{y}(x_0,x_1) = \\begin{bmatrix} x_0 \\\\ x_1^2 \\end{bmatrix}$, \n",
    "\n",
    "i.e. $\\mathbf{y}$ only squares the second component of its input, and\n",
    "$L(\\mathbf{y}) = |\\mathbf{y}|^2 = y_0^2 + y_1^2  \\  $ \n",
    "represents a simple squared $L^2$ loss.\n",
    "As starting point for some example optimizations we'll use \n",
    "$\\mathbf{x} = \\begin{bmatrix} \n",
    "  3 \\\\ 3\n",
    "\\end{bmatrix}$ as initial guess for solving the following simple minimization problem: $\\text{arg min}_{\\mathbf{x}} \\ L(\\mathbf{x}).$\n",
    "\n",
    "For us as humans it's quite obvious that $[0 \\ 0]^T$ is the right answer, but let's see how quickly the different optimization algorithms discussed in the previous section can find that solution. And while $\\mathbf{y}$ is a very simple function, it is nonlinear due to its $x_1^2$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 3 Spaces\n",
    "\n",
    "In order to understand the following examples, it's important to keep in mind that we're dealing with mappings between the three _spaces_ we've introduced here:\n",
    "$\\mathbf{x}$, $\\mathbf{y}$ and $L$. A regular forward pass maps an\n",
    "$\\mathbf{x}$ via $\\mathbf{y}$ to $L$, while for the optimization we'll need to associate values\n",
    "and changes in $L$ with positions in $\\mathbf{x}$. While doing this, it will \n",
    "be interesting how this influences the positions in $\\mathbf{y}$ that develop while searching for\n",
    "the right position in $\\mathbf{x}$.\n",
    "\n",
    "```{figure} resources/physgrad-3spaces.jpg\n",
    "---\n",
    "height: 150px\n",
    "name: three-spaces\n",
    "---\n",
    "We're targeting inverse problems to retrieve an entry in $\\mathbf x$ from a loss computed in terms of output from a physics simulator $\\mathbf y$. Hence in a forward pass, we transform from $\\mathbf x$ to $\\mathbf y$, and then compute a loss $L$. The backwards pass transforms back to $\\mathbf x$. Thus, the accuracy in terms of $\\mathbf x$ is the most crucial one, but we can likewise track progress of an optimization in terms of $\\mathbf y$ and $L$.\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "For this example we'll use the [JAX framework](https://github.com/google/jax), which represents a nice alternative for efficiently working with differentiable functions.\n",
    "JAX also has a nice numpy wrapper that implements most of numpy's functions. Below we'll use this wrapper as `np`, and the _original_ numpy as `onp`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.numpy as np\n",
    "import numpy as onp\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start by defining the $\\mathbf{y}$ and $L$ functions, together with a single composite function `fun` which calls L and y. Having a single native python function is necessary for many of the JAX operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting point x = [3. 3.]\n",
      "\n",
      "Some test calls of the functions we defined so far, from top to bottom, y, manual L(y), L(y):\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(DeviceArray([3., 9.], dtype=float32),\n",
       " DeviceArray(90., dtype=float32),\n",
       " DeviceArray(90., dtype=float32))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# \"physics\" function y\n",
    "def physics_y(x):\n",
    "    return np.array( [x[0], x[1]*x[1]] )\n",
    "\n",
    "# simple L2 loss\n",
    "def loss_y(y):\n",
    "    #return y[0]*y[0] + y[1]*y[1] # \"manual version\"\n",
    "    return np.sum( np.square(y) )\n",
    "\n",
    "# composite function with L & y , evaluating the loss for x\n",
    "def loss_x(x):\n",
    "    return loss_y(physics_y(x))\n",
    "\n",
    "\n",
    "x = np.asarray([3,3], dtype=np.float32)\n",
    "print(\"Starting point x = \"+format(x) +\"\\n\")\n",
    "\n",
    "print(\"Some test calls of the functions we defined so far, from top to bottom, y, manual L(y), L(y):\") \n",
    "physics_y(x) , loss_y( physics_y(x) ), loss_x(x) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can evaluate the derivatives of our function via `jax.grad`. E.g., `jax.grad(loss_y)(physics_y(x))` evaluates the Jacobian $\\partial L / \\partial \\mathbf{y}$. The cell below evaluates this and a few variants, together with a sanity check for the inverse of the Jacobian of $\\mathbf{y}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Jacobian L(y): [ 6. 18.]\n",
      "\n",
      "Jacobian y(x): \n",
      "[[1. 0.]\n",
      " [0. 6.]]\n",
      "\n",
      "Sanity check with inverse Jacobian of y, this should give x again: [3. 3.]\n",
      "\n",
      "Gradient for full L(x): [  6. 108.]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# this works:\n",
    "print(\"Jacobian L(y): \" + format(jax.grad(loss_y)(physics_y(x))) +\"\\n\")\n",
    "\n",
    "# the following would give an error as y (and hence physics_y) is not scalar\n",
    "#jax.grad(physics_y)(x) \n",
    "\n",
    "# computing the jacobian of y is a valid operation:\n",
    "J = jax.jacobian(physics_y)(x)\n",
    "print( \"Jacobian y(x): \\n\" + format(J) ) \n",
    "\n",
    "# the code below also gives error, JAX grad needs a single function object\n",
    "#jax.grad( loss_y(physics_y) )(x) \n",
    "\n",
    "print( \"\\nSanity check with inverse Jacobian of y, this should give x again: \" + format(np.linalg.solve(J, np.matmul(J,x) )) +\"\\n\")\n",
    "\n",
    "# instead use composite 'fun' from above\n",
    "print(\"Gradient for full L(x): \" + format( jax.grad(loss_x)(x) )  +\"\\n\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last line is worth a closer look: here we print the gradient $\\partial L / \\partial \\mathbf{x}$ at our initial position. And while we know that we should just move diagonally towards the origin (with the zero vector being the minimizer), this gradient is not very diagonal - it has a strongly dominant component along $x_1$ with an entry of 108.\n",
    "\n",
    "Let's see how the different methods cope with this situation. We'll compare \n",
    "\n",
    "* the first order method _gradient descent_ (i.e., regular, non-stochastic, \"steepest gradient descent\"), \n",
    "\n",
    "* _Newton's method_ as a representative of the second order methods, \n",
    "\n",
    "* and scale-invariant updates from _inverse simulators_.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient descent\n",
    "\n",
    "For gradient descent, the simple gradient based update from equation {eq}`GD-update`\n",
    "in our setting gives the following update step in $\\mathbf{x}$:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\Delta \\mathbf{x}_{\\text{GD}}\n",
    "&= \n",
    "- \\eta ( J_{L} J_{\\mathbf{y}} )^T  \\\\\n",
    "&=\n",
    "- \\eta \\big( \\frac{\\partial L }{ \\partial \\mathbf{y} } \\frac{\\partial \\mathbf{y} }{ \\partial \\mathbf{x} }  \\big)^T\n",
    "\\end{aligned}$$\n",
    "\n",
    "where $\\eta$ denotes the step size parameter .\n",
    "\n",
    "Let's start the optimization via gradient descent at $x=[3,3]$, and update our solution ten times with\n",
    "$\\eta = 0.01$:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GD iter 0: [2.94      1.9200001]\n",
      "GD iter 1: [2.8812    1.6368846]\n",
      "GD iter 2: [2.823576  1.4614503]\n",
      "GD iter 3: [2.7671044 1.3365935]\n",
      "GD iter 4: [2.7117622 1.2410815]\n",
      "GD iter 5: [2.657527  1.1646168]\n",
      "GD iter 6: [2.6043763 1.1014326]\n",
      "GD iter 7: [2.5522888 1.0479842]\n",
      "GD iter 8: [2.501243  1.0019454]\n",
      "GD iter 9: [2.4512184  0.96171147]\n"
     ]
    }
   ],
   "source": [
    "x = np.asarray([3.,3.])\n",
    "eta = 0.01\n",
    "historyGD = [x]; updatesGD = []\n",
    "\n",
    "for i in range(10):\n",
    "    G = jax.grad(loss_x)(x)\n",
    "    x += -eta * G\n",
    "    historyGD.append(x); updatesGD.append(G)\n",
    "    print( \"GD iter %d: \"%i + format(x) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we've already printed the resulting positions in $\\mathbf{x}$, and they seem to be going down, i.e. moving in the right direction. The last point, $[2.451 \\  0.962]$ still has a fair distance of 2.63 to the origin.\n",
    "\n",
    "Let's take a look at the progression over the course of the iterations (the evolution was stored in the `history` list above). The blue points denote the positions in $\\mathbf{x}$ from the GD iterations, with the target at the origin shown with a thin black cross."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc742a69430>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgk0lEQVR4nO3dfZBkdX3v8fd3Zh/Z3ZnVyM4uQVzIJuvyUAGhDA8pBGGDUSPgFU1pcuGaXEGwgl7Qy9YtH2IUohHFhAcpTSJ7TRnKKgsLn5YlBr2BFRWBSFw3WWUJq/sAqDOsMLMP871/dA9pZrtnpne653TPeb+qTu32Ob/T/Z2zvZ8+8zu//p3ITCRJ5dBTdAGSpJlj6EtSiRj6klQihr4klYihL0klYuhLUokY+pJUIoa+JJXInKILmGkREcARwNNF1yJJLbQE+FlO8o3b0oU+lcDfXnQRktQGRwI/nahBGUP/aYDHH3+cvr6+omuRpGkbGhrixS9+MUyhB6OMoQ9AX1+foS+pdLyQK0klYuhLUokY+pJUIqXt059IZrJ//34OHDhQdCkdo7e3lzlz5lAZ8SqpWxn64+zdu5cdO3bwzDPPFF1KxznssMNYsWIF8+bNK7oUSYfI0K8xOjrKo48+Sm9vL0cccQTz5s3zzJbKbz579+7liSee4NFHH+U3f/M36emxZ1Bqh8xsa+4UGvoR8Xbg7cDK6qp/Az6YmV+bYJ+LgL+o7vMfwP/OzK+2op69e/cyOjrKi1/8Yg477LBWPOWssXDhQubOnctjjz3G3r17WbBgQdElSbPGnpH9fGzDFu7evIt9B5K5vcG5awa4+rzVLJ7f2pgu+kx/O3ANlfAO4GLgSxFxUmb+2/jGEXE68HlgHfBl4M3AHRHxssx8pFVFeRZbn8dFar09I/t5/c33snX3HkZrJlBYv2kb9/34Sb54+RktDf5C/xdn5p2Z+dXM/I/M/PfM/D/AHuDUBrtcCXw9M/8qMzdn5nuB7wPvmKmaJamVPrZhy0GBDzCasHX3Hq7fsKWlr9cxp24R0RsRfwgsAjY1aHYacPe4dRuq6xs97/yI6BtbqExKJEkd4e7Nuw4K/DGjCRs372rp6xUe+hFxQkTsAUaATwEXZuYPGzRfDow/Aruq6xtZBwzWLE62JqkjZCb7Dkw4KSb7DySTTJzZlMJDH9gCnAj8DnALcFtEHNvC578O6K9Zjmzhc3eUnTt3cuWVV7Jq1SoWLFjAwMAAZ5xxBrfccstzQ1BXrlxJRBARLFy4kJUrV/LGN76Rb3zjGwVXL5VPRDC3d+KROnN6o6WjeQoP/czcm5lbM/OBzFwHPEyl776encDAuHUD1fWNnn8kM4fGFmZ4Hv1WfkJP5Cc/+QknnXQSd911F9deey0PPvggmzZt4j3veQ9f/vKXufvu/+oV++AHP8iOHTvYsmUL69evZ+nSpZx77rl8+MMfnpFaJf2Xc9cM0NMg03sC1q4ZH3nTU/TonXp6gPkNtm0CzgFuqFm3lsbXAAoxk8Ovxlx++eXMmTOH733veyxatOi59ccccwznn3/+8z58lixZwvLllR6xo446ijPPPJMVK1bwvve9jze84Q2sXr26LTVKOtjV563mvh8/edDF3J6AVcsWc9V5rf3/WOiZfkRcFxFnRsTKat/+dcBZwD9Ut6+vrhvzSeBVEXFVRLw0Ij4AnALcONO1NzI2/Gr9pm1s/8Wz7BoaZvsvnmX9pm28/uZ72TOyv+Wv+dRTT3HXXXdxxRVXPC/wa0326+GVV15JZvKlL32p5fVJamzx/Dl88fIzuPi0lRz5goUs71vAkS9YyMWnrWz5cE0o/kx/GbAeWEHlIuu/Audl5sbq9qOA0bHGmXlfRLwZ+BBwLZXx/Re0coz+dE1l+NX7X3dcS19z69atZOZBZ+gvetGLGB4eBuCKK67gIx/5SMPneOELX8iyZcvYtm1bS2uTNLnF8+fw/tcdx/tfd9zs/kZuZv7JJNvPqrPuC8AX2lXTdE1l+FWrQ7+R73znO4yOjvKWt7yFkZGRSdu3+80maXLt/j9Y9Jn+rNLM8KtW/sOuWrWKiGDLlud/ieOYY44BKlMoTOapp57iiSee4Oijj25ZXZI6T+Gjd2aTIoZfAfzar/0aa9eu5cYbb+RXv/rVIT3HJz/5SXp6erjgggtaWpukzmLot9hMD78ac/PNN7N//35OOeUUbr/9djZv3syWLVv43Oc+x49+9CN6e3ufa/v000+zc+dOHn/8cb71rW/xtre9jQ996EN8+MMfZtWqVW2pT1JniJkaR94pqlMxDA4ODh50Y/Th4WEeffRRjj766EOeRbLR5Eljw6/acTV+zI4dO7j22mv5yle+wvbt25k/fz7HHnssF110EZdffjmHHXYYK1eu5LHHHgNg3rx5LF++nFNPPZXLLruMs88+e8Lnb8XxkdR6Q0ND9Pf3A/RXv4/UkKFfo1WhtmdkP9dv2MLGzbvYfyCZ0xusXTPAVW0cpz8TDH2pMzUT+t2bQB1sJodfSVIz7NNvMwNfUicx9CWpRAx9SSoRQ7+Osl3cniqPi9T9DP0ac+fOBXhu7nk939hxGTtOkrqPo3dq9Pb2snTpUnbv3g3AYYcd5oVYKmf4zzzzDLt372bp0qXP+6KXpO5i6I8zNs/8WPDrvyxduvS54yOpOxn640QEK1asYNmyZezbt6/ocjrG3LlzPcOXZgFDv4He3l5DTtKs44VcSSoRQ1+SSsTQl6QSMfQlqUQMfUkqEUNfkkrE0JekEjH0JalEDH1JKhFDX5JKxNCXpBIx9CWpRAx9SSoRQ1+SSsTQl6QSMfQlqUQMfUkqEUNfkkrE0JekEjH0JalECg39iFgXEd+NiKcjYndE3BERqyfZ55KIyHHL8EzVLEndrOgz/VcANwGnAmuBucBdEbFokv2GgBU1y0vaWaQkzRZzinzxzHxV7eOIuATYDZwMfGviXXNnG0uTpFmp6DP98fqrf/58knaLI+KxiHg8Ir4UEcc1ahgR8yOib2wBlrSsWknqMh0T+hHRA9wA3JuZj0zQdAvwVuB84I+o/Az3RcSRDdqvAwZrlu2tqlmSuk1kZtE1ABARtwC/D/xuZk45mCNiLrAZ+HxmvrfO9vnA/JpVS4Dtg4OD9PX1TbNqSSre0NAQ/f39AP2ZOTRR20L79MdExI3Aa4Ezmwl8gMzcFxEPAqsabB8BRmpeazqlSlJXK3rIZlQD/0LglZn56CE8Ry9wArCj1fVJ0mxT9Jn+TcCbqfTPPx0Ry6vrBzPzWYCIWA/8NDPXVR+/D/g2sBVYCrybypDNz8xs6ZLUfYoO/bdX/7xn3Pr/AXy2+vejgNGabS8APg0sB34BPACcnpk/bFuVkjRLdMyF3JlSHbY56IVcSbNFMxdyO2bIpiSp/Qx9SSoRQ1+SSsTQl6QSMfQlqUQMfUkqEUNfkkrE0JekEjH0JalEDH1JKhFDX5JKxNCXpBIx9CWpRAx9SSoRQ1+SSsTQl6QSMfQlqUQMfUkqEUNfkkrE0JekEjH0JalEDH1JKhFDX5JKxNCXpBIx9CWpRAx9SSoRQ1+SSsTQl6QSMfQlqUQMfUkqEUNfkkrE0JekEjH0JalEDH1JKhFDX5JKpNDQj4h1EfHdiHg6InZHxB0RsXoK+10UET+KiOGI+EFEvHom6pWkblf0mf4rgJuAU4G1wFzgrohY1GiHiDgd+Dzwt8BJwB3AHRFxfNurlaQuF5lZdA3PiYjDgd3AKzLzWw3a3A4syszX1qz7NvBQZl42hdfoAwYHBwfp6+trUeWSVJyhoSH6+/sB+jNzaKK2RZ/pj9df/fPnE7Q5Dbh73LoN1fUHiYj5EdE3tgBLpl+mJHWnjgn9iOgBbgDuzcxHJmi6HNg1bt2u6vp61gGDNcv26VUqSd2rY0KfSt/+8cAftvh5r6PyG8TYcmSLn1+SusacogsAiIgbgdcCZ2bmZGfiO4GBcesGqusPkpkjwEjNa02jUklqncyc8UwqNPSj8tP+DXAhcFZmPjqF3TYB51DpChqztrpekjranpH9fGzDFu7evIt9B5K5vcG5awa4+rzVLJ7f/kgu+kz/JuDNwPnA0xEx1i8/mJnPAkTEeuCnmbmuuu2TwDcj4irgK1S6g04B3jajlUtSk/aM7Of1N9/L1t17GK0ZOLl+0zbu+/GTfPHyM9oe/EX36b+dSj/7PcCOmuVNNW2OAlaMPcjM+6h8ULwNeBh4A3DBJBd/JalwH9uw5aDABxhN2Lp7D9dv2NL2Ggo908/MSTuzMvOsOuu+AHyhHTVJUrvcvXnXQYE/ZjRh4+ZdvP91x7W1hqLP9CWpFDKTfQcm/jLs/gNJu78wa+hL0gyICOb2Tty5Mac32j6ax9CXpBly7poBehpkek/A2jXjR6O3nqEvSTPk6vNWs2rZ4oOCvydg1bLFXHXepJMMT1tHTbg2E5xwTVKR9ozs5/oNW9i4eRf7DyRzeoO1awa4ahrj9JuZcM3Ql6SCtOobud08y6YklUYR08IY+pJUIoa+JJWIoS9JJWLoS1KJGPqSVCKGviSViKEvSSVi6EtSiRj6klQihr4klYihL0klYuhLUokY+pJUIoa+JJWIoS9JbdZJ9y05tNu0SJImtGdkPx/bsIW7N+9i34Fkbm9w7poBrp7GHbJawdCXpBbbM7Kf1998L1t372G05iR//aZt3PfjJ/ni5WcUFvx270hSi31sw5aDAh9gNGHr7j1cv2FLMYVh6EtSy929eddBgT9mNGHj5l0zW1ANQ1+SWigz2Xdg4gu3+w9kYRd3DX1JaqGIYG7vxDc8n9MbhdwUHQx9SWq5c9cM0NMg03sC1q4ZmNmCal+/sFeWpFnq6vNWs2rZ4oOCvydg1bLFXHXe6mIKA6KTvjQwEyKiDxgcHBykr6+v6HIkzVJ7RvZz/YYtbNy8i/0Hkjm9wdo1A1zVhnH6Q0ND9Pf3A/Rn5tBEbQ19SWqzzGxrH34zoW/3jiS1WVEXbesx9CWpRFoW+hGxJiJ+0uQ+Z0bEnRHxs4jIiLhgkvZnVduNX5ZPq3hJKolWnunPA17S5D6LgIeBK5rcbzWwombZ3eT+klRKU76EHBEfn6TJ4c2+eGZ+Dfha9fmb2XV3Zv6y2deTpLJrZtzQlcBDQKMrw4unXc3UPRQR84FHgA9k5r0z+NqS1FC7R+pMVzOhvxX4RGZ+rt7GiDgReKAVRU1gB3AZ8D1gPvCnwD0R8TuZ+f0Gdc2vth2zpM01SiqZTp07v55mqvkecDJQN/SBBNr68ZaZW4DaOUnvi4jfAN4F/HGD3dYB729nXZLKq5Pnzq+nmQu5VwE3NNqYmQ9nZhFDQL8DrJpg+3VAf81y5EwUJakcOnnu/HqmHNKZuTMzH4uIsxu1iYhLW1NWU06k0u1TV2aOZObQ2AI8PWOVSZr1Onnu/HoO5cz86xHxVxExd2xFRLwoIu4E/rKZJ4qIxRFxYvV6AMDR1cdHVbdfFxHra9q/MyLOj4hVEXF8RNwAvBK46RB+Dkmalk6fO7+eQwn9s4ELge9GxLER8Roqo2j6qJx1N+MU4MHqAvDx6t8/WH28Ajiqpv084HrgB8A3gd8Gzs3Mf2r+x5Ck6en0ufPrafrqQmbeVz0z/xTwfSofHO8FPppNfpxl5j1McPE3My8Z9/ijwEebq1iS2ufcNQOs37StbhdP0XPn13OoF15/i8pZ+nZgP5VvyB7WqqIkqVt08tz59TQd+hFxDbAJ2AgcD7wcOAn414g4rbXlSVJnWzx/Dl+8/AwuPm0lR75gIcv7FnDkCxZy8WkrO264JhzCfPoRsQN4a3UKhbF1c4FrgT/LzPkNd+4AzqcvqZ2K+EZuM/PpH8pH0AmZ+WTtiszcB7w7Ir58CM8nSbNGJ120rafp7p3xgT9u2zenV44kdYdOGobZjM7qbJKkDtZNc+w00h1VSlLBum2OnUa8XaIkTUG3zbHTiKEvSVPQbXPsNGLoS9IkunGOnUYMfUmaRDfOsdOIoS9JU3DumoGDploY04lz7DRi6EvSJDKz6+bYaaTzxxdJUgHqjcl/xW8dzstXvpB7/v0J9h9I5vQGa9cMcJXj9CWpezUak//57/wnq5Yt5uvvPJNF83q7og9/PLt3JGmcqYzJ78bAB0Nfkg4yW8bk12PoS1KN2TQmvx5DX5JqzKYx+fUY+pJUIzNnzZj8ehy9I6n0xg/P7O0J+hbMYfDZ/dR24nTbmPx6DH1JpdZoeGZPQP/CuRw2v5fRUbpyTH493Vu5JLXARMMzh4b3ceFJv877/uDYru3DH88+fUmlNpXhmbMl8MHQl1Ris314Zj2GvqRSm83DM+uxT19SqYwfqfOrkX0N23b78Mx6DH1JpdFopE49s2F4Zj2GvqTSaDRSZ8yieb0sWTB31gzPrGd2/TSSNIGJRuoAvGDRPP7fe86eVX3443khV1IpjI6OTmmkzmznmb6kWWv8Rdsn94xM2H62jdSpx9CXNCs1c9EWZudInXrs3pE0K0120bbWbB2pU49n+pJmncyc9KJtb09w+OL5s3qkTj2F/oQRcSbwbuBkYAVwYWbeMck+ZwEfB44DHgc+lJmfbWedkjrf8/vvR3ni6Yn77w9fPJ/7rjmbnp5ydXgU/bG2CHgY+Dvgi5M1joijga8AnwLeApwDfCYidmTmhnYWKqlzNdt/D5WLtmULfCg49DPza8DXgKleMb8MeDQzr6o+3hwRvwu8CzD0pZJqpv8eynPRtp5u+5g7Dbh73LoN1fV1RcT8iOgbW4Al7SxQ0szb+MOdTQV+WS7a1lN0906zlgO7xq3bBfRFxMLMfLbOPuuA97e9Mkkz6rk+/B/uYsfg8IRte6LShz93Tk+pLtrWU4af+joqF37HLAG2F1SLpBZotg//iKULZ/30ClPVbaG/ExjfETcADDU4yyczR4DnLuP7jy51v2bH4K9dM+D//apu69PfRGXETq211fWSSmKyMfhjyt5/X0/R4/QXA6tqVh0dEScCP8/M/4yI64Bfz8z/Xt3+KeAdEfFRKsM8Xwm8EXjNDJYtqUBTucVhT8CK/gX83rHLS91/X0/RR+IU4J9rHo/1vd8GXELlC1tHjW3MzEcj4jXAJ4ArqfTN/6lj9KXyiIhJb3F4RP9C/uWaV85QRd2l6HH69wAN//Uy85IG+5zUtqIkdbxz1wywftO2ul08PQFrjy3nGPyp6LY+fUni6vNWs2rZYnrGnTLahz+5yJz9Nw2oVf2C1uDg4CB9fX1FlyPpEO0Z2c/1G7awcfMu9h/I0k2cVmtoaIj+/n6A/swcmqitoS+p62VmqYdkNhP6du9I6nplDvxmGfqSVCKGviSViKEvSSVi6EtSiRj6klQihr4klYihL0klYuhLUokY+pJUIoa+JJWIoS9JJWLoS1KJGPqSVCKGviSViKEvSSVi6EtSiRj6klQihr4klYihL0klYuhLUokY+pJUIoa+JJWIoS9JJWLoS1KJGPqSVCKGviSViKEvSSVi6EtSiRj6klQihr4klYihL0kl0hGhHxFXRMS2iBiOiPsj4uUTtL0kInLcMjyT9UpStyo89CPiTcDHgT8HXgY8DGyIiGUT7DYErKhZXtLuOiVpNig89IH/BXw6M/8+M38IXAY8A7x1gn0yM3fWLLtmpFJJ6nKFhn5EzANOBu4eW5eZo9XHp02w6+KIeCwiHo+IL0XEcRO8xvyI6BtbgCWtql+Suk3RZ/ovAnqB8Wfqu4DlDfbZQuW3gPOBP6LyM9wXEUc2aL8OGKxZtk+zZknqWkWHftMyc1Nmrs/MhzLzm8DrgSeASxvsch3QX7M0+nCQpFlvTsGv/yRwABgYt34A2DmVJ8jMfRHxILCqwfYRYGTscUQcWqWSNAsUeqafmXuBB4BzxtZFRE/18aapPEdE9AInADvaUaMkzSZFn+lDZbjmbRHxPeA7wDuBRcDfA0TEeuCnmbmu+vh9wLeBrcBS4N1Uhmx+ZqYLl6RuU3joZ+btEXE48EEqF28fAl5VMwzzKGC0ZpcXAJ+utv0Fld8UTq8O95QkTSAys+gaZlR12Obg4OAgfX19RZcjSdM2NDREf38/QH9mDk3UtutG70iSDp2hL0klYuhLUokY+pJUIoa+JJWIoS9JJWLoS1KJGPqSVCKGviSViKEvSSVi6EtSiRj6klQihr4klYihL0klYuhLUokY+pJUIoa+JJWIoS9JJWLoS1KJGPqSVCKGviSViKEvSSVi6EtSiRj6klQihr4klYihL0klYuhLUokY+pJUIoa+JJWIoS9JJWLoS1KJGPqSVCKGviSViKEvSSVi6EtSiXRE6EfEFRGxLSKGI+L+iHj5JO0viogfVdv/ICJe3Y667rzzToaHh+tuGx4e5s4772zHy0pS2xQe+hHxJuDjwJ8DLwMeBjZExLIG7U8HPg/8LXAScAdwR0Qc3+ra1q5dy6233npQ8A8PD3Prrbeydu3aVr+kJLVVZGaxBUTcD3w3M99RfdwDPA78TWb+ZZ32twOLMvO1Neu+DTyUmZdN4fX6gMHBwUH6+vomrW8s4C+99FIWLFhw0GNJKtrQ0BD9/f0A/Zk5NFHbQs/0I2IecDJw99i6zBytPj6twW6n1bav2tCofUTMj4i+sQVY0kyNCxYs4NJLL+XWW2/ll7/8pYEvqavNKfj1XwT0ArvGrd8FvLTBPssbtF/eoP064P2HWiBUgv/iiy/m5JNP5oEHHjDwJXWtwvv0Z8B1QH/NcmSzTzA8PMxtt93GAw88wG233dbw4q4kdbqiQ/9J4AAwMG79ALCzwT47m2mfmSOZOTS2AE83U2BtH/7SpUuf6+ox+CV1o0JDPzP3Ag8A54ytq17IPQfY1GC3TbXtq9ZO0P6Q1btoW9vHb/BL6jZFn+lDZbjm/4yIiyNiDXALsAj4e4CIWB8R19W0/yTwqoi4KiJeGhEfAE4Bbmx1YRs3bqx70XYs+Ddu3Njql5Sktip8yCZARLwDeDeVi7EPAX+WmfdXt90DbMvMS2raXwR8CFgJ/Afwnsz86hRfq6khm5LU6ZoZstkRoT+TDH1Js03XjNOXJM0sQ1+SSsTQl6QSMfQlqUQMfUkqEUNfkkqk6AnXCjM0NOGoJknqGs3kWRnH6f86sL3oOiSpDY7MzJ9O1KCMoR/AETQ58RqVefi3U5mls9l9i2btxbD2YpS19iXAz3KSUC9d9071gEz4SVhP5bMCgKcn+8Zbp7H2Ylh7MUpc+5TaeyFXkkrE0JekEjH0p24E+PPqn93G2oth7cWw9gmU7kKuJJWZZ/qSVCKGviSViKEvSSVi6EtSiRj6NSLiiojYFhHDEXF/RLx8kvYXRcSPqu1/EBGvnqla69Qy5doj4pKIyHHL8EzWW1PLmRFxZ0T8rFrHBVPY56yI+H5EjETE1oi4pP2V1q2jqdqrdY8/7hkRy2eo5LE61kXEdyPi6YjYHRF3RMTqKexX+Pv9UGrvlPd7RLw9Iv41Ioaqy6aI+P1J9mn5MTf0qyLiTcDHqQyXehnwMLAhIpY1aH868Hngb4GTgDuAOyLi+Bkp+Pm1NFV71RCwomZ5SbvrbGARlXqvmErjiDga+Arwz8CJwA3AZyLivDbVN5Gmaq+xmucf+90trmsyrwBuAk4F1gJzgbsiYlGjHTro/d507VWd8H7fDlwDnAycAnwD+FJEHFevcduOeWa6VIat3g/cWPO4h8p0Ddc0aH878OVx674NfKoLar8E+GXRx7xOXQlcMEmbjwCPjFv3j8DXu6D2s6rtlhZ9rMfVdXi1rjMnaNMx7/dDqL0j3+/V2n4O/MlMHnPP9IGImEfl0/fusXWZOVp9fFqD3U6rbV+1YYL2bXGItQMsjojHIuLxiGh4ttGBOuK4T9NDEbEjIjZGxBlFFwP0V//8+QRtOvW4T6V26LD3e0T0RsQfUvltcVODZm055oZ+xYuAXmDXuPW7gEb9rcubbN8uh1L7FuCtwPnAH1F5H9wXEUe2q8gWanTc+yJiYQH1NGMHcBnw36rL48A9EfGyogqKiB4qXWT3ZuYjEzTtlPf7c5qovWPe7xFxQkTsofKN208BF2bmDxs0b8sxL90sm4LM3ETN2UVE3AdsBi4F3ltUXbNdZm6hEkBj7ouI3wDeBfxxMVVxE3A88LsFvf50TKn2Dnu/b6FyLaofeANwW0S8YoLgbznP9CueBA4AA+PWDwA7G+yzs8n27XIotT9PZu4DHgRWtba0tmh03Icy89kC6pmu71DQcY+IG4HXAmdn5mQ3FuqU9zvQdO3PU+T7PTP3ZubWzHwgM9dRGQhwZYPmbTnmhj6VfwjgAeCcsXXVXx3PoXF/26ba9lVrJ2jfFodY+/NERC9wApXuh07XEce9hU5kho97VNwIXAi8MjMfncJuHXHcD7H28c/RSe/3HmB+g23tOeZFX73ulAV4EzAMXAysAW4FfgEMVLevB66raX86sA+4Cngp8AFgL3B8F9T+PuD3gGOoDPH8PPAscGwBtS+mEnwnUhmF8a7q34+qbr8OWF/T/mjgV8BHq8f9cmA/cF4X1P5OKv3Kq6h0S9xA5be0c2a47puBX1IZ/ri8ZllY06Yj3++HWHtHvN+r74czgZVUPnSuA0aBtTN5zGf0P0mnL8A7gMeoXGS5H/idmm33AJ8d1/4iKn10I8AjwKu7oXbgEzVtd1IZ935SQXWfVQ3M8ctnq9s/C9xTZ58Hq/X/GLikG2oH3gNsrQbOU1S+a3B2AXXXqzlrj2Onvt8PpfZOeb9TGW+/rVrHbiojc9bO9DF3amVJKhH79CWpRAx9SSoRQ1+SSsTQl6QSMfQlqUQMfUkqEUNfkkrE0JekEjH0pTbrlNs7SmDoS23VYbd3lJyGQZqOiDgc+AHw15l5bXXd6VTmUfl9KhN9vSYzj6/Z5x+p3DLxVTNfscrOM31pGjLzCSp3ZfpARJwSEUuA/0vlnsX/ROfeZlAl5Z2zpGnKzK9GxKeBfwC+R2Xq53XVzRPe3jG788Yv6mKe6UutcTWVk6iLgLdk5kjB9Uh1GfpSa/wGcASV/1Mra9bPtts7qsvZvSNNU0TMAz4H3E7lhhefiYgTMnM3lVvbvXrcLt18e0d1OUfvSNMUEX8FvAH4bWAP8E1gMDNfWx2y+QhwE/B3wCuBv6YyomdDQSWrxAx9aRoi4ixgI5XbHv5Ldd1K4GHgmsy8pdrmE8CxwHbgLzLzszNfrWToS1KpeCFXkkrE0JekEjH0JalEDH1JKhFDX5JKxNCXpBIx9CWpRAx9SSoRQ1+SSsTQl6QSMfQlqUQMfUkqkf8PBlq4U/dgIEYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
    "historyGD = onp.asarray(historyGD)\n",
    "updatesGD = onp.asarray(updatesGD) # for later\n",
    "axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='#1F77B4', label='GD')\n",
    "axes.scatter([0], [0], lw=0.25, color='black', marker='x') # target at 0,0\n",
    "axes.set_xlabel('x0'); axes.set_ylabel('x1'); axes.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No surprise here: the initial step mostly moves downwards along $x_1$ (in top right corner), and the updates afterwards curve towards the origin. But they don't get very far. It's still quite a distance to the solution in the bottom left corner."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton\n",
    "\n",
    "For Newton's method, the update step is given by\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\Delta \\mathbf{x}_{\\text{QN}} &= \n",
    "- \\eta \\left( \\frac{\\partial^2 L }{ \\partial \\mathbf{x}^2 }  \\right)^{-1}\n",
    "  \\frac{\\partial L }{ \\partial \\mathbf{x} }\n",
    "\\\\\n",
    "&=\n",
    "- \\eta \\ H_L^{-1} \\ ( J_{L} J_{\\mathbf{y}} )^T\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Hence, in addition to the same gradient as for GD, we now need to evaluate and invert the Hessian of $\\frac{\\partial^2 L }{ \\partial \\mathbf{x}^2 }$.\n",
    "\n",
    "This is quite straightforward in JAX: we can call `jax.jacobian` two times, and then use the JAX version of `linalg.inv` to invert the resulting matrix.\n",
    "\n",
    "For the optimization with Newton's method we'll use a larger step size of $\\eta =1/3$. For this example and the following one, we've chosen the step size such that the magnitude of the first update step is roughly the same as the one of GD. In this way, we can compare the trajectories of all three methods relative to each other. Note that this is by no means meant to illustrate or compare the stability of the methods here. Stability and upper limits for $\\eta$  are separate topics. Here we're focusing on convergence properties.\n",
    "\n",
    "In the next cell, we apply the Newton updates ten times starting from the same initial guess:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Newton iter 0: [2.        2.6666667]\n",
      "Newton iter 1: [1.3333333 2.3703704]\n",
      "Newton iter 2: [0.88888884 2.1069958 ]\n",
      "Newton iter 3: [0.59259254 1.8728852 ]\n",
      "Newton iter 4: [0.39506167 1.6647868 ]\n",
      "Newton iter 5: [0.26337445 1.4798105 ]\n",
      "Newton iter 6: [0.17558296 1.315387  ]\n",
      "Newton iter 7: [0.1170553 1.1692328]\n",
      "Newton iter 8: [0.07803687 1.0393181 ]\n",
      "Newton iter 9: [0.05202458 0.92383826]\n"
     ]
    }
   ],
   "source": [
    "x = np.asarray([3.,3.])\n",
    "eta = 1./3.\n",
    "historyNt = [x]; updatesNt = []\n",
    "\n",
    "Gx = jax.grad(loss_x)\n",
    "Hx = jax.jacobian(jax.jacobian(loss_x))\n",
    "for i in range(10):\n",
    "    g = Gx(x)\n",
    "    h = Hx(x)\n",
    "    hinv = np.linalg.inv(h)\n",
    "    \n",
    "    x += -eta * np.matmul( hinv , g )\n",
    "    historyNt.append(x); updatesNt.append( np.matmul( hinv , g) )\n",
    "    print( \"Newton iter %d: \"%i + format(x) )\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last line already indicates: Newton's method does quite a bit better. The last point $[0.052 \\ 0.924]$ only has a distance of 0.925 to the origin (compared to 2.63 for GD)\n",
    "\n",
    "Below, we plot the Newton trajectory in orange next to the GD version in blue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc7428c5bb0>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqEUlEQVR4nO3df3RU5b3v8fc3Qwg/QgIWIXiAoodK+WEPVGpVKqCQg9gef9badWwrt6dXLbrUHqot6/qrakVOW2utWi32FLGnHpf3enFp9SDoIVahglapXjCWFim0EKBKEqoJkHzvH3smDsnMJAMzs2eyP6+19gqz59kz3wyTT3ae/czzmLsjIiLRUBZ2ASIiUjgKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0BcRiRCFvohIhCj0RUQipE/YBRSamRlwDNAcdi0iIjk0CPiLd/OJ28iFPkHgbw+7CBGRPBgJ/DlTgyiGfjPAtm3bqKqqCrsWEZEj1tTUxKhRo6AHPRhRDH0AqqqqFPoiEjm6kCsiEiEKfRGRCFHoi4hESGT79DNxdw4ePEhbW1vYpQgQi8Xo06cPwWhbETkSCv1O9u/fz44dO3j//ffDLkWSDBgwgBEjRtC3b9+wSxEpaQr9JO3t7WzZsoVYLMYxxxxD3759dXYZMndn//797N69my1btvCxj32MsjL1Skrv5e3tWB7f46GGvpl9Hfg6MCa+6/8Bt7j7MxmOuRC4NX7M74FvufvTuahn//79tLe3M2rUKAYMGJCLh5Qc6N+/P+Xl5WzdupX9+/fTr1+/sEsSyal9Te/x5sPX8tE9q4l5G20WY+vQmUz68veorBqS0+cK+5RpO/Bt4ERgKvA88ISZTUzV2MxOBR4BfgZMAZYDy81sUi6L0plk8dH/ifRW+5reY89d0/nUrv/NCN/NMN5lhO9m6q7/w567prOv6b2cPl+oP0nu/qS7P+3uv3f3t939fwH7gJPTHHI18F/u/j133+TuNwC/Ba4sVM0iIrn05sPXMqptGzE7dMqcPtbOyLbtvPnwtTl9vqI5fTKzmJl9ERgIrE3T7BRgVad9K+L70z1uhZlVJTaCSYlERIrCR/es7hL4CX2sndF76nL6fKGHvpmdYGb7gFbgfuA8d9+YpnkN0NBpX0N8fzoLgcakTZOtiUhR8PZ2Yp55aHgfP4i3t+fsOUMPfaAemAx8GvgJ8JCZTcjh4y8CqpO2kTl87KKyc+dOrr76asaOHUu/fv0YPnw406ZN4yc/+UnHENQxY8ZgZpgZ/fv3Z8yYMXzhC1/g+eefD7l6keixsjLaLJaxTZvFcjqaJ/TQd/f97r7Z3V9194XABoK++1R2AsM77Rse35/u8VvdvSmxUeB59LuZ2jpn/vjHPzJlyhSeffZZbr/9dl577TXWrl3Lddddx1NPPcWqVR/2it1yyy3s2LGD+vp6li1bxuDBg5k9ezbf/e53C1KriHxo69CZHPTUUXzQy9g6dGZOn68Yx+mXARVp7lsLzALuStpXS/prAKHY13qQ76+oZ9WmBg60OeUxY/b44XxzzjgqK/Lzks+fP58+ffrwyiuvMHDgwI79xx13HOecc84hv3wGDRpETU3QIzZ69GimT5/OiBEjuPHGG/n85z/PuHHj8lKjiHQ16cvfY/tdLzOybTt97MNunINexrbYSCZ9+Xs5fb5Qz/TNbJGZTTezMfG+/UXATOA/4vcvi+9L+BFwppktMLOPm9nNBEM97yl07ensaz3I+fe9xLK177D9vQ9oaGph+3sfsGztO5x/30vsaz2Y8+f861//yrPPPssVV1xxSOAn6+5DZldffTXuzhNPPJHz+kQkvcqqIQy95gVeGXYBf7Fh7OIo/mLDeGXYBRx9zQs5H6cf9pn+MGAZMILgIuvvgDnuvjJ+/2ig41efu68xs38GbgNuJ/hw1rnu/mZBq87g+yvq2bxrH+2denXaHTbv2scPVtRz09kpP4Zw2DZv3oy7dzlDHzp0KC0tLQBcccUVLF68OO1jHHXUUQwbNox33nknp7WJSPcqq4Zw8hUPAh9+IveYPD1XqKHv7v/Szf0zU+x7DHgsXzUdqVWbGroEfkK7w8pNDTkP/XTWrVtHe3s7F198Ma2trd22d3dNOyESsnxOwQDhn+n3Ku7OgbbMF24PtnnOw3Xs2LGYGfX19YfsP+6444BgGoPu/PWvf2X37t0ce+yxOatLRIpP6KN3ehMzozyWOcz7xCznZ9Mf+chHqK2t5Z577uFvf/vbYT3Gj370I8rKyjj33HNzWpuIFBeFfo7NHj+csjSZXmZQO77ziNPcuO+++zh48CBTp07l0UcfZdOmTdTX1/OLX/yCt956i1jsw7HAzc3N7Ny5k23btvHCCy9w6aWXctttt/Hd736XsWPH5qU+ESkOVqhx5MUiPhVDY2NjY5eF0VtaWtiyZQvHHnvsYc/kmBi90/libpnB2GGVPD5/Wt6Gbe7YsYPbb7+dX/3qV2zfvp2KigomTJjAhRdeyPz58xkwYABjxoxh69atAPTt25eamhpOPvlkLr/8ck4//fS81JULufi/EemtmpqaqK6uBqiOfx4pLYV+klwFy77Wg/xgRT0rNzVwsM3pEzNqxw9nQR7H6fd2Cn2R9LIJfSVQHlRW9OGmsydy09kTNSJGRIqK+vTzTIEvIsVEoS8iEiEKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0BcRiRCFvohIhCj0e4l58+ZhZtxxxx2H7F++fHnePitw8803M3ny5Lw8tojkh0K/F+nXrx+LFy/mvffeC7sUESlSCv18K+DcRrNnz6ampoZFixalbfPiiy9y2mmn0b9/f0aNGsVVV13VMR3zPffcw6RJkzraJv5KuP/++w95juuvv56lS5fyne98hw0bNmAWTBe9dOlSAP70pz9xzjnnUFlZSVVVFV/4whdoaGjoeIzEXwgPP/wwY8aMobq6mi9+8Ys0Nxd0zXqRSFLo50NrMzx9Hdx1Atw5Pvj69HXB/jyKxWLcfvvt/PjHP2b79u1d7v/DH/7AmWeeyQUXXMDvfvc7Hn30UV588UWuvPJKAGbMmMHGjRvZvXs3AHV1dQwdOpTVq1cDcODAAdauXcvMmTO56KKLWLBgARMnTmTHjh3s2LGDiy66iPb2ds455xzeffdd6urqWLlyJX/84x+56KKLutSyfPlynnrqKZ566inq6uq6dE2JSO4p9HOttRkerIV1S2Dvn6B5R/B1/ZJgf56D/7zzzmPy5MncdNNNXe5btGgRF198Mddccw0f+9jHOPXUU7n77rtZtmwZLS0tTJo0iaOOOoq6ujoAVq9ezYIFCzpur1u3jgMHDnDqqafSv39/Kisr6dOnDzU1NdTU1NC/f3+ee+453njjDX75y19y4okn8ulPf5ply5ZRV1fH+vXrO2ppb29n6dKlTJo0idNOO40vf/nLPPfcc3l9bUREoZ97z90Ku+tJWs894O2wpx6evy3vJSxevJiHHnqITZs2HbJ/w4YNLF26lMrKyo5tzpw5tLe3s2XLFsyM6dOns3r1avbu3cvGjRuZP38+ra2tvPXWW9TV1fGpT32KAQMGpH3uTZs2MWrUKEaNGtWxb8KECQwePPiQesaMGcOgQYM6bo8YMYJdu3bl8FWQIxaxadejQlMr59rbz9Al8BO8HeqfhrmL81rC9OnTmTNnDgsXLmTevHkd+/ft28dll13GVVdd1eWY0aNHAzBz5kx++tOf8utf/5opU6ZQVVXV8Yugrq6OGTNm5KTG8vLyQ26bGe3taV43KZzW5uDE5e1noO0AxMrh+Lkw6waoGNT98VL0FPq55B78oGTSdiBol+cpl++44w4mT57MuHHjOvZ98pOfZOPGjRmXRJwxYwbXXHMNjz32GDNnzgSCXwSrVq3ipZdeYsGCBR1t+/btS1tb2yHHjx8/nm3btrFt27aOs/2NGzeyd+9eJkyYkMPvUHIu0TXZ+S/V9UtgywvwtZUK/l5A3Tu5ZBacGWUSK8974AOccMIJXHzxxdx9990d+771rW+xZs0arrzySl5//XV+//vf88QTT3RcyAX4xCc+wZAhQ/jlL395SOgvX76c1tZWpk2b1tF2zJgxbNmyhddff509e/bQ2trK7NmzO577t7/9LevWreMrX/kKM2bMYOrUqXn/vuUIFEHXpOSfQj/Xjp8LluZltTIYd1bBSrnlllsO6TL5xCc+QV1dHW+//TannXYaU6ZM4cYbb+SYY475sEQzTjvtNMyMz3zmMx3HVVVVMXXqVAYOHNjR9oILLuDMM8/k9NNP5+ijj+aRRx7BzHjiiScYMmQI06dPZ/bs2Rx33HE8+uijBfu+5TD1pGtSSp7WyE2Sk3VYE38i76kPflA6nrgMho7Tn8iHSWvk5pl7MLy4eUf6NoNGwL9uKshfqpKdbNbI1Zl+rlUMCoL9pEth8OjgB2Xw6OC2Al+KVRF1TUp+6UJuPlQMCkbozF1ckIu2Ijlx/Nzgoq2n6OIpcNek5I/O9PNNgS+lYtYNQRdk52tSia7JM64Ppy7JKYW+iATUNRkJ6t4RkQ+pa7LX05l+ClEb0VQK9H8SAgV+r6TQT5KYGuD9998PuRLpLPF/0nn6BhHJjrp3ksRiMQYPHtwx8deAAQPytuqU9Iy78/7777Nr1y4GDx5MLBYLuySRkqbQ76SmpgZAMz4WmcGDB3f834jI4VPod2JmjBgxgmHDhnHgQDeTp0lBlJeX6wxfJEdCDX0zWwicD3wc+ABYA3zL3eszHDMP+Hmn3a3untPP5sdiMQWNiPQ6YV/InQHcC5wM1ALlwLNmNjDjUdAEjEjaPprPIkVEeotQz/Td/czk2/Gz+F3AicALmQ/1nXksTUSkVwr7TL+z6vjXd7tpV2lmW81sm5k9YWYT0zU0swozq0psgD5WKCKRVTShb2ZlwF3AS+7+Zoam9cBXgXOALxF8D2vMbGSa9guBxqRte65qFhEpNUUzn76Z/QSYC3zG3XsczGZWDmwCHnH3G1LcXwFUJO0aBGxPNZ++iEgpymY+/aIYsmlm9wCfA6ZnE/gA7n7AzF4DUi786u6tQGvScx1JqSIiJS3U7h0L3AOcB5zh7lsO4zFiwAlAhiV/REQEwj/Tvxf4Z4L++WYzS3zkstHdPwAws2XAn919Yfz2jcBvgM3AYOBagiGbDxa2dBGR0hN26H89/nV1p/3/A1ga//doDl2teQiwBKgB3gNeBU519415q1JEpJcomgu5hZJpYXQRkVKkhdFFciliJ0bSu4XdvSNSnFqb4blb4e1noO0AxMqDhcNn3aBlA6WkKfRFOmtthgdrYXc9h1xOWr8Etryg9WKlpKl7R6Sz527tGvgA3g576uH520IpSyQXFPoinb39DF0CP8Hbof7pgpYjkksKfZFk7kEffiZtB3RxV0qWQl8kmVlw0TaTWHnQTqQEKfRFOjt+LliaHw0rg3FnFbYekRxS6It0NusGGDqua/BbWbD/jOvDqUskBxT6Ip1VDAqGZZ50KQweDYNGBF9PulTDNaXkaRoGke64qw9fipqmYRDJJQW+9CIKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0BcRiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IiIRotCX8ERshleRYtAn7AIkYlqb4blbg8XH2w4ESw8ePzdYuETz1IvknUJfCqe1GR6shd31QPuH+9cvgS0vaIESkQJQ944UznO3dg18AG+HPfXw/G2hlCUSJQp9KZy3n6FL4Cd4O9Q/XdByRKJIoS+F4R704WfSdkAXd0XyTKEvhWEWXLTNJFaupQlF8kyhL4Vz/FywNG85K4NxZxW2HpEIUuhL4cy6AYaO6xr8VhbsP+P6cOoSiRCFvhROxaBgWOZJl8Lg0TBoRPD1pEs1XFOkQMwjduHMzKqAxsbGRqqqqsIuJ9rc1YcvkgNNTU1UV1cDVLt7U6a2oZ7pm9lCM1tvZs1mtsvMlpvZuB4cd6GZvWVmLWb2hpmpM7gUKfBFCi7s7p0ZwL3AyUAtUA48a2YD0x1gZqcCjwA/A6YAy4HlZjYp79WKiJS4oureMbOjgV3ADHd/IU2bR4GB7v65pH2/AV5398t78Bzq3hGRXqVkundSqI5/fTdDm1OAVZ32rYjv78LMKsysKrEBulooIpFVNKFvZmXAXcBL7v5mhqY1QEOnfQ3x/aksBBqTtu1HVqmISOkqmtAn6NufBHwxx4+7iOAviMQ2MsePLyJSMopiamUzuwf4HDDd3bs7E98JDO+0b3h8fxfu3gq0Jj3XEVQqIpI77l7wTAo19C34bn8MnAfMdPctPThsLTCLoCsooTa+X0SkqO1rPcj3V9SzalMDB9qc8pgxe/xwvjlnHJUV+Y/ksM/07wX+GTgHaDazRL98o7t/AGBmy4A/u/vC+H0/AurMbAHwK4LuoKnApQWtXEQkS/taD3L+fS+xedc+2pMGTi5b+w5r/rCHx+dPy3vwh92n/3WCfvbVwI6k7aKkNqOBEYkb7r6G4BfFpcAG4PPAud1c/BURCd33V9R3CXyAdofNu/bxgxX1ea8h1DN9d++2M8vdZ6bY9xjwWD5qEhHJl1WbGroEfkK7w8pNDdx09sS81hD2mb4UkyL6oJ5Ib+PuHGjL/DN2sM3J9wdmw+7Tl7C1Ngdr1779TLByVaw8mPd+1g2a9VIkh8yM8ljmzo0+Mcv7aB6d6UdZazM8WAvrlsDeP0HzjuDr+iXB/tbmsCsU6VVmjx9OWZpMLzOoHd95NHruKfSj7LlbYXc9XRYr93bYUw/P3xZKWSK91TfnjGPssMouwV9mMHZYJQvmdDvJ8BFT6EfZ28/QJfATvB3qny5oOSK9XWVFHx6fP41LThnDyCH9qanqx8gh/bnklDEFGa4J6tOPLvegDz+TtgNa6EQkxyor+nDT2RO56eyJoXwiV2f6UWUWXLTNJFauwBfJozCmhVHoR9nxc7suUp5gZTBOC5KJ9DYK/SibdQMMHdc1+K0s2H/G9eHUJSJ5o9CPsopB8LWVcNKlMHg0DBoRfD3p0mC/xumL9DpFtVxiIWi5xAx00VakJJXycokSJgW+SK+n0BcRiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6pS5iH64TkSOjqZVLkZY4FJHDpNAvNYklDjuveLV+CWx5QXPmiBShMObNT0ehX2p6ssTh3MWhlCYiH9rXepDvr6hn1aYGDrQ55TFj9vjhfHPOuIKskJWO+vRLjZY4FCl6+1oPcv59L7Fs7Ttsf+8DGppa2P7eByxb+w7n3/cS+1oPhlabQr+UZLPEoYiE5vsr6tm8ax/tnX4U2x0279rHD1bUh1MYCv3SoiUORUrCqk0NXQI/od1h5aaGwhaURKFfarTEoUhRc3cOtGX+a/tgmxPWWiYK/VKjJQ5FipqZUR7L/Nd2n5iFNppHoV9qtMShSNGbPX44ZWkyvcygdvzwwhaURMslljotcShSdBKjdzpfzC0zGDusksfnT8vpsM1slktU6IuI5MG+1oP8YEU9Kzc1cLDN6RMzascPZ0Eexukr9DNQ6ItIoeX7E7laGF1EpIgUyxQMoNAXEYmUnIW+mY03sz9mecx0M3vSzP5iZm5m53bTfma8Xeet5oiKFxGJiFye6fcFPprlMQOBDcAVWR43DhiRtO3K8ngRkUjq8SVkM7uzmyZHZ/vk7v4M8Ez88bM5dJe77832+UqChmCKSB5lM27oauB1IN2V4cojrqbnXjezCuBN4GZ3f6mAz517WhRFpNcoprnzU8km9DcDP3T3X6S608wmA6/moqgMdgCXA68AFcDXgNVm9ml3/22auiribROKK0W1KIpIySvWufNTyaZP/xXgxAz3O5DXX2/uXu/uD7j7q+6+xt2/CqwBvpHhsIVAY9K2PZ81Zq0ni6KISNEq5rnzU8km9BcAd6W70903uHsYQ0DXAWMz3L8IqE7aRhaiqB7ToigiJa2Y585Ppcch7e473X2rmZ2ero2ZXZabsrIymaDbJyV3b3X3psQGNBessu5oURSRklfMc+encjhn5v9lZt8zs47VPMxsqJk9CdyRzQOZWaWZTY5fDwA4Nn57dPz+RWa2LKn9NWZ2jpmNNbNJZnYXcAZw72F8H+HToigiJa3Y585P5XBC/3TgPGC9mU0ws88SjKKpIjjrzsZU4LX4BnBn/N+3xG+PAEYnte8L/AB4A6gD/gGY7e7PZf9tFAktiiJSsop97vxUsg59d19DEO5vAr8F/i/wQ2Cmu2/N8rFWu7ul2ObF75/n7jOT2v+bu4919/7u/hF3P93d/zvb76GoaFEUkZJWzHPnp3K4F16PJzhL3w4cJPiE7IBcFRUpWhRFpKR9c844xg6r7BL8ibnzF8wZF05haWQ9tbKZfRv4DvBT4FqCkTMPE3TvfMnd1+a6yFwq+qmV9YlckZJTyLnzU8nrfPpmtgP4anwKhcS+cuB24Cp3r0h7cBEo+tAXkZIWxidyswn9w/kVdIK770ne4e4HgGvN7KnDeDwRkV6jmC7apnI4F3L3ZLiv7sjKEREpDcU0DDMbxTUphIhIESulOXbSKY0qexNdqBUpSYk5djpPubBs7Tus+cMeHp8/rSSCX8slFkJrMzx9Hdx1Atw5Pvj69HXBfhEpCaU2x046Cv18S0ydvG4J7P0TNO8Ivq5fEuxX8IuUhFKbYycdhX6+aepkkZJXinPspKPQzzdNnSxS8kpxjp10FPr5pKmTRXqNUptjJx2Ffj5p6mSRXsHdS26OnXSKf3xRqTt+bnDR1lN08WjqZJGilWpM/ozjj+akMUex+u3docyxkwtZz71T6go+905i9M6e+kODPzF1smbSFCk66cbkJ87qH58/jYF9Y0XTh5/N3Dvq3sk3TZ0sUnJ6Mia/WAI/WzrTLzR9Ilek6H1m8fNsf++DtPePHNKfF791RgErykxn+sVMgS9S1HrTmPxUFPoiIkl605j8VBT6+VKiZwEiUefuvWZMfiqlMcaoVLQ2B9MuvP1M8KGrWHkwZHPWDbpgK1LEOg/PjJUZVf360PjBQZJP30ptTH4qCv1cSQzN7DzPzvolsOUFjdQRKVKZhmdW9y9nQEWM9nZKckx+KqVbebHpycRqcxeHUpqIpJdpeGZTywHOm/J33PhPE0q2D78z9enniiZWEylJPZkyubcEPij0c0MTq4mUpN4+PDMVhX4uaGI1kZLVm4dnpqI+/VzRxGoiJaHzSJ2/tab/K73Uh2emotDPlVk3BKN00k2sdsb14dUmIkD6kTqp9Ibhmako9HMlMbHa87cFF20T4/THnRUEvoZrioQu3UidhIF9YwzqV95rhmem0ru+m7BVDAqGZc5drInVRIpQppE6AEMG9uXX153eq/rwO9OF3HzpxW8akVLU3t7eo5E6vZ3O9HNJZ/ciRaXzRds9+1oztu9tI3VSUegfKc23I1KUsrloC71zpE4qCv0jofl2RIpWdxdtk/XWkTqpqE//SPRkvh0RKTh37/aibazMqKnqx8gh/bnklDE8Pn9arxupk0qo36GZTQeuBU4ERgDnufvybo6ZCdwJTAS2Abe5+9J81plWT+bb0SRrIgVxaP99O7ubM/ffH11ZwZpvn05ZWbTOfcP+tTYQ2AD8O/B4d43N7FjgV8D9wMXALOBBM9vh7ivyWWgX2cy308svDImELdv+ewgu2kYt8CHk0Hf3Z4BngJ5eMb8c2OLuC+K3N5nZZ4BvAIUNfc23I1I0sum/h+hctE2l1H7NnQKs6rRvRXx/SmZWYWZViQ3I3ZXV4+cG0yykfGLNtyNSKCs37swq8KNy0TaVsLt3slUDNHTa1wBUmVl/d/8gxTELgZvyUo3m2xEJTUcf/sYGdjS2ZGxbZkEffnmfsl47vUJPReG7XkRw4TdhELA9J4+s+XZEQpFtH/4xg/v3+ukVeqrUQn8n0LkjbjjQlOYsH3dvBTou4+f8P71vpebbESmwbMfg144frsCPK7XQXwt07iivje8vHH0KVyRU3Y3BT4h6/30qYY/TrwTGJu061swmA++6+5/MbBHwd+7+lfj99wNXmtm/EQzzPAP4AvDZghWtT+GKhKonSxyWGYyo7sc/TqiJdP99KmG/ElOB/066neh7fwiYR/CBrdGJO919i5l9FvghcDVB3/zXCjpGvyefwtUHskTyxsy6XeLwmOr+vPjtMwpUUWkJdcimu692d0uxzYvfP8/dZ6Y4Zoq7V7j73xf807g9+RSuiOTV7PHDKUuT+2UGtROiOQa/J0ptnH64svkUrojkzTfnjGPssMouwa8+/O6F3b1TWvQpXJGiUFnRh8fnT+MHK+pZuamBg23eq5c4zCW9Mtk6fm5w0dZTdPHoU7giBVNZ0Yebzp7ITWdPxN01JLOH1L2TrVk3BJ+27Tz9gj6FKxIaBX7PKfSzlfgU7kmXwuDRMGhE8PWkSzVcU0SKnnnELjrGJ11rbGxspKqq6sgfUJ/CFZGQNTU1UV1dDVDt7k2Z2upM/0gp8EWkhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IiIRotAXEYkQhb6ISIQo9EVEIkShLyISIQp9EZEIUeiLiESIQl9EJEIU+iIiEaLQFxGJEIW+iEiEKPRFRCJEoS8iEiEKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0BcRiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIREhRhL6ZXWFm75hZi5m9bGYnZWg7z8y809ZSyHpFREpV6KFvZhcBdwLfAT4JbABWmNmwDIc1ASOSto/mu04Rkd4g9NAH/hVY4u4/d/eNwOXA+8BXMxzj7r4zaWsoSKUiIiUu1NA3s77AicCqxD53b4/fPiXDoZVmttXMtpnZE2Y2McNzVJhZVWIDBuWqfhGRUhP2mf5QIAZ0PlNvAGrSHFNP8FfAOcCXCL6HNWY2Mk37hUBj0rb9CGsWESlZYYd+1tx9rbsvc/fX3b0OOB/YDVyW5pBFQHXSlu6Xg4hIr9cn5OffA7QBwzvtHw7s7MkDuPsBM3sNGJvm/lagNXHbzA6vUhGRXiDUM3133w+8CsxK7DOzsvjttT15DDOLAScAO/JRo4hIbxL2mT4EwzUfMrNXgHXANcBA4OcAZrYM+LO7L4zfvhH4DbAZGAxcSzBk88FCFy4iUmpCD313f9TMjgZuIbh4+zpwZtIwzNFAe9IhQ4Al8bbvEfylcGp8uKeIiGRg7h52DQUVH7bZ2NjYSFVVVdjliIgcsaamJqqrqwGq3b0pU9uSG70jIiKHT6EvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IiIRotAXEYkQhb6ISIQo9EVEIkShLyISIQp9EZEIUeiLiESIQl9EJEIU+iIiEaLQFxGJEIW+iEiEKPRFRCJEoS8iEiEKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0BcRiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IiIRotAXEYkQhb6ISIQo9EVEIkShLyISIQp9EZEIKYrQN7MrzOwdM2sxs5fN7KRu2l9oZm/F279hZmflo64nn3ySlpaWlPe1tLTw5JNP5uNpRUTyJvTQN7OLgDuB7wCfBDYAK8xsWJr2pwKPAD8DpgDLgeVmNinXtdXW1vLAAw90Cf6WlhYeeOABamtrc/2UIiJ5Ze4ebgFmLwPr3f3K+O0yYBvwY3e/I0X7R4GB7v65pH2/AV5398t78HxVQGNjYyNVVVXd1pcI+Msuu4x+/fp1uS0iErampiaqq6sBqt29KVPbUM/0zawvcCKwKrHP3dvjt09Jc9gpye3jVqRrb2YVZlaV2IBB2dTYr18/LrvsMh544AH27t2rwBeRktYn5OcfCsSAhk77G4CPpzmmJk37mjTtFwI3HW6BEAT/JZdcwoknnsirr76qwBeRkhV6n34BLAKqk7aR2T5AS0sLDz30EK+++ioPPfRQ2ou7IiLFLuzQ3wO0AcM77R8O7ExzzM5s2rt7q7s3JTagOZsCk/vwBw8e3NHVo+AXkVIUaui7+37gVWBWYl/8Qu4sYG2aw9Ymt4+rzdD+sKW6aJvcx6/gF5FSE/aZPgTDNf+nmV1iZuOBnwADgZ8DmNkyM1uU1P5HwJlmtsDMPm5mNwNTgXtyXdjKlStTXrRNBP/KlStz/ZQiInkV+pBNADO7EriW4GLs68BV7v5y/L7VwDvuPi+p/YXAbcAY4PfAde7+dA+fK6shmyIixS6bIZtFEfqFpNAXkd6mZMbpi4hIYSn0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQsKecC00TU0ZRzWJiJSMbPIsiuP0/w7YHnYdIiJ5MNLd/5ypQRRD34BjyHLiNYJ5+LcTzNKZ7bFhU+3hUO3hiGrtg4C/eDehHrnunfgLkvE3YSrB7woAmrv7xFuxUe3hUO3hiHDtPWqvC7kiIhGi0BcRiRCFfs+1At+Jfy01qj0cqj0cqj2DyF3IFRGJMp3pi4hEiEJfRCRCFPoiIhGi0BcRiRCFfhIzu8LM3jGzFjN72cxO6qb9hWb2Vrz9G2Z2VqFqTVFLj2s3s3lm5p22lkLWm1TLdDN70sz+Eq/j3B4cM9PMfmtmrWa22czm5b/SlHVkVXu87s6vu5tZTYFKTtSx0MzWm1mzme0ys+VmNq4Hx4X+fj+c2ovl/W5mXzez35lZU3xba2Zzuzkm56+5Qj/OzC4C7iQYLvVJYAOwwsyGpWl/KvAI8DNgCrAcWG5mkwpS8KG1ZFV7XBMwImn7aL7rTGMgQb1X9KSxmR0L/Ar4b2AycBfwoJnNyVN9mWRVe5JxHPra78pxXd2ZAdwLnAzUAuXAs2Y2MN0BRfR+z7r2uGJ4v28Hvg2cCEwFngeeMLOJqRrn7TV3d23BsNWXgXuSbpcRTNfw7TTtHwWe6rTvN8D9JVD7PGBv2K95irocOLebNouBNzvt+0/gv0qg9pnxdoPDfq071XV0vK7pGdoUzfv9MGovyvd7vLZ3gX8p5GuuM33AzPoS/PZdldjn7u3x26ekOeyU5PZxKzK0z4vDrB2g0sy2mtk2M0t7tlGEiuJ1P0Kvm9kOM1tpZtPCLgaojn99N0ObYn3de1I7FNn73cxiZvZFgr8W16ZplpfXXKEfGArEgIZO+xuAdP2tNVm2z5fDqb0e+CpwDvAlgvfBGjMbma8icyjd615lZv1DqCcbO4DLgQvi2zZgtZl9MqyCzKyMoIvsJXd/M0PTYnm/d8ii9qJ5v5vZCWa2j+ATt/cD57n7xjTN8/KaR26WTQF3X0vS2YWZrQE2AZcBN4RVV2/n7vUEAZSwxsz+HvgG8OVwquJeYBLwmZCe/0j0qPYie7/XE1yLqgY+DzxkZjMyBH/O6Uw/sAdoA4Z32j8c2JnmmJ1Zts+Xw6n9EO5+AHgNGJvb0vIi3eve5O4fhFDPkVpHSK+7md0DfA443d27W1ioWN7vQNa1HyLM97u773f3ze7+qrsvJBgIcHWa5nl5zRX6BP8RwKvArMS++J+Os0jf37Y2uX1cbYb2eXGYtR/CzGLACQTdD8WuKF73HJpMgV93C9wDnAec4e5benBYUbzuh1l758copvd7GVCR5r78vOZhX70ulg24CGgBLgHGAw8A7wHD4/cvAxYltT8VOAAsAD4O3AzsByaVQO03Av8IHEcwxPMR4ANgQgi1VxIE32SCURjfiP97dPz+RcCypPbHAn8D/i3+us8HDgJzSqD2awj6lccSdEvcRfBX2qwC130fsJdg+GNN0tY/qU1Rvt8Ps/aieL/H3w/TgTEEv3QWAe1AbSFf84L+kBT7BlwJbCW4yPIy8Omk+1YDSzu1v5Cgj64VeBM4qxRqB36Y1HYnwbj3KSHVPTMemJ23pfH7lwKrUxzzWrz+PwDzSqF24Dpgczxw/krwWYPTQ6g7Vc2e/DoW6/v9cGovlvc7wXj7d+J17CIYmVNb6NdcUyuLiESI+vRFRCJEoS8iEiEKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0BcRiRCFvohIhCj0RfKsWJZ3FAGFvkheFdnyjiKahkHkSJjZ0cAbwN3ufnt836kE86jMJZjo67PuPinpmP8kWDLxzMJXLFGnM32RI+DuuwlWZbrZzKaa2SDgYYI1i5+jeJcZlIjSylkiR8jdnzazJcB/AK8QTP28MH53xuUdvTQXfpESpjN9kdz4JsFJ1IXAxe7eGnI9Iikp9EVy4++BYwh+psYk7e9tyztKiVP3jsgRMrO+wC+ARwkWvHjQzE5w910ES9ud1emQUl7eUUqcRu+IHCEz+x7weeAfgH1AHdDo7p+LD9l8E7gX+HfgDOBughE9K0IqWSJMoS9yBMxsJrCSYNnDF+P7xgAbgG+7+0/ibX4ITAC2A7e6+9LCVyui0BcRiRRdyBURiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IiIRotAXEYkQhb6ISIQo9EVEIkShLyISIf8felRNogzktxsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
    "historyNt = onp.asarray(historyNt)\n",
    "updatesNt = onp.asarray(updatesNt) \n",
    "axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='#1F77B4', label='GD')\n",
    "axes.scatter(historyNt[:,0], historyNt[:,1], lw=0.5, color='#FF7F0E', label='Newton')\n",
    "axes.scatter([0], [0], lw=0.25, color='black', marker='x') # target at 0,0\n",
    "axes.set_xlabel('x0'); axes.set_ylabel('x1'); axes.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not completely surprising: for this simple example we can reliably evaluate the Hessian, and Newtons's method profits from the second order information. It's trajectory is much more diagonal (that would be the ideal, shortest path to the solution), and does not slow down as much as GD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inverse simulators\n",
    "\n",
    "Now we also use an analytical inverse of $\\mathbf y$ for the optimization. It represents our inverse simulator $\\mathcal P^{-1}$ from the previous sections:\n",
    "$\\mathbf{y}^{-1}(\\mathbf{x}) = [x_0 \\ x_1^{1/2}]^T$, to compute the scale-invariant update denoted by PG below. As a slight look-ahead to the next section, we'll use a Newton's step for $L$, and combine it with the inverse physics function to get an overall update. This gives an update step:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\Delta \\mathbf{x}_{\\text{PG}} &= \n",
    "\\mathbf{y}^{-1} \\left( \\mathbf{y}(\\mathbf{x}) - \\eta\n",
    "  \\left( \\frac{\\partial^2 L }{ \\partial \\mathbf{y}^2 }  \\right)^{-1}\n",
    "  \\frac{\\partial L }{ \\partial \\mathbf{y} }\n",
    "\\right) - \\mathbf{x}\n",
    "\\end{aligned}$$\n",
    "\n",
    "Below, we define our inverse function `physics_y_inv_analytic`, and then evaluate an optimization with the PG update for ten steps:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PG iter 0: [2.1       2.5099802]\n",
      "PG iter 1: [1.4699999 2.1000001]\n",
      "PG iter 2: [1.0289999 1.7569861]\n",
      "PG iter 3: [0.72029996 1.47      ]\n",
      "PG iter 4: [0.50421   1.2298902]\n",
      "PG iter 5: [0.352947 1.029   ]\n",
      "PG iter 6: [0.24706289 0.86092323]\n",
      "PG iter 7: [0.17294402 0.7203    ]\n",
      "PG iter 8: [0.12106082 0.60264623]\n",
      "PG iter 9: [0.08474258 0.50421   ]\n"
     ]
    }
   ],
   "source": [
    "x = np.asarray([3.,3.])\n",
    "eta = 0.3\n",
    "historyPG = [x]; historyPGy = []; updatesPG = []\n",
    "\n",
    "def physics_y_inv(y):\n",
    "    return np.array( [y[0], np.power(y[1],0.5)] )\n",
    "\n",
    "Gy = jax.grad(loss_y)\n",
    "Hy = jax.jacobian(jax.jacobian(loss_y))\n",
    "for i in range(10):\n",
    "    \n",
    "    # Newton step for L(y)\n",
    "    zForw = physics_y(x)\n",
    "    g = Gy(zForw)\n",
    "    h = Hy(zForw)\n",
    "    hinv = np.linalg.inv(h)\n",
    "    \n",
    "    # step in y space\n",
    "    zBack = zForw -eta * np.matmul( hinv , g)\n",
    "    historyPGy.append(zBack)\n",
    "\n",
    "    # \"inverse physics\" step via y-inverse\n",
    "    x = physics_y_inv(zBack)\n",
    "    historyPG.append(x)\n",
    "    updatesPG.append( historyPG[-2] - historyPG[-1] )\n",
    "    print( \"PG iter %d: \"%i + format(x) )\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we obtain $[0.084 \\  0.504]$ as the final position, with a distance of only 0.51! This is clearly better than both Newton and GD.\n",
    "\n",
    "Let's directly visualize how the PGs (in red) fare in comparison to Newton's method (orange) and GD (blue)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc742ddba60>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxeUlEQVR4nO3de5RU5Z3v//e3+salaSBj21yUaZFIUMzB6BgvEaGBgzg5XpIYmTCTsHLmp4n6s5lDYkKMlxgjcXKDxFwckzNKosIvazkQjQ6CLWRUjMZEkwyklQBGlKtBoFGa7q7v749dhUV3VXUXVNWu6v15rVWr6V3Prvp2Uf3pXc9+9vOYuyMiItEQC7sAEREpHoW+iEiEKPRFRCJEoS8iEiEKfRGRCFHoi4hEiEJfRCRCFPoiIhFSGXYBxWZmBowC9oddi4hIHg0B3vBerriNXOgTBP7WsIsQESmAE4DXszWIYujvB3jttdeoq6sLuxYRkWO2b98+TjzxROhDD0YUQx+Auro6hb6IRI5O5IqIRIhCX0QkQhT6IiIREtk+/Wzcnc7OTrq6usIuRYCKigoqKysJRtuKyLFQ6Hdz6NAhtm3bxttvvx12KZJi0KBBjBw5kurq6rBLESlrCv0U8XiczZs3U1FRwahRo6iurtbRZcjcnUOHDrFr1y42b97Me9/7XmIx9UpK/xXv6iJWUVGwxw819M3ss8BngcbEpv8GbnP3x7LscwXw1cQ+rwBfcPdH81HPoUOHiMfjnHjiiQwaNCgfDyl5MHDgQKqqqnj11Vc5dOgQAwYMCLskkbzat20rm5o/yYA/vwFxIAYHTx7F2MVLqBt5Ql6fK+xDpq3AF4EzgbOAFmCFmZ2WrrGZnQc8CPwEOANYDiw3s4n5LEpHkqVH/yfSX+3btpU3LptBze/fwA8Y/o7hB4ya37/BG5fNYN+2/E4gEOpvkrs/7O6Puvsr7v6yu98ItAHnZNilGfhPd/+Gu29w95uA3wLXFatmEZF82tT8SXyvA927kg3f62xq/mRen69kDp/MrMLMZgODgXUZmp0LrO62bWVie6bHrTGzuuSNYFIiEZGSMODPb9Az8JMscX/+hB76Zna6mbUB7cCPgMvdfX2G5iOAHd227Uhsz2QBsDflpsnWRKQkxLu6gj78rI0S7fIk9NAHWoFJwAeBHwL3mdmpeXz8hcDQlFt+z4qUkO3bt9Pc3My4ceMYMGAADQ0NnH/++fzwhz88PAS1sbERM8PMGDhwII2NjXz84x+npaUl5OpFoidWUdF7CsfI62ie0EPf3Q+5+0Z3f8HdFwAvEfTdp7MdaOi2rSGxPdPjt7v7vuSNIs+j38vU1nmzadMmzjjjDB5//HHuuOMOfve737Fu3TpuuOEGHnnkEVavfrdX7LbbbmPbtm20trayZMkShg0bxvTp0/na175WlFpF5F0HTx4FZMoJT9yfP6U4Tj8G1GS4bx0wDViUsm0Gmc8BhKKtvZNvrmxl9YYddHQ5VRXG9AkNfG7meGprCvOSX3PNNVRWVvKb3/yGwYMHH94+duxYLr300iP++AwZMoQRI4IesTFjxjB58mRGjhzJzTffzMc+9jHGjx9fkBpFpKexi5fwxmUz0pzMdWyoMXbxkrw+X6hH+ma20Mwmm1ljom9/ITAFuD9x/5LEtqTFwEVmNt/M3mdmtxIM9byr2LVn0tbeyUd+8DRL1m1h65532LHvIFv3vMOSdVv4yA+epq29M+/P+eabb/L4449z7bXXHhH4qXq7yKy5uRl3Z8WKFXmvT0Qyqxt5AqOWr6L9/aOwwY4NdGyw0/7+UYxavirv4/TDPtI/HlgCjCQ4yfp7YKa7r0rcP4aU0xzu/oyZfQK4HbiD4OKsy9z9j0WtOotvrmxl48424t0+rcUdNu5s41srW7nlkrSXIRy1jRs34u49jtCPO+44Dh48CMC1117LnXfemfEx3vOe93D88cezZcuWvNYmIr2rG3kCk/6/4Lxav74i193/dy/3T0mz7efAzwtV07FavWFHj8BPijus2rAj76GfyXPPPUc8HmfOnDm0t7f32t7dNe2ESMgKGfgQ/pF+v+LudHRlP3Hb2eV5D9dx48ZhZrS2th6xfezYsUAwjUFv3nzzTXbt2sVJJ52Ut7pEpPSEPnqnPzEzqiqyh3llheX9aPpv/uZvmDFjBnfddRcHDhw4qsdYvHgxsViMyy67LK+1iUhpUejn2fQJDcQyZHrMYMaE7iNO8+MHP/gBnZ2dnHXWWSxbtowNGzbQ2trKz372M/70pz9RkfKRcf/+/Wzfvp3XXnuNX/3qV1x11VXcfvvtfO1rX2PcuHEFqU9ESoMVaxx5qUhMxbB37969PRZGP3jwIJs3b+akk0466pkck6N3up/MjRmMO76Wh645v2DDNrdt28Ydd9zBL3/5S7Zu3UpNTQ2nnnoqV1xxBddccw2DBg2isbGRV199FYDq6mpGjBjBOeecw2c+8xmmTp1akLryIR//NyL91b59+xg6dCjA0MT1SBkp9FPkK1ja2jv51spWVm3YQWeXU1lhzJjQwPwCjtPv7xT6IpnlEvpKoAKoranklktO45ZLTtOIGBEpKerTLzAFvoiUEoW+iEiEKPRFRCJEoS8iEiEKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0O8n5s6di5nx9a9//Yjty5cvL9i1ArfeeiuTJk0qyGOLSGEo9PuRAQMGcOedd7Jnz56wSxGREqXQL7Qizm00ffp0RowYwcKFCzO2eeqpp7jgggsYOHAgJ554Itdff/3h6ZjvuusuJk6ceLht8lPCj370oyOe48tf/jL33nsvX/nKV3jppZcwC6aLvvfeewH4y1/+wqWXXkptbS11dXV8/OMfZ8eOHYcfI/kJ4ac//SmNjY0MHTqU2bNns39/UdesF4kkhX4htO+HR2+ARafDtycEXx+9IdheQBUVFdxxxx1873vfY+vWrT3u//Of/8xFF13ERz/6UX7/+9+zbNkynnrqKa677joALrzwQtavX8+uXbsAWLt2Lccddxxr1qwBoKOjg3Xr1jFlyhSuvPJK5s+fz2mnnca2bdvYtm0bV155JfF4nEsvvZS//vWvrF27llWrVrFp0yauvPLKHrUsX76cRx55hEceeYS1a9f26JoSkfxT6Odb+3748Qx47h546y+wf1vw9fl7gu0FDv7LL7+cSZMmccstt/S4b+HChcyZM4d58+bx3ve+l/POO4/vfve7LFmyhIMHDzJx4kTe8573sHbtWgDWrFnD/PnzD3//3HPP0dHRwXnnncfAgQOpra2lsrKSESNGMGLECAYOHMgTTzzBH/7wBx544AHOPPNMPvjBD7JkyRLWrl3L888/f7iWeDzOvffey8SJE7ngggv4p3/6J5544omCvjYiotDPvye+CrtaSVnPPeBx2N0KLbcXvIQ777yT++67jw0bNhyx/aWXXuLee++ltrb28G3mzJnE43E2b96MmTF58mTWrFnDW2+9xfr167nmmmtob2/nT3/6E2vXruXv/u7vGDRoUMbn3rBhAyeeeCInnnji4W2nnnoqw4YNO6KexsZGhgwZcvj7kSNHsnPnzjy+CnLMIjbtelRoauV8e/kxegR+kseh9VGYdWdBS5g8eTIzZ85kwYIFzJ079/D2trY2rr76aq6//voe+4wZMwaAKVOm8G//9m/813/9F2eccQZ1dXWH/xCsXbuWCy+8MC81VlVVHfG9mRGPZ3jdpHja9wcHLi8/Bl0dUFEFp8yCaTdBzZDe95eSp9DPJ/fgFyWbro6gXYGnXP7617/OpEmTGD9+/OFtH/jAB1i/fn3WJREvvPBC5s2bx89//nOmTJkCBH8IVq9ezdNPP838+fMPt62urqarq+uI/SdMmMBrr73Ga6+9dvhof/369bz11luceuqpefwJJe+SXZPdP6k+fw9s/hX88yoFfz+g7p18MguOjLKpqCp44AOcfvrpzJkzh+9+97uHt33hC1/gmWee4brrruPFF1/klVdeYcWKFYdP5AK8//3vZ/jw4TzwwANHhP7y5ctpb2/n/PPPP9y2sbGRzZs38+KLL7J7927a29uZPn364ef+7W9/y3PPPccnP/lJLrzwQs4666yC/9xyDEqga1IKT6Gfb6fMAsvwsloMxl9ctFJuu+22I7pM3v/+97N27VpefvllLrjgAs444wxuvvlmRo0a9W6JZlxwwQWYGR/60IcO71dXV8dZZ53F4MGDD7f96Ec/ykUXXcTUqVOpr6/nwQcfxMxYsWIFw4cPZ/LkyUyfPp2xY8eybNmyov3ccpT60jUpZU9r5KbIyzqsyY/Iu1uDX5TDTxyD48brI/JR0hq5BeYeDC/evy1zmyEj4f9sKMonVclNLmvk6kg/32qGBMF+9lUwbEzwizJsTPC9Al9KVQl1TUph6URuIdQMCUbozLqzKCdtRfLilFnBSVtP08VT5K5JKRwd6ReaAl/KxbSbgi7I7uekkl2TTV8Opy7JK4W+iATUNRkJ6t4RkXepa7Lf05G+iKSnwO+XFPoiIhGi0BcRiRCFvogck6hd4FnudCJXRHLW1XaAXYsX0dbyJN7RgVVVUds0lfrmeVTUDu79ASQ0oR7pm9kCM3vezPab2U4zW25m43vZZ66ZebfbwWLVXKqSC6ObGdXV1YwbN47bbruNzs5OIDgau+eeezj33HOpq6ujtraW0047jebmZjZu3Bhy9VJOutoOsGX2bPbc/wAdr79O586ddLz+Onvuf4Ats2fT1XYg7BIli7C7dy4Evg+cA8wAqoDHzay3Q4V9wMiU298WsshjUcyPvhdddBHbtm3jlVdeYf78+dx666184xvfwN35xCc+wfXXX8/FF1/M448/zvr16/nJT37CgAEDuP12zZ4ofbdr8SIObdoE3dc/iMc5tGkTuxYvDqcw6ZNQu3fc/aLU781sLrATOBP4VfZdfXsBSzsmYX30rampYcSIEQB89rOf5T/+4z/4xS9+wUknncTSpUtZsWIFl1xyyeH2Y8aM4ZxzzlGfrOSkreXJnoGfFI/T1tICN36puEVJn5Van/7QxNe/9tKu1sxeJfik8lvgS+7+3+kamlkNUJOyqaCXFSY/+nY/Etpz/wMcWPcsjUuXFq3Pc+DAgbz55ps8+OCDjB8//ojAT2Uajy195O54R/aFgryzE3fX+6pEhd29c5iZxYBFwNPu/scsTVuBTwOXAv9I8DM8Y2YnZGi/ANibctuar5rTKYWPvu7O6tWrWblyJU1NTbz88stHrKAFMG/evMPr5J5wQqaXTuRIZoZVZZ+N0yorFfglrGRCn6BvfyIwO1sjd1/n7kvc/UV3Xwt8BNgFXJ1hl4UEnyCSt4ImXJ8++hbII488Qm1tLQMGDGDWrFlceeWV3HrrrWnb3njjjbz44ovcfPPNtLW1Fawm6X9qm6ZCLEN0xGLUNjUVtyDJSUl075jZXcCHgcnuntORuLt3mNnvgLQLv7p7O9Ce8lzHUmpvtYT60Xfq1Kn88Ic/pLq6mlGjRlFZGfz3vve976W1tfWItvX19dTX13P88cfnvQ7p3+qb53Fg3bM9P9HGYlSfPJb65ubwipNehT1k0xKBfznQ5O6bj+IxKoDTgSxL/hRH2B99Bw8ezLhx4xgzZszhwAf4h3/4B1pbW1mxYkVBnleipaJ2MI1LlzJ8zhyqRo+msqGBqtGjGT5nDo0PFu+clRydsI/0vw98gqB/fr+ZjUhs3+vu7wCY2RLgdXdfkPj+ZuBZYCMwDPg8wZDNHxe39PRqm6ay5/4H0nfxhPTRd/bs2Tz00EPMnj2bBQsWMHPmTBoaGnj11VdZtmwZFRUVRa9JyltF7WBG3PgluPFLOmlbZsLu0/8sQT/7GoIj9eTtypQ2YwjG4icNB+4BNgCPAnXAee6+vgj19qq+eR7VY8f27PMM8aOvmbFs2TIWLVrEo48+yrRp0xg/fjyf/vSnOfHEE3nqqaeKXpP0Hwr88qKF0VPka/HtYJz+YtpaWvDOTqyyktqmJuqbm/XR9yhpYXSRzHJZGD3s7p1+SR99+xktJiL9iEK/wBT4Zap9PzzxVXj5MejqgIqqYOHwaTdp2UApawp9ke7a98OPZ8CuViDlhPzz98DmX2m9WClrYZ/IFSk9T3y1Z+ADeBx2t0KLJqiT8qXQF+nu5cfoEfhJHofWR4tajkg+KfTTiNqIpnJQtP8T96APP5uujqCdSBlS6KeoSlxN+/bbb4dciXSX/D+p6uWK52NmFpy0zaaiSqN5pGzpRG6KiooKhg0bxs6dOwEYNGiQRt+EzN15++232blzJ8OGDSvO1cOnzApO2nqaLh6LwfiLC1+DSIEo9LtJLkKSDH4pDcOGDTv8f1Nw024KRunsbj0y+C0Gx42Hpi8Xpw6RAtAVuRl0dXXR0cuMmVIcVVVVxZ8fqH1/MEqn9dF3x+mPvzgIfA3XlBKTyxW5Cn2R3uiKXClxuYS+TuSK9EaBL/2IQl9EJEIU+iIiEaLQFxGJEIW+iEiEKPRFRCJEoS8iEiEKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0BcRiRCFvohIhCj0JTwRm+FVpBRoERUprvb98MRXg8XHk/PUnzIrWLikn8xT7+5acU1KlkJfiqd9P/x4BuxqBVJWpHr+nmClqn9eVbbB39V2gF2LF9HW8iTe0YFVVVHbNJX65nlU1A4OuzyRw9S9I8XzxFd7Bj4ESxLubg1WqipDXW0H2DJ7Nnvuf4CO11+nc+dOOl5/nT33P8CW2bPpajsQdokihyn0pXhefowegZ/k8WBpwjK0a/EiDm3aBPFuP1s8zqFNm9i1eHE4hYmkodCX4nAP+vCz6eooy5O7bS1P9gz8pHictpaW4hYkkoVCX4rDLDhpm01FVdktTejueEf2P2be2UnU1qKW0qXQl+I5ZRZYhrecxWD8xcWtJw/MDKvK/sfMKis1mkdKhkJfimfaTXDc+J7Bb7Fge9OXw6nrGNU2TYVYhl+lWIzapqbiFiSShUJfiqdmSDAs8+yrYNgYGDIy+Hr2VWU9XLO+eR7VY8f2DP5YjOqTx1Lf3BxOYSJpWNT6Gs2sDti7d+9e6urqwi4n2tzLrg8/k2Cc/mLaWlrwzk6sspLapibqm5s1Tl8Kbt++fQwdOhRgqLvvy9Y21NA3swXAR4D3Ae8AzwBfcPfWXva7Avgq0Ai8ktinT+P9FPpSaLoiV4otl9APu3vnQuD7wDnADKAKeNzMMh4amdl5wIPAT4AzgOXAcjObWPBqRfpAgS+lrKS6d8ysHtgJXOjuv8rQZhkw2N0/nLLtWeBFd/9MH55DR/oi0q+U05F+d0MTX/+apc25wOpu21YmtvdgZjVmVpe8AeV5tlBEJA9KJvTNLAYsAp529z9maToC2NFt247E9nQWAHtTbluPrVIRkfJVMqFP0Lc/EZid58ddSPAJInk7Ic+PLyJSNkpiamUzuwv4MDDZ3Xs7Et8ONHTb1pDY3oO7twPtKc91DJWKiORPGCO9Qg19C37a7wGXA1PcfXMfdlsHTCPoCkqakdguIlLS2to7+ebKVlZv2EFHl1NVYUyf0MDnZo6ntqbwkRz2kf73gU8AlwL7zSzZL7/X3d8BMLMlwOvuviBx32JgrZnNB35J0B10FnBVUSsXEclRW3snH/nB02zc2UY8ZeDkknVbeObPu3nomvMLHvxh9+l/lqCffQ2wLeV2ZUqbMcDI5Dfu/gzBH4qrgJeAjwGX9XLyV0QkdN9c2doj8AHiDht3tvGtlVmvS82LUI/03b3Xzix3n5Jm28+BnxeiJhGRQlm9YUePwE+KO6zasINbLjmtoDWEfaQvpaSELtQT6W/cnY6u7L9jnV1e8LUXwu7Tl7C17w/Wrn35sWDlqoqqYN77aTeV7ayXIqXIzKiqyN65UVlhBR/NoyP9KGvfDz+eAc/dA2/9BfZvC74+f0+wvX1/2BWK9CvTJzQQy5DpMYMZE7qPRs8/hX6UPfFV2NVKj8XKPQ67W6Hl9lDKEumvPjdzPOOOr+0R/DGDccfXMn/m+ILXoNCPspcfo0fgJ3kcWvs0W7WI9FFtTSUPXXM+nzq3kROGD2RE3QBOGD6QT53bWJThmqA+/ehyD/rws+nq6FcLnYiUgtqaSm655DRuueS0UK7I1ZF+VJkFJ22zqahS4IsUUBjTwij0o+yUWT0XKU+yGIy/uLj1iEjBKfSjbNpNcNz4nsFvsWB705fDqUtECkahH2U1Q+CfV8HZV8GwMTBkZPD17KuC7RqnL9LvlNRyicWg5RKzKIOTtlp0XKSnXJZL1OgdeVeJhmlX2wF2LV5EW8uTeEcHVlVFbdNU6pvnUVE7OOzyRMqKjvSlpHW1HWDL7Nkc2rQJ4inXFMRiVI8dS+PSpQp+ibxyXhhd5Ai7Fi/qGfgA8TiHNm1i1+LF4RQmUqYU+lLS2lqe7Bn4SfE4bS0txS1IpMwp9KVkuTvekf2qYe/sLPhUtCL9iUJfSpaZYVXZrxq2ykqN5hHJgUJfSlpt01SIZXibxmLUNjUVtyCRMqfQl5JW3zyP6rFjewZ/LEb1yWOpb24OpzCRMqXQL3f9vD+7onYwjUuXMnzOHKpGj6ayoYGq0aMZPmcOjQ9quKZIrjROvxxFeIlDXZEr0pOuyO3Pkkscdl/x6vl7YPOv+v2cOQp8KUeldLCi0C83fVnicNadoZQmIu9qa+/kmytbWb1hBx1dTlWFMX1CA5+bOb4oK2Rloj79cqMlDkVKXlt7Jx/5wdMsWbeFrXveYce+g2zd8w5L1m3hIz94mrb2ztBqU+iXk1yWOBSR0HxzZSsbd7YR7/arGHfYuLONb61sDacwFPrlRUscipSF1Rt29Aj8pLjDqg07iltQCoV+udEShyIlzd3p6Mr+abuzy0ObPkShX260xKFISTMzqiqyf9qurLDQRvMo9MuNljgUKXnTJzQQy5DpMYMZExqKW1AKXZxV7spgiUORqEmO3ul+MjdmMO74Wh665vy8DtvM5eIshb6ISAG0tXfyrZWtrNqwg84up7LCmDGhgfkFGKev0M9CoS8ixVboK3K1XKKISAkplSkYQKEvIhIpeQt9M5tgZpty3GeymT1sZm+YmZvZZb20n5Jo1/024piKFxGJiHwe6VcDf5vjPoOBl4Brc9xvPDAy5bYzx/1FRCKpz6eQzezbvTSpz/XJ3f0x4LHE4+ey6053fyvX5ysLZTQEs5SmixWRvsll3FAz8CKQ6cxw7TFX03cvmlkN8EfgVnd/uojPnX9ltChKV9sBdi1eRFvLk3hHB1ZVRW3TVOqb52kVKxFK/2Aol9DfCHzH3X+W7k4zmwS8kI+istgGfAb4DVAD/DOwxsw+6O6/zVBXTaJtUmmlaBktitLVdoAts2dzaNMmiL9b6577H+DAumdpXKrlCyWaSnXu/HRy6dP/DXBmlvsdKOifN3dvdfe73f0Fd3/G3T8NPAP8S5bdFgB7U25bC1ljzvqyKEqJ2LV4UY/AByAe59CmTexavDicwkRCVMpz56eTS+jPBxZlutPdX3L3MIaAPgeMy3L/QmBoyu2EYhTVZ2W0KEpby5M9Az8pHqetpaW4BYmUgFKeOz+dPoe0u29391fNbGqmNmZ2dX7Kyskkgm6ftNy93d33JW/A/qJV1psyWhTF3fGO7LV6Z2do08WKhKWU585P52iOzP/TzL5hZodX8zCz48zsYeDruTyQmdWa2aTE+QCAkxLfj0ncv9DMlqS0n2dml5rZODObaGaLgCbg+0fxc4SvjBZFMTOsKnutVllZ0iewRPKt1OfOT+doQn8qcDnwvJmdamZ/TzCKpo7gqDsXZwG/S9wAvp34922J70cCY1LaVwPfAv4ArAX+BzDd3Z/I/ccoEWW0KEpt01SIZag1FqO2qam4BYmErNTnzk8n59B392cIwv2PwG+B/wC+A0xx91dzfKw17m5pbnMT98919ykp7f/V3ce5+0B3/xt3n+ruT+b6M5SUMloUpb55HtVjx/YM/liM6pPHUt/cHE5hIiEq5bnz0znaE6+nEBylbwU6Ca6QHZSvoiKljBZFqagdTOPSpQyfM4eq0aOpbGigavRohs+ZQ+ODGq4p0fS5meMZd3xtj+BPzp0/f+b4cArLIOeplc3si8BXgH8DPk8wcuanBN07/+ju6/JdZD6V/NTKuiJXpOwUc+78dAo6n76ZbQM+nZhCIbmtCrgDuN7dazLuXAJKPvRFpKyFcTCUS+gfzZ+g0919d+oGd+8APm9mjxzF44mI9Bul/un3aE7k7s5y39pjK0dEpDyU0jDMXJTWpBAiIiWsnObYyaQ8quxPyuhErYi8KznHTvcpF5as28Izf97NQ9ecXxbBr+USi6F9Pzx6Ayw6Hb49Ifj66A3BdhEpC+U2x04mCv1CS06d/Nw98NZfYP+24Ovz9wTbFfwiZaHc5tjJRKFfaGU0dbKIpFeOc+xkotAvtDKaOllE0ivHOXYyUegXUhlNnSwi2ZXbHDuZKPQLqYymThaRzNy97ObYyaT0xxeVu1NmBSdtPU0XT5GnTtZcOSJ9l25M/oWn1HN243tY8/KuUObYyYec594pd0Wfeyc5emd365HBn5w6ucAzaXa1HWDX4kW0tTyJd3RgVVXUNk2lvnmeZsUUySDTmPzkUf1D15zP4OqKkjmIKuiEa+UulAnX2vcHo3RaHw368CuqgiP8pi8XPPC3zJ7dczHzWIzqsWNpXKrpkEXSufUX/82SdVvSDtGMGXzq3EZuueS04heWgUI/i9Bn2SziFbnbv/Y19tz/QPrFzGMxhs+Zw4gbv1SUWkTKyYfubGHrnncy3n/C8IE89YXSWSkul9DXidxiK+LHwbaWJ9MHPkA8TltLS9FqESkX/WlMfjoK/X7K3fGO7MNFvbOzbN+4IoXSn8bkp6PQL5SQw9TMsKrsw0WtsrJs37giheLu/WZMfjrlMcaoXLTvD6ZdePmxd0/YnjIrWPw8hLVua5umZu3Tr20qnT5JkTB1H55ZETPqBlSy951OUg/fym1Mfjo6kZsvyaGZ3efZKdLQzHSyjt45eawWMxch+/DMugFVDKqpIB6npMfkF3q5REmnLxOrzbqzqCVV1A6mcelSdi1eTFtLC97ZiVVWUtvURH1zswJfhOxTJu872MHlZ4zm5v91ar/pCtWRfr4sOj2YMjmTYWNg3h/y93xHQVfkivRUbsMz09GQzWIrk4nVFPgiR+rvwzPTUejngyZWEylb/Xl4Zjrq08+XEppYTUQy6z5S50B75k/p5T48Mx2Ffr5Muwk2/yrzxGpNXw6vNhEBMo/USac/DM9MR6GfLzVDgmGZIUysJiJ9k2mkTtLg6gqGDKgq6eGZx6p//TRhqxkSDMucdWdRJ1YTkb7Jtrg5wPDB1fzXDVP7VR9+dzqRWyj9+E0jUo7i8XifRur0dzrSzycd3YuUlO4nbXe3tWdt399G6qSj0D9WIcy3o4usRHqXy0lb6J8jddJR6B+LTPPtPH9PMJInj/PtaNlDkdz0dtI2VX8dqZOOQv9YFGm+nUwTp+25/wEOrHtWyx6KdOPuvZ60rYgZ9bU1/XqkTjqh/oRmNhn4PHAmMBK43N2X97LPFODbwGnAa8Dt7n5vIevM6OXH6BH4SR4Phm7mIfR3LV7Uc6ZMgHicQ5s2sWvxYi17KJF3ZP99nF37s/ff19fW8MwXpxKLRWs8S9g/7WDgJeDavjQ2s5OAXwJPApOARcCPzWxmgerLrIjz7WjZQ5Hskv33S9ZtYeued9ixr73Xbp3KCotc4EPIR/ru/hjwGPR5MrDPAJvdfX7i+w1m9iHgX4CVBSkykyLNt5PLsoc6uStRlUv/PUTnpG065fZn7lxgdbdtKxPb0zKzGjOrS96A/A2pOWVWMM1C2ifOz3w7WvZQpHer1m/PKfCjctI2nXIL/RHAjm7bdgB1ZjYwwz4LgL0pt615q2baTcG8Ot2DP8/z7dQ2TYVMH0O17KFEVFt7J7f+4r/50Ndb2Lb3YNa2MYOGITWcMHwgnzq3kYeuOT8SJ23TicJPvZDgxG/SEPIV/EWab6e+eR4H1j2bcdnD+ubmvDyPSLnIdQz+qGED+/30Cn1VbqG/HejeEdcA7HP3tEvfuHs7cPg0ft7/06trCz7fjpY9FDlSrmPwZ0xoUOAnlFvorwO6d5TPSGwvnhCuwq2oHRwMy7zxSzppK5HX2xj8pKj336cT9jj9WmBcyqaTzGwS8Fd3/4uZLQRGu/snE/f/CLjOzP4V+L9AE/Bx4O+LVnQRr8LNRIEvUdaXJQ5jBiOHDuB/njoiMhdd9VXYr8RZBGPuk5J97/cBcwku2BqTvNPdN5vZ3wPfAZoJ+ub/2d2LN1yzSFfhikh6ZtbrEoejhg7kqS9qgEM6oY7ecfc17m5pbnMT98919ylp9jnD3Wvc/eSiX43bl6twRaSgpk9oIJYh92MGM06N5hj8vii3IZvhKuJVuCKS2edmjmfc8bU9gl99+L0Lu3unvBT4KlydoBXpm9qaSh665ny+tbKVVRt20NnlkZs47WjplcnVKbOCk7aepovnKK7C1ZTJIkentqaSWy45jVsuOU0HTDkwj1hXRGIqhr179+6lrq4u9wdIjt7Z3Xpk8Cevws1h9E6mKZOJxageO1ZTJotIn+zbt4+hQ4cCDHX3fdnaqk8/V8mrcM++CoaNgSEjg69nX5XzcM2+TJksIpJPOtI/VsdwFe7GadPpeP31jPdXjR7NuCe6zy8nInIkHekX0zGctO3rlMkiIvmi0A+JpkwWkTAo9EOkKZNFpNgU+iGqb55H9dixPYNfUyaLSIEo9EOUnDJ5+Jw5VI0eTWVDA1WjRzN8zhwaH9RwTRHJP43eKSG6wEREjoZG75QpBb6IFJpCX0QkQhT6RRa17jQRKS2acK0INKmaiJQKncgtME2qJiKFphO5JUSTqolIKVHoF1hby5M9Az8pHqetpaW4BYlIpCn0C0iTqolIqVHoF5AmVRORUqPQLzBNqiYipUShX2CaVE1ESolCv8A0qZqIlBKN0y8yTaomIvmmcfolIt0fVAW+iIRJ0zDkmaZcEJFSpu6dPNKUCyISBnXvhERTLohIqVPo55GmXBCRUqfQzxNNuSAi5UChnyeackFEyoFCP4805YKIlDqFfh5pygURKXUK/TyKDR6kKRdEpKSVxMVZZnYt8HlgBPAS8P+6+3MZ2s4F/r3b5nZ3H1DQIjPIdDHWSStWEBs8SH34IlJSQg99M7sS+DbwGeDXwDxgpZmNd/edGXbbB4xP+T6UITGZLsbac/8DHFj3rC7GEpGSUwrdO/8HuMfd/93d1xOE/9vAp7Ps4+6+PeW2oyiVdqOLsUSk3IQa+mZWDZwJrE5uc/d44vtzs+xaa2avmtlrZrbCzE7L8hw1ZlaXvAFD8lW/LsYSkXIT9pH+cUAF0P1IfQdB/346rQSfAi4F/pHgZ3jGzE7I0H4BsDfltvUYawZ0MZaIlKewQz9n7r7O3Ze4+4vuvhb4CLALuDrDLguBoSm3TH8ccqKLsUSkHIUd+ruBLqCh2/YGYHtfHsDdO4DfAeMy3N/u7vuSN2D/MdR7BF2MJSLlJtTQd/dDwAvAtOQ2M4slvl/Xl8cwswrgdGBbIWrMRhdjiUi5CX3IJsFwzfvM7DfAcwRDNgeTGItvZkuA1919QeL7m4FngY3AMILx/X8L/LjYhSfXv921eDFtLS14ZydWWUltUxP1zc0arikiJSf00Hf3ZWZWD9xGcPL2ReCilGGYY4DUITLDgXsSbfcQfFI4LzHcs+gqagcz4sYvwY1f0vq3IlLytHKWiEiZ08pZIiKSlkJfRCRCFPoiIhGi0BcRiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IiIRotAXEYkQhb6ISIQo9EVEIkShLyISIQp9EZEIUeiLiESIQl9EJEIU+iIiEaLQFxGJEIW+iEiEKPRFRCJEoS8iEiEKfRGRCFHoi4hEiEJfRCRCFPoiIhGi0BcRiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IiIRUhKhb2bXmtkWMztoZr82s7N7aX+Fmf0p0f4PZnZxIep6+OGHOXjwYNr7Dh48yMMPP1yIpxURKZjQQ9/MrgS+DXwF+ADwErDSzI7P0P484EHgJ8AZwHJguZlNzHdtM2bM4O677+4R/AcPHuTuu+9mxowZ+X5KEZGCMncPtwCzXwPPu/t1ie9jwGvA99z962naLwMGu/uHU7Y9C7zo7p/pw/PVAXv37t1LXV1dr/UlA/7qq69mwIABPb4XEQnbvn37GDp0KMBQd9+XrW2oR/pmVg2cCaxObnP3eOL7czPsdm5q+4SVmdqbWY2Z1SVvwJBcahwwYABXX301d999N2+99ZYCX0TKWmXIz38cUAHs6LZ9B/C+DPuMyNB+RIb2C4BbjrZACIL/U5/6FGeeeSYvvPCCAl9EylboffpFsBAYmnI7IdcHOHjwIPfddx8vvPAC9913X8aTuyIipS7s0N8NdAEN3bY3ANsz7LM9l/bu3u7u+5I3YH8uBab24Q8bNuxwV4+CX0TKUaih7+6HgBeAacltiRO504B1GXZbl9o+YUaW9kct3Unb1D5+Bb+IlJuwj/QhGK75/5jZp8xsAvBDYDDw7wBmtsTMFqa0XwxcZGbzzex9ZnYrcBZwV74LW7VqVdqTtsngX7VqVb6fUkSkoEIfsglgZtcBnyc4GfsicL27/zpx3xpgi7vPTWl/BXA70Ai8Atzg7o/28blyGrIpIlLqchmyWRKhX0wKfRHpb8pmnL6IiBSXQl9EJEIU+iIiEaLQFxGJEIW+iEiEKPRFRCIk7AnXQrNvX9ZRTSIiZSOXPIviOP3RwNaw6xARKYAT3P31bA2iGPoGjCLHidcI5uHfSjBLZ677hk21h0O1hyOqtQ8B3vBeQj1y3TuJFyTrX8J0gr8VAOzv7Yq3UqPaw6HawxHh2vvUXidyRUQiRKEvIhIhCv2+awe+kvhablR7OFR7OFR7FpE7kSsiEmU60hcRiRCFvohIhCj0RUQiRKEvIhIhCv0UZnatmW0xs4Nm9mszO7uX9leY2Z8S7f9gZhcXq9Y0tfS5djOba2be7XawmPWm1DLZzB42szcSdVzWh32mmNlvzazdzDaa2dzCV5q2jpxqT9Td/XV3MxtRpJKTdSwws+fNbL+Z7TSz5WY2vg/7hf5+P5raS+X9bmafNbPfm9m+xG2dmc3qZZ+8v+YK/QQzuxL4NsFwqQ8ALwErzez4DO3PAx4EfgKcASwHlpvZxKIUfGQtOdWesA8YmXL720LXmcFggnqv7UtjMzsJ+CXwJDAJWAT82MxmFqi+bHKqPcV4jnztd+a5rt5cCHwfOAeYAVQBj5vZ4Ew7lND7PefaE0rh/b4V+CJwJnAW0AKsMLPT0jUu2Gvu7roFw1Z/DdyV8n2MYLqGL2Zovwx4pNu2Z4EflUHtc4G3wn7N09TlwGW9tLkT+GO3bUuB/yyD2qck2g0L+7XuVld9oq7JWdqUzPv9KGovyfd7ora/Av+7mK+5jvQBM6sm+Ou7OrnN3eOJ78/NsNu5qe0TVmZpXxBHWTtArZm9amavmVnGo40SVBKv+zF60cy2mdkqMzs/7GKAoYmvf83SplRf977UDiX2fjezCjObTfBpcV2GZgV5zRX6geOACmBHt+07gEz9rSNybF8oR1N7K/Bp4FLgHwneB8+Y2QmFKjKPMr3udWY2MIR6crEN+Azw0cTtNWCNmX0grILMLEbQRfa0u/8xS9NSeb8flkPtJfN+N7PTzayN4IrbHwGXu/v6DM0L8ppHbpZNAXdfR8rRhZk9A2wArgZuCquu/s7dWwkCKOkZMzsZ+Bfgn8Kpiu8DE4EPhfT8x6JPtZfY+72V4FzUUOBjwH1mdmGW4M87HekHdgNdQEO37Q3A9gz7bM+xfaEcTe1HcPcO4HfAuPyWVhCZXvd97v5OCPUcq+cI6XU3s7uADwNT3b23hYVK5f0O5Fz7EcJ8v7v7IXff6O4vuPsCgoEAzRmaF+Q1V+gT/EcALwDTktsSHx2nkbm/bV1q+4QZWdoXxFHWfgQzqwBOJ+h+KHUl8brn0SSK/Lpb4C7gcqDJ3Tf3YbeSeN2Psvbuj1FK7/cYUJPhvsK85mGfvS6VG3AlcBD4FDABuBvYAzQk7l8CLExpfx7QAcwH3gfcChwCJpZB7TcD/xMYSzDE80HgHeDUEGqvJQi+SQSjMP4l8e8xifsXAktS2p8EHAD+NfG6XwN0AjPLoPZ5BP3K4wi6JRYRfEqbVuS6fwC8RTD8cUTKbWBKm5J8vx9l7SXxfk+8HyYDjQR/dBYCcWBGMV/zov6SlPoNuA54leAky6+BD6bctwa4t1v7Kwj66NqBPwIXl0PtwHdS2m4nGPd+Rkh1T0kEZvfbvYn77wXWpNnnd4n6/wzMLYfagRuAjYnAeZPgWoOpIdSdrmZPfR1L9f1+NLWXyvudYLz9lkQdOwlG5swo9muuqZVFRCJEffoiIhGi0BcRiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IgVWKss7ioBCX6SgSmx5RxFNwyByLMysHvgD8F13vyOx7TyCeVRmEUz09ffuPjFln6UESyZeVPyKJep0pC9yDNx9F8GqTLea2VlmNgT4KcGaxU9QussMSkRp5SyRY+Tuj5rZPcD9wG8Ipn5ekLg76/KOXp4Lv0gZ05G+SH58juAg6gpgjru3h1yPSFoKfZH8OBkYRfA71Ziyvb8t7yhlTt07IsfIzKqBnwHLCBa8+LGZne7uOwmWtru42y7lvLyjlDmN3hE5Rmb2DeBjwP8A2oC1wF53/3BiyOYfge8D/xdoAr5LMKJnZUglS4Qp9EWOgZlNAVYRLHv4VGJbI/AS8EV3/2GizXeAU4GtwFfd/d7iVyui0BcRiRSdyBURiRCFvohIhCj0RUQiRKEvIhIhCn0RkQhR6IuIRIhCX0QkQhT6IiIRotAXEYkQhb6ISIQo9EVEIkShLyISIf8/DkIINNDELdQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "historyPG = onp.asarray(historyPG)\n",
    "updatesPG = onp.asarray(updatesPG) \n",
    "\n",
    "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
    "axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='#1F77B4', label='GD')\n",
    "axes.scatter(historyNt[:,0], historyNt[:,1], lw=0.5, color='#FF7F0E', label='Newton')\n",
    "axes.scatter(historyPG[:,0], historyPG[:,1], lw=0.5, color='#D62728', label='PG')\n",
    "axes.scatter([0], [0], lw=0.25, color='black', marker='x') # target at 0,0\n",
    "axes.set_xlabel('x0'); axes.set_ylabel('x1'); axes.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This illustrates that the inverse simulator variant, PG in red, does even better than Newton's method in orange. It yields a trajectory that is better aligned with the ideal _diagonal_ trajectory, and its final state is closer to the origin. A key ingredient here is the inverse function for $\\mathbf{y}$, which provided higher order terms than the second-order approximation for Newton's method. This improves the scale-invariance of the optimization. Despite the simplicity of the problem, Newton's method has problems finding the right search direction. For the inverse simulator update, on the other hand, the higher order information yields an improved direction for the optimization.\n",
    "\n",
    "This difference also shows in first update step for each method: below we measure how well it is aligned with the diagonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diagonal lengths (larger is better): GD 1.053930, Nt 1.264911, PG 1.356443 \n"
     ]
    }
   ],
   "source": [
    "def mag(x):\n",
    "    return np.sqrt(np.sum(np.square(x)))\n",
    "\n",
    "def one_len(x):\n",
    "    return np.dot( x/mag(x), np.array([1,1])) \n",
    "\n",
    "print(\"Diagonal lengths (larger is better): GD %f, Nt %f, PG %f \" % \n",
    "      (one_len(updatesGD[0]) , one_len(updatesNt[0]) , one_len(updatesPG[0])) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The largest value of 1.356 for PG confirms what we've seen above: the PG gradient was the closest one to the diagonal direction from our starting point to the origin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "<br>\n",
    "\n",
    "## y Space\n",
    "\n",
    "To understand the behavior and differences of the methods here, it's important to keep in mind that we're not dealing with a black box that maps between $\\mathbf{x}$ and $L$, but rather there are spaces in between that matter. In our case, we only have a single $\\mathbf{y}$ space, but for DL settings, we might have a large number of latent spaces, over which we have a certain amount of control. We will return to NNs soon, but for now let's focus on $\\mathbf{y}$. \n",
    "\n",
    "A first thing to note is that for PG, we explicitly map from $L$ to $\\mathbf{y}$, and then continue with a mapping to $\\mathbf{x}$. Thus we already obtained the trajectory in $\\mathbf{y}$ space, and not coincidentally, we already stored it in the `historyPGy` list above.\n",
    "\n",
    "Let's directly take a look what the inverse simulator did in $\\mathbf{y}$ space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc742e83a90>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAGDCAYAAAAyH0uIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhvElEQVR4nO3df3yedX3v8dcnd5sUmh/NYxRCRznWFrsjjuNUoG4eRZ1D0eOUopS6naH7oUydqGebrQ4Lx8lRfDhkcOAMmYVtlbpNQfEH8kMnKEqsOkUmEcaPIr9aW3OHpk1I8j1/3HfKnZikSZvkur/J6/l4XI803/t73dfnvpq+c/V7fa/ripQSkqT8NBRdgCTp4BjgkpQpA1ySMmWAS1KmDHBJypQBLkmZMsAlKVMGuCRlygCXpEwZ4JKUKQNckjJlgEtSpgxwScqUAa5sRMRLIyJFxOvHeG199bUXTrD+woj4YET8NCL2RcTPI+L2iHhFTZ/NEfFkRDwzIm6MiD0R8UhEnBcRMer9/ldEfKv6PnsjYltEnDHOtn8vIu6MiN6I2B0R34iI3xnV51URcVt1mz0R8cWIOH7qe0rzhQGunHwd2A68aYzX3gTcl1K6Y4L1NwEfBL4GvAP4a+Ah4Hmj+pWArwCPA38BbAPOry613gV8HzgP2AgMAP8cEa+u7RQRHwT+AXiq2veD1c/xspo+vw98EXgS+EvgfwPPBm6PiGdM8Jk0n6WUXFyyWYAPA/uAtpq2pVTCcdMB1v0BcMMB+mwGEnBJTVsANwB9wBE17YeNWnch8CPglpq2VcAg8FmgYVT/qH5tBnYDfzfq9aOAX4xud3EZXjwCV26uAZqA2qGKM4EFwD8eYN1fAMdHxHGT2M6lw39IKaXq943Ab9e07x3+c0S0A23AbYw8on8dlf/pXpBSGqrdQPV9AV4BLAE+HRFHDC9Ugv87wEsnUa/moQVFFyBNRUrpJxHRSWXI5Kpq85uAb6eU7j3A6ucB1wNdEXEXlWGSf0gp/XBUvyHgP0e1dVW/PmO4ISJeA3wAeC6VXyr7y6z588rq+909QV3Dv1BuHef18gTrah4zwJWja4BPRMQxVIJzDZUx7QmllL4RESuB3wV+B/gj4N0R8baU0ienUkBE/Hfg88A3gD8FHqUyjPNmYP1U3ounz0X9PvDYGK8PTPH9NE8Y4MrRtcDHgbOAw6gE59bJrJhS2gV8CvhURDRTCeBNQG2ANwDP5OmjboBnVb8+UP26lspY/Kkppb7hThHx5lGbvK/6fs+mMgY/lvuqX59IKd08mc8hgbNQlKGU0k7gy8DvURk++Uq1bUIR8Suj3udJ4F5GDn8M239EX50++A4qvyhuqTYPUhkqKdX0ewaVMe9a11EZQjkvIkb8e6uZlngjlWGSjRGxcIy6l074wTRveQSuXF0D/Ev1z381yXXujoivU5kWuAt4AZWToZeO6rcPeGVEXE3lJOKrgFcDH04p7aj2+SLwHuArEbEFOBJ4O5VfCCcMv1FK6d6I+OtqjbdFxGepzGY5EXgE2JBSKkfEOVSmGn4vIq4FdgDHVrf7TSYxRKT5J54+ES7lIyIaqYwXNwAdKaV9k1jn/cBrqQyHNAEPUgnNi1JKT1X7bKYS6icAlwMvAnqAKxg1kyQi3gK8j0rQ3g98hMpJzg+mlEZf9PNm4J1UhlJ6gR8CH6odMomIU6rvt6Za38+ozGq5NKW0bZK7RvOIAa4sRcQCKkewX0gp/eE0vu9m4IyUUvN0vac0UxwDV65eR+UCnmsKrkMqjGPgykpEnExleOOvgO+nlP6t4JKkwngErtycQ2Vs+gngfxZci1Qox8AlKVMegUtSpgxwScpU1icxq1eyLaMyT1eS5ooW4JF0gDHurAOcSng/XHQRkjQDjqFyMde4cg/wHoDt27fT2tpadC2SdMjK5TLLly+HSYws5B7gALS2thrgkuYdT2JKUqYMcEnKlAEuSZmaE2PgE0kpMTAwwODgYNGlFKpUKrFgwQKefoaApNzN6QDv7+/n0Ucfpbe3t+hS6sLhhx/O0UcfTWNjY9GlSJoGczbAh4aGuP/++ymVSixbtozGxsZ5e/SZUqK/v58dO3Zw//33c9xxx9HQ4OiZNNNSSjOaO3M2wPv7+xkaGmL58uUcfvjhRZdTuMMOO4yFCxfy4IMP0t/fz6JFi4ouSZrT0uAgD7zhjTzjnz9DlEoHXuEgzPnDMI80n+a+kGZP73e3sa+ri95tM/c0vDl7BC5Js21g926GeioXUO7avBkGBth19TUs7OgAoKGlhQXt7dO2PQNckqbJzsuvYPeWLZRaW0kDAwD0dnbywLqzGCyXaV+/no6NG6Zte/6fWpKmScfGDXRs2gSlEkPlMgBD5TKxoMTRF5w/reENBnhdOvvss4kIIoLGxkZWrVrFBRdcwED1N3pKiSuvvJIXvvCFtLa20tzczPHHH8+73vUu7r333oKrl+a39jPWsnjNySPaDj95DUtOP33at2WAT8FsPn7ula98JY8++ig//elPee9738umTZu46KKLSCmxfv16/uzP/ozTTjuNr371q9x9991cddVVLFq0iA996EOzVqOkX5ZSovfOTkrt7Sw64QRK7e30dnbOSH44Bj5JszElqFZTUxMd1RMf55xzDp/73Of4/Oc/z4oVK7j22mu5/vrree1rX7u//7HHHsuaNWtm9ZeMpF/W19XF0N69/OrffJzm3/otnrz9dn727vfQ19XFotWrp3VbBvgk1U4JWnzSSbO+/cMOO4yf//znfPrTn2b16tUjwrvWfL1YSaoXjStWsOrWWyg1NwPQ/KIXserWW4impmnflkMoExjYvZv+hx6i/6GHRkwJGm4b2L17xmtIKXHzzTdz44038rKXvYyuri5Wj/otfu6559Lc3ExzczPHHHPMjNckaXwNjY37w3tYqaWFhhm4hYUBPoGdl1/Bfae9mgfWnbV/Mv7wlKD7Tns1Oy+/Ysa2fcMNN9Dc3MyiRYt41atexZlnnsmmTZvG7Pv+97+fH/zgB5x33nk8+eSTM1aTpPriEMoEOjZuoOlZz2LHJz7B0K5dQGVK0IIjl3L0BefPyFnlYS996Uu5/PLLaWxsZNmyZSxYUPmrOu6447jnnntG9F26dClLly7lyCOPnLF6JNUfj8APYDanBNVavHgxq1at4thjj90f3gBnnXUW99xzD9dff/2Mbl9S/fMI/ABqpwQtXL6cp7Zv3z8lqIgThuvWreOzn/0s69atY8OGDZx66qkcddRRPPjgg2zdupXSLMyQkVQfCj8Cj4hfjYh/jIifR8TeiPhRRLyg6LqGDU8JWvaxi1jxma0su+ijDO3ZQ19XVyH1RARbt27l4osv5ktf+hIvf/nLWb16NW95y1tYvnw5t99+eyF1SZp9UeS84YhoB74PfA24HNgBHAfcl1K6bxLrtwLd3d3dv/RU+n379nH//fezYsWKQ7p16lB/P6m/f8RZ5cGeHqKpaUbOKs+k6donkmZOuVymra0NoC2lVJ6ob9FDKH8JbE8pvbmm7f6iihlLQ2MjjArqUktLQdVI0tOKHkJ5LfDdiPjniHgiIr4fEX88XueIaIqI1uEFMEklzVtFB/gzgXOAnwKnUhlGuSQi/mCc/huA7prl4dkoUpLqUdEB3gB8L6W0MaX0/ZTS3wFXAm8bp/+FQFvN4mWHkuatosfAHwXuHtX2H8DasTqnlPqAvuHvJzONz5s7Pc19Ic0tRR+BfxMYfXuuZwEPHuobL1y4EIDe3t5Dfas5Y3hfDO8bSXkr+gj8b4BvRcRG4DPAScCfVJdDUiqVWLJkCU888QQAhx9++Ly9U19Kid7eXp544gmWLFnixT7SHFFogKeUOiPi9VTGts+jMoXw3JTSP03H+w/fT3s4xOe7JUuW7N8nkvJX6IU8h2qiC3lqDQ4O8tRTT81eYXVo4cKFHnlLGcjpQp5ZUSqVDC9Jc07RJzElSQfJAJekTBngkpQpA1ySMmWAS1KmDHBJypQBLkmZMsAlKVMGuCRlygCXpEwZ4JKUKQNckjJlgEtSpgxwScqUAS5JmTLAJSlTBrgkZcoAl6RMGeCSlCkDXJIyZYBLUqYMcEnKlAEuSZkywCUpUwa4JGXKAJekTBngkpQpA1ySMmWAS1KmDHBJypQBLkmZKjTAI2JTRKRRy0+KrEmScrGg6AKAHwO/XfP9QFGFSFJO6iHAB1JKjxVdhCTlph7GwI+LiEci4j8j4p8i4tjxOkZEU0S0Di9AyyzWKUl1pegA/w5wNvBK4BxgBXBbRIwXzBuA7prl4VmoUZLqUqSUiq5hv4hYAjwIvCeldNUYrzcBTTVNLcDD3d3dtLa2zk6RkjSDyuUybW1tAG0ppfJEfethDHy/lNIvIqILWDXO631A3/D3ETFbpUlS3Sl6CGWEiGgGVgKPFl2LJNW7oueBfywiXhIRz4iI3wQ+BwwCny6yLknKQdFDKMdQCetfAXYAtwNrUko7Cq1KkjJQaICnlNYVuX1JylldjYFLkibPAJekTBngkpQpA1ySMmWAS1KmDHBJypQBLkmZMsAlKVMGuCRlygCXpEwZ4JKUKQNckjJlgEtSpgxwScqUAS5pv3p6Rq4OzACXBEAaHOSBtWeQBgeLLkWTZIBLAqD3u9vY19VF77ZtRZeiSSr6kWqSCjSwezdDPT0A7Nq8GQYG2HX1NSzs6ACgoaWFBe3tBVaoiRjg0jy28/Ir2L1lC6XWVtLAAAC9nZ08sO4sBstl2tevp2PjhoKr1HgcQpHmsY6NG+jYtAlKJYbKZQCGymViQYmjLzjf8K5zBrg0z7WfsZbFa04e0Xb4yWtYcvrpBVWkyXIIRZrnUkr03tlJqb2dhcuX89T27fR2dpJSIiKKLk8T8Ahcmuf6uroY2ruXZR+7iBWf2cqyiz7K0J499HV1FV2aDiBynrgfEa1Ad3d3N62trUWXI2VpqL+f1N9Pqbl5f9tgTw/R1ERDY2OBlc1P5XKZtrY2gLaUUnmivg6hSPNcQ2MjjArqUktLQdVoKhxCkaRMGeCSlCkDXJIyZYBLUqYMcEnKlAEuSZkywCUpUwa4JGWqbgI8It4XESkiLi66FknKQV0EeEScCLwV+GHRtUhSLgoP8IhoBv4J+GNgd8HlSFI2Cg9w4DLgiymlmw/UMSKaIqJ1eAG8YYOkeavQm1lFxDrgecCJk1xlA/DBmatIkvJR2BF4RCwHPgG8KaW0b5KrXQi01SzHzFB5klT3ijwCfz5wJPC9mqd+lIAXR8Q7gKaU0mDtCimlPqBv+HufFiJpPisywG8Bfn1U26eAnwAfGR3ekqSRCgvwlFIPcFdtW0TsAX6eUrpr7LUkScPqYRaKJOkg1NUj1VJKpxRdgyTlwiNwScqUAS5JmTLAJSlTBrgkZcoAl6RMGeCSlCkDXJIyZYBLUqYMcEnKlAEuSZkywCUpUwa4JGXKAJekTBngkpQpA1ySMmWAS1KmDHBJypQBLkmZMsAlKVMGuCRlygCXpEwZ4JKUKQNckjJlgEtSpgxwScrUtAZ4RPy3iBiczveUJI1tJo7AYwbeU5I0yoKpdI6Izx6gSxuQDr4cqb6klIjwmET1aapH4P8DWAR0j7M8Oa3VSQVKg4M8sPYM0qCjgqpPUzoCB/4D+NeU0lVjvRgRzwVec6hFSfWg97vb2NfVRe+2bSw+6aSiy5F+yVQDfBvwPGDMAAf6gIcOqSKpQAO7dzPU0wPArs2bYWCAXVdfw8KODgAaWlpY0N5eYIXS0yKlyQ9ZR0QTUEop9c5cSZMXEa1Ad3d3N62trUWXozngsQ9fyO4tWyi1tpIGBhgql2lobSUWLGCwXKZ9/Xo6Nm4oukzNYeVymba2NoC2lFJ5or5TGgNPKfWllHoj4taI+ODo1yOiPSJunez7RcQ5EfHDiChXlzsi4lVTqUmaTh0bN9CxaROUSgyVK/92hsplYkGJoy843/BWXTnYaYSnAO+IiOsiYnFNeyPwkim8z8PA+4DnAy8AbgWuj4jjD7Iu6ZC1n7GWxWtOHtF2+MlrWHL66QVVJI3tUOaB/zbQAXw7Ip5xMG+QUvpCSulLKaWfppS6UkrvpzKTZc0h1CUdkpQSvXd2UmpvZ9EJJ1Bqb6e3s5OpDDdKs+FQAvxRKkfbPwI6I+KUQykkIkoRsQ5YDNwxTp+miGgdXoCWQ9mmNJa+ri6G9u5l2ccuYsVntrLsoo8ytGcPfV1dRZcmjTClk5j7V6pcLn90SumJ6vcfAD4AfAT4QEqpNIX3+nUqgb2IytH3+pTSl8bpuwn4pbF3T2JqOg3195P6+yk1N+9vG+zpIZqaaGhsLLAyzQdTOYl5sAE+BHQMB3i1bS1wNXDYFAO8ETiWylWcZwB/BLwkpXT3GH2bgKaaphbgYQNc0lwxlQCf6jzwYSuAHbUNKaV/jYifUDkZOWkppX7g3uq32yLiROBdwFvH6NtHZa45gJc4S5rXDirAU0oPjtP+Y+DHh1RRZVy+6YC9JGmeO9gj8GkRERcCX6Zy9WYLsJ7KFMVTCyxLkrJQaIADRwLXAEdTuRnWD4FTU0o3FVqVJGWg0ABPKf1hkduXpJz5SDVJypQBLkmZMsAlKVMGuCRlygCXpEwZ4JKUKQNckjJlgEtSpgxwScqUAS5JmTLAJSlTBrgkZcoAl6RMGeCSlCkDXJIyZYBLUqYMcEnKlAEuSZkywCUpUwa4JGXKAJekTBngkpQpA1ySMmWAS1KmDHBJypQBrsKklIouQcqaAa5CpMFBHlh7BmlwsOhSpGwZ4CpE73e3sa+ri95t24ouRcrWgqIL0PwxsHs3Qz09AOzavBkGBth19TUs7OgAoKGlhQXt7QVWKOXFANes2Xn5FezesoVSaytpYACA3s5OHlh3FoPlMu3r19OxcUPBVUr5cAhFs6Zj4wY6Nm2CUomhchmAoXKZWFDi6AvON7ylKTLANavaz1jL4jUnj2g7/OQ1LDn99IIqkvLlEIpmVUqJ3js7KbW3s3D5cp7avp3ezk5SSkRE0eVJWSn0CDwiNkREZ0T0RMQTEXFdRKwusibNrL6uLob27mXZxy5ixWe2suyijzK0Zw99XV1FlyZlJ4q8mCIivgJcC3RS+d/Ah4HnAM9OKe2ZxPqtQHd3dzetra0zWqumx1B/P6m/n1Jz8/62wZ4eoqmJhsbGAiuT6kO5XKatrQ2gLaVUnqhvoUMoKaVX1n4fEWcDTwDPB75RRE2aWQ2NjTAqqEstLQVVI+Wt3sbA26pfd431YkQ0AU01Tf7LlzRv1c0slIhoAC4GvplSumucbhuA7prl4dmpTpLqT90EOHAZlfHvdRP0uZDKUfrwcsws1CVJdakuhlAi4lLgNcCLU0rjHlWnlPqAvpr1ZqE6SapPhQZ4VBL4b4HXA6eklO4vsh5JyknRR+CXAeuB3wV6IqKj2t6dUtpbXFmSVP+KHgM/h8pY9teBR2uWMwusSZKyUPQ8cAexJekgFX0ELkk6SAa4JGXKAJekTBngkpQpA1ySMmWAS1KmDHBJypQBLkmZMsAlKVMGuIDKw4Yl5cUAF2lwkAfWnkEaHCy6FElTYICL3u9uY19XF73bthVdiqQpKPp2sirIwO7dDPX0ALBr82YYGGDX1dewsKNyR9+GlhYWtLcXWKGkAzHA56mdl1/B7i1bKLW2kgYGAOjt7OSBdWcxWC7Tvn49HRs3FFylpIk4hDJPdWzcQMemTVAqMVQuAzBULhMLShx9wfmGt5QBA3weaz9jLYvXnDyi7fCT17Dk9NMLqkjSVDiEMo+llOi9s5NSezsLly/nqe3b6e3sJKXkA6OlDHgEPo/1dXUxtHcvyz52ESs+s5VlF32UoT176OvqKro0SZMQOV/AERGtQHd3dzetra1Fl5Odof5+Un8/pebm/W2DPT1EUxMNjY0FVibNX+Vymba2NoC2lFJ5or4OocxjDY2NMCqoSy0tBVUjaaocQpGkTBngkpQpA1ySMmWAS1KmDHBJypQBLkmZMsAlKVMGuCRlygCXpEwZ4JnK+RYIkqaHAZ4hn2EpCQzwLPkMS0ngzayy4TMsJY1WaIBHxIuBPweeDxwNvD6ldF2RNdUrn2EpabSih1AWA/8OvL3gOuqez7CUNFqhR+AppS8DXwYm9QiviGgCmmqa5tXNq9vPWEvvd75N+Qs37G/zGZbS/FX0EfhUbQC6a5aHiy1ndtU+w3LRCSdQam/f/wxLSfNPbgF+IdBWsxxTbDmzy2dYSqpVN8/EjIjEFE9izrdnYvoMS2nu85mYc5TPsJRUK7chFElSVdHzwJuBVTVNKyLiucCulNJDxVQlSXkoegjlBcDXar7/ePXr1cDZs16NJGWk6HngXwcOPAF8jkopTWr+uySNxTHwgnhHQUmHygAviHcUlHSoih4Dn1e8o6Ck6WSAzyLvKChpOjmEMou8o6Ck6WSAz7L2M9ayeM3JI9q8o6Ckg+EQyiyrvaPgwuXLeWr79v13FHRKoaSp8Ah8lnlHQUnTpW7uRngw6vVuhBMdTXtHQUkTmcrdCD0Cn2YHukCnobFxRHhD5Y6ChrekqTLAp5kX6EiaLZ7EnAZeoCOpCAb4NPACHUlFcAhlGniBjqQiGODTxAt0JM02h1AO0uipgl6gI2m2eQR+EMaaKugFOpJmm0fgB6F2quDik04CoHHFClbdesv+Od7NL3oRq269hWhqKrJUSXOYAT5JBzNVsNTSMttlSppHDPBJcqqgpHrjGPgkHbXhfU4VlFRXDPBJGD5pueT1r3OqoKS64RDKJOw/afndbU4VlFQ3DPBxjHXScudllzG4Zw8d799I8ymnsO/HP+Zn734PfV1dLFq9utiCJc073g98HI99+MIRJy2HymUaWlpgwQKGenr2n7T0Xt6SptNU7gdugE9g97/8Kzs+8QkGd+zY37bgyKUsPfdcx70lzQgf6DBNvL+JpHrmGPgEvL+JpHrmEfgEvL+JpHrmGPgEfACxpNk2lTFwh1Am0NDYCKOC2vubSKoXDqFIUqbqIsAj4u0R8UBE7IuI70TESTO9zfPOO2+mNyFJM6rwAI+IM4GPA+cDzwP+HbgxIo6cqW3u2LGDq666ip07d87UJiRpxhUe4MB7gCtTSp9KKd0NvA3oBd4yUxu88soreeSRR/jkJz85U5uQpBlXaIBHRCPwfODm4baU0lD1+xeO0b8pIlqHF2DSZxS3bNnCypUrWblyJZdccgkAF1988f62LVu2HOrHkaRZVfQR+BFACXh8VPvjQMcY/TcA3TXLw5Pd0Pr167n22ms54ogjePzxyuYef/xxli5dytatW1m/fv3B1C9JhSk6wKfqQqCtZjlmKiufeOKJvPGNbxzR9oY3vIEXvOAF01agJM2WogN8JzAIHDWq/SjgsdGdU0p9KaXy8AL0THWDt912GytXruSd73wnK1eu5LbbbjuowiWpaIUGeEqpH9gGvHy4LSIaqt/fMd3b6+3tZefOndx0001ccskl3HTTTezcuZO9e/dO96YkacYVfil9dRrh1cBbgTuBc4E3Ar+WUho9Nj563SldSj/8WWtvRDVWmyQVJatL6VNKWyNiKXABlROXPwBeeaDwPhhjhbTBLSlXhQc4QErpUuDSouuQpJwUfRJTknSQDHBJypQBLkmZMsAlKVMGuCRlygCXpEzVxTTCQ1UuTzjXXZKyMZU8K/xKzEMREb/KFO5IKEkZOSal9LOJOuQe4AEsY+o3tWqhEvzHHMS6c4374mnui5HcH0+b7X3RAjySDhDQWQ+hVD/chL+hxlJz+XzPge41MNe5L57mvhjJ/fG0AvbFpLbhSUxJypQBLkmZmq8B3gecX/0637kvnua+GMn98bS63BdZn8SUpPlsvh6BS1L2DHBJypQBLkmZMsAlKVNzNsAj4u0R8UBE7IuI70TESQfo/4aI+Em1/48i4rTZqnWmTWVfRMTZEZFGLftms96ZEhEvjogvRMQj1c/1ukmsc0pEfC8i+iLi3og4e+YrnXlT3RfV/TD65yJFRMcslTxjImJDRHRGRE9EPBER10XE6kmsV3hmzMkArz7p/uNUpv08D/h34MaIOHKc/r8JfBq4CvgN4Drguoh4zqwUPIOmui+qysDRNct/mek6Z8liKp//7ZPpHBErgC8CXwOeC1wMfDIiTp2h+mbTlPZFjdWM/Nl4YprrKsJLgMuANcArgIXAVyNi8Xgr1E1mpJTm3AJ8B7i05vsGKpfcv2+c/luBG0a1fRu4oujPUsC+OBv4RdF1z8J+ScDrDtDnI8Bdo9quBb5SdP0F7ItTqv2WFF3vLOyPpdXP+uIJ+tRFZsy5I/CIaASeD9w83JZSGqp+/8JxVnthbf+qGyfon4WD3BcAzRHxYERsj4jrI+L4GS61Xs3Jn4tD9IOIeDQiboqI3yq6mBnSVv26a4I+dfGzMecCHDgCKAGPj2p/HBhvvK5jiv1zcTD74h7gLcDvAr9H5WfkWxFxzEwVWcfG+7lojYjDCqinSI8CbwPWVpftwNcj4nmFVjXNIqKBylDZN1NKd03QtS4yI+u7EWr6pZTuAO4Y/j4ivgX8B/BW4K+KqkvFSindQ+WX+7BvRcRK4N3A7xdT1Yy4DHgO8KKiC5mMuXgEvhMYBI4a1X4U8Ng46zw2xf65OJh9MUJK6Sng+8Cq6S0tC+P9XJRTSnsLqKfe3Mkc+rmIiEuB1wAvTSkd6EExdZEZcy7AU0r9wDbg5cNt1f8WvZyaI8tR7qjtX/WKCfpn4SD3xQgRUQJ+ncp/oeebOflzMY2eyxz4uYiKS4HXAy9LKd0/idXq42ej6DO+M3QW+UxgH/AHwH8F/h+wGziq+vo1wIU1/X8TeAp4L/BrwCagH3hO0Z+lgH1xHvA7wDOpTDv8NLAXeHbRn2Ua9kUzldB5LpVZBu+u/vnY6usXAtfU9F8B7AE+Wv25+FNgADi16M9SwL44l8p5kVVUhhgupvK/u5cX/VmmYV/8X+AXVKYTdtQsh9X0qcvMKHznzeBfyjuAB6nc/vE7wMk1r30d2Dyq/xuojPH1AXcBpxX9GYrYF8Df1PR9jMo86N8o+jNM0344pRpWo5fN1dc3A18fY53vV/fHfcDZRX+OIvYF8BfAvVR+mf+cytz4lxb9OaZpX4y1H1Lt33W9Zoa3k5WkTM25MXBJmi8McEnKlAEuSZkywCUpUwa4JGXKAJekTBngkpQpA1ySMmWAS1KmDHBpCiLihIi4rfocxO0R8RdF16T5ywCXJikiWoGvUrlXzPOBPwc2RcSfFFqY5i3vhSLViIhnAGPdTvTfqDwH8a+BjlS5VS8R8X+oPE/y12atSKnKI3BppO2MfOr6b1C5+943qDzv8BvD4V11I7A6Itpnu1DJAJdqpJQGU0qPpZQeo3KP6Cuo3KR/E+M/BxHyf36qMuQzMaXx/T3QArwipTQUEUXXI41ggEtjiIgPAKcCJ6WUeqrN4z0Hcfg1aVY5hCKNEhFrqTxa7o0ppftqXroDeHFELKxpewVwT0pp92zWKIGzUKQRIuI5VB4793HgspqX+qk8A/IeKlMJP0Ll2ZB/D7w7pfR3s1yqZIBLtSLibOBTY7z0bymlUyLiBCrBfiKwE/jblNJHZrFEaT8DXJIy5Ri4JGXKAJekTBngkpQpA1ySMmWAS1KmDHBJypQBLkmZMsAlKVMGuCRlygCXpEwZ4JKUqf8P7WKg3gtvxL4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "historyPGy = onp.asarray(historyPGy)\n",
    "\n",
    "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
    "axes.set_title('y space')\n",
    "axes.scatter(historyPGy[:,0], historyPGy[:,1], lw=0.5, color='#D62728', marker='*', label='PG')\n",
    "axes.scatter([0], [0], lw=0.25, color='black', marker='*') \n",
    "axes.set_xlabel('z0'); axes.set_ylabel('z1'); axes.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this variant, we're making explicit steps in $\\mathbf{y}$ space, which progress in a straight diagonal line to the origin  (which is likewise the solution in $\\mathbf{y}$ space).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, neither GD nor Newton's method give us information about progress in intermediate spaces (like the $\\mathbf{y}$ space). \n",
    "\n",
    "For GD we're concatenating the Jacobians, so we're moving in directions that locally should decrease the loss. However, the $\\mathbf{y}$ position is influenced by $\\mathbf{x}$, and hence we don't know where we end up in $\\mathbf{y}$ space until we have the definite point there. (For NNs in general we won't know at which latent-space points we end up after a GD update until we've actually computed all updated weights.)\n",
    "\n",
    "More specifically, we have an update $-\\eta \\frac{\\partial L}{\\partial \\mathbf{x}}$ for GD, which means we arrive at $\\mathbf{y}(\\mathbf{x} -\\eta \\frac{\\partial L}{\\partial \\mathbf{x}})$ in $\\mathbf{y}$ space. A Taylor expansion with \n",
    "$h = \\eta \\frac{\\partial L}{\\partial \\mathbf{x}}$ yields \n",
    "\n",
    "$\n",
    "\\quad\n",
    "\\mathbf{y}(\\mathbf{x} - h) = \n",
    "\\mathbf{y}(\\mathbf{x}) - h \\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{x}} + \\mathcal{O}( h^2 )\n",
    "= \\mathbf{y}(x) - \\eta \\frac{\\partial L}{\\partial \\mathbf{y}} (\\frac{\\partial \\mathbf{y}}{\\partial x})^2 + \\mathcal{O}( h^2 )\n",
    "$.\n",
    "\n",
    "And $\\frac{\\partial L}{\\partial \\mathbf{y}} (\\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{x}})^2$ clearly differs from $\\frac{\\partial L}{\\partial \\mathbf{y}}$, which we would apply with GD when optimizing for $\\mathbf{y}$ directly.\n",
    "\n",
    "**Newton's method** does not fare much better: we compute first-order derivatives like for GD, and the second-order derivatives for the Hessian for the full process. But since both are approximations, the actual intermediate states resulting from an update step are unknown until the full chain is evaluated. In the _Consistency in function compositions_ paragraph for Newton's method in {doc}`physgrad` the squared $\\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{x}}$ term for the Hessian already indicated this dependency.\n",
    "\n",
    "With **inverse simulators** we do not have this problem: they can directly map points in $\\mathbf{y}$ to $\\mathbf{x}$. Hence we know eactly where we started in $\\mathbf{y}$ space, as this position is crucial for evaluating the inverse.\n",
    "\n",
    "In the simple setting of this section, we only have a single latent space, and we already stored all values in  $\\mathbf{x}$ space during the optimization (in the `history` lists). Hence, now we can go back and re-evaluate `physics_y` to obtain the positions in $\\mathbf{y}$ space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.asarray([3.,3.])\n",
    "eta = 0.01\n",
    "historyGDy = []\n",
    "historyNty = []\n",
    "\n",
    "for i in range(1,10):\n",
    "    historyGDy.append(physics_y(historyGD[i]))\n",
    "    historyNty.append(physics_y(historyNt[i]))\n",
    "\n",
    "historyGDy = onp.asarray(historyGDy)\n",
    "historyNty = onp.asarray(historyNty)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc7430c4b20>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAGDCAYAAAA26pu1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0/klEQVR4nO3de3xU1b338c8vkxuQCxHBQMVy89Cq8IiiYhXwLrVaK6Ag1B57ztP2eDkV6zltsR6KPK3U6ssqxeqjeIq0B6XtOYVKn9YLFAqKgnipihAv3KzcE5OQkOus5489iZOQhFxmZu+dfN+v134Ns2bPzC+T8M3O2nutZc45REQkvNL8LkBERLpGQS4iEnIKchGRkFOQi4iEnIJcRCTkFOQiIiGnIBcRCTkFuYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhJyCXEQk5BTkEhpmdqGZOTO7poXHZsQeO7eN52eY2Q/N7D0zqzKzQ2a23swujdtnsZkdNrNhZvasmVWY2cdmNsfMrNnr/ZuZvRR7nSNmttnMprby3l81s41mVmlmJWb2VzO7rNk+XzSzdbH3LDezP5rZqR3/pKSnUZBLmKwBdgMzW3hsJvCBc25DG8+fC/wQ+AtwK/BjYBdwRrP9IsCfgX3Ad4HNwN2xLd5twOvAHOBOoA74rZl9KX4nM/sh8CugNrbvD2Nfx0Vx+9wA/BE4DHwP+D/AKcB6MxvSxtckAs45bdpCswH3AFVAflxbf7yQnHuM574BrDzGPosBByyIazNgJVANHB/X3qvZczOAt4BVcW0jgHrgf4C0Zvtb7DYHKAEea/b4CcAnzdu1aWu+6YhcwmYJkAXEd2FMA9KBXx/juZ8Ap5rZye14n4UN/3DOudj9TOCSuPYjDf82swIgH1hH0yP8r+D95TvPOReNf4PY6wJcCvQFnjKz4xs2vF8ArwAXtqNe6cHS/S5ApCOcc1vNbBNeV8oTseaZwMvOufeP8fQ5wAqgyMzexus++ZVz7m/N9osCHzZrK4rdDmloMLMrgbuA0/F+uTSWGffv4bHX29JGXQ2/WFa38nhZG88VUZBLKC0BHjKzE/ECdBxen3ebnHN/NbPhwNXAZcD/Bm43s39xzi3qSAFmNh74A/BX4GZgD173zteBGR15LT49V3UDsLeFx+s6+HrSwyjIJYyeBh4Argd64QXosvY80TlXDPwS+KWZ5eAF8VwgPsjTgGF8ehQO8A+x2x2x2yl4ffWXO+eqG3Yys683e8sPYq93Cl4ffUs+iN3ud8690J6vQySe+sgldJxzB4E/AV/F61b5c6ytTWbWr9nrHAbep2m3SIPGI/zYZYe34v3CWBVrrsfrQonE7TcEr0883nK8rpU5Ztbk/1vc5YzP4nWf3GlmGS3U3b/NL0x6PB2RS1gtAX4X+/d/tPM5W8xsDd7lhMXAWLyTpgub7VcFTDKzJ/FONn4R+BJwj3PuQGyfPwLfAf5sZkuBAcAteL8YRje8kHPufTP7cazGdWb2P3hXv5wFfAzMds6VmdlNeJcovmZmTwMHgJNi7/si7eg6kp7LPj1xLhIeZpaJ15+cBhQ656ra8ZwfAF/G6ybJAnbihed9zrna2D6L8cJ9NPAIcD5QDjxKsytPzOyfgO/jBe524F68k6E/dM41Hzz0deBf8bpYKoG/AT+K70oxswtirzcuVt/f8a6CWeic29zOj0Z6IAW5hJKZpeMd0T7jnPvnBL7uYmCqcy4nUa8pkmzqI5ew+greQKAlPtch4jv1kUuomNk5eN0e/wG87pxb63NJIr7TEbmEzU14fdf7ga/5XItIIKiPXEQk5HRELiIScgpyEZGQC/XJztjIuEF41/mKiHQXucDHrp1936EOcrwQ/8jvIkREkuBEvEFhxxT2IC8H2L17N3l5eX7XIiLSZWVlZQwePBg60NMQ9iAHIC8vT0EuIj2WTnaKiIScglxEJOQU5CIiIdct+sjb4pyjrq6O+vp6v0sRIBKJkJ6ezqdrKohIV3XrIK+pqWHPnj1UVlb6XYrE6d27NwMHDiQzM9PvUkS6hW4b5NFolO3btxOJRBg0aBCZmZk6CvSZc46amhoOHDjA9u3bOfnkk0lLU++eL5wD/X/oNrptkNfU1BCNRhk8eDC9e/f2uxyJ6dWrFxkZGezcuZOamhqys7P9LqnnidbD4xfCN/4CaZFj7y+B1+0Ph3TEFzz6nvhs50uwbwvs2uB3JZIg3faIXETiVByC6lLv3xsehmitd5s3yGvLyoc+/fyrT7pEQS7SE6y7HzY+BtkFEK3x2na+CIsug6oSOPubMGm+vzVKp+lvXJGeYNJ8uPJBr0+8KnZkXlXq3b9qgUI85HwNcjPbYWauhe1hP+sKgr1793LbbbcxYsQIsrOzOeGEEzjvvPN45JFHGi+nHDJkCGaGmdGrVy+GDBnCddddx+rVq32uXgLpjBtg6PimbcMmwpiZ/tQjCeP3EflZwMC47dJY+299q6gNqVoW78MPP2TMmDE899xz3HPPPbz++uts2LCB7373u6xcuZIXXnihcd958+axZ88etm3bxpIlS+jbty+XXHIJP/7xj1NSq4SIc153Su9+8Jkzvdsd6712CTVf+8idcwfi75vZ94EPgMCtjF4fdXx54Xr+cOv5RNKSe/3tzTffTHp6Oq+++ip9+vRpbB82bBhXX311k18oubm5FBYWAnDSSScxYcIEBg4cyJw5c5g6dSojR45Maq0SIvvegZoKuHYxDL8I3l8Fv/u61154mt/VSRf4fUTeyMwyga8C/9naqhhmlmVmeQ0b3ioaKbFxezFF+8rZtKM4qe9z6NAhnnvuOW655ZYmIR7vWAObbrvtNpxzrFixIhklSlgdfzLMetsLcYARF8Ost7x2CbXABDnwFaAvsLiNfWYDpXFbUlcHKq6oYeehCnYequCJ9R9SW+94Yv32xrbiipqEv+f777+Pc+6oI+njjz+enJwccnJy+N73vtfmaxx33HEMGDCAHTt2JLw+CbH0LMhuNm9/dr7XLqEWpMsP/xn4k3Pu4zb2mQ88EHc/lySG+cLV77Fkw07ye2VQWx8F4JUPDzH5Fy9ReqSWr537WeZcdWqy3r6JjRs3Eo1GmTlzJtXV1cfc3zmnKQlEeohABLmZfRa4BJjc1n7OuWqgMcWSHVRzrjqVzw3M4/5nt3Goog6Asqo6BuRmMX/yKK4dOzjh7zlixAjMjG3btjVpHzZsGOANcT+WQ4cOceDAAYYOHZrw+kQkeILStfJ1YD/wR78Lae66sYM5d3jTEW/njeiXlBAH6NevH5deeikLFy6koqKiU6/x0EMPkZaWxle+8pXEFicigeT7EbmZpeEF+ZPOuTq/62nOOccr24s5rk8GJx3Xm13Flbz8YXFSuy5+8YtfcN555zF27Fjmzp3L6NGjSUtLY9OmTWzdupUzzzyzcd/y8nL27t1LbW0t27dv59e//jWLFi1i/vz5jBgxIin1iUiw+B7keF0qJwH/6XchLdm6t5zK6joennkG40/uz1+LDnDr0tfYureczw9MzoLPw4cP5/XXX+eee+5h9uzZfPTRR2RlZXHKKafwb//2b9x8882N+86ZM4c5c+aQmZlJYWEh48aNY9WqVVx44YVJqU1EgsdSNcglGWKXIJaWlpaSl9c0VKuqqti+fTtDhw7t0lSp1XX11NRFyc3OaGwrq6olKz2NrHRNAdoZifreiHRHZWVl5OfnA+Q758ra85wgHJEHWlZ65KjAzosLdRERvwXlZKeIiHSSglxEJOQU5CIiIacgFxEJOQW5iEjIKchFREJOQS4iEnIKchGRkFOQi4iEnII8gG688UbMjJ/85CdN2pcvX560ibrmzp3L6aefnpTXFpHkUpB3RArnpcnOzubee++lpKQkZe8pIuGkIG+vaD08NtG7TYFLLrmEwsJC5s+f3+o+69evZ/z48fTq1YvBgwfz7W9/u3EO84ULF3LaaZ8uqNtwNP/oo482eY+77rqLxYsXc/fdd/Pmm29iZpgZixcvBmDXrl1cffXV5OTkkJeXx3XXXce+ffsaX6PhSP5Xv/oVQ4YMIT8/n+nTp1NeXp7gT0REWqMgb6+dL8G+LbBrQ0reLhKJcM899/Dzn/+cjz46ejW7Dz74gEmTJjFlyhT+9re/sWzZMtavX8+tt94KwMSJE9myZQsHDhwAYO3atRx//PGsWbMGgNraWjZs2MAFF1zAtGnTuOOOOzj11FPZs2cPe/bsYdq0aUSjUa6++mqKi4tZu3Ytzz//PB9++CHTpk07qpbly5ezcuVKVq5cydq1a4/qFhKR5FGQt6XiEBR/6G0bHoZorXfb0FZxKKlvf80113D66afzwx/+8KjH5s+fz8yZM5k1axYnn3wyX/jCF1iwYAFLliyhqqqK0047jeOOO461a9cCsGbNGu64447G+xs3bqS2tpYvfOEL9OrVi5ycHNLT0yksLKSwsJBevXqxatUq3nrrLZYuXcqZZ57JOeecw5IlS1i7di2bNm1qrCUajbJ48WJOO+00xo8fzw033MCqVauS+tmIyKcU5G1Zdz8sPAsWXQa7XvLadr7o3V94lvd4kt177708+eSTvPvuu03a33zzTRYvXkxOTk7jdvnllxONRtm+fTtmxoQJE1izZg2ffPIJW7Zs4eabb6a6upqtW7eydu1azjrrLHr37t3qe7/77rsMHjyYwYM/XdbulFNOoW/fvk3qGTJkCLm5uY33Bw4cyP79+xP4KYhfwrxeQU+iIG/LpPlw5YOQFoGqUq+tqtS7f9UC7/EkmzBhApdffjmzZ89u0n748GG+9a1v8cYbbzRub775Ju+99x7Dhw8H4IILLmDNmjWsW7eOMWPGkJeX1xjua9euZeLEiQmpMSOj6fzsZkY0Gk3Ia4t/XH09O6ZMxdWn5ryQdJ6C/FjOuAGGjm/aNmwijJmZshJ+8pOf8Mwzz7Bhw6f982eccQZbtmxhxIgRR22ZmZnAp/3kv/3tb7ngggsAL9xfeOEFXnzxxcY2gMzMTOqb/Yf9/Oc/z+7du9m9e3dj25YtW/jkk0845ZRTkvcFSyBUvrqZqqIiKjdv9rsUOQYF+bE453Wn9O4HnznTu92xPqWXIo4aNYqZM2eyYMGCxrbvfe97vPTSS9x666288cYbvPfee6xYsaLxZCfA6NGjKSgoYOnSpU2CfPny5VRXV3Peeec17jtkyBC2b9/OG2+8wcGDB6muruaSSy5pfO/XXnuNjRs38rWvfY2JEycyduzYlH39kjp1JSXU7NpFza5dFC9eDHV1FD+5pLGtTpfDBpKC/Fj2vQM1FTBlEXxjNUx+HGoOe+0pNG/evCbdFaNHj2bt2rUUFRUxfvx4xowZw5w5cxg0aFDjPmbG+PHjMTPOP//8xufl5eUxduxY+vTp07jvlClTmDRpEhdeeCH9+/fnqaeewsxYsWIFBQUFTJgwgUsuuYRhw4axbNmy1H3hklIHH3mUD674EjumX994JF65aRM7pl/PB1d8iYOPPHqMVxA/aPHlY6mr9rbsuNevKoX0bEjP6vzr9mBafDnYSn733xx46CHqY5euAqQP6E//WbPoO3myj5X1DJ1ZfFlH5MeSntU0xAGy8xXi0m0VTJ1Cn3HnNGnrfc44hXiApftdgIgEi3OOyo2biBQUkDF4MLW7d1O5aRPOuaTN9SNdoyNyEWmiuqiI6JEjDLr/Pob+ZhmD7vsp0YoKqouK/C5NWqEjchFpInPoUEasXkUkJweAnPPPZ8TqVViWuhODSkEuIk2kZWZCbCxCg0jcyF0JHnWtiIiEnO9BbmafMbNfm9khMztiZm+ZmUabiIi0k69dK2ZWALwI/AX4InAAOBnQ8DERkXbyu4/8e8Bu59zX49q2+1WMiEgY+d218mXgVTP7rZntN7PXzewbre1sZllmltewAToDIyI9nt9BPgy4CXgPuBx4BFhgZv/Yyv6zgdK47eilc5IoVdMZNCy+bGZkZmYyYsQI5s2bR11dXWMdjz/+OOeeey55eXnk5ORw6qmnctttt/H++++npEYRCQ6/gzwNeM05d6dz7nXn3GPA48C/tLL/fCA/bjsxNWWmfm7mSZMmsWfPHt577z3uuOMO5s6dy3333YdzjhkzZvDtb3+bK664gueee44tW7bwxBNPkJ2dzY9+9KOU1CciweF3H/keYEuztneBKS3t7JyrBqob7qdyuHD83Mx9zj476e+XlZVFYWEhADfddBO///3v+cMf/sDQoUN5+umnWbFiBV/+8pcb9z/ppJMYN26cVnQR6YH8DvIXgZHN2v4B2OlDLUepKykhGlsNPn5u5oxYwKbl5pJeUJCSWnr16sWhQ4d46qmnGDlyZJMQj6e5MER6Hr+7Vn4GjDOzO81shJnNAL4JPOxzXUAw5mZ2zvHCCy/w7LPPctFFF1FUVMTIkU1/982aNatx3c4TT0xZb5OIBISvQe6c2wRcA1wPvA38BzDLOfdfftbVoPDO2RTOnQuRCNEyb1rgaFkZlh5h4Ly7Kbxzdtsv0AUrV64kJyeH7OxsvvjFLzJt2jTmzp3b4r4/+MEPeOONN5gzZw6HDx9OWk0iEkx+d63gnFsJrPS7jtYUTJ1C5SsvU/bMpyWmYm7mCy+8kEceeYTMzEwGDRpEerr3rTr55JPZtm1bk3379+9P//79GTBgQFJrEpFg8rtrJfDi52bOHj2aSEFB49zMydSnTx9GjBjBSSed1BjiANdffz3btm1jxYoVSX1/EQkPBfkxBG1u5unTpzN16lSmT5/OvHnzeOWVV9ixYwdr165l2bJlRCIRX+oSEf/43rUSdEGbm9nMWLZsGY8//ji//OUv+elPf0ptbS0nnngiF198MQ888IAvdYmIf7T4sqRcIL43zoEu1ZQA0uLLIu0RrYfHJnq3It2Aglx6np0vwb4tsGuD35WIJIT6yKVnqDgE1aXevzc8DNFa7zZvkNeWlQ99+vlXn0gXKMilZ1h3P2x8DLILIFrjte18ERZdBlUlcPY3YdJ8f2sU6aRu37US5pO53ZUv35NJ8+HKByEtAlWxI/OqUu/+VQsU4hJq3TbIMzIyAKisrPS5Emmu4XvS8D1KmTNugKHjm7YNmwhjZqa2DpEE67ZdK5FIhL59+7J//34AevfurZkBfeaco7Kykv3799O3b9/UD15yzutO6d0PCoZAyQ7YsV6XIkroddsgBxrn824IcwmGvn37Nn5vUmrfO1BTAdcuhuEXwfur4Hdf99oLT0t9PSIJ0m0HBMWrr6+ntrY2dYVJqzIyMvybRqCu2tuy435WqkohPRvS/RmpK9JcZwYEdesj8gaRSERzkIgX1s0DOzvfn1pEEqjbnuwUEekpFOQiIiGnIBcRCTkFuYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhJyCXEQk5BTkIiIhpyAXEQk5BbmISMgpyEVEQk5BLiIScr4GuZnNNTPXbNvqZ00iImEThPnI3wEuibtf51chIiJhFIQgr3PO7fW7CBGRsApCH/nJZvaxmX1oZv9lZie1tqOZZZlZXsMG5KawThGRQPI7yF8BbgQmATcBQ4F1ZtZaQM8GSuO2j1JQo4hIoAVq8WUz6wvsBL7jnHuihcezgPhFF3OBj461+LKISFiEfvFl59wnZlYEjGjl8WqguuG+maWqNBGRwPK7a6UJM8sBhgN7/K5FRCQs/L6O/H4zm2hmQ8zsC8DvgXrgKT/rEhEJE7+7Vk7EC+1+wAFgPTDOOXfA16pERELE1yB3zk338/1FRLqDQPWRi4hIxynIRURCTkEuIhJyCnIRkZBTkIuIhJyCXEQk5BTkIiIhpyAXEQk5BbmISMgpyEVEQk5BLokToLntRXoSBbkkRrQeHpvo3XZTQVqERSSeglwSY+dLsG8L7NrgdyVJ4err2TFlKq6++/6ikvDyexpbCbOKQ1Bd6v17w8MQrfVu8wZ5bVn50Keff/UlUOWrm6kqKqJy82b6nH223+WINKEgl85bdz9sfAyyCyBa47XtfBEWXQZVJXD2N2HSfH9r7IK6khKi5eUAFC9eDHV1FD+5hIzCQgDScnNJLyjwsUIRT6AWX+4oM8sDSrX4so9e+xWs/hEc3vtpW04hXDwHxsz0r64E2HvPfEqWLiWSl4erqyNaVkZaXh6Wnk59WRkFM2ZQeOdsv8uUbqYziy+rj1y65owbYOj4pm3DJoY+xAEK75xN4dy5EIkQLfP+P0XLyrD0CAPn3a0Ql8BQkEvXOOd1p/TuB58507vdsb7bXIpYMHUKfcad06St9znj6Dt5sk8ViRxNQS5ds+8dqKmAKYvgG6th8uNQc9hr7wacc1Ru3ESkoIDs0aOJFBRQuWmTLkWUQNHJTuma40+GWW9DduwcxYiLYdZbkJ7tb10JUl1URPTIET7zswfIOe88Dq9fz99v/w7VRUVkjxzpd3kigE52irQpWlODq6khkpPT2FZfXo5lZZGWmeljZdJddeZkp47IRdqQlpkJzQI7kpvrUzUiLVMfuYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhJyCXEQk5AIT5Gb2fTNzZvag37WIiIRJIILczM4CvgX8ze9aRETCxvcgN7Mc4L+AbwAlPpcjIhI6vgc58DDwR+fcC8fa0cyyzCyvYQM0xE5Eejxfh+ib2XTgDOCsdj5lNvDD5FUkIhI+vh2Rm9lg4CFgpnOuqp1Pmw/kx20nJqk8EZHQ8POI/ExgAPCamTW0RYAJZnYrkOWca7JkuXOuGqhuuB/3PBGRHsvPIF8FjGrW9ktgK3Bv8xAXEZGW+Rbkzrly4O34NjOrAA45595u+VkiItJcEK5aERGRLgjUwhLOuQv8rkFEJGx0RC4ikiSpWkpTQS4ikgT1UceVP19PfTT5Ya4gFxFJgo3biynaV86mHcVJf69A9ZGLiIRZcUUN5VW1ADyx/kNq6x1PrN/OwPxsAHKzMziuT2ZbL9Eplqo+nGSIzbdSWlpaSl5ent/liEgPN++Zd1iyYSf5vTKorY9SVlVHXnY6GZE0So/U8rVzP8ucq05t8zXKysrIz88HyHfOlbXnfRPatWJm/8vMNJBHRHqkOVedyj2TRxFJM8qq6gAoq6ojkmbMnzzqmCHeWcnoI9e4eRHpsa4bO5hzh/dr0nbeiH5cO3Zw0t6zQ33kZvY/x9glHwhvX0134xxoPhqRlHLO8cr2Yo7rk8FJx/VmV3ElL39YjHMuafNDdfSI/CogGyhtZTuc0Oqk86L18NhE71ZEUmbr3nIqq+t4aPoYlt9yPg9OG0NFdR1b95Yn7T07dLLTzP4GPOSce6KVx08HNjvnIokp75j16GRna7avg19dA19bDkPO97sakR6juq6emrooudkZjW1lVbVkpaeRlX7saOzMyc6OXn64GW8hiBaDHG+K2V0dfE1JlIpDUF3q/XvDwxCt9W7zBnltWfnQp1/rzxeRLstKjxwV2HlxoZ4MHT0izwIizrnK5JXUfjoib+bPs2HjY5BdANEaqCqF7HxIy4SqEjj7mzBpvt9Vikgbkn75oXOu2jlXaWarzeyoJdfMrMDMVnfkNSWBJs2HKx+EtIgX4uDdpkXgqgUKcZFuqrOXH14A3Gpmy82sT1x7JjCxy1VJ551xAwwd37Rt2EQYM9OfekQk6bpyHfklQCHwspkNSUw50mXOwc4XoXc/+MyZ3u2O9V67iHRLXQnyPXhH328Bm8zsgkQUJF207x2oqYApi+Abq2Hy41Bz2GsXkW6pU3OtxIbhD3TO7Y/dvwu4C7gXuEuXH/qortrbsuM+j6pSSM+G9Cz/6hKRdknF5YcNmgxPcs79yMzeBZ7s5OtJoqRnHR3Y2fn+1CIiKdHZIB8KHIhvcM79t5ltBcZ2uSoREWm3TgW5c25nK+3vAOqMFRFJIa0QJCIScgpyEZGQU5CLiIScglx8F+blBkWCQEEuvnL19eyYMhVXr3nTRTpLQS6+qnx1M1VFRVRu3ux3KSKh1dnryEU6ra6khGi5t1pK8eLFUFdH8ZNLyCgsBCAtN5f0ggIfKxQJF1+D3MxuAm4ChsSa3gHmOef+5FtRknQHH3mUkqVLieTl4eq8lcYrN21ix/TrqS8ro2DGDArvnO1zlSLh4XfXykfA94Ez8UaErgZWmNmpvlYlSVV452wK586FSIRomTeVRLSsDEuPMHDe3QpxkQ7yNcidc8845/6fc+4951yRc+4HeAs4j/OzLkm+gqlT6DPunCZtvc8ZR9/Jk32qSCS8AtNHbmYR4FqgD7ChlX2ygPgZoXJTUJokgXOOyo2biBQUkDF4MLW7d1O5aRPOOczs2C8gIo387lrBzEaZ2WG8hZsfBa5xzm1pZffZQGnc9lFqqpREqy4qInrkCIPuv4+hv1nGoPt+SrSiguqiIr9LEwmdTs1HntACzDKBk4B8YCrwv4GJLYV5K0fkH2k+8vCJ1tTgamqI5OQ0ttWXl2NZWaRlZvpYmYi/UjkfecI452qA92N3N5vZWcBtwLda2Lca78gdQH+Ch1haZiY0C+xIrnrKRDrD966VFqTR9KhbRETa4Pd15POBPwG78LpJZgAXAJf7WJa/nAP9pSEiHeD3EfkAYAmwDVgFnAVc7px73teq/BKth8cmerciIu3k6xG5c+6f/Xz/wNn5EuzbArs2wJDz/a5GRELC95OdPV7FIagu9f694WGI1nq3eYO8tqx86NPPv/pEJPAU5H5bdz9sfAyyCyBa47XtfBEWXQZVJXD2N2HSfH9rFJE2+T2Qze8+cpk0H658ENIiUBU7Mq8q9e5ftUAhLhJw9VHHlT9fT33UvzE5CvIgOOMGGDq+aduwiTBmpj/1iEi7bdxeTNG+cjbtKPatBnWtBIFzXndK735QMARKdsCO9boUUSSgiitqKK+qBeCJ9R9SW+94Yv12BuZnA5CbncFxfVI3QllBHgT73oGaCrh2MQy/CN5fBb/7utdeeJrf1YlIMwtXv8eSDTvJ75VBbX0UgFc+PMTkX7xE6ZFavnbuZ5lzVepm4/Z9rpWuMLM8oDT0c63UVXtbdtzXUFUK6dmQrkGuIkH0m1d3c/+z29hf3jhrCANys/j3y0dy7djBnX7dzsy1oj7yIEjPahriANn5CnGRALtu7GDOHd700uDzRvTrUoh3lrpWREQ6wTnHK9uLOa5PBicd15tdxZW8/GGxL5ci6ohcRKQTtu4tp7K6joemj2H5Lefz4LQxVFTXsXVvecprUR+5iEgnVNfVU1MXJTc7o7GtrKqWrPQ0stIjnX7dUM5HLiISRlnpkaMCOy8u1FNJXSsiIiGnIBcRCTkFuYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhJyCPJlCPNhKRMJDQZ4s0Xp4bKJ3KyKSRAryZNn5EuzbArs2+F1Ju4V5ugaRnkxBnkgVh6D4Q2/b8DBEa73bhraKQ35X2CpXX8+OKVNx9foLQqQtQTzgUZAn0rr7YeFZsOgy2PWS17bzRe/+wrO8xwOq8tXNVBUVUbl5s9+liARWEBZabokmzUqkSfNhwCmw+kdQWeq1VZVCTiFctSBwiynXlZQQLfem3CxevBjq6ih+cgkZhYUApOXmkl5Q4GOFIsESv9DyuGH9jv2EFFGQJ9oZN8D2tfDWbz9tGzYxcCEOcPCRRylZupRIXh6urg6Ayk2b2DH9eurLyiiYMYPCO2f7XKWIv4K20HJLNB95ojkHPzvFW4OzYAiU7PDW3rz9HUjxqiHtUfK7/+bAQw9Rf+BAY1v6gP70nzWLvpMn+1iZSDDMe+adJgstl1XVkZedTkYkLSkLLWvNziDY9w7UVMCURfCN1TD5cag57LUHUMHUKfQZd06Ttt7njFOIi8TMuepU7pk8ikiaUVbl/eVaVlVHJM2YP3lUQkO8s3wNcjObbWabzKzczPab2XIzG+lnTV12/Mkw620YfpF3f8TFMOstrz2AnHNUbtxEpKCA7NGjiRQUULlpUyDPzIv4JUgLLbfE7yPyicDDwDjgUiADeM7M+vhaVVekZ0F2s26e7HyvPYCqi4qIHjnCoPvvY+hvljHovp8SraiguqjI79JEAiN+oeXTB+dzXJ+MxoWWg8DXk53OuUnx983sRmA/cCbwVz9q6mkyhw5lxOpVRHJyAMg5/3xGrF6FZQXzF4+IHxoWWn545hmMP7k/fy06wK1LX2Pr3nI+P9D/83OBOtlpZiOA94BRzrm3W3g8C4hPmFzgo0Cd7BSRbidZCy23JNQnO80sDXgQeLGlEI+ZDZTGbR+lpjoR6cmy0iNNQhy8hZYTHeKdFZggx+srPw2Y3sY+84H8uO3EFNQlIhJogRgQZGYLgSuBCc65Vo+ynXPVQHXc81JQnYhIsPka5OYl8c+Ba4ALnHPb/axHRKSBcy40B4t+d608DHwVmAGUm1lhbOvlc10i0oMFdXKs1vgd5Dfh9XWvAfbEbdN8rKn9AnTFj4gkTvzkWGHg93Xk4fi7pSXRenj8QvjGXyAtGGeuRaTzwjA5VmsCcbIzlOJXABpyvt/ViEgXLVz9XpPJsQBe+fAQk3/xUlImx0okBXlHVByC6tg84/ErAOUN8tqy8qFPcOYoFpH2m3PVqXxuYB73P7uNQxWfTo41IDeL+ZNHBWZelZYoyDti3f2w8THILoBojdfWsAJQVQmc/U1vcQkRCaXrxg7mxfcPsuKNjxvbgjQ5Vmv8PtkZLpPmw5UPen3iVXErAKVFvBWAFOIioRb0ybFaoyDvqDNugKHjm7YFdAUgEemYhsmxHpo+huW3nM+D08ZQUV3H1r3lfpfWpkBNmtVRvqwQFLIVgESk/VI5OVZrQj1pVmiEbAUgETm2hgPaoE+O1RodkXdUXbW3xS8eUVXqHZWnePGIMA0hFgmq+qjjywvX84dbzyeS5v//Jx2Rp0JAVgBy9fXsmDIVV1+f0vcV6W7CNoqzJbr8MKQqX91MVVERlZs30+fss/0uRyRUwjyKsyUK8hCpKykhWu6dPS9evBjq6ih+cgkZhYUApOXmkl5Q4GOFIuEQ5lGcLVEfeYjsvWc+JUuXEsnLw9XVES0rIy0vD0tPp76sjIIZMyi8c7bfZYqEwm9e3c39z25jf3njEgcMyM3i3y8f6esAIPWRJ0tAftkV3jmbwrlzIRIhWuZ9f6NlZVh6hIHz7laIi3TAdWMHc+7wplNqhGEUZ0sU5McSrYfHJnq3AVAwdQp9xp3TpK33OePoO3myTxWJhEtDL0RYR3G2REF+LPGzHAaAc47KjZuIFBSQPXo0kYICKjdtCuUPn0iqxS8YEdZRnC3Ryc6WBHiWw+qiIqJHjvCZnz1AznnncXj9ev5++3eoLioie+RIX2oSCYv4Sw3HnNSXF79/UeMAoAn/0J/137+IrPTwHd/qZGdL/jy76SyHVaXeteJpmb7PchitqcHV1BDJyWlsqy8vx7KySMsMz+VSIqkSf6nh/1m5hRfe3c+lp5zAXV/6PBC8Sw07c7JTQd6a134Fq38Eh/d+2pZTCBfP0QRZIiEy75l3mlxqWFZVR152OhmRtEBeaqirVhJJsxyKhJ5zjjlXnco9k0cRSTPKqj5dMCKSZsyfPCpQId5ZCvLWOOctGtG7H3zmTO92x/rAXIooIm2LP7HZnS41bImCvDWa5VAk1OJPbHanSw1boqtWWnP8yTDr7U8nyBpxMcx6y5vlUEQCqbU5VCqq6zh8pJafXjuaK0YN4q9FB7h16Wts3VvO5weGf1S4TnaKSLfR2onN9EgaZc1ObKZ6wYj20snORArxLziRnqitE5vpLZzYDMOCEe2lIG9JCoflh/kvIpGg6EknNluiIG9Jiobla3EIkcToSSc2W6KTnQ18GJavxSFEOsc5R0llbY88sdkSX4PczCYA/w6cCQwErnHOLfelmHX3Nx2WD9515IsuS+iwfC0OIdI1DWtsnjWkgF+/vOuoxSE27yyhqi7KqztKuGLUoFDPodJefn9lfYA3gVt8rsML6SsfhLSIN7cKeLdpEbhqQcLmVjn4yKN8cMWX2DH9eio3bwagctMmdky/ng+u+BIHH3k0Ie8j0l01dKNMOm1gjzyx2RJfg9w59yfn3F3Oud/7WUejFAzL1+IQIh1XXFHDzkMV7DxU0aQb5ZyhxzH6xPwm+3b3E5stCVUfuZllAfHL1ecm9A3ih+UXDIGSHZ8OyzdL2NsUTJ1C5SsvU/bMysY2LQ4h0rq21tg8VFFDVrrx+YF57CqubDyxaQn8Pxt0fnetdNRsoDRu+yihr56iYflaHEKkY1q7PtwB2enGon88K/SLQ3RFYEZ2mpnjGCc7Wzki/yhhIzvrqr0tO+61qkq9YfnpWa0/r4Oqtm1j5w1fO2pxiM/++ldaHEKkDbc9/Tor3vi48f6X/9dAfnzNqMbFISC4IzbbqzMjO0PVteKcqwYal7xO+J9O6VlHB3Z2fsv7dkHm0KGMWL2qcXGInPPPZ8TqVVhW4n5ZiHQ38deHn3Rcb3YVV7JpRwk5WU1jLC8u1HuKsHWtdAtpmZlNVvgBiOTmaoUfkTZ0pzU2E83XrhUzywFGxO6+DnwH+AtQ7Jzb1Y7na9IskR6iuq6emrpot+pGaUkYu1bG4gV3gwdit08CN6a8mgTqaWfNRZItKz1yVGD3xG6Ulvh9Hfka55y1sN3oZ11dpTlURCSV1EeeBPFzqIiIJJvfXSvdhuZQERG/KMgT5OAjj1KydCmRvDxcnTdgoWEOlfqyMgpmzNDwexFJCnWtJIjmUBERvyjIE6hg6hT6jDunSZvmUBGRZFPXSgLFz6GSMXgwtbt3N86hoksRRSRZdESeQNVFRUSPHGHQ/fcx9DfLGHTfT4lWVFBdVOR3aSLSjQVm0qzO8HtkZ/Mj7WhNDa6mpsnw+/ryciwrS8PvRaRdOjOyU0fkndTSoB/NoSIiflCQd5IG/YhIUOhkZwdo0I+IBJGCvAM06EdEgkhdKx2gQT8iEkQK8nZquLpHg35EJGgU5O0Qf4WKFk4WkaBRkLdD/BUqGvQjIkGjAUGtiL9CZd/8n3D4L38h5+KL6f+d23G1tWSccELjFSoa9CMiidKZAUEK8lbs+fE9fPLUU41XqETLykjLy8PS03WFiogkjUZ2Joirr6dq82ZOmDNHV6iISOApyFvQ0CeeNeSzukJFRAJPA4JiWh61+SRH3niTtLw8Mk48kbo9ezQtrYgEjoI8psVRmy+/QrSiAtLS6D12LDkTxvP3279DdVER2SNH+lyxiIhHJzvjlPzuvznw0EPUHzjQ2BY5/ngGfOf2xu4UXaEiIsmkk51d1NKozT7nntukT1zT0opI0KhrJY6WahORMNIReRyN2hSRMFIfeRwt1SYifutMH7m6VuKkZWZCs8CO5Ob6VI2ISPsEomvFzG4xsx1mVmVmr5jZ2cl+zzlz5iT7LUREUsL3IDezacADwN3AGcCbwLNmNiBZ73ngwAGeeOIJDh48mKy3EBFJGd+DHPgO8Lhz7pfOuS3AvwCVwD8l6w0ff/xxPv74YxYtWpSstxARSRlfg9zMMoEzgRca2pxz0dj9c1vYP8vM8ho2oN0d2EuXLmX48OEMHz6cBQsWAPDggw82ti1durSrX46IiC/8PiI/HogA+5q17wMKW9h/NlAat33U3jeaMWMGTz/9NMcffzz79nlvt2/fPvr378+yZcuYMWNGZ+oXEfGd30HeUfOB/LjtxI48+ayzzuK6665r0nbttdcyduzYhBUoIpJqfgf5QaAeOKFZ+wnA3uY7O+eqnXNlDRtQ3tE3XLduHcOHD+df//VfGT58OOvWretU4SIiQeFrkDvnaoDNwMUNbWaWFru/IdHvV1lZycGDB3n++edZsGABzz//PAcPHuTIkSOJfisRkZTxfWRn7PLDJ4FvARuBWcB1wOecc837zps/t0MjOxu+1vh5U1pqExHxSyhHdjrnlplZf2Ae3gnON4BJxwrxzmgprBXgIhJ2vgc5gHNuIbDQ7zpERMLI75OdIiLSRQpyEZGQU5CLiIScglxEJOQU5CIiIacgFxEJuUBcfthVZWXtumZeRCTwOpNnvo/s7Aoz+wwdmAFRRCRETnTO/b09O4Y9yA0YRMcnz8rF+wVwYieeGwSq31+q3189of5c4GPXzoAOdddK7Its12+seHHD8svbO5dBkKh+f6l+f/WQ+jv0delkp4hIyCnIRURCrqcGeTVwd+w2jFS/v1S/v1R/M6E+2SkiIj33iFxEpNtQkIuIhJyCXEQk5BTkIiIh122D3MxuMbMdZlZlZq+Y2dnH2P9aM9sa2/8tM7siVbW2Uk+76zezG83MNduqUllvs3ommNkzZvZxrJavtOM5F5jZa2ZWbWbvm9mNya+01Vo6VH+s9uafvzOzwhSVHF/LbDPbZGblZrbfzJab2ch2PC8QP/+dqT9IP/9mdpOZ/c3MymLbBjP74jGe0+XPvlsGuZlNAx7Au8TnDOBN4FkzG9DK/l8AngKeAMYAy4HlZnZaSgo+up4O1R9TBgyM2z6b7Drb0Aev5lvas7OZDQX+CPwFOB14EFhkZpcnqb5j6VD9cUbS9HuwP8F1tcdE4GFgHHApkAE8Z2Z9WntCwH7+O1x/TFB+/j8Cvg+cCYwFVgMrzOzUlnZO2GfvnOt2G/AKsDDufhreUP7vt7L/MmBls7aXgUdDUv+NwCd+f+6t1OaArxxjn3uBt5u1PQ38OST1XxDbr6/f9bZQW/9YbRPa2CdQP/+dqD+wP/+x+oqBf07mZ9/tjsjNLBPvt+ELDW3OuWjs/rmtPO3c+P1jnm1j/6TpZP0AOWa208x2m1mrRwABFZjPv4veMLM9Zva8mZ3ndzEx+bHb4jb2CfLn3576IYA//2YWMbPpeH/hbWhlt4R89t0uyIHjgQiwr1n7PqC1PsvCDu6fTJ2pfxvwT8DVwFfxvq8vmdmJySoywVr7/PPMrJcP9XTUHuBfgCmxbTewxszO8LMoM0vD66Z60Tn3dhu7Bunnv1EH6g/Uz7+ZjTKzw3gjNx8FrnHObWll94R89qGe/VA8zrkNxP3GN7OXgHeBbwH/4VddPYVzbhtemDR4ycyGA7cDN/hTFeD1NZ8GnO9jDV3RrvoD+PO/De9cTz4wFXjSzCa2EeZd1h2PyA8C9cAJzdpPAPa28py9Hdw/mTpTfxPOuVrgdWBEYktLmtY+/zLn3BEf6kmEjfj4+ZvZQuBK4ELn3LEWXwnSzz/Q4fqb8Pvn3zlX45x73zm32Tk3G+/E+W2t7J6Qz77bBblzrgbYDFzc0Bb7E+1iWu+n2hC/f8ylbeyfNJ2svwkziwCj8P7kD4PAfP4JdDo+fP7mWQhcA1zknNvejqcF5vPvZP3NXyNoP/9pQFYrjyXms/f7jG6SzhJPA6qAfwQ+D/xfoAQ4Ifb4EmB+3P5fAGqBO4DPAXOBGuC0kNQ/B7gMGIZ3ueJTwBHgFJ/qz8ELstPxrji4Pfbvk2KPzweWxO0/FKgAfhr7/G8G6oDLQ1L/LLz+2RF4XQEP4v1VdbEPtf8C+ATvMr7CuK1X3D6B/fnvZP2B+fmP/WxMAIbg/TKZD0SBS5P52af8P0kKP9BbgZ14JxxeAc6Je2wNsLjZ/tfi9W1VA28DV4SlfuBncfvuxbsme4yPtV8QC8Dm2+LY44uBNS085/XY1/ABcGNY6ge+C7wfC49DeNfDX+hT7S3V7eI/zyD//Hem/iD9/ONdD74jVst+vCtSLk32Z69pbEVEQq7b9ZGLiPQ0CnIRkZBTkIuIhJyCXEQk5BTkIiIhpyAXEQk5BbmISMgpyEVEQk5BLiIScgpykQ4ws9Fmti62vuJuM/uu3zWJKMhF2snM8oDn8Ob1OBP4d2CumX3T18Kkx9NcKyJxzGwI0NLUqWvx1lf8MVDovOmGMbOf4K3p+bmUFSnSjI7IRZraTdPV2MfgzWj4V7x1FP/aEOIxzwIjzawg1YWKNFCQi8RxztU75/Y65/bizYv9KN4k/3NpfX1F8Hl9S+nZtGanSOv+E8jFm086amZ+1yPSIgW5SAvM7C7gcuBs51x5rLm19RUbHhPxhbpWRJoxsyl4y4dd55z7IO6hDcAEM8uIa7sU2OacK0lljSLxdNWKSBwzOw1vab0HgIfjHqrBW4dzG94liPfirc/5n8DtzrnHUlyqSCMFuUgcM7sR+GULD611zl1gZqPxAv4s4CDwc+fcvSksUeQoCnIRkZBTH7mISMgpyEVEQk5BLiIScgpyEZGQU5CLiIScglxEJOQU5CIiIacgFxEJOQW5iEjIKchFREJOQS4iEnL/H8VE4PHmwjdSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
    "axes.set_title('y space')\n",
    "axes.scatter(historyGDy[:,0], historyGDy[:,1], lw=0.5, marker='*', color='#1F77B4', label='GD')\n",
    "axes.scatter(historyNty[:,0], historyNty[:,1], lw=0.5, marker='*', color='#FF7F0E', label='Newton')\n",
    "axes.scatter(historyPGy[:,0], historyPGy[:,1], lw=0.5, marker='*', color='#D62728', label='PG')\n",
    "axes.scatter([0], [0], lw=0.25, color='black', marker='*') \n",
    "axes.set_xlabel('z0'); axes.set_ylabel('z1'); axes.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These trajectories confirm the intuition outlined in the previous sections: GD in blue gives a very sub-optimal trajectory in $\\mathbf{y}$. Newton (in orange) does better, but is still clearly curved. It can't approximate the higher order terms of this example well enough. This is in contrast to the straight, and diagonal red trajectory for the optimization using the inverse simulator.\n",
    "\n",
    "The behavior in intermediate spaces becomes especially important when they're not only abstract latent spaces as in this example, but when they have actual physical meanings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions \n",
    "\n",
    "Despite its simplicity, this example already shows surprisingly large differences between gradient descent, Newton's method, and using the _inverse simulator_.\n",
    "\n",
    "The main takeaways of this section are the following.\n",
    "* GD easily yields \"unbalanced\" updates, and gets stuck.\n",
    "* Newtons method does better, but is far from optimal.\n",
    "* the higher-order information of the inverse simulator  outperform both, even if it is applied only partially (we still used Newton's method for $L$ above).\n",
    "* Also, the methods (and in general the choice of optimizer) strongly affects progress in latent spaces, as shown for $\\mathbf{y}$ above.\n",
    "    \n",
    "In the next sections we can build on these observations to use PGs for training NNs via invertible physical models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Approximate inversions\n",
    "\n",
    "If an analytic inverse like the `physics_y_inv_analytic` above is not readily available, we can actually resort to optimization schemes like Newton's method or BFGS to obtain a local inverse numerically. This is a topic that is orthogonal to the comparison of different optimization methods, but it can be easily illustrated based on the inverse simulator variant from above.\n",
    "\n",
    "Below, we'll use the BFGS variant `fmin_l_bfgs_b` from `scipy` to compute the inverse. It's not very complicated, but we'll use numpy and scipy directly here, which makes the code a bit messier than it should be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BFGS optimization test run, find x such that y=[2,2]:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([2.00000003, 1.41421353])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def physics_y_inv_opt(target_y, x_ini):\n",
    "    # a bit ugly, we switch to pure scipy here inside each iteration for BFGS\n",
    "    import numpy as np\n",
    "    from scipy.optimize import fmin_l_bfgs_b\n",
    "    target_y = onp.array(target_y)\n",
    "    x_ini    = onp.array(x_ini)\n",
    "\n",
    "    def physics_y_opt(x,target_y=[2,2]):\n",
    "        y = onp.array( [x[0], x[1]*x[1]] ) # we cant use physics_y from JAX here\n",
    "        ret = onp.sum( onp.square(y-target_y) )\n",
    "        return ret\n",
    "    \n",
    "    ret = fmin_l_bfgs_b(lambda x: physics_y_opt(x,target_y), x_ini, approx_grad=True )\n",
    "    #print( ret ) # return full BFGS details\n",
    "    return ret[0]\n",
    "\n",
    "print(\"BFGS optimization test run, find x such that y=[2,2]:\")\n",
    "physics_y_inv_opt([2,2], [3,3])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonetheless, we can now use this numerically inverted $\\mathbf{y}$ function to perform the inverse simulator optimization. Apart from calling `physics_y_inv_opt`, the rest of the code is unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PG iter 0: [2.09999967 2.50998022]\n",
      "PG iter 1: [1.46999859 2.10000011]\n",
      "PG iter 2: [1.02899871 1.75698602]\n",
      "PG iter 3: [0.72029824 1.4699998 ]\n",
      "PG iter 4: [0.50420733 1.22988982]\n",
      "PG iter 5: [0.35294448 1.02899957]\n",
      "PG iter 6: [0.24705997 0.86092355]\n",
      "PG iter 7: [0.17294205 0.72030026]\n",
      "PG iter 8: [0.12106103 0.60264817]\n",
      "PG iter 9: [0.08474171 0.50421247]\n"
     ]
    }
   ],
   "source": [
    "x = np.asarray([3.,3.])\n",
    "eta = 0.3\n",
    "history = [x]; updates = []\n",
    "\n",
    "Gy = jax.grad(loss_y)\n",
    "Hy = jax.jacobian(jax.jacobian(loss_y))\n",
    "for i in range(10):    \n",
    "    # same as before, Newton step for L(y)\n",
    "    y = physics_y(x)\n",
    "    g = Gy(y)\n",
    "    y += -eta * np.matmul( np.linalg.inv( Hy(y) ) , g)\n",
    "\n",
    "    # optimize for inverse physics, assuming we dont have access to an inverse for physics_y\n",
    "    x = physics_y_inv_opt(y,x)\n",
    "    history.append(x)\n",
    "    updates.append( history[-2] - history[-1] )\n",
    "    print( \"PG iter %d: \"%i + format(x) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This confirms that the approximate inversion works, in line with the regular PG version above. There's not much point plotting this, as it's basically the same, but let's measure the difference. Below, we compute the MAE, which for this simple example turns out to be on the order of our floating point accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE difference between analytic PG and approximate inversion: 0.000001\n"
     ]
    }
   ],
   "source": [
    "historyPGa = onp.asarray(history)\n",
    "updatesPGa = onp.asarray(updates) \n",
    "\n",
    "print(\"MAE difference between analytic PG and approximate inversion: %f\" % (np.average(np.abs(historyPGa-historyPG))) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Next steps\n",
    "\n",
    "Based on this code example you can try the following modifications:\n",
    "\n",
    "- Instead of the simple L(y(x)) function above, try other, more complicated functions.\n",
    "\n",
    "- Replace the simple \"regular\" gradient descent with another optimizer, e.g., commonly used DL optimizers such as AdaGrad, RmsProp or Adam. Compare the existing versions above with the new trajectories."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}