rep_align_demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torchvision\n",
    "import torchvision.transforms as T\n",
    "import matplotlib.pyplot as plt\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing the Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform=T.Compose([T.ToTensor(), T.Normalize((0.5, 0.5, 0.5), (0.5, 0.5, 0.5))])\n",
    "train_dataset = torchvision.datasets.CIFAR10(root='./data', train=True, transform=transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binary_clf_dataset(data, labels, neg_class=0, pos_class=1, size=1000, seed=0):\n",
    "    rng = torch.Generator()\n",
    "    rng.manual_seed(seed)\n",
    "    labels_tensor = torch.Tensor(labels).long()\n",
    "    mask = torch.zeros(labels_tensor.shape, dtype=torch.bool)\n",
    "    idx = torch.nonzero(labels_tensor == neg_class)[:, 0]\n",
    "    idx = idx[torch.randperm(len(idx), generator=rng)[:size//2]]\n",
    "    mask[idx] = True\n",
    "    idx = torch.nonzero(labels_tensor == pos_class)[:, 0]\n",
    "    idx = idx[torch.randperm(len(idx), generator=rng)[:size//2]]\n",
    "    mask[idx] = True\n",
    "    binary_labels = labels_tensor[mask]\n",
    "    binary_labels[binary_labels == neg_class] = 0\n",
    "    binary_labels[binary_labels == pos_class] = 1\n",
    "    binary_labels = binary_labels.tolist()\n",
    "    binary_data = data[mask]\n",
    "    return binary_data, binary_labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_dataset.data, train_dataset.targets = binary_clf_dataset(train_dataset.data, train_dataset.targets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dataset size: 1000, number of classes: 2\n"
     ]
    }
   ],
   "source": [
    "print(f'dataset size: {len(train_dataset)}, number of classes: {max(train_dataset.targets)+1}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_loader = torch.utils.data.DataLoader(train_dataset, batch_size=64, shuffle=True, num_workers=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training the Network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 0\n",
    "torch.manual_seed(seed)\n",
    "torch.cuda.manual_seed(seed)\n",
    "torch.backends.cudnn.deterministic = True\n",
    "torch.backends.cudnn.benchmark = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sequential(\n",
      "  (0): Flatten(start_dim=1, end_dim=-1)\n",
      "  (1): Linear(in_features=3072, out_features=128, bias=False)\n",
      "  (2): ReLU()\n",
      "  (3): Linear(in_features=128, out_features=128, bias=False)\n",
      "  (4): ReLU()\n",
      "  (5): Linear(in_features=128, out_features=128, bias=False)\n",
      "  (6): ReLU()\n",
      "  (7): Linear(in_features=128, out_features=1, bias=True)\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "# Play around with the network's architecture. Change the depth, width, activation, etc..\n",
    "# For example try linear activation by replacing nn.ReLU() with nn.Identity().\n",
    "model = nn.Sequential(\n",
    "    nn.Flatten(),\n",
    "    nn.Linear(32*32*3, 128, bias=False),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(128, 128, bias=False),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(128, 128, bias=False),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(128, 1))\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using crossentropy loss with sigmoid activation and {0,1} labels.\n",
    "# It is equivalent to using the logistic loss in the paper with {-1,+1} labels.\n",
    "# See here: https://en.wikipedia.org/wiki/Loss_functions_for_classification#Logistic_loss\n",
    "# We will convert the labels to {-1,+1} later for computing alignment.\n",
    "critertion = nn.BCEWithLogitsLoss()\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch: 9, train loss: 0.055221106857061386\n",
      "epoch: 19, train loss: 0.017858062419691123\n",
      "epoch: 29, train loss: 0.00027390173181629507\n",
      "epoch: 39, train loss: 0.00011177455075994658\n",
      "epoch: 49, train loss: 6.065116349418531e-05\n",
      "epoch: 59, train loss: 3.754453257442947e-05\n",
      "epoch: 69, train loss: 2.5314709887425124e-05\n",
      "epoch: 79, train loss: 1.8092143818648765e-05\n",
      "epoch: 89, train loss: 1.443024348191102e-05\n",
      "epoch: 99, train loss: 1.0185814176111307e-05\n"
     ]
    }
   ],
   "source": [
    "num_epoch = 100\n",
    "model.train()\n",
    "for epoch in range(num_epoch):\n",
    "    train_loss = 0.0\n",
    "    batch_cnt = 0\n",
    "    for batch_X, batch_y in train_loader:\n",
    "        optimizer.zero_grad()\n",
    "        predict = model(batch_X)\n",
    "        loss = critertion(predict, batch_y.float().unsqueeze(1))\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        train_loss += loss.item()\n",
    "        batch_cnt += 1\n",
    "\n",
    "    if (epoch + 1) % 10 == 0:\n",
    "        print(f'epoch: {epoch}, train loss: {train_loss/batch_cnt}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting Representations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['input', '0', '1', '2', '3', '4', '5', '6', '7']\n"
     ]
    }
   ],
   "source": [
    "from torchvision.models.feature_extraction import get_graph_node_names\n",
    "all_nodes, _ = get_graph_node_names(model)\n",
    "# These are all the model's nodes.\n",
    "print(all_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torchvision.models.feature_extraction import create_feature_extractor\n",
    "# These are the nodes we extract the representations from.\n",
    "# They correspond to the flattened input and each layer's post activation hidden representation.\n",
    "# If you change the architecture you may need to change this array.\n",
    "nodes = ['0', '2', '4', '6']\n",
    "model.eval()\n",
    "feature_extractor = create_feature_extractor(model, return_nodes=nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([1000, 3, 32, 32]) torch.Size([1000])\n"
     ]
    }
   ],
   "source": [
    "input_data = []\n",
    "targets = []\n",
    "for batch_X, batch_y in train_loader:\n",
    "    input_data.append(batch_X)\n",
    "    targets.append(batch_y)\n",
    "input_data = torch.vstack(input_data)\n",
    "targets = torch.hstack(targets)\n",
    "print(input_data.shape, targets.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('0', torch.Size([1000, 3072])), ('2', torch.Size([1000, 128])), ('4', torch.Size([1000, 128])), ('6', torch.Size([1000, 128]))]\n"
     ]
    }
   ],
   "source": [
    "reps = feature_extractor(torch.Tensor(input_data.data))\n",
    "print([(key, rep.shape) for (key, rep) in reps.items()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing Alignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function computes the alignment curve efficiently.\n",
    "def get_alignment_curve(dot_products, y, s, grid):\n",
    "    cumu = np.zeros(len(grid))\n",
    "    cumu_ind = 0\n",
    "    for s_ind, this_s in enumerate(s):\n",
    "        while grid[cumu_ind] > this_s:\n",
    "            cumu_ind += 1\n",
    "            cumu[cumu_ind] += cumu[cumu_ind-1]\n",
    "        cumu[cumu_ind] += dot_products[s_ind]\n",
    "    cumu[cumu_ind+1:-1] = cumu[cumu_ind]\n",
    "    cumu[-1] = (y**2).sum()\n",
    "    return cumu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4bElEQVR4nO3dd3xc1ZXA8d+Zqi7ZsmRsy73b2LgXeo8xxQQIOIHQ4yxtIWQ3kOxmSbKbBNLAEELWCQYTA6YGG0KyOFRjjHtv2LhKuHfLllXm7h/vjawy0kjWjN7Mm/P9fOYzmjdv3jvjsXTm3HvfvWKMQSmllGqMx+kAlFJKJT5NFkoppaLSZKGUUioqTRZKKaWi0mShlFIqKp/TAcRDu3btTLdu3ZwOQymlksrixYv3GmMKIj3nymTRrVs3Fi1a5HQYSimVVERka0PPaTOUUkqpqDRZKKWUikqThVJKqahc2WehlFJOqaiooLi4mLKyMqdDaVBaWhpFRUX4/f4mv0aThVJKxVBxcTHZ2dl069YNEXE6nHqMMezbt4/i4mK6d+/e5NfFrRlKRKaKyG4RWVVjW1sRmS0iG+z7NvZ2EZEnRWSjiKwQkWE1XnOLvf8GEbklXvEqpVQslJWVkZ+fn5CJAkBEyM/Pb3blE88+i+eBcXW2PQy8b4zpDbxvPwa4DOht3yYBz4CVXIBHgNHAKOCRcIJRSqlElaiJIuxU4otbsjDGfALsr7N5AjDN/nkacHWN7S8Yy+dAnoh0AL4GzDbG7DfGHABmUz8Bxcz6pR/x8s1jWPv8M1Ts3Bmv0yilVNJp7T6L9saYHfbPO4H29s+dgO019iu2tzW0vR4RmYRVldClS5dTCu6rVQs5fdEhWPAkGx99El9hIelnDCZt4Ol427TBk5GOJz0dSbfuETmZoUWsm8eLeD3g9SJeL3g8tbN4zf189vM+H+LxnHyu5jFp2jcA8dnn8/kS/luNUiq+/vGPf3D//fdTVVXFnXfeycMPPxz9RVE41sFtjDEiErOVl4wxU4ApACNGjDil4/b82p28tPdpFhzP4mdtvk3Opn0cX7GcI7P/GaswW4fHYyUOTRrKIVkXXEDR5CecDiMlVVVVcc899zB79myKiooYOXIkV111FQMGDGjRcVs7WewSkQ7GmB12M9Nue3sJ0LnGfkX2thLg/DrbP4pXcDlZmVx2MIO/d/Pxs3aLeem7L+Hz+AiVllJ1tBRz/Bih48cJHTtG6HgZmJD1QmPAGIwxEAphqqqgKoSpqoRQjX1sJhSCkLGer6qy9686eQwAU/s1jTIGE7KOYSqrMJUV1vGUcsDRjz6mbOVKp8NIWQsWLKBXr1706NEDgIkTJzJz5sykSxazgFuAR+37mTW23ysiM7A6sw/ZCeX/gF/U6NS+FPhhvILLCvrYHirigUMl/JdnLTPWzeCmATfhyczEk5kZr9Mq5SpVBw9x9KOPnA4jIfz07dWs+epwTI85oGMOj1w5sMHnS0pK6Nz55HfvoqIi5s+f3+LzxnPo7MvAPKCviBSLyB1YSeISEdkAXGw/BngX2ARsBP4E3A1gjNkP/Dew0L79zN4WFz6vh62eLkw4WMJZHcby1NKn2FW6K16nU8qVxO/HlJc7HYaKsbhVFsaYbzbw1EUR9jXAPQ0cZyowNYahNarY3w1PpeE/elzD13cv4VcLf8Vvz/9ta51eqaQngQCmosLpMBJCYxVAvHTq1Int20+OCyouLqZTp4jjgppF54aqY3ea1c7X+eh+Jg2exHtb32NO8RyHo1IqeYjfT0iThWNGjhzJhg0b2Lx5M+Xl5cyYMYOrrrqqxcfVZFHH0czOVOCH3Wu4deCtdM/tzs/n/5zV+1ZTFdJOY6WikUAAKipODtZQrcrn8/H73/+er33ta/Tv35/rr7+egQNbXuHo3FB1ZKWnsd1bRI/d6wh4A/x4zI+ZNHsSE9+ZSLY/m+HthzPytJGMPG0kfdv2xSOab5WqSezJ6UxFhZU4VKsbP34848ePj+kxNVnUkZPuZyNd6LF7LQAjTxvJe9e+x4KdC1i4cyGLdi3io+KPAMhPy+fsTmdzTtE5jO04lpxAjoORK5UYwgnClJeDJgvX0GRRR3aaj3WhTlx6+GMoOwRpuRRkFHB5j8u5vMflAOwq3cWCnQuYUzKHD7d/yMwvZ+IVL2cUnMFZnc5ibIexDMgfgNfjdfjdKNX6qpOF9lu4iiaLOnLS/Kwo7wB+YM966Dyq3j7tM9tzZc8rubLnlVSGKlm5dyVziufwacmnPLX0KZ5a+hQ5gRxGdxjNmR3P5LLul5Hp1+s0VGqobobS4bOuosmijpx0H+tC9gUtu9dETBY1+Tw+hhYOZWjhUP512L+y7/g+5u+Yz7wd85j31Txmb53N7xb/jol9J3Jj/xvJT89vhXehlHO0snAnTRZ15KT5KTHtCPkz8Nj9Fs2Rn57P+B7jGd9jPMYYVu5dyXOrnuPPK//MC2te4OpeV3PrwFspyi6KQ/RKOU8rC3fSZFFHTrofg4cTeb1JP4VkUZOIMLhgMI9f8DibD23m+dXP88aGN3hl/St0zOxI7za9rVteb/q17Uf33O46Y6xKehI4ORpKuYcmizqy06x/kqO5fUjf8XHMjts9tzs/PfOn3H3G3by96W2+2P8FGw5uYG7JXCpNJQBFWUVc2u1SLu12KQPaDtDEoZKSVhbO2759OzfffDO7du1CRJg0aRL3339/i46pyaKOnDTrP/qBrJ4UlL4GpXshs13Mjt8+sz13Drqz+nF5VTmbD21m5d6V/HPrP5m2ehpTV02lKKuI8zufT3Ygu9brA94AYzqMYWD+QE0mKiHVGjqrHOHz+fjtb3/LsGHDOHLkCMOHD+eSSy5p0cyzmizqyEm3ksWe9O70Adi9FrqfE7fzBbwB+rbtS9+2fbmuz3UcLDvIB9s/4L0t7zFj3YzqqqOmyUzmtMzTuLjLxVzU5SKGFg7VYboqYXi0g9txHTp0oEOHDgBkZ2fTv39/SkpKNFnEUo7dDPVVoLu1Yc+6uCaLuvLS8rim9zVc0/uaiM8fOnGIj7Z/xD+3/pNX17/K9LXTyfZnk+5Lb9LxPR4PXvHi8/jwihePeOpVKBm+DL7e6+tc2fNKAl69qEo1jzZD1fD3h2FnjNf2OG0QXPZo9P1sW7ZsYenSpYwePbpFp9VkUUe23Qy1K9QG0nKt4bMJJDeYy4ReE5jQawKlFaXMKZnDop2LqAzVr0DqMhiqQlVUGftm/1zX9iPb+cm8n/CHZX/g5oE3c12f6/Q6EdVkOnQ2cRw9epRrr72WJ554gpycls0wocmijoDPQ5rfw+ETVVDQ32qGSlCZ/kzGdRvHuG7jYnpcYwzzdsxj6sqp/GbRb5iyYgoTek2gTbBNrf3y0vL4eq+v4/PofyN1klYWNTSjAoi1iooKrr32Wm688UauuSZyS0Vz6G95BDlpfg4fr4DC/rD6r9bypinUmSwinNnxTM7seCYr96xk6qqpTF8zHUP9WUQX71rMz8/6ufaZqGpaWTjPGMMdd9xB//79efDBB2NyTE0WEeSk+zlcVgFFA2Dxc3BkJ+R0cDosRwwqGMTjFzxOZaiy3pTTz69+nieXPknQG+SRsY/oDLwK0MoiEcydO5e//OUvDBo0iCFDhgDwi1/8okUz0WqyiCAnzceRskqrsgCr3yJFk0VYpKam7wz+DmVVZUxZMYWAJ8CPRv9Ih/Oq6soipMnCMWeffXbM1xPRZBFBdpqfg8fKTyaLPeugV73VYBVw75B7Ka8q5/nVzxP0Bvn+iO9rwkhx1WtYaDOUq2iyiCAn3c+2/cesi/EyCxJuRFQiEREeHP4gJ6pOMG3NNKpMFSNPG+l0WK2mY1ZH+rXt53QYCSXcDKWVhbtosoggJ81ndXCDVV0k8IioRCAiPDzqYcqrypm+djrT1053OqRWk+ZNY/6N87W/pgbt4HYnTRYRhDu4jTFIQX9YOh1CIfDoH4SGeMTDI2Mf4dsDvk15VWp8o3xr41u8tO4lKkIVBL1Bp8NJGOL1gsejHdwuo8kigpw0PxVVhhOVIdIK+0NFKRzaDm26Oh1aQhMReub1dDqMVrNw50LAmt9Lk0VtEghoZeEy+lU5gvDMs4ePV0DHodbGf/wQThx1MCqVaMJToZyoOuFwJIlH/H5MuSYLN9FkEUF4MsHDZRXQcQiMewy++Ds8Nw4OFTsbnEoY4Wqiokr/KNYlgYA2QzmsqqqKoUOHcsUVV8TkeJosIghPJnjouD3f0ph/gW+9Bge2wpQLoHiRg9GpROH3Wl8qtLKoT5uhnDd58mT69+8fs+NpsoigVmUR1vtiuGM2BDLgufGw8nWHolOJIlxZlIf0G3RdVjOU/rs4pbi4mL/97W/ceeed0XduIu3gjiC8ANKRsjozuRb2gzs/gFdugjfugC8/gEv/BzLaOhClclrAY/VZpMror+aQgF8rC+CxBY+xbv+6mB6zX9t+PDTqoUb3eeCBB/jVr37FkSNHYnZerSwiyKnZwV1XZj7cPBPO/h6seAV+PwKWz7AmG1QpJdzBrcmiPvFrn4VT3nnnHQoLCxk+fHhMj6uVRQQRm6Fq8gXg4p/AoG/A2/fDX78Ly1+Gy38H+akzdDTVVScLbYaqRysLS7QKIB7mzp3LrFmzePfddykrK+Pw4cPcdNNNTJ/esotltbKIIOjzEPB6OHw8yoJC7QfC7e/B+N9AyRL4w1j4+0NwqKR1AlWOqu6z0MqiHu2zcM4vf/lLiouL2bJlCzNmzODCCy9scaIAh5KFiHxPRFaLyCoReVlE0kSku4jMF5GNIvKKiATsfYP24432891aIT5y0n0caaiyqMnjgVHfgXsWWJXGwj/Dk0Pg7QfgwJY4R6qc5PdYFagmi/o8OnTWdVo9WYhIJ+BfgRHGmNMBLzAReAx43BjTCzgA3GG/5A7ggL39cXu/uMtO83O4bgd3Y3I6wNVPw31LYOi3YdmL8OQweHMSLHsJdqyASv3lcZNwZaFDZ+sTvw6dTQTnn38+77zzTkyO5VSfhQ9IF5EKIAPYAVwIfMt+fhrwE+AZYIL9M8DrwO9FREysJ2uvo9Zkgs3Rpitc8Ts499/gs6dg8TSrIxzA44OCflbzVbvekN/LurXtAQFd4zrZaAd3wySgzVBu0+rJwhhTIiK/AbYBx4H3gMXAQWNM+Kt8MdDJ/rkTsN1+baWIHALygb01jysik4BJAF26dGlxnNWr5Z3yATrCuF9aQ2v3fQk7V8CuVbBzFWyeczKBhOX3guv/Au0HtCxw1Wo0WTRMKwv3afVkISJtsKqF7sBB4DVgXEuPa4yZAkwBGDFiRIurjpw0P18dPN7Sw4DHCwV9rNug605uP3EU9m+C/V/Cvo0w/3/h9dvhOx9YF/6phKejoRqmHdzu40Qz1MXAZmPMHgAReRM4C8gTEZ9dXRQB4SFFJUBnoFhEfEAusC/eQVod3M3os2iuYBZ0GGzdADoOg+nXwHv/AVc8Hr/zqpgJX5SnfRb16XQf7uPEaKhtwBgRyRBr/c2LgDXAh0D4q/ctwEz751n2Y+znP4h3fwWEO7hb8T97r4vgzH+FRVNhzazWO686ZeHKQicSrE8rC/dp9WRhjJmP1VG9BFhpxzAFeAh4UEQ2YvVJPGu/5Fkg397+IPBwa8SZk+ajrCLEicqq1jid5cIfW1Oiz7oXDm5vvfOqU+IRDz6PTyuLCCQQIKSVhas4cp2FMeYRY0w/Y8zpxphvG2NOGGM2GWNGGWN6GWO+YYw5Ye9bZj/uZT+/qTViDF/FHdemqLp8Abj2WQhVWUNuq1rx3OqUBL1B7bOIQJuhnHfw4EGuu+46+vXrR//+/Zk3b16LjqdXcDcgPJngKQ2fbYn8nta0Ids+gzm/ad1zq2YLeAI6GioC8fuhogITCjkdSsq6//77GTduHOvWrWP58uUtnq5ck0UDctKtvv9WrSzCzrgBBk+Ejx+zLuZTCSvg1WQRiQSs/hxTqdWxEw4dOsQnn3zCHXdY1zYHAgHy8vJadEydSLAB2WlRJhOMt8seg/V/h09+BTe0fF4XFR8Bb0CboSIQv/X7Y8rLwU4cqWjnL37BibWxnaI82L8fp/3oR43us3nzZgoKCrjttttYvnw5w4cPZ/LkyWRmnvrFv1pZNOBkM5RD34zS82D0d2Ht27BrjTMxqKiC3qBWFhFUVxbab+GIyspKlixZwl133cXSpUvJzMzk0UcfbdExtbJoQLgZyrHKAmDMXfD5H6y+i+umOheHapDf49dkEUGtyiKFRasA4qWoqIiioiJGjx4NwHXXXdfiZKGVRQMc6+CuKaMtjLwDVr0Jezc4F4dqUNAb1KGzEVRXFimeLJxy2mmn0blzZ9avXw/A+++/z4ABLZtKSJNFAzICXrwecaaDu6ax94EvDeb8ztk4VETawR2ZBOzKQpuhHPPUU09x4403MnjwYJYtW8aPWljlaDNUA0SE7DSfs81QAFkFMOI2a+6o834Abbs7G4+qxe/1c6zimNNhJBxthnLekCFDWLRoUcyOp5VFI3LS/M42Q4WdeZ81IeGnOmdUogl69KK8SLSD2300WTQiJ93XvAWQ4hZIR3tBpZd0GpAEo81QkWll4T6aLBqRk+Zv2tKqreHsBwADcyc7HYmqQZNFZJ4UryxaYa7TFjmV+LTPohHZaT627E2Q9ui8LnDGN2HJNNixrPZzHj8U9LWmOz9tMBQO0DUxWknAG9DRUBGkcmWRlpbGvn37yM/Px5pYO7EYY9i3bx9paWnNep0mi0bktPY05dGc9xAc2weVZbW3VxyH1W/C4uesx+KBvK7gC9bez+uH9oOg0zDoNBzan25NXqhOmU4kGFkqD50tKiqiuLiYPXv2OB1Kg9LS0igqKmrWazRZNCInPUE6uMPyOsM3X478nDFwcBvsXGkt4bp3A5g606uXl8LG2bD8JeuxN2BVJN46SSW9DVw/TdcFbwKdSDCyVO7g9vv9dO/uvlGLmiwakZPmp7S8isqqED5vgnfviECbrtat/xUN72cMHNoOJUugZDHsWWdNiR52/ICVUHaugi6j4x93kgv3WRhjErLJwSmp3AzlVposGhGe8uPoiUryMlzSXCNi9X/kdYGBV9d/fs96eHqUVaVosogq4A1gMFSaSvzidzqchJHKlYVbJfjXZWdlOz2ZoBNyO1v3B7c6G0eSCNpNeNoUVVu4sghpZeEamiwakZOWAJMJtrZABmQWWJWFisrvsf4oarKoLVxZoJWFa2iyaER4adWE6uRuDXldNFk0Ubiy0OGztWll4T6aLBqR4/QCSE7RZNFkAa/1DbqiKsX+j0SRykNn3UqTRSNOrmmRQn0WYCWLQ9tB10+Oyu+1vlBoZVGbeL3g9WoHt4tosmhEdiKsaeGEvC5QVQ5HdzkdScILeuwObr0wrx7x+zHlKfa742KaLBqRHfQhkoqVRVfr/pBOWhhNuBlKO7jrk0BAKwsX0WTRCI9HyAr6Uq+yqB4+q/0W0WiyaJhVWei/i1tosogi4eaHag15eq1FU4WThfZZ1KeVhbtosogiO83n/NKqrS2QCRnttLJoguqL8rTPoh6tLNxFk0UUCTeZYGvR4bNNEvDo0NmGSECThZtosojCaoZKscoCNFk0kTZDNUybodxFk0UUOekp2MENdrLQay2iqe7g1maoerQZyl00WUSRkh3cYF9rcQJKdzsdSULTiQQb5vFrZeEmUZOFiNRbxSPStuYQkTwReV1E1onIWhEZKyJtRWS2iGyw79vY+4qIPCkiG0VkhYgMa8m5mysnzcfRE5WEQom9pm7Mha+10KaoRulEgg3TPgt3aUpl8UaEba+38LyTgX8YY/oBZwBrgYeB940xvYH37ccAlwG97dsk4JkWnrtZctL9GANHy1Os3yJPr7VoCu2zaJhoZeEqDS5+JCL9gIFArohcU+OpHKB5K33XPm4ucC5wK4AxphwoF5EJwPn2btOAj4CHgAnAC8YYA3xuVyUdjDE7TjWG5sipMeVH+OeUoBfmNYnP48MrXq0sItDKwl0aWymvL3AFkAdcWWP7EeA7LThnd2AP8JyInAEsBu4H2tdIADuB9vbPnYCa804U29tqJQsRmYRVedClS5cWhFdb9WSCxyuhTcwOm/iCWZCRr8miCcJLq6raxB/QZOEiDSYLY8xMYKaIjDXGzIvxOYcB9xlj5ovIZE42OYXPbUSkWZ0ExpgpwBSAESNGxKyDIVxNHEnVTm5NFlEFvAEdDRWBDp11l6aswb1RRH4EdKu5vzHm9lM8ZzFQbIyZbz9+HStZ7Ao3L4lIByA8DKcE6Fzj9UX2tlbRMS8dgC92HWF0j/zWOm1iyOsCu9Y4HUXCC3i0sohE/H5CFfrv4hZN6eCeCeQC/wT+VuN2SowxO4HtItLX3nQRsAaYBdxib7vFPi/29pvtUVFjgEOt1V8B0DU/g6I26Xz8xd7WOmXiCK9rYVJsJFgzaTNUZBIIgE5R7hpNqSwyjDEPxfi89wEvikgA2ATchpW4XhWRO4CtwPX2vu8C44GNwDF731YjIpzbp4CZS0sorwwR8KXQpSl5XaGyDI7uhuz20fdPUQFvQEdDRWBVFpos3KIpyeIdERlvjHk3Vic1xiwDRkR46qII+xrgnlid+1Sc16eAl+ZvY8m2A4xJpaaoPHugwMFtmiwaEfQGtc8iAu2zcJemfE2+HythlInIYRE5IiKH4x1YIjmzZz4+j/DxF3ucDqV15epU5U0R8AR0IsEIxO+HigqMThnjClGThTEm2xjjMcakGWNy7Mc5rRFcoshO8zOsSxs+SbVkoRfmNYk2Q0UmAeuCRa0u3KEp032IiNwkIj+2H3cWkVHxDy2xnNe3gNVfHWbPkRT6oxDMhvS2urxqFDp0NjJNFu7SlGaoPwBjgW/Zj48CT8ctogR1bu8CAOZsSLXqQq+1iEZHQ0UmfusaJb0wzx2akixGG2PuAcoAjDEHgEBco0pAAzvmkJ8ZSMGmKE0W0eh1FpFJwE4WWlm4QlOSRYWIeAEDICIFQMr1WHk8wjm92/HJhr2pNQNtOFnotRYN0j6LyMRvN0NpZeEKTUkWTwJ/BQpF5OfAp8Av4hpVgjq3TwH7S8tZ/VUKDQYLX2tRmmIVVTMEvUEdDRWBVhbuEvU6C2PMiyKyGOsaCAGuNsasjXtkCegcu9/ikw17GFSU63A0raTmtRZZhc7GkqD8Hj8nQlpZ1KV9Fu7S1MuRdwFzgM+A9NZegChRFGQHGdgxh4/Xp9C37Dy91iKaoDeofRYRVI+G0mThClErCxH5b6y1J77E7rew7y+MX1iJ67w+BUz5ZBNHyirIToX1LXRdi6h0NFRkHh066ypNqSyuB3oaY843xlxg31IyUYDVb1EZMnz25T6nQ2kdaTmQ3kaTRSMC3gBVpoqqUJXToSQUbYZyl6Yki1VYCyApYFiXNmQGvKk19UdeFzioF+Y1JLy0ql6YV5telOcuTZlI8JfAUhFZBVT34hljropbVAks4PNwZq92fPLFHowxiIjTIcVfXhfY84XTUSSsoDcIQHlVOem+dIejSRxaWbhLU5LFNOAxYCUpeH1FJOf2KWD2ml1s2H2UPu2znQ4n/vK6woZ/WtdapEJybCa/x/qjqP0WtWll4S5NaYY6Zox50hjzoTHm4/At7pElsK8NaE9W0Md//nUVValwgV5eF6g8bq1roeoJVxZ6YV5tWlm4S1OSxRwR+aWIjBWRYeFb3CNLYIU5afz0qoEs2LKfP378pdPhxF+RvfTI2lnOxpGgtM8isnBlEdJk4QpNaYYaat+PqbEtZYfOhl0zrBMfrN/N47O/4Jze7RhclOd0SPHTcRh0HAoL/gQj79SmqDoCHjtZaDNULdoM5S5NWc/iggi3lE4UYC23+ourB1GQHeSBGcs4Vl7pdEjxIwKjJsHe9bD5E6ejSTjVlYUmi1q0GcpdmrKeRVBEviUiPxKR/wrfWiO4RJeb4ee315/B5n2l/Pc7Lp8BZeA1kJEPC6Y4HUnCCScL7bOoTSsLd2lKn8VMYAJQCZTWuCngzJ7tmHROD15esI33Vu90Opz48afBsFtg/bt6gV4d4Q5unUywNq0s3KUpfRZFxphxcY8kiT14aR8+3biXh99cyemdcumY59Kx9iNuh7lPwMJn4ZKfOh1NwvB7rT+KWlnUJl4veL1aWbhEUyqLz0RkUNwjSWJBn5fJE4dSURnipmfns++oS/9o5HWGfpfDkheg4rjT0SSMoMe+KE9HQ9Ujfj+mXJOFGzQlWZwNLBaR9SKyQkRWisiKeAeWbHoVZvHsrSMpOXCcW55bwOEyl/6CjJoEx/fDqjecjiRhaAd3wyQQ0GYol2hKsrgM6A1cClwJXGHfqzpGdW/LH28azrodR7hz2iLKKlw4sVy3c6CgP8z/X109z6bJomESCGgzlEs0JVkciXD7Kp5BJbML+hXyuxuGsHDLfu5+cQkVVS6bIUUERn0Hdq6A7QucjiYh6EV5DbOaofTfxQ2akiyWAHuAL4AN9s9bRGSJiAyPZ3DJ6qozOvI/V5/OB+t28/1Xl7tvSpDBN0AwV4fR2mpOJKhqk4BfKwuXaEqymA2MN8a0M8bkYzVLvQPcDfwhnsElsxtHd+UH4/oya/lX3PrcAg6UuugPSTALht4Ia96C0r1OR+M4vYK7YVpZuEdTksUYY8z/hR8YY94DxhpjPgeCcYvMBe4+vxePXjOI+Zv2c8VTn7Kq5JDTIcXO4BsgVAkbZjsdieN8Hh+C6NDZCLTPwj2akix2iMhDItLVvv0A2CUiXnTK8qgmjurCq/8ylpAxXPvMZ7y2yCWLCJ02GLLaw4b3nI7EcSJiLa2qfRb1aGXhHk1JFt8CioC37FsXe5sXa8lVFcWQznm8c9/ZDO/ahn9/fQX/8deVnKhM8pFSHg/0ugS+fB+qXDwvVhMFPLoOdyQevw6ddYumTCS41xhznzFmqH271xizxxhTbozZ2BpBukF+VpAXbh/Fd8/rwYvzt3H105+xcfcRp8Nqmd6XQNkhKF7odCSOC3g1WUSizVDu0WCyEJEn7Pu3RWRW3VtLTywiXhFZKiLv2I+7i8h8EdkoIq+ISMDeHrQfb7Sf79bSczvF5/Xww8v68+wtI9h1uIwrnvqUlxdswyTr9Qo9LwDxalMUVrLQPov6tBnKPRqrLP5i3/8G+G2EW0vdD9ScqvUx4HFjTC/gAHCHvf0O4IC9/XF7v6R2Uf/2/OP+cxjRtS0/fHMld01fwsFjSfgLlZYLXcZqJzfW8FmdSLA+rSzco8FkYYxZbN9/HOnWkpOKSBFwOfBn+7FgLab0ur3LNOBq++cJ9mPs5y+y909qhTlpvHD7KH40vh/vr9vFZZPnsHVfEk7m2/sS2LUSDqf2dZp+r18riwi0snCPxpqhVtpzQdW9rRSR5S087xPADzg5miofOGiMCfeUFgOd7J87AdsB7OcP2fvXjXeSiCwSkUV79uxpYXitw+MRJp3bkzfvOotj5VX8y/QlyTdFSO9LrfsUry6CnqCOhopAKwv3aKwZKjwHVM3bVVgX45Wc6glF5Apgd7hyiRVjzBRjzAhjzIiCgoJYHjruBhXl8sQNQ1i74zD/NXOV0+E0T2F/yClK+X4L7eCOTPx+QhX67+IGjTVDbQ3fgLbAvcBHwM+Ad1twzrOAq0RkCzADq/lpMpAnIuH1NYo4mZBKgM4A9vO5wL4WnD8hXdCvkPsu7MWri4p5dWESXYshYjVFbfoIKlP3j4Imi8isWWe1snCDxpqh+ojIIyKyDngK2AaIvQb370/1hMaYHxpjiowx3YCJwAfGmBuBD4Hr7N1uwVqhD2CW/Rj7+Q9M0g4fatwDF/fhrF75/HjmKlZ/lURXe/e+FMqPwrZ5TkfiGL0oLzJthnKPxpqh1mF967/CGHO2MeYpIJ4N6g8BD4rIRqw+iWft7c8C+fb2B4GH4xiDo7weYfLEobTJCHD3i0s4dDxJfsm6nwveQEo3RQW9Qa0sItAObvdoLFlcA+wAPhSRP4nIRUBMRyEZYz4yxlxh/7zJGDPKGNPLGPMNY8wJe3uZ/biX/fymWMaQaNplBXn6xqGUHDjOv722PDmuwQhmQdezUrqT2+/xa7KIQAIBqKzEhHRmoGTXWJ/FW8aYiUA/rCaiB4BCEXlGRC5tpfhS0vCubfnh+P7MXrOLN5ac8liC1tX7Uti7Hg5scToSRwS9QR06G4H4rfXJtSkq+TVluo9SY8xLxpgrsTqel2I1Gak4uu3MbgzpnMdj/1jHkWRYojXFh9AGvAEqQknwObUyCVjTt2uySH5NmUiwmjHmgD1E9aJ4BaQsHo/w06sGsufICX7/QRJMwZXfE9p0T9lk4ffoRXmRVFcW2m+R9HzRd1FOOaNzHtePKGLq3M1cP7IzPQuynA6pYSJWdbH4eXjlptgdN6s9jHsMvIn9X1U7uCOTgCYLt0js30DFv3+tH39fuZOfvb2G528bSULPdDL0Rtj+Oez7MjbHqyqHtW9D9/NgwFWxOWachJuhQiaER5pVsLuaNkO5hyaLBFeQHeT+i3vzP39bywfrdnNR//ZOh9SwDmfAdz+J3fFCVTD5DFj456RIFgAVoYrqNbmVNkO5iX4FSgK3nNmNXoVZ/OydNcm/aFJzeLww/FbY/DHs3eB0NI0Kr8Ot/Ra1aWXhHposkoDf6+GRKwewdd8xnv10s9PhtK5hN4PHD4umOh1Jo8LVhPZb1KaVhXtoM1SSOKd3AZcOaM9T729k855S+p6WTZ/22fQ9LZvC7GBi92W0RFYhDJgAS1+EC/8TAplORxRRdTOUrmlRi0crC9fQZJFEfnLVQP7zrVV89MUeXltcXL09M+AlM+gj6PcQ9HkJ+jy0zQzw1DeHkpcRcDDiGBl5J6x6HVa9YVUaCSicLLQZqjatLNxDk0US6ZiXztRbRwKwv7ScL3YdYcOuI2zaW0pZRRUnKkKcqAxRWl7JR+v38Pbyr/j22G7OBh0LXcZA4UBY8CcY+m1rmG6CCScLnUywtuo+C00WSU+TRZJqmxlgTI98xvSotw4UAOOe+IS3lrkkWYjAyNvhb9+HksVQNMLpiOrRPovItIPbPbSD26WuGtKRxVsPsH3/MadDiY3BN0AgCxY+G31fB/g9VnOLJovatBnKPTRZuNSVgzsCMGu5S9bGDmbDGROtfotj+52Oph7ts4hMKwv30GThUp3bZjCiaxtmLitJjmnOm2LEHVB1ApZOdzqSesLNUDqZYG3hyiKklUXS0z4LF5swtBM/fmsV63YeoX+HHKfDabn2A6x1M+Y9DftO8SK9vK5w7r/FNi5ONkNpZVGbVhbuocnCxS4f1IGfzlrNzGVfuSNZAJzzILz9wKnNblt+DE4cghG3Q0bbmIalHdyRaZ+Fe2iycLG2mQHO6d2Ot5d/xQ++1hePJ/GGnDZbr4vhe6tO7bUrX4c37oCju2OeLKqHzmqyqOXk0FmtLJKd9lm43IQhnSg5eJxFWw84HYrzMgus+9I9MT+0JovITq6Up/8uyU6ThctdMqA96X4vM5clyfKs8ZRVaN2X7o75oXU0VGTi9YLXq5WFC2iycLnMoI9LBrTnbyt3UF4ZcjocZ1VXFntjfujqPgu9grseCQS0g9sFNFmkgAlDOnLwWAVzNsS++SWppLcF8Vh9FjEWHg2lEwnWJ36/dnC7gCaLFHBO7wLyMvzMXOaSC/ROlccDGe3i0mfhEY+uw90ArSzcQUdDpYCAz8PlgzrwxpJi/v215bWey033c9+FvcnN8DsUXSvLKoxLsgCr30KboerTysIdNFmkiBtHd2Xuxr3M3Vi7vX73kRN8unEvf7ljNAXZKbAcaGZ8Kguw+i10NFR9EtBk4QaaLFLEgI45fPTvF9TbPmfDHia9sJjr/3ce0+8cTae8dAeia0WZhbB/flwO7ff4NVlEIH6/NkO5gPZZpLhzehcw/c5R7D16gm888xmb9hx1OqT4yiyIy2gosJqhtM+iPgkEtLJwAU0WiuFd2zJj0hhOVIa4/n/nsXbHYadDip+sAqgohfLSmB866A3qRIIRePzawe0G2gylABjYMZdX/2UsN/15Plf9/lOygtH/a+RlBHj1u2OTq6+j5lXcMV7PW0dDRaYd3O6gyUJV61mQxet3ncnzczdzIsoFfFUhw0sLtjH9861875I+rRRhDGTaV3Ef3QNtusX00NrBHZkEAoSOuWQRrhTW6slCRDoDLwDtAQNMMcZMFpG2wCtAN2ALcL0x5oCICDAZGA8cA241xixp7bhTRae8dP7j8gFN2nfnoTJenL+Vu87vSZrfG+fIYiSznXUfp/mhNFnUp5WFOzjRZ1EJfN8YMwAYA9wjIgOAh4H3jTG9gfftxwCXAb3t2yTgmdYPWUVy+9nd2Xu0PLlW44vz/FCaLOqzLsrTf5dk1+rJwhizI1wZGGOOAGuBTsAEYJq92zTgavvnCcALxvI5kCciHVo3ahXJmT3z6ds+m6mfbk6e1fjiOfOsRy/Ki8SqLLSDO9k5OhpKRLoBQ4H5QHtjzA77qZ1YzVRgJZLtNV5WbG+re6xJIrJIRBbt2ZPicyC1EhHh9rO7sW7nEeZt2ud0OE3jC0Iw1+qziDHts4hMAgFCWlkkPceShYhkAW8ADxhjao3VNNbX1GZ9VTXGTDHGjDDGjCgoKIhhpKoxE4Z0om1mgKmfbnE6lKbLKohLZeH36kV5kejcUO7gSLIQET9WonjRGPOmvXlXuHnJvg83KpcAnWu8vMjephJAmt/LTaO78P66XWzZG/trF+IiMz7JIugN6tDZCLQZyh1aPVnYo5ueBdYaY35X46lZwC32z7cAM2tsv1ksY4BDNZqrVAK4aUxXfB7h+c+2OB1K08QpWehEgpFpZeEOTlQWZwHfBi4UkWX2bTzwKHCJiGwALrYfA7wLbAI2An8C7nYgZtWIwpw0rhzckdcWbedwWRL8UcgsiMuaFgGPjoaKRIfOukOrX2dhjPkUkAaevijC/ga4J65BqRa77azuvLm0hFcXbufOc3o4HU7jsgrh+H6oqgRv7H4FwkNnjTFYBbQCq7KgshITCiEenWEoWeknp2JiUFEuo7q15bm5W6ioSvDlW8MX5h2L7YSCAW8Ag6EyVBnT4yY78VtrpWhTVHLTZKFi5q4LelJy8DiPzFqd2NddhKf8iHG/ha7DHZkEAgDaFJXkNFmomLmgbyF3n9+Tl+Zv44V5W50Op2HhC/Ni3G8RXodbR0TVJgGtLNxAk4WKqX+7tC+XDGjPz95Zw5wNCXpxZPWUH7FthqquLLSTu5bqZiitLJKaJgsVUx6P8PgNQ+hdmMU9Ly5JzMWUqicTjG1lEfBazS2aLGqrbobSyiKpabJQMZcV9PGnm0fg83q4c9oiDh1LsD8SwRzwBmPeDKXJIjKtLNxB17NQcdG5bQZ/vGk4N/75c257fgFjeuTH9XzZaX4uH9SBLvkZ0XcWicvyqgGPnSy0g7sW7eB2B00WKm5GdW/Lo9cM5pFZq1lZciiu56qoMjz2j3WM7ZHP9SOLGDewA+mBRtbYyCqIeTOU9llEpkNn3UGThYqra4cXce3worif56uDx3ljcTGvLS7me68s57+CqzmvbwHpDSzKdPvhINn7tzH5teUxi2F/1TYAnvxgLfneU7sor+9p2Yl/UWMzebSycAVNFsoVOualc99Fvbnngl4s2LKfVxdtZ/6m/Q1e73FeRQajQl8wd2PsmqKq/KVQAKu+2ovvRPOPW1pexWuLi7lhZGey0/wxi8tp2sHtDposlKt4PMKYHvnR+0hmfwLz5vDZwxdafRgxsG7/Or7x9hP88roBXNSl3sw1UX24fje3PbeQVSWHGdszvn08rUk7uN1BR0Op1JRVCKEKKDsYs0NWd3CfYp/FoE65AKyKc/9Oa9PKwh00WajUVL28auyaoVo6dLZdVpCOuWmscFuy0MrCFTRZqNQUhyk/wsmiJdN9DCrKdW1lEdJkkdQ0WajUVF1ZxG5KkvDQ2YrQqTe3DOqUy+a9pcmxLkgT6dBZd9BkoVJTVuxnno3FRIKDivIAd/Vb6EV57qDJQqWm9LaAxDRZxGK6j3An98piFyYLrSySmiYLlZq8PsjIj2mfhc/jwyveFiWLtpkBOuWlx/2K99Z0soNbk0Uy02ShUldmQcwXQAovrdoSg4ty3ZksKrQZKplpslCpKytOyaKFEwme3imXrfuOJd5svadIvF7werWySHKaLFTqikNlEfQEY1JZAKz6ykXVRSCgHdxJTpOFSl2ZhXA0tsnC7/W3OFmc3tFKFivc1Mnt92sHd5LTZKFSV2Y7KD8CFcdjdsiAN9DiNbjbZAbo3DbddcNntbJIbposVOqKw7UWQW8wJosfDe6Ux4qSgy0PKEFIQCuLZKfJQqWuOFzFHfC0fDQUWJ3c2/cf5+Axd3wbF79fK4skp8lCpa5Mu7KIYb9FLIbOwslObrcMofUEAlpZJDlNFip1Zbaz7mN8FXdMKouO7koWaGWR9DRZqNRV3QwV25lnY9FnkZvhp2t+hmum/fD4tYM72WmyUKkrkAGBrNiuaeFp+WiosNM7uedKbh06m/w0WajUllkQ0/mhgt6WX5QXNrhTLsUHjnOgNPm/kevQ2eSnyUKlthhfxR3wBqiois036EEu6uQW7eBOej6nA2gqERkHTAa8wJ+NMY86HJJyg6xC2L8pZocLeAOcCMWuGQqsZHFun4KYHNMp4vdTvm0bJQ8+GHVf32kdyL3yCoL9+iEirRCdaoqkSBYi4gWeBi4BioGFIjLLGLPG2chU0stsB9vnx+xw4essjlUca/GxfF7o1s7Hsu27OFbRMQbROcd/5ii8GzdQtnZd4zsaQ/nsf7J/6lSCffqQO2ECOVdcgb99YesEqhqUFMkCGAVsNMZsAhCRGcAEQJOFapnMQqsZ6smhMTlcerCK42mG0S+NjsnxKIB9Bka/FJvDOWoieIwhaCBoIC0EflN/t4zjhiHrheFrvqDbr3/Nzt/8mj1twGiR0SS7u2dzw/QFMT9usiSLTsD2Go+LgVq/jSIyCZgE0KVLl9aLTCW3QdfBoe0QqozJ4a4PlZNduY9QTI4GpeWV7C8tJ8Lf1FMm5uTRpPrIsTxDZFUYTkj4FqJcDJUS4bxpsOkM69bmQBUD1lSQvz9W/6LuZ/Jy43LcZEkWURljpgBTAEaMGBH///nKHQr6wtf/GLvDATfH7GhKJY5kGQ1VAnSu8bjI3qaUUqoVJEuyWAj0FpHuIhIAJgKzHI5JKaVSRlI0QxljKkXkXuD/sIbOTjXGrHY4LKWUShlJkSwAjDHvAu86HYdSSqWiZGmGUkop5SBNFkoppaLSZKGUUioqTRZKKaWiEmPcd/2aiOwBtrbgEO2A2C1y4By3vA/Q95KI3PI+QN9LWFdjTMRZK12ZLFpKRBYZY0Y4HUdLueV9gL6XROSW9wH6XppCm6GUUkpFpclCKaVUVJosIpvidAAx4pb3AfpeEpFb3gfoe4lK+yyUUkpFpZWFUkqpqDRZKKWUikqTRQ0iMk5E1ovIRhF52Ol4WkJEtojIShFZJiKLnI6nOURkqojsFpFVNba1FZHZIrLBvm/jZIxN0cD7+ImIlNifyzIRGe9kjE0lIp1F5EMRWSMiq0Xkfnt7Mn4uDb2XpPpsRCRNRBaIyHL7ffzU3t5dRObbf8desZd1aPn5tM/CIiJe4AvgEqxlWxcC3zTGJOU63yKyBRhhjEm6C41E5FzgKPCCMeZ0e9uvgP3GmEftRN7GGPOQk3FG08D7+Alw1BjzGydjay4R6QB0MMYsEZFsYDFwNXAryfe5NPRerieJPhsRESDTGHNURPzAp8D9wIPAm8aYGSLyR2C5MeaZlp5PK4uTRgEbjTGbjDHlwAxggsMxpSRjzCfA/jqbJwDT7J+nYf1yJ7QG3kdSMsbsMMYssX8+AqwFOpGcn0tD7yWpGMtR+6HfvhngQuB1e3vMPhNNFid1ArbXeFxMEv4HqsEA74nIYhGZ5HQwMdDeGLPD/nkn0N7JYFroXhFZYTdTJXyzTV0i0g0YCswnyT+XOu8FkuyzERGviCwDdgOzgS+Bg8aYSnuXmP0d02ThXmcbY4YBlwH32E0irmCsttNkbT99BugJDAF2AL91NJpmEpEs4A3gAWPM4ZrPJdvnEuG9JN1nY4ypMsYMAYqwWkf6xetcmixOKgE613hcZG9LSsaYEvt+N/BXrP9IyWyX3dYcbnPe7XA8p8QYs8v+BQ8BfyKJPhe7XfwN4EVjzJv25qT8XCK9l2T+bIwxB4EPgbFAnoiEV0GN2d8xTRYnLQR62yMJAsBEYJbDMZ0SEcm0O+4QkUzgUmBV469KeLOAW+yfbwFmOhjLKQv/YbV9nST5XOzO1GeBtcaY39V4Kuk+l4beS7J9NiJSICJ59s/pWINz1mIljevs3WL2mehoqBrsoXJPAF5gqjHm585GdGpEpAdWNQHWOusvJdN7EZGXgfOxplreBTwCvAW8CnTBmn7+emNMQnceN/A+zsdq5jDAFuC7Ndr8E5aInA3MAVYCIXvzj7Da+pPtc2novXyTJPpsRGQwVge2F+uL/6vGmJ/Zv/8zgLbAUuAmY8yJFp9Pk4VSSqlotBlKKaVUVJoslFJKRaXJQimlVFSaLJRSSkWlyUIppVRUmiyUqkNE8kTkbvvn80XknTic43kRuS76ntX7d6s5e22d5z4SkRGxi06p+jRZKFVfHnB3c15gz1qslGtpslCqvkeBnvYEbb8GskTkdRFZJyIv2lcAh9cMeUxElgDfEJFLRWSeiCwRkdfsuYcQkUfttRNWiEjN6a/PFZHPRGRTuMoQy69FZJVY65HcUDc4EUkXkRkislZE/gqkx/nfQyl80XdRKuU8DJxujBkiIudjTZcwEPgKmAuchbV2AMA+Y8wwEWkHvAlcbIwpFZGHgAdF5GmsqSP6GWNMeHoGWwfgbKzJ32ZhTSt9DdZVxGdgXfm9UEQ+qRPfXcAxY0x/+yreJbF880pFopWFUtEtMMYU2xPMLQO61XjuFft+DDAAmGtXJLcAXYFDQBnwrIhcAxyr8dq3jDEhe4Gt8NTeZwMv2xPa7QI+BkbWiedcYDqAMWYFsCIWb1KpxmhloVR0NefVqaL2702pfS/AbGPMN+u+WERGARdhTe52L9biNHWPKzGLVqk40MpCqfqOANnNfM3nwFki0guqZ/7tY/db5Bpj3gW+h9W81Jg5wA32ojYFWFXEgjr7fAJ8yz7P6cDgZsaqVLNpZaFUHcaYfSIy1x6qehxrxthor9kjIrcCL4tI0N78n1iJZ6aIpGFVDw9GOdRfsdYkWI41++kPjDE77RXdwp4BnhORtVhTUi9u8ptT6hTprLNKKaWi0mYopZRSUWmyUEopFZUmC6WUUlFpslBKKRWVJgullFJRabJQSikVlSYLpZRSUf0/+9VZALGwB9EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Creating +-1 label vector\n",
    "y = targets.float() * 2 - 1\n",
    "\n",
    "for key, phi in reps.items():\n",
    "    # Adding a constant 1 feature to account for the bias unit in a target linear model\n",
    "    bias_feature = torch.ones(phi.shape[0], 1).float()\n",
    "    phi_bias = torch.hstack((phi, bias_feature))\n",
    "\n",
    "    # Normalizing each row\n",
    "    phi_norm = phi_bias / torch.linalg.norm(phi_bias, dim=1)[:, None]\n",
    "\n",
    "    u, s, vt = torch.linalg.svd(phi_norm)\n",
    "    dot_products = (u.T @ y)**2\n",
    "    \n",
    "    grid = np.linspace(0, 30, (30+1)*2)[::-1]\n",
    "    alignment_curve = get_alignment_curve(\n",
    "        dot_products.detach().numpy(), y.detach().numpy(), s.detach().numpy(), grid)[::-1]\n",
    "\n",
    "    plt.plot(grid[::-1], alignment_curve, label=key)\n",
    "plt.legend()\n",
    "plt.ylim(ymin=-50)\n",
    "plt.ylabel('Alignment')\n",
    "plt.xlabel('threshold')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}