{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "4c22782c-348c-409a-b2f6-ec4f5097a562",
      "metadata": {},
      "source": [
        "# Reusing Decompressors\n",
        "\n",
        "Copyright (C) 2026 Andreas Kloeckner\n",
        "\n",
        "<details>\n",
        "<summary>MIT License</summary>\n",
        "Permission is hereby granted, free of charge, to any person obtaining a copy\n",
        "of this software and associated documentation files (the \"Software\"), to deal\n",
        "in the Software without restriction, including without limitation the rights\n",
        "to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
        "copies of the Software, and to permit persons to whom the Software is\n",
        "furnished to do so, subject to the following conditions:\n",
        "\n",
        "The above copyright notice and this permission notice shall be included in\n",
        "all copies or substantial portions of the Software.\n",
        "\n",
        "THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n",
        "IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n",
        "FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n",
        "AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n",
        "LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n",
        "OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n",
        "THE SOFTWARE.\n",
        "</details>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "import matplotlib.pyplot as pt\n",
        "import scipy.linalg.interpolative as sli"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {},
      "outputs": [],
      "source": [
        "sources = np.random.rand(2, 800)\n",
        "targets = np.random.rand(2, 800) + 3\n",
        "\n",
        "all_distvecs = sources.reshape(2, 1, -1) - targets.reshape(2, -1, 1)\n",
        "dists = np.sqrt(np.sum(all_distvecs**2, axis=0))\n",
        "A = 1/dists\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7fd870e51010>]"
            ]
          },
          "execution_count": 31,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAGdCAYAAAAFcOm4AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAT9RJREFUeJzt3XlYVPXiBvB3ZoBhnWEblE0QXCBXciXQMlo0S9Osa6aVS9ly27y/uqGVebtFddu8lWnadUG7tHhTy7JcA1PRXFEEEUVQYGBYZhiWAWbO7w9kalQUcJgzA+/neebR+Z7DmXcimZdzvucciSAIAoiIiIgcmFTsAERERETXi4WGiIiIHB4LDRERETk8FhoiIiJyeCw0RERE5PBYaIiIiMjhsdAQERGRw2OhISIiIofnJHYAWzGZTCgsLISXlxckEonYcYiIiKgVBEFAVVUVgoKCIJW2vB+myxSawsJChIaGih2DiIiI2qGgoAAhISEtLu8yhcbLywtA038QhUIhchoiIiJqDZ1Oh9DQUPPneEu6TKFpPsykUChYaIiIiBzMtaaLcFIwEREROTwWGiIiInJ4LDRERETk8FhoiIiIyOGx0BAREZHDY6EhIiIih8dCQ0RERA6PhYaIiIgcHgsNEREROTwWGiIiInJ4LDRERETk8FhoiIiIyOGx0FynHVlqJP14EmV6g9hRiIiIuiyHKTSCIGDZsmUYM2YMJk2ahD179ogdCQDw9k9ZWJZ6BgfyysWOQkRE1GU5TKFZs2YNjh49ioULF2LMmDEYN24c6urqxI6FYeG+AID9ZytETkJERNR12bTQmEwmaDQaGAwtH54xmUzQ6/WXjT/44INYsmQJbrnlFjzzzDNwcXFBfX19R8ZtleE9mwoN99AQERGJxyaFRqPRICkpCb169YJKpcI333xz2TpGoxHz5s2DQqGAr68vevfuje3bt5uXu7i4mP/+xhtvYO7cuVAoFLaIf1XNheZEoRZVdQ0ipyEiIuqabFJoNm7cCK1Wa1FQLpWUlITk5GSkpaVBr9dj2rRpmDBhAgoKCszrCIKAl19+GdXV1fjnP/9pi+jXFKh0Q6ivG0wCcCi/Uuw4REREXZJNCs3s2bPx9ttvo2fPnldcLggCPvnkEzz99NOIiYmBi4sLFi5cCC8vL6xYsQIAUF9fj+nTp0OpVOKdd9655msaDAbodDqLR0dpnkdz4CwPOxEREYnBLiYF5+XlQa1WIz4+3jwmlUoRFxeH9PR0AMDSpUuxfv16fPXVVxg8eDAGDx6M3NzcFreZlJQEpVJpfoSGhnZY/uHmicEsNERERGJwEjsAAJSWlgIAVCqVxbhKpUJeXh6ApknBo0ePtlgeHBzc4jYTExMxb94883OdTtdhpaZ5Hs2R85WoazDC1VnWIa9DREREV2YXhUYikQBomhj8Z42NjZBKm3YiqVSqywrP1cjlcsjlcuuFvIqe/h7w93SBRl+PY+e15oJDREREtmEXh5yCgoIAAGq12mJcrVYjMDBQjEhtIpFIePo2ERGRiOyi0AQHByMiIgLbtm0zjxkMBqSmpmLUqFEiJmu9YZxHQ0REJBqbHHJqaGiAVqs1P9fr9dBoNHB1dYWnpycAYMGCBXj66acxbNgwxMTE4I033oCbmxvmzJlji4jXrbnQHDxXAaNJgEwqETkRERFR12GTPTR79+5FVFQUoqKi4Ofnh1deeQVRUVF49dVXzevMmjULixcvRlJSEhISEqDVarFz5074+PjYIuJ1iw5UwEvuBL2hESeLOu4UcSIiIrqcRBAEQewQtqDT6aBUKqHVajvsCsOPrtyPXdmlePXuGzA7/srX3CEiIqLWa+3nt13MoekszBODOY+GiIjIplhorKj5AnsH8srRRXZ8ERER2QUWGisaEKKEi5MUZdX1yC2tFjsOERFRl8FCY0VyJxliQr0B8Ho0REREtsRCY2XN82h4PRoiIiLbYaGxMhYaIiIi22OhsbIbe/hAJpXgQmUtLlTWih2HiIioS2ChsTIPuRP6BTWdJ8/Tt4mIiGyDhaYDNJ++vZ8Tg4mIiGyChaYDDOMF9oiIiGyKhaYDNN+oMqdEj/LqepHTEBERdX4sNB3A18MFvQOa7iLO69EQERF1PBaaDjKMp28TERHZDAtNBxnR84/7OhEREVHHYqHpIM3zaE4U6qA3NIqchoiIqHNjoekgQd5uCPFxg9Ek4NC5CrHjEBERdWosNB2o+Xo0POxERETUsVhoOlDzxOB0TgwmIiLqUCw0Haj5RpVHCiphaDSKnIaIiKjzYqHpQBH+HvD3dEF9owkZ57VixyEiIuq0WGg6kEQiMZ/txMNOREREHYeFpoMN48RgIiKiDsdC08Ga59EczKuA0SSInIaIiKhzYqHpYNGBCnjKnVBlaMTJIp3YcYiIiDolFpoOJpNKMCTMBwAPOxEREXUUFhobGM4bVRIREXUoFhobGP6nG1UKAufREBERWRsLjQ0MDFHCxUkKjb4eZzXVYschIiLqdByq0OzcuRN33303xo8fj7S0NLHjtJrcSYbBod4AeNiJiIioIzhMoSkuLsb999+P6dOnY8aMGbjvvvug0WjEjtVqzTeq3M+JwURERFZn00JjMpmg0WhgMBiuuo5er79sfNOmTZgwYQKmTp2KqVOn4q677sIPP/zQkXGtihODiYiIOo5NCo1Go0FSUhJ69eoFlUqFb7755rJ1jEYj5s2bB4VCAV9fX/Tu3Rvbt283Ly8sLER4eLj5eUREBC5cuGCL+FZxY5gPpBLgfEUtirS1YschIiLqVGxSaDZu3AitVmtRUC6VlJSE5ORkpKWlQa/XY9q0aZgwYQIKCgoAAAqFAlVVVeb1dTodFApFh2e3Fk+5E/oHKwFwLw0REZG12aTQzJ49G2+//TZ69ux5xeWCIOCTTz7B008/jZiYGLi4uGDhwoXw8vLCihUrAAAjRozADz/8gJqaGlRXV2PTpk0YOXJki69pMBig0+ksHmJrvq/TL5lqZJzXIkddhfyyGpTo6qCtaUBdg5GndRMREbWDk9gBACAvLw9qtRrx8fHmMalUiri4OKSnpwMA4uLiEBcXhx49egAAHnjgAQwbNqzFbSYlJWHRokUdG7yNhoX74ovdZ7H5WBE2HytqcT25k7Tp4SyDv6cc/546GL27edkwKRERkWOxi7OcSktLAQAqlcpiXKVSmZcBwIoVK5CZmYmTJ09iyZIlV91mYmIitFqt+dF86EpMt/RVISEqAOF+7ghUusLH3RnuLjJIJZbrGRpN0NU1orTKgJNFOqzdd06cwERERA7CLvbQSCRNn+hGo9FivLGxEVKpZecKCAho1Tblcjnkcrl1AlqJq7MMXzx6+V4lQRDQaBJQ12CEodFk/nNPbhle3XAcO7JL8LogmP87ERERkSW72EMTFBQEAFCr1RbjarUagYGBYkSyKYlEAmeZFF6uzvD3lCPExx2RKk/cd2MwXJykKCivxemSy09lJyIioiZ2UWiCg4MRERGBbdu2mccMBgNSU1MxatQoEZOJy93FCbERfgCA7VklIqchIiKyXzYpNA0NDdBoNOYr++r1emg0GosL6C1YsABLlixBSkoKsrOzMXv2bLi5uWHOnDm2iGi3EqKbDrHtOMlCQ0RE1BKbFJq9e/ciKioKUVFR8PPzwyuvvIKoqCi8+uqr5nVmzZqFxYsXIykpCQkJCdBqtdi5cyd8fHxsEdFujenbVGgO5legsqZe5DRERET2SSJ0kQuf6HQ6KJVKaLVah7ogHwDc8eGvOKXWY/HUwZg4OFjsOERERDbT2s9vu5hDQ1d3a1Q3AMAOzqMhIiK6IhYaB9A8j2ZXdikajSaR0xAREdkfFhoHEBPqDW93Z2hrG3C4oFLsOERERHaHhcYBOMmkuLlP01WUt/NsJyIiosuw0DiIW6Munr6dpb7GmkRERF0PC42DuLmPCjKpBKfUehSU14gdh4iIyK6w0DgIb3cXDAlruibPzmwediIiIvozFhoH0nzYifNoiIiILLHQOJCEi4Vmb24Zqg2NIqchIiKyHyw0DqRXgCdCfd1QbzTht9MaseMQERHZDRYaByKRSJBw8arBnEdDRET0BxYaBzPmT/NoushtuIiIiK6JhcbBjOjpC3cXGUqqDDhRqBM7DhERkV1goXEwrs4yxPfyB8CznYiIiJqx0Dig5ptV7uA8GiIiIgAsNA5pTN+mQnO0oBKlVQaR0xAREYmPhcYBBShcMSBYCYBnOxEREQEsNA7LfLNKzqMhIiJioXFUzfNo0nJKUd9oEjkNERGRuFhoHFT/ICX8PeWorjdi/9lyseMQERGJioXGQUmlEtwapQIAbM9Si5yGiIhIXCw0DuzWi7dB2JHFqwYTEVHXxkLjwOJ7+8NFJsW5shqc0VSLHYeIiEg0LDQOzFPuhBERvgB4thMREXVtLDQOrvn0bc6jISKiroyFxsE1F5oDeRXQ1jaInIaIiEgcLDQOLszPA70CPGE0CUjLKRU7DhERkShYaDoBXjWYiIi6OocqNOvXr8edd96JcePGYdOmTWLHsRvNhWZndgmMJp6+TUREXY/DFJq0tDR89dVX+Nvf/obp06djxowZyM3NFTuWXRgS5gOFqxMqahpwpKBC7DhEREQ252TNjZlMJpSXl8PLywtyubzFdWpqauDp6dmmbY8YMQKjRo0yP//www9RXc1rrwCAs0yKm/sG4PujhdiRVYIhYb5iRyIiIrIpq+yh0Wg0SEpKQq9evaBSqfDNN99cto7RaMS8efOgUCjg6+uL3r17Y/v27eblhYWFGDx48BUfhYWFcHFxMa+7bt069O7dGwMHDrRG/E6h+TYIP2YUQ6M3iJyGiIjItqyyh2bjxo3QarXYvn07IiIirrhOUlISkpOTkZaWhn79+uHNN9/EhAkTkJWVhdDQUPj7+2PVqlVX/Fp/f3/z35cuXYrU1FSsWbPGGtE7jVv6BMBL7oSzmmrc9sGvmH9XNO4fEgKJRCJ2NCIiog4nEax8EyCJRILk5GRMnz7dPCYIAgIDA/HEE0/g9ddfB9B06CkoKAhz587FokWLrrldQRCwYMEClJWV4bPPPoNUevWdSwaDAQbDH3sqdDodQkNDodVqoVAo2vfm7FzGeS1eWn8MJ4t0AICREb54a9IARKjadniPiIjIXuh0OiiVymt+fttkUnBeXh7UajXi4+P/eGGpFHFxcUhPT2/VNn766Se8/fbbSE9Px4033ojBgwdj586dLa6flJQEpVJpfoSGhl73+7B3A0KU+P6vcZh/VxRcnaXYd6YcYxen4ePtOahvNIkdj4iIqMNYdVJwS0pLmy74plKpLMZVKhXy8vJatY24uDgcOnTIYqxnz54trp+YmIh58+aZnzfvoensnGRSPD46EuP6B2LBhuNIPVWK97eewqajhUiaPABDwzlhmIiIOh+bFJrmeRxGo9FivLGx8ZqHjpoplUoMHjy41a8pl8tbPNOqKwj1dcfqmcOw6Wgh/vF9JnJK9JiydC+mjeiBv4+NgtLNWeyIREREVmOTQ05BQUEAALXa8gaKarUagYGBtojQJUkkEkwcHIztf7sZfxnatHfqy/R83PbBr9h8rAhWnj5FREQkGpsUmuDgYERERGDbtm3mMYPBgNTUVItry1DH8HZ3wTtTBiLl8ZGI8PdAaZUBT395CLNX/44SXZ3Y8YiIiK6bVQpNQ0MDNBoNNBoNAECv10Oj0UCv15vXWbBgAZYsWYKUlBRkZ2dj9uzZcHNzw5w5c6wRgVphZIQffnxuFJ5N6A1nmQQ7skowZeleFJTXiB2NiIjouljltO3U1FRMnjz5svEZM2bgww8/ND///PPP8emnn6KsrAwxMTF49913ER0dfb0v3yqtPe2rq8hRV2H26t+RX16Dbgo51s0ZgV4BXmLHIiIistDaz2+rX4fGXrHQXE6tq8P0FenIKdHD18MFa2YNR/9gpdixiIiIzOzqOjRkn7opXPHV3FgMDFGivLoeD36+DwfyysWORURE1GYsNF2cr4cL1s0ZgeE9fVFlaMSML9Lx66lSsWMRERG1CQsNwcvVGWtmDceYvirUNZgwZ/UB/JRRJHYsIiKiVmOhIQCAq7MMy2YMxfiBgWgwCnj6y0P45vcCsWMRERG1CgsNmbk4SfHvqTGYOiwUJgF48dtjWPXbWbFjERERXRMLDVmQSSVImjwAc+Kb7pP1+veZ+Hh7Dq8qTEREdo2Fhi4jkUiwYHw0XritDwDg/a2nkPRTFksNERHZLRYauiKJRILnbuuN1+6+AQDweeoZvLLhOEsNERHZJRYauqpZ8T3x7n0DIZUA69LzsfdMmdiRiIiILsNCQ9f0wLBQTB8ZBgBY9usZkdMQERFdjoWGWmVOfASkEuDXU6U4WaQTOw4REZEFFhpqlR5+7rhrQCCApvk0RERE9oSFhlpt7uhIAMCmo4U4X1EjchoiIqI/sNBQqw0IUSKulx+MJgH/2Z0ndhwiIiIzFhpqk+a9NCkH8lFZUy9yGiIioiYsNNQmo3r7IzpQgZp6I9buOyd2HCIiIgAsNNRGEokET9wcAQBYtScPdQ1GkRMRERGx0FA73DUgEMHebtDo67H+0Hmx4xAREbHQUNs5y6SYM6rp5pXLU8/AaOLtEIiISFwsNNQufxkWCm93Z+SV1eDnE8VixyEioi6OhYbaxd3FCQ+bb4eQy5tWEhGRqFhoqN0evikccicpjp7XYt+ZcrHjEBFRF8ZCQ+3m7ynH/UNDAADLUnNFTkNERF0ZCw1dl+abVu7K5k0riYhIPCw0dF3C/T0wrn/TTSuX86aVREQkEhYaum6Pj2660N6mo4W4UFkrchoiIuqKWGjoug0K9UZshB8aTQL+s/us2HGIiKgLYqEhq5h78XYI/92fD21Ng8hpiIioq3G4QlNbW4uRI0di0qRJYkehP7m5jwpR3b2ablqZzptWEhGRbTlcofn73/+Oe++9F7W1nKthTyQSiXkvzcrfzvKmlUREZFNWLTSVlZXYt28fNBpNi+ucP38eR44cQU1NTZu3/+233yIqKgpDhw69npjUQe4eGIQgpSs0+nr879AFseMQEVEXYpVCk5OTg9mzZyM6OhqxsbHYsmXLZevU1NRgwoQJ6NOnD6ZMmYJu3bph1apV5uX5+fnw9PS84iM/Px/nzp3D1q1b8dRTT1kjMnUAZ5kUs0c17aVZnsabVhIRke1YpdBkZGQgNjYWubktXy12/vz5OH78OPLy8nD69Gl88sknmDNnDk6cOAEACA0NRXFx8RUfoaGh+Pjjj7FmzRp4enrinnvuwbZt2zBu3DhrxCcrmjosFEo3Z5zVVOPdn7OQeqoU5ytqYGK5ISKiDiQRrHxXQYlEguTkZEyfPt08ZjQa4efnhwULFuDFF180j/fs2RNTpkzBv/71r2tut76+HvX19QCAXbt2YfHixdi0aRPc3NyuuL7BYIDBYDA/1+l0CA0NhVarhUKhaO/bo1Z47+dsfLLztMWY3EmKcD8PRKg80NO/6RGh8kCEvyd8PFxESkpERPZOp9NBqVRe8/PbyRZhcnNzodVqMWzYMIvx4cOH4/Dhw63ahouLC1xcmj74XF1dIZPJWiwzAJCUlIRFixa1PzS121NjIiGVACeLq3BWU41zZdUwNJqQra5CtrrqsvW93Z0RF+mPpPsGQOHqLEJiIiJydDYpNOXlTXdi9vX1tRj38/PDqVOn2ry9MWPGID4+/qrrJCYmYt68eebnzXtoqOO5uzhh3h19zc8bjSYUVtYhV6PH2dJqnNU0Pc6U6lGorUNlTQM2ZxQhxMcNiXdFi5iciIgclU0KjbNz02/dfz4EBDRdU6Z5WVvIZDLIZLKrriOXyyGXy9u8bbI+J5kUPfzc0cPPHWP6Wi6rrTdiy4kivPDVUazck4fpI8MQ6usuTlAiInJYNrkOTVhYGACgsLDQYrywsNC8jLomNxcZ7h0cjJsi/VDfaMJ7v2SLHYmIiByQTQqNv78/Bg4ciO+//948VllZidTUVCQkJNgiAtkxiUSC+XdFQyIBNh4pxNGCSrEjERGRg7HKISedTofMzEzz89zcXOzbtw8BAQGIiGi6Lslbb72FiRMnIjw8HDExMXjvvfcQHh6ORx55xBoRyMH1D1ZiUkww/nfoAt788SS+enwkJBKJ2LGIiMhBWOW07QMHDuCZZ565bHzixIlITEw0P9+6dSuWLl2KsrIyxMTEYP78+VCpVNf78q3S2tO+SDyFlbUY894uGBpNWP7wUNx+QzexIxERkcha+/lt9evQ2CsWGsfw7pYsLNmViwiVB35+fjScZQ53uzEiIrKi1n5+89OC7MqTt0TCz8MFZ0qrkbI/X+w4RETkIFhoyK54uTrj+dt6AwA+2paDqroGkRMREZEjYKEhuzN1eA9E+HugrLoeS39t+f5gREREzVhoyO44y6R4eVwUAGBF2lkUVtaKnIiIiOwdCw3Zpdtv6IbhPX1h4MX2iIioFVhoyC5JJBIsuHhfp+8OX8DxC1qRExERkT1joSG7NSjUGxMGBUEQgLd+PIkucoUBIiJqBxYasmsv3tkXLjIp9uSWYVd2qdhxiIjITrHQkF0L9XXHzLhwAE17aRqNJnEDERGRXWKhIbv31Jhe8HZ3Rk6JHl//fl7sOEREZIdYaMjuKd2c8eytTRfb+2DrKVQbGkVORERE9oaFhhzC9JFhCPNzh0ZvwLLUM2LHISIiO8NCQw7BxUmKl8c2XWzv89RcFGvrRE5ERET2xEnsAEStNbZ/dwwJ88HBcxWY+OluRKo8Eah0Q7C3KwK93RCodEWQtxuCvN3gKef/2kREXQl/6pPDkEgkeGV8NO5fuhdqnQFqnaHFdb1cnRCkdEOQtyv6BSnx5C2R8GDJISLqtCRCF7lamU6ng1KphFarhUKhEDsOXYeSqjqcKa1GkbYWhZV1KKysRZG26c/Cylro6i6fNDxlSAjeu3+QCGmJiOh6tPbzm7+yksMJ8HJFgJdri8v1hkYUVdaiUFuHHHUV3vzxJL49eB7j+ndHQnQ3GyYlIiJb4aRg6nQ85U7o3c0LN/dRYc6oCMyJ7wkASPxfBipr6kVOR0REHYGFhjq9v93RFxEqD5RUGfD6phNixyEiog7AQkOdnquzDO/fPwhSCbDhSCG2HC8WOxIREVkZCw11CTE9fDD35kgAwCsbMlBezUNPRESdCQsNdRnP39Ybfbp5QqOvx6sbj4sdh4iIrIiFhroMuZMM798/GDKpBJuPFeGHY4ViRyIiIithoaEuZUCIEk/f0nTo6dUNx1Fa1fLF+YiIyHGw0FCX89dbeyOquxcqahrwyoYMdJFrSxIRdWosNNTluDhJ8f4Dg+AkleDnE2psOspDT0REjo6FhrqkfkFKPJvQGwDw2sYTUOt4924iIkfGQkNd1pO3RKJ/sALa2gbM/x8PPREROTKHKjTV1dV4+eWXER8fj8GDByMjI0PsSOTAnGVSvH//YLjIpNieVYL1hy6IHYmIiNrJoW5Oed9998HZ2RlvvvkmlEolIiMjxY5EDq5vdy88f3tvvLslG4u+P4G4Xn4IVLqJHYuIiNrIqoXGZDKhvLwcXl5ekMvlLa5TU1MDT0/PNm07MzMTe/bswYULF+Dl5WWNuEQAgMdHReDnE2ocLajE39dnYPXMYZBIJGLHIiKiNrDKISeNRoOkpCT06tULKpUK33zzzWXrGI1GzJs3DwqFAr6+vujduze2b99uXl5YWIjBgwdf8VFYWIjTp0/jxhtvxKJFixAXF4cXX3wRNTU11ohPXZyTTIr37x8IFycpUk+V4oOtp3D8ghbVhkaxoxERUStZZQ/Nxo0bodVqsX37dkRERFxxnaSkJCQnJyMtLQ39+vXDm2++iQkTJiArKwuhoaHw9/fHqlWrrvi1/v7+8PLywuHDh/HYY49h8uTJePvtt5GUlIQ33njDGm+BurheAV548Y6+ePPHk/h4x2l8vOM0ACBQ6YoIlQciVZ6I8PdAhMoTkQGeCFS4QirlXhwiInshEax8aodEIkFycjKmT59uHhMEAYGBgXjiiSfw+uuvA2g69BQUFIS5c+di0aJF19xuRUUFhg4dimPHjsHDwwPvv/8+MjMz8cUXX1xxfYPBAIPhj6vA6nQ6hIaGQqvVQqFQXN+bpE7JaBLw2a7TSD2lQW6pHmVXuYGlq7MUEf6e6BekwIt39kWAwtWGSYmIug6dTgelUnnNz2+bTArOy8uDWq1GfHy8eUwqlSIuLg7p6emt2oaPjw9eeuklhIeHw9fXF3V1ddi8eXOL6yclJbWqKBE1k0kl+OutvfHXW5uuT1NZU4/c0mqcKdX/6U89zpXVoK7BhMwiHTKLdLhQWYt1c0Zw3g0RkYhsUmhKS0sBACqVymJcpVIhLy+v1duZO3cuHnroIRQXFyM8PBxOTi3HT0xMxLx588zPm/fQELWWt7sLhoS5YEiYj8V4g9GEgvIaZBVX4YWvjmBPbhm+O3wBk28MESkpERHZpNA0/+ZqNBotxhsbGyGVtm1esqenJ3r16nXN9eRyeYtnWhFdD2eZFBEqT0SoPJFXVo13t2Tjn5tPYkzfAPh4uIgdj4ioS7LJhfWCgoIAAGq12mJcrVYjMDDQFhGIOsRjoyLQp5snyqvr8fZPWWLHISLqsmxSaIKDgxEREYFt27aZxwwGA1JTUzFq1ChbRCDqEM4yKd6aNAAA8NXvBdh/tlzkREREXZNVCk1DQwM0Gg00Gg0AQK/XQ6PRQK/Xm9dZsGABlixZgpSUFGRnZ2P27Nlwc3PDnDlzrBGBSDRDw33x4PCm+Vnzv8tAfaNJ5ERERF2PVQrN3r17ERUVhaioKPj5+eGVV15BVFQUXn31VfM6s2bNwuLFi5GUlISEhARotVrs3LkTPj4+V9kykWP4+9go+Hu64HSJHsvTzogdh4ioy7H6dWjsVWvPYydqrw2HL+D5r45A7iTFz8+PRri/h9iRiIgcXms/vx3qbttE9mzi4CDE9/KHodGEVzceRxf5XYGIyC6w0BBZiUQiwT/v7Q8XJynScjTYdLRQ7EhERF0GCw2RFYX7e+CZMU3XSXrjh0xoaxpETkRE1DWw0BBZ2eM3RyBS5QGNvh7v/Mxr0xAR2QILDZGVyZ1k5mvTfJmej4PneG0aIqKOxkJD1AFGRPjhgaFN93aa/7/jaDDy2jRERB2JhYaogySOi4avhwuy1VVYkXZW7DhERJ0aCw1RB/HxcMGCu6IBAIu3n0JBeY3IiYiIOi8WGqIONPnGYMRG+KGugdemISLqSCw0RB1IIpHgn5P6w0Umxa7sUvyYUSx2JCKiTomFhqiDRao88eQtkQCAhZuOI0ddJXIiIqLOh4WGyAaevCUS0YEKaPT1+Mvn+3D8glbsSEREnQoLDZENuDrL8OWcERgYokR5dT0e/HwfDuTx+jRERNbCQkNkIz4eLlg3ZwSG9/RFlaERM75Ix6+nSsWORUTUKbDQENmQl6szVs8cjlv6qlDXYMKc1Qew5XiR2LGIiBweCw2Rjbm5yPD5jKEYPyAQDUYBT607hG8Pnhc7FhGRQ2OhIRKBi5MU/34wBg8MDYFJAP7vm6NYvSdP7FhERA6LhYZIJDKpBG9PHoiZceEAgIWbTuDTnad58T0ionZgoSESkVQqwWt334BnE3oDAP71czbe3pLFUkNE1EYsNEQik0gkmHd7H/N9n5b9egavbDgOk4mlhoiotVhoiOzEY6MjkDR5ACQSYF16PuZ9fQR1DUaxYxEROQQnsQMQ0R8eHN4DHnInzPvqCDYcKcSWE8UYEuaD2Ag/xEb6Y2CIEs4y/h5CRHQpFhoiOzNhUBA85TIk/i8Dap0Bv50uw2+nywCcgruLDMPCfREb6YfYCD/0D1ZCJpWIHZmISHQSoYvMPtTpdFAqldBqtVAoFGLHIbomQRBwukSPPbll2Jtbhn1ny1BZ02CxjperE0b09MXICD/c0jcAvQI8RUpLRNQxWvv5zUJD5CBMJgFZxVXYe6YMe3M1SD9bjqq6RvNyiaTpkNXf74yC0t1ZxKRERNbDQnMJFhrqbIwmAScKtdibW4bdpzVIy9EAAPw8XJB4VzTuuzEYEgkPRxGRY2OhuQQLDXV2e3PL8OrG4zhdogcADA/3xT8n9Uefbl4iJyMiar/Wfn7zdAmiTiI20g8/PjsKfx8bBTdnGfbnleOuxWlI+vEkqg2N194AEZEDc6hCs2nTJtx5550YNWoUkpKSeDVVoku4OEnx5C2R2DpvNG6/oRsaTQKWpZ7B7R/8ii3Hi/lvhog6LYc55KRWqxEVFYX//Oc/8PX1xbPPPosFCxbggQceaNXX85ATdUXbMtV4/fsTOF9RCwAY01eFRRP6o4efu8jJiIhaR5RDTiaTCRqNBgaD4arr6PX6Nm+7sbERPj4+GD9+PG6++WYMGDAADQ0N1/5Coi7sthu6YesLN+OvY3rBWSbBzuxS3P7hr/h4ew4KK2vRYDSJHZGIyCqssodGo9Fg+fLlWL58Oc6ePYvk5GRMnz7dYh2j0YgXX3wRn3/+Oerr6xEWFoalS5ciISEBAFBYWIi77rrritv/8ccfERQUhA8++ADvvPMO3N3dERMTg2+++QYymaxVGbmHhrq60yV6vLbxOPbklpnHJJKms6ICvFzRTSH/40+FK7opXBHgJUc3hSv8PV3gxCsUE5EIWvv5bZUrBW/cuBFarRbbt29HRETEFddJSkpCcnIy0tLS0K9fP7z55puYMGECsrKyEBoaCn9/f6xateqKX+vv74+srCwsXboUK1euhEKhwKJFi/Dll19ixowZ1ngLRJ1erwBPrJszApuOFmLx9hzkl9Wg0SRAo6+HRl+PzKKWv9bb3Rnrn7wJkSpeuI+I7JPV59BIJJLL9tAIgoDAwEA88cQTeP311wE0HXoKCgrC3LlzsWjRomtud8OGDfjwww+xc+dOSKVSLFy4EDqdDh9++OEV1zcYDBaHvnQ6HUJDQ7mHhugik0lAeU091Lo6lFQZUKKrg1pnuOx5qd4Ao0nA3QMD8cm0G8WOTURdjE330FxLXl4e1Go14uPjzWNSqRRxcXFIT09v1TZuv/12vPfeewgKCoKnpycMBgO2bNnS4vpJSUmtKkpEXZVUKoG/pxz+nnL0u8p6J4t0GLc4DZszivCcugq9eV0bIrJDNjkoXlpaCgBQqVQW4yqVyrzsWjw8PLB7924cOHAAW7ZswdmzZ9GvX8s/hhMTE6HVas2PgoKC9r8Boi4sOlCBO/t1gyAA/95xWuw4RERXZJNC03z5daPRaDHe2NgIqbRtEUJDQ9GrVy84OV1955JcLodCobB4EFH7PJvQGwDww7FCnC6pEjkNEdHlbFJogoKCADRdS+bP1Go1AgMDbRGBiK5DvyAl7rihaS/Nx9xLQ0R2yCaFJjg4GBEREdi2bZt5zGAwIDU1FaNGjbJFBCK6Ts17aTYdLTTfL4qIyF5YpdA0NDRAo9FAo2m6269er4dGo7G4gN6CBQuwZMkSpKSkIDs7G7Nnz4abmxvmzJljjQhE1MH6Bytx+8W9NJ/syBE7DhGRBasUmr179yIqKgpRUVHw8/PDK6+8gqioKLz66qvmdWbNmoXFixcjKSkJCQkJ0Gq12LlzJ3x8fKwRgYhs4Lk/7aXJLeVeGiKyHw5zL6frxSsFE1nHnNW/Y9tJNSbFBOPDvwwWOw4RdXKi3MuJiDq/529r2kuz8cgFnOFeGiKyEyw0RNQm/YOVuC06ACYB+IRnPBGRnWChIaI2ey6hDwBgw5ELOKupFjkNERELDRG1w4AQJRKimvbSfMwznojIDrDQEFG7PGeeS1OIPO6lISKRsdAQUbsMDPHGrVEBMJoEXj2YiETHQkNE7dZ8XZoNRy5wLw0RiYqFhojabVCoN8b0VcFoEvDJTu6lISLxsNAQ0XV57ramM56+O3wB58q4l4aIxMFCQ0TXZXCoN265uJfmU+6lISKRsNAQ0XVrnkuz/tAF5JfViJyGiLoiFhoium4xPXwwug/30hCReFhoiMgq/thLcx4F5dxLQ0S25SR2ACLqHIaE+WBUb3+k5WgwN/kgRvXxRy+VJ3oFND28XJ3FjkhEnRgLDRFZzQu398Hu0xpkFumQWaSzWNZd4WouN39++Hm4QCKRiJSYiDoLiSAIgtghbEGn00GpVEKr1UKhUIgdh6jTOn5Bi8P5FThdokdOiR6nS/QoqTK0uH6g0hX/mjII8b39bZiSiBxFaz+/WWiIqMNpaxtwukSP3BI9TpfqkaOuwulSPc5X1EIQALmTFCtnDsNNkSw1RGSJheYSLDRE9qfa0Ihn/nsYO7JK4OYsw6qZwzAiwk/sWERkR1r7+c2znIhINB5yJyx56EaM7qNCbYMRM1cdwMFz5WLHIiIHxEJDRKJydZbh8xlDENfLDzX1RjzynwM4UlApdiwicjAsNEQkOldnGVY8PAwjevpCb2jEjC/SkXFeK3YsInIgLDREZBfcXGT4z6PDMCzcB1V1jZj+RTpOFLLUEFHrsNAQkd3wkDth5czhiOnhDW1tA6avSEd2cZXYsYjIAbDQEJFd8ZQ7YfWs4RgYokRFTQMeWrEPp0tYaojo6lhoiMjuKFydkTxrBPoFKaDR1+PB5enILdWLHYuI7BgLDRHZJaW7M9bOHoGo7l4orTJg2vJ9yNNUix2LiOwUCw0R2S0fDxesnTMCvQM8odY1lRreyZuIroRXCiYiu1dSVYepn+/DmdJquMikGBLmg/je/ojv5Y/+wUrIpLy5JVFn5dC3Phg+fDhqamowaNAgrFu3zjxeXFyMF198EWfPnsWECRPw0ksvtXqbLDREjk2tq8Ps1Qdw/ILlXby93Z1xU6Qf4nupEN/LHz383EVKSEQdwaELzcmTJ5GZmYl//etf2Ldvn3l8zJgxGD58OO655x787W9/w1NPPYVHHnmkVdtkoSFyfIIg4KymGrtPa5CWo8G+3DJUGRot1unh64743v4Y1csfsZF+8HZ3ESktEVlDhxYak8mE8vJyeHl5QS6Xt7hOTU0NPD0927p5AMCRI0fwxBNPmAtNWVkZ+vTpg5KSEshkMmzevBlLlizB5s2bW7U9FhqizqfRaMLR81rsztHgt9MaHMqvQKPJ8kdasLcbIlQeiPD3QGSAJyL8PRGh8kCg0hUSCQ9VEdm71n5+O7VloxqNBsuXL8fy5ctx9uxZJCcnY/r06RbrGI1GvPjii/j8889RX1+PsLAwLF26FAkJCQCAwsJC3HHHHVfc/i+//IKgoKArLisqKkJISAhkMhkAoGfPnigsLGxLfCLqZJwuzqcZEuaD527rDb2hEelnypB2seDklOhxobIWFyprkZajsfhaN2dZU9FReSLC3wMRKg9EdVegd4AnpJyTQ+Rw2lRoNm7cCK1Wi+3btyMiIuKK6yQlJSE5ORlpaWno168f3nzzTUyYMAFZWVkIDQ2FSqVCSkrKFb9WpVK1+NoKhQI63R/HzrVaLZRKZVviE1En5yl3QkJ0NyREdwMAlFfX40ypHrmlepwprUZuaTXOlOqRX16D2gYjThTqcKLQck6Ol6sTYnr4YOjFojQ41Bse8jb9qCQiEbTpX+ns2bOvulwQBHzyySd4+umnERMTAwBYuHAhli1bhhUrVmDRokVwdnZG//792xw0JCQERqMRe/fuRWxsLNatW4fY2NgW1zcYDDAYDObnfy5DRNQ1+Hq4wNfDF0PDfS3GG4wm5JfXXCw5+oulpxoni3SoqmtE6qlSpJ4qBQBIJUB0oMK8J2hImA+Cvd14uIrIzlj11468vDyo1WrEx8ebx6RSKeLi4pCent7q7UyZMgVHjhxBYWEh+vfvj3/84x+YPHkyPv74Y4wbNw6enp5QqVT45ZdfWtxGUlISFi1adF3vh4g6J2eZFJEqT0SqPHE7upnHG40mZBVX4eC5Cvx+rgKHzlXgQmWteU/Omr3nAADdFa4YEu6Debf3QaSqffMEici6rFpoSkubfqO59NCRSqVCXl5eq7fzzjvvoLa21vw8ODgYADBx4kSUlJRArVYjNDT0qttITEzEvHnzzM91Ot01v4aIujYnmRT9g5XoH6zEIzeFAwCKtLU4eK7C/DhRqEOxrg6bjxXhZKEOW54fDRcnXqOUSGxWLTTNu2CNRqPFeGNjI6TS1v+Dj4yMbHGZi4tLq4qJXC5v8QwsIqLWClS64e6Bbrh7YNMJCzX1jThaoMUz/z2EM5pqrNmbhzmjrjynkIhsx6q/VjSfoaRWqy3G1Wo1AgMDrflSRESicHdxQmykH/52R18AwOLtOSjTG67xVUTU0axaaIKDgxEREYFt27aZxwwGA1JTUzFq1ChrvhQRkageGBqK6EAFquoa8cHWU2LHIery2lRoGhoaoNFooNE0Xc9Br9dDo9FAr9eb11mwYAGWLFmClJQUZGdnY/bs2XBzc8OcOXOsm5yISEQyqQQL77kBAPDf/fk4WcQzKYnE1KZCs3fvXkRFRSEqKgp+fn545ZVXEBUVhVdffdW8zqxZs7B48WIkJSUhISEBWq0WO3fuhI+Pj9XDExGJaWSEH8b17w6TALzxQybs8E4yRF2GXd7LqSPw1gdE1BEKymuQ8MGvqG80YdmMIbizX3exIxF1Kq39/Oa5hkRE1yHU1x2PjeoJAHjrx5MwNBqv8RVE1BFYaIiIrtNTt/RCgJcc58pqsPK3PLHjEHVJLDRERNfJQ+6El8ZGAQA+2XEapVU8jZvI1lhoiIisYHJMMAaFKKE3NOK9n7PFjkPU5bDQEBFZgVQqwWsXT+P++mABjl/QipyIqGthoSEispIhYb6YMCgIggAs+v4ET+MmsiEWGiIiK3p5XBRcnaU4kFeBzRlFYsch6jJYaIiIrCjI2w1zRzfdYDfpxyzUNfA0biJbYKEhIrKyJ26ORKDSFRcqa7E89YzYcYi6BBYaIiIrc3OR4eVxTadxL9mVi2JtnciJiDo/Fhoiog4wYVAQhoT5oLbBiHe3ZIkdh6jTY6EhIuoAEokEr93ddBr3/w5fwOH8CpETEXVuLDRERB1kUKg37rsxBADw+veZ0NY0iJyIqPNioSEi6kAvje0LdxcZjhZUYkTSNrz07VFknOdF94isTSJ0kSs/tfb240RE1pZ6qhRv/XgSWcVV5rFBIUpMHxmGewYFwdVZJmI6IvvW2s9vFhoiIhsQBAGH8iuQvPccfswoRr3RBABQujnjgaEheGhEGML9PUROSWR/WGguwUJDRPZCozfg698LsG5fPi5U1prHR/dRYcbIMNwaFQCZVCJiQiL7wUJzCRYaIrI3RpOAXdklSN53Dr+eKkXzT+Ngbze8NLYvJg4OFjcgkR1o7ee3kw0zERHRn8ikEiREd0NCdDecK6vGl+n5+Or3AlyorMVzKUdQpq/HrPieYsckcgg8y4mIyA6E+Xkg8a5o7EtMwJyLJeYfP2Ti052nRU5G5BhYaIiI7IirswwLxkfjhdv6AAD+9XM23tmShS4yO4Co3VhoiIjsjEQiwXO39caCu6IBAJ/tysXrm07AZGKpIWoJCw0RkZ16bHQE3pzUHxIJsHrvOby0/hiMLDVEV8RCQ0Rkxx4aEYYPHhgEmVSCbw+ex7P/PYz6RpPYsYjsDgsNEZGdmxQTgk+n3QhnmQSbM4rwxNqDqGswih2LyK6w0BAROYCx/btjxSPD4OosxY6sEsxceQDVhkaxYxHZDRYaIiIHcXMfFVbPHA5PuRP2ninDjC/Soa3lHbyJABYaIiKHMiLCD+vmjIDSzRmH8ivx4Of7UKY3iB2LSHR2eeuD4cOHo6amBoMGDcK6devM43v27EFSUhIqKysxefJkvPDCC63eJm99QESdyckiHWZ8kQ6Nvh49fN0xvKcvvN2c4e3uDKW7C7zdnOHj7tL0/OK4p9wJEgnvEUWOxaFvfbB69WpkZmbiX//6l3msoqICiYmJ+L//+z+4urpi7ty5CA8Px6RJk0RMSkQkjuhABb6eG4uHVqQjv7wG+eU11/wamVQCbzdn3NGvO16fcAPkTjIbJCWyjXYVGqPRiNLSUnh7e8PV1fWK6zQ2NkKv18Pb27vN24+OjobBYLkLValUYufOnZBKm46S3XrrrdDr9W3eNhFRZxGh8sQPz8Rj20k1yqrroa1pQGVNAypr61FR09D0/OLf6xtNMJoElFXX47/783G+ogbLZgyBu4td/l5L1GZt+j+5pKQEy5YtwxdffIFz584hOTkZ06dPt1jHaDTihRdewIoVKyAIAgIDA7F06VLccccdAIDCwkLz3y/1yy+/ICgo6IrLmosMAOzbtw+nTp3CJ5980pb4RESdjp+nHH8Z1uOa69U1GFFRU4+jBZWY9/VRpOVoMH1FOlY+OhxKd2cbJCXqWG0qNJs3b4bBYEBqairCwsKuuM6bb76JlJQUpKenIzo6GklJSbj33ntx8uRJhIWFQaVSISUl5Ypfq1Kprpnhl19+wbvvvouNGze2uHeIiIgsuTrLEKh0Q6DSDQEKV8xceQCH8ivxl8/3Ys2s4QhQ8OcpObY2FZqZM2dedbkgCPj000/x9NNPY8CAAQCABQsW4NNPP8UXX3yBf/zjH3B2dkb//v3bFXblypX4+uuvsWHDBnh6el51XYPBYHHYSqfTtes1iYg6mxt7+ODrubGY8UU6soqrcP+yvVg7ewRCfd3FjkbUblY9bTsvLw8lJSWIj4//4wWkUsTFxWH//v2t3s6UKVMwZcoUHDt2DP3798f//vc/nD9/HrNmzUJubi5GjhyJ/v37Y8mSJS1uIykpCUql0vwIDQ29rvdGRNSZ9O3uhW+fuAk9fN1xrqwG9322B6fUVWLHImo3q84GKy0tBQD4+/tbjKtUKpw9e7bV23nnnXdQW1trfh4cHAxPT09kZGRYrBcQENDiNhITEzFv3jzzc51Ox1JDRPQnPfzc8e0TsZjxxX5kq6vwwLK9WPnoMMT08BE7GlGbWbXQNF/foLHR8nLcjY2NkMlaf3pgZGTkFcfbcqhKLpdDLpe3en0ioq4oQOGKr+aOxMxVB3A4vxIPrUjH8oeHIq6X/7W/mMiOWPWQU3BwMABArVZbjKvV6hbPXiIiInF5u7tg7ewRGNXbHzX1RsxceQBbjheLHYuoTaxaaIKCghAZGYmtW7eax+rq6vDrr79i9OjR1nwpIiKyIg+5E1Y8MhTj+ndHvdGEp9YdxDe/F4gdi6jV2lRo6uvrUVxcjOLipuau1WpRXFxscQbRK6+8gs8++wxr167FiRMnMHPmTHh6emL27NnWTU5ERFYld5Lh4wdj8MDQEJgE4MVvj2FF2hkYTXZ3hxyiy7TpXk5paWm4//77Lxt/6KGH8P7775uff/HFF1iyZAnKysoQExODd955B3369LFO4nbivZyIiFpHEAS89eNJLE9rOpnDSSpBkLcbQn3dEOLtjlBfN4T6uiPEp+nvKk857xFFHaa1n992eXPKjsBCQ0TUeoIg4PPUM/hw2ynUNZiuuq6rsxQhPu4I8XHDTZF+mBXXE04yq85ooC6MheYSLDRERG1nNAkoqapDQXktCsprUFBR0/T3ihqcL69Bka4Ol36KDAv3weKpMQjydhMnNHUqLDSXYKEhIrK++kYTirS1KCivRVaxDh9ty4He0Ahvd2e8N2UQbruhm9gRycGx0FyChYaIqOOdK6vGX788jIwLWgDA7Pie+PvYKLg48RAUtU9rP7/5fxgREVlNmJ8Hvn0yFrPiegIAvth9FlOW7sG5smqRk1Fnx0JDRERWJXeS4bV7bsDyh4dC6eaMY+e1uPvfu/HDsUKxo1EnxkJDREQd4vYbuuHH50ZhSJgPqgyN+OuXhzH/uwzUNRjFjkadEAsNERF1mGBvN6Q8PhJP3RIJiQT4Mj0f9376G06X6MWORp0MCw0REXUoZ5kUL42NwuqZw+Hv6YKs4irc8/FufHvwvNjRqBPhWU5ERGQzJbo6PP/VEezJLQMAhPq6ITbCD7GRfoiN8Ed3pavICcne8LTtS7DQEBHZB6NJwJKdp/HxjtOoN1pehTjC3wMjI/0QG+GHkRF+UHnJRUpJ9oKF5hIsNERE9kVvaMSBvHLsyy3D3jNlOH5Bi0vvg9k7wPPi3pumguPj4SJOWBINC80lWGiIiOybtrYBB86WY8/FgnOySGexXCIBBgYrMaq3CqN6+yOmhw8v2NcFsNBcgoWGiMixVFTXI/1sGfZeLDin1JZnRnm4yBAb6Y/RffwxqrcK4X7uvOt3J8RCcwkWGiIix6bW1WF3jgapOaXYnaNBWXW9xfIQHzeM7qPC6N7+iI30h9LNWaSkZE0sNJdgoSEi6jxMJgGZRTqk5pQi7ZQGv58rR4Pxj48zqQRQujnD3cUJHnKZ5Z8uMrjLL/55cVzp5owgbzcEe7shyNsNrs4yEd8d/RkLzSVYaIiIOq9qQyPSz5Yh9ZQGaTmlyC29vntHqbzkCPZ2Q4iPG4J93BDi444Q76a/B3u7wUPuZKXkdC0sNJdgoSEi6jpKqwyorKlHdb0RNYbGpj/rG1FtsPxTb2hETb0RFTX1uFBRiwuVtaipv/atGQYEK/HQiB6YMDgI7i4sNx2JheYSLDRERHQtgiCgoqYBFypqcb6iBhcqa3G+oulxobIWFypqoKtrNK/v5eqE+24MwUMjeqB3Ny8Rk3deLDSXYKEhIiJrKK0y4H+HzuPL/fk4V1ZjHh/e0xfTR4ZhbL/uPJ3cilhoLsFCQ0RE1mQyCdh9WoO1+85h20m1+aKA/p4uuH9oKKYN74FQX3dxQ3YCLDSXYKEhIqKOUqStRcr+AqQcyIdaZwDQdCHAW/qoMHV4DwwN84GfJ2/j0B4sNJdgoSEioo7WYDRh+0k11qXnIy1HY7EswEuO6EDFxYcXbghUoKe/B5xkPDx1NSw0l2ChISIiWzqrqcaX6efwS6baYq7Nn8mdpOjTzQtR3b3MZSdC5QE/DxcWnYtYaC7BQkNERGLRGxqRXaxDZlEVThbpkFWkQ1ZxVYuniEslgL+nHN0UruimkCNA4YpuXs1/lyPAyxXdFK7w83CBVNq5b/fAQnMJFhoiIrInJpOA/PIanCzS4WTRH2WnSFt72V3HW+Iik2JEhC/uuKEbEqK7IcjbrWNDi4CF5hIsNERE5AiMJgFlegPUOgPUujqUVDX/WWceU+sMKKs24NJP8H5BCtwW3Q2339AN/YIUneJmnSw0l2ChISKizqTRaEJeWTW2nyzBtpNqHDxXYbFnJ1DpioToANwW3Q2xkX6QOznm/alYaC7BQkNERJ1Zmd6Andml2JapRmpOqcX8HA8XGUb3UWFQqDc6csrNxMHB6KZwteo2HbrQdO/eHXq9HsOHD8eOHTsuWz5jxgzs2bMHubm5rd4mCw0REXUVdQ1G7M0tw9aTamzLVKOkymCT1/3fUzfhxh4+Vt1maz+/7fKOWrm5uTh27BheeOGFy5YtXboUN954I7Zu3SpCMiIiIvvn6izDmKgAjIkKwD8n9sfxQi22nSzB+Yornz5uLb7uLh26/atpV6GprKxEVlYWevXqBX9//yuuc/78eWg0GvTp0wfu7m279LOHhwfc3C6fqX38+HFkZmbitddewzvvvNOe6ERERF2KVCrBwBBvDAzxFjtKh2rTVXtycnIwe/ZsREdHIzY2Flu2bLlsnZqaGkyYMAF9+vTBlClT0K1bN6xatcq8PD8/H56enld85Ofnt/jatbW1eOutt/Duu++2JTIRERF1AW0qNBkZGYiNjb3q3JX58+fj+PHjyMvLw+nTp/HJJ59gzpw5OHHiBAAgNDQUxcXFV3yEhoa2uN3169dj/fr18Pf3R3h4OEpKShASEtKW+ERERNRJtanQTJ48GXPmzGnxEJLRaMSqVavw5JNPIiAgAADwyCOPIDQ01LyXRiKRtLiH5mrny0+dOhVlZWUoLi7GiRMnoFKpkJ2d3eL6BoMBOp3O4kFERESdk1VvFJGbmwutVothw4ZZjA8fPhyHDx9u9XZuuukm3HTTTThw4AA8PT2xevVqODk5mYuPh4cHJBIJPDw8WtxGUlISlEql+XG1vT9ERETk2Kx6llN5eTkAwNfX12Lcz88Pp06davV2tm/fDqPxj/PnXV0tz2n39/fHmTNnrrqNxMREzJs3z/xcp9Ox1BAREXVSVi00zs7OAJoO9/xZbW2teVlrXOkMp0td68wpuVwOuVze6tckIiIix2XVQ05hYWEAgMLCQovxwsJC8zIiIiIia7NqofH398fAgQPx/fffm8cqKyuRmpqKhIQEa74UERERkVmbDjnpdDpkZmaan+fm5mLfvn0ICAhAREQEAOCtt97CxIkTER4ejpiYGLz33nsIDw/HI488Yt3kRERERBe16V5OBw4cwDPPPHPZ+MSJE5GYmGh+vnXrVixduhRlZWWIiYnB/PnzoVKprJO4nXgvJyIiIsfj0Den7AgsNERERI6ntZ/fVp1DQ0RERCQGFhoiIiJyeFa9Do09az6yxlsgEBEROY7mz+1rzZDpMoWmqqoKAHi1YCIiIgdUVVUFpVLZ4vIuMynYZDKhsLAQXl5eV70JZls131KhoKCgU0825vvsXPg+O4+u8B4Bvs/Opi3vUxAEVFVVISgoCFJpyzNlusweGqlUipCQkA7bvkKh6NT/8zXj++xc+D47j67wHgG+z86mte/zantmmnFSMBERETk8FhoiIiJyeCw010kul2PhwoWd/s7efJ+dC99n59EV3iPA99nZdMT77DKTgomIiKjz4h4aIiIicngsNEREROTwWGiIiIjI4XWZ69B0hIaGBvz222/QarUYOnQogoODxY5kdRs3bkRtba3FWL9+/TBgwACRElnP4cOHkZ2djVtvvRUBAQGXLTeZTEhPT4darUb//v3Rq1cvEVJeH0EQsGvXLqjVakyZMgVOTpb/5H/55ReUl5dbjPXs2RMjRoywZczrlp2djVOnTiEoKAgxMTFXvPhWZWUlfvvtNzg5OSEuLg6enp4iJG0/o9GIY8eO4cKFC4iMjER0dLTFcq1Wi59++umyr2vp/2971fw+CwoKEB4ejoEDB15xvZMnTyIrKwuhoaEYMmSIVS+Yagu1tbU4dOgQtFotoqKiEBERYbG8+efTn3l5eWH8+PG2jGk1WVlZOHLkCIYNG4bIyEiLZfX19fjtt99QVVWFYcOGITAwsF2vwULTTvn5+bj99tthNBrRo0cP7Nu3D++++y7++te/ih3NqubOnYuePXsiLCzMPCaVSh260Gzbtg2vvPIKysvLkZOTg507d172A1+r1WLs2LHIz8/HDTfcgD179uCZZ57B22+/LVLqtluxYgXeeecdCIKA3Nxc3H333Zd9iM+fPx8NDQ0WH44333yzwxSajIwMPP7446isrESvXr1w7NgxeHt7Y+PGjQgPDzev99NPP2Hq1Km44YYbUFdXhwsXLmDTpk0YOXKkeOHb4KeffsILL7wAd3d3BAcHY8+ePRg+fDjWr18Pd3d3AEBBQQEefPBBTJw4Ea6uruavHThwoMMUmt27d+PJJ5+Em5sbunfvjn379iEyMhI//PAD/Pz8ADSV9Llz5yIlJQWxsbE4fPgwBg8ejI0bN8LNzU3kd9A6ycnJWLhwIcLCwiCXy7F7927ce++9WLVqlfmXjtWrV+Obb77BqFGjzF/XvXt3hyw0FRUVuPvuu3HmzBksWbLEotCcOXMGt99+O2QyGYKDg7F//3589NFHeOyxx9r+QgK1y/jx44X4+Hihvr5eEARBSE5OFmQymZCVlSVyMuvq1q2bkJycLHYMq/ruu++EPXv2CAUFBQIAYefOnZet8/TTTwt9+vQRKisrBUEQhLS0NAGAsG3bNhunbb8VK1YIOTk5wnfffScAEKqqqi5bZ8iQIUJSUpII6azjt99+E9LT083P6+rqhNjYWGHcuHHmsaqqKsHPz09YsGCBeezRRx8VIiMjBaPRaNO87bVp0yYhNzfX/FytVguBgYFCYmKieSwjI0MAIBQVFYkR0Sp27NghnDt3zvxcq9UKISEhwksvvWQe++9//yvI5XIhIyNDEARBKCoqErp37y68/vrrNs/bXt999535Z4sgCEJWVpYgkUiEr7/+2jz23HPPCePHjxcjntVNnjxZeOONNwS5XC589tlnFssSEhKEhIQEoaGhQRAEQVi+fLng7OwsnDlzps2vwzk07VBWVoaffvoJzzzzDJydnQEA06ZNg0qlQkpKisjprC87OxsbNmzAoUOH0NjYKHac63bvvfciNja2xeWCIGDdunWYPXu2+XLb8fHxGD58ONauXWurmNdt9uzZrTpMlp+fj++++w779+9HXV2dDZJZz0033YThw4ebn8vlctx99904cuSIeWzLli2oqKjAc889Zx6bN28ecnNzsXfvXlvGbbd77rnH4pBEQEAA4uPjLd5ns927d+OHH35ATk6ODRNax5gxY9CjRw/zc4VCgcDAQPPNhQFg7dq1uOOOO9C/f38ATXstpk2b5lD/Nu+9916LS/mHh4dDLpdbvE+gaU/xpk2bsGvXLlRUVNg6plUsWbIEarUaiYmJly0rKirC9u3b8eyzz5r3TD366KPw8vLC119/3ebX4iGndsjMzITJZDL/gwKaDsP069cPGRkZIibrGF9//TUOHz6Mw4cPw9vbGykpKQ59yOlazp8/j8rKSovvLwAMGDDgih8gju7HH39Efn4+Tpw4AZPJhLVr11rs5nY027dvt/jeZWRkoFu3blCpVOaxfv36QSqVIiMjA3FxcWLEvC61tbXYs2cPpk2bZjEuk8nw3nvvQaFQYM+ePUhISMDatWvh5eUlUtK2a2xsxLfffou6ujrs2LEDtbW1eOmll8zLMzIy8NBDD1l8zYABA/DBBx+gtrbWYQ47lZSUYMeOHaiqqsJ///tfjBkzBg8++KDFOpmZmfjss8+gVquRk5OD999/H48//rhIidvu2LFjWLRoEdLT0yGTyS5bfvz4cQCw+Pfq5OSE6Ojodn2WstC0g1arBQD4+vpajPv5+aGsrEyMSB1m5cqVGDduHICmH6KTJ0/G1KlTkZGRcdW7njqyq31/KysrRUjUcf75z3/izjvvhEQiQWNjI+bMmYMHHngAp06dcqgPwWaLFy9GWloa0tLSzGNarfay76VUKoW3t7fDfj+feuopNDQ0YN68eeYxPz8/HDx4EIMGDQLQtOdt5MiRePnll/Hpp5+KFbXNjEYjNmzYAL1ejwMHDmD8+PHw8fExL7/S97N5fo1Wq3WYQqPRaLBhwwaUlZXhxIkTeOKJJ+Di4mJePmnSJLz11lvmOVKfffYZnnrqKdx4440YOnSoWLFbraamBlOnTsV7771nMZ/tz6z9s7ZzfiJ1sOZLNev1eotxvV5vMRmvM2guMwDg5uaG+fPnIzMzE6dPnxYxVcfqSt/fsWPHms8OcXJywsKFC1FcXIwDBw6InKztkpOT8eKLL2LNmjUWk5rlcvll30sAqK6udsjv50svvYTvvvsOmzdvRvfu3c3jgYGB5jIDAD169MBjjz2G77//XoyY7SaXy5GSkmI+bLZ//348++yzFsuv9G8TgEN9P2+44QakpKRg69at2L17Nz766CMsWbLEvPzmm282lxkAePLJJ9G9e3f8+OOPYsRts3fffRcNDQ1wdnZGSkoKUlJSYDKZ8Pvvv2Pz5s0ArP+zloWmHZpnaOfn51uMnzt37rJT7zqb5tu8l5aWipyk4/To0QNOTk78/jqQ5jlPK1euxNSpUy2WRUZGQq1Ww2AwmMeKi4thMBgc7vuZmJiIZcuW4eeff27Vb+kKhcLhvpd/plAoMHHiRIs9bpGRkVf8t+nj4wNvb28bJ7SO3r17Y+TIkRbv80q8vLwc5vsZFhaGIUOGYMOGDeaH0WjE0aNHsXXrVgAd8FlqjRnMXVHfvn2FJ554wvz8+PHjAgBhy5YtIqayrpKSEvPM82avvfaa4OrqajFD31Fd7SynO++8Uxg7dqz5eWlpqeDm5iYsXbrUhgmto6WznCoqKoTa2lqLsc8//1yQSCTC6dOnbRnxunz55ZeCi4tLi2fj5eXlCTKZTPj222/NYx999JHg7u4uaLVaW8W8bomJiYJCoRD27t17xeWFhYUWz00mkxAfHy+MGjXKFvGs4tL3IAiCMHbsWGH06NHm52+88YbQrVs3obq6WhAEQTAajcKQIUOE6dOn2yzn9TCZTEJxcbHFWE1NjRAaGirMmzfPPHbpf4tjx44JMplMWL16tU1ydoRLz3IymUxCeHi48Pzzz5vHfv/9dwGAsGvXrjZvn3No2umjjz7CPffcAxcXF0RGRmLx4sUYP3487rzzTrGjWU1GRgb+7//+D5MmTUJQUBB+++03rFu3Dh988IHFDH1Hc/bsWaSnp5svKLdjxw4UFxejf//+5slpb7/9NuLj4zFjxgzExsZixYoViI6OxqOPPipi8rY5dOgQTp06hd9//x0AsH79esjlcowePRpBQUG4cOEC/vKXv2Dy5MkIDw/H4cOHsXz5csyfP/+yC1/Zq19++QUzZszAhAkT4OTkZHGWYfOemrCwMPztb3/DY489htOnT6Ourg5vv/02kpKSzHuk7N0HH3yApKQkPP/888jLy0NeXh4AwMfHx/wzZ9myZThw4ABuu+02yOVyfP3118jMzMSWLVtETN42Tz31FLy9vTFixAgIgoAffvgBu3fvtngPzz77LNasWYM777wT06ZNwy+//IKzZ8/iq6++EjF56wmCgDFjxmDs2LHo168fdDodkpOTIZPJLOZEJSQkmM/mKiwsxMcff4yEhITLJg47MolEgo8++ghTpkyBVCpFjx498OGHH+K+++7DzTff3PbtCQLvtt1eBw8exJo1a6DVahEbG4tZs2aZT+PuLHJzc7F27Vrk5+ejR48eeOCBBy67Qqmj+fXXX/HZZ59dNj5lyhRMmTLF/PzUqVNYsWIF1Go1BgwYgCeffBIeHh62jHpdVq5ciZ9//vmy8b///e+IiYkB0HTa5OrVq5GTk4OgoCBMmDABw4YNs3XUdtu0aRO+/PLLKy679BIK3377LbZs2QInJydMmjTJoX75+Pe//409e/ZcNh4eHm5xscddu3bhhx9+QFVVFfr27YuZM2daTKi1dyaTCevXr8evv/4Ko9GIPn364OGHHzZP+m1WUVGBzz77DNnZ2QgJCcHcuXMtTve2dzU1NUhOTsbBgwfh5uaGQYMGYdq0aRbzRmpra7FmzRocPHgQSqUScXFxuPfee8ULbQUPP/wwHn74Ydx2220W4/v378fatWtRVVWFuLg4zJw584pnRV0LCw0RERE5PE4KJiIiIofHQkNEREQOj4WGiIiIHB4LDRERETk8FhoiIiJyeCw0RERE5PBYaIiIiMjhsdAQERGRw2OhISIiIofHQkNEREQOj4WGiIiIHB4LDRERETm8/wfRem62YJPiWQAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "k = 40\n",
        "U, sigma, VT = la.svd(A)\n",
        "pt.semilogy(sigma[:k])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 55,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(800, 40)"
            ]
          },
          "execution_count": 55,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Q = U[:, :k]\n",
        "Q.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "np.float64(8.928606558086162e-14)"
            ]
          },
          "execution_count": 36,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(Q@Q.T@A - A, 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 49,
      "metadata": {},
      "outputs": [],
      "source": [
        "def interp_decomp(A, k):\n",
        "    idx, proj = sli.interp_decomp(A, k)\n",
        "    P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]\n",
        "    return P, idx[:k]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 53,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "(800, 40)\n"
          ]
        }
      ],
      "source": [
        "P, J = interp_decomp(Q.T, k)\n",
        "print(P.T.shape)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 62,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "np.float64(2.3504495355553674e-15)"
            ]
          },
          "execution_count": 62,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(Q.T[:, J]@P - Q.T, 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1af2d697-74d0-43f4-bcd7-3c56bb635324",
      "metadata": {},
      "source": [
        "Recast this as $\\tilde P Q_J \\approx Q$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 66,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "np.float64(2.3504495355553666e-15)"
            ]
          },
          "execution_count": 66,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(P.T@Q[J] - Q, 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d7ff68b6-e601-4d0b-84c1-7a87b7d13bfd",
      "metadata": {},
      "source": [
        "Now try the same thing on $A$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 65,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "np.float64(1.5785357720000255e-13)"
            ]
          },
          "execution_count": 65,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(P.T@A[J] - A, 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "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.14.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}