tutorials/tutorial7/tutorial-7-dmdc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyDMD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial 7: Dynamic mode decomposition with control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we will show how to extend the dynamic mode decomposition to incorporate the effect of control (this technique has been introduced in the paper [Dynamic mode decomposition with control](https://arxiv.org/abs/1409.6358))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all we import the `DMDc` class from the pydmd package, we set matplotlib for the notebook and we import numpy and scipy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "from pydmd import DMDc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we create our dataset: since we want to add the control, the evolution of the complex system can be formally summarized as:\n",
    "$$\n",
    "\\mathbf{x}_{k+1} = \\mathbf{A}\\mathbf{x}_k + \\mathbf{B}\\mathbf{u}_k,\n",
    "$$where the operators $\\mathbf{A}$ and $\\mathbf{B}$ are the ones we will approximate using DMD. So, for a demostrative purpose, we create the original snapshots by using two random operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_system(n, m):\n",
    "    A = scipy.linalg.helmert(n, True)\n",
    "    B = np.random.rand(n, n)-.5\n",
    "    x0 = np.array([0.25]*n)\n",
    "    u = np.random.rand(n, m-1)-.5\n",
    "    snapshots = [x0]\n",
    "    for i in range(m-1):\n",
    "        snapshots.append(A.dot(snapshots[i])+B.dot(u[:, i]))\n",
    "    snapshots = np.array(snapshots).T\n",
    "    return {'snapshots': snapshots, 'u': u, 'B': B, 'A': A}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We got 10 snapshots of the evolving system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(25, 10)\n"
     ]
    }
   ],
   "source": [
    "s = create_system(25, 10)\n",
    "print(s['snapshots'].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can compute as usually the DMD algorithm on the data: the `fit` method in this version take as arguments the snapshots and the control input (the $\\mathbf{B}$ operator can be also passed). In this case, we do not perform any truncation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pydmd.dmdc.DMDc at 0x11ed66f70>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dmdc = DMDc(svd_rank=-1)\n",
    "dmdc.fit(s['snapshots'], s['u'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us visualize the original system and the reconstructed one: also because without truncation, the plots are the same!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4gAAAF1CAYAAACqI8QgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1rklEQVR4nO3deZhlVXnv8e/L1N3MYDODgEyCKHgBRYgGxQG9GmISicQQh9wQk2jU6yyJJjGDUzQmGpOOIpoYh6A4oghGIxokICCzzENDQ9M0NN009FTv/WPv1pLbQ1F1Tq2zVn0/z1MPdab9rjrd9K/es9baOzITSZIkSZI2KT0ASZIkSdJosEGUJEmSJAE2iJIkSZKkng2iJEmSJAmwQZQkSZIk9WwQJUmSJEmADaIKiIh3RMTHB/3cCRwrI2L/QRxLkiSNloj4XkT8n9LjkGpng6gpiYhXRMQVEbE8Iu6KiI9FxPYbek1m/nVmTugf8Efz3FESEbdExLNLj0OSNNr6vHgoIpb1OXpGRGxdelzrMswPWiNin/74mw3j+IPW/zn9ZelxSMNgg6hJi4g3Au8F3gxsBxwN7A2cGxFbrOc1VfzDL0nSNHpRZm4NHA48GXh72eFMjhkvtcEGUZMSEdsCfw68NjO/lZmrMvMW4CRgH+C3++f9WUScGRH/FhEPAK/o7/u3ccf6nYi4NSLujYg/HT/7Nv654z5dfHlE3BYRiyLitHHHeUpEXBAR90fEgoj4yPoa1XX8PK+IiJsiYmlE3BwRL4uILSJicUQ8cdzzdu5nS3eKiLkR8fW+3uKIOD8iNomIfwUeC3yt/0T4Lf1rj46I/+6f/5OIOG7ccb8XEX/ZP74sIr4WEY+JiM9ExAMRcVFE7DOJPypJUiUy8y7gHLpGEdhoduwYEZ+MiDsj4r6I+PK4x34vIm7o8+mrEbH7uMcyIl4dEdf3x/1oRET/2P4R8V8RsaTP2c/393+/f/lP+pz6zYg4LiLmR8RbI+Iu4JN9nv5g/M81fuYxIuZExN/2ub8kIn4QEXOAtce/vz/+0/rnvyoirul/vnMiYu9xx31ORFzbH+cjQKzvve1/R7i4z9S7I+KD/f3fiIjXPuK5l0fEi6PzoYhY2L/uiog4NCJOBV4GvGVtZvev2z0ivhgR9/S/S/zxuGP+WUT8R3S/Dy3tj3VgRLy9P/7tEfHc9Y1fmk42iJqsY4DZwJfG35mZy4CzgeeMu/tE4Exge+Az458fEYcA/0j3D+1udDORe2yk9i8BBwHHA++MiIP7+9cAbwDmAk/rH//Djf0gEbEV8PfA8zNzm/5nuywzVwKfo292eycD38nMe4A3AvOBnYBdgHd0b0GeAtxG/4lwZr4vIvYAvgH8JbAj8CbgixGx07hjvxQ4pf/59wMuAD7ZP/8a4F0b+1kkSfWKiD2B5wM39Lc3lh3/CmwJPAHYGfhQ/7pnAX9D96HtbsCtdHk23guBo4An9c97Xn//u4FvAzsAewL/AJCZz+gfP6zPts/3t3ftx7Y3cOoEfswPAEfQZe2OwFuAMWDt8bfvj39BRJxIl62/Rpe15wOf7X/GuXS/g/wJXe7fCBy7gbofBj6cmdvSZewX+vs/xbicj4jD6HL4G8Bz+3EdSPf7yUnAvZk5j+73mff1Y31RRGwCfA34Sf/644HXR8Ta9xXgRXR/ZjsAl9J9GLBJ//y/AP55Au+fNHQ2iJqsucCizFy9jscW9I+vdUFmfjkzxzLzoUc89zeAr2XmD/qG7J1AbqT2n2fmQ5n5E7p/iA8DyMwfZ+aPMnN1P5v5z8AvT/DnGQMOjYg5mbkgM6/q7/8UcPLaT1bpGrh/7b9fRRe8e/czqOdn5vrG/tvA2Zl5dv8+nAtcDLxg3HM+mZk3ZuYS4JvAjZl5Xv8e/wfdsiNJUnu+HBFLgduBhfz8A8H1ZkdE7EbXTL46M+/rc+i/+te9DDg9My/JzBV0S1af9oiVKO/JzPsz8zbgu/x81nIVXbO3e2Y+nJm/MBu4DmPAuzJzxToy/hf0TdSrgNdl5h2ZuSYz/7sf47q8GvibzLymz8K/Bg7vZxFfAFyVmWdm5irg74C7NlB+FbB/RMzNzGWZ+aP+/q8CB0bEAf3tU4DP97+TrAK2AR4PRD+OBes5/lHATpn5F5m5MjNvAv6F7sPftc7PzHPG5fpOdH8Oq+ga+H1iI+dxkKaDDaImaxEwN9a932C3/vG1bt/AcXYf/3hmLgfu3Ujt8QGwHNgaoF+q8fXoNvk/QBckc9d1gPEy80HgN+mCaEG/3OTx/WMX9jWO6+/bny5MAN5P9ynvt6Nbnvq2DZTZG3hJv5Tn/oi4n24mdLdxz7l73PcPreP2SJ60QJI0Zb/ar2A5jq4ZWZtdG8qOvYDFmXnfOo63O92sIfCz1T338osrdNaZpXQzegH8T0RcFRGv2sjY78nMhzf+IwLdzzWbbrZvIvYGPjzuZ1/cj20P/v/fH5IN/77xu3QzgddGt23jhf3rHgY+D/x238CeTP9BcGb+J/AR4KPAwoiYF90Wm/WNdfdH/Fm9g26F0VqPzPVFmblm3G0w6zUCbBA1WRcAK+iWffxMdGdeez7wnXF3b2hGcAHdEpa1r58DPGaSY/oYcC1wQL+E5B1sYD/CeP0nes+hC91r6T71W2vt8pNTgDPXBmFmLs3MN2bm44BfAf5vRBy/9pCPKHE78K+Zuf24r60y8z2T+1ElSa3pZwDPoFuGCRvOjtuBHdcz43QnXcMC/GwrxWOAOyYwhrsy8/cyc3fg94F/jA2fufSRefcg3bLXtbV3HffYIuBhuiWeGzsOdD/j7z/i55+Tmf9N9/vDXuPqxPjb6/i5rs/Mk+mW4r4XOLN/X6DL+ZfRLQtdnpkXjHvd32fmEcAhdA3mm9cz3tuBmx8x1m0y8wVIlbFB1KT0yyD/HPiHiDghIjbvl658gW5f3r9u6PXjnAm8KCKOie6EMn/GBJu6ddgGeABY1s/2/cFEXhQRu0TEiX1QrACW0S2ZWevfgBfTNYmfHve6F0a3mT+AJXR7INe+7m7gcY84xosi4nkRsWlEzI5uc/+eSJL0c38HPKffC7fe7OiXOn6TroHboc/htfv4Pgu8MiIOj4hZdCtqLuy3X2xQRLxkXDbdR9cIrS/b1uUnwBP62rPpch2AzBwDTgc+2J/QZdOIeFo/xnv6OuOP/0/A2yPiCf3YtouIl/SPfaOv82v9aqY/ptsPub6f67cjYqd+DPf3d4/147qg//5vGff7S0QcFRFPjYjN6RrfhzfwXvwPsDS6E/bM6X+2QyPiqI28X9LIsUHUpGXm++hm6T5A15hdSPcJ2vEb2E/wyGNcBbyWbu39ArrmbCFdo/ZovQn4LWAp3Qzg5zf89J/ZBPi/dJ+4Lqbbt/iz5jIzbwcuoQvJ88e97gDgvH7MFwD/mJnf7R/7G+BP+mUmb+qPsXaz/T1079Ob8f9BSdI42Z0E7dPAOyeQHafQ7ZO7li47X98f4zzgT4Ev0mXrfvziXrgNOQq4MCKW0W2peF2/nw66Zu9TfbadtJ7xX0d3wpXzgOuBR+5hfBNwBXARXea+F9ik32LyV8AP++MfnZln9Y9/rt86ciXdKiUycxHwEuA9dMtnDwB+uIGf6wTgqv7n+jDw0kfsmfw08ES6pnytbel+n7iPbsnuvXTbSwA+ARzSj/XL/VLRF9Lt5byZbrb043Qnt5GqEus/p4Y0/folqvfTLRO9ufBwfiYiTgfuzMw/KT0WSZI0WBHxO8CpmflLpccileYFTVVcRLyIbs9i0M1GXgHcUnJM4/VLZ38NzyIqSVJzImJLusti/WPpsUijYKPL2yJir4j4bkRc3Z/N6nX9/X8WEXdExGX9l5twNVkn0i3vvJNuichLN3C5iGkVEe+mW9Ly/lGa0ZQ0s5nN0mBEd53Ce+j2FP574eFII2GjS0yju87Obpl5SURsA/wY+FW6i4Uuy8wPbOj1kiRpsMxmSdKwbHSJaX+WrAX990sj4hp+8To6kiRpGpnNkqRheVRnUOz3Yj2Z7myVAK+JiMsj4vSI2GHQg5MkSRtmNkuSBmnCZzHtzy75X8BfZeaXImIXulP4JvBuuqUur1rH604FTgWIWZsfscWeOw1q7I9e4V1tEeW31WVO9hKDA6o/VrY+wJxZK4vWX3HfrKL1AQ7c9e6i9W++eeei9QEO2Pee0kPgkstXLsrMgf2j+LxnbpX3Ll4z6df/+PIV52TmCYMaj4ZvINk8e/MjZu0xdxpH/YtK55LZbDYDrLi/fDYftEvZbL7xlvLZfOA+i0oPgUsuXzHjs3lCDWJ/gdCvA+dk5gfX8fg+wNcz89ANHWf2/nvkPu/7/UkOdepWry57ybnNN5/8X45BWTNW9j1YuXSLovUBDtv/9qL1b/rSfkXrA5z7xvdv/ElD9Fuv+OOi9QG++el5pYfA7N1v/nFmHjmo4x1x2Oy88Jw9N/7E9dh8txsHOh4N16Cyec7+u+e+Hzh1OIOcgNWrNy1WG0Yjm1cVfg9WPjgC2fy4stl841fKZ/N331B26/BJryyfzed+6l9KD2HgWVhjNm90D2JEBN3FQK8ZH0ARsVu/BwLgxXRnepQkFZOsybHSg9A0MJslqRb1ZfNEroN4LHAKcEVEXNbf9w7g5Ig4nG4Zyy1AualBSRIJjJVeS6/pYjZLUgVqzOaJnMX0B3QXMH+kswc/HEmStDFmsyRpWCYygyhJqsQYdS1jkSSpdbVlsw2iJDUiSdZM8MzUkiRp+GrMZhtESWpIbfscJElqXW3ZXPaaB5IkSZKkkeEMoiQ1IoE1lX1KKUlSy2rMZhtESWpIbctYJElqXW3ZbIMoSY1IqG4jvCRJLasxm20QJakhdZ1IW5Kk9tWWzTOqQdxyzsqi9Zc/tEXR+gCr75lTtP6Lj72oaH2A7378qUXrf/hN/1S0PsBJf/CGovW/cMaHitYHeMZb3lR6CMAbSw9AKm7O7LLZ/ODyWUXrA6y+d3bR+i9+2sVF6wP85+lls/lDb5hXtD7Ar/3h64vW/4/Ty2fzMW81m0fBjGoQJallSVa3EV6SpJbVmM02iJLUioQ1dWWQJEltG3I2R8RewKeBXbpqzMvMD0/lmDaIktSIpL59DpIktWwasnk18MbMvCQitgF+HBHnZubVkz2gDaIkNSNYQ5QehCRJ+pnhZnNmLgAW9N8vjYhrgD0AG0RJkiRJaszciBh/Jql5mbnOsypFxD7Ak4ELp1LQBlGSGpHAmHsQJUkaGQPI5kWZeeTGnhQRWwNfBF6fmQ9MpaANoiQ1xCWmkiSNlmFnc0RsTtccfiYzvzTV49kgSlIjEhtESZJGybCzOSIC+ARwTWZ+cBDH3GQQB5EkSZIkTbtjgVOAZ0XEZf3XC6ZyQGcQJakhY+kMoiRJo2SY2ZyZP4DBTlHaIEpSI6ZhGcvAL8YrSVLLatz+YYMoSY1IgjXD3Tkw8IvxSpLUsmnI5oGbUQ3iA3duW7T+bvsuKlof4O675hStf9OyuUXrA9x3+Oqi9d976/OL1ge466ll/9f/wKKnF60PcNczy/49AODfBn/IIS9jGfjFeKUlCwpn8973Fq0PcNfdZbP5lgcfU7Q+wJInrypa/wO3Pa9ofYC7n1I2m//u3mOK1gdY+Kyyfw+Abp3MgNW2/aOudlaSNBIGdTFeSZI0WmbUDKIktWwA+xzmRsTF427Py8x5j3zSIC/GK0lSy9yDKEkqKFiTU1oYsigzj9xghQFfjFeSpLZNOZunnQ2iJDUigbEh7hwYxsV4JUlq2bCzeRhsECWpIUNexrL2YrxXRMRl/X3vyMyzh1lUkqSaucRUktSkYVyMV5IkjRYbRElqRGZ9+xwkSWpZjdlsgyhJDRlzgk+SpJFSWzbbIEpSI7pTadf1KaUkSS2rMZvrGq0kSZIkaWicQZSkZtS3z0GSpLbVl802iJLUiBqvtSRJUstqzGYbxGm06JJdSg+Bg4+5uWj9yy/dt2h9gE3q2ic8FLMXla3/9G1+WnYAwHd/dHTpIXDbEI65Jv0LrsqMlf07e8+l5bP50GNuKlr/skseV7Q+QBT+p2uTyLIDoHw2P2Oba8sOADj3h8eWHgK3DuGYtWWzDaIkNSKJ6jbCS5LUshqzua7RSpIkSZKGxhlESWrIWGUb4SVJal1t2WyDKEmNqPFaS5IktazGbLZBlKRGJFHdRnhJklpWYzbX1c5KkiRJkobGGURJakht11qSJKl1tWWzDaIkNSIT1lS2EV6SpJbVmM02iJLUjGCMuvY5SJLUtvqy2QZRkhqR1PcppSRJLasxm+sarSRJkiRpaJxBlKSG1HatJUmSWldbNs+oBvGoJ95YtP5FP923aH2Aq655bNH6z37aFUXrA5x3zeOL1n/OztcUrQ9w5j1l/x4cssXCovUB1swqPYLBS4Kxyq61JD31STcUrX/hCGTzFYWz+QXHXlq0PsC3rj2kaP3n7lQ+mz+/cO+i9Q/Z/N6i9QFWzyk9gsGrMZtnVIMoSa2r7VNKSZJaV1s22yBKUiMSGKtsI7wkSS2rMZvrGq0kSZIkaWicQZSkZgRrKrvWkiRJbasvm20QJakRNS5jkSSpZTVmsw2iJDWktk8pJUlqXW3ZXFc7K0mSJEkaGmcQJakRmVHdMhZJklpWYzZvdLQRsVdEfDciro6IqyLidf39O0bEuRFxff/fHYY/XEnShqzJTSb9pXqYzZJUj9qyeSJVVwNvzMxDgKOBP4qIQ4C3Ad/JzAOA7/S3JUmFJDBGTPpLVTGbJakCNWbzRpeYZuYCYEH//dKIuAbYAzgROK5/2qeA7wFvHcooJUkTEM4EzhBmsyTVor5sflSjjYh9gCcDFwK79AEFcBewy2CHJkmSNsZsliQN0oRPUhMRWwNfBF6fmQ9E/HzKMzMzInI9rzsVOBVgs7nbTW20U7To4S2L1t/kgfLnBNr3SfOL1j//m4cVrQ/w0hN/WLT+R773nKL1ATY/cXnR+id87k1F6wNs/sylpYcAHxns4bprLblUdCYZSDbvVDabFz60VdH6sbR8Nj/u0DuK1j/v7COK1gc4+VfOL1r/w+eXz+bNfqVsNh//+TcXrQ+w2S+PQDb//WAPV2M2T+hfxYjYnC6APpOZX+rvvjsidsvMBRGxG7BwXa/NzHnAPIDZ+++xzqCSJA3GGq9eNGMMKpvn7L+72SxJQ1RbNk/kLKYBfAK4JjM/OO6hrwIv779/OfCVwQ9PkjRRSTCWk/9SPcxmSapDjdk8kRnEY4FTgCsi4rL+vncA7wG+EBG/C9wKnDSUEUqSJmyssk8pNWlmsyRVorZsnshZTH8A6z3H6vGDHY4kSdoYs1mSNCzld2ZLkgYiE9a4VFSSpJFRYzbbIEpSQ9xLKEnSaKktm20QJakR3Ub4uvY5SJLUshqzua7RSpIkSZIAiIjTI2JhRFw5qGPaIEpSQ9YQk/6SJEmDN+RsPgM4YZDjdYmpJDUiqW+fgyRJLRt2Nmfm9yNin0Ee0wZRkppR3z4HSZLaVl822yBKUkPGXCoqSdJImWI2z42Ii8fdnpeZ86Y4pA2aUQ3iXUu2LVo/N82i9QGWrphdtP6agx4sWh/grBueVLT+JqvK/wJ/5F63Fa1/0c0HF60PsPN2S0sPgZ+WHsCjFBGnAy8EFmbmoaXHozbc/cA2ReuPQjY/uHJW2QEcvKxsfeDMGw4vWj9WlJ/heeretxatf+Gt5bN5lxHI5utKD+D/tygzj5zOgjOqQZSklk3DxXjPAD4CfHqYRSRJasU0ZPPA2SBKUkOGuc9hGBvhJUlq3TCzOSI+CxxHtxR1PvCuzPzEVI5pgyhJjeguxlvXp5SSJLVs2NmcmScP+pg2iJLUkNo2wkuS1LraTiBngyhJWmvaN8JLkqTRYoMoSY0Y9sV4JUnSo1NjNtsgSlJDatsIL0lS64aZzcNggyhJrcj6NsJLktS0IWfzMNTVzkqSJEmShsYZRElqRFLfmdIkSWpZjdlsgyhJDaltGYskSa2rLZttECWpETWeKU2SpJbVmM0zqkF8+MEtitafu9+9ResDrFi9adH6Ywu2LFofYPUOq4rWn7PXsqL1Aa5cuFvZAezzYNn6wB2XFX4PhqS2EJIeWjaraP2dH7e4aH2AlWvKZvPq27cqWh9gbMey2bzlXkuL1ge4atGuZQew9/Ky9YE7fmI2jwJPUiNJkiRJAmbYDKIktSyp71TakiS1rMZstkGUpIbUdqY0SZJaV1s22yBKUiuyvn0OkiQ1rcJstkGUpEbUeKY0SZJaVmM2e5IaSZIkSRLgDKIkNaW2TyklSWpdbdlsgyhJjajxTGmSJLWsxmy2QZSkhmRlISRJUutqy2b3IEqSJEmSAGcQJakptV1rSZKk1tWWzTaIktSIrPBaS5IktazGbLZBlKSG1LbPQZKk1tWWzTOqQZy15aqi9e+5fYei9QEOOPDOovWXrS7/Hhy094Ki9a+/5LFF6wNs/fhFRes/ePN2ResDbLLX8tJDGIL6zpQmzd5qZdH6d88vn0uPP+COovXvzx2L1gc4eJ+yv5/89NK9i9YH2OrAstk8dtecovUBNtnTbB4FnqRGkiRJkgTMsBlESWpdbctYJElqXW3ZbIMoSY1I6tsIL0lSy2rMZhtESWpFdmdLkyRJI6LCbLZBlKSG1HatJUmSWldbNnuSGkmSJEkS4AyiJDUjqW8jvCRJLasxm20QJakZ9V1rSZKkttWXzTaIktSQ2jbCS5LUutqy2T2IkiRJkiTAGURJakpt+xwkSWpdbdlsgyhJjcisL4QkSWpZjdk8oxrEhxduWbR+blZ+AfJti3coWj93f7hofYDVY2VXVq/ZZk3R+gBLLt2paP3DfvmGovUBLv3JvqWHMBS1bYSXHlq4VdkBbDZWtj5w6707Fq0fuz9UtD6U/7drzVYjkM2XlM3mI555XdH6ABf9ZL/SQxiK0n+/H60Z1SBKUutq2wgvSVLrastmT1IjSZIkSQKcQZSkptS2z0GSpNbVls02iJLUiCSqCyFJklpWYzbbIEpSQyrb5iBJUvNqy2YbRElqRYWn0pYkqWkVZrMnqZEkSZIkAc4gSlJbalvHIklS6yrL5o3OIEbE6RGxMCKuHHffn0XEHRFxWf/1guEOU5I0EZkx6S/Vw2yWpHrUls0TWWJ6BnDCOu7/UGYe3n+dPdhhSZImI3PyX6rKGZjNklSF2rJ5ow1iZn4fWDwNY5EkSRNgNkuShmUqJ6l5TURc3i9z2WF9T4qIUyPi4oi4eM2SB6dQTpK0IUl9y1g0cI86m1c/sHw6xydJM0qN2TzZk9R8DHg33c/8buBvgVet64mZOQ+YBzB7/z2KLmLado8HSpZn+UNbFK0PsPKhzYvWP2SvBUXrA9y8eMei9WMElvLtdNRdRetfetnjitYH2PaxS0oPYfASsNGbySaVzXP2373ov0rb7VH2/8XlD80qWh9g5cNlzxl42N53FK0PcN2inYrWH4Vs3vXoO4vWv+iy/YrWB9h2r7K/qw9Fhdk8qRnEzLw7M9dk5hjwL8BTBjssSdJk1LbPQYNjNkvSaBp2NkfECRHx04i4ISLeNtXxTqpBjIjdxt18MXDl+p4rSZpGOYWvCRh0CGlwzGZJGlFDzOaI2BT4KPB84BDg5Ig4ZCrD3eiahoj4LHAcMDci5gPvAo6LiMPphn0L8PtTGYQkafSNC6HnAPOBiyLiq5l5ddmRzTxmsySp9xTghsy8CSAiPgecCEw6mzfaIGbmyeu4+xOTLShJGpahb2gfeAhpcsxmSarFlLN5bkRcPO72vH4f+Vp7ALePuz0feOpUCpbdFS1JGqyp7SWc9hCSJKl5U8vmRZl55IBGMiE2iJLUimSqn1JOewhJktS0qWfzxtwB7DXu9p79fZM2lesgSpJmloGHkCRJmpKLgAMiYt+I2AJ4KfDVqRzQGURJaslwL1fxsxCiawxfCvzWUCtKklS7IWZzZq6OiNcA5wCbAqdn5lVTOaYNoiQ1ZXjLWIYRQpIktW+oS0zJzLOBswd1PBtESWrJkC94P+gQkiSpeUPO5kGzQZSkllQWQpIkNa+ybPYkNZIkSZIkYIbNIC5/aIui9VffO7tofYBj/9dPi9a/4EcHF60PsMOB9xatv3rHh4vWB7jj+p2K1p+124NF6wM8uHxW6SEMXgLDPZW2NHAPPVw2m1ctLp/NTz/82qL1f3jhCGTzAYuL1h+bu7xofYDbrt21aP05uy8rWh/K/64+FBVm84xqECWpdVnZMhZJklpXWzbbIEpSSyoLIUmSmldZNtsgSlJLKlvGIklS8yrLZk9SI0mSJEkCnEGUpKZEZctYJElqXW3ZbIMoSa1IqtvnIElS0yrMZhtESWpGVLfPQZKkttWXze5BlCRJkiQBziBKUlsqW8YiSVLzKstmG0RJakllISRJUvMqy2YbRElqSWUhJElS8yrLZhtESWpFUt1GeEmSmlZhNs+oBvHg3e4uWv+qO/ctWh9g2apZResf/ORbitYHuO6HZf8c5hx6X9H6AHMeu7Jo/SXztytaH2D7vZaUHoIk4OBdy2bzFSOQzcvXbFG0/mFH3FS0PsBV/7V/0fpbHra4aH2AOfuW/f1g8e3bF60PsP2eZvMomFENoiS1rraL8UqS1LrastkGUZJaUlkISZLUvMqy2esgSpIkSZIAZxAlqSm1LWORJKl1tWWzM4iSJEmSJMAZRElqS2Wn0pYkqXmVZbMNoiS1IqluI7wkSU2rMJttECWpJZWFkCRJzassm92DKEmSJEkCnEGUpKbUdqY0SZJaV1s22yBKUksqCyFJkppXWTbbIEpSSyoLIUmSmldZNs+oBvHq/9m3aP2Djrq1aH2Aa+7epWj9lQ9uUbQ+wOyDlhat/8C9WxWtD7DF1iuL1n/M3vcVrQ+w5MrHlB7CwEXWt4xFuuKixxWtf/ARtxStD3DVXbsWrf/wsllF6wNseciSovXvv2frovUBZm27omj9nR67uGh9gMVXzS09hIGrMZs9SY0kSZIkCZhhM4iS1LzKLsYrSVLzKstmG0RJaklly1gkSWpeZdlsgyhJDaltn4MkSa2rLZttECWpJZWFkCRJzassmz1JjSRJkiQJcAZRktpR4am0JUlqWoXZbIMoSS2pLIQkSWpeZdlsgyhJLakshCRJal5l2eweREmSJEkS4AyiJDWltn0OkiS1rrZsdgZRkiRJkgTMsBnENdusKVr/voe2LFofYPWKTYvWn3XrFkXrAxz8rNuK1r/0nn2L1gfYec+lRevPv2HnovUBjjj2htJD4MZhHLSyTymlsa3LZvP9K+YUrQ+wckXZX8dGIZufcPwtRetfdM/jitYH2GnbZUXrz79pp6L1AY465vrSQ+CmYRy0smyeUQ2iJDWtwlNpS5LUtAqz2SWmkiRJkiTAGURJaktln1JKktS8yrLZGURJaklO4WsKIuIlEXFVRIxFxJFTO5okSQ0plM2T5QyiJDUiKLrP4Urg14B/LjYCSZJGTOFsnhQbRElqSaEQysxrACKizAAkSRpVlTWILjGVJK01NyIuHvd1aukBSZKkyZns9g9nECWpFVM/lfaizFxvgETEecCu63jotMz8ypQqS5LUorKXuZjU9o+NNogRcTrwQmBhZh7a37cj8HlgH+AW4KTMvO/RjVeSNHBDDKHMfPbwjq5Hw2yWpIpUtv1jIktMzwBOeMR9bwO+k5kHAN/pb0uSSqvsTGmatDMwmyWpDlPL5mnf/rHRGcTM/H5E7POIu08Ejuu//xTwPeCtgxyYJKkeEfFi4B+AnYBvRMRlmfm8wsNqltksSTPGtG//mOwexF0yc0H//V3ALut7Yt/lngqw2dztJllOkjQRpfY5ZOZZwFllqqs3uWzeyWyWpGEaZjYPY/vHlE9Sk5kZsf4fOzPnAfMAZu+/R9FFTLPmPlSyPFttsaJofYCxB3csWv+QZ91QtD7AjV/cv2j9U175/aL1Ac769DOK1v/e695XtD7Ai//6LaWHMBwuFRWPLpvn7L970b81sx9TNptnb7aqaH2ANQ+VPWfgYSOQzdd+8cCi9V/28vOL1gf48r8/vWj9819bPptf9B6zeRRM9jIXd0fEbgD9fxcObkiSpEmZyh6HysJL62Q2S9KoKZjNEfHiiJgPPI1u+8c5E3ndZBvErwIv779/OeDpzSVpBERO/kvVM5slaQSVyubMPCsz98zMWZm5y0TPDbDRBjEiPgtcABwUEfMj4neB9wDPiYjrgWf3tyVJ0jQwmyVJwzKRs5ievJ6Hjh/wWCRJU+VM4IxgNktSRSrL5rK7oiVJA+VSUUmSRktt2WyDKEktqSyEJElqXmXZbIMoSa3wbKSSJI2WCrN5smcxlSRJkiQ1xhlESWpE9F+SJGk01JjNNoiS1JLKlrFIktS8yrLZBlGSGlLbmdIkSWpdbdnsHkRJkiRJEjDDZhDz6m2K1t/zubcWrQ9wI7sVrX//itlF6wM8sP+aovUvX7JH0foAY7PK1r961WPKDgDY9vbVpYcwHJV9SimNXVs4m59dPptvYNei9Uchm5ceVPbf5Mvu37NofYA1hf8Yrlg5Atl8q9k8CmZUgyhJzasshCRJal5l2WyDKEmtyPr2OUiS1LQKs9kGUZJaUlkISZLUvMqy2ZPUSJIkSZIAZxAlqSm1LWORJKl1tWWzDaIktaSyEJIkqXmVZbMNoiQ1pLZPKSVJal1t2eweREmSJEkS4AyiJLUjqW4ZiyRJTaswm20QJakllYWQJEnNqyybbRAlqRFBffscJElqWY3ZbIMoSS2pLIQkSWpeZdk8oxrEY593edH653/viUXrAxxx7A1F61/y4/2K1gfYdp8lRetfdt1ji9YH2P7oxUXrv+HSk4rWB1hxwgico+vrpQcglfe0Z19ZtP4PfvCEovUBjjy6bDb/+NIRyObHls3mq27ao2h9gO2PKpvNb7zsJUXrAzz83BHI5q+VHkB5M6pBlKTWRVb2MaUkSY2rLZttECWpFRWeKU2SpKZVmM02iJLUkNo2wkuS1LrastkGUZJaUlkISZLUvMqyeQR2gkqSJEmSRoEziJLUkNqWsUiS1LrastkGUZJaUlkISZLUvMqy2QZRklqR9X1KKUlS0yrMZvcgSpIkSZIAZxAlqS2VfUopSVLzKstmG0RJakRQ3zIWSZJaVmM22yBKUkuyshSSJKl1lWXzjGoQf3jOk4rW3/1pdxatD3DNwp2L1t/65vLbXj/zq6cXrf+b//zGovUB9n/8vUXrX/79A4rWB9jpsHtKD4Fbh3DM2j6llC4479Ci9fd66h1F60P5bN7mxvLZ/MUTP160/q/Me3PR+gD7HrC4aP0rfrB/0fowGtl8yxCOWVs2l/8XQZIkSZI0EmbUDKIkNS2pbiO8JElNqzCbbRAlqSExVnoEkiRpvNqy2QZRklpS2aeUkiQ1r7JstkGUpIbUthFekqTW1ZbNnqRGkjRlEfH+iLg2Ii6PiLMiYvvSY5IkSY+eDaIktSLprrU02a+pORc4NDOfBFwHvH2qB5QkqXpls3lSXGIqSQ0ptYwlM7897uaPgN8oMxJJkkZLbUtMbRAlqSVTC6G5EXHxuNvzMnPeJI7zKuDzUxqJJEmtsEGUJFVqUWYeub4HI+I8YNd1PHRaZn6lf85pwGrgM8MZoiRJGiYbRElqRDDcZSyZ+ewN1o94BfBC4PjMQhsnJEkaIcPO5mGwQZSkVhTc0B4RJwBvAX45M5cXGYQkSaOmYDZP1oxqEOccvrho/SUPzS5aH2D2FquL1l98wJqi9QFee91Li9Zf/tiyfwYA196zc9H6Bx5zc9H6ANdcsk/pIQxFwU8pPwLMAs6NCIAfZeari41G1ZjzpLLZ/MCK8tk8Z9aqovXvPbB8Np963clF6z+0V/lsvnHxY4rWf8KxNxatD3DFxY8rPYShKJXNEfF+4EXASuBG4JWZef/GXjejGkRJal65s5juX6ayJEkjrtyHt+cCb8/M1RHxXrpLUL11Yy/yOoiSJEmS1JjM/HZmrp0e/xGw50Re5wyiJDWkto3wkiS1bkSyecKXoLJBlKRWJDA2GikkSZIYRDZv8BrFw7gElQ2iJLXE/lCSpNEytWze4DWKh3EJKhtESWrIiCxjkSRJvYJnMZ3UJag8SY0kSZIktecjwDZ0l6C6LCL+aSIvcgZRklpS2cV4JUlqXqFsnuwlqKbUIEbELcBSYA2wekPrYyVJw+cSU5nNkjRaasvmQcwgPjMzFw3gOJKkqUg8SY3WMpslaRRUmM3uQZQkSZIkAVNvEBP4dkT8OCJOHcSAJEmTE0BkTvpLzTCbJWlE1JjNU11i+kuZeUdE7Ex3dpxrM/P745/Qh9OpAJvN3W6K5aZmyfyy9ffY756i9QHuvnRd19GcPqc8//sbf9KQnf33zyha/6hXXle0PsCV3zqwaP2v/eG/F60P8LzjDis9BG4exkHHhnFQVebRZfNOZnNpd/1kl6L1T3ne+UXrA3zjI2Wz+bBTbixaH+Da7+5XtP5Zv/e5ovUBnvcMs3kUTGkGMTPv6P+7EDgLeMo6njMvM4/MzCM33W6rqZSTJG1EbZ9SavAebTZvtu2W0z1ESZpRasvmSTeIEbFVRGyz9nvgucCVgxqYJOlRyil+qXpmsySNmAqzeSpLTHcBzoqItcf598z81kBGJUmSJsNsliRNyaQbxMy8CSi/UFiS1MtiF+PVaDCbJWnU1JfNg7gOoiRpRNR2MV5JklpXWzbbIEpSSyr7lFKSpOZVls02iJLUioSo7FTakiQ1rcJsntJlLiRJkiRJ7XAGUZJaUtkyFkmSmldZNtsgSlJL6sogSZLaV1k22yBKUkOisk8pJUlqXW3Z7B5ESZIkSRIww2YQZ+28vGj9O+/evmh9gCN/6bqi9f/ji88oWh/g7W/5j6L1//m0Xy9aH+CyD/9D0fpHv/U1ResDLPvSA6WHAC8+c/DHrOxTSql0Ni9YtF3R+gBHHVs2m7/w5fLZ/I43lc3mf3pn+Wy+5kMfK1r/6Le+umh9gAfOXFZ6CPDrZvOMahAlqWkJVHYqbUmSmlZhNtsgSlIjgqxun4MkSS2rMZttECWpJZWFkCRJzassmz1JjSRJkiQJcAZRktpS2aeUkiQ1r7JstkGUpFZUuBFekqSmVZjNNoiS1JDaNsJLktS62rLZBlGSWlJZCEmS1LzKstmT1EiSJEmSAGcQJakhWd2nlJIkta2+bLZBlKRWJNWFkCRJTaswm20QJakllZ0pTZKk5lWWzTOqQVyxaE7R+gc/fn7R+gAXXb5f0fpbPHFZ0foAlyzbu2j9LR5YXbQ+wGrWFK2/9Z0ri9YHePMTv1Z6CJxUegDSCFhxb9lsfsJBtxetD3Dh5fsXrT/7CUuL1ofy2TzrvvLZvCrLjmEksvmwr5QegtnMDGsQJal1tZ1KW5Kk1tWWzTaIktSSykJIkqTmVZbNNoiS1IoExuoKIUmSmlZhNtsgSlIz6juVtiRJbasvmzcpPQBJkiRJ0mhwBlGSWlLZp5SSJDWvsmy2QZSkllQWQpIkNa+ybLZBlKRWFNwIHxHvBk6kuxzwQuAVmXlnkcFIkjQqKjxJjXsQJakZCTk2+a+peX9mPikzDwe+Drxzyj+OJEnVK5rNk2KDKEmassx8YNzNreg+M5UkSZVxiakktWRq+xzmRsTF427Py8x5E31xRPwV8DvAEuCZUxmIJEnNcA+iJKmIqe9zWJSZR67vwYg4D9h1HQ+dlplfyczTgNMi4u3Aa4B3TWUwkiRVr8I9iDOqQYwVZVfUrlyzadH6ANtcV3YMbzn+G0XrA7zr7JOK1t/81GVF6wM84dw/LFp/1h+sKFof4K+ufUHpIQCXDP6QQ/yUMjOfPcGnfgY4GxtETcSqKFp++erNi9YH2Ob6stn89hHI5tO+WTibf2950foAB513atH6s/5gZdH6AO+57oTSQ6C2bB4G9yBKkqYsIg4Yd/NE4NpSY5EkSZM3o2YQJal55T6lfE9EHER3mYtbgVeXGogkSSOlshlEG0RJakYWC6HM/PUihSVJGmnlsnmybBAlqRUJjJW5ZpIkSVqHCrPZBlGSWlLZp5SSJDWvsmz2JDWSJEmSJMAZRElqS2WfUkqS1LxC2RwR76Y7s/gYsBB4RWbeubHXOYMoSc3I7mK8k/2SJEkDVjSb35+ZT8rMw4GvA++cyIucQZSkViRk1rURXpKkphXM5sx8YNzNrbrRbJwNoiRJkiSNprkRcfG42/Myc95EXxwRfwX8DrAEeOZEXmODKEktcamoJEmjZWrZvCgzj1zfgxFxHrDrOh46LTO/kpmnAadFxNuB1wDv2lhBG0RJaoknqZEkabQMMZsz89kTfOpngLOxQfxFf/q8s4rWf/dF/7tofYBjfvPqovU/cdsvFa0PMLb9qqL1X7TfVUXrA5z5w6OK1n945Zyi9QFe9/Tvlh4CfzToA2ZWdzFe6c+f86Wy9S9+YdH6AE8/6cqi9T85/5ii9QFyu9VF6//WQRdv/ElDdsZ/PaNo/YdXlT935VsP/XbpIfC7gz5gwWyOiAMy8/r+5onAtRN53YxqECWpec4gSpI0Wspl83si4iC6y1zcCrx6Ii+yQZQkSZKkxmTmr0/mdTaIktSQdImpJEkjpbZstkGUpGakS0wlSRop9WWzDaIktSLxMheSJI2SCrPZBlGSWpJ1LWORJKl5lWVz+fPZSpIkSZJGgjOIktSIBLKyZSySJLWsxmye0gxiRJwQET+NiBsi4m2DGpQkaRIyu2Usk/1SE8xmSRohFWbzpGcQI2JT4KPAc4D5wEUR8dXMvHpQg5MkPTq1fUqpwTKbJWn01JbNU5lBfApwQ2belJkrgc8BJw5mWJIkaRLMZknSlExlD+IewO3jbs8Hnjq14UiSpsSlojOd2SxJo6aybI6c5IUbI+I3gBMy8//0t08BnpqZr3nE804FTu1vHgpcOfnhNmEusKj0IArzPfA9AN8DgIMyc5tBHSwivkX3vk7Wosw8YVDj0fQzmyfNf498D8D3AHwPwGye0gziHcBe427v2d/3CzJzHjAPICIuzswjp1Czer4HvgfgewC+B9C9B4M8ns2dMJsnxffA9wB8D8D3AMxmmNoexIuAAyJi34jYAngp8NXBDEuSJE2C2SxJmpJJzyBm5uqIeA1wDrApcHpmXjWwkUmSpEfFbJYkTdVUlpiSmWcDZz+Kl8ybSr1G+B74HoDvAfgegO+BhsBsnhTfA98D8D0A3wPwPZj8SWokSZIkSW2Zyh5ESZIkSVJDpqVBjIgTIuKnEXFDRLxtOmqOkojYKyK+GxFXR8RVEfG60mMqJSI2jYhLI+LrpcdSQkRsHxFnRsS1EXFNRDyt9JimW0S8of//4MqI+GxEzC49pmGLiNMjYmFEXDnuvh0j4tyIuL7/7w4lx6iZx2w2m9cym81ms/ln95nNTEODGBGbAh8Fng8cApwcEYcMu+6IWQ28MTMPAY4G/mgGvgdrvQ64pvQgCvow8K3MfDxwGDPsvYiIPYA/Bo7MzEPpTqLx0rKjmhZnAI88zfXbgO9k5gHAd/rb0rQwmwGzeTyz2Ww2mztmM9Mzg/gU4IbMvCkzVwKfA06chrojIzMXZOYl/fdL6f7h2aPsqKZfROwJ/G/g46XHUkJEbAc8A/gEQGauzMz7iw6qjM2AORGxGbAlcGfh8QxdZn4fWPyIu08EPtV//yngV6dzTJrxzGazGTCbzeafMZs7ZjPT0yDuAdw+7vZ8ZuA/wGtFxD7Ak4ELCw+lhL8D3gKMFR5HKfsC9wCf7JfyfDwitio9qOmUmXcAHwBuAxYASzLz22VHVcwumbmg//4uYJeSg9GMYzaPYzabzZjNZnPHbMaT1EyriNga+CLw+sx8oPR4plNEvBBYmJk/Lj2WgjYD/hfwscx8MvAgM2zpQr+W/0S6QN4d2CoifrvsqMrL7nTSnlJaKsBsNpsxm83mdZjJ2TwdDeIdwF7jbu/Z3zejRMTmdAH0mcz8UunxFHAs8CsRcQvdUqZnRcS/lR3StJsPzM/MtZ9Qn0kXSjPJs4GbM/OezFwFfAk4pvCYSrk7InYD6P+7sPB4NLOYzZjNmM1gNoPZPJ7ZzPQ0iBcBB0TEvhGxBd2m169OQ92RERFBt7b9msz8YOnxlJCZb8/MPTNzH7q/A/+ZmTPq06nMvAu4PSIO6u86Hri64JBKuA04OiK27P+/OJ4ZdjKAcb4KvLz//uXAVwqORTOP2Ww2m82YzT2z+efMZrpp9aHKzNUR8RrgHLqzIp2emVcNu+6IORY4BbgiIi7r73tHZp5dbkgq5LXAZ/pfyG4CXll4PNMqMy+MiDOBS+jOIHgpMK/sqIYvIj4LHAfMjYj5wLuA9wBfiIjfBW4FTio3Qs00ZjNgNuvnzGaz2WweJ7rltZIkSZKkmc6T1EiSJEmSABtESZIkSVLPBlGSJEmSBNggSpIkSZJ6NoiSJEmSJMAGUZIkSZLUs0GUJEmSJAE2iJIkSZKk3v8DvJyijpp49j0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(16, 6))\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.title('Original system')\n",
    "plt.pcolor(s['snapshots'].real)\n",
    "plt.colorbar()\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.title('Reconstructed system')\n",
    "plt.pcolor(dmdc.reconstructed_data().real)\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, we built the approximation for $\\mathbf{A}$ and for $\\mathbf{B}$; we can now test the system with a different control input: differently by the other versions, we can pass as argument of the `reconstructed_data` method the control input we want to test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcMAAAFpCAYAAADgAUNWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa9UlEQVR4nO3df7Cc1X3f8c9HVz9AAlkggSyQQFBkEswYhBUMwSamkBQcN9ip45pOKIlpZLsmAZeOi3GmydTNjDNO3aRO6s5NIJAJxfXYYCjBxkCxqacGW/wwSAgDBgESEkIIIRBCQvd++8cu6/sQXe3V3bN7nnOf92tmR3f37j770SLu937POc95HBECAKDJpuUOAABAbhRDAEDjUQwBAI1HMQQANB7FEADQeBRDAEDjdS2GtpfYvsv2I7bX2L60/fgf295g+8H27QP9jwsAQHrudp6h7UWSFkXE/bYPlnSfpA9J+qikVyPiz/qeEgCAPpre7QkRsVHSxvbXr9heK+nIfgcDAGBQ9mvO0PZSScsl3dt+6BLbD9m+2vYhqcMBADAIXYdJO0+0D5L0fUl/EhE32F4oaYukkPQFtYZSP76X162UtFKSPGvmu2cccViq7FPK0NBo7ggVI6POHaG2Zk3fkztCxRuvzMwdoeMdhz2fO0LFk08fnjtCxfFHb8kdoeK+h3ZtiYjkP5T/2Vlz4sWtIz0d476Hdt0WEecmitTVhIqh7RmSbpF0W0R8eS/fXyrplog4cV/HmXXs4jjiP316klETq9nP+rfNey13hIpXXj0gd4TaOnZhvX6grb/zqNwROu785JdyR6j42Kc+kztCxfeG/zp3hIqhRY/fFxErUh/33ScdEPfetrinY8xY9LO+ZBtP1zlD25Z0laS1Ywuh7UXt+URJ+rCk1f2JCAAoS2gk6jXa1U3XYijpDEkXSnrY9oPtx66UdIHtk9UaJl0n6RN9yAcAKExIGlVZV0SayGrSH2jvg4q3po8DAMDgTaQzBABgv4xq6g2TAgAwYaHQSGEXjqcYAgCSK23OkI26AQCNR2cIAEgqJI0U1hlSDAEAyZU2TEoxBAAkFRILaAAAKOvEigYXwwPn7ModoWLbc3NzR6h6oz6bty4/6ancESqeuv643BEqLv+DG3JH6Ljg9y7LHaHiK//jv+WOUHHGZy7LHeEt/n3uALXR2GIIAOiPULCABgDQcCGNlFULKYYAgLRaG3WXhWIIAEjMGqnbRWO7YAcaAEDj0RkCAJIKSaPMGQIAmq60YVKKIQAgqdbepGUVQ+YMAQCNR2cIAEhuNMrqDCmGAICkShwmpRgCAJIKWSOFzcI1thju2Do7d4SKhUteyh2hYvNTh+aO0DHN9Vqjve099drk/cZNy3NH6NhwVr1+pNy4/ZTcESo2fXB37ghVX+vfoUsbJi2rdAMA0Af1+jUOAFA85gwBAJA1EmUNPFIMAQBJta5aQTEEADRcacOkZZVuAEDj2V5i+y7bj9heY/vSXo9JZwgASCqi73OGeyRdHhH32z5Y0n22b4+IRyZ7QIohACC50T4Ok0bERkkb21+/YnutpCMlUQwBAPXQOrViMLNwtpdKWi7p3l6OQzEEANTRAturxtwfjojhsU+wfZCkb0q6LCK29/JmFEMAQGJJ5gy3RMSKcd/BnqFWIbwuIm7o9c0ohgCApPp9nqFtS7pK0tqI+HKKYza2GA5tr9dfffvTh+WOUHHUmc/ljtDx8F3LckeomDa7XhuHb5w/N3eEjrk/q9e5ZYfP6GnkLLmFt87MHaFiXR+PPdLfjbrPkHShpIdtP9h+7MqIuHWyB6xXRQAAFK/fl3CKiB9IaZerctI9AKDx6AwBAMmNslE3AKDJBnmeYSoUQwBAUiH3ewFNcmWVbgAA+oDOEACQHNczBAA0WoS40j0AoOnc16tW9APFEACQVKi8zrCstAAA9AGdIQAgOc4zLMQ73/1U7ggVDz98dO4IFc88vjB3hI5fOuux3BEq7lu3JHeEil9fsiZ3hI7bX31f7ggV75v9eO4IFcOHlFUgJitkjRZ2nmFjiyEAoH/oDAEAjRYqb2/SstICANAHdIYAgMSsEc4zBAA0WYnDpBRDAEBypXWGZZVuAAD6gM4QAJBUhIsbJu2a1vYS23fZfsT2GtuXth8/1Pbtth9v/3lI/+MCAEowEtN6ug3aRN5xj6TLI+IESadJ+rTtEyRdIenOiFgm6c72fQBAw4Wk0faVKyZ7G7Suw6QRsVHSxvbXr9heK+lISedLen/7addK+p6k/9CXlACAgnhqX7XC9lJJyyXdK2lhu1BK0iZJ9dnMEgCA/TDhBTS2D5L0TUmXRcR2++dtbESE7RjndSslrZSkofnzegqb0uqnj8gdoeLATUO5I1TsXLYrd4SOx772jtwRKo758LO5I1R8/Vtn5o7Q8cYv7fXHQDYXfOXy3BEqdvzCaO4IA9E6z7CsUysmVAxtz1CrEF4XETe0H37e9qKI2Gh7kaTNe3ttRAxLGpakWccurtf/KQCAvihto+6JrCa1pKskrY2IL4/51s2SLmp/fZGkm9LHAwCU5s1LOPVyG7SJdIZnSLpQ0sO2H2w/dqWkL0r6uu2LJT0t6aN9SQgAKM5oYZ3hRFaT/kAad53r2WnjAAAweOxAAwBIKkIamYoLaAAA2B9TcjUpAAAT1VpAU9acYVlpAQDoAzpDAEBypV3PkGIIAEhqyu5AAwDAxJU3Z0gxBAAkl+MyTL1obDH0tHptk7prfr3yzJn7eu4IHa+eMZI7QsX25w7LHaHCB9fn384vLl+XO0LFY68dkztCxYJ/8mLuCBXrcgeokcYWQwBAf3DSPQAAEnOGAIBme/OqFSWhGAIAkittAU1ZfSwAAH1AZwgASIqT7gEAEAtoAABNF+UtoCmrdAMA0Ad0hgCApELlrSalGAIAkmOYFADQaG+uJu3l1o3tq21vtr06RebGdoYjO+v1V3/bsdtyR6iIGv1W56dn545QEYe9kTtCxdARr+WO0PHMtkNyR6jYs6Q+G85L0kurF+SOMDAD6AyvkfSXkv4uxcHoDAEAxYmIuyVtTXW8erVHAIDiJdqbdIHtVWPuD0fEcK8HHQ/FEACQXILVpFsiYkWKLBNBMQQApBXlrSalGAIAkipxb1IW0AAAimP7ekk/lHS87fW2L+7leHSGAIDk+t0ZRsQFKY9HMQQAJMWV7gEAUL027pgI5gwBAI1HZwgASI6rVgAAGi04zxAAgPLmDBtbDD1Sr/9QO346L3eEiqFjduSO0DHj1dwJqmYs25k7QkX8ZG7uCB2vLqvXVSJm/eyA3BEqXj98JHeEASlvNSkLaAAAjdfYzhAA0D8MkwIAGq3EvUkphgCAtKK1orQkFEMAQHKlnWfIAhoAQOPRGQIAkgqxgAYA0HjlnWdIMQQAJFfaAhrmDAEAjUdnCABIjjlDAECjRVAMizHjxXr91WMod4Kq3VtrtMHxonptbjxnWr0mQ3YdmDvBzx30QI3+3UiaefaW3BEqRu9bkDvCwLCABgDQeCygAQCgMHSGAIDkmDMEADRayBRDAAAKmzKkGAIAEivw1AoW0AAAGo/OEACQXmHjpF07Q9tX295se/WYx/7Y9gbbD7ZvH+hvTABASSLc023QJjJMeo2kc/fy+H+NiJPbt1vTxgIAlKy1Jdvkb4PWtRhGxN2Stg4gCwAAWfSygOYS2w+1h1EPGe9JtlfaXmV71cj2HT28HQCgBG9e6b6kYdLJLqD5qqQvqPV3/oKk/yLp43t7YkQMSxqWpFnHLq7NlOroMTtzR6gY2VmztUx76rMseuk7NuWOUPHM84fmjlARc0ZzR+iY9s5tuSNUbFs7P3eEimnHN6QhCEmFnVoxqZ/AEfH8m1/b/mtJtyRLBAAoXiM26ra9aMzdD0taPd5zAQANFD3eBqxrZ2j7eknvl7TA9npJfyTp/bZPVivyOkmf6F9EAAD6q2sxjIgL9vLwVX3IAgCYEtioGwCA4nagoRgCANJio24AAMpDZwgASI9hUgAAyhompRgCANKjMwQANF5hxZAFNACAxmtsZziybWbuCBXTX6nX7yULT36++5MG5LkfHpk7QoWXvZY7QsXo3DdyR+jY+ZN6bWI+ctSu3BEqRl6alTvCYDRlo24AAPaltI26KYYAgPQKK4b1GpsDAEwN4d5uXdg+1/ZPbT9h+4pe41IMAQBFsT0k6a8knSfpBEkX2D6hl2NSDAEAyTl6u3VxqqQnIuLJiNgt6WuSzu8lL8UQAJBWrxf27V4Mj5T07Jj769uPTRoLaAAAiU1s3q+LBbZXjbk/HBHDvR50PBRDAEAdbYmIFeN8b4OkJWPuL24/NmkMkwIA0uvvMOmPJS2zfYztmZI+JunmXuLSGQIA0uvjeYYRscf2JZJukzQk6eqIWNPLMSmGAID0+nzSfUTcKunWVMejGAIA0mJv0nIsOnZL7ggVW+5fmDtCxewZu3NH6Dj6vc/kjlDx7F1H5Y5QMfvd23JH6Nh9UH3+3UiSNs3OnaDigLfXa5N3/FxjiyEAoH8mcOJ8rVAMAQDpFVYMObUCANB4dIYAgORKGyalMwQANB6dIQAgPU6tAAA02sS2VKsViiEAIL3CiiFzhgCAxqMzBAAkV9pqUoohACA9iiEAoPEohmXYem+9Nsaef+rm3BEqfvbcYbkjdIxun5k7QsX0d+7IHaHi9Wfn5o7QMTR/V+4IFbNqtjH26NqDc0cYCEd5w6QsoAEANF5jO0MAQB9x0j0AoPEKGyalGAIAkittzpBiCABIr7BiyAIaAEDj0RkCANIq8NQKiiEAID2KIQCg8QorhswZAgAaj84QAJBcaXOGdIYAgMZrbGe46/CR3BEq5sys1wbHnlafDYVnvDSUO0LFipOfyR2h4oc7jssdoeOow1/MHaHiyXX12pD/zHPW5I5Q8fgf9vHghXWGjS2GAIA+KfDUCoZJAQCNR2cIAEivsM6QYggASI9iCABoMqu8OUOKIQAgvcKKIQtoAACNR2cIAEhrKp5aYftq25ttrx7z2KG2b7f9ePvPQ/obEwBQlOjxNmATGSa9RtK5b3nsCkl3RsQySXe27wMA0DLVimFE3C1p61sePl/Ste2vr5X0obSxAAAYnMnOGS6MiI3trzdJGncDQNsrJa2UpKH58yb5dgCAkpQ2Z9jzApqICHv8v3ZEDEsalqRZxy6uzcdz6JJtuSNUzJ7+Ru4IFTNm7skdoeM9Z6/u/qQBuue2E3NHqLj0I9/OHaHjv3/rvNwRKp763a/mjlBx4lc+lTvC4NTmp/3ETPbUiudtL5Kk9p+b00UCABSt1/nCOs4ZjuNmSRe1v75I0k1p4gAApgJHb7dBm8ipFddL+qGk422vt32xpC9K+lXbj0s6p30fAIAidZ0zjIgLxvnW2YmzAACmisLmDNmBBgCQXONWkwIA8I9QDAEAjZZpRWgvuGoFAKDxKIYAgKSc4NbT+9u/ZXuN7VHbKybyGoohACC9vCfdr5b0m5LunugLmDMEACSXczVpRKyVJHviPSbFEABQRwtsrxpzf7i913VfNLYY7rxnfu4IFVveW5+NsSXp9Y1zckfouEdH545QsXveaO4IFX/7xGm5I3TM3J47QdWfbj0ud4SKt/9oV+4IFWv7efDeO8MtETHufJ/tOyS9fS/f+nxE7PcWoY0thgCAPurzMGlEnJPyeBRDAEBamTbb7gWrSQEA6WVcTWr7w7bXSzpd0j/Yvq3ba+gMAQBTSkTcKOnG/XkNxRAAkFxpw6QUQwBAehRDAEDTldYZsoAGANB4dIYAgLQKvIQTxRAAkB7FEADQZFZ5c4YUQwBAehTDMiw9Z13uCBWP/ahem1G/89R1uSN0PPKTen02By55JXeEipdfOCh3hI65p2/LHaHi2kfrs4m5JO3+5zNyR6i6I3eA+mhsMQQA9I+jrNaQYggASIvVpAAAsIAGAIDiOkN2oAEANB6dIQAgOYZJAQCgGAIAGi3K6wyZMwQANB6dIQAgvcI6Q4ohACApNuoGAECS2I6tDOtvWpo7QsUBZ76cO0LF6seW5I7QseChek1tX/GB/507QsV//ovfzh2hY/p5e3JHqBj9f/NyR6jw8tdyRxiY0jrDev2UAQAgg8Z2hgCAPmGjbgAAJI/mTrB/KIYAgPToDAEATccCGgAACkNnCABIK8R5hgAAlDZMSjEEAKRXWDFkzhAA0Hh0hgCApNioGwCACBbQlGK0Zhtjzz3w9dwRKuYdvTN3hI4NQ/NzR6j4n5tOyx2h4uVTdueO0LFoaCR3hIoTf3N17ggV37//F3NHGBg6QwAACiuGLKABADQenSEAIDmGSQEAzRaSRsuqhhRDAEB6ZdVCiiEAIL3ShklZQAMAaDw6QwBAek066d72OkmvSBqRtCciVqQIBQAoW2nDpCk6w7MiYkuC4wAApoJQcQtomDMEADRer8UwJH3X9n22V6YIBAAoW+uqFdHTbdB6HSZ9b0RssH24pNttPxoRd499QrtIrpSkofnzeny7dHZsOih3hIo9h9WrSd+zrj6fz+m//GjuCBU/vbZemy0f9hsv5I7Q8eKPF+aOUPGHF/1d7ggVz522PXeEiqf7efDRfh48vZ5+AkfEhvafmyXdKOnUvTxnOCJWRMSKoblzenk7AEAhSusMJ10Mbc+xffCbX0v6NUn1ul4KAGDwIsGtB7a/ZPtR2w/ZvtH2vG6v6aUzXCjpB7Z/IulHkv4hIr7Tw/EAAEjhdkknRsS7JD0m6XPdXjDpOcOIeFLSSZN9PQBgqsp7pfuI+O6Yu/dI+ki317ADDQAguQQn3S+wvWrM/eGIGJ7EcT4u6X91exLFEACQXu+d4ZZ97Wpm+w5Jb9/Ltz4fETe1n/N5SXskXdftzSiGAIC0QnKfT62IiHP29X3bvyPpg5LOjuhemSmGAIApxfa5kj4r6Vci4rWJvIZiCABIL+9VK/5S0iy1NoORpHsi4pP7egHFEACQXsZaGBHH7e9rKIYAgORy7CLTi3ptiAkAQAbN7Qxr9kvL61sPzB2hYt47XsodoeOhb9VrY+xf+8S9uSNUPPC5U3JH6Pj74S/ljlBx8ac/kztCxTPX1Gz36ou+0b9jF9YZNrcYAgD6I1TcVSsohgCApKw8V57oBcUQAJBeYcWQBTQAgMajMwQApFdYZ0gxBACkxQIaAADKO+meYggASK+wYsgCGgBA49EZAgASi+I6Q4ohACCtEMUQAABWkxZi+stDuSNUHPKuF3JHqNi+6rDcETreOGln7ggVD710RO4IFTO2784doWN7zMwdoeLADTtyR6j4xIof5Y5QcWXuADXS2GIIAOgfTq0AAIBiCABotJA0SjEEADRaeadWcNI9AKDx6AwBAOkV1hlSDAEA6VEMAQCNxgIaAABCirK2oGEBDQCg8egMAQDpMWcIAGg05gzLsWfuSO4IFaNy7ggVs17KneDn/t0pt+WOUPGn3z4/d4SKaZe+ljtCx2/d9ancESqmf/aN3BEqXlh3Su4Ib/Ht/h26sM6QOUMAQOM1tjMEAPRRYZ0hxRAAkFh5e5NSDAEAaYWk0bLOM6QYAgDSK6wzZAENAKDx6AwBAOkV1hlSDAEAiQUn3QMAGi6kYKNuAADKQmcIAEiPYVIAQOOxgKYMv/O+/5s7QsU1D5yeO0LFGf/ykdwROm5+/qTcESpGD6nX5s+Xvut7uSN0/Pkd5+WOUDE6s14b8v/b476fO0LF7/XrwBGcdA8AQGmdIQtoAACNR2cIAEguGCYFADQbV60AADRdiFMrAAAQO9AAAJCP7S/Yfsj2g7a/a/uIbq+hGAIAkgpJMRo93Xr0pYh4V0ScLOkWSf+x2wt6Koa2z7X9U9tP2L6il2MBAKaIiNYwaS+3nt4+to+5O0et+rxPk54ztD0k6a8k/aqk9ZJ+bPvmiKjP1iUAgCwSdHc9sf0nkv61pJclndXt+b10hqdKeiIinoyI3ZK+Jun8Ho4HAMCbFtheNea2cuw3bd9he/VebudLUkR8PiKWSLpO0iXd3qyX1aRHSnp2zP31kt7Tw/EAAFNF76tJt0TEinEPH3HOBI9znaRbJf3Rvp7kmOSJkbY/IunciPg37fsXSnpPRFzyluetlPRmRT9R0upJveHUt0DSltwhaozPZ3x8NvvG5zO+4yPi4NQHtf0dtT73XmyJiHMn+f7LIuLx9te/L+lXIuIj+3pNL53hBklLxtxf3H6sIiKGJQ23Q63aV6VvMj6bfePzGR+fzb7x+YzP9qp+HHeyRSyhL9o+XtKopKclfbLbC3ophj+WtMz2MWoVwY9J+lc9HA8AgJ5FxL/Y39dMuhhGxB7bl0i6TdKQpKsjYs1kjwcAQC49bccWEbeqNTE5UcO9vN8Ux2ezb3w+4+Oz2Tc+n/Hx2bRNegENAABTBduxAQAabyDFkG3bxmd7ie27bD9ie43tS3NnqhvbQ7YfsH1L7ix1Y3ue7W/YftT2Wtun585UF7Y/0/5/arXt620fkDtTTravtr3Z9uoxjx1q+3bbj7f/PCRnxpz6XgzHbNt2nqQTJF1g+4R+v29B9ki6PCJOkHSapE/z+fwjl0pamztETf2FpO9ExC9IOkl8TpIk20dK+gNJKyLiRLUW+X0sb6rsrpH01lMerpB0Z0Qsk3Rn+34jDaIzZNu2fYiIjRFxf/vrV9T6YXZk3lT1YXuxpF+X9De5s9SN7bdJOlPSVZIUEbsjYlvWUPUyXdKBtqdLmi3pucx5soqIuyVtfcvD50u6tv31tZI+NMhMdTKIYri3bdv4Yb8XtpdKWi7p3sxR6uTPJX1WrZNnUXWMpBck/W17GPlvbM/JHaoOImKDpD+T9IykjZJejojv5k1VSwsjYmP7602SFuYMkxMLaGrC9kGSvinpsrdcfqSxbH9Q0uaIuC93lpqaLukUSV+NiOWSdqjBw1xjtee+zlfrF4YjJM2x/dt5U9VbtE4taOzpBYMohhPatq3JbM9QqxBeFxE35M5TI2dI+g3b69QaXv+ntv8+b6RaWS9pfUS8OZLwDbWKI6RzJD0VES9ExBuSbpD0y5kz1dHzthdJUvvPzZnzZDOIYtjZts32TLUmsW8ewPsWwbbVmvNZGxFfzp2nTiLicxGxOCKWqvXv5v9EBL/dt0XEJknPtvdglKSzJXE90ZZnJJ1me3b7/7GzxeKivblZ0kXtry+SdFPGLFn1tAPNRLBtW1dnSLpQ0sO2H2w/dmV7dx+gm9+XdF37F80nJf1u5jy1EBH32v6GpPvVWrH9gBq+24rt6yW9X63rBK5X65JGX5T0ddsXq7Wh9UfzJcyLHWgAAI3HAhoAQONRDAEAjUcxBAA0HsUQANB4FEMAQONRDAEAjUcxBAA0HsUQANB4/x8PMvQTi7hLAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "new_u = np.exp(s['u'])\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.pcolor(dmdc.reconstructed_data(new_u).real)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can also use a different timestep for the reconstruction, obviously passing the proper control input (the number of inputs we pass has to be the number of reconstructed snapshots we try to compute, except for the first snapshots). We continue halving the timestep and create a new random input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAFpCAYAAAD+yHPlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAij0lEQVR4nO3dfZBddZ3n8c83/UTSafIInSbBECAIAeQpw+DDuAwgAjtF1EEGakqDOBVnFnZ1d7ZWGKvQ0tHyabRqLEorI1lx1gFclSXjRCIgyowuDIENJBAwnZhAYkiT585jP9zv/nFP66W53X2S8725555+v6pO5T6c/tzfPTnd3/v7nXN/x9xdAABgbBPq3QAAABoFRRMAgJQomgAApETRBAAgJYomAAApUTQBAEhpzKJpZqea2eNm9qKZvWBmH08e/4yZbTWz1clyXe2bCwBA/dhY39M0sy5JXe7+rJl1SHpG0vsk3Shpv7t/teatBAAgB5rHWsHdt0naltzuNbN1kmbXumEAAOTNUR3TNLPTJF0k6ankodvN7HkzW2Zm06IbBwBAnow5PPu7Fc0mS/qFpM+7+4/MrFPSDkku6XMqD+HeWuXnlkhaIknW1npJS9dJmRs94bBlzhjik0phWU1NcVkD/U0xQZGzJJbitrvFbSp5a1zY9IkHQ3J27ZsckiNJbbvi3t+RGXHn/s2b+npY1pbDU2OCXm+JyZHkY47DpTetc19YVseEw2FZL63p2+Hu2f8oD/PeP273nbsGM2U88/yRle5+TVCTwqQqmmbWIunHkla6+9eqPH+apB+7+3mj5bTNm+Ndn/7Px9jU35vU3Zo5Y8jARfvDsqa2HwrL6tk+JSaoL+6PpB2Oy2o6GJflc+O2+00LVoXk3P/ou0JyJOmM++P20fUfbg/L+u6ffDMs687uD4Tk+LdODsmRpMPTgj64Srrhvz4SlnXV5BfDsi6du/kZd18YFpi45IIT/KmVczJltHRtqEnbshrzs5SZmaR7JK2rLJhm1pUc75Sk90taW5smAgAai2vQA4eTciTNAMQ7JX1I0hozW5089jeSbjazC1UeANwk6WM1aB8AoMG4pFLosaH8SHP27L9JqnYwa0V8cwAAyK/AQ90AAJSVNH6HZwEASM3lGkz5zYxGQ9EEAIQr6jFNJmwHACAlepoAgFAuabCgPU2KJgAgXFGHZymaAIBQLnEiEAAAaRXzCycNWjRbLtsVlnXoN1PDsnYdipuke2JvTE5rUI4k7Tsz7tfg3HduCMt67rl5YVnPvu/0kJwT/ixucvvb7v9RWNZfr7ohLOsL190YlmVnxVwk6W13rQ7JkaTJTUfCsh76wlVhWY93XxaWJX06MGt8aMiiCQDIL5dzIhAAAKm4NFjMmknRBADEKk/YXkwUTQBAMNNg1et8ND5mBAIAICV6mgCAUC6pxDFNAADSKerwLEUTABCqPPdsMYsmxzQBAEiJniYAIFzJi9nTpGgCAEIVeXiWogkACOUyDRb06F9DFs29W6fEhU0ZCIuae972sKzfbJsZkmMvTgzJkaT2zXG/BM0Xx80X0jrrYFjW9vfODsk5ODvu/X15w3vDspqb49q17vbpYVlTXm4Kyfnlb+Mm7z9t6u6wrD1nxvW6DnR2hGXpqbio4Yo6PFvMjwIAANRAQ/Y0AQD5xTFNAABSMw16MQcyKZoAgFDlq5xQNAEASKWow7PF/CgAAEAN0NMEAIRy55gmAACplQo6PEvRBACEKn/lpJg9zWK+KwAAaoCeJgAgGMc0AQBIhe9p5s1A3AHm6b9qDcs6uDdmsm9JmtwVs8PtXRA3IX3btMNhWc8+c0ZY1tSX4n45jwRdC2DinN6YIEmTWvrCsnpenhWWNfuZwbCs1y+KyXn3rFdjgiQ92zMnLGtiT1iUOrbG/U7X0mBBJ2xvzKIJAMitIl8arJjvCgCAGqCnCQAIV+JEIAAAxlbk72lSNAEAoVxW2BOBivlRAACAGqCnCQAIV9TvaRbzXQEA6sZdGvQJmZaxmNkyM+sxs7UjPG9m9vdm1m1mz5vZxRHvjaIJAAhmKmVcUviOpGtGef5aSfOTZYmkb2Z+W2J4FgAQzKWazz3r7k+Y2WmjrLJI0nfd3SU9aWZTzazL3bdleV16mgCAIpotqXJexS3JY5nQ0wQAhAv4nuZMM1tVcX+puy/NGppVQxbNW/7DE2FZK+efE5a1fceJYVml3piJ5Ns3xv0XH5neHpZVmhY36fRgW9yAycnPxkyOvml2R0iOJF1yetXzHI7J+tmdYVn+/1rCsiZlGjD7ve2H4rb7kf7AP49xzdL+Wfn/s+0ylbJ/T3OHuy/M8PNbJZ1acX9O8lgmDM8CAMINakKmJcBySR9OzqK9TNLerMczpQbtaQIA8stV+7lnzew+SZerPIy7RdKnJbVIkrt/S9IKSddJ6pZ0UNJHIl6XogkAaDjufvMYz7uk26Jfl6IJAAhmGkz3XcuGQ9EEAIQ6HsOz9ULRBACEK2pPs5gfBQAAqAF6mgCAUO5W2OHZMd+VmZ1qZo+b2Ytm9oKZfTx5fLqZPWJm65N/p9W+uQCARlDrq5zUS5qWDUj6a3dfIOkySbeZ2QJJd0h6zN3nS3osuQ8AGOdcOh5XOamLMYdnkxkUtiW3e81sncqT3i5S+YulknSvpJ9L+mRNWgkAaCCW695iFkf1rpLLsFwk6SlJnRVTEr0mKW5SSwAAcij1iUBmNlnSDyV9wt33mf2+++zubmY+ws8tUfkCoGqaMTVTY4f8zyffFZIjSU37m8KyOjbGDSkcmBOTc/DsIzFBkpp6YiaRl6RzvrE/LKt/+sSwrA23xHw69r5SSI4kPf7Vd4RlRX6y7bkkbn8f6IyZKH/iN08LyZGkri2Hw7I2fDBuf9g/Py5LNbpmSPl7mvkdYs0iVdE0sxaVC+b33P1HycPbhy7oaWZdknqq/WxyKZelktQ2b07VwgoAKJagSddzJ83ZsybpHknr3P1rFU8tl7Q4ub1Y0kPxzQMANJqhS4NlWfIqTU/znZI+JGmNma1OHvsbSV+U9H0z+6ikzZJurEkLAQANp1TQnmaas2f/TRrx/N8rY5sDAEB+MSMQACCUuzSY4yHWLCiaAIBweT4umQVFEwAQqnwiUDGPaRbzXQEAUAP0NAEA4Yp6PU2KJgAg1LifEQgAgPSKe0yTogkACJfny3tl0ZBF0/rj/jO8KW463EOdce1qPhCT46/HTbJemhU3+ftv/jTumuWRo0DtU/aG5BzYOjkkR5IOBO5Xh7ri9nfvOhSWpf6YCyfsXBDXu9k9f1JY1rQzdoZljXBtjGPySljS+NGQRRMAkF9MbgAAwFHgmCYAACkMXeWkiCiaAIBwRT0RqJj9ZwAAaoCeJgAgFJMbAABwFDgRCACANLy4JwIV86MAAAA1QE8TABDKVdyzZymaAIBwRR2epWgCAEJx9mzONB2IOxTbtidw8vfAI8SH5sdMjv6WU+Imin5t94lhWZHbamBe3MTh/fvbQnImHIl7g/vPGgjLssNx7Zq4ZmJY1slXbg3J2T09rk17X+8IyxpYMyMsa0Lc7lBTRS2anAgEAEBKDdnTBADkF3PPAgBwFDh7FgCANLy4xzQpmgCAUEU+e5YTgQAASImeJgAgXFF7mhRNAEAozp4FAOAoeEGLJsc0AQBIiZ4mACAc39MEACAF53uaAACkV9RjmuO+aJYCt0Dbrris0qaYq21s7jspJEeSJhxsCsvqer4UlnX41bgrW+y74kBIzkD7YEiOJE1fFbeTetx/oXrneVjWbzafHJIzZU1LSI4kzeyNe387L4jb3wda49pVO8U9e5YTgQAASGnc9zQBAPEYngUAIAXmngUAIC0vn0GbZRmLmV1jZi+bWbeZ3VHl+VvM7HUzW50sfxHx1uhpAgDC1fJ7mmbWJOluSe+RtEXS02a23N1fHLbqA+5+e+Rr09MEADSaSyV1u/tGd++TdL+kRcfjhSmaAIBQrvKJQFkWSTPNbFXFsqTiJWZLerXi/pbkseH+1MyeN7MfmNmpEe+N4VkAQLCQ72nucPeFGX7+nyXd5+5HzOxjku6VdEXWRtHTBACEq/GJQFslVfYc5ySPVby+73T3I8ndb0u6JOJ9UTQBAI3maUnzzWyembVKuknS8soVzKyr4u71ktZFvDDDswCAcLWc3MDdB8zsdkkrJTVJWubuL5jZZyWtcvflkv6LmV0vaUDSLkm3RLw2RRMAEKo8xFrbyQ3cfYWkFcMeu6vi9p2S7ox+3YYsmhY4X7HFzaOsQ51xWVHtmrI2bgLrtr1xG37P/LhfqEnb4trVuro9JKd0/qGQHEnad0bcr2n71rHXSWtCf9z/4ayfx7zHfXNDYiRJh2fGvb8T18dlNR0Ze520NsdFvUlRZwRqyKIJAMi3NLP6NCJOBAIAICV6mgCAcFzlBACAFFxG0QQAIK2CHtKkaAIAgh2Hr5zUCycCAQCQEj1NAEC8go7PjtnTNLNlZtZjZmsrHvuMmW2tuCL2dbVtJgCgkQRcGiyX0gzPfkfSNVUe/7q7X5gsK6o8DwAYp2p8lZO6GbNouvsTKk92CwDAuJblRKDbkytiLzOzaSOtZGZLhq68Pdh7IMPLAQAagau4w7PHeiLQNyV9TuVt8zlJfyfp1moruvtSSUslqW3enJBO90Bnf0SMJKl/etx/jg3EnYzcdCAmq3l/3PvrnRuXVQo8BW3X+XFZpdbBkJzWjRNDcqTYiwr0zosb9ypNitlWkrT77JgdomVfSIwkyQP30b0L++LCAveH8l/mGnBJOS58WRzTbuHu24dum9k/SPpxWIsAAA0vz8clszim7sywK2K/X9LakdYFAIxDnnHJqTF7mmZ2n6TLJc00sy2SPi3pcjO7UOW3tknSx2rXRAAA8mHMounuN1d5+J4atAUAUAj5PpknC2YEAgDEy/EQaxYUTQBALCZsBwAA9DQBAPEYngUAIK1iDs9SNAEA8ehpAgCQUkGLJicCAQCQUkP2NCfsagnLatkXN+7eticsSn0dMTn7z4qb3N5Kcdtq0ua4Xa9jU1y7dp8X8znSz4q7oo+91B6W1b4lblsdmZ6/Px8DcZtKHtilaNoR9zer+UADHCtkwnYAANIr6oTtFE0AQDyKJgAAKRV0eJYTgQAASImeJgAgnDE8CwBACjm/kHQWFE0AQDDjmCYAAOMdPU0AQDyGZwEASImiCQBAShRNAABSYO7ZfJlwysGwrCMnnBCW1bq3KSzrxE1BH9MsbqLoeVdsCsta19YVljUwsS0sK+rLZb4xbubw/hPjPrIfmVEKyzqhJ25/L7XG5PR3xG2rUvtgWFbX3J1hWa/vmRyWhaPXkEUTAJBvTG4AAEBaBS2afE8TAICU6GkCAMIVdXiWniYAACnR0wQAxOMrJwAApMBVTgAAOAoFLZoc0wQAICV6mgCAcEU9e5aiCQCIR9EEACAlimZ+tKyNmwy7fU9YlA7Mics6+AeHQ3IGD8X9F29eeVpYVvPUuN8oP+tAWNbgvpiZw0vNgX8xWuKymnfH7Q+luGsBqG960OTo7QMxOZJsQtx23/9wZ1jWiYfComrGvLjDs5wIBABoOGZ2jZm9bGbdZnZHlefbzOyB5PmnzOy0iNelaAIA4rllW0ZhZk2S7pZ0raQFkm42swXDVvuopN3ufqakr0v6UsTbomgCAOJ5xmV0l0rqdveN7t4n6X5Ji4ats0jSvcntH0i60swyT1NE0QQAhBs6rnmsi6SZZraqYllSET9b0qsV97ckj6naOu4+IGmvpBlZ31dDnggEAMi57CcC7XD3hQEtCUVPEwDQaLZKOrXi/pzksarrmFmzpCmSdmZ9YYomACBWxqHZFF9XeVrSfDObZ2atkm6StHzYOsslLU5u3yDpZ+6euf/L8CwAIF4Nv6fp7gNmdruklZKaJC1z9xfM7LOSVrn7ckn3SPpHM+uWtEvlwpoZRRMAEK/Gkxu4+wpJK4Y9dlfF7cOSPhj9ugzPAgCQEj1NAEA4ptEDAGCca8yeZuAnmL0X9Ydlndy1JyyrtTlm4umBwaaQHEl6fddJYVkeuOeVDsTNHN4+42BIzsHeE0JyJKnjubawrP1zS2FZb71sc1jWzkMxF2F4bUPm767/zoRDcX2K3rf1hWW1Twucsf1bcVFvUtCeZmMWTQBAfnGVEwAAQE8TABCvoD1NiiYAIB5FEwCAsZmKe0yTogkAiFfQosmJQAAApERPEwAQazx/5cTMlplZj5mtrXhsupk9Ymbrk3+n1baZAICG4hmXnEozPPsdSdcMe+wOSY+5+3xJjyX3AQAoG69F092fUPlaZJUWSbo3uX2vpPfFNgsAgPw51mOane6+Lbn9mqTOkVY0syWSlkhS04ypx/hyAIBGUtRjmplPBHJ3Nxt587j7UklLJalt3pyQzTh4wf6ImLIDrWFRu3snhWX174yZ8HtCn4XkSFLb6b1hWbed84uwrO9vvSQsq2/ZrJCc/vPjTkz/0K0/Dcu6f1Pcttr39beEZe25IObCAis/+pWQHEna53F/G265++NhWc2H4tpVUwUtmsf6m73dzLokKfm3J65JAICGlvV4Zo4L7rEWzeWSFie3F0t6KKY5AIAiMM+25FWar5zcJ+n/SnqrmW0xs49K+qKk95jZeklXJfcBACi0MY9puvvNIzx1ZXBbAABFkePeYhbMCAQACJfnIdYsKJoAgHgUTQAAUsj5GbBZcJUTAABSoqcJAAhlyVJEFE0AQLyCDs9SNAEA4Yp69izHNAEASKkhe5oDW9rDsjqfDYvSnvlxEym3nRszOfrAxskhOZI06ScnhmV9dce1YVknn7o7LKt/asznSG+K+5g9GHh0aOdvp4RlTS7FvccJfTE5d225PiZI0gQrhWWd8suDYVkTVv86LOv5sKQqCtrTbMiiCQDIOYomAAAp5HzS9SwomgCAeAUtmpwIBABASvQ0AQDhGJ4FACAtiiYAAOkUtafJMU0AAFKipwkAiFXgS4NRNAEA8SiaAACMzVTcY5oUTQBAPIpmfsx4646wrO0z4iYhb1/bFpY1uKYjJKf5/JiJ3yXpvHdvDsta/b/OC8sa/PeTwrIOXRs0UX5v3L5w7/ffE5Y1MSxJ2nlOXNb0K7aF5KztmRWSI0kDz8VNbm//Le73cNaUuP1dV8VFjRcNWTQBAPlmXsyuJkUTABCLs2cBAEiPE4EAAEiroEWTGYEAAEiJniYAIBzDswAApEXRBAAgBS9uT5NjmgAApERPEwAQr6A9TYomACBUvSdsN7Ppkh6QdJqkTZJudPfdVdYblLQmufuKu18/VjbDswCAeO7ZlmzukPSYu8+X9Fhyv5pD7n5hsoxZMKUG7Wn2PRw3YXFHS1iUes8cDMv6o4teCsnp3jszJEeSur8cN0P3ST0Hw7JeuXpSWNZFXb8NyVnbFDdx+OSNcTvpkSkWljXvpu6wrJ6Dk0NyJv5z3AUYWveVwrJ2lWIuwCBJW9vismqpzicCLZJ0eXL7Xkk/l/TJiGB6mgCAPJppZqsqliVH8bOd7j506ZzXJHWOsN4JSfaTZva+NMEN2dMEAORYzITtO9x94UhPmtmjkqoN6XzqDU1xd7MR+71z3X2rmZ0u6WdmtsbdN4zWKIomACCcxY1uV+XuI14N1My2m1mXu28zsy5JPSNkbE3+3WhmP5d0kaRRiybDswCAeJ5xyWa5pMXJ7cWSHhq+gplNM7O25PZMSe+U9OJYwRRNAEA482xLRl+U9B4zWy/pquS+zGyhmX07WeccSavM7DlJj0v6oruPWTQZngUAFIq775R0ZZXHV0n6i+T2rySdf7TZFE0AQCxXxHctc4miCQAIV9QJ2ymaAIB4BS2anAgEAEBK9DQBAKHqPWF7LVE0AQCxYiZdz6WGLJq9lx0KyyoN5HOE+qktc0NyJrb1heRI0v4PHgnLss1xk6wPToybemT1tlNCcvoPx/1qHb4gLEoWd00BrfnlmWFZMy+oOmHLUZt9y8aQHEl6ftPssCzfHzfpfuvMuL9/tURPEwCAtApaNPPZzQIAIIfoaQIAwjE8CwBAGi6pVMyqSdEEAMQrZs2kaAIA4hV1eJYTgQAASImeJgAgHpMbvJmZbZLUK2lQ0oC7L4xoFACgsRV1eDaip/nH7r4jIAcAUASuwp4IxDFNAABSylo0XdJPzewZM1sS0SAAQGMrX+XEMy15lXV49l3uvtXMTpb0iJm95O5PVK6QFNMlktQ0Y2rGlytrao6boHvGIyeEZfW3W1jWvjNaQ3L6eyaH5EiSdwbuyHMPhkWdd8prYVnb75kXktM3OW5fePjOr4RlXfF03Gfb2R94ISxr0xfeEZLzq1t+GJIjSd85+eSwrPtuvTYsq3lD3P6+Piypirg/07mSqafp7luTf3skPSjp0irrLHX3he6+sKmjPcvLAQAaRFF7msdcNM2s3cw6hm5LulrS2qiGAQAalAcsOZVleLZT0oNmNpTzT+7+cEirAADIoWMumu6+UVLg5XEBAMXgTG4AAEBaTG4AAEBa9DQBAEjBJeMrJwAAjG/0NAEA8RieBQAgpWLWTIomACBenmf1yYJjmgAApNSQPc3SK5PCsg52xk2s3bovLEql6X0hOVMW7AnJkaSDm2aEZZ192+awrO6/Oics67HPx0yO/vaHPxGSI0kfvuGvwrImnhc3gf+rPzg/LOtv3/ZPITkXfz5uW530zIGwrO13HgnLOvek/rAsvT0u6k0K2tNsyKIJAMgxV2GvckLRBACEMuX7SiVZUDQBAPEKWjQ5EQgAgJToaQIA4hW0p0nRBADE4kQgAADS40QgAADSKmjR5EQgAABSoqcJAAjmhe1pUjQBALFcFE0AAFLj7Nn8mPZiXNbuc+M+DfUvjJvg2fe2heQceLgzJEeSTt04EJb1+vvjJlk/cs6hsKy7tl0dkjP51y0hOZLU1L0xLGvgD94alnX9GWvCsj75yJ+F5Jzz8LaQHEk6fHrcBQoO98X9nelojpv8HUevIYsmACDf+MoJAABpUTQBAEjBJZUomgAApFDcr5wwuQEAACnR0wQAxCtoT5OiCQCIV9CiyfAsACDW0IlAWZYMzOyDZvaCmZXMbOEo611jZi+bWbeZ3ZEmm6IJAAjmkpeyLdmslfQBSU+MtIKZNUm6W9K1khZIutnMFowVzPAsAKBQ3H2dJJnZaKtdKqnb3Tcm694vaZGkUeeco2gCAOJlP6Y508xWVdxf6u5Ls4ZWmC3p1Yr7WyT94Vg/RNEEAMSKmdxgh7uPdjzyUUmzqjz1KXd/KOuLj6Qhi+buuLm+NTizPyzr0tmvjr1SSr/ac2ZIzrSX4t6ft4w61HFU5izeEJY1czBuN179rbeF5DRPComRJP36G28Jy+ro2BWW9cN1F4VlnfJ4zL614SNdITmSZGf3hmXppY6wqH9dfWFYlvTdwKxhanz2rLtflTFiq6RTK+7PSR4bFScCAQDGo6clzTezeWbWKukmScvH+iGKJgAgnnu2JQMze7+ZbZH0dkn/YmYrk8dPMbMV5eb5gKTbJa2UtE7S9939hbGyG3J4FgCQZ/Wde9bdH5T0YJXHfyvpuor7KyStOJpsiiYAIJZLKmX+rmUuUTQBAPGYRg8AgPGNniYAIF5Be5oUTQBAsOyTrucVRRMAEMslzz7pei5xTBMAgJToaQIA4jE8CwBASpwIlB9nXbo5LOvsE7eHZf3LhnPDsia+2hKSs/k/xh1XOOf8V8Kytu0/MSzrwC9OCsuyaTE5++cNxgRJat0wMSxrz7TWsCzFzd+vvfNijhS1BM6xvn9X3Ha//Kq1YVl7+uLa9dLfhkW9kTuTGwAAkFpBe5qcCAQAQEr0NAEA4ZzhWQAA0qjvVU5qiaIJAIjl4isnAACkxoxAAACMb/Q0AQChXJIXdHg2U0/TzK4xs5fNrNvM7ohqFACggbmXh2ezLDl1zD1NM2uSdLek90jaIulpM1vu7i9GNQ4A0Jjoab7ZpZK63X2ju/dJul/SophmAQCQP1mOac6W9GrF/S2S/jBbcwAAhZDjIdYsan4ikJktkbQkuXtk80fuyDxzcdx07dLD6VedKWlH4EsfTyFtj9zuR2Hcb/c6oe2Bvpt+1ePd9rm1CO3V7pWP+g9mZozJ1f/hkCxFc6ukUyvuz0keewN3XyppqSSZ2Sp3X5jhNeuGttcHba8P2l4fjdz2Su5+Tb3bUCtZjmk+LWm+mc0zs1ZJN0laHtMsAADy55h7mu4+YGa3S1opqUnSMnd/IaxlAADkTKZjmu6+QtKKo/iRpVler85oe33Q9vqg7fXRyG0fF8wLOhM9AADRmHsWAICUalI0x5pez8zazOyB5PmnzOy0WrTjaJnZqWb2uJm9aGYvmNnHq6xzuZntNbPVyXJXPdpajZltMrM1SbtWVXnezOzvk+3+vJldXI92Dmdmb63YnqvNbJ+ZfWLYOrnZ7ma2zMx6zGxtxWPTzewRM1uf/DtthJ9dnKyz3swWH79W/+71q7X9K2b2UrJPPGhmU0f42VH3r1oboe2fMbOtFfvFdSP8bF2n/Byh7Q9UtHuTma0e4Wfrut0xjLuHLiqfFLRB0umSWiU9J2nBsHX+k6RvJbdvkvRAdDuOse1dki5ObndI+nWVtl8u6cf1busI7d8kaeYoz18n6SeSTNJlkp6qd5tH2H9ekzQ3r9td0rslXSxpbcVjX5Z0R3L7DklfqvJz0yVtTP6dltyeloO2Xy2pObn9pWptT7N/1antn5H031PsU6P+TapH24c9/3eS7srjdmd541KLnmaa6fUWSbo3uf0DSVeamdWgLUfF3be5+7PJ7V5J61Se+agoFkn6rpc9KWmqmXXVu1HDXClpg7vXaS6Fsbn7E5J2DXu4cp++V9L7qvzoeyU94u673H23pEckHdfvs1Vru7v/1N0HkrtPqvyd69wZYbunUfcpP0dre/K370ZJ9x3PNuHY1KJoVpteb3jh+d06yS/rXkkzatCWY5YMGV8k6akqT7/dzJ4zs5+Y2bnHt2Wjckk/NbNnkpmYhkvzf1NvN2nkPx553e6S1Onu25Lbr0nqrLJOI2z/W1UejahmrP2rXm5PhpaXjTAsnvft/keStrv7+hGez+t2H5c4EagKM5ss6YeSPuHu+4Y9/azKQ4cXSPqGpP9znJs3mne5+8WSrpV0m5m9u94NOhrJJBnXS/rfVZ7O83Z/A3d3lf/QNRQz+5SkAUnfG2GVPO5f35R0hqQLJW1TeZiz0dys0XuZedzu41Ytimaa6fV+t46ZNUuaImlnDdpy1MysReWC+T13/9Hw5919n7vvT26vkNRiZlnnWAzh7luTf3skPajysFSlVFMf1tG1kp519+3Dn8jzdk9sHxrqTv7tqbJObre/md0i6U8k/XlS9N8kxf513Ln7dncfdPeSpH8YoU153u7Nkj4g6YGR1snjdh/PalE000yvt1zS0JmDN0j62Ui/qMdTcmzhHknr3P1rI6wza+j4q5ldqvI2rHvBN7N2M+sYuq3yyR3DJ8dfLunDyVm0l0naWzGkmAcjfuLO63avULlPL5b0UJV1Vkq62symJcOIVyeP1ZWZXSPpf0i63t0PjrBOmv3ruBt2TP79qt6mPE/5eZWkl9x9S7Un87rdx7VanF2k8lmav1b5jLVPJY99VuVfSkk6QeUhuG5J/y7p9HqfEZW0610qD6s9L2l1slwn6S8l/WWyzu2SXlD5DLwnJb2j3u1O2nV60qbnkvYNbffKtpvKFw7fIGmNpIX1bndF+9tVLoJTKh7L5XZXubBvk9Sv8vGxj6p8TP4xSeslPSpperLuQknfrvjZW5P9vlvSR3LS9m6Vj/kN7fNDZ7afImnFaPtXDtr+j8m+/LzKhbBreNuT+2/6m1TvtiePf2doH69YN1fbneWNCzMCAQCQEicCAQCQEkUTAICUKJoAAKRE0QQAICWKJgAAKVE0AQBIiaIJAEBKFE0AAFL6/zi7Cx8qPrZRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dmdc.dmd_time['dt'] = .5\n",
    "new_u = np.random.rand(s['u'].shape[0], dmdc.dynamics.shape[1]-1)\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.pcolor(dmdc.reconstructed_data(new_u).real)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recovering unstable dynamics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this subsection of the `DMDc` tutorial we are going to investigate the stabilizing effect of control in an unstable linear system. For simplicity we consider a 3D system with a known `B` operator. This tutorial is inspired by `Example 1` in the paper [Dynamic mode decomposition with control](https://arxiv.org/abs/1409.6358)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_instants = 5\n",
    "\n",
    "A = np.array([[2.1, 0, 0], [0, 0.5, 0], [0, 0, 6.9]])\n",
    "B = np.array([1.3 , 0, 3])[:, None] * 2\n",
    "\n",
    "snapshots = np.zeros((3, time_instants))\n",
    "snapshots[:, 0] = np.array([5, 6, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We construct our snapshots matrix using the usual formula: $\\mathbf{x}_{k+1} = \\mathbf{A}\\mathbf{x}_k + \\mathbf{B}\\mathbf{u}_k$, where `A`,`B` are specified above, and $\\mathbf{u}_k = -\\mathbf{x}_k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5.     -2.5     1.25   -0.625   0.3125]\n",
      " [ 6.      3.      1.5     0.75    0.375 ]\n",
      " [ 3.      2.7     2.43    2.187   1.9683]]\n"
     ]
    }
   ],
   "source": [
    "for i in range(1, time_instants):\n",
    "    previous_snapshot = snapshots[:, i-1][:, None]\n",
    "    snapshots[:, i] = (A.dot(previous_snapshot) - B*previous_snapshot).flatten()\n",
    "print(snapshots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We construct the control input matrix used in the formula, in order to feed it to the `DMDc` constructor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-5.     2.5   -1.25   0.625]\n",
      " [-6.    -3.    -1.5   -0.75 ]\n",
      " [-3.    -2.7   -2.43  -2.187]]\n"
     ]
    }
   ],
   "source": [
    "I = -snapshots[:, :-1]\n",
    "print(I)\n",
    "\n",
    "dmdc = DMDc(svd_rank=-1).fit(snapshots, I, np.diag(B.flatten()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now pause `DMDc` for one moment, and analyze the system resulting from the set of snapshots, as if we did not know that it is a control system. Therefore we use the standard DMD algorithm, and then plot the eigenvalues to assess the stability of its dynamics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFBCAYAAAA7XhdpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABjFElEQVR4nO2deVxU1/n/388MMIDggguioKACRkwwSFSMcYuNSdqoScy+mNY0TbN8v03zTZO0v2Zr0iZN0rTNZpYmMfu+mK3Z1bhGQFHRgEY0oogLLiAwwMz5/TEDmWGRAWbmXpzzfr3ui5l7zj3nM8/MPJw595znEaUUGo1Gowk8FqMFaDQaTaigHa5Go9EECe1wNRqNJkhoh6vRaDRBQjtcjUajCRLa4Wo0Gk2QMNThishzIrJXRDa2US4i8m8R2Soi60Uky6NsnohscR/zgqdao9FoOofRI9wXgDOPUX4WkOo+rgGeBBCROOBOYDwwDrhTRPoEVKlGo9F0EUMdrlJqKVBxjCqzgReVi1VAbxFJAGYCXyilKpRSB4EvOLbj1mg0GsMxeoTbHoOBnR7PS93n2jqv0Wg0piXMaAGBRkSuwTUdQVRU1NikpCSv8vDwcGw2G7W1tURFRVFZWdmijZiYGOx2OxEREdTX11NXV+dVHhERQXh4OHV1ddhsNqqqqo7ZRk1NDU6ns8NtxMbGUlNTQ2RkJHa7nfr6eq9ym82G1WqloaGB8PBwjh49esw2qqurW+horw0RISYmhurqaqKjo6mpqaGhocGrTmRkJCKC0+nEarVSXV3tVW6xWOjRo0dTG1VVVTTfYt5eG2IRwmxhNNgbiOkRQ8WRCnCCEyeCEGGNoNJZSaW9kjBLGCJC/4j+OJWTBqdLrxJFVHQUUi8cchwi0hnpahshzBJGg7OBaqpBgVM5SYhNoPJoJXWOOkQEQbCF24iOjuZQ5SFUuCLcEY4FC4I0aY2OjsbhcGCxWFBKUVtb6/VawsLCiIqKarJHa5/BjrbRmk179OhBfX09YWFhOBwO7Ha7V3nz78KRI0cQkS614Y/vU2VlZQsd/vhOdvT7tHbt2v1Kqf4tGuoAZne4uwBPD5noPrcLmNrs/OLWGlBKPQ08DZCVlaXy8/MDobND1NTUEBUVZbSMbqPjiP0Ia8vWkt4vnXpHPTNemsHOIzsZHDuYG8fdyI3jb+Rv3/6NMEsYg2IHkdInhYlJE6lz1BFuCW/xZe2sjkbqHHWUV5VTVlVGet90ahtqeey7xyirKqOsqozZ6bO5Zuw1ZDyRwcCYgYxNGMvU5KmcnXq2X3UEGq3DGxHZ0eVGlFKGHkAysLGNsp8DnwICTAC+c5+PA0qAPu6jBIhrr6/Ro0crM1BcXGy0BKWUOXUcrj2svin5RtU76tXX275WaY+mqej7olXOszlqxY8rVF1DndpQvkHVO+r9quHll19WeXl5fm2zvKpcfVL8ifrLkr+oB5Y9oJRS6op3r1DTF05Xt3x+i3p9w+vK6XS2uM6M74uRmEUHkKu66O9EGRgtTERewzVS7QeU41p5EA6glFogrqHJY7huiFUDv1RK5bqv/RXwR3dT9ymlnm+vv5NPPlmtXbvW3y+jwxw5coSePXsaLcMUOsoqy+ihevDW1rf4+4q/s+vILk6KP4l3LnyHcGs4ZZVlnND/BMIsgfsxtn//fuLj45k3bx7PPfdcwPoB2Hd0H3lleeTuzmX7oe08O+tZHln5CB8Wf8jYhLFMT5nOuH7j6Nunb0B1+IIZPh9m0iEieUqp7C61YaTDDTYnnXSSWr9+vdEy2L59O8nJyUbLMExH7u5c3ix8k0VFi2hwNvD5OZ/T0LMBe4M94M61NW6//XaefPJJlFLs3Lkz6F/uA9UHWLN7DWt2reGLbV/w+MTHcfZysnbPWn6e+nP69+jStGGnCfXPaXO0w+0g6enpqqioyGgZLF68mKlTpxotI2g6Ku2VfP7D5xTuK+SOKXfwz1X/5GDNQWalzyIrIYslS5YYZo/9+/eTnp5OcnIyFouFc889lz/+8Y/tXxhAli5dSs8+PVEOhUIhCBaL9424YKCU8nn++3jSYbVaGTZsWIt5Y384XLPfNNN0U47WHaVHRA9u+u9N/Gftf8hJymFO+hyUUvxuwu+MltfEww8/zIUXXkhBQQHTp0/nkUce4cYbbyQ2NtYwTbGxsSQMSKB///6ICA7lwCpWHMpBnaOOMAnDarFiEUtAHZHD4cBqtQasfTPqcDqd7Nu3j23btpGRkeH39rXD1fiN8qpynlv7HG9teouo8CiW/2o514y9hrun3U1Pm/FzcK0xatQoZsyYwfnnn8/AgQP561//2mKZW7CxWCz0798fi8W1TD5MXF9TK1ZsVluT442wRmDBglM5A+58Q4VG25eXlwekfe1wNV1CKcXOIzsZ0msIN3x6A71tvXlk5iOcOuRUAE7of4LBCo/NFVdc4fX817/+tUFKvGl0tp6ICFaxYsUK7gGfw+lyvgBhlrCmNceaztOa7f1FSDnc8PBwoyUAMHDgQKMlAF3TUWmv5JUNr/DEmieIiYhh+a+W8+bcNzv1ZTeLPfr2NX5lQEexWqxESiRO5aTeWY8Va9MSJKul8z/DrVYrqampTc/PO+88/vrXv3LRRRdx6623kpWVdYyr/Utrn6no6OgWm2G6AyHlcG02m9ESAEhJSTFaAtA5HWWVZSTEJnDLF7ewr3ofj8x8hOkp07s0qjKLPRISEoyWALTuYNqrbxVrk4N1OB3YHXbEIYRbwrFarB1u02azsXnz5hbXvfHGGx1qxx8cTyN2s8dS8CvNt0MaRXFxsdESAN911DnqeLPwTaa+MJWJz7l2cD358yd558J3OH3Y6V3+QpjFHjt37my/UhDo6Mqhp56CQYPAYnH9ffYZK1FhUYRbw2lQDSgUTuXEqZztN+ZB823fAOPGjePbb78F4J///CfJycmcdNJJXHzxxcyb54qSunv3bmbOnMno0aMZPXo0X3zxBQA333wzF154IePGjSMxMZH77rsPgOuvv57777+/qY+bb76ZO++8k8OHD5OTk0NGRgZpaWm8+uqrLfR88sknTJs2ren5vHnzePTRRwFYtmwZp5xyChkZGZx22mns2OHaKHbfffcxfPhw0tLS+MUvftEhm3SZru6c6E7H2LFjfdtSEmAcDofREpRS7euoqa9RSil1x9d3qCnPT1FvbnxT1TXUBV1HMMjJyVHffvut0TKUUkoVFBT4XHfBAqUiI5WCn47ISNd5T+ob6tXRuqOqpr5GNTga2m3XYrGo9PT0puPZZ59VSil1yimnqKVLl6qSkhI1aNAgVV5ermpra9XYsWPVlVdeqZRS6pxzzlGfffaZUsq1SywlJUUppdTvf/97NWbMGFVdXa12796tevXqpWpra9Xy5ctVdnZ2U9/Dhg1TW7duVXV1daqiokI5nU61e/dulZSU1PRZiYqKUkop9fHHH6upU6c2XXvllVeqf//736q2tlaNGTNG7dq1Syml1LPPPqsuuOACpZRS/fv3V9XV1Uoppfbt29fq62/tPcAPO81CakqhtUAaRrB06VJTrMNtS0dtQy2Pf/c4D654kLxr8rhjyh1dmg/srI5gs27dOiZNmmS0jA6NcO++G5r/cKutdZ3/zW9+OhdmdS0la3A2UOeoI1LcgXra+HVis9koLCxscznWsmXLmDBhAgMGDADg3HPPbfqlsmzZMq9fLVVVVRw+fBiAmTNnEhUVRVRUFHFxcezatYuJEydy4MABtm/fzp49e+jVqxfDhw/Hbrfzu9/9jpUrV2KxWNi7dy+7du2ieQCq1tiwYQNbtmxh+vTpgGu03qg1PT2dc889l9mzZ3PppZe225Y/CSmHq2mftWVrmf36bLISsvjqyq8Y3FNHvTQze/b4fl5ECLeGN61kqG2oRXCds4j/ZheVUuTn5xMdHd2izPM+itVqbYrQNWvWLF555RXKyso4//zzAXj66afZv38/BQUFREVFMXjw4BY3ysLCwrymPhqnDZ1OJyNGjGDdunUtNHzzzTf897//5YMPPuDBBx+kqKgoaDfUQ2oOV9M6Sine//591uxaQ2rfVN6Y+wbvX/w+GQP8v/Bb41/aWuBxrIUfjaNam9XW5HjrHfVtX9AKp556KqtWrWLfvn3U19fz/vvvN5WddtppXnOyK1eubLe9K664grfffpuPPvqIyy+/HIDDhw/Tv39/IiIi+Oijj9i9e3eL64YPH87WrVupqalh//79LFu2DICTTjqJiooKvvrqKwDsdjt5eXk4HA5++OEHfvGLX/DYY495jb6DgR7hhjhLti/htq9uo7q+mifOdi3xyknKMVqWxkfuvBN+9zvvaYXISNf59hBxxQ0Os4Q1zTE2OBsIs4Rht9u9dlqdfvrpPP74403PU1JS+P3vf092dja9evUiNTWVXr16AfDUU09x9dVXk5aWhsPhYMKECeTkHPszNXbsWI4ePUp8fDxDhw4FYP78+Zx11lmMHDmSzMzMVlezDB8+nFmzZjFy5EiSkpKaNEdGRvLGG29w4403UllZicPh4LrrrmP06NFceumlVFZWopTi6quvpl+/fu0by190dRK4Ox1paWmtTpAHm2+++cZoCaqmvkZ9/c3XavrC6erlgpeVw2ncjSsz2CMnJ0c9+uijRstQSim1bt26DtVfsECphASlRFx/m98w8xWn06lq62vV0bqjqq6hTjU0HPvm2qFDh5RSStXV1alp06apF1988Zj1O0t7OgKBvmnmB2JiYoyWANDuf/tAUnKwhD99/Sf2V+/nwws/5KupXxmmpREj7eHJ6NGjjZYAdHzd6W9+432DrCv92sJsOJwO6p31iEWabuC1pukPf/gDS5YswW63M3XqVC677LKui2iFQO78CjbHzyvxgebpQIyipKTEkH5fXv8y454dx8h+I3nnwncM09Ecs+goKyszWgLQ8XW4/sZqsRIZFokFCw3OBuwOe6treJ966im+//57SkpKeP755wPmGI22hz8JqRFuRESE0RIASExMDGp/JQdLGNxzMKlxqSz75TLS+6UboqMtzKKjf39j4s42xyw7q0SEMAlDORW1DbWEW8INidVgFnv4g5Aa4TZPEmcUgYpE1ByncvLo6kc55ZlTyN2dy/jE8U3ONpg62sMsOg4ePGi0BMA8IzqlVNONNZvVhkI1nQ+2juOFkHK4zTN7GkUwtpDWNtQyfeF0Xtv4Gst/tZyJSRMN0eELZtFhFsdvFjwdndViJcLq+oVY01BDvaM+aI7weHK4ITWlEAo4lZOCPQWcnHAyt556K2cMPyOgu8Q0oYWIYLPaqHPU4VAOV0xeP26aON7RljqO2HZwG9MXTufmz2/GqZyclXqWdraaDlNUVOQVmhF+CigDHjfV3I62ccnTt99+yy9/+UvAFVTmyy+/7FC/U6ZMYf/+/T7Xnzt3Li+88EKH+jAaPcI9Tvj8h8+57N3LuO3U2/jdhN/pUcdxxiXvXMLrG19vel7yvyWsLF3Jpe/8FAvgoZ89xM0Tb0bu/ukm07TkaXw972umL5zON9u/aTqv7uzaz/TGuV2AuoY6nDg5ddKpnHbaaQB89dVXxMTEMGPGDJ/bXLJkSYtzTqfTNKl+/EFIOVyzrFIYPNh/8Qkqaiqoqa8hKyHLawVCsHV0BbPoMMsqheZkPJFB4XWFvHb+a17nk3snc8noS1rUb82Zfj3v6w7329bqgHHjxpGdnc2yZcs4cuQITz31FNNmTGPRR4v41z/+xYIFC1i4cCEWi4U33niDf/3rX8ycObPp+sOHDzN//nwaM2j/6U9/Yt68eQwePJjc3FyOHDnCmWeeycknn8z69ev5+OOPefHFF3nrrbcQEWbMmOG16w1cAXNuuukmqquriYuL4+WXX27asWYmQmoYZJaMD/5aBlW4t5Bxz4zj3c3v0i+6X4ecrT91dBWz6GiMJmU0zR3dpn2b/Nr+wyse7pQOTxoaGli/fj0PPvggd999d9MWYXBF47ryyiv57W9/y/fff+/lbMGVlr5Xr14UFxdTXFzM2Wef3aL9HTt2cMMNN7B161bWr1/PJ598Ql5eHkVFRU1TG43Y7XZuvPFGPvjgAwoLC7nqqqu45ZZbfHqNwSakHK5ZVin4Y6H/x8UfM23hNO6Ycgc3jr/RMB3+wCw6WguOYgSBviv/f1/83zHLGx1tcx2eDviCCy4AXLsES0tLAZqmsRqcDU1LyFpjyZIl3HTTTU3PW/tlkZCQ0BRa8fPPP+eKK65oyqTc/B+jZyjGkSNH8sADD5jmvWxOSE0pmCXFTlpaWqevbfwSKBSLLlnEhMQJhujwJ2bR4Uuc1WDQfGTZPzq4Ux3x8fEcPnzYS0dFRYVX8JjISFc8XavVisPh8LreKu5UP8rRlFG4o3iGdmxv48OxQjGaDUNHuCJypogUichWEbmtlfJHRGSd+ygWkUMeZQ6PskW+9FdVVeVH9Z2nMYRcR6mpr+Hy9y7nxYIX+UXaL7rkbLuiw9+YRUfjnKLRNB9Z7r1lb1D779WrFwMGDODDDz909b93L998803TiLM9RIRePXtRVVnVlFHYk6lTp/LII480Pd+3b98x2/vZz37GSy+91JRAYO9eb3u0FYrRjBjmcEXECjwOnAWMAi4RkVGedZRSNymlxiilxgCPAu96FNc0limlZgVLt1GUHinltOddd4AvzLjQYDWaYHLzZzf7tb1Xz2+ZG6w5L774Ivfddx8jR45k8uTJ3HbbbYwaNard6xo5//zz+eTjTxgzegz//ey/NDgamsr++te/cujQIVJTU0lPT+fTTz89ZlvnnXceZ511FmPGjGHkyJHcc889XuWNoRhvvfVW0tPTycjIaHXFgynoarixzh5ADvCZx/PbgduPUX8F8DOP51Ud7bM7h2f89aJfqweWPaCcTqehOgKBGXSYOTwjd+F1fLv9W/Xt9m+9zv3+v79XSinV/+/9m86NenyUUkqpi9++2KtuycESn3T4Kyyiw+lQR+uOKnu9vVOfXx2e0T8MBjz3dJYC41urKCJDgRTAc21LpIjkAg3A/Uqp9wOk01BeKniJU4ecyoJfLNBra0OUttbMtna+temH185/rcWSsmBiEQtRYVHYHXbsDntTpolQpLvcNLsYeFsp5Tk7P1QptUtEhgFfi8gGpdQPzS8UkWuAa8B1d3Px4sVe5QMHDiQlJYXi4mIyMjJYunRpi85zcnIoKSkhMTGR8vLyFnv/Bw8eTGJiIiUlJaSlpbU6Jzlp0iSKi4tJSUnBbre30JGUlER8fDylpaWkpKSwcuVKSo+Ucrj2MMX7ivkx7EcmT55MYWEhaWlplJSUsKdZ4qrk5GTi4uIoLy8nMTGR1atXN7cFU6ZMoaCggIyMDGpra1voGDZsGLGxsVRUVBAfH09ubq5XeVhYGJMmTSI/P58xY8ZQWFjIgQMHvOqkpqZis9moqqoiLi6O/Px8r3KbzUZOTg55eXmMHTuWmpqaFjrS09OxWq3U1tYSGxtLQUGBV3lUVBTjx49vaiM/P58jR4541Rk1ahQOhwOn04nNZmPjxo1e5TExMWRnZzfN97Vmj9GjR2O327FYLFitVjZt8l6i1bNnT7Kyspp0rF69mpqaGq86mZmZVFZWEhkZicPhoKioyKu8T58+ZGZmNrUBtLgRZbFYvOZ2PR971nE6na3eyPJsQ0Q8fzk2ISKISFMbSqk2dfjSRqOeRsIlvGn1gsPhaGqrrTY8++msjsY6raV8b9TXmg6lFNu3b/f6PvmFrg6RO3vQgSkFYC0w8RhtvQDMba/PrKysNn5ABJf6+vp261z74bVq8vOT1YHqA4bqCAZm0JGTk6OWLFlitAylVMfSpAcSf05feeJwuKYYfEnXHkgdxyJQUwpG/kZdA6SKSIqIROAaxbZYbSAiI4E+wEqPc31ExOZ+3A84FWh3dbhZApB7ppBuToOzAaUU87Pm8+llnxIXFWeIjmBiFh1miVqmlGp1RGaEjkBgsViIsEZgd9hxOFuOxIOloy0CaXvDHK5SqgG4AfgM2Ay8qZQqFJF7RMRz1cHFwOvK2+onALkiUgB8g2sOt12Ha5atva0lwwOwN9iZ++ZcFhYsJHtQNtHhLdNMB0NHsDGLjkGDBhktAXCtbd23b5/hTjeQ86xhljAirBE0OBvarRvM+V6n08m+ffsCFrvB0DlcpdQnwCfNzt3R7PldrVy3Ajixo/2ZJQB5aWlpi2hMNfU1nPfmeUSHR3PpiZe2cWXgdRiBWXQ0X99pFFarlf379xsen1e55z8D3g8KVNuONVg6GrFarQwbNiwgbXeXm2Z+wSxbe3ft2tXCwfxj5T/oE9mHF899sWlPuhE6jMAsOtpbgB8sysvLmTp1qtEyWLx4cVB0LN2xlLlvzuX1ua8zPaXl5opg6QgGIeVwzYi9wc6eqj384dQ/YBGLjl+rCTkmD53M2xe+zdw35/LuRe8yacgkoyUFDL2w00DqHHVc+PaF3L/sfsKt4drZakKWyUMn88p5r1BprzRaSkDRI1yDaHA2cNm7l6GU4l9n/ctoORqN4fxs+M8A1/Ta1OSpZCVkGazI/4TUCNcsqxSSkpIo2u9a/P7WBW81Rc43QocZMIuO+Ph4oyUA5rGHUTqSeydz9itns3HvRkN1BIKQGuGaJQD55prNnD3sbN664C1DdZjFwZhFR58+fYyWAJjHHkbpOO+E87A32Jn58kxWzl9pGnv4g5Aa4ZphlcKC3AV8lv8Z1fXVRktpChxtNGbRYZZVCmaxh5E6LjnxEp455xn6RfczjT38QUg5XKMDkC/ZvoQ7F9/J9WdcT0xEjKFawDwbDsyiIyEhwWgJgHnsYbSOs1PPZn/1fhbuWBj03WaBIqQcrtEByN/d/C6vnPcKe77f037lILBy5cr2KwUBs+hoHtzGKMxiDzPoGNBjAEMrh3Lv0nuNluIXQsrhGkVVXRVF+4v411n/YsYw39NGazShTmRYJMPjhvN0/tO8t/k9o+V0Ge1wA4xTObnyvSt57LvHjJai0XRLwi3hvHfRewEN5BQstMMNMHcvvpvyo+U8dMZDRkvRaLot2YOymTx0Mg8se4D91fuNltNptMMNIBU1FXy05SPevfBdbGHmyBis0XRXRIRDtYeY++Zc6h3mCETVUeR4ufvnC9nZ2ap5BoNAUXqklEGxgxCkRaSjxijzRqN1/MTEiRP5+9//zqRJxu/jN4M9zKrD4XQw5405JPVM4omfPxFUHSKSp5TK7kobxlsziDRPfRIo9h7dy6TnJvHtjm9bDStXWFgYFB3toXV4U1JSYrQEwDz2MKMOq8XKK+e9Qt+ovj4FLzcbIeVwIyMjA95HvaOeuW/O5fKTLmdK8pRW66SlpQVchy9oHd6YZQupWexhVh09bT35y/S/ULivkBU7VxikqnOElMMNRoqd979/n8iwSO6Zdk+bdcwyktI6vCkrKzNaAmAee5hdR+mRUi5959JuFWEspBxuoDM+KKWYO2ou71/8/jFTmjfPtmsUWoc3zbMPG4VZ7GF2HWenns3pKadzyxe3BFlR5wkphxtI6h31TH9xOkUHigKei0yj0bj4x8x/sGHvBg7XHjZaik9oh+snHlzxIBHWCNL7phstRaMJGXpF9mLZL5fRI6KHKQJCtUdIhWcMFBv3buSRVY+Qd01eUJPdaTQa1/rcB759gJ1HdrLgFwuMlnNMQmqEG6hoYX0i+/DyuS8zpNcQn+onJycHREdH0Tq8MUu0MLPYozvpuGHcDXyy5RO+3PZl4AV1gZByuIHINf9W4VuICDNHzPT5mrg4c+wJ1zq86dmzp9ESAPPYozvp6BXZi2fOeYarF11NTX1w1tt3hpByuA0NDX5tb0P5Bq775LoOL8AuLy/3q47OonV4U1FRYbQEwDz26G46Zo6YybsXvUtUeFSAFXWekHK4/kyxU++o56oPruJvp/+NpF4dWzCfmJjoNx1dQevwZsCAAUZLAMxjj+6oIyshi4XrFvLVtq8CqKjzGOpwReRMESkSka0iclsr5VeJyD4RWec+rvYomyciW9zHPF/6O3r0qN+0ry9fT0rvFOafPL/D165evdpvOrqC1uGNWbaymsUe3VXHwJiBzF8035QbIgxzuCJiBR4HzgJGAZeIyKhWqr6hlBrjPp51XxsH3AmMB8YBd4pI0DIAHrEfISshi7cvfFuvStBoTMbMETOZMWyGKTdEGDnCHQdsVUptU0rVAa8Ds328dibwhVKqQil1EPgCODNAOlsw9825fFT8UbC602g0HeThMx7GqZymC3BjpMMdDOz0eF7qPtec80VkvYi8LSKNk6W+Xut3vtr2FdsObuPMEUHz7xqNpoP0iuzF0+c8zaHaQ0ZL8cLsGx8+BF5TStlF5DfAQmB6RxoQkWuAa8B1U2Tx4sVe5QMHDiQlJYXi4mIyMjJYunRpizZycnIoKSlh8ODBvLPiHR5Me5Dl3y5vKh88eDCJiYmUlJSQlpbGsmXLWrQxadIkiouLSUlJwW63t9CRlJREfHw8paWlpKSktJrAb/LkyRQWFpKWlkZJSUmLPebJycnExcVRXl5OYmJii7kvEWHKlCkUFBSQkZFBbW1tCx3Dhg0jNjaWiooK4uPjaR4/OCwsjEmTJpGfn8+YMWMoLCxsEYMgNTUVm81GVVUVcXFx5Ofne5XbbDZycnLIy8tj7Nix1NTUtNCRnp6O1WqltraW2NhYCgoKvMqjoqIYP358Uxv5+fkcOXLEq86oUaNwOBw4nU5sNluLJJExMTFkZ2eTl5cH0Ko9Ro8ejd1ux2KxYLVa2bRpk1d5z549ycrKatKxevXqFmFAMzMzqaysJDIyEofDQVFRkVd5nz59yMzMbGrj6NGjLXRkZWVRUVFBTEwMdrudLVu2eJX37duXjIwM1q1bR1ZWFsuWLWuxKic7O5vy8nLi4uKorKxk27ZtXuUDBgwgLS2NwsJCMjMzqaqqaqFj/PjxlJaWEh8fT0VFBdu3b/cq78j3KTExkfLycnbu3OlV3vz71JoOz+9TaWkpu3bt8ipPSkpiwIAB3P3+3Zw38Tzv4Zmbjn6f/IJSypADyAE+83h+O3D7MepbgcPux5cAT3mUPQVc0l6fY8eOVV2hwdGgXtvwmnI4HV1qR2M+cnJy1PLly42WofEzLxW8pHKezVFOp7PLbQG5qot+z8gphTVAqoikiEgEcDGwyLOCiHhu/ZkFbHY//gw4Q0T6uG+WneE+d0y6EoC8wdnAp1s/5aKMi44ZCcwXmo/WjELr8Kb5qNEozGKP40HHJaMvoaquig+LP/Sjos5j2JSCUqpBRG7A5SitwHNKqUIRuQfXf5JFwP+IyCygAagArnJfWyEif8HltAHuUUq1u2q9KwHIn1/7PK9tfI2fp/680200kpGR0eU2/IHW4U1KSorREgDz2ON40GG1WHnx3BfpH93fj4o6j6HrcJVSnyil0pRSw5VS97nP3eF2tiilbldKZSilMpVS05RS33tc+5xSaoT7eN6X/jobgLy6vpq7l9zN/TPu98sysOLi4i634Q+0Dm+azyUahVnscbzoGDNwDHWOOhZvX+wfQV0gpHaadTYA+asbXmVC4gTGDR7nFx179+71SztdRevw5uDBg0ZLAMxjj+NJx+7K3cx7fx61DbV+UNR5QsrhdpZfnfwr/jPrP0bL0Gg0neTUIaeSGZ/Jk2ueNFSHdrjt8I+V/2DFzhX0iuxltBSNRtMF/nr6X1mQt8DQzRDa4R6DXUd2cd+395HcO9loKRqNpouMHjCadb9Zh9Xi/zCtvhJSDrejAcjv+/Y+5p88n8Se/o2aNGzYML+211m0Dm8GDRpktATAPPY4HnVEhkVyyTuXsO/oPr+12RFCyuF2JAC5UoqyqjJuzrnZ7zpiY2P93mZn0Dq8iY42R/JPs9jjeNQhIkSHRfOftcbckwkph9vRAOTvXfQe8THxftdhlkDXWoc3lZXmCOdnFnscrzquO+U6FuQaM5cbUg7X1wDkTuVk0vOT+PHwjwHRER/vfyfeGbQOb/r0CVqEz2NiFnscrzrGDhpL9qBsth/a7td2fSGkHK6vAci/LvmaSnslST07lsnBV5oHhDEKrcOb77//vv1KQcAs9jiedbx94dsMjxvu93bbI6Qcrq88seYJrj/leh1cXKM5jjnrlbP4oeKHoPapHW4znMpJ36i+XHbSZUZL0Wg0AWR0/9EsyF0Q1D61w22Gw+ngmVnPEBMRY7QUjUYTQK7NvpYXCl4Ialr1kHK47U0R1DnqSH8snT1Ve45Zr6uEhZkj7rvW4U1Hlg0GErPY43jXMTxuOFecdAW7Kne1X9lPiCuubmiQnZ2tjjUB/2bhmzyZ+yTfzPsmiKo0ZmDixIk89NBDTJw40WgpmiDjVE6fYlyLSJ5SKrsrfYXUCLe6uvqY5U+seYLrsq8LuI7mKWeMQuvwpnnqG6Mwiz1CRce0hdPI250X0D4aCSmHe6ydRA6ng1OTTmXOyDkB1zFmzJiA9+ELWoc3qampRksAzGOPUNFx5vAzeWLNEwHto5GQcrjHSrFTWVfJfaffR7jVt80RXaGwsDDgffiC1uFNSUmJ0RIA89gjVHTMz5rPO5vf4WBN4OMhh5TDbWtrr8PpIP2xdMqryoOio3mWW6PQOrw5fPiw0RIA89gjVHQM6DGAO6bcwcHawDtcc9yGNJhVpasYGDMwIHETNBqN+fl9zu9xKmfA+wmpEW5bLCpaxKy0WUbL0Gg0BlFdX03yP5MDnoJHO1zghP4ncNHoi4yWodFoDCI6PJqhvYcGPNFkSDnc1tKk1zvqmZc5j9EDRgdNh1nuhmsd3iQlBSZYUUcxiz1CTcestFksKloU0D5CyuG2ttPs0e8e5bYvbwuqjo5mnggUWoc3vobvDDRmsUeo6Zgzcg79ovsFtI+QcrhOZ8tJ8UVFizht6GlB1VFVVRXU/tpC6/DmWMsGg4lZ7BFqOlL7pnLPtHsC2kdIOdzme+UPVB8gvyyf01NOD6qOuLi4oPbXFlqHNz179jRaAmAee4SijifXPMl9S+8LWPuGOlwROVNEikRkq4i0+F0vIr8XkU0isl5EvhKRoR5lDhFZ5z58mnhpvrX3YO1Bbpt0G1HhUV1+LR0hVLZM+opZdOitvd6Eoo6MARm8s/mdgLVvmMMVESvwOHAWMAq4RERGNau2FshWSp0EvA383aOsRik1xn10ak3X8D7D+eNpf+zMpRqN5jhkYtJEfjz8IzsP7wxI++06XBF5wJdznWAcsFUptU0pVQe8Dsz2rKCU+kYp1TgsXQX4LV95naOO1EdTOVrnW9odjUZz/BNmCeO32b8NWMhGX0a4P2vl3Fl+6Hsw4PlvpNR9ri3mA596PI8UkVwRWSUiczra+ZLtSxjQYwA9Inp09FKNRnMc85fpf2FC4oSAtN3m1l4R+S1wHTBcRNZ7FMUCywOipm0tlwPZwBSP00OVUrtEZBjwtYhsUEq1SFAkItcA1wAMGDCAxYsXA7DnyB7+J/5/+P7770lJSaG4uJiMjAyWLl3aov+cnBxKSkpITEykvLycnTu9f24MHjyYxMRESkpKSEtLY9myZS3amDRpEsXFxaSkpNDQ0NCko5GkpCTi4+MpLS0lJSWFlStXtmhj8uTJFBYWkpaWRklJCXv2eAdKT05OJi4ujvLychITE1m9enVzWzBlyhQKCgrIyMigvr6+hY5hw4YRGxtLRUUF8fHxLRL4hYWFMWnSJPLz8xkzZgyFhYUt9rqnpqZis9moqqoiLi6uxRyczWYjJyeHvLw8xo4dS11dXQsd6enpWK1WamtriY2NpaCgwKs8KiqK8ePHN7WRn5/PkSNHvOqMGjUKh8OB0+nEZrOxceNGr/KYmBiys7PJy3OF5mtNx+jRo7Hb7VgsFqxWK5s2bfIq79mzJ1lZWU06Vq9e3WK1Q2ZmJpWVlURGRuJwOFrMFffp04fMzMymNmpqalroyMrKoqKigpiYGOx2O1u2bPEq79u3LxkZGaxbt46srCyWLVvWInZIdnY25eXlxMXFUVlZybZt27zKBwwYQFpaGoWFhWRmZnL06NEWOsaPH09paSnx8fFUVFSwfft2r/KBAwf6/fvUmg7P71NpaSm7dnmPSLv6fSo6UMTUzKkkDEho+j75BaVUqwfQC0gGXgOGehxxbV3TkQPIAT7zeH47cHsr9WYAm4EBx2jrBWBue32OHTtWNXLnN3eq3F25SqNRSqmcnBy1fPlyo2VoTMK4Z8appduXep0DclUX/V6bUwpKqcO4fvKfrJTa4XFU+MfVswZIFZEUEYkALga8VhuIyMnAU8AspdRej/N9RMTmftwPOBXwHnq0gucqhbum3sXYQWP98To6TOOIymi0Dm/MkibdLPYIZR1jE8aSV+b/fo85h6uUcgBFIjLE3x0rpRqAG4DPcI1g31RKFYrIPSLSuOrgQSAGeKvZ8q8TgFwRKQC+Ae5XSrXrcBsDkG/cu5HffPgb/76gDjB2rDGOvjlahzcjR440WgJgHnuEso6JSRM5VHvI7+36Ep6xD1AoIt8BTbf0VSeXYnmilPoE+KTZuTs8Hs9o47oVwIkd7a9xbm1V6SpqGozbVVRQUEBmZqZh/WsdrbNlyxZT5DQziz1CWcflJ10ekHZ9cbh/DkjPBtB4EyFvdx5jE4z7733wYOADHfuC1uFNZWWl0RIA89gj1HX87r+/497p9xITEeO3Ntt1uEqpJX7rzSTsr9lP9qAuJd/UaDTHOSt2rmDdnnVMGjLJb236svFhgoisEZEqEalzb6k90t51ZuatC97i1CGnGi1Do9GYmLEJY/2ezdeXjQ+PAZcAW4Ao4GpcW3K7JVsObOGx7x4zWoZGozE52YOyWbtnrV/b9CmWglJqK2BVSjmUUs8DZ/pVRZCIjIxkyY4lfLfrO0N1pKenG9p/I1qHN0OG+H0xTqcwiz1CXcdlJ13GM+c849c2fXG41e51sutE5O8icpOP15kOETH8hhm0DBNpFFqHN1qHN6GuIzIsktc3vu7XeCu+OM4r3PVuwLUsLAk4328KgojT6SR/T75hGx4aqa0NbKI6X9E6vLHb7UZLAMxjD63DlRHGn9MKvqxS2OEe4SYD7wJFyhXdq9thtVr5/PLPiQ6PNlRHbGysof03onV407gxxmjMYg+twzWPm7c7z28rFXxZpfBz4Afg37huoG0VEX9ECws6R6uPkrs7l3CrsbmrmgdiMQqtw5utW7caLQEwjz20Dv9v8fVl48PDwDT3jTNEZDjwMd6hErsFTuXkX6v/xenDgptSR6PRdE9mj5xNTlKO39rzZQ63stHZutkGmGNLTgdRSpEQk2C0DI1G003oHdkboWW2787ii8PNFZFPROQqEZkHfAisEZHzROQ8vykJAgrFoNhBRsvQaDTdhJr6GrKfyW4MA9tlfJlSiATK+Sn49z5cGyDOARSuG2ndgjBrGBePvthoGURFBTdpZVtoHd7YbDajJQDmsYfWAbG2WCxi4YjdP5trfVml8Eu/9GQCoqOjGR433GgZjB8/3mgJgNbRnIyMDKMlAOaxh9bhIiEmgbKqMr+01S03MHSWfYf2UbDH+DuvoRzYuTXMokMHIPdG63Bx04Sb6BHun9yH4q+5ie5AeFK42rFph57H1bRg4sSJPPTQQ6aIh6sxJyKSp5TqUphBX9bhmmN/nx+Is8YxoMcAo2W0SKpoFFqHN82TOxqFWeyhdbi4b+l9PLr6Ub+05cuUwhYReVBERvmlRwPpEdaDMIsv9wkDS/PsskahdXhz9Kj/9sx3BbPYQ+twYQuzUXKoxC9t+eJwM4Fi4FkRWSUi14hIT7/0HmSM3mGm0Wi6H0G9aaaUqlRKPaOUmgjcCtwJlInIQhEZ4RcVQaLeUW+0BI1G080YHjecvlF9/dKWT3O4IjJLRN4D/olrq+8wXBsgPjnWtWZDETo3CDUajX+YkDiBx872T9ICXyY0t+BKRf6gO1tuI2+LyGS/qAgSDZYGoyUAMGqUOabDtQ5vUlJSjJYAmMceWoeLOkcd1398vV/aOqbDda9QeEEpdU9r5Uqp//GLiiBhhhtmAA6Hw2gJgNbRHK3DG63DRbglnOfWPeeXto45paCUcgC/8EtPJqCnzRz3+pxOp9ESAK2jOWZZk24We2gdLkSEcIt/brj7MuRbLiKPAW/gyvgAgFLKHIv0OoDdYY6I/mbZs691eBMebo5VLGaxh9bxE+HWcOx03X/4sixsDJAB3IPrhtnDwENd7hkQkTNFpEhEtorIba2U20TkDXf5ahFJ9ii73X2+SERm+tJfvd0cqxQ2btxotARA62jOtm3bjJYAmMceWsdP7PjdDr+048uysGmtHNO72rF7fvhx4CxgFHBJK5sr5gMHlVIjgEeAB9zXjgIuxvWP4EzgCV92xPkzrmV35pVXIDkZpk+fQnKy67lG08jx+vlofF0WC62+rgcffLBN575i54pWz3cUn4LXiMjPReQPInJH4+GHvscBW5VS29w50l4HZjerMxtY6H78NnC6iIj7/OtKKbtSqgTY6m5P0w6vvALXXAM7doBSwo4drufHy5dK0zWO18+H9+ui1dcVHR3NjBkzuPDCC1s43ms+vMYvOtqdwxWRBUA0MA14FpgLfOeHvgcDOz2elwLN47A11VFKNYjIYaCv+/yqZtcObq9Dq8XK+++/zyefGLt8ePfu3bz66quG9P3aa1Bd7X2uuhquvRaWLDFEkqH2aOSHH37g1VdfZf369YbqAP35aI4/7OHr6zrrrLPYvHkz2dnZpKens27dOkTEbyucfGllolLqJBFZr5S6W0QephvlMxORa4BrAAYMGMD+/fu9soD26tWLvn37snfvXhISEtiyZUuLNoYNG8aBAwfo3bs3lZWVVFRUeJX37t2bPn36sH//fuLj41tNRjhixAjKy8vp169fq1lI4+LiiI2N5dChQ/Tt27fV+cTU1FTKysoYMGAABw4c4PDhw17lffv2pUePHhw5coQ+ffpQUuK9/1tEqKpKa9VOVVWK2Nhi+vfvj81mo7q6mtjYWHbs8J67slqtjBgxgh9//JHExETKysqoqqryqjNgwADCw8Opra2lR48e/Pjjj17lYWFhDB8+nB07djB06FA2b96MxeL9Y2vgwIGICA0NDdhsNkpLS73KIyIiSElJaWrjxx9/pKamxqtOQkICSimcTifh4eHs2rXLq9xms5GcnMyOHTuIiYkhMTGxxXszePBg6uvrsVgsiAhlZd5bPKOiohgyZEiTjpKSEurqvJNaJyYmYrfbCQsLQynFnj17vMqjo6NJSkpqamPjxo0tbuANGTKEo0ePEhkZSX19PXv37vUqj4mJISEhgdLSUoYMGcLWrVtbLKcaOnQolZWVREdHY7fb2bdvn1d5bGwsVVUJ0MrUW+PnIyUlhYMHD9KzZ0+OHj3KgQMHvOoF4vtks9mIiIjwquP5fTp48CCHDh3yKm/+fXrmmbA2X1dmZkPT92n//v0cOnQIq9VKeno6lZWVlJeX08vai51e48NOopQ65gGsdv9dBQwCbLimAtq9tp12c4DPPJ7fDtzerM5nQI77cRiw3201r7qe9Y51pIxKUaHO0KFKuX5UeR9DhxqtzFhycnLU8uXLjZZhOMfr58OX1zV//nw1YMAAdf/996vKykqv61/f8LoCclUX/Z4vc7gfiUhv4EEgH9gOvNZlTw9rgFQRSRGRCFw3wRY1q7MImOd+PBf4Wiml3Ocvdq9iSAFS8WGaQ+rNcdPMyIDK990H0dHe56KjXeeNwugA043oAOTH7+fDl9d10UUX8cMPP3DrrbcSExPjVXfOyDld1gC0P8L1PHCNbnt11ct7tHc2rkhkPwB/cp+7B5jlfhwJvIXrpth3wDCPa//kvq4IOMuX/oae4PHvLIR5+WXXf3YR19+XXzZakfHoEe5PHK+fj668rl5/6+WXEa6vjnEicClwZePR1Y6NOFLTU323cABZs2aN0RKUUlqHJzk5OerZZ581WoZSyhz2UErr8KTHfT384nB9WaXwEjAcWAc0zsIr4MWujKyNQJQ5phSa32QyCq3Dm+Y33YzCLPbQOn6i3umfTVO+rFLIBkYppcyx0bwLNDjNES1Mo9F0L2anz+Yt3upyO77cNNsIDOxyTyZA6Xi4Go2mE7x5wZt+aceXEW4/YJOIfAc/RW9QSs3yi4IgchwM0jUaTZApPlDMwyse9ktbvjjcu/zSkwmos9S1XykIjB492mgJgNbRnGHDhhktATCPPbQOFyUHS9h2yD+Bjdp1uEopgzb0+Z9+Uf2MlgCA3W6OMJFahzf19eaIJmcWe2gdLsqqykiISfBLW206XBFZppSaJCKV4DX5KYBSSpkjmncH2Ht0b/uVgkDzbaxGoXV444qLZDxmsYfW4eJw7WGSeib5pa02Ha5SapL7b8uN/92UI/Yj1DnqiLBGtF85gFit7UaSDApahzdahzdah4v/nfC/APyVv3a5LV+y9sa1cpgjNH4H6RvRlz1Ve9qvGGA2bdpktARA62hO82A/RmEWe2gdLt4qfIttB/0zh+vLWD0f2IdrC+4W9+PtIpIvImP9oiJICMLuyt1Gy9BoNN2If67+p9/8hi8O9wvgbKVUP6VUX1wZGj4CrgOe8IuKIBEZFsn4wc1D7mo0Gk3b7K7c7bebZr443AlKqc8anyilPscVCnEVrmA23QancrK+3PgA0xqNpnuglKKssoyE2OA53DIRuVVEhrqPPwDl7hxi5sij7CP1qp63NnV9e15X6dnTHAs8tA5vevToYbQEwDz20Dpcu1MXX7WY6PDo9iv7gC8O91IgEXjffQxxn7MCF/pFRbAIxxRzuFlZWUZLALSO5qSnpxstATCPPbQOqK6vpm9UX7+150vW3v1KqRuVUie7jxuUUvuUUnVKqZa5ZExMmCPMFA7XLAG3tQ5vdAByb7QOWFW6ims/vtZv7fkSnrE/8AdcKckjG88rP6RKDzYxPWL48+Q/Gy2DsWPNsbhD6/Bm5MiRRksAzGMPrQM279tMalyq39rzZUrhFeB7IAW4G1eKnTV+UxBE7DV2UvqkUFNvbNzT1atXG9p/I1qHN4WFhUZLAMxjD60DcstyGZvgP4fvi8Ptq5T6D1CvlFqilPoV0O1GtwBOp5NL37mUb3/81lAdZgl0rXV4Y/Se/UbMYg+tA84ecTanDzvdb+354nAbI3qUicjPReRkIM5vCoJM9qBs8nabY25Ko9GYF6dyckHGBQzr478ocr443HtFpBdwM/B/wLPATX5TEGTGJowltyzXaBkajcbkrCpdxdQXpvq1TV/CM37kfngYmObX3g3gtKGnccR+xGgZGo3G5OTtzmNU/1F+bdOXVQopwI1Asmf97pjxITo6msSeifwm+zcopQwLx5eZmWlIv83ROrwZMWKE0RIA89gj1HXkluUyKWmSX9v0ZUrhfVwrEx4FHvY4uh0Ohyvp8HlvnMcX274wTEdlZaVhfXuidXhTXV1ttATAPPYIdR1Dew3l1CGn+rVNX1Ls1Cql/u3XXg2iMZDxsD7DyNudxxnDzzBER2RkZPuVgoDW4Y3NZo7QIGaxR6jruGfaPX5v05cR7r9E5E4RyRGRrMbD70qCQGMSSaNvnDWOtI1G6/BG6/AmlHWs2bWG3370W7+364vDPRH4NXA/P00nPNSVTt1BzL8QkS3uv31aqTNGRFaKSKGIrBeRizzKXhCREhFZ5z7G+NJvbW0tAOMTx/st3FpnKCoqMqxvT7QOb3788UejJQDmsUco61hVugqn8n9sLl+mFC4Ahiml/Jny9jbgK6XU/SJym/v5rc3qVANXKqW2iMggIE9EPlNKHXKX36KUersznQ/rM4zHzn6ss9o1Gs1xTiBumIFvI9yNQG8/9zsbWOh+vBCY07yCUqpYKbXF/Xg3sBfo7y8B/+/r/8e3O4zdcabRaMzJEfsRsgdl+71dXxxub+B7EflMRBY1Hl3sN14pVeZ+vAeIP1ZlERkHRAA/eJy+zz3V8IiIdPhuhyB8vOXjjl6m0WhCgPcueo+TE072e7u+TCnc2ZmGReRLYGArRX/yfKKUUiKiWqnX2E4C8BIwT6mmSZXbcTnqCOBpXNMRrd5SFJFrgGsA4uPjWbx4MQA5jhxKDpXw/fffk5KSQnFxMRkZGSxdurRFGzk5OZSUlJCYmEh5eTk7d+70Kh88eDCJiYmUlJSQlpbGsmXLWrQxadIkiouLSUlJwWKxNOloJCkpifj4eEpLS0lJSWHlypUt2pg8eTKFhYWkpaVRUlLCnj3eCTGTk5OJi4ujvLycxMTEFkE/RIQpU6ZQUFBARkYGQAsdw4YNIzY2loqKCuLj48nN9b65GBYWxqRJk8jPz2fMmDEUFhZy4MABrzqpqanYbDaqqqqIi4sjPz/fq9xms5GTk0NeXh5jx47F6XS20JGeno7VaqW2tpbY2FgKCgq8yqOiohg/fnxTG/n5+Rw54r2hZdSoUTgcDpxOJzabjY0bN3qVx8TEkJ2d3RT+r6GhoYWO0aNHY7fbsVgsWK3WFgkNe/bsSVZWVpOO1atXt9j7n5mZSWVlJZGRkTgcjhZzkn369CEzM7OpDbvd3kJHVlYWFRUVxMTEYLfb2bJli1d53759ycjIYN26dWRlZbFs2TIaGhq86mRnZ1NeXk5cXByVlZVs2+adGHHAgAGkpaVRWFhIZmYmNTU1LXSMHz+e0tJS4uPjqaioYPv27V7lAwcO9Pv3qTUdnt+n0tJSdu3a5VXele/T4drDOHFy8siTvb5PfkEpFfQDKAIS3I8TgKI26vXElcRy7jHamgp85Eu/Y8eOVY04nA51ytOnqAPVB5RGk5OTo5YvX260DI0JOOfVc9Qr619pcR7IVV30fW1OKYhIpYgcaeWoFJGu7o1dBMxzP54HfNBK/xHAe8CLqtnNMfeoF3FtFZuDa565XTwXtlvEwne//o64qODH4dGBnb0xiw4dgNybUNRRXV/N4u2LOWvEWQFpv02Hq5SKVUr1bOWIVUp1NcnQ/cDPRGQLMMP9HBHJFpFn3XUuBCYDV7Wy/OsVEdkAbAD6Aff60ml0tHdeou2HtnPzZzd38aV0HB3Y2Ruz6NAByL0JRR35ZflMSJxAn6gWK1X9gi83zfyOUuqAUup0pVSqUmqGUqrCfT5XKXW1+/HLSqlwpdQYj2Odu2y6UupEpdRopdTlSqkqX/o9evSo1/P4HvE8u/ZZKmoq/PsC26G1+SQj0Dq8aT7HaxRmsUco6pg0ZBKfXvZpwNo3xOEahdPpvZA5KjyKacnT+HRL4AzcGmYJdK11eFNX58+l5p3HLPYINR1O5eSWz2/BoQK3sy2kHG5rnJN2Dvll+e1X1Gg0xzVrdq3hk62fEGGNCFgfviwLO6751cm/MixMo0ajMQ+LihYxKy2wUWdDfoQrIvwn/z+sKl1ltBSNRmMgX5Z8yaz0wDrckBrhNl+l0Mjeo3sp2FDAhMQJQdGRlWWOYGtahzfp6elGSwDMY49Q07HkqiUBnU6AEBvhthXmbVb6LBYVLWoK3xhoKiqCuyqiLbQOb5rvVDMKs9gjlHR8tvUzCvYUYJHAusSQcriNAcibM6r/KGxhNkoOlQRFR0xMTFD6aQ+tw5uoqCijJQDmsUco6bh/+f2UVZW1X7GLhNSUQlsjWBGh4NoCIsOCE1k+1JbbtIdZdNTX1xstATCPPUJFR/GBYjbt2xSw3WWehNQItzEAeWvYrDZu+u9NQcno2zzwiFFoHd40D6JiFGaxR6joeHXDq/xqzK+whQU+xVJIjXCPhYhQWlnKy+tf5rpTrjNajkajCRJ/nvxnahvaHoz5k5Aa4bbHddnX8cSaJ4J280yj0RjLh0Uf8uW2L+kR0SMo/WmH68HU5Kn0iOjBtoPb2q+s0Wi6NUop7lpyV0Byl7VFSE0phIUd++WKCKvmrwr4zrO+ffsGtH1f0Tq86dWrl9ESAPPY43jXsWb3Gg7WHGTmiJkBab81QmqE68uyHxHhyveupKwycEtEGrMtGI3W4U1KSorREgDz2ON41/Hplk/5bfZvA7721pOQcrieAciPRVRYFM/mP9t+xU6ybt26gLXdEbQOb8xyV94s9jjeddw59U5+n/P7gLTdFiHlcNva2tuc6065jqfzn6bB2dB+5U4Qalsm28MsOvTWXm+OZx3PrX2ORUWLsFqsfm/7WISUw62q8ilOOZkDM8lJzAnYzbPWkkwagdbhTfNElUZhFnscrzqcysl9395HfI9jJgsPCCF106wjy73evODNgOlonk3VKLQOb9qKtRFszGKP41XH5z98Tu/I3owbPM6v7fpCSI1wO8q5b5xL0f6i9itqNJpuQ35ZPjeccoMhcbC1wz0GWQOzuGfpPUbL0Gg0fkIpxR9P+yO/PPmXhvSvHe4xuCnnJr4u+Zp1e9YZLUWj0fiBOW/MYXXpasP6DymH26NHx7bvxUTE8I8z/sGh2kN+1ZGdne3X9jqL1uGNWdKkm8Uex5uOL7d9yeZ9m8lKMG71RUg53M6E37vkxEs4bchpft0IUV5e7re2uoLW4c3BgweNlgCYxx7Hkw6lFLd9eRv3Tr+XcGu4H1R1jpByuO1t7W2LD4o+4Nw3zvVbUJu4uDi/tNNVtA5vYmNjjZYAmMcex5OOmoYaZqfPZu6ouX5Q1HlCyuF2dtnPnJFzqGmo4YOiD/yio7Ky0i/tdBWtwxtfdyIGGrPY43jRUe+op6yyjD9P+XNQt/G2hiG9i0iciHwhIlvcf/u0Uc8hIuvcxyKP8ykislpEtorIGyLiU+a3zkaOt4iFv53+N/741R9xOLu+VnPbNnNEI9M6vNm9e7fREgDz2ON40fHc2ue47hNzxLg2yt3fBnyllEoFvnI/b40apdQY9+GZv/gB4BGl1AjgIDA/sHLhrBFn8cw5zxj+H1Kj0fhOdX019yy9h3un3Wu0FMA4hzsbWOh+vBCY4+uF4lqtPB14uzPXdxYRYWLSRJ7Keypo0eE1Gk3XePy7x5mYNJFTBp9itBTAOIcbr5RqvO2/B2hrU3OkiOSKyCoRmeM+1xc4pJRq3O9XCgwOnNSfEBH+u/W/PP7d48HoTqPRdJH5WfP558x/Gi2jiYDFUhCRL4GBrRT9yfOJUkqJSFu3/4cqpXaJyDDgaxHZABzuoI5rgGsA4uPjWbx4sVf5wIEDSUlJobi4mIyMDJYuXdqijZycHEpKSkhMTOSWE24h9/tcvvrmK6ziijQ0ePBgEhMTKSkpIS0trdVgG5MmTaK4uJiUlBTCwsJa6EhKSiI+Pp7S0lJSUlJYuXJlizYmT55MYWEhaWlplJSUsGfPHq/y5ORk4uLiKC8vJzExkdWrvRd4iwhTpkyhoKCAjIwMLBZLCx3Dhg0jNjaWiooK4uPjyc3N9SoPCwtj0qRJ5OfnM2bMGAoLCzlw4IBXndTUVGw2G1VVVcTFxZGfn+9VbrPZyMnJIS8vj7FjxwK00JGeno7VaqW2tpbY2NgWgWWioqIYP358Uxv5+fkcOeKdAHTUqFE4HA6cTic2m42NGzd6lcfExJCdnU1eXh4ATqezhY7Ro0djt9uxWCxYrVY2bdrkVd6zZ0+ysrKadKxevZqamhqvOpmZmVRWVhIZGYnD4aCoyHu7eJ8+fcjMzGxqo76+voWOrKwsKioqiImJwW63twgl2bdvXzIyMli3bh1ZWVksW7asRQyC7OxsysvLiYuLo7KyssXc6IABA0hLS6OwsJDMzExqa2tb6Bg/fjylpaXEx8dTUVHB9u3bvco7+n0qLy9vkbyz+fepNR2e36fS0lJ27drlVZ6UlMSKihWMCBtB1qisFtdDx79PfkEpFfQDKAIS3I8TgCIfrnkBmAsIsB8Ic5/PAT7zpd+srCzlD656/yr1UsFLnb6+vr7eLzq6itbxEzk5OWrJkiVGy1BKmcMeSnVvHSUHS1TcA3Gq9HCp33QAuaqLvs+oKYVFwDz343lAi/VWItJHRGzux/2AU4FN7hf+DS7n2+b1rXGsNOkdYcHPF3D5SZd3OhdSYWGhX3R0Fa3Dm5KSEqMlAOaxR3fV4VRO5i+az62n3srgnkGZbfQZoxzu/cDPRGQLMMP9HBHJFpHGVAsnALkiUoDLwd6vlGr8PXcr8HsR2YprTvc/vnTqS4odX7CF2dhyYAsTnp3QqRtomZmZftHRVbQOb1JTU42WAJjHHt1Vx/ZD2+kR3oObc24OkKLOY4jDVUodUEqdrpRKVUrNUEpVuM/nKqWudj9eoZQ6USmV6f77H4/rtymlximlRiilLlBK+bTA1p8LuUfEjSCxZyJ3L767w9cuWbLEbzq6gtbhzdq1a42WAJjHHt1Rx6HaQyT3TmbRJcHP5uALelFpJxERnvz5kzy37jm+2/Vdh65Vftoi3FW0Dm+0Dm+6mw6ncnLeG+fx6oZXA6yo82iH2wXiY+J564K3GNJriNFSNJqQ56ncp6iur+aS0ZcYLaVNtMPtIpOHTsbhdLBw3cL2K2s0moBQVlnGHYvv4PnZz5tyKqER7XD9QJgljD98+QfW7FpjtBSNJiQZGDOQJVct4YT+Jxgt5ZiElMPtaAByX4mPieefM//JVR9c5dOqhfHjxwdER0fROrzJyMgwWgJgHnt0Fx3Pr32ej7d8zKj+o4KkqPOElMPtTAByX7l49MWMGzyODeUb2q1bWloaMB0dQevwZu/evUZLAMxjj+6go+RgCbd8cQvD+wwPoqLOE1Jp0jsbgNwXRITnZz+PUooD1QfoG923zbrx8W2FjgguWoc3Zgm4bRZ7mF1H4waHP5z6B9NPJTQSUiPczgYg7whLdyzltOdPo6quqs06FRUVAdfhC1qHN81jMRiFWexhdh2V9koy4zNNucGhLULK4XY2AHlHmDx0MjmJOVz53pVtbv1tHvDDKLQOb8rK/Je3riuYxR5m1rG6dDW1DbU8cuYjpl6V0JyQcrjBQER44udPUH60nNc3vm60HI3muGPbwW3Mfn02m/Ztar+yyQipOdxgYQuz8cmlnxBri6XOUUeE1acMQBqNph0q7ZXMfn02/2/y/2NayjSj5XQYPcINEL0ie1HnqCNzQSYFewrav0Cj0bRL7u5cpiVP4/pTrjdaSqcIKYcbHh7cfPSRYZHcNeUuZr8+m71Hf1pyNHBga3HZg4/W4U3fvm2vLAkmZrGH2XTk7c5jWso0/n3Wv3Fl2up+hJTDtdlsQe/zotEXcflJlzN/0U95LlNSUoKuozW0Dm8SEhKMlgCYxx5m0vFm4Zuc9+Z5VNrNkbq9s4SUw/VXAPKOcs+0e3jsrMeaor4XFxcboqM5Woc3zVO9GIVZ7GEWHcvXLuf6T67ng4s/INYWa7ScLhFSN838FYC8o1jEwtDeQ7n9y9tJ6pXEtWOvNURHc8yyldUsOswyojOLPcyi4+Xyl3ni7CcYM3CM0VK6TEiNcP0ZgLwzXJ11NXcvuZuPv/zYUB2NtJbgzwjMomPdunVGSwDMYw+jddgb7JRVlnFlryu5IOMCQ7X4i5ByuEYzPG44r573KiUHS9hdudtoORqNaVFKcd3H13Hn4juNluJXtMMNMqcPO50RcSNIiEkwTUR9jcZMKKW47cvbWFe+jofPeNhoOX5FO1wDiA6PZnflbk555hTKKs2xnVSjMQvbDm5jyY4lfH75593+JllztMM1iME9BzNn5BxmvDSDfUf3GS1HozEFy35cxrA+w1g5f+UxI+51V0LK4cbExBgtAYCcnBwA/t/k/8e5I8/lyvevNFSH0ZhFx+jRo42WAJjHHsHW8a9V/+Kq96/iiP2I18YGs9jDH4SUww1GtDBfKCkpaXr8l2l/YeGchdQ76jliD254QE8dRmIWHWaJFmYWewRTx4LcBTyy6hG+uvIrekX2MkxHoAkphxsRYY4gMomJiU2PRYQBPQbw6oZXmbZwGgeqDxiiw0jMoqN///5GSwDMY49g6VBKUbCngK+u/IqhvYcapiMYhJTDDWSKnY5QXl7e4tyVmVcyI2UG0xZO84q7EGwdRmAWHQcPHjRaAmAeewRDxxNrnmDbwW08+YsnGR7Xepocs9jDHxjicEUkTkS+EJEt7r99WqkzTUTWeRy1IjLHXfaCiJR4lI3xpd+6ujr/vpBO0toWUhHh/hn3c+7Ic3mx4EXDdBiBWXSY5YttFnsEUodSinuW3MO/V/+byLBIw3QEG6O29t4GfKWUul9EbnM/v9WzglLqG2AMuBw0sBX43KPKLUqpt4MjNziICHdPuxulFCt2riCxZyJDeg0xWpZG43fuWXIPb216i8VXLWZgjDmikgUDo6YUZgML3Y8XAnPaqT8X+FQpVR1IUWZBRFhbtpbJz0/2KQuwRtNdaHA2oJRixrAZIedswTiHG6+UarwlvAdoLz3oxcBrzc7dJyLrReQREQl+3MUAc/2467lv+n1Mf3E6n239zGg5Gk2XOVB9gDNeOoNXN7zKqUNOpV90P6MlBZ2ATSmIyJdAa/++/uT5RCmlRKTNPa4ikgCcCHh6ndtxOeoI4Glc0xH3tHH9NcA14ApkvHjxYq/ygQMHkpKSQnFxMRkZGa0G7MjJyaGkpITExETKy8tbzCkNHjyYxMRESkpKSEtLY9myZS3amDRpEsXFxaSkpGCz2VroSEpKIj4+ntLSUlJSUli5ciWDGczCMQthO3xT+g2TT5vMpk2bSEtLo6SkhD179ni1kZycTFxcHOXl5SQmJrJ69ermtmDKlCkUFBSQkZFBWFhYCx3Dhg0jNjaWiooK4uPjyc3N9SoPCwtj0qRJ5OfnM2bMGAoLCzlwwHtlRWpqKjabjaqqKuLi4sjPz/cqt9ls5OTkkJeXx9ixY7FYLC10pKenY7Vaqa2tJTY2loIC76wZUVFRjB8/vqmN/Pz8Fll3R40ahcPhwOl0YrPZ2Lhxo1d5TEwM2dnZ5OXlAa55xeY6Ro8ejd1ux2KxYLVa2bTJO49Wz549ycrKatKxevVqampqvOpkZmZSWVlJZGQkDoeDoqIir/I+ffqQmZnZ1EZDQ0MLHVlZWVRUVBATE4PdbmfLli1e5X379iUjI4N169aRlZXFsmXLaGho8KqTnZ1NeXk5cXFxVFZWsm3bNq/yAQMGkJaWRmFhIZmZmdjt9hY6xo8fT2lpKfHx8VRUVLRI8NjW96m2oZatFVuZHzefc1PP5fvvv/f5+9SaDs/vU2lpKbt27fIqb+371JzJkydTWFjo8/fJLzTGaA3mARQBCe7HCUDRMer+L/D0McqnAh/50u/JJ5+szEB1dXWHr/nz139WF7x5gaqyVxmqIxCYQUdOTo76+uuvjZahlDKHPZTyr47rPrpOvbjuRcN1dAUgV3XR9xk1pbAImOd+PA/44Bh1L6HZdIJ71Iu4tqPMATa2vKwlZlml0JmF3H887Y9EhUcx6flJ/Hj4R8N0BAKz6Ni92xwR3Mxij67qUErx9+V/5/v93/P4zx/niswrDNFhJoxyuPcDPxORLcAM93NEJFtEnm2sJCLJQBKwpNn1r4jIBmAD0A+415dOjUix0xppaWkdviYyLJIXZr/A5Sdezjub3jFMRyAwi46kpCSjJQDmsUdXdFTXV3PJO5fwzuZ36GnraZgOs2GIw1VKHVBKna6USlVKzVBKVbjP5yqlrvaot10pNVgp5Wx2/XSl1IlKqdFKqcuVUlW+9FtV5VO1gNPaHK8viAg3T7yZm3Ju4tMtn/Kf/P8YosPfmEXH+vXrjZYAmMceXdFx7hvnEmGNYMlVSxgUO8gwHWYjpHaaHU8M6zOMv6/4O//76f/S4Gxo/wKNJghs3LsRp3Ly9C+eZuGche1uagg1tMPtpqT3S2f11aspOlDE/376v0bL0YQ4SimezX+W6QunU3ygmKG9h3bbVOaBJKSSSB5v9I7szceXfsyBmgMcrj3Mkh1LmJU+y2hZmhCjtqGWi9++mK0VW1n6y6WM7DfSaEmmRY9wuzlWi5UBPQawu3I3//f5/3HJO5ewv3q/0bI0IYBSiq0VW4kMi+TCjAvJuyZPO9t2EBVCebXGjh2rGhe5G0lDQwNhYf7/cVFdX82fv/4zb29+m03XbaJHRA9DdHQUM+iYOHEi999/P5MnTzZUB5jDHu3p2FO1h2s/upbyo+Us/9VyLBK4sZtZ7CEieUqp7K60EVIjXLMEIC8uLg5Iu9Hh0Tw882FW/GoFPSJ68GLBi8cc7QZKR0cxiw6zRKUyiz3a0vFNyTdkLshk9IDRLJ63OKDO9lg6uiMh5XDNEoA8JSUloO0P7jkYpRTry9dz0pMn8d7m9wzR4Stm0TFoUNeWL/kLs9ijuY49VXsoqywjvV86H1/6MfdOvxdbWODXtpvFHv4gpByuWQKQl5aWBrwPEeGhMx7irQve4tYvb+Xl9S8bosMXzKJj797gBH5vD7PYo1GHUopX1r9C5oJMvi75mkGxg8ge1KVf1p3ScTwQUg7XLFt7mwfaCCSnDjmVddeu4/wTzue7Xd95jXaDqeNYmEXHvn3myJ5sFns06rjy/Sv527K/8fGlH3PZSZcZpuN4wPiZaE3AiQ6Pbnp865e38vy65/nr6X81UJHG7NTU13Cw1pVy6IZTbmDMwDFBmT443gmpEW6oM27wONb/dj3Tkqfx6w9/DYDD6TBYlcZMNDgb+E/+f0h7LI2Kmgpq6msYnzheO1s/oR1uiBEZFslNOTex4lcrADj9xdO56b83se+oOX5Oa4yhMXzg6xtf5+UNL/P2BW8zvM9wosKjjJZ2XBFSDtcsqxTMEJVKREhKSuKNuW9Q76xn5OMjuxwMp7OYwR4A8fHtJR4JDsG2xzcl3zDhPxP47IfPuPTES/n6yq8ZnzjeNO+LWXT4g5ByuOHh4UZLAMzzxY6Pjyc+Jp7Hzn6M767+jsyBmdgb7DyT9wx1juDdYDSLPfr0aZE82hCCZY+KmgpmvjyTqz+8mt+N/x1nDD8Di1iaYiCY5X0xiw5/EFIO1yyrFMyyzMVTx/C44WQPyuZg7UHe/f5dRj42klfWv4LTOzJmwHUYiVlWKQTaHlsrtrJ4+2J6R/bmypOuZPP1m7nkxEtabGAwy/tiFh3+IKQcrlkCkJtlIXdrOgbGDOTTyz7ludnP8diax1izaw1H645S21AbVB1GkJCQYLQEIHD2KNhTwK8X/ZoJz05g877NWMTCZSddRoS19ak2s7wvZtHhD0LK4ZolAHlrCe2M4Fg6piZPZcWvVjA+cTwfFX/EkEeGcOsXt7Lt4LY2rwmEjmDSPMmkUfjTHvYGOwV7XAk4H1r5EEN6DeH7G77nt6f8Nqg6uoJZdPiDkHK4mo7ROJd30eiLWDF/BQ7lYMKzE9hTtYequiq9pMzE7Di0gz9+9UeG/HMI9y+/H4CXzn2JP0/5c0imJzcL2uFqfGJE3AgeOuMhdv1+FwNjBvJU7lOMeHQEDyx7QC8pMwlO5WTjXtco/YHlD1BTX8PSq5by2vmvtXOlJljonWaaDhFuda30uHnizUweOpknc5/kpAUnse1/tuFUTqLDo3Wk/yBzoPoAz697ngW5CxgYM5Clv1zKEz9/wmhZmlbQDlfTaU4ZfAqnDD6F2oZaIsMi+d1/f8eSHUu44qQrOCftHFL7phot8bjlaN1R1uxew9TkqTy88mFKj5TyynmvMG7wOP0Pz8SEVADy7OxslZuba7QMnE4nFovxszn+1uFUTr7c9iVvb3qblaUrWfebdeTuzqXB2cCExAlYLdag6OgMEydO5O9//zuTJk0yVAcc2x6LihbxTP4zLNm+hAmJE/j40o+bfnUEU0cwMYsOHYC8g9TU1BgtAYDCwkKjJQD+12ERC2cMP4Onz3maDb/dgNViZWvFVq775DoSHk7glx/8kur66oDr6CwlJSVGSwB+skdjTON7l97LhW9dCLj+qV124mX8eNOPfH7F5wFztp46jMYsOvxBSE0pREaaI2VzWlqa0RKA4Oi47KTLuOyky9h+aDtfbvuSqLAo/rXqX3z2w2fMSp/FOWnnmMYeZthCWueog76ux/Pen8eyH5cxK30W12Zfi1KKOSPnBE2LWd4Xs+jwByE1wjVLih2zjKSCqSO5dzJXZ12NiHDVmKu4asxVLPtxGZkLMtlYvJHCvYW8sO4FNpRvoMHZEDRdnpSVlQW9T3uD6zP5+HePM+f1OcQ/FM/bK94G4LGzH+OH//mBf575T6anTA/63Gwofk4DjSEOV0QuEJFCEXGKSJtzIiJypogUichWEbnN43yKiKx2n39DRHyKSmOWjA979uwxWgJgnI5ekb24MONCXj7vZcr/r5zKA5VU1VXxxbYvuOCtC+h1fy+W/biMQ7WHeH7t86wvXx8UJ3zgwIGAtu9wOnA4HfxlyV+Y8/ockh5J4pcf/BKAmIgYLj3xUjZfv5lp/aYB0NPW09AbYKH+OQ0ERk0pbATOA55qq4KIWIHHgZ8BpcAaEVmklNoEPAA8opR6XUQWAPOBJwMvW+NvGm+kjU8czyuJrwBwxH6ECGsEe4/u5cuSL3lg+QPsPLKTR2Y+wjVjr+GNjW8wst9IUvqkEBsR6xenZLFY/HZjpt5Rj4hQUVPBU7lPkVeWR15ZHjeOu5E/nPoHFIpLT7yUh894mGF9hgEwb8y8puu/53u/6NCYD0McrlJqM9DeF2UcsFUptc1d93VgtohsBqYDl7rrLQTuQjvc44aetp4ADOk1hFfO+8kJNzgbaHA2sKh4EXctuYudh3dydurZvHnBm9z6xa1sP7ydQTGDGNxzMDfn3Myeqj0crD3IoNhB9LL1OubnbeHChWzfvv2YupzKiUUsbNy7kS0HtrC7cjcHag5wx5Q7eGfTO/xl6V/YXbmbQ7WH+PyKzxnVfxRVdVVceuKl/GPmP0jp7YoJcMeUO/xjKE23w8w3zQYDnnmrS4HxuG4pHFJKNXicHxxkbZog0+iEgSYnDDSFkbwg4wK2HNhCWVUZh2oPISJ8VfIV9y69l92Vu6l31rN43mL6Rvfl0ncuJcwSRpgljCtOuoJfj/019268l4mOidz/0v3UO+q5cdyNnD/qfMY/O56jdUfZU7WHCYkT+OjSj3ip4CWKDhSREJPAoNhBOJWTnKQcnp/9PAmxCfSP7t80cn/gZw8E11AaUxOwdbgi8iUwsJWiPymlPnDXWQz8n1KqxeJYEZkLnKmUutr9/ApcDvcuYJVSaoT7fBLwqVJqdBs6rgGucT8djWs6w2j6AfuNFoHW0Rytwxutw5t0pVRsVxoI2AhXKTWji03sAjzX6SS6zx0AeotImHuU23i+LR1PA08DiEhuVxcu+wOtQ+vQOrqnjq62YeZlYWuAVPeKhAjgYmCRcg3JvwHmuuvNAz4wSKNGo9H4jFHLws4VkVIgB/hYRD5znx8kIp8AuEevNwCfAZuBN5VSjVtObgV+LyJbcc3pGpOMS6PRaDqAUasU3gPea+X8buBsj+efAJ+0Um8brlUMHeXpTlwTCLQOb7QOb7QOb44bHSEVvEaj0WiMxMxzuBqNRnNccdw5XKO2DbfSfpyIfCEiW9x/W+TgFpFpIrLO46gVkTnushdEpMSjbEygdLjrOTz6WuRxPpj2GCMiK93v33oRucijrNP2aOu99ii3uV/bVvdrTfYou919vkhEZnb8lXdIx+9FZJP7tX8lIkM9ylp9fwKk4yoR2efR39UeZfPc7+EWEZnX/Fo/63jEQ0OxiBzyKPOnPZ4Tkb0i0uqSUXHxb7fO9SKS5VHWMXsopY6rAzgBSAcWA9lt1LECPwDDgAigABjlLnsTuNj9eAHw207q+Dtwm/vxbcAD7dSPAyqAaPfzF4C5frCHTzqAqjbOB80eQBqQ6n48CCgDenfFHsd6rz3qXAcscD++GHjD/XiUu74NSHG3Y+3k6/dFxzSP9/+3jTqO9f4ESMdVwGNtfEa3uf/2cT/uEygdzerfCDznb3u425oMZAEb2yg/G/gUEGACsLqz9jjuRrhKqc1KqaJ2qjVtG1ZK1QGN24YF17bht931FgJzOilltvt6X9uZi2sDR8uAsV2jozqaCLY9lFLFSqkt7se7gb1A/07210ir7/UxtL0NnO5+7bOB15VSdqVUCbCVzt2s9UmHUuobj/d/Fa415v7GF3u0xUzgC6VUhVLqIPAFcGaQdFwCBCQ5m1JqKa7BTlvMBl5ULlbh2geQQCfscdw5XB9pbdvwYPy7bTheKdUY728PEN9O/Ytp+YG6z/0T5hERsQVYR6SI5IrIqsZpDQy0h4iMwzXy+cHjdGfs0dZ73Wod92s9jOu1+3Ktr3S0rfm4RlWNtPb+BFLH+W5bvy2u3ZwdudafOnBPraQAX3uc9pc9fKEtrR22h5ljKbSJ+LBt2Ggdnk+UUkpE2lwO4v5veSKuNceN3I7LMUXgWo5yK3BPAHUMVUrtEpFhwNcisgGX4/EZP9vjJWCeUsrpPu2zPbo7InI5kA1M8Tjd4v1RSv3Qegtd5kPgNaWUXUR+g2v0Pz1AffnCxcDbSimHx7lg2sNvdEuHq8yzbbhNHSJSLiIJSqkytwPZeww9FwLvKaWaAvZ6jAbtIvI88H+B1KGU2uX+u01cMS5OBt4hyPYQkZ7Ax7j+ea7yaNtnezSjrfe6tTqlIhIG9ML1WfDlWl/xqS0RmYHrH9QUpVRTxPw23p/OOJh2dSilPAMDP4tr/r3x2qnNrl3cCQ0+6fDgYuD6Zhr9ZQ9faEtrh+0RqlMKwdg2vMh9vS/ttJifcjulxnnUOXQ+6E67OkSkT+NPdBHpB5wKbAq2PdzvxXu45sveblbWWXu0+l4fQ9tc4Gv3a18EXCyuVQwpQCrwnY/9dliHiJyMK0b0LKXUXo/zrb4/AdSR4PF0Fq6dnuD6BXaGW08f4Ay8f5X5VYdby0hcN6RWepzzpz18YRFwpXu1wgTgsHsA0HF7+OtOn1kO4Fxccyl2oBz4zH1+EPCJR72zgWJc/xX/5HF+GK4v1VbgLcDWSR19ga+ALcCXQJz7fDbwrEe9ZFz/KS3Nrv8a2IDLsbwMxARKBzDR3VeB++98I+wBXA7UA+s8jjFdtUdr7zWu6YhZ7seR7te21f1ah3lc+yf3dUXAWV38bLan40v3Z7bxtS9q7/0JkI6/AYXu/r4BRnpc+yu3nbYCvwykDvfzu4D7m13nb3u8hmtFTD0u3zEfuBa41l0uuJIh/ODuL9vj2g7ZQ+8002g0miARqlMKGo1GE3S0w9VoNJogoR2uRqPRBAntcDUajSZIaIer0Wg0QUI7XE23Q36KFLVRRD4Ukd6dbOcqEXnMz/Ia2/5jINrVdG+0w9V0R2qUUmOUK1NzBc12IRmJe3G8BdAOV9MC7XA13Z2VuAOGiMhwEfmviOSJyLfuXUqIyDniinO7VkS+FJH2AufcJSIviSs27xYR+bX7fIy44tTmi8gGEZntPp8srriuL+LamPEfIMo9Cn8lkC9e073olrEUNBoAEbECp/NTEtGnce0O2iIi44EncAVdWQZMUEopcQXT/gNwczvNn4Qr9mkPYK2IfIwr/sO5Sqkj7i2lq+Sn4NepuILtrHJru0ApNcZfr1VzfKAdrqY7EiUi63CNbDcDX4hIDK4tn2+5wi0ArsDh4Aoq8oY7RkAEUOJDHx8opWqAGhH5Blf81o+Bv4rIZMDp7r9xtLxDeQTb0WhaQ08paLojNe7R41Bc+9yvx/VZPuSe2208TnDXfxRXBoMTgd/gip3QHs33vCvgMlwB0ce6+y/3aOtoF16PJkTQDlfTbVGu7Aj/g2t6oBooEZELoOnmVaa7ai9+Cv3nax6u2SISKSJ9cYXgW+NuZ69Sql5EpuFy+G1RLyLhHXpBmuMe7XA13Rql1FpgPa4Ql5cB80WkAFe0q8aULXfhmmrIA/b72PR6XJGyVgF/Ua6UP68A2eIKzn4l8P0xrn8aWK9vmmk80dHCNJpmiMhduJIUPmS0Fs3xhR7hajQaTZDQI1yNRqMJEnqEq9FoNEFCO1yNRqMJEtrhajQaTZDQDlej0WiChHa4Go1GEyS0w9VoNJog8f8BmHMmp6cnR/QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pydmd import DMD\n",
    "dmd = DMD(svd_rank=-1)\n",
    "dmd.fit(snapshots)\n",
    "dmd.plot_eigs(show_axes=True, show_unit_circle=True, figsize=(5, 5), filename='eigs_ex1.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see the eigenvalues are all inside the unit sphere, therefore the corresponding dynamics are all stable. \n",
    "\n",
    "By contrast, if we use DMD with control we see that there are some unstable dynamics in the \"real\" `A` operator, which we lose when we use the original algorithm, but are recovered if we employ the modified one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAE9CAYAAACV7tCLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9hklEQVR4nO29e1yc1Z34//4ACRAhJJMLSQDDoAz5hhgQUMSSm7rrXb+t2m2y2rRr197WanXbrc1rf267Tb9tta1tbaut7larVavbq6ZatblIi2hgQ8zEhMQMSSYXciEXEJiE4fz+GIIZhkmGh+F5TvKc9+s1L8JznjnnnfPMfDjP5XyOKKUwGAwGt5LitIDBYDA4iQmCBoPB1ZggaDAYXI0JggaDwdWYIGgwGFyNCYIGg8HVpDktMBwmT56sCgsLk1LXsWPHGDt2bFLqSjbGzRrGzRq6uiXbq7Gx8YBSasrg7WdUECwsLGTt2rVJqau7u5vMzMyk1JVsjJs1jJs1dHVLtpeIbB9qu2tPh4PBoNMKcTFu1jBu1tDVzS4v1wbB3NxcpxXiYtysYdysoaubXV6uDYLt7e1OK8TFuFnDuFlDVze7vM6oa4LJpLW1lWTdZEk2xs0aurh1d3ezbds2wuHwwDalFOvXr3fQKj66uln1Sk1NpaioKOHria4NggbDaLFt2zYmT57MlClTSEmJnGyFw2FSU1MdNhsaXd2sePX19bF//362bdtGaWlpQu9x7emwwTBahMPhqABosI+UlBSmTJkSNQo/7XtG0cdgcC0mADrHcPvetUdq2rRpTivExbhZQ2c3EbG1vdTUVGbNmjXw+upXvwrAP/zDP9DU1OSoG8C4ceNOu49dXq69Juj1ep1WiItxs4bObnYHmvT0dDZt2hSz/bnnnovZ5kQQTAS7vFw7EmxpaXFaIS7GzRo6u/X19cUte/RRmDEDUlIiPx99dPQ8Lr74Yt544w0AHnroIQoLC5k7dy4f+9jHWLp0KQC7d+/myiuvZM6cOcyZM4dXX30VgHvvvZePfvSjXHzxxeTn57N8+XIAPv/5z/Otb31roI17772X+++/nyNHjlBTU8Ps2bPx+Xz86le/ivFZsWIFixYtGvh96dKl/OhHPwKgrq6Oiy66iNLSUubNm8f27ZEJH8uXL+e8887D5/Nx3XXXjbxTlFJnzKuyslIli3A4nLS6ko1xs4Yubs3NzTHb+vr6htz3kUeUyshQCj54ZWREto+ElJQUVVJSMvB67LHHlFJKXXTRRWrNmjUqEAioGTNmqLa2NtXT06MqKyvVxz/+caWUUtdff7165ZVXlFJKtbS0KK/Xq5RS6p577lHl5eWqq6tL7d69W+Xk5Kienh7117/+VVVVVQ20XVRUpLZu3aqOHTum2tvblVJK7d69WxUUFAwco8zMTKWUUi+99JJauHDhwHs//vGPqx/+8Ieqp6dHlZeXq127dimllHrsscfULbfcopRSasqUKaqrq0sppdT+/fuH/P8PdQyAtWqIuOLa0+E1a9awcOFCpzWGxLhZQ2e3vr6+IR/3+NrXoKcneltPT2T7pz9tvb14p8MnqKur45JLLmHq1KmEw2E+/OEPD4yk6+rqokbVnZ2dHDlyBIArr7ySzMxMMjMz8Xg87Nq1i0svvZSDBw/S2trK3r17ycnJ4bzzziMUCnH33XdTX19PSkoK+/btY9euXRQUFJzW/5133mHLli1cdtllQKT/pk6dCkBJSQkf/vCHufHGG1myZInlPjqBa4OgwaADe/cOb7sdKKVoamoa8uZFenr6wL9TU1M5fvw4ADfccANPP/00e/bs4aabbgLgZz/7GQcOHOCdd94hPT2dvLw8urq6oupLS0uLulTQ0/8Xoa+vj/PPP59169bFOKxcuZKXX36Z3//+9zzwwANs3ryZMWPGWP7/uvaaoMGgA/FuaI/2je4PfehDvPnmm+zfv5/jx4/zu9/9bqBs3rx5Udf46uvrT1vfbbfdxgsvvMCLL77IrbfeCsCRI0eYMmUK6enpvPjii+zevTvmfeeddx5bt26lu7ubAwcOUFdXB8DcuXNpb2/n9ddfByAUCtHY2Eg4HOa9997juuuu4+GHH44apVrF0ZGgiEwAHgPmAAr4J6XU6XvcYDhLuP9+uPvu6FPijIzI9pEQCoWYNWvWwO+XX345P/7xjwd+93q93HPPPVRVVZGTk0NxcTE5OTkAPProo3zqU5/C5/MRDoe55JJLqKmpOWV7lZWVvP/+++Tm5jJz5kwAbr/9dq6++mp8Ph9z584d8u79eeedxw033MCsWbMoKCgYmOWRkZHBs88+yxe+8AU6OjoIh8N87nOfY86cOSxZsoSOjg6UUnzqU59i8uTJI+usoS4U2vUCngA+1f/vscCEU+2fzBsjK1euTFpdyca4WUMXt6Euyvf29sbd/5FHlJo+XSmRyM+R3hRJlMOHDyullOru7laLFi1STz75pD0NJ8ip+ux0DOfGiCiHFl8XkRxgHVCkEpSoqqpSyUqqGgqFoq5v6IRxs4YubuvXr2fu3LlR25RS2j2P9+lPf5rVq1cTCoVYuHAhjz/+uFYzXUbSZ0MdAxFpVEpVDd7XydNhL7Af+G8RKQMagbuUUu+fvJOI3AHcATBjxgxWrVoVVcm0adPwer20tLRQWlrKmjVrYhqqqakhEAiQn59PW1sbO3fupKenh4yMDADy8vLIz88nEAjg8/kGrkucTG1tLS0tLXi9XoLBILt27YoqLygoIDc3l2AwiNfrHfI6yvz58/H7/fh8PgKBAHsHXf0uLCzE4/HwzjvvUFFRQUNDQ1S5iLBgwQKam5spLS2lpaWFffv2Re1TVFREdnY27e3t5ObmxmTiTktLo7a2lqamJsrLy/H7/Rw8eDBqn+LiYtLT0+ns7MTj8UTNMOjp6SEnJ4eamhoaGxuprKykubmZQ4cORdVRUlJCamoqPT09ZGdn09zcHFWemZlJdXX1QB1NTU0cPXo0ap/Zs2cTDofp6+sjPT2dDRs2RJVnZWVRVVU1UMdf//rXmC/xnDlzCIVCpKSkkJqaysaNG6PKx48fT0VFxUAdDQ0NdHd3R+1TVlZGR0cHGRkZhMNhNm/eHFU+ceJEysrKBupQSg05d/XkL/RQf/dTUlIG7iIP9f6UlJSBwKA+OJuKql9EEq7jpz/96Qejof46+/r6Buo44WPV48Q+Vus4UX6qOk54Dq5DKUVraysej4e2tjby8/Nj3j/QloMjwSrgTeBDSqkGEfkBcFQp9e/x3pPMkWBnZydZWVlJqSvZGDdr6OJ2powET6Crm10jQSfHvkEgqJQ6Mdx5Aaiwq/G2tja7mho2xs0aOrs5NdhIBF3d7PJyLAgqpfYCO0WkpH/T5cDGU7wlqezcudOupoaNcbOGzm66BhrQ180uL6cflr4TeFpExgLbgE867GMwGFyGo7eClFLrlFJVSqm5Sqn/q5Q6dPp3GQyGU7F582aKi4ujtp1IanAq3njjDT75ycg4ZMWKFbz22mvDanfBggUcOHAg4f1vvvlmfvGLXwyrjdHA6ZGgweAKbv3trTzrf3bg98BdAeqD9Sz5nw/mvj74dw9y76X3Il/74GbAosJF/GXpX7jsictY2bpyYLu6P/mnivPmzWPevHkAvP7662RlZXHFFVck/P7Vq1fHbOvr60MppWX6/hO4Ngjm5eU5rRAX42YNXd1Kf1LKO595h2dufiZqe+GEQhbPWRyz/1AB7i9L/5I0n4svvpiqqirq6uo4evQojzzyCFdddRUrVqzggQce4JFHHuGJJ54gJSWF5557jh/84AdceeWVA+8/cuQIt99++8AiSMuWLWPp0qXk5eWxdu1ajh49ylVXXcWFF17I+vXrWbFiBU888QTPP/88IsIVV1wRNXsFIkkbvvjFL9LV1YXH4+Gpp55KKNFCMnBtEDzVc0NOY9ysoavbxv0bk/oIynf/9l3uvfTeEdXR29vL+vXref755/n617/OVVddNVBWUlLC0qVLycrK4mtf+1rMe++77z5ycnIGMs3s378/Zp/t27fz+OOPc9lll/HCCy+wYsUKGhsbyc7Ojnm2NRQKceedd/LSSy8xY8YMHn/8cb70pS8NmQB2NNDn8XCbCQQCTivExbhZQ2e3ZN7p/NdX//WU5fEC7snbb7nlFiAykSAYDA6r/dWrV/PFL35x4PcpU6bE7DN9+vSBNFh//vOfue2228jOzgYYSIl1gpPTZs2aNYtvf/vb7N692zV3hx3D5/M5rRAX4zZ8duzYwbnnnuu0xpBMGTfF1oeRc3NzYzKrtLe3RyUwODFbKt7MkpGSyBoiJ4iXNuusf07QaYaaGqcLxm14KKW48MILufbaa51WGZJ9X9p3yvT6ySYnJ4epU6fyxz/+MdL+vn2sXLlyYGSWCNnZ2THTGE+wcOFCvv/97w/8PtTp8MlceeWV/PKXv6Sjo2PA52Tipc2yq89cGwQNZw9/+MMfOH78OG+//Tbt7e1O68Rw7ysju343mF/dFLtWx2CefPJJvvGNbzBr1izmz5/PV77yFWbPnp1wGzfddBMvvfQSs2bN4pVXXokq++Y3v8nhw4cpLi6mpKSEP/3pT6et6+qrr6a8vJxZs2bx9a9/Pao8IyOD5557jn/7t3+jpKSE0tLSIe80jxaOzR22QjLnDq9atUrbVOzGLXGUUlRUVAyMXG644YaYL5ndDJ63evIjLwBvfCKy0NG8X8wb2HbPJffw3Su/y9QHprK/KzKymj1lNv7P+Vn8P4t5dkP04zWFEwqT5hsOh7V8hGUkXmdKFhmDYcT84Q9/ACKPx1xwwQX85Cc/4e6778bj8Ths9gHqfjXkF3qoR2H2fWlfzLZnbnqGZ256Jma7ITmY02HDGc3u3bv53ve+h4gwceJE7r777mHf7TS4G9eOBGtra51WiItxS5zPfvazAPz85z+npKSE2267zWGjodEpWelgdHWzy0vP/70N6LxQt3Gzxp49e5xWGGDwnU2dr73r6mbVa7h3lV0bBIda9EUXjJs1Bj+E6xSpqans378/6suoY9LSE+jqZsWrr6+P/fv3D+uGimtPh4PBYEymDV0wbtbQ5fGYoqIitm3bFpXkVdfszaCvm1Wv1NRUioqKEt7ftUFw165d2n6ZjZs1dAmCmZmZA0tHnkC3R4tORlc3u7xcezpsMBgMYIKgwWBwOSYIGgwGV+PaIGhXwkYrGDdrTJo0yWmFuOjcb7q62eXl2iCYm5vrtEJcjJs1cnJynFaIi879pqubXV6uDYI6T60ybtbQ5e7wUOjcb7q62eXl2iCo80O/xs0aQ2U41gWd+01XN7u8XBsE6+vrnVaIi3GzxpYtW5xWiIvO/aarm11erg2CBoPBACYIGgwGl2OCoMFgcDUmCBoMBlfj2iA4f/58pxXiYtysMWvWLKcV4qJzv+nqZpeXa4Og3+93WiEuxs0auj7vBnr3m65udnk5HgRFJFVE/ldEXrSzXV0XEQfjZpXp06c7rRAXnftNVze7vBwPgsBdwLt2NxoIBOxuMmGMmzVOtwi4k+jcb7q62eXlaBAUkXzgWuAxu9veu3ev3U0mjHGzxuHDh51WiIvO/aarm11eTo8EHwK+DAxvZRSDwWBIEo6l1xeR64B9SqlGEVl4iv3uAO4AmDFjBqtWrYoqnzZtGl6vl5aWFkpLS1mzZk1MHTU1NQQCAfLz82lra2Pnzp10dnYO1JWXl0d+fj6BQACfz0ddXV1MHbW1tbS0tOD1egkGg+zatSuqvKCggNzcXILBIF6vd8gpP/Pnz8fv9+Pz+QgEAjF/6QoLC/F4PIRCIbq7u2loaBjcFyxYsIDm5mZKS0tpaWlh377oxbqLiorIzs6mvb2d3Nxc1q5dG1WelpZGbW0tTU1NlJeX4/f7OXjwYNQ+xcXFpKen09nZicfjoampaaCss7OT+vp6ampqaGxspLKykubmZg4dOhRVR0lJCampqfT09JCdnU1zc3NUeWZmJtXV1QN1NDU1cfTo0ah9Zs+eTTgcpq+vj/T0dDZs2BBVnpWVRVVVFY2NjQAcO3Ys5vMxZ84cQqEQKSkppKamsnHjxqjy8ePHU1FRMeDR0NBAd3d31D5lZWV0dHSQkZFBOBxm8+bNUeUTJ06krKxsoI76+npCoVDUPuFwmNbWVrKysgiFQjFT/CZNmkRpaSnr1q2joqKCuro6ent7o/apqqqira0Nj8dDR0cH27ZtiyqfOnUqPp8Pv99PWVkZq1evjlmxrbq6mmAwSG5uLu3t7bS2tkZ9F6x+n04mWd+ncDjMpk2bRvx9amtrIz8/P+b9AyilHHkB/w8IAq3AXqALeOpU76msrFTJYuXKlUmrK9kYt+GzePFitWzZMqc14qJrvymlr1uyvYC1aoi44tjpsFLqPqVUvlKqEPgY8Bel1K12tV9YWGhXU8PGuFlD5ywyOvebrm52eTl9TdAxPB6P0wpxMW7WyMrKclohLjr3m65udnlpEQSVUquUUtfZ2ebJa8LqhnGzxpEjR5xWiIvO/aarm11eWgRBJzjlhVKHMW7W0HVEA3r3m65udnm5NggOvvOqE8bNGlu3bnVaIS4695uubnZ5uTYIGgwGA5ggaDAYXI4JggaDwdW4NgiKiNMKcTFu1jBu1tDVzS4vUYOm1ehMVVWVGjwNzGAAWLJkCddddx1LlixxWsWgKSLSqJSqGrzdtSPBwXNZdcK4WWP79u1OK8RF537T1c0uL9cGwdLSUqcV4mLcrKHr826gd7/p6maXl2uDYEtLi9MKcTFu1tizZ4/TCnHRud90dbPLy7VBcHAKKp0wbtYYnIpLJ3TuN13d7PJybRA0GAwGMEHQYDC4HBMEDQaDq3FtECwqKnJaIS7GzRpTp051WiEuOvebrm52ebk2CGZnZzutEBfjZo3MzEynFeKic7/p6maXl2uDYHt7u9MKcTFu1ujs7HRaIS4695uubnZ5uTYI5ubmOq0QF+NmjZycHKcV4qJzv+nqZpeXa4OgznOQjZs1Bi9BqRM695uubnZ5uTYIGgwGA5ggaDAYXI4JggaDwdW4NgimpaU5rRAX42aN1NRUpxXionO/6epml5dJqmo4KzBJVQ2nwyRVHURTU5PTCnExbtYIBAJOK8RF537T1c0uL9cGwfLycqcV4mLcrDFz5kynFeKic7/p6maXl2uDoN/vd1ohLsbNGsFg0GmFuOjcb7q62eXl2iB48OBBpxXiYtysofO0OZ37TVc3u7xcGwQNBoMBTBA0GAwux7EgKCIFIrJSRDaKiF9E7nLKxWAwuBcnn5LsBe5VSjWJSDbQKCKvKqU22tF4cXGxHc1YwrhZY9q0aU4rxEXnftPVzS4vx0aCSqk9Sqmm/n93AO8CeXa1n56ebldTw8a4WWPMmDFOK8RF537T1c0uLy2uCYpIIXAh0GBXmzrfSTRu1ujp6XFaIS4695uubnZ5OT5pUESygP8B7lZKxSwcKyJ3AHcAzJgxg1WrVkWVT5s2Da/XS0tLC6WlpaxZsyamjZqaGgKBAPn5+bS1tbFz5076+vpobW0FIC8vj/z8fAKBAD6fj7q6upg6amtraWlpwev1EgwG2bVrV1R5QUEBubm5BINBvF4v9fX1MXXMnz8fv9+Pz+cjEAiwd+/eqPLCwkI8Hg8dHR10d3fT0BD9N0FEWLBgAc3NzZSWltLS0hKzNmtRURHZ2dm0t7eTm5sbk5MtLS2N2tpampqaKC8vx+/3xzyKUFxcTHp6Op2dnXg8nqgn9/v6+tizZw81NTU0NjZSWVlJc3Mzhw4diqqjpKSE1NRUenp6yM7Oprm5Oao8MzOT6urqgTqamppi1g2ePXs24XCYvr4+0tPT2bBhQ1R5VlYWVVVVNDY2AnD48OGYz8ecOXMIhUKkpKSQmprKxo3RV1vGjx9PRUXFgEdDQwPd3d1R+5SVldHR0UFGRgbhcJjNmzdHlU+cOJGysrKBOurr6wmFQlH7+Hw+WltbycrKIhQKsWXLlqjySZMmUVpayrp166ioqKCuro7e3t6ofaqqqmhraxv4jAzOnzh16lR8Ph9+v5+ysjJWr17N4Gmx1dXVBINBcnNzaW9vp7W1Neq7YPX7dDLJ+j6dc845bNq0acTfp7a2NvLz82PeP4BSyrEXMAZ4Bbgnkf0rKytVsli5cmXS6ko2xm34LF68WC1btsxpjbjo2m9K6euWbC9grRoirjh5d1iAx4F3lVLfc8rDYDC4m9MGQRH5diLbLPAh4DbgMhFZ1/+6Jgn1GgwGQ8IkMhL8uyG2XT3ShpVSdUopUUrNVUqV979WjLReg8FgGA5xb4yIyGeBzwHnicj6k4qygb+Otthoo+tjAWDcrGIekbGGrm52ecVNqioiOcBE4P8BXzmpqEMp5chCpSapqiEeJqmq4XQMO6mqUuoIsBO4UCm1/aSXnis1D5MTj1XoiHGzhs5JVXXuN13d7PI65TVBpVQY2Cwi59piYyOVlZVOK8TFuFnD6/U6rRAXnftNVze7vBK5MTIR8IvI6yLyhxOv0RYbbQY/vKsTxs0a27dvd1ohLjr3m65udnklMmPk30fdwgEGz3DQCeNmjffff99phbjo3G+6utnlddogqJRabYeIwWAwOEEiD0tfIiJvi0iniBwTkbCIxMzxNRgMhjORRK4JPgwsBrYAmcCngB+PppTBYDDYRUJzh5VSW4FUpVRYKfXfwFWjqzX6lJSUOK0QF+NmjenTpzutEBed+01XN7u8Erkx0iUiY4F1IvIdYA+a5CEcCampqU4rxMW4WSMlRd+Ppc79pqubXV6JfGpu69/vX4D3gQLgptGUsgOdE3AaN2scP37caYW46NxvurrZ5ZXI3eHt/SPBQuA3wGal1LHRFhttsrOznVaIi3GzRmZmptMKcdG533R1s8srkbvD1wLvAT8kcpNkq4iMOIuM0+j6gCgYNysolHlY2iK6utnllcjp8HeBRUqphUqpBcAi4Pujq2UwDE338W5e2foKy9cs5yPPfYTG3Y20HGzh2XeeJXA4wLjl41i+ZjkAz254lrd3vR2TYt5gOJlEgmBH/93hE2wDOkbJx2AYkvfa3+NY+BgvbXmJ5W8s53DPYT5a+lEKJxRS7Clm8QWL8U7wsu9L+/jsRZ8FoG5HHbf+9lZKf1LKr975lcP/A4OuJHJ3eK2IrAB+DSjgFuBtEfkIgFLqN6PoZ3A5O47s4M4/3Un9znpeve1Vbp59MzfPvjnu/lljswb+/fA1D6OUYvX21RwNHUUpxWNNj3Fb2W1kpGXYoW84A0hkJJgBtAELgIXAfiIPTV8PXDdqZqOMzhfRjVuE4NEglT+rpDqvmh1f3EHZtLJT7j927NiYbSLCwsKF3FByA13Hu3j5vZep/Fklb+16a7S0h8Qc0+Fjl1fcpKo6YpKquoPDPYdZ3bqaG2fdSPBokPzxp1gusZ9Ek6oqpXjO/xx3vXwXK5asoHKGnmmkDMln2ElVz3Z0TSQJ7nbb9/4+Fj2xiLodkbVqEwmAJ0gkqaqI8LE5H2Pj5zZSMb2C9W3r6VN9ln0Txc3H1CpaJFU9m9E1kSS41+1wz2H+/pd/zzXnX8N3/u47w37/cJKqTho3CYB7XrmHu1++e9TvILv1mI4EbZKqioiec2pGSFNTk9MKcXGrW29fL58o/wTfuOwbRJalHh7DTa8vIrzw0Rd4Y8cb/LDhh8Nubzi49ZiOBLu8EhkJbhGRB0Rk9qjb2MjRo/pmA3Oj24otKzgWPsbdl9xtKQACdHd3D/s9EzIm8Pwtz/Prjb8m3Be21G4iuPGYjhS7vBIJgmVAC/CYiLwpIneIyPhR9jK4iMChAEt/t5SjIWe+jOd7zqfuk3Uc7zs+qoHQoCenDYJKqQ6l1M+VUpcC/wbcD+wRkSdE5PxRNzSc9Xxuxef48qVfZtbkWY45iAj/suJfeLTxUcccDM6Q0DVBEblBRH4LPERkGl0R8EdgxejqGc52DnQdoPVwK3dfcrfTKtx+4e089OZDttwtNuhDQtcEgRuBB5RSFyqlvqeUalNKvQC8PLp6o8fs2fpe4nST2+Rxk/F/zs+Y1DEjrisvL29E778k/xLGjRnH69teH7HLYNx0TJOFXV6nDIL9d4Z/oZS6XSn1t8HlSqkvjJrZKBMO63vtxy1uncc6Wfw/i5P2eEpf38hGcCLCt6/4NhMzJybF52TcckyTiV1eiSy+fsZOjTsVI/3CjCZucXv1vVdp724nNSU5T2ElI5heef6VzJ4ym96+3iQYfYBbjmkyscsrkdPhv4rIwyIyT0QqTrxG3WyUSU9Pd1ohLm5xW7t7LZfkXZK0+saMGfkpNcC8/55H4+7kzlZwyzFNJnZ5JRIEy4FS4OtEbop8F3gwGY2LyFUisllEtorIV5JRZ6Js2LDBzuaGhVvcth/ZntS5uzt37kxKPRdOu5C1u5M7R90txzSZ2OWVyCMyi4Z4XTbShvuvN/4YuBqYDSw+2x7INpyapz7yFNf7rretvaefhsJCSEmJ/Hz66aH3q5pRReOe5I4EX3ttakJtn0309fXxzW9+k5aWFkfaH3y8YbJnqP0SmjssIteKyJdF5P878UqC48XAVqXUtv41S54lchfa4AJ6+3r5z9X/aXl2yHB5+mm44w7Yvh2Uivy8446hg9GlBZdy4bQLk9r2gw+WJNT22YSIkJaWxoc+9CE+/vGP2xoMhzrecO7MofY9bVJVEXkEGEckrf5jwM1AMpKx5QEnn78Egeok1Js0tm3bxre//W3b07Pv3r2bX/1Kz0zIyXI7Hj7OU+ufYueFyTmFfeutt3j33XfZt2/fkOXPPANdXdHburrgM5+B1auHrvOOx+9Iitszz0AoNLy27WS0P2/XXnstfr+fCy64gIqKCurr60etrRMsWxZ7vEGGHPQlkln6UqXUXBFZr5T6moh8F/jTSCUTRUTuAO4AmDFjBqtWrYoqnzZtGl6vl5aWFkpLS1mzZk1MHTU1NQQCAfLz82lra2Pnzp10dXUN1JWXl0d+fj6BQACfz0ddXSSN05EjR8jOzkYpxfnnn09bWxuTJ0/m0KFDHD58OKoNj8dDdnY2hw8fZtKkSWzbti3Go7i4mD179jB16lQOHjzIkSNHosonTZrEOeecQ3Z2Nl6vNyYhgIjg8/nYuXMnM2bMoK2tjY6O6JUOpkyZQnp6Ol1dXWRnZ8csPpSamsr555/Pjh07yM/PZ8+ePXR2dkbtM3XqVMaMGUNPTw/nnHMOO3bsGCjzer2MGzeO8847j+3btzNz5syB/hx8XESE3t5e0tPTCQaD0R5pqaTkpTBlyhRmzpzJjh07Yub+Tp8+HaUUfX19jBkzhl27dkWVp6enU1hYyPbt2/F4PAPH4GTy8vI4fvw4nZ0TgNhRZ2enoqpKBv4vgUCAzu5OOo914smMnD3l5+cTCoVIS0tDKcXevXuj6hg3bhwFBQUDdbz33nv09vae1IYvbtvZ2ZHRUVZWFtOnTycYDHLuueeydevWmEdEZs6cSUdHB+PGjSMUCrF///6o8uzsbHJzc9m9ezcFBQW0tLTE/AH3er0cOnSI8ePH8/7773Pw4EG8Xu9AQtqcnBwmTZrEvn37mD59Olu2bInxLioq4uDBg0yYMIGOjg7a29ujyidMmMDEiRM5cOAAubm5/O1vf+PQoUOMGTOG6dOns2rVKmpra2lpacHr9RIMBmOObUFBARkZGWzatAmv1ztk4Jw/fz5+vx+fz0cgEIg6Ljt2LBiyz4dEKXXKF9DQ//NNYAaQTuQ09rTvPU29NcArJ/1+H3Dfqd5TWVmpDGcHveFelfK1FBXuCyelvsWLF6unn346bvnMmUpFToyiXzNnxu77x81/VFc/dXVSvIbb9tlEOBxWH/3oR9X06dPVQw89pLq6umxre+g+r1RqiLiSyDXBF0VkAvAA0AS0As8kFmJPydtAsYh4+9c1/hjwhyTUmxC6JpIEd7ilpqSy9p+Tewf2VKm0li+HceOit40bF9k+mJaDLZw38bykeS1fDhkZ0aO6eG07wWh93lJSUvjkJz/Je++9x1133TXsdPkj8RrqeEOc+ZBDRcZ4LyKjwJzhvOc09V1DJEPNe8Cy0+1vRoJnF62HWlXgUCApdZ1uJKiUUk89FRkhiER+PvXU0Pvd99p96qnmOIUWSbRtQ/IY3OcweZsaIq4kck0QEbkUKKT/GqKIoJR6cvjxOSYAr8ChJAxr166lqipmuQEtcIvbU+uf4lDPIR78+6Q8djrkddiT+cd/jLxOxzcv/2ZSfE6mpGQtra1n/zFNJiP1Gny8RQ60D7VfIllkfknk4eha4KL+l349NkwG3wzQCbe4XZx38cBaIsmgp6dnxHXs7dzLV1//ahJsonHLMU0mdnklMhKsAmYrdQYtS2c4I1jkXUTw90HWt61nbu5cp3UAeLzpcfa/v//0OxrOGhK5MbIBmDbaIgb3kZaSxrM3P0te9shSYCWL3r5eHm18lM9f/HmnVQw2kshIcDKwUUTeAgYe+VRK3TBqVgbXUHtuLW/vepu0lDRyMnIcddmwbwO+ST7Kp5U76mGwl0SC4H+MtoQTzJkzx2mFuLjN7bGmx0iRFH563U9HVE9BQYHl9/apPsqnlfPqba+OyCEebjumycAur9MGQaWUBhN7kk9o8DwmjXCb23f+7jtc8NMLeG3ba1xRdIXleo4fP275vfe8cg+lU0r558p/tlzHqXDbMU0GdnnFvSYoInX9PztE5OhJrw4R0XONvmGQkqLvuvNuc8vJyOHn1/+c5W+M7Olhq8kY3tj+Br/2/5qP/J+PjKj9U+G2Y5oM7PKKOxJUStX2/8yOt8+ZTGqqvmvKu9HtyvOvZJF3EZ3HOhmTMob0tOEn1LTypdl5ZCf/8MI/8PPrf86kcZOG/f5EceMxHSl2eSXynKBniFdyUvg6yMaNG51WiItb3camjuWBvz7ALc/fQk/v8J/5GzwJ/3SE+8LMyJ7Bczc/x7W+a4fd3nBw6zEdCXZ5JfKnswnYT2R625b+f7eKSJOIJC8tsMEALJu/jHFjxnHN09fQ1tk2au38beffKH+0nPePv8+8mfNGrR2D/iQSBF8FrlFKTVZKTSKSCfpF4HPAT0ZTzuA+xqaO5emPPE1Nfg0vb3355HnmSaH7eDf/+ud/5aZf38Tyy5YzPn180uo2nJkkEgQvUUq9cuIXpdSfgRql1JtEEioYDEklNSWV5ZcvZ2n5Un676bfc/PzNbDqwaUR19qk+jvQcYX/Xfg50HWD9Z9ZzQ4l51NWQWBDcIyL/JiIz+19fBtr61wjRc62+BBg/Xt8RgHH7gGuKr6E8t5xFTyzi8icv53DP4bj7DpWq6WDXQR7824P4fuTjgb89wLk55/KL//sLppwzZRStYzHHdPjY5SWnO9UQkcnA/UQSKAD8FfgacAQ4Vym1dVQNT6KqqkqtXZvcHHSGM4Nj4WO8vPVlrvddz8NvPcxvNv2GyumVVE6v5IaSG/jEbZ9g7ry5zL9uPu8eeJdiTzGXeS9j5kMzWVi4kM9f9HkuzrvYtjVNDPohIo1KqZjkL4msNndAKXWnUurC/te/KKX2K6WO2RkAk40bEpeOBk65jU0dyw0lNyAi3Dr3Vu6rvQ9PpocX3n2BQz2H6DjWwfrN6/n3lf/OW7veQqEQEd77wns8+eEnqc6vdjQAmmM6fOzySmQkOAX4MpG1hzNObFdJWHZzuJiRoCEeS5Ys4brrrmPJkiVOqxg0xfJIEHga2AR4iZwGtxJJjX9G09DQ4LRCXIybNbZu1ffEROd+09XNLq9EguAkpdTjwHGl1Gql1D8Bto8Ck83glc10wrhZ49ixY04rxEXnftPVzS6vRLLInJiVvkdErgV2A0Ou5G4wGAxnGokEwW+ISA5wL/AjYDzwxVG1MhgMBptIJJXWi/3/PAIsGl0dg8FgsJfTBkER8QJ3ctJqc3DmZ5YuKytzWiEuxs0aM2fOdFohLjr3m65udnklcjr8O+Bx4I+cwTNEBtPR0cHEiROd1hgS42YNXS/wg979pqubXV6JBMEepdQPR93EZjIyMk6/k0MYN2uMGaNvhjed+01XN7u8EgmCPxCR+4E/E73QUtOoWdlAOBx2WiEuxs0afX36nqjo3G+6utnllchzghcA/wx8C/hu/+vB0ZSyg82bNzutEBfjZo09e/Y4rRAXnftNVze7vBIZCd4CFCml9H0S1WAwGCyS6OLrE0bZw2AwGBwhkZHgBGCTiLyNWXzdYDCcZSQSBO8fdQsH0PGRgBMYN2ucc845TivERed+09XNLi9HFl8XkQeA64FjwHvAJ5VSh5PdzqnQ9QFRMG5WMQ9LW0NXN7u8TrX4+uBF15O5+PqrwByl1Fwiq9jdN8L6ho2uiSTBuFklEAg4rRAXnftNVze7vE61+PqoLbrev1jTCd4Ebh6ttuJRWanvaqHGzRper9dphbjo3G+6utnllcjd4dHmn4A/2d1ofX293U0mjHGzxpYtW5xWiIvO/aarm11eidwYsYSIvAZMG6JomVLq9/37LAN6iWSvjlfPHcAdADNmzGDVqlVR5dOmTcPr9dLS0kJpaSlr1qyJqaOmpoZAIEB+fj5tbW3s3LmTzs7Ogbry8vLIz88nEAjg8/moq6uLqaO2tpaWlha8Xi/BYJBdu3ZFlRcUFJCbm0swGMTr9Q55AOfPn4/f78fn8xEIBNi7d29UeWFhIR6Ph6NHj9Ld3R2TWVdEWLBgAc3NzZSWltLS0sK+ffui9ikqKiI7O5v29nZyc3MZvBxBWloatbW1NDU1UV5ejt/v5+DBg1H7FBcXk56eTmdnJx6Ph6amDyYHdXZ2Ul9fT01NDY2NjVRWVtLc3MyhQ4ei6igpKSE1NZWenh6ys7Npbm6OKs/MzKS6unqgjqamJo4ejb7KMnv2bMLhMH19faSnp7Nhw4ao8qysLKqqqgZOm95///2Yz8ecOXMIhUKkpKSQmprKxo0bo8rHjx9PRUXFgEdDQ0PMHOSysjI6OjrIyMggHA7HPMQ7ceJEysrKBuqor68nFApF7RMOh2ltbSUrK4tQKBQTsCdNmkRpaSnr1q2joqKCuro6ent7o/apqqqira0Nj8dDR0cH27ZtiyqfOnUqPp8Pv99PWVkZq1evjlmzubq6mmAwSG5uLu3t7bS2tkZ9F6x+n04mWd+nrq4uNm3aNOLvU1tbG/n5+THvH+DE4tZ2v4BPAPXAuETfU1lZqZLFypUrk1ZXsjFuw2fx4sVq2bJlTmvERdd+U0pft2R7AWvVEHFl1EaCp0JEriKyeNMCpVSXEw4Gg8EAzl0TfBjIBl4VkXUi8ohDHgaDweU4MhJUSp3vRLsnU1FR4bRCXIybNXS+O6xzv+nqZpeXDneHHaG9vd1phbgYN2t0dnY6rRAXnftNVze7vFwbBLOyspxWiItxs4auyUFB737T1c0uL9cGwcGPMOiEcbPG8ePHT7+TQ+jcb7q62eXl2iCo84O1xs0ag58T0wmd+01XN7u8XBsEDQaDAUwQNBgMLscEQYPB4GpcGwQnTZrktEJcjJs1dL3LCXr3m65udnm5NgiWlpY6rRAX42aNU06Sdxid+01XN7u8XBsE161b57RCXIybNbZv3+60Qlx07jdd3ezycm0Q1HWqEBg3q5hpc9bQ1c1MmxtlhspxpgvGzRq6LiIOevebrm52ebk2CA5OWKkTxs0a4XDYaYW46NxvurrZ5eXaIGgwGAxggqDBYHA5JggaDAZX49ogWFVV5bRCXIybNYqKipxWiIvO/aarm11erg2CbW1tTivExbhZ48iRI04rxEXnftPVzS4v1wZBj8fjtEJcjJs1dJ42p3O/6epml5drg2BHR4fTCnExbtYYvF6wTujcb7q62eXl2iA4eOFqnTBu1hi8EL1O6NxvurrZ5eXaIGgwGAxggqDBYHA5JggaDAZX49ogOHXqVKcV4mLcrDF+/HinFeKic7/p6maXl2uDoM/nc1ohLsbNGtOnT3daIS4695uubnZ5uTYI+v1+pxXiYtysEQwGnVaIi879pqubXV6uDYJlZWVOK8TFuFlj5syZTivERed+09XNLi/XBsHVq1c7rRAX42aNd99912mFuOjcb7q62eXl2iColHJaIS7GzRrGzRq6utnl5WgQFJF7RUSJyGQnPQwGg3txLAiKSAHw98AOpxwMBoPByZHg94EvA3qOxQ0GgytwJAiKyI3ALqVUsxPtA1RXVzvV9GkxbtY4//zznVaIi879pqubXV5po1WxiLwGTBuiaBnwVSKnwonUcwdwB8CMGTNYtWpVVPm0adPwer20tLRQWlrKmjVrYuqoqakhEAiQn59PW1sbO3fuJBQKkZ6eDkBeXh75+fkEAgF8Pt+QS/3V1tbS0tKC1+slGAyya9euqPKCggJyc3MJBoN4vV7q6+tj6pg/fz5+vx+fz0cgEGDv3r1R5YWFhXg8Hvx+P+Xl5TQ0NAzuCxYsWEBzczOlpaW0tLTEZE4pKioiOzub9vZ2cnNzWbt2bVR5WloatbW1NDU1UV5ejt/v5+DBg1H7FBcXk56eTmdnJx6Ph6ampoGyUCjE+PHjqampobGxkcrKSpqbmzl06FBUHSUlJaSmptLT00N2djbNzdF/7zIzM6murh6oo6mpiaNHj0btM3v2bMLhMH19faSnp7Nhw4ao8qysLKqqqmhsbARg69atMZ+POXPmEAqFSElJITU1lY0bN0aVjx8/noqKigGPhoaGmJRcZWVldHR0kJGRQTgcjlnac+LEiZSVlQ3UUV9fTygUitpn0qRJZGdnk5WVRSgUYsuWLTHlpaWlrFu3joqKCurq6mJWW6uqqqKtrQ2Px0NHR0dMlpWpU6fi8/nw+/2UlZWxevXqmJsL1dXVBINBcnNzaW9vp7W1Neq7YPX7dDLJ+j51d3eTlpY24u9TW1sb+fn5Me8/gdh9Z0hELgBeB7r6N+UDu4GLlVJ7474RqKqqUoO/1FY5evSottOsjNvwWbJkCZdffjm333670ypDomu/gb5uyfYSkUalVEzOfttPh5VS7yilpiqlCpVShUAQqDhdAEw27e3tdjY3LIybNTo7O51WiIvO/aarm11ern1OsLW11WmFuBg3a+zfv99phbjo3G+6utnlNWrXBBOlfzRoMBgMjuDakaDBYDCACYIGg8HluDYITps21NM7emDcrDFhwgSnFeKic7/p6maXl2uDoNfrdVohLsbNGlOmTHFaIS4695uubnZ5uTYItrS0OK0QF+NmjT179jitEBed+01XN7u8XBsES0tLnVaIi3GzxqlmBTiNzv2mq5tdXq4NgkNNB9IF42aNTZs2Oa0QF537TVc3u7xcGwQNBoMBTBA0GAwuxwRBg8HgakwQNBgMrsa1QbCmpsZphbgYN2sUFxc7rRAXnftNVze7vFwbBAOBgNMKcTFu1tA5i4zO/aarm11erg2COj9TZtys4fF4nFaIi879pqubXV6uDYJtbW1OK8TFuFnjyJEjTivERed+09XNLi/XBsHB6yLohHGzxuC1UnRC537T1c0uL9cGQYPBYAATBA0Gg8sxQdBgMLga1wbBvLw8pxXiYtysofPdYZ37TVc3u7xcGwR1fSwAjJtVdA6COvebrm7mEZlRRtcHRMG4WWXfvn1OK8RF537T1c08LD3K+Hw+pxXiYtysMX36dKcV4qJzv+nqZpeXa4NgXV2d0wpxMW7DZ9y4cWzfvt1pjbjo2m+gr5tdXq4Ngoazi5/85CdceumlTmsYzkBMEDScFYwdO5aUFPNxNgwf86kxGAyuxgRBg8HgakQp5bRDwlRVVam1a9cmpa7e3l7S0tKSUleyMW7WMG7W0NUt2V4i0qiUqhq83bUjQV0XnAbjZhXjZg1d3c76xddF5E4R2SQifhH5jt3te71eu5tMGONmDeNmDV3d7PJyJAiKyCLgRqBMKVUKPGi3QzAYtLvJhDFu1jBu1tDVzS4vp0aCnwW+pZQKASilbJ/vtGvXLrubTBjjZg3jZg1d3ezycupqqA+YJyLLgR7gX5VSbw+1o4jcAdwBMGPGDFatWhVVPm3aNLxeLy0tLZSWlrJmzZqYOmpqaggEAuTn59PW1sbOnTvp7OwcqCsvL4/8/HwCgQA+n2/IJ9Vra2tpaWnB6/USDAZjDlBBQQG5ubkEg0G8Xi/19fUxdcyfPx+/34/P5yMQCLB3796o8sLCQjweD6FQiO7ubhoaGgb3BQsWLKC5uZnS0lJaWlpi5ssWFRWRnZ1Ne3s7ubm5DL6RlJaWRm1tLU1NTZSXl+P3+2MyMhcXF5Oenk5nZycej4empqaBss7OTurr66mpqaGxsZHKykqam5s5dOhQVB0lJSWkpqbS09NDdnY2zc3NUeWZmZlUV1cP1NHU1MTRo0ej9pk9ezbhcJi+vj7S09PZsGFDVHlWVhZVVVUDdXR1dcV8PubMmUMoFCIlJYXU1FQ2btwYVT5+/HgqKioG6mhoaKC7uztqn7KyMjo6OsjIyCAcDrN58+ao8okTJ1JWVjZQR319PaFQKGqfcDhMa2srWVlZhEIhtmzZElU+adIkSktLWbduHRUVFdTV1dHb2xu1T1VVFW1tbXg8Hjo6Oti2bVtU+dSpU/H5fPj9fsrKyli9ejWDb3xWV1cTDAbJzc2lvb2d1tbWqO+C1e/TySTr+xQOh9m0adOIv09tbW2nTsaglBqVF/AasGGI1439P38ECHAxEKD/TvWpXpWVlSpZrFy5Mml1JRvjZg3jZg1d3ZLtBaxVQ8SVURsJKqWuiFcmIp8FftMv9paI9AGTAX3XTDQYDGclTl0T/B2wCEBEfMBY4ICdAgUFBXY2NyyMmzWMmzV0dbPLy6lrgv8F/JeIbACOAUv7R4W2kZuba2dzw8K4WcO4WUNXN7u8HBkJKqWOKaVuVUrNUUpVKKX+YreDro8FgHGzinGzhq5uZ/sjMo6j6wOiYNysYtysoavbWf2wtA4MdctdF4ybNYybNXR1s8vLtUHQYDAYwARBg8HgckwQNBgMrsYEQYPB4Gpcm1S1r69P2zUpjJs1jJs1dHVLtpdJqjoIv9/vtEJcjJs1jJs1dHWzy8u1I8FQKER6enpS6ko2xs0axs0aurol28uMBAcRCAScVoiLcbOGcbOGrm52ebk2CA7OPaYTxs0axs0aurrZ5eXaIGgwGAxggqDBYHA5JggaDAZXc0bdHRaR/cD2JFU3GZsTuQ4D42YN42YNXd2S7TVTKTVl8MYzKggmExFZO9Ttch0wbtYwbtbQ1c0uL3M6bDAYXI0JggaDwdW4OQj+zGmBU2DcrGHcrKGrmy1err0maDAYDODukaDBYDCYICgid4rIJhHxi8h3nPYZjIjcKyJKRCY77XICEXmgv8/Wi8hvRWSCwz5XichmEdkqIl9x0uVkRKRARFaKyMb+z9ddTjsNRkRSReR/ReRFp11ORkQmiMgL/Z+zd0WkZrTacnUQFJFFwI1AmVKqFHjQYaUoRKQA+Htgh9Mug3gVmKOUmgu0APc5JSIiqcCPgauB2cBiEZntlM8geoF7lVKzgUuAz2vkdoK7gHedlhiCHwAvK6VmAWWMoqOrgyDwWeBbSqkQgFJqn8M+g/k+8GVAqwu3Sqk/K6V6+399E8h3UOdiYKtSaptS6hjwLJE/bI6jlNqjlGrq/3cHkS9ynrNWHyAi+cC1wGNOu5yMiOQA84HHYWCd8sOj1Z7bg6APmCciDSKyWkQuclroBCJyI7BLKdXstMtp+CfgTw62nwfsPOn3IBoFmhOISCFwIdDgsMrJPETkj2yfwx6D8QL7gf/uP1V/TETOGa3G0karYl0QkdeAaUMULSPy//cQOVW5CPi1iBQpm26Zn8btq0ROhR3hVG5Kqd/377OMyCnf03a6nWmISBbwP8DdSqmjTvsAiMh1wD6lVKOILHRYZzBpQAVwp1KqQUR+AHwF+PfRauysRil1RbwyEfks8Jv+oPeWiPQRma+430k3EbmAyF/DZhGByOlmk4hcrJSyJcnaqfoNQEQ+AVwHXG7XH4047AIKTvo9v3+bFojIGCIB8Gml1G+c9jmJDwE3iMg1QAYwXkSeUkrd6rAXREbzQaXUiVHzC0SC4Kjg9tPh3wGLAETEB4xFg4nkSql3lFJTlVKFSqlCIh+KCrsC4OkQkauInEbdoJTqcljnbaBYRLwiMhb4GPAHh50AkMhfsMeBd5VS33Pa52SUUvcppfL7P18fA/6iSQCk/3O+U0RK+jddDmwcrfbO+pHgafgv4L9EZANwDFjq8KjmTOFhIB14tX+k+qZS6jNOiCilekXkX4BXgFTgv5RSuqwc9CHgNuAdEVnXv+2rSqkVzimdMdwJPN3/h20b8MnRasjMGDEYDK7G7afDBoPB5ZggaDAYXI0JggaDwdWYIGgwGFyNCYIGg8HVmCBosB0RCYvIOhHZICJ/tJqFRkQ+ISIPJ1nvRN1fHY16DfphgqDBCbqVUuVKqTlAO/B5p4VOIBFSiExbNLgAEwQNTlNPf8IDETlPRF4WkUYReUNEZvVvv74/ycX/ishrIpJ7qgpF5D9E5JciUi8iW0Tkn/u3Z4nI6yLSJCLv9CepQEQK+/MRPglsIDLLI7N/tGrmRZ/luH3GiMFB+nMBXk5/yiQia0p8Rim1RUSqgZ8AlwF1wCVKKSUinyIyZe/e01Q/l0hijHOA/xWRl4B9wIeVUkf7k9S+KSInptgVE5kx9Ga/2y1KqfJk/V8N+mKCoMEJMvunkeURybH3an+mlUuB5/un4kFkah5EkiI8JyLTiczvDiTQxu+VUt1At4isJJJ38CXgmyIyn0j6qDzgxKhy+4kAaHAX5nTY4ATd/aOsmYAQuSaYAhzuv1Z44vV/+vf/EfCwUuoC4NNEsp6cjsHzQRXwj8AUoLK//baT6np/BP8fwxmMCYIGx+jPQPMFIqe2XUBARG6BgRsUZf275vBBeqylCVZ/o4hkiMgkYCGRbDM5RHLoHe9fWmHmKd5/vD8NluEsxwRBg6Mopf4XWA8sJjJSu11EmgE/H6TJ/w8ip8mNJJ7qbD2wkkj6//9USu0mkvy1SkTeAT4ObDrF+38GrDc3Rs5+TBYZw1mHiPwH0KmU0mrhLIOemJGgwWBwNWYkaDAYXI0ZCRoMBldjgqDBYHA1JggaDAZXY4KgwWBwNSYIGgwGV2OCoMFgcDX/P5A2FwkOT3uYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalue 6.900000000000039: distance from unit circle 5.900000000000039\n",
      "Eigenvalue 0.49999999999999867: distance from unit circle 0.5000000000000013\n",
      "Eigenvalue 2.099999999999993: distance from unit circle 1.099999999999993\n"
     ]
    }
   ],
   "source": [
    "dmdc.plot_eigs(show_axes=True, show_unit_circle=True, figsize=(5, 5))\n",
    "\n",
    "for eig in dmdc.eigs:\n",
    "    print('Eigenvalue {}: distance from unit circle {}'.format(eig, np.abs(1-np.linalg.norm(eig))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We seek a confirmation about this fact by computing the eigenvalues of the operator `A` which we used when we constructed the system, and we compare them with the eigenvalues of `dmd.atilde` which holds the approximation of `A` built by DMD with control."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues of A: [2.1 0.5 6.9] ; eigenvalues of A_tilde:  [6.9 0.5 2.1]\n"
     ]
    }
   ],
   "source": [
    "print('Eigenvalues of A:', np.linalg.eigvals(A), '; eigenvalues of A_tilde: ', np.linalg.eigvals(dmdc.atilde))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}