matlab/optimaltransp_2_benamou_brenier.ipynb

{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Optimal Transport with Benamou-Brenier Algorithm\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 details how to compute optimal transport maps by solving a space-time\n", "convex variational problem.\n", "\n", "\n", "This tour implements the algorithm described in:\n", "\n", "\n", "Benamou J.-D.; Brenier Y.,\n", "_A computational fluid mechanics solution of the Monge-Katonrovich mass transfer problem_,\n", "Numer. Math. 84 (2000), pp. 375-393\n", "\n", "\n", "The algorithm detailed in this paper is introduced as an\n", "augmented-Lagrangia method over the Legendre-Fenchel dual variational\n", "problem. In this tour, we introduce a proximal algorithm (Douglas-Rachford) directly over\n", "the primal algorithm. Both algorithm can be shown to be equivalent, but we feel that explaining a primal\n", "algorithm was conceptually simpler and fits nicely in the well developped framework of proximal splitting schemes.\n", "\n", "\n", "Special thanks to Jalal Fadili for advices for the implementation of the\n", "proximal solver."], "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_2_benamou_brenier')"]}, {"source": ["Usaful helpers."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["mynorm = @(a)norm(a(:));\n", "sum3 = @(a)sum(a(:));"]}, {"source": ["Optimal Transport of Densities\n", "------------------------------\n", "We consider transport $T : [0,1]^2 \\mapsto [0,1]^2 $\n", "with periodic boundary conditions.\n", "\n", "\n", "A valid transport should push forward the measure $ f_0(x) d x$\n", "onto $ f_1(x) d x $. In term of densities, this corresponds to\n", "the constraint\n", "$$ f_0(x) = f_1(T(x)) \\abs{\\det( \\partial T(x) )} $$\n", "where $\\partial T(x) \\in \\RR^{2 \\times 2}$ is the differential of $T$\n", "at $x$. This is known as the gradient equation.\n", "We call $\\Tt(f,g)$ the set of transport that satisfies this constraint.\n", "\n", "\n", "Optimal transport looks for a mapping that solves\n", "$$ W(f_0,f_1) = \\umin{T \\in \\Tt(f_0,f_1) } \\int C(x,T(x)) d x $$\n", "where $C(x,y) \\geq 0$ is the cost of assigning $x \\in [0,1]^2$ to $y \\in\n", "[0,1]^2$.\n", "\n", "\n", "We consider here the $L^2$ optimal transport, so that\n", "$C(x,y)=\\norm{x-y}^2$.\n", "\n", "\n", "The optimal transport geodesic $f(x,t)$ is defined as\n", "$$ f(x,t) = f_1( (1-t) \\text{Id} + t T(x))\n", "   \\abs{\\det( (1-t) \\text{Id}  + t \\partial T(x) )}. $$\n", "\n", "\n", "Benamou and Brenier showed that this geodesic solves the following convex\n", "problem over $f(x,t) \\in \\RR^+, m(x,t) \\in \\RR^2 $\n", "$$ \\umin{ (m,f) \\in \\Cc } J(f) =\n", "      \\int_{[0,1]^2} \\int_0^1 \\frac{\\norm{m(x,t)}^2}{f(x,t)} d t d x, $$\n", "where the set of constraints reads\n", "$$ \\Cc = \\enscond{(m,f)}{ \\text{div}(m)+\\partial_t f = 0,\n", "      \\quad f(\\cdot,0)=f_0, \\quad f(\\cdot,1)=f_1  }. $$\n", "\n", "\n", "Note that this convex program is very challenging because:\n", "1/ The functional $J$ tends to zero when $f(x,t)$ tends to\n", "infinity at some points, so that it is not coercive, which makes the\n", "proof of existence of minimizers non-trivial.\n", "2/ The functional $J(m,f)$ tends to infinity when $f(x,t)$ tends to\n", "zero at some points which makes the use of gradient descent methods impossible\n", "(its gradient is not Lipschitz).\n", "\n", "\n", "We propose to use a proximal method scheme (Douglas-Rachford algorithm)\n", "that can handle this kind of problem. It is equivalent to the algorithm\n", "initially proposed by Benamou and Brenier.\n", "\n", "\n", "Discrete Variational Problem\n", "----------------------------\n", "Densities is discretized on a spacial rectangular grid of size $N = n \\times\n", "n$. The time domain $[0,1]$ is discretized with $p$ points."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["n = 20;\n", "p = 20;"]}, {"source": ["Shortcut to generate Gaussian function."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["[Y,X] = meshgrid(linspace(0,1,n), linspace(0,1,n));\n", "gaussian = @(a,b,sigma)exp( -((X-a).^2+(Y-b).^2)/(2*sigma^2) );\n", "normalize = @(u)u/sum(u(:));"]}, {"source": ["Load two discretized densities $f_0,f_1$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sigma = .1;\n", "rho = .05; % minimum density value\n", "f0 = normalize( rho + gaussian(.2,.3,sigma) );\n", "f1 = normalize( rho + gaussian(.6,.7,sigma*.7) + .6*gaussian(.7,.4,sigma*.7) );"]}, {"source": ["Display $f_0,f_1$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAALMUlEQVR4nO3dz0tU\nfxvH4Rkdc9QMQ6MSF+1rEUTQqqC/o380aN0iighqk1GQWVnkaP4cdWaexbN74Nw92A9963Vtb26/\np+lML8/ifL7t0WjUAoA0Yyd9AQBwHAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAk\nAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyA\nSAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkY\nAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQS\nMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACI\nJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoAB\nEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQB\nAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAInVO+gL+una7\nfbxpp1N9OJOTk02jmZmZYnFqaqqYjo01/krR7/eLxZ2dnWK6t7fXNDo8PCwWh8NhMT0nRqPRSV/C\ncdR3PvzS6b/zPYEBEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARDp7J/EURsfH28adbvd\nYnF+fr5ptLi4WCxeu3atmBbHf/R6vWJxdXW1mK6trTWNtra2isX6nI7T/6I+cIZ5AgMgkoABEEnA\nAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDINLZP4mj3W4X0+IkjpmZmWJxaWmpaXTnzp1i8datW8V0\namqqafThw4di8fnz58X04OCgadTv94vFwWBw7CnAX+UJDIBIAgZAJAEDIJKAARBJwACIJGAARBIw\nACIJGACRzvuLzJ1O4ydQv8i8uLjYNLp9+3ax+ODBg2I6OzvbNHr16lWx2Ov1iunq6mrTaH19vVjc\n398vpgAnyBMYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRzv5JHLXinI7x8fFisdvt\nNo0uXbpULC4sLBTT4iSOy5cvF4vT09PFdGJiomk0Nlb9ElOfY1JMR6NRsQjw+zyBARBJwACIJGAA\nRBIwACIJGACRBAyASAIGQCQBAyDSeX+ReTgcNo0ODg6KxV6v1zT6+PFjsfjmzZtiOjMz0zRaXl4u\nFr99+1ZMd3d3m0ZHR0fFYv0+sreVgRPkCQyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQM\ngEhn/ySO+rSI4hyK7e3tYnFlZaVp9OzZs2Lxx48fxXRycrJp9Pnz52Lx7du3xXR9fb1pVB84UpxU\nAnCyPIEBEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARCpXR9UcQa02+1iOjbWmPDiUIxW\nqzU7O9s0WlhYKBbn5uaKaafTeDbK1tZWsVgf8NHr9ZpGe3t7xeJgMCimZ/7m+a/QP2Z95/OXHPtj\nP4W32Sm8pP/hCQyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkc77i8zFtHjHudVqTUxMNI3q\nN6AvXLhwvOs5PDwsFvv9fjEtdo+OjorFM397/D9CPwQvMv+O4utfnDbwy2lxL9Vf8BM5UuD03/me\nwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiHTeT+I4tuJF/foIj2OfDDIcDovFelr8\nLZ/5G+D3hX5ETuKo1d/T6enpptGVK1eKxfn5+WJ6cHDQNFpbWysWe71eMS1O8fidu/f03/mewACI\nJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiNQ56QtIVZx8UR+K8ZfORzj978xzztV3fnEu\nxvj4+PF+bP1NrH/s1atXm0b3798vFu/du1dMNzY2mkZPnjwpFl++fFlMi3M6BoNBsZjOExgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAInmR+V/zxjHnU6dT/Wtz8eLFptHc3FyxWLwBvbW1VSzW\nb/hev369afTw4cNi8dGjR8V0dXW1afT9+/di8d27d8V0c3OzaeRFZgA4dQQMgEgCBkAkAQMgkoAB\nEEnAAIgkYABEEjAAIgkYAJGcxAH8GcWhGK1Wa3p6upjeuHGjaXTz5s1icWJiomm0vLxcLK6srBTT\n4sScvb29YvHnz5/FdHt7u2nU7/eLxeFwWEzPLU9gAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnA\nAIgkYABEchIH8GfUJ3F0u91iurS01DS6e/fu8X7s0dFRsbi2tlZMv3792jR6/Phxsfjly5diurW1\n1TR68eJFsbi5uVlMB4NBMT3DPIEBEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJIXmYE/YzQa\nFdPhcFhM9/f3m0YbGxvFYvEi887OTrHY7/eLafEfffr0abH4+vXrYlq8W72+vl4s1n+W+rM9wzyB\nARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQyUkcwJ9RnwdRnyXx/v37plFxekWr1ep0\nGv8R+/TpU7FYH/Cxu7t7jFGr1frx40cxLQwGg2NPzy1PYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKA\nARBJwACIJGAARGqPRqOTvoa/q91un/QlkC30O3La7vzx8fFiOjk52TTqdrvF4thY42/h+/v7xWI9\nrY//OCdO/53vCQyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkbzIDL8Q+h3JuvOLqy1eVa4N\nh8NiGvrX+i+d/o/IExgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJE6J30BANWhD4PB\n4F9eCUE8gQEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEE\nDIBIAgZAJAEDIJKAARBJwACI1C7+T94AcGp5AgMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZA\nJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJ\nGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABE\nEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAA\niCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKA\nARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAk\nAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyA\nSAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkY\nAJH+Awh9ysUqSEb+AAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "imageplot({f0 f1});"]}, {"source": ["Boundary conditions, either periodic or Neumann."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["bound = 'per';\n", "bound = 'neum';"]}, {"source": ["We use first order finite differences with periodic boundary condition to\n", "approximate the spacial derivatives."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["if strcmp(bound, 'per')\n", "    dx = @(u)u([2:end 1],:,:)-u;\n", "    dy = @(u)u(:,[2:end 1],:)-u;\n", "else\n", "    dx = @(u)u([2:end end],:,:)-u;\n", "    dy = @(u)u(:,[2:end end],:)-u;\n", "end"]}, {"source": ["The adjoint operators are backward derivatives."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["if strcmp(bound, 'per')\n", "    dxS = @(u)-u+u([end 1:end-1],:,:);\n", "    dyS = @(u)-u+u(:,[end 1:end-1],:);\n", "else\n", "    dxS = @(u)[-u(1,:,:); u(1:end-2,:,:)-u(2:end-1,:,:); u(end-1,:,:)];\n", "    dyS = @(u)[-u(:,1,:), u(:,1:end-2,:)-u(:,2:end-1,:), u(:,end-1,:)];    \n", "end"]}, {"source": ["Check that |dxS| and |dyS| really implement $d/dx^*$ and $d/dy^*$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["Should be 0: 7.58e-15\n", "Should be 0: 7.76e-15\n"], "output_type": "display_data"}], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["fprintf('Should be 0: %.2e\\n', certify_adjoint(dx,dxS,[n n p]));\n", "fprintf('Should be 0: %.2e\\n', certify_adjoint(dy,dyS,[n n p]));"]}, {"source": ["Define spacial gradient and divergence, satisfying $\\text{div}=-\\nabla^*$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["grad = @(f)cat(4, dx(f), dy(f));\n", "div  = @(u)-dxS(u(:,:,:,1)) - dyS(u(:,:,:,2));"]}, {"source": ["Check that |div| really implements div."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["Should be 0: 5.04e-15\n"], "output_type": "display_data"}], "prompt_number": 13, "cell_type": "code", "language": "python", "metadata": {}, "input": ["fprintf('Should be 0: %.2e\\n', certify_adjoint(grad,@(v)-div(v),[n n p]));"]}, {"source": ["We use first order finite differences for the time derivatives.\n", "Note that zero is padded at the end to keep the same dimensionality."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 14, "cell_type": "code", "language": "python", "metadata": {}, "input": ["dt  = @(f)cat(3, f(:,:,2:end)-f(:,:,1:end-1), zeros(size(f,1),size(f,2)) );\n", "dtS = @(u)cat(3, -u(:,:,1), u(:,:,1:end-2)-u(:,:,2:end-1), u(:,:,end-1));"]}, {"source": ["In the following, we group the dicrete variables as $w=(m,f) \\in \\RR^{N\n", "\\times P \\times 2} \\times \\RR^{N \\times P}$, which is stored in a matrix\n", "|w| of size |(n,n,p,3)|, so that $m$ is stored in |w(:,:,:,1:2)|\n", "and $f$ is stored in |w(:,:,:,3)|.\n", "\n", "\n", "We thus introduce the following convex functional for $w=(m,f)$\n", "$$ J(w) =  \\sum_{x,t} \\frac{\\norm{m(x,t)}^2}{f(x,t)}. $$\n", "where $(x,t)$ runs over the discrete 3-D grid.\n", "\n", "\n", "The affine set of constraints reads\n", "$$ \\Cc = \\enscond{w}{A w = r_0}$\n", "      \\qwhereq A w = (\\text{div}(m) + \\partial_t f, f(\\cdot,0), f(\\cdot,1) ) \\in \\RR^{N \\times (P+2) }  $$\n", "and where the right hand size reads $r_0 = (0,f_0,f_1)$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 15, "cell_type": "code", "language": "python", "metadata": {}, "input": ["A = @(w)cat( 3, div(w(:,:,:,1:2))+dt(w(:,:,:,3)), w(:,:,1,3), w(:,:,end,3) );"]}, {"source": ["Its adoint reads, for all $r=(s, r_0,r_1)$ where $s \\in \\RR^{N \\times\n", "P}$ and $r_0,r_1 \\in \\RR^N$:\n", "$$ A^* r = ( -\\nabla s, \\partial_t^* s + U(r_0,r_1) ) $$\n", "where $U(r_0,r_1)$ is the zero-padding operator that puts zeros for\n", "intermediate times between $t=0$ and $t=1$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 16, "cell_type": "code", "language": "python", "metadata": {}, "input": ["U = @(r0,r1)cat(3, r0, zeros(n,n,p-2), r1);\n", "AS = @(s)cat(4, -grad(s(:,:,1:p)), dtS(s(:,:,1:p)) + U(s(:,:,end-1),s(:,:,end)) );"]}, {"source": ["Check that |AS| really implements $A^*$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["Should be 0: 3.30e-15\n"], "output_type": "display_data"}], "prompt_number": 17, "cell_type": "code", "language": "python", "metadata": {}, "input": ["fprintf('Should be 0: %.2e\\n', certify_adjoint(A,AS,[n n p 3]));"]}, {"source": ["Define the right hand side $r_0$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 18, "cell_type": "code", "language": "python", "metadata": {}, "input": ["r0 = cat(3, zeros(n,n,p), f0, f1);"]}, {"source": ["Proximal Operator of the $J$ Functional\n", "-----------------------------------------\n", "Given a convex functional $J(w)$, its proximal operator is defined as\n", "$$ \\text{Prox}_{\\la J}(w) = \\uargmin{q} \\frac{1}{2} \\norm{w-q}^2 + \\la J(w). $$\n", "Being abble to compute this proximate mapping is crucial to be able to\n", "use various proximal-splitting convex optimization methods.\n", "\n", "\n", "The discrete functional to be minimized in the Brenier-Benamou problem reads\n", "$$ J(w) = \\sum_{x,t} j(m(x,t),f(x,t)) $$\n", "where we introduced the following convex function\n", "$$ \\forall (m,f) \\in \\RR^2 \\times \\RR^+, \\quad j(m,f) ) = \\frac{\\norm{m}^2}{f}. $$\n", "\n", "\n", "Define the $J$ functional."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 19, "cell_type": "code", "language": "python", "metadata": {}, "input": ["J = @(w)sum3(  sum(w(:,:,:,1:2).^2,4) ./ w(:,:,:,3)   );"]}, {"source": ["The proximal operator of $J$ can thus be computed by applying the\n", "proximal operator of $j$ to each component, i.e.\n", "$$ \\text{Prox}_{\\la J}(w)(x,t) = \\text{Prox}_{\\la j}(w(x,t)). $$\n", "\n", "\n", "The proximal operator of $j$\n", "$$ (m,f) = \\text{Prox}_{\\la j}(m_0,f_0) $$\n", "is obtained by solving the following couple of equations (corresponding\n", "to the annulation of the derivative of the optimization problem defining\n", "the proximal map)\n", "$$\n", "      2 \\la \\frac{m}{f} + m-m_0=0 \\qandq -\\la \\frac{\\norm{m}^2}{f^2} +\n", "      f-f_0 = 0.\n", "$$\n", "Its solution is obtained by solving a cubic-equation\n", "$$\n", "      m = \\frac{m_0}{1+2\\la/f} \\qandq\n", "      P(f) = f^3 + (4\\lambda - f_0)f^2 + (4\\lambda^2-4\\lambda f_0)f -\n", "      (\\lambda \\norm{m_0}^2 + 4\\lambda^2 f_0) = 0.\n", "$$\n", "(only the largest real solution of this equation should be considered).\n", "\n", "\n", "Build the coefficients of the polynomial $P(f)$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 20, "cell_type": "code", "language": "python", "metadata": {}, "input": ["PolyCoef = @(m0,f0,lambda)[ones(length(f0),1), 4*lambda-f0, 4*lambda^2-4*f0, -lambda*sum(m0.^2,2) - 4*lambda^2*f0];"]}, {"source": ["Helper to compute the leading (largest) real root of a cubic polynomial."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 21, "cell_type": "code", "language": "python", "metadata": {}, "input": ["extract = @(A)A(:,1);\n", "CubicReal = @(P)real( extract(poly_root(P')') );"]}, {"source": ["Define the proximal operator of $j$. Note that it can operate in\n", "parallel over arrays |m| of size $k \\times 2$ and |f| of size $k \\times 1$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 22, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Proxj0 = @(m0,f, lambda)cat(2, m0 ./ repmat( 1+2*lambda./f, [1 2]), f );\n", "Proxj  = @(m0,f0,lambda)Proxj0( m0, CubicReal(PolyCoef(m0,f0,lambda)), lambda );"]}, {"source": ["Define the proximal operator of $J$ using the proximal operator of\n", "$j$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 23, "cell_type": "code", "language": "python", "metadata": {}, "input": ["ProxJ = @(w,lambda)reshape( Proxj( ...\n", "                   reshape(w(:,:,:,1:2), [n*n*p 2]), ...\n", "                   reshape(w(:,:,:,3  ), [n*n*p 1]), lambda ), [n n p 3] );"]}, {"source": ["Orthogonal Projection on the Constraints\n", "----------------------------------------\n", "The proximal operator of the indicator function $G=\\iota_\\Cc$ of $\\Cc$ is the\n", "projector, and does not depend on $\\la$.\n", "$$ \\text{Prox}_{\\gamma \\iota_\\Cc}(x)_i =\n", "      \\text{Proj}_\\Cc(w) = w + A^* (A A^*)^{-1} (r-Aw). $$\n", "\n", "\n", "Tolerance and number of iterations for the conjugate gradient."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 24, "cell_type": "code", "language": "python", "metadata": {}, "input": ["opts.epsilon = 1e-9; \n", "opts.niter_max = 150;"]}, {"source": ["Adapt conjugate gradient fucntion to handle variables that are not\n", "vectors."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 25, "cell_type": "code", "language": "python", "metadata": {}, "input": ["flat = @(x)x(:);\n", "resh = @(x)reshape(x, [n n p+2]);\n", "mycg = @(B,y)resh( perform_cg(@(r)flat(B(resh(r))),y(:),opts) );"]}, {"source": ["The operator $(A A^*)^{-1}$ can be computed using conjugate gradient."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 26, "cell_type": "code", "language": "python", "metadata": {}, "input": ["pA = @(r)mycg(@(s)A(AS(s)),r);"]}, {"source": ["Define the projection operator $\\text{Prox}_{\\la G}$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 27, "cell_type": "code", "language": "python", "metadata": {}, "input": ["ProxG = @(w,lambda)w + AS( pA(r0-A(w)) );"]}, {"source": ["Check that $\\text{Prox}_{\\la G}$ implements the projection on the\n", "constraint $Aw=y$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["Error before projection: 1.72e+03\n", "Error before projection: 1.47e-06\n"], "output_type": "display_data"}], "prompt_number": 28, "cell_type": "code", "language": "python", "metadata": {}, "input": ["w = randn(n,n,p,3);\n", "err = @(w)mynorm(A(w)-r0)/mynorm(r0);\n", "fprintf('Error before projection: %.2e\\n', err(w));\n", "fprintf('Error before projection: %.2e\\n', err(ProxG(w)));"]}, {"source": ["Douglas-Rachford Solver\n", "-----------------------\n", "The discrete optimal transport problem can be written as\n", "$$ \\umin{w} J(w) + G(w) \\qwhereq G=\\iota_{\\Cc}.  $$\n", "\n", "\n", "The Douglas-Rachford (DR) algorithm is an iterative scheme to minimize\n", "functionals of the form $J+G$\n", "where $J$ and $G$ are convex functions for which one is able to\n", "comptue the proximal operators.\n", "\n", "\n", "A DR iteration reads\n", "$$ \\tilde w_{k+1} = \\pa{1-\\frac{\\mu}{2}} \\tilde w_k +\n", "  \\frac{\\mu}{2} \\text{rPox}_{\\gamma J}( \\text{rProx}_{\\gamma G}(\\tilde w_k)  )\n", "  \\qandq w_{k+1} = \\text{Prox}_{\\gamma G}(\\tilde w_{k+1}). $$\n", "\n", "\n", "We have use the following shortcuts:\n", "$$   \\text{rProx}_{\\gamma G}(w) = 2\\text{Prox}_{\\gamma F}(w)-w $$\n", "\n", "\n", "One can show that for any value of $\\gamma>0$, any $ 0 < \\mu < 2 $,\n", "and any $\\tilde w_0$, $w_k \\rightarrow w^\\star$\n", "which is a minimizer of the minimization of $J+G$.\n", "\n", "\n", "To learn more about this algorithm, you can read:\n", "\n", "\n", "_Proximal Splitting Methods in Signal Processing_, Patrick L. Combettes\n", "and Jean-Christophe Pesquet, in: Fixed-Point Algorithms for Inverse\n", "Problems in Science and Engineering, New York: Springer-Verlag, 2010.\n", "\n", "\n", "\n", "Set the value of $\\mu$ and $\\gamma$.\n", "You might consider using your own values to speed up the convergence."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 29, "cell_type": "code", "language": "python", "metadata": {}, "input": ["mu = 1;\n", "gamma = 1;"]}, {"source": ["Define the rProx operators."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 30, "cell_type": "code", "language": "python", "metadata": {}, "input": ["rProxJ = @(w,tau)2*ProxJ(w,tau)-w;\n", "rProxG = @(w,tau)2*ProxG(w,tau)-w;"]}, {"source": ["Number of iterations."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 31, "cell_type": "code", "language": "python", "metadata": {}, "input": ["niter = 200;"]}, {"source": ["Initialization using linear interpolation of the densities."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 32, "cell_type": "code", "language": "python", "metadata": {}, "input": ["t = repmat( reshape(linspace(0,1,p), [1 1 p]), [n n 1]);\n", "f = (1-t) .* repmat(f0, [1 1 p]) + t .* repmat(f1, [1 1 p]);\n", "m = zeros(n,n,p,2);\n", "w0 = cat(4, m,f);"]}, {"source": ["Display the initialization."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAARzElEQVR4nO3dW29V\n5doG4FHa0pbNKlhQJEQlHhk4ICEmHmnikT/CP2risYkGY4x6gOAmVsBSSukGWvbrYOVb68v7TO1o\n6Szznr2uwyeds4M+s9wMc/uOiZcvX3YAkObI674AANgLAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAk\nAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQB\nBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEG\nQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZApKnXfQFDNzExseNkaqr9OczMzNS3On78\neDOZm5trJkeOtP8mePz4cTN5+PBhM9na2momT58+rRfw4sWLOhw1L1++PPhvul9b7rPibphbjlhx\ndwi2XFfc2fLocQcGQCQBBkAkAQZAJAEGQCQBBkCk8W8hVpOTk81kdna2mSwsLNQXnj9/vpmcO3eu\nmdQS1OrqajO5detWM1laWmomGxsb9QJqnWn0a0Kvy9623GfF3TC3bMW7Mrwt1xV3tjx63IEBEEmA\nARBJgAEQSYABEEmAARBp/FuI9bS02lyqZ6NduHChvtXVq1ebyeXLl5tJPVTtt99+aybXrl1rJk+e\nPGkm9dS1ruueP3++4+Rw2q8t91lxN8wtW/E/OMgtDzwD05ZHjTswACIJMAAiCTAAIgkwACIJMAAi\nHcYWYj3lrDaX6mlpXddduXKlmXzyySfN5OTJk83k+++/byZ9TlS7f/9+vYDt7e06pNu/LfdZcTfM\nLVvxPzjILdcVd7Y8etyBARBJgAEQSYABEEmAARBJgAEQafxbiFWfE9Xqc127rvvXv/7VTM6cOdNM\nannp9OnTzeTYsWPNZHp6upkcOTLg3xb1yuvEo13/Y29b7rPibphb7rPizpb/z/C2PLCFeJBbtuI+\n3IEBEEmAARBJgAEQSYABEEmAARDpMLYQX7x40UzqQ1TrEWdd1/3xxx/N5Keffmom9Si2GzduNJO7\nd+82k0ePHjWTZ8+e1QuoxSRVpb+zty33WXE3zC1b8a4Mb8t1xZ0tjx53YABEEmAARBJgAEQSYABE\nEmAARBJgAEQa/xp97afWVuvm5mYzWVxcrG/1zTffNJOVlZVmMjMz00xu377dTH7++edmUh86XtvA\n3aDSMP+xX1vus+JumFu24n9wkFuuK+5sefS4AwMgkgADIJIAAyCSAAMgkgADINLE2B8iWZ/VXZ/w\nXRtHAx8oXp87furUqWYyNdUWOzc2NppJrTzV80a3trbqBTx//ryZjOD6Xssl7deW+6y4G+aWI1bc\nHYIt1xV3tjx63IEBEEmAARBJgAEQSYABEEmAARDpMLYQ+3SZpqen61vVgtPRo0d3fPOnT582k8eP\nH+/4NfWQty6hFNSNTD9tb1vus+KBb76PW+6jXkA11OfWj/2WB/6E97bl2jmsJicnm0n9gwxUP0L1\noMU9L2v0/8JxBwZAJAEGQCQBBkAkAQZAJAEGQKTD2ELsY2AFqA771KJqKahPTSh3LyPST+ujLrTP\nigcO65+6z977XNLAGmRt1tUL6FOQy+qnHeSWB36vvW25vlWtQdbDGOfn5+sF1M5hfQB0PbCx7r2n\n0f9byB0YAJEEGACRBBgAkQQYAJEEGACRtBBfw5uP9888qJ+253fer73X95mdnW0mp0+frm9VS2u1\nn3bv3r1msr6+3kyy+mn7teW9nan4d8NGbSHWV509e7aZXLlypZlcvny5vnltGF67dq2ZXL9+fcdX\n9SnEdgl/U7kDAyCSAAMgkgADIJIAAyCSAAMg0tTrvoBso9/SYRjq83OnptpfpVpsq0XB+vmpp+Qt\nLCzUCzh//nwz2d7ebib1LMRHjx7t4ZLGT/0z1oUeO3asvvDEiRPNpG65/pBr5a/u9OrVq83ks88+\nqxewvLzcTFZXV5vJ4uJiM3n48OGOlxTKHRgAkQQYAJEEGACRBBgAkQQYAJG0EOGf1H5a13Vzc3PN\nZG/9tK2trWZSC3LPnz+vF1Dbg/XLxukx3/urHk5Yy5/nzp2rL7x48WIzqe3TWgK8e/fujpdUK6O1\nOtgN+gjt41O2E7kDAyCSAAMgkgADIJIAAyCSAAMgkgADIJIaPfxPz0fLz87ONpP5+fkdX1j7zbU8\n/eTJk2Zy7969egH1hbVYv7a2tuPXHKrK9X/VLdca/Ztvvllf+MEHH+z4wvr/M9TzdldWVprJ119/\nvePXdIO69devX9/xa8bm6N7KHRgAkQQYAJEEGACRBBgAkQQYAJG0EGHXatOsnqlaW4j1VbUeVouC\nDx48qBewubnZTGqfsF7SwHOBD6H6s6qLqD3PbtCPvbZGt7e3d/ya+j4//PBDM/n111/rBdQNrq+v\nN5N6QrQWIgCMFgEGQCQBBkAkAQZAJAEGQCQtRPif2k8b2NyrRa/6wnrgXu2n1c5h/XYDL6AW26rD\nec5hH/UnU1dz+/bt+sK6i8nJyWayvLzcTGrnsH67Oqn1wq5ff3KMO4eVOzAAIgkwACIJMAAiCTAA\nIgkwACJNjH1VqZbBGKrX8ok64C3Xcw5rG63qUxhL+X0csy3XhU5PT9cvq89frpdUC6J1knIo5eh/\nGt2BARBJgAEQSYABEEmAARBJgAEQSQuRfTZm/TQGGvstD/xefS6g/mRy/44d/St3BwZAJAEGQCQB\nBkAkAQZAJAEGQCRPZAZoDSzgjX4r77BxBwZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQB\nBkAkAQZAJAEGQCQBBkCkCcdTApDIHRgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACR\nBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEE\nGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQY\nAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgA\nkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACR\nBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEE\nGACRpl73BQzdxMTEjpMjR9ogn56erm81OzvbTI4ePbrjmz979qyZbG9vN5MnT57s+Kqu616+fFmH\no+a1XGT9sff5mrr3OpmcnOxzAS9evNhx8vz582YSsdCBRnbL7KPR/3y6AwMgkgADIJIAAyCSAAMg\nkgADINL4txCr2mWqZcITJ07UFy4sLDSTU6dONZNaWtvc3GwmKysrzWRtba2ZbG1t1QsYpxrbsO2t\na1o/CXUy8M1ra/Tx48fNpHZNLRRehTswACIJMAAiCTAAIgkwACIJMAAijX8LsRbGalGwHnJ45syZ\n+lbvv/9+M3nnnXeaSS2tLS0tNZNffvmlmdQ22tOnT+sF1OP1lNb6q3ufmZlpJrV9evz48T5vVVuj\ntX1al1Un9ZMA/B13YABEEmAARBJgAEQSYABEEmAARBr/FmJVz8SrbbQ33nijvvDixYvN5NKlS82k\nFhp///33ZvLw4cNmcv/+/WaysbFRL2BgNZGB+pyFWCujJ0+ebCbz8/P1zWsLsb5Vn9MR60LrZSua\nwt9xBwZAJAEGQCQBBkAkAQZAJAEGQKTD2ELc27N6u66bm5trJrW0Vr+mHq9XS4+11VYvkl2pP8A+\ne6+LGPhJmJpqf3Fqw7DPTm0ZXoU7MAAiCTAAIgkwACIJMAAiCTAAIh3GFmJ9rnE9kq4eV9h13fLy\ncjNZXFxsJrVheOfOnWaytrbWTGqHrV5k51i83ejz+ON6XOH29nYzGfhJ6PNE5idPnjST+rRlC4VX\n4Q4MgEgCDIBIAgyASAIMgEgCDIBIWohdN6h7du/evfrCmzdvNpPNzc1mUk/JW11dbSa3bt3a8X1q\nQa5TWgP4fw5jgHEY1LCvLfb6fy9sbGw0k4H/kqiH8NbSfC3W17eq/5bybxToz39CBCCSAAMgkgAD\nIJIAAyCSAAMg0sTYt576PEi+dt/n5ubqW508ebKZnDhxYsc3rx399fX1ZlJPjK2ttq5faa3PU+qH\nuvTX8onq86euq6ln8k5PTzeT+tkY+O1qxbGeEL2PLURbHns9/7Djt+VdcQcGQCQBBkAkAQZAJAEG\nQCQBBkCkw9hCrPr007p+FbXh9dMGvnm9gDqpK64XsI9Puw/qp9W99+msDlT/1H0ahnXSpyI7cGjL\nQfpseeDeD8OWd8UdGACRBBgAkQQYAJEEGACRBBgAkbQQd/GqPVfUGns+0vDYsWPN5OzZs81kYWGh\nmdRjFZeWlprJ6upqM6ntpoHXubev2Xf71U/b8/vUF/b5OdRX9VlxN8wt91xf9JYPWP1bYm+/yN3h\n2PKuuAMDIJIAAyCSAAMgkgADIJIAAyDSgOO26P6mftPnLLt6iGKfSe0lDjyM8a233momH3/8cTP5\n6KOPmsmDBw+ayZdfftlMvvvuu2ZSu0zdoGPWxkzPOmjPEwsbfbbcZ8XdMLc89iseqGfBuO6rvnC/\ntlxX3Nly4Q4MgEgCDIBIAgyASAIMgEgCDIBIWoi70OeByCdPnmwm8/PzzaQWnDY2NprJwJrQ22+/\n3Uw+/fTTZvL55583k1u3bjWT5eXlZnLz5s1msra2Vi9gbMpL/Q2sF544caKZnDp1qpnsbct9VtwN\nc8uHcMXdoC3XFXcHu+W64s6WC3dgAEQSYABEEmAARBJgAEQSYABE0kLchT7Pz3333XebyaVLl5rJ\n9PR0M7lx40YzWVxcrBdQT+rb2tpqJuvr681kc3OzmTx+/LiZ1APcDqc+D8/tuu69995rJvu15T4r\n7mz51fTZcl1xd7BbrivubLlwBwZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkNfpdqO3b2dnZZnLhwoVm\n8uGHH+74qmfPnjWTpaWlegF//fVXM/niiy+ayZ07d5pJPV3022+/bSZjfNznrvRZcTfMLfdZcWfL\nr2Zvv8jdwW65rriz5cIdGACRBBgAkQQYAJEEGACRBBgAkbQQd6EewVnPzdze3m4mDx48aCa1ufTw\n4cNmUs/oHPhWX331VTP58ccfm0mtRd2/f3/HCzhUR4L+V58Vd8Pccp8Vd7b8avb2i9wd7Jbrijtb\nLtyBARBJgAEQSYABEEmAARBJgAEQSQtxF2qZpxZ+fv3112ZSi0NTU+2P/c8//2wmtafUdd2jR492\nnKysrNQXNurZaGN8Wtqu9FlxN8wt91lxZ8uvZm+/yN3BbrnPirtDv2V3YABEEmAARBJgAEQSYABE\nEmAARJqox4KNmYmJieG9+eTkZDOZmZlpJvW0tPpA2Hrw2sCj2GoJagS9lk/U8LZcV9wNc8sRK+4O\nwZbrijtbHj3uwACIJMAAiCTAAIgkwACIJMAAiKSFOPRvV3tKVT2cLXcvY9ZP6/ntbPkAHOSWB34v\nWx417sAAiCTAAIgkwACIJMAAiCTAAIikhcg+G/t+Gp0tHw6jnw7uwACIJMAAiCTAAIgkwACIJMAA\niCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIg0/of5AjCW3IEBEEmAARBJgAEQSYABEEmAARBJ\ngAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmA\nARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYAB\nEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQ\nSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEOnfkgnVXY8AojYAAAAASUVO\nRK5CYII=\n", "output_type": "display_data"}], "prompt_number": 33, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sel = round(linspace(1,p,6));\n", "clf;\n", "imageplot( mat2cell(w0(:,:,sel,3), n, n, ones(6,1)) , '', 2,3);"]}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["\n", "ans =\n", "\n", "    0.1994\n", "\n", "\n", "ans =\n", "\n", "   2.8807e-07\n", "\n"], "output_type": "display_data"}], "prompt_number": 34, "cell_type": "code", "language": "python", "metadata": {}, "input": ["w = ProxG(w0,gamma);\n", "mynorm(A(w0)-r0)/mynorm(r0)\n", "mynorm(A(w)-r0)/mynorm(r0)"]}, {"source": ["__Exercise 1__\n", "\n", "Implement the DR iterative algorithm on |niter| iterations.\n", "Keep track of the evolution of the energy $J$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["[\bWarning: Input tol may not be achievable by CGS\n", "         Try to use a bigger tolerance]\b \n", "[\b> In cgs at 297\n", "  In perform_cg at 56\n", "  In @(B,y)resh(perform_cg(@(r)flat(B(resh(r))),y(:),opts))\n", "  In @(r)mycg(@(s)A(AS(s)),r)\n", "  In @(w,lambda)w+AS(pA(r0-A(w)))\n", "  In exo1 at 6\n", "  In pymat_eval at 38\n", "  In matlabserver at 27]\b \n", "[\bWarning: Input tol may not be achievable by CGS\n", "         Try to use a bigger tolerance]\b \n", "[\b> In cgs at 297\n", "  In perform_cg at 56\n", "  In @(B,y)resh(perform_cg(@(r)flat(B(resh(r))),y(:),opts))\n", "  In @(r)mycg(@(s)A(AS(s)),r)\n", "  In @(w,lambda)w+AS(pA(r0-A(w)))\n", "  In exo1 at 6\n", "  In pymat_eval at 38\n", "  In matlabserver at 27]\b \n"], "output_type": "display_data"}, {"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAQc0lEQVR4nO3d0Xai\nSBRAUWpW/v+Xax7oEAREUaDqVu39lI5Z04yxPd4LMSnnPABANP+VPgAA+ISAARCSgAEQkoABEJKA\nARCSgAEQkoABEJKAARCSgAEQkoABEJKAARCSgAEQkoABEJKAARCSgAEQkoABEJKAwU1SSpsfP/sa\nYJ+AQUVyzhoGbxIwuFtKaR6qxQcaBm8SMChPseADAgZljJPWOI0Nv2NZ6YOCSAQMrmW6gov8lD4A\naN90cqv0gUBTbC3gbvvbQrtEeJMVItzNdYZwCq/1AAjJBAZASAIGQEgCBkBIAgZASO0HzNVeAE1q\n/AeZx3qlNAyDjgEcUP816o0HbKbM96KSH0p1GA6j2mNwGNUeRulDeK35FeL8VwgWPAwATtZ8wIac\nh+nVTEoyBtCI9gM2qmAiB+BMVSxbr7PYJk/jV9P/0wDfquRU3L5eJrAFi0SA6PoK2Pz1hIYBhBbv\nMvrp4s7N8fblbw7MWboAWhByAss5b/5GpXFp+/KXLc0vSgQgqHgBO/e8ooYBBBVvhTjs/oj4eoW4\n+OLxJotEgIUQ774xFzJgY4QW9/X8os/5x88mtqlhKbmqHmDvpX+dgq0QL7pPI3ynAHgQbAKbX6Cx\nmLc2b3r1X5MugKiCBWzYitPLbeHuf80iESCkYCtEABgJmB8LAwhJwB5oGEAUAjYM3pweICAB+8ci\nESAWAdugYQD1E7A/FokAgQjYA4tEgCgEbEnDAEIQsD0aBlAtAdvgZBhA/QRsm0UiQOUE7CkNA6iZ\ngL1FwwBqI2B75ifDNAygKgL2goYB1EnAXtMwgAoJ2Fs0DKA2AvYuDQOoioAdoGEA9RCwYxYNkzGA\nUgTssJyNYgDlCdiHNAygLAH7nHUiQEEC9pXF+9ZrGMBtBOxb61NiMgZwAwE7h1EM4GYCdhqjGMCd\nBOxkRjGAe/yUPoDD0m8T8mMr0mMrFrfeafybp8MZPyh3OABtihew4TdOKaV5peYfpwoGHxkDuFS8\nFeLL0WoRtrLWG8UK2grQgpAT2KEBq/hqcTGKjR9XU1iAf2rYXR0SMmDTCnF903r8qmQas1EEKlfb\niZiXgq0QQ9ynO1xqD3CWYAHLOadf04uFqWpVnf3a4VJ7gO/FWyGuEzV9JkS9RjaKAF8KNoE1xq8W\nA/iYgJXnrBjABwSsCkYxgKMErCIaBvA+AauLhgG8ScCq45QYwDsErEZOiQG8JGD10jCAHQJWNQ0D\neEbAaqdhAJsELAANA1gTsBg0DGBBwMLQMIA5AYtEwwAmAhaMhgGMBCweDQMYBCwoDQMQsKj8+mag\ncwIW2NQwQxjQIQGLTcOAbglYeBoG9EnAmqJhQD8ErAUuSgQ6JGCN0DCgNwLWDg0DuiJgTdEwoB8C\n1hoNAzohYA3SMKAHAtYmbzQFNO+n9AEcln5nirx6kt65qUM5/xu/UtIzoEEhJ7Ccc845PW7Hxj9u\n3tQtb9IBNCxewPanq5RSSskENtEwoFXxVojDbFW4MHVr3rDFF/fcNrtEYEe43VXIgI0RevO+7rlY\no+lk2KBhwHPzZ8sQMQu2Qgxxn1bIhfVAe4JNYPMLNBYLw82bmJjDgMYEC9iwFafpM7q1T8OAlgRb\nIfIlu0SgGQLWHQ0D2iBgPdIwoAEC1ikNA6ITsH5pGBCagHVNw4C4BKx3GgYEJWBoGBCSgDEMGgYE\nJGD8o2FALALGHw0DAhEwHmgYEIWAsaRhQAgCxgYNA+onYGzTMKByAsZTGgbUTMDYo2FAtQSMFzQM\nqJOA8ZqGARUSMN6iYUBtBIx3LRomY0BZAsYBORvFgFoIGIdpGFADAeMTGgYUJ2B8SMOAsgSMz7ms\nAyhIwPiKyzqAUgSME2gYcD8B4xwaBtzsp/QBHJZ+nx3z/Cnz8abNW7lazn/pSmnwHQAuFS9gw2+c\nUkrrSulWWYuGDYOMAVeJF7D9RI1D2Pxr0uM+S+GuNt7BRjEIJ0Xb/scL2LB7L6+HM8UqwjoRwtl5\n6V+nkBdx5Jw3s6RVVXFZB3CpYAHbeVEQ4vVCb/ykM3CdYCvEnPP6KsRxYbh5E8U5JQZcJFjAhq04\nOd1VP6fEgNMFWyESl3UicC4B4z6LwUvDgG8IGLfy5r/AWQSMAqwTge8JGGVYJwJfEjCKsU4EviFg\nFKZhwGcEjPKcEgM+IGBUwToROErAqIiGAe8TMOqiYcCbBIzqOCUGvEPAqJFTYsBLAka9NAzYIWBU\nzToReEbAqJ03nQI2CRgBOCUGrAkYYWgYMCdgRKJhwETACMZlHcBIwIjHKTFgEDDi0jDonIARmIZB\nzwSM2DQMuiVghOeyDuiTgNECl3VAhwSMdmgYdEXAaIqGQT/iBSz92vmCO4+H2mgYdCJewIZhyDnn\nnDdDpV4MGgZ9iBewvPjtGjMppZ1b6YpLE6F5P6UP4BOHxqzFFytcP8Zv9fT9T2n5q8WAuXAbrMAj\ny2Le2gyVmYzhcYvo4QDvCPHkGWyFuPMCIf8ajFk8ckoMmhQsYOO1G6OpUuHGXu6nYdCeeOfA1tPV\n4jPGLzbl7HwYNCXYBAbfcGkitETA6It3TYRmCBg90jBogIDRKQ2D6ASMfmkYhCZgdE3DIC4Bo3cu\nTYSgBAxcmgghCRj8o2EQi4DBHw2DQAQMHmgYRCFgsOSyDghBwGCDyzqgfgIGT2kY1EzAYI+GQbUE\nDF7QMKiTgMFrLuuACgkYvMVlHVAbAYMDNAzqIWBwjIZBJQQMDnNKDGogYPAJp8SgOAGDz2kYFCRg\n8BUNg1IEDL7llBgUIWBwAqfE4H4CBqfRMLiTgMGZNAxuI2BwMqfE4B4/pQ/gsPT7fJDnzxOvboI7\njQ/AKV0pDR6ScLqQE1jOOeectl7Z7twEN7NOhEvFC9jOdGXwojYaBteJt0IcZqvCd25afEbkuFnO\nD7vEYbBOpFLhdlchAzZG6NkKcXGTYlGcU2KEMH+2DBGzYCvEQ7MXVMU6Ec4VbAKbX6AxvVhIKS2u\n3TB1USfrRDhRsIANW3GaPqNb1M86Ec4SbIUIbbBOhO8JGJShYfAlAYNivOkUfEPAoKTFCTANg/cJ\nGBTmd4nBZwQMqqBhcJSAQS2cEoNDBAwqYp0I7xMwqI6GwTsEDGpknQgvCRhUyhX2sE/AoF5OicEO\nAYPaaRhsEjAIwCkxWBMwiME6ERYEDCLRMJgIGARjnQgjAYN4XGEPg4BBUE6JgYBBYBpGzwQMYnNK\njG4JGIRnnUifBAwaoWH0RsCgHdaJdEXAoCnWifRDwKBBGkYPBAzaZJ1I8wQMmuUNO2ibgEHLnBKj\nYT+lD+Cw9PtPMC9eXu7eBD3L+S9d4wf+idCAkBNYzjnnnLZeTO7cBD2zTqQ98QK2M10ZvGCHdSKN\nibdCHGarwme3zku2+GKRo3PWiTwTbncVMmBjhNb39fiZRaIUCxbmDRuGISUNYxgeny1DxCzYCvHl\nfSpX8A7rRBoQbAKbX6AxtWrcGY6fdyEivM86kdCCBWzYKtP4GcWCD1gnElewFSJwOutEghIwYBi8\ndyIBCRjwj1GMWAQMeGAUIwoBA5a87xQhCBiwwTqR+gkY8JR1IjUTMGCPUYxqCRjwmlGMCgkY8BZX\ndlAbAQPeZZ1IVQQMOMY6kUoIGHCYUYwaCBjwIaMYZQkY8DlXdlCQgAFfsU6kFAEDTmCdyP0EDDiH\nUYybCRhwJqMYtxEw4GSu7OAeAgaczzqRGwgYcBXrRC4lYMCFjGJcR8CAyxnFuIKAAXdwZQenEzDg\nJut1oozxDQEDbmUU4ywCBtzNlR2cQsCAMqwT+VK8gKVfO19w5/EAHzOK8Y2f0gfwiZzzMAwppfy4\nTZcuiCjnv3SNHyzOk8GmeBNYfv7Qzjnv3ApUy5UdfCDkBHZo0lp8scJBncZ/mvNRzD/Wm4VbYoUM\n2LRCfP+LgRCsEwuaP1uGiFmwFWKI+xT4his7eFOwCSznPDVserGwvpoDiM4oxkvBAjZsrQQXnxEz\naMO8YYOzYqwEWyECXfH2iewQMKB2LrJnk4ABAbiygzUBA8KwTmROwIBIjGJMBAyIxyjGIGBAUK7s\nQMCAqFxk3zkBA2IzinVLwIDwXNnRJwEDGmGd2BsBA9phFOuKgAGtMYp1QsCABrmyowcCBrTJRfbN\nEzCgZUaxhgkY0DijWKsEDOiCUaw9Agb0wijWGAED+mIUa4aAAd0xirVBwIBOGcWiEzCgX0ax0AQM\n6J1RLCgBAzCKhSRgAP8YxWIRMIA/RrFABAxgySgWwk/pAzhT+n2U5cWjD+Cg8VlkStf4gaeWqrQ2\ngeWcc87J6yXgDEaxmjUVMIMXcDpnxarV1ApxmG0Rn31G5IAP5PzQrZQaXCeG2121FrCxT/Nvg2IB\np2j+rNj82TJEzNpZIYa4u4HobBTr0c4ENr92w9QFXKqHjWL92gnYoFvAjZrfKNavnRUiwP1cZ1+Q\ngAF8xVmxUgQM4ARGsfsJGMA5jGI3EzCAMxnFbiNgACczit1DwAAuIWNXEzCAC9koXkfAAK5lFLuI\ngAHcwSh2OgEDuIlR7FwCBnAro9hZBAzgbkaxUwgYQBky9iUBAyjJRvFjAgZQmFHsMwIGUAWj2FEC\nBlALo9ghAgZQFxl7k4AB1MhG8SUBA6iUUWyfgAFUzSj2jIAB1M4otknAAGKQsQUBA4jERnEiYADB\nGMVGAgYQ0noU6y1jAgYQ1WIUGzrbKAoYQGzdbhR/Sh/AYen3O5MXLzx2bwJoW84P3Ro/bvu5MF7A\nht8+pZTmoZr/cXETQA/Gp71Fxhp+Loy3QlQmgB39bBRDTmDD7oyVHr9Xiz/qH9CDDzaKKVro4gVs\nvIvXHco5TzfNvw2KBfTp6EZxcVLmugM7S7yADU+a5LwXwNoiYy1d3BEsYOOLgsXVhmO65oOXkgHM\nNXmNYrCAbZZp+qRuATzT3jWK8a5CBOBjLV2jKGDXquREqMOYcxhVHcPgMB7dcxhtZEzAADoV/e2A\nBQygX6HfDljAAHq33igOQ4CrOxr/2alKltoAQfwVof44NB4wAI6Kcnm9FSIAD0LUaxAwAIISMABC\nEjAAQgr2XoghrN9TuOC7DE9v0l/qGNbvvFzkMKa/uuBh7Hwv7jye2g5j/cfbDmPnLy17V9x8GC8f\nCdW+T7qAXWL+NvmLfyQ3P2Ou/947j2H9T7HgXVHwMOY/zrH+Xtz23Vn/VEmRB+r6MPbvnxuOodQD\ndfOuKPJA3Xkk3HkYR1khnq+Sb3A9D7WUUj0HU8T4635KH8XyMEod0uIwijw2Nr8j9z9Q63lglD6E\nD5nArtL5U/ZkvZQoZfMXduOBOvJAjfhIELDzLYbugsew+P2f3do81UGFD9RwT6DnKvVAreGR8BkB\nu0Txh8L6LEvZ46FOtT1Qyx5Mz4Le+QJ2ssXQM+64i1/DU+oY1n9vz4fx7GCKfHfqfKAOhe6NSh4h\nRQ7jnUdCDQ+MTV71ABCSqxABCEnAAAhJwAAIScAACEnAAAhJwAAIScAACEnAAAhJwAAIScAACEnA\nAAhJwAAIScAACEnAAAhJwAAIScAACEnAAAhJwAAIScAACEnAAAhJwAAIScAACEnAAAhJwAAIScAA\nCEnAAAhJwAAIScAACEnAAAhJwAAIScAACEnAAAhJwAAIScAACOl/EEzXyfef+oQAAAAASUVORK5C\nYII=\n", "output_type": "display_data"}], "prompt_number": 35, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo1()"]}, {"collapsed": false, "outputs": [], "prompt_number": 36, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Display the resulting density $f(x,t)$ for $t$ from 0 to 1."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAYO0lEQVR4nO3da29c\nZ7kG4JW2cewc7KQ59JC0lFZUgiJUgZDgC0j8Dv4ookJISAiqVlURB4WGtjk0ds4nxzm1zf6AxN56\n73fjt1M7mWdyXR8fzYyXZ63x3anuPGvP48ePJwCo5rmnfQAAMAsBBkBJAgyAkgQYACUJMABKEmAA\nlCTAAChJgAFQkgADoCQBBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgADoCQBBkBJAgyAkgQY\nACUJMABKEmAAlCTAAChJgAFQkgADoCQBBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgADoCQB\nBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgADoCQBBkBJAgyAkgQYACW98LQPYNft2bNn28kL\nL7Tvw759+/KlDhw40ExWVlaayXPPtf9N8ODBg2Zy9+7dZnLv3r1m8ujRozyAr7/+Oofz5vHjx0/+\nh+7du7eZPP/889s+Jk/f/v37t31M96W++uqrZrK1tdVMNjc3t31MnvfuSR95k/M6n1m+VPf63G2r\nq6vNJM/y0tJSMxk5p/l57/4FGPl037lzp5nkeR/5vA+e99nO8shfxal35PPGNzAAShJgAJQkwAAo\nSYABUJIAA6CkxW8hpmwuLS8vN5OjR4/mE1999dVm8vLLLzeTLDTeuHGjmXzxxRfN5NKlS80ku0xT\nr6r0VCp/JWRhLLtnJ06caCZ5il966aV88eyjZq9sfX29mZw/f76ZXL58uZlkm7HbRtuphuEONhWf\ngLza87OcpybPcn5sjx071kwOHTqUB5BvV366z50710zyvGd3cdBge3CGx+TnpYSSBw0AAgyAkgQY\nACUJMABKEmAAlLT4LcRs4Iw0l06dOpUv9ZOf/KSZ/PCHP2wmWXX77LPPmskHH3zQTB4+fNhMuj2l\nrKjlhH/LdYVra2vN5PXXX28m77zzTjN566238sWPHDnSTHKr4enTp7c9yJGtmF07tRNv0JyUFfNq\nz9Lv4cOHm0me5TfffLOZvPLKK80kT/HUW7R47dq1ZpJ/AfKc3r59e9vHdAvGc3Ii5odvYACUJMAA\nKEmAAVCSAAOgJAEGQEnPYgsxm0vZQsydeNM0vfvuu83kl7/8ZTPJFWoff/xxMxnZjnj9+vU8gPv3\n7+eQrjzL2UJ87bXXmskPfvCDbSdTb71ethCzsZabD/O857XRbaM94Tsyz618k7M9ePLkyWbyne98\np5nkdsTjx4/nj8s/FLdu3Wom2R/Os5yTNLgDM3cYzrYdsejV4hsYACUJMABKEmAAlCTAAChJgAFQ\n0uK3ENPIdsS8R/M0Taurq81k5EauWYvav39/M8nFfd0bpI50h57NezSPnNN9+/Y1k+wlZhstN+lN\nvdV5m5ubzWRjY6OZ5PWTJbrBEzrbWS7aNPuP2RrFeXf1PMvZVOzeiTs/3SON4nxWHnae0O6a05HO\n4WyTonwDA6AkAQZASQIMgJIEGAAlCTAASnoWW4i5ZCxviJzlommazp0710z+9re/NZMsQX3yySfN\nJHfi5Sa9L7/8Mg8gq0rPZucwZatq5Cznbsl8TFe+7fnEPKf54/IsD95qeeRKWLw2Wrea28j3Ifuo\nWQvMOnGuu5x6d1vOzYd5kHkpjtxavdtCzN9up3qJRfkGBkBJAgyAkgQYACUJMABKEmAAlCTAAChp\n8Wv0WTzN7nJuYj1//ny+1Pvvv99Mrl271kxyY+zFixebyT//+c9mcv369WbS7XN37zLO1KsF51m+\nfft2M7lw4UIzyVOT/y5i6v1DiDyD//jHP5rJ+vp6M8lifa557e6VHmld59Uy8o8uugXrOWldZyH+\n0aNHzeTWrVvNJD+k+Zj8uOUrT703+cqVK80kr407d+7kSzWyoJ9XwjT27y5mPssV+QYGQEkCDICS\nBBgAJQkwAEoSYACUpIU4TdN09+7dZpLVwalXVTp79mwzye5QdpCyFnXz5s1tf9Zkde83kS2yfJM/\n/fTTZpKrUbNeOE3Tiy++2Eyy4phLnPNm83kpZumxW4PMhmEWGu/du7ftj6slz06e5fxwnTlzppms\nrq42k2zlZVNx6n0As7aaPy5PxPHjx7eddJf55pLxnMzWVBxZlDyHSh40AAgwAEoSYACUJMAAKEmA\nAVDS4rcQ08j94LNUNvVuH56b0LLOlEWpfJ18zOANxfm3fGfyDcyuadYCt7a2mkmutpum6ciRI80k\nz+mlS5eaSRbb8mrJzmF3bV3+vrkyMXtl+Q50L7NC8oOTrbysBebnPbumb7zxRv64fEuzyJrN5L17\n9zaTt99+u5msra01k+55z6sxV25ubGw0k26leTH4BgZASQIMgJIEGAAlCTAAShJgAJSkhThNvS5W\nt++XS8aytJbdoVxbN3KrXIXDb2TkDRy5626WCbsNriyp5rUxctvfvA9v9txmbiHmVZ2/b25Q7F54\nc3I1jtyRPJuW+TvmetL8IOeVMPXe5OwcZqEx9xyePHmymRw9erSZ5O2nuy+VZzlXL+YhzXa37jnk\nGxgAJQkwAEoSYACUJMAAKEmAAVDSs9hCHNHtO+Ww2xDbVtHCzzwbeUtHmld5QrtttJGlgll+y6Zi\n3r87J91b5eaR5xOXlpa2neRv173yZ7vOd9zIGRxp5WUbM0uA+ayp9z7kw/LayM2ZMzt8+HAzef31\n15vJhQsXmkneizwvRXdkBoAnR4ABUJIAA6AkAQZASQIMgJK0EL8VfcK5NdupyQV03QJe1uGy/JaP\nyV2IWR3Mu/cOthBHDP52O/XjdtzILsSR1aP5JueNsPft25cvPrKpMtcq5k2ic09mFkS7uxAPHTrU\nTPLI8zH5UnlCR97bOeQbGAAlCTAAShJgAJQkwAAoSYABUJIWIotgtppclspmbiHmAeSL57NyFWGW\n37Kf1j3OnIwUI+ekXjhop7qX2dzLyeBuwCw05prBK1euNJNr1641k7W1tWaSaw+7Py73OmbZdeR9\nm5N1l9+Ub2AAlCTAAChJgAFQkgADoCQBBkBJAgyAktToeXZldTgn3QrybDX67HNn4zlvUd/tc2e3\nPtvb9+/f33aSh120Tv1f5Ns+sjE5/xFC1/LycjPJRv7ly5ebyeeff95MVlZWmkn3RNy8ebOZnDlz\nZtsfl7/LyFbiEnwDA6AkAQZASQIMgJIEGAAlCTAAStJChP+Vfa1s901jO3/zMVkdzB+XLcRuDTJ7\ndPmw3BSci18Xz0ibLt+ZO3fuNJOtra184v79+5tJLl9eXV1tJtevX28mGxsbzSSvhIsXL+YB5NV4\n6dKlZpJNxZGuqRYiADw5AgyAkgQYACUJMABKEmAAlKSFyGIaaVVlOyt7erk/sPvE7BzmvsSR9Xr5\nrOwlTr0eXR7SyMLGkZvNT3NcUcsD666ObGTDcH19feTHra2tNZO8YLLQmCcid2Bm5zBXGk6985Uv\nlZeZFiIAzBcBBkBJAgyAkgQYACUJMABK0kJkEYx0qGbrJWbLa+oVvbKFmC81Ug8bLApmsW3ESOdw\nnttoI925kePPdy/3B3bPe943Od/SrDjmS400VLsHMNga3ZYWIgA8TQIMgJIEGAAlCTAAShJgAJSk\nhchi2qlWVbfvN9IwHDFSKus+ZrbfrmjT7D92rzuXdzre3NzMh2XDMA9gthWUT30p5U71G58w38AA\nKEmAAVCSAAOgJAEGQEkCDICStBDhvxnsfe1eiat6dfAJ272tmNOsXdORF9/Vs5wvPnMNct74BgZA\nSQIMgJIEGAAlCTAAShJgAJQkwAAoSY2eRbBTLeRF6qw/4a72EzDbMt+detb0LZYv7xT/WqPhGxgA\nJQkwAEoSYACUJMAAKEmAAVDSnqI7HAF4xvkGBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgAD\noCQBBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgADoCQBBkBJAgyAkgQYACUJMABKEmAAlCTA\nAChJgAFQkgADoCQBBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgADoCQBBkBJAgyAkgQYACUJ\nMABKEmAAlCTAAChJgAFQkgADoCQBBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgADoCQBBkBJ\nAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgADoCQBBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQ\nkgADoCQBBkBJAgyAkgQYACUJMABKEmAAlCTAAChJgAFQkgADoCQBBkBJAgyAkgQYACUJMABKEmAA\nlCTAAChJgAFQkgADoCQBBkBJLzztA9h1a2trzeSFF9rfeu/evc1keXk5X2ppaWnbJ+bk+eefbyZf\nf/11M7l3714z2draygPI4cOHD5vJV1991UweP37cTPbs2ZMvPiKf+Nxz7X8D3b59e7YX/zb27dvX\nTEYONSf5rO57NfLiKU9EnqyRSfelcjKbwWsjL7wnYObrltns1EW1e3wDA6AkAQZASQIMgJIEGAAl\nCTAASlr8FmIWabIwtrKy0kyOHDmSL5XD1dXVZrJ///5mkr3EBw8eNJOrV682k/X19TyArH6NlMFG\nCnI72EucEyNFwTw12VDNEmn3pUbe5OwT5pUw0lQctIOnZm7PMs8438AAKEmAAVCSAAOgJAEGQEkC\nDICSFr+FmCWubJodOnSomZw4cSJfKoe5aPHgwYPNJHuJX3755bYH8OjRozyAXDN49+7dZjKyvmwH\nO4dz0k/L9ZLZHhzZgTmyynLqdQ7zYfmYvBTz3cvHdK+E+V9SB0+Ab2AAlCTAAChJgAFQkgADoCQB\nBkBJi99CzL5WNs1ypeHRo0fzpY4dO9ZMDhw40EyyT5gvnv20LMjduHEjDyAXJN68ebOZZB9vttsH\nD9YL56SFmEbe5JFdiN1fcKSFmC+VkzxZuR2xewCztRBHTtbcnlBIvoEBUJIAA6AkAQZASQIMgJIE\nGAAlLX4LMWUZLNcV5pLDaZoOHz7cTLKFmJ3DfFYeQG7A6x7A0tJSMxm5gW9W5kZWGnbbaIUqanmo\nWRQcmczcQhx5qdyK2V29CHT5BgZASQIMgJIEGAAlCTAAShJgAJS0+C3EkQ5ePmZ5eTlfaqRhmO3B\nnOQB5K2Wu220kc5hFttG3oHB9l2hFmIa+a1n3oU4slYxDb7tKZcoptnap1CIb2AAlCTAAChJgAFQ\nkgADoCQBBkBJi99CzKJX9vTu3bvXTPLGuFOv1pVLFA8ePNhMVlZWtv1xW1tbzWRzczMP4OHDh81k\npFGZ3cXssOVjBluIamzA07L4AcazYCRHR0J6ZCdvd5il+XzMyH835L+L6Dbm84mzPSb5LxIK8b8Q\nAShJgAFQkgADoCQBBkBJAgyAkha/hZj1sEePHjWTGzduNJONjY18qRdffLGZZI1+aWmpmWRrP3/c\nF1980Uxyve/U6+i/9dZbzSRLazdv3tz2xbOx1u2wzW2NPg8j34c8EVn5S/nPEqZew7D7sG0PIC/F\nfNtzTfDUu6oH64tUMfjJmq1rujB8AwOgJAEGQEkCDICSBBgAJQkwAEpa/BZilnmy+pU9vc8//3zk\nxXMt78mTJ5vJ8vJyM7l8+XIzOXv2bDPpdpB+9KMfNZNTp05t++Pyd/nwww+bSdYgc3HwPBtpIeZ5\nz9M3Uibs/riUB3D//v1mkr3ELJrmFTX16q/ZLM3LLHdGj/Qwp7npmi6MvMyyVpqTqdc5zKs6z+kC\nNxV9AwOgJAEGQEkCDICSBBgAJQkwAEpa/BZiNnCyHnb37t1m0u3gbW5uNpM7d+40k6wAHT16tJnk\nosVr1641k2PHjuUB/PznP28mP/7xj5vJ6upqMzl37lwzyYWNv//975vJpUuX8gBG7in8VGRNbqSv\nNXJtZFFwmqYDBw40kyx/5hMfPHjQTPJEvPHGG83k3XffzQN49dVXm8nFixebyZ/+9Kdmcvr06WaS\nl3T3hM7JWS4qO4dZIj1+/Hgzyb8bU+/vUn5Oc9XqyJVflG9gAJQkwAAoSYABUJIAA6AkAQZASYvf\nQhwpzmV1MLfGdYe5uS7rcPmYvJ9v1oSy1TZN0yuvvNJMXnvttWaS9cV81pUrV5rJmTNnmkkWI6fe\nNr/B5YG7beSO0jnJWmCeiK61tbVmknfrzhfPClm2Gb///e83k1/96ld5APmwrLbm75vnNLdBdt+B\nhd+FmL/g4J24R5Zw5rNeeumlZvKLX/yimfzsZz/LA8itre+9914z+eijj5pJ9hIHd2DOv7n46wMA\n35QAA6AkAQZASQIMgJIEGAAlaSFOU6+T0+1i5T1Sc/FgVv5ycvXq1WaSe/OybjRN061bt5pJbnHM\nYtuhQ4eayXe/+91mkrWo7g1h882cEyO73fLgR66NvXv35kvlXZK/973vNZO8RfJIhzNP34kTJ/IA\n8gzmdsTcfPjnP/+5meQGxe6VvzCr8/4/ebVneXiapsOHDzeTPIMj+1Hzb0J2TX/961/nAeSd00ca\nxfl3QwsRAJ4mAQZASQIMgJIEGAAlCTAASnoWW4gpS1bdTWjZOXz55ZebSd5ZNZtLuTcvD+Ds2bN5\nANkvevvtt5tJFpxyXVt2F3NxX7d4Nrc78XaqJpedwzxZ0zS9+eabzeSdd95pJrl4MCtk6+vrzSRr\ngd1bY2ehMc9gru7M7tnIxshpjs/7bEZukZy3xp56ZzkvmE8++aSZnD9/vpnkm5xLKbPFOvXuoJ3n\nfW6rwrvBNzAAShJgAJQkwAAoSYABUJIAA6AkAQZASYtfo5+tYJ13AZ+maWVlpZnk6tXcCpoV5CzW\nZ/n+448/zgN4//33m8mRI0eaSS5jzTvZ/+53v2smn3766bavM/X6x6UL1nnw2YpeW1vLJ+b5ypW7\ny8vLzSQ3Jp87d66Z/P3vf28mv/3tb/MActdz/tOIP/zhD81kY2OjmWSxvvQJHZSXcZ6sU6dO5RN/\n+tOfbvvEXMyd/xAiT8RvfvObZpL/xGLqbQr+8MMPm8kCr+5NvoEBUJIAA6AkAQZASQIMgJIEGAAl\nLX4LccRg8yrrfNevX28mly9fbia5JzSbirkxtnsn+1zmm4eU9cVc/Pqvf/2rmWQtqrsSNN+o0qW1\nPPiRLatT77znJKtfed7ThQsXmsl7772XD/vLX/7STLI1mruD8yDzLJc+oYPyLOf7kJ+aqVf+zBbi\nyKbsfJ0//vGPzeSvf/1rHkBWHPOc5gEs8Hpf38AAKEmAAVCSAAOgJAEGQEkCDICSnsUWYvascjda\nV96jPWuB2TTLgl8WxnInXjaXpt4NxU+fPt1MPvvss2Yy0rDKdtPgreUH37p5MNI5zFOTi+ymafro\no4+aSfbKUi6czBOaFbJuDTIX5Y38LoNnOS1YNTE/Efm258maem9gft6zR5rXRv4lycm1a9fyAFKW\nXRd482Eq89cHAP4vAQZASQIMgJIEGAAlCTAAStoz2w2LCzl06FAzmblSlY27ffv2NZO8gW/3lr6N\n7CnlrVenb9Ei29bMbbScXL16dUcO6RtZWlpqJiOdw5ErobuUMu/EffDgwW1ffKRzmBWyHfyEzvxS\n+bt0u5G7bffKkHkH9vxoT73Nh/k3ISu+I6Xf+TT/6eAbGAAlCTAAShJgAJQkwAAoSYABUNLitxBX\nV1ebyczr+0Y6eFlnykk+K5ezdRea5ckaOX072N0aqfYN7nDbWd3O2AwG36u8hPIspzzLI7fK3dVP\n6ODvm8fQvWHxbnuSKxm7P2vkT0ee07p/Y+f/yH0DA6AkAQZASQIMgJIEGAAlCTAASlr8OzKPNJcG\n200jLcTZmmYj7b7BJ85mwW65O2jm33r3mmY7uK5wtt2P/Fv3RDxTNzsuwTcwAEoSYACUJMAAKEmA\nAVCSAAOgpGexhThzL3G2lxrplQ0uOdypFtkOVt3m1lPvZz7hPXKzXWYz929hHvgGBkBJAgyAkgQY\nACUJMABKEmAAlCTAAChpz/zfNBoAkm9gAJQkwAAoSYABUJIAA6AkAQZASQIMgJIEGAAlCTAAShJg\nAJQkwAAoSYABUJIAA6AkAQZASQIMgJIEGAAlCTAAShJgAJQkwAAoSYABUJIAA6AkAQZASQIMgJIE\nGAAlCTAAShJgAJQkwAAoSYABUJIAA6AkAQZASQIMgJIEGAAlCTAAShJgAJQkwAAoSYABUJIAA6Ak\nAQZASQIMgJIEGAAlCTAAShJgAJQkwAAoSYABUJIAA6AkAQZASQIMgJIEGAAlCTAAShJgAJQkwAAo\nSYABUJIAA6AkAQZASQIMgJIEGAAlCTAAShJgAJQkwAAo6X8Aiwc/bogrkloAAAAASUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 37, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sel = round(linspace(1,p,6));\n", "clf;\n", "imageplot( mat2cell(w(:,:,sel,3), n, n, ones(6,1)) , '', 2,3);"]}], "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"}]}}}