predictive_coding.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zLygbz278oF1"
      },
      "source": [
        "# Supervised learning using predictive coding\n",
        "\n",
        "This notebook illustrates how the predictive coding library can be used to train a predictive coding network on a supervised learning task (MNIST)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tOeSG__s8oF3"
      },
      "source": [
        "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/Bogacz-Group/PredictiveCoding/blob/main/1_supervised_learning_pc.ipynb)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fJO9KSFm8oF4",
        "outputId": "a16b6ec8-b53d-4dab-88ba-f3510c43d014"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "using cpu\n"
          ]
        }
      ],
      "source": [
        "import torch\n",
        "import torch.nn as nn\n",
        "import torch.optim as optim\n",
        "from torchvision import datasets, transforms\n",
        "import torch.nn.functional as F\n",
        "import random\n",
        "import numpy as np\n",
        "import os\n",
        "import matplotlib.pyplot as plt\n",
        "from tqdm import tqdm\n",
        "import copy\n",
        "\n",
        "import predictive_coding as pc\n",
        "\n",
        "torch.manual_seed(42)\n",
        "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
        "print(f'using {device}')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "g72NMl0x8oF5",
        "outputId": "94fbdffc-644c-4dd2-f020-1e7a181226db"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "# train images: 60000 and # test images: 10000\n"
          ]
        }
      ],
      "source": [
        "# load data\n",
        "transform = transforms.Compose([transforms.ToTensor(), transforms.Lambda(lambda x: torch.flatten(x))])\n",
        "train_dataset = datasets.MNIST('./data', download=True, train=True, transform=transform)\n",
        "test_dataset = datasets.MNIST('./data', download=True, train=False, transform=transform)\n",
        "\n",
        "batch_size = 500\n",
        "train_loader = torch.utils.data.DataLoader(train_dataset, batch_size=batch_size, shuffle=True)\n",
        "\n",
        "print(f'# train images: {len(train_dataset)} and # test images: {len(test_dataset)}')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZveDwPLX8oF5"
      },
      "source": [
        "### Defining a model\n",
        "A model can be initalised in the same way as pytorch model, with the addition of `pc.PCLayer()` to include latent variables in the model.\n",
        "\n",
        "A `PCLayer()` contains the activities of a layer of latent variables under `pclayer._x`. A `PCLayer()` also contains the energy associated with that activity under `pclayer._energy` which is computed with `0.5 *(inputs['mu'] - inputs['x'])**2` where `inputs['x']` is the activity of that layer and `inputs['mu']` is the input to that layer.\n",
        "\n",
        "Check out the `PCLayer()` class in `predictive_coding/pc_layer.py` for more information."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KbcZEBnH8oF5",
        "outputId": "7fdc4413-3532-434d-b2eb-3e7ca41d727b"
      },
      "outputs": [
        {
          "data": {
            "application/vnd.code.notebook.error": {
              "message": "module 'predictive_coding' has no attribute 'PCLayer'",
              "name": "AttributeError",
              "stack": "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)\nCell \u001b[0;32mIn[4], line 10\u001b[0m\n\u001b[1;32m      4\u001b[0m activation_fn \u001b[38;5;241m=\u001b[39m nn\u001b[38;5;241m.\u001b[39mReLU\n\u001b[1;32m      5\u001b[0m loss_fn \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mlambda\u001b[39;00m output, _target: \u001b[38;5;241m0.5\u001b[39m \u001b[38;5;241m*\u001b[39m (output \u001b[38;5;241m-\u001b[39m _target)\u001b[38;5;241m.\u001b[39mpow(\u001b[38;5;241m2\u001b[39m)\u001b[38;5;241m.\u001b[39msum() \u001b[38;5;66;03m# this loss function holds to the error of the output layer of the model\u001b[39;00m\n\u001b[1;32m      8\u001b[0m model \u001b[38;5;241m=\u001b[39m nn\u001b[38;5;241m.\u001b[39mSequential(\n\u001b[1;32m      9\u001b[0m     nn\u001b[38;5;241m.\u001b[39mLinear(input_size, hidden_size),\n\u001b[0;32m---> 10\u001b[0m     pc\u001b[38;5;241m.\u001b[39mPCLayer(),\n\u001b[1;32m     11\u001b[0m     activation_fn(),\n\u001b[1;32m     12\u001b[0m     nn\u001b[38;5;241m.\u001b[39mLinear(hidden_size, hidden_size),\n\u001b[1;32m     13\u001b[0m     pc\u001b[38;5;241m.\u001b[39mPCLayer(),\n\u001b[1;32m     14\u001b[0m     activation_fn(),\n\u001b[1;32m     15\u001b[0m     nn\u001b[38;5;241m.\u001b[39mLinear(hidden_size, output_size)\n\u001b[1;32m     16\u001b[0m )\n\u001b[1;32m     17\u001b[0m model\u001b[38;5;241m.\u001b[39mtrain()   \u001b[38;5;66;03m# set the model to training mode\u001b[39;00m\n\u001b[1;32m     18\u001b[0m model\u001b[38;5;241m.\u001b[39mto(device)\n\n\u001b[0;31mAttributeError\u001b[0m: module 'predictive_coding' has no attribute 'PCLayer'"
            }
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "input_size = 28*28  # 28x28 images\n",
        "hidden_size = 256\n",
        "output_size = 10    # 10 classes\n",
        "activation_fn = nn.ReLU\n",
        "loss_fn = lambda output, _target: 0.5 * (output - _target).pow(2).sum() # this loss function holds to the error of the output layer of the model\n",
        "\n",
        "\n",
        "model = nn.Sequential(\n",
        "    nn.Linear(input_size, hidden_size),\n",
        "    pc.PCLayer(),\n",
        "    activation_fn(),\n",
        "    nn.Linear(hidden_size, hidden_size),\n",
        "    pc.PCLayer(),\n",
        "    activation_fn(),\n",
        "    nn.Linear(hidden_size, output_size)\n",
        ")\n",
        "model.train()   # set the model to training mode\n",
        "model.to(device)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YGR08KhU8oF5"
      },
      "source": [
        "### Defining a model trainer\n",
        "The predictive coding library is based around a `pc.PCTrainer()`.\n",
        "\n",
        "This trainer orchestrate the activity and parameter updates of the model to minimise the total error of the model. The total error is given by the sum of the energies in each pclayer as well as the loss functions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "zrTV_m728oF6"
      },
      "outputs": [],
      "source": [
        "# number of inference iterations where the latent states x are updated. Inference does not run till convergence but for a fixed number of iterations\n",
        "T = 20\n",
        "\n",
        "# options for the update of the latent state x\n",
        "optimizer_x_fn = optim.SGD          # optimizer for latent state x, SGD perform gradient descent. Other alternative are Adam, RMSprop, etc.\n",
        "optimizer_x_kwargs = {'lr': 0.01}   # optimizer parameters for latent state x to pass to the optimizer. The best learning rate will depend on the task and the optimiser.\n",
        "                                    # Other parameters such as momentum, weight_decay could also be set here with additional elements, e.g., \"momentum\": 0.9, \"weight_decay\": 0.01\n",
        "\n",
        "# options for the update of the parameters p\n",
        "update_p_at = 'last'                # update parameters p at the last iteration, can be set to 'all' to implement ipc (https://arxiv.org/abs/2212.00720)\n",
        "optimizer_p_fn = optim.Adam         # optimizer for parameters p\n",
        "optimizer_p_kwargs = {'lr': 0.001}  # optimizer parameters for parameters p, 0.001 is a good starting point for Adam, but it should be adjusted for the task\n",
        "\n",
        "trainer = pc.PCTrainer(model,\n",
        "    T = T,\n",
        "    optimizer_x_fn = optimizer_x_fn,\n",
        "    optimizer_x_kwargs = optimizer_x_kwargs,\n",
        "    update_p_at = update_p_at,\n",
        "    optimizer_p_fn = optimizer_p_fn,\n",
        "    optimizer_p_kwargs = optimizer_p_kwargs,\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "id": "thtspenN8oF6"
      },
      "outputs": [],
      "source": [
        "# get classification accuracy of the model\n",
        "def test(model, dataset, batch_size=1000):\n",
        "    model.eval()\n",
        "    test_loader = torch.utils.data.DataLoader(dataset, batch_size=batch_size, shuffle=False)\n",
        "    correct = 0\n",
        "    total = 0\n",
        "    for data, label in test_loader:\n",
        "        data, label = data.to(device), label.to(device)\n",
        "        pred = model(data)\n",
        "        _, predicted = torch.max(pred, -1)\n",
        "        total += label.size(0)\n",
        "        correct += (predicted == label).sum().item()\n",
        "    model.train()\n",
        "    return round(correct / total, 4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q7c9AIGY8oF6"
      },
      "source": [
        "### Train the model\n",
        "`trainer.train_on_batch()` is called for each batch of data. This function updates the activity of the latent states and the parameters for the given batch of data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 623
        },
        "id": "A0y1zzbP8oF6",
        "outputId": "48b3e7df-1ed6-4793-b724-7ea4ef5f3686"
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Epoch 1 - Test accuracy: 0.103: 100%|██████████| 120/120 [00:17<00:00,  6.70it/s]\n",
            "Epoch 2 - Test accuracy: 0.959: 100%|██████████| 120/120 [00:18<00:00,  6.53it/s]\n",
            "Epoch 3 - Test accuracy: 0.971: 100%|██████████| 120/120 [00:16<00:00,  7.27it/s]\n",
            "Epoch 4 - Test accuracy: 0.975: 100%|██████████| 120/120 [00:14<00:00,  8.55it/s]\n",
            "Epoch 5 - Test accuracy: 0.979: 100%|██████████| 120/120 [00:13<00:00,  8.91it/s]\n",
            "Epoch 6 - Test accuracy: 0.982: 100%|██████████| 120/120 [00:13<00:00,  8.84it/s]\n",
            "Epoch 7 - Test accuracy: 0.984: 100%|██████████| 120/120 [00:13<00:00,  8.77it/s]\n",
            "Epoch 8 - Test accuracy: 0.983: 100%|██████████| 120/120 [00:13<00:00,  9.00it/s]\n",
            "Epoch 9 - Test accuracy: 0.983: 100%|██████████| 120/120 [00:14<00:00,  8.24it/s]\n",
            "Epoch 10 - Test accuracy: 0.983: 100%|██████████| 120/120 [00:13<00:00,  8.67it/s]\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1RElEQVR4nO3de3xTdZ7/8XeStmkLbSnW3qBYVEZBARkQBtGZh9oR0WUXdQSRlQrj+NBBRLrOKqOAqFBxFkSFhRFR9yKKOt52cHC0oq6KgmAVf+tdaCvQlmtTWnpLzu+PNmlDCzSQ5NBzXs/HIw+Sk5P00zBD3n6/3/P9OAzDMAQAAGARTrMLAAAACCfCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsJQYswuINp/Pp507dyopKUkOh8PscgAAQCcYhqHq6mplZ2fL6Tz62Iztws3OnTuVk5NjdhkAAOA4lJWVqXfv3kc9x3bhJikpSVLzh5OcnGxyNQAAoDM8Ho9ycnIC3+NHY7tw45+KSk5OJtwAANDFdGZJCQuKAQCApRBuAACApRBuAACApRBuAACApZgabt5//32NHTtW2dnZcjgcevXVV4/5mnfffVc///nP5Xa7deaZZ+qZZ56JeJ0AAKDrMDXc1NTUaPDgwVq2bFmnzt+2bZuuvPJKXXzxxSouLtYdd9yhm266SW+++WaEKwUAAF2FqZeCjxkzRmPGjOn0+StWrFDfvn21aNEiSVL//v31wQcf6JFHHtHo0aMjVSYAAOhCutSamw0bNigvLy/o2OjRo7Vhw4Yjvqa+vl4ejyfoBgAArKtLhZvy8nJlZGQEHcvIyJDH49GhQ4c6fE1hYaFSUlICN1ovAABgbV0q3ByPWbNmqaqqKnArKyszuyQAABBBXar9QmZmpioqKoKOVVRUKDk5WQkJCR2+xu12y+12R6M8AABwEuhSIzcjR45UUVFR0LG33npLI0eONKkiAABwsjF15ObgwYP6/vvvA4+3bdum4uJi9ezZU3369NGsWbO0Y8cO/ed//qck6ZZbbtHSpUv1r//6r5o6dareeecdvfDCC1q7dq1ZvwIAdIphGGryGWryGmry+Vr+bL7v9RmB8/xNAdu2BvT3CXS0HG193PakY59z+Hsf/r5B73PY69s2K+zo9W17GXZ03HGEn9+ZJog4NsMwZBiS4b8fOC4Zan6u7ePW+63nG4EXtb7mSO+nwHP+lwSfHxfjVHpSfOR/8SMwNdx8+umnuvjiiwOPCwoKJEn5+fl65plntGvXLpWWlgae79u3r9auXauZM2fq0UcfVe/evfXkk09yGThgUf5A0NDkU6PXp4Ymnxpa/mz0Goc9bv3THxoavYa8PkNN/mP+QOH1qdFnyNs2ZLQ5p7ElcHQURDp+T1/HwSVQixEUYNCxIwWhdiHMcVhAahPo2oaxQGxqc/xw7f5WOvhrOvyQYbQ/qf05Hb2P0Ylz2h84Zhg5Cf28Tw+9/PtRpv18h9HR35KFeTwepaSkqKqqSsnJyWaXA5jOMIzmoOAPDy0Bob7N/QavT41NPtW3/NngbRs2jHbhouPQYaj+8PPanNt4hNBidTFOh5zO5i/ioH+M2/wXsRT8Rdb6X9FG0OPDzwMiwR8i2wbIQMhseW5Inx56/ubwLhkJ5fu7Sy0oBrq6YwWJtscbvF41NHX8nD8I1Lc9v93rDwsohz/f5vVdhdMhxcU4Fetyyt3yp/9xnMup2Bin4lwOxTidinE5FON0yOV0KtblkMvpUKzLqRino+U5Z8ux9ue4nM2v9d8//JwYZ/D7xLhazwk6ftg5zY9bjjub3ysa0zKBEBQITB08F3jc5rkOgpXanHf4iIL//YLe67BpC//PavtzjcOCnIzgeo74c9rW2smf09HHffihjv5O2p/T0fs4jnlOu9d09D6HHfSPXrUboXK0Tv0dPmIVFDaO8vzh04lt36/tuR3VdTIj3MDW/GHjUKNX9Y1e1TX6dKjRq7qW26GWY/VNXh1qaDne5Gu+3+RVfWPr/ebzfapreS//a/3v4w8aJzuX06G4ltAQF+Nsve9yKjam+blYV/BzbUNGc+hwBB2P6zCEOBUX41CcyxV0frvQ0ub1LmfX+cf1ZBKY4unw4+MzhfUQbtAlGIYhT12TKj11qvDUq+pQY5vw0RwcWsOHV4cafC3ho4OQ0eZ+XaNXZi6F6ChIuGPaB4vYNue4Owoeh70+9gjnuGOcinO52r2m7TkECABdHeEGpjvU4FWFp04VnjqVe+pU6alvflxdHzhe4alTXWNkRz2cDik+1qWEWJfiY11yxzoD9+NjnYqPcSk+ztX85+HPBe77X9/2WPN9d4xT7hgXQQIAIoxwg4hpaPJp98GWgFJVFxRYKj31Km8JLdV1TZ1+z5SEWGUku9UjMa4lSLSEiBiXEuKaA4n/fnxM83MJcS652wWS1hDibrkf64rO+gcAQGQRbhAyr8/Q3pr65oBSVaeK6uaposrA6Evz/b01DZ1+z4RYlzJT4pWe5FZmSrwykpvvZyQ3389Mjld6slvxsa4I/mYAACsg3CDAMAxVHWpURZtRFf8al9bpoXrtPljf6T07Yl0OpSfFKyO5Nag039xBf3Z3xzBqAgAIC8KNjRmGoUV//1Yf/7g3MPrS2at5nA4prXtrQElvGV3x389oCTSpiXFysq4EABBFhBsb+3FPjZau/77d8dTE2OZpoeR4ZfinhlLa3E+OV1r3OMW4ulRrMgCATRBubGzb7hpJUt+0bvrTbwYpIzlepyaxrgUA0LURbmysZF+tJKl/VpKG5fY0uRoAAMKDeQUbK93bPHLTp2c3kysBACB8CDc2tn1v88jNaackmlwJAADhQ7ixsdJ9hBsAgPUQbmyqyevTT/v94YZpKQCAdRBubGpXVZ0avYbiXE5lJsebXQ4AAGFDuLGpkpb1Njk9E2jeCACwFMKNTZXsa75SiikpAIDVEG5syj9y06cni4kBANZCuLGpkpY9bnK5UgoAYDGEG5sq2cuVUgAAayLc2JBhGIE9bvowcgMAsBjCjQ3tPliv2gavnA6pd2qC2eUAABBWhBsb8k9JZaUkyB1DB3AAgLUQbmyohJ5SAAALI9zYkL8bOIuJAQBWRLixIbqBAwCsjHBjQyUtV0qxxw0AwIoINzbkn5bq05NpKQCA9RBubKbqUKP21zZKYo8bAIA1EW5sprRlvU1ad7e6u2NMrgYAgPAj3NhMazdwRm0AANZEuLGZwB43dAMHAFgU4cZmStjjBgBgcYQbm2F3YgCA1RFubMYfbrhSCgBgVYQbG6lr9KrcUydJymVaCgBgUYQbGylr2Zk4yR2j1MRYk6sBACAyCDc2sr3NlJTD4TC5GgAAIoNwYyP+K6WYkgIAWBnhxkZK97GYGABgfYQbG9nOBn4AABsg3NhIKRv4AQBsgHBjE01en37af0gSG/gBAKyNcGMTOw/UqclnKC7GqczkeLPLAQAgYgg3NuHvBt6nZ6KcTi4DBwBYF+HGJlhMDACwC8KNTfgXE3MZOADA6gg3NuFvmMkGfgAAqyPc2ATdwAEAdkG4sQHDMAK7EzNyAwCwOsKNDeyurtehRq+cDqlXjwSzywEAIKIINzbgv1Iqu0eC4mL4KwcAWBvfdDZAN3AAgJ0QbmyAbuAAADsh3NgAG/gBAOyEcGMDdAMHANgJ4cYGSlqmpegGDgCwA8KNxVXVNupAbaOk5qaZAABYHeHG4vzdwE9NcqubO8bkagAAiDzCjcWVsJgYAGAzhBuLK6EbOADAZgg3Fkc3cACA3RBuLI4rpQAAdkO4sbjAtBRrbgAANmF6uFm2bJlyc3MVHx+vESNGaOPGjUc9f8mSJTrrrLOUkJCgnJwczZw5U3V1dVGqtmupa/SqwlMviWkpAIB9mBpu1qxZo4KCAs2dO1dbtmzR4MGDNXr0aFVWVnZ4/urVq3X33Xdr7ty5+uqrr7Rq1SqtWbNGf/zjH6Ncedfg7ymVFB+jHomxJlcDAEB0mBpuFi9erN/97neaMmWKBgwYoBUrVigxMVFPPfVUh+d/9NFHGjVqlK6//nrl5ubqsssu08SJE4852mNX2/f42y4kyuFwmFwNAADRYVq4aWho0ObNm5WXl9dajNOpvLw8bdiwocPXXHDBBdq8eXMgzPz444964403dMUVVxzx59TX18vj8QTd7KI0sJiYKSkAgH2YtmXtnj175PV6lZGREXQ8IyNDX3/9dYevuf7667Vnzx5deOGFMgxDTU1NuuWWW446LVVYWKh58+aFtfaugg38AAB2ZPqC4lC8++67WrBggf793/9dW7Zs0csvv6y1a9fqgQceOOJrZs2apaqqqsCtrKwsihWba3vLlVIsJgYA2IlpIzdpaWlyuVyqqKgIOl5RUaHMzMwOXzN79mzdcMMNuummmyRJAwcOVE1NjW6++Wbdc889cjrbZzW32y232x3+X6AL8E9LsTsxAMBOTBu5iYuL09ChQ1VUVBQ45vP5VFRUpJEjR3b4mtra2nYBxuVySZIMw4hcsV1Qo9enn/YfksQGfgAAezG1TXRBQYHy8/M1bNgwDR8+XEuWLFFNTY2mTJkiSZo8ebJ69eqlwsJCSdLYsWO1ePFiDRkyRCNGjND333+v2bNna+zYsYGQg2Y7DxyS12fIHeNURlK82eUAABA1poabCRMmaPfu3ZozZ47Ky8t13nnnad26dYFFxqWlpUEjNffee68cDofuvfde7dixQ6eeeqrGjh2r+fPnm/UrnLT8i4n79EyU08ll4AAA+3AYNpvP8Xg8SklJUVVVlZKTk80uJ2L+a8N2zX7t/ymvf7qezD/f7HIAADghoXx/d6mrpdB5gcvAuVIKAGAzhBuLohs4AMCuCDcWRTdwAIBdEW4syDCMwB43bOAHALAbwo0FVVbXq67RJ5fToV6pCWaXAwBAVBFuLMjfDTy7R7xiXfwVAwDshW8+CyphSgoAYGOEGwsqbbOBHwAAdkO4sSB/N3AuAwcA2BHhxoJK97GBHwDAvgg3FtS6OzEjNwAA+yHcWMyB2gZVHWqUxJobAIA9EW4sxj9qk57kVmKcqU3fAQAwBeHGYugpBQCwO8KNxZTs8feUYjExAMCeCDcW07qBHyM3AAB7ItxYTGADP8INAMCmCDcW49/Aj9YLAAC7ItxYSG1Dkyqr6yWxoBgAYF+EGwvx70ycHB+jHolxJlcDAIA5CDcW4t/jJjeNKSkAgH0RbiyEbuAAABBuLIVu4AAAEG4shW7gAAAQbiwl0A2caSkAgI0Rbiyi0evTjgOHJDFyAwCwN8KNRezYf0hen6H4WKfSk9xmlwMAgGkINxbh7ynVp2einE6HydUAAGAewo1FlOylGzgAABLhxjICG/hxGTgAwOYINxYRuFKKcAMAsDnCjUUEpqW4UgoAYHOEGwvw+YzABn5MSwEA7I5wYwGV1fWqb/LJ5XQou0eC2eUAAGAqwo0F+HtK9eqRoFgXf6UAAHvjm9ACSllMDABAAOHGAkr20Q0cAAA/wo0FbA/sccOVUgAAEG4swD8t1Ydu4AAAEG66OsMwAguK6QYOAADhpss7UNuo6romSYzcAAAgEW66PH838IxktxLiXCZXAwCA+Qg3XZy/7cJpdAMHAEAS4abLo2EmAADBCDddHOEGAIBgIYebuXPnqqSkJBK14DjQDRwAgGAhh5vXXntNZ5xxhi699FKtXr1a9fX1kagLnVRCN3AAAIKEHG6Ki4u1adMmnXPOOZoxY4YyMzN16623atOmTZGoD0dR29Ck3dXN4ZIFxQAANDuuNTdDhgzRY489pp07d2rVqlX66aefNGrUKA0aNEiPPvqoqqqqwl0nOuBfb5OSEKuUxFiTqwEA4ORwQguKDcNQY2OjGhoaZBiGUlNTtXTpUuXk5GjNmjXhqhFHULKXKSkAAA53XOFm8+bNuu2225SVlaWZM2dqyJAh+uqrr/Tee+/pu+++0/z583X77beHu1YcpnQfi4kBADhcyOFm4MCB+sUvfqFt27Zp1apVKisr00MPPaQzzzwzcM7EiRO1e/fusBaK9vzdwE+j7QIAAAExob5g/Pjxmjp1qnr16nXEc9LS0uTz+U6oMBxbKXvcAADQTsjhZvbs2ZGoA8ehZB/dwAEAOFzI01LXXHONFi5c2O74ww8/rGuvvTYsReHYGpp82rH/kCRGbgAAaCvkcPP+++/riiuuaHd8zJgxev/998NSFI5tx4FD8hlSfKxT6Ulus8sBAOCkEXK4OXjwoOLi4todj42NlcfjCUtROLa23cAdDofJ1QAAcPI4rqulOtrD5vnnn9eAAQPCUhSOjYaZAAB07LgWFF999dX64YcfdMkll0iSioqK9Nxzz+nFF18Me4HoGOEGAICOhRxuxo4dq1dffVULFizQSy+9pISEBA0aNEhvv/22fvWrX0WiRnSADfwAAOhYyOFGkq688kpdeeWV4a4FIdhO6wUAADp0Qr2lYA6fz1DpPv/uxIzcAADQVsgjN16vV4888oheeOEFlZaWqqGhIej5ffv2ha04dKyiuk4NTT7FOB3K7hFvdjkAAJxUQh65mTdvnhYvXqwJEyaoqqpKBQUFuvrqq+V0OnXfffdFoEQcbvue5lGb3qkJinEx+AYAQFshfzM+++yzWrlypf7lX/5FMTExmjhxop588knNmTNHH3/8cSRqxGFYTAwAwJGFHG7Ky8s1cOBASVL37t1VVVUlSfqHf/gHrV27NrzVoUMldAMHAOCIQg43vXv31q5duyRJZ5xxhv7+979LkjZt2iS3O/Q2AMuWLVNubq7i4+M1YsQIbdy48ajnHzhwQNOmTVNWVpbcbrd+9rOf6Y033gj553Zl7HEDAMCRhRxurrrqKhUVFUmSpk+frtmzZ6tfv36aPHmypk6dGtJ7rVmzRgUFBZo7d662bNmiwYMHa/To0aqsrOzw/IaGBv3617/W9u3b9dJLL+mbb77RypUr1atXr1B/jS6NbuAAAByZwzAM40Te4OOPP9ZHH32kfv36aezYsSG9dsSIETr//PO1dOlSSZLP51NOTo6mT5+uu+++u935K1as0J/+9Cd9/fXXio2N7dTPqK+vV319feCxx+NRTk6OqqqqlJycHFK9JwPDMDTovr+rur5Jf5/5S/0sI8nskgAAiDiPx6OUlJROfX+HNHLT2NioqVOnatu2bYFjv/jFL1RQUBBysGloaNDmzZuVl5fXWozTqby8PG3YsKHD17z++usaOXKkpk2bpoyMDJ177rlasGCBvF7vEX9OYWGhUlJSArecnJyQ6jzZ7K9tVHV9kySpD2tuAABoJ6RwExsbq7/85S9h+cF79uyR1+tVRkZG0PGMjAyVl5d3+Joff/xRL730krxer9544w3Nnj1bixYt0oMPPnjEnzNr1ixVVVUFbmVlZWGp3yz+buCZyfGKj3WZXA0AACefkNfcjBs3Tq+++moESjk2n8+n9PR0PfHEExo6dKgmTJige+65RytWrDjia9xut5KTk4NuXZl/MXEfFhMDANChkHco7tevn+6//359+OGHGjp0qLp1C17Uevvtt3fqfdLS0uRyuVRRURF0vKKiQpmZmR2+JisrS7GxsXK5Wkcs+vfvr/LycjU0NCguLi7E36brKaGnFAAARxVyuFm1apV69OihzZs3a/PmzUHPORyOToebuLg4DR06VEVFRRo3bpyk5pGZoqIi3XbbbR2+ZtSoUVq9erV8Pp+czuZBp2+//VZZWVm2CDYSV0oBAHAsIYebtouJT1RBQYHy8/M1bNgwDR8+XEuWLFFNTY2mTJkiSZo8ebJ69eqlwsJCSdKtt96qpUuXasaMGZo+fbq+++47LViwoNOBygoC01IsJgYAoEMhh5twmjBhgnbv3q05c+aovLxc5513ntatWxdYZFxaWhoYoZGknJwcvfnmm5o5c6YGDRqkXr16acaMGbrrrrvM+hWirnVaipEbAAA6EvI+N8faqO+pp546oYIiLZTr5E82NfVNOmfum5Kkz+deppSEzu31AwBAVxfK93fIIzf79+8PetzY2Kgvv/xSBw4c0CWXXBLq2yEE/lGbHomxBBsAAI4g5HDzyiuvtDvm8/l066236owzzghLUehYKYuJAQA4ppD3uenwTZxOFRQU6JFHHgnH2+EI6AYOAMCxhSXcSNIPP/ygpqamcL0dOrCdPW4AADimkKelCgoKgh4bhqFdu3Zp7dq1ys/PD1thaM8/LdWHaSkAAI4o5HDz2WefBT12Op069dRTtWjRomNeSYUTE5iWYuQGAIAjCjncrF+/PhJ14BgamnzaeeCQJMINAABHE/Kam23btum7775rd/y7777T9u3bw1ETOvDT/lr5DCkxzqVTu7vNLgcAgJNWyOHmxhtv1EcffdTu+CeffKIbb7wxHDWhAyX7WtsuOBwOk6sBAODkFXK4+eyzzzRq1Kh2x3/xi1+ouLg4HDWhAyV7/HvcMCUFAMDRhBxuHA6Hqqur2x2vqqqS1+sNS1Fozz9ywwZ+AAAcXcjh5pe//KUKCwuDgozX61VhYaEuvPDCsBaHVqV0AwcAoFNCvlpq4cKF+uUvf6mzzjpLF110kSTpf//3f+XxePTOO++EvUA02763eVqKbuAAABxdyCM3AwYM0BdffKHx48ersrJS1dXVmjx5sr7++mude+65kajR9nw+Q2X7uQwcAIDOCHnkRpKys7O1YMGCcNeCIyj31KmhyacYp0NZKfFmlwMAwEkt5JGbp59+Wi+++GK74y+++KL+4z/+IyxFIZh/SiqnZ6JiXGFrBwYAgCWF/E1ZWFiotLS0dsfT09MZzYkQFhMDANB5IYeb0tJS9e3bt93x0047TaWlpWEpCsG201MKAIBOCzncpKen64svvmh3/PPPP9cpp5wSlqIQzN8NnD1uAAA4tpDDzcSJE3X77bdr/fr18nq98nq9eueddzRjxgxdd911kajR9gLdwJmWAgDgmEK+WuqBBx7Q9u3bdemllyompvnlPp9PkydP1vz588NeoN0ZhtEabpiWAgDgmEION3FxcVqzZo0efPBBFRcXKyEhQQMHDtRpp50Wifpsb19Ngw7WN8nhaL5aCgAAHN1x7XMjSf369VO/fv0kSR6PR8uXL9eqVav06aefhq04tPaUykyOV3ysy+RqAAA4+R13uJGk9evX66mnntLLL7+slJQUXXXVVeGqCy1K9tINHACAUIQcbnbs2KFnnnlGTz/9tA4cOKD9+/dr9erVGj9+vBwORyRqtLXWxcRcKQUAQGd0+mqpv/zlL7riiit01llnqbi4WIsWLdLOnTvldDo1cOBAgk2EBDbwY+QGAIBO6fTIzYQJE3TXXXdpzZo1SkpKimRNaINu4AAAhKbTIze//e1vtWzZMl1++eVasWKF9u/fH8m60KJ0H5eBAwAQik6Hmz//+c/atWuXbr75Zj333HPKysrSP/3TP8kwDPl8vkjWaFsH65u052CDJKalAADorJB2KE5ISFB+fr7ee+89bd26Veecc44yMjI0atQoXX/99Xr55ZcjVact+a+U6tktTsnxsSZXAwBA1xBy+wW/fv36acGCBSorK9N///d/q7a2VhMnTgxnbbZHN3AAAEJ3QvvcSJLT6dTYsWM1duxYVVZWhqMmtChhvQ0AACE77pGbjqSnp4fz7WyvdQM/rpQCAKCzwhpuEF50AwcAIHSEm5MY3cABAAgd4eYkVd/k1c6qQ5KYlgIAIBQhh5vTTz9de/fubXf8wIEDOv3008NSFKSf9h+SYUiJcS6ldY8zuxwAALqMkMPN9u3b5fV62x2vr6/Xjh07wlIUgi8Dp28XAACd1+lLwV9//fXA/TfffFMpKSmBx16vV0VFRcrNzQ1rcXZGTykAAI5Pp8PNuHHjJEkOh0P5+flBz8XGxio3N1eLFi0Ka3F2xmJiAACOT6fDjb9/VN++fbVp0yalpaVFrCi07nFDTykAAEIT8g7F27Zta3fswIED6tGjRzjqQQv/7sRMSwEAEJqQFxQvXLhQa9asCTy+9tpr1bNnT/Xq1Uuff/55WIuzK6/P0E/7mi8Dp68UAAChCTncrFixQjk5OZKkt956S2+//bbWrVunMWPG6A9/+EPYC7SjXVWH1OD1KdblUHaPBLPLAQCgSwl5Wqq8vDwQbv76179q/Pjxuuyyy5Sbm6sRI0aEvUA78l8GnpOaKJeTy8ABAAhFyCM3qampKisrkyStW7dOeXl5kiTDMDrc/wah86+3YTExAAChC3nk5uqrr9b111+vfv36ae/evRozZowk6bPPPtOZZ54Z9gLtiD1uAAA4fiGHm0ceeUS5ubkqKyvTww8/rO7du0uSdu3apd///vdhL9CO2u5ODAAAQhNyuImNjdWdd97Z7vjMmTPDUhDYwA8AgBNxXF3B/+u//ksXXnihsrOzVVJSIklasmSJXnvttbAWZ0eGYQQ28KMbOAAAoQs53CxfvlwFBQUaM2aMDhw4EFhE3KNHDy1ZsiTc9dnO3poG1TR45XBIOT25DBwAgFCFHG4ef/xxrVy5Uvfcc49cLlfg+LBhw7R169awFmdH/imprOR4uWNcxzgbAAAcLuRws23bNg0ZMqTdcbfbrZqamrAUZWdMSQEAcGJCDjd9+/ZVcXFxu+Pr1q1T//79w1GTrbGYGACAE9Ppq6Xuv/9+3XnnnSooKNC0adNUV1cnwzC0ceNGPffccyosLNSTTz4ZyVptoZQN/AAAOCGdDjfz5s3TLbfcoptuukkJCQm69957VVtbq+uvv17Z2dl69NFHdd1110WyVltgAz8AAE5Mp8ONYRiB+5MmTdKkSZNUW1urgwcPKj09PSLF2REb+AEAcGJC2sTP4Qhu4piYmKjERL6Ew6W6rlF7axokseYGAIDjFVK4+dnPftYu4Bxu3759J1SQnfkXE5/SLU5J8bEmVwMAQNcUUriZN2+eUlJSIlWL7bGYGACAExdSuLnuuutYXxNBgcvAWW8DAMBx6/Q+N8eajsKJYwM/AABOXKfDTdurpRAZbOAHAMCJ6/S0lM/ni2QdUOuaG8INAADHL+T2C4iM+iavdlYdksS0FAAAJ+KkCDfLli1Tbm6u4uPjNWLECG3cuLFTr3v++eflcDg0bty4yBYYBWX7DskwpG5xLp3SLc7scgAA6LJMDzdr1qxRQUGB5s6dqy1btmjw4MEaPXq0Kisrj/q67du3684779RFF10UpUojq+1iYhZvAwBw/EwPN4sXL9bvfvc7TZkyRQMGDNCKFSuUmJiop5566oiv8Xq9mjRpkubNm6fTTz89itVGDouJAQAID1PDTUNDgzZv3qy8vLzAMafTqby8PG3YsOGIr7v//vuVnp6u3/72t8f8GfX19fJ4PEG3kxEb+AEAEB6mhps9e/bI6/UqIyMj6HhGRobKy8s7fM0HH3ygVatWaeXKlZ36GYWFhUpJSQnccnJyTrjuSKAbOAAA4WH6tFQoqqurdcMNN2jlypVKS0vr1GtmzZqlqqqqwK2srCzCVR6fUnYnBgAgLEJqvxBuaWlpcrlcqqioCDpeUVGhzMzMduf/8MMP2r59u8aOHRs45t9/JyYmRt98843OOOOMoNe43W653e4IVB8+Xp+hsv1MSwEAEA6mjtzExcVp6NChKioqChzz+XwqKirSyJEj251/9tlna+vWrSouLg7c/vEf/1EXX3yxiouLT9opp2PZeeCQGr2G4lxOZaUkmF0OAABdmqkjN5JUUFCg/Px8DRs2TMOHD9eSJUtUU1OjKVOmSJImT56sXr16qbCwUPHx8Tr33HODXt+jRw9Jane8K/EvJu7dM0EuJ5eBAwBwIkwPNxMmTNDu3bs1Z84clZeX67zzztO6desCi4xLS0vldHappUEhoxs4AADh4zBs1hHT4/EoJSVFVVVVSk5ONrscSVLhG1/pz+//qBsvyNV9/3iO2eUAAHDSCeX729pDIl0EG/gBABA+hJuTQAndwAEACBvCjckMwwjqKwUAAE4M4cZkew42qLbBK4dD6p3KZeAAAJwowo3JSvc1j9pkpyTIHeMyuRoAALo+wo3Jtu9hvQ0AAOFEuDEZi4kBAAgvwo3JSlsWE/fpyWJiAADCgXBjsu0te9zkMnIDAEBYEG5M5u8rRTdwAADCg3BjIk9do/bVNEhijxsAAMKFcGOi0pYpqbTuceruNr2HKQAAlkC4MZG/p1QfuoEDABA2hBsTbW+5UiqXKSkAAMKGcGMi/7QUi4kBAAgfwo2JSvb5G2YSbgAACBfCjYn8a264UgoAgPAh3JikrtGrck+dJOk0FhQDABA2hBuT/LS/VoYhdXfHqGe3OLPLAQDAMgg3JmnbDdzhcJhcDQAA1kG4MQndwAEAiAzCjUnoBg4AQGQQbkxCN3AAACKDcGMSuoEDABAZhBsTNHl9+mk/e9wAABAJhBsT7KqqU6PXUFyMU1nJ8WaXAwCApRBuTODfmTgnNUFOJ5eBAwAQToQbE7T2lGJKCgCAcCPcmKC1pxSLiQEACDfCjQlKWva4oacUAADhR7gxQWDkJo1pKQAAwo1wE2WGYQT2uGHkBgCA8CPcRNnug/WqbfDK6ZB6pxJuAAAIN8JNlJW2TEll90hQXAwfPwAA4ca3a5Rt50opAAAiinATZXQDBwAgsgg3UVayj27gAABEEuEmypiWAgAgsgg3Uca0FAAAkUW4iaKqQ43aX9soiZEbAAAihXATRf7LwNO6u9XNHWNyNQAAWBPhJopau4EzagMAQKQQbqKIbuAAAEQe4SaKWruBs5gYAIBIIdxEESM3AABEHuEmigg3AABEHuEmSuoavSr31EmSTjuFaSkAACKFcBMlZS1tF5LcMUpNjDW5GgAArItwEyWBtgtpiXI4HCZXAwCAdRFuooQrpQAAiA7CTZSUtkxL9WExMQAAEUW4iRL/tFQu4QYAgIgi3EQJ3cABAIgOwk0UNHl9+mn/IUlSbhojNwAARBLhJgp2HqhTk89QXIxTGUnxZpcDAIClEW6iwN8NvE/PRDmdXAYOAEAkEW6ioITFxAAARA3hJgpKWEwMAEDUEG6igIaZAABED+EmCvwb+BFuAACIPMJNhBmG0WbkhmkpAAAijXATYbur63Wo0SunQ+rVI8HscgAAsDzCTYT52y70Sk1QXAwfNwAAkca3bYTRDRwAgOgi3EQY3cABAIiukyLcLFu2TLm5uYqPj9eIESO0cePGI567cuVKXXTRRUpNTVVqaqry8vKOer7Z6AYOAEB0mR5u1qxZo4KCAs2dO1dbtmzR4MGDNXr0aFVWVnZ4/rvvvquJEydq/fr12rBhg3JycnTZZZdpx44dUa68c+gGDgBAdDkMwzDMLGDEiBE6//zztXTpUkmSz+dTTk6Opk+frrvvvvuYr/d6vUpNTdXSpUs1efLkY57v8XiUkpKiqqoqJScnn3D9x3Le/X/XgdpG/W3GReqfFfmfBwCAFYXy/W3qyE1DQ4M2b96svLy8wDGn06m8vDxt2LChU+9RW1urxsZG9ezZs8Pn6+vr5fF4gm7RUlXbqAO1jZLYwA8AgGgxNdzs2bNHXq9XGRkZQcczMjJUXl7eqfe46667lJ2dHRSQ2iosLFRKSkrglpOTc8J1d5a/G/ipSW4lxsVE7ecCAGBnpq+5OREPPfSQnn/+eb3yyiuKj4/v8JxZs2apqqoqcCsrK4tafYGdiXsyagMAQLSYOpyQlpYml8ulioqKoOMVFRXKzMw86mv/7d/+TQ899JDefvttDRo06Ijnud1uud3usNQbqsAeN7RdAAAgakwduYmLi9PQoUNVVFQUOObz+VRUVKSRI0ce8XUPP/ywHnjgAa1bt07Dhg2LRqnHhW7gAABEn+kLQQoKCpSfn69hw4Zp+PDhWrJkiWpqajRlyhRJ0uTJk9WrVy8VFhZKkhYuXKg5c+Zo9erVys3NDazN6d69u7p3727a79GRErqBAwAQdaaHmwkTJmj37t2aM2eOysvLdd5552ndunWBRcalpaVyOlsHmJYvX66Ghgb95je/CXqfuXPn6r777otm6cfEtBQAANFn+j430RatfW7qGr06e/Y6SdJns3+t1G5xEftZAABYXZfZ58bK/D2lkuNj1CMx1uRqAACwD8JNhGzf0zol5XA4TK4GAAD7INxECN3AAQAwB+EmQkroBg4AgCkINxGy3X+lFN3AAQCIKsJNhDAtBQCAOQg3EdDo9WnH/kOSpFz2uAEAIKoINxGw88AhNfkMuWOcSk8yp68VAAB2RbiJAP9i4j49E+V0chk4AADRRLiJgNaeUkxJAQAQbYSbCCgJbODHYmIAAKKNcBMBdAMHAMA8hJsIoBs4AADmIdyEmWEYgT1uTuvJyA0AANFGuAmzyup61TX65HI61Cs1wexyAACwHcJNmPm7gffqkaBYFx8vAADRxrdvmLGYGAAAcxFuwqx0L+EGAAAzEW7CjG7gAACYi3ATZnQDBwDAXISbMPP3laIbOAAA5iDchNGB2gZVHWqU1Nw0EwAARB/hJoz8ozbpSW4lxLlMrgYAAHsi3ISR/zJwpqQAADAP4SaM/N3AWUwMAIB5CDdhVEJPKQAATEe4CaPABn5pTEsBAGAWwk0YtW7gx8gNAABmIdyEyaEGryqr6yXRegEAADMRbsLEvzNxSkKseiTGmVwNAAD2RbgJk301DUpJiGXUBgAAk8WYXYBVjDzjFH0+9zLVNXrNLgUAAFtj5CbM4mPZmRgAADMRbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKXEmF1AtBmGIUnyeDwmVwIAADrL/73t/x4/GtuFm+rqaklSTk6OyZUAAIBQVVdXKyUl5ajnOIzORCAL8fl82rlzp5KSkuRwOML63h6PRzk5OSorK1NycnJY3xut+Jyjg885Ovico4fPOjoi9TkbhqHq6mplZ2fL6Tz6qhrbjdw4nU717t07oj8jOTmZ/+NEAZ9zdPA5Rwefc/TwWUdHJD7nY43Y+LGgGAAAWArhBgAAWArhJozcbrfmzp0rt9ttdimWxuccHXzO0cHnHD181tFxMnzOtltQDAAArI2RGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEmzBZtmyZcnNzFR8frxEjRmjjxo1ml2Q5hYWFOv/885WUlKT09HSNGzdO33zzjdllWdpDDz0kh8OhO+64w+xSLGnHjh3653/+Z51yyilKSEjQwIED9emnn5pdlqV4vV7Nnj1bffv2VUJCgs444ww98MADnepPhCN7//33NXbsWGVnZ8vhcOjVV18Net4wDM2ZM0dZWVlKSEhQXl6evvvuu6jVR7gJgzVr1qigoEBz587Vli1bNHjwYI0ePVqVlZVml2Yp7733nqZNm6aPP/5Yb731lhobG3XZZZeppqbG7NIsadOmTfrzn/+sQYMGmV2KJe3fv1+jRo1SbGys/va3v+n//u//tGjRIqWmpppdmqUsXLhQy5cv19KlS/XVV19p4cKFevjhh/X444+bXVqXVlNTo8GDB2vZsmUdPv/www/rscce04oVK/TJJ5+oW7duGj16tOrq6qJToIETNnz4cGPatGmBx16v18jOzjYKCwtNrMr6KisrDUnGe++9Z3YpllNdXW3069fPeOutt4xf/epXxowZM8wuyXLuuusu48ILLzS7DMu78sorjalTpwYdu/rqq41JkyaZVJH1SDJeeeWVwGOfz2dkZmYaf/rTnwLHDhw4YLjdbuO5556LSk2M3JyghoYGbd68WXl5eYFjTqdTeXl52rBhg4mVWV9VVZUkqWfPniZXYj3Tpk3TlVdeGfS/a4TX66+/rmHDhunaa69Venq6hgwZopUrV5pdluVccMEFKioq0rfffitJ+vzzz/XBBx9ozJgxJldmXdu2bVN5eXnQvx8pKSkaMWJE1L4Xbdc4M9z27Nkjr9erjIyMoOMZGRn6+uuvTarK+nw+n+644w6NGjVK5557rtnlWMrzzz+vLVu2aNOmTWaXYmk//vijli9froKCAv3xj3/Upk2bdPvttysuLk75+flml2cZd999tzwej84++2y5XC55vV7Nnz9fkyZNMrs0yyovL5ekDr8X/c9FGuEGXdK0adP05Zdf6oMPPjC7FEspKyvTjBkz9NZbbyk+Pt7scizN5/Np2LBhWrBggSRpyJAh+vLLL7VixQrCTRi98MILevbZZ7V69Wqdc845Ki4u1h133KHs7Gw+ZwtjWuoEpaWlyeVyqaKiIuh4RUWFMjMzTarK2m677Tb99a9/1fr169W7d2+zy7GUzZs3q7KyUj//+c8VExOjmJgYvffee3rssccUExMjr9drdomWkZWVpQEDBgQd69+/v0pLS02qyJr+8Ic/6O6779Z1112ngQMH6oYbbtDMmTNVWFhodmmW5f/uM/N7kXBzguLi4jR06FAVFRUFjvl8PhUVFWnkyJEmVmY9hmHotttu0yuvvKJ33nlHffv2Nbsky7n00ku1detWFRcXB27Dhg3TpEmTVFxcLJfLZXaJljFq1Kh2Wxl8++23Ou2000yqyJpqa2vldAZ/1blcLvl8PpMqsr6+ffsqMzMz6HvR4/Hok08+idr3ItNSYVBQUKD8/HwNGzZMw4cP15IlS1RTU6MpU6aYXZqlTJs2TatXr9Zrr72mpKSkwNxtSkqKEhISTK7OGpKSktqtYerWrZtOOeUU1jaF2cyZM3XBBRdowYIFGj9+vDZu3KgnnnhCTzzxhNmlWcrYsWM1f/589enTR+ecc44+++wzLV68WFOnTjW7tC7t4MGD+v777wOPt23bpuLiYvXs2VN9+vTRHXfcoQcffFD9+vVT3759NXv2bGVnZ2vcuHHRKTAq12TZwOOPP2706dPHiIuLM4YPH258/PHHZpdkOZI6vD399NNml2ZpXAoeOf/zP/9jnHvuuYbb7TbOPvts44knnjC7JMvxeDzGjBkzjD59+hjx8fHG6aefbtxzzz1GfX292aV1aevXr+/w3+P8/HzDMJovB589e7aRkZFhuN1u49JLLzW++eabqNXnMAy2aQQAANbBmhsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAtudwOPTqq6+aXQaAMCHcADDVjTfeKIfD0e52+eWXm10agC6KxpkATHf55Zfr6aefDjrmdrtNqgZAV8fIDQDTud1uZWZmBt1SU1MlNU8ZLV++XGPGjFFCQoJOP/10vfTSS0Gv37p1qy655BIlJCTolFNO0c0336yDBw8GnfPUU0/pnHPOkdvtVlZWlm677bag5/fs2aOrrrpKiYmJ6tevn15//fXI/tIAIoZwA+CkN3v2bF1zzTX6/PPPNWnSJF133XX66quvJEk1NTUaPXq0UlNTtWnTJr344ot6++23g8LL8uXLNW3aNN18883aunWrXn/9dZ155plBP2PevHkaP368vvjiC11xxRWaNGmS9u3bF9XfE0CYRK3/OAB0ID8/33C5XEa3bt2CbvPnzzcMwzAkGbfcckvQa0aMGGHceuuthmEYxhNPPGGkpqYaBw8eDDy/du1aw+l0GuXl5YZhGEZ2drZxzz33HLEGSca9994beHzw4EFDkvG3v/0tbL8ngOhhzQ0A01188cVavnx50LGePXsG7o8cOTLouZEjR6q4uFiS9NVXX2nw4MHq1q1b4PlRo0bJ5/Ppm2++kcPh0M6dO3XppZcetYZBgwYF7nfr1k3JycmqrKw83l8JgIkINwBM161bt3bTROGSkJDQqfNiY2ODHjscDvl8vkiUBCDCWHMD4KT38ccft3vcv39/SVL//v31+eefq6amJvD8hx9+KKfTqbPOOktJSUnKzc1VUVFRVGsGYB5GbgCYrr6+XuXl5UHHYmJilJaWJkl68cUXNWzYMF144YV69tlntXHjRq1atUqSNGnSJM2dO1f5+fm67777tHv3bk2fPl033HCDMjIyJEn33XefbrnlFqWnp2vMmDGqrq7Whx9+qOnTp0f3FwUQFYQbAKZbt26dsrKygo6dddZZ+vrrryU1X8n0/PPP6/e//72ysrL03HPPacCAAZKkxMREvfnmm5oxY4bOP/98JSYm6pprrtHixYsD75Wfn6+6ujo98sgjuvPOO5WWlqbf/OY30fsFAUSVwzAMw+wiAOBIHA6HXnnlFY0bN87sUgB0Eay5AQAAlkK4AQAAlsKaGwAnNWbOAYSKkRsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAp/x+QZI5CXLXiVgAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "epochs = 10\n",
        "\n",
        "test_acc = np.zeros(epochs + 1)\n",
        "test_acc[0] = test(model, test_dataset)\n",
        "for epoch in range(epochs):\n",
        "    # Initialize the tqdm progress bar\n",
        "    with tqdm(train_loader, desc=f'Epoch {epoch+1} - Test accuracy: {test_acc[epoch]:.3f}') as pbar:\n",
        "        for data, label in pbar:\n",
        "            data, label = data.to(device), label.to(device)\n",
        "            # convert labels to one-hot encoding\n",
        "            label = F.one_hot(label, num_classes=output_size).float()\n",
        "            trainer.train_on_batch(\n",
        "                inputs=data,\n",
        "                loss_fn=loss_fn,\n",
        "                loss_fn_kwargs={'_target': label}\n",
        "            )\n",
        "    test_acc[epoch + 1] = test(model, test_dataset)\n",
        "    pbar.set_description(f'Epoch {epoch + 1} - Test accuracy: {test_acc[epoch + 1]:.3f}')\n",
        "\n",
        "plt.plot(test_acc)\n",
        "plt.xlabel('Epoch')\n",
        "plt.ylabel('Test Accuracy')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wIlwOcPP8oF6"
      },
      "source": [
        "The trained model achieves a classification accuracy of above 95% on MNIST which is comparable to a backpropagation trained model with the same architecture."
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "gpuType": "T4",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.11.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}