Lesson7_BP-GPR_MLP.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/cc-ats/mlp_tutorial/blob/main/Lesson7_BP-GPR_MLP.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Lesson 7: Behler-Parrinello Gaussian Process Regression (BP-GPR) for Machine Learning Potentials\n",
        "\n",
        "\n",
        "$\\Delta$MLP with PyTorch for the Claisen Rearrangement reaction\n",
        "\n",
        "For this tutorial, we will be combining the Gaussian Process Regression (GPR) from Lesson 2 and the symmetry functions from the Behler-Parrinello and ANI models from Lesson 4 to train a $\\Delta$ Machine Learning Potential ($\\Delta$MLP) model to reproduce the energy and forces for the Claisen Rearrangement reaction. With a BP-GPR, we can use the BP symmetry functions for feature extraction and the GPR kernel for predicting observables. The goal of this model is to train with data that is from semiempirical (PM3) and DFT (B3LYP) levels of theory. The BP-GPR model will be used to correct the semiempirical values to obtain DFT level accuracy, that makes it a $\\Delta$MLP model."
      ],
      "metadata": {
        "id": "6ZczLhOzG8Wc"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Importing PyTorch Lightning and Libraries\n",
        "\n",
        "We will first install PyTorch Lightning and import libraries needed to train our machine learning models"
      ],
      "metadata": {
        "id": "v208qXsE8h4F"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%capture\n",
        "!pip install pytorch-lightning > /dev/null\n",
        "!pip install gpytorch > /dev/null\n",
        "\n",
        "import math\n",
        "from typing import Sequence, Tuple\n",
        "\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "from torch import Tensor\n",
        "import gpytorch\n",
        "\n",
        "from torch.utils.data import TensorDataset, DataLoader, random_split\n",
        "import pytorch_lightning as pl\n",
        "from pytorch_lightning import loggers as pl_loggers"
      ],
      "metadata": {
        "id": "uuyTEvnmwHBc"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Defining Symmetry Functions\n",
        "\n",
        "Now we will define the radial symmetry functions:\n",
        "\n",
        "pairwise_vector: Distance function for calculating distances between each atom.\n",
        "\n",
        "symmetry_function_g1: BP radial symmetry function.\n",
        "\n",
        "symmetry_function_g2: BP angular symmetry function.\n",
        "\n",
        "symmetry_function_g2ani: ANI angular symmetry function"
      ],
      "metadata": {
        "id": "ng2talt98vmg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def pairwise_vector(coords: Tensor) -> Tensor:\n",
        "    num_batches, num_channels, _ = coords.size()\n",
        "    rij = coords[:, :, None] - coords[:, None]\n",
        "    mask = ~torch.eye(num_channels, dtype=torch.bool, device=coords.device) # remove self-interaction\n",
        "    rij = torch.masked_select(rij, mask.unsqueeze(2)).view(num_batches, num_channels, num_channels - 1, 3)\n",
        "    return rij\n",
        "\n",
        "def symmetry_function_g1(rij: Tensor, Rcr: float, EtaR: Tensor, ShfR: Tensor) -> Tensor:\n",
        "    dij = torch.norm(rij, dim=3)\n",
        "    fij = (torch.cos(dij / Rcr * math.pi) + 1) * 0.5 * (dij <= Rcr)\n",
        "    g1 = torch.sum(torch.exp(-EtaR.unsqueeze(dim=2) * (dij.unsqueeze(dim=-1) - ShfR.unsqueeze(dim=2))**2) * fij.unsqueeze(dim=-1), dim=2)\n",
        "    return g1\n",
        "\n",
        "def symmetry_function_g2(rij: Tensor, Rca: float, Zeta: Tensor, EtaA: Tensor, LamA: Tensor) -> Tensor:\n",
        "    c = torch.combinations(torch.arange(rij.size(2)), r=2)\n",
        "    #print(c)\n",
        "    rij = rij[:, :, c]\n",
        "    r12 = rij[:, :, :, 0]\n",
        "    r13 = rij[:, :, :, 1]\n",
        "    r23 = r12 - r13\n",
        "    d12 = torch.norm(r12, dim=3)\n",
        "    d13 = torch.norm(r13, dim=3)\n",
        "    d23 = torch.norm(r23, dim=3)\n",
        "    f12 = (torch.cos(d12 / Rca * math.pi) + 1) * 0.5\n",
        "    f13 = (torch.cos(d13 / Rca * math.pi) + 1) * 0.5\n",
        "    f23 = (torch.cos(d23 / Rca * math.pi) + 1) * 0.5\n",
        "    cosine = torch.einsum('ijkl,ijkl->ijk', r12, r13) / (d12 * d13)\n",
        "\n",
        "    g2 = torch.sum(2**(1 - Zeta.unsqueeze(dim=2)) * (1 + LamA.unsqueeze(dim=2) * cosine.unsqueeze(dim=-1))**Zeta.unsqueeze(dim=2) * torch.exp(-EtaA.unsqueeze(dim=2) * (d12**2 + d13**2 + d23**2).unsqueeze(dim=-1)) * (f12 * f13 * f23).unsqueeze(dim=-1), dim=2)\n",
        "    return g2\n",
        "\n",
        "def symmetry_function_g2ani(rij: Tensor, Rca: float, Zeta: Tensor, ShfZ: Tensor, EtaA: Tensor, ShfA: Tensor) -> Tensor:\n",
        "    c = torch.combinations(torch.arange(rij.size(2)), r=2)\n",
        "    rij = rij[:, :, c]\n",
        "    r12 = rij[:, :, :, 0]\n",
        "    r13 = rij[:, :, :, 1]\n",
        "\n",
        "    r23 = r12 - r13\n",
        "    d12 = torch.norm(r12, dim=3)\n",
        "    d13 = torch.norm(r13, dim=3)\n",
        "    f12 = (torch.cos(d12 / Rca * math.pi) + 1) * 0.5\n",
        "    f13 = (torch.cos(d13 / Rca * math.pi) + 1) * 0.5\n",
        "    cosine = torch.einsum('ijkl,ijkl->ijk', r12, r13) / (d12 * d13)\n",
        "    cosine = torch.cos(torch.acos(cosine).unsqueeze(dim=-1) - ShfA.unsqueeze(dim=2))\n",
        "    g2 = torch.sum(2**(1 - Zeta.unsqueeze(dim=2)) * (1 + cosine)**Zeta.unsqueeze(dim=2) * torch.exp(-EtaA.unsqueeze(dim=2) * (0.5 * (d12 + d13).unsqueeze(dim=-1) - ShfZ.unsqueeze(dim=2))**2) * (f12 * f13).unsqueeze(dim=-1), dim=2)\n",
        "    return g2"
      ],
      "metadata": {
        "id": "HGM0pzjQxA6J"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Defining Feature Extraction Classes"
      ],
      "metadata": {
        "id": "5BVuAjaQ83gz"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Defining the BP Feature Class\n",
        "\n",
        "Now we will create the feature class, which takes input information and uses functions to extract characteristics unique to the called input features. Here we will be using the previously defined symmetry functions for feature extraction."
      ],
      "metadata": {
        "id": "JdnQPFmMsRo-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "class Feature(nn.Module):\n",
        "    def __init__(self, Rcr: float, EtaR: Tensor, ShfR: Tensor, Rca: float, Zeta: Tensor, EtaA: Tensor, LamA: Tensor) -> None:\n",
        "        super().__init__()\n",
        "        assert len(EtaR) == len(ShfR) == len(LamA)\n",
        "        assert len(Zeta) == len(EtaA)\n",
        "        self.Rcr = Rcr\n",
        "        self.Rca = Rca\n",
        "        self.EtaR = torch.Tensor(EtaR)\n",
        "        self.ShfR = torch.Tensor(ShfR)\n",
        "        self.Zeta = torch.Tensor(Zeta)\n",
        "        self.EtaA = torch.Tensor(EtaA)\n",
        "        self.LamA = torch.Tensor(LamA)\n",
        "\n",
        "    def forward(self, coords: Tensor, atom_types: Tensor) -> Tensor:\n",
        "        num_batches, num_channels, _ = coords.size()\n",
        "        rij = pairwise_vector(coords)\n",
        "        EtaR = self.EtaR.to(device=coords.device)[atom_types]\n",
        "        ShfR = self.ShfR.to(device=coords.device)[atom_types]\n",
        "        Zeta = self.Zeta.to(device=coords.device)[atom_types]\n",
        "        EtaA = self.EtaA.to(device=coords.device)[atom_types]\n",
        "        LamA = self.LamA.to(device=coords.device)[atom_types]\n",
        "        g1 = symmetry_function_g1(rij, self.Rcr, EtaR, ShfR)\n",
        "        g2 = symmetry_function_g2(rij, self.Rca, Zeta, EtaA, LamA)\n",
        "\n",
        "        return torch.concat((g1, g2), dim=2)\n",
        "\n",
        "    @property\n",
        "    def output_length(self) -> int:\n",
        "        return len(self.EtaR[0]) + len(self.EtaA[0])"
      ],
      "metadata": {
        "id": "FPztptTTxIFH"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Defining the ANI Feature Class\n",
        "\n",
        "We will create a similar feature class for ANI that uses the ANI angle symmetry function.\n"
      ],
      "metadata": {
        "id": "jLAoPesH87am"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "class FeatureANI(nn.Module):\n",
        "    def __init__(self, Rcr: float, EtaR: Tensor, ShfR: Tensor, Rca: float, Zeta: Tensor, ShfZ: Tensor, EtaA: Tensor, ShfA: Tensor) -> None:\n",
        "        super().__init__()\n",
        "        assert len(EtaR) == len(ShfR)\n",
        "        assert len(Zeta) == len(ShfZ) == len(EtaA) == len(ShfA)\n",
        "        self.Rcr = Rcr\n",
        "        self.Rca = Rca\n",
        "        self.EtaR = torch.Tensor(EtaR)\n",
        "        self.ShfR = torch.Tensor(ShfR)\n",
        "        self.Zeta = torch.Tensor(Zeta)\n",
        "        self.ShfZ = torch.Tensor(ShfZ)\n",
        "        self.EtaA = torch.Tensor(EtaA)\n",
        "        self.ShfA = torch.Tensor(ShfA)\n",
        "\n",
        "    def forward(self, coords: Tensor, atom_types: Tensor) -> Tensor:\n",
        "        num_batches, num_channels, _ = coords.size()\n",
        "        rij = pairwise_vector(coords)\n",
        "        EtaR = self.EtaR.to(device=coords.device)[atom_types]\n",
        "        ShfR = self.ShfR.to(device=coords.device)[atom_types]\n",
        "        Zeta = self.Zeta.to(device=coords.device)[atom_types]\n",
        "        ShfZ = self.ShfZ.to(device=coords.device)[atom_types]\n",
        "        EtaA = self.EtaA.to(device=coords.device)[atom_types]\n",
        "        ShfA = self.ShfA.to(device=coords.device)[atom_types]\n",
        "        g1 = symmetry_function_g1(rij, self.Rcr, EtaR, ShfR)\n",
        "        g2 = symmetry_function_g2ani(rij, self.Rca, Zeta, ShfZ, EtaA, ShfA)\n",
        "\n",
        "        return torch.concat((g1, g2), dim=2)\n",
        "\n",
        "    @property\n",
        "    def output_length(self) -> int:\n",
        "        return len(self.EtaR[0]) + len(self.EtaA[0])"
      ],
      "metadata": {
        "id": "uLYEFnZexP4g"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Creating the BP-GPR Model\n"
      ],
      "metadata": {
        "id": "gyFHUEgp9CZy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Defining the BP-GPR Class\n",
        "\n",
        "This class will combine the symmetry functions with the GPR model."
      ],
      "metadata": {
        "id": "dL5HsXiEspio"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "class BPGPR(gpytorch.models.ExactGP):\n",
        "    def __init__(self, descriptor: Tensor, Y: Tensor, likelihood, learning_rate=5e-4) -> None:\n",
        "        shape = descriptor.shape\n",
        "        X = descriptor #.reshape(shape[0], shape[1]*shape[2])\n",
        "        super(BPGPR, self).__init__(X, Y, likelihood)\n",
        "        self.mean_module = gpytorch.means.ConstantMean()\n",
        "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
        "\n",
        "    def forward(self, x):\n",
        "        mean_x = self.mean_module(x)\n",
        "        covar_x = self.covar_module(x)\n",
        "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)"
      ],
      "metadata": {
        "id": "lt8r5J9GxVpK"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Importing Data and Seeding\n",
        "\n",
        "We will now import our semiempirical and DFT calculation data."
      ],
      "metadata": {
        "id": "-7Tv-R5NJSuF"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "ds = np.DataSource(None)\n",
        "\n",
        "qm_coord = np.array(np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/qm_coord.npy\", \"rb\")), dtype=\"float32\")\n",
        "atom_types = np.loadtxt(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/qm_elem.txt\", \"r\"), dtype=int)\n",
        "\n",
        "elems = np.unique(atom_types).tolist()\n",
        "atom_types = np.array([[elems.index(i) for i in atom_types]])\n",
        "atom_types = atom_types.repeat(len(qm_coord), axis=0)\n",
        "\n",
        "energy = np.array((np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/energy.npy\", \"rb\")) - np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/energy_sqm.npy\", \"rb\"))) * 27.2114 * 23.061, dtype=\"float32\")\n",
        "energy = energy - energy.mean()\n",
        "qm_gradient = np.array((np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/qm_grad.npy\", \"rb\")) - np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/qm_grad_sqm.npy\", \"rb\"))) * 27.2114 * 23.061 / 0.529177249, dtype=\"float32\")\n",
        "\n",
        "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
        "\n",
        "qm_coord = torch.from_numpy(qm_coord).to(device)\n",
        "atom_types = torch.from_numpy(atom_types).to(device)\n",
        "energy = torch.from_numpy(energy).to(device)\n",
        "qm_gradient = torch.from_numpy(qm_gradient).to(device)"
      ],
      "metadata": {
        "id": "QTw4Jjy8tC6m"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Using ANI for Feature Extraction\n",
        "\n",
        "The ANI model is set to be used and the ANI parameters are defined.\n"
      ],
      "metadata": {
        "id": "s01tP30n9J1q"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "ani = True\n",
        "\n",
        "if ani:\n",
        "  # From TorchANI\n",
        "  Rcr = 5.2000e+00\n",
        "  Rca = 3.5000e+00\n",
        "  EtaR = [1.6000000e+01]\n",
        "  ShfR = [9.0000000e-01,1.1687500e+00,1.4375000e+00,1.7062500e+00,1.9750000e+00,2.2437500e+00,2.5125000e+00,2.7812500e+00,3.0500000e+00,3.3187500e+00,3.5875000e+00,3.8562500e+00,4.1250000e+00,4.3937500e+00,4.6625000e+00,4.9312500e+00]\n",
        "  Zeta = [3.2000000e+01]\n",
        "  ShfZ = [1.9634954e-01,5.8904862e-01,9.8174770e-01,1.3744468e+00,1.7671459e+00,2.1598449e+00,2.5525440e+00,2.9452431e+00]\n",
        "  EtaA = [8.0000000e+00]\n",
        "  ShfA = [9.0000000e-01,1.5500000e+00,2.2000000e+00,2.8500000e+00]\n",
        "  EtaR, ShfR = np.array(np.meshgrid(EtaR, ShfR)).reshape(2, -1)\n",
        "  Zeta, ShfZ, EtaA, ShfA = np.array(np.meshgrid(Zeta, ShfZ, EtaA, ShfA)).reshape(4, -1)\n",
        "  EtaR = np.repeat([EtaR], 3, axis=0)\n",
        "  ShfR = np.repeat([ShfR], 3, axis=0)\n",
        "  Zeta = np.repeat([Zeta], 3, axis=0)\n",
        "  ShfZ = np.repeat([ShfZ], 3, axis=0)\n",
        "  EtaA = np.repeat([EtaA], 3, axis=0)\n",
        "  ShfA = np.repeat([ShfA], 3, axis=0)\n",
        "  descriptor = FeatureANI(Rcr, EtaR, ShfR, Rca, Zeta, ShfZ, EtaA, ShfA)\n",
        "else:\n",
        "  Rcr = 6.0\n",
        "  Rca = 6.0\n",
        "  ShfR = np.zeros((3,6)) # H, C, O\n",
        "  EtaR = [0.0, 0.04, 0.14, 0.32, 0.71, 1.79]\n",
        "  EtaR = np.repeat([EtaR], 3, axis=0) # H, C, O\n",
        "  Zeta = [1, 2, 4, 8, 16, 32, 1, 2, 4, 8, 16, 32]\n",
        "  Zeta = np.repeat([Zeta], 3, axis=0) # H, C, O\n",
        "  EtaA = [0.0, 0.04, 0.14, 0.32, 0.71, 1.79, 0.0, 0.04, 0.14, 0.32, 0.71, 1.79]\n",
        "  EtaA = np.repeat([EtaA], 3, axis=0) # H, C, O\n",
        "  LamA = [1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1]\n",
        "  LamA = np.repeat([LamA], 3, axis=0) # H, C, O\n",
        "  descriptor = Feature(Rcr, EtaR, ShfR, Rca, Zeta, EtaA, LamA)"
      ],
      "metadata": {
        "id": "lAfRrZy1xj1y"
      },
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## GPR Hyperparameters"
      ],
      "metadata": {
        "id": "LUjFyU0UwYmx"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Reshaping and Initializing\n",
        "\n",
        "The atom types and coordinates are reshaped with the energies for use in the regression kernel. Then the marginal likelihood function is initialized."
      ],
      "metadata": {
        "id": "1A--F-On9OFk"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "nskip = 1\n",
        "\n",
        "# Reshape Matrix inputs to vector\n",
        "X = descriptor.forward(qm_coord, atom_types)[::nskip,:,:]\n",
        "shape = X.shape\n",
        "X  = X.reshape(shape[0], shape[1]*shape[2])\n",
        "Y = energy[::nskip]\n",
        "\n",
        "# Initialize likelihood and model\n",
        "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
        "model = BPGPR(X, Y, likelihood)"
      ],
      "metadata": {
        "id": "daVvWWblxpIT"
      },
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Training the GPR Hyperparameters\n",
        "\n",
        "Now we can begin training to optimize the GPR hyperparameters."
      ],
      "metadata": {
        "id": "0yjHzzd29UWG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "training_iter = 500\n",
        "\n",
        "# Find optimal model hyperparameters\n",
        "model.train()\n",
        "likelihood.train()\n",
        "\n",
        "# Use the adam optimizer\n",
        "optimizer = torch.optim.Adam(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters\n",
        "\n",
        "# \"Loss\" for GPs - the marginal log likelihood\n",
        "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
        "\n",
        "for i in range(training_iter):\n",
        "    # Zero gradients from previous iteration\n",
        "    optimizer.zero_grad()\n",
        "    # Output from model\n",
        "    output = model(X)\n",
        "    # Calc loss and backprop gradients\n",
        "    loss = -mll(output, Y)\n",
        "    loss.backward()\n",
        "    if (i+1) % 100 == 0:\n",
        "      print('Iter %d/%d - Loss: %.3f   outputscale: %.3f   lengthscale: %.3f   noise: %.3f' % (\n",
        "        i + 1, training_iter, loss.item(),\n",
        "        model.covar_module.outputscale.item(),\n",
        "        model.covar_module.base_kernel.lengthscale.item(),\n",
        "        model.likelihood.noise.item()\n",
        "    ))\n",
        "    optimizer.step()\n",
        "\n",
        "# Print Trained hyperparameters\n",
        "print('outputscale: ', model.covar_module.outputscale.item())\n",
        "print('lengthscale: ', model.covar_module.outputscale.item())\n",
        "print('noise: ', model.likelihood.noise.item())"
      ],
      "metadata": {
        "id": "WqSWoTHBxz7b",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "6385577b-ed8b-4c61-c755-8be6bba8bbeb"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Iter 100/500 - Loss: 1.700   outputscale: 2.948   lengthscale: 3.301   noise: 0.851\n",
            "Iter 200/500 - Loss: 1.516   outputscale: 4.564   lengthscale: 3.560   noise: 0.192\n",
            "Iter 300/500 - Loss: 1.460   outputscale: 5.940   lengthscale: 3.748   noise: 0.074\n",
            "Iter 400/500 - Loss: 1.400   outputscale: 7.171   lengthscale: 3.976   noise: 0.033\n",
            "Iter 500/500 - Loss: 1.394   outputscale: 8.245   lengthscale: 4.165   noise: 0.018\n",
            "outputscale:  8.256046295166016\n",
            "lengthscale:  8.256046295166016\n",
            "noise:  0.01762973517179489\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Evaluating the Model's Accuracy"
      ],
      "metadata": {
        "id": "OkKS46hBwoEi"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Evaluation of the Predicted Energy (RMSE)\n",
        "We will now calculate the RMSE for the energy and forces predicted by the BP-GPR mode"
      ],
      "metadata": {
        "id": "DHcSnXUqvmit"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "model.eval()\n",
        "likelihood.eval()\n",
        "\n",
        "# Prepare testing set\n",
        "qm_coord = torch.autograd.Variable(qm_coord, requires_grad=True)\n",
        "X_test = descriptor.forward(qm_coord, atom_types)\n",
        "shape = X_test.shape\n",
        "X_test = X_test.reshape(shape[0], shape[1]*shape[2])\n",
        "\n",
        "# Make predictions\n",
        "y_preds = likelihood(model(X_test))\n",
        "\n",
        "y_mean = y_preds.mean\n",
        "y_var = y_preds.variance\n",
        "y_covar = y_preds.covariance_matrix\n",
        "\n",
        "def eval(ref, pred):\n",
        "    rmse = torch.sqrt(torch.mean((ref-pred)**2))\n",
        "    q2 = 1 - torch.sum((ref-pred)**2)/torch.sum((ref-torch.mean(ref))**2)\n",
        "    return rmse, q2\n",
        "\n",
        "print(torch.mean(y_var))\n",
        "# Evaluate performance on energy predictions\n",
        "y_rmse, y_q2 = eval(energy ,y_mean)\n",
        "print(\"Energy Prediction RMSE: \", y_rmse.item(),\"Q^2: \", y_q2.item())"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "OLPc-bFlx6-t",
        "outputId": "29838205-0f71-40ca-a0e7-a3cfff35a1f3"
      },
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.10/dist-packages/gpytorch/models/exact_gp.py:284: GPInputWarning: The input matches the stored training data. Did you forget to call model.train()?\n",
            "  warnings.warn(\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "tensor(0.0675, grad_fn=<MeanBackward0>)\n",
            "Energy Prediction RMSE:  0.050121136009693146 Q^2:  0.9998199343681335\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Evaluation of the Predicted Forces (RMSE)"
      ],
      "metadata": {
        "id": "76P2oVOG9h1f"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "auto_grad, = torch.autograd.grad(y_mean.sum(), [qm_coord])\n",
        "rmse = torch.sqrt(torch.mean((auto_grad - qm_gradient)**2))\n",
        "q2 = 1 - torch.sum((qm_gradient-auto_grad)**2)/torch.sum((qm_gradient-torch.mean(qm_gradient))**2)\n",
        "print(\"Force RMSE: \", rmse.item(), \"Q^2\", q2.item())"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "UMcC9Uv4JhG1",
        "outputId": "d8678f81-9c96-482d-d9b7-f76870020191"
      },
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Force RMSE:  4.071394920349121 Q^2 0.8092787265777588\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Plotting RMSD for the BP-GPR Model\n",
        "\n",
        "Here we plot the RMSD for the predicted and reference energy and forces for the BP-GPR model."
      ],
      "metadata": {
        "id": "LelFo-t29mo0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "fig, ax = plt.subplots(1,2,figsize=(10,5))\n",
        "\n",
        "e1 = energy.cpu().detach().numpy() + np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/energy_sqm.npy\",\"rb\")) * 27.2114 * 23.061\n",
        "e2 = y_mean.cpu().detach().numpy() + np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/energy_sqm.npy\",\"rb\")) * 27.2114 * 23.061\n",
        "ax[0].plot(e1, e2, linestyle='none', marker='.', color='springgreen')\n",
        "ax[0].plot([np.max(np.concatenate((e1,e2))), -np.max(np.concatenate((e1,e2)))], [np.max(np.concatenate((e1,e2))), -np.max(np.concatenate((e1,e2)))], color=\"k\", linewidth=1.5)\n",
        "ax[0].set_xlabel(\"Reference Energy (kcal/mol)\", size=14)\n",
        "ax[0].set_ylabel(\"Predicted Energy (kcal/mol)\", size=14)\n",
        "ax[0].annotate('RMSD: %.3f' % np.sqrt(np.mean((e1 - e2)**2)), xy=(0.05, 0.95), xycoords='axes fraction', size=14)\n",
        "\n",
        "f1 = -qm_gradient.cpu().detach().numpy().reshape(-1) - np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/qm_grad_sqm.npy\",\"rb\")).reshape(-1) * 27.2114 * 23.061 / 0.529177249\n",
        "f2 = -auto_grad.cpu().detach().numpy().reshape(-1) - np.load(ds.open(\"https://github.com/cc-ats/mlp_tutorial/raw/main/Claisen_Rearrangement/qm_grad_sqm.npy\",\"rb\")).reshape(-1) * 27.2114 * 23.061 / 0.529177249\n",
        "\n",
        "ax[1].plot(f1, f2, linestyle='none', marker='.', color='springgreen')\n",
        "plt.plot([-np.abs(np.max(np.concatenate((f1,f2)))), np.max(np.concatenate((f1,f2)))], [-np.max(np.concatenate((f1,f2))), np.max(np.concatenate((f1,f2)))], color=\"k\", linewidth=1.5)\n",
        "ax[1].set_xlabel(r\"Reference Force (kcal/mol/$\\AA$)\", size=14)\n",
        "ax[1].set_ylabel(r\"Predicted Force (kcal/mol/$\\AA$)\", size=14)\n",
        "ax[1].annotate('RMSD: %.3f' % np.sqrt(np.mean((f1 - f2)**2)), xy=(0.05, 0.95), xycoords='axes fraction', size=14)\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.savefig('rmsd.png', dpi=300)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 527
        },
        "id": "NbC3Mh98aEnQ",
        "outputId": "cd1e75db-b772-44e5-b10e-bde161f7b193"
      },
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1000x500 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAHqCAYAAADyGZa5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAD1w0lEQVR4nOzdd3xN9x/H8de9ETMSIsQIsUepqhV7kyqhNjXTKG1RRGu0qE2rtVraIrVj14iqFasIDWr8WtQmtYUbQobc8/vjJJfIzbq5I+PzfDw8Hsm533vO56bl5H2+S6MoioIQQgghhBBCCCEsQmvrAoQQQgghhBBCiMxMgrcQQgghhBBCCGFBEryFEEIIIYQQQggLkuAthBBCCCGEEEJYkARvIYQQQgghhBDCgiR4CyGEEEIIIYQQFiTBWwghhBBCCCGEsCAJ3kIIIYQQQgghhAVls3UB6Z1er+fWrVvkzZsXjUZj63KEEEJkAoqi8OTJE4oWLYpWK8/ArUHu50IIIcwtNfdzCd7JuHXrFsWLF7d1GUIIITKhmzdv4ubmZusysgS5nwshhLCUlNzPJXgnI2/evID6w3R0dLRxNUIIITKDsLAwihcvbrjHCMuT+7kQQghzS839XIJ3MuKGozk6OsqNWgghhFnJkGfrkfu5EEIIS0nJ/VwmlgkhhBBCCCGEEBYkwVsIIYQQQgghhLAgCd5CCCGEEEIIIYQFSfAWQgghhBBCCCEsSIK3DV27dg2NRhPvj729PcWKFaNr164cP37c6Pv69etnaP/DDz8kev5u3boZ2i1dujTea4qisHLlSpo1a0aBAgXInj07rq6uvP3223zyySccOHAgXvsJEybEq9POzo58+fJRvnx5unTpwpIlSwgPD0/zz+RVt2/fxsfHhyJFipAzZ04qVKjA1KlTiY6OTvW5Vq1aRe3atcmTJw/58+enbdu2nDx50mjbkiVLJvjvEvenSZMmRt8TGRnJpEmTKFeuHDlz5qRo0aIMGDCAe/fupbpWIYQQGYvcz1MuKCgIOzs7NBoNM2bMSPX7d+7cSePGjcmbNy+Ojo40bdqUwMDABO3279+f6L087k+ZMmXivefy5ctMmDCBdu3aUaxYMTQaDSVLljT1owohRDyyqnk6UKZMGXr16gVAeHg4J06cYP369WzevJk9e/bQqFEjo+/Lli0bv/zyC4MHD07wWmhoKFu2bCFbtmy8ePEiwesffPABS5cuNYTQYsWK8fz5c06fPo2fnx9hYWE0btw4wfs6depElSpVAHX5/GvXrrF//342bNjA+PHjWbFiRaLhNDXu3LmDh4cHISEhdOjQgXLlynHgwAHGjh3Ln3/+yebNm1O8GvDUqVMZO3Ys7u7ufPTRRzx58oQ1a9ZQr149AgMDqV+/foL3ODk5MWzYsATHjd2A9Xo97du3Z+fOndSpU4dOnTpx8eJFFi9eTGBgIEePHqVgwYKp/REIIYTIYOR+nrRnz57Rt29fcuXKZVK4X7lyJb1796ZgwYL069cPgLVr19KyZUvWrVtH586dDW1LlizJV199ZfQ8e/bs4fDhw3h6esY7/scffzBx4kTs7OyoVKkSd+7cSXWNQgiRKEUkSafTKYCi0+nMfu6rV68qgOLp6ZngtenTpyuA0qhRowSv9e3bVwEULy8vBVBOnTqVoM3cuXMVQGnXrp0CKEuWLDG8dvDgQQVQqlWrZvRzPXr0SDl8+HC8Y1999ZUCKKtXr07QPiIiQpkxY4ai1WqVPHnyKKdPn07Jx09Snz59FED58ccfDcf0er3SvXt3BVD8/f1TdJ5///1XyZYtm1K+fHnl8ePHhuN//fWXkiNHDqVSpUpKTExMvPe4u7sr7u7uKa71l19+UQClR48eil6vNxz/8ccfFUAZMGBAis8lhMgaLHlvEcbJ/Vxl7fv5qwYPHqw4OTkpU6ZMUQBl+vTpKX5vaGioki9fPsXFxUW5efOm4fjNmzcVFxcXxcXFRQkLC0vRuapUqaIAyokTJ+Idv3z5shIUFKQ8e/ZMURRFyZEjR6p+HxBCZD2pubfIUPN0ysfHB4ATJ04k2qZv377Y2dnh5+eX4LUlS5ZQqVIl6tatm+C1oKAgw/uN7WWaL18+6tWrl+Jac+TIwahRoxg/fjzh4eGMHj06xe815smTJ6xdu5bSpUszcOBAw/FXh6UtWrQoRedasmQJL1684Msvv8TJyclwvFq1avTo0YNz585x6NChNNUbV8v06dPj9cIPHDiQ0qVLs2rVKp4/f56mawghhMiYsvL9/FX79u1j/vz5zJo1i2LFiqX6/evXr+fx48cMGTIENzc3w3E3NzcGDx7MgwcP2LRpU7Ln+fPPP/nf//5HtWrVqF69erzXSpcuTZ06dciVK1eq6xNCiORI8E7nsmVLfDZAsWLFaNWqFf7+/kRFRRmOnzx5klOnTuHt7W30fQUKFADg33//NWutI0aMIHfu3OzcuROdTmc4HjfPKqVD1oKCgoiMjKRly5YJhpO7u7tToUIFDh8+TExMTLLn2r9/PwCtWrVK8FrcELPX57+BOmd76dKlTJs2jR9++IFjx44ZPX9ERATHjh2jQoUKuLu7x3tNo9HQsmVLwsPDE53fJ4QQImvIivfzOE+ePMHb25tWrVrxwQcfmFSTqffz1y1evBiA/v37m1SHEEKYSoJ3OhV3Y2jQoEGS7Xx8fHj48CFbtmwxHPPz8yNbtmz06dPH6Hs8PT1xdHTkp59+omfPnmzYsIHr16+nuWYHBwdq1KiBXq9P8sl+ci5evAhAuXLljL5erlw5oqKiUlTzxYsXcXBwoHDhwkbP8+r1XnXnzh28vb358ssvGTJkCHXq1KF27dpcvnw5XrvLly+j1+uTrDWxawghhMj8svL9PM7w4cN59OhRikerGZPU7wYpvdeGh4ezZs0acuXKRc+ePU2uRQghTCGLq6UDly5dYsKECcDLxVj27duHq6srM2fOTPK97dq1w8XFhV9++YUuXboQERHB6tWradOmDa6urkbf4+bmxsaNG/Hx8cHf3x9/f38AChYsSNOmTRk4cCDNmjUz6bMULVoUgAcPHhiO1a5dm3PnzpE7d+4UnSPu6fqrQ8NfFTec7tWn8Emdq1ChQqk6j7e3Nw0bNqRKlSo4ODjw77//MmvWLFasWEHz5s05e/YsefPmNXutQgghMja5nyf0+++/4+fnx88//0zx4sVNqgWSvt+m9F67fv16njx5Qq9evciXL5/JtQghhCkkeKcDly9fZuLEifGOFS5cmD/++IOyZcsm+V57e3t69erFvHnz+O+//zh48CCPHj1KdihXixYtuHz5Mvv37+fgwYOcOHGCQ4cOsW7dOtatW8eYMWOYNm1amj8bQO7cualYsaJZzmUNr6+CWq1aNZYvXw7AihUrWLRoEb6+vrYoTQghRDom9/P4Hj16RP/+/WnevDkDBgwwSw1pETeHPm7evRBCWJMMNU8HPD09URQFRVG4d+8eM2fO5N69e7Rr146nT58m+/4PPvgAvV7P0qVL+eWXXyhcuDDvvvtusu/Lli0bLVq0YNKkSfz22288fPiQRYsWkS1bNqZPn57oPtdJuXXrFkCats+Ke5qd2JPrsLCweO2SO5c5zgMYFno7fPiwRWoVQgiRscn9PD5fX190Op1huH1aJHW/Tcm99sKFCxw6dIiyZcuafZs0IYRICQne6UzBggX57LPP+OKLLzh37hxjx45N9j1vvvkmtWrVYv78+ezdu5c+ffokuYhLYrJly0b//v15//33AXUF0tR4+vQpJ06cwM7OLsFKoamR3Fytixcvkj17dkqUKJGicz19+tToXpzJzSV/nYuLC0C8vUdLly6NVqtNstbUXEMIkX6FoGMfVwlBpo6I5Mn9HP766y/Cw8MpVaoUGo3G8CdusbgxY8ag0WgYNmxYsudK6neDlNxrpbdbiKwtPdzDJXinU1988QVFixZlwYIFXLt2Ldn2H3zwAbdv30av15u8YmgcBwcHk9733Xff8ezZM1q3bp2mHt46deqQPXt2du/ejaIo8V67fv06Fy5coH79+in6ZaRx48YA7Nq1K8FrO3fujNcmOXErm5csWdJwLFeuXNSuXZsLFy4kWNBGURR2795Nnjx5qFmzZoquIYRIn/w4SQllNs02jKJEzCz8SH0PosiasvL9vGPHjvj4+CT406hRIwBq1aqFj4+P0a3SXpeW+/mLFy9Yvnw52bJlo1+/fiZ+GiFEepSSQO3HSdyZQzOW4c4cFisnWL9+PXq93oqVSvBOt3LlysWoUaOIjo5m8uTJybbv1asXmzZt4vfff6dChQpJtt2xYwdbtmzhxYsXCV67dOkS69evB5JfgTVOZGQk33zzDZMmTcLBwYHp06fHe/3Zs2ecP3+eGzdupOh8jo6OdO/enStXrvDzzz8bjiuKwpgxYwD48MMP471Hp9Nx/vx5bt++He+4t7c32bJlY+rUqfGGp506dYrVq1dTqVKleJ/z/PnzPHv2LEFN58+fZ9SoUQCGHoQ4cfPWxowZE+9Bwc8//8yVK1fo2bOn7AkqRAYWgo4BBKBM2g9d1qP0/pUBylbp+RYpkpXv5+PHj2fx4sUJ/sT1eHfs2JHFixfTrVu3ZK/RtWtXnJyc+P777wkJCTEcDwkJ4YcffsDFxYUOHToYrWPbtm3cvXuXd9991+guJ0KIjOn1QG3soXjcPVyP+ju6XtEzYPggunbtyuDBg61aryyulo4NGDCAr7/+muXLl/PFF19QpkyZRNs6ODjw3nvvpei858+fZ/jw4bi4uNCoUSPKlCmDoihcunSJ7du3ExUVxccff4yHh0eC927YsIHz588D6lC0q1evcvDgQR48eEDx4sVZuXIlVapUifeeP//8k6ZNm9K4cWPDPpzJmTFjBvv27eOTTz5hz549lC1blgMHDnD06FG8vLzo3r17vPabNm3C29ubvn37snTpUsPx8uXLM2HCBMaOHctbb71Fp06dePLkCWvWrAFg0aJFaLUvnz+tWbOGWbNm0ahRI9zd3cmTJw///vsv27dvJzo6mjFjxhie1Mfp27cva9euZfXq1Vy9epXGjRtz6dIlfv31V0qVKsWUKVNS9JmFEOnTRULRTzsIE2L3CK5VFL0GLhGKG7J+g0heVr6fp1Zi18ifPz8//PADvXv3pnr16oawvnbtWh4+fMjatWsNO468Lm6YeXJ7dz948IDPPvvM8H10dDQPHjyI10v+7bffGqaeCSGsJwQdFwmlHM644ZQwUKMwkAA8KRPv3nyRUEMbFAVG7EKZq45irVGjhlU/gwTvdCxnzpyMGTOGIUOGMHHiRMPK2mnVs2dPHBwc2LlzJ2fPnmX37t1ERETg4uJCq1at6NevH506dTL63o0bN7Jx40a0Wi0ODg4UKlSIJk2a0KZNG7p27ZqqLUaSUqRIEY4dO8bYsWP57bffCAgIwN3dncmTJzNy5Eg0Gk2Kz/Xll19SsmRJ5syZw48//kj27Nlp2LAhkydPTjB3rWnTppw7d46//vqLP/74g2fPnuHi4sK7777LJ598QqtWrRKcX6vVsmXLFmbMmMGKFSuYPXs2zs7O+Pj4MGXKlDQtTCOEsL3Ab1fBl3vVb2Y0h+F1sUNDWZxtW5jIMLLy/dycevXqhYuLC9OmTWPJkiVoNBpq1KjB2LFjadGihdH33Lp1i99//52iRYsmu1Dd06dPWbZsWbxj4eHh8Y5NmDBBgrcQVubHSUPI1qJhIV6UJv/LQB0rBiXBQ/FyOKNFg17Rw8jdMPsoADN+nmP1NR80yuuTaEU8YWFhhpWx4/aJFEIIkTXMnTvXsPCTZnJTlLGNsEPDz3jhg+mLTsm9xfrkZy6EEBlPCDrcmRMvZNuhIQgf6uCX4Pg1hiUYjbZYOcGAMUNRvlZ3Jur14xhWfGSebRZTc2+ROd5CCCGEEQsWLDCE7nHjxnFj7Cb20ZdrDEtT6BZCCCFEysQbKh4rBoVwolmIF3aoo2DjHoq/HroVReHq2F8NoXvyDzPNFrpTS4aaCyGEEK9ZtGgRgwYNAmD06NFMnDgRDRqZ0y2EEEJYkWGo+Gs922Vxpgml8KQMlwilbOzc71cpisL48eOZNk0N2vPmzWPIoCFWrf9V0uMthBBCvGLp0qUMHDgQAF9fX6ZNm5aqdSWEEEIIYZrXtwdzwynJnm03nGhCKaMPxidOnGhY5HjOnDkMGWK70A3S4y2EEEIYrFq1ig8++ABFURgyZAjffvuthG4hhBDCCowtouZDdXyonmTPtjGTJ09m4sSJAHz33XcMHTrU0uUnS3q8hRBCCGDdunX06dMHRVH46KOPmDt3roRuIYQQwgoS2x7s1Z7vxHq2Xzd16lTGjx8PwMyZM/H19bVc4akgwVsIIUSWt2nTJt5//330ej0+Pj7Mnz9fQrcQQghhJYktonaJ0FSdZ8aMGYwdO9bw9WeffWa2GtNKgrcQQogsLSAggG7duhETE0OfPn1YuHAhWq3cHoUQQghrOcGtBMfiFlFLqZkzZzJmzBhA7fUeNWqU2eozB/nNQgghRJb1+++/07lzZ6Kjo+nRowe//PKLhG4hhBDCikLQMYo9CY5Pp0WKdxOZNWsWI0eOBGDSpEl88cUXZq3RHGRxNSGEEFnS7t276dChA1FRUXTu3Jnly5djZ2dn67KEEEKITC8EHRcJpRzORoeZA9SiaIrONWfOHEaMGAHAhAkTGDdunFlrNRcJ3kIIIbKcffv20a5dOyIjI2nfvj3+/v5kyya3RCGEEMLSXl+9/GtaJLpXd3LmzZvH8OHDARg3bhxfffWVxepOKxlPJ4QQIkv5448/aNu2LREREbRp04a1a9dib29v67KEEEKITM/Y6uWj2cMMWiS6V3di5s+fb9gm7IsvvjBsH5ZeyeN9IYQQWUZQUBDvvvsuz549w9PTkw0bNpAjRw5blyWEEEJkCYmtXl6LolxjWIr36v7xxx8ZPHgwAKNHj2bKlCnpfjcSCd5CCCGyhODgYN555x2ePn1K8+bN2bRpEzlz5rR1WUIIIUSG9up87eQCczmcEx1W7oZTgvcbO/fChQv55JNPAPj888+ZNm1aug/dIMFbCCFEFnDy5ElatWpFWFgYjRs3ZuvWreTKlcvWZQkhhBAZ2qvztTXACOoxFI8kA/RCvBhIADEoSQ4rf30u+EK8UBafZODAgQD4+vry9ddfZ4jQDaBRFCXhEnLCICwsDCcnJ3Q6HY6OjrYuRwghRCqdOXOGpk2bEhoaSv369dmxYwcODg42rUnuLdYnP3MhhDCvEHS4MyfB0HEtGnypawjgrwfoGbTAHSc0QF2KGw3dxs6t+eUv6B+AoigMHTqU2bNnJxq6U9MLnxapubdkqsXVZsyYgUajYdiwYYZjERERDBo0iAIFCuDg4ECnTp24e/eu7YoUQghhNX///TfNmzcnNDQUDw8Ptm/fbvPQLYQQQmQGR7hpdBswPQrfcoQSzGYmhxMspjaS3XRjA93ZyE4uGz13grngy06h9N+KoigMGTIkydDtx0ncmUMzluHOHPw4mfYPawaZJngHBwfz888/U7Vq1XjHhw8fTkBAAOvXr+fAgQPcunWLjh072qhKIYQQ1nL+/HmaN2/OgwcPqFGjBjt27JCeTiGEEMIM/DhJDzYm2UYBRrPbaDgHNYQPJIAQdAlei5sLDsCK0+C9RT3hoFpUnduX/zRh7ONqgvcaWzU9sWtYW6YI3k+fPqVnz54sWrSI/PnzG47rdDr8/PyYNWsWzZo1o0aNGixZsoQjR45w9OhRG1YshBDCki5evEizZs24e/cu1apVY9euXeTLl8/WZQkhhBAZXjAhfMjWRAP1q/TwMkAbEYPCJUIJQRcvSLvhxEK80PifhX6xofujGvB9awZqtlGC2UZ7tBNbNf0SoSZ9VnPKFMF70KBBtGnThhYtWsQ7fuLECaKjo+Mdr1ixIiVKlCAoKMjaZQohhLCwEHSsunKAxs2acvv2bapUqcLu3btxdna2dWlCCCFEhufHSergZzRyG4vXdmjoRKVEz2eHhuPcMjo0PM+af6H3ZtArMKA6zG8DGnVF9Ljrv96jHa+n/JVrlMX2vwdk+FXN16xZw8mTJwkODk7w2p07d8iePXuCXg5XV1fu3Llj9HyRkZFERkYavg8LCzNrvUIIISzDj5N8eH0FSrMlEKKjSKVSBAYG4uLiYuvShBBCiAzv9WHcr7JDQxA+rOMfZnEEfeyx6bRgNHsSPWd/qjOKPQmGhketO8Pgnj4oej30rw4/tgWt8Z7zuB7tuO3IUrpqurVl6OB98+ZNhg4dyu7du822F+v06dOZOHGiWc4lhBDCOkLQ8WHISpRmy+C6DsoX4G5gR6IK5bB1aUIIIUSmYGwYN6hDqH/Gi1q4UQs3huJBEDfj9UobowHewjXh0PANfzP4/cno9Xq8vb2p+/MnfKz9zRCkX+3xhoQ92j5Ux5MyXCLUsD94epChg/eJEye4d+8e1atXNxyLiYnh4MGD/PDDD+zcuZOoqCgeP34cr9f77t27FC5c2Og5x4wZg6+vr+H7sLAwihcvbrHPIIQQIm1C0LHo1n41dF95BGXyw94+6Is4GJ6ACyGEECJt4oZxvxqUtWg4ig+1cDMc28nleHt7a8Bo9FaAQWyP//qv56DHRvQxeur1bcuiRYuw09rRmnKGIL2Ty8n2aMf1fqcnGTp4N2/enLNnz8Y75u3tTcWKFRk1ahTFixfH3t6ewMBAOnXqBMCFCxe4ceMGdevWNXrOHDlykCOH9JAIIURGMJZApt79HZovg4uhUDIf7O0LxRzTzZwuIYQQIiN6fS/sxIZxvxq6Xx+OntzyawpqMNcC+s3nodsGeKGH3lU56ledjXbnqBe713dckE6vPdrJydDBO2/evFSpUiXesTx58lCgQAHDcR8fH3x9fXF2dsbR0ZEhQ4ZQt25d6tSpY4uShRBCmEkX1rHhfjA0Xw7nH0BxR9jbB0o4pas5XUIIIURG48dJQ4DWomEGLahJUTwpQxA+HOImDSgeL3RD4sPRk6IAXQOysbbrejV0v/8mLGmP3k5DNzagRcNCvPDh5Sjn9NijnZwMHbxTYvbs2Wi1Wjp16kRkZCSenp4sWLDA1mUJIYRIg6X8xYaHx6HFCvj7PhTNq/Z0l8rPbDzpzBsZ7oYshBBCpAfG9sIeyW7g5crlCiQIxCHouE94okPLE/Xbv2zotB6i9dC9Cix7D+xebr4Vt+CaJ2Uy9L090wXv/fv3x/s+Z86czJ8/n/nz59umICGEEGbVllX89ugMtFoJZ+5CYQe1p7usMxqQ0C2EEEKkQVK91q8efTUQp3RedwK/X4SO64iJjoEub8CKDpAt4Y7Xr65cnlFlin28hRBCZA1tWMVvurPguRJO3oaCuSGwD1RQtwz7goYZ+qYshBBC2JqxvbATE4NCEDcTzOvWomE8jZJ+885L0GEtRMVAp0qwqiNk0xq9dmZYt0WCtxBCiAxhGxfY/uR/0HoVBN+CArnU0P1GQUDt6Z5CcxtXKYQQQmRscYuo2aUgfNuhQSHhlmExKFShUOIBfvdlaL8GImOgQ0VY3Qns7SD2XJ9Rz3D9zLJuS6Ybai6EECLzGcJ2fgg/BG38ISgE8ueEPX3gTVcAltCefrxt4yqFEEKIzCFu5fAgbtKVDUbbaIHh1KUU+RJsM2aHhroUj7cKukHgFWgXG7rbV4A1nQ2hO+69Q/FgKB4ZbuXypEjwFkIIka7V4meOP7sBXqvhjxvglAN29YZqhQGoQzEJ3UIIIUQqvL5VWGKvJzZPuwklOch1vuUIswiiN1VZyZkEe2u/uvXXOe7zyd756v084gV4lYd1XdBkt0MTG9xf793ODIE7jgRvIYQQ6dZ49nI84ia8twb2XYO82WFnL6hZ1NBmPV1tV6AQQgiRwby+VdjrW3W9+npiC6Xt55rhaz0KKzlDED6EE23ooX413DehFPcP/ANt/eH5C2hTDtZ3gex2sXt5K6yjM3Vj9+zOjGSOtxBCiHSpKguYHLkXOq2D3Vcgjz383hM8Xu4ZOpOWmfYGnZkcPHgQLy8vihYtikajYfPmzfFeVxSF8ePHU6RIEXLlykWLFi24ePFivDahoaH07NkTR0dH8uXLh4+PD0+fPrXipxBCiIzP2FZhAwkgBJ3R15XYP8nN9o5B4RqPDQHdj5O4M4dmLMOdOYz6w49+73ZVQ/c7ZWFDV8jxsg9YDxQkT6a+p0vwFkIIke40Zylno25D1w2w/SLkyga/vQ/1SxjafEFDPqO+DasUKRUeHs5bb72V6Nae33zzDfPmzeOnn37i2LFj5MmTB09PTyIiIgxtevbsyd9//83u3bvZtm0bBw8eZMCAAdb6CEIIkSkY2yosbmXyfVwlgAtGtxJLydZg3dhgCNofsvVluD90nW9af8yzZ89o1aoVf2z6HW1O+3jvzQyrlidHhpoLIYRIV7Zxgb3Rl6HHRth6AXJmg4Ae0Likoc1MWkrozkBat25N69atjb6mKApz5sxh7NixtG/fHoDly5fj6urK5s2b6d69O+fOnWPHjh0EBwdTs2ZNAL7//nveffddvv32W4oWLWr03EIIIeKL2yrs1XCtAbqzMdG9u1Mq7t3xznPkprobSXg0lVrWptnm0eTImTPeomuZZdXy5EiPtxBCiHRjLIF4vVgFvTfBr+cgux1s7gbNSxvaBNBDQncmcvXqVe7cuUOLFi0Mx5ycnPDw8CAoKAiAoKAg8uXLZwjdAC1atECr1XLs2DGr1yyEEBnV61uFxYXBtIZuo46GwDsr4WkUNCvFuc0tGJ3rALVZzB9c5xrD2EdfrjEs3hzzzEp6vIUQQqQLbVjF9pgL4L0F1v4N9lrY2BU8yxradOYN2lLBhlUKc7tz5w4Arq6u8Y67uroaXrtz5w6FChWK93q2bNlwdnY2tHldZGQkkZGRhu/DwsLMWbYQQmRYr640fo9wuiWyXVia/PkfeK6EJ1HQtKQ6ci33y+HlyziNB8XwokKm7+mOIz3eQgghbO5TtrNd/y98GAArz0A2LazrAm3LG9p8SUNZwVyk2PTp03FycjL8KV68uK1LEkKIdMMNJ5pQinoUR5vI0mka4EfeTXZhtQSO3yJ7q9UQFslbjT0YFfBDvNAd5xO2484c/DiZ6vozIgneQgghbOoDNvO9/hh8vA2WnAI7DazuBO9VBNQb/5/0ZwrNbVqnsIzChdX92O/evRvv+N27dw2vFS5cmHv37sV7/cWLF4SGhhravG7MmDHodDrDn5s3b1qgeiGEyNheH3r+KiX29VQNQj95G1quIEoXTsOGDTm0bQ+d8ryVaPPXV1XPzCR4CyGEsJmqLGCJ8hd8+jssPAlaDSzvAJ3fMLRZRDtq4ZbEWURGVqpUKQoXLkxgYKDhWFhYGMeOHaNu3boA1K1bl8ePH3PixAlDm71796LX6/Hw8DB63hw5cuDo6BjvjxBCZGUh6NjH1QQh14fqfI/xBTC9WJ3yC/x1G1osh8cRUL843X6bioODA7Vwoy+Jh+8YFC4RmvLrZFAyx1sIIYRNfMBmzip3wXcnzA9Wu7aXtIf33wSgEi7soneWmfuVmT19+pRLly4Zvr969SqnTp3C2dmZEiVKMGzYMKZMmUK5cuUoVaoU48aNo2jRorz33nsAVKpUiXfeeYcPP/yQn376iejoaAYPHkz37t1lRXMhhEgBP04a9ufWomEhXuZd0Oz0HWixAh5FQF03+L0nQ/PupT1VccOJpXRgELX4jX+ZyMF4b80KW4mBBG8hhBA28Cnb1Z7uUXtgTuyq1Iu8oM/LJ+ISujOP48eP07RpU8P3vr6+APTt25elS5cycuRIwsPDGTBgAI8fP6ZBgwbs2LGDnDlzGt6zatUqBg8eTPPmzdFqtXTq1Il58+ZZ/bMIIUR6FoKOi4RSDmfDPTQEnSF0w8vh3VUpxFOiKYczBcht+kXP3IXmyyH0OXgUgx29IG8OQ092XB21cKMWbhQnHwPYih51+HVW2EoMJHgLIYSwsjasYjsXYfw+mHlEPfhjG/B5+eT9G1pmiZtwVtGkSRMUJfFZghqNhkmTJjFp0qRE2zg7O+Pv72+J8oQQIlNIrFf7CDcTbBcWg4IHi1EALRrG0MC0i/7vnhq6Hz6HWkVhZy9wzAEk15OtQZ1Fnuql2zIsmeMthBDCaj5gsxq6Jx2AKX+oB+e9Ax+93J95CLX5XPbpFkIIIVIssV7tbzlMDzYafU9cFNejMJU/Un/Rv+9Bs2Xw4BnULEqrXV9h55QLUEO3sZ7sYEKM1pkVFleTHm8hhBBWYejpnv4HfLVfPfhdKxjycnGsppRkHu/apkAhhBAig7pIqNFe7VHsSXDcLM7dh2bL4f4zqF4EdvXCL183AC4RStnYoe6vDn3fyWU+ZGuCal4fkp5ZSfAWQghhcYae7u+OwBd71YMzmoNvXUObmhRhL/1sU6AQQgiRgZXDGS2aeCFbC5YJ3RceqKH7XjhUKwy7e0P+XFwilCaUMgTo14e+KyhGq8kqi6vJUHMhhBAWVY/FLOEUzDsGn+1WD05qAqNezieriDPBDLRJfUIIIURG9/p+3HZomEFLtOaeQ/3vQ2i6DO48hbdcYU9vcM6VIDwbG/puLHRrExmSnhlJj7cQQgiLKclsrqODH4Nh6A714NiGMK6xoU1zSrJHerqFEEKIJBlbsfxVPlTHkzLxhno7k4uBBBBjjp7vi7Gh+/ZTeLMQ7OkDBdTV0GfQIl5Nxoa+v04LHMWHWrilvbYMQIK3EEIIi6jE92roXnwSPtmuHhxVHya93FbKmZwSuoUQQohkpHQfbjec4gXguDDem1/Zz3XTC7gcqobuW0+gckEI7AMuL7cgq0nReM2NDX1/VdxnyCqhG2SouRBCCAuYxRHO8xCWnYIBAerB4XVgenPQqMPe8pOTh4y2XZFCCCFEBpDYiuUpXQncDSe+oaXpBVx5pIbu/55AJRc1dBfME69JHuwTXPPVoe8JKXhSxvSaMiAJ3kIIIcyqGUsZwS7wPwveW9T9SgbXUlcwjw3dXXmDUAndQgghRLISW7H8EqGAukXXLI4QTIjh9RB07OOqIZzXwo0uvJH6i197rIbum2FQ0QX29gVXh4TNeJzgmA/VucYwZuGZ4DU9GOrPKmSouRBCCLOpxPdqT/f6v6H3JjV0D6wB81obQvcXNGQqzW1baBawb98+AgMDOXz4MCEhITx48IDcuXNTsGBB3nzzTRo3bkzbtm0pXLiwrUsVQgiRBGPDtuMWM+vHJpZx2nC8L2/REPcEw9I9KUNFXFJ34euPoclSuKGD8gVgbx8onDB0A3RjA2FExRv+HjcnvQHFE60/K9EoimKBNeYzj7CwMJycnNDpdDg6Otq6HCGESLc+YLO6evmmc9BlPcQo8EE1WNQOtGro9qYav/CeLctMFyx1bwkPD2fevHksWrSI69evE3eLz5kzJ87Ozjx//hydToderwfA3t4eLy8vhg8fTv369c1WR3ok93MhREbmx0nDIml2sSuBV6UQtVmcoK0G4vWPxw32TlXou6FTQ/fVx1DOGfb3g6J5sUPDKBownT8SnE8D3GA4bjglmJPem6qs5Ey8+o3NUc9oUnNvkR5vIYQQadaNdazjHwi4AN02qKG7d1VY6GUI3WXJL6Hbgn766ScmTpzI3bt3qVq1KpMnT6Zu3brUrFmTvHnzGtopisLFixc5duwYu3btYsuWLWzatIn27dvz3XffUapUKRt+CiGEEMYYW7F8PHuNtn09EKe6lzUkTB1efvUxlHWGfX2hqHof8aSs0dAdd50gblIXEsxJX8kZgvAhnGhD/VmNWYK3Tqcjd+7c2NvbJ99YCCFEplKEmdwhHHZcgs7rIVoP3avAkvZgpy4l4kxOLjLUxpVmbkOGDKFHjx6MHDmSKlWqJNpOo9FQvnx5ypcvT+/evXn+/DmrVq1i+vTprFixgvHjx1uxaiGEECn1+orlOS3Rh/pfmNrTfeURlM6vhu5iL3tyt3Mx2VMkNic9nGiakHUf7qb6v9bTp09Zv369Yd7YrVu3ePHiBQB58+blzTffpEmTJrRv356aNWuavWAhhBDpRy1+VkP3nivw3hqIioFOlWBFB0Po1oCsXm4Ff//9N+XLl0/1+3LlykX//v3x9vbmxo0bFqhMCCFEWiS2f7fZ50jfeqL2dF9+BKXyqaHbLeVTc7RoqEtxbvNE5nQbkeLgffPmTaZMmcLq1at5+vQpAPnz56dMmTKGeWOhoaEcPXqUw4cPM23aNKpVq4avry89e/a02AcQQghhG5+ynePchv3XoN1qiIyB9hVgdSfI9nLTDD0TbFZjVmJK6H6VnZ2dDDMXQoh0JrH9u+NWK399PrfJbj+BZsvgYii4O6mhu0TKh4PHzdveyeV4w8xffS0rDi9/VYqC9+jRo5k3bx4xMTG0bt2arl27UrduXaM36GfPnnHixAl27dqFv78/vXv3Zvbs2SxatIi3337b7B9ACCGE9c3kMN/zJxy6AW394fkLaFMO1nYGeztDO0VCtxBCCJEqcT3cDtjzIVsNEVaPwodsJZTnjGZPguHcJi2iBnDnKTRbDhceqmF7X19wz5fit39FI96gEKXIRx384tWlBYLwoRZuqa0q00nRqub58uXj008/ZdiwYTg7p26IwJ49e5g0aRItWrTIkPPGZBVUIYSIbxsX8GI1BN2EVivhaRS0KgNbukNO9XmuFoiR0J0oW95bLl68yKlTpwx/Tp8+TUhISPJvzODkfi6EyAhe7eFOTW+2BviGljwigmn8kfIL3gtXh5f/cx+KO6qrl5fOn+DcJFKLBtDEDit/fXh5nH30zbRzu82+qvnVq1fJnz9/8g2NaNGiBS1atODRo0cmvV8IIUT60YaVbOcSBP8H76xSQ3ezUrC5myF0N8CNP+hv40qzhhMnTjB+/HguX75MzZo1mTlzJkWKFAEgMjKSM2fOxAvZZ8+eJTw8HMCw1ZhGo0n0/EIIIawnBF28Ydqp6blWgFFGesGTdD9cHV7+z30ollft6S6dMPMpwDo6U5A87OYKM/gDPepDdgXirV7+Opnb/VKKgrepodvc5xBCCGE7zVnGXq7CX7fVnu6wSGjkDlu7Qy51V4scaCV0W8mVK1do2rQpjRs3ZurUqSxZsoTmzZvz9ttvc+rUKS5evEhMTAxAvP28FUWhcuXK+Pj4ULly5SRXQBdCCGE9xlYDT41UvffBM2i+HP6+r24Vtq8vlDEekLVAXYoDatB+jwqEE809wunGBqPt9cjc7tfJPt5CCCGS9SWBaug+cxdarIDHEVCvOGzrAXmyA+pwswgy3pSijGrbtm2Eh4ezcuVKnJycaN68Oc7Ozly4cAGtVkvZsmWpUqUKVapU4c0336RKlSqUKVMGe3t76taty7Bhw2z9EYQQQryiHM6JDtc2q4exofvsPSjioIbucgUSba4AY9nLCs7EW+TNkzJGVy/P6vt1JyZFwfvgwYMmX6BRo0Ymv1cIIYTtbeOCOl/sn/vQYjmEPgePYvB7T8ibA4DS5OMyw2xbaBZTtGhRFEXhyJEjtG7dmiNHjgDq0PGhQ4cyceJE8uTJY+MqhRBCpJQbTgyilrp4qaWEPlcfoJ+5C655YG9fKJ946AY1eC/jtOF7PQoDCeAaw1iIFwMJIAbF0MMtC6kZl6Lg3aRJE5PngMUNcxNCCJHx9GOTerO98ECdB3b/GdQoAjt6gWMOQ7sDeNuwyqypU6dOtGvXDi8vL9zd3VEUhW+++YbvvvuO2bNn4+/vz6RJk/Dx8ZF53EIIkc4ktjd3SfJZ7qKPnkPLFXDqDhTKo/Z0V3Qx6VQxKFwiFB+q40kZLhEqPdzJSFHwHj9+vNy0hRAii5nFETV0XwpVtxm5Gw5vucKu3pAvp6HdYtrJjdYGNBoNmzdvJjg4mLt379KwYUOcnJz46KOPmDp1KnPmzGHgwIHMmzePb7/9llatWtm6ZCGEEMC3HDYshPbqsO2LhFKepHufTfY4Qg3dJ29Dwdxq6K5U0OTTvbpomhtO8ntACqQoeE+YMMHCZQghhEhPqrKAs9yDq4/UbUZuPYEqhWBPH3DOZWi3hPb0420bVipq1aoV73sHBwemT5/ORx99xMiRI1m/fj2tW7emVatWzJw500ZVCiGEAJjJYUay2/C9HoUBBKCgoABaLNDZqYuAVivgxG1wya0OL3/D9NCtAVk0zQRaWxcghBAifWnOUjV0X3+shu6QMHUo2p7e6g07Vk2KSOhOx9zd3Vm7di2HDh2iZs2a7Ny5k7fffltGsAkhhI2EoIsXuuPoY0N33NdmFRapbv8ZfAsK5ILAPuqD9DTQosGTMmYqMOtI86rmhw8f5tSpU4SFheHo6Ei1atWoX7++OWoTQghhZcGEsJdrathuthyu66CcM+ztA64Ohnb9qMYS3rNZnSLl6tWrx7Fjx1i1ahVjxowhJCSENWvWUK5cOYYPH062bLLBiRBCmJuxOdz92WrdIp5Ewjsr4WiIOlptTx+o6prm08bN75Ye79TRKHGbe6bSkSNH8Pb25tKlS4C6R2jcU/Ry5cqxZMkS6tata75KbSQsLAwnJyd0Oh2Ojo62LkcIISxKwwS4/QQaL4WLoVA6PxzoB24v//2bSUs+Qx6wpoWt7i0RERHMnDmTmTNnEh4eToUKFZg7dy4tW7a0Wg22IvdzIYS1+HGSAQTEm8NdlULUZrH1ingSCa1XweGbkD+n2tP9dhGznNoODdcYJsGb1N1bTAref//9Nx4eHjx79oyWLVvStGlTihQpwp07d9i3bx+7du3CwcGBo0eP8sYbb5j8QdIDuVELIbKCEHQUZzbcfQpNlsH5B1Aynxq6S7y8sf5Jf9kmxAxsfW+5c+cOX3zxBcuXLwfgxYsXVq/B2mz9MxdCZA0h6HBnToK9rb+gIZNJuEVzHdzoSCVGsdt8g8yfRsG7q+CPG+piqHt6Q42iZjl13JZhPlQ3y/kyutTcW0waXzZp0iSioqLYvn0777zzTrzXRo0axY4dO2jXrh2TJk1izZo1plxCCCGElfhxUh3+dj8cmi9XQ3dxR3V4+SuhezHtJHSncwsWLKB9+/YUK1YsyXaFCxfml19+YciQIYwYMcJK1QkhROZ3kdAE87RjUMiVSOw6SgjHCDFf6A6Pgrb+auh2ygG7epkUujXAVnrwnGhA3eYsnGjZMiwNTOrxdnV1pXnz5vj7+yfa5v333ycwMJC7d++mqUBbkyfkQojMzNDTHfpc3af79F0omlft6S7rbGgnPd3mZal7i1arRaPRUL16dd577z3atWvHm2++abbzZ2RyPxdCWIOxHm8tGsOq5Rb1LFoN3fuugWMO2N0baif9IDYxGuAGwyVkJyM19xaTVjXX6XSUKlUqyTalSpVCp9OZcnohhBBWUpzZL/f2PH0XXPOoPd2vhG7p6c44goOD+eKLL4iKimLcuHFUq1aNMmXKMGLECA4cOIBer7d1iUIIkam54cRCvLCL3RbMzlqh+3k0tFuthu682WFnL5NDN4ACXCLUbOUJE4N30aJFOXr0aJJtjh07RtGi5plLIIQQwvw0TFC3GfFcCSdvQ8HYvT0ruBjaLKG9zOPKQGrUqMHkyZM5ffo0V65c4dtvv6VEiRLMmzePZs2a4erqire3N5s3b+bZs2e2LlcIITIlH6pzjWGsozNf0NA6obv9Ggi8Cg7ZYUcvqJO2B+Z2aCiLc/INRYqZFLzbtWvH/v37GTduHBEREfFei4iI4KuvvmLfvn20b9/eLEUKIYQwLw0TXq54+ud/L/f2fKOgoY3s052xlSxZkuHDh7Nv3z7u3r3LkiVLaNiwIRs2bKBjx464uLjg5eXFL7/8wr1792xdrhBCZBoh6JjLMbqz0eiCamYV8QLeWwu7r0Aee/i9J9QrbtKp4oJh3AJqMszcvEya4/3w4UM8PDy4evUqBQoUoHbt2ri6unL37l2Cg4O5f/8+pUuX5s8//8TZOWM/KZE5YUKIzMaeCbwIj4J3/eHgdXXF073xtxmRfboty5b3lsjISHbv3s2WLVvYtm0bd+/eRavVUqdOHQ4dOmTVWqxJ7udCCGuYyWFGsts6F4t8AR3Wwu+XIHds6G7kbtKp7NAQhI8soJZKFl/VvECBAhw9epSRI0eyZs0atm/fbngtZ86ceHt78/XXX2f40C2EEJlNTibx4lk0tFujhu64xVdeCd1DqM083rVhlcKScuTIQdu2bWnbti2KohAUFMSWLVvYunWrrUsTQogMK5gQxhBIIFetc8HIF9BpnRq6c2WD3943OXQD9KKqrOdiYSb1eL8qOjqa8+fPExYWhqOjIxUrVsTe3t5c9dmcPCEXQmQWGiaoQ9Lar4Fdl9V5YLt7x5sH9i7l+I2etisyi5B7i/XJz1wIYSn92MQyTlvvglEx0HkdBPyrhu5t70OzpBe+To4dGq4xTHq6U8niPd6vsre3l61KhBAindMw4eXT8V2XX84DeyV0N6eUhO4M7IMPPjDpfRqNBj8/PzNXI4QQWUMwIdYP3V3Xq6E7ZzYI6JHm0A3qXuOXCJXgbUFpDt5CCCHSNw0TYm/UG2D7xZdPxxuUMLRxJid76Gu7IkWaLV261KT3SfAWQgjT/cEN610sOga6b4AtFyCHHWzpDs1Lm+XUWpBVzC3M5OB96NAhvv32W06fPs2tW7d48eJFgjYajcbocSGEENahYQK80MP7G2HrBfXp+NYe0KSkoc1QPJhDa5vVKMzj6lUrzSsUQohMKgQdFwmlXCoWF2tIieQbmUN0DPTYCJvOvwzdrcqY7fS+1DN8ZlN+DiJ5JgXvFStW0K9fPxRFoXTp0tSuXZts2aTzXAgh0pMifKOG7t6bYOM5yG4Hm7pBi5dPx7vwhoTuTMLd3fRFdYQQIqvz4yQDCECPghYNC/HCh+rJvq8WbnhShp1ctlxxL/TQ89f493LPsmY7vR0ahuIBmP5zEMkzKS1PnjyZ/Pnzs337dmrXrm3umoQQQqRRXqbwNCYKPtgCa/4H9lrY2BXeeXmjli3DhBBCCLWHNy5sAuhRGEgAnpRJUY/vYtpRnNmWKS7uAfr6f17ey1uXS/NpNYBC/D270/pzEEkzKXjfvHkTHx8fCd1CCJEOaZgAegU+DIAVZyCbFtZ1gbblDW0q4iyhO4s4fPgwS5cu5dSpU4YdSN5++2369OlDgwYNbF2eEELY3EVCDWEzzuuLjRkbfh137Di3DEHWrGL00HfzywfoG7rGu5ebQgvMoCU9qMIlQuPt2Z2Sn4MwnUnB293dnaioKHPXIoQQIo3smQiKAh9vgyWnQKsB/47wXkVDm5zYcY5PbVeksJrhw4czb9484nYO1Wg0KIrCiRMn8PPzY+jQocyaNcvGVQohhG2Vwxktmnih0w6NYbExY8OvgXi9w2YXowfvLeB/9uUD9HYVTDqVBviGltSkKHmw5ynRADRBXQ097gGCA/ZJ/hxE2mhNedOHH37Itm3bCA0NNXc9QgghTGTPBF4oevj0d1h4Ug3dKzpAl8qGNs0pyXPG2bBKYS3Lli1j7ty5lCtXjlWrVhkWQr19+zb+/v6UL1+euXPnsnz5cluXKoQQNuWGEwvxwg4NkPzw6wFstXzo9tmqjlqz08DazvEeoKfGVzTmBsP5jPpc5hF18KMZy3BnDn6cxI+TuDOHZiyjDn70pqrRn4NIO40S9xg8lQYNGsTevXsZO3Ysb731VqIbhpcoYaWV/iwkNZuiCyGErWiYoPZ0j9gFs4+qj7eXtIe+1QxtvKnGLzK8PF2wxr2lTp063Lp1i7Nnz+LklPCXJp1Ox5tvvknRokU5evSoRWpIT+R+LoRITgg6w/BrUIde3yecbmywXhF6BfpvVUet2WlgTWfo/IZJp9Ki4Sg+1MKNEHS4MydBb7YeJd7jAzs0BOFDONHxhqEL41JzbzF5KfLq1avj7+9Pnz59Em0j24kJIYTlZY8bXj56jxq6ARZ5xQvdhckjoTuL+fvvv+nfv7/R0A3g5OREp06dWLx4sZUrE0KI9MkNJ9xwije0XAOWmb9tjF6BgQEvQ7d/J5NDN6i983XwYyFelCa/0fnbr4tBIZxowzB0YT4mBe/vv/+eYcOGYW9vT9OmTSlSpIhsJyaEEDZgxwT0AOP3wTdH1IM/tgGfl1t/aIHbfG6L8kQ6p9FobF2CEEKkK68PLVcAq/xLqY9dn2XxX+pUsZUdoWvl5N+X3GljVyYPwsfo/G1jPd4yp9syTErLs2fPplixYhw5cgQ3Nzdz1ySEECIF8jJVDd2TD8CUP9SD896Bj2oa2uQiG88Ya5P6hG1VrlyZjRs3MnnyZBwcHBK8/uTJEzZu3Ejlymn/xU4IITILYyt7W7y3W1Fg8PaX67Msfw+6VzHb6eN6sRfixUACiEExzN8GEhyT4eWWYVLwvnPnDgMHDpTQLYQQNuLENHVV0hmHYPx+9eB3rWCIh6HNUDyYQ2vbFChsbuDAgfj4+FC3bl0mTJhA48aNcXFx4cGDB+zfv5+JEycSEhLCpEmTbF2qEEKkCyHouE94gqHlWlAfdFuCosCQ3+HH42rX+tL20LOqWS8R14vdhFJ4UibBNmLGjgnzMyl4ly1blsePH5u5FCGEECnhw2bCiILvjsCYQPXg9ObgW9fQpiZFJHRncd7e3vz111/88MMPdO3aFQCtVoter/76qCgKQ4YMoW/fvrYsUwgh0oXE5nXboaEDFdnAOfNfVFFg2A6YH/xyUdTeb5n1Eq/3YsfNY3+VsWPC/EwK3sOHD2fEiBFcv34dd3d3c9ckhBAiEVWZz1nuw/fH4LPd6sGJTWB0A0MbB+wJZqBN6hPpy7x58+jSpQtLly7l1KlThIWF4ejoyNtvv03fvn1p2LChrUsUQgibMzavW4uGcTTCiRx8xi7zX1RRwHcnzPtT/X5xu3iLoqZUM9ypSwmm8ke843ZoWE0n6lJcQnU6YVLwLlOmDI0bN6ZmzZoMGzYsye3EGjVqlKYChRBCqEoxm2vo4Kfj8OkO9eDYhjC+saFNTYpI6BbxNGzYUAK2EEIkwdi8bj0KkzhgmQsqCny+G+YcU79f5AUfvG3SqbpQmY+ojRM5Gc1u9Lzs5e6C+eaJi7QzKXg3adIEjUaDoiiMGzcuyVVRY2JiTC5OCCGEKheTiSAG/E7Cx7+pB0fWg0lNDW2+oxW+1LNRhUIIIUTG5IB9ghW/LUZRYNQe+C5I/f7nttC/etLvSUJbKgDwOfXpQRWZq52OmRS8x48fL1uQCCGEldTiJzV0LzsFHwaoB4fXgRktIPbfYnecJHQLo2JiYggJCeHWrVtER0cbbSOj04QQWUkIOi4SSjmc2cnleMPMLUpR4ItAmBm7/eeCd2FADZNP14Zy8QK2zNVO30wK3hMmTDBzGUIIIYwpxAzuEwH+Z8F7izrxbFAtdQXz2NCdm2xcY7htCxXpjl6vZ9q0acydO5fQ0NAk29p6dNqECROYOHFivGMVKlTg/PnzAERERDBixAjWrFlDZGQknp6eLFiwAFdXV1uUK4TIwF5fRM0KcVulKDBuH8w4rH7/Q2v4uFaaTrmdiwQTQi1kp6mMIMXBe+TIkXTo0IG6desm31gIIUSa2TFB3b5k/d/QZ5P628GA6jCvtSF0FyQn9xhtyzJFOjVmzBhmzpxJoUKF8Pb2pkiRImTLZtLzdquoXLkye/bsMXz/aq3Dhw/nt99+Y/369Tg5OTF48GA6duzI4cOHbVGqECKDMraImtVM2A9TYxdAm/sODKqd5lMqQB38WIgXPpg+XF1YR4rvwLNmzeK7776jYMGCeHl50b59e1q2bEmOHDksWZ8QQmRJeZmqhu7N5+H9XyFGAe9q8GNb0KqhWwMSukWili1bRoUKFQgODsbBwcHW5SQrW7ZsFC5cOMFxnU6Hn58f/v7+NGvWDIAlS5ZQqVIljh49Sp06daxdqhAigzK2iJpVTNwPkw6qX8/2hE89zHZqPQoDCcCTMjLMPJ3TprTh7du3WbhwIbVr18bf35/27dvj4uJCp06dWL58OQ8fPrRknUZNnz6dWrVqkTdvXgoVKsR7773HhQsX4rWJiIhg0KBBFChQAAcHBzp16sTdu3etXqsQQqSUF6t4SjRs+xe6rocXeuhdVV319JXQrWeCTesU6dvTp09p06ZNhgjdABcvXqRo0aKULl2anj17cuPGDQBOnDhBdHQ0LVq0MLStWLEiJUqUICgoyFblCiEykBB07OOqYRE1q5pyECbEro7+bUsYZv6HhTEoXCLpKUXC9lIcvAsWLIiPjw9bt27lwYMHbNiwgU6dOnHw4EH69etHkSJFaNKkCXPmzOHKlSuWrNngwIEDDBo0iKNHj7J7926io6Np1aoV4eHhhjbDhw8nICCA9evXc+DAAW7dukXHjh2tUp8QQqRWF9ayjYuw4xJ0WgfReuheBZa0Bzv1n2w7JHSL5FWtWpVbt27ZuowU8fDwYOnSpezYsYMff/yRq1ev0rBhQ548ecKdO3fInj07+fLli/ceV1dX7ty5k+g5IyMjCQsLi/dHCJH1+HESd+bQjGXUwY/eVMXOWuF7+h/qvG6Ar1vACMssgmqHhrI4W+Tcwnw0iqKkabyFXq/n0KFDbN68ma1bt3LlyhU0Gg1vvPEGHTp0oF27dtSsWdNc9Sbp/v37FCpUiAMHDtCoUSN0Oh0FCxbE39+fzp07A3D+/HkqVapEUFBQioanhYWF4eTkhE6nS3SvciGEMIdgQqjNYthzBdr6Q2QMdKoEazpDNjV0Z0NDNF/ZuFKRVta4t/z222906dKFQ4cOUb16xpr79/jxY9zd3Zk1axa5cuXC29ubyMjIeG1q165N06ZN+frrr42ew9iCbYDcz4XIQkLQ4c6ceMPL7dDwPa1xITc/c4JArlrm4l8fgtGB6tfTm8PoBha5TNye3TLH2zZScz9P8yorWq2WRo0a0ahRI2bNmsXff//N5s2b2bJlC1OnTmXq1KkULVqUmzdvpvVSydLpdAA4O6tPfJIbnibzwoQQ6UltFsP+a9ButRq621eA1Z0MobswebjN57YtUmQYbdq0YenSpbRu3Zp27drx1ltvJfpLQZ8+faxcXdLy5ctH+fLluXTpEi1btiQqKorHjx/H6/W+e/eu0TnhccaMGYOvr6/h+7CwMIoXL27JsoUQ6YyxOd0xKHzCdrRoaB+7B7bZfXvkZeie0tTsoVsDTKEZ5XCmLsVlbncGYfblTStXrkzlypX58ssvuXXrFlu2bGHr1q3mvkwCer2eYcOGUb9+fapUqQJg0vC0yMjIeE/VZWiaEMLSDD3dh26oPd3PX8C75WBtZ7C3A8ABewndIlUiIyMJCAjgwYMH+Pn5AaDRxB9eqSgKGo0m3QXvp0+fcvnyZXr37k2NGjWwt7cnMDCQTp06AXDhwgVu3LiR5E4rOXLkkAVghcjiyuGMFo3RBdX0KGzivPkvOjsIPt+tfj2xCXzZyKynj/tX/Ev2okUjK5pnIBbdV6Ro0aJ8/PHHfPzxx5a8DACDBg3if//7H4cOHUrTeaZPn250aJoQQlhCW1bxGxfhaAi8uwrCo6FVGdjYFXKo/0SXJR8XGWbbQkWG4+vry6pVq6hatSqdO3dO19uJffbZZ3h5eeHu7s6tW7f46quvsLOzo0ePHjg5OeHj44Ovry/Ozs44OjoyZMgQ6tatKyPXhBCJCkHHRUL5mhaMZg8x1ljNfO5R8N2lfv1VYxjf2Gyn1gIfUoNFnDQ8SJAVzTOWFN2BDx48aPIFGjUy71MeYwYPHsy2bds4ePAgbm4vN5AvXLhwqoenydA0IYS1NGcpe7kGx2+B50p4EgXNSsHmbpBT/ec5J3YSuoVJ1q9fT40aNQgKCkq3gTtOSEgIPXr04OHDhxQsWJAGDRpw9OhRChYsCMDs2bPRarV06tSJyMhIPD09WbBggY2rFkLYUlywLodzgtDpx0nDft1aNLxDWbZz0bIFfX8Mhu1Uvx7bUA3eZuJJGRbTjouE8jMn4r0Wt6K5BO/0L0V34iZNmiQYnpZSMTExJr0vJRRFYciQIWzatIn9+/dTqlSpeK+bMjxNhqYJIazhSwLV0P3XbWi5AsIioZE7bO0OuewN7Z4zznZFigwtIiKCpk2bpvvQDbBmzZokX8+ZMyfz589n/vz5VqpICJGevR6sXx1uHYLO8BqovcIWD90LguHTHerXYxrApKZgYnYyZieXuc0To0PnZUXzjCNFd+Px48ebHLwtadCgQfj7+7Nlyxby5s1rmLft5ORErly5ZHiaECJd2sYFpvEHnL2rhu7HEVCvOGzrAXmyG9opsmWYSIMaNWpw6dIlW5chhBBmZSxYvzrc2tiCahb183EYtF39emQ9mNrMrKE7zmFuMoy6LMSLgQQQg2JY0Vx6uzOGFAXvCRMmWLgM0/z444+A2iP/qiVLltCvXz9AhqcJIdKXfmxiGafhn/vQfDk8fA61i8H29yHvy9E2ErpFWk2bNo3mzZuzbds22rZta+tyhBDCLBJbqTyIm7gQyh6uWK+YRSfgo9/Urz+rCzNaWCR0A9RHnfrqQ3U8KcMlQilrZJi9SL/S//izJKRkC3IZniaESC+CCVFD94UH0GwZ3H8G1YvAzl7glNPQTkK3MIfdu3fTpEkT2rdvT7NmzRLdTkyj0TBunExpEEJkDMaGW2uA7my0bk+330kYsE39engd+KalxUJ3X96iFi/XsXLDSQJ3BqRRUpJes7DUbIouhBBJ0TABLoVC46Vw6wm85Qp7+4JzLkMbCd1ZgzXuLVqtNkXtNBqNRddjSS/kfi5E5vHqHG+bWHoKPtgCCjDUA2Z7mj10a4FxNKYN5eKFbpG+pObeYnKP95MnT/jhhx/Ys2cPt27dirf3dRyNRsPly5dNvYQQQmQaGibA1UdqT/etJ1C5IOzuLaFbWMy+fftsXYIQQliEJ2VsF7qXn34ZugfXskjoBujEG0ygqdnPK2zHpOB9//596tWrx+XLl3F0dDQk/aioKJ4/fw6oe3jb29sncyYhhMjcQtBRnNlwQwfNlsPNMKjoAoF9oGAeQ7ubDLdhlSIzatzYfFvZCCFEenKEm7a58Moz0G+zGro/qQnzWltsePmvnCMEnQwpz0RSNg7tNRMmTODy5cssX76cR48eATB8+HDCw8M5duwYtWvXpmTJkvz9999mLVYIITISP06qofu/MGi6DK49hnLOsLcPuDoY2s2kpdxYhVksWLCA//77z9ZlCCGESULQsY+rhKCzdSkJ+Z+FvpvV0D2wBnz/rsVCN7zcn1tkHiYF7+3bt9O8eXN69eqVYJuxWrVq8fvvv3Pt2jUmTpxoliKFECKj2cYF+rMVbj9RQ/eVR1A6vzqnu0heQ7tvaMln1LdhpSIzGTx4MCVKlKBWrVpMnTqVs2fP2rokIYRIET9O4s4cmrEMd+bgx8lE29ajOFbd6Hjt/6D3JtAr0P9tWNAGtJatQPbnznxMCt63b9/m7bffNnxvZ2dnGGIOkD9/flq3bs26devSXqEQQmQwbViFF6vh7lN1ePnFUHB3Unu63V4uvHGT4XwuoVuYUXBwMF988QVRUVGMGzeOatWqUaZMGUaMGMGBAwfQ6/W2LlEIIRJIbG/uuJ5vYz3hA6hhneLW/w09f1VD9wfV4Gcvi4XuuLPK/tyZk0lzvJ2cnIiOjjZ8nz9/fkJCQuK1cXR05O7du2mrTgghMpi2rGI7F+HBM2ixAs4/UMP23r7gns/QThZSE5ZQo0YNatSoweTJk7l27RqbNm1i69atzJs3jzlz5uDs7Ezbtm1p3749rVq1Infu3LYuWQghktybO4zLhlCuRUNvqrKCM9ZZXG3jP9BjI8Qo0K8aLGpn0Z7utXSmIHlkf+5MyqTtxOrWrUvhwoXZtGkTAJ6enpw+fZq///6bAgUK8Pz5c9566y20Wi3nz583e9HWJNuPCCFSKpgQarMYQp9D8+Vw6g4UcYCD3lBWHS6WDYiW0J3lWfveEhoayrZt29i8eTO7d+8mPDycnDlz0rx5czp06EDbtm0pVKiQxeuwJbmfC5F+haDDnTkJwnRcxLXJ+uWbzkHXDfBCD72rwpL2YGfSYOEUsUPDNYZJ4M5gUnNvMen/nlatWhEYGMizZ88AGDhwIPfu3eOtt96iS5cuVKlShcuXL9OvXz9TTi+EEBlSbRbD4whotUIN3a55YF9fQ+jOjkZCt7AJZ2dn+vTpw6+//sqDBw/YunUrPXv25Pjx4/Tv35+iRYvSoEEDW5cphMii3HBiIV4JgomCjUL3lvMvQ3fPN80eujXAR9TALvbRggwtzxpM6vG+ffs2Bw8epHnz5ri4uADw3XffMWXKFHQ6Hbly5eKTTz5hxowZ2NnZmb1oa5In5EKIlNAwAcIi1dB97D8omBv294M3Chra/El/auFmsxpF+pFe7i2KohAUFMSWLVvYunUr586ds1ktlpZefuZCCONC0LGQE0zmoG0LCbgAndZBtB56VIEVHSzS072PvpTFmUuEytDyDCw19xaTgndiYmJiePDgAYUKFUqw2nlGJTdqIURyNEyAp1HguRKO3ATnXGpPd1VXQ5u+vMVSOtiuSJGuyL3F+uRnLkT69S2HGcWeZOdta7BwD/j2i9BhLUTFQLfKsLIjZDN/6JZh5ZlHau4tJi2ulhg7OztcXV2TbyiEEJmAYU53eBS08VdDd76csKd3vNC9hPb04+0kziRE2n3wwQcmvU+j0eDn52fmaoQQImVmcpiR7E62nRaw6L4MOy69DN1d3jBr6C5Lfi7zCAUZVp6VmRS8Dx8+zMaNGxk5ciSFCxdO8Prt27eZOXMmXbt2pU6dOmkuUggh0pt+bGIZp+F5NLRbAwevg2MO2NUL3i5iaPcu5SR0C6tYunSpSe+T4C2EsJUQdIxKQegGC4funZfgvTVq6O5UCVaZt6d7H/0AZFh5FmdS8J41axZnzpxh1qxZRl8vUqQI27Zt47///mPt2rVpKlAIIdKbYELU0B3xAt5bC3uvgkN22NETahUztKtJEX6jpw0rFVnJ1atXbV2CEEKkykVCjQ4d16JBQbHOwmp7rqj38sgYeK8irO4E9uZdo2onl/GhugTuLM6k4B0cHEzz5s2TbNOoUSN2707ZEywhhMhIarMYIl9A53Ww6zLktoffe0Ld4oY2b1KIYAbasEqR1bi7u9u6BCGESJVyOKNFk2Bud0cq8Q5lGUgAMZaM33uvgtdq9UF6uwqwtrPZQzfAQALwpIwE7yzOpDEU9+7do1ixYkm2KVy4MPfu3TOpKCGESK/smADRMdBtA/x2EXJlg9/ehwYlDG0C6MEZPrFdkUIIIUQG4IYTX9MiwfFNnKMqhfCnE8PwsMzF91+Dtv5q6G5bHtZ1huyW2Y0pBoVLhFrk3CLjMKnHO1++fNy4cSPJNtevX8fBwcGkooQQIj3SMgHlhR7e/xW2XIAcdrC1BzQpaWgjW4aJ9Chu15HIyEijr5coUcLocSGEsLQaFE1wLAYFDxZbrq/7wDV1UdTnL+DdcrChC+Qw65rT8WiBsjhb7PwiYzDp/7A6deqwadMmbt68SfHixRO8fuPGDTZv3kyzZs3SXKAQQthaCDqKMxti9NB7E2z4R30qvrk7tChtaNeFNyR0i3TlxIkTfPHFFxw8eJCoqCijbTQaDS9evLByZUKIrCwEHRcJpRzOlMPZ6DZhFgvdf1xXQ/ezaHinLGzsatHQDTCDljLMXJg21NzX15dnz55Rv359li9fzu3btwF1NfNly5ZRv359nj9/zogRI8xarBBCWJsfJ1+Gbu8tsOZ/YK9Vn46/U9bQzptqrKOrDSsVIr5Tp07RsGFDgoKCaNWqFYqiULVqVVq1aoWLiwuKotC4cWN69+5t61KFEFmIHydxZw7NWIY7c+jGeussogZw+Aa86w/h0dCqDGzqBjkt2dOtYSYt+Zz6FruGyDhM+j+tUaNGzJo1ixEjRuDt7Q2oT8wVRf1ro9VqmTt3Lo0aNTJfpUIIYWUh6OjPVtArMCAAVpwBO426+IpXBUO7erjxC+/ZrlAhjJg8eTIAx44do1KlSmi1Wjp06MD48eMND8c3bNjAL7/8YuNKhRBZRQg6BhBgWExNj8IRQqxz8aCb8M4qeBqljlbbbJnQrQG+piW1KCpbh4l4TN6gbujQoZw8eZKBAwdSvXp1SpcuTY0aNfj444/566+/GDRokDnrFEIIqyvObFAU+OQ3+OUUaDXg3wk6VDK0KUM+DtPfdkUKkYhDhw7Rrl07KlV6+f9r3APyXLly8cMPP1C0aFG++OILW5UohMhiLhKaYAVzqzgWAp4r1dDdrBRs6Q657C1yKQUYw554oTsEHfu4Sgg6i1xTZAxpesxTtWpVFixYYK5ahBAi3dAwQQ3dn/4OP59QH2Gv6ABdKxvaVKUQp2X1cpFO6XQ6Spd+uQaBvb09T58+NXyv1Wpp0qQJq1evtkV5QogsKLHtwyzqz/+g1Up4EqUuhrq1u7oNqAXFrWLuhhN+nDT08mvRsBAvfKhu0euL9MmkHu8jR46kqN3s2bNNOb0QQtiUIXSP2AU/BKuhe0l7eP9NQ5sSOEroFulaoUKFePTokeH7woULc/HixXhtIiIiePbsmbVLE0JkQXELqs2gBXZoALBDQz1LLkp6/Ba0WgFhkdDIHbb1gDzZLXe9WHZoKIuz0aH1AwmQnu8syqTg7eXlxblz55JsM2fOHD777DOTihJCCFvJFhe6xwTC7KPqwYVe0LeaoU1BcnEdX5vUJ0RKvfHGG1y4cMHwff369dm1axdBQUEAnDt3jnXr1lGxYkVblSiEyCJeXVBtNHsYTQPG04jvac0YGlKYPOa/6Mnb0HIF6CKhQQn47X2zh+5pNGUfffmGlvEeJvyMF244GR1aL3t6Z10mDTXPlSsX77zzDkeOHKFYsWIJXp83bx6+vr7Uq1cvzQUKIYS1ODGNGICv9sPXh9WDC96F/i+HhOXGnnuMskV5QqRKmzZtGD58OLdv36ZIkSKMGjWKTZs20aBBA5ydnXn06BF6vV7meAshLMpYr+9U/rDsRU/dgRbL4XEE1CsO298HB/P3dLegDLVwowml6EEVLhEab263saH1cb3hIusxqcd7x44dhIWF0apVK0JD4z+x+f777xk2bBh16tRhx44dZilSCCEsLTdTCCMKJh+AyQfVg3PfgY9rGdoUw4FwvrRRhUKkzkcffcR///1HgQIFAHjrrbcIDAzknXfewcXFhRYtWhAQEECHDh1sXKkQIjOz+oJqp+9A8+XwKALquMHvPSFvDotcKpxow9duONGEUvFWMXfDiYV4Ge0NF1mPST3eVapUISAgAE9PT9q2bUtgYKBhhdShQ4dSu3Ztdu7ciYODg7nrFUIIs3NiGs95AV8fgvH71YPftoRPPQxtHLAnBJk+IzIOe3t7XF1d4x2rV68ev/32m40qEkJkNSHouE+49RZUO3tXDd2hz6F2MdjRExwtE7q1Key59qE6npRJ0Bsush6TVzVv0KABq1evplOnTnTu3JlWrVoxfPhwatasyc6dO8mbN6856xRCCIvIySQi0cOsIBgdqB6c1gxGvJwqUxInrjLcRhUKIYQQGY8fJ/mQrYa4rQHLRu+/76mh++FzqFkUdvYCp5wWu9zXtEhxiHbDSQK3MH0fb4B27drx008/8fvvv+Pr60uNGjXYs2cPTk7yP5YQIv0rxxw1dH9/TF3BHGBiExjT0NDGmVwSukWGdPjwYXx9fblz547R12/fvo2vry9Hjx61cmVCiMwuBB39XwndoIZuDeBuiQD6z31othzuP4MaRWBXL8hnmdCtBb6hJZ9R3yLnF5lXinq8Dx48mOhr5cqVo127dhw6dIixY8dy6tSpeK83atQoTQUKIYQl5GAiUSjw03H4NHY9ii8bwriX/2Y5kp2HspCayKBmzZrFmTNnmDVrltHXixQpwrZt2/jvv/9Yu3atlasTQmRmUzCeHRTgurm30jr/AJotg3vh8HZh2NUb8ucy7zViadFwFB9qWXILNJFppSh4N2nSBI1Gk2QbRVHo2LFjguMxMTGmVSaEEBbiwBQ1dP/yF3wcO991ZD2Y3BRi/63TAjpktWeRcQUHB9O8efMk2zRq1Ijdu3dbqSIhRFYQgo6fOWGdi114AE2Xwd1weMsVdvcGZ8uEblBXZH91QTUhUiNFwXv8+PHJBm8hhMgICvE14byA5aeh/1b14DAPmNHCELodyMYTxtqwSiHS7t69e0a3/HxV4cKFuXfvnpUqEkJkBXM5Zp0LXXyohu47T+HNQrCnDxTIbdFLaoA82Fv0GiLzSlHwnjBhgoXLEEIIyyvCTO7zHFafBe8t6pi3QbVglqchdGdHI6FbZAr58uXjxo0bSba5fv267EAihDCbEHTMIsjyF7oUqobu20+hSiEI7AMulg3doP7aUAc/FuKFD9Utfj2RuaRpcTUhhMgoavMTdwiH9X9D702gV2BAdZjX2hC6NUAkX9m2UCHMpE6dOmzatImbN28aff3GjRts3ryZevXqGX1dCCFSyyp7dl+ODd3/PYE3Cqqhu2Aey17zFXoUBhJAiLnnqotML0XBOzw8PM0XMsc5hBDCFFomEMwd2Hwe3v8VYhTwrgY/tgXty2k0eibYrEYhzM3X15dnz55Rv359li9fzu3btwF1NfNly5ZRv359nj9/zogRI2xcqRAiMwhBhx8nLXuRq4/U0B0SBpVcYG8fKGT+0N2LNwmgB+NpxDA8Erweg8IlQs1+XZG5pWioealSpfj888/55JNPyJMndf9zBwUFMXHiROrXr8+4ceNMKlIIIUxViBnqs/dt/0LX9fBCD72qwiKveKFbkdAtMplGjRoxa9YsRowYgbe3NwAajQZFUXujtFotc+fOld1HhBAmC0HHRUIJ5ArT+MOyfd3XHquh+2YYVHSBvX3B1TJTZVZxFn/+hx4FLZoEe5DboaEszha5tsi8UhS833//fcaOHcukSZPo2LEjnTt3pk6dOhQsWDBB2xcvXnD69Gl2797NqlWr+Oeff3jjjTdo3bq12YsXQoiklGUu94mAnZeg0zqI1kO3yrCkPdi9HPAjoVtkVkOHDqVp06b89NNPBAcHo9PpyJcvH7Vr1+ajjz6iSpUqti5RCJFB+XGSAQRYfmg5wPXHaui+roPyBdSe7sKWW59CAZTYz/UyfCvoUUP3z3jhZon9yEWmplHiHn0n4+LFi0yePJn169cTFRUFQLFixXB1dSVfvnxEREQQGhrK1atXiYyMRFEU3njjDXx9fenXrx9abcacTh4WFoaTkxM6nQ5HR0dblyOESKGyzOUyjyDwCrRdDREvoFMlWN0J7O0M7SR0C1tIT/eWmJgY7Ozskm+YwaWnn7kQGV0wIXiw2BqRG27qoPFSuPoYyjrD/r5QzPp/h9fRmYLkoSzOErqFQWruLSkO3nFCQ0Px9/cnMDCQI0eOcP/+fcNr9vb2VKlShcaNG9OhQwcaNmxo2idIR+RGLUTGk4cpPOMFHLgGrVfB8xfQrgKs7wLZJXQL27PGvWX+/PkMGjQoyTYxMTF0796d9evXW6SG9ETu50KYh1V7ukPCoMlSuPwIyuSH/f3AzbJ/f7XE9Xi/ZIeGawyTwC0SsGjwfl10dDQPHz4kV65cODllvv8Z5UYtRMZShJnq6uWHb4DnSgiPhnfLwa9dIYc6uyY7ECmhW9iQNe4tdnZ2rF27ls6dOxt9Xa/X061bN3799VdiYmIsUkN6IvdzIdIuBB3uzLFO6P4vTB1efjEUSudXe7qLWy5raIAR1GMoHuzkMgMJIAbFMLRctg8TxqTm3pLm8d/29vYULlw4U4ZuIUTGUo/Faug+GqL2dIdHQ8vSsPFl6LZDQrfIGho0aEDv3r3Zt29fgtf0ej09evRg48aNyfaKCyGynhB07ONqgi2zrLJdGMDtJ9BsuRq6S+aDfZYN3aD2cM8miNs84TwPDJ/TKp9XZAkZc+K1EEK85lO2E0QIHL8F76yEJ1HQtCRs7g451dCdDXghoVtkEQEBAZQvX54OHTpw6tQpw3G9Xk/Pnj1Zv349n3zyCfPmzbNdkUKIdMePk7gzh2Ysw5058bYIK4czWjRJvNsM7jxVe7r/fQglnNTQXcI6HXwxKHiwmG85YojbCsi+3cIsJHgLITK8IWzne/6EU3eg1QrQRULDEhDQA3LbA1CY3ERL6BZZiKOjIzt37iR//vy0bt2aK1euoCgK77//PmvXruWjjz7ihx9+sHWZQoh0JARdvPnbehQGEsCP/Mk6/sdq/mdY7dsi7j6FZsvgwkMo7qgOLy+Zz3LXM8LYp5N9u4U5pGg7MSGESK9q8TPHuQ1n70KL5fAoAuq6wW/vQ57shna3GWnDKoWwjcKFC7Nz504aNGhAq1atePvtt9m4cSMffvghCxYssHV5Qoh0xthQ8hgUPmG75S9+L1wdXn7ugbqA2v5+UCq/5a8bSwvoE3lN9u0W5iA93kKIDKs436mh+5/70Hw5PHwOtYvB7z0hbw4gbnXSCTatUwhbKl++PNu3b+fevXv8+uuv9O/fn59//tnWZQkhrCyxeduvvuaAvaUHkht3P1y9j/9zH4rlVYeXl7Z86NYCM2nJPvpylP5Gh9HLvt3CXKTHWwiRIS3lL0J4AhceqDfr+8+gehHY0ROcchraxUjoFlnEpEmTkny9du3anDp1iqJFi8Zrq9FoGDdunKXLE0LY0KtbgGnRsPCVVbpffU2D8aHWFvXgGbRYAf+7B0Uc1NBd1vK9y1o0HMWHWrgZji3Ey7CauRbwjV3lXEK3MIc0byeW2cn2I0KkP0v5C2+2wKVQaLwUbj2Bt1whsA8UyG1oJz3dIr2yxL1FqzVtEJtGo5HtxITIxIxtARa3LzVgve3BjHn4TH14fvouFHZQ53RXcLHY5bRo0CezRVgIOi4RSlmcJXCLZKXm3mJSj/exY8fw8PAwqTghhEiLksziOmFw9ZG6AMutJ1C5IOzubQjdubDjGdKDJ7IWY9uGCSFEYvO2LxGKgg23ywp9Di1XqKHbNY/a023B0B1AD6pRONlQ7YaTBG5hESYF77p16/Lmm2/y4Ycf0qtXL/Lly2fmsoQQIqEP2KyG7hs6dQGWm2FQ0UXt6S6YBwB7kNAtsqTGjRvbugQhRDoUtwXY6z3ecYuFvf6aVTx6ru5C8tcdKJQH9vZV7+cW4kkZ2lIBQEK1sBmTxqX16tWLS5cu8emnn1K0aFH69OnDH3/8Ye7ahBDCYBZHWMIp+C9M7em+9hjKOcPePuDqYGgXJcPLhRBCCAM3nFiIF3axC4e9uljY669ZZWG1xxHguRJO3IaCudWH528UtNjltGiYTNNEF5YTwlpMCt7Lly/n1q1bfP/991SsWJGVK1fSpEkTKlasyHfffceDBw/MXacQIgsryWxGsAtuP1F7ui8/Ulc73dsXiuQ1tLvJcBtWKYSwhvnz51OyZEly5syJh4cHf/75p61LEiJDeHVv7lf5UJ3RNACssLCaLjZ0B98Cl9jQXaWQxS6nRUNvqlIHP5qxDHfm4MdJi11PiKSYvJ2Yk5MTgwYN4uTJkxw/fpwBAwZw9+5dPv/8c9zc3OjWrRt79uwxZ61CiCyoG+u4jk7d37P5cvj3Ibg7qT3dbi8XsVhMOxk+JrK0d955h+DgYJPeGx4ezowZM5g/f76ZqzKvtWvX4uvry1dffcXJkyd566238PT05N69e7YuTYh0KwQdAwgwhGoFGEiAofd3JoeZihVGroZFwjur4M//wDkX7OkNb7pa7HI/8i5H8WEFZ+I9dHj1swthTWbZx7t69er8+OOP3Lp1i6VLl+Li4sKGDRvw9PSkdOnSfPPNNzx58sQclxJCZCGfsp11/BO71chyOPdADdt7+4J7PkO7mww3ujKpEFnJ/fv3qVOnDk2bNmXJkiXodMn/Ynn06FEGDx6Mu7s7kydPxtXVcr8Em8OsWbP48MMP8fb25o033uCnn34id+7c/PLLL7YuTYh0K6nF1ULQMYrdli/iSSS0XgVHQyB/TrWn+63CFr1kRQrylOhEP7sQ1ma2fbwfPXrE8uXLWbx4Mbdu3UKj0VC/fn3OnTvH6NGjmTNnDlu2bKFWrVrmuqQQIhPrwlo2cO7lqqdnY/f33NtHHWYeS7YME0J14sQJli1bxsSJE/Hx8eHDDz+kQoUK1KhRA1dXV/Lly0dERAShoaFcuHCB48eP8+TJE+zs7OjevTtTpkyhRIkStv4YiYqKiuLEiROMGTPGcEyr1dKiRQuCgoIStI+MjCQyMtLwfVhYmFXqFCK9SWpxtYuxK5tb1NMoeNcfjtyEfDlhTx+oZtnQndTica++JoQ1pTl479u3j0WLFrF582YiIiIoWLAgn3/+OQMHDqR06dJERkbyyy+/MHLkSIYMGcLRo0fNUbcQIhMLJkQN3Y8j1FVPT91RtxrZ2xfKFQDUBWD0ErqFiKdv37706dOH7du3s2TJEvbv38/KlSsTtNNqtVStWpUOHTrQv39/ihQpYoNqU+fBgwfExMQk6JV3dXXl/PnzCdpPnz6diRMnWqs8IdKtuAXUBhJAzCt7WMdNz7LoqubhUdDGHw7dAKcc6vDy6pb990aLhhm0MHy+pD67ENZkUvC+e/cuS5Yswc/PjytXrqAoCo0bN+ajjz6iY8eO2NvbG9rmyJGDjz/+mEuXLqX7uWNCiPShNotj54LFrnoatwDLK1uN3JCF1IQwSqPR0KZNG9q0aQPAuXPnCAkJ4eHDh+TKlYuCBQtSuXJlnJwy9y+eY8aMwdfX1/B9WFgYxYsXt2FFQthGCDpKk58gfAgnOt4e1m44MYYGlpnj/Swa2q6Gg9fBMQfs7g01ipr/OrEaUYI/uIEehVHsIT+58KE6PlTHkzIJ9u8OQcdFQimXxJ7eQpiTScHbzc0NvV5P/vz5GTZsGAMGDKBChQpJvqdgwYJERUWZVKQQImsIQUdxZscOS1sFx15ZgKXyy1VPZSE1IVKuUqVKVKpUydZlpJmLiwt2dnbcvXs33vG7d+9SuHDCYas5cuQgR44c1ipPiHTJj5MMIAA9Clo0fE0LQ9+2G04EE0Ioz81/4WfR4LUa9l+DvNlhVy+oVcz813nFQW4Yvo5bRM2TMoZt0179veH1n8tCvGStGGFxJi2u5uHhwbJly/jvv//47rvvkg3dAKNHj0av15tyOSFEFvAth9XQHR4Fbf3hcOxcsN294y3AEkAPuTkKkQVlz56dGjVqEBgYaDim1+sJDAykbt26NqxMiPQpbjXzV1f0/pzdhm216rOY2izmR46b98LPo6H9Gth7FRyyw85e4OFm3mukQGKLqBn7uchK58IaTOrxPnTokLnrEEJkYZ+yne/58+XN+kDssLRdveLNBevMG7Ql+Qd9QojMydfXl759+1KzZk1q167NnDlzCA8Px9vb29alCZHuGFvNPI4ehSOEmP+iES/gvbWw5wrksYcdPaGubaZ4JLaIWlKrvMtoOmFJZlvVXAghTNGGVWznonqz7rAWAmOfkO/oGW9Y2hBqM493bVipEMLWunXrxv379xk/fjx37tyhWrVq7NixI91vgyaELRhbzdyi4u7juy6rofv3nlDfOjsl2KGhF1VZyZlkF1FLapV3ISzJpODdrFmzZNtotVocHR2pUKEC7733Hh4eHqZcSgiRiY0lUA3dUTHQeR3svAy57WH7+/GekA+W0C2EiDV48GAGDx5s6zKEyBB8qctsgoixdPiOfAGd1sGOS7H38Z7Q0N2y10QNzKvpRF2K44YTU2iWYBG11yW3yrsQlqJRFCXVfxO1WnVquEajwdjbXz+u0Wjw9vZm8eLFaSjVNsLCwnByckKn0+Ho6GjrcoTINAwLqUXHQNcNsPk85Mqm3qyblDS0a0M5ttHTdoUKYQFyb7E++ZmLrOTVxcM0wAjqUYg8jGaP+XvAI19A5/Ww7V+j93FLGkgNxtLIpNAcgi7ZkC5EclJzbzFpcbXnz5/j5eVFpUqV8Pf35/r160RERHD9+nX8/f2pXLky7dq14+bNm+zatYvq1auzZMkSfvzxR5M+kBAi8ynObHihh56/qqE7hx1s6R7vZt2VNyR0CyGEEKnw+uJhCjCbIHREmD90R8U+PN/2L+TMBgE9rBa6AX7mBO7MwY+TqX6vG040oZSEbmE1JvV4jx49mnXr1nH27Fny5MmT4PXw8HDefPNNunbtyowZM3j8+DEVK1akRIkS/Pnnn2Yp3FrkCbkQ5qdhAsToofcmWP0/yG4Hm7tB63KGNhUpwDmG2K5IISxI7i3WJz9zkVXs4yrNWGb5C706Yi0udLcobfnrGmGHhmsMkxAtrM7iPd7+/v507NjRaOgGyJMnDx07dmT16tUA5MuXj3feeYdz586ZcjkhRCYRgu5l6P5gqxq6s2lhQ5d4obsZpSR0C2EGL168YPbs2dSuXRtHR0eyZXu5tMupU6f45JNP+Pfff21YoRDC3OIWD7Oo6BjovjH+iDULhm4PirGYdtgl8rkS2zpMiPTEpMXV7t+/T3R0dJJtXrx4wb179wzfFylShJiYGFMuJ4TIBPw4SX+2gl6Bgdtg+Wmw08DazuD1couwcTRiEskv4CiESNrz589p1aoVR44cwcXFBUdHR8LDww2vlypViiVLluDs7MyUKVNsWKkQwpx2chnFkoupxU0T+/Vc7Ii17tCqjOWuB2ygK2444UkZgrhJNzbE+4SyKrnICEzq8S5Tpgzr16/n4cOHRl9/+PAh69ato0yZl38Jb926hbOz/IUQIisKQaeGbkWBQb+B31+g1YB/J+hYydCuMLkldAthJtOmTePw4cNMnz6dO3fu0L9//3ivOzk50bhxY3bu3GmjCoUQ5hY3v9tisfuFHnr9Cuv/UUP3pm7wTllLXQ0NsJh2hiHkbjjRhSp8Q0tDr76sSi4yCpOC95AhQ7hz5w7Vq1dn3rx5nDhxgps3b3LixAnmzZtH9erVuXv3LkOGqENF9Xo9e/fupVatWmYtXgiRMRRnthq6h+6An06od9Ll70HXyoY2buTlNiNtVqMQmc3atWtp2rQpI0eORKPRoNEkHKJZunRpbty4YYPqhBApFYKOfVwlBF2ybS8Sarl9u1/ooc8mWPs32GthY1d4t1zy7zORBjhGf3yoHu+4HycZFbs6uwaYTosEbYRIj0waaj5w4ED+++8/pk+fzvDhw+O9pigKWq2WMWPGMHDgQABCQ0P57LPPqFevXtorFkJkGIYtwxQFPtsF3/+p3kl/aQ89qxravUlBzjDIdoUKkQnduHGDDh06JNkmb9686HTJ/zIvhLCNV7cF06JhIV74UJ0QdFwklHKvbYVVzlLDrWP00G/zy7VZ1neBtuUtc61YChBO/KmtxlZsH8MeelBFerxFumdS8AaYNGkSvXv3xt/fnzNnzhAWFoajoyNvvfUW3bt3p3z5l38ZXVxcGDp0qFkKFkJkDIY53YoCXwTCrKPqCz+3hX7VDO1K4CihWwgLyJs3b7y1Voy5fPkyBQsWtFJFQojUeD1k6lEYSAChPDfsx/1qGLeYGD14b4FVZ9XQva4ztK9ouevFipu3/epDBmM9+nELq0nwFumdScF7+fLluLq64unpyVdffWXumoQQGZxhTjfAhP0w47D69fx34cMahnZFyMN1fK1foBBZQJ06dQgICODx48fky5cvwes3b95k+/btyfaKCyFsI7GQOYrdhqN6FAawlbxkpxT5WMVZ8xahV6B/AKw4oy6IuqYTdKiU/PvSKG7e9k4ux+vx/5oWaNHE+7nIwmoiozBpjrePjw87duwwdy1CiEyiOLPVL6YchEkH1a/neMInL9d5cCQ7t/jcBtUJkTV8/vnnPHr0iObNm3P48GFevHgBwLNnzwgMDMTT05MXL17g6ysPv4RIj4xtC6ZFk2AGtx7oxgZqs5i5HDNfAXoFPtwKS0+poXt1J+j0hvnOnwgNEIQPnpRJ0OM/mj3MoIVhWzFZWE1kJCYF7yJFihhu4BnF/PnzKVmyJDlz5sTDw4M///zT1iUJkSlpmKB+8fUhGLdP/XpmSxhax9DGAXt0fGH94oTIQho1asQPP/zA2bNnadSoEdOmTQPUIeitWrXi0qVLLFiwgBo1aiRzJiGELbjhxEK84oXMr2lh6R26VXoFPtoGv5xSdyFZ1RG6VE72beYQN7c7sR7/WhTlGsPYR1+uMUwWVhMZhklDzdu1a8fu3buJjIwkR44c5q7J7NauXYuvry8//fQTHh4ezJkzB09PTy5cuEChQoVsXZ4QmYYhdM8OgtGB6tfTmsFnLxdWzE02nvCl9YsTIgv6+OOPadKkCT/99BPHjh0jNDQUR0dHPDw8+OSTT6hc2Tq/SAshTONDdTwpwyVCDcOpR7LbshfVK/DJb7DopBq6V3SAblUse81XvDp0PLFh5W44SS+3yHA0iqKkes8BnU5Hs2bNKFy4MN988026v3F7eHhQq1YtfvjhB0Dd3qx48eIMGTKE0aNHJ/nesLAwnJyc0Ol0ODo6WqNcITIkQ+j+4U8Y8rv69YTG8FUTQ5uc2PGccVavTYj0Ru4t1ic/c5EZ7OMqzVhmuQsoCgzaDj8ej936swP0qprs28zpG1ryOfUBdaHWgQQQg2IYVi493CI9Sc29xaQe77fffpvIyEhOnTrFjh07yJkzJ4UKFUqwR6hGo+Hy5cumXMJsoqKiOHHiBGPGjDEc02q1tGjRgqCgIBtWJkTmkZvJ6hc/H38Zur9oAOMbG9qUIR+XGGb94oQQQohMIm7et0X26lYU+PT3l6F76XtWD90AtShq+Pr1Hn/p5RYZmUlzvPV6PdmzZ6dEiRKUKFHCMFxbUZR4f/R6vVmLNcWDBw+IiYnB1dU13nFXV1fu3LmToH1kZCRhYWHx/gghEmfHBJ4TA7/8BR/9ph78vB5MaQaxD+MKkktCtxBWdvjwYXx9fY3e6wBu376Nr68vR48etXJlQghTvT7v26Rf5I1RFBi+E34IVkP3L+2hz1vmOnuKaSHBCuVuONGEUhK6RYZnUo/3tWvXzFxG+jF9+nQmTpxo6zKEyBDsmIAeYMVp6B+7fdhQD/i6hSF0Z0PDPUbZrEYhsqpZs2Zx5swZZs2aZfT1IkWKsG3bNv777z/Wrl1r5eqEEKbyoTpXeMQ0/sAsXVyKAiN2wdzYFdEXt4N+1cxx5mRpwNB3rwEW0k4Ctsi0TAreGYmLiwt2dnbcvXs33vG7d+9SuHDhBO3HjBkTb2uVsLAwihcvbvE6hchotExQb5Zr/gf9tqh3zk9qwmxPQ+h2wF4WUhPCRoKDg2nevHmSbRo1asTu3RZeqEkIkWoh6LhIKOVeG14dgo4pHORnTpjnQooCI3fD7NiRLwvbwgdvm+fcyYjbl9s99vPVpbiEbpGppTl4//PPP5w/f57w8HB69+5tjprMKnv27NSoUYPAwEDee+89QB0qHxgYyODBgxO0z5EjR4ZYqV0IW8rBRDV0b/gHev0au9dndfj+XUPodsOBm3xm0zqFyMru3btHsWLFkmxTuHBh7t27Z6WKhBAp4cdJw/7VWjQsjF1QzI+TfMhW883uVhQYvQe+jV3z6Mc28KH1theM25f7GsMkcIssweSpIcHBwVSrVo0333yTLl260K9fP8NrBw8eJHfu3GzdutUcNaaZr68vixYtYtmyZZw7d46PP/6Y8PBwvL29bV2aEBlONiYShQJbzkOPjRCjqEPSfmqrbjuCOrxcQrcQtpUvXz5u3LiRZJvr16/j4OBgpYqEEMkJQWcI3aCG0wFsZRS76G/u0P3lXvjmiPr9/Hfho5rmOnuKxaBwiVCrX1cIWzApeP/99980a9aMq1evMnz4cFq3bh3v9YYNG+Li4sL69evNUmRadevWjW+//Zbx48dTrVo1w2rsry+4JoRImgNTiEGB3/6FLuvhhR56vgmLvQyhOzsaovnKxpUKIerUqcOmTZu4efOm0ddv3LjB5s2bqVevnpUrE0Ik5iKhCVYs1wPfcMR8F1EUGL8Pph9Sv/++NXxSy3znT4VX9+wWIrMzKXh/9ZX6S/WJEyf49ttvqVUr/l9WjUZD3bp1CQ4OTnuFZjJ48GCuX79OZGQkx44dw8PDw9YlCZGhFOEbwnkBOy9Bx3UQrYduldXtRuxe/lNyWVYvFyJd8PX15dmzZ9SvX5/ly5dz+/ZtQF3NfNmyZdSvX5/nz58zYsQIG1cqRNYRgo59XCUEndHX47YLs6iJB2DKH+rXczxhcG3LXi8RcftyyzBzkVWYNMf7wIEDdOrUibJlyybapkSJEuzYscPkwoQQ6YcLM3hIBARegffWQlQMdKwEKzpAtpehe7GsRipEutGoUSNmzZrFiBEjDFOrNBoNiqL2pmm1WubOnUujRo1sWaYQWUZic7dfX0htIV4MJEAdYWZukw6owRtgVisYWsf810ghi+xFLkQ6ZlLwfvLkiWHv7sQ8f/6cmJgYk4oSQqQfOZhEFHo4cA28VkPEC/AqD6s7gb1dbBstlxgqoVuIdGbo0KE0bdqUn376ieDgYHQ6Hfny5aN27dp89NFHVKlSxdYlCpElGJu7PZAAHvGcUexJEMYf8ZzPMfOOA1MPwlf71a9ntoThdc17/lRSgIEE4EkZ+f1BZAkmBe/ixYtz9uzZJNucPHmSMmXKmFSUECJ9cONbNXQfvgFt/OH5C2hdFtZ3gexq6LYDIhhv20KFEImqWrUqCxYssHUZQmRpxuZux6AYQje8DONVKcRIc4fuGYdg7L7Yr5vDZ5ZZ20EDvIkrZ7ibbFt4ubiaBG+RFZg0x7tt27bs2rWLPXv2GH193bp1HD161LB9lxAi4ynCTP7jKRwLgdarIDwaWpaGX7tBDvWZnQP2vGCCbQsVQhhlZ2dHz549bV2GEALjc7e1JBxuHYOCP2fNOwj7m8MwJlD9elozGNXAnGc3qEJBjtGfVpRJ8Tx1WVxNZCUmBe8vvviCokWL8u677/Lhhx9y/PhxABYsWEDv3r15//33KVmyJL6+vmYtVghhHd1Yxx3C4cQt8FwJT6KgaUnY3B1yvgzdT/jStoUKIRLl6OhI8eLFbV2GEAIMc7ftYgOpHRpm0NJIGNcwh2Pmu/B3R2BUbEfZ5KYwpqH5zv2achTAg8V8y5FE529r0RjChyyuJrIak4aaFyxYkAMHDtC7d2/8/PwMxwcPHgyAh4cHq1evxslJ/iIJkdGMJZB1/AOn7kDLFaCLhIYlIKAH5LYHoB5uHKa/jSsVQiSldu3anD592tZlCCFi+VAdT8pwiVDKxi6k5kwuw0JqWjTmXXBsdhB8FjtkfUJjGGvZhRQ3cz7J6uOC9us/AyGyCpOCN0Dp0qU5fPgwp06d4ujRo4SGhuLo6IiHh0eC7cWEEBnDTA4zlT/gf/egxXJ4FAF13eC39yFPdgAKkVtCtxAZwIQJE2jSpAnLly+nT58+ti5HCIHa8/162Hx1jrfZzDsGvrvUr8c1gq+amO/cRmgg0eq1wBo6U5fihs8ugVtkRSYH7zjVqlWjWrVqZihFCGFLwYSoi7mcuw/Nl8PD51CrKPzeE/LmMLS7y0gbVimESKndu3fTpEkTvL29+f7776lVqxaurq5oNPGHtmo0GsaNG2ejKoXIuuJWOjf7plrz/4ShsVv6ftkQJjYx9xUSGExt5hOc4OFBXC93F2QHBSE0StyGnsKosLAwnJyc0Ol0ODo62rocISzCj5P0Zyv8+xAaL4U7T6F6EdjTG/LnMrS7yXB5Si2EGVjj3qLVpmwZF41GkyW2/5T7uUhv9nGVZiwz70l/DIZPtqtfj64P05qDJmULnaXFOjoTRtQrw+bBl3oMxUN+bxCZWmruLSb3eN+/f58lS5YQHBzM48ePjd60NRoNgYGBpl5CCGEFwYSooftyKDRbpobuqq6wq1e80P0lDeXmKUQGsm/fPluXIIRIQtxK52YbYr7wxMvQ/Xk9q4VugO5sZCFeXGOYzN8WIhEmBe8zZ87QrFkzHj16RFId5q8PZxNCpC+Gnu5rj6HZcvjvCbxRUO3pLpDb0K4zlZhCc9sVKoRItcaNG9u6BCHEa0LQcZFQysVuoVWW/PxLaNpPvPgkDNymfj2iLnzdwmqhG17uQX6NYTShlNWuK0RGYlLwHjFiBKGhoYwdOxYfHx/c3Nyws7Mzd21CCAsKQaeG7ps6aLoMbuigQgEI7AMF8wBQkFz8Rk9q4WbjaoUQQoiMKS5sB3KF6RxCj5LkYmSp9stf8GGA+vUwD5jZ0qqhO04MCpcIlZ5uIRJhUvAOCgrivffeY9KkSeauRwhhJRX4Hv4LU0P3tcdQ1hn29oXCDgB05Q3W0tW2RQoh0uzw4cMsXbqUU6dOERYWhqOjI2+//TZ9+vShQYMGti5PiEzNj5MMICDBcHKzhe5lp6D/VvXrT2vDLE+bhG5QF1IrG9uTL4RIyKTgnT17dsqUKWPuWoQQVuLAVJ7deawOL7/8CErlg719oGheANxxktAtRCYwfPhw5s2bZ5gWptFoUBSFEydO4Ofnx9ChQ5k1a5aNqxQic4pbtdys24S9asVp8N6ipvhBtWDOOzYL3drY1cult1uIxKVsydPXNG7cmOPHj5u7FiGEFdgzkfB7j9Utw/59CCWcYF9fKK7eLHNhxzWG27ZIIUSaLVu2jLlz51KuXDlWrVrFrVu3ePHiBbdv38bf35/y5cszd+5cli9fbutShciULhJqudC96gz0iw3dH9eE71tbJXRrYv+8SgscxQcfqlv8+kJkZCYF72+//Zb//e9/fPvtt+auRwhhQVom8OJBOLRYDv/cBzdHNXS75wPUm+kzZD9fITKDH3/8ETc3N44dO0aPHj0oXLgwGo0GV1dXunfvztGjRylWrBgLFiywdalCZCoh6NjHVTZxzjIXWH0W+mwGvQIDqsMP71qtp1sBRlAPu9j4bYeGhbSTtWCESAGThppPnTqVKlWqMGrUKH766SeqVatmdN8yjUaDn59fmosUQqRdbiajhD6Hlivg7D0o4qAOLy+dHwBHsqPjCxtXKYQwl7///pv+/fvj5GR86KeTkxOdOnVi8eLFVq5MiMwrsTndZrPub+i1SQ3d/d+GH9uC1rrDy13JI9uGCWECk4L30qVLDV9fuXKFK1euGG0nwVuI9MGJaTx/HA6eK+HUHSiUR11IrVwBAPJgL6FbiCxItv0UwnwsPqd7wz/w/kY1dHtXg5+9rB66AUazh+5UkW3DhEglk4L31atXzV2HEMJCsjGRmLAIaL0Kjt8Cl9xqT3dFF0Cdb/KUL21bpBDC7CpXrszGjRuZPHkyDg4OCV5/8uQJGzdupHLlyjaoTojMx6Jzun89B903QIwCfd+Cxe1sErpBtg0TwlQmBW93d3dz1yGEsIAizCTmaSS8uwqOhoBzLtjTGyoXMrS5LgupCZEpDRw4EB8fH+rWrcuECRNo3LgxLi4uPHjwgP379zNx4kRCQkJka1AhzKQczmjRmD98bz4P3WJDd++q4Ge70A2ybZgQpjIpeAsh0r8CzCD02RNo6w+Hb0K+nLC7N7xV2NDmSxrKE2shMilvb2/++usvfvjhB7p2VbcH1Gq16PV6ABRFYciQIfTt29eWZQqRabjhxEK8GEgAMeYK31svQNf18EIP778JS9qDnUlrI5uFnWwbJoTJUvw3t1mzZgm2HDl27Bjz5s0z2v7rr7+mQIECaatOCGGSoswk9PkTaL8GDlwHxxywqxdUL2Jo04xSTKG5DasUQljavHnzOHDgAP369aNatWqULFmSatWq4e3tzYEDB5g7d66tSxQi0whBR2ny8w0tzXPCbf9C53UQrYfuVWDZe1YN3Vo0zKQl++jLn/RnH325xjDZNkwIE6W4x3v//v00adIk3rEdO3YwadIkPv300wTtIyIiePz4cVrrE0KkUhFmcidCBx3XwZ4r4JAdfu8JtYoZ2jSlJIFIL5cQWUHDhg1p2LChrcsQIlMz+2rm2y9Cp9jQ3bUyrOgA2azb062gkJ9csoiaEGZiu7EqQgizK8TX3IkKgy7rYcclyG0P29+HesUNbbryBnvpZ7sihRAWExYWRlRUlK3LECJL2cYFPmSr+UL3zkvQcS1ExUCnSrDS+qEb1D27BxJACDqrX1uIzEiCtxCZhBeruB/9VF31dNu/kDMbbOsBDV8uhliTIqylqw2rFEJYUv78+fn666/jHUtqWpgQIm36sQkvVptvObVdl9VpYpEx0KEirO4E9nbmOnuqxa1gLoRIOwneQmQC9VnMthcXoOevsOk85LCDrd2h6cvhYd/RimAG2rBKIYSlKYqCosSPADt27GD4cNm9QIi0CkHHPq4aeoCDCWEZp813gT1XXobu9hVgTWebhm6QFcyFMCdZ1VyIDM6d2dyIeQR9N8P6fyC7HWzqBi3LGNq0oRy+1LNdkUIIIUQG9uocbi0aFuKFjgjzXWDvVWi3GiJegFd5WNdFvZ9b2WfUYzZBxKDICuZCmJkEbyEysJLM5ob+MfhsBf+z6hywDV2gdTlDmzcpxDZ62q5IIYQQIgMLQRdv4TQ9CgMJYDHtzHOB/dfUrT+fv4A25WC9bUK3FhiKB0Px4BKhlMVZQrcQZpSq4L1t2zbu3Llj+P748eMAfPLJJwnaBgcHp7E0IURSKjKP6/rHMDAAlp0GOw2s7QxeFQxtipGXMyT8+ymEEEKIlLlIaIKF02JQ+JVzaT/5wevQJjZ0ty4LG7tCDuv3i2mAhbQzBG0J3EKYX6r+Zh8/ftwQtl/1008/GW2v0WhMq0oIkaR8TEOnRMLg7bD4L9BqYFVH6FjJ0GYoHsyhtQ2rFEIIITK+cvy/vfsOj6J62zj+3QQSCKTQexMQUZAOgggklFAFRBBFmnTpIM1CsaIoCFiQIuUHr6CoaASE0BFCkY4CAlIMVSkJNSHJef9YsrBkU0myKffnuvaSnTk788zsmtlnz5nn5MYFi13y7YKFAP56uA3/dhqaL4Kbd8C/NPzwQqon3RbgC5rTknJKtkVSWIL/716/fn1KxiEiCZSNCYSZKBjyK3z5u/WqOb8NvFDB1qYE3kq6RTKphQsXsm3bNtvzY8eOAdC8eXOH7S0WC8uXL0+V2ETSq2HU5hO22lLvh546bOs/0GwR3LgDjR+x1mbJlrpJtwvWXu4eVE3V/YpkVgn+P7x+/fopGYeIJEB23rEm3SMCYdoO68KvW8PLT9raeJCFk6iCsUhmdezYMVuyfb9ff/3VYfu0MDqtZMmSnDp1ym7ZBx98wOjRo23P9+/fT//+/dm5cyf58uVj4MCBjBw5MrVDlUzm/qJqyWZbMDRdCNfDoWEpWNYRsmdNvu0n0E+8SEvKxd9QRJKFiquJpBPuTCDcRMHra+GTIOvCmS2hW2Vbm6Lk5B9ec06AIuJ0J06ccHYISfb222/Tq1cv23NPT0/bv0NDQ2nSpAmNGjVixowZHDhwgFdeeQUfHx969+7tjHAlE3iwqFqy2HEG/BfCtXDwLQk/vwgeqZ90A+TEzSn7FcmslHiLpAPZeYdwDEzYCBO3WBd+3hx6VbO1yUlWJd0imVyJEiWcHUKSeXp6UrBgQYfrFi1aRHh4OF9//TVubm488cQT7N27l8mTJyvxlhTjqKjaQ/n9LDT5H4SGQf0SEJA6SXdbHmMZh+2ORPNzi6Q+F2cHICJxK8On3CYS3ttkTbwBPvWHV2vY2pTGh2u84aQIRUQe3sSJE8mTJw9VqlRh0qRJRERE2NYFBQVRr1493Nzu9dD5+/tz5MgRrly54nB7YWFhhIaG2j1EEiMnyZgU7zoLjf8HIWHwTHH45SXIkfI9zi0oy08ciZF0a35ukdSnHm+RNKwYnxDMNfhoC7x5t8DhpMYw+ClbmxJ4cYwhzglQRCQZDBo0iKpVq5I7d262bt3KmDFjOHfuHJMnTwbg/PnzlCpVyu41BQoUsK3LlStXjG1+8MEHTJgwIeWDl3QvmBCOcpmyD8xbfZ07ybODPeesSffV2/B0MVjRCXKmfNJdl2Ks5NgD1dghiB7UoGiK719E7KnHWySNysNEa9I9JQhGrbEufM8PXqtja1OQHJxkmJMiFBGJ3ejRo7FYLHE+Dh8+DMCwYcNo0KABTz75JH379uWTTz5h+vTphIWFJXn/Y8aMISQkxPb4559/kuvQJIMIJoQRrKY4U/BjPiX4lDnstq2Pnkbsoew7D43+B1duQ+2isDJ1km6Al6gYY6h8FHAjuX5QEJFEUY+3SBr0NHO4zG34fAcMW21dOK4+vP6MrY0XbpxjhJMiFBGJ2/Dhw+nWrVucbR555BGHy2vVqkVERAQnT56kXLlyFCxYkAsXLti1iX4e233h7u7uuLu7Jz5wyRQcVSuPwtCbn3mS/NSgKEXxpjNPMp99SdvJ/gvQcAFcvgW1isCvL4Nn6nwmu1KJVpRjACvtjlH3dos4jxJvkTRmCCvZyj8wcxcMWGld+Hpda+J9V06yEsLrTopQRCR++fLlI1++fEl67d69e3FxcSF//vwA1K5dmzfeeIM7d+6QNav1vtvAwEDKlSvncJi5SFziqlYeBTzFHGbSCn9KsyCpSffBi9ak+9ItqFkEVr0MXimbdLtgoT2PM5zatqHkM2lFHwKIxOjebhEnU+Itkobk4gOuEgZz90CfX6wLX6sN7/rB3bl2PciiQmoikmEEBQWxfft2fH198fT0JCgoiKFDh/Lyyy/bkuqXXnqJCRMm0KNHD0aNGsXBgweZOnUqU6ZMcXL0kh7FV608CkMfAphOs6TVNP/jIvjNh/9uQvXC1qTbO1uS402IVjzKW9SjEJ4c5TLBhFAUb3pQFX9Kc4zLlHngHnYRSV0JSrxjGwoWH4vFwvHjx5P0WpHMpjSfWpPuhfuhx8/WhYNrwUeNbUm3K3CDN50XpIhIMnN3d2fx4sWMHz+esLAwSpUqxdChQxk27F79Cm9vb1avXk3//v2pVq0aefPmZezYsZpKTBLkweJp0fdux5V8R2LYxdnE7+zQv+C3AP69CVULweqXwSdlk26AAP7iF/4CwGDt/Z5JK3pQlaJ4K+EWSQMsxph4f8wrWbIkFot9cYnw8HDOnTsHQJYsWciTJw+XLl2yTf9RqFAh3NzcOHHiRAqEnXpCQ0Px9vYmJCQELy8vZ4cjGVRJpnCKEFh8EDr9AFEG+lW3ztV93/97hvHOC1JEko2uLalP5zxzuv9e7vuT0Tnstg3BdsGCwTz8jN2H/4MG8+DCDahcENZ2gdzZk+EoksYVCycZoqRbJAUl5tqSoB7vkydP2j2/evUqjRo1omzZsrz33nvUrl0bFxcXoqKi2Lp1K2+++SY3btxgzZo1ST4IkcwiDxOthdS+/xNevpt096wCnynpFpGEe+WVV5L0OovFwpw5c5I5GhHne/Be7ugh5P6UpgdVeZL8/MZp6lKcCWxkOUeTvrO/LlmHl1+4AZUKwJrOTk26wdprf4zLSrxF0ogk3eM9atQobt++zfbt23F1dbUtd3FxoW7duqxZs4ZKlSoxatQoZsyYkWzBimQ02XibMKLgp8PQ8XuINNCtMnzVClysSbcFiFLSLSLxmDdvnsPlFosFR4Pbopcr8ZaMytG93NHJ6CqO25JyCzxcb/fRS+A7H85dh4r5YU0XyOPxMFtMFqpgLpK2JGke759++omWLVvaJd33y5IlCy1btuSnn356qOBEMrK6zLYm3SuOQvvvICIKOlWE2feSboDTDHVilCKSXpw4ccLucfz4cVq2bEmePHl455132LBhA4cOHWLDhg28/fbb5MmTh1atWnH06EP08omkYY7m4XbFQg6y2vWEP1TSfeyyNek+ew0q5LcOL8+bekm35e4j+t/Rx6sK5iJpT5J6vENDQwkJCYmzTUhISLxtRDIrP+axhWBYfRyeWwJ3oqDDEzCvDbje+z1sNs/qoikiCVKiRAm75xMnTmT79u3s27ePQoUK2ZaXK1eOevXq0b17d6pUqcLSpUsZOXJkaocrkuKK4h1jOq2h1OYEV+MsrJZgf1+xJt1nrsHj+axJd74cD7/dBIg+lsHUArBVLb//3/r+IJK2JKi42oNq1KjBsWPH+P333yldunSM9UePHqVGjRo8+uij7NixI1kCdRYVY5Hk1pJF1vvI1p2AFv8HtyPgufKwuB1ktY4icQFOMVQXTZEMKjWuLWXLlqVp06ZMnz491jYDBgxg1apVmaLXW9fzzCuYEKaynU/YmhzpttXJq1B/HpwOgcfywoauUCBncm09ThbgZ16kJeVSZX8iErtkL672oDfffJO2bdtSpUoVevToQd26dcmfPz8XL15k8+bNfP3119y4cYM339S0RyL3m8xWa9K96RS0+saadLd6FL65l3S7YSGMcU6OVETSu+DgYLJli3sao2zZshEcHJxKEYk4z2SCki/pPnXVWr38dAg8mgfWdUm1pBusQ+Nbs9hWod2RB6dQExHnS1KPN8CCBQsYOHAg165ds5tqzBiDl5cX06ZNo0uXLskWqLPoF3JJLk8zm60Ew9Z/oMn/4MYdaFoGlr0A7vd+A1P1cpGML7V6vI0xHDx40GECfvPmTSpWrIiLi4t6vCXDcJRwrucEfsxPnh2cDrEm3SeuQtncsKEbFPZMnm0nUmzThcU2hZqIJL/EXFuSVFwNoEuXLgQHBzN//nyGDh3KK6+8wtChQ5k/fz6nT5/OEEm3SHKZzFZr0r3jDDRdaE26Gz0CP3RQ0i0iKaJnz578/fffPP300/z0009cunQJgEuXLrFs2TLq1q3LyZMn6dWrl5MjFUkeH7OFEnyKH/MpwafMYTfguMhakvwTYr2n+8RVKJMb1nd1WtIN1grtU9lutyy2KdSCUd0lEWdLco93ZqFfyOVhNWQ+6zgBu89Z5/gMCYMGJWH5S+CR1dZOSbdI5pEa15aoqCh69erF3LlzbSPTXFxciIqKAqwj1Lp3787s2bPtRq5lVLqeZ2yT2MJIAu2W3d8j3IFv+Y4/k76DM6HWe7qPX4FHcsHGblDU+Z+jB2vCxNa7v56uNKBUKkcnkvGl+D3e97t+/Tp//fUXN27c4JlnnnnYzYlkKOWZzmEuwb7z0GiBNemuWxwCXlTSLSIpysXFhTlz5tClSxfmz5/P/v37CQkJwdvbm0qVKtG5c2caNGjg7DBFHlowIYx6IOmGe3N2n+MaSx8m6T57zdrTffwKlPKx9nSnQNJtAbpQiYXsJ/LuMPH4qq9HAUH8Q/u7iXd07/79r9N83iJpQ5IT75MnTzJ48GBWrFhBVFQUFouFiIgIALZs2UKvXr344osvdFGXTOsVllmT7oMXodH/4MptqF0UVrwEOd1s7ZR0i0hKql+/PvXr13d2GCIp5iiXHaanLljYydkYPeGJcu5u0n30MpT0sSbdxVOuWFlzytKfGtzgDjnIylPMiTf5foGlhBJOD6o6nEJN83mLpA1JSrxPnz7NU089xaVLl2jdujXnz58nKCjItr5WrVr8999/fPPNN0q8JVP6mC3MZS8c+hcaLoD/bkKNwrCyE3i629op6RYREXk4jnp5AcZQ12FPeIKdvw5+C+CvS9Zke31XKOHzcMHGwWBNoqMLoj1CrgTNN26APgTgT2mK4k0PquJPac3nLZLGJKm42rhx47hy5QobN25k6dKlNG7c2G59lixZeOaZZ9iyZUuyBCmSngQTwggCrRdqvwVw8QZUKQirXgbve5WFlXSLSEqLiIhgypQp1KxZEy8vL7Jkufd7+969e3n11Vf566+/nBihyMOL7uV1vVtAzQXoQzWepEDSpxC7cN36w/nh/6CYlzXpLumTTBHHLbogWk6yJrgoXPSw+mhF8aYBpZR0i6QhSUq8V61aRdu2balTp06sbUqUKMGZM2eSHJhIelWMKXD8srWQ2vnrUDE/BHaGXNltbZR0i0hKu3XrFr6+vrz22mucOnUKLy8v7q+nWqpUKebOncuCBQucGKVI8uhBVU4yhNeogwG+YhcvsDRpG7t4w5p0//kvFPG0Jt2P5ErWeOMTieEGd+x+UHDFwmvUYQc9YyTkuo9bJO1LUuJ9+fJlSpYsGWcbYwxhYWFJ2bxIuvQLR7AwHk5etfZ0n7kGj+eDNV0gj4etnZJuEUkN77//Plu2bOGDDz7g/Pnz9OzZ0269t7c39evXZ9WqVU6KUCT5TSYo6b3cYL01rNEC+ONf61RhG7pB6dRPaKMT6egfFNbTlZMMYRJNqEHRGAm57uMWSfuSdI93gQIFOHr0aJxtDhw4QPHixZMUlEh6U50Z7OK8dY5Pv/lwOgTK5YG1XSB/Dls7Jd0iklqWLFmCr68vI0eOBHA4ZdgjjzzCnj17Ujs0kRRxlMsJuic6VpduWnu6D1yEQjmtPd1lnJN0359IF8U7RlKt+7hF0p8k9Xg3btyYX375hf379ztcv3nzZtatW0fz5s0fKjiR9OAVllmT7jOh1p7uE1etF+p1XaFgTls7Jd0ikppOnz5N9erV42zj6elJSEhIKkUkkrLKkjuBd0Q7cPmWdQaS/Res1+51XeHRPMkZXoJ9Qzt6UDXedrqPWyR9SVLi/eabb5I9e3bq1avHe++9x7FjxwBYuXIlb731Fk2bNiVv3ryMGDEiWYMVSWsmRVcvP3+3CMuxy9Y5Ptd1sQ5Ru0tJt4ikNk9PTy5evBhnm+PHj5MvX75UikgkZRXFmw9pHH/DB125ZR1evvc8FMhhvYY/ljf5A0wAVyzUpphT9i0iKStJQ81LlizJqlWr6NixI2+99RYWiwVjDC1btsQYQ/HixVm6dCmFChVK7nhF0oxgQqxzg0YXYTlyd7qRdV2h2L1fn5V0i4gzPPXUUwQEBHD16lV8fHxirP/nn39YsWIFbdu2Tf3gRJJBMCEc5TJl7xtqPYKnsQAjCUzYoPOrt6Hx/2DPecjnYb2Gl0+dH6NcgKj7nutebZGMLUmJN1jn6j569CgBAQFs376dy5cv4+XlRa1atWjdujVubm7JGadImlOMKdb7wRrdV/l0XRe76Ub+YajzAhSRTG3EiBH4+vrSsGFDpk2bRkREBAA3b94kKCiIgQMHEhERwbBhw5wcqUjizWE3vQkgCoMF+JDGjOBpADpSwTqtZ3xCbkOT/8Guc5D3btL9eOok3a5YCKIHN7hDDrJygzu6V1skg0ty4g3W+brbtm2rX8slUwkmxJp0X7ll/ZX8/iIs91U+fZ1ndAEVEaepV68en332GYMHD6ZevXq25Z6e1ttgXF1d+eKLL6hWrZqzQhRJkp0E04ufbT3aBmsP93Eu8yb1eJ5v499IaBj4L4SdZyFPdusP5xXyp2TYNi53e7ZrUDRV9iciaUOS7vH28/OLd97PhQsX4ufnl6SgRNKqOey2Jt0ht60X7D3nrVXL13WFsveKsLSgLO/R0ImRiohAv3792LdvHwMGDKBGjRqULl2aKlWq0LdvX/bs2RNjijGRtG4Ou3mKOQ6HkX/FLooxhe2ciXsj18Kg6ULYfgZyZ7fOQFKxQIrE65jBn9KpuD8RSQuS1OO9YcMGGjRoEGebU6dOsXHjxqRsXiRNCiaEnvx894K9yPoreV4P6wX7viIsHXicJXRwYqQiItaq5j4+PpQvX56pU6c6bHPt2jWuXLmi6T8lXQgmxDa8PMmuhUGzRRAUDLmywZrOUKlg8gWZAFHAMS5rVJxIJpOkHu+EuHHjBlmzZk2pzYukumJMgevh0Pz/YFuw9VfyNZ3thqZVJL+SbhFJE0qVKhVrwh1t2rRplCpVKpUiEnk4Dz1P9/VwaPF/sOUf8MkGgZ2hSuoXAnbFQhlSf35wEXGuBPd4nz592u751atXYywDiIyM5J9//uH777+nZMmSDx2gSFpgYTzcvAOtvoHfTt+7YN/3K/lAajINzV0vImmDMQZj4k5S4lsvkpaUfZhk9cbdpHvzafB2t17DqxVOvuDiED23uEGVy0UyswQn3iVLlsRisf7psFgsTJ06Nc5f0o0xTJo06eEjFHEyC+Ph1h1ovRg2nARPN1j1MlS99yv5DnqqSIqIpDvBwcG2YmsizuBoSrDYDGNV0nYS/cP5plPg5Q6rO0P1lE+6XYBh1GEwtQDr8HJVLhfJvBKceHfp0sU2X/eCBQuoVKkSlStXjtHO1dWV3Llz4+fnR9OmTZMzVpFUZ2E8hEXAc9/Cmr8hpxv8+jLULGJrM4CaSrpFJE14++237Z5v2LDBYbvo0WmLFy/mqaeeSoXIRGK6f0owFyzMpBU9qOqw7U6C+Y4/E7+T6KR7/cl7P5zfdw1Pbi5YmEgjalA4RpKthFskc7OYJIwzK1WqFEOHDmXQoEEpEVOaEhoaire3NyEhIXh5eTk7HElFFsZDeCQ8/y0E/AUeWWFlJ6hXwtamNLk4xmDnBSki6VJKXVtcXO6Vbon+sTwuhQsX5scff6RGjRrJFkNapet52hJMCCX41O6ebVcsnGSIXYIa3SP+E4eZyvbE7eTWHXh28b0fzle9DHWKJdchOOQCnGKokmyRTCIx15YkVTU/ceJEkgITSS8sjIc7kdBxqTXpzpYFAl60S7pL4K2kW0TSlPXr1wPW2738/Pzo1q0bXbt2jdEuenTaY489Zpesi6QWR4XSIjF21b7v7xFPtNsR0GaJNenOcfeH8xROukEVy0UkdklKvP/880/WrFnDiy++SL58+WKsv3jxIosXL6Zx48aUL1/+oYMUSU0WxkNEFLz8I/x4GNxd4aeO4Hev8m958vAnA50XpIiIA/Xr17f9e9y4cTRo0MBumUhaUZbcuGCJ0eMdXe37oaYOC4uA55bA6uPW0WorOkHd1JkyTxXLRSQ2SfqZe+LEiXz44YfkyZPH4fo8efIwadIkPvroo4cKLi4nT56kR48elCpViuzZs1O6dGnGjRtHeHi4Xbv9+/fzzDPPkC1bNooVK5aiMUn6Z2E8REZB12Xw7R+Q1QV+eAGalLa1yU02Jd0ikuaNGzdOSbekWUXxZiatcL1b8/vBat9b+SfpSXe7b2HlsbtJ90t2o9VSkiqWi0hcktTjvXnzZho2bBjr8DRXV1caNmzIpk2bHiq4uBw+fJioqCi++uorypQpw8GDB+nVqxc3btzg448/Bqxj7ps0aUKjRo2YMWMGBw4c4JVXXsHHx4fevXunWGySPlkYD1EGevwM/3cAsrjA0g7QvKytTSFycJYRzgtSRCSBtm7dytKlSxk5ciQFCxaMsf7cuXNMmjSJDh06qMCaOEUPquJP6RjVvqOHmCdaeCS0/w6WH4XsWeCXF6F+yeQNOhZT8Od5HlfSLSKxSlLiff78eYoVi/s+mSJFinDu3LkkBZUQTZs2taua/sgjj3DkyBG+/PJLW+K9aNEiwsPD+frrr3Fzc+OJJ55g7969TJ48WYm32MnKBGvS3ScA5u8DVwssbgfPlrO1KYE3JxnqxChFRBLuk08+Yf/+/UyePNnh+kKFCvHLL79w5swZlixZksrRiVgVxTtGMbUkDTEPj4QO39nXZfEtFf/rkoErFiXdIhKvJA01z5EjBxcvXoyzzcWLF8mWLVuSgkqqkJAQcue+d19NUFAQ9erVw83NzbbM39+fI0eOcOXKlVSNTdIuC+OJMFEwYAXM3gMuFlj0HLR73NamGgWVdItIurJz507q1q0bZ5t69eqxbdu2VIpIJH6Oiq7FK7oY6k9HrHVZfu4IDR9JmQAdmEgjJd0iEq8kJd5Vq1Zl2bJlXL161eH6K1eu8OOPP1K1quO5GFPCsWPHmD59On369LEtO3/+PAUKFLBrF/38/PnzDrcTFhZGaGio3UMyLgvjwRgYugq+/B0swPw28EIFW5s6FON3+jorRBGRJLl48SJFisQ9X3HBggXj/SFdJDUEE8J6TpCTrLjcve87Qe5Ewovf2xdDbVw6/tclo+oUTtX9iUj6lKTEu3///ly6dAlfX98Y93Fv3LgRX19frly5woABAxK97dGjR2OxWOJ8HD582O41Z86coWnTprRv355evXol5ZBsPvjgA7y9vW2P+IbUS/qVlQnWpHtkIEy9OzfonGfh5SdtbR4jD1vo4aQIRUSSzsfHh9OnT8fZ5tSpU+TMmTOVIhJxbA67KcGn+DGfp5hDM8ok7IURUdDpB/j+ELi5wo8vgH8CX5sElruP+6mKuYgkVJLu8W7dujVDhw5lypQp+Pr64u7uTsGCBTl//jxhYWEYYxgxYgRt2rRJ9LaHDx9Ot27d4mzzyCP3hg+dPXsWX19f6tSpw8yZM+3aFSxYkAsXLtgti37uqNAMwJgxYxg2bJjteWhoqJLvDMiT96zDy99YBx8HWRd+1RK6V7G1yYYrh1S9XETSqaeeeooff/yRf/75x+F17PTp0yxbtgw/Pz8nRCdi7eXeyj9293RHYVjO0fhfHBEFnX+E7/68OwNJB2hWNv7XPaQPacwY1hCJURVzEUmUJCXeYC3a4uvryxdffMHOnTsJDg7Gx8cHPz8/+vfvT7NmzZK03Xz58jmcG9yRM2fO4OvrS7Vq1Zg7d26MKuu1a9fmjTfe4M6dO2TNmhWAwMBAypUrR65cuRxu093dHXd39yTFLunDWNZxnTswYSN88Jt14WfNoHc1Wxsv3AjhdSdFKCLy8IYNG0ZAQABPP/007777Lo0bN6ZQoUKcO3eO1atX8+abb3Lr1i2GDx/u7FAlE4quXJ6kKcOip/1cfNCadH/fAVo8muwxPsgANSjMSYbEqMQuIhIfizEmCX/xnO/MmTM0aNCAEiVKMH/+fFxdXW3ronuzQ0JCKFeuHE2aNGHUqFEcPHiQV155hSlTpiS4qnloaCje3t6EhITg5eWVIsciqWcQK5jODnhvE7y53rpwij8MuTeVTnUKsZM+sWxBROThpda1ZerUqQwfPpzoS73FYrH928XFhU8//ZT+/fun2P7TEl3P045gQijBp0lPurv9BAv33532sz20fiz5g3TAFQsnGaJkW0RsEnNtSXKPt7MFBgZy7Ngxjh07RtGiRe3WRX+p8Pb2ZvXq1fTv359q1aqRN29exo4dq6nEMqnqfMUuzsGkLfeS7o8a2SXd3anM17RxToAiIsls8ODB+Pr6MmPGDHbu3ElISAg+Pj7UrFmTvn37UqFChfg3IpLMklS5HKxJd4+frUm3qwWWPJ+iSbf1nm4LURpWLiLJIN32eKcW/UKeMUxmK8NZDZ9us1YwB3jXF96oZ2vjS0nW0c05AYpIpqJrS+rTOU87ktTjHWWg588wd6816V78PDz/eLwvexgBvEhlCmpYuYjEKjHXlgRVNXdxcSFLliz89ddftueurq7xPrJkSbcd6pKBdONHa9L9+Y57Sfe4+nZJdx2KKekWERFJBUXxZiatcE3otGFRBnoH3Eu6/69diifdAP9jP0XxpgGllHSLyENLUGZcr149LBYLHh4eds9F0rpfOMJ89sHMXTBgpXXhmLrWxPuu2hTVlGEiku698sorWCwW3n//fQoUKMArr7ySoNdZLBbmzJmTwtGJ2OtBVZ4kPzWZHXfDKAP9foE5e8DFAgufgw5PJGpfdSlGQXIQhWEZR4hK4OuW8gfBNFHSLSLJQkPN46GhaenXHHbTk59h7h545Wfrwtdqw0eN4e4PR0+Sn3286sQoRSQzSolri4uLCxaLhUOHDvHoo4/GmOkjNhaLhcjIyGSJIS3T9TztWc8J/JgfewNj4NXlMGOXNen+X1t4qWKS9uWChZm0ogA5aMU3iYixKw0olaR9ikjGlymKq4nEZSfB1qR74X5rIRaAQTXtku6C5FDSLSIZxokTJwAoUqSI3XORtCiYEA7zb+wNjIEBK6xJtwWY3ybJSTdY5wfvxc9YHhje7oKFYdSmA49Ti9l2d527YqEMuZO8TxGR+ynxlgzH1tO95KB1nk8D9KsOnza1Jd3teZxv6eDUOEVEklOJEiXifC6SVsQ7h7cxMPhX+OJ3a9I9tzW8/ORD79cA5r59umBhGz2ogXV2nFk8Sx8CiFQVcxFJAQlKvBN6n9iDdN+YpLZgQqxJ9w+HoNMPd6ugVoHPmtuS7oHUZBrNnRypiIhI5hFMCEe5TE6yxp90D1sF03dYn89+FrpWTpGYojDc4I7teQ+q4k9pVTEXkRSRoMR73rx5DpdbLBYc3SIevVyJt6SmYEKYyS74+Qi8sBQiDXStBF+1st4bBjzP40q6RUREUtH9PdwWiH0SMWPgtdXw6Xbr81mt4JUqCdpHR55gMX8kKi5HQ8mL4q2EW0RSRIIS7wfvE4uKimLw4MFs27aNwYMH88wzz1CgQAEuXLjApk2bmDZtGrVr12bKlCkpErTIg2wX9RV/wfPfQkSU9V6wOc/aku4AXqQl5ZwcqYhIytDoNEmLggmx6+GOM+ketQYmb7M+/6ol9Kya4P0kJenWUHIRSU0JSrwfvE9s4sSJbN++nX379lGoUCHb8nLlylGvXj26d+9OlSpVWLp0KSNHjkzeiEUeYLuorz4Gzy2BO1HQ/nFrIRZXa1Xfj2ispFtEMjSNTpO06CiXHQ4rj+75dsVCpImC19fCpK3WlV80h97VUiQeF2AYdRhMLSXdIpKqEjbXyAPmzJlDhw4d7JLu+xUpUoQOHTowa9ashwpOJCGOcpmodX9D68UQFgltH4NFz0EWF1ywMInGjOBpZ4cpIpKiTpw4Yfc4fvw4LVu2JE+ePLzzzjts2LCBQ4cOsWHDBt5++23y5MlDq1atOHr0qLNDlwysLLlxeaCSOIAFC69Rhx/NC5R/8wBM3GJd8Vkz6FfjbpvYuQBjqZeoWFyAbfRkkubmFhEnSFJV8+DgYLJlyxZnm2zZshEcHJykoEQS4/Lmw9DqG7gdAS0fhcXP45o1C9/QjtoU08VVRDIFjU6TtCC6iFrZu8XJiuLNTFrRi5/t+r2jMHzCVj4e9x68vwmAMtNeYnj/weTBg9oU4xzXqMlsh/vpRTV6UY132Rxnj/r9osCumFpijkNE5GElqce7aNGi/Pjjj9y+fdvh+ps3b/Ljjz9StGjRhwpOJD5bt26lW/MOcPMONC0DS9vj6paFr2hFeyroYikimVZaGZ323nvvUadOHTw8PPDx8XHY5vTp07Ro0QIPDw/y58/PiBEjiIiIsGuzYcMGqlatiru7O2XKlIl1aL04RzAhjGA1JfgUP+ZTgk+ZxBbWc4Ir3HL4GjNhA7xjTbqZ4s+xgWWpRmHb9fsEV2Pd3yx2s4rjfEijGOtcsPAzL8boaU/ovNxz2G13HHPYHe9rRETik6TEu2fPnvz99988/fTT/PTTT1y6dAmAS5cusWzZMurWrcvJkyfp1atXsgYrcr8dO3bQrFkzrl+/TqNGjfjrh82sd+/BSYbQg4QXZBERyYjSyui08PBw2rdvT79+/Ryuj4yMpEWLFoSHh7N161bmz5/PvHnzGDt2rK3NiRMnaNGiBb6+vuzdu5chQ4bQs2dPVq1alaKxS8LMYTfFmcLHbLX1PkdhGEkgfsxnBIEx+6Tf2QjjN1r//XFjGPIUAN9w0NbkEjdj3WcUhj4E0JEKfERj2xdaVyzMpBUtKcdMWuF6N/lOaDG1B4vBRe8nmJCEnQwRkVhYjKOKK/GIioqiV69ezJ07F8vduZFdXFyIiooCwBhD9+7dmT17tm19ehUaGoq3tzchISF4eXk5Oxy5a/fu3TRs2JCrV6/SoEEDli9fjoeHh7PDEhFJkNS4tpQtWxZjDAcPHnSYgN+8eZOKFSvi4uKSKvd5z5s3jyFDhnD16lW75StXrqRly5acPXuWAgUKADBjxgxGjRrFv//+i5ubG6NGjWL58uUcPHgvKevYsSNXr17l119/TdD+dT1PGcGEUIJPY5+X25H3N8Mb66z//qgRjLhXh8UCzOJZgLjn+75rPV1pQCmCCXE4//b9y4F4h4+v5wR+zI91PyIi90vMtSVJPd4uLi7MmTOH9evX07VrV6pUqULJkiWpUqUK3bp1Y926dcyZMyfdJ92SNu3bt4/GjRtz9epV6tatS0BAgJJuEZEHpJfRaUFBQVSsWNGWdAP4+/sTGhrKH3/8YWvTqJH9kGJ/f3+CgoJSNVaJKbaq5bH68Ld7SfcHDe2SbrDel92bnxOUdN8/dLwo3jSgVIyEOnr5Ko4naPi4o2JwCR2iLiISlyQVV4tWv3596tevn1yxiMTr4MGDNGrUiMuXL/PUU0+xYsUKcubM6eywRETSnBEjRvDXX38xd+5cnnvuOcDx6LQRI0Y4M0zOnz9vl3QDtufnz5+Ps01oaCi3bt0ie/bsMbYbFhZGWFiY7XloaGhyhy7cS1QTlHx/vBVGr7X++11fGF3XYbO7n9AYy+tTgt84TSQmUfNwxzZ83J/SDhP1mbSiDwGJ3o+ISFweKvEWSU2HDx+mYcOG/Pfff1SvXp1ff/0VT09PZ4clIpImRY9O69KlC/Pnz2f//v2EhITg7e1NpUqV6Ny5Mw0aNEjStkePHs2HH34YZ5tDhw7x2GOPJWn7yeGDDz5gwoQJTtt/ZvFgoho9T3YBcjCaNURGJ9CTg2BEoPXfbzeAN2KfCsw6HDNmMr+ZU2yjJze4E2NIeVwc9cpHYjjGZYfb6EFV/CntcOi6iEhSJTnxjoiIYPr06XzzzTccPnyYmzdv2iqQ7t27l5kzZzJkyBAeffTRZAtWMq+jR4/i5+fHxYsXqVy5MqtXr8bbWxdCEZH4pMTotOHDh9OtW7c42zzyyCMJ2lbBggXZsWOH3bILFy7Y1kX/N3rZ/W28vLwc9nYDjBkzhmHDhtmeh4aGUqxYsQTFJIkTW6LakQp8x58M+3QCDF9tbTyuPrwV++cxuof5MP/xMVvt1kVPB5bYe60d9crHN3w8eio0EZHkkqTE+9atWzRp0oStW7eSN29evLy8uHHjhm19qVKlmDt3Lrlz5+bdd99NtmAlc/r777/x8/Pj3LlzVKxYkTVr1pArVy5nhyUikmnly5ePfPnyJcu2ateuzXvvvcfFixfJnz8/AIGBgXh5efH444/b2qxYscLudYGBgdSuXTvW7bq7u+Pu7p4sMUr8HCWqRfHm+vTNMPRu9fk3n7Em3rGYgj/P8zhF8SaYED5hq10/dVLvtdbwcRFJC5JUXO39999ny5YtfPDBB5w/f56ePXvarff29qZ+/fqa5kMe2qlTp/D19SU4OJjHH3+cNWvWkCdPHmeHJSKSLkRERDBlyhRq1qyJl5cXWbLc+7197969vPrqq/z1118pGsPp06fZu3cvp0+fJjIykr1797J3716uX78OQJMmTXj88cfp3Lkz+/btY9WqVbz55pv079/fljj37duXv//+m5EjR3L48GG++OILvv32W4YOHZqiscvD+fzzzxk7aJT1yet14W1fiKXwrisWW9IN1mR5Fs8mejqw2PSgKicZwnq6atpREXGKJPV4L1myBF9fX0aOHAngsHr5I488wp49ex4uOsnU/vnnH3x9fTl9+jTlypVj7dq1tt4QERGJW1oZnTZ27Fjmz783PVOVKlUAWL9+PQ0aNMDV1ZVffvmFfv36Ubt2bXLkyEHXrl15++237WJdvnw5Q4cOZerUqRQtWpTZs2fj7++fYnHLw5kxYwYDBgywPhn1NLzrFyPptmAtoRZbUp3c91pr+LiIOFOSEu/Tp0/Ttm3bONt4enoSEhKSpKBEzp49i5+fHydOnKBMmTKsW7fOdq+fiIjEL3p02sSJExkxYgQTJkzgnXfesa2/f3RaSibe8+bNY968eXG2KVGiRIyh5A9q0KCBftBPI4IJcTgfdvTyHTN/ZnS/IQD0eW0gMz/IjXHQ0b2E58lHjjiT6vuT5dj2KyKSHiQp8fb09OTixYtxtjl+/Hiy3f8lmcuFCxfw8/Pj2LFjlCpVinXr1lG4cGFnhyUikq5odJqkhDnstk3N5YKFmbSiB1XvLZ+zC/oEADB06FA++egTyli2MoJAu+24YqE2xRKcQMe2XxGR9CJJ93g/9dRTBAQEcPXqVYfr//nnH1asWEG9erFPFSHiyL///oufnx9HjhyhePHirFu3TlVoRUSS4PTp01SvXj3ONhqdJokR23zYOwm2Lp+7G3pZk27L4FoM/WQsZyyhVKMwb/CM7UtnYu/Xjm2/weizKyLpR5J6vEeMGIGvry8NGzZk2rRptmnEbt68SVBQEAMHDiQiIsJuGg+R+Fy6dIlGjRrx559/UqRIEdatW0fJkiWdHZaISLqk0WmS3GKbD/tjgoiavwd6/Gy9aXtgTcwUf6ZZdjCZIFsv9UQaU4PCib5fO7b9fscftOcJDTsXkXQhSYl3vXr1+Oyzzxg8eLBdr7anpycArq6ufPHFF1SrVi15opQM78qVKzRu3Jj9+/dTqFAh1q1bR+nSpZ0dlohIunX/6DQfH58Y66NHp8VXs0UkWk6yOlz+7cJvoPtP1qT71eowtSkuFost6QZrL/UY1nCSIYlOlB3Nww0wjNW8RqCGnYtIupCkoeYA/fr1Y9++fQwYMIAaNWpQunRpqlSpQt++fdmzZ0+MKcZEYhMSEoK/vz979uwhf/78rF27lkcffdTZYYmIpGsjRozgypUrNGzYkC1bttiNTlu7di3+/v4anSaJcp07MRf+3wHousyadPetBtObg8XC8zzhsJf6GJcTtc/ogmoTaWSbWux+GnYuIulFknq8N23ahJeXF5UrV2bq1KnJHZNkIteuXaNZs2bs3LmTPHnysHbtWsqXL+/ssERE0j2NTpPktouz9gsWH4TOP0KUgV5V4fMW4GLBBXiN2izlT7vk2xULZcid4P09WFDtw7vJ9zBW27WLTug15FxE0rIk9Xj7+voyc+bM5I5FMpkbN27QvHlzgoKCyJUrF2vWrKFChQrODktEJMPQ6DRJLsGEMIo19xZ89we8/IM16e5RBWa0tCXdM3mWGhRlJq1svdTJUVBtNGuoS3FcHuj5TmxCLyLiDEnq8c6fPz/ZsmVL7lgkE7l58yatWrXit99+w9vbm8DAQCpXruzssEREMpzy5ctrdFomk5zzXUdv619u3Ou9/v5PePF7iDTQrTLMbAUu1mR4Mc/THuuP6D2oij+lOcblZCuodoM7zKQVfQggEpPohF5ExFmSlHg3btyYDRs2YIxxOC+oSFxu375NmzZtWL9+PZ6enqxatUpDHUVEkpmrqysdO3Zk0aJFzg5FUlFyznf94LYsgPnxEHS8m3R3qQSz7yXdLnfn5r5fUbztkuKE/ijgqKBadM92A0olOaEXEXGWJA01nzhxIpcuXaJ3795cvpy4IhmSuYWFhfHcc88RGBhIjhw5+PXXX6lVq5azwxIRyXC8vLwoVqxY/A0lw0jO+a4dbYufjkCHpRARBS8/CV8/C673vkp+SKM4k+A57KY4U/BjPsWZwhx2x9q2KN5xDlUvijcNKKWkW0TSjST1eL/88sv4+Pjw9ddfs3DhQkqVKkWBAgVi9H5bLBbWrl2bLIFK+hceHk779u1ZuXIlHh4erFixgjp16jg7LBGRDKlmzZrs27fP2WFIKopteHZSCo/F2FbAEUz7byEiCr8XW9Fk3mu84bqeSAwuwEQa8xpP220jmBC28g8ApfChFz/btmiAXvyMP6Vjje1hhqqLiKQ1SUq8N2zYYPt3WFgYhw8f5vDhwzHaaRi6RLtz5w4dO3YkICCAbNmyERAQYFdlV0REktf48eNp0KABCxYsoEuXLs4OR1JBXMOzH2pby/+C57+DO1G0euE5fliwhCyuWehEpViT4jnstku0LfDATwLW50H8Q/s4EuoHh6qLiKRXSUq8o6KikjsOycAiIiLo3LkzP/74I25ubvz000/4+fk5OywRkQwtMDCQBg0a0L17d6ZPn06NGjViHZ321ltvOSlKSU7Rw7OTo/BY9LZ6r5xM1HPfQngk1ds34oeFS8iSJYutjaNtBxNil3RDzKRbRCSzsRhj9LcwDqGhoXh7exMSEoKXl5ezw0l3IiMj6dq1K4sWLSJr1qwsW7aM5s2bOzssERGnSo1ri4tLwsq4WCwWIiMjUySGtCQzXc+DCUmW4dmrVq2idevWhIWF0bzdsyz7ZilZs2aN93XrOYEf8+Nt54KFUwxRj7aIpFuJubYkqsc7KCiIN954g507d2KxWKhVqxbvvvuuimOJQ1FRUfTs2ZNFixaRJUsWvvvuOyXdIiKpZP369c4OQZwkOYZnBwYG0qZNG8LCwmjbti1LvlmSoKQbrMPUHQ0tt9x9RJH4eb1FRNK7BCfeBw4coGHDhty+fdu2bO3atWzdupUdO3bwxBNPpEiAkj5FRUXRt29f5s2bh6urK4sXL6Z169bODktEJNOoX7++s0OQdOLBKb7Wrl3Ls88+y+3bt2ndujWLFy9OcNIN1sR/Fs/aDTePntpMxdJEJLNK8HRiEydO5Pbt27zxxhucP3+e8+fP89Zbb3Hr1i0+/PDDlIxR0hljDAMHDmTWrFm4uLiwcOFC2rVr5+ywREQyhaCgIPz8/PD09MTLy4vGjRuzY8cOZ4cladQcdlOCT/FjPiX4lBHrv6Jlq1bcvn2bRi39+fbbb3Fzc0v0dntQldMM5Vue51ue5xRD6EFVTQMmIplWgu/xLl68OCVLlmTTpk12y+vXr8/Jkyc5depUigTobJnpnrDkYIxh6NChTJ06FYvFwoIFC3j55ZedHZaISJqSUteWAwcOUKtWLbvRaQDZs2fP9KPTdD2PKZgQSvDpvSroG09C8/+Dm3egeVksP7zALPe29KBqrK+/v6dcRCSzScy1JcE93hcuXOCpp56KsbxWrVpcuHAh8VFKhmOMYeTIkUydOhWAOXPmKOkWEUlFGp0miWE3V/fmU9DibtLdtAx83wHj7kofAggmJMZrH+wpn8Nu27pgQljPCYevExHJrBJ8j/edO3fImTNnjOU5cuTgzp07yRqUpD/GGN58800+/vhjAL766iu6d+/u5KhERDKXzZs3U7duXd555x3bsgkTJrBhwwY2btzoxMgkLbLN1b3lFDRbBDfuQJPS8OMLkM36FTESwzEu2/VoBxNCbwJsSXsUhj4E4E9pVnHcti76vu7YesxFRDKTBPd4i8Tl7bff5v333wfgs88+o3fv3k6OSEQk89HoNEmMongzKqgENL2bdDd6BJbdS7rBWn28DLntXmfXU35XJIYg/nGYkKvnW0QkkdOJLVy4kG3bttktO3bsGIDDaaIsFgvLly9/iPAkPXj//fcZP348AJMnT6Z///7ODUhEJJPS6DRJjG3btvGZ/yC4Hk4Vv9os+el7NmU/Rx8CiMTEOuWXraf8vuTbFQsGHCbkD/aYi4hkRolKvI8dO2ZLtB/066+/xlhmsViSFpWkGx9//DFvvPEGAB9++CFDhw51ckQiIiISnx07duDv78+1a9do0KABywOW4+HhQVkKxTvlV1G8mUmrGAl6HYo5TMgf7DEXEcmMEpx4nzhxIiXjkHRo6tSpjBgxAoB3332XkSNHOjkiERHR6DSJz++//06TJk0IDQ2lXr16fPXLIrZ7XLBVJ49+xKUHVR0m6I4ScvV2i4gkYjqxzErTjzj2xRdf2IaUjx07lgkTJjg5IhGR9COlri0uLokv3WKxWIiMjEy2GNIqXc+tdu/eTcOGDbl69Sp169al48oPGJRzbbIWQwsmJM4ecxGRjCIx15ZEDTUXAZg1a5Yt6R49erTt/m4REXEujU6TuOzZs4dGjRpx9epV6tSpw6wV3/BEzjkOq5M/TMKckB5zEZHMRom3JMq8efPo06cPAMOHD+f999/XvfwiImlEiRIlnB2CpFH79u2jUaNGXLlyhdq1a7Ny5Up2eV5SMTQRkVSixFsSbNGiRbzyyisYYxg0aBCTJk1S0i0iIpLGHThwgIYNG3L58mVq1qzJypUr8fLyouzd4eUqhiYikvI0j7ckyLfffkuXLl0wxtCvXz8+/fRTJd0iIiJp3MGDB/Hz8+PSpUvUqFGDVatW4e1t7c2Ork7uivV6fn8xtGBCWM8JzcEtIpJM1OMt8frhhx946aWXiIqKomfPnnz22WdKukVERNK4P//8Ez8/P/777z+qVavG6tWr8fHxsWvjqDr5HHbTm4BkLbgmIpLZqcdb4hQQEMALL7xAZGQkXbt25auvvkpS1VwRERFJPYcPH8bPz49///2XKlWqOEy6oxXFmwaUsvV0RyfdcK/gmnq+RUQejjIoidXKlSt5/vnniYiI4KWXXmLOnDlKukVERNK4I0eO4Ovry4ULF6hcuTJr1qwhd+6E3bd9lMuxFlwTEZGkUxYlDgUGBtK2bVvCw8Np37498+fPx9XV1dlhiYiISBwuXryIr68v58+f58knn0xU0g1Qlty4YH87mQquiYg8PCXeEsP69et59tlnCQsLo02bNixatIgsWVQOQEREJK3Lly8fnTp1okKFCqxZs4Y8efIkqlBaXAXXREQk6SzGGBN/s8wrNDQUb29vQkJC8PLycnY4KW7z5s00bdqUmzdv0rJlS77//nvc3NycHZaISIaS2a4taUFmOufGGK5fv46np2eSC6UFE2JXcE1ERGJKzLVFPd5iExQURPPmzbl58yZNmzZl6dKlSrpFRETSGYvFgqen50MVSru/4JqIiDw8Jd4CwM6dO2natCnXr1+nYcOG/PDDD7i7uzs7LBEREUkiFUoTEUk7lHgLu3fvpkmTJoSGhlK/fn1+/vlnsmfP7uywRERE5CGoUJqISNqhxDuT279/P40bN+bq1as8/fTT/PLLL3h4eDg7LBEREXlIKpQmIpJ2qFR1JvbHH3/QsGFDLl++zFNPPcWKFSvImTOns8MSERGRZNKDqvhTWoXSREScTIl3JnX48GEaNmzIf//9R/Xq1Vm5cmWGr/IqIiKSGRXFWwm3iIiTaah5JnT06FH8/Py4cOEClStXZtWqVfj4+Dg7LBERERERkQxJiXcm8/fff+Pn58e5c+eoWLEigYGB5M6tIisiIiIiIiIpRYl3JnLq1Cn8/PwIDg6mfPnyrFmzhrx58zo7LBERERERkQxNiXcmERwcjJ+fH6dOneLRRx9l7dq15M+f39lhiYiIiIiIZHhKvDOBs2fP4ufnx99//03p0qVZt24dhQoVcnZYIiIiIiIimYIS7wzuwoULNGzYkKNHj1KyZEnWrVtHkSJFnB2WiIiIiIhIpqHEOwP7999/adiwIYcPH6ZYsWKsX7+e4sWLOzssERERERGRTEWJdwZ16dIlGjVqxB9//EGRIkVYv349JUuWdHZYIiIiIiIimY4S7wzoypUrNG7cmP3791OwYEHWrVtH6dKlnR2WiIiIiIhIpqTEO4MJCQnB39+fPXv2kC9fPtatW8ejjz7q7LBEREREREQyrQyReIeFhVG5cmUsFgt79+61W7d//36eeeYZsmXLRrFixfjoo4+cE2QquHbtGs2aNWPnzp3kyZOHdevWUb58eWeHJSIiIiIikqlliMR75MiRFC5cOMby0NBQmjRpQokSJdi1axeTJk1i/PjxzJw50wlRpqwbN27QokULgoKCyJUrF2vWrKFChQrODktERERERCTTS/eJ98qVK1m9ejUff/xxjHWLFi0iPDycr7/+mieeeIKOHTsyaNAgJk+e7IRIU87Nmzdp1aoVmzdvxtvbm8DAQCpXruzssERERHjvvfeoU6cOHh4e+Pj4OGxjsVhiPBYvXmzXZsOGDVStWhV3d3fKlCnDvHnzUj54ERGRZJKuE+8LFy7Qq1cv/ve//+Hh4RFjfVBQEPXq1cPNzc22zN/fnyNHjnDlyhWH2wwLCyM0NNTukZbdvn2bNm3asH79ejw9PVm1ahXVqlVzdlgiIiIAhIeH0759e/r16xdnu7lz53Lu3Dnbo02bNrZ1J06coEWLFvj6+rJ3716GDBlCz549WbVqVQpHLyIikjyyODuApDLG0K1bN/r27Uv16tU5efJkjDbnz5+nVKlSdssKFChgW5crV64Yr/nggw+YMGFCisSc3MLCwmjXrh2BgYHkyJGDlStXUqtWLWeHJSIiYhN9TY2vh9rHx4eCBQs6XDdjxgxKlSrFJ598AkD58uX57bffmDJlCv7+/skar4iISEpIcz3eo0ePdjjk7P7H4cOHmT59OteuXWPMmDHJuv8xY8YQEhJie/zzzz/Juv3kEh4eTocOHVixYgXZs2dn+fLlPP30084OS0REJEn69+9P3rx5qVmzJl9//TXGGNu6oKAgGjVqZNfe39+foKCg1A5TREQkSdJcj/fw4cPp1q1bnG0eeeQR1q1bR1BQEO7u7nbrqlevTqdOnZg/fz4FCxbkwoULduujn8f2q7q7u3uMbaY1d+7c4cUXX+Tnn38mW7ZsBAQEUL9+fWeHJSIikiRvv/02fn5+eHh4sHr1al599VWuX7/OoEGDAOsotegRa9EKFChAaGgot27dInv27DG2GRYWRlhYmO15Wr91TEREMrY0l3jny5ePfPnyxdtu2rRpvPvuu7bnZ8+exd/fnyVLltiGW9euXZs33niDO3fukDVrVgACAwMpV66cw2Hm6UFERASdO3fmhx9+wM3NjWXLltGwYUNnhyUiIpnI6NGj+fDDD+Nsc+jQIR577LEEbe+tt96y/btKlSrcuHGDSZMm2RLvpEhPt46JiEjGl+YS74QqXry43fOcOXMCULp0aYoWLQrASy+9xIQJE+jRowejRo3i4MGDTJ06lSlTpqR6vMkhMjKS7t27s2TJErJmzcoPP/yge9tERCTVJXR0WlLVqlWLd955h7CwMNzd3WMdwebl5eWwtxust44NGzbM9jw0NJRixYolOSYREZGHkW4T74Tw9vZm9erV9O/fn2rVqpE3b17Gjh1L7969nR1aokVFRdGzZ08WLlxIlixZ+O6772jRooWzwxIRkUwooaPTkmrv3r3kypXLdutX7dq1WbFihV2bwMBAateuHes20sOtYyIiknlkmMS7ZMmSdoVYoj355JNs3rzZCREln6ioKPr27cu8efNwdXVl8eLFtG7d2tlhiYiIxOv06dNcvnyZ06dPExkZyd69ewEoU6YMOXPmJCAggAsXLvDUU0+RLVs2AgMDef/993nttdds2+jbty+fffYZI0eO5JVXXmHdunV8++23LF++3ElHJSIikjgZJvHOqIwxDBo0iFmzZuHi4sL//vc/2rVr5+ywREREEmTs2LHMnz/f9rxKlSoArF+/ngYNGpA1a1Y+//xzhg4dijGGMmXKMHnyZHr16mV7TalSpVi+fDlDhw5l6tSpFC1alNmzZ+t2KxERSTcsxlE3sdiEhobi7e1NSEgIXl5eqbpvYwzDhg3j008/xWKxMH/+fDp37pyqMYiISPJz5rUls9I5FxGR5JaYa0uam8dbrIwxjBo1ik8//RSA2bNnK+kWERERERFJh5R4p0HGGN566y0mTZoEwIwZM3jllVecHJWIiIiIiIgkhRLvNOidd97hvffeA2D69On06dPHyRGJiIiIiIhIUinxTmM++OADxo0bB8DkyZMZMGCAkyMSERERERGRh6HEOw355JNPeP311wH48MMPGTp0qJMjEhERERERkYelxDuNmDZtmm3O0nfeeYeRI0c6OSIRERERERFJDkq804Avv/ySwYMHA/DWW2/x5ptvOjkiERERERERSS5KvJ1s9uzZvPrqqwCMHj2aCRMmODkiERERERERSU5KvJ1o/vz59O7dG4Bhw4bx/vvvY7FYnByViIiIiIiIJCcl3k6yaNEiunfvjjGGgQMH8vHHHyvpFhERERERyYCUeDvBd999R5cuXTDG0LdvX6ZOnaqkW0REREREJINS4p3KfvzxR1588UWioqLo0aMHn3/+uZJuERGRdCSYENZzgmBCnB2KiIikE1mcHUBmEhAQwAsvvEBkZCRdunRh5syZuLjotw8REZH0Yg676U0AURhcsDCTVvSgqrPDEhGRNE5ZXypZt24dzz//PHfu3OHFF1/k66+/VtItIiKSjgQTYku6AaIw9CFAPd8iIhIv9XinkieeeIKyZctSvnx5FixYgKurq7NDEhERkUQ4ymVb0h0tEsMxLlMUbydFJSIi6YES71RSoEABNm3ahKenJ1my6LSLiIikN2XJjQsWu+TbFQtlyO3EqEREJD3QWOdUlDt3brJmzersMERERCQJiuLNTFrhirUoqisWvqKVertFRCRe6noVERERSaAeVMWf0hzjMmXIraRbREQSRIm3iIiISCIUxVsJt4iIJIqGmouIiIiIiIikICXeIiIiIiIiIilIibeIiIiIiIhIClLiLSIiIiIiIpKClHiLiIiIiIiIpCAl3iIiIiIiIiIpSIm3iIiIiIiISApS4i0iIiIiIiKSgpR4i4iIiIiIiKQgJd4iIiIiIiIiKUiJt4iIiIiIiEgKUuItIiIiIiIikoKUeIuIiIiIiIikoCzODiCtM8YAEBoa6uRIREQko4i+pkRfYyTl6XouIiLJLTHXcyXe8bh27RoAxYoVc3IkIiKS0Vy7dg1vb29nh5Ep6HouIiIpJSHXc4vRz+1xioqK4uzZs3h6emKxWJwdjkOhoaEUK1aMf/75By8vL2eHky7pHCYPnceHp3P48NLDOTTGcO3aNQoXLoyLi+76Sg1p5XqeHj6fKSkzH7+OXceuY894EnM9V493PFxcXChatKizw0gQLy+vDPuhTi06h8lD5/Hh6Rw+vLR+DtXTnbrS2vU8rX8+U1pmPn4du449s8nox57Q67l+ZhcRERERERFJQUq8RURERERERFKQEu8MwN3dnXHjxuHu7u7sUNItncPkofP48HQOH57OoaRlmf3zmZmPX8euY89sMvOxO6LiaiIiIiIiIiIpSD3eIiIiIiIiIilIibeIiIiIiIhIClLiLSIiIiIiIpKClHhnEGFhYVSuXBmLxcLevXvt1u3fv59nnnmGbNmyUaxYMT766CPnBJkGnTx5kh49elCqVCmyZ89O6dKlGTduHOHh4XbtdA7j9/nnn1OyZEmyZctGrVq12LFjh7NDSrM++OADatSogaenJ/nz56dNmzYcOXLErs3t27fp378/efLkIWfOnLRr144LFy44KeK0b+LEiVgsFoYMGWJbpnMozvbee+9Rp04dPDw88PHxcdjGYrHEeCxevNiuzYYNG6hatSru7u6UKVOGefPmpXzwDykhx3769GlatGiBh4cH+fPnZ8SIEURERNi1SY/H7kjJkiVjvM8TJ060a5NRv2tkhu8H48ePj/H+PvbYY7b1Gel6tGnTJlq1akXhwoWxWCwsW7bMbr0xhrFjx1KoUCGyZ89Oo0aNOHr0qF2by5cv06lTJ7y8vPDx8aFHjx5cv349FY/COZR4ZxAjR46kcOHCMZaHhobSpEkTSpQowa5du5g0aRLjx49n5syZTogy7Tl8+DBRUVF89dVX/PHHH0yZMoUZM2bw+uuv29roHMZvyZIlDBs2jHHjxrF7924qVaqEv78/Fy9edHZoadLGjRvp378/27ZtIzAwkDt37tCkSRNu3LhhazN06FACAgL47rvv2LhxI2fPnuW5555zYtRp186dO/nqq6948skn7ZbrHIqzhYeH0759e/r16xdnu7lz53Lu3Dnbo02bNrZ1J06coEWLFvj6+rJ3716GDBlCz549WbVqVQpH/3DiO/bIyEhatGhBeHg4W7duZf78+cybN4+xY8fa2qTXY4/N22+/bfc+Dxw40LYuo37XyEzfD5544gm79/e3336zrctI16MbN25QqVIlPv/8c4frP/roI6ZNm8aMGTPYvn07OXLkwN/fn9u3b9vadOrUiT/++IPAwEB++eUXNm3aRO/evVPrEJzHSLq3YsUK89hjj5k//vjDAGbPnj22dV988YXJlSuXCQsLsy0bNWqUKVeunBMiTR8++ugjU6pUKdtzncP41axZ0/Tv39/2PDIy0hQuXNh88MEHTowq/bh48aIBzMaNG40xxly9etVkzZrVfPfdd7Y2hw4dMoAJCgpyVphp0rVr10zZsmVNYGCgqV+/vhk8eLAxRudQ0pa5c+cab29vh+sA8+OPP8b62pEjR5onnnjCbtkLL7xg/P39kzHClBPbsa9YscK4uLiY8+fP25Z9+eWXxsvLy3a9Te/Hfr8SJUqYKVOmxLo+o37XyCzfD8aNG2cqVarkcF1Gvh49+PcrKirKFCxY0EyaNMm27OrVq8bd3d188803xhhj/vzzTwOYnTt32tqsXLnSWCwWc+bMmVSL3RnU453OXbhwgV69evG///0PDw+PGOuDgoKoV68ebm5utmX+/v4cOXKEK1eupGao6UZISAi5c+e2Pdc5jFt4eDi7du2iUaNGtmUuLi40atSIoKAgJ0aWfoSEhADYPne7du3izp07duf0scceo3jx4jqnD+jfvz8tWrSwO1egcyjpS//+/cmbNy81a9bk66+/xtw302tQUFCMz7e/v3+6/xwHBQVRsWJFChQoYFvm7+9PaGgof/zxh61NRjr2iRMnkidPHqpUqcKkSZPshtVnxO8ame37wdGjRylcuDCPPPIInTp14vTp00Dmuh6dOHGC8+fP2x2rt7c3tWrVsh1rUFAQPj4+VK9e3damUaNGuLi4sH379lSPOTVlcXYAknTGGLp160bfvn2pXr06J0+ejNHm/PnzlCpVym5Z9EXu/Pnz5MqVKzVCTTeOHTvG9OnT+fjjj23LdA7j9t9//xEZGWn35Qms5+jw4cNOiir9iIqKYsiQITz99NNUqFABsH6u3NzcYtwXWaBAAc6fP++EKNOmxYsXs3v3bnbu3Bljnc6hpBdvv/02fn5+eHh4sHr1al599VWuX7/OoEGDAOtn2dHf19DQUG7dukX27NmdEfZDi+24otfF1SY9HvugQYOoWrUquXPnZuvWrYwZM4Zz584xefJkIGN+18hM3w9q1arFvHnzKFeuHOfOnWPChAk888wzHDx4MFNdj6KPx9F7fv//1/nz57dbnyVLFnLnzp3hzseD1OOdBo0ePdphsZX7H4cPH2b69Olcu3aNMWPGODvkNCeh5/B+Z86coWnTprRv355evXo5KXLJbPr378/BgwdjFFOSuP3zzz8MHjyYRYsWkS1bNmeHI5lIUq4vcXnrrbd4+umnqVKlCqNGjWLkyJFMmjQpBY8g6ZL72NO7xJyPYcOG0aBBA5588kn69u3LJ598wvTp0wkLC3PyUUhyaNasGe3bt+fJJ5/E39+fFStWcPXqVb799ltnhyZpiHq806Dhw4fTrVu3ONs88sgjrFu3jqCgINzd3e3WVa9enU6dOjF//nwKFiwYo2pi9POCBQsma9xpSULPYbSzZ8/i6+tLnTp1YhQyyaznMKHy5s2Lq6urw3Ok8xO3AQMG2IqKFC1a1La8YMGChIeHc/XqVbtfyHVO79m1axcXL16katWqtmWRkZFs2rSJzz77jFWrVukcSopI7PUlsWrVqsU777xDWFgY7u7usV6DvLy8Ur3HNzmPvWDBgjGqWz94bU1Lx+7Iw5yPWrVqERERwcmTJylXrlyG/K6Rmb8f+Pj48Oijj3Ls2DEaN26caa5H0cdz4cIFChUqZFt+4cIFKleubGvzYHG9iIgILl++nOHOx4OUeKdB+fLlI1++fPG2mzZtGu+++67t+dmzZ/H392fJkiXUqlULgNq1a/PGG29w584dsmbNCkBgYCDlypVLl8OWEiqh5xCsPd2+vr5Uq1aNuXPn4uJiPxAks57DhHJzc6NatWqsXbvWVok3KiqKtWvXMmDAAOcGl0YZYxg4cCA//vgjGzZsiDG8sFq1amTNmpW1a9fSrl07AI4cOcLp06epXbu2M0JOcxo2bMiBAwfslnXv3p3HHnuMUaNGUaxYMZ1DSRGJub4kxd69e8mVK5ftR/XatWuzYsUKuzaBgYFO+Rwn57HXrl2b9957j4sXL9qGnQYGBuLl5cXjjz9ua5NWjt2Rhzkfe/fuxcXFxXbsGfG7Rmb+fnD9+nWOHz9O586dM9U1vVSpUhQsWJC1a9faEu3Q0FC2b99um+Ggdu3aXL16lV27dlGtWjUA1q1bR1RUlC1/ybCcXNxNktGJEydiVDW/evWqKVCggOncubM5ePCgWbx4sfHw8DBfffWV8wJNQ4KDg02ZMmVMw4YNTXBwsDl37pztEU3nMH6LFy827u7uZt68eebPP/80vXv3Nj4+PnbVauWefv36GW9vb7Nhwwa7z9zNmzdtbfr27WuKFy9u1q1bZ37//XdTu3ZtU7t2bSdGnfbdX9XcGJ1Dcb5Tp06ZPXv2mAkTJpicOXOaPXv2mD179phr164ZY4z5+eefzaxZs8yBAwfM0aNHzRdffGE8PDzM2LFjbdv4+++/jYeHhxkxYoQ5dOiQ+fzzz42rq6v59ddfnXVYCRLfsUdERJgKFSqYJk2amL1795pff/3V5MuXz4wZM8a2jfR67A/aunWrmTJlitm7d685fvy4WbhwocmXL5/p0qWLrU1G/a6RWb4fDB8+3GzYsMGcOHHCbNmyxTRq1MjkzZvXXLx40RiTsa5H165ds/3/DJjJkyebPXv2mFOnThljjJk4caLx8fExP/30k9m/f79p3bq1KVWqlLl165ZtG02bNjVVqlQx27dvN7/99pspW7asefHFF511SKlGiXcG4ijxNsaYffv2mbp16xp3d3dTpEgRM3HiROcEmAbNnTvXAA4f99M5jN/06dNN8eLFjZubm6lZs6bZtm2bs0NKs2L7zM2dO9fW5tatW+bVV181uXLlMh4eHqZt27Z2PwhJTA8m3jqH4mxdu3Z1+P/6+vXrjTHWKXQqV65scubMaXLkyGEqVapkZsyYYSIjI+22s379elO5cmXj5uZmHnnkEbu/FWlVfMdujDEnT540zZo1M9mzZzd58+Y1w4cPN3fu3LHbTno89gft2rXL1KpVy3h7e5ts2bKZ8uXLm/fff9/cvn3brl1G/a6RGb4fvPDCC6ZQoULGzc3NFClSxLzwwgvm2LFjtvUZ6Xq0fv16h/9vd+3a1RhjnVLsrbfeMgUKFDDu7u6mYcOG5siRI3bbuHTpknnxxRdNzpw5jZeXl+nevbvtR7mMzGLMfXNWiIiIiIiIiEiyUlVzERERERERkRSkxFtEREREREQkBSnxFhEREREREUlBSrxFREREREREUpASbxEREREREZEUpMRbREREREREJAUp8RYRERERERFJQUq8RURERERERFKQEm9Jc44ePUrbtm0pVKgQLi4u+Pj4ODskyaCMMVSrVo0mTZrYLW/QoAEWi8VJUaW9OBLCUaxr1qzBYrGwYsUKJ0UlIiIikjYo8ZYEO3nyJBaLxe6RNWtWihQpQocOHfj9998feh+RkZG0adOGFStW0KJFC8aOHcvo0aOTIfrMYfz48THeowcf48ePd3aYacaCBQvYvXs3b7/9trNDSZKlS5disVj45ZdfnB2KQ40aNaJu3bqMHDmSyMhIZ4cjIiIi4jRZnB2ApD+lS5fm5ZdfBuDGjRvs2rWL7777jmXLlrFmzRrq1auX5G2fOHGCP//8k169ejFz5szkCjnTadeuHRUqVHC4rkGDBqkbTBoVFRXF+PHjeeaZZ3jqqaecHU6S/PTTT+TIkYNGjRo5O5RYjRw5kmeffZbFixfTqVMnZ4cjIiIi8YiKiuKHH37gzJkztGzZktKlSzs7pAxBibckWpkyZWL0mk6cOJExY8bw1ltvsXHjxiRv++zZswAULlz4YULM9J5//nk6duzo7DDStJUrV3Ly5EneeOMNZ4eSJBERESxfvhx/f3+yZcvm7HBi1bRpU/LmzcuMGTOUeIuIiKQDzz33HD/99BMAb7zxBuvXr6dGjRpOjir901BzSRY9evQAYNeuXTHWhYeHM3nyZKpWrUqOHDnw9PTkmWee4eeff7ZrV7JkSerXrw/AhAkTHA6NTui2ALp164bFYuHvv//mk08+4fHHH8fd3Z1u3brZ2ly8eJGhQ4dSpkwZ3N3dyZs3L+3atePgwYMxtleyZElKlizJ9evXGTx4MIULF8bd3Z0nn3ySpUuXOjwv4eHhTJkyhRo1auDp6UnOnDl5/PHHGTZsGFeuXLFrm5hYksu8efOwWCzMmzeP1atXU6dOHTw8PMiTJw9du3bl0qVLDl+3f/9+OnbsSKFChXBzc6NEiRIMHDgwRvvo2xO6devGoUOHaNu2LXny5MFisXDy5EkAbt68yciRIylWrBjZsmWjQoUKzJo1iw0bNti9/yEhIeTIkYMnnnjCYUxRUVGULFmSXLlycevWrXiPfe7cuVgsFtq1a5fg87VkyRLc3d2pVKkS586dsy03xjB37lyeeeYZfHx88PDwoGzZsvTp04fTp0/b2u3atYsBAwZQoUIFvL29yZ49OxUrVmTixIncuXMnwXEAbN68mStXrtC6dWvbsuj7rMPCwnj99dcpXrw42bNnp1q1aqxZswawnsf+/ftTuHBhsmXLRu3atdmxY4fDfRw8eJAOHTqQP39+3N3dKVWqFEOGDIn1c+FI1qxZadOmDb/99hvHjh1L1DGKiIhI6tqxYwcBAQEEBARw+fJl6tevz9ixY50dVoagHm9JVlmy2H+kwsLCaNq0KRs2bKBy5cr06NGDO3fusHz5clq3bs306dMZMGAAAEOGDGHv3r3Mnz+f+vXr24ZER/83Mdu638CBA9m2bRstWrSgVatW5M+fH4Djx4/ToEEDgoODadKkCW3atOHixYt8//33rFq1irVr11KrVi27bd25c4cmTZpw5coV2rVrx82bN1m8eDEdOnTg119/tSvSdevWLRo3bsyWLVsoW7Ys3bt3x93dnaNHj/LVV1/RpUsXcuXKleRYktPPP//M8uXLadWqFXXq1GHTpk0sWLCA48eP89tvv8Vo26FDB1xcXGjdujXFihXjzz//5LPPPmPVqlVs377ddlzRjh07xlNPPUXFihXp1q0bly5dws3NjcjISFq2bMn69eupWLEiL730EpcvX2b48OExhsR7e3vTsWNHvv76a7Zu3UqdOnXs1gcGBnLq1Cn69+9P9uzZ4zxeYwzr16+nXLlyMWKNzfTp0xk8eLDthx5vb2/AmvC/8MILLF26lCJFivDiiy/i5eXFyZMn+fbbb2nWrBnFixcHYNasWQQEBFCvXj2aN2/OzZs32bBhA2PGjGHnzp18//33CYoFYNmyZbi6utKyZcsY61544QUOHDjAs88+y61bt1i0aBEtW7Zky5Yt9O7dm/DwcNq3b8+///7LkiVLaNq0KSdOnLAdE8Bvv/2Gv78/4eHhPP/885QsWZKgoCCmTp3KL7/8wrZt28ibN2+CYq1duzazZ89m3bp1lClTJsHHKCIiIqnrwoULFC9e3Pb9onv37rz++utOjiqDMCIJdOLECQMYf3//GOvef/99A5gWLVrYLX/99dcNYN566y0TFRVlWx4aGmqqV69u3NzczJkzZ2zL169fbwAzbty4GPtI7La6du1qAFO0aFFz6tSpGNurU6eOcXV1Nb/++qvd8iNHjhhPT09TsWJFu+UlSpQwgGndurUJCwuzLV+zZo3D8zJ8+HADmM6dO5uIiAi7dVevXjXXrl1LciyxGTdunAFMu3btzLhx4xw+zp07Z2s/d+5cA5gsWbKY3377zbY8IiLCNGjQwAAmKCjItvy///4zXl5epkiRIubkyZN2+/7mm28MYAYMGGBbFv2ZAczYsWNjxDt79mwDmGbNmtmdoz/++MNky5Ytxmdh+/btBjDdunWLsa3nn3/eAGbv3r3xnqc//vjDAKZTp04O19evX9/c/+cx+rPXtm1bc+vWLbu206dPN4Bp2LChuXnzpt26mzdvmkuXLtmenzp1KsZnISoqyrzyyisGsHsPHMVxv5IlS5oGDRo4bF+3bl1z/fp12/IlS5YYwPj4+Jj27dubO3fu2NZ9+OGHBjCffPKJbVlkZKQpXbq0AWJ8JkeMGGEA88orryQ41n379hnAdOnSxeF6EclYwsPDzbhx40yZMmWMm5ubAcyPP/7o7LAkDYuKijJVq1Y1jRs3tlse1/fC1BT9faZr165OjSMhHjbWK1eumBIlSpiAgABz+fJl07x5c9OxY0e7NocPHzaurq7m888/T4aIMw8l3pJg0f8jly5d2pbEvfbaa8bX19cApkCBAubPP/+0tY+MjDS5cuUypUuXtkuUo/38888GMNOnT7cti+0PbFK2FZ14T506NUb73bt3O0weog0bNswA5sCBA7Zl0Yn333//HaN9iRIlTO7cuW3P79y5Yzw9PY23t7e5fPmyw308TCyxiU6843rs2bPH1j468XaUEEWvmzZtmm3Z5MmTDWAWLFjgcP9Vq1Y1efPmtT2P/swULFjQ7seKaNHJ/e7du2Os6927t8PPQpUqVUyOHDlMSEiIbdnFixeNm5ubqVGjRqzn5n6rVq0ygBk2bJjD9dFJZEREhOnRo4cBTK9evWIkzcYYU758eePq6mr++uuvBO3bkV27dhnAjB8/3mEcD9q7d68BzJQpUxy237hxo93yyMhIkzVrVgPE+BHq9OnTMT4DmzZtsv0g8qBr166Z3Llzm2zZstm9p3El3ufPnzeA8fPzc3wCRCRV3P9jaPQjS5YspnDhwqZ9+/Zm586dybKfDz74wACmXr165vXXXzfjxo0zhw4dSpZtZ3SO3qMHH1euXHF2mMlu3rx5MX7sNyb9JN47d+6M8T3UWZIj1o8++sj2eStXrpwJDg6O0ebll182+fPnN6GhockVeoanoeaSaMePH2fChAl2ywoWLMjmzZvthpEeOXKEK1euULhw4RjtAf79918ADh8+HO8+H2ZbNWvWjLFs27ZtgHU4jaPptaK3c/jwYbvq4D4+PpQqVSpG+6JFixIUFGT3+mvXrtGoUaN4hzInNZa4fPPNN4kqrlatWrUYy4oWLQrA1atXY8S6fft2jh8/HuM1t2/f5r///uO///6zG4ZcqVIl3NzcYrTft28fOXLkoEqVKjHWPf300w4r2/fp04e+ffvyf//3f/Tt2xewTgsWHh5Or1694jlSq+h7lOObI75du3b89NNPvPHGG7z77rsx1l+/fp1Dhw5RpkwZypYtG+9+w8PD+eyzz1i8eDGHDx/m+vXrGGNs66OLC8YnuuDJ/fd3369y5cp2z11cXMifPz83b960DXuPVqhQoRj73rNnD+C4An7OnDmpXr06q1ev5siRI1SsWDHeeHPnzg3Af//9F29bEUl5KTk7CcAvv/xCzpw5CQwMdPi3X+J3/3v0oLRcUDMpMsosIxD7dTktiS/Wa9eu8cknn2CxWDDGUKVKFYoUKRKj3ciRI1m4cCHTpk1Lt4VqU5sSb0k0f39/fv31V8Ca8M6fP59Ro0bx7LPPsmPHDnLmzAnA5cuXAfjjjz/4448/Yt3ejRs34t3nw2yrQIECsW5v+fLlLF++PMHbu/8e2PtlyZKFqKgo2/OQkBAAh3+okiuW5OTl5RVjWfT9+vfPvxwd6+effx7n9m7cuGGXeDt6DwBCQ0MpVqyYw3Wxveall17itddeY/bs2bbEe86cOeTMmZMXX3wxzriiRd8Dfvv27Tjbbdq0iWzZstG8eXOH6xPzPoO12nxAQACPPvooL7zwAvnz5ydr1qxcvXqVqVOnEhYWlqDtLFu2jIoVKzr8EQhifz/jep/vL+4WGhoKxP4eRCfr0e3iE13szsPDI0HtRSRlpeTsJGD9IS9PnjxKuh+Co/coo0rvs4yANZmtWrVqrN9p0pL4Yn3//fe5cOECQ4cOZdq0aezdu9dhu4oVK/Lkk08ya9YsxowZg4uLanbHR2dIHkq+fPl47bXXeP311zl06BBvvvmmbV30l/x27dphrLc1OHzMnTs33v08zLYsFkus25s+fXqc2+vatWuSzkt0T+qZM2cSfGwpFUtyio71wIEDccZaokQJu9c5eg+itxc9WuFBFy5ccLjc09OTTp06sWvXLvbu3cuWLVs4dOgQHTt2tP3oE598+fIB935IiM3atWvx8PCgadOmbN26Ncb66B9iEvI+79y5k4CAAPz9/fnzzz+ZNWsW7733HuPHj0/U6ITTp0+zZ88e2rRpk+DXJFb0+xzbe3D+/Hm7dvGJPs/R511E0p64ZicB6w+RrVq1Im/evLi7u1O2bFnefPNNbt68aWszfvx4LBYLJ06c4NSpU7bZSUqWLJnobQF2s1ts3bqVJk2a4OPjE+OakpTt/f777zRu3BhPT0+8vb1p27atbbYNR8fepk0bChQogLu7O8WKFeO5556LUXw0MbEkp7lz51KrVi1y5sxJzpw5qVWrFvPmzYvRLjHnMzWPNymzjISEhFC/fn1cXFyYPn16ouMPDw9n+vTp+Pv7U6xYMdzd3cmfPz/PPfecbdRXQp04cYIDBw7Y9SA/eK59fX3x9PQkX758vPrqq7YfpJcvX07t2rXJkSMHBQoUYOTIkURERDjcT0Lf58TGer+TJ08yZcoU8ufPz4QJEyhTpgxHjx6NdbaYDh06cOrUKdavX5+oODIrJd6SLF5//XUKFy7MF198YbtwlS9fHi8vL37//fdET5X0oOTcFmCrEH7/8PDkVK5cOby8vNi5c2eMacNSO5bklNyxVqpUiRs3bjj8NdVRohutT58+gLVK+OzZswESPMwc4IknnsDFxYUjR47E2a5KlSqsW7cONzc3mjZtypYtW+zWR08Pd+LECY4ePRrntqKH5rdo0QJXV1e7dZs3b05w7NFT56XkcLboof8bNmyIse7GjRv8/vvvZM+enXLlyiVoe9HnOSHD0kXEuR6cnQTgyy+/pEGDBmzZsoUWLVowaNAgihYtynvvvUfjxo0JDw8HrLenjBs3Dm9vb7y9vRk3bhzjxo1jyJAhid7W/bZu3WqbLrF379688MILD7W9nTt3Uq9ePdzc3OjTpw/Vq1dn2bJlNGrUKMZIqKlTp9KgQQMCAwNp3Lgxw4cPx8/Pj3379sWYSjQpsTysQYMG8corr3DmzBl69OhBjx49OHPmDN27d2fw4MEOXxPX+Uzt4zVJmGXk3Llz1KtXj23btvHNN98wcODARMd/+fJlhgwZQlhYGM2bN2fo0KE0aNCAFStWUKdOHXbu3JmgWMA6Cg0cX5e3b99Ow4YN8fb2pk+fPhQvXpwvv/ySXr16sWTJEp5//nlKlChBnz598PHxYdKkSbz//vsxtpOU9zmxsQKMGjWKsLAw3n77bTw9PalQoQKRkZHs37/fYfvatWsD1o4KSYCUvIFcMpa4qpobY8zUqVNjFAkbNWqUAcygQYNMeHh4jNccOHDAXLhwwfY8riIaid1WdHG1EydOOIy3Vq1axmKxmMWLF8dYFxkZaTZs2GC3rESJEqZEiRIOt+WosFR09eeEVDVPbCyxiS6u9s033ySofXQBtblz58ZY5+i9uHjxovH09DT58uUzBw8ejPGaGzdu2BVGia/Ax6xZs2zV8CMjI23LDx065LCq+f1q1KhhvL29TY4cOcyTTz6ZoOO9X+XKlY2np6fdfqM9+H7u27fP5M2b1+TMmdNs3rzZru3nn39uANOoUaMYVc1v3bplq2q+detWA5gOHTrYtTl48KDJlSuXw/Pk6HPVqFEjU7RoUYfHFFeBs7g+v4CpX7++7fn9Vc0DAwPt2o4ePTrRVc0nTJhgALNq1SqH60UkdSRldpI//vjDZMmSxVSqVMn8999/duuii6h9/PHHdstj+3uT2G1FX4cA8/XXXyfr9h683nbu3DnG9XPv3r3GxcXFFC5cOMZ3iaioKLuZVJJynhxxVMj2/sf919iNGzcawJQvX95cvXrVtvzy5cvm0UcfNYDZtGmTw+N3dD6dcbzxzTLy4HeRI0eOmJIlSxpPT88Y16fExH/79m2HBcMOHjxocubMaRo1amS3PK7vM/Xr1zclS5Z0GDdgli1bZlseHh5unnzySWOxWEzevHnNjh07bOtCQ0NN/vz5Te7cue2+5yb2fU5srNG2bNliAPP444/bvrdGf6+cMWOGw9eEhITYCilK/JR4S4LFl3jfunXLFC5c2GTJksUcO3bMGGP9w9a4cWPbRaR79+5m1KhR5uWXXzaVKlWKUcEyrsQ7sduKL/H++++/bZXKn3rqKfPqq6+a4cOHm/bt25uiRYsad3d3u/aJTbxv3bplnnnmGQOYsmXLmkGDBpkRI0aYdu3amRw5cthVF09sLLFJyHRi93+pSGzibYwxv/zyi8mePbtxdXU1LVq0MMOHDzcDBgwwLVu2NJ6ennafj/gS74iICFOvXj0DmCeffNKMGjXK9OrVy3h6eppWrVoZwEyYMMHha+fMmWO7qCWlimh0Mrhly5YY6xy9n/v37zf58uUzOXPmtLvARUVFmQ4dOhjAFClSxPTr18+MGjXKvPjiiyZ37ty2KXQiIiJMzZo1DWCeeeYZM2LECPPCCy+Y7Nmz26ZCiy/xvnLlismaNat59dVXHR5TciXexhizefNm4+HhYbJmzWpeeuklM2bMGFsV+tKlS5uLFy8meN9169Y1uXLlijEVm4ikrsTOTmKMMYMGDYrxxT5aZGSkyZcvn6lWrZrd8tj+3iR2W9HXoapVqzo8nqRuz1GSEL3u/tku+vXrF2uS+rCxxCa+qub3z2YRPRXlkiVLYmxn0aJFMX4kje98OuN445tl5P7vIjt27DD58uUz+fLlM7///vtDxR+XVq1aGTc3N7vkN7bvM//9959xdXU1gwcPdhi3r69vjO2//fbbBjDdu3ePsS76Pb1/Bp3Evs+JjdUY63eZ6O8oK1assC1funSpAUzfvn1jvCZatmzZzCOPPBLrerlHxdUk2WTLlo0xY8YwcOBAJkyYwIIFC3B3d2flypXMmTOHBQsW8P333xMWFkaBAgV4/PHH6du3b4KHnybntgBKlSrFnj17mDx5MsuWLWPu3Lm4urpSqFAh6tWrx/PPP5/UUwFYz0dgYCCfffYZCxcuZNasWbi6ulK8eHH69u1rd89bcsfy/fff8/333ztc17p160TdU/ygFi1asGfPHiZNmsSaNWsIDAwkR44cFC1alO7du8dahdURV1dXVqxYwbhx4/jmm2/49NNPKV26NJ988gm5c+cmICAg1vuIO3bsyKuvvoqLi0ui9hmtZ8+evPPOOyxcuJA6derE275ixYqsW7eOhg0b0qxZM1asWEG9evWwWCwsXryYJk2aMHv2bBYsWIAxhiJFitChQwdbxXhXV1d++eUXRo8eza+//srOnTspW7YsH3/8Mc2aNYsxhM+RFStWcOfOnRS9vzta3bp12bZtG2+//TarV68mJCSEwoULM3jwYN5880274nlxOXnyJFu2bGHw4MEZrhKvSHqV0NlJ4N5sFqtWrXI4nDRr1qwJmp3kYbZVo0aNZN1eQmfy2LFjBwBNmjRxuP/kiCU29xeyjU1cM1D4+voCOLyVK7bz6YzjTegsI5s3b+aTTz4hX758rFq1yuFMIomJH6zn5qOPPuK3337j/PnzMW5l/O+//2zFRGOzfPlyIiMjEzzLCNwrUBrXurNnz9oKqCb1fU5MrP/3f//Hjh07aNy4Mc2aNbMtj55NJ67t586dW7OWJJSzM38REUfeeOONGL+83i96HsrOnTsneR8vv/yyyZUrV7qZg7JDhw7G29vb4a0WadUbb7xhsmbNahsFIyLO42jk2sWLF82kSZOMi4uLKV++vN1tUMYYU6ZMmTh7YKMf94utxzux24ruNRw7dqzD40nq9hyNqnPUS1imTBljsVgc3pL0sLHEJr7RhfcrXbq0cXFxMVFRUTHWRUZGGovFYsqUKWNblpDzmdrH++OPPxrAvP766w7XR8ccfUvWc889F+P2vaTEv2XLFuPu7m7c3d1Nq1atzLBhw8zYsWPNuHHjbKMo7x8xGVsvctu2bU3u3LnNnTt3HMbt6LMW12jD6JGL69evty1L7Puc2Fhv3rxpihUrZlxcXMy+ffvs1kVERJhs2bKZHDlyxHpec+XKZQoWLOhwndhTcTURcapz587FWPbnn38ybdo0fHx8HP7CCzBp0iQA+vXrl+R9v/vuu9y6dStGRdS0KDw8nF9//ZVmzZqRNWtWZ4eTIFeuXGH69On069eP0qVLOzscEXEgrtlJ4N7sBaGhoXHOZpEQSd1WXDNjJFdsjvj4+GCMcXidSu1YYttnVFSUw9lBLl68iDHG4aix2M6nM443obOMDBgwgB49evDDDz/w0ksvOaz8nZj433vvPcLCwlizZg0///wzn3zyCRMmTGD8+PEULFgw3teDdUrS1atX06JFC4dFCZNLUt/n+8UV68cff8w///xDVFQUlSpVss1GYLFYyJIlC7dv3+bGjRsOi8hGRUUREhKiWUsSSIm3iDhVv379qFy5Mr1792bUqFE8//zzVKlShevXrzN58mTbnNtgnUpr4sSJdO7cmW+//RZ/f39bRc2kKFGiBPPnz8fT0zM5DiVFubm5ERISwjfffOPsUBLsxIkTDB06lLFjxzo7FBGJh6PZSeDebBbRQ4sfRnJuKyW296CaNWsCsHr1aqfH4khcM1BEL3M0nDk2zjjehM4y4uLiwqxZs+jVqxfffvstnTp1ipF8Jyb+48ePkzt3burWrWu3/ObNm+zevTtBsa9Zs4YbN26k6CwjkDzvc2yxnjt3jo8++ogcOXLYqqU/+HjiiScAx8PNjx49SlRUlGYtSaiU60wXEYnfwoULTd26dU2ePHlMlixZjI+Pj2ncuLH59ddfY7SNHrqVM2dO06pVK3Pu3DknRCwikj4lZXaSAwcOmCxZsphy5cqZU6dOxXjNlStXzO7du+2WxTbUPLHbimu4bnJvz9Hw3P379xtXV1dTuHBhc/LkSbv2D1bJTsp5ciQxQ82jq10//vjjJiQkxLb86tWr5rHHHjOA2bhxo215fOfTGcdrTNyzjDwYc1RUlOnTp48BzPPPP283bDox8Tdp0sRYLBa7GVoiIiJsBdpIwFDznj17Gnd39xi3ZziK+36JHWqe2Pc5MbF2797dAOaDDz6IEUu0yZMnG8CMHj06xrr58+cbwHz11Vexvl7uUXE1EXGqTp060alTpwS1bdCgQbIP1RMREavevXvz4YcfsmDBAl5//XVKly5NhQoV+OKLL+jXrx/lypWjefPmlC5dmmvXrvH333+zceNGunXrxowZM+LdfnJuKyW296CKFSvy6aefMmjQIJ544gnatGlDiRIlOH/+PJs2baJFixZ8+umnqRKLI/Xq1WPgwIFMnz6dChUq0K5dO4wxfP/99wQHBzNo0CDq1auX5o+3bdu2jBs3jm3btsVb7NRisfDll1/i4uLCl19+iTGGxYsXkyVLlkTFP3DgQFavXk3dunXp0KED2bJlY8OGDZw5c4YGDRo47F2+X1RUFAEBATRs2JCcOXMm5PQm2cO+z7HFumfPHubPn8+jjz7KsGHDYn19XD3egYGBZMmShZYtWyb9ADMT5+b9IiIiIpIaEtKbOn36dIeFK3fs2GE6duxoChcubLJmzWry5s1rqlatakaPHm0OHTpk1zau6QsTs634emiTc3txTX+5fv1607JlS5M7d27j5uZmihYtatq1a+dwOsrEnCdHEtPjHe3rr782NWrUMB4eHsbDw8PUqFHD4ZRaCT2fqXm8xhhz5swZkyVLFtOvX78ExxwVFWX69+9vK7h2f9HRhMa/dOlSU7VqVePh4WHy5s1rOnToYI4fP+5wOtoHPx/Rc17PnDkz1nMY27lObI93tIS+zwmNNXqK0JUrVzo8hmhnzpwxQIwCajdu3DA5c+Y0bdq0ifP1co/FGHUfiYiIiIiIc3Tu3Jnly5dz6tSpdFF3ZdSoUUyaNImzZ88muBibs6RUrLNnz6ZXr15s3LgxUSMrMjMl3iIiIiIi4jSnTp3iscce46233uL11193djjxeuyxx8iVKxdBQUHODiVeKRFrREQEjz76KBUrVuSnn35Ktu1mdLrHW0REREREnCZ6lpELFy44O5QEOXz4sLNDSLCUiPX06dN06dKFzp07J/u2MzL1eIuIiIiIiIikIM3jLSIiIiIiIpKClHiLiIiIiIiIpCAl3iIiIiIiIiIpSIm3iIiIiIiISApS4i0iIiIiIiKSgpR4i4iIiIiIiKQgJd4iIiIiIiIiKUiJt4iIiIiIiEgKUuItIiIiIiIikoKUeIuIiIiIiIikoP8HvHr9rfYfUIgAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Tabulating the RMSE for the BP-GPR Model\n",
        "\n",
        "Finally, we can tabulate the RMSE for the BP-GPR model."
      ],
      "metadata": {
        "id": "GnnzS2tH9yy9"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "qm_coord_train = torch.autograd.Variable(qm_coord[::nskip,:,:], requires_grad=True)\n",
        "X_train = descriptor.forward(qm_coord_train, atom_types[::nskip,:])\n",
        "shape = X_train.shape\n",
        "X_train = X_train.reshape(shape[0], shape[1]*shape[2])\n",
        "\n",
        "# Make predictions\n",
        "y_preds = likelihood(model(X_train))\n",
        "\n",
        "y_mean = y_preds.mean\n",
        "y_var = y_preds.variance\n",
        "y_covar = y_preds.covariance_matrix\n",
        "\n",
        "def eval(ref, pred):\n",
        "    rmse = torch.sqrt(torch.mean((ref-pred)**2))\n",
        "    q2 = 1 - torch.sum((ref-pred)**2)/torch.sum((ref-torch.mean(ref))**2)\n",
        "    return rmse, q2\n",
        "\n",
        "# Evaluate performance on energy predictions\n",
        "y_rmse, y_q2 = eval(Y ,y_mean)\n",
        "print(f'Energy Prediction RMSE: {y_rmse.item():6>.4f}\\t Q^2: {y_q2.item():6>.4f}')\n",
        "\n",
        "\n",
        "# Evaluate Forces\n",
        "auto_grad, = torch.autograd.grad(y_mean.sum(), [qm_coord_train])\n",
        "rmse = torch.sqrt(torch.mean((auto_grad - qm_gradient[::nskip,:,:])**2))\n",
        "q2 = 1 - torch.sum((qm_gradient[::nskip,:,:]-auto_grad)**2)/torch.sum((qm_gradient[::nskip,:,:]-torch.mean(qm_gradient[::nskip,:,:]))**2)\n",
        "print(f'Force Prediction  RMSE: {rmse.item():6>.4f}\\t Q^2: {q2.item():6>.4f}')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "wHcFRm-p3b_H",
        "outputId": "8cba9a06-1611-47a4-b264-99b4bb12be12"
      },
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Energy Prediction RMSE: 0.0501\t Q^2: 0.9998\n",
            "Force Prediction  RMSE: 4.0714\t Q^2: 0.8093\n"
          ]
        }
      ]
    }
  ]
}