{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Skeletonization using Proxies"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "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"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEF9JREFUeJzt3W9oXNeZBvDn0UiuJdvrrSW128TxKNCQbEiIS4eQxV9SOxTHTVu2sJCu3ISmIFLtLg60lBZBFwfEsv1QnA+tsyJ1u9vMNoRtQrdxtq3jpoRCSTpKnTpeN2AS2/W6YFlu/ihO/Ed698PM+I88Gs3onrnnXr/PD4Q0M3fOebE1j84999x7aWYQEb+6YhcgInEpBEScUwiIOKcQEHFOISDinEJAxLnuEI2QPAzgHQCzAM6bWSlEuyLSeUFCoOYTZnYyYHsikgLtDog4xxArBkm+AeDPAAzAv5nZRINtRgCMAMCKFSs+ftNNNyXuV0Qam5ycPGlmg61sGyoErjGz4yQ/BGAPgH8ysxcW2r5UKlmlUkncr4g0RnKy1bm5ILsDZna89v0EgKcB3B6iXRHpvMQhQHIFyVX1nwF8EsCrSdsVkXSEODrwYQBPk6y3959m9rMA7YpIChKHgJm9DuC2ALWISAQ6RCjinEJAxDmFgIhzCgER5xQCIs4pBEScUwiIOKcQEHFOISDinEJAxDmFgIhzCgER5xQCIs4pBEScUwiIOKcQEHFOISDinEJAxDmFgIhzCgER5xQCIs4pBEScUwiIOKcQEHFOISDinEJAxDmFgIhzwUKAZIHk70g+E6pNEem8kCOBbQAOBmxPpDXlMjAwAJDVr4GB6nNXS38dFiQESK4F8CkAj4VoT6Rl5TLwwAPA9PTF56angS9+sfkHs1wGhoaArq7q99HRi48HBoBVqy5+yFeurD5Xf+2++67s74EHchsENLPkjZD/BeBfAKwC8FUzu6fZ9qVSySqVSuJ+xblyGbj/fmB2duFtCoXq68UisGUL8OyzwJEj1Q93gN/9yxSLwOHDYdtcIpKTZlZqZdvEIwGS9wA4YWaTi2w3QrJCsjI1NZW0W/GuXAZGRpoHAHDx9SNHgJ07q9+B8AEAAEePhm8zBSF2BzYA+AzJwwCeALCR5OPzNzKzCTMrmVlpcHAwQLfi2rZtwOnTsau43Lp1sStYksQhYGbfMLO1ZjYE4F4AvzSzrYkrE2mkPil36T55VszM5HJeoDt2ASItq+8CZG0EUDc9Xa0PAIaH49bShiATg+3SxKAsydDQxX36LMvABGGqE4MiqcnLxFte6qxRCEh+rFkTu4LW5GyCUCEgEtr4eOwK2qIQkPw4dSp2Ba3J0aQgoBCQPMnZMDsvFAKSH3kYZhcKsStom0JA8iMPw+w774xdQdsUApIv/f2xK2ju0KHYFbRNISD58sgj1VN6sypnawQAhYDkzfAw8MEPxq5iYTmcvFQISP5k+VBhHiYv51EISP7k8K9tlikEJH/Gx4G+vthVNDY2FruCtikEJH+Gh4GJierZelmjiUGRlAwPV0/XXbEidiWXy+GuikJA8qtcBt59t/339fcDy5aFrweoXsw0Z3RlIcmn+lWGlqKTlyZ79tnOtd0hGglIPo2NZfMyY5oTEElJVi8zpjkBEcf6+rRYSMStYrF62DIPZzrOo4lByaf67cWy4MtfBr773dhVLJlGApJPSz0y0AmPPprLm47UKQQkn7L0l9esemPUnAaBQkDyKWsfuNnZ6ugka3W1QCEg+ZTFE3VOn85mXYtQCEg+ZXVRTlbrakIhIPmU1UU5Wa2ricQhQHI5yZdIvkLyAMntIQoTaSqL1xRwvFjoDICNZnYbgPUANpO8I0C7IgurX1MgK1cfLhRyu1gocQhY1UztYU/tK/37nYs/w8PAyZPA44+3Hwb9/WGvRTA3l8sAAALNCZAskNwH4ASAPWb2YoNtRkhWSFampqZCdCtSVQ8DsrXti8Xq9jMz1QAJsVuRw7mAuiAhYGazZrYewFoAt5O8pcE2E2ZWMrPS4OBgiG5FLtfKB3H+fnuz3YqentYuPpLTuYC6oEcHzOxNAL8CsDlkuyItWWiysH6zkoVO8rl0t6JYrI4oikXg+98Hdu26+Fx//8WwqN9zMMcnDtXRLNnuO8lBAOfM7E2SvQB+AeBfzeyZhd5TKpWsUqkk6lekoXK5umDn6NHqyGB8PNcf0KUiOWlmpVa2DXEW4UcA/DvJAqojiyebBYBIRw0Pu/zQJ5E4BMzs9wA+FqAWEYlAKwZFnFMIiDinEBBxTiEg4pxCQMQ5hYCIcwoBEecUAiLOKQREnFMIiDinEBBxTiEg4pxCQMQ5hYCIcwoBEecUAiLOKQREnFMIiDinEBBxTiEg4pxCQMQ5hYCIcwoBEecUAiLOKQREnFMIiDinEBBxLnEIkLyO5PMkD5I8QHJbiMJEJB0h7kp8HsBXzOxlkqsATJLcY2b/G6BtEemwxCMBM/uTmb1c+/kdAAcBXJu0XRFJR9A5AZJDqN6m/MUGr42QrJCsTE1NhexWRBIIFgIkVwL4MYCHzOzt+a+b2YSZlcysNDg4GKpbEUkoSAiQ7EE1AMpm9lSINkUkHSGODhDA9wAcNLNvJy9JRNIUYiSwAcAXAGwkua/2tSVAuyKSgsSHCM3s1wAYoBYRiUArBkWcUwiIOKcQEHFOISDinEJAxDmFgIhzCgER5xQCIs4pBEScUwiIOKcQEHFOISDinEJAxDmFgIhzCgER5xQCIs4pBEScUwiIOKcQEHFOISDinEJAxDmFgIhzCgER5xQCIs4pBEScUwiIOKcQEHEu1K3Jd5E8QfLVEO3FMrp7FN0Pd4Pbie6HuzG6ezR2SSIdF2ok8AMAmwO1FcVd/3EXdlZ2YtZmAQCzNoudlZ0Ng6C8v4yhHUPo2t6FgW8NYOBbA+ja3oWhHUMo7y8vGCYLPX9pe/U2RNJCMwvTEDkE4Bkzu2WxbUulklUqlSD9JjW6exSPVh6FofG/Q4EFnP/m+QuPy/vLGPnpCE6fO91WP9esvAbHZ45f8fzywnLM2izOzZ278FxfTx8mPj2B4VuH2+pDpI7kpJmVWtk2tTkBkiMkKyQrU1NTaXXb1OjuUeys7FwwAIDqiKD+l5rbia1PbW07AAA0DAAAeH/2/csCAABOnzuNsb1jbfchshSuRwLdD3dfGP4vhCB6e3qX9MFPgiDm/nku1T7l6pHJkUAWLRYAAPCBwgdSDwAAWLd6Xep9ik+uQ6DAwqLbvD/7fgqVXGnm7IwmCCUVoQ4R/gjAbwDcSPIYyS+FaLfTRj4+EruEBU2/N42Rn44oCKTjgs0JtCMrcwIAwO2MXUJTxdVFHH7ocOwyJGc0J3AVOfrW0dglyFXOfQj09/bHLqEpTRBKp7kOgfL+MmbOzsQuY0EEMb5pPHYZcpVzGwL1lX9nZs/ELmVBD5Ye1KpB6Ti3ITC2dyzK8f9WdbELG9ZtiF2GOOA2BLI+4TZnczpEKKlwGwJ5mHDTOQSSBrchML5pHH09fbHLWFTWRyySf25DYPjWYUx8egLF1UUArS0hXgwRfuFRHkYskm/dsQuIqT7zvpTrA8xHsOkpyUvR19OnQ4TScW5HAnWhjhIYLMhooq6/t18XFpFUuB4JAGH3uVs5NbkVm67fhOfuey5IWyKLcT8SyOI+99439upag5Ia9yHw0TUfjV1CQ0feOqJ1ApIK1yEwunsUe9/YG7uMBWmdgKTBdQhMTE6k2t9SDiFqnYB0musQCDWR16offu6HF9YltCqLcxZydXF/dCBN9cN9W5/a2tL2WicgaXA9EkhTuyMArROQtLgeCfT39mP6vemO99PT1YMtN2zB0I4hHHnrSEvvWblspQJAUuE6BNJiMDz28mNX3GmoGU0ISlpch8Cp906l0s/5ufOLbzSPJgQlLa7nBLL6QdOEoKTJdQhk9ZoCvd29sUsQR1yHwPCtw7j/tvtjl3EF3X1I0uQ6BADgyQNPxi6hIS0ZlrS4D4E0DhEulY4QSBpC3ZB0M8nXSB4i+fUQbQqwpndN7BLEgcQhQLIA4DsA7gZwM4DPk7w5abtpWdGzInYJC3r7zNuaF5COCzESuB3AITN73czOAngCwGcDtNtx5f1lnJ09G7uMBZ2bO6d5Aem4ECFwLYA/XvL4WO25zBvbO9bWKr4YNC8gnRYiBBqdJH/FZXdJjpCskKxMTU0F6Da5pB+wUJcYL64uLniCUVYXNMnVI0QIHANw3SWP1wI4Pn8jM5sws5KZlQYHBwN0m9xSPmD1KwoXVxfxYOnBxIuNlhWWYXzTeMOFS1o5KGkIce7AbwHcQPJ6AP8H4F4Afx+g3Y4b3zR+xT0Herp6sKywDO+ee/fCc/29/Xjk7kcantW3Yd0GjO0dw9G3jmLd6nUXPrSNntv2P9suOyTZqN3579OZhNJpNEt+wwySWwDsAFAAsMvMmv75KpVKVqlUEvcbQnl/WR88ueqQnDSzUkvbhgiBdmUpBESuRu2EgPsVgyLeKQREnFMIiDinEBBxTiEg4pxCQMQ5hYCIcwoBEecUAiLOKQREnFMIiDinEBBxTiEg4pxCQMQ5hYCIcwoBEecUAiLOKQREnFMIiDinEBBxTiEg4pxCQMQ5hYCIcwoBEecUAiLOKQREnFMIiDinEBBxLlEIkPw7kgdIzpFs6eaHIpItSUcCrwL4HIAXAtQiIhF0J3mzmR0EAJJhqhGR1CUKgXaQHAEwUnt4huSrafXdggEAJ2MXcYms1QNkrybV09yNrW64aAiQfA7AXzV4aczMftJqR2Y2AWCi1mbFzDIzh6B6Fpe1mlRPcyQrrW67aAiY2V3JyhGRLNMhQhHnkh4i/FuSxwD8DYDdJH/e4lsnkvTbAapncVmrSfU013I9NLNOFiIiGafdARHnFAIizkULgawsOSa5meRrJA+R/HqsOmq17CJ5IitrKEheR/J5kgdr/1fbMlDTcpIvkXylVtP2DNRUIPk7ks/ErgUASB4muZ/kvlYOFcYcCURfckyyAOA7AO4GcDOAz5O8OVY9AH4AYHPE/uc7D+ArZvbXAO4A8A+R/30A4AyAjWZ2G4D1ADaTvCNyTdsAHIxcw3yfMLP1raxdiBYCZnbQzF6L1X/N7QAOmdnrZnYWwBMAPhurGDN7AcCpWP3PZ2Z/MrOXaz+/g+ov+rWRazIzm6k97Kl9RZvdJrkWwKcAPBarhqS8zwlcC+CPlzw+hsi/5FlFcgjAxwC8GLeSC8PvfQBOANhjZjFr2gHgawDmItYwnwH4BcnJ2nL9pjp67kCoJccd1OjMJx0znYfkSgA/BvCQmb0dux4zmwWwnuRfAnia5C1mlvo8Csl7AJwws0mSd6bdfxMbzOw4yQ8B2EPyD7VRZkMdDYEcLDk+BuC6Sx6vBXA8Ui2ZRLIH1QAom9lTseu5lJm9SfJXqM6jxJhM3QDgMyS3AFgO4C9IPm5mWyPUcoGZHa99P0HyaVR3excMAe+7A78FcAPJ60kuA3AvgP+OXFNmsHqO+PcAHDSzb8euBwBIDtZGACDZC+AuAH+IUYuZfcPM1prZEKq/O7+MHQAkV5BcVf8ZwCexSEDGPES41CXHwZjZeQD/CODnqE56PWlmB9Kuo47kjwD8BsCNJI+R/FKsWmo2APgCgI21w037an/1YvoIgOdJ/h7VEN9jZpk4NJcRHwbwa5KvAHgJwG4z+1mzN2jZsIhz3ncHRNxTCIg4pxAQcU4hIOKcQkDEOYWAiHMKARHn/h/Ywd6NIp4rQQAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "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\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {},
      "outputs": [],
      "source": [
        "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)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {},
      "outputs": [],
      "source": [
        "def numerical_rank(A, eps):\n",
        "    _, sigma, _ = la.svd(A)\n",
        "    return np.sum(sigma >= eps)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check the interaction rank:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "9"
            ]
          },
          "execution_count": 29,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "numerical_rank(interaction_mat(targets, sources), eps)"
      ]
    },
    {
      "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",
      "execution_count": 30,
      "metadata": {},
      "outputs": [],
      "source": [
        "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)])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAE3hJREFUeJzt3X2MHOV9B/Dv914c7mzXjc+XNMHcHhIISqEQ3QlRuX8QjCLjkERFqkR6LihEOsG1kiOlapOe1MpIVtX8EcEfxPRESNrcKgiloDY2DQGHiCJFIXsEYqhBsuDOcR3Jh11ejMEvd7/+Mbu+t32ZnZ2Z59l9vh9pdbd7szO/nb397swzzzNDM4OIhKvLdQEi4pZCQCRwCgGRwCkERAKnEBAJnEJAJHA9acyE5CyA9wEsALhgZqNpzFdEspdKCJR91szeTnF+IpID7Q6IBI5p9Bgk+RaA/wNgAP7FzKaqTDMOYBwA1q9fP3L11Ve3vFwJ26lTwNwcsLi49FhXF1AoAJs3137ezEztv42MpL88F2ZmZt42s8E406a1O7DNzI6T/ASAZ0i+bmbPL5+gHAxTADA6OmqlUimlRUuohodXfiCB6P7iIlDv32t4OPowr1YoNH5ekuW5QLLKK6wuld0BMzte/nkCwJMAbkxjviL1HD3a3OMVe/cC/f0rH+vvjx7PYnm+azkESK4nubHyO4DPAXi11fmKNDI01NzjFWNjwNRU9M1PRj+npqLHs1ie79LYEvgkgBdIvgLgRQAHzOwnKcxXpK6k3+hA9IGfnY025WdnGwdAq8vzWcshYGZvmtn15dsfmVmbrxJxpViM9ru7uqKfxWL96ZN+oyeVdHnNvq7cmVnut5GRERNZbnrarL/fDFi69fdHj7czV68LQMlifh5TOUTYLB0dkNXqtdjPzuZdTXpcvS6SMxaz5646C4kXOrXlvR1el0JAvNCpLe/t8LoUAuKFTm15b4fXpRAQL+Td0p+XdnhdahgU6UBqGBRveH+M3DMu1lea5xMQWaFYBMbHgTNnovtzc9F9wK/NYV+4Wl/aHZDMdOqx/6ykub60OyBeaIdj5D5xtb4UApKZdjhG7hNX60shIJlph2PkPnG1vhQCkpl2OEbuE1frSw2DIh1IDYMiEptCQCRwCgGRwCkERAKnEJDENC7AjbTXu8YOSCIaF+BGFutdhwglEY0LcCPuetchQsmcxgW4kcV6VwhIIhoX4EYW610hIIloXIAbWax3hYAkonEBbmSx3tUwKNKBnDQMkuwm+WuS+9Oap4hkL83dgd0ADqc4P5F4ikVgy5Zo+5iMfs+y51Ley8tYKiFAciuAzwN4JI35icRWLAL33AOcPLn02MmTwFe+Uv+Dubrb3cTE0v0tW4CNG5c+5Bs2RI9V/nbXXWuXd889bRsEqbQJkPwRgH8CsBHA35jZ7fWmV5uApKJYBO6+G1hYqD1Nd3f090IB2LkTeOqpqLcNGV0kOE0e9ZTKtU2A5O0ATpjZTIPpxkmWSJbm5+dbXayErtJ/tl4AAEt/n5sD9u1b6m6XRYN4m/aUSmN3YBuAL5KcBfAYgFtITq+eyMymzGzUzEYHBwdTWKzkybvBQrt3L3Wg94XjnlLL3yPgj6+L+7yWQ8DMvmlmW81sGMCdAH5mZrtana/4o/KlOzcXfYFWBq04CYJKo9zyfXJfnD7tLB1Xv0dA77q4z1VnIWlocnLtl+6ZM9Hjuar8p/sYAEBUl6N0rPYexaXOQtJQV1f1XWgSWFzMsZBaQ+h846CBcO17NAqzEmM9N5uSpJN4M1ioXRreHNSpAUSSKW8GC23enPMCE3LQQFjtPYpLISANabBQkxwMpVz9HgHnz8V9rtoEpH3UapzwjQc16sxC0pl0xpJMKASkfbTDGUu6u11X0DSFgLSPdmiEuPlm1xU0TSEg7WVgwHUF9R054rqCpikEpL08+GClc7yf2qUvwzIer02RKsbGgI9/3HUVtbVh46VCQOrybvQgAJw65bqC2hw3Xk5MAD09ADAyEvc5CgGpyavRg8u14bdtHiYmolMmNDrFwmrqLCQ1eXupsdUX5POJw5XT07M8ADSASFLg7aXGlveR9Y3DldPsFkCFQkBq8mb0YDVjY9E37vr1ritZyeHKSdpPSSEgNXkzerCWYhH44IPmnzcwAKyLfeKd5uzcmc18Y6hcorxZCgGpyevRg5V2gSROngTOxR5k15ynnspmvjF85zvAffc1v0WghkFpT76eZSj30y3VKkOjCKXT+RgAgCcNJs1RCIikxasGk/gUAiJp8KrBpDk9rgsQSaRyeTEf3Hdf1CrXprQlIO0p6ZGBLDz8sAd9qZNTCEh78umb1yy6MGqbBoFCQGLxbjSh8wJWWVhwProq6bUIYWa530ZGRkzax/S0WX+/WfSVF936+6PHnSkUVhbky61QcLI61r5HI2YxP4/qLCQNeTma0NfTjzvqLLT2PdIoQkmRl6MJfe2U46iuVt6LlkOA5CUkXyT5CsnXSO5pdZ7iFy9HE7Zy3a2sOOws5PpahGcB3GJm1wO4AcAOkjelMF/xhJejCSujm3w5+3B3t9POQk6vRVhulzhdvttbvnm4syZJeTuacGwMePttYHq6+TAYGEj3XASLi05XiPNrEZLsBjAD4AoAD5nZ31WZZhzAOAAMDQ2NzPk6AETaV9zGwuUtmmmdqsz5OddWyn0UoZktmNkNALYCuJHktVWmmTKzUTMbHRwcTGOxIivF2TFevR9Tb7eitzfeyUec7xu1JtWjA2b2DoCfA9iR5nxFYqm1Y1y5WEmt/ZjluxXL93m+9z3g0UeXHhsYWAqLypk7vNk3Sq7l3QGSgwDOm9k7JPsA/BTAP5vZ/lrPUT8ByUyxCExORsfMhoaiYGjjD2hSzewOpDGK8FMA/rXcLtAF4PF6ASCSqbGxID/0rWg5BMzsNwA+k0ItIuKAegxKYt4NKgpE2utdJxWRRFYfWatcogzQ1niWsljvGkAkiXg5qCgAcde7zjYsmfNyUFEAsljvCgFJxMtBRQHIYr0rBCQRLwcVBSCL9a4QkES8HVTU4bJY72oYFOlAahgUkdgUAiKBUwiIBE4hIJlS1+LmuFhf6jYsmVHX4ua4Wl86OiCZUdfi5qS5vnR0QLygrsXNcbW+FAKSGXUtbo6r9aUQkMyoa3FzXK0vhYBkRl2Lm+NqfalhUKQDqWFQRGJTCIg3OrVjke+vS52FxAud2rGoHV6X2gTEC53ascjV61KbgLSdTu1Y1A6vSyEgXujUjkXt8LoUAuKFpB1l8m50a3Z5bdFhysxyv42MjJjIatPTZoWCGRn9nJ5uPH1/vxmwdOvvb/y8VupLsrxmX1caAJQs5uex5Q80gMsAPAfgMIDXAOxu9ByFgKShUFj5gazcCoXGz03ywWxleXlrJgTSOER4AcDXzewlkhsBzJB8xsz+J4V5i9SUtNEt6WG7dmjkS6LlNgEz+52ZvVT+/X1EWwSXtjpfkUaSNrpNTi4FQMWZM9HjWSzPd6k2DJIcRnSZ8l9W+ds4yRLJ0vz8fJqLlUAlbXRL+o3eFo18ScTdb2h0A7ABwAyAOxpNqzYBSUve+/YuGvmSQBNtAqn0GCTZC2A/gKfN7NuNplePQXFpdZsAEH2jd9Iw51x7DJIkgO8COBwnAERc03kOVmp5S4DknwL4bwCHACyWH/57M3uq1nO0JSCSrWa2BFo+RGhmLwBgq/MRETfUbVgkcAoBkcApBEQCpxAQCZxCQCRwCgGRwCkERAKnEBAJnEJAJHAKAZHAKQREAqcQEAmcQkAkcAoBkcApBEQCpxAQCZxCQCRwCgGRwCkERAKnEBAJnEJAJHAKAZHAKQREAqcQEAmcQkAkcAoBkcApBEQCl0oIkHyU5AmSr6YxP1cmDkyg5/4ecA/Rc38PJg5MuC5JJHNpbQl8H8COlOblxK3/div2lfZhwRYAAAu2gH2lfVWDoHioiOEHhtG1pwtbvrUFW761BV17ujD8wDCKh4o1w6TW48vnV5mHSF5avjT5xRmRwwD2m9m1jab16dLkEwcm8HDpYRiqr4duduPCP1y4eL94qIjxH4/jzPkzTS3n0xs+jeOnj695/JLuS7BgCzi/eP7iY/29/Zj6whTGrhtrahkiFc1cmjy3NgGS4yRLJEvz8/N5LbauiQMT2FfaVzMAgGiLoPJNzT3Erid2NR0AAKoGAAB8tPDRigAAgDPnz2Dy4GTTyxBJIugtgZ77ey5u/tdCEH29fYk++K0giMV/XMx1mdI5vNwS8FGjAACAj3V/LPcAAIChTUO5L1PCFHQIdLO74TQfLXyUQyVrnT53Wg2Ekou0DhH+EMAvAFxF8hjJr6Yx36yNj4y7LqGmkx+exPiPxxUEkrnU2gSa4UubAABwD12XUFdhUwGzX5t1XYa0GbUJdJCj7x51XYJ0uOBDYKBvwHUJdamBULIWdAgUDxVx+txp12XURBB7t+91XYZ0uGBDoNLz7+zCWdel1HTv6L3qNSiZCzYEJg9OOjn+H1cXu7BtaJvrMiQAwYaA7w1ui7aoQ4SSi2BDoB0a3DSGQPIQbAjs3b4X/b39rstoyPctFml/wYbA2HVjmPrCFAqbCgDidSFuhEi/41E7bLFIe+txXYBLlZb3JOcHWI1g3SHJSfT39usQoWQu2C2BirSOEhgsla2JioG+AZ1YRHIR9JYAkO4+d5yhyXFsv3w7nr3r2VTmJdJI8FsCPu5zH3zroM41KLkJPgSu2HyF6xKqmnt3Tv0EJBdBh8DEgQkcfOug6zJqUj8ByUPQITA1M5Xr8pIcQlQ/Acla0CGQVkNeXD+44wcX+yXE5WObhXSW4I8O5KlyuG/XE7tiTa9+ApKHoLcE8tTsFoD6CUhegt4SGOgbwMkPT2a+nN6uXuy8cieGHxjG3LtzsZ6zYd0GBYDkIugQyIvB8MhLj6y50lA9ahCUvAQdAqc+PJXLci4sXmg80SpqEJS8BN0m4OsHTQ2CkqegQ8DXcwr09fS5LkECEnQIjF03hruvv9t1GWvo6kOSp6BDAAAef+1x1yVUpS7DkpfgQyCPQ4RJ6QiB5CGtC5LuIPkGySMkv5HGPAXY3LfZdQkSgJZDgGQ3gIcA3AbgGgBfJnlNq/PNy/re9a5LqOm9s++pXUAyl8aWwI0AjpjZm2Z2DsBjAL6UwnwzVzxUxLmFc67LqOn84nm1C0jm0giBSwH8dtn9Y+XHvDd5cLKpXnwuqF1AspZGCFQbJL/mtLskx0mWSJbm5+dTWGzrWv2ApXWK8cKmQs0BRr52aJLOkUYIHANw2bL7WwEcXz2RmU2Z2aiZjQ4ODqaw2NYl+YBVzihc2FTAvaP3ttzZaF33OuzdvrdqxyX1HJQ8pDF24FcAriR5OYD/BXAngL9IYb6Z27t975prDvR29WJd9zp8cP6Di48N9A3gwdserDqqb9vQNkwenMTRd49iaNPQxQ9ttcd2/9fuFYckq8139fM0klCyRrPWL5hBcieABwB0A3jUzOp+fY2OjlqpVGp5uWkoHirqgycdh+SMmY3GmjaNEGiWTyEg0omaCYHgewyKhE4hIBI4hYBI4BQCIoFTCIgETiEgEjiFgEjgFAIigVMIiAROISASOIWASOAUAiKBUwiIBE4hIBI4hYBI4BQCIoFTCIgETiEgEjiFgEjgFAIigVMIiAROISASOIWASOAUAiKBUwiIBE4hIBI4hYBI4BQCIoFrKQRI/jnJ10gukox18UMR8UurWwKvArgDwPMp1CIiDvS08mQzOwwAJNOpRkRy11IINIPkOIDx8t2zJF/Na9kxbAHwtusilvGtHsC/mlRPfVfFnbBhCJB8FsAfVPnTpJn9R9wFmdkUgKnyPEtm5k0bguppzLeaVE99JEtxp20YAmZ2a2vliIjPdIhQJHCtHiL8M5LHAPwJgAMkn4751KlWlpsB1dOYbzWpnvpi10Mzy7IQEfGcdgdEAqcQEAmcsxDwpcsxyR0k3yB5hOQ3XNVRruVRkid86UNB8jKSz5E8XH6vdntQ0yUkXyT5SrmmPR7U1E3y1yT3u64FAEjOkjxE8uU4hwpdbgk473JMshvAQwBuA3ANgC+TvMZVPQC+D2CHw+WvdgHA183sDwHcBOCvHK8fADgL4BYzux7ADQB2kLzJcU27ARx2XMNqnzWzG+L0XXAWAmZ22MzecLX8shsBHDGzN83sHIDHAHzJVTFm9jyAU66Wv5qZ/c7MXir//j6if/RLHddkZna6fLe3fHPWuk1yK4DPA3jEVQ2tCr1N4FIAv112/xgc/5P7iuQwgM8A+KXbSi5ufr8M4ASAZ8zMZU0PAPhbAIsOa1jNAPyU5Ey5u35dmY4dSKvLcYaqjXzSMdNVSG4A8O8AvmZm77mux8wWANxA8vcBPEnyWjPLvR2F5O0ATpjZDMmb815+HdvM7DjJTwB4huTr5a3MqjINgTbocnwMwGXL7m8FcNxRLV4i2YsoAIpm9oTrepYzs3dI/hxRO4qLxtRtAL5IcieASwD8HslpM9vloJaLzOx4+ecJkk8i2u2tGQKh7w78CsCVJC8nuQ7AnQD+03FN3mA0Rvy7AA6b2bdd1wMAJAfLWwAg2QfgVgCvu6jFzL5pZlvNbBjR/87PXAcAyfUkN1Z+B/A5NAhIl4cIk3Y5To2ZXQDw1wCeRtTo9biZvZZ3HRUkfwjgFwCuInmM5Fdd1VK2DcBfArilfLjp5fK3nkufAvAcyd8gCvFnzMyLQ3Oe+CSAF0i+AuBFAAfM7Cf1nqBuwyKBC313QCR4CgGRwCkERAKnEBAJnEJAJHAKAZHAKQREAvf/u5bR9+ycJbkAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "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\")"
      ]
    },
    {
      "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",
      "execution_count": 32,
      "metadata": {},
      "outputs": [],
      "source": [
        "target_proxy_mat = interaction_mat(targets, target_proxies)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check its numerical rank and shape:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "24"
            ]
          },
          "execution_count": 33,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "numerical_rank(target_proxy_mat, eps)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(200, 25)"
            ]
          },
          "execution_count": 34,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "target_proxy_mat.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now compute an ID (row or column?):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "metadata": {},
      "outputs": [],
      "source": [
        "idx, proj = sli.interp_decomp(target_proxy_mat.T, nproxies)"
      ]
    },
    {
      "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",
      "execution_count": 36,
      "metadata": {},
      "outputs": [],
      "source": [
        "target_skeleton = idx[:nproxies]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check that the ID does what is promises:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "6.738945834835623e-15"
            ]
          },
          "execution_count": 37,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "P = np.hstack([np.eye(nproxies), proj])[:,np.argsort(idx)]\n",
        "tpm_approx = P.T @ target_proxy_mat[target_skeleton]\n",
        "\n",
        "la.norm(tpm_approx - target_proxy_mat, 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the chosen \"skeleton\" and the proxies:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHJVJREFUeJzt3XtwHFeVBvDvzEPvly3JDrZeKZIimwVskICAFwriACZAwlJFVYKzBlcKlTFbFSDAwirsYhIXbJaQQGEcVCFgx4IQdqF2iZeEYDCBKh6RQoIJTqiQWImsgGQ5lh96jDRz9o+j9tjSjDSPnuke9fermpKnp6f7quV7+va593aLqoKIgivkdQGIyFsMAkQBxyBAFHAMAkQBxyBAFHAMAkQBF3FjIyJyBMApAHEAs6ra5cZ2iajwXAkCc96iqsdc3B4RFQEvB4gCTtwYMSgizwJ4EYAC+Iaq9qZYpxtANwBUV1d3XnLJJXnvl4Lt+HFgcBBIJJLLQiGgvR1YuTL99wYG0n/W2en+/rwwMDBwTFWbM1nXrcuBDao6LCKrADwkIk+q6sPnrjAXGHoBoKurS/v7+13aNQVVR8f5FRKw94kEsNh/r44Oq8zztbcv/b1c9ucFEUnxG6bmyuWAqg7P/RwB8EMAr3Vju0SLee657JY7du4EqqrOX1ZVZcsLsT+/yzsIiEi1iNQ6/wbwNgB/zHe7REtpa8tuuWPzZqC31878Ivazt9eWF2J/fudGS2A1gF+JyOMAfgdgv6o+4MJ2iRaV6xkdsAp/5Ig15Y8cWToA5Ls/P8s7CKjqM6q6bu7196pa4oeEvNLXZ9fdoZD97OtbfP1cz+i5ynV/2f5eRaeqRX91dnYq0bn27VOtqlIFkq+qKlteyrz6vQD0a4b10ZUuwmyxd4DmWyxjf+RIsUvjHq9+LxEZ0AxH7nKwEPnCcs28l8LvxSBAvrBcM++l8HsxCJAvLNfMeyn8XgwC5AvFzvQXSyn8XkwMEi1DTAySb/i+j9xnvDhebt5PgOg8fX1AdzcwMWHvBwftPeCv5rBfeHW8eDlABbNc+/4Lxc3jxcsB8oVS6CP3E6+OF4MAFUwp9JH7iVfHi0GACqYU+sj9xKvjxSBABVMKfeR+4tXxYmKQaBliYpCIMsYgQBRwDAJEAccgQBRwDAKUM84L8Ibbx51zBygnnBfgjUIcd3YRUk44L8AbmR53dhFSwXFegDcKcdwZBCgnnBfgjUIcdwYBygnnBXijEMedQYBywnkB3ijEcWdikGgZ8iQxKCJhEfm9iNzv1jaJqPDcvBy4AcBhF7dHlJm9e4HWVhs9s3YtcOedwMxM4fY3M2P7WLvW9tnaamUoUa4EARFpAfBOAHe5sT2ijO3dC2zbBgwN2fM+h4eBj33MLpQXCwR791pKPRSyn9/6FjA+DoyNAceOAbt2JSv5ihVAY2Oywl9/ve1jeNj2OTRkZSjVQJDpk0sXewH4LwCdAN4M4P6l1udTickVe/aohsPnP/LXeYVCqiKqa9eq7tqleuKE6pkz9nPXLtXKyvPXr6hQ3b1bdWxM9dZb7X2q7S72am31+oichSyeSpx3S0BE3gVgRFUHllivW0T6RaR/dHQ0391S0DktgHg89eeJhFXNo0eBG28Evv1t4PnngVjM+tMmJ89ff2oKuOUWG497++32PltDQ9l/xwfcuBzYAOAqETkC4F4Al4vIvvkrqWqvqnapaldzc7MLu6Vi8t1koZ6ehRU5nakp4AtfAMJh4PRp4IUXUq83PAz87W/AX/+aW5laWnL7nkvO/RsBr3xFpt/LOwio6mdUtUVVOwBcA+Bnqnpdvtsl/3AmrQwO2snVmbTiSSBwkoDZnnVHRoA//AH47W+BpqbU67zkJZZHyPUkdeaMZ3mB+X8jIFqW6Xc5WIiW1NOTnLXmmJiw5UV1bhIwW3V1wMmTdplw1VVA2bw6Ul5uyb6ZGRt5M//zc5WVAW9/O1Bbe/7y48c9SxCm+htlytUgoKoHVfVdbm6TvOebyUI33bT4JUBZGbBx48IKHIkAGzZY5n9oCHjd66yir1xpn69aZQHgNa+xy4X3vAf4yEesxSAC1NRYhRexZd3dwKc+BVRWLizD5KSVs8jy+VvwfgK0pLa21NNXiz5ZaLEWQFMTcM01wMUX2+t73wNefNEqb1cXcOGFFhzGx4FnnwU6O4Err7SKfOKEXQKsWmUtgePHgauvtu3F41axq6ttP5OT9pqetkuMbMtZIOn+RplgEKAl7dx5/o0sAI8mCzU0WMWer7oa+NKXLPE3PAy8/vXAG99oldl5VVTYmby21pr+09P2vqzM3s/M2C84Pe1cVNv2Tp2ydeJxW7+yEohG7frf2c58HiQIU/2NMsUgQEtyJqf09Fizs63N/tMVfbKQSOrliYT1AJSXA/X1VnGjUXvf0GCfx2I2CGhqyipyc7NdDkxPJyt7VZW9ysutxRCLWeUvL7faNTtr68fjwNe+ljoAiFhXY5HN/xupzsQy/S4nEFHpCIWSZ+n5du+27L5zJj950kb5nThhzfaJCTvbT0/bCMAVKyxgNDRYSyISsUoeiVjQGB21Zn9DA7B6tX02Pm7bq6gA3vc+W28+kdTLiyybCURsCVDpaGmxAT/z1ddbxW9stDP5iRN2dq+rs4qcSFjOoLo6GURE7H1TkwWMhgar+KOjlhxsbLR1y8vtc+c7TqBJV9E9OKnmi12EVDpuucUq5bmiUeuuU7V/h8PJvv5o1FoPF1xgTX/n+r+uzs7mTsZf1QYIxWK2vKnJtuGsn0gkg084bEEhlKbqhMOFPQYFwCBApWPLFmDHDsviO911W7cCb3mLnbnDYausiYSdsWtr7a4bzc1WaUMhWxaJWIshErFWRFWVLVe1y4aRERs5GItZS0LEfo6MWOJxYsKSj6l88INFPSRuYBCg0nL99cAvfmGtAhGbLfjZzwKPPGIZ+9lZO9NHIpYEdCr07KwlBMfHbbmILYvHrcVQVWXN/BMnksHh9Gl7nTljYwXicUsuvviiDTjq7EwmK0MhC1K33ebt8ckBcwJUWurrge98B7j55uQkH2fqb0uLtQqmpqziqtp1/+SkVebJSTuLr1mTPPNPT1uTPxazqcP19cmkYDxu/y4vt8/Gxy2gOHMQtm4FbrjBLiGqqqx1MDvr7fHJAYMAlZZoFPjP/1w4y29qyvotGxrsTP7Sl1pFPXnSxglMTtqAoQsusLP9sWNWuVetsuz/0aM2otAZGDQ1ZfuamrIAUFdnAaSx0fIF5eXJewxMTNglx6lTtv8Sw8sBWpTvZg8CVmFTOXbMzu6jo8Bjj9nPsjKroM6AoETCzuSrVlmycHbWKn40aus7A4rq6215KJT8mUjYWT8UsnzEwYPAhz9seYBrrgHuvXfhrYCLbPt2i4F2rZIZBgFKy1ezB8+VbkReU5N95gwFnpmxSj41ZZ/FYtacr621IFBVZZV6etqu8+vrk70LThJxxQr73Mk1TE3ZZcHBg8Ddd1srA7Dtfu5zwHe/W6yjsMD27TZcIt0tFtJhEKC0fDN7cL5bblk4eaesDHj/+62CiljT/sQJq8TV1VaJnUrd3GxRrazMmvTl5VbpIxHL/j/zjP2itbW2zurVyfkCdXUWaL7//YW3L5ua8mTykKO3N7fvMQhQWr6ZPTjfli12o8/WVqvwK1daE+UNb0gOEW5qskARClnFdSp0ebnlCuJxq9wzM1bBx8eTXYc1NVbhnVGE9fW234YG63p0Rg+m4uHdhbJtATgYBCgtXz9qbMsWi0ZHjwL79gHr1iXv9jc7mxwN6Jy9Gxrsur+mxj5bsyaZB/jLX+zfTm4gkbBA8be/WaWPxSxYVFdb0HEGFKXi4d2Fch2nxCBAafn+UWPOWXztWjvLO7MAV6+2Zv/MjE0bHhy0wLB6ta3b2Gjfe/xxCyIvvpicF+CMNVi9OjnUuL7egsHERHKcwbp1qcu0aVNxj8E5nEeUZ4tdhJSWb2YPpjIzYxUyHLYK/8pX2tnbuWfAmTP2WWWlLQuF7N+nT1tr4fhxu4xYudLO+M778vLkZYFzc5JIxL577FiyZfH736cu1wMPFO8YzPP1r9vP3t7sLg3YEqBFbd5sz71PJOynLwIAYGfle+8FLrvMhga/9a2WrBsZsSb92rXJ0YDRqJ3dx8asL//o0WRQmJy0noKpKUsIOl2E4bD1idbU2HorViR7E5wbj6Ti8R2Hv/51Z7zSwKJ3/z4XgwCVpj17gE9+0ib+qNpZevduyw88/bQFg3A4eTZ3piGrJm8UMj1trYa6OmtJOEnDSCQ54SgcTk4aikRsnZoau6RIxeM7DueCQYBK0xe/uHDU4MwM8OMfW7MlGj2/Ce8sq621Sjw0BDz6qAWMo0ft9HnppcBFF1nld3oTJiftWujJJ4GnnrLblZ86ZYOD5s9orKz05IYi+WIQoNKU7tkAJ09aNn942ILCmTNWaZub7ewdj1ur4OmnLThUVNi2BgftLL5mjbUWZmasNVFebpcR4+O2bGzMrosuuwz40IcsgShi3ZV33mm9FiWGiUEqTWvWpB4+XF9v+YJVq+zMX1dnQSAWs16AkRE7qycS9qqstODgrOdU6PJyax0cPWrdi4ODydmHdXUWDNavB975TuBNb/J8uHA+2BKg0nTzzQtvLR6N2izCUMiShW1tNmnIucHI7Kw17c+cSSb/Jift+r+21oJCXZ0lAWdnk5cD8bjlChIJCw4VFfZZS4ut99e/FvYpyAXGIEClaetWu2W4c4cfEWuir1tnZ3fV5E1AVq60inr0qFV850aiIlZ5//xnO9MfO2YtBmestJNLmJy0PMLx45YTmJ62FseBA3Z/g4susmBTok8lZhCgjPhuNuHevdYn79zrTxX45S9tmvHAQHKcQHm5fVZba62BWMya9Q8/bAMgPvEJu1vRwYN2iTE9bUFgbMy6DBMJCxjOKMHKSss7PPII8NWvWqBxHnzq8ePJPXsWIS1/vpxNeNNNqZ8cPDFhBdu/30YAJhLJKcQzM9a0/9OfgAcftM8BywXcd5+d2Z27EZWXWxCZnrbA4TzCrKbG7knQ17fwluMePX0I4LMIqcB8OZtwsUE5sRhwzz2WxX/+eTtLnzyZvCvxAw8svIaPxezuRCdPJu8ONDtrr9ZW4JJLgFe/2vIFk5M25yDbchVQPs8iZO8ALcmXswnT3X7cceyYnZWPH7eK+/a3Wx5h7dpkC2C+kREbEjk2Zmf7rVutZXDPPba9pibrDYhEkjcZSVUuD+T1t1DVvF4AKgD8DsDjAJ4AsGOp73R2diqVjvZ2Z6jd+a/2dg8LtWePamVl6oKlepWVqX7846q/+Y3q6tWZfSccVo1EFi4Lh1OvX1lp5fLAwr9Rp2qGddiNy4FpAJer6joA6wFsEpHLXNgu+YQvZxM69xRYsSKz9WMxS9o1Ntpw44qKpb8Tjy+8cWg8nnp2Tjjs6WChVH+jTOUdBOYC0em5t9G5V+k9hoXS2rzZZqa1t1uvWnu7vfd8MtGWLdbc37Mns2b42JhV/quuAm6/3S4NnKcK5SuR8HS04Py/EVDkZxGKSBjAAICLAOxS1X9JsU43gG4AaGtr6xzM9TnKROks9qxCwLoAn3jCaolztyBnSvL69elvYJqJ1lYf3HIpKZtnEbrSO6CqcVVdD6AFwGtF5OUp1ulV1S5V7WpubnZjt0TnW6w1UFFhDymZnT2/3RyNWkC46aaFlwjONOT5y+bfwqdEJw45XO0iVNUTAA4C8O72KhRcqW5ACtjY/9tuA6691ip8qoq9bRvwjW8k71vY2grcdZe9nGUtLTZA6KtfTV5KlPDEIUfelwMi0gxgRlVPiEglgJ8A+A9VvT/dd/hociqYvXvtrD40ZJX2lltKuoLmqtiPJn8JgD1zeYEQgPsWCwBEBbVlSyArfT7yDgKq+gcAr3KhLETkAQ4bppz5blJRQLh93DlsmHLiTFhxxqs7k4oAH4wfWMYKcdxdGSeQLSYGS19Hh/0HnK+93ebtUGFketyLPk6AgseXk4oCoBDHnUGAcuLrR5QtY4U47gwClBNfTioKgEIcdwYByolvJxUtc4U47kwMEi1DTAwSUcYYBIgCjkGAKOAYBKigOLQ4O14cLw4bpoLh0OLseHW82DtABcOhxdlx83ixd4B8gUOLs+PV8WIQoILh0OLseHW8GASoYDi0ODteHS8GASoYDi3OjlfHi4lBomWIiUEiyhiDAPnGch1Y5Pffi4OFyBeW68CiUvi9mBMgX1iuA4u8+r2YE6CSs1wHFpXC78UgQL6wXAcWlcLvxSBAvpDrQJliJ92y3V9JDJhS1aK/Ojs7lWi+fftU29tVReznvn1Lr19VpQokX1VVS38vn/Llsr9sfy83AOjXDOtj3hUaQCuAnwM4DOAJADcs9R0GAXJDe/v5FdJ5tbcv/d1cKmY++yu2bIKAG12EswBuVNVHRaQWwICIPKSqf3Jh20Rp5Zp0y7XbrhSSfLnIOyegqi+o6qNz/z4FaxGszXe7REvJNenW05MMAI6JCVteiP35nauJQRHpgD2m/LcpPusWkX4R6R8dHXVztxRQuSbdcj2jl0SSLxeZXjcs9QJQA2AAwHuXWpc5AXJLsa/tvUjy5QJZ5ARcGTEoIlEA9wN4UFW/vNT6HDFIXpqfEwDsjL6cpjkXdcSgiAiAbwI4nEkAIPIa73NwvrxbAiLyDwB+CeAQgMTc4n9V1f9L9x22BIgKK5uWQN5dhKr6KwCS73aIyBscNkwUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcAwCRAHHIEAUcK4EARG5W0RGROSPbmzPK9v3b0fk8xHIDkHk8xFs37/d6yIRFZxbLYFvA9jk0rY8ccXeK7C7fzfiGgcAxDWO3f27UwaCvkN96LijA6EdITTd2oSmW5sQ2hFCxx0d6DvUlzaYpFt+7vacbRAVS96PJj+7IZEOAPer6suXWtdPjybfvn877uy/E4rUxyEsYcz+2+zZ932H+tD9o25MzExktZ81NWswfHp4wfKKcAXiGsdMYubssqpoFXrf3YvNr9ic1T6IHNk8mrxoOQER6RaRfhHpHx0dLdZuF7V9/3bs7t+dNgAA1iJwztSyQ3DdD67LOgAASBkAAGAqPnVeAACAiZkJ9BzoyXofRLkIdEsg8vnI2eZ/OgJBZbQyp4qfD4Eg8e+Jou6Tlg9ftgT8aKkAAADl4fKiBwAAaKtvK/o+KZgCHQTCEl5ynan4VBFKstDp2GkmCKko3Ooi/C6AXwN4mYgMicj1bmy30Lo7u70uQlpjk2Po/lE3AwEVnGs5gWz4JScAALJDvC7Cotrr23Hko0e8LgaVGOYElpHnxp/zugi0zAU+CDRWNnpdhEUxQUiFFugg0HeoD6djp70uRloCwc6NO70uBi1zgQ0Czsi/6fi010VJa1vXNo4apIILbBDoOdDjSf9/pkISwoa2DV4XgwIgsEHA7wm3hCbYRUhFEdggUAoJN84hoGIIbBDYuXEnqqJVXhdjSX5vsVDpC2wQ2PyKzeh9dy/a69sBZDaEeCkC9wcelUKLhUpbxOsCeMnJvOdyf4D5BLLolORcVEWr2EVIBRfYloDDrV4ChbrSmnA0VjbyxiJUFIFuCQDuXnNnMjU5Exsv3IifbvmpK9siWkrgWwJ+vOY+8OwB3muQiibwQeCilRd5XYSUBscHOU6AiiLQQWD7/u048OwBr4uRFscJUDEEOgj0DvQWdX+5dCFynAAVWqCDgFuJvEzd8957zo5LyJQfcxa0vAS+d6CYnO6+635wXUbrc5wAFUOgWwLFlG0LgOMEqFgC3RJorGzE2ORYwfcTDUVx5cVXouOODgyOD2b0nZqyGgYAKopAB4FiUSjuevSuBU8aWgwTglQsgQ4CxyePF2U/s4nZpVeahwlBKpZA5wT8WtGYEKRiCnQQ8Os9BSojlV4XgQIk0EFg8ys24wPrPuB1MRbg04eomAIdBADgvifu87oIKXHIMBVL4INAMboIc8UeAioGtx5IuklEnhKRp0Xk025sk4CVlSu9LgIFQN5BQETCAHYBeAeASwFcKyKX5rvdYqmOVntdhLROTp9kXoAKzo2WwGsBPK2qz6hqDMC9AK52YbsF13eoD7F4zOtipDWTmGFegArOjSCwFsDz57wfmlvmez0HerIaxecF5gWo0NwIAqkmyS+47a6IdItIv4j0j46OurDb/OVbwdy6xXh7fXvaCUZ+HdBEy4cbQWAIQOs571sADM9fSVV7VbVLVbuam5td2G3+cqlgzh2F2+vbsa1rW96DjcrCZdi5cWfKgUscOUjF4MbcgUcAXCwiFwI4CuAaAO93YbsFt3PjzgXPHIiGoigLl+HMzJmzyxorG/GVd3wl5ay+DW0b0HOgB8+NP4e2+razlTbVsht+fMN5XZKptjv/e5xJSIUmqvk/MENErgRwB4AwgLtVddHTV1dXl/b39+e9Xzf0HepjxaNlR0QGVLUro3XdCALZ8lMQIFqOsgkCgR8xSBR0DAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRwDAJEAZdXEBCR94nIEyKSEJGMHn5IRP6Sb0vgjwDeC+BhF8pCRB6I5PNlVT0MACLiTmmIqOjyCgLZEJFuAN1zb6dF5I/F2ncGmgAc87oQ5/BbeQD/lYnlWdzLMl1xySAgIj8FcEGKj3pU9X8y3ZGq9gLondtmv6r6JofA8izNb2VieRYnIv2ZrrtkEFDVK/IrDhH5GbsIiQIu3y7CfxSRIQCvB7BfRB7M8Ku9+ey3AFiepfmtTCzP4jIuj6hqIQtCRD7HywGigGMQIAo4z4KAX4Yci8gmEXlKRJ4WkU97VY65stwtIiN+GUMhIq0i8nMROTz3t7rBB2WqEJHficjjc2Xa4YMyhUXk9yJyv9dlAQAROSIih0TksUy6Cr1sCXg+5FhEwgB2AXgHgEsBXCsil3pVHgDfBrDJw/3PNwvgRlX9OwCXAfiIx8cHAKYBXK6q6wCsB7BJRC7zuEw3ADjscRnme4uqrs9k7IJnQUBVD6vqU17tf85rATytqs+oagzAvQCu9qowqvowgONe7X8+VX1BVR+d+/cp2H/0tR6XSVX19Nzb6NzLs+y2iLQAeCeAu7wqQ76CnhNYC+D5c94PweP/5H4lIh0AXgXgt96W5Gzz+zEAIwAeUlUvy3QHgE8BSHhYhvkUwE9EZGBuuP6iCjp3wK0hxwWUauYT+0znEZEaAP8N4KOqetLr8qhqHMB6EWkA8EMRebmqFj2PIiLvAjCiqgMi8uZi738RG1R1WERWAXhIRJ6ca2WmVNAgUAJDjocAtJ7zvgXAsEdl8SURicICQJ+q/sDr8pxLVU+IyEFYHsWLZOoGAFeJyJUAKgDUicg+Vb3Og7KcparDcz9HROSHsMvetEEg6JcDjwC4WEQuFJEyANcA+F+Py+QbYnPEvwngsKp+2evyAICINM+1ACAilQCuAPCkF2VR1c+oaouqdsD+7/zM6wAgItUiUuv8G8DbsESA9LKLMNchx65R1VkA/wzgQVjS6z5VfaLY5XCIyHcB/BrAy0RkSESu96osczYA+CcAl891Nz02d9bz0ksA/FxE/gAL4g+pqi+65nxiNYBficjjAH4HYL+qPrDYFzhsmCjggn45QBR4DAJEAccgQBRwDAJEAccgQBRwDAJEAccgQBRw/w8c/w9Q3nUpegAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "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\")"
      ]
    },
    {
      "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",
      "execution_count": 16,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "1.8776453951593488e-09"
            ]
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "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)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Source Skeletonization"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {},
      "outputs": [],
      "source": [
        "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)])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAE9pJREFUeJzt3X2MHOV9B/Dv9/bOwWs7bnx3SQPm9pCCoAgE0a0Qlf8hGEXGIYkaKRLpOqAQaQXXVo6Uqkp0UitHOlXljwj+iE1OhKTNbYOiBtSCaRLjkKJKUcgegRjqIFlwd3EdyW/l5Tjwy92vf8yu78X7MrM7u8/MPt+PtDrveG6eZ8fe784+83tmaGYQEX/1ue6AiLilEBDxnEJAxHMKARHPKQREPKcQEPFcfxwbITkL4F0ASwAumlk+ju2KSOfFEgIVnzKz0zFuT0S6QF8HRDzHOCoGSb4J4P8AGIDvmtlUjXWKAIoAsGnTprHrr7++7XZFpLaZmZnTZjYcZt24QuBKMztB8qMADgH4GzN7od76+XzeyuVy2+2KSG0kZ8KOzcXydcDMTlR+ngTwFIBb49iuiHRe2yFAchPJLdU/A/g0gFfb3a6IdEccZwc+BuApktXt/auZ/TSG7YpIF7QdAmb2BoCbY+iLiDigU4QinlMIiHhOISDiOYWAiOcUAiKeUwiIeE4hIOI5hYCI5xQCIp5TCIh4TiEg4jmFgIjnFAIinlMIiHhOISDiOYWAiOcUAiKeUwiIeE4hIOI5hYCI5xQCIp5TCIh4TiEg4jmFgIjnFAIinlMIiHhOISDiudhCgGSG5G9JPhPXNkWk8+I8EtgL4GiM2xMJp1QChoYAMngMDQXLeqW9DoslBEhuB/AZAI/FsT2R0Eol4P77gTNnVpadOQN85SuN35ilEjA6CvT1BT/Hx1eeDw0BW7asvMk3bw6WVf/u3nsvb+/++1MbBDSz9jdC/huAfwSwBcDfmtndjdbP5/NWLpfbblc8VyoB990HLC3VXyeTCf4+lwN27waefRaYmwve3DH8318jlwNmZ+PdZotIzphZPsy6bR8JkLwbwEkzm2myXpFkmWT51KlT7TYrviuVgGKxcQAAK38/NwccOBD8BOIPAACYn49/m10Qx9eBHQA+R3IWwBMA7iA5vX4lM5sys7yZ5YeHh2NoVry2dy+wuOi6F2uNjLjuQUvaDgEz+6aZbTezUQD3APiFme1pu2citVQH5VZ/J0+KhYVUjgv0u+6ASGjVrwBJOwKoOnMm6B8AFApu+xJBLAODUWlgUFoyOrrynT7JEjBA2NWBQZGuScvAW1r6WaEQkPTYts11D8JJ2QChQkAkbpOTrnsQiUJA0uPsWdc9CCdFg4KAQkDSJGWH2WmhEJD0SMNhdibjugeRKQQkPdJwmH377a57EJlCQNJlcNB1Dxo7dsx1DyJTCEi6PPJIMKU3qVJWIwAoBCRtCgXgIx9x3Yv6Ujh4qRCQ9EnyqcI0DF6uoxCQ9Enhp22SKQQkfSYngWzWdS9qm5hw3YPIFAKSPoUCMDUVzNZLGg0MinRJoRBM1920yXVP1krhVxWFgKRXqQS891703xscBDZsiL8/QHAx05TRlYUknapXGWpFJy9N9uyzndt2h+hIQNJpYiKZlxnTmIBIlyT1MmMaExDxWDarYiERb+VywWnLNMx0XEcDg5JO1duLJcGDDwL797vuRct0JCDp1OqZgU549NFU3nSkSiEg6ZSkT16z4MaoKQ0ChYCkU9LecEtLwdFJ0voVgkJA0imJE3UWF5PZryYUApJOSS3KSWq/GlAISDoltSgnqf1qoO0QIHkFyRdJvkLyNZL74uiYSENJvKaAx8VC5wDcYWY3A7gFwC6St8WwXZH6qtcUSMrVhzOZ1BYLtR0CFlioPB2oPLp/v3PxT6EAnD4NTE9HD4PBwXivRbC8nMoAAGIaEyCZIfkygJMADpnZr2usUyRZJlk+depUHM2KBKphQIZbP5cL1l9YCAIkjq8VKRwLqIolBMxsycxuAbAdwK0kb6yxzpSZ5c0sPzw8HEezImuFeSOu/97e6GvFwEC4i4+kdCygKtazA2b2FoBfAtgV53ZFQqk3WFi9WUm9ST6rv1bkcsERRS4HfP/7wOOPrywbHFwJi+o9B1M8caiKZu19fSc5DOCCmb1FciOAnwP4JzN7pt7v5PN5K5fLbbUrUlOpFBTszM8HRwaTk6l+g7aK5IyZ5cOsG8cswo8D+GeSGQRHFj9uFAAiHVUoePmmb0fbIWBmvwPwyRj6IiIOqGJQxHMKARHPKQREPKcQEPGcQkDEcwoBEc8pBEQ8pxAQ8ZxCQMRzCgERzykERDynEBDxnEJAxHMKARHPKQREPKcQEPGcQkDEcwoBEc8pBEQ8pxAQ8ZxCQMRzCgERzykERDynEBDxnEJAxHMKARHPKQREPNd2CJC8muTzJI+SfI3k3jg6lialEjA6GtwBe3Q0eC6SFnHclfgigK+b2UsktwCYIXnIzP4nhm0nXqkEFIvA4mLwfG4ueA7o5riSDm0fCZjZH83spcqf3wVwFMBV7W7XhVY+0ScmVgKganExWN6J9kRiZ2axPQCMApgH8OEaf1cEUAZQHhkZsaSZnjbLZs2AlUc2GyxvhFz7O9UH2Zn2RMIAULaQ71sG67eP5GYA/wVg0syebLRuPp+3crkcS7txGR0NDuXXy+WA2dnk/J5IGCRnzCwfZt1Yzg6QHADwEwClZgGQVPPz0ZZXTU4C2ezaZdlssLwT7YnELY6zAwTwPQBHzezb7XfJjZGRaMurCgVgair4BCeDn1NTzQcFW21PJG5xHAnsAPBlAHeQfLny2B3DdtsSddCt1U90IHjDz84Cy8vBzzBnBVptT4OJEruwgwdxPsbGxjozGlLR6qDb9LRZLhcM6uVynR+ki9qeBhMlLLgYGIyi0wODvTro1quvS+LX9YHBpOnVQbdefV3iVk+GQK8OuvXq6xK3ejIE2hnkS7JefV3iVk+GQKun7ZKuV1+XuNWTA4MivvN+YFBEwktVCKhQJjrtM2kmjusJdIXm7UenfSZhpGZMQIUy0Wmf+asnxwRUKBOd9pmEkZoQUKFMdNpnEkZqQkCFMtFpn0kYqQkBFcpEp30mYaRmYFBEwuvJgUER6QyFgIjnFAIinktMCKi81Q3td0lE2bDKW93QfhcgIWcHVN7qhvZ770rd2QGVt7qh/S5AQkJA5a1uaL8LkJAQUHmrG9rvAiQkBFTe6ob2uwAJGRgUkXilbmBQRNyJ69bkj5M8SfLVOLbnyvjBcfR/qx/cR/R/qx/jB8ddd0mk4+I6EvgBgF0xbcuJO//lThwoH8CSLQEAlmwJB8oHagZB6UgJow+Pom9fH4YeGsLQQ0Po29eH0YdHUTpSqhsm9Zav3l51GyLdEtuYAMlRAM+Y2Y3N1k3SmMD4wXE8Wn4Uhtr7IcMMLv79xUvPS0dKKD5dxOKFxUjtXLn5SpxYOHHZ8isyV2DJlnBh+cKlZdmBLKY+O4XCTRqhk9YkckyAZJFkmWT51KlT3Wq2ofGD4zhQPlA3AIDgiKD6Sc19xJ4n90QOAAA1AwAAPlj6YE0AAMDihUVMHJ6I3IZIK7oWAmY2ZWZ5M8tnMsOJmLQyNTPVdB2CKD5dxNzbNeprO2j+bbdle5pY5A8nE4jm5oDl5ZU/u5q0Uv3+38iHMh9q6ZO/XSNb3ZXtaWKRX5ycIqwGQNXiIjDh4Og3w0zTdT5Y+qALPbncwvkFZwOEExMrAVDl6t9IOi+uU4Q/AvArANeRPE7yq1G34WLSSnGs2P1GQzrz/hkUny46CQJNLPKLk4pBMm/A2rMDrqavch+732gEua05zH5ttqttaopx+iXy7MCaRte1qkkr9bkYINTEIr84CYFcLjmTVgY3DrppOCQXA4SaWOQXJ2cHtm0DklArVDpSwsL5BdfdqIsgJne6+fgtFPSm94W3E4iqlX/nls657kpdD+QfUNWgdJy3ITBxeMLJ+f+w+tiHHSM7XHdDPOBtCLiuyGtm2ZadnSIUv3gbAi4r8sLSHALpBm9DYHLnJLID2eYrOpb0IxZJP29DoHBTAVOfnUJuaw5AuBLiZoj4C4/ScMQi6eYkBGZmkjEzrXBT4dIRQZjJRI0QbDgluRXZgayzU4RV4+NAf39QL9DfHzyX3uLsSKA6M811EMR1lsBgsRxNVA1uHHR+YZHxceDAAWCpko9LS8FzBUFvcT53wHU9et++vtg/wdu185qdeO7e51x3A/39KwGwWiYDXLx4+XJJjsTPHVjN9cy0JH7nPvzm4URca7BWADRaLunkPARc3/LqE9s+4bYDdcy9Pee8TiBT59tNveWSTk5DwPXMtPGD4zj85mF3HWjCdZ1Asc7lFuotl3RyFgJJmJkW5hqDcWrlFKLLOoH9+4EHH1z55M9kguf79zvrknSAkxAYGwsGA13PUmv3tGBUP/zCDy/VJYTlesxi//5gENAs+KkA6D1OphL7qnq6b8+Te0Ktn4Q6Ael9zgcGfRH1CCAJdQLiB6+PBAY3DuLM+2c63s5A3wB2X7sbow+Phr5/weYNmxUA0hVeh0C3GAyPvfTYZXcaakQTh6RbvA6Bs++f7Uo7F5ejl9e5HhAUf3g9JpDUN5oGBKWbvA6BpF5TYGP/RtddEI84CYGzZ5Nxs8vCTQXcd/N9bhpvwOXdh6p0Q1J/OJlFmMnkbXl55Zrj2ay76sGhh4a6coagFS7uPgRcfkNSwO2/kUQXZRah86nEVboN2eUIYvkflpuvGDPdhiz9uj6VmOQukq+TPEbyG61sw/WU4iTatnGbk3Z1Q1K/tB0CJDMAvgPgLgA3APgSyRuibsfVlOJNA5vcNBzCO+fecTIuUO/fwvW0b+mMOI4EbgVwzMzeMLPzAJ4A8PmGjSbkhqSlIyWcXzrf/YZDurB8wclUYt2Q1C9xhMBVAP6w6vnxyrK6knJD0onDE5Gq+FxwUTmoG5L6JY6KwVoja5eNNpIsAigCwMjISM2Bp25r9w0W1xWGq5OLas0rcFXQpBuS+iOOI4HjAK5e9Xw7gBPrVzKzKTPLm1l+eHg4hmbb18obrHpF4dzWHB7IP9B2sdGGzAZM7pysWbikykHphjiOBH4D4FqS1wD4XwD3APjLGLbbcZM7J1F8urjmkuMDfQPYkNmA9y68d2nZ4MZBPHLXIzVn9e0Y2YGJwxOYf3seI1tHLr1pay3b+59719Qk1Nru+t/TTELptFjqBEjuBvAwgAyAx82s4cdXPp+3crncaJWuKR0p6Y0nPSfxxUJJCgGRXpSq+w6IiFuJCQFNWHFD+10ScVGR9RNWqvcpBHSaqpO03wVIyJiAJqy4of3eu1I3JqAJK25ovwuQkBDQhBU3tN8FSEgIaMKKG9rvAiQkBDRhxQ3tdwESMjAoIvFK3cCgiLijEBDxnEJAxHOpCgGVuEanfSbNJKJsOAyVuEanfSZhpObsgEpco9M+81dPnh1QiWt02mcSRmpCQCWu0WmfSRipCQGVuEanfSZhpCYEVOIanfaZhJGagUERCa8nBwZFpDN6NgR6tUimV1+XuJOaYqEoerVIpldfl7jVk2MCvVok06uvS+Ln/ZhArxbJ9OrrErd6MgR6tUimV1+XuNWTIdBqkUy3B92itqfiH+kIM+v6Y2xszDptetoslzMjg5/T083Xz2bNgJVHNtv899rpXyvtRX1d4icAZQv5fmzrzQzgiwBeA7AMIB/297oRAlHlcmvfkNVHLtf8d1t5Y7bTnkgzUUKg3VOErwL4AoDvtrkd51oddGv1tJ0G+SQp2hoTMLOjZvZ6XJ1xqdVBt4mJlQCoWlwMlneiPZG4da1YiGQRQOUzEudIvtqttkMYAoaWgZEcwFXBaMtzc/Nz5Omz9X91bKzW0rk5gJyZadDktgbt9QE4HfE1dNoQktUn9aex68Ku2DQESD4H4E9r/NWEmf172IbMbArAVGWbZQtZyNANQX9OJaw/ydk/QPL6pP40RjJ0NV7TEDCzO9vrjogkWU/WCYhIeG2FAMm/IHkcwJ8DOEjyZyF/daqddjtA/WkuaX1SfxoL3R8nE4hEJDn0dUDEcwoBEc85CwGSXyT5Gsllks5OrZDcRfJ1ksdIfsNVPyp9eZzkyaTUUJC8muTzJI9W/q32JqBPV5B8keQrlT7tS0CfMiR/S/IZ130BAJKzJI+QfDnMqUKXRwLVkuMXXHWAZAbAdwDcBeAGAF8ieYOr/gD4AYBdDttf7yKAr5vZnwG4DcBfOd4/AHAOwB1mdjOAWwDsInmb4z7tBXDUcR/W+5SZ3RKmdsFZCCSk5PhWAMfM7A0zOw/gCQCfd9UZM3sBQIPqxO4ysz+a2UuVP7+L4D/6VY77ZGa2UHk6UHk4G90muR3AZwA85qoP7fJ9TOAqAH9Y9fw4HP8nTyqSowA+CeDXbnty6fD7ZQAnARwyM5d9ehjA3yGYSZsUBuDnJGcq5foNdXTuQFwlxx3EGst0znQdkpsB/ATA18zsHdf9MbMlALeQ/BMAT5G80cy6Po5C8m4AJ81shuTt3W6/gR1mdoLkRwEcIvn7ylFmTR0NgRSUHB8HcPWq59sBnHDUl0QiOYAgAEpm9qTr/qxmZm+R/CWCcRQXg6k7AHyO5G4AVwD4MMlpM9vjoC+XmNmJys+TJJ9C8LW3bgj4/nXgNwCuJXkNyQ0A7gHwH477lBgkCeB7AI6a2bdd9wcASA5XjgBAciOAOwH83kVfzOybZrbdzEYR/N/5hesAILmJ5JbqnwF8Gk0C0uUpwlZLjmNjZhcB/DWAnyEY9Pqxmb3W7X5UkfwRgF8BuI7kcZJfddWXih0Avgzgjsrpppcrn3oufRzA8yR/hyDED5lZIk7NJcTHAPw3yVcAvAjgoJn9tNEvqGxYxHO+fx0Q8Z5CQMRzCgERzykERDynEBDxnEJAxHMKARHP/T/zehTkb4x1fwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "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\")"
      ]
    },
    {
      "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",
      "execution_count": 19,
      "metadata": {},
      "outputs": [],
      "source": [
        "source_proxy_mat = interaction_mat(source_proxies, sources)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(25, 200)"
            ]
          },
          "execution_count": 20,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "source_proxy_mat.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now compute an ID (row or column?):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {},
      "outputs": [],
      "source": [
        "idx, proj = sli.interp_decomp(source_proxy_mat, nproxies)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {},
      "outputs": [],
      "source": [
        "source_skeleton = idx[:nproxies]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "5.2523961564651685e-15"
            ]
          },
          "execution_count": 23,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "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)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the chosen skeleton as well as the proxies:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X2QHHd95/H3d3qen/ZZlrTaXfmCMQG7bKI9w50JZzvE+BxjQyiesmAHSAkjciWe6uocCVMEKzkqBNspsGXFtkBoE5dNIATjBHyxeUgA2ytsAcYOMUYPq6d9fpzn7t/90buSVvug2d3Z7Znt76tqajW9M92/6VV/pvvbv1+3GGNQSvlXwOsGKKW8pSGglM9pCCjlcxoCSvmchoBSPqchoJTPBSsxExE5BIwDNlAyxnRWYr5KqZVXkRCYcrUxZqCC81NKrQI9HFDK56QSPQZF5DfAMGCA+4wxe+Z4zVZgK0Aikdjyqle9atnLVUrN7cCBAwPGmJZyXlupENhojDkuIuuAx4H/ZYz5wXyv7+zsND09PcterlJqbiJyoNzaXEUOB4wxx6d+9gHfAK6oxHyVUitv2SEgIgkRSU3/G7gW+MVy56uUWh2VODtwAfANEZme398ZY/6lAvNVSq2CZYeAMeZl4LIKtEUp5QE9RaiUz2kIKOVzGgJK+ZyGgFI+pyGglM9pCCjlcxoCSvmchoBSPqchoJTPaQgo5XMaAkr5nIaAUj6nIaCUz2kIKOVzGgJK+ZyGgFI+pyGglM9pCCjlcxoCSvmchoBSPqchoJTPaQgo5XMaAkr5nIaAUj6nIaCUz2kIKOVzGgJK+VzFQkBELBF5VkQerdQ8lVIrr5J7AtuBFyo4P6XK090Nzc0g4j6am91pa2V5K6wiISAim4A/AO6vxPyUKlt3N3zgAzA4eGba4CC8//0Lb5jd3bB5MwQC7s9t2848b26GVOrMRp5MutOmf3fzzbOX94EP1GwQiDFm+TMR+Rrwl0AK+KQx5oaFXt/Z2Wl6enqWvVzlc93dcMstYNvzv8ay3N93dMD118Njj8Hhw+7GXYH/+zN0dMChQ5Wd5xKJyAFjTGc5r132noCI3AD0GWMOnOd1W0WkR0R6+vv7l7tY5Xfd3bB168IBAGd+f/gw3Huv+xMqHwAAR45Ufp6roBKHA1cCN4rIIeAh4BoR2X/ui4wxe4wxncaYzpaWlgosVvna9u2QyXjdipna271uwZIsOwSMMbcZYzYZYzYD7waeMMa8d9ktU2ou00W5s4/Jq8XERE3WBYJeN0Cpsk0fAlTbHsC0wUG3fQBdXd62ZREqUhhcLC0MqiXZvPnMMX01q4IC4aoWBpVaNbVSeKuVdk7REFC1o7HR6xaUp8YKhBoCSlXarl1et2BRNARU7Rga8roF5amhoiBoCKhaUmO72bVCQ0DVjlrYzbYsr1uwaBoCqnbUwm72VVd53YJF0xBQtaWpyesWLOyll7xuwaJpCKjacvfd7pDealVjfQRAQ0DVmq4uaGjwuhXzq8HipYaAqj3VfKqwFoqX59AQULWnBr9tq5mGgKo9u3ZBPO51K+a2Y4fXLVg0DQFVe7q6YM8ed7RetdHCoFKrpKvLHa6bSHjdkplq8FBFQ0DVru5umJxc/PuamiAcrnx7wL2YaY3RKwup2jR9laGlWMlLkz322MrNe4XonoCqTTt2VOdlxrQmoNQqqdbLjGlNQCkfi8e1s5BSvtXR4Z62rIWRjufQwqCqTdO3F6sGH/4w3HOP161YMt0TULVpqWcGVsLu3TV505FpGgKqNlXTN68x7o1RazQINARUbaq2Dc623b2TamtXGTQEVG2qxoE6mUx1tus8NARUbarWTjnV2q4FaAio2lStnXKqtV0LWHYIiEhURJ4WkYMi8ryIfKYSDVNqQdV4TQEfdxbKA9cYYy4DLgeuE5HXV2C+Ss1v+poC1XL1Ycuq2c5Cyw4B45qYehqaeqz+/c6V/3R1wcAA7N+/+DBoaqrstQgcpyYDACpUExARS0SeA/qAx40xT83xmq0i0iMiPf39/ZVYrFKu6TAQKe/1HR3u6ycm3ACpxGFFDdYCplUkBIwxtjHmcmATcIWIXDLHa/YYYzqNMZ0tLS2VWKxSM5WzIZ573L7QYUUoVN7FR2q0FjCtomcHjDEjwPeA6yo5X6XKMl+xcPpmJfMN8jn7sKKjw92j6OiAvXvhwQfPTGtqOhMW0/ccrOGBQ9PEmOUdvotIC1A0xoyISAz4LvA5Y8yj872ns7PT9PT0LGu5Ss2pu9vtsHPkiLtnsGtXTW+gSyUiB4wxneW8thKjCDcAXxERC3fP4uGFAkCpFdXV5cuNfjmWHQLGmJ8Br61AW5RSHtAeg0r5nIaAUj6nIaCUz2kIKOVzGgJK+ZyGgFI+pyGglM9pCCjlcxoCSvmchoBSPqchoJTPaQgo5XMaAkr5nIaAUj6nIaCUz2kIKOVzGgJK+ZyGgFI+pyGglM9pCCjlcxoCSvmchoBSPqchoJTPaQgo5XMaAkr5nIaAUj6nIaCUzy07BESkTUSeFJEXROR5EdleiYbVku5u2LzZvQP25s3uc6VqRSXuSlwCPmGM+amIpIADIvK4MeaXFZh31evuhq1bIZNxnx8+7D4HvTmuqg3L3hMwxpwwxvx06t/jwAtA63Ln64WlfKPv2HEmAKZlMu70lVieUhVnjKnYA9gMHAHSc/xuK9AD9LS3t5tqs3+/MfG4MXDmEY+70xciMvM90w+RlVmeUuUAekyZ2624r18+EUkC3wd2GWO+vtBrOzs7TU9PT0WWWymbN7u78ufq6IBDh6rnfUqVQ0QOGGM6y3ltRc4OiEgI+Aeg+3wBUK2OHFnc9Gm7dkE8PnNaPO5OX4nlKVVplTg7IMADwAvGmC8sv0neaG9f3PRpXV2wZ4/7DS7i/tyz5/xFwaUuT6lKq8SewJXA+4BrROS5qcf1FZjvsiy26LbUb3RwN/hDh8Bx3J/lnBVY6vK0mKgqrtziQSUfW7ZsWZlqyJSlFt327zemo8Mt6nV0rHyRbrHL02KiKhdeFAYXY6ULg2u16LZWP5eqvFUvDFabtVp0W6ufS3lrTYbAWi26rdXPpby1JkNgOUW+arZWP5fy1poMgaWetqt2a/VzKW+tycKgUn7n+8KgUqp8NRUC2lFm8XSdqfOpxPUEVoWO2188XWeqHDVTE9COMoun68y/1mRNQDvKLJ6uM1WOmgkB7SizeLrOVDlqJgS0o8zi6TpT5aiZENCOMoun60yVo2YKg0qp8q3JwqBSamVoCCjlcxoCSvlc1YSAdm/1hq53VRXdhrV7qzd0vSuokrMD2r3VG7re166aOzug3Vu9oetdQZWEgHZv9YaudwVVEgLavdUbut4VVEkIaPdWb+h6V1AlhUGlVGXVXGFQKeWdSt2a/EER6RORX1Rifl7Z++xeNn1hE4HPBNj0hU3sfXav101SasVVak/gy8B1FZqXJ+4/cD/bvr2NY+PHMBiOjR9j22Pb5gyCol1kNDfKYGaQ3c/spu3ONgKfCdB+Zzv7Du6bN0zmm77v4D7a72yfMQ+lVkvFagIishl41BhzyfleW001gb3P7mXnkzs5Pn58zt+3plrp/Xjv6efTARAMBHn4+YfZ/p3t5Eq5078PShARoegUT0+LBqO869Xv4qHnHyJv509Pj1gR3vnqd/K1F75GtpQ9PT0WjLH7ht3cfNnNlfyoykcWUxNYtRAQka3AVoD29vYth+fqqrbK9j67l22PbZuxEZ9LEL781i+z84md9I71sj65no+//uO8+5J387oHXjdveJwrIAEc45Q9vS3dxpGPaa8dtTRVWRg0xuwxxnQaYzotq6UqBq186slPLRgAAHWROm599FaOjh3FYDgxcYJPfe9TPPL8I5wYP1H2suba0Bea3jvWO+f01aIDi/zDkwFEhw+D45z5t1eDVs73LR6xIojIjF11gFwpx+d//HkuSFzAycmTZS1rsXsCm9KbyprvStCBRf7iySlC55z/85kM7Nix+u3YmNo47+9a4i38xdV/wUhuZM7fn5g4wXsvfS9hKzxjuiUWQZmZrRErwpv/y5uJWJEZ08NWmLde/Faiweis+U8WJj0rEO7YcSYApnn1N1Irr1KnCP8e+DFwsYj0isgHFzsPLwatfPbqz87aYAGCgSAf2fIRfu+3fo8LkhfM+d7meDOXtlzKzZfeTFOsCXCDY+cbdrL9iu2np62Lr+O2N9zG7f/jdj56xUdZF1/nvj/WzG3//Tb+6k1/xR1X30FdpG7G/IdyQ9z66K2eBIEOLPIXT3oMinQamHl2wKvhqw2fa5jz274p1sQjb3+E7x/+Pp/70efI2WdqBxErwsde9zGubL8S27GZLEySCCXIFrMkIgkKxQKhYIh4KE5rupWiU2RgcoCR/AihQAjbsdmY2khDvAEMjOZHuf7vrqc/0z+rHV4UCHWIce1bTGHQk5pAIDDzkMDLQSujudE5pw9mB3HE4arNVxG0gtx74F5OjJ/ggsQF3HL5Lad37y3Lone0FxEBIGpFKdpFxrJjZEtZHByCBMmVckStKOlIGsdxyJfyjOfGGc4PUygWGMgMzNkOLwqEu3bNrAmADixayzwJgY4ONwSOHHGHre7a5V3BqTXdOueG1hJvoVgsMlQY4g1tb+DqzVcTCUZIhVMYDH0TfRSdIhOlCVLRFBO5Cfdb3xQxxmBjEzIhjowcwbZtgpa7qqPBKBErwnhhnJMTJ6mP1pMtZWmKNTGQnR0EXhQIp/8WO3ZUx99IrSxPQqCxEaqhr9C+g/uYyE/Mmh4KhHjLK97CYG6QaDBKMpxkPD9OvpSnf7KfdCjNaH6UsIQRSygWi4zlx0AgW8xSKBWIh+OErTDHJ45jbEMoFCIRStA32YcgBCRAyS5BAJKRJF2XdnHvgXsp2IXT7YgFY9xxzR2ruUpO6+rSjd4vquIag17Yd3Aftz5666zTf8lwkvdd8j4ua7mMoARJh9MMZNzj+dZkK/FQnKHsEBYWeSePbdvEgjE2JDdwavIUQQmyIb2BXDHHaGGUxmgjpyZOkc/nSYaTDIwP0JfroyHSQKaYoeSUaEw0ctXmqwDo/nk3g9lBNiQ3cMc1d2ivQbXifBsCO5/YOSsAAOLBODf99k38ZvA3DOWHcIyDFbBoTbYSC8XIlXI0xBpIR9IMZgYJW2EioQi5Yo5sKUvBLlBySiQiCUZzo+RKOYxjGMuNYRx3j2AyN0nABMDAyfGTxMIxkuEkxjmrSCtgBaxVXCPKr3w7lHi+glt/pp9sLksilKAh4m7sAFk7y1h+jIJdYDQ3ynBumLH8GIlIgsZoI/FwnLa6NhLhBKPZUQISoCXZgjGGTXWbSMVSFJ0ih4cOEw6GsYIWsUgMIwbbsXnk+UfY/exuBrIDGAzHx497dopQ+YtvQ2C+gltLvIVYJMaFDRcSD8aZKE4QkAC5Qo5MMcNgZpCSU8JxHCJWhJHMCP2Zfop2kfpoPR11HaRjaXLFHGO5MeLhOOtS62iMNZIIJ2hINFAfrScdStMYa8QxDiLCt1761ox6AEC2lGXnEztXY3UoH/NtCNxxzR3EgrEZ0yJWhFsuu4X6aD2tda00J5pJhVLYxmasMEbIcot74UCYVCRFU6IJxzjki3nqInUU7AKTxUle2fhKOuo7aIo30VbXhiCII0zkJwgHwuTtPKFgiGwx6z4v5avqFKHyF9+GwM2X3czuG3bTlm4D3O6+eTvP/p/v59+P/juRYISwFSYcCpMOp2mJtxC2wthiEwvHSEfSxENx0rE0EnCHDoesEE3xJiQgiAipSIoAARqiDVzYfCHRiNuHIFfMMTQ5hG1s6qJ1JMNJWuItc7bTyzEEyh88CYEDB6pjZNrNl918eo/ANjbgjgn4syf+jK8e/CpFp0hYwoSsEEErSNgK0xhudF83foLR3CgWFslQkmgw6u7uxxoIBoKkwinioTi2sTkxeYJ8Kc/6+HqiVpRUJEV9rJ6oFcWyLNrSbWx/3fZZYwi8PEU4bds2CAbdC5EGg+5ztbZ42m04Hvf+6rbtd7ZzdOzorOnNsWYe+sOHCFkhoqEoJafEeHacw2OHiVpRmuPNhK0wRgytqVZaEi1kS1m3o5CxKdgFekd7Gc4OY9s2GTtDIpigRIliqUhDvAHbsQkFQ0SsCLlSjqePPc3e5/ZyYuIErelWdl2zy9NThNu2wb33zp7+4Q/DPfesfntU+Ty5qMhinD12wOv+6IHPBDDMXgeC8JMP/ATHcRjJj2AFLGzHZjw/TsEpEAwEiVpRAoEA65LrWBdfR6aYoX+yn2wxSyQUoXe0l6AEaU40M5gdpD5Sz2RxkmwhS2tdKwj0TfSxIbmBklMiZIWIh+Nc1HgRdbG6OVq7uoJBsO3Z0y0LSqXVb48qX9WPHTib1yPTNqU3zb0nEG/m0PAhik7RHfCT3kjBLtCcaAYgX8xjWRbpcJp80e1JOFGYYDg7zHh+3D3VWMwSC8UIWSFiwRjGGKLBKPFQnHwpz8mJk0QCESYKE8TCMX545Id8seeLnJo4xab0Js87C80VAAtNV7XJ88Kg17e8uv2Nt88a5x+xItz6O7diG5uQFQJxry14dOwo4/lxSk6JExMnGM4MM5Qb4qXBlzg6dpR8KU/ezjOSH8ExDi3JFuLBOC8NvQS43YOnOxfFQ3GKpSJZO8tobpR/fflf+ewPP8vJiZMYDEfHjnreT8Cap6/SfNNVbfI0BLwemZYpZPjd9t/l02/8NOvi6xCEplgTn3z9J7my/UrWJ9e7RUPbZjAziCD8euDXnBo7Rd94H73jvfSO9FK0iySCCWxj8/LwyxhjKNklTo2dQgJCY6yRgATYmNpILBRjU2oT7fXttNa30ppqpTXdyt8++7czLkIK3vcTmL6aULnTVW3yLASq4ZZXA5kBvvniN7n76bvpz/TTFG/imo5rSEVSHBs7RqaUwcbGsiyGckNuIGBzbPwYIStEJBAhEU5QH6/HwWE0N4oxhoJTcL/1Q2GigSgITOQmCAQCbK7fzMb0RkpOiaJd5Nj4MY6PH5/zWgLgbT+Be+5xi4DT3/yWpUXBtciTmsCWLdUxivChXzzE7d+//fQ38EBmgH/81T8SDoTp3NDJSG6EVDhFXbQOEWGyMEksHCMWjJEKp8gWs4RCIeoj9RwdOUrRKZItut2Lm+PNJENJDIbGSCN1sTra6trIFDOM58Y5MnqEZNg9tTgwOUBDtIHh3PCsNnrdT+Cee3SjX+s8rwl46a6n7pq1C150ivzzy/98urNQySlRckpYWNjGJkDAva+AXSQVSWEcQzgQZjQ3SrFUJB1JYxdtcsUc6UiahngDDYkGXtH4Crf6H4pzbOIYhZI70ChTzPDsqWcpObPL7dXQT0CtfZ6fHfDSyYm5rxQ8lBtyrwAUcBjJjDAwOUC2kCUWiVEoFcgUMgxkBmiONVMfqyccCBMNRmmMNVJv6mlNtzKZnyRXylEv9VzUcBH18frT8y+UCtRF6xjLj/Gj3h/xxZ4vzho30BBt4K7r7tKhxGrF+ToE5ruqUFOsiUjQvfrPxtRGJguTNMeb+c+h/yRfyBMOhUlGktiOjYhwfOI4DfEG0tE0iVCCnJ1jMjxJ3s7Tmm4lHAwTCoQYzY0ylB1CEIKBIPWxevb9bN+sAAD3ugYaAGo1+PpwYMcbdszqqhsOhLnxohvJlrJYYlEyJUbzowxODjKWd0cF1oXrSAQTNCQbaK9rpz5az+a6zdRF6hCEoewQE4UJjGOwxDp96zJjDCJCPBjnxYEX6R/v14FDynO+3hN4x2veQUAC/PkP/pzj48dpibfw1le9lYvqLuLw6GGaYk3Ewm4nn2wxS320nrpIHUYMGPeioolwAskIk8VJJguTBK0gTdEmbGOfHidQdIoYDMlIEgwMZAcISpATEydoiDYwlBua1TavC4LKP3wdAsFAkHe+5p1c+1vXMpobRRAOjRxyv7mzoxgMuWKOfCmPWELIDjGSGyEWipEpZSjZJSLBCM3xZoKBoBsEk5OsS65zBxBF4pyaPMVYboz1qfWkIikyhQwnxk+QjqRJR9N88PIPcvczd1fNtQWV//g6BOKhOKO5Ub568Kt8qedL9E320RBtoOvSLq5YfwWnsqdoCbUQDoQxGGzLva5AIpSgPlLP8OQwJ8dOkoqkSAfTrE+tZ3BikOHMMCWnxERxgqgVZSw3RtF271I8lB2iJd5CwSng4HDDK28gbIX5m2f+hvHCOMCcdyRSaqV4UhMYGqqOm12GrBDf+tW32PVvuzg1eQqDYSg3xH0H7uPRXz+KGHFvLmJPMl4YpzHRyMVNF9OYaHSvJBSvpzHe6N6OPOAW+5oTzdiOzUR+gkKpwFDG3dWPhWIcGznGweMH6R3rpW+yj6AEiYTcLstnn6oczg173mVYb0jqH56MIrSsTuM4Z3oLeTmkuO3OtjmLcOlwmo/+14/Slm4jFU0xnBtmIDPAhY0Xkg6lCchUfgbcYccAJdv99h/ODp8+719ySjTGGhnNjzKcHWY4M0zcirO+bj3RcJSACfChxz7EYHZwdts8uj35uTckBW//Rmrxamoo8TSvhhTPN5QY4O5r7yYZTvL0saf52otfYzA7SGO0kT+65I+48ZU3gkCumKM53oxjHCZLbnHQcRyMMRCAYsm9Gcnx8eMEA0GSoSQFCiStJIlIAsc4vOvr75p3OLPz6blvXb6S9DZktW/VhxKLyHXA3YAF3G+M+b+LnYdXQ4o3pDbMe4vy7d/dTkhCGAwl436zD+WG2PPsHuKhONe/4nqe6H2CB597kL5MH83xZj702g9Rcko8ePBB+jP9NEYbecsr3sJr1r2Gnxz/Cd/59XcYyY/QGG3kbRe/jRsvvnHeW5w3RBtW9LPPR29I6i/L3hMQEQv4FfD7QC/wDPAeY8wv539P9ewJ/PE3/piv/Owri37fuvg6/uTyP+Gvn/rrGcfzlliIyIxuwKFAiM71nfSc7KHoFE9PD1thPvG6T9AQa2DHkztm/A4gJCHuv+n+Ve80pHsCtW8xewKVKAxeAbxkjHnZGFMAHgJuWnCh5yzVqyHF+w7u4+FfPryk9/Zn+nng4AOzxh7Yxp41DqDoFHnq+FOzNvKCXeDBgw9y48U3kgglZi2jaIqeDCXetcv9m5zN62HfauVUIgRagbMvzdM7NW1eHR3uQ8TbIcXz3YWoHOuT6+mb7Cv79Q5zH9v3TfZRF61jND/33ZG96DnY1eX+Tarhb6RWXiVCQOaYNusYQ0S2ikiPiPTYdj+HDrl3Jj50yLv/XEvdwGLBGLe/8XY2pDaU/Z7TZxPOcfpCI/P0EPSq52BXF1XxN1IrrxIh0Au0nfV8EzCr0maM2WOM6TTGdLa0zH2N/dVW7gYmCPXRegShNdXK53//87zn0vdwx9Wzb2ASlCChQGjGtGgwyttf9XaiVnTW9NuuvI14KD7nzVC056BaDZUIgWeAi0TkQhEJA+8G/qkC811xc214IQnNOD6vj9Rzz/X30PfJPpxPO/R+vJdtV2yjLlrH+1/7/tM3MBGEtnQbD9z0APffeP+MaffdcB/db+/mzjffycbkRgRhY3Ijd157J1s7txKyQjNuhjL9vt037NaRhGrFVaSfgIhcD9yFe4rwQWPMgiWkzs5O01MNlxbCLQ7ufGInvWO9VXGFX6Uqoeo7C1VTCCi1Fq32KUKlVA2rmhDQASve0PWuqmIo8bkDVg4fPnNtez01tXJ0vSuokpqAdlP1hq73tavmagI6YMUbut4VVEkIzHc/Qq/vU7jW6XpXUCUhoANWvKHrXUGVhIAOWPGGrncFVVIYVEpVVs0VBpVS3tEQUMrnNASU8rmaCgHt4rp4us7U+VRFt+FyaBfXxdN1pspRM2cHtIvr4uk68681eXZAu7gunq4zVY6aCQHt4rp4us5UOWomBLSL6+LpOlPlqJkQ0C6ui6frTJWjZgqDSqnyrcnCoFJqZazZEFirnWTW6udS3qmZzkKLsVY7yazVz6W8tSZrAmu1k8xa/Vyq8nxfE1irnWTW6udS3lqTIbBWO8ms1c+lvLUmQ2CpnWRWu+i22OVp5x+1Iowxq/7YsmWLWWn79xvT0WGMiPtz//7zvz4eNwbOPOLx879vOe1byvIW+7mUPwE9psztcVkbM/AO4HnAATrLfd9qhMBidXTM3CCnHx0d53/vUjbM5SxPqfNZTAgs9xThL4A/BO5b5nw8t9Si21JP22mRT1WLZdUEjDEvGGP+o1KN8dJSi247dpwJgGmZjDt9JZanVKWtWmchEdkKTH1HkheRX6zWssvQDM0OtHeAnBWMxjl8+MhhkYGh+d+6ZctcUw8fBpEDBxZYZOMCywsAA4v8DCutmepqk7ZnYReX+8LzhoCI/D9g/Ry/2mGM+Wa5CzLG7AH2TM2zx5TZkWE1uO3pr7L2VM/6geprk7ZnYSJSdm+884aAMeZNy2uOUqqarcl+Akqp8i0rBETkbSLSC/w34Nsi8p0y37pnOctdAdqe86u2Nml7FlZ2ezwZQKSUqh56OKCUz2kIKOVznoWAiLxDRJ4XEUdEPDu1IiLXich/iMhLIvJ/vGrHVFseFJG+aulDISJtIvKkiLww9bfaXgVtiorI0yJycKpNn6mCNlki8qyIPOp1WwBE5JCI/FxEnivnVKGXewLTXY5/4FUDRMQCvgT8T+DVwHtE5NVetQf4MnCdh8s/Vwn4hDHmt4HXAx/xeP0A5IFrjDGXAZcD14nI6z1u03bgBY/bcK6rjTGXl9N3wbMQqJIux1cALxljXjbGFICHgJu8aowx5gfAAr0TV5cx5oQx5qdT/x7H/Y/e6nGbjDFmYuppaOrhWXVbRDYBfwDc71UblsvvNYFW4OhZz3vx+D95tRKRzcBrgae8bcnp3e/ngD7gcWOMl226C/jfuCNpq4UBvisiB6a66y9oRccOVKrL8QqSOabpOdNziEgS+Afgo8aYMa/bY4yxgctFpB74hohcYoxZ9TqKiNwA9BljDojIVau9/AVcaYw5LiLrgMdF5MWpvcw5rWgI1ECX416g7aznm4DjHrWlKolICDcAuo0xX/e6PWczxoyIyPdw6yheFFOvBG4Ukese2NaEAAAA8UlEQVSBKJAWkf3GmPd60JbTjDHHp372icg3cA975w0Bvx8OPANcJCIXikgYeDfwTx63qWqIiAAPAC8YY77gdXsARKRlag8AEYkBbwJe9KItxpjbjDGbjDGbcf/vPOF1AIhIQkRS0/8GruU8AenlKcKldjmuGGNMCfhT4Du4Ra+HjTHPr3Y7ponI3wM/Bi4WkV4R+aBXbZlyJfA+4Jqp003PTX3reWkD8KSI/Aw3xB83xlTFqbkqcQHwbyJyEHga+LYx5l8WeoN2G1bK5/x+OKCU72kIKOVzGgJK+ZyGgFI+pyGglM9pCCjlcxoCSvnc/wdQ3gLfx5Tj+QAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "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\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check that it works for 'our' targets:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "5.966880531704515e-09"
            ]
          },
          "execution_count": 25,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "imat_error = (\n",
        "    interaction_mat(targets, sources[:, source_skeleton]) @ P\n",
        "    -\n",
        "    interaction_mat(targets, sources))\n",
        "\n",
        "la.norm(imat_error, 2)"
      ]
    },
    {
      "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."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Hierarchical Skeletonization"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7fb821315fd0>]"
            ]
          },
          "execution_count": 18,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABKrklEQVR4nO3dfXRU9Z0/8PfNJJkEzARQyeMItiDPCIorAaNxoeVUugeaH5Qiij2LntWFLdHtA+xxq9VqXF0r9GxFq0dxSykKG6RFBSkVioJbEdgFba0PAQImwbVhJgkxD3e+vz++vZmH3JncO0/33pn365w5Se5878y988nc+7nf+31QhBACRERERBbJsXoDiIiIKLsxGSEiIiJLMRkhIiIiSzEZISIiIksxGSEiIiJLMRkhIiIiSzEZISIiIksxGSEiIiJL5Vq9AUYEAgF8+umnKCoqgqIoVm8OERERGSCEQHt7O8rLy5GTE73+wxHJyKeffgqv12v1ZhAREVEcmpqaUFlZGfV5RyQjRUVFAOTOeDwei7eGiIiIjPD7/fB6vf3n8WgckYxot2Y8Hg+TESIiIocZrIkFG7ASERGRpZiMEBERkaWYjBAREZGlmIwQERGRpZiMEBERkaWYjBAREZGlmIwQERGRpZiMEBERkaUcMegZEWUmVQUOHACam4GyMqC6GnC5rN4qIko30zUjZ8+exS233IKLL74YhYWFmDJlCg4fPhxznX379uGqq66C2+3GmDFjsHHjxni3l4gyREMDMHo0cOONwM03y5+jR8vlRJRdTCUjbW1tmD17NvLy8vDaa6/h/fffx+OPP47hw4dHXaexsRHz58/HjTfeiGPHjqGurg633347du/enfDGE5EzNTQAixYBZ86ELz97Vi5nQkKUXRQhhDBaeM2aNXjrrbdw4MABw2/wgx/8AK+88gpOnDjRv+xb3/oWzp8/j127dhl6Db/fj+LiYvh8Ps5NQ+RwqiprQCITEY2iAJWVQGMjb9kQOZ3R87epmpFf//rXmDFjBhYvXoyRI0di+vTpeOaZZ2Kuc+jQIcydOzds2bx583Do0KGo63R3d8Pv94c9iCgzHDgQPREBACGApiZZjoiyg6lk5JNPPsGGDRswduxY7N69G3fddRe+853v4IUXXoi6TktLC0pKSsKWlZSUwO/3o6urS3ed+vp6FBcX9z+8Xq+ZzSQiG2tuTm45InI+U71pAoEAZsyYgYcffhgAMH36dJw4cQJPPfUUbrvttqRt1Nq1a3HPPff0/+33+5mQUNbI9B4mZWXJLecEmR5TokSZqhkpKyvDxIkTw5ZNmDABp0+fjrpOaWkpWltbw5a1trbC4/GgsLBQdx232w2PxxP2IMoG2dDDpLpatglRFP3nFQXwemW5TJANMSVKlKlkZPbs2fjggw/Clv35z3/GqFGjoq5TVVWFvXv3hi3bs2cPqqqqzLw1UcbLlh4mLhewfr38PTIh0f5ety4zag6yJaZECRMm/OEPfxC5ubnioYceEh9++KH45S9/KYYMGSI2bdrUX2bNmjXi1ltv7f/7k08+EUOGDBHf+973xB//+Efxs5/9TLhcLrFr1y7D7+vz+QQA4fP5zGwukWP09QlRWSmEbL458KEoQni9slym+K//GrjPXq9cngmyMaZEkYyev03VjFxzzTXYvn07fvWrX2Hy5Ml48MEHsW7dOixbtqy/THNzc9htm8svvxyvvPIK9uzZgyuvvBKPP/44nn32WcybNy9Z+RSR42VjD5PaWuDkSeCNN4DNm+XPxka5PBNkY0yJ4mV6OPivf/3r+PrXvx71eb3RVWtqanD06FGzb0WUNbK1h4nLBdTUWL0VqZGtMSWKByfKI7KBbOxhkukYUyLjmIwQ2UC29TDJBowpkXFMRohsIJt6mGQLxpTIOCYjRDZRWwts2wZUVIQvr6yUyzOlYWc2YUyJjDE1UZ5VOFEeZROO1pl5GFPKVkbP36Z70xBRamVyD5NsxZgSxcbbNERERGQpJiNERERkKSYjREREZCkmI0RERGQpJiNERERkKfamIcdid0kia/E7SMnCZIQcqaEBWL06fFbUyko54iUHkiJKPX4HKZl4m4Ycp6EBWLRo4PTsZ8/K5Q0N1mwXUbbgd5CSjSOwkqOoKjB69MCDoEZR5NVZYyOri4lSgd9BMsPo+Zs1I+QoBw5EPwgCgBBAU5MsR0TJx+8gpQKTEXKU5ubkliMic/gdpFRgMkKOUlaW3HJEZA6/g5QKTEbIUaqr5f1oRdF/XlEAr1eWI6Lk43eQUoHJCDmKyyW7DgIDD4ba3+vWseEcUarwO0ipwGSEHKe2Fti2DaioCF9eWSmXc4wDotTid5CSjV17ybE4+iORtfgdpMEYPX9zBFZyLJcLqKmxeiuIshe/g5QsvE1DRERElmIyQkRERJZiMkJERESWYjJCRERElmIyQkRERJZiMkJERESWYjJCRERElmIyQkRERJYylYzcf//9UBQl7DF+/Pio5Tdu3DigfEFBQcIbTURERJnD9AiskyZNwm9/+9vgC+TGfgmPx4MPPvig/28l2lSPRERElJVMJyO5ubkoLS01XF5RFFPliYiIKLuYbjPy4Ycfory8HF/60pewbNkynD59Omb5jo4OjBo1Cl6vFwsWLMB777036Ht0d3fD7/eHPYiIiCgzmUpGrr32WmzcuBG7du3Chg0b0NjYiOrqarS3t+uWHzduHJ577jns2LEDmzZtQiAQwKxZs3DmzJmY71NfX4/i4uL+h9frNbOZRERE5CCKEELEu/L58+cxatQo/OQnP8GKFSsGLd/b24sJEyZg6dKlePDBB6OW6+7uRnd3d//ffr8fXq930CmIiYiIyD78fj+Ki4sHPX+bbjMSatiwYbjiiivw0UcfGSqfl5eH6dOnD1re7XbD7XYnsmlERETkEAmNM9LR0YGPP/4YZWVlhsqrqorjx48bLk+pp6rAvn3Ar34lf6qq1VtERJR6PPbZi6lk5Lvf/S7279+PkydP4uDBg/jGN74Bl8uFpUuXAgCWL1+OtWvX9pd/4IEH8Prrr+OTTz7BkSNHcMstt+DUqVO4/fbbk7sXFJeGBmD0aODGG4Gbb5Y/R4+Wy4mIMhWPffZj6jbNmTNnsHTpUnz++ee49NJLcd111+Htt9/GpZdeCgA4ffo0cnKC+U1bWxvuuOMOtLS0YPjw4bj66qtx8OBBTJw4Mbl7QaY1NACLFgGRLYbOnpXLt20Damut2TYiolThsc+eEmrAmi5GG8CQMaoqrwKidWpSFKCyEmhsBFyutG4aEVHK8NiXfkbP35ybJgsdOBD9ywjIK4amJlmOiChT8NhnX0xGslBzc3LLERE5AY999sVkJAsZ7czETk9ElEl47LMvJiNZqLpa3heNNmehogBeryxHRJQpeOyzLyYjWcjlAtavl79Hfim1v9etYwMuIsosPPbZF5ORLFVbK7uwVVSEL6+sZNc2IspcPPbZE7v2ZjlVlS3Hm5vlfdLqal4VEFHm47EvPdIyNw05n8sF1NRYvRVEROnFY5+98DYNERERWYrJCBEREVmKyQgRERFZiskIERERWYrJCBEREVmKyQgRERFZiskIERERWYrJCBEREVmKyQgRERFZiskIERERWYrJCBEREVmKyQgRERFZiskIERERWYrJCBEREVmKyQgRERFZiskIERERWYrJCBEREVmKyQgRERFZiskIERERWYrJCBEREVmKyQgRERFZiskIERERWSrX6g3INqoKHDgANDcDZWVAdTXgclm9VUREZKVsPzeYqhm5//77oShK2GP8+PEx19m6dSvGjx+PgoICTJkyBa+++mpCG+xkDQ3A6NHAjTcCN98sf44eLZcTEVF24rkhjts0kyZNQnNzc//jzTffjFr24MGDWLp0KVasWIGjR49i4cKFWLhwIU6cOJHQRjtRQwOwaBFw5kz48rNn5fJs+qdzJFUF9u4FfvEL+VNVrd4i+7LzZ6WqQG+vvbbJjuwcwwzDc4OkCCGE0cL3338/Xn75ZRw7dsxQ+SVLlqCzsxM7d+7sXzZz5kxMmzYNTz31lOGN9Pv9KC4uhs/ng8fjMbyeXaiqzHIj/9k0igJUVgKNjdlVLZc0qgrs2wd8+ilQXg7U1CT3g9y6Fairk6+vKS8H1q0DFi+2dtvsQlWBQEAeOe+5J/7PKhAAcnKS/xn19QFdXUB3d/A93G6gsBDIHeRudTbEMPSzb2iI//+dTMmGc4Ph87cw4b777hNDhgwRZWVl4vLLLxc333yzOHXqVNTyXq9XPPHEE2HLfvjDH4qpU6fGfJ8vvvhC+Hy+/kdTU5MAIHw+n5nNtY033hACGPzxxhtWb6kDvfSSEOXl4R9keblcnqzXV5SBwVIU+Yj1PqneNjvo7RXC7xfis8+EeOaZ+D6r0NdobZU//X65PFnb+Pnn8rXb2oTw+eTP1la5PNb7ZHoMIz/7Z5+N//+dTMuGc4PP5zN0/jZ1m+baa6/Fxo0bsWvXLmzYsAGNjY2orq5Ge3u7bvmWlhaUlJSELSspKUFLS0vM96mvr0dxcXH/w+v1mtlM22luTm45+qutW4ElS8Kv4AD5QS5ZIp/XY7SqXlXlFaJe5aG27O679V8n3m1zkr4+wO+XNQ6KAvzwh+Y/q9DXyM0FCgrkz64uubyvT/+9zdxu6eqSr1NYKC8vFUX+LCwM1pjoyfQYRn72eXnAv/5rfP/vFBeeG4JMJSNf+9rXsHjxYkydOhXz5s3Dq6++ivPnz+Oll15K6katXbsWPp+v/9HU1JTU10+3srLkliPElyj09QHt7UBbW/DR3h79hKdVzUcjhLyxu29f4tvmRKEn+UOHYh8xo31WZhMFszFUVXlrJj9f//n8fPl8ZCyyIYaRn/1bb8UXQ4obzw1BCY0zMmzYMFxxxRX46KOPdJ8vLS1Fa2tr2LLW1laUlpbGfF232w2PxxP2cLLqannfT1H0n1cUwOuV5cggs4lCPFfgsV4/Vrl4k5h4WdHYMPIkb/TSLfRzMZsoxBPDQCDYFkJPTo6MRyAQvjzdMQTS27hW77OPJ4aUEJ4bghJKRjo6OvDxxx+jLEraVlVVhb1794Yt27NnD6qqqhJ5W8dxuYD16+Xvkf902t/r1jm3gZIlzCYK8VTVl5cbe4/IcvEmMfHYuhW47DJg7lxg+XL587LLUn8LIfIkb/TSLfSzMpsoxBPDnBz5iEw2QrdBUQZuQzpjaLa2Jxn0Pvt4YkgJ4bkhyFQy8t3vfhf79+/HyZMncfDgQXzjG9+Ay+XC0qVLAQDLly/H2rVr+8uvXr0au3btwuOPP44//elPuP/++3H48GGsWrUquXvhALW1wLZtQEVF+PLKSrm8ttaa7XIsM4lCvFX1NTVy/ViXLRUVsly825YIK9s0RJ7kr7tOnszMfFZCBGsD9IQmCvHG0OWSvWZ6evTX6+mRz0ce7dMVw3jbzCRKL0kzG0N2/00KnhskU8nImTNnsHTpUowbNw7f/OY3cfHFF+Ptt9/GpZdeCgA4ffo0mkOq+mbNmoXNmzfj5z//Oa688kps27YNL7/8MiZPnpzcvXCI2lrg5EngjTeAzZvlz8bG7PlnSyoziUK8VfUul7ws0V4v8vUB4IknBp7I4k1izLC6TUPkSd7lAh55RP4+2Gel1QRoJ+HPPtOvCQhNFOKNIRDsvtvVJT8PLQnSEoDCwoHrpCOGQPyNaxOll6SZiaFVNXIZiucGk+OMWMXp44xQimg1A0D4SVk7cL74ohwXQVVl1Xdurn59p6rKA//w4frP640zUlEhD8zRxl0wum3x2rtXngAG89vfAnPmxP8+sWhX9X19smZCG6Ni7VogtMdc6GcVuU4gAPh88qRbUAAMGyZfp6dHxsvjkT8TjWHoOCNCyDgMNs5IqmOY6D4lSi9+gYC8HP/hD8PbkITGUPtcIk8dyfpcKKMYPX8zGSFnM5ootLfLk5HeVbC2vKgo+vvEM/BVPEmMUb/4hbwiHcx//idw662JvVcseif53FzgD38AWlsHflZ6cdBew+eTCUlRkX6ikGgMAfMDq6Uyhr29MhkpKNCvgRFCfq7Dhslut6kQLUnLzwfefHPg/7uqyhqQaG1lFEWWP3UqOxo60KCYjFD2MJIoRLsKjLwCt2Lb4mGHmpFQRk7yg9UE9PTIx8UX67cNybQYWl0zEvleRpI0u/3fke0ZPX9z1l5yPpdr8AOfdrKKvAosLDQ2JHgqty0eWpuG5mb9diPaFWqibRqMcrkGP2EO1u4jLy/YaFVPpsVQa7cRrbanpyfYliTVjMQPSG8vI8oqTEbI+Yxeuebmymr8IUNSNwdKumiNa5cskSdkvTYNeo1rrRTag0Nvu3p7ZU1ArMraTIohIJON3l6ZkOjV9uglKVZKVy8jyjoJjTNCZLl4WvW7XPIq3MknMUC2WXjxxYHjQ5SX27MRYbRutlrvmnPngl1aBxtnI1NiqNX2aL1nuruDvWtSddspEenqZURZh21GyLnYql9y0qyy8famyQapnLU4mVLdy4gyChuwUmZjq37nCu3B4ffLxEOrHQhNPIz2kKH0S2UvI8oobMBKmc3M3CFs1W8vWrsPt1u2lygu1u89o42qOmQIE0q7WbxYjsjllBo5sj0mI+RMbNXvfNpoo9HG0AgdVZUnOftJVS8jykpswErOxFb9zhfvJHZElHH4LSdnYqt+54t3EjsiyjhMRsiZ4p3EjuwlnknsiCjjMBkh53LaOBs0kNPG2SCilGDXXnI+J42zQdE5ZZwNIjKMXXspe7hcMgHhiYyIyJGYjJCzhQ6gpSUjetPPk30xhkRZj990cq5oU8p3dcnBtNjmwP4YQyICkxFysq6uYGNHjcsl/+7qkg8OJW5vHR0yTkOHBm+vMYbOwfZalCRMRsiZVFVW6+sNIw5wKHG76+uTiUhrq6wN6e2VMQu9NcMY2pve/DTl5bLLPXuykUlMRsiZAoFg+wI9HErcvrRbM11dMk6FhTJWX3wRfmuGMbSvaDNmNzcHZ/RlQkImcJwRcqZUDCWuqvJkqKrJ2UbSp91eGzpUzksjhEw2CgqCA54BjKFdqaqsEdEbFUJbdvfdjAGZwmSEnCmZQ4n39QHt7UBbW/DR3i6XU3KF3l5zueTP0Bjm5cm/VZUxtCszM2YTGcTbNORchYXyKrirK7wnRk+P8aHE2ZsjvSJvr0XGUFFkLDo75XOMof1wxmxKAdaMkHMlYyjx0B45LldwWnvtNbVbBpQckbfXQmOoqrLdiKrKWziMoT1xxmxKAV4ukLPl5squn0OGmB+BlT1y0k+7vdbVFaz10GKoqrJGZOhQYNgwY6/HGKafNmN2c7N+uxFFCXbzJTKINSOUGVwu2d7AzAnHTI8cSp5oM/X29MjnLrrI+GsxhunHGbMpBZiMUPZKRY8cGlwyZ+plDK3BGbMpyXibhpwtkREg9W4ZhNKu1HmFl3yRt9eEkElD5JX2YBhD6yxeDNTWcgRWSgomI+RcyRgBMhk9cih+2mBniUySxxhax+UC5syxeisoAzAZiUJVgQMHZButsjKgupoJv60kawRI7ZaBNmusdoWudStll9DUSVaXXMaQKCqnnMsUIfSaQ9uL3+9HcXExfD4fPB5Pyt+voQFYvRo4cya4rLISWL9e1kqSxVQVuOyy6OMYaK35T50y961TVfM9cih+7e3Rb69oy81OkscYEvWzw7nM6Pk7oVZdjzzyCBRFQV1dXdQyGzduhKIoYY+CgoJE3jalGhqARYvCgwfIAQUXLZLPk8VSNQJkPD1yKD5Gu+SaHVKcMSQC4LxzWdzJyDvvvIOnn34aU6dOHbSsx+NBc3Nz/+PUqVPxvm1KqarMImNNuVBXxykXLMcRIJ2PXXKJUsaJ57K4kpGOjg4sW7YMzzzzDIYPHz5oeUVRUFpa2v8oKSmJ521T7sCBgVlkKCGApiZZjizEESCdj11yiVLGieeyuL7pK1euxPz58zF37lxD5Ts6OjBq1Ch4vV4sWLAA7733Xszy3d3d8Pv9YY90aG5ObjlKEW0EyGjdQBUFqKjgCJB2lsyJDokojBPPZaaTkS1btuDIkSOor683VH7cuHF47rnnsGPHDmzatAmBQACzZs3CmRhpW319PYqLi/sfXq/X7GbGJXL8nkTLUYpwBMjMEG0k1q4udsklSoATz2WmetM0NTVhxowZ2LNnT39bkZqaGkybNg3rtJPDIHp7ezFhwgQsXboUDz74oG6Z7u5udHd39//t9/vh9XpT3ptGVYHRo2UDn2hTLlRWAo2NPM/Zgt44IxUVMhHhCJD2pw0Br43AqnXJNTvOCFmHvZdsyU7nMqO9aUwlIy+//DK+8Y1vwBWy9aqqQlEU5OTkoLu7O+y5aBYvXozc3Fz86le/MvS+6ezaq7VABsKDqF1wb9vG7r22ksgIrGQNbSbd0IHOcnNlEpKfz/g5gV4MmUTail3OZSnp2jtnzhwcP34cx44d63/MmDEDy5Ytw7FjxwwlIqqq4vjx4yizU/1QiNpaGaSKivDllZVMRGxJGwHy1lvlT57I7E0b6Ey7FVNQIH/29Mhl9h/2iKLFsKsrOIgdWc5p5zJTKWxRUREmT54ctmzo0KG4+OKL+5cvX74cFRUV/W1KHnjgAcycORNjxozB+fPn8dhjj+HUqVO4/fbbk7QLyVdbCyxY4IxR64gcpasrOCmexuWSf3d1yYfZgc4ovRhDx3DSuSzp9WmnT59GTkh3vLa2Ntxxxx1oaWnB8OHDcfXVV+PgwYOYOHFist86qVwudsYgSiqjA50NGWLPoyUxhg7klHMZh4MnovTo7QXa2mS1vl63bCHkiWzYMDmKKtkPY0gmpWU4eCIiwzjQmfMxhpQi/I8hovTgQGfOxxhSijAZIaL04UBnzscYUgqwQzgRpU9uLuDxBMeo0AY6KyzkGBVOwRhSCvC/hojSKzdXdv0cMoSjdzoVY0hJxmSEiKzhcvEE5nSMISUJ24wQERGRpZiMEBERkaV4m4aI0oszvTobJ6ekFGAyQkTpwZlenW/rVqCuTiYimvJyYN06YPFiq7aKMgCPAESUetpMr319cv4SbRTPri45xLjHw4TE7rZuBZYsGTizcnOzXA4wIaG4sc0IEaVe6EyvLpccl0Kb6VWrMSH7UlVZI6I3lZm27O67ZTmiODAZocygqvIKmwdD+zE602tPD2NoV1obkWiEAM6eleWI4sB6UXI2tkOwv0AgGJtoz/v9MhFxuRhDO4qViMRTjigCa0bIubR2CNqcGAUFwTkztPYJZL1YM7329QE+n6wVyc9nDO2qvDy55YgiMBkh52I7BGeINdNrV5d8eDwyGWEM7ammRiYaiqL/vKIAFRWyHFEcmIyQMxlth8D2B/agN9NrT4+sFSko0J/plTG0D5dLdt8FBiYk2t9PPMHxRihuTEbImQZrh5CTI094ercGKP20mV61Gg+twarbDQwbpt82hDG0l8WLgRdfBMrKwpeXl8vl7NZLCWDrMHKm0HYIeldjgYC8YouWrFD6Rc70KoRsFxKrYStjaC+LFwO1tRyBlZKOyQg5k9YOoatLv4q/pyfYloTsJXSmV8bQeVwuYM4cq7eCMgwvOci59NohqGqwd43eCY7shTEkIrBmhJxMa4egjTMihKzWLyzkGBVOwRgSEZiMkNNFtkPgTLDOwxgSZT0mI5QZQtshkDMxhkRZi8kIEaWHqrLmw8kYP0ohJiNElFqcP8jZGD9KA/4nEVHqaPMH9fXJEVW1sWG6uuTEeB4PT2h2xvhRmrBrL2UWVeU09HYSz/xBjKF9cP4nShOmtJQZWJWcXkbaDxidP2jIEPkajGF6DRZDs/EjSgC/4eR8kVXJQsjfOzrkFfbQocHh4zPhoJloQ8JE1jeTMESbP0h7f0UJzj2jDQ2vF8MvvgAuukguz4T4Ac6IoZn4ZUpcyDJZm4yoKnDgANDcLOd9qq7m98mxtKrkvDzgwgU5jLh2ID13TiYjHk9iV9p26EmQaM1BMtY3034gcv4g7f21+Aghl3s8cllkDLVHVxdQXAyMGJFYTQljaC6GZuJHjtXTAzz5JPDxx8CXvwz84z9GrwxLpYTajDzyyCNQFAV1dXUxy23duhXjx49HQUEBpkyZgldffTWRt01YQwMwejRw443AzTfLn6NHy+XkMFpVck6OPMh+8UVwvIoLF+Tj/Hl5FacNO64djI3o6wPa24G2tuCjvd34+sminUS0YdILCsztT6LrA+bbD2jzB2mJRmh83O5gwuj3A52d4TEUQq7X2yvfp6MjOEy8mfhp+84YSmZiaDR+nZ3p/ywpKb7/fXmX7e67gf/4D/lzyBC5PN3iTkbeeecdPP3005g6dWrMcgcPHsTSpUuxYsUKHD16FAsXLsTChQtx4sSJeN86IQ0NwKJFwJkz4cvPnpXLmZA4jHZA7O6WJ6uCAnmw7O6Wz48YIU9sFy6Yb3gXefDPy5Ov1dFh/oSYqEQbEia6vtH2A1qjU60Ran6+/Oz+8hd5UnO75fNffCGfGz5cLr9wITyGvb0yrkOHytFZVVWuY7bhJGMYFE8Mc3Pl+8SKHxuyOtL3vw889tjAduKqKpenOyFRhBDC7EodHR246qqr8OSTT+LHP/4xpk2bhnXr1umWXbJkCTo7O7Fz587+ZTNnzsS0adPw1FNPGXo/v9+P4uJi+Hw+eBKoElRVWQMSmYhoFAWorAQaG3nLxjFUFWhtBfbulQfMsjLgb/5GXvlq1c09PbLNwcUXy8CqqjyADh8eO9Dt7fIgm5cXXj2dkyNfY8QI+RpGtjHR9gFtbXJ/ojU0jLU/ia4PyOSgrU0mCooy8HnthOTxyJN96G0EIFg7lZMjf+bnB28t9PTIGCpK8PW1+GnxCo2h9n6DxQ9gDEPFiqGqBhPI4mL5WloMtfXy8mTcIuNn9PtEttHTI2tAYnVY0yqXE71lY/T8HVfNyMqVKzF//nzMnTt30LKHDh0aUG7evHk4dOhQ1HW6u7vh9/vDHslw4ED0RASQx9CmJlmOHKKhAZgxA1i+HKirA5YsAWbOBHbvlieOnh558NXufQPyd63hXTTRbv+43fJnIAD83/8Fa2D0JOv2QLSGhJrB9ifR9bUyoZ+hRtvHv/wF+Pxz+QX77DP5nHYbQbva9njkiW7YMFnboSWLiiIThdZW+Vm3tcmaC+29ImNoZHsBxlCvTGQMtf07fz4Yw+bmYDJYUBBMQHJzZU1VZPyMxoNs48knB+85r6qyXLqYTka2bNmCI0eOoL6+3lD5lpYWlJSUhC0rKSlBS0tL1HXq6+tRXFzc//B6vWY3U1dzc3LLkcW2bpXJR2TAWluB730P2LFDHjDd7uBVOSCv9Pr65AE0mmi3f7Sq8SFD5Ot0duqvn4z7+5poiUDotobuX7LXB8LbD2hC2xEEAsHeLkLIz0VV5d9Dh8rnv/hCXl2HXj339QE+n/xstEu1nBz5uWtJSWQMtdhoV/PRMIbhImMY2Q5EW/+LL+Tnos2gnJcXjGFf38DaDyPvTbby8cfJLZcMpv57mpqasHr1avzyl79EQUFBqrYJa9euhc/n6380NTUl5XXLypJbjiykqrImJFZC8fjjwZlgte6i7e3yyl07oUS7ytWu9rQTaKRAQB7Yow3OlczBovQSgVDavfxoVeSJrq/RquW7uuQ+a71ehAjePigoGLiPLpe8xaIlKKG6uuTj0kuB0lJZVushcuGCLFNUFB7Dzz+Xz2m1KIxhfDHs6Aj2YOrpkevm5cnPOzJ+Wl19T8/Az8roe5NtfPnLyS2XDKaSkXfffRfnzp3DVVddhdzcXOTm5mL//v346U9/itzcXKg6X+jS0lK0traGLWttbUVpaWnU93G73fB4PGGPZKiulm1C9G55A3K51yvLkc3t2wd8+mnsMufOAW+8IX/PzZUnrvPn5YFz2LDYV7nagVmr5o+kd/tHY7ahoBGRiYAQwd4lubnyNWONWjrY+oWFg29Dbq681VJYKLe/vV3uf2Fh+Fgu2j6Gnri057WEROst4/MFE5iiIhkXt1u2DbnkErlPFy4EX/vcORmroqLBayoYw4G0GObny89Ve43IGEbGr7Aw2O5GqzUx+95kG//4j4Pnji6XLJcuppKROXPm4Pjx4zh27Fj/Y8aMGVi2bBmOHTsGl87eVVVVYe/evWHL9uzZg6qqqsS2PA4uF7B+vfw9MiHR/l63jgm+IwyWiGg6OmTtSGenPHkMGyYf2u2CWFe5Q4cGx73QO/hH3v7RJOP+fqTQREBrXKhd1WoNPmO1aYi2fmGhuflFcnNlIlBcLH+OGCF/ho5ZobePOTnyfYYODb6/dkWtJYbaNmoDnF10UbCbttsdrOEYOTL42TOG8cVQG3tn+PBgG5DQGEbuX26ujEdBQfDWVzzvTbaQnw/cc0/sMvfck97xRkz9BxUVFWHy5Mlhy4YOHYqLL764f/ny5ctRUVHR36Zk9erVuOGGG/D4449j/vz52LJlCw4fPoyf//znSdoFc2prgW3bgNWrwxuzVlbKRKS21pLNIrPKy42VmzRJnjA//zx4RRgp2rDWbre8Om9rC/YYUBR5AC4slFeIWhV+qMjBoiLFe49dSwS0W0+hYzwYGYQscv1EBv/KywueyIFgVb7W/iByH3t6wrvpau0PgPBbbdo2hvbsGD5clvf5ot8OYAzNyckJJuTa+qEx1BKk0P0LBORnqb23ELK1/6efyu9jTQ2v5Bzk0Uflz5/8JLwyzuWSiYj2fLokPZ09ffo0ckL+gWfNmoXNmzfj3nvvxb/8y79g7NixePnllwckNelUWwssWMARWB2tpkYeAJub9duNKErwAKmdUPTaDQDhV4GR/wRFRcGuqtrVu6LIk2u06mntSr6rS//5nh79E6BR2glEu3oOfQ+tpkBri1FUFH39ROjto3Zy/+ILeXQbMkQuj6zK19pxdHfL53y+4BW/duJ1uYL9D7XbF4M18GQMzb2/3v5pMWxvl++rdYEO/axcLtmLra4uvIayvFxe0S1eHP92UVo9+ijw4x/bYwTWuMYZSbdkjTNCGUbrTQOEJyTa1fqLL8oDY6JjNIQOwa31MBhsCO5ow25rB/VEq7aTMe5EovT2UWsHEgjIan2tBkX7rCLX0Wo8urrkLYBhw4KvE/o5MYbxv0800fbvwgWZkGi1UJGflfa9izx1RH7viGD8/M1khJxt69aBV2gVFcATT4QfELXBr/SucrXlelegocwOfBXPCdCowQYh02oCtPYxqRJtH7WTW+RnpRcH7TW0xqxFRfqfE2OYfLH2T1EGflaqClx2WfQ2W1qN5KlTrGomAMbP32x1RM62eLG876b1rol271qrfu7q0r/KNdIbwGzVeDLbaERKVZsGs8zsY7QeKtpraF1Xo7XtYQyTb7D9i9yuwXqxCSHn1ti3D5gzJxVbTBmKyQg5n8s1+IFPq1aPvArUGjKmsjdAMtpo6L1mKts0xLM9g73XYD1U8vKCJ2A9jGHqGN0/o73YjJYj+ismI5QZjFS/p/Iq1wrJqClIp8FqAoyMjMsYWstoLzaj5Yj+iskIOVvoPW/t5DTYPf1UXOVawcqagnhEqwkIbTPidstGlYyhPWNophcbkQk2+08nMiFab4BoYzRkIqfVFETWBETrTcMYWr1V+lwu2X13yRKZeOj1YnviCftuP9kWZzYi50rm3CFOp42jYveTQOQooufPB3uMDB9ubGTcTOWUGC5eLLvvRk7iVV7Obr0Utwy/5KCMZXTukMgROcl6ob1ntFFWzYyMS9Yz2ouNyCAmI+RMZuYO4QHSnrSarHhGxiXrGenFRmQQb9OQM4X2zNCTrnE2KH6MIRH9Fb/l5Exaz4yeHv3ntRlheUVtX4whEf0VkxFyLq3rY1eX/vTwdhujgQZiDIkIbDNCTua0MRpoIMaQiMBkhJzOSWM0kD7GkCjr8TYNERERWYo1I+Rs8QwHT/bCGBJlPX7Tybk4HLzzMYZEBN6mISfjcPDO19ERnKeGMXQmVZWJo6pavSXkYLzkIGficPDO1tcnE5HWVlkb0tsrYxZ6a4YxtDfeXqMkYs0IOZOZ4eDJXrRbM52dMk5azdYXXwRv2QCMoZ1pMdTGgykoCI4XExpDIoOYjJAzpWIocVY3p4d2e23oUDkvjRAyGSkoCA54BjCGdsZbpJRkrEsjZ9KGEu/q0h+ls6cneKAcDKub0yf09prLJX+GxjAvT8ZOVRlDu+ItUkoB1oyQcyVjKHFWN6dX5O21yBgqivzMOzsZQ7viLVJKASYj5FzaUOJa1XB3d7Dq2GiXUFY3p1fk7bXQGKqqbDeiqvIWDmNoT5xtmVKA9ZfkbIkMJc7q5vTTu72mxVBVZY3I0KHAsGHGXo8xTL9k3iIl+iumrpQZXC7Z3sDMAZDVzdaIdntNO4lddJHx12IMrcHZlinJWDNC2Su0ulkviWF1c2okc6ZextAanG2Zkoz/MeRsqgrs2wd8+ilQXg7U1BivHWF1s3Uib69pJzNFMfc6jKF1ONsyJRGTEXKurVuBujqZiGjKy4F164DFi429RmGhHJdCG5Jcu8ru6WF1czoIIRutJtIllzG0lsvFJIQSxmSEnGnrVmDJEnkyC9XcLJcDxhISVjdbJ1mT5DGGRI6nCBF5NLcfv9+P4uJi+Hw+eDweqzeHrKaqwGWXhdeIhFIUWUNy6pS5KzZVZXVzOrW3R7+9oi0vKjL3mowhka0YPX+batW1YcMGTJ06FR6PBx6PB1VVVXjttdeilt+4cSMURQl7FBQUmHlLooG0NiLRCAGcPSvLmRFPjxyKj9EuuWaHdWcMiRzJVP1lZWUlHnnkEYwdOxZCCLzwwgtYsGABjh49ikmTJumu4/F48MEHH/T/rZhtoGYRVQUOHJC1/mVlQHU1j2+2ESsRiaccpZ+ZLrn84hHFzSnnMlPJyN/93d+F/f3QQw9hw4YNePvtt6MmI4qioLS0NP4ttEBDA7B6NXDmTHBZZSWwfj1QW2vddtFflZcntxylH7vkEqWck85lcX/TVVXFli1b0NnZiaqqqqjlOjo6MGrUKHi9XixYsADvvffeoK/d3d0Nv98f9kiXhgZg0aLw4AGy1n/RIvk8WaymRiYa0WrZFAWoqJDlyJ60Lrk9PfrP9/TI5+14CUfkAE47l5lORo4fP46LLroIbrcbd955J7Zv346JEyfqlh03bhyee+457NixA5s2bUIgEMCsWbNwJvLTiVBfX4/i4uL+h9frNbuZcVFVmUXqNenVltXVcXZyy7lcsvsuMDAh0f5+4gmeyOyOo3gSpYQTz2Wme9P09PTg9OnT8Pl82LZtG5599lns378/akISqre3FxMmTMDSpUvx4IMPRi3X3d2N7u7u/r/9fj+8Xm/Ke9Ps2wfceOPg5d54gxfdtqA3zkhFhUxEjI4zQtbRhoDXJjjUuuSaHWeErJHIgIOUUnY6lxntTWP6256fn48xY8YAAK6++mq88847WL9+PZ5++ulB183Ly8P06dPx0UcfxSzndrvhdrvNblrCmpuTW45SbPFieeOTB0Rn0WbSDR3oLDdXJiH5+YyfEyRjwEFKGSeeyxK+9AgEAmG1GLGoqorjx4/jpptuSvRtU6KsLLnlKA1cLmDOHKu3goyKNtBZT4/8mZdn9RbSYJI14CCljBPPZabajKxduxa///3vcfLkSRw/fhxr167Fvn37sGzZMgDA8uXLsXbt2v7yDzzwAF5//XV88sknOHLkCG655RacOnUKt99+e3L3Ikmqq2VL41jtIr1eWY6I4tDVJRMRbb4YRZE/CwuDNSZkX6oqa0RiNUa4+257NUbIQk48l5lKRs6dO4fly5dj3LhxmDNnDt555x3s3r0bX/nKVwAAp0+fRnNIvU9bWxvuuOMOTJgwATfddBP8fj8OHjxoqH2JFVwu2eUJiN4uct061iITxSVVA51R+qRqwEFKKieeyzgcvA69vtlerwye3fpmEzlGby/Q1gYUFOhfsgkhk5Fhw3i7xq5+8Qtg+fLBy/3nfwK33pr67aGY7HAuS1kD1mxQWwssWOCMUeuIHIMDnTkfBxx0FCedy1gzQkTpk4rJ8Sh9tEkqm5v1243EO0klZayUTJRHRJQQDnTmbBxwkFKEyQgRpU9uLuDxBHvPaAOeFRbK5RzozP4WLwZefHFgv9Dycrmc3XopDrxNQ0TWUNXgoGe8knYejsBKBrABKxHZm8vFk5eTccBBSiLepiEiIiJLMRkhIiIiS/E2DRGlF9uKOBvbilAKMBkhovTQm63X7Q529yX742y9lCI8AhBR6kWbrberSw4Tz2699sfZeimF2GaEiFKPs/U6G2frpRRjMkLOp6rA3r1yEq+9e3lAtBujs/X29MhaEsbPfjhbL6UY60XJ2XgP2/4CgWAbkWjP+/0yEXG52JbEjmIlIvGUI4rAbzo5F+9hO0Os2Xr7+gCfT9aKFBcDeXlsS2JHnK2XUoy3aciZeA/bOVwuWdPR0zPwua4u+fB45O0atiWxp5oamWhETo6nURSgokKWI4oDkxFyJt7Ddha92Xp7emStSEGB/my9WlsSJpTW42y9lGJMRsiZeA/bWfRm6+3pkTUmw4bp34rJyZFJSyCQ9s0lHZytl1KIN2PJmXgP23lyc4GiImDIEJlgCCEbrsZq2Koo0Z+n9Fu8GKit5QislHRMRsiZtHvYzc367UYUJXigJHsJna3X7Za3bvRu0/T0BMclIfvgbL2UArzkIGfiPezMoNeWRFXl37m5+kkKEWUcJiPkXLyH7Xx6bUm0kVrZrZcoayhC6NVx24vf70dxcTF8Ph88Ho/Vm0N2w1lEMwNn8yXKOEbP37zsIOfjPezMENqWhIiyCpMRIkoP1nw4G+NHKcRkhIhSSxtJtbs7eDLj3DPOwfhRGvA/iYhSp69PjiXS1ydHVNXmqOHcM87A+FGasDcNZRZV5TT0dtLVFewd43IZm3uGMbSPeOJHFAemtJQZWJWcXkbaD6iqjEd+vv7z2twzQ4bI12AM02uwGJqNH1EC+A0n5wutSna5gP/+b9nN95JLgOuvl1XJ2jT2Tj9oJqMbcyINEc0kDIFAsIze+ytKcO4ZbWh47XaAEPJ3nw/Yu1f+rKzMjG7bicYw0fWNxlAvfqoKvPmmHPm4tBS45hpZxukxIcsxGUkzVQUOHJDf5bIyoLqa3+OEaVXJu3YBa9bID1dTUgL867/KAdDivdK2yzgmW7cCdXXhk/+Vl8uRaI0M8JZozYPZ9gNaAqidrLT37+kJJiAul1yvp0c+n5cHXLgg/965E/jxj4HW1vj2N5IdeoMkGsNk/A8YjWFk/LZvH/j9Ki2V771kyeDvTTFl/blBmPDkk0+KKVOmiKKiIlFUVCRmzpwpXn311ZjrvPTSS2LcuHHC7XaLyZMni1deecXMWwohhPD5fAKA8Pl8pte1k//6LyEqK4WQR2H5qKyUyylOfX1CfPaZEM8/L4SihH+4oY/nnhOirU2I1lYhPv9ciN5eY6//0ktClJeHv1Z5uVyeTi+9pL9/iiIfg21Pb6/c79ZW+Tn4fOY/D79flvf7Bz605dHW+ctfhDh5UoiPPxbi9GkhzpwR4s9/ln+fOydEc7OMo1bmyScT29/Ifff75eu3tsqffr/x/4FkSTSGia4vhPkYastfeCF58aABMvncYPT8bWoE1t/85jdwuVwYO3YshBB44YUX8Nhjj+Ho0aOYNGnSgPIHDx7E9ddfj/r6enz961/H5s2b8W//9m84cuQIJk+ebDhhyoQRWBsagEWLBs7ppk2jsm2bnAyTTOrtBf7v/4Crrw6/YotUVga8/7681NAmZisqiv3aW7fKK75oQUvXkPOqClx2WfjVcOT2lJcDp05Fv5Rqb48+IZ2Rz0NVgbY2edUcrX1BXx8wfHj489qVeFubvMLWZuzt6QkOBX/hgqytcbvlz7w84G/+BmhpiX9/I98/siYg9P3T0R4l0Rgm438gnhj29cl1pk5NTjxogEw/Nxg9fyc8HPyIESPw2GOPYcWKFQOeW7JkCTo7O7Fz587+ZTNnzsS0adPw1FNPGX4PpycjqgqMHg2cOaP/vKLI2+GNjfwum6aqwI4dwP/7f4OX/c1vgBtuiH7ijHzdRA/+ybJ3LzB37uDlfvtb/ZFo400kQvX2ytcoKBg4MSEgj6Td3cCwYTKZCNXdLT9H7RaJosjEQLs91NMjb8VoE+MdOmSs2j/a/oZKNAlLlkRjmOj6QPwxfP11YN68xN6bdGXDucHo+Tvurr2qqmLLli3o7OxEVVWVbplDhw5hbsQXaN68eTh06FDM1+7u7obf7w97ONmBA9H/2QB5DGhqkuXIJJcL+PxzY2W1mpOcnGDDyWi0NiLRCAGcPSvLpVqs7TBSLlpDUo2RzyO0/UC091AU/ffIyQGGDgVGjACKi+XJrqgoWCORlxecuTcnJ/oVeKTBPhejvUHS0YU40Rgmuj4QfwxD2+zE+96ki+eGINPJyPHjx3HRRRfB7XbjzjvvxPbt2zFx4kTdsi0tLSgpKQlbVlJSgpZBDjb19fUoLi7uf3i9XrObaSux7h7EU44ijBplrJw2u2+sE6cmGQf/ZCkvT6xcIomExuWSt1F6evSf7+mRz+tdvmnvrygy8YgsEwjIZMXtlrdsRo6Mvh2hBvtckpGEJUuiMUx0fSD+GCbjvUkXzw1BppORcePG4dixY/jv//5v3HXXXbjtttvw/vvvJ3Wj1q5dC5/P1/9oampK6uunW+QM94mWowhz5sgDoV7VMxC8pXLddfLvWCdOjZ0OwDU1g+9fRYUspyeRRCKUdlulq0vWJgghf3Z1BW+xxPv+Ho/s+eRyye6iERcxYQbbX00ykrBkSTSGia6viSeGyXpvGoDnhiDT38L8/HyMGTMGV199Nerr63HllVdi/fr1umVLS0vRGlHF19raitLS0pjv4Xa74fF4wh5OVl0t7/vF+i57vbIcxcHlkt0LgYEfsvb3ww/Ln4OdODV2OgAb2b8nnoidTMSbSITSGnxqo292dwdH5xysIaiR9y8qkrdwhgwBHnpI7lu8+wskLwlLhkRjmIz/ASC+GCbrvWkAnhuCEr4kCAQC6O7u1n2uqqoKe/fuDVu2Z8+eqG1MMpXLBWj5WrTv8rp1/C4nZPFi2bsl8hKirAx45hngppuMnzgB+x2Ao+1febmxXj2JJBKRr1NUJBu7Dhsmf4a2/0jk/bUyF10ELFgg4xZZQ2J0fzXJSMKSJdEYJrq+Jp4YJuu9KQzPDUGmetOsXbsWX/va13DZZZehvb29v6vu7t278ZWvfAXLly9HRUUF6uvrAciuvTfccAMeeeQRzJ8/H1u2bMHDDz+clV17AdmFa/Xq8AZLXq/8Z3Ny1y1b0RugDIh/sCu9QaYqKmQiYsUB2OoRWJPB6FDy2sBoBw4ktr+hg70JIY/yVg4zb/UIrImwywCAGSaTzw0p6dq7YsUK7N27F83NzSguLsbUqVPxgx/8AF/5ylcAADU1NRg9ejQ2btzYv87WrVtx77334uTJkxg7diweffRR3HTTTSnZGSfI+lH2nIgH4MxgdRJGFEOmnhvSNs5IOmRSMkJERJQtUj7OCBEREVEyMBkhIiIiSzEZISIiIksxGSEiIiJLMRkhIiIiSzEZISIiIksxGSEiIiJLMRkhIiIiSzEZISIiIksxGSEiIiJLMRkhIiIiSzEZISIiIksxGSEiIiJLMRkhIiIiSzEZISIiIksxGSEiIiJLMRkhIiIiSzEZISIiIksxGSEiIiJLMRkhIiIiSzEZISIiIksxGSEiIiJLMRkhIiIiS+VavQFkLVUFDhwAmpuBsjKguhpwuazeKiKi1OKxz16YjGSxhgZg9WrgzJngsspKYP16oLbWuu0iIkolHvvsh7dpslRDA7BoUfiXEQDOnpXLGxqs2S4iolTisc+eFCGEsHojBuP3+1FcXAyfzwePx2P15jieqgKjRw/8MmoURV4lNDay2pKIMgePfeln9PzNmpEsdOBA9C8jAAgBNDXJckREmYLHPvtiMpKFmpuTW46IyAl47LMvJiNZqKwsueWIiJyAxz77YjKShaqr5X1RRdF/XlEAr1eWIyLKFDz22ZepZKS+vh7XXHMNioqKMHLkSCxcuBAffPBBzHU2btwIRVHCHgUFBQltNCXG5ZJd2ICBX0rt73Xr2ICLiDILj332ZSoZ2b9/P1auXIm3334be/bsQW9vL7761a+is7Mz5noejwfNzc39j1OnTiW00ZS42lpg2zagoiJ8eWWlXM6+9kSUiXjss6eEuvZ+9tlnGDlyJPbv34/rr79et8zGjRtRV1eH8+fPx/s27NqbQhyFkIiyEY996WH0/J3QCKw+nw8AMGLEiJjlOjo6MGrUKAQCAVx11VV4+OGHMWnSpKjlu7u70d3d3f+33+9PZDMpBpcLqKmxeiuIiNKLxz57ibsBayAQQF1dHWbPno3JkydHLTdu3Dg899xz2LFjBzZt2oRAIIBZs2bhTIzO3vX19SguLu5/eL3eeDeTiIiIbC7u2zR33XUXXnvtNbz55puorKw0vF5vby8mTJiApUuX4sEHH9Qto1cz4vV6eZuGiIjIQVJ6m2bVqlXYuXMnfv/735tKRAAgLy8P06dPx0cffRS1jNvthtvtjmfTiIiIyGFM3aYRQmDVqlXYvn07fve73+Hyyy83/YaqquL48eMo46gyREREBJM1IytXrsTmzZuxY8cOFBUVoaWlBQBQXFyMwsJCAMDy5ctRUVGB+vp6AMADDzyAmTNnYsyYMTh//jwee+wxnDp1CrfffnuSd4WIiIicyFQysmHDBgBATUQT5Oeffx7f/va3AQCnT59GTk6wwqWtrQ133HEHWlpaMHz4cFx99dU4ePAgJk6cmNiWExERUUZIaJyRdOE4I0RERM5j9PzNuWmIiIjIUkxGiIiIyFJMRoiIiMhSTEaIiIjIUkxGiIiIyFJMRoiIiMhSTEaIiIjIUkxGiIiIyFJxTZRHZAeqChw4ADQ3A2VlQHU14HJZvVVE2YPfQUoWJiPkSA0NwOrVwJkzwWWVlcD69UBtrXXbRZQt+B2kZOJtGnKchgZg0aLwgyAAnD0rlzc0WLNdRNmC30FKNs5NQ46iqsDo0QMPghpFkVdnjY2sLiZKBX4HyQzOTUMZ6cCB6AdBABACaGqS5Ygo+fgdpFRgMkKO0tyc3HJEZA6/g5QKTEbIUcrKkluOiMzhd5BSgckIOUp1tbwfrSj6zysK4PXKckSUfPwOUiowGSFHcblk10Fg4MFQ+3vdOjacI0oVfgcpFZiMkOPU1gLbtgEVFeHLKyvlco5xQJRa/A5SsrFrLzkWR38ksha/gzQYo+dvjsBKjuVyATU1Vm8FUfbid5CShbdpiIiIyFJMRoiIiMhSTEaIiIjIUkxGiIiIyFJMRoiIiMhS7E1DZDPsLpl5GFOi2JiMENlIQwOwenX4rKiVlXLESw4k5UyMKdHgeJuGyCYaGoBFiwZOz372rFze0GDNdlH8GFMiYzgCK5ENqCowevTAk5ZGUeTVdGMjq/edgjElMn7+Zs0IkQ0cOBD9pAUAQgBNTbIcOQNjSmQc24wQ2UBzc3LLOUUmN+zM1pgSxcNUzUh9fT2uueYaFBUVYeTIkVi4cCE++OCDQdfbunUrxo8fj4KCAkyZMgWvvvpq3BtMlInKypJbzgkaGuRtjBtvBG6+Wf4cPTpz2lFkY0yJ4mUqGdm/fz9WrlyJt99+G3v27EFvby+++tWvorOzM+o6Bw8exNKlS7FixQocPXoUCxcuxMKFC3HixImEN54oU1RXy/YDiqL/vKIAXq8slwmyoWFntsWUKBEJNWD97LPPMHLkSOzfvx/XX3+9bpklS5ags7MTO3fu7F82c+ZMTJs2DU899ZSh92EDVsoG2gkakO0JNNrJbNu2zOgKmk0NO7MlpkTRpKUBq8/nAwCMGDEiaplDhw5h7ty5YcvmzZuHQ4cORV2nu7sbfr8/7EGU6Wpr5cmpoiJ8eWVlZp20sqlhZ7bElChRcTdgDQQCqKurw+zZszF58uSo5VpaWlBSUhK2rKSkBC0tLVHXqa+vx49+9KN4N43IsWprgQULMrdRJ5B9DTuzIaZEiYo7GVm5ciVOnDiBN998M5nbAwBYu3Yt7rnnnv6//X4/vF5v0t+HyI5cLqCmxuqtSJ1sbNiZ6TElSlRcyciqVauwc+dO/P73v0dlZWXMsqWlpWhtbQ1b1traitLS0qjruN1uuN3ueDaNiGxOa9h59mx4OwqN1maEDTuJsoepNiNCCKxatQrbt2/H7373O1x++eWDrlNVVYW9e/eGLduzZw+qqqrMbSkRZQSXS87LAgzsaaL9vW4db2MQZRNTycjKlSuxadMmbN68GUVFRWhpaUFLSwu6urr6yyxfvhxr167t/3v16tXYtWsXHn/8cfzpT3/C/fffj8OHD2PVqlXJ2wsichQ27CSiUKa69ipROsw///zz+Pa3vw0AqKmpwejRo7Fx48b+57du3Yp7770XJ0+exNixY/Hoo4/ipptuMryR7NpLlJkyeQRWIjJ+/uZEeURERJQSnCiPiIiIHIHJCBEREVmKyQgRERFZiskIERERWYrJCBEREVmKyQgRERFZiskIERERWYrJCBEREVmKyQgRERFZKq5Ze9NNGyTW7/dbvCVERERklHbeHmywd0ckI+3t7QAAr9dr8ZYQERGRWe3t7SguLo76vCPmpgkEAvj0009RVFQUdbK+ePj9fni9XjQ1NWXsnDeZvo/cP+fL9H3k/jlfpu9jKvdPCIH29naUl5cjJyd6yxBH1Izk5OSgsrIyZa/v8Xgy8h8sVKbvI/fP+TJ9H7l/zpfp+5iq/YtVI6JhA1YiIiKyFJMRIiIislRWJyNutxv33Xcf3G631ZuSMpm+j9w/58v0feT+OV+m76Md9s8RDViJiIgoc2V1zQgRERFZj8kIERERWYrJCBEREVmKyQgRERFZKmOTkfr6elxzzTUoKirCyJEjsXDhQnzwwQeDrrd161aMHz8eBQUFmDJlCl599dU0bG184tnHjRs3QlGUsEdBQUGatticDRs2YOrUqf0D8VRVVeG1116LuY6T4md2/5wUOz2PPPIIFEVBXV1dzHJOimEkI/votDjef//9A7Z3/PjxMddxUgzN7p/T4gcAZ8+exS233IKLL74YhYWFmDJlCg4fPhxznX379uGqq66C2+3GmDFjsHHjxpRuY8YmI/v378fKlSvx9ttvY8+ePejt7cVXv/pVdHZ2Rl3n4MGDWLp0KVasWIGjR49i4cKFWLhwIU6cOJHGLTcunn0E5Ch7zc3N/Y9Tp06laYvNqaysxCOPPIJ3330Xhw8fxt/+7d9iwYIFeO+993TLOy1+ZvcPcE7sIr3zzjt4+umnMXXq1JjlnBbDUEb3EXBeHCdNmhS2vW+++WbUsk6MoZn9A5wVv7a2NsyePRt5eXl47bXX8P777+Pxxx/H8OHDo67T2NiI+fPn48Ybb8SxY8dQV1eH22+/Hbt3707dhoosce7cOQFA7N+/P2qZb37zm2L+/Plhy6699lrxD//wD6nevKQwso/PP/+8KC4uTt9GJdnw4cPFs88+q/uc0+MnROz9c2rs2tvbxdixY8WePXvEDTfcIFavXh21rFNjaGYfnRbH++67T1x55ZWGyzsthmb3z2nx+8EPfiCuu+46U+t8//vfF5MmTQpbtmTJEjFv3rxkblqYjK0ZieTz+QAAI0aMiFrm0KFDmDt3btiyefPm4dChQyndtmQxso8A0NHRgVGjRsHr9Q56JW4Xqqpiy5Yt6OzsRFVVlW4ZJ8fPyP4BzozdypUrMX/+/AGx0ePUGJrZR8B5cfzwww9RXl6OL33pS1i2bBlOnz4dtawTY2hm/wBnxe/Xv/41ZsyYgcWLF2PkyJGYPn06nnnmmZjrWBHDrEhGAoEA6urqMHv2bEyePDlquZaWFpSUlIQtKykpQUtLS6o3MWFG93HcuHF47rnnsGPHDmzatAmBQACzZs3CmTNn0ri1xh0/fhwXXXQR3G437rzzTmzfvh0TJ07ULevE+JnZP6fFDgC2bNmCI0eOoL6+3lB5J8bQ7D46LY7XXnstNm7ciF27dmHDhg1obGxEdXU12tvbdcs7LYZm989p8fvkk0+wYcMGjB07Frt378Zdd92F73znO3jhhReirhMthn6/H11dXanZ0JTVudjInXfeKUaNGiWamppilsvLyxObN28OW/azn/1MjBw5MpWblxRG9zFST0+P+PKXvyzuvffeFG1ZYrq7u8WHH34oDh8+LNasWSMuueQS8d577+mWdWL8zOxfJLvH7vTp02LkyJHif/7nf/qXDXYLw2kxjGcfI9k9jpHa2tqEx+OJejvRaTGMNNj+RbJ7/PLy8kRVVVXYsn/6p38SM2fOjLrO2LFjxcMPPxy27JVXXhEAxIULF1KynRlfM7Jq1Srs3LkTb7zxBiorK2OWLS0tRWtra9iy1tZWlJaWpnITE2ZmHyPl5eVh+vTp+Oijj1K0dYnJz8/HmDFjcPXVV6O+vh5XXnkl1q9fr1vWifEzs3+R7B67d999F+fOncNVV12F3Nxc5ObmYv/+/fjpT3+K3NxcqKo6YB2nxTCefYxk9zhGGjZsGK644oqo2+u0GEYabP8i2T1+ZWVlA2pbJ0yYEPNWVLQYejweFBYWpmQ7MzYZEUJg1apV2L59O373u9/h8ssvH3Sdqqoq7N27N2zZnj17Yt7Dt1I8+xhJVVUcP34cZWVlKdjC5AsEAuju7tZ9zmnx0xNr/yLZPXZz5szB8ePHcezYsf7HjBkzsGzZMhw7dgwul2vAOk6LYTz7GMnucYzU0dGBjz/+OOr2Oi2GkQbbv0h2j9/s2bMHDPnw5z//GaNGjYq6jiUxTEl9iw3cddddori4WOzbt080Nzf3P0KrmG699VaxZs2a/r/feustkZubK/793/9d/PGPfxT33XefyMvLE8ePH7diFwYVzz7+6Ec/Ert37xYff/yxePfdd8W3vvUtUVBQYPjWQDqtWbNG7N+/XzQ2Nor//d//FWvWrBGKoojXX39dCOH8+JndPyfFLprIWxhOj6GewfbRaXH853/+Z7Fv3z7R2Ngo3nrrLTF37lxxySWXiHPnzgkhnB9Ds/vntPj94Q9/ELm5ueKhhx4SH374ofjlL38phgwZIjZt2tRfZs2aNeLWW2/t//uTTz4RQ4YMEd/73vfEH//4R/Gzn/1MuFwusWvXrpRtZ8YmIwB0H88//3x/mRtuuEHcdtttYeu99NJL4oorrhD5+fli0qRJ4pVXXknvhpsQzz7W1dWJyy67TOTn54uSkhJx0003iSNHjqR/4w34+7//ezFq1CiRn58vLr30UjFnzpz+E7UQzo+f2f1zUuyiiTxROz2GegbbR6fFccmSJaKsrEzk5+eLiooKsWTJEvHRRx/1P+/0GJrdP6fFTwghfvOb34jJkycLt9stxo8fL37+85+HPX/bbbeJG264IWzZG2+8IaZNmyby8/PFl770pbDzSiooQgiRunoXIiIiotgyts0IEREROQOTESIiIrIUkxEiIiKyFJMRIiIishSTESIiIrIUkxEiIiKyFJMRIiIishSTESIiIrIUkxEiIiKyFJMRIiIishSTESIiIrIUkxEiIiKy1P8HPXiPueUv7vYAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "gathered_skeletons = np.concatenate(\n",
        "    [targets[:, target_skeleton],\n",
        "     targets[:, target_skeleton] + np.array([0, 1]).reshape(-1, 1),\n",
        "     targets[:, target_skeleton] + np.array([1, 0]).reshape(-1, 1),\n",
        "     targets[:, target_skeleton] + np.array([1, 1]).reshape(-1, 1),\n",
        "    ], axis=1)\n",
        "\n",
        "pt.plot(gathered_skeletons[0], gathered_skeletons[1], \"ro\", alpha=0.05)\n",
        "\n",
        "target_proxies = 4 + 2 * np.array([np.cos(angles), np.sin(angles)])\n",
        "pt.plot(target_proxies[0], target_proxies[1], \"bo\")\n",
        "\n",
        "parent_proxy_mat = interaction_mat(gathered_skeletons, target_proxies)\n",
        "\n",
        "idx, proj = sli.interp_decomp(parent_proxy_mat.T, nproxies)\n",
        "parent_skeleton = idx[:nproxies]\n",
        "\n",
        "pt.plot(gathered_skeletons[0, parent_skeleton], gathered_skeletons[1, parent_skeleton], \"ro\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.4+"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}