COT_1b_monotone_rearrangement.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# MATH3349 / 4349 Special Topics: Computational Optimal Transport\n",
    "\n",
    "## Assignment 1, part B: Optimal Transport with monotone rearrangment\n",
    "\n",
    "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n",
    "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n",
    "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n",
    "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n",
    "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n",
    "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n",
    "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n",
    "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n",
    "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n",
    "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$\n",
    "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n",
    "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n",
    "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n",
    "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n",
    "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n",
    "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n",
    "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n",
    "$\\newcommand{\\CC}{\\mathbb{C}}$\n",
    "$\\newcommand{\\RR}{\\mathbb{R}}$\n",
    "$\\newcommand{\\EE}{\\mathbb{E}}$\n",
    "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n",
    "$\\newcommand{\\Ww}{\\mathcal{W}}$\n",
    "$\\newcommand{\\Vv}{\\mathcal{V}}$\n",
    "$\\newcommand{\\Nn}{\\mathcal{N}}$\n",
    "$\\newcommand{\\NN}{\\mathcal{N}}$\n",
    "$\\newcommand{\\Hh}{\\mathcal{H}}$\n",
    "$\\newcommand{\\Bb}{\\mathcal{B}}$\n",
    "$\\newcommand{\\Ee}{\\mathcal{E}}$\n",
    "$\\newcommand{\\Cc}{\\mathcal{C}}$\n",
    "$\\newcommand{\\Gg}{\\mathcal{G}}$\n",
    "$\\newcommand{\\Ss}{\\mathcal{S}}$\n",
    "$\\newcommand{\\Pp}{\\mathcal{P}}$\n",
    "$\\newcommand{\\Ff}{\\mathcal{F}}$\n",
    "$\\newcommand{\\Xx}{\\mathcal{X}}$\n",
    "$\\newcommand{\\Mm}{\\mathcal{M}}$\n",
    "$\\newcommand{\\Ii}{\\mathcal{I}}$\n",
    "$\\newcommand{\\Dd}{\\mathcal{D}}$\n",
    "$\\newcommand{\\Ll}{\\mathcal{L}}$\n",
    "$\\newcommand{\\Tt}{\\mathcal{T}}$\n",
    "$\\newcommand{\\si}{\\sigma}$\n",
    "$\\newcommand{\\al}{\\alpha}$\n",
    "$\\newcommand{\\la}{\\lambda}$\n",
    "$\\newcommand{\\ga}{\\gamma}$\n",
    "$\\newcommand{\\Ga}{\\Gamma}$\n",
    "$\\newcommand{\\La}{\\Lambda}$\n",
    "$\\newcommand{\\si}{\\sigma}$\n",
    "$\\newcommand{\\Si}{\\Sigma}$\n",
    "$\\newcommand{\\be}{\\beta}$\n",
    "$\\newcommand{\\de}{\\delta}$\n",
    "$\\newcommand{\\De}{\\Delta}$\n",
    "$\\newcommand{\\phi}{\\varphi}$\n",
    "$\\newcommand{\\th}{\\theta}$\n",
    "$\\newcommand{\\om}{\\omega}$\n",
    "$\\newcommand{\\Om}{\\Omega}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This project we solve the discrete optimal transport problem \n",
    "problem (with measures that are sums of Diracs) using linear\n",
    "programming."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.optimize\n",
    "import scipy.sparse\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: Discrete distributions in 1 dimension\n",
    "## Using monotone rearrangement\n",
    "\n",
    "Now we're going to look at computing the monotone rearrangement for an arbitrary 1 dimensional point cloud.\n",
    "\n",
    "We're doing this _after_ having considered the 2 dimensional problem as we can test our monotone rearrangement agains the calculation from the linear programming method.\n",
    "\n",
    "Let's define two point clouds in $\\RR$ and two vectors of weights. We'll then explicitly build the transport coupling using what we know. Again we have\n",
    "\n",
    "$$\n",
    "\\alpha = \\sum_{i=1}^n a_i \\delta_{x_i} \\quad\\text{and}\\quad \n",
    "\\beta = \\sum_{j=1}^m b_j \\delta_{y_j}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 8\n",
    "m = 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(2)\n",
    "\n",
    "# Generate normally distributed, or sum of normally distributed points. \n",
    "X = np.random.randn(int(n))\n",
    "Y = np.concatenate((1. * np.random.randn(int(m/2)) + 3, 1.*np.random.randn(int(m/2)) + 5))\n",
    "\n",
    "a = np.random.random(n)\n",
    "a /= a.sum()\n",
    "b = np.random.random(m)\n",
    "b /= b.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A crucial step is to order the points spatially - which will make the calculation of inverse cumulatives easier. We label the cumulative density functions (CDFs) as ```a_cdf``` and ```b_cdf```, noting that we use the ```np.cumsum()```. We also make use of ```np.argsort```, which gives us the order to sort the ```X``` and ```Y``` vectors, which we can apply to ```a``` and ```b```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_sort = np.argsort(X)\n",
    "Y_sort = np.argsort(Y)\n",
    "\n",
    "X = X[X_sort]\n",
    "Y = Y[Y_sort]\n",
    "\n",
    "a = a[X_sort]\n",
    "a_cdf = a.cumsum()\n",
    "\n",
    "b = b[Y_sort]\n",
    "b_cdf = b.cumsum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's have a look at the distributions and the corresponding cumulative density function (CDF)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3AAAAJOCAYAAAD27eW+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABQjUlEQVR4nO3df5hV1Z3n+/e3qLKkqDJCUNqAFhjRiEyPKCZRO44dmbZCt9EQb2PUTjeTjNOXZHTozNgmPTedcZ6naTvTTSe3ZSbGhNvTwQmOoXOZG8FoiHe8oi2lmDTIkCiClr/4VcaiCuvnun/sQgqogoJTdc7Zdd6v56mnzl5rn7O/B45Sn1p7rRUpJSRJkiRJ5a+q1AVIkiRJkobHACdJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJUkWLiI9FxLYiXetrEfG9YlxLkjQ2GeAkSWUhInZExIGIaIuItyNiQ0T8YUSM6r9VKaUnUkoXHFHHvNG85nCUSx2SpPJigJMklZPrUkoNQCPw58AfA98pbUmSJJUPA5wkqeyklH6VUloDLAR+PyJmA0REbUT8p4h4JSLeioj/EhHj+/uujoiWiPhSROyKiDciYtHB14yI+RHxQv8I32sR8W8HPq//8d8B5wD/IyL2R8SdEfGjiPjXA+uLiJ9HxKeOrDsipkdEiojbIuL1/hr+7VDvMyI+GRFb+kccH4+IC4eqo7A/UUnSWGGAkySVrZTSM0AL8LH+pj8HzgcuBs4DpgJfHfCUXwPe19/+OeDeiJjY3/cd4F/1j/DNBtYPcr3fA14hGwmsTyn9BfC3wK0Hz4mIf9r/+j86Rum/CcwEfgv448FuhYyI84H/Bvwb4AzgYbLAdsoQdUiSZICTJJW914FJERHAbcCSlNK+lFIb8GfATQPO7QbuTil1p5QeBvYDFwzomxURp6WUWlNKzw3z+muA8yNiZv/x7wGrUkpdx3jOf0gptaeU/hFYAXxmkHMWAj9KKT2aUuoG/hMwHrhimHVJkiqQAU6SVO6mAvvIRqnqgGf7bzl8G1jX337Q3pRSz4DjDqC+//GngfnAzoj4fyPi8uFcPKX0LrAKuLV/QZXPAH93nKe9OuDxTuADg5zzgf6+g9fp63/e1OHUJUmqTAY4SVLZiojLyALN/wfsAQ4AF6WUTu//el9Kqf6YL9IvpbQxpXQ9cCbwQ+DBoU4dpO1vgVuAa4COlNJTx7nc2QMen0M2inik18kWawGgf4TxbOC1Y9QhSapwBjhJUtmJiNMi4neA7wPfSyn9Y/8I1beBZRFxZv95UyPi2mG83ikRcUtEvK//dsV3gL4hTn8LOHdgQ39g6wP+kuOPvgH8HxFRFxEXAYvIRvCO9CDw2xFxTUTUAF8COoENQ9UhSZIBTpJUTv5HRLSR3Ur4J8BfkQWgg/4YeBF4OiLeAR7j0By34/k9YEf/8/6QbERtMEuBf99/m+bAFST/K/BPgOFsxP3/9tf5E+A/pZR+fOQJKaVtZIuj/J9ko4vXkS1acnBu3VB1SJIqWKTkHRqSJB1PRHwWuC2l9BvHOGc68DJQc8RcPEmSRoQjcJIkHUdE1AGLgftKXYskqbIZ4CRJOob+OXa7yeakPVDiciRJFc5bKCVJkiQpJxyBkyRJkqScqC51AUeaPHlymj59eqnLkCRJkqSSePbZZ/eklM4YrK/sAtz06dNpbm4udRmSJEmSVBIRsXOoPm+hlCRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJkiRJyomyW4VSklQ5Oro7WLF5Bau2reLtzrc5vfZ0Fl6wkEWzF1FXU1fq8iRJKjsGOElSSXR0d3DzwzfT0tZCZ28nAK2drazYsoJHX3mUB+Y/YIiTJOkI3kIpSSqJFZtXHBbeDurs7aSlrYUVm1eUqDJJksqXAU6SVBKrtq06Krwd1Nnbyaptq4pckSRJ5c8AJ0kqibc73y6oX5KkSmSAkySVxOm1pxfUL0lSJTLASZJKYuEFC6kdVztoX+24WhZesLDIFUmSVP4McJKkklg0exHTGqYdFeJqx9UyrWEai2YvKlFlkiSVLwOcJKkk6mrqeGD+Ayy6aBHVke1qM7F2IosuWuQWApIkDcF94CRJJVNXU8cX5nyB5reaAVjR5NYBkiQdiyNwkiRJkpQTBjhJkiRJygkDnCRJkiTlhAFOkiRJknLCACdJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJkiRJygkDnCRJkiTlhAFOkiRJknLCACdJkiRJOWGAkyRJkqScGFaAi4imiNgWES9GxF2D9P9RRLwQET+PiJ9EROOAvt6IeL7/a81IFi9JkiRJlaT6eCdExDjgXuCfAy3AxohYk1J6YcBpm4C5KaWOiPjfgb8AFvb3HUgpXTyyZUuSJElS5RnOCNyHgRdTSttTSl3A94HrB56QUvppSqmj//BpYNrIlilJkiRJGk6Amwq8OuC4pb9tKJ8D1g44PjUimiPi6Yi4YbAnRMRt/ec07969exglSZIkSVLlOe4tlCciIm4F5gL/bEBzY0rptYg4F1gfEf+YUnpp4PNSSvcB9wHMnTs3jWRNkiRJkjRWDGcE7jXg7AHH0/rbDhMR84A/AT6ZUuo82J5Seq3/+3bgcWBOAfVKkiRJUsUaToDbCMyMiBkRcQpwE3DYapIRMQf4Fll42zWgfWJE1PY/ngxcCQxc/ESSJEmSNEzHvYUypdQTEV8EHgHGAd9NKW2JiLuB5pTSGuDrQD3w3yMC4JWU0ieBC4FvRUQfWVj88yNWr5QkSZIkDdOw5sCllB4GHj6i7asDHs8b4nkbgH9SSIGSJEmSpMywNvKWJEmSJJWeAU6SJEmScsIAJ0mSJEk5YYCTJEmSpJwwwEmSJElSThjgJEmSJCknDHCSJEmSlBMGOEmSJEnKCQOcJEmSJOWEAU6SJEmScsIAJ0mSJEk5YYCTJEmSpJwwwEmSJElSThjgJEmSJCknDHCSJEmSlBMGOEmSJEnKCQOcJEmSJOWEAU6SJEmScsIAJ0mSJEk5YYCTJEmSpJwwwEmSJElSThjgJEmSJCknDHCSJEmSlBMGOEmSJEnKCQOcJEmSJOWEAU6SJEmScsIAJ0mSJEk5YYCTJEmSpJwwwEmSJElSThjgJEmSJCknDHCSJEmSlBPDCnAR0RQR2yLixYi4a5D+P4qIFyLi5xHxk4hoHND3+xHxy/6v3x/J4qWxqKO7g3s33ctV37+KX//bX+eq71/FvZvupaO7o9SlSZIkqcSOG+AiYhxwL/AJYBbwmYiYdcRpm4C5KaVfBx4C/qL/uZOAPwU+AnwY+NOImDhy5UtjS0d3Bzc/fDMrtqygtbOVRKK1s5UVW1Zw88M3G+IkSZIq3HBG4D4MvJhS2p5S6gK+D1w/8ISU0k9TSgd/snwamNb/+Frg0ZTSvpRSK/Ao0DQypUtjz4rNK2hpa6Gzt/Ow9s7eTlraWlixeUWJKpMkSVI5GE6Amwq8OuC4pb9tKJ8D1p7IcyPitohojojm3bt3D6MkaWxatW3VUeHtoM7eTlZtW1XkiiRJklRORnQRk4i4FZgLfP1EnpdSui+lNDelNPeMM84YyZKkXHm78+2C+iVJkjS2DSfAvQacPeB4Wn/bYSJiHvAnwCdTSp0n8lxJmdNrTy+oX5IkSWPbcALcRmBmRMyIiFOAm4A1A0+IiDnAt8jC264BXY8AvxURE/sXL/mt/jZJg1h4wUJqx9UO2lc7rpaFFywsckWSJEkqJ8cNcCmlHuCLZMFrK/BgSmlLRNwdEZ/sP+3rQD3w3yPi+YhY0//cfcB/JAuBG4G7+9skDWLR7EVMa5h2VIirHVfLtIZpLJq9qESVSZIkqRxUD+eklNLDwMNHtH11wON5x3jud4HvnmyBUiWpq6njgfkPsGLzCu7/x/vpST1MrJ3IwgsWsmj2Iupq6kpdoiRJkkpoWAFOUvHU1dTxhTlfoPmtZgBWNLl1gCRJkjIjugqlJEmSJGn0GOAkSZIkKScMcJIkSZKUEwY4SZIkScoJA5wkSZIk5YQBTpIkSZJywgAnSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSZIkKScMcJIkSZKUEwY4SZIkScoJA5wkSZIk5YQBTpIkSZJywgAnSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSZIkKScMcJIkSZKUE9WlLkBjX0d3Bys2r2DVtlW83fk2p9eezsILFrJo9iLqaupKXZ4kSZKUGwY4jaqO7g5ufvhmWtpa6OztBKC1s5UVW1bw6CuP8sD8BwxxkiRJ0jB5C6VG1YrNKw4Lbwd19nbS0tbCis0rSlSZJEmSlD8GOI2qVdtWHRXeDurs7WTVtlVFrkiSJEnKLwOcRtXbnW8X1C9JkiTpEAOcRtXptacX1C9JkiTpEAOcRtXCCxZSO6520L7acbUsvGBhkSuSJEmS8ssAp1G1aPYipjVMOyrE1Y6rZVrDNBbNXlSiyiRJkqT8McBpVNXV1PHA/AdYdNEiqiPbtWJi7UQWXbTILQQkSZKkE+Q+cBp1dTV1fGHOF2h+qxmAFU1uHSBJkiSdDAOcJEmSTkznftjwTdh4P3Tsg7pJcNnn4Yrboba+1NVJY5oBTpIkScPXuR/unwetL0PPu1lbx1548hvwwhr4/GOGOGkUDWsOXEQ0RcS2iHgxIu4apP+qiHguInoi4sYj+noj4vn+rzUjVbgkSZJKYMM3Dw9vB/W8m7Vv+GZp6pIqxHEDXESMA+4FPgHMAj4TEbOOOO0V4A+ABwZ5iQMppYv7vz5ZYL2SJEkqpY33Hx3eDup5FzZ+p7j1SBVmOLdQfhh4MaW0HSAivg9cD7xw8ISU0o7+vr5RqFGSJEnlomPfsfsP7C1OHVKFGs4tlFOBVwcct/S3DdepEdEcEU9HxA2DnRARt/Wf07x79+4TeGlJkiQVVd2kY/ePf39x6pAqVDH2gWtMKc0Fbgb+OiI+eOQJKaX7UkpzU0pzzzjjjCKUJEmSpJNy2eeh+tTB+6pPhcs+V9x6pAoznAD3GnD2gONp/W3DklJ6rf/7duBxYM4J1CdJkqRycsXtMHHG0SGu+tSs/YrbS1OXVCGGE+A2AjMjYkZEnALcBAxrNcmImBgRtf2PJwNXMmDunCRJknKmtj7bKuDKO6CqJmurm5wdu4WANOqOu4hJSqknIr4IPAKMA76bUtoSEXcDzSmlNRFxGfD3wETguoj4Dymli4ALgW/1L25SBfx5SskAJ0mSlGe19fCbX4EdT2bHi35U2nqkCjKsjbxTSg8DDx/R9tUBjzeS3Vp55PM2AP+kwBolSZIkSRRnERNJkiRJ0ggwwEmSJElSThjgJEmSJCknDHCSJEmSlBMGOEmSJEnKCQOcJEmSJOWEAU6SJEmScsIAJ0mSJEk5YYCTJEmSpJwwwEmSJElSTlSXugBJkiRJZaBzP2z4Jmy8Hzr2Qd0kuOzzcMXtUFtf6urUzwAnSZIkVbrO/XD/PGh9GXrezdo69sKT34AX1sDnHzPElQlvoZQkSZIq3YZvHh7eDup5N2vf8M3S1KWjGOAkSZKkSrfx/qPD20E978LG7xS3Hg3JACdJkiRVuo59x+4/sLc4dei4DHCSJElSpaubdOz+8e8vTh06LgOcJEmSVOku+zxUnzp4X/WpcNnniluPhmSAkyRJkirdFbfDxBlHh7jqU7P2K24vTV06igFOkiRJqnS19dlWAVfeAVU1WVvd5OzYLQTKivvASZIkScpC2m9+BXY8mR0v+lFp69GgHIGTJEmSpJwwwEmSJElSTngLpSRJkqTK0bkfNnwz27y8Y1+2hcJln88WasnBXD8DnCRJkqTK0Lkf7p8HrS9Dz7tZW8deePIb8MKaXCzY4i2UkiRJkirDhm8eHt4O6nk3a9/wzdLUdQIMcJIkSZIqw8b7jw5vB/W8Cxu/U9x6ToK3UJaxju4OVmxewaptq3i7821Orz2dhRcsZNHsRdTV1JW6PEmSJClfOvYdu//A3uLUUQBH4MpUR3cHNz98Myu2rKC1s5VEorWzlRVbVnDzwzfT0d1R6hIlSZKkfKmbdOz+8e8vTh0FMMCVqRWbV9DS1kJnb+dh7Z29nbS0tbBi84oSVSblS0d3B/duupervn8Vv/63v85V37+Kezfd6y9BJEmqRJd9HqpPHbyv+lS47HPFreckGODK1Kptq44Kbwd19nayatuqIlck5Y8j2ZIk6TBX3A4TZxwd4qpPzdqvuL00dZ0AA1yZervz7YL6JTmSLWkUde6Hn/4Z/MW58LXTs+8//bOsXVL5qq3Ptgq48g6oqsna6iZnxznYQgCGGeAioikitkXEixFx1yD9V0XEcxHRExE3HtH3+xHxy/6v3x+pwse602tPL6hfkiPZkkbJwX2knvxGtn8U6dA+UvfPM8RJ5a62Hn7zK3D2R6DxN+DOl7LjHIQ3GEaAi4hxwL3AJ4BZwGciYtYRp70C/AHwwBHPnQT8KfAR4MPAn0bExMLLHvsWXrCQ2nG1g/bVjqtl4QULi1yRlD+OZEsaFWNgHylJ+TWcEbgPAy+mlLanlLqA7wPXDzwhpbQjpfRzoO+I514LPJpS2pdSagUeBZpGoO4xb9HsRUxrmHZUiKsdV8u0hmksmr2oRJVJ+eFItqRRMQb2kZKUX8MJcFOBVwcct/S3DcewnhsRt0VEc0Q07969e5gvPbbV1dTxwPwHWHTRIqoj265vYu1EFl20iAfmP+A+cNIwOJItaVSMgX2kJOVXWWzknVK6D7gPYO7cuanE5ZSNupo6vjDnCzS/1QzAiiYXXJBOxKLZi3j0lUePWsjEkWxJBamb1D/3bQg52EdKUn4NZwTuNeDsAcfT+tuGo5DnSlJBHMmWNCrGwD5SkvJrOCNwG4GZETGDLHzdBNw8zNd/BPizAQuX/Bbw5ROuUpJOkiPZkkbcFbfDC2uOXsgkR/tIScqv447ApZR6gC+ShbGtwIMppS0RcXdEfBIgIi6LiBbgfwO+FRFb+p+7D/iPZCFwI3B3f5skSVI+jYF9pErG/fOkgg1rDlxK6WHg4SPavjrg8Uay2yMHe+53ge8WUKMkSVJ5ObiP1I4ns+NFPyptPXlwcP+8gSOXB/fPe2GN4VcapmFt5C1JkiQVxP3zpBFhgJMkSdLoc/88aUQY4CRJkjT63D9PGhEGOEmSJI2+uknH7nf/PGlYDHCSJEkafe6fJ40IA5wkSZJG3xW3Z/vkHRni3D9POiEGOEmSJI0+98+TRsSw9oGTJEkaFZ37s+XjN96fLXJRNym71e6K2/2Bfixy/zypYAY4SZJUGm7srLHGX0ioCLyFUpIklYYbO2ssOfgLiSe/kf0ignToFxL3z8v6pRFggJMkSaXhxs4aS4r5C4nO/fDTP4O/OBe+dnr2/ad/ZkisEAY4SZJUGm7srLGkWL+QcKSv4hngJEmVw99alxc3dtZYUqxfSHjrccUzwEmSKoO/tS4/buyssaRYv5Dw1uOKZ4CTJFUGf2tdftzYWWNJsX4h4a3HFc8AJ0mqDP7Wuvy4sbPGkmL9QsJbjyueAU6SVBn8rXV5Orix89kfgcbfgDtfyo4Nb8qbYv1CwluPK54beUuSKkPdpP65b0Pwt9aSCnXwFxI7nsyOF/1o5K9xxe3ZRvdH3hLurccVwxE4SVJl8LfWksYCbz2ueI7ASZIqg7+1ljRWFGOkT2XLEThJUmXwt9aSpDHAEThJUuXwt9aSpJwzwEmSxp7O/dm+bhvvz1afrJuUzYG74nZH2iRJuWaAO4aO7g5WbF7Bqm2reLvzbU6vPZ2FFyxk0exF1NXUlbo8SdJgOvfD/fMOn+vWsRee/EY2B+7zj5W2PkmSCmCAG0JHdwc3P3wzLW0tdPZ2AtDa2cqKLSt49JVHeWD+A4Y4SSpHG7559EIlkB23vpz1S5KUUy5iMoQVm1ccFt4O6uztpKWthRWbV5SoMknSMW28/+jwdlDPu7DxO8WtR5KkEWSAG8KqbauOCm8HdfZ2smrbqiJXJEkalo59x+4/cIzNvCVJKnMGuCG83fl2Qf2SpBKpm3Ts/vHvL04dkiSNAgPcEE6vPb2gfklSiVz2+Wxz7sFUnwqXfa649UiSNIIMcENYeMFCasfVDtpXO66WhRcsLHJFkqRhueJ2mDjj6BBXfWrWfsXtpalLkqQRYIAbwqLZi5jWMO2oEFc7rpZpDdNYNHtRiSqTJB1TbX22VcCVd0BVTdZWNzk7/vxj7gMnSco1txEYQl1NHQ/Mf4AVm1dw/z/eT0/qYWLtRPeBGyPc408a42rr4Te/AjuezI4X/ai09UiSNEKGFeAiogn4BjAOuD+l9OdH9NcC/xW4FNgLLEwp7YiI6cBWYFv/qU+nlP5whGofdXU1dXxhzhdofqsZgBVNbh0wFrjHnyRJkvLquLdQRsQ44F7gE8As4DMRMeuI0z4HtKaUzgOWAfcM6HsppXRx/1duwpvGLvf4kyRJUl4NZw7ch4EXU0rbU0pdwPeB648453rgb/sfPwRcExExcmVKI8c9/iRJkpRXw7mFcirw6oDjFuAjQ52TUuqJiF8BBzfamRERm4B3gH+fUnriyAtExG3AbQDnnHPOCb0B6US5x9/xOUdQkiSpPI32KpRvAOeklOYAfwQ8EBGnHXlSSum+lNLclNLcM844Y5RLUqVzj79jOzhHcMWWFbR2tpJI780RvPnhm+no7ih1iZIkSRVrOAHuNeDsAcfT+tsGPSciqoH3AXtTSp0ppb0AKaVngZeA8wstWiqEe/wdm3MEJUmSytdwAtxGYGZEzIiIU4CbgDVHnLMG+P3+xzcC61NKKSLO6F8EhYg4F5gJbB+Z0qWT4x5/x+YcQUmSpPJ13ACXUuoBvgg8QrYlwIMppS0RcXdEfLL/tO8A74+IF8lulbyrv/0q4OcR8TzZ4iZ/mFLaN8LvQTohB/f4W3TRIqojmwY6sXYiiy5a5BYCOEdQkiSpnA1rH7iU0sPAw0e0fXXA43eB/22Q5/0A+EGBNUojzj3+hnZ67em0drYes1+SJEmlMdqLmEjKGecISpIklS8DnKTDOEdQkiSpfBngJB3GOYKSJEnla1hz4CRVFucISpIklScDnCRJx9O5HzZ8EzbeDx37oG4SXPZ5uOJ2qK0vdXWSpApigJMk6Vg698P986D1Zeh5N2vr2AtPfgNeWAOff8wQJ0kqGufASZJ0LBu+eXh4O6jn3ax9wzdLU5ckqSIZ4CRJOpaN9x8d3g7qeRc2fqe49UiSKpoBTpKkY+nYd+z+A3uLU4ckSTgHTpKkY6ublM15G8r49xevFknlZ7BFjmrq4LRppa5MY5QBTpKkY7ns89mCJYPdRll9Klz2ueLXJKk8DLXIUbRm3zv3u8iRRpy3UEqSdCxX3A4TZ2RhbaDqU7P2K24vTV2SSm+oRY5SX9bmIkcaBQY4SZKOpbY+2yrgyjugqiZrq5ucHbuFgFTZjrXIUepzkSONCm+hlCTpeGrr4Te/AjuezI4X/ai09UgqDy5ypBJwBE6SJEk6GXWTjt3vIkcaBQY4SZIk6WRc9vmj58ceFFUucqRRYYCTJEmSTsZQixxFVdbmIkcaBQY4SZIk6WQMtcjRaVPh1/6pixxpVBjgJEmSpJN1cJGjsz8Cjb8Bd74EpzdC1bhSV6YxygAnSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSZIkKScMcJIkSZKUEwY4SZIkScoJA5wkSZIk5YQBTpIkSZJywgAnSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSZIkKSeGFeAioikitkXEixFx1yD9tRGxqr//HyJi+oC+L/e3b4uIa0ewdkmSJEmqKMcNcBExDrgX+AQwC/hMRMw64rTPAa0ppfOAZcA9/c+dBdwEXAQ0Acv7X0+SJEmSdIIipXTsEyIuB76WUrq2//jLACmlpQPOeaT/nKciohp4EzgDuGvguQPPG+p61dXVqb6+vqA3NdI6ejoAqKuuq8jrj5RyfB/lWNNB5VBbOdQwUsbSexmLRu3vp6s9+37KhOG1n+zrqXB5/bPNa90jxf+WDhn4nor1/kbrOmPx72cwZfw+f/WrXz2bUpo7WF/1MJ4/FXh1wHEL8JGhzkkp9UTEr4D397c/fcRzpx55gYi4Dbit//EwSiquUv/AV+rrj5RyfB/lWNNB5VBbOdQwUsbSexmLRu3vZ6h/lE/2H+sy/Ed+zMjrn21e6x4p/rd0yMD3VKz3N1rXGYt/P4PJ6fscToAbdSml+4D7AObOnZuam5tLXJEkSZIklcaxBrWGs4jJa8DZA46n9bcNek7/LZTvA/YO87mSJEmSpGEYToDbCMyMiBkRcQrZoiRrjjhnDfD7/Y9vBNanbHLdGuCm/lUqZwAzgWdGpnRJkiRJqizHvYWyf07bF4FHgHHAd1NKWyLibqA5pbQG+A7wdxHxIrCPLOTRf96DwAtAD/CFlFLvKL0XSZIkSRrTjrsKZbE5B06SJElSJYuIIVehHNZG3pIkSZKk0jPASZIkSVJOGOAkSZIkKScMcJIkSZKUEwY4SZIkScqJsluFMiJ2AztLXccImQzsKXUR0ijws62xys+2xio/2xqrxupnuzGldMZgHWUX4MaSiGgeavlPKc/8bGus8rOtscrPtsaqSvxsewulJEmSJOWEAU6SJEmScsIAN7ruK3UB0ijxs62xys+2xio/2xqrKu6z7Rw4SZIkScoJR+AkSZIkKScMcJIkSZKUEwa4URYRX4+I/xURP4+Iv4+I00tdk1SIiGiKiG0R8WJE3FXqeqSREBFnR8RPI+KFiNgSEXeUuiZpJEXEuIjYFBH/T6lrkUZKRJweEQ/1/6y9NSIuL3VNxWCAG32PArNTSr8O/AL4conrkU5aRIwD7gU+AcwCPhMRs0pblTQieoAvpZRmAR8FvuBnW2PMHcDWUhchjbBvAOtSSh8C/ikV8hk3wI2ylNKPU0o9/YdPA9NKWY9UoA8DL6aUtqeUuoDvA9eXuCapYCmlN1JKz/U/biP7IWBqaauSRkZETAN+G7i/1LVIIyUi3gdcBXwHIKXUlVJ6u6RFFYkBrrj+BbC21EVIBZgKvDrguAV/yNUYExHTgTnAP5S4FGmk/DVwJ9BX4jqkkTQD2A2s6L89+P6ImFDqoorBADcCIuKxiNg8yNf1A875E7JbdFaWrlJJ0rFERD3wA+DfpJTeKXU9UqEi4neAXSmlZ0tdizTCqoFLgP+cUpoDtAMVMTe/utQFjAUppXnH6o+IPwB+B7gmufGe8u014OwBx9P626Tci4gasvC2MqW0utT1SCPkSuCTETEfOBU4LSK+l1K6tcR1SYVqAVpSSgfvlniICglwjsCNsohoIrtt4ZMppY5S1yMVaCMwMyJmRMQpwE3AmhLXJBUsIoJsHsXWlNJflboeaaSklL6cUpqWUppO9v/s9YY3jQUppTeBVyPigv6ma4AXSlhS0TgCN/r+BqgFHs1+PuDplNIflrYk6eSklHoi4ovAI8A44LsppS0lLksaCVcCvwf8Y0Q839/2lZTSw6UrSZJ0HP8aWNn/S+XtwKIS11MU4R19kiRJkpQP3kIpSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSZIkKScMcJIkSZKUEwY4SZIkScoJA5wkSZIk5YQBTpIkSZJywgAnSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSZIkKScMcJIkSZKUEwY4SZIkScoJA5wkSZIk5YQBTpIkSZJywgAnSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSRUjIm6OiOaI2B8Rb0TE2oj4jYj4WkR0R0Rb/9cvIuJvIuKsAc+9OiL6+p978Ot/lPL9SJIqjwFOklQRIuKPgL8G/gyYApwDLAeu7z9lVUqpAZgEfAr4NeDZgSEOeD2lVD/g67qivQFJkjDASZIqQES8D7gb+EJKaXVKqT2l1J1S+h8ppX838Nz+9i3AQmA38KUSlCxJ0qAMcJKkSnA5cCrw98N9QkqpF/i/gY+NVlGSJJ0oA5wkqRK8H9iTUuo5wee9TnZL5UEfiIi3B3z97siVKEnS8VWXugBJkopgLzA5IqpPMMRNBfYNOH49pTRtZEuTJGn4HIGTJFWCp4BO4IbhPiEiqoDrgCdGqSZJkk6YI3CSpDEvpfSriPgqcG9E9AA/BrqBecBvAh0Hz42IamAm8DWylSj/qugFS5I0BEfgJEkVIaX0l8AfAf+ebHXJV4EvAj/sP2VhROwHfgWsIbvt8tKU0uvFr1aSpMFFSqnUNUiSJEmShsEROEmSJEnKiYICXER8NyJ2RcTmIfojIr4ZES9GxM8j4pJCridJkiRJlazQEbj/C2g6Rv8nyCaCzwRuA/5zgdeTJEmSpIpVUIBLKf1PDt8f50jXA/81ZZ4GTo+Iswq5piRJkiRVqtHeRmAq2SpfB7X0t70x8KSIuI1shI4JEyZc+qEPfWiUy5IkldpbHW+x58CeIfsbaho4rfa0E37d6qpq6mvqCylNkjRWvfM67H9r6P5T3wennl60coby7NaX96SUzhisryz2gUsp3QfcBzB37tzU3Nxc4ookSaPtB7/4AfdsvIcDPQeO6htfPZ67PnwXC2YuKEFlkqQx69m/hXV3QXfH0X01dfCJe+CSzxa/riNExM6h+kZ7FcrXgLMHHE/rb5MkVbimGU1UDfHPUBVVNE0/1hRrSZJOwuwFEENEoKiCi8r/F4ejHeDWAJ/tX43yo8CvUkpvHO9JkqSxb0LNBJbPW86E6gnUVNUA2cjbhOqsva6mrsQVSpLGnNoGuOUhOKUexp2StdXUZce3PAS15X8LfkEbeUfEfwOuBiYDbwF/CtQApJT+S0QE8DdkK1V2AItSSse8P9JbKCWpsnR0d7D0maXsPbCXeY3zaJreZHiTJI2uzv2w9k7YvxtmXZeNvJVReIuIZ1NKcwfrK2gOXErpM8fpT8AXCrmGJGlsq6up49IplwJww3k3lLYYSVJlqK2Hxiuzx3NuKW0tJ2i0b6GUJEmSJI0QA5wkSZIk5YQBTpIkSZJywgAnSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSZIkKScMcJIkSZKUEwY4SZIkScoJA5wkSZIk5YQBTpIkSZJywgAnSZIkSTlhgJMkSZKknDDASZIkSVJOGOAkSZIkKScMcJIkSZKUEwY4SZIkScqJ6lIXIEmSJKlMdLbB5tWwbztMOhdmL4DahlJXNfI622Dnk9C+B/p6cvU+DXCSJEmSYOdTsPJGSH3Q3QE1dfDIV+CWh6Dx8lJXN3J2PgWrboXTpsLEc+CXP4af3A0Lv5eL92mAkyRJkipdZ1sW3rr2H2rr7si+r7wRvrQNautLU9tI6mzLwtunvw3vvJG1zbkFXlqftd/xs7J/nwY4SZIkqdJtXp2NvA2mtwvW3gmNVxa3ptGw88ls5O2dN2D/W1A/JWv/4MfhnI/CltVwyWdLW+NxuIiJJEmSVOn2bT804nak3i7Yv7u49YyW9j3ZbZOQhbcpsw71nXUxtO4oRVUnxBE4SZIkqdJNOjeb8zZYiKupg1nXZbca5l1fTzbnbbD38sbzcP61RS/pRDkCJ0mSJFW62QsghogGUQUXLShuPaNl9gJ45elszttAL63P2nPwPh2BkyRJkipdbUO22uTKG7NbJnu7spG3qMray3xhj2GrbchWm1x1azbn7ayLs5G3V57O2nPwPg1wkiRJkrIl9L+0LVuwZP/u7LbJixbkItSckMbLs9Umt6zO5rydfy186lu5eZ8GOEmSJEmZ2vpDq02OhTlvQ6mtL/vVJofiHDhJkiRJygkDnCRJkiTlhAFOkiRJknLCACdJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJkiRJygkDnCRJkiTlhAFOkiRJknLCACdJkiRJOWGAkyRJkqScKCjARURTRGyLiBcj4q5B+s+JiJ9GxKaI+HlEzC/kepIkSZJUyU46wEXEOOBe4BPALOAzETHriNP+PfBgSmkOcBOw/GSvJ0mSJEmVrpARuA8DL6aUtqeUuoDvA9cfcU4CTut//D7g9QKuJ0mSJEkVrZAANxV4dcBxS3/bQF8Dbo2IFuBh4F8P9kIRcVtENEdE8+7duwsoSZIkSZLGrtFexOQzwP+VUpoGzAf+LiKOumZK6b6U0tyU0twzzjhjlEuSJEmSpHwqJMC9Bpw94Hhaf9tAnwMeBEgpPQWcCkwu4JqSJEmSVLEKCXAbgZkRMSMiTiFbpGTNEee8AlwDEBEXkgU475GUJEmSpJNw0gEupdQDfBF4BNhKttrkloi4OyI+2X/al4B/GRE/A/4b8AcppVRo0ZIkSZJUiaoLeXJK6WGyxUkGtn11wOMXgCsLuYYkSZKkI3S2webVsG87TDoXZi+A2oaRed2dT0L7HujrGbnX1YgpKMBJkiRJKrKdT8HKGyH1QXcH1NTBI1+BWx6CxssLe91Vt8JpU2HiOfDLH8NP7oaF3yvsdTWiDHCSJElSXnS2ZeGta/+htu6O7PvKG+FL26C2/uRed9Wt8OlvwztvZG1zboGX1mftd/zs5F5XI84AJ0mSJOXF5tXZyNtgertg7Z3QeBIzmHY+mY28vfMG7H8L6qdk7R/8OJzzUdiyGi757MnXrREz2vvASZIkSRop+7YfGnE7Um8X7D/JBd/b92S3TUIW3qbMOtR31sXQuuPkXlcjzhE4SZIkKS8mnZvNeRssxNXUwazrslsfT1RfTzbnbbDnvvE8nH/tib+mRoUjcJIkSVJezF4AMcSP8FEFFy04+dd95elszttAL63P2k/2dTXiHIGTJEmS8qK2IVttcuWN2S2TvV3ZyFtUZe0nu9BIbUO22uSqW7M5b2ddnI28vfJ01u4CJmXDACdJkiTlSePl2WqTa+/M5rzNui4bISs0ZDVenq02uWV1Nuft/GvhU98yvJUZA5wkSZKUN7X1h1abPJk5b8d6XVebLGsGOEljXnt3O+teXsfOtp00NjTSNKOJCTUTSl2W+rV3t9P8ZjOtna309vX69yNJ0jG4iImkMe25t57jmgev4Z6N97Bi8wru2XgP1zx4Dc+99VypSxPZ38/81fP5ResvOKXqFJ547Qnmr57v348kSUNwBE7SmNXe3c7ixxbT3tP+XtuBngMALH5sMet/dz11NXWlKq/itXe3s+TxJSz92FJ2dewC4IbzbmDD6xtY8vgS1i5Y69+PJElHMMBJGrPWvbyOPvoG7evq62LpM0u5dMqlRa5KBzW/2cyUuins6tjFngN7mDx+MgBXfOAK5pw5h3U71rFgpstWS5I0kLdQShqzdrbtfG/E7Ujdfd3sPbC3yBVpoNbOVqbWTwVg8vjJzJw4872+CyddSEtbS6lKkySpbDkCJ2nMamxoZHz1+EFD3Pjq8cxrnMcN591Q/MIEQG9fL0+89sSgfwdb923lqmlXFb8oSZLKnCNwksasphlNVA3xv7kqqmia3lTkijRQ04wmNu3axIbXNxzWvuH1DWzatcm/H0mSBuEInKQxa0LNBJbPW87ixxbT1ddFd18346vHU0UVy+ctd4GMEptQM4FlVy9jyeNLmHPmHC6cdCFb921l065NLLt6mX8/kiQNIlJKpa7hMHPnzk3Nzc2lLkPSGNLR3cHSZ5ay98Be5jXOo2l6k+GgjHR0d7Buxzpa2lqY1jDNvx9J+dbZBptXw77tMOlcmL0AahtG5zpr/xja98CHfnv0rqOSiIhnU0pzB+0zwEmqBD988YcAznmTJI2enU/Byhsh9UF3B9TUQVTBLQ9B4+Uje51Vt8JpU2HiOZASvPI0LPzeyF5HJXOsAOctlJIkSVKhOtuy8Na1/1Bbd0f2feWN8KVtUFs/MtdZdSt8+tvwzhtZ25xb4KX1WfsdPxuZ66hsGeAkSZKkQm1enY28Daa3C9beCY1XFn6dnU9mI2/vvAH734L6KVn7Bz8O53wUtqyGSz5b+HVUtlyFUpIkSSrUvu2HRtyO1NsF+3ePzHXa92S3TUIW3qbMOtR31sXQumNkrqOy5QicJEmSVKhJ52Zz3gYLcTV1MOu67FbHQvX1wC9/PPhrvfE8nH9t4ddQWXMETpIkSSrU7AXZgiWDiSq4aMHIXeeVp7M5bwO9tD5rH6nrqGw5AidJkiQVqrYhW21y5Y3ZLZO9XYevQjlSC4vUNmSrTa66NZvzdtbF2cjbwVUoXcBkzDPASZIkSSOh8fJstcm1d2Zz3mZdl42IjXSoarw8W21yy+psztv518KnvmV4qxAGOEmSJGmk1NYfWm1yJOa8Hes6rjZZkZwDJ0mSJEk5YYCTJEmSpJwwwEmSJElSThjgJEmSJCknDHCSJEmSlBMGOEmSJEnKCQOcJEmSJOWEAU6SJEmScsKNvCUNqr27nXUvr2Nn204aGxppmtHEhJoJpS7rpLR3t9P8ZjOtna309vXm+r1IkqTK5gicpKM899ZzXPPgNdyz8R5WbF7BPRvv4ZoHr+G5t54rdWkn7Lm3nmP+6vn8ovUXnFJ1Ck+89gTzV8/P5XuRJElyBE7SYdq721n82GLae9rfazvQcwCAxY8tZv3vrqeupq5U5Z2Q9u52ljy+hKUfW8qujl0A3HDeDWx4fQNLHl/C2gVrc/NeJEllrLMNNq+GfdvhtA9ASzO8+yvo64HZC6C2odQVagwxwEk6zLqX19FH36B9XX1dLH1mKZdOubTIVZ2c5jebmVI3hV0du9hzYA+Tx08G4IoPXMGcM+ewbsc6FsxcUOIqJUm5tvMpWHkjpD7o7sjaIuCcy+GXP4af3A0LvweNl5e2To0Z3kIp6TA723a+N+J2pO6+bvYe2Fvkik5ea2crU+unAjB5/GRmTpz5Xt+Fky6kpa2lVKVJksaCzrYsvHXtPxTeAFKClmfhU9+CT38bVt0KnftLV6fGFEfgJB2msaGR8dXjBw1x46vHM69xHjecd0PxCzsJvX29PPHaE4PWu3XfVq6adlXxi5IkjR2bV2cjb4NKsGU1XPJZOOejhx5LBSpoBC4imiJiW0S8GBF3DXHO70bECxGxJSIeKOR6kkZf04wmqob4X0MVVTRNbypyRSevaUYTm3ZtYsPrGw5r3/D6Bjbt2pSr9yJJKkP7th8+8jZQbxfs3Z49PutiaN1RrKo0xp30CFxEjAPuBf450AJsjIg1KaUXBpwzE/gycGVKqTUiziy0YEmja0LNBJbPW87ixxbT1ddFd18346vHU0UVy+ctz9WiHxNqJrDs6mUseXwJc86cw4WTLmTrvq1s2rWJZVcvy9V7kSSVoUnnQk3d4CGupg7ef272+I3n4fxri1qaxq5IKZ3cEyMuB76WUrq2//jLACmlpQPO+QvgFyml+4f7unPnzk3Nzc0nVZOkkdPR3cHSZ5ay98Be5jXOo2l6U24DT0d3B+t2rKOlrYVpDdNy/V4kKfcGrtg46dx8r9LY2QZ/+aFsDtyRTqmHL22DlmfgB/8S7vgZ1NYXv0blUkQ8m1KaO1hfIXPgpgKvDjhuAT5yxDnn9xfwJDCOLPCtG6TA24DbAM4555wCSpI0Uupq6t5bbTIvc96GUldT52qTklQOjlyxsaYOHvkK3PJQPldprG2AeV+Dh/9df0P/wEiMg4sWwN//K3jl6WwVSsObRshoL2JSDcwErgamAf8zIv5JSuntgSellO4D7oNsBG6Ua5IkSVKxDVyx8aCDtx6uvDEbrcpbyOlsg8f/HG5aCf/rR7B/N5z/W7B7Gzz/APzz/5CtRJm396WyVkiAew04e8DxtP62gVqAf0gpdQMvR8QvyALdxgKuK4157d3trHt5HTvbdtLY0EjTjCYm1Ewoeg3NbzbT2tlKb19vSWqQJI0hx1qxsbcL1t4JjVcWt6ZC7XwSTpsKB96G95+X1T/nlqyv7Q2orjW8acQVsgrlRmBmRMyIiFOAm4A1R5zzQ7LRNyJiMtktldsLuKY05j331nNc8+A13LPxHlZsXsE9G+/hmgev4bm3nitqDfNXz+cXrb/glKpTeOK1J5i/en5Ra5AkjTHHW7Fx/+7i1jMS2vfAxP7pP/VTYMqsQ32uPKlRctIjcCmlnoj4IvAI2fy276aUtkTE3UBzSmlNf99vRcQLQC/w71JK+dkFWCqy9u52Fj+2mPae9vfaDu7Htvixxaz/3fWjvvhGe3c7Sx5fwtKPLWVXxy4gmwO34fUNLHl8CWsXrHUBEEnSiTveio2zrjs0epUXfT3wyx8PXrcrT2qUFDQHLqX0MPDwEW1fHfA4AX/U/yXpONa9vI4+Br+9pKuvi6XPLH1vYZHR0vxmM1PqprCrYxd7Duxh8vjJAFzxgSuYc+Yc1u1Y54IgkqQTN3tBtmDJYKIqW/Qjb2YvgJ/cDS+thw9+/FD7S+uzxUs+9a3S1aYxa7QXMZF0Ana27XxvxO1I3X3d7D0w+gPYrZ2tTK2fCsDk8ZOZOXHme30XTrqQlraWUa9BkjQG1TZkq02uvDG7ZbK3Kxt5i6qsPY9zxWobshUmV90K53w0u23yjeddeVKjygAnlZHGhkbGV48fNMSNrx7PvMZ5o76kf29fL0+89sSg19m6bytXTbtqVK8vSRrDGi/PVptce2c2523WddnIW56DTuPl2R5vW1Znc97Ov9aVJzWqClnERNIIa5rRRNUQ/1lWUUXT9Kai1LBp1yY2vL7hsPYNr29g065NRalBkjSG1dZnqzVedANc8tmxEXRq67P3cs1Xx857UtlyBE4qIxNqJrB83nIWP7aYrr4uuvu6GV89niqqWD5veVEWD5lQM4FlVy9jyeNLmHPmHC6cdCFb921l065NLLt6mQuYSJIklVBk64yUj7lz56bm5uZSlyGVVEd3B0ufWcreA3uZ1ziPpulNRQ9OHd0drNuxjpa2FqY1TCtJDZKkMWrTyux73ladlIokIp5NKc0drM8ROKkM1dXUvbfa5GjPeTtWDa42KUmSVF6cAydJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJkiRJygkDnCRJkiTlhAFOkiRJknLCACdJkiRJOVFd6gIkSZKUQ51tsHk17NsOk86F2QugtmF4z9v5JLTvgb6e4T9PEmCAkyRJ0ona+RSsvBFSH3R3QE0dPPIVuOUhaLz82M9bdSucNhUmngO//DH85G5Y+L1jP0/SewxwkiRJGr7Otiy8de0/1NbdkX1feSN8aRvU1g/+vFW3wqe/De+8kbXNuQVeWp+13/GzwZ8n6TAGOEmSJA3f5tXZyNtgertg7Z3QeOXRfTufzEbe3nkD9r8F9VOy9g9+HM75KGxZDZd8dvTqlsYIFzGRJEnS8O3bfmjE7Ui9XbB/9+B97Xuy2yYhC29TZh3qO+tiaN0xklVKY5YjcJIkSRq+Sedmc94GC3E1dTDruuzWyCP19WRz3gbre+N5OP/aES9VGoscgZMkSdLwzV4AMcSPkFEFFy0Y+nmvPJ3NeRvopfVZ+1DPk3QYR+AkSZI0fLUN2WqTK2/Mbpns7cpG3qIqax9qIZLahmy1yVW3ZnPezro4G3l75ems3QVMpGExwEmSJOnENF6erTa59s5sztus67IRtOOFsMbLs9Umt6zO5rydfy186luGN+kEGOAkSZJ04mrrD602Odi8tmM9z9UmpZPmHDhJkiRJygkDnCRJkiTlhAFOkiRJknLCACdJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJkiRJygkDnCRJkiTlhAFOkiRJknLCACdJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwoKcBHRFBHbIuLFiLjrGOd9OiJSRMwt5HqSJEmSVMlOOsBFxDjgXuATwCzgMxExa5DzGoA7gH842WtJkiRJkgobgfsw8GJKaXtKqQv4PnD9IOf9R+Ae4N0CriVJkiRJFa+QADcVeHXAcUt/23si4hLg7JTSj471QhFxW0Q0R0Tz7t27CyhJkiRJksauUVvEJCKqgL8CvnS8c1NK96WU5qaU5p5xxhmjVZIkSZIk5Vp1Ac99DTh7wPG0/raDGoDZwOMRAfBrwJqI+GRKqbmA60ojor27nXUvr2Nn204aGxppmtHEhJoJpS4LyGprfrOZ1s5Wevt6y6o2SZIklU4hI3AbgZkRMSMiTgFuAtYc7Ewp/SqlNDmlND2lNB14GjC8qSw899ZzXPPgNdyz8R5WbF7BPRvv4ZoHr+G5t54rdWk899ZzzF89n1+0/oJTqk7hideeYP7q+WVRmyRJkkrrpEfgUko9EfFF4BFgHPDdlNKWiLgbaE4prTn2K0il0d7dzuLHFtPe0/5e24GeAwAsfmwx6393PXU1dSWrbcnjS1j6saXs6tgFwA3n3cCG1zew5PElrF2wtmS1SZIkqfQKuYWSlNLDwMNHtH11iHOvLuRa0khZ9/I6+ugbtK+rr4ulzyzl0imXFrmqTPObzUypm8Kujl3sObCHyeMnA3DFB65gzplzWLdjHQtmLihJbZIkSSq9UVvERCpXO9t2vjfidqTuvm72Hthb5IoOae1sZWp9tpjr5PGTmTlx5nt9F066kJa2llKVJkmSpDJQ0AiclEeNDY2Mrx4/aIgbXz2eeY3zuOG8G4pfGNDb18sTrz0x6PW37tvKVdOuKn5RkiRJKhuOwKniNM1oomqIj34VVTRNbypyRYc0zWhi065NbHh9w2HtG17fwKZdm0pamyRJkkrPEThVnAk1E1g+bzmLH1tMV18X3X3djK8eTxVVLJ+3vKSLhEyomcCyq5ex5PElzDlzDhdOupCt+7ayadcmll29zAVMJEmSKlyklEpdw2Hmzp2bmpvdaUCjr6O7g6XPLGXvgb3Ma5xH0/SmsglIHd0drNuxjpa2FqY1TCur2iRJ/TrbYPNq2LcdJp0LsxdAbUOpqyqezjZY+8fQvgc+9NuV9/6lURQRz6aU5g7aZ4BTJfvhiz8EKNmcN0lSTu18ClbeCKkPujugpg6iCm55CBovL3V1o2/nU7DqVjhtKkw8B1KCV56Ghd+rjPcvjbJjBThvoZQkSToRnW1ZeOvaf6ituyP7vvJG+NI2qK0vTW3F0NmWhbdPfxveeSNrm3MLvLQ+a7/jZ2P7/UslZoCTJEk6EZtXZyNvg+ntgrV3QuOVxa2pmHY+mY28vfMG7H8L6qdk7R/8OJzzUdiyGi75bGlrlMYwA5yKor27nXUvr2Nn204aGxppmtHEhJoJJa+p+c1mWjtb6e3rLYuaJEk5sG/7oRG3I/V2wf7dxa2n2Nr3ZLdNQhbepsw61HfWxdC6oxRVSRXDAKdR99xbz7H4scX00ceBngOMrx7P1zd+neXzlnPJlEtKVtOSx5cwpW4KU+un8sRrT/DNTd9k2dXLSlaTJCknJp2bzXkbLMTV1MGs67JbCseqvh745Y8Hf49vPA/nX1v0kqRK4j5wGlXt3e0sfmwx7T3t722cfaDnAO09WXvHUL/BHOWaljy+hKUfW8rNF97MPzv7n/HXv/nXLP3YUpY8vqQkNUmScmT2gmzBksFEFVy0oLj1FNvsBdmCJS+tP7z9pfVZ+1h//1KJOQKnUbXu5XX0Mfg8ga6+LpY+s5RLp1xa1Jqa32xmSt0UdnXsYs+BPUwePxmAKz5wBXPOnMO6HetYMNN/fCRJQ6htyFabXHljdstkb9fhq1CO9QU8ahuy1SZX3ZrNeTvr4mzk7eAqlGP9/UslZoDTqNrZtvO9kbcjdfd1s/fA3iJXBK2drUytnwrA5PGTmTlx5nt9F066kJa2lqLXJEnKmcbLs9Um196ZzXmbdV028lQp4aXx8my1yS2rszlv518Ln/pW5bx/qYQMcBpVjQ2NjK8eP2iIG189nnmN84q+B1tvXy9PvPbEoNfdum8rV027qqj1SJJyqrb+0GqTY3nO21Bq611tUioB58BpVDXNaKJqiI9ZFVU0TW8qckVZTZt2bWLD6xsOa9/w+gY27dpUkpokSZKk4XAETqNqQs0Els9bzuLHFtPV10V3Xzfjq8dTRRXL5y2nrqauJDUtu3oZSx5fwpwz53DhpAvZum8rm3ZtYtnVy0pSkyRJkjQckVIqdQ2HmTt3bmpubi51GRphHd0dLH1mKXsP7GVe4zyapjeVPCh1dHewbsc6WtpamNYwrSxqkiTlzKaV2fdKvIVS0qiJiGdTSnMH63METkVRV1P33mqTxZ7zNpS6mjpXm5QkSVKuOAdOkiRJknLCACdJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJkiRJygkDnCRJkiTlhAFOkiRJknLCACdJkiRJOVFd6gIkSVKF62yDzath33aYdC7MXgC1DaWu6vg622Dnk9C+B/p68lO3pFwzwEmSpNLZ+RSsvBFSH3R3QE0dPPIVuOUhaLy81NUNbedTsOpWOG0qTDwHfvlj+MndsPB75V23pNwzwEmSpNLobMvCW9f+Q23dHdn3lTfCl7ZBbX1pajuWzrYsvH362/DOG1nbnFvgpfVZ+x0/K8+6JY0JBjhJklQam1dnI2+D6e2CtXdC45XFrWk4dj6Zjby98wbsfwvqp2TtH/w4nPNR2LIaLvlsaWuUNGa5iIkkSSqNfdsPjbgdqbcL9u8ubj3D1b4nu20SsvA2ZdahvrMuhtYdpahKUoVwBE6SJJXGpHOzOW+DhbiaOph1XXZrYrnp68nmvA1W2xvPw/nXFr0kSZXDEThJklQasxdADPGjSFTBRQuKW89wzV4ArzydzXkb6KX1WXu51i1pTHAETpIklUZtQ7ba5Mobs1sme7uykbeoytrLdSGQ2oZstclVt2Zz3s66OBt5e+XprL1c65Y0JhjgJElS6TRenq02ufbObM7brOuyEaxyD0GNl2erTW5Znc15O/9a+NS3yr9uSblngJMkSaVVW39otclynPM2lNp6V5uUVHTOgZMkSZKknDDASZIkSVJOFBTgIqIpIrZFxIsRcdcg/X8UES9ExM8j4icR0VjI9SRJkiSpkp10gIuIccC9wCeAWcBnImLWEadtAuamlH4deAj4i5O9niRJkiRVukJG4D4MvJhS2p5S6gK+D1w/8ISU0k9TSgd353wamFbA9SRJkiSpohUS4KYCrw44bulvG8rngLWDdUTEbRHRHBHNu3fvLqAkSZIkSRq7irKISUTcCswFvj5Yf0rpvpTS3JTS3DPOOKMYJUmSJElS7hSyD9xrwNkDjqf1tx0mIuYBfwL8s5RSZwHXkyRJkqSKVsgI3EZgZkTMiIhTgJuANQNPiIg5wLeAT6aUdhVwLUmSJEmqeCcd4FJKPcAXgUeArcCDKaUtEXF3RHyy/7SvA/XAf4+I5yNizRAvJ0mSJEk6jkJuoSSl9DDw8BFtXx3weF4hry9o725n3cvr2Nm2k8aGRppmNDGhZkKpyzph7d3tNL/ZTGtnK719vbl9H5IqWGcbbF4N+7bDpHNh9gKobSh1VWNDZxvsfBLa90Bfj3+2knQMkVIqdQ2HmTt3bmpubi51GWXhubeeY/Fji+mjjwM9BxhfPZ4qqlg+bzmXTLmk1OUN23NvPceSx5cwpW4KU+unkkhs2rWJZVcvy9X7kFTBdj4FK2+E1AfdHVBTB1EFtzwEjZeXurp82/kUrLoVTpsKE8+BlOCVp2Hh9/yzlVSxIuLZlNLcQfsMcOWpvbudax68hvae9qP6JlRPYP3vrqeupq4ElZ2Y9u525q+ez9KPLWVXRzYN8obzbmDD6xv48hNfZu2Ctbl4H5IqWGcb/OWHoGv/0X2n1MOXtkFtffHrGgs62+AbF8Onvw3vvJG1zbkFXloPP/iXcMfP/LOVVJGOFeAKuoVSo2fdy+voo2/Qvq6+LpY+s5RLp1xa5KpOXPObzUypm8Kujl3sObCHyeMnA3DFB65gzplzWLdjHQtmLihxlZJ0DJtXZyNvg+ntgrV3QuOVxa1prNj5ZDby9s4bsP8tqJ+StX/w43DOR2HLarjks6WtUZLKTFH2gdOJ29m2kwM9Bwbt6+7rZu+BvUWu6OS0drYytT7b333y+MnMnDjzvb4LJ11IS1tLqUqTpOHZtz27bXIwvV2wf3dx6xlL2vdkt01CFt6mzDrUd9bF0LqjFFVJUllzBK5MNTY0Mr56/KAhbnz1eOY1zuOG824ofmEnqLevlydee2LQWrfu28pV064qflGSdCImnZvNeRssxNXUwazrstv+dOL6euCXPx78z++N5+H8a4tekiSVO0fgylTTjCaqhvjrqaKKpulNRa7o5DTNaGLTrk1seH3DYe0bXt/Apl2bcvM+JFWw2QuyBUsGE1VwkbeBn7TZC7IFS15af3j7S+uzdv9sJekojsCVqQk1E1g+bzmLH1tMV18X3X3dh61CmZeFPybUTGDZ1ctY8vgS5pw5hwsnXcjWfVvfW4UyL+9DUgWrbchWm1x5Y3bLZG/X4atQusjGyattyFabXHVrNuftrIuzkbeDq1D6ZytJR3EVyjLX0d3B0meWsvfAXuY1zqNpelMuQ09Hdwfrdqyjpa2FaQ3Tcvs+JI0BJ7ufW9ub8OBnYf+uLGzM+xo0/Nqol1sROvdnC5a07oCJ07ORN8ObpArmNgI598MXfwiQizlvklTWTnY/N/cqkyQVkdsISJLU2ZaFt4H7uR1cmGTljUPv59bZloW3wfYqW3Wre5VJkorKACdJqgwnu5+be5VJksqIq1BKkirDye7n5l5lkqQy4gicJKkynOx+bu5VJkkqI47ASZIqw8nu5+ZeZZKkMuIInCSpMpzsfm7uVSZJKiMGOElS5Wi8PFttcu2d2Zy3WdcNb8+xxsuz1SYP7lV2/rXwqW8Z3iRJRWeAkyRVltr6Q6tNDjav7VjPc7VJSVKJOQdOkiRJknLCACdJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJkiRJygkDnCRJkiTlhBt5H0d7dzvrXl7HzradNDY00jSjiQk1E4p6/eY3m2ntbKW3r7fo15ek3Opsg82rYd92mHQuzF4AtQ1Z+84noX0P9PUcapckKQcipVTqGg4zd+7c1NzcXOoyAHjuredY/Nhi+ujjQM8BxlePp4oqls9bziVTLinK9Zc8voQpdVOYWj+VRGLTrk0su3pZUa4vSbm18ylYeSOkPujugJo6iCqY9zV4/M/htKkw8RxICV55GhZ+DxovL3XVkiQBEBHPppTmDtpngBtce3c71zx4De097Uf1TaiewPrfXU9dTd2oXn/+6vks/dhSdnXsAuCG825gw+sb+PITX2btgrWjen1Jyq3ONvjLD0HX/kE6A25aCQfezg7n3AIvrYcf/Eu442dQW1/MSiVJGtSxApy3UA5h3cvr6KNv0L6uvi6WPrOUS6dcOmrXb36zmSl1U9jVsYs9B/YwefxkAK74wBXMOXMO63asY8HMBaN2fUnKrc2rs5G3ofyvH8H7z4P6KdnxBz8O53wUtqyGSz5bnBolSTpJLmIyhJ1tOznQc2DQvu6+bvYe2Duq12/tbGVq/VQAJo+fzMyJM9/ru3DShbS0tYzq9SUpt/Ztz26bHFSC/buz8DZl1qHmsy6G1h1FKE6SpMI4AjeExoZGxlePHzTEja8ez7zGedxw3g2jdv3evl6eeO2JQa+xdd9Wrpp21ahdW5JybdK52Zy3wUJcjINZ12W3Tg70xvNw/rVFKU+SpEI4AjeEphlNVA3xx1NFFU3Tm0b9+pt2bWLD6xsOa9/w+gY27do06teXpNyavSBbsGQwqQ/q3n9420vrs4VMLvK2dElS+XMEbggTaiawfN5yFj+2mK6+Lrr7ug9bhXK0FxCZUDOBZVcvY8njS5hz5hwunHQhW/dtfW8VShcwkaQh1DbALQ9lq1D2dmVfA1ehXHM7PP9AdtvkG88fWoXSBUwkSTngKpTH0dHdwdJnlrL3wF7mNc6jaXpTUcNTR3cH63aso6WthWkN04p+fUnKrc79sPbObM7brOuyEbba+qx9y+psztvE6YfaJUkqE65CWYC6mrr3VpsczTlvx7q+q01K0kmorYfGK7PHA+e81da72qQkKbecAydJkiRJOWGAkyRJkqScMMBJkiRJUk4Y4CRJkiQpJwxwkiRJkpQTBjhJkiRJyomCAlxENEXEtoh4MSLuGqS/NiJW9ff/Q0RML+R6kiRJklTJTjrARcQ44F7gE8As4DMRMeuI0z4HtKaUzgOWAfec7PUkSZIkqdIVspH3h4EXU0rbASLi+8D1wAsDzrke+Fr/44eAv4mISCmlIV+1Yx9sWllAWYf7xdu/YH/X/oJeY2JXGw2nNEBb+whVJUkqiv1vQf2UUlchSdKIKSTATQVeHXDcAnxkqHNSSj0R8Svg/cCegSdFxG3Abf2H++OSW7cVUFc5mcwR71UaI/xsK2duHe6JfrY1VvnZ1lg1Vj/bjUN1FBLgRkxK6T7gvlLXMdIiojmlNLfUdUgjzc+2xio/2xqr/GxrrKrEz3Yhi5i8Bpw94Hhaf9ug50RENfA+YG8B15QkSZKkilVIgNsIzIyIGRFxCnATsOaIc9YAv9//+EZg/THnv0mSJEmShnTSt1D2z2n7IvAIMA74bkppS0TcDTSnlNYA3wH+LiJeBPaRhbxKMuZuC5X6+dnWWOVnW2OVn22NVRX32Q4HxCRJkiQpHwrayFuSJEmSVDwGOEmSJEnKCQPcKIuIr0fE/4qIn0fE30fE6aWuSSpERDRFxLaIeDEi7ip1PdJIiIizI+KnEfFCRGyJiDtKXZM0kiJiXERsioj/p9S1SCMlIk6PiIf6f9beGhGXl7qmYjDAjb5HgdkppV8HfgF8ucT1SCctIsYB9wKfAGYBn4mIWaWtShoRPcCXUkqzgI8CX/CzrTHmDmBrqYuQRtg3gHUppQ8B/5QK+Ywb4EZZSunHKaWe/sOnyfbLk/Lqw8CLKaXtKaUu4PvA9SWuSSpYSumNlNJz/Y/byH4ImFraqqSRERHTgN8G7i91LdJIiYj3AVeRrXpPSqkrpfR2SYsqEgNccf0LYG2pi5AKMBV4dcBxC/6QqzEmIqYDc4B/KHEp0kj5a+BOoK/EdUgjaQawG1jRf3vw/RExodRFFYMBbgRExGMRsXmQr+sHnPMnZLforCxdpZKkY4mIeuAHwL9JKb1T6nqkQkXE7wC7UkrPlroWaYRVA5cA/zmlNAdoBypibv5Jb+StQ1JK847VHxF/APwOcE1y4z3l22vA2QOOp/W3SbkXETVk4W1lSml1qeuRRsiVwCcjYj5wKnBaRHwvpXRrieuSCtUCtKSUDt4t8RAVEuAcgRtlEdFEdtvCJ1NKHaWuRyrQRmBmRMyIiFOAm4A1Ja5JKlhEBNk8iq0ppb8qdT3SSEkpfTmlNC2lNJ3s/9nrDW8aC1JKbwKvRsQF/U3XAC+UsKSicQRu9P0NUAs8mv18wNMppT8sbUnSyUkp9UTEF4FHgHHAd1NKW0pcljQSrgR+D/jHiHi+v+0rKaWHS1eSJOk4/jWwsv+XytuBRSWupyjCO/okSZIkKR+8hVKSJEmScsIAJ0mSJEk5YYCTJEmSpJwwwEmSJElSThjgJEmSJCknDHCSJEmSlBMGOEmSJEnKif8fG4p/h+m5fpoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Nice colours\n",
    "c1 = '#2ca02c' # Forest green\n",
    "c2 = '#ff7f0e' # Burnt orange\n",
    "\n",
    "def nice_stem_plot(x, weight, col, alpha=None): \n",
    "    markerline, stemlines, baseline = plt.stem(x, weight, linefmt='-', markerfmt='o', basefmt='k-')\n",
    "    plt.setp(markerline, 'markerfacecolor', col, 'color', col, 'linewidth', 2, markersize=8)\n",
    "    plt.setp(stemlines, 'color', col)\n",
    "    baseline.set_xdata([0,1])\n",
    "    baseline.set_transform(plt.gca().get_yaxis_transform())\n",
    "    #plt.axis('off')\n",
    "\n",
    "def nice_cdf_plot(x, y, col, alpha=None):\n",
    "\n",
    "    plt.step(x, y, where='post', c=col, alpha=0.5)\n",
    "    \n",
    "    plt.scatter(x, y, color=col, s=50, zorder=2)\n",
    "    plt.scatter(x[1:], y[:-1], color='white', s=50, zorder=0, edgecolor=col)\n",
    "    \n",
    "    x_start = plt.xlim()[0]\n",
    "    x_end = plt.xlim()[1]\n",
    "    \n",
    "    plt.step([x_start, x[0]], [0, y[0]], where='post', c=col, alpha=0.5, zorder=0)\n",
    "    plt.step([x[-1], x_end], [y[-1], 1], where='pre', c=col, alpha=0.5, zorder=0)\n",
    "    \n",
    "plt.figure(figsize=(15, 10))\n",
    "\n",
    "ax = plt.subplot(2, 1, 1)\n",
    "ax.set_title('Density plot')\n",
    "nice_stem_plot(X, a, c1)\n",
    "nice_stem_plot(Y, b, c2)\n",
    "xlims = plt.xlim()\n",
    "\n",
    "ax = plt.subplot(2, 1, 2)\n",
    "ax.set_title('CDF')\n",
    "nice_cdf_plot(X, a_cdf, c1)\n",
    "nice_cdf_plot(Y, b_cdf, c2)\n",
    "\n",
    "plt.ylim([0., 1.])\n",
    "plt.xlim(xlims)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construct the monotone rearrangement plan\n",
    "\n",
    "Remember, it is possible, given two smooth densities $f$ and $g$ on $\\RR$, to define an optimal Monge map \n",
    "\n",
    "$$\n",
    "T(x) = G^{-1}\\circ F(x).\n",
    "$$\n",
    "\n",
    "However, we saw earlier in the course that two arbitrary discrete measures usually don't posess an optimal Monge map. Rather, they have a Kantorovic style coupling.\n",
    "\n",
    "The coupling takes a form that is somewhat similar to the map $T$ above. First define the cumulative density vector $A_i = \\sum_{k=0}^i a_k$ and similar for $B_i$. Note that we are indexing from $0$, to make things easier when translating to Python. We can also define a cumulative density _function_ \n",
    "\n",
    "$$\n",
    "A(x) = \\sum_{x_i \\le x} a_i = \\int_{-\\infty}^x \\sum_{i=0}^n a_i \\delta_{x_i} \\, dx\n",
    "$$\n",
    "\n",
    "and a generalised inverse\n",
    "\n",
    "$$\n",
    "A^{-1}(t) = \\inf \\{ x : A(x)  > t \\} = \\min \\{ n : t < \\sum_{k=0}^i a_k \\}.\n",
    "$$\n",
    "\n",
    "Through these generalised inverse it's possible for find a monotone rearrangement coupling that satisfies the notion that $B^{-1}(A(x_i)) = y_i$ and that the marginals hold.\n",
    "\n",
    "First we define a _cumulative_ coupling $Q$ of $P$, \n",
    "\n",
    "$$\n",
    "Q_{i, j} = \\sum_{k, \\ell = 0}^{i, j} P_{i, j}.\n",
    "$$\n",
    "\n",
    "Then as it happens, the cumulative density of the optimal coupling is simply \n",
    "\n",
    "$$\n",
    "Q^\\star_{i, j} = \\min \\{ A_i, B_j \\}\n",
    "$$\n",
    "\n",
    "from which we can find the optimal coupling $P^\\star$. For now I just leave this defined. As an exercise, you could check that this defines a monotone coupling for which the marginals hold.\n",
    "\n",
    "### TASK 5: Find $P^\\star$\n",
    "\n",
    "Implement this method, or something similar, below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAAFlCAYAAAA+gTZIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAb4klEQVR4nO3df7DddX3n8eeLIKCoaJPoKAGTGbJi1K26V7Rr13VltaG1pDMbuuDUQsXJ7kzpqluni7tOdFk7s3Y6YmcWO2YEy4IV2VRnMzWVsqKD7Vg2CbJq+LFGQEjEkgBigwMYeO8f94Q9Xm/uOffe873ne7/n+Zg5k3O+5/vjnTMwvHh/Pt/PN1WFJEmSmnPcuAuQJEnqOgOXJElSwwxckiRJDTNwSZIkNczAJUmS1DADlyRJUsOOb+Kkz3ve82rlypVNnHpkVq1aNe4SJElacvfeey+HDh3KuOsYJMli1626oao2jqSYEWgkcK1cuZKtW7c2ceqRefe73z3uEiRJWnJTU1PjLmFoycJzYVW1qrPikKIkSVLDGulwSZIkLdYiO1wjrGTxDFySJKmVFhO42sbAJUmSWicJxx238JlPTz311AirWTzncEmSJDXMwCVJklopyYJfQ55/Y5K7kuxLcuks3785ya1JjiTZ3Lf9NUm+kWRvkm8l+deDruWQoiRJaqUm53AlWQFcAbwN2A/sSrKjqm7v2+0+4CLgAzMO/wnw21X13SQvBfYkuaGqfnSs6xm4JElSKzU8af4sYF9V3d271nXAJuCZwFVV9/a+e7r/wKr6v33vf5DkQWA18KNjXczAJUmSWqnhwHUqcH/f5/3AG+Z7kiRnAScA35trPwOXJEnqolVJdvd93lZV20Z5gSQvAa4BLqyqp+fa18AlSZJaZz6T34/hUFXN9RyjA8BpfZ/X9LYNJcnzgS8B/6mq/m7Q/gYuSZLUSotZh2sIu4D1SdYxHbTOB945zIFJTgC+CPz3qto+zDEuCyFJklqpyWUhquoIcAlwA3AHcH1V7U1yWZJze9d/fZL9wHnAp5Ls7R3+m8CbgYuS3NZ7vWau69nhkiRJrdT0o32qaiewc8a2rX3vdzE91DjzuGuBa+dzLTtckiRJDRsqcA1aiVWSJGmUFjOc2MaHXg8cUhxyJVZJkqSRamNwWqhh5nANXIlVkiRp1LoUuIYZUpxtJdZTZ+6UZEuS3Ul2Hz58eFT1SZIkLXsju0uxt3rrNoC1a9fWqM4rSZImU8PrcC2pYQLXolZilSRJmq+2Tn5fqGEC14JXYpUkSVqoiQpcVXUkydGVWFcAV1XV3gGHSZIkLcpEBS6YfSVWSZIkDcdH+0iSpFaauA6XJEnSUjNwSZIkNSjJxC0LIUmStOS61OHqTnSUJElqKTtckiSplbrU4TJwSZKkVjJwSZIkNahrj/ZxDpckSVLD7HBJkqRW6lKHy8AlSZJayXW4JEmSGmaHS5IkqUFOmpckSdK8NNLhWrVqFRdddFETpx6ZrVu3jruEgd71rneNu4Q5rV+/ftwlSJI6rEsdLocUJUlSKxm4JEmSGuZdipIkSQ1y0rwkSZLmxQ6XJElqpS51uAxckiSplZzDJUmS1LAudbi6Ex0lSZJayg6XJElqnSQOKUqSJDWtS0OKBi5JktRKBi5JkqQGdW1IsTt/E0mSpJaywyVJklrJIUVJkqSGdWlI0cAlSZJaaaI6XEmuAt4BPFhVr2q+JEmSNOkmcdL8nwEbG65DkiSpswZ2uKrq5iRrl6AWSZKkZ3RpSLE7vTpJktQpSRb8GvL8G5PclWRfkktn+f7NSW5NciTJ5hnfXZjku73XhYOuNbJJ80m2AFsATj/99FGdVpIkTaCm53AlWQFcAbwN2A/sSrKjqm7v2+0+4CLgAzOO/QXgw8AUUMCe3rGPHOt6I/ubVNW2qpqqqqnVq1eP6rSSJElNOAvYV1V3V9WTwHXApv4dqureqvoW8PSMY38FuLGqHu6FrBsZMN/dZSEkSVIrLXIO16oku/s+b6uqbX2fTwXu7/u8H3jDkOee7dhT5zpgmGUhPge8henC9wMfrqorhyxIkiRpQRY5pHioqqZGVctiDXOX4gVLUYgkSVK/hu9SPACc1vd5TW/bsMe+ZcaxX5vrAO9SlCRJrbOYOxSHDGq7gPVJ1iU5ATgf2DFkeTcAb0/ywiQvBN7e23ZMBi5JkjRxquoIcAnTQekO4Pqq2pvksiTnAiR5fW861XnAp5Ls7R37MPBfmA5tu4DLetuOyUnzkiSplZp+tE9V7QR2zti2te/9LqaHC2c79irgqmGvZeCSJEmt1KWV5g1ckiSpdbr28GoDlyRJaqUudbi6Ex0lSZJayg6XJElqpS51uAxckiSplZzDJUmS1KB5LGC6LHQnOkqSJLWUHS5JktRKDilKkiQ1rEtDigYuSZLUOi58KkmStATscA2h7an0sssuG3cJA1XVuEuY0/vf//5xlzDQE088Me4SBvrkJz857hIkSQ2zwyVJklrJDpckSVLDDFySJEkNcuFTSZIkzYsdLkmS1Epd6nAZuCRJUisZuCRJkhpm4JIkSWpYlwKXk+YlSZIaZodLkiS1js9SlCRJWgJdGlI0cEmSpFYycEmSJDWsS4GrO4OjkiRJLWWHS5IktU7XnqVo4JIkSa3UpcA1cEgxyWlJvprk9iR7k7x3KQqTJEmT7WiXayGvthmmw3UE+P2qujXJ84A9SW6sqtsbrk2SJKkTBgauqnoAeKD3/h+S3AGcChi4JElSY9rYqVqoec3hSrIWeC1wSyPVSJIk9Uxk4EryXOAvgPdV1Y9n+X4LsAXg9NNPH1mBkiRp8rR1LtZCDbUOV5JnMR22PltVX5htn6raVlVTVTW1evXqUdYoSZImUJcmzQ9zl2KAK4E7qurjzZckSZLULcN0uN4EvAt4a5Lbeq9fbbguSZI04brU4RrmLsW/AdpXuSRJ6rQ2BqeFcqV5SZLUSgYuSZKkBrV1aHChhrpLUZIkSQtnh0uSJLVSlzpcBi5JktRKXQpcDilKkqRWanpZiCQbk9yVZF+SS2f5/sQkn+99f0vvEYckeVaSq5N8O8kdST446FoGLkmSNHGSrACuAM4BNgAXJNkwY7eLgUeq6gzgcuBjve3nASdW1auBfwL8m6Nh7FgMXJIkqZUa7nCdBeyrqrur6kngOmDTjH02AVf33m8Hzu49gaeAk5McDzwbeBL4uedM9zNwSZKk1llM2BoycJ0K3N/3eX9v26z7VNUR4FFgJdPh6zHgAeA+4I+r6uG5LuakeUmS1EqLnDS/Ksnuvs/bqmrbIks66izgKeClwAuBryf5X1V197EOMHBJkqRWWmTgOlRVU3N8fwA4re/zmt622fbZ3xs+PAV4CHgn8OWq+inwYJK/BaaAYwYuhxQlSdIk2gWsT7IuyQnA+cCOGfvsAC7svd8M3FRVxfQw4lsBkpwMvBG4c66L2eGSJEmt1OQ6XFV1JMklwA3ACuCqqtqb5DJgd1XtAK4ErkmyD3iY6VAG03c3fibJXiDAZ6rqW3Ndz8AlSZJaqemFT6tqJ7Bzxratfe8fZ3oJiJnHHZ5t+1wMXJIkqXW69vBqA1eLtf0ftMsvv3zcJQz00EMPjbuEge68c85h/1Y488wzx12CpAnU9v8OzoeT5iVJkhpmh0uSJLXSccd1py9k4JIkSa3UpSFFA5ckSWqdrk2a706vTpIkqaXscEmSpFbqUofLwCVJklrJwCVJktQwA5ckSVLDuhS4nDQvSZLUMDtckiSpdbq2LISBS5IktZKBS5IkqWEGLkmSpIZ1KXA5aV6SJKlhdrgkSVLrJOG447rTFxoYuJKcBNwMnNjbf3tVfbjpwiRJ0mTr0pDiMB2uJ4C3VtXhJM8C/ibJX1XV3zVcmyRJmmATFbiqqoDDvY/P6r2qyaIkSZK6ZKjB0SQrktwGPAjcWFW3zLLPliS7k+w+ePDgiMuUJEmT5ujipwt5tc1Qgauqnqqq1wBrgLOSvGqWfbZV1VRVTa1evXrEZUqSpEkzcYHrqKr6EfBVYGMj1UiSJLG4sLUsA1eS1Ule0Hv/bOBtwJ0N1yVJkiZclwLXMHcpvgS4OskKpgPa9VX1l82WJUmS1B3D3KX4LeC1S1CLJEnSM9rYqVooV5qXJEmtZOCSJElqmIFLkiSpQV17lmJ3/iaSJEktZYdLkiS1kkOKkiRJDTNwSZIkNaxLgcs5XJIkSQ2zwyVJklqnrY/oWSgDlyRJaiUDlyRJUsMMXJIkSQ3rUuBy0rwkSVLD7HBJkqRW6lKHy8ClTlu5cuW4SxjoZS972bhLGOg973nPuEsY6BOf+MS4S5A0Qt6lKEmStAS69PBqA5ckSWqlLnW4uhMdJUmS5iHJxiR3JdmX5NJZvj8xyed739+SZG3fd/84yTeS7E3y7SQnzXUtA5ckSWqlo/O4FvIa4twrgCuAc4ANwAVJNszY7WLgkao6A7gc+Fjv2OOBa4F/W1WvBN4C/HSu6xm4JElS6ywmbA05FHkWsK+q7q6qJ4HrgE0z9tkEXN17vx04O9Mnfzvwrar6PwBV9VBVPTXXxZzDJUmSWmmRk+ZXJdnd93lbVW3r+3wqcH/f5/3AG2ac45l9qupIkkeBlcA/AirJDcBq4Lqq+qO5ijFwSZKkLjpUVVMNnft44JeB1wM/Ab6SZE9VfeVYBzikKEmSWqnhIcUDwGl9n9f0ts26T2/e1inAQ0x3w26uqkNV9RNgJ/C6uS5m4JIkSa3UcODaBaxPsi7JCcD5wI4Z++wALuy93wzcVFUF3AC8OslzekHsnwO3z3UxhxQlSVLrNL3SfG9O1iVMh6cVwFVVtTfJZcDuqtoBXAlck2Qf8DDToYyqeiTJx5kObQXsrKovzXU9A5ckSWqlphc+raqdTA8H9m/b2vf+ceC8Yxx7LdNLQwzFIUVJkqSG2eGSJEmt5LMUJUmSGtT0HK6lZuCSJEmt1KXA1Z1enSRJUkvZ4ZIkSa3UpQ7X0IGr91Tt3cCBqnpHcyVJkiRN7qT59wJ3AM9vqBZJkiSge5Pmh4qOSdYAvwZ8utlyJEmSpjX8aJ8lNWyv7hPAHwBPH2uHJFuS7E6y++DBg6OoTZIkqRMGBq4k7wAerKo9c+1XVduqaqqqplavXj2yAiVJ0mTqUodrmDlcbwLOTfKrwEnA85NcW1W/1WxpkiRpkrUxOC3UwA5XVX2wqtZU1Vqmn5J9k2FLkiQ1KQnHHXfcgl9t4zpckiSplbrU4ZpX4KqqrwFfa6QSSZKkjrLDJUmSWmliO1ySJElLxcAlSZLUoKOT5ruiO38TSZKklrLDJUmSWskhRUmSpIYZuCRJkhpm4JIkSWpYlwKXk+YlSZIaZodLkiS1TteWhTBwSZKkVurSkKKBS5IktZKBS5IkqWFdClzdGRyVJElqKTtc0pgdPnx43CUM9Nhjj427hIHuu+++cZcw0Ite9KJxlzCnk046adwlSM9w0rwkSdIScEhRkiRJQzNwSZIkNcwhRUmS1EpdGlI0cEmSpFYycEmSJDUoiYFLkiSpaV0KXE6alyRJapgdLkmS1Epd6nAZuCRJUisZuCRJkhrWpcDlHC5JkqSGGbgkSZIaZuCSJEmtc3QdroW+hrzGxiR3JdmX5NJZvj8xyed739+SZO2M709PcjjJBwZdy8AlSZJaqcnAlWQFcAVwDrABuCDJhhm7XQw8UlVnAJcDH5vx/ceBvxrm72LgkiRJrdRwh+ssYF9V3V1VTwLXAZtm7LMJuLr3fjtwdnonT/IbwD3A3mEuNlTgSnJvkm8nuS3J7mGOkSRJarFTgfv7Pu/vbZt1n6o6AjwKrEzyXOA/AP952IvNZ1mIf1FVh+axvyRJ0oItclmIVTOaRNuqatsiSzrqI8DlVXV42Bpdh0uSJLXSIgPXoaqamuP7A8BpfZ/X9LbNts/+JMcDpwAPAW8ANif5I+AFwNNJHq+q/3asiw0buAr46yQFfGqECVGSJOnnzOduwwXaBaxPso7pYHU+8M4Z++wALgS+AWwGbqqqAv5ZX50fAQ7PFbZg+MD1y1V1IMmLgBuT3FlVN/fvkGQLsAXg9NNPH/K0kiRJs2sycFXVkSSXADcAK4CrqmpvksuA3VW1A7gSuCbJPuBhpkPZggwVuKrqQO/PB5N8kemZ/TfP2GcbsA1gamqqFlqQJEnSUqiqncDOGdu29r1/HDhvwDk+Msy1Bt6lmOTkJM87+h54O/CdYU4uSZKk4TpcLwa+2GvrHQ/8eVV9udGqJEnSxOvSw6sHBq6quhv4xSWoRZIk6RkTFbgkSZLGoUuBy0f7SJIkNcwOlyRJap0lWIdrSRm4JElSKxm4JEmSGtalwOUcLkmSpIbZ4ZIkSa3UpQ6XgUuSJLVSlwKXQ4qSJEkNs8MlSZJap2vLQtjhkiRJapgdLkmS1Epd6nAZuCRJUisZuCRJkhpm4JIkSWqYgUvSRDn55JPHXcJAy6HGj370o+MuYU6bN28edwkDnXnmmeMuQVoQA5ckSWodl4WQJEnSvNjhkiRJrWSHS5IkSUOzwyVJklrJDpckSZKGZodLkiS1Upc6XAYuSZLUSl0KXA4pSpIkNcwOlyRJah0XPpUkSdK82OGSJEmt1KUOl4FLkiS1UpcCl0OKkiRJDTNwSZIkNcwhRUmS1EoTN6SY5AVJtie5M8kdSX6p6cIkSdJkO7o0xEJebTNsh+tPgC9X1eYkJwDPabAmSZKkThkYuJKcArwZuAigqp4Enmy2LEmSNMna2qlaqGGGFNcBB4HPJPlmkk8nObnhuiRJkjpjmMB1PPA64E+r6rXAY8ClM3dKsiXJ7iS7Dx48OOIyJUnSpOnSHK5hAtd+YH9V3dL7vJ3pAPYzqmpbVU1V1dTq1atHWaMkSZpAExW4quqHwP1JXt7bdDZwe6NVSZIkdciwdyn+HvDZ3h2KdwO/01xJkiRJ3VqHa6jAVVW3AVPNliJJktRNPtpHkiS1UtNzuJJsTHJXkn1JZrsh8MQkn+99f0uStb3tb0uyJ8m3e3++ddC1DFySJGniJFkBXAGcA2wALkiyYcZuFwOPVNUZwOXAx3rbDwG/XlWvBi4Erhl0PQOXJElqncV0t4bscJ0F7Kuqu3uLul8HbJqxzybg6t777cDZSVJV36yqH/S27wWeneTEuS5m4JIkSV206uj6oL3Xlhnfnwrc3/d5f2/brPtU1RHgUWDljH3+FXBrVT0xVzHD3qUoSZK0pBZ5l+Khqmr0hr8kr2R6mPHtg/Y1cEmSpFZqeFmIA8BpfZ/X9LbNts/+JMcDpwAP9WpbA3wR+O2q+t6gizmkKEmSJtEuYH2Sdb11Rs8HdszYZwfTk+IBNgM3VVUleQHwJeDSqvrbYS5m4JIkSa3U5KT53pysS4AbgDuA66tqb5LLkpzb2+1KYGWSfcC/5/8/S/oS4Axga5Lbeq8XzXU9hxQlSdJEqqqdwM4Z27b2vX8cOG+W4z4KfHQ+1zJwSZKkVpq4R/tIkiQtpfmsGL8cOIdLkiSpYQYuSZKkhjmkKEmSWqlLQ4oGLkmS1EoGLknSvH3oQx8adwlz+vrXvz7uEgZ6xSteMe4SBlq3bt24S5jTgQMzF1PXUnAOlyRJUsPscEmSpFbq0pCiHS5JkqSG2eGSJEmt48KnkiRJmhcDlyRJUsMcUpQkSa3kkKIkSZKGZodLkiS1kh0uSZIkDc3AJUmS1DCHFCVJUis5pChJkqShGbgkSZIa5pCiJElqJYcUJUmSNLSBgSvJy5Pc1vf6cZL3LUFtkiRpQh19ePVCX20zcEixqu4CXgOQZAVwAPhis2VJkiR1x3yHFM8GvldV32+iGEmSpC6ab+A6H/jcbF8k2ZJkd5LdBw8eXHxlkiRJHTF04EpyAnAu8D9m+76qtlXVVFVNrV69elT1SZKkCdWlOVzz6XCdA9xaVX/fVDGSJEldNJ91uC7gGMOJkiRJo9bGTtVCDdXhSnIy8DbgC82WI0mS1D1Ddbiq6jFgZcO1SJIkdZKP9pEkSa00cUOKkiRJWjgDlyRJUsMcUpQkSa3kkKIkSZKGZodLkiS1TltXjF8oO1ySJEkNM3BJkiQ1zCFFSZLUSg4pSpIkaWh2uCRJUivZ4ZIkSVrmkmxMcleSfUkuneX7E5N8vvf9LUnW9n33wd72u5L8yqBrGbgkSdLESbICuAI4B9gAXJBkw4zdLgYeqaozgMuBj/WO3QCcD7wS2Ah8sne+YzJwSZKkVjq6FtdCXkM4C9hXVXdX1ZPAdcCmGftsAq7uvd8OnJ3pk28CrquqJ6rqHmBf73zHZOCSJEmT6FTg/r7P+3vbZt2nqo4AjwIrhzz2ZzQyaX7Pnj2Hknx/hKdcBRwa4fkmkb/haPg7joa/4+L5G47GyH/He+65Z5Sna8LLxl3AMPbs2XNDklWLOMVJSXb3fd5WVdsWW9dCNRK4qmr1KM+XZHdVTY3ynJPG33A0/B1Hw99x8fwNR8Pfsb2qamPDlzgAnNb3eU1v22z77E9yPHAK8NCQx/4MhxQlSdIk2gWsT7IuyQlMT4LfMWOfHcCFvfebgZuqqnrbz+/dxbgOWA/877ku5jpckiRp4lTVkSSXADcAK4CrqmpvksuA3VW1A7gSuCbJPuBhpkMZvf2uB24HjgC/W1VPzXW95RK4xjbm2iH+hqPh7zga/o6L5284Gv6OE6yqdgI7Z2zb2vf+ceC8Yxz7h8AfDnutTHfGJEmS1BTncEmSJDWs1YFr0JL7GizJaUm+muT2JHuTvHfcNS1XSVYk+WaSvxx3LctVkhck2Z7kziR3JPmlcde0HCV5f+/f5+8k+VySk8ZdU9sluSrJg0m+07ftF5LcmOS7vT9fOM4a1W2tDVxDLrmvwY4Av19VG4A3Ar/r77hg7wXuGHcRy9yfAF+uqjOBX8Tfc96SnAr8O2Cqql7F9GTf88db1bLwZ0w/gqXfpcBXqmo98JXeZ6kRrQ1cDLfkvgaoqgeq6tbe+39g+j9wc66Gq5+XZA3wa8Cnx13LcpXkFODNTN/1Q1U9WVU/GmtRy9fxwLN76wI9B/jBmOtpvaq6mem7zPr1P7blauA3lrImTZY2B655L5uvufWecv5a4JYxl7IcfQL4A+DpMdexnK0DDgKf6Q3NfjrJyeMuarmpqgPAHwP3AQ8Aj1bVX4+3qmXrxVX1QO/9D4EXj7MYdVubA5dGKMlzgb8A3ldVPx53PctJkncAD1bVnnHXsswdD7wO+NOqei3wGA7hzFtvntEmpgPsS4GTk/zWeKta/nqLWXrbvhrT5sA172XzNbskz2I6bH22qr4w7nqWoTcB5ya5l+mh7bcmuXa8JS1L+4H9VXW0w7qd6QCm+fmXwD1VdbCqfgp8AfinY65pufr7JC8B6P354JjrUYe1OXANs+S+BkgSpufM3FFVHx93PctRVX2wqtZU1Vqm/zm8qarsKMxTVf0QuD/Jy3ubzmZ6lWbNz33AG5M8p/fv99l488FC9T+25ULgf46xFnVca1eaP9aS+2Muazl6E/Au4NtJbutt+4+91XWlpfZ7wGd7/xN1N/A7Y65n2amqW5JsB25l+i7kb+Jq6QMl+RzwFmBVkv3Ah4H/Clyf5GLg+8Bvjq9CdZ0rzUuSJDWszUOKkiRJnWDgkiRJapiBS5IkqWEGLkmSpIYZuCRJkhpm4JIkSWqYgUuSJKlhBi5JkqSG/T8hE005ewiHoAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "P = np.zeros((n, m))\n",
    "\n",
    "# TODO: Calculate the optimal coupling P here\n",
    "\n",
    "plt.figure(figsize = (15,6))\n",
    "plt.imshow(P, cmap='Greys')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing the monotone rearrangent coupling\n",
    "\n",
    "There are two ways we want to test that ```P``` is correct:\n",
    " 1. That the marginals are correct, i.e. $ P \\mathbb{1}_m = a $ and $P^T \\mathbb{1}_n = b$\n",
    " 2. That we obtain the same results with the linear programming approach\n",
    " 3. (Optional) again, test for cyclic monotonicity\n",
    "\n",
    "### TASK 6: Test the marginals of $P$\n",
    "\n",
    "In the box below, test that the marginals of $P$ agree with $a$ and $b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TASK 7: Implement the LP method for this 1d example\n",
    "\n",
    "Now we take to implementing the LP method in this case. Much of the code will be the same as in the 2 dimensional case, except for the distance matrix, which is given for you below.\n",
    "\n",
    "In the following cells, you have to build the ```A_eq``` and ```b_eq``` matrix as before, with the new dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_matrix_1d(x, y, p=2):\n",
    "    return np.abs(x.T[:,np.newaxis] - y.T)**p\n",
    "\n",
    "C = distance_matrix_1d(X,Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 564,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Compute the S and T matrix with the dimensions n and m\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Build the arguments for the scipy.optimize.linprog routine and execute it\n",
    "\n",
    "results = scipy.optimize.linprog('''TASK: Fill in the arguments here ''')\n",
    "\n",
    "P_lp = results.x.reshape((n, m))\n",
    "\n",
    "primal_min = results.fun\n",
    "print(f'Optimal value / W_2 distance: {results.fun}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TASK 8: Compare the results of the monotone rearrangement and the LP method\n",
    "\n",
    "For example, you could compare the plans graphically, or you could compute the $\\ell^2$ norm of the difference of the two plans."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is the metric we calculate in these examples?\n",
    "\n",
    "As has been hinted before, what we are computing is the **Wasserstein metric**, which is what we call optimal transport in the specific case that $c(x, y) = \\| x - y \\|^p$, or in a more general metric space, $c(x, y) = d(x, y)^p$.\n",
    "\n",
    "Assuming that $c(x, y) = \\|x - y\\|^2$, we are computing the **Wasserstein-2 metric**:\n",
    "\n",
    "$$\n",
    "W_2^2 (a, b) = \\min_{P \\in U(a, b)} \\langle C, P \\rangle_F\n",
    "$$\n",
    "\n",
    "If we take our coupling/plan $P^\\star$ from either the explicit monotone rearrangement or the LP solver, we have $W^2_2 (a, b) = \\langle C, P^\\star \\rangle$, so both couplings if they are equal should give the same value, and furthermore their values should be the same as ```results.fun```, if we have ```results = scipy.optimize.linprog(...)```  above.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Monotone rearrangement: W^2_2(a, b) = \\t {((C * P).sum())}\\n')\n",
    "\n",
    "print(f'LP method: W^2_2(a, b) = \\t\\t {((C * P_lp).sum())}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}