{
  "metadata": {
    "name": "",
    "signature": "sha256:bcd8ac50f3861f1e3988a696e15e04380d827e0a32513c3fbf83eeaf0468bbbc"
  },
  "nbformat": 3,
  "nbformat_minor": 0,
  "worksheets": [
    {
      "cells": [
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "# Skeletonization using Proxies"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "import numpy as np\n",
            "import numpy.linalg as la\n",
            "import matplotlib.pyplot as pt\n",
            "\n",
            "import scipy.linalg.interpolative as sli\n",
            "\n",
            "eps = 1e-7"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 1
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "sources = np.random.rand(2, 200)\n",
            "targets = np.random.rand(2, 200) + 3\n",
            "\n",
            "pt.plot(sources[0], sources[1], \"go\")\n",
            "pt.plot(targets[0], targets[1], \"ro\")\n",
            "\n",
            "pt.xlim([-1, 5])\n",
            "pt.ylim([-1, 5])\n",
            "\n",
            "pt.gca().set_aspect(\"equal\")"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "display_data",
              "png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAEACAYAAABccqhmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX10VdWZ/z8PCXkxAURAo+hMJI7VivWl1Z9vC6/j4mWM\nr/VnFeuy45SO7TjA2P7Gzg9E43SYztjf6mjQVqdqp5ai1qli8Y6YIIZofakdQUVwqoH4AkYUUZIQ\nEpLs3x/PObnnntwLgXtzb+46z2etrHvuPvvssw/c/d3Pfvaz9xHnHIZhRJNR+a6AYRj5wwTAMCKM\nCYBhRBgTAMOIMCYAhhFhTAAMI8IUZ1qAiLQCO4E+YI9z7vRMyzQMIzdkLACAA2LOuU+zUJZhGDkk\nW0MAyVI5hmHkkGwIgANWicgfRORbWSjPMIwckY0hwNnOuQ9FZBLQKCJvOeeey0K5hmEMMxkLgHPu\nQ+/zYxF5HDgdGBAAEbHFBoaRJ5xzex2eZzQEEJGDRGSMd1wBzADeSFGJEfV366235r0OhVCnkVov\nq9PQ/oZCphbAYcDjIuKX9SvnXEOGZRqGkSMyEgDn3Gbg5CzVxTCMHBPJSMBYLJbvKgxiJNYJRma9\nrE7ZQ4Y6VjjgG4i44b6HYRiDERHccDoBDcMobEwADCPCmAAYRoQxATCMCGMCYBgRxgTAMCKMCYBh\nRBgTAMOIMCYAhhFhTAAMI8KYABhGhDEBMIwIYwJgGBHGBMAwIowJgGFEGBMAw4gwJgCGEWFMAAwj\nwpgAGEaEMQEwjAhjAmAYEcYEwDAijAmAYUQYEwDDiDAmAIYRYUwADCPCZEUARKRIRNaKyIpslGcY\nRm7IlgUwH9gA2EsADaOAyOj14AAiciRwAbAY+G7GNTKMIdAcj9NQX09xdze9paXMmDePabW1WS0r\nm/cYqWQsAMC/AX8PjM1CWYaxT5rjcZ6eP5/FLS0DaQu94301XP9cx9atvPfuu2zbtYuJfX2cDPQC\nRwB3P/ccvzz0UHZu28YNXV1MC9xj/SuvsPXFFynu7uaDnTspAQ4dO7ZwBcI5d8B/wIXA3d5xDFiR\nIo8zjGyycMYM52DQ380zZ7o1Tz7pFtTUJKUvqKlxa558MuW5OeDWeMdrwF0fKvN6cH8NbqF3/mvl\n5QN5F4Ty+vcZKXhtb69tOFML4CzgYhG5ACgDxorIg865a4OZ6urqBo5jsRixWCzD2xpRpri7e+C4\nGXgQ6AD6V63ijRde4Lvt7Un5F7e0sGjJEpxzSVYDwM+ARcA0oAG4GrgZNY17ve+NwA+AhcDBXV00\nA3cDj4Tq5d8nX1ZAU1MTTU1N+3VNRgLgnFsALAAQkXOB/xNu/JAsAIaRCc3xOOtefVWPgV8AVcB9\nAH190N4+4IiaFriuaPfutGUWeZ8fA0+jziyfhcAn3vFi4BIvz8EkC8UM7357u89wE+5cb7vttn1e\nkw0fQBCbBTCGDX/sP7G9nYXADuBw4J/882gvPhbtodcDW9GGvfV3v6O0rCxluX3e52fAvaFzi4Gr\nAt9LgMnA28CkQPov/LLS3GOkkjUBcM6tAdZkqzzDCNNQX8/ilhbqgEOBzWiDBG384d7728CX0F7p\n3t5emjs6WBjKcz3wde/48DT3rQocHwc8CZxCQnhALYV/KSnhH+bO3d/HyivZtgAMY9jwx/69aM8+\nBe2JQXv+xaH89wBXkhir+0OCRcB7wKfAX6A+hJ8AXWnuO8b7vBFoB8oBQUXHL3Mx8L9LSgpuFsBC\ngY2Cobe0FNDx9nveZwfa+6brycpD36ehDr2jgZOAqcBhwMPA97yyglwPbEMDXXaivobfoL3/06gI\n+FSMHr3/D5VnTACMgmHGvHksrKkBtOGv9tJnAhvTXLMdddY1h9L7UAH5V7T39v0Hn6BWw2zUUvg6\nMAEoBe73rm32yhyN+hr8siuPPvqAny1f2BDAKBim1day/pVXWHb77TzWpQZ7M7AMuAEGje8XoCb+\nemAJUI967/uBv/TSx5DafzAH7fkbgVnA/3jpqfJ+B/jnoiIuvuiiLD1p7hCNFxjGG4i44b6HEQ1+\nUlfHU4sX8+Xe3qSpt2Z0DD8O+Bw1a2tQp947wI8DZSxEhw9noo37cOAQkh16PovQ4cIc1DJYjvb8\n6fL219Qw8847R4wfQERwzsne8tgQwCgIflJXx+uLF7Oit5c6ksfg01Dv/L3oWL4YbbhbSW78oD13\nNfAUOpPQBmxKc88/okOBTcAXgG+S3mQuQgOBGpcsOYCnyx8mAEZBsOauu7intzcpbTHaizcDL6Nm\n/ZVob30z8EGasorQcf37wPHosCAVxwIPAXWoZfEN0vsa/FiCfAYCHQjmAzAKgvJQ4/fZhvoAngqk\nXY8ODx4kdbRen3fdQSQcgKniA07yjhvQKUWfVL6GWd5xZAOBDGM46SpO/VN9H/ivUNq9qFOwjMHB\nOv+B9uYXAs966cH4gCLU4fcd1LqA5Ebi570BjUSsQRv/NOCbInyjwAKBbAhgFATVM2bwrVDaHNSU\nT8UO4K5Q2mJgC7p7zd+QHCPgxwfUAcd439/zvm8keRpxGjr9d5D3fTUqHm0VFSPGAThUzAIwCoLR\n27dzChqXXwbsRhviU2nyf4423qDpDzpT4B/PQXewCToKg+b8n3hlAAPi41/718C1ge/fAGq/9739\ne6gRgAmAURAUd3fzN2jP7dMMxFFz/aeB9DnA90k0Tj+6bxqwC/ULBEXhBjQsuAJdW7Aa7eHPDZT5\nMzQa8DHUMuj28vwKnU1oHzWKqaedlulj5hwTAKMg8MOA/Yi9j9F4/LiXtgh4F50B+AeSlwIv9s6v\nBG7yzn0HdRJeC3SNGkVLfz+novEEPgu971vRhiLoWoDHQnmmAz/o78/rXgAHivkAjIJgxrx5XHzQ\nQSxDHXuTSHjm/fH7g+gPelqK698j4awDtRgOR5fxTrrySv6ksjLJ0w8qHGu8+9WhYlNFsj/An4qE\nwpsCBBMAo4Ao6e0daKTpTNd0k3A9qOUQbLxFaHz/uw0NVKSZZTg+9H0xKh7B9QX+hiKFNgUINgQw\nCoSG+nqm9vQMfE8dFaA/6CvRhus7AFeiJv801GT3Nwp5H23I7Z9+Sk9lZcry+lKkfQG1CBYG8txY\nVcVlBTYFCCYARoFQ3N2d1OhnMDgg52vAkSR79f1NQXzTfyYaOBQ097/tHF9qb+f6oiLu7Us0+eu8\nvzB+jsXAZehUZKGudjEBMAqC3tLSpEbvN+jLULP/GNTMD8f+34M6AH3CUX1+niuBc/v6uGrCBI6b\nOpW+sjIOPeQQfv7QQ0k+heA0Id697wNoaytIJ6AJgFEQzJg3j6dbWpjZ0jIQsbcOtQS2op74d9Jc\nWxQ4TveDPx4NEhpXVUVdYGfdy19+mcs2bWIsGhcQdCQC7AnepwCdgCYARkHg96yNS5ZQtHs3G9ev\nZ8b27XyILgT6CA3LTUVwHP/aXvIsBq5qa0tKP/Goo/hg0yY6gA9Jbvx/hU4LXonGDBSiE9D2AzAK\nkuunTmXSm2/yMTrF97/Q1X9jSR4G+Pv4HYk28h406OdngTy+WT8NuLqigmO/8hW27dxJD9D5xz+y\nq7MTP8av0SvvQzSAyBeEbwGHzJ7Nvy5blv2HPUCGsh+AWQBGwdEcjyObNjED3aQjGA78TbRhTgLW\novv8+Y3Uf6HHwUAtcBoqCkGzvqazk7o1urn1QtRpiHfd4UAlKiLhl4L8DLiqoSE7D5hDLA7AKDga\n6uu5p6uLBgY7/e5Hd/ipQx10wcb/NNpwv44G9OxBA4j8PAtQX4LPYi//cu/zDjQo6PM09ZKudPsK\nj1zMAjAKDn978GISocHBNf9FqCXQ733eT/K24Q1emh9CHFwCHI4i3IFOGxK6V3g9AUBHIE6hUDAB\nMAqK5nicjevXcz36ToCtwL8Hzn8XNf0PQX/cE1EnXX8gj/+jn0aiAd9MYn9Bv5Fv88qv8467SewM\nDMmLjBYAYydPzvj5co0NAYyCoTkeZ/mcOdywfTsT0TfT/nsoz4+Bc9C9+x9BzfVzgY5RiZ96qijC\nGWgg0dMkYv9/ApwN/Ll3nGodwE9RgTkSmHLccRk9Xz4wATAKhocXLeLHbW0D5vzeNuj0uQdttKOL\nigZ6bD+gCBJ7/N+H9vDhtwsFF/sEj32+gApNAzDdQoENY/jo3LwZSPxo060HCMfvlwFuzx4+AS5F\nd/LxLYMa4AEvX12a8orSHAfvVTR6dMFFAYJZAEYB0S06pe03/GBP7hP25IPuHXAD6v0/AXXqxb3j\nBwL5hiIowePgvfpEaI7H9/kMI42MBEBEykTkZRFZJyIbROSH2aqYYYTZM3483ybR8Keh8/SL0I09\nZqF7+Qc989ehe/81oPsFBE38zlD5+xKUOajj8Vrvnn78wI3Ad3t6Cu6dAJDhEMA5t1tEznPO7RKR\nYuB5ETnHOfd8lupnGIA6AP901y4uRcfhnwAXoe/sOwxtnNOAy730SnRM/3ckBOE7JL/Rtzt0Dz/d\nDxL62MtzH+pcnIFuSdbs1eEB7/My79rVUVwL4Jzb5R2WoEOkTzMt0zDCNNTX82MvTj/Yw9+ARun5\nnAi8ji7PDW7dBeqxXxS4vpLBS4p/RfJ+gsH7+PsR+tOH/qvDfApxLUDGAiAio4BXUX/KT51zGzKu\nlWGE8IN/wkwKfe9DI/y+lKacdwPH16Jbh/sbiPjj+1Rbiu0IfZ/jXe+zoKaGWVGcBXDO9QMni8g4\n4GkRiTnnmoJ56urqBo5jsRixWCzT2xoRw98UNEzYKdeGCkA6h57vEJzkXRtDdwGu887fnOa68SSi\nBjeWl1N96aU0fvopq3fvpq+sjFlz5+Z9FqCpqYmmwFLmoZDV1YAisgjocs79v0CarQY0MqY5Hufp\n+fNZ3NIykDanpISOnh5K0QU6FWivfH9NDcWdnVS1tSW/8rukhFN7ethCstlfi84KQJpXhY8ezYel\npfxZdTVjJk9m+gho7ENhKKsBMxIAEZkI9DrnPhORcvTf7jbn3DOBPCYARlZojscH9gPoKyvj8DPO\nYMvSpUmisKCmhll33gnAg4sW0dnaSglQefTRnHDRRWxZupSZLS00ktyb73n22QEfQzNQX1LCkcce\nW1ANPkwuBOBEdGflUd7fL51zPwrlMQEwho2wKOyrsabLv7/lFALDLgBDrIQJgGHkgaEIgEUCGkaE\nMQEwjAhjAmAYEcYEwDAijAmAYUQYEwDDiDAmAIYRYUwADCPCmAAYRoQxATCMCGMCYBgRxgTAMCKM\nCYBhRBgTAMOIMCYAhhFhTAAMI8KYABhGhDEBMIwIYwJgGBHGBMAwIowJgGFEGBMAw4gwJgCGEWFM\nAAwjwpgAGEaEMQEwjAhjAmAYEcYEwDAiTEYCICJHicizIvKmiKwXkXnZqphhGMNPpq8HrwKqnHPr\nRKQS+G/gUufcxkAeezuwYeSBYX87sHOuzTm3zjvuADYCR2RSpmEYuSNrPgARqQZOAV7OVpmGYQwv\nxdkoxDP//xOY71kCSdTV1Q0cx2IxYrFYNm5rGEaApqYmmpqa9uuajHwAACIyGngSeMo5d0eK8+YD\nMIw8MBQfQKZOQAF+AWx3zt2YJo8JgGHkgVwIwDlAM/A64Bf0f51zKwN5TAAMIw8MuwAMsRImAIaR\nB4Z9GtAwjMLGBMAwIowJgGFEGBMAw4gwJgCGEWFMAAwjwpgAGEaEMQEwjAhjAmAYEcYEwDAijAmA\nYUQYEwDDiDAmAIYRYUwADCPCmAAYRoQxATCMCGMCYBgRxgTAMCKMCYBhRBgTAMOIMCYAhhFhTAAM\nI8KYABhGhDEBMIwIYwJgGBHGBMAwIowJgGFEGBMAw4gwGQuAiDwgIh+JyBvZqJBhGLkjGxbAz4FZ\nWSjHMIwcU5xpAc6550SkOvOqDD/xxjj1y+rpdt2USinzrp5H7fTaQXkWLVnE5m2bkT6h+tBqLj73\nYl7c+OJer9ufewTzbfloC20ft3H44YdzxIQj9lq2YWSbjAWgUIg3xpl/93xaTmkZSGu5W4/9Bhdv\njDPnR3Nom9wGnwGjYMdbO1i7dS1cyqDrgKTGfubxZ7L0haV7vceg+7QBlbB953bWf7ae1//pde7j\nPhMBIyeIcy7zQtQCWOGcOzHFOXfrrbcOfI/FYsRisYzvub/MvG4mDdUN0Aq0oIOfTqjsquTLX/ky\npVLKx9s/Zu2EtXr+fO/CZwLHAcp/XU5fZR89F/Qk0h4vp+uyrkF5T331VC46+yLuevQueot66fis\ng74z++A1oDJU/jNwyphTeDX+anYe3IgMTU1NNDU1DXy/7bbbcM7J3q7JiQBk4x6ZEvvLGGtkTXLj\nBngCKAIqQD4SXLmDSwLnnwXOS1Hgw8AkVEj6gRpgc5q8v/byBPN3AF3AVYOzj39yPJ++8ul+PZ9h\nhBGRfQpAZIYApVI6uPG3ktQDOxw8CfweaEcb6ydpCnQM6rkZ3PkrfcCh+5G/KE26YWSZbEwDPgS8\nABwrIu+LyHWZVyv7zLt6HmW7ypITw4IAcCHwrpd+HnAasCKUZwVwZijtfLShLw+lPwF0p7jP+UBv\n6roefdjRqU8YRpbJxizA7GxUZLipnV7L8UuOZy1rE4np5G9i4Lja+1wKHAQcjPb+fnorST4FdgFx\n4HO0Jx8FlHj5/Gt8ioCVJE2iVj1fxT/e9I9DfCrDyIxIRQL+YO4PqFlbk0joT5Mx7LKoBiZ4x38O\nVHjHrSSsiPOAq9Fx/jbgMOBK4AovvcXLH6QUOA5YDRUrKpj57kzuu8lmAIzcESkBqJ1eyzVnXcOE\n+ATGPTWO8m3lFK8IGUHLgSkpLg6KxRjgUcB31LcGzl0CjCa1yb8p8H0l6gh8CZ1y3ANzZ8+1xm/k\nlMg4AUHn35e+sJTttdsH0kb/ZjTlvy6nZ3QPfeP6tFG2kGyur0LH60XomL4S7dl9nvE+/WtK0lTg\nc3RW4RPgT0kSik46mfOjORYDYOSUSFkA9cvqk4J0APZcvocu6aKvpw+2oo2yBliNNtY4OrbfA+xG\nRWBfvXsPqRmHikA5CUdjgLaz21jy0JL9fi7DOFAiZQF0u+7UJyaSmL9/FO3Jq0ke44P29Okk059t\nXYWKSBwIduTLgZ3ouL8UFZTfov6EGgash939u4f4NIaROZESgFIpTX0i6PTzp/0uYvA0YQ3wSprC\nP0TF4zRUDGpQEehBLYd2NBYgEFLMM14+3yiphrJRoalKwxhGIjUEmHf1PMpXlicnriLZ6VeNCsKj\nIDtCQVTV6Nj9ieRkeULgLLTxt5Bo1CcAY4EqdArx0uTrOB94k4EhRNXzVcydPfeAns0wDoSshALv\n9QYjJBTYZ+oFU3nzkzc1OKcPbbTVoUyrYczOMdAP7V9tH1zIE+hMwOcwRsbQ2dlJ/9e9aYJW1B/w\nAToUqEDvNQq4PEWFHgNOBfm9sOKOFeYANLLGUEKBI2UBAEw+bLL2zL6n/iWSp/FWAh9Ae0877UXt\nOk4Pshx1BH4G7IL2He30S7+KQisqJlOA8Wi0YAU6Y3BwmgodDGyCIldkjd/IOZHyAQBMKJkAz6Hj\n8XLUXH8Nil8oppde/Re51svcCqwB/gMVDF9LzyRhNQTH8a95aS2oDyG4krAGXWdwYaAyq4BjgFeh\ndFQa/4RhDCOREoB4Y5zlry/XCD2fZ4CToPzVctrHtScabCvakM9lsDMwOO/vOwYnog6/56GorIg+\n+pLtq2pUIFajQuLQxl+t1xx26GHZekzDGDKRGgLUL6una1ZoCZ7ngGvva9cQ3mdInv5LtWDIn/f3\n812BTiNeAYyFvo4+zRcONT4JFYnz0JDianTIcQ7sqthFvDGe8TMaxv4QKQHY8tEWbeDPkmjooD3y\nBDR+32/0fjBPun+hDrTnHxUq62J0PYA/NHgmdN1OEkFGq1EHIRoEdMtdtxzgkxnGgRGZIUC8Mc7b\n29+GCwKJfuN8Dx3jjybRcF/0zqVaMNSKika6cOAeNHT4TXTN/1JUYDqAr6Uob7Vet2HLBuKNcXMG\nGjkjMhZA/bL6pO27AO3tm4AydK5+PAmH3k40sKeTwfsBvIQ6+cJl+eHA49BFQWUk1hD4I4/WFJXz\nnIu7K3ZbKLCRUyJjAaQMA25FG31wCzDfAmhDHXv96Jz/w2gwTzfawFMh6Jj+OO97Ddrrh8uH5NgD\nx8CMgIUCG7kkMhZAyjDgFpIbJyR68iNQZ935aJhvJRo4BOm38vqIgTH9Psv3eQJdhOTNCFgosJFL\nIiMA866eR9XvqpITd6TJ7E/TQWLfwAuBy9Bxfy/a0wdZhcYVXEKigadbFfg58Bt0rQBouHA1lK8s\nt1BgI6dEZghQO72W+7iPW+66hXXvrKO/r19N+1R8iA4N/I07Lwyd/xrac4fn9Dd75wUVjs/SlN8N\nnEFiGLAMypeXc9Psm8wBaOSUyK0FAKg8qZLOr3YOXu4L2rBPItE4HwO+mqKQVNuFr0bn9x9Fl/ue\nlaL8FcCJJPkAKn5bwSP/8og1fiOr2Lbgaeju9wbq1V6C35NvAc4h2UFXRGIfAH///2o0aCiIH9a7\nCl1g9HKK8h0pNxQtGV1C/bJ6ABMBI6dEUgDKisrooEO/VJNokI+Q3Phb0c07wmHA/42a8Y+iwtCN\njve70YVGkPAhBMvHu8Yv+y1gFuxgBw008PqP7LVgRm6JjBPQJ94YVydeOEJvFckefFDzPfze4/PR\nDT4mow7By9FdfyeScCq+go7xw/dYju4nsBqNJQiVbdGARq6JlAVQd3sdty+7na7SLthOsmnuv4nL\n3w0I0stjMalfKLIaSjaU0L+nn95q760fwXt0oA7G4MrCEJs/2pz6hGEMA5ERgHhjnNsfvz3x8s5W\ntKcuRbfrOtvL+BJqphez91d9pULQaMNfemUUoY19KjoMiKMvDtlB+p2D05VtGMNAZARg0ErAau/z\n92hj3Iw6+aagPXU3Ol//ODr/77MK9fCnwh/3V5A8c7ASWOdd109il6Dwm4dXwfjy8fv1XIaRCZER\ngJShwC+hjTEY1/8EiV74KHRoEJ7vP8bLdwmJqcTPSbwCLDxzcBy6MKifxAKiVlR8QsOQcdXp4owN\nI/tERgCSQoFbgfVowzuIxAKd19DGGwzffQyd2w/zAmrml4TyL0cFJDxz0EfytmDV3qe/mUgbcAyM\nnTR2iE9kGJmTjbcDzxKRt0TkbRH5fjYqNRzMu3qevhewFe2xL0R3BjofeB143ssYjvpL1yGPBg5h\ncKz/pWgUYZDzUT/DFPT+/p4ELagP4jzvPqfbWgAjt2QkACJSBNyFTmh9EZgtIsdno2LZpnZ6LXfe\ncCcT3pww2IN/Mbr4p4jBy3X9vfyCPIH2/Oli/VN5+H1bK/gy0fO9vMuBE2wtgJF7MrUATgfecc61\nOuf2oItmw33iiKF2ei1Tvzg19UlBe39/IU8r2lNvRntv37O/Gm345ei4PxWpIp8rUXM/LD4XobMN\nm+DQ4kMtCMjIKZkKwGTg/cD3D7y0Ecs+3w70GYNf+z0b3earFPUHXIG+4DNdsM8nobRVaIRgZZpK\nVWm5Hd0dQ34Ow8gGmQrAyFrlMwQGfAFBgm8HKiV1Tz2L5EjBchK7Agf3+BuF2kV+2qMkdv9NN/fv\nbQhSNbEqTQbDGB4ynQXYgk6W+RyFWgFJ1NXVDRzHYjFisViGtz1wfBP7lrtuYcOWDeyu2D3QQMtX\nljNl4hRat7XSSefgiwM9eHGP9x6BapJj/b39/QZ2/C0JnK+B4v8qpveC3kT+Fej04ElwpByZ+QMa\nkaWpqYmmpqb9uiaj5cAiUgz8D9pfbkVntmc75zYG8oy45cA+8cY4Sx5awu7+3ZSNKmPu7LnUTq9l\n5nUzaahuGHyBt9y35tUarjn7Gpa+sDTpdePjnxlP58ed9BzWo736FKhqqeKI8UcwZvwYykaVccZx\nZ7DidysS4jMFqNYy7/zbO80HYGSNoSwHzng/ABH5C+AO1Id+v3Puh6HzI1YA0hFvjDP/7vlJjbt8\nZTlTxk/hyKojB4QilYAAKUUl1T2Gks8wDpScCMAQKlFwAgDWQI3CxwTAMCKMvR3YMIy9YgJgGBHG\nBMAwIowJgGFEGBMAw4gwJgCGEWFMAAwjwpgAGEaEMQEwjAhjAmAYEcYEwDAijAmAYUQYEwDDiDAm\nAIYRYUwADCPCmAAYRoQxATCMCGMCYBgRxgTAMCKMCYBhRBgTAMOIMCYAhhFhTAAMI8KYABhGhDEB\nMIwIYwJgGBHGBMAwIswBC4CIXCEib4pIn4icms1KGYaRGzKxAN4ALgOas1SXnNHU1JTvKgxiJNYJ\nRma9rE7Z44AFwDn3lnPuj9msTK4Yif9ZI7FOMDLrZXXKHuYDMIwIU7y3kyLSCFSlOLXAObdieKpk\nGEauEOdcZgWIPAt8zzn3aprzmd3AMIwDxjknezu/VwtgP0h7k31VwDCM/JHJNOBlIvI+cAYQF5Gn\nslctwzByQcZDAMMwCpeczAKMpKAhEZklIm+JyNsi8v181sWrzwMi8pGIvJHvuviIyFEi8qz3f7Ze\nROaNgDqVicjLIrJORDaIyA/zXScfESkSkbUiMmIc4yLSKiKve/X6fbp8uZoGHBFBQyJSBNwFzAK+\nCMwWkePzWSfg5159RhJ7gBudcyegQ7wb8v3v5JzbDZznnDsZ+BJwnoick886BZgPbABGkjntgJhz\n7hTn3OnpMuVEAEZQ0NDpwDvOuVbn3B7gYeCSfFbIOfccsCOfdQjjnGtzzq3zjjuAjcAR+a0VOOd2\neYclQBHwaR6rA4CIHAlcANzHXpzheWKf9YlaINBk4P3A9w+8NCMNIlINnAK8nN+agIiMEpF1wEfA\ns865DfmuE/BvwN8D/fmuSAgHrBKRP4jIt9JlytY0YKEEDY0kE23EIyKVwH8C8z1LIK845/qBk0Vk\nHPC0iMScc035qo+IXAhsc86tFZFYvuqRhrOdcx+KyCSgUUTe8qzNJLImAM656dkqaxjZAhwV+H4U\nagUYIURkNPAbYKlzbnm+6xPEOfe5iMSBrwBNeazKWcDFInIBUAaMFZEHnXPX5rFOADjnPvQ+PxaR\nx9Hh7yCXvOtWAAAA3klEQVQByMcQIJ/jpD8AfyYi1SJSAlwJ/DaP9RmRiIgA9wMbnHN35Ls+ACIy\nUUQO9o7LgenA2nzWyTm3wDl3lHPuaOAqYPVIaPwicpCIjPGOK4AZqCN+ELmaBhwRQUPOuV7gb4Gn\nUa/tI865jfmoi4+IPAS8ABwrIu+LyHX5rI/H2cA1qKd9rfeX75mKw4HVng/gZWCFc+6ZPNcpzEgZ\nYh4GPBf4t3rSOdeQKqMFAhlGhInaLIBhGAFMAAwjwpgAGEaEMQEwjAhjAmAYEcYEwDAijAmAYUQY\nEwDDiDD/H5h6FLMi52vMAAAAAElFTkSuQmCC\n",
              "text": [
                "<matplotlib.figure.Figure at 0x7fd7ac83d160>"
              ]
            }
          ],
          "prompt_number": 2
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "def interaction_mat(t, s):\n",
            "    all_distvecs = s.reshape(2, 1, -1) - t.reshape(2, -1, 1)\n",
            "    dists = np.sqrt(np.sum(all_distvecs**2, axis=0))\n",
            "    return np.log(dists)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 3
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "def numerical_rank(A, eps):\n",
            "    _, sigma, _ = la.svd(A)\n",
            "    return np.sum(sigma >= eps)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 4
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Check the interaction rank:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "numerical_rank(interaction_mat(targets, sources), eps)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 5,
              "text": [
                "9"
              ]
            }
          ],
          "prompt_number": 5
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "**Idea:**\n",
            "\n",
            "* Don't want to build whole matrix to find the few rows/columns that actually matter.\n",
            "* Introduces \"proxies\" that stand in for\n",
            "    * *all sources outside the targets* or\n",
            "    * *all targets outside these sources*"
          ]
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "## Target Skeletonization"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "nproxies = 25\n",
            "\n",
            "angles = np.linspace(0, 2*np.pi, nproxies)\n",
            "target_proxies = 3.5 + 1.5 * np.array([np.cos(angles), np.sin(angles)])"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 38
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "pt.plot(sources[0], sources[1], \"go\")\n",
            "pt.plot(targets[0], targets[1], \"ro\")\n",
            "pt.plot(target_proxies[0], target_proxies[1], \"bo\")\n",
            "\n",
            "pt.xlim([-1, 5])\n",
            "pt.ylim([-1, 5])\n",
            "\n",
            "pt.gca().set_aspect(\"equal\")"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "display_data",
              "png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAEACAYAAABccqhmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXt8HNWV57/HkvVAtsHYgHntCCtDIHGGRwILgY9phrXk\nIJ5hA5jkQ4aJEyYQy5tkh+z4gZUhmuyQ/WRiyTBMAsmEeA2ECTGxe7AlY2RBIIQsJsQY8hCIgMEY\njMGybEuWdPePU6WuLnXr1a3uFnW+n48+XV1169atLt3fvffcc26Jcw7DMKLJpHwXwDCM/GECYBgR\nxgTAMCKMCYBhRBgTAMOIMCYAhhFhijPNQEQ6gL1AH3DIOXd2pnkahpEbMhYAwAEx59y7WcjLMIwc\nkq0hgGQpH8Mwckg2BMABm0TkNyLyxSzkZxhGjsjGEOA859ybInIU0CIiLznnHs9CvoZhjDMZC4Bz\n7k3v820R+TlwNjAgACJiwQaGkSecc0MOzzMaAojIYSIy1duuAKqB36UoREH9rVixIu9lmAhlKtRy\nrVixgurqpejoM/mvpmbZkOeO9bzhzi3E32kkZGoDOAZ4XESeA54G1jvnmjPM0zCGpbs7def14MGi\nIc+rq6umqmpp0r6qqiUsWjRv3K5ZyGQ0BHDOvQKcnqWyGMaIKS3tTbm/rKxvyPNqa+cC0NS0nIMH\niygr62PRovkD+8d+zYkpAtkwAk44YrFYvoswiEIsE+SmXPF4G42NzXR3F1Na2ktdXfWQFTIWi3HW\nWZNob19Ke3vDwH5tyecPe73a2rkjqvBh6uqq016zoqI/7Xmjvb9cIiMdK4z5AiJuvK9hTFzi8TYW\nL94YqlRLWbmyZthKEo+30dTUEmjJ5417xRrtNTO5v0wREdwwRsBcGCKcYaSjunqpAzfor6ZmWb6L\nlhXyeX9e3RuyflowkJFXPoiGtSCFfn8mAEZeGasxb6JQ6PdnAmDklUym5SYChX5/ZgQ08k4+jHm5\nJF/3NxIjoAmAYXxAGYkA2BDAMCJMJB2BjPGnkJ1fCp1c/nYmAEbWSeX80t6uhjATgaHJ9W9nQwAj\n6zQ2Nif9AwO0tzfQ1NSSpxJNHHL925kAGFmn0J1fCplc/3YmAEbWKXTnl0Im17+dCYCRdQrd+aWQ\nyfVvZ34AxrjwQXfuGU+y9duZI5BhRBhzBDIMY0hMAAwjwpgAGEaEMQEwjAhjrsDGqDE///yTrWdg\nAmCMCvPzzz/ZfAY2BDBGhfn5559sPgMTAGNUmJ9//snmMzABMEaF+fnnn2w+AxMAY1SYn3/+yeYz\nMFdgY9SYn3/+GckzyFksgIgUAb8BXnfOXRo6ZgJgGHkgl7EAi4Ht6AvTDcOYIGTsByAiJwAXAw3A\n1zIukWGMgLZ4nObGRoq7u+ktLaW6ro65tbVZzSub1yhUsuEI9C/A3wPTspCXYQxLWzzOxsWLaWhv\nH9i31NseruL6x/a98QZ/fvVVdu3fz8y+Pk4HeoHjgDsef5yfHH00e3ft4uYDB5gbuMa2Z57hjaee\nori7m9f37qUEOHratIkrEMO9PXSoP+AS4A5vOwasS5Emmy88NQy3tLp68Ot2wS2rqXFb1q93S6qq\nkvYvqapyW9avT3lsIbgt3vYWcDeG8rwR3JfALfWOX11ePpB2SSitf51CgRG8HTjTHsAngctE5GKg\nDJgmIvc6564PJqqvrx/YjsVixGKxDC9rRJni7u6B7TbgXmAf0L9pE7978km+1tmZlL6hvZ3lTU04\n55J6DQA/AJYDc4Fm4DpgGdo17vW+twC3AUuBIw4coA24A3ggVC7/OvnqBbS2ttLa2jqqczISAOfc\nEmAJgIhcAPzPcOWHZAEwJg6FGPTTFo/z3LPP6jbwY2AWcDdAXx90dg4YooIlLTp4MG2evv/c28BG\n1JjlsxR4x9tuAC730hxBslBUe9cb6jrjTVfXJJ56qnfgeY2EbAcD2SzAB4RCDPrxx/4zOztZCuwB\njgW+5R9HW/FpaAu9DXgDrdhv/PKXlJaVpczX9597D/i30LEG4NrA9xLgeOCPwFGB/T/280pzjfEm\n1fNKlrI0DDdGyPQPswFMSKqrl6YaZruammV5K5M/9l8B7g5wC8AtDozfw2PyG710S4ZI86WADWBx\nqhsO7V8G7lOpxv/gPlVSkjcbQOrnNf42AOMDSiEG/fhj/160ZZ+NtsSgLX+4vbsLuIbEWN3vtywH\n/gy8C3wKtSHcCRxIc92p3udXgU6gHBC0x+Hn2QD895KSvI3/0z2v4bBYACMlhRj001taCuh4+8/e\n5z50nJ7u37889H0uatA7CTgNmAMcA9wPfN3LK8iNwC7U0WUvamv4GTrs2IiKgE/F5Mmjv6ksMdIx\nfxgTACMlhRj0U11Xx9KqKkAr/mZvfw3wYppzdqPGurbQ/j5UQP4Zbb19+8E7aK9hAdpT+CwwAygF\n7vHObfPynIzaGvy8p5x00pjvLVNSPa+RYEMAIyW+oa+paXkg4GR+XmcB5tbWsu2ZZ1hz++08dEA7\n7G3AGuBmtPUODgOWoF38bUAT0Iha7/uBv/H2T/XyCFv/F6ItfwswH/i9tz9V2i8D/1RUxGWXJoXB\n5JRUz2vjxuHPs2hAY8JwZ309jzQ08PHe3qSptzZ0DH848D7ara1CZwj+BHw3kMdSdPhwLlq5jwWO\nJDGTEGQ5OlxYiPYM1qItf7q0/VVV1KxcWTDegPZiEOMDw5319Tzf0MC63l7qSR6DzwVOQafw7ke7\ntbehhsLvhvJpACqBR4CjgZ3Ay2mu+Qd0KPAy8GHgC6TvMhehjkAtTU1juLv8YQJgTAi2rFrFXb3J\nhq4GtBVvA55Gu/XXoK31MuD1NHkVoeP614BT0WFBKk4G7gPq0Z7F50lva/BNo/l0BBoLZgMwJgTl\nvamt3LtQG8AjgX03osODe0ntrdfnnXcYCQNg2H5wIzpLAGocvCtwLJWtYb63nS9HoLFiAmBMCA4U\np/5XfQ34z9C+f0ONgmUkj9eXAv+OtuaXAI95+4P+AUWowe/LaO8CkiuJn/Zm1BOxCq38c4EviPD5\nRYtGflMFgA0BjAlBZXU1XwztW4h25VOxB1gV2tcA7EBXr7mJZB8B3z+gHviQ9/3P3vcXSZ5GnItO\n/x3mfd+MisfOioqCMQCOFOsBGBOCybt3cwbql18GHEQr4iNp0r+PVt5g1x90psDfXoiuYBM0FAa7\n8//FywMYEB//3C8B1we+fx6o/frXR3dTBYAJgJGSQosELO7u5ia05fZpA+Jod/1fA/sXAt8gUTl9\n95i5wH7ULhAUhZtRt+AKNNhnM9rCXxDI8weoN+BDaM+g20vzf9HZhM5Jk5hz1lmZ3mZWqK+/k1Wr\ntoworQmAMYhCjAT03YB9j723UX/8uLdvOfAqOgPwv0gOBW7wjm8AbvGOfRk1El4PHJg0ifb+fs5E\n/Ql8lnrf30AriqCxAA+F0swDbuvvz+taAD719XfS0PA8vb0PAD8d/oThooUy/cOiAScchRgJuGX9\nenfpYYcNrNizNE3kXm2a/dcHov6CkX1/C+6WBQvcVVOmpDzv6hRRf6nyceBWXHBB3n4fnxkzrh5V\nNKAZAY1BFGIkIEBJb+/AdFy6rmu6SbgetOcQNOYVof79rzY3U5FmluHU0PcGdLgRjC/wf5VCmALs\n7Q2HPw2NDQGMQRRiJGBzYyNzenoGvqeLfStGnYFOJWEA3IB2+eeiXXZ/oZDX0Irc+e679EyZkjK/\nVHf8YdQ4uDSQ5quzZnFlAUwBFhenC2pOjfUAjEEUYiRgcXd3UqWvZnDo7tXoqr4PwIC78BrgBBI2\ngRrgee/Yj73PE53jv3V2cmNRcg/nBnR8H8YXhQZ0Sew30VDhQuArX7mA4uK/G3F66wEYgyjESMDe\n0tKBSt9AokJfiXb7P4R288O+/3ehBkCfsFefn+Ya4IK+Pq6dMYNT5syhr6yMo488kh/dd1+SQTE4\nTYh37bsBdu4sECPgTcCdrFp1Lbt3D5/eBMBISW3t3LwvABqkuq6Oje3t1LS3D3jsPYf2BN5AW+o/\npTk32K6n+4c/FXUSOnzWLOoDK+te9fTTXPnyy0xD/QJ8rz+fQ8HrFEgcQH39TdTX34RIeN3iwZgA\nGBMCv2VtaWqi6OBBXty2jerdu3kTDQR6C3XLTUVwHP/bIdI0ANfu3Jm0/2MnnsjrL7/MPrSrH6z8\nf4tOC16D+gwUghFwtNh6AMaE5MY5czjqhRd4G3XM+a9o9N80kocB/jp+J6CVvAd1+vlBII3frZ8L\nXFdRwcmf+AS79u6lB+j6wx/Y39WF7+PX4uX3JupA5AvCF4EjFyzgn9esyf7NjpGRrAdgPQBjwtEW\njyMvv0w1ukhH0B34C2jFPArYiq7z51dS/4UeRwC1wFmoKAS79VVdXdRvUS+6pajREO+8Y4EpqIiE\nO9c/AK5tbs7ODeYQmwUwJhzNjY3cdeAAzQw2+t2DrvBTjxrogpV/I1pxP4u+TOQQGgDkp1lCstW/\nwUu/1vv8Hjpr8H6acsmB0U3BFQLWAzAmHP7y4MUkXIODMf9FaE+g3/u8h+Rlw5u9fb4LcTAEOGz2\n3INOJRK6VjieAGBfwE9homACYEwo2uJxXty2jRvRdwK8AXw/cPxraNf/SPSfeyZqpOsPpPH/6eeS\nqMDLSKwv6FfyXV7+9d52N4mVgSE5yGgJMO344zO+v1xjQwBjwtAWj7N24UJu3r2bmeibab8fSvNd\n4Hx07f4H0O76BcC+SYl/9VRehNWoI9FGtJtfjwYCnQf8tbc9i2RXYt8t+BrUyDj7lFMyur98YAJg\npCUeb6OmZhmxWD01NcuIx8Or6+eW+5cv57s7dw5054daoNPnLrTSTi4qGmixg16E/hr/d6MtfPjt\nQv66g+Ftnw+jQtMMzMuzK3D4eY0EGwIYKSnEkOCuV14BEv+06eIBwv77ZYA7dIh3gCvQlXz8nkEV\n8EMvXX2a/IrSbAevVTR5cl69AMf6clDrARgpaWxsDv0zQXt7A01N4TYwd3SLTmn7FT9VPEDYkg+6\ndsDNqPX/o6hRL+5t/zCQbiSCEtwOXqtPhLZ4fNh7GC9SPa+RkJEAiEiZiDwtIs+JyHYR+XYm+RmF\nQyGGBB+aPp2/I1Hx56Lz9MvRhT3mo2v5B/snN6Br/zWjC4AEq0hXKP/hBGUhani83rum7z/wVeBr\nPT15fSfAWF8OmtEQwDl3UEQudM7tF5Fi4AkROd8590Qm+Rr5p9BCgtvicf5i/36uQMfh7wCXou/s\nOwatnHOBq7z9U9Ax/f8gIQhfJvmNvt2ha/j7fSeht700d6PGxWp0SbI2rww/9D6v9M7dnMdYgLG+\nHDSbK/8cBjwDfCS0f1xXQDHGh/Xrt7iqqiVJC+RUVf2DW79+S17Ks7S6OuWKPTeFvq8AdyW4a9Ks\nDLQssP0Fb4Wf4PEvpVjxJ9V1wnk5cMtqavLy2ziX+nkxghWBMjYCisgk4FnUnvKvzrntmeZp5J9C\nCwn2nX/CHBX63od6+P1VmnxeDWxfjy4d7i8g4vdtUt3hntD3hd75Pkuqqpifx1mAsb4cNGMBcM71\nA6eLyOHARhGJOedag2nq6+sHtmOxGLFYLNPLGjmgkEKC/UVBw4SNcjtRAUjXIfYNgkd558bQVYDr\nvePpJs+mk/AafLG8nMorrqDl3XfZfPAgfWVlzF+0KO9rAVRU9HPOOb6Npij3bwcWkeXAAefc/wns\nc9m8hhFN2uJxNi5eTEN7+8C+hSUl7OvpoRQN0KlAW+V7qqoo7upi1s6dya/8LinhzJ4edpBsDKxF\nZwUgzavCJ0/mzdJS/rKykqnHH8+8AqjsI2Ek0YAZCYCIzAR6nXPviUg5+tt90zn3aCCNCYCRFdri\n8YH1APrKyjj2nHPYsXp1kigsqapi/sqVANy7fDldHR2UAFNOOomPXnopO1avpqa9nRaSW/NDjz3G\nd721ANqAxpISTjj55AlV4cPkQgA+hi6tNsn7+4lz7juhNCYAxrgRFoXhKmu69KPNZyIw7gIwwkKY\nABhGHhiJAJgnoGFEGIsFMEZNob03MIpk6xmYABijohCDhKJGNp+BDQGMUVGIQUJRI5vPwATAGBWF\nGCQUNbL5DEwAjFFRaEFCUSSbz8AEwBgVhfjewKiRzWdgfgDGqInH22hqagkECc0zA2COGckzMEcg\nw4gw5ghkGMaQmAAYRoQxATCMCGMCYBgRxlyBjXHB4gXGTi5/OxMAI+tYvMDYyfVvZ0MAI+tYvMDY\nyfVvZwJgZB2LFxg7uf7tTACMrGPxAmMn17+dCYCRdSxeYOzk+rczV2BjXLB4gbGTrd/OYgEMI8JY\nLIBhGENifgBG3vmgOw0V8v2ZABh55YPuNFTo92dDACOvfNCdhgr9/kwAjLzyQXcaKvT7MwEw8soH\n3Wmo0O/PBMDIK5k4vsTjbdTULCMWq6emZhnxeNt4FXPM1yx0p6iMjIAiciJwL3A04IDvO+cas1Ew\nIxr4hrCmpuUBx5f5wxrI8mFcG8s1x3p/OcM5N+Y/YBZwurc9Bfg9cGoojTOMbFNdvdSBG/RXU7Ns\n2HPXr9/iqquXugsuWOGqq5e69eu3jPs184FX94aswxn1AJxzO4Gd3vY+EXkROA54MZN8DWM4xmpc\ny6TnUOgGvbGQNRuAiFQCZwBPZytPw0jHWI1rmUzLFbpBbyxkxRFIRKYA/wEsds7tCx+vr68f2I7F\nYsRisWxc1ogwdXXVtLcvTarMalybP+R5mbTiY71mrmhtbaW1tXVU52QcDCQik4H1wCPOue+lOO4y\nvYZhpGIsUXM1Nctobv5Wiv3L2bDhtnG5Zr4Y92hAERHgx8Bu59xX06QxATAKhlQ2gKqqJaxcWUCW\n+SyRCwE4H2gDnkenAQH+wTm3IZDGBMAoKCZSK54Jth6AYUQYWw/AMIwhMQEwjAhjAmAYEcYEwDAi\njAmAYUQYEwDDiDAmAIYRYUwADCPCmAAYRoQxATCMCGMCYBgRxgTAMCKMCYBhRBgTAMOIMCYAhhFh\nTAAMI8KYABhGhDEBMIwIYwJgGBHGBMAwIowJgGFEGBMAw4gwJgCGEWFMAAwjwpgAGEaEMQEwjAhj\nAmAYEcYEwDAiTMYCICI/FJG3ROR32SiQYRi5Ixs9gB8B87OQj2EYOaY40wycc4+LSGXmRRl/4i1x\nGtc00u26KZVS6q6ro3Ze7aA0y5uW88quV5A+ofLoSi674DKeevGpIc8bzTWC6Xa8tYOdb+/k2GOP\n5bgZxw2Zt2Fkm4wFYKIQb4mz+I7FtJ/RPrCv/Q7d9itcvCXOwu8sZOfxO+E9YBLseWkPW9/YClcw\n6DwgqbKfe+q5rH5y9ZDXGHSdncAU2L13N9ve28bz33qeu7nbRMDICeKcyzwT7QGsc859LMUxt2LF\nioHvsViMWCyW8TVHS80NNTRXNkMH0I4OfrpgyoEpfPwTH6dUSnl799tsnbFVj1/knfhoYDtA+U/L\n6ZvSR8/FPYl9Py/nwJUHBqU989kzufS8S1n14Cp6i3rZ994++s7tg98CU0L5PwpnTD2DZ+PPZufG\njcjQ2tpKa2vrwPdvfvObOOdkqHNyIgDZuEamxP4mxhbZkly5AR4GioAKkLcEV+7g8sDxx4ALU2R4\nP3AUKiT9QBXwSpq0P/XSBNPvAw4A1w5OPn39dN595t1R3Z9hhBGRYQUgMkOAUikdXPk7SGqBHQ7W\nA78GOtHK+k6aDB2DWm4GN/5KH3D0KNIXpdlvGFkmG9OA9wFPAieLyGsickPmxco+ddfVUba/LHln\nWBAALgFe9fZfCJwFrAulWQecG9p3EVrR14b2Pwx0p7jORUBv6rKedMxJqQ8YRpbJxizAgmwUZLyp\nnVfLqU2nspWtiZ3p5G9mYLvS+1wNHAYcgbb+/v4OkmwK7AfiwPtoSz4JKPHS+ef4FAEbSJpEnfXE\nLP7xln8c4V0ZRmZEyhPwtkW3UbW1KrGjP03CsMmiEpjhbf81UOFtd5DoRVwIXIeO83cBxwDXAJ/x\n9rd76YOUAqcAm6FiXQU1r9Zw9y02A2DkjkgJQO28Wj73yc8xIz6Dwx85nPJd5RSvC3WC1gKzU5wc\nFIupwIOAb6jvCBy7HJhM6i7/y4HvG1BD4K/QKcdDsGjBIqv8Rk6JjBEQdP599ZOr2V27e2Df5J9N\npvyn5fRM7qHv8D6tlO0kd9c3oeP1InRMPwVt2X0e9T79c0rSFOB9dFbhHeAvSBKKLrpY+J2F5gNg\n5JRI9QAa1zQmOekAHLrqEAfkAH09ffAGWimrgM1oZY2jY/tDwEFUBIZr3XtIzeGoCJSTMDQG2Hne\nTpruaxr1fRnGWIlUD6Dbdac+MJPE/P2DaEteSfIYH7SlTyeZ/mzrJlRE4kCwIV8L7EXH/aWooPwC\ntSdUMdB7ONh/cIR3YxiZEykBKJXS1AeCRj9/2u9SBk8TVgHPpMn8TVQ8zkLFoAoVgR6059CJ+gIE\nXIp51Evnd0oqoWxSaKrSMMaRSA0B6q6ro3xDefLOTSQb/SpRQXgQZE/IiaoSHbs/nLxbHhb4JFr5\n20lU6o8C04BZ6BTiFcnncRHwAgNDiFlPzGLRgkVjujfDGAtZcQUe8gIF4grsM+fiObzwzgvqnNOH\nVtrKUKLNMHXvVOiHzk93Ds7kYXQm4H2YKlPp6uqi/7PeNEEHag94HR0KVKDXmgRclaJADwFngvxa\nWPe9dWYANLLGSFyBI9UDADj+mOO1ZfYt9b8ieRpvA/A6dPZ00lnUqeP0IGtRQ+B7wH7o3NNJv/Sr\nKHSgYjIbmI56C1agMwZHpCnQEcDLUOSKrPIbOSdSNgCAGSUz4HF0PF6Odtd/C8VPFtNLr/4i13uJ\nO4AtwL+jguFr6bkkeg3BcfxvvX3tqA0hGElYhcYZXBIozCbgQ8CzUDopjX3CMMaRSAlAvCXO2ufX\nqoeez6PAaVD+bDmdh3cmKmwHWpEvYLAxMDjv7xsGZ6IGvyegqKyIPvqS+1eVqEBsRoXEoZW/Us85\n5uhjsnWbhjFiIjUEaFzTyIH5oRA8zwDX2depLryPkjz9lypgyJ/399N9Bp1G/AwwDfr29Wm6sKvx\naahIXIi6FFeiQ47zYX/FfuIt8Yzv0TBGQ6QEYMdbO7SCP0aiooO2yDNQ/32/0vvOPOl+oX1oyz8p\nlNdlaDyAPzR4NHTeXhJORptRAyHqBHTrqlvHeGeGMTYiMwSIt8T54+4/wsWBnX7l/DM6xp9MouI+\n5R1LFTDUgYpGOnfgHtR1+AU05n81KjD7gKtT5LdZz9u+YzvxlrgZA42cEZkeQOOaxqTluwBt7VuB\nMnSufjoJg95e1LGni8HrAfwKNfKF8/LdgQ9Hg4LKSMQQ+COPjhSF84yLBysOmiuwkVMi0wNI6Qbc\ngVb64BJgfg9gJ2rY60fn/O9HnXm60QqeCkHH9Kd436vQVj+cPyT7HjgGZgTMFdjIJZHpAaR0A24n\nuXJCoiU/DjXWXYS6+U5BHYcg/VJebzEwph82f5+H0SAkb0bAXIGNXBIZAai7ro5Zv5yVvHNPmsT+\nNB0k1g28BLgSHff3oi19kE2oX8HlJCp4uqjA94GfobECoO7ClVC+odxcgY2cEpkhQO28Wu7mbm5d\ndSvP/ek5+vv6tWufijfRoYG/cOcloeNXoy13eE7/Fe+4oMLxXpr8u4FzSAwD1kD52nJuWXCLGQCN\nnBK5WACAKadNoevTXYPDfUEr9mkkKudDwKdTZJJqufDN6Pz+g2i47ydT5L8O+BhJNoCKX1TwwP9+\nwCq/kVVsWfA0dPd7A/VKb4ffku8AzifZQFdEYh0Af/3/StRpKIjv1rsJDTB6OkX+jpQLipZMLqFx\nTSOAiYCRUyIpAGVFZexjn36pJFEhHyC58negi3eE3YD/H9qNfxAVhm50vN+NBhpBwoYQzB/vHD/v\nl4D5sIc9NNPM89+x14IZuSUyRkCfeEtcjXhhD71NJFvwQbvv4fceX4Qu8HE8ahC8Cl31dyYJo+Iz\n6Bg/fI216HoCm1FfglDe5g1o5JpI9QDqb6/n9jW3c6D0AOwmuWvuv4nLXw0I0stjMalfKLIZSraX\n0H+on95K760fwWvsQw2MwcjCEK+89UrqA4YxDkRGAOItcW7/+e2Jl3d2oC11Kbpc13lewl+h3fRi\nhn7VVyoE9Tb8iZdHEVrZ56DDgDj64pA9pF85OF3ehjEOREYABkUCVnqfv0Yr4yuokW822lJ3o/P1\nP0fn/302oRb+VPjj/gqSZw42AM955/WTWCUo/ObhTTC9fPqo7sswMiEyApDSFfhXaGUM+vU/TKIV\nPhEdGoTn+z/kpbucxFTi+yReARaeOTgFDQzqJxFA1IGKT2gYcnhlOj9jw8g+kRGAJFfgDmAbWvEO\nIxGg81u08gbddx9C5/bDPIl280tC6deiAhKeOegjeVmwSu/TX0xkJ/AhmHbUtBHekWFkTjbeDjxf\nRF4SkT+KyDeyUajxoO66On0vYAfaYl+Crgx0EfA88ISXMOz1l65BngwcyWBf/ytQL8IgF6F2htno\n9f01CdpRG8SF3nXOtlgAI7dkJAAiUgSsQie0PgIsEJFTs1GwbFM7r5aVN69kxgszBlvwL0ODf4oY\nHK7rr+UX5GG05U/n65/Kwu/3tYIvE73IS7sW+KjFAhi5J9MewNnAn5xzHc65Q2jQbLhNLBhq59Uy\n5yNzUh8UtPX3A3k60Jb6FbT19i37m9GKX46O+1ORyvN5CtrdD4vPpehsw8twdPHR5gRk5JRMBeB4\n4LXA99e9fQXLsG8Heo/Br/1egC7zVYraAz6DvuAznbPPO6F9m1APwSlpCjVL893XvW/E92EY2SBT\nASisKJ8RMGALCBJ8O1ApqVvq+SR7CpaTWBU4uMbfJLRf5O97kMTqv+nm/r0FQWbNnJUmgWGMD5nO\nAuxAJ8t8TkR7AUnU19cPbMdiMWKxWIaXHTt+F/vWVbeyfcd2DlYcHKig5RvKmT1zNh27Ouiia/DJ\ngRa8uMd7j0Alyb7+3vp+Ayv+lgSOV0HxfxbTe3FvIv06dHrwNDhBTsj8Bo3I0traSmtr66jOySgc\nWESKgd/+IoF0AAAEmElEQVSj7eUb6Mz2Aufci4E0BRcO7BNvidN0XxMH+w9SNqmMRQsWUTuvlpob\namiubB58ghfuW/VsFZ8773OsfnJ10uvGpz86na63u+g5pkdb9dkwq30Wx00/jqnTp1I2qYxzTjmH\ndb9clxCf2UCl5rnyKyvNBmBkjZGEA2e8HoCIfAr4HmpDv8c59+3Q8YIVgHTEW+IsvmNxUuUu31DO\n7OmzOWHWCQNCkUpAgJSikuoaI0lnGGMlJwIwgkJMOAEAq6DGxMcEwDAijL0d2DCMITEBMIwIYwJg\nGBHGBMAwIowJgGFEGBMAw4gwJgCGEWFMAAwjwpgAGEaEMQEwjAhjAmAYEcYEwDAijAmAYUQYEwDD\niDAmAIYRYUwADCPCmAAYRoQxATCMCGMCYBgRxgTAMCKMCYBhRBgTAMOIMCYAhhFhTAAMI8KYABhG\nhDEBMIwIYwJgGBFmzAIgIp8RkRdEpE9EzsxmoQzDyA2Z9AB+B1wJtGWpLDmjtbU130UYRCGWCQqz\nXFam7DFmAXDOveSc+0M2C5MrCvFhFWKZoDDLZWXKHmYDMIwIUzzUQRFpAWalOLTEObdufIpkGEau\nEOdcZhmIPAZ83Tn3bJrjmV3AMIwx45yToY4P2QMYBWkvMlwBDMPIH5lMA14pIq8B5wBxEXkke8Uy\nDCMXZDwEMAxj4pKTWYBCchoSkfki8pKI/FFEvpHPsnjl+aGIvCUiv8t3WXxE5EQRecx7ZttEpK4A\nylQmIk+LyHMisl1Evp3vMvmISJGIbBWRgjGMi0iHiDzvlevX6dLlahqwIJyGRKQIWAXMBz4CLBCR\nU/NZJuBHXnkKiUPAV51zH0WHeDfn+3dyzh0ELnTOnQ78FXChiJyfzzIFWAxsBwqpO+2AmHPuDOfc\n2ekS5UQACshp6GzgT865DufcIeB+4PJ8Fsg59ziwJ59lCOOc2+mce87b3ge8CByX31KBc26/t1kC\nFAHv5rE4AIjICcDFwN0MYQzPE8OWJ2qOQMcDrwW+v+7tM9IgIpXAGcDT+S0JiMgkEXkOeAt4zDm3\nPd9lAv4F+HugP98FCeGATSLyGxH5YrpE2ZoGnChOQ4XURSt4RGQK8B/AYq8nkFecc/3A6SJyOLBR\nRGLOudZ8lUdELgF2Oee2ikgsX+VIw3nOuTdF5CigRURe8nqbSWRNAJxz87KV1ziyAzgx8P1EtBdg\nhBCRycDPgNXOubX5Lk8Q59z7IhIHPgG05rEonwQuE5GLgTJgmojc65y7Po9lAsA596b3+baI/Bwd\n/g4SgHwMAfI5TvoN8JciUikiJcA1wC/yWJ6CREQEuAfY7pz7Xr7LAyAiM0XkCG+7HJgHbM1nmZxz\nS5xzJzrnTgKuBTYXQuUXkcNEZKq3XQFUo4b4QeRqGrAgnIacc73AV4CNqNX2Aefci/koi4+I3Ac8\nCZwsIq+JyA35LI/HecDnUEv7Vu8v3zMVxwKbPRvA08A659yjeS5TmEIZYh4DPB74rdY755pTJTRH\nIMOIMFGbBTAMI4AJgGFEGBMAw4gwJgCGEWFMAAwjwpgAGEaEMQEwjAhjAmAYEeb/AxH1A1REiTAK\nAAAAAElFTkSuQmCC\n",
              "text": [
                "<matplotlib.figure.Figure at 0x7fd7ac14ebe0>"
              ]
            }
          ],
          "prompt_number": 39
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Construct the interaction matrix *from* the target proxies *to* the targets as `target_proxy_mat`.\n",
            "\n",
            "**A note on terminology:** The `target_proxies` are *near* the targets but *stand in* for far-away sources."
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "target_proxy_mat = interaction_mat(targets, target_proxies)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 61
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Check its numerical rank and shape:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "numerical_rank(target_proxy_mat, eps)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 62,
              "text": [
                "24"
              ]
            }
          ],
          "prompt_number": 62
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "target_proxy_mat.shape"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 42,
              "text": [
                "(200, 25)"
              ]
            }
          ],
          "prompt_number": 42
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Now compute an ID (row or column?):"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "idx, proj = sli.interp_decomp(target_proxy_mat.T, nproxies)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 43
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Find the target skeleton as `target_skeleton`, i.e. the indices of the targets from which the remaining values can be recovered:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "target_skeleton = idx[:nproxies]"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 44
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Check that the ID does what is promises:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "P = np.hstack([np.eye(nproxies), proj])[:,np.argsort(idx)]\n",
            "tpm_approx = P.T.dot(target_proxy_mat[target_skeleton])\n",
            "\n",
            "la.norm(tpm_approx - target_proxy_mat, 2)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 45,
              "text": [
                "5.7389031370755866e-15"
              ]
            }
          ],
          "prompt_number": 45
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Plot the chosen proxies:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "pt.plot(sources[0], sources[1], \"go\")\n",
            "pt.plot(targets[0], targets[1], \"ro\", alpha=0.05)\n",
            "pt.plot(targets[0, target_skeleton], targets[1, target_skeleton], \"ro\")\n",
            "pt.plot(target_proxies[0], target_proxies[1], \"bo\")\n",
            "\n",
            "pt.xlim([-1, 5])\n",
            "pt.ylim([-1, 5])\n",
            "\n",
            "pt.gca().set_aspect(\"equal\")"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "display_data",
              "png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAEACAYAAABccqhmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmUG1ed7z9XJZVUkrrb7T2Lx+14HgPDkoSBrCbpxCsx\nmGUITDI5zPAOMJMZ7DDvDVnshDisIcyD2A4MzGGZgUyABIiX9OB47Th7AjEE4mSGmDRJHLvddq9S\nay3d90dJ3ZJasrtbsqRO/T7n9HHp1q1bt6p8v3Xv7/5+t5TWGkEQ3Imn3hUQBKF+iAAIgosRARAE\nFyMCIAguRgRAEFyMCIAguBhvpQUopbqAQcAGUlrr8yotUxCE2lCxAAAaaNda91ahLEEQaki1hgCq\nSuUIglBDqiEAGtillPqlUuoTVShPEIQaUY0hwMVa68NKqVnATqXUC1rrh6tQriAIp5iKBUBrfTj7\nb49S6n7gPGBEAJRSEmwgCHVCa33C4XlFQwClVFAp1ZTdDgHLgN+WqERD/d166611r8NUqFOj1uvW\nW29l2bJ1OKPPwr/ly28+4bGTPe5kxzbifRoPldoA5gAPK6V+DTwJPKC13lFhmYJwUhKJ0p3XeNw4\n4XFr1ixj4cJ1BWkLF65l9eqlp+ycjUxFQwCt9UvAOVWqiyCMG78/XTI9ELBPeNzKlZcAsGnTLcTj\nBoGAzerVK0bSJ3/OqSkC1TACTjna29vrXYUxNGKdoDb16ujYx8aNO0gkvPj9adasWXbCBtne3s47\n3+nh4MF1HDz4xZF0502+4qTnW7nyknE1+GLWrFlW9pyhUKbscRO9vlqixjtWmPQJlNKn+hzC1KWj\nYx/XXfdgUaNax4YNy0/aSDo69rFp0868N/nSU96wJnrOSq6vUpRS6JMYAWthiNCCUI5ly9Zp0GP+\nli+/ud5Vqwr1vL5s2zth+5RgIKGuvB4Na/k0+vWJAAh1ZbLGvKlCo1+fCIBQVyqZlpsKNPr1iRFQ\nqDv1MObVknpd33iMgCIAgvA6ZTwCIEMAQXAxrnQEEk49jez80ujU8t6JAAhVp5Tzy8GDjiFMRODE\n1PreyRBAqDobN+4o+A8McPDgF9m0aWedajR1qPW9EwEQqk6jO780MrW+dyIAQtVpdOeXRqbW904E\nQKg6je780sjU+t6JH4BwSni9O/ecSqp178QRSBBcjDgCCYJwQkQABMHFiAAIgosRARAEFyOuwMKE\nET//+lOtZyACIEwI8fOvP9V8BjIEECaE+PnXn2o+AxEAYUKIn3/9qeYzEAEQJoT4+defaj4DEQBh\nQoiff/2p5jMQV2Bhwoiff/0ZzzOoWSyAUsoAfgm8qrV+b9E+EQBBqAO1jAW4DjiA88F0QRCmCBX7\nASilzgSuAL4I/J+KayQI42Dv1q3s3LgRbyJB2u9n6Zo1XLZqVVXLquY5GpVqOAJ9HfgM0FyFsgTh\npOzdupWdn/40X3rppZG0tX/4A8BJG25unxGPk/R6mXH22Rz/+c/5yquvjpR104svsv+RR+i5916+\n/Mc/Fpzjt089xZEnnsCbSHB0cJAEMK+5eeoKxMm+HnqiP+A9wDey2+3AthJ5qvnBU0HQNy1ePPZz\nu6DXLlmi92zZom9asKAg/aYFC/SeLVtK7vuk368fKlHWldOmjUl7CPTfBQKF58ym55+nUWAcXweu\ntAdwEbBKKXUFEACalVI/0Fp/ND/T+vXrR7bb29tpb2+v8LSCm/EmEqV3DA6y/Wtf4yt5PQOAL730\nEus2bUJrXdBrAPh2IsEtQPEcRiCVGlP8DuBb8XhB2hdh5PjceerVC+js7KSzs3NCx1QkAFrrtcBa\nAKXUpcA/Fzd+KBQAYerQiEE/e7du5YXnnmM9kAaWMdp4tWliDg+XPM4oargF+0qkxYyxqfmNZR+O\nIHiB32d/X3KS85xqolEPjz+eHnle46HawUAyC/A6oRGDfnJj/3v7+kbScu4wHaedxsXvfz+P3Htv\nyWMTHg+o0jNiB0wTksmR39fPncv5y5Zx444d3H7kyGg+vx8SCfYBD+K8/YvrYQcCE7+wKlDqeRXW\nsDTiCCSUZPnym9mx4wsl0m9h+/bP16FGsHbJEr60e/eY9A83NfF3N97Iueecw/7nn2fn177G7a+9\nNrL/xnnzWH7bbSRsm87PfY7bX3lldN/cuZy+dCmHDxxAx2KkvF4Wf/jDXHLeeezbu5fOBx7AC9jB\nIHPe/na6t2zB99prjL0z8OFAgGt/8pO6DAFKP6+T+wFIOLBQkkYM+ik39v+z+fM5/y1vgVSK89/2\nNlL//M9cf++9BGybhM/Hio9/nEuWLMGORjHWrWPd3Xejh4dJeb2868orufD880EpkvE4KpHA9HjA\n6+Uv3vpWLjn3XPpTKaYtXIgRDrPnDW/gRzffDNHomHrMmT+/buP/cs/rZIgACCVpxKCftN9fMj3p\n9UI6jWmaRCIR3nH22Zz7rncxc84cYvE44VCI5PAwpsfDpZdfztL3vpdINIoJ2JkMRiaDaRjEBgfp\ne/ll0Bo7kyHg9xMbHqbJshjo66PFslixfDmdd98Nv/rVmHqETzvtFN+B8ox3zF+MBAMJJWnEoJ+l\na9awdsGCgrTPnHYaF3/oQ0SCQSLJJEYoRGtzM0Y8TuzoUSIDA0R6eoj195McHsbOZACws428t7ub\nZF8fkdxfPE4iHod4HDMYBNMEyyIcCmFmMthK0f43f8MNRY39pvnzufzaa2t2L4op9bzGg/QAhJLk\nDH2bNt2SF3Cyoq6zAJetWkUmk+GmDRvwxuOkvV4WX3UVi97xDidDIIAVDpO0bTxeL8lEgkw0Siyd\nJtbfT6ypCXPmTJIeD4bXix2JQDSKrRSxaBR8PpqCQfypFJHhYUylGM5kiMfjeAwD2zCwmpq4/PLL\nUbbNtd/6FpFDh0gqhRUK4THNut2bUs/rwQdPfpwYAYUpw+7Nm3nwa1/DjMdJ+ny0f+xjXHD++Zge\nD5GhISzLwh4exradYYqtNcbQECZAMslwOk1w5kxioRDTZ81ioLsbIx7H7u/HiMcZSCaxLIt4Xx+W\nYeCxLJTf70wvNjVhhsPYSmH4/Tzy+OPsue02vpznQbh2wQKW3nlnw3gDjicYSHoAwpRg9+bN7Pr0\np7kjzzX3xq4uWLuWiy+5BLO5GcvjwbQsIoODGEoRO34cw+ulxbIAnC68zwepFJHBQaxQiGQ8TuzY\nMaYrhW9gAE8gAD4fhmWRjscxtSbm94PXSyQaxfJ6MZRi1/e/X+A+DPV3BJoMYgMQpgQ7N2wo8MsH\nuP2113j4hz/E0Jrevj4i0SiR4WFs08SOx0lGIqTyxv2GYWBaFmY4jGGaRHp7SXR1YcXjJIeGSPb3\nM9jdTTKRwE4mR/wGFGAaBi0zZmAGg6QSCejvL1nPejoCTQYRAGFK4CszBWhkDXbTLIuW2bMxLIvk\n8DAx28ZqaSEdixE7fpyBwUHw+UhqDcEgA8eOkezpwRgYIKwUyeFhpoXDTPP5SPX3o5NJzFCIlmnT\n8KTThHMzDdOmERsYwMj2KoqplyPQZBEBEKYEqTJTgHGvF8MwsBMJ7KxIzJo+3WmomQwpr5fhRIJE\nIsGxY8eIALG+PqbPnk0ymcQfDDIYj2O1tDCg9YjHoGFZmMEgMa2xstOBht9PMpkk0NzMolWruGHO\nnIK63DBvHktWrz7Vt6KqiA1AmBLMvfBC/u7hh/m2PeqH8EmPhz9ZsIDk4CC+pibsRIJkNAqGQTIW\nw+v1Ys2bhx2PMzwwQCAUIq0105uaMDIZrDlzoK8Pn9ZgGLQEg+hYDJ9tYzY30zs4iHfGDLRSWIEA\npmliK0UmEmHJBRcAsLajAyORIBUMcum1106p8T+IAAhThKNPPcVf2za34ATv2MA1mQwdBw8C0BIK\nEYnHifb3k8pkGBgaYu7MmZBOYwSDhKdNwwwG6R0exvB6MTMZrKYm0BrDMEj09WGYJoNKYTY3E0+l\naJk+HTM77u8ZGnKmGA8dorW5me0PPMDj//VfeFMp0oEA5115JUs+8pF63qJJIQIglKTRIgG9iQSX\nMDZsd3s6TTI7hmdoiGafj0g0yqymJsxEArTmaE8PwVmz4NgxBrTG0Jpkdk7f9vux/X76UilagkHC\npgmpFCaQVIr40BDE43hbWug9fBjTMNje2cmj3/0uXzl8eKQeN3R307pwYUP0ANav/yZ33fXQuPKK\nAAhjaMRIwHJuwFgW4Wz3347HSSpFePZs7Hic2NAQ9uAgls8H3d2Eg0FSSpHq7iaSTpNpaaH5jDOw\nLQuruRlzcBArECA5PIxh2xCJkEylaJ01CwYGGBocJGOaPHr//QWNH+ArL7/cEFOA69d/ky9+8VnS\n6Z8ApSMj8xEBEMZQ/tNTt9RNAJauWcPHnn2WM3p68OKsBXCotZVVy5Zhp9OgNVYwSG80StgwMIJB\nopEIvmCQaCxGwOMh5vEQSCQwTZOAz0dMKexEAqutDaOvD8s0MT0eDMOA4WFsw8C2bYxIBJJJvEeP\nOvsGB0vWsRGmAO+666Fs4x8fIgDCGBoxEhCgRamCMNzrgHBzM9g2diZD38AAw7bNrFgMIxCAYJC+\nWMzp5ns8xGMxWk3TceltbQWvl3BLC5GhIQzDwAiFHB8Ay3Ks/l6v06hjMexYjGbTxI5GSXtKT541\nwhRgOl16erIcIgDCGBoxEnDnxo3cefRoQdqGvj5u+tnPOO/yyyGTobmpiUAmQ2R4GCIRhjMZmpqb\nicViTJ82jZhSGOm082aPx4kpBT09xAIBzOZmkpkMVjDoFK4Ug8PD2B4Ptm1j+v3EYjEMpVi0ZAl/\n++qrnDkwMNobmTWLjzbAFKDXG5tQfvEDEMbQiJGA5dYC8KVSGF4v4XCYWCZDMpEA0yRjGGjDIBaJ\nYLS2EsvGo9i2jZ1IEEsksDwewraNEYthmSbJTIaY1iSVgnAYPWcOvmnTMEyT2NAQ6UiEdDKJ3+9n\nWrY3sh74AjCtTK+g1nzqU5fi9f79uPNLD0AYQyNGApYzAsZ9PiLJJMmhIUzDwAqFnDX+/H58QDgY\npCeRIDlnDtFkktjAAC1+P4HWViyvl97BQcxZs0j292NNm0bStjGCQQylaJk+nd4XX+T4iy/ijccx\nDYOk1uzaupU7i1yBv97d3SBGwH8Avsldd/0Vx4+fPL8IgFCSlSsvqfsCoPksXbOGtX/4Q8Gqvtef\ncQbLP/lJwqaJBRzp7iaeydCsNYbPx3A67RgIDYNUdzczZ80i2dyMHY0SjcdJeTw0z5yJ6fORtG2I\nxTCyUX85Ws48k5ht0xeJMFNrWgIBgmWiWxvBCAiOCKxf/w8odXJjoAiAMCXIvVnXbdqEEY+T8vlY\n/olP8K7LLiNy/DjJ/n5awmFiR49CJsOx48cZymQIHDuGZVlY06ePxAJYM2eSevllPKkUyaEhaGoC\nv98JKy4aapg+H8asWdh/8idw/Di2YZAqY+xrBCPgRJH1AIQpSWxgACMWw06n6X35ZaxkkmhvL/He\nXk6bNm0kEjAQChHzejFDIYzWVpLhMPbgIC0eD3Y8TjgYpDeVwpo/H7xekj4fVnPzSARhMhLB7unB\nHh6G/n5MpXjo0Ud55N/+jTvyjJI3zJvHkjvvZOkHP1ivWzIGWQ9AeF1i2zZ2NIpp29iDg4STSYxo\nFAMINTUxMDyMSqfxmSbhUMhx621tJQkMHDlCSzjMQCxGKholkkhg+P309fbSPHOmIyqZjNPgvV7H\n7TgWwz58mLDPRzKT4W1vfjMHFi3iyp07sWybuM/Hog98gAsXLarznZk4jWG6FIQJYMfjWMGgM1ef\nTGL5fFjhMJFkEqulBbOpCU84jBEOY/t8+FpaSGYbtd+2MTIZjHicZtNkejhMi1JEurqIHDmCaRgY\n8TiWbZN87TWSx48Ti0RQtk10aAgjk+HRRx+l/5FHuG9oiB8MD3PvwACH77+ffb/4Rb1vzYQRARCm\nHtkAHiMUIpaduhtIJDCmTXPG8paFbZoQDBJJpTADAad7bxjEAgFiqRTTW1uxQiEig4OgNYFMBiud\nxo5GifT1EevpwYzFSB47Rks6TbNl4WlqYjiR4IGf/hTv0aOsB27G+SrQl195hYe+97363pdJIEMA\nYUqxd+tWHvz61/HF48SV4l2rVrHiXe+ClhZiySTHhoYY6uvD9PkAMIJBbK/XcfI5/XRmzphB5H/+\nB8PjwU4kCAeDJE2T8Ny5zic+o1GSR45geb0M9PbyyuHDzA6F8Hu9pAMBfvP887QcO1bgkZjzmDDT\nk1uau56IAAhThr1bt7Lj05/m9rypwBu7urC05rzzzsMOBEgNDhJSivCMGZDJEPV4CM+aheHxYIXD\njo9AUxMD2dWAUYrwzJkYiQS9PT3M0hrT4yHa20uqu5tgfz9NQCKdxtvczO7Nm/lWUUPPfSA0k/Mi\nnEKIAAhlabSQ4B0bNvDloq/73n74MJ+5/37eftFF9Pf0EEqnsSwLI7t8V5PXSwScNK0xtWb69OlE\nUimMQADLsiCdJpJOEzMMYrEYg0NDMDCApRQzmpsZGBzEFwrhAbxl3vJdpsn/rrMrcPHzGg8iAEJJ\nGjEk2Bsr7eduplJMmz2b4b4+pmcX+TCBSH8/1owZxPr6aDFNbL+fgUjE6eo3NTmLgwwOElYKMxQi\n2NdHKrsakBoawptK0dTczEAwiAKGtS67NFlm7ty6egFO9uOgYgQUSlI+JHhnnWoEqTIf3kgoxfHf\n/55Yby+RWIxk1i/fzC4UEotEnFBfrbFsG0trWpqaaAkGCXq9GH4/6cFBPIDKZCCVwhMK4fX5GBoe\nZt/zz3P7N7/JNzdsoGdoiH/I8xQEWAuYsRh7t249xXegPKWe13ioqAeglAoADwF+wAS2aK1vqqRM\noTFoxJDg2eefz98/+STfynO5veG007h4xQq86TQhyyJp25haE0mn2fvYY+zetg2/x4MnHObya65h\n6eWXY9s2kZ4eIseOMTcUcoJ/Uim8to0NxKJR0rEYgUSCp59/nl9u3cq/9PSMnPNKw+AfgVk4S5Ot\nAC7p6alrLEBdPg6qtY4rpS7TWg8rpbzAI0qpRVrrRyopV6g/jRYSvHfrVnruu4+r4/GRdQEPmCaz\n/+zPeHzXLp7etg1tmlz47ndz4Tvfyd59+/jVD3/InXkr99zY1cVz+/dzZP9+PPE4KdPkig99iLe9\n7W088cQTPLptGyqRIO71ct7ixSw65xwef+IJ/l9e4we4L7s24fqiOtYzFmCyHwet2AagtR7Obpo4\nz6W30jKF+rNmzTIOHlxX0K10QoJX1KU+OzduHDEA5iwQ+5JJ/vPxx7kjz3//htdeI+338/TevXy1\naNmuKw4d4j///d/5djI5knbToUM8vWgRvQ8/zO3d3SPpN/b04PH7sbylm0ipflA9YwFKPa/xULEA\nKKU8wDPAQuBftdYHKi1TqD+NFhJcaj2AHcC3i9K/cuQI13d0YGZ9+cfkz2v8AF8+coQrt2/nvkik\nIP32I0e4fs8ePGWm9g74/ZB37psWLGBZHWcBJvtx0Gr0ADLAOUqpFuBBpVS71rozP8/69etHttvb\n22lvb6/0tEINaKSQ4FLrAZT7z+u1bTKh0LjzWyXEAsBn2yy59lrWvvoqX3rllZH0G+fP56Krr2bd\n0087HxcNBFi2enXd1wIIhTJccEGub2LURgByaK0HlFIdwDuAzvx9+QIgCJOh1HoABwIBKDHu1s3N\nXP6JT7D2j3/kS3nfEyyXP5b1GiwmEw5z8bvfTQy48Qc/wJdKkbEsljdAYy9F8cv1tttuO+kxFYUD\nK6VmAmmtdb9SygIeBG7TWu/OyyPhwEJV2Lt1K7uy6wHYgQBzzj+fI/fcUyAKNy1YwLI77+SSlSvZ\n87Ofsedf/xVvIoEdDDL3ggvG5L9x/nzmXn013T/+cYGT0Q3z57Pkjju4dOVKjEDAWQ14ijGecOBK\nBeCtwH/g+BN4gB9qrb9alEcEQDhlFIvCkpO8ncvln2g5U4FTLgDjrIQIgCDUgfEIgHgCCoKLkVgA\nYcI0WpCQG6nWMxABECZEIwYJuY1qPgMZAggTohGDhNxGNZ+BCIAwIRoxSMhtVPMZiAAIE6LRgoTc\nSDWfgQiAMCEa8buBbqOaz0D8AIQJ09Gxj02bduYFCS0VA2CNGc8zEEcgQXAx4ggkCMIJEQEQBBcj\nAiAILkYEQBBcjLgCC6cEiReYPLW8dyIAQtWReIHJU+t7J0MAoepIvMDkqfW9EwEQqo7EC0yeWt87\nEQCh6ki8wOSp9b0TARCqjsQLTJ5a3ztxBRZOCRIvMHmqde8kFkAQXIzEAgiCcELED0CoO693p6FG\nvj4RAKGuvN6dhhr9+mQIINSV17vTUKNfnwiAUFde705DjX59IgBCXXm9Ow01+vWJAAh1pRLHl46O\nfSxffjPt7etZvvxmOjr2napqTvqcje4UVZERUCk1D/gBMBvQwL9prTdWo2KCO8gZwjZtuiXP8WXF\nSQ1k9TCuTeack72+mqG1nvQfMBc4J7sdBv4beFNRHi0I1WbZsnUa9Ji/5ctvPumxDzzwkF62bJ2+\n9NJb9bJl6/QDDzx0ys9ZD7Jt74RtuKIegNb6CHAkux1RSj0PnA48X0m5gnAyJmtcq6Tn0OgGvclQ\nNRuAUqoNOBd4slplCkI5Jmtcq2RartENepOhKo5ASqkw8FPgOq11pHj/+vXrR7bb29tpb2+vxmkF\nF7NmzTIOHlxX0Jgd49qKEx5XyVt8suesFZ2dnXR2dk7omIqDgZRSPuAB4Bda6ztL7NeVnkMQSjGZ\nqLnly29mx44vlEi/he3bP39KzlkvTnk0oFJKAf8BHNda/1OZPCIAQsNQygawcOFaNmxoIMt8laiF\nACwC9gHP4kwDAtyktd6el0cEQGgoptJbvBJkPQBBcDGyHoAgCCdEBEAQXIwIgCC4GBEAQXAxIgCC\n4GJEAATBxYgACIKLEQEQBBcjAiAILkYEQBBcjAiAILgYEQBBcDEiAILgYkQABMHFiAAIgosRARAE\nFyMCIAguRgRAEFyMCIAguBgRAEFwMSIAguBiRAAEwcWIAAiCixEBEAQXIwIgCC5GBEAQXIwIgCC4\nGBEAQXAxFQuAUup7SqlupdRvq1EhQRBqRzV6AN8HVlShHEEQaoy30gK01g8rpdoqr8qpp2NnBxvv\n2UhCJ/ArP2uuXsPKpSvH5Lll0y28dPQllK1om93GqktX8fjzj5/wuImcIz/foe5DHOk5wmmnncbp\nM04/YdmCUG0qFoCpQsfODq77xnUcPPfgSNrBbzjbuQbXsbODj3/14xw54wj0Ax7oe6GP/a/th/cz\n5jigoLFf+KYLufuxu094jjHnOQKE4fjgcX7X/zue/cKzfIfviAgINUFprSsvxOkBbNNav7XEPn3r\nrbeO/G5vb6e9vb3ic06U5R9bzo62HdAFHMQZ/EQhHAvzF+/4C/zKT8/xHvbP2O/sX5w9cHfedh7W\nvRZ22CZ5RXI07X6L2AdiY/K+/Zm3896L38td991F2kgT6Y9gX2jDb4BwUfm74dymc3mm45nqXLjg\nGjo7O+ns7Bz5fdttt6G1Vic6piYCUI1zVEr737bzkHqosHEDbAEMIASqW6EtDe/L278XuKxEgT8G\nZuEISQZYCLxUJu+92Tz5+SNADPirsdlbH2il9+neCV2fIBSjlDqpALhmCOBX/rGNv4uCN7BGwwPA\nU8AQTmM9VqZAzZg3N2Nf/g42MHsC+Y0y6YJQZaoxDfgj4DHgDUqpV5RSH6u8WtVnzdVrCAwHChOL\nBQHgPcAfs+mXAe8EthXl2QZcWJS2GKehby5K3wIkSpxnMZAuXdcFcxaU3iEIVaYaswBXVaMip5qV\nS1fypk1vYj/7RxPLyd/MvO227L93A0FgGs7bP5feRYFNgWGgAxjAeZN7ADObL3dMDgPYTsEk6txH\n5vK56z83zqsShMpwlSfg51d/noX7F44mZMpkLDZZtAEzstuXA6HsdhejvYjLgKtxxvlHgTnAR4Ar\ns+kHs/nz8QNvBPZAaFuI5X9czneulxkAoXa4SgBWLl3JNRddw4yOGbT8ogXrqIV3W1EnaDNwVomD\n88WiCbgPyBnqu/L2vQ/wUbrL/4e839txDIFP4Ew5pmD1Vaul8Qs1xTVGQHDm3+9+7G6Orzw+kub7\nmQ/rXoukL4ndYjuN8iCF3fVdOON1A2dMH8Z5s+fYnf03d4xZpgIDOLMKx4D5FAhFlCgf/+rHxQdA\nqCmu6gFsvGdjgZMOQOovU8RUDDtpw2s4jXIhsAensXbgjO1TQBxHBE72dk9SmhYcEbAYNTTmceTi\nI2z60aYJX5cgTBZX9QASOlF6x0xG5+/vw3mTt1E4xgfnTV9OMnOzrbtwRKQDyH+RbwYGccb9fhxB\n2YpjT1jISO8hnomP82oEoXJcJQB+5S+9I9/ol5v2ey9jpwkXAk+XKfwwjni8E0cMFuKIQBKn5zCE\n4wuQ51LM7my+XKekDQKeoqlKQTiFuGoIsObqNVjbrcLEXRQa/dpwBOE+UH1FTlRtOGP3LYXJaouC\ni3Aa/0FGG/WbgWZgLs4U4vsLj2Mx8BwjQ4i5j8xl9VWrJ3VtgjAZquIKfMITNIgrcI63XPEWnjv2\nnOOcY+M02raiTHugabAJMjD0waGxhWzBmQkYgCbVRDQaJfPX2WmCLhx7wKs4Q4EQzrk8wF+WqNDP\ngbeDekqx7c5tYgAUqsZ4XIFd1QMAOGPOGc6bOWepf4LCabztwKswlBxiyBhyxun5bMYxBPYDwzDU\nN0RGZRxR6MIRk7OAVhxvwRDOjMG0MhWaBvwBDG1I4xdqjqtsAAAzzBnwMM543MLprv8GvI95SZN2\n7shHs5m7gIeAf8cRjJyWXshoryF/HP+bbNpBHBtCfiThQpw4g/fkVWYX8KfAM+D3lLFPCMIpxFUC\n0LGzg83PbnY89HLsBs4G6xmLoZah0QbbhdOQL2WsMTB/3j9nGJyJY/B7BIyAgY1d2L9qwxGIPThC\nonEaf5tzzJzZc6p1mYIwblw1BNh4z0ZiK4pC8LIGuCF7yHHh3U3h9F+pgKHcvH8u35U404hXAs1g\nR2wnX7EVWWA3AAAKpUlEQVSr8dk4InEZjktxG86QYxEMh4bp2NlR8TUKwkRwlQAc6j7kNPC9jDZ0\ncN7IM3D893ONPufMU+4ORXDe/J6islbhxAPkhga7i44bZNTJaA+OgRDHCeizd312klcmCJPDNUOA\njp0d/P747+GKvMRc43wZZ4zvY7ThPp7dVypgqAtHNMq5AydxXIefw4n5vxtHYCLAh0uUt8c57sCh\nA3Ts7BBjoFAzXNMD2HjPxoLluwDnbd8JBHDm6lsZNegN4jj2RBm7HsATOEa+4rJy7sAtOEFBAUZj\nCHIjj64SlcsaF+OhuLgCCzXFNT2Akm7AXTiNPn8JsFwP4AiOYS+DM+f/YxxnngROAy+FwhnTvzH7\neyHOW7+4fCj0PdCMzAiIK7BQS1zTAyjpBnyQwsYJo2/y03GMdYtx3HzDOI5DUH4pr25GxvQnLT/H\nFpwgpOyMgLgCC7XENQKw5uo1zH10bmFiX5nMuWk6GF038D3AB3DG/WmcN30+u3D8Ct7HaAMvFxU4\nAPwMJ1YAHHfhNrC2W+IKLNQU1wwBVi5dyXf4Dp+967P8+sVfk7EzTte+FIdxhga5hTvfU7T/wzhv\n7uI5/Zey+xWOcPSXKT8BXMDoMOAesDZbXH/V9WIAFGqK62IBAMJnh4l+MDo23Bechn02o43z58AH\nSxRSarnwPTjz+/fhhPteVKL8bcBbKbABhLaG+MntP5HGL1QVWRa8DIlMdqDelk3IvckPAYsoNNAZ\njK4DkFv/vw3HaSifnFvvLpwAoydLlK8puaCo6TPZeM9GABEBoaa4UgACRoAIEedHG6MN8icUNv4u\nnMU7it2Af4XTjb8PRxgSOOP9BE6gEYzaEPLLJ3tMruwXgBXQRx872MGzX5XPggm1xTVGwBwdOzsc\nI16xh94uCi344HTfi797vBhngY8zcAyCf4mz6u9MRo2KT+OM8YvPsRlnPYE9OL4ERWWLN6BQa1zV\nA1h/x3ruuOcOYv4YHKewa577ElduNSAoL49eSn9QZA+YB0wyqQzptuxXP/LPEcExMOZHFhbxUvdL\npXcIwinANQLQsbODO+6/Y/TjnV04b2o/znJdF2czPoHTTfdy4k99lULheBv+MFuGgdPY34IzDOjA\n+XBIH+VXDi5XtiCcAlwjAGMiAduy/z6F0xhfwjHynYXzpk7gzNffjzP/n2MXjoW/FLlxf4jCmYPt\nwK+zx2UYXSWo+MvDu6DVap3QdQlCJbhGAEq6Aj+B0xjz/fq3MPoWnoczNCie7//TbL73MTqVOMDo\nJ8CKZw7eiBMYlGE0gKgLR3yKhiEtbeX8jAWh+rhGAApcgbuA3+E0vCCjATq/wWm8+e67P8eZ2y/m\nMZxuvlmUfzOOgBTPHNgULgvWlv03t5jIEeBPoXlW8zivSBAqpxpfB16hlHpBKfV7pdQN1ajUqWDN\n1Wuc7wJ24byx34OzMtBi4FngkWzGYq+/ci9kHzCdsb7+78fxIsxnMY6d4Syc8+fWJDiIY4O4LHue\n8yQWQKgtFQmAUsoA7sKZ0Ppz4Cql1JuqUbFqs3LpSjb84wZmPDdjrAV/FU7wj8HYcN3cWn75bMF5\n85fz9S9l4c/1tfI/Jro4m3cz8GaJBRBqT6U9gPOAF7XWXVrrFE7QbPE7sWFYuXQlb/nzt5TeqXDe\n/rlAni6cN/VLOG/vnGV/D07Dt3DG/aUo5fkcxunuF4vPe3FmG/4As72zxQlIqCmVCsAZwCt5v1/N\npjUsJ/06UD9jP/t9Fc4yX34ce8CVOB/4LOfsc6wobReOh2C4TKXmOuVGEpFxX4cgVINKBaCxonzG\nwYgtIJ/8rwP5Kf2mXkGhp6DF6KrA+Wv8eXD6Rbm0+xhd/bfc3H92QZC5M+eWySAIp4ZKZwEO4UyW\n5ZiH0wsoYP369SPb7e3ttLe3V3jayZPrYn/2rs9y4NAB4qH4SAO1tlucNfMsuo52ESU69uC8N7g3\nmf2OQBuFvv7Z9f1GVvw18/YvBO9/eUlfkR7Nvw1nevBsOFOdWfkFCq6ls7OTzs7OCR1TUTiwUsoL\n/DfO+/I1nJntq7TWz+flabhw4BwdOzvY9KNNxDNxAp4Aq69azcqlK1n+seXsaNsx9oBsuO/CZxZy\nzcXXcPdjdxd8brx1dyvRnijJOUnnrX4WzD04l9NbT6eptYmAJ8AFb7yAbY9uGxWfs4A2p8wNn9og\nNgChaownHLji9QCUUu8G7sSxoX9Xa/3lov0NKwDl6NjZwXXfuK6gcVvbLc5qPYsz5545IhSlBAQo\nKSqlzjGefIIwWWoiAOOoxJQTAJAGKkx9RAAEwcXI14EFQTghIgCC4GJEAATBxYgACIKLEQEQBBcj\nAiAILkYEQBBcjAiAILgYEQBBcDEiAILgYkQABMHFiAAIgosRARAEFyMCIAguRgRAEFyMCIAguBgR\nAEFwMSIAguBiRAAEwcWIAAiCixEBEAQXIwIgCC5GBEAQXIwIgCC4GBEAQXAxIgCC4GJEAATBxUxa\nAJRSVyqlnlNK2Uqpt1ezUoIg1IZKegC/BT4A7KtSXWpGZ2dnvaswhkasEzRmvaRO1WPSAqC1fkFr\n/T/VrEytaMSH1Yh1gsasl9SpeogNQBBcjPdEO5VSO4G5JXat1VpvOzVVEgShViitdWUFKLUX+L9a\n62fK7K/sBIIgTBqttTrR/hP2ACZA2ZOcrAKCINSPSqYBP6CUegW4AOhQSv2ietUSBKEWVDwEEARh\n6lKTWYBGchpSSq1QSr2glPq9UuqGetYlW5/vKaW6lVK/rXddciil5iml9maf2e+UUmsaoE4BpdST\nSqlfK6UOKKW+XO865VBKGUqp/UqphjGMK6W6lFLPZuv1VLl8tZoGbAinIaWUAdwFrAD+HLhKKfWm\netYJ+H62Po1ECvgnrfWbcYZ4/1jv+6S1jgOXaa3PAd4GXKaUWlTPOuVxHXAAaKTutAbatdbnaq3P\nK5epJgLQQE5D5wEvaq27tNYp4MfA++pZIa31w0BfPetQjNb6iNb619ntCPA8cHp9awVa6+HspgkY\nQG8dqwOAUupM4ArgO5zAGF4nTloftzkCnQG8kvf71WyaUAalVBtwLvBkfWsCSimPUurXQDewV2t9\noN51Ar4OfAbI1LsiRWhgl1Lql0qpT5TLVK1pwKniNNRIXbSGRykVBn4KXJftCdQVrXUGOEcp1QI8\nqJRq11p31qs+Sqn3AEe11vuVUu31qkcZLtZaH1ZKzQJ2KqVeyPY2C6iaAGitl1arrFPIIWBe3u95\nOL0AoQillA/4GXC31npzveuTj9Z6QCnVAbwD6KxjVS4CVimlrgACQLNS6gda64/WsU4AaK0PZ//t\nUUrdjzP8HSMA9RgC1HOc9Evgfyml2pRSJvARYGsd69OQKKUU8F3ggNb6znrXB0ApNVMpNS27bQFL\ngf31rJPWeq3Wep7WegHwV8CeRmj8SqmgUqopux0CluEY4sdQq2nAhnAa0lqngU8BD+JYbX+itX6+\nHnXJoZT6EfAY8Aal1CtKqY/Vsz5ZLgauwbG078/+1Xum4jRgT9YG8CSwTWu9u851KqZRhphzgIfz\n7tUDWusdpTKKI5AguBi3zQIIgpCHCIAguBgRAEFwMSIAguBiRAAEwcWIAAiCixEBEAQXIwIgCC7m\n/wOqOI+0WlJ0lgAAAABJRU5ErkJggg==\n",
              "text": [
                "<matplotlib.figure.Figure at 0x7fd7ac0e4828>"
              ]
            }
          ],
          "prompt_number": 46
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "What does this mean?\n",
            "\n",
            "* We have now got a moral equivalent to a local expansion: The point values at the target skeleton points.\n",
            "* Is it a coincidence that the skeleton points sit at the boundary of the target region?\n",
            "* How many target proxies should we choose?\n",
            "* Can cheaply recompute potential at any target from those few points.\n",
            "* Have thus reduce LA-based evaluation cost to same as expansion-based cost.\n",
            "\n",
            "Can we come up with an equivalent of a multipole expansion?"
          ]
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "----------------\n",
            "\n",
            "Check that this works for 'our' sources:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "imat_error = (\n",
            "    P.T.dot(interaction_mat(targets[:, target_skeleton], sources))\n",
            "    -\n",
            "    interaction_mat(targets, sources))\n",
            "\n",
            "la.norm(imat_error, 2)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 49,
              "text": [
                "6.1074574232342807e-09"
              ]
            }
          ],
          "prompt_number": 49
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "## Source Skeletonization"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "nproxies = 25\n",
            "\n",
            "angles = np.linspace(0, 2*np.pi, nproxies)\n",
            "source_proxies = 0.5 + 1.5 * np.array([np.cos(angles), np.sin(angles)])"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 54
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "pt.plot(sources[0], sources[1], \"go\")\n",
            "pt.plot(targets[0], targets[1], \"ro\")\n",
            "pt.plot(source_proxies[0], source_proxies[1], \"bo\")\n",
            "\n",
            "pt.xlim([-1, 5])\n",
            "pt.ylim([-1, 5])\n",
            "\n",
            "pt.gca().set_aspect(\"equal\")"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "display_data",
              "png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAEACAYAAABccqhmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX90XVWZ9z9Pk+YH/QGlBcqvmdA4CFJFUBgQVwmy2lQi\nv/RFKLJQX+ugQtNR39GxPyCIeZ3BdzmSFsQRcAZrERmh2N6hTUqbBkQBpVhKy6iB8KuEQlsgTZu0\nSfb7x3NO7rk396ZJ7s099+Y8n7Xuyrnn7rP3Pvdmf/ezn/3sfcQ5h2EY0WRc2BUwDCM8TAAMI8KY\nABhGhDEBMIwIYwJgGBHGBMAwIkxxphmISBvwHtALHHTOnZ1pnoZh5IaMBQBwQJVzbncW8jIMI4dk\nawggWcrHMIwckg0BcMB6EfmDiHw5C/kZhpEjsjEEOM8594aIHAU0icgLzrnHspCvYRijTMYC4Jx7\nw/v7log8BJwN9AuAiNhiA8MICefcoMPzjIYAInKYiEzyjicAc4DnUlQir1433XRT6HUohDrla72s\nTkN7DYVMLYBjgIdExM/rF865xgzzNAwjR2QkAM65l4APZ6kuhmHkmEhGAlZVVYVdhQHkY50gP+tl\ndcoeMtSxwogLEHGjXYZhGAMREdxoOgENwyhsTAAMI8KYABhGhDEBMIwIYwJgGBHGBMAwIowJgGFE\nGBMAw4gwJgCGEWFMAAwjwpgAGEaEMQEwjAhjAmAYEcYEwDAijAmAYUQYEwDDiDAmAIYRYUwADCPC\nmAAYRoQxATCMCGMCYBgRxgTAMCKMCYBhRBgTAMOIMCYAhhFhTAAMI8JkRQBEpEhENovI6mzkZxhG\nbsiWBbAQ2AbYQwANo4DI6PHgACJyAnARUA98I+MaGcYQaInFaGxooLi7m57SUubU1jKrpiareWWz\njHwlYwEA/g34J2ByFvIyjEPSEouxbuFC6ltb+88t9o4P1XD9z/bu2MErL7/Mzn37mNbby4eBHuA4\n4PbHHuPnRx/Nezt3cv3+/cwKlLH16afZ8bvfUdzdzWvvvUcJcPTkyYUrEM65Eb+ATwG3e8dVwOoU\naZxhZJPFc+Y4BwNeS6qr3aY1a9yiysqE84sqK92mNWtSfjYf3CbveBO465LyvA7cP4Bb7H3+2fLy\n/rSLktL65eQLXtsbtA1nagF8DLhERC4CyoDJInKvc+7aYKK6urr+46qqKqqqqjIs1ogyxd3d/cct\nwL3AXqBv/Xqee+IJvtHRkZC+vrWVpcuW4ZxLsBoAfgosBWYBjcDVwBLUNO7x3jcBtwCLgSP276cF\nuB24P6lefjlhWQHNzc00NzcP65qMBMA5twhYBCAi5wP/J7nxQ6IAGEYmtMRiPPvMM3oM/CcwHbgL\noLcXOjr6HVGzAtcVdXWlzbPI+/sWsA51ZvksBt72juuBS700R5AoFHO88gYrZ7RJ7lxvvvnmQ16T\nDR9AEJsFMEYNf+w/raODxcAe4Fjge/7naC8+Ge2htwI70Ia947e/pbSsLGW+vd7fd4CfJH1WD1wV\neF8CHA/8BTgqcP4//bzSlJGvZE0AnHObgE3Zys8wkmlsaKC+tZU64GjgJbRBgjb+5N77K8CH0F7p\nJz09tOzdy+KkNNcBn/OOj01T7vTA8SnAGuAM4sIDain8S0kJ/7xgwXBvK1SybQEYxqjhj/170J59\nBtoTg/b89Unp7wSuJD5W94cES4FXgN3AJ1Efwh3A/jTlTvL+fh3oAMoBQUXHz7Me+F8lJQU3C2Ch\nwEbB0FNaCuh4+xXv7160903Xk5UnvZ+FOvROAk4HZgLHAL8EvunlFeQ6YCca6PIe6mv4Ndr7r0NF\nwGfC+PHDv6mQMQEwCoY5tbUsrqwEtOFv8M5XA9vTXLMLdda1JJ3vRQXkX9He2/cfvI1aDfNQS+Fz\nwFSgFLjbu7bFy3M86mvw85540kkjvrewsCGAUTDMqqlh69NPs/LWW3lwvxrsLcBK4HoYML5fhJr4\nW4FlQAPqve8DvuCdn0Rq/8F8tOdvAuYC/+OdT5X2q8D/LSrikosvztKd5g7ReIFRLEDEjXYZRjS4\no66OR+rr+UhPT8LUWws6hj8ceBc1aytRp95fgR8G8liMDh/ORRv3scCRJDr0fJaiw4X5qGWwCu35\n06Xtq6yk+rbb8sYPICI452SwNDYEMAqCO+rq2FJfz+qeHupIHIPPQr3zP0HH8sVow91BYuMH7bkr\ngEfQmYR24MU0Zf4ZHQq8CLwf+BLpTeYiNBCoadmyEdxdeJgAGAXBpuXLubOnJ+FcPdqLtwBPomb9\nlWhvvQR4LU1eRei4/lXgVHRYkIqTgfuAOtSy+DzpfQ1+LEGYgUAjwXwARkFQntT4fXaiPoBHAueu\nQ4cH95I6Wq/Xu+4w4g7AVPEBp3vHjeiUok8qX8Nc7ziygUCGMZrsL079r/oq8N9J536COgXLGBis\n8x9ob/4pYKN3PhgfUIQ6/L6KWheQ2Ej8tNejkYiVaOOfBXxJhM8XWCCQDQGMgqBizhy+nHRuPmrK\np2IPsDzpXD3wOrp7zddIjBHw4wPqgPd571/x3m8ncRpxFjr9d5j3fgMqHu0TJuSNA3ComAVgFATj\nd+3iDDQuvwzoQhviI2nSv4s23qDpDzpT4B/PR3ewCToKg+b833h5AP3i41/7D8C1gfefB2q++c3h\n3VQeYAJgFATF3d18De25fVqAGGqu/zhwfj7wbeKN04/umwXsQ/0CQVG4Hg0LnoCuLdiA9vDnB/L8\nKRoN+CBqGXR7aX6BziZ0jBvHzLPOyvQ2c44JgFEQ+GHAfsTeW2g8fsw7txR4GZ0B+GcSlwLXe5+v\nBb7lffZV1El4LbB/3Dha+/o4E40n8Fnsvd+BNhRB1wI8mJRmNnBLX1+oewGMFPMBGAXBnNpaLjns\nMFaijr2jiHvm/fH7veg/9KwU179C3FkHajEciy7jPerKK/mbiRMTPP2gwrHJK68OFZvpJPoD/KlI\nKLwpQDABMAqIkp6e/kaaznRNNwl3ALUcgo23CI3vf7mxkQlpZhlOTXpfj4pHcH2Bv6FIoU0Bgg0B\njAKhsaGBmQcO9L9PHRWg/9BXog3XdwCuRU3+WajJ7m8U8irakDt27+bAxIkp8+tNce79qEWwOJDm\n69Onc3mBTQGCCYBRIBR3dyc0+jkMDMj5LHACiV59f1MQ3/SvRgOHgub+V5zjQx0dXFdUxE96403+\ni94rGT9FPXA5OhVZqKtdTACMgqCntDSh0fsN+nLU7H8fauYnx/7fiToAfZKj+vw0VwLn9/Zy1dSp\nnDJzJr1lZRx95JH87L77EnwKwWlCvLLvAmhvL0gnoAmAURDMqa1lXWsr1a2t/RF7z6KWwA7UE//X\nNNcWBY7T/cOfigYJHT59OnWBnXU/8+STXP7ii0xG4wKCjkSAg8FyCtAJaAJgFAR+z9q0bBlFXV1s\n37qVObt28Qa6EOhNNCw3FcFx/J8GSVMPXNXennD+gyeeyGsvvshe4A0SG///RqcFr0RjBgrRCWj7\nARgFyXUzZ3LU88/zFjrF9/fo6r/JJA4D/H38TkAb+QE06OengTS+WT8LuHrCBE7+6EfZ+d57HAA6\n//xn9nV24sf4NXn5vYEGEPmC8GXgyHnz+NeVK7N/syNkKPsBmAVgFBwtsRjy4ovMQTfpCIYDfwlt\nmEcBm9F9/vxG6j/Q4wigBjgLFYWgWV/Z2UndJt3cejHqNMS77lhgIioiyQ8F+SlwVWNjdm4wh1gc\ngFFwNDY0cOf+/TQy0Ol3N7rDTx3qoAs2/nVow/0cGtBzEA0g8tMsQn0JPvVe+lXe3x+hQUHvpqmX\n7E+3r3D+YhaAUXD424MXEw8NDq75L0ItgT7v790kbhve6J3zQ4iDS4CTowj3oNOGJJWVvJ4AYG8g\nTqFQMAEwCoqWWIztW7dyHfpMgB3Avwc+/wZq+h+J/nNPQ510fYE0/j/9LOINeAnx/QX9Rr7Ty7/O\nO+4mvjMwJC4yWgRMPv74jO8v19gQwCgYWmIxVs2fz/W7djENfTLtvyel+SHwcXTv/vtRc/18YO+4\n+L96qijCOWgg0Trisf93AOcBn/COU60D+DEqMCcAM045JaP7CwMTAKNg+OXSpfywvb3fnB9sg06f\nO9FGO76oqL/H9gOKIL7H/11oD5/8dKHgYp/gsc/7UaFpBGZbKLBhjB6dL70ExP9p060HSI7fLwPc\nwYO8DVyG7uTjWwaVwD1euro0+RWlOQ6WVTR+fMFFAYJZAEYB0S06pe03/GBP7pPsyQfdO+B61Pt/\nGurUi3nH9wTSDUVQgsfBsnpFaInFDnkP+UZGAiAiZSLypIg8KyLbROT72aqYYSRzcMoUvkK84c9C\n5+mXoht7zEX38g965r+I7v3XiO4XEDTxO5PyP5SgzEcdj9d6ZfrxA18HvnHgQME9EwAyHAI457pE\n5ALn3D4RKQYeF5GPO+cez1L9DANQB+Df7tvHZeg4/G3gYvSZfcegjXMW8Bnv/ER0TP+PxAXhqyQ+\n0bc7qQz/vB8k9JaX5i7UuTgH3ZKsxavDPd7fy71rN0RxLYBzbp93WIIOkXZnmqdhJNPY0MAPvTj9\nYA9/PRql5/NBYAu6PDe4dReox35p4PqJDFxS/AsS9xMMluPvR+hPH/qPDvMpxLUAGQuAiIwDnkH9\nKT92zm3LuFaGkYQf/JPMUUnve9EIvw+lyeflwPG16Nbh/gYi/vg+1ZZie5Lez/eu91lUWcncKM4C\nOOf6gA+LyOHAOhGpcs41B9PU1dX1H1dVVVFVVZVpsUbE8DcFTSbZKdeOCkA6h57vEDzKu7YK3QW4\nzvt8SZrrphCPGtxeXk7FZZfRtHs3G7q66C0rY+6CBaHPAjQ3N9McWMo8FLK6GlBElgL7nXP/L3DO\nVgMaGdMSi7Fu4ULqW1v7z80vKWHvgQOUogt0JqC98t2VlRR3djK9vT3xkd8lJZx54ACvk2j216Cz\nApDmUeHjx/NGaSl/V1HBpOOPZ3YeNPahMJTVgBkJgIhMA3qcc++ISDn63d3snHs0kMYEwMgKLbFY\n/34AvWVlHHvOOby+YkWCKCyqrGTubbcBcO/SpXS2tVECTDzpJE67+GJeX7GC6tZWmkjszQ9u3Njv\nY2gBGkpKOOHkkwuqwSeTCwH4ILqz8jjv9XPn3A+S0pgAGKNGsigcqrGmSz/cfAqBUReAIVbCBMAw\nQmAoAmCRgIYRYUwADCPCmAAYRoQxATCMCGMCYBgRxgTAMCKMCYBhRBgTAMOIMCYAhhFhTAAMI8KY\nABhGhDEBMIwIYwJgGBHGBMAwIowJgGFEGBMAw4gwJgCGEWFMAAwjwpgAGEaEMQEwjAhjAmAYEcYE\nwDAijAmAYUQYEwDDiDAmAIYRYUwADCPCmAAYRoQxATCMCJORAIjIiSKyUUSeF5GtIlKbrYoZhjH6\nZPp48OnAdOfcsyIyEfgjcJlzbnsgzZh8OnAs1kJDQyPd3cWUlvZQWzuHmppZYVfLMPoZytOBizMp\nwDnXDrR7x3tFZDtwHLB90AsLnFishYUL19HaWt9/rrV1MYCJgFFQZGQBJGQkUgFsAk5zzu0NnM9b\nC2CkvXh19RIaG7+X4vxS1q69ZVTKNIzhMuoWQKCgicB/AQuDjd+nrq6u/7iqqoqqqqpsFJsRmfTi\n3d2pv7aurqJRK9MwDkVzczPNzc3Du8g5l9ELGA+sA/4xzecuH5kzZ7EDN+BVXb1k1K7NpEzDGC5e\n2xu0/WY6CyDA3cA259yPMskr14y0FweorZ1DZeXihHOVlYtYsGD2qJVpGKNBpkOA84BrgC0istk7\n9x3n3NoM8x11Skt7Up4vK+s95LW+ub5s2VK6uoooK+tlwYK5hzTjMynTMEaFQ5kImb7I0yHAmjWb\nXGXlogRTvLLyO27Nmk1jqkwjujCEIUDWZgHSkatZgJF412OxFpYtawr04rNH3Rk30jJt9sAYLkOZ\nBRgTFkDqnnXRmOlZx/r9GaMDo+0EzBcaGhoTptYAWlvrWbasKaQaZZexfn9GeIwJARjr3vWxfn9G\neIwJARjr3vWxfn9GeIwJARjpvHyhMNbvzwiPMTULkGuPfi4Z6/dnZJ+hzAKMGQEwDCORoQjAmBgC\nGIYxMrKyGnA0sMCXkWPfnTFU8lIAbNnsyLHvzhgOeTkEsMCXkWPfnTEc8lIALPBl5Nh3ZwyHvBQA\nC3wZOfbdGcMhLwXAAl9Gjn13xnDI2zgAC3wZOfbdGWCBQIYRaSwQyDCMQTEBMIwIYwJgGBEmlEhA\nC1UNH/sNDAhBACxUNXzsNzB8cj4EsFDV8LHfwPDJuQBYqGr42G9g+ORcACxUNXzsNzB8ci4AFqoa\nPvYbGD6hRAJaqGr42G8w9rFQYMOIMDkJBRaRe0TkTRF5LtO8DMPILdnwAfwMmJuFfAzDyDEZBwI5\n5x4TkYrMqzL6xJpiNKxsoNt1Uyql1F5dS83smgFpli5byks7X0J6hYqjK7jk/Ev43fbfDXrdcMoI\npnv9zddpf6udY489luOmHjdo3oaRbfJyU9DRINYUY+HtC2k9o7X/XOvteuw3uFhTjPk/mE/78e3w\nDjAO9rywh807NsNlDLgOSGjs5556LiueWDFoGQPKaQcmwq73drH1na1s+d4W7uIuEwEjJ2TFCehZ\nAKudcx9M8Zm76aab+t9XVVVRVVWVcZnDpfqL1TRWNEIb0IoOfjph4v6JfOSjH6FUSnlr11tsnrpZ\nP7/Qu/DRwHGA8l+V0zuxlwMXHYife6ic/ZfvH5D2zGfO5OLzLmb5A8vpKeph7zt76T23F/4ETEzK\n/1E4Y9IZPBN7Jjs3bkSG5uZmmpub+9/ffPPNuZkFOJQA5MMsQNUXqtgkmxIbN8DDQBEwAeRNwZU7\nuDTw+UbgghQZ/hI4ChWSPqASeClN2l95aYLp9wL7gasGJp+yZgq7n949rPszjGSGMguQkyFAdfWS\n0FedlUrpwMbfRkIP7HCwBngK6EAb69tpMnQM6LkZ2PkrvcDRw0gfYkSurRKMFhkLgIjcB5wPTBWR\nV4EbnXM/C6ZpbPxe/3FYq85qr66l5TstdNEVP5ksCACfAh4ArvDetwGrgYsDaVYD5yZddyFqTawi\nwV/Aw0B3inIuBFakrutJx5w0yJ2MHrZKMHpkPA3onJvnnDvOOVfqnDsxufEnE9aqs5rZNZx63KmJ\nJ9Pd/bTAcQXwQbSxPghsQHv/Cu/zNrQ33wh0ohZDDFgJ3A8cAEq8dMkUAWsTT01/fDrfveG7h7yf\n0cBWCUaPUHYECmvV2S0LbqFyc2X8RF+ahMkuiwpgqnf8CWCCd9xG3Iq4ALgaHefvBI4BrkQtiau9\ndG1J+ZYCpwAbYMLqCVS/XM1d3wpvBsBWCUaPUAQgrFVnNbNruOZj1zA1NpXDHzmc8p3lFK9O+qdf\nBcxIcXFQLCahwwTfUd8W+OxSYDypTf4XA+/Xoo7A36NTjgdhwbwFoU7/2SrB6JHzOABddRZO4GCs\nKcaKJ1awq2ZX/7nxvx5P+a/KOTD+AL2H92qjbCVu4gOsB3pQk/1h1HF4ReDzR72//jUlaSrwLjpU\neBv4WxKEopNO5v9gfqgxALW1c2htXZwwDAjz9zJGn5wsBqquXpIXq876YwGSeQDt4Q+ijfLv0d5a\ngH1o4z+IevMnoI7CZDagwwNQf8E1adK8h8rufhJFxK/jy9WsvWftwA9yhK0SHDvkzTTg2rW35KKY\nQ9LtulN/MI34/P0DaE9eQeIYH7SnTzdo8r/m9aiIxIBgR74Kbfyl3usg8BtUUCrptx66+gKzFCFQ\nUzPLGnyEiEwoMHixAKkIGkFnEZ/2S54mrASeTpP5G6h4nIWKQSUqAgeALjSu4GgSpwgf9dL5kcMV\nUDaubIh3YxiZE6nnAtReXUv52vLEk+tJdPpVoILwAMieJOupAh27P5x4Wh4W+Bja+FuJN+rTgMnA\ndOAwEhs/qLg8T7+DcPrj01kwb8GI7s0wRkLkNgSZedFMnn/7eQ3O6UUbbUVSog0w6b1J0Acdn+4Y\nmMnD6EzAuzBJJtHZ2Unf57xpgjbUf/AaOhSYgJY1DvhMigo9CJwJ8pSw+kerbRGQkTXs2YApOP6Y\n47Vn9j31vydxGm8t8Bp0HOigo6hDx+lBVqGzAe8A+6BjTwd90qei0IaKyQxgChotOAF19h2RpkJH\nAC9CkSuyxm/knEj5AACmlkyFx9DxeDlqrv8Jip8opoce/Uau9RK3AZuA/0AFw9fSc4lbDcFx/J+8\nc62oDyG4krASXWcQnEFYD7wPeAZKx6XxTxjGKBIpAYg1xVi1ZZVG6Pk8CpwO5c+U03F4R7zBtqEN\n+XwGOgOD8/6+Y3Aa6vB7HIrKiuilN9G+qkAFYgMqJA5t/BV6zTFHH5Ot2zSMIZOTIUB19RJisZZc\nFDUoDSsb2D83aQme54Dr6O3QEN5HSZz+S7VgyI/q89NdgU4jXgFMht69XuRccqjx6ahIXIDGDFSg\nQ46Pw74J+4g1xTK+x2xRV3cH06ZdyRFHfIFp066kru6OsKtkjAI5cQKCo7JyMbfdVh3qHPPMi2by\n/N7nE9fwV6DRef7X8AlUBA4AnyT9fgC/QZ1705LyAp3+KyM+NAhaFU951/hWQAcqDBW6ccgfH/5j\nNm41I+rq7qC+fgs9PXf2nysu/gqLF3+IurqvhVgzYzjklRMw7FVlsaYYf9n1l/jCHb93bwNeAd4E\ndhMf07/jXZhqwVAb2oD9nj+YF6h47EWn+PagkYGPoEOFz6Ii41sBl9K/RmDb69vywgpYvnxTQuMH\n6Om5k+XLw7fijOyS01mAMFeVNaxsSNi+C9CG24z21tNRz73fa7+HBvZ0ooFBQX5P4v4Afl7+Yp/D\n0YZdRnwNgT/yaEtROU+juyZ0sey+ZUO/qVGip6c8zXkLUhpr5NQJGOaqspRhwG1oow9uAeZbAO3E\nzftJ6BZgh6Fm/+FpChF0TH+K974StQSS84fE2ANH/4xA2KHAAMXFqbcqKi4Ov25GdsmZBRD2s+dS\nhgG3ktg4Id6TH0fcvH8DXQHo61e6rbzeRAViKPn7PAzsoH9GIB9CgW+44XyKi7+ScK64+DpuuMHW\nCIw1crQn4FIWLJgbqgOw9upatvxgC+3ntcdP7kmT2HfQwYB9AwHd5HMtiY9DWY/GFVyKTvVVoL6A\nVLwL/BodIoCGC1dA+dpyFnw3/FBgdfTdwfLlV9HTU0ZxcRc33DDLHIBjkEiFAseaYty4/Eae/euz\n9PX2qWmfamnvSnRoUIL29qnS+OHAvljMIL4r8EbgJOB3wLwU1/4SOIf4MGAllB9WzrfmfYu6b9WN\n8O4MIxF7OGgaJp4+kc5Pdw5c7gvasL1pOUBj9T+dIpNU04P+ngAPoMt9P5Yi/9XoHoMV8VMTfjOB\n+//lfgsFNrJK3uwHkG9093kD9QrvhB+d9zrwcRIddEXE9wEIzvfvTMrUD+tdjy4wejJF/o6BG4q2\nQsn4EhpWNgCYCBg5JZICUFZUxl726psK4g3yfhIbfxu6eUdyGPAfUWffA6gwdKPj/W50oRHEfQjB\n/PGu8fN+AZgLe9hDI41s+YE9FszILZFbDRhriunc/KNJH6wn0YMPar4nb4d3IbrBx/FoINBn0F1/\npxF3Kj6NjvGTy1iF7iewAY0lSMq7/bx2blx+47DuxzAyIVIWQN2tddy68lb2l+6HXSSa5v6TuIIP\nAUknj8WkfqDIBijZVkLfwT56KrwddoNl7EWnFIMrC5N46c2Xhn1fhjFSIiMAsaYYtz50a/zhnW1o\nT12KxuOf5yX8PWqm+xt3piJdPJOg0YY/9/IoQhv7THQYEEM3Gd1D+p2DbQduI4dE5tmAA1YCVnh/\nn0Ib40uok28G2lN3o/P1DwGXBzJaj3r4U+GP+yeQOHOwFnjWu66P+C5ByU8eXg9TyqcM676yjT0b\nMFrkRADy4dmAKUOBf482xmBc/8PEe+ET0aFB8hr+93npLiU+lfgu8UeAJc8cnIIuDOoj8ZmDTzFg\nGHJ4Rbo449HHng0YPXI+BNBVgUtz/g+VEArcBmxFG95hxBfo/AltvMHw3QeJ7/cf5AnUzC9JSr8K\nFZDkmYNeErcFq/D++puJtAPvg8lHTR7iHWWf9M8GzP3vZeSGjGcBRGSuiLwgIn8RkW8P5ZowVgXW\nXl2rzwVsQ3vsT6E7A10IbAEe9xImR/2l65DHA0cyMNb/MjSKMMiFqJ9hBokPE21FfRAXeOWcHe5a\nAHs2YPTISABEpAhYjk5ofQCYJyKnDn5VOKsCa2bXcNv1tzH1+akDPfiXoIt/ihi4XNffyy/Iw2jP\nny7WP5WH329bwYeJXuilXQWc5q0FCHFbcHs2YPTI1AI4G/irc67NOXcQjXJP7hMTCHNVYM3sGmZ+\nYGbqDwXt/f2Vem1oT/0S2nv7nv0NaMMvR8f9qUgV+TwRNfeTxedidLbhRTi6+OhQg4Bqa+dQWbk4\n4VzYqziN0SVTH8DxwKuB96+hT9ZLoLp6aeBZc+GuCjzk04HeIfUagbVow/f9ASvRWP9kT/4qBloG\n69EIwe1pKjUduAD2xvYO4Q5GD/93WbYsf34vY3TJVACGtMonX54NCOoLaL29ldYzWuMn/Th+0DH5\n0wx8cOdc4mG8oBZAhXcc9OSPQ+0i/9zbxB8+EigyAW9DkOnTpo/gjrKLPRswWmQqAK+jk2U+J6JW\nQAJ1dXX9x1VVVVRVVWVY7MjxTewbl9/Itte30TWhq38zjvK15cyYNoO2nW100jnw4onxw+ID3nME\nKkiM9ff3AqhArYaSwOeVUPzfxfRcFBhrr0anB0+HE+SEzG/QiCzNzc00NzcP65qMlgOLSDHwP6gR\nvAOd2Z7nnNseSJN3y4F9Yk0xlt23jK6+LsrGlbFg3gJqZtekf4y4t9y38plKrjnvGlY8sSLBkpjy\n6BQ63+rkwDEH+vcImN46neOmHMekKZMoG1fGOaecw+rfro6LzwygQvO87YbbbCGQkTVysh+AiHwS\n+BHqQ7/bOff9pM/zVgDSEWuKsfD2hQmNu3xtOTOmzOCE6Sf0C0UqAQFSikqqMoaSzjBGim0IkgHW\nQI1CxwTMTjknAAAE5ElEQVTAMCJMXj0YxDCM/COU5cC24ix87DcwIAQBsBVn4WO/geGT8yFA+hVn\n4T03MGrYb2D45FwAbMVZ+NhvYPjkXABsxVn42G9g+ORcAGzFWfjYb2D4hBIHEIu1sGxZU2DF2Wxz\nPuUY+w3GPhYIZBgRxgKBDMMYFBMAw4gwJgCGEWHy9slAFqo6cuy7M4ZKXgqAhaqOHPvujOGQl0MA\nC1UdOfbdGcMhLwXAQlVHjn13xnDISwGwUNWRY9+dMRzyUgAsVHXk2HdnDIe8jQS0UNWRY9+dARYK\nbBiRxkKBDcMYlLyMAxgJYz34ZazfnxEOY0IAxnrwy1i/PyM8xsQQYKwHv4z1+zPCY0wIwFgPfhnr\n92eEx5gQgLEe/DLW788IjzEhAGM9+GWs358RHiOOAxCRK4A64BTgLOfcM2nS5SQOYCTBL2F41kda\npgX3GMNlKHEAOOdG9EIb/snARuDMQdK5fGPjxo1uzZpNrrJykQPX/6qsXOTWrNk0auUOVubGjRtH\nrdxMyMd6WZ2Ghtf2Bm3HIx4COOdecM79eaTXh0lzc3PGnvVYrIXq6iVUVdVRXb2EWKzlkNcMVmZz\nc/OQ659L8rFeVqfsMSbiAEZCJp71kc7LmzffyDcGtQBEpElEnkvxujhXFRwtMvGsj9R6MG++kW9k\nvBhIRDYC33SDOAEzKsAwjBHjDuEEzNYQIG0hh6qAYRjhMWInoIhcLiKvAucAMRF5JHvVMgwjF4z6\nfgCGYeQvOYkEFJErROR5EekVkTNzUeYgdZkrIi+IyF9E5Nth1sWrzz0i8qaIPBd2XXxE5EQR2ej9\nZltFpDYP6lQmIk+KyLMisk1Evh92nXxEpEhENovI6rDr4iMibSKyxavXU+nS5SoU+DngcuDQk+Wj\niIgUAcuBucAHgHkicmqYdQJ+5tUnnzgIfN05dxo6xLs+7O/JOdcFXOCc+zDwIeACEfl4mHUKsBDY\nBuSTOe2AKufcGc65s9MlyokA5FHQ0NnAX51zbc65g8AvgUvDrJBz7jFgT5h1SMY51+6ce9Y73gts\nB44Lt1bgnNvnHZYARcDuEKsDgIicAFwE3MUgzvCQOGR9xsRioGFwPPBq4P1r3jkjDSJSAZwBPBlu\nTUBExonIs8CbwEbn3Law6wT8G/BPQF/YFUnCAetF5A8i8uV0ibIWCSgiTcD0FB8tcs7ly9gon0y0\nvEdEJgL/BSz0LIFQcc71AR8WkcOBdSJS5ZxrDqs+IvIpYKdzbrOIVIVVjzSc55x7Q0SOAppE5AXP\n2kwgawLgnCuEtamvAycG3p+IWgFGEiIyHvg1sMI5tyrs+gRxzr0rIjHgo0BziFX5GHCJiFwElAGT\nReRe59y1IdYJAOfcG97ft0TkIXT4O0AAwhgChDlO+gPwdyJSISIlwJXAb0KsT14iIgLcDWxzzv0o\n7PoAiMg0ETnCOy4HZgObw6yTc26Rc+5E59xJwFXAhnxo/CJymIhM8o4nAHNQR/wAcjUNmBdBQ865\nHuAGYB3qtb3fObc9jLr4iMh9wBPAySLyqoh8Mcz6eJwHXIN62jd7r7BnKo4FNng+gCeB1c65R0Ou\nUzL5MsQ8Bngs8F2tcc41pkpogUCGEWGiNgtgGEYAEwDDiDAmAIYRYUwADCPCmAAYRoQxATCMCGMC\nYBgRxgTAMCLM/wd3SfE+A0GaywAAAABJRU5ErkJggg==\n",
              "text": [
                "<matplotlib.figure.Figure at 0x7fd7ac867b00>"
              ]
            }
          ],
          "prompt_number": 56
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Construct the interaction matrix *from* the sources *to* the source proxies as `source_proxy_mat`:\n",
            "\n",
            "**A note on terminology:** The `source_proxies` are *near* the sources but *stand in* for far-away targets."
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "source_proxy_mat = interaction_mat(source_proxies, sources)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 73
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "source_proxy_mat.shape"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 74,
              "text": [
                "(25, 200)"
              ]
            }
          ],
          "prompt_number": 74
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Now compute an ID (row or column?):"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "idx, proj = sli.interp_decomp(source_proxy_mat, nproxies)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 75
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "source_skeleton = idx[:nproxies]"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 76
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "P = np.hstack([np.eye(nproxies), proj])[:,np.argsort(idx)]\n",
            "tsm_approx = source_proxy_mat[:, source_skeleton].dot(P)\n",
            "\n",
            "la.norm(tsm_approx - source_proxy_mat, 2)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 78,
              "text": [
                "6.3367197566876269e-15"
              ]
            }
          ],
          "prompt_number": 78
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Plot the chosen proxies:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "pt.plot(sources[0], sources[1], \"go\", alpha=0.05)\n",
            "pt.plot(targets[0], targets[1], \"ro\")\n",
            "pt.plot(sources[0, source_skeleton], sources[1, source_skeleton], \"go\")\n",
            "pt.plot(source_proxies[0], source_proxies[1], \"bo\")\n",
            "\n",
            "pt.xlim([-1, 5])\n",
            "pt.ylim([-1, 5])\n",
            "\n",
            "pt.gca().set_aspect(\"equal\")"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "display_data",
              "png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAEACAYAAABccqhmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYXEW99z91Tu/dM5NlspAECAwoviCil02WMBiZCQQS\nUDFBr6JXFAQSBOUuCXkZr+a9710elAkCiuK9oCxeF26SkSySTAZEEa6iLFEuEwgh+zKTmd67T9f9\no7on3bNkmaW7Z/r3eZ5+5nR1nao650x9T9WvflWltNYIglCZWKUugCAIpUMEQBAqGBEAQahgRAAE\noYIRARCECkYEQBAqGNdQE1BKvQ10AQ6Q0lqfO9Q0BUEoDkMWAEAD9VrrA8OQliAIRWS4ugBqmNIR\nBKGIDIcAaOBXSqmXlFJfHIb0BEEoEsPRBbhQa71TKTUJWK+U+rPW+tlhSFcQhBFmyAKgtd6Z/btX\nKfUL4FygRwCUUjLZQBBKhNb6sN3zIXUBlFIBpVRV9jgINACv9FOIsvrcfffdJS/DaChTuZZLynR0\nn6NhqC2AKcAvlFK5tH6stV43xDQFQSgSQxIArfVbwFnDVBZBEIpMRXoC1tfXl7oIfSjHMkF5lkvK\nNHyoo+0rDDoDpfRI5yEIQl+UUuiRNAIKgjC6EQEQhApGBEAQKhgRAEGoYEQABKGCEQEQhApGBEAQ\nKhgRAEGoYEQABKGCEQEQhApGBEAQKhgRAEGoYEQABKGCEQEQhApGBEAQKhgRAEGoYEQABKGCEQEQ\nhApGBEAQKhgRAEGoYEQABKGCEQEQhApGBEAQKhgRAEGoYEQABKGCEQEQhApmWARAKWUrpf6glFo1\nHOkJglAchqsFcBvwOiCbAArCKGJI24MDKKVmAFcAy4E7hlwiQTgK2lpaWNfcjCuRIO310rB4MbPm\nzh3WtIYzj3JlyAIAfAu4E6gehrQE4Yi0tbSw9rbbWN7e3hO2NHt8pIqb+y28YwfvbN3KnmiUWsfh\nLCANTAO+8+yzPDp5Ml179nBLLMasvDxeffFFdvzmN7gSCd7t6sIDTK6uHr0CobUe9Ae4EvhO9rge\nWNVPHC0Iw8nShgatoc/nrsZGvWn1ar2krq4gfEldnd60enW/v90AelP2eBPoG3uleSPoL4Femv39\nk35/T9wlveLm8ikXsnXvsHV4qC2AC4B5SqkrAB9QrZR6RGv92fxITU1NPcf19fXU19cPMVuhknEl\nEj3HbcAjQBjI/OpXvPL889zR3V0Qf3l7O8tWrEBrXdBqAHgIWAbMAtYBnwLuwjSN09nv64FvAEuB\ncbEYbcB3gCd7lSuXT6laAa2trbS2th7TOUMSAK31EmAJgFLqEuBrvSs/FAqAIAyFtpYWXv79780x\n8B/AVOD7AI4D3d09hqhZeefZ8fiAadrZv3uBtRhjVo6lwL7s8XJgfjbOOAqFoiGb3+HyGWl6v1y/\n/vWvH/Gc4bAB5COjAMKIkev713Z3sxToAI4Dvpn7HfMWr8a8oV8FdmAq9o5f/xqvz9dvuk72byfw\n3V6/LQcW5n33ANOB/wEm5YX/Ry6tAfIoV4ZNALTWm4BNw5WeIPRmXXMzy9vbaQImA29hKiSYyt/7\n7X0TcCbmrfTddJq2cJilveLcCHw6e3zcAPlOzTs+DVgNfJBDwgOmpfD/PR7+ftGiY72skjLcLQBB\nGDFyff805s1+MuZNDObNv7xX/AeBBRzqq+e6BMuAd4ADwOUYG8L9QGyAfKuyf28HugE/oDCik0tz\nOfAJj2fUjQKIK7Awakh7vYDpb7+T/RvGvH0HepP5e32fhTHonQR8ADgDmAI8AXw1m1Y+NwJ7MI4u\nXRhbw88wb/+1GBHIEXS7j/2iSowIgDBqaFi8mKV1dYCp+Buy4Y3A5gHO2Y8x1rX1CncwAvLPmLd3\nzn6wD9NquA7TUvg0MBHwAj/IntuWTdONsTXk0g6ddNKgr61USBdAGDXMmjuXV198kcf+5V/4ecw0\n2NuAx4BboE//fgmmif8qsAJoxljvM8DnsuFV9G8/uAHz5l8PzAH+kg3vL+6Xgf9n28y76qphutLi\noYy/wAhmoJQe6TyEyuD+piaeXr6cv0qnC4be2jB9+BrgIKZZW4cx6r0J3JOXxlJM9+HDmMp9HDCB\nQoNejmWY7sINmJbBU5g3/0BxM3V1NN57b9nYAZRSaK3V4eJIF0AYFdzf1MSfli9nVTpNE4V98FkY\n6/x3MX15F6bi7qCw8oN5c88EnsaMJOwCtgyQ5xuYrsAW4L3AFxi4yWxjHIHWr1gxiKsrHSIAwqhg\n03338WA6XRC2HPMWbwNewDTrF2De1ncB7w6Qlo3p128D3ofpFvTHe4DHgSZMy+J6BrY15HwJSukI\nNBjEBiCMCvy9Kn+OPRgbwNN5YTdiugeP0L+3npM9L8AhA2B//gEfyB6vwwwp5ujP1jAne1yxjkCC\nMJLEXP3/q24Dftkr7LsYo6CPvs46/455m18JbMyG5/sH2BiD35cxrQsorCS5uLdgPBHrMJV/FvAF\npbh+lDkCSRdAGBXMbGjgi73CbsA05fujA7ivV9hyYDtm9ZqbKfQRyPkHNAGnZL+/k/2+mcJhxFmY\n4b9A9vsGjHjsCgbLxgB4tEgLQBgVuPfv54MYv3wfEMdUxKcHiH8QU3nzm/5gRgpyxzdgVrDJNxTm\nN+dPyKYB9IhP7twvAZ/N+349MPerXz22iyoDRACEUYErkeBmzJs7RxvQgmmuP5AXfgPwdxyqnDnv\nvllAFGMXyBeFWzBuwUHM3IINmDf8JXlpPoTxBvw5pmWQyMb5MWY0oduyOOOcc4Z6mUVHBEAYFeTc\ngHMee3sx/vgt2bBlwFbMCMDfUzgVeHn29zXA32Z/+zLGSPhZIGZZtGcyfAjjT5Bjafb7DkxFUZi5\nAD/vFecy4BuZTEnXAhgsYgMQRgUNixczLxDgMYxhbxKHLPO5/vsjmH/oWf2c/w6HjHVgWgzHYabx\nTlqwgBNCoQJLPxjh2JTNrwkjNlMptAfkhiJh9A0BggiAMIrwpNM9lXSgputAg3BJTMshv/LaGP/+\nrevWERxglOF9vb4vx4hH/vyC3IIio20IEKQLIIwS1jU3c0Yy2fO9f68A8w+9AFNxcwbANZgm/yxM\nkz23UMg2TEXuPnCAZCjUb3pOP2HvxbQIlubFuX3qVK4ZZUOAIAIgjBJciURBpW+gr0POJ4EZFFr1\nc4uC5Jr+jRjHofzm/k1ac2Z3NzfaNt91DlX5z2c/vcnFWA5cgxmKHK2zXUQAhFFB2ustqPS5Cn0N\nptl/CqaZ39v3/0GMATBHb6++XJwFwCWOw8KJEzntjDNwfD4mT5jADx9/vMCmkD9MSDbv7wPs2jUq\njYAiAMKooGHxYta2t9PY3t7jsfcypiWwA2OJf3OAc+2844H+4d+HcRKqmTqVpryVdT/+wgtcs2UL\n1Ri/gHxDIkAqP59RaAQUARBGBbk36/oVK7DjcTa/+ioN+/ezEzMRaDfGLbc/8vvxfzxMnOXAwl27\nCsLff/zxvLtlC2FgJ4WV/28ww4ILMD4Do9EIKOsBCKOSG884g0mvvcZezBDfeZjZf9UUdgNy6/jN\nwFTyJMbp56G8OLlm/SzgU8Eg7zn7bPZ0dZEEIm+8QTQSIefjtz6b3k6MA1FOEL4ITLjuOv75sceG\n/2IHydGsByAtAGHU0dbSgtqyhQbMIh357sBfwFTMScAfMOv85SppbkOPccBc4ByMKOQ36+siEZo2\nmcWtl2KMhmTPOw4IYUSk96YgDwEL160bngssIuIHIIw61jU382Asxjr6Gv1+gFnhpwljoMuv/Gsx\nFffTGIeeFMaBKBdnCcaWkGN5Nv5T2b/fxjgFHRygXCo20LrC5Yu0AIRRR255cBeHXIPz5/zbmJZA\nJvv3BxQuG74uG5ZzIc6fAtzbi7ADM2xIr7x6zycACOf5KYwWRACEUUVbSwubX32VGzF7AuwAvpf3\n+x2Ypv8EzD93LcZIl8mLk/unn8WhCnwXh9YXzFXyPdn0m7LHCQ6tDAyFk4yWANXTpw/5+oqNdAGE\nUUNbSwtP3XADt+zfTy1mZ9rv9YpzD3ARZu3+JzHN9UuAsHXoX70/L8IGjCPRWg75/t8PXAh8JHvc\n3zyABzACMwM4+bTThnR9pUAEQBg1PLFsGffs2tXTnD/cAp05HsRUWrdt97yxcw5FcGiN/+9j3vC9\ndxfKn+yTf5zjvRihWQdcJq7AgjByRN56Czj0TzvQfIDe/vs+QKdS7AOuxqzkk2sZ1AEPZ+M1DZCe\nPcBxfl622z3qvABBWgDCKCKhzJB2ruLnv8lz9Lbkg1k74BaM9f90jFGvJXv8cF68oxGU/OP8vByl\naGtpOeI1lBtDEgCllE8p9YJS6mWl1OtKqX8aroIJQm9S48dzE4cq/izMOP0yzMIeczBr+edb5j+P\nWftvHWa9gPwmfqRX+kcSlBswhsfPZvPM+Q/cDtyRTI66PQFgiF0ArXVcKXWp1jqqlHIBzymlLtJa\nPzdM5RMEwBgAT4xGuRrTD98HXIXZs28KpnLOAj6eDQ9h+vRf4ZAgfJnCHX0TvfLIheechPZm43wf\nY1xswCxJ1pYtw8PZv9dkz91QiXMBtNbR7KEH00U6MNQ0BaE365qbuSfrp5//hr8F46WX4/3AnzDT\nc/OX7gJjsV+Wd36IvlOKf0zheoL5+eTWI8wNH+a2DssxGucCDFkAlFIW8HuMPeUBrfXrQy6VIPQi\n5/zTm0m9vjsYD78zB0hna97xZzFLh+cWEMn17/tbUqyj1/cbsufnWFJXx5xKHAXQWmeAs5RSNcBa\npVS91ro1P05TU1PPcX19PfX19UPNVqgwcouC9qa3UW4XRgAGMujlDIKTsufWY1YBbsr+ftcA543n\nkNfgZr+fmVdfzfoDB9gQj+P4fMxZtKjkowCtra205k1lPhqGdTagUmoZENNa/1temMwGFIZMW0sL\na2+7jeXt7T1hN3g8hJNJvJgJOkHMW/kHdXW4IhGm7tpVuOW3x8OHkkm2U9jsn4sZFYABtgp3u9np\n9XLqzJlUTZ/OZWVQ2Y+Go5kNOCQBUErVAmmtdadSyo+5d1/XWj+TF0cEQBgW2lpaetYDcHw+jjv/\nfLb/6EcForCkro45994LwCPLlhF5+208QOikkzj9qqvY/qMf0djeznoK3+apjRt7bAxtQLPHw4z3\nvGdUVfjeFEMA3o9ZWdnKfh7VWv9rrzgiAMKI0VsUjlRZB4p/rOmMBkZcAI6yECIAglACjkYAxBNQ\nECoYEQBBqGBEAAShghEBEIQKRgRAECoYEQBBqGBEAAShghEBEIQKRgRAECoYEQBBqGBEAAShghEB\nEIQKRgRAECoYEQBBqGBEAAShghEBEIQKRgRAECoYEQBBqGBEAAShghEBEIQKRgRAECoYEQBBqGBE\nAAShghEBEIQKRgRAECoYEQBBqGBEAAShghEBEIQKZkgCoJQ6Xim1USn1mlLqVaXU4uEqmCAII89Q\ntwefCkzVWr+slAoB/w1crbXenBdnTO4O3NLSRnPzOhIJF15vmsWLG5g7d1apiyUIPRzN7sCuoWSg\ntd4F7Moeh5VSm4FpwObDnjjKaWlp47bb1tLevrwnrL19KYCIgDCqGFILoCAhpWYCm4DTtdbhvPCy\nbQEM9i3e2HgX69Z9s5/wZaxZ840RyVMQjpURbwHkZRQCfgrcll/5czQ1NfUc19fXU19fPxzZDomh\nvMUTif5vWzxuj1iegnAkWltbaW1tPbaTtNZD+gBuYC3wlQF+1+VIQ8NSDbrPp7HxrhE7dyh5CsKx\nkq17h62/Qx0FUMAPgNe11t8eSlrFZrBvcYDFixuoq1taEFZXt4RFiy4bsTwFYSQYahfgQuCvgT8p\npf6QDfsHrfWaIaY74ni96X7DfT7niOfmmusrViwjHrfx+RwWLZpzxGb8UPIUhBHhSE2EoX4o0y7A\n6tWbdF3dkoKmeF3dP+jVqzeNqTyFyoWj6AIM2yjAQBRrFGAw1vWWljZWrFif9xa/bMSNcYPNU0YP\nhGPlaEYBxkQLoP8365Ix82Yd69cnjAyMtBGwXGhuXlcwtAbQ3r6cFSvWl6hEw8tYvz6hdIwJARjr\n1vWxfn1C6RgTAjDWretj/fqE0jEmBGCw4/KjhbF+fULpGFOjAMW26BeTsX59wvBzNKMAY0YABEEo\n5GgEYEx0AQRBGBzDMhtwJBDHl8Ej9044WspSAGTa7OCReyccC2XZBRDHl8Ej9044FspSAMTxZfDI\nvROOhbIUAHF8GTxy74RjoSwFQBxfBo/cO+FYKFs/AHF8GTxy7wQQRyBBqGjEEUgQhMMiAiAIFYwI\ngCBUMCXxBBRX1dIjz0CAEgiAuKqWHnkGQo6idwHEVbX0yDMQchRdAMRVtfTIMxByFF0AxFW19Mgz\nEHIUXQDEVbX0yDMQcpTEE1BcVUuPPIOxj7gCC0IFUxRXYKXUw0qp3UqpV4aaliAIxWU4bAA/BOYM\nQzqCIBSZITsCaa2fVUrNHHpRRp6Va1bS/HgziUwCr+Vl8XWLmTdnXp849z52LzEnhkd5uPmTN+Nx\nebjvyfsOe96x5AHgOA7xZJyW9S088JMHSKkUPst32LQFYbgpy0VBR4KVa1bylfu/wlt/9VZP2Jb7\ntwD0VLj+4rz67VexEhZ7G/b2OQ8oqOznn3Y+j/3mscPmAabyH4wcZNX6VXztnq+xz95n2mIZeGX5\nK33iC8JIMSxGwGwLYJXW+v39/Kbvvvvunu/19fXU19cPOc9j5aPXf5RnTn6mT3j9m/Ws+t4qFIp5\nN85jQ92GvidvAD5SGHT686cTdofZes7WnjDfGh/x0+IwszDu7C2zuXXBrdz72L0kdAI7Y3P9x67n\nnu/dw+bEZpidF/kZODN4Jn98+o+DvlahMmltbaW1tbXn+9e//vXijAIcSQDKYRTg4s9czHOnPNcn\n/OzNZ/PUd57Csi3m3zSfF9/3Yt+TNwKXFgZV/VcV3fO7+8btRyxOaTuFuDfOu+e92xM2/YXpHNh2\ngNgnYn2SqFlVQ+dLnUdxVYIwMEczClCULkBj410ln3Xmtbz9hruVm5gTw2t7sRnAFbYf/crYmf7j\n9nO793TsoWt+V0HY9vO2Y70zgA22hB65MkuwshiyACilHgcuASYqpbYB/1dr/cP8OOvWfbPnuFSz\nzhZft5g3v/MmW88+1GSf8cIMPvM3nyHtpEkn0iy4cgFvPfoWuy/Y3RNn3MZx6ITmIAd7wqb9dhre\nGi9v8Ra98XZ4SZDo+T79henYE2y66OoT16M8xIn3CZ85eeZgL3NIyCzByqMojkC9X6GNjctYs+Yb\nI5pvf/xk5U944KcPEM/EcePmuquuY/bs2cQiMfxBP+FYmDVr1/CTX/6ElErh0i4+duXHsFIWT294\nmgQJAq4AV116FVEryoP/+SC7zz8kFlOen8I1H7qGV995lUg6gnIUcz8yl5+v/Tmvfei1PuU5eePJ\nHLQOsv+S/T1hk5+dzENfe6gkRsDGxrsKxPpQeGmelzA0yqYL0JtSzTr7+NyP0zi7Ectl4TgOuzt3\nk0ln8HpM90CnNQ2XNXDpJZey/+B+MjqD3+dnavVUPv+JzxOOh1GWYm/HXvbG93J9/HrWbFiDg4MH\nD9d+8lrmN8zHSTvEiLGzYyeOy+Fg8iC7n97Nvgv29ZRlyvNTuPn6m/F7/Tze8jhpnSbkCbHoa4tK\nNgIgswQrj5IIQKlmndm2zcZNG3us8W7cXD/vei655BJ0RjM+OJ6tO7eivApvyIvX5SXoCWLZFgpF\nwB1gX3gfHrcHT9zDmWedyQfO/gDjfePxe/wE3UECgQDbd24nkokQjUeJJqJccNEFpHWa9c+sx7Ec\n/Laf+fPm42QcHl35KCmdIuQKsWhh6So/yCzBSqToAmBmnZXGcXDlmpXc8eAdZpz+baAdXviXFzjl\n30/h9s/dzkUXX8TEmol0JjrJJDNkyOANerGUMdbZto3P7UP7NLjAFXYRT8axHIuACjCxaiLxVJwM\nGRzbwV/lx1EOGZ3hrLPP4oz3n0HIF8LlcvHyiy9zz3/ew+4PH+pCvHP/O0DpfAAWL26gvX1pgQ2g\nlM9LGHmKYgNobLyrLGad9fgCvA20UzD+PuOFGdy54E4uuuAifH4fGk08HcfSFn6vHzIQ6Y7g4ODY\nDi7bRTKZxBfw4VIuqt3VdMY6wQP79u4jZafQLk0sFiNNGsdxcJIOOq1JZVLcc+89bP7Q5r5lfOuj\nrP/30q3MI7MExw5lYwMoFwNSIpO1zveq/ADvnvcuT6x9ggtnXYhWmupgNUEnSOfBThzt0NHRQcbO\nEE/HqQnUYCmLKk8VtmXjsl1EU1EyZIh2R3HbbrSjsV02/qCfnft3knbShA+GsfwWGZUhoRN9ygcQ\nd/qOChSTuXNnSYWvICrGFRjyfAEGGH5PZBIEAgFi0RgJdwKtNf6An+5wN/5qPxkrQ/RglM5IJxOr\nJ5IhgztjKnsilcC2bWon1JJIJkimk1gZC6UUdtCmM9pJwpPAV+PjQOcBXLr/W++zfSN09YLQl4ra\nF2DxdYs58aUTYQAfHq/yYtu2afI7kElmsByLUDCEZVvYts240DhIQSKeAAeqglVYlsWkCZOYWDMR\nNHg9Jh3LsrCxmTxlMuND45k2YRpWytzyCz58AbXP1xbkP/OlmSxauGikb4Mg9FBRLYB5c+YRT8b5\nx/v/kTeffpPE5Yea4TNemMENX7gBK2NhKYugKwguiDgREokEGTIoS+HxeKitqcVWNiFPCNuxURlF\nPBVHofC7/aTSKVwZF/v27yMYCHKg+wDxZJxEOkE0GiVUE+LC2RcSqgrxzPpnsF0247zjuPOWO2US\nkFBUKm5FoEgswupnVvONFd9g275taFszbcI0bv/87VzReAXd3d1kkhls2yYYCBJOhFEexcHwQSwr\nOxqgbDKpDF5lRghSOkVapQmFQiit8Ll97OvYh9fj5UD4AFbAorurm65EF1u3b2X8lPG4lRuP7aHa\nXc20cdMY7x7PydNPLvHdEcYSZWMELCfWbljLnQ/cybYLtvWERX8Xxa3dhDvDYIM74Ea5Fd3JblKJ\nFPFwnHQmTTqZxrZtPB4PPo+PJEmUS+F1e0mEE0TCEUKBEN1d3QS8ARLpBMHqINFkFJ/fR0dXB36f\nn/D+MCdMPYGAO8DEmolorfF5pO8vFJ+KEgDHcWh+vJlt524rCN927jYebXmUK+ZcQcbKEIlHSDtp\n4uk4uEFlFJbbIpVIUVNVg9ty4/a4yTgZkk6Szv2deFweotEoAVeAgCeA2+Pm4IGDONohk86QSCYI\n+ALYlo1CEXQF8bl8pFNpbG3j9/hLdFeESqYoRsDGxrtoaWkrRlaHJZ6MF0zUyacz1sn2vdvZ27GX\nVDJFIplAuRSxeIwECXwBH4HqAMqtSOgEkXiEVCrF/u79OC4HfOCp8tAV7yIRSxCOh4nEjN9ASqdQ\nbkW1v5oZtTOodldjWRZpncZreZlQMwGtNI5TPh53TU33U1u7gHHjPkdt7QKamu4vdZGEEaAoLYB1\n675ZFrPKWta3sHnzZji1729el5c4cfxuP7FoDBxwBVwcCB/A7/dDxAzRacz4fldHF53xTroz3djK\nZrI9GRcuEk6CuIrjtb2EgiHe3f0uylI4yiFgBYiEI/jdfrTW/PbXv2XVr1ah3Aq/7efma2/mk/M+\nWfwb04umpvtZvvxPpNNP9oQtX34TcD9NTTeXrmDCsFPU2YClnFW2cs1KFt+3mK2TtvZxBJr6/FQW\nXbuIM888E7fPbTwBoxrLYxFJRAiGggR8ATKpDDqt0Wje2fkOjsdBuRUel4d0PE21q5oJ4yeQjCTx\n+r1YLotYLMYbb7+B7bfxur3MmDKDrlgX//27/+aBnz7AzvN39pTjxJdOpPmW5pKPBNTWLmD//if7\nhE+cuJB9+54oQYmEwVCUZcGPhVLOKmt+vNks3zUTqMOs3LMRQk+FWNiwkJlnzKQr2QUK4rE4nQc6\nIQEBO0A8GjfNcwVaa7oj3UyunUzADuByXASsAJPGTQIFLp+LuBMnaSXpjnVzMHqQ6onVhGpCKJei\nM9JJLBnjiZVPFFR+gK1nb2XFEytKcHcKSaf7t0ek02KoHGsU1QhYylllPW7AYERgpjmc+uJUTj3n\nVMJOGJ/Xx96uvUwMTiTqRJnom4ilLKy0Rcf+DsYHxuPz+bCx8QQ8KKXYfXA3B8MHCflDpHWaWHcM\nv8+Pg0Ncx+lMdWJZxq9ApRUd4Q4sZeFY/d+LUrsCA7hcfZcpM+GlL5swvBStBVDqvecGWhLM4/IQ\n6Y6QJk1CJ3BcDu072rGrbbRbE6gKEPQFCQVDhGNhLG3hsT10RDrAC26fm1g6xo59O0jH0uiUxuf3\nEfAESCaTPX4BHreHoD+I3+cnlUnhUZ5+y1MOrsC33noJLtdNBWEu143ceqvMERhrFGlNwGUsWjSn\npAbAxdctZsv9WwqW7J78/GQa5zZSHaomlUihlLH614RqqAnVEA1HIQPdsW4yTgY7aePgcKDzAFEr\niuW2yNgZqqurzapCfj9JnSTRYTwHI90R3B43LlxEuiNU+atwMg7RriiXX3o5e3+5l13n7+opz4kv\nnsiiW0vvCmwMffdz330LSad9uFxxbr11lhgAxyAV5QmY27SjK9FF555OkiQJjA+gHMUll17Ch8/7\nMG7lpmNfBxMmTMDGxuvzkrJSZlQg7WJS7SR27d+Fy+Vif8d+vH4vXreXqqoqUrEUyqWIRqMcN/E4\ntu/ZTsbKEPAF8Ckfu/fvBg/UeGpweVz8/sXfs3L9StKkqfZUc9t1t3H1FVeX+jYJYwTZHHQAHvnp\nIyz94VLePffQMt2Tnp/EwtkLObnuZCZPmkxXogtv0EsqmmL8hPGQgumTppNMJemOdpPRZkZR2klT\nO7EWy7KId8XRliYcDVPjryFUFSKaiJJKpfApH5FUhEwmw+RJk01rIxbDj59p46YxZcIUbFuW3hKG\nDxGAAbj4Uxfz3Hv77hFw6m9O5e6ld+PyueiOd4Myi4BY2qK2uhaf34fX9uL1eNnfsR+3x004Hqa6\nppqQK0QzXIRfAAAMCUlEQVSVv4p4Ko7P9pGIJbC8Zu1BK2NBBpRX4bJdeD1eNJpNmzbxk1U/ARsC\nroBsCyYMKzIXYADS9L/2nfIoguODRLojeCwPPp8PUpBOp/GP8xOLxYikIoScEAFPgLSVJh6J41Zu\naifW4sZtKrPfGAD9QTOcZmUsIl0RaoI1JBLGJfiZjc/wb4/9G9vP3d6Tf3/biAnCSFJR6wGAmQ/g\nVf2PCPgsY4G33TY+vw8n7RAIBQgFQvgtP9GuKH6fn0QqgeWxGD9uPKfVnUbQGyRDhkwqwwT/BGKx\nGONqxpFOpAmHw6SSKXxeH10dXfjdfqKRKA8/+XBB5Qd466/eovnx5hG/B4KQo6IE4KlfPsXs62ez\na+8uvE8XisC030zjE1d+glQyBbaZN5BOm3H98Z7xBF1BTpp0ElV2FR7lwe1247N9ONphSu0Uqqqq\ncLldeG0vk3yTsNM28YhZMNTGxuPykElnsCzLLCGWifZbxmi6/3BBGAkqpguwcs1KvvLAV8zOQKcC\nb4P3515mTJ3BuNA4Fly3gHPOO4dEKkF3optwJIxLuThh0glUV1UTi8WoClRhuS1UShHwBHDSDqFA\nCBcuAr4Alm0RCAR459132BfbZyYA+arxB/xEY1HGTzTrAETj0QH9AAbyVxCEkaBi9gZsfry5YFsw\nZkJiZoLazbXcfefdaFuDArflJmgHGTdxHK6MywwLdnQwfvx4s/hHKoWFhaUt0pk0kUiEkD9EOBym\n2l3Ntn3b2BXdRSQTIa3T7N+3n1p/LUFvEJ/tI+2kicajzG+Yz85f7CxYFnzab6dx04039S18EZG9\nASuLos0GzFGqWYEFrsB5RNNRbJcNLkgmknhtL37lJ1gVxOPy4LJc+DI+otEotm2WAaueXE13dzfx\nZBzlNsN5XpeX7Z3bSbvTZFwZ4tE4KVJknAxux42tbWzLdAUsn8V5F59HzImxbuM6HBy8lpdPL/g0\nl3/08qLel3xkb8DKo+g2gPb25axYUfx17wdqWr+95W3++Kc/4na7SWfSZkGQVARtaTx+D9Gk6ZOH\nqkI4GYeacTVUhapwu9yMC47DiZtJQqFQiKQryYGDBwjHw/ir/Hi8HvxVfuPfb0E8Hsfr8VIdrOY3\nm37D2g1re9yCP/aRj/Gpqz+FZZfOLNPcvK6g8kPpnpdQHIZjd+A5wLcxm1p/X2v9z0c6pxSzAvtz\nBeYp6HZ1s/RbSzl+wvEs/NhC6i+tN9N7U8aCfyByAMtlUROoweV10RUxMwYtt4VlW0ycNJF4NI7K\nKEhhdgIigxNzcHvd/O6537Fm3RpQUOOt4TPzPkPCSfDQfz3Erg8fcgPes3IPJ08/mStnX1n0e5ND\n9gasPIYkAEopG7gP+CiwHXhRKbVSa913y5s8SjErMDe2/rm7PkdHqAO6AB8wB1Kk2MIWHl71MCF/\niIZZDTi2Q1ekC5/fRzQcJWpHcWXMgh9u3FQFqlizYQ2Pr36cpE7iUz6uabiGfeF9/HLjL3Esh0hH\nhIPRg3Rdfmhr8G0/3kbACRRUfjDLkj345IN84vJPFPGuFCJ7A1YeQ20BnAu8qbV+G0Ap9QQwHxhQ\nAEq519y8OfM4/cen89wpz8Ez9NkdaNf5u3hqzVPMb5jPmg1reGzlY6RVGjtjc+0V1zL74tkk0gm0\nS/Pcr5/jX3/8r+w4f0fP+W98/w20V9NR33Eo0WcwW5HNNF+3nbsN/3/1P99+d8fukroDy96AlcdQ\nBWA6kL/C5rvAeb0jNTYuy9trrrSzAo+0O5CDQ9tzbdzzxD0Fjjq7f7abmmANH7nkI8QTcR7+2cMF\nlR/ggOsA1PdKcDZm8ZGZh4ISsf4Nkrt37+43vFjknsuKFeXzvISRZagCcFRO/uWyNyDk2QIyb/X7\ne8gT4pFfPNLHS2/7udt59KlHuWL2FWY1X9VPc3kg+10vb2yv8hJ7JlbYAvkVTJk45egvZISQvQEr\ni6EKwHbg+Lzvx2NaAQU0NTX1HNfX11NfXz/EbAdPzhaw7N5l/OXpvxTsDnTiiydy85du5ltPfKvf\nc9OkqQnUAOBT/SzcMcCWY/kyecKLJ1B1fBWvTXvNtAxU9vdTYLqefuwXJAhZWltbaW1tPaZzhjQb\nUCnlAv6CeZftAH4HXJdvBCzH2YA5Vq5ZyYonVhB3zAy+RQsXMW/OvEPbiPeivr2elQ+tJJPOsHHT\nRu548I6CUYVJ6yehfIo9F+/pCZv87GSmuqdSXVuNz/ZxyydvIZlO8rcP/W2BY9JJL53Et2/5tkwE\nEoaNEZ8NqLVOK6VuBdZihgF/cKQRgHJi3px5/Va4/oYMT3zxRL78pS9jZSwCvgBXX3E1lmUVCshd\nZjWfgrCvLeo3D5/HVxjvlv7jCcJIUpHrARwNA7UOBGG0IAuCCEIFU3b7AgiCUF6UZDqwzDgrPfIM\nBCiBAMiMs9Ijz0DIUfQugMw4Kz3yDIQcRRcAmXFWeuQZCDmKLgAy46z0yDMQchRdABYvbqCubmlB\nWKn3Daw05BkIOUriB9DS0saKFevzZpxdJsanIiPPYOwjjkCCUMGII5AgCIdFBEAQKhgRAEGoYMp2\nZyBxVR08cu+Eo6UsBUBcVQeP3DvhWCjLLoC4qg4euXfCsVCWAiCuqoNH7p1wLJSlAIir6uCReycc\nC2UpAOKqOnjk3gnHQtl6Aoqr6uCReyeAuAILQkUjrsCCIByWsvQDGAxj3fllrF+fUBrGhACMdeeX\nsX59QukYE12Ase78MtavTygdY0IAxrrzy1i/PqF0jAkBGOvOL2P9+oTSMSYEYKw7v4z16xNKx6D9\nAJRS1wJNwGnAOVrr3w8Qryh+AINxfimFZX2weYpzj3CsHI0fAFrrQX0wFf89wEbgQ4eJp8uNjRs3\n6tWrN+m6uiUadM+nrm6JXr1604jle7g8N27cOGL5DoVyLJeU6ejI1r3D1uNBdwG01n/WWr8x2PNL\nSWtr65At6y0tbTQ23kV9fRONjXfR0tJ2xHMOl2dra+tRl7+YlGO5pEzDx5jwAxgMQ7GsD3ZcXqz5\nQrlx2BaAUmq9UuqVfj5XFauAI8VQLOuDbT2INV8oN4Y8GUgptRH4qj6MEXBIGQiCMGj0EYyAw9UF\nGDCTIxVAEITSMWgjoFLqGqXUNuB8oEUp9fTwFUsQhGIw4usBCIJQvhTFE1Apda1S6jWllKOU+lAx\n8jxMWeYopf6slPofpdTflbIs2fI8rJTarZR6pdRlyaGUOl4ptTH7zF5VSi0ugzL5lFIvKKVeVkq9\nrpT6p1KXKYdSylZK/UEptarUZcmhlHpbKfWnbLl+N1C8YrkCvwJcAxx5sHwEUUrZwH3AHOD/ANcp\npd5XyjIBP8yWp5xIAbdrrU/HdPFuKfV90lrHgUu11mcBZwKXKqUuKmWZ8rgNeB0op+a0Buq11h/U\nWp87UKSiCEAZOQ2dC7yptX5ba50CngDml7JAWutngY5SlqE3WutdWuuXs8dhYDMwrbSlAq11NHvo\nAWzgQAmLA4BSagZwBfB9DmMMLxFHLM+YmAx0DEwHtuV9fzcbJgyAUmom8EHghdKWBJRSllLqZWA3\nsFFr/XqpywR8C7gTyJS6IL3QwK+UUi8ppb44UKRh8wRUSq0Hpvbz0xKtdbn0jcqpiVb2KKVCwE+B\n27ItgZKitc4AZymlaoC1Sql6rXVrqcqjlLoS2KO1/oNSqr5U5RiAC7XWO5VSk4D1Sqk/Z1ubBQyb\nAGitR8Pc1O3A8Xnfj8e0AoReKKXcwM+AH2mtnyp1efLRWh9USrUAZwOtJSzKBcA8pdQVgA+oVko9\norX+bAnLBIDWemf2716l1C8w3d8+AlCKLkAp+0kvAacqpWYqpTzAAmBlCctTliilFPAD4HWt9bdL\nXR4ApVStUmpc9tgPXAb8oZRl0lov0Vofr7U+CVgIbCiHyq+UCiilqrLHQaABY4jvQ7GGAcvCaUhr\nnQZuBdZirLZPaq03l6IsOZRSjwPPA+9RSm1TSn2+lOXJciHw1xhL+x+yn1KPVBwHbMjaAF4AVmmt\nnylxmXpTLl3MKcCzefdqtdZ6XX8RxRFIECqYShsFEAQhDxEAQahgRAAEoYIRARCECkYEQBAqGBEA\nQahgRAAEoYIRARCECuZ/AQEB96VcovUCAAAAAElFTkSuQmCC\n",
              "text": [
                "<matplotlib.figure.Figure at 0x7fd7ac113e80>"
              ]
            }
          ],
          "prompt_number": 79
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Check that it works for 'our' targets:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "imat_error = (\n",
            "    interaction_mat(targets, sources[:, source_skeleton]).dot(P)\n",
            "    -\n",
            "    interaction_mat(targets, sources))\n",
            "\n",
            "la.norm(imat_error, 2)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 81,
              "text": [
                "1.5507578112941745e-08"
              ]
            }
          ],
          "prompt_number": 81
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "* Sensibly, this is just the transpose of the target skeletonization process.\n",
            "    * For a given point cluster, the same skeleton can serve for target and source skeletonization!\n",
            "* Computationally, starting from your original charges $x$, you accumulate 'new' charges $Px$ at the skeleton points and then *only* compute the interaction from the source skeleton to the targets."
          ]
        }
      ],
      "metadata": {}
    }
  ]
}