msm-vampe-convergence.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:29.006887Z",
     "start_time": "2021-09-30T11:28:28.999970Z"
    },
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "# Papermill params\n",
    "ratio = 0.9          # Train-Test split ratio\n",
    "attempt = 0        # Number of times to run\n",
    "width = 512\n",
    "depth = 2\n",
    "learning_rate = 1e-2\n",
    "dropout = 0.05\n",
    "regularization = 1e-7\n",
    "epsilon = 1e-7\n",
    "attempts = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "colors = sns.color_palette(\"husl\", 8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural network\n",
    "\n",
    "In this notebook we evaluate convergence of the individual ensembles by comparing timescales for limited numbers of trajectories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:35.439766Z",
     "start_time": "2021-09-30T11:28:29.011591Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/thomas/anaconda3/envs/tf/lib/python3.5/site-packages/pyemma/__init__.py:130: UserWarning: You are not using the latest release of PyEMMA. Latest is 2.5.7, you have 2.5.4.\n",
      "  .format(latest=latest, current=current), category=UserWarning)\n"
     ]
    }
   ],
   "source": [
    "%run model.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:35.445767Z",
     "start_time": "2021-09-30T11:28:35.441629Z"
    }
   },
   "outputs": [],
   "source": [
    "def statdist(X: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    Calculate the equilibrium distribution of a transition matrix.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X\n",
    "        Row-stochastic transition matrix\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    mu\n",
    "        Stationary distribution, i.e. the left\n",
    "        eigenvector associated with eigenvalue 1.\n",
    "    \n",
    "    \"\"\"\n",
    "    ev, evec = eig(X, left=True, right=False)\n",
    "    mu = evec.T[ev.argmax()]\n",
    "    mu /= mu.sum()\n",
    "    return mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:35.508600Z",
     "start_time": "2021-09-30T11:28:35.447609Z"
    }
   },
   "outputs": [],
   "source": [
    "def test_split(data: MaybeListType[np.ndarray], lag: int, p: float=0.1, mask: Optional[np.ndarray]=None):\n",
    "    data = make_list(data)\n",
    "    lengths = np.array([len(d) for d in data])\n",
    "    nframes = lengths.sum()\n",
    "    \n",
    "    inds = np.empty((nframes, 4), dtype=np.int)\n",
    "    inds[:, 0] = np.repeat(np.arange(len(data), dtype=np.int), lengths)\n",
    "    inds[:, 1] = np.concatenate([np.arange(n) for n in lengths])\n",
    "    inds[:, 2] = np.arange(nframes, dtype=np.int)\n",
    "    inds[:, 3] = np.zeros_like(inds[:, 0]) if mask is None else mask\n",
    "    inds = unflatten(inds, lengths=[lengths])\n",
    "    \n",
    "    # Local (frame) shuffling\n",
    "    shuf_traj_inds = [np.random.choice(\n",
    "        d[:, 1], size=d.shape[FRAMES], replace=False) for d in inds]\n",
    "    \n",
    "    # Sort out too short trajectories, split out lagged part\n",
    "    n_pairs = 0\n",
    "    xt, xttau = [], []\n",
    "    for i, traj in enumerate(inds):\n",
    "        n_points = traj.shape[FRAMES]\n",
    "\n",
    "        # We'll just skip super short trajectories for now\n",
    "        shuf_traj_inds[i] = shuf_traj_inds[i][shuf_traj_inds[i] < (n_points - lag)]\n",
    "        if n_points <= lag:\n",
    "            continue\n",
    "                \n",
    "        n_pairs += n_points - lag\n",
    "        xt.append(traj[:n_points - lag][shuf_traj_inds[i]])\n",
    "        xttau.append(traj[lag:n_points][shuf_traj_inds[i]])\n",
    "        \n",
    "    # Shuffle externally\n",
    "    shuf_full_inds = np.random.choice(\n",
    "        np.arange(n_pairs, dtype=np.int), size=n_pairs, replace=False)\n",
    "    xt_shuf = np.vstack(xt)[shuf_full_inds]\n",
    "    xttau_shuf = np.vstack(xttau)[shuf_full_inds]\n",
    "    \n",
    "    # These are the entries for the test set\n",
    "    n_frames_test = int(xt_shuf.shape[FRAMES] * p)\n",
    "    inds_t = xt_shuf[:n_frames_test]\n",
    "    inds_ttau = xttau_shuf[:n_frames_test]\n",
    "    data_flat = np.vstack(data)\n",
    "    test_xt, test_xttau = data_flat[inds_t[:, 2]], data_flat[inds_ttau[:, 2]]\n",
    "    \n",
    "    # Mask out unwanted frames with NaNs\n",
    "    mask_pair = xt_shuf[:, 3] | xttau_shuf[:, 3]\n",
    "    test_xt[mask_pair[:n_frames_test]] = np.nan\n",
    "    test_xttau[mask_pair[:n_frames_test]] = np.nan\n",
    "    \n",
    "    # We can't just remove our test frame pairs, as the training set\n",
    "    # would then be out of sync! So we replace the test samples with\n",
    "    # NaNs instead, we can check for those later in the DataGenerator.\n",
    "    data_flat[np.union1d(inds_t[:, 2], inds_ttau[:, 2])] = np.nan\n",
    "    data_train_valid = unflatten(data_flat, lengths=[lengths])\n",
    "        \n",
    "    return data_train_valid, (test_xt, test_xttau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "### Trajectories\n",
    "Trajectories were acquired in five rounds of 1024 simulations each, totalling 5119 runs (one simulation failed to run) at 278 K in the $NVT$ ensemble. Postprocessing involved removing water, subsampling to 250 ps timesteps, and making molecules whole."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:35.552417Z",
     "start_time": "2021-09-30T11:28:35.510575Z"
    }
   },
   "outputs": [],
   "source": [
    "sim_names = (\"holo\", \"control\", \"apo\")\n",
    "top, trajs = {}, {}\n",
    "trajs = {k: sorted(glob(\"trajectories/{0}/r?/traj*.xtc\".format(k))) for k in sim_names}\n",
    "top = {k: \"trajectories/{0}/topol.gro\".format(k) for k in sim_names}\n",
    "KBT = 2.311420 # 278 K\n",
    "nres = 42\n",
    "traj_rounds = {\n",
    "    \"apo\": [1024, 1023, 1024, 1024, 1024],\n",
    "    \"holo\": [1023, 1024, 32],\n",
    "    \"control\": [1024, 1023]\n",
    "}\n",
    "\n",
    "# This is only really necessary for the residues in the plots\n",
    "topo = md.load_topology(top[\"apo\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use minimum distances as features for the neural network:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:56.293417Z",
     "start_time": "2021-09-30T11:28:35.554112Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='Obtaining file info'),), layout=Layout(max_width='35%', min_width='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='Obtaining file info'),), layout=Layout(max_width='35%', min_width='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='Obtaining file info'),), layout=Layout(max_width='35%', min_width='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r"
     ]
    }
   ],
   "source": [
    "allpairs = np.asarray(list(itertools.combinations(range(nres), 2)))\n",
    "inpcon = {}\n",
    "for k in sim_names:\n",
    "    feat = pe.coordinates.featurizer(top[k])\n",
    "    feat.add_residue_mindist(residue_pairs=allpairs)\n",
    "    inpcon[k] = pe.coordinates.source(trajs[k], feat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:56.321923Z",
     "start_time": "2021-09-30T11:28:56.295360Z"
    }
   },
   "outputs": [],
   "source": [
    "lengths, nframes = {}, {}\n",
    "for i, k in enumerate(sim_names):\n",
    "    lengths[k] = sort_lengths(inpcon[k].trajectory_lengths(), traj_rounds[k])\n",
    "    nframes[k] = inpcon[k].trajectory_lengths().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:56.340584Z",
     "start_time": "2021-09-30T11:28:56.323963Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t\tholo\t\tcontrol\t\tapo\n",
      "Trajs: \t\t2079\t\t2047\t\t5119\n",
      "Frames: \t1225868\t\t1114503\t\t1259172\n",
      "Time: \t\t306.467 µs\t278.626 µs\t314.793 µs\n"
     ]
    }
   ],
   "source": [
    "print(\"\\t\\t\" + \"\\t\\t\".join(sim_names))\n",
    "print(\"\\n\".join((\n",
    "    \"Trajs: \\t\\t\" + \"\\t\\t\".join(\"{0}\".format(len(trajs[k])) for k in sim_names),\n",
    "    \"Frames: \\t\" + \"\\t\\t\".join(\"{0}\".format(nframes[k]) for k in sim_names),\n",
    "    \"Time: \\t\\t\" + \"\\t\".join(\"{0:5.3f} µs\".format(inpcon[k].trajectory_lengths().sum() * 0.00025)\n",
    "                           for k in sim_names)\n",
    ")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VAMPNet\n",
    "VAMPNet[1] is composed of two lobes, one reading the system features $\\mathbf{x}$ at a timepoint $t$ and the other after some lag time $\\tau$. In this case the network reads all minimum inter-residue distances (780 values) and sends them through 5 layers with 256 nodes each. The final layer uses between 2 and 8 *softmax* outputs to yield a state assignment vector $\\chi: \\mathbb{R}^m \\to \\Delta^{n}$ where $\\Delta^{n} = \\{ s \\in \\mathbb{R}^n \\mid 0 \\le s_i \\le 1, \\sum_i^n s_i = 1 \\}$ representing the probability of a state assignment. One lobe thus transforms a system state into a state occupation probability. We can also view this value as a kind of reverse ambiguity, i.e. how sure the network is that the system is part of a certain cluster. These outputs are then used as the input for the VAMP scoring function. We use the new enhanced version with physical constraints[2], particularly the ones for positive entries and reversibility.\n",
    "\n",
    "[1] Mardt, A., Pasquali, L., Wu, H. & Noé, F. VAMPnets for deep learning of molecular kinetics. Nat Comms 1–11 (2017). doi:10.1038/s41467-017-02388-1\n",
    "\n",
    "[2] Mardt, A., Pasquali, L., Noé, F. & Wu, H. Deep learning Markov and Koopman models with physical constraints. arXiv:1912.07392 [physics] (2019)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data preparation\n",
    "We use minimum residue distances as input ($\\frac{N(N-1)}{2}$ values, where $N$ is the number of residues) and first normalize the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:28:56.346392Z",
     "start_time": "2021-09-30T11:28:56.342412Z"
    }
   },
   "outputs": [],
   "source": [
    "for k in sim_names:\n",
    "    filename = \"intermediate/mindist-all-{0}.npy\".format(k)\n",
    "    if not os.path.exists(filename):\n",
    "        print(\"No mindist file for {0} ensemble, calculating from scratch...\".format(k))\n",
    "        con = np.vstack(inpcon[k].get_output())\n",
    "        np.save(filename, con)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:30:25.009002Z",
     "start_time": "2021-09-30T11:28:56.348055Z"
    }
   },
   "outputs": [],
   "source": [
    "idx = np.triu_indices(nres, k=1)\n",
    "mat = np.zeros((nres, nres), dtype=np.int)\n",
    "full_flat, full_data = {}, {}\n",
    "for k in sim_names:\n",
    "    raw = np.load(\"intermediate/mindist-all-{0}.npy\".format(k))\n",
    "    mat[idx] = np.arange(raw.shape[1])\n",
    "    redinds = mat[np.triu_indices_from(mat, k=3)]\n",
    "    full_flat[k] = ((raw - raw.mean(axis=0)) / raw.std(axis=0))[:, redinds]\n",
    "    full_data[k] = unflatten(full_flat[k], lengths[k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Neural network hyperparameters\n",
    "To allow for a larger hyperparameter search space, we use the self-normalizing neural network approach by Klambauer *et al.* [2], thus using SELU units, `AlphaDropout` and normalized `LeCun` weight initialization. The other hyperparameters are defined at the beginning of this notebook.\n",
    "\n",
    "[2] Klambauer, G., Unterthiner, T., Mayr, A. & Hochreiter, S. Self-Normalizing Neural Networks. arXiv.org cs.LG, (2017)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:30:25.017939Z",
     "start_time": "2021-09-30T11:30:25.011136Z"
    }
   },
   "outputs": [],
   "source": [
    "activation = \"selu\"                 # NN activation function\n",
    "init = \"lecun_uniform\"              # NN weight initialization\n",
    "lag = 20                            # Lag time\n",
    "n_epoch = 100                       # Max. number of epochs\n",
    "n_epoch_s = 10000                   # Max. number of epochs for S optimization\n",
    "n_batch = 10000                     # Training batch size\n",
    "n_dims = full_data[k][0].shape[1]   # Input dimension\n",
    "nres = 42                           # Number of residues\n",
    "epsilon = 1e-7                      # Floating point noise\n",
    "dt = 0.25                           # Trajectory timestep in ns\n",
    "steps = 6                           # CK test steps\n",
    "bs_frames = 1000000                 # Number of frames in the bootstrap sample\n",
    "n_tries = 1                         # Number of training attempts for each model, we pick the best scoring one\n",
    "\n",
    "# outsizes = np.array([4, 2, 3, 5, 6])\n",
    "outsizes = np.array([8, 10, 12, 14, 16])\n",
    "lags = np.array([1, 2, 5, 10, 20, 50, 100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bound ensemble\n",
    "We filter out unbound frames to avoid issues with potential non-equilibrium conditions, as we see essentially no unbinding events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:30:28.660804Z",
     "start_time": "2021-09-30T11:30:25.019896Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='Obtaining file info'),), layout=Layout(max_width='35%', min_width='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='Obtaining file info'),), layout=Layout(max_width='35%', min_width='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r"
     ]
    }
   ],
   "source": [
    "inter_pairs = list(itertools.product(range(nres), (nres,)))\n",
    "intercon = {}\n",
    "for k in (\"holo\", \"control\"):\n",
    "    feat = pe.coordinates.featurizer(top[k])\n",
    "    feat.add_residue_mindist(residue_pairs=inter_pairs)\n",
    "    intercon[k] = pe.coordinates.source(trajs[k], feat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:30:30.247119Z",
     "start_time": "2021-09-30T11:30:28.662613Z"
    }
   },
   "outputs": [],
   "source": [
    "inter_con = {}\n",
    "for k in (\"holo\", \"control\"):\n",
    "    filename = \"intermediate/mindist-inter-all-{0}.npy\".format(k)\n",
    "    if not os.path.exists(filename):\n",
    "        print(\"No mindist file for {0} ensemble, calculating from scratch...\".format(k))\n",
    "        inter_con[k] = np.vstack(intercon[k].get_output())\n",
    "        np.save(filename, inter_con)\n",
    "    else:\n",
    "        inter_con[k] = np.load(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:30:30.326644Z",
     "start_time": "2021-09-30T11:30:30.249073Z"
    }
   },
   "outputs": [],
   "source": [
    "cutoff = 0.5  # Compromise given the above plot\n",
    "masks = {k: (inter_con[k] < cutoff).any(axis=1) for k in (\"holo\",)}\n",
    "masks.update({k: np.full(full_flat[k].shape[0], True) for k in (\"apo\", \"control\")})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T11:31:23.505924Z",
     "start_time": "2021-09-30T11:30:30.328502Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading existing input file for holo ensemble...\n",
      "Reading existing input file for apo ensemble...\n",
      "Reading existing input file for control ensemble...\n"
     ]
    }
   ],
   "source": [
    "input_data, input_flat, test_data = {}, {}, {}\n",
    "for k in (\"holo\", \"apo\", \"control\"):\n",
    "    filename = \"intermediate/input-mask-{0}.npz\".format(k)\n",
    "    if not os.path.exists(filename):\n",
    "        print(\"No input file for {0} ensemble, resplitting...\".format(k))\n",
    "        input_data[k], test_data[k] = test_split(full_data[k], lag=lag, mask=~masks[k])\n",
    "        input_flat[k] = np.vstack(input_data[k])\n",
    "        input_flat[k][~masks[k]] = np.nan\n",
    "        np.savez(filename, data=input_flat[k], test_t=test_data[k][0], test_ttau=test_data[k][1])\n",
    "    else:\n",
    "        print(\"Reading existing input file for {0} ensemble...\".format(k))\n",
    "        raw = np.load(filename)\n",
    "        input_flat[k], test_data[k] = raw[\"data\"], (raw[\"test_t\"], raw[\"test_ttau\"])\n",
    "        input_data[k] = unflatten(input_flat[k], lengths=lengths[k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run\n",
    "We evaluate the Koopman model with varying numbers of trajectories and compare the implied timescales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T15:24:41.888424Z",
     "start_time": "2021-09-30T14:46:54.648034Z"
    },
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00232: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00818: early stopping\n",
      "Lag loss: -1.9996654987335205\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00427: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00389: early stopping\n",
      "Lag loss: -1.9883699417114258\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00070: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00302: early stopping\n",
      "Lag loss: -1.9999778270721436\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00051: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00452: early stopping\n",
      "Lag loss: -1.9951941967010498\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00537: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00134: early stopping\n",
      "Lag loss: -1.9904100894927979\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00090: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00344: early stopping\n",
      "Lag loss: -1.9924297332763672\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00283: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00193: early stopping\n",
      "Lag loss: -1.9824985265731812\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00210: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00431: early stopping\n",
      "Lag loss: -1.991247296333313\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00051: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00442: early stopping\n",
      "Lag loss: -1.9851287603378296\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00271: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00239: early stopping\n",
      "Lag loss: -1.987882137298584\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00121: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00613: early stopping\n",
      "Lag loss: -1.9918204545974731\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00569: early stopping\n",
      "Restoring model weights from the end of the best epoch\n",
      "Epoch 00498: early stopping\n",
      "Lag loss: -1.9835622310638428\n"
     ]
    }
   ],
   "source": [
    "n = 2\n",
    "for k in (\"holo\", \"control\"):\n",
    "    filename = \"intermediate/k-conv-{0}-{1}-1.npy\".format(k, n)\n",
    "    n_trajs = np.array([1024, 2048])\n",
    "    k_conv = np.empty((len(n_trajs), attempts, n, n))\n",
    "    for j, nt in enumerate(n_trajs):\n",
    "        for i in range(3):\n",
    "            generator = DataGenerator(input_data[k][:nt])\n",
    "            print(\"Analysing trajs={0} n={1} i={2}...\".format(j, n, i), end=\"\\r\")\n",
    "            koop = KoopmanModel(n=n, network_lag=lag, verbose=0, nnargs=dict(\n",
    "                width=width, depth=depth, learning_rate=learning_rate,\n",
    "                regularization=regularization, dropout=dropout,\n",
    "                batchnorm=False, lr_factor=2e-2))\n",
    "            koop.load(\"models/model-ve-{0}-{1}-{2}.hdf5\".format(k, n, i))\n",
    "            koop.generator = generator\n",
    "            k_conv[j, i] = koop.estimate_koopman(lag=50)\n",
    "            del koop\n",
    "            gc.collect()\n",
    "    np.save(filename, k_conv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Concatenate partial runs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T15:29:47.828962Z",
     "start_time": "2021-09-30T15:29:47.823050Z"
    }
   },
   "outputs": [],
   "source": [
    "n = 2\n",
    "k_conv = {}\n",
    "for k in (\"holo\", \"control\"):\n",
    "    k_conv[k] = np.concatenate([np.load(\"intermediate/k-conv-{0}-{1}.npy\".format(k, n))[:, :2],\n",
    "                                np.load(\"intermediate/k-conv-{0}-{1}-1.npy\".format(k, n))[:, :3]], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T15:30:19.277607Z",
     "start_time": "2021-09-30T15:30:19.271271Z"
    }
   },
   "outputs": [],
   "source": [
    "its_conv = {}\n",
    "for k in (\"holo\", \"control\"):\n",
    "    n_rounds = len(k_conv[k])\n",
    "    its_conv[k] = np.empty((n_rounds, 5, n - 1))\n",
    "    for j in range(n_rounds):\n",
    "        for i in range(5):\n",
    "            its_conv[k][j, i] = -50 * dt / np.log(np.sort(np.linalg.eigvals(k_conv[k][j, i]))[::-1])[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-30T15:37:29.028133Z",
     "start_time": "2021-09-30T15:37:28.519232Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALoAAAEBCAYAAADPdL/zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFK1JREFUeJztnXuUHUWdxz9fw2uJhOURjBhwCD5W8YFkFAgsA3oU1vWBgBKPIiqPI4eHgeDjQCAjRzQIrkBYZUFNEDgGDBvWuBEESUggLBhXTkR5BTK6QJCESAIEAoHf/lF1SedO33u778yd3Dv1+5zTp6er61ddNfO9daurevorM8Nxhjuv29wVcJyhwIXuJIEL3UkCF7qTBC50Jwlc6E4SuNCdJCgsdEl9kqzG9mSNmAmS5klaLWmdpKWSJkkaUec6H5O0QNIaSc9JulvSsQ3qdqyke2L+NTH+Y0Xb5gx/tiiZfw1wcU76c9UJkj4J3AC8CFwHrAY+DvwAOAD4dE7MKcB04GngGuAl4ChgpqR3m9mZOTEXAZOBx4Arga2AicBcSaea2WUl2+gMR8ys0Ab0AX0F844CngLWA92Z9G2AxYABE6tiuggfiqeBrkz6DsCyGLN/VcyEmL4M2KGqrKdjeV1F6nzooYea07E0/Pu2aox+FDAamGVmSyqJZvYiMCUenlQV82Vga+AyM+vLxPwd+E48/EpVTOX4/JivEtMH/Hss70tFKrxq1aoi2ZwOpazQt5b0eUlnSfqqpENqjLc/GPc35ZxbCKwDJkjaumDMr6vyDCTGSZCyY/QxwNVVacslfcnMbs+kvT3uH6ouwMw2SFoO7AWMA+4vELNC0vPAWEnbmtk6SSOBNwHPmdmKnLo+HPdvK9IwZ3hTpkefAXyIIPaRwLuB/yCMh38t6b2ZvNvH/ZoaZVXS/7GJmO2r9mWusQmSTpS0RNKSlStX1spWmt7eXiQ1vfX29g5aXZxIkYF8vQ24iHBDOCeT9lBMe0uNmMoN6X6ZtJdi2hY1Yp6I58fE413j8WM18m8Zz79YpB3jx48f+C1RAXp6eqynp2dIrpUQQ3IzenncH5RJq+59qxlVla9MzNqC+Rv1+E5CDIbQn4r7kZm0B+O+3/hY0hbAHsAG4NGCMW+M5T9mZusAzOx54HHg9fF8NW+N+35jfic9BkPo+8d9VrS3xf1hOfkPArYFFpvZ+oIx/1KVZyAxToIUErqkvSTtmJP+ZqCy8nhN5tRsYBUwUVJ3Jv82wLfj4Y+qiptBWGA6RVJXJmYH4Kx4eHlVTOX47JivEtMFnBzLm1G3cU4SFJ1e/DTwTUnzgeXAs8CewL8SVjvnEW5KATCztZJOIAh+gaRZhEcAPkGYRpxNeCyATMxySV8DLgWWSLqOjY8AjAW+b2Z3VcUslvRvwBnAUkmzCY8AHA3sCJxqmcUnJ12KCn0+QaDvIwxVRgLPAHcQ5tWvNtv0v6zN7EZJPcDZwJGED8Qygigvrc4fY6ZL6gPOBL5A+Mb5MzDFzK7Kq5iZTZa0FDgFOBF4Ffhf4EIz+1XB9jnDnEJCt7AYdHvDjP3j7gQ+WjJmLjC3ZMxVQO4HwXHAn0d3EsGF7iSBC91JAhe6kwQudCcJXOhOErjQnSRwoTtJ4EJ3ksCF7iSBC91JAhe6kwQudCcJXOhOErjQnSRwoTtJ4EJ3ksCF7iSBC91JAhe6kwQudCcJXOhOErjQnSRwoTtJ4EJ3ksCF7iSBC91JAhe6kwRNC13SMRmL9ONr5Gm53bmkEdF2famkF6Id+zxJE5ptmzP8aEroknYjWJn3s0bP5DmF8FbcdxFMAq4kGGzNjLbmeTEXATOBN8b81xDc7+bG8qrzC5hFsF3fimBKMIfgqrEw2rQ7TnlXOkDArcAjwIUE57fjq/J0MQR258BnY8ydwDaZ9PcT3C6eArYr0i53petoWuJKdxrBjflLwPM18gyV3XnFZn2KBfv1SszvCI4aowmOGU7ilBK6pHcA04BLzGxhnawttzuP9uoTCHbri0pcx0mQwkKPtolXA39lo3lWLeranRO+CcZK2jaW3Yzd+VuAEcCjZrahYIyTKGV69HMJHkZfNLMXGuQdCrvztrVId9qPovaLHyD04v2c4ZpEcd/PsKsBZfI3vIaZXWFm3WbWPXr06JJVcTqJhkLPDFkeAs4pWO5Q2J03Y8PuJEqRHv31hHHuO4AXM4tEBkyNea6MaRfH46GwO18GvAKMix/GIjFOohSxX1wP/KTGuX0I4/Y7COKuDGtuAw4gWJdXD3Xq2Z0fE2Oq3Z77xZjZekmLgX+O2/yC13FSpMhke60N6CV/wWgP2mfBaFSRtviCUUfT8O9b1Dm67IdnqOzOZwFHxHL/IGkusFOMGQGcYGZrcZKnJUKHobE7NzOT9FlgMWE19lRCz78Q+LaZLR70hjkdiczKzvANT7q7u23JkiUtv87BBx8MwIIFC1p+rYRQowz+PLqTBC50Jwlc6E4SuNCdJHChO0ngQneSwIXuJIEL3UkCF7qTBC50Jwlc6E4SuNCdJHChO0ngQneSwIXuJIEL3UkCF7qTBC50Jwlc6E4SuNCdJHChO0ngQneSwIXuJIEL3UkCF7qTBC50Jwlc6E4SuNCdJCjjSneBpN9K+r+MFfkfJE2VtFONmAnRrny1pHXRxnySpBF1rtNyW3UnPcr06KcTLFluAS4BrgU2EMwAlkbb9NeI9uQLCXblcwimuFsR7Mxn5V1gKGzVnUQp4hYQXy29TY308wmuEz/MpI0iuE2sB7qzZRDeZW7AxKpyuhgCl4xamztelGPq1KkWf/dNbVOnTh3M6gyeRbplLMiruD7u35pJO4pgTz7LzF576XgsY0o8PIlNGSpbdWcQ6O3trSmqnp4eenp66gqvt7d3SOs7GDejH4/7pZm0enbnCwm25hOizXmRmEGxVXfSpbS1i6QzCZaM2wPdwIEEkU/LZKtnkb5B0nJgL2AccH+BmBWSXrNVN7N1TdqqV7flRIKFDLvvvnutbM4woBkPozOBN2SObyLYpmc9xltleT4y5lvX5DU2wcyuAK6AYO1SK5/T+ZQeupjZGDMTMIbgCDeO4Ai3T4limrFIHwpbdWeY0vQY3cz+ZmZzgI8QLA9/ljndjH35UNiqO4ky4JtRM/sLwVZxL0k7x+R6FulbEAx3NwCPZk4Nha26kyiD9QjArnH/StxXbMkPy8l7ELAtsNjM1mfS68XUs1UvG+MkSCGhS/onSWNy0l8n6XxgF4JwK3PZs4FVwERJ3Zn82wDfjoc/qipuBmGB6RRJXZmYHYCz4uHlVTGV47NjvkpMF3ByLG9GkTY6w5uisy6HARdKWgg8Qlh1fAPQQ7gZfRI4oZLZzNZKOoEg+AWSZgGrgU8QphFnA9dlL2BDZ6vuJEhRod9KmIY7AHgvYcruecL492rgUjNbnQ0wsxsl9QBnA0cSlv+XEUR5qVl/y2obAlt1J00KCd3M7iMMBUphZncCHy0ZM5fwYFeZmKuA3A+C44A/j+4kggvdSQIXupMELnQnCVzoThK40J0kcKE7SeBCd5LAhe4kgQvdSQIXupMELnQnCVzoThI08xaAZNj3P68d9DLvX/W3lpUNcPcRn2tJuZ2O9+hOErjQnSRwoTtJ4EJ3ksCF7iSBC91JAhe6kwQudCcJXOhOErjQnSRwoTtJ4EJ3ksCF7iRB0ddG7yTpeElzJC2LztFrJN0h6ThJueW4c7TTLhTt0T9NcGXeF7gbuBi4geDw/GPgeknKBrhztNNOFH0e/SHCu83/28xerSRKOgu4h/Ba6CMI4kfSKILoXgEOrpjqSjqH4EBxlKSJZjYrU1YXcBHhPerdlfeaSzoP+B0wWdIN2XekS5oATCa8s/39FSMCSRcCvwcukvQrf0e6U6hHN7PbzGxuVuQx/Uk2uk4cnDnlztFOWzEYN6Mvx/2GTJo7RzttxYCEHh3mvhAPs2Kr6xwNLCcMm8YVjFlBcNgYK2nbeO0BO0c76TDQHn0a4cZxnpndnElvlXN0Nt+AnaMlnShpiaQlK1eurJXNGQY0LXRJpxFuBB8AjikbHveb1TnazK4ws24z6x49enTJYp1OoimhSzoZuIRgpHVItVEX7hzttBmlhS5pEnAZcB9B5E/mZHPnaKetKCV0Sd8gLPjcSxD5UzWyunO001YUFnpc7JlGWIj5kJmtqpPdnaOdtqLQymh81uQ8wkrnIuC0qhV/gD4zmwnuHO20H0UfAdgj7kcAk2rkuZ3wzAngztFOe1HUOboX6C1buDtHO+2Cv2Q0YVr1otNWvki12Zeo+j9eOEngQneSwIXuJIEL3UkCF7qTBC50Jwlc6E4SuNCdJHChO0ngQneSwIXuJIEL3UkCF7qTBC50Jwlc6E4SuNCdJHChO0ngQneSwIXuJIEL3UkCF7qTBC50Jwlc6E4SuNCdJHChO0ngQneSwIXuJEFRi/SjJE2XtEjSWkkm6ZoGMW6P7rQNRXv0KYTXMu9NsFOpi9ujO+1GUaGfTvAWGkV/x+dNyLFHP87Mvkb4kNxFtEeviuliU3v0k83sdOA9BPvzyZL2r4rJ2qO/x8xON7OTgfGxnIuyzhlO2hS1SJ9vZg/nvbw/B7dHd9qOVtyMuj2603a0Quhuj+60Ha0QekfYo4NbpKfE5phHbwt7dHCL9JRohYdRs/boO8eYp+vEdIQ9+mPX3cAT18+pm+eeIz9f89yun/kUY48+crCrlTStEPqDQDdhfPz77IkG9ug7x5i7qmJy7dElPQ68SdIbc8bpm9UefezRR7pQ24xWDF3cHt1pO1rRo88GLiDYo0+vzKUXsEf/OsEefUZlLr2APfoxBHv0Gytz6W6PPjR02vCsqEX64cDh8XBM3O8vaWb8eZWZnQluj54KnTY8K9qj7w1UP1w1jo1z4X8h2JoDbo/utB8ts0h3e3SnnfDn0Z0kcKE7SeBCd5LAhe4kgQvdSQIXupMELnQnCVzoThK40J0kcKE7SeBCd5LAhe4kgQvdSQIXupMELnQnCVzoThK40J0kcKE7SeBCd5LAhe4kgQvdSQIXupMELnQnCVzoThK40J0kcKE7SeBCd5LAhe4kwbAQuqSxkn4q6QlJ6yX1Sbo4vl/dcVpiBDCkSNoTWAzsAvwX8ADwAeCrwGGSDjCzPF8kJyGGQ4/+Q4LITzOzw83sm2b2QeAHBOOB8zdr7Zy2oKOFLmkc8BGgj2CLnmUqwYj3mGjA6yRMRwudjRbovzGzV7MnzOxZ4E6COdh+Q10xp73o9DF6TWv1yMOEHv9twG+rT0o6kWAJA/CcpAcHvYZDjKhtkDUcqNG+m8wsz53wNTpd6AOySjezK4ArBrtSTvvR6UOXRjRrre4MMzpd6M3YsTsJ0ulCr4yp31bj/Ga1SnfaB+VYfnYMcbFoGWF6cc/szIuk7YAVhA/zaDN7frNU0mkLOrpHN7NHgN8AXQRb9CzfAkYCP2tG5JKOkjRd0iJJayWZpGsaxEyQNE/SaknrJC2VNEnSiJy8e0vqlXSnpBWSXpL0uKSfS9qnYB13lvRkrNsdJdu3k6TjJc2RtEzSC5LWSLpD0nGScrVRpo05sZJ0S6yvScqdDJH0bknXZur1uKT5ko6uVa+GmFlHb8CewN8IN5w3At8FbovHDwI7NVnuvbGMZ4H748/X1Mn/SWAD8BzwE+BCwuMIBvwiJ///xHNLCItdFwA3x7SXgU8VqOMNsX4G3FGyfV+JcU8A18bf20+BZ2L6bOI3frNtzLnmqTH+hRizRU6ej8f2rwd+AUwDfgw8HWOubOrvubmFOkhi3w2YQRiqvESwbL8E2HEAZR5CGOMLOLie0Ak3vU/FP053Jn0bwnM4BkzM+aO/Jaesz8X8q4Ct6tTvCzHfSU0K/YNRVK+rSh8D/DWWeeRA2lhV7tuBdVG4fXWE/qd4rienXpUObfckhd7qrYDQvxzPX1VDUAbcXuJ6D8WY8TXO7x573h8Thm2lhd7g+mfFMqcPRhsJ6zX3RBFv3UDoLwBrapTzy3q/l3pbR4/R24jKowg35ZxbSOjJJkjaumB5L8f9huoTkgTMJEyZnlGumoXJu/5A2jgFeB9wrJmtb3DtPwGjJB2YTZS0C+Gp1CeAPzcoox8u9MGh5qMIZrYBWE7o1cY1KkjSvsA7gceB+3KyTCJ8wxxnZmubrG+9629BGBbBpqJuqo2S3g+cDUwzsyUFqnA6sBa4VdJ1kr4r6UrCB+BZ4HAze6FEk4DOfwSgXRjQowgV4j+KXB0PzzCzV6rOvxP4DnC5md3aZF0bMQ14FzDPzG7OpJduo6R/ILTnz8B5RS5uZosk7Q9cD3wmc+pZwn3YH4uUU4336ENDw0cR4qPEvyTcAH/PzK6vOr8lQTQrgK+3pJLSacBkwkzKMWXD4z7bxu8Revhjzezl/iG5dfgwsIjwjTaeMEW8J+F+5Hzgt7WmJeuyuW/0OmGj8c3o76h/83hfPP+OGudHAgtinu/XyHMu8Cr9ZyO6GISbUcI6hBGGCGMG2kagJ9b33Jy8feTcjAI7AquBx4Btc+LmxLgvlm2f9+iDQ81HEWLvswfhxu7RnPPbAb8mCON7Zja5xjX2IfSaCzILLkYYGwMcENOeKVt5SZOAywhiPcTMnszJVraN74v1/Va2vrHOb455Xo5pe8fjCcAOwN1mti6nDvPjfny5FvoYfbC4jTD/fRjw86pzBxH++WOhVc04SNqecMO3H3C+mU2pc41bCHPr1bweOJowx/wrwuxHYSR9gzAuvxf4sJnlXQPKt/E+wqJSHkfHev+U0ENX/qe3MmMzukZcJf2lGudrs7mHBZ2wUWzBaCXlFox2YONwoN/Xe4m6ddHk0AU4h42rs3UX15ppY52y+sgfuuxKmNp8BfhI1bndCAtWBny0bFs7+qGuViLpcODweDgGOJTwtbwopq0yszOr8s8GXgRmEcaanyBMy80GPmOZX7ak+YQP0CNArWdobjSzexvUs4swfLnTzA6sl7cq7ljCfPwrwHTyZ1P6zGxmJqZUG+tcu48wfNnSwtRk9ty5hOeUXiV8Qz1A+P0fQfgWmGNmRxRt52ts7t6yXTegl9B71Nr6cmIOAOYBfyes8P2RMC88ok6vVm/7YoF6dtHcIwCN2mfAgoG0sc61K23vtzIaz3+ScN+ykjDuX0v41jipzHWym/foThL4rIuTBC50Jwlc6E4SuNCdJHChO0ngQneSwIXuJIEL3UkCF7qTBP8P9ToL9nVyxp4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 144x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALoAAAEBCAYAAADPdL/zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFFNJREFUeJztnXuwHUWdxz9fgoCgYQGj0YUYg49VfCDEB8Hl4ZbCuqCsRIHahag8VoqHILhaiHi1RFFQQVARVECxDBgKV6yIohACxFVSygYEVCBXBYMkgAkkIZDw2z+6Lwwn5zHn3EfuSX8/VVOT6elfT3fOd/r0dM89X0UExmzsbLKhK2DMWGChmyKw0E0RWOimCCx0UwQWuikCC90UQW2hSxqUFC22+1vEzJA0V9JDklZJWiTpBEkT2lxnP0nzJC2X9KikX0ma1aFusyT9OudfnuP3q9s2s/GzaZf5lwNnN0l/tDFB0ruAK4DHgMuAh4D9gS8DuwPvaRJzLHAu8CBwKfA4MBO4WNJrIuLkJjFnAScB9wIXApsBBwNXSTouIs7rso1mYyQiam3AIDBYM+9E4AFgDTC9kr4FsAAI4OCGmKmkm+JBYGolfRvgrhyzW0PMjJx+F7BNQ1kP5vKm1qnzPvvsE6Zv6fj5jtYYfSYwCZgdEQuHEiPiMeDUfHh0Q8wHgM2B8yJisBLzMPDZfPjBhpih49NzvqGYQeCrubz316nwsmXL6mQzfUq3Qt9c0n9KOkXShyTt3WK8/da8v7rJufnAKmCGpM1rxvykIc9wYkyBdDtGnwx8tyFtsaT3R8T1lbRX5P0fGguIiLWSFgM7AdOAO2rELJG0Ethe0pYRsUrSVsA/Ao9GxJImdf1j3r+8TsPMxk03PfpFwL+QxL4V8BrgG6Tx8E8kva6Sd+u8X96irKH0f+ghZuuGfTfXeAaSjpK0UNLCpUuXtspmNgJqCz0iPhUR10bE3yJiVUTcFhEfBL4EPBsY6OK6Gip2lGPa5o+ICyJiekRMnzRpUpfFmn5iJB5Gz8/7PSppjb1vIxMb8nUTs6Jm/k49vimIkRD6A3m/VSXt93m/3vhY0qbAS4C1wD01Y16Yy783IlYBRMRK4D7gOfl8Iy/L+/XG/KY8RkLou+V9VbTX5v2+TfLvAWwJLIiINTVj/rUhz3BiTInUmWwnzZBs2yT9xaTZjQBOqaRPBJbS3YLRS9iAC0a77rrrMNcszAak4+dbd3rxPcDHJF0HLAYeAXYE/i2Ldy5wVuXmWSHpSGAOME/SbNIrAO8kTSPOIb0WUL3hFkv6CPAVYKGky3j6FYDtgS9GxC8bYhZI+hLwYWCRpDmkVwAOArYFjovK4pMpmDp3A7An8H3gTuDvwBOkHvsa4DBALeJ2J90EDwOrgVuBE4EJba61P3A96WZaCdwMzOpQv1k538ocdz2wX522DW3u0fuajp+vIvwrAADTp0+PhQsXds5oxiPqlMHvo5sisNBNEVjopggsdFMEFropAgvdFIGFborAQjdFYKGbIrDQTRFY6KYILHRTBBa6KQIL3RSBhW6KwEI3RWChmyKw0E0RWOimCCx0UwQWuikCC90UgYVuisBCN0VgoZsisNBNEVjopggsdFMEPQtd0qEVi/QjWuQZdbtzSROy7foiSauzHftcSTN6bZvZ+OhJ6JJ2IFmZr2eNXslzLHAV8GqS3fmFwItIdudntYg5C7gYeGHOfynJ/e6qXF5jfgGzSbbrmwHnAVeSXDXmZ5t2Y+pbpA9tpJ/o/TlwN3AmyXHiiIY8UxkD9wrgkBxzE7BFJf0NJLeNB4Dn1mmXfx+9r+n4+fbSox9PcmN+P+mH95sxVnbnQzbrp0ayXx+KuZnkqDGJ5JhhCqcroUt6JXAGcE5EzG+TddTtzrO9+gyS3foNXVzHFEhtoWfbxO8CfwZO6ZC9rd056Ztge0lb5rJ7sTt/KTABuCci1taMMYXSTY9+GvB64H0RsbpD3rGwO7dFuqlNLaFLeiOpF1/PGa5HRtzuvJdrhC3Si6Gj0CtDlj8An6hZ7ljYnfdiw24KpU6P/hzSOPeVwGOVRaIAPpnzXJjTzs7HY2F3fhewDpiWb8Y6MaZQ6hjqrgG+1eLcLqRx+40kcQ8Na64leYzuW0kbop3d+aE55qJOMRGxRtIC4J/zdl3N65gSqTPZ3moDBmi+YDQmdufUWzCaWKctXjDqazp+vnUt0ru9ecbK7nw28O5c7m8lXQVsl2MmAEdGxApM8YyK0AEi4lxJg8DJJBv1TYDbSauYl7SIOUnSIuBY4CjgSeA3wJkR8eMm+UPSIcAC0mrscaSefz7wmYhYMOINM32JLdIztkjva2yRbgxY6KYQLHRTBBa6KQIL3RSBhW6KwEI3RWChmyKw0E0RWOimCCx0UwQWuikCC90UgYVuisBCN0VgoZsisNBNEVjopggsdFMEFropAgvdFIGFborAQjdFYKGbIrDQTRFY6KYILHRTBBa6KYJuXOk+L+kXkv5SsSL/raRPStquRcyMbFf+kKRV2cb8BEkT2lxn1G3VTXl006OfSLJkuQY4B/gesJZkBrAo26Y/RbYnn0+yK7+SZIq7GcnOfHazC4yFrboplDpuAfmnpbdokX46yXXia5W0iSS3iTXA9GoZpN8yD+DghnKmMgYuGa02O170NSNnkR4VC/IGLs/7l1XSZpLsyWdHxFM/Op7LODUfHs0zGStbdVMgI/Ewun/eL6qktbM7n0+yNZ+Rbc7rxIyIrbopl66tXSSdTLJk3BqYDryFJPIzKtnaWaSvlbQY2AmYBtxRI2aJpKds1SNiVY+26o1tOYpkIcOUKVNaZTMbAb14GJ0MvKByfDXJNr3qMT5aludb5XyrerzGM4iIC4ALIFm7tMpn+p+uhy4RMTkiBEwmOcJNIznC7dJFMb1YpI+FrbrZSOl5jB4Rf4uIK4G3kywPv1M53Yt9+VjYqptCGfbDaET8iWSruJOk5+Xkdhbpm5IMd9cC91ROjYWtuimUkXoF4EV5vy7vh2zJ922Sdw9gS2BBRKyppLeLaWer3m2MKZE6k+3APwGTm6RvwtMLRjdV0icCS+luwWhMbNVbbV4w6mtGzCJ9X+BMSfOBu7OIXgDsSXoYvR84snLzrJB0JDAHmCdpNvAQ8E7SNOIc4LKGG26sbNVNidS5G0jvnnwVuAVYRhpfLwduJr3rsm2LuN2BucDDwGrgVtI7MxPaXGt/4HrgEWBlvsasDvWblfOtzHHXA/vVadvQ5h69r+n4+doiPWOL9L7GFunGgIVuCsFCN0VgoZsisNBNEVjopggsdFMEFropAgvdFIGFborAQjdFYKGbIrDQTRFY6KYILHRTBBa6KQIL3RSBhW6KwEI3RWChmyKw0E0RWOimCCx0UwQWuikCC90UgYVuisBCN0VQS+iStpN0hKQrJd2VnaOXS7pR0uGSmpZj52gzXqjbo7+H5Mr8JuBXwNnAFaRf2f0mcLmkZ/zQo52jzbiizk/ukrw69wc2aUifDPyZ9GP8B1bS7RxtxpKOn2+tHj0iro2IqyLiyYb0+4Hz8+FelVN2jjbjipF4GH0i79dW0uwcbcYVwxJ6dpg7LB9WxdbWORpYTDLznVYzZgnJzWJ7SVvmaw/bOdqUw3B79DNID45zI+KnlfTRco6u5hu2c7SkoyQtlLRw6dKlrbKZjYCehS7peOAk4E7g0G7D836DOkdHxAURMT0ipk+aNKnLYk0/0ZPQJR0DnEMy0t07Ih5qyGLnaDOu6Frokk4AzgNuI4n8/ibZ7BxtxhVdCV3SR0kLPreQRP5Ai6x2jjbjizqT7ZEWYT5BGu8upIWvaCWvnaPNWNLx863lHJ3fNfk0sA64ATi+YcUfYDAiLs43j52jzfiizt1AcoeODtu8JnF2jjZjQcfP187RGTtH9zV2jjYGLHRTCBa6KQIL3RSBhW56YmBgAEk9bwMDA2NaX8+6ZDzrMnLstddeAMybN2+sLulZF2PAQjeFYKGbIrDQTRFY6KYILHRTBBa6KQIL3RSBhW6KwEI3RWChmyKw0E0RWOimCCx0UwQWuikCC90UgYVuisBCN0VQ6yfpzMbJjNn/Nyrl/u6BR0et/AUHv66nOPfopggs9FGg3/5CvgQ8dBkFBgYGWop1A/yFvKG+RfpMSedKukHSCkkh6dIOMbZHN+OGukOXU4FjgZ1JdiptsT26GW/UFfqJJG+hiazv+PwMJE0kiW4dsFdEHB4RHyHdJL8EZko6uCFmKnAWySxgekQcExEnAq8F7gZOkrRbQ8wMkive3cBrI+LEiDgG2DWXc1Yu15jaFunXRcQfo97Petke3Yw7RmPWxfboZtwxGrMube3RJS0GdiLZo99RI2aJpKfs0SNile3RNzx/mfN17r3iG23z/PKQnVue2/7A/2KHmW1HwSPKaAh9tOzRt8r5VvV4jfWQdBRwFMCUKVPaZTUN7DDz6DEV6nDZEAtG48IeHWyRXhKj0aP3ao/+vBzzYJuYMbVHH413NUbzPRDo/V2QjZ3R6NFtj27GHaMhdNujm3HHaAh9DrAMOFjS9KFESVsAn8mHX2+IuYhkp35sdZFH0jbAKfnw/IaYoeOP53xDMVOBY3J5F/XeDLMxUdci/QDggHw4Oe93k3Rx/veyiDgZbI8O/Tf1VgJ1H0Z3JtmQV5mWN4A/AScPnYiIH0raE/g4cCCwBXAXSZRfabbCGhHnShrM5RxG+ra5HTg1Ii5pVqmIOEnSItJ7OEcBTwK/Ac6MiB/XbNuI029TbyVQS+gRMQAMdFNwRNwEvKPLmKtIL3Z1E3MJ0PRGMGYI/+GFKQIL3RSBhW6KwEI3RWChmyKw0E0RWOimCCx0UwQWuikCC90UgYVuisBCN0VgoZsisNBNEVjopggsdFMEFropAgvdFIGFborAQjdFYKGbIrDQTRFY6KYILHRTBBa6KQIL3RSBhW6KwEI3RbBRCF3S9pK+LemvktZIGpR0dvV3003ZjIaH0ZgiaUdgAfB84H+AO4E3Ah8C9pW0e0Q080UyBbEx9OhfI4n8+Ig4ICI+FhFvBb5MMh44fYPWzowL+lrokqYBbwcGSbboVT4JrAQOzQa8pmD6Wug8bYH+s4h4snoiIh4BbiKZg715rCtmxhf9PkZvaa2e+SOpx3858IvGk1XnaOBRSb9vzNNv6JANXYPRpUX7ro6IZu6ET9HvQh+WVXpEXABcMNKVMuOPfh+6dKJXa3WzkdHvQu/Fjt0USL8LvaW1esZW6QYANbH87BvyYtFdpOnFHaszL5KeCywh3cyTImLlBqmkGRf0dY8eEXcDPwOmkmzRq3wK2Ar4Ti8ilzRT0rmSbpC0QlJIurRDzAxJcyU9JGmVpEWSTpA0oUnenSUNSLpJ0hJJj0u6T9L3Je1Ss47Pk3R/rtuNXbZvO0lHSLpS0l2SVktaLulGSYdLaqqNbtrYJFaSrsn1DUlNJ0MkvUbS9yr1uk/SdZIOalWvjkREX2/AjsDfSA+cPwQ+B1ybj38PbNdjubfkMh4B7sj/vrRN/ncBa4FHgW8BZ5JeRwjgB03y/28+t5C02PV54Kc57Qng32vU8YpcvwBu7LJ9H8xxfwW+l//fvg38PafPIX/j99rGJtc8LsevzjGbNsmzf27/GuAHwBnAN4EHc8yFPX2eG1qoIyT2HYCLSEOVx0mW7ecA2w6jzL1JY3wBe7UTOumh94H84UyvpG9Beg8ngIObfOgvbVLWf+T8y4DN2tTvsJzv6B6F/tYsqk0a0icDf85lHjicNjaU+wpgVRbuYBuh/y6f27NJvYY6tClFCn20txpC/0A+f0kLQQVwfRfX+0OO2bXF+Sm55/0madjWtdA7XP+UXOa5I9FG0nrNr7OIN+8g9NXA8hbl/Kjd/0u7ra/H6OOIoVcRrm5ybj6pJ5shafOa5T2R92sbT0gScDFpyvTD3VWzNs2uP5w2ngq8HpgVEWs6XPt3wERJb6kmSno+6a3UvwK3dyhjPSz0kaHlqwgRsRZYTOrVpnUqSNKbgFcB9wG3NclyAukb5vCIWNFjfdtdf1PSsAieKeqe2ijpDcDHgTMiYmGNKpwIrAB+LukySZ+TdCHpBngEOCAiVnfRJKD/XwEYLwzrVYQh8h+KfDcffjgi1jWcfxXwWeD8iPh5j3XtxBnAq4G5EfHTSnrXbZT0bFJ7bgc+XefiEXGDpN2Ay4H3Vk49QnoOu7VOOY24Rx8bOr6KkF8l/hHpAfgLEXF5w/lnkUSzBPjvUamkdDxwEmkm5dBuw/O+2sYvkHr4WRHxxPohTevwNuAG0jfarqQp4h1JzyOnA79oNS3Zlg39oNcPG50fRm+m/cPjbfn8K1uc3wqYl/N8sUWe04AnWX82Yioj8DBKWocI0hBh8nDbCOyZ63tak7yDNHkYBbYFHgLuBbZsEndljntft+1zjz4ytHwVIfc+LyE92N3T5PxzgZ+QhPGFiDipxTV2IfWa8yoLLkEaGwPsntP+3m3lJZ0AnEcS694RcX+TbN228fW5vp+q1jfX+cU5zxM5bed8PAPYBvhVRKxqUofr8n7X7lroMfpIcS1p/ntf4PsN5/Yg/fHH/GiYcZC0NemB783A6RFxaptrXEOaW2/kOcBBpDnmH5NmP2oj6aOkcfktwNsiotk1oPs23kZaVGrGQbne3yb10EN/0zs0YzOpRdxQ+uMtzrdmQw8L+mGj3oLRUrpbMNqGp4cD6329d1G3qfQ4dAE+wdOrs20X13ppY5uyBmk+dHkRaWpzHfD2hnM7kBasAnhHt23t65e6RhNJBwAH5MPJwD6kr+UbctqyiDi5If8c4DFgNmms+U7StNwc4L1R+c+WdB3pBrobaPUOzQ8j4pYO9ZxKGr7cFBFvaZe3IW4WaT5+HXAuzWdTBiPi4kpMV21sc+1B0vDlWZGmJqvnTiO9p/Qk6RvqTtL//7tJ3wJXRsS767bzKTZ0bzleN2CA1Hu02gabxOwOzAUeJq3w3UqaF57Qpldrt72vRj2n0tsrAJ3aF8C84bSxzbWH2r7eymg+/y7Sc8tS0rh/Belb4+hurlPd3KObIvCsiykCC90UgYVuisBCN0VgoZsisNBNEVjopggsdFMEFropgv8HIho3fk5SjUsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 144x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i, k in enumerate((\"holo\", \"control\")):\n",
    "    fig = plt.figure(figsize=(2, 4))\n",
    "    ax = fig.add_subplot(111)\n",
    "    its_conv_m = its_conv[k].squeeze().mean(axis=1)\n",
    "    its_conv_l, its_conv_h = np.percentile(its_conv[k].squeeze(), (2.5, 97.5), axis=1)\n",
    "    ax.bar([0, 1], its_conv_m, yerr=[its_conv_m - its_conv_l, its_conv_h - its_conv_m],\n",
    "           capsize=8, color=colors[i + 4])\n",
    "    ax.set_xticklabels([0, 1024, 2048], fontsize=20)\n",
    "    ax.set_ylim(0, 5100)\n",
    "    ax.set_yticks(np.arange(0, 6000, 1000))\n",
    "    ax.tick_params(axis=\"x\", length=0)\n",
    "    ax.tick_params(axis=\"both\", labelsize=20)\n",
    "    sns.despine(ax=ax)\n",
    "    plt.savefig(\"figs/conv-bar-{0}.pdf\".format(k), bbox_inches=\"tight\", transparent=True)\n",
    "    plt.savefig(\"figs/conv-bar-{0}.svg\".format(k), bbox_inches=\"tight\", transparent=True)"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}