BPGAP_ClaisenRearrangement.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 35,
      "metadata": {
        "id": "STOpB3BJoPPp"
      },
      "outputs": [],
      "source": [
        "#@title Import data from GitHub (mlp_class) { display-mode: \"form\" }\n",
        "#@markdown Install GpyTorch.\n",
        "\n",
        "#@markdown Files from GitHub:\n",
        "#@markdown - **qm_coord.npy** (2100, 14, 3)\n",
        "#@markdown - **qm_elem.txt** ([8, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1])\n",
        "#@markdown - PM3\n",
        "#@markdown  - **energy_sqm.npy** (2100,)\n",
        "#@markdown  - **qm_grad_sqm.npy** (2100, 14, 3)\n",
        "\n",
        "#@markdown - B3LYP/6-31+G*\n",
        "#@markdown  - **energy.npy**  (2100,)\n",
        "#@markdown  - **qm_grad.npy** (2100, 14, 3)\n",
        "\n",
        "%%capture\n",
        "!rm *py*\n",
        "!rm qm_elem.txt\n",
        "!rm -r sample_data\n",
        "!wget https://github.com/cc-ats/mlp_class/raw/main/Claisen_Rearrangement/energy.npy\n",
        "!wget https://github.com/cc-ats/mlp_class/raw/main/Claisen_Rearrangement/energy_sqm.npy\n",
        "!wget https://github.com/cc-ats/mlp_class/raw/main/Claisen_Rearrangement/qm_grad.npy\n",
        "!wget https://github.com/cc-ats/mlp_class/raw/main/Claisen_Rearrangement/qm_grad_sqm.npy\n",
        "!wget https://github.com/cc-ats/mlp_class/raw/main/Claisen_Rearrangement/qm_coord.npy\n",
        "!wget https://github.com/cc-ats/mlp_class/raw/main/Claisen_Rearrangement/qm_elem.txt\n",
        "\n",
        "!pip install pytorch-lightning > /dev/null\n",
        "!pip install gpytorch > /dev/null"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Import libraries \n",
        "\n",
        "#@markdown - math, typing (Sequence, Tuple)\n",
        "\n",
        "#@markdown - Torch (nn, nn.functional, Tensor. TensorDataset, DataLoader, random_split)\n",
        "\n",
        "#@markdown - gpytorch\n",
        "\n",
        "import math\n",
        "from typing import Sequence, Tuple, Optional\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"
      ],
      "metadata": {
        "id": "qe5Z8mkdoYC4"
      },
      "execution_count": 36,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "qm_coord = torch.from_numpy(np.array(np.load(\"qm_coord.npy\"), dtype=\"float32\")) #.cuda()\n",
        "atom_types = np.loadtxt(\"qm_elem.txt\", dtype=int)\n",
        "elems = np.unique(atom_types).tolist()\n",
        "atom_types = torch.from_numpy(np.array([elems.index(i) for i in atom_types])) #.cuda()\n",
        "atom_types = atom_types.repeat(len(qm_coord), 1)\n",
        "\n",
        "energy = torch.from_numpy(np.array((np.load(\"energy.npy\") - np.load(\"energy_sqm.npy\")) * 27.2114 * 23.061, dtype=\"float32\")) #.cuda()\n",
        "energy = energy - energy.mean()\n",
        "qm_gradient = torch.from_numpy(np.array((np.load(\"qm_grad.npy\") - np.load(\"qm_grad_sqm.npy\")) * 27.2114 * 23.061 / 0.529177249, dtype=\"float32\")) #.cuda()"
      ],
      "metadata": {
        "id": "g3bUic8_ojz0"
      },
      "execution_count": 37,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "atom_types = np.loadtxt(\"qm_elem.txt\", dtype=int)\n",
        "elems = np.unique(atom_types).tolist()\n",
        "atom_types = torch.from_numpy(np.array([elems.index(i) for i in atom_types])) #.cuda()\n",
        "atom_types = atom_types.repeat(len(qm_coord), 1)"
      ],
      "metadata": {
        "id": "EDwyLwXCokwh"
      },
      "execution_count": 38,
      "outputs": []
    },
    {
      "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_g2(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\n",
        "    g2 = torch.sum(torch.exp(-EtaR.unsqueeze(dim=2) * (dij.unsqueeze(dim=-1) - ShfR.unsqueeze(dim=2))**2) * fij.unsqueeze(dim=-1), dim=2)\n",
        "    #g2 = torch.sum(torch.exp(-EtaR.unsqueeze(dim=1) * (dij.unsqueeze(dim=-1) - ShfR.unsqueeze(dim=1))**2) * fij.unsqueeze(dim=-1), dim=2)\n",
        "    return g2\n",
        "\n",
        "def symmetry_function_g3(rij: Tensor, Rca: float, Zeta: Tensor, EtaA: 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",
        "    g3 = torch.sum(2**(1 - Zeta.unsqueeze(dim=2)) * (1 + 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 g3\n",
        "\n",
        "def symmetry_function_g3ani(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",
        "    #cosine = torch.cos(torch.acos(cosine).unsqueeze(dim=-1) - ShfA.unsqueeze(dim=1))\n",
        "    g3 = 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",
        "    #g3 = torch.sum(2**(1 - Zeta.unsqueeze(dim=1)) * (1 + cosine)**Zeta.unsqueeze(dim=1) * torch.exp(-EtaA.unsqueeze(dim=1) * (0.5 * (d12 + d13).unsqueeze(dim=-1) - ShfZ.unsqueeze(dim=1))**2) * (f12 * f13).unsqueeze(dim=-1), dim=2)\n",
        "    return g3"
      ],
      "metadata": {
        "id": "WewQkZN4onp6"
      },
      "execution_count": 39,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "class Feature(nn.Module):\n",
        "    def __init__(self, Rcr: float, EtaR: Tensor, ShfR: Tensor, Rca: float, Zeta: Tensor, EtaA: Tensor) -> None:\n",
        "        super().__init__()\n",
        "        assert len(EtaR) == len(ShfR)\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",
        "\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[atom_types].to(device=coords.device)\n",
        "        ShfR = self.ShfR[atom_types].to(device=coords.device)\n",
        "        Zeta = self.Zeta[atom_types].to(device=coords.device)\n",
        "        EtaA = self.EtaA[atom_types].to(device=coords.device)\n",
        "        g2 = symmetry_function_g2(rij, self.Rcr, EtaR, ShfR)\n",
        "        g3 = symmetry_function_g3(rij, self.Rca, Zeta, EtaA)\n",
        "\n",
        "        return torch.concat((g2, g3), dim=2)\n",
        "\n",
        "    @property\n",
        "    def output_length(self) -> int:\n",
        "        return len(self.EtaR[0]) + len(self.EtaA[0])\n"
      ],
      "metadata": {
        "id": "vGbznStaoxUa"
      },
      "execution_count": 40,
      "outputs": []
    },
    {
      "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[atom_types].to(device=coords.device)\n",
        "        ShfR = self.ShfR[atom_types].to(device=coords.device)\n",
        "        Zeta = self.Zeta[atom_types].to(device=coords.device)\n",
        "        ShfZ = self.ShfZ[atom_types].to(device=coords.device)\n",
        "        EtaA = self.EtaA[atom_types].to(device=coords.device)\n",
        "        ShfA = self.ShfA[atom_types].to(device=coords.device)\n",
        "        g2 = symmetry_function_g2(rij, self.Rcr, EtaR, ShfR)\n",
        "        g3 = symmetry_function_g3ani(rij, self.Rca, Zeta, ShfZ, EtaA, ShfA)\n",
        "\n",
        "        return torch.concat((g2, g3), dim=2)\n",
        "\n",
        "    @property\n",
        "    def output_length(self) -> int:\n",
        "        return len(self.EtaR[0]) + len(self.EtaA[0])"
      ],
      "metadata": {
        "id": "xhpIkbU8o2Yi"
      },
      "execution_count": 41,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The Gaussian Approximation Potential (GAP) formalism is motivated by considering an observable, such as molecular energy (E), as a linear combination of unobservables, such as atomic energies (ε): \n",
        "\\begin{equation}\n",
        "E = \\sum_{i}^{\\textrm{Atoms}}ε(\\textbf{r}_{ij})\n",
        "\\end{equation}\n",
        "where $\\textbf{r}_{ij}$ is the relative position between atoms i and j.\n",
        "\n",
        "To use Gaussian processes to determine the underlying function of atomic energies when only the linear combination of functions are avalible, as observations, an NxM matrix $\\textbf{L}$ describes the relationship between N unobserved atomic energies and M observed energies:\n",
        "\\begin{equation}\n",
        "\\textbf{y}=\\textbf{L}\\textbf{y}'\n",
        "\\end{equation}\n",
        "The covariance of the observations can then be given by:\n",
        "\\begin{equation}\n",
        "\\textbf{K}_{MM}=\\textbf{L}^{\\textrm{T}}\\textbf{K}_{NN}\\textbf{L}\n",
        "\\end{equation}\n"
      ],
      "metadata": {
        "id": "hXJ4cge-5PuO"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "class SOAP_kern(gpytorch.kernels.PolynomialKernel):\n",
        "    def __init__(self, natms):\n",
        "        super().__init__(power=2.0)\n",
        "        self.noise = torch.nn.parameter.Parameter(torch.ones((1)))\n",
        "        self.natms = natms\n",
        "    def K(self,\n",
        "        x1: torch.Tensor,\n",
        "        x2: torch.Tensor=None,\n",
        "        diag: Optional[bool]=False,\n",
        "        last_dim_is_batch: Optional[bool]=False,\n",
        "        **params,\n",
        "        ) -> torch.Tensor:\n",
        "        if x2 is None:\n",
        "            x2 = torch.clone(x1)\n",
        "        kern = super().forward(x1,x2,diag=False,last_dim_is_batch=False)+self.noise*torch.eye(int(x1.size()[0]),dtype=torch.float)\n",
        "        L = torch.zeros(int(x1.size()[0]), int(x1.size()[0]/self.natms), dtype=torch.float)\n",
        "        j, k = 0, 0\n",
        "        for i in range(int(x1.size()[0])):\n",
        "            L[i,k] = 1.\n",
        "            j += 1\n",
        "            if j == self.natms:\n",
        "                j = 0\n",
        "                k += 1\n",
        "        LT = L.T\n",
        "        return torch.matmul(torch.matmul(LT,kern), L)\n",
        "    def Ks(self,\n",
        "        x1: torch.Tensor,\n",
        "        x2: torch.Tensor,\n",
        "        diag: Optional[bool]=False,\n",
        "        last_dim_is_batch: Optional[bool]=False,\n",
        "        **params,\n",
        "        ) -> torch.Tensor:\n",
        "        kern = super().forward(x1,x2,diag=False,last_dim_is_batch=False)\n",
        "        L = torch.zeros(int(x2.size()[0]), int(x2.size()[0]/self.natms), dtype=torch.float)\n",
        "        j, k = 0, 0\n",
        "        for i in range(int(x2.size()[0])):\n",
        "            L[i,k] = 1.\n",
        "            j += 1\n",
        "            if j == self.natms:\n",
        "                j = 0\n",
        "                k += 1\n",
        "        LT = L.T\n",
        "        return torch.matmul(kern, L)\n",
        "    def forward(self, x1, x2=None):\n",
        "        if x2 == None:\n",
        "          x2 = torch.clone(x1)\n",
        "        x1_eq_x2 = torch.equal(x1, x2)\n",
        "        if x1_eq_x2:\n",
        "            res = self.K(x1)\n",
        "        else:\n",
        "            res = self.Ks(x1, x2)\n",
        "        return res"
      ],
      "metadata": {
        "id": "M6qYTJyEo3mR"
      },
      "execution_count": 42,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "class BP_GAP_GPR(gpytorch.models.ExactGP):\n",
        "    def __init__(self,train_x, train_y, likelihood, natm):\n",
        "        super(BP_GAP_GPR, self).__init__(train_x, train_y, likelihood)\n",
        "        self.natm = natm\n",
        "        self.mean_module = gpytorch.means.ZeroMean() #ConstantMean()\n",
        "        self.covar_module = SOAP_kern(natm)\n",
        "    def forward(self,x):\n",
        "        mean_x_ = self.mean_module(x)\n",
        "        covar_x = self.covar_module.K(x)\n",
        "        mean_x = torch.zeros(covar_x.size()[0])\n",
        "        j,k = 0,self.natm\n",
        "        for i in range(int(covar_x.size()[0])):\n",
        "            mean_x[i] = torch.sum(mean_x_[j:k])\n",
        "            j += self.natm\n",
        "            k += self.natm\n",
        "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n"
      ],
      "metadata": {
        "id": "WprC2EpppRvP"
      },
      "execution_count": 43,
      "outputs": []
    },
    {
      "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 = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] # H, C, O\n",
        "  EtaR = [[0.0, 0.04, 0.14, 0.32, 0.71, 1.79], [0.0, 0.04, 0.14, 0.32, 0.71, 1.79], [0.0, 0.04, 0.14, 0.32, 0.71, 1.79]] # H, C, O\n",
        "  Zeta = [[1, 2, 4, 8, 16, 32], [1, 2, 4, 8, 16, 32], [1, 2, 4, 8, 16, 32]] # 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], [0.0, 0.04, 0.14, 0.32, 0.71, 1.79]] # H, C, O\n",
        "  descriptor = Feature(Rcr, EtaR, ShfR, Rca, Zeta, EtaA)\n"
      ],
      "metadata": {
        "id": "qgZUGmHwpVAq"
      },
      "execution_count": 44,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "nskip = 5\n",
        "\n",
        "# Prepare training data \n",
        "X = descriptor.forward(qm_coord, atom_types)[::nskip,:,:]\n",
        "shape = X.shape\n",
        "X  = X.reshape(shape[0]*shape[1],shape[2])\n",
        "X = torch.nn.functional.normalize(X) # Preprocess Inputs\n",
        "Y = energy[::nskip]\n",
        "\n",
        "# Initialize likelihood and model\n",
        "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
        "model = BP_GAP_GPR(X, Y, likelihood, 14)\n"
      ],
      "metadata": {
        "id": "py4KHB-IpZ9O"
      },
      "execution_count": 45,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "training_iter = 50 \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",
        "    print('Iter %d/%d - Loss: %.3f' % (\n",
        "        i + 1, training_iter, loss.item(),\n",
        "    ))\n",
        "    optimizer.step()\n",
        "\n",
        "# Print Trained hyperparameters\n",
        "for param_name, param in model.named_parameters():\n",
        "    print(f'Parameter name: {param_name:42} value = {param.item()}')\n",
        "\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "6OS0oZUkpgq8",
        "outputId": "c2664795-4ab6-4c8f-e6e1-925de0b0bb97"
      },
      "execution_count": 46,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Iter 1/50 - Loss: 2.559\n",
            "Iter 2/50 - Loss: 2.531\n",
            "Iter 3/50 - Loss: 2.502\n",
            "Iter 4/50 - Loss: 2.473\n",
            "Iter 5/50 - Loss: 2.445\n",
            "Iter 6/50 - Loss: 2.419\n",
            "Iter 7/50 - Loss: 2.400\n",
            "Iter 8/50 - Loss: 2.398\n",
            "Iter 9/50 - Loss: 2.421\n",
            "Iter 10/50 - Loss: 2.430\n",
            "Iter 11/50 - Loss: 2.415\n",
            "Iter 12/50 - Loss: 2.399\n",
            "Iter 13/50 - Loss: 2.392\n",
            "Iter 14/50 - Loss: 2.391\n",
            "Iter 15/50 - Loss: 2.394\n",
            "Iter 16/50 - Loss: 2.397\n",
            "Iter 17/50 - Loss: 2.400\n",
            "Iter 18/50 - Loss: 2.401\n",
            "Iter 19/50 - Loss: 2.401\n",
            "Iter 20/50 - Loss: 2.399\n",
            "Iter 21/50 - Loss: 2.396\n",
            "Iter 22/50 - Loss: 2.393\n",
            "Iter 23/50 - Loss: 2.390\n",
            "Iter 24/50 - Loss: 2.387\n",
            "Iter 25/50 - Loss: 2.385\n",
            "Iter 26/50 - Loss: 2.385\n",
            "Iter 27/50 - Loss: 2.387\n",
            "Iter 28/50 - Loss: 2.389\n",
            "Iter 29/50 - Loss: 2.389\n",
            "Iter 30/50 - Loss: 2.388\n",
            "Iter 31/50 - Loss: 2.386\n",
            "Iter 32/50 - Loss: 2.384\n",
            "Iter 33/50 - Loss: 2.383\n",
            "Iter 34/50 - Loss: 2.383\n",
            "Iter 35/50 - Loss: 2.384\n",
            "Iter 36/50 - Loss: 2.385\n",
            "Iter 37/50 - Loss: 2.385\n",
            "Iter 38/50 - Loss: 2.385\n",
            "Iter 39/50 - Loss: 2.384\n",
            "Iter 40/50 - Loss: 2.384\n",
            "Iter 41/50 - Loss: 2.383\n",
            "Iter 42/50 - Loss: 2.382\n",
            "Iter 43/50 - Loss: 2.382\n",
            "Iter 44/50 - Loss: 2.383\n",
            "Iter 45/50 - Loss: 2.383\n",
            "Iter 46/50 - Loss: 2.383\n",
            "Iter 47/50 - Loss: 2.383\n",
            "Iter 48/50 - Loss: 2.382\n",
            "Iter 49/50 - Loss: 2.382\n",
            "Iter 50/50 - Loss: 2.382\n",
            "Parameter name: likelihood.noise_covar.raw_noise           value = -0.8300074338912964\n",
            "Parameter name: covar_module.raw_offset                    value = 3.820683002471924\n",
            "Parameter name: covar_module.noise                         value = 0.3428843915462494\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "model.eval()\n",
        "likelihood.eval()\n",
        "\n",
        "# Break predictions into sections to limit memory usage \n",
        "nsections = 5\n",
        "for i in range(nsections):\n",
        "  nstart = int(len(qm_coord)/nsections * (i))\n",
        "  nend = int(len(qm_coord)/nsections * (i+1)) \n",
        "  qm_coord_test = torch.autograd.Variable(qm_coord[nstart:nend,:,:], requires_grad=True)\n",
        "  X_test = descriptor.forward(qm_coord_test, atom_types[nstart:nend,:])\n",
        "  shape = X_test.shape\n",
        "  X_test = X_test.reshape(shape[0]*shape[1],shape[2])\n",
        "  X_test = torch.nn.functional.normalize(X_test)\n",
        "  y_preds = model(X_test) #likelihood(model(X_test))\n",
        "  mean_y_test = y_preds.mean\n",
        "  auto_grad_test, = torch.autograd.grad(mean_y_test.sum(), [qm_coord_test])\n",
        "  if i == 0:\n",
        "    pred_test = mean_y_test.detach()\n",
        "    grad_test = auto_grad_test.detach()\n",
        "  else:\n",
        "    pred_test = torch.concat((pred_test, mean_y_test.detach()))\n",
        "    grad_test = torch.concat((grad_test, auto_grad_test.detach()))\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",
        "y_rmse, y_q2 = eval(energy ,pred_test)\n",
        "print(\"Energy RMSE: \",y_rmse, \"Q^2: \", y_q2)\n",
        "f_rmse, f_q2 = eval(grad_test, qm_gradient)\n",
        "print(\"Force RMSE: \",f_rmse, \"Q^2: \", f_q2)\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "07XSLYsupkIJ",
        "outputId": "d3661927-c3b4-4db4-828b-d33e9181d253"
      },
      "execution_count": 47,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Energy RMSE:  tensor(2.2788) Q^2:  tensor(0.6278)\n",
            "Force RMSE:  tensor(7.6583) Q^2:  tensor(-1.0244)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "print(energy.shape, pred_test.shape)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "S28V1SeQDGBb",
        "outputId": "50da4435-90b8-44de-de12-b8a0077e0c9f"
      },
      "execution_count": 48,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "torch.Size([2100]) torch.Size([2100])\n"
          ]
        }
      ]
    },
    {
      "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(\"energy_sqm.npy\") * 27.2114 * 23.061\n",
        "e2 = pred_test.cpu().detach().numpy() + np.load(\"energy_sqm.npy\") * 27.2114 * 23.061\n",
        "ax[0].plot(e1, e2, linestyle='none', marker='.', )\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(\"qm_grad_sqm.npy\").reshape(-1) * 27.2114 * 23.061 / 0.529177249\n",
        "f2 = -grad_test.cpu().detach().numpy().reshape(-1) - np.load(\"qm_grad_sqm.npy\").reshape(-1) * 27.2114 * 23.061 / 0.529177249\n",
        "\n",
        "ax[1].plot(f1, f2, linestyle='none', marker='.', )\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(\"Reference Force (kcal/mol/A)\", size=14)\n",
        "ax[1].set_ylabel(\"Predicted Force (kcal/mol/A)\", 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": 369
        },
        "id": "Ebm0q3mBBlDm",
        "outputId": "fe3da39b-9ea2-4492-bd79-63782d1b8ace"
      },
      "execution_count": 49,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 720x360 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAFgCAYAAACmDI9oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3iUVfbA8e+dmSQECCRAEsBQFVBAZQEFG1hQUUEELKCyKrq6K+rys6HSNCi2XRfdtbMWVAxCIoiFBaSIu9ISUOlIJCGUNCYhpE27vz+mMJPMJBPIpHE+z5MH8s4773tnV25OznvuuUprjRBCCCGEEMLJUN8DEEIIIYQQoiGRAFkIIYQQQggvEiALIYQQQgjhRQJkIYQQQgghvEiALIQQQgghhBdTfQ8gFNq1a6e7du1a38MQQggfqampeVrr2PoeR32QeVkI0RAFmpebZIDctWtXNm/eXN/DEEIIH0qpjPoeQ32ReVkI0RAFmpelxEIIIYQQQggvEiALIYQQQgjhRQJkIYQQQgghvEiALIQQQgghhBcJkP24++67UUqhlMJkMtG5c2f+8pe/YDabfc7r2rUrSik++eSTStcYNGgQSin+9re/eY79/vvv3HnnnSQkJBAREUHHjh254YYb2LJlS6VrKqVo1qwZnTp1YvTo0SxdurTGn8NqtTJlyhTOO+88WrRoQYcOHbj99tvJzMys8n0pKSlcc801xMbGEhUVxaBBg/jqq698zrn88ss94/T+6tOnj8/9ExMTOfPMM2nWrBnnn38+y5Ytq/HnEEKIpjIvA37nTqUUkyZNqvJ9WmvmzJnD2WefTUREBB06dOCpp57yvL5mzRq/1921a5fPdV5//XXOPvtsIiMjSUhIYNKkSRw/fvykPosQTZUEyAEMGzaMw4cPs3//fubOncvSpUt58MEHK53XqVMnPvjgA59j27ZtY9u2bbRt29ZzzGq1cvXVV5Obm8sXX3zBnj17SE5O5sILL+To0aM+758xYwaHDx9mz549JCUl0bVrV0aPHs1DDz1Uo89QUlJCWloaU6dOJS0tjSVLlnDgwAGGDx+OzWYL+L61a9dy5ZVX8s0337Blyxauv/56Ro8ezbp16zznpKSkcPjwYc/X/v37iYqK4tZbb/WcM23aNN5++23eeOMNduzYwZ///GdGjx7t84NHCCGC1RTmZcBn7jx8+LAn0PaeP/157LHHeOutt3j55ZfZuXMn3377LUOGDKl03vbt232u36NHD89r8+fP58knn2Tq1Kns3LmTefPm8e233/LXv/61xp9DiCZNa93kvgYMGKBPxV133aVvuOEGn2OPPvqobtOmjc+xLl266ClTpuiIiAi9b98+z/HJkyfrP/7xj7pLly761Vdf1VprvWXLFg3ovXv3Vnlv7/d4e/fddzWgV61adbIfS2ut9fbt2zWgf/nllxq974ILLtCPPvpowNc//fRTbTQadWZmpudYhw4d9Jw5c3zOGzNmjL7jjjtqNmghmghgs24Ac2R9fMm8HNh9992ne/bsWeU5u3bt0iaTSe/YsSPgOatXr9aAzs3NDXjOpEmT9JAhQ3yOzZgxQ/fp06dmgxaiiQg0L0sGOQjp6eksW7aMsLCwSq+1a9eOkSNH8uGHHwJgsVj49NNPuffee33Oi42NxWAwkJycXGX2NpB7772XmJgYkpOTPccuv/xyLr/88hpd59ixYwDExMTU6H1FRUVVvuf9999n+PDhdOrUyXOsvLycZs2a+ZwXGRnJjz/+WKN7CyFERU1lXj5+/DhJSUn86U9/qvK8JUuW0L17d5YtW0b37t3p2rUrd911Fzk5OZXOHThwIB06dOCqq65i9erVPq9deumlbN26lfXr1wOQmZnJV199xfXXXx/0mIU4HUiAHMCyZcto2bIlkZGRnHnmmezYsYMpU6b4PXfixIl8/PHHOBwOvvrqK6Kjoys99jrjjDN44403SExMJDo6mqFDhzJ9+nS2b98e1HiMRiM9e/YkPT3dc6xz58507tw56M9ksVh47LHHGDlyJAkJCUG/78033yQrK4sJEyb4fX3Pnj2sXbu20gR/7bXXMmfOHHbv3o3D4WDFihWe0gwhhKippjgvz58/H4vFwl133VXleenp6WRkZJCUlMRHH33EJ598wq5duxg5ciQOhwOADh068Pbbb5OcnExKSgq9evXiqquu8imPGzduHLNnz2bIkCGEhYXRpUsXzj33XF5++eWgxyzE6UAC5ACGDBnC1q1b2bhxIw8//DDXX389jzzyiN9zr732WrTWrFixgn//+99MnDjR73mTJk3iyJEjzJ8/n0svvZQlS5bQr18/v4tJ/NFao5TyfD9v3jzmzZsX1HttNht33nknBQUFnqxKMJKTk3niiSeYP38+Xbp08XvO+++/T4cOHbjhhht8jr/++uv06tWL3r17Ex4ezkMPPcQ999yDwSD/2Qkhaq6pzcvgnD9HjRpFbGzVO5A7HA7Ky8v55JNPGDJkCJdddhmffPIJGzduZNOmTQD06tWLP//5zwwYMICLLrqIt956i+HDh/Pqq696rrN27VpmzZrFW2+9RVpaGikpKaxZs4aZM2cGPWYhTgcSqQTQvHlzzjrrLM4991zeeOMNSkpKmDVrlt9zDQYDd911F7Nnz2bVqlVVZgKioqK48cYbeeGFF/j555+54oormD59erXjsdvt7Nmzh+7du9f4s9hsNsaPH88vv/zC999/77NIpSqLFi1iwoQJzJs3j5EjR/o9x2Kx8PHHH3PPPfdgMvnuXB4bG8vixYspLi4mIyODXbt20bJly5P6DELUh9QMM2+u/o3UDHP1J4uQa0rzMsDWrVvZvHlzteUV4MwOm0wmevbs6TnWo0cPjEZjlZ2JBg0axN69ez3fT5s2jfHjx3Pfffdx7rnnMnr0aGbPns0rr7xyUmUmQtSlupyTJUAO0syZM3n55Zc5dOiQ39cnTpzIunXruPrqq+nYsWNQ11RKcfbZZwfVXmfu3LkUFBRw880312jcVquV2267jV9++YXVq1fTvn37oN73xRdfMGHCBD766KMq77l48WLy8vIq1fZ5a9asGWeccQY2m43k5GRGjRpVo88gRH1IzTBzx9z1/H35bu6Yu94zIR86dIiJEydKW6wGoLHOy27vvfce3bp1Y9iwYdWee8kll2Cz2di3b5/nWHp6Ona7PeDTPXAG4R06dPB8X1JSgtFo9DnHaDTiXKskRMMVaE4G+N///seUKVM85Ua1wVT9KQKcCy969+7N888/z1tvvVXp9e7du5OXl0dkZKTf92/dupWZM2cyYcIET8nB2rVr+eCDDxg/frzPuUVFRRw5cgSr1cqBAwdYuHAh//znP3nooYcYOnSo57w//vGPAAEf59lsNm655RY2bdrE0qVLUUpx5MgRAFq3bu0Za8XrJCUlMWHCBP72t78xZMgQz3vCw8Np06aNzz3ee+89rrrqKr8ZlA0bNnDw4EH69evHwYMHefbZZ3E4HDz55JN+xytEQ7I+PR+LzYFDg9XmYH16PgnNLFx11VVkZWXxyCOP0K9fv/oe5mmtMc7LbiUlJXz22Wc8+eSTPiUabk8//TQbN27k+++/B5wt7vr378/EiROZM2cOAJMnT2bQoEEMHDgQgDlz5tC1a1f69OnjWZi4ePFin0WEI0eO5LXXXmPgwIEMGjSI3377jenTpzNixIhKTwGFaEj8zckDusSwdetWrr/+emJjY5kyZUqlOOVkyb+GGnjssce45557mDJlit/f2Kv6PyUhIYHu3buTmJjI/v37cTgcdO7cmccff9yn0TtAYmIiiYmJhIeHExsby8CBA0lJSeHGG2/0Oa+6DT+ysrJYsmQJAAMGDPB57cMPP+Tuu+/2e5133nkHm83G5MmTmTx5suf40KFDWbNmjef79PR0Vq1aRVJSkt/7l5WVMW3aNNLT02nZsiXXX389n3zyCdHR0VWOW4iGYHD3toSbDFhtDsJMBs6OdgYpmZmZLFu2TILjBqKxzctuCxYsoLi4mHvuucfv64cPH/bJFhsMBr7++mseeeQRhgwZQmRkJFdffTWvvfaaZ12HxWLhiSeeICsri8jISPr06cM333zj06Fi2rRpKKWYPn06WVlZno4fL7zwQlDjFqK+VJyTB3dvy+7du7nmmmuIiopi5cqVtRYcA6im+Fhl4MCBevPmzfU9DCFEI5eaYWZ9ej692xp59K4x7Nq1i2+++YYrr7zypK6nlErVWg+s5WE2CjIvCyFOlXtOHty9Le04xqWXXkp5eTnr1q2jV69eJ3XNQPOyZJCFECKAAV1i6BFjZNiwYezYsYMlS5acdHAshBDi1AzoEsOALjEcOXKEyy4bRlFREWvWrDnp4LgqEiALIUQAx48f57rrrmPLli2kpKQwfPjw+h6SEEI0Kd5Z4QFdqt/E7OjRo1xzzTUcOnSIFStWhKzcTQJkIYTwo6SkhJEjR7JhwwYWLFgQsNWhEEKIk+PuTGGxOQg3GfjsvsFVBslFRUVcd9117N69m2+++YaLL744ZGOTNm9CCFFBWVkZN910E2vXrmXevHmMHTu2vockhBBNjr/OFIGUlZUxatQoUlNTWbBgQVDtEU+FBMhCiAavLpvDWywWbr75ZlasWMEHH3zA7bffHvJ7CiHE6cjdmcKo8HSm8MdqtXLrrbeyevVqPvzwQ2666aaQj01KLIQQDVpNH8GdCqvVyrhx4/jmm2945513PK0QhRBC1L4BXWL47L7BPjXIFWuS7XY7d911F0uXLuXNN99kwoQJdTI2CZCFEA1aoObwtc1mszFhwgS+/PJL3njjDR544IFav4cQQghf7s4UUDkh8um9g3j/xaf5/PPPefHFF3nwwQfrbFwSIAshGjR/zeFrm8PhYOLEiSxYsIBXX32Vhx9+uNbvIYQQouquFd4JEYvVzlNTnuT7L/7NU089VWnznlCTAFkI0aB5P4KLaR7uWcRRW1lkh8PBAw88wCeffMKsWbN4/PHHa+W6QgjR2NW0BVt176uuZM47IXLspy/4fu0nPPjgg8yePbtWP1cwJEAWQtSLmky87tdruxZZa83DDz/M3LlzmTZtGtOmTTul6wkhRFNxsus/qnpfcloW5VYHGv8lc+6EyEt/e41Faz/hzjvv5J///CdKqVB9zICki4UQos65J9C/L9/NHXPXB9WdoibtgIKhteaxxx7jrbfe4vHHHycxMfGUrieEEE3Jyc65gd6XmmFmUWoW2nWe0aD8lsz9unoJi/71PDfddBMffvghBkP9hKoSIAsh6tzJTLzBtgMKhtaaqVOn8o9//INHHnmEV155pV4yFEIIUd8CtdE82TnX+31Gg+JQQanniaHN7gBAAbcM7FQpI52cnMy9997LsGHDSEpKwmSqv0IHKbEQQtS5k1l4568d0MmaNWsWL774Ig888ABz5syR4FgIcVqqqhziZOdc9/uS07JYlJrF5xszSU7LYsaIPj7z/pj+CT7v+89//sP48eMZNGgQixcvJiIiotY/b01IgCyEqHPBTrwV65QrtgM6mWD55ZdfZubMmdx999289dZbEhwLIU5b1bXR9J5za2JAlxhPxth9bXOJJeC8v27dOkaPHk3v3r359ttvadGiRa18vlMhAbIQol5UN/FWldk42cUj//jHP3jqqacYP348c+fOrbfaNiGEaAhC2UbT37X9zftpaWmMGDGCTp06sXz5cqKjo2ttDKdCAmQhRIMUqE55fXo+BwtKa7x5yFtvvcWjjz7K2LFjmTdvHkajsS4+hhBCNFi1Wbrmzf2E7+6LurL98DGu69uh0rVTM8wsXr2Bfz12J9HR0axcuZK4uLhauX9taHABslLKCGwGDmqtRyilugFJQFsgFZigtbbU5xiFEKFXMfsQ0zzckzU2GQ2YDAq7Q3sWgczfkIm5xOJ3u9KtK1OYNGkSI0eOZP78+fW68EMIIRqSky2jCKbXsUM7F+Rt2n8UwDNHA9zy6mIyP34ctOaD5GQ6depUK5+ntjTEnxJ/BXYCrVzfvwz8Q2udpJR6B7gXeLu+BieEqBsVMxveGWW73cG4CzuTW1TO97tymL8hEw0YFISbDMwY0YfEr7djsTko3bmG7K/+zvDhw1m4cCHh4eH1/dGEEKJRcwfB5VYHRoPivku7ERUZVmmuBtCAxepgxpJtOLQm3GTg6q7hHPjsGbTNQsfbX+KQbhhlFd4aVICslEoAbgBeAB5VztUzVwK3u075GHgWCZCFOC14ZzZ2HynCoBSgCTMZ6NOxNTOWbMPunoXBU3Lx3bbDWGwOinasI2/pa/TsN5iUlJR6XxUthBBNwfr0fM+GHzaH5p0f0gEwuYLlcJPBEyQbAINB4dAah4ayY2aSnpuBvbiADuNfoEXH7rVa+1xbGlSADMwBngSiXN+3BQq01jbX91nAGf7eqJS6H7gfoHPnziEephCiLqVmmEn8erunpGLGiD6YSyw+wbGbwaC4rm8HVv3na/KWvkpkwjn8+7MviIyMrIeRCyFE0zO4e1uMBoWtwhxsc2jm/vg7iaP6Yi6xENM83PNn4tfbKS8+zuGFM9H5Wbz98RfY2veu1drn2tRgAmSl1AggR2udqpS6vKbv11q/B7wHMHDgwMo/NYUQjVaK1/akWmtPHVuYUWGx+/5zdzg0uTvXk7P4ZbqefS4ffbGES89J8H9hIYQQ1fLXcjNxVF/PUzzvWdjhmqMnXXGWzzW6tDZy922jseXu58uUFEaOHFm3H6KGGkyADFwC3KiUuh5ohrMG+XUgWillcmWRE4CD9ThGIUQtq2qhx/r0fGKah7Nw84ET25MaT7QLumVgJz7bkOlzveL0LTz6t1n0PLs3P61b02BaBgkhRGMUqK3m7YM606t9FOvT8ykqtTL3x989NcaDu7f1mcNzC4/zxYt/Zd+vm3n+9ffIbHkOqRnmSguqG1ImucEEyFrrp4GnAVwZ5Me11ncopRYCN+PsZHEXsKTeBimEqFWBJt7UDDPj3/sJq11jUHiCYwXcPCDBM4mO6Z9Asld2uSzzV3JSnicspiO2a55mX6FmgMTHQghx0rwX3VmsDuas3MPkYT19gtmr+7Tn6j7tPYHu7iNFzFiyDZtDox128pe+SvGuH7n3qReZl3MGluW7Ky2orklP+7rQYALkKkwBkpRSzwNbgH/X83iEEKfAO1vgPfGWWx28s3Yf/TpFs/VAgad0wq7BqJxfYSYDY722Jx3QJYYZI/owY8k2ijN3kLPoOUyt44m77XnsYS2D6o8shBDCl/c87W65abE6cAA/7s3jp3353HdpNz76ab9PcDvpirNIzTAz3V16oR3kL/sXxbt+pM0VEyk/6wosv+U5g22bg/d+2FfjnvZ1pUEGyFrrNcAa19/TgQvrczxCiNpRMWM8Y0QfTAZnHbEGVuzIZsWObAwVdn++6px4zu8UTUzzcM+GIe5J1FxiofTQbrIXzsTYsg1x457H2CIaBxDTXFq6CSFEsFIzzKSkZbFw8wFsDu0JfGeM6MN7P+wjI7/E07nivXXOzhUVg9uUtCxXcKwxr/o3xb+uIPri24i95Gau69uBTfuPeoLi/fklgLNFp7vf/Zurf/Mpt6ivEowGGSALIZqmirvjmUss3DKwk6ePsZtDOydM958FJRZ+PlDAmj252Oy+j+JiSg9xZMEMjJFRnHH7i6iWbTw9kc0lsqeQEEIEw7u3sXs+ttocpKRlkZyWRZnV4XO+1mA0KJSr9aY7uM0tKgeg8L/zKdq8hC6XjeXxZ2dx0ZntGNAlhl7to5izcg8/7s1zztXAJWe147q+HSqVWwB+y/DqggTIQog64/2oTilFTPNwBndv63fyHdAlhrTMAuwOzcb9Zp/XLDZnHdyITjYe+eMYYtu0ZsyMueTpVmzOMIMGk2sxnxBCiOq5Exjeaz6MRgPbDhZSXmF+BmfG99mRfXzauFlsDkwGxbFNiyn87+e06DuMsEvv9gTH4JzbJw/ryab9Rz07pU4e1rNSAsX9tLC+SjAMdXIXIYTgRM2wu2n8s0u3k5yWxd0XdcXoVVZhMip6xkfh8NPnWOHMLK9av4U7Ro/AYTDxzvwlrDzgDKQd2rmoz+GoPKEL/5RSHyilcpRS27yOtVFKrVBK7XX9GeM6rpRSbyilflNK/aKU6l9/IxdC1BZ3AsPo2pH06t7xoDW/ZBVSeSZ2Lpi+fVBnJl1xFuYSiyeQPZq2DPOquTTveTFtr3sYq0MxZ+UeUjNOJDrcO6U+ek0vT1bY+/5hrk4Y/o7VFckgCyFCzruGzFxi8eyoZLE5mL8h01NO4danQyv6dGxNmGs3JjeDgo6tm/H77+lkJ01FK2h+07N8m6Gx2nwDYruDBrXgo4H7CPgXMM/r2FPA91rrl5RST7m+nwJcB/RwfQ3CubPpoDodrRAiaN7t1tw95P3Ni+6g1XsR9fId2X6vaTQ4kxXuoPdQQSkmg6Jg21ryl/2LZt36027kEyiDEYB1e/PY8PtRPv/TYJ9Msvc4Kt7f/Zq/Y3VBAmQhREh5t2wLMyqevbEvJqNv4OsOjp0bScOvBwvZnb2dZ0f2Yc3uHLYfKuRwYZlzUUdGpjM4ttuIHz8bU5tOfL8zG5NRYfXaNCTMqKTEIkha6x+UUl0rHB4FXO76+8c4F05PcR2fp7XWwHqlVLRSqoPW+nDdjFYIESzvhdEO7ZxjI8IC1/K6g9bUDDNrd+dUcWXF5xszWbj5ACiFze6gNH0TuV//nYiEc4gd/QzKFObzDournrmqILdi0BzoWF2QAFkIEVLJaVmelm0Wu+bNNb8FLH8Id2WM3dnlD/77O+k5x3GfbTuWR/bnT6PLS4gfP5vw2K6AsxXcVb3iiI2KILeonHZREYztnyDZ41MT7xX0HgHiXX8/AzjgdV6W61ilAFkpdT9wP0Dnzp1DN1IhhF/eu5CCMwHhXctbsUOEdxeLiruUerO7shrOpISmNPNXclJeJDy2K3E3z8QQ1szv+xrTNscSIAshQqpCxzYOmksDnlteIav8W85xz/e240fJTnoGe+kx2o97gfD4M33eu2pXNgseuFiC4hDQWmulVI1/tmmt3wPeAxg4cGBj+tkoRKOXmmH22YUUnAvPwrx2unNnl01GA5f3jGXN7hysdh10IBtmVJQd3ktOciJhrePpMG4WRLTwvK440ZEozKh8+tg3dBIgCyFqhXcmAvD8fUz/BJI2HfBkHE6GvbiAnKSp6GIzM/71KdZ2PViw+QA2rwyHzeHMVkuAXGuy3aUTSqkOgPt560Ggk9d5Ca5jQogGZH16PjbXvKuAq3s7+8m7s8Vvrv7N54ldoHrjqlx/ho2P3nqONm3b8smX32EmyrNJCIDRqLiiVxxxURGMcQXHFfscN1QSIAshTlnFTARaY3NoDEqROKovs0b1ZfriX6niiV1A9tJjZC+Yhq0wh6lvzOPZP43xZEYqqpitFqfkK+Au4CXXn0u8jj+klErCuTivUOqPhWh43B0g3K3UHhh6pk9QOrh720rrQWrCaj7Eu+88TUzzZvz441q6d+8OQK/2UaSkZZFbVM6a3Tl8vzObcJOBPh1bN9htpf2RNm9CiFNWsX+l1e7sUmFzaGYs2Uav9lF88eeLOT+hdY2u6yg7Ts4XM7AePUj82GksPBTlyVTbKmSkw00GT4ZC1IxS6nPgJ6CXUipLKXUvzsD4aqXUXmCY63uAb4F04DfgfeDBehiyEKIa7raaF5/Vjhkj+lQKRncfKarU/ce58Uf1bMfyyE6ahkHbmfj8vzEbfbtRvDD6XM7vFI3NoT0/F77bdthvn+OGSjLIQohqVbfVZ1GpFXBmcMNMBux2hydb7HBoUtKy6BgdyW0XdGbnke1BZSwc5SVkL5yJJWc/cWOmEtH1D5RbHSSnZTG2f4InM2I0KG4Z2IkxsijvpGmtxwd46So/52pgUmhHJIQ4VakZZk/GdtP+o/RqH+WzffP0xb9WqjVu3yqCQwVlVV7XXlJI9oJpGK3FxI5/gfl7HSz6fX2ljHDFDLZ7m2n39w29y5AEyEKIKnmXT/h7LPbStzt554d0z/cTL+5K57YtmLFkGw6HxmQysHDzAWwOjcmg6NymOb/lFld5z9hmDvYtfB7rkb3Ej36GZmdeADhXQC9KzaJvx9aM7Z+ABulWIYQ4bVWVvFifnu/pYGGxnsjYrk/P52BBqd+St4N+gmN3+0048VTPfiyHHne/SGnbs3BoPMmL6voa92ofVS89jU+GBMhCiIBSM8zMWbnHM8lW3OozNcPMe+vSfd6z/fAxnrr+HE8d2ob0fE9AbLFrz9+9J11vDms525KeoyxzG7Ejn6DZWYN9XrfZHcxwLQIxGhR9O7Zu8BOtEELUtuqSFzHNw08EtsDe7CL+uWovFpsDFUQdRb+E1vxysNDTp95hKSNnUSKW3AzixkyjpG0vz7nu5EV1CYv66ml8MiRAFkL45T35apwBrVKKmObhnnPWp+ejK0S51/Xt4OmlmbQx02+WQgHxrSI4cqzc57i2Wcn9cjZlGb/SdsSjND/nMp/XDYBBKU/9sXeNc2OZdIUQojZUXPtRcedQc4nFJxHx1c+HANfGTEEsmN6aVej5u3NufoHyQ7tod+OTRJ45sNL5NrvvGOZvyPQkM6ranKShkgBZCOGX9+TrDI7BoTWJX2/3BKSDu7clIsxAudWZkbj/Mucq5lvf/anKtm4a6NymuU+ArO1Wcpe8SNnvqbQZ/ggt+1zheU0puKBLDGfFR9G3Y2tmLNnmCZIdWjfZLaWVUucA44GhQFcgEsgF0oDvgGStdXnACwghmhx3WUVRqRWDUqB1pZre1AwzBwtKcb0MOANjo8HrQJC0w07u0lco27+Fttc9QouzL/V7nkPjSaCkZph92r25Szwa0zwtAbIQwi/vBRZKKewOXanMwl1jlpKW5UlITF38a1Dz76GCUk8DeRx28r56ldLfNtLm6r8Qdf41PudqDZv2m/nlYCFj+yeQOKqvs8ZZa8IbwWKPmlJK9QdeAS4F/gv8D1gElAJtgL7AC8A/lVKvAHMkUBai6au4dTSAUcGMEX0AZ4/hmObhJH69nTJr5cXQ/TtHs2m/2edYi94gJNIAACAASURBVHAjxRa73/tp7SD/u9cp3fMTMVf+iZbnXeP3PHA+4TOXWADnDn7eSRJlUI1unq42QFZKGYDL8Z/BWK61rtyMVAjR6HkvsHBPuIFWHy9MzapxL80s92IQh53i5a9Tsud/xFx5H1H9b/B7vndwPumKsxrVYo+T8CXOAPkWrbU50ElKqYuA/wMexxkwCyGaGO+FeN5P9tzsGhZvyeKXg4WeeTjQA7ye8VFsyTTjPV0HDo415pXvUbxtFa0vvYNWF4wKuHbEAISHnfjZUPGcK8+Oa3TzdMAAWSkVCTyKs8dlG2ArcAhnBqMbMBJ4Vym1HEjUWq8P/XCFEHXJe0FFoID03bX7TrrRvDM78U+Kt60ieuhdtLrgpirPNxpPTMCNabHHSeihtbZUd5LW+ifgJ6VUeHXnCiEaH3fGuNzqbGl536XdCDcZPAun3dxZ4eoe3vXp2JrEUefy5uq9HCooq/L8gnWfUpT2NVEX3ETri8cBcF5Ca7q1a8GSnw+hXSUbf7q0G1GRYT4/G/p29O15f0WvuBp+8vpXVQZ5L87G8ffjzBRbK56glOoC3A4sUEo9r7V+PzTDFELUN38BaWqGme935QR4R9W01hxd/hbF21bS+pLbaT34lmrfc3nP2KYcFHsEExwrpVoA47XWc4M5XwjR+Hi3arM5NHN//J3EUX3ZfqiQzzdmejLFwVYVuzPNFQPsigo3LOLYTwtoed41xFxxL8rV9mLbwUK2HToGGkwG506ptw/qXOn95hKLp4TOoE6UXjQmVe2kN1xrfYvW+ht/wTGA1jpDa/0i0ANYE4oBCiEaLmcXi5rvH621xvz9exzfuoxWg2+m9SWB9qmo8L4a36npUUpdpJSaCxwB5tT3eIQQoTO4e1vnwjoXh9aYSyy8MPpcnr/pXIzBbHvnZeN+M2VewbHCWcN8TvsozzlFW76lYM1HND/7MtpcO8kTHIOznMO9HsXu0Gw/VIg/7jUsRkWjXScSMIOstd4W7EVc2Yu9tTIiIUSD5K8hvXsS9LcYJBCtNQVrPqQodSlRA0cRPeQunwkY4IKuMUQ3Dyc1w8zR4hOZh7ioiKDG1dQopdoCdwH3Ab2Ab4B7ga/rc1xCiNBwz2sxzcO58uw4Vu7MRmswGSsEm8p/Vwqjwm+LzYo0zvN2ZxcBULxjDUeXv01k94G0G/EoymCs8r1JGzPp07F1pSyyv01CGpuqapD7B3sRrXVa7QxHCNEQea+cNhkN3Dwggb4dW2MusXD3RV1594f0oLO7hT9+xrGNKbT8ww3EXHlfpeAYYMuBAme/Y/uJwDvcZGBM/4SA4/LXKL+xU0pdizMoHgFsAF4D3gae0lrvqM+xCSFql7t/fE5ROWv35GKzn2iz6Z5f7XYHu484g9kZXm3UKgomOPbm0FCydwN5X79GRKc+tLvpaZQxrNr32TUBe9E39nUiVdUgbwbP/gBV0UDgXzGEEI2Kv4ys98ppi83B/A2ZnvMN+C998LfaufB/Cyj8XxItz7uGNlc/4Dc4BrDZNQrnYzyDgkvOasfkYT39bqVaVaP8xkwptR8oAz4BntBa73cdf7sehyWECIHUDDPj31/vd8Gz9zxq1zB98a9cdU68pxd8bSjdv5XcJS8R3v5M4sbOwBBW+WldIDZH0+xFX1WA3K3ORiGEaBACZWTdpRT+FnY4cAbJFaf1SsHxhhQK1n1Ciz5XEDt8EhgMAVsRGZWzY4Xd7mwrN3lYT8DZ49NfiUeg9nONXHtgCc4OQtJOU4gmLDkt+FaZdg3Ld2TX2r3LD+4iN+V5wmI6EnfLcxgimtf4GjGukrjGXFJRUVU1yBl1ORAhRP1KzTAzZ+UevxlZdz3Zy9/tZNN+c6Xgt1WkiYJSW8BrH0tdSsGaD+gx+GoemvUG7aKaM/OrbTgCPAecddO5Pm3lAL+Be1Ooc6tCZ5x1x38DPlBKJQGfImsVhWhSUjPMLErN8jlmAJpHGDle7r9HcW2x5KSTs3AmxhYxxN02C2NkqxpfQwHbDhWS+PX2JlXuFvROekqpeGAS0BvnBL0DeEtrXXu/xggh6oV3r013WUPFjOzuI0Vs3O9/z4rSKhbpFW1dhnnlu0T2GIzxiof5flce3dq1wBogOI5tGe6pZ3NPsG+u/i1gKUVjr3MLRGudA7wKvKqUugyYCKzGOW8/oJR6T2u9vT7HKIQ4devT8z3rLRTQtmU4ecctIQ+OrUcPkr1gBiqsGfHjnsfUsk3AcxVgNCoMOLtXKINCa+ei63CTgbyics/Pj6ZS7hZUgKyUugRYBmTj7I0McAfwf0qpa13N6oUQjZS7llfjzFy4a37hRFnDd9sOB3x/eYBHg8d//Z6j//kXkd0HEnvjFEodBn7OKuTnLP+tgQByj1u4Y+56nwxEEy+lqJbWeh2wTin1CM65dyLwsFJqt9b6nPodnRCiOlWVH3jPb0ajgfzjoe8ZbDuWQ3bSNNAO4se9iKl1fMBze3eI4obzOnrmXe8ne+5OG89+tc3zaMtYsdNGIxVsBvlvwOfAn7XWDvBsQf0O8Hfg4tAMTwhRFyoGoO7g2Lus4e6LurJub17Q1yzesZb8716nWdc/EDv6GZSp+hXRbv6yxE24lMIvpZTSFZpMa62LcM677yilzsPZ6k0I0UC5O1Ms3HwAm0N7yg8An/nMPb8dLCj1WQQdCvZiM9lJ03BYSmg/fjZhbTtVeb7F5mDSFWd5Po+b++ndm6t/8ywYVMDNAxKaxBwdbIDcD7jbHRwDaK0dSqnXgC0hGZkQos4M6BLDjBF9+G7bYa7r28Ez6bnLGsqtDvblFXNO+yh2uloMVaV493/J+/rvRCT0JnbMVJQp+J2QFc7tSytmIJpqKUUVcpRSS3Eu1FuutS71flFr/Qvw13oZmRCiWhVL18D5y39yWhYprkV57oB595Eilm8/QoTJQLhRYalpn7Yg2cuOk71gOvbj+cTdOovw+DOrfU+bFuE+nyfQIm53gmWsVzvOxrxwL9gAuRBnV4vdFY53AwpqdURCiDo3f0Omp6fmpv1H6dU+isHd22IyOCdqDazYkV1tz0eAkt82kPfVK0R07EXczTMxhDULehwGd8/7AO3fTjMjgRuB2cBnSqmVOIPlpVrr4FP5Qoh64b1NNDh/+Q8zGVDgs6binbX7WFGLXSkCcVhKyVk4E+vRLOLGzqBZQu+g3tcj3rnLXqC2moGe8DX2PvVVbTXtLQn4t1LqDqVUN9fXncBcnKUXQohGJDXDzJurfyM1w8z8DZlMXfwrNtf2oRbriYnvloGdfILi6nIapemp5C5+kfD47sTd8iyG8Mgqz1c4JyGDApNrO1WNsxn++vT8k/+ATYDWer3W+hmtdR/gD8CPOGuPDyqlflBKPaaUOqt+RymECCSmebjPnDmsdzyf3TeYMf0TMBmU52lZzrGykI9F2yzkpszCcngvsSOfJLJbcHvBGQ3Ks0HT4O5tMRkNrgV7vnXGA7rEMOmKs3wCYH8BdWMSbAb5SZw/yz7weo8V145OIRiXECJEKu6KZ7M7fHYqVQrPxDemfwLJaVl++x9XVJrxM7lfvkBY287E3ToLQ0SLaseiXV8mpbjv0m589NP+03YhXlW01ntxrgX5m1IqFmd2eRSQqJT6HZiitf6mPscohPBlLrFgUM5d6gwK+nWKZkCXGGcdr3JtpaQU3dq1qHLh8qnSdhu5S16mLOMX2t7wfzTvVf2yMYVzzFedHVfhYtr3zyo09sXVQQXIWmsL8Fel1NOAu2Bln9a6JGQjE0KERMXf6itOc1edE19pcdyclXuqXKBXdmAbucmJmKI7EH/bLIzNWtZoTHaH5li57bRbiHcytNa5OJMVHyilmgPXAOX1OyohREUVu1McLCj11OTa7A7P07L84tB1rdDaQf63cyj9bQNtrv4zLfteVe17eneIwmJzkHG0hBU7slm1K4fEUX0xl1g8TxrtDk1yWlaV83VjX1wddB9kAFdA/GuIxiKEqAMVJ22Hw4G7S5vRAA8MPbPSworJw3oGDJDLD+4iZ9FzGKNiiR/3PMbmrYMeizu7ooFFqVmM7Z/gWS19OlNKBW5IWoHWenEoxyKEODnuANHdxSJpYyaLNh+gX6dowJmlVQZFmTU0/Y611hxd/jbFO9YQPeSPRPUfEdT7dhz2XYhtc2hmLNlG4qi+J352GBSLUrOw2auuL27Mi6uD7YMcATwIXAHEUaF2WWt9Ye0PTQgRChV/q1+x/Qjv/pCOc22c4uXvdpKaWeBpAO+e+Nq0CONosdXnWuVHfiN74UyMzaOdwXGL4CdCBQw7J54VO7J96o4b62Ray/KovuTb9YwWY+iHI4SoSqBuDQO6xDgzxg6NQ4PFrn02XLLZNZsCbMB0KrTWFKz9iONbv6PVoLG0GnzLKV3PoTXmEovnZ8ehglI+35jpd/OmpiLYDPL7wAicK6h3IFudCtGouX+rT80wM/fH3z3/oG0VJu9yq4OUtCwGdInh1gGdeOeHdM9rlpx0chZMxxDRgvjxL2CKalfjcVzeK44f9uY22hq1ELqivgcghAhOdd0a3E/tglnLUVuOrV/IsQ3JtOx3HdFD70adYmcgk2tRnvfPjuS0rCY9dwcbIN8IjNJarw3lYIQQdWt9ej52R+ApWwNJGzPp07E1ndu2oFUzE8fKbFjyMsleMB0VFkH8+NmYWsUFvEYgSuGTkWiMNWqhInOtEI1HoPZn3i7rEcvvecX8nnccu/+NR2vNsdSlFPwwjxa9L6fNNX855eBYAZf3jPV0oaiqtVtTEmyAnIPzkZ8QogkZ3L0tEWEGLFYHuHoQVwyX7RqmfvnriUb3Rw+SkzQVpQzEj3uBsOj2Qd/P0+cYCDf5ZiREYK4ytzuA3jj/L9oOfK61lsV5QtSzQN0a3LvoLdiU6VnnYVBwTvsorHYHbVqEk5pprtWA+fi27zGvfJfIHoNpe/1knJsenxqlYNWubFbuzPbJkDf1uTvYAPkZYLZS6m6tde0Xywgh6oz70ZjC2cbts/sG8+7afXy/MxsHYFTODMHmDDPu5LInOC44Qvbnz6C1g/jxLxLW5oyg73tWbAsyzaXY7A4MSjFjRJ8mPbnWFqVUb2AZ0IoTi6T/BDynlBqutd5Zb4MT4jTjr9bYXzbVXXZRZvWNfh2aE7uR5hbX6thK9vyP/G9fp1mX84m98UmUsUZ9GAJyaDw/C5pqvbE/wf6vtxx4AOfWp0dw9kD20Fp3r+2BCSFqX2qGmfHv/eTZxnRhahbPjuzDql05uHc2dWjnzklpmScCZADbsRxncGyzED9+NuHtOtfs5kphszsfQyqcCz5EUF4HtgATtNbHAJRSrYBPgTnAtfU4NiFOG1XVGlfMpqakZVUKjkOp9Pc0cr96hfAOPYgdMw1lCq/1e7h3AmyK9cb+BBsgz8P5aG8OkI0s0hOiUVqfno/VfuKfr9Xm4Ltth33qkDWwIT3f80gQwFaUT/bnU3GUFxM/7gXC47rV+N6lFhsmowG7veku6giRS4AL3MExgNb6mFJqKrC+/oYlxOnFu9bYYnMwZ+UeruvbAXOJpdIWy0kbM+tsXGVZO8hNeYGwNgnE3fJctTuYBnJN73jW7M7B7tAYDQoHzoXbAGFGxS0DOzG2f8JpkT2G4APkq4ErtdYbQjkYIURoDe7eljCj8mSQw0wGruvbgZ/2OdsQuf3m9ejPXmwmO2kq9pIC4m+dRUT7E32KW0eaKCy1VbqPgRO75LkdLCjDqJwbkTww9MzTZpKtBWVAtJ/jrV2vCSHqgHurZasrSF63N491e/MwKHwyyslpWdjrKI1oyd5HzsJnMUa1PalNmtwUcH6naB4YeqanVATwKcc73ebsYAPkTGSnJiEavQFdYvj8/ov8TnozlmzzCZIB7CWFZCdNw16US9wtzxFxxtk+rx/zExwDGIyK+y7pxvs//u6TnbZrWLUrhweGnun3fcKvpcD7Sqk/cSJjfBHwLvBVvY1KiNOR1pUeoXtnlCcP68mp9YwInjX/ANnuVps17ENfkVL4XTR9ugXF3oJd3vh/wCtKKdniSogm4IzoSJ/g+PZBnVnwwEVc0zvec4697DjZC6ZjKzhM7NgZNOvUt9J1AiVJbHbn1tFfPHARdwzqjNHrJ4ZDa0+7IBGUvwJ7gXU4M8ZlwFpgDzC5HsclxGnFveGHPw4NP+7N44656+nTsTWGEEfJtsJsspOmgTIQP+75k2q1KaoWbAZ5IRAB7FZKlQM+aSOtdavaHpgQ4uQF2tWpukUm53eKZvmObBzlxeR8MR1rfiZxY6YT2eX8Go8hr6jck43o07E1M5Zsw+HanU/qj4OntS4ARimlegDuFP5OrfVv9TgsIU4rqRlmDhaUOnsK6xNBskHh0+3HanOw/VBhSLPItuNHyU6ahraWEX/7SzXqJhSIQ8NjX2zl/iFncvugGi7AbqKCDZAfRhbmCdEoeAfBJqOBmwckeBZWVGxon5yW5RNIF5VacVhKyVn4LJbsdGJHP0Nk9wEnNY41u3NIzXB2hTSXWEgc1bfSYhYRPK31XpyZZCFEHfKeUw0GhdEVIyuFT62xu8uDhpDVINtLj5GzYBr2YjPxtz1/UgumA9mfX8IzX/5KZn4xUZFhp/1cXWWArJQ6T2v9i9b6ozoajxDiFFVcaf35hkxS0rL47L7BPg3tjQbFF5sPYLNrjApm3XQuSzank5OcSPmh3bQbNYXmZw2q0b0VJ36Ttjs0KWlZJKdlBdyCVQRHKTUa5/bTcVQojdNa31ovgxKiCfN+Cuc9pyqHZtyFndHAgk2ZPqnDs9tHcazMytKtB0MyJkd5CTkLZ2I1Hybu5pmV1oTUlnd/SAcgIuz0nrOrq0H+QSmVrpT6h1JqqKqNLVmEECHlDoLdj/jcj/3cdb9j+ydw5TnxxLZq5mnhY9fw9MJUtn44nfLMbbQb8Rgtel1So/t6P1I0cCKTUnELVlEzSqm/AwuAc12H7BW+hBC1yJ0x/vvy3dwxdz0xzcMxGRQKMBoUY/oncEZ0ZKUd8HYeKeJgQRlF5bX/z9JhLScnORHLkX3EjnqKyK79av0ebu4ORBbr6T1nV1diEQtcCdyIsyl9pFLqG2AJsExrXRLi8Qkhasi9q1NKWhYLNx/A7tCEmQwUlVq57d2f/C4y0XYruV/Opmz/FtpeP5kWvYfW+L7uqyrgkh7tmDysJ+BsmF9xC1ZRI3cBt2itl9T3QIQ4HVQsRdt2qNBZT4GrrgKIaR7u88QslLTdSt6Slyg/sJ12Ix+jeY+aPdkLxGgAUGjt7HvctkU42cfKPZ/JYFCn9ZxdZYCstbYC/3F9TVJKXYAzWH4O+FQp9T3OYPkrrXXOqQxEKdUJ54Yk8Tj/m3tPa/26UqoNzuxJV2A/cKtsdy1E1dyL48b0T2B9ej4xzcP9tnED0HYbuUtepjR9M22ufYiW5w6r9vpGg8LhqNzuyP3a5GE9PY/lKm7BKmqsBNhV34MQ4nThfgpnsTpAKTam52OzO5y1xXbn2o2Fmw/UTXDssJP39WuU7ttEm2sn0aL35bVyXQXMGnUuvdpH+fQ99tRaK0XiqL6n9Zxdo426tdabgE3AdKXUmTiD5QnAm0qpR7XWb57CWGzAY1rrNKVUFJCqlFoB3A18r7V+SSn1FPAUMOUU7iPEacMdKL+5+jcc2k9w7LCT9/XfKd27nphhDxDVb3hQ11VoDAbl0+PY/fix4qRasa+mqLGXgCeVUg9orf03ng4hpdR+oAhnOYdNaz1QEheiKRvQJYYZI/owfck27A7t2TjJvQgvr6jcZ0fSUNFac/Q/b1Kyax3Rl99NVL/rau/awOrdOdw+qLPP/CwJjRNqFCB701rvA/4B/EMp1RZocyoD0VofBg67/l6klNoJnAGMAi53nfYxsAYJkIWokcHd22IyKM+krnEGx/nfznFNvhNpNWBk0NdzaOfk7aaAS11lFaf7pBoC7wMjgYNKqT2A1ftFrfWVdTCGK7TWeV7fP4UkLkQTZi6x4Kj4xE3BjBF9nCUXIaa1xrz63xz/ZTmtLrqV1oNurvV75BwrY/6GTL7bdpjr+nbwBMsyhzsFDJCVUmOCvIbWWn8J1Folt1KqK/AHYAMQ7wqeAY7gLMHw9577gfsBOneWHn5CVOKqoTMYFDa7naP/eZPi7auJvmwCrQcF+8/dpcLPjYplFaJWvQNcCiwDsmkYLTclcSGatMHd2xJmMmCxnViJpzUkfr29cuAcAoX/S6Jo02Ki+o8g+rIJIblHt3YteObLXwHnttmA9ED2UlUGeVGQ19CAsRbGAoBSqiWQDEzWWh9T6sTaeK21Vkr5/S9Ta/0e8B7AwIEDG8IPECHqlXebomTXQjl3Dd3RFe9y/JfltL7oNlpffFuNr909riVZ5hKpVasbtwGjtdYr6un+Gljumnvfdc21krgQTZZ77px4cVe+SM3iaLHF81qZ1VHFO2vHsc1LKPzxM1r0vYqYYffjHQfVljCj4ve8Yp9j3207LAGyl4ABsta6zlu6KaXCcAbHn2mtU1yHs5VSHbTWh5VSHYBTWgwoxOlg/oZMZrjq58KMCrt2LqjTWmNeNZfjW76h1YVjaH3ZnSd1/YmXdKNX+yhS0rLQQK/2UbU6fuEjFwhNY9XgXKq1PqiUigNWKKV8FgxK4kI0Je4Wb+VWR708qjn+y3LM379P854X0/a6RwhVd12HQxPfqhlwolzkur4dQnKvxuqka5Brm3L+ivRvnFuovub10lc42xy95PpTWh0JUYXUDLNPxwqLu+5Yawp++JiizUuIGjCS6MvvqXFmon2rCG7qdwbmEgu7jxR5NgFxb0QiWeSQmAkkKqXu1lofr+uba60Puv7MUUp9CVyIJC5EE+POGh8sKMViqxwct2pm4lhZaNfIFu/6kfxl/6JZ1z/QbuQTKMOpP5z3bkVnAEwmA3a7s+3mA0PP5PJecT41yOKEoANkpZQJ58TYGQj3fk1rPa8WxnIJzo4YvyqltrqOPYMzMP5CKXUvkAHIrlFCVGF9er7fjhWF//2cY+sX0bLfcGKuqvljO6OCR67qSeLX2z2lFXZXqzf3JiASIIfEEzi7RWQrpTKpvEjvvFDdWCnVAjC4Fk63AK4BEpHEhWgC3EFxTPNwz7xmMhowGZQnseB2dvsoNu4PXaOW0n2byVv6NyI6nk3s6KkoU9gpX/PCrjFszSrE4to5ddaovj5t3dwL8iQw9i+oAFkpdTawFOiG8xcSu+u9VqAcZ//iU6K1/hHfzbi8XXWq1xfidOHuWOGe4BVQsH4hhf+dT4u+w2hzzYMnVdOmDIrthwo9DfTB2Vxeay2bgIRWsOtBQiEe+NL134sJmK+1XqaU2oQkLkQj5i6lqPjLvsXmoE3zMM6Ma8mWAwXY7RqlYHMIg+OyA9vIXTyb8NguxN08A0N4s1q5bo/4KDZnOMftcGi2HSqs1NZNBBZsBnkOkAr0w7kgox/QGngbmBaaoQkhquOdATGXWIhpHs7q3Tk+PToLNy2mYO3HNO89lLbXPRx0Tds1vePRwMod2a62cM4fIOEmg2dnvBkj+mAusUjPzBBQSsVorc1a6+fqawxa63TgfD/H85HEhWjE1qfnexbcVXzidrTEytH9ZtpHRZBdVE4om1aUH95LzqLnMLaKI+7WRAzNWtbataMiTJ5kiQYWbj7A2P4JMlcHKdgA+QJgqNa6WCnlAEyuDT2eBP4JhOwRnxDClzsoLiq1MvfH3z2ZD3+K0r7BvGouzXteTLsbHg26ps1kVDww9EwA1u3N9QTEY/snMNa1O58ExSF3RCn1X5zlC4u11hn1PSAhmoqiUmu15xwpKg/pGCx5meQsnImhWRTxtz2PsXnrWr3+T+n5XN4rjuU7sgGw2jUpaVkybwcp2ABZ4dzuFJwrqs8AdgNZwFkhGJcQwo+arLAu+nk5R1e8TeRZF9LuxuoXfJzTPopdR4o82eL16flMuuIsvzsryQRbJ7rg3K10FPCSUmo3zmB5idY6rV5HJkQjt/3wsXq9v7XgCDkLpqEMRuLHPY+pVbtav8eOw8eIa+VbriGtZIIXbIC8DedjtnRgIzBFKWUH/gT8FqKxCSEqWJ+e73eFdUXHt63i6LJ/0qzbAGJHPY0yVr/g41iZlTBX+YTReKKmWHZWqh9a6yM4W6S95+oPPxxnsLxSKVWEc6HcEmC11tpefyMVovFwP4Gr/c7CwbMV5ZGdNBVtsxJ/+4uExXQMyX3sDk1cVAThRucuqmFGxdj+CSG5V1MUbID8AtDC9fdpwDfAaiAPWZwhRJ2puADPn+Kd68j/dg7NupxL7Ohngl4NbbE5nFtFwYk/RYPgau+2CFiklDLi3MXuRmAu0Eop9bDW+rN6HKIQDZ73wrw62AzPL3tJITkLpuMoPUb8uBcIj+0asnsZlGJM/wT6dGztaeUmyY7gBRUga63/4/X3dOAcpVQbwKy1/CQVItRSM8wkp2WhgH6dogO2GyrZ8xN5S18l4oxziB0zA0NYRND3yC+2oLXzEZzVrkmWWrUGyZUt/t719VelVD/g1HtCCdHEuZ/A1Vdw7CgvJueLGdgKs4m75TkiOvSs9XsYlDO/oRTcd2k3AE8Lu037j9KrfZTM60EKts1be5wL87Lcx7TWR5VSCUopq9Y6O2QjFOI0l5ph5rZ3/4fNtcNpoEeDJfs2kbvkZcI79CDu5pkn1SrIoMDuCpIXpWbJiud6opTqH+SpWmu9JaSDEaIJmL8hk+Xbj2BQ1EuA7LCWkbPoOSy5+4kbM41mnc8NyX2Ucv6McGj46Kf9HCu3eX4pkH71NRNsicWnwALg/QrHrwVuw9k8XghRy1IzzCQu3e4JjsH/IovS37eQ++VswuO6En/LcxgimtfoPgblbN82pEesZ8WzzS6TaT3ajPP/6upKJTVw6tttCdGEzd+QyTNf/ur5vnmYgRKro4p31C5ts5L7jet0CgAAIABJREFU5WzKs3bS7sYniDzzgpDdy+H6WO4NnBS+rTmlX33wgg2QBwKT/BxfB7xae8MRQri56+XKqpnIyzJ/JTflecLanHHSfTQvOasdk4f1ZPeRIk+A7NAQ0zy8mneKEOlW3wMQoqlYsCnT5/s6DY4ddvKWvkrZ72m0Gf4wLc4ZEtL7mYzKuZubw7mB05j+CYyR1pwnJdgA2QT4K2ZsFuC4EOIUuevlqlKWtZOcRc9hah3v7KMZ2arG9wk3KiYP68mALjGe1d0aMADmEstJjV2cGul5LMSpcXerGNy9LeGm4DZHqm1aO8j/7p+U7PkfMVfeR9T514b0fgq4ZWAnv73qJTCuuWAD5A3AX1xf3iYBm2p1REKcpuZvyPSsNL59UGfPxB6o53H54T3kLJyJsWUb4sY9j7FF9End99kb+3omz8Hd2xIRJo/jGiKlVEegM+CT1tda/1A/IxKi4UnNMJOSlsXCzQewOTQGpegZX3u70wVLa435+/cp3raS1peMp9UFN4X8nu42btKas3YEGyBPBVYppc4DVrmOXQn8ARgWioEJcTrxrpFbtzePzPxiisptXNYjlpxjZfycVehzviU7nZwF0zFGRhE/bjamlm1O6r4K3yzxgC4xfjcGEfXHFRjPB4Zwoi7Z+3cmqUEWAv8bKTm0ZsfhojofS+GPn1GUupSogaNofcntIb+fO3ssc3btCbbN23ql1EXAE8AY1+EtwINa659DNTghThffbTvs8/07P6QHPNeSu5/sBdNQEc2JHz/7lHZgCjOqSlliyT40OHMAO9Ab5xO74UA8kAj8Xz2OS4gGJdiNlEKtcEMKhf9LosW5VxNz5X0oFfptSYwGZ89jUXuCbfNmcAXCd/p5rZXWun73bBSikbuubwfW7c2r9jxr/gGyk6ahjCbix72AqXV8je7TtW1zMo+W4NCScWhEhgI3aK13KaU0kKu1/q9SqhyYBayo3+EJUf9SM8xsPVBQbz2O3Yq2LqNgzQc073UpbYc/FJLgODrSREGpzfO9UUHiqL4yl9eyYEss/g3cU/GgUqo1sBwYVJuDEuJ0kpphZvuhQtpHRXCkqDzgeVbzIbKTpoKC+HEv1Hh7UoOC+4ecSeLX2z01xpJxaBQice5aCnAUiAP2ADuA8+prUELUJ3etsQb6dmxN4tfbKa/D7hT+FO9Yy9H/vEmz7gNoN/IxlCE01U+FXsGxAsZd2JnbB3UOyb1OZ8EGyIOVUq9qrZ9wH3AFxyuAwsBvE0JUJTXDzPj311fbrcJWmE3251PRdhvx42cT1rZTje91/2XduX1QZ3q1j5Ia48ZlF3A2sB/YCvxZKXUA5yLpg/U4LiHqRcV50+ja/KM+k8clv20k75vXiOjUh9ibnkYZQ7e5pffnNBgUGuf/JjKf165gA+RrgR+VUvla65eUUtGcCI5HhGx0QjRx69PzsVYXHB/LI/vzZ9CWEuLHv0h4bNca36dVMxNX92kPSI1xI/Q60N7190RgGTAeKAfuqq9BCVFfKs6bDo1nhzyDgvbRkRw0l9bZeMoyfiF38YuEx3UnbuwMDGE138W0Joxe20krBUkbM0lJy+Kz+wbL3F6Lgl2kl6mUuhb4QSllBcYBBfw/e3ceH1V9Ln788z2ThQRCCIEEMCyiiBJcCqhY9wWXFkRBBLWbrYr36u1i722rbMqitrVqb8ttXdra/hRZBEVxx+IuoKGoLFIRCYQlC0xCyDbLeX5/nJkwSSbJBGbJ8rxfL15kzpyc840v/M6T5zzf5wsTRKT5Z8JKqRaNHZpNcpLVbAbZd/ggxYvvwV9TSe60+aTkDo3ouo23Uz1U6+OGxz/i2dvO0Qm0gxGRZ0K+3mCMGYKTUd4lIq0XrivVSQR7G2elp5DkMnj9RyY5CfzxC3ENjuv2bqNkxTySe/Un5/q272LaVkkWXHJyLn0yUjHAs+t36TbSMRJpBhkR2WqM+TawGvgIuFqDY6UiEzqxu6s99Z0j1u44wA+/OYQn3tuBv9HzQX9VOSWLZ+A/fJCc6+eR2v+kiO8nAifm9GB7yeH6Y16/6ATaARljUgBLRGoBRKQa2GCM6WaMSRER3c1FdXrBFm4en02Sy8JvNwqOE1Bf4Sn52ulFn55JztR5R7VRU1sZY1i9tZiUJIvZ4/N1G+kYajZANsZ8TviSHh8wBPg4uDpTRHShiFLNaNyb0zJOSx6/X7AJPC5r9D3+mkMUL5mJr6KEnCn30i3vlDbdU4DLTs5h18Hq+ux0uJZuqkNYBrwDPNzo+O3ARUDsdyBQKsGCLdxsodU1G/HgPbiH4qWzMEkp5EydT1JGfObWYNbc67NxV3u0b30MtZRBfi5uo1CqE2vcm9MWsEPSxY0zx3btYUqWzMJ7cA85182h26BTj+q+O8qqePbWsfUrvYM7LKkO51yczZoaexO4J85jUSohxg7NJskyeBpPmAngO1RK8ZKZYNvk3vgAyb36tf5NUWIMWFCfMdY1JbHTbIAsIvfFcyBKdVbBLaM9XpvW8h52XTXFS+fgKS0kZ9IM0oaccdT3fWtrMdMvPIEF1x5dgK3ajXScJ3eN2UBGnMeiVEKMHpzFRcNzeGNLcULH4a8qp3jJTOzaKqejUJ+2dxQ6FtPPH0pGWrJmjOPAau4N08bu1m09X6muIrh98w1nD2r+fzjA9tRS8ty9eIq303fiL0k74cw23cfC6YkZJOJkr1WH9xlO14rGbgQ2xXksSiXEonW7WL01scGxXXuY4qWz8R8qI2fKHFL7nRjze7oswzVnDOD0vEwuH5HLuPx+3HHxiRocx0FLn9dfGGO+Y4xJbekCxphTjDFPAL+K7tCU6hyCC/RKKuuazSDb3jpKV8ylbs8X9Bn/36SfdE6r17UMXD4ilxSXwWUgJdli+gVDSbIMFs5rrTnuFOYC9xhjnjHG/CjwZxHOnKtP+lSnUVDoZuGa7RQUuhscX7RuFzOe/zyhu+TZnhpKlt2Lt2wXfa+9h255+TG/5yn9Mlg6/Ry+e84QthVXsnprMVMf+4hF63bF/N6q5Rrk24DfAH80xrwFfALsBWqBLGAEcB5wEvC/wB9jO1Sl2r9gMBx8/BW68rq5Vdbi81L6/P3UFn5O9vi76H7K+U3OaZAZDvydZBmmX3gC0y88ocE9x+X300UbnYiIvGKMmQDMxJlrAf6F00no1cSNTKnoCZ0rU5Ks+p6+BYVuZr3weUI3ARGfh9IVC6jb92/6TPwlaUNHx+W+NV4/owdnsXDN9voFirYIs1duYni/DJ3fY6ylGuR3gLONMd/Eebw3FRjMkW1P/wX8FXhaRMrjMFal2rVwE3zoyutwxO+ldOUD1H5dQO8rf0yP/IubnGOAcSNy+ecXJfhCLuQLtG1r/LhNF210PiLyGs4GIQ0YY5JFxJuAISkVVaFzZZ3XZsWGIkYPzmLtjgNNFjLHk9h+Sl/8DbWFG8n+1s/oPvzcuN37jIG9AGcdi2UMdiDLYtvasjMeWiqxAEBEPhSR/xKRM0QkS0S6iUieiEwQkT9qcKyUI3SC9/qcCX7j7vLmg2PbT9mLv6Vm+3p6j/sPMk6/vNlr98lIrZ8cg2wgKz0lij+Bao+MMfOaOZ4CLI/zcJSKibFDs0lyOSGJAMs+2U1BoTvQvSIxYxKxOfDKo9R8uZasy6bT49RL43ZvAwzLddbgjh6cxdyJI7V8Ls4i3ihEKdWyYLcKr8/G5bJY8vEummvXKbafslUPU/3vD8m65FYyRn272esKMHJAJilJVn0vZXBqkN3VukdEF/AjY0ypiATLK4LB8QogL3HDUip6Rg/O4rrReTy7bpezI54tLN9Q5JSXGUP4bRliR0Q4+OZjVG1eQ+b536Hn6AlxvX9SoG99aNnekunnaPlcHGmArFSUjB6cxezx+by6aR9llXVs3V8Z9jwRmwOv/oHqre/Q68If0PPMia1ee/PeCp65ZSzLNxTxXEERfr/unNSFXAWsMcYcEJFnAsHx8zjB8SWJHZpS0TN5VB4rNhQ5SQbLsOyT3Q22k46n8nf/weF/vUzPsyaRec7UuN9/RP+ebNtfyeyVm/DbQmqyU7Z3x8Wx75yhHBogKxUlBYVu5q7aTK23+W7HIsLBN/6Pqk2ryTz3RjLHXhfRtZd9sptJo/K4/9pTmTwqT7MIXYiIfGqMmQisMsbUAjcDxwGXiIj28VOdRmiSoc7rZ/1Od+vfFAMVa5dxaO0yepx+Jb0uuplEdLE9Z2g2s1duql934vHaWnccZxogKxUFBYVuHl39b+paCY7dbz3O4Y2v0XPsFDLPDdfaNjyv33ncGFyAp5Nk1yIi7xljbsTJHG/GCY4PJnhYSkVVMMnQUtefWKvc8DLl7/yd9BEX0vvy/4h7cGyA6Rc4m4GErjuxLKNPDONMA2SljkJoXdibm/fz+Hs7EGm+Sk5EKH/7b1QWvETGmIn0uuB7bZp4BXh2/S56pibxq2+dEpWfQbVfxpgXm3mrDKgCngr++xGRq+M1LqViqbWuP7F2ePMaDr75J9JOPIs+3/oZxnLFfQwLrj2VG88eREGh29mB1WdjGcPciSM1MRJnEQXIxpgXgCeBV0Sktd1ylerUQtu5uSwTUY1cxXtPc2j9Cnp849tkXXLLUWUlRODP7+7gq7Iqbr/wBJ0sO7fmSidej+solIqjrPQULGMQkbj3Pa7+90ccePkRUgedRt+Jv8K44p8/vHxELsP7HelcEWwVquV0iRHpv4AqYAlQYYx5CvibiHwZs1Ep1Y6FZjnsCILj8g8XU/HREnqcdjm9x00/5kd2b24p5r0vS+sb6avOR0RuTvQYlIqXgkI3yzcUsfST3Q16vcdLzc6NlL74a1L6DSNn0kxMUvzbZ7osw+qtxbwbMrdrOV1iRdRdUERuAvoD84DLgG3GmHeNMd8zxqTFcoBKtTeh/TpbU7FuBRXvPU33/IvpfeWdGBOdhp5en7NgQymlOrLgE7lF63bhS0DHitqirZSumEdy7zxyptyLlZoe1/tbxqk79ttS30Nf5/b2IeJPaxE5JCJ/EpGzgFOBAuAxYJ8x5jFjjBZGqi5h9OAsTs/LbPW8QwUvUf72X0k/+Xyyv/XTiIJjV6ARPDgTp8sEWoDi/G3hTKYuXbDRqRljVhtjzovgvF7GmBnGmP+Kx7iUipaCQjcL12xn+YaiFhc3x5KneAclz92Lq0dvcq+fhystI673dwXm9uCvBQa0fWc70uYiG2PMAGAiMB7w4ezkNBD4zBhzt4g8FN0hKtW+FBS6KShsuf1Q5cbXcK9+jLSTzqHP+J9HvNjjpJwejD99AFnpKbirPewtr+HZ9U7jfCNgLIPYciRqVp3V/wOeNcZUA6uAT4C9QC2QBYwAzgOuBF4AfpGgcSoVVuhC5sZlAo3XcSRiTZ73QBHFS2dhpaSTO3UBrh7xLWU4sW93fnjeUOau2lzf93nKmIFMGpWnZRXtRKSL9JJxguIfAuOAfwG/AZ4VkcOBc64G/gFogKw6teUbilpcZX3489UcfP2PpJ1wJn2v/kWbFnts3V/Jd89J4cazBwFHavO8PhtjDHZg8Yrfrz0xOzMR+bsxZhEwBbgBZ+4NPrYQYAvOgr1viMi2xIxSqfBCA+CUJKvJeom2ruOINl9FCcVLZgKQO3UeSZk5cR/DLncNw/tl6EK8dizST+59ONn/RcCvROSzMOe8CySmq7dSMRbMhmSlp7D0k93Nnle15W0OvPJ7ug35Bn2vuRvjSm7zvV7dtK8+QA5dyZyVnlKfbdDHcJ2fiHhx5txFAMaYTCANOBB4T6l2KTQADq2pDQaCwXUcHl/8Syv8h90UL5mB7amh3w0PkJydmN3ag0mOOy4+UQPjdirSAPlnwDIRqW3uBBEpB46PyqiUakdCsyGWMc2usq7a9gFlqx4mddBI+k6acdQroa8a2b/B69CVzMP7ZWi2oYsSkQqgItHjCDLGXAn8HnABT4rIgwkekmonxg7NJiXJqv9lPis9pUlG+brReSxatyuu4/LXVFK8ZCb+wwfJuX4+KblD43r/U/pl8FVZFX6/Jjk6gogCZBH5f7EeiFLt1dodB6jz2gjOhh8uA42fClZvX0fZi78hdcBwcibPxkru1qZ7JFmQPyCTqWcOqs8eh6Ntf1R7YIxxAQtxSu6KgI+NMS+KyJbEjky1B417+IbLKE8elcdzBUVxyyLbddWULJuD172HnOvupVte/PsKbCuuZP41p+Ku9miSowOItAZ5DeE3CROcRSPbgb+LyIYojk2phAqWVVTWeOv/8QuQmZbMweojT7hrdhRQ+sIDpOQOddoEpUTe+dAAN5w9iMm6MEN1LGcB20VkB4AxZjHOOhUNkBXQ9Jf5JJeTUcYY3ti8n8oaL2fkZbJ+Z+wrM8XnoWTFfDz7t9P32ntIG3JGzO8Zji2waW8F9197akLur9om0hKLrcCNwH5gfeDYmUA/nBXU5wP/aYy5UkTeivoolYqzxmUVoRoExzs3Uvr8ApKzB5Fz/Tys1O7NXtMAlmU4KacH2/ZXIkCyy2hwrDqi44DQYvwi4OwEjUW1c9v2V+L1OU/h/LbwaVEFnxbFp1pI/D5KVz5I3a7PyB7/c9KHjY3LfZtTVlnHwjXbNYPcAUQaINcCT4nIT0MPGmN+B4iIjDLG/B6YD2iArDq80EeCIuFrjmt3b6J0xTySevUnd+o8XN16tHhNy4BtC1+WHHbKNXA+LLQbheqsjDG3AbcBDBrUfOmQ6rwWrdvFjOc/T0grN7H9lL38CDXb19N73H/QI//iuI+hX0YqJYfrEIEkl+HtbSWs3loctruHal8i3Sjk+zj1Zo09BgS3RH0CpzenUh3e2KHZJFkGQ/jaoro9X1Dy3H24MvqSO20+rvSWNw45MadHfVBsi+CyDC6jTeFV5IwxY4wxU40x3QOvuxtj2tzLPkr24PS/D8oLHGtARB4XkTEiMqZv375xG5xqHxat28WMFxIUHItw8I0/Ub31HXpd+H0yRn07AaOA3MxuzL/mVM4b1oeLh+fg0x3zOoxIJ1cD5ANfNjo+IvAegAdIzHY4SkVRQaGbB1/ditcv4YPj/dspXjYHV3ovcqctwNW99QzAZSfn8JS7un5V9+zx+bpQQ0XEGJMLrMSp+xVgGLADeBjn6d5PEjCsj4FhxpjjcQLjaThleErV929f8vFumnkAF1MiQvnbf+Pwp6/Rc+wUMsdOif8gAjbtPcTWFzfhs4Uky5DksrSLRQcRaYD8d+AvxphhOBMjODXIvwSeCry+ENgU1dEpFWcFhW6mPvYhzS2s9pTsoGTJLKxuPci9YQFJGa1PcBaQkZasDeHV0XoEKAaygdC+WMuAPyRiQCLiM8bcibNZiQv4q4hsTsRYVPsSXL8R7PyTCBUfLeHQ+hVkjPo2vS74XoJG4fDbgs2RkrqpZw3kuF5p+jnQAUQaIP83zgT9M5yFeeAs2PstR3bOex14NaqjUyqOFq3bxf++9e/mg+OyXRQvnolJTiV32gKSera8+5LBqTsOZgq0RZs6SpcCl4qI2zRcMPoVkLDCXhF5BXglUfdX7VNw/UaiguNDn7xIxXtP0z3/YrIum06j/2fizuAsxvbbQnKSpYuyO5BWA+RAjdttwD9E5EFjTE8AETkUep6IxLfjt1JHIdi6rfFv79/7yzre/bKs2e/zHtxDyeIZGMtF7rQFJPfq1+y5AD1SXNzz7RFaRqGiIQ2nhK2xvjglFkq1G8H1G54EbCF9+PPVuN96nLRhY8n+1k8xJtJlVrFjGbj36pH6WdABtfqvR0R8OJni5MDrQ42D43gwxlxpjNlmjNlujPlVvO+vOr7go7/fvbGNm55cS0Gh03/zwVe2thwcl++n+Nl7ELHJmbaA5N7HtXqvwx4/c1dt1glRRcO7wA9CXktgo45fol2DVHuUgKxt1bYPOPDq/9Jt8Bn0vfqXGMsV9zH07ZHCiTk96NOj4S6q7mqPbindAUVaYrEWGA0UxnAszdJdm1Q0NN7NafmGItbuOMALG5ssvq/nO1TiBMc+D7k33E9Kn8ifaAdXKeukqI7RL4B3jDFnAqnA73AWTWcC5yZyYEo1tmJD/HbHC6rZUUDZi78ldcBw+k6aiUlKjuv9g/r0SGX86QPISk9h7qrN9YuydTFexxRpgPwE8JAxZhBQAFSFvhmHHfR01yZ1zMYOzSYlydnNyWUZnisowue3m11l7asso/jZGdh1VeROW0BKzvFtup9OjCoaRGSLMeZU4D+AOqAbzgK9hSKyL6GDUypEQaGbJZ/sbv3EKKrdvYnS5+8nuc8gcq6bg5XSLa73D/VFcSXb3thGinYq6hQiDZAXBf5+OMx7grOKOZZa3bVJG9Kr1owenFXfSWJveQ3Prt+F3Uxw7K9yU7x4Jv7qcnKnzie134ltutflI3KZfuEJOjGqqBCR/cCcRI9DqZas3XEAXxxrj+v2b6fkubm4evYl9/q5WK1s1hRz4vS69frs+rIK1XFFGiC3LXWWACLyOPA4wJgxYxK1gFa1c8FOEovWtRAcV1c4wXFlKTnXzyV1wPA23eP2C4byq2+dEoXRKgWBdmrlIvJ0o+PfAXqKyP8lZmRKHVFQ6GZPeQ3GEJfex96y3ZQsnY3Vrbuzk2n3XrG/aRgWcGpeJucMzeapj3ZqWUUnElGALCIJqT0OEdGuTUpFatPeirDH/TWVFC+Zha98H32vm0O3vPw2XfeaMwZocKyi7afAj8Ic3wn8DdAAWSVMQaGbx975ire2FhOv5LG3fD/FS2aAZZE7dT5JPeO/S6MBXJZh7sSR3Hi289R6XH4/7XXfiUS8Takx5irgDmAocIWI7DbG3AJ8LSKxXkmtuzapqAq3xtquq6Jk2Wy8B3aRM2kWaYNPb/U6vdOTOVjtrb/msNyM6A5UKSchEC5JURR4T6m4CW2VCXDD4x/FtaWbr/IAJUtmOgunb3wwoq5CsXByvwzmX3sqAAvXbNde951QRAGyMeYm4M/AkzhN64NLRF04K6xjGiDrrk0qGkIn9ozUhv/07bpqSpbdi6d4B32vvYe0oaMjumbvHqlUe/36WE3F0n7gDJyMcahRQPP9CZWKskXrdjHrhc+xxdn8YsqYgXjjGBz7qysoWTILf3UFuVPnk9J3SNzu3dgXxZVs21/J3FWb8fhsUpIsnrllrAbInUikGeRfALeKyOJA1jhoLTA3+sNqSndtUsci2APZ47OxLNNgIYntraVk+Vzq9m6jz8Rfkn7i2S1cqaEfnns8w/tl6GM1FUuLgP81xlQBbweOXQw8CjyTqEGprqWg0M2sFz6vL6Pw+IW3t5XEbcc8J4kxB2/5PnKvv6/Na0OiTQRe3bSvQetQbevZuUS6zcww4KMwxw8DPaM3HKViY+2OA9R5nYksNDgWn4fSFQuoK9pCn/E/p/vwyNrKGpzFeMHaM6ViaA7wAc4TtOrAn1eBD4FZCRyX6kLW7jjQpMZ4T3l8NnIMJjE8JV/T95q76TbotLjctyUuy3DVyP6kJFm4jLb17IwizSDvBU6iaR3cBcBXUR2RUjGQlZ7SJNMhPi+lz99P7c5/kf2tn9F9xIURXy85yWJcfr8GmWl9xKaizTh75Q4DbsUJhr8ReGujiHyZsIGpLmfs0GxSXPHfQlr8XkpfeIC63ZvpM+G/ST/xrLjePxyXZZgXWJynTxA7r0gD5MdxHvEFyysGGmPOB34D3BuLgSl1tEJrjYMTlrvag2Wob+0mfh+lL/6amh2f0PuKO+lx6qVtukfwcRqgj9hULAmwERghItuB7Qkej+qiRg/O4tnbzuHP73zF6i3FcSmtENtP2Uu/o3ZHAb2vuLNNSYxoMxypu540Kq9+nteFeZ1XpG3efmOMyQTexNnFaQ3Ojk4PicjCGI5PqYgVFLpZvqGofoe80Izu2KHZJLksPD7bmXRX/Y6aL9eSddl0Ms64stlrXj4il4uG57Dk4118WnSkNZzgZKWH98uo351PH7GpaBMRMcZsA/qiwbFKsNGDszhjYC9WbymO+b1EbA689keqt71P1sU/bHGejrUkl+H6MQOZHBIYq84v4jZvIjLDGLMAGIFTu7xFRA7HbGRKtUGw1KHOa9dnNuq8No+98xWnD+xFVnoKtu0ExwdeeZTqL96j10U/pOfoCc1e86whWTz+vTEADO+XwdTHPsIXSEEbnKx06O58+ohNxcgvgIeMMXcAn4rEYxsGpRyNn8iNHZqNZYHfjt09RQT3P/9C1edvkvnNafQ8a1LsbtaCIdnpfPPEPhoYd1ERB8gAIlINfBKjsSh11NbuOOBkh0OOCfDGlmLe2FKMAWyxOfj6Qqo2r6HX+d8l8+yWJ93QnsajB2cxd+JIZq/chN8WXJYhKz2l/j2dPFUMLcV5clcA+IwxdaFviogulFYx0XiNxQ/OGcLqrcUxDY4BKj5YROUnK8kYPYHM826K7c1aMLB3OvcHeh2rrifSPsjdgJ/g9EDOoVH3CxFJ/JJS1aWNHZpNSpKFJ5BBbpxis0U4+OZjHP7sDTK/OY3Mb05t8XrGONcoKHTXB7/D+2Vwyck5vPVFCbYIc1dtZni/DA2OVazdmegBqK4lmDXeU15Tv8aizmvz53d3xPzeh9Y/T8UHz9J95GVkXXorxoTb1ik+8vvr755dWaQZ5P8DrgWW4bQW0kd8ql0ZPTiL2ePz6zO8oZzHdU9y+F8v0/OsSRFlJFwGFq/fxYoNRTxzy1iAJiUcuihPxYOI/D3RY1BdR2jWOMllkWQZfLZgx+FTv/LT13Gv+Qvpw88l+6r/wmnikhgGyEhLbvU81XlFGiBfA0wRkdWxHIxSLQnXnSKUu9qDLdKwzEKE8nf/Xv+4rtdFN7eakThzSBYFhe4GnSmABiUcBu17qeLHGJMK3ISzBkSAzcCzIlLX4jcq1UbBcjVbwO+3mXbWIDbtqWiwSDkWqra+y8HX/ki340fTZ8J/Yyw9N84hAAAgAElEQVRXTO/XEstAis7vXV6kAXI1sDuWA1GqJZH0G/6yuLJJlqPig2c5tPY5epxxJVmX3tZicGwMTD9/KOPy+zn38toYY5p0q3BZTVv9KBUrxpgRwGs4mzJ9Hjh8K3CfMeZKEdmasMGpTidYrhbszJM/IJNlBUUxvWf1Vx9Ttup3pOaNoO+1d2Ncicvc3n7BUDLSknXRtYp4J73fAHeZRBYDqS4tNKsRmtUNWrRuFy9s3NvgWMXaZVR8sIjup15G78v/s9XMsQhU1vnqyzUsy9TXGgM8c8tY7rp8OM/edg4Lrj1VJ08VL78H/gUMEpHzReR8YBDwKc5200pFTbAzz12XD+eZW8ayeW8FHl/sVuXV7vqcshceIKXvEHKum42V3C1m92rN5SNyGZffL2H3V+1LpBnkccD5wJXGmC2AN/RNEbk62gNTKlTjrEZWegoL12yv/y3/1U37Gpx/6OMXKH/n76SPuJDsKyOvZQsmoIPlGqEB+R0Xn6hBsUqEc4EzReRQ8ICIHDLGzADWJm5YqjMK9pM3wLb9lSz5eFfM7lW379+ULJ9LUmYuOdfPxUrtHrN7tcQAqckWFw3P0Z1RVb1IA+Qy4PlYDkSpcELrjoP9hitrvMxeuQlbpH4Sy+/fk/e+LAOgcsPLuP/5JOnDz6XPt++KuJYtyYLJo/KApgG51qKpBKoFeoU5nhl4T6moKCh0M/WxDwkmjEN3H402T+lOSpbOwZXWk5yp83ClZ8bmRq0wwHnD+vDTy04K+6RSA+SuK9Kd9G6O9UCUaixc3fHYodlc/9hH9Z0qar02P3rqYypqnIcalZ++zsE3/0TaiWfTZ8L/RBQcn5jTg7OO792gGbxuAKLakZeAJ4wxt3IkY3wO8BjwYsJGpTqNYCLi093lhFZTxCo49rr3UrJkFiYpmZxpC0jK6BObG0XAsgxXjexfP8drYkQFtRggG2NOAr5sbucmY0wycL6I/DMWg1NdW+Pf5pdvKGLznoombdzKA8Hx4U3/rF8F3XfirzCuyB6QnH18bxaEaQavG4CoduInwN+B9wB/4JiFExz/NFGDUp1DaCIiHouMfIfKKF48E7H95N74IMm94lPz2ycjhbJKT5PjfrthT3tNjKig1iKIrUB/oATAGLMLJyAuDLzfG3gTSFw/FtVphZY5GMuw5OPdTYLjoKqt73LglUfpNvg0+l57Dyap+VXQfXqkUFHjxecXkpMsJgXKKpRqj0SkHJhojDkROCVweKuIbE/gsFQnEZqIiDV/dQXFS2Zi11aSe8MDpPQZFPubBhyq8XHBsD68GyjFCxVaTqGJERXUWoDc+BfKLJoGw9rZQsXMpFF5lFXW8dbWYvzNTODV//6QspceIvW4U+g7aRZWcmqL1zxU4+W+q0firvZolkC1W8aYh4FZIlIVOFQLvCwiMd7oV3UlY4dmk2QZPM1NsFFi1x6mZOls/IdKyLl+Lqn9Tozp/Rrz+GxqvX7Gjcjl67Iqeqcns7GoAr9fyylUeJEu0muJ7qqnoq6g0M0Nj3+E1y8tLhSp/upjSlf+hpT+w8i5bg5WSvgWQYYj/1D9tuCu9nDHxfGdoJVqo58ADwLBAHkLcAYQ+/1+VddiQmfI6LM9tZQ8NxdPaSE5k2bSbeDImN2rJet3uus7Vvx6srNDqpZTqOZEI0BWKuqWbyiqz2j4xVlN3bgSvubrf1H6/P2k5Awhd8p9WKnpzV7vzCFZfLanQhdfqI6k8dM5fVqnomrRul08/u5XeGPY51h8XkqfX0Dd3i/oc/UvSDthTMzuFdF40NadKjKtBcgCZBljfCGvexljegde9w7/bUq1XWhLt7LKhjvoNg6Oa3d9TumK+ST3Po6c6+dhdevR4rXrfDazx+drWYVSSuEEx/c8/3nrJx4Dsf2UvvQbanf+i+yrfkL3k8+L6f3C6ZHiotrjR3ACGAs0SaIiEkkN8pZGrz9u9FpLLNQxC11JbVmmyWK80Fe1RVsoee4+kjJzyZ22AFdaRqvX/3xPBduKN2vjd9WRtJagcA6KHIz7yFSHE5qAGD04q8kGIC5Ds+s8joaIzYFXf0/Nvz8i69Jb6XHauOhdvA0Oe/wkuwxTxgxk5IBMTZKoiLUWIF8cl1GoLi90JbXdwixdt3cbJcvm4OrRm5xp81ttLt+vZyollXXa+F11RJEmKLSLkGpRuJ7yqUkNdxeNbnAsuFc/TtWmf5J53k30HDMxehc/Cl6/YIAbz45f1wzV8bUYIIvIO/EaiOragi3d6rx2s48kPMVfUbJ0Nq70THKn3U9Sj9YrfE7L68W7X5Zq7bHqiDRBoaIiNAHh8drMfWkzW/Ydav0bj1L5e09TuWEVGWdeQ+Y3p8XsPm2hj7pVW+kiPdVuTBqVxzv/LmWPu6bJe57SnRQvmYVJTXeC454t77xkcOrMpl94AtMvPEFXKqsORxMUKlqCCQiP18YGPiuqiFnAWLHuOQ59tIQep11O1sU/wpjEry1Nchkma7971UYaIKuEC23pFm7S9h7YTfHimRhXkhMcZ+Y0ey0DjBuRy+kDezUIiDUwVkp1VcEd4h5d/W8+2F4Ws01BKv/1CuVvP0X6yefT+4o7EhYcjxuRy8XDc9i0twKDk3zRzwDVVhogq4QpKHSzYkMRb28rabZJvde9l+LFM8BA7rT7Sc7qH/Y8g9PKMyWQNdbJUCmljhg9OIufXnYSH+88iNdnIwLRbO5WteVtDr7xJ9JOOJM+43+OsRJXGl/r9Wu9sTpmGiCruAqupK6s8fLE+183u3U0gK+imOJnZyB+H7k33E9ydvOPyCwD084aVJ8paLxiWymlOqvG811zr7PSU5g8Ko/Syjo2761gT3ltVO5f/eU6ylY9TOqgkfSZ+CuMK7GhxVUjwydSlGoLDZBV3ARLKSLZ0tR3qJTiZ+9BvDXkTruflL5DWjzfFhjQK63+w6Dxim0NkpVSnVHj+W72+Hzmrtpc//oH5wxpkIyIdm/Wmp0bKV35ICn9TiBn0iys5NQoXr1tDDD9gqGaPVZR0WyAbIz5a6QXEZEfRmc4qrMqKHTzy+WfRRYcHz5I8eIZ+GsqyZ02n5Tcoa1+jwBZ6SmAswtfsBuGtnZTSnVmoR0qvD6bVzftO9Kxwmfz2Hs7Gmy0FM3guG7PVmfDpqwB5LSym2ksGeC8YX24amR/3NUeCgrdOuerY9ZSBrlvo9cX4JQsBbfeGYmzKc27MRiX6kQKCt3c8IST4WiNv6qcksUz8B8+SM7180jtf1Kz5/ZKT+ZQjRc7sBV1cGJ8rqCo/kPAZRlt7aY6nPaQoDDG3AvcCpQGDt0jIq8E3rsb+BHgB34sIq/HYgyqdcEOFcFWlleN7M+6r506Y2MMdoxW5HlKdlCy7F5c3bPImToPV1rPmNwnEqnJzs8dmjnXJ4fqWDUbIIvIhODXgcmwBrhZRKoCx7oDf+FIwKxUWMEMR2v8NYcoXjITX0UJOVPupVveKc2e67IM08YM5Mn3vwaElECP47U7DuDzO/cywJQxA3WSVB1Re0lQPCIiD4UeMMaMAKYB+cAAYLUx5iQR8cd4LCqMYIeKYM0xQDBlbAy4LPA3mn6PtczCe3APxUtmY1LSyJ02P6Ke9LGQZMHUM521J40z6frkUB2rSGuQfwxcGgyOAUSkyhgzD3gLWBCLwamOK3SRSGWNt9Xz7drDlCyZhffgHnKum0O3Qac2e64BLjk5h6c+2onfFlyWYfb4/PrJMDSbMkl7X6oOqJ0nKCYCi0WkDvjaGLMdOAv4KAFjUThBcnD+W7hmOz7baZlp284Oco0dS3DsO1RC8eKZgJA7dR5JmbnHcLWjY4DT8jKZPSG/QRAcOvfrk0N1rCINkHvgZAq2NDreH0hM0ZFqtxovGsnrldbi+XZdNcVL5+ApLSRn0gzShpzR4vmpyRY5Gal4fE6dsYjgrvYATbMpmkFQnUAiExR3GmO+B3wC/FxE3MBxwNqQc4oCx1Q7EFpyYYzBF8USC3+Vm+LFM7E91fS74X6SswdG7dqRsMyRTaAaB8c696toizRAXg78zRjzPxyZGMcCvwZWxGJgquMKfdRV57XZXlrV7Lm2p5aS5+7FU7ydvhN/RdoJZ7Z4bZeB2ePzGd4vg+UbisJmC0KzKUp1AjFLUBhjVgP9wrw1A/gTMA8n4TgP+B3QpnpnY8xtwG0AgwZpZ4F4CA0Us9JTmL1yU1SCZH/tYYqXzMJ/+AA5188jJfeEKIy2eS4Djdd0uyzDlDEDmdzMxh8696toijRA/g+cyfEpIDlwzIfziO+/oz8s1ZGNHZqNyzLYzeyMF2R76yhdMZe6PV/Q5+pfkH7SOa1e2xZnMZ5mC1QXErMEhYhcFsl5xpgngFWBl3uA0NRhXuBYuOs/DjwOMGbMmFjtbqwCQkvbstJTWLjmS1KTLHyeYysPtz01lCybg/dgETmT59Atb0SURty8/pndKGrUp9lvC8cF2nkqFWsRBcgiUgP8Z2CCDv7a+FXoIz+lQrW2clp8XkqfX0Bt4edkj7+L7iefF9F1jaE+W6zZAtVFJCRBYYzpLyL7Ai+vBTYFvn4RWGSMeRgnsz0MWB+rcajIhJa2AVHbTlp8HkpXzMOz70v6XvMr0o7/RnQu3AJjYMSATPYdqm2wwNCWI+08lYq1tm4Ukhb4szGwQEOpJtbuONDi5Cx+L6UrH6D26w1kX/VjeuRfHPG1bzt/qAbFqktJYILiN8aYM3BKLHYC0wPj2WyMWYpT8uED7tAOFokXWtoWLeL3Ubry19QWfkb2t39G+knfjN7FW7qvwJtbips8gbSgfr2J7paqYi2iANkYkwH8FZiMM1kOA3YYY/4M7BeRe2M2QtUhFBS6WbHB6T/cMzWp2T5CYvspe/G31GxfT+/L/5Mep10e0fWNgennD+VX32q+9ZtSnVxcExQi8t0W3luAdi9KuNAgcezQbCxjsCU6EbKIzYFXHqVm+zp6j7udHiMvjcp1I75/o9cGSEl21pvobqkqHiLNIP8a51HaKOD9kOOrcCbJe6M7LNXehU7MQEQbgYjtp2zVw1T/+0OyLrmVjG98K+L7WUBGWnKr5ynV2WiCQoUTLkicO3Eks1duckrcTMMyi4xUF5V1kSX6RYSDb/yJqi1v0+uC75ExanyMforIJLsaLs5buGa79jxWMRdpgHw1cK2IbDTGhP5itxVofR9g1ak0npgnjcrD21pwLDYHXv0D1VvfodeFP6DnmRMjvl+wrY/2tVRdlCYoVBPhNsa44+ITAXh10z6yu6fwwsa99ee3JTguf+cpDm98lZ5nX0fmOdfHZPwt6dsjhdLDnvrXti0NOlc03j1QPxtULEQaIGcBB8Icz8DZalR1IY0nZgO4XAZf4548ASLCwdf/j6pNq8k87yYyx17X7LVTkyxu/uYQnvpoJ16fjctlcd3ovGbb+ijVBWiCQjURLkhctG4Xs1duqt9A6WgcWruMQ+uW0+Mb36LXhd+P8qgjk5mW3CBA9gsNssTaxUjFQ6QB8sc4k/SjgdfBSXo68GG0B6Xat8YTc/6ATCxTRLiiYxHB/dbjHP70NXqOnULmN6e1eG2Pz+apj3Yye3w+7mqPTn5KaYJChRFui+lZL3xe3zv4aHofHyp4ifJ3/0H3ERfRe9ztGHN0QfaxOGtIFtd8I497nj+ySWSSRZMssXYxUrEWaYB8D/C6MSY/8D13Bb4+C7ggVoNT7dfkUXmUVtbRJyOVTXsrwtYfiwjlb/+NyoKXyDjzGnpd8L1WJ1zByUq7qz31jwuV6uI0QdGFhevWEHosOE/e9o9Pmmys0RaHN72Fe/VjpA0bS/a3fooxVjSG3yKDswA7GMsnuwy/vOqU+p9zyce7yOnZjdsvPEGDYRV3kfZB/tAY802cnptfAZcCG4BzROTzFr9ZdSrB+uM6r93iJiAAFe89zaH1K8gY9W2yLv5RxNkIyzJaU6bUEZqg6KLCLcQDwh57c2vxUd+netuHHHjl93QbfDp9r/4FxtXWDrBHb/41p7JpbwUGyB+QydodzsOSG88exI1n6+6LKnEi/r8gEAgnpiBJtRvB+uPWguPyDxdT8dESepx2OVmXTW/To7rWNhlRqivRBEXXFW4hHhD22NF2d6v5egOlL/6G1P4n0XfSTExS/DbiEJwFhT+97CSgaeCvWWOVSJH2QfYD/UWkpNHxbKBERFyxGJxqf8YOzW49c7xuORXvPU33/IvpfeWdrT6qs6BBSyIBbdujVAhNUHRNY4dmk2QZvH5n0V3wyVrjxXnb9lce1fVri7ZQumIByX0G0nfKvVgpacc8ZhP403JfoyPe/7KMj3ceZPKoPG3dptqVSIuMmkv/pQKeZt5TndC2/ZUtZioOffIi5W//jfSTz4+4ju2Gswcx/5pTSbIMlnEmfy2xUMphjPEbY3LCHM8OJC9UZxZ8+hb4O7g4767Lh9dnWd/eVtLCBcLzFH9FybJ7cfXsQ+71c3F16xG1IV99xoBmgwZwao1Pz8us30/KG3gqmZJk4TLa1lO1Dy1mkI0xdwW+FOB2Y8zhkLddwPnAFzEam2ongrvklVTW8XlRebPnVW58Ffdbj5N20jn0Gf9zjBXZg4X8AZncePYghvfL0LY9SjWlCYouau2OA/j8TvDo9x8ppwidJxet28UbW9pWf+w9sJviJbOwUruTO3Ueru7Rm28FePHTvXRLtqjxHskjJ1mGYTk9GDU4i0mj8gC44fGP6rPjk0c57Tz1M0C1F62VWPxX4G8D3ELDlkIeYCdwe/SHpdqLgkJ3RLvkHf58NQdfX0jaCWe2aZGHAdzVzme8tu1R6ghNUKjGLTWz0lMa1On+4Jwh/PndHW26pq+imOLFM8GyyJ02n6SeTR5OHDNbaBAcA/ht4esDVcy/9lRGD86ioNCNQP0f0M8A1b60GMWIyPEAxpg1wCQRccdiEMaY3wITcILur4CbRaQ88N7dwI9wgvMfi8jrsRiDCm/tjgOt7pJXteVtZwX0kG/Q95q7Ma7It4TWR2lKNUsTFF1csJxixYYiBHh7W0l9ByGP1+axtgbHhw9SvHgm4q0l98YHSe59XEzGHY7gLC4M1hav2FCEN9CXzusXVmwo0uBYtSuRdrG4kjCP+Ywx3QBbRI71Md+bwN0i4jPG/Bq4G/ilMWYEMA3Ix9lqdbUx5iQR0bq7OMlKT8GyDP5mOktUffE+ZaseJnXQSPpOmtHiCmjLgDENr3XdaN0hT6lw4pWgUO3f8g1FTVprhvYPjoS/5hAlS2bir3KTO3U+KTnHR32crbHF+UyBpttKae8i1d5EukhvKfCfYY7fHnjvmIjIGyLiC7xcC+QFvp4ILBaROhH5GtiO0/tTxUFBoZu5qzZj2xK2CLL6y3WUvfRbUgecTM7k2VjJ3Vq8nmXg1vOOdxbjAd2SLSaPymvxe5RSXAnUND5ojOlmjIlfTy4VUwWFbhau2U5BYcPfg9buONA0OAYuPSWXSLtn2nXVlCybg9e9j76TZ5F63MlRG3fomMKNp/GhzXsr6n/GJJfBACkuo58Fqt2JNIN8LjAjzPE3cZrYR9MPgSWBr4/DCZiDigLHmjDG3AbcBjBokDYXj4ZwE3NQzY4CSlc+QEruCeRE2B5IBDLSklky/RxdiKFU5JYC7wAPNzp+O3ARcE28B6SiK9yGIMG5sbLG22QOTk22GNqne0S9j21vHSXL5+LZ/xV9r51B2uDTo/8DABPPGEC/nt0a1ESf0i+DHt2SKCh012e7l3yym2Wf7MZnC0mWYerZg5g0Sp8kqvYn0gA5HfCFOW4DGZFcwBizGugX5q0ZIrIycM6MwH2eiXBc9UTkceBxgDFjxujTmmNUUOjmnW0l4YPjnRspfX4BKX0Gk3P9XKzU9FavZ3Gk3lgXYijVJvFMUKgECN0QpM5r19fjFhS6efL9rxucawycP6wvyzcUtXpd8XspW/kgdbs302fCz0kfdnasfgSG5WYwdmg2O8qq2FFWxc4DVXyxv7LJZ4jfL/hxSir8tjCgV5p+Hqh2KdIA+TPgBmBOo+M3ApsiuYCIXNbS+8aYHwDjgUtF6n8v3gMMDDktL3BMxVBLnStqd2+idMU8knr1JyfC3pkGOHdYH3562Uk6ESrVdsecoFDt29ih2SS5rPpdShev30X+gEzc1Z4m6z9E4M0I2rqJ7ads1cPUfPUxva+4g+4jLorN4HFauIV22LCMwRYJm2BJTrJABL8tukhbtWuRBshzgZXGmBOBfwaOXQpMAa491kEYY64EfgFcKCLVIW+9CCwyxjyMs0hvGLD+WO+nWtZc54q6PV9Q8tx9uDL6kjttPq70zAbvp6e4qPX6sSWkHk0gJdnS4Fipo3fMCQrVvo0enMV1o/NYtG4XAH6BWSs3ccnJOSS5TH23h0iJCAdfX0j1F+/R66KbyTjjqlgMG3CC47kTR+Ku9tRnwREJW498el4msyfkA2iZnWr3IgqQReQVY8wEYCbwv4HD/wKuFpFXozCOP+I0vX/TOP9XrRWR20VkszFmKbAFJ4Nyh3awiL2xQ7NJTrIaZJDr9m+neNkcknpkkTNtQdjG8tUePykuw5QxA+sbweskqNQxi2mCQrUPk0flsfRjpzYXnPKD1VuKsSJciBckIrjX/IXDn71Bz3Omknn25KMeU3Cnu5beH39af9zVHrLSUxr0bL5gWN8GG5i4LMPsCfn1nwX6maDau0gzyIjIa8BrsRiEiJzYwnsLgAWxuK8Kb/TgLJ69dSxzX9rMp0UVeEp2ULJkJla3HuRMXUBSRvOPxBrXlOkkqNSxiUOCQrUDowdnMXfiSGav3ITflvoNNGxxOgBF2tKt4sPFVH78AhmjJ9Dr/O8c05hau6UAL2zci2WcbaJnj8/HXe2pL5t498tS6rw2JtDBSD8PVEcSaZs31cWMHpzF1DMH4Tuwi+LFMzHJ3cidtoCknn2b/R6DbvyhVCyIyGsicp6IdA/8OU+D487nxrMHMXfiSE7LyyTJcj6gAxULETn0yUoq3n+G7iMvJevSWzGR9oE7RraA12fjrvZwx8Un1i/Enj0+H1cgBf7URzubtLBTqj1rNoNsjDkEDBWRMmNMJS38MikiPWMxOJU4BYVuZv79dfY/OwNjucidtoDkXkeakNx+wVB2lFXx1tZiBKcOLVhaoVkCpZRqu2DveY/PJsllkT8gg8+KKiLaROPwZ2/gfusJ0k/6JtlX/Rhj4pf/ai454q72YIvUB9DBXfSU6ghaKrH4L6Ay8PWdcRiLakdWvf8vdj19NyI2uTc80GBLUoPTz/jx742hoNCtdcZKRZkmKLqm0HZvfr9N/nGZbNp7qNmdTIOqvnifA6/9kW5DvkGfCf+DsVzHPBYDnJaXyaY9FTReI9gj1cXhOn/9eec106Vo7NDsBnXJ+nRRdSTNBsgi8vdwX6vOb9UHn/HIz76L+Dz0v+F+uuUMrl84As5ii+BEpz2NlYoJTVB0QcF2b16fjbEMG3e5Ww2Oa776hLKXHiJ1wMn0vXYGJik5KmMRIKdnN+yiiibvfefswTz10c76wLe5LkWjB2fxzC1jNYmiOqSIF+mpruHVdZuZNOFK/NWHOO6mB/jO+AsZOSDTeezntbECLX10olMqdjRB0TVt21+JL9AL2ecXtuyrbPH82t2bKH3hflL6DiZnyhyslG5RG4sJ/HFZpkGC5KwhWfzqW6cwLr9fRIGvJlFUR9VSDbJN64tYARCRY3+eoxJu//793DzlanxV5eROnU9yzlCO65XGjWcPYni/DM0CKKVUjBQUupm1chNNO9CHV7fvS6cvfc+cwI6m3Y/63n16pDC0T3e2lxymvMaLLc6H/9vbSrjlvON54v2v6zPZG3eXU1Do1sBXdXotZZCv50iAnIvTi/N54KPAsXOAa2javF51QGVlZVx22WVUlO0nb9pckgcMb1AzppOhUvGjCYquZ/mGolbLKYI8pYWULJuDldaT3KlNN21qK3e1h/U7PYDTOSPY/9hvCxlpyUw9cyDPrtuFAF6/sDywFbZSnVlLNcjPBb82xrwI3C0iT4Sc8ldjzHqcIPn/YjdEFWsHDx5k3LhxfPXVV7z88stknvANzRYrlViaoOhCCgrdLP1kd7Pvp7gMnsBKOW/5fkqWzsK4ksidOp+knn2O6d6De6dTePDIBrY24Ap0hwtNkjz3yW48fqc/83MFRUzWjkWqk4u0BvkS4K4wx9cAj0ZvOCreKioquOKKK9iyZQsvvvgil1xyCaAbfCiVSJqg6FrW7jiAv4XtpIPBsa+yjOLFMxCfl9wbHyQ5q/8x3zs0OA4S4NJTcrn9whPqPwumjBnIokAW2e/Xlm2q84u0UWIZcF2Y49cBpdEbjoqnyspKrrrqKjZu3Mhzzz3HFVdckeghKaWaugQnGdHYGuCi+A5FxcLYodm4XA039Wi8xbS/uoKSJbOwaw6Rc/19pPQdHLPx2AJrvihpcGzSqDxSky1cRjeEUl1DpBnk2cDfjDEXc+QR31jgMuBHsRiYiq3q6momTJjA+vXrWbp0KRMmTEj0kJRS4QUTFA82Oq4Jig4utI/8xcNzeHNLcf17/TO7sae8FgC7roqSpbPxVRSTM+U+UvufFJPxhG5pbYs0yBJryzbV1UQUIIvIP4wx24AfA1cHDm8FzhWRdbEanIqN2tpaJk6cyHvvvcczzzzDpEmTEj0kpVTzNEHRCRUUurnpybX1u+b1TGv4cbw3GBx7ayl57j48pTvJmTSTboNOjcl4LAOXnZLLP78owRYhJUyWWBdrq64k4j7IgUD4phiORcVBXV0dkydP5q233uJvf/sb06ZNS/SQlFIt0ARF5xS6a57HZ1NW6WnwvgDi81L6/P3U7fmCPhP+h7QTzjyqexnAGEiyDBcNz2HNthJ8fiHJZTA43SqSkyymX3gC0y88QbPEStGGANkYkwt8FxgKzA5sgXousFdEvo7VAKIJFVQAACAASURBVFX0eL1epk6dyiuvvMJjjz3G97///UQPSSkVAU1QdD5Z6SkYY0DCL84T20/ZS7+l9usNDL7mZzD8/CbnjOifwdZ9la32Axw3IpfTB/aqD3qDpR1Z6Sls3luBQIOuFBoYKxVhgGyMGQ28BXwN5AMP4dTFjQNOAm6M1QBVdPh8Pm666SZWrlzJH/7wB2677bZED0kpFSFNUHQui9btYtYLn9Nc4woRmwOv/oHqf39I70tuYdy1NzaoTw76ztgh7DpQxWPv7mg2SE5xGaaHdKMI2lNew+/f+hKf3yYlyWLyqLxj/KmU6lwizSA/BPxeROYYY0L3vnwduDn6w1LR5Pf7ufnmm1m2bBkPPfQQd955Z6KHpJSKkCYoOpeCQjezV25qITgW3G89QdWm1WSeewM9z7wmbHAMzk53b31REjY4dhmnVVswOA7NGs9dtZk6r13/fV6ftm1TqrFIA+TRhF8Msg+nib1qp2zbZvr06Tz99NPMnz+fn//854keklKqbTRB0Yms3XGgxR3zKt5/hsqCl8gYM5HMc29sNjNsGXhzS3Gz79sC735ZyvQLT6Cg0M0Nj3+E1y9YJlDfHDjPoG3blAon0gC5Bgj3q+XJQEmY46odEBHuvPNO/vKXvzBr1ixmzJiR6CEppdpOExQdWGgrN4C95TW4LPDZTc+tWLeCig8X0+O0y8m65BanRrkZoUFuc+8HM8N7ymvqNxvxC7gsgwvBZRmmjBnIJN0VT6kmIg2QVwJzjDFTAq/FGDME+DWwPAbjUsdIRLjrrrv405/+xC9+8Qvuu+++RA9JKXV0NEHRQYVmbl0WWJaFz28HAt+G4W3lxtcof/uvpJ98PtlX3IFlDGFiaMDpRhGahW56NSfDHMwMr9hQ1OC90YN6ceHwHO1UoVQLIg2Q/xt4BacpfTrwPk7m4gNgZmyGpo6WiHDPPffw6KOP8pOf/IQHH3ywxUyEUqpdi2mCInDde4FTgLNE5JOQ9+7GyV77gR+LyOuB41cCvwdcwJMi0ngTEwUs31B0ZJtoG7CdkNdCSLIMti3YQNWWdzj4+kK6DR1Nzvi7OG94LoN6p/Ps+l2EVmO4LJh25iDyB2Qyd9VmPD4byxjmThzJ8H4ZrNhQhAAjB2TirvY0CICXfLIbX2AsG4sq+OVVp2hwrFQLIg2QfThbml4AjMLZonqDiKyO0bjUMbjvvvt48MEHuf3223nkkUc0OFaqY4t1gmITMAl4LPSgMWYEMA1nYeAAYLUxJriF20KcRYJFwMfGmBdFZEsUxtKpNJ55XZbT1i05yeLK/H5s3F1Oz9LPePmVR0gbmE/utXeT2i2Vn17m/GdevqEIr8/G5bK4bnReg1Zsw/tlNOlX3FzAO3pwFlPHDGTRul0I4PfrojylWtNqgGyMcQEVwOki8k/gnzEflTpqDzzwAPfddx8333wzCxcu1OBYqY4vpgkKEdkKhJsrJgKLRaQO+NoYsx04K/DedhHZEfi+xYFzNUBuZNKoPJYVOEFucpLFvRPycVd7qKzx8ud3d1Bb+BnFy+bQZ9AwHnhyMXUmtUHAO3t8Pq9u2sdVI/vXB8RwZEe7YHeKhWu2t1ouMWlUXn3ArYvylGpdqwGyiPiNMYVAShzGo47BI488wj333MNNN93EE088gWVZiR6SUuoYJDhBcRywNuR1UeAYwO5Gx8+O16A6ktGDs3j21rFNMr3f/cs66vZuo2TFPJKz+tNtwiweWrOLZ24ZW39OQaG7voxi3Y4DYEx9z+LgeaHbVTc+Hi67/MwtY+vLMJRSLYu0xGIe8KAx5jsiUhbLAamjs3DhQu666y6uu+46nnrqKVwuV6KHpJQ6RtFKUBhjVgP9wrw1Q0RWHsu1W7nvbcBtAIMGDYrVbdq1YIAazP4CVO3bQcmyObjSM8mZOh9XWs8mvYhDt6L2+p2eFaGdKUYPzmp4TuA4EDZoDlq+oQiPz2bFhqIm7ymljmjLIr3jgT3GmCKgKvRNETkt2gNTkXvyySe58847ufrqq1m0aBFJSRHvIK6Uav+OOUEhIpcdxbftAQaGvM4LHKOF443v+zjwOMCYMWO6ZOIyNMubZBk87r0U/eN/MMmpnPWfj1AsGYg07UU8dmg2KUkWHp+NAazAor7Q84LnhJZNhAuawwbdujmIUi2KNJJaTsstF1WC/OMf/+C2227jyiuvZOnSpSQnJyd6SEqp6EpUguJFYJEx5mGcRXrDgPU4a8+GGWOOxwmMp6G7+TURLHPYU15TH5TWuEvY98w9iG0z4IYHuHBMPuD8B23ci3j04Cxmj893dt2zBRcw7axBDc4Llk00LqdoHDQHhQuolVLhRRQgi8i9MR6HOgpLlizh5ptv5pJLLmHFihWkpqYmekhKqeiLaYLCGHMt8AegL/CyMWajiFwhIpuNMUtxFt/5gDtExB/4njtxdvJzAX8Vkc2xGl9HU1DoZvmGIp4rKMLnt0lyWU7muNJN8ZKZSF0V/W+4n7ScQfXnpCRZTBqV1+Ra7moPtjilFX5bGNArrUnGN7hYL/R1uKC5tfeUUg21GCAbY9KB3wLXAMnAapxemFqHnGDPP/88N910E+eeey4rV64kLS0t0UNSSsVArBMUIvI88Hwz7y0AFoQ5/gpO6zkVYtG6XcxeuQlfSPNiv99m4oheLL3v51jVB7n79/+PHUl5pCW7WL21GFvA00y5w9FmfBsHzZG+p5Q6orUM8n3AD4BngFrgBuBPwJQWvkfF2Msvv8zUqVM588wzefnll+nevXuih6SUijJNUHQsBYXuJsExgPHX8dajP6Nk11c8/OQi/vhlOh5fGS7L1G8CYgtkpTddh9lSxjdcpwqlVPS0FiBPAn4kIosBjDFPAx8YY1zBR20qvt544w0mT57MaaedxquvvkpGRkaih6SUig1NUHQga3ccwJaGwbH4POxdPo/aXZ/x4MK/YgaegWfrNmwB8R8518IppwgnXMa3ufZuSqnoaa1R7kDgveALEVmPU4s2IJaDUuG9/fbbXHPNNQwfPpw33niDXr16JXpISqnYCSYobhORHwPfBq4J9EZWMRDcdKOg0N3m7w2WQ1gGLANi+yl98TfU7NxI9pU/4a97cshKTyElycJlnK4VKS6Dy0BKctsWzDXX3k0pFT2tZZBdQONfa30RfJ+Ksg8++IDx48dz/PHHs3r1anr37p3oISmlYqtJgsL8//buPN6qqv7/+OuNgAOKYaCkCIJTipYKX7/4K02ZRCwQU3PIITHCrDQ1S4vCkl9W9rVoIHHI4UtSgRqoKcgV8/tVNHECJwTCZBIFJwSRYX3/WOvi5nDOne/d517fz8djP+49e/zsdc5d93PWXnttqbKB4tWSW1md1LdVNtsdYuftWnPB14az9uVZdOj/Ndod3I/1Gzbx5poPij4dr7bdJDwahVnjqy7RFfDfktZl5m0HXC9pTeWMEMKQxgjOoscff5zjjjuOPfbYgxkzZtCpU6e8QzKzxucGiibUEGME9+rWgRACF1/0Ld6ZU8FRp3+DZd0Hs3FjTGQ77NB289Px/rloFT/8fE+WvrWWyU8u3rx9KYV9jj0ahVnjqq6ivaXIvP9ujECsuKeeeopjjz2WTp06MWPGDDp3LvYwLDNrgdxA0YQaolV29itvcuwZI1n1yF/pcMQX+eWYK5G0OZHNJuEfbNjEqDTGMcCkJ17l9hFHFE12S7VuOzE2azxVJsghhK80VSC2tTlz5jBgwADat29PRUUFXbpsPU6mmbVYbqBoQvVpla1s3Z10029Z9chf2fGQQbQ/8hzueGoJ/3/YwVvsqzIJl7TFiBfrN4aSrdZ+Ap5Z0/OlujL14osv0r9/f7bddlsqKiro1q1b3iGZWRNyA0XTq0urbGXr7huPT2XltHHseODn2GXA+SAxafZivljiyXcddmjL6KmxuwVAm220Vat1ZeJdeXOf+xybNR0nyGVo/vz59O3bF0lUVFSw99575x2SmZkVMWvhSlY9PYOV08axwz6Hc/SI0cxd9l58+t3GrVt7s0n4/p134o4nFxNgi0Qatu5W8cPP9+TNNR+4z7FZE3GCXGYWLVpE3759Wb9+PTNnzmT//ffPOyQzMyvh/Zdn8fo917Jdt0+xx0lXcFqfHvz47udq1NpbVYt1YbeKN9d8wAXH7NNYp2FmBZwgl5FXX32Vvn37snr1aioqKujZs2feIZmZWUZ2NIlV857gim8O56BPH8pXxtzI53p2pVe3DnUevi3LQ7mZ5csJcplYtmwZ/fr1Y+XKlcyYMYNDDjkk75DMzCwj2+1h47IXWfGXH7L//vszc8a0Lcamb4gRJjyUm1m+nCCXgRUrVtCvXz+WLl3K9OnT6d27d94hmZlZgcpuD+8vX8hrE39Eu5134Rc3/qXRHtxU6jHTTprNGp8T5JytXLmS/v37s2jRIu677z6OOOKIvEMyM7Mi+vT4OLy9lNf+Mgq13YGdT/wxF075FxN23a1JktX6Pu3PzGquVd4BfJS99dZbDBw4kHnz5jFlyhSOOuqovEMyM7MSOvIOa+4aTdttWtH51KvYZuddN49LXJ3Zr7zJ7x6cz+xX3qzz8YuNh2xmjcMtyDl55513GDRoEHPmzOGuu+6if//+eYdkZmYlLF++nP79+7Pmvfc4ffQNPLyq3eZHSFd3A11Dtfz6xj2zpuMEOQerV6/m+OOPZ/bs2UyaNInBgwfnHZKZmZWwatUqBgwYwJKly+h0yo958PXtad0qcOrhXTmxYPziYhrqSXi+cc+s6ThBbmJr165lyJAhPPLII0ycOJGhQ4fmHZKZmZXw7rvvctxxxzFv3jxGjBnPPSs7sinAxk2B3T+2fY2S1IZs+W2IETLMrHpOkJvQ+++/z7Bhw5g5cya33XYbJ598ct4hmZlZCe+//z5Dhw5l9uzZTJ48mS6HHMX0G2bVOtF1y69Z8+MEuYl88MEHnHLKKdx///3ceOONnHHGGXmHZGZmJaxfv55TTjmFmTNncuutt26+2lfXRNctv2bNixPkJrBhwwZOP/10pk6dyu9//3vOPffcvEMyM7MSNm7cyNlnn83UqVMZN24cX/7ylzcvyya6HpPYrOUqq2HeJF0iKUjqmF5L0lhJ8yU9K+mwvGOsrY0bN3LmmWcyefJkrr32Ws4///y8QzIzsxJCCJx//vncfvvtXH311YwcObLoepUjU/xy2kucccOseg3fZmblp2wSZEl7AgOBf2dmHwfsm6YRwLgcQquzTZs2MXz4cCZOnMjVV1/NRRddlHdIZmZWQgiByy67jOuvv57LL7+c7373uyXX9ZjEZi1b2STIwLXAZUDIzBsK3BqiWcDHJH0il+hqqbIV4pZbbuHKK6+ssqI1M7P8jRkzhmuuuYYLLriAE0dcWuWDPSpHpthGeExisxaoLPogSxoKLAkhPCMpu2gP4NXM68Vp3rIi+xhBbGWma9eujRdsDYQQuPDCCxk/fjxXXHEFo0aNyjUeMzOr2tixYxk1ahRnnnkmZ19yJV++8bEqH+zhkSnMWrYmS5AlPQB0LrLo+8AVxO4VdRZCGA+MB+jdu3eoZvVGU3mJ7je/+Q0XX3wxV111FQVJv5mZlZGbb76ZCy+8kGHDhnHTTTdx3cOLavRgD49MYdZyNVmCHEIo+ixlSQcD3YHK1uMuwJOSDgeWAHtmVu+S5pWtUaNGbb5Ed8011zg5NjMrY5MmTWL48OEMGDCA22+/ndatW/uRzmaWfxeLEMIcYNfK15IWAb1DCG9ImgJ8Q9JE4D+Bt0MIW3WvKBdXXXUVY8aM4bzzzmPs2LFOjs3Mytwdd9xBnz59uPPOO9l2220Bd58wszJIkKtxLzAYmA+sAb6Sbzil/eIXv9jcf+26666jVatyuv/RzMyKufXWW1m7di3t2rXbYr67T5h9tJVdghxC2CvzewAuyC+amhk7diyXXXYZX/rSl7jpppucHJuZNROtW7dmp512yjsMMyszzuTq6brrrtt8c8dtt91G69Zl953DzMzMzGrBCXI9/PGPf2TkyJEcf/zxTJw4kTZt2uQdkpmZmZnVkxPkOvrTn/7E8OHDGThwIJMmTaJt27Z5h2RmZmZmDcAJch1MmjSJs846i8997nPceeedbLfddnmHZGZmZmYNxAlyLU2ZMoXTTjuNPn36MHXqVHbYYYe8QzIzMzOzBuQEuRbuu+8+Tj75ZA477DDuvfdedtxxx7xDMjMzM7MG5gS5hmbMmMGwYcPo2bMn9913H+3bt887JDMzMzNrBE6Qa+Dhhx9myJAh7LvvvkybNo0OHTx4vJmZmVlL5QS5Go8++iiDBw+ma9euTJ8+nY4dO+YdkpmZmZk1IifIVXjiiScYNGgQnTt3ZsaMGey22255h2Rm1qAknSzpOUmbJPXOzN9L0lpJT6fpD5llvSTNkTRf0lhJyid6M7PG4QS5hGeeeYaBAweyyy67UFFRwe677553SGZmjWEucCLwjyLLFoQQDknTyMz8ccBXgX3TNKjxwzQzazpOkIt4/vnn6d+/P+3ataOiooI999wz75DMzBpFCOGFEMJLNV1f0ieA9iGEWSGEANwKnNBoAZqZ5cAJcoF58+bRr18/2rRpQ0VFBd27d887JDOzvHSX9JSkhyQdmebtASzOrLM4zTMzazFa5x1AOVmwYAF9+/Zl48aNPPTQQ+y77755h2RmVm+SHgA6F1n0/RDC30pstgzoGkJYKakXcJeknrU87ghgBEDXrl1rs6mZWa6cICevvPIKffv2Ze3atTz44IMccMABeYdkZtYgQgj967DNOmBd+n22pAXAfsASoEtm1S5pXrF9jAfGA/Tu3TvU5vizX3mTWQtX0qfHx+nVzUNrmlnTcoIMLFmyhL59+/L2229TUVHBpz71qbxDMjPLlaROwKoQwkZJPYg34y0MIayS9I6kPsBjwFnAbxry2LNfeZMzbpjFBxs20bZ1Kyac18dJspk1KfdBBjZs2ED79u25//77Oeyww/IOx8ysyUgaJmkxcARwj6T706KjgGclPQ1MAkaGEFalZV8HbgDmAwuAvzdkTLMWruSDDZvYFGD9hk3MWriyIXdvZlYttyAD3bp1Y/bs2bRq5e8LZvbREkK4E7izyPzJwOQS2zwBHNRYMfXp8XHatm7F+g2baNO6FX16fLyxDmVmVpQT5MTJsZlZeejVrQMTzuvjPshmlhsnyGZmVnZ6devgxNjMcuNmUzMzMzOzDCfIZmZmZmYZTpDNzMzMzDKcIJuZmZmZZThBNjMzMzPLcIJsZmZmZpbhBNnMzMzMLMMJspmZmZlZhhNkMzMzM7MMJ8hmZmZmZhkKIeQdQ4OT9DrwSt5xJB2BN/IOohG19PODln+OPr+m0y2E0CnvIPLQRPVyOb3XNeWYm4ZjbjrNLe6i9XKLTJDLiaQnQgi9846jsbT084OWf44+P2spmuN77ZibhmNuOs017kLuYmFmZmZmluEE2czMzMwswwly4xufdwCNrKWfH7T8c/T5WUvRHN9rx9w0HHPTaa5xb8F9kM3MzMzMMtyCbGZmZmaW4QTZzMzMzCzDCXIjknSJpCCpY3otSWMlzZf0rKTD8o6xriT9QtKL6TzulPSxzLLL0zm+JOnYPOOsD0mD0jnMl/S9vOOpL0l7SnpQ0vOSnpN0YZq/i6Tpkl5OPzvkHWt9SdpG0lOS7k6vu0t6LL2Xf5bUNu8Yre4knZw+w5sk9c7M30vSWklPp+kPmWW9JM1Jn4GxklQOMadlRevMcquDJI2WtCRTvoMzy8q23i+3cixF0qL0GX1a0hNpXlnVz5JukrRC0tzMvKIxNvecxwlyI5G0JzAQ+Hdm9nHAvmkaAYzLIbSGMh04KITwKWAecDmApAOBU4GewCDg95K2yS3KOkox/474nh0InJbOrTnbAFwSQjgQ6ANckM7pe8CMEMK+wIz0urm7EHgh8/pnwLUhhH2AN4HhuURlDWUucCLwjyLLFoQQDknTyMz8ccBX+bAOHtT4YW6haMyl6swyroOuzZTvvVDe9X4Zl2Mpx6SyrfwSVW71881s/bdTKsZmnfM4QW481wKXAdm7IIcCt4ZoFvAxSZ/IJbp6CiFMCyFsSC9nAV3S70OBiSGEdSGEfwHzgcPziLGeDgfmhxAWhhA+ACYSz63ZCiEsCyE8mX5/l5hA7kE8r1vSarcAJ+QTYcOQ1AU4HrghvRbQF5iUVmn25/hRF0J4IYTwUk3XT/Vs+xDCrBDvTL+VJv4MVBFzqTqzOdVB5VzvN6dyLKas6ucQwj+AVQWzS8XYrHMeJ8iNQNJQYEkI4ZmCRXsAr2ZeL07zmrtzgb+n31vKObaU8yhK0l7AocBjwG4hhGVp0XJgt5zCaii/In453ZRefxx4K/OFrkW9l7aV7ql7zUOSjkzz9iC+75XK6TNQqq4p1zroG+ly+U2Zy/3lGiuUd2yFAjBN0mxJI9K85lA/l4qxOZX9VlrnHUBzJekBoHORRd8HriB2r2jWqjrHEMLf0jrfJ166n9CUsVndSdoRmAxcFEJ4J9sVM4QQJDXbsR8lfR5YEUKYLenovOOxuqtJ/VPEMqBrCGGlpF7AXZJ6NlqQBeoYc1mp5n/bOOAnxETuJ8AviQ0k1jA+G0JYImlXYLqkF7MLm0P93BxirCknyHUUQuhfbL6kg4HuwDMp8egCPCnpcGAJsGdm9S5pXlkqdY6VJJ0DfB7oFz4cULtZnWMVWsp5bEFSG2JyPCGEcEea/ZqkT4QQlqXLXyvyi7DePgMMSTcPbQe0B35NvLTXOrUit4j3sqWrrv4psc06YF36fbakBcB+xPe7S2bVRvkM1CVmqq5rmrwOquk5SLoeuDu9LOf6spxj20IIYUn6uULSncTuIc2hfi4VY7Mp+2LcxaKBhRDmhBB2DSHsFULYi3hJ4bAQwnJgCnBWurOzD/B25rJEsyJpEPEy9pAQwprMoinAqZK2ldSd2Dn/8TxirKd/Avsqjn7QlngDypScY6qX1Bf3RuCFEMJ/ZRZNAc5Ov58NNIuWrmJCCJeHELqkv71TgYoQwhnAg8BJabVmfY5WmqROlTeHSepBrH8Wpnr2HUl90t/BWZTPZ6BUnVl2dVBB/9FhxBsPobzr/bIrx2IktZO0U+XvxKvQc2ke9XOpGJt3zhNC8NSIE7AI6Jh+F/Fu2gXAHKB33vHV47zmE/sWPZ2mP2SWfT+d40vAcXnHWo9zHEwcoWMB8RJp7jHV83w+S7w0+mzmfRtM7KM7A3gZeADYJe9YG+h8jwbuTr/3IP7Dng/8Fdg27/g81eu9HUZsfFgHvAbcn+Z/EXgufbafBL6Q2aY3MeFYAPyW9CTZvGNOy4rWmeVWBwG3pf9dzxKTn09Udw7lMJVbOZaIsQfwTJqeq4yz3Opn4HZiV6b16fM8vFSMzT3n8aOmzczMzMwy3MXCzMzMzCzDCbKZmZmZWYYTZDMzMzOzDCfIZmZmZmYZTpDNzMzMzDKcILdwkjpLmibpvZbydJuPEkk/knRT5vVMSb9txOMdLSlI6thYx6ji2OdIWp15fYGkqU0dh1ljcX1cfiT9UdIPM68XSbq0EY+3RT3XlCSNljS3+jWr3c/Bkpak8ZpbLCfIZUrSzSlRCZI2SPq3pHGSOtRyV5cCuwOHAJ+oZt1mKf3RhyLT8rxjq4/0uNFLgKvKIJY+kt6ofAhDE7kB6CXpyCY8ptlWXB/XXBX18Ql5x1YoPfn2BOBXZRDLqZKezunY35N0V8G8EyVtlDShcP0QwhxgFnBxU8WYByfI5e0BYiW6F3Ae8AXg97Xcxz7A7BDCyyE+za/WJLVq4sSoLl4illV2OrgxD9gE5XIe8HgIYWEjHqOmhhIfurGxqQ4Y4mOD/wR8q6mOaVYF18c1V6w+/ntddpSeftdYvglMDiG804jHqKmhwF3VrtV0xz4P+DlwQokvgn8EzpfUurGDy4sT5PK2LoSwPISwOIQwDfgz8fGTm0n6iqTnJb0vaZ6kb0tqlZYtIn7wz0rf4G9O83eWNF7SCknvSnpIUu/MPs+RtFrS4HQ55gPgAEltJf1M0mJJayT9U9Kxme0qL8/3k/RYWucJSYcVxNxHUkW6zPh2+n33tEySLpO0QNJaSXMkfbkGZbUhlVV2ej1zzEWSfiDpOknvpHP4TkFcdS2X3SRNSfG+kt6TuZJGp+1uknR3wbFapVaoqr6Bnw5U2cUglfVbkkZm5p2dym2dpNck3ZJZdrGkZ1PZL5F0g6SPVVmy0VDS40NTa9rdkr4raXl6D69O5zQ6ld9ySd8tiLWrpDtT2b4r6Q5JXao57hRgiKQdahCjWWNyfVy/+nhd2ufBkh5I+1uV6pOdM/Fk65fFxKe1IWl3SRMkrUzn8rSkYzLbfUHS7FT2/5I0RlUk14pfMk6h+jr2y4r/M4ZkyuQSSS+nOnaxpJ9m1r9a0kvp/BZJ+rmk7ao5RhvgOD6sY2cqXqH4ZSqj1yVdqPgo798p1vn/lnRmwX6qLNsSx+5MfMrk3Zl5XYBjgGuILcVnFNl0GrAL8YmlLVPej/LzVHwCbiY9Jje97kF8/OTyzLyvEh/5eBLQndiisRz4RlreCZhOrMg7AzsTH/34P8A9wOHEFo2fAO+QHhsKnANsAB4FPgPsB+wETCD+sRyV4vkGsbL+dNruaOKjjB8n/nF9ErgfeAE2P7Xx08BaYDzxMuMBwNeArmn5GGLrw6B0TqcD7wHHV1FWo4G51ZTnImBlinkfYstBAI5Iy+tTLvcRHw96RDqnGcC7wOi03RFpu+xjWY9NZdepRLy7AJuAzxTMnwn8Nv1+UorvlMzyrwHvEy997Q/0Ar6TWX4R0JfYCvY54iNjb8ssr3wPO2bm7QusAdplPpvvAH9I7/FpKdb7gJ+mchmZ9tMrbdMKeAp4hFgZ9yZ+lp7IfDbOAVYXnO8OwEagX95/k54+uhOujxukPgbaAUuJM0jioQAACXlJREFUrZUHpzpoHrEVN1vW76bzOyit1474GOP/BY4E9gZOBI5J2xybyuwradkxKe5rqojz0FQ+exTMXwRcmn6/EHgTOCqz/KfAW8C56f06Avh6Zvmo9D7tRXzE9b+Bn2SWn8PW9dxA4JXM65npfEYT699LUqx/TzFVfk7WZT4nNSnbrd4bYATwj4J5o4C/ZeJ9qkQZzsqeW0ubcg/AU4k3JlYSG4DVqQILafp2Zp1/A2cWbHcR8Hzm9d3AzZnXfdM+ty/Y7mngsvT7OWSSmzRvb2IS1LVgu7uA36ffj07bHZtZ/pk0r0t6PQF4tMQ5t0vnemTB/F8B91ZRVqOJSdTqgun2zDqLsq/TvJeBH9SzXPZP8/pk5u2Z4hmdmTcX+F7m9Z+BSVWc0yFpv90L5s8EfpsqtbeBgQXLFwNX1+JzNohYybYqeA+zCfJ3gCkFn81XgW0y854AninY9yI+/EczIJXJXpnlPdJnqn+mfFcXiXEVMDzvv0lPH90J18fZ+XWpj59Ly76a6q2dMutXxrlPpqxfB7bNrPNVYtLcscQx/wGMKph3Qjq2SmxzQirDVgXzFxH7iv8EeA04NLNsR2IDxMhafHZGAvMzr7eq54DfAWMzr2dm3xfiF6nX2bIebkP8QnRSLcp2NFsnyPcAFxcca2FmvzsSvxT1LnJud5BpYGlpU4vtO9JC/IOYCG1P/PDvDYwFkNSJmIhdJ2lcZpvWxA94Kb2IrXKvS1ustl3af6UNxEq60mFpv88XbLctUFFwjGczvy9NP3clJm+HAneWiO3AFMd92vIO7zbESqsqC4jf1rMK7xR+tuD10hQX1L1cPkmsZJ+onBFCeFXSUrZ0PfB14GpJuxAvtQ6r4ny2Tz/fL7LsBGIrz1EhhEcrZyre1LcHsQW7KEl9gcuJLUU7A9sAbYktWoUxVxpK7G+W9XzYsj/ya8RWFQrmVZbvAcDSEMKiyoUhhIWpnA4k9u8sZS0flodZXlwfR3Wpj9ennwcAz4YQ3s0se4RYhx4IzE/z5obUJSM5NG33Ronj9QIO15bduloR36vOxJb9QtsD60MIm4osu5DYSv8fIYSXM/MPJJZxVXXsScQvRvsQk8tt0lSVIcTEOWvz+xZCCJJWAHMy89ZLepMt69ialG021h2Bfmx5n0c/oAOp60kIYbXiDXzDyfyfS1p03ewEubytCSFUfqi/JelB4qWP0XzYf3wk8Y+gploRE5diIwNkb1RYV5AAtSJ+E/0PPqzsKq0teJ1dXlmx1qS/e+U6XyC2xpTaZzEfZMqqlMJ9hMwx61ouNXUb8DNJnyVW9q8TL3eWUvmPoANbV+7PEC+hDZc0K6Sv8tWR1I3YWnA98ENil5PDgNuJSXKxbXYF/pN4OTOrWFlWVb5VqS7+XYjlZZYn18fF91lMTerjQtl64L1abtsKuBL4a5FlpeqON4C2knYIIawpWPY/xKtrpwE/rmkQkvoAE1Ms3yY2Ggwh9uUttU1vYiL9UMGipqhjBxFbtxdk5p0HfAx4L/PlS8C7ki4pKKtdqP7LUrPlBLl5uRL4u6TxIYSlqfVt7xDCrbXYx5PAbsCmULvREZ4i/pF0DiE8WIvtiu2nb4llzxMv93cLIRS2gjS2upbLi8QKqhfwGGy+wWH37EohhFWS7iD2WzsUuKVEy0WlBcR/kAcSyyXrX8Q+1DOB8ZJGhGiFpCXEFoDpRfbZm5gIf7vyn62kz1dzfl8AHgshrKhmveq8AOwuaa/KVmRJPYjlVHh+m0nam9iK9WQ9j2/W0Fwf194LwLmSdsq0dP4/Yh36QjVxnimpY4lW5CeBT9YyKa9skT+QrVtGZwP/BUyXFEIIP8nEv45Yx77M1j4DLMmsX9kwUZWhwD0hhA21iL2YupTtFqNXpKubJwBns3WdO4PYvz77+T6I2M2iRfIoFs1ICGEmsdL6QZr1I+AyxTul95d0kKSzJF1exW4eIN7o8DdJx0nqLukISVeqivFmQwjziP3VbpZ0kqQeknpLulRSYetiVX4BHKp41/anU9znSeqa/qivAa6RdK6kfSQdImmkpBHV7Le14iD8W0y1iKuu5fISsSX4D4p3gx9C7I6whq2/tV9PvBv408BNVCElzw8Any2xfCHxRpRBxMu6lV/1xwAXpc/Efqn8LknLXib+zV+Uzu804qXAqmwevaKeHiBeMpyQPje9iZ+nJ9n6knDWkcDCgsucZrlzfVwnE4h1462KIy4cBVwH3FFNcvsnYAWxnI5M5ztEH45i8WPgdEk/TuX+yVQuPy+1wxBHOXqS0nXsP4k3z10i6Qdp3rvAr4GfKo5YsrekwyWdnzabB+wh6YwU4/nEVuiqNFQdW6uyVRye7fiCY59JbJiZEEKYm52IifB5me33Inbpm9YAsZclJ8jNzy+Jl9a7hRBuILZInkm87P4wsY/cv0ptnC7HDyYmJdcT7/T9C/Fms1J9UCt9hZj8/ZzYcno38Q7qV2oafAjhaaA/se/uLGKr66l8eOmo8pLlpcS7xKcDX6zqnJL9iV0RtphUwzEa61ku5xD7880kDks2gViZF/Yfnlm5Xg1bi8YDX1KJMU/TZbGjicMDXafY1DEOuIDYR3IucWSJnmn9Z4l96y4m/mM/j1jORSkOrdafBqi8U/kOJV7ufDBNy4ETqukichrx/TArR66PayFdnj8WaE8cXeNvxNE5zq1mu/eIozIsJvaNnUtswQ9p+f3EZO+YtN/Hge+xddeQQuMpPoRZ5XEfJybJl1YmycR7OH5GLJsXgMlAl7T+VOKXjl8RGwQGELuzFZWuou1HrKfrpQ5lexQxoc62ng8H7irRjfCvwJGS9kuvTwOmhRBq/HlrblT1/yYzqy3FxzQvBU4LIUzOzN8eWAJ8M4Sw1dOJSuzrUeJd6bc1SrBVH3sYMCaEcGBTHzsd/yDiZb39Qghv5xGDmbVciuMTv0gcfeThHI7/bWBACKHwBvOmOPavgTYhhK/XYdttiVckTwsh/G+DB1cm3AfZrJ4UR4bYiXiH8a7Ebg5vkFoFFB8U0JHYeruW2EJUU18j9lnOw3vAd6tdq/HsDpzl5NjMGkMI4X1JZxFvNsvDEuK4ynl4jtjCXBfdiI0nLTY5Brcgm9Wb4tOrriGO67uGeKny25X9vlJfrX8RLw8OD/EpXGZmZlamnCCbmZmZmWX4Jj0zMzMzswwnyGZmZmZmGU6QzczMzMwynCCbmZmZmWU4QTYzMzMzy/g/Kde8wxgC2foAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Evaluate Training Set \n",
        "y_preds = model(X) #likelihood(model(X_test))\n",
        "mean_y = y_preds.mean\n",
        "y_rmse, y_q2 = eval(Y ,mean_y)\n",
        "\n",
        "print(y_rmse, y_q2)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "HrmpCRWn-Y0X",
        "outputId": "de9465dc-6798-4cf7-89e5-27d968ebe4a1"
      },
      "execution_count": 50,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/gpytorch/models/exact_gp.py:275: GPInputWarning: The input matches the stored training data. Did you forget to call model.train()?\n",
            "  GPInputWarning,\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "tensor(2.1751, grad_fn=<SqrtBackward0>) tensor(0.6254, grad_fn=<RsubBackward1>)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Check manual calculation\n",
        "print(mean_y)\n",
        "test = torch.matmul(torch.matmul(model.covar_module.Ks(X, X), torch.linalg.inv(model.covar_module.K(X) + model.likelihood.noise*torch.eye(420))), Y.reshape((Y.shape[0],1)))\n",
        "test_y = torch.zeros(int(2100/nksip))\n",
        "j,k = 0,14\n",
        "for i in range(int(2100/nskip)):\n",
        "    test_y[i] = torch.sum(test.flatten()[j:k])\n",
        "    j += 14\n",
        "    k += 14\n",
        "print(test_y)\n",
        "y_rmse, y_q2 = eval(Y ,test_y.flatten())\n",
        "print(y_rmse, y_q2)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "ly85oFNGq-1q",
        "outputId": "02eab6ec-385a-444a-b646-6fe5bc40bf99"
      },
      "execution_count": 51,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "tensor([-1.0088, -0.3525,  2.4512, -0.2988,  0.8789,  1.0293,  0.8486, -0.1289,\n",
            "         1.8066, -1.1396, -0.8477,  0.4609, -2.6621, -1.7822, -1.4375, -1.0449,\n",
            "        -2.1738, -0.8906, -2.1475,  1.1299, -1.5518,  0.7676,  1.6260,  0.9355,\n",
            "        -0.2910,  1.1924, -2.5234, -0.9785, -2.2725,  3.7461,  0.0605, -2.5674,\n",
            "        -1.8682, -1.4727, -0.3369,  2.5508,  0.6270, -1.2188, -0.6670, -2.7598,\n",
            "        -1.8555, -1.5068,  1.7529,  0.5898,  1.1699, -1.2422,  0.9922, -1.8311,\n",
            "        -2.8984, -2.0586, -1.0156, -2.5947, -1.0176,  1.5889,  0.6523,  1.5625,\n",
            "        -0.0566,  0.5225, -1.1895, -1.4531, -0.4570,  0.0273, -2.2266,  1.9453,\n",
            "        -1.6006, -0.6016, -1.1357, -0.5029, -1.4785,  0.1514, -0.3018,  0.0312,\n",
            "         0.2148, -0.2314, -1.6084, -0.2461,  1.7998, -1.8926,  0.7568,  0.0293,\n",
            "         0.2246,  3.0596,  0.9229,  2.5273, -0.3311, -0.2002,  1.4561, -1.6250,\n",
            "         2.0068, -0.6602,  0.9609,  0.0596,  1.4805,  0.0967, -1.4727,  1.4570,\n",
            "         0.3730, -0.0322,  0.6631,  1.0254, -0.2803,  2.0283,  2.2021,  1.6104,\n",
            "         2.1631,  0.8623,  0.4883,  2.2305,  1.7656, -3.1211,  0.1426,  1.8125,\n",
            "        -0.1650,  0.4512, -0.5811, -1.2695,  0.5820,  0.6533,  0.6973,  1.4092,\n",
            "         3.4111,  1.2754,  2.0078,  0.9004,  2.0498,  2.8672,  0.0713,  2.1689,\n",
            "        -0.2900,  2.0234,  3.5430,  1.7705,  3.3086,  5.8379,  2.1758,  1.5586,\n",
            "         3.2227,  0.5029, -0.9834,  3.3281,  2.0625,  3.5117,  3.5605,  2.8428,\n",
            "         1.3174,  0.3867,  0.1475,  0.1475,  1.6895,  3.1064,  2.5898,  3.4102,\n",
            "         1.9746,  0.9199,  2.4072,  1.6621,  2.9980,  2.2500,  2.5615,  3.5596,\n",
            "         2.4902,  4.9785,  0.7295,  5.6377,  6.6650,  2.2832,  2.1182,  0.7207,\n",
            "         2.6514,  1.6885,  1.7080,  2.8379, -0.9746,  3.6943,  3.6787,  3.1680,\n",
            "        -0.2695,  0.8418,  3.0186,  3.3164,  2.8242,  4.6484,  1.6113,  2.2354,\n",
            "         1.5693,  4.3750,  2.9590,  3.0508,  1.2041,  4.7168,  5.9082,  2.8037,\n",
            "         1.9414,  2.0332,  2.9043, -0.2197,  2.0088,  0.2900,  1.6963,  2.8301,\n",
            "         1.0312,  0.8799,  5.1084,  1.9805,  5.2188,  1.4453,  2.5166, -1.1943,\n",
            "         0.7432,  1.4316,  1.0176,  3.7803,  1.4062, -0.1982,  3.1211,  0.3877,\n",
            "        -0.3184,  0.3281, -0.7617,  4.5098, -4.7012, -3.3047,  0.5361, -2.1719,\n",
            "         1.4580, -3.2373,  5.3789, -0.6758, -2.6074, -1.0547, -1.0459,  1.5527,\n",
            "        -3.6299, -3.9268, -1.9521, -1.1582, -1.5381,  0.9375, -1.7432, -2.7158,\n",
            "         1.7012,  1.6641,  0.1465,  0.6592, -2.7656, -2.3291, -0.2783,  0.1182,\n",
            "        -0.1738, -1.9854, -1.4268, -0.2461,  1.6914, -0.7715,  1.8379, -1.7188,\n",
            "         2.4082, -0.2295,  1.4990, -0.9482, -0.4990,  2.4766, -0.2480, -1.4434,\n",
            "        -0.4990,  4.0723,  0.9395,  1.9746,  2.8613,  1.6611, -1.9561,  2.1885,\n",
            "        -0.2246, -0.3350,  2.1289,  1.5762, -1.3047,  0.4551, -0.3086,  0.3213,\n",
            "         0.7051,  1.8789,  0.3184,  2.1445,  0.0078,  3.3027, -0.4912, -0.9131,\n",
            "         2.7705,  0.8213, -0.1055, -1.7109, -1.8232, -1.1992,  2.1777,  0.9424,\n",
            "         1.1748,  4.4658,  1.8643,  0.2578, -1.9258, -0.7754, -0.7500, -2.3389,\n",
            "         0.3965, -0.3799, -0.8848,  2.0840, -2.6904,  0.5391, -1.7422,  2.3799,\n",
            "        -1.0391, -2.0527,  1.1719,  0.0264, -0.0713, -1.8604,  1.5938, -1.4941,\n",
            "        -2.5000, -2.7520, -4.4941, -2.2695,  0.4082, -1.8799, -3.4639, -0.6875,\n",
            "        -3.0625,  0.2451,  0.8779,  0.9717, -0.0654,  2.4355,  1.0410, -0.3477,\n",
            "        -0.5361,  0.6377,  2.8125,  0.4355, -1.7588, -0.6592, -1.5244, -1.7383,\n",
            "        -4.6504, -2.3574, -0.7061, -0.1592, -1.2031, -3.2607, -2.4121, -0.1035,\n",
            "        -0.6279, -0.2012, -3.3291, -2.7383, -2.5410, -0.7607, -2.1533,  0.1396,\n",
            "        -2.5938, -3.4512, -0.6738, -0.0527, -3.7197, -2.2627, -0.7363, -3.1934,\n",
            "        -0.2930, -2.8145, -1.0547, -3.4092, -1.4863, -4.1436, -0.2178, -1.2197,\n",
            "        -1.1377, -6.7275, -3.5049, -3.5078, -2.4551, -4.9053, -4.9238, -0.7920,\n",
            "        -2.4980, -4.4707, -4.9551, -4.3750, -5.0967, -4.0615, -6.3721, -6.3887,\n",
            "        -6.5400, -3.8047, -2.8555, -4.8018, -7.1113, -4.5830, -6.9639, -5.7031,\n",
            "        -2.5557, -6.8311, -2.7891, -5.7812, -3.4180, -2.4014, -5.9375, -3.0977,\n",
            "        -3.8525, -4.4697, -5.1514, -5.1504, -3.8164, -6.6895, -2.1865, -5.8682,\n",
            "        -4.2930, -4.7031, -7.5371, -4.2793], grad_fn=<ViewBackward0>)\n"
          ]
        },
        {
          "output_type": "error",
          "ename": "NameError",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-51-7de8c647d60f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmean_y\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mtest\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatmul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatmul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcovar_module\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mKs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcovar_module\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlikelihood\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnoise\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meye\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m420\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mtest_y\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2100\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnksip\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0mj\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mk\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m14\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2100\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnskip\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mNameError\u001b[0m: name 'nksip' is not defined"
          ]
        }
      ]
    }
  ]
}