Wasserstein_Kantorovich.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    ""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some Visualizations for the Earth Mover's Distance\n",
    "see [Wasserstein GAN and the Kantorovich-Rubinstein duality]()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P_r:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAFkCAYAAACq4KjhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAB69JREFUeJzt3TFuVFcUgOHzLAovgDotuGU9bCW7IVvIMtIiqkg0NBQp\nqfxSDEUEP0GDPH5z7e9r3I1O+Wt85p5t3/cBAPjWzdEDAADXSSQAAEkkAABJJAAASSQAAEkkAABJ\nJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkk\nAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQA\nAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAA\nSSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJ\nJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkk\nAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQA\nAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAA\nSSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJ\nJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkk\nAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQA\nAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAA\nSSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAADpxdEDwHO2bdtvM/Py6Dl+wed93z8ePQRw\nWSIBDnIKhJsPM/e3R89yvpsv27a9EgrwtIkEOM7LUyC8m5m7o2c5w/uZeXs7p29ARAI8YSIBDnc3\nM2+OHgLgOxYXAYAkEgCAJBIAgCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQS\nAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIA\ngCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQSAID04ugBgGW93rbt6Bl+xed9\n3z8ePQSsQCQAZ/o0py8h7/84epJfc/Nl27ZXQgF+TiQAZ/pnZu5n5t3M3B08y7nez8zb25l5OTMi\nAX5CJAC/6G5m3hw9BHBBFhcBgCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQSAIAkEgCAJBIAgCQS\nAIAkEgCAdLUHnrZt+21Ol9pW41Y9AE/C1UbCzM2Hmfvbo6c4n1v1ADwNVxwJ97fr3at3qx6Ap+OK\nI2HGvXoAOI7FRQAgiQQAIIkEACCJBAAgiQQAIIkEACCJBAAgiQQAIIkEACCJBAAgiQQAIIkEACCJ\nBAAgiQQAIIkEACCJBAAgiQQAIIkEACCJBAAgiQQAIL04egCAA7zetu3oGc71ed/3j0cPwfMiEoBn\n5NOcvkC9/+PoSc5382XbtldCgcckEoBn5J+ZuZ+ZdzNzd/As53g/M29vZ+blzIgEHo1IAJ6hu5l5\nc/QQcPUsLgIASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAABJJAAASSQAAEkkAADJgafL\nWPFW/Yx79QD8h0h4UCvfqp9xrx6A/xIJD2rVW/Uz7tUD8C2RcBFu1QOwPouLAEASCQBAEgkAQBIJ\nAECyuAiwDm+w8KhEAsDV8wYLxxAJAFfPGywcQyQALMMbLDwui4sAQBIJAEASCQBAEgkAQBIJAEAS\nCQBAEgkAQBIJAEASCQBAEgkAQBIJAEASCQBAcuCJb614r96teoALEAl8tfK9erfqAS5BJPDVqvfq\n3aoHuBSRwDfcqwfgxOIiAJBEAgCQRAIAkEQCAJAsLgLwGLzBsiCRAMAFeYNlZSIBgAvyBsvKRAIA\nj8AbLCuyuAgAJJEAACSRAAAkkQAAJIuLAPBjK77vMPu+//UQnyMSAOA7K7/vMDMzD1I2IgEAvrPq\n+w4zpzceHoZIAIAfet7vO1hcBACSSAAAkkgAAJJIAACSxUWeihV/y/z66AEA/o9IYHHL/5YZ4GqJ\nBBa38m+Z/5yZ348eAuCHRAJPxIq/ZX64B08ALsHiIgCQRAIAkEQCAJBEAgCQrnxxcbXFrr+//l1t\n7pl1Z1917pl1Z1917pl1Z1917pl1Z1917pnTzA+zyL3t+/4gHwQAPC3+3QAAJJEAACSRAAAkkQAA\nJJEAACSRAAAkkQAAJJEAACSRAAAkkQAAJJEAACSRAAAkkQAAJJEAACSRAAAkkQAAJJEAACSRAAAk\nkQAAJJEAACSRAAAkkQAAJJEAACSRAAAkkQAAJJEAACSRAAAkkQAAJJEAAKR/AWYtuEsiApcGAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1112ae2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P_t:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAFkCAYAAACq4KjhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAACCpJREFUeJzt3bFOnEcYheEzkYu9gK0oaOF6uEe33ILvIL3ljnqlNC7o\nJgVJpDgnipYsnh14nsYdOuUrPP/HmHMGAOBHv6weAABcJ5EAAFQiAQCoRAIAUIkEAKASCQBAJRIA\ngEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkE\nAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQi\nAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACV\nSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBA\nJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIA\nUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEA\nAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEok\nAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKAS\nCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCo\nRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAA\nKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIA\ngEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkE\nAKASCQBAJRIAgEokAADVp9UD4P8aY9wmOa7e8UqnOefT6hEAjUhga2OM24x8zcxh9ZZXGXkeY9wJ\nBeAaiQR2d8zMIQ/Z73cJpySPOeRluUgAro5I4H04JrlZPQLgffFwEQCoRAIAUIkEAKASCQBAJRIA\ngEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkE\nAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAqk+rBwC5\nH2Os3nCu05zzafUI4G2JBFjle5KRZObz6ilnG3keY9wJBXjfRAKs8pxkJnlIcly85RynJI855GW1\nSIB3TCTAasckN6tHAPyTh4sAQCUSAIBKJAAAlUgAACqRAABUIgEAqEQCAFCJBACgEgkAQCUSAIBK\nJAAAlUgAACqRAABU/gokfxlj3GavP1qcJPerBwC8VyKBJH8EwsjXzBxWbwHgOogE/nTMzCEP2et3\nCd+SfFk9AuB9Egn83THJzeoRZzitHgDwfnm4CABUIgEAqEQCAFCJBACgEgkAQCUSAIBKJAAAlUgA\nACqRAABUIgEAqEQCAFCJBACgEgkAQCUSAIDKn4oGXut+jLF6w2uc5pxPq0fADkQCcJ7vSUaSmc+r\np7zKyPMY404owH8TCcB5npPMJA9Jjou3nOuU5DGHvCwXCfAfRALwOsckN6tHAG/Jw0UAoBIJAEAl\nEgCASiQAAJVIAAAqkQAAVCIBAKhEAgBQiQQAoBIJAEAlEgCASiQAAJVIAAAqkQAAVCIBAKhEAgBQ\niQQAoBIJAEAlEgCASiQAAJVIAAAqkQAAVJ9WDwBY4H6MsXrDuU5zzqfVI/hYRALwcXxPMpLMfF49\n5Wwjz2OMO6HAzyQSgI/jOclM8pDkuHjLOU5JHnPIy2qRwE8jEoCP55jkZvUIuH4eLgIAlUgAACqR\nAABUIgEAqEQCAFD5uuHCxhi32evjqj/drx4AwHURCRc0xrjNyNfMHFZvAYD/SyRc1jEzh+0OtSTJ\ntyRfVo8A4JqIhLew46GW0+oBAFwbDxcBgEokAACVSAAAKpEAAFQiAQCofN0AsI/7McbqDa9xmnM+\nrR7B+UQCwLX7nmQkmfm8esqrjDyPMe6Ewn5EAsC1e04yky0PtZ2SPOaQl+UiYTMiAWAXOx5qY2se\nLgIAlUgAACqRAABUIgEAqDxcBOBn2PHGw4e/7yASAHg7O994cN9BJADwhna98eC+QxKRAMDP4MbD\nljxcBAAqkQAAVCIBAKhEAgBQiQQAoPJ1AwD8ux2PQGXO+eslfo5IAIAf7XwE6sVFykYkAMCPdj0C\nlbwcgroQkQAA/+aDH4HycBEAqEQCAFCJBACgEgkAQCUSAIDqar9uGGPcZr8PT+5XDwCAS7naSMjI\n18wcVs8AgI/qeiNh5rDdEYtvSb6sHgEAl3G9kZDsd8TigleuAGA1DxcBgEokAACVSAAAKpEAAFTX\n/XBxt4eAv/3x7267k32377o72Xf7rruTfbfvujvZd/uuu5OLbh5zzsv9NADg3fDfDQBAJRIAgEok\nAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKAS\nCQBAJRIAgEokAACVSAAAKpEAAFQiAQCoRAIAUIkEAKASCQBAJRIAgEokAACVSAAAKpEAAFS/AxeS\nofs3TBWcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1116073c8>"
      ]
     },
     "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": 69,
   "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": [
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EMD:  0.706237424547 \n\nGamma:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVgAAAFdCAYAAABGoXXzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAABbRJREFUeJzt2yGOlmcUhuH/NLML0mRW0Q2wgAaBxCMaXF01npSEBdQ2\nYR2YSnRFWQVJT0VxEH4z97yEuS79Js+n7hzzze5eALh7P5z+AIDvlcACRAQWICKwABGBBYgILEBE\nYAEiAgsQuakHfpw59ifDP6eGL5fLnwe3/zq4/fLgNtyn3Z1rb1ywABGBBYgILEBEYAEiAgsQEViA\niMACRAQWICKwABGBBYgILEBEYAEiAgsQEViAiMACRAQWICKwABGBBYgILEBEYAEiAgsQEViAiMAC\nRAQWICKwABGBBYgILEBEYAEiAgsQmd1tB2baga/4/dTw5XJ5c3D7/cFteCh2d669ccECRAQWICKw\nABGBBYgILEBEYAEiAgsQEViAiMACRAQWICKwABGBBYgILEBEYAEiAgsQEViAiMACRAQWICKwABGB\nBYgILEBEYAEiAgsQEViAiMACRAQWICKwABGBBYgILEBEYAEis7vtwEw78I365eD264Pb8FDs7lx7\n44IFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEF\niAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEZnfbgZl2gM88Orj9\n4eA23KfdnWtvXLAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICI\nwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWIDIzekP\n4O49P7j92+6x7Zk5tg1f4oIFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxAR\nWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICI\nwAJEBBYgcnP6A7h7T2/Pbc/Mse0/ji3/79nhfb49LliAiMACRAQWICKwABGBBYgILEBEYAEiAgsQ\nEViAiMACRAQWICKwABGBBYgILEBEYAEiAgsQEViAiMACRAQWICKwABGBBYgILEBEYAEiAgsQEViA\niMACRAQWICKwABGBBYgILEBkdrcdmGkH4JPHh/dfHNz++eD2Q7W7c+2NCxYgIrAAEYEFiAgsQERg\nASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASIC\nCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxCZ3W0HZtoB+OT28P7fB7d/Pbj96uD2vwe3P+7O\ntTcuWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxAR\nWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQGR2tx2YaQeAo54c\n3H57cHt359obFyxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAi\nsAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAyu9sO\nzLQDwIP108Htd7tz7Y0LFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGAB\nIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgIL\nEJndPf0NAN8lFyxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAi\nsAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkT+A3i+Tfoq\nG0pLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110c93860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVgAAAFdCAYAAABGoXXzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAABb5JREFUeJzt271t22AUhlEq4BqBeya9YXgox6N4OfYCB1EGyI9c5OEV\ngnNqAu914Qdfo8vtdlsA+Pe+TB8A8L8SWICIwAJEBBYgIrAAEYEFiAgsQERgASJrPXC5XMZ+ybBN\nDQ9vfxvcnvy7v38dHF+WZXseHH99mtt+fpvbfnmf216Wy70PvGABIgILEBFYgIjAAkQEFiAisAAR\ngQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWI\nCCxARGABIgILEBFYgIjAAkQEFiCy1gNbPfAX++A2A47pA+Zsy3Vw/WNwe9DL+91PvGABIgILEBFY\ngIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjA\nAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisACRtR7Y6oEHtU8fwPmO6QNmbMt1cP1j\nbvrl/e4nXrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJE\nBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWIDIWg98qwf4\nxT59AOc7pg+YsS3XufEf9z/xggWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgIL\nEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFY\ngIjAAkTWemCrB3go+/QBnO+YPmDGZ9rmBQsQEViAiMACRAQWICKwABGBBYgILEBEYAEiAgsQEViA\niMACRAQWICKwABGBBYgILEBEYAEiAgsQEViAiMACRAQWICKwABGBBYgILEBEYAEiAgsQEViAiMAC\nRAQWICKwABGBBYgILEBkrQe+f60X/uIY3OZ0+/QBnG/wf3z7xDdesAARgQWICCxARGABIgILEBFY\ngIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjA\nAkQEFiAisAARgQWICCxARGABIgILEBFYgMhaD2zP9cKDOqYP4Gz79AE8HC9YgIjAAkQEFiAisAAR\ngQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWI\nCCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxAZM0XXp/yiT/ZluvY9qhj+gDOtk8fwG95wQJE\nBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYg\nIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASJrvvD8lk/82cfY8rZcx7ZH\nHdMHcLZ9+oAH5gULEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWI\nCCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWIrPnC\ny3s+8Zg+xpa35Tq2PeqYPoCz7dMH3OEFCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICI\nwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJEBBYgIrAAEYEFiAgsQERgASICCxARWICIwAJE\nBBYgIrAAEYEFiFxut9v0DQD/JS9YgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAi\nsAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAARgQWICCxARGABIgILEBFYgIjAAkQEFiAisAAR\ngQWI/ASfXirGU5OFuAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110cb0a90>"
      ]
     },
     "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": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dual EMD:  0.706237424547\nf: \n [ 1.  0. -1. -2. -1.  0.  1.  2.  1.  0.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAFkCAYAAABxWwLDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X2UZHV54PHvw8v4gsuYYGQEWY0BlSTOYPcIEmFE8SXq\nKoN6NK0GBWI0xEAmm6OMkVUaAxGjjNGgHllFY+yEzTkocpRRNMOsr2S6YYgGJBshGsGBYHbYHYXF\n4dk/qhpqiunuquq6dV/q+zmnz5muvrf6V2eqqp+qb92qyEwkSZL6sU/ZC5AkSfXjACFJkvrmACFJ\nkvrmACFJkvrmACFJkvrmACFJkvrmACFJkvrmACFJkvrmACFJkvrmACFJkvpW6AARERsj4tqIuDsi\ndkTE5RHx5B72OyEiZiPinoi4OSJeX+Q6JUlSf4p+BuJ44IPAMcDzgP2BL0XEIxbaISKeCFwJfAVY\nA3wAuCQinl/wWiVJUo9ilB+mFRGPAe4A1mXm1xbY5j3AizJzdcdpM8DKzHzxaFYqSZIWM+rXQDwa\nSOAni2zzTODqrtM2A8cWtShJktSf/Ub1iyIigE3A1zLznxbZdBWwo+u0HcCBEfGwzLy363wPAl4I\n3ArcM7wVS5LUeA8Hnghszsy7+tlxZAMEcDHwq8Czhny+LwT+esjnKUnSOHkt8Jl+dhjJABERHwJe\nDByfmbcvsfmPgYO7TjsYuLv72Ye2WwE+/elPc+SRRy53qZWwYcMGLrroorKXMTRNujxNuizg5amy\nplyWCy+Ez34W7r13A6effhFnnFH2ioajKf8/N954I6973eug/be0H4UPEO3h4STg2Zn5gx52+Sbw\noq7TXtA+fW/uATjyyCOZmJgYeJ1VsnLlysZcFmjW5WnSZQEvT5U14bJccw387d/Cpk3w4Q+v5NJL\nJzjjDKj5xQKa8f/Tpe+XABT9PhAX03pa5DXArog4uP318I5tzo+IT3bs9hHgSRHxnoh4SkScAbwS\neH+Ra5UkDc+uXXDaaXDccfAHfwCHHw5Pexq84Q3w//5f2avTMBR9FMabgQOBLcBtHV+v6tjmccBh\n899k5q3AS2i9b8T1wAbg9MzsPjJDklRRGzfC7bfDJz4B++zT+rr0UrjxRjjvvLJXp2EoNGFk5pID\nSmaeupfTtgKThSxKklSoa66BD36wlS4OP/zB09esgXPOgelpOPnkZqSMceZnYVTQ1NRU2UsYqiZd\nniZdFvDyVFldL0t3upg3f3k2bmxGyqjr/88wjfSdKIsQERPA7OzsbNNe0CJJtXPmmXDJJXDDDXs+\n+9Bp+3ZYuxbOPtucUba5uTkmJycBJjNzrp99fQZCkjQU8+niggsWHh7gwZRxwQUw19efLFWJA4Qk\nadkWShcLaUrKGGcOEJKkZes+6mIp++/vURl15wAhSVqWXtNFN1NGvTlASJIG1m+66GbKqC8HCEnS\nwPpNF91MGfXlACFJGsig6aKbKaOeHCAkSX1bbrroZsqoHwcISVLflpsuupky6scBQpLUl2Gli26m\njHpxgJAk9WzY6aKbKaM+HCAkST0bdrroZsqoDwcISVJPikoX3UwZ9eAAIUlaUtHpopspo/ocICRJ\nSyo6XXQzZVSfA4QkaVGjShfdTBnV5gAhSVrQqNNFN1NGdTlASJIWNOp00c2UUV0OEJKkvSorXXQz\nZVSTA4Qk6SHKThfdTBnV4wAhSXqIstNFN1NG9VTgaiFJqpKqpItupoxqcYCQJD2gaumimymjOhwg\nJEkPqFq66GbKqI4KXj0kSWWoarroZsqoBgcISVLl00U3U0b5Ch0gIuL4iLgiIn4UEfdHxMuW2P7Z\n7e06v3ZHxGOLXKckjbuqp4tupozyFX01OQC4HjgDyB73SeAIYFX763GZeUcxy5Mk1SVddDNllGu/\nIs88M68CrgKIiOhj1zsz8+5iViVJmle3dNFt40a4/PJWyti2DVasKHtF46OKT1QFcH1E3BYRX4qI\n3yh7QZLUVHVLF91MGeWp2tXlduBNwCuAlwM/BLZExFGlrkqSGqiu6aKbKaMckdnrSxOW+Ysi7gfW\nZ+YVfe63BfjXzHz9Aj+fAGbXrVvHypUr9/jZ1NQUU1NTA65Ykppr1y5YvRoOOaQ1SNTx2YdO990H\nRx8Nu3ebMhYyMzPDzMzMHqft3LmTrVu3AkxmZl/jVx0GiAuBZ2Xmsxb4+QQwOzs7y8TExBBWKknN\nd+aZcMklcMMN9X72odP27bB2bSvLTE+XvZp6mJubY3JyEgYYIOowcx5FK21IkoagKemi23zKOP98\nU8YoFP0+EAdExJqO1zA8qf39Ye2fXxARn+zY/qyIeFlE/EpE/FpEbAKeA3yoyHVK0rio+1EXS/EN\npkan6Gcg1gLXAbO03t/hfcAccG7756uAwzq2X9He5gZgC/A04MTM3FLwOiVpLNT9qIuldB6V8e53\nl72aZiv6fSCuYZEhJTNP7fr+vcB7i1yTJI2r+XSxaVOz0kW3+ZQxPQ3r14MvjytGA+dPSVK3pqeL\nbqaM4jlASNIYaHq66GbKKN4YXI0kabw19aiLpXhURrEcICSpwcYtXXQzZRTHAUKSGmzc0kU3U0Zx\nxvDqJEnjYVzTRTdTRjEcICSpgcY9XXQzZQyfA4QkNdC4p4tupozh82olSQ1jutg7U8ZwOUBIUoOY\nLhZnyhgeBwhJahDTxeJMGcPj1UuSGsJ00RtTxnA4QEhSA5gu+mPKWD4HCElqANNFf0wZy+fVTJJq\nznQxGFPG8jhASFKNmS6Wx5QxOAcISaox08XymDIG59VNkmrKdDEcpozBOEBIUg2ZLobLlNE/BwhJ\nqiHTxXCZMvrn1U6SasZ0UQxTRn8cICSpRkwXxTJl9M4BQpJqxHRRLFNG77z6SVJNmC5Gw5TRGwcI\nSaoB08VomTKW5gAhSTVguhgtU8bSvBpKUsWZLsphylhcoQNERBwfEVdExI8i4v6IeFkP+5wQEbMR\ncU9E3BwRry9yjZJUZaaLcpkyFlb0MxAHANcDZwC51MYR8UTgSuArwBrgA8AlEfH84pYoSdVluiiX\nKWNh+xV55pl5FXAVQERED7v8HvD9zHxr+/vvRcRxwAbgy8WsUpKqaT5dbNpkuijTfMqYnob162Fi\nouwVVUPV5tlnAld3nbYZOLaEtUhSaUwX1WLKeKiqDRCrgB1dp+0ADoyIh5WwHkkqhemiWkwZD9WY\nq+V995W9AkkaDo+6qCaPythToa+BGMCPgYO7TjsYuDsz711sx1e9agNr1qzc47SpqSmmpqaGu0JJ\nKpDpoto2boTLL2+ljG3bYMWKslfUu5mZGWZmZvY4befOnQOfX2QueXDEUETE/cD6zLxikW3+DHhR\nZq7pOO0zwKMz88UL7DMBzO6zzyzXXjvB5OSwVy5Jo3PmmXDJJXDDDT77UFXbt8Pata1hYnq67NUs\nz9zcHJOtP5yTmdnX8ypFvw/EARGxJiKOap/0pPb3h7V/fkFEfLJjl4+0t3lPRDwlIs4AXgm8f6nf\ndfjhvrhFUr2ZLurBlNFS9Gsg1gLXAbO03gfifcAccG7756uAw+Y3zsxbgZcAz6P1/hEbgNMzs/vI\njIc491y46SY477xhLl+SRmM+XRx/vOmiDjwqo+ABIjOvycx9MnPfrq/T2j8/NTOf27XP1syczMxH\nZOYRmflXvfyuJz8Z3vGO1uQ+O1vEpZGk4swfdfHxj3vURR14VEaDjsIAePvbnQgl1Y/pop7GPWU0\naoCYnwhNGZLqwnRRb+OcMho1QEBrIjRlSKoL00W9jXPKaOTV1ZQhqQ5MF80wrimjkQOEKUNS1Zku\nmmUcU0YjBwgwZUiqNtNFs4xjymj01daUIamKTBfNNG4po9EDhClDUtWYLpptnFJGowcIMGVIqhbT\nRbONU8oYi6uvKUNSFZguxsO4pIyxGCBMGZLKZroYL+OQMsZigABThqRymS7GyzikjLG6GpsyJJXB\ndDGemp4yxmqAMGVIGjXTxXhrcsoYqwECTBmSRst0Md6anDLG8upsypA0CqYLQXNTxlgOEKYMSUUz\nXahTE1PGWA4QYMqQVCzThTo1MWWM9dXalCGpCKYL7U3TUsZYDxCmDEnDZrrQYpqUMsZ6gABThqTh\nMl1oMU1KGV69MWVIGg7ThXrRlJThAIEpQ9LymS7UjyakDAeINlOGpOUwXagfTUgZXs07mDIkDcJ0\noUHUPWU4QHQwZUjql+lCy1HnlOEA0cWUIakfpgstR51Thlf3vTBlSOqF6ULDUNeUUfgAERG/HxG3\nRMTPIuJbEfGMRbZ9dkTc3/W1OyIeW/Q6O5kyJC3FdKFhqmPKKHSAiIhXA+8D3gk8HdgObI6Ixyyy\nWwJHAKvaX4/LzDuKXOfemDIkLcZ0oWGqY8oo+mq/AfhoZn4qM28C3gz8FDhtif3uzMw75r8KXuOC\nTBmS9sZ0oSLULWUUNkBExP7AJPCV+dMyM4GrgWMX2xW4PiJui4gvRcRvFLXGpZgyJHUzXahIdUoZ\nRT4D8RhgX2BH1+k7aKWJvbkdeBPwCuDlwA+BLRFxVFGLXIopQ1In04WKVKeUUamrf2benJkfy8zr\nMvNbmXk68A1aKaQ0pgxJYLrQaNQlZUSrKhRwxq2E8VPgFZl5RcfplwIrM/PkHs/nQuBZmfmsBX4+\nAcyuW7eOlStX7vGzqakppqamBrwEe9q+HdauhbPPNmdI42jXLli9Gg49FLZs8dkHFeu+++Doo2H3\nbti2DVasWP55zszMMDMzs8dpO3fuZOvWrQCTmdnXuFLYAAEQEd8Cvp2ZZ7W/D+AHwF9k5nt7PI8v\nAXdn5isX+PkEMDs7O8vExMSQVr53557bGh6+/W2YnCz0V0mqmDPPhEsugRtu8NkHjcb8A9eNG2F6\nupjfMTc3x2TrD1rfA0TRM/T7gTdGxCkR8VTgI8AjgUsBIuKCiPjk/MYRcVZEvCwifiUifi0iNgHP\nAT5U8Dp7YsqQxpPpQmWYfw1eVVNGoQNEZl4G/DEwDVwHrAZemJl3tjdZBRzWscsKWu8bcQOwBXga\ncGJmbilynb3yqAxp/HjUhcpU5QeuhVe8zLw4M5+YmY/IzGMzc1vHz07NzOd2fP/ezDwiMw/IzF/K\nzBMzc2vRa+yHR2VI48WjLlSmKh+V4c1hAFWeCCUNj+lCVVDVlOEAMQBThtR8pgtVSRUfuDpADMiU\nITWb6UJVUsWU4c1iGao4EUpaPtOFqqhqKcMBYhlMGVLzmC5UZVV64OoAsUymDKlZTBeqsiqlDG8e\nQ1CliVDS4EwXqoOqpAwHiCEwZUj1Z7pQnVThgasDxJCYMqR6M12oTqqQMryZDFEVJkJJ/TNdqI7K\nThkOEENkypDqx3ShOivzgasDxJCZMqR6MV2ozspMGd5cCmDKkOrBdKEmKCtlOEAUwJQhVZ/pQk1S\nxgNXB4iCmDKkajNdqEnKSBnebApkypCqyXShJhp1ynCAKJApQ6oe04WabJQPXB0gCmbKkKrFdKEm\nG2XK8OYzAqYMqRpMFxoHo0oZDhAjYMqQyme60DgZxQNXB4gRMWVI5TJdaJyMImV4MxohU4ZUDtOF\nxlHRKcMBYoRMGdLomS40zop84OoAMWKmDGm0TBcaZ0WmDG9OJTBlSKNhupCKSxkOECUwZUjFM11I\nDyrigasDRElMGVKxTBfSg4pIGd6sSmTKkIphupAeatgpo/ABIiJ+PyJuiYifRcS3IuIZS2x/QkTM\nRsQ9EXFzRLy+6DWWxZQhDZ/pQlrYMB+4FjpARMSrgfcB7wSeDmwHNkfEYxbY/onAlcBXgDXAB4BL\nIuL5Ra6zTKYMabhMF9LChpkyir55bQA+mpmfysybgDcDPwVOW2D73wO+n5lvzczvZeZfAn/XPp/G\nMmVIw2G6kJbWmTJuvHHw8ylsgIiI/YFJWs8mAJCZCVwNHLvAbs9s/7zT5kW2bwRThrR8pgupd/MP\nXN/1rsHPo8hnIB4D7Avs6Dp9B7BqgX1WLbD9gRHxsOEur1pMGdLymC6k3s0/cL3llsHPY7+hraZk\nGzZsYOXKlXucNjU1xdTUVEkr6t/b3w6f/WwrZczOwooVZa9IqoetW1vpYtMm04W0kJmZGWZmZvY4\n7dBDd/KDHwx2ftGqCsPXThg/BV6RmVd0nH4psDIzT97LPtcAs5n5Rx2nvQG4KDN/YYHfMwHMzs7O\nMjExMdwLUYLt22HtWjj7bHOG1Itdu2D1ajj0UNiyxWcfpH5ce+0cxxwzCTCZmX0d3FnYTS0z7wNm\ngRPnT4uIaH//jQV2+2bn9m0vaJ8+FtasgXPOaaWMIj49TWoa04U0uP2W0SGKvrm9H3hjRJwSEU8F\nPgI8ErgUICIuiIhPdmz/EeBJEfGeiHhKRJwBvLJ9PmNj40aPypB64VEXUnkKHSAy8zLgj4Fp4Dpg\nNfDCzLyzvckq4LCO7W8FXgI8D7ie1uGbp2dm95EZjdZ5nK4ZQ9o7j7qQylX4iygz82Lg4gV+dupe\nTttK6/DPsTafMqan4eSToQEv75CGaj5dbN5supDK4M2uwkwZ0t6ZLqTyOUBUmClDeijThVQNDhAV\n51EZ0p486kKqBm9+NWDKkFpMF1J1OEDUgClDMl1IVeMAUROmDI0704VULd4Ma8SUoXFlupCqxwGi\nRkwZGkemC6maHCBqxpShcWO6kKrJm2MNmTI0LkwXUnU5QNSQKUPjwHQhVZsDRE2ZMtR0pgup2rxZ\n1pgpQ01lupCqzwGixkwZaiLThVQPDhA1Z8pQ05gupHrw5tkApgw1helCqg8HiAYwZagJTBdSvThA\nNIQpQ3VnupDqxZtpg5gyVFemC6l+HCAaxJShOjJdSPXkANEwpgzVjelCqidvrg1kylBdmC6k+nKA\naCBThurAdCHVmwNEQ5kyVHWmC6nevNk2mClDVWW6kOrPAaLBTBmqItOF1AwOEA1nylDVmC6kZvDm\nOwZMGaoK04XUHIUNEBHxCxHx1xGxMyL+IyIuiYgDltjnExFxf9fXF4pa47gwZagKTBdSsxT5DMRn\ngCOBE4GXAOuAj/aw3xeBg4FV7a+pohY4TkwZKpvpQmqWQm7GEfFU4IXA6Zm5LTO/AfwB8FsRsWqJ\n3e/NzDsz8472184i1jiOTBkqi+lCap6iHgccC/xHZl7XcdrVQALHLLHvCRGxIyJuioiLI+IXC1rj\n2DFlqAymC6mZihogVgF3dJ6QmbuBn7R/tpAvAqcAzwXeCjwb+EJEREHrHDumDI2a6UJqpv362Tgi\nLgDetsgmSet1DwPJzMs6vv1uRPwj8C/ACcDfL7bvhg0bWLly5R6nTU1NMTXlSyi6bdwIl1/eShnb\ntsGKFWWvSE01ny42bTJdSGWbmZlhZmZmj9N27hz8VQKRmb1vHHEQcNASm30f+G3gzzPzgW0jYl/g\nHuCVmfm5Pn7nHcCfZObHFvj5BDA7OzvLxMREr2c79rZvh7VrW8PE9HTZq1ET7doFq1fDIYe0Bgmf\nfZCqZ25ujsnJSYDJzOzreem+noHIzLuAu5baLiK+CTw6Ip7e8TqIE4EAvt3r74uIx9MaWG7vZ51a\n2nzKmJ6G9evB2UvDNp8uNm92eJCaqJCbdWbeBGwGPhYRz4iIZwEfBGYy88fz27VfKHlS+98HRMSF\nEXFMRDwhIk4EPgvc3D4vDZlHZagoHnUhNV+RjwteA9xE6+iLK4GtwJu6tjkCmH/hwm5gNfA54HvA\nx4B/ANZl5n0FrnNsdR6V8e53l70aNcX8URfHHedRF1KT9ZUw+pGZ/xt43RLb7Nvx73uA3yxqPdo7\nU4aGzXQhjQdv3jJlaGhMF9L4cICQKUNDYbqQxosDhIAHU8b55/sGUxrMfLr4xCdMF9I48GauB5gy\nNCjThTR+HCD0AFOGBmG6kMaTA4T2YMpQv0wX0njy5q6HMGWoV6YLaXw5QOghTBnqhelCGm8OENor\nU4aWYrqQxps3ey3IlKGFmC4kOUBoQaYM7Y3pQhI4QGgJpgx1M11IAgcI9cCUoXmmC0nzHCC0JFOG\nwHQhaU8OEOqJKUOmC0mdvBtQz0wZ48t0IambA4R6ZsoYT6YLSXvjAKG+mDLGj+lC0t54d6C+mTLG\nh+lC0kIcINQ3U8Z4MF1IWowDhAZiymg+04WkxXi3oIGZMprLdCFpKQ4QGpgpo5lMF5J64QChZTFl\nNI/pQlIvvHvQspkymsN0IalXDhBaNlNGM5guJPXDAUJDYcqoP9OFpH54N6GhMWXUl+lCUr8KGyAi\n4u0R8fWI2BURP+ljv+mIuC0ifhoRX44I785qwpRRT6YLSYMo8hmI/YHLgA/3ukNEvA14C/C7wNHA\nLmBzRKwoZIUaOlNG/ZguJA2isLuLzDw3Mz8A/GMfu50FnJeZV2bmd4BTgEOA9UWsUcUwZdSH6ULS\noCrzeCMifhlYBXxl/rTMvBv4NnBsWetS/0wZ9WC6kLQclRkgaA0PCezoOn1H+2eqEVNG9ZkuJC3H\nfv1sHBEXAG9bZJMEjszMm5e1qgFs2LCBlStX7nHa1NQUU1NTo16K2jZuhMsvb6WMbdtgha9kqYz5\ndHHRRaYLaVzMzMwwMzOzx2k7d+4c+PwiM3vfOOIg4KAlNvt+Zv68Y5/XAxdl5i8ucd6/DPwLcFRm\n3tBx+hbguszcsMB+E8Ds7OwsExMTvV0Qjcz27bB2bWuYmJ4uezWCVrpYvRoOOaQ1SPjsgzS+5ubm\nmJycBJjMzL6eL+7rGYjMvAu4q599+jjvWyLix8CJwA0AEXEgcAzwl0X8ThVvPmVMT8P69eCMV775\ndHHVVQ4PkgZX5PtAHBYRa4AnAPtGxJr21wEd29wUESd17LYJeEdEvDQingZ8Cvg34HNFrVPF86iM\n6phPF+efD0ccUfZqJNVZkY8/poE54J3Ao9r/ngMmO7Y5AnjghQuZeSHwQeCjtI6+eATwosz0z06N\neVRGNXQedXHmmWWvRlLd9ZUw+pGZpwKnLrHNvns57V3Au4pZlcpiyiif6ULSMHk3opExZZTHdCFp\n2BwgNDKmjHKYLiQVwQFCI+UbTI3efLr4+MdNF5KGx7sTjZwpY3RMF5KK4gChkTNljIbpQlKRHCBU\nClNG8UwXkork3YpKY8oojulCUtEcIFQaU0YxTBeSRsEBQqUyZQyf6ULSKHj3otKZMobHdCFpVBwg\nVDpTxnCYLiSNkgOEKsGUsXymC0mj5N2MKsOUMTjThaRRc4BQZZgyBmO6kFQGBwhViimjf6YLSWXw\n7kaVY8ronelCUlkcIFQ5pozemC4klckBQpVkylia6UJSmbzbUWWZMhZmupBUNgcIVZYpY+9MF5Kq\nwAFClWbKeCjThaQq8O5HlWfKeJDpQlJVOECo8kwZLaYLSVXiAKFaMGWYLiRVi3dDqo1xThmmC0lV\n4wCh2hjXlGG6kFRFDhCqlXFMGaYLSVVU2N1RRLw9Ir4eEbsi4ic97vOJiLi/6+sLRa1R9TROKcN0\nIamqinw8sz9wGfDhPvf7InAwsKr9NTXkdanmxiVlmC4kVdl+RZ1xZp4LEBGv73PXezPzzgKWpAaZ\nTxnT07B+PUxMlL2i4ZtPF1ddZbqQVD1VvFs6ISJ2RMRNEXFxRPxi2QtSNTU5ZZguJFVd1QaILwKn\nAM8F3go8G/hCRESpq1IlNTVlmC4k1UFfA0REXLCXFzl2fu2OiCcPupjMvCwzr8zM72bmFcB/AY4G\nThj0PNVsTTwqw6MuJNVBZGbvG0ccBBy0xGbfz8yfd+zzeuCizBwoRUTEHcCfZObHFvj5BDC7bt06\nVq5cucfPpqammJryNZhNd999cPTRsHs3bNsGK1aUvaLBXXMNnHACXHQR/OEflr0aSU0yMzPDzMzM\nHqft3LmTrVu3AkxmZl8Pw/oaIAaxnAEiIh4P/CtwUmZeucA2E8Ds7OwsE018JZ16sn07rF3bevQ+\nPV32agazaxesXg2HHNIaJHz2QVLR5ubmmJychAEGiCLfB+KwiFgDPAHYNyLWtL8O6Njmpog4qf3v\nAyLiwog4JiKeEBEnAp8FbgY2F7VONUMTUobpQlKdFHk3NQ3MAe8EHtX+9xww2bHNEcB8d9gNrAY+\nB3wP+BjwD8C6zLyvwHWqIep8VIZHXUiqmyLfB+JU4NQlttm349/3AL9Z1HrUfPNHZaxd2zoqoy4p\nw6MuJNWRT5SqUeqYMkwXkurIuys1Tp1ShulCUl05QKhx6vIGU6YLSXXmAKFGqkPKMF1IqjPvttRY\nVU4ZpgtJdecAocaqasowXUhqAgcINVoVU4bpQlITePelxqtSyjBdSGoKBwg1XlVShulCUpM4QGgs\nVCFlmC4kNYl3YxobZaYM04WkpnGA0NgoK2WYLiQ1kQOExkoZKcN0IamJvDvT2BllyjBdSGoqBwiN\nnVGlDNOFpCZzgNBYGkXKMF1IajLv1jS2ikwZpgtJTecAobFVVMowXUgaBw4QGmtFpAzThaRx4N2b\nxt4wU4bpQtK4cIDQ2BtWyjBdSBonDhASw0kZpgtJ48S7OaltOSnDdCFp3DhASG2DpgzThaRx5AAh\ndRgkZZguJI0j7+6kLv2kDNOFpHHlACF16TVlmC4kjTMHiAqamZkpewlDVcfLs1DK6LwsTUgXdfy/\nWUyTLk+TLgt4eZqokLu9iHhCRFwSEd+PiJ9GxD9HxLsiYv8e9p2OiNva+305Ig4vYo1V1rQrZl0v\nz95SxvxlaUq6qOv/zUKadHmadFnAy9NERT1ueioQwBuBXwU2AG8G/nSxnSLibcBbgN8FjgZ2AZsj\nYkVB65QWtFDKMF1IUkEDRGZuzszTM/MrmXlrZl4J/Dnw8iV2PQs4LzOvzMzvAKcAhwDri1intJS9\npYwmpAtJWq5R3v09GvjJQj+MiF8GVgFfmT8tM+8Gvg0cW/jqpAV0pow772xGupCk5dpvFL+k/TqG\ntwB/tMhmq4AEdnSdvqP9s4U8HODGG29czhIrZefOncwN66MhK6AJl+dtb4PXvQ7uv38nRx01x3HH\nDe/TO8vUhP+bTk26PE26LODlqaqOv50P73ffyMzeN464AHjbIpskcGRm3tyxz6HAFuCrmfmmRc77\nWOBrwCGUQoMoAAAHL0lEQVSZuaPj9L8F7s/MqQX2ew3w1z1fCEmS1O21mfmZfnbo9xmIPwc+scQ2\n35//R0QcAnwV+Npiw0Pbj2m98PJg9nwW4mDgukX22wy8FrgVuGeJ3yFJkh70cOCJtP6W9qWvZyD6\nOuPWMw9fBf4B+O3s4RdFxG3AezPzovb3B9IaJk7JzP9RyEIlSVLfinofiENoZYt/Bd4KPDYiDo6I\ng7u2uykiTuo4aRPwjoh4aUQ8DfgU8G/A54pYpyRJGkxRL6J8PvCk9tcP26cFrddI7Nux3RHAyvlv\nMvPCiHgk8FFaR238T+BFmdnnhytLkqQiFZYwJElSc/k2OJIkqW8OEJIkqW+1HyAi4vcj4paI+FlE\nfCsinlH2mgYREcdHxBUR8aOIuD8iXlb2mgYVERsj4tqIuDsidkTE5RHx5LLXNaiIeHNEbI+Ine2v\nb0TEb5a9rmGIiLPb17f3l72WQUTEO9vr7/z6p7LXtRwRcUhE/FVE/Hv7QwW3R8RE2esaRPu+ufv/\n5/6I+GDZa+tXROwTEed1fEjk/4qId5S9ruWIiEdFxKaIuLV9mb4WEWt73b/WA0REvBp4H/BO4OnA\ndlofvvWYUhc2mAOA64EzaL3YtM6OBz4IHAM8D9gf+FJEPKLUVQ3uh7TeQG0CmKR1ePLnIuLIUle1\nTO1h+3dp3W7q7Du03i9mVfvruHKXM7iIeDTwdeBe4IXAkcB/Bf6jzHUtw1oe/H9ZResF9glcVuai\nBnQ28CZa99FPpXWE4Vsj4i2lrmp5/jtwIq33Uvp14MvA1RHxuF52rvWLKCPiW8C3M/Os9vdB687+\nLzLzwlIXtwwRcT+wPjOvKHstw9Ae6O4A1mXm18pezzBExF3AH2fmUm+sVkkR8ShgFvg94Bzgusxc\n7K3mKyki3gmclJm1fITeLSL+DDg2M59d9lqKEBGbgBdnZu2ekYyIzwM/zsw3dpz2d8BPM/OU8lY2\nmIh4OPB/gJdm5lUdp28DvpCZ/22p86jtMxARsT+tR4OdH76VwNX44VtV82hajzoW/DC1umg/jflb\nwCOBb5a9nmX4S+DzmfnVshcyBEe009+/RMSnI+Kwshe0DC8FtkXEZe38NxcRv1P2ooahfZ/9WlqP\neuvoG8CJEXEEQESsAZ4FfKHUVQ1uP1pvq3Bv1+k/o8dn8UbyYVoFeQytC7+3D996yuiXo71pPyu0\nidbbmde2TUfEr9MaGOan9pMz86ZyVzWY9gB0FK2nl+vuW8AbgO8BjwPeBWyNiF/PzF0lrmtQT6L1\nrND7gD8Fjgb+IiLuzcy/KnVly3cyrff9+WTZCxnQnwEHAjdFxG5aD8D/JDP/ptxlDSYz/29EfBM4\nJyJuovW38zW0HoD/cy/nUecBQvVwMfCrtCb1OrsJWEPrDvCVwKciYl3dhoiIeDytge55mXlf2etZ\nrszsfP/+70TEtbTeAfdVLP25PVW0D3BtZp7T/n57e3h9M1D3AeI04IuZ+eOyFzKgV9P6A/tbwD/R\nGsI/EBG31Xi4ex3wceBHwM+BOeAztJ7dX1KdB4h/B3bTevFUp4NpfTCXShYRHwJeDByfmbeXvZ7l\nyMyf8+AHxV0XEUcDZ9F6tFgnk8AvAXPtZ4eg9UzeuvaLwR7Wy+fWVFVm7oyIm4HDy17LgG4Hbuw6\n7Ubg5SWsZWgi4j/TekH1+rLXsgwXAhd0fC7TdyPiicBGajrcZeYtwHPaL3A/MDN3RMTf0PGhmIup\n7Wsg2o+eZmm9ghR44OnyE2m1KpWoPTycBDwnM39Q9noKsA/wsLIXMYCrgafRevS0pv21Dfg0sKbO\nwwM88OLQw2n9Ia6jr/PQBPsUWs+q1NlptJ4ir+vrBaD1uqfdXafdT43/js7LzJ+1h4dfoHX0z2d7\n2a/Oz0AAvB+4NCJmgWuBDbT+ky8tc1GDiIgDaN3xzT8qfFL7RTo/ycwfLrxn9UTExcAU8DJgVzz4\nIWo7M7N2H7keEecDXwR+APwnWi8EezbwgjLXNYj26wL2eC1KROwC7srM7ke+lRcR7wU+T+sP7KHA\nucB9wEyZ61qGi4CvR8RGWoc6HgP8DvDGRfeqsPYDuzcAl2bm/SUvZzk+T+vDHv8N+C6tw7o3AJeU\nuqpliIgX0Pqb8z1an011Ia37h0t72b/WA0RmXtY+RHCaVrq4HnhhZt5Z7soGshb4e1pHKyStF1FB\n6wVHp5W1qAG9mdZl2NJ1+qm0PmG1bh5L6//hccBO4AbgBQ05ggHq/b4jj6fVbA8C7gS+BjwzM+8q\ndVUDysxtEXEyrRfsnQPcApxV1xfqtT0POIx6vial01uA82gdwfRY4Dbgw+3T6molcAGt4fsnwN8B\n78jM7mda9qrW7wMhSZLKUft2I0mSRs8BQpIk9c0BQpIk9c0BQpIk9c0BQpIk9c0BQpIk9c0BQpIk\n9c0BQpIk9c0BQpIk9c0BQpIk9c0BQpIk9e3/AzHh6muOTX+yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1125dcc88>"
      ]
     },
     "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": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAFkCAYAAABxWwLDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XuUpHV54PHvIzAacRkTXWdEiUogSjYy2I2jxCBj8BJ1\nFTAmptWgQLzE29ibHA0El8AkTMQLo0ajR1bQGDvLeo43jjgGDcN6QZxuwGhAshGiERgQcobdUchk\nePaPqoKeYvpS1fXWe/t+zulzmHfet+ZX3Vr1dH3ftyoyE0mSpEE8qOwFSJKk+nGAkCRJA3OAkCRJ\nA3OAkCRJA3OAkCRJA3OAkCRJA3OAkCRJA3OAkCRJA3OAkCRJA3OAkCRJAyt0gIiI0yPiqoi4KyJ2\nRMRnIuKXl3HchoiYjYi7I+KGiHhVkeuUJEmDKfoViGOBDwBPA54NHAB8OSJ+bqEDIuLxwCXAV4B1\nwPuACyLiOQWvVZIkLVOM88O0IuKRwG3AMzPzawvs807g+Zl55LxtM8DqzHzBeFYqSZIWM+5zIB4O\nJHDnIvs8Hbisb9tW4JiiFiVJkgaz/7j+oYgIYAvwtcz8x0V2XQvs6Nu2AzgoIh6cmff03e4jgOcB\nNwF3j27FkiQ13kOAxwNbM/OOQQ4c2wABfAj4FeAZI77d5wF/M+LblCSpTV4BfGqQA8YyQETEXwIv\nAI7NzFuW2P1WYE3ftjXAXf2vPnTdBPDJT36SI444YqVLrYTp6WnOP//8spcxMk24P7Oz8NrXwn77\nTfOEJ5zPJz8JBxxQ9qpWrgk/m/madH+adF/A+1NV1113Ha985Suh+1w6iMIHiO7wcAJwXGb+cBmH\nfBN4ft+253a378vdAEcccQQTExNDr7NKVq9e3Zj7AvW/P7t2wW//Nhx7LESs5hvfmOCSS2DTprJX\ntnJ1/9n0a9L9adJ9Ae9PDQx8CkDR7wPxITovi7wc2BURa7pfD5m3z7kR8fF5h30YODQi3hkRT4yI\nNwAvBd5b5FqlhZx+OtxyC3zsY7B6NZx5JmzeDHNzZa9MkspT9FUYrwcOAi4Hbp739Tvz9nk0cEjv\nD5l5E/BCOu8bcQ0wDZyWmf1XZkiF27YNPvCBzsBw2GGdbWecAU9+Mrz61fDv/17q8iSpNIUmjMxc\nckDJzFP2se0KYLKQRUnLtGsXnHpqJ128+c33bz/gALjoIjj66E7GaELKkKRB+VkYFTQ1NVX2Ekaq\nrvdnfrp4UPf/Kb37sm5dM1JGXX82C2nS/WnSfQHvTxON9Z0oixARE8Ds7Oxs005oUYm2bYMNG2DL\nFti4cd/77N4N69fDnj2wfTusWjXWJUrSis3NzTE5OQkwmZkD/TrkKxBSn4XSRb9eyrjuOjOGpPZx\ngJD67CtdLKQpKUOSBuUAIc2zr6suluJVGZLayAFC6lpuuuhnypDURg4QUtcg6aKfKUNS2zhASAyX\nLvqZMiS1iQOEWm/YdNHPlCGpTRwg1HorSRf9TBmS2sIBQq02inTRz5QhqQ0cINRao0oX/UwZktrA\nAUKtNcp00c+UIanpHCDUSkWki36mDElN5gCh1ikqXfQzZUhqMgcItU6R6aKfKUNSUzlAqFXGkS76\nmTIkNZEDhFpjXOminylDUhM5QKg1xpku+pkyJDWNA4RaoYx00c+UIalJHCDUeGWli36mDElN4gCh\nxiszXfQzZUhqCgcINVoV0kU/U4akJnCAUGNVJV30M2VIagIHCDVWldJFP1OGpLqr2MOqNBpVTBf9\nTBmS6swBQo1T1XTRz5Qhqc4cINQ4VU4X/UwZkuqq0IfXiDg2Ij4fET+OiHsj4sVL7H9cd7/5X3si\n4lFFrlPNUYd00c+UIamOiv797EDgGuANQC7zmAQOB9Z2vx6dmbcVszw1SV3SRT9ThqQ62r/IG8/M\nLwFfAoiIGODQ2zPzrmJWpabqpYutW6ufLvr1UsamTXDSSTAxUfaKJGlxVXyYDeCaiLg5Ir4cEb9W\n9oJUfXVMF/1MGZLqpGoDxC3A64DfAl4C/Ai4PCKOKnVVqrS6pot+pgxJdVJowhhUZt4A3DBv05UR\n8UvANPCqxY6dnp5m9erVe22bmppiampq5OtUtdQ5XfQzZUgqyszMDDMzM3tt27lz59C3F5nLPbdx\nZSLiXuDEzPz8gMedBzwjM5+xwN9PALOzs7NM+GjbOtu2wYYNsGULbNxY9mpGY/duWL8e9uyB7dth\n1aqyVySpqebm5picnASYzMyBLiavw+9rR9FJG9JempIu+pkyJNVBoQkjIg4EDqNzYiTAoRGxDrgz\nM38UEZuBgzPzVd39NwI3At8DHgK8BngW8Jwi16l6alK66GfKkFR1RT/sHg1cDczSeX+H9wBzwNnd\nv18LHDJv/1Xdfb4DXA48GTg+My8veJ2qmSZcdbEUr8qQVGVFvw/ENhYZUjLzlL4/vwt4V5FrUv01\nNV3066WMo4/uvBJhzpBUJQ174VdtUKfPulgpPytDUlU1/OFXTdOGdNHPlCGpihwgVBttSRf9vCpD\nUhU5QKg22pQu+pkyJFVNyx6GVVdtTBf9TBmSqsQBQpXX1nTRz5QhqUocIFR5bU4X/UwZkqqi5Q/H\nqjrTxQOZMiRVgQOEKst0sW+mDElV4AChyjJdLMyUIalsPiyrkkwXSzNlSCqTA4Qqx3SxPKYMSWVy\ngFDlmC6Wz5QhqSw+PKtSTBeDM2VIKoMDhCrDdDEcU4akMjhAqDJMF8MzZUgaNx+mVQmmi5UzZUga\nJwcIlc50MRqmDEnj5ACh0pkuRseUIWlcfLhWqUwXo2fKkDQODhAqjemiGKYMSePgAKHSmC6KY8qQ\nVDQftlUK00XxTBmSiuQAobEzXYyHKUNSkRwgNHami/ExZUgqig/fGivTxfiZMiQVwQFCY2O6KIcp\nQ1IRCh0gIuLYiPh8RPw4Iu6NiBcv45gNETEbEXdHxA0R8aoi16jxMV2Ux5QhadSKfhg/ELgGeAOQ\nS+0cEY8HLgG+AqwD3gdcEBHPKW6JGgfTRflMGZJGqdABIjO/lJn/PTM/B8QyDvkD4AeZ+bbM/H5m\nfhD4NDBd5DpVLNNFNZgyJI1S1V5IfjpwWd+2rcAxJaxFI2K6qA5ThqRRqdrD+VpgR9+2HcBBEfHg\nEtajFTJdVI8pQ9IoVG2AUIPs2gWnnWa6qBpThqRR2L/sBfS5FVjTt20NcFdm3rPYgdPT06xevXqv\nbVNTU0xNTY12hVq200+Hm2+GSy81XVRNL2Vs2gQnnQQTE2WvSFLRZmZmmJmZ2Wvbzp07h769yFzy\n4oiRiIh7gRMz8/OL7PMXwPMzc928bZ8CHp6ZL1jgmAlgdnZ2lgkfBStj2zbYsAG2bIGNG8tejfZl\n925Yvx727IHt22HVqrJXJGnc5ubmmJycBJjMzIHOjCr6fSAOjIh1EXFUd9Oh3T8f0v37zRHx8XmH\nfLi7zzsj4okR8QbgpcB7i1ynRsurLurBlCFpJYp+Yflo4Gpgls77QLwHmAPO7v79WuCQ3s6ZeRPw\nQuDZdN4/Yho4LTP7r8xQhXnVRX2sWwfveIdXZUgaXKHnQGTmNhYZUjLzlH1suwKYLHJdKk7vqost\nW7zqoi5OPx0+85nOVRmmDEnL5e+HGhnTRT2ZMiQNwwFCI2O6qC9ThqRB+TCvkfANo+rv9NN9gylJ\ny+cAoRUzXTSDKUPSIBwgtGKmi+YwZUhaLh/utSKmi+YxZUhaDgcIDc100UymDEnL4QChoZkumsuU\nIWkpPuxrKKaL5jNlSFqMA4QGZrpoB1OGpMU4QGhgpov2MGVIWogP/xqI6aJ9TBmS9sUBQstmumgn\nU4akfXGA0LKZLtrLlCGpn08DWhbThUwZkuZzgNCSTBcCU4akvTlAaEmmC/WYMiT1+HSgRZku1M+U\nIQkcILQI04X2xZQhCRwgtAjThRZiypDk04L2yXShpZgypHZzgNADmC60HKYMqd0cIPQApgstlylD\nai+fHrQX04UGZcqQ2skBQvcxXWgYpgypnRwgdB/ThYZlypDax6cJAaYLrZwpQ2oXBwiZLjQSpgyp\nXQofICLijRFxY0T8LCKujIinLrLvcRFxb9/Xnoh4VNHrbDPThUbFlCG1R6FPFxHxMuA9wFnAU4Br\nga0R8chFDkvgcGBt9+vRmXlbketsM9OFRs2UIbVD0b9vTgMfycxPZOb1wOuBnwKnLnHc7Zl5W++r\n4DW2lulCRTBlSO1Q2AAREQcAk8BXetsyM4HLgGMWOxS4JiJujogvR8SvFbXGtjNdqCimDKn5inza\neCSwH7Cjb/sOOmliX24BXgf8FvAS4EfA5RFxVFGLbCvThYpmypCarVK/d2bmDZn50cy8OjOvzMzT\ngG/QSSEaEdOFxmF+yvizPyt7NZJGbf8Cb/snwB5gTd/2NcCtA9zOVcAzltppenqa1atX77VtamqK\nqampAf6pduili61bTRcqVi9lnHMOnHgiTEyUvSKpvWZmZpiZmdlr286dO4e+veicllCMiLgS+FZm\nbuz+OYAfAu/PzHct8za+DNyVmS9d4O8ngNnZ2VkmfHRa0rZtsGEDbNkCGzeWvRq1we7dsH497NkD\n27fDqlVlr0hSz9zcHJOTkwCTmTnQGUtF//75XuA1EXFyRDwJ+DDwUOAigIjYHBEf7+0cERsj4sUR\n8UsR8V8iYgvwLOAvC15nK/TSxa//uulC42PKkJqpyIRBZl7cfc+Hc+iki2uA52Xm7d1d1gKHzDtk\nFZ33jTiYzuWe3wGOz8wrilxnW5guVBZThtQ8hSaMcTBhLI/pQmUzZUjVU+WEoQowXagKTBlSszhA\ntEAvXVx4oelC5eqljHPP9Q2mpLrz6aThfMMoVY1vMCU1gwNEg5kuVEWmDKkZHCAazHShqjJlSPXn\n00pDmS5UdaYMqd4cIBrIdKE6MGVI9eYA0UCmC9WFKUOqL59eGsZ0oboxZUj15ADRIKYL1ZEpQ6on\nB4gGMV2orkwZUv34NNMQpgvVnSlDqhcHiAYwXagJTBlSvThANIDpQk1hypDqw6ebmjNdqGlMGVI9\nOEDUmOlCTWTKkOrBAaLGTBdqKlOGVH0+7dSU6UJNZ8qQqs0BooZMF2oDU4ZUbQ4QNWS6UFuYMqTq\n8umnZkwXahtThlRNDhA1YrpQG5kypGpygKgR04XaypQhVY9PQzVhulDbmTKkanGAqAHThWTKkKrG\nAaIGTBdShylDqg6fjirOdCHtzZQhVYMDRIWZLqQHMmVI1VD4ABERb4yIGyPiZxFxZUQ8dYn9N0TE\nbETcHRE3RMSril5jVZkupH0zZUjlK/RpKSJeBrwHOAt4CnAtsDUiHrnA/o8HLgG+AqwD3gdcEBHP\nKXKdVWS6kBZnypDKVfTvtdPARzLzE5l5PfB64KfAqQvs/wfADzLzbZn5/cz8IPDp7u20hulCWpop\nQypXYQNERBwATNJ5NQGAzEzgMuCYBQ57evfv59u6yP6NZLqQlseUIZWnyKenRwL7ATv6tu8A1i5w\nzNoF9j8oIh482uVVk+lCGowpQxrez342/LH7j24Z5Zqenmb16tV7bZuammJqaqqkFQ3OdCENrpcy\njj4aNm3qfEl6oJmZGWZmZvbaNju7c+jbK3KA+AmwB1jTt30NcOsCx9y6wP53ZeY9i/1j559/PhMT\nE8OsszJ66eJLXzJdSIPopYxzzoGTToKaPxRIhej/pfqKK+C44+bonG0wuMKepjJzNzALHN/bFhHR\n/fM3Fjjsm/P373pud3uj9dLFuefC4YeXvRqpfkwZ0vLt2gWnnAJHHTX8bRT9e+57gddExMkR8STg\nw8BDgYsAImJzRHx83v4fBg6NiHdGxBMj4g3AS7u301jz08Vb3lL2aqR68qoMafnOOKPzivdZZw1/\nG4UOEJl5MfBHwDnA1cCRwPMy8/buLmuBQ+btfxPwQuDZwDV0Lt88LTP7r8xolF66+NjHTBfSSnhV\nhrS0K66A97+/8/+TX/zF4W8nOldW1ldETACzs7OztTwHYts22LABzj8f3vrWslcj1d/u3bB+PezZ\nA9u3w6pVZa9Iqo5du+DII+HggzvPP9dcM8fk5CTAZGYONHb7+26JTBfS6JkypIX10sUoXvF2gCiR\n6UIqhilDeqD56WIUJ+v7tFUSr7qQiuVVGdL9elddjPIVbweIEpgupOKZMqT7jTJd9DhAlMB0IY2H\nKUMafbro8elrzEwX0niZMtRmRaSLHgeIMTJdSONnylCbFZEuehwgxsh0IZXDlKE2Kipd9Pg0Niam\nC6lcpgy1SZHposcBYgxMF1L5TBlqkyLTRY8DxBiYLqRqMGWoDYpOFz0+nRXMdCFViylDTTaOdNHj\nAFEg04VUPaYMNdk40kWPA0SBTBdSNZky1ETjShc9Pq0VxHQhVZspQ00yznTR4wBRANOFVH2mDDXJ\nONNFjwNEAUwXUj2YMtQE404XPT69jZjpQqoXU4bqrIx00eMAMUKmC6l+TBmqszLSRY8DxAiZLqR6\nMmWojspKFz0+zY2I6UKqN1OG6qTMdNHjADECpgup/kwZqpNeurjwwvJe8XaAGAHThdQMpgzVQS9d\nbN4Mhx1W3jp8ulsh04XULKYMVdn8dPHmN5e7FgeIFTBdSM1jylCVVSFd9DhArIDpQmomU4aqqCrp\nosenvSGZLqRmM2WoSqqULnocIIZgupCaz5ShKqlSuugpbBkR8fMR8TcRsTMi/i0iLoiIA5c45sKI\nuLfv64tFrXFYpgupHUwZqoKqpYueIp/+PgUcARwPvBB4JvCRZRx3KbAGWNv9mipqgcMwXUjtYspQ\nmaqYLnoKGSAi4knA84DTMnN7Zn4DeDPwuxGxdonD78nM2zPztu7XziLWOAzThdQ+pgyVqYrpoqeo\n5RwD/FtmXj1v22VAAk9b4tgNEbEjIq6PiA9FxC8UtMaBmS6kdjJlqAxVTRc9RT0NrgVum78hM/cA\nd3b/biGXAicDvwG8DTgO+GJEREHrXDbThdRupgyNU5XTRc/+g+wcEZuBty+yS9I572EomXnxvD9+\nLyL+AfhnYAPw94sdOz09zerVq/faNjU1xdTUyk+hMF1I6qWMo4/upIxzzil7RWqyXrrYunV0r3jP\nzMwwMzOz17adO4c/SyAyc/k7RzwCeMQSu/0A+D3g3Zl5374RsR9wN/DSzPzcAP/mbcCfZOZHF/j7\nCWB2dnaWiYmJ5d7sQN7yFrjgArj2Wl99kNru7LNh0ya46ioo6CFHLXfFFXDccbBlC2zcWOy/NTc3\nx+TkJMBkZg4U6AZ6BSIz7wDuWGq/iPgm8PCIeMq88yCOBwL41nL/vYh4LJ2B5ZZB1jlKvXRx/vkO\nD5I6vxl+9rOdlLF9O6xaVfaK1CR1SBc9hZwDkZnXA1uBj0bEUyPiGcAHgJnMvLW3X/dEyRO6/31g\nRJwXEU+LiMdFxPHAZ4Eburc1dqYLSf28KkNFqvJVF/2KXN7LgevpXH1xCXAF8Lq+fQ4Heicu7AGO\nBD4HfB/4KPBt4JmZubvAdS7Iqy4k7cu6dXDmmV6VodGq+lUX/QY6B6KKijoHYts22LChky7e+taR\n3aykhti9G9avhz17TBlauV274Mgj4eCDO88/4/qldSXnQPh79T6YLiQtxZShUapTuuipyTLHy3Qh\naTlMGRqFuqWLHp8e+/iGUZIGccYZvsGUhlenqy76OUDMY7qQNChThlaijumip2bLLZbpQtIwTBka\nRl3TRY9Pk12mC0krYcrQIOqcLnocIDBdSFo5U4YGUed00VPTZY+W6ULSKJgytBx1Txc9rX+6NF1I\nGiVThhbTSxfHHlvfdNHT6gHCdCFp1EwZWkwvXTThFe+aL39lTBeSimDK0L40JV30tPZp03QhqUim\nDM3XpHTR08oBwnQhqWimDM3XpHTR05C7MRjThaRxMGUImpcuelr39Gm6kDROpox2a2K66GnVAGG6\nkDRupox2a2K66GnY3Vmc6UJSGUwZ7dTUdNHTmqdR04WkMpky2qXJ6aKnFQOE6UJS2UwZ7dLkdNHT\n0Lu1N9OFpCowZbRD09NFT+OfTk0XkqrElNFsbUgXPY0eIEwXkqrGlNFsbUgXPY2+e6YLSVVkymim\ntqSLnsY+rZouJFWZKaNZ2pQueho5QJguJFWdKaNZ2pQuehp5N00XkurAlNEMbUsXPY17ejVdSKoT\nU0a9tTFd9DRqgDBdSKobU0a9tTFd9BR2dyPijIj4ekTsiog7BzjunIi4OSJ+GhF/FxHLfkHIdCGp\njkwZ9dTWdNFT5NPsAcDFwF8t94CIeDvwJuC1wHpgF7A1IlYtdezsrOlCUn2ZMuqlzemip7ABIjPP\nzsz3Af8wwGEbgU2ZeUlmfhc4GTgYOHGpA88+23Qhqb5MGfXS5nTRU5m7HRFPANYCX+lty8y7gG8B\nxyx1/E9+0u4fpKT6M2XUQ9vTRU+Vnm7XAgns6Nu+o/t3i3rjG00XkurPlFFtvZP125wuevYfZOeI\n2Ay8fZFdEjgiM29Y0aqG8O1vT/PiF6/ea9vU1BRTU1PjXookDa2XMo4+GjZt6nypOk4/HW6+GS69\ntH6veM/MzDAzM7PXtp07dw59e5GZy9854hHAI5bY7QeZ+R/zjnkVcH5m/sISt/0E4J+BozLzO/O2\nXw5cnZnTCxw3AczOzs4yMTGxvDsiSRV39tmd4eGqq8CHtmrYtg02bIAtW2DjxrJXMxpzc3NMTk4C\nTGbmQOFsoFcgMvMO4I5Bjhngtm+MiFuB44HvAETEQcDTgA8W8W9KUlWdcQZ89rOdlLF9O6xa8lo0\nFcl08UBFvg/EIRGxDngcsF9ErOt+HThvn+sj4oR5h20BzoyIF0XEk4FPAP8KfK6odUpSFR1wAFx4\nYeeqDDNG+XyfoQcq8ttwDjAHnAU8rPvfc8DkvH0OB+47cSEzzwM+AHyEztUXPwc8PzM9lUhS6xx1\nVOeqjM2bvSqjTL2PSGj7VRf9BjoHooo8B0JSk+3eDevXw549powy7NoFRx4Jj3kMXH558159WMk5\nEA37VkhSs5gyymW6WJjfDkmqOFNGOUwXi3OAkKQa8A2mxsurLpbmACFJNWDKGC/TxdL8tkhSTZgy\nxsN0sTwOEJJUI6aMYpkuls8BQpJqxJRRLNPF8vntkaSaMWUUw3QxGAcISaohU8ZomS4G5wAhSTVk\nyhgt08Xg/DZJUk2ZMkbDdDEcBwhJqjFTxsqYLobnACFJNWbKWBnTxfD8dklSzZkyhmO6WBkHCElq\nAFPGYEwXK+cAIUkNYMoYjOli5fy2SVJDmDKWx3QxGg4QktQgpozFmS5GxwFCkhrElLE408Xo+O2T\npIYxZeyb6WK0HCAkqYFMGXszXYyeA4QkNZApY2+mi9Hz2yhJDWXK6DBdFMMBQpIarO0pw3RRHAcI\nSWqwtqcM00Vx/HZKUsO1NWWYLorlACFJLdC2lGG6KF5hA0REnBERX4+IXRFx5zKPuTAi7u37+mJR\na5SktmhbyjBdFK/Ib+sBwMXAXw143KXAGmBt92tqxOuSpFZqS8owXYxHYQNEZp6dme8D/mHAQ+/J\nzNsz87bu184i1idJbdT0lGG6GJ8qvrCzISJ2RMT1EfGhiPiFshckSU3R9JRhuhifqn17LwVOBn4D\neBtwHPDFiIhSVyVJDdLUlGG6GK+BBoiI2LyPkxznf+2JiF8edjGZeXFmXpKZ38vMzwP/FVgPbBj2\nNiVJD9S0lGG6GL/9B9z/3cCFS+zzgyHX8gCZeWNE/AQ4DPj7xfadnp5m9erVe22bmppiaspzMCWp\nXy9lPPWpnZRR95zRSxdbt5ouFjIzM8PMzMxe23buHP40w8jMla5p8X8g4lXA+Zk58LkMEfFY4F+A\nEzLzkgX2mQBmZ2dnmZiYWNliJallzj67MzxcdRXU9SF02zbYsAG2bIGNG8teTb3Mzc0xOTkJMJmZ\nAwWtIt8H4pCIWAc8DtgvItZ1vw6ct8/1EXFC978PjIjzIuJpEfG4iDge+CxwA7C1qHVKUpvVPWWY\nLspT5As95wBzwFnAw7r/PQdMztvncKDXHfYARwKfA74PfBT4NvDMzNxd4DolqbXqflWGV12UZ9Bz\nIJYtM08BTllin/3m/ffdwG8WtR5J0r71rsrYtAlOOqk+KaN31cWWLV51UQbnNUlS7VKG6aJ8DhCS\npNqlDNNF+fy2S5KA+rzBlG8YVQ0OEJKk+1Q9ZZguqsMBQpJ0n6qnDNNFdfjtlyTtpaopw3RRLQ4Q\nkqQHqFrKMF1UjwOEJOkBqpYyTBfV449BkrRPVUkZpotqcoCQJC2o7JRhuqguBwhJ0oLKThmmi+ry\nxyFJWlRZKcN0UW0OEJKkJY07ZZguqs8BQpK0pHGnDNNF9fljkSQty7hShumiHhwgJEnLVnTKMF3U\nhwOEJGnZik4Zpov68McjSRpIUSnDdFEvDhCSpIGNOmWYLurHAUKSNLBRpwzTRf34Y5IkDWVUKcN0\nUU8OEJKkoa00ZZgu6ssBQpI0tJWmDNNFffnjkiStyLApw3RRbw4QkqQVGzRlmC7qzwFCkrRig6YM\n00X9+WOTJI3EclOG6aIZHCAqaGZmpuwljFST7k+T7gt4f6qsrvdloZTRuz9NSRd1/fmMUiEDREQ8\nLiIuiIgfRMRPI+KfIuJPI+KAZRx7TkTc3D3u7yKidfNp0/6H2aT706T7At6fKqvrfVkoZfTuT1PS\nRV1/PqNU1I/vSUAArwF+BZgGXg/8+WIHRcTbgTcBrwXWA7uArRGxqqB1SpJGbKGUYbpolkIGiMzc\nmpmnZeZXMvOmzLwEeDfwkiUO3QhsysxLMvO7wMnAwcCJRaxTklSM/pSxZ08z0oXuN84XkB4O3LnQ\nX0bEE4C1wFd62zLzLuBbwDGFr06SNDL9KeO665qRLnS//cfxj3TPY3gT8N8W2W0tkMCOvu07un+3\nkIcAXHfddStZYqXs3LmTuVF+Rm7JmnR/mnRfwPtTZU25L6eeCueeC/feu5M//MM57rprtB8BXpam\n/HzmPXc+ZNBjIzOXv3PEZuDti+ySwBGZecO8Yx4DXA58NTNft8htHwN8DTg4M3fM2/4/gXszc2qB\n414O/M3BoNN1AAAHHklEQVSy74QkSer3isz81CAHDPoKxLuBC5fY5we9/4iIg4GvAl9bbHjoupXO\niZdr2PtViDXA1YsctxV4BXATcPcS/4YkSbrfQ4DH03kuHchAr0AMdMOdVx6+Cnwb+L1cxj8UETcD\n78rM87t/PojOMHFyZv6vQhYqSZIGVtT7QBxMJ1v8C/A24FERsSYi1vTtd31EnDBv0xbgzIh4UUQ8\nGfgE8K/A54pYpyRJGk5RJ1E+Bzi0+/Wj7ragc47EfvP2OxxY3ftDZp4XEQ8FPkLnqo3/DTw/M4f4\nlHlJklSUwhKGJElqLq/GlSRJA3OAkCRJA6v9ABERb4yIGyPiZxFxZUQ8tew1DSMijo2Iz0fEjyPi\n3oh4cdlrGlZEnB4RV0XEXRGxIyI+ExG/XPa6hhURr4+IayNiZ/frGxHxm2WvaxQi4o+7/3t7b9lr\nGUZEnNVd//yvfyx7XSsREQdHxF9HxE+6Hyp4bURMlL2uYXQfm/t/PvdGxAfKXtugIuJBEbFp3odE\n/p+IOLPsda1ERDwsIrZExE3d+/S1iDh6ucfXeoCIiJcB7wHOAp4CXEvnw7ceWerChnMgcA3wBjon\nm9bZscAHgKcBzwYOAL4cET9X6qqG9yM6b6A2AUzSuTz5cxFxRKmrWqHusP1aOv+/qbPv0nm/mLXd\nr18vdznDi4iHA18H7gGeBxwB/CHwb2WuawWO5v6fy1o6J9gncHGZixrSHwOvo/MY/SQ6Vxi+LSLe\nVOqqVuZ/AMfTeS+lXwX+DrgsIh69nINrfRJlRFwJfCszN3b/HHQe7N+fmeeVurgViIh7gRMz8/Nl\nr2UUugPdbcAzM/NrZa9nFCLiDuCPMnOpN1arpIh4GDAL/AHwDuDqzFzsreYrKSLOAk7IzFr+ht4v\nIv4COCYzjyt7LUWIiC3ACzKzdq9IRsQXgFsz8zXztn0a+GlmnlzeyoYTEQ8B/i/wosz80rzt24Ev\nZuZ/X+o2avsKREQcQOe3wfkfvpXAZfjhW1XzcDq/dSz4YWp10X0Z83eBhwLfLHs9K/BB4AuZ+dWy\nFzICh3fT3z9HxCcj4pCyF7QCLwK2R8TF3fw3FxG/X/aiRqH7mP0KOr/11tE3gOMj4nCAiFgHPAP4\nYqmrGt7+dN5W4Z6+7T9jma/ijeXDtArySDp3fl8fvvXE8S9H+9J9VWgLnbczr22bjohfpTMw9Kb2\nkzLz+nJXNZzuAHQUnZeX6+5K4NXA94FHA38KXBERv5qZu0pc17AOpfOq0HuAPwfWA++PiHsy869L\nXdnKnUTnfX8+XvZChvQXwEHA9RGxh84v4H+SmX9b7rKGk5n/LyK+CbwjIq6n89z5cjq/gP/Tcm6j\nzgOE6uFDwK/QmdTr7HpgHZ0HwJcCn4iIZ9ZtiIiIx9IZ6J6dmbvLXs9KZeb89+//bkRcRecdcH+H\npT+3p4oeBFyVme/o/vna7vD6eqDuA8SpwKWZeWvZCxnSy+g8wf4u8I90hvD3RcTNNR7uXgl8DPgx\n8B/AHPApOq/uL6nOA8RPgD10Tp6abw2dD+ZSySLiL4EXAMdm5i1lr2clMvM/uP+D4q6OiPXARjq/\nLdbJJPCfgbnuq0PQeSXvmd2TwR68nM+tqarM3BkRNwCHlb2WId0CXNe37TrgJSWsZWQi4hfpnFB9\nYtlrWYHzgM3zPpfpexHxeOB0ajrcZeaNwLO6J7gflJk7IuJvmfehmIup7TkQ3d+eZumcQQrc93L5\n8XRalUrUHR5OAJ6VmT8sez0FeBDw4LIXMYTLgCfT+e1pXfdrO/BJYF2dhwe47+TQw+g8EdfR13lg\ngn0inVdV6uxUOi+R1/V8Aeic97Snb9u91Ph5tCczf9YdHn6eztU/n13OcXV+BQLgvcBFETELXAVM\n0/khX1TmooYREQfSeeDr/VZ4aPcknTsz80cLH1k9EfEhYAp4MbAr7v8QtZ2ZWbuPXI+Ic4FLgR8C\n/4nOiWDHAc8tc13D6J4XsNe5KBGxC7gjM/t/8628iHgX8AU6T7CPAc4GdgMzZa5rBc4Hvh4Rp9O5\n1PFpwO8Dr1n0qArr/mL3auCizLy35OWsxBfofNjjvwLfo3NZ9zRwQamrWoGIeC6d55zv0/lsqvPo\nPD5ctJzjaz1AZObF3UsEz6GTLq4BnpeZt5e7sqEcDfw9nasVks5JVNA54ejUshY1pNfTuQ+X920/\nhc4nrNbNo+j8HB4N7AS+Azy3IVcwQL3fd+SxdJrtI4Dbga8BT8/MO0pd1ZAyc3tEnETnhL13ADcC\nG+t6ol7Xs4FDqOc5KfO9CdhE5wqmRwE3A3/V3VZXq4HNdIbvO4FPA2dmZv8rLftU6/eBkCRJ5ah9\nu5EkSePnACFJkgbmACFJkgbmACFJkgbmACFJkgbmACFJkgbmACFJkgbmACFJkgbmACFJkgbmACFJ\nkgbmACFJkgb2/wFhgKTVb1IlNQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111507f60>"
      ]
     },
     "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": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhUAAAFdCAYAAAC0B5/iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAACFZJREFUeJzt3TFu5PYdhuEf5UWgMoWaNG5X2/o60VV8G/sKPkZaOZUB\nNWlUpLMQGMMUkm0kWDjR7Gtz/6PnadQshG+BGc4rgkNu+74PAMCnujp6AABwGUQFAJAQFQBAQlQA\nAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAk\nRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAV\nAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAA\nCVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACRE\nBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUA\nQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJ\nUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQF\nAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBA\nQlQAAAlRAQAkRAUAkHh39IDCtm1fzszN0TvO9Ljv+8PRI4DLsvBx0TFxYdu+70dv+GTb9sWPM6fr\no3ec5+pp5vTemwioPAfF1d/XPC46Jq7sIs5UPL9xvpmZD0cPeaX7mbm7nue/JryBgMrNmsdFx8TV\nXUhUzDy/cb46egTAZ8RxkT+WCzUBgISoAAASogIASIgKACAhKgCAhKgAABKiAgBIiAoAIHFBN78C\naC38/IzbowfwNokKgI9Y+/kZcAxRAfBxiz4/Y2bmu5n5+ugRvEGiAuA3rfj8jPujB/BGuVATAEiI\nCgAgISoAgISoAAASogIASIgKACAhKgCAhKgAABKiAgBIiAoAICEqAICEqAAAEqICAEiICgAgISoA\ngISoAAASogIASIgKACAhKgCAhKgAABLvjh4AAP/ldtu2ozec43Hf94ejRxxJVADwmfjHPJ9AP317\n9JLzXD1t2/b+LYeFqADgM/HPmTnNzDcz8+HgLa91PzN31zNzMzOiAgA+Dx9m5qujR3AGF2oCAAlR\nAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQOKCbn51f/SAM/yyecX73F/N863vVvSnmfnX0SPO\nsOrumTVfL7fPP1Y8tvzw8nO17avunln5eL7v+9+q37Xt+179rsNs2xc/zpyuj95xpp9mybi7+mnm\ntODumTU/32bW3T0zXucHWPX1surumVn0db7ve1ZBy/3nP+70fp7vt76iFd9Bt88P/Fnx/vzfzczX\ns972VXfPvDwT4d3M/HVmvj94zCudVnx/vjgtemZr1d0zs+bxPHURUfHyRLg3+wCXP9qvp/ZWvD//\nz6coV9u+6u7/8H15mhX4/LhQEwBIiAoAICEqAICEqAAAEhdxoSa83mrfg//hf/8TgIOJCt6YP8/z\nCbq7o4cAXBxRwRvzl3n5Gvlq90y4nZlvjx4B8FtEBW/VUvdMWO22v8Db5EJNACAhKgCAhKgAABKi\nAgBIiAoAIOHbH3yC1W4gNbPmZoA1iArO8Thz9TRzd330kPNcPc2cHo9eAXBpRAWvtu/7w7Zt72fm\n5ugt5zk97vv+cPQKgEsjKjjLy4eyD2YAfuFCTQAgISoAgISoAAASogIASIgKACAhKgCAhKgAABKi\nAgBIiAoAICEqAICEqAAAEqICAEiICgAgISoAgISoAAASogIASIgKACAhKgCAhKgAABKiAgBIiAoA\nICEqAICEqAAAEqICAEiICgAgISoAgISoAAASogIASIgKACAhKgCAhKgAABKiAgBIiAoAICEqAICE\nqAAAEqICAEiICgAgISoAgISoAAASogIASIgKACAhKgCAhKgAABKiAgBIiAoAICEqAICEqAAAEqIC\nAEiICgAgISoAgISoAAASogIASIgKACAhKgCAhKgAABKiAgBIvDt6APAa90cPOMOKm4FziApYw+PM\n1dPM3fXRQ85z9TRzejx6BfD72vZ9P3oD8H/Ytu3Lmbk5eseZHvd9fzh6BPD7EhUAQMKFmgBAQlQA\nAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAk\nRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAV\nAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAA\nCVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACRE\nBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUA\nQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJ\nUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQF\nAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQOLfh4u1zA6kXBUAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bf5be48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhUAAAFdCAYAAAC0B5/iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAACGBJREFUeJzt3bFuHNcZhuH/CELAIiUrFWrF61Gb61PLW/AdpLbhTjUB\nNwFCBEYmBaUgQeQAXL328Cyfp2FHfAKGy1fLszPrOI4BAPheb84eAABcB1EBACREBQCQEBUAQEJU\nAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACTenj0A/mhrrfczc3v2jgs8\nHMfx+ewRAL9FVPCqrLXez5qf5pibs7c825rHtdYHYQG8VKKC1+Z2jrmZj7PXexUPM3M/N/O0WlQA\nL5Ko4HW6nZl3Z48AuC4OagIACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAk\nRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAV\nAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAA\nCVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACRE\nBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUA\nQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJ\nUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQOLt2QOAZ7lba529\n4RIPx3F8PnsE8PsSFbCDv83MmpljPp095SJrHtdaH4QFXDdRATt4nJljZj7OzO3JW57rYWbu52ae\nlosKuGKiAnZyOzPvzh4B8G0OagIACVEBACREBQCQEBUAQMJBTS6y1no/+30OYWbm7uwBANdKVPBs\na633s+anOebm7C0AvByigkvczjE3W94z4eeZ+eHsEQDXSVRwuR3vmfBw9gCA6+WgJgCQEBUAQEJU\nAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEA\nJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQ\nFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQA\nAAlRAQAk3p49oLDWej8zt2fvuNCbmfnn2SOe6W5mZh5OXnGJX7583W37rrtn/nPz3VrrxCEX2fHn\n86s/zcw/zh5xgV13z2x6vRzH8dfqe63jOKrvdZr1Zv19jrk5e8eFfp0d427Nr3NsuHtmZs3Mjpf9\nrrufuM7/aLteL7vufrLldX4cR1b72/3jv+mYm/k4+71X8TAz9/N2Zv4yMz+evOZ5jj2LfGZmjk3/\nJ7Tr7ic7Xi93c8ynLV9bfp6ZH2a2277r7pm9X89D1xEVM08X4LuzR1zsx/LtJ+D7/ftPNTu+tnz9\nk9Nu23fd/d9e9eu5g5oAQOJ63qnY+wAbAGzvOqJizePcb3pQc83jHPICgP1dR1Qc82H2O9bz5JiH\n4zg+nz0DAL7XVUTFl1/KfjEDwIkc1AQAEqICAEiICgAgISoAgISoAAASV/HpDwCuyI537tlx8+9A\nVADwMtzM01NK788eciE3MxQVALwQf56vjz3f80mfbmYoKgB4cV71kz535qAmAJAQFQBAQlQAAAlR\nAQAkRAUAkBAVAEBCVAAACVEBACREBQCQcEdNgP9nxyc5/HL2AF4rUQHwbQ+z5nHu5+bsIbALUQHw\nDcdxfF5rfZiZ27O3XOBuZj6dPYLXR1QA/IYvT5zc7qmTa62zJ/BKOagJACREBQCQEBUAQEJUAAAJ\nBzUBrtVu99jYbS//Q1QAXJ9977Gx5nEOebGrdRzH2RsAiK213s+e99h4+PJRXjYkKgCAhIOaAEBC\nVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEB\nACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQ\nEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJU\nAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEA\nJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQ\nFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQA\nAAlRAQAkRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAk\nRAUAkBAVAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAV\nAEBCVAAACVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAA\nCVEBACREBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACRE\nBQCQEBUAQEJUAAAJUQEAJEQFAJAQFQBAQlQAAAlRAQAkRAUAkBAVAEBCVAAACVEBACT+BfkXtTkY\ny78wAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11152e128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "emd:  0.706237424547\n"
     ]
    }
   ],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAFkCAYAAACq4KjhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAADFNJREFUeJzt3U2InWcZx+H7mabJJE5tGyfaaDtN/SyYEFJBa60WtRQr\nWkSC4NemCxFERFxY9yq46qaYjUZwYTcBkSBFEemXuLGJQW3aJJJ4mubEZpJmTJrJxySvizFIyx2d\n92Qmz3nnXNfqrA7/zQw/znnP85SmaQIA4I3Gag8AAIaTSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEA\nSIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAl\nEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAlEgCAlEgA\nAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAlEgCAlEgAAFIiAQBI\niQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUS\nAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAA\nUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJ\nBAAgJRIAgJRIAABSIgEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIA\ngJRIAABSIgEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABS\nIgEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkE\nACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAlEgCA\nlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAlEgCAlEgAAFIi\nAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACAlEgCAlEgAAFIiAQBIiQQA\nICUSAICUSAAAUiIBAEitqD0ARlkp5YsR8ZHaOwbwh6ZpHq89AlhapWma2htgJJVSpiLiUESUylMG\n0UTEhqZperWHAEvHJwlQz2RElC2ffX/cMPmm2lsW7NT0a7F7599KzO8XCbCMiQSo7NPfui+mNr+j\n9owF6+15OXbv/FvtGcA14MFFACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIA\ngJRIAABSIgEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABS\nIgEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASK2o\nPQBGXX//sdoTWunaXmBwpWma2htgJJVS7i5j5Y/Npc7+De6IiH/WHjGAPzRN83jtEdAFIgEqKaXc\nFRHPPfS9+2Ny6ubacxas95d+/O7Hz9aecTWaiNjQNE2v9hAYdr5ugMo2ffK9MbX5HbVntLZ165ZY\nt+6G2jNaOXbsVOzYsbtExGREiAT4P0QCMJBHHnkw7rprqvaMVnbt6sWOHbtrz4DO8OsGACAlEgCA\nlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASA3tBU+llKmYv6mta6Zd\nQQvAcjC0kRBRXoxoxmuvaK+cLaW8TygA0HVDHAnN+IYND8fq1etrD1mw2dl+HDq0fTzcVQ/AMjDE\nkRCxevX6WLOmW/fVA8By4cFFACAlEgCAlEgAAFIiAQBIiQQAIDXUv26Yne3XntBK1/YCwP8yxJFQ\n4tCh7bVHDKBERLOy9goAuFpDHAlNTE5+LK6//qbaQxbswoWTMT39dETE+dpbAOBqDXEkRKxb99FO\nHaZ05kzvciQAQOd5cBEASIkEACAlEgCAlEgAAFIiAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIA\ngNRQ390AsEQeK6Wcqj2ipd83TfOj2iMYLSIBGEUfrj1gAA+UUh5vmqZXewijQyQAI+eO2ydjcu1E\n7RkL9tqZc/H8i/2IiMmIEAlcMyIBGDkPf/neePD+TbVnLNgL+/rxla//pPYMRpAHFwGAlEgAAFIi\nAQBIiQQAICUSAICUSAAAUiIBAEiJBAAgJRIAgJRIAABSIgEASIkEACDlgqel8alSyp21RwxgX9M0\nf6o9AoDhIBIW0YULM5df/qDmjqswV0p5l/vqAYgQCYtqbu5MRESM3/vtuO7GWyuvaefizOE4++yj\nK8J99QD8h0hYAivfeV+suGVj7RmtzB39a5x99tHaMwAYIiIBGDlHjp6MF/b1a89YsIO96doTGFEi\nARg527Y/Gdu2P1l7RitjYyUuXWpW1t7BaBEJwMjZ+tAHYvPG22rPWLAjR09ejprzlacwYkQCMHI2\nb7wtHrx/U+0ZC/bCvn7nPvlgeRAJwED27u3Od/qXdXHzGzxWSjlVe8QAft80zY9qj6A9kQC0subN\n4zE2VuKrX/1Z7SkDGRsrMTd3sfaMVqZPnL788sM1d1yFB0opjzuDpXtEAtDKjW+7IS5damLr1i2x\neXN3vtePiDhy5NXYtu2ZWLHiutpTWjl1+mxERNy+/i0xvvL6ymvaOXv+QvyjfzzCGSydJBKAgXzu\nc1viS1/6YO0Zreza1Ytt256pPWNgX37wnphaP1l7Riu9/nT8cPvO2jMYkAueAICUSAAAUiIBAEiJ\nBAAgJRIAgJRIAABSfgIJ0BH94zO1J7TWxc38l0gAGHI3TIxHKSV+9quna08ZSCklmsYNll001JEw\nO9utc9bPnXPnO7D4JtdORNM0seXO2+O2t66tPaeVV0+fiWd2vRjhBstOGuJIKHHo0PbaIwCGxpb3\nTsUHN7279oxWev3py5FABw1xJDSxfv1DsWpVd44gPX36QExPd/PjQAB4oyGOhIibbtoUa9ZM1Z7R\nikgAYLkY6kigiodLKZ+qPaKl55qm+U3tEQDLjUjgjb5Re8AALpVS7nBXPcDiEgm83j33RNxxR+0V\nC3fiRMQTT4yFu+oBFp1I4PXuvjvi4x+vvWLhDhyIeOKJ2isAliXHMgMAKZEAAKREAgCQEgkAQMqD\niwBcC78spZyrPaKlJ5um+VrtETWJBACuhW4dnzvvPaWU74/yGSwiAYAlN7F6daxa1Z3boufmLsbM\n6dMRI34Gi0gAYMl94YFPxIc2vb/2jAXr9Y/GD37689ozqvPgIgCQEgkAQEokAAApkQAApDy4CABX\n9utSyvnaI9pqmub2xXgfkQAAV3ZL7QE1iQQAuIKJiYkYX7Wq9oxW5ubmFu29RAIAXMF3v/XN+PT9\nn6w9o5W9+/Yv2nuJhCVw8fjfa09orYubAVhaImERNc3FiDIWr/36O7WnDKaMRUxM1F4BwJAQCYuo\nlOsimksxfu+347obb609p5ULLz8X5//8i4ibb649BYAhIRKWwMp33hcrbtlYe0Zr5//8i9oTABgi\nIoHl4rullGO1R7S0rvYAgP9FJLBcfKH2gEHN/PNU7QkAKZHAsjBx962xakO3nqeYOzEbM789EGf+\ndbb2FICUSGBZeMvnN8baz9xZe0YrZ55/JWZ+e6D2DIArcsETAJDySQJARxybOR29/nTtGa30j8/U\nnsBVEAkAHbHzqd2x86ndtWe0VkqJNePduv+AeSIBoCM2v+e2ePtb19ae0cqRV07Env0vxY1Oc+2k\noY6E2dl+7QmtnDs3/zFgF+9BuDhzeP7F4cN1h7T1n73nDs/EmedfqTymnbN/PxEREdO9V6O35+XK\naxauv3/+OIqDB4/Frl29ymva2bt3/n/KkaMn44V93fn/crA3/79lz/6XYs/+lyqvGUx/+njtCa1c\n3nukf3RRL0y6Fg7+Y/H+LkvTNIv2ZgDA8uHXDQBASiQAACmRAACkRAIAkBIJAEBKJAAAKZEAAKRE\nAgCQEgkAQEokAAApkQAApEQCAJASCQBASiQAACmRAACkRAIAkBIJAEBKJAAAKZEAAKREAgCQEgkA\nQEokAAApkQAApEQCAJASCQBASiQAACmRAACkRAIAkBIJAEDq3yBO6U7Tf1MKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111191d68>"
      ]
     },
     "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": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAFkCAYAAACq4KjhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAADE9JREFUeJzt3UGIXed9xuH/NyPkdJAtCRRVasMYpY7slJi4IrQEx46L\nQ3ECNaUNWZRmUUMXbqBp6SK4NCtjFy8Suig1ZGFDS5tCN00g0CzcJHK8iKnVlrYosuMqGmyP7MxY\nUjweS7JmThdmWmRetXNGI39zR8+zuqvLq4sEP937nXPaMAwFAPBuU70HAABbk0gAACKRAABEIgEA\niEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgA\nACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQS\nAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJ\nBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABE\nIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAA\nkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkA\nQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQC\nABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKR\nAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBI\nJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAg\nEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEA\niEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgA\nACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBoR+8BcLVaaw9W1T29d2zQd4dheLz3CICkDcPQewNs\nWGtttqpO9d5xlW4ehmGu9wiAd/NNApNuX1XVvffeWnv3zvTeMsqZM8v11FMnqt75M4gEYMsRCWwL\njz32W3XkyGzvGaMcOzZXH/vYo71nAFyRg4sAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAg\nEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEA\niEQCABCJBAAgEgkAQCQSAIBIJAAAkUgAACKRAABEIgEAiEQCABDt6D0AqL9orb3Re8RI/zQMw2O9\nRwDXlkiA/j7ee8AG/Fpr7evDMMz1HgJcOyIBOjtw6/tr9/4be89Yt/NLF+rUv7xcVbWvqkQCbGMi\nATr7zB/eU7/y2Tt6z1i3uX97uR751F/2ngG8BxxcBAAikQAARCIBAIhEAgAQiQQAIBIJAEAkEgCA\nSCQAAJFIAAAikQAARCIBAIhEAgAQiQQAIPIUSP5Ha+3Pq+pXe+8Y6Wd6DwDYrkQCVVXVWputqi/2\n3gHA1iESWLOvqqr276/aubPzlBHefLPqzJneKwC2JZHA5R56qOqWW3qvWL/vfKfqq1/tvQJgW3Jw\nEQCIRAIAEIkEACASCQBAJBIAgEgkAACRSAAAIpEAAERupsS2cPz4fO8Jo03iZuD6IhKYbCsrNTXV\n6vOff7L3kg2Zmmo1c9P7es8AiEQCk216ulZXh7r/oU/Vvtm9vdeM8qNnT9XRJ5+t3T97Y+8pAJFI\nYFu4/d7DNfvRn+89Y7SjTz7bewLAFTm4CABEIgEAiPzcAGzUA621+3qP2IDnhmH4du8RMAlEAjDK\nuVffWHv5hZ47rsJqa+3QMAxzvYfAVicSgFGWf3q+qqruOnJr7d0103nNOGeWluvpYyemqmpfVYkE\n+H+IBGBD7rrjcM0e3Nd7xihz8wv19LETvWfAxHBwEQCIRAIAEIkEACASCQBAJBIAgEgkAACRSyCB\nDZlfPNd7wmiTuBl6EgnAKCsXV6q1Vk9+42jvKRvSWqthGHb23gGTQCQAo0zvnK5hGOrBB+6pnzuw\np/ecUV45fbYef+K7VVUXO0+BiSASgA2585dvqdsOH+w9Y5QfPj+/FgnAOji4CABEvkngci+91HvB\nOK+91nsBwLYlElizs1qr+spXeu8AYIsQCay5WMNQ7/vEH9X07g/03rJub7/8XF3817/tPQNgWxIJ\nXGbnBz9ZOw58pPeMUUQCwLXh4CIAEIkEACASCQBAJBIAgEgkAACRSAAAIpEAAETukwBcj77UWvtJ\n7xEjPTMMw9d7j+D6IhK4zMrii70njLJybsKeNUFXC68vrb38XM8dG/T7rbVnhmGY6z2E64dIYM3O\nalP15rf+uPcOuGbeWDpfVVX33v3h2rtnpvOa9TtzdrmeOnq8VdW+qhIJvGdEAmsu1rBaBw/eXzfc\nsK/3lnVbWvpRLSwc7T2DCfO7v31n3Xb4YO8Z6/bD5+frqaPHe8/gOiQSuMyePbfXzMxs7xmjiIQ+\nTs4t9J4w2iunz/aeABNFJACjrFxcqampVl9+9B96TwGuMZEAjDK9c7pWV4d6+OFfr0OH3t97zijP\nPPNCPf74071nwMQQCcCGfPrTt9eRI5P101RViQQYQSSwLcy/MGmXvFctzJ3pPQHg/yQSNllr7WtV\ndU/vHRtwQ+8BGzEMK9WmWj3x4N/3ngKw7YiETdRam62q3+u942q8/fa53hNGaW26htWhDvzBx+uG\nD+zuPWeUpWMv1+Lf/XvvGQBXJBI21zs3GNi1v2p6Z+cpI61crFp6rS5dWu69ZEN2332oZn5xf+8Z\no4kEYCsTCdfATb/5tdpx4CO9Z4xy6fR/1E//6jd6zwBgC/EUSAAgEgkAQCQSAIDImQS2hfMvvt57\nwmgXXpqsK0mA649IYKINw0rVVKtTX/rH3lMAth2RwERrbbpqdag9e36pduy4sfecUc6ff7WWlk70\nnsFkeaC1dl/vERvw3DAM3+49gvFEAtvCwYOfmbhHXC8u/kAksC4Lry+tvfxCzx1XYbW1dmgYhrne\nQxhHJABscW8sna+qqruO3Fp7d810XjPOmaXlevrYial652ZzImHCiASACXHXHYdr9uC+3jNGmZtf\nqKeP+cZsUrkEEgCIRAIAEIkEACASCQBA5OAiAO+F+1prt/UeMdLzwzD8c+8RPYkEAK6Zc0vLay8f\n6bljgy611n7her6/g0gA4Jp56/zFqqq6/5OfqH17dndes34LZ8/VN7/3/R11nd/fQSQAcM3dfssH\na/bggd4z1m1u/nR983vf7z2jOwcXAYBIJAAAkUgAACKRAABEIgEAiFzdAABXdltrrfeG0YZhOLYZ\n7yMSAOBdzi0t1dTUVK2urv5N7y0btCllIxIA4F2Wz1+o1dXVevhPHqpDN8/2njPKyVObd+8nkQAA\nV3Do5tn68OEP9Z7RjYOLAEAkEgCASCQAAJFIAAAikQAARFv26obW2l9X1d29d4y0s/cAANgsWzYS\nqup3eg8AgOvZVo6Emp7eVa1N956xbqurF2t19a3eMwBgU2zpSDh8+Is1MzM5d7paXPxB/fjHT/Se\nAQCbwsFFACASCQBAtKV/bgBge5hfWOw9YZSFs+d6T9gSRAIA18ylldVqrdUT3/hW7ykbsplPVHyv\nXDdPgXzrrfneE0a5cGGhqqpWFl/svGS8tc0XLizU8vLk/KNY+8wn7e9K1f9un3/hJ52XjLMwd6aq\nqo4fn7zP/OTJdz7rk3MLnZeM88rps1VVNb84ef+7PbO0XMMw1J0f/VDdtGum95x1e+nVxfrP/3ql\nvvzon/WesiF/+sijm/I+bRiGTXkjAGB7cXARAIhEAgAQiQQAIBIJAEAkEgCASCQAAJFIAAAikQAA\nRCIBAIhEAgAQiQQAIBIJAEAkEgCASCQAAJFIAAAikQAARCIBAIhEAgAQiQQAIBIJAEAkEgCASCQA\nAJFIAAAikQAARCIBAIhEAgAQiQQAIBIJAEAkEgCA6L8ByNXUNzo3uEgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1125c9ba8>"
      ]
     },
     "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()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}