notebook/md_compress_adaptive.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "\n",
    "from sklearn.decomposition import TruncatedSVD\n",
    "from sklearn.utils.extmath import svd_flip, randomized_svd\n",
    "\n",
    "from scipy.sparse.linalg import svds\n",
    "from scipy.misc import imfilter\n",
    "import scipy.ndimage as ndi\n",
    "\n",
    "import scipy.io\n",
    "\n",
    "from heapq import heappush, heappop\n",
    "\n",
    "from random import uniform\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "from pprint import pprint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def svd_wrapper(Y, k, method='svds'):\n",
    "    if method is 'svds':\n",
    "        Ut, St, Vt = svds(Y, k)\n",
    "        idx = np.argsort(St)[::-1]        \n",
    "        St = St[idx] # have issue with sorting zero singular values\n",
    "        Ut, Vt = svd_flip(Ut[:, idx], Vt[idx])\n",
    "    elif method is 'random':\n",
    "        Ut, St, Vt = randomized_svd(Y, k)\n",
    "    else:\n",
    "        Ut, St, Vt = np.linalg.svd(Y, full_matrices=False)\n",
    "        # now truncate it to k\n",
    "        Ut = Ut[:, :k]\n",
    "        St = np.diag(St[:k])\n",
    "        Vt = Vt[:k, :]\n",
    "        \n",
    "    return Ut, St, Vt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class WeightedReservoirSampler:\n",
    "    def __init__(self, reservoir_size):\n",
    "        self.m = reservoir_size\n",
    "\n",
    "        self.reservoir = []\n",
    "        \n",
    "    def add(self, weight, data):\n",
    "        rw = uniform(0.0,1.0)**(1.0/weight)\n",
    "        if(len(self.reservoir)< self.m):\n",
    "            heappush(self.reservoir, (rw,data))\n",
    "        elif rw > self.reservoir[0][0]:\n",
    "            heappop(self.reservoir)\n",
    "            heappush(self.reservoir, (rw,data))\n",
    "        #print(len(self.reservoir))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MDTrSampler:\n",
    "    def __init__(self, n_atoms, n_dim = 2, conv_size = 50, n_samples=1000, batch_size=100, manifold_size=64):\n",
    "        self.n_dim = n_dim\n",
    "        self.conv_size = conv_size\n",
    "        self.n_atoms = n_atoms\n",
    "        self.batch_size = batch_size\n",
    "        self.l = manifold_size\n",
    "        self.n_samples = n_samples\n",
    "        self.bq = np.zeros((batch_size,n_atoms,3)) # batch queue\n",
    "        self.total_samples = 0\n",
    "        self.bq_index = 0\n",
    "        self.Btp = np.zeros((n_atoms*3,self.l))\n",
    "        self.strm_smplr = WeightedReservoirSampler(n_samples)\n",
    "        self.last_Vt = None\n",
    "    \n",
    "    def traj_char(self, c, s):\n",
    "        x = c.shape[0]\n",
    "        # scaled data <-- not so much meaningful as we're doing l2 normalization later\n",
    "        ps = np.mat(c[:,:self.n_dim]) * np.mat(np.diag(s[:self.n_dim]))\n",
    "        # gradient or changes of each data point\n",
    "        psd = ps[0:(x-1),:self.n_dim] - ps[1:x, :self.n_dim];\n",
    "        \n",
    "        # convoluted gradient --> smoothed gradient\n",
    "        psdm = ndi.convolve(psd, np.ones((self.conv_size,2))/self.conv_size*2);\n",
    "        # L2 normalized smoothed gradients --> now we focus on the angle of gradient only as atoms may move different speed \n",
    "        npsdm = np.divide(psdm, np.mat(np.sum(np.abs(psdm)**2,axis=-1)**(1./2)).T  * np.mat(np.ones((1,n_dim))))\n",
    "        # Angle smoothing with the assumption that the overall angle can not be radically changed\n",
    "        nmpsdm = ndi.convolve(npsdm, np.ones((self.conv_size,2))/self.conv_size*2);\n",
    "        # Smoothed angle changes (accelerations)\n",
    "        psdd = nmpsdm[0:(x-2)] - nmpsdm[1:(x-1)]\n",
    "        # results for gradient of changes\n",
    "        psdu = np.sum(abs(psd),axis=1)\n",
    "        # results of smoothed normalized acceleration\n",
    "        psddu = np.sum(abs(psdd),axis=1)\n",
    "        prob_dist = (abs(psddu) / np.sum(psddu))\n",
    "        return nmpsdm, psddu, prob_dist\n",
    "    \n",
    "    def strmML(self, t2):\n",
    "        n_t = t2.shape[1]\n",
    "        Ct = np.concatenate( (self.Btp, t2), axis = 1)\n",
    "        Ct = Ct[:, ~(Ct==0).all(0)]\n",
    "\n",
    "        Ut, St, Vt = svd_wrapper(Ct, self.l) #SVD_l(matrix)\n",
    "        Ut_l = Ut[:, :self.l]\n",
    "        St_l = St[:self.l] - St[self.l-1] # to be adaptive, added singular substraction\n",
    "        Vt_l = Vt[:self.l, -n_t:]\n",
    "        self.Btp= np.dot(Ut_l, np.diag(St_l))\n",
    "\n",
    "        return Ut_l, St_l, Vt_l\n",
    "    \n",
    "    def batch_sampling(self, trace):\n",
    "        (x,y,z) = trace.shape\n",
    "        self.sampling_rate = float(self.n_samples) / float(x)\n",
    "        t2 = trace.reshape((x, y*z),order='C')\n",
    "        c,s,v = svd_wrapper(t2, 2, method='random')\n",
    "        nmpsdm, psddu, prob_dist = self.traj_char(c,s) \n",
    "        total = 10**-10 + prob_dist[0]\n",
    "        sampling_entries = int(x * self.sampling_rate)\n",
    "        target = np.zeros((sampling_entries, y, z))\n",
    "        time_stamps = np.zeros(sampling_entries)\n",
    "        #output initialization: the first one should be added always\n",
    "        target[0,:,:] = trace[0,:,:]\n",
    "        time_stamps[0] = 0;\n",
    "        t_idx = 1;\n",
    "        for i in range(2,x):\n",
    "            total = total + prob_dist[i-2];\n",
    "            if(total > 1 / float(sampling_entries - 1)):\n",
    "                target[t_idx,:,:] = trace[i,:,:]\n",
    "                time_stamps[t_idx] = i\n",
    "                total = total - 1 / float(sampling_entries-1);\n",
    "                t_idx = t_idx + 1;\n",
    "        return target, time_stamps\n",
    "    \n",
    "    \n",
    "    def adaptive_sampling_update(self):\n",
    "        t2 = self.bq[0:self.bq_index,:,:].reshape((self.bq_index, self.n_atoms*z),order='C').T\n",
    "        Ut, s, VTt = self.strmML(t2)\n",
    "        if self.last_Vt is None: # for the first time\n",
    "            c = VTt.T\n",
    "            self.strm_smplr.add(10**10,[0, self.bq[0,:,:]])\n",
    "            nmpsdm, psddu, prob_dist = self.traj_char(c, s)\n",
    "            sidx = self.total_samples - self.bq_index\n",
    "            for i in range(2,self.bq_index):\n",
    "                self.strm_smplr.add(psddu[i-2], [sidx+i, self.bq[i,:,:]])\n",
    "        else:\n",
    "            c =  np.concatenate( (self.last_Vt, VTt.T), axis = 0)\n",
    "            nmpsdm, psddu, prob_dist = self.traj_char(c, s)\n",
    "            sidx = self.total_samples - self.bq_index\n",
    "            for i in range(0,self.bq_index):\n",
    "                self.strm_smplr.add(psddu[i], [sidx+i, self.bq[i,:,:]])\n",
    "                \n",
    "        self.last_Vt = VTt.T[-2:,:] # for last two to start compare against\n",
    "        \n",
    "        \n",
    "\n",
    "    def adaptive_sampling_step(self, dataframe_t):\n",
    "        self.total_samples = self.total_samples + 1\n",
    "        self.bq[self.bq_index] = dataframe_t\n",
    "        self.bq_index = self.bq_index+1\n",
    "        if self.bq_index == (self.batch_size):\n",
    "            self.adaptive_sampling_update()\n",
    "            self.bq_index = 0\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat= scipy.io.loadmat('../data/collision_0-5k.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = mat['trace']\n",
    "(x,y,z) = trace.shape\n",
    "t2 = trace.reshape((x, y*z),order='C')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ut, St, VTt = svd_wrapper(t2, 2, method='random')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = Ut\n",
    "s = St"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5000, 2)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_dim=2\n",
    "step=1000\n",
    "c.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X901PWd7/HnW6JQaQWCoOBAQRNdoALWJHgt7hUxSt0Q/MFKXLcStcWtukJTW/HuhS7QCnqEwln1dj03INuqoXJbAUURdN0DreWXghZZSxRaAlSRAIoKCLzvH58JmYSEJGSS+ZHX45w5853vfGbmnRyYV77fz4+vuTsiIiJVTkt0ASIiklwUDCIiUoOCQUREalAwiIhIDQoGERGpQcEgIiI1pGwwmNkIM3vPzMrNbGKi6xERSReWivMYzKwd8CcgH6gA1gK3uPu7CS1MRCQNpOoRQx5Q7u4fuPthoAwYleCaRETSQqoGw3nA9pjHFdF9IiLSTKkaDFbHvtQ7JyYikoRSNRgqgF4xjyPAzgTVIiKSVlI1GNYC2WbW18zOAIqAxQmuSUQkLaToqKSzHbrQvv123J2zzz6bHj16JLosEZGkdfbZZ7Ns2bJl7j6iobYZrVFQ/PUB1lFeDpFIomsREUkNZnZ2Y9ql6qkkAHr1ariNiIg0TUoGw8UXV2+3a5e4OkRE0lFKBsMZZ0Bubtg+dgyKihJbj4hIOknJYABYswZOPz1s/+Y3ia1FRCSdpGwwAMyYEe6//BK+9rXE1iIiki5SOhhKSuCOO8L2gQPQtWti6xERSQcpHQwApaXV4VBZCf37J7YeEZFUl/LBACEcqjqjN2+GyZMTW4+ISCpLi2CA0BmdmRm258xJbC0iIqksbYIB4OWXw/0nn8AVVyS2FhGRVJVWwZCbC0OHhu1Vq+CFFxJbj4hIKkqrYABYuRK6dQvb118PFRWJrUdEJNWkXTAAfPRR6G84ehSGDUt0NSIiqSUtgwHgnnvCfXm5lswQEWmKtA2GqVOrh7AuWAATJiS2HhGRVJG2wQBhCOugQWF7zhx46qmEliMikhLSOhgANmyArKywPW6cOqNFRBqS9sEAsGULnHtuWGzvm99MdDUiIsmtTQQDwK5dYRjr7t0weHCiqxERSV5tJhgAJk4M9xs3wqxZia1FRCRZtalgKCmp7ox+8EH1N4iI1KVNBQOEzuhIBA4f1uQ3EZG6xCUYzGyEmb1nZuVmNrGO59ub2YLo86vNrE90f56ZbYjeNprZDfGopyFvvAEdOoTJb1psT0SkpmYHg5m1Ax4Hvg30B24xs9qXy7kT2OvuWcDPgYej+/8I5Lj7YGAE8O9mltHcmhoSicDPfha2V62CO+9s6U8UEUkd8ThiyAPK3f0Ddz8MlAGjarUZBcyPbi8EhpuZufvn7n4kur8D4HGop1FKSmDkyLA9d64mv4mIVIlHMJwHbI95XBHdV2ebaBDsB7oCmNkQM9sEvAP8U0xQ1GBm48xsnZmt2717dxzKhsWLq5fN+P731RktIgLxCQarY1/tv/zrbePuq919AJALPGhmHer6EHd/0t1z3D2nW9W62nGwZk2YGX3wIFxzTdzeVkQkZcUjGCqAXjGPI8DO+tpE+xA6AZWxDdx9M/AZ8I041NQkzzwT7jdvhvz81v50EZHkEo9gWAtkm1lfMzsDKAIW12qzGBgb3R4NvObuHn1NBoCZfR24CNgWh5qaJDcXZs4M2ytWqDNaRNq2ZgdDtE/gXmAZsBn4tbtvMrOpZlYYbVYKdDWzcqAEqBrSOhTYaGYbgN8Cd7v7x82t6VSUlMAdd4TtuXNh8uREVCEiknjm3moDgeImJyfH161b1yLvXVgIS5aE7fHjYfbsFvkYEZFWZ2br3T2noXZtbuZzQxYvhjFjwvacOTqtJCJtj4KhDmVlNU8racE9EWlLFAz1KC2tngA3ebLmOIhI26FgOInFi8NqrJ99Bvffn+hqRERah4KhAYXRcVVLl8LatYmtRUSkNSgYGjB1apjn8OmnYfLbCy8kuiIRkZalYGiENWtg6FDYvz8cQejIQUTSmYKhkVauDEcO7jBtWqKrERFpOQqGJigqCvcvvqijBhFJXwqGJigpCaeUjh2De+5JdDUiIi1DwdBEzz4LX/taOGLQrGgRSUcKhiaKRGDChLD91FM6pSQi6UfBcAqmTq0+pfS97yW6GhGR+FIwnKKVK8Os6I0bdXEfEUkvCoZmqDqltGKFFtoTkfShYGiG4uLqJbqnTVN/g4ikBwVDM5WVwdVXw759MHy4VmEVkdSnYIiD5cur11N69NFEVyMi0jwKhji5/PJw//rrCS1DRKTZFAxxMnt29SilqqUzRERSkYIhjqpGKS1cqI5oEUldCoY4Ki4OHdFHj2oFVhFJXQqGOLv11nC/ZIku6iMiqUnBEGfFxTByZNh++OGEliIickoUDC3giSegUydYtUod0SKSehQMLSASgcmTw/aCBTqlJCKpRcHQQkpKqk8p/epXia1FRKQpFAwt6MYbw72OGkQklSgYWlBsR/ScOQktRUSk0RQMLWzSJMjI0NLcIpI6FAwtLDcXbrstbGsdJRFJBQqGVnDDDeF+yZJwnWgRkWSmYGgFBQXVF/QpLU1sLSIiDVEwtJKyMhg6VJPeRCT5KRha0Z13hnsNXxWRZKZgaEVaR0lEUkFcgsHMRpjZe2ZWbmYT63i+vZktiD6/2sz6RPfnm9l6M3snen9VPOpJZrHrKBUWJroaEZETNTsYzKwd8DjwbaA/cIuZ9a/V7E5gr7tnAT8Hqv5e/hgY6e4XA2OBXza3nmQXiYQlMk4/XaOURCQ5xeOIIQ8od/cP3P0wUAaMqtVmFDA/ur0QGG5m5u5vufvO6P5NQAczax+HmpJaQQF85zthW6OURCTZxCMYzgO2xzyuiO6rs427HwH2A11rtbkJeMvdD9X1IWY2zszWmdm63bt3x6HsxCotrR6llJeX6GpERKrFIxisjn3elDZmNoBweumu+j7E3Z909xx3z+nWrdspFZpsVq4MM6PXroUrrkh0NSIiQTyCoQLoFfM4Auysr42ZZQCdgMro4wjwW+A2d38/DvWklN/8prozumo4q4hIIsUjGNYC2WbW18zOAIqAxbXaLCZ0LgOMBl5zdzezzsCLwIPu/rs41JJyIhGYPTtsz51bfYEfEZFEaXYwRPsM7gWWAZuBX7v7JjObamZVAzJLga5mVg6UAFVDWu8FsoBJZrYheuve3JpSTXExjB8ftqdN00glEUksc6/dHZD8cnJyfN26dYkuI+6KisKs6M6d4Z13wtGEiEi8mNl6d89pqJ1mPieRsjK4+mrYtw/uvjvR1YhIW6VgSDIPPVQ9+W3ChERXIyJtkYIhyeTmwsRoD8ycOVpsT0Ran4IhCU2dqmtFi0jiKBiS1BNPQJcu4VrRWmxPRFqTgiFJRSKwbBl07Rr6G/LzE12RiLQVCoYklpsLL70Uhq/qyEFEWouCIcnl5sIvf1k9UmnWrERXJCLpTsGQAgoKquc1TJsGFRWJrUdE0puCIUXMnl09+e2WWxJdjYikMwVDClm+vPoaDlqmW0RaioIhxTz7bPUy3RqpJCItQcGQYiKRcORQNcdB4SAi8aZgSEG5ufD226HPQeEgIvGmYEhRkQjMm1d95KCrv4lIvCgYUlgkUj2vYd48hYOIxIeCIcUVF8PMmWEC3Ny5CgcRaT4FQxooKYHp08P23LlaqltEmkfBkCZKSmDMmLCtpbpFpDkUDGnk0Uc1jFVEmk/BkEYikROHsa5dm+iqRCTVKBjSTNUEuKpwGDpUK7KKSNMoGNLU8uVwxx1w+DD88IfqkBaRxlMwpLHS0uoO6YcfTmwtIpI6FAxprqysekXWCRMSXY2IpAIFQxswbFi4f+opXeRHRBqmYGgDpk4NRw3798Pttye6GhFJdgqGNuLZZ6vnOBQW6shBROqnYGgjquY4jBwJS5bAgAEKBxGpm4KhDYlEYPHicFrpk0907WgRqZuCoQ2KvTxoUVGiqxGRZKNgaIMiEfjVr8JS3QsW1L1U9x133EH37t35xje+cXxfZWUl+fn5ZGdnk5+fz969ewFwd+677z6ysrIYOHAgb7755vHXzJ8/n+zsbLKzs5k/f36L/2wi0nwKhjaqoABmzAjbc+eGoayxiouLefnll2vsmzFjBsOHD2fLli0MHz6cGdE3eOmll9iyZQtbtmzhySef5Pvf/z4QgmTKlCmsXr2aNWvWMGXKlONhIiLJS8HQhsUu1f300zWf+9u//VsyMzNr7Fu0aBFjx44FYOzYsTz//PPH9992222YGZdddhn79u1j165dLFu2jPz8fDIzM+nSpQv5+fknhI2IJB8FQxsXu1R3Q1d/+/DDD+nRowcAPXr04KOPPgJgx44d9OrV63i7SCTCjh076t0vIsktLsFgZiPM7D0zKzeziXU8397MFkSfX21mfaL7u5rZf5rZATN7LB61SNPEXjd67lyYPLnp7+HuJ+wzs3r3i0hya3YwmFk74HHg20B/4BYz61+r2Z3AXnfPAn4OVC3pdhCYBNzf3Drk1BUXw/jxYXvatPpXYj3nnHPYtWsXALt27aJ79+5AOBLYvn378XYVFRX07Nmz3v0iktziccSQB5S7+wfufhgoA0bVajMKqBqSshAYbmbm7p+5+ypCQEgCzZ4dJr9B/ZcGLSwsPD6yaP78+YwaNer4/v/4j//A3fnDH/5Ap06d6NGjB9deey2vvPIKe/fuZe/evbzyyitce+21rfHjiEgzZMThPc4Dtsc8rgCG1NfG3Y+Y2X6gK/BxYz/EzMYB4wB69+7dnHqlHk88EeY2rFgB5513C8eOvc7HH39MJBJhypQpTJw4kZtvvpnS0lJ69+7Nc889B8B1113H0qVLycrK4swzz2TevHkAZGZmMmnSJHJzcwGYPHnyCR3aIpJ8rK7zwE16A7O/B6519+9GH38HyHP3f45psynapiL6+P1omz3Rx8VAjrvf25jPzMnJ8XXr1jWrbqnbCy/ADTfAkSPhQj+lpYmuSETixczWu3tOQ+3icSqpAugV8zgC7KyvjZllAJ2Ayjh8tsRZQUH1RX3mztWV30TaongEw1og28z6mtkZQBGwuFabxcDY6PZo4DVv7qGKtJjY+Q319TeISPpqdjC4+xHgXmAZsBn4tbtvMrOpZlYYbVYKdDWzcqAEOD6k1cy2AbOAYjOrqGNEkyTAD38IGRmNm98gIuml2X0MiaA+htYxa1YICICZM8ORhIikrtbsY5A0VVJSPb9h4sQT11MSkfSkYJCTmj079Dd8+SV873vqjBZpCxQM0qCysjD57cgRuPVWXflNJN0pGKRRYq/8dvfdia5GRFqSgkEa7YEHwv2SJepvEElnCgZptIKC6vkNJSWwdm1i6xGRlqFgkCYpK4Orr4a9e+GaaxQOIulIwSBNtnx5CId9+yA/X53RIulGwSCnZPny0Bm9fz/cfnuiqxGReFIwyCmbNUvLZoikIwWDnLLc3JorsZ7KZUFFJPkoGKRZYpfNmDat+vrRIpK6FAzSbLNnh4v6gNZUEkkHCgaJi9LS6jWVvvtdrakkksoUDBI3VXMcjh6Ff/xHDWMVSVUKBomr2sNYNQFOJPUoGCTuVq4MRw4rVsBllykcRFKNgkFaRNWRw7FjcM89ia5GRJpCwSAtZuXKMNdh7VrIy0t0NSLSWAoGaVFr1tQMB3VIiyQ/BYO0uNhwGDhQQ1lFkp2CQVrFmjXVy3WPGqUOaZFkpmCQVqMOaZHUoGCQVvXss/C1r4UjhsLCRFcjInVRMEirikTgmWfg9NPDtaMVDiLJR8Egra6gAJ58MlzLYckSKCpKdEUiEkvBIAlRXFx9LYcFC7Qiq0gyUTBIwpSUhBVZq7YVDiLJQcEgCVW1IuvevWHRvQkTEl2RiCgYJOGWL6++0M+cObp+tEiiKRgkKZSWwqRJYXvuXIWDSCIpGCRpTJ16YjhobSWR1qdgkKQydSrMnAlnnBHC4dJL1Skt0toUDJJ0Skpg1aowYumjj9QpLdLaFAySlHJzw4gldUqLtD4FgyS12p3SmiUt0vLiEgxmNsLM3jOzcjObWMfz7c1sQfT51WbWJ+a5B6P73zOza+NRj6SX2E7pBQu0vpJIS2t2MJhZO+Bx4NtAf+AWM+tfq9mdwF53zwJ+DjwcfW1/oAgYAIwAnoi+n0gNseGgxfdEWlY8jhjygHJ3/8DdDwNlwKhabUYB86PbC4HhZmbR/WXufsjdtwLl0fcTOUHViKWqxfeuuCLRFYmkp3gEw3nA9pjHFdF9dbZx9yPAfqBrI18LgJmNM7N1ZrZu9+7dcShbUlFJCfz2t3DWWWHk0hVXaK6DSLzFIxisjn3eyDaNeW3Y6f6ku+e4e063bt2aWKKkk4ICePppaNcuhMPFF+tSoSLxFI9gqAB6xTyOADvra2NmGUAnoLKRrxU5QUEBvPFGWIBv3z4YNkxHDiLxEo9gWAtkm1lfMzuD0Jm8uFabxcDY6PZo4DV39+j+ouiopb5ANrAmDjVJG5CbGxbgGzQIPvsMbrxRRw4i8dDsYIj2GdwLLAM2A792901mNtXMqsaOlAJdzawcKAEmRl+7Cfg18C7wMnCPux9tbk3StrzwQvV1pC+/HGbNSnRFIqnNwh/uqSUnJ8fXrVuX6DIkiVRUwN13h9FKEO4LChJbk0iyMbP17p7TUDvNfJa0EInA4sUwdGh4/PDD4UhCRJpOwSBpZeXKEA6rVsHIkVpfSeRUKBgk7axcWb343ty5miUt0lQKBklLsYvvLVmixfdEmkLBIGlr6tTqI4cFCxQOIo2lYJC0VlpaMxx0WkmkYQoGSXu1TysVFmoinCSHo0ePcskll1AQHVu9detWhgwZQnZ2NmPGjOHw4cMAHDp0iDFjxpCVlcWQIUPYtm3b8feYPn06WVlZXHTRRSxbtiwudSkYpE2ovTLr0KG6XKgk3pw5c+jXr9/xxw888AA/+MEP2LJlC126dKG0tBSA0tJSunTpQnl5OT/4wQ944IEHAHj33XcpKytj06ZNvPzyy9x9990cPdr8OcIKBmkzSkrg978Pp5YOHw6XC83PT3RV0lZVVFTw4osv8t3vfhcAd+e1115j9OjRAIwdO5bnn38egEWLFlFePhYz+Id/GM1zz72KmTNgwCLefruIs89uzw039CUrK4s1a5q/qlBGs99BJIXk5obbeefBtGmwYgXk5UEc/i+JNMmECRN45JFH+PTTTwHYs2cPnTt3JiMjg8JCePXVCF98sYMOHeDQoR1UrTfqXrUO6R5gB3AZBw7Axo0AEV55ZQfnngsffghdusCBA3DoUNNq0xGDtElVp5ZOOy30NwwerH4HaT0vvPAC3bt359JLLwXgT3+C/HznL38Bs3C68/PPwd2iX+rVSxeddlpoE65aUNeSRsZf/wruUFkZjo7NoEOHxtenYJA2q6QE/vCHsDrrxo3wrW9pprS0jilTfscvfrGY007rw1VXFbFly2ts2DCBL7/cBxyJtqrArCft20NGRoTOnbfjDocOHSEzcz/HjmXy0EMRHnpoO9u3w7nnhtdAT849tyo8qoWAuWRwY+pTMEiblpsLGzaE5TO+/DLMlL7iCq2zJC2nfXtYt2467hW4b8O9DLgKeJqOHYeRkbEQgLvums9jj43i4EGYPbuQMWPC1ZEXLlzIVVddhZlRWFhIWVkZ3bod4ve/30rfvls4ciSPXbvg2LFw1NCuHXTrFj4X3trQmBoVDCKEBfiqRitVrbOkCXESD/n54TTO4MHhr/joCFQgXKL2nHPCaDl3ePvth7nkkllccEEWe/bs4c7oIeydd97Jnj17yMrKYtasWcyYMQOAAQMGcPPNN9O/f39GjBjB448/Trt27Wp8/pEj8NFHcPBg42vWstsitRQWVi/fnZurjmk5Nfn5YXBDfRLx1atlt0VO0eLFYUJcVcd0Xh489VSiq5JUUlcoDBoUbvPmJSYUmkJHDCL1WLsW7rmnerTSmDFQVpbYmiQ1nHZazS//ZPma1RGDSDNVnUa6+urweMGCcPQg0pDYIGjKMNFkoWAQacDy5TB+fOg4XLsW+vfXnAepX+1hol98kZg6mkPBINIIs2eHvoevfAU2b4bLL9ecBzlR1641HyfLKaSmUjCINFJBQZihevXVYQig5jxIbZWV1dupGgqgYBBpkkgknFqaNClMHKqa86DrPEg6UTCInIKpU+GNN0IoQJj3MHiwjh7asti+hXnzEldHPCgYRE5Rbm7od5g5Ezp2DOstjRwZTi9J2zN+fPV2qvc/KRhEmqmkBP77v6uPHlatgr59NSmurZkzp3o7DtfKSSgFg0gcRCLV6y316QPbtsHtt4ejh4qKRFcnLS2dTiOBgkEkrgoKYOvW6klxq1bB3/yNLiOazmJDYdAgKC5OWClxo2AQaQHLl4ejh0GD4LPPwmmG/v3VOZ1uai1kyoZGLWqd/BQMIi2koCB8UUyaBP36hYlxI0eGjknNnE59ZuGaB1VSed5CbQoGkRY2dSq8+24IBwgT4y6/XHMfUlns6aOsrPQKBVAwiLSad98NRw8jR4aZ01VzH9T/kJwOHjxIXl4egwYNYsCAAfzkJz8BwGwrMATIJjt7DJs2hSvvHDp0iDFjxpCVlcWQIUPYtm3b8feaPn06WVlZXHTRRSxbtqz1f5imcveUu1166aUukspmznQfNMg9/K3pnpXlvmZNoquSWMeOHfNPP/3U3d0PHz7seXl5Dm84/L3Ds96nj/tdd93lTzzxhLu7P/74437XXXe5u/uzzz7rN998s7u7b9q0yQcOHOgHDx70Dz74wM8//3w/cuRIQn4mYJ034jtWRwwiCVBSEvofxowJj8vL4bLLwuVE1UGdHMyMr371qwB8+eWXrFnzJWDAa3zlK6PZuhXGjh3L888/D8CiRYsYO3YsAKNHj+bVV1/F3Vm0aBFFRUW0b9+evn37kpWVxZokvyyggkEkgcrKwiml3NzQkblggWZPJ5OjR4/SpctgOnbsDuQDF9C+fWc+/zwDgEgkwo4dOwDYsWMHvXr1AiAjI4NOnTqxZ8+eGvtrvyZZKRhEEqygIFwQaN68EBBQPXt68uTE1tbWXXllO/bt2wBUAGu4/PLNRCI121i0J9rr6IE2s3r3J7NmBYOZZZrZcjPbEr3vUk+7sdE2W8xsbMz+n5nZdjM70Jw6RNJBcXEIiPHjoX37MHt62rRw9DBrlpbYaG0TJoSADjpzySVXcv31f2Dfvn0cOXIEgIqKCnr27AmEI4Ht27cDcOTIEfbv309mZmaN/bVfk6yae8QwEXjV3bOBV6OPazCzTOAnhG78POAnMQGyJLpPRKJmzw59DmPGQKdO4cvphz+sXmJDAdHyunaFOXN2A/sAuOeeL+jYcQX9+vVj2LBhLFy4EID58+czatQoAAoLC5k/fz4ACxcu5KqrrsLMKCwspKysjEOHDrF161a2bNlCXrJfI7YxPdT13YD3gB7R7R7Ae3W0uQX495jH/w7cUqvNgaZ8rkYlSVuxfbv7pEnuY8ZUj2CCMKJp5sxEV5d+hg51z8io+j1vdBjsXbte7AMGDPApU6a4u/v777/vubm5fsEFF/jo0aP94MGD7u7+xRdf+OjRo/2CCy7w3Nxcf//994+/709/+lM///zz/cILL/SlS5cm5Gdzb/yoJPNmzMwws33u3jnm8V5371Krzf1AB3f/afTxJOALd380ps0Bd/9qA581DhgH0Lt370v//Oc/n3LdIqnoqafg6adhxYrqfUOHwrBhYRKdNE+7djVnMnfokJrXaz4ZM1vv7jkNtWvwVJKZrTCzP9ZxG9XYWurY1+Q0cvcn3T3H3XO6devW1JeLpLzi4pprMEE4zTRtWpgop47qpquogLy8mstbnH463HFH+oVCUzQYDO5+tbt/o47bIuBDM+sBEL3/qI63qAB6xTyOADvjUbxIW1S1BtO8eeGIAcJFgqZNCyOZZs1KbH2p4ooroFevmutWZWXB4cNQWpq4upJBczufFwNVo4zGAovqaLMMuMbMukQ7na+J7hORZiguhpUrwxXkqpb53rYtdFT36NH8gNi+fTvDhg2jX79+DBgwgDnRK9FUVlaSn59PdnY2+fn57N27Fwj9lffddx9ZWVkMHDiQN9988/h7zZ8/n+zsbLKzs4930CZCURFkZ0P37rEjjiAzM/QqbNmSsNKSS2M6Iuq7AV0Jo5G2RO8zo/tzgP8b0+4OoDx6uz1m/yOEI4pj0ft/bcznqvNZ5ERLlrjn5p7YST1ypPu8eU1/v507d/r69evd3f2TTz7x7Oxs37Rpk//oRz/y6dOnu7v79OnT/cc//rG7u7/44os+YsQIP3bsmL/xxhuel5fn7u579uzxvn37+p49e7yystL79u3rlZWVcfmZG2PMmNCp3Llzzd8NuHfo0GplJAUa2fncausbxfOmYBCp35o17v36uZ95Zs0vwX793MePDwFyKgoLC/2VV17xCy+80Hfu3OnuITwuvPBCd3cfN26cP/PMM8fbV7V75plnfNy4ccf3127XErZvDz/vWWedGAaZmSEs2uLaVI0NBs18FkkzublhJdf33gtzIapmU2/eHC4YNHJk6Ky+887GX3Z027ZtvPXWWwwZMoQPP/yQHj16ANCjRw8++ih0Lda39ENrLQnxwguhI7lv39B3sHkzfPJJeK5z5/B7mDkT9uwJS5FU/V7kRBmJLkBEWkYkEr4AIXxpzpkD+/eHztaNG8Pt6afh+utDm7y8sLhfbQcOHOCmm25i9uzZnHXWWfV+ntcx9L01loQoKoKXXqoOgSrt24fLqvbuHa7HLY2nYBBpAwoKwg3CsNY33gizq7dtCwv3Qbj/t3+DUaPCEUVxcVhV9KabbuLWW2/lxhtvBOCcc85h165d9OjRg127dtG9e3eAepd+iEQivP766zX2X3nllc3+mfLza87pAMjIgD594NxzQ8e8nJpmTXBLlJycHF+3bl2iyxBJeUVF0LFjmBuxe3fN577+defo0bGYZfKXv8w+vv9HP/oRXbt2ZeLEicyYMYPKykoeeeQRXnzxRR577DGWLl2vbSHaAAAGEklEQVTK6tWrue+++1izZg2VlZVceumlx0cpffOb32T9+vVkZmY2us6KCrj/fli9Olzk6NChmvW2bx9GYm3d2qxfR9pr7AQ3HTGItGFVp5og/AV+8GD40l27Fv78598BvwQu5rTTBnPmmdCz50NkZk6kY8ebKS0tpXfv3jz33HMAXHfddSxdupSsrCzOPPNM5s2bB0BmZiaTJk0iN3pSf/LkyY0OhdNOg4EDw2mv2jIywqS0f/7nsL6UxI+OGETkBJMnw8svh0X8ap+ugTDu//rr4bPPYP16+Lu/g1/+MoRLfX0VjdG9e/jrf/fuEFCxzjgDevYMRwzf+hY8+ignLIEtJ9fYIwYFg4icVEVFWI/p3HNDv8Rf/9rwawYNCh3d3bqFcNm9O2xDOCopL4dzzglf9nv3hraVlXD06Inv1a1bWMdo1674/lxtUVoHg5ntBlpiFb2zgY9b4H1bkmpuHar5uEiP8Hf9sWPhZM9XOsKn+6FzJhw+BGeedEHM+u0Guh6DI1+GE0UO7P0Y/tzIQbUJkUr/Lj4GcPcRDTVMyWBoKWa2rjFpmkxUc+tQzS0v1eqF1Ky5MTTBTUREalAwiIhIDQqGmp5MdAGnQDW3DtXc8lKtXkjNmhukPgYREalBRwwiIlKDgkFERGpI+WAwsxFm9p6ZlZvZxDqeb29mC6LPrzazPtH9Xc3sP83sgJk9Vs97LzazP8Y8/nsz22Rmx8wsp1bbB6Of8Z6ZXZvsNTfmvZKw5nwzW29m70Tvr0qBmvPMbEP0ttHMbkj2mmOe7x19v/uTuV4z62NmX8T8nn9RX73JUnP0uYFm9kb0+XfMrMPJ6m5VjbloQ7LegHbA+8D5wBnARqB/rTZ3A7+IbhcBC6LbHYGhwD8Bj9Xx3jcCzwB/jNnXD7gIeB3IidnfP/rZ7YG+0ZraJXnNJ32vJK35EqBndPsbwI4UqPlMICO6XXVd9Ixkrjnm+f8HPAfcn8z1An1i26XIv+UM4G1gUPRxV+r5zkjELdWPGPKAcnf/wN0PA2XAqFptRgFVF5ldCAw3M3P3z9x9FXCw9pua2VeBEuCnsfvdfbO7v1dHHaOAMnc/5O5bCZcwzUvmmk/2Xklc81vuvjP6cBPQwczaJ3nNn7v7kejDDoTpvPVJipqjr7ke+IDwe076epsgWWq+Bnjb3TdG2+1x9zoWBEmMVA+G84DtMY8rovvqbBP9D7qfkM4nMw2YCXwexzqa0rY1am6KZKz5JuAtdz9Uz/NJU7OZDTGzTcA7wD/FBEVS1mxmHYEHgCkNNE2KeqP6mtlbZvZfZnZFCtR8IeBmtszM3jSzHzfyda0i1YOhrstA1f6LrDFtqhubDQay3P23ca6jKW1bo+amSKqazWwA8DBw18maNaKeVqnZ3Ve7+wAgF3jwJOeSk6XmKcDP3f1AA+2Spd5dQG93v4TwV/szZlbfpeaSpeYMwmmpW6P3N5jZ8Ca8vkWlejBUAL1iHkeAnfW1MbMMoBNQeZL3/B/ApWa2DVgFXGhmr8ehjmSruSmSpmYziwC/BW5z9/dToeYq7r4Z+IzQP5LMNQ8BHom+ZgLwv8zs3mStN3oKd090ez2hD+HCeponRc3Rz/gvd//Y3T8HlgLfbOA1rSbVg2EtkG1mfc3sDEJHUe2ruy4Gxka3RwOvebS3py7u/n/cvae79yEk+Z/c/coG6lgMFEVHM/QFsoE1SV5zUyRFzWbWGXgReNDdf5ciNfeNfrlgZl8ndERuS+aa3f0Kd+8Tfc1s4CF3r2sUTlLUa2bdzKxddPt8wv+/D+ppnhQ1A8uAgWZ2ZvTfx/8E3m3gNa3Hk6AHvDk34DrgT4S/Ev4lum8qUBjd7kAYWVFO+LI+P+a12wh/CRwgJHjt0Ql9qDnC4IZou0PAh8CymOf+JVrDe8C3U6Tmk75XstUM/G/CX9wbYm7dk7zm7xA6cDcAbwLXp8K/jZg2/0o9o5KSpV5Cf9MmwgijN4GRqfA7Bv4xWvcfgUdO5fuvpW5aEkNERGpI9VNJIiISZwoGERGpQcEgIiI1KBhERKQGBYOIiNSgYBARkRoUDCIiUsP/B/MZFikpN0KaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(c[:,0],c[:,1], marker=\".\", c='b', s=0.5)\n",
    "ax.set_xlim(min(c[:,0]), max(c[:,0]))\n",
    "ax.set_ylim(min(c[:,1]), max(c[:,1]))\n",
    "for i in range(0,x,step):\n",
    "    ax.annotate(str(i),(c[i,0],c[i,1]))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mds = MDTrSampler(y, n_dim = 2, conv_size = 50, n_samples=100, batch_size=100, manifold_size=64)\n",
    "\n",
    "nmpsdm, psddu, prob_dist = mds.traj_char(c,s)     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGaxJREFUeJzt3X+QXWV9x/H3d3eTTQKEH826WXsTQvihidUibtJSja3WKGX4UR3qAJYSwMlYp1RaWpXSQsEyU3GM+KMdyzTIFAnUSh2SWAeSYkVmWtmNruQXgQiBXLObLGObYCJJlnz7x3Oud3/vvXvP7rnnuZ/XzJ1z7zkn5z7Pzdnvee73ec5zzd0REZF4NGVdABERSZcCu4hIZBTYRUQio8AuIhIZBXYRkcgosIuIREaBXUQkMgrsIiKRUWAXEYlMSxZvOm/ePF+0aFEWby0ikltbtmx5xd3bJtovk8C+aNEiuru7s3hrEZHcMrOXKtlPqRgRkcgosIuIREaBXUQkMgrsIiKRUWAXEYmMAruISGQU2EVEIpPJOPZadXXB1VfDwAC0JDW49VZ4y1vgwx8O6wYG4L3vhbVrYcUK6Osr7z8wAG1t8PTTUCyG7aXjDAzA978PhUL5/W67DR54YOg+peOMthxvn/nzw/ErsXQp/OIXI4/z1rfC+vXpfJYiEh/L4jdPOzs7fbI3KN1/P1x3XTrlaG6G118ffVuhALNmwYEDcOhQOu9XMmMGnHnm+BeGF16AEyfGPsasWeEiMdkLTFr7DN/25S/DJZek+3mJSGBmW9y9c8L98hTYu7pg+fLwvLkZFiwIAeWll+D48fL6s86CffvgyJHyv21rg5NOCvsfPAj9/eVtTU2weHF4PlZAbW0N7weTD4SHDw9930qULjCl40zFhSZtbW1w6qmVXRjWrYNly7Itr0heVBrYc5WKWbYspE+uvRZ27Bi6bflyOHYMenrK6264AZ54Am6/HVatGrr/xo1w440we/bIYw1PgaSZ+rj/frjjjsouDMNTQiW33QYPPTT1rfFqW+x79oTn/f2VX8CWLx958Rq+/MIX9C1A6suaNeHb6Vh/o1nLVYtd6t/y5SGoV3Jh2Lt37FTYcKWLushUu+yy8Rtya9bAzTeXX+/dO3pwv+wy2Lq13NeXhihTMRKfpUvLabTRgn+xGJYAixbBiy9mVlRpAEuXws6dcOmlYwf3lSth8+by67lzQ3p3sOXLQ+q4ZP58OPnk0L82PENQjUoDu4Y7SqZ27IDnnw+PF18cuTx+POTsIaR6brop0+JKxNasCUEd4EMfGn2frq6Qfmlqgq99Law7dCi0zks2biwH9ebmsOzrg927w/GXLp2a8g+mwC5178ABWLIkPH/kkWzLIvFasSIEYrPQD1YsjtxnwwY4ejT02a1aBddfH9b/7/+W9ymlDN/1rvBt89d/Pbw+7bRwHtfSYq+UArvkwo4dIY9ZLEJHR9alkRjt3x9GxH34w/Dyy0MHYpSURuV98INheccdYQTYU0+Vz8svfjEs3/GOsOzpgb/5mxD8pyOogwK75Mj73x+WfX0K7jI13EPwbm2F9vaR2597bujrQqE8XLevL7T2Dx0Kefd77invd+edU1fm0SiwS26sXVv+WtvXFzqxRNI0Zw6cccboo7W6uuCv/zoMzR0c9DdtKp+XJdu3T205J6LALrnS01P+I9q8eejIA5HJKhZDuuS888JNiseOjZz6Y/9+eO01uOuukTfV9fSEztRzzhl7+ON0UmCX3OnpKXemPvhgtmWROPT2wrPPwmc+E1IxZqHlPlh7e0jRrFgx+jFWrQqjubIO6qDALjn10Y+G5dy52ZZD4rBhQ2iNP/00nH9+aDh87nNDR8bs3x9GxOzfn105K6XALrm0YkXIdT7wwOjD0kQmY/Hi0OL+y78MI1iGj4zJ4H7OSVFgl1xatiz88e3ZM/QuQJFqdXXB3XePvu2VV8rP29tHdpzWKwV2ya3SmOJ587Ith+Rfc3PoPC1NFlg6pwafW9u3w8yZ0160ScnV7I4iIlOhqSnMD1MyvHXe1QUf/zj84z/mY5pptdgl18yyLoHkXUcHLFw49Ka3/fvDPEUbNoTXGzaEqbxfeCGbMlZLgV1yq7093FAyOA8qUq3e3tAB39tbXnfJJfBXfxXGrG/cWE77lZb1ToFdcqujI/zE4I036kYlmbyOjjAaZvg0FZdeGnLv7e3l+V/yQoFdcqtQgM9+dvzfhhWpxGidohs2hHTMb/1WGHlVCvJ5UHNgN7MFZvZdM9tpZtvN7BNpFEykEnn5Q5P8Kv3Qy3vek4+OU0hnVMwAcLO7/9DMTgG2mNkmd5+mCSpFRNJXmpHxoYfCFLwPP5xteapRc4vd3Xvd/YfJ81eBncCv1npcEZGs3XlnmP8lT0EdUs6xm9ki4O3AD0bZttrMus2su7/Sn7AXEZkGra1ZlyBdqQV2MzsZeAS4yd0PDd/u7ve6e6e7d7aVfsRSJCV5mJhJ6lOhAJ//fH3MypiWVAK7mc0gBPUH3f3f0zimSKWamuBTn9JkYDI5XV1w661xnT9pjIoxYC2w093X1F4kkcotWwZf/nL4jcrBN5iIVKJYhJtvDjciqcU+1DuBa4D3mllP8rg4heOKVOQtb9FYdqlNbL+hW/NwR3d/CtCMHSKSOzHm10F3nkokmnQmyySUUjEx5ddBgV0iMNrsfCKVevXVrEuQPgV2yb3RZucTqURvb5iKN7ZzR4Fdcm+s2flEJtLREfLrsZ07CuwShbz8ZJnUl95e2LdPLXYRkWgsWwYPPJCfWRsrpcAuIg2rWIS779aoGJG61Noa39dpkclSYJfcK91kEuN4ZJlaukFJpI7FNqpBpkexGN8EYKDALpGIteUlU6tQgPvui++8UWCXKMR6a7jIZCiwSzSOHs26BJI3xSJcf318DQIFdhFpWIVCfHOxgwK7RCS2362UqafOU5E6ps5TmQy12EXqmDpPZTLUYhcRiYyGO4rUMaViRMoU2CUKsX6llqml4Y4idSzWTjCZWkrFiNQxdZ6KlCmwi0jDUipGpI6p81QmQ6kYkTqmVIxMRrEYX1AHBXYRaVCxpmFAgV0ioVSMVCvWNAykFNjN7D4zO2Bm29I4nki1lIoRKUurxX4/cFFKxxIRmXJKxUzA3Z8EfpbGsUQmQ6kYqZZSMSkws9Vm1m1m3f39/dP1ttIgSlMKdHVlXRKR7E1bYHf3e929090729raputtpUGUphTQfDFSKaViRHJg2TLNFyOVUypGJAc0MkaqFWNQh/SGOz4E/DfwJjMrmtkNaRxXRGQqxdoISGtUzFXu3uHuM9y94O5r0ziuSDU0MkaqoRy7SA4oFSPVUI5dJCeOHs26BJInMQZ1UGCXyLS2Zl0CyZNYv90psEs0lGOXaijHLpIDyrFLNZRjF8kJ5dhFFNglMsqxS6WUihHJAeXYpRpKxYjkQGmGxxhbYDI1YgzqoMAuESnN8BjrH6ukL9ZGgAK7REOjYqQayrGLiERGOXaRHFDnqVQr1nNFgV2ioc5TqVas54oCu0RDnadSDeXYRXJALXaphnLsIjmgFrtIoMAu0VCLXaqhVIxIDqjFLtVQKkYkB9Ril2rFGNRBgV0iUmqBiVQq1kaAArtEJ9a8qaTo8GF+sXgJx85+E8Wdr2ZdmtS1ZF0AkTQpzy4TOn4czjmH2X19LAb4/U7YtSvrUqVKLXaJivLsMqG//Vvo6yu//sIXMivKVFFgl6gUTj/Mo72dFH64Hq69Vr+VJxM7eBBOnMi6FKlSYJd4uMOFFzJ72xa4/HL4l3+BVauyLpXUmz/8QwBeWXM/f7fk63D11fDkkxkXKl0K7BKP556DrVuHrtu6NeRURUp+9CMA5rXPYPVtHWGdWuwjmdlFZrbLzHab2afTOKZIVV59FS64AIDDp7SV12/fDgsXwuHDGRVM6op7+XlHR3blmGI1B3Yzawb+Afg9YClwlZktrfW4IhV77TWYPx+OHAFg68BSDjz8BFx8cdje1wednUP/qKUxnDgBmzfDjTfCsWOwezd85CPw9a9TPPu3+fu7BsJ+X/pS2B6JNIY7Lgd2u/sLAGb2MHA5sCOFYw914gRs2ABXXAFm4XVzc/iDPXECmppGX9bLPvVSjtjqU1oPsGQJzV/9Dm9492x48pvlc+fZZ6GlJZw3+lwbp6yvv14+N77ylfL5UCxSWNjEHec+ANuARx+F1laYOXPqyjpjBrS3wzPPwNy5qYfHIdy9pgdwBfDPg15fA3xllP1WA91A98KFC31Svvc99/Bx6aHHyMeb3+x79wz4ypXue/e6+9Gj7h/7mPucOdmXTY/6erzjHSGmPP749L/34sWTi38hkHZXEpctCbqTZmZ/AHzA3T+avL4GWO7uN471bzo7O727u7v6N1OLvb72qZdyzJkDy5bBxo0wezZdXeHlLx08CAsWhJRN1mXN0+caQ1khtI7nzoV9+0LcePOb4amn4KSTKL58gsKOx+Ev/iLcpNTUNHVlTaHFbmZb3L1zov3SSMUUgQWDXheAfSkcd6SmpjCMTaMcZAylG5SGzNp36qlw6FCm5ZL6UyzC9R9t4r77LqKw7aKsi5OqphSO0QWca2ZnmdlM4EpgfQrHFalazFOxSrpiPldqbrG7+4CZ/QnwGNAM3Ofu22sumYiITEoaLXbc/T/c/Tx3P9vd70rjmCKTEfOv4ki6Yj5XUgnsIvVCsztKpWJOxSiwS1Q0u6NUI8agDgrsEpmYW2GSvlgbAArsItKQlGMXyYmY/1glXTF/u1Ngl6io81SqEet5osAuUVHnqVQj1vNEgV2iEvPXa0lXzGk7BXYRaUgxNwIU2CUqMbfCJH0xBnVQYJfIqPNUKhXzxV+BXaKizlOpROzf7BTYJSox500lPbGfJwrsIiKRUWCXqMT+FVvSEft5osAuUVHnqVRCqRiRHFHnqVQq1qAOCuwSGbXYpVIxX/wV2CUqarFLJZRjF8kRtdilEsqxi+SIWuxSqViDOiiwS2TUYpdKxXzxV2CXqKjFLpVQjl0kR0q5U5HxKMcukkMxt8ZEJqLALtGJvTUmtVMqRkQkMrFf/GsK7Gb2B2a23cxOmFlnWoUSqUXsrTGRidTaYt8GfAh4MoWyiKRCQx5lIrFf/GsK7O6+0913pVUYkTRoyKNMRKmYlJjZajPrNrPu/v7+6XpbaUBqsUslYj4/JgzsZrbZzLaN8ri8mjdy93vdvdPdO9va2iZfYpEJqMUulYj5/GiZaAd3f990FEQkLWqxy0RKOfZY0zEa7ijRKRbh5pvjbpFJbZRjH4eZfdDMisCFwLfN7LF0iiUiIpNV66iYb7l7wd1b3b3d3T+QVsFEJqtQgM9/Pt7WmNROwx1FckadpzIRpWJEckadp9LoFNglOuo8lYkoFSMiEhmlYkRyRp2n0ugU2CU6SsXIRJSKERGJjFIxIjmjVIw0OgV2iY5SMTIRpWJERCKjVIxIzigVI41OgV2ioykFZCJKxYjkjKYUkIkoFSOSM+o8lUanwC4iDUepGJGcKRRg3bqsSyH1TKkYkZyKuUUmMh4FdomSOlBlPErFiOSQhjzKeJSKEckhtdilkSmwS5TUYpfxKBUjkkNqsct4Yj8/FNglSmqxy3hiPz8U2CVKsbfIpDaxnx8K7BKl2FtkUpvYzw8FdolS7C0yqY2GO4rkkCYCk/EUi/EGdagxsJvZ58zsWTN7xsy+ZWanpVUwEZGpEPtQR6i9xb4J+DV3fxvwHHBL7UUSqZ1+RUnG0ghpupoCu7s/7u4Dycv/ASL+qCRPYu8ck8lrhHMjzRz79cB3xtpoZqvNrNvMuvv7+1N8W5GRGqFVJpMTe8cpVBDYzWyzmW0b5XH5oH1uBQaAB8c6jrvf6+6d7t7Z1taWTulFxqDOU2lkLRPt4O7vG2+7mV0LXAL8rrt7WgUTEZkKpc7TmFvtEwb28ZjZRcCngN929yPpFEmkduo8lbE0Qpqu1hz7V4BTgE1m1mNmX02hTCI1UypGxtIInac1tdjd/Zy0CiIiMh3UeSqSU0rFSCNTYJcoKRUjY9GdpyIikVEqRiSnlIqRRqbALlFSKkbGolSMSI4dPZp1CaQeKRUjkmOtrdDbm3UpRKafArtEqZRjVzpGhlMqRiTHOjqUjpGRlIoREZHcUWAXkYaiVIxIzh07lnUJpN5odkeRHOvthRde0MgYGaoRZndUYJdodXSEVllHR9YlkXqizlORHOvthT171GKXxqPALtHavx+OHAlLkZJiEa6+WqkYkdxqbobnnsu6FCLTS4FdonXJJXD33XD77dDVlXVppF40wsyfNf00nki9O+88OHxY6ZjRrFwZ+iAABgagpWXs5bvfDWvXjn6cc8+tbv+slWb+XLcu3uCuwC5Ra2+HWbPCMmtnnTU0+MHYgXS8bWns09tb3XQLu3fDww/DG9849DgvvwwnToy+/ze+AfPmhW9Mq1aN3GfjRvizPwvHeetbYf36yssj41Ngl+i5j73tsstg586RgbBQgO9/P7y+6Sb49rdrC6gvvTR+ObJy2mkh+I5Xn1degZ//PHRE79498hhNTbBwYXn/AwfCvj//eXhcdx188pPhotDTA1deCd/7HvT1lY+xZw+cdFL4djWe88+HgwfL7/Xii9XXWakYkZzbvx9eew22b4dly4Zumz07bBvNnj0wc2bY59ChdMpiBmefXR8t9pYWuOoquPPOysp+ww3wxBOjH+f550fuf+WVsGVLuKAdPw79/eFhNnS/lhaYMyd8xkeOhM7uxYtHL3OxGJaDNTXBmWeOrPO6dUP/vwd/O7jmGnj8cfjmN+MN7uYZNCM6Ozu9u7t72t9XGs+aNSGf2tISWneFAixfDt3d5Rb03LlwxhnloPCzn4XHYLNmwfz5tQXU0QJgI1ixIrT0B7fQTz45XOR6esLrm26CL36xsuPNnx/+/U9+Mv63oEIB3v52WL0aLr106LZZs+DJJ0de7OudmW1x984J91Ngl9h1dISgMjjgQmg9/umfwj33jPw3998Pd9wR/s2SJcr/pqHUah6rY7VYDBeB8S6Yw1vio3XcvvgivP76yOPPmhX6A44dC8d4+umpq+tUUWAXSRSLsGDB0HVtbSEXLHFauhR++tOhabSvfS104q5cCZs2ZVa0mlQa2DWOXaJXKMDevXDOObBoEWzYoKAeux07QifrJz4R/t9LQR3yG9SroRa7iEhOTEuL3cw+Y2bPmFmPmT1uZm+s5XgiIlK7WlMxn3P3t7n7+cBG4LYUyiQiIjWoKbC7++ARvicBdXgLhohIY6n5BiUzuwv4I+Ag8J5x9lsNrAZYuHBhrW8rIiJjmLDz1Mw2A/NH2XSruz86aL9bgFnufvtEb6rOUxGR6lXaeTphi93d31fhe64Dvg1MGNhFRGTq1Doq5txBLy8Dnq2tOCIiUquaxrGb2SPAm4ATwEvAx9z9pxX8u/5k/8mYB7wyyX+bV6pzY1Cd41drfc9097aJdsrkBqVamFl3JTmmmKjOjUF1jt901VdTCoiIREaBXUQkMnkM7PdmXYAMqM6NQXWO37TUN3c5dhERGV8eW+wiIjKOXAV2M7vIzHaZ2W4z+3TW5amFmd1nZgfMbNugdWeY2SYzez5Znp6sNzP7UlLvZ8zsgkH/5tpk/+fN7Nos6lIJM1tgZt81s51mtt3MPpGsj7nOs8zsaTP7cVLnO5L1Z5nZD5Ly/6uZzUzWtyavdyfbFw061i3J+l1m9oFsalQ5M2s2sx+Z2cbkddR1NrM9ZrY1mem2O1mX3bnt7rl4AM3AT4DFwEzgx8DSrMtVQ33eDVwAbBu07m7g08nzTwOfTZ5fDHwHMOA3gR8k688AXkiWpyfPT8+6bmPUtwO4IHl+CvAcsDTyOhtwcvJ8BvCDpC7fAK5M1n8V+OPk+ceBrybPrwT+NXm+NDnfW4Gzkr+D5qzrN0Hd/5xwN/rG5HXUdQb2APOGrcvs3M78A6nig7sQeGzQ61uAW7IuV411WjQssO8COpLnHcCu5Pk/AVcN3w+4CvinQeuH7FfPD+BRYGWj1BmYA/wQ+A3CDSotyfpfntfAY8CFyfOWZD8bfq4P3q8eH0AB+E/gvYTpvK0B6jxaYM/s3M5TKuZXgb2DXheTdTFpd/degGT5hmT9WHXP5WeSfN1+O6EFG3Wdk5RED3AA2ERoef6fuyc/qT2k/L+sW7L9IPAr5KzOwD3AJwl3pEOoQ+x1duBxM9uSzGQLGZ7bNU/bO41slHWNMqRnrLrn7jMxs5OBR4Cb3P2Q2WhVCLuOsi53dXb314Hzzew04FvAktF2S5a5r7OZXQIccPctZvY7pdWj7BpNnRPvdPd9ZvYGYJOZjTdv1pTXOU8t9iIw+LfmC8C+jMoyVfabWQdAsiz95PJYdc/VZ2JmMwhB/UF3//dkddR1LnH3/wP+i5BTPc3MSo2qweX/Zd2S7acCPyNfdX4ncJmZ7QEeJqRj7iHuOuPu+5LlAcIFfDkZntt5CuxdwLlJ7/pMQkfL+ozLlLb1QKkn/FpCHrq0/o+S3vTfBA4mX+0eA95vZqcnPe7vT9bVHQtN87XATndfM2hTzHVuS1rqmNls4H3ATuC7wBXJbsPrXPosrgCe8JBsXQ9cmYwgOQs4F3h6empRHXe/xd0L7r6I8Df6hLt/hIjrbGYnmdkppeeEc3IbWZ7bWXc6VNlBcTFhNMVPCD/0kXmZaqjLQ0AvcJxwpb6BkFv8T+D5ZHlGsq8B/5DUeyvQOeg41wO7k8d1WddrnPq+i/C18hmgJ3lcHHmd3wb8KKnzNuC2ZP1iQpDaDfwb0Jqsn5W83p1sXzzoWLcmn8Uu4PeyrluF9f8dyqNioq1zUrcfJ4/tpdiU5bmtO09FRCKTp1SMiIhUQIFdRCQyCuwiIpFRYBcRiYwCu4hIZBTYRUQio8AuIhIZBXYRkcj8P8Gh5liQhyRhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_b, ax_b = plt.subplots()\n",
    "ax_b.scatter(range(0,nmpsdm.shape[0]),nmpsdm[:,0], marker=\".\", c=\"b\", s=0.5)\n",
    "ax_b.scatter(range(0,psddu.shape[0]),psddu[:],marker=\"x\", c=\"r\", s=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampling_rate = 0.005\n",
    "target, time_stamps = mds.batch_sampling(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt0lNW9//H3JoSES1QgaVApBE6lNIQ4QKAgoLEqUO1Sy9GFVJGLPWDVCvhbIh6P4daug9QW23P0KD0C1lpFkdYWuagcUy1iMdiAASEqF+XgCeFOuMgl398fMxlzzyQzyWSefF5rzcrMM/vZ892TJ9/s2c+e/TgzQ0REvKNVtAMQEZHIUmIXEfEYJXYREY9RYhcR8RgldhERj1FiFxHxGCV2ERGPUWIXEfEYJXYREY9pHY0XTU5OtrS0tGi8tIhIzNq0adMBM0upq1xUEntaWhp5eXnReGkRkZjlnNsTSjkNxYiIeIwSu4iIxyixi4h4TFTG2EVi3dmzZ9m7dy+nT5+OdijiQYmJiXTt2pX4+PgG7a/ELtIAe/fuJSkpibS0NJxz0Q5HPMTMOHjwIHv37qVHjx4NqkNDMSINcPr0aTp37qykLhHnnKNz585hfRpUYhdpICV1aSzhHltK7CIiHhOxxO6ci3PO/cM5tzJSdYpIzeLi4vD5fMHb/PnzQ943NzeXH/zgBw1+7dr2T0tL48CBAwBcccUVDX6N8vLz81m1alW999u3bx+33HJLg15z6dKl7Nu3r9775eTk8NZbbzXoNSMlkidPpwIfAxdEsE4RqUHbtm3Jz8+Pdhi1eu+99yJST35+Pnl5eVx//fVVnjt37hytW1efyi655BKWL1/eoNdcunQpGRkZXHLJJSHvc/78eebOndug14ukiPTYnXNdgRuA/45EfSLScGlpafzrv/4rQ4YMISsriw8//JCRI0fyT//0Tzz99NPBcseOHeOHP/wh6enp3H333ZSWlgLwxhtvMGTIEPr378+tt95KSUkJAGvWrKF3794MGzaMFStWBOs5ePAgI0aMoF+/fkyZMgUzCz7XoUMHwN/Dz87O5pZbbqF3797cfvvtwXKrVq0K1nv//fdX+SRw5swZcnJyWLZsGT6fj2XLljF79mwmT57MiBEjuPPOO9m9ezfDhw+nf//+9O/fP/gPZffu3WRkZAD+pPvggw8ycOBAMjMzeeaZZ4KvsWDBAvr27cvll1/OzJkzWb58OXl5edx+++34fD5OnTrFunXr6NevH3379mXSpEl89dVXwfd77ty5DBs2jFdeeYUJEyYE/5ls2rSJq666igEDBjBy5Ei+/PJLAH7zm9+Qnp5OZmYmt912W7i/8qrMLOwbsBwYAGQDK+sqP2DAABOJZdu2bYt2CNaqVSu7/PLLg7eXXnrJzMy6d+9uTz31lJmZTZs2zfr27WvHjh2z/fv3W0pKipmZvf3225aQkGCfffaZnTt3zq699lp75ZVXrLi42IYPH24lJSVmZjZ//nybM2eOnTp1yrp27WqFhYVWWlpqt956q91www1mZvbTn/7U5syZY2ZmK1euNMCKi4vNzKx9+/bB17vgggvsiy++sPPnz9vgwYPt3XffDda7c+dOMzO77bbbgvWWt2TJErv33nuDj2fNmmX9+/e3kydPmpnZiRMn7NSpU2ZmVlhYaGU5ZteuXdanTx8zM3vmmWds3rx5ZmZ2+vRpGzBggO3cudNWrVplQ4YMsRMnTpiZ2cGDB83M7KqrrrIPPvjAzCwY544dO8zMbNy4cbZw4cLg+/3YY48FYxs/fry98sordubMGRsyZIjt37/fzMxeeuklmzhxopmZXXzxxXb69GkzMzt8+HC1v9/qjjEgz0LIyWEPxTjnfgDsN7NNzrnsWspNBiYDdOvWLdyXFYlJRUVFpKamRqSu2oZibrzxRgD69u1LSUkJSUlJJCUlkZiYyJEjRwAYNGgQPXv2BGDs2LH87W9/IzExkW3btjF06FDA31seMmQI27dvp0ePHlx22WUA3HHHHSxatAiAd955J9iDv+GGG+jYsWO1MQ0aNIiuXbsC4PP52L17Nx06dKBnz57B+dpjx44N1luXG2+8kbZt2wL+L4zdd9995OfnExcXR2FhYZXyb7zxBlu2bAn2po8ePconn3zCW2+9xcSJE2nXrh0AnTp1qrLvjh076NGjB7169QJg/PjxPPnkk0ybNg2AMWPGVLtPQUEB1113HeD/xHDxxRcDkJmZye23387NN9/MzTffHFJ76yMSY+xDgRudc9cDicAFzrnfm9kd5QuZ2SJgEUBWVpZVrUbE24qKirjnnnt46qmnIpbca5KQkABAq1atgvfLHp87dw6oOqXOOYeZcd111/Hiiy9WeC4/P7/WKXihTM8rH0dcXBznzp2rMGxTX+3btw/eX7hwIampqWzevJnS0lISExOrlDcz/uM//oORI0dW2L5mzZo6468rzvKxlN+nT58+bNiwocpzr7/+Ou+88w5//vOfmTdvHlu3bq3xPEFDhD3GbmYPm1lXM0sDbgP+p3JSFxFITU1tkqQeqo0bN7Jr1y5KS0tZtmwZw4YNY/Dgwaxfv55PP/0UgJMnT1JYWEjv3r3ZtWsXn332GUCFxH/llVfywgsvALB69WoOHz4ccgy9e/dm586d7N69G4Bly5ZVWy4pKYnjx4/XWM/Ro0e5+OKLadWqFc8//zznz5+vUmbkyJH813/9F2fPngWgsLCQEydOMGLECBYvXszJkycBOHToUJXX7N27N7t37w6+L88//zxXXXVVrW379re/TXFxcTCxnz17lq1bt1JaWsoXX3zB1VdfzYIFCzhy5EjwPEakaB67SBOKZFI/depUhemOM2fOrNf+Q4YMYebMmWRkZNCjRw9++MMfkpKSwtKlSxk7diyZmZkMHjyY7du3k5iYyKJFi7jhhhsYNmwY3bt3D9Yza9Ys3nnnHfr3788bb7xRr6HWtm3b8tRTTzFq1CiGDRtGamoqF154YZVyV199Ndu2bQuePK3snnvu4bnnnmPw4MEUFhZW6EGX9cZ//OMfk56eTv/+/cnIyGDKlCmcO3eOUaNGceONN5KVlYXP5+Pxxx8HYMKECdx99934fD7MjCVLlnDrrbfSt29fWrVqxd13311r29q0acPy5ct56KGHuPzyy/H5fLz33nucP3+eO+64g759+9KvXz+mT5/ORRddFPJ7FgoXzkehhsrKyjJdaENi2ccff8x3vvOdaIfhCSUlJXTo0AEz49577+Wyyy5j+vTpEal706ZNPPDAA/z1r3+NSH1NqbpjzDm3ycyy6tpXPXYRiarf/va3+Hw++vTpw9GjR5kyZUpE6s3Ly2Ps2LFMnTo1IvXFEvXYRRpAPXZpbOqxi4hIkBK7iIjHKLGLiHiMEruIiMcosYvEqJ///Of06dOHzMxMfD4ff//732stP3v27OAc7Wgqv6xvY3j66af53e9+V2uZhi4DHCt0zVORGLRhwwZWrlzJhx9+SEJCAgcOHODMmTPRDqtZqOuLQ1D7MsBeoB67SFNYswbKphab+R+H4csvvyQ5OTm4/kpycnJw3fDyPeK8vDyys7OD+23evJnvfe97XHbZZfz2t78N1nXllVfi8/nIyMjg3XffBeAnP/kJWVlZ9OnTh1mzZgXrCGVZ4NzcXK688spqlwUu7/e//z2DBg3C5/MxZcqUapcCSEtL46GHHmLQoEEMGjQo+LX+PXv2cM0115CZmck111zD559/DlT8ZJKdnR3ct1evXrz77rvVLgPsNUrsIo1tzRr4/vdh+nR/Up8+3f84jOQ+YsQIvvjiC3r16sU999wT8jcrt2zZwuuvv86GDRuYO3cu+/bt4w9/+AMjR44kPz+fzZs34/P5AP9QT15eHlu2bOGvf/0rW7ZsCdbzzW9+kw0bNjB8+PDg+uPvv/8+OTk5wTIbN27kl7/8JR999BGfffZZhTXcwT9Pe9myZaxfvz64KmPZmjOVXXDBBWzcuJH77rsvuKLifffdx5133smWLVu4/fbbuf/++6vd99y5c2zcuJEnnniCOXPm0KZNG+bOncuYMWPIz8+vdmXGWKfELtLYRo6EqVPh17+GVq38P6dO9W9voA4dOrBp0yYWLVpESkoKY8aMYenSpXXud9NNN9G2bVuSk5O5+uqr2bhxIwMHDmTJkiXMnj2bjz76iKSkJABefvll+vfvT79+/di6dSvbtm0L1lN+WeDvfve7JCUlkZKSUu2ywHFxccFlgctbt24dmzZtYuDAgfh8PtatW8fOnTurjXvs2LHBn2WLam3YsIEf/ehHAIwbN65K/WVGjx4NwIABA4KLjXmdxthFGptzsHChP6GXWbjQvz0McXFxZGdnk52dTd++fXnuueeYMGECrVu3Dg57nD59ulIoVZfqvfLKK3nnnXd4/fXXGTduHA8++CDDhw/n8ccf54MPPqBjx45MmDChQl0NXRa4PDNj/Pjx/Pu//3udbS2/b01L7Na0vSy+sqWCWwL12EUaW9nwS3llwzINtGPHDj755JPg4/z8/OCKi2lpaWzatAmAV199tcJ+r732GqdPn+bgwYPk5uYycOBA9uzZwze+8Q3+5V/+hbvuuosPP/yQY8eO0b59ey688EKKiopYvXp1vWOsblng8q655hqWL1/O/v37Af9yuXv27Km2rrJx8GXLljFkyBDAf6Hsl156CYAXXnihSv21qWsZ4FinxC7S2Nau/Xr4pbT062GZtWsbXGVJSQnjx48PXjdz27ZtzJ49G/Avozt16lSGDx9OXFxchf0GDRrEDTfcwODBg3n00Ue55JJLyM3Nxefz0a9fP1599VWmTp3K5ZdfTr9+/ejTpw+TJk0KXlGpPqpbFri89PR0fvaznzFixAgyMzO57rrrgtcEreyrr77iu9/9Lr/+9a9ZuHAh4L9u6JIlS8jMzOT555/n1+U/EdWhrmWAY50WARNpgHovArZmjX9M3Tl/T33tWhg1qvECjLLc3Fwef/xxVq5cGXZdaWlp5OXlkZycHIHIYkc4i4BpjF2kKZRP4s55OqlL9Cmxi0jElZ3UjYSWMpMlkjTGLtJA0RjGlJYh3GMr7MTunEt0zm10zm12zm11zs0Jt06R5i4xMZGDBw8quUvEmRkHDx4kMTGxwXVEYijmK+B7ZlbinIsH/uacW21m70egbpFmqWvXruzdu5fi4uJohyIelJiYSNeuXRu8f9iJ3fxdlpLAw/jATd0Y8bT4+Hh69OgR7TBEqhWRMXbnXJxzLh/YD7xpZrWvHyoiIo0mIondzM6bmQ/oCgxyzmVULuOcm+ycy3PO5enjq4hI44norBgzOwLkAlUm6ZrZIjPLMrOslJSUSL6siIiUE4lZMSnOuYsC99sC1wLbw61XREQaJhKzYi4GnnPOxeH/R/GymYX/PWIREWmQSMyK2QL0i0AsIiISAfrmqYiIxyixi4h4jBK7iIjHxPTqjkOHDuXEiRPEx8dz9uzZ4Pbyj+t732v7xMfH06FDB95++22uvvpqSkpKqi23du1aUlNTQ3vjRaRZi9nEPnToUN57771ohxEzaroeZJkuXbqQnp5OfHx8k/wDys/Pj1zjRKSCmE3s69evV489hPunTp2isLAwuH+vXr1o27ZthXI7duzg9OnTFa5C39gSEhLIzMwMqW3jx49neuVrhopEwcKFC1m+fDnr16+Pdii10qXxWoDc3FymTZvGE088UePFD8r+SULj/wPasmULZ86cqVcbunTpwpNPPsno0aPrtZ9IKIqKihg/fjxr1qypsczChQt54IEHAGjbti0nT56stpzP5wNg9OjR5OTkUFBQwC9/+UuWLFkSdpyhXhpPiV2iwufzhfTPYPfu3Rw9ejRY5le/+pV67xJRRUVFfOtb36KkpISRI0fWmNz79OlT4VNtUlISx44dq1CmU6dOHD58OPi4c+fOHDx4EIAJEyaEndxDTeyYWZPfBgwYYCKhmjBhgrVr187wLwdtH330UbRDEg959dVXg8fWnDlzqi0zZ84cAyw9Pd0++ugjS0pKMsCysrKCZW666SYDrGPHjnbhhRcG6wSsc+fOEYkVyLMQcqymO0qzt2TJEk6cOEHZ4nFXXHEFixcvjnJU4jWXXnops2bNYsWKFVWey8jwL1g7b948MjIyeOKJJwDIy8vD5/MxY8YMXnvtNQCmTZvGkSNHmDBhAllZWVxxxRUcOHCg6RqC5rFLDHn55ZcBOH78OHfddRdz586NckTiJZMnT6Zt27b06tWrynN/+ctfADhy5AgAkyZNIj09HYDNmzfzi1/8AoBu3bqRk5MD+DskH3zwQVROtCqxS8zIzs7m7bffpn379gDMmjWLgoKCKEclsa5Tp04kJSUB/llkGzdurPD84sWLWbp0KRMmTGDSpEnB7Vu3biUr6+vh7o4dO7Jnz56mCboOSuwSU7KzsykpKaFLly4ATJkyJcoRSSwrKipi0qRJtGnTJpjcy5+sB+jZsyfgn3Zb2QcffMDbb7/NTTfdxKFDhxo/4BApsUtMevHFFwF47733qh0TFQnFM888w65duxgzZgw/+tGP6Nu3b7WfBFu1qjlVZmdn86c//amxQ60XJXaJSdnZ2cyZMyfaYUiMK5u+2L59e1JTU5k4cSLHjx+vMBxz6NAhSktLm1WPvC5K7BKzymYqiDREbm4uf/zjHxkzZgwLFiwAoHv37gBcdNFFwXJ79uzhoosuqvakanOlxC4xq2yGwvvvvx/lSCQWJScn07ZtW/7t3/4tuK1Tp04Vhl1WrFjBAw88QE5OTkx1JJTYJWZNmjSJBx98kIULF5KbmxvtcCQGVTd23qpVKx555BGKioqCY+3Hjx9v6tDCosQuMW3w4MGcO3cupsY/pfmonNizs7N59NFH2b59O3/4wx9Yu3YtEHvDfmGv7uic+ybwO6ALUAosMrNfh1uvSCh69epF+/btY2r8U5qHlJQUvvWtbwW/0VympKQEILjgV+vWrWPu+IpEj/0c8P/M7DvAYOBe51x6BOoVqdOBAwf46quvmvwr2xL7iouL2blzJ8XFxRW2Dx48uMLjRx99NOZ67GEndjP70sw+DNw/DnwMXBpuvSKhKDsBlpycHO1QJMakpKTQu3fvKj320aNH8+qrr3LFFVfw7LPPBpcIiCURvdCGcy4N6Af8vZrnJgOTwb+egkgklPXY33zzzZjrVUn0JSQkVLt99OjRMb32f8ROnjrnOgCvAtPM7Fjl581skZllmVlW5f+QIg2VnZ3N/PnzmTFjhmbGiAREJLE75+LxJ/UXzEzf75Ym1b17d82MESkn7MTu/FdJfhb42Mx+FX5IIvXTqVMn4uLi6NSpU7RDkRhSXFxMmzZtoh1Go4hEj30oMA74nnMuP3C7PgL1ioQkOTmZDh066ASqhKyoqIiZM2fy2GOPkZqaGu1wIi7sk6dm9jfARSAWkQarbfU9keocO3asyowYr9Bfg3iCErvUR3FxMdu2basyh90r9NcgMS8lJYV+/fpFOwyRZkOJXWJeamoqCxcuZNq0aRQVFUU7HIkBKSkppKenayhGpLnbsmWLZz9aS+TV9OUkL1BiF0+o6evhIjXx6lRHUGIXkRYoNTXVs1MdQYldPKK4uJjt27drKEZCUjaP3avnZJTYxRM0FCP1ZWbRDqHRKLGLZ3j5ZJhEloZiREQ8pqioiFmzZmkoRqS5U49dQpWamsqcOXPUYxdpzlJTU5k/f75n/1AlstRjF4kBRUVFPPTQQ579Q5XIUo9dJEacOXMm2iFIjFCPXUTEY9RjFxHxGPXYRWKEZsVIqNRjF4kBqampPPvss9EOQ2KEeuwhcM4tds7td84VRKI+kYaaNGmSZ/9YJXJSU1N56qmn1GOvw1JgVITqEmkwzYwRiVBiN7N3gEORqEtEpLEVFRVxzz33ePbTXZONsTvnJjvn8pxzeVpaVUSiTUMxEWBmi8wsy8yytLSqiERLWW/dyzQrRkRaFK+fOAUldhERz4nUdMcXgQ3At51ze51zd0WiXhGRSPP6iVOA1pGoxMzGRqIeEZHGpqEYERGJOUrsItKitIShGCV28RQtBCZ10VCMSAzR5fFE/JTYxTN0eTwJhYZiRGKMFgGTunh9LXZQYheRFsbra7GDErt4jE6eSl108lQkhujkqYifErt4RtlH7IICXchLaqaTpyIxpOykmNfHTyU8GooRiTEZGRmen/EgUhcldvEUzWWXumgoRiQGOeeiHYI0YxqKEYkxmhkjosQuHqOhGKmLhmJEYpCGYqQ2WlJAJMZoKEbqoiUFQuScG+Wc2+Gc+9Q5NzMSdYo0REv4o5Xw6ORpCJxzccCTwPeBdGCscy493HprtWYNmPnvm/kfS/PQWL+byvX+7Gdf1z9uHPh80LMnqd/+Ni//5S+kXnYZJCZC+/ZwwQWQmekvt3o1zJvn/1lbjM3hGGsOMVRWn5jqKhtqXfUtt2YNrFrl/z2Xlvp/dukCrVv7b9/4Bp179ADn/LcePfxl5s71HxerV/v3X7Pm6+PFzH9bvfrrY698PDXFGK3foZmFdQOGAGvLPX4YeLi2fQYMGGANtnq1/y2eOtWstNT/E/zbJboa63dTud6bbvI/Hjy47M+tYbdVq6qPsTkcY80hhnBiqqtsqHXVt1zZsVF269IlvGOk7HbTTWb331/xcfl4qotx3ryI/w6BPAslL4dSqNYK4Bbgv8s9Hgf8Z237hJXYK7+Z5d84ia7G+t1UV6/PF5k/2OpibA7HWHOIIZyY6iobal3hlGusW+Vj7/77Kyb9shjPn4/47zDUxO78ZRvOOXcrMNLMfhx4PA4YZGY/rVRuMjAZoFu3bgP27NnT8Bc1g1blRpFKS/0fqST6Gut3U7ne8+chLi78eqH6GJvDMdYcYqisPjHVVTbUuhparrFUPvZKS/0/q4sxwr9D59wmM8uqq1wk3oW9wDfLPe4K7KtcyMwWmVmWmWWlpKQ0/NXMYPr0itumT/dvl+hqrN9NdfUOGBBeneVVjrE5HGPNIYbK6hNTXWVDrSucco2l8rE3bZr/Vt706f4kHq3fYSjd+tpuQGtgJ9ADaANsBvrUto/G2D1KY+yR0xxiCCcmjbHH9hi7/7W4HigEPgMeqat8WIndzP/GlB+rU1JvPhrrd1O53nnzvq7/jjvMLr/crEcPO3fBBXYK7FyHDmYJCWbt2pklJZn17esvt2qV2dy5/p+1xdgcjrHmEENl9YmprrKh1lXfcqtXm73+uv/3fP68/2dqqllcnP+WkmJn27b9OiGnpfnLzJnjPy5WrfLvv3r118dLaan/tmrV18de+XhqijHCv8NQE3vYY+wNkZWVZXl5eU3+utIyFBQUkJGREe0wpJkqW1IgFueyN+UYu0izoS8oSV20pIBIjGkJf7QSnpbwz1+JXTylJfzRSni0pIBIjFGPXUSJXTxGPXapS0tYj12zYsRzNCtG6lJUVBSTn+o0K0ZaJF1BSUSJXTxIV1CS2mgoppFoKEYak4ZipC4aihGJITp5KqLELh6j6Y5Sl5YwFKPELp6iHrvURV9QEokxZX+0IjWJ1fH1+lBiF0+aNGmSeu1SRUsYhgEldvEoTXmU6rSUczBK7OI5qampzJ8/3/N/vFJ/LeUcjBK7eI6+fSo1aQknTkGJXTxKQzHSkoWV2J1ztzrntjrnSp1zdX4bSqQpaChGaqKTp6EpAEYD70QgFpGIaCnjqFJ/OnkaAjP72Mx2RCoYkUhoKX+8Un8t5Z++xtjFc3TyVGrSUk6etq6rgHPuLaBLNU89YmavhfpCzrnJwGSAbt26hRygSEPo5Km0ZHX22M3sWjPLqOYWclIP1LPIzLLMLCslJaXhEYvUQSdPpSY6eSoSozQUIzVpKedfwp3u+EPn3F5gCPC6c25tZMISCc+ZM2eiHYI0Qy3l5GmdY+y1MbM/An+MUCwiIo2qpZw81VCMeFJCQkK0QxCJGiV28RydPJWa6OSpSIzSyVOpiU6eisQwnTyV6rSUk6dK7CLSYqjHLiLiMeqxi8SwhIQEiouLox2GNDOa7igSo8pmxbSEnplIdZTYxZMyMjJaxFiq1I+mO4rEME15lOro5KlIjNOUR6lMJ09FYpyWFZDK1GMXiWFaVkCqox67SAzTGLtURz12kRiny+NJZS2lx+7MrMlfNCsry/Ly8pr8daVlKSgoICMjI9phSDMTy8eFc26TmWXVVU49dvEkDcVIdYqKipg5c6bnjwsldvEsTXeU6kRjlKKphXvN018457Y757Y45/7onLsoUoGJhEvTHaWy1NRUHnvsMZ08rcObQIaZZQKFwMPhhyQSPk13lOq0lJOnYSV2M3vDzM4FHr4PdA0/JJHwaYxdqqPpjvU3CVgdwfpEwqIxdqmspfTYW9dVwDn3FtClmqceMbPXAmUeAc4BL9RSz2RgMkC3bt0aFKxIfWiMXarTEtZjD3seu3NuPHA3cI2ZnQxlH81jl6YQy/OVJfLKluyN5cTeJPPYnXOjgIeAG0NN6iJNQWPsUllLGV+H8MfY/xNIAt50zuU7556OQEwiEaExdimvpYyvQwhj7LUxs29FKhARkcakHruIiMe0pB67Ert4VkJCAsXFxdEOQ5oJ9dhFYlzZN091AlXKqMcu4gEpKSk6gSpB6rGLeMRXX30V7RCkmVCPXcQDiouL2b59u8bZBfD32GP5y0n1ocQunpWSkkLPnj1JSUmJdigiTUqJXTyruLiYnTt3qscugH8o5q677tJQjEgsU49dKmsJV08CJXbxsOLiYgoLC9VjF6DlXD0JlNjFww4cOMCxY8coLCyMdijSDGhWjIhHOOeYPXt2i/hjltppHruIB2RnZ7NgwQK2bdvGxx9/HO1wPKuoqAifz8eMGTOiHUqtioqKmDlzZov4Jx/W6o4izV337t05f/48hw4dinYoNVqxYgVz584lPj6es2fPAlS4X/lxKPebap+zZ8+ybds2ADZv3szSpUu55JJLaN++PevXr6+13aNGjSIzM5MFCxbU7w0LQ0v5wpoSu3hap06daNWqFZ06dar2+RkzZvDGG28AFRPW0KFDefLJJwGYO3cuK1asCO4TyYR5/Phxdu7cGbEGvrB0AAAHhUlEQVT2RltxcXHwZHW7du3o1asXAKNHjyYnJ4dRo0Zx8OBBtm3bxsmTJ1m7di1Lly5lwoQJtSb4hQsX8txzzwUf33///UyaNKlxGxPDwr40XkPo0njSVHJzc7n22msZN24cS5YsCW6fMWMGixcv5uDBgzXu27FjR1q3bt0ks2ratWtHenp6zPXYgWDv/N5772Xjxo01nrCOi4vj/PnzNb4HnTt3pkePHlVep6ioiP/7v/+rUv7SSy8lOTm5SmyXXnopf/rTnwD/5RHvuOMOANLS0jh06BCvvPJKzI6zh3ppPCV28bSCggIGDhzI6dOnSUlJoXv37uzatatCQu/cuTNdu3YNJoh9+/ZVSeYdO3YMXoQ90gmzc+fOrFmzppHegegoKChg4sSJnD17ls8//5zDhw8Hn+vZsydJSUnk5OSwbt06li1bVus/2DIXXnghaWlp7Ny5k+PHj9datn379nTv3j04TFQmPj6eDz/8MGavhavELhJw880389prr1XZ3rFjR3784x9XOwRQfvil/LCMNMy9997L+vXrycnJYfTo0TU+X9M/wH79+lX4xHXbbbexffv2KuUOHz7M559/XqHuxMREOnfuzP/+7//SqlUr1q1bR3Z2dmM0s9EpsYuUM3HiRP7xj38Ek4CStXctXryY3/zmN8THx9OmTZvgSdzFixfTs2fPmE3q0ESJ3Tk3D7gJKAX2AxPMbF9d+ymxi4jUX6iJPdx57L8ws0wz8wErgZww6xMRkTCFldjN7Fi5h+2BlrHCjohIMxb2PHbn3M+BO4GjwNW1lJsMTAaCswtERCTy6hxjd869BXSp5qlHzOy1cuUeBhLNbFZdL6oxdhGR+gt1jL3OHruZXRvia/4BeB2oM7GLiEjjCWuM3Tl3WbmHNwLbwwtHRETCFe4Y+3zn3LfxT3fcA9wdfkgiIhKOsBK7mf1zpAIREZHI0HrsIiIeo8QuIuIxSuwiIh6jxC4i4jFK7CIiHqPELiLiMUrsIiIeo8QuIuIxUbmCknOuGP83VcOVDByIQD2xoiW1tyW1FdRer4tUe7ubWUpdhaKS2CPFOZcXykpnXtGS2tuS2gpqr9c1dXs1FCMi4jFK7CIiHhPriX1RtANoYi2pvS2praD2el2Ttjemx9hFRKSqWO+xi4hIJTGZ2J1zo5xzO5xznzrnZkY7noZyzi12zu13zhWU29bJOfemc+6TwM+Oge3OOfebQJu3OOf6l9tnfKD8J8658dFoSyicc990zr3tnPvYObfVOTc1sN2TbXbOJTrnNjrnNgfaOyewvYdz7u+B2Jc559oEticEHn8aeD6tXF0PB7bvcM6NjE6L6uaci3PO/cM5tzLw2Mtt3e2c+8g5l++cywtsax7HspnF1A2IAz4DegJtgM1AerTjamBbrgT6AwXlti0AZgbuzwQeC9y/HlgNOGAw8PfA9k7AzsDPjoH7HaPdthraezHQP3A/CSgE0r3a5kDcHQL344G/B9rxMnBbYPvTwE8C9+8Bng7cvw1YFrifHjjOE4AegeM/Ltrtq6HND+C//vHKwGMvt3U3kFxpW7M4lqP+5jTgzRwCrC33+GHg4WjHFUZ70iol9h3AxYH7FwM7AvefAcZWLgeMBZ4pt71CueZ8A14DrmsJbQbaAR8C38X/RZXWge3B4xlYCwwJ3G8dKOcqH+PlyzWnG9AVWAd8D1gZiN2TbQ3EVl1ibxbHciwOxVwKfFHu8d7ANq9INbMvAQI/vxHYXlO7Y/L9CHz07oe/F+vZNgeGJvKB/cCb+HugR8zsXKBI+diD7Qo8fxToTOy09wlgBv5rIIM/dq+2FcCAN5xzm5xzkwPbmsWxHO7FrKPBVbOtJUztqandMfd+OOc6AK8C08zsmHPVNcFftJptMdVmMzsP+JxzFwF/BL5TXbHAz5htr3PuB8B+M9vknMsu21xN0ZhvazlDzWyfc+4bwJvOue21lG3S9sZij30v8M1yj7sC+6IUS2Mocs5dDBD4uT+wvaZ2x9T74ZyLx5/UXzCzFYHNnm4zgJkdAXLxj69e5Jwr61SVjz3YrsDzFwKHiI32DgVudM7tBl7CPxzzBN5sKwBmti/wcz/+f9qDaCbHciwm9g+AywJn29vgP/Hy5yjHFEl/BsrOjI/HPw5dtv3OwNn1wcDRwEe9tcAI51zHwBn4EYFtzY7zd82fBT42s1+Ve8qTbXbOpQR66jjn2gLXAh8DbwO3BIpVbm/Z+3AL8D/mH3j9M3BbYCZJD+AyYGPTtCI0ZvawmXU1szT8f5P/Y2a348G2Ajjn2jvnksru4z8GC2gux3K0T0A08KTF9fhnVHwGPBLteMJox4vAl8BZ/P+578I/zrgO+CTws1OgrAOeDLT5IyCrXD2TgE8Dt4nRblct7R2G/2PmFiA/cLveq20GMoF/BNpbAOQEtvfEn6w+BV4BEgLbEwOPPw0837NcXY8E3ocdwPej3bY62p3N17NiPNnWQLs2B25by/JQczmW9c1TERGPicWhGBERqYUSu4iIxyixi4h4jBK7iIjHKLGLiHiMEruIiMcosYuIeIwSu4iIx/x/8QS+/q8zUQsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_c, ax_c = plt.subplots()\n",
    "ax_c.scatter(range(0,nmpsdm.shape[0]),nmpsdm[:,0], marker=\".\", c=\"k\", s=0.4, label='Embedding trajectories')\n",
    "ax_c.scatter(time_stamps, np.zeros(time_stamps.shape),marker=\"x\",c=\"r\", s=25.0, label='Subsample point')\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.ylim(-3.5, 4.5)\n",
    "plt.savefig(\"pre-print-MDtrj-smpl.pdf\", format=\"pdf\", bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mds2 = MDTrSampler(y, n_dim = 2, conv_size = 50, n_samples=64, batch_size=512, manifold_size=256)\n",
    "for i in range(x):\n",
    "    mds2.adaptive_sampling_step(trace[i,:,:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGiZJREFUeJzt3X2QXFWZx/Hvk5kQXqIJ2YxhYAgB\nAxbEoqJM4jsaN/FtA5ZbqybB9QWUBcqUCKUrldpsIcu6omAULRE3mAIFdHUpEV8gKFtquZJMNLiB\nSBLzIiNDmJSEoEjCJM/+cW47PUl3T7/c7tv39O9T1XW7b9+595w7t5977nNO3zZ3R0RE4jEh6wKI\niEi6FNhFRCKjwC4iEhkFdhGRyCiwi4hERoFdRCQyCuwiIpFRYBcRiYwCu4hIZLqz2Oj06dN91qxZ\nWWxaRCS3NmzYsMfde8ZbLpPAPmvWLAYGBrLYtIhIbpnZrmqWUypGRCQyCuwiIpFRYBcRiYwCu4hI\nZBTYRUQio8AuIhIZBXYRkchkMo49LfPnw1NPwcgIdHfDOefAnXfC3Lnw5z+H+T09sG4dXH45fPe7\nYbnC8iMj8LOfhXUtXQqDg6Pzly+HK64Y3db69fCudx359/VMJ06ERx6pvo7Dw0eup69vtOwiIsUs\ni9887e/v90a/oHTMMfDccykVqIypU2H6dPjTn+CJJ9Jdd1cXnHxy5RPAzp3heTkTJ8Ipp1Q+kUA6\nJ6Napt/6Fsybl+7+EhEwsw3u3j/ecrlssc+fPxrUe3rguONgz54QgAtOOAEOHgyt3YKjjw4t3UIQ\n2r4dDh0afX/atPD4/e/hwAHYuzc8Crq64NRTGwt6IyPhyuDgwRC4q9HTA1OmjK7nj38Mj+efh23b\n6t6NTTN/ftj/kydXt0927Mi6xCJxyWVgX7cuBI/LLoP3v390/kUXwU9/Cp/7HCxeHObdc09Iq7zi\nFSFNc7izzgoBcsWKsetatCgEzULwmTIFNm5Mp/zr18OyZdUFvRtvHK1LsTVr4Oqrx/97aG2LfceO\ncNKq5QpnwoRw5VFpvQ88EE7KIu2g8Bm+9NKxKdt2kdtUjLSvuXPh6aerOxH87ndQ7SH47neXPjmL\npOn88+Huu8u/v359aFgW/Mu/wCc/eeRyhX692bNh7dp0ylZtKkaBXTJ3+umVTwC///1oymzevHDF\nJtIMixbB/ffDwoXlg/GSJfDNb46d99hjY68oL78cPv/50deF/rrCgI16rz6rDewa7iiZ27o1pHDK\nTQ8eDPl6CK2lG27ItrwSp8FB2LQpPN+4Mbwu5bjjwvRrXxs9Lt/0prHL3HdfmHZ1henevSG1u3Nn\naMiUW3daFNglF555Bgq38P/P/8y0KBKpoaFwnL32tWEwxv33H7nMmjVwyy1w1FEwZw7ccUeYv3lz\n6M8DWLkyvD7zzNBCX7gwzJ86NRzDW7c2v79IgV1y48Ybw3Tz5pDHF0nT7t3w7LOwYEF4PX36kcsU\n5n3qUyEtuHgxnHdemHfBBWH61a+GaaEVv3ZtyMM/9VS4Am3FIAAFdsmNxYtHx8c/9JCCu6TPPbSs\njzkGZswov9wZZ4w+P+20MN23D8xGR4QVH5+lOlebSYFdcmXdujBGHkJwX7ky2/JIPGbMCN91mTYt\n9OuUsmVLmO7ZMzpv1Sr4yEdCUC+48MKxw6dbTYFdcmdoaDS433ZbtmWRuEyYEDpODxyAm24a+97g\nIKxeHYL/nDlj31u1KozcmjUrdKquXt2yIpekwC65NDQUPkQ7d4ahZSKN+t73wvSCC0KLe82a0Q5R\nCMfcrl3w5S+Xv2XGjh3ZttQLFNglt5YvD9OZM7Mth+TfmjVwzTWj35e45JIw8qU4z17oXC3Vqdpu\nFNglt4o7sETS8LGPhdb47t3hflS7d499P4Pvc9ZFgV1ya8aMMHphxYrwxSWRem3fPvZ1oXO0uJO0\n0LlaabRMu1Bgl9yaNw/+7d9Kt6xEanHaaTBp0uiY9OnTx45yKZiQk4iZk2KKlKZ0jDRq/Xr46EfD\nKJhCp+jcueGbo5/5zOjX/wudq3mgwC65V6plJVKt3l448cTRr/5D+Hboxz4WvuW8ceORnavtToFd\ncm3GjHAJnYe8p7SnoaHQKh8aGju/ePRL4fm11+bj18EU2CXXdu8OP5Si33+VevX2hhx7b+/Y+YVG\nw5YtozcEy0vqT4Fdcm3xYrjuuvBrUhoZI/U66qgj5xWGPF555ei91fPyWwANB3YzO9nMHjCzzWb2\nsJl9JI2CiVTrda+r/KPfIvUo3Bem4IQTWn8zr3ql8ZunI8CV7v4rM3sBsMHM1rr7IymsW0QkE4Xf\nMi38HGNeWuuQQovd3Yfc/VfJ82eAzcBJja5XRCRrV1wRAnqegjqknGM3s1nAy4AHS7x3sZkNmNnA\n8PBwmpsVEWnIpElZlyBdqQV2M5sMfAe43N33Hf6+u9/s7v3u3t/T05PWZkWA/IwvlvbT1wfXX9+a\nXzZqlVQCu5lNJAT1b7j7f6exTpFq7d4N+/frtgJSn8HBMPKl2T8w3UppjIoxYDWw2d31+/HScjNm\nhF+LL75hk0gt9u/PugTpSqPF/hrgH4E3mtnG5PG2FNYrUpV58+ALXwg/T6ax7FKP2HLsDQ93dPef\nA7pbh2RqzhyNZZf6KMcu0sbycktVaS/KsYu0qd7e8BN5h9/vQ2Q8Q0PKsYu0pXJ36BOpZHAQPvzh\nrEuRvjRuKSCSuXJ36BMZz4EDcPPNyrGLtKVSd+gTqWRoCHbuzLoU6VNgF5GONW8e3HZbPn48oxYK\n7CLSsQYHw/38YxoRAwrsIiLRUWAXkY4V45eTQIFdIhLb18Kl+QYHYcUKpWJE2lJfH9x+e9alkLzp\n64NbblGLXaStLVsWX+tLpFYK7BKV2L4aLs01OAgXXhhfY0CBXUQ6Vl8fXHutUjEibU0dqFILdZ6K\ntLlYh65J86jFLtLmYryvtjSXWuwiIpHRcEeRNqdUjEigwC7RiPWyWppHwx1F2lysHWHSPErFiLQ5\ndZ6KBArsItKxlIoRaXPqPJVaKRUj0uaUihEJFNhFpGMpFSPS5pSKkVopFVOBmd1iZk+a2aY01idS\nD6ViRIK0WuxrgLektK7KRkbgttvAPUxHRlqyWQD27oXTT4cTT4TZs2HqVBgebk05Rkbgy1+GM86A\nM88M03PPhT/9qbnbBXjuObjkEti3D046Cf78Z3jjG+vfdq3/w+eeC9fLCxbAhz4Er3992HZhPc8/\nDx/8IFM//iHe/PgtTH/nAvjLXxr/v4yMwOrV8IEPhH1+0UWj226m4npdckmYtvpYLy5HvZ+1Rv5+\nZATWrIFbbw2PNWvqq//h6/niF8NnePZs6OpiZNp0ppwyhT9e86X6t9GO3D2VBzAL2FTNsuecc47X\n7dZb3cOhEh633lr/ump15pljt93Kchxe78JjwYLmbtfd/Z/+Kd1t1/o/LLX9BQvK75O0/i9Z7fNy\n223lsV6qHLVuv5G/L7UP6ql/NcdIlvu4RsCAVxFjLSzbODObBdzj7i8t8/7FwMUAM2fOPGfXrl31\nbcgdJhRdaBw6BGb1ratWIyMwcWLp95pdjsPrXVymrq7mbRdC3Upto95t1/o/LLX9kZGwjlL75PC/\nrff/ktU+L7fdVh7rpcpR6/Yb+ftS+6Ce+pfbl+W0eh/XyMw2uHv/eMu1rPPU3W9293537+/p6al/\nRV//euXXzXT22eXfa3Y5yq1/0aLmbhfgssvS3Xat/8NS21+0qLp93sj/Jat9Xm67rTzWS22v1u03\n8vellq2n/o2WOae6sy5AzZYuDdP3vCf8EwqvW+EXv4B580KO+dhjYc8e2LoVfvSj5pdj6dKQ2/3c\n50Jr8eBBOOEEuPvu5m4XYNWqML3uupBr3rIFzjuv/m3X+j9ctQoOHIAdO0JudMuWsO2jjw7vL1kC\nl14KZmw/4VWc9vPb4Ac/gG9/u7H/y9KlYbs//zn88pfw6lfDtm3N3+eFMi9ZAsuXw403wp13tvZY\nLy5HvZ+1Rv5+6dJwjBda24cO1Vf/w9ezbx98/vOhJb9jBxx/PAf3P0/Xf/w7HHdc6/dxk7QsFVOs\nv7/fBwYGUtmuSMHgICxbBrffHt/wNRFocSrGzO4A/hd4iZkNmtlFaaxXRKTZYhwem0pgd/el7t7r\n7hPdvc/dV6exXpFa9PWF1rpItfTNU5GcWLYsvg+qNIe+eSqSE/v3Z10CyZPYgjoosEuEJk3KugSS\nF7Fe2SmwS1R0IzCpVqz5dVBgl8joRmBSrVjz66DALhFSjl06nQK7REc5dqmGUjEiOaEcu1RLqRiR\nnBgchBUr4myFSfpiDOqgwC6R6euDa6+N9wMr6Yq1AaDALlHRqBiplnLsIiKRUY5dJCfUeSq1iPU4\nUWCXqKjzVGoR63GiwC5RUeepVEs5dpGcUItdqqUcu0hOqMUuosAukVGLXaqlVIxITqjFLtVSKkYk\nJ9Ril1rEGNRBgV0ioxa71CLWBoACu0RFLXaplnLsIjmhFrtUSzl2kZxQi11qEWNQBwV2iUzMrTBJ\nX6wNAAV2EelIyrGL5ETMH1ZJV8xXd6kEdjN7i5k9ambbzOwTaaxTpB6FD6tIJ2s4sJtZF/Al4K3A\nWcBSMzur0fWKNEKtdhlPzFd3abTY5wPb3H27ux8A7gTensJ6ReqiIY9SDaViKjsJeKzo9WAybwwz\nu9jMBsxsYHh4OIXNipSmIY9SrRiDOqQT2K3EPD9ihvvN7t7v7v09PT0pbFakNLXYpVqxnvzTCOyD\nwMlFr/uAx1NYr0hd1GKXaijHXtl64HQzO9XMjgKWAHensF6RusScO5X0xHycdDe6AncfMbMPA/cC\nXcAt7v5wwyUTEZG6pDKO3d1/4O5nuPuL3f3aNNYpUq+YL7ElPTEfJ/rmqURHnadSjZhTMQrsEh11\nnkq1YgzqoMAuEYq5JSbpivXkr8AuIh1JOXaRHIn5AyvpifnKToFdoqPOU6lWrMeIArtER52nUq1Y\njxEFdolOzJfYkp6YU3YK7CLSkWJuACiwS3RibolJumIM6qDALhFS56lUK9aTvwK7REedp1KNmK/s\nFNglOjHnTiU9MR8nCuwiIpFRYJfoxHyJLemJ+ThRYJfoqPNUqqFUjEiOqPNUqhVjUAcFdomQWuxS\nrVhP/grsEh212KUayrGL5Iha7FIN5dhFckQtdqlWjEEdFNglQmqxS7ViPfkrsEt01GKXaijHLpIj\narFLNZRjF8kRtdil0ymwS3RibolJepSKERGJTMwNgIYCu5m908weNrNDZtafVqFEGhFzS0zSMzgY\nZ1CHxlvsm4C/B36aQllEUqHOUxlP7Cf/hgK7u29290fTKoxIGtR5KuOJOQ0DLcyxm9nFZjZgZgPD\nw8Ot2qx0oMKHVqSSWIM6VBHYzex+M9tU4vH2Wjbk7je7e7+79/f09NRfYpEqxXypLY2L+djoHm8B\nd1/YioKIpEl5dqmkkGOPNR2j4Y4SpcFBuPLKuFtlUj/l2Csws3eY2SDwKuD7ZnZvOsUSEZF6NToq\n5i5373P3Se4+w93fnFbBRBrR1wfXXx9vi0wao+GOIjmkIY9SiVIxIjmkzlPpZArsEiV1nkolSsWI\niERGqRiRHFLnqXQyBXaJklIxUolSMSIikVEqRiSHlIqRTqbALlFSKkYqUSpGRCQySsWI5JBSMdLJ\nFNglSrqlgFSiVIxIDumWAlKJUjEiOaTOU+lkCuwi0nGUihHJIXWeSiVKxYjkkDpPpZMpsEuU1Hkq\nlSgVI5JDarFLJUrFiOSQWuzSyRTYJUpqsUslSsWI5FDhUluklNiv6BTYJWoxt8qkfrFf0SmwS7Ri\nb5VJ/WI/NhTYJVqxt8qkfrEfGwrsEq3YW2VSPw13FMkp3QhMOlVDgd3MPmNmvzWz35jZXWY2Na2C\niYg0i4Y7VrYWeKm7nw1sAa5qvEgi6dCNwKSc2NN0DQV2d7/P3UeSl78EIt1Nkkexd5BJ/WI/NtLM\nsV8I/LDcm2Z2sZkNmNnA8PBwipsVKS32VpnUr+M7T83sfjPbVOLx9qJlVgAjwDfKrcfdb3b3fnfv\n7+npSaf0IhWo81Q6Vfd4C7j7wkrvm9n7gMXA37q7p1UwEZFmKXSextpqHzewV2JmbwH+GXi9uz+b\nTpFE0qHOUykn9jRdozn2LwIvANaa2UYzuymFMomkQqkYKSf2ztOGWuzuPjutgoiItErHd56K5JVS\nMVLO4GDcx4UCu0RLqRgpJfZvnYICu4h0mNjTMKDALhFTKkY6lQK7REupGClFqRiRnNu/P+sSSLtR\nKkZERHJHgV1EOopSMSIikVEqRiTnJk2CoaGsSyHSWgrsEq3CcEeNjJFiSsWI5FxvLzzzTNalkHYS\n+50dQYFdIjc0BNu3Kx0jo2K/syMosEvkentDy6y3N+uSSLtQ56lIzg0Nwc6darFLZ1Fgl6jt3g3P\nPhumIhBSMMuWKRUjkmtmWZdApLUU2CVqM2bAscfCnj1Zl0TaRSfc9bOhn8YTaXe9vXDKKbB8OcyZ\nA/PmZV2i9rF+fUhJjIxAd/f4061bS69nyRLYsCEsB2HZiRPhkUdaV5daFO76efvt8QZ3BXaJWl8f\nfPrT8O53Z10SWLkS7rij+kDazOnISO055gkT4MUvHrueJ58MfRildHfDySfD7Nmwdm3pZebPh6ee\nCuvasaO28kh5CuwSvRkz4NCh8u/Pnw/Dw6PBCsLzc8+F1avD60WLYNu2+gPp3r3h0W66ukLwHa/8\nu3aBe9gHpUyeDNOnh+eFk8bBg2FE0s6d4dYOM2eG8eNz5sC73gV/+AM8//zoOiZMgM9+Fq64onx5\ni68Ourth6VL45Cdrq3MnpGLM3Vu+0f7+fh8YGGj5dqUz3XMPnHcefO97sHjx2PnveMdoMC9l8mQ4\ncCA80jBpUkgPZd1iryddcuqppddTfAIsKKR59uwZ/4TW0xOWK4SiadPghS88clvlrg6mTg2P4uX/\n7u9g1aqxy82dC08/Ha4gnnkGvv3t/AV3M9vg7v3jLqfALrFbuRKuuSZ8iB97LMx70YtCK72gpwem\nTBkN8qWCyNSpoVVabyCdNat8SiJmN9wAN944toXe1RVOFCtWwPvfH+ZNmlTdCbRwdTA0VPmHVF74\nQjjppHDyOv74sScYM3jwwfz1uSiwiyTWrw/pFgit1JGR0dZhdzfcddfYlnxB8WX/8uWVUwRSnfnz\nQ/DeuLH0+697HTzxRPmT4+FXB4sWhVRP8XJPPAHPPVd6/cceO3rCvv76/P1PFdhFihRa7QVm0N8P\n69ZlVyZpnvPPh1//emwH8cKF4Ypp5cpw9ZW3oA4K7CJHWLMGrr668tA9iUsh13/OOXDnnVmXpnEK\n7CIikak2sDf0zVMzu8bMfmNmG83sPjM7sZH1iYhI4xq9pcBn3P1sd58L3AOsTKFMIiLSgIYCu7vv\nK3p5HND6vI6IiIzR8DdPzexa4L3A08CCCstdDFwMMHPmzEY3KyIiZYzbeWpm9wMnlHhrhbt/t2i5\nq4Cj3f1fx9uoOk9FRGpXbefpuC12d19Y5TZvB74PjBvYRUSkeRodFXN60cvzgd82VhwREWlUQ+PY\nzew7wEuAQ8Au4BJ3/0MVfzecLN+I6UAn/XxCp9UXOq/OnVZfUJ1rdYq794y3UCZfUEqDmQ1Uk2uK\nRafVFzqvzp1WX1Cdm0U/jSciEhkFdhGRyOQ5sN+cdQFarNPqC51X506rL6jOTZHbHLuIiJSW5xa7\niIiUkLvAbmZvMbNHzWybmX0i6/I0wsxuMbMnzWxT0bxpZrbWzLYm0+OT+WZmX0jq/Rsze3nR37wv\nWX6rmb0vi7pUw8xONrMHzGyzmT1sZh9J5sdc56PNbJ2ZPZTU+epk/qlm9mBS/m+a2VHJ/EnJ623J\n+7OK1nVVMv9RM3tzNjWqjpl1mdmvzeye5HXs9d1pZv+X3Ol2IJmX3XHt7rl5AF3A74DTgKOAh4Cz\nsi5XA/U5F3g5sKlo3nXAJ5LnnwA+nTx/G/BDwIBXAg8m86cB25Pp8cnz47OuW5n69gIvT56/ANgC\nnBV5nQ2YnDyfCDyY1OVbwJJk/k3Apcnzy4CbkudLgG8mz89KjvdJwKnJ56Ar6/pVqPcVhG+j35O8\njr2+O4Hph83L7LjOfIfUuPNeBdxb9Poq4Kqsy9VgnWYdFtgfBXqT573Ao8nzrwBLD18OWAp8pWj+\nmOXa+QF8F1jUKXUGjgV+BbyC8AWV7mT+X49r4F7gVcnz7mQ5O/xYL16u3R5AH/Bj4I2E23lbzPVN\nylcqsGd2XOctFXMS8FjR68FkXkxmuPsQQDJ9UTK/XN1zuU+SS+6XEVqwUdc5SUtsBJ4E1hJan3vd\nfSRZpLj8f61b8v7TwN+QrzqvAj5O+EY6hPLHXF8Ityy/z8w2JHeyhQyP64Zv29tiVmJepwzrKVf3\n3O0TM5sMfAe43N33mZWqQli0xLzc1dndDwJzzWwqcBdwZqnFkmmu62xmi4En3X2Dmb2hMLvEolHU\nt8hr3P1xM3sRsNbMKt03q+l1zluLfRA4ueh1H/B4RmVplt1m1guQTJ9M5pere672iZlNJAT1b7j7\nfyezo65zgbvvBf6HkFedamaFhlVx+f9at+T9KcAfyU+dXwOcb2Y7gTsJ6ZhVxFtfANz98WT6JOHk\nPZ8Mj+u8Bfb1wOlJD/tRhM6WuzMuU9ruBgq94e8j5KEL89+b9Ki/Eng6uby7F3iTmR2f9Lq/KZnX\ndiw0zVcDm939hqK3Yq5zT9JSx8yOARYCm4EHgH9IFju8zoV98Q/ATzwkXO8GliSjSE4FTgfWtaYW\n1XP3q9y9z91nET6fP3H3C4i0vgBmdpyZvaDwnHA8biLL4zrrToc6OineRhhN8TvCj31kXqYG6nIH\nMAQ8TzhbX0TIL/4Y2JpMpyXLGvClpN7/B/QXredCYFvy+EDW9apQ39cSLi1/A2xMHm+LvM5nA79O\n6rwJWJnMP40QqLYB/wVMSuYfnbzelrx/WtG6ViT74lHgrVnXrYq6v4HRUTHR1jep20PJ4+FCXMry\nuNY3T0VEIpO3VIyIiIxDgV1EJDIK7CIikVFgFxGJjAK7iEhkFNhFRCKjwC4iEhkFdhGRyPw/yjoH\n6XeqkQAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbb626995f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "adaptive_samples = sorted([ mds2.strm_smplr.reservoir[i][1][0] for i in range(64)])\n",
    "fig_d, ax_d = plt.subplots()\n",
    "ax_d.scatter(range(0,nmpsdm.shape[0]),nmpsdm[:,0], marker=\".\", c=\"b\", s=0.5)\n",
    "ax_d.scatter(adaptive_samples, np.ones(len(adaptive_samples)),marker=\"x\",c=\"r\", s=10.0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0,\n",
      " 514,\n",
      " 516,\n",
      " 517,\n",
      " 522,\n",
      " 533,\n",
      " 535,\n",
      " 538,\n",
      " 539,\n",
      " 540,\n",
      " 558,\n",
      " 1025,\n",
      " 1038,\n",
      " 1045,\n",
      " 1046,\n",
      " 1048,\n",
      " 1056,\n",
      " 1059,\n",
      " 1068,\n",
      " 1537,\n",
      " 1544,\n",
      " 1570,\n",
      " 1572,\n",
      " 1576,\n",
      " 1871,\n",
      " 2060,\n",
      " 2067,\n",
      " 2077,\n",
      " 2084,\n",
      " 2090,\n",
      " 2171,\n",
      " 2192,\n",
      " 2217,\n",
      " 2224,\n",
      " 2250,\n",
      " 2560,\n",
      " 2572,\n",
      " 2577,\n",
      " 2591,\n",
      " 2596,\n",
      " 2598,\n",
      " 3064,\n",
      " 3092,\n",
      " 3622,\n",
      " 3880,\n",
      " 4129,\n",
      " 4145,\n",
      " 4435,\n",
      " 4442,\n",
      " 4486,\n",
      " 4493,\n",
      " 4498,\n",
      " 4501,\n",
      " 4505,\n",
      " 4512,\n",
      " 4513,\n",
      " 4518,\n",
      " 4524,\n",
      " 4528,\n",
      " 4530,\n",
      " 4533,\n",
      " 4535,\n",
      " 4541,\n",
      " 4555]\n"
     ]
    }
   ],
   "source": [
    "pprint(adaptive_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}