svhn.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d0380a2f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "778ae8ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import torchvision\n",
    "import torchvision.transforms as transforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "3d805a7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set data path\n",
    "def set_data_path():\n",
    "    return \"../data/\"\n",
    "#     raise NotImplementedError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "9169089d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pre_process(torchset,n_samples,num_classes=10):\n",
    "    indices = list(np.random.choice(len(torchset),n_samples))\n",
    "\n",
    "    trainset = []\n",
    "    for ix in indices:\n",
    "        x,y = torchset[ix]\n",
    "        ohe_y = torch.zeros(num_classes)\n",
    "        ohe_y[y] = 1\n",
    "        trainset.append(((x/np.linalg.norm(x)).reshape(-1),ohe_y))\n",
    "    return trainset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "777dbd11",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using downloaded and verified file: ../data/train_32x32.mat\n",
      "Using downloaded and verified file: ../data/test_32x32.mat\n"
     ]
    }
   ],
   "source": [
    "# load svhn data\n",
    "transform = transforms.Compose([\n",
    "    transforms.ToTensor()\n",
    "])\n",
    "\n",
    "data_path = set_data_path() ## set this data path\n",
    "\n",
    "trainset0 = torchvision.datasets.SVHN(root=data_path,\n",
    "                                    split = \"train\",\n",
    "                                    transform=transform,\n",
    "                                    download=True)\n",
    "testset0 = torchvision.datasets.SVHN(root=data_path,\n",
    "                                    split = \"test\",\n",
    "                                    transform=transform,\n",
    "                                    download=True)\n",
    "\n",
    "trainset = pre_process(trainset0,n_samples=5000, num_classes=10)\n",
    "train_loader = torch.utils.data.DataLoader(trainset, batch_size=128,\n",
    "                                          shuffle=True, num_workers=2)\n",
    "\n",
    "\n",
    "testset = pre_process(testset0,n_samples=5000, num_classes=10)\n",
    "test_loader = torch.utils.data.DataLoader(testset, batch_size=128,\n",
    "                                         shuffle=False, num_workers=2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "5e7f0612",
   "metadata": {},
   "outputs": [],
   "source": [
    "# rfm import \n",
    "from rfm import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "5218e22f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaders provided\n",
      "Round 0 MSE:  0.05582288420917185\n",
      "Round 0 Acc:  0.6542\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████| 2500/2500 [00:50<00:00, 49.26it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Round 1 MSE:  0.048552343881175526\n",
      "Round 1 Acc:  0.693\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████| 2500/2500 [00:36<00:00, 69.14it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final MSE:  0.047642564293349135\n",
      " Final Acc:  0.6844\n"
     ]
    }
   ],
   "source": [
    "# run rfm\n",
    "M, _ = rfm(train_loader, test_loader, iters=2, loader=True, classif=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae02aed4",
   "metadata": {},
   "source": [
    "We have run three steps of RFM (the first iterate is the original laplace kernel), returning the M matrix of the final iterate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "7a730cd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize M matrix\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "id": "66b0f905",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsUAAAEWCAYAAACOieoEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJOklEQVR4nO3deXxU5b0/8M+ZfbJNdpJA2FUERVtFQa8LqEVraWur3FpbwfZSW4XWerUtbX+C2hZrba1ar167iLdV26rUWqugVnBDREtBQEG2QFhC9pkks885vz9oRsLM9zsJSyZmPu/Xa14vmCdne+bM8zw5ec7nGJZlWSAiIiIiymG2bO8AEREREVG2cVBMRERERDmPg2IiIiIiynkcFBMRERFRzuOgmIiIiIhyHgfFRERERJTzOCgmIiIiopzHQTERERER5TwOiomIiIgo53FQTB9JhmFg4cKFx3QbixcvhmEYqKurS753/vnn4/zzzz+m282mwX58RJTb+qPvoI8uDooHie4BnGEYeP3111PKLctCbW0tDMPApz71qSzsYXorV67EwoUL0d7enu1dISKiXli/fj0uv/xyjBgxAh6PB0OHDsVFF12E++67L9u71u9Gjhw5oPpUOjKObO8AHV0ejwePPfYY/uM//qPH+6+88gp2794Nt9udpT1Lb+XKlbj11lsxe/ZsFBcX93q5UCgEh6P/T98XXnih37dJRDRQrFy5ElOnTsXw4cMxZ84cVFVVob6+HqtWrcI999yDefPmZXsXiQ4bB8WDzCc/+Uk88cQTuPfee3sMGh977DGcdtppaG5uzuLeHRnTNBGNRuHxeODxeLKyDy6XKyvbJSIaCH784x/D5/Ph7bffTrmQ0djYmJ2dIjpKOH1ikLnyyivR0tKCF198MfleNBrFk08+iS9+8Yu9Xk/3n4RWrFiB008/HV6vFyeffDJWrFgBAFiyZAlOPvlkeDwenHbaafjXv/7VY/l3330Xs2fPxujRo+HxeFBVVYWvfOUraGlpSf7MwoULcfPNNwMARo0alZz+0T2H1zAMzJ07F48++igmTJgAt9uNpUuXJsu654WFQiGMGzcO48aNQygUSq6/tbUV1dXVOOuss5BIJNTj3bhxI6ZNmwav14thw4bhRz/6EUzTTPm5Q+fcRqNR3HLLLTjttNPg8/mQn5+Pc845B8uXL09ZtqWlBV/+8pdRVFSE4uJizJo1C+vWrYNhGFi8eHGPn3355ZdxzjnnID8/H8XFxfjMZz6D999/v8fPLFy4EIZhYOvWrckr7T6fD9dccw2CwWCPn3344Ycxbdo0VFZWwu12Y/z48XjggQfUOiEiOtS2bdswYcKEtH/Zq6ys7PH/3rY7R9rfzJ49GwUFBdi+fTumT5+O/Px81NTU4LbbboNlWRmPac+ePfjKV76CIUOGwO12Y8KECfjd737X+0o5SF1dHQzDwF133YX7778fo0ePRl5eHj7xiU+gvr4elmXh9ttvx7Bhw+D1evGZz3wGra2tPdbx17/+FZdeeilqamrgdrsxZswY3H777Wn7se5teL1enHHGGXjttdfS3hsSiUSwYMECjB07Fm63G7W1tfjOd76DSCRyWMc5WPFK8SAzcuRITJkyBY8//jguueQSAMDzzz8Pv9+PL3zhC7j33nt7va6tW7fii1/8Iq699lp86Utfwl133YUZM2bgwQcfxPe//31cd911AIBFixZh5syZ2Lx5M2y2A79nvfjii9i+fTuuueYaVFVVYePGjXjooYewceNGrFq1CoZh4HOf+xw++OADPP7447j77rtRXl4OAKioqEjuw8svv4w///nPmDt3LsrLyzFy5MiU/fR6vXjkkUdw9tln4wc/+AF+8YtfAACuv/56+P1+LF68GHa7XTzOhoYGTJ06FfF4HN/73veQn5+Phx56CF6vN2MdBQIB/OY3v8GVV16JOXPmoKOjA7/97W8xffp0rF69GqeeeiqAA1e5Z8yYgdWrV+Mb3/gGxo0bh7/+9a+YNWtWyjpfeuklXHLJJRg9ejQWLlyIUCiE++67D2effTbWrFmTUgczZ87EqFGjsGjRIqxZswa/+c1vUFlZiZ/+9KfJn3nggQcwYcIEfPrTn4bD4cDf/vY3XHfddTBNE9dff33G4yQiAoARI0bgzTffxIYNG3DSSSepP9uXdudI+hsASCQSuPjiizF58mTceeedWLp0KRYsWIB4PI7bbrtN3Mf9+/dj8uTJyYswFRUVeP755/HVr34VgUAAN9xww2HV06OPPopoNIp58+ahtbUVd955J2bOnIlp06ZhxYoV+O53v4utW7fivvvuw0033dRjEL548WIUFBTgxhtvREFBAV5++WXccsstCAQC+NnPftajfufOnYtzzjkH3/72t1FXV4fPfvazKCkpwbBhw5I/Z5omPv3pT+P111/H1772NZx44olYv3497r77bnzwwQd4+umnD+sYByWLBoWHH37YAmC9/fbb1q9+9SursLDQCgaDlmVZ1hVXXGFNnTrVsizLGjFihHXppZdmXN+IESMsANbKlSuT7y1btswCYHm9Xmvnzp3J9//3f//XAmAtX748+V73tg/2+OOPWwCsV199Nfnez372MwuAtWPHjpSfB2DZbDZr48aNacsWLFjQ47358+dbNpvNevXVV60nnnjCAmD98pe/zHisN9xwgwXAeuutt5LvNTY2Wj6fL2XfzjvvPOu8885L/j8ej1uRSKTH+tra2qwhQ4ZYX/nKV5LvPfXUUyn7k0gkrGnTplkArIcffjj5/qmnnmpVVlZaLS0tyffWrVtn2Ww26+qrr06+t2DBAgtAj+1YlmVddtllVllZWY/30n0e06dPt0aPHt3jvUOPj4joYC+88IJlt9stu91uTZkyxfrOd75jLVu2zIpGoyk/29t250j7m1mzZlkArHnz5iXfM03TuvTSSy2Xy2U1NTUl3z+07/jqV79qVVdXW83NzT326Qtf+ILl8/nSHsOh+35wn7pjxw4LgFVRUWG1t7cn358/f74FwDrllFOsWCyWfP/KK6+0XC6XFQ6Hk++l2+a1115r5eXlJX8uEolYZWVl1qRJk3qsb/HixRaAHu3473//e8tms1mvvfZaj3U++OCDFgDrjTfeUI8xl3D6xCA0c+ZMhEIhPPvss+jo6MCzzz7bp6kT3caPH48pU6Yk/3/mmWcCAKZNm4bhw4envL99+/bkewdfZQ2Hw2hubsbkyZMBAGvWrOn1Ppx33nkYP358r3524cKFmDBhAmbNmoXrrrsO5513Hr75zW9mXO65557D5MmTccYZZyTfq6iowFVXXZVxWbvdnpxnbJomWltbEY/Hcfrpp/c4zqVLl8LpdGLOnDnJ92w2W8rVkn379mHt2rWYPXs2SktLk+9PnDgRF110EZ577rmUffj617/e4//nnHMOWlpaEAgEku8d/Hn4/X40NzfjvPPOw/bt2+H3+zMeJxERAFx00UV488038elPfxrr1q3DnXfeienTp2Po0KF45plnevxsX9qdI+lvus2dOzf57+4rv9FoFC+99FLaY7EsC0899RRmzJgBy7LQ3NycfE2fPh1+v79P/dXBrrjiCvh8vpT9/tKXvtTjfp8zzzwT0WgUe/bsSb53cL11dHSgubkZ55xzDoLBIDZt2gQAeOedd9DS0oI5c+b0WN9VV12FkpKSHvvyxBNP4MQTT8S4ceN6HOO0adMAIO10v1zF6RODUEVFBS688EI89thjCAaDSCQSuPzyy/u8noMbIgDJL3htbW3a99va2pLvtba24tZbb8Uf//jHlJsv+jIIGzVqVK9/1uVy4Xe/+x0mTZoEj8eDhx9+GIZhZFxu586dyQbrYCeccEKvtvvII4/g5z//OTZt2oRYLJZ233fu3Inq6mrk5eX1WHbs2LEp+yJt+8QTT8SyZcvQ1dWF/Pz85PuHfk7dDWJbWxuKiooAAG+88QYWLFiAN998M2W+sd/v79F4ExFpJk2ahCVLliAajWLdunX4y1/+grvvvhuXX3451q5dm7yQ0Zd250j6G+DARYbRo0f3eO/4448HgB5Z8wdrampCe3s7HnroITz00ENpf+Zwbx48kuPZuHEjfvjDH+Lll1/ucXED+LD/7O4rDu1DHA5HyhS7LVu24P333+8xNfFgvEHyQxwUD1Jf/OIXMWfOHDQ0NOCSSy7pU9xZN2kervS+ddANDTNnzsTKlStx880349RTT0VBQQFM08TFF1+c9gY2SW/m9R5s2bJlAA5cnd6yZUufBtWH4w9/+ANmz56Nz372s7j55ptRWVkJu92ORYsWYdu2bcd0290yfR7btm3DBRdcgHHjxuEXv/gFamtr4XK58Nxzz+Huu+/u0+dBRNTN5XJh0qRJmDRpEo4//nhcc801eOKJJ7BgwYI+tztH0t8cru59+NKXvpT2/g7gwF/pDsfhHk97ezvOO+88FBUV4bbbbsOYMWPg8XiwZs0afPe73z2s9to0TZx88snJ+20OdehAPZdxUDxIXXbZZbj22muxatUq/OlPf+rXbbe1teEf//gHbr31Vtxyyy3J97ds2ZLys725kttb7777Lm677TZcc801WLt2Lf7rv/4L69evz3gVdMSIEWn3bfPmzRm3+eSTT2L06NFYsmRJj2NZsGBByjaWL1+OYDDY42rx1q1bU35O2vamTZtQXl7e4ypxb/ztb39DJBLBM8880+PqBf9kRkRHy+mnnw7gwBQwoP/bHdM0sX379uTVYQD44IMPACDtDdrAgb+qFhYWIpFI4MILLzwm+9VXK1asQEtLC5YsWYJzzz03+f6OHTt6/Fx3X7F161ZMnTo1+X48HkddXV2PwfyYMWOwbt06XHDBBUe1zx2MOKd4kCooKMADDzyAhQsXYsaMGf267e7fhA/9Tf6Xv/xlys92D/CO9Il2sVgMs2fPRk1NDe655x4sXrwY+/fvx7e//e2My37yk5/EqlWrsHr16uR7TU1NePTRRzMum+5Y33rrLbz55ps9fm769OmIxWL49a9/nXzPNE3cf//9PX6uuroap556Kh555JEedbJhwwa88MIL+OQnP5lxn3qzj36/Hw8//HCf10VEuW358uVpr9J23+/QPfUrG+3Or371q+S/LcvCr371KzidTlxwwQVpf95ut+Pzn/88nnrqKWzYsCGlvKmp6ZjtqyRdvUWjUfzP//xPj587/fTTUVZWhl//+teIx+PJ9x999NGUqSUzZ87Enj17evQ/3UKhELq6uo7mIXyk8UrxICb9OehYKyoqwrnnnos777wTsVgMQ4cOxQsvvJDymy4AnHbaaQCAH/zgB/jCF74Ap9OJGTNm9Plq6I9+9COsXbsW//jHP1BYWIiJEyfilltuwQ9/+ENcfvnl6mDyO9/5Dn7/+9/j4osvxre+9a1kJNuIESPw7rvvqtv91Kc+hSVLluCyyy7DpZdeih07duDBBx/E+PHj0dnZmfy5z372szjjjDPw3//939i6dSvGjRuHZ555JplPefBv7z/72c9wySWXYMqUKfjqV7+ajGTz+XzJbOa++MQnPgGXy4UZM2bg2muvRWdnJ37961+jsrIyeVWHiKg35s2bh2AwiMsuuwzjxo1DNBrFypUr8ac//QkjR47ENddcA6D/2x2Px4OlS5di1qxZOPPMM/H888/j73//O77//e+Lc2kB4I477sDy5ctx5plnYs6cORg/fjxaW1uxZs0avPTSSykZwsfaWWedhZKSEsyaNQvf/OY3YRgGfv/736f8IuJyubBw4ULMmzcP06ZNw8yZM1FXV4fFixdjzJgxPfqUL3/5y/jzn/+Mr3/961i+fDnOPvtsJBIJbNq0CX/+85+xbNmy5JX+XMcrxXRMPPbYY5g+fTruv/9+zJ8/H06nE88//3zKz02aNAm333471q1bh9mzZ+PKK6/s82/na9aswU9+8hPMnTu3x5+Rvve972HSpEmYM2eOeiW6uroay5cvx8SJE3HHHXfgl7/8Ja6++mp861vfyrjt2bNn4yc/+QnWrVuHb37zm1i2bBn+8Ic/pDQwdrsdf//73/Gf//mfeOSRR/CDH/wANTU1ySvFBz+h78ILL8TSpUtRVlaGW265BXfddRcmT56MN95447DmSJ9wwgl48sknYRgGbrrpJjz44IP42te+1qvjIyI62F133YWpU6fiueeew4033ogbb7wRq1evxnXXXYe33noref9Kf7c7drsdS5cuRUNDA26++Wa8/fbbWLBgAW6//XZ1uSFDhmD16tW45pprsGTJEsydOxf33HMPWltbe2S995eysjI8++yzqK6uxg9/+EPcdddduOiii3DnnXem/OzcuXNx7733YteuXbjpppvw2muv4ZlnnkFxcXGPPsVms+Hpp5/GHXfcgfXr1+Omm27Crbfeirfffhvf+ta3ekw5yXWGdTRmqxPRYXn66adx2WWX4fXXX8fZZ5+d7d0hIvrImT17Np588skef53LVaZpoqKiAp/73OfSTpcgHa8UE/WTgx9BDRx4AtN9992HoqIifPzjH8/SXhER0UdROBxOmVbxf//3f2htbU15zDP1DucUE/WTefPmIRQKYcqUKYhEIliyZAlWrlyJn/zkJ32OniMioty2atUqfPvb38YVV1yBsrIyrFmzBr/97W9x0kkn4Yorrsj27n0kcVBM1E+mTZuGn//853j22WcRDocxduxY3HfffT2ewkRERNQbI0eORG1tLe699160traitLQUV199Ne64447kk1apbzinmIiIiIhyHucUExEREVHO46CYiIiIiHIe5xRnYJom9u7di8LCQj4ekWgAsSwLHR0dqKmpgc3G3+8pu9hXEA1MfekrOCjOYO/evaitrc32bhCRoL6+HsOGDcv2blCOY19BNLD1pq/goDiDwsJCAMCTv/sy8vNS7+Z88Cm7uGxpkVwGAKccL98dumxVSCzTbG7Qr1BcNS1PLGtojohlXrd8qowfJd+ruWaTqe6P2yXX0aih8rE893pQLHtrj7yvHTH9vtLzhyfEsqEV8r6ur5OPM5bQP5OJI+Sy8hL5HFm9Sf68RlXK2/R36XWwt00uK82Xl23skH8DL/Lo50FZkby/sXj6bcZjUSx//pHkd5Qom7rPw+Ipt8FweFLKE44Ccdm4ofcVlsspF7rcYlF5VZVYVuArUbfpdReJZQ6n3MYaNrmNsNtj8vZcGa6uW3GxSL34p+yPdkU/AX1/Egn5M4vKu4qEdQR/1bLLK3Y45PXa7cr5JXd5iMf0OrDUZl1eNm7Kx2G36/WTMOVzyIb0xxkPd+Cft36sV30FB8UZdH9p8vNcaQfFTqd8srlcevV6vfKAx+lUzlSF3aGfxB6P3IBqCS5uZVDs9cqNjsuVYVDsluvP65GPxeGUv1SGXelATH1AqK1X+zztDvnzMjP8KdWldAZut/yhOBzysThdcsPijOqfiV05bR1OpcNTGmWHU9+m06nUkaF/ZvxTNQ0E3eeh4fDAcKTmjqd778Nl9UExnErj7JTbdJsrXyyzu/UBgt0jlzucchtr2OTvujYodrgPf1CsjqMOc1BsZBgUG8qgWBnzwTCPYFDs6N9BMZR1AoCV0NpmZVllYJtpUGwcxqA4uWwv+goOig8RiUQQiXx4BS4QCGRxb4iIaCBiX0E0+PDulEMsWrQIPp8v+eIcMSIiOhT7CqLBh4PiQ8yfPx9+vz/5qq+vz/YuERHRAMO+gmjw4fSJQ7jdbrjd8hwtIiIi9hVEgw8Hxb307Ot5aRvA17fJk75PrZHLAKCqPPUO5W7Kzb244Up5sv2Pf6dPNF/8QpdYVpEv3yCR746KZU1t8jZb/PoNVg3tclnClDuc+lb5jxxDvPKdA3kZJvGXF8nlhXnycZ44TD7O6nK943y/Tk4aaesIi2XDK+SbBsLKzXRDSvVzZMt+edl/bJP3B5FWsaimXP/T8m6/fMPG8RXpP8+4cjMLUbYk4ISB1BvRXPk+cRmPR056AABLuUG6sExOkfAUyDfLebzyjX8A4HTK5U7tBqyEnAzkMOS2JdOvF27l/mnLlJN4rLhcZlpyv2Zod8sBsCspEk5laGVayoHYlDIA8ZjSnyr74zDkzzJhyPsa1+MlYBpKHSiDmETk8IIEACBhKTdHCjfSWRlumjwYp08QERERUc7joJiIiIiIch4HxURERESU8zgoJiIiIqKcx0ExEREREeU8DoqJiIiIKOcxku0IDcmXo0U27NOjr/yhTrHMoTyvPRyVP7ZASHsWOeBUfg0aWiYX7mmRo1k27pTLxlTpUSj5bnnZF9fIcTlFXvk4bYZcZpcT6QAALqe8v20BOaJne4O8zbYOOXINAJwOeZuxuLxeLeqt1CdH+2zcplfC0GL5M+mI5Ytley25fva2N6vbhE0+luMr9OgoooGktLwYNmdeyvu2omJxGSNPjmsDAJtXjlbzVVSIZXElViyW0GOx7Kbc/ubb5T7IpbRnTshtS2GGTDaPXT4Wl12JenPK/bDLJbctNpseR2Yp0WAJU+5L4wm5bQ4qUZoA4O+UP5P2oFwWicptfsKReq4m2eTYWACwDPlDCyfkSFoLcr9mJfQ6UBZFQuj7E9pCh+CVYiIiIiLKeRwUExEREVHO46CYiIiIiHIeB8VERERElPM4KCYiIiKinMdBMRERERHlPA6KiYiIiCjnZSWnOBwOIxqVM/UAwOVywePRM/L60/Y9ETidqVl3LiWKeHdY/53D1SGXRZUIyWt/oazTrufxfWykXBZXttnYKR+oxyFvc8f+DPszVs5sXP2BnEsZliMQ4bLL+ZEn16q7g3e2yftbVSTnJwaj8jYnjpaPEQA6uuSKLymSl32/LiyWrd8n111Fnv6ZtATlz7o1puROu4rlMj2uGogFxSKHsDtWhjhLomyoHV0Ghzs1z7vTlPuzuNOlrtPIk5d1F8jL2kz5i2fF9C+QR8np9brl9sULuV0qcMrLlXr1XP/KYjkjvaRQHspUlMrLlZbIdefJkJusxTyHlMzgaET+TPa3ys8uAIBGv5wpvHNfm1hW1ySv1690pqZb7ytMqXEGYCqLKqcl7Jmu1VpyxVuWsFHp/TT6fVAcDocxakQRGhqVUQ2Aqqoq7NixY0ANjImIiIhocOr3QXE0GkVDYwy73jkVRYXpf8sIdCQw/PS1iEajHBQTERER0TGXtcc8F+UbKMoXrqFr19aJiIiIiI6yrA2KLUue5tGH6R9EREREREcsi4NiGywr/YRqcbI0EREREdExkLVBccK0I2Gm33yCd5UTERERUT8aoFeKB158cnWZHS5XanW1d8kpGiPy9NF9MC7PnR5ZLMeOrNojx73UFHjVbb72gXwVPm4d3lzuApd8nOVF+rLRuLw/44bK+9PQKm+zS4m8icT0v0IUKjE8/2qQvy5Tx8qf1/e/skfd5mNLx4hlf3tDjtLJU1KcTqiQo486QvrnPKpUPpaCTrn+Clxy2e4OPW6pISp/nm1d6fc3rnx/iLLltI+PgDuvMOX9HU0RcZmGgN4uRW1KuSH3QY40MaLdTCXaCgBsSkwiIvI28/Pk/ntYmdw/jaiQo9MAYGytT96m0hYWp34USV5lOSXZ80C5UmaVyitWqg4jKkrUbQbktDtsr1LqZ8tesWxzg5wN2xqX+xEAiBpyeVyJVktAOS8Nva9Q8z2FMYzZh7FN1gbFpmmHaaY/eFMLuCMiIiIiOsp4pZiIiIiIcl6fRp8PPPAAJk6ciKKiIhQVFWHKlCl4/vnnD2vDpmWoLyIiIiKi/tKnQfGwYcNwxx134J///CfeeecdTJs2DZ/5zGewcePGPm/YtOzJKRQpLyvTnBIiIiIioqOnT9MnZsyY0eP/P/7xj/HAAw9g1apVmDBhQp82bFkGLOGKsPQ+EREREdGxcNhzihOJBJ544gl0dXVhypQpfV6ec4qJiIiIaKDo86B4/fr1mDJlCsLhMAoKCvCXv/wF48ePF38+EokgEvkwiiYQCAAATNOAaaYf/JoD8DHPwbCJuJkaG+Vxysu0RfXBfblbjqFa1ah8NE65bK8S2QIApyhxW9U+OfWjrkU+luagXHbiMD1JxOuWl+3okutnSIm8XGuHvFymv0K4HPL+zjxdXtbllCN4fvSbGnWbDS1y7NrwCnmb5cXyyffGRjn3pyum14FXiVYr8shlk06Q9yfxnpJDBKAyzyOWffy49OuNRk288w91tUT97pyzRiK/MDWLsmx7SFzm3e2t6jr3tsvfu0hCXq+hJDmZkS51m2ZMLi9UssxGVsrRYOOHl4plo4bIbQAADK2Qy7T4NK0XVnvoNP19z23KdWtTYsXcyuzQNKdND2VKeVGRHHdXUDhC3mZdm1i2aY8SywdgT5sS/anEbNog9xWWoY+bTCXOzabFtfVSny/JnnDCCVi7di3eeustfOMb38CsWbPw3nvviT+/aNEi+Hy+5Ku2thbAgcGJdJMdp08QERERUX/q86DY5XJh7NixOO2007Bo0SKccsopuOeee8Sfnz9/Pvx+f/JVX18P4N832ikvIiIiIqL+csQ5xaZp9pgecSi32w23O/UxYQP1RjtpugcREVE39hVEg0+frhTPnz8fr776Kurq6rB+/XrMnz8fK1aswFVXXdXnDXcPiqVXtkjTPYiIiLqxryAafPo0KG5sbMTVV1+NE044ARdccAHefvttLFu2DBdddFGfN2yaNvWVLdJ0DyIiom7sK4gGnz5Nn/jtb3971DY8UKdPSNM9iIiIurGvIBp8jnhO8eFKmDYkhCvC0vvZ1NxuwuFMjTOzK7s6slCOKwGAArccLeKyy2Xjh8plwbAegdYelH/hcLvksvaIfKDVBXLM27/q9F9wxg+Vo7oSSiJOZakc6VJRIm+zfn9U3Z+WLvk4P9gt7+sFpxeKZa+t7VC32d4l729hnrzcB7vkY1m7a59Yds7x1er+rNkv1+34UvmcfnebvD/+sH7zbKEST7hhR/r1xmN6zBtRNtSUAgVporPaQ3Jk1k756woA2Nckx6NZCSXeSlmnPSHfCwQAxXlyWzh6WLFYdupxVfJylXJbV5Gv7g6UplAN4tJGE2rvlGEYosd/KbGgyoojUbkvBQC3Q25Hy5VEO+dweajncMtZd+GYHNcGAC2t+8WyZr98frkKy8SyTI+pMJUfsIQz3rL0cdHBsjYoHqhXiomIiIgo93BQTEREREQ5L2uDYtMyxCfXmRwUExEREVE/yt6VYhiwhDk50vtERERERMcCp08QERERUc7L3qDYNMQ8YkuYVkFEREREdCxkcU7xgZdURkRERETUX7I3KFaeXJfNJ9pJtjTZYEuTETh5rDyCbwroo/vtbXL1R5SM3uEhOSt2Z4tedyPL5RV7XPKyJ1XJ2yzOl6/smxnyAavK5IzNhhY5h7Zun1w2boQc2Dh6qBLmCKChPSyWaYeyuS4olv1zj/41O3uUnE35+ha5bmt98jovPbVGLJt8UoG6Px2vdoplQ5RtNgbkfS1wKSc0ACXaFI1+dVGiAcWL9Jm6IyrlZTYXutR1NufL7VZDi5I3bMjfyaJ8/cEj1T552fGjh4hlw5Us4kql6XHpEb3wKFHnCWVZmx6RLrIyBOYmEnKHYLPJdaD9HTzPpe9s1JLbUYeSf5ynbHSEHBmMpmqlwQewZcsescxtyX203ZQ/sFBcz892uOS8byjr7S3OKSYiIiKinJfFQbF85a0PDx8hIiIiIjpiWZxTbBMf16c9xo+IiIiI6Gjj9AkiIiIiynlZvNFOeaIdI9mIiIiIqB9xTjERERER5TxOn+gluwGkS1nZ3ypHpNQF9OqNKlfECxzybwbFBfKc60hMj77a0ihHvgyNyHEm0bi8r3VKDNyEofr+7G9VotVG5YtlEWVfK8vkel+/RY5OA4DjauRjefk9ebmqEvnzmjxcr4P/vEiORjJekuNpAl3yOpsDcpka4QRgd6d8jhR55Wi+kZVy3W3ao/+m+9YOedkFV6WPowqFHFj5grpaon7n/vfrUFosVkWBHE0JAEUOeeFWU25DtW9dvl3fZqVPjr6q9Mnf12K52YZXSRxTEkEBAIZyMGnSUpO01ldJVUNUbuoA6A8ZsykbdR9mRBwA2Cx5p5RTBC7In3WBU16wIk//UGp88ofdEZJjBtvjcgXZjMOvIOkc0c6dQ2VtUJywDCSEwa/0PhERERHRsZC1QTGUK8XgoJiIiIiI+lEWI9mUG+04KCYiIiKifsQ5xURERESU85g+QUREREQ5L7vTJ4Qrwpw+QURERET9idMneml0hQmnMzUG7PkdSqxNXI++KvHIkSXnHqfHeEmaO/UIle2d8mX4oT55uQKPvNzwCrmsulw+RgBoapPjhJpa5frzd8qRbDd8cbtY9v+2D1P3JxKTj+XSj8lRMU6HXO+F+fpn+cCSsFjmVk6vt/bKhV1y9SCmRDgBQHNUibRrlMscNjkuqDWkx+wM88l1dM3nX0n7fqAjjnk3q6sl6neef78OVaj0tjUlBeo6t7nlKEmH0s8klGwwt9JmAUBpidy+lJTIy7mU41Tj0bR8NAAJZWEL8pihuU2un8b2DrEsrHffsEOuvzyv3O9VlspRd+Vlejtps8nbtCA3+jZLbvNdhryvvjx1d1BSIMeJFitxbh1++bO2G/r4TztLLGGagfR+Opw+QUREREQ5j495JiIiIqKcx+kTRERERJTzeKMdEREREeU8zikmIiIiopzHQTERERER5bws3mgn31BnHl4a2TH13j4bbI7UuJTJlXIMVXtIr95NHfI0kSa/XAk7WuTYlk9NUjK8AFRtl+O/8tzy/qQ59KRt++XfYkKRqLo/eR55m9UV6YKNDnA75fVuqqsRy2Lyx3Vgf9xyjExBnvx5doXkFQ+tlI8DAPa3yssGgnLd/sdwOWZnX7tcr/ku/bfOKVXy/hhKIE5bUNmmU9/mLr98gk3/+ulp34/HogD+qa6XqL/Zkb5jdSmzAsu0vDYARW75+1GYJ8diRSw5pqsgX2+XiosLxTKPsuhhX9OyZ5g2qSTI+QNy2Qe79otlu3Y3iWXBkJJrCSAWk8vLivPFshNHV4plJaVV6jbthnwe2JRINrsSWaeFwLn0hDg4bMqnrQzkHDZ5fwxlXwEAlrxeacm+nJNZGxQf6F6lQ+CcYiIiIiLqP1m+UiyXERERERH1F0ayHSISiSAS+fBRNoGA8ncZIiLKSewriAYf/TmPx5Bp6a9sWbRoEXw+X/JVW1ubvZ0hIqIBiX0F0eCTxUGxob6yZf78+fD7/clXfX191vaFiIgGJvYVRINP9qZPQIlk69c96cntdsPtlu/mJSIiYl9BNPhwTnEvhS0DtjQRclqm8sdH6eu0tssRKg1++SJ+QtlmOKLHyMSV4q6wvOLTT5QjZhrag2LZqScUqPvT2CpHxGma2+U4svd2yPsaiXap69X+SrFrf0gsqx0if5Wa2/RYuqoyednOkHycLUqk37vt8jqnuPVcuoRyo+uoIfJ5uWqbvD87Q/ofpYa55Y1KUYExJdaHKFsMpM9P0tKttMhLQE8rM5TYNe2u9VhC7ytCMbndCptyG+tRvupaK6AHxAGGUgctnfKx7GuR2/z97XL/EwpGxDIA6OrsFMuCQbmsotQrliXMDJFsyolgg0spUwYNSlEswxXKuDLXNa6cezZD/mXSlmH8p06vlQr7MCeX6RNERERElPOyOH1CzimW84uJiIiIiI4+XikmIiIiopyXxTnFyo122bzTjoiIiIhyDm+0IyIiIqKcl73pE8pDOrL58A4iIiIiyj1ZGxRDmT6R1aBiIiIiIso5WbxSLD+5LptPtJMMzTdhd6RmIXZF5X11OfXjqCqSsxXfanD2fucOsmGHnGsL6FmP5T65cG+TnOdYVSyv09+hZ/R2BpV8SWWbNZVyzuHWXXIdFOTpeblF+fJXYneTnFtZUih/XmOG6SGkr6yRc56HVcjL5nnlso535SziN/dn+Np37BKLgtGhYlmhS/ltVj5EAEBbTP5c4kKealwL7Cb6COnq0NvtWExuJ50OOfc2rnxFEqacawsA7QH5jvcWv7KgslpT6dYSGbr9pv3ywWzb2SiW1e9rF8vaAnKecCymfyYJpdxmlz8TLWs4wybVrGtDDbNWtqkEG4T17hsh5QRLGHKbbikZ8/Y0z4PotaMw+4A32hERERFRzuOgmIiIiIhyXhZzig2YwmVy6X0iIiIiomMhi0+0k++n44ViIiIiIupPnD5BRERERDmP0yeIiIiIKOf1aVC8aNEiLFmyBJs2bYLX68VZZ52Fn/70pzjhhBP6vOGP2vSJsya44HanxoD9caUcfbV5txyjAwAbm+V8mnKnnJNSUySXBZWIOADY2ymHunz2XDlGZvWGLnW9kkCXXD8AEIkpkS6mXH9alNsZJ+WJZfua9VN++245O8yjxAlt3B4Sy04cla9usyhf/kz2t8r1p9VdtU85f4qVDB4A21uGiWUbGlrlBT3lcllEWQ5Al6tELNvWmD7aJxHX4/WIskHq27R4tJYWORoMAMJyOiU6OuUVh0z5O9KlrRSAv7VNLOv0y1ld5R45k83nlsuclh5d2dUh58C1d7SLZVGb3HD7yuQ2y2HT+2+vQ673kVVlYtnwoZVimZJUlpHWqmvJaiHlvGzu1Pvv5k653wtqG7XJkarqgQB6vtpRmH7Qp17llVdewfXXX49Vq1bhxRdfRCwWwyc+8Ql0dfV9wGRagGkKr4E4KiYiIiKiQatPV4qXLl3a4/+LFy9GZWUl/vnPf+Lcc8/t04Y5p5iIiIiIBoojmlPs9x/4c0Zpaan4M5FIBJHIh08CCwQCADgoJiIiIqKB47An5ZmmiRtuuAFnn302TjrpJPHnFi1aBJ/Pl3zV1tYeWP7fj3mWXkRERERE/eWwB8XXX389NmzYgD/+8Y/qz82fPx9+vz/5qq+vP1BgZXgREREREfWTw5o+MXfuXDz77LN49dVXMWyYfKc6ALjd7rSpDZYl31DH6RNERERE1J/6NCi2LAvz5s3DX/7yF6xYsQKjRo067A1/1OYUB0NxxBOpkTEji+XYlv0d+oV4nxK7Vp53eLFrBW698nxueb37W+RjSSiRIF1huay8WK+DkcVyRM/40fJxvvx2RCzL88j74yvU92dImRzfE4rGxLKyInm9b23QsmmAcET+TLTYtbomeZ3Dy+S627Jfn55U4JS3mZ8vRxh1de4Xyy4eL993AABLN+wRyzZ3CJ+1KX8eRNmS+PfrUHEl4autM0MbEZPbF1/JELHMmdDaMz0GTtunut1yFFe7S47xKnDLZWZEj0CLBANiWWGx3KaNHDlCLKutLhbLiuSEUgBAgVP+TMoL5Hi5YmW9tgxjH0PpvrTa01LOlFQ1tAaUQgDtfjl5zN8lfyaeokKxzDD0/smAXG6JR5op5+1DfRoUX3/99Xjsscfw17/+FYWFhWhoaAAA+Hw+eL0ZzqBDfNQGxUREREQ0ePVpTvEDDzwAv9+P888/H9XV1cnXn/70pz5v2LT0FxERERFRf+nz9ImjhVeKiYiIiGigOKKc4iPBQTERERERDRRZGxRr0yQ4fYKIiIiI+hOvFBMRERFRzsvelWLTgGmmj9aQ3s+mrpCFaCI11mN1k1yFY/P1GJBRpXKISqVProM9rfJvDeEMKVVuOSkGja1yzNmkCQVi2cp1crRPUb5+ihUVyOUNLfJynSG5bjdslWN/OkN67M/eZnlZLVIpz6tUbIbf8k46Lk8sGz5E/kB//1xYLNvUIO9PoRLLBwCxhHzuVXvlZessORqqTon7A4Ain5x3HpAWTcjnK1G2xACkCzOr2ycvs98vtzsA0BqWvz9hS05+0taaMPW22eUskbcZle/Rb43L7V1rp7xHdpt+378Vl9ulPIfc3o0/uUIsG1Ikby9Deic8yjDFrZSZUbl+nE597BNNKH2JQ142piT+vbe+XizbtUuPZIuE5P4pHpH3Jx6Xz4NoXK94p9sjljkc6cvsRu/jO7N3pRjyg+t4oZiIiIiI+hOnTxARERFRzsvijXYWTOGOOjOLo+JIJIJI5MM/ywYC8lN0iIgoN7GvIBp8+vTwjqPJyvDKlkWLFsHn8yVftbW1WdwbIiIaiNhXEA0+2RsUm4ApvKzeP6b6qJs/fz78fn/yVV8vT0InIqLcxL6CaPDJ2vSJgXqnndvthtvtzt4OEBHRgMe+gmjwYfpEL7ldgMuVGjFyZoUcLVLg0Y/kld1Osay6Vb5c7nPL693gV6LBAHxytBxNsqdZXm5fS4dYVlkib7M1oEehhKNy1FBCiZ/Rpp2Ho3LdNbbp0WDFBfIfT97fLW90rPLEGW1/ACDQJS+7u1E+R9q75Ei2fJe8zbPHu9T92bBDjjpr6pDrZ2Kp/F1Y05yhqVGSiD41Jv16Y7EElr2rr5aov3WFAVuar1iLX257Wjv0drK1Q87UMtz5YpmlfLGsDMmnNrv8nTUNuSyq9OAJS24jjAxRrNGQXAdlcbmdVBLQYFM2abPp/bdTqVvtT/DWMfr7vJJyho6AXOhvl+fCB9qD6jYPnlN/qLw8OdIvpnTgmc5LjWWkX1h6P50s5hQfeEllRERERET9hVeKiYiIiCjnZTGn2IIlXEKX3iciIiIiOhY4fYKIiIiIcl720ic4gYKIiIiIBogsXilWnmin3MlPRERERHS0ZXFOsRyrxSnFRERERNSfmD7RS1v3WXA4Uic7R5RswMIM0XjHF8q5lR1ReeGYFrUbl3MDAcDpkDOFffnyNusa5U/FG5QngUdj+qd5yvFyxmZnUM7uHFkjh+ZXlcmZlXubO9X92bFf3t+Pj5EDJj82Ts7+3bxTz0ZubZfzN599U66DKSfK29zXLC+3erO8PQBo7pLPkcvP9ohl//eKcu5F5SxMAKgsKBTL3qlPvz9mPGsP5CQSNbcA6SJ1m1rkXPFOPQ4WYS1sF0qGMeS2x4B+847NJndu8YSyXiUT1rK0tlDfnxjkNs20y22htlZTaUJMLTwdUPZGH8MklEf2Wgm9TbM75H3SBnN5eXJpdVWpWLa3U68DZ4t84kYdci3YtMcp6N2lGsRgCmVmH25Uy96g2LRgCdMkpPeJiIiIiI4FXikmIiIiopzHG+2IiIiIKOfxRjsiIiIiynlZzSm2OIGCiIiIiAaA7F4pFm4I5JViIiIiIupPvNGul44faoPLlZojsrlejq2pa9FyR4DxNXJMyLb98nKnjJTLOjbL0TQA8O5uuXbPOk5e7rTj5JizUcPkbQb0BDS0BeQ4oTHD5Ni1f22So2CaWuV1RjJExAXCciSO2ymXxeJydE1nUM+YCYbl8+DCj8kRaE5lf15YK2/z7OPV3cHqd5Q4t/fk9doMuTkZVlSgbvPjw+X17mtL/34ibqJRXStR/9uxuwN5aeIt9zTKbVZQT9KECbkvcdrltsdmyN8rp6G3hU6bvGw8JjfslpLiZdjl49CitgDAbpPbyVhU3tdWv7zOErl5hTNP3R1oiW125VDcSixqpuCwww2hLFSO5YTjq8WyDpsclQkAW5VzOqREtSqpfcg0AtQj2dJv86MRyWZZ4sFl+nIQERERER1NWUyfOPCSyoiIiIiI+gunTxARERFRzmNOMRERERHlvMOdt01ERERENGjw4R1ERERElPM4faKX3tlqwe5IvQOwPSJX4alD5bg2APC65Qv12zrl2JaJlhw/c0qNvs3WTjkLZfU2uWxYsZwZJD+EBbDb1OwVdIWUKK5muez80+W4tj++IMfEjK6Ro+UAIJ6Q48hiCfnz6gzJxzlulJL7A6AoT76ztLFd3ubSNzvEsu2d8nmZtyPDnawhOQ9wQ2OtWOZS/u7k1NMJ0RyQz6GOSPoVJ5QYPKJs2byjGW5vanu5q1Vum/1KFCQAxJUu0RaXIyiddrkNLS7Uu39vmv6uWyQsbzOWkMtgl7cZ0Q4SAAx5va3NIbHs3fXbxbJAkxw5VlUi9zEAUFUmL5vnko+lolD+rPXeCUowHxCNyn20zSUfixbXVl6qR2nm58mfp0NotwEgrozxjqRVl7bYlykRvNGOiIiIiHIeI9mIiIiIKOfxSjERERER5bysDYp5px0RERERDRTZmz5hHXhJZURERERE/SWLV4qVC8IcFBMRERFRP2JOcS+5HRbsjtQdu/xUOSQlYepBIHa7XB5OyFE6j6+Tl/tYuV55WlSMZcnrLS+SyyqU6Jqde+WoHABwK7F0FSVy3QbD8nEMHyKf1qOG6vFoHUG53kuLXGJZi19errxYzyN7+R25jlr8Slxbh1x344vl5coL9DtZS0rk2DWPTa730WVyHby6S99m3JTDiEKJ9CE9pvA+UTbtagzA6Uk93/0xue1J2PX4L0eavidZZpO/dyUFchtx/PASdZs1pfJ3MhHrFMviphxrmTDktrkjpLcRW7ftEstam9vEss2b6uTl9sn1XlqkB6TV1sj1V1wgH+dJo4aKZUOK5T4GAFxKV+J2ZAp0S8+pNKN5+u6gwCPXn0PpKwxT3qiVYQBoKZG08piy94PKLE+fEHKKB+CgmIiIiIgGL14pJiIiIqKcxxvtDhGJRBCJfPhkmEAgkL2dISKiAYl9BdHg05en3x19lvDKokWLFsHn8yVftbXy/EoiIspN7CuIBp+sDYq7p09Ir2yZP38+/H5/8lVfX5+9nSEiogGJfQXR4MPpE4dwu91wu/U7gYmIKLexryAafHijHRERERHlvCxeKbaUSLaBNypuCRqwOVJnm3ywK5Lmpw9YWa/nBtbkRcWyMyrlOtgbkGe9/KtZz8SdWCpn/BWnydbsVl0hXxFpbJXrIBzVsyeLlDzH2iFy3uXDfwuLZWeeVCiWtbTLxw8AZ5+SJ5a9t13+vLpC8nqdDq+6zZJCuQ7er5frYOwQ+RzZ0STnQO5u18+RYiXLuj4kn3vFXfJyo4r0beYpOaz5QnuQiFtoVddK1P+CCQecidTvtOmQc4rtLj0/3W3GxbJypf0YWSGXnXpcgbrN45Up0k7j8K6Qa61vmx5pD3u0SyzbFJLb5vY2ebl2U+67goGguj9dnfIO+/Ll5YpcSjZyYY26TW3Cq8N2eLNhtZa5MMPH7MuXD9Rpk+vHsA4/Y17LHE4IRYmPQk4xrxQTERER0UDBQTERERER5bw+X29/9dVXMWPGDNTU1MAwDDz99NOHteHuG+2kFxERERFRf+nzoLirqwunnHIK7r///iPa8ECNZCMiIiKi3NPn6ROXXHIJLrnkkiPesPacDo6JiYiIiKg/ZW9OsQmYQjCBpQcWEBEREREdVcd8UCw9H/6jdqW4PN+C3ZE6Wq9rkWegeGz6kdT45NF/RE7gwSglVq0yqkedvNsqf+RTquT4r7q9cnRNebES6pIheSU/TzsF5UqoLJG3OXGsHM+TMPUdWrdF/jx3N8r10xSQ1xkM69E+1RUuscxQdlfbZldMXrBQiVwDgMagXAd5yoSrPGW9lYX6b7r17fKKq4Rl4zETW9S1EvU/u8OA3ZH6/XOlifT8sFBfZywstz1OpZ9x2+TvXUGGbfqUci2pS5uTKR8FYM/Q8deUyrF1rUVyNFgiLNeBR4lHMxN6ux2JyP1wIKFEvYXlMpueXImEkmmXiMr9nmWT+wPTLkfHOjOMEPO9coSpyyn331Zc+7APP67NEE4+6f10jvljnqXnw5um/iIiIiIi6i/HfFAsPR+eN9oRERER0UBxzKdPSM+H/6hNnyAiIiKiwavPg+LOzk5s3bo1+f8dO3Zg7dq1KC0txfDhw3u9Hm2aBKdPEBEREVF/6vOg+J133sHUqVOT/7/xxhsBALNmzcLixYt7vR5eKSYiIiKigaLPg+Lzzz8f1lGY9MsrxUREREQ0UGQvpxgfrSvF4ZgBu5UaFVJVJI/gjxuiH0l9ixw94o/I90COKZdzWRo69DiTSiU2K67ElRV4lUgX5bnc/9qh7g4i0S6xbOsueZtup1z2Qb2cJfTmuk51f84/XY6nqS6Tvy4el/yZbG/Qz4PmQFguC8oZPeMq5W22heRtluTpv3UGlTi3oR55m1va5H0dUahkCUE/36UYPTN+zO8TJuqz2ppSuPMKU95vV87xnXsb1HVGIyGxrMAhx2LFiuUyr5xwBiBDSlxcjtvSYiTddrkNLcoQETdySLFY1rSvQyzzt7SJZbGQHI8WV2LVACASkvuuypEVYplDyV1LaJl1AGxKHXny5MJoQm7zDSUGLipXKwAgFJLPS0M5ERIx+fxxGFrgHwCb/D0yE+nXa/XhSmv2BsVKygTTJ4iIiIioP2VtUGwCkC4wcvYEEREREfUnXikmIiIiopyX1UGxdKWYg2IiIiIi6k/Zmz5hAYY0fYKDYiIiIiLqR1m8UmzASpPm0F1GRERERNRfeKW4l2IJIJEmYmTFHjnCa4RXv2VwVKkcU2Wzycv68uVfGkqUKC4AyHfK29zWKmezVBTJyw3xysuV5MnRKwDQ1inv74VnFIhl+1vk7Jr3tgXFsspSJX8GwKtr5Mi2poC83AWnyfv63m45ugcAzhwpn0N72+X629MuR9OU58v16pY3BwBw2eVlXUr1FStxf54M2yxXviu1penL4jETdfpqifrd+DGl8BYUpby/zy9/l0NdLeo6m9vkNsSA3DZHInKsWFdQ7yviysWpPIc8dNCCErUtSuOBblVlXrFsTG2xWBYNKXFt7XKkmKlElAJAXqHc5h8/pkYsGzlsiFjmyjAiM5TKjSv1l7Dk9jVuySvtkrtSAEAwKNdfMCj30XZDjgrMGLWgHIvtKIwpOSgmIiIiopzH9AkiIiIiynlZGxQnlEfaJTgoJiIiIqJ+xCvFRERERJTzsnylOP1Edl4pJiIiIqL+xCvFRERERJTz+ES7Xvr8f+TB7XanvB/ojIrL7G7U48iGVqSur1tbh7zs7iY5gqe8UN0kyn1yptbJasXLsS37W+V9zfRZnjLWI5aVFMgLtwXk/ekIypEtY2rlGB0A2N8qR7JVl8jLba6TI5NK8vSImbc3y9E1gBwLVFMsrzehbDKibQ7ApvYMkTiCccVKDJxDPxGGFMnlQSFVKp7hOIiyYfJEFwqKXCnv72xJfa+bmZBjugBg2y65HWhtbRXL2rvkvqK1I6xusyMmR6DZ5ENBhvRFUZrE0x5KfHLZxHGVYll1Rb5Y1tkhR9ZligbzFshHWlEu72yJ0gU5j+ARDVr9GUqEXlQ+RRCO6u12JCr3/fG4Mi7QcvsOr/s5sF5hwCG9n05Wp09I+8lINiIiIiLqT3yiHRERERHlvKwNiuPW0Xn6yNEWiUR6PAUoEFAeY0ZERDmJfQXR4KPN7DimTEt/ZcuiRYvg8/mSr9ra2uztDBERDUjsK4gGHw6KDzF//nz4/f7kq76+Pns7Q0REAxL7CqLBJ3vTJ2DAJswdNpU77o81t9udNmWCiIioG/sKosGHkWxERERElPOyF8kGOYF1II6Jd+4LwuVKDfSLx+W9fXWbnAkMADUNcsZxNCFfLd/UIecjjsqQiTvJK+/vRuWvf5VFSvZvjXwaVRSru4Mt9XJOpMct19/mupBY5g/Kx7hhq5xDDAAeJX9Ty/7d0yRv84xx+tWkDTvkOtir5DGfNFrOeG7vkEN8F6/SsjkBOOR8UkTkTNTKQjkk+5WdeqjwsAL5WEaVKkGaRAOM1wbkpfna5ivdgd0KquuMx+W+IhaXG6bOiFzWqGb0Ao0BuR2wF8vL5SmjCsOU98dmyzCbUxkYFCvZv8X5ck4xTKUsw+5oucBq4r/2h3BLb+sShrxmUxnOaU9M0OKqWzOcI0Elx1ibBmsewShPzxyWzq/ehx9nbVAcMeWTyjqC8GYiIiIior7K4pViA4ZwrdjK4pxiIiIiIso92bvRzoQ8f4JXiomIiIioH2VtUHxg1o00N2QgziomIiIiosEqe4NijomJiIiIaIDI4qDYlO+05J12RERERNSPsjgotuRA4gEYVFzgtcHlSs1oCSmRN6dUaUEoQExJX2nulG82LFKifbR1AsAuJTqsQEkOsyvxNC+ul9dZXaj/glNdIpc1t8lxMKccL2fwtAXk5er369Fg2qmnRcwUKVF3My/S45b2Pi5/DUs75fpraZePRU03ciqRawCmDZPX+/JOOXbNq8TZwZmnbnNUqbzNkoL034VYjDfk0sDj/fcr5X2l3fY6tC8P4HTKMZyuPLkt7ErIHUJ9c4e6zcoyOa7M65H3x66knNmVRtRh6H2FQ7sBX2m4HUpj6FDayUytixbRqUaOKbFqcWSIZLMOb4e12LU9LXLZ3sYudX/aO+UxTldYriCjQD4OM0PFq4l2R2FMmcVBcRyQPmBLH0wSERERER1NWZ4+IfwmwekTRERERNSPspg+YeJoPH2EiIiIiOhIZW9QbCYAQ5g/Y/KxrkRERETUf3ilmIiIiIhyHucUExEREVHOy+KgOCGnTEj5xVl0zsdsyPOmpmW8+JYc9ZHn1rNFXtkhV/8Yn5zA4bHJ2zyuQq+7pg45CiWoRFydMkZeLp6Q47Q2NOqnWHtI/gVoVFSuA8OQM2baAvJyhXlaVhmQ55HLI1F5Xxvb5Hr/0wt6HNmYYXIdtQTkSJw12+Xz4MKPyRFPM8bKkXUA0K6k8Eyuls8RrxLpp1QrAGBbq1wHUyvTH2csykg2GpjSnZmJqPzz4bCSmQUgqkR/Gg65fensDIhlu/b71W2WFMhtSLlviFjmkxPi4HDIuXQZmgi9XIlzsysxZ07I+6Mkpx2gXLyzlIUNbZuGHs0XVfYpplxL3N8ml+3YLUeG7mnS+4rOiNwGx5XhpV0deurtuhi7ppRpyxyKOcVERERElPOyeKNdHDCE3whM5hQTERERUf/J4o121r9fUhkRERERUf/IbvqEOCeHN9oRERERUf/J3qA4HgFsUk4xp08QERERUf/p90Gxy+VCVVUVGna9qP5cVVUVXC79TkwiIiIioqOh3wfFHo8HO3bsQDSq5NPgwODZ4/H0015l9tLqBFyu1Cvbo4d6xWWa/Z3qOgMd+8WykydWiGW1nXLEjMuhx5kUeOT52hv2ysuu3yFfvR9aKi93EvSr/k45nQZ72uT1tnbK6813y8fYpUQbAcBl58nxRq+u6RDLRlQ5xbKWgBxZBwAb6+TYmxUfNItl+YXVYtn6bXLEUzBDlNk5E+U6+GBXSCzb2SSv84wqvQ72B+TApXxv+pMkas8U4kTU/ywLMNM0QWGly4vHMpzLDvkCUahLbgsjptzFt3XobfPOva1iWUWh3EYUOAvFslIlndKXobvX7jSya4Vac6d0B5kel2BTQuIMU144qhyJadPb5i4lIa2xXS77YJvcd23aJue1NcqJfgCAWCJfLDOUsYgBub+0MswUMLSvihTNp0T2HSor0yc8Hs+AGvASERERUW7jpRYiIiIiynkcFBMRERFRzuOgmIiIiIhyHgfFRERERJTzOCgmIiIiopzHQTERERER5bwsPub5o2XSeDe8XnfK+y+91SUuk5DjhAEAx1WWi2XhqJxl2N4ll42p1j9SZ1TO66sulHe43i8HCpcVyOsMZcjEbYrI5Q4lfNJuk8saO+Tf9ZxqoCXwt9flPMeOsLyvvgK57ur269tMmPJ6x1VVimUeh5znWN8u10Fc2R4A/O0tOePYpeRKB9WsVT0n0uvUPuv0+yu9T5RN7UEgnqYZbg/Ky3TFlC8WANNI7Xu6BSNyBrhll5eLG/r3p6FVDsV9f7ucsW+Ly9/12go51792iP6wrhI5Ehc2Q68/iRZbnyna1lQ2mbDJbaHStSMgx8ADAPY0ywvXN8hjkfVb5RD5bXvkPi9oyp8XAISVIWTCkivIroRHGxnOSy2x2rLS98NWptDpg/BKMRERERHlPF4pPkQkEkEk8uFvyIFAhke6EBFRzmFfQTT48ErxIRYtWgSfz5d81dbWZnuXiIhogGFfQTT4cFB8iPnz58Pv9ydf9fX12d4lIiIaYNhXEA0+nD5xCLfbDbdbvjmBiIiIfQXR4MMrxURERESU83iluJcaWmxwu1N/h2jvlKM+Sgr13zlOHtb3GCoAmD65UCyrb5BjdABg4y45OuyKqfJ63/1AzhOaMLZALHtqRae6P5VFcv0VeuU6KPPJp67bKS/X3C7HmAHAuFFyHazeKEfXvLZZ/qyH+fRItsoiuey8Uzxi2evr5fyeMVXyNodX6dFH7R1yHRXmy/W+fpsc5dYR0mN2Sgvl/T1uePqrcWF5c0RZs3RlMzx5qe3w/i657W3M8P0IQb4iXViaJ5aZpty+2iw5yu3AsnI7sadVbiMCgX1i2QaHXAcFhr4/kz82XiwrK5L3NU+5mJ9QNqm32kBC+cg6lbapuSMqlm3apk/BaQvK/XujXy7zd8n9UwhOsSzu0IeIhl3un/JsSsU7Di9CDwBgk89pw5l+m3FDrvOU1fd5h4iIiIiIBhkOiomIiIgo53FQTEREREQ5j4NiIiIiIsp5HBQTERERUc5j+kQGlnXgHtSDH+d5sHhMvqsxFtV/54ib8v2tNku+tTUclj+2aFRPn0jE5buGw2F5WW294bB8J2kirt/1GY/JdRBzyHUQVerWUOouFpPvfgbkzxnQP2szLu9PPK7fxxxT7oCORLT1Kueeuk51dxCNyudIRLl7XKufRFy/u147D6Tzsvv97u8oUTZ1n4fhYPqUmkhQ/u5EwxnaiKjSplly+2sp6RNGhvQJw5QbipihHItDPhat/XBkSJ/o6pQfo+025PQJU+mCjlX6RJfSxgY75R2Szp1ukZC8bFTpv2NhpR8xlfQJm/6ZJOxyuWHTKl7r1/RtqukTQspEPHwgBas3fYVhsUdR7d69m4/vJBrA6uvrMWzYsGzvBuU49hVEA1tv+goOijMwTRN79+5FYWEhDMNAIBBAbW0t6uvrUVSkBMzmKNaPjvWj60v9WJaFjo4O1NTUwGbjTDDKLvYVfcP60bF+MuttHfWlr+D0iQxsNlva3yyKiop4oipYPzrWj6639ePz+fphb4gyY19xeFg/OtZPZr2po972Fby8QkREREQ5j4NiIiIiIsp5HBT3kdvtxoIFC+B2K8/1zmGsHx3rR8f6ocGC57KO9aNj/WR2LOqIN9oRERERUc7jlWIiIiIiynkcFBMRERFRzuOgmIiIiIhyHgfFRERERJTzOCgmIiIiopzHQTERERER5TwOiomIiIgo53FQTEREREQ57/8DI6GiiemdQnsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, axarr = plt.subplots(1,2,figsize=(10, 3))\n",
    "axarr[0].axes.xaxis.set_ticklabels([])\n",
    "axarr[0].axes.yaxis.set_ticklabels([])\n",
    "axarr[1].axes.xaxis.set_ticklabels([])\n",
    "axarr[1].axes.yaxis.set_ticklabels([])\n",
    "\n",
    "pcm = axarr[0].imshow(np.mean(np.diag(M).reshape(3,32,32),axis=0),cmap='cividis')\n",
    "axarr[0].set_title(\"M matrix diagonal\")\n",
    "f.colorbar(mappable=pcm, ax=axarr[0], shrink=0.8,location=\"left\")\n",
    "axarr[1].imshow(torch.moveaxis(trainset0[3][0],0,2))\n",
    "axarr[1].set_title(\"Sample Image\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b8bb8a2",
   "metadata": {},
   "source": [
    "Here we plot the diagonal (averaged across channels) of the M matrix to see which coordinates are being highlighted by RFM. The center of the image is highlighted (where the digits appear). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "451d553b",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}