Earth Mover's Distance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P_r:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAA+xJREFUeJzt3dFNxDAQQEEbXf8tmwoQcKdjyctMA87X00rZOPucswDo+ph+AADeS+gB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiDuMf0A8JW915k495y1J86FdzHRA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1A3GP6AeC/2XudqbPPWXvqbLpM9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPEHfJKxB8og7wcyZ6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIu+SPR6Bq6qc6fqjTZqIHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4VyD8kk/Ugasx0QPECT1AnNADxAk9QJzQA8TZugHGtsnWslH2F0z0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8Q5wqEi/CJOvAsEz1AnNADxAk9QJzQA8QJPUCcrRtg1NRG2Z22yUz0AHFCDxAn9ABxQg8QJ/QAcbZugFu60/1RJnqAOKEHiBN6gDihB4jzMpZvTb60Al5nogeIE3qAOKEHiBN6gDihB4jb51ioACgz0QPECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8R9As2BJAbtmiimAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P_t:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAA+9JREFUeJzt3UFqw0AQAMFs8L+NXz65B3JIQEzUrvrA6tQMaNg9M/MBQNfn9gcAcC2hB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoC4x/YH8P+d15mNc+c5Z+NcqDHRA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxD22PwB+cl5nNs6d55yNc+EqJnqAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4d93cxNa9L8D9megB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gzqVm8M3mBXLznLN1Nl0meoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDgvTME/svW6lZet2kz0AHFCDxAn9ABxQg8QJ/QAcbZufmlrKwLgr0z0AHFCDxAn9ABxQg8QJ/QAcbZugNVtMvfsXM9EDxAn9ABxQg8QJ/QAcUIPEGfrBljlVa3rmegB4oQeIE7oAeKEHiBO6AHihB4gznol8Jbe6SI3Ez1AnNADxAk9QJzQA8QJPUDcLbduNv+WA9yNiR4gTugB4oQeIE7oAeKEHiDuzFhgASgz0QPECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8R9AYdZKIsiKnskAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "l = 10\n",
    "\n",
    "P_r = np.array([12,7,4,1,19,14,9,6,3,2])\n",
    "P_t = np.array([1,5,11,17,13,9,6,4,3,2])\n",
    "P_r = P_r / np.sum(P_r)\n",
    "P_t = P_t / np.sum(P_t)\n",
    "\n",
    "plt.bar(range(l), P_r, 1, color='blue', alpha=1)\n",
    "plt.axis('off')\n",
    "plt.ylim(0, 0.5)\n",
    "# plt.savefig(\"discrete_p_r.svg\")\n",
    "print(\"P_r:\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(range(l), P_t, 1, color='green', alpha=1)\n",
    "plt.axis('off')\n",
    "plt.ylim(0, 0.5)\n",
    "# plt.savefig(\"discrete_p_t.svg\")\n",
    "print(\"P_t:\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A: \n",
      " [[1. 1. 1. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " ...\n",
      " [0. 0. 0. ... 1. 0. 0.]\n",
      " [0. 0. 0. ... 0. 1. 0.]\n",
      " [0. 0. 0. ... 0. 0. 1.]] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "D = np.ndarray(shape=(l, l))\n",
    "\n",
    "for i in range(l):\n",
    "\tfor j in range(l):\n",
    "\t\tD[i,j] = abs(range(l)[i] - range(l)[j])\n",
    "\n",
    "A_r = np.zeros((l, l, l))\n",
    "A_t = np.zeros((l, l, l))\n",
    "\n",
    "for i in range(l):\n",
    "\tfor j in range(l):\n",
    "\t\tA_r[i, i, j] = 1\n",
    "\t\tA_t[i, j, i] = 1\n",
    "\t\t\n",
    "A = np.concatenate((A_r.reshape((l, l**2)), A_t.reshape((l, l**2))), axis=0)\n",
    "print(\"A: \\n\", A, \"\\n\")\n",
    "\n",
    "b = np.concatenate((P_r, P_t), axis=0)\n",
    "c = D.reshape((l**2))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/dist-packages/scipy/optimize/_linprog_util.py:704: OptimizeWarning: A_eq does not appear to be of full row rank. To improve performance, check the problem formulation for redundant equality constraints.\n",
      "  warn(redundancy_warning, OptimizeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EMD:  0.7062374245472838 \n",
      "\n",
      "Gamma:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAA8hJREFUeJzt3bGJlFEYhtH7yzZgqggWITZgAWIXBqZmxuaiYBOCWIAFCGJovIHahOA1MdhgzfZjmWfOCSd4meThwgTfHHvvBXTcue0vANwsUUOMqCFG1BAjaoi5mBh9cBwjP6n/mBhda30Y2Pw2sLnWWq+Hdjk9e+/jus+91BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xx8R/aR1D10TfToyutd4PbH4f2ISrXBOFMyFqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTEndU10youBzXcDm3CVa6JwJkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTGihhhRQ4zDg0PuD+3+HNrl9Dg8CGdC1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiLm77C1Q9H9p9NXD9da21juPaw5ScIC81xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTGihhhRQ4xrokPuDu1OXf389Whkdt37OrPL/3mpIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmKOvffNjx7HzY+emCdDu5+Hdqd8Gtp9OrR7Svbe156W9VJDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xLgmylprrYdDu5dDuy8HNt8MbK611p+h3d+uicJ5EDXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMQ4Pwj/PhnY/Du1uhwfhPIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTGihhhRQ4xrojDs8dDuF9dE4TyIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUPMyDVR4PZ4qSFG1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooaYv40rTsz3JZkeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAA9VJREFUeJzt3EFKm2EUhtGb0m2I8+g8iIuyXYqb61y6kHQDOijkKv/jOcMMLkF8+CCD93S9Xgfo+PHVXwC4LVFDjKghRtQQI2qI+blx9HQ6rfykft44unT3YeHmzN7f4PFu5+75snN3nu9vf/PycvubMzNPv3fuzpze+9BLDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BCzsia6tXj5Z+kuM/P3q7/A/znP28LV14Wbiz5YKfVSQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcQcak10i5XSRQdaKd1ZKJ1ZWym1Jgrfg6ghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiFkZHnzYOHowxgyXGTSc+fX+x15qiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTGihpiVNdHzxlFmxkrpqgMtlM583JmXGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFmZU308W7j6hxu7fFIrJQuWvq/tSYK34SoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSsrImeLxtXF1kpXWOl9PN5qSFG1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiVtZE5/l+5ex53lburrBQuspK6ce81BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpidoYHLy8rZ2deV64aNGSmM2bopYYYUUOMqCFG1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qI2VkTffq9cnbP7VdKD7VQOmOldNFnr5R6qSFG1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKghRtQQI2qIETXEiBpiTtfr9au/A3BDXmqIETXEiBpiRA0xooYYUUOMqCFG1BAjaogRNcSIGmJEDTGihhhRQ4yoIUbUECNqiBE1xIgaYkQNMaKGGFFDjKgh5h/LBSn87ajzIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import linprog\n",
    "from matplotlib import cm\n",
    "\n",
    "opt_res = linprog(c, A_eq=A, b_eq=b, bounds=[0, None])\n",
    "emd = opt_res.fun\n",
    "gamma = opt_res.x.reshape((l, l))\n",
    "print(\"EMD: \", emd, \"\\n\")\n",
    "\n",
    "plt.imshow(gamma, cmap=cm.gist_heat, interpolation='nearest')\n",
    "plt.axis('off')\n",
    "# plt.savefig(\"transport_plan.svg\")\n",
    "print(\"Gamma:\")\n",
    "plt.show()\n",
    "\n",
    "plt.imshow(D, cmap=cm.gist_heat, interpolation='nearest')\n",
    "plt.axis('off')\n",
    "# plt.savefig(\"distances.svg\")\n",
    "print(\"D:\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dual EMD:  0.7062374245472838\n",
      "f: \n",
      " [ 1.  0. -1. -2. -1.  0.  1.  2.  1.  0.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VPW5x/Hvk4SEhEDCErYECISwJRDEFEEUF1ABFRAh1d5qvdqLXrW1AqLgjoIKgrXV2lK1tbdWDTsoiyAqbiABCSSEJUSWhG3YZZGQ8Nw/5gwOmJBlzpwzM+f7fr3mRTI5mfNzDDyZmXPmI6oKIiJynjC7F0BERPbgACAicigOACIih+IAICJyKA4AIiKH4gAgInIoDgAiIofiACAicigOACIih4qwewEX06RJE01OTrZ7GUREQWPNmjUHVDWhOtsG9ABITk5GTk6O3csgIgoaIrKjutvyKSAiIofiACAicigOACIih+IAICJyKA4AIiKH8nkAiEgrEflURDaKSL6IPFTBNiIifxKRQhFZLyI9fN0vERH5xozDQMsAjFbVtSJSH8AaEVmqqhu9thkIINW4XAbgDeNPIiKyic+PAFR1j6quNT7+AUABgMQLNhsC4F/qthJAvIi08HXfRBT8fjxTjuycXThZWmb3UhzH1NcARCQZwCUAVl3wpUQAu7w+L8bPh4TnNkaKSI6I5LhcLjOXR0QB6MVFmzB25npMWLCx6o3JVKYNABGJBTALwB9U9Vhtb0dVp6tqpqpmJiRU62xmIgpSXxcewD+/3o7E+Gi8v3oXlm/aZ/eSHMWUASAideD+x/9dVZ1dwSYlAFp5fZ5kXEdEDnXsxzN4ZOZ6tGtSDwsfuhKdmtfHo7M24PCJUruX5hhmHAUkAN4CUKCq0yrZbD6AO42jgXoBOKqqe3zdNxEFrwkLNmLP0VOYmpWBuOg6mJqVgSMnS/HEvDy7l+YYZjwC6APgDgDXisg64zJIRO4TkfuMbRYCKAJQCODvAO43Yb9EFKSWbtyHmWuKcf/V7XFJ64YAgLSWcXioXyo+Wr8H83N327xCZ/D5MFBV/RKAVLGNAnjA130RUfA7ePw0xs1ejy4tGuD3/VLP+9p9V6VgWcF+PDk3D5e1bYRmDeratEpn4JnARGQZVcUTc/Nw7FQZXvlld0RGnP9PUER4GKZlZeB0WTkenbUe7t8dyV84AIjIMvPW7caivL0YdX0HdGxev8Jt2iXE4rEBnfDZZhfeX72rwm3IHBwARGSJPUdP4al5echs0xD/c2W7i257Z+9kXJ7SGM99uBE7D560aIXOwwFARH6nqhg7cz3OlCteHpGB8LCLvmyIsDDBlBEZCBfBmBm5KD/Lp4L8gQOAiPzu36t24outBzD+xs5IblKvWt+TGB+Np27ugm+3H8LbX37v5xU6EwcAEfnV9gMnMOmjAvTtkIBfX9a6Rt87/NIkXNelGaZ8vBlb9v3gpxU6FwcAEflN+VnFmBm5qBMumHxrN7jPG60+EcELw7oiNioCo7LX4Uz5WT+t1Jk4AIjIb/7+RRFydhzGhCHpaB5Xu2P6m8RGYdIt6cgrOYY/Ly80eYXOxgFARH6xae8xTPt4CwamN8eQ7i19uq0B6S1wyyWJeP3TQuTuOmLSCokDgIhMV1p2Fg9/kIsG0RF4fmh6jZ/6qcgzg9OQEBuFUdnr8OOZchNWSRwARGS6P32yFQV7juGFYd3QODbKlNuMi66DKSO6YZvrBKYs2WzKbTodBwARmeq7nYfxl88KMcI4gsdMV6Ym4I5ebfD2V9/jm20HTb1tJ+IAICLTnCotx+jsXLSIcx/D7w/jBnVCm0YxGDMjFz/8eMYv+3AKDgAiMs1Lizeh6MAJTBnRDfXr1vHLPmIiIzA1KwN7jp7C8x8W+GUfTsEBQESm8OQd77o8GZenNPHrvi5t0wgj+6bggxxmJH1hVhLybRHZLyIVpnxE5GoROeoVjHnKjP0SUWA49uMZjJmRi3YJ9fDogE6W7PPh61KZkfSRWY8A/glgQBXbfKGq3Y3LBJP2S0QBYMKCjdj3w2lMy+qO6MhwS/YZFRGOaVndmZH0gSkDQFVXADhkxm0RUXD5Ke+Ygu6t4i3dd5eWDfCH/h2YkawlK18D6C0iuSKySETSLNwvEfmJJ++Y1rIBfndtatXf4Af39m2HS1rH48m5edh37Edb1hCsrBoAawG0UdUMAH8GMLeyDUVkpIjkiEiOy+WyaHlEVFOqisfnuPOO07J+nne0SkR4GKaOcGckx85kRrImLPk/pqrHVPW48fFCAHVEpMLDBFR1uqpmqmpmQkKCFcsjolqYu64Ei/P3YvRF8o5W8WQkP9/iwnvfMiNZXZYMABFpLsabgYhIT2O/PI2PKEi58475+EVyQ/y2iryjVe7snYw+7Rvj+Y+Ykawusw4DfQ/ANwA6ikixiNwjIveJyH3GJsMB5IlILoA/AbhN+TiNKCh58o7lZ6uXd7RKWJhgynBmJGsiwowbUdXbq/j6awBeM2NfRGQvT97x+aHpaNO4enlHq7SMj8bTg9MwZkYu3vqyCCP7pti9pIDGM4GJqNq8847/VcO8o1Vu7ZGI67o0w8tLtjAjWQUOACKqlvKzitE+5B2tci4jWZcZyapwABBRtUxfUYQ1Ow7juaG1zztaxZ2R7MqMZBU4AIioSpv2HsMrS7dgUNfmGJzhW97RKgPSm2MYM5IXxQFARBf1U96xDp4f2jVgn/qpyNOD09C0PjOSleEAIKKL8uQdXxzWFY3qRdq9nBqJi66DycPdGcnJi5mRvBAHABFVyjvv2N/kvKNVmJGsHAcAEVXIiryjVcYN6oTkxsxIXogDgIgqZEXe0SrujGR3ZiQvwAFARD/zlZF3/O8+/s87WuXSNg1x71XujOQnBcxIAhwARHSBYz+ewSMW5x2t8of+P2UkDzEjyQFAROfzzjvWrWNN3tEqnozk0VOleHJunuPbARwARHTOx/l7MXNNMR6wIe9olXMZyQ3MSHIAEBEAd95x/JwNSGvZAA/alHe0iicj+dS8fEdnJDkAiChg8o5WiQgPw7Ss7o7PSIb2/2UiqpZAyjtapW2Tehg3sDM+3+LCf77dafdybGFWEextEdkvInmVfF1E5E8iUigi60Wkhxn7JSLfBWLe0Sp39GqDPu0bY+JHBdhx8ITdy7GcWY8A/glgwEW+PhBAqnEZCeANk/ZLRD4I1LyjVZyekTRlAKjqCgCHLrLJEAD/UreVAOJFpIUZ+67I7iOncOJ0mb9unihk/HvlDnyx9QAev7FzwOUdrdIyPhrPDE7D6u2H8daXRXYvx1JWvQaQCGCX1+fFxnWmO3KyFANf/QKTFvJ0b6KL2X7gBCYt3ISrOiTgVz0DM+9olWE9EnG9AzOSAfcisIiMFJEcEclxuVw1/v74mEhkZSbh3VU78fmWmn8/kRN45x1fCuC8o1VEBJOGdUX9uhF4+IN1KC1zRkbSqgFQAqCV1+dJxnU/o6rTVTVTVTMTEhJqtbPR13dEatNYjJ2Zi6Mn+c5/RBcKpryjVZrERmHiLV2Rv/sYXlu+1e7lWMKqATAfwJ3G0UC9ABxV1T3+2lndOuF45ZfdcfB4KZ6eX+GBSUSOVbAn+PKOVjmXkfxsmyMykmYdBvoegG8AdBSRYhG5R0TuE5H7jE0WAigCUAjg7wDuN2O/F5OeGIffXZuKuet2Y+EGv80aoqBSWnYWo7KDM+9oFSdlJM06Cuh2VW2hqnVUNUlV31LVv6rqX42vq6o+oKopqtpVVXPM2G9V7r8mBd2S4vD4nA3Y/4NzT/cm8nj1ky0o2HMML90afHlHq8RF18GU4RmOyEgG3IvAZqoTHoZpWRk4UVqO8bM3OPZ0byIAWLvzMN74bBuyMpPQr3Nw5h2tckVqE9zZ252R/HrbAbuX4zchPQAAoH3T+nh0QCcsK9iPGWuK7V4OkS1OlZZjjJF3fPKm4M47WuWxge6M5CMz1odsRjLkBwAA/PflyejVrhEmLNiI4sMn7V4OkeVCKe9oFe+M5HMfbrR7OX7hiAHgOd0bAMbMyMVZh53uTc7myTve3adtyOQdrXJpm4a476oUZOcUY9nG0MtIOmIAAECrRjF48qbOWFl0CP/8ervdyyGyhCfvmJJQD2MHdLR7OUHpISMj+djs0MtIOmYAAEBWZitc26kpXlq8CYX7j9u9HCK/e3Z+6OYdrRIV4T6v6OipUjwxN7QOJnHUABARvDisK6IjwzF6Ri7Kyp1xujc508f5ezFrrTvvmBGieUerdG7hzkgu3LA3pDKSjhoAANC0QV1MHNoVubuO4I3Pttm9HCK/OHD8NMbNdkbe0SqejOSTc/Ow92honFfkuAEAADd2a4HBGS3x6idbkVdy1O7lEJnKnXfcgB9+dEbe0SqejGRp+Vk8Ois0MpKO/cmYMCQNjepFOuJ0b3KWOd+VYEn+Poy5wTl5R6u0bVIP4weFTkbSsQMgPiYSLw3vhi37juOVpVvsXg6RKXYfOYWn5+ejZ3Ij3HOFs/KOVvn1ZW1wRfsmIZGRdOwAAIBrOjbF7T1bY/oXRVi9/WJBM6LAp6p4dJZz845WCQsTTB7eDeFhwZ+RdPQAAIAnbuyMVg1jMDo7lxlJCmreecfWjWPsXk5IaxkfjWdudmck3/wieDOSjh8A9aIi8PKIDOw6fBITmZGkIPX9gROYuLCAeUcLeTKSUz/egs17gzMj6fgBAAA92zbC/1zZDv9ZtROfbd5v93KIaqT8rGJ09jpEhocx72gh74zkqOzgzEhyABhGXdcBHZrF4tFZ65mRpKDytxXbsHbnEeYdbdAkNgqThgVvRtKsItgAEdksIoUi8lgFX79LRFwiss64/NaM/Zqpbp1wTMtyZySfYkaSgoQn73hj1xbMO9rkhrTmGNbDnZFcF2QZSZ8HgIiEA3gdwEAAXQDcLiIVveH4B6ra3bi86et+/cGTkZy3bjc+Ws+MJAW202XlePiDdYiLjsRzQ9P51I+Nnr45ODOSZjwC6AmgUFWLVLUUwPsAhphwu7a4/5oUZCTF4Ym5zEhSYHt12VZs2vsD844BwJORLHKdwEuLN9m9nGozYwAkAtjl9Xmxcd2FbhWR9SIyU0RaVXZjIjJSRHJEJMflcpmwvJqpEx6GqVndcbK0HONmhdY7/1HoWLPjMP76OfOOgeSK1Cb4Te82+MdX24MmI2nVi8ALACSrajcASwG8U9mGqjpdVTNVNTMhIcGi5Z2vfdNYjB3QCZ9s2o8ZOcxIUmA5WVqGMTOYdwxEjw3sjLZN6gVNRtKMAVACwPs3+iTjunNU9aCqnjY+fRPApSbs16/OZSQ/3Ihdh5iRpMDx0qJN+P7ACbw8IoN5xwATHRmOqVkZQZORNGMArAaQKiJtRSQSwG0A5ntvICItvD4dDCDgz7jyzkg+MpMZSQoMX249gHe+2YG7+7RF75TGdi+HKtCjdfBkJH0eAKpaBuBBAEvg/oc9W1XzRWSCiAw2Nvu9iOSLSC6A3wO4y9f9WqFVoxg8dVMXrCw6hH8wI0k2O3rqDB6ZybxjMAiWjKQprwGo6kJV7aCqKao60bjuKVWdb3w8TlXTVDVDVa9R1aB5mXxEZhL6dWqKycxIks2eXZCP/cw7BoVgyUjyTOAqiAheuLUrYiLDMTp7HTOSZIsl+Xsxe20J845BpHOLBnj4usDOSHIAVEPT+nXx/NCuyC0+ir8wI0kWO3D8NMYz7xiU7u2bgh4BnJHkAKgmT0byT8xIkoXO5R1Pl+GVXzLvGGzCwwRTs7rjTLlibABmJPnTVAOejOTDHwTX6d4UvGavNfKO13dAh2bMOwajtk3qYdygTlixxYV3VwVWRpIDoAbiYyIxeXg3bN1/HNOYkSQ/233kFJ5h3jEkeDKSkxYGVkaSA6CGru7YFL+6rDX+/kURvv2eGUnyj7NnFWNnrke5Mu8YCrwzkqOzAycjyQFQC48Pcmckx8xgRpL849+rduDLQuYdQ0nL+Gg8OzgNOTsCJyPJAVALzEiSP31/4AQmMe8Ykm65JBE3pAVORpIDoJa8M5KfMiNJJikrP4tR2esQFRHOvGMIEhFMusWdkXz4A/szkhwAPjiXkZy5HkdOBu7p3hQ8/raiCN/tPIIJQ9KYdwxRjY2M5MY9x/BnmzOSHAA+8GQkD50oxVPz8u1eDgW5jbuP4Y/LmHd0ghvSmuPWHkn4i80ZSQ4AH6UnxuH3/VIxP5cZSaq902XlGJXNvKOTPD24C5rZnJHkADDB/Vd7ZSSPBd7p3hT4mHd0ngZ162DKCHszkhwAJojwykg+Njtw3/mPAtOaHYeYd3SoPu3tzUhyAJjEk5Fcvmk/snN2Vf0NRHDnHUdnM+/oZN4ZyWMWZyRNGQAiMkBENotIoYg8VsHXo0TkA+Prq0Qk2Yz9BppzGckFzEhS9by4aBO2HzzJvKODnZeRXGBtRtLnASAi4QBeBzAQQBcAt4vIhb/K3APgsKq2B/AKgJd83W8gCgsTvDwiAyKCMTOYkaSL+3LrAfyLeUeCOyP5v1enYMYaazOSZjwC6AmgUFWLVLUUwPsAhlywzRAA7xgfzwTQT0L0MIekhu6M5KrvmZGkyjHvSBd6qF8HdG7RwNKMpBkDIBGA95PexcZ1FW5jNISPAgjZX3k8GcmXFm9C4X77T/emwPPsfOYd6XyREWGYlpWBo6dK8fgcaw4mCbgXgUVkpIjkiEiOy+Wyezm14slI1osMx6jsXJxhRpK8LM7bi9nfMe9IP+fJSB4+WYpTFpwbYMYAKAHQyuvzJOO6CrcRkQgAcQAOVnRjqjpdVTNVNTMhIcGE5dmjaf26mHhLV6wvPoq/fMqMJLkdOH4aj8/ZgPRE5h2pYvf2TcF/ftsLMZERft+XGQNgNYBUEWkrIpEAbgMw/4Jt5gP4jfHxcADL1QEHyw/q2gJDurfEn5dvxYZiZiSdTlUxfrY77zgti3lHqlh4mCDMov6Dzz+BxnP6DwJYAqAAQLaq5ovIBBEZbGz2FoDGIlIIYBSAnx0qGqomDE5H49hIW0/3psAwe20JPt7IvCMFDgnkX8QzMzM1JyfH7mX47LPN+3HXP1ZjZN92GD+os93LIRuUHDmFAa+sQOcWDfDeyF4sfJHfiMgaVc2szrZ8DGoB74zkqqIKX/qgEObOO+Yy70gBhwPAIucykjNzcZwZSUf5v5U78FXhQTxxYxfmHSmgcABYpF5UBKZmZaD48ClM/IgZSacoch3HC4sKcHXHBNzes1XV30BkIQ4AC/0iuRFGXtkO733LjKQTlJWfxegZucw7UsDiALDYw8xIOoZ33rFZA+YdKfBwAFjMOyP5JDOSIYt5RwoGHAA2SE+Mw0P9UrEgdzc+XL/b7uWQyTx5x/gY5h0psHEA2OR/jfeBeWJuHjOSIeaPzDtSkOAAsElEeBimjsjAKWYkQ8qaHYfwt8+34ZeZrXBtJ+YdKbBxANiofdNYPMqMZMg4WVqGUUbe8YmbeMY3BT4OAJvddXkyerdrzIxkCHhh4SbsYN6RgggHgM3CwgRTRnRjRjLIfbHVhf9buQP3XMG8IwUPDoAAkNQwBk/d7M5Ivv3V93Yvh2ro6KkzeGTGerRvGotHbmDekYIHB0CAGHFpEvp3borJSzYzIxlknp2fD9fx05iWlcG8IwUVDoAAISKYNIwZyWCzOG+PO+94TXt0S2LekYILB0AAYUYyuLh+OI3xc/KQntgAv7u2vd3LIaoxnwaAiDQSkaUistX4s2El25WLyDrjcmEukrwM6toCQ5mRDHiqivFzNuC4kXesE87fpSj4+PpT+xiAT1Q1FcAnqDz1eEpVuxuXwZVsQ4ZnB6ejSWwUM5IBbNbaEizduA+PXN+ReUcKWr4OgCEA3jE+fgfAUB9vjwDExdTBS8O7Yev+45j68Wa7l0MXKDlyCs/Oz0fP5Ea4+4q2di+HqNZ8HQDNVHWP8fFeAJWd+15XRHJEZKWIcEhUw1UdEvBfl7XGm19+z4xkADl7VvHIDOYdKTRUOQBEZJmI5FVwGeK9nbrfzKays5jaGJHiXwH4o4ikXGR/I41hkeNyuWry3xJyxjMjGXD+9c12fL2NeUcKDVUOAFXtr6rpFVzmAdgnIi0AwPizwsyVqpYYfxYB+AzAJRfZ33RVzVTVzISEhFr8J4WOelERmMaMZMAoch3Hi4s3Me9IIcPXp4DmA/iN8fFvAMy7cAMRaSgiUcbHTQD0AbDRx/06RmZyI4zsa2QkNzEjaZey8rMYlc28I4UWXwfAiwCuE5GtAPobn0NEMkXkTWObzgByRCQXwKcAXlRVDoAaGHVdB3RsVh+PzmJG0i5/W1GEdbuO4Lmh6cw7UsjwaQCo6kFV7aeqqcZTRYeM63NU9bfGx1+raldVzTD+fMuMhTtJVEQ4pmZlMCNpk/zdR915x27MO1Jo4dkrQYIZSXucLivHqA9yER8TieeHpNu9HCJTcQAEEWYkrffK0q3YvM+dd2zIvCOFGA6AIBIRHoZpWRn48QwzklZYs+MQpq/Yhtt+wbwjhSYOgCCTkvBTRvKD1cxI+suJ0+68Y8v4aDxxUxe7l0PkFxwAQeg3vd0Zyec+ZEbSX15YVICdh9x5x9ioCLuXQ+QXHABBiBlJ/1qxxYV/r9yJu/u0Ra92zDtS6OIACFJJDWPwNDOSpjt68gzGzmTekZyBAyCIDb80Cf07N8PkJZuxdR8zkmZ4ZgHzjuQcHABBTETwwrCuiI2KYEbSBIs27MGc70rwIPOO5BAcAEEuoX4UJg5Nx4aSo3j900K7lxO03HnHDeiaGIcHmXckh+AACAEDjYzka8sLmZGsBVXFuNkbcKK0HNOyMph3JMfgT3qIYEay9mauKcaygn0Ye0NHpDLvSA7CARAi4mLqYDIzkjVWcuQUJizYiJ5tG+HuPsw7krNwAISQvh0S8OtezEhWlyfveFYVU0dkIIx5R3IYDoAQM35QZ7RuFIPRM5iRrMq5vONNXdCqEfOO5DwcACEmJjICU0dkoOTIKUz8iN2dymxzHccLi9x5x9t+wbwjOZNPA0BERohIvoicFZHMi2w3QEQ2i0ihiDzmyz6papnJjXBv3xS89+0uZiQr4Mk7RkeGYzLzjuRgvj4CyAMwDMCKyjYQkXAArwMYCKALgNtFhG+v6GcPX5eKTs2ZkazIXz/fhtxdR/DckHQ0Zd6RHMzXJGSBqlZ1yElPAIWqWqSqpQDeBzDEl/1S1TwZycMnmZH0lr/7KF79ZCtu6tYCNzPvSA5nxWsAiQC837i+2LiuQiIyUkRyRCTH5XL5fXGhLK3lTxnJBbnMSHrnHZ9j3pGo6gEgIstEJK+Ci19+i1fV6aqaqaqZCQkJ/tiFo9x3VQq6t4rHk/OYkZy2dAs27/sBk2/txrwjEaoxAFS1v6qmV3CZV819lADwPswiybiOLBARHoapRkby0VnrHZuRzNl+CNNXFOH2nq1wTaemdi+HKCBY8RTQagCpItJWRCIB3AZgvgX7JUNKQiweG9AJn252OTIjeeJ0GUbPyEVSw2g8fiOPPyDy8PUw0FtEpBhAbwAficgS4/qWIrIQAFS1DMCDAJYAKACQrap8VdJid/ZOxuUpzsxInss7Dmfekcibr0cBzVHVJFWNUtVmqnqDcf1uVR3ktd1CVe2gqimqOtHXRVPNuTOSGQgTwWgHZSQ/N/KO9/Rpi8uYdyQ6D88EdpDE+Gg8dXMXfOuQjKQ775iL9k1jMYZ5R6Kf4QBwGCdlJJ+en4cDx0uZdySqBAeAwzglI7lowx7MXbcbv7uWeUeiynAAOFBC/ShMuiV0M5LeeccHrmHekagyHAAONSC9BW65JDHkMpLMOxJVH/92ONgzg9PQJDYKD4dQRpJ5R6Lq4wBwsLhod0aycP9xvLwk+DOSxYdP4lnmHYmqjQPA4fp2SMAdvdrgra++x8ogzki6847ut7pg3pGoejgACOMGdUKbRjEYE8QZyXe+2Y5vig7iSeYdiaqNA4DcGcmsDOwO0ozkNtdxvLhoE67t1BS/ZN6RqNo4AAgAcGmbRhgZhBlJ77zji8O6Mu9IVAMcAHSOJyM5dtZ6HD4RHBnJNz5j3pGotjgA6BxPRvLIyVI8OS/P7uVUKa+EeUciX3AA0HnSWsbhD/074MP1ewI6I3m6rByjs3PRqB7zjkS1xQFAP3Nv33a4pHVgZyQ9eceXmHckqjUOAPqZiPAwTB0RuBlJ5h2JzOFrEWyEiOSLyFkRybzIdttFZIOIrBORHF/2SdZoF6AZSeYdiczj6yOAPADDAKyoxrbXqGp3Va10UFBgCcSM5KSFzDsSmcXXJGSBqgb/m8hQhQItI/n5FhfeXcW8I5FZrHoNQAF8LCJrRGTkxTYUkZEikiMiOS6Xy6LlUWUS46Px9OA02zOSnrxjKvOORKapcgCIyDIRyavgMqQG+7lCVXsAGAjgARHpW9mGqjpdVTNVNTMhIaEGuyB/ubVHIq7rYm9G8un5eTh4vBTTsroz70hkkioHgKr2V9X0Ci7zqrsTVS0x/twPYA6AnrVfMlnN7ozkQiPv+OC17dE1Kc7SfROFMr8/BSQi9USkvudjANfD/eIxBZEmsT9lJF9bbl1Gcv8PP+Jx5h2J/MLXw0BvEZFiAL0BfCQiS4zrW4rIQmOzZgC+FJFcAN8C+EhVF/uyX7LHuYzkp4VYX3zE7/tTVYxn3pHIb3w9CmiOqiapapSqNlPVG4zrd6vqIOPjIlXNMC5pqjrRjIWTPZ4ZnIaE2CiMys71e0ZyxppiLCvYz7wjkZ/wVyqqkbjoOpgywv8ZyeLDJzFhwUZcxrwjkd9wAFCNXZnq34ykd97xZeYdifyGA4BqxZ8ZSeYdiazBAUC14p2RfP5D8zKShfuZdySyCgcA1ZonI/n+6l1Yvmmfz7dXVn4Wo2cw70hkFQ4A8oknI/norA0+ZyQ9ecfnhzLvSGQFDgDySVREOKZldfc5I+nJO96c0RI3dWPekcgKHADksy4tG5zLSM6vRUbyxzPeecc0P6yQiCrCAUCmOJeRnJuHfTXMSL7ilXeMj2HekcgqHABkCk9G8nRZzTKSq7cfwvQvmHcksgMHAJmmXUIsxg3sjM82u/Ax/9O8AAAE3klEQVR+NTKSJ06XYXQ2845EduEAIFPd0asN+rRvjOerkZGctLAAuw6fxNQR3Zl3JLIBBwCZKixMMGV41RlJT97xt1e0Rc+2jSxeJREBHADkBy2ryEh65x1HX8+8I5FdOADIL7wzklsuyEg+xbwjUUDwNQgzRUQ2ich6EZkjIvGVbDdARDaLSKGIPObLPik4nJ+RXHcuI7lwwx7MY96RKCD4+ghgKYB0Ve0GYAuAcRduICLhAF6HOwjfBcDtIsJDPhzAnZHsirySY3hteeG5vGO3JOYdiQKBT4deqOrHXp+uBDC8gs16AihU1SIAEJH3AQwBYN5bSFLAGpDeHMOMjOTnW1zMOxIFEDP/Ft4NYFEF1ycC8D4ovNi4jhzi6cFpaFo/Cut2HcHYGzqifVPmHYkCQZWPAERkGYDmFXzpcVWdZ2zzOIAyAO/6uiARGQlgJAC0bt3a15ujABAXXQd/u+NSLN+0n3lHogBS5QBQ1f4X+7qI3AXgJgD9tOLz/0sAeJc9kozrKtvfdADTASAzM7N67ydAAa9bUjy6JVV4jAAR2cTXo4AGABgLYLCqVnba52oAqSLSVkQiAdwGYL4v+yUiIt/5+hrAawDqA1gqIutE5K8AICItRWQhAKhqGYAHASwBUAAgW1XzfdwvERH5yNejgCo8lk9VdwMY5PX5QgALfdkXERGZi8fiERE5FAcAEZFDcQAQETkUBwARkUNxABAROZRUt91qBxFxAdhRy29vAuCAicsJZrwvzsf743y8P34SCvdFG1VNqM6GAT0AfCEiOaqaafc6AgHvi/Px/jgf74+fOO2+4FNAREQOxQFARORQoTwAptu9gADC++J8vD/Ox/vjJ466L0L2NQAiIrq4UH4EQEREFxFyA4AB+p+ISCsR+VRENopIvog8ZPea7CYi4SLynYh8aPda7CYi8SIyU0Q2iUiBiPS2e012EpGHjb8neSLynojUtXtN/hZSA4AB+p8pAzBaVbsA6AXgAYffHwDwENxvS07AqwAWq2onABlw8P0iIokAfg8gU1XTAYTD3S4JaSE1AOAVoFfVUgCeAL0jqeoeVV1rfPwD3H/BHdtjFpEkADcCeNPutdhNROIA9AXwFgCoaqmqHrF3VbaLABAtIhEAYgDstnk9fhdqA4AB+kqISDKASwCssncltvoj3AW7s3YvJAC0BeAC8A/jKbE3RaSe3Yuyi6qWAHgZwE4AewAcVdWP7V2V/4XaAKAKiEgsgFkA/qCqx+xejx1E5CYA+1V1jd1rCRARAHoAeENVLwFwAoBjXzMTkYZwP1vQFkBLAPVE5Nf2rsr/Qm0A1ChA7wQiUgfuf/zfVdXZdq/HRn0ADBaR7XA/NXitiPzb3iXZqhhAsap6HhHOhHsgOFV/AN+rqktVzwCYDeBym9fkd6E2ABig9yIiAvdzvAWqOs3u9dhJVcepapKqJsP9c7FcVUP+N7zKqOpeALtEpKNxVT8AG21ckt12AuglIjHG35t+cMCL4j41gQONqpaJiCdAHw7gbYcH6PsAuAPABhFZZ1w33mg0E/0OwLvGL0tFAP7b5vXYRlVXichMAGvhPnruOzjgrGCeCUxE5FCh9hQQERFVEwcAEZFDcQAQETkUBwARkUNxABARORQHABGRQ3EAEBE5FAcAEZFD/T+pK6Ffg+7fCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g: \n",
      " [ 1.  0. -1. -2. -1.  0.  1.  2.  1.  0.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8FPX9x/H3JwnhJhwJZ4Bw3yRqRLxR0QJSQCBRfj1sa0v9qa0WBDm8QEVRDrVaK7W2tvWnhktROQQv8MAakJCLI0SOhCtc4QgkJPn8/tgNLpiQZHd2Znbm/Xw88iDZnex82Qfkk93ZnZeoKoiIyH3CrF4AERFZgwOAiMilOACIiFyKA4CIyKU4AIiIXIoDgIjIpTgAiIhcigOAiMilOACIiFwqwuoFXEx0dLTGxcVZvQwiopCxYcOGQ6oaU5NtbT0A4uLikJqaavUyiIhChojsqum2fAqIiMilOACIiFyKA4CIyKU4AIiIXIoDgIjIpQIeACLSXkQ+FZEsEckUkfsr2UZE5EURyRGRzSJyaaD7JSKiwBjxMtBSABNVdaOINAawQURWq2qWzzZDAXTzflwB4BXvn0REZJGAHwGo6j5V3ej9/ASAbADtLthsJIB/qcd6AE1FpE2g+6bQkX/sND7YvNfqZRCRD0PfCCYicQAuAfDNBVe1A7DH5+s872X7KrmN8QDGA0CHDh2MXB5ZpLi0DHf981ts2X8CAsGt/Tn7iezAsIPAItIIwGIAD6jqcX9vR1UXqGqiqibGxNTo3cxkcy+s2Y4t+0+gXdP6ePjddBw8ccbqJRERDBoAIlIHnh/+b6rqkko2yQfQ3ufrWO9l5HAbdh3FXz/fgdsT2+ON3wxAUUkZpi5Oh6pavTQi1zPiVUAC4O8AslV1XhWbLQPwS++rgQYCKFTVHz39Q85SVFKKiSmb0CaqPh4e3gtdWzbC5CE98fGWg0hJ3VP9DRBRUBlxDOBqAL8AkC4im7yXTQPQAQBU9a8AlgMYBiAHQBGAXxuwX7K5Z1Zswc7DRXjrdwPRuF4dAMCvr4rD6qz9mPl+Fq7qEo32zRtYvEoi9wp4AKjqFwCkmm0UwL2B7otCxxfbD+FfX+/Cb67uhCu7tDh3eViY4Lmx8Rj6wjo8uDANb/1uIMLCLvrPh4iChO8EJsMVnj6LSYvS0CWmISYP6fGj69s3b4BHh/fGN98fwT++2mn+AokIAAcABcGM9zNx8EQx5t+egHp1wivdJikxFjf1bIlnV25BzsETJq+QiAAOADLYyoz9WLIxH/fe0BX9Y5tWuZ2I4Okx/dAgMhwTUtJwtqzcxFUSEcABQAY6dLIY05emo2+7JvjDjV2r3b5l43p4clQ/bM4rxF8+3WHCConIFwcAGUJVMW1JOk4Ul2JecgLqhNfsn9at/dtgRHxb/PmT7UjPKwzyKonIFwcAGWLJxnx8lHUAk27pge6tGtfqe2eO7IMWjSIxIWUTzpwtC9IKiehCHAAUsL3HTuPxZZkYENccv7mmU62/v2mDSMwe0x/bD57EvNXbgrBCIqoMBwAFpLxcMWlRGspUMScpHuF+vqZ/UI+W+J8rOuBv63LxTe5hg1dJRJXhAKCA/Hv9LnyZcxgP39obHVoE9q7e6cN6oX2zBnhwURpOFpcatEIiqgoHAPktt+Aknl6RjUE9YjBuQPvqv6EaDetGYE5SPPKOnsZTH2YbsEIiuhgOAPJLaVk5Ji5MQ92IcMwe0x+ecwIGbkCn5hh/bWe89d/d+HTrQUNuk4gqxwFAfnl1bS6+230MT4zqi1ZN6hl623+6uTu6t2qEhxZtxrGiEkNvm4h+wAFAtZa19zieX7Pt3Gv4jVavTjjmJSfgyKkSPPJepuG3T0QeHABUK8WlZZiQsglNG0TiyZF9g7afvu2i8MebuuH9tL1sCRMFCQcA1crz3rzj7DH90KxhZFD3dc+gLoiPjcLD72bg4HFmJImMZlQS8nUROSgiGVVcP0hECkVkk/fjUSP2S+basOsIXvXmHW/s2Sro+4sID8Pc5AScLinDlCXMSBIZzahHAP8EMKSabdapaoL3Y6ZB+yWTePKOaWjb1JN3NEvXlo3w0JCe+IQZSSLDGTIAVHUtgCNG3BbZ0zMrtmDXkSLMSYo/l3c0y6+uisOVnVtg5vtZ2HOkyNR9EzmZmccArhSRNBFZISJ9TNwvBWjd9oJzeceBnVtU/w0GCwsTPJfkea/BgwvTUF7Op4KIjGDWANgIoKOqxgP4M4B3q9pQRMaLSKqIpBYUFJi0PKpK4emzmLRwM7q2bIRJP/lx3tEssc1+yEi+/uX3lq2DyElMGQCqelxVT3o/Xw6gjohEV7HtAlVNVNXEmJgYM5ZHFzFjWSYKThZjXnJ8lXlHs5zLSK7ayowkkQFMGQAi0lq85woQkQHe/fKUjza3MmM/lnxXfd7RLBUZyYbMSBIZwqiXgb4F4GsAPUQkT0TuEpG7ReRu7yZjAWSISBqAFwHcoXxNn63VNu9olpaN6+Gp25iRJDJChBE3oqrjqrn+JQAvGbEvCj5VxVRv3vHtWuQdzTKsXxuMTPBkJG/s2RL9YqOsXhJRSLLX/2yyhcUb87Ham3fsVsu8o1lmjujLjCRRgDgA6Dz5x05jRgB5R7NENahzLiM596OtVi+HKCRxANA55eWKyQbkHc0yqEdL/OyKDnjti++ZkSTyAwcAnVORd3xkeOB5R7NMY0aSyG8cAATg/LzjHZcHnnc0S8O6EZibXJGRzLJ6OUQhhQOAUFpWjgkpxucdzXJ5XEVGcg8+3cKMJFFNcQAQXl2bi017gpN3NMu5jORiZiSJaooDwOWCnXc0CzOSRLXHAeBiZuUdzdK3XRTuZ0aSqMY4AFzMzLyjWf53UBfEt2/KjCRRDXAAuFRF3vGOy83JO5olIjwMc5PicbqkDA8t3syMJNFFcAC4UFFJKSacyzv2tno5hqvISH66tQDvfMuMJFFVOABc6OnlW7Dbm3dsVNeQ8wHaTkVG8okPmJEkqgoHgMus216Af6+3Lu9oFmYkiarHAeAihUX2yDuaJbZZAzz6U2YkiarCAeAij79vn7yjWZIui8XgXp6M5PYDzEgS+TKqCPa6iBwUkYwqrhcReVFEckRks4hcasR+qeZWZuzD0u/ycZ9N8o5mERHMGs2MJFFljHoE8E8AQy5y/VAA3bwf4wG8YtB+qQYKThRj2tIM9GsXhftslHc0S0VGMj2/EC9/mmP1cohsw5ABoKprARy5yCYjAfxLPdYDaCoibYzYN12cqmLa0nScLC7FvOR42+UdzTKsXxuMSmiLlz7JQXpeodXLIbIFs34atAPg+4LsPO9lFGShkHc0y4wRfRHdqC7+xIwkEQAbHgQWkfEikioiqQUFBVYvJ6Sdyzt2snfe0SxRDepg9tj+yDl4EnNWMSNJZNYAyAfgWxmJ9V72I6q6QFUTVTUxJibGlMU5UXm5YtLCNJSrYm4I5B3Ncn33GPzsig74+5ffYz0zkuRyZg2AZQB+6X010EAAhaq6z6R9u9K/vt6Jr3YcxsPDe6N989DIO5rlXEZyITOS5G5GvQz0LQBfA+ghInkicpeI3C0id3s3WQ4gF0AOgL8BuMeI/VLlcgtO4pmVW0Iu72iWioxk/jFmJMndDDkRjKqOq+Z6BXCvEfuiiwv1vKNZLo9rjvHXdcarn+filt6tcUPPllYvich0tjsITIFxQt7RLBNu7o4erRrjocWbcfQUM5LkPhwADpK5txDPr9mG4SGedzRL3YhwzE2O92YkK30TO5GjcQA4RHFpGSa8k4amDSLxhAPyjmapyEh+sHkf3k9jRpLchQPAIeav3o6tB07g2TH9HZN3NEtFRvKR95iRJHfhAHCADbuOYMFaT96RBzNrLyI8DPOSmZEk9+EACHGnip2ddzRLl5hGmDKUGUlyFw6AEPf0imzH5x3NcueVzEiSu3AAhLC12wrwn/W7cZfD845m8c1ITmRGklyAAyBEFRadxeRFnrzjgy7IO5qlIiP5X2YkyQU4AEKUG/OOZvFkJFsxI0mOxwEQglakuzPvaBYRwdOj+6FR3QhmJMnROABCjCfvmO7avKNZYhrXxVOj+iI9vxAvfcKMJDkTB0AIUVVMXZKOUyVlrs47mmVoRUby0xxszjtm9XKIDMefICFk0YY8rMk+gMk/Yd7RLDNG9EVMo7qYkJLGjCQ5DgdAiMg/dhoz38/y5B2vZt7RLMxIkpNxAISAC/OOYcw7mur67jH4+UBmJMl5jCqCDRGRrSKSIyJTKrn+VyJSICKbvB+/NWK/bsG8o/WmDeuFDs2ZkSRnCXgAiEg4gJcBDAXQG8A4EanspDTvqGqC9+O1QPfrFjsKTuLpFVtwA/OOlmoQGYG5SZ6M5JMfMCNJzmDEI4ABAHJUNVdVSwC8DWCkAbfrehV5x/qRzDvaQaI3I/n2t3vwyZYDVi+HKGBGDIB2AHxPn5jnvexCY0Rks4gsEpEqf5UVkfEikioiqQUFBQYsL3T99fMdSNtzDE+M7IuWzDvawg8ZyXRmJCnkmXUQ+H0AcaraH8BqAG9UtaGqLlDVRFVNjImJMWl59pO5txAvfLwdw/u3wU+Zd7SNuhHhmHd7PI4VMSNJoc+IAZAPwPc3+ljvZeeo6mFVLfZ++RqAywzYr2Mx72hvfdr+kJFcxowkhTAjBsC3ALqJSCcRiQRwB4BlvhuISBufL0cAyDZgv441b/U25h1t7u7ruyChfVM88m4GDjAjSSEq4AGgqqUA7gOwCp4f7CmqmikiM0VkhHezP4pIpoikAfgjgF8Ful+nSt15BAvW5mLcAOYd7SwiPAxzk+NRXMqMJIUusfM/3MTERE1NTbV6GaY5VVyKYS+uQ7kqVtx/HQtfIeAfX36PGe9n4enR/TBuQAerl0MEEdmgqok12ZbvBLaRc3nHscw7hoo7r4zDVV1a4ElmJCkEcQDYxOc+eccrmHcMGZ6MZDzCRDAxJQ1lzEhSCOEAsAFP3jGNeccQ1a5pfU9GcucRvP4FM5IUOjgAbOCxZRk4fLIE85MTmHcMUWO9GcnnPtqKbcxIUojgALDYivR9eHfTXtx3Y1f0i42yejnkp/MzkpuYkaSQwAFgId+84703MO8Y6ioykhn5x5mRpJDAAWAR5h2daWi/NrjtknbMSFJI4E8dizDv6FyPj+jDjCSFBA4AC+QdLcIM5h0dK6p+HTzrzUg+x4wk2RgHgMk8eUfPqQOYd3Su67wZydeZkSQb4wAw2Rtf78TXuYfxCPOOjseMJNkdB4CJdhScxDMrtuDGni1xO/OOjleRkdzLjCTZFAeASXzzjs+M7se8o0t4MpJdmJEkW+IAMMkrnzHv6FZ/urkberZmRpLshwPABBn5zDu6Wd2IcMxN9mQkH2ZGkmzEkAEgIkNEZKuI5IjIlEqurysi73iv/0ZE4ozYbygoLi3DxJQ0NG/IvKObVWQkP2RGkmwk4AEgIuEAXgYwFEBvAONEpPcFm90F4KiqdgUwH8DsQPcbKiryjrOZd3Q9ZiTJbox4BDAAQI6q5qpqCYC3AYy8YJuRAN7wfr4IwE3igqOg3zLvSD6YkSS7MWIAtAOwx+frPO9llW7jbQgXAnB09eRUcSkeXJiG2Gb1Mf3WCx8QkVt1iWmEKUN64rOtBXjrv3uq/wZynXe/y8fUJemmnEbEdgeBRWS8iKSKSGpBQYHVy/HbrOXMO1LlflmRkfwwC7sPMyNJP9hXeBqPvJeB7QdOmHKCSCP2kA/A911Nsd7LKt1GRCIARAGo9P3xqrpAVRNVNTEmJsaA5Znv820FePOb3fjtNcw70o9VZCTDRfDgQmYkyUNVMXnRZpSWKeYkxSPchNPEGDEAvgXQTUQ6iUgkgDsALLtgm2UA7vR+PhbAJ+rQJ0Ar8o7dWjbCxFuYd6TKMSNJF/rPN7uxbvshTL+1F+KiG5qyz4AHgPc5/fsArAKQDSBFVTNFZKaIjPBu9ncALUQkB8AEAD96qahTVOQd5zHvSNUYe1ksbu7NjCQBOw+dwqwPs3Fd9xj87IoOpu3XkCeZVHW5qnZX1S6q+pT3skdVdZn38zOqmqSqXVV1gKrmGrFfu1nOvCPVAjOSBABl5YqJC9NQJ1zw7Jj+pp4mxnYHgUPVwRNnMJ15R6ql6EZ1Mes2T0byz8xIutKCtbnYsOsoZo7si9ZR5p4mhgPAAKqKad684/zbmXek2hnS15ORfPnTHKTtYUbSTbbsP475q7dhaN/WGJlg/mli+JPKAAs35GFN9kFM/kkPdG3JvCPV3g8ZyU3MSLpESWk5/vROGprUj8CTo/pacoZgDoAA5R0twsz3s3AF844UgKj6dfBcUn/sKDjFjKRLvPjxdmTvO46nR/dHi0Z1LVkDB0AAfPOOc5h3pABd2y0GvxjYEa9/+T2+3sGMpJNt3H0Uf/ksB0neV4JZhQMgAP/8inlHMtbUYT3R0ZuRPHHmrNXLoSA4XVKGB1PS0CbK814QK3EA+Cnn4EnMXsm8IxmrQWQE5ibHY1/haTz5QbbVy6EgmL1yC3IPncJzSf3RuF4dS9fCAeCH0rJyTFzIvCMFx2UdPRnJd1KZkXSaL3MO4Z9f7cSvrorDVV2irV4OB4A/KvKOT45i3pGCgxlJ5zl+5iwmLUxD55iGeGhIT6uXA4ADoNYq8o4/jW+L4f2Zd6TgYEbSeWYsy8KBE8WYl5yA+pH2OE0MB0AtnDnrm3fsY/VyyOH6tI3CA4O7MyPpAB9l7sfijXm4Z5CnCmcXHAC1ML8i7zi2P5o2YN6Rgu/313XGJR2YkQxlh08WY9rSdPRp2wR/uLGb1cs5DwdADX278wgWrMvFuAEdcEMP5h3JHBHhYZib5MlITl7EjGSoUVVMX5qB46dLMS85AZER9vqRa6/V2NSp4lJMTKnIO/ayejnkMp29GcnPtzEjGWre3ZSPlZn7MfGW7ujR2n6nieEAqIFZy7Ox52gR5iYlMO9IlvjllXG4uiszkqFk77HTePS9TFwe1wy/vbaz1cupFAdANXzzjgM6Nbd6OeRSYWGC58YyIxkqVBUPLd6MsnLz8o7+CGgAiEhzEVktItu9fzarYrsyEdnk/bgwF2lbzDuSnbRtWh+PjejDjGQI+M/6XVi3/RCmDeuFji3MyTv6I9BHAFMAfKyq3QB8jKpTj6dVNcH7MaKKbWznUW/ecf7tzDuSPYy5tJ0nI7mKGUm72nnoFGYt32J63tEfgQ6AkQDe8H7+BoBRAd6ebSxP34f3Nu3FH27shr7tmHckeziXkazHjKQdlZUrJqRssiTv6I9AB0ArVd3n/Xw/gKrOa1pPRFJFZL2I2H5IVOQd+8dG4Z4buli9HKLzeDKS/ZiRtKFX1+7Axt3H8MQo8/OO/qj2JS0isgZA60qumu77haqqiFR1ZKqjquaLSGcAn4hIuqruqGJ/4wGMB4AOHcx/+OSbd5yXzLwj2dOQvq0x2puRvKlnS8Tb6N2lbpW9z5N3HNavNUbEh8ZpYqr96aaqg1W1byUf7wE4ICJtAMD758EqbiPf+2cugM8AXHKR/S1Q1URVTYyJifHjrxSYirzjQ0N6Mu9ItvbYiD5o2ZgZSTsoKS3HhJQ0RNWPxJOjQucMwYH+ersMwJ3ez+8E8N6FG4hIMxGp6/08GsDVALIC3G9Q7DniyTsO7Nwcv74qzurlEF1UVP06eHasJyP57EpmJK30wsfbkL3vOJ4Z3Q/NG4bOaWICHQDPALhZRLYDGOz9GiKSKCKvebfpBSBVRNIAfArgGVW13QAoL1dMWpQGAHhuLPOOFBqYkbText1H8cpnO5B0WSwGW5h39IfY+dwiiYmJmpqaasq+Xv/ie8z8IAuzx/TD7Zfb+6VbRL6KSkox7IV1OFumWPnAtZZXptykqKQUt774BUpKy21z34vIBlVNrMm2PMKJ8/OOyYnMO1JoYUbSOrNXbMH3Nsk7+sP1A6C0rBwTUzahQWQ4nhkTOgdviHxd1rE5fn89M5Jm+jLnEN74ehd+fbU98o7+cP0A+MtnO5CWV4gnR/VDy8b2f90uUVUeGOzJSE5elI4jzEgGlR3zjv5w9QDIyC/Eix9vx4j4tri1fxurl0MUkLoR4ZiXnIDC0yV45N0MtgOCyDfvGMqniXHtADhztgwTUjahecNIzGTekRyid9smnoxkOjOSwbLKpnlHf7h2AMxfvQ3bDpxk3pEcpyIj+eh7mcxIGuzQyWJMW2LPvKM/XDkAmHckJ4sID8O85ARmJA3myTum48QZe+Yd/RH6f4Naqsg7tm/WAA8z70gO1Sm6IaYO7cWMpIGWfpePVZkHbJt39IfrBsBT3rzjnKR4NGTekRzsFwM7nstI7jp8yurlhLS9x07jsWX2zjv6w1UD4LOtB/F/3+zG767tzLwjOR4zksYoL1dMXmT/vKM/XDMACovO4qHFm9G9VSNMuLm71cshMkXbpvXx+Ig++HbnUfz9i1yrlxOS/vPNLnyRcwjTb7V33tEfrhkAFXnHUH/dLlFtjb60HW7p3QpzVm1jRrKWvj90CrOWZ+P67jH4nwHOO0eYKwbAh5s9ecc/3sS8I7mPiGDW6H5ozIxkrZSVKyambEJkeBhmh0De0R+OHwAHT5zBw++mIz42CvcMYt6R3Cm6UV08xYxkrYRa3tEfjh4Aqoqpi9NRVFKGuckJiGDekVzMNyOZtueY1cuxtay9oZd39IejfyIuTM3Dx1sOYvKQnujaspHVyyGyHDOS1Ssu9ZwmJtTyjv4IaACISJKIZIpIuYhUGSAQkSEislVEckRkSiD7rKk9R4ow8wPmHYl8MSNZvRfWbMeW/SdCLu/oj0AfAWQAGA1gbVUbiEg4gJcBDAXQG8A4Eekd4H4vinlHoqpd2y0Gv7ySGcnKbNh1FH/9fAeSE0Mv7+iPgAaAqmaranW/RgwAkKOquapaAuBtACMD2W91/vHVTqzPPYJHh/dG++YNgrkropA0ZWhPxLVogAcXpuHEmbNWL8cWikpK8eDCNLSJqo9Hhgf1d1TbMOMYQDsAvicjyfNeVikRGS8iqSKSWlBQUOudHSsqwdyPtuKmni2RlBhb+9USuYAnI5mAfYWn8cQHWVYvxxZCPe/oj2pPhiMiawC0ruSq6ar6ntELUtUFABYAnih8bb+/aYNI/PuuK9C+eX1HH7whCtRlHZvh99d3wSuf7cAtvVu74imPqnyxPfTzjv6odgCo6uAA95EPwLe0Huu9LGgu69gsmDdP5BgPDO6GT7ccxJQl6fioYzPHH/SsTOHps5i0KA1dQjzv6A8zngL6FkA3EekkIpEA7gCwzIT9ElE16kaEY/7t7s5Izng/EwcdkHf0R6AvA71NRPIAXAngQxFZ5b28rYgsBwBVLQVwH4BVALIBpKhqZmDLJiKj9Grj3ozkqsz9WLIxH/cO6oL4EM87+kPsPPETExM1NTXV6mUQOV5pWTmSXv0auQWnsOqB6xx76gNfh04W4yfz16J1VD0svedqRxS+AEBENqhqle/L8uWMvzERBcQ3I/nQYudnJFUV05Y4K+/oD3f+rYnoRzpFN8S0YZ6M5P/9d7fVywmqJRvz8VGWs/KO/uAAIKJzfn5FR1zTNRpPfZjt2Izk3mOn8bgD847+4AAgonPCwgTPju2P8DBnZiTP5R1VMTcpwVF5R39wABDRedo2rY/Hf+rMjKRv3rFDC54mhgOAiH7ENyO5db8zMpJOzzv6gwOAiH7kwoxkSWloZyRLy8oxweF5R39wABBRpaIb1cWs0f2Qufc4Xvpku9XLCcira3PxncPzjv7gACCiKv2kT2uMvrQdXv5sR8hmJLP2Hsfza7bh1n5tHJ139AcHABFd1GM/Dd2MpG/e8YlRffnUzwU4AIjooqLq18FzY+NDMiNZkXecPcb5eUd/cAAQUbWu6RaNO70Zya92HLJ6OTWyYdeRc3nHm3q5t3VwMRwARFQjU4b2Qqfohpi0cLPtM5JFJaWYmOKuvKM/OACIqEbqR4ZjbnJ8SGQkn1mxBTsPF2FOUrxr8o7+4AAgohq7tEMz3H19F6Sk5mFN1gGrl1OpL7Yfwr++3oXfXN0JV3ZpYfVybI0DgIhq5f7B3dCzdWNMWZKOI6dKrF7OeXzzjpOH9LB6ObYXaBEsSUQyRaRcRKoMEIjIThFJF5FNIsLCC1EI881IPvxuuq3aAW7OO/oj0EcAGQBGA1hbg21vUNWEmpZqiMi+erVpgj/d3B3L0/fbJiO5MsPdeUd/BDQAVDVbVUPrhcFEZIjfX9cFl3ZoikfezcD+wjOWruXQyWJMX5qOPm2b4L4bu1m6llBi1jEABfCRiGwQkfEX21BExotIqoikFhQUmLQ8Iqqt8DDB3OQEnC1TSzOS5/KOxaWYf7t7847+qPaeEpE1IpJRycfIWuznGlW9FMBQAPeKyHVVbaiqC1Q1UVUTY2JiarELIjJbp+iGmDqsp6UZyYq844O3dEf3Vu7NO/ojoroNVHVwoDtR1XzvnwdFZCmAAajZcQMisrmfX9ERq7MO4KkPs3FN12h0bNHQtH1X5B0HxDXHXde4O+/oj6A/VhKRhiLSuOJzALfAc/CYiBzANyM5McW8jGR5uWLSojSUqWJOUrzr847+CPRloLeJSB6AKwF8KCKrvJe3FZHl3s1aAfhCRNIA/BfAh6q6MpD9EpG9tImqjxkj+iB111G8ts6cjOS/1+/ClzmHmXcMQLVPAV2Mqi4FsLSSy/cCGOb9PBdAfCD7ISL7u+2SdliVuR9zP9qGQT1aokfr4D0fn1twEk+vYN4xUDxcTkSGEBHMui34GcnSsnJMXJiGuhHheHYs846B4AAgIsO0aFQXTwc5I1mRd5w5sg9aNWHeMRAcAERkqFv6tMaYS2Px8mc7sMngjCTzjsbiACAiwz02ojdaGZyRZN7ReBwARGS4JvXq4LmkeOQWnMLslVsMuc3nmXc0HAcAEQXF1V09Gcl/fLkz4Izkhl1H8OrnO3B7YnvmHQ3EAUBEQTNlaC90DjAj6Zt3fHh4L4NX6G4cAEQUNPUjwzEnwIwk847BwwFAREF1aYdm+N9B/mUk120vYN4xiDgAiCjo7r+pO3q1aVKrjGTh6bOYtHAz845BxAFYyqNMAAAFDUlEQVRAREEXGRGGecnxtcpIzliWiYKTzDsGEwcAEZmiV5smmHBzjxplJFdm7MeS7/Jx7w1dmXcMIg4AIjLN+Os647KOzS6akazIO/Zt1wR/uLGrySt0Fw4AIjJNeJhgblI8zpYpJleSkVRVTPXmHeclJ6BOOH9EBRPvXSIyVVx0Q0wb1hNrtxXgzW/Oz0gu3piP1cw7mibQIMxzIrJFRDaLyFIRqfTJOhEZIiJbRSRHRKYEsk8iCn0/H9gR13aLxqzl2dh1+BQAIP/Yacxg3tFUgT4CWA2gr6r2B7ANwNQLNxCRcAAvwxOE7w1gnIj0DnC/RBTCRM7PSJaWlWMy846mC2gAqOpHqlrq/XI9gNhKNhsAIEdVc1W1BMDbAEYGsl8iCn1toupj5khPRnLc39bjy5zDePjW3sw7msjIYwC/AbCiksvbAdjj83We9zIicrlRCe0wpE9rfLvzKAb1iMG4Ae2tXpKrVNsEFpE1AFpXctV0VX3Pu810AKUA3gx0QSIyHsB4AOjQga1PIicTEcwa3Q8doxvgrms68Rz/Jqt2AKjq4ItdLyK/AjAcwE1a+dv78gH4jvVY72VV7W8BgAUAkJiYWP3bBYkopDVvGImpQ3mWTysE+iqgIQAmAxihqkVVbPYtgG4i0klEIgHcAWBZIPslIqLABXoM4CUAjQGsFpFNIvJXABCRtiKyHAC8B4nvA7AKQDaAFFXNDHC/REQUoGqfAroYVa30fdqquhfAMJ+vlwNYHsi+iIjIWHwnMBGRS3EAEBG5FAcAEZFLcQAQEbkUBwARkUtJTdJsVhGRAgC7/Pz2aACHDFxOKON9cT7eH+fj/fEDJ9wXHVU1piYb2noABEJEUlU10ep12AHvi/Px/jgf748fuO2+4FNAREQuxQFARORSTh4AC6xegI3wvjgf74/z8f74gavuC8ceAyAiootz8iMAIiK6CMcNAAbofyAi7UXkUxHJEpFMEbnf6jVZTUTCReQ7EfnA6rVYTUSaisgiEdkiItkicqXVa7KSiPzJ+/8kQ0TeEpF6Vq8p2Bw1ABig/5FSABNVtTeAgQDudfn9AQD3w3NacgJeALBSVXsCiIeL7xcRaQfgjwASVbUvgHB42iWO5qgBAAboz6Oq+1R1o/fzE/D8B3dtj1lEYgHcCuA1q9diNRGJAnAdgL8DgKqWqOoxa1dluQgA9UUkAkADAHstXk/QOW0AMEBfBRGJA3AJgG+sXYmlnoenYFdu9UJsoBOAAgD/8D4l9pqINLR6UVZR1XwAcwDsBrAPQKGqfmTtqoLPaQOAKiEijQAsBvCAqh63ej1WEJHhAA6q6gar12ITEQAuBfCKql4C4BQA1x4zE5Fm8Dxb0AlAWwANReTn1q4q+Jw2AGoVoHcDEakDzw//N1V1idXrsdDVAEaIyE54nhq8UUT+Y+2SLJUHIE9VKx4RLoJnILjVYADfq2qBqp4FsATAVRavKeicNgAYoPchIgLPc7zZqjrP6vVYSVWnqmqsqsbB8+/iE1V1/G94VVHV/QD2iEgP70U3AciycElW2w1goIg08P6/uQkuOCgeUBPYblS1VEQqAvThAF53eYD+agC/AJAuIpu8l03zNpqJ/gDgTe8vS7kAfm3xeiyjqt+IyCIAG+F59dx3cMG7gvlOYCIil3LaU0BERFRDHABERC7FAUBE5FIcAERELsUBQETkUhwAREQuxQFARORSHABERC71/2AfCmYGt8GCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# linprog() can only minimize the cost, because of that\n",
    "# we optimize the negative of the objective. Also, we are\n",
    "# not constrained to nonnegative values.\n",
    "opt_res = linprog(-b, A.T, c, bounds=(None, None))\n",
    "\n",
    "emd = -opt_res.fun\n",
    "f = opt_res.x[0:l]\n",
    "g = opt_res.x[l:]\n",
    "\n",
    "#print(dual_result)\n",
    "print(\"dual EMD: \", emd)\n",
    "\n",
    "print(\"f: \\n\", f)\n",
    "plt.plot(range(l), f)\n",
    "# plt.savefig(\"f_function.svg\")\n",
    "plt.show()\n",
    "\n",
    "\n",
    "print(\"g: \\n\", f)\n",
    "plt.plot(range(l), g)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAA8tJREFUeJzt3dFtwzAMQEGr6P4rszs0SVno3S0g/QQPRGj7zMwDQNfX9gUA2CUEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHHf2xf4lHOe2Th35jkb58Kn+U3dy0QAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHHXPlAGN9p6qIu7mQgA4oQAIE4IAOKEACBOCADihAAgTggA4oQAIE4IAOKEACBOCADihAAgTggA4oQAIE4IAOKEACDOh2mAf23zYzwzz9k6+y+ZCADihAAgTggA4oQAIE4IAOKEACBOCADizszaii5vtrlvXVPZL6fBRAAQJwQAcUIAECcEAHFCABDn7aO8bGuDxpYUvIeJACBOCADihAAgTggA4vxZfBGvPQB+w0QAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcWdmtu8AwCITAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxAkBQJwQAMT9AP4vHO/GpHZoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAA9pJREFUeJzt3QFqw0AMAMFeyb+DX67+oa0R8c584M6QsAjk5MzMFwBd39sXAGCXEADECQFAnBAAxAkBQJwQAMQJAUCcEADECQFAnBAAxL22L8DnO9dZ+Z2Sec/ZOBeexkQAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABD32r4A/Na5zmydPe85W2fDfzMRAMQJAUCcEADECQFAnBAAxNkaepDNLRrgc5kIAOKEACBOCADihAAgTggA4oQAIE4IAOKEACBOCADihAAgTggA4oQAIE4IAOKEACBOCADihAAgTggA4oQAIE4IAOKEACBOCADihAAgTggA4oQAIE4IAOKEACBOCADihAAgTggA4oQAIE4IAOKEACBOCADizsxs3+EW5zorDzbvORvn0rD1ua6qfJ9NBABxQgAQ99q+wF0qIx3AX5kIAOKEACBOCADihAAgTggA4h67NQQ8gw3A+5kIAOKEACBOCADihAAgTggA4oQAIE4IAOKEACDOC2XwQTZfrvLvaM9lIgCIEwKAOCEAiBMCgDghAIgTAoC4M2MjDKDMRAAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABAnBABxQgAQJwQAcUIAECcEAHFCABD3A3j/JW9UkDDbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "emd:  0.7062374245472836\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "plt.bar(range(l), np.multiply(P_r, f), 1, color='blue', alpha=1)\n",
    "plt.axis('off')\n",
    "plt.ylim(-0.5, 0.5)\n",
    "# plt.savefig(\"p_r_f-scaled.svg\")\n",
    "plt.show()\n",
    "\n",
    "plt.bar(range(l), np.multiply(P_t, g), 1, color='green', alpha=1)\n",
    "plt.axis('off')\n",
    "plt.ylim(-0.5, 0.5)\n",
    "# plt.savefig(\"p_t_f-scaled.svg\")\n",
    "plt.show()\n",
    "\n",
    "#check sum\n",
    "emd = np.sum(np.multiply(P_r, f)) + np.sum(np.multiply(P_t, g))\n",
    "print(\"emd: \", emd)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAABMZJREFUeJzt3TFuVUcUgOH7ECsIUFFRpjVIpDRyyRrSUJhsIEJChAKjSOwAjESXHVBQYBEpZUSK9FmBYQtcehorGHP8fn/fBs5Uv0bz5s3drOu6ANB1aXoBAJwtoQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAuMvTC+D8e//u8Tox9+adJ5uJuVBjRw8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8Q5/VKTrSze316CcAp2NEDxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8T5Zyzn1uGH39eJuftXH24m5sJZsaMHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeI89bNlti8fj3y7suyLMunu1OTgW/Bjh4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4vxhihO9/PhxegnAKdjRA8QJPUCc0APEOaOHL7x4dG/sAbn7T19tpmbTZUcPECf0AHFCDxAn9ABxQg8QJ/QAca5Xwhf2D25PLwG+KTt6gDihB4hzdLMlfvjrx7nhP/07Nxs4NTt6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDj36OEcef/u8cjXrW7eeeLLVmFCvyVuHB3NDf91bjRweo5uAOKEHiBO6AHinNFzoufHxyNzf7l2bWQu1Aj9/3TlwX8jtyJuTAwFEhzdAMQJPUCc0APECT1AnB9jOdGln4c+Y/jmw8xciBF6YHnx6N7IbbJlWZb7T195Z+eMOboBiBN6gDihB4gTeoA4oQeIc+sGGHX424ORGz/7B88uzG0fO3qAOKEHiBN6gDihB4gTeoA4oQeIc70SuJD++fNo7CG3nd2973q1044eIE7oAeKEHiBO6AHihB4gbitv3dy6dTj2a/mytzc2GuBr2NEDxAk9QJzQA8Rt5Rk90LGux9NLyBN64EL6++0fY7N3dr/vpY7Nus5dYAHg7DmjB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiPsM5o9ABjqw+kwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "current_bottom = np.zeros(l)\n",
    "r = range(l)\n",
    "\n",
    "for i in r:\n",
    "\tplt.bar(r, gamma[i, r], 1, color=colorMap.to_rgba(i), bottom=current_bottom)\n",
    "\tcurrent_bottom = current_bottom + gamma[i, r]\n",
    "\n",
    "plt.axis('off')\n",
    "plt.ylim(0, 0.5)\n",
    "# plt.savefig(\"earth_move_p_t.svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAABFpJREFUeJzt3bFtFUEUQNFd5A7AFdCBcQEgUldB8E1KZMkCAkBILsE/cBVOERRgSInoAGoYV4AQSN/Dvz6ngdno6kn7ZncdYywAdD2Y/QAA7JbQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1AnNADxAk9QJzQA8QJPUCc0APECT1A3MHsB4Df2f78OGacu3l0vs44F3bFRA8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPECf0AHEHsx8A/jdjXI5ZZ6/r6TrrbLpM9ABxQg8QJ/QAcUIPECf0AHFCDxAn9ABxQg8QJ/QAcUIPELeXn0A4Pt5Ou6J+c7NxRR3YKyZ6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeI28sfj0DV189vp/xU58mzd36oE2aiB4gTeoA4oQeIE3qAOKEHiBN6gDihB4gTeoA4oQeIE3qAOJ9A+EsPz35MuaL+6+KxK+rAPzHRA8QJPUCc0APECT1AnNADxNm6AZbL1y+mbJMty7KcfriyUbZjJnqAOKEHiBN6gDihB4gTeoA4oQeIE3qAOKEHiBN6gDihB4jzCYQ9sV5fT7uiPk5OXFGHPWaiB4gTeoA4oQeIE3qAOKEHiLN1A0y1fXM2ZaNs8/7i3myTmegB4oQeIE7oAeKEHiBO6AHibN0A99K3L5+mfT/q6OnzO934MdEDxAk9QJzQA8QJPUCcl7H80dH3V1NeWr08PJxxLOSY6AHihB4gTugB4oQeIE7oAeLWMabdAgbgDpjoAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4oQeIE7oAeKEHiBO6AHihB4gTugB4m4BS+UvEEZAUWkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.colors as colors\n",
    "\n",
    "cNorm = colors.Normalize(vmin=0, vmax=l)\n",
    "colorMap = cm.ScalarMappable(norm=cNorm, cmap=cm.terrain)\n",
    "\n",
    "current_bottom = np.zeros(l)\n",
    "\n",
    "for i in r.__reversed__():\n",
    "\tplt.bar(r, gamma[r, i], 1, color=colorMap.to_rgba(r), bottom=current_bottom)\n",
    "\tcurrent_bottom = current_bottom + gamma[r, i]\n",
    "\n",
    "plt.axis('off')\n",
    "plt.ylim(0, 0.5)\n",
    "# plt.savefig(\"earth_move_p_r.svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}