DeepRetinotopy_explain.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "DeepRetinotopy_explain.ipynb",
      "private_outputs": true,
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/felenitaribeiro/explainability_CorticalSurfaceGDL/blob/main/DeepRetinotopy_explain.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oM2bI7F9ygV9"
      },
      "source": [
        "# Overview\n",
        "\n",
        "* Setting up the programming environment \n",
        "  * Integration with google drive\n",
        "  * Models download from Open Science Framework\n",
        "  * Installation of miniconda, PyTorch and PyTorch Geometric\n",
        "  * Setting up modules and helper functions\n",
        "\n",
        "* Explainability for a surface-based geometric deep learning model of retinotopy\n",
        "  * Imports and helper functions\n",
        "  * An occlusion experiment with fixed neighborhood size\n",
        "  * Generating a saliency map\n",
        "  * Visualizing predicted maps with intact and occluded features\n"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Setting up the environment"
      ],
      "metadata": {
        "id": "1Kw65i_b2PtY"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N35YN_yJqZHz",
        "cellView": "form"
      },
      "source": [
        "# @title Setup google drive integration\n",
        "\n",
        "from google.colab import drive\n",
        "import os\n",
        "\n",
        "google_drive_dir = '/content/drive/My Drive/deepRetinotopy_explain'\n",
        "local_scratch_dir = '/content'\n",
        "drive.mount('/content/drive')\n",
        "\n",
        "if not os.path.isdir(google_drive_dir):\n",
        "  os.mkdir(google_drive_dir)\n",
        "os.chdir(google_drive_dir)\n",
        "\n",
        "\n",
        "# # Select where to store the data - a free google account is sufficient to store \n",
        "# # all data in this example. \n",
        "# data_directory = google_drive_dir\n",
        "\n",
        "# However, if you need more storage, you can store the data on the local \n",
        "# scratch diretory:\n",
        "data_directory = local_scratch_dir\n",
        "\n",
        "# Define working directory\n",
        "work_dir = os.path.join(data_directory,'osfstorage')\n",
        "\n",
        "# Define checkpoint directory\n",
        "checkpoint_directory = os.path.join(work_dir,'checkpoints')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BlurtYp3qoVO"
      },
      "source": [
        "!pwd "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EKoLOaDOrJnt",
        "cellView": "form"
      },
      "source": [
        "# @title Setup OSF integration\n",
        "\n",
        "# Install the client code of Open Sience Framework (OSF)\n",
        "!pip install osfclient\n",
        "\n",
        "# Download data from OSF\n",
        "data_exists = False\n",
        "print('testing if ',work_dir, ' exists ...')\n",
        "if os.path.isdir(work_dir):\n",
        "  for element in os.listdir(work_dir):\n",
        "    if ('B1Map') in element:\n",
        "      print('data already exists. To trigger a re-download: Delete ', work_dir)\n",
        "      data_exists = True\n",
        "  \n",
        "if not data_exists:\n",
        "      print('data does not exist yet ...')\n",
        "      print('dowloading data to: ',data_directory)\n",
        "      os.chdir(data_directory)\n",
        "      !osf -p f4dez fetch Models/deepRetinotopy_PA_LH_model5.pt ./deepRetinotopy_PA_LH_model5.pt\n",
        "      !osf -p f4dez fetch SaliencyMaps/meanErrorVSnodes_dorsalEarlyVisualCortex_10neighbor_model5.npz ./meanErrorVSnodes_dorsalEarlyVisualCortex_10neighbor_model5.npz"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!ls ."
      ],
      "metadata": {
        "id": "KERIvmZ7KZ2Q"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Miniconda Installation\n",
        "%%bash\n",
        "\n",
        "MINICONDA_INSTALLER_SCRIPT=Miniconda3-py37_4.10.3-Linux-x86_64.sh\n",
        "MINICONDA_PREFIX=/usr/local\n",
        "wget https://repo.continuum.io/miniconda/$MINICONDA_INSTALLER_SCRIPT\n",
        "chmod +x $MINICONDA_INSTALLER_SCRIPT\n",
        "./$MINICONDA_INSTALLER_SCRIPT -b -f -p $MINICONDA_PREFIX"
      ],
      "metadata": {
        "id": "iDf-rouwnNpk",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Cuda 10.1\n",
        "%cd /usr/local/\n",
        "!ls\n",
        "!rm -rf cuda\n",
        "!ln -s /usr/local/cuda-10.1 /usr/local/cuda\n",
        "!stat cuda"
      ],
      "metadata": {
        "id": "DeRI8rhPY3Ya",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Updating conda\n",
        "%%bash\n",
        "conda install --channel defaults conda python=3.7 --yes\n",
        "conda update --channel defaults --all --yes"
      ],
      "metadata": {
        "id": "2lgQg_zjlCTC",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Correct path to packages\n",
        "import sys\n",
        "sys.path.insert(0,\"/usr/local/lib/python3.7/site-packages\")\n",
        "sys.path"
      ],
      "metadata": {
        "id": "RDybrZ7xahCc",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title PyTorch installation\n",
        "!yes | conda install pytorch==1.6.0 torchvision==0.7.0 cudatoolkit=10.1 -c pytorch"
      ],
      "metadata": {
        "id": "YJ2o5qQKa7A-",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title PyTorch geometric installation\n",
        "%%bash\n",
        "pip install scipy\n",
        "pip install --no-index torch-sparse -f https://pytorch-geometric.com/whl/torch-1.6.0+cu101.html\n",
        "pip install --no-index torch-scatter -f https://pytorch-geometric.com/whl/torch-1.6.0+cu101.html\n",
        "pip install --no-index torch-cluster -f https://pytorch-geometric.com/whl/torch-1.6.0+cu101.html\n",
        "pip install --no-index torch-spline-conv -f https://pytorch-geometric.com/whl/torch-1.6.0+cu101.html\n",
        "pip install torch-geometric==1.6.3"
      ],
      "metadata": {
        "id": "mWpdU1_wvbfD",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import torch\n",
        "print(torch.version.cuda)\n",
        "import torch_geometric"
      ],
      "metadata": {
        "id": "FiLrXDUprus2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "s_xYtCkk6aqg",
        "cellView": "form"
      },
      "source": [
        "# @title Other requirements\n",
        "\n",
        "%cd $data_directory\n",
        "!git clone https://github.com/Puckett-Lab/explainability_CorticalSurfaceGDL.git\n",
        "%cd ./explainability_CorticalSurfaceGDL/\n",
        "\n",
        "!pip install scikit-learn==0.22.2\n",
        "!pip install scipy==1.1.0\n",
        "!pip install matplotlib==3.2.1\n",
        "!pip install -r ./requirements.txt",
        "\n",        
        "!pip install git+https://github.com/felenitaribeiro/nilearn.git\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Explaining deepRetinotopy"
      ],
      "metadata": {
        "id": "oDdbESlcIY-v"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DKi8zNkm0k3_",
        "cellView": "form"
      },
      "source": [
        "# @title Imports\n",
        "\n",
        "import os.path as osp\n",
        "import scipy.io\n",
        "import torch\n",
        "import torch.nn.functional as F\n",
        "import torch_geometric.transforms as T\n",
        "import sys\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import pandas as pd\n",
        "import os \n",
        "\n",
        "sys.path.append('../..')\n",
        "\n",
        "\n",
        "from torch_geometric.data import InMemoryDataset\n",
        "from torch_geometric.data import DataLoader\n",
        "from torch_geometric.nn import SplineConv\n",
        "from Retinotopy.functions.def_ROIs_DorsalEarlyVisualCortex import roi as roi2\n",
        "from Retinotopy.functions.def_ROIs_WangParcelsPlusFovea import roi\n",
        "from Retinotopy.functions.error_metrics import smallest_angle\n",
        "from Retinotopy.functions.plusFovea import add_fovea\n",
        "\n",
        "from nilearn import plotting\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wk1E7ms8rJ9w",
        "cellView": "form"
      },
      "source": [
        "# @title Helper functions\n",
        "\n",
        "class Net(torch.nn.Module):\n",
        "    def __init__(self):\n",
        "        super(Net, self).__init__()\n",
        "        self.conv1 = SplineConv(2, 8, dim=3, kernel_size=25)\n",
        "        self.bn1 = torch.nn.BatchNorm1d(8)\n",
        "\n",
        "        self.conv2 = SplineConv(8, 16, dim=3, kernel_size=25)\n",
        "        self.bn2 = torch.nn.BatchNorm1d(16)\n",
        "\n",
        "        self.conv3 = SplineConv(16, 32, dim=3, kernel_size=25)\n",
        "        self.bn3 = torch.nn.BatchNorm1d(32)\n",
        "\n",
        "        self.conv4 = SplineConv(32, 32, dim=3, kernel_size=25)\n",
        "        self.bn4 = torch.nn.BatchNorm1d(32)\n",
        "\n",
        "        self.conv5 = SplineConv(32, 32, dim=3, kernel_size=25)\n",
        "        self.bn5 = torch.nn.BatchNorm1d(32)\n",
        "\n",
        "        self.conv6 = SplineConv(32, 32, dim=3, kernel_size=25)\n",
        "        self.bn6 = torch.nn.BatchNorm1d(32)\n",
        "\n",
        "        self.conv7 = SplineConv(32, 32, dim=3, kernel_size=25)\n",
        "        self.bn7 = torch.nn.BatchNorm1d(32)\n",
        "\n",
        "        self.conv8 = SplineConv(32, 32, dim=3, kernel_size=25)\n",
        "        self.bn8 = torch.nn.BatchNorm1d(32)\n",
        "\n",
        "        self.conv9 = SplineConv(32, 32, dim=3, kernel_size=25)\n",
        "        self.bn9 = torch.nn.BatchNorm1d(32)\n",
        "\n",
        "        self.conv10 = SplineConv(32, 16, dim=3, kernel_size=25)\n",
        "        self.bn10 = torch.nn.BatchNorm1d(16)\n",
        "\n",
        "        self.conv11 = SplineConv(16, 8, dim=3, kernel_size=25)\n",
        "        self.bn11 = torch.nn.BatchNorm1d(8)\n",
        "\n",
        "        self.conv12 = SplineConv(8, 1, dim=3, kernel_size=25)\n",
        "\n",
        "    def forward(self, data):\n",
        "        x, edge_index, pseudo = data.x, data.edge_index, data.edge_attr\n",
        "        x = F.elu(self.conv1(x, edge_index, pseudo))\n",
        "        x = self.bn1(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv2(x, edge_index, pseudo))\n",
        "        x = self.bn2(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv3(x, edge_index, pseudo))\n",
        "        x = self.bn3(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv4(x, edge_index, pseudo))\n",
        "        x = self.bn4(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv5(x, edge_index, pseudo))\n",
        "        x = self.bn5(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv6(x, edge_index, pseudo))\n",
        "        x = self.bn6(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv7(x, edge_index, pseudo))\n",
        "        x = self.bn7(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv8(x, edge_index, pseudo))\n",
        "        x = self.bn8(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv9(x, edge_index, pseudo))\n",
        "        x = self.bn9(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv10(x, edge_index, pseudo))\n",
        "        x = self.bn10(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv11(x, edge_index, pseudo))\n",
        "        x = self.bn11(x)\n",
        "        x = F.dropout(x, p=.10, training=self.training)\n",
        "\n",
        "        x = F.elu(self.conv12(x, edge_index, pseudo)).view(-1)\n",
        "        return x"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%cd /content/explainability_CorticalSurfaceGDL\n",
        "sys.path.append('..')"
      ],
      "metadata": {
        "id": "BAQQH7HSUuw1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Data instances are PyTorch Geometric data objects. If want to know more about setting up your dataset, check out the [documentation](https://pytorch-geometric.readthedocs.io/en/latest/notes/create_dataset.html).\n",
        "\n",
        "Here we load our test set:"
      ],
      "metadata": {
        "id": "l8vvaiMdvR2P"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from Retinotopy.dataset.HCP_3sets_ROI import Retinotopy\n",
        "\n",
        "path = osp.join('/content/explainability_CorticalSurfaceGDL', \n",
        "                'Retinotopy', 'data')\n",
        "norm_value = 70.4237\n",
        "pre_transform = T.Compose([T.FaceToEdge()])\n",
        "\n",
        "test_dataset = Retinotopy(path, 'Test', transform=T.Cartesian(max_value=norm_value),\n",
        "                           pre_transform=pre_transform, n_examples=181,\n",
        "                           prediction='polarAngle', myelination=True,\n",
        "                           hemisphere='Left') # Change to Right for the RH\n",
        "test_loader = DataLoader(test_dataset, batch_size=1, shuffle=True)"
      ],
      "metadata": {
        "id": "zCMaOZndLHom"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Running our 'gray patch' on PyTorch Geometric data"
      ],
      "metadata": {
        "id": "EfuzUdU2NYOy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "To exemplify the usage of our approach, we can first find out about our functions by using the help() method, as follows:"
      ],
      "metadata": {
        "id": "JlG8xPV4wd4K"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from Explainability.neighborhood import node_neighbourhood\n",
        "\n",
        "help(node_neighbourhood)"
      ],
      "metadata": {
        "id": "xIbeY3x8VZ-8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "We need to specify the vertex's index and neighborhood size defining the patch we want to occlude.\n",
        "![occludedNeighborhood.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAJcAAACXCAYAAAAYn8l5AAAACXBIWXMAAC4jAAAuIwF4pT92AAAAGXRFWHRTb2Z0d2FyZQB3d3cuaW5rc2NhcGUub3Jnm+48GgAAIABJREFUeJzlvXeUHPd1JvpVdVV3dc7d09MzmIAMDBEIESAIEiIVGUwzycrv+ay93H3Pz3vOOq336L0Nlr3HT7uWLdnPOiItL23ZEmkFM0kUZYkiCQaRIhgADAlgkGYwoWc6x+qu6lDvj577Q3WcnghQ+50zh+BMd1V11e37u797v/tdDv8LIxaL2cvl8rZ8Pj+cTqeHc7ncULFYDBWLRV8kEvHmcjlvpVIxcxwHTdOc5XKZ5ziuZjAYMoqioFwuFyVJSjgcjoTFYom7XK6Iy+WaCgaDk263e9LpdE74/f7c1f6cVwvc1b6AjUI6nXZPT0/flEgkDqZSqT3FYnFPPp8fqVarDfdA0zQUCgWUSiUAgKqqKBQKqNVqqNVqqFarqFarqNVq0DSNvY/jWm+lyWTSvF7vJb/ff3JoaOhkOBx+fXh4+FW3251e5497TeCX1rhisZh9Zmbmw2fPnr39zJkzt9RqtZ02m40zm80QBIEZR6FQQLFYBAA4HA6YzWZ2jEWPBQDsv/S+SqWCSqUCVVVhMBjA8zxKpRKKxSJqtRpMJhM4joMsy9A0DYqioFKpwGw21xwOx+mtW7ce27Vr14/D4fBzgUAgv/F3aP3xS2VcsVis/9SpU59477337p2YmDhSq9WMNpsNXq8XFouFGUs8Hkc8Hke5XEalUoEoiujr64PH44EgCA3HrFQqKBQKUBQFBoMBoihCEASIosiMSu+12nkwTdOY16vVapBlGYqiQJZlVKtVdXh4+OVdu3Y9sWfPnu/5/f7Iut+oDcL73riSyaTz7Nmznz1+/PinU6nUzYIg8OQ1NE1DuVyG0Whk3kNVVSSTSciyzI5hNBoRCoVgsVhgMplgMBgA1JfESCSCZDLJPBfHceA4DjzPQxAEGI1GSJIEu90Oh8PRYmzVahWpVArpdBqapsFut8PpdMJsNoPjOKiqinw+j2QyiUqlUhsYGDh28ODBx44cOfIox3HZjb2ba4v3rXGdP3/+ppdeeunBM2fOfNLpdFpcLheq1SrS6TQSiQSKxSIqlQp4noemaUvGR0DdyKxWKwRBQLVaRS6XQ7lc7voegtlsxqZNm2CxWFoMDKgvp6VSCel0GrlcPcZ3OBzM0HieR7lcRjKZRCKRAMdxhb17937nyJEjD2/duvW11dyrq4X3lXFpmsa/9tprd73wwgtfmJmZudHn88Hr9YLneSSTScRisQaP1A1LGcty38NxHARBgMPhQDgchiiKHV9fq9VQLBYxOzuLQqEAq9UKk8kEk8kEj8cDURRRKpVQKBSQTqfh8/nevvXWW79y8ODBb3McV1n2hV8lvC+MS9M00/PPP/9vjx079ruFQmHI4XDAZrPBYDAgl8thfn5+SaMyGAwwGAyo1WpIpdJQVbVhqeN5HmazxAJ6Wv40TWPLX6VSYbvFThAEAQMDA3A6nSwm6/CZoKoqstksUqkUZFlGrVYDz/NwOp3w+/2wWq0AgFwuh2g0CqPROPmxj33sz44ePfoNjuOUldzLjcQ1bVyaponHjh37zPPPP/9fVVUdcbvdMBgM7KGn02ksLCxAVVVmDHqIooh8vgBVVWEymaCqZUxNXcapU+NQlPqzqb+Fg9EoYvPmUQSDAQwODiAcDsNut6NarcJkMsFoNKJWqyGbzbJz1mq1ttfN8zxEUYTVaoXH44HT6Wy4tmKxiGQyCUVR4HA4YLfbUalUkEgkkMlkUC6XwXEcbDYb/H4/e38ul8PCwgKsVuv0Bz/4wS/feuutX7+WjeyaNa4333zzvmeeeeZ/lMvlzRaLhRkV/aiqikqlglKpxB6Mqqool8solRRmWK+//gtkMhnUahrbtTWCAz13juPgdDpx3333YP/+fTAY2sdO5XIZqqqiVCohkUh09ZpWqxWBQAAOhwOiKLJjyLKMZDKJbLYes7tcLvj9fvA8j0QiwT4PUI/n+vr64HK5AACKoiAej8NsNp+7/fbb/2D//v1PrsEtX3Ncc8Z18eLFPc8888xX5ubmbiOjogC5+YfnedRqNczNRfAv//ITzM7OQZaLKJfLqFarDUH8UiAjMhgM2LJlM8Lhflx//X6MjAx3jJ00TUMqlUI2m0WxWISqqm2XTI7jEAgEEA6HW5bJcrmMXC4HjuNgsVggSRKAegoklUohFouhVCpB0zRYrVaEQiG2K1UUBZlMBqFQ6Kcf+chHfmfTpk3jPX/gDcA1Y1yappmfeOKJP3z99de/YDAYxEwmA0EQ4Pf7W7xWNpvF6dNnEI3GkEqlMDcXQSKRXPU16I3IaDTiV37lDhw9eguLs4xGIwRBaDEQTdNQrVahKAoKhQJyuRxkWUa5XIYgCDCbzXC5XCxYX851kJGlUinmMR0OB3w+H6xWK0RRpBxaZffu3V87evToFziOK6z6ZqwBrgnjGh8fv+3RRx99WFGULaVSCalUCl6vt8WwFEXB22+fwM9+9jyy2fUq2V1ZJsfGduPw4UPguPpSJggCfD4fgsEgy4Wxd+kMgvJrsiwz49K/vpedavNrKJ2SSCRYzCcIAvr6+uDz+djxrVbrucOHDz84MjLy4krvwFrhqhqXpmmmb3/723/005/+9A8A8LIsw2QyIRQKwWq1MqOan1/A5OQkzpyZwKVLl1Cttg+k1xK0THm9HoRCfdi2bQtMJhOcTicGBgbY8tX8nuUcfylomsa8JJWRLBYLVFVFNBpFOp1GtVqFxWKBx+OB2+2GIAgwGAza1q1b/+rmm2/+D1cz4L9qxjU1NbXra1/72rcmJib20YMMBAKwWq3MW5XLZbzxxnEcO/YySiVlWTHUWqEeL/lx+PAh+HxeGAwGmM1mWCwWWK1WmM1mGI1GlqwlD7eU8ej/LssycrkcqtUqnE4nNE1DsVhkns/tdqNUKmF+fp55QZ/PB0EQoKoq5ubmIMsyjEYj+vr64HA4YDQa4fF43r355ps/GwwGT673fWr7Ga/GSZ9//vnPPfLIIw9ns1mLJElwOp3weDzgeR6FQgGaBsTjcRw//ibOnTu/uMO7ek7WYrHg4MED2Lx5tKWOKAgCJEmC1WpFpVJBsViE1WqF3W5nMVEnQ6PfFwoFxONxZiyFQqEhDyeKIjNc+p0gCAgEAnC5XCyJvLCwwIzb5XKhv78fNputdODAgd/avXv3I+t8m1o/30aeTNM04Zvf/OafPPXUU38IAH19fcxTVSoVLCws4PLlaZw8OY5sNoNKpXnndXUMzGDgMTo6iqGhQXg8HtjttiWz9TzPw2QywWazweVyQZ9OaX4tUI/pOI5DrVZDJpNBLBZb/KJpDa9rfi8leMPhMABgZmaGsTxMJhOCwSC8Xi927Njx8NGjR/8dx3HqmtyUHrBhTyubzXq/+tWvPnH8+PGb6WZQ3EKlkHPnzuPdd0+jUulW4bh6HsxoNOK663bjuut2t7AnmqE3Bp7nYbFY4HQ64XQ6YTKZGnac7QyH6qTRaJQZSzfwPI9gMAhBELCwsMBqonTdQ0ND2LFjx4s333zz/U6nc/Vb6x6wIU/q8uXLm7/61a/+8Pz589s5jkN/fz8cDgcrDs/NzWF8/F1MTl7umPVuxNUxMI7jMDw8hEOHPgCbzdbT65v/32g0wuFwwOPxsCJ3u9cCdW9WKBRYDbKXOK5doR4Ai9N27tx54cMf/vCdfr9/YskPsEp0//qtAc6dO3f4i1/84tOxWMzLcRz8fj/sdjtyuRymp2dw+vQZzM/PI5e79vlyVHLKZLI9GZemaS0pCkVREIvFkMlkMDQ0BIfD0fa95XIZCwsLyOVyrFTVfLx256tWqx094cLCAvL5/OZoNPry9PT03YODg68v+SFWgfZV1TXCyZMnP/jHf/zHP47H414KMo1GIxYWFjA5OYU33jiO8+cvrMCwNn7XSJDlIkolBU6nE5IkMc8jiiKMRmPP6YhqtcrSCwBaPI0gCLDZbMwTUR6r+XW0c9Wj2666UCjg7Nmz/oceeui5iYmJj/R0sSvEuq0vJ0+evOOhhx76frVaNXMcB0mSEIvFkM/nUSjIeOedE4jF4qs8y8YvjwYDjxtu+ADuu+8emEwmljIwm80wmUwoFArIZDKQZZnFjgaDgZWFBEGAxWJhBWsiDbJP1GSc5XIZ2WwWBoMB2WwWyWSywTuJooj+/n7IsoxEItGwJC616bDZbMrnP//5T+7du/epNb1JdI71OOh7771352OPPfa4oihGAIxpKcsyisUijh9/C8lkao3OtnEGxvM8BgcHccMNB3Dw4A2w2+tLY/NyRQX1bDaLSqUCm80GWZahqiqcTidsNhvj2LdDNx4Y7SSJn1+tVlk8ZTQaGY2HNgEdGkcgCAJKpRICgYDyqU996t4dO3Y8u9r70/I51vqAZ8+evfmxxx57tlgsWjVNQyaTYfkXWZbx1lvvrEkdsBEbY2A8z+PIkcO466474HQ6l3w95aYMBkPDv1dS/tGDdoL5fB7z8/MoFovgOA4OhwP9/f3geR4zMzPIZrPM8Gn1oEYSr9cLs9mMSqWC/v7+4tGjR+/cunXrCz3fjB6wpjHXzMzM4ccee+xHxWLRWq1WMTs7i/n5eUbxffPN9TCsjUOtVsP8/Dwymd6o7RzHsViJOPe9xmTd4iZRFCGKItxuN0ZGRuD3++Hz+VAoFHDp0iUoioLh4WH4/f6GDiaXy4XR0VGEw2Gk02mYzWZ4PB5Uq1Xzz3/+86cikcjBni6uR6yZcS0sLIw++uijTxSLRZuiKJibm0Mul2OMy3feOYlkcr0Ma+MC/ESingnv9PDL5TJSqRSSySTjzKdSKcbNArobDu0ok8kkUqlUV9YrAEiShHA4jIGBAQwNDQEApqamkM1mEQ6HEQqFGDUpGo0im83C5/Ohv78f8/PzjM6kqqr9+eef/2E0Gt2ygtvSFmuSishms96/+7u/+1E8Hg/E43Ekk0mWryqXyzhx4hQWFqJrcaou0LARy2M+n8fMzCz27LkOJpOp4W/FYhHRaBSZTAbVahWSJDFqtNlshtvthsVigSzLrGREydhqtYpSqYR8Ps8oO7VaDS6XC8FgsGVHqAftWF0uFyRJQiKRYH/r6+uDJEmYmZmBqqq4fPkySqUS+vr6YLPZGnarsiz7jh079nQmkzmyFonWVRuXpmnGhx566ImZmZlt2WyW7ViA+jJy5swEIpH51Z6m16vBehtYuVxBJFLPyzWnHnieR7VaZTGRnqGaz+dZcZm6kux2OzweD+sFyGazjFlLSCaTbCPRiY+vhyRJCIVCDbRvt9sNk8mEmZkZ5PN5LCwsoFQqYWhoqIFftkh+3PHKK698X9O0j3EcV+50nl5gWPol3bFly5a/fPnllx/gOA6RSKThm3Dp0iTOnbuwwWyG9fdeHAeMjAzD6/WyeiDRlguFQsfyFXVqU3BfLBaRzWaRz+eRyWSgKErbCkWtVoMkST3n0dq1tomiCKfTyc5Lu1mLxQKj0djw2lwuN1woFOyPPPLIj3u+Ke2uYzVvfuWVVz7/9NNP/xYF7/obk0ymdIyGjcT6GjIZUzKZYktZLBbDzMxMz3VA/ZeNmBTd7lOpVMLU1BTi8fiS8Vrz8fWgzqSBgQEYDAYUi0VcuHCB8fj1eOONN/79G2+88b8v+WG6YMXGNT8/f93DDz/8sCiKyOVyLTcnm82iXF66xa5d187qsX4G5vf7cNttt2LHjm3IZrOYmZlBJBJBPp9f1hdpud6c8ledAvxyucyaOkg2oB2oBDc6OgpRFFEul3Hp0qWGcIbwT//0T1+bm5vbsawL1WFFxqVpmunrX//6PxSLRXO5XGYfhOd5RqE5d+5C2w/I8zwkyQSfz4fNm0dZE8TaY30MrFgssdaxeDyOdDq9BIujM5ZrYFTkb+cdBUFAsVjExYsXMTc3h0gk0rBDbYbdbkdfXx+Auvckz0jgOA7FYtH6rW9969uaphk7HacbVhRzjY2N/dlPfvKTe2u1GruxdrsdgUAAFosFL774UsOFAnWj6u8PYdOmARw4cD0OHvwAdu/ehWAwiEgksgaloHZY+/irVqvB7XZh8+bNUFUFiqJAkiS2XC77CpfhuRd3dKjVarDb7S20HeLY5/N55PN5FItFmM3mjk0hoigik8mwOFBVVUbaJGmBZDIZqlarpscee+yny/1syzau8+fP3/aXf/mXX6tWqxx98+x2O/r7+yGKIl5++VWcPn2m4T2CIGBwcBB7945hZGQEIyMjsFgs0DQNMzMzuHDh4ro2XKwlSB5p69YtsNvtMBqNjBFKjbOULSfNCbpP9a5uM6s1Go1GuFwu2Gw21nxL71nqGigQb96t5vN55rGo9cxgMLAvgB5knPl8njEqSC9DFEXIsoxSqYTJycnD//iP//jCX/3VX00t514tKxWhaZr5C1/4wsOqqrKrpPYvnucRiczj1KnxBnfvcrmwZcsoAoE61YZeW61WEY1Gce7cBczPLyznMpaJtU9PkNhJKLQbDodD33kDj8eDUqnEZAFyuRxSqRQqlQpTuOE4DplMBpIkwWazQRRFuFwudt9mZmaY+Fw7lEolXL58mbFM9cwMm82GfP4Ky0RVVUxPT6NUKiEYDDZ4MYq/NE3D3NwcarUaEokEk5IKBoPk2fjvfe97f6Np2l6O4zpfWBOWZVxPPvnkH50+fZplcOniKJlYXyK1RfqtAWazBWNjuxAMBhAMBtm3NpvNsqYEnl+PgL4Za2dggUAAN9xwgJVW9EtTPZ6UGrwEdQzR8kkPl9gQZJhms5ndk15iMer2djgc7P6TBIAgCA1xIHG5iKjZzMLwer1Ip9MsF5fNZuF2u2G1WuFwOJDJZDA9Pb3tBz/4wf8N4D/1eq96Nq6FhYU9v//7v/87+g8uCAITywDAuk9MJiMcDiccDhtIfM1utwOof+tjsRiLT6xWK4xGccVB8UbDYKgLhYhiffnqhe5sNBpZLoliG9phOxwO1nVNSdReYzcyML/fz4yW8mHN91PTNORyOVQqlZYapyAIcDqdjLOfSCTgcrlYpzgVwJ977rk/jEQij4ZCofd6ulc9fQoA4XD4H86ePbtV/zur1couAqjf+HC4Hw6HHTablYmi0VJI5RH9B6+74mRP+aHVY/XeS5YLuHhxEtPTMzCZjMxTtMuet+NmZTIZpgVB1JhsNsuSqL1wsQhEg67VarBareB5nt3ndvoVZMz6pKmeF0bGV6lU4HK5mMFS0F+pVAzJZHLHM888881e7lVPqYhTp07d/8orr3y04Y08D5/Pxy6OdjLUfyeKIgKBADOsSqWCWCzW0DhQ/8DakgHs2mH16YlaTUM+n8fFi5dw+vQZzM7OdoyPmpc3ykUlk0mWj6IMffM96DVNUavVWBBP9BoytGZUKhUmPNd8HpPJhJGREUiSxNQQAbRQhMbHxz88Pj5+dy/XtqRxaZomfve73/3vzUZBRVICuWiSM/J4PA07qFQq1eKdKpUKJieneqawrA3WJv9VLqtIJlPI5/MNHqflbLrfS5IEr9fLhESWvNIeDaxUKiESiTDyoCAIHZfrXC7Xkv8i9qokSUxJh7wV0dMJtVoNP/zhD/9M07QlQ6olX/DKK6/85smTJzdbLBYmtrHYzcteU61WEY/H2bfPbDbD6XQydWOq9jffrPaSRu8PVKs1ZDL1KgQ1q3ZTw6Hg3+PxwGq1MkoNKdh0wlJNGQDYl7dUKsHhcDCKdbv35nI5TE5OssC+VCrB6XQyY/T7/SgWi4yZ4XA44Ha7EY1G2fO9ePHitp///Oe/DuBvu11X16+Ppmmmxx9//AtOpxNer5dZvNvtZuuxoiiIRqMoFK4Iq1ArfjQaZZ0u7Za+el1N7vkbunZYm/MViyWWq1rKE+k/IzWrDgwMgATten1v9+spIhaLsdCk03tzuRympqYwPT2NmZkZTE9Ps2criiKTespkMgDqHldP+dE0Dc8+++x/1jTN1HJwHbrekVdfffXfzM7ODgaDQaRSKWiaBpPJxNwkeSxSKiaUy2XE4/GuBVlN0xCJzK8hl365WL2B5fN5RKOxls/YySPr7xHP83A4HBgYGEB/fz8TXul4tW16Eduh03mb30uaq5VKBclkEpOTkyx2NJlMcDgcSCQSiMViAOoORf8Z5ufnN7322mu/0e1aOhqXpmnC008//fs+nw+lUontPjweDzjuioSiPmFHoK12t5thtVoXJZLWvXWyC1ZnYOVyGXNzc5ifv9LhXKlUkMlkEI/H28Zizf8viiJ8Ph8GBgba6uC3XPEqvHw3A6XlUpZlcByHcDgMq9XKurftdnuLh33xxRf/Q7fYq6M/fuCBBz793HPP/YbNZsP8/DwTMvP7/SAac7FYXPaHNRgMMJlM8Hq9kGV50SWvipO2Blh5iiKVSmN6egbZbG6xeCwzbVNFUViRu3nZbE5k6ssupIxIoQf1Lep1NQRBWHU40S5VIssy7HY7605KJBKw2+2wWCwNDboAkE6nXYODg+9+/etff7fd8Tsa1/Dw8DfK5XI4k8kwvg/dhGZOeK+g9AVpxkciEczNda/ebwxWblyk/Dc1dRnz8wuslqcf10JTOtoVkPUP2GAwwGKxsEoGtfwD9WXJ5/MxlWjS4tKr4awUzW1xFosFFosFoiiiWCyiWq0y8WE992sx/TTy7LPP/k2747Y1roWFhZv++Z//+T9bLBbWvUOg/ruVgFT2VFVFKpXCe++dYWv61cfqEqzVahXJZAqxWBylkgqjUYTZLLHOJ7PZ3NIAy87cxouR/hc9aBqGYDKZ2N8ogauqKmsbW02lQ59QdTqdTEBlfn6e5c/000QAIJlM9j/++OPP/sVf/MVs8/HaGteBAwe+GI/H91FJQo9u35LmWlsziGIryzJkuYhIJLKObIiVYPUZfFUtszRDrabBbJZYQVkftJfLZeTzedbTSLOF6O9EVeZ5nonLEeh1xJ8zm81wOBxM3IWkxjuVgpa8CxwHt9vNMv6FQgGFQgE+n6+BogNcUT/8/ve//3TzcVqCsWQy6fzKV77ySUEQ2Fa0l4uhb5UgCIwF0A7kBWlIE8cBV0EwsAPWpsBdZ3zEkM8XkEymsG/fnpZcGElP0uSMQqHAxrU0c9rbgQR97XY7q9sCYCRGTdMY7YckLntdPklynYwrHA5jfn6eyWY2l5beeeedT2ua9ntc06yiFuO6ePHiZwuFgqWu594bu8JsNiMUCkEQhLaU53ZIp1NYWIheQ4a19pBlGbFYnAnANcdXlMqhIL2uo1H3EHp5pXag/KIsy/B6vcwgqftH0zRWoKbCeSKR6MmLVSoVKIoCk8kETdMgiiLrf7TZbC1E0HQ6bXvzzTc/CeAb+t+3XP3x48c/Q0bSTuOpOW6wWCxwu93sw6XT6bbGRfQSUm0pFkvXKBNiba2djKZSqbTEVrQk0n2g2uPMzMySnsZms8HpdCKTyWBubq5BhdBgMLCUhqZpjK0SCAR6Log3OxYqzrcz+lqthjfeeOPTzcdpeFWhUAhduHDhiKZpLbEW8X4omQaAFafJLesnrDZfmMfjgd/vZ1wnu90Gs7lVEfnaQHcD4zhA6JFPksvl8LOfvYCf/vRnmJq6zJonSqVS2wetaVcGjOpb0ZpB976vrw+VSgWRSKQhTdB8TGpb66WmSc267c5rNBpbeh0B4OzZs7fm8/mg/rUNy+KpU6c+WalUeBpyqYckSSxrS5l60qjSa3k2ey1ypV6vl7W6A3WPZ7fbr7GAXo/28ZfNzOPeO4/ALPH4/lPHkMx1N0RN0xCLxXHs2MuYmZnFoUM3wGKx4Ny587BaLXA6HW0Tp6qqIpPJoFAosN1icxKThlAZjcYWtkMzaKQLPdduHoyO2+41pPXabMj5fN5w4sSJBwB8jX7X8KlOnDhxj9FobNF0EEURwWCQ7XBoN6G/AH0WX/8+amHneb6hIl+raS30m2sZHAfs3TWIj3/0VlhcAwB4/Pb/uRk//tEP8eZ7C2jRBm6Coig4e3YCyWSSbe8dDge2bBnFyMgITKYrQby2qGBIE2VJGok6pxuv68pYl25eiQyYGjk68cZoIFanJZk2b+02e++99969aGdcmqbZf+d3fueIKIoNRWgATOFudnaWTV6l39HNaNdTZ7fb2ZZWURTGaNS0+ihg8mLXLurey+WQ8Gv3fRzhoS0Ab2aLJmcO4Vfu+zT27zuJ7z75ApJZreuCWqvVGjQzisUSI1jqjQuo7ziJoqSn1JBWf6dmi06gdAiRAunY7YyIErR0HVTiI9koWq2a3zs+Pn5U0zQbx3F5QGdcExMTH1ZV1VitVhuWROJl07eIkml6V14ul1sMUhAE1v6kaXWC3RVPxYHnDeB5A4BrMaivg+OAo4f34aYjh2Gyetv8nUOVsyO0+TD+r387hGef+QHeOhNf0ovpUe/aVmCz2Vjuipo7isUim+1DVGhKqi6378BkMqG/v5/lyDpRfej5U9qkWq2yJCo5FP2YZz2KxaLp3LlzHwTwQ0BnXOPj47cDdZJYc/eOJElsqXS5XHC73Q0frlAoNCxx5FqJTEjSi1dcMeBw1IcAdApCrzYGwkHcc/fH4A8OAtwSlBjw4Cxh3P3AZ3H9xRP4pyeOIV3QekqzpFIZnDr1LjZvHsXIyBArjwH1Z2GxWBiD1GazweFwsAe7HAOjpCrtHkkXnyamORwOFvKQc6EEKfH+FUVhlQa3292gpgPUPfPp06dvR7NxTUxMHBUEoWGXKAgCvF4v290YDAa4XK6GdZ+22vptsNfrZd8uitOaS0ZmswSXy4l0Og1AQ612bSS8RFHEXXd8CLt27YJRWlqxmVD3Yg6EttyE3/4/hvCjZ36At88ksdSYIlVVMTV1mQXkoVCIKUVLksR2bSR1qQ/q9Vlyis1Ii6sTR4ziZTIuknzyeDxwOBwsUaqHIAjI5/OYnZ3FyMgIe8bNxqVpGiYmJm5h7wOAdDrt/uIXv7idrBSor+FEZKMBRvVOnSuxgaqqbFjoCxM1AAAgAElEQVQ5XYTZbIbdbocgCIww2K6jpX6BdTZrLpdbp47r5WHfvutwy02H4A2EVnwMDQYYbJtw7yc+iz1n38H3n34FGbn7F6dWqyEeT+AXvziOUqmEQ4cOYmhoE9vyU9DeDrSrpCVTkiQUi0W43e62jbAAGoro/f39NIyK9U7q4zeaXAuAFbENBkODIesxPz9/XTKZdHo8ngwZ102VSoXXB+RGoxFutxuyLC96F7Qk0Ei5mFy02+1mxVT6ezdel9/vg9frwdTUZSQSyatGefb5vLjrzo9j06ZNMAhLz0PsBRXeheGdR/HvBobxgx88g1PnlvZi2WwWr7/+BuLxBA4dugG7d+9aUu+eOqoURUEul2NfalmW4ff72TLaCaIooq+vj8kKNL9WVVWWBVic6wigcWY4cMWLZjIZPplMHgbwrAAA0Wj0oCRJbPem91rE4aEiKYGCdDoo6aZTnEVdKd2yzPTtqX/LTJDljWgvuwKDwYDbbj2K/dfvhdVqX/oNy0QNPAT7MB745Oew5/Q7+N7TL0NWunuxcrmMiYlzSKVSiMfjOHDgAAIBf9vdINUrSTKA2sKAKx6H6DLEnGi3XNLvqOBNSzEVrPXno7jLYDDA6XS2yDppmoapqalDAJ41AMAdd9zx7+fm5naoqopiscjKBZQyIA0Bfes4lSr0za36+lmpVGqQr+wGm80KRSmvo2ZqIzgO2L9nJx74xD3YsXMXjMauVPBVQ+MleAOD+MDeYUQjU0hllK4pC3qws7NzyGQyMJnq8al+0IGiKExvtdPOjzyYqqpIJpMol8stcRu9rh6aXBlmpSgKG/bVHNtRQwflRJvP7XA4kk8++eR3BACIxWJ7ALBRblQ/1O9Imq2+VCoxF9muV65QKPTUj0hlpZ07tyMSiSyZaV4teB74V5+7B4MjOwF+bZbAXlCDAUbHED77+f8NJ996DU/9y3EsJV8myzLefvsEkskkMpkbMTY2BrNZYkqEuVwO5XK56+qgKApzAhQzBQIBmEwmNsUjnU4jm82yVFMmk2Gph2aQrKYoimxpbH7Oc3NzewDAoGma/ZlnnvlSuVzmSA2FdnqLGk1QFAVer5cteZqmsW8MADaeV594SyQSPWXgRVGEx+NhRdhodH3Jgz6XGbd96CPgRfO6nqcjeAl94WFU8jOYmk0v+fJ621gas7NzzDhIIqnXZmIyPjIwWv5SqRRTeNanHyjX1Uxjp38TJZsMq7mPQpZl1xtvvPHnQiaT2a4oCovimqm4FosFxWKxoQGW1l6CyWRqeF+lUunJsDiOY3wm2hC0S86tJQ7esAcG09rHV8sBDw1KcXnzjhKJJF555VVMTV3G6OgI+vqCEHqtnuugb1Cmru9OaPccNE1DNBplvalULdC/VlVVPp1ObxEymcxwt4dpNpsbhC6AulGIosg8V/OSSFtV+oZQ5plExYhr5HA42HJbKBQQicyv646R57FOKobLgwBlRVPaisUSLly4CFmWYbWaGxgqy0Gd8rTyzVO5XGa10U75tGg0OizEYrHhbr2F5Fn0IF4PGZl+Fwlc2TkSr4jyKNRNUi6XIUkSa6WqJ1uvVOvXy3NZJQEul29djr0cGDgVxdLK+hCurBpXb6gpALZDJTp28zOLxWIjQjQaHWpXStAPQmqXJ7HZbMyA2m2TaSoXnZwugDRUjUYjM6x6E2kVPp8X09Mzq/pWdcO+67ZBkFxLv3CdwdUUqKsoqZbLZUQi84vNGlcndiRjMplMsFqtLRuxVCo1LGQymb5m4yKxMNJtaodOOpuEdhwlfVylT9iRjjsFm+uFXbt2ANy6jpjsDbXSkjvFbshksjh1ahz5fAFbtozC7Xa1dQDrCdLDNRqNCAaDTMrpyjVmgryiKL5sNsv+QA+6VCp1Hdu2UpCna+6VU9UyVLX7tno1MAqA23v1l0QA0Koq1MrqPmc+X8DExAROnhzH/PzCBspQ1UGlPwAtLJn69eV9fC6X85EOBACWTKMq+kZAFEU4HHaEQn0rDlKXwrYtAzBZPEu/cANQKa/OcxEURcWlS5N4++0TuHhxcsMZJrTKtGspzGQyPkFVVTcJZxgMBpRKJdZD1xzIrxdoDo7f78PAQD/y+XzPnUe9Yv/ePeBWmDQtGAz4dn8YNQCfmZuDo7o6y1CVAirVtfHQ9TF9C0yqamRkGDabdam3rQn0PZYkmUVQVdXDV6tVs74gSe6tm4CYHrStJcnDlYLUcwYGwgiH+5eUFVoOBAMQXBT0XwkuWKyYNUmImCScW4MHJ+cSqKxhxqVOi85gfPw9jI+/17Yksx7QU6Wbey0rlYpZIPYpGRf912QyNbg68m4UxOmDcVIU9Pl8HafNLwUqN/T390OSzCgWi5ibi6zoWM0I+pyw2le+JG6RCxgqFVEDh21tVH2WAwNXQyoRx3qk82RZxsTEOdRqVRw4sH/dVx4iNVC9Up+SqFQqRqFcLhsp5Q9c6XPTUzX0hL9yudygZ05lAj1bspf2pXYgTlFdnjqLTCbbQp9eCW644Xpwwsq37JZqFf9qZnrV1wEAIqcgkcyu23SieudWlQ2/stvtEMX1kaki1gUAlpai1atWq5l4juNEAKyy7nA4WA6LQIp1qVSKFTcJ+gI2pRPagbyeqqotMwX1oHG927ZtxcjIMHh+9bvVwYHwqo+xVuChIJ4urmunealUwpkzE3jjjTdx7tz5lmfSbmTeSkGOpI1DEYVarVYGYMzlcoy92Jz3UhSlISbTd/3o6R5EeSaRXWoAtVgsTECfJqc6nc4GIdd2qJeKVvfhrRYTHO5rIwUBALymIFNY33hobi6CaDQKVS0vatsr2LZtC5NvdzqdLE5eTWxGJT3gCrdMh7LAcZwKwKiqKlRVhSiKDZRaEtkl0NJIOgLN21+9lnm1WgXP8wgEAigUCozR2suE+kwmsyaB6Q0H9kC8yoVqPbSqguwStOfVQh9DZ7M5vPfeaShKCQcOXI8tW7bAZrOxdrVexWbagXoryJnoYTAYFMFoNLIil6IoLXVCmuqgPyBdOC1zepTL5YYsO7VJ6Y2QdBI6QdM0JJNJ5POrj7c2bx5Z9THWEpWygnxpY5tRFEXB+fMXYTJJCAaDjLPXqcpCu7+lhOWq1SpT5yFevo4GrfImk4ktyM1rMzW76o2FZBj1/98JgiBAkiTWd6e/KD2HqBnlchmxWHzVSUGjaID3GsnKE0rFIopXoZtOVVWMj4/jxz/+F5w5c7bFCeghiiL6+/s71pX1oHRV82okCEJRsNvtCQBhAC1iFVRo1oOq8mSpNputrYVT+5EkSS2zm6kPrhMbo1QqdVTLWQ7sdsuabAjWEnm5tOrSz0pRLlfw3ntnUCopOHz4UMeityiKbNxfc0+qHsSKARpXNACwWq1x3u12s+azZgPQ7wT1r6Flkpop9S1QtNujank74wLqbMZOS6MkmREO97MldaVIJHP4+7/7e2TjlwBtY2tvnRCNJVC+ipdSrVZx4cJFvPDCiy1yCpQK6uvrYx6p2xdcb1z6OA8A7HZ7grfZbHH9iwmapjU0YOihz8aSZrnT6WQyScFgkBEM9Tk0/Xmak7T689ZqVQwNbUIwGOj4wXpFJJbHXz/0j5gYfxWotg5b2lhomJy++oJ35GXoy019kaFQCIODg6zRpl7z7ZwUJ0cCoIXLb7Va44LD4VigXjfa5RmNRsadb3fAZlaq2+1mxDEyGL2h2u12loag1qVOjItyuYwTJ04iFouvaiejh1IG/unxF3DT1DRu/fBHYTD71+S4ywVXUzATyWA9BowuF8ViCbOzcxgdHUEgEGBd8vovfDuiqB6iKHY0LrfbHRWcTueUJEnsj+RlSFi1GXprJSxVg6StbzKZZJrmeu9HpSVKtCqKiomJc2sqIV6tAS8dv4CLlyP41Cfvg9M3BG1lI75XBEM1jfjCNObj14YeWTKZwokTpyCKRmzfvr1tp4+++77dCqaf21jXt71S/nG73Zd4t9s9SXktPXWC8ljtDrjcojJVzUljnaaZEUiYLBaLQZZl7Ny5HUNDQ+tCgJuNyvj633wbF8ZfBldZfaqjGzhUIZQXEJ98Bf/499/A1x95AkWFHtLV18ZIpzN4++0TePPNt1sGTVAi3ePxdPReiqI00G708Pl8k4Lb7b6kT5rSizrJ5ABg3UDdjKwuDVRiXHtqqm2XQCW9AxrnVtc4WN0okm6QFQ3ffvwYbtx/GR/6yO3gpbVdJg0oQSst4OLESTz/8knMJWodgvirvzwmEgmk02m2q6aNXD6fZ5PpOhER9ZWc5lE0fX19lwSLxTLhdru1qakpDrhiXG3S+ewg0WgUbre7a/mGunwpWKTG2XYXSGkJou/E44l1FyapVIFXjk9ieubv8YkH7obDvxna8kZ+t0CspaHk5vDG6y/hrTNxRNO1JfUhrraBUZhTLBaZYmA8Hmc6tjQWpl2IQn0QABp6HB0OR81qtZ4XOI7LPfTQQ5cEQRhtHg/cDmTZmUymQXREj0KhwIrY6XS64SKaUavVWJ5MVcu4ePHS4hCE9d/ZaQCm5ot46G+/i/vvugmjO28AhOWVigyogC8vILVwEa+8+jrOTpeQLmjrQqlZD9RqNRw//ubikPsA0uk0KpUKgsEgeJ6HwWCA3+9HoVBoMTDa8Te3qvX395/nOK4gAMDAwMBJg8EwWi6XmZqdvmOnHYrFIvL5fIv30rT6oG4S383lcpAkqUUwjkATHwDyijFWg9wo5Esavv3EKzh06RI+9JE7YbAsLaEkaAVUSzFcOPMWXj1+BtPxGopLiIx0xtX1XvVlMAeD4cpOX69eaLVa2UDPZnugzZi+RDg4OHgKWNTnCoVCJ81m87369IPFYmFyPO2gaRoymUzLqDQq9egbPtLpdMtASIJeOZoo1lcDlSrw87fnMD3zTXzi/jvgCu5AjWueZKFBqCSRz0Rw/LWXcepcAslcp3hqubh6BkZTygiiKDb0T5DqUTqdbiExtBPuHRoaumJcAwMDr9ntdhYjAVc0Tbspz7TjBVE6QQ8Sw2g3uYt2qJqmLY4cuXoTzGoacDmq4uG/fwp3f2QGO/behJrBCQNXRVmOYm76Al489hpmYyUUr/agtTWEz+dh7GKr1YpAINDSnEPt+7FYjBkU5cWIWgXU7SYcDv8cWDQuURR/HgqFatPT0zwZC0lCd0O7nV9zGYBAlBu9DBNdtNvtXpxcKqO6dAS87sgXNXzvmTdxeHYO+z5wE955+y2ceG8S+aLWQ4C+Gmy896rXgH0wmUzweDzw+XxtE9yUTkokEg3sU6DuwSge83g8VZvN9jqwOEHD7XanBwcHT9O4EEIvFfHm11QqlY6jcJvpO3SOuoZEBbJ8NeZdt0e5Crz0dgR//Y3v49gbl5AprLdhETb281utFoyMDGHTpk2L/QvtpS4BNIQ7ejFAvcjfwMDASY/HkwF041m2bt16zGq1Nmwpu3F+gKUz881o199GzIv64Mu1bSdbLTQNG2RQLWfesDMNDAxg586dPQ2B1zsNvW3oHcbWrVtfon+zow0MDPzY6/U2MCFIPa4T2iVRSXC3HSqVSosKHsktlkpFqOq1KRv+ywqDwYDR0ZGOYr7N0Hu1dqrSgiBg586dP6bfM+Pq7+//aTAYVPQvbj6IHp3YpDR8sh1IFE6fw6Jtb6VSXVedCP356hKbTrb1vjax/t7LarVg8+bRnsts+rBF/18Sf3O73aWBgYEX6PXs7nIcV9i5c+crRqOxISDvxkbs5Ea7NceSxgBtaWnjYDAIG6LmLIoihoY2YdeunRga2rRhXeUrw/oamN/vR19ffdAY7daX4m8B9eRpKBQCx3EN6aOxsbEXOY5jnqPBOsbGxp5wOp0NxkXzlNudqBMfa6nOa31gCFyZL7MRO8X63CIrhoY24ZZbbsaBA9d3TPBeG1gfA+M4Dlu3bmFTTObm5rCwsNDVuKg3keSxADRMqtu/f//j+tc3WEcgEPjupk2bqvoEKAm1tbu4dmhHDmwH/YdQFHVJWfG1giSZYDabYTQa4fP5cODAfhw6dAObhPq/CsxmCVu2bEYul8Pk5CSi0ShEUex6D4xGIyRJYlUZGoMIAD6fr7pp06bOxmWz2eb37dv3EsVMdZGx9u1H1GTRbBDEy+oGVVURj8fZLqPdOL31AC3BNNScirJjY7tx8803Ydu2revWnbw6rO2Xjud5bNq0CT6fF3Nzc8jn80w6q1uFhLwWrWS0QQOAsbGxn9nt9qj+9S1munPnzsfozTTKrlOg3U5RmGpNS6HOfoijrn0vr3saQhAENvGLBIL11YiBgfDiMIGNIxAuD2tjYIIg4KabbsR9993T0nfaTlBED47j4PP52EaOcpqCIOCGG254rPn1LcY1PDz8WCAQKAAra/sWBKHtPMBm0C4jHo8jnc6sKeu0GQaDAYODA7juut0Ih+uzbqhpl1BXoK6sqbrOtQifz4vbbvsg+vtDTZPkuLbTaJuhz3tSvNXf35/fvn37d5tf22JcHMdl9u/f/x2gfsO7naydAfE8D5fL1bP8Eg0AX480BGWRPR43HI66/pfdbkOlUkEqlWpoXzObzdiz5zpcf/3+nvr1rg5W7722bt0Cn8+PYrHYEIos1YzRDjTh7tChQ49yHNfC324bve3du/dvyELbPXSKXTo1WWia1nNwXp/iUFyXNITJZML+/XvxoQ/dtihqcuXj1mq1Bl1WjuPgcjmxdetmbNu2BV6v5xoN8FduYGazGWNjuyEIhpYeCepT7BU0ydZiseADH/jAN9q9pu3dczgcP7/++utfp6C35U08D6/XywhleixX4IKabNdypyiKAvx+PwYGwvD7/RgertfO9LJQRHrUl7tIGn14eAi7du3EwEAYknQt5sFWdq+8Xg8CAT8URWkY2kqzLpdDd4rFYqhWqzhw4MAbPp/vF+1e03HtOnDgwF+89dZbj0mS1DJ+w2azdWxqrVarSKVSPQvAVirVNY+3zGYLDhzYj+HhIVaykCQJJpOJpTyIyGgwGBgViOd5JtxPLIGFhSjOnp2AXjf2/YpSqYS33noHo6PDrHXPZDIhFotBVRU2ZodmKXYCkRCMRiNuueWWL3V6XUfjGhwc/N727dsvZjKZUZpkRVBVta1h0UNbzkNox/9aLajDSM+SpWYCGudGNJFmHX2e59lU1r4+HsPD9dHA4+PjmJmZYcMYrj6WT8+JxxN48cWXcPHiRSQSSaZgVE8r5XDp0hR27tyBTZsGYbfbWbc8DasgUOpo165dl0ZHR5/odL6OxsVxXPXcuXN//u677/5/NH1Uf/B2FGdZlpedDDWbJWzaNLjIm1/9cAOqVVqtjcVYalWnwQulUgm5XK7tt1S/JRcEATt2bIPT6cD4+LuYmDi3ITm53rB8A8tmszhzpoBardrQ+X3mzAQEQUA8nsDOndtx881HFsV85+F0Ohs2OdFoFLVaDXfeeeeXOI7ruER1jVi3bNnyN1u2bJlsx5PXC8IRarVa28QqoXm4AVD3FMFgYLEE0+1qekNdW9XatqpAM4koXRIIBGC321tkDPTMSnpfONyPw4cP4eDBDyAQCFyju8neUB+t0vg7SohevHgRoljX21pYWEAqlUKhUGCft1gsIplMYufOnVNbt259pNt5uuYLOI5Tp6am/vTSpUsPJZPJhoCPJpEGAgH2LaclplO85XA4YDabUa1WkUwmG+Y11omHPDRtdbtGg4FnMeESn63t0q6qKlKpFMtG01JPinx79+5hceVGsDiWxtqyV3negNnZWRiN9aQzDQEjxGIxGAwG3H777X+0KBzY+VhLnWzTpk2PjI2NnWse01KtVpHL5ZiLBOqxjsfjaflWU0BNdSmi09LDFUURW7ZsRjjcv2rJI0EQ2050WArksagDimYRRiIRNsG+UqmwXebVmsfdHmu30agPvY/BaDSxOJWYI5VKBdlsFvv27Tuza9euf1jqWEs+AY7jyvF4/A/OnDnzBE0TFUWR7fByuRwSiQR8Ph9r/6Zvtl5rgIQuaFl0u93geZ5NlOU4Drt374QgGHD58syyx42QJzKbzbDZrMvOUZEmWKlUYk26+l0yTemy2x3Yvn07isXSmutZrA5r48GMRiOuu24X47oJgsB6KRYWFsBxHO65557f5ThuyZ1NT0/A5/M9eeONN/44FArB4XAwYV5COp1mN5nnefh8vhZhi2YBEyo39PX1wev1MuWb+jK7/EkXZrMZBw5cj8OHD8HrbfWeS4GW9OYeAI7jYLfb4Xa7F6m9AkKh+jW/n+OuTpAkEyTpSmOO3+9nI4uTySSOHj36TCgU+lEvx+r5633kyJHfD4fD5VAoBKfTiYGBAcbpIV0IgsFgQDAYbKucogd5G/oBaMSabdkPzmQyob8/hM2bRyEIwrL7H0nLwu/3s9wYXZ/FYmHGBdSXjmw2g0rlWoi59Fj98mizWWE0ikzYLxiskwmj0ShcLpd6xx13/F6vx+rZuOx2+/jHP/7xL1PWXhRFBINBpoZDCjW0nNHfKWfUDSaTiX1D3G4XrrtuF0KhvmUtbSZTnadVLpeRTCYRjUaX3U1EIwADgQATtNMv54R6GiPPMtvXFlZ+TSaTEXv37oHD4UA4HMbwcL1kRhSpe++9908dDseZXo+3rKj3uuuu+68nTpy45913390JXNHOJFpOJBKBJEkIBAIQRRFGoxGhUIgJw3UCBY3kcRwOB8bGdkGSJEQi84xO2y2INpslSJKJ9U3m83lUKhW43W6mcEgJVOBK/VPfjELeymq1srppO2aIJEnYvHkUxWIR8/Pz10Sv5WrB8zxuueVmHDp0EE6ns4GSlEqlMDY2dvbgwYP/73KOuSzj4jhOSaVSvzU9Pf1cNpvlm8fQUta7Wq2iv7+/JbPbDc3FbrPZjB07tqO/v84QjcVimJyc6sj7qsdDRpTLV4R8S6USotEoMxp9KoTUjMvlcgOTo1qtwmw2d71uSZJw3XVj8Pl8OHnyFCYmzl01GYL2WF5wz3EcDh8+hLvuuqNllVm8R7XPfOYz/5rjuGWR7pZNu3S73S+88847X/nOd77zu7TLa1566sPQI6ybuleOVPNxRFFg7ITh4SEMDQ3hzTffgiwXYTSKMJlMLMPudDohCAaUSo27TP2uM5VKMQ5SszGrqgqO41Cr1diS2G05J4JhNpvF5OTUNWZcwHIMbOvWLfiVX7mrbW5wYWEBd95553/3eDwvL/cKVsTp3bt37xcuXrz4oR/84Af7aJfVvGTl83nk83lGy9CrPrcDJVL15Rg9W0KSJOzZcx36+/uhqiokyaSTZuLA83XvRCmDduhGwdanFJLJJEqlEtxuN6zWzmkN8tbNE7veTwiH+/H5z38WVqulRVRElmVs27btrSNHjvyXlRx7RcbFcZyiKMrnT548+Yt33nmna0dlPp9HoVBAoVBAKBTqaGCiKCIQCDAd/PoMwTSKxSLbTdZHvXRWAWwnbb4SaJrWkH/rtOvVNA0+nxejo8M4e3YChcK1UnMkdPdemzYN4jOf+SQcjvquX29ciwnU/P333/+ZpTLxnbBiNpzJZHr3wQcf/NcDAwNLfmPpYc3OznadAEvdJUajke0gBwYG0N/f37aLm5Y2+umlOWS56LYRoRpkfTm9Vr1W63VxHIfNm0fw67/+eZZq0IcJi//WPvGJT/ymw+GYWOmZV0W1DIfDj/7Gb/zGX3aixzb3NtJovV69CyVe9bQY2gnS4KlEIsF+crncmngug8HAxgJ2W8qp7ubxeBhN5dpMrDYaWDjcj4985EOQJInJKzQZFg4fPvzloaGh76zmrKu+E5qmic8+++xPvvWtb32wmW5D5SBKC1SrVZhMJgwODnZ9aN1QLpexsLDQUJrRC7+uNu6hOmggEIDZbO5qLCSiks1mIctFZLNZnD59BqlUCn6/H9FoPfdnt9uRy2WhKFezVMQtzhAQcOONB7F58+hi48ogbDYbcwQ8z2P37t3P3Xbbbbf3UuLphlU36XEcV9Y07f5UKvXqc889tz0evyKUS4xFn88HnueRSCTYjqxXNGsT0FLUTj5xGdfMfkj3kxLDRqOR/XcpL8Rx9SkTJINut9vA8/UlMhzuRzweh6Ko8Hg8mJycxOnTZxankmz8Emo2S9ixYzs8Hjf8/jrV2Wyuj3umbi1N0zA6Onr6tttu+7XVGhawBsYFABzHJYvF4h2pVOq11157LUBdIQBNIIstduF4uorvNoP0VYvFIhspUq2ujhZNVGZJktiSq98wNPPNOl0XjQNUFIVNGymXy4sGZ0A+n0cwGGTynz6fF319QZw8OY7JySnd/CS0cKvWGk6nAwcOXA+/39cwNExVVZZ+0DQNoVAo/uEPf/hXOY5LdTter1iz9mKz2XxJluVfnZqa+kkul7PrPUmlUkEmk0E+n2e8KH2zBHBF1VnfUUSZ+1wux7qCVwvKylO5qBdDIsMul8tQFIX90Izndl3nhUKBjaOhzprR0VH4fD6cO3cBp0+fhs1mhckk4fz585BlGRzHs3TKWqDeWT2IsbHdbEKZPoSgJmFN0+D1erN33333HZIknV+Tk2MdNBJnZ2ePfOlLX/rx1NQUUy/heR5msxmKoqBSqTAKjsfjaQjUk8kkW2b0D524Y4lEYs0oLvU6pps1ZLQDLev5fJ5l87vl0QhkUFQRaPbU9b7JNKxWCzKZDMbH38XcXAROpxO1Wg1zcxEUi0W4XK6W/sJeIAgCBEHA9u1bMTo60lDaouvjeR4jIyNwOp1wu93yvffee6ff739xWSdaAuuytYlGo7f/yZ/8yRORSMREM4XoQS4sLLAlwW63NzBZKZbq9LAVRUEymUQ2m10TLyYIApuc2qn/UpZlJJNJNDepdANpt9P0r06gZT+VSiObzcBgMKBWqyGRSKJW07B9+zak02lcvHgJMzOzUFUFg4ODKJUUJJOJxb4GDqFQH6rVutaD1+tGKHzDEt0AAAmsSURBVBSCJNWnyekHfenjzEAggL6+PrjdbuXXfu3XftXr9f7Lim5iF6zbvjkSidz25S9/+am5uTlboVCfsUNj8yg5WiqVIEnS4s1o5by3Q7FYRCQSWXa5pd1uUhRF9PX1LUkNymaziEajPdOaqUWrm1ckUOxD6RWS3+Z5nslXiaKITCaLSqWMUCgEVS0zpqzRaMKWLaOsk0cfRtDnbvZafr8f/f398Hq98n333XdfMBhcc8MC1lk6OBqN3vLwww//YHx83EHdQ2azmQ2LrOtEpOFyuRAMBnvKEfViXDzPw263NywH9EPDGagBtq+vb0nV6kKhwBQRe/VeRqMRfr+/5xZ5am6JRqMs5qI2N1JYpjY82oC004CneC+RSDBaNr2GGMDhcBihUCh///333+Pz+X7W0wWuAOue8ZNl+eDf/u3fPv38888H9PQWWo5yuRwqlcqSIhgkXCLLMnK5XEcvQio2VquV5W/0oKHfQP3h0a6xh8/RsKQvBUEQWHdRr7y0bDaLVCoFVVWZ7jtxyWRZZk0hFosFfr+/4/2qi6rUewESiQT7rIFAAH6/H5s3b164//7777JarW/2dGErxLqLUVksll8Ui8UbBUH44bPPPrsTqH/4SCSCUqnEuPfdvBYF9LIsd21do2PrH0ozqKtnuWie4bwUKpUKYrEY46d10tXQgxpRK5UKCxNkWUY6nUahUGClJpvNtuTgetJPo/cEg0F4vV6MjY29d9ddd91lNpsne/4wK8SGKJ2ZzeZLmqYd8Xq9//z973//VlVVIYoikskkZFlGX19fx5iLyhGFQgH6/FknEGtSEIQlY6leQZn45aYIyuUyEokECoUCXC7XkjEYNePyPI98Ps9iU/10Cq/X21WFh8bhxONxZlhE0z569OhzH/rQhz7BcdyGDFfaMBk9juNSmqZ9dGho6E/++q//+g8BwOPxIJVKYXp6Gl6vt0GbVC8U0k2ArhlWq7Wn7Loe+tpaJ53X5jG7y0GpVEKpVFrSY5K3yWQybYfO11WoO+vFEzlSHxvWm46DeOCBBx7eu3fvb3Mct2HE/w3VaFwsKfzH2dnZd772ta99Y2ZmxhoOh1lBW1EU9PX1sWWyWq0iHo+39RjE5adyEM0qMhqNPRWQKYCmzvFqtcoGKNFxKXah+GU1dUtZlpFIJNjS1+76ZFlGPB5vq51BxfROn0uWZczOzja8VxRF3HjjjfnPfe5zvxkKhVZVhF4JrloJX1GUHY8++ui3fvKTn1xPPPe5uTlIkgSfz8d474lEoqE7m+B0OuHz+QBcqTdSTbAXVKtVRKPRlvF7tMPUs1X1RmUwGCBJElRVZQlh2hAQmxVoP6+SGkA8Hk/LJoJyXmTwzS1udF5K/jbPCM/n84hGoyzJbDabcffddx9/4IEHPidJ0oppM6vBVeWHaJomHD9+/P955JFH/pMsy7zT6Vws9irwer1M6KRYLLaoD1osFoRCoWUJljWdm+2maBlZyjPRAyZSI9XmRFFkgrWiKLLykCRJLLtutVqZinQ3Zqs+79XuS+V2uxEIBBqOQYZJubi+vj7twQcf/KsDBw78wUqJfmuBa4J8lMvljn7zm998+OWXX95OqYFsNgu73c6aT2VZRiwWY8sT1R0pWF2OkemNqFarMQUfmmzbblmixCglNZfKvANoYICQ0dCwVIopSaGxUqkwuaJcLseYsHqIooiBgYGWzQ+RA4rFIo4ePXr6c5/73L9ZCed9rXFNGBdQ54WdOnXqdx9++OEvRiIRo8vlWpT1vtIOVi+N1ImBegPp1YtR0Zo6xCmpaDQa2d9yuRzbaenBcXVtL+J59fB5WNacxipTkE49B5SJJwHgZnpR8/npy6Y3LlmWMTc3B6/XW37wwQf//MCBA/+F47hrolvkmjEugqIou370ox/9xZNPPvkxVVXR19fXQixUVRXZbBa5XI6VS0gSqVPHDn279e1vAJiiILEf9L2MBGoeoXMsVdKhMcylUqnBaFYLq9WKcDjMzk9Ge/vtt//o7rvv/l2TydRzw+pG4JozLsL8/PzdTz311J8dO3Zsm8vlapuzIj1Voup0imco9xOLxZZV8KYir9VqhdfrZTvUbqhUKlhYWOiak6PSzXKuhZZlKoarqoqxsbEzd9xxx++FQqFnej7QBuKaNS4A0DSNn5+ff+Dxxx//07Nnz24mQYxMJsMeei85LUrCkufqJXgH6jsur9e7ZJNs87kURUGhUGDFY1IFqs8dsjH6UfMEt07geR4ej4cZ+MjIyOWPfvSj/23Tpk3/cy0Yo+uFa9q4CJqmGWdnZ3/zpZde+oMzZ86MUFKSdlM2m60nQ6O8FQm8US6rmwcxGo2M3AhcoTb3AmpNIwqxngWaSqV69qR2ux0DAwPYt2/fhUOHDv2P4eHh/7mRydCV4n1hXARN0/hUKnXXq6+++h9PnTp1UyaTYV1A5XIZoigy+SbSl6C+R30XERkgxWFUyO4EfReTIAhwOBywWCyoVqsst7Yc0RQaC9iNl0aMiIGBARw+fPitQ4cOfTUQCHz7WvZUzXhfGZcesiwfevfddx985513PnX58mVbnXSXgizLsNls8Hq9MBgMLPtP6YZmcZGVBtzkwaibXD82kEiGpVIJgiAsdohLzOgpDUHJ0kKhsEhz5phE+djYWO7AgQOP7dmz5xsWi6Wtzvu1jvetcRE0TbMvLCx86uTJk59+9913b52fnzcoisIeFHBltJ+qqoxh0KtB6dmb1ClE3spms7Ut5VC2XFtU1iH6T7vNQK1WYzSeUChUvf7665/fu3fvo8Fg8Dscx+Vb3vA+wvveuPTQNC04Pz//wNmzZ++dmJi4ZWFhQdIbEnmNbDaLRCLRMlS0HaxWa4OkkMlkWnIAAOlHdNoEUPfQooRBad++fce2bdv2RDgc/h7HcbGV34FrC79UxqWHpmmWVCp12+zs7McnJiY+ePr06bHJyUmeCIJkIPl8HrlcrisJkPJcbrcbLperxQNRbZMESTrx8avVKvx+f3VwcHB8x44dLw4MDPzY7/e/wOlG9/4y4VqcXLkmWHxgP1z8gaZpjkKhcNPCwsKhc+fO7bl8+fKehYWFUVEUefJolKZoBrXt026UJkjIsoxisdhASxZFER6PB36/H16vt9bf338xGAyeCIfDJ4PB4OtWq/XnHMctTUz7JcAvrefqBZqmWQuFwtZSqTScSqVGEonEcCaTCebzeX86nfYWCgWvpmnWxZKNXRRFQRCECsdxuUVmaN5sNicdDkfc7XbH3G531Ol0XgoEApOSJE1ardZzHMd134r+EuP/B4w95R/fKByeAAAAAElFTkSuQmCC)"
      ],
      "metadata": {
        "id": "5Cx11F5xxVIZ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from Explainability.neighborhood import node_neighbourhood\n",
        "\n",
        "data = test_dataset[0]\n",
        "vertex = 1982\n",
        "neighborhood_size = 10\n",
        "\n",
        "new_data, neighbors = node_neighbourhood(data, vertex, neighborhood_size)"
      ],
      "metadata": {
        "id": "RgAAM49mPuWt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "After running our 'node_neighborhood' function to occlude a patch of the input data, you can verify the new (occluded) input features at the specified location (see 'neighbors' dictionary), and compare with the intact features."
      ],
      "metadata": {
        "id": "n6ui7V7rXbq0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "new_data.x[neighbors['level_10']]"
      ],
      "metadata": {
        "id": "1u4Ti_Q8VsZq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "test_dataset = Retinotopy(path, 'Test', transform=T.Cartesian(max_value=norm_value),\n",
        "                           pre_transform=pre_transform, n_examples=181,\n",
        "                           prediction='polarAngle', myelination=True,\n",
        "                           hemisphere='Left')\n",
        "print(test_dataset[0].x[neighbors['level_10']])"
      ],
      "metadata": {
        "id": "n9aGcrtZW_Pp"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Occlusion experiments - `Where` are the most `important vertices`?"
      ],
      "metadata": {
        "id": "03VnJ54ANmjc"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Here is a sample of the code on our github repository to determine saliency maps by occluding target vertices' neighborhoods. You can run it yourself after changing your runtime to GPU. Alternatively, we made available the predicted maps generated with one selected model (model 5, as in our manuscript) after occluding each target vertices within the early visual cortex. So, you can skip this cell entirely and go straight to the plots. "
      ],
      "metadata": {
        "id": "QLtJS0BCYDC2"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import time\n",
        "from Explainability.neighborhood import node_neighbourhood\n",
        "\n",
        "nodes = np.load('/content/explainability_CorticalSurfaceGDL/Models/'+\n",
        "                'nodes_earlyVisualCortex.npz')['list']\n",
        "# neighborhood_sizes = np.arange(1, 21, 1)\n",
        "neighborhood_sizes = [10]\n",
        "\n",
        "for neighborhood_size in neighborhood_sizes:\n",
        "    # for i in range(1):\n",
        "      i = 4\n",
        "      device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
        "      model = Net().to(device)\n",
        "      model.load_state_dict(\n",
        "          torch.load(\n",
        "              '/content/deepRetinotopy_PA_LH_model' + str(i + 1) + '.pt',\n",
        "              map_location=device))\n",
        "\n",
        "      # Create an output folder if it doesn't already exist\n",
        "      directory = '/content/explainability_CorticalSurfaceGDL/Models/testset_results'\n",
        "      if not osp.exists(directory):\n",
        "          os.makedirs(directory)\n",
        "\n",
        "      for node in nodes:\n",
        "          def test():\n",
        "              model.eval()\n",
        "              MeanAbsError = 0\n",
        "              y = []\n",
        "              y_hat = []\n",
        "              for data in test_loader:\n",
        "                  new_data, _ = node_neighbourhood(data, node,\n",
        "                                                    neighborhood_size)\n",
        "                  pred = model(new_data.to(device)).detach()\n",
        "                  y_hat.append(pred)\n",
        "                  y.append(data.to(device).y.view(-1))\n",
        "                  MAE = torch.mean(\n",
        "                      abs(data.to(device).y.view(-1) - pred)).item()\n",
        "                  MeanAbsError += MAE\n",
        "              test_MAE = MeanAbsError / len(test_loader)\n",
        "              output = {'Predicted_values': y_hat, 'Measured_values': y,\n",
        "                        'MAE': test_MAE}\n",
        "              return output\n",
        "\n",
        "          start = time.time()\n",
        "\n",
        "          evaluation = test()\n",
        "\n",
        "          end = time.time()\n",
        "          print(end - start)\n",
        "\n",
        "          torch.save({'Predicted_values': evaluation['Predicted_values'],\n",
        "                      'Measured_values': evaluation['Measured_values']},\n",
        "                      osp.join('/content/explainability_CorticalSurfaceGDL/Models',\n",
        "                              'testset_results',\n",
        "                              'testset-node' + str(\n",
        "                                  node) + '_neighborhood' + str(\n",
        "                                  neighborhood_size) + '_model' + str(\n",
        "                                  i + 1) + '.pt'))"
      ],
      "metadata": {
        "id": "CidjluPILxQM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Saliency map"
      ],
      "metadata": {
        "id": "6ilx-NBvmEjR"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Considering the time required to download all predictions, we also made available the final saliency map, which has been previously generated using ./Manuscript/stats/left_hemi/MeanDeltaTheta_vertices_PA_neighbor.py\n",
        "\n",
        "Following script can be used to plot the saliency map on a spherical surface model."
      ],
      "metadata": {
        "id": "ezH6SAunHr1h"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Background settings\n",
        "background = np.zeros((32492, 1))\n",
        "threshold = 1  \n",
        "\n",
        "# Early visual cortex\n",
        "early_visual_cortex = ['V1d', 'V1v', 'fovea_V1', 'V2d', 'V2v',\n",
        "                       'fovea_V2', 'V3d', 'V3v', 'fovea_V3']\n",
        "V1, V2, V3 = add_fovea(early_visual_cortex)\n",
        "cluster = np.sum(\n",
        "    [np.reshape(V1, (-1, 1)), np.reshape(V2, (-1, 1)),\n",
        "     np.reshape(V3, (-1, 1))], axis=0)\n",
        "cluster[V3 == 3] = 3\n",
        "cluster[V2 == 2] = 2\n",
        "cluster[V1 == 1] = 1\n",
        "\n",
        "label = ['Early visual cortex']\n",
        "\n",
        "# ROI\n",
        "visual_hierarchy = ['ROI']\n",
        "final_mask_L, final_mask_R, index_L_mask, index_R_mask = roi(\n",
        "    visual_hierarchy)\n",
        "\n",
        "# Selecting nodes\n",
        "nodes = cluster[final_mask_L == 1]\n",
        "nodes = np.where(nodes != 0)[0]\n",
        "\n",
        "error_map = np.zeros((32492, 1))\n",
        "\n",
        "# Loading error\n",
        "induced_error = np.load('/content'\n",
        "                        '/meanErrorVSnodes_dorsalEarlyVisualCortex'+\n",
        "                        '_10neighbor_model5.npz')['list']\n",
        "\n",
        "tmp = error_map[final_mask_L == 1]\n",
        "tmp[nodes] = induced_error\n",
        "error_map[final_mask_L == 1] = tmp + 1\n",
        "error_map[final_mask_L != 1] = 0\n",
        "\n",
        "# Location of highest and lowest saliency score\n",
        "print('Most important vertex is: ' +\n",
        "      str(nodes[np.where(induced_error == np.max(induced_error))[0]]))\n",
        "\n",
        "print('Least important vertex is: ' +\n",
        "      str(nodes[np.where(induced_error == np.min(induced_error))[0]]))"
      ],
      "metadata": {
        "id": "lyzlmVoBREIv"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Downloading surface gifti file\n",
        "!osf -p 95w4y fetch data/raw/S1200_7T_Retinotopy181.L.sphere.32k_fs_LR.surf.gii /content/surface/S1200_7T_Retinotopy181.L.sphere.32k_fs_LR.surf.gii"
      ],
      "metadata": {
        "cellView": "form",
        "id": "YVNuMgEqZFbu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Predicted map\n",
        "view = plotting.view_surf(\n",
        "    surf_mesh='/content/surface/S1200_7T_Retinotopy181.L.sphere.32k_fs_LR.surf.gii',\n",
        "    surf_map=np.reshape(error_map[0:32492], (-1)), bg_map=background,\n",
        "    cmap='inferno', black_bg=False, symmetric_cmap=False,\n",
        "    threshold=threshold, vmax=15)\n",
        "view"
      ],
      "metadata": {
        "id": "2ZT3P12HWFit"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Make sure you select the ``Custom_LH`` view in the drop down menu to see the map."
      ],
      "metadata": {
        "id": "cW3slrZFbT6h"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Downloading extra data\n",
        "\n",
        "If you want to check out the saliency maps generated with 4 other models, you just need to uncomment the following cell:"
      ],
      "metadata": {
        "id": "aG0sLSNYXkBV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# !osf -p f4dez fetch SaliencyMaps /content/SaliencyMaps"
      ],
      "metadata": {
        "id": "jXwDUBX4VTsu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fdjWB-mZ2lyJ"
      },
      "source": [
        "## Polar angle maps\n",
        "\n",
        "We can also visualize the predicted polar angle maps when the most important vertex is occluded and compare with predicted maps generated with intact features."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Downloading predicted maps\n",
        "\n",
        "!osf -p f4dez fetch OcclusionExp/testset-node1982_neighborhood10_model5.pt /content/testset_results/testset-node1982_neighborhood10_model5.pt\n",
        "!osf -p f4dez fetch OcclusionExp/testset-node2562_neighborhood10_model5.pt /content/testset_results/testset-node2562_neighborhood10_model5.pt\n",
        "!osf -p f4dez fetch Predictions_IntactFeat/testset-intactData_model5.pt /content/testset_results/testset-intactData_model5.pt"
      ],
      "metadata": {
        "cellView": "form",
        "id": "qQfLrCYbhF34"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from Retinotopy.functions.def_ROIs_WangParcelsPlusFovea import roi\n",
        "\n",
        "participant_index = 0\n",
        "\n",
        "label_primary_visual_areas = ['ROI']\n",
        "final_mask_L, final_mask_R, index_L_mask, index_R_mask = roi(\n",
        "    label_primary_visual_areas)"
      ],
      "metadata": {
        "id": "YOQpbsCMir7G"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Ground truth maps\n",
        "\n",
        "predictions = torch.load('/content/testset_results' +\n",
        "                         '/testset-intactData_model5.pt',\n",
        "                          map_location='cpu')\n",
        "\n",
        "measured = np.zeros((32492, 1))\n",
        "measured[final_mask_L == 1] = np.reshape(np.array(predictions['Measured_values']\n",
        "                                              [participant_index]), # Change the index for different participants\n",
        "                                     (-1, 1))\n",
        "\n",
        "# Rescaling polar angle valuess to match with the right visual field\n",
        "measured = np.array(measured)\n",
        "minus = measured > 180\n",
        "sum = measured < 180\n",
        "measured[minus] = measured[minus] - 180 + threshold\n",
        "measured[sum] = measured[sum] + 180 + threshold\n",
        "measured = np.array(measured)\n",
        "measured[final_mask_L != 1] = 0\n",
        "\n",
        "# measured maps\n",
        "view = plotting.view_surf(\n",
        "    surf_mesh=osp.join(work_dir,\n",
        "                       '/content/surface/S1200_7T_Retinotopy181.L.sphere.32k_fs_LR.surf.gii'),\n",
        "    surf_map=np.reshape(measured[0:32492], (-1)),\n",
        "    cmap='gist_rainbow_r', black_bg=False, symmetric_cmap=False,\n",
        "    threshold=threshold, vmax=361)\n",
        "view"
      ],
      "metadata": {
        "cellView": "form",
        "id": "klrl37bwkaxQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Predictions generated with `intact` data\n",
        "\n",
        "predictions = torch.load('/content/testset_results' +\n",
        "                         '/testset-intactData_model5.pt',\n",
        "                          map_location='cpu')\n",
        "\n",
        "predicted = np.zeros((32492, 1))\n",
        "predicted[final_mask_L == 1] = np.reshape(np.array(predictions['Predicted_values']\n",
        "                                              [participant_index]), # Change the index for different participants\n",
        "                                     (-1, 1))\n",
        "\n",
        "# Rescaling polar angle valuess to match with the right visual field\n",
        "predicted = np.array(predicted)\n",
        "minus = predicted > 180\n",
        "sum = predicted < 180\n",
        "predicted[minus] = predicted[minus] - 180 + threshold\n",
        "predicted[sum] = predicted[sum] + 180 + threshold\n",
        "predicted = np.array(predicted)\n",
        "predicted[final_mask_L != 1] = 0\n",
        "\n",
        "# Predicted maps\n",
        "view = plotting.view_surf(\n",
        "    surf_mesh=osp.join(work_dir,\n",
        "                       '/content/surface/S1200_7T_Retinotopy181.L.sphere.32k_fs_LR.surf.gii'),\n",
        "    surf_map=np.reshape(predicted[0:32492], (-1)),\n",
        "    cmap='gist_rainbow_r', black_bg=False, symmetric_cmap=False,\n",
        "    threshold=threshold, vmax=361)\n",
        "view"
      ],
      "metadata": {
        "cellView": "form",
        "id": "0VW0cm6Ff2Er"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Predictions generated with `occluded` data \n",
        "# @markdown `Most` important node\n",
        "\n",
        "predictions = torch.load('/content/testset_results' +\n",
        "                         '/testset-node1982_neighborhood10_model5.pt',\n",
        "                          map_location='cpu')\n",
        "\n",
        "predicted = np.zeros((32492, 1))\n",
        "predicted[final_mask_L == 1] = np.reshape(np.array(predictions['Predicted_values']\n",
        "                                              [participant_index]), # Change the index for different participants\n",
        "                                     (-1, 1))\n",
        "\n",
        "# Rescaling polar angle valuess to match with the right visual field\n",
        "predicted = np.array(predicted)\n",
        "minus = predicted > 180\n",
        "sum = predicted < 180\n",
        "predicted[minus] = predicted[minus] - 180 + threshold\n",
        "predicted[sum] = predicted[sum] + 180 + threshold\n",
        "predicted = np.array(predicted)\n",
        "predicted[final_mask_L != 1] = 0\n",
        "\n",
        "# Predicted maps\n",
        "view = plotting.view_surf(\n",
        "    surf_mesh=osp.join(work_dir,\n",
        "                       '/content/surface/S1200_7T_Retinotopy181.L.sphere.32k_fs_LR.surf.gii'),\n",
        "    surf_map=np.reshape(predicted[0:32492], (-1)),\n",
        "    cmap='gist_rainbow_r', black_bg=False, symmetric_cmap=False,\n",
        "    threshold=threshold, vmax=361)\n",
        "view"
      ],
      "metadata": {
        "cellView": "form",
        "id": "cro5LdWwi5ZS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# @markdown `Least` important node\n",
        "\n",
        "predictions = torch.load('/content/testset_results' +\n",
        "                         '/testset-node2562_neighborhood10_model5.pt',\n",
        "                          map_location='cpu')\n",
        "\n",
        "predicted = np.zeros((32492, 1))\n",
        "predicted[final_mask_L == 1] = np.reshape(np.array(predictions['Predicted_values']\n",
        "                                              [participant_index]), # Change the index for different participants\n",
        "                                     (-1, 1))\n",
        "\n",
        "# Rescaling polar angle valuess to match with the right visual field\n",
        "predicted = np.array(predicted)\n",
        "minus = predicted > 180\n",
        "sum = predicted < 180\n",
        "predicted[minus] = predicted[minus] - 180 + threshold\n",
        "predicted[sum] = predicted[sum] + 180 + threshold\n",
        "predicted = np.array(predicted)\n",
        "predicted[final_mask_L != 1] = 0\n",
        "\n",
        "# Predicted maps\n",
        "view = plotting.view_surf(\n",
        "    surf_mesh=osp.join(work_dir,\n",
        "                       '/content/surface/S1200_7T_Retinotopy181.L.sphere.32k_fs_LR.surf.gii'),\n",
        "    surf_map=np.reshape(predicted[0:32492], (-1)),\n",
        "    cmap='gist_rainbow_r', black_bg=False, symmetric_cmap=False,\n",
        "    threshold=threshold, vmax=361)\n",
        "view"
      ],
      "metadata": {
        "cellView": "form",
        "id": "4ARJR9ppkPaJ"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}