saw_time_series.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Saw time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from torch.autograd import Variable  # CHECK: do I need this?\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "import time\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 425,
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.manual_seed(2012)\n",
    "\n",
    "N = 1000\n",
    "theta = torch.tensor([3.4159])\n",
    "a = 0.75\n",
    "z = torch.empty(N)\n",
    "x = torch.empty(N + 1)\n",
    "x[0] = 0\n",
    "for n in range(N):\n",
    "    z[n] = torch.randn(1) + x[n]\n",
    "    x[n + 1] = torch.randn(1) + a * (theta + x[n])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.plot(x, alpha = 0.5)\n",
    "plt.plot(z, alpha = 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distributions for ELBO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 427,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Vectorized for Monte Carlo estimation of the ELBO\n",
    "def log_joint(x, theta, z):\n",
    "  return -0.5 * torch.sum((z - x[:N])**2, dim = 1)\\\n",
    "    - 0.5 * torch.sum((x[1:] - a * (z + theta))**2, dim = 1)\n",
    "\n",
    "def log_q(theta, z, nu_mean_theta, nu_sd_theta, nu_mean_z, nu_sd_z):\n",
    "    log_q_theta = torch.sum(- torch.log(nu_sd_theta) - 1 / (2 * nu_sd_theta**2)\\\n",
    "                    * (theta - nu_mean_theta)**2, dim = 1)\n",
    "    log_q_z = torch.sum(- torch.log(nu_sd_z) - 1 / (2 * nu_sd_z**2)\\\n",
    "                          * (z - nu_mean_z)**2, dim = 1)\n",
    "    return log_q_theta + log_q_z\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SANITY CHECK\n",
    "if False:\n",
    "    N_mc = 10\n",
    "    theta = theta.to(device).repeat(N_mc).reshape(N_mc, 1)\n",
    "    z = z.to(device).repeat(N_mc).reshape(N_mc, N)\n",
    "    x = x.to(device)\n",
    "    nu_mean_theta = torch.randn(1).to(device)\n",
    "    nu_sd_theta = torch.exp(torch.randn(1)).to(device)\n",
    "    nu_mean_z = torch.randn(N).to(device)\n",
    "    nu_sd_z = torch.exp(torch.randn(N)).to(device)\n",
    "\n",
    "    print(log_joint(x, theta, z))\n",
    "    print(log_q(theta, z, nu_mean_theta, nu_sd_theta, nu_mean_z, nu_sd_z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model class for VI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 429,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Inference_NN(nn.Module):\n",
    "    def __init__(self, input_dim, hidden_dim):\n",
    "        super(Inference_NN, self).__init__()\n",
    "\n",
    "        self.FC_input = nn.Linear(input_dim, hidden_dim)\n",
    "        self.FC_input2 = nn.Linear(hidden_dim, hidden_dim)\n",
    "        self.FC_mean = nn.Linear(hidden_dim, 1)\n",
    "        self.FC_sd_log = nn.Linear(hidden_dim, 1)\n",
    "\n",
    "    def forward(self, x):\n",
    "        h_ = F.relu(self.FC_input(x))\n",
    "        h_ = F.relu(self.FC_input2(h_))\n",
    "        mean = self.FC_mean(h_).reshape(N)\n",
    "        sd_log = self.FC_sd_log(h_).reshape(N)\n",
    "\n",
    "        return mean, sd_log\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 430,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Model(nn.Module):\n",
    "    def __init__(self, input_dim = 1,\n",
    "                 use_avi = False, hidden_dim = 0, const_z = False, N_mc = 100):\n",
    "        super(Model, self).__init__()\n",
    "        self.use_avi = use_avi\n",
    "        self.const_z = const_z\n",
    "        self.input_dim = input_dim\n",
    "        self.hidden_dim = hidden_dim\n",
    "\n",
    "        # variational parameters for q(theta)\n",
    "        self.nu_mean_theta = torch.nn.Parameter(torch.randn(1))\n",
    "        self.nu_sd_theta_log = torch.nn.Parameter(torch.randn(1))\n",
    "\n",
    "        # variational parameters for q(z)\n",
    "        if (const_z and (not use_avi)):\n",
    "            self.nu_mean_z = torch.nn.Parameter(torch.randn(1))\n",
    "            self.nu_sd_z_log = torch.nn.Parameter(torch.randn(1))\n",
    "        elif (use_avi):\n",
    "            self.inference_nn = Inference_NN(input_dim, hidden_dim)\n",
    "        else:\n",
    "            self.nu_mean_z = torch.nn.Parameter(torch.randn(N))\n",
    "            self.nu_sd_z_log = torch.nn.Parameter(torch.randn(N))\n",
    "        \n",
    "        self.N_mc = N_mc\n",
    "\n",
    "    def reparam(self, nu_mean_z, nu_sd_z, nu_mean_theta, nu_sd_theta, N_mc):\n",
    "        epsilon = torch.randn((N_mc, N)).to(device)\n",
    "        z = nu_mean_z + nu_sd_z * epsilon\n",
    "        \n",
    "        epsilon_theta = torch.randn((N_mc, 1)).to(device)\n",
    "        theta = nu_mean_theta + nu_sd_theta * epsilon_theta\n",
    "\n",
    "        return z, theta\n",
    "    \n",
    "    def variational_z(self, x):\n",
    "        if (self.use_avi and self.input_dim == 1):\n",
    "            nu_mean_z, nu_sd_z_log = self.inference_nn(x[1:].reshape(N, 1))\n",
    "        elif (self.use_avi and self.input_dim == 2):\n",
    "            x_input = torch.transpose(torch.stack([x[:N], x[1:]]), 0, 1)\n",
    "            nu_mean_z, nu_sd_z_log = self.inference_nn(x_input)\n",
    "        elif (self.const_z):\n",
    "            nu_mean_z = self.nu_mean_z.repeat(N)\n",
    "            nu_sd_z_log = self.nu_sd_z_log.repeat(N)\n",
    "        else:\n",
    "            nu_mean_z = self.nu_mean_z\n",
    "            nu_sd_z_log = self.nu_sd_z_log\n",
    "        return nu_mean_z, torch.exp(nu_sd_z_log)\n",
    "    \n",
    "    def compute_elbo(self, x):\n",
    "        nu_mean_z, nu_sd_z = self.variational_z(x)\n",
    "        z, theta = self.reparam(nu_mean_z, nu_sd_z,\n",
    "                                self.nu_mean_theta, torch.exp(self.nu_sd_theta_log),\n",
    "                                self.N_mc)\n",
    "        Elbo = log_joint(x, theta, z)\\\n",
    "            - log_q(theta, z, self.nu_mean_theta, torch.exp(self.nu_sd_theta_log),\n",
    "                    nu_mean_z, nu_sd_z)\n",
    "        return torch.mean(Elbo)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run VI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 499,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_vi(seed, n_iter_optimizer, use_avi=False, input_dim=1, hidden_dim=0, const_z=False,\n",
    "           print_output=False):\n",
    "    torch.manual_seed(seed)\n",
    "\n",
    "    model = Model(input_dim, use_avi, hidden_dim, const_z).to(device)\n",
    "\n",
    "    optimizer = torch.optim.Adam(model.parameters(), lr = 1e-3)\n",
    "    loss_saved = torch.empty(n_iter_optimizer)\n",
    "\n",
    "    start_time = time.time()\n",
    "    for i in range(n_iter_optimizer):\n",
    "            loss = - model.compute_elbo(x.to(device))\n",
    "            optimizer.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "            loss_saved[i] = loss.data\n",
    "            if (i % 1000 == 0 and print_output):\n",
    "                print(\"loss\", loss_saved[i])\n",
    "    end_time = time.time()\n",
    "    run_time = end_time - start_time\n",
    "\n",
    "    return loss_saved, model, run_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('seaborn-white')  # seaborn-v0_8-white\n",
    "\n",
    "params = {'axes.labelsize': 18,'axes.titlesize':25, 'legend.fontsize': 15, 'xtick.labelsize': 15,\n",
    "          'ytick.labelsize': 15}\n",
    "plt.rcParams.update(params)\n",
    "\n",
    "linewidth = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# manual examine a few tuning parameters for A-VI.\n",
    "if (False):\n",
    "    seed = 2012\n",
    "    n_iter_optimizer = int(30e3)\n",
    "\n",
    "    n_method = 9\n",
    "    loss_all = torch.empty((n_method, n_iter_optimizer))\n",
    "\n",
    "    loss_all[0,:], model_fvi = run_vi(seed, n_iter_optimizer)\n",
    "    loss_all[1,:], model_const = run_vi(seed, n_iter_optimizer, const_z=True)\n",
    "    loss_all[2,:], model_avi_4 = run_vi(seed, n_iter_optimizer, use_avi=True, hidden_dim=4)\n",
    "    loss_all[3,:], _ = run_vi(seed, n_iter_optimizer, use_avi=True, hidden_dim=6)\n",
    "    loss_all[4,:], model_avi_10 = run_vi(seed, n_iter_optimizer, use_avi=True, hidden_dim=10)\n",
    "    loss_all[5,:], model_av2_1 = run_vi(seed, n_iter_optimizer, use_avi=True, input_dim=2, hidden_dim=1)\n",
    "    loss_all[6,:], model_av2_2 = run_vi(seed, n_iter_optimizer, use_avi=True, input_dim=2, hidden_dim=2)\n",
    "    loss_all[7,:], model_av2_3 = run_vi(seed, n_iter_optimizer, use_avi=True, input_dim=2, hidden_dim=3)\n",
    "    loss_all[8,:], model_av2_4 = run_vi(seed, n_iter_optimizer, use_avi=True, input_dim=2, hidden_dim=4)\n",
    "    loss_x2k6, _ = run_vi(2013, n_iter_optimizer, use_avi=True, input_dim=2, hidden_dim=6)\n",
    "    loss_x2k10, _ = run_vi(2013, n_iter_optimizer, use_avi=True, input_dim=2, hidden_dim=10)\n",
    "    loss_x1k20, _ = run_vi(2013, n_iter_optimizer, use_avi=True, input_dim=1, hidden_dim=20)\n",
    "    \n",
    "    algo_names = ['F-VI', 'const', 'n=1, k=4', 'n=1, k=6', 'n=1, k=10', 'n=2, k=4']\n",
    "\n",
    "    plt.plot(loss_all[0,:] + min_elbo, color = 'black', label = 'F-VI', linewidth = linewidth)\n",
    "    plt.plot(loss_all[1,:] + min_elbo, color = 'red', label = 'const', linewidth = linewidth)\n",
    "    plt.plot(loss_x1k20 + min_elbo, alpha = 0.5, label = \"n=1, k=20\", linewidth = linewidth)\n",
    "    plt.plot(loss_all[8,:] + min_elbo, alpha = 0.5, label = \"n=2, k=4\", linewidth = linewidth)\n",
    "\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.legend()\n",
    "    plt.yscale(\"log\")\n",
    "    plt.xscale(\"log\")\n",
    "    plt.ylabel(\"Elbo\")\n",
    "    plt.grid(which = 'major', visible = 'true', c = 'grey', alpha = 0.25)\n",
    "\n",
    "    plt.savefig(\"to_discuss/k_vae_time_series.pdf\", bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 530,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_experiment(seed, n_iter_optimizer, nn_widths, nn_inputs):\n",
    "    loss_all = np.empty((n_iter_optimizer, 2 + nn_widths.size))\n",
    "    saved_times = np.empty(2 + nn_widths.size)\n",
    "\n",
    "    loss_all[:,0], _, saved_times[0] = run_vi(seed, n_iter_optimizer)  # FVI\n",
    "    loss_all[:,1], _, saved_times[1] = run_vi(seed, n_iter_optimizer, const_z=True)\n",
    "\n",
    "    for i in range(nn_widths.size):\n",
    "        loss_all[:,i+2], _, saved_times[2] = run_vi(seed, n_iter_optimizer, use_avi=True,\n",
    "                                                    input_dim=nn_inputs[i],\n",
    "                                                    hidden_dim=nn_widths[i])\n",
    "    \n",
    "    np.save(\"deliv/timeseries_loss_\" + str(seed), loss_all)\n",
    "    np.save(\"deliv/timeseries_time_\" + str(seed), saved_times)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 504,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seed  2013\n",
      "seed  2014\n",
      "seed  2015\n",
      "seed  2016\n",
      "seed  2017\n",
      "seed  2018\n",
      "seed  2019\n",
      "seed  2020\n",
      "seed  2021\n",
      "seed  2022\n"
     ]
    }
   ],
   "source": [
    "n_iter_optimizer = int(30e3)\n",
    "\n",
    "nn_widths = np.array([4, 6, 10, 20, 2, 4, 6, 10])\n",
    "nn_inputs = np.array([1, 1, 1, 1, 2, 2, 2, 2])\n",
    "\n",
    "init_seed = 2013\n",
    "for i in range(10):\n",
    "    seed = init_seed + i\n",
    "    print(\"seed \", seed)\n",
    "    run_experiment(seed, n_iter_optimizer, nn_widths, nn_inputs)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analyze results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 531,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 2013\n",
    "loss_all = np.load(\"deliv/timeseries_loss_\" + str(seed) + \".npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 532,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAGgCAYAAACg6sNQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAACJD0lEQVR4nO3dd3gU5fo38O9sL+mFQEISakKA0DuiAipVpAj8FCxgwfqiePRgQRCP4jnHdvQcO6hgoWOhNwFRFFCalFBDSEJCerJ9d2beP4aEhGxNdne23J/rypVknmdm780ku3eeyvA8z4MQQgghhAAAJGIHQAghhBASSCg5IoQQQgiph5IjQgghhJB6KDkihBBCCKmHkiNCCCGEkHooOSKEEEIIqYeSI0IIIYSQeig5IoQQQgipRyZ2AMGoT58+sFgsSExMFDsUQgghhLippKQECoUCBw8edFqPWo6awGw2w2az+eTavrouaTq6J4GH7klgofsReOie2Gez2WA2m13Wo5ajJmjRogUAYMeOHV69LsdxyMvLQ1paGiQSylsDAd2TwEP3JLDQ/Qg8dE8cGz58uFv16KdGCCGEEFIPJUeEEEIIIfVQckQIIYQQUg8lR4QQQggh9VByRAghhBBSDyVHhBBCCCH1UHJECCGEEFIPJUeEEEIIIfVQckQIIYQQUg8lR4QQQggh9VByRAghhBBSD+2tRgghhIQQjudh43iU6c0oN9hwoUyPvApDs6/bMkqF4hozeJ4HAOh0erz11psuz+N5HpzVDKlCVXdMKpVhzJgxKCsrwy+/7G1QP6tzZ0wbMxz3De8l2t5wlBwRQgghAY7leNSYbVj55yWXdXmeR1VVFaILOTAM47UYiqpN+O2337Fly2aPzmMYpkFiBAAsa8MPP3xvt/7JEyfw0okT+GRFOr7/1xy0iItpashNRskRIYQQEoBMVhZL918UO4w6r7zyil8fLy/vImZ/uBbfvjjTr48LUHJECCGEBBS9xYavD+SJHUYD/k6MahUV5IvyuJQcEUIIIQHiWGEV9l0oEzuMBsRKjACgc+t4UR6XkiNCCCEkAHz+Wy6sLCd2GA3YbKxoj92xYwZe/38zRHlsSo4IIYQQEfE8j09/vSB2GHa99to/6r6uyjkAbetMyLRRHl3DqqsCa9JBlZACuVwOlUoNuVyOsWPHoqSkBJs2bWxQ/7bbRmD04J64c0Anrw4o9wQlR8Svfv/9d9x7770u60VGRuLgwYN2y0pKSjB06FBoNBr8/PPPUCqVTq81fvx4nDx5EsuXL0fPnj0xd+5crFu3Dh999BGGDh3apOdBCCHe4q3EqENiBG7skACG53Hp0iWkpaU1ayq8xWLBrG1LGxzT552s+7qiogIxMTFNvn6dp6c2/xpeRskREUWPHj0wY4bj5lK5XO6wLDExEbfeeis2btyITZs2Yfz48Q7rHj16FCdPnkTnzp3Rs2fP5oRMCCFet/VkscfnDG6XgM4tIx22qnCcd7rmnP3j+ddff3knMQpQlBwRUSQlJWHkyJFNPn/atGnYuHEjli9f7jQ5Wr58OQDg7rvvbvJjEUKIL5TrLcgt17tdf1rfNGgVgfG23aVLF7FD8KnA+CmHoe3bt+Pzzz/H+fPnGxw3m80uu4nE0K5dO8yYMQO33HKL2KEAAPr06YPMzEwcOnQIp0+fRkZGRqM6Op0OmzZtQnR0NMaOHStClIQQYh/L8Vh92P1p6vf1T4dSJvVhRA1VVFQ4LLt4MXDWXvIVSo5EsH37dowaNQo2m03sUNz222+/YeXKldi8eTOGDx8udjgAhNag+fPnY8WKFZg3b16j8u+//x4GgwEzZsyAWq0WIUJCCLHvYJ7j5ON6Dw1q6/eByXFxcQ7L0tLS/BiJOGjjWRF8/vnnQZUY1bLZbFiyZInYYdQZN24cIiMj8f3338NoNDYqX7lyJRiGwV133SVCdIQQYh/P8zhSUOlW3YcHtxNtxlY4o+SIiMJqtaK6utrhh8HgepNEjUaD8ePHo6amBhs3NpwKeuTIEZw6dQo33HAD0tPTffU0CCHEY+7OTntwUFsfR+I5vd79MVLBjLrVRDBjxgysXLky6FqPZDIZZs70zh43O3fuRN++fR2WDx8+HO+++y50Ol2D43K5HJGRkXXf33XXXVi2bBlWrFiBSZMm1R2vHYg9ffp0r8RLCCHeYHZzUcUZA9pAIlKL0e7dux2WaTQaP0YiHkqORHDLLbdg8+bNWLJkSVANyJ45c6bXxhv17dsXTz75pMPy2NhYrF+/Hs8//3yD4/369cOyZcvqvm/fvj0GDhyIffv24dSpU+jUqRNqamqwadMmtG7dGjfeeKNX4iWEkOaysRy+/N31YGalTAq5VLyOnZtvvlm0xw4UlByJZPjw4Y0SDY7jkJeX1+yFu4JBXFwc+vfv77RObGwsli5tuABZVFTjlVmnTZuGffv2Yfny5ViwYAF++OEHGI1G3HXXXSH/cySEBI+dp0vcqndf/8AcCjBkyBCxQ/AbSo5IwEpMTERiYqLLesOGDUOrVq2wceNGvPjii/juu++gVCobdLMRQojY3FnTaFhGCz9E0jS7du0SOwS/oX+rSdCTSqWYOnUqqqqqsHz5chw9ehRjxoxBbGys2KERQggAoKCy8YxaezokRvg4EudqamocloVTS3z4PFMS0qZMmQK5XI4333wTAK2ITQgJHCYriw3HL7usd1unJD9E49zbb78tdggBgZIjEhLi4+MxYsQImEwmdOvWDdnZ2WKHRAghAICl+91bUbpNvNbHkbi2YMECsUMICDTmiPhV//79kZOT45Nrv/XWW3jrrbdc1nvjjTfwxhtv+CQGQgipj+N5t+pN6J7i40iaZ+XKlWKH4FfUckQIIYT4yGduLviYGBF4S7jUd+edd4odgl9RckQIIYT4QKXR6la9hwJkJWyO4xyWhdsWJpQcEUIIIT6w9WSxW/UCJfHYsmWL2CEEDEqOCCGEEB+oNFpc1pkYQGONRo8eLXYIAYOSI0IIIcTLzpXqXFcCkBDgY40AoFWrVmKH4HeUHBFCCCFetiPniss6abHBsYnr1q1bxQ7B7yg5IoQQQkRwawAs+liL5zjIGEB6dfgTg2tfd8nKAs/awPM8eNbW4Byn1+R58FZr3dd1xzhO+Lj+GMvWHRMbrXNECCGEeNHxy9Vu1ZNKPBuIzfM82OJCmI4eBFtS5LQeU1WF6uhotwd78zzwZAf73WfVyz/zKE5vip42S5THpeSIEEII8RIry+GX86Uu603p2brRMZ7nAasF1au+8EFkzpWXlfn9Md2h37kB2mFj/P64lBwBOHXqFN58800cO3YMLMuiU6dOePrpp9G7d2+xQyOEEBJEcssMbtWL0SgAAGx1JXQ/rvBlSG5Zv2G92CHYZbucL8rjhv2YI71ej/vuuw+dO3fGrl27sGfPHmRmZuKhhx6CXq8XOzxCCCFB5Kczrgdit46UQ7d5Haq+/jggEiMAyM3NFTuEgBL2yZHZbMbf/vY3PPnkk1Cr1dBoNJg8eTL0ej0KCwvFDo8QQkiQ+KuwymUdtqwEmYe3gC1znUQFgluG3yJ2CKIIyuRo7dq16N27NzIzM5Gfb7/JzWq1YsmSJRg3bhy6d++O3r17Y/r06dixY0eDenFxcZg8eTLkcjkAoLi4GJ9++imys7PRtm1gLOlOCCEk8P16wfG4Hd5shunoH7AW5CFewvoxqubpP6C/qI8fMXaqKI8bVGOOysvLMW/ePOzYsQNqtdpp3aeffhrbtm3D0KFDcf/998NsNmPlypV47LHHsGDBAtx1110N6tfU1GDgwIGwWq0YMGAAPvroI8hkQfXjIYQQIpKiapPDMsvZHHAGYVHIm9U6BMhuIW4R830wcsI0SDQRojx2UL37T5o0CVarFZ9++ik++eQT7N+/3269bdu2Ydu2bRg7dizeeuutuuPjx4/HHXfcgX/961+47bbbEB8fX1cWGRmJv/76C8XFxfjf//6HKVOm4LvvvkNUVJTPn1e42rt3L9auXYvDhw+jrKwMMpkMiYmJ6NGjB8aPH48BAwaIHaJd58+fx4YNG/Dkk0+KHQohJADoLTb8cMz+MAzT0T8afJ8ud72liLsUmdlQtMuAJDa+wZR9juNQmZeHqLQ0SCSuO4isVivemf6I3bK3RZpKL7agSo569OiBefPmIS4uDp988onDeuvWrQMAzJw5s8FxtVqNqVOn4l//+hc2bdqE6dOnNzo3KSkJ8+fPR58+fbBlyxZMnjzZu0+CwGAwYO7cudiyZQtSU1Nx++23o23btrBYLDh9+jR+/PFHrFu3DnfeeSfmz58PhUIhdsgNbNu2Df/9738pOSKEAAB+Od+4O43neJj/+rPBsXsjyz2+tnbYGEhbtAQYCRg3Ep2mWL16tU+uG8yCKjl655133Kp36NAhqFQqZGVlNSrr2bMnAODPP//E9OnTsWvXLixcuBCbN29u8CbM8zx1q/nI3//+d2zduhVTp07FvHnz6sZ71XriiSfw+OOPY/Xq1YiPj8ecOXNEitS+I0eOiB0CISSA5JY1nNlsLzEC4FZ3mnbEeMgS/LtytrPGhnAVcu/+Op0O5eXlSE9Pt9ucmJycDAC4ePEiAKBbt24wGo1YtGgRnnnmGUgkErz//vuQSqUYNGiQ08fiXCyd7imO4+o+QtXevXuxdetW9OzZEy+//DIkEkmj5xsVFYW33noLL774Ilq1atWgvKCgAB9++CF+/fVXlJaWQqVSISMjA3fddRfGjLm2UJjFYkH37t3Rt29fvPvuu/jXv/6Fn3/+GUajEWlpaXjggQdw++23N3jcH374AStXrkRubi5qamoQHx+PPn36YOzYsWjdujUKCgpwyy3XZm5kZmYCAE6ePOmLHxVxIBz+ToJJON8Pnufx2b7c644BFjuJ0f9FlDvdGkPepiNUA24CY+c10VOe3pNdu3bZPZ6ZmRmW9xUIweSodm0ijcb+hn61x3U6YXBcXFwclixZgrfeegvDhw8Hx3Ho1KkTPvvsMyQlOc7ebTYb8vLymhynau9eRKxeDdl110i0WmG7riUlENjS0qC7806YbrihWdf56quvAADjxo1zONOw1rx58wCg7udcUFCAp556CiaTCWPHjkXHjh1RUlKCbdu24W9/+xsOHz6Me+65B8C1fXyqq6sxbdo0ZGRkYMaMGdDpdFi1ahWee+45WCwW9O3bFwCwatUqfPbZZ+jRowemTZsGtVqN/Px8/Pjjj9i5cyc++OADxMTE4KWXXsL777+PqqoqvPTSSw3iI/7BcRyuXBGmQbsznoL4VjjfjxMlRlRVGa8d4HngzAm7dY1cJYz2CuQK8LeOF5qVXLwmustb9+SVV14Judc3m83mVq9QyCVHrtS+adYfvJaVlYXPPvNs7xiZTIa0tLSmBbF9O5gZM8DYbI2KVE27ou8dOgTthg3gN24Ehg9v8mVOnz4NABg7diwiIjybhfDGG2+guroa7777LkaMGFF3/NFHH8XEiRPx7bffYubMmWjV6tr+QDk5OZg9ezYeeeTaYMP27dvjueeewx9//IFJkyYBAH7++WdotVp89dVXkEqldXVHjRqFt956CyaTCRkZGcjIyMCSJUtQVVWFadOmNelnQJqn9j/ZNDcHmxLfCuf7sSX/AqKjrw3HMB8/DKjsv4pHR0Y3OqYdMwXS6Bivx+WtezJmzBiHDQ3Byt3hMiGXHEVGRgKAw9Wta4/X1muOJv/SffklYCcxCnSMzQbmiy+AW29t8jXKy8uh0Wg8ngVoNBqxZ88eJCYmYtSoUQ3KoqKiMHbsWHz44YfYvXs37r777msxMwzuu+++Bveqa9euAIQ1rWqPy+VyGAwGHD9+HD169Kir27t3b7zxxhsNXmRqE+tweyMIJBKJpO6DiC9c70f9f7LZmmrAQRdUttLYaAPYiDGTIY2J81ls3rgnnv4DG0pC7jdZo9EgMTERRUVFYNnGC23VduW0adPGz5ERQMjanfW7O5KbmwuWZZGRkWG3vEOHDgCACxcuNDiekJAArVbb4JhSqQQgNK/WeuSRRyCRSHD33XfjnnvuwUcffYSjR4+GbX87IcS5+q9jnMkE64UzDut2UTRcAynq7od9mhiR5gu55AgA+vbtC4vFgmPHjjUqq10bqV+/fv4O65oZM4BgnAknkwHXLY/gqaSkJBiNRpSXezaltbbFT+Wgybr2uNHYsFff3WUAbr31VqxZswbjx4/H+fPn8c4772Dy5MkYNmwY1q8PzA0ZCSHiybmiq/vacvq4w3rt5GYomWuJlGbIbY1akcTk6WtxuAjCd2jXpkyZgo0bN2Lx4sV4//33645XV1djxYoViImJwciRI8UL8JZbgM2bgSVLgPPn6w7zEPZ6UyqVCJw/navatRMSo2aMNwKEbqrz589jz549GD9+vNO65eXliIsT/ruqbd6tqKiwW9dgEHbCvr6VyBNZWVl4/fXXwfM8cnJysGvXLixbtgzvv/8+0tPTG81uI4SEp0qDBXvOlgAAzMedL+0xRN1wiIcspYljVX1kw4YNYocQkIImOSooKGjQElSb7e7Zs6fuDTQlJQXZ2dkYOHAgJk2ahDVr1mDWrFkYOXIkDAYDvvrqK5SWluLdd9/1ypijZhk+vFGiwXMcivLykJaW5rPFvsQ2YcIErFq1Ch9//DFGjRpV18V1PZ1Oh0mTJiElJQVffPEF2rRpA7lcjgsXLtidbVA70Lu2e605GIZBp06d0KlTJwwfPhxjx47F5s2bKTkihAAAjhZWg+cB87E/XFeuR3vbHWDqTfgIBF988YXd4zExMX6NI9AETXL0+++/4/nnn290/JVXXqn7esKECXjjjTcAAP/4xz/QuXNnrFy5EvPnz4dcLkf37t2xcOHCuunbxP969+6NqVOnYsWKFXjqqafw5ptvNmrtqaqqwuOPP47CwkLcc889kMlkkMlkGDZsGLZs2YKNGzdi3LhxdfVramrwww8/QKFQYNiwYR7HVFhYiIcffhjDhg1rtOBkbfJWPxmrnc1mMpkcdvMRQkLXibwrsJ5xvb7ZbZrquq8jx90FSWTgbUe1c+dOu8c//PBDP0cSWIImOZo4cSImTpzodn2JRILp06fb3SKEiOvFF18Ey7JYvXo1brvtNowbNw4ZGRmw2Ww4ffo0vvvuO+j1ejzzzDOYMWNG3Xm10+/nzZuHnJwcZGVloaysDMuXL0dxcTFefvnlBvvluSs5ORktW7bExx9/jIsXL2LQoEHQaDQoKirCqlWrIJVKG2xUnJqairy8PMybNw+ZmZm4/fbbna6JRQgJHRfLDbBdvADejckaLaXCpI+oux8OqHFG7ujfv7/YIYgqaJIjEjqUSiVee+01jB8/HmvXrsXOnTuxfPlySCQSJCcnY8KECZgyZUqjLrLWrVtj9erV+N///of169fjiy++gEajQXZ2NubOnYubbrqpyTF98MEHWLx4MbZs2YJ9+/bBaDTWbYL77LPPNhjA//TTT6OoqAhbt27FkSNHGqyaTQgJbZuPXgJnMris10VhAsMAjFIVdIkRQDO6Gb4p86rD3PCrY4V27Njh1etyHIe8q2OOwm29kEBF9yTw0D0JLOFyP3iex56zpTh+/Axslwtc1r8vShgXGzFqEqRxCb4OrwFP7omjxC1UUwN337+p5YgQQghx4UyJDidzi91KjLQSDoxciYgxd0KiDdyFFEM1AfIGSo4IIYQQF/afyoflrHubTI/VVCFqysM+jqj5Dhw4IHYIASt020AJIYQQL+AMOpQf/8vt+vEjgmPZj6VLl4odQsCi5IgQQghxgDPoULb2G7frt5VbIGvRynXFAPC///1P7BACFiVHhBBCiAOWMyfxk9H9cUOj75/qw2j844EHHhA7BNFRckQIIYQ4kHPoLxTZ5G7VHZCVCmkw7pt5nVdffVXsEERHyREhhBBiR0lFNXZ70GrUIqOjD6Pxn1atgqNb0JcoOSKEEEKuw/M8vl2+xe36Em0kUlvG+TAi4k+UHBFCCCH1cCYjSr7+1KNz+g7uBakkeFbCZllW7BACGiVHhBBCyFU8a0PNmqX4xah1Xfkqdbc+6N++hQ+j8r6TJ91bsylcBf/IMUIIIcQLeKsV1SuX4IxFiUs2hVvnyJJTcU+/NB9H5n3ffvut2CEENEqOCCGEhD2e51GzdhlW6mJg5NzvVLnjpu5QyaU+jMw31qxZY/f4ww8H/sre/kDdaoQQQsKe+fDvOGtkPEqMlF16IDla5cOofCcnJ8fu8fHjx/s3kABFLUeEEELCGltViQOHcnDI7NkmsYxU6nBX+2A1ZMgQsUMICNRyRPxq7ty5yMzMxO+//y52KHYNGzYMmZmZYodBCPGjQweO4JBZ49E58rS26J4S45uARBQR4VmCGKqo5Yj41bRp03DzzTejY8fmL5ZmsVjwzjvv4IsvvkCrVq2wc+dOL0ToOzabDd999x1WrlyJixcvwmAwICUlBSNGjMDDDz8Mrbbh7BiDwYBPPvkEmzZtQmFhIdRqNbp164ZHH30UvXv3FulZEBJajuRX4tczVzw+TxIdhz5psT6IiAQCSo6IX2VnZyM7O7vZ1zlx4gSee+45FBQUgOM4L0TmWyzLYtasWdi7dy8GDhyIp556ChKJBNu2bcNHH32Effv24ZtvvoHs6tYDFosFM2bMwOHDh3HHHXfg0UcfRVVVFb766ivcc889+N///oehQ4eK/KwICW6czYafdzWtFXtKr9ZBta5RfcHwmik2So5I0Dl//jymTJmCTp064bvvvsNtt90mdkguff/999i7dy9uu+02vP/++3XHp06digceeAB79+7F9u3bMXLkSADAN998g8OHD2PWrFmYM2dOXf0xY8Zg9OjRWLBgAQYPHgyFwr3pxoSQxgp27QRXU+3xedKYeMRqgvdvj9Y4co2SI5HYLufDcj4HnK6m7hjPcWDKy6E/GQdGEljDwSQRkVC0y4SsVetmXWfu3LlYt24dli5div79+8NiscBkMrk8TyqV1nU7mUwmzJgxA//v//0/yOXubQjZXC+88ALWrVuHxx57DLNnz4ZOp3Prvy+NRgOZTAaZTIbRo0fj/vvvb1Tn5ptvxt69e3Hq1Km65GjdunVgGKZR/YSEBIwdOxZff/01fvnlF2o9IqSJeIsZRZcKAbi/2CMAQCKBND7BJzH5y2+//Wb3eGpqqp8jCVyUHInAdjkf+p82AXzDN1ee54GqKrC8LeBmQLClxbBePA/t0FHNTpDqW79+PZ5//nmX9fr164dly5YBADp37ozOnTt7LQZXlixZgnXr1uHee+/F7NmzAQDjxo1DQUGBy3Nrk8Bx48Zh3LhxdutUVwv/uUZHRwMQxhrl5OQgPT0dcXGN92rq2bMnvv76a/z555+UHBHSRNWrvsBpa5RH5zAyGeRtOmBwl+Bb9LG+f/7zn3aPP/jgg36OJHBRciQCy/mcRolRUOA5WM7neDU5GjJkCJYuXeqyXlSUZy9i3rJs2TKsWLECEydOxAsvvFB3/O2334bZbHZ5flZWltNyg8GANWvWQC6X49ZbbwUAXLp0CTzPIzk52e45tTtmX7x40d2nQQipx3TkAACgjHX/LVASEQV5245gGCArKdJXofnFmTNn7B6nafzXUHJERJWYmIjExESxw7Br48aNWLRoEYYMGYKFCxc2aM3r0aNHs69vsVgwZ84cFBQUYM6cOWjdWkg69Xo9AKFLzp7a7sXaeoQQ97HVlTD/9ScsvGet84p2wgzbMV1aQSYNrGEP3jJ48GCxQwgYlByJQNEuE9aL54Ov9YiRQNEuPNYA2rdvH5577jl0794dc+fOhVTq3e0BysvL8dhjj+HQoUO4//77MWvWLLfP5Xneq7EQEi54qwXVP6zAfrMWORal2+cpu/So+7plVHCuiO0OmuBxDSVHIpC1ag3t0FF2B2SDkUEaF7oDsq9nsVig0+lc1pPL5YiM9F9T9uOPPw6r1YozZ86guLgY7dq1a1BeVVUFlmVdXicyMrLRoPFTp07hscceQ1FREZ5//vlGg65rn6ejlqHa4/78eRASCqpXfo4jFrVHiREjkYK5+s/RhO4pQTt9n3iGkiORyFq1bpRocByHsrw8aNPSIAmw5MhXmjIg2x+GDBmC0aNHY/bs2Xjttdewdu1aqFTX/mOcMGGCRwOya+3fvx+zZs2CXC7Hp59+arcZOzU1FRKJxOH18/PzAQBt2rTx8FkREr4Mv/4EADhqVnt0nrzdtQVrEyPcT6oClcViETuEoEDJERFVoA7I/s9//gMAmDlzJhYvXozXX38dCxcurCtvyoDsI0eOYNasWYiNjcXixYvRtm1bu+eoVCpkZ2fjyJEjuHLlClq0aNGgfP/+/QDQIOkihDhmzbsA64XTHp8nTUiCROPhVP8At2jRIrFDCAqUHBFR+WNAdmFhIYxGI1JTUz3uU3/qqafw66+/YsWKFRgwYABGjx4NwPMB2ZWVlXjyySehUqmwdOnSusHXjkyePBlHjhzBkiVLMHfu3Lrjly5dwqZNm9C2bVv069fPoxgICUdsRRkMP28FywM/6KPdPk+ekt5gPaOerWN8EJ3/bdiwQewQggIlRyToHDt2rFGXk9FoxObNm+u+z87ORkpKCgDg73//O/bv348ff/wRGRkZHj2WTCbDCy+8gCeeeAIvvfQSunTpgvT0dI9j/vjjj1FcXIzx48fjr7/+wl9//dWoTlxcXF3CM3HiRGzatAmff/45ysvLMXjwYJSVleGLL74Ax3F4/fXXvT5InJBQY83PhWH3FgDAVzWN1wxzRhLTsH6v1NDYR+3AgQN2jz/88MN+jiSwUXJEgs7XX3+NdevWNThWXl5et0AjIDQdT5w4sUGdpo7jatGiBV577TU88cQTmD17NlasWAGl0rOxB7XJ0HfffYfvvvvObp3646qkUik++ugjfPrpp/jhhx+wceNGaLVa9O3bF0888QQ6derUpOdCSLiwXblclxj9ZrK/LIYjkogoMFen6zMMg2EZiSE/ELv+Om4EYHiaF+yx4cOHAwB27Njh1etyHIe8vDykhdGAbH/gOA69evXCtm3bPO7Co3sSeOieBJZAvB88x6H6208BADpOgjW6GI/OV3buDkYmw5D2CUiP00CjCK52BGf3xNHuCyzLBsz98yV337+D644T0gS//fYbWrVqFbCLTRJCvKegtBqb1myCgfOsG62WqltvAED3lBhktRRnZX4xhENi5AlKjkjIq6ysxMsvvyx2GIQQH6uurMaaVVsANP+Nvn+bpiVXgYw6itxHyREJebUzzAghoYmtLIf55FEcOXkRgGfji+qTxsYDAAa0ifdSZIGluLhY7BCCBiVHhBBCghan10G//UfwZhMOmZve2sPIFZC1EDZ1zk4Oze60d999V+wQggYlR4QQQoISz/Mw7tsF1mTCFbbpb2fyNh0g0UbUbRPiaNBysHv77bftHr/lllv8HEngo+SIEEJI0LFdKYJ+2/c4b1XgZ2PTW4yk0XGQRl1bHLJzCA/Ctlqtdo8/8cQTfo4k8NHwdEIIIUGFt9lg2LUJVh742RjRrGvJ0q5t46NVyJCd7P4q2qGCxmU2Ri1HhBBCgoo17zysFgu+8XDV6+vJU9ugtgftpg6JSIvVQK0Iv5Xn5XK52CEEHEqOCCGEBAWe53GssBo/788H28zESNEuA5KISADAQ4Pahuw4I9I0lBwRQggJCr+cL8OJompwuppmXUfWqnVdYtQ9JSYsEqMLFy6IHUJQoTFHhBBCAp7FxuFEUTV4lgVvtTTrWrLEpLqvM1o0b8xSsPh//+//iR1CUKHkiBBCSEDjeB6bTxYBAMzHDzfrWoqMznVfR6rkiNUomnW9YLF+/XqxQwgq1K1GCCEkoFQYLDhXqoeE59CypghHL5bhUrUVbGV5s67LyORglOq672/NbNHcUIPeK6+8InYIAYmSI+JXc+fOxbp167B06VL0799f7HAaGTZsGAoKCpCTkyN2KISEHRvHYevJK8ivNIDnAcvxQ+A5zmvXV2R1q5udlhSpQkKE0mvXDlZz5swRO4SARMkR8atp06bh5ptvRseOHZt1ndzcXHz88cc4ePAgLl++DK1Wi27duuHBBx8MyKTLntLSUowZMwaVlZUBmywS4k+nimqQX2kAp9fBcs77/6DUH3c9qnNLr18/GEVEhMeYK09RckT8Kjs7G9nZ2c26xq+//opHHnkECoUC06ZNQ9u2bZGXl4evvvoK999/P95///2gWA5/wYIFqKysFDsMQgLGrxfKwJlNsF446/VrS2MabiarkNGQW+IYJUckqPA8j3nz5oFlWSxduhSdO18bXHnjjTdi6tSpePfddwM+OdqwYQO2bduGzp0748SJE2KHQ4hf2VgORdUmcACSo1SwcTzkUiFZYa8UgedYrz+mLDW97utwWwW7tLRU7BCCDiVHIsmvNOJ0cQ2qzba6YzzHoby8GnGVhWAkgfVfTZRShoykSLSOUbuu7MT1Y44sFgtMJpPL86RSKbRaLWpqajBgwABEREQ0SIwAoEePHoiJicHZs2dhs9kgk3nv1/uFF17AunXr8Nhjj2H27NnQ6XTg3BgLodFoGsVRVlaGV199FUOGDEH37t0pOSIhjed5nCyuQW6ZHmq5DB0TtdiXW44Kg+W6egBbdgVsRZnXY5CnXlvkMUatQLeU8EqOpk6dKnYIQYeSIxHkVxqx6UQReJ5vcJzneVQZbbDpzAG3KNmVGuBcmR6jOrdsdoJU3/r16/H888+7rNevXz8sW7YMUVFReO211+zWqU207CUkzbFkyRKsW7cO9957L2bPng0AGDduHAoKClyea28s0SuvvAKr1YpXX30Vq1ev9lqchASiowVV+P3itVlmZ0rsL+BoK8wDW1bi9cdnZHJIIoXNZEd2bolWUaq6VqpwsXPnTrFDCDqUHIngdHFNo8QoGPA8j9PFNV5NjoYMGYKlS5e6rBcV5Xqn7FWrVsFkMmHKlCneCA0AsGzZMqxYsQITJ07ECy+8UHf87bffhtlsdnl+VlZWg+83b96MLVu2YOHChWjVqpXX4iQkUB0vqnZazvMAV17im8RIpYaiTQcwMhmSIlVIi9V4/TGC2QcffCB2CAGLkiMiqsTERCQmJjb7Ovv27cM///lPtGzZEk899VTzAwOwceNGLFq0CEOGDMHChQsbtOb16NHD4+uVl5dj4cKFGDhwIDVzk7Chqz90gAfA2sDUa9m1XbrQ7PWLHFFeXfCRAYM+abE+eYxgdt9994kdQsCi5EgEGUmROFemD7rWI4ZhkJEUKXYYjaxZswbz589HbGwsFi9ejPj4eNcnubBv3z4899xz6N69O+bOnQuptPk7db/66qswGo149dVXm30tQgKdjeVwKL+y7nu2shy2wnzwNiskai3k6e1gu5wPtqrCJ4+vaNsRGS0iEaWSoU2cFnHa8FgJ2xMajUbIWANsjGudL74ARErgKDkSQesYNUZ1bml3QLbMIkNchDJkB2Rfz2KxQKfTuawnl8sRGdkwMWNZFv/+97/x+eefIysrCx9++KHXuqoef/xxWK1WnDlzBsXFxWjXrl2D8qqqKrCs6xk1kZGRkMvl2LZtGzZu3Ih58+YhNTXVKzESEmjK9BaU6c1oEaHErxfKkV9pAABwZjOsl3KvNh0BnFEP86ljPomBUSghS05FRttk3NQhIeDGb/qb04kjgf6zuf9+QCYDpk3z+0NTciSS1jHqRokGx3HIy7MhLS0ZkgBLjnzF0wHZtWw2G2bPno3t27dj1KhRWLRoEdRq746FGj16NGbPno3XXnsNa9euhUqlqiufMGGC2wOyMzMzsWDBAnTt2hXDhw9HUVFRXXltYlheXo6ioiLExcVBoaD/cElw0VlYfPrrBYeJCFtWUpcY+ZJEGwlF+wykx2kwNCMx7BMjADh7tvGaUU8BeMfvkTTR9OmUHJHw05QB2RzH4dlnn8X27dvxyCOP4KmnnvL6i+B//vMfAMDMmTOxePFivP7661i4cGFduScDsk+ePInS0lKUlpbi5ptvtluvdpwUrZRNgo3RymLTmWpERTueNMGV+2edHXmb9hiW0QIdEmnV51p33nlng+9XAPDelJXQRckREVVTBmQvXrwYGzduxAMPPICnn37aZf3CwkIYjUakpqZ63Crz1FNP4ddff8WKFSswYMAAjB49GoBnA7IzMjLw0Ucf2S1bv3491q9fjzlz5iAjIwMZGRkexUeIGDiex+H8ShzMqwDP8+DhuFWI5+GTRR2vd1P7eHQe3D7spum7cuzYte7LW0GJkbsoOSJBpaysDB988AE0Gg26dOmCzZs3263Xr18/xMXFAQD+/ve/Y//+/fjxxx89Tj5kMhleeOEFPPHEE3jppZfQpUsXpKenuz6xntjYWAwdOtRu2V9//QVASLaoxYgEAp7nUW6wQquQQiVvPBGB43lsOF6Ey1VGx9ewWACpFIxUCt6o92W4AIC4lGR0u3UQdaO5sFXsAJpi925RHpaSIxJUzp49C4NBGOTpbDdpe91TTR3H1aJFC7z22mt44oknMHv2bKxYsQJKJe3mTUJPpdGKjccvQ2e2gQGDXqkx6JUa0yDpKKwyOUyMeJsN1tyz4Ax6gGEgjYnzyYrXtSSR0ZCnt8eInq0pMXIhuOZGX/Xrr8DAgaI8NMMH23zyADB8+HAAwI4dO7x6XWFAdh7S0tLCZkC2P3Ach169emHbtm0ed+HRPQk8dE985/ujhSiuabidz7jsZCREKKAz2RClluO7I4Uo1ZvBWyzgaqrAy+TQsRyiY2Nhu5TrszWL7MkaMhD92ichWi3322MGg/p/I1KpFDcB2OXqpLNngfbtfR+cyNx9/6aWIxLyfvvtN7Rq1cori00SEkoMFhtOX9FBb7GhZZSqUWIEAD8cK4RUwoDleKjlUhitLDhdDay5Z8HXThNnJOCjov2aGHXq0x23dm3tt8cLRr///jsANxIjaiNphJIjEvIqKyvx8ssvix0GIQHFYLHhqwN5dd8fv+x4mw+WE948jVZhYLWtqPBaYgQARgMsJw77JM5GJBKouvZE/660Xpgrn3zyCd5zVclmc1UjLFFyREJe7QwzQoiA4/kGiZE7eI4HW14CXq8HZ3C9cKuvMHIFWkapqCvNDV988QU+d1bhs88AL6z+H4ooOSKEkDBTXN24+8wZnmVhPnEkILpfMlITMLxLS7HDCArZrio88IA/wghKlBwRQkiYOXip0u26vNUKy5kTAZEYRUg43DqsHw3Ed9NRZ4VGx0sxEEqOCCEk7HjScsSWl4IPgHEpCVIbbh3QBRLqBnILU1PjvEK97ZBIY5QcEUJIiKkxWVFmsCBBq8SRgkrklhmgVcrQOzUGBgsLzk4rEM/xYCSN1wqyFRf6I2SHlFndkKpgMapvBzBy2nfQXendujkutFj8F0iQouSIEEJChJXl8ONfl1Gqa7zvn95iw6YT1zY95nmAqyyHrfASeFZoGZJERkOWmATOaACjUEASGe232O1Rdu0JRiLBjb3bg5HT25XXyGkwuyv020YIIUGG53kcyq/E6Ss6yKUSZCdHIS1Wg28OXoKt/hR7J9jiQtiuXG5wjKupgqWmyhche0waGw9GIsEtmUmIUNJblSecdT3+uWoVevkxlmBFv3GEEBJkjhVW42BeRd33u86UeHQ+z/ONEqNAwjCAunUauqfGom28RuxwQkqH224TO4SgQMkRIYQEmSMF9lt3OL0ObFUFGIlUaHmxswcgZzbDlp/r4wibrq/KgK6D+0PToT2kdsZAERfOn3dY9A6Ap6Oi/BdLEKPkiBBCgozR2nj2GFtVCWve+bop92x5KRTtMxskSLzVCuv50+CtgTcgV8oA7TUc+t55J2QardjhBC8n+6PNAfC0/yIJapQcEUJIEDFd3cLjemzplQZrEfE2K9iqcshatAJvtcBWfBlseam/wnSbIrMLpvVqDaVMAkVEhNjhBLcAWIsqVFByRAghQeRcqd7ucU7feF0bW1Eh2NIS8Darr8PymDKrGxi5HA8NaguGoe4zr3CyOCbNT/MMJUeEEBJELKx7s9FqBWJiJGWAVvGRuLFDIiVGfmIDsGHDBrHDCBqUHBFCSBAxOuhWCxbTI8vBAIjNThY7lNCit9+iCAB/v/p51KhR/oklBFByRAghQeRkUXWjY8Ey1GSougZSBpBE0Iwpr3MyXutfVz9TK537KDkihJAAxPO83Tczlqs36JoHeIMOlnM5/gytSW5U65AmF7r4tLfdIXI0hDhHyREhhAQQS+5ZmI8eBG8yQta6DdT9hoCRCS/Vp4prwJaVgK0oB8/awJvd30DWn9LlFpSyMug5YYBwlsKENjJh+QB5WjtI1LSwo1dt3eqwiBZFaBpKjgghJEDYSq/A+OvOun4y64XTYCQSqAfchCMFVfj1j1OwFuSJHKVz90WVAwA4HqjgpIjQqqE0GSCJjIa8XQaUWd1FjjAEjRjhsMhw9fO///1v/8QSIig5IoSQAGE5c6LRACLLuVOwFFzCL3wa2JrG440CSaz02mBxCQMk9+wJVXZvESMitR599FGxQwgqlBwRQogf8TwPrqYKvF4HaUISmHo7pFvP2x87dEVvhk1f7K8QmyRCwmKYuuFaSxJa6dr35s51WFR/xJpWS/fCE5QcEUKIF/EsC9vlS+BqqiFrmQJpbPy1Mp6H6cAvsJw5XndM0bEL5GltASc7qRs4x4v7ieGuyApIwUMCoIKVoMymQ1uNGrL6e6FJpJCnOd7KgnjJP/8pdgQhiZIjQgjxEp5lYdizFbbCq+OCGAbqQcOgaNMBAGC7nN8gMQIAy5njjY5dL1AmYMdIWYxQV0PBXOv6i5WykEgskDLqBnVVvQc2aBUjPuBkDYd3/BhGKKLkiBASEjijAdZLF8BbzJBGx0GWnArGSWuML1jOnryWGAEAz8N8eD/k6e0BjoPhp41NvLL4Cxl1UxrRWWGCkmkcC9+jP6IG3gi+pgpseSmk8YmQRsX4P8hwc+ONDovm1Pv6008/9X0sIYaSI0JI0OMMOui3/QhO13DAsmbIbUKXlbcfT6+D7XI+GLUGspbJYKQyWM6egungL3bq1oCrKINu05omP14FJ/5LdReFqUGLUR1GArRIAQBIo2MhjY71c2RhbO9et6rNnDnTx4GEHvH/4gghpJms5880SowAwLB3OyIn3QOJUuXWdXieB1tWAt5kgCwpGYxcce0x8nNhOXkMtiuFDc6RtUiG5qbbYPrjV4fXbU5ixPPAEbPadUUfSpVZoGB4yFLSAY6D7fKlujJFRheAus/8z0mXWpfrvpc42ZCW2EfJESEk6JmO7LdfwHOwnj8NZVY3h+faSothvXheGPfDXpuKzsgV0A4fC0lMHEyH98Ny6qj9868UwvTnbz7b4LWK82/X4PXS5RYMUgn7dqm69YEkOgbW82fAVpZBmtgS0tS2wKVLLq5CvM5JwnPCj2GEKkqOCCEhzVZc6DA5shZchGH3VoBvvNM9b7XAeGAvJBotrJcuOH0My7lTXonVnlNW91q9fGGQWo+OcjMAQJacBklsPBiGgaJjVl0djmv8syMk2FFyRAgJbU7evC0nj9lNjGqxZVfAlvkiKPexIo3F7iA3o4NMSIzU/YZA3i6DNi4NFFbHrZTJ133/888/+zaWEEXJESEktDlJfmzFBX4MpGn83S7TT6VHqsyKmFatIIltD2Xn7rQXWqBRKBwWXb7u+8GDB/s2lhBFyREhJKTx9cYRNTjuZEBrIDHxzRtM21NpQAkrRyknhZbhkCE3I9emwGVb40HUA1R6dE6JhXboqAaD0Unwota+pqHkiBAS0tiSItSsXwl170GQtWoNzmiA8bfdDdcjCmCFdpIYdw1R69BObgFganC8g9yMdUn9oTt7pu6YlOHRVmmDduhoWrwxkDlJ6tP9GEaoo/l9hJCgxhn0rutUVUC/ezM4gw6GXZuCJjFqjiSZDakyi90yCQOM7p+JyPg4AICK4TFUrUNklx6UGAW6ESMcFoX+b7X/UMsRISSoGfZsda8iy8KS8xfY8lLfBuRFZt7zLpEshQmtZFa0klohc3C6IrMroqPUuG/Szai8cAFqfQXkLVpB2jKlmRETn9u2ze2qJ07QpP6mouSIEBK0OF0N2LIrbtc3nzjiw2i8rykbzmbKTYiWOh/GrejYGQAgkUoR16FDk2IjgaWnnWOZmZl+jyNUULcaISRosRUiz7P3sYu25g+KVnTqBlmr1gAASWQ0NENupS0+gpWT7UIO+y+KsEAtR4SQoMTzPKwXz4odhk+dtDR/AUhlVjYkmgjwHAeGtpEIbkOGiB1B2KC/FEJIULIcPQjrxXNih+FTliaMOYqQNOxSk2giAIASoxA2w86xKVOm+D2OUEJ/LYSQ4MOxsJw+LnYUomOkQuM/IxNmmHVUmCGtl09JIqPFCIv4QlGRw6Iv7Bx74403fBZKOKBuNUJI8KmuBG+1hPUCd7IWrTBuaE+czC9DDSdFK7YGmSd3N6hTO/CahIBWrTyqHs5/G97gteSIZVns3r0bv/32G86fP4+amhowDIPo6Gh07NgRQ4YMQf/+/b31cISQcMaG92anEm0EpIlJaB2nRes4bd1xtkMSLDnHwRv1kKW2haJdhohREhK8vJIcHT9+HM888wwuXrwIoPGy/Hv27MHixYvRvXt3vPnmm2jdurU3HpYQEq6c7JcWSiSaCEhiYsGbTZCotWCu7nHGqNSw1zAgjY6Fut8Nfo6SiMnekp1jxozxexyhptnJUXFxMWbOnImqqiooFAr06tUL7dq1Q1RUFHieR1VVFc6cOYMjR47g8OHDmDlzJtauXYuIiAhvxE8IISFHmtgSysyukOfxDgdSJ0er/RwVEc2HHzosstk59uWXX6KmpsZ38YSBZidHS5YsQVVVFW6++Wa8/vrriIuLs1vv8uXLeO6553Dw4EEsXboUjz32WHMfmhASrrjQbTnS3joOshbC+BIm/7zDeqkxlByFDQ/fL2NjYyk5aqZmz1bbu3cvIiIi8PbbbztMjACgVatWeO+996BUKrF9+/bmPiwhJJwFabeaNC4BjEwOaUJS40JGgshJ99YlRtcPT7heZlKkL0IkQeQOsQMIYc1uOSosLET37t2h0Whc1o2NjUV2djbt90IIaZ7SYrEj8Ji8TQdoBg8Hz/NgGAZseSkMu7eAM+jAyJVQD7gJEtW11iDWSXIUq1FAJZf6I2wiNictQD/YOUaz1Lyj2ckRx3GQSt3/I1UqlbDZ7PWSEkKIe5jcM0B0cK3hox44FMC1Ny9pXAIixt8N3qAHo9Y0GlvkrOcwLdb1P6MkRERFeVR961Y3N2ImTjU7OUpKSsLx48dhtVohl9sbN3+N1WrFiRMn0KJFi+Y+LCGEBCxFRhfwNhu4mirIklOh7NzD7sBqhmHAaO1PTmE5xy1HnIsuNxK+hg8f7rJLlrjW7DFHAwcORFVVFRYuXAiWZR3Ws9lsWLBgASoqKjB48ODmPiwhJEzxbDNbnr3U7aAdNqbBCtQSTQTUA25G5B13Q933BmgG3oyI2+6AqmuvJm3dUWG0OCxzkjeRMJHp4Dh1q3lHs1uOHnzwQfzwww9YvXo1fvnlF9xyyy3o2LEjoqOj66by5+TkYPv27bhy5QrUajUeeughb8ROCAkznNkEw/b1zbqG5sbbYNi9pVnXUA+4CbJWrREx5k6wJcWAVAppfAuv7l+WU6xzWEYtA2Hi8GGHRaf9F0VYanZylJqaivfffx/PPPMMCgsLsWzZMrv1eJ5HfHw83n77baSkpDT3YQkhYYTneZiP/QHzsT+anRjIUtIhS06DrTCvwXFGpQZvMrp1DUX7TsI5UhlkLX3zenamxPFAXEqNwkTPnh5VHzdunI8CCT9eWSH7hhtuwObNm7F69Wr8/PPPOHfuHKqrqwEA0dHRyMjIwI033oiJEyciysPBZYQQYsn5C+ZjfzT7OoqOXcAwDNSDhsK0/2fYLhdAEhkFVc/+dUkOp6sBW1UB3qCHcf+eRtdQ9R7U7Diai8YcEXu+/fZbsUMIGV7bWy02NhYPPfQQdZkRQryKMxlh+uPXZl+HUWmg6NQVACBRqqAZcmvdtPr6JBGRkEQIawjZigtgvXju2jXkSsjT2jU7luai3CgMOLnJKgfH3VlSh7jHa8kRIYT4gjX3bJPPlbfpCGl8IgAG8rQ2kGgazgxzNXhVPXAoJNpI2AovQRIRCWV2b0g0Wqfn+ENihFLsEIivzZjhsMjsxzDCldeSo3PnzmH9+vU4cuQI8vPzodfrwTAMIiIikJ6ejt69e2Ps2LFITk721kMSQsJAc1qNJBGRUHbKbvL5jFQKVc/+QM/+Tb5GU7gaV9UhUfwEjfjYl1+KHUFY80py9Oabb+Lzzz8Hx3GN/qhLS0uRm5uLPXv24L///S9mz56NBx54wBsPSwghTslbtxE7hCaxsI5XgGyXEAGljFbHDldfODi+efNmf4YR8pqdHK1atQqfffYZACA9PR2DBw9GWloatFrhPxudTofc3Fz88ssvyM/Px5tvvonk5GSMGjWquQ9NCCEOSaLjIIlLEDuMJsktMzgsa59ArUYhz0nLoaPOthEjRvgmljDV7OTo22+/BcMwePrpp/Hwww87rMfzPN577z18+OGH+PLLLyk5IoT4hEQbCWl8C6h6DwzaBfH2ni91WBakT4l4QkljysTW7OTowoULSE9Pd5oYAcLAx9mzZ2PDhg3Iyclp7sMSQohdkePvFjuEZnO2dYhC6r2FJkmAslrFjiDsNfuvTCaTIT093e36bdq08WijWkJI+OKd7b4appIiHU3kJqHuPgfHP//8c7/GEQ6anRx17NgRhYWFbtcvLi5Ghw4dmvuwhJAwYD5yQOwQAgrDMJBKqF8tpJWUOCxa6uD4ffc5SptIUzU7Obr33ntx5swZ7NnTeCXZ6+3btw85OTm4++7AbfbesGEDMjMzsXbtWrFDISTsWc57toOUMqu7jyIJDN2So11XIsGtRQuPTwnWsXWBrNljjkaOHImXX34Zzz33HCZNmoQRI0YgIyMDKpXQ9GuxWHD+/Hls3boVy5Ytw1NPPRWw+79UVFRg0aJFtMooIQGAt9nAmxzP2rqeLDkNyuzePoyIEBIuPEqOevXq5bDMYrFgyZIlWLJkCQBALpeDYRhYLJa6OgqFAl999RW++eYbt1qa/O21117DyJEjsXPnTrFDISTs8WaT23UjRk68uhI2IaEp1sFxiYQG6PuCRz9Vg8Hg8MNms4Hn+boPi8UCs9nc4JjZbEZpaSlKnPSpumPt2rXo3bs3MjMzkZ+fb7eO1WrFkiVLMG7cOHTv3h29e/fG9OnTsWPHDrv1d+/ejT///BNPP/10s2IjhHgHbzK6rMMoVdDccEvYJEbUexLinMzkrnRwvH4DBPEej1qOli51NBzMfSzLwmxu2s4w5eXlmDdvHnbs2AG1Wu207tNPP41t27Zh6NChuP/++2E2m7Fy5Uo89thjWLBgAe666666ujqdDvPnz8err75at3glIURcnJOWI75NR0QMHQGpWkPjLUjo6NTJ41No9rdveJQc9evXz1dxuGXSpEmwWq349NNP8cknn2D//v12623btg3btm3D2LFj8dZbb9UdHz9+PO644w7861//wm233Yb4+HgAwL///W/0798fQ4YM8cvzIIS45rRbLbEVJCp1SCZGrvZVI4T4XlB1Vvbo0QM//PCDyyRm3bp1AICZM2c2OK5WqzF16lQYDAZs2rQJAHDgwAHs3LkTzz//vG+CJoQ0CW+0Pxhboo0EEjyf0RMsnKVGoZcKEnd0FDuAMORRy9F///tfrz3wE0884fE577zzjlv1Dh06BJVKhaysrEZlPXv2BAD8+eefmD59OtauXYuqqqoG+9JUV1fj1VdfxbZt2/Dhhx86fBzOywvUcRxX90ECA90T8diqKuy2omhGTgB3uShk7wnL8Q5bj6JVsoB73vQ34iWrVjlsrTjr4PjJkyft/tzpnjSfx8mRt5qxm5IcuUOn06G8vBzp6el2R/EnJycDAC5evAgAmDt3LmbPnt2gztSpUzFjxgynSw7YbDbk5eV5MXLhF/rKlSsAaAZCoKB7Ih7mwjmgqqrhQU0Eyi8XhfQ9qTDZUFVVbbdMqpcgz1ju54ico78R72jzf//n8Tkqlcru+xDdE8dsNhtkMtepj0fJUd++fZsckL/o9XoAcLhWUe1xnU4HAIiOjkZ0dMOF1aRSKaKiohAXF+fwcWQyGdLS0rwRcp3aLD8tLY1+oQME3RPx1ByUg7/ub1OWnAbV1b+7UL0nRRfLER1t/5/Qtm3c36rJX+hvxLe2Oilz9B5E98QxdxIjwMPkaNmyZU0KJpDUNlc7awFzd50jX/zSSSSSug8SGOieiIOxWRvNXZdqIxrcj1C8J2dLDQ5fnwL1+Yby/fALJ9Px73BymrOfN92T5gm5n1pkZCSAay1I16s9XluPEBJ4eI4Dz9oaF8jk/g/GzwwWO8+bhDal0mGRozmbJ06c8E0sBEATkiOdTgeTyf2Va6+3evVqLFq0qMnnu6LRaJCYmIiioiKwLNuovHbRyDZt2vgsBkJIM9lLjAAwYZAcEeIOexOOiPd4nBz17dsXzz77rMPyHTt24OjRow7Lf/rpJ68sJulM3759YbFYcOzYsUZltWsjib1mEyHEMd5qtXuckTd7O8iA5myNozbxtEBtuPlNoRA7hLDlcXJUuxWII48//jg+/fTTZgXVXFOmTAEALF68uMHx6upqrFixAjExMRg5cqQYoRFC3MDb7CdHod6tVmV08LwBRKtC+7mHrYICh0UDaWsQ0fjk3zBfrPBaUFDQoCWovFyYzrpnz566WWUpKSnIzs7GwIEDMWnSJKxZswazZs3CyJEjYTAY8NVXX6G0tBTvvvsujTkiJIDZCi/ZPR7q3WoW1vFrZ4QytFvNwlbr1h6fsnWrszlsxBuC5q/t999/t7uK9SuvvFL39YQJE/DGG28AAP7xj3+gc+fOWLlyJebPnw+5XI7u3btj4cKFQbEkASHhijOZYPrjV7tloZ4cSZwsI9c+gbrViODWW28VO4SQFzTJ0cSJEzFx4kS360skEkyfPh3Tp0/3YVSEEG+z5p5xWMbIQzs5ckYlpw1Gw0l2RARwdT0+4n8hN5WfEBLcHLUaAQDcXMCNkKDwyCMOi/6ixEhUlBwRQoJGqHereX+0JgloH3/s8Slr1671QSDkepQcEUICBltW4rSc0YT2uBuWo/SIODdhwgSxQwgLlBwRQgKG6fDvTsslSpWfIhGHlaVd1MOGg7W8AMD+zqDEnyg5IoQEBN5qga3I8ZovyqzufoxGHBUG+2+Y7RIi/BwJ8TknCzwa/RgGsa9JoxtLSkqwY8eOJpWXlDhvNieEhB+eZaHfudFpHUXnHv4JRkS/5ZbZPd4tOdrPkZBAtHGj878R4j1NSo6OHj2KJ554wm4ZwzBOywkh5HrslctgS4sdlmuHjYFEFdpdas4Wz41S0Sy9kGJzvLmw6fXXgRdesFs2atQoX0VErtOkv7jmroDNME5WOiOEhB3L+dNOy6VJyX6KRDzOtg6hNY5CjJONz9UOEiPiXx4nR6dOnfJFHISQMMbpaxyWyVqmgJGE/vBInYUVOwTiL072UyOBIfRfcQghAc/ZLDR5egc/RiIembO9Q0jYy83NFTuEsELJESFEdIxC6bgsTLYModEGYeKnnxwWnfjxR4dl6enpvoiGOEDJESFEdDznZH0fueMpz4QEnWHDHBZ1uf12PwZCnKHkiBAiPpvjwcihvmUIIa6kpqaKHULYoeSIECI63klyJImI9GMk4jFaaUB2yDOZHBbx8+Y5LLtw4YIvoiFOUHJECBGdo5WxZS1TIAnx/dRq5ZYZxA6B+Jpa7bBoxL59DsukUlrKwd8oOSKEiIoz6B2WaYbc6sdIxHX6iuPlDEjo27Z9u9ghkHooOSKEiMp64Yzd45LoWKez2AgJKs4mHTgZUzR48GAfBENcoeSIECIqW1G+3eMMzVIjocTJRrPF+/c7LPv55599EQ1xgZIjQoiouOpKu8dliS39G4iIWM7xlkxaBe2rFhJYxwPuW7Zq5bCMttsSByVHhBDR8BwHzsEMHllaWz9HI57LVUaHZRktIvwYCfGJZu5HSvyPkiNCiGh4ixngGv9HLYmKgSwhSYSIxHHwUqXDsk5J4bGUQUgbMsRhEe9kW5CfnKymTXyLkiNCiGgsp4/bPa7I7OrnSMRVpjc7LItU0SKYQe+XXxwW3TlnjsOym2++2QfBEHdQckQIEQXP8zAf+8NumUTleD2YUES9LiHMyVgjAFi7dq2fAiGeoOSIECIKrqrCYRmjVPkxEvE5So7axIXHApghzdl+aVeuOCxq5WSQNvE9So4IIaLgrRaHZWGXHMF+dkQTlULApk0Oi345fdphWX6+/SUuiH9QckQIEYeT7gYmzLrVHJFQdhTcDM63hLnhhhsclkkk9PYsJvrpE0JE4WyzWVoZW0C5UZDTOukWraz0WxjEc5QcEUJEwVsdJEcSKRj6rxkAwICyo1BV5aQsLy/Pb3EQ++gViBAiCt5of8NZVY9+fo6EEB9wMUstJibGYVmqk73WiH9QckQI8Tu2qhKmQ7/bLZOnt/NzNIFLIaOX6KAlc7LtS0mJ/+IgTUJ/eYQQvzP/ZX99IwBgZLToYa0uraLEDoH4gCXK8X2dPn26HyMhjlByRAjxK57nYc0967hCmCVHBovNYVmMOrx+FiFjyxanxUql4wkHn3/+ubejIU1AyREhxK94g/2xRgDASGVhNxj7aIGzobkkKI0c6bhMp3N6qsxZdxzxm/B6FSKEiM566YLDMkeLIYayo4WUHIUTXqMROwTiBkqOCCF+Zb1wxnEhx/kvEEJ84dlnnRZnZGQ4LOPo9z9gUPsdIcRvOLMJbLnjmTry9PZ+jIYQH3jzTcdlHIezTrqNGVr1M2BQyxEhxG/YogKHZYxSBVW3vn6MhhAvc7Sw6VUWF+UkcFDLESHE53ieh/nQ7zCfPGK/glQK7YjxkETS1HUSxBQKx2UbNzqdpUZdaoGFWo4IIT7Hll1xnBgBiBg1CdLIaD9GRIifjRrltJi61AILJUeEEJ+z5p5zWi7RRvopkuCRGkuzmoLKqVNOi81ms8MyrbMNaokoKDkihPicJeeYwzJGpQETpmu71Jgcj0FJinTcBUMCUFaW47KdO6FSqRwWV1XRcg6BhpIjQojPOdsSRJ6S5sdIAsu+3HKHZVlJNP4qaBQXOy8fOtRpsVQq9WIwxBsoOSKE+J6T5EjZo78fAwkspTqLwzK1gt4wg0bLlo7LbrgBly5dclJ8gw8CIs0Vnm3ZhBC/4XkevMlgt0zeug0kTrobQl/4rQgedn7+GWlOBlvv2bPHj8EQd1HLESHEp2z5Fx2WKTp392MkgYdzkBslRNB4o6DhYhaaKzRLLTBRckQI8SnrxbN2j8vT2kOW6KQ7IgywDrIjmYTeMIPG5s2Oy3Q6fPvttw6L33rrLR8ERLyBkiNCiE+xlfYHHasH3OTnSAKP2cbaPS6h1oTgcPy483KtFnfffbfD4jlz5ng5IOItlBwRQnyKNxkbHZO1SgUjdzxIO9zRNP4g0bWr47J//APl5Y5nI5LARgOyCSE+w3MceEvjxe9omxCg0uB4plqnJFoUM+CVON5AGQDw4ouId9ICuH//fi8HRLyJWo4IIT7DW8wA33hcDaMM5xlqAp3FfpcaAKjlNI0/4LVo0azT+/alTZYDGSVHhBCf4XTVdo8zKrWfIwk8zkYVSWhAdmA7dMh5ucmEXbt2+SUU4huUHBFCfIJnbdBv+c5umYRajmC0Om45ogHZAa5XL+flSiWGOlkV22hsPA6PBBZKjgghPmE5c9JhGbUcAaeKa8QOgTRFjYv7ZrG4HIjtbJ81EhgoOSKE+IT1wmmHZTTmCCisotaDoBTlYjKBXI74+HiHxfPnz/dyQMQXKDkihHgdz3FgqyrsF0qkkGgj/BsQId5w5Ijz8i1bXF5iwYIF3omF+BQlR4QQr+NqqgDW/pgaeXp7MHKFnyMixAt69HBefttt+PTTT/0SCvEtSo4IIV5nOviL3eOS6Fio+9/o52iCS7SaFscMSJcvOy9/7jkAwMMPP+ywCg3EDh6UHBFCvIq3mGErtv9GorlpBBgpreHjTGqsRuwQiD3Jyc7L//EPnD9/3mkVGogdPCg5IoR4FVtZAfBco+OMQglJBK2MDQAmJ9P4s2h17MBz2vHkAgDAnj2AXI727ds7rPLHH394OSjiS5QcEUK8ijPo7B5XdMgCQ+v3AADOleodlsVQt1pg4XkgM9N5nSFDcOLECadVerlaG4kEFEqOCCFeZcuz37Wg6JDl50gC1y/nSx2WUQIZYJ54wnn50qUAgC5duvghGOIvlBwRQrzG9OdvsF660LiAYcBotP4PiJDm0OuBDz5wXmf6dJcDrW02mxeDIv5AyREhxCvY8lKYT9pfB4ZRaWggNgk+ES7W4/r+e4BhoNE4H0Qvpd/9oEPJESHEKyxnHI+5oEUf3ROrofWfAobZ7LrOuHGwWq1Oq7AO1vsigY2SI0JIs7GV5bCcdbyXGiVH7hnSPkHsEEgtV9PuS0oAAAqF84RWIqG32WBEd40Q0mymA/YXfawlb9vRT5EENxqLHSAuXnRdJyHB5QazrspJ4KLkiBDSLLbSYtiuFDosV2b3hiw5zY8RBTYb13gNKBJArFagTRvnddauBQCnG8wCQGxsrJeCIv4mEzsAQkhws+aetXuckcoQMWEaJEpaFbi+y1UmsUMgzrha0wgAJkwAz/NOqxw9etRLARExUMsRIaRZ2HL7a/Zox9xJiZEd+ZWOp30zoH41UdXUABfsLEVR39WWP1djibKzs70VFREBtRwRQjzGcxzY4kLwHAfOTnIkS06DNDJahMgC3/HL1WKHQOzheSDKxfY2Y8YADAODweC02ldffeXFwIgYKDkihHiEt5ih37EBbHmJwzrSmDg/RhQ6ZBJqORKNO7PKfvwRAKDVOl/QdNq0ad6IiIiIutUIIR4x5/zlNDECAAklRw5xTsaqxGpoXzVRVFW5rpOXBzAMampqnFbbvHmzl4IiYqLkiBDiEcsZx+sZ1ZJG0ywde8w2xwsCDmwbT/uqicFiAWJiXNdLTQXHcYhy0fU2YsQI78RFREXJESHEbTzHgTc63lEeABi5AhJKjuwqqHQ8Uy2jBS2UKQql0nWd338HAAwYMMBptVWrVnkjIhIAKDkihLiNq3E9mFiWnEr7qDlQUOV4ppqEWo387513XNfZswfo1w8cx+HAgQNOq955551eCoyIjZIjQojbuMoyl3XkqW39EElwcrbGkZQGY/tXfj4wZ47rekOGAHC9eez+/fu9ERUJEJQcEULcwttsMPyyw2kdSUQUZK3b+CegIFRptDgso5YjP+J5IDXVvXoAODdWNe/bt29zoyIBhJIjQohLPM9Dv2N93ZuFI9pbb6cutSaQ0eak/uXOz/uzz+q+dNVqZLPZmhsRCTD0F0kIcYm9chlsabHTOpobR0CioUHFTZEepxE7hPDhbiL6wAMAgGeffdZlVVfJEwk+lBwRQlyylRQ5LVd17wdZ63Q/RRN6OtJMNf9Yu9Zl6yeAujrHjx/Hm2++6aKqG9cjQYdWyCaEuMSWXnFYpup7A5QZXfwYTehJiVaLHULo27MHmDTJdb0TJ+q+7Nq1q9OqL7/8cnOjIgGKWo4IIU5Z8y7AVnDRYbk0Jt6P0QQvlnPcwkAz1XyM44CbbnJd79FHgawsAMA7bkzzf+WVV5obGQlQlBwRQhziTCYYf9vtuIJECmksbRXiDitrf8YTDcb2MZYF3B0T9L//AQDMZjPmuJjmX1TkvKuZBDf6qySEOGT+6w/wVrPDckX7TDByhR8jCl4Xy+3v5N69dbSfIwkzMjdHj3AccHU5BZVK5bJ6UlJSc6IiAY7GHBFCGuFZFuajB2HJ+ctpPVXvQX6KKPidLtHZPR6lpM1mfYLn3Z+ZtnZtXWL0f//3f25cmgZhhzpKjgghdXieB28xw/znb7Ccz3FaV9WjP61p5IHLDrYOidVQcuQT7iZGr78OTJgAACgvL8eKFSucVl+/fn1zIyNBgJIjQggAwHY5H4Z9u1xuLFuLcaPrgQhMVtZhWbyWuiW96r33gNmz3av79dfA3XcDAHQ6HeLjXU8uGDNmTHOiI0GCkiNCwhhvtcJWmAfObILpwF6PzpXGt/BRVKHnjIMuNQBgaNsQ71m2zP3ECKhLjEwmEyIjI11Wd2cbERIaKDkiJExxBj30238EV1Pl8bmyFsmQxtAsNXcdK6wWO4TQd+QIcO+97tevN25IrXa9ztSGDRsokQ0jlBwREqbMxw95lBgpOnUDV1UOaUw8lNm9fRhZ6NGZrWKHENq+/bauFcgt1mv3Y926dW6dMnr0aE+jIkGMkiNCwpQt3/HCjteTJraEuvdAH0YTuioNFrFDCG1Wq2eJ0aZNddP7jx49iokTJ7o8hbrTwg+tc0RIGOJ0NeAMjsfBXE8aS6tgN9X+ixUOy5Qymu3XLKtXAwoPBrRXVgIjRwIAbDYbunfv7vKU/Px86k4LQ5QcERKGXG0kez1pXKKPIgl9ueWOZ/8NaZ/gx0hCzGOPAZMnu18/NxeIvrbgplzuegmFpUuXIiUlpQnBkWBHyREhYYi9ctntuoxMDllKug+jCV2uFgtMi6MNZ5vk0UeBDz90v35VFZAu/A5XVla63RJ0zz33NCU6EgJozBEhYYYtL4Xl7Em36jJyBdQDboKE1jRqEmebzQK0r5rHSkuBRA9bMW22BnurxcbGunUarYId3ig5IiRM8DwP04G9sJw54VZ9zc0jIWvZmlbBbgaW3mC9w2YDfvwRcGPwdAPvvdcgMerUqZNbp126dMmzxyEhh5IjQsKE9VyO24kRo9ZSYuQF+ZX2twwhHsjPB1JTPT9v3DjgyScBAHq9Hu+88w5ycpxviQMAX375JVq3bu3545GQQskRIWGAMxlhOvSbe5UZCVQ9+lFi5AU7cq44LLuxAw1yd0mna1pitH8/0LcvAKCkpAQtWri3mvuqVatw5513ev54JORQckRIGDAfOQDeYnZYrrnhFkAmA2/QQ5qQRFP3/aBdvFbsEALX9u3Arbc27dzLl4GWLeu+dTcxAkCJEalDyREhIY6troTlnOPuBHmbDpCltaO1XLzM2YBejUIGhYwGY9v18cfAI4807dzrfuajRo1y+1Sbzda0xyQhiZIjQkIUZzbBdjkfxl92OKyjuWkEZCnplBj5gIV1vKryjR1ofaNGli3zbG+0+tRqoQuuHk9+p6urqyGlbmRSDyVHhIQgtrwU+p82gjc5GBAskSLy9qmQRLjeiZw0TW6ZwWFZ6xha36jOX38B2dlNP//SJaDeAGqe5yHxYImEyspKREbS3wFpiJIjQoIQz/OwnDoG64XTYJQqKLK6Q54sDFzlLWbot/8I3up4Ty9Fh06UGPnYgTz724Z0SoqChFrqBFu21G3n0SQHDjRIjIxGIzQajdun19TUICIioumPT0IWJUeEBCFLzl8w/bmv7ntb8WVobhgOW0kxLKeOOj2XUamh7NrL1yGGPYPF/hiWtNgwbzUyGoExY4CffmredQoLgVat6r61WCweJUbffPMNJUbEIUqOCAkyvM0G89E/rjvIwfDzNpfnMkoVtDePgkTt/psI8S6ZNEwHYufkANOnAwcPNu86paVAfMPZlBcuXEC7du3cvsRvv/2G/v37Ny8OEtIoOSIkyFjzc8FbHU/Ld0QalwD1oGGQRru3fQJpOmebzcqlYdil9tNPwLBhzb+OxQJct2Hsiy++iNdff93tS7As69GYJBKeKDkiJMhYz57y+BztiPGQJST5IBpiz9aTxQ7L4rUKP0Yikt9/B/7zH+DiReDXX5t/vZtuAnbtanAoNzcXbdu29egyxcXFlBgRt1ByREgAq10rp3Zasq2kCLbiArfPZ+QKaIbcSomRH5msrNPykN5s1mgEbrnFOwlRrV27hOSonj/++AN9+vTx6DImkwlKpdJ7cZGQRskRIQHKkntWWNnaaIAsJQ3y9A4w7nNvEKs0IQmKzK6Qp6SBkYdBS0UA2X22xGFZdnK0HyPxs+pqwM0d790ilwNmM3DdzL5PPvkEs2bN8uhSZrMZCgX9HRD3UXJESADgeR6WnL9gzT0LRi6HNC4B5pNH61b8teadhzXvvMvrSCKioOzaE/J2mbSwowh4nsfFcsfrG/VsHeO/YHztwAHg22/BbN6MNidPevfat94KbNjQIDH65ZdfcMMNN3h8KRpjRJqCkiNCAoD5+CGYjxyo+95W5H7XmTQ2HvL2nSCNioE0KRkMvRGIptxgdVqukgfxKswcB9T+bt19N/DttwAAr6bgiYnCWKUJE+oGXhuNRvTs2RM5OY63wHGEEiPSVPRbQ4jI2JoqmI/92aRzJVEx0N5yO5SZXSFr1ZoSI5FtOlHksGxk55YOywLamTPAjTcCWq2w4CLD1CVGXnP0aivplSvAXXcBKhUAICcnBxqNxuPE6LnnnqPEiDQLtRwR4idseSmMB38BV1MFRq6AJDIa0uhYWM6eBDjng3jtkbduA/XAm8EoaJBpIDBZWYcLPwJAqyiVH6NpJrMZKCgA2rdveLzA/RZNt4wYAXz2WYNVrgGhe3LkyJHYunWrx5f8+uuvcffdd3srQhKmKDkixA/YqsoGW3rwJiO4mirYCvOadD1ZSjrUQ26llqIAkl/pYB+7q+SBvPhjTg6wYAGwfLl/Hq9VK2Ga/3VrFgHCAo0DBw5s0mWtVitkMnpbI80XwH+thIQGnudh2v+z073OPCGNjYdm8HBKjAJMYZXJYVn3lBj/BeKKqV6cNTXA3LlAp07+S4x27RK2/rCTGHXv3r3JiZHNZqPEiHgNvboS4gU8a4PlXA7Mfx0CW17aoMx6/jRsVwqbfG15altIouPAqLVQtMuEZvhYMHbeWIi4lDLHL6d90wNgVfJTp4A+fQC1Whg3xDBAVBTwz3/65/EnTQIuX260ZpHZbAbDMGAYBkePOt8X0J5vvvkGHMdBKg3iwe4k4FCaTYgDbGU5LJcuABfOw3DuOGAygjPqwcgVULTvBEVWNhipDGxNFQy7NoOrrhROPLIfisyuUPXsD95qg+nQbx49rrxtBsDaAJ6HLK0dFG06eP/JEa87UlDpsEwi1rIKa9cC//qXsGK1WKqrgcjIRod5nsddd92FFStWNPnStLAj8RVKjgi5DltVCfPRA7DmnQfP82CqqmCLjq5bN4g3m2A6sh+WC6ehyOgC89GD4C0N9zqz5PwF9spl8CwL3uygu4WRADzX4JA0sSXUA26iLrMg8/O5Uodlo305S620VGgR6tABePNNYe2hPXuA9HRhTI8Y2rcHpk4F7r8f6NixQZHZbMYjjzyCL774otkPU1NTQ4kR8RlKjgi5itPrYD72Byznc+oWX3Rav7oSpoO/OCxnK8oclmluGgFZchrY4kJY83PB1VRDGp8IZVY3SoyCDM/zOFlU7bC8ZVNmqVmtwvT2du2EmWO//w506SKsQP33vwOLFzs/X6zEaPJkYeySnd/ho0ePonv37s1+iDFjxuDHH3+kRU6JT1FyRMIez3HCIozHDwGs51PqPSVPbQt56zYAAFmr1pC1au38BBLQcq7onJbL3JmlZjBcGwv0xx/C2KAgYunUCbJff4XEzvYhPM9jxIgR2LZtW7MeY8SIEVi9ejUiIiKadR1C3EHJEQkJPMcJs8EsFvDWqx9Xv2ZkMkgTkiDRaBudx5mMMP6yw6MVqZuDkSug6jPIL49FfK+4xoQ9TvZS0yjqvcSyrLDIYVwcUNsdVF4OZGcLs7eCyY8/AmPHApWV4KxWFOr1SIu+tm/c0aNH0b9/f5hMjmfweSI/Px8pKSleuRYh7qDkiAQlnufBll2B7VIurAUXhcHQzrrCGAayVqlQtO8EWUoaGKkUtitFMOzdDt6od/5gjAQSbSQkGi14qwVcVYXjqnKFwyn7Ek0E1IOGQqKh/3xDQX6lERuPX7ZfyPHAoUO4ofwskNcReOABoKqqYZ1nnwX+/W/fB9pcrVoB//2vkMRxnNDVVztbMiZGOKbXw2q14r777mvWAOvrrV27FhMmTPDa9QhxFyVHIYjnefAWM3iz6dpns7nh12YTeMvVcpMJPMeCUaoh0Wgh0UaA0Wgh0URAotGC0UQIx0ReQ4RnWbBXLsN6KRfW/FzXSU2Dk3nYCvNgK8wDo1JD1jIF1ovnGw2IbkAihaJDFvioeER0zIBEIgHPcbCcPg7zkQPgbQ330VJ27gFl154wHfxVGLd0lTQuAYpO3SBPaweGphv7VkWFMDPK1e9qcTGwdCmg1wM9egB33NFo93e7/voLrNGIkvZZ2HjyinDMaAR27waKrwBpqUCvXsC776J1WSFSD+10/DsWDIlR//7A5s1CEuTAt99+i+nTp3v1YV944QW8+uqrtP0HEQ0lRwCGDx+O4uLiBgP8kpOTsWXLFhGjEpIc2Gz1khohmeHM9ZIai7lhudkktFy4MaC40eOZjOCqyh2WMwrl1cQp4lripK1NoLTCZ2nzfqWE52wFZzSANxnBGw3gjEawZcWwFVwCbzW7voirxzAZYc0967gCw0DRNgPK7N6ARgvkXVvFmpFIoOyUDXlqW5iOHhSSLaUaquxekKcLWy2oB94MRZceYMtKII2OhSQ2PnwGj27ZAvz0E5CcDDz0kDCOxpf0euDdd4U38L17hWPR0cCiRcCjj16rt28f8MsvwqDm338HXnml4XVuvx34/vtrCVJenpBgbd4stPpcZZSrsLH3LSiLjBMO9OwF5F0Eyq4Ovs+9IMwWAzDy0A5ImvB3GBBkMuCvv4DMzAaHjUYjduzYAY7jcPjwYcyfP9/rD63T6aDVNu4CJ8SfKDkCUFVVhffeew/Dhg1z+xzeZIT+p01ejYPneTDFxdAd0wJWC3iLyS8DhN3FW8xgLWbAySwsRqEEpDJhxpVUCkYiFWau1PuakUoBiVQ4xjBCwnc1EeJNRvCs4/2pfE2iiYB6yC2QJSQBADjO/n/9Em0ENANvvnbgwgVg9mxhkbuRIyG9/35Io2J8H7Cv1e7Efv488OuvQnLRs6f9um+9Bfztb9e+X7ZMSEqa2uJYUABERAjJTi2zGSgvB2M2C7HZG5xbVQU89hjQrx/Qu7ewyOHcuc4f68cfhX2+XAwaPtKmy7XECAAO2d8wePjR3cGVGD33HDBtGtCtW6MinuexYMECLFy40KchGAwGqH2dTBPiprBPjliWhU6nQ6ydWRbO8DZbk/fFcnhNngeqqsCx0UHb0iCs92NGEL0t1JG1SoV60DBIVB5Ovb58Wdi1PD9f+H7VKmDHDqHbxhvdaDwPlJQInxMT7U6T9rrTp4HHHxe2erBdl6wuWABc32JgNAL/+EfDYwcPAl9/Ddx3X8PjublCK07//sKaOMXFwtT11FShXKcD7rxTaIWqdeONQtJz5AgkAFJatIDkyhXnz6FPH+EN/+uv3XvObsymOtqmi1uXalcs0lR6d/35p5AISSQOuxNZlsXvv/+OwYMH+zQUGmxNAlFQJkdr167Fa6+9Bp1Ohx07dqB168ZToa1WK5YtW4bvvvsOFy9ehEwmQ1ZWFmbMmIHhw4fX1auqqgLP8/jyyy/x1FNPwWQyoUePHnjhhReQnp7uz6dFmoiRySFLToUkOhaMQglGLgcjFz5DLgdbVAjLuVPgdA7WomEYKLv1gbJLz4ZJ6alTwKZN0DIMMH06kJBg//znnruWGNX65huhxWTJkqYnSHv3Ct0/+/YJ3UeA0GqzcqWw6J8rRqPwht+yJdC9O3DmjNDS8vPPwnO77TZhKwebDdi6Fbh0SUjC8vOBjz92fN0FC4Q31tqBsjwvJIOVlY3rzp0rdLP98IMwHsiVQYOEFqrrXe2qqiVzlRjVcjcxcoNe6V6rRv/TBxGQ/9r83/8BTz0ltKbZac3jeR4Mw6CmpgZRUVE+D2fOnDn497//TeOKSEAKquSovLwc8+bNw44dO1w2vz799NPYtm0bhg4divvvvx9msxkrV67EY489hgULFuCuu+4CIKzY2qVLF2RlZWH+/PlgWRavvPIKZsyYgQ0bNoRmMy8jAaNUComEUiV8lsrAGfXg9TpwRoPzgcpiYzkwNhvkSi1kvASyKguY3D+F/4ClQncdZLK6r2UaDRRt2oBtHwlLVSmsly7UdVcyKg00g4Y2XGuouBh4/nng888hAZAIgH/+eaEV4vHHG3Yr/for8NVX9uNculSI49NP3WvtWblSeDOPiQFGjhTG7OivG3R+6BBw113A/v3C8718GfjPf4CdO4UEaMECICUFOHdOaJkpc9wFitdfB95/X1i07xfHi1naNXGikLBotcKb7pkz9usVFQFffun+de0lRgHi6xsnu1Wv28UTPo7E3oN2E8Z3/f678Dv69tvCYpGffSZ0P95/P/Deew0S9W3btuG1117DwYMHob/+98yHxo0bh8WLFyPB0T8bhAQAhueDp2N86NChsFqtWLRoET755BPs37/fbsvRtm3b8MQTT2Ds2LF466236o4bjUbccccdKCkpwfbt2xEfH2/3caqrqzFgwAC88847GDFiRKPy4cOHgyvMx3dR3t38kwfAAZBoNWBUKkClBtTOPzORUWAiI4UkR1kv2VGqwCgaHpNcrQO5wmm3Hc9xwuBsgw68QQ9OrwNn0IM3CJ85vQ68yeD+oG+rDdDVCDuAX/+h0wkv3rXx1G6IWft17WeOg6RaB9mVcsgLr0BqNDftv/PISHAdOoBtnw6kpUHWpTuYjAxhm4PYWCFZWLhQ2A/KkUGDhNaisWOFBOSPP5w/5sMPAx9+2DBB+vln4J13hPVuZs0Czp4VkiF3rV8vjL+ZNUvYQqJWVhbw22/AqFEBnWgEk0pNFFYNugO8G13dN/+1FxmXz/s2oFdfFVoyn30WOHwYGDhQ+F2y83qWe+YMftm7Fxldu6J169aYNWsWfvzxR9/GZ8fJkyfRqVMnvz9uuOI4Dnl5eUhLS6OWuevU9hzt2LHDab2gajnq0aMH5s2bh7i4OHzyyScO661btw4AMHPmzAbH1Wo1pk6din/961/YtGmTw+mnUVFRiImJQUmJ48XdGJaFpNiNbgIPSXlAwnJgOA4Me/Wj/tcsB4bjr30NCK0TMTHCm7s7H/XrRkU1GnPASCRgVGpIwACMDJAqAbkGUEQASiOgMoLX68FXV4HTVYOrqQavqxGSqaoqoKIcKC8HX1kBVFaCNxiuJj0Q3mAYBjyDa4kQx0FiE56PxMZee242VvhZ1H7mvJDH19RAcugQJIcONS6TyRqPr7Hn11+B8eOBtm2FgdiufPKJcO3//ld4vm++Kbyx1Vq+3O3w64wfbz/WkyeFaemUGHnFmgFjGw7AdqFdca5vApkzR9iao1cvQKEQjq1a1agaz/PYtm0b3njjDfz000++icVDtDksCUZBlRy98847btU7dOgQVCoVsrKyGpX1vNol8ueff2L69Ok4efIkVq9ejRdffLEuwy4vL0d5eTlSaweI2sFYWURedJw8+ZXNJrQe1G9BcBMvkQjJUmQkYLEI41SMRmE2kBPM1Y+Q+p/EncSoPncSo1offAB+3z4gPR3Md9959jj2OIt1167mXz/MVasjsPyGiR6d0+fsIcgczG6sxb33HqDXgyksBPP++3br8HPmgJ88Gdi8GUxlJfhJk4D6g6I5DmazuUHC8c033+Cee+7xKF5f6tChA5YsWVI3mNvRrE/iGxzH1X2Qpgmq5MgdOp0O5eXlSE9Pt9ucmJycDAC4eHVjxri4OKxduxZqtRpPPPEE9Ho95s2bh44dO/p8lkYgYDhO2MKg3PH6RsR9tpYtIS0vB2NpvEo2c+iQMGaIBLQD7XvgULvGU9pd6XXhmMMyTq1G0fLlsNSfKj9nTt2X0pISSCorYe3QAVabDSaTCZH333+t7tW1tvbs2YP7rp/9FwAefPBBzJ07t24wd+1rb169NcKI/3AchytXJy1Qt1pDNpsNMjeWFwm55Kh2YKFGo7FbXntcpxM2i0xKSsLixYvxzjvvYPDgwVCr1ejXrx8+++wzt36AJLRxU6eiRqFA1Lp1YHTONxgFAMnHH4OXSoGJE+0mSCRw8G3agMnNrftep9LiXFK63cSInzETSG0NVFeDeffdRuVjspPBvfqb0NV6/ryw7ADPg1m9WugGmzEDLXv3dhxMWhp4nsewYcOw57qZeYFs9uzZeOaZZ2C1Wml8SwCpbTGie9KYu+/rYffuXzv+vP6A5F69emHZsmWeXUirFVbV9SKe52GqqIDKZAJTWSlMfa6sFAYsh6KoKGHfJpVKGNzt7INhhKn0LVsK57Rs2fDrFi2EOizb8KO2y/HsWWFGVf3PzmbodO8uDM4ePBgVeXmIfO89MMuWCQsKFjjYpHbUKEhuv12IY/VqYNIkYf2e5urRQ/hde/XV5l/r//0/4WfywQf2y9u0AebNE9Yo6tBBGOdS2728Z4+w0OXhw54/7qFDQsJw4QIwY0bj8shIYTbVm28Cx49fO/7SS8DMmcIg+U2bAJUKvEYDnDkDxt1u0OJiYTD8uHHCsggAEBMDZu1acO3bw7T4C/xZUIUTnXoDWVcHDfM8cPqM8LNKTwNzdbHJC+UVKGjXHgMvXICU5yC32dBLyyHllRevLUjZogUA4R+wn3kecrkcN2RnQ3XdmxTHcZg2bRqWN2XMmUjy8vLAcRySk5Mhv7q/Wu3gX4lEQm/EAaT2ftA9aZqQS44iIyMBwOHU1NrjtfWaLClJWLvFi3iOQ/HVGQZM7S80xwmzuioqHH/UJlL2PsRYYVutFqaTJyc7/mjVyv7qxr5y440Nv+d54U3z+oSJZYU30enThWnPtYlpVBTw5JPAgw8KicWiRQ2nycfHC1Pqa5Pu228XpuZPnuz5WKb6EhKENZPatBGSg4MHr5X16QM8/bQw027IEGH2miPDhwsz2J5+WnjuAwcKaxBdn+g984yQjFw3mQGA8DOs3y1oMAjrG2k0wvU0GuH7W25peN7DDwsJXu01Fi0SFpms77PPgClThAUjz58Xfuc7dRISZwD4/PO6qryjmTiPP9446duypS5Zwe7dMG7YhANXjNgkicRfK3+GQvE7amrMqKw0IfH0aQxNSYFKpcLp06eh0+nQvn17JF79PT19+gyWL18Onufwh8mAmw9sxOPz5yHryScBrRbV1dXQ6/VISEhAYWEhbr75ZuTWa5kKZl9//TXuuuuuoF2clhBPhVxypNFokJiYiKKiIrAsC+l1C/DlX12sr02bNiJE1wQSibB9QnS08AbpCZ4Xpsq7k0Tp9cK0crX62odK1fB7d8oiIuzOgAs4DHOt9WnIEPfPU6uFBOKhh4Tp+Tt3CongnDlCklLf+PHAunVCS8n1g+WVSmEV6G++abgkglQqrF6cny8kCCNHXts+Y9cu4XqlpcCAAcIyArU/5++/B15+WSizWoXr2GzArbcCb7xxLUGoNX26kIwsXCi0kPG8kFw88oj7PwuNpnHr6bBhwvpJ8+cLcYwfD7z22rVyhhHK775bmAAACEsjTJp0rbx9+0YPpbfYcLSgChaWQ1qMGiYrB5bjIZEAVpaD2cZB+eJLkF+8CGzYAE4qhXX+AvA3DUVppREMgBaRSnwZ3xnLt63AmTOnGz1GQX4+cnJyEB0djaKiy/VuiRQymQzmq8knz3E4fWgHTpgM+OD554U1sULUq6++ipdeeknsMAjxu5BLjgCgb9++2LhxI44dO4Yetf+xXrV//34AQL9+/USIzM8YRuiuiIwE0tLEjia0REUJi+z9/e/O640dCxQWAidOAFeuCCtQA0JLS9u2QjIxa5aQyEilQgtKt25297iCViskNfaMGCF8eEKhELb8uH7bj+ZgGCFZ+NvfAJNJ+N273qRJQtfcjh1AXBwwdWqDxQkrjVYcuFiOapMNiRFKdEuJxndHC2CxCa14p4qqUVVVhYRSBplJkThXoofJxkKjkOHGpStwOf8KDpdZcCE/H7olG9G6dWvExcWC54Fvly/HWUcLVgIwGg0wGg0NjrEsC/ZqC2zJb+th1VU0aWPnYDJz5ky88cYbSExMFDsUQkQRksnRlClTsHHjRixevBjv15suW11djRUrViAmJgYjR44UMUISVuRyYQyTPTNnAqNHA8eOCStvN3PV4OunePsTz/O4UGbApUoDIhQyZCdHQ3FdHY7ncbnKBFnrtmjx+OPQmW04eK4cZ0pqAAAqmRQm27Wu4DK9GaeKhQU5WZbDmTNnwLIsWrRIRI3eiM3785CYmAiVSoX8oitY9PtBqFQqHD12FKdOnqy7Tt++/RAZGek0MXLFUl0Ga01ozeqMiIjAnDlz0Lt3bzAMgzZt2qBr167UfUbCXtAkRwUFBTh27NpU2fKrU8/37NmDuDhhkbaUlBRkZ2dj4MCBmDRpEtasWYNZs2Zh5MiRMBgM+Oqrr1BaWop33323+WOOCGmGNWvW4LPPPoPNZsO0adNw7733gmEYFF2+DLlc7vHWCl9++SWeeeYZlJeXY+zYsfj8888drgDfXDZO6NJSSCW4dOkSZs+ejYMHDyJj8EjcPOUBJCUJXXg5V2owsXsKairL8ebb7+BsuRktsgejR4/uYBgGcRoFLpdWYP2mLbhw4Tz0ej0GDR6MAf0HIDLy2ng0q9UKo9GIb7/9FkVFRU2K+cCB/c1+3vq8U82+hpi+/vprREdHo7KyEnq9Hl26dMHAgQNpwC4hdgTN9iFr167F8y769idMmIA33ngDgDCD4ptvvsHKlSuRm5sLuVyO7t2749FHH0Xfvn2bFYu7y497ipZ8t0+v1yM/Px9t2rTxa6tIZWUldu/ejaKiInTu3BlJSUlITExETExMo/+sbSyH/Eoj9BYWKdEqxGiubzO5ZtWqVZg6dSoc/elJpVI88sgjeP/99936D/7LL7/E/fXXxAFw++23o0+fPvj6668RGxuLBQsWoH379li9ejVMJhOysrKwceNG7N+/H5GRkXjppZdwxx13AABqTFbs2X8IF0+fwKbv16C6uhqdO3dG586dcaHCDGV6Ntq074j8/Et4b8Fz4G0WaFpnQKpQITo6Bo8++kjdfbLZbPjg9Xm4UmOCPCIGAJCRkQlthDCA+dzZsw5/BkNuvBG7AmCVZ1PJJRguX4DJV6tfe5FarUZ2djb69OmDoUOHYuzYsVAqlT5tCaLXrcBD98Qxd9+/gyY5CiSUHPne6dOnsWHDBmzYsAF79uyB1WqFVqvF1KlTMXPmTAwaNMinL/irV6/GQw89hEo7O80nJibirrvuwvPPP4+4hBY4WVSNo4XVMFqFWWkMAP2ZAzj3x16YzWZwHAepVIq+ffti4sSJ6Ny5M4qLi13GsGjRIjAMg9LSUkyYMAGDBg1qVOfQoUMYNGgQTCaTw+vIImIhUahgrbwCnrvaZcVIIFNHwGYS9raTKFT413sfIqJ1Bj76dl3dAnKsyQBD4VkwUjnkUXFQxia5/uGFAJ5lUXX6IAz5OWKH4tAHH3wAlUqFTp06YcCAAaJ1hdHrVuChe+IYJUc+FGzJEcuyKC8vR2lpKUpKSuo+1/+6uroaCoUCSqUSKpWqyZ81Gg1iY2MRGxsLtVrtdoxmsxl79uypS4jOOmhRqJWRkYEZM2bg3nvvrVv13F02mw2XLl2CQqFAcnJy3ZsKx/Mwma14+aUXGmxYbA8jUyCuQ3fcMnUmBgweAo1GeK5GoxErV63ChfPnUfbHNlgqXHQDSSSIaNMV8ohYWKvLob90Cjxrf22k+fPnY8GCBdAZzbhUZcaRY8fx1AN3o/hSrt362tRO0KZ3hkx9rYvKVFoA1qiDulU7SGRy8CwL5roZneGu/OiegGklat26dd0M21pvv/02nn76aZEiaozeiAMP3RPHKDnyIbGTI6PRaDfRKSktxZWSMpSUV6CsogpllVWoqNahWqcHpHJIZApI5AowMuFrRqaARCYHI1OAkUjBmvSwGaphM1QJn/XVYI06gG/aIpRKpRKxsbGIi4urS5jqfx0bGwuJRILt27dj+/btDtemckYikWDkyJGYOXMmMjMzYbVaG3xYLBYUFRXh1KlTdR9nz56F9erijMnJyRg+YjQybxgJeUIqvvthPU4f/h3VZ/6Atbqs8ePJVdCmZ0Gb2gkSmbAInkKhROcundG+XXv89NNPKC8XzuOsFpTs3wjWUG03dkYiRXyfEVBEXxtfZNNXo+zwTuEchoEiOhHyyFjwPA/OagFnNiCm6w2Qa6LAX70v1ppymMuLoYhOACOVwWbUQRGTCKlC5fHPM9yxFhOu7F0Lnm3G2lRN8NZbb2HUqFHo0KED9u7di5MnT+LGG29E165dsXz5cqxbtw6tWrXC448/jo7XLxkhMnojDjx0Txyj5MiHhg8fDqPRiJdf8cKKxfWwLIvcixfBQYLSikqUlFeirLIa5ZVVqKzWoUqnR5XOAIuNvZrkXEtuhGRH7vWmdZ7nwRpqriZLV5Omqx+c2ejx9RipHFKlGhKlGlKlBhKFCkJHFF9vejR/9ctrx3iOhbWmAjZdhZeeGYRWm7TOiGibXZfo1OJ5Hvq8U6g5dxg8a4UsIhba1hnQJLcHI/VsHoOh8BxYswGWqhKYS661AsRm3wh1yzaN6nNWC8qP7II2rRPULWgJhubgOQ76SzmISM8CZ7OBMxtgKi2Aseg8rNVlkMhViMroDUVMEmyGKlTnHIDNQTLbVElJSRg6dCgOHjxY1yKqVCqxdu1ajB492quPJQZ6Iw48dE8cczc5CprZaoGmpEqHB/691I+PqAU0Wsg0/r1pDMNApo2CTBsFJLZuUMbZrNdamGoTJovpWuKj1ECqUAmflWpIlWqPE4vr2Yw6mK5cgqkkD5bKK07Xm5FqoqCIFlpQeM4GjmXBszbwnA0SuQqR7bo16HK6/nlHpGdBnZQO1mxo0LrjKU3ytUUNjVfyUHFsDyLSOttNjABAIlcgoc9tTX68UGWpLoMiyvkMPFNJPmQRMZCpI8Cajag8sQ/m0nxUnz5gtz5nNaHy+C9eiS8lJQWPPPIIoqKiMHXqVCQlJTXa5LKoqAgnTpxAnz59EBUV5ZXHJYR4HyVHpMkkMjkUUfEu37Cag2EYtG6diviEeOScOgUjgIj0LKElwGKGqTQfxisXYSkvhkwbBUVMCyhik6CIaeGVbiWFNgqDbxuFhIQEGAwG5F3Kw4n6e3/Vw5qNkCobjrMS/mtjwF0dCK1ukQb18GsLOfbuLcwqKi4uhl6vx9q1a5ods7+Yy4tgqSpBZNtsr1+7+JfvoG6RBpkmCpaqEhgKzwqJsEQCZVwyZNooRLbtBolcmBVoKslH+dFdwnYvDANGKhMGn9vZl/Cee+7BmjVrYDAYGpV16tQJp05dm7KfkJCAL774AlOmTKmrP3HiRHzxxReIjIzEhQsXsG/fPkRHR2P06NGNWm6v3+SyZcuWaNmyZXN/PIQQH6PkiEClUkGlUoPjONhsNrCs7epnEfZlgzAduX2HDsjomIH27dvXDXa2jR6DnJwcHDp8COfPnYNEoYQmuX2DlhlviomJweTJk6HRaBAdHQ2GYdC/fz8UDxmCnT/9hNM5wkwmU2khdBeOwlpdhvi+I+uSxeSUFPzf1P+DRCLB8uXLkZ9/qcH1lUolhg69GVqtBu3atQUAHDp8CBfOn/c41l69eqNbt2wsXbqsLhFrrri4eKSkpODYsaNQMhzaRssxZPit2PD9OpRdOInbBvRBnz7jsHv3brRIScett4/HhBFDIdNEQ6aJhDwqARK5Am8vehUdkhOwavUa/Pffr0EeFQ95ZDxshiqYyy7jo02/g4mIg8ViQXFhAQamqBGnehT/+c9/UFlZic2bf70WFMfhpl6d0a1bN7z59ttQRCUgMyMDb7z2AkaPFprJeZ7H+fPnIZPJkJ6eDo7jYLFYIJPJ6pKVpUuX4o8//kCfPn3qLj1x4kR89tlnKC0txZo1ayCRSDBlyhS0adMGer0eFosFCkXDJRratm2Ltm3beuXnTQgJHDTmqAmGDx+OvKISGFJ6ih1KIxKJFBqNBlqtBmqNBlqNFhqNBhrtta8jtBpER2oRGxWJmKgIaBTCeJsqkw16sw08ro7z4XmwLAubjW2QNDX8nr16TPiwWCwwmUwwmowwGo0wGYWvTUYjjEYTTCZj3R5V9SUltUTHjI7I6NgRKSmtIZE4HztVXV2NI0eO4tChQ6iocG/VYoVCiYSEBCQkJCA+Pl5YZ+fc2QbT9XmeB8Mw6NCxIyZOmACVSoWqqqq65Kg+rrIYf25dgx9Xfl23OKFEqUZivzHo1rsf7rjjDsjlwpux1WrDxo0bcfjwYQA8FAoF7rrrrkZ7/JWUlGDx4iUwm2un5jNISUlBeXkZjEYjoqKikJCQKHR3ymRIaZ2CzlmdER8fB4VMgrPnL2LlypXQ63V11xzYtzc+e+FhHDx6AhazCd8u/Rw5x49i2IgxaD90Evb9eRQA0CWzI4ovnsWF0mqkp6WjY8cOdc/5tk5JaBOvdfkzLi0txdtvv42cnBz069cPs2fPhkp1rQXvgQcewJIlS+q+nzx5MlasWIFLlUZYbBxSotVQK5zPnvP2eIqysjLYbDYkJYXHMgXeRuNbAg/dE8doQLYP+TM5UiiU1yU6akRo1IjSqhEdoUVMZARioyMRHxOF+NgYxERFQCmTQiFloJBJoJRKIL/6WSGTQCGVQCphHA7ctnEcqo02VJmsqDRaUVXvo/62Ds3BcZyQQBlNsFjMiIyMQkSE6zdee3geKCwswJUrJeB5HlKpFBKJpMFnuVyO+Pg4RERENtoPl+eF1daLLuQg/8AOmPQ16DZiCtp07QWGYcDzfIPkSC6VoH1CBLq2ikKcVmhFYFkWe/fuxaZNm1BWVoYRY+9AZEYfFFebYb26mnStiooKVFRUomXLJGg0mrrj/dLjcLZEh3KDBcXFV3Dw4EHYWBu6dumK9u3bgWVZqORSWK/2EjEMgxvaxSNeqwTH84hUyaCRC0nF7uMX8cnarSivqkFGm9aYd9/tSIy69lj1WWwcCqqM4HgeabEayKUSbD5RhLyKa11OCVolxndPhsQLg/1rF2c9cuQIMjMzMXPmTI9fvOmFP7DQ/Qg8dE8co+TIh4YPH45qvRH3zZnn1evyHAe9XoeM9m2RlBCPxPg4RGnVkEslUF5NbORSx4mNr5ms7LWEyVQvcTJZGyQA9alkUqgVUmgVUqjlMmgU0msfcik0ChmkV1uJeF5os+J5CG1X9b6vMlmRW2bAxQp93Qak3hChlKNfeizaJ2jrfq48z+NsiR6/XyyH3mxFVVUV2rZKROdW0eiYGAGFzLMXmzK9BRuPX4bR2ji5lEoY9E+PQ9fkaJTozPjuaKHdlbPHd0tBlEqGc6V6mKwsUmLUaBnleEyVjeNgsnLQKqQe/75YWQ6H8itRqjMjWi1Hr9axLltz/Ile+AML3Y/AQ/fEMZqt5mNRWjWemDLKq9cM9F9olVyKlnJpozdlnuehM9tQZbTCxvNQy6XQKmRQySWQeel5xGkVaBuvBcsl4HK1CRfL9cgtM0BvabweDQMG8VoFWkar0DJShYQIBXgesHE8rCwHG8fDxnKIUMoQp1U0ahFhGAYdW0SgfaIWV6pNKC0GstonQ9rExRLjtQrc3rUVNhwvahBvvFaJoR0T61qgEiOU6J4SjcP5lQ3ObxOvRYtIYTuOLq3cm+Ekk0gQoWzaz14ulaBfelyTziWEkFBAyRFpNoZhEKmSI1Ild125maQSBq1j1Ggdo8agtjxK9RbklRtgtLJQyiRoFaVCi0iVx6079kgYBi0ilTBVeN76cr0YjQKTeqTgaGEVakw2JEerkNEisq7VrFaftFiYbRxOFglr7bSMUuGmDk1fRoAQQojnKDkiQYthGCRGKJEY4b/NaJtDJZe6bJGRMAyGtE9Av/RYWFkeEUr6EyWEEH+jV15CApBSJgXlRYQQIo7AG9hCCCGEECIiSo4IIYQQQuqh5IgQQgghpB5KjgghhBBC6qHkiBBCCCGkHkqOCCGEEELqoeSIEEIIIaQeSo4IIYQQQuqh5IgQQgghpB5KjgghhBBC6qHkiBBCCCGkHtq9qQmuXLkClmUxfPhwr1/bZrNBJqPbEkjongQeuieBhe5H4KF7Yt/ly5chlUpd1qOWoyZQKpU++6WjX+bAQ/ck8NA9CSx0PwIP3RP7ZDIZlEqly3oMz/O8H+IhhBBCCAkK1HJECCGEEFIPJUeEEEIIIfVQckQIIYQQUg8lR4QQQggh9VByFIR0Oh1eeOEF9O/fHz169MD999+PCxcuiB1WWDt37hwmTZqEzMxMsUMhAIqLi/HMM89g4MCB6Nu3L+69916cOHFC7LDC2qlTp/Dggw+if//+6NOnD6ZPn44//vhD7LAIgA0bNiAzMxNr164VO5SAQclREJo/fz5OnDiB1atXY8+ePUhOTsZDDz0Ei8UidmhhadOmTbj33nuRnp4udijkqkcffRQGgwEbNmzA7t27kZ6ejocffhgmk0ns0MKSXq/Hfffdh86dO2PXrl3Ys2cPMjMz8dBDD0Gv14sdXlirqKjAokWLoNFoxA4loFByFGQqKiqwceNG/P3vf0dqaiqioqLw4osvori4GD///LPY4YUlnU6H5cuX+2RRUOI5nU6HTp064YUXXkBcXBw0Gg0efPBBlJSU4Ny5c2KHF5bMZjP+9re/4cknn4RarYZGo8HkyZOh1+tRWFgodnhh7bXXXsPIkSMRGxsrdigBhZIjP1i7di169+6NzMxM5Ofn261jtVqxZMkSjBs3Dt27d0fv3r0xffp07Nixo0G9EydOgOM4ZGdn1x3TarVo3749jh075tPnESq8eT8AYPLkyUhNTfV12CHNm/ckIiICr7/+eoN7UlBQAIZhEB8f79PnEUq8eU/i4uIwefJkyOVyAEK356effors7Gy0bdvW588lFHj7dQsAdu/ejT///BNPP/20L0MPSpQc+VB5eTkef/xxvPDCC+A4zmndp59+Gv/85z+RnJyM+fPn429/+xv0ej0ee+wxfPvttw2uKZVKERER0eD8mJgYlJeX++R5hApf3A/SPP64JxUVFViwYAGmTp2Kli1bevsphBxf3pOamhp07doVN954I0pLS/HRRx/RSs4u+Op+6HQ6zJ8/H6+88gq0Wq0vn0Jw4onP3HzzzfzgwYP5PXv28NOnT+czMjL4S5cuNaq3detWPiMjg58zZ06D4waDgb/11lv5Hj168KWlpTzP8/wPP/zAd+7cudE17r33Xn7evHm+eSIhwhf3o77169fzGRkZPos/FPn6nly8eJG/7bbb+CeffJK3WCw+ex6hxNf3hOd5vqioiJ83bx4/dOhQvqqqyifPI1T46n68/PLL/HPPPVf3/dChQ/k1a9b47okEGWo58qEePXrghx9+wJAhQ5zWW7duHQBg5syZDY6r1WpMnToVBoMBmzZtAgAkJCTAZrOhpqamQd3y8nIkJCR4MfrQ44v7QZrHl/fkzz//xOTJk3HLLbfg3XffrevSIc754+8kKSkJ8+fPR0VFBbZs2eKdwEOUL+7HgQMHsHPnTjz//PO+CToEUHLkQ++88w7i4uJc1jt06BBUKhWysrIalfXs2ROA8EIPAF26dIFcLseRI0fq6lRXV+P8+fPo0aOHdwIPUb64H6R5fHVPTpw4gVmzZuFvf/sbnn32WUgk9FLnLl/ck127dmHYsGGNZtTyPE/dai744n6sXbsWVVVVGDFiBPr374/+/fvj8uXLePXVV/Hoo4969wkEKfqtFJlOp0N5eTnS09PtvoAnJycDAC5evAgAiIqKwoQJE/DPf/4T6enpiIyMxMKFC9GuXTsMHjzYr7GHIk/vB/E9T+8Jy7KYO3cuZsyYgcmTJ/s11nDh6T3p1q0bjEYjFi1ahGeeeQYSiQTvv/8+pFIpBg0a5NfYQ5Gn92Pu3LmYPXt2gzpTp07FjBkzMG7cON8HHAQoORJZ7RofjtaYqD2u0+nqjr344ot4/fXXMWnSJFitVvTr1w8ff/wxpFKp7wMOcU25HyNGjEBhYSF4ngeAupmEr776KsaPH+/DaMODp/fk0KFDyMnJwfnz5/Hhhx82qEv3xDs8vSdxcXFYsmQJ3nrrLQwfPhwcx6FTp0747LPPkJSU5J+gQ5in9yM6OhrR0dEN6kilUkRFRbnVShUOKDkKcLVvuAzD1B1TqVRYuHAhFi5cKFZYYcve/aAxE+K6/p706dMHOTk5YoYU9uz9nWRlZeGzzz4TK6SwZu9+XG/nzp3+CicoUEe8yCIjIwHA4Sqxtcdr6xHfovsReOieBB66J4GF7of3UXIkMo1Gg8TERBQVFYFl2UbltYt9tWnTxs+RhSe6H4GH7kngoXsSWOh+eB8lRwGgb9++sFgsdle43r9/PwCgX79+/g4rbNH9CDx0TwIP3ZPAQvfDuyg5CgBTpkwBACxevLjB8erqaqxYsQIxMTEYOXKkGKGFJbofgYfuSeChexJY6H54Fw3I9pGCgoIGGXzt1h579uypmw2QkpKC7OxsDBw4EJMmTcKaNWswa9YsjBw5EgaDAV999RVKS0vx7rvvUl9xM9H9CDx0TwIP3ZPAQvdDPAxfO4ydeNXatWtdrj46YcIEvPHGGwAAjuPwzTffYOXKlcjNzYVcLkf37t3x6KOPom/fvv4IOaTR/Qg8dE8CD92TwEL3QzyUHBFCCCGE1ENjjgghhBBC6qHkiBBCCCGkHkqOCCGEEELqoeSIEEIIIaQeSo4IIYQQQuqh5IgQQgghpB5KjgghhBBC6qHkiBBCCCGkHkqOCCGEEELqoeSIEEIIIaQeSo4IIYQQQuqh5IgQQgghpB5KjgghhBBC6vn/+Upb5S0xOooAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.plot(loss_all[:,0], color = 'black', label = 'F-VI', linewidth = linewidth)\n",
    "plt.plot(loss_all[:,1], color = 'red', label = 'Const', linewidth = linewidth)\n",
    "plt.plot(loss_all[:,5], alpha = 0.5, label = \"i=1, k=20\", linewidth = linewidth)\n",
    "plt.plot(loss_all[:,7], alpha = 0.5, label = \"i=2, k=4\", linewidth = linewidth)\n",
    "\n",
    "plt.legend()\n",
    "plt.yscale(\"log\")\n",
    "plt.xscale(\"log\")\n",
    "plt.gca().invert_yaxis()\n",
    "plt.ylabel(\"Elbo\")\n",
    "plt.grid(which = 'major', visible = 'true', c = 'grey', alpha = 0.25)\n",
    "\n",
    "plt.savefig(\"to_discuss/k_vae_time_series.pdf\", bbox_inches = 'tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 533,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_window(loss, index_center, window):\n",
    "    return np.mean(loss[(index_center - window):(index_center + window)])\n",
    "\n",
    "def sd_window(loss, index_center, window):\n",
    "    return np.std(loss[(index_center - window):(index_center + window)]) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 534,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_iter_optimizer = int(30000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 535,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FVI:  553.0837867431641  +/-  0.06985303834506525\n",
      "Const:  2933.555895751953  +/-  6.7674346535909455\n",
      "k = 4:  891.9649091796875  +/-  0.06985303834506525\n",
      "k = 6:  892.8204157714844  +/-  6.7674346535909455\n",
      "k = 10:  886.1229639892579  +/-  2.5709975260837434\n",
      "k = 20:  886.220542541504  +/-  2.5966791410048455\n",
      "k = 2:  2933.6476901855467  +/-  2.588699204945348\n",
      "k = 4:  553.0766517333984  +/-  2.5385311546707383\n",
      "k = 6:  553.0869779052734  +/-  7.0208861080178\n",
      "k = 10:  553.116260925293  +/-  0.10228968027303144\n"
     ]
    }
   ],
   "source": [
    "window = 500\n",
    "index_center = n_iter_optimizer - window\n",
    "print(\"FVI: \", mean_window(loss_all[:, 0], index_center, window), \" +/- \",\n",
    "      sd_window(loss_all[:, 0], index_center, window))\n",
    "\n",
    "print(\"Const: \", mean_window(loss_all[:, 1], index_center, window), \" +/- \",\n",
    "      sd_window(loss_all[:, 1], index_center, window))\n",
    "\n",
    "for i in range(nn_widths.size):\n",
    "  print(\"k = \" + str(nn_widths[i]) + \": \", mean_window(loss_all[:, i + 2],\n",
    "        index_center, window), \" +/- \",\n",
    "        sd_window(loss_all[:, i], index_center, window))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 536,
   "metadata": {},
   "outputs": [],
   "source": [
    "def iter_convergence(loss, rel_tol = 0.01, init_iteration = 10, window = 10,\n",
    "                     grid = 10, max_iteration = 20000, oracle = True, set_target = 0):\n",
    "    iteration = init_iteration\n",
    "    tol = 0\n",
    "    diff = tol + 1\n",
    "    while(diff > tol and iteration < max_iteration):\n",
    "        iteration += grid\n",
    "        if (oracle):\n",
    "            target = set_target\n",
    "        else:\n",
    "            target = mean_window(loss, iteration - grid, window)\n",
    "\n",
    "        tol = rel_tol * abs(target)\n",
    "        diff = abs(mean_window(loss, iteration, window) - target)\n",
    "        # print(\"diff:\", diff)\n",
    "        # print(\"tol:\", tol)\n",
    "    return iteration\n",
    "\n",
    "# Q: should we report convergence to F-VI's solution or for each individual solution?\n",
    "def iter_conv_wrapper(loss, loss_benchmark, rel_tol = 0.01, window = 10):\n",
    "    return iter_convergence(loss, oracle = True, rel_tol = rel_tol,\n",
    "                            set_target = mean_window(loss_benchmark, 19500, window))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 537,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get results across multiple seeds\n",
    "init_seed = 2013\n",
    "num_seed = 10\n",
    "init_manifold = False\n",
    "n_algorithms = 10\n",
    "num_tol = 3\n",
    "rel_tol = np.array([0.5, 0.1, 0.01])\n",
    "\n",
    "iter_to_conv = np.empty((num_seed, n_algorithms, num_tol))\n",
    "Elbo_final = np.empty((num_seed, n_algorithms)) \n",
    "\n",
    "for i in range(num_seed):\n",
    "    loss_all = np.load(\"deliv/timeseries_loss_\" + str(init_seed + i) + \".npy\")\n",
    "\n",
    "    for j in range(n_algorithms):\n",
    "        Elbo_final[i, j] = mean_window(loss_all[:, j], n_iter_optimizer - window, window)\n",
    "        for k in range(num_tol):\n",
    "            # Use F-VI as the benchmark for the asymptotic ELBO\n",
    "            iter_to_conv[i, j, k] = iter_conv_wrapper(loss_all[:, j], loss_all[:, 0], rel_tol[k])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 538,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean iter:\n",
      " [[14405. 15925. 17595.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [ 5405.  6870.  9085.]\n",
      " [ 4780.  6305.  7780.]\n",
      " [ 3360.  5025.  7235.]\n",
      " [ 3210.  4650.  7145.]]\n",
      "min iter:\n",
      " [[14240. 15750. 17330.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [ 3330.  5070.  7500.]\n",
      " [  240.  2220.  3620.]\n",
      " [  310.  3240.  5460.]\n",
      " [  430.  2160.  3470.]]\n",
      "max iter:\n",
      " [[14530. 16060. 17800.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [20000. 20000. 20000.]\n",
      " [15060. 17490. 20000.]\n",
      " [ 6750.  8050.  9200.]\n",
      " [ 6520.  7820.  9510.]]\n"
     ]
    }
   ],
   "source": [
    "print(\"mean iter:\\n\", np.median(iter_to_conv, axis = 0))\n",
    "print(\"min iter:\\n\", np.min(iter_to_conv, axis = 0))\n",
    "print(\"max iter:\\n\", np.max(iter_to_conv, axis = 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 539,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 553.07710651, 2933.55018848,  891.96490918,  891.98230627,\n",
       "        886.06733527,  885.16987256,  553.04929572,  553.06468109,\n",
       "        553.07584314,  553.06945428])"
      ]
     },
     "execution_count": 539,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nanmin(Elbo_final, axis = 0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "5b7e68a6523d9654d119e3ab45bb2c5c51f02effa6643d52b99addcee48c066f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}