{"metadata": {"signature": "sha256:68319bb264d7e2249c56d0363ee296fb7893f537fcb29b9821e0100a95a3d661", "name": ""}, "worksheets": [{"metadata": {}, "cells": [{"metadata": {}, "cell_type": "markdown", "source": ["Dijkstra and Fast Marching Algorithms\n", "=====================================\n", "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n", "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n", "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n", "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n", "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n", "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$\n", "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n", "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n", "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n", "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n", "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n", "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n", "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n", "$\\newcommand{\\CC}{\\mathbb{C}}$\n", "$\\newcommand{\\RR}{\\mathbb{R}}$\n", "$\\newcommand{\\EE}{\\mathbb{E}}$\n", "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n", "$\\newcommand{\\Ww}{\\mathcal{W}}$\n", "$\\newcommand{\\Vv}{\\mathcal{V}}$\n", "$\\newcommand{\\Nn}{\\mathcal{N}}$\n", "$\\newcommand{\\NN}{\\mathcal{N}}$\n", "$\\newcommand{\\Hh}{\\mathcal{H}}$\n", "$\\newcommand{\\Bb}{\\mathcal{B}}$\n", "$\\newcommand{\\Ee}{\\mathcal{E}}$\n", "$\\newcommand{\\Cc}{\\mathcal{C}}$\n", "$\\newcommand{\\Gg}{\\mathcal{G}}$\n", "$\\newcommand{\\Ss}{\\mathcal{S}}$\n", "$\\newcommand{\\Pp}{\\mathcal{P}}$\n", "$\\newcommand{\\Ff}{\\mathcal{F}}$\n", "$\\newcommand{\\Xx}{\\mathcal{X}}$\n", "$\\newcommand{\\Mm}{\\mathcal{M}}$\n", "$\\newcommand{\\Ii}{\\mathcal{I}}$\n", "$\\newcommand{\\Dd}{\\mathcal{D}}$\n", "$\\newcommand{\\Ll}{\\mathcal{L}}$\n", "$\\newcommand{\\Tt}{\\mathcal{T}}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\al}{\\alpha}$\n", "$\\newcommand{\\la}{\\lambda}$\n", "$\\newcommand{\\ga}{\\gamma}$\n", "$\\newcommand{\\Ga}{\\Gamma}$\n", "$\\newcommand{\\La}{\\Lambda}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\Si}{\\Sigma}$\n", "$\\newcommand{\\be}{\\beta}$\n", "$\\newcommand{\\de}{\\delta}$\n", "$\\newcommand{\\De}{\\Delta}$\n", "$\\newcommand{\\phi}{\\varphi}$\n", "$\\newcommand{\\th}{\\theta}$\n", "$\\newcommand{\\om}{\\omega}$\n", "$\\newcommand{\\Om}{\\Omega}$\n"]}, {"metadata": {}, "cell_type": "markdown", "source": ["This numerical tours details the implementations\n", "of Dijkstra and Fast Marching algorithms in 2-D.\n", "\n", "\n", "The implementation are performed in Matlab, and are hence quite slow."]}, {"prompt_number": 1, "metadata": {}, "cell_type": "code", "outputs": [{"output_type": "stream", "stream": "stdout", "text": ["Populating the interactive namespace from numpy and matplotlib\n"]}, {"output_type": "stream", "stream": "stderr", "text": ["WARNING: pylab import has clobbered these variables: ['pylab']\n", "`%matplotlib` prevents importing * from pylab and numpy\n"]}], "input": ["from __future__ import division\n", "from nt_toolbox.general import *\n", "from nt_toolbox.signal import *\n", "from nt_toolbox.graph import *\n", "from nt_solutions import fastmarching_0_implementing as solutions\n", "%pylab inline\n", "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload 2"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Installation\n", "------------\n", "*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) for details about how to install the toolboxes."]}, {"metadata": {}, "cell_type": "markdown", "source": ["Navigating on the Grid\n", "----------------------\n", "We use a cartesian grid of size $n \\times n$, and defines operators to navigate in the grid.\n", "\n", "\n", "We use a singe index $i \\in \\{1,\\ldots,n^2\\}$ to index a position on\n", "the 2-D grid.\n", "\n", "\n", "Size of the grid."]}, {"prompt_number": 2, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["n = 40"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["The four displacement vector to go to the four neightbors."]}, {"prompt_number": 3, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["neigh = array([[1, -1, 0, 0], [0, 0,  1, -1]])"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["For simplicity of implementation, we use periodic boundary conditions."]}, {"prompt_number": 4, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["boundary = lambda x: mod(x,n)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["For a given grid index |k|, and a given neighboring index k in ${1,2,3,4}$,\n", "|Neigh(k,i)| gives the corresponding grid neighboring index."]}, {"prompt_number": 5, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["ind2sub1 = lambda k: [int( (k-np.fmod(k,n))/n ), np.fmod(k,n)]\n", "sub2ind1 = lambda u: int( u[0]*n + u[1] )\n", "Neigh = lambda k,i: sub2ind1(boundary(ind2sub1(k) + neigh[:,i]))"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Check that these functions are indeed bijections."]}, {"prompt_number": 6, "metadata": {}, "cell_type": "code", "outputs": [{"output_type": "stream", "stream": "stdout", "text": ["[13, 27]\n", "134\n"]}], "input": ["print( ind2sub1( sub2ind1([13, 27]) ) )\n", "print( sub2ind1( ind2sub1(134) ) )"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Dikstra Algorithm\n", "-----------------\n", "The Dijkstra algorithm compute the geodesic distance on a graph.\n", "We use here a graph whose nodes are the pixels, and whose edges defines the\n", "usual 4-connectity relationship.\n", "\n", "\n", "In the following, we use the notation $i \\sim j$ to indicate that an\n", "index $j$ is a neighbor of $i$ on the graph defined by the discrete\n", "grid.\n", "\n", "\n", "The metric $W(x)$. We use here a constant metric."]}, {"prompt_number": 7, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["W = ones( (n,n) )"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Set $\\Ss = \\{x_0\\}$ of initial points."]}, {"prompt_number": 8, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["x0 = [n/2, n/2]"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Initialize the stack of available indexes."]}, {"prompt_number": 9, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["I = [sub2ind1(x0)]"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Initialize the distance to $+\\infty$, excepted for the boundary conditions."]}, {"prompt_number": 10, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["D = ones( (n,n) ) + inf\n", "u = ind2sub1(I)\n", "D[u[0],u[1]] = 0"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Initialize the state to 0 (unexplored), excepted for the boundary point $\\Ss$ (indexed by |I|)\n", "to $1$ (front)."]}, {"prompt_number": 11, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["S = zeros( (n,n) )\n", "S[u[0],u[1]] = 1"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Define a callbabk to use a 1-D indexing on a 2-D array."]}, {"prompt_number": 12, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["extract   = lambda x,I: x[I]\n", "extract1d = lambda x,I: extract(x.flatten(),I)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["The first step of each iteration of the method is\n", "to pop the from stack the element $i$ with smallest current distance $D_i$."]}, {"prompt_number": 13, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["j = np.argsort( extract1d(D,I)  )\n", "if np.ndim(j)==0:\n", "    j = [j]   # make sure that j is a list a not a singleton\n", "j = j[0]\n", "i = I[j]         \n", "a = I.pop(j)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["We update its state $S$ to be dead (-1)."]}, {"prompt_number": 14, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["u = ind2sub1(i);\n", "S[u[0],u[1]] = -1"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Retrieve the list of the neighbors that are not dead and add to I those that are not yet in it."]}, {"prompt_number": 15, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["J = [] \n", "for k in np.arange(0,4):\n", "    j = Neigh(i,k)\n", "    if extract1d(S,j)!=-1:\n", "        # add to the list of point to update\n", "        J.append(j)           \n", "        if extract1d(S,j)==0:\n", "            # add to the front\n", "            u = ind2sub1(j)\n", "            S[u[0],u[1]] = 1\n", "            I.append(j)"], "collapsed": false, "language": "python"}, {"prompt_number": 16, "metadata": {}, "cell_type": "code", "outputs": [{"metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD9CAYAAACY9xrCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAAvtJREFUeJzt2EEKwjAARUEr3rvpyevShQER1KpvZllKyOYR+KcTAAAA\nAAAAAADwtOVdB48x9nedDTw2xrjr+3zERYBjCB5CBA8hgocQwUOI4CFE8BAieAgRPIQIHkIEDyGC\nhxDBQ4jgIUTwECJ4CBE8hAgeQgQPIYKHEMFDiOAhRPAQIngIETyECB5CBA8hgocQwUOI4CFE8BAi\neAgRPIQIHkIEDyGChxDBQ4jgIUTwECJ4CBE8hAgeQgQPIYKHEMFDiOAhRPAQIngIETyECB5CBA8h\ngocQwUOI4CFE8BAieAgRPIQIHkIEDyGChxDBQ4jgIUTwECJ4CBE8hAgeQgQPIYKHEMFDiOAhRPAQ\nIngIETyECB5CBA8hgocQwUOI4CFE8BAieAgRPIQIHkIEDyGChxDBQ4jgIUTwECJ4CBE8hAgeQgQP\nIYKHEMFDiOAhRPAQIngIETyECB5CBA8hgocQwUOI4CFE8BAieAi5HH0BvtO6rtPv27Z9+Ca8khce\nQgQPIYKHEMFDiOAhxErPdJFflmX6777v0+/W+9/ghYcQwUOI4CFE8BAieAix0jNd2K3x/8kLDyGC\nhxDBQ4jgIcRox5Rx7j954SFE8BAieAgRPIQIHkIEDyGChxDBQ4jgIUTwECJ4CBE8hAgeQgQPIYKH\nEMFDiOAhRPAQIngIETyECB5CBA8hgocQwUOI4CFE8BAieAgRPIQIHkIEDyGChxDBQ4jgIUTwECJ4\nCBE8hAgeQgQPIYKHEMFDiOAhRPAQIngIETyECB5CBA8hgocQwUOI4CFE8BAieAgRPIQIHkIEDyGC\nhxDBQ4jgIUTwECJ4CBE8hAgeQgQPIYKHEMFDiOAhRPAQIngIETyECB5CBA8hgocQwUOI4CFE8BAi\neAgRPIQIHkIEDyGChxDBQ4jgIUTwECJ4CBE8hAgeQgQPIYKHEMFDiOAhRPAQIngIETyECB5CBA8h\ngocQwUOI4CFE8BAieAgRPAAAAAAAAAAAAAAAAAAAAAAAAAAAADdXMhYWTM7DMbMAAAAASUVORK5C\nYII=\n", "text": ["<matplotlib.figure.Figure at 0x10a604f50>"]}], "input": ["imageplot(S)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Update neighbor values.\n", "For each neightbo $j$ of $i$, perform the update,\n", "assuming the length of the edge between $j$ and $k$ is $W_j$.\n", "$$ D_j \\leftarrow \\umin{k \\sim j} D_k + W_j. $$"]}, {"prompt_number": 17, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["DNeigh = lambda D,k: extract1d(D,Neigh(j,k))\n", "for j in J:\n", "    dx = min(DNeigh(D,0), DNeigh(D,1))\n", "    dy = min(DNeigh(D,2), DNeigh(D,3))\n", "    u = ind2sub1(j)\n", "    w = extract1d(W,j);\n", "    D[u[0],u[1]] = min(dx + w, dy + w)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["__Exercise 1__\n", "\n", "Implement the Dijkstra algorithm by iterating these step while the\n", "stack |I| is non empty.\n", "Display from time to time the front that propagates."]}, {"prompt_number": 18, "metadata": {}, "cell_type": "code", "outputs": [{"metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAATkAAAD/CAYAAACKJ6HsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADslJREFUeJzt3c+LJdUZxvHHMBBhwEUgEJwkoGMSHQ0MOpjFCAoRA1kk\nCFlk55/S9J+SXRaBoIuEiAEFXSgzIiSZMVFHSMYgCFkIgoJgFlPn1Hv7vl116nfVW9/PIl2pvt19\nR6+nn/c+VWckAAAAAAAAAAAAAAAAAAAAAAAAAMCM7pvuW598M933xjCnE/573wNe2+t1/Nr+1hJP\nAwDmwiIHIDQWOQChscgBCI1FDkBoLHIAQmORAxAaixyA0C4s/QQALOdEp/n4VCcLPpPpkOQAhMYi\nByA0xlVgh9KY+vzB2ZijK0kOQGgkOWAnbMmQEtwzPzAP+I999Gn1v9tPdCQ5AKGxyAEIjXEVCM4r\nGfKY+ow596Z5QB5dt19GkOQAhEaSAwJqLRlSgrvuf31OdQHKCJIcgNBY5ACExrha4I6J/g9vLKpj\nX0pLhjymnjOu5q8NUEaQ5ACExiIHIDTG1QZpTH3oBXPuVUZXrEuvJrU6/uTqd/KpS/pf489palzv\nHa3zvweSHIDQSHJn2JIhJ7hf1Z9/yD62SnUkOszNS29SecmQEtwNXatPXr2RD5tSnV9GSGu9jo4k\nByA0FjkAoTGuVrySIY+pL/pfk0ZXygjMpfE6OKm4ZEhj6g095f+ganQtLiOk1V5HR5IDENquk1xb\nyZAS3BsPPp1PPfviO0ffhzICUyq+REQqLhlSgrtpiwdPYRkhrfemfpIcgNBY5ACEtstxtbRkSGPq\na/p5ffLB+rBpdKWMwFCdb7Y3x20lw818rmVctTZaRpDkAIS2myTXp2RICe5V/cL/plWqKy0jJFId\nmg25D1UqLxnS5z9744f1uWcLn+TGygiSHIDQWOQAhBZ6XHVHVKm4ZEhj6jv/MDn+cecHFZYREtfR\nwTdkR98+JUMeU1+vv89nmm50XbKMIMkBCC1kkmu8REQqLhlygvtD/aXvqCzVeYlO4hIT1Mba7LJP\nyZAT3OtypVRXnOik7peYzLT5JkkOQGgscgBCCzOull4HJ3UoGdKY+nv/Z+bRlTIChabY0bdPydA2\nruavnb2MkMa+jo4kByC0zSe5rvehSuUlQ05w739sztkMVn1tYRkhcb/rXk252WWvkiEfv23O/ezc\n5y/NVEZIo19iQpIDEBqLHIDQNjmuDrnZXupQMuQx9Xfm3Ev1cdPo6o2tUueb+hlbY5hyR99+JUMa\nU//sP+GG0XXKMkIa/6Z+khyA0FjkAIS2yXHVjnBprDseHM+MhA86DzAj5UFDmqRx1I6oj5qf9Nvq\n42/qU08//oYk6QX9JZ97Xn+tn9N/q+f0R/NzXrn34eNX61OMqbG8lg7MCHbQKDryWGdGvTZpfLQj\npcsdR8255858NMffffbf+dQ11c/tqer4mm4eff7Se2ZEfevMR0ky/yzerP4ZvWY+TbsKAOfYZJKz\nUiXwkklBQ1JdY6KT6vQm5QSX0ptUJzg3vUl1gnulPpUSnKk3EEydROrSrDTVHbxRX5jqbCHQmOpe\nd9KbPTbnUoLz0tu98zePPp8TnE1t6dhJb1Kd4LjjAQAKsMgBCG3z46o3BqTR1RtbJTO69ikjCksG\nd0SV3JIhjalL/i3jmMfhv+Pj0XXRMuK542OvZPBGVPv50pLBG1El9pMDgE42n+QS7zfkFGVEccng\npDfJLxlIcPvkproly4jn6sOmksFLb1L3kmHK9GaR5ACExiIHILQw46rVVEZILaNrSxlRXDI4I6pE\nyQBf03V0xWWE1Hl0tWNracngjqhS55Jhrv8GSHIAQguZ5JLz6vqmS0zayog+96FSMqDUkDJC6n6J\niS0jBt2Hao8XLBk8JDkAobHIAQgt9LhqDbqOzv5FND1utmdERR9zlxGDbraXVlEyeEhyAELbTZKz\nut7velBGcB8qZjZXGbG2zS7HQpIDEBqLHIDQdjmuJn3KCEoGLGnKHYbXtqPvWEhyAEJjkQMQ2q7H\nVau0caVJxRpMscPw2nb0HQtJDkBoJLkz2soISgaszWg7DG/gZvs+SHIAQmORAxAa42oDr4zYQjzH\nfg25qX8LN9v3QZIDEBpJrsDWfnMBvcqIACWDhyQHIDQWOQChMa4CwZWWERFKBg9JDkBoJDlgJ9rK\niAglg4ckByA0FjkAoTGuAju0p7t5SHIAQiPJATsWNb1ZJDkAobHIAQiNRQ5AaCxyAEJjkQMQGosc\ngNBY5AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABg5+6b7luffDPd98YwpxP+\ne9+Dbb22T3QqSXrSnEvH3/9Rfe7uB/Xxu2c+StKpTiZ4dmM7fm3zt3UBCO3C0k9gC/5U/SaUpF9u\n4rcZ9u7EvGafPPNRMgnuUXPOfoMPdCR9z20kuhpJDkBoJLkGKcFd/7Y59xWpDutV+v5bTnCP+d8n\np7qGRCdtI9WR5ACExiIHIDTG1TNsyZDG1Aeu1p+//p55bDW6MrZiSX1KhjymnjOu5q+1/2ejoytJ\nDkBoJLmKVzLkBGeS3APma1Kqo4zA3Lz0Zo9bS4bCJGc1lRH2Oa0t0ZHkAITGIgcgtF2Pq20lQx5T\nf+p/fRpdKSMwN29ElTqUDNXxl0/Up+4v/NltZcTakOQAhLbLJFdaMuQEZ885KCMwl5edOxr6lAwp\nwX108eF87vITd/LxkFT3spmQfr2C1z5JDkBoLHIAQtvNuNqrZKjOfWke1xbjm8oIidEV3b3sXBPn\njqhSccmQxtQP9Uh98mJ9mEbX0rFV8q+jS899ybGVJAcgNBY5AKGFHle9EVUqb1LTmHrrYp39r1y9\nnY+borzXuEpcR4dyjU2qN6La45YmNY2pH+my/8Or0TVC40qSAxBayCTXeB2cVFwypAR3S1fqk+bN\n2ZTqSssIievo0Ky4ZPDSmzluKxlSgjsoHjwBygiSHIDQWOQAhBZmXC2+Dk4qLhnSmPpP/dj/oVWU\nLy0jJG7qh69zyXDOuFpaMrQWD56NlhEkOQChbT7Jdb7Z3pxvKxlSgrttiwcPZQR6GFQyOOlNKi8Z\niosHj1NGSD1SnVNGSOOnOpIcgNBY5ACEtslxdcjN9lJ5yZDG1Ftt46pFGYEG3ogqDdvRt0/JkM7d\n+cQUD5dK/gRnDLiO7rwdhse+jo4kByC0TSW5ITv69ikZ0udvf2KSXOlvO6eMkLrf70oZEcOUO/r2\nKRlygnuvfkXe0Tipbm2XmJDkAITGIgcgtNWPq2Pt6NunZMhj6o37jx4nadDoyg7D8c21o2+vkiGN\nqX+3z9gZXWcuI6Txb+onyQEIjUUOQGirHFen2NG3V5OaxlQzMtrgnb++S6TveB0dOwxvz9w7+vZp\nUvOYejCuWvceG6FxJckBCG1VSW7KHX37lAw5wd047xnff/D9JM1SRkhcR7c2S+7o26tkaE1yyfbL\nCJIcgNBY5ACEtvi4OteOvr1KhhtnPp5r3jJC4qb+tVjDjr69SobicdXaQhlxjCQHILTFktzcO/r2\nKhnS8aefm3P27X8PZUR0a9vRt1fJkI6/vmvOHeSkBiPd1D/BDsMekhyA0FjkAIQ267i65I6+vUqG\nPKa+Zc5dN1/TNLqOc1M/ZcT6rG5H3z4lQx5T3zXn7NcMGF1XssNwQpIDEBqLHIDQZhlXV7FteZ8m\nNY+pb8uVRtfCxvXgeS66jTqGWNu25b2a1DymmnHVSqNr8dgqreE6OsZVALszS5JLb3rbNJESRnkG\nMknmovvQY+a3yMGb/02/F2wqsyVDZs59r3rsNfPpdGwT6rUv8+Fjl25Jkq7oVn2uOv6J/pXPpc9f\n+cKkN1uUpOO/HZ/73Dzura/ufTwsHk6F/u6atND17XmbTopfx1b1mj5IS6U5x6ayr70HmErlQvVY\nU5TkY++cJF299zp/+NJH+dQj+lCSdFll5yTp8hdV8WDT6O0zH+3x+8d/EoskByA0FjkAoc16ndzB\n2/fVGGWv6SodXe0b8MNG15aYn0bXT50RVapHUzuupjHVGVGlegx9zIyraUy1I2waU90RVarHVHMu\njalpRJXOrUwwwMFb9dXo2v2qsoGjq3k9Dxpd7dh6wfwpmkZTZ0SV6jH1YPRU4bkvTMmQxlRvNLXn\nqjH1Lrd1AdizWZPcqXnz+yS9+W1SxyxlhJR/C/YqIwpLBi+9ST1KBi+9mfNeyWDT2yl3OozOvehi\nQBkhmVQ3oIyQbKrrUUaUlgxOepO6lwxuepOKS4aU4M65CCYjyQEIjUUOQGiL7SeXxqgTe83WFsqI\nwpLBG1GlHiWDM6JKzSUDI+q03LddrI5lhLSS6+gKSwZvRJW6lwzuiGqPW0qGlns2MpIcgNAW/zse\nSssIqTnVzVZGFJYMXnqzny8uGZz0JlEyrIU7kSRbKyMKSwYvqdnzxSWDl97scUvJkI5PW+7mIckB\nCI1FDkBoi4+rVlMZIXW/jm6SMmLAzfZS95LBG1ElSoa1CVFGDLjZ3p4fdLO9VFwylL72SXIAQmOR\nAxDaqsbV5Nzo3/E6OrdxlTqPrrZxHXKzvdS9ST3vZnvG1PXabOM64GZ7e37IzfZSlya1DEkOQGir\nTHLWKm7qN78B597Rl/S2XaVlhLSOHYaX3NF3rJLBQ5IDEBqLHIDQVj+uWmu4qX/uHX0ZUWNoLCOk\nVewwvOSOvmOVDB6SHIDQNpXkkiXLCHb0xRBTlhHSsEtMltzRd6ySwUOSAxAaixyA0DY5rlpzlxHs\n6IuxrO2m/iV39J3ytU+SAxDa5pNcMtcOw+zoiyms4X7XZXf0nQ5JDkBoLHIAQgszrlpT7jDMjr6Y\n0pJlxNp29B0LSQ5AaCxyAEILOa4mU+wwzI6+mMvcjevadvQdC0kOQGihk5w11k393GyPuc21w/Da\ndvQdC0kOQGgscgBC2824ag26qZ+b7bGgKXcYXtuOvmMhyQEIbZdJLulTRlAyYA2mKCPWtqPvWEhy\nAEJjkQMQ2q7HVau0jKBkwNqMdlP/Bm6274MkByA0ktwZbWUEJQPWbND9rhu4D7UPkhyA0FjkAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAENf/ARmJQ5FDbouVAAAAAElFTkSu\nQmCC\n", "text": ["<matplotlib.figure.Figure at 0x10a860ad0>"]}], "input": ["D = solutions.exo1(x0,W)"], "collapsed": false, "language": "python"}, {"prompt_number": 19, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["## Insert your code here."], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Display the geodesic distance map using a cosine modulation to make the\n", "level set appears more clearly."]}, {"prompt_number": 20, "metadata": {}, "cell_type": "code", "outputs": [{"metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD9CAYAAACY9xrCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACRJJREFUeJzt3VGO3cYRBdB6QZYQ7yRenITI8cSCvTh5J9Ieks8EYTeM\nBlmt9txzPokHukTrgkCpWF0FAAAAAAAAAAAse3Xd+J9V/x5d//nv12ufvnwY3uPzD79eL377efwf\n/Nv4+oevP12u/frjL8Pf/vz7+Najmj9++TT87W8//Ot6cbHmf3z9eLn2y4+/jWtbqLnzOf8Za37v\nfzdev1/z/Zfxfw14jwQeggg8BBF4CPLXrhuPmhlV4ybMsAFTNW5oLDRgqsZNmJUGTNW4CTNswFQt\n1TxqwFSNmzCrNXc951mTaXjvr+NbrDTzZr8d3Xvpz1eT/4eTmleaebPfju699Peo1p5zva7Pzhse\nggg8BBF4CCLwEKStadc1IbXSnKsaN1VWmnNV9yekVppzVWs1dz3nWc0rjbHpbxeaebNm5fDZTe67\nUvO0ibbQzJvVPPw7OrnvSjNv+ucrTTuIJvAQROAhiMBDEIGHIG1d+q5x2e3fK1fdHpfd/l14VVvN\nS53wlfHQyb27xnCndayM4U7u3TWGO61j9q8vA97wEETgIYjAQxCBhyBtTbuucdnt37JX3R6X3f4t\ne1VbzZ3fp9/9pv571Hz3m/qlMdzJved/n98ul7zhIYjAQxCBhyACD0EEHoL0dembxmV3L6+ouj8u\nu3t5RVVfzbsXUnSN4U5/v3mJxuqo+NISjcFBct7wEETgIYjAQxCBhyACD0HauvRd8/G7l1dU3Z+P\n3728oup+zbPnvNJV3r5EY/NZdtM6DjnLbvQn8YaHIAIPQQQeggg8BGlr2nWNy+5eXlF1f1x29/KK\nqsbn/MByh66FFNvPspvc+5Sz7Ea84SGIwEMQgYcgAg9BBB6CtHXp28ZlNy+CqHpgXPaQmp94zp3L\nHW4vpNh8lt3094ecZfc2qNkbHoIIPAQReAgi8BCkr2nXNS67+bvwqgfGZb9DzV3Pefe33iefZVe1\nv+aV51yv63P2hocgAg9BBB6CCDwEEXgI0tal7xqX3b28our+uOzu5RVVfc/5lOUOJ5xlN6vjmLPs\nSpceogk8BBF4CCLwEETgIUhbl75rPn778oqqI2refabe9jPWJvc++iy7yb1POctuxBseggg8BBF4\nCCLwEKStadc1erp9eUXV9ppPOFNv9xlrVfeXaHyPmu8u0eit+e1yyRseggg8BBF4CCLwEETgIUhf\nl75p9HT38oqqvppPPlNvVvMxyx1OOMuu6vYSjdaz7F7Xa97wEETgIYjAQxCBhyBtTbuu0dPd37JX\n3a9597fsVX01n/J9+gln2VU98E1941l2nwfXvOEhiMBDEIGHIAIPQQQegrR16bvGZXcvr6i6Py67\ne3lFVWPNBy+k+DPW3HmW3Yg3PAQReAgi8BBE4CGIwEOQti5923z85kUQVQ/Mx7+jmo9eSPGOan5i\nicbboGZveAgi8BBE4CGIwEOQvqZd17js5kUQVQ+My76jmk9eSLH7LLuq/TWvLNGo17UOb3gIIvAQ\nROAhiMBDEIGHIG1d+q5x2d2LIKruj8u+p5qPXkix+Sy7qjNqnv62dOkhmsBDEIGHIAIPQdqadl3j\nstu/C69S8/86+Pt0Nf8xb3gIIvAQROAhiMBDEIGHIG1d+q7R0+2LIKrU/Ae/raojFlKo+f98e7tc\n8oaHIAIPQQQeggg8BBF4CNLXpV9Yc/zEOVqzbudwrnnSRV3phKfW3HkO4IevP12urcylr6xGrwqo\necAbHoIIPAQReAgi8BCkrWm3tPX0gXO0Vhpjs6bKUmMstOYnFpuMGl1Va6Ono2bXUqOr6t3XbLQW\nwgk8BBF4CCLwEETgIUhbl75r6cDKeOisjqWR1sm9u0Zaqw6veWH0dKWzXXV/9DS15um/iL2u17zh\nIYjAQxCBhyACD0HamnZL2zgfOEfr6O/T31HNXd+FVz0wehpQ88o5gJ8H17zhIYjAQxCBhyACD0EE\nHoL0dekn44Wjbmfn2V8nLKRYPRfu5Jq7FkFUPTB6unl5RVVfzU+cAzjiDQ9BBB6CCDwEEXgIIvAQ\npK1LP5snXuky3p67rzpjIcXmc+Gq+mo+5ty0A5ZXVN2veaUbX7VW89vgt97wEETgIYjAQxCBhyBt\nTbuu5Q4rY7jTOjYvpNh+Ltzk3k/UfMq5aUcsr6i6PS67uthkpeZ6Xe/tDQ9BBB6CCDwEEXgIIvAQ\npK1L37XcYWkMd3Lv3Us0dp8LV9VX88lnvZ1S8xPLKx6puXTpIZrAQxCBhyACD0H6mnaT8cK733qv\njiKe8E399nPhqtpqPvmst93fslfdH5dtrXnAGx6CCDwEEXgIIvAQROAhSFuXfjYSeXu5wwNnrG1f\norH5XLiqvud88llvu5dXVN0fl+2sub69XS55w0MQgYcgAg9BBB6CCDwEaevSdy136DxjrWuJxvZz\n4Sb3fuQ5H3zW2+7lFVVn1Dz9F7HX9Zo3PAQReAgi8BBE4CFIW9OubbnD5jPWpr8/+Cy7qr7nfPRZ\nb5uXV1Ttr3nlOY/+JN7wEETgIYjAQxCBhyACD0H6uvST8cK7yx12n7FWdX+Jxu6z7Kr6nvPRZ71t\nXl5R1VfzE895xBseggg8BBF4CCLwEKStaffxy6fh9ZWmwwlnrFU98E395rPsZnU88ZxPPutt97fs\nVfdr7twZ8Db4rTc8BBF4CCLwEETgIYjAQ5C2Lv20az7o/q50JLefsTa599Fn2U3u/cQY7innph2x\nvKLq9rjsE8959i9i9bre2xseggg8BBF4CCLwEETgIUhbl37WVR52zRfmtnefsVZ1fyHF7rPspr9/\noOaTz3o7pebOJSGjjvz073Pp0kM0gYcgAg9BBB6C9DXtJuOFo+ZH1xju9Pebl2hsP8uuqq3mk896\n2728our+uOxKc65q8RzAAW94CCLwEETgIYjAQxCBhyBtXfrZSORKZ/TuGG7VIUs0Np9lN63jgZpP\nPutt9/KKqvvjskvd+Kqlmuvb2+WSNzwEEXgIIvAQROAhSFvTrutb75Ux3Kozvqnffpbd5N6PjOEe\nfNbb7m/Zq+6Py67uZlh5zp9f12ve8BBE4CGIwEMQgYcgAg8AAAAAAAAAAAAAAAAAAAAAAAAAAMB/\n/Qcr6ab1f7vvAwAAAABJRU5ErkJggg==\n", "text": ["<matplotlib.figure.Figure at 0x10ad19e10>"]}], "input": ["displ = lambda D: cos(2*pi*5*D/max(D.flatten()) )\n", "imageplot(displ(D))\n", "set_cmap('jet')"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Fast Marching\n", "-------------\n", "The Dijstra algorithm suffers from a strong metrization problem, and it\n", "actually computes the $\\ell^1$ distance on the grid.\n", "\n", "\n", "The Fast Marching algorithm replace the graph update by a local\n", "resolution of the Eikonal equation. This reduces significantly the grid\n", "bias, and can be shown to converge to the underlying geodesic distance\n", "when the grid step size tends to zero.\n", "\n", "\n", "Over a continuous domain, the distance map $D(x)$ to a set of seed\n", "points $ \\Ss $ is the unique solution in the viscosity sense\n", "$$ \\forall x \\notin \\Ss, \\quad \\norm{\\nabla D(x)} = W(x)\n", "      \\qandq\n", "    \\forall y \\in \\Ss, \\quad D(y) = 0. $$\n", "\n", "\n", "The equation is then discretized on a grid of $n \\times n$ pixel, and a\n", "solution $ (D_{k,\\ell})_{k,\\ell=1}^n \\in \\RR^{n \\times n} $ is found by using\n", "an upwind finite difference approximation, that is faithful to the\n", "viscosity solution\n", "$$ \\forall (k,\\ell) \\notin \\tilde \\Ss,  \\quad\n", "     \\norm{ (\\nabla D)_{k,\\ell} } = W_{k,\\ell}$\n", "  \\qandq\n", "  \\forall (k,\\ell) \\notin \\tilde \\Ss, \\quad D_{k,\\ell}=0,\n", "$$\n", "where $\\tilde \\Ss$ is the set of discrete starting points (defined here\n", "by |x0|).\n", "\n", "\n", "To be consisten with the viscosity solution, one needs to use a\n", "non-linear upwind gradient derivative. This corresponds to computing\n", "the norm of the gradient as\n", "$$ \\norm{ (\\nabla D)_{k,\\ell} }^2 =\n", "      \\max( D_{k+1,\\ell}-D_{k,\\ell}, D_{k-1,\\ell}-D_{k,\\ell}, 0 )^2 +\n", "      \\max( D_{k,\\ell+1}-D_{k,\\ell}, D_{k,\\ell-1}-D_{k,\\ell}, 0 )^2.\n", "$$\n", "\n", "\n", "A each step of the FM propagation, one update $ D_{k,\\ell} \\leftarrow d $\n", "by solving the eikonal equation with respect to $D_{k,\\ell}$ alone.\n", "This is equivalent to solving the quadratic equation\n", "$$ (d-d_x)^2 + (d-d_y)^2 = w^2 \\qwhereq w=W_{k,\\ell}. $$\n", "and where\n", "$$ d_x = \\min(D_{k+1,\\ell},D_{k-1,\\ell}) \\qandq\n", "   d_y = \\min(D_{k,\\ell+1},D_{k,\\ell-1}). $$\n", "\n", "\n", "The update is thus defined as\n", "$$\n", "      d = \\choice{\n", "          \\frac{d_x+d_y+ \\sqrt{\\De}}{2} \\quad\\text{when}\\quad \\De \\geq 0, \\\\\n", "          \\min(d_x,d_y)+w \\quad \\text{otherwise.}$\n", "      }$\n", "      \\qwhereq\n", "      \\De = 2 w^2 - (d_x-d_y)^2.\n", "$$\n", "\n", "\n", "Note that in the case where $\\De<0$, one has to use the Dijkstra\n", "update.\n", "\n", "\n", "Once the Dijstra algorithm is implemented, the implementation of the Fast\n", "Marching is trivial. It just corresponds to replacing the graph udpate"]}, {"prompt_number": 21, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["D[u[0],u[1]] = min(dx + w, dy + w)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["by the eikonal update."]}, {"prompt_number": 22, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["Delta = 2*w - (dx-dy)**2\n", "if (Delta>=0):\n", "    D[u[0],u[1]] = (dx + dy + sqrt(Delta))/ 2\n", "else:\n", "    D[u[0],u[1]] = min(dx + w, dy + w)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["__Exercise 2__\n", "\n", "Implement the Fast Marching algorithm.\n", "Display from time to time the front that propagates."]}, {"prompt_number": 23, "metadata": {}, "cell_type": "code", "outputs": [{"metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAATkAAAD/CAYAAACKJ6HsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnU1sXNd1x89EpEiR1pCekWcoUrKHDGl9WFJSy4iCpAWC\nNgiycpYJEGSbhTcGmo0XBggBWWiTAt544W02XhQFnIURGE1RFEnQGGrrWmosRYo8jiiaHJs0PyRL\n1IfZxTvn3v+bOSMy5tfw8v/b3Ms7X28e77z3P/fcc44IIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC\nCCGEEEIIIYSQbaSwdW89ubp17002xvkt/L/vBTi3O5fWuf2VnTgMQgjZLniRI4QkDS9yhJCk4UWO\nEJI0vMgRQpKGFzlCSNLwIkcISRpe5AghScOLHCEkaXiRI4QkDS9yhJCk4UWOEJI0vMgRQpKGFzlC\nSNLwIkcISRpe5AghSdO10wdACNlaJuW8iIhUYaykbRHGlqA/r+0sjJ2XyU0/tu2ASo4QkjRUcoQk\nxAVVbeMwVtP2CIxVVMoV+uPY6p3Yb6iUm4LXnNb3vg5jr+wCdUclRwhJGl7kCCFJQ3OVkF3Om2pG\nioic1nZiID7efVQ7FXhRWVt4XmEx9qtz2jbi2JmbWXsNnleDz/5Rh5quVHKEkKThRY4QkjQ0V9fB\nayDJX+5QSU72Bj/XufgNGDvXE/vFCe2MwhPGtB2GMdsoB95VAe9q2Cg3HYe6b2TtyQ/j2JFrsf/O\nSnZs78LbvNoBvxcqOUJI0lDJKW/oHfIMjI3vy9pyOY79eC6quuuPsvZ9eM1PO+DORdLiAlgS39H2\nG6X4eOE0PPm4thMw5ik5dUKstnE8iDkcQMmJOTDgfYrw2/juJR2bj2Ndeuw7uZ+OSo4QkjS8yBFC\nkmZPm6sXwQw4a3uITsATbPEWJHl5Dvq6AHvugzj2fCN7zxdotpINYk6G78DYOTMVvw6Dp6Bvpuvx\nOLSq5up0Kdq4sxquvywHw9jB0nLoV0ez0Pzh+Wh7FsxcxUh/NHfViXHuPXgczd0dgkqOEJI0e0bJ\nvQGq7UVtq+fgCS9oi4u4puRwwRbvTOZKvxSHzl7M2pk/xM/7FbyEjgmyXmybCDoZgoJ7wRkTCapu\nBkIe6hqifzN4DtooOQElp0mWjpZuhrFaqS4iIkMD4KHALSi9zne4l7XvzMffw/e2+TdAJUcISRpe\n5AghSZO0uYom6g9xV/g/aOfb8ORvZs0qSP/rpSwD1yxENldPxYjl8fks21YB9ySpaVsFZ8UPfwOP\n665wmq3EA4PtLZIhtw/OnAxoooLpeuPokIiI/EmOhbGr2kdzdVon6oIMhrFBWQj9YV2XmYa1Gnvu\nsxNXw9hY70zrl7gXuwWNojj3n3HsTf0NbFdAP5UcISRpklZyL0I/qDcRke9pC2MfnXpKREQugefB\n7oZ4NxsGz8OzpeyOdvrvo+fhmconWQcWZDGP/otvZ+1P1/UNyF7BohpQtIU4VNgOEp4A20ZMvYmI\nXNKYHZzHpuTqIUdwVHVzd6LJUe6P+6OOSuZwmIX9IuakWBEwi6I4lLF7quowBlZ9FEXYenX6ctZi\nJMdWRkRQyRFCkoYXOUJI0iRprlokQ24fHDoZ1Ey9cuqZMPR7+ZaIiPwPrOia9McFW5PxIiKnNTS/\nAZL+W6d+LyIix+Wj+Hkg3y3j6kXYR8foCGKFZzCjb9iniY4tNV1xHxw6GcxMRXP1j3JSRESuzsXn\nPbiiiyhQqeb2kadCf/q4OibK0TFxX/a3HHePrIR+38RdEREZWoR9dFbTEPaXTuhPaBwTAmwhVHKE\nkKRJRsm94cWh4q7wb8auORlMvYmI/If8XcvYtatfyzqX42tvnHou9KePHRaR/K5x48Cpz0P/mcYn\n8QG9o52FxINvNLitZC8yCXO2pm03LOSHFEljccjiUNGJcNXZLmLqTUTk8kdqnVzsjm9kc7oOnxff\nUh7MZErv8guwV0UNn/1yPwxhlIRtQamO/W8YK5jhEw0g6da5XwMlh+dis4tYU8kRQpKGFzlCSNIk\nY65iRt+QLgk2HWEkgy3KopPBzNRrv/tafOKvtb0I7w0m8LXv63PBqWHyvRJWXEWe/vq/hX7BcuKD\nCXwGyr6RvQNmLArV7bFs4HBTKzFdEjrDsG9mLDoZgpn6W3hvm9NXYAz344Xgh2jiXn0ie8/BcoyM\nqMI8tz2kGNQ/Mjzf8h3sOx6BITwXmw2VHCEkaZJRclaPQUSi6x0qFlkcqkh0uV8C/RecDL+Or5F/\n1vbKUhyrY/yCvvZQVH/Dxz7Wj66HsWOlP4X+xOhUy7GN2x32Uctbk4TBDEoV+6PsPAHUnUUgYCQC\nRuSYqgtbRESi1YAWSVB1UKHktzl7KCPuIJEHQ9l73vy2v6XKO7aRynz+u4iE71iBsRLUhdhsqOQI\nIUnDixwhJGl2vblqhZ+xbGCQ/LDYiemSTN7jgq0r6YOZ+k8w9o+xf1FNAnBq3Dx2NPcZzZ89MTyV\nP0Y49tcaLGK9l8CFD6uPgDUTLMkDlg20PZm4NxPTJYWAe4hkCCsn6GQIZuq/+Ad3RU3XcRibavoM\nEVnoj5/tHZsdewEzCDtjRZqrhBDy5eBFjhCSNLveXDWzDivbh7KBEBSMGX3Dfh7wDIVwLQwFM08q\nmqjHwciw50JuL3tPzDtXFdgIZ8OQX2tO+1uo2EkHAj57WdUkDrkq9s6YlQ30wqlEYk44DLYP4Vq4\nD87zpOJuU3tuDR4+kv+M5s+2Y8JjC8fu5JhbhTE8F5sNlRwhJGl2vZIzrsMeMyv6LBAEb/UYRCCj\nL+wRsmD7EMWAXHTUm4jI97Nm4lgMSLb3fFZiHnz8bHGOzY59swOTSWeDyr2hf1TnnCeAIWBFn71I\nA5FoSViqJJEYbA+iK3LFUW8irpXSfXwp9xnNn23HhMcWjh2/rH7HBoxtpRVDJUcISRpe5AghSZOM\nuQrBKXLuA+1AZXssG2iFZxpOkQ4Mtg/hWk6FcpFopn5Lfh/G/kbeyz4DPrzwHrzehj+IQ3jsZO8A\nRl3Y1lbFZA3TTa2IDM9nhh0Gwef3ZGZzGjP6xpxwkE/OHsZ9cDXo2zx/4UEYOla+qk+rhzE0XYPT\nbR6MT+c7mAmLW/nwXGw2VHKEkKRJRslhVt3nNXLgLEYvQPSDlQ20egwIur8t2N6iGET8Gg+m3kSi\nqnvmMmQDhsK6FlHxX3DHZkbgvQk6mk5r5M4ZzKB7QzsQmFPQfq1UD2MY8WAWSa4eg2b0tVRJIjHY\nPienIPeRORlMvYmInJQ/ZmPgVMO+KbzCDYncaGpF5IF+xzo8bSudblRyhJCk4UWOEJI0yZiriJX4\nm4Gyf1UM4NfAYCwbaIVnMKOv5YTDhV3cF2R74dDJEMzU38Dn/S52Z/+QP0ZCRESua3sNohtO2l5K\nzKqrvrKhgfjEZyeiyZirbq9Y4RnM6Gs54TDYHiMZbFkGnQzHnPmO+0GH7OAxEYBlwoZ9ofa067I9\nUMkRQpImSSVn/Ar6PwRlFeIXIHbOygZiPQbL6IupkjAO1SIZcltEzMkA6m0JPhuPiRDjFVX2NSjN\nd0RVUBGtEEu7BGmKxnpn4h/qmMCiz+ZMw0gEU2qYKgnjUL34blNyqN7GbsJnW7oy2Lplqm7pWhyy\nh1/ZJmuGSo4QkjS8yBFCkiZpczW3/2wlmgEvvp21uWBo9ScUQFZb0ZmQzReeJyJxMRXlue6DMweD\nSN5E5Z448jh+BPPjHZ2z38XIHTNTe/3Xj93LzMe+ibthzMxQdJpZZARm8cU9ombaorlqTogh9I5A\nac2wXRTGVvXY/xCt59x33A6o5AghSZO0kkNQQf1U24uwxeSsqTK8M1nZQFz4RfVnr4E4VItk4BYR\nslHe1RbrH5x7z3niPeirM21oMaqt6lgWY43xrutVchiHGiIZcIsIWjH224FjfHc+/112Aio5QkjS\n8CJHCEmaPWOueqBJ+YYG9Z+BwHmrbI/lDufAXLWMvpgqiY4Fslm8qnOpC/bOme/gG2CiFpz6CZi7\nqKBW6shwND2tsj2WO8zVl7DfATrazFwF5xyaruZkeBfM63/X9tUd/F1QyRFCkoYXOUJI0uxpcxXx\nzMzJR5mZUAITFgtusPAM2Q4w/OnnarouoccV8hUWnXKcYasbBvqXsiZX2R7NXnt/z1yFYHsM17K9\ncOhJ3Ukz1aCSI4QkTWHr3npydevem2yM81v4f98LdNbcfhMcE6e1nQCHQrdlFo55JuLeT3ieoOPB\nFCFYMZbRFwMecJvcdkcy+LTObSo5QkjS8CJHCEkaOh4I2eWgmXhBTddxMClr2oc6NVJxHA+r4Hiw\n6vZY56auLWb03a6ccBuBSo4QkjRUcoQkhKesJlXdVWGspEoNg/+X4HEbxqLPu3XLFJUcISRpeJEj\nhCQNzVVCEme3mpmbBZUcISRpeJEjhCQNL3KEkKThRY4QkjS8yBFCkoYXOUJI0vAiRwhJGl7kCCFJ\nw4scISRpeJEjhCQNL3KEkKThRY4QkjS8yBFCkoYXOUIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC\nCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEII2eMUtu6tJ1e37r3Jxji/hf/3vQDn\ndufSOrdZrYsQkjRdO30AhJCdY1LOh/55mdzBI9k6qOQIIUlDJUdIglwAhVbStgiPWx8vAG/Da5aa\nWhGReW1f2WWKj0qOEJI0vMgRQpKG5iohu5zXwcwc0bYEj7vmak/WdsEV4OHD2F9a0RZeY+bqW/B5\nt+DxlzrUjKWSI4QkDZWcYq50vNsd1PYAjN2F/rK2eLdL1Q1POgNTbTUYG4F+VdsKSLnCgHYG4Imq\n5KQXxu7FblGVnCzGsVXtN+bj2Cy83BwXdRjrBHVHJUcISRpe5AghSbMnzdXXVFZXYSzIfBgz0/XA\nvjh291Hsm5nagNec0PdGGf9yB0h2svu44Jim49rWYE6WcSLbBEbPg2eu9muLVwBwPMgdbcFcLWi/\nCuZqFSb/iE76CvxG3nRM2O3eZ0clRwhJmqSV3C/A1T0O4zVtR/BuWNYO3gGf0LYHxlagfztrTsKd\nbW4ua2/B3Qxd7te1/RnVHXF4DebKCW1zc1fnZ+FpGKw4fVR3ZpJ4Sq6N48FTcsF0QTMFlFxZP7v0\nlzhW0t8G/qzsO26XhUMlRwhJGl7kCCFJk6S5+kuVwydgbAJMzuKwdjyZX4ax/qZWJMp47M/FoXIj\n34qI1KZjf2Qlf4wiIj+h6brned2Zs9avjsKgmanDMIb9w9riPDZbcaPmqi3LHIaxj6Gv87wAnzOq\npmvfh3GsW1uM1NjK/XRUcoSQpOFFjhCSNMmYq+jBPK3tKMr4o9D3JL9nrppXqp2kN28TmKvB2wQm\nahG8TWdvZm0JHrdjvwRv8ypN2OR53ZmzZ3BZxczUMXiRzWP0rnqmKy7F6Jx+ADGLywPZpL4PWwf2\nw9aBg4vZRO/GmEWb57gx1DOLMTZSzeIqfK8Darp2wU6FrTRdqeQIIUmz65XcO3oHOANj1QntTMDg\nWou3laZWJC7UtnM82KIs3tkcJefdVUfhPfuuZS0mArDd7rstCyt5PN4+OJGo4Io4T23+opKzx1HJ\noZWic3qp2h2GPt2XTbrlkHIi9tsquYHlXCsicmg4k3LF2Qfx87zICvy9mBUEe1JN6J0BZ8RDUHWb\nvY+OSo4QkjS8yBFCkmZXmqvvgOR/XmVw2dtgtJa5CpL/ge79mRuIm3xM0n8OhmQfZJQ7qBnlyotx\nM1G37RsCZ8Nakt4WZZ//AJ4HYWFk93PhMfvgRMBMxTn7rLaeCQtz9x6Yq9P9QyIi0oC4rgUZFBGR\nz7QVEbkd5nZfGOuTz0P/CZ3bT8pCGGvsy/qV4RjXNTwwE/q9FgaJ5mpXUwugf+LEldg3YxiL8Wxk\n2YZKjhCSNLtKyb3lOBmCgjsNg8e1xbsi9Fd1IXe6FCXWrN75GuAlWF7jbmdKrjIQPQ/VgewuN1yJ\nUfsFVHJ2l8Ogf/0voDf+zOWsxa0xP6ATYtdS0xaD7XORDGNNrUhUcM48nh+N+5o+Bg/atPanYWxO\nZxYqOVN37RwPg6rgUMmVdQ8JOjDu9sffxuHjmbet1IP7rBwspRNYK1VwPIyrQwLyXmwIKjlCSNLw\nIkcISZqON1d/6ewKr6J8N3P1uDMGK7t3JuL1vN5TExGRm7DByOT9p2A0LsiTIiKyHBLLiRy0JHIi\nMiifiYjIIQh5GNYNcgulm2Gs1l8P/f7eL7KOd+YhM6vJ99PX4hiD+ncXrzv5DGu4dOHtdfOcDDDf\nPxnN5iLOXa8/C8su5oSYg7ltputdWIo5AEsxZqaWYW5XNJEcmqu4lLOipu+j0Tj3n4LfS8DmuZcQ\nQERq6sebB3vVzuWXiYagkiOEJE1HKrlftNkVHmJRPYeCo+6WTsRd39f3xSXfP8tXRcRXct6CbTvH\ng93lhiG8wRZ08W630hMXd8dPZLmBiwK7xu3OhlmH9c42Cne4eYiisHPEDMOdSw37utWp0C7m1NnW\nZH10MticrcO7Y98ex3nsOdWCklsBJdfTquQqEM5jc/s2zO37sj/0H2FYg7JvNPMilFZAtt1pakVy\ndT2tlkQN0jw1NrClikqOEJI0vMgRQpKmI81V3EuEGX3dxVlnL5E5GdBEvRq2j4v8WT9hLZlvToiF\nxbi/aHAg7hsyh8NczlmRPRdN3Ico47V7bOJqGOq/p84IL/gf0jhNQP8Wmrako7BFcqxsX/ZqXnop\nkiB6wSIZcB+cZ65+mDNdsx/ENKTvtTndmI1REF/M6IbNOJ0Fprk0hrLJ+Fk1DpoDru3cVvbBBrge\nXYPpO3ojjPWaL8LLOiwS5jyWWhzRpZovk5KJSo4QkjS8yBFCkqajzNXXnIrhxXV6oFYhHMb2wZkX\nNetH0/W6juOYyf+pW2AvTKlX69M4NHMobnSaOZId3MJIlPQm5T0ZLyLSpVK+pyfamyfHMilfQMlu\nfchVV8TiOBr6gvnJtquOJXk8ZqZi6VM3XyEWhBluaiUG20875iruDKjD+o3N45sr8fHF60P2YMTi\n6sFchagv+WIoM2dnavGHdXc8m9v3e3ANKWJze7/cD2O2G+FAf/Tcjg3rhzsmqojEuq5wrqpqruIS\nwHqhkiOEJE1HKTm782Fle3eh1lF3GGzv3+1qoW8K7qocC2Of/J++EaR8CXe+GRgbgn4tU3pTx6PX\nY+W51rtcFyzE2p3NgvtFRAZL2e105Gm4tZlq87IXi8iIpnKqMiVTR4BpgWwmVjC6wfoo7zAjg/5v\nMaOvRSrk97xVWsbQyWAKbvEKTFSb09fh86a0BStFDkH/iLYQsLD4UN8Toov2g0Vic/sJmNs2z3G+\nH6pqhuEG7BXFc2HnCPaF2rkswU9kvdmzqeQIIUnDixwhJGl23FydBJlvKrXsyPhcH8Yso+8s2AFe\niJZnugYTVUTkPW0vw+eZvJ+CsSPQN6kPESufqP3c81yrjBeJeboGYcXXQmcqh6MW7/YWqqFv56gK\nzgg7l+fpgNh2PMsUK8mHjNBF54kiwVyzojMicc8l7sM0E9ab7yLgZMBllyvSOlbXtt1SjE1PSBoR\nPqMrPnH6OcyvmNm2OLctPGwBvBr2HYtl+HA8F3aO4PzZuURzFV/yOKjkCCFJs+NKDm9sQag4d7hc\nHxSN1WTA4GOLVJh2MqaKwDYRvLOZgnvPGWun5JwsMlaCbWowKsfBkdbsqhWJefLt2OcG4krrUGXR\nXhBxCvlWcIuJczhkeyh6fU/JeWMSCz9jYofPnNoMpupwvmMkQ1Bo6GSwee5ZKfeiEpOZGMngzm27\nWsTMY9I4FD/7yWo2zw+BN8P7DvYdHxSjkute57lyz/MaUMkRQpKGFzlCSNLsuLl6EPpBfoIczpU3\ncxYkTfqizLeMvrhgixl/QyRDHd7b5DtK+imT8v8NY8/DE/paj9f2Gh2JOcA+HYmfbcdkx9juOwwN\nqLmKmhzPxROtD+O5JNtLzoyyrZK4ZbK/qW3qLw9k8wXngOVtw0V7z/wLwfYi0ZGASyx1bdGEvWfL\nJe/DGJSIuq5mKM5t+0hwUOBnWzC/VzAHc9CF+T4QfyOlfvDeeedKz2URzmlxnUkqqOQIIUnDixwh\nJGl23Fw9gH0L5/JkvkjwXOKYVbfHHFde3ivMCRecP7hHaKqpFZFopv6rf/BTf9v6GntPCJfBz/58\noC93jHicn+PZsO8YFb0r3w9ACNwBhnhtO7Y3EX9IXfYH/u+8Mehb/VOsg2rzAses8AymLc8F2Vsf\nw7VsTqInNZipvxOXe9/W18Ln2Ht6nwfHdLcnvsaOHX+L3neVXjBX7bzgSdWxLhjrUnMV99qel1ao\n5AghSbPjSu4u9k2JOAVdRCRGFsBYn75DvrL97ZYxzOgb0iXhDu8jTa1Ik5PBgDHvNfaeEOyMnx0D\n9G+3jPXh2bDv2KZsm52ju6De8FyS7cGiS94GDfHQK7nnjUHfqtdjFXubFzhmZQOx6AwaKWHNH4Pt\nbU6iKkMnQwDGevvyr8X39D4PjglLG9qx9zlj+L3c8/KwdewhjFk3H+HTquWo5AghScOLHCEkaXbc\nXF2Gfii9iCElXm1GKIDh5auyyvZY/Rur3FtGX8sHJyJxUTUXzqKS3RwMInnT9JS2WHmnZs+L+hs/\n247JjrHddwjfEepR5s7F7daH8VyS7QX/D0tqhRW9ZZc7zpiIHFzM5svBgfhftLxsXsD7kzBmRWdE\nYkbf3Dy1p+Lctn1w5mAQiSaqSJzTNXiNvSeYsF+Bz/aOzY7dzTG36NRixb6zPLME5xTP+eOgkiOE\nJM2OKzm8Glus+cl2ud/nmp4oIuXFTPJUBuLgIaeyfa5soNZkwIy+uYVPw3Z5tAvQt7vdKRjTrKlH\nRm6GITwO66O6s1RL9l1EJH5H7/uLhPz4cCrWfWcjm8+S18eSe4uPGRORbn0RKjlPGcUED/E/j2UD\nQ00GVG1OuqQwt9EZkct6rS1kAfbUXaUaE02EeQwT1fsOpuS68aSt81y553kNqOQIIUnDixwhJGl2\n3FzFPS4ndI/LHJhlZbTHGk2tiHR/nLXVgSibzSRccAKFReLuayw6Yxl9czvSbV9Q20I22oKkf+q5\nv+hD9TB2VFpNVzRhq5pbzr6LiLjfFft2jmbhYWYE3jm8apKrYG4VPEeSsyxzaDhO/sa+zMQr55Y2\nsv/4Qi4/W4yesbKBoegMgr92e/lahWzQqabzfGA8/iC8eYy5EqOjLZqrhx7p93GWX0TEdTDaufTO\n81pQyRFCkmbHlRxi1/9bsIs/p+Smm1oRES3NN1yJ1/WFUqacMG0Nxs55hZ+tJgNm9A3pktre7TJv\nBToZTMF9FfLaeKoup+7m53PfRUT87wrnws4RKjmyc2BZvLfUImmA1KhaH/9hWFxa/7fF2VimrzKc\nPXnZSVN0u83cDoWfwboINRkwbZIJvTbFpcPjtThkCu5oD1omH0M/m6xVmKg2huoufEf8bXvFpbFC\n53zL0JqlCA0qOUJI0vAiRwhJmo4yV19W+fkWBNnWwFwrmjmHZfo01r4AxW9q/XUREVnpiY4Fz0T1\nKttj0RnL6IupkjDY3tuPZ2Yomqhoun5V/pw9vhIfL9zQDpqrf2kdW4JzYa9+mc6GjuOWtmiZVj1H\nEjqabBsnzOPhgcw8vNsfzdGYfml/GPPmNla2t7KBWHQmZPRtY65aJAPug4tOs3jgNfkw9N2lGHvN\nHfDeeUsxazjdwlKW/PVQyRFCkoYXOUJI0nSUuWpgvY0RCMg9ayoY649aURvImtvf+4WIiIyfgHcC\nRW9mqlfZHvckWSiYZfNtfk3ZMVetj5LdTFQRkfFH2TH1X/siHtA1baPyj/2b8DQ4F3iOSGfxkpNj\nbkTtrTIutaC5ZmYqFGnqVW/o4ePxiSuaTfeRY6KK+HPbchdaXVSRGAqGGYYxR52FYWH4mO2D85Zn\nROISDY4d1uf2wjwOD+P3d3YRzIG9b2bqS19ieYZKjhCSNB2p5H4GV+tfwt2wpFf7Ua+sG9aF0G9V\nlLjn6NjE1dDv0UXZfHomu3PF24eVDcQd5ZjR19IlYbC9p+TQyRAU3AdwvNeaWuh/CHc4fMnP6HDo\neOrQr6iPqwSOpAJWiLeahs7cLvXE7BGPRlESZewDB9p+uS8i+dRGNre9yvZYjwEz+j7pWDb228B9\ncDjPPcdD6UM99jWcaqjkVnW8Dvtl6/LloZIjhCQNL3KEkKTpSHMV+YkTLtMHZl3VzFTvm0Aerf57\ncaH/5Fi2MW2wFBdibYG1AZvwLIQGw2byC7rLudeKwOLsfAxACfvgRHzT9Iq2YI/O6uOX4Gk/oYm6\nq8BF8jd17pYgLmkUzTWvBKUzp5/S5ZJ9o9FL1eMUv8GlGDM9vcr2WBYQC8sMPiaXnedoE4lOhmCi\nivjz3XGqoela13OEzrUv43AwqOQIIUnT8UoOMVWDBamfV/WDu0qCgmtT2rCgd4qRp+NttXI4688N\nxDtTVHLxE7FsoN0tMaNvSJfkLbSKxLsY3tn0O8yBkrPSvz+gekuCurYQ0CB9sGUoWCS4M+Qxv87S\nSlRLfUejqXCgv1XJmWpDJXd7DSvFHBdevQZ0zmEkQ9gmgvPd2x51o6kVkVl43BRcXTYHKjlCSNLw\nIkcISZpdZa6+qqbbBaySrXtpzlyOQ1UzU7GkGRbFMCsVMwyrv2GoEp84NKB93LvkvaeXvdfJeSci\nrrlqTob34Wnfo5maFJb77DWYu93w+AGdF0VxeOj0YR72QtGaseHMfDxUjfvbPt2XLeZ4eenaOR68\nMpmW0Rdz3uXmuWeu2nx3nG9LYKLiHlDrrzdf3FpQyRFCkmZXKTkDr/A/1zvjXXj8tN4pRlF1YeZR\nU1vDMGY7R9CDYbdVdOtj6ULLRe+9dzslp3c7jGQwhwqdDOmDqbFeB1XXpSLqDKiboOpg53+Yf2tY\nKcVGVFvFcqbuHhSjk2B5IJvUbZWcFn7OlQ10SoK68afOfEclZwrufXAMopLb7PRhVHKEkKThRY4Q\nkjS70lzMv7ivAAACY0lEQVRFXnWkrQX1z4OUngCTsmhyG9PeeOZqf1MrkjcTrO+ZqyDpMaOvpUtC\nec5Ihr3JS47p+hBMuBMaCVP19nuuUdowN491c143JAQo9avd2wvrL7gUY5/jOexwvmNG38cUmpp1\nnAz4G9hIRMNaUMkRQpKGFzlCSNLsenPVw8y/X4D36hZI/ppK5xHwApWdQiIhjRzmqkPTwfYngblg\nle1vtcmFZSErzAdHkJecfXTmHx0HU6+m5mMBzUivZmkVxsxNi/nrvIQAa5mrZiJjhR5nj+jqGsH2\nZqZuVxEmKjlCSNIUtu6tJ1e37r03ht0t8WZnffRF2A3wAARN3wWFZjc2vJnNNrUinVg28PwW/t/3\nAtszty2ypwZj49rWYE6WcSLbBEaLZKCpFYlKDm05jKzwlJz10dEBk99qMmBGX1NwP9q230Dr3KaS\nI4QkDS9yhJCk2ZPmqsekmgYYIG3hzJi/DsPHLHQZtyyd7zjT1IPm6sbYubn9umPCjkA/LLuAuVrw\nzFVzprVzPJiDDczVVe03wFzFZRkrG1iHsa3c/+ZDc5UQssegktuTUMltjM6a2xjob6oO/Q7WRyul\nqEquCxwPD8HxsKRKzgusQL/DLehvv2rzoJIjhOwxeJEjhCRNkhEPhOwlPDMRs2e75qqao10QwYPb\n5JaaWpFopm5Wxt7tgkqOEJI0VHKEJMh61dYkKL7dsf3pr4dKjhCSNLzIEUIIIYQQQgghhBBCCCGE\nEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYSQdPl/Lqkd\nyBZ7pYsAAAAASUVORK5CYII=\n", "text": ["<matplotlib.figure.Figure at 0x10a6047d0>"]}], "input": ["D = solutions.exo2(x0,W)"], "collapsed": false, "language": "python"}, {"prompt_number": 24, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["## Insert your code here."], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Display the geodesic distance map using a cosine modulation to make the\n", "level set appears more clearly."]}, {"prompt_number": 25, "metadata": {}, "cell_type": "code", "outputs": [{"metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD9CAYAAACY9xrCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHTtJREFUeJztnW9wluWZxc8jrKAEQSsIFWugyBIICuIgHaOGVgV2CW0X\nlD++DkFBoKJEtKIk4qsSFSsYEMufgIYh8kehVcIWUFqjxhEYEZBAWKAktlgQrIJEBRe894M7sx+u\nc+2ss/t+2fv8Pp5ccz/3+7zPNZk57/WcGxBCCCGEEEIIIYQQQgghhBBCCCGEEEJ8b5JMLfwIEJie\nzrVa+c4UXePOgUutuH6Wc8UeXK640UjvjepJS/v22kH1GdutdrGzi9TNVvvg5RxaOwDrqX70xh9Z\nceOHzhU3OvoXRLvAqb2MaJc6tR85+j6ifebUnufoNxDpclrZ5o2/GG09BtDaK2+po3rlK1b72NnZ\nFP7IYNO2K4z2kyXkgQGAwjeIuJPXDphM5YXrbjPa2B6VtDaptf19Fr+aEOL/I2p4ISJCDS9ERKjh\nhYiIjJl2IZebdk/vnGi0KT2e44vUpok4gZa2Okkvh2PvtzdaTR6/3FtcRnFHqx09kEVrLztlzavj\n2e34wocrnCs2EC2fl2Y7+oNW6jCOGWtAf2wwWg/HTNrpmKMb0N9oBxcwMxDAU1xGQzURmQYA2VZq\nV0grWzUcpvq+ZnZ/bTo10trSer6L64mWV8NrW191yGjHm3stOI/LuWkjzdh5Ny2dksyVaSdEzKjh\nhYgINbwQEaGGFyIiMmbaLQwp6qLd2YNMz1FzDgAeMUpe4JNl7zxxE9XnFVutpXO11AyuFz7wW6Mt\nGcjNQ6xnU0/edBo3WzDXTqKl7iqnpbPAJ7LaPEbMp+X8crv3WM3bsTd/160rEUfw2qPTuOE5GXaK\nsvL5sXyRiWyS0DF/vV0PsBOeo9Zxs6zi6V9RvXKK1U44u5hQarVrp75Oa2sSMnUIAHjUSsTIAwDU\nJjLthIgZNbwQEaGGFyIi1PBCRIQaXoiIyJhLjwGBz7quTxPRuvEAMCisMlrV2Fto7YxFfBvM6+xt\nX6UGACTffMv/0PkIEZ3RR9xqpQo+Ylo2ahzVJy1YaLQvJvGrVZziOnsT3XuHn72tn9OE19adcXSi\nee+We2/lFzaz2nmzee3scXcarWjJAmdhPlIMvEQ059eX/W2pHM62/zO3kjgDgCcXTBnDawvKX6b6\n2mQoUYlz/50ul16ImFHDCxERanghIkINL0REZM60w0xu2hFTKy/wgMZ3xtpx2VLHnHO8D3we7Fhl\nzswGXnz/SmeVT6yUfw+tHPWmNfMqKvlY5oc2jxAA8DuidXZ2liKhoACAqVbaMYKbh6tgjaAtuJrW\n9sFmqg+FNVivWO6YZU9wubLWavt5Kf6FaJeTqW0AKEzZ0WgAWNKPGHTVc5wrXsTlZ4YZqe6+bFp6\nfmIHlp3HGcXOA31tuR3FrUl40CfQTqadEDGjhhciItTwQkSEGl6IiFDDCxERGXTpX6cufauTNvX0\n2EybLAsAM0h4RaFztfpgj/wBgJ9MIMf+zPecWCdldcVAI9UO+zEt7d7pgNG8xFOWGQEAQwqt9sGL\n/Liq4U6qxb57yf0ocy5InXdvHNW5R8zVL3JWeJYf6bWCJGZcOZofE7W6wmokxwMATx0GgF0HOhkt\nd+WfefHwdc7q5D6N57/gvDfPnlfVMeH3osK52hQSotH6PpuGCwDHm7eXSy9EzKjhhYgINbwQEaGG\nFyIi1PBCRETmXPoKHoAROttLznPOeutDtBZkNh4AciY08EXms+xp67oDwFmHufN+Zr+NVa529vwu\n0YqtOQsA2LXNusQAkPvxLiv2bM4X+bSC6/R8Oi8C4xqidXNqdzs6++ReBEY2ly8stNr2k7S09uLu\nRuvey/5CAgCl5IcagH/qfOdcuCad+Zlz37Zjrr7j6I+3mdZ187Jp6Zdk7h4AthBtgrPnJM/2t/7D\nCxERanghIkINL0REqOGFiIiMmXbvhSuoaXeajBI2OGukSLps8rKTq3G/Ny6bb5Q2oTWtPLKIG4Ll\n5Hizc5yrpeZbrWCck0Da62a+yHaWiMvOUgPc8+lKzjVSzuMf0NIxJIYhH9W0tprcTwBYRCJI6h6+\nku9t+ldcp2fD2XP2AAA9bXjFoG2v0NKqBTzpuHK81b52djaWH+2HtmOsuXY0OeasUm2lZ/gYbriF\nt2YlScTNdq52Lelv/YcXIiLU8EJEhBpeiIhQwwsREWp4ISIiYy596Alqp7Mxx2I2/QogGULOeuvM\nHW/PzT3r8HVGO1NlR2UBYA5x4wGAhQDn7+W1yTfkY+d646jeZ5lspfX88y3uP5Lqt1eSYIxp/Gpv\nkIAOLx7ai8u+kYVMPMZrX0jZoAsAuGPDMisO8H6dmEU07sajlo8Jh7Pt41/dhS/Bg9SBe4h736TA\nG8N9m6jO59vPP0tYbf9Hl9qJXQBAiVx6IeJGDS9ERKjhhYgINbwQEZEx0+4pcNNuCjF3Cg84Z38l\n7AQx5/y3Fc6IYgf7Ecudd9m9PNa+xIM5Z4Uz4juG7e8zXjucnG0GoGS5PRju8SeepLU1JNkXADYS\nzXsb/gaidXTe4a933i1n1/PehmfXA4A8ksj68NSHaO30EeSAuhVsJBkALuDyInsu3NfDeUts4j4v\nzfYd672ffpA8M8O9kXC7NwAYFezJgxWd+NmFSb1MOyGiRg0vRESo4YWICDW8EBGhhhciIjLm0i91\nXPr+wdqdbQee4IusT1stn2gAat/kibNHE5tkepBfDSlvXPZt8lHGOIkIuMhKqwbTykNDeBDHeVnH\njfb0l/xq+c4u8idZbU8ZD/h4jMzcVn1ZQGsLWlRRfRqZo+1axJNXq2dTmUZuPNCC137R2Mpo7Vc7\nwRND13Adn1hpEZ+vDtc5gRRkFLeDc7U2gZxl1885y646zfUBVj+yriUtbZs0yqUXImbU8EJEhBpe\niIhQwwsREWp4ISIicwEYN3OXvvXSQ0Y73pxNYgPA1UYZFXitN09cSsIdikmUNAAkec58fC6bj3eC\nqt+3jnw4zW/zur58iTqiTeamOdau+SnVC979oxUH8TVwjJ2FxnYBADlcbk3O61vLS6uu+RnVBw3+\nk9Fm8R8F6C4GbuK1SVPne72KufdOUHUtn20PNfa7LSXx1wBQ/L3eI/HeONhslFYnee3x5u3l0gsR\nM2p4ISJCDS9ERKjhhYiIjJl2W0MOdUp6tycJrodJ8gEAVNh0h9CEb7n0Nr5EMQlyKNjmnPWWdOeL\n4C0rreLhFSxwo9Ix5/7Budowm3GAbr/cSmvrWjrntzWyKGDvirdaKZeMCANALRlHBQC8RLR/56VZ\nPGY154Q9+27373vT2pUkG8W5GlKemccCKYZ6IRrXU3VQ2GW0ql48cZYmNi919nbGMRoLSa+0c1JQ\nDicy7YSIGTW8EBGhhhciItTwQkSEGl6IiMiYS98mfERtxqOJHZ8EWBw1UBZ+bbR+yUJay+KCAaAr\nCx3o5YQObE9zfbjVDy3n4RXbEhte4Z2ONmwb15NT5Nb1dbKP2S8IANCcOLd7uI+9+NJRRrv9r+Rs\nOgAvXOKcC/fREit2dX4VOOn8KsOc8E08Uzw0s4/uyl58VX4qH9ArkBCNEU6Ixoo013tavXYbD2PZ\nQ8JYvGj0N8OdVC9KfkNU8rMOAGC0XHohYkYNL0REqOGFiAg1vBARkTHTDjcEPhu4MW21uUQDEJra\n7aWdd43ThVxPppP3mzu8yIvZiCmAkvCU0R7K4me9zSfpspMdTyX5oTM+2Xc1EZ33tKenqPxqcX+j\n/Xzk67R2HfHn7JDrdziDvBhIvLzXlt1Ea39RuoEvUlJJRCd3YNMQI4W/8cd5FveEMZ4k4j7Z6Jxl\nlzzIF2EjxQdH08pQYj9LuoKvmvYyG06TZ2ZimhfjUZl2QsSMGl6IiFDDCxERanghIkINL0REZM6l\nxw7Hgv7MKKnAB2Ofb27HC9ec4qt2CzxNtXcbErjxKQuHALCeBzOEreTXAidzIE3SZbutccIrkq/4\nImiw0lruxodu/CvcaieK4QTAwg7WAh3H8Nr6RVwng7VwgnbR206YAgCS3eSRGcScewDINkpOOJdW\n7h7MQzTS5IaknanfpLfzOA8gz9KF/DnaerSb3VvC04EHN+OXu+ukHS2vTLwB3X5y6YWIGTW8EBGh\nhhciItTwQkSEGl6IiMigSz+T25rZk410pL4lLV2eNBrtnlx+tS47SQYwgH0JS5ngMcKLA7emuyR2\n2Pw03wYagz3rraAlOecNABodS3g6ieceyb+q1cSNB/jkfaqM19476Qmjla3mM+VFQ/g7BM/Onmq0\nyiJ+PWc6HkOIe58sc9zxEnLvsvhPJ1Un+Fl2WSSMpamzt73BCf5I2M8WPAb9smATOvb2IDnqAObU\n8n2MCFlGa9vxBC9uUEy1EFGjhhciItTwQkSEGl6IiMigafcId1vmp40UPuHbSD9CtGX8askWx9wp\ns9dDCdEAhH909kHOrUtPcvZxM9lHnjPK25yPYL76tQ2v6JDw8Ao+mAmkiFeZfOnco7wPiegloTpp\nEjWXGym04Pez0kmXZcPRB4MTonEOCdE46dznGmdk+hUyMj2bL5H2zoD7N3JPp6d5cZHVQx/nmRvp\n7ONRsoeLnO91vEw7IaJGDS9ERKjhhYgINbwQEaGGFyIivEnC/wMuoGqHcSTsgmdX4GKi7RjhvOzv\nuJpslZzHnRBmZ0yV7WNP2aW8+HwmOmesOWe9sTjpNF8BaWdcljryeRXOKmS0ebxzxfksQhtAnr2n\nSQ13j0OZ40yTUdy0E61N7122c58HOUt8br/Di2d/xIuncTnngP3cddPZEwOAfFc7An+eLx7pnJZI\nIsU71PHagyTSXf/hhYgINbwQEaGGFyIi1PBCREQGR2v/QB2bO8KfjTY5uZuuYPNtgQ3BO/vr584+\nrCE1M8yllT2SeVTvTLTisJjWLk/aE5WfyLY43Ev19uT9+660EpgT7LvsAFCW/DNRublzwem+Rvv7\n0g609ge3HaT6Z003EZUbUkXhX6l+T2Lfqd9DK4FD5P30O5JnnWpu0o4Ih4xWmtxBa/c7K+8ME4x2\nXzLRqbbvrZeE12hl/4TnDjArfFZ4jtYuTu7WaK0QMaOGFyIi1PBCRIQaXoiIUMMLEREZdOl3UZe+\nLNiEgS6JPS8LAPo0sdrI06/S2tcTJ7kT9ky2rcGe8QUAm51zviaQYNGWNUdobWMWSUrItUm9ABD+\nwG9/6Y+sVuyc9ZYMcMIPhqat5ozLhqtJEMRovmz6RWcfm8k+SNgJAGCVs4/1dh+lzll2xX8he/gn\n517UzqJyVqNNNjmR15bWzuOhyLianGnYOyHnGQIA7Dl5NwWe2Lys6S+ovuWM1fYGewYjABQlC+XS\nCxEzanghIkINL0REqOGFiAg1vBARkcEADB4ksBM9jHajs0IdcST7YDOtfR3eLL11TKuRTyt7OKHP\n9cShLWhRRWuXszSP2k9o7QuX8PPKriQpB/WOW11Uzmeuy1ictBNe8YO5dj7+7/h+s/QYzdbmkdbe\n+XT1Q63G30Jw7p1zn72EFfYdsu8a4O9TAN6z5Ln09t2CPuCz9OzZBwAWu7Ka9JSH/sMLERFqeCEi\nQg0vRESo4YWIiAyO1s6mc44dwkCj/TWnC12hnKQf9HFSPnsme519lBslJ/Smlbs7cb283mrXBp5a\nm3N+gxWP8dFONPDgj/DQ2UZLk7RS4L9Jrb3qf5taO4SXeqm1JNwBNYW0Mrz/PVJrua+J5MlvrJjN\ngyDQmo82132ebbR3Em42j+3Il+52YKtdN7Haf65ilO2BP/tbEh5WMpYkoVxSx5/9g0kXjdYKETNq\neCEiQg0vRESo4YWICDW8EBGRwdFaFjINHFxAXHbHif34Eatdsdw5c4s4vACAso+NVPewdUsBAI85\n+7A5Ceha5JxBtpZoefwMOXTlZ6G99vVNRitYzs9Yq3Q+d9hmnXDvrDfkfWg1L7zCGZdFzeV2Dy24\nG+/tuYBory2z9wIAcA67d859Zt8J+He4gpe6z0bdw2z4dw0vJp/be55/7+2D9ArtKQf9hxciItTw\nQkSEGl6IiFDDCxERGRytnckdomw75niknid3Lk8ajXZPLr9al538ReZ9yTai3kJrFwceDduFnPV2\nmm8DjeGnRito+UenuJTr04uNFEbyr2p1J77E10RLOWO4906y59OVreZn+Hnvsj87254L55lz53AZ\nQw5YLVnmGI0l5N5l2fsGAFUnfkb1rORPRvNc7L3kLDsAuCNhQQUv09rLQi+7bg8SiQxgTi3fx4iQ\nZbS2HZ3E5oZEo7VCxIwaXoiIUMMLERFqeCEiQg0vRERk0KXf4dirduQ2Ffh44fPN7ZlZa07xVbuR\nM74AoHcbkiD66Qy+yPopVA5bydlr3BBGmsyHdlvDAxHqkq/4Imiw0lp7Rh4AhG78K9xK3HueswuM\nIlpH5yw7Lz13CdHYqCwA9CZuPAAku8kjM8iex/Yd2UbJCefSyt2DebBJmtyQtPPDSdLbeZwHkGfp\nQv4cbT1qzzTc7ZxnOLgZv9xdJ+05jJWJN1rbTy69EDGjhhciItTwQkSEGl6IiFDDCxERmXPpbwjc\n1tyYttpcogEITYk7Pp5fLl3I9WQ6mSrv8CIvxq1ULQlPGe2hLD5TPv9Lq03+nbO3HzrOb18WBc2m\n4wFM5+79q8X9jfbzkTxEYx2JwP6AX809620gGTX3wit+UbqBL1LCHHln8n6TjdEOf+OP8ywns2N8\nC6s92cjfIZiePMgXwUtWOjiaVoYS+1nSFXzV9HyuJ6fJMzMxzYvxqFx6IWJGDS9ERKjhhYgINbwQ\nEZEx065N+Ig6UkdJ6ICXhFoWfm20fokdLQQAJ8sWXYOdMc3t9WdevD3N9eFWP7S8NS3dlhw32hfO\n3oaxbA4AySly6/rWOKu8xeXmZPZ3D091XXypHa69/a/8MLsXLnGCID4iw7VOKi9OOvOruN5Km/Jo\nZWhmH92VNl8CAHCec7VeoZXR2o84xotXpLne0+q1235MS/ckdqbYG4p9M9ixcgAoSn5DVMcVxmiZ\ndkLEjBpeiIhQwwsREWp4ISJCDS9ERGTMpd8acqhL37s9CaQ47Li2FSSuuQnfcik5/w0AikkKcME2\nHiO8NunOF2FO+KoJtDJ0sPur7MtXdTxsDCOma7dfOiEaLZ1h10YW8uFdkYwU517ES2s/cdYgI6be\nWW9ZPCAi54Qd6N39ex5esZL8sONcDalNXE8Okkd06DxnFfILAoBBYZfRqnrxGPRSkqRevNTZ2xln\n7LqQ9Eo7J43lsGKqhYgaNbwQEaGGFyIi1PBCRETGTLtwM6jr0HrpIaMdb77RWeVqo4wKvLai06+o\nXlpvtWLvXeM8xyjJXUlE5z3t9wcbKZzmt3mdY+axHNPJTgTs2jX2LDsAKHiXnGc3iK+BY+v+h7sA\nAJ4OjNYDrbaWl1Zdw896GzTYjl3PcqJ22S4GeuZcU+d7vWoNEZ3cgdphVA419rstdTIbijtarfDA\nb2ntkuQGvgg2G6XVSV57vHl7mXZCxIwaXoiIUMMLERFqeCEiQg0vRERkzKVfCu7S9w9ZRms78ARf\nZH3aavlEA1D7Jg8dOEpCBw7yqyG1l+vJ2+SjjCl3ViEjqauscw8Ah4bwEI3zsmyIxtMkDRcA8p1d\n5E+y2p6yS2ntY5hmtKov+c8CBS24bT4Njxmta9FHtLZ6NpVRTbQHSLIsAHzRSMIrVjvhFUOZGw8A\nZEx40VhaGa7jrVLZxWodnKu1YWEs/Zwwluo01wdY/ci6lrS0bdIol16ImFHDCxERanghIkINL0RE\nqOGFiIiMufRPOS79lO81T8ziq9lcO4AV91CZBVKU8+RjNzK4b6PVzlnhzGePYfv7jNcO5yEaJcun\nGu3xJ/hZdjVO9gF74+BiXgo2id2RBIcAQD0JcfCu9/H3uB4A5JFsh4enOme9jXjCiiu88IoLuLzI\nzsd/PZy3xCb74xIAHo8+1kkUp4Ebw+fwYvDZ/VHBpqN475Ek9ba/9R9eiIhQwwsREWp4ISJCDS9E\nRGQuAKMnN+1ocicLWAWQDPnWip154qx3gthZh68z2pkq7sDM4VOVuJxo+d4Y7jfkY+eSpF4AgPdZ\nJltpPf98i/uPpPrtleRsODtBCwB4g4SE7Hd21tnRbyRmLJm2BQC8kHLOp9uwzIoDvJP5ZhGNp8Wi\nthuVw9n28a8mo7IA8KGzi3vIhHWTAuLyAvi23dtEdT7ffv5Zwmr7P7qUhwCjhPS3/sMLERFqeCEi\nQg0vRESo4YWICDW8EBGRMZf+vXAFdelPJzuM1uCskfqL1ZKXnZHW+70RxXyjtAk8eOLIIh4QUU7c\neyekGikSgV0wzjnLrtfNfJHtbETUc6vv5nLJuUbKedye3QYAY7DIaPk0jgKodiI3FmGM0eoeds69\nm/4V1/Ec0fivE+hpx5IHbXuFllYt4I53JYmTdkKqMdbJO2k7xoZ8HE2cIA52T59xRsJvcQI3fmS1\nbOdq18qlFyJu1PBCRIQaXoiIUMMLEREZM+1QEai7FjrbS85z3k/vQ7QWgRtrORMa+CLz2dwuOQcN\nwFmHefLtmf12FLfa2fO7RCt23i3ftc2mmAJA7se7rNizOV/k0wquUyvUeyP+GqLxcVTAGxNmn9x7\nIz6byxcWWm37SVpae3F3o3XvZROKAT7ODfBPne+8y96kszcuy1Jn2Vl9AMbbGdi6edm09MuEJ/5u\nIdoE7/37PJl2QkSNGl6IiFDDCxERanghIkINL0RENM3YyoVvULn1yUNGO1bantbOIImshY57+V7g\nVvhPQCza+XwM99t23FVOSEJtbeCOfnEn6xR7LnFXcu4dAIRCO7j7wdEcWjscJOgCwL57C61YxvcB\nbCZapVPrZfuS+eMiZ4Vn7Xg1AKwgvwxcObqO1q6usNqrzs6KWTgHgF0H7K8kyUrnrLc8x3lnubXE\njQeA9+bZZ/R853l+zbnaFJLs2/oq21PfYftK/+GFiAg1vBARoYYXIiLU8EJEhBpeiIjI3Cw9ZjpJ\nFbcaJS/wEOB3xt5ktFKb1QAAJH7hOz4ns/c5Mxt48f3OuXX4xEr5PLhg1Js2vKKikp/99eFt/Gr2\n9DA/HjqV6/zBHk+HHSO4w74KQ422BVfT2j7U0QeGYpXRrljOTl4DQI6FA4DKWqt5cdns1MHLl/La\nwpRzdmE/crZftRekchGXn7FnwNXdl01LmSPvPM4odh7oa8tfN1pNwoLUAaCdZumFiBk1vBARoYYX\nIiLU8EJEROZMuwE8AAPr00R8hJYOCtYIqhrLE0hnOO7HDUTrTdJwASD5hpxlBwCdjxCRJcsCzJRE\nBTfLykaNo/qkBQuN9sUkfrWKU1z/jGhe/AUb2s1pwmvrzjg60bz4iwscvbCZ1c6bzWtnj7vTaEVL\nFjgLO+YhXiIaMfIAYH9bKoez7f/MrSRZFgA2Em2KY84VlDtJx4k1WIFH+SJ4VKadEDGjhhciItTw\nQkSEGl6IiFDDCxERGXPpF4YUdenv7EHmH2vTzirWvc8LzOsE3nnCjuECwDwSotHSuVqKJVoDKHzA\njmYuGei4uetZcAQPOXDPhZtrz1NL3cUPN5uFyVRv8xiJVeZZGdi9x2rejnlIONCtKxFH8Nqj02zs\nNwBMxiyjVT5PgjUAYCI7a4+dTQe4ux6QMtKodfzXl4qn+Xh0Jcm6OOHsYgIJr7h2qh2VBYCahP2+\nBFBHPjfNS2sTufRCxIwaXoiIUMMLERFqeCEiImOmXcgFNe2e3jnRaFN6OGYLNfO4WdbqJJ/kPfa+\nTe6scc6Fe4vLNPX06AFuPF12yo5xHs9uxxc+XOFcsYFo+bw029EftFKHcXzEtD82GK0HdtLanehB\n9Q3ob7SDC5yE26e4jIZqIjINoOfTtSukla0aDlN9XzO7vzad+BlypfV8F9cTLc85642lyx5v7rWg\nM7pNDLoZO7n5OyWZK9NOiJhRwwsREWp4ISJCDS9ERKjhhRBCCCGEEEIIIYQQQgghhBBCCCGEEEII\nIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEL8F/8BdtP0lqErnowAAAAASUVORK5CYII=\n", "text": ["<matplotlib.figure.Figure at 0x104b28490>"]}], "input": ["imageplot(displ(D))\n", "set_cmap('jet')"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Computation of Geodesic Paths\n", "-----------------------------\n", "We use a more complicated, non-constant metric, with a bump in the\n", "middle."]}, {"prompt_number": 26, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["n = 100\n", "x = linspace(-1, 1, n)\n", "[Y, X] = meshgrid(x, x)\n", "sigma = .2\n", "W = 1 + 8 * exp(-(X**2 + Y**2)/ (2*sigma**2))"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Display it."]}, {"prompt_number": 27, "metadata": {}, "cell_type": "code", "outputs": [{"metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD/CAYAAAA6/dD3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADthJREFUeJzt3dli2zgMhWF2m67v/6Td17kpG+QEgEDJchvj/27k2PLa\n4RyQoqgxAAAAAAAAAAAAAAAAAAAAACx6cvLr/zr59QHE3Pb99NqfAsDfRaMHmqHRA83Q6IFmaPRA\nMzR6oBkaPdAMjR5ohkYPNEOjB5qh0QPN0OiBZmj0QDM0eqAZGj3QDI0eaIZGDzRDoweaodEDzdDo\ngWZo9EAzNHqgGRo90AyNHmiGRg80Q6MHmqHRA83Q6IFmaPRAMzR6oBkaPdAMjR5ohkYPNEOjB5qh\n0QPN0OiBZp7/7Q+A63jy5MlV3ufXr19XeR/sR9IDzZz9v3/+t38F10rxS6EauBr3PwySHmiGRg80\nQ3n/SFyqhD+rK3Cpkp3S/6Io7wGQ9P+kvWl8VjXwt1Oc9N+NpAdA0v8TVhK6su9ZlcKZSb3y2iR/\nGUkPgKT/q46kdvbcM6sBdTTFo8cuXR00RdIDoNED7VDeX8nRklsfWy3vL1Hy7y25V8r7PV2B1X0a\nobwHQNKfbk/CVpL6zH0qKgl9zX22ntMUSQ+ApD/NkUNtWUJXHjsz4dVKUus226f6WPa+W89pgKQH\nwBp5F3V0RN1L6pV9ssTf099Xq/3tKOG9ffQzrPbbo+fZ79Q89f8g6YFmSPoL2Ht8/UiKP336dHOf\nyutUv8cY6/31StL//PkzfGzFfF5WKaxUEbeMpAeaodEDzVDeH7CnrN87AKdlvd1H7/P22TPYp/YO\n0mkJP/+275mV+ZVDfrovZX6MpAeaIel3WJl4s2dwzRukm/fpdu8+KwN6lQE8TWp7n27t55r3RYlv\n71uhie99j66JT9IDzZD0Cy6R8F7SR31ye1u3z54929zHe51Kvz+y2l+Pkt7u8+PHj3vvrYnvvf9K\n8q9M4OmS+CQ90AxJf0Bl4k2W9CsJPZPdS3p9bG6fP3/+4HV0e6mk91L8+/fv7mMz3e176mPe59KE\nX+3rV0b2OyDpgWZo9EAzlPcbVla+WT0cVxmkm7ejrX3eLOfnttIFsJ9nRVTW29J9fo5Z5s+t9xto\nWW9fZ+XzVERlfZcz8kh6oBmSPrB6bny0TzRxxrvPS3FNbU1z775K0tvn6z4RL311sM7uM29r0s+t\nvZ1NFlqhiZ8ldtepuiQ90AxJv2DrJJq9fXpN8Szp//vvv3v3jzHGixcv7m2zakCriWyq7pRNjZ1p\nrqlub8/tt2/fwt9AP8PWfcqbdrv13K6H8Eh6oBmSXuzpr1em2Hp9+mhEfib2GHcJrWk+E997TLf2\n9nw9r28f9au9U1q1L69pbm/PrTdm8fXr13vvOff1/h1WVsPVE3kye/v9jxVJDzRDoweaobzfcKnJ\nOVl5Hx1yG+Nhqf7y5ct7f9v75naW/l55H5X59vNE5b09HBeV9V55P0v4L1++PHj9I+V3trpOdgiQ\nyTkAWiHpfzsygOe9RpTwlSm2XkK/evVqjHGX4jPV7e1sH01/r6pYOWQXHY6bqT7GXbJrVZGd5++J\n1s+rrKPnpXf0Xl0G9Eh6oBmSPnA0+Svnymuf3uuLaz9dU93ejrb2eZr09r1W+vQz2TXpZ7rb99DX\nrUyYye7zVtCZv2tWDVTW/tN9byHZFUkPNEPSb8hG5PX+bEqrNzElGr33Jt5owr958+bPPq9fv763\n9ZJe+/vedF5vau4Y/sk0M+FnH35uvYlFmvjZ6bwrq+rapI9W57XvlfX31S0m/ETSA83Q6IFmWpf3\nR8/drgzkRWW+vS+aXz/Gw4k3s0yfpfwYd6X+3Gq5b5+nZb53Jl404OWdQTfL+s+fPz/47NFAnhUt\noW3L69mt0AlB9rfUufY6sGcfOzJIdwsTeEh6oJnWSe+55ll2OtCVJb0mtZf08753797du98+pknv\nJXM0kGeTXifjeK8XrcTjpbgmvbcCz9zOz+kdsptb7yy76N+r26E7kh5ohqQPZOkd7bM6OSc60cam\npR5i8w7HacK/ffv23naMu9TPpuquJP2chKOH6my6R1WTd+hPV96xJ+7M+/R3sp8n+p29y2LvWV3n\nlpD0QDMk/Q5biW9vZ5ea0gk7XtJHk3O8pJ9p7iX9vB317e37V5JeR+2jkf8xHk7jrayjZ3+DaOUd\n+1vO184uxR39u91immdIeqAZGj3QTMvyfu+knCNLYFcWxszK+8pZdlGZb2/rgN7RgbzoOWM8LOsr\nq+t4hy0ry3dHZf3KuROrZf5j7R6Q9EAzLZPes5L+0b7ZQF52sYvoohdjxMta2wG46Dz6ylRdm/Te\n+4/x8LCa3SeqCuzzorXysu+Xnf3nVU3R7539u64k9WNNdQ9JDzRD0m+oJMXePr0exvP6qtEEHpv0\nOoFHJ/LY23oyjk36mbJRetu++Eq/P/p82ffLVhiqrLW30qfXx8e4jUSPkPRAMyT9Dpc6KUfT0pvA\nE51+650SqxN5bKLqKbq6tc/XE2V04oulCW9Xw52vrRN4vDGL7HJb0e9UOZkmc/TU6seKpAeaodED\nzVDeB7zSb2sAaOWMPCtbIlof8xbYjAbDsuvTZ1e21fI+K+tnOZ+9Z7RAZvb9vN9SP9fq77317+et\nsnOLA3okPdAMSS8uPbiz8nrZAFU2yaeySo8OkGWX2YqWqPZWqqkcRlv57JUBz8xZ/363lPgkPdAM\njR5ohkYPNEOfXlQuebTn9Sq8Nd+j9eC9Szpll33SSTTeqrNR31lXo81ez3vPlc/uXXjS7r/l0n3v\nW+rLTyQ90AxJH/ASP/q/vpdO0eWVvdfQtd+zdeG91NV15nRrb+uCFtkxeP0M3uIXuu5dtv6dd4pu\n9P2839I+L9pH/872ye6/xYSfSHqgGRo90Azl/Q6V0k/Ly8oAlXdOut43y2qvjJ6P6YUo7H3Zmnbe\nNd3t/ba8n6+nW+899TN7l8fSAUHvN9C/vdI960apWy7hMyQ90AxJv8E7CUMfywaNosNp9rYOcHnr\nzOnWJuq8rdt5HvsYtTXt9OQZ/Qx2IG2m+KdPn+69l33P6HPZzx59P/u5dLAw+y1X/k1Ul+Qn6YFm\nSPrfViblVJKiMqlGE8u77JOuFa8ry45xl656uehs9Rmv77xn3fv53h8/fhxj3CW/fUy39rNH38++\nl1ZAXtJXJvmolWS/pSqApAeaaZn0e6fa6vOO9ul1oo13meZozXhv8QtdBMP7ftmadkcuYDkT/sOH\nD3/20fT3+v066u9dBUfHFLyxjyN9+r0p/ljTn6QHmqHRA820LO+PWpnnrWW9N3BWKe91ck22/l12\njfj5HvP1jl7AUsv8WdKPcVfqR2W+va0TeSrlvTeBJxvIu1RZ/9iR9EAzJH0gS4HKQF50eKky8cam\nnE6Y8S4yqavEet9Bz5Sbh/dWJvB4hxKjAb0x7hL+/fv39x7zkl4n7nhJn03giX7nlYG8yhl5t4Ck\nB5oh6UXlcF7lZJrKajZRwtvknn1dXS3W7hN9VjttVpN5Jr03NhBNN/ZO8tGktyk+k31uKxN4tG9v\nb68kfeWQXdeTckh6oJnWSZ+dTLPy/KxPPxM6W5NuJpe3Fn1l3Xv9PNl0Xr18tDc2EPXpvdV6spN8\nNMU1+e1jeopulvTe54mmO6+ccFNxC8lP0gPN0OiBZlqX9xUrh+7swNJ8TEtR73JNWt7b+fBZOT9F\nZ+l55b2egWfn8EeTe7yFKbXk1u7DGPFZdpXJOdk5997ioDofP1tuW/++RNn/mJD0QDMkfWDl0J33\n90wYvSCjd5ZddJnlrfv0fbOBvJmcmvR2IC+qKrxBsehwo7dGng7y2WpA9/EG8uZj2QpDRybneG45\n9Ul6oBmS/rc9yT55lzPWPuXcx/ZDtQqYW5ty0WfIUs6bzqvn3HtTblf69NF7VdYCsEmvYwHe62gf\nPrvMVmXlnO6TdEh6oBmSfkNlAk8lTTTxx3h4maat94neU5NYR+rtbU14r0+/tervGPE69V7SR9sx\n4hVzsstjHV05J3JLaZ4h6YFmaPRAM5T34siAnrePXpIpe91sxZts0kllgU09H9+bZ79yyE4nBHmH\nCaPDet658rPM914nGsDjLLt9SHqgGZJ+QbQEdrav/q3JX3mud5+XcjMBNdXtbW8AT/eJeAOP2Zlv\n0bTgyiBdlvTZxS6uuQT2Y0XSA82Q9IGsb7+n319J+Mrn8SqGmYo6tTZLet2uiqa9VpI+G4+onEyz\ncgHLzEqldktIeqAZkn7D3sk5kdXEjxLevs5McZ0o413AUrf2O21NCtq7BmClGtDHKifT7O3TV77f\nLSPpgWZo9EAzlPcHrBzC82ipXznMNMvybBHOufUW4YwW2tTbnqy8r5z9l125d2WQrnKufOUsu65I\neqAZkn5BdKjOu3/P4F52eDBLME1v73BclPCVVXpWziL0KpBsNZuVfVYm3qwkfLfkJ+mBZkj6HbYS\n3z62J0W8VXWzyTl6QY2svx5tve+jVvvQlX5/tM/e91qZYtst4SeSHmhm/VpOa1r8r3TldNksWbP0\n1a3XX6+k+EpfXq2uN1dJ6D39de/zkPAu9x+UpAeaodEDzTCQdwErZ+TtLS2z19kzSHfp8t77+1ID\ncCuDdJT120h6oBkG8k4SpWZl0C9L4cpjlRSvLLNdsZK+lYkzq49l77v1nAYYyANA0p+ukqh7+ttn\n7lOxp3995j5bz2mKpAdA0l/NnsTPHltdP//o+4+xP2FX+uJHR99J+HtIegAk/V91JH1Xk37PPhVH\n0/fIqDupvomkB0CjB9qhvP8HrJTaZ5bulxjI2/u8ldemrC+jvAdA0v+Tzkrqva9zqWQ9s1KAi6QH\nQNI/Gmel+KX87WoALpIeAEl/E85K77OQ5ldD0gOg0QPtUN43ca0uAKX7P4XyHgBJD9wykh4AjR5o\nh0YPNEOjB5qh0QPN0OiBZmj0QDM0eqAZGj3QDI0eaIZGDzRDoweaodEDzdDogWZo9EAzNHqgGRo9\n0AyNHmiGRg80Q6MHmqHRA83Q6IFmaPRAMzR6oBkaPdAMjR5ohkYPNEOjB5qh0QPN0OiBZmj0QDM0\neqAZGj3QDI0eaIZGDzRDoweaodEDAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAM73P1VN\nODaxSAeKAAAAAElFTkSuQmCC\n", "text": ["<matplotlib.figure.Figure at 0x10a8e5450>"]}], "input": ["imageplot(W)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Starting points."]}, {"prompt_number": 28, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["x0 = [round(.1*n), round(.1*n)]"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["__Exercise 3:__ \n", "\n", "Compute the distance map to these starting point using the FM algorithm. _Important:_ use symetric boundary conditions."]}, {"prompt_number": 29, "metadata": {}, "cell_type": "code", "outputs": [{"metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD/CAYAAAA6/dD3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvWucVcWVNv60QEBAYLqB7kATGqQVEQWFCBEydAxRY4zx\nOiZKIv7DTC46iYlmNJlcyP2miZMxiVHzE0eTaLzfNRJtJqBiUFERMYAcIhgQ6AEEBLn0+2E9T9U6\n1bs7OO97vvxPrS+1zz51ate+nP2sWutZawFZsmTJkiVLlixZsmTJkiVLlixZsmTJkiVLlixZsmTJ\nkiVLlixZsmTJkiVLlixZsmR5m1JTycHbL0c7AGy/8AAAwPU9zw/f3YiPAwCemjfNdtzMLx5yA5R2\nc+Ml7WC71XXqx7aJ7WH82CN2OZHtR605Zto8AMDHcWPocv6u6wEAfa7aZzvmWDN/SRzmqeToOuKH\nu1lbF08P+Lo1Pxn2GQDAt3bZji0XNcQ+V7dz4zdsV7CtdQN9mOcwwtorrfnyoV8PPb7b9m0AQM2l\n9nnVddbe60ZpY3sE2zOay+cJAI/MmAoA+BJ+DAB47puT7YvZ6nGTG1FzrWd7rjUz+4Ueg67/KwDg\nx/gSAOC8R39vX3wljvLgQmvZhLv5YXek5lnW7uS5f6nPDwEAV837t9jp02yXzefGPLbd3UjTrZk8\nwdqrrfnqOJvQt1/9fux6sTVP32rtXO7e03E0TNJkf2TNd0ZfHPp87eXLy+fX+jw3HnQjaVT+F0ZP\nLZsfAFw4zQb/yRa7yT0usv1L58Q+ut8a7aud/L+7F+38fyX6s/+852cBANdhVvhu+Q3jbGMOd7Tq\nm4WIoj/7erZvdnE09X3GmtJh8aurJ1m7zJqnZtrF/Z/zBoQuO3oeCAC44MJfAAD67LI//9Tr4jBv\nrCo7AtayXbDX2lPmxb64x5oTLngYAPB4z2MBALefOCP2aeU9Waa5asQ33ECc9JP807da8/ShE0KP\nZ2vt90ePt2swYojtr30tjqI//etst9r/Ef1ein3GYGlZ+9xY/ukb2WFNfeyM1RpJX1qzYkzoseHl\ndwEASoc22Q6eAobEUQazPZCt7vD62AXNPI9evP7Dxr5qG027Y6dGvuSXaY4a0QMEr8K68in/bZxN\naOuQCBT9BtvYev3qG/8E6pqGi8pJ14+Os+/buAEAsG3gIO7RiAciikbifd+I8nkC2ISB1vbvDwBo\nqNti83SjaMQ2dC0H/J3vs2TJ8v8zqSjSS50Xwi//+bj4pVSXJXortrItuRH0XtV7tugtqT5CR/3e\nYwX3tbZYu9HQYPm2OJ/rLrA59u5p410489f2xeuhC6YR9Tdtt1YKrtqly2PfMX+09vBprwAA3j/W\nFMTHTn5f6NN231DbWHYk97yQnAMQ4GjzDmsX9QYAPPupo0KPxRgPICI9qLo3OqRPcbm0y9oj3ZyH\nrjKMOHTEy7ZjLL8YpakMdfPSoyNlklrKmoj0UlJWEenbRvQCANS+a2fokuJyW9ICAHQe1E6GjP0b\nAGDQ8L+FLhsa38Wt9BnxSM9tISkv7XrOYn23qMn0q7Mv63pytF0dR9uhjS1s+awMdg9NXZ9NAIqQ\n3i0/0/ltRvk8AWxCHb8y7bShnx20tk/s02N7wZAFkpE+S5Yqk/ynz5KlyqSi6r0s9MFo56yRWEKr\nTLCylgpGkLol/VIGL69myvgli5SUba/e6zsuBZbQSnr1iNBjeV+b443n2ZyHDDGd8vSzopW1N1Xh\nIx+wNhjF2L6AKGMWcYNq/rFjH7e22+Ohz30tZ3GDqt46nddaRJGiy3NYbAa8DQvfFXq8MMls8rvH\n/g4A0IPq/VBnWEyVXR3hyFWxD3h+B49YCQAYdKjp0xuadKxG11kmJF1nznON61Ky5lUMK2trh8Q1\nRX1/zm8LysSr0djElhd8CPX9enePNzR0pt574f3fxo/rNOxgHqYu9GwebCfST/NzyzxJMLdqsrwE\nA4J+DtRRR189UHt6//35acXk1Hup9WpRZwu23k69PzCr91myZCmSiiJ98MHP4Y4lHn1ThJdhY7Tr\nQ3/lwPqyj2VdlhF95rPdKKPYfN8pOZbm0ztuz7FjPNVkc75tmqHdmClLQ5fRx9vb9UhC+nJ6joTw\nHp+X8rsx9ECOe9XQbeKwRaHPfZPPtI3xdN09JI3G6wwJNq+gq87xB16edCgAYGX/4TbPETbPpv6x\nz0FbykcTeMIZ+0DUb+LGCLYR6b2D6CC2uqcceU977LLGzutv9NGpHeeQvgeBmXapcA29a2w3h+5B\ntBV6DvbanOgP3Xkt93SBpDLBbbT7v5HusI0O6YPC0AXSy2G4g3PvzXYA/if0CagfvcOUovntLm83\nR2PfG7zeAel1bx3S7++fOSN9lixVJhVF+sCya0W6gY4IL66Y42JN51ubTLoDTrZFyxH1EQlfWG+/\n23cfX3k3UyuYe7o7lrhKqUvMzaf17LI5PzD5JADApJ6RLDT6hGtsY4E1hxHNU3cYEC0LY3RIAvz4\nYYtDn+aDjZ21fCxtHg/JxuAZecmCcTORdFkkW63AwWWtkL7H4NAFtQmSaj26dVPs04+orzXzO6UG\naCnv0WqzR30goqiDxHV0he21dmM3LmzdvAKSdjIaALRx7vWc60DqKQPhJq81c1/Nryuk59lvNqTf\nvJfr5W7/0GFeQtKuRnuTHsiI9HFNf5COpWun+W0r+utpMc97vi1qHkJ6tQHhe8ZfF82xSDLSZ8lS\nZVJZpA88eqFlqaCTFuhE+DMdXZj84qlTHmEP47ZOCix4YGH9MQCAez95CgBg/ugP2BcD3Di3SXvQ\nmqmIBMM5PmSU3S132SLxsbNbQo8TRhuldvREQ9IxNEs8RXDzSK/V5g6uk3vzkGNOizYC0V0D0gsN\nNnuk1xzTtX20pL+6hdbx/tZCS3CHqP0cCQdwdNddrg+BvX6XnVB9T55Yl0gvTS1BKSCQoDZvtB9u\nqidy+dMjoqZ6g+e4a8R6IekuQ9IBPSOihrmFa1iEexqVa3r+/I3NRNG6vrHrfiC9RnuTNGzQet7b\n6SkHyVXQAemdPSmIns/EywBgx17rv6Nb77J5oVfsU0T3KZKM9FmyVJnkP32WLFUmlVXvO4TGevOM\niDf0w8lod1HsccYUC+WcBSO9n7ic+vSjsc+040zVP7LZ9OfrpphifTtcNNtmjj1XPr8kLMrPsXS0\nta2mLC08e1LosRC2LfVeHPehHK7kRguRAFT5xlDNH9UW2SsjavkLeeqkRm/20WxpfBfnvC6q9ztL\npi+/xmixEMXm1ftkNN2ZIhJMn9ctwrBuWLCclbcAUEoV3oIRqT7v22i66Ob6xN3kJqbR0vkB7qnh\n0H222PwOGuyiEVP1ucOIflSOSAPc7m3W9426g2LXxFCmUYpU6DA/qfe7AisfvXtyW2p4UMe7MuTt\nLpsfAOzgHHf0P7BsXt6Ql9X7LFmyFEplkb5DPLx/FxHeRLyhW05GO8Ah/KNEeMW2+7h1bp84ixvH\n8YhTIlrO/yiNe4u5b6Og1cdyaY6iuxrJZ93TI0OPpRMYQaboM3rYhtKFVxT7J51iDCPEav4a+wyr\nVVw4X+mNhIEl3tKVRponMddAoJNuGmdQvJMI38uTc9imsXE+cl9IpaixAcNIMhlI5Onr719nIzpt\nToYoGcz0G4/0fctHSUcrG1HIx3keGOPcIsIHJO0K9xJD2Tbr+5aHzcRQVvRHSag0YX7v2Lkv9HlH\nz13l8woD+fmlc91TNh4AvLXT5vZW/57l4zlDXibnZMmSpVAqjPQltnpXewRj8AyX2SLeyC0HuDU8\nEX6NxZOE1EUAMJ37tMI9cZj95vnmI0Kfx0+2rDWBwHOXAndc2piA+pzzCtJ5l8UepQlNAIC1g+08\nhg6x34hCUZQLZWu6w3FX6riIrm2wtm2AAm68A0u3KFmPOm+VtkXRfKOPwV6v/tHnk67pdUf82jlF\nermeevS1dnffInRKXXYOo3eWt7vwDhvHr0MVr47OJcxxV3nbE2/FTh2QNNUdCuaYfNT8yn7GVGhd\n6Q3hjOm66+HcoGGO+7WmT3QHdyn37enGQ3Qrm1do/84cvWSkz5KlyqTCSL81+ezf50Q1cnNErfXE\nm2ClJ+AL4UtuFO2bqXU+fzOpOY6jsZ8bPbn82IX4kmRXcXnKlGFFgRlD6wy+azsJDwU6Wnbh+oii\neVA3a9s6WJ6Bju9vrmPdek/bbxGpduj3PSPSa5T0hu8p+sC2JyG1W3fbsbvsx+lIHVa4Hcbby9/s\ndT/dX3TabylcM/8d2dPFdxX+h+y37OlW/rmLef29M89InyVLlUn+02fJUmVSYeWls0gsIHDIGQ+v\naDlx6YFIvJF6P73IkBc6s6XLbiHiOBo7GuWKoraTOYuI4tLUK1NLiO5ioyiwotGCoi4XkHNXhcip\nvXRlBW28LJo8GZHca+eq0faB/F3gfjuDkpxbqSZbqLH30s/pJqK7qGxJ0WGkAvpKYrzqxt/0dPNK\nbWpdSmK82uutWB3siOl189K9rIkf98YuBXbJ/Z1fu7uoYY7/m/H8zem+l4fYWz6Om3LBAqtQMtJn\nyVJlUmGkb2Ir15ingnAfM97InaZoOSBSa0W8kVtupifnCOFZR+OhZttxL+I4wVUXkukUzUeYzDmL\nv+Oy9DTRhDj0dfrfGJVWVIpDeBd0HXkrHTVWOdna1tHpF9xw3gCqV7pG5Dx9xBu1kmAY3E6VwRkN\nUxddMppJn/JWmsg+RYRt850LnX7lIwrp+2p+NkCNz+W2q3yUIiAMI2o8zm8HXKSa5ha0kSJo1ePe\no/xjX8tR8A7vAkzcg/uhN4SNXc4lKW0pzCtMx89LoyfUYW/w7GVzC0iv+b0dzYGSkT5LliqTCiO9\nSDAqBFVy3zG3jHLaMeNNiIdHDJ7ROl3EGx9wo++E8CqsMX9BHCdk8NkoTFZeG//+lhuPc7b6EWiY\n8Eroofj3kJ+OQTRFSC+ED8DOGPf2mMQ2ZIdFiRAWYnE8PVij6lZx/e+DX2h3kM2hlwhALrGMdBoB\ng0YrQ3qxjAbr59yxkQFLnhAUtJEuRgwx7u2cOWfhXZud2EP8g5naRXYnmggAR6nVjiILS6LfUAM5\noK9ZPHp7Wq+0kV3lo/kn5sCkFai/1SvaNYL7NCEqdV2ijb9xdpuevXaVzzHRRID9B/2M9FmyVJlU\nFulVOTYUk/ShrNrmQls57VzGG4XHKnhG1FpPvJGVXmv4gPBXukPNVYZWLep1bI9zSbXbFmsmuYKa\nYVsJbZmNxmfBlQiDmrROZnDOitoYErtChgMpHgHp/XUStmhEBg25FPR9m6xIonLbqfyTp/xKd0hX\nj54YrQ9tgw1iQnZYzcsH+QTdoYsRqY30bdjIjxzAaSBC+nQ0h+HRLkLPx+b+BtGbvWFDc+uA9B6b\nNWq5XeSgAbSFeBvP9vK2CEX15wmWBWoOb3SLs38TiT0kzK/ISpBoSy6Rz4F93uQ3O8rH8eG3BSMW\nSUb6LFmqTPKfPkuWKpPKqvcnslV9+DJDntxmYswwTXVIYomQ8Ubx8IqW6zoFNr+Y64ouhBTYLmQO\nQHQpAmBWHM25/6lGun+fS5M9ehkz5lC9X0r1uageuKL5ezO7jjJ8L0Ws6voyDrENGQaV3rpwRCm5\nNDiOit809SlZG+tImTj1Po2CCIsFVyxBGXdUlCIYGgvV+85GdOo9DYyD+9hEQjVXX2CjrXg0zx8P\nlVkTA6MvQxXmFuaYjuhH5Ryp3ivBZp1fd2iTanSR6qzRDtL8eIu8gTEsQWQE/V8YGgHnjtUSJDE0\nAtmQlyVLlk6kskjPbDgBYFUfHkB805XYCr2dgUM57ZjxRmgeo+Xc2LLRbUwMhGWdNHYTWzcfbXLO\nJ/W0KpUn4OHYR5tPWiNdRZjiSccyXYYsOxN5KvIFAli+8kjtpMii14XuIEOnIw0dzN+plYFxq0PU\ndETNtfcQt5PGxlW8PtGlyO/3eONTiqRhxLjLhsEwqh5qPdK3byoezRvyeifEpteoiawP+hSiNhLg\nritSNOdITURlsuq8KpMkRCjLMJSOJqSXh9NpIAHpNfTbmZ+zU6pUVkiprQvmiE5dOQG9ZKTPkqXK\npKJIf8w0I9M8NZNc2Y3uzaxy0UFKbH3xRq4BldMuZLx5u6WqJU1seeyxbj4zy+d8Jm4DAIxesDr2\n+YM1zyflrCR+Vs0ESSH8c8Nsca+MugCA+UTtUIyyxNa/s/X25+jSHKLCgDE895GvMfifSF9y6z0B\ng1aN4cyHxT7K7ruSBoOV63ndgymk5Dqn2Ec49tWs+XMVxFQbbA4A1hJRU5Sq8x80WRKbhPSvwakp\nwd2pxXjRmp76g9bKgdT0elkLIDw+KvvVZTCVgqhqNYPInFJxzKhEyDpQ5LLjiPpXOgKWSmWF4pic\n1w5HdPp7gTaSjPRZslSZVBTpP44bAQD/c54tTpZvGxe/vJoLyFAuupVtyY0gtNYCS2helHe2M7II\n0GENL4T/dOzRfN5zZXM+/bUH7YtbY58dZAE/z8/pWj5m5UNAeLzfmj/hHwEAj+86NvaR2WGd3v6l\n5ByAzujBjeNinaojNCPZBgpIQ7pKmutwfdEc+2w9yvQAeRj2LeZiVcpT4YiJJuK8Chhr3ohD8RcA\nQO0yMkkis7lDQU2N5mtcBm2Ej0yJ93PDy47T3IHYpPl5PwD1B4VLN2l4Uz3q9zrtkHaHTQkN10uH\nYCoqHq+72b++ndshA5Oe5S5G1FreIX2Hop0iNeU1fZYsWf6e5D99lixVJhVV78/fdT0AYEdPU9qu\nu2BW+G55X6r6c6hqqz6847p3LJZR5MqSJGpmdJohEG9a+HGmNVLpgVhYQ3PGHH4RM3JjHlUpGfCk\nOEqjHeNUZan1L461YhmtPPiWh1wqnlZtyHip8/MqKc+ngfvorZwQAgCAiXjaNrhrBy9bLKAVRYrn\nCNnAjorfPd/NFijBrahDBPXej6gliEbkyTsD48iDLSoxRCfqNF0FXY2Ycu69PRCsN7Kh2SxwL+NQ\n2+G5Vh2WIGn0n5urBh9lyw+p9/1eccsqqvcy7RUtFgINSZeA19S7EretGWQbXar3geZjjdR6dxEG\n8z8weEtb2cR8GEOaeaEzyUifJUuVSUWRvs9VVt7nggt/AQDo3TO+3W487+MAgKea6D4TffYh59JS\nMcmA+CW23h0jc0oT2yRaDoh0YBJv5JaT0Q6ICK85g/n45q+Kw2gWKcVnilK1HeemdZo19+MkAMDD\n20+wHfe5Pis6K/DpC1iSwKOiINNN3TgWj4cezcuIl0Tml6iReL1IV0NzDq4/d7mf4odF62mFfFLf\n6CI4l1YQwpE0kYnxm4mc0Hg8aztoaNzqrmk6os68foTbSSKSDIxLdY+XuD4dDHm6tj6OkHOlatZw\nsE0kJTUBAP5mTRqd6AlYwVxHhN860q7Bq94PWmLbAemLMvpwrlIGGyOVXO7EHpyXIN47TtPMC51J\nRvosWapMKkvDnWNNn12GnhfO/HX4asgQWzTdNm0lAOCByYaIW+7ya16iB4tJhlJTPvBD6x8trLWm\nbIldFDwjaq2IN8Et5+YqhF9IFHkm9gj6hbBD+XbruH7H6bHvQ0PKc/Vtu5lru4fcgMF+kRb4bIpd\nBhD7mPa3pf4xAOWBQKGgJ5G0q8LgRyp/G20Dr4yP1/txsPzXXLrqAtJr8VykYRGGifDSRICYf2D4\n4g3l83OkodSK0aQvvEkmoTC/sCmlLwMR6lO7j3P+deecqeXIlTgK9gyWIT1zEvg1M1BOD24UKYc2\nh1e7GcKX/P3rkCtBIxbF+dOlSIWktinylTtQmAuCvYqc1UWSkT5LliqTiiL9fKLlVJWYdgu4088y\nlB0zxSy7k3oaKjx2dkvos/BsW2OGctECHFdqKqx/CDjKaecz3ggVFTwTqLWOeCMrvdbwQnj/JhW2\naRk8RuvXs6xZfnw0t96GM228pz+gHSZrPD04xWR5HmLOfiH8Aadu58c/AgDevcotaG1XoAcXZfKR\nItRPc6b9oRXvC30e28ttaSMbhUoeAiWk9/QaUTZPaSIA0CJtRJrIoo6jCZ2Cz0WaiPMCrB1vupUo\nzLvn9ysbz0SQmtJ8nAlcQUocW16FMbuWdpjYmk7Cpr2FoEdCGlImpFVFSL9T6/Oi0B09WdTH+PNh\n3SJfubOsSF2QjTuVjPRZslSZ5D99lixVJhVV75W+8g2qzNOui9/1pio1+nhTtUefcA0A4ITRMX5d\n6tzSCeaqUX14T35Q2mdljZHK5tX7kPFGQzNaTlx6IBJvuoqRl9I9SWoi1fqt55sZ6rc4J/S9ZTvJ\nRnO4IxjwPPlIirhUUbkbXYzZqdacUG+TPwn32w5HGsICaxQTUOT4myS1mR7S1dPMsPgAXYoA0Hbz\n0GSuuoNScv3VoDtVdcXONF69zz9w9DJeTS4/nqZK6pcfqQEvLD9cEOZ/470AoqExhmn47Egltilp\nyFkEOXb/ybY+VMxCn8V0074Uu2q0LhYLIepPh1iJg61d7wIQAmlIIxYRzMoz+WgtNsyFI4asSFTv\nd3SRtSktJpdKRvosWapMKor0QksZxTa5iKAjH2AraibRavTEGKUetuliWTvY3og+M4mijkKpKdm3\nvJFH23RByeD1vOuio6aUDsddiQj/Mfb9lLXXdzvfWpwf+m67ii46kY5CSJ2DkyBNbMnAOTl+M+hc\ne7WfhjsBAOMWUUX6Q+yzkDYeIajQ00f99RNykjR0FzfuXH9a7BSKgog9k87VoeZAWq9IePrI0Ds5\n/F2xzx3W7P5va3WrvStRBrwpyj7Dea47Plb6/CPVidV/5A1o1TfeoZrmTyAmD3C4JwpzT3sgUvqy\nN+TpWqbPwyhfJp6XY8NYowd3iE4ECujBaXQiEK5CEz/yNEfFH6NpC31+vDVrt5eP5qW2YJ+XjPRZ\nslSZVBTpm9jqHbfCfSfv3XKi7mFsx/jilApgIagMHWJoPrTOrWTkVRJpQSDl3trKWivcEqp7d0e6\nttT6fYyjlWoNL4S/rv951rKU1upfu8R1c9hupDsooFIRZWaKNVofz4w9zib8niGfH92MW/8Y+6Rr\nea0op/r8dyQOzZ9oa3HZH/Zd51ApUISpdoUrJDyeHvtaHRKMPPdFAMC5+C0AoPkBF5Rj4I9Wxn6n\nFgwgWAbQO9FE7kXMiiwqc1AiFmst7zURnb1wjnrOVNdluv1OFOZ3v061kEi/ymX0SfUGWQjqRrqd\n1ECF8C/omJ401KE8elE5D2olfHz6jjUy0yF4OfTooXH+Wj5aUY4lb8spkoz0WbJUmVQU6T/M9c8C\nVtf1SC8M0TpP6PuUI/AM5fZQAo9W8rVxuYc2oojezGq9hThJbBrEWzmFjgqeCdTas9ycaaXXGl4I\nv+SGd1uHq9yACoIJqFlkAac+MYpv+pnWHD/h7tBDCFo7h1lnaLVvdVRWnaswboq+cMv1rbNs7nNo\nd3jqfkKr86gAj7DVndJcW6w50QUxXWSIpZDks5bf12G854mg0nGEcd7WMCHxhDw50UKuRW4CgHU3\nEF6DuUDXtMgPQPVwAK+pU07ec3BCYZa2xHn6SF3draSAeVkWYnkDnsYEAMCzm44qG8+ks7Dw+o7b\n1BwO6WMIL5owgKjUJEVTPZlX93+o7A57USgZ6bNkqTKpKNLX0Zh9CtfpS906Wwivd/XWpAWiZzMt\nCXygywD6ZidtkSQ1YsoQJyTAUHgs18CeWqt1sKz0YQ0vhF/sV4Jz2aYr2aNjl76GECGpx7nlefoA\nYPICJvpQIBDhyK9mhXEaecTx3HA5AK/uZoaIX6+mQeJyflHyPgxhskakzjCWF+erseeFw38KAPjC\ndlYKvcKaNXfGPq1sU1vDh715mZ6QtTNtp67tHxZ+JPaZw3aNDDZFfoCm8jkrnPrMWOFRHIL3tT1h\nO/hcrufSvoSOIjw+WnN2Np7VY81Ds4g7Az24LBBID30aWuRyJyfBYspuHJKPAOGGb/pr+WheNMU6\n2XJeLeiEjPRZslSd5D99lixVJpWNp/86Wxqfxjg30xgVgaQKUlSiQhRIKXFdZciT0iQV3jtEpKKJ\nWpIWogAQctrJ+KV4eG9QErU2EG/m8ItgtJuLKKVkZowB94H+TBnY/5J1/GhWsBmv3h770DC2imQc\nkXi9YqtlylRlw6Faf+vYyPL5Kb5gG1/lfFp1pVvdSBqVKnIjqUnfsea8Kb8MPb+ObwMAen3DPm/i\nPO9FFN2vJran0ztYE1MlYvulhjs/w+cAANes/qx9caUbqFUpwuXP1dwLaMFjuY+koZOHxhmdBGOE\n1YjCzOGKIiq1GJPmXaOoP0cPFi04FDAR/8pbrDvQg/U0uoSKHLvHZFvcHsVMQw3L3DqW6n2Jxrm0\neIkfWfH9Wb3PkiULgAoj/U+GfQYAcMIFZkA5fJqrckDUH0PoGkPbzA6XP61EqqG8eHq7FRES9M4X\niaLJcU5CueikmGRAd8SstSKCKONNiIcHIrIHuqoMLV25kITDH7RmputyiRmZPtXzVwCAC3ZZLkHv\n+ttBMo50iBQ9AeAMnTQRfv5phnrfxb+HPuu+xNf/TULNO5IRgYCWA3nORNsZH7kWAPBj/FvoOegS\nK6S4/ory0bympjmeSxdrjwus3fqtiE/f7WZz/NHaS22HNJGbi0qNl9gWGEUHUpOaac1hHzH8Phu3\nhC7vXkSLHZF+OY2iMrV595eKgUySe5gknw3TYv1o5SJ4ZeHhtiPUTHUGuA5ZAHmzujuXHZ/Ho+oM\n4cfLEuhdfypVxo8pFcnPOSgRnujmJCN9lixVJhVF+m/tskX94z1t7fP+sXHNe+xYo0KOe5WvML7V\nerv6lWOI+mOULUSg5AJ3IETXK0/hjj6bqsA2KSapUlNAzEuvrLUhp91tiBJCTvVKT1eDnmCqNTwR\nnrRVfDUi2IVDfwYA+Bys7XMFQzznxFHu4LlKhxA+nN4z9hHCL7vA3vXfwDcBAM9915X0losOv2Er\nTPZoadoNrrbmk2eYyvEfuz5v87twX+i6imt4eeikhXk36BmyndjPsfpiu6bfcxrINS/yS7kD70o1\nESAukhNVSTW2AAAgAElEQVSE79USu9BOMOhie1hmwrIbn9vm7CPU0HY/aq3sI0W0Kd29HmI60Q06\n16mHj8k+o+ci5BT0DtU0SJuWJW9P4jATGAA0cS//DM71t4l/k1Rv8EgfahmMRZeSkT5LliqTiiL9\nlossgd3tJxrMPXZyzMd2bDdD+onDmBt9mL3WxpwW10Oj2swqXpPkBYMzaoYywVwqtRPpV9RGUo0C\nIpRNVdZWX0wyVJ5R0Ine3mU57YQNaW47vcUdaur1PZMfifCfOfinoccXYNtD/4NYQ/S83b3OhXF6\noysMpZ8j3qz9pn17Gb4PAHj057TaOzJNVB+kM3CuDad06HLxCWauv/zVr9kOgvFCR7zxfgog+iRa\npridNAE8cootiH8AW7c/eruLHf4B20Uy5qSkJiBeX9KW+3KBfWHs0X+2eUAuwM8BAJ/ea3aSmqvd\nMFQe5vL56axSEQAcKU2RCL98ij1PD+BDoc/yeazSFDRAWQeKshSKjGPPog8EGj7JjAsKBOo3n9aF\nZ2Ofl2i17zR7MBDX8uPRpWSkz5KlyiT/6bNkqTKpLDnnahqtWmsAAG33Rb7xfS0WVnXfZCO/NB9s\nHHDPNx5RWwIADKs1lkEdg+cPcmmE3yANZxNj8FRSaIVT1l7GIQCA5Stpnplfw9bNtZVtKDWVqvJA\nR7VNZjWpnRPiVyKg0C0no51UegAY+XPm8jaNFPdQw3W2zKDY0hyIRo679Yro9voSfgwAuPvXJLIH\ntXeOG6nElkSSURzx5tjjexOMwPPlBfTVXcR5JdFyQFxuqL5HIyn9+Fbs8/3BNsAPd5lav+UyLqF8\nNOIeRfal+QYcN12Lh1HUX7m0abw4BnN8FubuvGCvXcx+V/A+xjAGtPL6KtpAd1pPynTn5g2pAxn9\np/j+B/bGnILByBsMeF0lTqcBT0Y2F/0ntf5YuX7ZtDtDXimZcwHFJ6j12ycLy/ehSDLSZ8lSZVJZ\npJd7aBnfcsuOjF/dR6Qab6i7fOy4shZAfAU3GVrWNhDpuzmk32tI37aO0falXtZ6KqTy5i1OPq/b\n4ToJX9Py2EVUoCa2tFol8fBApNaKeCO3XDDaARHhCVjCCe/4EyA00+W3kyD8uW7/Efr87ob/zzYC\nvXUO25IbiXOdSOLNXXZNbxw6I/SYcRPdWxdbc22SvbbMHSfE+oo1j37sPQCAbwXuNTDvBoa6yVi3\nTMjszYC6zjprWQIdKYoZgaXBTHu/Wc5mueD9Ga9x7rymikp8xJG9OtMlWtj2diCu6L+7h5glT3Ts\ntjlOAwmZhqQVdmXA47PPG9pwQiSqvQ8W5z9yMTU/ag4vuEclddVp1PqC+P5FPeUPfApFkpE+S5Yq\nkwojveBWb0C3Wl2nHOuE84foIxkQu4Qk442G3m0D7DdtfV2fbWw3s12TtACwWYQYzaeUtEBchwkH\nCuKeQ8gO/S1JTjvFwwMRhUStDcQbl1lGa/gU4T8Yu+BIIk47gf3CPv8JALjhN5+JnTpFeLdwnG5z\nHvSI+T9v5WJ12jciGqzievz3KJdz2TZfEPdt/Q+7Ppd1Mxj/5cIv2hfeTThXKK5CoZqXDxNJaMqT\ne1vrgnIaP2kawjnMIvQJ/BcA4PAHHK37Jmt2kGI7N6ljAHREeF2dRiF8VHrw5BTTOG+ElVR/Yh4T\nLTgbCEo6P1kJ0hKnQHCNiid1qo4do8/eL81Hu57uOPdOKD7AUa4TlaSFIcNjRvosWbKg4kivN57W\n4CX3XYr+7LvZvSU300a5RPv0nvOrXr2/9S4UYntSTVvS6jc+xELoI7toE1tXTFKVZ8QtmWmNctr5\njDchPFaW6jnWeOKNzjxF+CNdXr52/v7TtbaY//VvLiw7NgBgz7Xc0DXlSGfGrP3jbrWF4t0MJBp+\njmVcfeR3cRiFDcmUMkPrRWoZ1x4fofCL238CANg2S/n9ZR/xekKJrc5Q8zkxdpleU3Y+w881osqZ\njv+soJkQMKPCo67Kz/MMnklX116n0HnJWlAvphMzPD1zSszrr/yHt788QztMylhJcv+kR3P1AfrS\nvs5npnmaaYO+EtDIhVzLkx78fFLHwIuCakJ4uCNDPTcipZd7N0mUjPRZslSZ5D99lixVJhVW76U/\nKbmwt66lianVllwfqUtSD7sn+4Goou9hW6S6p6LxvJEuca0MoJrvbGEywqjUlApRKE11SGIJBHVQ\n8fCKlvOeRC1WdIhgtHNamdT6a24hAT5o2NfGTkERJFVmlp3DcdcGnxIe2mL3ogep4zct6DgfLS8m\n0WC3/CqzpJ6HGwAAT3ztuNj5O9qgBa1D2mw3Yl+q9TO52xnpjhtnc5Q6r7JYDQ+4AAup8VR/lyYu\nTqBjWTLNwrsZp3OVWKNwA7PR4cnjzGh3nZtYSCCqDD46zZAmHEhS5CAuIFriLi0FZ9jMTmFugFCI\nFAAT+qBdZde42yeJ7WDAk2HQZfKRWv/43hhTUiQZ6bNkqTKpLNKfSDfck2w3ezKMHBJCqaKSFGli\n667QW+ivU/JZ8vSelJFOqO4MLg38vaKfCL8HnBqD91UuWsUkVWoqFKJwRjHltEtjxrwzR8gq4o3c\nckJ3wCH8R7UnNdoBwD9Zc5FFcP3TTw2Zb1kyM3ahW+pnzJumK/rvrnBID1NY8J2TjJ3ztYcvLz/2\nZhebHsyQuqZW4guTXSID0njHnv1nANHl5o10zfOo/RHthObrXdYYIV9RHkWJ7rBsjzIZNvqoP7nm\niPB3DzPizfVUQe5+8WOxr27BHO1oZetJ0qkTkAeb7DRR3tuTh5sGcyqfndpbY2pu2fQW8C9QZMCT\nDjFGl5cIr4xPNsMWAOV09yLJSJ8lS5VJZZFeb8tWtot6x+8WMzhlBdtAoPHvOb3T07y4XmPQmFqn\nC+F92SC++QbQPSQKqY875hrpgOmG7C31Ro30JAqtw0K56MR1tNDVRUqzsmg2sSxjDJ4RtVbEm+CW\nAwrW8Lo+58Y+l5mr5pPfN2PAdfP+FQCwI4Z+46dUWJRO73NEirbWXqHPcaSDPncaL8ZdOguRa/w1\nJRlnhmHsAZfbAb5QHxb7gaA0+h6uuBWP73K3LU2CjNLiJ15SEnST++5oZhLqJwaq1rouXcCLk0Zy\nGpbyWNTa5+7n+foSX6GEltbwmqGfWVKAtICOPfZDpuVIu5m6hJYI525cw8ucZmnwVzvYJjpk8olG\np8f2Ml9FiO8vloz0WbJUmVQU6b98qAVfPH0oC/x9KnIGNyxkihsFvywjCq+IGW+wjtsb+VlUW7cc\ngoBK9N2BbN0wYUGkBR8RvnFcDM2cwCR9CnNUkcN3r1oS+oS3M9frKhetYpJFlM8mtspp5zPeKDxW\nwTOBWjvTDdQZwn81BlV+5ttGlPnFPbYWX8OKUB64lNPnFBql777aoOLUBZEkErO5/Cw5C87r0xF7\nmn9pnoof4Mt2fvOoDfyrOz9er3t4fbQm7wrFZfNw1pZw+5qlphRliJE1m0j4zAgbwee0exiW//DR\nlVSBbuYzJ2rtEj+zVrZd4S8P1sBsODP58VORHvxxUobP2ULSkuw+f4ijpIQiXQt/DUboXO0U8GTz\nuLJzAoC226jRpmmNEslInyVLlUn+02fJUmVSUfX+u21W+ujZWlNUFjt97IVJZpp4edKhAIAVOBgA\n8OqWYaHPzhKVPVKTu1TvpdYzOUvfpg2hS1OfEgDgYCqYqgp6RHAGARMZ2tSsElUyNrlSXCKnLyQ3\nWr8u4nnL8JIWolASSyBmvOkQD7+niHhDtxyNdlLpAeAXd5pav4rcnBu4v8WN0kLv2zcvtmyVs7/5\nQ9sx+2nXS0UlVPnVWOojH3wRAPAbp0pOvoJEpB9ZcwtjCjxdJc3yInXVE2YadbubO2n9NgdY12x+\nRiU8BYAXOKoqyCrSbPnTLj+DDFxSf1v1ha6Bp/ukjrPELQdEtZ73VjkUzncZiz7BeIwe2sUlT2tB\n8lNJk47k3Kky3O1k2TUVZXl4fbwnIb5/RVeu7Yz0WbJUnVQU6WtYqejo8S+VtQCwe6xZNFb2t7ih\ngPT9I9K/Ns6y928aZzC+mda6t/CO0OdAGliUN6+ebr4heC30aSK1V0g/8jWqDr6O+KKk5XfPuyKA\nKbKnpYV8Auy0mKQKUShNNeBy2nWIh/coQ/gm8UZuORntgI4IH+i0Lj/cx2dcAwC46cP/bDvuE690\ndeyE2daQDvr7SeZgPOtHBiFrLo09pWfI9CUU/1dXhKMf0SnY0uhGWzc+Qtg8mKZXIr4px+FrGBL6\nxH3v5GczAm9YGZ8VLKFRTnbX5D7aQeQWFqKnBLGusiS1WDPaGfJmWlN7if1+Vjcznc5yFOmGOaQT\n0737NOcXdcyOJkLFdfZwrGc9Bvf2MV6vyq7tu9kl9guuugXoSjLSZ8lSZVJRpFfpo1Bux63TenB7\n9IjVZW0oSwVAL/udXBe/0cdS5uxw8fS9hfTbLYVOL62V/oooQuvlnbQAdvCl/1KSccVjbppTR64k\nrfJGHO86J8UkVWoqFKIACrLWlti63DkMnhG1VsQbueWAAoRnUYeTTou02QePIlQsljuO13BqLDF1\n/J8sL8DDiyyyqJ3pA35YUCOcFgaMUBZcagG3jojnJ/LLYyz0uGEBb+6v3EAltqlbdrPrk2ZHUrvR\n9dmTkrvU+iAvPRw6Ia19i7IkKZ8jCb3iwLjsOsPPMzaW3HIiIw2/KdqTlGJhKcFXuWyK8uVKUxwh\nDvHpsc8zk0yXErHouXkiULmBNupcvfO4o2Skz5KlyqSiSC9bcC2X141xmY2htI43qYyxrNyDYx9t\n9+qvlq/8nttiHxI/QqmrTWx9+lBub+XxS/yNR3G9I9PcOl605pL1eao0GFpU4Yg3t441xFO56FBM\nsrDUVIkt4cRlvFF4rIJnRK31xJsWtlrDC+EfbHZQsYLWeuknl5wBALjyx58KXT5/oa37b2FGWVmV\nZXKodyWi7v6UqTUn0Hy//IpxHSe2THRpUVj/zPYN1ykNjS6SzqzRewr67Ek+F4mwVXlomtg6q8wo\n0ruluDDoaOqkGFqrjD7nMOtz7a86Bl4t53OuVXaa+xeIz9MkacK8bW0zIkVaWtNtm6wNob6tbqBA\n8ykKSYqSkT5LliqT/KfPkqXKpKLqfZqO0juHQkwc1fJatv2ccU1K2EHJZ0+CkQIpdVyKozeUpCkz\n9blIhU/TY/pKppPEn1cEF9WwrbPsV1d3i6ryT2ElotZ9ifHOoT78HDdiiS1NgUxTrSSWQMx4o1hw\nRcuV1cfl2HLLBaNdUOmBYJi6ysx9f7zAsqscd9QTocd36d5SiMK/M5T/jivtN2e8+EDoG4yyK1q5\noWMVub3SxKY+q0AqRbkTtC9dAvgnIc2upM91ro+O28SWd1mcMVdJVrk7Gz9kD+RJDPg/DbF074mr\nqLtLnedXz7tcAJ3x6l3ZF3wgXSbOtOa3LpLyFpwNANh9Ha9hMOD52mwpzadYMtJnyVJlUlGkl4FC\nNjUfw5RmxivKFpLmwkk/Ax3NNvtjxtE4Ppub7IdNbI9M47OBGKPNN/L8iYa3c5hDOeRVA4Cv8ig3\nSRdhia+yM6XBjqWmVIjibhcErpx2ynijeZ7iDiVqbSDeBLec01OI8E9fYCSfo+vNrTPbGTxlsGv8\ni7WHNBv0Lz+MRrplN8XOAVV0l+mzc9mtg5tLRCVRpf0NlPtNVOuiYiVrkj76jbPnBmq2xtaD4KMt\nSdEOl0Xzoo21+fCY41DRlsqnoJTVg+5xB5WlmoGK83mPPJlX2mWK8B/0BmshPG/ATYPNyKpCGwDw\nym8O55fcsVEqcVEmHx8J2FEy0mfJUmVSUaQ/g+u+rSTKyFUGRLyThy3NjQN0nue2SFItwLtEUpuA\n3oPDXZ9AIEre/nBUyNXTrLDDXXw1/xbnAACeup8qwOWxL1rlNiFTJrhRvFuIdBoWk1SpKRWiAGLW\nWp27Mt4oHh5wwTOBWsuzp1sOiGv4FOFnR7YrnllrBJBhv2G58EMUiDKbrSMNnWksldqb7U5e2s1I\nQ8oMDABDFxHnVERSp+VVLPJ1VF755Z5Gy/XBNNpeKap2AVV38y6jaL+x2e72vm3mcuvfGN1XdT3t\naRNFOwZg2fke5Ti7yq9Q+yhVCAVeuaw/6xWAxc9p5L0/VelDHxDCu4ImcvXe0WzXVzn7nnrYpbqV\nK3SJ/jGpExCIT70PaeooGemzZKkyqank4O03wriRegU6y3x4+5Mws5UvsPVOG9B6X1pAEY0jRfbU\nXgwA9YxJ6C1gSMM5gVgIkMvsV8bbArCVFFIAeIAm9DvXG9Lvu44D6y1c8mEUrWy1qiPCN7ikbQyF\nvHGCES5mfMNINY98K3bR+1xZa9/YbISNugUOT6YKkRlxQWrtlX9yxJujzLI/m2AmhH9obUSTD36B\nc75yDvfwGBMtc84xf44w9yAX77XnGxKu509a46w6FPtOqxgA0UahW3GEctN7sBqftNTGdrv792p/\nu1+vU49TcFY37HXHshkNpH45pM3uTU3R8ligT0v8Vn5+xj2faVFziV9RS6+bpCy24ku53P93jzat\n7Vew+/XgPHbymuN9evqVvsklZAyii6ZnrKbw/52RPkuWKpP8p8+Spcqkour9H9qntgPRUDJ0laPM\nSKVK1HwXBh+tfNuTtki/V1ixuBhev5daLxVLBsajIrnj+W6mGj1F/f5xmOErpBUG0HYzo7CUSDFU\njRIf2ztrpNRSwRtIlct5vb53ghF4vnyT5cBeRQ/NDbFLjJxjFaTxJxmZ5rmaya7X7LL2+HZGy114\naujxXfLpWZICr7eb0W7CF5bGYa6Uq6/Jmqtszg9e0AIAOPH0qN7fQiKKVFzFp3nS0GG8JwfS/fk6\nb3/J9dFjkEa0e9pNugRQ7H69q6sRojO1dFMNCJ9lSc+P4jT0OBZEZi5PMgGtSboCcbmZRt77a9As\nl29CvPndkBgmKZfvHxZwn1LHx5ogiGq91iA6uqeP8Rk7kebDh7J6nyVLFlQY6ce1P1GG9Ifi5fDd\nwVgJAGgi1MuNUr8rvkv7vL7PNvRmVuuRXoFIQnrCQtvgGKH0N0L9Kr6LV/Lt6N1Cyt+3aL29mvfN\n5YC+cIC2N0o9kZlNeODxiflPBrIaOiPULj4jFoO4fMHXbIN2mx/y1H1U9wwWk/zOVSw1dRqtO3f9\nzPUi9j1px2rvZrf1lnfHHlKIziDxpuYpxp/PmOPGoXr0kCHGW5NtnPVML+4D6MS7mUqDVDu9hjfX\nRgR7oYxsGp+D97GoBgAMfYBwq8IhZPrOd5AqbEsNZt4gKMUuLVzmvYNpFP2byX5PHnujYB9QrIFI\n85giI+T7XSeCryLmRK2VuxcAnniYfmFFMQaK7YOIIiNxWkrLuVEnkol0CT9/NCN9lixZUGGkx+x2\ngxMRXsbGrwYdaguoEUT6dwrp3aqpjov6AfgfADFLTk9Ev8ku2ILxDVJwNnFRv9EFWojMoXblekP6\nfYtdfjEFSTyZtBs3xT4h74lWssIB4YnLlNpI3x/XZ588gyWnXnXVIIjw1/LYQpxLR8cuy1+yt/ch\nD5PjeaKoIK2x04lGgf39gxacc2adGRu+50woCp455EpSa2ukLj0VO11lrrm3Zthj8cyA8iNd6k7v\nufm2wv4AU8tu+DgX1TfHPtij60N61Cg+bmfGLgdcZAvtWfWmR3yKaXWO/oPL/sI8Ae3UAlTo0eeH\n6SwPQlFR8/2R1L2Y0rQB4Gi5grVuF5HLeWWfGW96gOLhFTgTaLVAjMGaqxkK4YvKpwjhW6wZ5fyW\nLBjafIHRiZfXjM9InyVLlkojPYj0CnrwhsampG1MWiAGaAy0N2CPvva269Y9Lurf2mlIL9olNvKU\nigI2SmzFa/CRiGF7VdLJM4rSZHHJW3ese+ty6T7jI5YZ9epdxrXs87F9ocs9tIDL5q9AymZXaerY\n443/+cQ/EEY2z+Y3nwl9RrZbsriVPzJV6oeMfZnpZrqgneGxhxEul3GcibNDnwf/bOcxtsas9AoR\nupTr9kuu/Xboe8VZTAF0m1BJGohfaWs7LTnuXSt2TIzi+v8ya8Z9MoYXf45lts5vsxjWGmX5dZG+\nu3n4Z2j3KXG/v2Mp5zylbvv1f2obCPn5fa0pOVCU5fc4Y1D5UlP3smRpIHTNoXrgY5cCtbaVrZ49\nr5skefebynPuA0DDpVZO63OwQqhfqbkyI32WLFkqjvQ30kScvjeBCOl8v6oA5QDXRUjfN2n9K1t+\n2DRjqs+Uqu09enOW2BZlySsKBJYIC/S6pw37RGKGy3933pRfAgB+DAt7HfTPNsGFzgQujAzVcGip\nv/aqmHL1X75GWPuOMttyfp+OSP/EL83z0Fhja7mA0C6nXc1U3oqxghjTSo5pj47shae3AABmUwOZ\nTVC5ZL4h/BXN7gRDgg7e0wbO+ZLYJdBmVVFG577RJ35wESwAApINbIm7OGbjpaZ1KeusctQBwOgF\nTNGiwpBSFLyippz8jkoLIHp+fEUZLeLl75edxYVarx1vGst/470AImXbI/3q+/lD+dxlmd/sJyYv\nkJ5HPaf+/9JiTSO1SSJ87WXxGf5sNyNjXIBfAADeWbMlI32WLFnynz5LlqqTCqv336BOWZTzprMM\neN6ccmDSp6vcOTISpbF5fl/6XVH0fjo/H3UvSi35n9LCL7LfXjj8p6Hn12Eq8aBLTK1fdYXtj9Hm\n8aw+R1fm1sV2fkN3RpVtW99B3JrN9osAgOb2VaHPX64wPfonVIPF+Hy+Pcbcnyrr4AqOc6a1m26N\nhrdHakzVl1J5ULupkuPPIqPnNp9zj0ubmycAAP7zbLP2Xbjk17ELI9O2fszO67vdLPrvR09/I/aZ\nyXaJVH6p+/4et1gjNyiLhPS/aF3ocXZPU/WVy+59e40A1O/ZAkddao+V2u/Ve1K31w42Ff5llt96\nwcWqLyRlW+0rC+iGC/RsRELXYhXj0Pl5d5xoR3q+9cy1xC6juBymUTUtpQUAnyK7Z+Qcuy415xf/\nvzPSZ8lSZVJhpP8FX29C2KLiBKn0KNiXZskrkreTJU9SROSUgZEGmF4uqkPcUxY+GHmulXCWYekL\n268MXXsRzNYT4WVc87qFgl8aCf+f/Zh1/uXHvhg73SzDG7WTGWbAu/3Gk0KX0+vNJBiy4TCM/pCr\nYyaY5TX/wy1zD9XuMX/TplnRRzp7DluiyeBrzTi2oebFst/avEzNKZ1tFq/hR1hanJuWxC4yiaog\n41RN+frYZ+ZgMzrd8CsaJoMLqigfnwzCvBGemCKta6bd96nDWwEA4/Fs6NITb9l5EVlF9lKZtDfR\nO/QVuSsQukTdXh+p2/uepAVQSora6G1EpM92Vv4UiM+hDMQt1kx0Wu9MaxouMLecnrnz3cUceRM1\nH+6qeTQjfZYsWVBxpN9CpJc7zLvIOstC31WponQ/0BH9i7SC1G4gVPcuEaJGA38n18x01+VMW/N+\nZKj5tJQP7qzlXMRdEbtu4lLr90zcolVbixuuhYj86NXvAQC8f6FlYMXkHa7Xj9iaXnDAOkPWvcv6\nhh63+EEBnG1ggJo72uPOS7QeN+bOD9uNDnxezVWhi4I3ZxHQD/g8f3/TbGsbZoe+//k3ruGPsDX8\nbCK8564oxFRJY3WHz3F9RjB+6NaLrX7UuZtMJ9p9qkO5+coMJN+fRvJsrxZrmhrLPsJHIMsd3MS2\nL8rFZ9eVm1ePrpQNn7BGWk1JO4Tm3h2nag+pC7gwv441yig8M/Y47GyjcH2cnGQhfCiFDQS68vJH\nrT2kk/93RvosWapMKpoNFzP5tl7BddCauB4Kb9BQYrir7PipFlCE9F1lydM212xFtGARSYjwB0y3\nQJCW+hgGqtznp5Fh0fwAT4KoviYWPwnoJoSXzbfFBa3gW2q+bhuB+/J714nXcLLZFr5QT36vi9sR\nCP0rrdChXLSPhdU1JIpIS2l1PYQ3ITw2BM8ou27sKyv9TQnCn/3d2Gf1V8zz8MV/tvV+K+fjk4RM\n4ZhnLTFtafz1ZgH/6J9iFchnvsTSM5drhrKQeB41EbXEmzqHV3yO1z0IfLr/nggGlCfcEOp3IHZ5\nfvfqZN/+ELv00LlUGw2ckHKezLTGF8tUkcxP7DI473Md6dy+WCY5Pq0FR/eSkT5LliqT/KfPkqXK\npKLq/aDrLWZ+w8uMtfZGkBLbNVS51tGwsdEZONLyRVK/ijLnqJXKNtD1UTmjJrahrFE0dI082LK6\nTGRg/SRGjbU4ZenoZSRUqH5FUrAw9oxmSh3yDEVp/Vvs8/3BFgA97wbq3HO1GCi5kZgZhbHSctVs\nvSf2kOLZj1wcxW7H+vBAUNFpmFQhCp+N5sP0QN2ujDchHp6/He860xsohZZlNoNKDwBNL9q3377W\ndPivTjJL5+B/jsNoIbN1jrVnvG6q8oLfvzf0Oe/H9uXvx9PJeRn9c2sWIkpackKq/1zXh+exhsu9\nNZ4IlkoaGVhE+ioqtgmUE8xkLNYCj0tcbyCmWt9rht2TM/tbnMWZLkneR15lUIEiDJlp6Pnole1Q\nLLMzyUifJUuVSUWR/sf4EgCgdGgTAGAVWyCSHpS/bv1eFinYGK0r+zYSeoT4XSG93C8DDL37NsQw\nu8F9DHGGwbLPKC/fofhL6KP8bSJzDF/MGkw+CIyljXb/t7Wt9JYoHt7TLXSm54raycw1j5wS6yH/\ncBcD33+gPYq7c6ShvubGGXv2nwEAo+8x49E9LlIs6EbMzfZYKNDhqzcQ3ZS9aFXHOStrbRS5OtnL\ngebW75gB9ZiPG8rJcCmjHRAR/mtfM7/c5m/bvb0cXwt9/omoL8Tfwxj5sz8Uj3XL9TMBAE3nlgAA\nPxpP5tNlk2Kn+0QB0iSF+B6ZpX/9PSzcX0nz6hRZiOkKlgtY7riYqBhHTzBWjwzEKod9+LxXYqc0\nhyDvX1fFMjuTjPRZslSZVBTpz3uU728yWdtGxAy1KdJv7GaL8E31Mbfd5noWJSTi7MI7AAB73bS7\nEfXVHbMAACAASURBVPYP4sL/IL7ZB7qA+sEoRvraZc5HI1BcnLSLYpenXy/vmkY/+/f76VRSejBG\nfvXFttb9gUo6A9hyGY0Ny0TmKLF1CDbTmnMUqkM7gndWBacUM7hsWKAk8H92vbjOlK2DgOwpTMpL\nL60n5LRbQaR3LsDvft+CZ3540mwAwDwiUKvrozW8EP6KX5tPctisV0Ofz79m5bamE7yF07c7DesM\n1gP44ZV2rKaJdv++d++/hz5rriCi3sTsw4u1aPbazvqk9VpAZ9KVK1jr9SZrBtAl7G0fUuw4nZHT\nXuTHqDYpSOikLbZu7yFUd3abHSTczGPu/jRLI9CxWGZnkpE+S5Yqk8oWsJzEApaqMPMu9+WQpNWy\nyL9I+5e3u7nm3Ov0k55c29ak1Ut8Etu0eo6Axi2ZxJzcyrXSS7vKdgPoGCqRhkl82M29hkErW79n\nWPqlbj8GAFxz++djJwbuYM8vuaFV2Vdin8V2i/4yzjSj5pFm3Z4dI2tDLcTu7Xah3vkrGkE+PcfN\nnpNbbKla21fbuL+MaeoDXaSx3fo2fpkX8Qc/S3oAWGQQ1j7MxllFw4In3vwT2zGWJhD/MetfAABf\nXP+T0Gde/T8CAKaeY6vTe0g28WvVkFlI62KaBB495T2hj4o//n41Sb63UYfxSXrkPSqx3SnvTREF\nXHeX40hJbXBdmtgqq44QPpptcNjhdibHwijW8gZ5pG94lA+tcv4xCnq5C16SZV50IGmX/u+iuzOV\nnqKaVzMNN0uWLMh/+ixZqk4qash7kDqJNHcfV1RPlb2H9JNElQcQLRN0x/Wgel8WVyeNTC6stAQW\nELTmdmqra/m5q7SYaQkkf9y0WOMEqXcfi323X2rvU2WLueZFqvU/iH2wR9zqhKE/PWplx40zTnrz\nPJvhUqfWS5SeeR6zu0RuT0FOAbk/udQa7L7SUmYSS02pEMW+H7Twm3tj55mmw858weLh51z+WQCR\nSw9EN5zccjLaTfh6tI5+D3Z97v7lGQCAU7jk2up4N8EUx5+dQaLSccufCH0mXWg/mD7c1OaHL7bk\nlAsuPjb0Wfccs1zKSCti2GYa4IoKo8oVLLXep2in5bbxULty0e0bGTMieUm9b1jAB/OPbhwaLXfw\nnBckxjqgY7KfJra+WOaRcseexdYlKPKSkT5LliqTiiJ9Wv6grAwCX3j9tpR/V5Qhr3vy2UsadyeT\nTFcZ8or6pJiYojoQ365TVM6ILjK9WdfOjGaVn+FzAIAfraWLThF0izxUy1ylMyPldmbsEaiYNPK8\nkPwCQOB/lDTD4K1M8xAAYHKV7ZMP4E9j8Y1QLpEuo1knmf/tmlHUUlY4pGdOO2W8+dDFNkFFywGR\nWivEl1tu6rJophPC/7S/Jb77t29ZfP+M8+OhrqUBVuf+Ji/h6Q7J+iy28/jnsyzjzhmn2HVTyXEA\neGGc0YuXjjMq7Gt4JwBgM7nbbyGyk7rBEiH0hlGZVWKt3hGXm6hSHUwHqpD+yLZo/q1RFh39Gfi5\n3bmCn0rKdKWuYCDeb7mFZSusP851Ypm0nTP5OSN9lixZgAojvVA7LQnspStCZFHu286kq9w6nYlH\nSxFOZXdoYnuYo6b2k8tICM+0s09OHAcAuB4Rnq5ZbWtcfJVncZeCX3wAiK4Mg+wn29py+LkxMknU\nTJCcoevlXTVCehGewrq9THgsGi9e7mnr/6Nq48rxDt0gahUqJnnNZUT6WT4ZABeinzbMOfdMi/dW\nPDwQg2dErdWZt7kYcK3hhfCPHm9uuOO+Etfr51JZuoFrXRGT9BkAjmFKvQlE1No7jXh18pRHQ5+T\nx9v29iMM617raf5iFT3dVYb09kQpj56KqJaVUv8rtSQBuy6lc7Vpewe/e2F7+U+Ajq7gIi1TRb8n\nycWtRMenxT6rTzEC2G9DgbSYs9FLRvosWapMKor0CrfUKsgjfbq+fjsZ8roSndBBbl+PZJ+Ivt5y\nLaNsvZLfinHjKZUqVHi8uRhUnFChrH9Y6JguesneLAKI4nG9bqN3OamjZNn4kMqGB8zosZ5rQF03\nnw9GSP+amE6FSM+rSaRfyhDPo4+ISN+D4D2fYDaV5aJVTPK5y1ricBu1SDWI3X2qhbv6jDcKj1Xw\njKi1nngjK73W8EL4Vy6ILJiRr5sh4lym+buDnpqSG0fp5V8ihB6p1hW51HXq02wI3TzELkZzf14U\nV7m8U69QEemLZK92iyTHCvegl5K2s5LaQNTelOP3mG7xuzrl+hPCk/k0b/Qxoc9tMPvILYH1lZE+\nS5YsyH/6LFmqTiqq3jdTXW1Oue9AUJN2U19to/rko4ZS1b8r41zqFvSuv1q52KQ/yVo3zHXqpDqp\nKpMCsTrpHxkydT+sesO6G/jjOW68Vhnu5OaS+cnPrMUaxlY3ftJ0Ul+NVZFWz6NcfERfB0NeSOXs\nr1i5IU/qPcbfHnoMpvot19hUZmn53PHGvf/kJS6fdFD16ehjmuqQxBIx443i4RUtF4q0umPJLSej\nnVR6ANjwTWPIDNppJ3YeI/nucGr0iqTVImphtLuhkdtDeXw9Bgfx+eiYTwB4k+r9mwzIbNsbv9Ph\n0+Vr0TJWoqWmN8Q2sVWMwQg9g9NcJxYK2XCKXYt0aQkAD66kVe+6rkNqMtJnyVJlUlGk30k7Qi/x\nUf7qvuRbtwfbeiJ/vafP6jWpsPe0rjgAyNihKKiuao2rVbSfq1i1odneoELAxbTgLXSx7SJ6rP4j\nX8WqNa52jSfeKCC8xFY6SDS8hLJMLMiomPl3L3I+H3qchGAapdlbIWnVE9mk2JBHyxQBdCUOtg1n\nqJQBqZVtO41g59PH9o1Lvxn6rrmK12WNzHJ0yF0eiaHKaaeMN4qHP+Mzblo0UArx5YaT0Q6ICL/9\nW4ZRffqbIe6MG2Ofp5eVj7M2af12qhX24DG7OxegJCV/eaPy29FAhexNbL0htjHRLpUBSegOAI8O\nNlemEP4eWLTkK3dHF2lIWX4XupSM9FmyVJlUFOm/1Mde18PGmk9jyNi/he+GcIFfR87oQC7yB+yK\nMNVnC8kPegN3hfRE+N1Kq9c/1iwS+UIuLbWBtopYingp38EvbDI6xO75bg3emrSh/LAWqUVZWvSu\nJwL2dcHWLNY47f3mcPoE/st23Bq7LF1ePlqYsavduK7Z1JpXpcIUIj3Kvgvr/7HxK6FPK9sFXJxq\nbT/r8zEtzuwLCcWXKUPNb5IWIWutctop481nvhaj7hU8I2qtNJo73L3WGl4Iv+ErXOOPiHWoJvCa\nHS23YFv5eEDHfDlFZLG3I2khdT0pPrBMTtlRsisVuYLFeeIa/s8j7KbMdSlzH4YFEM17mUn25NWN\n3l1gsXSO1i7nnZE+S5Yqk4oi/VXzmOS9yd5Ag4ZHpFfggsoGB6TvGWHqoMH2Tj6QQQ8qNexlL6F+\nB0tWKZ/eZlezSEi/nu9gIX3Ixw/ErCpaThfkyENJyK51bBoi4SkXeu8T4Xu1WHth7NF4scG4ctkf\n/gA5qS43mieyAA7pHVLIDrFhJdG7MOCGdmOCo67BbqcxiJh0oLIHcf9UZaj9fPQq/PSiLwAAtlzH\nAVbIn+CwVXnpmbVWOe0OPSVmIVZ4rIJntKYvuZnLSq81vBD+mY/FlfHRk222Nbx2E4j4E5x5ZBMv\nb4kW+JQoU0SYScuhet9LSt0eyoV7jc8QpesrjYr3rd0pfAtrjcYt71ArsxnP3RSRfvddPLLimcRG\n2ukJvWl+nWLJSJ8lS5VJ/tNnyVJlUtmqtTRUodGUow2NUe/Z0MBtUayVmtlXEtW2bHK90FGkwW5L\n2o2uT1prPG2BqJWGfWlhciAqnTIJpXFRTa4v1fqBjI8iUan/7Eg6+Sws68yM10iQYaTY8678lxS1\nsFgQgcTxZF4QrWMJSRmhEnBBRDbdn5t32cV9tb/juL/L5lZL9T5kD6LWOHpBVBvPnmKq/jUzGIE3\nuyWZMRDUTRaiUJrqX138qdBDGW8UD69ouYcQRXdAbrlgtJsc4wZuGGFk9CM+b8bUo0/hd0+GLqij\nql8nz6rIYnIT+6q1IuHIUKxn0Ov3CuJQ5JvIXm7JpFvzygi7zotxFABgUfDPRVfw45us3f0QD+Iv\nQivbNSL/P8XWq/f7Z5rMSJ8lS5VJZZF+GfMPL5Opo6BmfHeiUyhL5bqkCN89aYGI9HpLp4gPuBrj\n2qG3pS/fuDbZp7dmV/l10rgoF28+lm/rmdYMutiYSRfg56HLBXu5zWYHjVC+LKOOKP5GiOl3hwqo\nEYxWOpeCpG9s3thsZqjX66ODaeQwQ3rt0SjPEAkn/SGOdtIUs+5dM5N5A+YwTrHkCcJCYhJ3WIji\n92eeE3oop50y3ige/iUHYEJ6OUTllqtxBk8h/C9g8xk2wtzEE0ZES6wy2wSKr8hiutWenKNLp2dN\nZC+P9LxQbcPsAS2R7RWIT4iuYGljz9KSt/zlcXGc1qSVdlLyBUhT47H4xf6Z1OSGoyvJSJ8lS5VJ\nZZE+UFELs+RZs4ft5qQt65/m0PH5cFOCZJEDJo3eTz/736do7ueTFirkgm0AP5/oujKk+bCP2Bt6\nJq4HAHx6769Cl35X8FgMQZ9LpPHUUb27AxmYBI5nRkR31UJRewOodZFljae5b5u5OFU6DEAIOhKh\nRNhS0hzc+vh9ex8DAEwd3goAmN/CnABzfFElYTTRSaWmbov3T1lrldNOGW+OdEifUmpFvJFbDohr\neCH8z2H1xJTxBohIP2rwSgDAkMEiiJnmp3x4QMyRp1Jqb9Il7F3BG2mIEv1ZSO9JX6tXMkMxi5aE\ne+RdwdrerOdS4VVF1KKu7Eh6JlxZtALJSJ8lS5VJhZFew6f5aPdXeiSfu5ru/yZLnh9fY2udnqI6\nEN6kA4i/SXFCnBnNvycPtZBahcme22YW+pqr3XAkmzySkGH8rISbjVrDW5wF5oaoDGD501wfhnTr\n3i0hqSs7nf6Nhhzd4GJFOX1pF5pHsAk79O33rF3f8ROttPf8yUJ6H0qirHiyi3BVPv/I0EN56WXB\nVk47n/FG4bFCeuGfJ95oHaw1vBB++a/i2nl5d27LQNJkJ9x/oBHCDuwZkb47r8semu/f3GVIv8WV\nUsc6Gpt0uUvJBIFI+lIbbo23uqeuo6L1uu6GLC6ynThUb+KdO5mfr0KhZKTPkqXKJP/ps2SpMqmw\nei+9VwpiV3lxukqNiYLvUkmNfEWpMdO4qDrXJ1Hnu1NVGu26yF0mYsx0I8G852Azap2gcqOINcdD\nbLxine8IXdBKtV4GM10B7/SartUGY6tfnGTWNkVdAYgkjhAb4NLFBOFA5OLU9TTj1WDvtqQhMS3p\nFR6SgswyIR4iaL0+a4uud5JjxpGPVGpKhSiUptoTXBoT9V4z3uSqDot4I2Od2qDSA4CWVlqFNZh6\nvmWAXZQtMTCzo3TlCl6XtHu8q02zLiWz9/eos/w6Pl6viS1JXwP43XTXhWp9/4/aRLZk9T5LlixA\npZF+8gRr9Qb01NjwxtRb8Y3kM9C1FiBJEb7IPSi4ZKHCtCghEO11gtkkKgoA+k+2E5nQ04xFKkr4\nPrIq3tcWCzQE4gjb3QSwuS4zkBwzOrskIbaNQ8MdmF/uTlY3eHTlh2KnUD9DOoOQ1TNJmqzh+Smf\nwUCf05lz69QE6vMZ8BBBU+DwZXbPNbruCUaXXB8aH1VqSoUolKYaiDntdEf1pJScDVLUWhFv5JYr\n09SE8EtoMFsihNV1KiqcJkndvX4m2lek0apPV3l2dNwOpVZil778TsZjuYdPjeNNG25VMU/C/QCA\nSwuOBGSkz5Kl6qSySK81VFGAS0B/oq/KBfusLym1tiugT0sLFwXuNCRtk+szytbDDQcbZBwKi/nW\n2hAAjiA2T4Rlfn3361xIquywI4toeznXr0WRznpHC+G1PKv/sOtEhL97mFU5CNlPb3Zr51ZtyOmn\nkT0dk0hBDUZFF4e0uSANbqZx5anj1IvKPhVe9zX9kt5EvZ3tro+dhwguKjWlQhRu5h0sBGXhJQqe\nIbVWxBu55QCENXxEeGU8KoqkT6WIvNUZoatINPuiwu1NbInsDbziLqhKiZOF8Mccag/Y+13N6+lU\n+Y57zTTOjPRZsmQBUGGk/+q4rwAA/jbO3t7r3dvtdb7xRGXcvNcgQoEgALB7G9+O2/jm04u1II4E\nfQ09DuhrNoGDBsRAGWXjUT6+elpOh6kekdsWAo6CrQnH7IpIr/DPDlRKtusdWUSra1EwhEoeNWU+\naGHbqOynrkzzk8epOOZMAMBz9/P1fzOcPM1Wa+cCiqZsE/y5NJgazxF5vXyukvCQ+LJPtOTv6HId\n3Jk4ZN1cTm9VlqNQagoxL32PJFttGT7LVsLJi1or4g0QrfT7V1q1MykidKV2JO85kodI+lyTNb2c\nFiT7UQfvUNRS3jPU7EfvxZ8AAO+DeYxatkf1spcITS4wqkgy0mfJUmVSUaT/9qvfBwBsHWJvx/Xd\nItLrjb6R7eZu/wAAeKMuOkvfqLM35luEFQU/eBFd8h30Fyto4iAXEqu3for09Xujj7rfK1yXCfnS\nFohLZu5bRUVBbueS65r6uvVeL/PBK0OqEN6Sx+KZU6LV9jpm37j7xY9ph8kSbyFOPf1ClZiDXlbf\n5sOfAwAcJbO5T+DL9XDqNQ6Y5GsJ9NcRaYuR3cUnouhU9nTYjPeYKoTTKlR5Js1LX4b0Oi776Dnw\n1Nroh+9MO/HFoeV56JG0RQXO0zy4DsVT+5G8CT4bLhF+0JS/8mO5d8hvH7vd1uu9tJSPlbiDiWK9\nD+YpkIz0WbJUmeQ/fZYsVSaVddldbE2/wabk9quLxpnmwdyWFtU/aYGo4qktmm1aR3x70gIxUU5a\ncdAX1NS2snRT1V3j2JIltp3l2PEiJVBOM8WVHelKaYWyRdTcn5xiRjup9ADw69XMJ6dS46FkUasb\nKC3YxKON6h270NUjNXGCrI+LY5flPFctjDqQQX0pLeaF0/IsEK88TbVTV5hTkXlPFe3XTTe0K8Z1\nkYios6d8vO4+ivDvis/sRFZUX14FuSJ9nsb9cgWzpVo/aJw9WEWuYC25dG/GvebWlrLXyffLqMKt\n7v49w/+AX5EWSUb6LFmqTCqK9E8za6nen3UuYKNfiuwpqgMxwENvV2UmLcqRtytpPdILfWih2krk\n3+RopQJ0oXaqHAARATvLrePpFiHamefXQ/Hwx7tOZ1lz9xDbeSOZOLe/PCP2EcLP0Y5H2MZMsFGa\ndFRrTo7fNH7I3v/TSeaofZSWL2f0UVR3Sk4N5q2Rse/awXZXQ3ksKXGeat0B6dNMSAikHhneAtnH\n3RuVi07Bv+zhTZ4NGX33hC/2R/y8uC2lS+40T91mBucDGuxhG1xvut+7nCu4ifrhwXQBH4qXAZQj\nvdzCfRZ14hIGgka2ibf9BSowBR7Xv0sVykifJUuVSUWRXnEgwdnh3t4Hvp60KG/974oy46Wit1tR\naE4a0lOURe9/U3ZYS9wmtke7JWGNXDIKkCDCL58SI1JUdljU2ifmHWdfxDqRIRd+XMPLx+ZnLyxu\nsUZ83o/GHgr1DeG/dPn4NWGabycEGyufuwteUZbXldJppCbs8Vcwdf7xCvp1sUJ9VdZMOe1cYNKb\nVErSe1P2PMgdl7gSlfEGgLM3dBbA5UbUep0IP/bcPwMAJrlcxYPl+k1KtBWRvoa/tsF2yL/rFbWX\nkpZ9ljubky5viW3qEvaz947HIslInyVLlUlFkT7NT/t/myFvf2R/Qh+6OpbWsaJdeHtuiL7lMrFO\na1wh4ETXmVlrN0wzCFJOuwcQQ2If2Gvm+7Y5fDeLWhtCZYG4hhfC6yp6CwINBqO5jyaBqZMeCT1O\nw50AgEH3EO5U0rnArtGBxCuukDs/5XFfut5CYmNiDK8v+JoBQJrIA0C4qELL+l2chYv4beP6NcXl\nfkUf2IrWW5bTLiC9rqGelgLijbQRzlUIP8upYZPajOhUo2QeqQcIiLn1V3XyGcAa/q6E8nZ/0mx4\nVG9iezSfz3/pxHGRkT5LliqT/KfPkqXKpKLqfZohz+fESd1fXeXGeTsquxQ0f2Kd5dTx6mGHxNdy\ntQ1znVR/U+puEh21euyg0FUpnVVr/DEa2ZbPcznbbmOrmuMl0X3mu4PKhKMrJLXe1bVqoIo905rh\n55murfTbAHDiKurzlpkb68nT9vakNL7/aLlPGfW1dnxc7CykW3Dfk+wUIgx9xgDNWVecIze5LqO0\nqwQA6PNXuq2cESutI18UyxaC2Xh5FL0Z0lTbToqevrR2lRtRqwIOI6OdVHoAqFEOOrnWZL9zevkm\n3tK1VLV1h31BtbSA2tvJrXOE80jWHcUNZcWOFdTKJCN9lixVJhVF+knKAKM3n3PDyDKxg/vklnnT\nGR9SN9r+5MLVG9E5akI8dm8hlwhB3koniFfZYRnpPG2WCL9hrBnnlsIQ9mlYLkBfflhI+MrCw22H\nMtb68sNPxt4myprnC1tJhL9EeKE7EEqC115iv/s4/gsAcA5+E/uwdJY8djqiRxxdOykyvXU6NEr+\nN94b+ur8glJS0jfekKc7lzg3neuv8VCjlyiPQWCbRK9XmKNGEx57U6bum4pJKhNPYbamDoa8At2P\nBjwRb2RorHEZeIXwa6g9SWvy1zTNotdVjp7OXMJATA58mDTQNAYfCPdpNz2/GemzZMkCoNIBNz9i\nW5TqmwuZ3lvK2zL6rLYVK13kgkipuaLu+hzmeoFr3adX6RDXh9tbR5oF4NVutphf4SLgVYJYCB/K\nD2+yxdTu+Q4p5idtQHVPnFQcvJBdeFCUX4cIP4o+rpmxR/9LDMJmdTN3ktxKtb9ywe13cjpE0K5I\nvFN0fYgY644zeJF9AgBeWUANJpgfpKUU5dznnAdQ/3Kx5Gmeek2s/a+xTxrQFOjBXlOj7SUWkRyh\nHVEC0qcjFhCpifSi1oasvz5IK7mWaVYDL525hIH4OHbqEgY6ZmemnWX7tIjbj/c0O9Jj4T59u2Am\nGemzZKk6qSjSf2e0xdbWj9bbMqLAAKa9Fe1Sn3u792TvXWbvf8dOs+j28HnXKe0KsCDCv9XL3qlv\ndIvv0jf4XlUY6CaaZF93qybl71MAicoNr3Km5pXrDXX3LaZxQBRWWW8dpTUWMVRghfDAr9dTxNE7\n3xeBbLFmMrFipjUNn4qLy/MZjTML1wIAht9EyufvQhc8zzmmaOTXxcqzU6PamIwuVTWdsqo68jh0\n0GA8jSTJF5RSkwGM50U7so2/pxdgRcckvUE055p3uZ1c9EobC+WifTHJUHkmpbqkGXAQ1B4FzwRq\nrSfe8HFOi0h7KrksMfXJ50YffKbzCAt3tgXZdVaPNg/RU7SpyEvkt59a+Y/ck5E+S5YsyH/6LFmq\nTiqq3n/t5csBAH0bTd2s6xMJ1UpSKbVeiSwPcqlXejOp4Tt6ml4fiiU62UtLnhIqKiXzm85pJx62\nWhE3Xt8e1ftta0isKXHHiqQt2hf45lLwvJFOA2lJUxR5kEasU8Hu66o3KiaefPqxH7JoL7nlAOAT\nLHTfMIfWUNa9X+6Kb8hFlxa8cqkz0Sz3D9X6Z8abnqlowNX3O19bcD3KgOdJORKdF8+Hav1hhz8T\neojTXqNlAtX7khtFVy6h+JQVuaRNNUT/rV5prS+WGZdWaUwALbw+JiAhDQX13hkYRbxJR/N8+Ba2\njVLVpbr7uY8u/65tvLkdF+Oo0GUxdX25hYNL+LnD4ziK2WhFl5KRPkuWKpPKuuxIGtk2cFBZCwCr\nSW8MdMe0BWKkU9p6EWNnZ9L6XG2qd7Axade5Pmm5YZE6fAmmDlny1BYVWkrJlB3y0CC+9pnTTkUO\nXMYbzLBxTh5uyfHOJHf3nC2/D116zOEGMxUtJcV2AaKkWfQIjJjkyUen8Sxm2IVWnP+d6/nFba7v\nYl0XIX2ab8cdRRoEedlFqZ2liuygvbOEjiIz5yiRrMbG714ZYTAtNyoWs+xXGdJr1P2gBxN9lfEm\nxMO76DhRa1NHqzeOBoRXMiQGWa4dHf2NcgGnruBwLihwC0szCoZUt73T04M6Skb6LFmqTCqL9K1C\ngdqkBcqJsohkGk+qEbJ3Tz57SUtdFSF9hwKYct14ZE7RWvaHNwr6dJZnx5Nq9L6X3UCo7uizQipF\nJp1qTfO0GNRxCiNkTiW7ZuoSroedO07lsJ/mevgp7i6i2Cor7wc0rdNdp5nW/BbnAgBuwdkAgH1z\nCK13ub4hPWtKGfblPHiuzMQ7ctqLAIAWt+hsWEA7BFHqBRKy/J1J4/t7F7i0tP59VjvlRi1LB5QS\nejunBytrrXLaBY3BrelTV12RLhduOxH+7tGWQumBkAoZWEQa99Itdr12LqotPwe/LbdwsDUVVWMp\nonFHyUifJUuVSWWRHg+yLciCmga6buNUtnkNIM2O19V00zw9HoXTrHhpQK/f3p/yw2lwZ0q9AIJ5\ntju/07rWEVOE8A0nvMKPlrgu5LEDcBLuBwDU3koVhqjuixS20kEgvUoo6UkiQvgPCtyYidel2Mfv\nhnwEAPBbnAMAeOU3tAwrT9/mIlRJaT7OH5BoMCqlPN2nBmKuvnYimY7g70wHyhLBvN1dS1m1l788\nTjtQPiLQkSLMOauY5H7Qg/2aPl0562nwxBs9BlrDC+GvWx8v/L7bqEmJ0qwwZW+P2JN6iEpsu8qv\nUywZ6bNkqTLJf/osWapMKqzeS1V+O/W/i2R/UmT+b1NipsdIs6gU5dcJDGq2VBO9DUuqotT6FmuG\nT4o6m9xVqjX+fqq9Ixc6X6JqjlPjX0PXVkzE3DG3jmZ8hOvzAUUU0vsmd6oMSwAwB+cDAJ54+Djt\nMFkio6Z3Akrd1CKCan1DTPEttf7oCaa3Kg13w6MusQLtgU/xESkyQTWxbZShjQGHC2tjFqLAQW/l\njqDee3ZVSvPhyAWx6So1pUIUUu/XFGT0kejp6OdjAjhnueNktAsqPQBczXaJDMxaS5TcQGk5Ylv6\nVAAACRxJREFUlv3Jr1MsGemzZKkyqTDSM5VHh/whQOd5cfwbLM2Vsz/GNcnbzZKXJr1Okq4BkTjU\nxFbIk8Y6A+gx2c71qLpnAQAT8DSAlJhiyDlyMZG9oOZ4O11ZC/iC78opk5rSJnmqpxCe9qM7mj8I\noLxY5h8WmCEvIM9cXe/WgqPqmsq8xkw6p8YevWa08dDm6ztpC62PD8Q+O6iypGZB79wNGksXmXwe\n35Qg/WY9a2kJDyBeqcPKxlV9eCAWkwylpgpIQ3qq9TQFHcdfdx4iEHDolitLgxgQXoZvUZqLDHMp\nBagov46u2PkFv89InyVL1UllkX40fSqivW523wUQF5oUZRHrSgtIpTM099t0B+qsPeVXtGAFXei1\n3eT6JGWH+441auYhfVSUMOajOQqG8IoXn7jXFpn95rtz0BJZ9ElTBvC8WzemWfOKYrY1xWPYjlA2\n1ALizU2DzwAAXM8dj85znF9ldw0kHCGP7BD++utitFhzYvlxAODM/rcDiIU2epAm7DySWEAyjs5P\nd9GD5QhpVIzz//MIU618Jp/dDyX01MJMPgnNp4GfSX+e6NgwoZS3dvESlNxouhqBHqwsTi4dgoJn\nRKkNxJslcKLnJkX4osh8PZg8SHeniSZZdcL9TCQjfZYsVSaVRXqtDWWM9mWMA/rzbbuNa2hPn02D\naIrS4XZPWlF1PZ1X20ku89AC8QXaaIEktU0Gt8O6dSxGOIoW4UNI0TwUfwFQXn64YRkt1Cl98tl4\nyHbue6GL9Xq6qktWowCAKcqQqiyoRHgFzgCRWqty2E89zIWxRwNlwwkIrxmlGfGBgPATibAzrfGl\ntBQcdPg8ZvkhsWi5Q7nOjnCMrzAt0xCZq3PJ+pm7aXrso1DfktbHstp77SS5ekLE6faAeXvLuNdI\nghFpiNpXUQZArapDTrsCenAInkk0B06abYrwPvCZLotGaquyb/hsuFSsG6bZ9V6XkT5LlixA/tNn\nyVJ1UlH1/sJplgNbiSg3BTdYzGKjpJVqd+yN3Psd20zNeWunkZn37aHOt8fpft0tqLlHL8uq07OX\nZdk5sE80CCorT5qEc6ArjaoUx/VU4KTKD3E5j5VFpWmLuYF6pLXGi2qOU0vcxPYll8a7xHZ/cuvI\nbCZnzAifNFH8Gqr1z0wy9VXx8ECMmAt8ehVfbfXqb6rW6xpK6W6JXZNU3IedbdF/vsDGR16li44G\nvB10RXpiURonoBiBusmuE9X6Rwe/B0BM0Ln7LudybdWGsvKIPOSJXU3W9K0vO533DDW13qv3IYiQ\nSzAtFvw90pyD0bGA5KOMN4qHD+p94NIDHRcNSWETAGjh/0K39ERbho47OF5NZSFS+0kUS0b6LFmq\nTCqK9D/ZcikAYFN/szRtdj6yNG9dQPpuDun7E+n7G9LvhbfumHRjBQy1vRkrf6Arl5nm3xPilyH9\nFsOcHkpxLIB38dN/r9b4Jte3REQv8XNKogQ6z60z3O0LUfjKcKOXv8tGvZPEm3v7mPvtTjJxbtsU\nkX73dRxdEXOBWtvqjpa65hLEaXSONHJ6Gi4wo9HHmZjvE7tujH20STLOPLrnirLpSZOZJLpwZAdj\nwylmiVWuvnkv0z94X+yDNTqfzqL/gGAYU3Qeh3kv/gQAOHb7E7ErAXRTF5l8gllQZdKofSlNNRCj\n/0LGm5Am3Uf/6amQVkItqtFFnPJWDr/A7pHcoD43gTSVQQvtOc9InyVLFgAVRvoeF1nbUGfuq4Z+\nLtCiju97vSX7JC0QS1TJ85SWsAKiG29X0nrXn8pjaUEmUPAFNV9Pvns9aQHs4PbahFCi1ZlH8c4K\nFhbl1hGehkSpvjy21rZyWxEBn2yOwSb3c9F7L9PYPjePP7oJUUS42SiEETPIOwhThG+xpokI72Lv\nVSxTJbTOx/UAgD7X7YuduJafT00odc8BkVgUQuOF8P8U+wjh71GaXuXq88VAQ74g3TBdaRcFNYBX\nnAh/zKG2cFfAU68/xq4hkw81thSL/dxDMcnJmsmk0CcU+hRpKMT/lNxIuiLF9GCbs63hhfC67sEd\nCkT6ts+bVyAZ6bNkqTKpKNIvnWOt1qq1DsV7p8ieonrRvo5L+ljUMkV8V7sx7CNCqzz2G65YZpoR\nT292b61NM+R1likPiIiQZgdscn20Qq4XzVRpzp3RVgj/4lhjfoiY4ktMPbzetvfdnOSya/UzUoTH\nC2yLrNsJtVYWeiK80B2IxTKF8CHnvi+lxfWr7Om6lj6YRvSTehGLaJ+YN/qY0EdeiFfuPlw7THb6\ndbG2dTea2EbUjZl8rM/7GeHUsp2mehfotHVx8aiennSEnkcisopJ+lJTIS99h/JfXdB80mxLiFZ6\nreEDwl8X+2j+T/simwWSkT5LliqT/KfPkqXKpKLq/b1sQ1pMp04fuL18AkWpM9NcNkX5c1LVWlr+\njoJ9XaXF7CqtZmeiORcl+JZJRu63EXJF+fCxNLsO1frnRsROf4JVIG2lyv3YXossa7vNKZpyXcmw\ntVEx5F3l19HsCyrkJnx6ueVmOV1Sav3ImxhYoVJaLrnO/pTSOlJGMBKLVp9i7q7bcEbo8+BK6vw3\nc8di3Z0imo+OwvNqcgQeBhROG25qvRJ09lJ8v5v7M1wSpmbBJnfEOi3HuARTfXiv3occoEG9lzmz\nC5pPQXVfEW4CgUhLEbckefC18iN0Jhnps2SpMqko0v+/ypCXikf8/9vMeOmYqeZxkOsTDJJsg8uN\nBp26Ia6zIq4E2gXZdbZPtnfuop4qSmjGK6E7ADy+11Cj7T4iu9DcZZHGCl0FQZXe9UXljTRrEXod\nnCQx8c1nW9GN85ksT+gOFBTLJOK0uiN1VkprqndJMhX3Th4zFtr4aOxzHUtUhTh/HaWI5iPdigY8\nly6g/0dNK1Fa8eNeIxmHbOH1rrhEZwa8o70xmYfYTSPkY4zvj/Xh3VRDqamuch7xCtH1p2g5O5Qh\nvYg30hy80a5IhyiSjPRZsmTJkiVLlixZsmTJkiVLlixZsmTJkiVLlixZsmTJkiVLlixZsmTJkiVL\nlixZsmTJkiVLlixZsmTJkiVLlixZsmTJkiVLlixZsmTJkiVLlixZsmTJkiVLlixZsmTJkqXy8n8A\nxMbTCRYAPO4AAAAASUVORK5CYII=\n", "text": ["<matplotlib.figure.Figure at 0x104b28490>"]}], "input": ["D = solutions.exo3(x0,W)"], "collapsed": false, "language": "python"}, {"prompt_number": 30, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["## Insert your code here."], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Once the geodesic distance map to $\\Ss$\n", "has been computed, the geodesic curve between any point $x_1$ and $\\Ss$\n", "extracted through gradient descent\n", "$$ \\ga'(t) = - \\eta_t \\nabla D(\\ga(t)) \\qandq \\ga(0)=x_1 $$\n", "where $\\eta_t>0$ controls the parameterization speed of the resulting\n", "curve.\n", "\n", "\n", "To obtain unit speed parameterization, one can use $\\eta_t =\n", "\\norm{\\nabla D(\\ga(t))}^{-1}$ (one need to be careful when\n", "$\\ga$ approaches $\\Ss$ since $D$ is not smooth on $\\Ss$).\n", "\n", "\n", "Compute the gradient $G_0(x) = \\nabla D(x) \\in \\RR^2$ of the distance map.\n", "Use centered differences."]}, {"prompt_number": 31, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["G0 = grad(D)"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Normalize the gradient to obtained $G(x) = G_0(x)/\\norm{G_0(x)}$, in order to have unit speed geodesic curve (parameterized\n", "by arc length)."]}, {"prompt_number": 32, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["d = sqrt(sum(G0**2, axis=2))\n", "U = zeros((n,n,2))\n", "U[:,:,0] = d\n", "U[:,:,1] = d\n", "G = G0 / U"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["The geodesic is then numerically computed using a discretized gradient\n", "descent, which defines a discret curve $ (\\ga_k)_k $ using\n", "$$ \\ga_{k+1} = \\ga_k - \\tau G(\\ga_k) $$\n", "where $\\ga_k \\in \\RR^2$ is an approximation of $\\ga(t)$ at time\n", "$t=k\\tau$, and the step size $\\tau>0$ should be small enough.\n", "\n", "\n", "Step size $\\tau$ for the gradient descent."]}, {"prompt_number": 33, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["tau = .8"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Initialize the path with the ending point."]}, {"prompt_number": 34, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["x1 = round(.9*n) + 1j*round(.88*n)\n", "gamma = [x1]"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Define a shortcut to interpolate $G$ at a 2-D points.\n", "_Warning:_ the |interp2| switches the role of the axis ..."]}, {"prompt_number": 35, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["Geval = lambda G,x: bilinear_interpolate(G[:,:,0], imag(x), real(x) ) + 1j * bilinear_interpolate(G[:,:,1],imag(x), real(x))"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Compute the gradient at the last point in the path, using interpolation."]}, {"prompt_number": 36, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["g = Geval(G, gamma[-1] )"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Perform the descent and add the new point to the path."]}, {"prompt_number": 37, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["gamma.append( gamma[-1] - tau*g )"], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["__Exercise 4__\n", "\n", "Perform the full geodesic path extraction by iterating the gradient\n", "descent. You must be very careful when the path become close to\n", "$x_0$, because the distance function is not differentiable at this\n", "point. You must stop the iteration when the path is close to $x_0$."]}, {"prompt_number": 38, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["gamma = solutions.exo4(tau,x0,x1,G)"], "collapsed": false, "language": "python"}, {"prompt_number": 39, "metadata": {}, "cell_type": "code", "outputs": [], "input": ["## Insert your code here."], "collapsed": false, "language": "python"}, {"metadata": {}, "cell_type": "markdown", "source": ["Display the geodesic curve."]}, {"prompt_number": 40, "metadata": {}, "cell_type": "code", "outputs": [{"metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAQ0AAAEACAYAAACpjCPWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFvlJREFUeJzt3V2vG8d9x/HfkudBRxIsGQaCKK5toJd+doECKWDkKjaa\n20SV5YcGfQnKjdO+ivQiaVMgdm6TVnUaoC6Qi15YRpIiqSzXTVIU6EMkx1acBq2ObEnnkDzcXuzO\n4XA4sztLLrnL5fcDEMsll+Ty2Br+9j+zsxIAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAVJY0vQMt8nlJF5reCWBF\nLkm6O88LaTQmLkn6WtM7AazI/ZJuzfPCXs07AqDjaDQAVEKjAaASGg0AldBoAKiERgNAJTQaACqh\n0QBQCY0GgEpoNABUQqMBoBIaDQCV0GgAqGSr6R3okj1Jn5N0XtLDkm5IuizpiqR7De4XUCcajZp8\nR9LvS/odSTvW41+W9CtJP5V0sYH9AupGo1GDPWUNxu96ntvJH0/y7UgcWHfUNGrwOWUJo8iD+XbA\nuqPRqMF5TR+S+Ozk2wHrjkajBg/XvB3QZjQaNbhR83ZAm9Fo1OCypEHJNoN8O2Dd0WjU4IqybtUi\nH+TbAeuOLtca3FM2DiOR9KC2taPh8XMDZQ3GT0R3K7qBRqMmF2VGhH5L53VFD+s/dUM/1WXdYUQo\n0FGXJKWL3z5MpTS/vVHD+3HjtpTbWc2Jmkbt7BEbaWN7ASwLjUbtzOUx9yV9pckdAZaCRqN2pkx0\nRtKfN7kjwFLQaNTOPjz5rLLGA+gOGo3aXbXun5P0V03tCLAUdLnW7oKykRmnJN2S9Gqzu9OAJEmW\n9t5pSnG5aSSN2u1Lup3fPyvqGugaksZSdL/bdZlpYt7PJYWsBkljKeh2RXeRNJbiuqSHlPWc/FjS\no8oakPWzSKKoM43EpIjQ55FA6kXSWIrb1v3PiB4UdAlJYyle0rr2oFRJB3WmkLI0sEgtw34tqWNx\nJI2loAcF3UXSWJr16EGJTQtF21VNHL7tzWNVkoDZtkoto2rCwSwajaW5Kul5ZYnjtLKiaLuKofM0\nBEX/4Of5nLJt5/lHHWpMihoRGo94HJ4szQVlU+/cp6zx+HazuwPUhKSxNPvKJgA0tpvakWN1JIuY\npLFI8nB/8ZeVPJbxOZuCpLFUd637ZfOVA+uBpLFU/6w21DWqJIwq64u8toz7S2+nBTc5LFI8Db2W\nWkcYSWOpqGuge0gaS7Wv7BBlL19f7a/WPD0goXTge7zsNWX7ESOUOOz77jImJcT2sJA4ZpE0ls5M\nymMfogDri6SxdBck3dT0IcoXl/qJ8/SEuKmhLHH4kkYdyaMoWZj1UMJwl+7nVUkevvcgbWRIGktn\nDlEM/sfDeiNprMRdSQ9o2fNrVE0YMUnDXfZ6veP12NfUkTTsZVnCGI/H3sdjFCUO6hsZksZKXM+X\nZn4N6hpYXySNlfDNr3GxtnefN2FUSQt2wjDroeeqJg9bTL0ilCjMunl/X+Ioq3+4+0HimEXSWImX\nJN3J76/X/BqAi6SxEmZ+jVOazK9RTw9KlXNBYpKGL1H4lkmSFD7nWy8aTVpWy7BTg7nvLs3nFSUO\ncz+WbySq+x02LXHQaKyM78pr8w8pX7SxMMvYRsIs+/3+8XrZtjGNh6usqDkej4ONhlkeHR15P89u\nMNz3L5OmaWF3rO/xruLwZGW48hq6gaSxMsu78to8CUOKSwt2srDX+/3+zLZbW1uF77VI0rCXo9HI\n+5ybMHyP+1KHb71o3zYtWbhIGiuzL+lafv+s6HrFuiJprNTiXa9Fv9JVu099ScNOEkXLXq93nCzM\n0reNbxmjqF5hPs8kDrN0P8dNHmY95nOLbPrp9CSNlaLrFeuPpLFS83e9VjnN3X2uqGcklDDc9OCm\niq2treikYZ43zPM+bhpw6xdHR0fH27hJw13G1FBC3MRRZVLiricOksbKrcelDYAQksbK1XNpg3lO\nNvMljrJkYdZ3dnamHt/e3tb29vbUY7HJwx4fYoSGhJtUYacLN1EMh8OZ71X2twoJTc7je92m9qaQ\nNFaOKQCx3kgaK1ft0gZlv4wxJ50VjcFw04CbHsy6WZrEYSeN0NIdt2HW7X023PEZbi3DThUmWZil\n+x0Gg8HU9zfb2Z9ZdqKa4Q5Jj3mN0dUEQtJoBJc2wPoiaTRi/ksbLDJOw5c0QvUINzXs7u7OrJvH\n7PThW7qJo9/vB5OGqV34EoZZmvsmURweHnq/r09ZwghNUhyTUrqaLFw0Go0onzc05rDELMuGkRcN\n4AodnpjliRMnJE0aBtNQ7O7uBp9zGxG3QYophIaKnYPB4LiRiB22XvQ5sUv7LNfQf5uixqRLDQmH\nJ41g3lCsL5JGY6p1vc5zWBJzMpr5tXbTgZsazLpJFydOnJi6by/dQxn3kCfm8MQkCzdpHB4eznTp\nlh2W+K6V4q67c3L4EkfZHCBG1w9TSBqNoesV64mk0Zj5rr5WNN9m2exbvpPP3HpDKHGYFHHy5ElJ\n0t7envb29qaec5fmvdyBYfZp9YY7mMskDFPsNMvt7e3gILLQCXExs36Fah12AgkVSV1dTRgGSaNR\npq6x3EsbAHUiaTTquqSHNLm0waNKktvFL/Eoq2mEEod9enuoi9VNGiZdnDx5cip12MtQ4rATQqg+\n4PaamIRxcHBwvH+hmoYRGpKepmmwS9e8lzuYy65xuH/fKomiS3UOkkaj3Pk1qGug/UgajXpJ0keS\ndvP18OnXIYsMI7dPbw8lDTc12EnD3D99+vTxY/Y2btKwe1NCNQ3fuAz3PUKn1rs9ML6T39wT4Mz3\nd3tPfDObh/6+m9aLQtJo1L6kt/L7WddrmjIFINqNpNE4d3To65K+NLNVTK9J1XEavol03BPTzNKX\nOEzCOHXq1NTSJI7QiNGYpGFGfdq9JpJ/jIfh9ry468Ph8Pi++71DUwbaiSM0HmTTelFIGo1zR4cC\n7UbSaAV7dOgppekZJUn1E9hiE4fd++BOulOltmEShZs0zDJU24hJGm6viX3eiuHWMELT/pn6yPb2\ndvB0evP+5r18l10IJYyuJooQkkYruKNDX292d4ACJI1WcEeHlptnur+iSXjKJt3xJQ67J0Uqr23M\nU9Nwz2SVwuepuEu3HuIb4xG63ILvb1hW04hJHF1IJzQarXFX0gPKGpBL0TNo+7aLHeTlm4081Hi4\njYjvhDW3EXG7YE2jYQ9fN9zipa+xkKavsBaaV6NoXo+ya7KE/mZFf+eyBoBT47Ek1/NlNjqUrle0\nFUmjNaZHh6bpa0qS8zNbxSSLKl2woZO+QjN62V2xMd2y9tI+ZT50eOLOLO4+PxqNjhOF+/m+4er2\n0jf7eux1Z4sOT4wuHHrEIGm0xsuSDpveCaAUSaMlkmRfafqW7Il5Yrpe65gz1L4fmkHcN3Vf2YQ9\noaVvKLjb1Wn4umLN+7jdsaF9tb+L+/5FycJe+sTWnbqGpNEqL2jS9fqc0vS1hvcHmEXSaJEsbfi7\nXmOGLpf9SvoeD02Z5z7uGwwVqh24dRJfb4abNEIJw/SMFF1DNtSN6vsu7mNFfxv38dhaRtdrGySN\n1rmaLzmBDe1E0midFyR9qMkhyuvq9WZ7UaqIOfYOHdsXjV+IOSHOXreTQNnUfKH3sMeWlI2t8H2n\n0Oe66qhXdDVxkDRaJit82iewfZa0gVah0Wilq9b9T2s8/mZjewK4aDRa6QVJd/L7t5Qkrza5M8AU\nahotlPWiXJP0rKSzStMfKU2fiD5d3hVzTF02hb9vvexyAO6V30NjMeznQhP/2u8duuxA2bq9z2Xq\nqEN0rZZh0Gi01sfW/XM6OvpL9fsvFhbXQtchDT0vzc6nGZpn0/1HPRqNgvNXuCef+YaGu/94Q2es\n+ubKKPu80L6naTrzWNHfxn287OLRofWu4fCkpZLkZXGIgjYiabSU/xDlcU0nkOJftbIrodu/9qHD\nAfc0dPvX3b6+qjR7inrR6e1lJ6yZ93CX9lXjffvke9z+Lm7Cib2KvE/XE0UISaPVpg9RxuOf0f2K\nxpE0WixJXlaa/lzSufyRcxqN/kJbWy9558qMTRb20q0HuAkjdMWzwWAwc51V9wSyoqHhZZPwmDRx\n7969qfc+ODiY+Vx3Gdr38Xg8UyOJvaarXdOIrW10FUmjxZJkX0nymOzaRq/31SZ3CSBptNXk5Kfp\n2sZ4/EMlyRMyhy4xvSihrki7+9J3jRDf0k4VviunSeWnt1eZI9QkjLt3s1Gy9+7dm0od9tLsW2jf\nR6PRTKKK7b4t+juX6VoCIWmsBbf79VuN7QlA0mgJ97RqW1bbuKnJNV+nt5+3pmFf2zSUMIom7bWv\neubb99C1TOa57ompbdy5c2cqddjbmKXbm2N/p9CYjkVqGlWSRBdSB0ljDWQjQa/ka5MLKgFNIGms\niV7vosbjf5P0KUnPaTR6T1tbTypNJxMSxyYMu7ZQljTcMRf2MpQwzOe69YlFruVq1zbu3LlzfF8q\nTxwxSSPUm1I0IrQLqWEeJI01kaUN++zXcxqN3iNxYOVIGmsiTdPguI1+/4+Pt7GXMSeUhUZTutPs\nuVP/u9Pl2Z/rnkdielnscRxlScPXWyNlqcIkjE8++eT4MXsbtzfFrm0UjeEoWsbUNGJGkXYBSWON\nJMm+er3HZY/b6Pf/tMldwgYiabSc26uSHaZMxm2MRm+r13s6H88RPxJUmk4a7tXUTT3AnTIvVMeQ\nZhOG+YX3XcSorA4S6kU5ODg4ThZmWaW2ETpfJWZkaNVE0dXEQaPRMkVdrxPT4zYGg3e1vf2UkmR6\nwJc736ZvXgvzj8adx7Ns3k37c0LFVPdKaPahTejwxH0v31B1t1FwGxG3sShqNNzPLRrktcjhR5ca\nEA5P1lCv94qkm9Yj5zQc/gtFUawESWPNZL9Yt9TrPa7x+Geyi6LD4de1vf3K8a+jSQfur6mdGtyk\nYX7RfcnC5jv5y00a7vByezBY6PDEPcTxpZfQMPKyLtfBYFA6cU9R4tj0AqhB0lhT/qLonzW5S9gQ\nJI014at1uEXR4fDKVG3DTRy+5FF25bZQbcU3dZ5bh/Bd/9W8Z6imESqI2mkhNPDLN2GPvT4cDmdO\nny87gW2emkbXEwdJY+1NF0WpbWDZSBotVdaLYp7v9V6ZqW0cHl7Tzs7TwcRh1zbcFGJ+0UOf5zut\nPlTLCJ0yH1PTCL233QPiJgm3huEbRh6qYcSeuOZ7bNOSB0ljzSXJvvr9J+T2poxG32hql9BxJI01\nE6pt9PtP6OjoX5Uljtsaj89qPL5vatBXKHHYyq5hav+6uj0dZQnDrmmEkoZb0/AljtAkO24C8U3C\nE+o9qXJqfNnfputIGh2RJPva2npS0m+UXTz68zo8vEZ9A7UjabRcbG0j2+aWsjNhv5A/ktU3dnef\nmbo6m3+iH/8vf2iEZNEEPu5JbWbpu+K8UfQ59rLoJDt3bIdJHPZrY2sZDCMPI2l0THbG63R9g8SB\nOpE01pRvur/MLae+IZmG48SJ31OS7Bdez7Rsohn7l9j8SpclC/dyBb7T6g23zuI7RyQ0ErXsco2+\npBFTyzDLOqb76wIajTUROkzJ5tlwB2S5hVFJOqeDg3e0s/O03Ku0xXyue9gwGo1mBm2FGg13GSN0\nuFDUaIQOaexllcFc9tJnU69/wuFJR00Ko9OHKoPBuxyqYCEkjTVTVBh1f+GSxH+oEpM4yubk6Pf7\n3jlHpXDCcAeSzfO59lXhfCnEXvpOe6+aMOhynUXS6DgSB+pG0ugAN31USRymOFp2zRSTFuwE4CYJ\ne3Ife3/qSBr2ellaKKpbxF5Jrah7dVMSRQhJY0OEEsfBwTskDlRC0lhTvtrGIonDHQBW9otrD9AK\n1S6KEkZoMFno8+1EEHsdVt/zVZNFTMLYtORB0tgwocRxePieDg+/R+pAKZLGmvON06ieOG5L+pTS\n9A91795V7e4+o17v9tRrfZMVuxMXu4mibBn6PmXLsiHuRfWKeRKGbz302CYgaWwokziS5LKkn1jP\nMOwcxUgaHVA0WtR+3J84XlSanvHUOX4p6W3t7r4y1bti96K40wfGJouYpOGuV0kLi1z5nYRRjqSB\nQJ1jV9np9dQ6MI2k0SHzJo7sOV/PyqTWYUaRmqvU+6YKjE0W8yQNe32RukRsoiBhhJE0cGySOP5O\n0vfl1joYRQqJpNFJiyaOra0/yp+frXUcHv6XpHeUJJ9oZ+dPlCT7CyWL2O/iW4/t8Sh6TdnnlT2+\niWg0Oqys8bCf8/+j8B2ynJT0rNJUOjy8ll9nZd/7OYs0HjGNxTzrMYchVZ7fRByeoNDsbOe3rGfP\naTj8pQaDf9Rw+H0OXTYESWMDhBKH77lQ4siuHftN9Xpf1Xj8Q/mSx2Dwrra2npwZGOaqUggtej72\nUGKeoiYJI4ykgShZ4nhRvd77hcljNLpO8ug4ksYGmR3cNXuym/tcWfJIkleVpj+SL3kMh/8t6ZqS\n5GP1+18+rn1U+RUv2rbOYibJIh5JA3OJSx6nlDUgX9Bo9CsNh2+SPjqARmODlf3y2rfQ49lzWfKQ\n/kZJ8oz8DciupOc0Gt3QcHhTw+GbGo/vm3kvM1GOufk+r+xW9h183xPxaDRQCzt5+BsQKRthuifp\nAUnP6ejofY1Gv9Zo9A8kkDVCTWPDFdU5irazt5197pZ6vYv5/ceUpq/l909Lei6/f1vZ5SP3lDUg\n1yVdk/SxkuTlqQmBqn6HuraFH40GliobMXpeUjbC1N+ASKb+kW33vtJ0vgYEy0ejgSmxycO3rct3\nEafZBuQPJH1aWf3jbL6l3YDcVJpekXRh7saDdFEvahpoRJLsq9c7ryR5VNJfS3paRQVU6UOl6f8o\nTa8rTa8oTf+eOggad0lSyi3uliRJ7TfpbCp9N5UeTqUPUinNb/vWffv2SSq9nUpv5q9t/u+yRjcT\n6yrj8ARzKYr8857Vmh1+XMzf/zFJr+fPnJL0fH5/X9IZ6/Fn8/u/lvQbSdeVFVlfyrdF3Wg00EpZ\nA/IlScoPQ0wDcknSjyV9RtN1kF1JD+U3SfpAWW8MDQiWh8OTtbmdSaXvpLOHMbes++5hzIep9NtU\n+kH++qa/Q+O3uQ9PKIRiDZnDmBuSHpX0PUlvSHpS0of5Nu5w9nPKBpU9ryyFvC3pTU0OdYDqSBqd\nuPlSyP/lS1KIdSNpABk3hXxX0lPK0kjopDo7hdyU9FtJPxApBGVIGp2/laUQt2u30ymEpAGUK0oh\nb0j6J2d7UgiKkTQ2/nYmld5Is3RRlkLuplkCuZFOBpitVRKZO2lggkaDW36zD2PeSKW/TbPDk7IR\nqgfpGh3GzN1ozH9Biu65JOlrTe8E2uqMpG8r+/d2Wv4RqsY9SXclXZV0QS0dWHa/pivC0RgRCkTZ\nl/TF/L7dgHxFkxGq0vQ8IaYOcje/RQxx35b0iLKSy5l8s1/kLx3W+oXmRqMBVGY3IFL2L9yXQuwG\n5AFNhrh/JOktzaSQ85IezF/St97+qfytPpB0uc7vMR8aDWBhoRQSOozZ1UwK2f536cFXpPs/mn37\nvrKDCSlLIg0nDhoNoFbzHMY8ID3yc+m+/y1+6/uUHbr8R+07XQmNBrA0sYcx+9Kjl6V+SYTo529B\nowFsioIUcuaRuLdowZgyGg2gEU4K2b8Z3HLmZQ1jGDnQBr8YSkcl2xwp635tGI0G0AZmCEeR2/l2\nDePwBGiDobJxGNLsOI0jTcZptGCAF40G0BaXxYhQABUNlXWpNtytWoSaBoBKaDQAVEKjAaASGg0A\nldBoAKiERgNAJTQaACqh0QBQCY0GgEpoNABUQqMBoBKuezJxQtmEjcAm2Jc0bnonAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAkv0/WzwTNDHnsDcAAAAASUVORK5CYII=\n", "text": ["<matplotlib.figure.Figure at 0x10a601f10>"]}], "input": ["clf\n", "imageplot(W) \n", "set_cmap('gray')\n", "h = plot(imag(gamma), real(gamma), '.b', linewidth=2)\n", "h = plot(x0[1], x0[0], '.r', markersize=20)\n", "h = plot(imag(x1), real(x1), '.g', markersize=20)"], "collapsed": false, "language": "python"}]}], "nbformat_minor": 0, "nbformat": 3}