Multiple_shooting_v_1.ipynb

{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "name": "Multiple shooting v.1.ipynb",
   "provenance": [],
   "collapsed_sections": [],
   "authorship_tag": "ABX9TyPcCUzC8Lwvh2GhuNsJgYbU",
   "include_colab_link": true
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "view-in-github",
    "colab_type": "text",
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/pkhamchuai/Optimization/blob/main/Multiple_shooting_v_1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Multiple shooting for IVP\n",
    "\n",
    "*   Original paper: https://www.sciencedirect.com/science/article/pii/S1474667017612059\n",
    "*   This is the version taught by Prof Dr. Guido Kanschat in Numerical analysis for ODEs class at Universität Heidelberg in 2016\n"
   ],
   "metadata": {
    "id": "bl01MHIlaST_",
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import plotly.graph_objects as go"
   ],
   "metadata": {
    "id": "FPvJsfGB9-A2",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": 1,
   "outputs": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "4C7UQCmYaRRc",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "t0 = 0.0 # time\n",
    "tn = 5.0"
   ]
  },
  {
   "cell_type": "code",
   "source": [
    "eq_nr =2 # number of ODEs\n",
    "node_nr = np.arange(10, 20, 1)\n",
    "\n",
    "ns = node_nr[0] # number of shooting nodes\n",
    "tm = np.linspace(t0, tn, ns+1) # time at shooting nodes/boundaries\n",
    "u0 = np.array([1.0, 0.0]) # initial values of u, u(0)\n",
    "v = np.zeros([eq_nr, ns+1]) # initial values of shooting nodes vk\n",
    "v[:, 0] = u0 # v0 = u0"
   ],
   "metadata": {
    "id": "_nBSrZpOaTQs",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": 3,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Define functions"
   ],
   "metadata": {
    "id": "oX8g3rWXY-To",
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def rhs_evaluate(t, u):\n",
    "  f = lambda t, u: np.array([t*u[1], 4*(u[0]**(3/2))]) # pre-defined ODEs RHS << change\n",
    "  return f(t, u)\n",
    "# --------------------\n",
    "\n",
    "def expl_euler(u0,  t0, tn, h):\n",
    "  u = u0\n",
    "  t = t0\n",
    "  epsilon = 1e-7\n",
    "  y = u0\n",
    "\n",
    "  while True:\n",
    "      f_val = rhs_evaluate(t, u)   # evaluate values of rhs at t\n",
    "      if abs((t+h) - tn) < epsilon:    # stop computation when we reach tn\n",
    "          un = u + (h*f_val)           # un = new u\n",
    "          y = np.vstack([y, un])\n",
    "          break\n",
    "      else:\n",
    "          if t+h > tn:    # reduce step size if we are clse to tn\n",
    "              h = tn - t\n",
    "              un = u + (h*f_val)\n",
    "              y = np.vstack([y, un])\n",
    "              break\n",
    "          else:            # normal step size h\n",
    "            un = u + (h*f_val)\n",
    "            t = t + h\n",
    "      y = np.vstack([y, un])\n",
    "      u = un\n",
    "  return np.transpose(y)"
   ],
   "metadata": {
    "id": "datViejvGlNA",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": 4,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "def integrator_ee(u0, t0, tn, h, y, f):\n",
    "  # Explicit Euler method integrator for Fundamental Matrix\n",
    "  # y is values of u(t,u), needed for evaluation of Fundamental Matrix\n",
    "  # f is the Jacobian/Gradient\n",
    "  u = u0\n",
    "  t = t0\n",
    "  epsilon = 10e-7\n",
    "  i = 0\n",
    "\n",
    "  #print(f'integrator_ee y: {y} {y.shape}')\n",
    "  while True:\n",
    "      \n",
    "      f_val = np.dot(f(t, y[0, i]), u) \n",
    "      # print(f'integrator_ee: {f(t, y[0, i])} {u}')\n",
    "      # print(f'fval: {f_val}')\n",
    "      if abs((t+h) - tn) < epsilon: # stop computation when we reach tn\n",
    "          un = u + (h*f_val)\n",
    "          break\n",
    "      else:\n",
    "          if t+h > tn:    # reduce step size if we are clse to tn\n",
    "              h = tn - t\n",
    "              un = u + (h*f_val)\n",
    "              break\n",
    "          else: # normal step size h\n",
    "              un = u + (h*f_val)\n",
    "              t = t + h\n",
    "      u = un;\n",
    "      i = i+1;\n",
    "  return un\n",
    "\n",
    "def compute_fundamental(a, b, d, y, Gradf_u):\n",
    "  # d = dimension of fundamental matrix\n",
    "  # need all y solved by IVP to evaluate GradF\n",
    "\n",
    "  h = 1/30\n",
    "  U = np.eye(d)\n",
    "  for i in range(d):\n",
    "    U[:, i] = integrator_ee(U[:, i], a, b, h, y, Gradf_u)\n",
    "  return U"
   ],
   "metadata": {
    "id": "cMQKXX8EKubs",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": 5,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Perform multiple shooting"
   ],
   "metadata": {
    "id": "W_eNep_QZB0t",
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def MultipleShooting(v, t0, tn, ns):\n",
    "  \n",
    "  h = 1/30 # step size for Explicit Euler\n",
    "  g = v[:, 0] # boundary condition of original problem\n",
    "  tm = np.linspace(t0, tn, ns+1) # time at shooting nodes/boundaries\n",
    "\n",
    "  Ba = np.array([[1, 0], [0, 0]]) # Ba = dr(v1, vm+1)/dv1 << change\n",
    "  Bb = np.array([[0, 0], [1, 0]]) # Bb = dr(v1, vm+1)/dvm+1 << change\n",
    "\n",
    "  u = np.zeros([eq_nr, ns+1])\n",
    "  GradF = np.zeros([2*(ns+1), 2*(ns+1)])\n",
    "\n",
    "  # analytical Jacobian (not used in AutoDiff)\n",
    "  GradF_func = lambda t, u: [[0, t], [6*np.power(u, (1/2)), 0]]\n",
    "\n",
    "  # --------------------\n",
    "  # compute u(i) and GradF(v)\n",
    "  for i in range(0, ns):   # for every shooting nodes <<< check 1st index!!!\n",
    "      # print(f'Node: {i}')\n",
    "      y = expl_euler(v[:, i], tm[i], tm[i+1], h)\n",
    "      u[:, i] = y[:, -1]      # store values of last y for uk(tk)\n",
    "      \n",
    "      if i == 0:\n",
    "          ua = y[:, 0]        # get ua from 1st subinterval, needed for r <<< check index!!!\n",
    "      if i == ns-1:\n",
    "          ub = y[:, -1]      # get ub from last subinterval, needed for r\n",
    "      \n",
    "      # to compute GradF: put G, -I, Ba, Bb in their positions\n",
    "      # G is computed by integration of U_prime over interval of shooting nodes\n",
    "      # done in compute_Variational func\n",
    "      \n",
    "      # switch between analytical and automatic differentiation (AD)\n",
    "      \n",
    "      GradF[2*i:2*i+eq_nr, 2*i:2*i+eq_nr] = compute_fundamental(tm[i], tm[i+1], eq_nr, y, GradF_func)\n",
    "      # GradF[2*i-1:2*i, 2*i-1:2*i] = compute_fundamental_AD(tm[i], tm[i+1], 2, y);\n",
    "      GradF[2*i:2*i+eq_nr, 2*i+(1*eq_nr):2*i+(2*eq_nr)] = -np.eye(eq_nr)\n",
    "      GradF[2*ns:2*ns+eq_nr, 0:eq_nr] = Ba\n",
    "      GradF[2*ns:2*ns+eq_nr, 2*ns:2*ns+eq_nr] = Bb\n",
    "  print(GradF)\n",
    "  # --------------------\n",
    "  r = np.dot(Ba, ua) + np.dot(Bb, ub) - g     # boundary condition function r\n",
    "  F = compute_F(u, v, r, ns)     # compute F\n",
    "\n",
    "  # print(f'u1(0) = {ua(0)}')\n",
    "  # print(f'u2(0) = {ua(1)}')\n",
    "\n",
    "  # solve linear system GradF*v = -F for v (update of shooting nodes)\n",
    "  deltav = np.linalg.solve(GradF, -F)\n",
    "\n",
    "  # --------------------\n",
    "  # Globalization with line search (with Armijo & Backtracking)\n",
    "  # reduce Newton step size with t, t = (0,1]\n",
    "  # t which minimizes F(u, v+t*deltav, r, ns)\n",
    "  # return deltav = t*deltav\n",
    "  deltav = backtracking(deltav, u, v, r, ns)\n",
    "  # --------------------\n",
    "\n",
    "  return deltav, ua"
   ],
   "metadata": {
    "id": "hw8LG7HVALrQ",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": 9,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "def compute_F(u, v, r, ns):\n",
    "  # to computes F, need u, v, r and ns\n",
    "  F = np.zeros((eq_nr, ns+1))\n",
    "  for i in range(0, ns): # compute F(v)  << beware index!!\n",
    "    F[:, i] = u[:, i] - v[:, i+1]\n",
    "\n",
    "  F[:, -1] = r  # last 2 elements of F(v)\n",
    "  F = np.reshape(F, (-1, 1))\n",
    "  return F\n",
    "\n",
    "def backtracking(deltav, u, v, r, ns):\n",
    "  # compute a new deltav with backtrackng\n",
    "  # find t = (0,1] which minimizes F(u, v+t*deltav, r, ns)\n",
    "  # first check if t = 1 satisfies Armijo (mu = 0)\n",
    "  # if not, reduce t and check again\n",
    "  # returns deltav = t*deltav\n",
    "  # needs u, v, r, ns for compute_F func\n",
    "\n",
    "  t = 1\n",
    "  for i in range(18):     # step from t = 1 back to some t that is not so small\n",
    "      v = np.reshape(v, (-1, 1))\n",
    "      v_new = v + (t*deltav)\n",
    "      v_new = np.reshape(v_new, (eq_nr, ns+1))\n",
    "      v = np.reshape(v, (eq_nr, ns+1))\n",
    "      \n",
    "      # compute norm of F(x+td) and F(x)\n",
    "      normF = np.linalg.norm(compute_F(u, v_new, r, ns))\n",
    "      normF_old = np.linalg.norm(compute_F(u, v, r, ns))\n",
    "      \n",
    "      # if F is \"better\", return deltav, if not, reduce t and repeat\n",
    "      # if there is no t satiefies the condition, use t = 1\n",
    "      if normF < normF_old:\n",
    "          deltav = t*deltav\n",
    "          break\n",
    "      else: t = t - 0.05\n",
    "  return deltav"
   ],
   "metadata": {
    "id": "Ymt3aZC_Mu13",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": 10,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Somehow the Jacobian always explodes, still need to fix\n",
    "The Matlab version runs fine."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "source": [
    "print(f'Nodes = {ns}')\n",
    "iteration = 0\n",
    "while True:\n",
    "  deltav, ua = MultipleShooting(v, t0, tn, ns)\n",
    "  v = np.reshape(v, (eq_nr*(ns+1), 1))\n",
    "  v = v + deltav\n",
    "  v = np.reshape(v, (eq_nr, ns+1))\n",
    "\n",
    "  print(f'Iteration {iteration}')\n",
    "  print(f'Norm v: {np.linalg.norm(deltav)}')\n",
    "  # iterate until deltav is small enough\n",
    "  if np.linalg.norm(deltav) < 1e-7 or iteration > 2:\n",
    "    break\n",
    "  iteration = iteration + 1;\n",
    "\n",
    "print(f'nr. nodes:  {ns}  nr. iterations: {iteration}')\n",
    "print(f'ua: {ua}')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(tm, v[0, :].real, '-', label=\"u_1(t)\")\n",
    "plt.plot(tm, v[1, :].real, '-', label=\"u_2(t)\")\n",
    "plt.legend()\n",
    "plt.xlabel('Time t')\n",
    "plt.ylabel('u(t)')\n",
    "plt.title('Multiple Shooting: {:d} nodes'.format(ns))\n",
    "plt.ylim([-4, 1.5])\n",
    "plt.xlim([0, 5])\n",
    "plt.show()\n",
    "\n",
    "  "
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "hMIzRhYb-SIz",
    "outputId": "e6b32290-aae9-4863-fde8-62a1460eaa03",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": 13,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Nodes = 10\n",
      "[[nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan\n",
      "  -1. -0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan\n",
      "  -0. -1.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  nan nan -1. -0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  nan nan -0. -1.]\n",
      " [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  1.  0.]]\n",
      "Iteration 0\n",
      "Norm v: nan\n",
      "[[nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan\n",
      "  -1. -0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan\n",
      "  -0. -1.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  nan nan -1. -0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  nan nan -0. -1.]\n",
      " [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  1.  0.]]\n",
      "Iteration 1\n",
      "Norm v: nan\n",
      "[[nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan\n",
      "  -1. -0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan\n",
      "  -0. -1.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  nan nan -1. -0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  nan nan -0. -1.]\n",
      " [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  1.  0.]]\n",
      "Iteration 2\n",
      "Norm v: nan\n",
      "[[nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -1. -0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan -0. -1.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan\n",
      "  -1. -0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. nan nan\n",
      "  -0. -1.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  nan nan -1. -0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  nan nan -0. -1.]\n",
      " [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "   0.  0.  1.  0.]]\n",
      "Iteration 3\n",
      "Norm v: nan\n",
      "nr. nodes:  10  nr. iterations: 3\n",
      "ua: [nan nan]\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAXQklEQVR4nO3de5RdZZ3m8e9DEikuoVWIIAQJXrgIItOWtPcrERAQRXuU8Q40g0tEpxVvzFLUQW1kcGzQxaKVRoERFU2LUbloe0FAMQgEEGRQoAmCQBRMWrkov/njvJUqQlWqKlTVrqS+n7Vq5eyz99n7t3elznPe991n71QVkiRt0HUBkqTpwUCQJAEGgiSpMRAkSYCBIElqDARJEmAgaAIkqSRPXsP8a5K8aIzruinJnhNW3OB6f5jk0Ile7wjbGvP+ro+SvCjJsq7r0PgZCDNYe/O9P8kWqz1/eXuTX7AW6zwtyf8a+lxV7VJVP3xExY5t2x9McmOSlUmWJfnKFGyzy/19fJJzkvx2uN9Xkg2TnJrkj0luT/KPk12T1m0Ggm4EDhqYSPI0YOPuylk7Sd4MvBHYs6o2BfqB73db1aR7EDgXePUI848BngJsB7wYeG+SvaemNK2LDASdDrxpyPSbgS8NXWD17pYkb0nyk9VXlOQw4PX03nhWJvlWe35VN1CSY5KcneQrSVYk+UWSpw9XWJINkrw/ya+TLE/y1SSPHWE/ngmcV1W/Bqiq26vqlNWW2S7JRW275w9tGSV5Revqubvt785D5u3cnru7LfOKce7vV5N8qW33miT9Q9b9t61FtiLJ19pxeUiLYyRV9buq+hzw8xEWeTPwsar6Q1VdC/wL8JbhFhz4nSY5PskfWktrnyHzt26tkd8nuSHJPwyZt1FrKf0hyS/p/S5Y7bVfT3JnW++RQ+btkWRJa8X8LskJY9l3TQ4DQT8FNmtverOA1wFnrM2K2hvwmcBxVbVpVe0/wqIHAF8DHgv8X+DfkswZZrl3AK8EXghsDfwB+Owa9uNNSY5K0t/2ZXX/DXgr8DjgUcB7AJLsAHwZeBcwD/gO8K0kj2p1fQs4v73uHcCZSXYcx/6+AjgLeDRwDnBS2+6jgEXAae1YfBl41dAXthB63gjrHVGSxwCPB64c8vSVwC5reNnfAb8CtgCOA76QJG3eWcAyer+H1wAfT/KSNu/DwJPaz170gmigjg3oHb8rgW2AlwLvSrJXW+QzwGeqarP2+q+Od181cQwEwWArYSFwLXDrJG/vsqo6u6oeAE4A+oBnDbPc4cDRVbWsqu6j1wXymiSzV1+wqs6g92a9F/Aj4I4k71ttsX+tquur6s/03nh2b8+/Fvh2VV3Qajoe2Ah4TqtrU+CTVXV/Vf07sJgh3Wxj8JOq+k5V/ZXesR5oET0LmA38c1U9UFXfAC5dbb8eXVUPa42Nwabt33uGPHcPMHcNr7m5qv6l1flFeoGyZZJtgecC76uqe6vqCuDzDLYs/ytwbFX9vqpuAf55yDqfCcyrqo+24/cbei2V17X5DwBPTrJFVa2sqp+uxb5qgjzsD0sz0unAj4HtWa27aJLcMvCgqh5M74yUrYdZbjtgUZIHhzz3V2BLhgmtqjqT3qf3OfRaFmcmuaKqzmuL3D5k8T8x+Ka5NXDzajXdQu8T7V+AW6pqaA03t3ljtfp2+1qobQ3cWg+9wuQtTIyV7d/NgHuHPF4xljqr6k+tcbApsDnw+6oa+tqb6Y3TQG8/bllt3oDtgK2T3D3kuVnAhe3xIcBHgeuS3Ah8pKoWr3nXNFlsIYiqupne4PLLgW8Ms8h/8tCB5q3WtLoxbHLbgQetS2E+8NthlrsF2Kd9Sh746auqNbZg2qftrwFLgV3HUM9v6b1xDdSUVuOtbd62rc4BT2AwkB7J5YJvA7YZ0i0DQ47NI1FVf2jrHzo+83TgmrVY3W+BxyYZ2roYegxu46F1P2HI41uAG1f7Hc6tqpe3Ov9fVR1Erzvun4Czk2yyFjVqAhgIGnAI8JKq+s9h5l0BHJhk4/S+b3DIGtbzO+CJo2zrGUkObJ+S3wXcR28MYHUnA8cm2Q4gybwkBwy3wjYoum+SuW0weh96/eU/G6UW6HUf7Zvkpa118e5W08Xt9X+iN3A8J73vF+xPr099rPs7kkvotXiOSDK77dse41lBkj5gwza5YZse8CXgfyZ5TJKdgH+gN14xLq0b6GLgE0n6kuxG7//AwFjTV4EPtO3Mp9d1N+BSYEWS97XB51lJdk3yzFb/G5LMay2wgVbE0NaYppCBIACq6tdVtWSE2Z8G7qf35vdFegOpI/kC8NQ2GPpvIyzzTXr99n+gd6roga3vfnWfoTcIe36SFfRC4+9GWOcfgQ8C/0HvjeU44G1j6X+vql8BbwBOBO6i94a/f+vzvr9N79PmfQ54U1VdN479HWm79wMH0ntzvbvVsJheGAHQzl56/hpW82cGu4eua9MDPgz8ml4Xzo+AT1XVueOpcYiDgAX0WguLgA9X1ffavI+0bdxIb/D99IEXtfGI/eiN19xI7xh+HvibtsjewDVJVtL7fb+ujfGoA/EGOZpKSY4BnlxVb+i6lukoyc+Ak6vqX7uuRTOPLQSpQ0lemGSr1mX0ZmA3el82k6acZxlJ3dqRXh/8JsBvgNdU1W3dlqSZyi4jSRJgl5EkqVmnuoy22GKLWrBgQddlSNI65bLLLrurquaNttw6FQgLFixgyZKRzoyUJA0nyc2jL2WXkSSpMRAkSYCBIElq1qkxBEkayQMPPMCyZcu49957R194PdXX18f8+fOZM2e424uMzkCQtF5YtmwZc+fOZcGCBTz0ArIzQ1WxfPlyli1bxvbbb79W67DLSNJ64d5772XzzTefkWEAkITNN9/8EbWQDARJ642ZGgYDHun+GwiSJMBAkCQ1BoIkdeBrX/sau+yyCxtssMHDrsBw+eWXc8ghvRsT/vCHP+Tiiy9eNe+kk07i1FNPnZSaDARJ6sCuu+7KN77xDV7wghc8bN7HP/5xjjzySODhgXDwwQdz4oknTkpNnnYqab3zkW9dwy9/+8cJXedTt96MD++/yxqXuemmm9hvv/24+uqrATj++ONZuXIlxxxzzMOW3XnnnYddx4oVK1i6dClPf/rTuemmmzj55JOZNWsWZ5xxBieeeCLPf/7zWbBgAZdeeil77DGuW3CPykCQpGlkyZIl7LrrrkDvgp6HH344m266Ke95z3tWLdPf38+FF15oIEjSaEb7JD+d3Xbbbcybt+YrVT/ucY/juuuum/BtO4YgSRNk9uzZPPjgg6um1+ZLYhtttNGor7v33nvZaKONxr3u0XQaCElOTXJHkqu7rEOSJsKWW27JHXfcwfLly7nvvvtYvHjxuNex8847c8MNN6yanjt3LitWrHjIMtdff/2qbqWJ1HUL4TRg745rkKQJMWfOHD70oQ+xxx57sHDhQnbaaacRl120aBHz58/nkksuYd9992WvvfYCYKedduKee+5ZFQL7778/ixYtYvfdd+fCCy8E4KKLLmLhwoUTXn+qasJXOq4CkgXA4qoaNe76+/vLO6ZJGs6111474pk765pPf/rTzJ07l0MPPfRh8y6//HJOOOEETj/99GFfO9xxSHJZVfWPtt2uWwijSnJYkiVJltx5551dlyNJk+5tb3sbG2644bDz7rrrLj72sY9Nynan/VlGVXUKcAr0WggdlyNJ4/L2t7+diy666CHPvfOd7+Stb33riK/p6+vjjW9847DzJqOraMC0DwRJWpd99rOf7bqEMZv2XUaSpKnR9WmnXwYuAXZMsizJIV3WI0kzWaddRlV1UJfblyQNsstIkgQYCJLUiaOOOoqddtqJ3XbbjVe96lXcfffdq+Z5PwRJmkEWLlzI1VdfzdKlS9lhhx34xCc+sWqe90OQpIny3ffD7VdN7Dq3ehrs88k1LjKe+yG87GUvW/X4Wc96FmeffTbg/RAkaUY79dRTee1rXwt4PwRJmlijfJKfTo499lhmz57N61//eqDb+yEYCJI0QcZ7P4TTTjuNxYsX8/3vf58kwAy+H4IkrU/Gcz+Ec889l+OOO45zzjmHjTfeeNXzM/l+CJK03hjP/RCOOOIIVqxYwcKFC9l99905/PDDgW7vh2CXkSRNoCOPPHLVKaNrMrQVsLqDDz6Yr3zlKxx66KHssMMOLF26dNW8yy+/nF122YXNN998QuodyhaCJE0z3g9BktZD3g9BkjpQVavO1pkupvJ+CI/0lsh2GUlaL/T19bF8+fJH/Ka4rqoqli9fTl9f31qvwxaCpPXC/PnzWbZsGTP53ut9fX3Mnz9/rV9vIEhaL8yZM4ftt9++6zLWaXYZSZIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAjoOhCR7J/lVkhuSvL/LWiRppussEJLMAj4L7AM8FTgoyVO7qkeSZrouWwh7ADdU1W+q6n7gLOCADuuRpBmty0DYBrhlyPSy9txDJDksyZIkS2byvVIlabJN+0Hlqjqlqvqrqn/evHldlyNJ660uA+FWYNsh0/Pbc5KkDnQZCD8HnpJk+ySPAl4HnNNhPZI0o83uasNV9ZckRwDnAbOAU6vqmq7qkaSZrrNAAKiq7wDf6bIGSVLPtB9UliRNDQNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSUBHgZDk75Nck+TBJP1d1CBJeqiuWghXAwcCP+5o+5Kk1czuYqNVdS1Aki42L0kaxrQfQ0hyWJIlSZbceeedXZcjSeutSWshJPkesNUws46uqm+OdT1VdQpwCkB/f39NUHmSpNVMWiBU1Z6TtW5J0sSb9l1GkqSpMa5ASPKYJLskeWKStQ6TJK9Ksgx4NvDtJOet7bokSRNj1C6jJH8DvB04CHgUcCfQB2yZ5KfA56rqB+PZaFUtAhaNv1xJ0mQZyxjC2cCXgOdX1d1DZ7Qvlb0hyROr6guTUaAkaWqMGghVtXAN85YASya0IklSJ8Y8DpDk+2N5TpK0bhrLGEIfsDGwRZLHAANfL94M2GYSa5MkTaGxjCH8d+BdwNbAZQwGwh+BkyapLknSFBvLGMJngM8keUdVnTgFNUmSOjDqGEKS5wGMFAZJNkuy60QXJkmaWmPpMnp1kuOAc+l1GQ18D+HJwIuB7YB3T1qFkqQpMZYuo/+R5LHAq4G/p3fBuj8D1wInV9VFk1uiJGkqjOnidlX1+yRbAze2H4ACXgoYCJK0HhjP1U5XDnncB+xHr5UgSVoPjDkQqup/D51OcjzgRekkaT3xSC5/vTEwf6IKkSR1a8wthCRX0Rs3AJgFzAM+OhlFSZKm3njGEPYb8vgvwO+q6i8TXI8kqSPjGUO4eTILkSR1y1toSpIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAjoKhCSfSnJdkqVJFiV5dBd1SJIGddVCuADYtap2A64HPtBRHZKkppNAqKrzh9x+86fA/C7qkCQNmg5jCAcD3x1pZpLDkixJsuTOO++cwrIkaWYZ8z2VxyvJ94Cthpl1dFV9sy1zNPAX4MyR1lNVpwCnAPT399cklCpJYhIDoar2XNP8JG8B9gNeWlW+0UtSxyYtENYkyd7Ae4EXVtWfuqhBkvRQXY0hnATMBS5IckWSkzuqQ5LUdNJCqKond7FdSdLIpsNZRpKkacBAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRLQUSAk+ViSpUmuSHJ+kq27qEOSNKirFsKnqmq3qtodWAx8qKM6JElNJ4FQVX8cMrkJUF3UIUkaNLurDSc5FngTcA/w4jUsdxhwGMATnvCEqSlOkmagVE3Oh/Mk3wO2GmbW0VX1zSHLfQDoq6oPj7bO/v7+WrJkyQRWKUnrvySXVVX/aMtNWguhqvYc46JnAt8BRg0ESdLk6eoso6cMmTwAuK6LOiRJg7oaQ/hkkh2BB4GbgcM7qkOS1HQSCFX16i62K0kamd9UliQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgR0HAhJ3p2kkmzRZR2SpA4DIcm2wMuA/+iqBknSoC5bCJ8G3gtUhzVIkprZXWw0yQHArVV1ZZLRlj0MOKxN3pfk6smubx2xBXBX10VMEx6LQR6LQR6LQTuOZaFUTc4H9CTfA7YaZtbRwAeBl1XVPUluAvqratRfXJIlVdU/sZWumzwWgzwWgzwWgzwWg8Z6LCathVBVew73fJKnAdsDA62D+cAvkuxRVbdPVj2SpDWb8i6jqroKeNzA9HhaCJKkybOufQ/hlK4LmEY8FoM8FoM8FoM8FoPGdCwmbQxBkrRuWddaCJKkSWIgSJKAdSQQkuyd5FdJbkjy/q7r6VKSU5PcMdO/j5Fk2yQ/SPLLJNckeWfXNXUlSV+SS5Nc2Y7FR7quqWtJZiW5PMnirmvpUpKbklyV5IokS0ZdfrqPISSZBVwPLASWAT8HDqqqX3ZaWEeSvABYCXypqnbtup6uJHk88Piq+kWSucBlwCtn4v+L9M7f3qSqViaZA/wEeGdV/bTj0jqT5B+BfmCzqtqv63q6Mt6zONeFFsIewA1V9Zuquh84Czig45o6U1U/Bn7fdR1dq6rbquoX7fEK4Fpgm26r6kb1rGyTc9rP9P6kN4mSzAf2BT7fdS3rmnUhELYBbhkyvYwZ+oev4SVZAPwX4GfdVtKd1kVyBXAHcEFVzdhjAfwfetdJe7DrQqaBAs5Pclm7DNAarQuBII0oyabA14F3VdUfu66nK1X116rand43//dIMiO7E5PsB9xRVZd1Xcs08byq+ltgH+Dtrct5ROtCINwKbDtken57TjNc6y//OnBmVX2j63qmg6q6G/gBsHfXtXTkucArWt/5WcBLkpzRbUndqapb2793AIvodcGPaF0IhJ8DT0myfZJHAa8Dzum4JnWsDaR+Abi2qk7oup4uJZmX5NHt8Ub0TsC4rtuqulFVH6iq+VW1gN57xb9X1Rs6LqsTSTZpJ1yQZBN6959Z49mJ0z4QquovwBHAefQGDr9aVdd0W1V3knwZuATYMcmyJId0XVNHngu8kd4nwCvaz8u7Lqojjwd+kGQpvQ9QF1TVjD7dUgBsCfwkyZXApcC3q+rcNb1g2p92KkmaGtO+hSBJmhoGgiQJMBAkSY2BIEkCDARJUmMgSE2SzYecwnp7klvb45VJPjcJ23tlkqdO9HqlteVpp9IwkhwDrKyq4ydxG6cBi6vq7MnahjQethCkUSR50cB19ZMck+SLSS5McnOSA5Mc1645f267nAZJnpHkR+2iYue1y3UPXedzgFcAn2qtkCdN/Z5JD2UgSOP3JOAl9N7QzwB+UFVPA/4M7NtC4UTgNVX1DOBU4NihK6iqi+ldguWoqtq9qn49lTsgDWd21wVI66DvVtUDSa4CZgEDlwO4ClgA7AjsClzQu+QSs4DbOqhTGhcDQRq/+wCq6sEkD9TgQNyD9P6mAlxTVc/uqkBpbdhlJE28XwHzkjwbepfpTrLLMMutAOZOaWXSGhgI0gRrt3p9DfBP7UqTVwDPGWbRs4Cj2s3gHVRW5zztVJIE2EKQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1Px/BLDg7s6HiAUAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "db0dGVRUU9OJ",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}