physgrad-hig-code.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2tT-NooLyyGW"
      },
      "source": [
        "# Coupled Oscillators with Half-Inverse Gradients\n",
        "\n",
        "In this notebook, we'll turn to a practical example comparing the half-inverse gradients (HIGs) to other methods for training neural networks with physical loss functions. Specifically, we'll compare:\n",
        "\n",
        "1. Adam: as a standard gradient-descent (GD) based network optimizer,\n",
        "2. Scale-Invariant Physics: the previously described algorithm that fully inverts the physics,\n",
        "3. Half-Inverse Gradients: which locally and jointly inverting physics and network.\n",
        "\n",
        "\n",
        "\n",
        "## Inverse problem setup\n",
        "\n",
        "The learning task is to find the control function steering the dynamics of a coupled oscillator system. \n",
        "This is a classical problem in physics, and a good case to evaluate the HIGs due to it's smaller size. We're using two mass points, and thus we'll only have four degrees of freedom for position and velocity of both points (compared to, e.g., the $32\\times32\\times2$ unknowns we'd get even for \"only\" a small fluid simulation with 32 cells along x and y). \n",
        "\n",
        "Nonetheless, the oscillators are a highly-non trivial case: we aim for applying a control such that the initial state is reached again after a chosen time interval. We'll use 24 steps of a fourth-order Runge-Kutta scheme, and hence the NN has to learn how to best \"nudge\" the two mass points over the course of all time steps, so that they end up at the desired position with the right velocity at the right time.\n",
        "\n",
        "A system of $N$ coupled oscillators is described by the following Hamiltonian:\n",
        "\n",
        "$$\n",
        "  \\mathcal{H}(x_i,p_i,t)=\\sum_i \\bigg( \\frac{x_i^2}{2}+ \\frac{p_i^2}{2} +  \\alpha \\cdot (x_i-x_{i+1})^4+u(t) \\cdot x_i \\cdot c_i\\bigg),\n",
        "$$\n",
        "\n",
        "which provides the basis for the RK4 time integration below.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5NOyhZ-MOnoZ"
      },
      "source": [
        "\n",
        "## Problem statement\n",
        "\n",
        "More concretely, we consider a set of different phyiscal inputs ($x_i$). Using a corresponding control function ($u_i$), we can influence the time evolution of our physical system $\\mathcal P$ and receive an output state ($y_i$).\n",
        "\n",
        "$$y_i=\\mathcal P(x_i;u_i)$$\n",
        "\n",
        "If we want to evolve a given initial state($x_i$) into a given target state ($y^*_i$), we receive an inverse problem for the control function $(u_i)$. The quality of between desired target state ($y^*_i$) and received target state ($y_i$) is measured by a loss function ($L$).\n",
        "\n",
        "$$\\text{arg min}_{u_i}\\; L(y^*_i,\\mathcal P(x_i;u_i))$$\n",
        "\n",
        "If we use a neural network ($f$ parameterized by $\\theta$) to learn the control function over a set of input/output pairs ($x_i,y^*_i$), we transform the above physics optimization task into a learning problem.\n",
        "\n",
        "$$\\text{arg min}_{\\theta}\\; \\sum_i L\\big(y^*_i, \\mathcal P(x_i;f(x_i,y_i; \\theta)) \\big)$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "We6eG_ZWOspS"
      },
      "source": [
        "Before we begin by setting up the physics solver $\\mathcal P$, we import the necessary libraries (this example uses TensorFlow), and define the main global variable `MODE`, that switches between _Adam_ (`'GD'`), _Scale-invariant physics_ (`'SIP'`), and _Half-inverse gradients_ (`'HIG'`)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 95,
      "metadata": {
        "id": "A0cyHl7tyyGg"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import tensorflow as tf\n",
        "import time, os\n",
        "\n",
        "# main switch for the three methods:\n",
        "MODE = 'HIG'  # HIG | SIP | GD"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yCN4F_AjyyGk"
      },
      "source": [
        "## Coupled linear oscillator simulation\n",
        "\n",
        "For the physics simulation, we'll solve a differential equation for a system of coupled linear oscillators with a control term. The time integration, a fourth-order Runge-Kutta scheme is used.\n",
        "\n",
        "Below, we're first defining a few global constants:\n",
        "`Nx`: the number of oscillators,\n",
        "`Nt`: the number of time evolution steps, and\n",
        "`DT`: the length of one time step. \n",
        "\n",
        "We'll then define a helper function to set up a Laplacian stencil, the `coupled_oscillators_batch()` function which computes the simulation for a whole mini batch of values, and finally the `solver()` function, which runs the desired number of time steps with a given control signal."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 96,
      "metadata": {
        "id": "hevmSaVvyyGl"
      },
      "outputs": [],
      "source": [
        "Nx = 2\n",
        "Nt = 24\n",
        "DT = 0.5\n",
        "\n",
        "def build_laplace(n,boundary='0'):\n",
        "    if n==1:\n",
        "        return np.zeros((1,1),dtype=np.float32)\n",
        "    d1 = -2 * np.ones((n,),dtype=np.float32)\n",
        "    d2 = 1 * np.ones((n-1,),dtype=np.float32)\n",
        "    lap = np.zeros((n,n),dtype=np.float32)\n",
        "    lap[range(n),range(n)]=d1\n",
        "    lap[range(1,n),range(n-1)]=d2\n",
        "    lap[range(n-1),range(1,n)]=d2\n",
        "    if boundary=='0':\n",
        "        lap[0,0]=lap[n-1,n-1]=-1\n",
        "\n",
        "    return lap\n",
        "\n",
        "@tf.function\n",
        "def coupled_oscillators_batch( x, control):\n",
        "    '''\n",
        "    ODE of type:    x' = f(x)\n",
        "    :param x_in:    position and velocities, shape: (batch, 2 * number of osc) order second index: x_i , v_i\n",
        "    :param control: control function, shape: (batch,)\n",
        "    :return:\n",
        "    '''\n",
        "    #print('coupled_oscillators_batch')\n",
        "    n_osc = x.shape[1]//2\n",
        "\n",
        "    # natural time evo\n",
        "    a1 = np.array([[0,1],[-1,0]],dtype=np.float32)\n",
        "    a2 = np.eye(n_osc,dtype=np.float32)\n",
        "    A = np.kron(a1,a2)\n",
        "    x_dot1 = tf.tensordot(x,A,axes = (1,1))\n",
        "\n",
        "    # interaction term\n",
        "    interaction_strength = 0.2\n",
        "    b1 = np.array([[0,0],[1,0]],dtype=np.float32)\n",
        "    b2 = build_laplace(n_osc)\n",
        "    B = interaction_strength * np.kron(b1,b2)\n",
        "    x_dot2 = tf.tensordot(x,B, axes=(1, 1))\n",
        "\n",
        "    # control term\n",
        "    control_vector = np.zeros((n_osc,),dtype=np.float32)\n",
        "    control_vector[-1] = 1.0\n",
        "    c1 = np.array([0,1],dtype=np.float32)\n",
        "    c2 = control_vector\n",
        "    C = np.kron(c1,c2)\n",
        "    x_dot3 = tf.tensordot(control,C, axes=0)\n",
        "\n",
        "    #all terms\n",
        "    x_dot = x_dot1 + x_dot2 +x_dot3\n",
        "    return x_dot\n",
        "\n",
        "@tf.function\n",
        "def runge_kutta_4_batch(x_0, dt, control, ODE_f_batch):\n",
        "\n",
        "    f_0_0 = ODE_f_batch(x_0, control)\n",
        "    x_14 = x_0 + 0.5 * dt * f_0_0\n",
        "\n",
        "    f_12_14 = ODE_f_batch(x_14, control)\n",
        "    x_12 = x_0 + 0.5 * dt * f_12_14\n",
        "\n",
        "    f_12_12 = ODE_f_batch(x_12, control)\n",
        "    x_34 = x_0 + dt * f_12_12\n",
        "\n",
        "    terms = f_0_0 + 2 * f_12_14 + 2 * f_12_12 + ODE_f_batch(x_34, control)\n",
        "    x1 = x_0 + dt * terms / 6\n",
        "\n",
        "    return x1\n",
        "\n",
        "@tf.function\n",
        "def solver(x0, control):\n",
        "    x = x0\n",
        "    for i in range(Nt):\n",
        "        x = runge_kutta_4_batch(x, DT, control[:,i], coupled_oscillators_batch)\n",
        "    return x\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VtY4-wbXyyGn"
      },
      "source": [
        "## Training setup\n",
        "\n",
        "The neural network itself is quite simple: it consists of four dense layers (the intermediate ones with 20 neurons each), and `tanh` activation functions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 97,
      "metadata": {
        "id": "et1j8lYRyyGp"
      },
      "outputs": [],
      "source": [
        "act = tf.keras.activations.tanh\n",
        "model = tf.keras.models.Sequential([\n",
        "    tf.keras.layers.InputLayer(input_shape=(2*Nx)),\n",
        "    tf.keras.layers.Dense(20, activation=act),\n",
        "    tf.keras.layers.Dense(20, activation=act),\n",
        "    tf.keras.layers.Dense(20, activation=act),\n",
        "    tf.keras.layers.Dense(Nt, activation='linear')\n",
        "    ])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zilE67gYyyGq"
      },
      "source": [
        "As loss function, we'll use an $L^2$ loss:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 98,
      "metadata": {
        "id": "km83dkWgyyGq"
      },
      "outputs": [],
      "source": [
        "@tf.function\n",
        "def loss_function(a,b):\n",
        "    diff = a-b\n",
        "    loss_batch = tf.reduce_sum(diff**2,axis=1)\n",
        "    loss = tf.reduce_sum(loss_batch)\n",
        "    return loss"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZyS31DHdyyGr"
      },
      "source": [
        "And as data set for training we simply create 4k of random position values which the oscillators start with (`X_TRAIN`), and which they should return to at the end of the simulation (`Y_TRAIN`). As they should return to their initial states, we have `X_TRAIN=Y_TRAIN`.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 99,
      "metadata": {
        "id": "_zz7VSIpyyGs"
      },
      "outputs": [],
      "source": [
        "N = 2**12\n",
        "X_TRAIN = np.random.rand(N, 2 * Nx).astype(np.float32)\n",
        "Y_TRAIN = X_TRAIN # the target states are identical"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1KChFddqyyGs"
      },
      "source": [
        "## Training\n",
        "For the optimization procedure of the neural network training, we need to set up some global parameters. The next cell initializes some suitable values tailored to each of the three methods. These were determined heuristically to work best for each. If we try to use the same settings for all, this would inevitably make the comparison unfair for some of them.\n",
        "\n",
        "1. _Adam_: This is the most widely used NN optimizer, and we're using it here as a representative of the GD family. Note that the truncation parameter has no meaning for Adam.\n",
        "\n",
        "2. _SIP_: The specified optimizer is the one used for network optimization. The physics inversion is done via Gauss-Newton and corresponds to an exact inversion since the physical optimization landscape is quadratic. For the Jacobian inversion in Gauss-Newton, we can specify a truncation parameter.\n",
        "\n",
        "3. _HIG_: To obtain the HIG algorithm, the optimizer has to be set to SGD. For the Jacobian half-inversion, we can specify a truncation parameter. Optimal batch sizes are typically lower than for the other two, and a learning rate of 1 typically works very well.\n",
        "\n",
        "The maximal training time in seconds is set via `MAX_TIME` below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 100,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PYjx7Sk8yyGz",
        "outputId": "9ee462ad-499c-4f7c-d3a3-0fb11ab36abe"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Running variant: HIG\n"
          ]
        }
      ],
      "source": [
        "if MODE=='HIG': # HIG training\n",
        "  OPTIMIZER = tf.keras.optimizers.SGD\n",
        "  BATCH_SIZE = 32 # larger batches make HIGs unnecessariy slow...\n",
        "  LR = 1.0\n",
        "  TRUNC = 10**-10\n",
        "\n",
        "elif MODE=='SIP': # SIP training\n",
        "  OPTIMIZER = tf.keras.optimizers.Adam\n",
        "  BATCH_SIZE = 256\n",
        "  LR = 0.001  # for the internal step with Adam\n",
        "  TRUNC = 0   # not used\n",
        "\n",
        "else: #  Adam Training (as GD representative)\n",
        "  MODE = 'GD' \n",
        "  OPTIMIZER = tf.keras.optimizers.Adam\n",
        "  BATCH_SIZE = 256\n",
        "  LR = 0.001\n",
        "  TRUNC = 0 # not used\n",
        "\n",
        "# global parameters for all three methods\n",
        "MAX_TIME = 100 # [s]\n",
        "print(\"Running variant: \"+MODE)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4VFJ4pmEyyG0"
      },
      "source": [
        "The next function, `HIG_pinv()`, is a crucial one: it constructs the half-inverse of a given matrix for HIGs. It computes an SVD, takes the square-root of the singular values, and then re-assembles the matrix. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 101,
      "metadata": {
        "id": "c3B4qtUPyyG1"
      },
      "outputs": [],
      "source": [
        "from tensorflow.python.framework import ops\n",
        "from tensorflow.python.framework import tensor_shape\n",
        "from tensorflow.python.ops import array_ops\n",
        "from tensorflow.python.ops import math_ops\n",
        "from tensorflow.python.util import dispatch\n",
        "from tensorflow.python.util.tf_export import tf_export\n",
        "from tensorflow.python.ops.linalg.linalg_impl import _maybe_validate_matrix,svd\n",
        "\n",
        "# partial inversion of the Jacobian via SVD:\n",
        "# this function is adopted from tensorflow's SVD function, and published here under its license: Apache-2.0 License , https://github.com/tensorflow/tensorflow/blob/master/LICENSE\n",
        "@tf_export('linalg.HIG_pinv')\n",
        "@dispatch.add_dispatch_support\n",
        "def HIG_pinv(a, rcond=None,beta=0.5, validate_args=False, name=None):\n",
        "\n",
        "  with ops.name_scope(name or 'pinv'):\n",
        "    a = ops.convert_to_tensor(a, name='a')\n",
        "\n",
        "    assertions = _maybe_validate_matrix(a, validate_args)\n",
        "    if assertions:\n",
        "      with ops.control_dependencies(assertions):\n",
        "        a = array_ops.identity(a)\n",
        "\n",
        "    dtype = a.dtype.as_numpy_dtype\n",
        "\n",
        "    if rcond is None:\n",
        "\n",
        "      def get_dim_size(dim):\n",
        "        dim_val = tensor_shape.dimension_value(a.shape[dim])\n",
        "        if dim_val is not None:\n",
        "          return dim_val\n",
        "        return array_ops.shape(a)[dim]\n",
        "\n",
        "      num_rows = get_dim_size(-2)\n",
        "      num_cols = get_dim_size(-1)\n",
        "      if isinstance(num_rows, int) and isinstance(num_cols, int):\n",
        "        max_rows_cols = float(max(num_rows, num_cols))\n",
        "      else:\n",
        "        max_rows_cols = math_ops.cast(\n",
        "            math_ops.maximum(num_rows, num_cols), dtype)\n",
        "      rcond = 10. * max_rows_cols * np.finfo(dtype).eps\n",
        "\n",
        "    rcond = ops.convert_to_tensor(rcond, dtype=dtype, name='rcond')\n",
        "\n",
        "    [ singular_values, left_singular_vectors, right_singular_vectors, ] = svd(\n",
        "        a, full_matrices=False, compute_uv=True)\n",
        "\n",
        "    cutoff = rcond * math_ops.reduce_max(singular_values, axis=-1)\n",
        "    singular_values = array_ops.where_v2(\n",
        "        singular_values > array_ops.expand_dims_v2(cutoff, -1), singular_values**beta,\n",
        "        np.array(np.inf, dtype))\n",
        "\n",
        "    a_pinv = math_ops.matmul(\n",
        "        right_singular_vectors / array_ops.expand_dims_v2(singular_values, -2),\n",
        "        left_singular_vectors,\n",
        "        adjoint_b=True)\n",
        "\n",
        "    if a.shape is not None and a.shape.rank is not None:\n",
        "      a_pinv.set_shape(a.shape[:-2].concatenate([a.shape[-1], a.shape[-2]]))\n",
        "\n",
        "    return a_pinv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bj90LjpyyyG2"
      },
      "source": [
        "Now we have all pieces in place to run the training. The next cell defines a Python class to organize the neural network optimization. It receives the physics solver, network model, loss function and a data set, and runs as many epochs as possible within the given time limit `MAX_TIME`.\n",
        "\n",
        "Depending on the chosen optimization method, the mini batch updates differ:\n",
        "1. Adam: Compute loss gradient, then apply the Adam update.\n",
        "2. PG: Compute loss gradient und physics Jacobian, invert them data-point-wise, and compute network updates via the proxy loss and Adam.\n",
        "3. HIG: Compute loss gradient and network-physics Jacobian, then jointly compute the half-inversion, and update the network parameters with the resulting step.\n",
        "\n",
        "The `mini_batch_update()` method of the optimizer class realizes these three variants."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 102,
      "metadata": {
        "id": "peqXxMNpyyG5"
      },
      "outputs": [],
      "source": [
        "class Optimization():\n",
        "    def __init__(self,model,solver,loss_function,x_train,y_train):\n",
        "        self.model = model\n",
        "        self.solver = solver\n",
        "        self.loss_function = loss_function\n",
        "        self.x_train = x_train\n",
        "        self.y_train = y_train\n",
        "        self.y_dim = y_train.shape[1]\n",
        "        self.weight_shapes = [weight_tensor.shape for weight_tensor in self.model.trainable_weights]\n",
        "\n",
        "    def set_params(self,batch_size,learning_rate,optimizer,max_time,mode,trunc):\n",
        "        self.number_of_batches = N // batch_size\n",
        "        self.max_time = max_time\n",
        "        self.batch_size = batch_size\n",
        "        self.learning_rate = learning_rate\n",
        "        self.optimizer = optimizer(learning_rate)\n",
        "        self.mode = mode\n",
        "        self.trunc = trunc\n",
        "\n",
        "\n",
        "    def computation(self,x_batch, y_batch):\n",
        "        control_batch = self.model(y_batch)\n",
        "        y_prediction_batch = self.solver(x_batch,control_batch)\n",
        "        loss = self.loss_function(y_batch,y_prediction_batch)\n",
        "        return loss\n",
        "\n",
        "\n",
        "    @tf.function\n",
        "    def gd_get_derivatives(self,x_batch, y_batch):\n",
        "\n",
        "        with tf.GradientTape(persistent=True) as tape:\n",
        "            tape.watch(self.model.trainable_variables)\n",
        "            loss = self.computation(x_batch,y_batch)\n",
        "            loss_per_dp = loss / self.batch_size\n",
        "        grad = tape.gradient(loss_per_dp, self.model.trainable_variables)\n",
        "        return grad\n",
        "\n",
        "\n",
        "    @tf.function\n",
        "    def pg_get_physics_derivatives(self,x_batch, y_batch): # physics gradient for SIP\n",
        "\n",
        "        with tf.GradientTape(persistent=True) as tape:\n",
        "            control_batch = self.model(y_batch)\n",
        "            tape.watch(control_batch)\n",
        "            y_prediction_batch = self.solver(x_batch,control_batch)\n",
        "            loss = self.loss_function(y_batch,y_prediction_batch)\n",
        "            loss_per_dp = loss / self.batch_size\n",
        "\n",
        "        jacy = tape.batch_jacobian(y_prediction_batch,control_batch)\n",
        "        grad = tape.gradient(loss_per_dp, y_prediction_batch)\n",
        "        return jacy,grad,control_batch\n",
        "\n",
        "    @tf.function\n",
        "    def pg_get_network_derivatives(self,x_batch, y_batch,new_control_batch): #physical grads\n",
        "\n",
        "        with tf.GradientTape(persistent=True) as tape:\n",
        "            tape.watch(self.model.trainable_variables)\n",
        "            control_batch = self.model(y_batch)\n",
        "            loss = self.loss_function(new_control_batch,control_batch)\n",
        "            #y_prediction_batch = self.solver(x_batch,control_batch)\n",
        "            #loss = self.loss_function(y_batch,y_prediction_batch)\n",
        "            loss_per_dp = loss / self.batch_size\n",
        "\n",
        "        network_grad = tape.gradient(loss_per_dp, self.model.trainable_variables)\n",
        "        return network_grad\n",
        "\n",
        "    @tf.function\n",
        "    def hig_get_derivatives(self,x_batch,y_batch):\n",
        "\n",
        "        with tf.GradientTape(persistent=True) as tape:\n",
        "            tape.watch(self.model.trainable_variables)\n",
        "            control_batch = self.model(y_batch)\n",
        "            y_prediction_batch = self.solver(x_batch,control_batch)\n",
        "            loss = self.loss_function(y_batch,y_prediction_batch)\n",
        "            loss_per_dp = loss / self.batch_size\n",
        "\n",
        "        jacy = tape.jacobian(y_prediction_batch, self.model.trainable_variables, experimental_use_pfor=True)\n",
        "        loss_grad = tape.gradient(loss_per_dp, y_prediction_batch)\n",
        "        return jacy, loss_grad\n",
        "\n",
        "\n",
        "    def mini_batch_update(self,x_batch, y_batch):\n",
        "        if self.mode==\"GD\":\n",
        "            grad = self.gd_get_derivatives(x_batch, y_batch)\n",
        "            self.optimizer.apply_gradients(zip(grad, self.model.trainable_weights))\n",
        "\n",
        "        elif self.mode==\"SIP\":\n",
        "            jacy,grad,control_batch = self.pg_get_physics_derivatives(x_batch, y_batch)\n",
        "            grad_e = tf.expand_dims(grad,-1)\n",
        "            pinv = tf.linalg.pinv(jacy,rcond=10**-5)\n",
        "            delta_control_label_batch = (pinv@grad_e)[:,:,0]\n",
        "            new_control_batch = control_batch - delta_control_label_batch\n",
        "            network_grad = self.pg_get_network_derivatives(x_batch, y_batch,new_control_batch)\n",
        "            self.optimizer.apply_gradients(zip(network_grad, self.model.trainable_weights))\n",
        "\n",
        "        elif self.mode =='HIG':\n",
        "            jacy, grad = self.hig_get_derivatives(x_batch, y_batch)\n",
        "            flat_jacy_list = [tf.reshape(jac, (self.batch_size * self.y_dim, -1)) for jac in jacy]\n",
        "            flat_jacy = tf.concat(flat_jacy_list, axis=1)\n",
        "            flat_grad = tf.reshape(grad, (-1,))\n",
        "            inv = HIG_pinv(flat_jacy, rcond=self.trunc)\n",
        "            processed_derivatives = tf.tensordot(inv, flat_grad, axes=(1, 0))\n",
        "            #processed_derivatives = self.linear_solve(flat_jacy, flat_grad)\n",
        "            update_list = []\n",
        "            l1 = 0\n",
        "            for k, shape in enumerate(self.weight_shapes):\n",
        "                l2 = l1 + np.prod(shape)\n",
        "                upd = processed_derivatives[l1:l2]\n",
        "                upd = np.reshape(upd, shape)\n",
        "                update_list.append(upd)\n",
        "                l1 = l2\n",
        "            self.optimizer.apply_gradients(zip(update_list, self.model.trainable_weights))\n",
        "\n",
        "\n",
        "    def epoch_update(self):\n",
        "        for batch_index in range(self.number_of_batches):\n",
        "            position = batch_index * self.batch_size\n",
        "            x_batch = self.x_train[position:position + self.batch_size]\n",
        "            y_batch = self.y_train[position:position + self.batch_size]\n",
        "            self.mini_batch_update(x_batch, y_batch)\n",
        "\n",
        "    def eval(self,epoch,wc_time,ep_dur):\n",
        "        train_loss = self.computation(self.x_train,self.y_train)\n",
        "        train_loss_per_dp = train_loss / N\n",
        "        if epoch<5 or epoch%20==0: print('Epoch: ', epoch,', wall clock time: ',wc_time,', loss: ', float(train_loss_per_dp) )\n",
        "        #print('TrainLoss:', train_loss_per_dp)\n",
        "        #print('Epoch: ', epoch,' WallClockTime: ',wc_time,' EpochDuration: ',ep_dur )\n",
        "        return train_loss_per_dp\n",
        "\n",
        "    def start_training(self):\n",
        "        init_loss = self.eval(0,0,0)\n",
        "        init_time = time.time()\n",
        "        time_list = [init_time]\n",
        "        loss_list = [init_loss]\n",
        "\n",
        "        epoch=0\n",
        "        wc_time = 0\n",
        "\n",
        "        while wc_time<self.max_time:\n",
        "            \n",
        "            duration = time.time()\n",
        "            self.epoch_update()\n",
        "            duration = time.time()-duration\n",
        "\n",
        "            epoch += 1\n",
        "            wc_time += duration\n",
        "\n",
        "            loss = self.eval(epoch,wc_time,duration)\n",
        "            time_list.append(duration)\n",
        "            loss_list.append(loss)\n",
        "\n",
        "        time_list = np.array(time_list)\n",
        "        loss_list = np.array(loss_list)\n",
        "        time_list[0] = 0\n",
        "        return time_list, loss_list"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3hCSV9KcyyG7"
      },
      "source": [
        "All that's left to do is to start the training with the chosen global parameters, and collect the results in `time_list`, and `loss_list`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 103,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "6BazTjB5yyG8",
        "outputId": "f4cf85e6-e7d3-4ed1-8c67-b442e1653bfe"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Epoch:  0 , wall clock time:  0 , loss:  1.4401766061782837\n",
            "Epoch:  1 , wall clock time:  28.06755018234253 , loss:  5.44398972124327e-05\n",
            "Epoch:  2 , wall clock time:  31.38792371749878 , loss:  1.064436037268024e-05\n",
            "Epoch:  3 , wall clock time:  34.690271854400635 , loss:  3.163525434501935e-06\n",
            "Epoch:  4 , wall clock time:  37.9914448261261 , loss:  1.1857609933940694e-06\n"
          ]
        }
      ],
      "source": [
        "opt = Optimization(model, solver, loss_function, X_TRAIN, Y_TRAIN)\n",
        "opt.set_params(BATCH_SIZE, LR, OPTIMIZER, MAX_TIME, MODE, TRUNC)\n",
        "time_list, loss_list = opt.start_training()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dPAJ6CGiyyG9"
      },
      "source": [
        "\n",
        "## Evaluation\n",
        "\n",
        "Now we can evaluate how our training converged over time. The following graph shows the loss evolution over time in seconds."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 104,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 316
        },
        "id": "YzHtLIqgyyG9",
        "outputId": "0388cf17-63e8-4587-a168-362f903e19d5"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5xU9bnH8c+znS3ssvSySJVexd6wRNFw1cQGmqq5xsQaS2JuipqbqLElKuQaE40pgl1jbGii2BtSpAkCIiydhV1Y2GXbc/+YA44rOyyws2dn9vt+vea1c87MOec5HJ3v/H6/M+eYuyMiItKQlLALEBGRlk1BISIiMSkoREQkJgWFiIjEpKAQEZGYFBQiIhKTgkL2i5kdbWaLwq5jT8zMzazfPi4738zGNnFJLc6ejqWZ9Qr+HdOasy4Jn4JCGsXMlpvZifXnu/sb7j4gjJqai7sPcffp8Vi3mU03s+/VmzfWzIobeo+Z5ZnZncEx2WZmK8zscTM7dH9qqX8sGzrmjWVmPczsCTPbaGZlZjbPzL4T9XqGmf3SzBYF+7HKzF4ws5Pq1VBhZlvNrNTM3jazi81Mn13NSN8MJCGZWZq714RdR3Mzs0zgFaAUGA8sBLKAU4LHe+FV9yV/B+YABwA7gGFAl6jXHwe6A98CZgXzjge+CrwU9b7/cvd/m1k+cCxwF3Ao8N24Vi+7KJVlv+zm2+9yM7vGzD4KvkU+YmZZUa+PN7PZUd8Oh0e9dp2ZLQ2+PS4ws69FvfYdM3vLzH5nZiXADbupJdXM/idqHR+aWVHUW040s0+CbU82MwuW62tmr5hZSfDt9yEzK6i3TycGz28ws0fN7G/BNuab2Zgm+udsjG8CPYAz3H2eu9e6+zZ3f9zdb9jdAmb2VzO7OnjePeg+uiSY7mtmm8wsJfpYmtnfgZ7Av8ys3Mx+HLXK84NWzEYz+1mMWg8GHgzqq3H3We7+QrD+E4GvAKe7+3vuXhU8XnT3K3a3Mncvc/dngHOBb5vZ0Mb/s8n+UFBIPJwDjAN6A8OB7wCY2SjgAeD7QHvgj8AzwbdkgKXA0UA+cCPwDzPrGrXeQ4FlQGfgN7vZ7lXAROBUoC1wAbA96vXxRD68hgc1nhzMN+BmoBswCChiN0EU5TTgYaAAeAaYFOO9Te1EYJq7b9uLZV4DxgbPjyXyb3hM1PQb7l4XvYC7fxNYQeTbfK673xr18lHAAOAE4JdmNqiB7b4LTDazCWbWczf78Z67F+9muZjc/X2gmMh/K9IMFBQSD3e7+2p33wT8CxgZzL8I+GPwDbLW3f9KpEviMAB3fyxYrs7dHwE+AQ6JWu9qd78n+HZasZvtfg/4ubsv8og57l4S9fot7l7q7iuAV3fW5e5L3P1ld9/h7huAO4l8gDbkTXd/3t1riXSvjNjLf5/67g5aOaVmVgo8G+O9HYC1OyfMbGSw3JYYA9GvAUcF/frHALcCRwavHRu8vjdudPcKd59DpGupof0/G3gD+AXwadCSPLiB/SgM9qPMzCobUcNqoHAv65Z9pKCQeFgb9Xw7kBs8PwC4ut6HYhGRb/KY2beiuqVKgaFEPlB2WrmH7RYRaZXsVV1m1tnMHg4GU7cA/6i33T2tJ8t2cyZQ0A1WHjzujbG+y929YOeDSMunISXArlaWu88Olvk6kLm7Bdx9KbCNSDAeTSSIVpvZAPYtKBo6vvW3u9ndr3P3IURagbOBp4Muv/r7sSnYj4Ma2o96ugOb9rJu2UcKCmlOK4HfRH8ounu2u081swOAPwGXAu2DD415RLqFdtrTpY5XAn33oa6bgnUPc/e2wDfqbXefuPtNQbdNrrtfvL/rC/wHOMnMcvZyudeAs4AMd18VTH8baEfkA3x3muzS0u6+EbidyJeCQiL7cbCZ9djbdQWtku7Am01Vn8SmoJC9kW5mWVGPvT1r7k/AxWZ2qEXkmNlXzSwPyCHywbQBwMy+S6RFsTf+DPyvmfUP1j/czNo3Yrk8oBwoM7PuwLV7ud3m9DdgDfCUmQ0NBvCzgD0NqL9GJIRfD6anB9NvBl1ou7MO6LOvhZrZb4Ma04Jj/ANgibuXuPtLRLr/ng7+e8gws3SCbsgG1tfWzMYTGR/6h7vP3dfaZO8oKGRvPA9URD1u2JuF3X0G8N9EBn83A0sIBrrdfQFwB/AOkQ+oYcBbe1nfncCjRE6t3ALcD7RpxHI3AqOBMuA54Mm93G6zcfdK4DhgAZFatwCLiAzSnxNj0deIBOLOoHgTyI6a3p2bgZ8HXYHX7EO52cBTRE7lXUak6/G0qNe/RqQb7B/Bez4Fzufzkwx2+peZbSXSYvwZkeOsU2ObkenGRSIiEotaFCIiEpOCQkREYlJQiIhITAoKERGJKSkvCtihQwfv1atX2GWIiCSUDz/8cKO7d6w/PymDolevXsyYMSPsMkREEoqZfba7+ep6EhGRmBQUIiISk4JCRERiUlCIiEhMCgoREYlJQSEiIjEpKEREJCYFRZR/L1jHozP2dBM1EZHWJSl/cLcv3J0p76/gzU82MqRbW4Z0yw+7JBGRFkEtioCZcfvZI2iXk86lU2ZRvqMm7JJERFoEBUWUwpwM7p4wis9KtvGzp+aimzqJiCgovuTQPu350YkH8s/ZqzVeISKCgmK3fnhcP47q14Hrn5nPorVbwy5HRCRUCordSE0x7jx3BLmZ6VwyZSbbqzReISKtl4KiAZ3ysrhrwkiWbijn+n/OD7scEZHQKChiOLJfBy47rh+PfVjMkzOLwy5HRCQUCoo9uPyE/hzSu5CfPz2PJevLwy5HRKTZKSj2IC01hbsnjCIrPZVLp8yksro27JJERJqVgqIRuuRnccc5I/h47Vb+99kFYZcjItKsFBSNdNyATnz/2D489N4Knv1oddjliIg0GwXFXrjmpAGM7lnAdU/M5bOSbWGXIyLSLBQUeyE9NYW7J44iNcW4dMosdtRovEJEkp+CYi/1aJfNbWcNZ+6qMm5+/uOwyxERiTsFxT44aUgXvntkLx58eznT5q8NuxwRkbhSUOyj604ZyLDu+Vz72ByKN28PuxwRkbhRUOyjzLRUJp03Cne4bOosqmvrwi5JRCQuFBT74YD2Odx85jBmrSjl9mmLwi5HRCQuFBT7afzwbpx/aE/++PoyXv14fdjliIg0OQVFE/jF+MEM7JLHVY/OZk1ZRdjliIg0KQVFE8hKT2Xy+aPZUVPHFVNnU6PxChFJIgqKJtK3Yy6/+dpQ3l++ibv+80nY5YiINBkFRRP62qgenDOmB5NeXcKbn2wMuxwRkSahoGhiN5w2hH4dc7nykVms31oZdjkiIvtNQdHEsjPSmHz+aMp31HDlw7OprfOwSxIR2S8tPijMrI+Z3W9mj4ddS2Md2DmPX502lLeXljD51SVhlyMisl/iGhRm9oCZrTezefXmjzOzRWa2xMyui7UOd1/m7hfGs854OHtMD84Y2Y3f/3sx7y4rCbscEZF9Fu8WxYPAuOgZZpYKTAZOAQYDE81ssJkNM7Nn6z06xbm+uDEzfv21YfRqn8MVD8+ipHxH2CWJiOyTuAaFu78ObKo3+xBgSdBSqAIeBk5397nuPr7eI6F/6pybmcY9541i8/Zqrnp0DnUarxCRBBTGGEV3YGXUdHEwb7fMrL2Z3QuMMrOfxnjfRWY2w8xmbNiwoemq3U9DuuXzi/GDeW3xBu57Y1nY5YiI7LW0sAvYE3cvAS5uxPvuA+4DGDNmTIv66v6NQ3vyztKN3DZtEQf3asdBBxSGXZKISKOF0aJYBRRFTfcI5iUtM+OWM4fTrSCLy6fOpnR7VdgliYg0WhhB8QHQ38x6m1kGMAF4JoQ6mlXbrHQmTRzN+q2VXPPYR7i3qEaPiEiD4n167FTgHWCAmRWb2YXuXgNcCkwDFgKPuvv8eNbRUowoKuC6Uwbx74Xr+Mtby8MuR0SkUeI6RuHuExuY/zzwfDy33VJdcGQv3llaws0vLOSgA9oxoqgg7JJERGJq8b/MTjZmxu1nD6djbiaXTp3JlsrqsEsSEYlJQRGCguwM7jlvFKtLK7nuCY1XiEjLpqAIyUEHFHLNSQN4fu5a/vHeirDLERFpkIIiRN8/pg/HHtiR/312AfNXl4VdjojIbikoQpSSYtx5zgjaZadz2ZRZlO+oCbskEZEvUVCErH1uJndNGMXykm38/Km5Gq8QkRZHQdECHNanPVeeeCBPz17NYzOKwy5HROQLFBQtxCXH9eOIvu355TPzWLxua9jliIjsoqBoIVJTjN9PGEluZhqXPDST7VUarxCRlkFB0YJ0ysvi9+eOYsmGcm54plVc1UREEoCCooU5qn8HLhnbj0dnFPPULI1XiEj4FBQt0JUn9ueQXoX87Kl5LN1QHnY5ItLKKShaoLTUFO6aOJLMtBQueWgmldW1YZckIq2YgqKF6prfhjvOGcHHa7fy6+cWhF2OiLRiCooW7PiBnbnomD78490VPPvR6rDLEZFWSkHRwl178gBG9Szgp0/M5bOSbWGXIyKtkIKihUtPTeHuCaMwg0unzGJHjcYrRKR5KSgSQFFhNredPYK5q8q45YWPwy5HRFoZBUWCOHlIF75zRC/+8tZyps1fG3Y5ItKKKCgSyE9PHcjQ7m259rE5FG/eHnY5ItJKKCgSSGZaKpMmjqbO4bKps6iurQu7JBFpBRQUCaZXhxxuOXMYs1aUcvtLi8IuR0RaAQVFAho/vBvnHdqTP762jFcXrQ+7HBFJcgqKBPXL8YMZ2CWPqx+dw9qyyrDLEZEkpqBIUFnpqUw6bzSV1bVcPnUWNRqvEJE4UVAksH6dcvn1GUN5f/km7v7PJ2GXIyJJSkGR4L4+ugdnHdSDe15dwpufbAy7HBFJQgqKJPCr04fQt2MuVz4ym/VbNV4hIk1LQZEEsjPSmHzeaLZWVvOjR2ZTW+dhlyQiSURBkSQGdMnjxtOG8NaSEv7w6pKwyxGRJKKgSCLnHlzEaSO68bt/L+a9ZSVhlyMiSUJBkUTMjJu+Poyehdlc/vAsSsp3hF2SiCQBBUWSyc1MY9J5o9m8rZqrH5tDncYrRGQ/KSiS0NDu+fx8/CCmL9rAn95YFnY5IpLgFBRJ6puHHcApQ7tw27RFfPjZ5rDLEZEEpqBIUmbGLWcOp0t+FpdPnUXp9qqwSxKRBKWgSGL5bdKZdN5o1m+t5NrHP8Jd4xUisvcUFEluZFEBPxk3kJcXrOPBt5eHXY6IJCAFRStw4VG9OWFgJ256fiEfFZeGXY6IJBgFRStgZtx+9gg65mZy6ZRZbKmsDrskEUkgCopWol1OBndPHMWq0gp++uRcjVeISKO1+KAws0Fmdq+ZPW5mPwi7nkQ2plchV590IM99tIaH3lsRdjkikiDiGhRm9oCZrTezefXmjzOzRWa2xMyui7UOd1/o7hcD5wBHxrPe1uDiY/pyzIEd+dWzC1iwekvY5YhIAoh3i+JBYFz0DDNLBSYDpwCDgYlmNtjMhpnZs/UenYJlTgOeA56Pc71JLyXFuPOcERS0SefSKTPZtqMm7JJEpIWLa1C4++vApnqzDwGWuPsyd68CHgZOd/e57j6+3mN9sJ5n3P0U4PyGtmVmF5nZDDObsWHDhnjtUlLokJvJXRNGsbxkGz9/ep7GK0QkpjDGKLoDK6Omi4N5u2VmY83sbjP7IzFaFO5+n7uPcfcxHTt2bLpqk9Thfdtz+Qn9eWrWKh77sDjsckSkBUsLu4A9cffpwPSQy0hKlx3fn/eWbeKX/5zHqKIC+nfOC7skEWmBwmhRrAKKoqZ7BPOkmaWmGHdNGElORhqXTJlJRVVt2CWJSAsURlB8APQ3s95mlgFMAJ4JoQ4BOrXN4nfnjmTxunJueGZ+2OWISAsU79NjpwLvAAPMrNjMLnT3GuBSYBqwEHjU3fUJFaJjDuzID8f25ZEZK3l6lhp3IvJFcR2jcPeJDcx/Hp3q2qJc9ZUDef/TTfzsqbkM75FPn465YZckIi1Ei/9ltjSPtNQU7p44ivS0FC6ZMovKao1XiEiEgkJ26VbQhjvOHsHCNVv4zXMLwy5HRFoIBYV8wQmDOvPfR/fm7+9+xvNz14Rdjoi0AAoK+ZJrTx7IiKICfvL4R6wo2R52OSISskYFhZnlmFlK8PxAMzvNzNLjW5qEJSMthUkTR4HBpVNnUlVTF3ZJIhKixrYoXgeyzKw78BLwTSIX/JMkVVSYzW1nDeej4jJueeHjsMsRkRA1NijM3bcDXwf+4O5nA0PiV5a0BOOGduXbhx/AA299yssL1oVdjoiEpNFBYWaHE7l663PBvNT4lCQtyU9PHcSQbm255rE5rCqtCLscEQlBY4PiSuCnwFPuPt/M+gCvxq8saSmy0lOZfN5oauucy6bMpLpW4xUirU2jgsLdX3P309z9t8Gg9kZ3vzzOtUkL0atDDjd9fRgzV5Ryx0uLwy5HRJpZY896mmJmbc0sB5gHLDCza+NbmrQkp43oxsRDirj3taVMX7Q+7HJEpBk1tutpsLtvAc4AXgB6EznzSVqR6/9rCAM653HVo3NYW1YZdjki0kwaGxTpwe8mzgCecfdqQPfPbGWy0lOZfP4oKqpqufzhWdRovEKkVWhsUPwRWA7kAK+b2QHAlngVJS1Xv055/O8ZQ3n/003c/cqSsMsRkWbQ2MHsu929u7uf6hGfAcfFuTZpoc46qAdnju7BPa98wltLNoZdjojEWWMHs/PN7E4zmxE87iDSupBW6lenD6FPhxyueHg2G7buCLscEYmjxnY9PQBsBc4JHluAv8SrKGn5cjLTmHTeaLZWVvOjR2ZTV6chK5Fk1dig6Ovu17v7suBxI9AnnoVJyzeoa1uu/68hvLlkI3+YrvEKkWTV2KCoMLOjdk6Y2ZGArucgTDykiP8a0Y07X17M+59uCrscEYmDxgbFxcBkM1tuZsuBScD341aVJAwz46avDaWoMJvLp85i07aqsEsSkSbW2LOe5rj7CGA4MNzdRwHHx7UySRh5WelMPm80m7ZVcfWjGq8QSTZ7dYc7d98S/EIb4Ko41CMJamj3fH721UG8umgDf35zWdjliEgT2p9boVqTVSFJ4VuHH8DJQzpz64uLmLlic9jliEgT2Z+gUP+CfIGZceuZI+iSn8VlU2ZRtr067JJEpAnEDAoz22pmW3bz2Ap0a6YaJYHkZ6dzz8RRrNtSybWPz8Fd3ydEEl3MoHD3PHdvu5tHnrunNVeRklhG9WzH1ScN4KUF65i+aEPY5YjIftqfrieRBl14VG96tc/m5hcW6iqzIglOQSFxkZGWwo/HDWTxunKemFkcdjkish8UFBI3pwztwqieBdzx0mK2V9WEXY6I7CMFhcSNmfE/pw5i/dYd3P/Gp2GXIyL7SEEhcXVwr0JOGtyZe19bqsuRiyQoBYXE3U9OGUhlTR13/+eTsEsRkX2goJC469sxl4mHFDHl/RUs3VAedjkispcUFNIsrjjhQLLSUrj1xY/DLkVE9pKCQppFx7xMvn9sX6bNX8cHy3XfCpFEoqCQZvO9o3vTKS+Tm55fqEt7iCQQBYU0m+yMNK76yoHMWlHKC/PWhl2OiDSSgkKa1VkH9eDAzrnc+uLHVNXo0h4iiUBBIc0qLTWF604ZyPKS7Ux577OwyxGRRlBQSLM7bkAnDu/TnrtfWcKWSt2zQqSlU1BIs9t5aY9N26q4d/rSsMsRkT1o8UFhZmPN7A0zu9fMxoZdjzSNYT3yOX1kN+5/81PWlFWEXY6IxBDXoDCzB8xsvZnNqzd/nJktMrMlZnbdHlbjQDmQBeh61UnkmpMG4A63T1scdikiEkO8WxQPAuOiZ5hZKjAZOAUYDEw0s8FmNszMnq336AS84e6nAD8BboxzvdKMigqzufDo3jwxs5iX5ut0WZGWKq63M3X3182sV73ZhwBL3H0ZgJk9DJzu7jcD42OsbjOQ2dCLZnYRcBFAz54996NqaU5XntifNz/ZyLWPf8TQ7vl0K2gTdkkiUk8YYxTdgZVR08XBvN0ys6+b2R+BvwOTGnqfu9/n7mPcfUzHjh2brFiJr8y0VO6ZOIqa2jqueHiWbpsq0gK1+MFsd3/S3b/v7ue6+/Sw65Gm16tDDjd9fRgfLN+sS5GLtEBhBMUqoChqukcwT1qx00d25+yDenDPq0t4e+nGsMsRkShhBMUHQH8z621mGcAE4JkQ6pAW5sbTh9C7Qw5XPjybknLdDU+kpYj36bFTgXeAAWZWbGYXunsNcCkwDVgIPOru8+NZhySG7Iw0Jk0cTWlFNVc/Noe6Ol1hVqQliGtQuPtEd+/q7unu3sPd7w/mP+/uB7p7X3f/TTxrkMQyuFtbfvHVQUxftIEH3vo07HJEhAQYzJbW5xuHHcDJQzrz2xc/5qPi0rDLEWn1FBTS4pgZt545gk55WVw2dRZbdeFAkVApKKRFys9O564JIyneXMHPnpqnO+KJhEhBIS3WmF6F/OjE/jwzZzWPfajLfImERUEhLdoPxvbjiL7tuf6f81myfmvY5Yi0SgoKadFSU4zfnzuS7IxU/vtvH+r3FSIhUFBIi9epbRb3fesg1pRV8N0HP6B8R03YJYm0KgoKSQgHHVDIpImjmb96Cz/4x4dU1ejigSLNRUEhCePEwZ25+evDeOOTjVyjX26LNJu43o9CpKmdM6aIkvIqfvvix7TPzeCX4wdjZmGXJZLUFBSScC4+tg8btu7ggbc+pWNeJj8c2y/skkSSmoJCEo6Z8fOvDqJk2w5ufXERHXIzOWdM0Z4XFJF9oqCQhJSSYtx21gg2bavip0/OpTA7gxMHdw67LJGkpMFsSVgZaSnc+42DGNqtLZdMmcmM5ZvCLkkkKSkoJKHlZKbxwHcOpntBGy548AMWr9Ovt0WamoJCEl773Ez+esEhZKWn8q3732flpu1hlySSVBQUkhSKCrP56wWHsL2qhq/94W3mrNR9LESaioJCksagrm154gdHkJWewrn3vcOL89aEXZJIUlBQSFLp3zmPp354JAO7tOUHD83kvteX6l4WIvtJQSFJp2NeJg9fdBinDu3KTc9/zM+enkd1ra4NJbKv9DsKSUpZ6ancM3EUB7TP5g/Tl7Jy03Ymnz+atlnpYZcmknDUopCklZJi/HjcQH575jDeWVrCWf/3NsWbdUaUyN5SUEjSO/fgnvz1gkNYU1bJGZN1RpTI3lJQSKtwZL8OPKkzokT2iYJCWo36Z0TdPm2RBrlFGkFBIa3KzjOizhrdg0mvLuHse9/hs5JtYZcl0qIpKKTVyUpP5bazRzDpvFEs3VDOqXe9wRMfFuv3FiINUFBIqzV+eDdevPIYhnTL5+rH5nDFw7PZUlkddlkiLY6CQlq17gVtmHrRYVz9lQN5bu4aTvn9G7pcuUg9Cgpp9VJTjMtO6M9jFx9Oaopxzh/f4XcvL6ZGA90igIJCZJfRPdvx3OVHccbI7tz1n0849753dclyERQUIl+Ql5XOneeO5K4JI1m8divjfv86f5i+hMrq2rBLEwmNgkJkN04f2Z3nrziaw/t24NYXF3HCHa/xrzmrdWaUtEoKCpEGFBVm8+dvj2HK9w6lbZt0Lps6i7PufYfZugSItDIKCpE9OKJfB5697Ch+e+YwPivZzhmT3+KKh2exqrQi7NJEmoWCQqQRUlOMcw/uyfRrx3Lpcf14cd5ajr99One8tIhtO2rCLk8krhQUInshNzONa04ewCvXjGXc0C7c88oSxt4+nb+/s1yBIUnLknFwbsyYMT5jxoywy5BWYOaKzfz62QXMXFFKbmYaZ47uzjcOO4D+nfPCLk1kr5nZh+4+5kvzFRQi+8fdmbmilH+8+xnPfbSGqto6DutTyDcP68VJQzqTnqqGuyQGBYVIMygp38GjM4p56L3PKN5cQce8TCYe0pOJhxTRNb9N2OWJxJSwQWFmRwPnE7m/92B3P2JPyygoJGy1dc5ri9fz93c+Y/riDaSYceKgTkw4uCdH9GtPZlpq2CWKfElDQZEW540+AIwH1rv70Kj544C7gFTgz+5+S0PrcPc3gDfM7Azgg3jWK9JUUlOM4wd25viBnVm5aTsPvbeCRz5YwbT568jLTOP4QZ04eUgXjj2wIzmZcf3fUGS/xbVFYWbHAOXA33YGhZmlAouBrwDFRD78JxIJjZvrreICd18fLPcocKG7b93TdtWikJZoR00tby3ZyLR563h54To2basiMy2FYw7syMlDunDioE4UZGeEXaa0YqG0KNz9dTPrVW/2IcASd18WFPYwcLq730yk9fElZtYTKIsVEmZ2EXARQM+ePfe/eJEmlpmWuquV8ZvaOj5Yvplp89cybf5aXl6wjtQU47A+hYwb0oXjBnaiR7vssEsWAZphjCIIimejWhRnAePc/XvB9DeBQ9390hjruBGY5u5vN2abalFIInF3Piou48X5a5k2by3LNkZuzVpU2IbD+7Tn8L7tObxPB7rkZ4VcqSS7UFoUTcXdrw+7BpF4MTNGFBUwoqiAH588gCXry3njk428u6yEafPX8eiMYgB6d8jhsD6FHBaER6c8BYc0jzCCYhVQFDXdI5gn0uqZGf0759G/cx4XHNWb2jpn4ZotvLushHeWlvDsnDVMfX8lAH075jC0ez492rWhR7tsitpl06NdG7oVtCEjTb/dkKYTRlB8APQ3s95EAmICcF4IdYi0eKkpxtDu+Qztns/3ju5DTW0dC9Zs4Z2lJby7rISZKzbz7EdrqK37vAvZDDrnZVFUGAmQHu3a7AqRHu2y6VqQpR8Byl6J91lPU4GxQAdgHXC9u99vZqcCvydyptMD7v6bptyuxiikNamprWPd1h2s3LSd4s0VFG/ezspNkb/FmytYU1ZBVI6QYtA1vw3d27X5QogUFUb+dmmbRZqCpFVK2B/c7QsFhcjnqmvrWFtWycogOIp3BUoFKzdvZ+2WSqI/BtJSjK4FWfQoyN7VKolunXTOyyIlxcLbIYmbhB7MFpF9l56aQlFhNkWFuz/dtqqmjtWlFZ+3RoJAWblpO9MXbWD91h1feH9GagrdCrLoVtCGgux08ttkkN8mnYLsdArapJPfJp387PRgXgYFbdLJzkjFTOGSqC9TtnMAAArMSURBVBQUIq1cRloKvTrk0KtDzm5fr6yuZVV0kATdWmvKKlm8rpzS7dWUVVRRXdtw70RailGQnU7bNpEwKciOhMvOgIn+m98mY1foFGRnkKrWS+gUFCISU1Z6Kn075tK3Y26D73F3Kqprg9Co3hUenz+vprSimrLg+fqtlSxet5Wyimq2VjZ8H48Ug8KcDNrnZNI+N4P2uZl0yM2gQ24m7XMi0+1zM+gQvK6WS3woKERkv5kZ2RlpZGek0a1g766SW1Nbx9bKGkorqindHgmXsopqNm2rYtO2KjaWV1FSvoOSbVXMLS6lpLyKrQ3cJCorPYX2OZEwiQ6VjnmZu/7ufN42K02h0kgKChEJVVpqCu1yMmiXkwHsvvurvsrqWkq2BQFSXsXGIEh2TW+rYt2WSuavLqOkvIqaui93i2WkpdAxN5MOeZl03BkiuRlfCJP2uZkU5mS0+lBRUIhIwslKT6V7QRu6N6L1UlfnlFZUs7F8Bxu2Rh67ngd/izdvZ/bKUkq27WB3J4Kmp9oXu8ByMiiMer6zC2zn85wk6wJTUIhIUktJiXzIF+ZkcOAeblFbW+ds2la1K0Q2bYu0UHa2VnZ2hX1Wsp2S8h1sq6rd7XrapKfSrSCL7u2y6V6QRfeCyC/md/7tkp9YP3pUUIiIBFJTbFfXU2N8oQtsWxUl5VVs2raDdVt2sLq0glWlFSxYXcbG8qovLJdi0LltFl3zs+iaHwmOrvlZwd82dM3PolNeZov54aOCQkRkHzW2C2znKcarSytYtTnyt7i0grVllSxcs4X/fLyOyuq6LyyTYtAxL5PObbPokPv5AH37nMg4ys5usA65mbTLTo9rqCgoRETibE+nGLs7WypqWF0WCY81ZZWsLauI/N1Sydqy2APzZtAuO4MOuRnc980xDf4mZl8pKEREQmZmkV+zZ6czqGvbBt9XV+dsqazedcrwxvIqSrbtiJreQW5W03+sKyhERBJESopFLouSnUG/Tg3/ALLJt9tsWxIRkYSkoBARkZgUFCIiEpOCQkREYlJQiIhITAoKERGJSUEhIiIxKShERCQm891dUzfBmdkG4LN9XLwDsLEJy2nJWtO+QuvaX+1r8orn/h7g7h3rz0zKoNgfZjbD3ceEXUdzaE37Cq1rf7WvySuM/VXXk4iIxKSgEBGRmBQUX3Zf2AU0o9a0r9C69lf7mryafX81RiEiIjGpRSEiIjEpKEREJCYFRRQzG2dmi8xsiZldF3Y9TcnMiszsVTNbYGbzzeyKYH6hmb1sZp8Ef9uFXWtTMbNUM5tlZs8G073N7L3g+D5iZhlh19gUzKzAzB43s4/NbKGZHZ7kx/VHwX/D88xsqpllJcuxNbMHzGy9mc2LmrfbY2kRdwf7/JGZjY5XXQqKgJmlApOBU4DBwEQzGxxuVU2qBrja3QcDhwGXBPt3HfAfd+8P/CeYThZXAAujpn8L/M7d+wGbgQtDqarp3QW86O4DgRFE9jkpj6uZdQcuB8a4+1AgFZhA8hzbB4Fx9eY1dCxPAfoHj4uA/4tXUQqKzx0CLHH3Ze5eBTwMnB5yTU3G3de4+8zg+VYiHybdiezjX4O3/RU4I5wKm5aZ9QC+Cvw5mDbgeODx4C1Jsa9mlg8cA9wP4O5V7l5Kkh7XQBrQxszSgGxgDUlybN39dWBTvdkNHcvTgb95xLtAgZl1jUddCorPdQdWRk0XB/OSjpn1AkYB7wGd3X1N8NJaoHNIZTW13wM/BuqC6fZAqbvXBNPJcnx7AxuAvwTdbH82sxyS9Li6+yrgdmAFkYAoAz4kOY/tTg0dy2b7zFJQtDJmlgs8AVzp7luiX/PIudIJf760mY0H1rv7h2HX0gzSgNHA/7n7KGAb9bqZkuW4AgT986cTCchuQA5f7qpJWmEdSwXF51YBRVHTPYJ5ScPM0omExEPu/mQwe93O5mrwd31Y9TWhI4HTzGw5kS7E44n04xcE3RWQPMe3GCh29/eC6ceJBEcyHleAE4FP3X2Du1cDTxI53sl4bHdq6Fg222eWguJzHwD9g7MnMogMkD0Tck1NJuijvx9Y6O53Rr30DPDt4Pm3gX82d21Nzd1/6u493L0XkeP4irufD7wKnBW8LVn2dS2w0swGBLNOABaQhMc1sAI4zMyyg/+md+5v0h3bKA0dy2eAbwVnPx0GlEV1UTUp/TI7ipmdSqRvOxV4wN1/E3JJTcbMjgLeAObyeb/9/xAZp3gU6Enk0uznuHv9wbSEZWZjgWvcfbyZ9SHSwigEZgHfcPcdYdbXFMxsJJFB+wxgGfBdIl8Ck/K4mtmNwLlEzuSbBXyPSN98wh9bM5sKjCVyKfF1wPXA0+zmWAZBOYlI19t24LvuPiMudSkoREQkFnU9iYhITAoKERGJSUEhIiIxKShERCQmBYWIiMSkoBARkZgUFJKQzOx3ZnZl1PQ0M/tz1PQdZnZVjOUfNLOzgufTzWxMI7c7dudly/eh5vI9vF5gZj+Mmu5mZo/HWmYf67jBzFaZ2a9ivKevmc3eU83SOigoJFG9BRwBYGYpRH6gNCTq9SOAt0Ooa38UALuCwt1Xu/tZMd6/P37n7r9s6EV3X+ruI+O0bUkwCgpJVG8DhwfPhwDzgK1m1s7MMoFBwEwz+6WZfRDc5Oa+4NesjWJmB5vZ22Y2x8zeN7O8eq8XmtnTwU1j3jWz4cH8XDP7i5nNDV47s95yHczsHTP7ar1N3gLs/CZ/m5n12nkDGzP7TrCtl81suZldamZXBVeMfdfMCoP39TWzF83sQzN7w8wGNmI/jw22OTtYX96elpHWJW3PbxFpedx9tZnVmFlPIq2Hd4hcxuFwIpeenuvuVWY2yd1/BWBmfwfGA//a0/qD6309Apzr7h+YWVugot7bbgRmufsZZnY88DdgJPALItfdGRasq13UejsTuUbPz9395Xrruw4YuvObfHA5+GhDiVwePgtYAvzE3UeZ2e+AbxG5/Mx9wMXu/omZHQr8gchFEWO5BrjE3d8Kri5cuYf3SyujoJBE9jaRkDgCuJNIUBxBJCjeCt5znJn9mMgNbgqB+TQiKIABwBp3/wBg5yXZ6zVIjgLODF5/xczaB4FyIpGLERK8tjl4mk7kDmWXuPtre7uzwKvBTae2mllZ1H7MBYYHH/JHAI9F1ZnZiPW+BdxpZg8BT7p78T7UJklMXU+SyHaOUwwj0vX0LpEWxRHA22aWReQb9VnBt/s/Efk2HpYaIjfZOXkfl4++yF1d1HQdkS99KURu4DMy6jFoTyt191uIXFivDfBWY7qrpHVRUEgie5tIV9Imd68Nro5aQCQs3ubzUNgYfNvem4HhRUBXMzsYwMzyou53sNMbwPnB62OBjUHL42Xgkp1viup6cuACYKCZ/WQ329wK7PP4QLDtT83s7GC7ZmYj9rScmfV197nu/lsil9tXUMgXKCgkkc0lcrbTu/Xmlbn7xuDe0X8i0tqYRuRDsFGC+6afC9xjZnOIfPjXb43cABxkZh8RGYjeec+AXwPtggH0OcBxUeutBSYCx0efChu8VkLkG/08M7utsbXWcz5wYbDd+TTuvu9XBtv8CKgGXtjHbUuS0mXGRVoZM7sBKHf32xvx3nJ3z41/VdKSqUUh0vqUAxc15gd3RG6eI62cWhQiIhKTWhQiIhKTgkJERGJSUIiISEwKChERien/ASdLu6N130vZAAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.plot(np.cumsum(time_list),loss_list)\n",
        "plt.yscale('log')\n",
        "plt.xlabel('Wall clock time [s]'); plt.ylabel('Loss')\n",
        "plt.title('Linear chain - '+MODE+' with '+str(OPTIMIZER.__name__))\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xnr9BvlDU8H6"
      },
      "source": [
        "For all three methods, you'll see a big linear step right at the start. As we're -- for fairness -- measuring the whole runtime, this first step includes all TensorFlow initialization steps, which are significantly more involved for HIG and SIP. Adam is much faster in terms of initialization, and likewise faster per training iteration.\n",
        "\n",
        "All three methods by themselves manage to bring down the loss. What's more interesting is to see how they compare. For this, the next cell stores the training evolution, and this notebook needs to be run one time with each of the three methods to produce the final comparison graph."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 105,
      "metadata": {
        "id": "bZMVQaJiyyG-"
      },
      "outputs": [],
      "source": [
        "path = '/home/'\n",
        "namet = 'time'\n",
        "namel = 'loss'\n",
        "np.savetxt(path+MODE+namet+'.txt',time_list)\n",
        "np.savetxt(path+MODE+namel+'.txt',loss_list)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KS1lL5LsV9JV"
      },
      "source": [
        "After runs with each of the methods, we can show them side by side:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 106,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 316
        },
        "id": "OsKQ0MyX6Zmx",
        "outputId": "c76d6002-24e3-4f7d-e3a2-042ce15f9478"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3hUZdrH8e89k94rJCGEXgWkg6iAooJiW3VFV9e+9tXVtZe1rRV27Yu9r93VV7GLIoKoNAXpHUII6b1nnvePMwlDSCWTnGRyf67rXDNzzpkz98lAfnme5xQxxqCUUko1xGF3AUoppTo2DQqllFKN0qBQSinVKA0KpZRSjdKgUEop1SgNCqWUUo3SoFBtRkSOFJENdtfRFBExItL/IN+7RkSmermkDktEzhGRr+yuQ7Uv0fMoVGuJyHbgEmPMN3bXcjBExAADjDGb7a5FqY5IWxTK54iIn901+CL9uXZdGhSqzYjIVBFJ9Xi9XURuEJFVIpIvIu+ISJDH8hNF5FcRyRORH0VkhMeyW0Rki4gUishaEfmDx7ILRGSxiDwqItnA3fXU4hSR2zy2sVxEenqscoyIbHJ/9tMiIu739RORb0UkW0SyROS/IhJVZ5+OcT+/W0TeFZHX3J+xRkTGNvLzOUREvhaRHBHZKyK3uecHishjIpLmnh4TkUDPn6mI3CQiGSKyR0ROFZETRGSje1u3eXzG3SLyvvtnXSgiK0Tk0IP9ubrnLXIvF/eyDBEpEJHVIjLMvSzS/XPIFJEdInKHiDg8trtIROaISK6IbBOR4xv6OSn7aVCo9nYmMAPoA4wALgAQkVHAS8BlQCzwLPBxzS9IYAtwJBAJ3AO8ISKJHtudAGwFugP31/O51wNnAycAEcBFQInH8hOBce6azgSmu+cL8CCQBAwBelJPEHk4GXgbiAI+Bp6qbyURCQe+Ab5wb7s/MN+9+HZgIjASOBQYD9zh8fYEIAjoAfwDeB44FxiD9TO6U0T6eKx/CvAeEAO8CXwkIv7uZa35uR4HTAYGut9/JpDtXvake15fYApwHnBhne1uAOKAR4AXa8JZdUDGGJ10atUEbAeOqWf+VCC1znrnerx+BHjG/XwucF+d928ApjTwmb8Cp7ifXwDsbKLGDTXr17PMAEd4vH4XuKWBdU8FVta371gB8o3HsqFAaQPbOdtzO3WWbQFO8Hg9Hdju8TMtBZzu1+Hu+id4rL8cONWjpp88ljmAPcCRB/Nzdc9b5H5+NLARK9QcHus4gQpgqMe8y4AFHtvY7LEsxL0PCXb/W9ap/klbFKq9pXs8LwHC3M97AX93d/3kiUge1l/vSQAicp5Ht1QeMAzrr9Eau5r43J5Yv4BbVJeIdBeRt0Vkt4gUAG/U+dymthPUQN9+Y/UkATs8Xu9wz6uRbYypdj8vdT/u9Vheyr6fK3j8bIwxLiAVL/xcjTHfYrWYngYyROQ5EYlwv9+/nn3o4fE63WM7NS07z5pVB6JBoTqKXcD9xpgojynEGPOWiPTC6l65Gog1xkQBv2N1C9Vo6vC9XUC/g6jrAfe2hxtjIrC6eLzRRbILq1umPmlYwVkjxT3vYNWOxbjHCZKBNG/8XI0xTxhjxmC1ngYCNwJZQGU9+7C7FfugbKRBobzFX0SCPKaWHiHzPHC5iExwD5KGishMd19+KNYvrEwAEbkQ6y/flngBuE9EBri3P0JEYpvxvnCgCMgXkR5Yvwi9YR6QKCJ/cw9eh4vIBPeyt4A7RCReROKwxiHeaMVnjRGR09zfyd+AcuAnWvlzFZFx7u/LHygGygCXu7XzLnC/e796YY0RtWYflI00KJS3fIbV5VEz3d2SNxtjlgF/werKyAU24x7oNsasBf4FLMHqYhkOLG5hff/G+uX1FVAAvAgEN+N99wCjgXzgU+B/LfzcehljCoFjgZOwumE2AUe5F/8TWAasAlYDK9zzDtb/AbOwfq5/Bk4zxlR64ecagRXwuVhdS9nAbPeyv2KFx1ZgEdYg+kut2AdlIz3hTikfJiJ3A/2NMefaXYvqvLRFoZRSqlEaFEoppRqlXU9KKaUapS0KpZRSjfLJi3zFxcWZ3r17212GUkp1KsuXL88yxsTXne+TQdG7d2+WLVtmdxlKKdWpiMiO+uZr15NSSqlGaVAopZRqlAaFUkqpRvnkGIVSSnlDZWUlqamplJWV2V2KVwUFBZGcnIy/v3/TK6NBoZRSDUpNTSU8PJzevXvjK/dVMsaQnZ1Namoqffr0afoNaNeTUko1qKysjNjYWJ8JCQARITY2tkWtJA0KpZRqhC+FRI2W7pMGhYcFuxbw4aYP7S5DKaU6FA0KN2MM7218j/t+uo/1OevtLkcppWrt3buXP/3pT/Tt25cxY8Zw2GGH8eGHH7JgwQIiIyMZNWoUgwYNYvLkycybN8/rn69B4SYi/PPwfxIdGM0N399AcWWx3SUppRTGGE499VQmT57M1q1bWb58OW+//TapqakAHHnkkaxcuZINGzbwxBNPcPXVVzN//nyv1qBB4SE6KJqHJz/MrsJd3LvkXvTKukopu3377bcEBARw+eWX187r1asXf/3rXw9Yd+TIkfzjH//gqaee8moNenish/XpBSQFHcKVh17JU78+xYTECZw24DS7y1JKdQD3fLKGtWkFXt3m0KQI7jrpkEbXWbNmDaNHj272NkePHs3s2bObXrEFNCg87Hj1MrqVbOLI2GSWxsXy4JJ7GZ6bzoDuoyAiCcITwS/A7jKVUl3YVVddxaJFiwgICKg3ENqiJ0SDwsPEQw8hfVU6zqwN3JWfx5+Tornhtyd4Ky2dkJoffmg3KzQiergfEz2e97DCJCDE3h1RSnldU3/5t5VDDjmEDz74oPb1008/TVZWFmPHjq13/ZUrVzJkyBCv1qBB4SFyxh1EzriD33blcdu8tWTv/oHslDe48ZCTeKrPMUjhHijYDQV7IHc77FgMZXkHbig4el94hHsGiUfABEW0+/4ppTqfo48+mttuu425c+dyxRVXAFBSUlLvuqtWreK+++7jhRde8GoNGhT1OLRnFO9cfhhfre3Lnd9nsVC+YOavo5l9/BUckhS5/8oVxVZwFOyGgrR9jzWhkrYSijMP/JCAcHdrJKn+IInoYQWOD57so5RqPhHho48+4rrrruORRx4hPj6e0NBQHn74YQB++OEHRo0aRUlJCd26deOJJ55g2rRp3q3BF4/sGTt2rPHWjYtKKio4/cPzSC3ZRPH2qzlt2GhuOG4QCZFBzd9IVbk7ONI8wsQzXNKgKB2Ma//3+QU13CKpmUK7gUMPXlOqLaxbt87r3TgdRX37JiLLjTEH9Glpi6IJIQEBvHLi45zxyRlEDfiAj3+LZd6qNC49si+XTelHaGAzfoR+gRDd25oaUl0FRXvdrZG0/VsnBWmw6ycrXFyV+7/P4ecOk3qCJLzmMQGczbtKpFJK1aVB0QzdQ7vzwBEPcOX8K/nDMcspTjuFJ77dzJu/7OLvxw3kzLE9cTpa2UXk9IPIHtbUEJcLSrLqtEw8nu9ZBRu+gKrSOm8UCOvecKukJlT8W9BKUkp1GRoUzXRk8pFcOOxCXv79ZWZPOYyLjpjE/Z+u49b/reaVxdu5beYQpgw84J7k3uVwQFg3a0oaWf86xlgD7AeEifsxezNs+wHK8w98b0js/kdv1dflFRjWtvuolOpwNCha4K+j/sqKvSu4+8e7ee/E93j/8sP44vd0HvpiPee/9AtHDojj9plDGJxg4xFNItYgeHA0dG/kcL7ywgPHSTxDJXUplGQf+L7ASI/DgusMvteESlCUDsIr5UN0MLuF0orS+OMnfyQ5PJnXj3+dAGcAFVUuXluynSe/3UxhWSV/HNOT648bSPeITt6VU1lazyB82v5T0V6gzr8h/xCPcZMGBuJD4nQQXnV4Opht0RZFCyWFJXHf4fdx7XfX8u/l/+aW8bcQ4OfgkiP7csaYZJ78djOvLdnOx7+lcdmUvlw6uS8hAZ30x+wfDDF9rakh1ZVQmL4vSGqDxR0qOxZb81xV+7/P4V/nZEWPIKkZhA/rbo3dKKVspf8LD8LRKUdz7pBzeWPdG4xLGMe0FOuY5aiQAO48cSjnHdaLh79Yz2PfbOLNn3dyw/RB/HFMsk/eAAWnP0T1tKaGuKqtc0lqDwuu0zrZvQLWzYPq8v3fJ06Yfj9MvKJt90GpDu7+++/nzTffxOl04nA4ePbZZ7n55puZM2cOY8eOpXfv3oSHhyMiJCQk8Nprr5GQkOC1z9egOEjXjbmOFRkruHPxnQyJGUJSWFLtsl6xofznnDEs35HDffPWcdP7q+gVE8KEvrE2Vmwjh9M6RDc8ARo6qMsYKMnZ/9Dg7x6EnT9pUKgubcmSJcybN48VK1YQGBhIVlYWFRUVB6z33XffERcXx2233cYDDzzAE0884bUatJP4IAU4A5gzeQ7GGG5ceCOVdc9vAMb0iuHJs0cBsCOn/lPulZsIhMZCwnAYOB3GXmR1eZXm2F2ZUrbas2cPcXFxBAYGAhAXF0dSUlKD60+ePJnNmzd7tQZtUbRCz4ie3DXpLm78/kaeXPEk14+9/oB1YsOsq81mFx34F4BqQkgM5O2yuwqlLJ/fAumrvbvNhOFw/EONrnLcccdx7733MnDgQI455hhmzZrFlClTGlx/3rx5DB8+3KtlaouilWb0nsGZA8/k5TUvszB14QHLQwL8CPJ3kFNcXs+7VaNCYuo/RFepLiQsLIzly5fz3HPPER8fz6xZs3jllVcOWO+oo45i5MiRFBQUcOutt3q1Bm1ReMFN42/i18xfuX3R7bx30nskhO4/iBQbGqgtioMRHGN1PRmj52Uo+zXxl39bcjqdTJ06lalTpzJ8+HBeffXVA9apGaNoC9qi8IJAZyBzpsyhvLqcmxfeTFWdQ0G7RQSyt7DMpuo6sZBYqCqDSh3fUV3Xhg0b2LRpU+3rX3/9lV69erVrDRoUXtInsg93TryTFRkrmPvb3P2WJUeHkJpb9/pLqkkhMdZjiQ5oq66rqKiI888/n6FDhzJixAjWrl3L3Xff3a41aNeTF53U7ySWpi/l+VXPM7b7WA5LOgyA5Ohgvvh9D9Uu0/qLB3Ylwe6gKM1p/DwNpXzYmDFj+PHHHw+Yv2DBgtrn27dvb9MatEXhZbeMv4W+kX255YdbyCrNAqygqKw2ZGj3U8uEuM870QFtpWylQeFlIf4hzJkyh5LKEm5ZeAvVrmqSo617aGv3Uwtp15NSHUKHDwoR6SsiL4rI+3bX0lz9o/tz24Tb+Dn9Z55f/TzJ0cEApObqoGyLBGtQKNURtGlQiMhLIpIhIr/XmT9DRDaIyGYRuaWxbRhjthpjLm7LOtvCqf1PZWbfmcz9bS7p5WsA2JWjLYoWCY62HvXsbKVs1dYtileAGZ4zRMQJPA0cDwwFzhaRoSIyXETm1Zm6tXF9bUZEuHPinaSEp/CPH28jLrJSWxQt5fSDoEhtUShlszYNCmPMQqDu//LxwGZ3S6ECeBs4xRiz2hhzYp0poy3ra2uh/qHMnjKbvPI8nN3eZldusd0ldT4hsTqYrZTN7Bij6AF4XsAnlYavKYqIxIrIM8AoEWnwvHQRuVRElonIsszMTO9V20qDYwZz07ibKPFbw5byT+0up/MJjYfiTv33glKtEha2/+2HX3nlFa6++moA7r77bubMmVO77N///jeDBw9m+PDhHHrooVx//fVUVh54wdKW6vCD2caYbGPM5caYfsaYBxtZ7zljzFhjzNj4+Da+d3ULnTnoTFICJ1ISNo/l6SvtLqdzCU+0LjuulGrUM888w1dffcVPP/3E6tWrWbp0Kd26daO0tPVjo3YExW7A8+ypZPc8nyUinJ5yPaYyipsW3kR+eb7dJXUeET2smx354C17lfKm+++/n7lz5xIVFQVAQEAAt9xyCxEREa3eth1nZi8FBohIH6yAOAv4kw11tKv+cXGU7j6b7MBnuWPxHTxx1BO+ecc7b4tIhMpiKC+wBraVssnDvzzM+pz1Xt3m4JjB3Dz+5kbXKS0tZeTIkbWvc3JyOPnkk/dbp6CggKKiIvr06ePV+mq09eGxbwFLgEEikioiFxtjqoCrgS+BdcC7xpg1bVlHR5AUFYSrrCfTEy9mwa4F/Hfdf+0uqXMIT7QetftJdVHBwcH8+uuvtdO9997b5Hu+/PJLRo4cSe/eveu9/EdLtWmLwhhzdgPzPwM+a8vP7mjiw4IA6B94PFN7buBfy//FyG4jGRY3zObKOrgI93EOBWnQbYi9taguram//O0UERFBWFgY27Zto0+fPkyfPp3p06dz4okn1nvb1Jbq8IPZviIi2I8APwdZxRX88/B/Ehccxw3f30BhRaHdpXVsEe4WReEee+tQqoO79dZbueKKK8jLywPAGENZmXeuL6dB0U5EhPiwQDILyokMjGT25NmkF6dz1493YXSgtmHa9aRUs1xxxRVMmzaNCRMmMGLECA4//HBGjRrFqFGjWr1t8cVfUmPHjjXLli2zu4wDnPr0YsKD/Hj94gkAvLj6RR5b8Rh3TLiDWYNn2VxdB/ZIXxhyMpz0mN2VqC5m3bp1DBnim12e9e2biCw3xoytu662KNpRfHggmYX77p194bALObzH4Tyy9BGvH03hUyKStOtJKRtpULSjbnWCwiEOHjjiAaICo7jx+xsprtRLfNQrPAkKfPpUG6U6NA2KdhQfHkhOSQWV1a7aeTFBMTw0+SF2Fu7kvp/u0/GK+kQkWifdKWUDX/w/2dJ90qBoR7FhgRgDuSX7H642LmEcVxx6BZ9u/ZSPNn9kU3UdWEQPKMmCqvKm11XKi4KCgsjOzvapsDDGkJ2dTVBQULPfo/fMbkdRwf4AFJRW0i18/y/pL8P/wrL0ZTzw8wMMjxtO/+j+dpTYMYV7HCIb3dvWUlTXkpycTGpqKh3pQqPeEBQURHJycrPX16BoR5HuoMgrOfBqjk6Hk4cmP8TpH5/ODd/fwJsz3yTEP6S9S+yYas6lKNCgUO3L39+/zS6L0Zlo11M7igppOCgA4oLjeOjIh9iav5WHfnmoPUvr2MKTrMdCPZdCKTtoULSjmhZFfmnD14c/LOkwLhl+CR9u/pBPtnzSXqV1bGHuGx0W6w2MlLKDBkU7igoOACCvkaAAuHLklYzuNpr7frqPbfnb2qO0ji04GsShNzBSyiYaFO0oPMgPkcZbFAB+Dj8envwwgc5Abvj+BsqqvHO9lk7L4bRuiVrsWwOKSnUWGhTtyOEQIoL8yS9p+mqOCaEJ3H/E/WzM3cicZXOaXN/nhcZDcZbdVSjVJWlQtLOoEP8mu55qTE6ezAWHXMA7G97hi+1ftHFlHVxonLYolLKJBkU7iwz2b7LrydM1o69hRPwI7vnxHnYV7GrDyjq40HgNCqVsokHRziKD/Rs8PLY+/g5/Hpn8CCLCDQtvoKK69Tch6ZS060kp22hQtLOWtigAeoT14L7D72Nt9loeXf5oG1XWwYXGWffNruziA/tK2UCDop1FhbQ8KACmpUzjnCHn8Ma6N5i/c34bVNbBhcZbjyXaqlCqvWlQtLOo4ADySipwuVp+kbHrx1zPkJgh3Ln4TtKKuthZyjVBoeMUSrU7DYp2FhXij8tAUUVVi98b4AxgzpQ5uIyLGxfeSKWr5S2TTis8wXpMW2lvHUp1QRoU7az2Mh4tGND2lBKRwt2T7mZV5iqeXPmkN0vr2BJHQfI4+PZ+KM21uxqluhQNinYWFeK+jMdBBgXAjN4z+OPAP/Ly7y/zQ+oP3iqtY3M4YOa/oTQHvv2n3dUo1aVoULSz2ivIlrbuMNebxt3EgOgB3L7odvYW7/VGaR1f4ggYdwksewlyd9hdjVJdhgZFO4tqxhVkmyPIL4g5U+ZQVl3GTQtvosrV8jGPTmnSX8EYWPmG3ZUo1WVoULSzyCbuSdESfSP7cufEO1mRsYJnfnum1dvrFKJSoN/RsPJ1KC+yuxqlugQNinbWnHtStMRJ/U7ilH6n8Nyq51iStsQr2+zwJt8Ahenw+U12V6JUl6BB0c4C/ZwE+zvJa8YVZJvrtgm30SeyD7f+cCtZpV3ghLRek6yw+PW/sOpdu6tRyudpUNggKsSfXC90PdUI8Q9hzpQ5FFUWccsPt1DtqvbatjusKbdAz4kw7zrI2mx3NUr5NA0KGyRFBbMrp8Sr2xwQPYBbx9/Kz3t+5oXVL3h12x2S0w/OeBGcAfDueVBRbHdFSvksDQobDOwezsa9hRjT8st4NOa0AadxfJ/j+c9v/2FZ+jKvbrtDikyG05+HzHXw/sXQFVpSStlAg8IGw3tEkltSyfr0Qq9uV0S467C7SA5L5uaFN5NTluPV7XdI/Y+B4x+BjZ/D5zdbh84qpbxKg8IG0w/pTqCfgyfmb/L6tkP9Q5kzZQ655bncvuh2XMbl9c/ocMb/BQ67GpY+D0uesrsapXyOBoUNYsMCuWbaAD7/PZ3nF271+vaHxA7hxnE3smj3Il5d86rXt98hHXsfDD0FvroD1nxkdzVK+ZQmg0JEDheRUPfzc0Xk3yLSq+1L822XTe7LzOGJ3P/ZOp75fovXxyvOGnQWx/Y6lidWPMGvGb96ddsdksMBf3gWksfDh5fB9sV2V6SUz2hOi2IuUCIihwJ/B7YAr7VpVV2An9PBY2eNZOaIRB76fD1X/ncFmYXlXtu+iHD3pLvpHtqdmxbeRH55vte23WH5B8PZb0NkT3jjdFj/qd0VKeUTmhMUVcb6c/cU4CljzNNAeNuW1TX4Ox08dfYobj1+MPPXZXD0vxbw+pLtVB/ETY3qExEQwezJs8kszeTOxXd6vdXSIYXGwoWfQ/wgePtP8PU/oLqLXAdLqTbSnKAoFJFbgXOBT0XEAfi3bVldh4hw2ZR+fP63IxmRHMmd/7eGk59axHcbMrzyi314/HCuG30d3+36jjfXv+mFijuBsHi46EsYexEsfhzeO1/DQqlWaE5QzALKgYuNMelAMjC7TavqgvrFh/HGxRN44uxRFJRVcuHLSznjmSX8uKX1l+T489A/MyV5CnOWzWFN1hovVNsJ+AfBiY/CjIdg/TwrLErz7K5KqU5Jmvqr1T2QXWaMqRaRgcBg4HNjTIe9D+fYsWPNsmWd94SziioX7y3fxZPzN5NeUMaolCguPqIPMw5JwM95cAeq5ZXl8cd5f8RP/Hj3pHcJD+hCvYdL/mMdDRUUAYddBeP+AsFRdlelVIcjIsuNMWMPmN+MoFgOHAlEA4uBpUCFMeactijUGzp7UNQoq6zmnaW7eGnxNnZkl5AUGcQ5E3vxxzHJdIsIavH2Vmas5MIvLuSYXscwe/JsRKQNqu6g0n6FBQ/Cxi/AGQgDp8OIM6H/sVbrQynVqqBYYYwZLSJ/BYKNMY+IyG/GmEPbqtg6nz8EuBaIA+YbY+Y29R5fCYoa1S7Dt+szeGnRNpZszcbpEKYN7saZY3syeWA8AX7Nb2W8sPoFHl/xOHdOvJMzB53ZhlV3UHt+g1/fhN8/gOJMCIyEwTPhkD9A36ngF2B3hUrZpjVBsRK4EngUa5xijYisNsYMb8aHvgScCGQYY4Z5zJ8BPA44gReMMQ81Y1sO4DVjzLlNretrQeFpa2YR7yzdxfvLU8kuriAiyI/jhyVy0qFJTOwb02TXlMu4uPKbK1mavpQ3Z77JoJhB7VR5B1NdBdu+h9XvW4fRludDUCQMPskdGlPAqcdsqK6lNUExBev8icXGmIdFpC/wN2PMNc340MlAEdYv+GHueU5gI3AskIrVlXU2Vmg8WGcTFxljMkTkZOAK4HVjTJOH7vhyUNSoqHKxeHMWH/+Wxldr0imuqCYuLICZw63QGJ0SjcNRf9dSdmk2f/zkj4T6h/LOie8Q4h/SztV3MFXlsHUBrPnQHRoFEBQFQ060QqOPhobqGg46KDw2EAZgjGnR/SdFpDcwzyMoDgPuNsZMd7++1b3duiFR37Y+NcbMbGDZpcClACkpKWN27NjRkjI7tbLKar5bn8Enq9KYvy6D8ioXCRFBTD+kO9OHJTC+94EtjaXpS7nkq0s4oc8JPHDEA11rvKIxVeWw5Vt3aHwGFYUQHA2DT4Shp0Kfydo9pXxWa1oUw7HOxI4BBMgEzjPGNOs4y3qC4gxghjHmEvfrPwMTjDFXN/D+qcBpQCCwyn3CX6O6QouiIUXlVXyzdi+f/76H7zdmUlbpIjrEn2OGdGfGsAQO7x9HkL8TgLm/zuU/v/2Heyfdyx8G/MHmyjugyrJ9obHhM6gossY0Bk6HISdB/2kQEGp3lUp5TUNB4deM9z4LXG+M+c69oanA88Akr1bYAGPMAmBBe3yWLwgL9OPUUT04dVQPSiqqWLgxky9+T+eL39N5b3kqQf4OJvWL46jB3Zg56M8s27uMB35+gBHxI+gX1c/u8jsW/yAYfII1VZZZ3VPrP7FaGqvfBb9gKyyGnGSFR3C03RUr1SaaExShNSEB1i/umosEHqTdQE+P18nuecrLQgL8mDEskRnDEqmocvHjliy+W5/Btxsy+HZ9BgD9Ek7BxKznyq//xvsnv014oP6FXC//IBg0w5pOrIKdP8K6ebDuE+uEPocfpBxmBcaA6RA3ALQ7T/mI5nQ9fQisAF53zzoXGGOMaVZfRT1dT35Yg9nTsAJiKfCn5nZlNUdX7npqDmMMWzKLrdBYn8HyjJ8I7PkiFI7n6LirOWZId6YMiiciSAdwm+RyQdpKq6Wx8UvIWGvNj+5tBcbA46DXEXquhuoUWjNGEQ3cAxzhnvUD1mB0bjM+9C1gKtY5EHuBu4wxL4rICcBjWEc6vWSMub8F+9IkDYqWKSir5I4Fj/Dd3rdxZp1DXuZw/J3ChD6xTBvSjWOGdKdnTBc/Mqq58nbCpq9g41ewbSFUlYJ/iHXk1MDjYMBx1i1cleqAWn3UU2eiQdFyVa4qLvryIjbkbODOUc+xekcAX6/dy9bMYgAGdAtj2pDuTBvSjVE9ow76UiJdSmUpbF9ktTQ2fWmFCED3YdYtXPtOgZ4TIUBDWHUMLQ4KEfkEaDBFjDEne68879KgODjpxemc8ckZJIQk8N+Z/yXQGcj2rGLmr8/g28tg4/oAACAASURBVPV7+XlrDlUuQ1SIP5P6xTKpXxyT+sXSJy5UD69tijGQucEKjI1fwa6fwFUFzgDoOcE67LbPFOgxWs/ZULY5mKCY0tgGjTHfe6k2r9OgOHjf7/qeq7+9mlmDZnHHxDv2W1ZQVsmiTVl8uz6DHzdnkZZfBkBSZBCH9YtjdK8oRvSIYmBCGIF+TjvK7zzKi2DnEuvs8K3fQ/pqwIB/KPQcB70Oh16ToMcY64ZMSrUD7XpSzTZn6RxeXfsq/5ryL47rfVy96xhj2J5dwqLNWfy4OYslW7PJK7EuKOznEHrFhtC/Wxj94q2pf7cw+nULIyywOQfadUElObD9B6urascS2Ps7YMDhb4VFr8Os8Og53rrUiFJtQINCNVtldSXnf3E+2/K38e5J79IzvGeT7zHGsCunlNW781mTls+WzCI2ZxSxI7uEKvcd+5wO4ahB3ThrXE+mDorXcY7GlObCrl9gx2LY8aN1ZJWrCsRhjXH0mmRNKZOsGzUp5QUaFKpFUgtTOfOTM0mJSOH141/H/yD7zSurXezMKWFzRhErduTywYrdZBWVExcWyMmHJjF1UDzjescQHKBdVY2qKIHUpVZ31Y7FsGupdUQVQOwASJkISSMhYQR0P0TPGFcHRYNCtdg3O77hugXXce6Qc7l5/M1e2WZltYv56zL434pUFmzIpKLaRYDTwZhe0YzvE8PY3tGMSonWLqqmVFVYl0zf+aPV4tj1s9UKAUAgtj8kDIfEEdZjwggI62Zryarja815FPUd/ZQPLAOeNcaUea1KL9Gg8J4Hfn6At9a/xRNHPcFRKUd5ddslFVUs3Z7L4s1ZLNqUxbr0AowBh8CQxAjG9opmTO8YxvaKJilKB3QbZQzkp1qD4umrrMc9qyB/5751QuIgbqB11njcQIgfZD2P7AkObdGp1gXF40A88JZ71iygACs8Iowxf/Zyra2mQeE95dXl/PmzP7O7aDfvn/Q+iWGJbfZZhWWVrNyZx7IduSzbnsOvu/IoqagGoEdUMGN6RTMiOZK+8aH0jg2lZ0wI/jrO0bjSXHd4rIbM9ZC1yTpMtzRn3zp+QVYLpCZAaqbY/nqORxfTmqBYaowZV988EVljjDnEy7W2mgaFd+0s2MmZ886kf1R/Xp7xMv6O9jnOv6raxbo9hSzbkcOyHbks355LesG+BqzTISRHB9MnzgqOvvGhHN4/jn7xYe1SX6dWnA1ZGz2mTdZj3g4wrn3rRaZ4tEA8QiQ0Xq9l5YNaExTrgOnGmJ3u1ynAl8aYISKy0hgzqk0qbgUNCu/7fNvn3LTwJi4adhHXjbnOlhqMMeSWVLItq5jtWcVszy5ma83zrGKK3a2PwQnhnDgikZkjkugTp4O6LVJZBjlb64SIO0gqS/atFxTp0foYAHGDrOfRvcGp40udVWuC4gTgGWAL1v0o+mDdGnUB8BdjzGNer7aVNCjaxt0/3s0Hmz5g7jFzOaLHEU2/oR0ZY0jLL+OrNenMW7WH5Tusgd0hiRGM7x3N8OQohveIpE9caIvuMa7cXC4oTNu/9ZG1ETI3QlH6vvUc/hDT1+q2inU/xvSD2H4QnqitkA6uVUc9iUggMNj9ckNHHMD2pEHRNsqqyjj707PJLs3mvZPeo3tod7tLalBaXimfrd7D12v3snp3fu1Yh9Mh9IwOpm98GCkxIfum2BB6RofoYboHoywfsjbv3wLJ3mK1TKrL963nH+IOjb77wqMmSELjNEQ6gNYGxSSgNx73rzDGvObNAr1Jg6LtbM3bylmfnsXQ2KG8cNwL+Dk6fjdDtcuwNbOINWkFbM0sYktmMVsyi9iVU1LbXQXW76leMSEMSYyonQYnhNMjKrjB+4+rRrhcUJDqDo0t1mPN89zt1gmENQIjPFoi/TyCpJ/eEKodtabr6XWgH/ArUPO/yhhjrvF6lV6iQdG2/m/z/3HH4ju4/NDLuWrkVXaXc9Bqxjx25pSwM6eEbZnFrE8vYN2eArZn7+uPD/J30DcubN8lSbqF0i8+jD5xobW3lVUtVF1lDZznbHUHyOZ9YZK/a/8B9eCYfeERmQwRSRDRY99jSIy2RryktYPZQ00nOjNPg6Lt3b7odj7Z8gnPHfccExMn2l2O1xWVV7EhvYD16YVsdbdANmcUsTuvlJr/CSKQHB1cez2r2mtaxYcSExqgV9Q9WFXlVoujtiWy2f18GxTuAVO9//rOwDrh4fk80Qqa4Gir1eLQ8anGtCYo3gOuMcbsaavivE2Dou2VVJZw1qdnUVBewPsnv09ccJzdJbWL0opqtmVZwbGlphsro4itWUWUVe77KzgqxN8dHqH7gqRbGD2jg/UaV63hqoaiDChIg4LddR49nrsqD3yvOKyjtYKjrSkk1gqRkFgIifZ4HuOxLAb8Att/P23SmqD4DhgJ/ALUjkzp/SjUhpwNnPPZOYzqNopnj30Wh3TdX4Aul2F3Xum+8MgsYkuG9TyraN+AboDTQe+4EPrEhRIfHkhsaCBxYQHEhgUSFeyPv58Df6cDf6cQ4HQQ5O8kJMBJaKCfdnM1l8sFJdlWaBTusU46LM1zP+ZaJxuW5lpX7C3JsV5XFDW8vYAwCIpyB0yUNQW5H4Oj6yyL3rcsMLLTtWBaExT13pdC70ehAN7b+B73LrmXv476K5eOuNTucjqk/JJKNte2QIrYkmGdA5JdVE5uST1/+TYgKsSfntEhJEcHkxwdTO+4UAZ0C6d/tzBiQgPacA+6gKpyd3BkW8FR93lpHpTl7QudmudVjRwAKg6rZRIabx3VFRrfwOReFmj/iaJ6UUDVJowx3LzwZr7c8SUvTX+JMd3H2F1Sp1JV7SKnpIKc4grySiqpqjZUVruoqHZRXuWirLKasspqCsuqSMsrJTW3lF25JaTmllJRta+rKzY0gAHdwxjbK4aJfWMZlRJFqF5Yse1VljYcIqW5UJwFxZkeUxaUF9S/Lf+Q/QMlogckj4OUCRDdp10G7A/mDneLjDFHiEgh+18UULCOeopom1JbT4OifRVVFHHmvDMpry7n/ZPeJzpID2dsay6XYU9BGZv2FrI5wxpoX7ungDVpBVS7DCLQNy6UxMhg4sIC6B4ZxLTB3RnbK1oP9bVbZRmUuAOkKPPAIKl5nrt9X6hE9YKB02HAdOh9BPgHtUlp2qJQbWpt9lrO/excJiZO5KlpT3Xp8Qo7FZVXsWx7Dr/tymftnnwyCsvJKipnb345FdUukiKDGNs7htBAJyUV1RSXV1NUXklxeTUGQ4DTQYCfg0A/J4F+DoIDnAT5OQnyt8ZLAv2t58H+TgI95/tZj0H+jtr5gX7W2EpIgLVcjwJrIVe1dSHHHT/C5m+sW+ZWlVqHCV+zok0+srUn3DmB7ux/wt3Oht9hLw0Ke7y57k0e/OVB/j7m71ww7AK7y1Eeisur+HrtXuatSmNTRhHF5dWEBjoJCfAjPNCP0EAnDhHKq1xUVLkor6qmrNJFWVW1u/vL6gYr9+juagk/h9QOyocEOAkL9CMkwK+2htBAP0IDnIS4H0MDPZa516t5b6h7/S53KZbKUti+2OraGn5Gm3xEawaz/wrcBewFav6VGGPMCK9X6SUaFPYwxnDdguv4ftf3vHL8Kxwaf6jdJSkvc7lM7dhJ7RhKVTXllTXPXZS7H8sqqimpqKK45rG8muLyKqslU1FFcbk1r3ad8qr9zpRvir9T3EHiDpGasAnwIyzQWee1HyGBVsgEBzgJdh9Ntu+5H8H+1usuF0AeWhMUm4EJxpjstirO2zQo7JNfns+sebMwxvDuSe8SGRhpd0mqE3G5DKWVVpCUlFdT5BEsJe6gKa5wz3Mvs9bZFzpF5dX7vS5pQfiA1fqpCY0DQ8UKmhCP5SH1rutX+7zmEOea5x25G66hoGjOYRG7sO5op1STIgMjeWTyI5z/+fn8Y/E/eOyoxzrsfwrV8Tgc4u528oNw72zT5TKUVFotlqLyKkorqymtqKa0spqSiv2fl1XuC5eyepbnFpfWvr+kwtpWZXXLxnkdgjuIrK602lCqCRZ3+NR01cWFBdItIpBu4UF0Cw+ke0RQu1+8sjlBsRVYICKfsv8Jd/9us6pUpzYifgRXj7qax1Y8xg+7f2By8mS7S1JdmMMhhAVa3U9tcdfwymoXpZXV7q42d7jsF0ZV9YTR/sFUWumitKKK9ILK/d5bVFZFlWv/IDp1ZBKPndW+twFqTlDsdE8B7kmpJp039Dz+t+l/PLr8USYlTeoUV5lV6mBYZ9I7iAjy/p0fXS5DXmklGYVl7C0oJ6OgjB423D++yf+9xph72qMQ5Vv8nf5cO/pa/v793/l4y8ecNuA0u0tSqtNxOISY0ABiQgMYnGBfHQ0GhYg8Zoz5m4h8wv4n3AEd+1pPqmM4ttexjIgfwVMrn2JG7xmE+IfYXZJS6iA01qJ43f04pz0KUb5HRPj7mL9z/hfn8/ra17ns0MvsLkkpdRAaDApjzHL3Y4e9+J/q+EZ3H83RPY/mpd9f4vSBp3eZy5Er5UuaPLNERAaIyPsislZEttZM7VGc8g1/G/M3yqvLeea3Z+wuRSl1EJpzCuLLwFygCjgKeA14oy2LUr6lT2Qfzhh4Bu9vfJ9t+dvsLkcp1ULNCYpgY8x8rLO4dxhj7gZmtm1ZytdcfujlBDoDeXzF43aXopRqoeYERbmIOIBNInK1iPwBsP8OG6pTiQuO48JhFzJ/53xW7G2bK18qpdpGc4LiWiAEuAYYA5wLnN+WRSnfdN7Q84gPjudfy/+FL17eXilf1WhQuC8vPssYU2SMSTXGXGiMOd0Y81M71ad8SIh/CFeNvIpVmav4esfXdpejlGqmBoNCRPyMMdXAEe1Yj/Jxp/Q/hf5R/Xl8xeNUVjf/ftFKKfs01qL4xf24UkQ+FpE/i8hpNVN7FKd8j5/Dj+vGXMfOwp28u/Fdu8tRSjVDc8YogoBs4GjgROAk96NSB+XIHkcyPmE8z/72LIUVhXaXo5RqQmNB0U1Ergd+B1a7H9e4H39vh9qUjxIRrh97Pbnlubz0+0t2l6OUakJjQeHEOgw2DOsWImF1pnYhIlNF5AcReUZEprbX56q2dUjsIZzQ5wReX/s66cXpdpejlGpEYxcF3GOMubc1GxeRl7C6qTKMMcM85s8AHscKoxeMMQ81shkDFGF1gaW2ph7VsVwz+hq+3vE1T658kvuPuN/ucpRSDWisReGN+1e+AszYb6PWIbdPA8cDQ4GzRWSoiAwXkXl1pm7AD8aY44GbAb03hg/pEdaD84aex8dbPubbnd/aXY5SqgGNtSimtXbjxpiFItK7zuzxwGZjzFYAEXkbOMUY8yCND5LnAoENLRSRS4FLAVJSUlpRtWpPV468kiV7lnDn4jsZGjuUhFAb786ilKpXgy0KY0xOG31mD2CXx+tU97x6uQ/HfRbr/hhPNbSeMeY5Y8xYY8zY+Ph4rxWr2laAM4DZk2dT5ari5oU3U+WqsrskpVQdzTk81lbGmP8ZYy4zxswyxiywux7lfSkRKfzjsH+wImOFXopcqQ7IjqDYDfT0eJ3snqe6sJl9Z3Jq/1N5btVz/LLnl6bfoJRqN3YExVJggIj0EZEA4CzgYxvqUB3MreNvpVdEL2754RZyytqq51Mp1VJtGhQi8hawBBgkIqkicrExpgq4GvgSWAe8a4xZ05Z1qM4hxD+EOVPmkF+ez+2LbsdlXHaXpJSijYPCGHO2MSbRGONvjEk2xrzonv+ZMWagMaafMUYPoFe1BsUM4sZxN7Jo9yJeX/u63eUopegEg9mq65k1aBbTUqbx2IrHWJOljU2l7KZBoTocEeGeSfcQHxzPjQtvpKiiyO6SlOrSNChUhxQZGMnDkx8mrSiNe3+6V++Ip5SNNChUhzWq2yiuHHkln2/7nI82f2R3OUp1WRoUqkO7eNjFTEiYwIO/PMjWvK12l6NUl6RBoTo0p8PJg0c+SLBfMNd8d42eX6GUDTQoVIcXHxLP40c9zt7ivVz5zZUUVxbbXZJSXYoGheoURnYbyewps1mfs57rvruOyupKu0tSqsvQoFCdxtSeU7nrsLtYsmcJty/WM7eVai+N3Y9CqQ7nDwP+QE5ZDo+teIzYoFhuGncTIt64x5ZSqiEaFKrTuWjYRWSVZvHGujeIDY7lkuGX2F2SUj5Ng0J1OiLCjeNuJKcsh8dXPE5sUCx/GPAHu8tSymdpUKhOySEO/nn4P8krz+OeJfcQHRTN1J5T7S5LKZ+kg9mq0/J3+vPo1EcZEjOEG76/gZUZK+0uSSmfpEGhOrUQ/xCePuZpEkMTuWr+VWzO3Wx3SUr5HA0K1enFBMXwzLHPEOQM4rJvLmN3kd5ZVylv0qBQPqFHWA/mHjOX0qpSzvn0HH7P+t3ukpTyGRoUymcMihnE68e/TpBfEBd+cSHf7PjG7pKU8gkaFMqn9IvqxxsnvMHA6IFcv+B6Xvn9Fb2XhVKtpEGhfE5ccBwvTn+RY3sdy7+W/4v7frqPSpdeG0qpg6XnUSifFOQXxOwps0lZmcILq19gd9Fu5kyZQ3hAuN2lKdXpaItC+SyHOLh29LXcM+keftnzC+d9fh5pRWl2l6VUp6NBoXzeaQNOY+6xc9lbvJc/ffonPSJKqRbSoFBdwsTEibx+gh4RpdTB0KBQXUbdI6KeWPGEDnIr1QwaFKpLqTki6pT+p/D86ue54PML2FWwy+6ylOrQNChUlxPkF8R9h9/H7Cmz2Za/jTM+OYOPt3ys51so1QANCtVlzeg9gw9O/oDBMYO5fdHt3PzDzRRWFNpdllIdjgaF6tISwxJ5afpLXD3yar7a/hVnfHyGXq5cqTo0KFSX53Q4uezQy3j1+FdxiIMLvriA//z6H6pcVXaXplSHoEGhlNuh8Yfy3knvMbPPTOb+NpcLv7iQ1MJUu8tSynYaFEp5CAsI44EjH+DhIx9mc95mTv/4dF5Y/QLl1eV2l6aUbTQolKrHCX1P4P2T32d84ngeX/E4J394Ml9s+0KPjFJdkgaFUg3oEdaDJ49+kheOe4HwgHBuXHgj531+HqszV9tdmlLtSoNCqSZMSJzAOye+wz2T7mFX4S7+9NmfuHnhzewp2mN3aUq1Cw0KpZrB6XBy2oDT+PS0T/nL8L8wf+d8TvroJJ5c+SQllSV2l6dUm9KgUKoFQv1DuWb0NXxy6idMS5nGc6ueY+aHM3l7/dsaGMpniS8Ozo0dO9YsW7bM7jJUF/Bb5m/MXjqb3zJ/I9Q/lJP7ncysQbPoF9XP7tKUajERWW6MGXvAfA0KpVrHGMNvmb/xzoZ3+HL7l1S6KhmXMI5Zg2ZxdMrR+Dv87S5RqWbRoFCqHeSU5fDhpg95b+N77C7aTVxwHGcMPIPTB5xOQmiC3eUp1ahOGxQiciRwDtb9vYcaYyY19R4NCmW3alc1i9MW8/b6t1m0exEOcTC151ROG3AaExMnEuAMsLtEpQ5gS1CIyEvAiUCGMWaYx/wZwOOAE3jBGPNQM7Z1KtDdGPNsU+tqUKiOJLUwlfc2vsf/Nv2PvPI8wvzDmJw8mWkp0ziixxGE+IfYXaJSgH1BMRkoAl6rCQoRcQIbgWOBVGApcDZWaDxYZxMXGWMy3O97F7jYGNPkdaA1KFRHVFFdwU97fmL+zvl8t/M7cstzCXQGMilpEtNSpjG151QiAyPtLlN1YbZ1PYlIb2CeR1AcBtxtjJnufn0rgDGmbkh4biMFuNMY85dG1rkUuBQgJSVlzI4dO7y1C0p5XZWripUZK/lmxzfM3zmfvSV7cYqTsQljOSblGCYnTyYpLMnuMlUX05GC4gxghjHmEvfrPwMTjDFXN7KNe4AvjTE/NucztUWhOhNjDGuy19SGxvaC7YB1CZHxCeMZlzCO8Qnj6R7a3d5Clc9rKCj87CimpYwxd9ldg1JtRUQYFjeMYXHDuHb0tWzN38qStCUsTV/K/J3z+XDzhwD0iujF2O5jGZ8wnvGJ44kLjrO5ctVV2BEUu4GeHq+T3fOU6vJEhH5R/egX1Y9zh55Ltauajbkb+SX9F5amL+XL7V/ywaYPAOgT2YchMUPoEdaDpLAkeoT1oEdYDxJDE/F36rkbynvs6HrywxrMnoYVEEuBPxlj1njrM7XrSfmqKlcVG3I21AbH1vytpBenU22qa9cRhPiQeJLDkkkKSyIpLGm/5wmhCXoSoKqXXUc9vQVMBeKAvcBdxpgXReQE4DGsI51eMsbc783P1aBQXUmVq4rMkkxSi1JJK0ojrShtv+fpJem4jKt2fYc46B7Sfb9WiOfzbiHd8HN0il5p5WWd9oS7g6FBodQ+la5K9hbvJa0ojd1Fu9ldtHu/5xklGRj2/R7wEz+6h3avN0SSwpLoFtINh+j1RH1Rpx7MVkodPH+HP8nhySSHJ9e7vLK6kj3Few4IkN1Fu1m0exGZpZkHbC8xNJHE0EQiAiOIDIwkIsB6jAyItOZ5PEYGRhLsF4yItMfuqjagQaFUF+fv9CclIoWUiJR6l5dXl9d2Y3m2SNKL09mSt4X88nzyK/KpclU1+Bl+4kdEYMS+QKknXGpee64TGRCJ0+Fsq11XzaRBoZRqVKAzkD6RfegT2afBdYwxlFaVUlBRQH55foOPNc8zSzLZkreFgvICCisbvtiCQxxEBUYRExRDbFCs9RhsPdZOwTG1y7Xl0jY0KJRSrSYihPiHEOIf0uKr5Fa5qiiqKCK/In+/UMkrzyOnLMeaSq3HNdlryCnLoaiyqN5tBTmD9guQ2KBYYoNjiQuO2/95cCzh/uEaKs2kQaGUspWfw4+ooCiigqKa/Z6yqjJyy3LJKcshuyz7gEDJKcshsyST9dnrySnLococ2C0W4AioDQ3PMIkLjqudHxMUQ3RQdJcPFQ0KpVSnE+QXRGJYIolhiU2u6zIu8svzyS7NJqssi6zSLLJLs62pLJus0izSitJYnbmanLKc/Y4Aq+Hn8CMm0GqlRAdG1z7GBsdar2u6wNzrhPiF+FSwaFAopXyaQxxEB0UTHRRNf/o3um61q5rc8lwrVEqzyCnLqW25eD7fVbiLnLIcSqrqv096sF8wCaEJJIVaJzgmhSXVHimWGJZIt5BuneqkRw0KpZRyczqctV1PgxjU5PqeXWCeYZJRmkF6cTppRWmsy1lHTlnOfu9ziIP44HgSQhPoHtKd7qHdax8TQhJICE0gLjiuw5z42DGqUEqpTqi5XWBlVWXsKd5jTUV7ap/vLd7LxtyNLExdSFl12X7vcYiDuKA44kPirXGUoH1He3ke+RUbHEtUYFSbhooGhVJKtbEgv6BGDzE2xlBQUUB6cTp7S/ZaU7H1mFGSQUZJRqMD84IQFRhFbHAsjx/1eIPnxBwsDQqllLKZiNSeZDgopuEuL5dxUVhRaB3pVXrgEV/ZZdmE+od6vT4NCqWU6iQc4qgNlL6Rfdvvc9vtk5RSSnVKGhRKKaUapUGhlFKqURoUSimlGqVBoZRSqlEaFEoppRqlQaGUUqpRGhRKKaUaJcYceEndzk5EMoEdB/n2OCDLi+V0ZF1pX6Fr7a/uq+9qy/3tZYyJrzvTJ4OiNURkmTFmrN11tIeutK/QtfZX99V32bG/2vWklFKqURoUSimlGqVBcaDn7C6gHXWlfYWutb+6r76r3fdXxyiUUko1SlsUSimlGqVBoZRSqlEaFB5EZIaIbBCRzSJyi931eJOI9BSR70RkrYisEZFr3fNjRORrEdnkfoy2u1ZvERGniKwUkXnu131E5Gf39/uOiATYXaM3iEiUiLwvIutFZJ2IHObj3+t17n/Dv4vIWyIS5CvfrYi8JCIZIvK7x7x6v0uxPOHe51UiMrqt6tKgcBMRJ/A0cDwwFDhbRIbaW5VXVQF/N8YMBSYCV7n37xZgvjFmADDf/dpXXAus83j9MPCoMaY/kAtcbEtV3vc48IUxZjBwKNY+++T3KiI9gGuAscaYYYATOAvf+W5fAWbUmdfQd3k8MMA9XQrMbauiNCj2GQ9sNsZsNcZUAG8Dp9hck9cYY/YYY1a4nxdi/TLpgbWPr7pXexU41Z4KvUtEkoGZwAvu1wIcDbzvXsUn9lVEIoHJwIsAxpgKY0wePvq9uvkBwSLiB4QAe/CR79YYsxDIqTO7oe/yFOA1Y/kJiBKRxLaoS4Ninx7ALo/Xqe55PkdEegOjgJ+B7saYPe5F6UB3m8rytseAmwCX+3UskGeMqXK/9pXvtw+QCbzs7mZ7QURC8dHv1RizG5gD7MQKiHxgOb753dZo6Ltst99ZGhRdjIiEAR8AfzPGFHguM9ax0p3+eGkRORHIMMYst7uWduAHjAbmGmNGAcXU6Wbyle8VwN0/fwpWQCYBoRzYVeOz7PouNSj22Q309Hid7J7nM0TEHysk/muM+Z979t6a5qr7McOu+rzocOBkEdmO1YV4NFY/fpS7uwJ85/tNBVKNMT+7X7+PFRy++L0CHANsM8ZkGmMqgf9hfd+++N3WaOi7bLffWRoU+ywFBriPngjAGiD72OaavMbdR/8isM4Y82+PRR8D57ufnw/8X3vX5m3GmFuNMcnGmN5Y3+O3xphzgO+AM9yr+cq+pgO7RGSQe9Y0YC0++L267QQmikiI+990zf763HfroaHv8mPgPPfRTxOBfI8uKq/SM7M9iMgJWH3bTuAlY8z9NpfkNSJyBPADsJp9/fa3YY1TvAukYF2a/UxjTN3BtE5LRKYCNxhjThSRvlgtjBhgJXCuMabczvq8QURGYg3aBwBbgQux/gj0ye9VRO4BZmEdybcSuASrb77Tf7ci8hYwFetS4nuBu4CPqOe7dAflU1hdbyXAhcaYZW1SlwaFUkqpxmjXk1JKqUZpUCillGqUBoVSSqlGaVAopZRqlAaFUkqpRmlQKKWUapQGheqURORREfmbx+svReQFj9f/EpHrG3n/KyJyhvv5AhEZ28zPnVpz2fKDK7zGhgAAA8FJREFUqLmoieVRInKlx+skEXm/sfccZB13i8huEbm3kXX6icivTdWsugYNCtVZLQYmAYiIA+sEpUM8lk8CfrShrtaIAmqDwhiTZow5o5H1W+NRY8w/GlpojNlijBnZRp+tOhkNCtVZ/Qgc5n5+CPA7UCgi0SISCAwBVojIP0RkqfsmN8+5z2ZtFhEZJyI/ishvIvKLiITXWR4jIh+5bxrzk4iMcM8PE5GXRWS1e9npdd4XJyJLRGRmnY98CKj5S362iPSuuYGNiFzg/qyvRWS7iFwtIte7rxj7k4jEuNfrJyJfiMhyEflBRAY3Yz+nuD/zV/f2wpt6j+pa/JpeRamOxxiTJiJVIpKC1XpYgnUZh8OwLj292hhTISJPGWPuBRCR14ETgU+a2r77el/vALOMMUtFJAIorbPaPcBKY8ypInI08BowErgT67o7w93bivbYbnesa/TcYYz5us72bgGG1fwl774cvKdhWJeHDwI2AzcbY0aJyKPAeViXn3kOuNwYs0lEJgD/wbooYmNuAK4yxix2X124rIn1VRejQaE6sx+xQmIS8G+soJiEFRSL3escJfL/7d29axRRFMbh3ysEFQz4UVlbaAo/QGyChYrgH6Aikk7BJo2VtiksFMXGMrWVYCsiKCJZAlbJmsLKJmCzKiHbaXIszh0cF5lMFkQ3+z6wsLt3Z+6dZs79mL1Hd8gENweBFVoECuAo8Dki3gNUW7IPDEjOApdL+WtJh0pAuUhuRkgp+1beTpAZymYj4u12LxZ4U5JOrUtaq11HFzhRbvLTwLNaO3e3OO8C8FjSU+B5RKwO0TbbwTz1ZKOsWqc4Tk49LZIjimmgI2kP2aO+Unr382Rv/F/5QSbZuTTk8fVN7jZrnzfJTt8uMoHPqdpraquTRsR9cmO9vcBCm+kqGy8OFDbKOuRU0teI2Ci7o+4ng0WHX0GhV3rb21kY/ggclnQGQNJkLd9B5R0wU8rPAb0y8ngFzFY/qk09BXADOCbp7h/qXAeGXh8odX+SdLXUK0kntzpO0pGI6EbEA3K7fQcK+40DhY2yLvm00+LAd2sR0Su5o+fJ0cZL8ibYSsmbfg14ImmJvPkPjkbmgNOSlsmF6CpnwD3gQFlAXwLO1867AVwHLtQfhS1lX8ge/QdJD9u2dcAMcLPUu0K7vO+3S53LwHfgxZB12w7lbcbNxoykOaAfEY9a/LYfEfv+fqvsf+YRhdn46QO32vzhjkyeY2POIwozM2vkEYWZmTVyoDAzs0YOFGZm1siBwszMGv0Ep6TsxqAlFHAAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "from os.path import exists\n",
        "\n",
        "# if previous runs are available, compare all 3\n",
        "if exists(path+'HIG'+namet+'.txt') and exists(path+'HIG'+namel+'.txt') and exists(path+'SIP'+namet+'.txt') and exists(path+'SIP'+namel+'.txt') and exists(path+'GD'+namet+'.txt') and exists(path+'GD'+namel+'.txt'):\n",
        "  lt_hig = np.loadtxt(path+'HIG'+namet+'.txt')\n",
        "  ll_hig = np.loadtxt(path+'HIG'+namel+'.txt')\n",
        "\n",
        "  lt_sip = np.loadtxt(path+'SIP'+namet+'.txt')\n",
        "  ll_sip = np.loadtxt(path+'SIP'+namel+'.txt')\n",
        "\n",
        "  lt_gd = np.loadtxt(path+'GD'+namet+'.txt')\n",
        "  ll_gd = np.loadtxt(path+'GD'+namel+'.txt')\n",
        "\n",
        "  plt.plot(np.cumsum(lt_gd),ll_gd,   label=\"GD\")\n",
        "  plt.plot(np.cumsum(lt_sip),ll_sip, label=\"SIP\")\n",
        "  plt.plot(np.cumsum(lt_hig),ll_hig,   label=\"HIG\")\n",
        "  plt.yscale('log')\n",
        "  plt.xlabel('Wall clock time [s]'); plt.ylabel('Training loss'); plt.legend()\n",
        "  plt.title('Linear chain comparison ')\n",
        "  plt.show()\n",
        "else:\n",
        "  print(\"Run this notebook three times with MODE='HIG|SIP|GD' to produce the final graph\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XifTlo4xWPDb"
      },
      "source": [
        "This graph makes the significant differences in terms of convergence very clear: Adam (the blue `GD` curve), performs a large number of updates, but its rough approximation of the Hessian is not enough to converge to high accuracies. It stagnates at a high loss level.\n",
        "\n",
        "The SIP updates don't outperform Adam in this scenario. This is caused by the relatively simple physics (the _linear_ oscillators), and the higher runtime cost of SIP. If you run this example longer, SIP will actually overtake Adam, but start suffering from numerical issues with the full inversion.\n",
        "\n",
        "The HIGs perform much better than the two others: despite being fairly slow per iteration, the half-inversion produces a very good update, that makes the training converge to very low loss values very quickly. The HIGs reach an accuracy that is around four order of magnitudes better than the other two methods."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UiaLlPVz-clP"
      },
      "source": [
        "## Next steps\n",
        "\n",
        "There's a variety of interesting directions for further tests and modifications with this notebook:\n",
        "\n",
        "* Most importantly, we've actually only looked at training performance so far! This keeps the notebook reasonably short, but it's admittedly bad practice. While we claim that HIGs likewise work on _real_ test data, this is a great next step with this notebook: allocate proper test samples, and re-run the evaluation for all three methods on the test data.\n",
        "* Also, you can vary the physics behavior: use more oscillators for longer or shorter time spans, or even include a non-linear force (as employed in the HIG paper). Be warned: for the latter, the SIP version will require a new inverse solver, though."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "physgrad-hig-code.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3.8.10 64-bit ('tf250': conda)",
      "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.10"
    },
    "orig_nbformat": 4
  },
  "nbformat": 4,
  "nbformat_minor": 0
}