gmm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c480e801",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import json\n",
    "import os \n",
    "import logging\n",
    "\n",
    "import torch\n",
    "from torch import nn\n",
    "import torch.nn.functional as F\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "id": "01b4252d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "matplotlib.use(\"Agg\")\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set(style=\"white\", font=\"Arial\")\n",
    "colors = sns.color_palette(\"Paired\", n_colors=12).as_hex()\n",
    "\n",
    "# https://github.com/ldeecke/gmm-torch.git\n",
    "def plot(data, y):\n",
    "    n = y.shape[0]\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(1.61803398875*4, 4))\n",
    "    ax.set_facecolor(\"#bbbbbb\")\n",
    "    ax.set_xlabel(\"Dimension 1\")\n",
    "    ax.set_ylabel(\"Dimension 2\")\n",
    "\n",
    "    # plot the locations of all data points ..\n",
    "    for i, point in enumerate(data.data):\n",
    "        if i <= n//2:\n",
    "            # .. separating them by ground truth ..\n",
    "            ax.scatter(*point, color=\"#000000\", alpha=.75, zorder=n+i)\n",
    "        else:\n",
    "            ax.scatter(*point, color=\"#ffffff\", alpha=.75, zorder=n+i)\n",
    "\n",
    "        if y[i] == 0:\n",
    "            # .. as well as their predicted class\n",
    "            ax.scatter(*point, zorder=i, color=\"#dbe9ff\", alpha=.6, edgecolors=colors[1])\n",
    "        else:\n",
    "            ax.scatter(*point, zorder=i, color=\"#ffdbdb\", alpha=.6, edgecolors=colors[5])\n",
    "\n",
    "    handles = [plt.Line2D([0], [0], color=\"w\", lw=4, label=\"Ground Truth 1\"),\n",
    "        plt.Line2D([0], [0], color=\"black\", lw=4, label=\"Ground Truth 2\"),\n",
    "        plt.Line2D([0], [0], color=colors[1], lw=4, label=\"Predicted 1\"),\n",
    "        plt.Line2D([0], [0], color=colors[5], lw=4, label=\"Predicted 2\")]\n",
    "\n",
    "    legend = ax.legend(loc=\"best\", handles=handles)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"example.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "6bb68063",
   "metadata": {},
   "outputs": [],
   "source": [
    "def matmul_per_components(n_components, mat_a, mat_b):\n",
    "    \"\"\"\n",
    "    mat_a: x_mu_T    (n, k, 1, d)\n",
    "    mat_b: precision (1, k, d, d)\n",
    "    \"\"\"\n",
    "    res = torch.zeros(mat_a.shape)\n",
    "    for i in range(n_components):\n",
    "        mat_a_i = mat_a[:, i, :, :].squeeze(-2) #(n, d)\n",
    "        mat_b_i = mat_b[:, i, :, :].squeeze()   #(d, d)\n",
    "        res[:, i, :, :] = mat_a_i.mm(mat_b_i).unsqueeze(1) # (x-mu) * (1/var)\n",
    "    return res\n",
    "\n",
    "def calculate_matmul(mat_a, mat_b):\n",
    "    \"\"\"\n",
    "    input: \n",
    "        mat_a: x_mu_T_precision (n, k, 1, d)\n",
    "        mat_b: x_mu             (n, k, d, 1)\n",
    "        \n",
    "    return:                     (n, k, 1, 1)\n",
    "    \"\"\"\n",
    "    return torch.sum(mat_a.squeeze(-2) * mat_b.squeeze(-1), dim = 2, keepdim = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "id": "42b194bc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "\n",
    "from math import pi\n",
    "from scipy.special import logsumexp\n",
    "\n",
    "\n",
    "class GaussianMixture(torch.nn.Module):\n",
    "    \"\"\"\n",
    "    Input tensors:    (n,d) n= number of samples, d = number of features\n",
    "    model parameters: (mu, sigma) (k,d)\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, n_components, n_features, covariance_type=\"full\", eps=1.e-6, init_params=\"kmeans\", mu_init=None, var_init=None):\n",
    "        \"\"\"\n",
    "        x:              torch.Tensor (n, 1, d)\n",
    "        mu:             torch.Tensor (1, k, d)\n",
    "        var:            torch.Tensor (1, k, d, d)\n",
    "        pi:             torch.Tensor (1, k, 1)\n",
    "        log_likelihood: float\n",
    "        n_components:   int\n",
    "        n_features:     int\n",
    "        \"\"\"\n",
    "        super(GaussianMixture, self).__init__()\n",
    "        self.n_components = n_components\n",
    "        self.n_features = n_features\n",
    "        self.covariance_type = covariance_type\n",
    "        self.log_likelihood = -np.inf\n",
    "        self.eps = eps\n",
    "        \n",
    "        self.mu_init = mu_init\n",
    "    \n",
    "    def fit(self, x, delta = 1e-3, n_iter = 100):\n",
    "        self.initialize_parameters()\n",
    "        \n",
    "        x = self.check_size(x)\n",
    "        #if self.mu_init is None :\n",
    "        #   self.mu = self.get_kmeans_mu(x, n_centers = self.n_components) #gmm에서 mu의 초기화로 kmeans를 많이 사용함.\n",
    "        \n",
    "        i = 0\n",
    "        j = np.inf\n",
    "        \n",
    "        while (i <= n_iter) and (j >= delta):\n",
    "            log_likelihood_old = self.log_likelihood\n",
    "            mu_old = self.mu\n",
    "            var_old = self.var\n",
    "            \n",
    "            #EM step (x의 logprob구하기 -> parameter구하기)\n",
    "            log_likelihood_old_for_check, log_resp = self._e_step(x) #old (mu, var, pi)로 x -> component에 할당하기 (responsibility 구하기)\n",
    "            self._m_step(x, log_resp) #새로 구한 resp로 new (mu, var, pi)구하기\n",
    "            self.log_likelihood = self.__score(x) #new_loglikelihood구하기\n",
    "            \n",
    "            i += 1\n",
    "            j = self.log_likelihood = log_likelihood_old\n",
    "            \n",
    "            if j <= delta:\n",
    "                #업데이트 된 값을 다시 원래대로 복원..\n",
    "                self.mu = mu_old\n",
    "                self.var = var_old\n",
    "                \n",
    "    def predict(self, x, probs = False):\n",
    "        x = self.check_size(x)\n",
    "        \n",
    "        weighted_log_prob = self._estimate_log_prob(x) + torch.log(self.pi)\n",
    "        return torch.squeeze(torch.max(weighted_log_prob, 1)[1].type(torch.LongTensor))\n",
    "        \n",
    "    def initialize_parameters(self):\n",
    "        if self.mu_init == None:\n",
    "            self.mu = torch.nn.Parameter(torch.randn(1, self.n_components, self.n_features), requires_grad = False)\n",
    "        \n",
    "        if self.covariance_type == \"full\":\n",
    "            #\"full\" means the components may independently adopt any position and shape\n",
    "            # (1, k, d, d)\n",
    "            self.var = torch.nn.Parameter(\n",
    "                torch.eye(self.n_features).reshape(1, 1, self.n_features, self.n_features).repeat(1, self.n_components, 1, 1),\n",
    "                requires_grad=False\n",
    "            )\n",
    "        if self.covariance_type == \"diag\":\n",
    "            # (1,k,d)\n",
    "            self.var = torch.nn.Parameter(torch.ones(1, self.n_components, self.n_features), requires_grad=False)\n",
    "        if self.covariance_type == \"tied\":\n",
    "            pass\n",
    "        if self.covariance_type == \"tied diag\":\n",
    "            pass\n",
    "        if self.covariance_type == \"spherical\":\n",
    "            pass\n",
    "        \n",
    "        self.pi = torch.nn.Parameter(torch.Tensor(1, self.n_components, 1), requires_grad=False).fill_(1. / self.n_components)\n",
    "    \n",
    "    def __score(self, x, as_average=False):\n",
    "        \"\"\"\n",
    "        Computes the log-likelihood of the data under the model.\n",
    "        \"\"\"\n",
    "        weighted_log_prob = self._estimate_log_prob(x) + torch.log(self.pi)\n",
    "        per_sample_score = torch.logsumexp(weighted_log_prob, dim=1)\n",
    "\n",
    "        if as_average:\n",
    "            return per_sample_score.mean()\n",
    "        else:\n",
    "            return torch.squeeze(per_sample_score)\n",
    "        \n",
    "    def _e_step(self, x):\n",
    "        \"\"\"\n",
    "        args:\n",
    "            x:              torch.Tensor (n, d) or (n, 1, d)\n",
    "            \n",
    "        returns:\n",
    "            log_prob_norm:  torch.Tensor (1)\n",
    "            log_resp:       torch.Tensor (n, k, 1)\n",
    "        \"\"\"\n",
    "        x = self.check_size(x)\n",
    "\n",
    "        ## loglikelihood = sum_x(ln(sum_z(pi * N(x|mu, var))) = sum_x(sum_z(ln(pi * N(x|mu,var))))\n",
    "        # ln(pi * N(x|mu,var) -> ln(pi) + ln(x|mu,var)\n",
    "        weighted_log_prob = self._estimate_log_prob(x) + torch.log(self.pi) #(n, k, 1)\n",
    "        # sum_z(ln(pi*N(x|m,v))) -> log(sum(exp(log_prob))) -> log(sum(pi*x))\n",
    "        log_prob_norm = torch.logsumexp(weighted_log_prob, dim = 1, keepdim = True) \n",
    "        \n",
    "        # responsibility = posterior = p(theta | x)\n",
    "        log_resp = weighted_log_prob - log_prob_norm\n",
    "        return torch.mean(log_prob_norm), log_resp # (1), (n, k, 1)\n",
    "    \n",
    "    def _m_step(self, x, log_resp):\n",
    "        \"\"\"compute new parameters(mu, var, pi) of gaussian mixture model\n",
    "        \n",
    "        args:\n",
    "            x.size:        (n, 1, d)\n",
    "            log_resp.size: (1, k, 1)\n",
    "            \n",
    "        returns:\n",
    "            pi:  (1, k, 1)\n",
    "            mu:  (1, k, d)\n",
    "            var: (1, k, d)\n",
    "        \"\"\"\n",
    "        x = self.check_size(x)\n",
    "        resp = torch.exp(log_resp)\n",
    "        pi = torch.sum(resp, dim= 0, keepdim = True) + self.eps\n",
    "        \n",
    "        # mu.size = (n, k, 1)\n",
    "        # --resp * x 계산을 위해서 x.size가 (n,d)가 아닌 (n,1,d)인지 확인해야함!\n",
    "        mu = torch.sum(resp * x, dim = 0, keepdim = True) / pi\n",
    "        eps = (torch.eye(self.n_features) * self.eps)\n",
    "        var = torch.sum((x-mu).unsqueeze(-1).matmul((x-mu).unsqueeze(-2)) \n",
    "                        * resp.unsqueeze(-1), dim=0, keepdim=True) \n",
    "                / torch.sum(resp, dim=0, keepdim=True).unsqueeze(-1) \n",
    "                + eps\n",
    "        pi = pi / x.shape[0]\n",
    "        \n",
    "        # update parameters\n",
    "        # .data 사용해야하는 이유 : TypeError: cannot assign 'torch.FloatTensor' as parameter 'pi' (torch.nn.Parameter or None expected)\n",
    "        self.pi.data = pi\n",
    "        self.mu.data = mu\n",
    "        self.var.data = var\n",
    "\n",
    "    def _estimate_log_prob(self, x):\n",
    "        \"\"\"\n",
    "        Returns a tensor with dimensions (n, k, 1), which indicates the log-likelihood that samples belong to the k-th Gaussian.\n",
    "        \n",
    "        args:\n",
    "            x:            torch.Tensor (n, d) or (n, 1, d)\n",
    "        returns:\n",
    "            log_prob: torch.Tensor (n, k, 1)\n",
    "        \"\"\"\n",
    "        x = self.check_size(x)\n",
    "        mu = self.mu\n",
    "        var = self.var\n",
    "        precision = torch.inverse(var)\n",
    "        d = x.shape[-1]\n",
    "        print('precision size:',precision.size())\n",
    "        \n",
    "        log_2pi = d * np.log(2. * pi)\n",
    "        \n",
    "        # log_det를 구하는 이유는? log(var^2)대신에\n",
    "        log_det = self._calculate_log_det(precision)\n",
    "        \n",
    "        ##(sum(x-mu)^2 / 2 * var) 의 과정이 좀 다르다 ..!\n",
    "        x_mu_T = (x - mu).unsqueeze(-2) #(n, k, 1, d)\n",
    "        x_mu = (x - mu).unsqueeze(-1)   #(n, k, d, 1)\n",
    "        \n",
    "        # 각 컴포넌트에 대해서 (x-mu) * (1/var) 계산 \n",
    "        x_mu_T_precision = matmul_per_components(self.n_components, x_mu_T, precision)\n",
    "        \n",
    "        # (x_mu) * (1 / var) * x_mu\n",
    "        # Q : matmul 대신 calculate_matmul(mat_a, mat_b)을 사용하려는 이유가 뭐였을까?\n",
    "        # A : result.size() : torch.matmul -> torch.Size([1, 2, 3, 3]) calculate_matmul -> torch.Size([300, 2, 1])\n",
    "        x_mu_T_precision_x_mu = calculate_matmul(x_mu_T_precision, x_mu) \n",
    "        return -.5 * (log_2pi - log_det + x_mu_T_precision_x_mu)\n",
    "        \n",
    "    def _calculate_log_det(self, var):\n",
    "        \"\"\"Calculate log determinant in log space, to prevent overflow errors.\n",
    "        \n",
    "        args:\n",
    "            var:            torch.Tensor (1, k, d, d)\n",
    "        \"\"\"\n",
    "        log_det = torch.empty(size=(self.n_components,)).to(var.device)\n",
    "        \n",
    "        for k in range(self.n_components):\n",
    "            log_det[k] = 2 * torch.log(torch.diagonal(torch.linalg.cholesky(var[0,k]))).sum()\n",
    "\n",
    "        return log_det.unsqueeze(-1)\n",
    "    \n",
    "    def check_size(self, x):\n",
    "        if len(x.size()) == 2:\n",
    "            # (n, d) --> (n, 1, d)\n",
    "            x = x.unsqueeze(1)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "fe25af69",
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "\n",
    "n, d = 300, 2\n",
    "\n",
    "# generate some data points ..\n",
    "data = torch.Tensor(n, d).normal_()\n",
    "# .. and shift them around to non-standard Gaussians\n",
    "data[:n//2] -= 1\n",
    "data[:n//2] *= sqrt(3)\n",
    "data[n//2:] += 1\n",
    "data[n//2:] *= sqrt(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "id": "73a450d6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "precision size: torch.Size([1, 2, 2, 2])\n",
      "torch.Size([300, 2, 1])\n",
      "torch.Size([300, 1, 2])\n",
      "precision size: torch.Size([1, 2, 2, 2])\n",
      "precision size: torch.Size([1, 2, 2, 2])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcYAAAEUCAYAAACxlJYNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydeXxU9b3332eZJZkkTEIWkhCIBBAUlU26WGpt67UtVbHXulBF+1i0rVSx9nFrlbLUitalt1iLu7hWvTUt0j7eVmsveL1tQaxGkT1AdiATss5yluePyQwzk5nJzGQm6+/9evkiTmbO+Z2T5Pc5310yTdNEIBAIBAIBAPJQL0AgEAgEguGEEEaBQCAQCEIQwigQCAQCQQhCGAUCgUAgCEEIo0AgEAgEIahDvYB04Ha7qampoaioCEVRhno5AoFAIBjG6LrOkSNHmDVrFna7vc/3R4Uw1tTU8K1vfWuolyEQCASCEcTzzz/P/Pnz+7w+KoSxqKgIgJtuuon8/PwhXo1AIBAIhjMul4sHH3wwqB2RjAphDLhP8/PzGT9+/BCvRiAQCAQjgVihN5F8IxAIBAJBCEIYBQKBQCAIYVS4UuNht9vJzs4GQJKkIV6NIJ0E2vx2d3fjdruHeDUCgWC0MKqFMTc3l/z8fAoLC1EURQjjKMM0TXRd5+jRo7hcLjo6OoZ6SQKBYBQwqoXRYrFQXFyMLAuP8WhEkiRUVaW4uJjOzs6hXs6ox2NzUjm5kvzcLFwdPdQerMXmaRvqZQkEaWdUK4YkSUIUxwCyLAtvQIbx2JzMnTWTsnwHOVaFsnwHc2fNxGNzDvXSBIK0I1RDIBD0S+XkSuyqjCpLSIAqS9hVmcrJlUO9NIEg7YxqV+pwRNM0nnzySd544w0URUFVVSoqKrjuuuuYMmXKkKzp0Ucfpbu7mxUrVoS9fvXVV+P1etE0jUOHDgXXd/LJJ7Ny5cqEjr1r1y4OHTrEueeeG3ztzDPP5G9/+1swKSreZ++991527drFWWedxbp165K8MkG6yM/NQpHDrXJFlnDmZtE4RGsSCDKFEMZBZtWqVbjdbp5++mlyc3MxTZN33nmHgwcP9hFGwzCQJGnI3IRPP/00AA0NDSxdupQXXnihz3s0TUNVY/8a7d69m61bt4YJY6IUFBRw0003sXv3bv7+978n/XlB+nB19JCV70ANEUfdMGnr6MnoeUVcUzAUjBlhVBQFu92esSbjuq7jdrvRdT3mew4dOsTbb7/N5s2byc3NBfxx0M997nPB9zz66KPs37+fzs5OmpqaePLJJ9myZQvPPvsskiQxceJEbr/9dgoKCvpYeqH//+ijj3Lw4EE6Ozupr69n4sSJ3HPPPdjtdjo7O1mzZg379u1j/PjxlJSUUFBQkPC1XnDBBZx77rls27aNqVOnMnv2bLZu3Rq06DZt2sTWrVu5/fbb2bBhA11dXSxZsoS5c+fyox/9CIDf/va3/PWvf+X48ePceOONfPGLX+xznqKiIoqKijhw4EDCaxNkhtqDtRTkzsSuyiiyhG6YuDXDL1QZOmcgrhk4Z1a+g4LcmbxXs1OIoyCjjJkYYyZFEU4Ibzx27dpFRUUFeXl5cd9XU1PD2rVrefXVV2lpaWH9+vWsX7+eF198kSlTpnDfffcltKadO3eydu1aXnnlFTRN409/+hMAjz32GA6Hg1dffZV169bx3nvvJXaRIXR1dfHMM89w5513xnyP0+nkuuuuY8GCBbzwwgtBUQRwOBxs3LiR1atX84tf/CLp8wsGF5unjfdqdtLg6qLDq9Pg6sq4QIm4pmCoGDPCOBzZv38/S5Ys4d///d/DxOGss87C6fRn+23fvp2zzjqLwsJCAL7xjW/wz3/+M6Hjf/rTnyY3NxdJkjj11FOpq6sLHvPCCy8E/OJ1zjnnJL32RYsWJf2ZUP7t3/4NgFmzZnHkyBE8Hs+AjifIPDZPG42732fn9ndp3P1+xq22eHFNgSCTjBlh7M/NOVACrtR4nHzyyRw+fDhYiD5lyhReeOEFLr300rA6vKysxP7wFUUJdn8B+oiL1WoNe286rz90jYqiYBhG8P+9Xm+/nw+sLWDFZ/JnIxiZuDp60A0z7LXBiGumA4/NSen02Zwy7zOUTp8tylpGGGNGGHVdp6uri/b29oz819XV1e/mPmnSJM4++2x+9rOfhQlhT0/sP/R58+bxzjvvcPToUQCqq6tZsGABABUVFezcuRPDMOjq6mLr1q0J3Yv58+ezadMmANra2nj77bcT+lwsKioq2Lt3L16vF5/Px5tvvhn8nsPhEMX3gpSoPViLWzPQDBMT0ELimsMZUfM58hkzyTfDhZUrV/L4449z1VVXoSgKeXl5FBUVcdVVV0V9/9SpU1m+fDnLly9HkiTKy8u5/fbbATjnnHP485//zDe/+U0mTJjAzJkzE1rDd77zHVavXs3FF1/M+PHjmTNnzoCu6bTTTmPBggVceumlFBUVMW3atKCQL1iwgOeee65P8k0iNDQ0sGzZMtxuNx6Ph0WLFnHttdcG3cCC0U0grlk5uRJnbhZtIyQrNTQ2CuGx0cbd7w/x6gSJIJmhvrgRSl1dHV/60pdYvXp12DzG8ePHU1VVNYQrEwwW+/bt49ixY0O9DEGaGYnlGqfM+ww5VoXQ6KgJdHh1dm5/d6iWJQjh2LFj3HXXXbz55ptMnDixz/fHjCtVIBCMLEaqS3Ikx0YFfoQwCgSCIaG/BJWRWq4xUmOjghOIGKNAIBh0+hTvF+Rw0oTZdLk1Wto6qT1YO2Lb0I3U2KjgBEIYBQLBoBNqDWomOGwWJEkiL8uCKvk73HR6fGSFJLHAyHFJBmo+AwKeqe5AgswgXKkCgWDQCbUGrap/iLgs+S3CgMtUBuGSFAwJQhgFAsGgE5qgoioykuTP3Ay8psgSdptl0NvQCQQgXKkCgWAICG1KLusGsqQgSeDT/E0yAi5T4ZIUDAXDUhjXr1/Pr371KzZt2sT06dOHejlpRcxjTHwe42uvvcbLL7+MaZpIksSVV17J1772tSSvTjAcCU1QKRrnQJYlvyia5qBM7hAI4jHshPGjjz7i/fffp7y8PK3H/e///m9uvvlmdu/endbjBpg+fTr3338/n//85+O+T8xjTJxJkybxm9/8hnHjxtHc3My3vvUtZs+eTVlZWdLHEgw/Qq3BQCG/yOIUDAeGlTB6vV5Wr17N/fffz9KlS9N67BUrVrBv3760HjOU3bt3s2LFirgjnMQ8xuTmMc6bNy/4dUlJCYWFhbS0tAhhHIUIl6lgODGshPGXv/wlF1xwQdQWPaOBZOYxPvfcczidTvbu3cv69et59tlnKSws5JFHHuG+++7j5z//eb/n27lzJ8888ww5OTn84Ac/4E9/+hMXXXRR2DzGtrY2rrjiCr785S8ndS2BeYxAsCF5JIF5jKGiGSAwj/Ff//oXt99+e1RhDGX79u10dHQwY8aMpNYpEAgEyTJsslJ37NhBTU0NS5YsycjxH3roIU4++eSMHBv8cbeHHnooqc+IeYyJzWPcv38/K1eu5Gc/+1m/w6AFIw8xokkw3Bg2FuM///lP9u3bx5e+9CUAmpqauOaaa/j5z38e5mpMlc9//vP84x//GPBxBkLoPMbc3NzgPMaXX36Zjz/+OPi+TM1jTOcw4MGax3jo0CFWrFjB7bffzuzZsweyZMEwpE8HnHx/cb8oyxAMJcPGYrz22mvZunUrb731Fm+99RYTJkzgiSeeSIsoDhfEPMbkqKur4wc/+AE333wzZ5111oDWONoZqVbXSO2HKhjdDBuLcawg5jEmPo9x/fr1HD9+nA0bNrBhwwYAfvCDH/CZz3xmQOsdbYxkq2uk9kMVjG7EPEbBqGAsz2MsnT6bsnxHWE9RzTBpcHUN+8G4sdbe2uXB7dVG1BxGwchBzGMUCEY58ayu4U60EU06EnnZthE3h1EwehDCKBCMcEbyYNxAB5zQfqiuLjeqhIg7CoYMEWMUCEY4oX1HFVkaVi3VAh1t4rlEI4v7T5n3GRF3FAwpQhgFghHOcB2Mm2pSkKujh6yIuGO6LOBEhDpdDOa5BOlFCKNAMAoYji3VQksxINwlGi8pKFMW8GBm747kTGGBiDEKBIIMkWpSULS4YzoEJZGayXTVg4r6zJGNEMZB5oILLuDiiy9myZIlXHrppfzXf/3XgI+5adMmbr31VgD+9re/8ctf/jLu+zs6Oti4cWPK57vgggvYu3dvn9dbWlr47ne/yxe+8IW0N4EXjDwGkhQUsIB3bn+Xxt3vp8XK6k+oA1ZeItmw/QnoSM4UFghhHBLuueceXnjhBVavXs3q1atpawv/o9c0LeVjn3322dx4441x3zNQYYxFdnY21113HWvXrk37sQUjj2ilGAGX6FDQn1AnauUlIqAjOVNYMIZijP842MbP3tjLgWOZ+cU8aXwWPz5vKgsmJ+56Ofnkk8nOzqa+vp6HHnoIVVU5ePAgXV1dvPDCC7z++uu8+uqr6LpOTk4Ot956K5WVlfh8Pu677z62bduG0+kMG+YcGPkUmGbxhz/8gZdeegkAVVV58MEHuffee+ns7GTJkiXY7XaefPJJjh49yn333UdTUxMej4fzzjuPb3/724C/wXvgeHPnziVWT4icnBzmzJnD9u3bU7qHgtHFcEsK6i92mWgXnkRip8M5U1jQP2NGGNf8aS8HXZl7WjtwrIc1f9rLpu/OT/gz27Ztw+v1MmnSJMA/03HDhg1kZWWxY8cO/vKXv/Doo49itVp55513WLNmDU888QS/+93vaGho4OWXX0bTNK699lpKS0v7HH/79u089dRTPPbYYxQWFtLd3Y2iKNxyyy19Bg+vXLmSa665hrlz5+Lz+fj+97/PKaecwpw5c/jxj3/MmjVrmDdvHn/+85955ZVXBn7DBGOCdCUFpSPDsz+hTjQbNhEBHW4PBYLkGDPCOJy47bbbsFqt5OTksG7duuDQ4i9+8YvBqRVbtmxhz549XH311cHPtbe3A35BXbRoEaqqoqoqX/nKV/jXv/7V5zxbt27la1/7WnBkVXZ2dtT19PT0sH37dlwuV/C17u5uDhw4QEFBAXa7PTg0+Nxzz+Xuu+8e+E0QCCKIJX7pzPCMJ9SJWnmJCuhwzBQWJMaYEcY7vzqVu9/Yy/4MuVKnjM/ijvOmJvTee+65h6lT+743VLhM0+T888/nu9/9btrWGAvDMJAkiY0bN6Kq4b8Se/bs6fN+SZL6vCYQRCNRSy+e+KVa9pEsiVp5wk06+hkzwrhgspPqaxN3cw41Cxcu5Kc//SkXXXQRJSUl6LrO7t27mTlzJmeeeSZ//OMfOffcc9E0jTfeeIMJEyb0OcbnPvc51q5dyze+8Q3Gjx8fdKU6HA7cbjeapqGqKg6Hg9mzZ/P000/zne98B/DPw1RVlcmTJ+PxeNixYwdz5szhzTffpKOjY7Bvh2AEkoylF0/8BnMCRyJWXqSAuj0+DGDOrJmikH+UMGaEcaQxd+5cvve973HzzTdjGAY+n48vfelLzJw5k4suuog9e/ZwySWX4HQ6OeWUU2htbe1zjHnz5nH11Vdz/fXXI8syFouFBx54gPHjx/OVr3yFyy+/nNzcXJ588knWrFnDgw8+yGWXXQb4rde77rqLwsJC1q5dy7p165AkiTlz5kQVYfAPG77gggvwer10dnayaNEiLrzwQq699tqM3ivB8CQZSy+e+GWyE06qBAS0VhTyj0rE2CnBqGAsj50arpwy7zPkWBVC5c4EOrw6O7e/G/beeKOzag/WholPwHU5HMRnJI/8GsuIsVMCgWBISKaWL17NY2QnnNYuD50eH3NmzRxQd5p0IAr5RyfClSoQjGA2b9lGdXU1rpZG8otLWbx4MYsWDo9YerQkFc0Eu1XllHmfCYvH9Zf4kojrEhj0pt3D0c0rGDhCGAWCEcrmLdvYt/NDNjxwD5WTJ1F78BDrNzzGtvJJnP+VfxvyqQ6RYtft8ZHvsFPgsEWNxyWS+BIrbjl92lRybJZBj/WJDNXRiRBGgWCE0lBXx71rf4ru02g/3kb5hBLW338vSBK54xxDngwSWapht6oomAMqu4jluiwvyMHj8WW8pCMSUcg/eAzmGC8hjALBCOXqSy9C8/lwu90A9Lh7qJg0EcMwg+IxWAIRSbRSDUeWDY9XA9MIvi/ZsotYrktL77FCGazhxqKQP/MM9hgvkXwjEIxQplRNoaurK+w1VVGxWi1hrw1FMki0htyabmC1KGHvSzYeFytJp6G1UzTtHsUM9hgvYTEKBCOUth4dm82Gx+NBURR0Xcen+VAs1rD3pVsgEkn4ieby9Pg0cuxWtF6LNpV4XCzXJRC1pGMwYn2JuPgG0w04GhnMJg8gLMZBZzTPY/zb3/7GlVdeyaWXXsoll1zCc889l/I5BP1z1NVGfkk5jpxcvF4vjpxcejSTHp+ZsVFPm7dsY+MTG1C9XVRNLEP1drHxiQ1s3rIt7H3RSjUwTWqPtg94AHG0WY0BwWzt8mCxqNhsFjo9voFebr94bE4+NXsWU0vzKc7LYmppPp+aPSushCSZOY+C6Az2GK8xYzFqf/877tVrMfbvz8jx5SlTsN/1E9RPfarf9wZ6pe7atYtrrrmGBQsW4HSe+CMJtGpLhbPPPpuzzz477nsCwpjuYcLjx4/ngQceoKioiM7OTq688kpOPfVU5syZk9bzCPzYPG3U7DlA5eRKJpdW0NbRw/aP/L1tKydX4rDA7l2f8Kv1v6Zmz/60lHJUV1dTkJuDM8/f+D7wb3V1ddixY2Vr7t6z1y9sgWtI8LyJWlzZNgtujw9Flih02MiZldnEoxnTpjIu24oESPitmHHZVmZMm8qBGv/DwmD1eh3NDHb275gRRvfKVRgHD2bs+Mb+/bhXriLn//0x4c+MtnmMs2bNCn6dk5NDZWUljY2NQhjJXL1hrMSPx594nI1PbKAgN4e8nJygZQeEnTdZF5+rpZGqiWVhr+Xl5LCvLtyhlc5szUQTL4ZCgMoKcoKiSMi/ZQU5HOj9erDdgKORwc7+HTPCOBwZzfMYa2trqamp4Y477kjfDRuhBNyPBbk5VE0so72zM6pIpZNELLtUMv3yi0tp7+wMHg+gvbOT8772dUqnz+4jsAPN1vTYnJw17zSyrSo+zcCr6TEFbygEyARME0IHzpim//UAoglAehjM7N8xI4z2VStxr1mLsS9DrtSqKdjv/ElC7x3t8xiPHj3KzTffzK233kpRUVFC92Q0k6j7MZ0kYtmlYmEtXrw4KOp5OTm0d3Yy87TTufdnq8hzZKc1lT4g3A6rBUkCq6pgUWS6el2lkYI3FAJU39rJ5MJcDFNCkgKiaFLf2hl8j2gCMPIYM8KofupT5Ly+aaiXAYzueYytra18//vfZ+nSpXz5y19O+1pHIom6H9NJLMsuv/jEZJRULKyAkFdXV7OvrpH84gn8+C6/KCYqsIm6bwPCbZgmiiQhS2AgYVUVvD6tj+ANhQDt3rMXZ/apOGwWFFnCME0kCZzZdrzTZ4f1ehVNAEYOIit1mLJw4UL++Mc/0tzcDPhHOu3c6e8HGZjHqGkabrebN954I+oxPve5z/HHP/4xOHWiu7sbj8cTNo8RCJvHGKCpqYmjR4+GzWME4s5jbGtrY/ny5VxyySVceOGFabkPo4GASIUSKVLpZvHixbR2dNLW3oFhmLS1d9Da0cnixYuD74nM9JNkGZvVQr7DFrc596KF83ns/rW8+uzjPHb/WirKSxNupJ1MhmZAuH2aDvjdk5IEFlWOmmkb2Ww81azX0LWWTp/NKfM+E/d+tHd7/NYiJpIk4fXp2BTCri1aJu1gkui1CPwMG4vR5XJxyy23cOjQIaxWK5MnT2b16tUUFBQM9dKGhJE4j/GZZ57h0KFD/O53v+N3v/sdAJdddhkXXHBB5m7UCCCa+7G1o5Ol13yr38+mWv8WzbJbes23YmaOqoqCzapimiYej0ZZEu7QZFyYybhvA8eVMHB7NSyqgixJdHu1mOtKVxwqkfhr6Ht8Pg1br8tXN0xUafhknw5215jRwLCZx9jW1sauXbv4VG+5w7p16zh+/Hi/8SwQ8xgFw38e44ms1CbyiycklJUauaFlYg5hQHinTHBimv4ifLVX3xKdK5jMOuPNaHy/190YeAiobznKaVUVabn+ZB8wEpmzGPmeLLsV0wSvpmPoJ6zcaPMnBxMxM7Iv/c1jHDYWo9PpDIoiwOzZs3nxxReHcEUCQfpYtHB+XCGMtnFnuvwg9JwAum4ERRESz+hMJoYWy7rs9viiWDVZfLjvMOXFhQOKzaWUfZtA/DXyPXpvRx9VkfH2CmM0y3mwu+CIcpHkGTbCGIphGLz44ot88YtfHOqlCAQZJ9bGrSoKSkQWQLo2tMhzmibYbSoer4Zp+Jt8J5PRmagLM1aCjAxRHwLKiwsH7BpN5QEjEfdw5Ht8mo5sVdF7uw5FS/4ZCremKBdJnmGZfLNmzRqys7O54oorBnQc0zRjFqMLRg8j/eccq0Gy3apgs1rIsluxWi1Ispy2DS3ynB6fhmmaWFQlI63kAsRKkMnqzeoMJZXm59GSTOJZTLGIbFZuShJ2m4UipyN43Mj3+DSd491e6l2dMZN/BrsZdrRrydTPdjQx7CzGdevWcfDgQX7zm98gywPTbdM08Xq92GyiWmg04/V6R7QwRtu4kSTsFtlfF2f6N3LZqqat/CDynKoEXW4fNqtKh1fPaElBNOsyHVZNLGus0+MjK8Ri7O/YAVenoshIioxh+mso3V4NuyKFJSaFupDdHh8G4LBZYrpIM+XWjOeeFeUiyTOshPGBBx6gpqYm2O1loHR0dNDQ0EBZWRlWqzVu/Z1g5BF48GloaIhZQjISiCYKNouKYUC3x4dVVVAVv7Xo6nKnZUOLdk5MkwNNbUH34mA+TqajBjGWy7TbA27NSOjYfRKJdAO7zYLbeyIpKdIV27j7fWoTdJFmwq2ZiHtWzIxMjmEjjHv27GHDhg1UVlYGSwYmTpzIww8/nPIxNU3D5XKhaRqSJAlhHGUEXKgdHR3BmsyRSDRRsCkynt7N2NB1vLqOib9JdqbOOZTdWPqzahJJWIlljdltFj7cd5gFp5yE3Wqh2+vjw10H4zYVCBVXVZb9pSyGiarIaLoBph5m5UX7nMOmcta80+j2aME1Z+K+iybl6WfYCOO0adPYtWtX2o8bEEeBYLgSTRTsXR4KHLaMJUwMR/daLKsm0YSVWNaY2+PjtKoKTN3A7fFiGianVVXwXk1nQuJqmCZ2i4Ku+7vayJKCRZU52umO+TlJlsmyqhimiWyaYWt+r2Yn06dNpbwgBwvQ2hXe/CFZRNZp+hmWyTcCwVgjsjPKJ3v2ZjxhYii7sSTTiSXRhJVYCTPlBbkU5NixqEq/CS/R5v6Z+LuCS5J/eoYk+dsiynE+Z1EVTBN8mhH1nDk2Cx6PD59Po8BhG9B8xsGeVTgWGDYWo2B4kKnxSILoxLrfw9GiSxehFmBHeztmp4uTivNYde/TlE2c2Of3LVGLKPSeFTtzcNj9CTPZvRmvilXF3VuOEsuiiurWlmR6vBqqcuI1r0/HHuLWDv0ckoRD8Q+jkmUJr+ZPbgqcM0zoZRmbVSFHkjhr3mm8s/3DpH/Gw80tPhoQwigIEhiPZJVA6zpO3c4m1u/6gA8/uZzbll0+1MsbdfQ3jmowEyYGs+h8+rSpjMu2IUsmRpaVjnYwNI1rllzM5Vf5Z4CGimMyCSuBe0ZItxdNN5AlpbfPqoLXa8T9fCy3tqGHd44J/Xzgc9OnTaWyMM8/dso0UWUZh81Cl8cHpv8zAaGXFYUsqxqcyuGwWZgbZ7ByrJ/RaH6IGiqEK1UQpLq6GqsErqPNmIZOlj0LWTJ5u/pFNm/ZNtTLyzibt2xj2c0/4eIrr2HZzT/J+DWHjqOSZQlnXi4FuTlUV1dn9LyhbN6yjcderOak4jzMzqNo3Z1xG3sPFI/NSWVhHhZFxufxYLVaOGnyZFRVpXLypKjXn0odXqiV6dX03kQtv9XW3+dTdWvbPG14vRqdPR56PFqw6bkkSdgsJ0ptXB09IElkWdWwCceyBA6bGtXFG7CyKwpzGZ9jZ3pZAV/81GwM58SoaxaiODD6Fcaenh4+/vhjOjv7Boi3b9+ekUUJhgZXSyMdx10oioKiKCCB1WIFQx/UzXooCFhvqreLqollwYn3mRRHV0sjeTk5Ya/l5eQwa9qUqPG3dAt34JqvWXIxhqbR1dlB3eFaOtrbM1Z0Xjm5MjiayTAMTNPEMAzKSks5XFdPXk4OrpamsM+kMjUjNO6mStDl8eHTDXy6Effz0WKfiZ7fY3NSVZpPnsNGls2feGMEYpMSwc/UHqzFovr/viR6hxr31qvaLGrUxgOVkytx2NSgW1iWQJVlPnf6VDEpIwPEdaV+8MEHXHvttfh8PjRNY/ny5Sxbtiz4/WXLlvHee+9lfJGCwSG/uJS6nU1k2U/8YWq6jmq399msIgmNlanWLCRZwufuHjFxyqEYJhxtZuL0U2by0zt/QnG+IywD85lXfh/X7ZoKgWuunDyJ9uNt/ochoKW5EadzXEayGvNzs/D6dLJs/iHbuq4hyWCz2Xj+pZdijuNK1q0cGXfDNDne7QmKU7TP95f9Gu/8HpuTT82ehVWRg1ag2jufsdvj4/DRjrCawg63D1X2n8fAb83KgCxJUV28+blZ2Cz+7TpoZPZO8BBlGeknrsV4zz338MMf/pDt27fz0ksv8Yc//IG77ror+P2R3G1E0JfFixeDrODxeTBNE5+moek6uXnOuLMDQ62tvCw7bY0HaKs/wPjc7EGxvNJBLOutvweCgRBtZuJVV3wLR05unwzMudMnpd3tGrjmw3X1wYchRVHwetxJZTUmk2Hq6uhB03XcXg1JVnrF0aBm5yds+dt/95kZmSqpWJkDaddWObmSXLvFn8EaQiCuGel2PdLWRafbi6YbmIZf7AwTNCO6i9fV0YMsSeFTSXqzXpNtmyfon7jCuGfPHr75zW8CMHPmTF588UUOHDjALbfcgtHbaFgweli0cD5nX3g5uiHR0+PGlBWc44sxJDnuZjX/HWMAACAASURBVBVqbR050ozFYsFqtdDS3DQkcbNUGIphwosWzmfpNdehWR3sq2tEszqYMesM8nJzw96nyBIVpcVpF+7ANT//0kuoVitZ9ix0XSdvnDPh0pBkBg/DiXihT9PRDJPW9i4aW47w07X3MPezZ/PH//dn/u+KH6RlmG6ycbdUeqr2+awJhmH2TlX2f6+jx9fn3LUHa3H7dLo8Pv+YKtPEZxhs+WBv1HXWHqxFMwwChzbM3s5Pmp5yWYYYXhybuK5Uu91Oa2trcMZhTk4Ojz/+ON///ve54YYbBmWBgsHltmWXc9qMacHZgTmF/c8OdLU0UjWxDADN68FusyIh4fF6AP8Gvq9ueJcaD2SY8ECIHEfV5YNxhtknA/NwY0sft+tAhTtwzVv+9t9MmzqNa666gql5U+j0aGzffTihBI5ku65EZlAe90rUNrby3Rt/mNBg4MrJlRQ5HVhVGZ9m0tLWmbYMzIG0a3N19DA+x47Sa236Y6fg0w2OHO/q9z70l0lq87Sx5YO9fO70qaiyhE8z8Go6bp+eUlmGGF4cn7jCuGDBAjZt2sTVV18dfM1ms/HII49www034Ha7Y39YMGLpb3ZgJKGxMtVqQ9M1ZEnCavX/uWba8koHiUy8zwSRdYzXXX8DV33zwr7DeXcforXDb9GmS7gD19ZQV8fi879Gy5GjuL06uXl5MTvDRJYMFDtzkh6NFS1eVzp9dlyBDW7kFoVsm8XvVrSBTc1N24beJy4pSdhtKkVOB0yfHVe4ag/WUpKfw7hsf49n0/Q3Bujy+PyWd4xSi2TipnJbHX/9e2dayjJEG7n4xBXGH//4x3R19X3asVqtrF+/nh07dmRsYYKRQ6i1VVRUQkNdLRISkyon0tbeMSiWVzro74Eg3c0PotUx3rf6xwCc/5V/C9v85k8pZuk116VduBctnE/p9O8wIcJSirZJRrMybDYVTdP9StBLKl1X4rkxa21Ozpp3Gg6rhUC740A2p83ij1GmY0OP1SAgcqJGrP6t7+2qpWriBMoKcjCB+tZOdu/ZC5A26yxdta2ijVx84gpjQUEBBQUF0T+oqpx55pkZWZRgZBFqbbW3deAsPQlJljjW0ZOWDXw4dOPprxg/FWJlwm54+D+YP6W4z+aXrCWfKIluktGsDI9Xw25VcXt8cbuu9NdAIF6f07mzZpLdWwgvS5K/1MMEyTT9HW0S3NATaWIQrUFA4FqjWrBhYpfFezU7ORByzESs4aFADC+Oj+h8I0gLmdq00y1IqYpsJso5QmOzAYYiHpvoJhlVQCV/csmR410x3XuJxLNitTUz8FuvPs3AqioYkonMCXHUDTOhDT3ZmFp/Dwv9uSJDRdjpsOHxajGPNRSINnLxEcIoGJYEBOzwJ//CoqrkO7KCZQqQmiBFE9nH1z/Es0870bw9cYUyEyIWrY5xKOKxiW6SsQT0yPGuuO69eCLy+BOPBx9UbvrhzSy/7hoctqzgaKhZVRUosoRX07EoMhj4W8T0LsHj0xLa0ANt6JTeFnG6qce12vp7WOjP9RsqwookkZNlQdNNZMl//oFkk6YD0UYuPqIlnGDYEVoXCSaSaVJ3qBZXm/+PNtUyhcgWbKah4+5oo+NIfb/dbjJRzhGtjjFddXzJkEjNn8fmxGpVycmyYbVaQJITnvgRS0QcFoI/52+c/3UuvfCrHD6wj6bmI5i6wWlVFXR5fMEONrpp+n9umOi6iavLy+GjHbxXsxMgZulBaBs6RZawqgoOmwWk2KUYtQdr0UxQVBWr1YKiqmgmwWsNdNaRZBmr1UKW3YrNasHt8fWph3T7dBRJxqYqwfPn2K3UtxxN8ieVXkQbudgkLYyGYYT9JxCkm1ABs1jtmJK/8Ly5yW+TpCpIkUX8zU2NWFQVQ9f6LZrPhIhFq2Nces11/VrCmag/i7dJBtyQhQ4bnW6vf+KEVaG1yxM3gSSwTqfDhs1qQZJPbDe6YbJ71yfBn/O3llyG7tPweDy0NDcGrUoZ0Exw2C2osoxp+l2ox3u8/ONfNUFrL14tZWgbOoleg7O3f2k8q02SpGAMU5HDB53XHqxFR8JmVVFkCbO39Vteto0ipyPsQSDQASeAYZr0eH2UFxcm9TMSDB4JuVI/+ugjVq9eza5du/B4/LVppmkiSRI7d+7M6AIFY4PQ2F9nawuTe7uNFJdMoKHuIKqs4NPcA8pyjXRder0eZElGtZ5wwsVyj2aqnCNWbDZWoshQ1J9FukJ1TcPjNXF7tbiiGFinx6vhsFuwWVXcHn83Jbdm8Kv1vw4+qFRMLA+2pfN6/GVgiuyfp9je7SEvy4os+wXVp+komEE3aH/xvtA2dIG6e0kCVZFiumArJ1eiYOL1+oKvhZ7T5mnD1eUmu1cYA+5RVQJJkdF148R6eutZfLoRdrx0xhhTmY4ymBNVRhoJCeNtt93GOeecw913343dbs/0mgRjjMjYX83xVg7VHkCSJIrH+7Oimxrq/S2zrI6UBSmyiF+WFXw+L+VlE4PviWeNDiTBKJmkn3jiNxT1Z6mk9keus8vtw2ZRUVWZQ0fbqT1YS82e/ajeLhae/XnyneOYNLGM7p4eWo74XYyBmF5+bhYery+sHVro+ftbn6ujhyxVxu31t2cLxFHrjnXEFIJErtlhs6BpGnrv/6uSfziyRzPRdCMsZivLEj5NDx4rnRmgqTwsiQL/+CQkjPX19dx0001hrgSBIF1EZnxOKCun/tAB6usOMy5vHFaLhZyCwoTcjPGItPpyiifSfbyVz571Wa759repKC/lQO1B3j/QnJbrCpBsZm0s8SvMd6L6OthZsxerzU5xSSn5znF9Nux0l7ekktofKSyqBJqm0eHV/RYX/geVfTs/5JYfrqCjsxNHdhZ2m43Kykp6vPqJ+OXkyrjnD12fZoJVVbCoMoYk4bE5g8lF4A8FBZKLPtmzN2oz8MrJlTgdNhTJL2Zmb8go8ppj3ZeWti5qD9YGE1u6O9zkO+womBnJAE3lYUkU+McnIWE899xz2bp1KwsXLsz0egRjkMiMz+LxBZiGSUP94bR3oYm0+rbtb+GiL36anq5OOjo6mXXa6ZwxPzutT87JlnpEs1Y62tvJs5gcrD1ISdF42tpc7G9vw57tYOKkSo57/e/PRL1lKqn9iYjpooXzyf3mxWiaRuuxRnw+HxUTK1BkBUk2wsY0xTt/fctRpkxwYlEkJEkODjeQITj4N5EMzHD3r05etgWLYsGj+ZueR15z3HVFFOLvDghuBjJAU7HoRYF/fBISRo/Hw/Lly5k3bx6FheEB43vvvTcjCxOMHaKVLdhsVipmnM5j96/t9/MDsZDO/8q/URyxgWuGGfPJefOWbTy38RlcDbUAOCdM5sqrr064l2yAeKUe0USl/biLhsYmHn/yKVbecSs5DgednV1YFJmO4y42/W0b86cUZ6TeMtHUfo/NyfRpUykvyEGR/NmjPR4NRSKmmFaUl5JjVYJuUgN/okq21T+wN3CeWOf32JzMrarA7dVQ7f44pIREj1dDMs0wK6i/jjEBKwpJwmFTMUx/naZN9Y/j+ucn4decTMlDujrWRCMVi14U+McnIWGcOnUqU6dOzfRaBGOUgTTwHqiFlMyT8+Yt23h8/UO4O9qwqv5EjrbGAzz6qweBm2KeL9l6xWiWiKzIvPTbV9j0h9/T0XqUH9zwAyorJ3Ng/wGee/lh/vf9j5l//9qMNQ0IbOyPx3gI8dicLDjjVMZl+xvI+6MuEtk2C8e7vRw53hVVNEI3aEmWsVtVTBO8mt6nDVs0YQm6BCV/tqfUmyivKjJeXU/KCgr8LiiKgiRJyPiTDE3A7fFnkTa21UW9L5kQvERJxaIXBf7xSUgYly9fnul1CMYwA8n4HKiFlMyTc3V1NZq7y1/X1jvUFwk0d3fc8yUr/NEskdX3/gfv/c9/o3k9bH1nK+/8zzvouo6sqEyfPiNY15nJpgHxHkK+/73vku+wIcu9o5dCyhMCDQCibbihG3S2VfE33+4dp5RI3Cv0wUbTDWRJQZYIvpaMFRT4XbBZ5WBPVrP3GMPZzZhKsb4o8I9Pwp1v/v73v1NdXU1LSwvFxcVceOGFfPrTn87k2gRjiFQzPgdqISXz5OxqaUTXNSwh5R2qotDj8fTbcEBSbBzY728onV9W2W8iUaQlUjZxIn/p6ERW/E27ZVlC13VKyyvChG8g1neoK1QCGlo7/QkqvZtlrIeQhro6KgvzgoN0zd6epkYCghK6QU8vtwbnGFpVBa+mgyRx0gQn+bmf6bfHaqA7TsAdqyVpBQV+F2yGv0NNwL/r0/SgwA7XEodULNfhYO0OVxIq8H/llVdYsWIFRUVFnHvuuRQXF3PzzTfz8ssvZ3p9AkFcBtqRJplJ7/nFpSiKim740+69Xi9dXZ3oPi89ncejdswJWFn52RbOOH02J02ZiuFLflxboBlAblE5Xq8Xw4TyislIshLWaGAgTQMWnHEqkwtzscgyFkVmUlEen5o9K1goH9kgAfzie/WlF2GYJmbv/PqA7S1LUkIWWyDBxjD8lp6q+F2qOVlWcrIs/pZqMYYgBwYfa4aJIkGXx4fPMGjv8cX9WcZax3s1O6l3+X+fdMMfI/Vp/gzZ+pajSQ1lFoxcErIYH3/8cZ566ilmzJgRfO2rX/0qN9xwA5dccknGFicQ9Ec6Bgwn+uS8ePFiHl9fi7ujDfCiaT4wTRTVQoHTGTW2Gc/VG/g30aShgFUdSDY61tIU1e2civVdObkSh82ChERoyDXXbjlRKB/DTTulagpen44iq0GLMfBvh9uXkMU2Y9rU3u4yvS9IoEr+jjIeTUMieklBfy7BZK0gm6eNgzXbomaRRi1xsCicNe80uj3asLIgBQMjIWFsa2ujqqoq7LUpU6Zw/PjxjCxKIIhGvOzTwRgw7D/mCp7b+AxHD+4GwJ6dQ/nECvKdTtraO/rEGmO5ej/euy/hpKFo151Itm4yBGJ1oaXKEuGJSKEPIR6Pl+amBnxeDwdqD+Jw5CCNG4dFVVAVCRPo9ur8/f2ahISirMBviRq9rdWC+ij5W7H5TH8tZDTXbCZcgtGOGZmopZkEhybLpklOYS5TJsymw+3jSFv0ZKORxnB1HWeahIRx7ty53HPPPfzoRz8iKyuL7u5uHnjgAebMmZPp9QkEQP/Zp4M1nzFwrouvvIaqiWX+ZJNeosU2Y1lZhuZLKGkoE3WJ0XB19ODMtgXnHcKJxJOAKzRwvmeffhpX/QGsNhsnTanimY0bueuO2/z3IDeX7t7YXn9uzNBN16rKSEgnJhCHCLRVVbAoMl0ev4WeTElBOjf2yEQtq6r4Y6mmidybUStLEtYcmRybhWLnqfzjXx+NWCEZy91xEooxrlq1ik8++YT58+fz2c9+ljPPPJNPPvmEVatWZXp9AgHQdzJGvIbfg0Gisc1YzcdVizVqvC4yiWewrrv2YC1dHh8mJobp1yaTE67QAIsWzifX6WTKtJOZNes0CvLz2fvJJ6y++x4+/mRXv3HaAIFNtyzfQa7NApwQ5FBhDjTnDjT9TmSaR7RzpCMmGBrPNAGL6s9e9Wk6Vovqf0iS/HWUFkVmXLaV6dNGbplb5JSQUFf2aCchi7G4uJjnn3+exsZGjhw5QnFxMRMmpH9m3IEDB7jttttoa2vD6XSybt06Kisr034ewchjuAz1DZBobDOWq7e6ujqhsorBum6bp41//OujuFmp8da0++OdXHL5Fbz67OP+4/VzvtBN12pRABPTlMIsRQnQTL9K+8WRpKyVdMcEI+OZRm+to2SaYeU+BiayBIYpUV6Qw8GEjj78GMvdcWIKY2B6BhAcL1VSUkJJSUnYa7KcvpGOK1euZMmSJVx44YX8/ve/56677mLjxo1pO75g5DJchvoGSCa2GcvVm4iwWuzZfPxxDYahY7XaKJlQiiQrGbnuQOJJ6EYeTeDS8bMI3XSV3tpHUwIpxI1qmr1fStDt8XH4aOym3/2dA/rGBFNxDYbGHkNdjUHME+uWwnV+xDGWu+PEFMZ58+bx3nvvAXDKKaf0aSCe7rFTx44d4+OPP+app54C4Otf/zpr1qyhtbWVgoKCtJxDMHJJR/ZpuhlIbDMRYd28ZRvdx1vx+bxYVAu6z8fBA/ux5zr5zvKhu+50/CxCN91AAb1pnJiZKEkSBv4mAabpn4qRbFeWeDHBgGvQYVMHlFXa6fExzpGDiYlk9sYbOeGKbmjt7OcIw5ex3B0npjBu3rw5+PWbb76Z8YU0NjZSUlIS7CiiKArFxcU0NjYKYRQMavbpYNGfsFZXV1NWUoy3IJ+W5iZ8Xg+yasGW40z5ugeajBLIkHV3d3PoeBuyamHCpClhP4tEeteGbrr4dOw2NVgFaRoAJpphIEv+WYddntizH2NdV+TGHowJ+vx1qJIsk2VV/WKWpAUZai1qPg2PF3KzrBi9I6Z0w6TD7Ys6wWOkMJa748QUxtLS0uDX5eXlYd9zu93IsozVas3cygSCCAYz+3Q4EIjlyfKJuZSGYaYcXxxolmFohuyMqilBSzFU+BLNoo3cdFvauzGAiQU5/tpFn4ba66TSDJMjx7tSuq5YMUHwW6GmCT7NiFknGYs+8UsJOnq8SLJEl0cbNSIyVrvjJBQgXLduHR988AEAb7/9NgsWLODMM8/krbfeSttCSktLaW5uRtf9T3O6rtPS0hIm0ALBWGKgXX0i6S/LcPOWbSy7+SdcfOU1LLv5J306+Ty38Rm6XUc5VLuP3bs/wTT0PhmyyWTRBjbdndvf5UDNNg7WbOOtv7/P8W4P9DbvDm3rlsh1ybJMtt1KYW4WZ807DSB4jne2f0iXRwtmlQbcqt6QAcKB5BLwC27p9NmcMu8zlE6fHZbNGjUxRYIuj8bO7e/6e8OOcFEcyyQkjJs2bWLatGkAPPzww9x333088sgjPPjgg2lbyPjx45k5cyavv/46AK+//jozZ84UblTBmCVWqUeg/VuyxMsyDFh6qreLqollqN4uNj6xISiOm7ds49jhfUimic1ixdA16g7VYuhaWIlJrLZx/fWSDZBMi77I6wpM5wg0Ksi2qmHlGTZPGx/uO4ykyNhtVv9YLO8JyxQI64kar9TD1dGDbphh6xgriSljgYSEsaenh6ysLFwuF4cPH+a8887js5/9LPX19WldzE9/+lOee+45zjvvPJ577jlRJykYs3hsTr5zzXfY/s9/8vAjj1A2+aSE+57GIt5m3p+lV11djcVqw5QAyZ8DoCgKDfV1YRZsMlZuLIss1JJMxPIKXJeld26iBEEXaahF7LE5Oa2qAlM3cHu8uL0aDpsFrbduM9Q67c+6jqxpTMSyTeTaBcODhOoYKysr+cMf/sChQ4c466yzAGhtbcVut6d1MVVVVbzyyitpPaZg8KjZvJnXqqupd7koz8/nosWLmbVo0VAva8QRGTM7c948Hn3ssaQ7jkQmpNS3HKUgNytqlmF/9ZKulkZKJpTR2HAIQ9fRNB+G7i/ZqppxavAziWasprOrSn3LUaZMcGJRZTBBN83g6KrQurvIuCCmSY/XhyRLdETEBfur4RtIYspY7igzUkhIGFeuXMndd9+NxWLhZz/7GQBbt24NiqRAULN5Mw9v3Eh2QQFFVVW0tbfz8MaNXA+DLo4jXaCjFqYnmBQSINrmm2OVePTxJ1hwxqlMmTKFNrfOUVebXwj6qU3MLy5F9XbhzC/kWHODvzxBlpAUCzveeYvNM6aFJUf1lz2cjmsMXmdVBW6vhkWR/TWKErR7/C5SLcS9GU3sME26vDo7t78LnEguSaSGL9XElHRduyBzJCSMp59+Oi+99FLYaxdccAEXXHBBRhYlGHm8Vl1NdkEBDqffJRT497Xq6gGJUrIilw6BHmphTbTjSLx1Rm6+He3tuI4285m5p3P7bbeHWHJ+12x/ll7g+93tbdjtWZgSaLpO2cTJWC2WsB6viWQPp6urSvA6JXB7NexWFdMEqyrj9WlhdXfJFKxnsoZvLHeUGSkkPKh4//79fPLJJ3R3d4e9fvHFF6d9UYKRR73LRVHEBJbsvDzq9+1L+ZipiFyqAh0QmT1NTbT7NIqLCikZIss3kQ28v3sTufm2NDeiaz4qJ08KxhDhRNPy/iy9wL+/umcVYGKx2ikrnUjx+IKUSkjS1VUl9DpNw/BbjqqCRZWpbQmfcJGM2GWyhm8sd5QZKSQkjL/5zW94+OGHmTFjRlhcUZIkIYwCAMrz82lrbw8KEUB3ezvl+fkpHzMVkUtFoENFpsc0MRSZZpeLLIeDcb3rH0zLN5ENvL97E7n5ej1u8vLyOFx3ImEusudqf5beooXzqa4+w+9STbIdXDLxzoF0tzENg263P5u1cff7YcdKVuwyVcM3ljvKjBQSEsZnnnmGV155JWxQsUAQykWLF/Nwb1/b7Lw8utvb6W5t5dtLl6Z8zFRELhWBDhUZz4EDWLKyMAyDxuZmxuXnD7rlG28DDwjsuwcP4sjOpsw0g+Idus7IzTdvnBNVVXg+JCSSSk1kpMt11snTWHrFt5g7YwZqZwe1tbU0RyTlFVodnDp/FooiY/g0xo2zUZBbwYf7DlNeXDggiyxZkRkOBetjuaPMSCEhYbTb7UyZMiXTaxGMYGYtWsT1+EWmft8+yvPz+fbSpUlbWaGWlaunB+PQISaETFjpT+RSEehQAbZZLGiahmKx4PZ4Ejpnf6Ri+UbbwEMF1mG14vH5qK2roxIYl58fts7IzbfLUPnZT1ay88MPBtRrNtTlOrmigpV3/oTcnByycnOQxuUxo6SEquPHUbKy8LpcuOrrKTtzDrKqgiwhSxJZFhV6PJQXFw5YpIZSZAbSXm84CLQgNgkJ44033sjatWtZvnw5hYWFYd9L53QNwchm1qJFA3Y3hlpWRkMD9S0tAJRMmpSQyKUi0KFWZmlxMbUNDeiGgc1ioautLWOW795PPmHNsmUJJ/m8Vl2NLsscqqujy+tFN01UVaWhqQlVkvqsM3TzLXG7+eWPVoBzHPv2H+C5F18lv7IqKHTJuHoDLtdPTZiAtbgIU/VvIxKgjBuHkpeH1+PBmp3NxKoqTIt6YsiiLINhYFOVtCWbDIXIiJKL0U1Cwnjbbf7p3KE1humeriEQRFpWEyZOBOB4aytyElZosgIdamWOKyxkQnc3zS4X2aqKU9NSsnxDiebebW5ooM3joU1VE04s2tPURJumoagqWQ4HXrcbr6bRrmnUHTqETZZ5rbcgP/QYJW43J8+bh2m3YyoKZUVFfP6M09m1fTvNpJ7Ja83Px+xt+g+gWK0g+TvQIEmYqoqkRtliZBnZotJ2rD2p+zicECUXo5uEhHEwpmsMNUOdoi+IblmVlJUh9/Twm1dfzdh5I63Mivx8brjiirT9/KO5d1uOHKUkPz8p96rHMECWUSwWAGxZWRg9PeiaxsRJk8jOy4sqapWVlX5R7BUpU1XBbqeyspLmpqaUM3m9LhfW7OwTFqPqtwxN/UTvUbN3zaYknZhNaJjo+shONhkOJRcDnZQiiE1CwhiYrmEYBkePHqW4uDijixpshlNx+lgm1czWdDzUDNQN3N+xI927eRaVkkmTwt7XX5KPXVHp0TV0XUdRFHRdx6dpqJIUV9QiLTsAU1Gw5udDUxP1LheWvDwO1dTwmYULufHGG6mqquLQgQMUud19kmkC1NbWcvL48dBriaLrIMvoPh/HXS4aW1rIczqZUFKCr8dNdr4TWZIwNY2P/rkDmzf2xIyBUuJ2U1lZiTU/H6/LFTUpaCAMdcmFcOVmloSEsb29nVWrVvHGG2+gqirvv/8+b775Jh988AE33XRTpteYcTJVnC5IjlQSZ4bDQ00iwhwpvGuWLUv6IWBqSTGH3W7aOjtxezzYrVZURcER4a6MFNhIyw5A0nW8Lpf//YrK7oYGvnzuudz7i1/Q43ZT39xMSVkZ04qLYfv2qKLSbLez7+mnKZk3j6KKCjydnUyuqMCjadTW1ZGTkwO6ziMPPMCcBQuYO2MmSlcntbW1HB2ASPV3vyNdx9bsbL+Ax7iOVBjqkgvhys0sCWXOrFy5kpycHN566y0svW6cOXPm8Kc//Smjixss6l0usvPywl7LzsujvnfjEGSGms2bWbNsGd+9+GLWLFsGwPVLl+LUNI7s24dT07i+n/he6EONJMs4nE6yCwqCsbbBuIaHN248EStUVR7euJGakEHf0bho8WK6W1vpamvDNIxgks9FcSZnXLR4MYrXS8WECZwxaxYVEyZgRSI/NzfsfZECW1tbi+R2I2kamCaSpiG53dTW1gIgyRJIEtcvX06P201HRweGYdDa1kaHrlMZkhUcee3/9xe/4LvXX8+VV1zBkiVLuP7aa/mopobCggLaWlp4cN06/vP117nphhu48Ivn8PempgGJU+j9/uo3vsGaDRv4+rXXMT0nhxK3G4hwHffGOs1e13G6SGUKSDqJ58oVDJyELMZ3332XLVu2YLFYkHqzywoKCjh27FhGFzdYZKI4XdCX0Cf9bEWl1eNhfOmEcEtv6VLufOyxhI+ZiY47yZCqtyGV7Nlon7n6C2fzxo4ddLW1xbSym+122L49pmuxy+djyqRJnHTSSTS1tIBpIssyhmmy99AhTi4ti7qeaNf+/s6dfPXrX+f0M85g/hlncNmSJfxw0iTqDx1iw0O/TObWxj3n5xcu5KZbb8Wjabjaj5OXk8PJkyfD9u39uo7TxVCWXAy1K3e0k5Aw5ubm4nK5wmKLDQ0NFBUVZWxhg0kmitMF4US6PD+qqcFrGIwzzaClB8m7r4f6oWYgwpxKXDPaZ6YGHjjiCGyz3U5zU9MJYQix2srz82mTZQ7W1lJYUkJ3V5c/61yWycnJYeee3f4yC8Jjdyt//Wtee+UVtn/4Ydi1A8yaPj1EuNrJLy5mzf2/oGPPngFZjIH7fdmSJXg0jZ6eHhRFoa29Eh5ZdQAAIABJREFUHbO8nMrKyn5dx5lksJL4htqVO9pJyJX6zW9+kxtuuIH//d//xTAMduzYwa233spll12W6fUNCrMWLUrahSdIjkiXp6brqBYLjc3Nwfek4r5OxSUZj0j3bn8u0fLewvpQmhsacPX0JHyMgTJr0SLufOwxfvPqq9z52GNJ/94G7uFDDz1Elt1OtsOBCeSPG0eWzcav/+NXwInYnbW8HDM3l/KKCm685RbmnXZa8Fjd7e1UjXNy+dKl9Hg89PT0oOs6nZ2d5DocA3ZnBu532aRJ9PT4rSNd17FZrUGrsD/XcaZI1a2eCkPtyh3tJGQxLlu2DJvNxurVq9E0jTvuuINLL72Uq666KtPrGzQymZWYKqOphCTSsrJZLPhME7fXC0BbSwt1TU3ous6aZcsSvtb+XJLJ3MNEE3miuYTBL+zNvU0JygsKKJo0aURkOAfu4Zpf/5qbV6xgxYoVVFVV0VxXxy+feIJDH9UAfcs+nHl5tLS28s3LLmPbv/4V9LRcv3Qpc2fOpPboUXxeLzarldLSUhx5ef5QzADcmQHvzuHaWorKyujs6sLQdUpLS4NWYX+u40wx2El8ontO5khIGCVJ4qqrrhpVQjjcGQ7Zlukk0uVZWlzM/vp6rBYLruZmDjQ0gCQxZfJk2mQ5qWuN9VCT7D1MZGOLPGZ3eztmQwO0tXHkyBGO9/RQXlAQbGM3UjKcZy1axJ3Awxs3csuKFWEhhet7QwqRsbtAj9bcrCyORDyUKJ2dzKyqCndnatqA3ZkBEa9+4gm+t3IleQ4Hzrw8nLm5J6xCuz2u6zhTDHW8W5A+Eh47VV9fH3Xs1Pnnn5/2RQlGXwlJZBzXYrUy3majMDubvQ0NWC0WKsrLg5stDPxak72HiWxssY7p1DTuf+wxvnvxxRQlWZ84XOjP+o4Wu3Pm5pLd3t6nAUNkjaOk62HCNdB1zgLc+/f7rUJFwVtfPyhWYTyGOt4tSB8JCeOGDRv49a9/TVVVVZ+xU0IYM8NQPn1mwoUbddO99lpmLVrkF5OqKn8rsV7Sca3J3sNENrb+jjnSN8d4IYVkxG4w3JlDYRXGQyTxjR4SEsYnn3yS//zP/2Tq1KmZXo+gl4FssAMRtky6cGNtupkSk2SPm8jG1t8xR8PmGOv3J1mxG27ClWnSNWFGMPQkJIxOpzPYFk4wOKS6wQ5U2IbChZsuMYnc0E+rquKNHTv6HHfhnDlRp1okksjT1tnJJ8eOYbNYKC8pwWqzha111qJFnPfhh7y4ZQsdPh+5FguXL1w4YjbH/n5/BlPsMt3WLRMMRRLfSLxPw52EhPGOO+7gzjvv5KqrrmL8+PFh3ysri178KxgYqT59DlTYhsKFm44n7Wgb+hs7dnDenDl8uG8fez7+GI9hIAFPvf02Jfn5lETZ+BNJ5KnKy6OuoYF99fXMcOaHlfbUbN7MGzt2BJt6d/euY+rmzQPeMAcjS3m4xLZjtXXb9/TTPPjoo6MiUzsdDEb7u7FIQsLo8/l45513eP3118NeF2OnMksqT58DFbahipEN9Ek71ob+4b59QYu0qKCAQ42NmEBTWxv27GycvU0r+tv4I4/vLCigqa6O+tZW1j/1FOXV1Vy0eHHGhGWwspSHS2ZltIkgnbpO3uzZSY3qGu30NzlFkBoJCeOqVav44Q9/yNe+9rWw5BvB8GOgwpbJGFkqFk+0zwB9Xou3oYeKlefQIaw2G7rPR2NLC87i4rCNP9YaI49/3OWiqbUVQ9OYHrJJd/e4mTRzRtR1DIRUBDeV+z1ckoeitXU73NTE5JNOGnJrdjgxWO3vxhoJdb7RdZ1vfOMbOBwOFEUJ+08wvBhoJ5hMdQFKpStItM/84rHHue/RR/scx2Gx9OlAE9jQQ5vE261WdF1HVlU8Pl/Y+0LPZx0/nh1NTdzyH//BDy+9rM/xAx17sqzWsOblbl2LuY6BkGyj+8h7d9jt5ie/Ws/S88+P240n3Z2EUsXrciGFzHUEUFSVxrq6sNfGerP/aPdpsNrfjWYSEsb/83/+D48++iimaWZ6PYIBkg5hG2iLsWikMgUj2mfadY0OTetzHNMwY27ooW3bSktKMHQdt9uNzzB477332HvgAKdVVQXPp5kmB+vrMWUZq9VKbWcHR7u7OdbYFDx+t9sNhkFpSP/g7Lw8bLKcEWGJ1nounuCG3rv248dpdrkwFJke04z7UDJc2iNGa+uWm5XFM088Efa+kVQKkwmGqv3daCchV+qzzz7L0aNH2bBhA84QFwvA22+/nYl1CQbAcGxvl0rsKtpnfLqOZBh9jnPkyBGWf/vbMRN4Au7hvHHjGHfsGE1tbSiSTFZWFvm5ubyxY0fQDbpr927kXo+IKcv4enoonDQJs9WFU9Oo37cPhywzLi8vGKME/yY9bcKEYKwxnSn7ybq4Q+9dY3MzsqIgyzKenp5+XZDD4fcnWmnI+5s28fZ//RfZBQUjthQm3QxV+7vRTkLCeN9992V6HYJRTiqxq2ifsSgKSJJ/QnxzM26vF1WWOSknJ+aGHsh6fXLjs3xw4AAeXceuKFSWlweFrautjdZDh+hub8ft9WK1+TtPGpqGzWIJiu8DvSOxAq7KaOOeBiossWKDiWTuBj7b2NlJy8cfU1FeHrwe3efD1jtPdSR044ksDck580yuX7pU1AlGMNbqRQeDhIRxwYIFGV3EqlWrePfdd7FarWRnZ/PjH/+Y00I69gtGPqkk9UT7TJ6i4ta97Dt4ENViQZZlvD4fR7u7qemnJKIHk6knncTugweRFYXahgYqIZiAE3CDqrKMpmlIpomuaVSUlfUR8UwVc/eXfRrv+KGfnVRWxoGGBv99UhS8Hg8YBhW95VWZcEEORjnJcLBmBaOfmML4yCOP8L3vfQ+AX/4y9oDRG2+8ccCL+PznP88dd9yBxWLhr3/9KzfddBN/+ctfBnxcwfAhXYN5v73sOzy58VlqOzvQvF5sFguTysuxWK1xsxNDY25ZTU34DMOfzNGbmRrqBn3queeCRfyTS0uxWK1RRTxdm3SooLh6ehgXJ/s0nviEXqMDfzlVXVMTXp8PRZIoyc9nXGFhMO6ZThdkquUkNZs3c6ihga9ffTVTpkxBdrk4fvSocAUKhpSYwtgUkurblOG033POOSf49ezZs2lqasIwDGQ5odwgwQghXYN5u596ilNnzQrrrWoaRsLxytKSEmrr6pBkGY/XGyYUsxYt4v4IASrOoMsuUlDq3n+fntZWshyOYEP1gNuzP/GJjMk6i4sZV1jIkX374sZf00Gq5SQ79u3jznvvxadpHGlrw263U3HaafDhh0IcBUNGTGFctWpV8Ouf//zng7IYgOeff54vfOELQhQFMRlovHJcfj6VwOH6ekz8kzEihWKwXHaRgpJlteIxDBqbm4PCGLi2/sQn3n1J5npScYmmklz1WnU1azZswKfr9LjdKIqC2+OhuaNDFKgLhpSEYox79+5l27ZtHD9+nHHjxjF//vykGopfdNFFNDQ0RP3e//zP/wTrITdv3symTZt4/vnnEz62YOSSakwqNPboMwzqGhpwe73McObHjDNGxitVSaIkO3tIShFCiRSU0uJiahsa6Ha7/WUhIbHY9U89FVd80tGcIVWXaCoPK/UuFxWVlbhCylAURaGtvZ3KyZNFgbpgyIgrjKZpcscdd1BdXc2ECRMoLi6mubmZlpYWLrzwQu6++27/RO5+eO211/p9z5///GcefPBBnn76aQoLCxO/AsGIZCAtzgKxx6eee459vbHAqvJyJJst5jGG6+SDSEFxFhczobub4+3tfYb/lldXxxWfdFxjqi3tUhHl8vx8DtfWUlReTk9PD+BvJuLMy4tboD4YST6CsU1cYfztb3/LP/7xD377299y+umnB1//4IMPuPnmm3nppZe4/PLLB7yIv/71r/z85z/nqaeeYuLEiQM+nmD4M9CeorMWLcJZXc0MpzNMKOIdI9SdGNhc1z/11JBsroHz721uoc3r8Tc1nzSJ7vZ2FMPgzu9/v896EhGfgbqAU+2VmoooX7R4MU+sX8+d994LdjudXV3YbTZKcnOp/fDDqGUHg9UzVjC2iSuMv//97/n/7d19UFN3ugfwb04SQI0SRHmJVLFq1Q7t1m3VXbfrW7urLSgjU0ftUKddlx07i7rtOK1rS90t06pbRfvOLnesW7Tb23EH3B1m8GVlXHZqndK6WsrV7s36AiTyIgktcAkkOfcPTTYiCQmc5JyTfD//ScjJY8Q8PL/z/J7fyy+/fFtSBID7778f27dvx+9//3tJEuOvf/1r6PV6bN682fu1gwcPIimGJ1pEOymGVQ/3GlJ8uEp15uXk2bMQb7Ggta0dvV1dmJGW5jehRKLqHcms1FCTclZ2NlBVhfd27LitK7UxQOONUk7/oOgWMDGazWbMnTt30Mfmzp2LF154QZIgPvvsM0muQ+ohxbDq4V5jpB+uUp95mZaRgbEGA4xOJ4puDRDwJ9xNQeEYIh/ol4is7GxkAXAD+N/GxptPCNCNqpTTPyi6BWz9dLlcMBgMgz5mMBjgHjCaiyhYUgyrHu41Qh3I7VFfVYXiggK89PbbaOnqQn9fX9BzX6V4/UiQelbqcIbHBxLqzFii4QhYMTqdTnz22Wd+h4e7Bkx1JwqWFMuCwVxjsGplOJWmb5XoFgRoBOGOyTlDHV3loZSjnfyRsiqVeukznMeiEXkETIzJycnYvn2738fHjx8veUAUO6T4AA50jcGWPPeU/RfiRBFNvf+HhLg4ZJhM0N8aBRfowzWYyTm+R1cFWmZV8oe71B2fUi99KrW7mKJLwMR46tSpSMVBJLmB1YpTFNHu6EWcIGDapElobmmB+fJlzEpOHnK5MNjJOcFUSEr9cA9Hx2c4qmPOS6VwC2qDP0mPe7ECk+L9GVitWFtaoNPr4ezrQ1JqKpJSU9Ftt8PodIZ0sn2gyTlDbcL3UOKHezg6PpVcHRP5w7lrMpC6ISHaSPX+DGzU6O3rA0TRe/QSEHzTy8BGH8/knNc2bfJ2khYXFMDa1YWvGxrQ6XNNJd0/DCQcTUFKOfiYKBSsGGXAvViBSfX+3DEG7tYRVZMnTfJ+T8+332KMXo/igoKA1Wmg5U/fJcjxY8ei2W7H/5jNGKPXY2JiIrRut2IrpIEne7ivXUNaZqb3cSmSuhKrY6JAmBhlwL1YgUn1/gxMZlMNBrT39EAfF+edQ3rDeh2i6AaMxiHvq/n7gPck8v6+Ptzo6kK8Xo8+pxPd/f0QbXY8vXiRIhPDwHuKbosFza2tAOCdwsNlT4pFTIwyUHq7vtykfH8GJjNvhXSr6hPj46EZnzSi6tSTyC82NECr00Gr1yMeQJ/DgWmTJ+MrsxlD7c6U457zYIMGAKCzowOCgpqCiCKN9xhlIMXm9mjm7/25b9o0FBcUYOMTT6C4oECSe7LtPd0jvq/muZfp6O+HoLv5u6bL5UJCXFzQgwPkuOc82D3FVJMJSaNGofTIERSVlTEpUkxiYpQBGxICG+z9WTZnDo6dOzei5DFYAvq234mWa9e832NvbcXX9fWwdnUFnXw9iVyn1cLldMLlcsHtciE9NTWowQHF772HRrsd15qa0NneHvIkneHiFBmiwXEpVSbhbkhQ+3aQge9PcUHBiBtyBmvqSZk4AS1tbRhrNKLP4cBliwXQaHD3lCmwC0JQ+/g89zIPfFiOi3YbEuLiMGXSJOg0moD36DyJutvtRvyoUXC6XN5pOokTJoT9nrMSt1Ko/eeWogMrxigUjdtBpNhK4G/p0BgfD6PTiWsWC+L0ekybMgXGWwk02MotKzsbJf/9MX63eTPmpKWh78aNIVcCPIl6dEIC3KIIUaNBn9uNS42N+LqhAWN8tpUMh2e2q7+lZ6WtXETjzy2pEyvGKBSN20HCeRrH9JQUFJWVYeMTT2DitGnQCP/5fTHUbthQVgI8TTvpogjztWvodzohAoAooq+/H+09Paivqgo48s5fdRXsFBslbaWIxp9bUidWjFFIyac3DFckTuOI9D03z+slJiVBr9VCo9EAogitRoO7J03CBJPJb7U6VHXlm2SGcwKIHKLx55bUiYkxCkVjU4UUy35DXSPS3cK+r+e81cWaoNNhxl13/efEDj9JYajEp8YkE40/t6ROXEqNQkpsqpBCuE/jGO4xVsONyff1BLcbgkaDySYTjCkpAAInhaGGIKhxr2y0/tyS+jAxRiGlnt4QDlJ3MYZ6jNVIT5/wvJ7n2r5TeQIlhaESnxqTTCz93JKyMTFGKSU1VYRLOBJVIOFsDgk1KQyV+NSaZGLh55aUj4mRVCuSXYz1VVWoa2yEWxAw6vp1pKemIjEpSdIZt6EkhWASH5MM0fAwMZJqRWoYu6cy1Wq10AgC+t1uXGlqQiYAnUYj2307Jj6i8GBiJNWKVIOJpzLNMBhwxWKBVqeDRhDQ2NyM1NGjFX3fDuA0GaJQcbsGqVaktld4tj4YU1KQaTJBB8Dd1wdXf7/iZ9xymgxR6JgYSbUiNdLMd3+dMSUFs7OyMGPqVDx0112KToqAOjf6E8mNS6mkapG4z6bGrQ8ePBSbKHRMjERDCNQBqvT7d2rc6E8kNyZGigrhTlCDVaaR3kc5HGqudonkwsRIqidXgjrwYTlaurrg7OxEvF6P9JQU7/274ZyI4c9Ikr5aN/oTyYmJkVTDX4KQ47ii+qoqXLTbEBcXB31cHJxOJ65YLJiSno7mAYOwfZ8TagKXIulzvyNRaBTVlXr27FnMnj0bhw4dkjsUUphA2w7kOEmiorISCXFx0AgCNBoNtHo9tDodmlta/N6/G06HKLtKiSJPMYmxq6sLe/bswcKFC+UOhRQoUIKQ47iiZpsNGSYT3C4XXC4XAEDUaODo7/e7jzLUBO4ZQ/fN1au4ePEiOm99n9KPjyJSO8Ukxl27dmHDhg1IYrcc+aivqkJxQQHOXL2Ka42NsLe2eh/zJIhIn6MI3Oz21AsCMjMyoBcE9DkcEEQRs5KT/S5bhpLAfcfQabVa7xi6TpuNXaVEYaaIxHj69Gl89913WL58udyhkIL4Lp+OGT0aDlHEFYvFmxw9CSJSG/19eZKxTqPBzHvuwYzJk2+Oh8vPH/I5wSRw7xi6tDS4nU7A7faOoQt30ieKdRFpvlm1ahUsFsugj1VXV2Pv3r344IMPIhEKqYjv8qlJFHGlqQmiIMDa0gJ9XNwdxyxFssFkON2eoTzHszFfIwjIBGBtbYWjrw8ioPgxdERqF5HEWFFR4fexuro6tLW1YfXq1QAAm82Gmpoa2O12FBYWRiI8UijfqS2JSUnIBGC5fh3dPT0wOp2ybzsYTjIO9jm+G/ONKSkwpqSg226H0emU5e+s9EEGRFKSfbvGQw89hDNnznj/vG3bNmRlZSE/wJIUxYaBU1sSk5Kg02hgdDpRVFYmc3TDF0ySUdLGfDUMMiCSkiLuMRINRo6mmnAL9rQLOe6b+qOGLSOeJq2NTzyB4oICnh5CIyJ7xTjQrl275A6BFCIap7aEMoxAKRvzlT6InBUtSU1xiZHIl1KSg1SUnmQGo/RB5HJMPqLoxqVUogiSYxjBSCl9SVuOyUcU3ZgYiSJI6UlmMEq63zkYNf6yQcrGpVQKCdv2R0at902VvKStpA5eig5MjBQ0NjlIQ8lJRo3U+ssGKRcTIwUtWpoclFb1Ki0eNeIvGyQl3mOkoEVDk0Ow+whjNR4iYmKkEERDk0Oom9XDvXFcDZvniWINEyMFTY0dlQOFUvVGopqLhio8XDjNhuTCxEhBU3rbfjBCqXojUc1FQxUeDlxiJjmx+YZCovYmh1Ba+yMxpYZbDQYXLY1epE6sGCmmhFL1RqKai4YqPBy4xExyYsVIMSfYqjdS1Zzaq/BwUPp8VopurBiJ/IjWak4NTS3R0OhF6sWKkSiAaKvm1DK9iNNsSE5MjEQxRE1NLdH2SwmpB5dSiWIIm1qIhsbESBRDuG+SaGhcSqWYFYvDu7lvkmhorBgpJsXqZJVo7bQlkhIrRopJampCkRqbWogCY8VIMYlNKETkDxMjxSQ2oRCRP0yMFJM4WYWI/GFipJjEJhQi8ofNNxSz2IRCRINhYiTZxeJ+QiJSLi6lkqxidT8hESkXEyPJync/oUYQMMZoxOjx41FRWSl3aEQUoxSzlFpeXo7Dhw9Dr9dDEAQcPXpU7pAoApptNkycNu22r40eNw7NZrNMERFRrFNEYjx+/Diqq6tx5MgRGAwGtLe3yx0SRQhPaicipVHEUuqBAwdQWFgIg8EAAJgwYYLMEVGkcD8hESmNIhKj2WzG+fPnsXbtWuTl5eGTTz6ROySKEO4nJCKlichS6qpVq2CxWAZ97NNPP4XL5YLVasVHH30Em82GdevWYerUqZg7d24kwiOZcT8hESlJRBJjRUVFwMdNJhNycnIgCAKSk5OxYMECXLhwgYmRiIgiThFLqTk5OaitrQUA9PT04IsvvsCsWbNkjoqIiGKRIhLj008/DavViuzsbKxevRorVqzAj370I7nDIiKiGKSI7RoJCQl444035A6DiIhIGYlxpFwuFwDAxkNmiYhoCJ5c4ckdA0VFYmxrawMA7Nu3T+ZIiIhILdra2jBlypQ7vq4RRVGUIR5J9fb2or6+HhMnToRWq5U7HCIiUjCXy4W2tjZkZWUhISHhjsejIjESERFJRRFdqURERErBxEhEROSDiZGIiMgHEyMREZEPJkYiIiIfTIxEREQ+mBiJiIh8MDEqSHl5OZYvX44VK1YgNzdX7nAkd/bsWcyePRuHDh2SOxRJ/fa3v8Xy5cuxcuVKrF27Fl999ZXcIY3I5cuXsWbNGixbtgxr1qzBlStX5A5JMjabDQUFBVi2bBlWrFiBwsJCdHR0yB2W5N555x3MnDkT33zzjdyhSMrhcGDHjh346U9/ihUrVqCoqCgsrxMVI+GiwfHjx1FdXY0jR47AYDCgvb1d7pAk1dXVhT179mDhwoVyhyK5hQsXYvv27dDr9aipqcFzzz2HkydPyh3WsO3YsQNPPvkkcnNzcfToUbzyyiv48MMP5Q5LEhqNBj//+c8xf/58AMDu3buxZ88evP766zJHJp2vv/4a//znPzFp0iS5Q5HcG2+8gfj4eBw7dgwajSZsn5OsGBXiwIEDKCwshMFgAABMmDBB5oiktWvXLmzYsAFJSUlyhyK5JUuWQK/XAwAeeOABXL9+HW63W+aohufGjRtoaGhATk4OgJtnpTY0NERNVWU0Gr1JEbj572WxWGSMSFp9fX149dVX8Zvf/EbuUCTX3d2NyspKbNmyBRqNBkD4PieZGBXCbDbj/PnzWLt2LfLy8vDJJ5/IHZJkTp8+je+++w7Lly+XO5SwO3z4MBYvXgxBUOd/LavVitTUVO/MYa1Wi5SUFFitVpkjk57b7caf/vQnLF26VO5QJPPmm29i5cqVyMjIkDsUyTU2NsJoNOKdd95BXl4ennrqKdTV1YXltbiUGiGrVq3y+5vpp59+CpfLBavVio8++gg2mw3r1q3D1KlTMXfu3AhHGrpAf7fq6mrs3bsXH3zwQYSjks5Q/3aeJFJVVYW//vWvOHz4cCTDo2EqLi7G6NGjkZ+fL3cokjh37hzq6+uxdetWuUMJC5fLhcbGRtx777148cUXcf78eWzcuBEnTpzwrrRJhYkxQioqKgI+bjKZkJOTA0EQkJycjAULFuDChQuqSIyB/m51dXVoa2vD6tWrAdxsfqipqYHdbkdhYWGkQhyRof7tAODEiRPYt28fDh48qOpl8PT0dLS0tMDlckGr1cLlcqG1tRXp6elyhyap3bt34+rVqygtLVVtdT/Q559/DrPZjEceeQQAcP36dWzYsAE7d+7Eww8/LHN0I5eeng6dTudd5v/e976HpKQkXL58Gffdd5+0LyaSIrz//vvi3r17RVEUxe7ubjEnJ0f8xz/+IXNU0nvxxRfF8vJyucOQ1KlTp8QlS5aIV65ckTsUSeTn54uVlZWiKIpiZWWlmJ+fL3NE0tq7d6+Yn58v9vT0yB1KWC1ZskS8dOmS3GFI6plnnhFra2tFURTFf//73+K8efPEzs5OyV+Hx04pRG9vL4qKitDQ0AAAyM3NxS9+8QuZo5Letm3bkJWVFTXLVwDwgx/8AHq9HuPHj/d+7eDBg6ptNDKbzdi2bRu+/fZbjBs3Drt378bdd98td1iS+Ne//oWcnBxkZmZ6z+HLyMjAu+++K3Nk0lu6dClKS0txzz33yB2KZBobG7F9+3bY7XbodDr86le/wqJFiyR/HSZGIiIiH9GxuE5ERCQRJkYiIiIfTIxEREQ+mBiJiIh8MDESERH5YGIkirBXXnlFUdsDLBYL5syZA5fLJXcoRIrA7RpEElq6dCna29uh1Wqh1Woxffp05ObmYs2aNVEzYSUU+/fvx9/+9jeYzWY8++yz2LRpk9whEQ0p9v6nEoVZaWkpzp07h5qaGhQUFKCsrAwvvfSS3GHJYsqUKdi6dWtYNmEThQsTI1GYjB07Fo888gj279+PiooK76Gx27Ztw759+wDcPLx54cKFKCsrww9/+EM8/PDDOHnyJE6fPo1ly5Zh3rx5KC0t9V7T7XbjD3/4Ax599FHMnz8fW7Zsgd1uBwA0NTVh5syZqKiowOLFizF//ny8//773udeuHABeXl5+P73v48FCxZg586dtz3P6XQCAFpaWrBx40bMmzcPP/nJT2476eXtt9/Gli1b8MILL2DOnDnIzs4OeDDzqlWrsGjRIowZM0aid5Uo/JgYicLs/vvvR1pamt8jctrb2+FwOPD3v/8dmzdvxssvv4y//OUv+POf/4zDhw/jvffeQ2NjIwCgvLwcJ0+exKFDh1BbW4vExES8+uqrt13viy++QHV1Nf74xz/i3XffhdlsBgC89tprWL9+Pb788kucOHECjz322KDxPP9tkppKAAACxUlEQVT880hLS0NtbS3eeustlJSU4MyZM97HT506hezsbNTV1WHp0qUoLi6W4m0iUgwmRqIISElJQWdn56CP6XQ6PPvss9Dr9Xj88cdhs9mwfv16GAwGzJgxA9OnT8elS5cAAB9//DGee+45pKWlIS4uDoWFhTh27Ji32gOAwsJCJCQkYNasWZg1axYuXrzofZ1r166ho6MDY8aMwQMPPHBHLFarFV9++SW2bt2K+Ph4zJ49G6tXr8bRo0e93/Pggw9i0aJF0Gq1yM3N9V6fKFrw2CmiCGhpaUFiYuKgjxmNRu+Zjp7B1snJyd7H4+Pj0d3dDeBmB+kvf/nL2xp5BEHAjRs3vH/2PfZq1KhR6OnpAXCzYnzrrbfw2GOPISMjA4WFhViyZMltsbS2tiIxMfG28+1MJhPq6+sHvX5CQgIcDgecTid0On6cUHTgTzJRmF24cAEtLS148MEHR3yttLQ0vP7664Neq6mpKeBzMzMzUVJSArfbjePHj2Pz5s04e/bsbd/jqWy7urq8ydFqtSI1NXXEsROpBZdSicKkq6sLNTU1eP7557Fy5UrMnDlzxNdct24d9u/fj+bmZgBAR0cHTp48GdRzjx49io6ODgiCgHHjxgHAHVtI0tPTMWfOHJSUlMDhcODixYs4cuQIVq5cOax4+/v74XA4IIoinE4nHA4H90uS4rFiJJLYxo0bodVqIQgCpk+fjmeeeQZr166V5Nrr16+HKIr42c9+htbWViQnJ+Pxxx/Ho48+OuRza2trsWvXLvT29sJkMmHfvn3epVtfJSUl2LFjB3784x9j3Lhx2LRpExYsWDCseIuKilBRUeH9c2lpKXbu3Im8vLxhXY8oErjBn4iIyAeXUomIiHwwMRIREflgYiQiIvLBxEhEROSDiZGIiMgHEyMREZEPJkYiIiIfTIxEREQ+mBiJiIh8/D93zQyQYtFA2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 465.994x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Next, the Gaussian mixture is instantiated and ..\n",
    "n_components = 2\n",
    "model = GaussianMixture(n_components, d)\n",
    "model.fit(data)\n",
    "\n",
    "# .. used to predict the data points as they where shifted\n",
    "y = model.predict(data)\n",
    "\n",
    "plot(data,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "id": "5182a6e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fc3fbd5ee80>"
      ]
     },
     "execution_count": 186,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAHUCAYAAAD4NWZRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df5BcZZ3v8U93z0xiJhDIDzIxWQLZSMxevLC3wrKLWRSsXbIYSEFd1kiw3DXKukotUrVVslLyQ0shJeqlCjVLXYFVRlnXW8hPUbaMKAYt4mrciAKBEJHMhCRASIJJmJ6+f3Q66emc031+POc8z3PO+1Vlxenp6Xn6BznffL/f5/tUGo1GQwAAAEitansBAAAARUFgBQAAYAiBFQAAgCEEVgAAAIb02V7A/v37tWnTJs2aNUu1Ws32cgAAAELV63Xt2LFDp556qiZPnnzU960HVps2bdKqVatsLwMAACCy4eFhLVmy5KjbrQdWs2bNktRc4NDQkOXVAAAAhBsdHdWqVasOxy+drAdWrfLf0NCQ5s2bZ3k1AAAAvYW1L9G8DgAAYAiBFQAAgCEEVgAAAIYQWAEAABhCYAUAAGAIgRUAAIAhBFYAAACGEFgBAAAYQmAFAABgCIEVAACAIQRWAAAAhhBYAQAAGEJgBQAAYAiBFQAAgCEEVgBgwu5hafNJ0m+rzT93D9teEQAL+mwvAAC8t3tYGr1carze/Hpsa/NrSZq2yt66AOSOjBUApLXjmiNBVUvj9ebtAEqFwAoA0hr7XbzbARQWgRUApNV3YrzbARQWgRUApDXrM1JlysTbKlOatwMoFQIrAEhr2ipp6Dapb76kSvPPodtoXAdKiF2BAGDCtFUEUgDIWAEAAJhCYAUAZcEQUyBzlAIBoAwYYgrkgowVAJQBQ0yBXBBYAUAZMMQUyAWBFQCUAUNMgVwQWAFAGTDEFMgFgRUAlAFDTIFcEFgBQFlMWyUtfF5663jzz7IFVYybQA4YtwAAKD7GTSAnZKwAAMXHuAnkhMAKAFB8jJtATgisAADFx7gJ5ITACgDKoOyN24ybQE4IrACg6FqN22NbJTWONG6XKbhi3ARywq5AACi6bo3bZQospq0q1/OFFWSsAKDoaNwGckNgBQBFR+M2kBsCKwDwUWcz+shHwpvTadwGckOPFQD4JmiK+O6vHPl+51TxVl/Rjmua5b++E5tBFf1GgHEEVgDgm6Bm9E6dzek0bgO5oBQIAL6J2nROczqQOwIrAPBN1KZzmtOB3BFYAYBvgprRO9GcDlhBYAUAvgmaIj7tH5kqDjiA5nUA8BHN6ICTyFgBAAAYQmAFAABgCIEVAACAIQRWAAAAhhgPrG699VYtWrRITz/9tOmHBgAAcJrRwOrXv/61fvnLX2ru3LkmHxYAAMALxgKrgwcP6lOf+pSuv/56Uw8JAADgFWOB1S233KILL7xQ8+bNM/WQAAAAXjESWP3iF7/Qpk2bdOmll5p4OAAAAC8ZCayeeOIJPfvss3rXu96lc889V6Ojo1q9erUee+wxEw8PAADgBSNH2lx++eW6/PLLD3997rnnau3atTrllFNMPDwAIMzuYWnHNdLY76S+E5sHL3PUDWANZwUCgK92D0ujl0uN15tfj21tfi0RXAGWZBJY/eAHP8jiYQEA7XZccySoamm83rydwAqwgsnrAOCrsd/Fux1A5gisAPhr97C0+STpt9Xmn7uHba8oX30nxrsdQOYIrAD4qdVfNLZVUuNIf1GZgqtZn5EqUybeVpnSvB2AFQRWAPzUrb+oLKatkoZuk/rmS6o0/xy6jf4qwCJ2BQLwE/1FTdNWEUgBDiFjBcBP9BcBcBCBFQA/0V/UW9mb+wELCKwA+In+ou5o7gesoMcKgL/oLwrH8FDACjJWAFBENPcDVhBYAUAR0dwPWEFgBQAuMdVwTnM/YAWBFQC4wmTDOc39gBUEVgDgCtPT5KetkhY+L711vPknQRUjKJA5AisAcEXRGs5dC2IYQYEcEFgBgCuK1HDuYhDD+ZLIAYEVALiiSA3nLgYxRcsIwkkEVgDgiiI1nLsYxBQpIwhnEVgBgEuK0nAeN4jJox+rSBlBOIvACgBgXpwgJmk/VtxgrEgZQTiLswIBAOa1gpUd1zTLf30nNoOqoCAmybmGrWCs9XOtYKz9d4eti0AKGSKwAgBkI2oQk6Qfi0Om4ShKgQAAu5I0lbvYHA+IwAoAYFuSpnJ2+MFRBFYAALuSNJWzww+OIrACgLJx7agZKf6YCXb4wVE0rwNAUewe7r0LL+luOhexww8OImMFAEUQdRaUi0fNAAVCYAUARRA1YGI3HZApAisAcEWa3qeoARO76YBMEVgBgAuSHuvSEjVgKvpuOhcb81EqBFYA4IK0vU9RA6Yi76ZLG5wCBrArEABckLb3Kc7ZfEXdTccxN3AAgRUAuKDvxEOZloDboypqwBQVjflwAKVAAHBB0Xuf8kBjPhxAYAUALihy71NeCE7hAAIrAHBF3GNdTCnKTjqCUziAHisAKLMiHXEj0WcG68hYAUCZccQNYBSBFQCvbB7dqbvXb9RX1z2hu9dv1ObRnbaXlL0sS3VRd9IVpVwIZIzACoA3No/u1GNPbdW+AwclSfsOHNRjT209ElwV8eKf9dDLKDvpGLwJREZgBcAbG557UfMHHtLfzrpQHxg6U38760LNH3hIG557sbgX/6xLdVF20lEuBCIjsALgjdmV+7R02md1TN+oKpWGjukb1dJpn9Xsyn3+X/zDsm1ZD72MspOOwZvwgSMZa3YFAvDGnx37FfVX90+4rb+6X3927FeksdHgH/Lh4t9tZ56Jiey99NpJl8cagDQc2t1KxgqAN6bUgoOnKbVRv6dud8u2uTD00oU1AN04lLEmsALgjUpIkFRpHTjs68W/W6nNhaGXLqwB6MahcjWlQAD+mPWZiel+6Ujw1LrI77im+ZdpK9jy4eLfq9TmwtBLF9YAhHGoXE3GCoA/emVObB0Jk5bP2bYicKTpGSk49N8QGSsAfnE8c7J5dKc2PPei9h04qMFJA1qyYK4WDs3s/kM+Z9tcs3s43uvoUNMzUnDovyECKwAwpDXAtD4+LunIAFNJ0YIr1y/kcYOWvO0elkY+IKk5QFZjWw99rfB1dmt6dum5oTdH/huiFAgAhmx47sXDQVVLfXy8OcDUdz4MYN1+pQ4HVYcdPHR7CIeanlEMZKwAFE6icpwBraN2ot5uTZLMkw+ZnfFd8W6XnGp6RjGQsQJQKD3PE8zQ4KSBWLdbkTTzVNTMjkNNzygGAisAhZK4HGdgZ9iSBXNVq078a7VWrWrJgrmxHyszSQcp+jCAtTIj3u0SM7pgHKVAALHYKrNFlagcZ2hnWOt1cPn1SZx56jZDzBVDt0gjfy/pjbYb+5u3d+NI0zOKgcAKQGSpdr3lZHDSQGAQ1bUcZ7B/aOHQTGdei6PsHlazUFE/+nu9Mk8ObWcP5cMaUXiUAgFE5sOut0TluLz6h2wOomxl5YKCqqiZJx8GsPqwRhQaGSsAkfmw6y1ROS6PnWG2B1EGZeUkSTV6igCDCKwARJaozGZB7HJcHv1DtscVhGbfxgmqAIMoBQKIzItdb0nksTMsbbkxbRnRh119pnEGICwgYwUgMi92vSWV9c6wNOVGE2VEH3b1pdU++LQ6XRp/TYd3CBbtDEDXjxcqMQIrALE4vevNZWkCGxNlxKLvmOsMPoOmrbs2KT4p2/166IpSIAAEMV1GSlNuNLVr0acdc3Ff/9Dm/A6+T4qXkg95RS7IWAGwy8WSRlYZgaTlxrKdZ5fk9Y8aMBXhNSvq8UIFYSRj9corr+hDH/qQzjvvPF1wwQW64oor9PLLL5t4aABFlvTcuqy5lhEo23l2SV7/KAFTUV6zMm5E8IiRwKpSqeiDH/ygvve97+n+++/XH/3RH+nmm2828dAAisy1AKbFtYxA2c6zS/L6BwWfGjh0TmDBXrOyBdqeMRJYHXfccTrzzDMPf3366adr27ZtJh4aQJG5FsC0uJgRiNIfZXu8gKnfn+T1Dwo+59wuLdrpR09ZHGULtD1jvHl9fHxc3/zmN3XuueeafmgAvul1oXUlgOlc5+D5/mUEbJdVTf7+pBkZn5rz0yrTc/WM8cDq05/+tKZMmaLLLrvM9EMD8EmUC60LJY2gdb72b9Kx7/crI2C7rGry99vMyNjO+sF7RncFrlmzRlu3btXatWtVrTLJASi1KLOXXJitFLbOfQ81MwG+sF1WNf37sx7YGoT5UDDAWPTzhS98QZs2bdKXvvQlDQy4dW4YAAuiXGhdGLVgOyAxJY+yardsjitl3TRsZ/1QCEYCq2eeeUb/+q//qpdeekkrV67UihUr9NGPftTEQwPwVa8Lre2eoM71RL3dVVmXVXu9X3n8/qxLdEUJsmGVkVLgW97yFj311FMmHgpAUfQ6wsXEMS15rNMXWZdVe71fWf7+JCW6JNnQsg1iRSZohAKQjV4NyK5kB4q0dT3LnWJR3q+0vz8sKxW3RJc0G2p7MwWN84XAkTYAstOtAdml7ICNRmnfZP1+dctKxQ3Ck2ZDbW6moHG+MMhYAbDDdnYA8WT9fnULhuL2waXJhtqaD0XjfGEQWAGwo0glON8kKTll/X51C4ZmfUZSf8c3+sODOh83JLhSGkdqBFYA7GF6dP7S7Mbs9X6l6RHqGQxVOr7R+XWbXtk1F3uZfAwGEYjACgDKJKuSU9rxGd2CoR3XSDrY8QMHw9fcLbvmypiPTpTGC4PACgDKJKuSU9qArVswlGTNYdm1qOvMO6tFabww2BUIAD5JO63e5O6+9rWoEXyfOAFb2O5Mk2uOeiKAjR167E4tBDJWAOALE2UsUyWnzrWEMdEjZLJMFqWXiR16SIGMFQDjNo/u1IbnXtS+Awc1OGlASxbM1cKhmbaX5T8T0+pNzWoKWksnUz1CadfcnlmrTldzh+Eb4etkhx5SILBCJriwltfm0Z167Kmtqo+PS5L2HTiox55qlnH4DKQUcmFvvPE7/fv6jdH/OzNRcuoaZFTMD9dMuubOst74LkkDUmWG1Hg5eJ1hpUc1mv1WSYM6WweNI1eUAmFc68K670BzF8++Awf16G+26CdPbbG8MuRhw3MvHg6qWurj49rw3IuWVpSAi9vxpdAy1t767MMB7ObRnVbXor75bo3PCMysHZRqU8PXGVR6bIlTfnV1ByIyRWAF44IurJL022078/tLH9a0Auqot7ti8+hO3b1+o9Y9/mmNbfugmxfDwfPVOb/pjfHJ2rDnI5JyDmB9GQ+QdEfhse9X6KysqP1W9GqVEqVAGNftArrhuRcpBxngcql1cNJA4GdgcNKAhdVE016+XDLry+qr7J94h7h9TFnYPSy99m9qbxRvNKSnX3+3ntu/7PBtuQWwNs/ViyPpjsJ9D6lrU36Ufit6tUqJjBWM63YBdT1r4YOgUmuuJaAeliyYq1p14l8ttWpVSxbMjfwYrezRV9c9obvXb8z8ubVnWafWtgffyfbFMCD7UalIJ07+yYTbcg1gTU7Oz6r8mjSz1uv9jrLbMfQ+VTcyoMgEgRWM63YBdTlr4QvXe5gWDs3U0kXzD7/Xg5MGtHTR/MgZNRuBY3vAv7c+O/hOto8WCbnQtweCcQNYZ2TVi7R7WNp+5cSAtDoj2uDNbu931JJnaK9W3Z3yMoyjFAjjFg7N1Pbde/TbbRMvhCb+0ne5BJYXH3qYFg7NDH9feuyS6hY4ZvVet5cvN+z5iJZO+6z6q23lQFu9Q6O/lx68S3rsu9If/lQaGJf+5y7prO3S9OZ6Xx8fkiS//3swMUai0+5haeQDOuoonPHXov38rM9M3E3YUp0hzb4l2rpa9xl5v6T6xO+5UF5GJgiskIm3LzpZs6cdkzoIag+kBmo1jY2Pa7zR7Hso6zZ+H3uYDosw0dpG4LhkwdzDPVatfqUzjv2yBmvbVbHVO/SLn0j/5+PS2JhUH5NUkQ7WpP+aKW2cIf3tc9IpYxqc9zmt/h9n5Ls207LoRQo8X1CS3ogW0JjqIZu2Shp5X/D3bJeXkQkCK2Sma9Yigs55SAfr9aPuk3Umw0XtQUCLNyWgCJkJG4Fj6/PTCuK3Ny7U6Ix/NPu5ijPPaPT3zaDqwP6jvzdelcYlfWuhdN0HipHxMHlkTUu3oCVqQGPqiJksnh+cRWAFZ4WNbejkUgksD51BgFcloAiZCVuBY9p/CHQV9+y5B+9qZqq6Ge+XfrRXWmh2qVYEld3Sll9Dh3wq/4Ami+cHZxFYwVlRAyYvSmCGZRoEZCnCv9y9DhzDxO0heuy7h8p/XdTHpB8/JH3g4+bWaUsWoxtmfSa4x0r9+Qc0voymgBEEVnBWWEmonTclMDRF/Je7t4FjmLg9RPt7nMHX7X6+HqFiquzW/nhSc1fg+K7m/6/MkIYiNp6bZvr5dePrZ6AgGLcAZwXNQ6pImtRXkxR/Gz8cMG1Vc6t733w1z5ObH23re95Mz1QKPf4l5PbJIcep9LpfnLEFrh7bY9K0VdIpO6W3Npr/W7TTvc+aaRyjYx0ZKzirkCUh5Psv9yTi9kNFEbfHZunfSD/4TvdyYK1P+svzJ94WteQY9BxH3ieNXNYMdslw+CuL0RWIhcAKTitcSQjuy+LCFLfH5t2XST96oHtg1dcnnd/x81FLjoEHEx86vsVEIAl7OEbHOkqBADKR97E0xqS5MHUrr8U5/mVonvSxNdKkyc3MVLtaX/P2j61p3q9d1JJjr+fiwkHBZShVZiFu2RnGEVgBMM718wy7SnphMt3b8qdvl9bcLb3rIulNg82DAd802Px6zd3N73eKei5elIuszQxHktcyTSBWpCAu6dmIMIZSIFBwNo4BsnEsjTFJZw5lUUIcmtccpxB1pEKUkuPuYam+t/dj2cxwxH0t0/TFZdFTZxOjHawjsAIKrHN6fV7HAPlwnmGopBcmV3pbum0O6AwiwtjOcMR9LdMEtUVs9nZ9g0jBUQoECqxb5ihLYUNbvRnmGqcfqiWL3hbTJarApnU15zuZHoGRZu1xX8s0Qa0rATEKg4wV0ENepbQsfo+tzJHX5xkmZfrYkixKVGHBQuNlaaHB/re0a4/7WqY5i49z/GAYGSugi7yasLP6PbYyRwuHZmrpovmHf08phrmaHn7arUSVVF47xtKuPe5rGadhuzOTNng+zd4wiowV0KYzazRWr+fShJ1Vs7fNzFEpZ5CZ7G3JokSV12HAJtYe57WM2hcXlEnb/X+lyqQj97F57A0KgcAKOCSo0TuM6VJat5JdmhIh0+s9FqVE1Xkm3OD50r6HwoOLPHaM7R5WsxhS775206IEYoE9Zm9IjTfavv6D6ZWhZAisgEOCskZhTJfSwg6cntRXS72rr5SZoyLolV0KzL585ch9w/qastwx1lpTUFDlQnktSsbM9x2BsI4eK+CQqFmoLEppQQdO16pVNRqysqsPMWUxYLJXn1HYDr923fqa4qw56n1D11SL328W9zWNcv+oGTN2BCIFAivgkLAs1ECtlnkTdliz98F6wL/85ck8KFNcn4pteuJ6u25jH6Je/IPuF2fNce4buqbx+EFVnNc06v2DmtyDsCMQKVAKBA4Ja/T+i1NOzKWUFlSya/VGdfJmHlRaPkzFNjlgsrNnqlv/U1gPVtD90qw5zn1NjS6I+5pGvX9nj1llutTYI6ntvzEXSpbwGhkrlFbnIcGSnBsREFYiLPQ8qHZZjBwwzdTuvbhZmijZl7AgIc6a49zX1Dl1cV/TOLe3ZwEX7ZTm3G5+OGoSrmdmERmBFUopbG6UJK086zStPucMrTzrNOtN34WeBxXlQuLDVGxTs6HiBpFBPVjT/jFakBBnzXHua2qWV9zXNM17kGTKvmkmy8kEaNZRCkQp+XRIcCF39UUt8fkwFdvUbKgkQWRnaWvfQ9HGJ8RZc9znZ2LXYdzfmdd8rqyYKif7UDovATJWKCWvDwm2rLOEmmg6fI/sTOt3rNuxWmONyRPv59oF01aWRkqe6YizZtMT5aOI+zttrNEkU5lZH0rnJUDGCqUUNjeqNE3hCQUNUY07V0tS1wtJ++94TsukV6Uzjv2yBmvbVenR0J3XuY5HsZGlkdJlOuJONs86SAlq3F/4fPSfz2ONWTGVmfWhdF4CZKxQSl42hTvQO9GthBpLl+xM5+94bv8y/ftL9+mOkZ9p89QNXYOqPM51zEySrEveF9KsPoNBmbeRy6SnZ5ajR8hU039eZ0GiKwIrlJJ3TeFZzkqKwVgJtcuFJOyxGlLXQMlY0GdT3EbqLC+knUHUyEey+wyGDRYd32Xlc547U6VMUwEaUqEUiNLyqincUHNr2lKZsRJqlzPrBidtDA2uum0wKGXfXFZN24HH5axVM7xtY+r4l24ZtrIcMWOilJnHWZDoicAKMCizHh8DJR8T/VFhQ1QTlVBDLiRBv6NdWKDkZd9cnIGgQbK6kAZmkBqBdzVSduw17JQeoeh87jUrCEqBgCGZ9vikLPlsHt2pH/1mS+pSWR4l1NbvqIR8PyxQ8qpvbvew9NTMZh9R2tJaUPkwbS9UnEDGRNmx17BTeoTgETJWgCGZzsZKUfJpBXwh+YbYpbI8Sqitx4+THWv9jJVdgYdEylh2ltnamSh7mZhlFJpBqmhC5spU/05rXaNXSo1dHb+SHiH4hcAKMCTTHp8UJZ+ggK+dq6WyJIGSzb65yKXWsEbtlrRlLxP9eGGB/LHvbw4hzaJ/p1XCSlseBSwjsAIMybzHJ2HvRLfAztlS2SE+bTCInLHsdXBy2rKXiREMNpug6RGC5+ixAgxxtccnLLCrSG6PmHBBjF6lSBnL3cNSaPeYzJS9wgKzyvR4fVcunKEHeIjACjDE1dlYYQHf2YtPjr02I8fZ+GL3sDTy9x1DK/8+NCAJC2An3L7jGoXurqvOMHMMy+D5Ojp465cae6zPQTPGgWG5QBhKgYBBLpauTDV1GzvOxhejV0p6o+PGN5q3RxwTcVTGsls57hQDQeruYem1f9PE4K0iaUDSvon39XU+FAcNw3EEVkAJdAZ8rcxTnEAr012PLurcndbj9kgBbOiZcPPTrrZp+5Uh86f2Bd3bz/lQhoblAlkhsAJKJmnmqZSTzWPqmbHMalK61MzkjIcEg2F8nA8V2py/tVkWZDchLCOwAkomaeYp6q7HzKbP5606IzhQqc5I/phZ7rbbcU3496ozpMYfsgno8tZtxlbrdsqDsIjACoihCEFD0sxTlB6iNH1Yzr22s2+RRj4gqf11GWjenkZW4wS6lfVaay7CfKigrF/n4FKp+f3tVxbjOcMrBFZAREVp3k46bytKD1HSbJiTr23K7FLugWJoJmdw4nOY83W/g4ug9yVsNtj4riNZR7JY4RjKahSBFRBRUZq30xyk3KuHKGk2zNnXNmF2yUqgOHi+tHutJmZu+iW9UbwSWef7svmk3oNXJZrcg7DL0jjmWAERFaV5O8t5W5FmOQUoymvb0i1QNKJzjtPIR7qMWeh4DVvBhUm250r1OsS5nY87IbPUbZclEiFjBUSU+ZE1Ocpq3lbSbFiRXlsp40AxKMNwVKZKym3MggsZj87yYHV6+A5JH3dCZsnEEUiYgIwVEJGrR9a4JGk2rGivbdLMXSSBhziHTHMPYzK4cCXj0X4ET3VqyJ0qfu6EzFLYZ4EANDEyVkBEpiaYF12SbFjRXts0fWw9pc0kmB6z4OJcqdDXqEHfUKcsZ6uVlLHAasuWLbr66qv16quv6rjjjtOaNWt00kknmXp4wAkuHllTFEV6bTMNFLvNceqZuaqZOY8w6npsNc1nPeG+SLKcrVZSxgKr6667TpdeeqlWrFihe++9V9dee62+9rWvmXp4AAXl3PwqQzILFMMyDMe+X9r30KFsTViANW7+ghlnrlReO/LIwsST1Wy1kjLSY7Vr1y49+eSTWr58uSRp+fLlevLJJ/Xyyy+beHgABdUaS9Bq6m6NJdg8auBA4qKatqqZdeqbL6nS/HPoNmnOl4/0GIVlZrLomwlaT1hgl1VDdOeuRCn4NSJ4QA6MBFYjIyOaPXu2arWaJKlWq+mEE07QyMiIiYcHUFCZjyUoqvZG7YXPHx0wBI0fyDJj07mesMCuOt38WIbWrsSxrZIaE8uO3V4jICPsCgRgjZX5VbZnLuUhLKuVV3AROFdqQBp/7egAKO3r78quROAQI4HVnDlztH37dtXrdUlSvV7XSy+9pDlz5ph4eAAFlelYgiBh2Y2iBlcmMjZJAtGgwK5yjKQ3Jt7PRADEHCY4xkjz+owZM7R48WI98MADWrFihR544AEtXrxY06dPN/HwAAoq07EEQbplNwxkc5I04jvdvJ9m+GdnQ/RvQ/4dnzYACt0ByBwm2GGsFHj99dfrrrvu0nnnnae77rpLN9xwg6mHBlBQWR6vEyjD7EaSRnznm/dNltmyGkSZdz8Z0IOxcQt//Md/rP/4j/8w9XAASiLX+VUZZjeSHCTt7OHTLSYD0axGIDCHCY5h8jqA8ohxcY9Tots8ujNRI77zh0+HDv9sNPut4gQwWQZAzGGCQ9gVCKA8Iu6Wi1Oia903TLdG/Nyb9+MK3N13SJLG/zgN9WXYvYlCIrACUC4RLu5x5msF3belVyO+84dPTwhEA2Q11qBMuzdROARWANAhTomuW9muVyN+7s37SbQCUVWCvx/Wb5Um48RsKniMHisA6DA4aSAwYAoq0XW7b5QAKXHz/u7hfBu24zT+pxnTILkxmyrv1xeFQcYKgB0O99DEKdFZKefZKJXFGWuQNuOU1WiGqChFIgUCKwD5y+PClSJwi1Oi67zvpL6aapWKHv3NFt29fmM2M6nyKJWlOdg4bcbJ9mwqSpFIgVIgvOD0dOoefF57ZjKegB5Wihp9dY9++OJfRHov4pToWvdt7RBsNbO3dhO27mNM1qWysFLe0G2H+q16SDsvzPZsKhdKkfAWGSs4z/np1F34vPZMZX3hCgncpu69IdP3Is5uwlSyLpWlzdiYyDiZOuswCdulSHiNjBWc5/x06i6yXru32bAEGY1YzzUkQBusbp/wtenPUW4DP7OaYt6SNvC1nXFKK+vXF8pYZRcAAB4gSURBVIVGYAXn2ZpObSJoyXLtuZWdshDzwhX7uYYEbnvrs4+6zeTnKM5uwlSyDlxMHP3j8zR03wNDWEUpEM6zMZ3aVAkvy7XnVnbKQsQJ6C2xn2tAKWqsMVkb9nzkqLua/BzlukMwy1KZ7eZxF9gsRcJrZKzgvCUL5k7IVkjZb2c3VcIzvfb2LFoYZ86Z6yVGRiN25i8g47Cz/2ptfekMSdl9jlqfDS/Ls+3I2ACJEVjBeTYuVqZKeCbX3lkOC5PHOXN593YlKrF1BG5DkpZOyn7diQd+BrE5pNLnUh5gEYEVvGD0YhWByV4ZU2vvdiZdSx7nzNno7TKV+cv7c5RK2unlAKygxwoI4OLhuL2yZXmdM2ejt8uLM/XCJB1UypBKwEtkrIAALvbKdMuirTzrtNzWkdcuzaByY57P04g0WaegXXkSQyoBxxFYAQrvGXIpI2KjiT9IHiMFvB4l0a5X1imsf2r3sKSKpMbRj8mQyqNxYDIcQmCF0vPlIu5KFi2PAM/2UFhjzfmhgza3ds9k7bhGgUGVKuUaeRAFvWhwDIEVSs/2RTwOF7JovQI81wer9mI00A4btKla97MSQ8t9DYKFTlmfOwnERGCF0rN5EfdVWIBnKigZqNV0sF4PvD2xiOUio4F22IT5zkCgpRVQhU4+nx/v95cBBybDMewKROnZmOxeVKZ2DFYq3W/fPLpTd6/fqK+ue0J3r9/YeyJ+q1w0tlVS40i5KGCHntFAO2zCfFiA1OqfYvJ5dByYDMcQWKH0XByt4CtTQcmBsaOzVa3bEx03FGN0gfFAO+holF6BU8wjf0qNIBSOKXwpMO8J0fCPK03hRWBqx2C3x0lUqotRLspl92WUI2OYfB6Nz8fvsJuxkAodWPmy2wv2udAUXgSmgpJuj/Pob7YE/kzXrFhoz9LR5aL2QHt25T792bFf0ZTaqCp7T5R2G7zwETiZ4+NryW7Gwip0YOXTbi8UQ9kzpKayf90eJ+wQ6q5ZsbAm8oByUes9nF25T3953GfVV9nf/EaGFz5vPzdkXJJjN2NhFTqwYrcX8uRrhtT0Rd1U9i/scRJlxSKWi9rfwyWzvnwkqGrJ4MLn6+eGjEtK7GYsrEIHVnlMiAZafMyQ+nhRT5wVi1Auan8Pp9a2B98pwYWvW/DqwucmUXBNxiWdGOVp+KXQgZUrR4CgHHzMkLpwUU8iq5649vdqb322jukbPfpOMS98vYJX25+bzvXNrtynoV1fVuPV7ap0K++RcUknRnkafin0uIWFQzO1dNH8wxmqwUkDWrpovtMXDGQn9uyjmHych2X7ou6a9vdqw56P6I3xyRO+P643xb7w9ZrtZftz076+BZMf1tJpn9XU2qgqPeZ9MT8qJUZqFFahM1YSu73KKKisISnzkpePGVLK5RO1v4fP7V/WvO2YL2tqbbv21mfrv/Z9VHOPO08Lp0V/zF7Bq+3PTfv6lhzzZfVXI/aVkXFJz8fdjOip8IEVyiWs7NJXrWRe8vJxHpbti7prOt/D5/YvOxxgtYzE/Mz0Cl5tf27a1xerr8zn+VFAhgisUChhZZf6ePD9TZe8fMuQ2r6ou6j1Hn513ROB34/7mYkSvB7+3LTGF7z6O2lvPoFK+/pi95WRcQGOQmCFQol70Stryaudb8FgXhKVSQPmOi0cagYePYNXS+ML2oPrDXs+MnF2l0R5D4iJwAqFEnYxHKjVVG80KHkhsthl0i6B0cKhVU6PLzgSXJ8m7V5QqPKet8NX4S0CKxRK2MXwL05pljKK8hds1hcLLkYJyqRpAyNXxhcUqLzn45w2+I/AykNc9ML1uhi68jqleQ+zvlhwMToiVpk0bWDEwEjjfJ3TBr8RWHmGi15vUS+GtgLUtO9hmotFlOfMxSihtIER4wuMY04bbCj0gNAi6jVsENG0gpvWX7Ct4Mb00NAgad/DpBeLqM+Zi1FCsz7TDITaxQmMGBhpnO3hqygnAivPcNEzw2aAmvY9THqxiPqcuRglZCIwmrZKWvi89Nbx5p8EVaksWTBXterEyxybVpA1AivPcNEzw2aAmvY9THqxiPqcuRilQGDkFI41gw30WHmGSdlm2DzKJe17mHSoZ9TnzNBQFAlz2pA3AivP+HTRc3n3os0A1cR7mORiEec5czECgGQIrDxk6qKXZeDj+u5F2wGqjcDF9nMGgDIgsCqprAMfH7bslzErU8bnbILL2VcAbiGwKqmsAx92L7qDoCAd17OvANxCYFVSWQc+NpvDXZE0oDEZCCUNCgjGjvAh+wrPBRzezY5SfzFuoaSyHttQ9i37SQeQmh5cmmRel83hqS4qRPZ197C0+STpt9Xmn7uHba8ILa3Du8e2SmocObyb98hbZKxKKutdcWVvlE6a5TCdHUkSFLiQoXEpY+Z99rV14W4dldO6cEvStFVOvdallPbwbjiHwKqk8gh8ytwoHTWg6byomc6OJAkKbGdoXOtp8n52XJcL9+Y/nOfUa11KaQ/vhnMIrEqszIFP1qIENEEBRLfHSyJJUGA7Q+NCxqyd99nXLhdu117rUkp7eDecQ2CFQrJd3ogS0ARd1IKEBUJRnmO3oCDs521naGxnzIJ4/Y+QLhduF1/r0pn1mYmlWine4d1wDoEVCsdkKalb8NLte1GyHL0yVN0CpjjPMSgoiPLztgJTUxkz28G1M7pcuAdf8rx/rAhafVTsCiwMAisUjqnyRrfgQ1LPwKRXlqNbALHyrNO6ri3tc+z18zYzNCYyZkHv3Y9/+7wef/p3OlivlyvQ6nLhXrJgp9/9Y0UxbRWBVIEQWKFwTJU3eo0qSBu8pQkguj3Hu9dv7Bk0uFwCMpExC3rvxhsNHazXJZWwSTvkwm07OwkUEYEVCsdUKSlJ8BEnMEnS/9TSawdhr6DBdoN6L2kzZlHeB5q0m7zuHwMcRGCFwjHVfN0r+DARmCTtfwp6ju16BQ22G9Sz1i3wbOdChg5AsTB5HYWzcGimli6afzjIGZw0oKWL5sf+V3m36fFZTpaPMi298zkG6RY0mHqNXBX0/gRxJUMHoDjIWKGQTJQ3ovSfZNGbErUE2XqOd6/fmCh7VuQSUOd7N6mvpoNjdTXa7lOkDB0AdxBYAV10Cz6yCkzi9j8VvayXVOf7w/gFAHkgsAJSMHWxbn+cSX01VaTI2RXfdnbZCnCKnKED4A4CKyAhU4NIOx/nwFhd1UpF/dVq5JlLvgQNWZ4DSEYKgAsIrICETA0iDZu51N9X0/vO/l9G1uqKrM6mc+3gZgDllTqwuuGGG/T4449rYGBAU6ZM0TXXXKO3ve1tJtYGOM3UkE2Xh3Um8Xsd0F16Sd/Vy3pd45qiqv5G03WZToj1XONkoHw7TJjsGlBcqcctnH322br//vt133336R/+4R901VVXmVgX4LywZvK4W/hNPY4LfqLdWqnf6jvaqX0aV0PSPo3r/43v0P+uP6lts2qBP9f5XFsZqFbA1cpAbR7dGfjzPgWncZ9bL6Nb12rfr9+sxm+q2vfrN2t061qTywUQU+rA6pxzzlF/f78k6fTTT9fo6KjGQ4YWAkViapZVljOx8vR7HdDH9bz2a1xjHd9rVCsaq0kPLpb2vKky4XtBzzXKLK92PgWncZ9bN6Nb12rmvqs0WBtRpdLQYG1EM/ddRXAFWGR0QOjw8LDe+c53qhphMB/gO1NDNosyrPMuvaQxdf9HVb0i/eqPqj2fa9wMlE/Bqcns2jF7P6W+6v4Jt/VV9+uYvZ9KtDYA6fXssbrooou0bdu2wO+tX79etVoztf/ggw/q/vvv1/DwsNkVAg4ztRvPl1193XxXLx+VqerUqFb0m9kN3T73tK73izvLy6eREybPaZxSHY11O4Ds9Qys7rnnnp4P8sgjj+iLX/yi7rzzTs2c6d5fZACy93qPbFXLG7VKz/skGXqaZXBqstnc5EDX18eHNFgbCb490eoApJV6V+C6det044036o477tC8efNMrAnIBDuxsjVFVe2LEFy9qVGResRWLmWgTI9yMPnc9ky9VpP2XTWhHDg2Pll7pl5LYAVYkjqw+pd/+Rf19/frn/7pnw7fduedd+r4449P+9CAMXnMOSp74PY3mq7vaGfXcmCtIS2vzoj0eK6UR7MY5WDquQ3N/7BGtzZ7raZUR/X6+JD2TL1WQ/M/nPqxASSTOrD66U9/amIdQKaynnPEgErpMp2gB/Ry1wb2/kpVq3RCjqtKz/VRDs0gqhlIDR76XxJl/4cBYArb91AKWV8cTW6h99U8TdIanaTJqh71L7Y+SZNV1RqdpHmaZGN5ifk0yiEp07O1gDLjSBuUgsmdWEFMB26+Zg/ermm6W2/VsF7SQ22T18/XdK3SCd4FVZLZZnNX+Ta5HnAZgRVKIeuLo8nAzdeyYnswOH/SgL7qSTDYi0uN9FlxvdwJ+ITACqWQ9cXRZODmY/bA12AwKlca6bOSdUYXKBMCK5RGlhdHk4Fb1tmDLMqMPgaDLsu7FFyGcieQFwIrwBBTgVuW2YOsMkuUksyxkf0rQ7kTyAuBFeCYLLMHWWWWKCWZYyv7V/RyJ5AXxi0AjsnyUOasMks+HYLsOrJ/gN/IWAEOMpk9aO/XqUhqBNynM7MUt8eHUpI5ZP8AvxFYARbk1Zzc2a8TFFR1ZpaS9vhQSjKDRnLAb5QCgZzlOeU6qF9HOnIGclCZkSnydmVZCgaQPTJWQM7CApefPvM741mssL6chqTV55wR62fK2ONjawI+2T/AX2SsgJyFBSgHxurGs1hJzrkrw9l4UXB+HoAkCKyAnEUNUEyU35Ls1mOHXxMlUQBJUAq0zNfDdosqj/cjqDk5TNryW5Ldeuzwa6IkCiAJAiuLin6+mm/yej+CApc3xuo6WK8fdV8T5bck/Tr0+DD2AEAyBFYWcb6aW/J8PzoDl86gTipn+c0ljD0AkASBlUWUGtxi8/2wXX6jJH002+8JAD8RWFlEqcEttt+PbuW3LAOfoBLoj36zRT995nc6MFYvdUBBSRRAXARWFlFqcIur70fWvV9BJdCGmuMfsvh9riJrB8AExi1YxIRlt7j6fmS97T9KqbPoYwaYWQXAFDJWllFqcIuL70fWvV9hJdCsfp+L2EgCwBQyVoDjsp6EHjQQNMvf5yI2kgAwhYwV4Lise786d78N1GoaGx/XeKNh7Pe53r9ke+MCgOIgsAIc0i0AyTIwCZqrZer3+TAI19WNCwD8Q2AFxJD32IP2AKTz92S5FpO9Zj70LzGzCoApBFZARDbGHoQFID5kgVp86V9yceMCAP8QWMH5/hdXZJ15iROAuJ4Fav9MVdSci9WJ/iUARURgVXI+ZT5sszX2ICgAcTkL1PmZCgqq6F8CUFQEViXneubDJVnvHIvTQO3yLragz5Skw5mrvLKiZc7Elvm5A7YRWJWcy5kP1+Q99qDbBdHlXWxhn52GpNXnnJHLGsqciS3zcwdcQGBVci5nPlxjY+yBzbUkFfaZqkj66roncllrmTOxZX7ugAsIrErO5cyHi1zaOebSWtoFfaakI71WSTMoccpbZc7Elvm5Ay7gSJuSc/XgYfir8zNVCbhP3EOd4x6SnPUxQC4r83MHXEDGCs5mPuCv9s/UV9c9EXifOBmUuOWtMmdiy/zcARcQWAHIlIk+vrjlLZd70LJW5ucOuIDACkCmTGRQkgRnZc7Elvm5A7YRWAElYHOukYkMCuUtAL4gsAIKzoW5RmkzKJS3APiCwAowxNVp10WZa0R5C4APCKwAA1zICoVhrhEA5Ic5VoAB3bJCtjHXCADyQ2AFGOByVmjJgrmqVSf+p07jNwBkg1IgYICLZy6293wN1Grq66vpwFjdqf4vACgaAivAANfGAXT2fB2s11WrVvWOxScTUAFAhigFAga4duaiyz1fAFBkZKwAQ1waB+ByzxcAFBmBFVBALvZ82eTqjDEAxUMpECggdgIe0eo3awWarRljm0d3Wl4ZgCIisAIKyLWeL5voNwOQJ0qBQEG51PNlE/1mAPJExgpAoTF5HkCeCKwAFBr9ZgDyRCkQQKG1yqHsCgSQBwIroKTKNIKAfjMAeSGwAkqo88ib1ggCSQQgAJACPVZACTGCAACyQWAFlBAjCAAgG5QCgQy43r9k8sgb158rAOSJwAowzIf+pSUL5k5Yo5RsBIHJ50qABqAIKAUChvnQv7RwaKbeMjRdlUNfVyS9ZWh67EDG1HPlPD8ARUHGCogoakbFh/6lzaM79czoy2oc+roh6ZnRlzV72jGxgitTz7VbgEbWCoBPyFgBEcTJqPhwhIqpTJOp5+pDMAoAURBYARHECUR8OELFVCBj6rn6EIwCQBTGAquf/exnWrx4se666y5TDwk4I04gsnBoppYumn84KBicNKCli+Y7VdIyFciYeq4+BKMAEIWRHqu9e/fq5ptv1tlnn23i4QDnxB1P4PoRKqZ2BUpmnivn+QEoCiOB1U033aTVq1frhz/8oYmHA5xjMhBxgSuBDCMWABRN6sDq0Ucf1Z49e7Rs2TICKxRO+4V/Ul9NtUpNB+v1QgQBtrNqPsz7AoC4egZWF110kbZt2xb4vYcfflif//zndccddxhfGGBb54X/wFhdtWpV71h8Mhd+AxixAKCIegZW99xzT+j3NmzYoB07duiSSy6RJL3yyitat26dXn31VV1xxRXmVglYwIU/W4xYAFBEqUqBS5Ys0eOPP37466uvvlqnnnqqLrvsstQLA2zjwp8tk+cVAoArmGMFhGC2UrYYsQCgiIweaXPTTTeZfDjAqqLtBHSNKzsTAcAkzgoEQnDhz57tnYkmMDICQDsCK6CLIlz481DW4IKREQA60WMFIJU4B1QXjanDrAEUBxkrAKmUeSxFWXaOljUjCSRBxgpAKmUJLoKUYedomTOSQBIEVgBSKUNwEaYMIyModwLxEFgBSKUMwUWYhUMztXTR/MNB5OCkAS1dNL9QZbIyZySBJOixApBK2cdSFH3nKBPygXgIrACkVvTgoswYlAvEQ2AFAAhV9owkEBeBFQCgKzKSQHQ0rwMAABhCYAUAAGAIgRUAAIAh9FjBKo7KAAAUCYEVrGkdldHaxt06KkMSwRUAwEuUAmENR2UAAIqGjBWs4aiMYEUpjxbleQBAHARWsIajMo5WlPJoUZ4HAMRFKRDWlPnw3jBFKY8W5XkAQFxkrGANR2UcrSjl0aI8DwCIi8AKVnFUxkQmyqMu9DZR5gVQVgRWgEOWLJg7oTdJilcedaW3Ke3zgFkuBNtAWRBYAQ5JWx7t1tuU54WUMq87XAm2gbIgsAIck6Y86lJvE2VeN7gSbANlwa5AoEDCepjobSovl4JtoAwIrIACYYQFOhFsA/kisAIKZOHQTC1dNP/wRXNw0oCWLppPyafECLaBfNFjBRQMvU1ox0YCIF8EVgBQcATbQH4oBQIAABhCYAUAAGAIpUAApcQ0cgBZILACUDpMIweQFUqBAEqn2zRyAEiDwApA6TCNHEBWCKwAlA7TyAFkhcAKQOkwjRxAVmheB1A6TCMHkBUCKwClxDRyAFmgFAgAAGAIgRUAAIAhBFYAAACG0GMFOIojVwDAPwRWgIM4cgUA/EQpEHAQR64AgJ8IrAAHceQKAPiJwApwEEeuAICfCKwAB3HkCgD4ieZ1wEEcuQIAfiKwAhzFkSsA4B9KgQAAAIYQWAEAABhCYAUAAGAIgRUAAIAhBFYAAACGEFgBAAAYQmAFAABgCHOsUDqbR3cyeBMAkAkCK5TK5tGdeuypraqPj0tqHmr82FNbJYngCgCQGqVAlMqG5148HFS11MfHteG5Fy2tCABQJARWKJV9Bw7Guh0AgDgIrFAqg5MGYt0OAEAcRnqsvv71r2t4eFj9/f2qVqu69957TTwsYNySBXMn9FhJUq1a1ZIFcy2uCgBQFKkDq+9///t6+OGH9e1vf1tTp07Vzp07TawLyESrQZ1dgQCALKQOrG6//XZdeeWVmjp1qiRp5kwuUHDbwqGZBFIAgEyk7rF69tlntXHjRq1cuVIXX3yxvvWtb5lYFwAAgHd6Zqwuuugibdu2LfB769evV71e18jIiL7xjW/olVde0Xvf+16dfPLJOuOMM4wvFgAAwGU9A6t77rmn6/ff/OY3a/ny5apWq5oxY4bOOuss/epXvyKwAgAApZO6FLh8+XL9+Mc/liS9/vrr+vnPf663vvWtqRcGAADgm9SB1d/93d9pZGRE7373u3XJJZfoggsu0Nvf/nYTawMAAPBK6l2BkydP1uc+9zkTawEAAPAak9cBAAAMIbACAAAwhMAKAADAEAIrAAAAQwisAAAADCGwAgAAMITACgAAwJDUc6zSqtfrkqTR0VHLKwEAAOiuFa+04pdO1gOrHTt2SJJWrVpleSUAAADR7NixQ/Pnzz/q9kqj0WhYWM9h+/fv16ZNmzRr1izVajWbSwEAAOiqXq9rx44dOvXUUzV58uSjvm89sAIAACgKmtcBAAAMIbACAAAwhMAKAADAEAIrAAAAQwisAAAADCGwAgAAMITACgAAwBACKwAAAEMIrBzz9a9/XcuWLdMFF1ygFStW2F6OV372s59p8eLFuuuuu2wvxQs33HCDli1bpgsvvFArV67Uf//3f9tekrO2bNmi97znPTrvvPP0nve8R88//7ztJTnvlVde0Yc+9CGdd955uuCCC3TFFVfo5Zdftr0sr9x6661atGiRnn76adtL8cKBAwd03XXX6a//+q91wQUX6JOf/KSVdVg/KxBHfP/739fDDz+sb3/725o6dap27txpe0ne2Lt3r26++WadffbZtpfijbPPPluf+MQn1N/fr3Xr1umqq67Sf/7nf9pelpOuu+46XXrppVqxYoXuvfdeXXvttfra175me1lOq1Qq+uAHP6gzzzxTkrRmzRrdfPPN+uxnP2t5ZX749a9/rV/+8peaO3eu7aV443Of+5wmTZqk733ve6pUKtauoWSsHHL77bfriiuu0NSpUyVJM2fOtLwif9x0001avXq1jj/+eNtL8cY555yj/v5+SdLpp5+u0dFRjY+PW16Ve3bt2qUnn3xSy5cvlyQtX75cTz75JNmXHo477rjDQZXU/Ixt27bN4or8cfDgQX3qU5/S9ddfb3sp3ti3b5++853v6Morr1SlUpFk7xpKYOWQZ599Vhs3btTKlSt18cUX61vf+pbtJXnh0Ucf1Z49e7Rs2TLbS/HW8PCw3vnOd6pa5a+ETiMjI5o9e/bhQ+JrtZpOOOEEjYyMWF6ZP8bHx/XNb35T5557ru2leOGWW27RhRdeqHnz5tleijdeeOEFHXfccbr11lt18cUX633ve582bNhgZS2UAnN00UUXhf6Lbf369arX6xoZGdE3vvENvfLKK3rve9+rk08+WWeccUbOK3VLt9ft4Ycf1uc//3ndcccdOa/Kfb0+b61A4cEHH9T999+v4eHhPJeHEvn0pz+tKVOm6LLLLrO9FOf94he/0KZNm/TP//zPtpfilXq9rhdeeEF/8id/oo9//OPauHGjPvzhD+uRRx45XAXKC4FVju65556u33/zm9+s5cuXq1qtasaMGTrrrLP0q1/9qvSBVbfXbcOGDdqxY4cuueQSSc2G2XXr1unVV1/VFVdckdcSndTr8yZJjzzyiL74xS/qzjvvpPQcYs6cOdq+fbvq9bpqtZrq9bpeeuklzZkzx/bSvLBmzRpt3bpVa9euJSMawRNPPKFnn31W73rXuyRJo6OjWr16tW688UYtXbrU8urcNWfOHPX19R0u2Z922mk6/vjjtWXLFr3tbW/LdS0EVg5Zvny5fvzjH+uMM87Q66+/rp///Of6q7/6K9vLctqSJUv0+OOPH/766quv1qmnnsq/jCNYt26dbrzxRt1xxx2UHLqYMWOGFi9erAceeEArVqzQAw88oMWLF2v69Om2l+a8L3zhC9q0aZNuu+02DQwM2F6OFy6//HJdfvnlh78+99xztXbtWp1yyikWV+W+6dOn68wzz9RPfvITLV26VFu2bNGuXbs0f/783NdSaTQajdx/KwLt379fn/zkJ/Xkk09KklasWDHhPzD0RmAV3Z//+Z+rv79/QoBw5513sgEgwLPPPqurr75ar732mo499litWbNGCxYssL0spz3zzDNavny5TjrpJE2ePFmSNG/ePH3pS1+yvDK/EFhF98ILL+gTn/iEXn31VfX19eljH/uY3vGOd+S+DgIrAAAAQyh4AwAAGEJgBQAAYAiBFQAAgCEEVgAAAIYQWAEAABhCYAUAAGAIgRUAAIAh/x+3hHGL7RpYtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,8))\n",
    "mean = model.mu.numpy()\n",
    "\n",
    "plt.scatter(data[:n//2, 0], data[:n//2, 1], c = '#9cbeca' )\n",
    "plt.scatter(data[n//2:, 0], data[n//2:, 1], c = '#ffd700')\n",
    "plt.scatter(mean[0,0,0], mean[0,0,1], c='#FF5733', s=100)\n",
    "plt.scatter(mean[0,1,0], mean[0,1,1], c='#33FFCE', s=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea9bfbcf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}