{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Optimal Transport with Linear Programming\n", "=========================================\n", "\n*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_matlab/) for details about how to install the toolboxes.\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}$\n"], "metadata": {}, "cell_type": "markdown"}, {"source": ["This numerical tours details how to solve the discrete optimal transport\n", "problem (in the case of measures that are sums of Diracs) using linear\n", "programming."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 2, "cell_type": "code", "language": "python", "metadata": {}, "input": ["addpath('toolbox_signal')\n", "addpath('toolbox_general')\n", "addpath('solutions/optimaltransp_1_linprog')"]}, {"source": ["Optimal Transport of Discrete Distribution\n", "------------------------------------------\n", "We consider two dicretes distributions\n", "$$ \\forall k=0,1, \\quad \\mu_k = \\sum_{i=1}^{n_k} p_{k,i} \\de_{x_{k,i}} $$\n", "where $n_0,n_1$ are the number of points, $\\de_x$ is the Dirac at\n", "location $x \\in \\RR^d$, and $ X_k = ( x_{k,i} )_{i=1}^{n_i} \\subset \\RR^d$ for $k=0,1$\n", "are two point clouds.\n", "\n", "\n", "We define the set of couplings between $\\mu_0,\\mu_1$ as\n", "$$ \\Pp = \\enscond{ (\\ga_{i,j})_{i,j} \\in (\\RR^+)^{n_0 \\times n_1} }{\n", "     \\forall i, \\sum_j \\ga_{i,j} = p_{0,i}, \\:\n", "     \\forall j, \\sum_i \\ga_{i,j} = p_{1,j} } $$\n", "\n", "\n", "The Kantorovitch formulation of the optimal transport reads\n", "$$ \\ga^\\star \\in \\uargmin{\\ga \\in \\Pp} \\sum_{i,j} \\ga_{i,j} C_{i,j}  $$\n", "where $C_{i,j} \\geq 0$ is the cost of moving some mass from $x_{0,i}$\n", "to $x_{1,j}$.\n", "\n", "\n", "The optimal coupling $\\ga^\\star$ can be shown to be a sparse matrix\n", "with less than $n_0+n_1-1$ non zero entries. An entry $\\ga_{i,j}^\\star \\neq 0$\n", "should be understood as a link between $x_{0,i}$\n", "and $x_{1,j}$ where an amount of mass equal to $\\ga_{i,j}^\\star$ is transfered.\n", "\n", "\n", "In the following, we concentrate on the $L^2$ Wasserstein distance.\n", "$$ C_{i,j}=\\norm{x_{0,i}-x_{1,j}}^2. $$\n", "\n", "\n", "The $L^2$ Wasserstein distance is then defined as\n", "$$ W_2(\\mu_0,\\mu_1)^2 = \\sum_{i,j} \\ga_{i,j}^\\star C_{i,j}. $$\n", "\n", "\n", "The coupling constraint\n", "$$\n", "     \\forall i, \\sum_j \\ga_{i,j} = p_{0,i}, \\:\n", "     \\forall j, \\sum_i \\ga_{i,j} = p_{1,j}$\n", "$$\n", "can be expressed in matrix form as\n", "$$ \\Sigma(n_0,n_1) \\ga = [p_0;p_1] $$\n", "where $ \\Sigma(n_0,n_1) \\in \\RR^{ (n_0+n_1) \\times (n_0 n_1) } $."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["flat = @(x)x(:);\n", "Cols = @(n0,n1)sparse( flat(repmat(1:n1, [n0 1])), ...\n", "             flat(reshape(1:n0*n1,n0,n1) ), ...\n", "             ones(n0*n1,1) );\n", "Rows = @(n0,n1)sparse( flat(repmat(1:n0, [n1 1])), ...\n", "             flat(reshape(1:n0*n1,n0,n1)' ), ...\n", "             ones(n0*n1,1) );\n", "Sigma = @(n0,n1)[Rows(n0,n1);Cols(n0,n1)];"]}, {"source": ["We use a simplex algorithm to compute the optimal transport coupling\n", "$\\ga^\\star$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["maxit = 1e4; tol = 1e-9;\n", "otransp = @(C,p0,p1)reshape( perform_linprog( ...\n", "        Sigma(length(p0),length(p1)), ...\n", "        [p0(:);p1(:)], C(:), 0, maxit, tol), [length(p0) length(p1)] );"]}, {"source": ["Dimensions $n_0, n_1$ of the clouds."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["n0 = 60;\n", "n1 = 80;"]}, {"source": ["Compute a first point cloud $X_0$ that is Gaussian.\n", "and a second point cloud $X_1$ that is Gaussian mixture."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gauss = @(q,a,c)a*randn(2,q)+repmat(c(:), [1 q]);\n", "X0 = randn(2,n0)*.3;\n", "X1 = [gauss(n1/2,.5, [0 1.6]) gauss(n1/4,.3, [-1 -1]) gauss(n1/4,.3, [1 -1])];"]}, {"source": ["Density weights $p_0, p_1$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["normalize = @(a)a/sum(a(:));\n", "p0 = normalize(rand(n0,1));\n", "p1 = normalize(rand(n1,1));"]}, {"source": ["Shortcut for display."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["myplot = @(x,y,ms,col)plot(x,y, 'o', 'MarkerSize', ms, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', col, 'LineWidth', 2);"]}, {"source": ["Display the point clouds.\n", "The size of each dot is proportional to its probability density weight."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3d65Kq\nuAIGUDk17//KnB/2ViQXwj0Ja9WuqR6bRkTMZ0IuwziOLwBozf/uPgAA2EKAAdAkAQZAkwQYAE0S\nYAA0SYAB0CQBBkCTBBgATRJgADRJgAHQJAEGQJMEGABNEmAANEmAAdAkAQZAkwQYAE0SYAA0SYAB\n0CQBBkCTBBgATRJgADRJgAHQJAEGQJMEGABNEmAANEmAAdAkAcaJhmG4+xCAbg3jON59DPQmmluu\nNOBYamAcLFXrUhsDjvXf3QdAtz4VLsEFnEENjFOMiZ8BjiLAAGiSAAOgSQKMUwyJnwGOohs9x9ON\nHriAGhjHC7NKegGHUwMDoElqYAA0SYAB0CQBBkCTBBgATRJgADRJgAHQJAEGQJMEGABNEmAANEmA\nAdAkAQZAkwQYAE0SYAA0SYAB0CQBBkCTBBgATRJgADRJgAHQJAEGQJMEGABNEmAANEmAAdAkAQZA\nkwQYAE0SYAA0SYAB0CQBBkCTBBjtGYbh7kMA7jeM43j3McCyTGi5huGZBBi1K6xvuZLhaQQYVZul\n1+xinSWbixkexT0w6jVNrzFIr/BB98bgUQQYDchXrFS74JkEGJX6VKdK8umzjUoYPIcAA6BJAoyq\nlTcPakiEpxFg1GhPS6BWRHgIAQZAk/67+wBoQ1ityQ+6mo/fMkILOJoaGAuGYYg2yqUef8XSTrMe\ncDgBRs5i8OQ32Fzt2lNjU9uDh9CEWJ2FSLipdA6f9XOUw5CbkGwMJnxaZShOQVU8eBoBVpfCGs81\nMZYfSjwWZNghoVKSYc9Mr7U3JqEzmhBrEd5SGn//zTa+7MAyJWLqV2ExuqFg3VYWP6cEX3tjEvoj\nwGpUw8S1O/c/juMnSzaHyucPh3Qda/qrh6dX4W+hG5oQq1A+719Jw92V8re49h/eOI6fk/MNqtiT\n3n4qLpO5WgQXj6IGdr9Vs9a+rq2H1SBMJun1St+YnG0GHVMDa9LOrn0HuuYw3vkULZSfE12F6rk2\n4GwC7GZrq18fZ5dT07a7SsgqYEoTYvNO782x6VecSpLDSw2sEnWWR59KWHQk1vC75WVHtcE04ys/\n1ELl47uhYwKMIqd2NTxPdFbGmg94P9VinkMTIjmLZX3NYVA+13BDpmPj4OHUwFiQ6gFYc3S9Zs2G\nnwcnvz3v+C9bSmbWkNhQoy779dcwvoEAq0L9tzTa/YSMvz+fXXG5oNEyOrh7tsGBT0eFZpdZ9w3j\nKZoQm/fMC3ezU0/WZY2Tuen/XQ+9a7oN/FhqYDfL9/TLcAlfZkNbzazR8owvyN+bYRcuUEBV3m/5\neddY/dTAmiS9Nlt76sK2mgMP5hDTeZN5Dm/5S4DVYE+/MiVXSvSsXhw+D5wmn8sMk/8+lgCry+Ll\n+MzVQ3YaggVZTjp1Jeu/wFH0OxVgVZhefKmy75oiuCepU7Tq1M023jAwzjvFsVxjH0+871et8vsr\n3rULXDacC9hGgFUnH2Per2tYugXqpxt9dWbdo2cPcjHLa12puQlfuJcaGET8Da76/O+/H3xeTqLh\ngQ3UwJ5o1ffcJ38pno0uf84Lv9j8dmNQ633mKF0WCbDHWTUUN7rxo6Z+uLH98CG9SKLTLod1Xwjp\nRv9c47SYCOabeJtt3GcJGnN7T+XopMBXHsD1ouc3dYnCSw3ssaZfdaeLjOQ3nm7ffavOja8us5JZ\nr+e8z1fFydTAnm5VeqUe4SRjUPftrCJS8nJcb6QIsIfKTKT2qKbCmnkXQrOWbR5OE+JzhcWAEpOa\nfZfx/PdDrw2qFFIDe5z9n3lfgK/xhPNccjVmzoPa2MOpgT3Rz9zBny+zs20mzYy+5d4lUkvutM5R\ncplNNxBcvNTAKDH8/nvrtSStwSHz6Lcivx7e93qb/dVnA5WwB1MDIykzDWCXJWlVxnF8yEDmKVnE\nKgKMnHeROStWnlCS1uA55zlMayghwFj2qYo9p0jlYu9LS4yxigB7us+X38xddIXK20ltek+eLnkm\nfOHfrvM6ExHQiYOv/F30hzujcpAalmu47sdPj9nJDzoT8RJgvIIyIlpMvJQU/xx1FhYjSoa9pa5P\n6HZuUNayomDesQ190TVEvr896Fl6oqGVkADjK5phrpCPo+6B5dPrb5vdzwLdE2Dwes1CpWB96kOq\nX/ldyDDI0wsRIut5XjAXhlCCnQQYT5fqB3hSvae8a0ZmJpTprqLHWVihhKYJMPgTnXakTtN8mmVt\nGJCZnFv1RB8SkUoIMGjeJ8OO7Xmf2tsDJ2mkTsaBwZ/U8NhjU6G8uF981syOxsm/v72tfBWzrpKZ\ndbqHf1btH/ZTA+PpUjPJzprpPhun9rN2oNL+uZGSawXs2+0ru0pc5nnPu3EIUS643mje2exzr2ix\nMjE7q2vHgB8yDizaTePvJazcVWbn0T8I1+gy8ppbCLB+ZIrRJ7/LG2ZwSAVMtJguaTpbm2Gb8yAa\nPJsDLLP1kH6WVU8Ee7gH1ol8MfrY+xMbpsqdZkBqFeDo/sfgX+Ywwvn9wjWvMy9h0RD8cLY999tg\nGwHWg0wxGt3mIeatqdnfFkqV0dEaR75M39Cgt3x46Xlv1YrojwBrXqbG8Hpwhs1Oyzj5uUTJZoUz\nQuWeZRwLc2VbhuUfXH7SlY+/VMK4ll6I/Vuc0OFpPick7DV3RrG7eP5nQ7hSXTDK+/ipbPEQamCd\nyJdYTyvPNteNTi36C2+8zdzy3n27qPym789QuYuPCQICrG0aam53QTl+e1REe5ek3H60PIcmRKpT\nw0S0qVHG0y7p7+MsH4/cTcteJy+D9qmBta2bMrEe09azmUfVdj/nYdYFZvGCe9RZ4l4CrBML625c\ndBS1yITQ67cWld9PanjWvOvH0k52quHt80WJCgkwOjcfF1zwJ/lg+1ZNCqq/a7OnhqyaqfCQ4M1U\nUj0onbnuYe/12ikKF/9qeenI2a8Kni6/n3tHIm/ozPnYi41bqIF1ZdrYtarnWJcyZeiG9Fr87au4\nyXHxkLbt4XD5ltjQky82bqEG1olttY3nKFmVeE9d6qiZlCtcAbmwHnZ73PJAAqwrNRR/NXSC3+CQ\nJU6m+2noteeVzPoovbiFcWBdUXbst/MM9vcWzBZIK5kWC64hwDhSKzODpI5zcUauNl7e0WazNeY3\ng8toQmSj1NLPlTchLtwsXPzzz5b1vbTLWPWbSqiBsU4qAKb3fkp6TNwidTvnmfWqzSp8Z3km3ehZ\noaSFsIn0mk2JVN2xAgXUwCi12E/vZ92N4sWrrlHSybBcVS8NHksNjCIlATCr1tTToaM8vUwpCQ0R\nYCybDmVdrHrUXDfJHNvi/TDpBbURYCzYUJH6hkEFlbA9E/pFH9F+CJVwD4xS3Rfb02FeqeCVXlAP\nNTCKrC22q6qEFVo80MwUwG+HHxKQoQbGI6xdR3j8fVDFCyqkBkbOQ2oV36AKHikn5OBiamA8wnDm\nPTzRBbdQA4NIxUsiQf3UwKjR9TMChzPNuwF2hsrneqYtamDk3FLEzG687bwP950mf3HLyT/OEL6z\nD7nJykkEGEXWFjOtV19aP/4KLa5jAGsJMPpXXglbtRnlZrORqeZyCAHGgrWl/6vK6sv0VZTP2FvP\n8XepxdHuVEWAsWxVhu0visZxnCbHUSky3U80xgbpdTInlGPphcg67yJ+cT2w1+4AOCM/PqtFv2Wy\nVnpB/QQYRd4F+qf03zxt4O1mLyS1AWc4dTg5D1TXsrnUb/F2hSuKULimqNZa9hNgbBQmmWuJDLVe\nDifAgIv40sOxBBg1in5bd60CUwKM6uRvs7ligTe9EKnLb3rtWZ8L6JwaGBWZpFdupJmLllPN2gBc\nb9VSA6Mh4ZoncJj8dMNirEKmkqJCCyWFqfM43OJFZfGXCqmBAU+XuPP6/f11h8IaamDAoy2l12u6\n/ItKWFUEGBUyWxW3yF9XMqw6AoyGKDg42FLHV6omwKjIpGoVrtglvbidkKuLAKNmw+TfH+2H3E4r\nYiX0QqQu+fW6HptenxPy2DNwMme1SQKMGn2KaWNIX79x7oScw1qbTRJgVE1JrbWqQi7LSrgHBq0w\nGul2TntdBBg0wVf+U/x2fKUxAgyaoHi9QP4k60dTHQEGrRiUoWf4PZnRDHPmKyXAoGpKzAsEGTZE\nxyBSGwta0rZpd4a+L2Yd6M+W7xrjzFdIN3pogwL0bLPRh+Hj1EaAAfyQWK3QhEjz+mhbM3sWrCXA\n4GaLo5J9SCFKgMGdYun1/kjOH/dRhRkBBrcpWMz+tWcpGZ0R6JtxYHCP37WAM9GyMXXCut0wDCZR\npCcCDG6wMkhWz+Gb2VKG0Q3d6OFehRWscduUEIuzJEG71MDgRqubB0vqT9/lm3c+GdRNDQzWyUSI\nXhJwJQHGMkX222LtZxiKuvW6CwWH0ITIgnxp+5yyODFga96BMNrT76Duf6XfFb5z+s0OY/cRQFWM\nA+tHtIjc8/6GO8wsXtv3hbQ0YCt55iPn8OfxVSdtxZJUKs08gSbEHuT7TG8rsH6WKQl+G58r4hFm\nJ2N1b/VN78jqMz2O43MqxzyWGljzgnIqEi4b3uVUT7bIljuepQm/I47nv4z9KjU7VP5XuUP4+YM1\n5/k566XxQGpgbUs3bf3E2Npv/eXpNfurhxWRqWCbjdk66pxsGQr2sHeEZ9GJoxvRcurbxeC8BqW+\nC8hs9SsjtfHf45NcKX9fNAnCDwHWh1NCpO9kOs7G81ScYXIL4jQhNqy4UrWu6WlDXW3jNEedGNZk\nWOo8FZ0/7YEwJcCg3LbmxAPCvTy61o6m6GM9a55JgEGhIfHzK1YJy4VWoo97skWxMF0WZ6CXUnTG\nPTDmtvS5P+M4GhOPt9nJHMcxSK9wOo/5I6sm8E3v9kHTpvAQamANKx6sqtg6Sn6wV7zmNMuw7KQe\nGwekZ/e5MKBCtYx2qYH1IRNR29NL7mUXTc4vo1wSDCtaHYstdN+Hngiwti11xd44C0NqNtio7mfi\nKCv9F/LsLT2qLN+oOPvbwn2GR5jbDzRHE2JPTpm/Nd9JvOP0KivopyMIDpm0aXYzTNhAkhpY8xbL\nym2F6fSvUpU7hevEqvTa2NCn8gRTamA9yK/ccYh8wdlf9Wuzo0/FsZUwVTq6IsD6cUuKiK6Pc07F\nsXkjveiKJkTigg54RZ0UepIOpGHSgLpikckCAgZWUAMjLzqo6G14Vb+EykHNqtOOLPEBy8edh89z\nzfol5nro/HuZhR1uoBMCjIh7Owt8VyPbFwnRV7F1UqVh8QbSvuma4r0Zz1DzFw5YRYBRizBv9iRZ\nbGqlLatBTvaTqRKt3Xn5BPY/tbHUefithL2CnVuUmT4JMKqQr/OtrdwkxvZ+w2D2dJsqZLP9j/va\nVMO4XSedtdAtnTiIuPJ7+jAMicltr5uLdrbn4JA+B7b4yIJNi1iWVkPzG7znEV44PmiKGhh5qcau\nM7Ik+kRhj4b9O0zeG0s/MqZ3ON9zvhIW63Mxrc+FC7Ws8J0DbNBsSP8EGItyN2yO62dxQDYE+6zf\n55bVT6f839++XuvPs9DiCTQhEvdbAk6HPQ3ljVrHHtHfoRTcLdv3LKePewvO25ETMcNzCDCS8uXs\n/lK1uPq1/RmyD0ZXSLlI4R0pN64gQxMiOak1M+ssVWOHOmv/jHbNuK3JMTWJ5avsDGeqm3W+QXAs\nAcaCmorC8rBZHBdcPgzr9HjbdoaPHXgALdKESH8WewxuaDlcsbTn2bER6+V/9cADqIEaGA35K46L\ny+Vo/SmcRGO6//xcghVWaFKj077nSj2MXgkw+hb2Q1k70GrWGrkwKm7Vvat9U2QtDjx4NTHhMmym\nCZHbFE9LsVaYUjufJTqiIP9cSXua9c5oEhz+mT1y+BPB4dTAaEWq3e+1PplS/RLzt80i0zv9bHHH\nqLiCzTYOANcNhPqpgXGn9XMDRscXhzP2Zna4eezXwkyDpXs5IBIunanysueCtdTAqMfCUiBZ0dmY\nhoI8K3fYoKvL62q51x6MRZtXT91Co1pqYNwsmGxi+P23ZZfFDx7jmmrKJc8iqGiJAKMKZd/xDyle\nU70Q8x3oo+OuJlu03PEhfeTyjKppQqQW4VIgwSPlI7Hy7WbhOiaFO0/dPzs3um6Pxt3rf8IptG7T\ngLLBT7Ntot0L9yRB/pNyyjQcieha9RRFB/bviUpvQCo3qIEmRPqQSa9X/o5a2YzvpRscWFtK7+q8\nCtn0LJ26bA0cQIDRgN/e9tEp5zfu9rPn2qoUwXrQ29YnKz0zieXfwgMYz4hq2MY9MFpU3ukgbDxM\n3h4LMyzRsJZy2NT1QXpFNlmVZ4Vrjx3RYgnXEWC0IV28fjf590P0htm8FK6typWQWcxscerhf7vY\nNMg6e7avWGUGFmlCpAPjUnpNnV7yHte2thg8h8/6CC1RA6MZsUpYpoIS1UDJXpB/s44q1/cVbOA0\n8gRqYLQkKJFTc3YcsOhJ9A+P2OwoCy/huPQSV1RKDYz7peoc0SI4mLsvvsvM786++3Xa/qfzPY6/\njxz/1L/13cjQ70ZuItIzA5m502Jz2drr88wdFg1kXnyK6RGmtkz3fry6tTBzPhUd3E6AcZv0Da15\noblv5eKNe4jtZ7ljSOaJomEQ3T6bmqU7OUp42AoNKiHAuEfZUKd/v771Kl061MhxhtWsVVWZhk4O\n3EgnDm5QVq2pZdKH7GovG4JkeVqN4BmBCJ04uFcbtYd8LWpNHWg2tnp47zbcQ9CHIndg8EwCjKsV\nd4v4bJMs5S82ayEsP55tNcjFyUf23xrcuTe4lwBjNatDLb7k9JxMW+YwPKobRT4O16Yy3E6AsU5Y\nCCr4VlpdG9t/bmPRFe/zWUNNFwoJMDaLrrx10hMlnyVcvpklY+J//86kDKMVAowVYrevLpiY/Fuw\nzn4xa6mrvtgNz9XGO3wlo6FjG+e3HGUYbRFgNGxVOX7s0xU+b3AD7O+PVu0kfxj5sNnWZQaaYBwY\nV/tdXnn1X0cHUZ2dXsMwpHpAZHpGpH/1GUC2Jy3Om3Hx9bp77B2UEGCs8Js90WngDzetQMwGUV3U\nxpXtAfG3Qapvy++fhFPpbz6Anz9P5evK6teHxkPaoAmRA6ytAE1GOOW7lZesjHXuKLH0pCG5Piyr\n/mppSt+fn/OthW5c8ShqYKwTnTNi3y5T1biS6abOVTDl1eJ8V6m/WvXUmUMaC45hI62IVE4NjNUO\n+ZofTDORKStL+86daVsni8xf/XT5e+XGPqduGYa94VPnQbWMPqmBcZtxHBtp8lo8yPM6PoTNj+HP\nP2LHoCJFnwQYNxv/CR+865Bu8QmebAXu0nPytLeA5mhCpBbp4nJnR4+NcwleeAdoQyvoujkVoUtq\nYHRuqX950vX1j32RGemmuHXInSZH2iDAqFdB+bvQuTzWZ2/BMFF6oL/KYqNkVo783axp780DI0d6\n0QwBRiuGdPFdYpz9EE7IlA6txSc6b0LhxY6aP+chO8hhXTK5AUb9BBhVC4rRyPQfBUXtwnSFBZWt\nzaV/GLRDpvoV/PkrfMmpHi4F5yGf+sOZSQzHM3R/7qjFAznWhol0J38Vn8MptuRxuMPFsczxQr8g\nFOd1o5I9HLeU5Swj9z4FXE+AfR2+fDu3S3RJL58yI7p9fBKp6BWSvqgiLXv5PZyzrGVwWK5z2iHA\n/sw+29F+zc5Vc4o7YpQ0vqX/eMWKXLOnW0ivw/mWRk8E2Ov1+6nONCFtW7TJGb7Xmqa8hT3F/7js\n/ZUccDgB9npNwya1wb8ftpVWTvK9wrdj6zoje7+UHHhbCzATx5dSpFfHJcTe1VuqyiqNBLRON/pT\nKBEa4W16lp1D1KmNACuyp5wTZgBn0IRYZMMXNrnVqWuWH7vC+17gAy/UB77kXgmwr5Pm9/7p4uiT\n07xO0uvtaRfk015v9wTYj2iGbetG/4p1OQvX3qV13lO4iwB7vX77Vc8ybPP37cyN4me221TJqlqv\nl+6IWVpQaibA/swyLLpB+d5SI6O7an56MAXZE2hBqZ9eiF+Z63LzJTsm/ldH3nsduNJjam542pXv\nau/DWw81sB9KokcqbEjstp3tgleUKvQrP5nxVQzcBaiGGtiJXOA1C9bcyvOle7u2ajOZieXG2Gbc\nSICdyAVeubIMW1jymLzZ/eDPv+gG9Ui9zd7+qgiwS9X4SX22xNrH03+pjVk2rc2k7ge/asqweo6E\nEgLsFJ+SbloEDrENuF3Je6GnxgaLizzUmWGLXAf10InjLJl++YrC2ny/cATFqDfrVP1MzMUdBNiJ\nfted+j54y8FQwruzUzhSqq0TGv3MzkjcegiwcykQeY5vk8OaXuZ1VsLM0dIE98AAIqKxunlmVM6g\nBtaGRseBwqLaql/h3es664i8Xi/jyRuw2Ch/xpto2jd2WuyF+LfZvx+qutgyH7qqjvPh1MBqVzIv\n8FET26RmL/17dp9bjlZtzeZ9teuEVTk1sKql0uu7weTnnW9l4UAcFwzl8pWwGkZGWkpmqrl6pwCr\n12J6/W322WbrWzm7ajNLeu58Ih6opAnh9vT6HEa0BH/IBX/LrYqdBFi9Cm8hvPZlWGFMvuoocWhR\nvmS88VpaNf1Hx9d8/itszd9fdaPnz+JV2ejEP9zujJX2DpF69nDG4efIz8FfGzWw2yz2j/jrB1iy\nq+x+So6h/M/q7DMGm737QGU+C31f89e09JxEDawKDVVoLCpNZ/LpRc10o79HDR3WfWihRMcDmVcV\nAp/zUM+C1GpgNZqGWYWfnCquXODxBNg9Fr+/lDfQ7WyY3pNGWhGBGwmw28zyJuz4NF0VE7jLd03a\nYejyS1v5t9jaWl8E2J2idabFOQumD1bYLwiak/myGHmkuxgrfzG1vWwBdrMweFLjhYfg3yFquyLh\nFtMMK/mUdZZhjdILsRapT0NqUtHZBrdQ7eMCl138meWYo5Ng1dMZ77EEWKXmd8hiMXbbDHK3PCtP\nUli/OXzRnzDDojNTdPMR+LzekhWoK7xh4RvEbWqYOnPbULAKr2N6Ep9Rd7ZNdg+revl+CvHZ18QW\nlzHbpuQl1/l63QO7R8kXzCsb2Tfcxa3qOqYbkSW4EtWgXGmb/uyEXTC+6y8XT2zdq+htvwPvuB9O\ngN1g9jmZ/UtteYZZJ5FF1V7H9CofJPkM2/DxeWbXjNmX0Uw3ltq+troHViR5a/eE5UveD172MZo2\n+ufbwS2nwsVWzW803X7x4xPdMtOJI/+8Hai5s1iGe2ALDr9TVefcz5kb1+Hr3/6SXWws2XtfNnxk\nVr343X+4WeEqEH23pbfymVUDy4k2xw+xbQrf5mrnz51998z36S8X3myo/PNAJS64SrQoZLRyQtwD\nS4reqZr+fPYCjxcvXDKOY37hwcyiZeWeeY+BQpsvj8wyjNsGn+SPw0VciW5rYJHK05rvFIX9kaYV\nsvLqRc3fbTacpdQLT7XVwHlKxjPN2htS13zhrtYcHQfrM8Ci3+O2tV8tX8Hpovn6Jb5q1tlNb5qW\nb2zIdGva3/B4YI8wOmxCzA8BWbWHwgvq4ra+eoz/pH77/uHYyRthlS1fW9NzkO60v3Riqs8a2Fu+\nK92qPz9WSdNEr3zN5Bo70yDfpX7n+Jm/nezrEcarvwBLVZ7Obr/K73/nFdnuBV14swGOcuAVNutY\n/zruAh6DH1S+tuktwPY7sCK/83Lv45oWWpRbNbfs1JD4+ZDLb/9OFu9KTDs6GW1SrsN7YDUrX2S5\nj/SCi1X4wdkw+tP9sEK9BVgqITZfDkPwb6eSDDPEkscq/5L3Ed0yP66xTo0dbgV6bkLcFjbTRozU\nbo+6zhbncGvuEwgHKvms7Z/qjHZ1GGC5ZVU3XdlhuuRHhyzvMJg3M/UVMnwwvlSSTywdKZxguq1v\ne1UfXLN6a0J8i17KO9PrdfQluHbepszyENtWjoBqpUZivRKN+ZWnV7lOXsZVdHeJy9x3HYLHL5iX\numQQyakHALe4ZeHybdMnLu6tZC99T3J/uA6bEK8R/VSdtAZB+UKxOuDSmcw6VSdd6uFzGV9cLQG2\nTmE73YFBkvr6lpmA0SeNzlx2Sd/bFO82wFp93gPbL9qXN9JkF/z7bnzEfankTa/gGIADhR/q7Uu9\nrB8Y4DtoIQG2LDUCLHqJnXTd5Qfwv4LsPOcooGeZiegOe4qCX0mvcgIsKddLcGmRsMfOTw+Eop0q\no//LKu6B5YSzeW6w/6aU72PQunB8ql7E+6mBFZleWGsXCTvVEPwA7JSfiG6YWLXb/OxW0msDNbAm\njaILjhbOABJ+uOa1qPXtK4LqQGpg66y69A6fMjG/c/04YKfZzarZr6IfK5+1G6mBNcygZjhcGFSz\nz9HOpd45kABrw7ZpuYENyufR/vzKl8VbCLB1Vq2lcsjiYYUNFLOhzcIMVlmcrerwZS7YT4DVbrY+\nWfSDctTqnfBA+arVK8innZ8vKyIdSCeOBqRGQRoCCTvNb3clpoU76emmj+sMsoGm2xXuXRbhlnUl\noG/5D/X0U7yz5lSYTz7FqwiwdQoz7LxpzbQ/wFFKPs752NkWYKlVBtfuE/fANkrdkTq7FcDFDfUo\n7H+4mJQ6Xm3jHtg64e3c2b/MxkCFyieHO/UOmRnAN1ADWy2zROxsG6BX5XWmVffOWUUNbKPUvJz5\n+TqB2hSt0RX77eGfcwXHWmpgu8gqwEwcd1EDA55u500n6XUXAQY8VGbu+RKH97XQeWMtAQY8VzjN\nzSvWqXhPtKifnUeAAY9WONXhNNKmP68Yy1z2W4FXToABT5fKjLBT8bbZRz87yaQjG+g8A7AsHPq5\nqvCM/Pn78R37RIABXMF83IcTYADXqWc+7u8UIc2mgAADeJZ6QnQnnTgAHiSzqObFR7KfAAN4olPX\nnr6GAAN4iujU+O01Hf4jwABokgADoEkCDOApFucEaYsAA3ii6czFb831pBdgAA+ycxGZqhjIDPA4\nsx7zjQaBAAOgSZoQAWiSAAOgSQIMgCYJsNoNw9DcBGUAF2LbABcAAAEwSURBVNCJo1Kp0PJ+Abyp\ngdUoU+VSGwN4E2DVmUbU2MWSBwBnEGB1maXX1DTGZBiAAKuUO10AeQKsRpn0EmwAbwIMgCYJMACa\nJMBqpIcGwCIB1hjZBvAmwOqSX/BbegF8/Hf3AZD0jqtx8vOHCaUA1MCqMwunQXoBxJjMt17hdBve\nLIAPAQZAkzQhAtAkAQZAkwQYAE0SYAA0SYAB0CQBBkCTBBgATRJgADRJgAHQJAEGQJMEGABNEmAA\nNEmAAdAkAQZAkwQYAE0SYAA0SYAB0CQBBkCTBBgATRJgADRJgAHQJAEGQJMEGABNEmAANEmAAdAk\nAQZAkwQYAE0SYAA0SYAB0CQBBkCTBBgATRJgADRJgAHQJAEGQJMEGABNEmAANEmAAdAkAQZAkwQY\nAE36P5dadnDr0jsyAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "for i=1:length(p0)\n", "    myplot(X0(1,i), X0(2,i), p0(i)*length(p0)*10, 'b');\n", "end\n", "for i=1:length(p1)\n", "    myplot(X1(1,i), X1(2,i), p1(i)*length(p1)*10, 'r');\n", "end\n", "axis([min(X1(1,:)) max(X1(1,:)) min(X1(2,:)) max(X1(2,:))]); axis off;"]}, {"source": ["Compute the weight matrix $ (C_{i,j})_{i,j}. $"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["C = repmat( sum(X0.^2)', [1 n1] ) + ...\n", "    repmat( sum(X1.^2), [n0 1] ) - 2*X0'*X1;"]}, {"source": ["Compute the optimal transport plan."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gamma = otransp(C,p0,p1);"]}, {"source": ["Check that the number of non-zero entries in $\\ga^\\star$ is $n_0+n_1-1$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["Number of non-zero: 139 (n0+n1-1=139)\n"], "output_type": "display_data"}], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["fprintf('Number of non-zero: %d (n0+n1-1=%d)\\n', full(sum(gamma(:)~=0)), n0+n1-1);"]}, {"source": ["Check that the solution satifies the constraints $\\ga \\in \\Cc$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["Constraints deviation (should be 0): 5.38e-16, 1.01e-17.\n"], "output_type": "display_data"}], "prompt_number": 13, "cell_type": "code", "language": "python", "metadata": {}, "input": ["fprintf('Constraints deviation (should be 0): %.2e, %.2e.\\n', norm(sum(gamma,2)-p0(:)),  norm(sum(gamma,1)'-p1(:)));"]}, {"source": ["Displacement Interpolation\n", "--------------------------\n", "For any $t \\in [0,1]$, one can define a distribution $\\mu_t$ such\n", "that $t \\mapsto \\mu_t$ defines a geodesic for the Wasserstein metric.\n", "\n", "\n", "Since the $W_2$ distance is a geodesic distance, this geodesic path solves the\n", "following variational problem\n", "$$ \\mu_t = \\uargmin{\\mu} (1-t)W_2(\\mu_0,\\mu)^2 + t W_2(\\mu_1,\\mu)^2. $$\n", "This can be understood as a generalization of the usual Euclidean\n", "barycenter to barycenter of distribution. Indeed, in the case that\n", "$\\mu_k = \\de_{x_k}$, one has $\\mu_t=\\de_{x_t}$ where $ x_t =\n", "(1-t)x_0+t x_1 $.\n", "\n", "\n", "Once the optimal coupling $\\ga^\\star$ has been computed, the\n", "interpolated distribution is obtained as\n", "$$ \\mu_t = \\sum_{i,j} \\ga^\\star_{i,j} \\de_{(1-t)x_{0,i} + t x_{1,j}}. $$\n", "\n", "\n", "Find the $i,j$ with non-zero $\\ga_{i,j}^\\star$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 14, "cell_type": "code", "language": "python", "metadata": {}, "input": ["[I,J,gammaij] = find(gamma);"]}, {"source": ["Display the evolution of $\\mu_t$ for a varying value of $t \\in [0,1]$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3dbZKj\nOLOGYXzi3dH0mqbXNLOm6TVxflClkpVSKvUBVsJ9xcREtY0xBsxjffLa930DAMCb//v0BgAA0IMA\nAwC4RIABAFwiwAAALhFgAACXCDAAgEsEGADAJQIMAOASAQYAcIkAAwC4RIABAFwiwAAALhFgAACX\nCDAAgEsEGADAJQIMAOASAQYAcIkAW8Lr9VKeOly5PZil48hyxN0pHazqoeQoD/rfpzcAmtfrte+7\n/Bve6Uf2+CdH3LXql5f0GkcJ7POO85iz+X76jiyh5UvfUebXyRSUwD5v3/fjbE6+A5zf3o0cWS5w\nXoTjGx9ljt01CLCFcNLfVdORPa6DnAy+dBzi8H+OdTcCbCGUwO6q9chy6N1pOsQ0bM/C7luC8qM7\nfDE4Uh7ZjyzVyH6VjrL88iaJRYANYvcBAFyiFyIAwCUCDADgEgEGAHCJAAMAuESAAQBcIsAAAC4R\nYAAAlwgwAIBLBBgAwCUCDADgEgEGAHCJAAMAuESAAQBcIsAAAC4RYAAAlwgwAIBLBBgAwCUCDADg\nEgEGAHCJAAMAuESAAQBcIsAAAC4RYAAAlwgwAIBLBBgAwCUCDADgEgEGAHCJAAMAuESAAQBcIsAA\nAC4RYAAAlwgwAIBLBBgAwCUCDADgEgEGAHCJAAMAuESAAQBcIsAAIOP1er1er09vBTSvfd8/vQ0A\nsBCZW1wn10QJDAB+ZEtdFMXWRAkMAL7EQfX39s/xx7/b7+MPrparoQQGAKmQXsnfWAoBBgBwiQAD\nALhEgAFAKrR7JX9jKXTiAIAfSodDrparoQQGAHWk14IIMAD4Uip+kV5rIsAAYNsKg8CwMtrAAOAn\nvZLoKo1i7phuKnkJ195x//v0BgD+MFeeL8fxmnWMStWMyruUpqfitBlEgAFtuBi5IA9T/Ej3wapO\nitgUlmefNrf/pUWAAQ3CFeGv7b/jjz/br/DUza4OfllipuNg6Y1k8XCxeP3ZyslZY8uUiHrCL61b\nfRjgbDLAtijD+DZ9XHLV1mNmez9kSdKUMkbv3yHbzPTWtZPKgvu+lz7Onc5SSmBAszi9sKZSzOjh\ndBgpHv29/XO8PCnr0K3xDHSjB3ATpbKOFBZobRgL+ZeEnHwkWf/c+ajiTxr+yy5wb1QhAlbxRUG2\ngW33qpxxx55egb2LfJNQCMs+pb+10Ui72p3OUqoQgR5xbh1GrgvhenTSxWVuP/Lbkzvq2IFKMsWU\nZZQWuOniZrBr3vF6lMAAq2qz+ax1TvxWnrHNawoZ0/QqY0+KbN+QkEayXFVaMjblF4/yFkkvkjsd\n64ASGNDmqDyUJbCJpvR1ttSD3ebSNljpN2WHl9JryxXdjrebUkKqFgpvcHxLCDCgzZ/t11/bf9k2\nMEmvG1Q65d9svM7NxJmhlPn+3X6Xikcd03lAIsDOwom4uJHbPh0ZtpXTq2kaiDM65Vf7mxyPn1qO\nfJRSdWJf58Pq1SPbvpUs0PG+7tCNfr7X6xXOrde3z24SYtUjYrmnxp/tlz29mhY4b5BZsmbSS3pF\ntsL3t3tEl+zsrj/VHVEPSa+NEth0Ss3AsmfVoybJTj5sdkKNrfwTuPrLN1v0kW+hnA/kSofqcanS\nC82xo2IwFK2yZSylPax1q6pBNWWOR6fWvap69H7GHzt26XPraXdP19PlkOSHsh9k8HesXy48ZWIq\nSxXi+LssZWQcWAdjl/p4+e53vM0xmo4Amym6asR79dwhPt2eVoee7TFRYkyR0j5UKuuMLU+D+99Y\nHLnTUZ4SYEksle4NNk7Js1KXRUgE2EyFANu6M+zUcYjVa9ydzg1L2SjRWhJSJqofNHGE2fhqV9aU\nYUlIyLFc1YmAS+IEapqVQ3nH+x2sKejEcZ3x9Co9OMP+/d/NnT0Pb1N6xd3xFX1dgfTz7d4XxGrS\n9JWllFxMnioNbc7O/KSMg4aOALtAT+S8F+bmR0uuuU5f5g6a0msk6prKXvY3ao2xUkrdNb3iz2Wf\nz2nLJVNTn4tjYT2H+trAmMNeRy/Ek8y69O/vf3917T3zArTP2/ibsOzwbK5U5+wI48kUciVNJ0Dc\nUe2uuRWLu9KMFGjCGOSOhS1z4VPYmoIS2EyzRmZcXvR5Zf++zfXug0VJmV6y2jC7jP6qjurE2xzN\nKsuHVYZkKfRSXTWTKE5NR4Dlvcr0F5a+OasNZxbb+fr+Dz2yx92STKXyWXZJbqRpp8SYHiSlmkCl\nYtDiePkRnBS/ZqEKMWWZRsHePN4948P3U68ndKw42/go11gybjRZczZmQtd5YwjJVItfzmSJdvYv\nY5YlaeLGLT3POvo9QkcJ7E3uFM/0oRgvTulreH/qJf6YU557ci+1DtkiePYolOIn/M1cGx8hT2lL\nX0FdnD0d95asLklZTcePuB+WjnkdTURyohdx1Qvr6cmkM0a83u+saBrFHKyQNHG/fNlH/35H6lTx\n2S67BY4M/Op4lb6SZFUc6CwC7IstvX4W/1qufe/p00114zjqOgJMn/MpmVxDviqecSP5W75cmSc+\nu8EEWDe96qJUipINV5amLHtzF+nVhzYwyXKudPZoH04v+V4/095zlitCY1V3K9RWyJLkwWwaZR8v\nLa8/jkF6g+jcyaLsTWjoQ4BtmzYF1Hn6yl6yT8dOhrWqZlh3eilvRxqtQ35NLMWyUOWo3Gd5KwSS\nbDuo5hbFLws6cXT7Oqu6OlPsU/usD/Wzeo7kbl6lxUp1gyP09ApFtCnvhQ72kMgGTxxp4UGZbfaJ\nQkgvI0pgsXPPlShdlJjJduioFteYPsMkrj6qzhAfJ0pHumQLXrIB7FhMLkmJ7WKlqkWlqKTfEqz6\nLnob20Z6GRBgn2XvfxiHX/Kq/fv/F0w05V62I6hlxHFWtTayo72N+sbPSuaS71iDnmfZpKS7fB8C\nbKTabbDck335K/d49l0ocvUrjGfoETKsNPyrdNPnrCS96IJ4PTkLoj3PjDlkGVnPwbXg1/q2/VzI\nWndFfkyYYeDXNjX80ptncky7KZcV+yxQs1Tv6cyBnksfJWZ5KssyG4gcUs3BtSDAtm1GgG3N8xUl\nCZR9oRJy2jhrjumI0kHUC0PZ5fV5Di09O5R35CifoaNQnmSPEkVTboWDGAG2bf3d6EfqHsPLZYOW\nZbUE2ImmZJh9yLO+/tJ7cZRPkhx9vWNhqVIxW5AKa65OOsXBNSLAtq0zwJTyU1qnFz2YfZWlNJYg\nwK6TXNHsbVpJSaupMSx5SbI8h/g81dFa1SaxOJNk8UufsIoj24ROHNs2Ovu7MgeVsYhWrSekv8Yn\nJZXD9orBUh9FmUalMl/2vbjGnSocjvEbhrUeKY5sKwJMsmSYvdd76ZFsb0Nj/8P6RMOYS++1KEeM\nheCRxS9LR8TSYlzgLlBtqQr97Kvd5WduFnKYiePL+9mmnMHKDBqtBSaZlJb0sq8fk8UnSevwr+Mu\nKjKW4uJXfNtluQYuiF486v7Xn0WA/RAZlmSVjK7sOdqULt1RFCo8X7lNxVnCSaLPFq88bsfkUouL\ng2qPlBY+/ijd6JnM60AV4ptca4ReGvug4rsfH4Hvw0SDQ56PLvUjd/PiaF4mXATkiObDYEOXXAO6\nEWCpgTkaqu1V2eHMlpUoitMkMqfULOPptbW0gUkcx3tQeutwiPtwjRulXt068qlpho5Sf/2vvzm4\n42QfeuP4rWrvednBOlkPY5Y/KNuZflZ/dzlZD/oQYKNqASZHmA0Of95y65SPk2ETxAe3OrNGrDrV\nfXJc9DfSJypDh+Rr23qHMPb8IgiwIbXKpbkBtonyGQF2rrhs1DRmWZc9KNXKJXtNJgddoexGY4yx\ne9dBgA0RU3hM76ZYpRTL6MoxpJReHRNqBPrhmHg/Uo57FnPA3wydOPqNTUDVpNhT49M9IR8nqS0s\nlcz6Qs4yH3R2zbL8RxceqTRHVNyyxX7zhXFgnd4vNGEYVt/AZF3T14nv3olCTlQzKR62HEeOnk9N\nkwjLx+Nx0NX3erK4U4bs0b7mfnu9Xmtu2GcRYD0KZ1LpwekDjXfxn9wGzvUhlskPLewZpl+e7IOa\nyTBdtkthYqn9FkcXMZagCrHZ+wmU9KrIznBoly3AWdrV9vfoohzWT14gLHdJls/G0yEeZaOkiu/4\nw1Jh1Zqg1S4nD5Qc1ji9stWJcV3i9YO3qiV1ZdseVQVKCazNe7vX/v6HhbKkfKqpD0imHPaoU3mK\n7IVDzsyryC6j1Dq2/qBOZrtvei2kZK4NOfWGcoCmlIde3/R3/Hv7J962Y5nsko8qpVECGxSavvZC\ngSm7vDQlaUYKf3hT6haf3CGlFCGlnh3Zx0cyrIRCmMURCaX7exnvbNnd7yNZv5wBLkx7n33EXkq7\nMQJsxKy0sJS04o6ISiUhI8D6VSsPD90Veh0v7HgVBsnYOBwPHjGWJN9EMjWVDNvU22M+IcOoQhxh\nOTmqFYyy3WtCj4/bn7jTJReObB9Ci1lTyCvrSXobbrkqSmIvkZ0MPmRSXEGnx1I2PIxl6Li2ML5B\nc1I9WHo7Kfuqjvtw+kUJbFB2HLEcH7ZHlY3J8jHla6BMwJFb+gE/viZSbogc0sIYCZYSWzJxVHZ5\neTfn7D8JKgulB0eVMdiatkGeb0qZr/SgElShlHb768DNP94ZcuOXS/0SZccK2UylT1efpR8yqhDb\nZOfSLU1mqE97aJHcwVIuUB0ZrbwkSTXOgc1QPCo1bskF4lrE5CXVXV3ajFKNZZVe0nrIbcYIsGaG\n6eeb5j9sCjDLwWIGqTbVyeAVI30llEwaybAYp4GcfaMjLYyFsO4M61OtKjQmq2u0gTUr33Q1e/OU\nbINWqWGsqbUsi16IbarFL111ht8kdcKZYxxMNuLeV65W4XJvbyKa3pjUcUSUbchOI/I0lMCGiEHN\nP8+UX6TfFcVYVitu0deiHFYbGWCDrUpKmaxaeRiLawWVCky5fHgv4wbfmNIPPlaNgWzRTVb9tZbA\nslWXHU1u+m2j730m0IljyPvsq3Hv9u4hWcoLq1NskF7NLPPnBtWQi2++3JEoxltllt76eDlHfytU\n1mVbsCxKWRJnWNNul9tQbQlTusv/u/1+VM/DGFWIo95P3OrMh3vub9mto1THaOymiGZxKScbG6UH\nw+OlsCnXOec3IHmQfoaDsn0uwt/H1b8vAJpq8OJ+89V1Jp31q93lky15TtUiJbAJmn7FJy99n79j\nLzwbSyoeHz0T2hTh8DVNnyHJju+HWV2ZZ40wuz0ZFdk6wLgEY1ltspK+9LLQxyknmyQ3I37k9lcD\n2sDmS+ZpLbeTfS1uXK0lJjmaffT7mCi33co+IgMvPi7He7V2x7e3mXEOKGWdwaJJqR+jvfXL2Pc9\neUdd6UM94UwgwK5Qy7CtGmMfnBj7OZJ9a2zxqpbPkkOj9OPQ30hHgAXHHlZqDo1m9Z7fzAFWervW\n2DNu1Q3QBnYFw8lUXCBpQdkj2QXQLdmNet+/KyfCsNw/8/iDMyGY2A6k5Mf0HT4StMFzTgPawC4S\nVQDKzoSveLGmdc7YNPxI6mmz0yFmS2b2YWRJk5uxZavUwFZ9OyiUMpZe6LkgtzqGXT/wgkAJ7CNe\nW66/4gPPvwXJEm12nlzlFpfVTIrXb79VZqmbYna1T5adt1ey18t9qlPfczoTdqMN7Go0Ynlh7DlW\nLRiVjuzciYU4f2KWfSuLX/ZBxK17u9oGNiVrH3gOUAK7mvyBTyPWmvacZJmRru36QW9aM+dPorRD\nktumNM3GNDJY2FgonP6+t0cJDGhmnBHKODVG0+iI7MJ8i0vkZL5bbcp5ZZlD397Obkn1HY0z9o5s\nmGuUwIBmc68U+tqy5fUHXqr6xDtKuZNWUymnr+43uyVV9iWfeUpQAgN6lMYjB31d25NR8JaF+QpX\n6ZHTMd65e58bG+csvSJbpxK+JQIM6FGtRWRs1lKU5Kje0DIxeEDtHUzoQ19FgAE94suQMqcG36/V\njHf+nHJM5WZ035r58MwzjQADOpUyjLFZjmTzTO/QMeuYlvrWd8TYY08zAgzop/+c58vlguW+l2d0\n9iu9r3Hqxfjljz3TCDBgCCPT7yE7BXDsjPHCI/eXIMA25kIEBj322oFxe+aOS29P6a+dO5mLR4wD\nA/B0TVNPTWeZ80V/4Ukbtj4CDMDTWTKAqXUXRIABwBdSyhcCDAAqCLY1EWAAoE0Y//AJc1dGN3oA\n2LaxTu34CEpgALBtjbcFwAoogQHAm6QoxkVyWQQYAMAlqhABAC4RYAAAlwgwAIBLBBgAwCUCDADg\nEgEGAHCJAAMAuESAAQBcIsAAAC4RYAAAlwgwAIBLBBgAwCUCDADgEgEGAHCJAAMAuESAAQBcIsAA\nAC4RYAAAlwgwAIBLBBgAwCUCDADgEgEGAHCJAAMAuESAAQBcIsAAAC4RYAAAlwgwAIBLBBgAwCUC\nDADgEgEGAHCJAAMAuESAAQBcIsAAAC4RYAAAlwgwAIBLBBgAwCUCDADgEgEGAHCJAAMAuESAAQBc\nIsAAAC4RYAAAlwgwAIBLBBgAwCUCDADgEgEGAHCJAAMAuESAAQBcIsAAAC4RYAAAlwgwAIBLBBgA\nwCUCDADgEgEGAHCJAAMAuESAAQBcIsAAAC4RYAAAlwgwAIBLBNjHvF4v5amD8XEsq+Moh6fO3C6c\nRT9wHNa5/vfpDUDq9Xrt+579e9u245/x4/CodJTjf3KU74ToOgMlsM84zuaOc/r4bc51zYXuowyn\nlCO+7ztf2+kogX3Gvu9HDiXnun6Kl36zY019R3l7L23Dke4jjj4E2IdxZj+B/SgrVYtwhAN3DQLs\nw/il9gQc5afhiF+Dn3gfo1QThbM/acxPHsf6OMpPo1f/Uqqei70JAHCJXogAAJcIMACAS3TiAIAv\ncecLmlfWRxsYABRHH1+/JbCjChHA05VmS2EWlcVRhQgAX0KBi+BygRIYAGxblF7x3xTCVkaAAXg0\nIsovAgzAo9FTwy8CDAC27b3dK/xNvK2MThwA8IXKRF8ogQF4ulIxi+LX4hjIDABwiRIYAMAl2sAA\nKybKA5ZCFSJQVxoqxNcH+CBKYABQxCS/K6MEBlSES9h/21/hwV/bn+MPvkF3VZ2hg0P/cXTiAHrE\nYYb7keklw4o5qD6OAAOAN2+9db7/S/6WS+J6BBjQI1Qh4sZKVYRk2CJoAwMq4ivUUXMYpxffoNUo\niWI5WOHl1UWZL/HjCDCgjm70LlgKQ9VDdqzEeFzJsM+iChGoy16euGYtJUmvPddkdSym5Fzpqdf3\nf1gKJTAA7iXdLvLLvP9z3/d6R/nsC3Pr5EL6EQxkBuCbJb02kUYf6XyRFhOJvTFUIc70+vbpDQGe\nyJIG1WVkreMs8srAtWIQJbA5khPxqx2Yn1dP8tN7jeN+IXunwYZ1vv9zj6sKe1vCiq1rL9px+hFg\nE3BqPpnys5qjvyY9hLJP7eqzre+evBcXim4E2Ez/bH8ff/ze/v3sluAC1fofLkyXOWMvZ+aOsr1Q\n/oI5o5iIjV6I48KpGdLrEDKMPXxLSXolUyMm83RwDpynadjWz6tsi+U7HBbWcxzlzCSKUXfHtE//\n+2uvlN3Oi7dhHCUwPE7rDTKUH9RbYVZfOWEHznB2J4hj7XrPxmq3xgV7amQ3yWOFAb0Q8SBKH9Hs\nU8mDchl9Tvrw7IKXsHs46YKbrPZVKH4VXy66Mv789HlfbfbtrnRel8trEGCjwskXt3tRf7iabD7J\n+IkXq/6gnn5HFYZhfNy+78d3NvwxHZN6TOSvzLigwclDcTZLsUm2WsmgylYJymez7WH6mcBtf7u1\n9o+oJofc7XKSKn1tTV0Wm47yrPE5yU5boSzYhwCbgwvQyuxlJhlReu+MsEzyeNO9m/kBNKIpwKx9\nN3K7PftG3ePDWo/s3D4Xt5mcmk4cc8TN++5OgntrqvGTUbSd3Bcj2x+E3h92oaz8qmWYvGaX+gRm\nv8XxG2XX05RhTReKUp+LsGHmt93CS+iFiAyPJ8GNnddeFfcznJJ58RaGHL3499DbpIIOz+Rqhumq\nCWSZ/1dZ+aH19ScNN4w7/Xs81gc6cQBvLHWG1VcpXpHmjTuTnA5ttS0sia+/pS2W3Qiz1+z6TInZ\n2sWkp+K37Gp/kqx99xYn2u+O1dP6qlyDABv1emdZZs3rFyziWr6k14aeYdlRruucA9WRbYtLMiyO\nk6TL34S5oN4Nr+9N9RKRzcLHogpxSHbk0FYY8Vp6uetfQC5MvCIrBTK9F2JTVlX7fZwnhLG7djhZ\nv1fqItg6EaL9fZuO11vo2gZBZ5vfVvkR9AkEWEbaa7a989hRK50sEM81FQ8ao9/HasLlu9rmIaPR\n2FVkE1k4HhgTu8K6S69DaTKnrMEGs+R9v/p3FH6Sxu+VbNw6pXCPuHT+qJ5J6fj83CyI2Wl8k2kS\nswvbY5JDZnfsvY4SWGnwVvi1IevcjpdUm9CyPeyzOkYIda/nNv2qS+JO8PFH1Sc2tK82Ybrds+2W\n0KUlkkS8zZFqQoB9Mf4OkqV+GU5JhpXSK1lYH3qibAYUfb0QjbmSPTrGcWP2d7FIPqZcc1+G3ek0\nS0ZxdQxnVtZ5qmyGyfLcnQ6WHVWI6VmYzZsQM62tVnp6HQscK5cVifGGHeuJN+Oa89V7v+ozhHMg\n+wv61/anNDdH+Ge1gq66q/XslC1YlhMm+Tg3O9zJcLHxpqPp6VUclPa+tdV5QB6FEljxfihSUlqK\nS2DVG4DFtwoz3nhFr6K84MBNbFP5lI4SWCldZCooLS5JYikjlJvm7Ag6Cn++DtwZ3uI5PJhbsqn4\nZSnSybzMFKHkW7w/pbzFYw/u0wPMnl6HOD/i13bcwVLGUjbASlWUBJhRU4bpsx0my1QnzUvWYF+5\nUhZPXluawirOTstkjPdWrIo/nk0ebKw8TKcTbOwrGJ7Kvqslw558ZKlCnKDv/svyVXHfxTNOyqZA\nktHu/TbTceVeaQH5YLafYWmajFKzvDJTcFa2F2tWkmFK5eRjO7tqvYXFIx27SHYJkXEli1/yrS29\nItMeKI88oLFHD2RuLX7FSxovB/9sf9tXHm9S+COJje4U0QdZ963ThfgwTZ9WI/t2YXxr9gyxnDbG\n9IoHU2fFH+SZF7umc3uw5+HbqvQ1tL8ks/wjD2ji0QE2Tg7FN55VR7C11ltOKQNl31fpgXa8qevi\nV5Jh2akLS8OQB9+0dD6Ex7PlJ+US+d/2l3ELuSt0YhchUZ1TQ06KkW1Iy/5Tvnt2e+784/F8D61V\nOHSUwDZDN4rqD7RSJ474WT0wWo+aZciasUXH7wlj/DFe6n8Rs/S2GBkh1LSFpQWSZfweuBH2jhsJ\ny5jopmJWvSOGuhL6zWdRAmtLr47l9ZeXgkp5l5NOXMswINffmdKv7DPos7Q0reoocsnOismDspRJ\nemUZ9/5g1XrrK43phQQBNp+8WJRqCzvq5ZRi3xmtWdU6Fnf2snix7KS9yT8tRXC9us9eH7gVOtyH\nB5VG2XscuHGv9gJQdYWD776Xn6Vq0YIAa06RZBrDrKYVJu1hv7d/j//iZ+WrSu+ov7tsThssUN5D\nvA+zlYf2Mo1Mr3h8sYyrbD/4UqplH0/KYXowI7FHbWPx36WFE9kvW3f2vKL/lDfd+FHyjQDrZw+q\nOJD6ekPEMVPNSL0yMN4YbLUeE6e+tUyvcIfMpgw73Ls3aR/9Ql96LhtjpYVD3sjsKW1SvFVNkckB\nTjAO7AOO/FDm5jiUCl6l1ca9P+y9/K+c2mNBpWFbcbFpM/flS9ZWGlhGz8ArlXrNtA66kq+Nn+2b\nJuMVFf7i9Tzxq9jl0b0Qt4GZOCRlgqi+n8bVSe7l8sr09sk2yNlDnnkm6IemNAVUqct1/EJl1ii9\nA+ERn0o3SGXa+2cexCZfEwUYF/7+ozTzkzKiS+8uqHTHr28MR/kbVYhfLAkRL5M0Ten5Z2+NSOr6\n7NsmF0suzcm7k16xOAz0pq+DpeEzWadS9Sd7Z1BEO0lxOH+hArB1PLJeH/i25PstpKuoPMx6eoBl\nAyPLXvZSHqzGWPIF62uviqcLqb477fwxvQ+FXdzTvdoRUXnHaqcPjLN0/2vND/tQs+ryaZ8Ovq0R\n2sDeqviS1qn4wVjr7FBvZ+qMlvbs/ISlNrPkjOcLIMU3Oqn29EteEi9jSRd5I+aOmfItdx3DXK/3\nuTOUukT7Nzy++CQ9D0srMTZvPwT74ktTrsj8UIZ5KekVN0QljVLVGe6rGfbw3hlGluNeCphSH0Ll\nhZaksQShfYJ8xMLhjneTMjxrE7lSeklYILO2rpvjKGvmKAdPr0IMptekGe8QVmpIqxatqs8ywMui\netCNbVdZchC0ZZOqZbKzO/ejW1/tit5GHg9T+3oXxkt8I8DetI79LN1qJHurlOxp19rERSxNpxzu\nalRY6htlT5CwQDjZlKlAkgeVUh0/zPvs739bdmJ1sb4jwfDzVrSB5YUTKO5n39HLYy9PCZpdvu9d\ndFSaV8WHKTRLZKdu2kS0xKOP5YOSTKBwdMJbH4Ww7jtwoiTs4ZfIrSy9pCNbv5hy92KUwDTJKDF7\n14l4DTK95Ev0lRwjnVtvLYZWpV++stuFvp5sz/u4R6Jy25SwAa39QbhQTmesp9uj/07dDEJRIsA6\nyTFhTenS2vQlJ0js2GZYrNPAoN+vkvTq81O5Uuu8bnxQWWzkuCSd7Fc5KRdDFWI/OSNU/KxSGRhK\nVNW3qNYobrWekLhAqDDMdkRsLcAl1c5UEk6XdF7Xp+hNDkd2eeWFs7Zz4mrvhBKYJpwrSookZaND\nX5ZUe4LIlVuaxDjjm3TsLmM5SV+yYzPie6kYV4uDLN/kZ4KPmieT5RMnDTeOq7Xp3CHRvF/RcW+U\n8LexJ71xYfuqkhVyiFs1deLYxPjiuSmVPQPltIoc5Q76t1vuUuPVgGNxGaoQK7Kl+Fg8B7xFXCs4\n5bYm1ByeRM5Jb6THWOts9LJDLOPAZpGTDOjZo3QqlivEBR+HWTIAAA8FSURBVAiwumqGxTkU3xul\nYzr5LVdyit+9ddZ8vlFnK03vJGOM4FmZ/ZvCd2odVCF2OkIlW/wqZYxxOuBs9jRlGOk1zljc6cin\nkUq/0g8pDjSeiU4co2QnQGPFoDJ9ohRfoao9SizvDl08HqtU49fRQYMuhQhKN+KBHSWwTieVwJQ7\nUm4tPUo4rONKvSf0EDLO/DtygN7uhciBdih/h2gOZTvawIbELV7d7N1Aqm3IAdNHjcvubZleyrCt\nUtoNHhqOrGvFm2rynW3HLuuUzDJV1VGzpx+a7Ncg7uLIkZ3I8pPZdFMMDsp6fu6xcv7ReSs9hwej\nBThDmlACu5o+uUZ4tnoex30j6Ul/NstlJdtrtGkN2ddyRTuJ3NWWnvSzJFMJ0w7WhxJYP3t+KLM9\nleajah3iWuqCD0cYYHSZanH5pB3+lZHZp85837siwIZYMsw+db3eg0PfgASHdaJrmtzjd2GijVNl\n6/G+ngqPE2AeEGBDkktbkk9NEyRW06tUuUSl03lapxqa8l7yfphnvN2T/TR6ZZ/9/uOMHV56a26V\n0ocAG9U0jMMy9suYXsrCmEUvFW1T97+cgDF5O471LEoxaLsqwOINoBNHNwJsDr31YnD81pXlAATK\nZBxnZBgBdg29+PW1zLHAmbWIWRziVvRCnEAPmGN4x5RG44mz12NE64S8QFC6GpBeHSiBjVJmKazO\naiifVdZf6sHIETyDZS7EI8PmlsC28m1cONCz6FWI21X9Ka4cf3ZXlMCmkQETBnXFY+zP6L3GF+BO\nzpi8A1Lp3sqXDcnimI5jMt8hDCXGLMrljCsdkEWAneu8YGPM8qmqe/WMNrDsPeM5vtP9zJnyXt56\nndwFEdMRYKsLX6T4Vin04LjM9Z019ncXv/sDvUSSsdu9oPlkSLUKcVY5iXFg18sOAtvoWHEj9Ab0\njgAbpWfYxIo+Zty43qdmzANgQYCNau1GD18o+wLLIsAm4Hf6EzBqB1gNATYH08MAwMUIsGloowKA\nKxFgAACXGAcGAHCJAAMAuESAaZpuVgkAuBJtYBmMzweA9RFgbxjRBQBeEGA/LHNqbGQYAKyBG1pm\ncHOvO6FCGLgrSmBfSnPylm5cwn5zgRlSgBujF6JGue2Wiw6Kr9fr2M6X8OlNu0LyMeNboshnAbhD\nFWJRnF5xsSw8/notXX4NF+jslXrxjR+XvZvX8Ue4odftd8JdMbEyDpTA6pJKRaWFTBZ0rvmZb3+X\neONLr7pZWS0peCWP3OMzPkdyvG5zlqIPAWb1e/s3lL1khinfolO/XbKSMDyV/XH6e/v3n+3vUgZn\nP4XHC0TYYJlev7Y/8S2Vw/LuPuPT6MeIw/dMVCFW6D0Ss1+b4yVxDeSx2PTqDv1LK0cFyCa9uA7N\nsjbvNTYyumL3qFGMj1RyTP1+urdfZslT0TJ+PyD6EGBf9n0/viRHGSU8Hv5p7FsfFgt/KD1BOsiY\nSd7o+BrLxfTNkGlnfKFTSclMDzZHbt/wKbd+J8MejCrEjOolW1lAPhXy4KRajmy5qrrMQf5Olzkd\nVzneo6JGbxVzrXr59ngEf7psFBaIH6dV7FEogeXFV3xZJouXPKmkdZJ/tr+V7bz3CO6kmHX8U/at\nlxkgr4af+plv2ZJs+bv7jVwUaLKf9lGlsXvU8Hd40DHWjXznZYbJJAhPDe7wpg6HeqYmrV/V9DrW\n5uWEkf04QoAlnemTf8YfUN/bF++K0sYoifvf9lf8MY+/lc1WPu8Hj7teAtO/D15O10HJD47n5Bkl\nsG0r16EdZS9LjeJlZZds43xWdbP7fqK6+2H7a/uTrSEsVRvafyJctitGfl3FGTbyLmse97DF9Ow4\nJGMMtrvHGAGmtQCF7ht91YNLVSpmOyLeuKkgrkmTlYdNnTiyffGPPz54fcwmU3JA7Z1T4heWdo6j\nMNhrJbObUX7XOjpqHZ4eYINX8DjbkvrD89JL73NR2oylypHXUFqD4gzr6IJoL9ZM1/3WyquU9Ere\n8SNXw3AcX4VaxPDg6/2fN2AfCKEPkrtrht32g1mUDnn2Om4PpCQ5kn+O7/Cw2fp2Kr1LlI7ySob5\nagMLWn+jlNrMpGyz2RmMTXGWT5rdVGXcd3DZh80qNYMliRX/86d20XDRX/Oszh7Q0qYu1V57med2\noz+j9kyZ5GIWPb0k2Wcy21G+tHz18ZX1daeu5lZwWef7fd/1q/CUjuP6x/nsSIOfkFbrBndb8Uvu\nsQU735c2xvQzZaAY+srpXdnpHlqFqPTa2N4r01oLXuuoFrzkU/LjJ6ty9Dsu+dY1NXqV+n181qz+\n8Y4OYiz++MleKFUtfj37Pr+ai74qH8kMvRJyW/Lr/9AAC2SvDXujUbaRSWmU2uadAXoUZbe8KV9D\nhnkseG21Rp3N1pK02vQcSiQrm5oMEpDr6diMT13IlPh5iYrE/GI3nZIq2S0dB9jpnllug66hjHwy\nXrIt1W7SrEFgxjFbsb4qx5iXU6WaXofkoq9EmtIL8XBxZ/rsh1I23pjE1UJndRjZR1TzOGkS2/Ti\nWljm0x+zaUxefSeYp2jRB9t9fLckntsGtl1bvNCbMaaTbV0jH/bijR/R3Y+8O73km57B0s8iLFAt\nnHVUkK5WHt0a2zh9nL4RS6vn6/tOFN2rCiv8WrK0BrHkIh4dYNJ45Vti/zawUW9rsy8sP0tThu0R\n+6s+K/vtyt4/ZbNdxFfu16ArVZyGeIv/PiwYUVmtPQuarrjLXqlLLNtp3F3693zNq8DT28Cq/TVK\n6ZWto1uks4Pseei0KatJtuawdFE2XqyVIc8hBo5HTm0esBe/AllBWloyrmMs9V75bB/6oDtU9F4e\ny2oaIJF5+QNGcz80wOLW4NLFPRtdeteGy6LifoOOJ4ovweHqfFyaw/+Vl+jTdnxwCLPR4LjsZI7j\ndT5schGXgXTXi7XSvWJr+dRJb8zs48b1rFMr89wqxNIxkAUyOTRYH0p1DXtSdmTqOieokX0obscV\nOVRCxrWRx9+LFE22WhmrdW3JpzusMJ/vVh7klL1V2M8ayoutLG6dMn5q/fGw2lCv2HpYP36qxx5a\nAjskE4jpPeB1l9XRJWXHpuDs6LuIYJ2ySInsMd/0WuWF7moO+4fx9r7wU7QxA5aXr1Sc6vDcElgQ\njp9yM5TSgyWntgDHJ9xgT40pL7yBpAKtKj4Eq/VzGexasufM2rZW2d4H+mQc8TKtX8J0fPSnJ6Go\ndg6M7YUSZ/xs/J98F/2jrhntvuN3luQ0jZuyjJ04ktavC/Zq91fLMvrN41mRrUUs9T+sBlWyzMd3\nyPHpLOHU3eVynerQQBsLJRc2r1Yf6pQuvEDFqWngmnhkE0+VXltd+G0U3TKnx0aAxbLflr7BzheP\nbC1RupxY5srydW7IiSqqY3tLy8guDB/fFfaOiNkbddoD7OOfNChdvrOtWck45WqYlda5iSTwFWD2\n9Eperrxk2fTaCLBEKRJa5924eK/KiUWyJcjWSkJ354ae6NVbMyeLJQuv0w5knF7k9gEW7IUE2sTF\nt9SVceUAq769ctKbpjY2b88650ZAgBV119Fdv0tlP5TkDimt82M5vXPKljtqljmikut+dvmP7w09\nw1rvyVl6+cc/ZlC6fDeVrsZbbj4eYNtAhjUFmD7T8TonRowAq+jpAfW5ANsMdYOyXjGcuDcIsC13\noa82C1n6bqywN7p/VOkBtlRrX9AdYHO5CLDtvS7xZX5V/Nrs4LB1zocsAsykYcq1D+3PpgxLlH55\n+T03kgyrTkvhJcA2WzVpUz+ONdNra6lC7FD6pKXRvh/RF2DBSIB5QYA1yIbEtkwXvsHCYto84PzE\nyO4N/SKeRN1qHRETstPK8Uf11ioxLx9Q6WJnkW/rWuzzZhkzbKQKcSPAHkK/m8kKEyEOVmEv9fNz\nnDHD9D4dh2X3huxe3zeQec0P2PqbLNv/UOnT0bdVF6tmmL6P7KU3LzskRoBZ6el1WCHDtrGyVN/s\nMivLdusIf2drF2UGLLtDZPf6JI+TB7NkQXydz9t6w5RSgGUeX+Yz6pRJIJVO8NnBBvn1hyWd7JDY\no6eSuquRE9HjSayTU3orV3Nfxa+54gvlOjGWNNBWx3hlF1hzFgmjZA9UP4vrD9uKqaTwFKXLsT5h\n0goX8T5N82OtfPurtwJiaZnyy9f9YGbVCb32li4bMdfFr40SGB4lKY0lX1qP3+Hw8zzc9kXeGyXI\nViomn3pkRuDzZDvKth4tywwdK5Pnp6xdzJY++xrPXCDA8ETJtcD1KII4w7bvWy1nE8gSS/EyS01V\nrg+zDZQlbnC9jnUMc/beHCgtdIKub/1eiOjgbvYBqfOOvbW7Gi748Wfd1nLBj9ZKBljnGAPPu4I2\nsB4Tb2KCzyqNpso+uyy9jWQvSJY5fzMnSLb8FV21X+YWHS8fttse/XdvlMDaKDdeObA/HVFmF1x5\nhG9Vdx9CRwMB9d8W+/u9apPH70EpgZXmLslWIbreJwRYsxvUOOGg36Nkwfltkbh3RFUl00VmA0zp\nt+l3/HJAgPXw0mAAheUOW2QYVlYthFlqwEtNoS7OeXoh9nBxaAE8hOwr39F4m/wuX6oPagmdOADA\nJdmfpXlsnJinJrZ+JyYCDMBtrX8JHpSOaKwtYHzKC6oQAdyNnNdxu8X1OqvU5TI8pU9A45qDWk7g\nJPRCvJ9qkYujWeKxSydViEBmjqXVJgOERTJv/UMG8z4ZJTA8Gj/Y70SfHtD7zOuQKIHh0R4+59Cd\nVCe35VjeD5048HSyJzG5dW8uRjjBggADvnBRA3yhChEA4BIBBgBwiQADcCulfqX0QrwfAgzAHeix\ndPMZpZ6K3jgA7iMZ2LeL6OKKdyeUwADchz65Lel1MwQYgFsppRTpdT9UIQIAXKIEBgBwiQADALhE\ngAEAXCLAAAAuEWAAAJcIMACASwQYAMAlAgwA4BIBBgBwiQADALhEgAEAXCLAAAAuEWAAAJcIMACA\nSwQYAMAlAgwA4BIBBgBwiQADALhEgAEAXCLAAAAuEWAAAJcIMACASwQYAMAlAgwA4BIBBgBwiQAD\nALhEgAEAXCLAAAAuEWAAAJcIMACASwQYAMAlAgwA4BIBBgBwiQADALhEgAEAXCLAAAAuEWAAAJcI\nMACASwQYAMCl/wdXFiiQH4dOswAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 15, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "tlist = linspace(0,1,6);\n", "for i=1:length(tlist)\n", "    t=tlist(i);\n", "    Xt = (1-t)*X0(:,I) + t*X1(:,J);\n", "    subplot(2,3,i);\n", "    hold on;\n", "    for i=1:length(gammaij)\n", "        myplot(Xt(1,i), Xt(2,i), gammaij(i)*length(gammaij)*6, [t 0 1-t]);\n", "    end\n", "    title(['t=' num2str(t,2)]);\n", "    axis([min(X1(1,:)) max(X1(1,:)) min(X1(2,:)) max(X1(2,:))]); axis off;\n", "end"]}, {"source": ["Optimal Assignement\n", "-------------------\n", "In the case where the weights $p_{0,i}=1/n, p_{1,i}=1/n$ (where $n_0=n_1=n$) are\n", "constants, one can show that the optimal transport coupling is actually a\n", "permutation matrix. This properties comes from the fact that\n", "the extremal point of the polytope $\\Cc$ are permutation matrices.\n", "\n", "\n", "This means that there exists an optimal permutation $ \\si^\\star \\in \\Sigma_n $ such\n", "that\n", "$$ \\ga^\\star_{i,j} = \\choice{\n", "      1 \\qifq j=\\si^\\star(i), \\\\\n", "      0 \\quad\\text{otherwise}.\n", "  } $$\n", "where $\\Si_n$ is the set of permutation (bijections) of\n", "$\\{1,\\ldots,n\\}$.\n", "\n", "\n", "This permutation thus solves the so-called optimal assignement problem\n", "$$ \\si^\\star \\in \\uargmin{\\si \\in \\Sigma_n}$\n", "      \\sum_{i} C_{i,\\si(j)}. $$\n", "\n", "\n", "Same number of points."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 16, "cell_type": "code", "language": "python", "metadata": {}, "input": ["n0 = 40;\n", "n1 = n0;"]}, {"source": ["Compute points clouds."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 17, "cell_type": "code", "language": "python", "metadata": {}, "input": ["X0 = randn(2,n0)*.3;\n", "X1 = [gauss(n1/2,.5, [0 1.6]) gauss(n1/4,.3, [-1 -1]) gauss(n1/4,.3, [1 -1])];"]}, {"source": ["Constant distributions."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 18, "cell_type": "code", "language": "python", "metadata": {}, "input": ["p0 = ones(n0,1)/n0;\n", "p1 = ones(n1,1)/n1;"]}, {"source": ["Compute the weight matrix $ (C_{i,j})_{i,j}. $"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 19, "cell_type": "code", "language": "python", "metadata": {}, "input": ["C = repmat( sum(X0.^2)', [1 n1] ) + ...\n", "    repmat( sum(X1.^2), [n0 1] ) - 2*X0'*X1;"]}, {"source": ["Display the coulds."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAULElEQVR4nO3d23Li\nuhYFUHxq//8v6zykQ4wlC1/ky7LGqFQXAQJOGjRZkiwNKaUXAETzv6sPAAC2EGAAhCTAAAhJgAEQ\nkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAA\nhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABC+u/qA4DHGobh50JK6dojgUcavLWguXd0TXi7\nQUO6EKGxufSq3wSspQsRWhpH1LjaElzQnAoMDpFmvlWEQSsCDNoz0gUnEGBwEqkGbQkwOImuQ2hL\ngEF7sgpOIMDgEMPMt04Fg1YEGLQ0zqdh9AU0J8CgsUqNpfyChiwlBUexFiIcSoABEJIuRABCEmAA\nhCTAAAhJgAEQkgCDeIZhsKo9mIUIYRRDy1uYbgkwiKFecnkj0yE7MkMANnqGnDEwiMRGz/AmwCAM\nvYQwJsAACEmAARCSAIMwDHPBmACDSGz0DG8CDAKw0TPkBBjEYKNnmLASBwTzcVKz9y8dE2AAhKQL\nEYCQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCS\nAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARDSf1cfAFENw5BfmVI6/0iAPqnA2KKYXpXr\nAZoTYKw2Tqn0+1W8FeA4Aozt0sxlGQacQICxUT7YZfgLOJMAYx3VFXATAgyAkAQY65goD9yEAGOj\nvCfxfY2QA07gRGZWSyn9jIQZDQMupAKjMeUXcI5Bc8NmkxmJXkvAmQQYACHpQgQgJAEGQEgCDICQ\nBBgAIQkwAEISYACEJMAACMlSUnBT7/PEnawJRU5khtspbrrmrQoTAgxupL5fqHcrjOlCJLxxo/+Y\nJn78a7x/vWHwiRP+CDACy+uVn2uCtvJ/g16f1yc710CJWYhEVeltq3fERRQykOFgKjDCi9Xblofr\nquO87e8F51OBEVKlt+3OiqXhqnpResGbAONpbtvA7+nzfFqXKLSgCxHO8DFVcnz96A4ppZ+7DTP3\nqT/s3+Or0uiDCgxOtbDPcxh9/d35M5nm6rbnzWGBIgHG09y58Z7E1ZAVWymlufqpkl7p96t4KzyV\nLkRC2tDbditD6fJc6tS7BNPn5fvPw4RWBBiXabVYbbHVv2fbPaw/K3l5er2viZLisJMA4wKTUmPb\n8hk/9w83i+FjTCu75iWBYDFjYJyt7dSD8aBR+rX94A4zOao0k16TW19Gs2CeCoxT1WeTbx62uWdo\nVeShtK3wGkq9iNAJFRjXiLWCRhOr+0g3Pcs7BcOFOqwlwDjP3PpPxfs8Wz7xffUj/ObT8CqfNAaP\nJ8DgFtpmj/KLHggwOM+4bPpqc6RJLzohwDjPkua7n8a32OOXRrf+u2bmD9JJXytUCDCuMVmKop/G\nuJ7QC4eyxqOJ1pGiWxab4WyV5rWrV+PXmPlae+U3m4JIVwQYFwi3fMat1CdzyjD64URmLjBZBSp0\nU5uHcehfBwIRYFwmekNfWRMr+q8GIZjEAQ2cOY2iMplT/yFdEWCwxWQe4PtyfodjD2PmMvRAgMF2\nxe24znjeUYGVz7xXftEJAQYhzaWU9KIfJnFAVE+azAkbCDCITW7RLQHGM50zlz3fTzLETArnrvEM\nTljhUc5c4yPfXTrETApLefEYJnHwHJUziw9/6pjpddOjhGVUYDxEXg+9TgmVuXS84Tur+Cf6d9P7\n+vsdNswxBsbTpM/LRxdff+tifCbZz7f3zIM7HhOspwuRJ6gs0J6y+xx6AAuvvxuRRkQCDPaadM3Z\nXhLOIcCgmUnv5TnFXxN3Pz4oEWA8wZKhpqOHoyr7I79er+HXocew0C0OAnYTYDxKZYeRk80974UZ\nNlkC+H3BFESCEmA8xHiXrPHX5NZz5Bl1k1GxfBl7iEuAQRtDdmHiJtVNMctTSsovwnEeGM8xWZ19\nfOU58qUR/x3D6MIdih5ZxTNYiQMaqPcKfvQfvq/01oN9dCFCA5U0Ktdk0gt2E2DQRp5JKUuvO/Qf\nwmMYA4NmJusizg2JAU0IMDhQseTSfwhN6EKE9mpDYtILGjELEQ504Zx+eDwBBkBIuhABCEmAARCS\nAAMgJAEGQEjOAyOY4qqD5iJBh1RgRDK3Zu5NdjoGziTACGOcUilbaVCGQW8EGPGkmctAV4yBcYG8\nWvo6iPX+kfx+N9klEjiZAONsxb6+YWiwKEyTB+GGNnzioQe6EDlVZaTKIBZFc594zj8S7kaAcZ5x\nN2A6YCKGT+XP4xMPFQKM6y2Jnb+9IrObNGNPdfQnHqIzBsbZ9ldJdonsnGk7/BBghJFS8ol7s7jz\nIGIcJVcQYETy0+baJXKVyvIl/nSEJsA427D7M7Vmd7OfP9w70GQYoZnEwXmKEzGG0bca0+Ym8yDe\nl/M73Nbdj4/rCDCuMXxG10t6Ham4fMn9+cRDnS5ETlUcxHppifgmr8O8ZhBgXEDTw0I+8VAhwOD5\n9k+cuZa4osgYGPTIzAgeQIDBk43nQQyfFyZ3gHAEGHRkUnhJL0JzGiN0wTwInkeAARCSLkQAQhJg\nAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAE\nGAAhCTAAQhJgAIQkwAAISYABENJ/Vx8APNkwDONvU0rva1JKVxwRPMfgXQRHmETXHG9A2EwXIrS3\nML1W3ROYEGDQ2LJM+iu8ZBhsI8DgKjoPYRcBBoeqp5QMg+0EGBwnfV4WV9CSAINzGOiCxpwHBpeQ\nZ7CXCgyOM7xew6jnUGhBSwIMjvY9t8ykhw0EGNyCDIO1BBg0tnh1qGRqIuwhwKC9aobluSXDYAsB\nBodIv15W7IVjCDA4lvSCgwgwuJzpG7CFAIMzjOqwSVxJL9jIShxwvkJo6WmEtezIDOdZcrKXtyQs\nJMDgbGIMmtCFCJezWCJsoQLjgfIS5yav8/naa3x4/+5zk2OG2zILkacphsQwDJcvNlg9ALUXrCbA\neJQsJA4vYoZfa34oZRdeowxTeMEiuhB5js8Umbyw2/fLFUNr7vFHd06j4xlf/vnWDHtYSgXGI51R\neM1dv6YaG8fYa64j8fLOT7gnAUYnWkZaVurt2RilEk5/DyjDICfAYI9U/HZB3hRHv/I76DyEWQIM\nNpuE0EcOzfclLpmsYbcw+M6JzDzSMDeJo8FDT2NpmLlcvPPkbsVZGxILFjELkUeZmYi4awri+DF/\nHuHIEaniEQ4vcxEhowLjUVJKo3TZGzN5UM1EV3n6eyOmb0CZMTCeplippJTWVjCLy6zlD1se2fo8\nNnEFS6nAeKDWvW31GYP1efAbbnUuMyyiAoOCbOGM18y3ufodZistEQVrCTBYq5402/sAW3V+Qid0\nIULFcclh2UPYSwUGFTunVFQGzEzWgL0EGKy1OXuUXNCSAIOC0rz2YbJe1P7skV6whzEw+GrVLieL\nTmoWXbCfpaRg1o5zmT8KtWHwRoP2VGAw63Nhqld1geBa1EkvOIIxMFioEkKmFMIFVGCwxMIS6udu\n6yZ65Kvdr5L3cyr46ISueaj5jYfZXU6+qrzF5sbYlr8rK6N03to8ni5EWKJ9J2ElexZOHsnG5yQW\nfRFgsFBlI+b1j1XIno8E+pphn8sNp/zykbtuwi0IMKj57Igb8tOZKz+68Bk2/VSTp4bYBBh8sXUw\nSQEExxJg8N2yDGs8CqUPEOpMo4dFxhn2Of40XiwxTZaS2jMVcMfPSj66IMBgs/T776JNUg6rqAaD\nXvRJgMFmQynD/pmp2IoP8hol0NKQ+1zm6v0gLRfLh5szBgZNTOfBV++Z27JXS3afj/AbhsEoGs9m\nJQ7YshrT5zDY39X5j8/c81Uptla9K7+mlPc4T6UCo3fFAFhTu6w6OWyskCsppbV583n/LedEQ1DG\nwOhaabeU4X3TTzYUF9vNb8rvM7m6chh7iqT5Cm/R1poQlwCjX1n8TCfH5/n0c00xxtaE0IpcWf/g\nYzKMJxNgkPu+DuE4UVamy5Yey/cd5p/LQBfdMQYGXzUZWBpG/243/7wqLbojwOD1GVFzq+tumKbx\nen3UTHM//nU9383ZKdV4MgEGuXy++86ZinWTiST5pIzKppqVgJRePJwAo1Nb42dLMbRhcvz8s9cV\n5vQ7D4ynEmDw+mzxJ/VWvm/kdMGLRU/wbQfLzIrUqUSU9OLBzEKEt0oUHRcDczPd1y3R+w6qfdPu\nIRIVGJ3avU3lnoTYPDr1/QfbdVfC3anAiG3NWhh19bmCzSdEFB9wsra9SRlQowIjsLnxp+VzK8Y/\ntGbuw7rFdhdPGKnM1DcpA6asRk9UpWUMXxsa+noNN/MsS5+o+UK6m38p73SeR4AR1ZJtSpq8vLMM\n+8iG+lN8HuSuMFv+u1RS0/udJzEGRnTFdSuOGy5akV7jO646pEaL08PDGQMjpDOb6bl5fWs2vVyn\n0W9XWB9LvPEkKjBYpF3nW16NNSsZ57tV4YFUYIQUYiznZgd5q4OBBlRgRJevWPF166zzDfNnd01W\n771pF9+k7/FOf1v6JcCIKqX026retNH/FHUfr+KwmQWruANdiDzTTdrWNYeRmtaOeepsScr6pA9T\nQriWCozAfhr6tb1b7VafWmRUKdY1CIPP5xp3ru598NlT7eA6TmSmL/c5yXfuSPYfRsPf8V9XYfGm\nTQ8IDQkwOtJq9akdTzp9ovetbZ+6SZX5d2zFW7c+LLSiC5E+pc/LR3WJVZYb/mn3D2r9hQo9MImD\nXlRP8j2kuZ/Zhbl4K7CaAIM/DUNlfhn7MLXRXz9ndpP+Q+5AgMHRTir4DjXMXIYLCTD4c5OJiPcx\n/oNMd/xUfnE1AUYvRq1tm5N8OzGXUtKLy5mFSJ/+ncV0UXTdcLXGmijHSW9UYHQka4gPXKC2WvAB\nDajA6Mu21ad2O3XxqijG/wv+GmxgJQ440H1WrrqV45bRoisCDA5nM62x+tzLzv84rCLAgFMVl1g0\nO58NTOIAzjO3QLDIYgMBBtyCDGMtAQacTVbRhGn0tJSPzxvPqOvzLzbIMFpQgdHGMAzF2WX3X+7v\nKv5iE53+2uwgwGgvTSaY9doir9HLVmFzW7Q89hfmSKbR08CSBtcrbay6u+aulRKL/xe3+uM7D4xW\njIGx3YMLhaDm/keG4UYfVYureY1vgoUEGBtJrxaK7fXGNfI/lhZ8Xzm69VbxMI6xWx0YgRgDY4sN\n6SXwFmuWXq/bT/ZLKUkvNhNgrLY5imRYpv0fJE8D+cBTCTBOJcO+OfDv44/PwwgwuMDndpfD54XJ\nHYAyAcZGa9vXvwZbHfB6vaYRdeB+K38T9ls87DCy/9FgDwHGRqtaL9VEUTFRtsXM3AnCDeWhJca4\nlmn0cKUjugqLkbLziSpBdbcJ+vRDBcZqa1srbds5jkuRyRz991fxDnAaFRgHKg7y+LR+nOIiFw3/\n4PkulIKLCwkwtqisBjQ2aOOucNBHBJ87uBsBxnZzDeU42Cbppfx6GB9QuJAAo72jO7K4RHEXSunF\nhQQYR5FYwKHMQgS++LoLpQ8rXEIFBqygz5D7UIEB31VqLOUXV3EKPbCOXSi5CQHGLh9rNHgtAScS\nYGw0dxazVxQny1+KXoSdEGBsUV+Dw4uK01Reil6HjyfA2OLdasxtaeV1xXGKoTXeIfTvSq/DRzON\nntWK6fWyqhDHm+24/rzsddgJAUZL2g4OsmrHFq/DTjgPjEPYIIqD5LuRvebjyuvw2QQYhzD2QENz\nvdZ0ToDRko+73IoPUs8mwFhtvLTrOLGkF2eq9CJ6KXbCJA62SCm9e3XyxmLPx14npbKHOfRdUYGx\n0Vzr0Da9fq40FE/RwqkcPJUTmdmr1dKuH8sqvq98X+OF2rev586PebV0QoBxC5VpZjLskdYuA231\nMnK6ELk7LdPD5H3Ck2vKK0WldESvNaGZxEFjk9ZH48JYpZAqZtjk9fPz7dyt9EaA0czcFIyXhobX\n6zUzzPmqTr4YhsIwh5cTP3Qh0sZxEwVNLXueev5M5xaag8oMAUZj6fPrx/I2SFv1ePW57+n3X0UW\nXwkwGqjMIVzYDI07hYbPr/wO9EYRRpEA4yRf26BKREmvPvlfp06AcSPFoJJeD6OYohWzEDnJwhwS\nVz0YZqorS0OxigqMBsbr009og3grjnSu/UF4E2A0ZlcLKupRVJy/A3OshUgzlWkaXmbkxie5e/Gw\ngQCjpblV7M4/EsKxFRxrCTAAQjIGBkBIAgyAkJwHxncGJ4AbMgZGjW1wgdvShchSFggHbkWAMWu8\nxvw7urZtkgLQnAADICQBxhdLFl0FOJ8AAyAkAcYXxWEuY1/A5QQYACE5D4wa54EBt6UCo6YSUdIL\nuJYA44tiUEkv4HK6EAEISQUGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQ\nBBgAIQkwAEISYACEJMAACEmAQV+GYahvtA1R2A8MelHMLS0AcQkweL56yaURIChdiNCRNPp606NI\nUAIMHu6dT5M6S9lFdAIM+iXDCE2AARCSAAMgJAEG/TJ5g9AEGDzce5b8JK6kF9H9d/UBAOcphpbz\nwAhKBQbPV4ko6UVcVuKAjkzOWfb2JzQBBkBIuhABCEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEIS\nYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQ\nBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMg\nJAEGQEgCDICQBBgAIQkwAEISYACE9H90tQVuk+dQYAAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 20, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "myplot(X0(1,:), X0(2,:), 10, 'b');\n", "myplot(X1(1,:), X1(2,:), 10, 'r');\n", "axis equal; axis off;"]}, {"source": ["Solve the optimal transport."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 21, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gamma = otransp(C,p0,p1);"]}, {"source": ["Show that $\\ga$ is a binary permutation matrix."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAH2ElEQVR4nO3dwU7b\nQBhGUVL1/V85XVSCqsImgdieO3POtpuQhF79iw/f7vf7GwDU/Lr6BQDAdwgYAEkCBkCSgAGQJGAA\nJAkYAEkCBkCSgAGQJGAAJAkYAEkCBkCSgAGQJGAAJAkYAEkCBkCSgAGQJGAAJAkYAEkCBkCSgAGQ\nJGAAJAkYAEkCBkCSgAGQJGAAJAkYAEkCBkCSgAGQJGAAJAkYAEkCBkDS76tfwFput9vVLwEedb/f\nt/7JN/kIO284n3KBAZAkYAAkCRgASQIGQJKAAZAkYAAkCRgASQIGQJIhM/C5nbWyjTMjcIEBkCRg\nACQJGABJAgZAkoABkCRgACQJGABJAgZAkiEzD9karlqtrsnnzghcYAAkCRgASQIGQJKAAZAkYAAk\nCRgASQIGQJKAAZBkyMxDBh+u2lnDglxgACQJGABJAgZAkoABkCRgACQJGABJAgZAkoABkGTIzFG2\nxsVvB+yLDZZhQS4wAJIEDIAkAQMgScAASBIwAJIEDIAkAQMgyQ6Mo7x8m3XmsAwYnwsMgCQBAyBJ\nwABIEjAAkgQMgCQBAyBJwABIEjAAkgyZ+TD4UniE1wCMwwUGQJKAAZAkYAAkCRgASQIGQJKAAZAk\nYAAkCRgASYbMfLAUhlcZ/M8CzMEFBkCSgAGQJGAAJAkYAEkCBkCSgAGQJGAAJAkYAEmGzDA/o9rz\neWNP4AIDIEnAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJENmfsRCNsFnwZRcYAAkCRgASQIGQJKA\nAZAkYAAkCRgASQIGQJKAAZBkyMyPrLaQNdyGcbjAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJDsw\neIKx1zhs8nCBAZAkYAAkCRgASQIGQJKAAZAkYAAkCRgASQIGQJIhM5D02rWyWXSRCwyAJAEDIEnA\nAEgSMACSBAyAJAEDIEnAAEgSMACSDJnnZJUJT/F7UeQCAyBJwABIEjAAkgQMgCQBAyBJwABIEjAA\nkgQMgCRD5jlZZcIJ/MWAa7nAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJAEDIMmQmQuYfzIHX9dr\nucAASBIwAJIEDIAkAQMgScAASBIwAJIEDIAkOzAuYD2zpq39n+8D3+MCAyBJwABIEjAAkgQMgCQB\nAyBJwABIEjAAkgQMgCRDZjiWp3e+W+3n5WguMACSBAyAJAEDIEnAAEgSMACSBAyAJAEDIEnAAEgy\nZIZjrbbeNdzmNC4wAJIEDIAkAQMgScAASBIwAJIEDIAkAQMgScAASDJkhnlsjYjPXBBbK3MaFxgA\nSQIGQJKAAZAkYAAkCRgASQIGQJKAAZAkYAAkGTLDPKIjYg9x5ntcYAAkCRgASQIGQJKAAZAkYAAk\nCRgASQIGQJIdGHCxwcdeZmrDcoEBkCRgACQJGABJAgZAkoABkCRgACQJGABJAgZAkiEzLM1K90ve\nh2G5wABIEjAAkgQMgCQBAyBJwABIEjAAkgQMgCQBAyDJkBliXjs9ttKlywUGQJKAAZAkYAAkCRgA\nSQIGQJKAAZAkYAAkCRgASYbMYZ6luyYfLvzlAgMgScAASBIwAJIEDIAkAQMgScAASBIwAJIEDIAk\nQ+Ywg1amZ63PDhcYAEkCBkCSgAGQJGAAJAkYAEkCBkCSgAGQdNuZWXAmoxZYnP+Nn+UCAyBJwABI\nEjAAkgQMgCQBAyBJwABIEjAAkgQMgCQPtDyVtXKChyhCggsMgCQBAyBJwABIEjAAkgQMgCQBAyBJ\nwABIEjAAkgyZ4X/WynWm6ItwgQGQJGAAJAkYAEkCBkCSgAGQJGAAJAkYAEkCBkCSITMwm5218tbG\n2cC5yAUGQJKAAZAkYAAkCRgASQIGQJKAAZAkYAAkCRgASYbMy/GwWlbmSz4TFxgASQIGQJKAAZAk\nYAAkCRgASQIGQJKAAZAkYAAkGTIvx5Azwd4cvuQCAyBJwABIEjAAkgQMgCQBAyBJwABIEjAAkuzA\nRmH3w7986PAlFxgASQIGQJKAAZAkYAAkCRgASQIGQJKAAZAkYAAkGTKPwnD1CObhMDEXGABJAgZA\nkoABkCRgACQJGABJAgZAkoABkCRgACQZMjOzM9fKW6Npi2k4iAsMgCQBAyBJwABIEjAAkgQMgCQB\nAyBJwABIEjAAkgyZ4TUMlnnQzoPCeYoLDIAkAQMgScAASBIwAJIEDIAkAQMgScAASBIwAJIMmQM8\n6hdmsvWba+D8LBcYAEkCBkCSgAGQJGAAJAkYAEkCBkCSgAGQZAcWYO91hJ3NjTccElxgACQJGABJ\nAgZAkoABkCRgACQJGABJAgZAkoABkGTIzKJGWCsbU8NPuMAASBIwAJIEDIAkAQMgScAASBIwAJIE\nDIAkAQMgyZCZsSy17Z3vJ4IzucAASBIwAJIEDIAkAQMgScAASBIwAJIEDIAkAQMg6bazGwWAYbnA\nAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJAED\nIEnAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyA\nJAEDIEnAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJAEDIEnAAEgSMACSBAyAJAEDIEnAAEgSMACS\n/gBnzYsPQSvn6QAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 22, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "imageplot(gamma);"]}, {"source": ["Display the optimal assignement."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAdeUlEQVR4nO3d23Kj\nuhYFUHNq//8v+zyk42AuQhICtGCM6upKx1zkpKOZJYQY3u/3CwCi+d/VDQCAGgIMgJAEGAAhCTAA\nQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIM\ngJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASP9d3QC4rWEYfj54v9/XtgRuafCjBc19omvC\njxs0ZAgRGltLr/RLQClDiNDSOKLG1ZbgguZUYHCI98o/FWHQigCD9lzpghMIMDiJVIO2BBicxNAh\ntCXAoD1ZBScQYHCIYeWfbgWDVgQYtDTOp2H0B2hOgEFjiRpL+QUNWUoKjmItRDiUCgyOIrfgUAIM\ngJAEGAAhCTAAQhJgAIQkwCCeYRisag+eBwZhTELr55/mOvJYKjCIYa3kUo3xWCowCMCDnmFOBQaR\neNAzfAgwCMPFLhgTYACEJMAACEmAQRguc8GYAINIPOgZPgQYBOBBzzAnwCAGD3qGCTcyQxg/QfV1\nU7Po4sEEGETySS/RBYYQAQhJgAEQkgADICQBBvG4AAYvAQZAUAIMgJAEGIThoV8wJsAACEmAARCS\nAAMgJAEGwZhDDz8EGAAhCTAAQhJgEIM59DAhwAAISYABEJIAAyAkAQaRmEMPHwIMgJAEGAAhCTAI\nwBx6mBNgcDjxA0f47+oGENVip2yKAXAaFRg11koKpQZwGgFGsXFKvX//LL5KWwpcGBNg1HuvfCzD\ngBMIMCrNawHVAXAmAUYZ1dX5fM1hkQADICQBRhnzCIBOCDAqzUe1Pp8RcsAJ3MhMsff7/XNVxpWZ\nM/m1ACZUYDSmnwXOIcCo8X6/F4NKegGnMYRIvXFcmep9EF9YWKMCoyW9LXAaAQZASAKMNj7DiYow\n4BwCDAIwOwbmBBgAIQkwmjGKCJxJgEG//CoACQKMllyqaUh6QZoA4xA63z2GYZh8AX09YU6AQUfm\n0TV+6eTGQOcsJUVjn7XqTzM+3W3GMMdv4/P2hmG4zRuE/VRgHOWEGFscaotbqXxaPskokQWLBBhR\nJYIqboatkWEwZwiR9v6eeHnKkFes0bZ5uBa1s9v3BedTgRFS0NG2xdKwqF6UXvAhwLibbjv4PWOe\ndxsShRYMIXKIk0cR+/c1VXL8+dEGf1+0lW3Sh/07vi84z6ACg1NljnkOoz9/G38nkzvGeDgBxlGu\nqgN67rwnX5FhVmy93++1r1sivd6/fxZfhbsSYBzuiM70b+X7ybman+kY49JqPIq4eB9bItVek/wb\nn0KGcXeugXGZv5mE+2q1xX66z+tAP/VWUbCk38j8tdLjQ1wqMA601vlO6oy65TPyR9u68nVNayWB\ngBwCjDNM4mpzm3zp4bV+TBr5ia7l8nEUY0YCYY0A41SLUw/2d9ZRYuxjYYZh7XHgsQQYx1od5Uv+\ns+25OlHavJ2R1vlXA/YTYJxkfKEr0bM2GTHrf9htPvG9+AijeZiLN43B7Qkw6ELb7FF+8QQCjMOd\n2Zl23nGv3b62qDrSOv8iQCsCjAskuuZWnW//o4iLI37v0av/PrN+K8IRrYJABBjXGL4/fk5nnE7o\nzEtZ46uJ1pHisQQYZ1jstTcXq214rqCWv24xn4UGzQkwznba8hmdFyKf9/se2XvM3w86f+/QhLUQ\nOcnnYVev37570skGfXLY4tq7rQ4FJKjAONs4xg5aQeOcIFxbwrFVDkWMcziTAON6f5PLW5cgp9U0\nZ06jSMzFtwwHjyLAOM+detXJPMDPx/MNjm3GysfwBAKMCySuG8Uaf6t4GMrnDe55p+N3d9BkTuif\nAOPm7jozIuKz0KAtAcapEpXWcVfCHiLcM2VgJwHGbT2qN3/Um4UfAoyONCzCcp4BvV9iHmDPhpmr\nWwQ1BBhnO3SocNwjnzaX/XOaq2ZSFJ1r7d41MUY4Aoy+7Im3C7vg5vMAN49Q92Ynexl2JDQBxgWO\nqE7GXfM5C7QnJk1snu6SS1bzL9Hr+wulCCMWAcaVFnvMnWOM75WPD7IWYz0Pyim8uAcBxh2sPWHk\nFae8uLZ5Io2IBBjXSI+hxbonLD16CRxEgHGxEBGVaTJ6eVyMtV/1uO3h4BQCjE4VFWGZcyKOnVg/\nP933qRevip0wX3FOXHEPAozLHDIXMeMz5+gwJCZLAH8+8BAWghJgXG+tMCrqTyd3Fs/XaD/T/Lyd\nJMN8GXuIS4ARwP6hvxNmhSyux/HVhuxDHVoJLR7cQsBE9N/VDeDR3u/3Z9mnVleDjrjOlGlYSanx\nEy/352esR6bBcVRgdK2icnqPlO5bZ/HaUtqVq17daNonDyfAeIqr7i1brslUP7CbAONim7myJ3iu\nzYn5rWDNk7P6DUpQbkCA8QiXjJt1OFRn/JA7EWBcb7MaaDX6d/4o4uKc/j3VjwSCDwFGRw566Mlx\nB5+fJYRYrYU1AowYYi3v+yq8s+qcNxXlSweZBBhdOK0mODkIx+/ra7b9+tmVR5BJgNGXnJ69NHvC\nVW9H+HtkmoDkLgQYnGR/fE6OIIp4OAFGL3LqpP1F2J6DAF0RYDzLM0PL+CG3ZDFfohpHUWm//FlE\nGIhLBUZHNgf3Sj+/eXAxBnEJMML4KrlmKw3ePoqq36DxQ+5KgNGXnE72vfJx0cH76c1vH71wENfA\nuEDOMycnj7j8KyNmR6t4SuT44E2epdlKfkv6aTNcRYBxtsWCo0mKdBVFa3a2sHT6yT3GDy98yjY9\nM4TIqRKd7zldbZ+FFwlrv/Gc3xJ6I8A4z3gY8J0xEWPP3cp1DTtZxWq/T+u4c37j4bEEGNfL6cX/\nKqfZS6Xd2HMKr+jjh6W/8fA0AoyzFcXVolZPiTSKGJTvGT9M4qBr44w5dPmMtmHW1ZPAfhoTdx5E\njFZyBRUYkSSeElmUBEHX8x1n+ebG48tm5kFwSwKMs2X2momMeX/b2Z7+C5GGSeMaEnciwDjP4kSM\nndexxnZ2xzt37zxKx9EVK8N6bx/XEWBcYzIFY4+6fj/iKGJ+U6ePvpxt0HvV+Xq9jv+Nh+gEGKdK\nD/qtPXF4s+PemUb6wT59PTdg9huP7xoCjAu8Z65qxuQzPVdjFV+l0F18+nsR+q3Rimn0dOEzRX4y\nnT1/6vzaEXJclVsHnXd+2CHImOHHnqeV8hwqMLqzf853RTCc00uu3ZKVcND9W/1WmtKLbAKMXiw+\nPKX6CKW73G+xwZ+3Np4HMXx/MN6sH9KLfAKMjix2WHW92G1yqFT6jU9e6y0hpBdFBBh9ScwnzMmk\nPUVYzz3mnmW01t5X5+/36iYQgACjO63u0Crd/dpRxLZd9vLDrDuY+ZkQfe18zifA6NpPp1YUaXfq\n/qozOBzpRQUBRo9a9WL5U/CbnK5nPb9H6UUdAUan9gwk7p/Q2G0pk9PF73/7Z+q/hXRLgNGvcYYd\nvXRhz7/7Z7YtYhKYdsgeAowY9vTOEXv2JjofmpNe7CTA6Fp1v7a/Q7w89r6Wss1uTJQkkF7sJ8Do\nXZMld28/fTFx21yH70t60YQAI4CdT/z6cXlFxQ/pRSsCDFadlnmJfryoi59PdektIaQXDQkwYqgr\np0r3OrRL3X/wtSNEKS6lF20JMMLov8tr2MImmdTtV6zbhhGLACOS/U8lLi3C+i9u5qOFfba521FN\n4hJgRFXdTR/Xv5/2hOVwpBdHEGAEU7FO0s5O88L8qGh5h1HRYZO4BwFGPEffpByon+0/G25QPtIt\nAUZst1xvfv/oaCfv17RDDiXACO+I9eY7nMoRLgCkF0cTYIRUETD360Mn9VZX5Zf04gQCjOe6pK66\n37K8c9KLcwgwompShCX2OnoU8aCnfF0eGNKL0wgw7qOHK1U/Pi3JeSRKXVKujR9eS3pxJgHGreQv\ntJGzbsWeLngYhsQjTtLnvQHpxQkEGIHN0+hVXsHk6zByupq+0UkzeA4Bxh0Mw7BnmcQrkmkjcSdv\np8Np/RPSi/MJMO4jM5DGm23ucmR33ObIPSRHD23ggQQYsa0lUDfr6qY79IJZkQs795EWHZaDPIQA\n41YyL4adVYS9vz+uX5l37ZOXlz6mHXIhAcbd9NqN3rBMkV5cS4AR3ryEKpoin1OEjTUaMRt+/6Ta\nljDe5pLwkF5cToBxc6V5kzMnMP9gr9cwGjlskHyd3EkmveiBAOMO5n1ozsWwxSLsALlV3RUzIWtI\nLzohwLiVPX1r/kDiEaVP6fzDq6ZvSC/6IcC4s/wZhvlrULXzLp2aePmEdelFVwQYN7HWnxZN0Dhl\nbY55bi0vupFow+Wz56UXPRBg3E1Ovz+Wfyt00YqL71/pzSpcfsuX9KITAoz7K+pwu+2d549oObOp\n0osOCTDuI9G3psflJq+evshv7/c4Sy/6JMC4oaKhwiIVz22ZGZL/TJ0x5/PNXT5zBNYIMJ4ocX3r\nlCJseK2sxLEWS1eliGmH9EyAcSuZQ4VrGxSd4ogduyp3pBedE2A8S/51sqLp7JkvZWbYj8Q2J8SJ\n9KJ/AozHqbh5+eDCaOOO5vPLMulFCAKMuxmnzmY1M3kps7Ou7tNXWvL5ZNbCHEcnivQiCgHG3aw9\nBHKtX04/NHKzCMsvj5Jb9nLpS3oRiADjVmYhsX3Fq84kAgtH+d6zD17fddjfWepqxDrSi1gEGPfx\n3dd/huP+xuXWwmCx4y56PObmq98NS/zza279VZMSpRchCDBuqbj/zZz1N1/PaZ+fg+Q+8fKc8kt6\nEYUA4yG217Oove71fpU/GGV8sPWXLHUIKQKMpyudeThLr9f6P1PHG328mGHvyWbHRYv0IigBBtXX\nvSYhNN0yY9J8zpHPmzQPsQgwbmneKTe4Z3mlQBlHV36G/fxZPOAFI4cv5RcBDf7Xcicr43tZQ2Tz\n3vyK6mQtI9sHjPQiOhUYt/LdES+v+J63b9lpa3fMcUiISi9uQIBxN4vd8fv9zummJ8tQZZ4wv2mL\nH3+37YyaT3pxD/9d3QBob0+nPF/8YmvGYHqWR8WrlcvYZ5Je3IYKDNJKJ8qnN1ittM7JEunFnQgw\nmNrq2dOv1o8B7hn8bHIuiEWAQcJxvfzy+OEkrtrGjBuWuRkBBgk7p1QkLpjVPOV5V1OkF7cjwKBU\ndcAcOzsjdWLpxR0JMFiwNK99elfZ/jBIH6FV2FgsirsyjR42FT2LOT11/nejlXBqHjamHXJjAgyW\nLd0Qtrrt7J9fsTEMxWu2NQkb6cW9CTBYNcuw+dOTFz9eOE7DVmWSXtyea2CQKZEBbcb9Gk61kF48\ngQoMcmRmwM9mZeGx87pXenfpxY0JMGho4XlgiQjZP2VDevFkhhAhR/uZ6InsyQy22fW5r7iSXtye\nAINMQ/KfhcdayJ6vBNrMsNEG4x3/Pnb7F7cnwCBl5QmZOdlQdNmsdK8mp4bYBBhsqB2LUwDBsQQY\nbMvLsOlVqJ2MAUKaWYiQZZxh39efxoslvufLcDQ5YyHJxyMIMKj2/v17+yEprwMrqsFFL56peIk2\neLJZCM0Xrf99YbliW7NwnM2fzaXGWICDB3ENDHIlc2g6Dz655cKxK57VMtvmq3nDMLiKxr2pwGAh\nmSY/F1vR9bfhfPfvC2av+cYLRyz5qdxMKT/j3JUKjKdbDIDxJ7cSoujmsLGFXHm/3/sevFJzTzQE\nZRIHj7b0tJRh6aWpn9hY3GYlgY56+PJ6hZf1aE2IS4DxXLP4WZwcv+CTN+MYKwmhglwpP/iYDOPO\nBBjMba9DOE6UwnTZGLFcO+PWuVzo4nFcA4NNTS4sDaO/662fV6XF4wgweH1H1NrquhXTNF6vr5pp\nbffN9Xyrs1OqcWcCDObm891rxv2yTSaSzCdlDEut+t1oNSClFzcnwHio2vipKYYqJsevnz1teU6/\nmfTckgCD13ePP+n958+NrMmGzSdYzhQEXk46WpiD+zELkcepmgdx3By/tZnuqSV653PrPx8vTrv/\nvOXJ/EkITQXGs+yuQtbuGq47yK4dF4uqxeHKwpWFIQa/jhFb/loYWx335lzB+fab92bln32zDdu3\nJBetXu8HnxsQYASWSIWS1XhXj7C+19fFsPQPUfmp9y6fkWjPpDF+/AlNgBHV0jKGr8UnkqRzLlHD\nrT/9a/lEW408RDJovzb7fKwU4x4EGFFVPKZkU0bmTcuj7PKrTV1VEcbzg7xkGLcgwIjqtwteez5k\nyryvzynXFo+T3mDUyMr0qrqel3U0GUZ0AoyQ1suvf6/nHyoxQy89ea9k0sSu8qs2w1YL07W01hsQ\niwAjpCYBlhld1VYCbB5mC/E2qRE3M2xWU863Xw1gpRhBuQ+MkPb3s+Pu+6BCZOdxJvmaLgEz7lCe\nvjr8cpcYQanAiKp6EsfRhddXU7ZHEZeLs/yLVSuRs1qYFs31WDuLfoMeCDACKy0XzoyuxbNsWb05\nujzDtkdWN5MsZ95jYnc4mgAjqur0qpiOUWH//cvpe7Ez7tTOKkw3UypNB8KFBBghVU92X9ux7Q9C\ncnpFgXYZlnWnc6rY+v7n36CnDoTrCDDiyc+DoukJrX4WMif15T9LLP/4RV+Z/CVIEsOROhAuJMAI\npqLf35rm0HIS+Txd0ssNf14tGuGsy7C127fXWvjXtqWjKcK4nAAjkgPS699W8x0rZKbX5rlyltjP\nvOc65wayxEn/bbPYyOQp4AQCjBjaRdcrPcW8uGVLp1tKr+KCb08pVpdhq+ddbF7y+HACAUYAe1Zb\nfy3ESapDrviJyBiXqyz4Tsuw+Tbpy2DSix5YiYOuLa5AMXf+DV6Lp3u/30uny/nMssW4Sg8ezi+q\nvVZ+A5i0dn7Yr1fHbchsPRxMgNGvzOhaWxRqJU42jlbdwv2lW36TMounyRdn8/gbq1X9/oFOGEKk\nRwdVXZurTxX9OGSvjtFgxDLz9rVzhhNzWgInEGB0p+iGraJlNWYbbyx+kVASCfNX9+blWKsMy9ym\nqDFwKAFGRw693NXwRubSgqbVrP2KDBu/Wpph6RbmrKMIhxJg9GJP4VU0nW/tft5M5em1rPk1s5xS\nrG2GlTYpsa+OiAoCjOsdN2Y436XVnV6Lh5pvUHcjc34b0seszrD8zeZbVtyFnXkimBBgXOmEKfJH\npFdOdOW/uqclcwdlWP6WS7fBfR3E2COtCDAuU5RedRlwQnpVVIT77R9uPTTDXuvfuK9txrvnnQXG\nBBgXKLpKVJ0Q56fXmT9N45u9SkftijKsyZY5y1NZ3YNSAoyz7UmvhhMFdx7nwuiatyExapeeSV89\nb7764tlfmxf3zWgJfAgwzpMZXa/jJwruOc4lY4Zr9mTYa/Z1PiHDfl5a21OAUUSA0VKic9+8dN9q\nZlqT9Mq8orP/x2dnHGbOnmiVYYtn3NzstfQFVIGxnwCjjfzqamJz0lp1M5qk15nRVXfwzDLxtR5j\nezIsvX0iX10DYz+L+dLee32MaLplsmjoJL2GYTj+itfX16zot4HFWNq8w3qyV3rCxfyMmdvPN/uL\n2Enb0qeEJSowGqgov5pfUjouvfa3LXmi1VKkug5LyxxObHiPXc6dzjnHgQkBRr26YcON2dXdpNeh\nVVfzAHu1qNtehw0npi+85Z8OxgQYlZqnV9tLVnsOMnHEz8jv6VLT8XaNoI4PtKJVhs13z99sMq4I\nRQQYNXamV8PiZn96XTKilQywvUOIk93S36qcaaKHZhhUM4mDYtUTDodfR7SkeXpVzCKpOH/zI85b\nnH4PORe9cr5lk3ouZ1ppw/8JPJMA42IXjhwmOtnrioNj+/T07MTEG8/MsJx8yow62CTAuFKf6VXX\npCKjswy/oTWM0+u4Zqxl2GYplpk0OfmUPxEfElwDo9hnNaCiXme+fW/pdf7PQsMQLb1BePPUpcOG\na+1Jb+mSGHsIMIrV3PX1s+P4Mz2l14U/BQ2jdC3DEitcFC3YsWbnlA0ZRjUBRrH9ASa9FhXdRJzY\nfVHpwOCZpdjxC51wT66BUay0fzmiN6obZJt3lF31lUWZkdi9yV5rafoe/ZlvvHl8l8RoSAVGjdzr\n+eNdPp/cXX7tuUS0sw1H21mETQ5SdKhEKVZ6gW3z+IYTaUKAUS9rcKnR3I3qnr3nMcNFTTJs56kX\npW66Lr/ZeW0XGUY+AUZ7zZdlapVeIf63Xxhgr/T3bm2Xn1er6jyXxNhDgHGU5ksdVi+ttLMB57s2\nw16J4cT5lp+Xqhe+MpxILZM4OMq40zkzvRZvntUDFtl/lTH/+ImbnesOznMIMI6yv4yoS6/5J8Ol\nVw8NXh7cS/+zNsPW9pVhpBlC5BANb9jaM4ku7n/vy0cR5y3JtOc6pUtiFBFgtHdyel27quFB+gmw\nV94t0nu+6WYnUscQIo1JryZ23tTcVs4Xc89wX86NzIYTmRNg7DKMTF46Yj2k+akXX4qeXh1KrFqy\n+AtHzXpjMoxChhCpVLfyXv4xi5bXu+v4UlcDiT+K1k48Ys2U9He/+rxEJMCocVV6LXZVd02vV5cB\n9sq4Ub3tmss5MzsSjeHGBBg1FhfH2/O0lIr0und0ffSZYa+tUuzQDFseYPy8OmsMdyXAKNZkadfF\nA67tuDis9IT0enUcYK+tUqxths2P+fXSZMd95yUKkzhoqaK3SHdz85kaj0qvV2fTESfS1yn3z7nI\nmZ24vGPFyQhIgHGIzL5mM73G//x0Z89Jr/4dnWHzU7xf06eRvWZrguw/KSH8d3UDuKfSW7jSFznW\n+kHp1YP5bxUfk09OIq3mXHt25nZUYLSU/+vuWg6tjRmm97q3nkcRP9LfjvFbOPNdPOc/yTMJMIr9\ndUbfidUkvSYnWuv4dEwdOmE48et0k1OsfMyNmYVIpc2bgXJ2zKktHjJdflPP0xEnMvMp/42sTXyd\nnOZtDv3DqMCotNY7lKbX5piS9Aon89u0vw7LnMrBXanA2OunG6qbtbF5W49fqMcCFWGv1nXY5r3z\nFcckOrMQ2avi1+38O1Lffq0OK/07ykfp1MTFw0msZzKEyBnm06k/HyeWOf+3wVGNiifEdMSJ6V1c\nKyutjN/R2lL0daPW3JgKjMbyl/956XqeYfx/IPGfYTHDFvMv1lAqx3ENjGaKyoK1O5cX/zsOS7s8\nVv5Fx958Vd7jzyf3ivhOOYchRNooWKdua8xweuSq9txeoFHEufS3fzq3MPI75VCGEGls406drega\nXPS6u817uT5/Cy7SVGA0kBgA/HwmXXh9LdPw/We+wcNFnMqx03PeKUUEGCfZ7IMy443n8F0nTYDR\nkcWgkl5zoYuweC2mV66BcZLMHBJXT7B2pdPSUBRRgdHAeH36CX3QQSIWYYtXOkt3hA8BRmOeakHC\nxlrPS/N3YI0bmWlgzwQN9oh7U/Pru/GWa6GCAGMv6XWhOy2qNP+PdIM3xaFM4mAXD+uiFf9/KOUa\nGPWk1+UiTuWAVgQYlaQXcC3XwNiW/u3ef6HL3elKGORzDYwUEzSAbhlCJNd7tlCC9OqEK2E8kwBj\n1XiN+b9F5Ucf6y6BCwkwuANFGA8kwNiQs+gqwPkEGNyES5I8jQBjw+KAlFGqnhlF5CEEGAAhuZGZ\nFLcwh+OmZp5DBUZKohPUPwLXEmBsWAwq6dUt8+l5DktJsU1cAR1SgcHdKMJ4CAEGQEgCDG7IqC9P\nIMDgzowicmMCDICQBBjck6kc3J4AAyAkAQa3pQjj3gQYACEJMLgzRRg3JsAACEmAwc0pwrgrAQZA\nSAIMgJA8kRkeYTJ+6AefG/A8MHiiT55JMuIyhAj3l5i+YWYHcQkweJD36M+HDCMoAQY39zda+P15\nQ4dEJ8DguWQYoQkwAEISYACEJMDguUzeIDQBBjf3txbi9+elF9G5kRkeZDG03MtMUCowuL9EREkv\n4rIWIjyIFRG5EwEGQEiGEAEISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIA\nAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQk\nwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAh\nCTAAQhJgAIT0f04XIpMBAk1YAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 23, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "[I,J,~] = find(gamma);\n", "for k=1:length(I)\n", "    h = plot( [X0(1,I(k)) X1(1,J(k))], [X0(2,I(k)) X1(2,J(k))], 'k' );\n", "    set(h, 'LineWidth', 2);\n", "end\n", "myplot(X0(1,:), X0(2,:), 10, 'b');\n", "myplot(X1(1,:), X1(2,:), 10, 'r');\n", "axis equal; axis off;"]}], "metadata": {}}], "nbformat": 3, "metadata": {"kernelspec": {"name": "matlab_kernel", "language": "matlab", "display_name": "Matlab"}, "language_info": {"mimetype": "text/x-matlab", "name": "matlab", "file_extension": ".m", "help_links": [{"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md", "text": "MetaKernel Magics"}]}}}