penalty_matrix_for_psplines

penalty_matrix_for_psplines.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/avivajpeyi/b5329ca13d51c0d6b77465f96400d748/penalty_matrix_for_psplines.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Penalty matrix for P-splines"
      ],
      "metadata": {
        "id": "ymAGJYlClECF"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 651
        },
        "id": "lIYdn1woOS1n",
        "outputId": "a546efc8-721e-4bb5-eeb6-49e3af0c659a"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Shape: (8, 8)\n",
            "Frobenius norm ||P_quad - P_skfda|| = 0.2300470611495532\n",
            "Max abs diff = 0.11499252631071144\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1400x600 with 4 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAABSgAAAJoCAYAAABlbg9qAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdeVJREFUeJzt3XmcT3X///HnmRkzY5vBmLFkDFKESKgmlSWRS0pElIylHSXVVaqvtUzLdUlFxFVIaY9WhMuSUOJyIZKKKFljxjrDfN6/P/zmczmzMJ+ZM3M+c+Zxv93O7dbnfM55n9c5M3j1ei/HMsYYAQAAAAAAAIALQtwOAAAAAAAAAEDJRYESAAAAAAAAgGsoUAIAAAAAAABwDQVKAAAAAAAAAK6hQAkAAAAAAADANRQoAQAAAAAAALiGAiUAAAAAAAAA11CgBAAAAAAAAOAaCpQAAAAAAAAAXEOBEgCAYmbBggXq16+fLrzwQkVFRSkiIkLVqlXTddddpxdffFH79u1zO8RiZfv27bIsS7Vq1XI7lKCwYsUKtW/fXpUqVVJISIgsy9L06dPPeV7r1q1lWZZtCw8PV7Vq1dS5c2d99tlnhR98Ici8ryVLlrgdikaOHOl/trGxsTp58mSux/75558KCwvzH//WW28VYaTnFkzPFQAAuI8CJQAAxcT+/ft13XXXqX379po+fbpOnjypNm3aqFu3brrooou0YsUKDR06VHXq1NG3337rdrgohnbt2qVOnTpp4cKFatSokXr37q2kpCTVrVs3z200adJESUlJSkpK0o033qgyZcro888/14033qgHH3ywEKMvWpnFwpEjR7py/f379+vTTz/N9fsZM2YoIyPD8etSWAQAAIUhzO0AAADAuaWkpOiqq67Sli1bVL9+fU2ZMkVXX3217Zi0tDTNmDFDI0aM0J9//ulSpMXPeeedp82bN6tUqVJuh+K6r776SocOHdJtt92mt99+O19tdOnSxVa08/l8euqpp5ScnKyXX35ZN910k9q2betQxCVT8+bN9f333+uNN95Qt27dcjxm2rRpioiIUL169bR+/foijvDc3nzzTR07dkw1a9Z0OxQAABAEGEEJAEAxMHjwYG3ZskW1atXSN998k604KUkRERG6++67tW7dOl100UUuRFk8lSpVSvXr19f555/vdiiu27FjhyTpggsucKzNkJAQjRkzRnXq1JEkvf/++461XVI1adJEl156qebPn69du3Zl+/7rr7/WTz/9pC5duqhixYouRHhuNWvWVP369VWmTBm3QwEAAEGAAiUAAEHu119/1axZsyRJ48aNU6VKlc56fJUqVVSvXr1s+999911de+21qlSpkiIiIpSQkKD+/fvrp59+yrGdWrVqybIsbd++XXPnzlXr1q0VHR2tihUr6oYbbtCGDRv8x86aNUuJiYkqX768KlSooK5du+qXX37J1uaSJUtkWZZat26tY8eO6YknnlDdunUVGRmp6tWra8CAAfrjjz9yjGfhwoUaPHiwLrnkElWuXFkRERGqUaOGbr31Vq1evTrHc86chrtjxw4NGDBA8fHxKlWqlPr27Svp7GtQbt26Vf3791ft2rUVERGhcuXKKSEhQZ06ddK0adNyvOb8+fN1ww03KC4uTuHh4apevbpuvfVWff/99zkef+aU2XXr1qlr167++2vQoIH++c9/yhiT47nnktef+fTp02VZlkaMGCFJGjVqlH/tQifW5gwNDdUll1wi6fTzPtOuXbs0dOhQXXTRRSpTpozKly+vFi1aaMKECTp16lS2tvr27etfF3Pbtm264447VLVqVUVEROj888/XU089pbS0tGznHT58WFOnTlXXrl11wQUXqGzZsipbtqwuvvhiPfnkkzp06FCe78eyLI0aNUqS/VlZlqW+ffsqNTVVUVFRCgsL086dO3Nt529/+5ssy9Krr76a52tn6t+/vzIyMjRjxoxs373xxhv+Y3IT6PPI/LO7dOlSSVKbNm1s9525TumZf54yMjI0btw4NW3aVOXKlZNlWf72zjVV/N///re6d++uGjVqKCIiQrGxsWrRooVGjBihAwcOZDv+p59+0j333KPzzz9fkZGRio6O1jXXXJPr2pspKSl66qmndPHFF6ts2bKKiIhQ9erV1bJlSw0fPvys63sCAIBCYAAAQFB76aWXjCRToUIFc+rUqYDP9/l8pk+fPkaSCQsLM23btjU9e/Y0F154oZFkypQpY+bOnZvtvISEBCPJPP7448ayLNOyZUvTo0cP/3kVKlQwP//8s3n00Uf97d5yyy0mPj7eSDLVq1c3f/31l63NxYsXG0kmMTHRXHHFFaZMmTLmb3/7m+nevbupVq2akWSqVq1qfvrpp2zxnH/++SY8PNw0bdrU3HjjjaZr166mQYMG/vv68MMPs50zYsQII8ncdtttplKlSqZq1aqmW7dupmvXrubhhx82xhizbds2I8kkJCTYzt2wYYOJiooykky9evVM165dTffu3U1iYqIpV66cadKkSbbrPfXUU0aS/3n16tXLXHLJJUaSCQ0NNa+//nq2c1q1auV/zuHh4eaiiy4yPXv2NK1atTKhoaFGknnwwQfP8hPOLtCf+ddff22SkpJMkyZNjCTTpEkTk5SUZJKSkvzP6Vwy72PEiBE5ft+uXTsjydx4443+fUuXLjUVK1Y0kkytWrXMjTfeaDp06ODf1759e5Oenm5rJykpyf9MoqKiTEJCgunRo4dp166dKV26tJFkunTpku36X3/9tZFkYmNjzVVXXWVuvfVW0759exMTE2Mkmbp165r9+/fnel+LFy+2xZDTs0pKSjJTp041xhgzePBgI8k88cQTOT6Pn3/+2ViWZaKioszhw4fP9XiNMf/7fR4wYID566+/TGRkpLngggtsx6SmppqyZcuamjVrmoyMDH/8M2fOLNDz2Lx5s0lKSjJVqlQxkkyHDh1s9/31118bY/7356lmzZrmxhtvNOHh4ebaa681vXr1Mo0bNz7rc82U+ewkmUsuucT07NnTdOzY0dSpUyfHc95//30TGRlpJJn69eubm2++2bRt29aULVvWSDL9+vWzHX/06FHTqFEj//137tzZ9OzZ07Ru3dpUrVrVSDIHDx7M088EAAA4gwIlAABB7o477jCSTNu2bfN1/qRJk4wkU7lyZfOf//zHv9/n8/kLHhUqVDB79+61nZdZoIyIiDALFy707z916pTp3r27kWQaNWpkYmJizLp16/zfHz161Fx55ZVGknn66adtbWYWKDMLIL/99pv/u+PHj5tu3boZSeaKK67Idh+zZ8/OVvDM3B8WFmZiYmLMsWPHbN9l3p8k07t3b3PixIls5+dWoOzXr1+O92CMMceOHTNLly617Zs7d66RZCIjI81XX31l++5f//qXkWRKlSplNm7caPsus1AjyUyePNn23aJFi4xlWSY0NNTs3LkzWxy5ye/PPPO73IqMZ3O2AuWuXbtM+fLljSQzfPhwY4wxf/75p4mJiTGWZZlXX33VZGRk+I/fv3+/adu2rZFkRo0aZWsrs0ApyTz55JO2ov2GDRv8RakVK1bYztu5c6dZuHCh7TrGnP59zSzm3n///bneV9ai2Lme1U8//WQsyzJxcXE5/t49/PDDRpIZPHhwjufn5MwCpTHG9OrVy0gyy5Yt8x8zdepU23POrUDp9PPIlPnnSZKpUaOG2bJlS47H5dbOyy+/bCSZmJgY8+9//zvbed9++63ZsWOH//P69etNRESEiYyMNB999JHt2O3bt5uLL77YSDIzZszw758xY4aRZDp27JitAJ6RkWGWLFli0tLScowbAAAUDqZ4AwAQ5Pbt2ydJiouLy9f5//jHPyRJw4cP90+zleSf0tu4cWMdOnRIU6dOzfH8Bx54QNdee63/c2hoqIYNGyZJ2rhxo0aPHq0mTZr4vy9TpowefvhhSdKiRYvOGteZL8iIjIzUq6++qjJlymjVqlVasWKF7fjc1tPr0qWLunfvrgMHDmjx4sU5XqtSpUqaMGGCIiIico0nqz179kg6PQ03q9KlS+uaa67Jdj+SdP/99+u6666zfTdgwADdcMMNOnnypF566aUcr9e1a1fdc889tn1t27ZVhw4dlJGRkeu95aSgP3OnHD16VMuWLdMNN9ygw4cPq2zZsrrzzjslSePHj9eBAwc0cOBA3XfffQoJ+V9aGhMTozfffFOlSpXShAkTcpzi3qxZM40ZM0ahoaH+fY0aNdIdd9wh6fSSAGeqUaOGrr32Wtt1pNO/r5MmTVJYWJg++OADx+79ggsuUMeOHbV3795s7R4/flxvvPGGLMvSwIED832NzCncmVO6Jen111+XZVnq16/fWc8tiucxduxYXXjhhXk+/tSpUxozZowkacqUKWrTpk22Yy677DLFx8f7Pz/zzDNKS0vT008/ra5du9qOTUhI0Ouvvy5Jevnll/37M/9sX3fdddlejhUSEqJWrVopPDw8z3EDAICCo0AJAICH/f777/61IJOSkrJ9f2YhI7cCWE4FujNfonK273N6gYckVahQQTfeeGO2/XFxcbr++uslKce16Xbt2qWpU6fq4Ycf1p133qm+ffuqb9+++uGHHyRJW7ZsyfF67dq1U3R0dI7f5eayyy6TJN13332aP3++Tpw4keuxp06d0jfffCNJ/rUtsxowYICk3J9z586dc9yf+cKj3NbmzMqJn3lBnLkmY7ly5dSqVSutXbtWcXFx+uSTT/zFpS+++EKSdOutt+bYznnnnacLLrhA+/bt09atW7N9f8MNN9jWNMx0rue1YsUKPffccxo4cKD69eunvn376v7771d4eLj27dungwcP5uu+c/Lggw9KkiZMmGDbP2vWLB08eFDt2rXLcb3YvLr22muVkJCgDz74QEeOHNHmzZu1atUqtWnTJs9rhxbm88jtDeO5WbNmjfbt26fKlSvr5ptvPufxPp9Pc+fOlZT771Hz5s1Vrlw5/ec///H/GW7RooUk6fnnn9ebb76pv/76K6A4AQCA88LcDgAAAJxdbGysJGnv3r0Bn5tZpImJiVFUVFSOx2S+vTq3gs6ZoxwzlStX7qzfly9fXpJyLeplvoAnJ7Vr15Z0utB2plGjRumZZ54568srUlNTc71eoB599FEtX75cCxcu1PXXX69SpUqpSZMmuuaaa9SzZ09/kUOSDhw44L/XzPizys9zluT/uZ2tQHomJ37mBdGkSRP/qM1SpUqpUqVKatasmTp37qzSpUv7j/v1118lKcc30me1b9++bCPxAn1ee/fuVbdu3bR8+fKzXis1NdWxN19fd911uuiii/Ttt99qzZo1atasmSRp4sSJkqRBgwYVqP3Ml/KMGjVK7733nn788UdJZ385TqbCfh5xcXEBv6H7t99+kyTVq1cv178fznTgwAH/n/kzR1We7fjzzjtPrVu31mOPPaYXXnhBSUlJsixLF1xwgVq2bKmbbrpJnTt3zjayFAAAFC4KlAAABLlmzZpp5syZWrt2rTIyMmxTWovCuf5HvbD+R/7Mab0ff/yxRo4cqXLlymnChAlq27atqlevrtKlS8uyLD3xxBNKTk7O9W3XZxbG8qpMmTJasGCBVq9erXnz5mnFihVasWKFvv/+e40bN07333+/v9DkBK8URLp06aKRI0ee8zifzydJuuWWW1S2bNmzHhsTE5NtX6DP684779Ty5cuVmJioUaNGqUmTJqpYsaJ/im/16tX1559/5vuN6TmxLEuDBw/W/fffrwkTJmjatGlauXKl/vOf/6hWrVq64YYbCnyNfv36afTo0ZoyZYp+++03RUdHZ5vqnJPCfh75+TMXqMzfISnn0cJZnbnEw7PPPqt7771Xn332mZYvX65vvvlG06ZN07Rp09SiRQstXrz4nL+XAADAORQoAQAIcjfccIOGDh2qQ4cO6dNPP83T1MdM5513nqT/jTTKaURd5ki2zGOLwvbt28/5XY0aNfz73n//fUmn15u7++67s52T0xRgp7Ro0cI/WvLUqVOaM2eO+vTpo1dffVW33HKL2rRpo5iYGEVERCgtLU2//vqrGjdunK2donrOwfozzyo+Pl5bt27VY489pubNmxfqtY4ePaovv/xSISEh+vLLL1WhQoVs3+/evbtQrt2nTx898cQTevfdd/WPf/zDP90767qb+ZWQkKC2bdv613u99957z1kcdPN5nE3mqNiffvpJxphzjqKsXLmySpcurePHj+sf//iHKleuHND1atWqpcGDB2vw4MGSpNWrV6t3795avXq1nn/+eY0aNSp/NwIAAALmja56AAA87Pzzz1evXr0kSQ8//PA510vbu3evfy3GGjVq+KfzTp8+Pduxxhj//pxeSFFYDh06pM8++yzb/n379mnevHmSpNatW/v3Z95zQkJCtnP27t2rBQsWFE6gWYSFhemWW25Rhw4dJEnr1q3z77/qqqsk5fycpf+9yKSwn3Ow/syz6tixo6T/FZ8LU0pKijIyMhQVFZWtGCdJb731VsAjBTNfonLq1KmzHle2bFkNGDBAJ06c0NixY/Xhhx8qMjLSvyapE+6++27FxMQoJiYmT+0W5Hnk9b7zo3nz5qpcubL27dunOXPmnPP40NBQ/wupnPg9atGihe6//35J//uzDQAAigYFSgAAioFXXnlFdevW1bZt23TVVVfluG5cenq63njjDTVt2lSbN2/273/kkUckSWPGjNF///tf/35jjJ5++mmtW7dOFSpU0F133VX4N3KGhx9+2LbOZFpamgYOHKijR4/qsssuU8uWLf3fZb74ZMqUKUpPT/fvT0lJUVJSklJSUhyP79VXX83xpTu7d+/W999/L8leMM18c/mkSZOyvb18+vTp+vTTT1WqVCn/i1MKU7D+zM/06KOPqkKFCho3bpz++c9/2n6umbZt26a33nqrwNeqUqWKKlasqEOHDmnmzJm271atWuV/K30gMkf4Zr6g6WwGDRqkkJAQjRs3Tunp6erVq1eO09bzq0ePHtq/f7/279+fp9GoBXkegdx3oMLCwvTkk09KOl10XbZsWbZjVq9ebft7Y8SIEQoPD9ejjz6qGTNm2KZ9Z9q4caM+/vhj/+fZs2dr2bJl2Y49efKkv4Mkp84QAABQeJjiDQBAMVCxYkV98803uvXWW7VkyRJdffXVql27tho3bqwyZcpoz549+u6773TkyBFFRUWpevXq/nPvuecerVixQjNnzlTz5s3VqlUrxcXFae3atdqyZYtKly6tWbNm+V/GUxQSExPl8/lUr149tW3bVmXKlNHy5cu1a9cuxcXF6c0337QdP2TIEL355pv68ssvVadOHV1xxRU6efKkli5dqjJlyqh///7+EYpOmTJligYOHKjatWurUaNGioqK0r59+/T111/r+PHjatu2re1N5B07dtRTTz2lp59+Wtddd51atmypmjVr6scff9TatWsVGhqqyZMnq2HDho7GmZNg/JlnVaNGDX3yySfq1q2bHnnkET3//PNq1KiRqlWrppSUFG3evFm//PKLLr/8cvXu3btA1woNDdXw4cP10EMPqU+fPpo4caLq1KmjHTt2aMWKFerdu7eWLVvmf0lLXnTo0EFly5bVnDlzdNVVV+mCCy5QaGioWrZs6X9LeqZatWrpxhtv9I8KLOjLcQqqIM+jW7dumjZtmv7+979r4cKFiouLk2VZ6t+/v6688soCx/bggw9qy5Ytmjx5slq1aqWmTZuqXr16Sk1N1Y8//qhff/1Vixcv9hdKL730Ur311lvq27ev+vbtq6eeekoNGjRQbGys/vrrL23YsEG///67br31Vv/anEuXLtVLL72kypUrq2nTpoqLi9Phw4e1atUq7d27V+edd57+/ve/F/heAABA3lGgBACgmIiLi9PixYs1b948vfPOO1qxYoUWLVqktLQ0xcTEKDExUZ06ddIdd9yhSpUq+c+zLEtvvvmmOnbsqClTpmjNmjU6evSoqlatqr59++rxxx9XvXr1ivRewsPD9cUXX2jUqFH68MMP9ccff6hixYrq27evRo8ene2NvLVr19Z//vMfPfXUU/r666/1+eefq2rVqurVq5dGjhypSZMmOR7jM888oy+++EKrVq3SqlWrlJKSori4OF1++eXq16+fevXqpbAweyo1ZswYtWzZUq+88oq+/fZbrVq1SpUrV1b37t31yCOP6LLLLnM8zpwE4888J9dcc41++OEHTZgwQV988YVWr16ttLQ0xcXFqWbNmurdu7e6devmyLWGDBmi2rVr6/nnn9emTZv0ww8/qH79+po4caLuvffeXN++npsqVapo7ty5Gj16tNasWaOVK1fK5/Pp1KlT2QqU0umC5pw5c5SYmKhLL73UkXsqiPw+j06dOmnq1KmaNGmS/v3vf+vYsWOSpKuuusqRAqVlWZo0aZJuuukmTZ48WatWrdLGjRtVoUIF1a5dW0lJSdnWeO3evbtatGihl19+WQsWLNA333yjjIwMValSRXXr1tWgQYN0yy23+I/v27evSpcureXLl2vTpk1aunSpoqOjVbNmTQ0ZMsQ/ZR4AABQdyzj5qkIAAICzWLJkidq0aaNWrVppyZIlbocDFJmrrrpK33zzjWbNmuVfUxYAAACnsQYlAAAAUIjmzp2rb775RjVr1rSN5AMAAMBpTPEGAAAAHHbgwAE99thjOnjwoL788ktJ0vPPP69SpUq5HBkAAEDwoUAJAAAAOOzw4cN6/fXXFRYWpjp16ujhhx/Wrbfe6nZYAAAAQYk1KAEAAAAAAAC4hjUoAQAAAAAAALiGAiUAAAAAAAAA11CgBAAAAAAAAOAaCpQAAAAAAAAAXEOBEgAAAAAAAIBrKFACAAAAAAAAcA0FSgAAAAAAAACuoUAJAAAAAAAAwDUUKAEAAAAAAAC4hgIlAAAAAAAAANdQoAQAAAAAAADgGgqUAAAAAAAAAFxDgRIAAAAAAACAayhQAgAAAAAAAHANBUoAAAAAAAAArqFACQAAAAAAAMA1FCgBAAAAAAAAuIYCJQAAAAAAAADXUKAEAAAAAAAA4BoKlAAAAAAAAABcQ4ESAAAAAAAAgGsoUAIAAAAAAABwDQVKAAAAAAAAAK6hQAkAAAAAAADANRQo4VkjR46UZVluh4Espk+fLsuytH37drdDcUVh3f/zzz+v+vXry+fzOdpuYZk8ebJq1qyptLQ0t0MBAAAAALiMAiUc8cMPP6h3794677zzFBERoerVq6t3797atGmT26E5bsWKFRo5cqQOHTrkdijZZBa/MrfIyEhdeOGFGjRokPbs2eN2eLkqrGd65vNYvnx5tu+NMYqPj5dlWbrhhhsCbj9YfhdSU1P13HPP6bHHHlNIiP2v9fnz59t+J0JDQ1WrVi099NBDOnLkiEsRS3379lV6erpee+0112IAAKAoZHaa79+//5zHrl69WldeeaXKli0ry7K0bt26c7ZbktDRTUf3mejwBryFAiUK7OOPP9all16qRYsWqV+/fnr11Vc1YMAA/fvf/9all16qTz75xO0QHbVixQqNGjXK9aLU2YwePVozZ87UhAkTdOWVV2rSpElKTEzUsWPH3A4tR4X9TCMjIzVr1qxs+5cuXarff/9dERER+Wo3P3HfcccdOn78uBISEvJ1zZy88cYbOnXqlHr16pXtu//+97+SpBdffFEzZ87Ua6+9poYNG2r8+PF6+OGHHYshUJGRkUpKStK4ceNkjHEtDgAAgsXJkyfVvXt3/fXXX/5/t53MF86lOHZ008ldcLl1dAdrJ/eZ6PAGvIUCJQrkl19+0R133KE6depo/fr1evrppzVgwACNGTNG69evV+3atdW7d29t27bN7VBzdfToUbdDkORsHB07dlTv3r115513avr06RoyZIi2bdvmuWJxXv3tb3/TBx98oFOnTtn2z5o1S82aNVPVqlWLJI6jR48qNDRUkZGRjo54mDZtmm688UZFRkZm+279+vUqW7asHnjgAf/vxKeffqqEhATXfx969Oih3377TYsXL3Y1DgAAgsEvv/yi3377TY888ojuvvtu9e7dWxUrVizyOIpTR3dJ6uSWirajO1g7uc9EhzfgLRQoUSAvvPCCjh07pilTpig2Ntb2XeXKlfXaa6/pyJEjeuGFF/z7+/btq1q1amVrK+s0ld9++03333+/6tWrp9KlSysmJkbdu3fPcUrD8uXL1aJFC0VGRur888/PtRct8xqbNm3SbbfdpooVK+qqq67K8/VGjhypRx99VJJUu3Ztf4/i9u3b83xf54pDkv744w/1799fVapUUUREhBo2bKg33ngjx3vKi7Zt20qSrVCc12tkxvrzzz+rb9++qlChgqKjo9WvXz9bohrIzytr+zk902nTpsmyLM2ePTvbObNmzZJlWVq5cmWe7r9Xr146cOCAFixY4N+Xnp6uDz/8ULfddlu24wv6u3Dmc8v6M846Nef48eOqX7++6tevr+PHj/vb/+uvv1StWjVdeeWVysjIyPXetm3bpvXr16tdu3Y5fv/f//5XjRs3tvWIh4aGKi4uTocPHz7nsytMzZo1U6VKlVwvlAIAEAz27t0rSapQoYKrcdDR/T/B1MktqUg7uoO5k/tMdHgD3kGBEgXy2WefqVatWrr66qtz/P6aa65RrVq19NlnnwXc9urVq7VixQr17NlTL7/8su69914tWrRIrVu3thXGNmzYoPbt22vv3r0aOXKk+vXrpxEjRuRY2MrUvXt3HTt2TGPHjtVdd92V5+t17drV37uY2Zs4c+bMbMXZvMopjj179uiKK67QwoULNWjQIL300kuqW7euBgwYoPHjx+frOr/88oskKSYmJt/X6NGjhw4fPqzk5GT16NFD06dP16hRo/zf5/XnlVVuz7RHjx6Kj4/X22+/ne2ct99+W+eff74SExPzdP+1atVSYmKi3nnnHf++uXPnKiUlRT179sx2vJO/Czn9jM9UunRpzZgxQz///LOefPJJ//6BAwcqJSVF06dPV2hoaK73tmLFCknSpZdemu279PR0bdmyRU2bNrXt37Nnj3744Ycczylql156qb755hu3wwAAIN8OHz6sIUOGqFatWoqIiFBcXJyuu+46rV27NtdzfvvtN9WtW1eNGjXSnj171LdvX7Vq1UrS6dzBsiy1bt3af3xeOuPz21l8Lvnt6Hark3v79u1avHixIx3dwdTJLdnXoCxoJ7d09o7uYO7kPhMd3oB3hLkdAIqvlJQU7dq1SzfddNNZj2vcuLE+/fRTHT58WOXLl89z+506ddItt9xi29e5c2clJibqo48+0h133CFJGj58uIwx+vrrr1WzZk1JUrdu3XTxxRfn2naTJk2yTdfIy/UaN26sSy+9VO+88466dOmS44jJQOQUx5NPPqmMjAxt2LDBX1C899571atXL40cOVL33HOPSpcufdZ2U1JStH//fp04cULffPONRo8erdKlS/vXyMnPNZo2barXX3/d//nAgQN6/fXX9dxzz0nK+88rq7M90969e2vcuHFKSUlRdHS0JGnfvn366quvbMW8vLjttts0bNgwHT9+XKVLl9bbb7+tVq1aqXr16tmOdfJ3IevPePr06dmOufzyy/X3v/9dzz33nG6++Wbt2bNH7777rsaPH68LL7zwrPf1448/Sjqd4Ga1adMmnTx5Uueff77279+v9PR0/fDDD3riiSeUlpamESNGnLXtolCnTh3NnDnT7TAAAMi3e++9Vx9++KEGDRqkBg0a6MCBA1q+fLk2b96cY2fgL7/8orZt26pSpUpasGCBKleurHvuuUfnnXeexo4dqwceeEAtWrRQlSpVJP2vMz42NlYjR47UqVOnNGLECP/3mc7sYK1Ro4a2b9+uSZMmqXXr1tq0aZPKlCmTr/vLraPbsiwNGjRIsbGxmjt3rgYMGKDU1FQNGTLEdn6PHj1Uu3ZtJScna+3atfrXv/6luLg4fw6Z37i7du2qn376Se+8845efPFFVa5cWZIUGxurhIQEf0f3zTffbDsvkI7uMzu5O3bsKMneyf3yyy/bjs/LvZwt7jN1795dF1xwgcaOHZvj9OXMTu6WLVvqySef1Lhx4yTlvZNbyr2jO7OTO2vnejB1cp+JDm/AIwyQTzt37jSSTO/evc963O23324kmT/++MMYY0xSUpJJSEjIdtyIESNMbr+S6enpZv/+/Wbfvn2mQoUKZsiQIcYYY06dOmVKly5tevbsme2cv/3tb9nay7zG0qVLzxpzbtczxpgXXnjBSDLbtm2znRPIfeUWh8/nMxUqVDB333232bdvn22bNm2akWSWL1+ea9yZx2TdEhISzLx58/J1jcxYv/vuO9u1xo0bZySZlJSUgJ5f5jXOfH65PdPNmzcbSeZf//qXf98rr7xiJJmtW7fm+hyyXmv16tVm7969JiwszLz//vsmNTXVlC5d2kydOtUYY0xCQoLp1KlTjm3k53fBmNx/xjndvzHGpKWlmYsvvtjUrl3bxMbGmlatWhmfz3fOe7zvvvtMWFhYjt/NmDEjx9+HevXqmS+++OKcbTutVatWZubMmbZ9jz32mJFkjh49WuTxAADghOjoaDNw4MBcv8/MCfbt22c2b95sqlevblq0aGH++usv23GLFy82kswHH3xg29+lSxcTGRlpfvvtN/++TZs2mdDQUFuOeezYsWzXXrlypZFk3nzzzXPeR2aOsnDhQrNv3z6zc+dO8+6775qYmBhTunRp8/vvvxtjjBkwYICpVq2a2b9/v+38nj17mujoaH8cmffdv39/23E333yziYmJyVfcWfOos+Viw4YNMxEREebQoUP+fZn54IgRI/L0LFavXm0mTJhgypcv74+ze/fupk2bNsaY7DlkXu8lLzlkr169co3rzPOGDRtmQkJCzLJly8wHH3xgJJnx48ef9f4yPfXUU0aSOXz4sG3/f/7zHyPJ/POf/zT79u0zf/zxh/nqq69M8+bNTWhoqFmwYEGe2i8qd999tyldurTbYQAoIKZ4I98yR0Oea4j/4cOHZVmWv3cwr44fP67hw4crPj5eERERqly5smJjY3Xo0CGlpKRIOj2a7vjx47rggguynV+vXr1c285ptFlerue0rHHs27dPhw4d8q/peebWr18/Sf9bn+hsJk6cqAULFmjx4sXatGmTfv31V3Xo0KFA18gcnZopc9H2gwcPSiqc51e/fn21aNHCNs377bff1hVXXKG6desG1FZsbKzatWunWbNm6eOPP1ZGRka2UZKZnLyXnH7XchIeHq433nhD27Zt0+HDh/1rcBbEf//7X4WFhemrr77SggULtHTpUv3666/68ccf9be//a1AbTvF/P8RAU6upQQAQFGqUKGCvv32W+3ateusx23cuFGtWrVSrVq1tHDhwjy9ACcjI0Pz589Xly5dbLnYRRdd5M/tMp05++XkyZM6cOCA6tatqwoVKpx1unlW7dq1U2xsrOLj49WzZ0+VK1dOs2fP1nnnnSdjjD766CN17txZxhjt37/fv3Xo0EEpKSnZrnXvvffaPl999dU6cOCAUlNTHY07qz59+igtLU0ffvihf997772nU6dOqXfv3nlup0ePHjp+/Lg+//xzHT58WJ9//nmO07udvpeszy03I0eOVMOGDZWUlKT7779frVq10gMPPJCncw8cOKCwsDCVK1fOtn/9+vWSpIcfflixsbE677zz1L59ex0+fFiffvpprmufZ7IsS7///nuO323evFkXX3yxypcvr48++ijb961bt9Zbb72Vp/gzVaxYUcePHw/KFzkByDumeCPfoqOjVb16df8/YLlZv369atSoofDwcEm5FyKyrpEyePBgTZs2TUOGDFFiYqKio6NlWZZ69uwpn89XoNhzmiJd0Ovl9b7OFkfmdXr37q2kpKQcz2ncuPE5Y7nsssvUvHnzHL/L7zVymyKSWWAqrJ9Xnz599OCDD+r3339XWlqaVq1apQkTJuSrrdtuu0133XWXdu/erY4dO+a6CL2T93Ku6fhnmj9/viTpxIkT2rp1a56KmzExMTp16lSOSyisX79edevW1XXXXRdQzEXp4MGDKlOmTEDPCQCAYPL8888rKSlJ8fHxatasmf72t7+pT58+qlOnju24zp07q0qVKpo/f362glBuztUZ/+WXX/o/Hz9+XMnJyZo2bZr++OMP27TgQDpYJ06cqAsvvFBhYWGqUqWK6tWr51+H8MyO7ilTpuR4ftaO7rN1ckdFRTkWd1ZndnQPGDBAUv46us/s5D527Ng5O7mdupdAO7kz1yh1spP7yy+/lGVZCg8PV3x8fJ5jOpsXXnhBnTt31tixYwvcViY6vAFvoECJAuncubNee+01LV++3PYW6kxff/21tm/frqFDh/r3VaxYUYcOHcp27G+//Wb7/OGHHyopKUn//Oc//ftOnDhhOzc2NlalS5fW1q1bs7W3ZcuWgO4lL9eTcv+HL6/3dTaxsbEqX768MjIyztkzmV+FdY28Pr+cnC2Z6Nmzp4YOHap33nlHx48fV6lSpXTrrbfmK8abb75Z99xzj1atWqX33nsv1+MK+ruQH+vXr9fo0aPVr18/rVu3Tnfeeac2bNjgX3szN/Xr15d0epHzrIXl9evX5/oCqzOtXr1aAwYM0Pbt23Xbbbdp06ZNuvvuu9W7d29ZlqWdO3eqRo0akk73at95553q3bu3xo4dq9dee00HDx5Uw4YN9dprr9liyNpuTsXdbdu26aKLLjpnjAAABKsePXro6quv1uzZs/XVV1/phRde0HPPPaePP/7Yv26hdHqN9BkzZujtt9/WPffc43gcTnWwOt3R7VYnt+RcR7cXO7ml3Du6C7OTe8eOHf4XLzmFDm/AG5jijQJ55JFHVKZMGd1zzz06cOCA7bu//vpL9957r6KiojRo0CD//vPPP18pKSm2kZd//vlntrfshYaGZlsQ+pVXXrGNSAwNDVWHDh00Z84c7dixw79/8+bN/n+o8yov15OksmXLSlK2YlVe7+tcMXTr1k0fffSRNm7cmO37ffv25bmtor5GXp9fTnJ7ppJUuXJldezYUW+99ZbefvttXX/99QEvF5CpXLlymjRpkkaOHKnOnTvnelxBfxcCdfLkSfXt21fVq1fXSy+9pOnTp2vPnj166KGHznlu5gLv33//vW3/7t27tXfvXjVo0OCs56enp6tr16667777dODAATVs2NC/YPq51K9fX99//70OHDig6667Tn369Am43bVr1+rKK6/M0/UAAAhW1apV0/333685c+Zo27ZtiomJ0TPPPGM75oUXXtCAAQN0//33Z3tJYm4C6Yw/s4P1lltu0XXXXaerrrqqwHlK1njO7OjOaYuLiwuozYLEfa7O4p49eyo0NFTvvPOO3n777Xx3dN98880KCQnRqlWrcp3eLeX9Xgqrk7tp06a688478zxa88yO7qxtNmzY8Kzn+nw+PfDAA6pcubIqVKigFi1aaP/+/dmOmz17ts4//3xt3bpVHTt21OLFi3XnnXeqXLlyOnDggFavXq3GjRsrKipK9957r62QO3bsWCUkJCgqKkqJiYm5ztyjwxvwBkZQokDq1q2rN998U7169dLFF1+sAQMGqHbt2tq+fbtef/11HTx4UO+++66tF69nz5567LHHdPPNN+uBBx7QsWPHNGnSJF144YW2tVluuOEGzZw5U9HR0WrQoIFWrlyphQsX+t8gmGnUqFGaN2+err76at1///06deqUXnnlFTVs2PCc08/PlNfrNWvWTNLpN2H37NlTpUqVUufOnfN8X+fy7LPPavHixbr88st11113qUGDBvrrr7+0du1aLVy4UH/99Vee2yrKa+T1+eUkt2eaWQDs06ePfyrNmDFjAo7tTLn19p+poL8LmXHn1dNPP61169Zp0aJFKl++vBo3bqzhw4frqaee0i233HLW9SLr1KmjRo0aaeHCherfv79//3//+19JOmdyuXLlSoWFhem+++6TJA0aNEjPP/98nuLu2rWr/7+feOIJPf300zpy5IjKlSuXp3bXrFmjv/76SzfddFOergcAQLDJyMjQkSNHbDMe4uLiVL16daWlpdmOtSxLU6ZM0eHDh5WUlKRy5crpxhtvPGv7WTvjM6dL59QZX5DO4rzK7OieNWuWNm7cqEaNGtm+37dvX7a3UeelzcLo5JbsHd0nTpzId0d3Zif39u3bg7qTe9u2bWrRooUeeughvfHGG+c8/8yO7syRr3nt5P7qq6+0YsUK/frrrypbtqz++9//KjIy0nbMu+++q//7v//TwoULVbt2bc2dO9c2GyezQ/uJJ57QnXfeqcmTJ+tf//qX7r77bkn/6wyvUKGCxowZoz59+mjdunXZYlm7dq1uv/32vDwyAEGMAiUKrFu3blq7dq2Sk5P1r3/9S3v37pXP51NkZKTWrFmT7R+3mJgYzZ49W0OHDtXf//531a5dW8nJydq6dautkPfSSy8pNDRUb7/9tk6cOKGWLVtq4cKF2RYEb9y4sebPn6+hQ4dq+PDhqlGjhkaNGqU///wzoAJlXq/XokULjRkzRpMnT9a8efPk8/m0bds21apVK0/3dS5VqlTRd999p9GjR+vjjz/Wq6++qpiYGDVs2FDPPfdcntsp6mvk9fnlJLdnmpm8de7cWRUrVpTP5ztnIu+Egv4uBFKgXLt2rcaOHatBgwapTZs2/v2PP/64PvnkE91111364Ycfcp1KJEn9+/fX8OHDdfz4cf/Ulszf/XMll3/++ad/+rZ0+n+ezvx8NlOnTtX48eP1+++/y7IsGWN04MABlStXLk/tfvDBB6pZs6bj03wAACgqhw8fVo0aNXTLLbeoSZMmKleunBYuXKjVq1fblorJFBISorfeektdunRRjx499OWXX57z38G8dsYXpLM4EE53dBdmJ7fkXEe31zq5pZw7uvPayV2qVCkdPnxYP/74o1q0aKFLL73U9v2bb76p6dOna+HChUpISMixjXN1aJ+tMzwTHd6AhxT5e8NRIsyYMcNYlmXuuOMOt0OBB5w8edLExsaa/v37ux1KUDp06JCpVKmS+de//hXwuUuWLDG1atWy7atRo4aZOXOmMcaYMmXKmJ9++sn/Xf369c3MmTPNtm3bTJkyZczq1atNRkaGOXbsmLEsy2zbti1P7Z44ccJUrVrVjB8/PuCYAQAIFmlpaebRRx81TZo0MeXLlzdly5Y1TZo0Ma+++qr/mBEjRhhJZt++ff59x44dM61atTLlypUzq1atMsYYs3jxYiPJfPDBB9mus3TpUtOsWTMTHh5u6tSpYyZPnuxvN9PBgwdNv379TOXKlU25cuVMhw4dzI8//mgSEhJMUlLSOe9l2rRpRpJZvXr1OY/ds2ePGThwoImPjzelSpUyVatWNddee62ZMmXKWe/7zOtk5gyBxJ31XGOMGTNmjDnvvPNMSEhItu+MOf0zqlixoomOjjbHjx8/570F8iwSEhJMp06d/J8DuZfc4s7tuWW9/zVr1piwsDAzePBg2zGnTp0yLVq0MNWrVzcHDx48572OGzfOlCtXzhw7dswYY8zzzz9vJJn169fn6dwmTZqYKlWqmIcfftikp6cbY4yRZGJjY01ycnK2c1q1auXPB9955x1z1VVX2b6/4oor/N9PmTLFNGjQwERFRZno6GgjyWzfvt12/GOPPWZq1qxpfD7fOeMFENwoUKLQPPvss0aSGTZsmNuhoJj74IMPjCSzZMkSt0MJWs8++6ypV6+eycjICOi8tLQ0c95555nXXnvNpKenmwkTJpjQ0FB/YpiYmGhGjx5tTp06ZWbMmGHCwsLMzJkzzYYNG0y5cuXM9u3b/f9zdmaB8lztTpo0ycTHx5sTJ044+hwAAADOREf32RWkozvTjh07TKNGjcy0adOMMacLlMuWLTPx8fHm3XfftR17ZoHybB3a5+oMN4YOb8BreEkOCs1jjz0mY4zGjh3rdigopr799ltNnTpVQ4cOVdOmTdWqVSu3Qwpajz32mH788UeFhAT213p4eLg++ugjvfLKK4qJidH69ettL6158cUX9fbbb6tSpUpas2aN/7tGjRrpnnvuUePGjVWrVi3Vrl1b4eHheW733nvv1Y4dOxQREVHAOwcAAMjdnDlztG/fPtvL/PA/0dHR+vvf/64XXnghoDeNf//991q9erVOnTql8uXLq1SpUrY3tteuXVvz58/XQw89pC+++CLHNhITE3Xy5ElNmTJFJ0+e1MSJE/Xnn39Kko4cOaKQkBDFxsbq1KlTGjFiRLbzp02bplKlSunee+8N8K4BBCPLmCyr+AJAkOjbt6/eeustXXLJJZo+fXq2hdhROM5cvBwAAKA4+vbbb7V+/XqNGTNGlStXDmhNeJzbokWLNGTIEP/66z169ND48eMVGhoqy7K0c+dO1ahRQ2vWrFGnTp303nvvqVWrVtnyzG+//VZ33nmnfvvtN/Xq1UubN2/W3Xffrd69e+uRRx7R1KlTVbZsWf3f//2fHnroIf3444+qVauWuzcPoFBQoAQA2FCgBAAAxR0d3QBQvPAWbwAAAACAp0yfPl3Tp093OwwAQB4xghIAAAAAAACAa4p8BKXP59OuXbtUvnx5WZZV1JcHAACFxBijw4cPq3r16gG/sAkIJuSrAAB4E/lq8CryAuWuXbsUHx9f1JcFAABFJHNhfKC4Il8FAMDbyFeDT5EXKMuXLy9JmmbVVhmr5FWr3x/8qdshuObOXhXcDsE1VSL2ux2Ca0KU4XYIAIrIkSNHdNXV1/j/rQeKq5Ker3744Gduh+CafrdWcDsE11SN2Od2CK4hXwVKDvLV4FXkBcrMaTJlrBCVsUKL+vKuKxVRcv8QlC0X5XYIrikfmeZ2CK4h4QNKHqbEorgjXyVfLYnKR55wOwTXkK8CJQ/5avApeV3CAAAAAAAAAIIGBUoAAAAAAAAArqFACQAAAAAAAMA1FCgBAAAAAAAAuIYCJQAAAAAAAADXUKAEAAAAAAAA4BoKlAAAAAAAAABcQ4ESAAAAAAAAgGsoUAIAAAAAAABwDQVKAAAAAAAAAK6hQAkAAAAAAADANRQoAQAAAAAAALiGAiUAAAAAAAAA11CgBAAAAAAAAOAaCpQAAAAAAAAAXEOBEgAAAAAAAIBrKFACAAAAAAAAcA0FSgAAAAAAAACuoUAJAAAAAAAAwDUUKAEAAAAAAAC4hgIlAAAAAAAAANdQoAQAAAAAAADgGgqUAAAAAAAAAFxDgRIAAAAAAACAayhQAgAAAAAAAHBNvgqUEydOVK1atRQZGanLL79c3333ndNxAQAAAPlGvgoAAFB8BFygfO+99zR06FCNGDFCa9euVZMmTdShQwft3bu3MOIDAAAAAkK+CgAAULwEXKAcN26c7rrrLvXr108NGjTQ5MmTVaZMGb3xxhuFER8AAAAQEPJVAACA4iWgAmV6errWrFmjdu3a/a+BkBC1a9dOK1euzPGctLQ0paam2jYAAACgMJCvAgAAFD8BFSj379+vjIwMValSxba/SpUq2r17d47nJCcnKzo62r/Fx8fnP1oAAADgLMhXAQAAip9Cf4v3sGHDlJKS4t927txZ2JcEAAAA8ox8FQAAwF1hgRxcuXJlhYaGas+ePbb9e/bsUdWqVXM8JyIiQhEREfmPEAAAAMgj8lUAAIDiJ6ARlOHh4WrWrJkWLVrk3+fz+bRo0SIlJiY6HhwAAAAQCPJVAACA4iegEZSSNHToUCUlJal58+a67LLLNH78eB09elT9+vUrjPgAAACAgJCvAgAAFC8BFyhvvfVW7du3T8OHD9fu3bt1ySWXaN68edkWIgcAAADcQL4KAABQvARcoJSkQYMGadCgQU7HAgAAADiCfBUAAKD4KPS3eAMAAAAAAABAbihQAgAAAAAAAHANBUoAAAAAAAAArqFACQAAAAAAAMA1FCgBAAAAAAAAuIYCJQAAAAAAAADXUKAEAAAAAAAA4BoKlAAAAAAAAABcQ4ESAAAAAAAAgGsoUAIAAAAAAABwDQVKAAAAAAAAAK4JczsAAAAAJ504cULp6emOtRceHq7IyEjH2gMAAEDJRr6aHQVKAADgGSdOnFD10uV0UBmOtVm1alVt27at2Cd9AAAAcB/5as4oUAIAAM9IT0/XQWVoRmQdlXFgJZtj8ilp969KT08v1gkfAAAAggP5as4oUAIAAM8poxCVsUIL3pApeBMAAABAVuSrdhQoAQCA51hhlkIsq+DtmIK3AQAAAGRFvmrHW7wBAAAAAAAAuIYRlAAAwHOsUiGyrIL3w1rGI3NmAAAAEFTIV+0oUAIAAM8JCbUUElLw6S4hPm9MmQEAAEBwIV+1Y4o3AAAAAAAAANcwghIAAHiOVcqS5UCPtOWRHmkAAAAEF/JVOwqUAADAc0LCmDIDAACA4EW+ascUbwAAAAAAAACuYQQlAADwHKbMAAAAIJiRr9q5VqB8f/CnKhVR3q3Lu2bIsm5uh+CaanfOcDsE1+w7VcXtEFxTIeyQ2yG4ypLP7RCAEikk1FJIqANTZjK8kfAhfz588LMSma8+sKSr2yG4plp/8tWSiHyVfBVwA/mqHVO8AQAAAAAAALiGKd4AAMBzrFBLlgM90pa80SMNAACA4EK+ascISgAAAAAAAACuYQQlAADwHMfW9PFIjzQAAACCC/mqHQVKAADgOVaIQ29FNN5I+AAAABBcyFftmOINAAAAAAAAwDWMoAQAAJ5jhYbICi14P6wl40A0AAAAgB35qh0FSgAA4Dms6QMAAIBgRr5qxxRvAAAAAAAAAK5hBCUAAPAcy3Jo0XGfN3qkAQAAEFzIV+0oUAIAAM+xQuXIlBnLG0v6AAAAIMiQr9oxxRsAAAAAAACAaxhBCQAAPMcKtWQ50iPtjSkzAAAACC7kq3aMoAQAAAAAAADgGkZQAgAAz7FCQmSFFLwf1ok2AAAAgKzIV+0oUAIAAM+xQhx6K6IDbQAAAABZka/aeaPMCgAAAAAAAKBYYgQlAADwnJBQSyEOLDoe4pFFxwEAABBcyFftKFACAADPYcoMAAAAghn5qh1TvAEAAAAAAAC4hhGUAADAcyzLobciWvTlAgAAwHnkq3YUKAEAgOcwZQYAAADBjHzVzhtlVgAAAAAAAADFEiMoAQCA5zj2VkSfN3qkAQAAEFzIV+0YQQkAAAAAAADANYygBAAAnsOaPgAAAAhm5Kt2FCgBAIDnWCEOvRXRgTYAAACArMhX7QK+i2XLlqlz586qXr26LMvSnDlzCiEsAAAAIH/IVwEAAIqXgAuUR48eVZMmTTRx4sTCiAcAAKDAMqfMOLGh+CFfBQAAwc6tfDU5OVktWrRQ+fLlFRcXpy5dumjLli22Y06cOKGBAwcqJiZG5cqVU7du3bRnzx4nbz+bgAuUHTt21NNPP62bb765MOIBAAAoMBK+ko18FQAABDu38tWlS5dq4MCBWrVqlRYsWKCTJ0+qffv2Onr0qP+Yhx56SJ999pk++OADLV26VLt27VLXrl2dfgQ2hb4GZVpamtLS0vyfU1NTC/uSAAAArshM+Fq0aKFTp07piSeeUPv27bVp0yaVLVtW0umE74svvtAHH3yg6OhoDRo0SF27dtU333zjcvQlF/kqAAAoKebNm2f7PH36dMXFxWnNmjW65pprlJKSotdff12zZs1S27ZtJUnTpk3TRRddpFWrVumKK64olLgKvUCZnJysUaNGFfZlAAAA/Jx+K2LWglVERIQiIiKyHR+sCR/OjnwVAAAUNbfy1axSUlIkSZUqVZIkrVmzRidPnlS7du38x9SvX181a9bUypUrCy1fLfRX/QwbNkwpKSn+befOnYV9SQAAAEfFx8crOjravyUnJ+fpvEATPriDfBUAABR3+clXfT6fhgwZopYtW6pRo0aSpN27dys8PFwVKlSwHVulShXt3r27MEKXVAQjKPNasQUAAHDK6R7pgvfDZvZI79y5U1FRUf79ecltginhw9mRrwIAgKIWDPnqwIEDtXHjRi1fvrzAcRRUoRcoAQAAipoVYikk1IEpMxmn24iKirIlfHkRTAkfAAAAgovb+eqgQYP0+eefa9myZapRo4Z/f9WqVZWenq5Dhw7ZOtX37NmjqlWrFjje3ARcqj1y5IjWrVundevWSZK2bdumdevWaceOHU7HBgAAUCxlJnyLFy/ONeE7U2EnfCUN+SoAAEDOjDEaNGiQZs+erX//+9+qXbu27ftmzZqpVKlSWrRokX/fli1btGPHDiUmJhZaXAGPoPz+++/Vpk0b/+ehQ4dKkpKSkjR9+nTHAgMAAMgvpxcdzytjjAYPHqzZs2dryZIlZ034unXrJqloEr6ShnwVAAAEO7fy1YEDB2rWrFn65JNPVL58ef8yQ9HR0SpdurSio6M1YMAADR06VJUqVVJUVJQGDx6sxMTEQn2hY8AFytatW8sYUxixAAAAOMIKCXFoTZ/A2gjWhK+kIV8FAADBzq18ddKkSZJO50tnmjZtmvr27StJevHFFxUSEqJu3bopLS1NHTp00KuvvlrgWM+GNSgBAAAcEqwJHwAAACApT524kZGRmjhxoiZOnFgEEZ1GgRIAAHiOm1O8z8WNhA8AAADBxa18NVhRoAQAAJ5DwgcAAIBgRr5qV/DJ7gAAAAAAAACQT4ygBAAAnuPWouMAAABAXpCv2nnjLgAAAAAAAAAUS4ygBAAAnsOaPgAAAAhm5Kt2FCgBAIDnMGUGAAAAwYx81c4bdwEAAAAAAACgWGIEJQAA8B7LOr050Q4AAADgNPJVGwqUAADAcyzLoTV9PJLwAQAAILiQr9oxxRsAAAAAAACAaxhBCQAAPIdFxwEAABDMyFftKFACAADPsUIcmjLjQBsAAABAVuSrdt4oswIAAAAAAAAolhhBCQAAPIcpMwAAAAhm5Kt23rgLAAAAAAAAAMUSIygBAIDnWCHOrMdj0ZULAACAQkC+akeBEgAAeA6LjgMAACCYka/aeaTOCgAAAAAAAKA4YgQlAADwnpCQ05sT7QAAAABOI1+1oUAJAAA8x7IsWZYDU2YcaAMAAADIinzVzhtlVgAAAAAAAADFEiMoAQCA51ghIbIcmO7iRBsAAABAVuSrdq4VKO/sVUFly0W5dXnXVLtzhtshuKbSvi1uh+CaY3Hl3A7BNSdNKbdDcFUp66TbIbjGks/tEACgQPrdWkLz1f7kqyUR+WrJRb4KIBgwghIAAHiOFWLJCnFgTR8H2gAAAACyIl+1o0AJAAC8x3LorYiWN6bMAAAAIMiQr9p44y4AAAAAAAAAFEuMoAQAAN7j0JQZeWTKDAAAAIIM+aoNBUoAAOA5lhUiy4HpLk60AQAAAGRFvmrnjbsAAAAAAAAAUCwxghIAAHhPiOXMdBePTJkBAABAkCFftaFACQAAPMcKCZHlwFsRnWgDAAAAyIp81c4bdwEAAAAAAACgWGIEJQAA8BzLobciOvJmRQAAACAL8lU7RlACAAAAAAAAcA0jKAEAgPdYlmQ50A9reaNHGgAAAEGGfNWGAiUAAPAcpswAAAAgmJGv2jHFGwAAAAAAAIBrGEEJAAC8JyTk9OZEOwAAAIDTyFdtKFACAADPsSxLlgPr8TjRBgAAAJAV+aqdN8qsAAAAAAAAAIolRlACAADvsRyaMuPEmxUBAACArMhXbShQAgAAz+GtiAAAAAhm5Kt23iizAgAAAAAAACiWGEEJAAC8xwpxZrqLR6bMAAAAIMiQr9p44y4AAAAAAAAAFEuMoAQAAN4TYp3enGgHAAAAcBr5qg0FSgAA4DmWFSLLgekuTrQBAAAAZEW+aueNuwAAAAAAAABQLDGCEgAAeA9TZgAAABDMyFdtKFACAADPsUJCZIU4MGXGgTYAAACArMhX7QK6i+TkZLVo0ULly5dXXFycunTpoi1bthRWbAAAAEBAyFcBAACKn4AKlEuXLtXAgQO1atUqLViwQCdPnlT79u119OjRwooPAAAgcJbl3IZihXwVAAAUC+SrNgFN8Z43b57t8/Tp0xUXF6c1a9bommuucTQwAAAAIFDkqwAAAMVPgdagTElJkSRVqlQp12PS0tKUlpbm/5yamlqQSwIAAJxbiCU5sR6PRxYdL8nIVwEAQFAiX7XJ95Pw+XwaMmSIWrZsqUaNGuV6XHJysqKjo/1bfHx8fi8JAACQN0yZgchXAQBAECNftcl3gXLgwIHauHGj3n333bMeN2zYMKWkpPi3nTt35veSAAAAQJ6RrwIAABQP+ZriPWjQIH3++edatmyZatSocdZjIyIiFBERka/gAAAA8sMKCZHlwJQZJ9qAO8hXAQBAMCNftQuoQGmM0eDBgzV79mwtWbJEtWvXLqy4AAAA8s8KOb050Q6KFfJVAABQLJCv2gRUoBw4cKBmzZqlTz75ROXLl9fu3bslSdHR0SpdunShBAgAAADkFfkqAABA8RNQmXXSpElKSUlR69atVa1aNf/23nvvFVZ8AAAAgbOs//9mxAJuHll0vCQhXwUAAMUC+apNwFO8AQAAgp1lhchyYLqLE22gaJGvAgCA4oB81c4bdwEAAAAAAACgWMrXW7wBAACCWuaUFyfaAQAAAJxGvmrDCEoAAAAAAAAArmEEJQAA8B4r5PTmRDsAAACA08hXbShQAgAA77EceqOhR96KCAAAgCBDvmrjjTIrAAAAAAAAgHNatmyZOnfurOrVq8uyLM2ZM8f2vTFGw4cPV7Vq1VS6dGm1a9dOW7duLdSYKFACAADvCQlxbgtAMCZ7AAAACEIu5auSdPToUTVp0kQTJ07M8fvnn39eL7/8siZPnqxvv/1WZcuWVYcOHXTixImC3nWuKFACAADvyVzTx4ktAMGY7AEAACAIuZSvSlLHjh319NNP6+abb872nTFG48eP11NPPaWbbrpJjRs31ptvvqldu3Zl63x3EmtQAgAAOKRjx47q2LFjjt9lTfYk6c0331SVKlU0Z84c9ezZsyhDBQAAgIekpqbaPkdERCgiIiLgdrZt26bdu3erXbt2/n3R0dG6/PLLtXLlykLLWRlBCQAAvCfEcm7T6YTvzC0tLS3gkM6V7AEAAKAEcThfjY+PV3R0tH9LTk7OV1i7d++WJFWpUsW2v0qVKv7vCgMjKAEAgPdYVr6mu+TYjk4nfGcaMWKERo4cGVBTbiV7AAAACEIO56s7d+5UVFSUf3d+Rk+6iQIlAADAORT3hA8AAADeFhUVZctX86tq1aqSpD179qhatWr+/Xv27NEll1xS4PZzwxRvAADgPZbl3Kb/JXyZW34KlGcme2fas2eP/zsAAACUEA7nq06pXbu2qlatqkWLFvn3paam6ttvv1ViYqKj1zoTBUoAAIAi4FayBwAAAJzpyJEjWrdundatWyfp9Frp69at044dO2RZloYMGaKnn35an376qTZs2KA+ffqoevXq6tKlS6HFxBRvAADgPSEhpzcn2gnAkSNH9PPPP/s/ZyZ7lSpVUs2aNf3J3gUXXKDatWvr//7v/wo92QMAAEAQcilflaTvv/9ebdq08X8eOnSoJCkpKUnTp0/X3//+dx09elR33323Dh06pKuuukrz5s1TZGRkwePNBQVKAADgPU5NdwmwjWBM9gAAABCEXMpXJal169YyxpylSUujR4/W6NGjCxJZQChQAgAAOCQYkz0AAAAg2FGgBAAA3mOFnN6caAcAAABwGvmqDQVKAADgPZZDa/p4JOEDAABAkCFftfHGXQAAAAAAAAAolhhBCQAAvMfFRccBAACAcyJftaFACQAAvIc1fQAAABDMyFdtvHEXAAAAAAAAAIolRlACAADvYcoMAAAAghn5qo1rBcoqEftVPjLNrcu7Zt+pKm6H4JpjceXcDsE18X+ucjsE1+ysdoXbIbjqlCm5/UBh1im3Q3CNJZ/bIQBwQNWIfSofecLtMIoc+WrJRL5acpGvlkzkqwg2JfdvIgAA4F0hIac3J9oBAAAAnEa+akOBEgAAeI6xLBkHprs40QYAAACQFfmqnTfKrAAAAAAAAACKJUZQAgAA77EsyXKgH9YjPdIAAAAIMuSrNhQoAQCA91ghDiV8TDYBAABAISBftfHGXQAAAAAAAAAolhhBCQAAPIdFxwEAABDMyFftGEEJAAAAAAAAwDWMoAQAAN7Dmj4AAAAIZuSrNhQoAQCA91iWM2809MiUGQAAAAQZ8lUbb5RZAQAAAAAAABRLjKAEAADeExJyenOiHQAAAMBp5Ks2FCgBAIDn8FZEAAAABDPyVTtvlFkBAAAAAAAAFEuMoAQAAN7DWxEBAAAQzMhXbShQAgAAzzFWiIwDyZoTbQAAAABZka/aeeMuAAAAAAAAABRLjKAEAADeY1mnNyfaAQAAAJxGvmrDCEoAAAAAAAAArmEEJQAA8Bwjh9b0oS8XAAAAhYB81Y4CJQAA8B6mzAAAACCYka/aeKPMCgAAAAAAAKBYYgQlAADwHsuSHJgy45UeaQAAAAQZ8lUbCpQAAMBzjGXJOJCsOdEGAAAAkBX5qh1TvAEAAAAAAAC4hhGUAADAe6wQh6bM0JcLAACAQkC+akOBEgAAeI6RJSMHpsw40AYAAACQFfmqXUBl1kmTJqlx48aKiopSVFSUEhMTNXfu3MKKDQAAAAgI+SoAAEDxE9AIyho1aujZZ5/VBRdcIGOMZsyYoZtuukn/+c9/1LBhw8KKEQAAICDGCpFxYLqLE22gaJGvAgCA4oB81S6gAmXnzp1tn5955hlNmjRJq1atIuEDAACA68hXAQAAip98r0GZkZGhDz74QEePHlViYmKux6WlpSktLc3/OTU1Nb+XBAAAyBsWHYfIVwEAQBAjX7UJuEC5YcMGJSYm6sSJEypXrpxmz56tBg0a5Hp8cnKyRo0aVaAgAQAAAmEsS8ZyYNFxB9pA0SNfBQAAwY581S7gMmu9evW0bt06ffvtt7rvvvuUlJSkTZs25Xr8sGHDlJKS4t927txZoIABAACAsyFfBQAAKF4CHkEZHh6uunXrSpKaNWum1atX66WXXtJrr72W4/ERERGKiIgoWJQAAAABYNHxko18FQAABDvyVbt8r0GZyefz2dbsAQAAcJ1lnd6caAfFHvkqAAAIOuSrNgEVKIcNG6aOHTuqZs2aOnz4sGbNmqUlS5Zo/vz5hRUfAAAAkGfkqwAAAMVPQAXKvXv3qk+fPvrzzz8VHR2txo0ba/78+bruuusKKz4AAIDAOTRlxitvRSxJyFcBAECxQL5qE1CB8vXXXy+sOAAAAIACI18FAAAofgq8BiUAAECwMbJkVPD1eJxoAwAAAMiKfNWOAiUAAPAc3ooIAACAYEa+aueNuwAAAAAAAABQLDGCEgAAeI8lyXJguos3ZswAAAAg2JCv2lCgBAAAnmMUIuPARBEn2gAAAACyIl+188ZdAAAAAAAAACiWGEEJAAA8x1iWjANTZpxoAwAAAMiKfNWOAiUAAPAc3ooIAACAYEa+aueNuwAAAAAAAABQLDGCEgAAeI6RJePAKw2daAMAAADIinzVjhGUAAAAAAAAAFzDCEoAAOA5rOkDAACAYEa+akeBEgAAeA5vRQQAAEAwI1+180aZFQAAAAAAAECxxAhKAADgOSw6DgAAgGBGvmpHgRIAAHgOa/oAAAAgmJGv2nnjLgAAAAAAAAAUS4ygBAAAnsOUGQAAAAQz8lU7CpQAAMBzjByaMsNkEwAAABQC8lU7b9wFAAAAAAAAgGKJEZQAAMBzmDIDAACAYEa+ascISgAAAIdNnDhRtWrVUmRkpC6//HJ99913bocEAAAA+AVbvkqBEgAAeI6xLBkrxIEt8B7p9957T0OHDtWIESO0du1aNWnSRB06dNDevXsL4U4BAABQHJGv2lGgBAAAnpM5ZcaJLVDjxo3TXXfdpX79+qlBgwaaPHmyypQpozfeeKMQ7hQAAADFEfmqnWtrUIYoQyHKcOvyrqkQdsjtEFxz0pRyOwTX7Kx2hdshuCb2+A63Q3BVSpkqbofgmpMm3O0QXBNqlbx/3yTJp1C3Qyg0qampts8RERGKiIjIdlx6errWrFmjYcOG+feFhISoXbt2WrlyZaHHCWeRr5Y85KslE/kq+WpJRL7qPcU9X2UEJQAA8JzTU2ac2SQpPj5e0dHR/i05OTnH6+7fv18ZGRmqUsX+P3tVqlTR7t27C/2+AQAAUDyQr9rxFm8AAOA5xlgyxoG3Iv7/Nnbu3KmoqCj//px6owEAAIC8Il+1o0AJAABwDlFRUbaELzeVK1dWaGio9uzZY9u/Z88eVa1atbDCAwAAQAlX3PNVpngDAAAPCpFxYAs0VQoPD1ezZs20aNEi/z6fz6dFixYpMTHR4XsEAABA8UW+eiZGUAIAADho6NChSkpKUvPmzXXZZZdp/PjxOnr0qPr16+d2aAAAAEBQ5qsUKAEAgOcYWTJyYE2ffLRx6623at++fRo+fLh2796tSy65RPPmzcu2EDkAAABKLvJVOwqUAADAc9xM+CRp0KBBGjRoUIGvDwAAAG8iX7VjDUoAAAAAAAAArmEEJQAA8By3e6QBAACAsyFftaNACQAAPIeEDwAAAMGMfNWOKd4AAAAAAAAAXMMISgAA4DnGWDLGgR5pB9oAAAAAsiJftaNACQAAPIcpMwAAAAhm5Kt2TPEGAAAAAAAA4BpGUAIAAM+hRxoAAADBjHzVjhGUAAAAAAAAAFzDCEoAAOA59EgDAAAgmJGv2lGgBAAAnmPk0FsRPZLwAQAAILiQr9oxxRsAAAAAAACAaxhBCQAAPMcnSz4HepOdaAMAAADIinzVjgIlAADwHNb0AQAAQDAjX7VjijcAAAAAAAAA1zCCEgAAeI4xDi067kAbAAAAQFbkq3YUKAEAgOcYOTPdxRQ8FAAAACAb8lU7pngDAAAAAAAAcA0jKAEAgOcwZQYAAADBjHzVjhGUAAAAAAAAAFzDCEoAAOA5RpZDa/p4o0caAAAAwYV81a5AIyifffZZWZalIUOGOBQOAABAwWVOmXFiQ/FGvgoAAIIR+apdvguUq1ev1muvvabGjRs7GQ8AAADgCPJVAACA4iFfBcojR47o9ttv19SpU1WxYsWzHpuWlqbU1FTbBgAAUJiMJJ8DmynqwOEY8lUAABDMyFft8lWgHDhwoDp16qR27dqd89jk5GRFR0f7t/j4+PxcEgAAIM+YMgPyVQAAEMzIV+0CLlC+++67Wrt2rZKTk/N0/LBhw5SSkuLfdu7cGXCQAAAAQF6RrwIAABQvAb3Fe+fOnXrwwQe1YMECRUZG5umciIgIRURE5Cs4AACA/OCtiCUX+SoAACgOyFftAipQrlmzRnv37tWll17q35eRkaFly5ZpwoQJSktLU2hoqONBAgAABMKp6S5emTJTkpCvAgCA4oB81S6gAuW1116rDRs22Pb169dP9evX12OPPUayBwAAAFeRrwIAABQ/ARUoy5cvr0aNGtn2lS1bVjExMdn2AwAAuIUpMyUX+SoAACgOyFft8vUWbwAAAAAAAABwQkAjKHOyZMkSB8IAAABwjs+c3pxoB8Uf+SoAAAg25Kt2BS5QAgAABBumzAAAACCYka/aMcUbAAAAAAAAgGsYQQkAADzHGEvGONAj7UAbAAAAQFbkq3YUKAEAgOcYc3pzoh0AAADAaeSrdkzxBgAAAAAAAOAaRlACAADP8cmSz4EFw51oAwAAAMiKfNWOEZQAAAAAAAAAXMMISgAA4DksOg4AAIBgRr5qR4ESAAB4DouOAwAAIJiRr9oxxRsAAAAAAACAaxhBCQAAPMfIknFgwXAn2gAAAACyIl+1o0AJAAA8x2dOb060AwAAADiNfNWOKd4AAAAAAAAAXMMISgAA4D0OvRVRHnkrIgAAAIIM+aoNBUoAAOA5vBURAAAAwYx81Y4p3gAAAAAAAABcwwhKAADgOT5Z8jnwRkMn2gAAAACyIl+1YwQlAAAAAAAAANcwghIAAHgOa/oAAAAgmJGv2lGgBAAAnmMceiuiI29WBAAAALIgX7VjijcAAAAAAAAA1zCCEgAAeI7PnN6caAcAAABwGvmqHQVKAADgOazpAwAAgGBGvmrHFG8AAAAAAAAArmEEZRGz5HM7BNeUsk66HYJrTpmS+0ctpUwVt0NwVeSpo26H4Bor1CNdefmQZiLdDsEVGUH0d52RJSMHFh13oA2guCFfLZnIV0su8tWSiXzVfeSrdsHzkwEAAHCITw6t6VPwJgAAAIBsyFftmOINAAAAAAAAwDWMoAQAAJ7DouMAAAAIZuSrdoygBAAAAAAAAOAaRlACAADPoUcaAAAAwYx81Y4CJQAA8ByfseQzBX+joRNtAAAAAFmRr9oxxRsAAAAAAACAaxhBCQAAPIcpMwAAAAhm5Kt2FCgBAIDnkPABAAAgmJGv2jHFGwAAAAAAAIBrKFACAADPMUbyObAVZo/0M888oyuvvFJlypRRhQoVcjxmx44d6tSpk8qUKaO4uDg9+uijOnXqVOEFBQAAgCJRHPLVokSBEgAAwAXp6enq3r277rvvvhy/z8jIUKdOnZSenq4VK1ZoxowZmj59uoYPH17EkQIAAKCkKqpOddagBAAAnmOMJWMsR9opLKNGjZIkTZ8+Pcfvv/rqK23atEkLFy5UlSpVdMkll2jMmDF67LHHNHLkSIWHhxdabAAAAChcxSFflf7XqZ6YmKjXX3892/eZnepVq1bVihUr9Oeff6pPnz4qVaqUxo4dm+frMIISAAB4Tuai405skpSammrb0tLSCv0eVq5cqYsvvlhVqlTx7+vQoYNSU1P1ww8/FPr1AQAAUHiKS746atQoPfTQQ7r44otz/D6zU/2tt97SJZdcoo4dO2rMmDGaOHGi0tPT83wdCpQAAADnEB8fr+joaP+WnJxc6NfcvXu3rTgpyf959+7dhX59AAAAFB9u5KuSc53qTPEGAACek7louBPtSNLOnTsVFRXl3x8REZHj8Y8//riee+65s7a5efNm1a9fv+DBAQAAoNhyK191mlOd6hQoAQCA55w53aWg7UhSVFSULeHLzcMPP6y+ffue9Zg6derk6dpVq1bVd999Z9u3Z88e/3cAAAAovtzKV6Xg7FSnQAkAAOCQ2NhYxcbGOtJWYmKinnnmGe3du1dxcXGSpAULFigqKkoNGjRw5BoAAAAoeYKxU50CJQAA8Byne6QLw44dO/TXX39px44dysjI0Lp16yRJdevWVbly5dS+fXs1aNBAd9xxh55//nnt3r1bTz31lAYOHFhkU3YAAABQONzMV4OxU50CJQAA8Byn1/QpDMOHD9eMGTP8n5s2bSpJWrx4sVq3bq3Q0FB9/vnnuu+++5SYmKiyZcsqKSlJo0ePLrygAAAAUCSKQ74qFV2nOgVKAAAAF0yfPl3Tp08/6zEJCQn68ssviyYgAAAAIIui6lSnQAkAADynOEzxBgAAQMlVXPLVoupUDynQ2QAAAAAAAABQAIygBAAAnuPznd6caAcAAABwGvmqHQVKAADgOcVlygwAAABKJvJVO6Z4AwAAAAAAAHANIygBAIDn0CMNAACAYEa+ahfQCMqRI0fKsizbVr9+/cKKDQAAIF98knzGgc3tG0HAyFcBAEBxQL5qF/AIyoYNG2rhwoX/ayCMQZgAAAAIHuSrAAAAxUvA2VpYWJiqVq2a5+PT0tKUlpbm/5yamhroJQEAAAJijJFxYL6LE22g6JGvAgCAYEe+ahfwS3K2bt2q6tWrq06dOrr99tu1Y8eOsx6fnJys6Oho/xYfH5/vYAEAAPIic00fJzYUP+SrAAAg2JGv2gVUoLz88ss1ffp0zZs3T5MmTdK2bdt09dVX6/Dhw7meM2zYMKWkpPi3nTt3FjhoAAAAICfkqwAAAMVPQFO8O3bs6P/vxo0b6/LLL1dCQoLef/99DRgwIMdzIiIiFBERUbAoAQAAAmB8ks+BFcONV1YdL0HIVwEAQHFAvmoX8BTvM1WoUEEXXnihfv75Z6fiAQAAABxDvgoAABD8ClSgPHLkiH755RdVq1bNqXgAAAAKjDV9kIl8FQAABCPyVbuACpSPPPKIli5dqu3bt2vFihW6+eabFRoaql69ehVWfAAAAAHzGec2FC/kqwAAoDggX7ULaA3K33//Xb169dKBAwcUGxurq666SqtWrVJsbGxhxQcAAADkGfkqAABA8RNQgfLdd98trDgAAAAc49R0F69MmSlJyFcBAEBxQL5qF1CBEgAAoDgwPiPjwHwXJ9oAAAAAsiJftSvQS3IAAAAAAAAAoCAYQQkAADzHqQXDPdIhDQAAgCBDvmrHCEoAAAAAAAAArmEEJQAA8BwWHQcAAEAwI1+1o0AJAAA8x+cz8jkw38WJNgAAAICsyFftmOINAAAAAAAAwDWMoAQAAJ7DlBkAAAAEM/JVOwqUAADAc0j4AAAAEMzIV+2Y4g0AAAAAAADANYygBAAAnuMzRj4HupOdaAMAAADIinzVjgIlAADwHOM7vTnRDgAAAOA08lU7pngDAAAAAAAAcA0jKAEAgOcYGRkHprsYeWPKDAAAAIIL+aodIygBAAAAAAAAuIYRlAAAwHOMT/Kxpg8AAACCFPmqHQVKAADgOcY4NGXGI29FBAAAQHAhX7VjijcAAAAAAAAA1zCCEgAAeI7PnN6caAcAAABwGvmqHQVKAADgOcZnZBzI1pxoAwAAAMiKfNWOKd4AAAAAAAAAXMMISgAA4DnGnN6caAcAAABwGvmqHQVKAADgOT6fkc+B6S5OtAEAAABkRb5qxxRvAAAAAAAAAK5hBCWKjCWf2yG4Jsw65XYIrjlpwt0OwVVWqDd6s/Ij+uifbofgmm2RjdwOwRVHfRluh+BnjJFxYL6LE20AKD7IV0sm8tWS+28d+WrJQ74avBhBCQAAAAAAAMA1jKAEAACeY3ynNyfaAQAAAJxGvmpHgRIAAHiOzxj5HJju4kQbAAAAQFbkq3ZM8QYAAAAAAADgGkZQAgAAz2HRcQAAAAQz8lU7CpQAAMBzfD4jn8+BKTMOtAEAAABkRb5qxxRvAAAAAAAAAK5hBCUAAPAcY05vTrQDAAAAOI181Y4CJQAA8BxjjIwD0128sqYPAAAAggv5qh1TvAEAAAAAAAC4hhGUAADAc4wx8vFWRAAAAAQp8lU7RlACAAAAAAAAcA0jKAEAgOcYn0Nr+jjQBgAAAJAV+aodBUoAAOA5JHwAAAAIZuSrdkzxBgAAAAAAAOAaRlACAADP8ZnTmxPtAAAAAE4jX7WjQAkAADyHKTMAAAAIZuSrdkzxBgAAAAAAAOAaRlACAADPMcbIGAd6pB1oAwAAAMiKfNWOEZQAAAAAAAAAXMMISgAA4Dk+n+RzYD0en8+BYAAAAIAsyFftKFACAADPYcoMAAAAghn5qh1TvAEAAAAAAAC4hhGUAADAc4zPyDgwZcaJNgAAAICsyFftKFACAADPIeEDAABAMCNftWOKNwAAQBHbvn27BgwYoNq1a6t06dI6//zzNWLECKWnp9uOW79+va6++mpFRkYqPj5ezz//vEsRAwAAAIWHAiUAAPAcn4x8xoFNhdMj/eOPP8rn8+m1117TDz/8oBdffFGTJ0/WE0884T8mNTVV7du3V0JCgtasWaMXXnhBI0eO1JQpUwolJgAAABSdYM9XpaLtVA+4QPnHH3+od+/eiomJUenSpXXxxRfr+++/D/jCAAAAhSVzyowTW2G4/vrrNW3aNLVv31516tTRjTfeqEceeUQff/yx/5i3335b6enpeuONN9SwYUP17NlTDzzwgMaNG1coMXkJ+SoAAAh2wZ6vSkXbqR7QGpQHDx5Uy5Yt1aZNG82dO1exsbHaunWrKlasGNBFAQAAipPU1FTb54iICEVERDh6jZSUFFWqVMn/eeXKlbrmmmsUHh7u39ehQwc999xzOnjwIPlXLshXAQAAnHH99dfr+uuv93+uU6eOtmzZokmTJukf//iHJHunenh4uBo2bKh169Zp3Lhxuvvuu/N8rYAKlM8995zi4+M1bdo0/77atWsH0gQAAEChM8bIGAcWHf//bcTHx9v2jxgxQiNHjixw+5l+/vlnvfLKK/5ET5J2796dLc+qUqWK/zsKbjkjXwUAAMWB0/lqUXSoS4XXqR7QFO9PP/1UzZs3V/fu3RUXF6emTZtq6tSpZz0nLS1Nqamptg0AAKA42blzp1JSUvzbsGHDcjzu8ccfl2VZZ91+/PFH2zl//PGHrr/+enXv3l133XVXUdyOp5GvAgCAkig+Pl7R0dH+LTk52fFrZHaq33PPPf59u3fv9neiZzqzUz2vAhpB+euvv2rSpEkaOnSonnjiCa1evVoPPPCAwsPDlZSUlOM5ycnJGjVqVCCXAQAAKBDjM/I5sB5P5po+UVFRioqKOufxDz/8sPr27XvWY+rUqeP/7127dqlNmza68sors63TU7VqVe3Zs8e2L/Nz1apV8xJ+iUS+CgAAigOn89WdO3fa8tWzjZ58/PHH9dxzz5213c2bN6t+/fr+z4XdqR5QgdLn86l58+YaO3asJKlp06bauHGjJk+enGvCN2zYMA0dOtT/OTU1Nds0KQAAACc5tWB4oG3ExsYqNjY2T8f+8ccfatOmjZo1a6Zp06YpJMQ+sSUxMVFPPvmkTp48qVKlSkmSFixYoHr16jG9+yzIVwEAQHHgdL6a1w51KTg71QMqUFarVk0NGjSw7bvooov00Ucf5XpOYc15BwAAKK7++OMPtW7dWgkJCfrHP/6hffv2+b/LTORuu+02jRo1SgMGDNBjjz2mjRs36qWXXtKLL77oVtjFAvkqAADA2QVjp3pABcqWLVtqy5Yttn0//fSTEhISAmkGAACgUDm96LjTFixYoJ9//lk///yzatSokeM1o6Oj9dVXX2ngwIFq1qyZKleurOHDhwf0NsSSiHwVAAAUB8Ger0pF26keUIHyoYce0pVXXqmxY8eqR48e+u677zRlypRswzsBAADcZHw+GZ/PkXYKQ9++fc85rUaSGjdurK+//rpQYvAq8lUAAFAcBHu+KhVtp3pABcoWLVpo9uzZGjZsmEaPHq3atWtr/Pjxuv322wO6KAAAAFAYyFcBAACcUZSd6gEVKCXphhtu0A033FCgiwIAABQmn0NvRXSiDRQ98lUAABDsyFftAi5QAgAABLvisKYPAAAASi7yVbuQcx8CAAAAAAAAAIWDEZQAAMBzjM/IODDdxYk2AAAAgKzIV+0YQQkAAAAAAADANYygBAAAnkOPNAAAAIIZ+aodBUoAAOA5PvnkMz5H2gEAAACcRr5qxxRvAAAAAAAAAK5hBCUAAPAc43NmuosDndoAAABANuSrdhQoAQCA57CmDwAAAIIZ+aodU7wBAAAAAAAAuIYRlAAAwHOMMTLGgR5pB9oAAAAAsiJftWMEJQAAAAAAAADXMIISAAB4js/nk89X8BXDnWgDAAAAyIp81Y4CJQAA8BwWHQcAAEAwI1+1Y4o3AAAAAAAAANcwghIAAHiOMT4ZU/DpLk60AQAAAGRFvmpHgRIAAHgOU2YAAAAQzMhX7ZjiDQAAAAAAAMA1jKAEAADe41CPtDzSIw0AAIAgQ75qQ4ESAAB4js/45HNgPR4n2gAAAACyIl+1Y4o3AAAAAAAAANcwghIAAHgOi44DAAAgmJGv2jGCEgAAAAAAAIBrGEEJAAA8xxifjK/g6/EYj6zpAwAAgOBCvmpHgRIoApa88RdGfoRaGW6H4Ko0E+l2CK7ZFtnI7RBcc/7RtW6H4IrUo0fdDsGPKTMAEBjy1ZKLfLVkIl91H/mqHVO8AQAAAAAAALiGEZQAAMBzjPE5Mt3FK1NmAAAAEFzIV+0oUAIAAM/x+SSfA9NdHFgWCAAAAMiGfNWOKd4AAAAAAAAAXMMISgAA4DnG59BbEb3SJQ0AAICgQr5qR4ESAAB4Dm9FBAAAQDAjX7VjijcAAAAAAAAA1zCCEgAAeA5vRQQAAEAwI1+1YwQlAAAAAAAAANcwghIAAHgOa/oAAAAgmJGv2lGgBAAAnsNbEQEAABDMyFftKFACAADPyTh1NKjaAQAAAM5EvmpHgRIAAHhGeHi4qlatqu8X9XCszapVqyo8PNyx9gAAAFByka/mjAIlAADwjMjISG3btk3p6emOtRkeHq7IyEjH2gMAAEDJRb6aMwqUAADAUyIjI4t9ggYAAADvIl/NLsTtAAAAAAAAAACUXBQoAQAAAAAAALiGAiUAAAAAAAAA11CgBAAAAAAAAOAaCpQAAAAAAAAAXEOBEgAAAAAAAIBrKFACAAAAAAAAcA0FSgAAAAAAAACuoUAJAAAAAAAAwDUUKAEAAAAAAAC4hgIlAAAAAAAAANdQoAQAAAAAAADgmoAKlLVq1ZJlWdm2gQMHFlZ8AAAAQJ6RrwIAABQ/YYEcvHr1amVkZPg/b9y4Udddd526d+/ueGAAAABAoMhXAQAAip+ACpSxsbG2z88++6zOP/98tWrVytGgAAAAgPwgXwUAACh+AipQnik9PV1vvfWWhg4dKsuycj0uLS1NaWlp/s+pqan5vSQAAACQZ+SrAAAAxUO+X5IzZ84cHTp0SH379j3rccnJyYqOjvZv8fHx+b0kAAAAkGfkqwAAAMVDvguUr7/+ujp27Kjq1auf9bhhw4YpJSXFv+3cuTO/lwQAAADyjHwVAACgeMjXFO/ffvtNCxcu1Mcff3zOYyMiIhQREZGfywAAAAD5Qr4KAABQfORrBOW0adMUFxenTp06OR0PAAAAUGDkqwAAAMVHwAVKn8+nadOmKSkpSWFh+X7HDgAAAFAoyFcBAACKl4ALlAsXLtSOHTvUv3//wogHAAAAKBDyVQAAgOIl4C7l9u3byxhTGLEAAAAABUa+CgAAULzk+y3eAAAAAAAAAFBQFCgBAAAAAAAAuIYCJQAAAAAAAADXUKAEAAAAAAAA4BoKlAAAAAAAAABcQ4ESAAAAAAAAgGsoUAIAAAAAAABwDQVKAAAAAAAAAK6hQAkAAAAAAADANRQoAQAAAAAAALiGAiUAAAAAAAAA11CgBAAAAAAAAOAaCpQAAAAAAAAAXEOBEgAAAAAAAIBrKFACAAAAAAAAcA0FSgAAAAAAAACuoUAJAAAAAAAAwDUUKAEAAAAAAAC4hgIlAAAAAAAAANdQoAQAAAAAAADgGgqUAAAAAAAAAFxDgRIAAAAAAACAayhQAgAAAAAAAHBNWFFf0BgjSTpy5EhRXxqAC3wKdTsEV2WYIv9rNmgc9WW4HYJrUo8edTsEVxw+ekzS//6tB4or8lWgZCFfJV8tichXyVeDTZH/TXT48GFJ0lVXX1PUlwYAAEXg8OHDio6OdjsMIN/IVwEA8Dby1eBjmSIuG/t8Pu3atUvly5eXZVlFeWmlpqYqPj5eO3fuVFRUVJFe223ce8m8d6lk3z/3zr1z70XLGKPDhw+revXqCglhFRkUX+Sr7inJ98+9c+/ce8nBvZOvIrsiH0EZEhKiGjVqFPVlbaKiokrcXwKZuPeSee9Syb5/7p17L2ncvHd6ouEF5KvuK8n3z71z7yUN9869FzXy1eBEuRgAAAAAAACAayhQAgAAAAAAAHBNiSpQRkREaMSIEYqIiHA7lCLHvZfMe5dK9v1z79x7SVOS7x3wipL+57gk3z/3zr2XNNw79w6cqchfkgMAAAAAAAAAmUrUCEoAAAAAAAAAwYUCJQAAAAAAAADXUKAEAAAAAAAA4BoKlAAAAAAAAABcQ4ESAAAAAAAAgGtKVIFy4sSJqlWrliIjI3X55Zfru+++czukQrds2TJ17txZ1atXl2VZmjNnjtshFZnk5GS1aNFC5cuXV1xcnLp06aItW7a4HVaRmDRpkho3bqyoqChFRUUpMTFRc+fOdTssVzz77LOyLEtDhgxxO5RCN3LkSFmWZdvq16/vdlhF5o8//lDv3r0VExOj0qVL6+KLL9b333/vdlhFolatWtl+9pZlaeDAgW6HBiBAJTFflUpuzkq+Sr4qlax8VSJnLak5K/kqzqXEFCjfe+89DR06VCNGjNDatWvVpEkTdejQQXv37nU7tEJ19OhRNWnSRBMnTnQ7lCK3dOlSDRw4UKtWrdKCBQt08uRJtW/fXkePHnU7tEJXo0YNPfvss1qzZo2+//57tW3bVjfddJN++OEHt0MrUqtXr9Zrr72mxo0bux1KkWnYsKH+/PNP/7Z8+XK3QyoSBw8eVMuWLVWqVCnNnTtXmzZt0j//+U9VrFjR7dCKxOrVq20/9wULFkiSunfv7nJkAAJRUvNVqeTmrOSr5KslMV+VyFlLYs5KvopzsYwxxu0gisLll1+uFi1aaMKECZIkn8+n+Ph4DR48WI8//rjL0RUNy7I0e/ZsdenSxe1QXLFv3z7FxcVp6dKluuaaa9wOp8hVqlRJL7zwggYMGOB2KEXiyJEjuvTSS/Xqq6/q6aef1iWXXKLx48e7HVahGjlypObMmaN169a5HUqRe/zxx/XNN9/o66+/djuUoDBkyBB9/vnn2rp1qyzLcjscAHlEvnpaSc5ZyVfJV72er0rkrOSsp5GvIqsSMYIyPT1da9asUbt27fz7QkJC1K5dO61cudLFyFCUUlJSJJ1OfEqSjIwMvfvuuzp69KgSExPdDqfIDBw4UJ06dbL9uS8Jtm7dqurVq6tOnTq6/fbbtWPHDrdDKhKffvqpmjdvru7duysuLk5NmzbV1KlT3Q7LFenp6XrrrbfUv39/kj2gGCFfhUS+Sr5acpCzluyclXwVOSkRBcr9+/crIyNDVapUse2vUqWKdu/e7VJUKEo+n09DhgxRy5Yt1ahRI7fDKRIbNmxQuXLlFBERoXvvvVezZ89WgwYN3A6rSLz77rtau3atkpOT3Q6lSF1++eWaPn265s2bp0mTJmnbtm26+uqrdfjwYbdDK3S//vqrJk2apAsuuEDz58/XfffdpwceeEAzZsxwO7QiN2fOHB06dEh9+/Z1OxQAASBfBfkq+WpJQc5Kzkq+ipyEuR0AUBQGDhyojRs3lpi1TSSpXr16WrdunVJSUvThhx8qKSlJS5cu9XzSt3PnTj344INasGCBIiMj3Q6nSHXs2NH/340bN9bll1+uhIQEvf/++56fKuXz+dS8eXONHTtWktS0aVNt3LhRkydPVlJSksvRFa3XX39dHTt2VPXq1d0OBQAQAPJV8tWSgpyVnJV8FTkpESMoK1eurNDQUO3Zs8e2f8+ePapatapLUaGoDBo0SJ9//rkWL16sGjVquB1OkQkPD1fdunXVrFkzJScnq0mTJnrppZfcDqvQrVmzRnv37tWll16qsLAwhYWFaenSpXr55ZcVFhamjIwMt0MsMhUqVNCFF16on3/+2e1QCl21atWy/c/MRRddVGKmC2X67bfftHDhQt15551uhwIgQOSrJRv5KvlqSc1XJXLWkpazkq8iNyWiQBkeHq5mzZpp0aJF/n0+n0+LFi0qUWuclDTGGA0aNEizZ8/Wv//9b9WuXdvtkFzl8/mUlpbmdhiF7tprr9WGDRu0bt06/9a8eXPdfvvtWrdunUJDQ90OscgcOXJEv/zyi6pVq+Z2KIWuZcuW2rJli23fTz/9pISEBJcicse0adMUFxenTp06uR0KgACRr5ZM5Kt25KslL1+VyFlLWs5KvorclJgp3kOHDlVSUpKaN2+uyy67TOPHj9fRo0fVr18/t0MrVEeOHLH1RG3btk3r1q1TpUqVVLNmTRcjK3wDBw7UrFmz9Mknn6h8+fL+9Zuio6NVunRpl6MrXMOGDVPHjh1Vs2ZNHT58WLNmzdKSJUs0f/58t0MrdOXLl8+2blPZsmUVExPj+fWcHnnkEXXu3FkJCQnatWuXRowYodDQUPXq1cvt0ArdQw89pCuvvFJjx45Vjx499N1332nKlCmaMmWK26EVGZ/Pp2nTpikpKUlhYSXmn3fAU0pqviqV3JyVfJV8NVNJyVclctaSnLOSr+KsTAnyyiuvmJo1a5rw8HBz2WWXmVWrVrkdUqFbvHixkZRtS0pKcju0QpfTfUsy06ZNczu0Qte/f3+TkJBgwsPDTWxsrLn22mvNV1995XZYrmnVqpV58MEH3Q6j0N16662mWrVqJjw83Jx33nnm1ltvNT///LPbYRWZzz77zDRq1MhERESY+vXrmylTprgdUpGaP3++kWS2bNnidigACqAk5qvGlNyclXyVfDVTSclXjSFnLck5K/kqzsYyxpiiK4cCAAAAAAAAwP+UiDUoAQAAAAAAAAQnCpQAAAAAAAAAXEOBEgAAAAAAAIBrKFACAAAAAAAAcA0FSgAAAAAAAACuoUAJAAAAAAAAwDUUKAEAAAAAAAC4hgIlAAAAAAAAANdQoAQAAAAAAADgGgqUAAAAAAAAAFxDgRIAAAAAAACAa/4fYwDGzqp71jkAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "! pip install scikit-fda -q\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from skfda.misc.operators import LinearDifferentialOperator\n",
        "from skfda.misc.regularization import L2Regularization\n",
        "from skfda.representation import basis as skfda_basis\n",
        "\n",
        "\n",
        "def penalty_skfda(knots, degree):\n",
        "    basis = skfda_basis.BSplineBasis(knots=knots, order=degree + 1)\n",
        "    operator = LinearDifferentialOperator(1)\n",
        "    reg = L2Regularization(operator)\n",
        "    P = np.asarray(reg.penalty_matrix(basis))\n",
        "    return P, basis.n_basis, basis\n",
        "\n",
        "\n",
        "def penalty_quadrature(basis, knots, n_grid=8000):\n",
        "    K = basis.n_basis\n",
        "    xg = np.linspace(knots[0], knots[-1], n_grid)\n",
        "    dx = xg[1] - xg[0]\n",
        "\n",
        "    # Evaluate derivative of basis\n",
        "    Bp = basis(xg, derivative=1)\n",
        "    Bp = np.asarray(Bp)\n",
        "\n",
        "    # Reshape to (n_grid, K)\n",
        "    if Bp.ndim == 3:\n",
        "        if Bp.shape[1] == 1:        # (M,1,K)\n",
        "            Bp = Bp[:, 0, :]\n",
        "        elif Bp.shape[2] == 1:      # (K,M,1)\n",
        "            Bp = Bp[:, :, 0].T\n",
        "        else:\n",
        "            raise RuntimeError(f\"Cannot reshape Bp of shape {Bp.shape}\")\n",
        "    elif Bp.ndim != 2:\n",
        "        raise RuntimeError(f\"Cannot reshape Bp of shape {Bp.shape}\")\n",
        "\n",
        "    assert Bp.shape == (n_grid, K)\n",
        "\n",
        "    return (Bp.T @ Bp) * dx  # (K,K)\n",
        "\n",
        "\n",
        "def plot_matrix(ax, matrix, title, fig):\n",
        "    vmax = np.max(np.abs(matrix))\n",
        "    cax = ax.matshow(matrix, cmap='coolwarm', vmin=-vmax, vmax=vmax)\n",
        "    fig.colorbar(cax, ax=ax, shrink=0.8)\n",
        "    ax.set_title(title, fontsize=12)\n",
        "\n",
        "    n = matrix.shape[0]\n",
        "    ax.set_xticks(np.arange(n))\n",
        "    ax.set_yticks(np.arange(n))\n",
        "    ax.set_xticklabels(np.arange(n))\n",
        "    ax.set_yticklabels(np.arange(n))\n",
        "    ax.xaxis.tick_bottom()\n",
        "\n",
        "\n",
        "# 4. MAIN: compare and plot\n",
        "xi = np.array([0.00, 0.07, 0.21, 0.25, 0.60, 1.00])\n",
        "degree = 3\n",
        "\n",
        "# Compute both penalties\n",
        "P_skfda, K, basis = penalty_skfda(xi, degree)\n",
        "P_quad = penalty_quadrature(basis, xi)\n",
        "\n",
        "# Print difference metrics\n",
        "fro_err = np.linalg.norm(P_quad - P_skfda)\n",
        "max_err = np.max(np.abs(P_quad - P_skfda))\n",
        "\n",
        "print(\"Shape:\", P_quad.shape)\n",
        "print(\"Frobenius norm ||P_quad - P_skfda|| =\", fro_err)\n",
        "print(\"Max abs diff =\", max_err)\n",
        "\n",
        "\n",
        "fig, axes = plt.subplots(1, 2, figsize=(14, 6))\n",
        "fig.suptitle('Comparison of Penalty Matrices', fontsize=16)\n",
        "\n",
        "plot_matrix(axes[0], P_quad, r'Quadrature Penalty Matrix ($P_{\\rm quad}$)', fig)\n",
        "plot_matrix(axes[1], P_skfda, r'skfda Penalty Matrix ($P_{\\rm skfda}$)', fig)\n",
        "\n",
        "plt.tight_layout(rect=[0, 0, 1, 0.95])\n",
        "plt.show()\n"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}