example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "env: FLOX_STRICT_TYPING=True\n"
     ]
    }
   ],
   "source": [
    "%env FLOX_STRICT_TYPING=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:using flox with strict typing\n",
      "WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    }
   ],
   "source": [
    "import flox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "\"\"\" core system definitions \"\"\"\n",
    "\n",
    "\n",
    "from functools import partial\n",
    "from typing import Callable, TypeVar\n",
    "\n",
    "import jax\n",
    "from jax import numpy as jnp\n",
    "from jax_dataclasses import pytree_dataclass\n",
    "from jaxtyping import Array, Float  # type: ignore\n",
    "import jax.scipy as jscipy\n",
    "from flox.flow import Transformed\n",
    "\n",
    "from flox.util import key_chain\n",
    "\n",
    "T = TypeVar(\"T\")\n",
    "\n",
    "Scalar = Float[Array, \"\"] | float\n",
    "VectorN = Float[Array, \"N\"]\n",
    "\n",
    "Potential = Callable[[T], Scalar]\n",
    "\n",
    "def double_well_potential(\n",
    "    x: Scalar,\n",
    "    a: Scalar = 2,\n",
    "    b: Scalar = 5,\n",
    "    c: Scalar = 0.5\n",
    ") -> Scalar:\n",
    "    \"\"\" The simple double potential \"\"\"\n",
    "    return jnp.sum(a * x ** 4 - b * x ** 2 + c * x ** 3)\n",
    "\n",
    "\n",
    "def harmonic_potential(x: Scalar) -> Scalar:\n",
    "    \"\"\" The simple harmonic potential \"\"\"\n",
    "    return jnp.sum(x ** 2)\n",
    "\n",
    "\n",
    "@pytree_dataclass\n",
    "class State:\n",
    "    \"\"\" State holding the reaction coordinate\n",
    "        and auxiliary noise.\n",
    "    \"\"\"\n",
    "    pos: VectorN\n",
    "    aux: VectorN\n",
    "\n",
    "\n",
    "def joint_base_potential(s: State) -> Scalar:\n",
    "    \"\"\" combine both potentials into one \"\"\"\n",
    "    return jnp.sum(harmonic_potential(s.pos) + harmonic_potential(s.aux))\n",
    "\n",
    "\n",
    "def joint_target_potential(s: State) -> Scalar:\n",
    "    \"\"\" combine both potentials into one \"\"\"\n",
    "    return jnp.sum(double_well_potential(s.pos) + harmonic_potential(s.aux))\n",
    "\n",
    "\n",
    "def sample_base(key: jax.random.PRNGKeyArray) -> Transformed[State]:\n",
    "    chain = key_chain(key)\n",
    "    pos = jax.random.normal(next(chain), shape=(1,))\n",
    "    aux = jax.random.normal(next(chain), shape=(1,))\n",
    "    logprob = -jnp.sum(\n",
    "        jscipy.stats.norm.logpdf(pos)\n",
    "        +\n",
    "        jscipy.stats.norm.logpdf(aux)\n",
    "    )\n",
    "    return Transformed(State(pos, aux), logprob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "\"\"\" init random state. \"\"\"\n",
    "\n",
    "seed = 42\n",
    "chain = key_chain(seed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "\"\"\" plotting tools \"\"\"\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "\n",
    "def plot_potential(potential: Potential, left=-2, right=1.5, res=100):\n",
    "    \"\"\" plot the potential \"\"\"\n",
    "    xs = jnp.linspace(left, right, res)\n",
    "    us = jax.vmap(potential)(xs)\n",
    "    plt.plot(xs, us)\n",
    "\n",
    "\n",
    "def density(potential: Potential, left=-3, right=3, res=1000) -> tuple[Array, Array]:\n",
    "    \"\"\" get the density for plotting \"\"\"\n",
    "    xs = jnp.linspace(left, right, res)\n",
    "    us = jax.vmap(potential)(xs)\n",
    "    pis = jnp.exp(-us)\n",
    "    Z = jnp.trapz(pis, xs)\n",
    "    ps = pis / Z\n",
    "    return xs, ps\n",
    "\n",
    "\n",
    "def plot_samples(\n",
    "    s: State,\n",
    "    potentials: tuple[Potential, Potential] | None = None,\n",
    "    weights: Float[Array, \"N\"] | None = None,\n",
    "    bins: tuple[int, int] = (50, 50),\n",
    "    range: tuple[tuple[float, float], tuple[float, float]] = ((-3, 3), (-3, 3)),\n",
    "    label: str | None = None\n",
    "):\n",
    "    \"\"\" plot samples \"\"\"\n",
    "\n",
    "    mosaic = [\n",
    "        [\"pos\", \"aux\"],\n",
    "        [\"joint\", \"joint\"]\n",
    "    ]\n",
    "    if weights is not None:\n",
    "        mosaic.append(\n",
    "            [\"joint_weighed\", \"joint_weighed\"]\n",
    "        )\n",
    "    fig, axes = plt.subplot_mosaic(mosaic, figsize=(8, 8 + (1 if weights is not None else 0)))\n",
    "\n",
    "    for i, (ax, data, title) in enumerate(zip((\"pos\", \"aux\"), (s.pos, s.aux), (\"pos\", \"aux\"))):\n",
    "        axes[ax].hist(data[:, 0], bins=bins[i], density=True, histtype=\"step\", range=range[i],  # type: ignore\n",
    "                      label=\"samples\", color=\"blue\")\n",
    "        axes[ax].set_title(title)  # type: ignore\n",
    "\n",
    "    if weights is not None:\n",
    "        for i, (ax, data, title) in enumerate(zip((\"pos\", \"aux\"), (s.pos, s.aux), (\"pos\", \"aux\"))):\n",
    "            axes[ax].hist(data[:, 0], bins=bins[i], density=True, histtype=\"step\", range=range[i],  # type: ignore\n",
    "                          label=\"samples (reweighed)\", color=\"green\", weights=weights)\n",
    "    if potentials is not None:\n",
    "        for ax, pot in zip((\"pos\", \"aux\"), potentials, strict=True):\n",
    "            xs, ps = density(pot)\n",
    "            axes[ax].plot(xs, ps, color=\"red\", label=\"density\")  # type: ignore\n",
    "    \n",
    "    for ax in (\"pos\", \"aux\"):\n",
    "        axes[ax].legend()  # type: ignore\n",
    "    \n",
    "    axes[\"joint\"].hist2d(s.pos[:, 0], s.aux[:, 0], bins=bins, range=range)  # type: ignore\n",
    "    axes[\"joint\"].set_title(\"joint\")  # type: ignore\n",
    "\n",
    "    if weights is not None:\n",
    "        axes[\"joint_weighed\"].hist2d(s.pos[:, 0], s.aux[:, 0], bins=bins, range=range, weights=weights)  # type: ignore\n",
    "        axes[\"joint_weighed\"].set_title(\"joint (reweighed)\")  # type: ignore\n",
    "    if label is not None:\n",
    "        fig.suptitle(label)  # type: ignore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sample from target potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jonas/miniconda3/envs/flox-debug/lib/python3.10/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "\"\"\" use numpyro sampler \"\"\"\n",
    "\n",
    "import numpyro as npy  # type: ignore\n",
    "from jaxtyping import PyTree  # type: ignore\n",
    "\n",
    "\n",
    "def sample(\n",
    "    key: jax.random.PRNGKeyArray,\n",
    "    potential_fn: Potential[T],\n",
    "    num_samples: int,\n",
    "    init: T,\n",
    "    num_chains: int,\n",
    "    num_warmup: int\n",
    ") -> T:\n",
    "    \"\"\" samples from potential using numpyro's no-u-turn sampler (HMC). \"\"\"\n",
    "    mcmc = npy.infer.MCMC(npy.infer.NUTS(potential_fn=potential_fn), num_warmup=num_warmup, num_samples=num_samples, num_chains=num_chains)\n",
    "    mcmc.run(key, init_params=init)\n",
    "    return mcmc.get_samples()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 10100/10100 [00:12<00:00, 789.94it/s, 21 steps of size 1.72e-01. acc. prob=0.96] \n"
     ]
    }
   ],
   "source": [
    "num_chains = 1\n",
    "\n",
    "init = State(\n",
    "    jax.random.normal(next(chain), (num_chains,)),\n",
    "    jax.random.normal(next(chain), (num_chains,))\n",
    ")\n",
    "\n",
    "npy.set_host_device_count(num_chains)\n",
    "data = sample(\n",
    "    next(chain),\n",
    "    joint_target_potential,\n",
    "    num_samples=10_000,\n",
    "    init=init,\n",
    "    num_chains=num_chains,\n",
    "    num_warmup=100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqQAAALjCAYAAADEGxd9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACcmElEQVR4nOzdeXxU1fnH8c9kkkwSQhIgEBYTUEA2ZQuLqChqFFFRaq1Ua0FUahHc0P5qXEBqa2zdaBVFrYAWKVjrUquiFotURZAEFBFQZAkiCYQlgSRkmbm/P24yEAiQO8zMneX7fr3mlZnJuXOfyXLz5JzznOMwDMNARERERMQmMXYHICIiIiLRTQmpiIiIiNhKCamIiIiI2EoJqYiIiIjYSgmpiIiIiNhKCamIiIiI2EoJqYiIiIjYSgmpiIiIiNhKCamIiIiI2EoJqYhIEzz44IM4HA67wxARiUhKSEVEAuiZZ55hzpw5dochIhLSlJCKiASQElIRkeNTQioiIiIitlJCKiJymE8++YSBAweSkJBA586dee65545oM3v2bM4//3zatGmDy+WiZ8+ePPvssw3adOrUiTVr1vDxxx/jcDhwOBwMGzYMgN27d3P33Xdz+umnk5ycTEpKCiNGjODLL78MxlsUEQkpDsMwDLuDEBEJFatXr2bw4MG0bt2aCRMmUFtby9NPP01GRgZfffUV9ZfMQYMG0atXL/r06UNsbCxvv/02H3zwAU8//TQTJ04E4M033+TWW28lOTmZ++67D4CMjAwuvPBCVqxYwc9//nN+9rOfcfLJJ1NcXMxzzz3H/v37+eabb2jfvr1tXwMRkWBTQioicoif/OQnLFy4kPXr15OVlQXA2rVrOf3003G73d6EtLKyksTExAbHXnzxxXz33Xd8//333udOO+000tPTWbx4cYO2VVVVxMXFERNzcKBq8+bNdO/enfvuu48HHnggQO9QRCT0aMheRKSO2+3m/fffZ9SoUd5kFKBHjx4MHz68QdtDk9HS0lJKSko499xz2bhxI6Wlpcc9l8vl8iajbrebXbt2kZycTLdu3SgoKPDTOxIRCQ9KSEVE6uzcuZPKykq6du16xOe6devW4PGnn35KTk4OzZo1Iy0tjdatW3PvvfcCNCkh9Xg8PPnkk3Tt2hWXy0V6ejqtW7fmq6++atLxIiKRRAmpiIhF33//PRdccAElJSU88cQTvPPOO3z44YfceeedgJlsHs/DDz/M5MmTOeecc5g7dy7vv/8+H374Ib169WrS8SIikSTW7gBEREJF69atSUxM5Lvvvjvic+vXr/fef/vtt6mqquJf//pXg6H9//73v0ccd7TdnV577TXOO+88XnzxxQbP7927l/T0dF/fgohIWFIPqYhIHafTyfDhw3nzzTcpLCz0Pr927Vref//9Bu0ADq0JLS0tZfbs2Ue8ZrNmzdi7d2+j5zq8pvQf//gH27ZtO9G3ISISdtRDKiJyiGnTprFw4UKGDh3KLbfcQm1tLU899RS9evXiq6++AuCiiy4iPj6ekSNHcvPNN7N//35eeOEF2rRpw/bt2xu8XnZ2Ns8++yy///3v6dKlC23atOH888/nsssu43e/+x3jxo3jzDPPZPXq1bzyyiuccsopdrxtERFbadknEZHDLFmyhMmTJ7N69WpOOukk/u///o/t27czbdo0b6/m22+/zf3338+3335L27ZtmTBhAq1bt+aGG25g06ZNdOrUCYDi4mJuvPFGlixZwr59+zj33HNZvHgxVVVV3HfffcybN4+9e/fSv39/HnvsMe655x6AI5aJEhGJZEpIRURERMRWmkMqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiHvwQcfxOFwsG7dOq6++mpSUlJo1aoVt99+OwcOHPC2q62t5aGHHqJz5864XC46derEvffeS1VVVYPXW7FiBcOHDyc9PZ3ExEROPvlkbrjhhmC/LRGRoNmyZQu33HIL3bp1IzExkVatWvGzn/2MzZs3N2hXf7093Jw5c3A4HN72H330ETExMUyZMqVBu3nz5uFwOHj22WcD9VYkQsXaHYBIU1199dV06tSJvLw8Pv/8c/7yl7+wZ88eXn75ZQBuuukmXnrpJa666iruuusuli1bRl5eHmvXruWNN94AYMeOHVx00UW0bt2ae+65h7S0NDZv3szrr79u51sTEQmoL774gs8++4yf//znnHTSSWzevJlnn32WYcOG8c0335CUlGTp9c4//3xuueUW8vLyGDVqFP3792f79u3ceuut5OTk8Otf/zpA70QiliES4qZOnWoAxuWXX97g+VtuucUAjC+//NJYtWqVARg33XRTgzZ33323ARgfffSRYRiG8cYbbxiA8cUXXwQtfhERu1VUVBzx3NKlSw3AePnll73P1V9vDzd79mwDMDZt2uR9rry83OjSpYvRq1cv48CBA8all15qpKSkGFu2bAnIe5DIpiF7CRsTJ05s8PjWW28F4N133+Xdd98FYPLkyQ3a3HXXXQC88847AKSlpQHw73//m5qamkCGKyISMhITE733a2pq2LVrF126dCEtLY2CggKfXjMpKYk5c+awdu1azjnnHN555x2efPJJsrKy/BW2RBElpBI2unbt2uBx586diYmJYfPmzWzZsoWYmBi6dOnSoE3btm1JS0tjy5YtAJx77rn89Kc/Zdq0aaSnp3PFFVcwe/bsI+aZiohEksrKSqZMmUJmZiYul4v09HRat27N3r17KS0t9fl1zzrrLCZMmMDy5csZPny45uOLz5SQSthqbOJ9Y88d/vnXXnuNpUuXMmnSJLZt28YNN9xAdnY2+/fvD1SoIiK2uvXWW/nDH/7A1VdfzauvvsoHH3zAhx9+SKtWrfB4PN52R7uGut3uRp+vqqpi8eLFAHz//fdUVFT4PXaJDkpIJWx89913DR5v2LABj8dDp06d6NixIx6P54g2xcXF7N27l44dOzZ4/owzzuAPf/gDK1as4JVXXmHNmjXMnz8/4O9BRMQOr732GmPHjuXxxx/nqquu4sILL+Tss89m7969Ddq1aNEC4Ijn60eZDjd16lTWrl3LY489xqZNm7jnnnsCEb5EASWkEjZmzJjR4PFTTz0FwIgRI7jkkksAmD59eoM2TzzxBACXXnopAHv27MEwjAZt+vbtC6BhexGJWE6n84hr31NPPXVEz2fnzp0BWLJkife58vJyXnrppSNec9myZTz22GPccccd3HXXXfzmN7/h6aef5uOPPw7AO5BIp2WfJGxs2rSJyy+/nIsvvpilS5cyd+5crr32Wvr06QPA2LFjef7559m7dy/nnnsuy5cv56WXXmLUqFGcd955ALz00ks888wz/OQnP6Fz587s27ePF154gZSUFG9SKyISaS677DL+9re/kZqaSs+ePVm6dCn/+c9/aNWqVYN2F110EVlZWdx444385je/wel0MmvWLFq3bk1hYaG33YEDBxg7dixdu3blD3/4AwDTpk3j7bffZty4caxevZpmzZoF9T1KmLO7zF/keOqXIfnmm2+Mq666ymjevLnRokULY9KkSUZlZaW3XU1NjTFt2jTj5JNPNuLi4ozMzEwjNzfXOHDggLdNQUGBcc011xhZWVmGy+Uy2rRpY1x22WXGihUr7HhrIiJBsWfPHmPcuHFGenq6kZycbAwfPtxYt26d0bFjR2Ps2LEN2ubn5xuDBw824uPjjaysLOOJJ544YtmnO++803A6ncayZcsaHLtixQojNjbWmDBhQpDemUQKh2Ec1ocvEmIefPBBpk2bxs6dO0lPT7c7HBEREfEzzSEVEREREVspIRURERERWykhFRERERFbKSGVkPfggw9iGIbmj0pEmTFjBp06dSIhIYHBgwezfPnyY7bfu3cvEydOpF27drhcLk499VTvlrkiIuFOyz6JiATZggULmDx5MjNnzmTw4MFMnz6d4cOHs379etq0aXNE++rqai688ELatGnDa6+9RocOHdiyZQtpaWnBD15EJADCosre4/Hw448/0rx58+NuDSki4gvDMNi3bx/t27cnJiawg0eDBw9m4MCBPP3004B5jcvMzOTWW29tdKebmTNn8uijj7Ju3Tri4uJ8OqeuoyISaCdyHQ2LhPSHH34gMzPT7jBEJAps3bqVk046KWCvX11dTVJSEq+99hqjRo3yPj927Fj27t3LW2+9dcQxl1xyCS1btiQpKYm33nqL1q1bc+211/Lb3/4Wp9PZ6Hmqqqoa7D62bds2evbs6ff3IyJyOF+uo2ExZN+8eXPAfIMpKSk2RyMikaisrIzMzEzv9SZQSkpKcLvdZGRkNHg+IyODdevWNXrMxo0b+eijj/jFL37Bu+++y4YNG7jllluoqalh6tSpjR6Tl5fHtGnTjnhe11ERCZQTuY6GRUJaP7yUkpKiC6mIBFQoDmd7PB7atGnD888/j9PpJDs7m23btvHoo48eNSHNzc1l8uTJ3sf1fyh0HRWRQPPlOhoWCamISKRIT0/H6XRSXFzc4Pni4mLatm3b6DHt2rUjLi6uwfB8jx49KCoqorq6mvj4+COOcblcuFwu/wYvIhIgWvZJRCSI4uPjyc7OZtGiRd7nPB4PixYtYsiQIY0ec9ZZZ7FhwwY8Ho/3uW+//ZZ27do1moyKiIQbJaQiIkE2efJkXnjhBV566SXWrl3LhAkTKC8vZ9y4cQCMGTOG3Nxcb/sJEyawe/dubr/9dr799lveeecdHn74YSZOnGjXWxAR8SsN2Yv4idvtpqamxu4w5CgOH/K20+jRo9m5cydTpkyhqKiIvn37snDhQm+hU2FhYYMlUzIzM3n//fe588476d27Nx06dOD222/nt7/9rV1vQSQgdB0NbYG8jobFsk9lZWWkpqZSWlqqyfgScgzDoKioiL1799odihxHWloabdu2bXTCfaRfZyL9/Ul403U0fATqOmq5h3TJkiU8+uij5Ofns337dt54440Ga+kdy6effsq5557LaaedxqpVq6yeWiQk1V9E27RpQ1JSUkhWaUc7wzCoqKhgx44dgFkkJCKhQ9fR0Bfo66jlhLS8vJw+ffpwww03cOWVVzb5uL179zJmzBguuOCCI6pLRcKV2+32XkRbtWpldzhyDImJiQDs2LGDNm3ahMzwvUi003U0fATyOmo5IR0xYgQjRoywfKJf//rXXHvttTidTt58881jtj18h5GysjLL5xMJhvq5TklJSTZHIk1R/32qqalRQioSInQdDS+Buo4Gpcp+9uzZbNy48agLOB8uLy+P1NRU703bhkqo0/BSeND3SSR06fczPATq+xTwhPS7777jnnvuYe7cucTGNq1DNjc3l9LSUu9t69atAY5SREREROwS0GWf3G431157LdOmTePUU09t8nHaYSTwCguhpKThc+npkJVlTzyRqLGvcaBEyvfu+uuvZ+/evced1iMiwWH33wpdR60L1+toQBPSffv2sWLFClauXMmkSZMAc0cSwzCIjY3lgw8+4Pzzzw9kCNKIwkLo0QMqKho+n5QEa9dGxi+k3Y72NQ4Ufe9ExN/s/luh62h0CWhCmpKSwurVqxs898wzz/DRRx/x2muvcfLJJwfy9HIUJSXmL/jcueYvO5i/hNddZ35Ov4wnrrGvcaDoeycigWD33wpdR6OL5Tmk+/fvZ9WqVd51RDdt2sSqVasoLCwEzPmfY8aMMV88JobTTjutwa1NmzYkJCRw2mmn0axZM/+9E7GsRw/o39+8BfqXPVod+jUO1M3X791rr73G6aefTmJiIq1atSInJ4fy8nK++OILLrzwQtLT00lNTeXcc8+loKCgwbEOh4PnnnuOyy67jKSkJHr06MHSpUvZsGEDw4YNo1mzZpx55pl8//333mMefPBB+vbty3PPPUdmZiZJSUlcffXVlJaWHjVGj8dDXl4eJ598MomJifTp04fXXnvN+/k9e/bwi1/8gtatW5OYmEjXrl2ZPXu2b18QETEdOADz59P+qVwe4Hdk1y6jfz/Dtr8Vuo5Gx3XUckK6YsUK+vXrR79+/QBzT+Z+/foxZcoUALZv3+5NTiVMfPstne77BbfxZwj9jbvED7Zv384111zDDTfcwNq1a1m8eDFXXnklhmGwb98+xo4dyyeffMLnn39O165dueSSS9i3b1+D13jooYcYM2YMq1atonv37lx77bXcfPPN5ObmsmLFCgzD8E7VqbdhwwZeffVV3n77bRYuXMjKlSu55ZZbjhpnXl4eL7/8MjNnzmTNmjXceeedXHfddXz88ccAPPDAA3zzzTe89957rF27lmeffZb09HT/f8FEosU778App8A119B2ziP8jql0v/4MGDECfvzR7uhCiq6jfmaEgdLSUgMwSktL7Q4lIuTnGwaYHw232zD69TOfAOPbGR/YHV5YqaysNL755hujsrKywfMNvsYB5su58vPzDcDYvHnzcdu63W6jefPmxttvv+19DjDuv/9+7+OlS5cagPHiiy96n/v73/9uJCQkeB9PnTrVcDqdxg8//OB97r333jNiYmKM7du3G4ZhGGPHjjWuuOIKwzAM48CBA0ZSUpLx2WefNYjnxhtvNK655hrDMAxj5MiRxrhx45r8vo/2/TKMyL/ORPr7Ez+YOdP7t8DIzDR2/OwW4++MNtxx8eZzWVnG129+5/drm66jpmi/jgZlHVIJYcuWwcqV3oetX33axmAkWPr06cMFF1zA6aefzs9+9jNeeOEF9uzZA0BxcTHjx4+na9eupKamkpKSwv79+48Y+ejdu7f3fkZGBgCnn356g+cOHDjQYGOLrKwsOnTo4H08ZMgQPB4P69evPyLGDRs2UFFRwYUXXkhycrL39vLLL3uHsCZMmMD8+fPp27cv//d//8dnn33mh6+OSBT6178wJkwAYOdPb2blq9/xwRUzuIb5rJ3/FZx6KhQW0vm2S0hjj83BhgZdR/0roEVNEgbeeQeAA526k7B5Hc2XL4LqaoiPtzkwCSSn08mHH37IZ599xgcffMBTTz3Ffffdx7Jly5gwYQK7du3iz3/+Mx07dsTlcjFkyBCqq6sbvEZcXJz3fv1CyY095/F4fIpx//79ALzzzjsNLr6Ad1m4ESNGsGXLFt59910+/PBDLrjgAiZOnMhjjz3m0zlFotIPP+AeOw6nYfAsv+aWfz4D/zR/f5OSoPmAbrBkCQwaRELhd7zAeOC1Y79mFNB11L/UQxrtPv8cgOJfTGYHrXFWlsPy5TYHJcHgcDg466yzmDZtGitXriQ+Pp433niDTz/9lNtuu41LLrmEXr164XK5KPHTQoCFhYX8eMg8tM8//5yYmBi6det2RNuePXvicrkoLCykS5cuDW6H7t7WunVrxo4dy9y5c5k+fTrPP/+8X2IViQqGAePH49y7my8YQNqcP5Of7yA/H/LzD1kGKSMD3ngDwxnLVfyTtEX/tDvykKDrqP+ohzSaGQZ8+SUAld378SVncDlvm0P4Z59tc3ASSMuWLWPRokVcdNFFtGnThmXLlrFz50569OhB165d+dvf/saAAQMoKyvjN7/5DYmJiX45b0JCAmPHjuWxxx6jrKyM2267jauvvpq2bdse0bZ58+bcfffd3HnnnXg8Hs4++2xKS0v59NNPSUlJYezYsUyZMoXs7Gx69epFVVUV//73v+mhJSNEjqt+wfmUJe/QZeFC3LHx/KL2FeafHk///kc5qH9/iq7/Le1e/AMn/elWmDjC7EKNUrqO+pcS0igWV7LdvCLFxFB5Si++pI+ZkNYlqXLi1q4NzXOkpKSwZMkSpk+fTllZGR07duTxxx9nxIgRtG3bll/96lf079+fzMxMHn74Ye6++26/xNqlSxeuvPJKLrnkEnbv3s1ll13GM888c9T2Dz30EK1btyYvL4+NGzeSlpZG//79uffeewGIj48nNzeXzZs3k5iYyNChQ5k/f75fYhWJVPULzldV1LKGuwB4rPYOtiWdyvGKq4tueoADL77CySWb4S9/gXvuCXi8uo42FLHXUctlUDZQdah/1VcTfvfUe+adHj2M/HzDuIpXzccDB9odYtg4WrXhli2GkZR0sGA10LekJPOcoWzq1KlGnz59bI1BVfaR+/6k6er/BnxyyyuGAUZ1Wrqx8uPSJl1D8vMN4zpeNl8gNdUw9uw54Xh0HW26SL6Oqoc0irk2rzPv9OoFwJf0AcDz1des+sINTicQOfv7BlNWlvkft/ZgFpHQZND/g0cAiLvrdvqek9LkI+dxLc+f8giJG7+Bv/4V/NTzdzhdR6OLEtIo5vpxs3nnlFNIT4ftiZ2pqozHVVXJTwb9QCEdAe3v66usLH3NRCQ0Xco7JG5YDcnJMHGipWM9ONlx3V10/N2N5rD97bfDIZXh/qTraPRQlX0Ui/9xk3mnUyeysmDNOidGx04AfPjs9+Tnm3sIV1QE7z9UiVwPPvigd8thEbHXbfzFvPPrX0OLFpaP/6zTtdS0bANbt7LpsX+iDRqDI5Kvo0pIo1h8fQ9pp06A+V9oQs/OAJzq/F573IuIRKD4HzZyER+aD46x5WRj0tPNUbNrb0jgd7vNLS233fs0PXqgpFROiBLSqGUcHLKvS0gBcw9jgLodHEREJLKkv/lXAMrOuAhOPtnSsfXzOvPz4WcLb8JwOjmbT+lQ8a1G0uSEKCGNUmnsxVletxVZx44HP9HZ7CFl48bgByUiIoFVU0Orf80CoOTKX/n0EllZ0L8/9B7eDsfFFwNwPXP8FaFEKSWkUSqLurGV+vGXevUJqXpIRUQiz6JFxO0qZget2XvO5Sf+euPGATCWl8DtPvHXk6ilhDRKtaXIvNO+fcNP1G8l9sMPwQ1IREQCr27B81e52j+V8SNHUpvaig78SPMvPjrx15OopYQ0SmVQbN45fKuxk04yP+7YAdXVwQ1KREQC58ABeOMNAObzc/+8Znw8e3KuAiBt0Wv+eU2JSkpIo5S3hzQjo+En0tMhPt68/+OPwQ1KbDds2DDuuOOOoJzrwQcfpG/fvkE5l4gACxdCWRnVbTrwGWf67WX3XlCXkH70OtTW+u11w5Wuo75RQhqlvAnp4T2kDgd06GDe37YtuEFJVLn77rtZtGiR9/H111/PqFGj7AtIJNLVDdfvuXA0hh///O/LHkYJrYjbWwJLlvjtdeX4Iuk6qoQ0Sh01IYWDCanmkUoAJScn06pVK7vDEIkOBw7Av/8NwJ6LRvv3tWNjeYOfmPf/8Q//vrYcUyRdR5WQRqljJqT180jVQxrRysvLGTNmDMnJybRr147HH3+8weerqqq4++676dChA82aNWPw4MEsXrzY+/k5c+aQlpbG+++/T48ePUhOTubiiy9m+/bt3jaLFy9m0KBBNGvWjLS0NM466yy2bNkCNBxqevDBB3nppZd46623cDgcOBwOFi9ezPnnn8+kSZMaxLVz507i4+Mb9AqIyHF8/DGUl0P79lT0Guj3l38Nc9ie11+Pqmp7XUf9R3vZRylvUdPhc0hBQ/YnyjDM/VaDLSnJnHLRRL/5zW/4+OOPeeutt2jTpg333nsvBQUF3ovbpEmT+Oabb5g/fz7t27fnjTfe4OKLL2b16tV07doVgIqKCh577DH+9re/ERMTw3XXXcfdd9/NK6+8Qm1tLaNGjWL8+PH8/e9/p7q6muXLl+NoJMa7776btWvXUlZWxuzZswFo2bIlN910E5MmTeLxxx/H5XIBMHfuXDp06MD5559/gl8wkSjy9tvmx0svtXSdaKqPOJ/a5FRid+yAL76AM844sRfUdTT6rqNGGCgtLTUAo7S01O5QIkJ+vmGU0NIwwDC+/vrIBk8+aX7u6quN/Hzzbn5+0MMMC5WVlcY333xjVFZWHnxy/37zixbs2/79TY573759Rnx8vPHqq696n9u1a5eRmJho3H777caWLVsMp9NpbNu2rcFxF1xwgZGbm2sYhmHMnj3bAIwNGzZ4Pz9jxgwjIyPD+3qAsXjx4kZjmDp1qtGnTx/v47FjxxpXXHHFEV/fFi1aGAsWLPA+17t3b+PBBx9s8ns9/PWO+H7VifTrTKS/PzkGj8cwOnY0rxP/+pffr+v1r7f7wqvNO/ffb+l4XUd1HTUMw9CQfRRyVFfRit3mg2PNIVUPacT6/vvvqa6uZvDgwd7nWrZsSbdu3QBYvXo1brebU089leTkZO/t448/5vtDNk1ISkqic/1mCkC7du3YsWOH9/Wuv/56hg8fzsiRI/nzn//cYBiqKRISEvjlL3/JrFnmzjIFBQV8/fXXXH/99b6+dZHo8/XXsGULJCTABRcE7DSlZ11i3nn33YCdI5ToOupfGrKPQrG7zR90T2wcMS1aHNmgfrF8Lfvkm6Qk2L/fnvP6yf79+3E6neTn5+N0Oht8Ljk52Xs/7rCFtR0OB4ZheB/Pnj2b2267jYULF7JgwQLuv/9+PvzwQ86wMJx300030bdvX3744Qdmz57N+eefT8dDt7sVkWOrK2biggv8ep04XNlZI8w7BQWwfTu0a+f7i+k66n0cLddRJaRRKG63OX+0tkVr4mMa6SSvn1daXBzEqCKIwwHNmtkdxTF17tyZuLg4li1bRlZWFgB79uzh22+/5dxzz6Vfv3643W527NjB0KFDT+hc/fr1o1+/fuTm5jJkyBDmzZvX6IU0Pj4edyPFEKeffjoDBgzghRdeYN68eTz99NMnFI9I1KmfP3rZZQE9TW3LNjBwoDmH9L334IYbfH8xXUcbiIbrqIbso5Cz1Byur01Lb7xBfUJaUUFMZXmQopJgSk5O5sYbb+Q3v/kNH330kXf4JqbuH5RTTz2VX/ziF4wZM4bXX3+dTZs2sXz5cvLy8njnnXeadI5NmzaRm5vL0qVL2bJlCx988AHfffcdPXr0aLR9p06d+Oqrr1i/fj0lJSXU1NR4P3fTTTfxyCOPYBgGP/nJT078CyASLfbsgWXLzPuXXBL48116qfmxideJcKbrqH8pIY1CsWVmQupOadl4g+RkSEw02+5SL2mkevTRRxk6dCgjR44kJyeHs88+m+zsbO/nZ8+ezZgxY7jrrrvo1q0bo0aN4osvvvD2BBxPUlIS69at46c//Smnnnoqv/rVr5g4cSI333xzo+3Hjx9Pt27dGDBgAK1bt+bTTz/1fu6aa64hNjaWa665hoSEhBN74yFixowZdOrUiYSEBAYPHszy5cuP2nbOnDneZVzqb5HydZAAW7wYPB7o1g2a+Lvrq7VrYd3J5rB97QeLKNwY+bs26TrqRz6VWAWZqkP9a8tvZxgGGLvPv/LojeoqMtfOWaoq+2M4VrWh+M+mTZuMmJgYI/8EfxBDpcp+/vz5Rnx8vDFr1ixjzZo1xvjx4420tDSjuLi40fazZ882UlJSjO3bt3tvRUVFls6p62iUmjjRrB6fONH7lL+r7LdsMYykJPM1Y6g1dpNmGGAMdS0ztmw5/vG6jgZHqF9HLfeQLlmyhJEjR9K+fXscDgdvvvnmMdu//vrrXHjhhbRu3ZqUlBSGDBnC+++/70vuLH5y3B5S8A7bx6mHVGxUU1NDUVER999/P2eccQb9+/e3OyS/eOKJJxg/fjzjxo2jZ8+ezJw5k6SkJG8VbGMcDgdt27b13jIaW0P4EFVVVZSVlTW4SRT6z3/MjwGsrs/KMntH8/Phi3wnjmHDADirahElJQE7rTRRuFxHLSek5eXl9OnThxkzZjSp/ZIlS7jwwgt59913yc/P57zzzmPkyJGsXLnScrDiH87SXQDUph5ju7G6P3axu5WQin0+/fRT2rVrxxdffMHMmTPtDscvqquryc/PJycnx/tcTEwMOTk5LF269KjH7d+/n44dO5KZmckVV1zBmjVrjnmevLw8UlNTvbfMzEy/vQcJEz/8AOvXY8TE8GWLYRQUmAXwa9f6/1RZWdC/v3lL+6mZ/F7AItauxXveggIoLPT/ueXYwuU6arnKfsSIEYwYMaLJ7adPn97g8cMPP8xbb73F22+/Tb9+/ayeXvwgtq6oyZ16jB7SNm0AiKtbIkrEDsOGDWuw/EkkKCkpwe12H9HDmZGRwbp16xo9plu3bsyaNYvevXtTWlrKY489xplnnsmaNWs4qX6r38Pk5uYyefJk7+OysjIlpVFm16uLaAV84clm8HkNl/hLSoL0o9S1nrC63tiz+JQW1x2gioPzFZOSzIQ4wNNZ5RDhch0N+rJPHo+Hffv20bLl0ZOhqqoqqqqqvI811ORf9UP2tU0YslcPqYj9hgwZwpAhQ7yPzzzzTHr06MFzzz3HQw891OgxLpfLu02gRKmPzH3Km428gPwHG34qPT2ASWH37tCuHYnbt/PVs5+xf5C5PeXatXDddVBSooRUjhT0hPSxxx5j//79XH311Udtk5eXx7Rp04IYVXRxljWhh7R+Dql6SJskHP77lND4PqWnp+N0Oik+bJ3f4uJi2ja2c1oj4uLi6NevHxs2bAhEiBIJDIPmy8z5o/GX5NArmNMGHQ6zl3TuXE7dugh+3bT90kPh91OOL1Dfp6Au+zRv3jymTZvGq6++Spu6IeHG5ObmUlpa6r1t3bo1iFFGvvoh+2P2kNZ9f9RDemz1O2xUVFTYHIk0Rf336fCdUYIpPj6e7OxsFi1a5H3O4/GwaNGiBr2gx+J2u1m9ejXtTmQnHIls69YRX7KdA7jY3/vM4J+/vojqkJ/zo9F1NLwE6joatB7S+fPnc9NNN/GPf/yjwWT+xmioKbDqe0ibUtSkKvtjczqdpKWlefcdTkpKwuFw2ByVHM4wDCoqKtixYwdpaWlHbOMXbJMnT2bs2LEMGDCAQYMGMX36dMrLyxk3bhwAY8aMoUOHDuTl5QHwu9/9jjPOOIMuXbqwd+9eHn30UbZs2cJNN91k59uQULZkCQCfcSZpCYnBP//5db2iK1aYW4AeslXm4XQdDQ+Bvo4GJSH9+9//zg033MD8+fO5tH4XB7GHYTStqKl+DukeDdkfT/0wa/3FVEJXWlpak4fFA2n06NHs3LmTKVOmUFRURN++fVm4cKG30KmwsNC72wuY2xGOHz+eoqIiWrRoQXZ2Np999hk9e/a06y1ICCosxLvMUqe3/kdL4H8MZaQdwWRlmbfCQnOnqOMsO6XraPgI1HXUckK6f//+BvOWNm3axKpVq2jZsiVZWVnk5uaybds2Xn75ZcAcph87dix//vOfGTx4MEVFRQAkJiaSmprqp7chTbZvHw63uXtGk4bsS3cTSw1g3xBnqHM4HLRr1442bdo02KZNQktcXJztPaOHmjRpEpMmTWr0c4sXL27w+Mknn+TJJ58MQlQSrgoLoUcPqB/13sQntARWuM5mXKCq6Y9n6FB45RX43/+Om5DqOhoeAnkdtZyQrlixgvPOO8/7uH5ZkbFjxzJnzhy2b99O4SELjT3//PPU1tYyceJEJk6c6H2+vr0E2W6zd7SSBIxjDeO0bAlOJ7jdtGEH0CE48YUxp9MZUgmPiESPkhIzGZ07F3q32EqnS7dgOJ08s3IImXZVtJ999sGEtIl0HY1elhPS461ndXiSefh/+mKzuoR0N8foHQWIiTHXBSkuphW7UEIqIhL6evSA09d/AoCjXz8yexx97mbADR1qfvz8c6jRSJscW1Cr7CUE1CWkuzhGQVO9VmabdLT3m4hI2KjvkTz7bHvj6NHDHG2rqADtzijHoYQ02jS1hxS823goIRURCSOfmD2k3h5Ku8TEHEyKLQzbS3RSQhptdpn72DcpIa3rITWH7EVEJNQ5y/bA11+bD846y95gQAmpNJkS0mizZ4/5gRbHaYh6SEVEwkyzLz8Dw4BTT/Uu32er+l7aTz4Bj8feWCSkKSGNNmVl5gdSjt+2LiFVD6mISHhIXlnXE2n3cH29/v0hMRF27SJhy3q7o5EQpoQ02uzbBzQxIVVRk4hIWEleVTd/1O6Cpnrx8TBwIABJqz+3ORgJZUpIo01dQrqP5sdvqyF7EZGwEUc1SWtXmA/OtGH/+qM54wwAmikhlWNQQhpt6obsm5SQqqhJRCRs9OYrYqqrzKWWuna1O5yDlJBKEyghjTZWhuzVQyoiEjYGs8y8M2gQOBz2BnOowYMBSPz+a5qx3+ZgJFQpIY02Vobs1UMqIhI2zqCuB7KuRzJktG8PmZk4PB4GsMLuaCREKSGNNlaG7Ot6SFPYh6O6KpBRiYjICfL2kNb1SIaUuiTZmzSLHEYJabSxMmSfmooRY/6IxJaql1REJFQ5S3dzKt+ZDwYNsjeYxighleNQQhptrAzZx8RQm2oO2zv3KiEVEQlVzdYsB+BAx1PNoqZQU5eQDmaZuXC/yGGUkEYTj8daQgrUppnD9rF7VdgkIhKq6ivYy08LweF6gH798MTG0Y4i4osK7Y5GQpAS0mhSXu6926Qhe8Bd10OqIXsRkdCV9LU5fzRkE9LERCq79QW0/JM0TglpNKnrHTWcTg6Q0KRD1EMqIhLiDINm9Qnp6SFWYX+I+mT5wMefU1CA91aoDlNBCWl0qauwdyc1B5q2Rl39HFIlpCIiIeq774gt20MlCVR26W13NEcVM8RMlvcsXEZ2Nt5bjx5KSkUJaXSp6yF1N2vacD0c7CF1asheRCQ0LTN7RwvoD3FxNgdzdOmX1VXax+VT8Hk1+fkwdy5UVECJ+jyinhLSaFKXkHqaNa2gCaA2TT2kIiIhrS4h/ZzQHa4H4JRToEULYmqq6Rf3Nf37m72jIqCENLo0GLJvGu8cUvWQioiEpvx8AFYwwOZAjsPhgAF1Ma7Qjk3SkBLSaOLtIW36kL1bRU0iIqHL7YYvvwTqhuxDnRJSOQolpNGkvofUypC9ippERELX+vVQWYk7KZnv6Gp3NMenhFSOQglpNKnvIbUyZO9NSDVkLyIScgoKAKg8tS9GOPxJz842P379NRw4YG8sElLC4KdX/MaXKvv6rUPLy6CmJiBhiYiIj1auBKCiWz+bA2mirCxITzf/nqxebXc0EkKUkEaBwkLzn+gd35tD9jsqm95D6k5pgad+zdLduwMRnoiI+Kquh7SiexjMHwUVNslRKSGNcIWF5rIa2dnw7gKzh/RvbzUnKcn8J/W4nE72kmbe36VhexGRkGEY3h7Syu5h0kMKSkilUUpII1xJibno8Ny5cMX5ZkI68Z4U1q41R06aYhfmsL16SEVEQsimTVBaCvHxVJ7c0+5omk4JqTTCckK6ZMkSRo4cSfv27XE4HLz55pvHPWbx4sX0798fl8tFly5dmDNnjg+hyono0QNaxJhD9lm9mjc5GQXYTUvzjnpIRURCR91wPb17h/QOTUeoL2xaswZHZYW9sUjIsJyQlpeX06dPH2bMmNGk9ps2beLSSy/lvPPOY9WqVdxxxx3cdNNNvP/++5aDlRNUV9RE86bPIYVDekiVkIqIhI76hLRfGA3XA3ToABkZ4HaT9N2XdkcjISLW6gEjRoxgxIgRTW4/c+ZMTj75ZB5//HEAevTowSeffMKTTz7J8OHDrZ5eTkR9QprS9Cp7UEIqIhKS6uaP0j9MCprq1Rc2vfMOSd+sAIbYHZGEgIDPIV26dCk5OTkNnhs+fDhLly496jFVVVWUlZU1uIkf1H8dfewhLVqzi4IC85/ywkJ/ByciIk1mGOHbQwreeaRmQioShIS0qKiIjIyMBs9lZGRQVlZGZWVlo8fk5eWRmprqvWVmZgY6zOjgw5B9ejqUxZkJ6dsv7SI725z+06OHklIREdts3w47doDTac4hDTf1Cem6fJsDkVARklX2ubm5lJaWem9bt261O6TwZxg+DdlnZcGkKWZCetX5u8nPNyv2KyrMCn4REbFBfe9ojx6QmGhvLL6oK2xK2LSWZuy3ORgJBQFPSNu2bUtxcXGD54qLi0lJSSHxKL9ELpeLlJSUBjc5MTEHKsDjMR9YHLJv2cWssm/h3kX//ub1T0ROzIwZM+jUqRMJCQkMHjyY5cuXN+m4+fPn43A4GDVqVGADlJC2979mQrqrYz8KCmDtWpsDsqpdO2jfHofHQz9W2h2NhICAJ6RDhgxh0aJFDZ778MMPGTJEk5iDKaa8rnfU4YBmzawd3EpFTSL+tGDBAiZPnszUqVMpKCigT58+DB8+nB07dhzzuM2bN3P33XczdOjQIEUqoaiwEJb82UziHnqnP9nZcN11NH3Dk1BR10vanwKbA5FQYDkh3b9/P6tWrWLVqlWAuazTqlWrKKybUJibm8uYMWO87X/961+zceNG/u///o9169bxzDPP8Oqrr3LnnXf65x1IkzgrDpk/6nBYO1gJqYhfPfHEE4wfP55x48bRs2dPZs6cSVJSErNmzTrqMW63m1/84hdMmzaNU045JYjRSqgpKYE+bjOJu+X5fuTnQ34+ljY8CQl1xVjqIRXwISFdsWIF/fr1o1/dD9LkyZPp168fU6ZMAWD79u3e5BTg5JNP5p133uHDDz+kT58+PP744/z1r3/Vkk9BFlPuW4U90DAhNQz/BSUShaqrq8nPz2+w+khMTAw5OTnHXH3kd7/7HW3atOHGG29s0nm0Wknkcu7dRUfMv7OnXt2X/v3NlZ/CKhkF73JVSkgFfFiHdNiwYRjHSEoa24Vp2LBhrFypHzg7Oct9WxQfOJiQVleb1UxYHPIXEa+SkhLcbnejq4+sW7eu0WM++eQTXnzxRe/IVFPk5eUxbdq0EwlVQlTSevPv6YHMLiSkptoczQmo69jqxRpWVx0AEuyNR2wVklX24n/eIXtfCsSaNTu4LZ2G7UWCat++ffzyl7/khRdeIN3CBEGtVhK5ktaZw/WV3cJw/dFDZWZSm9qSWNwkfv+13dGIzSz3kEp4itl/AkP2DofZS1pUVJeQhtu4kEjoSE9Px+l0Nrr6SNu2bY9o//3337N582ZGjhzpfc5Tt2JGbGws69evp3Pnzkcc53K5cLlcfo5eQkFiXQ9pRff+tLA5lhPicFDRrR8pyxeRuG4lMMDuiMRG6iGNEg2KmnyhwiYRv4iPjyc7O7vB6iMej4dFixY1uvpI9+7dWb16tbeYdNWqVVx++eWcd955rFq1ShuHRKH6HtKKcO8h5WAvb/00BIle6iGNEjEnMmQPDRPSlv6JSSRaTZ48mbFjxzJgwAAGDRrE9OnTKS8vZ9y4cQCMGTOGDh06kJeXR0JCAqeddlqD49PS0gCOeF6iwL59uAq/A6Cye/gnpBXdzcKmRCWkUU8JaZRwnsiQPTRMSLv6JyaRaDV69Gh27tzJlClTKCoqom/fvixcuNBb6FRYWEhMjAawpBFffonDMNjKSdS2bGN3NCesvpc36dsvwe02t0KVqKSENEr4bch+927/BCQS5SZNmsSkSZMa/dzixYuPeWxjq5lIlKjbMnQl/TjJ5lD8oSqrK+Uk0ayqAtavh5497Q5JbKJ/waPECQ/Zt6wbp9ccUhER+9QtoVhAf5sD8ROnky/pY97X8pBRTQlplPDrkL2IiNjjkB7SSOF9L0pIo5oS0igRoyp7EZHwduAArFkDRFAPKYe8FyWkUU0JaZQ4oYXxQQmpiIjdvv4a3G5q0tL5ISJmkJoa9JBqe+qopYQ0SmjIXkQkzBUcukOTw95Y/GgNvTCcsbBnDxQW2h2O2EQJaZTw25C9quxFROyx8uAOTZGkGheVnXuZD+qSbok+SkijxAkP2ddX2e/ZY64VJyIiwVXfQxoBC+Ifrn7HJs0jjV5KSKOCceJD9vUJqWHg3LfXL1GJiEgT1dbCV18BkddDCoe8JyWkUUsJaRRwUYXDXWs+8DUhjY/3HhtbqnmkIiJBtW6dWWXfvDlVJ3W2Oxq/q1APadRTQhoFmrPv4IPkZN9fqG4eqRJSEZEgq59b2bcvROC2spWn9gGHA7Ztg5077Q5HbBB5P9VyhBTqhuubNTuxfYLrElJnqQqbRESCqr7nsH/kDdcDeJo1hy5dzAfqJY1KSkijgLeH1Nfh+np180jVQyoiEmT1PaT9Iq+gyav+vanSPiopIY0C3oTU1wr7ehqyFxEJPo8n4ntIgYPvTT2kUUkJaRTwDtmfaA+pd8heCamISNBs3Aj79oHLBd272x1N4PRTYVM0U0IaBfw2ZF/fQ7pXCamISNDUD2H37g1xcfbGEiBr18JXzrqE9LvvWPW/fdq0KcooIY0CGrIXEQlfpYvNhHRnVn8KCszkLVKkp0NSElx3HfTJac0PdADg1nNW0aOHdhKNJkpIo4Dfh+zLVGUvIhIMhYXwxXPmEPb9/+xHdraZvCUlmclcuMvKMhPs/HzzlnyOOY/0yV+upKICSkpsDlCCJtbuACTwVGUvIhKeSnYa9PGYPaR3vtyfm+u2fE9PN5O5SJCVdch7GdYPlrzNyXs1jzTaKCGNAn4fstccUhGRoIjbsY3WlGA4nXT/2emQYHdEAVZX2JS4XglptNGQfRRQlb2ISHhKWmf2jh44uSckRHo2ysGE9Ps1xFNlczASTEpIo4C/q+ydBypwceAEoxIRkeOpT0grukfw+qOHysqCFi1wuGvpxRq7o5EgUkIaBfw2ZJ+a6t16tBXqJRURCbT6oeuKbhG8Q9OhHA7vAvn90Y5N0cSnhHTGjBl06tSJhIQEBg8ezPLly4/Zfvr06XTr1o3ExEQyMzO58847OXBAPWzB4rche4fDW9jUElXai4gEmreHtEeU9JCCd9i+H5pHGk0sJ6QLFixg8uTJTJ06lYKCAvr06cPw4cPZsWNHo+3nzZvHPffcw9SpU1m7di0vvvgiCxYs4N577z3h4KVp/DZkD96EVD2kIiIBtnMn8cU/AFDZtY/NwQSREtKoZDkhfeKJJxg/fjzjxo2jZ8+ezJw5k6SkJGbNmtVo+88++4yzzjqLa6+9lk6dOnHRRRdxzTXXHLdXVfzHb0P24J1HqoRURCTA6rbQ/JaueJL9cP0OF3UJaR++BLfb5mAkWCwlpNXV1eTn55OTk3PwBWJiyMnJYenSpY0ec+aZZ5Kfn+9NQDdu3Mi7777LJZdcctTzVFVVUVZW1uAmvvPbkD0oIRURCZa6LUMLiKLheoBTT8WdkEQzKtj2328pKMB7085NkcvSOqQlJSW43W4yMjIaPJ+RkcG6desaPebaa6+lpKSEs88+G8MwqK2t5de//vUxh+zz8vKYNm2aldDkGPw6ZK+EVEQkOOp6SFfSj1NtDiWonE5qe/bBWbCUv/92JX+nh/dTSUnmzk6RsimAHBTwKvvFixfz8MMP88wzz1BQUMDrr7/OO++8w0MPPXTUY3JzcyktLfXetm7dGugwI5ajppqE+rXc/JiQqqhJRCTAorWHFHCdYQ7bP/nLld5tRefORduJRjBLPaTp6ek4nU6Ki4sbPF9cXEzbtm0bPeaBBx7gl7/8JTfddBMAp59+OuXl5fzqV7/ivvvuIybmyJzY5XLhcrmshCZHEVO+7+ADFTWJiISHsjLYsAEwe0ijTt080owfV5IRffl4VLLUQxofH092djaLFi3yPufxeFi0aBFDhgxp9JiKioojkk5n3VqWhmFYjVcsclaYCanHlQBxcSf+ghqyFxEJvFWrAKjOyGQX6fbGYoe6hJSVK0G5QlSwvJf95MmTGTt2LAMGDGDQoEFMnz6d8vJyxo0bB8CYMWPo0KEDeXl5AIwcOZInnniCfv36MXjwYDZs2MADDzzAyJEjvYmpBE59D6k7qbl/5mcoIRURCbyCQ3ZoKj5O20h02mkQGwu7d5uVTB072h2RBJjlhHT06NHs3LmTKVOmUFRURN++fVm4cKG30KmwsLBBj+j999+Pw+Hg/vvvZ9u2bbRu3ZqRI0fyhz/8wX/vQo7KWW5W2LubpeCH/tEGCWmlP15PRESOVFfQVNG9H3xscyx2cLmgVy/48kvza6GENOJZTkgBJk2axKRJkxr93OLFixueIDaWqVOnMnXqVF9OJSfIO2TfzA/zR6FBQvqDf15RREQOV9dDWhkte9g3pl+/gwnpqFF2RyMBpr3sI9yhQ/Z+cWiVveb1iIj4X2WlubYRUbSHfWMOnUcqEU8JaYSrH7L3NPPTLh91VfaxuInZrw0LRET8bvVqc4ei1q2padPB7mjsU5+Q1vUWS2RTQhrhYuqG7N3+GrJPTMTjSgQgtlSFTSIiflefgPXrBw6HvbHYqU8f8+O2bbBzp72xSMApIY1wzrohe4+/huyB2jRz2F4JqYjvZsyYQadOnUhISGDw4MHe7ZUb8/rrrzNgwADS0tJo1qwZffv25W9/+1sQo5Wgqh+i7h/F80cBUlKga1fzvobtI54S0gh3aJW9v9SmKiEVORELFixg8uTJTJ06lYKCAvr06cPw4cPZsWNHo+1btmzJfffdx9KlS/nqq68YN24c48aN4/333w9y5BIU9T2k0Z6QguaRRhElpBHO70P2HExInXuVkIr44oknnmD8+PGMGzeOnj17MnPmTJKSkpg1a1aj7YcNG8ZPfvITevToQefOnbn99tvp3bs3n3zySZAjl4CrqYGvvjLv94vigqZ6SkijhhLSCBeIIXt3fQ9pmfazF7Gqurqa/Px8cnJyvM/FxMSQk5PD0qVLj3u8YRgsWrSI9evXc8455xy1XVVVFWVlZQ1uEgbWroXqanOr51NOsTsa+6mwKWooIY1wgRmyr6u015C9iGUlJSW43W7vZiL1MjIyKCoqOupxpaWlJCcnEx8fz6WXXspTTz3FhRdeeNT2eXl5pKamem+ZmZl+ew8SQIcO18foT7Q3If3uO+8yhhKZ9NMe4WL8vTA+B3tInUpIRYKmefPmrFq1ii+++II//OEPTJ48+YiNSA6Vm5tLaWmp97Z169bgBSu+y883P2r+qKlNG+hgLn2V+O2XNgcjgeTTTk0SPpz+XhifQ4qaNIdUxLL09HScTifFxQ03KC8uLqZt27ZHPS4mJoYuXboA0LdvX9auXUteXh7Dhg1rtL3L5cLlcvktbgmOqqUFuIDNLfuzu8C7Pn5069cPtm0jaf1K4Gy7o5EAUQ9phItRlb1ISImPjyc7O5tFixZ5n/N4PCxatIghQ4Y0+XU8Hg9VVVWBCFFsUrjJTW3+KgBGPJBNdjZcdx0kJUF6ur2x2apu2D5pneaRRjL1kEY4v+9lzyFV9kpIRXwyefJkxo4dy4ABAxg0aBDTp0+nvLyccePGATBmzBg6dOhAXl4eYM4HHTBgAJ07d6aqqop3332Xv/3tbzz77LN2vg3xs/3568miglpXEq/871Rwms+np0NWlr2x2aouIU1cr0r7SKaENML5fS97wF1f1KQqexGfjB49mp07dzJlyhSKioro27cvCxcu9BY6FRYWEnNIQUt5eTm33HILP/zwA4mJiXTv3p25c+cyevRou96CBEB9D+CB7n3pP9BpczQhpD4h/X4N8VQBmooSiZSQRrLaWpwHKoAADdlrDqmIzyZNmsSkSZMa/dzhxUq///3v+f3vfx+EqMROiXUJaUX3bJJtjiWkdOwILVrg2LOHXqwBVPAViTSHNJLt3++9G5Ah+/IycxFnERE5YUnehFQJVwMOh7eXtB8ato9USkgj2T5zuL6aOIx4/w1xuFNa4MFhPtitYXsRkRPm8XgT0soeSkiPUJeQ9keFTZFKCWkkq9uZpQz/DdcD4HSylzTz/i4N24uInLDvv8dZvo9KEqg8uafd0YQe9ZBGPM0hjWR1PaT78N9wfb1dtKIle9RDKiLiD3U7NH1Fb+Ji9af5CHUJaR++5K01brxLEKBVCCKFfuojWQAT0t2YlfbqIRUR8YO6HZoK6M9gm0MJSd264UlIpNmBCn435jvW0937qaQkcwMBJaXhTUP2kSxQQ/aYPaTmHSWkIiInrK6HNJ9smwMJUU4nMX37APDeH1aSn2/m8HPnQkUFlJTYHJ+cMCWkkSzAQ/bmHSWkIiInxDC8CWmBljQ6urph+5P3FNC/P/TvDz162ByT+I0S0kimhFREJPRt2QJ79uCJjWMNveyOJnTVJaSsVGFTJFJCGsk0ZC8iEvrq5o8e6HI61dqF6Oj61/UeFxSYvcoSUZSQRrIgFDXt2bibggLz+lBY6PfTiIhEvgItiN8kp50GcXGwZw9s2mR3NOJnSkgjWYAS0vR02B9v9pB++dEusrMhO9ucy6OkVETEIiWkTeNyQR+zsIkVK+yNRfxOCWkkC9CQfVYW/Gl2OgCDO5eo0lFExFeG4R2yV0LaBAMGmB+/+MLeOMTvlJBGsgAO2Wec1gaAxH07VekoIuKrbdtg505wOqns2tvuaELfwIHmR/WQRhyfEtIZM2bQqVMnEhISGDx4MMuXLz9m+7179zJx4kTatWuHy+Xi1FNP5d133/UpYLEggAkpbcyElJIScLv9//oiItGgrneUHj0wEhLtjSUc1PeQ5ueDx2NvLOJXlhPSBQsWMHnyZKZOnUpBQQF9+vRh+PDh7Nixo9H21dXVXHjhhWzevJnXXnuN9evX88ILL9ChQ4cTDl6OI4BV9qSbQ/Z4PKq0FxHxVf3Qc33Pnxxbz56QmGh2uHz7rd3RiB9ZTkifeOIJxo8fz7hx4+jZsyczZ84kKSmJWbNmNdp+1qxZ7N69mzfffJOzzjqLTp06ce6559KnfmKyBE4ge0hjY6FV3dJPR/lnREREjkMJqTWxsQfXI9U80ohiKSGtrq4mPz+fnJycgy8QE0NOTg5Lly5t9Jh//etfDBkyhIkTJ5KRkcFpp53Gww8/jPsYw7xVVVWUlZU1uIkPApmQAmRkmB+LiwPz+iIikcwwDs6FVELadJpHGpEsJaQlJSW43W4y6hOROhkZGRQVFTV6zMaNG3nttddwu928++67PPDAAzz++OP8/ve/P+p58vLySE1N9d4yMzOthCn1AjlkDwfnkaqHVETEuo0bYfduiI+H3ipoajJV2kekgFfZezwe2rRpw/PPP092djajR4/mvvvuY+bMmUc9Jjc3l9LSUu9t69atgQ4z8ng8sH8/EMAeUiWkIiK+q0+o+vY1k1Jpmvoe0pUrobbW3ljEb2KtNE5PT8fpdFJ82BBtcXExbdu2bfSYdu3aERcXh9Pp9D7Xo0cPioqKqK6uJr6RX0KXy4XLpe3TTkh5ufduwIfslZCKiFin+aO+6doVUlKgrIzEjWsA1aREAks9pPHx8WRnZ7No0SLvcx6Ph0WLFjFkyJBGjznrrLPYsGEDnkOWZ/j2229p165do8mo+EndcL0RE0MlAVpKpL6HVHNIRUSsU0Lqm5gYc3tAIOkbzSONFJaH7CdPnswLL7zASy+9xNq1a5kwYQLl5eWMGzcOgDFjxpCbm+ttP2HCBHbv3s3tt9/Ot99+yzvvvMPDDz/MxIkT/fcu5Eh1BU3u5FTAEZhzaMheRMQ3tbUH1yBVQmpd3TzSpG80jzRSWBqyBxg9ejQ7d+5kypQpFBUV0bdvXxYuXOgtdCosLCQm5mCem5mZyfvvv8+dd95J79696dChA7fffju//e1v/fcu5Eh1PaTuZikQqEUKlJCKiPhm7Vpzv+XkZOjWze5owk9dEq8e0shhOSEFmDRpEpMmTWr0c4sXLz7iuSFDhvD555/7cirxVV1C6mkWoAp70BxSERFf1Q/XZ2fDITUW0kR1PaSJ331FPFWA6k7Cnfayj1SH9pAGiuaQioj4pj4hHTTI3jjCVadO0KoVMbU19OYru6MRP1BCGqnqE9LkICSkFRXEVJYfu62IiBykgqYT43B4e0kHonmkkUAJaaQKRg9pcrK5pzAQu1vD9iIiTVJVBV/V9eopIfVd3dduAJpHGgl8mkMqYaB+DmlSgNYgBfM/1DZtYMsWYncXAycH7lwiIpHiyy+hpgZ3y3S+3NURdptPr11rb1hhpy4hHcwyqmwORU6cEtJIFYweUvAmpHHqIRURaZq64fr/7B3IxQMaLsuXlATp6XYEFYYGDwagB2v5al8pkGpvPHJClJBGqmDMIQXvPFIN2YuINNHy5QB87hnI3LnQo8fBT6WnQ1aWTXGFm4wMqjqcjGvbJpqtWQ7nXmh3RHIClJBGqmAs+wTepZ9i9yghFRFpkqVLAficMxjZA/r3tzmeMFZ+2hlmQrr6c0AJaThTUVOkCuaQPRC3S0s/iYgc165d8N13gJmQyokp721uW24mpBLOlJBGqmAlpHU9pHG7igJ7HhGRSFC3ScyBTt3ZSwubgwl/5aebSX2zrz8Hw7A5GjkRSkgjVbAS0vbtAYgt2R7Y84iIRIK64fr9dT17cmIqT+1DJQnElu729jxLeFJCGqnq55AGuqipLiGN2/ljYM8jIhIJ6hLS+p49OTFGXDz5ZJsPtEV5WFNCGqmC1UParh0A8Tt/BDRcItJUM2bMoFOnTiQkJDB48GCW11VeN+aFF15g6NChtGjRghYtWpCTk3PM9hKi3G5vhX356eoh9RfvXNy6ZF/CkxLSSGQYQU9IY6oqSaEssOcSiRALFixg8uTJTJ06lYKCAvr06cPw4cPZsaPx1SoWL17MNddcw3//+1+WLl1KZmYmF110Edu2bQty5HJCvv4a9u+H5s05cEpPu6OJGN6EVD2kYU0JaSSqqoKaGiAICWlSEqSlAdAeDduLNMUTTzzB+PHjGTduHD179mTmzJkkJSUxa9asRtu/8sor3HLLLfTt25fu3bvz17/+FY/Hw6JFi4IcuZyQ+h68wYPB6bQ3lgiylLre5q++gvJye4MRnykhjURlB3sqPUnJgT9fXS+pElKR46uuriY/P5+cnBzvczExMeTk5LC0iUOOFRUV1NTU0LJly6O2qaqqoqysrMFNbFb//R2i4Xp/+pEOVGecBB6PdxcsCT9KSCNR/R+e5OTg/BdeV9ikhFTk+EpKSnC73WTULZlWLyMjg6Kipi2f9tvf/pb27ds3SGoPl5eXR2pqqveWmZl5QnGLHyghDZjtHc2v6bZ/fk5BARQW2hyQWKaENBLVJ6QpAR6ur1eXkLZDSz+JBNojjzzC/PnzeeONN0hISDhqu9zcXEpLS723rVu3BjFKOUJJycFlic5Qhb2/pKebM8f+vNz8mn7x9OdkZ5vbsSopDS9KSCNRsBNSDdmLNFl6ejpOp5Pi4oa7mxUXF9O2bdtjHvvYY4/xyCOP8MEHH9C7d+9jtnW5XKSkpDS4iX12vH1wQfyCTS1Yu9bmgCJEVhasXQs3zzIT0ktbLmXu3wwqKsz/ASR8KCGNMIWF8P2qfQCUO1OCc9HTkL1Ik8XHx5Odnd2gIKm+QGnIMYZy//SnP/HQQw+xcOFCBgwYEIxQxU8KC2H2r8yEdN7mIWRnw3XXmT176ek2BxcBsrKg2zX9IS6OuN076NN8o90hiQ9i7Q5A/Kew0BymGFVRxivAZ2tSgnPR05C9iCWTJ09m7NixDBgwgEGDBjF9+nTKy8sZN24cAGPGjKFDhw7k5eUB8Mc//pEpU6Ywb948OnXq5J1rmpycTHJyEAoX5YSUlMDA2s8AuOC+IeRfaT6fnm4mU+IHCQkwcCB89hnJK/8HdLY7IrFICWkEKSmBigq4/foymAMDzk8h/9EgXPQOGbIvDeBpRCLF6NGj2blzJ1OmTKGoqIi+ffuycOFCb6FTYWEhMTEHB7CeffZZqqurueqqqxq8ztSpU3nwwQeDGbr4wFFTzRDMgqaOvzibjj1sDihSDR1qJqSr/gdcb3c0YpES0gjUobk5h7RFVgot+gfhhIcM2ZcaBuAIwklFwtukSZOYNGlSo59bvHhxg8ebN28OfEASMElr80nkADVp6cR17253OJFr6FD44x/rekgl3GgOaQRylttT1JTIAZz79gbnnCIiYaI+QSrvNxQc+oc9YM46CxwOEgq/I4OmLaEmoUMJaQSKCXZCmphIbfM0AOJKNI9URORQyQVLANjfb6jNkUS4tDQ4/XQAzuYTe2MRy5SQRiDn/iAnpEBNa3PYPm6nKu1FRLzcbpqtMpOjff3PsTmYKDDUTPqHomH7cKOENAI599eVFjVvHrRzehPSHduCdk4RkZD39dfE7i9lH8lUdu1jdzSRTwlp2FJCGoGc+/ead1q0CNo5a1p3ACBupxJSERGv/5mJ0aecBbGqIw64uoS0D18SUz9aKGFBCWkE8hYWpaUF7ZzVbc11peKLtFebiIjXEnP+6P/Q/NGgaN+eqg6n4MTD7n9/RkEB3pu2Eg1tPiWkM2bMoFOnTiQkJDB48GCWL1/epOPmz5+Pw+Fg1KhRvpxWmkgJqYhICDAMbw/pEjR/NFhqhpjJ/9JH/0d2Nt6b9rcPbZYT0gULFjB58mSmTp1KQUEBffr0Yfjw4ezYseOYx23evJm7776boUP1X2KgeRPSIA7ZKyEVETnMhg1QVIQnLp4vGGh3NFEjeYSZ/N/W73/k50N+Psydi/a3D3GWE9InnniC8ePHM27cOHr27MnMmTNJSkpi1qxZRz3G7Xbzi1/8gmnTpnHKKacc9xxVVVWUlZU1uEnTOPAcXIdUPaQiIvap6x2t6DWIKhJsDiaKnGMmpM3XLKN/9wr69zd7RyW0WUpIq6uryc/PJycn5+ALxMSQk5PD0qVLj3rc7373O9q0acONN97YpPPk5eWRmprqvWVmZloJM6qlUIbDMMwHqalBO29NW/N75CzfB6XaQFREhLodt/Zruafg6twZTjoJqqvhs8/sjkaayFJCWlJSgtvt9u63XC8jI4OiosZ3Rfjkk0948cUXeeGFF5p8ntzcXEpLS723rVu3WgkzqqWx17yTmAguV9DO60lsRgmtzAeapCMiUaywEAryDaoXLgLg6zbn2xxRlHE44IILzPuLFtkbizRZQKvs9+3bxy9/+UteeOEF0tPTm3ycy+UiJSWlwU2axpuQBnG4vl4hWXV3lJCKSHQqLDSHh68dsJ74nT9yABcjHzmTpCSw8GdQTtT5df8EfPSRvXFIk1laFC09PR2n00lxcXGD54uLi2nbtu0R7b///ns2b97MyJEjvc95PB7zxLGxrF+/ns6dO/sStxyF3Qlpf1YqIRWRqFVSYhbP/G3sIngJqgeexaczE0lPh6wsu6OLIvUJ6YoVsHcvkGZjMNIUlnpI4+Pjyc7OZtEhXeAej4dFixYxZMiQI9p3796d1atXs2rVKu/t8ssv57zzzmPVqlWaGxoAdiek5h0lpCIS3boWmn8nU0ZdQP/+SkaD7qST4NRTwePxrgUroc3ythGTJ09m7NixDBgwgEGDBjF9+nTKy8sZN24cAGPGjKFDhw7k5eWRkJDAaaed1uD4tLpE6fDnxT+UkIqI2CsGN8n5i80H52v+qG3OPx++/dacRzr2crujkeOwnJCOHj2anTt3MmXKFIqKiujbty8LFy70FjoVFhYSE6MNoOyihFRExF59WUVs2R5ISYEBA+wOJ3pdcAHMnGnOIx1rdzByPD5trDtp0iQmTZrU6OcW1y1zcTRz5szx5ZTSRKGQkFZtKGRNgfmc5k2JSLS5gLppbeeeq/3r7TRsmPnx66+J3VUMZByrtdhMXZkRxq6END0ddiZ2BCC26AcGZ9doqzYRiUrehLR+6SGxR3o69OkDQPMV/7U5GDkeJaQRxq6ENCsL/ru2LR5XIk48rHpzi7ZqE5Go46iuYijmDk2aPxoC6v4paL5c65GGOiWkEcbOIfusjg5iuphbw/ZK+F5btYlI1Gm26lOSqKSmVQaoeNd+9Qnp5x8Ahr2xyDEpIY0wdiakgLllG8D339tzfhERG6V+9h4AZWcMN3cMEnsNGwYuF66iQnqw1u5o5BiUkEYYJaRhxuOBf/4T7r4b5s6F2lq7IxKRE5BSn5CeNcLmSASApCRvcdMI3rM3FjkmJaQRRglpGKmshMsvh6uugscfh1/+EkaOhAMH7I5MRHyxdSuJ36/BTQxlZ1xkdzRSb4T5z4ES0tCmhDTCtGCPeUcJaei7/XZ45x1ISIDrrzf/k1+4EP7v/+yOTER88Z6Z8CxjMO7UljYHI151CelQ/kdMxX6bg5GjUUIaSWpqSKXMvN+qlT0x1CekGzeCoQnkR/Xee/DCC+Ycs7fegtmz4R//MD/39NOwapWt4YmID+oS0vfQcH1I6dqVqg6n4KKa5l98ZHc0chRKSCNIbNluAAyHw74e0o4dISYGKiqILSmyJ4ZQ5/Ec7AW94w64qG5o75JLYPRoM5H//e9tC09EfFBdbW5RiRLSkONwUFo3pzflUw3bhyolpBEktnQXAO6UFuB02hNEfLx3aybXDxq2b9Srr8LXX5v/NEyZ0vBzDzxgfnz9ddi6NeihiYiPPvsM9u2jpmUbCuhvdzRymLIz6xLSz97T6F2IUkIaQZx1CWltqk3D9fXqhu2VkDbCMOAPfzDv33XXkT3ZvXqZ2w0aBrz8ctDDExEf1Q3Xlw0ZjqE/rSFn34DzOIAL1/YtsG6d3eFII/RbE0Fi94ZYQrr1O3vjCEWffGL2jiYlwaRJjbcZN878OGeO/pMXCRdvvw1A2VmX2ByINMZITGIxw8wH//qXrbFI45SQRpDYveYenW67E9K6LZoSNmkR4iM884z58dprjz7P96qrIDERNmyAr74KWmgi4qP162HtWoiLo/RMzR8NVW8yqu7Om3aGIUehhDSCxIbKkH1dQpq48Rt74wg1xcXmIvgAEyYcvV2zZgcLnd56K/BxiciJqU9wzj8fT/NUW0ORo/sXl5t3Pv+c1e//SEEBFBbaG5McpIQ0gngT0jSbE9KePQFwbd1ALDX2xhJK5s2DmhoYOBD6H6fo4YorzI9KSEVCX31COmqUnVHIMaSnQ2lSe5ZyBgAzLv4X2dlm/4mS0tCghDSC1Bc12T5kf9JJkJyMw11LFzbYG0soWbDA/PjLX3qfKiyEgoKGt8JC4NJLzTVKCwqgSMtniYSs7dvh88/N+5dfbm8sclRZWeasiqxbRwHwpzPeYO5cqKiAkhJ7YxOTEtIIEjJFTQ6Hd9i+Jxq2B2DTJli2zFyj9Wc/A8zEs0cPyM5ueOvRAwoPtIG+fc1jP9JCziIhq75AZvBgaN/e3ljkmLKyoMPEUQCkrPiI007aa2s80pAS0ggSMnNIwTts3wMVNgEHe0eHDYO2bQHzv/KKCpg7F/LzzVuD/9gvuMA8RgmpSOjScH146dbN/K+/tpaUT961Oxo5hBLSCBIy65CCekgPN3+++fHnPz/iUz16mFNK+/f3ftlM559vfqzb/UVEQsyePQd/P5WQho+671Xaf9+wNw5pQAlpBPHu1GR3UROoh/RQ69bBl19CbCxceWXTjxs61Dxm82ZzyF8iyowZM+jUqRMJCQkMHjyY5cuXH7XtmjVr+OlPf0qnTp1wOBxMnz49eIHK0b3xhlmoePrp0L273dFIU/3kJwCkfPouSZTbHIzUU0IaKQyD2FJzL/tQ6iHtzjpwu20Oxmb1w/UXXQStLHxvkpPNeWkAixf7PSyxz4IFC5g8eTJTp06loKCAPn36MHz4cHbs2NFo+4qKCk455RQeeeQR2tZN+ZAQcIyRDwlhAwbAKafgPFDB5WiR/FChhDRSlJXhcNcCIZKQnnwy7oQkEjmAqzCKd2wyjBP7o3X22ebHzz7zX0xiuyeeeILx48czbtw4evbsycyZM0lKSmLWrFmNth84cCCPPvooP//5z3G5XEGOVhq1YwdG3XD9171Ge1fJWKtBodDncMA11wBwDX+3ORipp4Q0Uuwyh+srSMRISLQ5GMDppPLUPgAkrSuwORgbffWVOWTvch1cW9SKM880PyohjRjV1dXk5+eTk5PjfS4mJoacnByWLl3qt/NUVVVRVlbW4Cb+s/u5f+DweFjOQE4f1dm7SsZ115k7A6en2x2hHFNdQnoxC3HWjS6KvZSQRoq6hHQXIdA7Wqeiu7n4e9L6lTZHYqP63tFLL4WUFOvHDxlifvzmG7OAQsJeSUkJbrebjIyMBs9nZGRQ5Mc1Z/Py8khNTfXeMjMz/fbaAvGvm7/bzl9c410lo/62dq25xJCEsF69qOxyOvHUkPbR63ZHIyghjRzFxQDsoI3NgRxU2a0fAInR2kN66HD96NGWDl27tm6h/K2tOZDV1XyyfvFtkSbIzc2ltLTUe9u6davdIUWOrVtJXvUJHhzEX3e1d5WM+puS0fCw++JrAWi5cJ7NkQhArN0BiJ/UJaRFtKWdzaHUa9BDahjmvJ1o8sUXZoV8s2Zs7X0pOw/Lyxuba5aebg73XXfdwedmcybX8x2l731G6ogRAQ1ZAi89PR2n00lx3e9sveLiYr8WLLlcLs03DZS5cwH4H0Np3qaDzcGIr/YM/zkdns4lOX8xbNsGHfS9tJN6SCNF3VBfMRnHaRg8Bzr3opo4Ysv2wJYtdocTfHW9o+UXXE737GZH7MjU2Fyz+u3tDh3+63GDOY80ZumndrwL8bP4+Hiys7NZdMj6sh6Ph0WLFjGkfoqGhC7DgLris9mMszkYORHV7TvxP87GYRjw8st2hxP1fEpIrayf98ILLzB06FBatGhBixYtyMnJOWZ78VFdb0soJaRGXDxr6GU+WBll80jdbu9yT0XDfn7EjkzHmmuWldVw+K/5cDMhTfp6GdTWBvudSABMnjyZF154gZdeeom1a9cyYcIEysvLGTfOTHDGjBlDbm6ut311dTWrVq1i1apVVFdXs23bNlatWsWGDRvsegvR65NPYMMG3EnJvMZVdkcjJ2gWN9TdmWX+syG2sZyQWl0/b/HixVxzzTX897//ZenSpWRmZnLRRRexbdu2Ew5eDhGCCSlAAeawPfn59gYSbP/7H/z4I6SmUjZkONBwRyYrc80OnNKTvaTiPFBhVu1L2Bs9ejSPPfYYU6ZMoW/fvqxatYqFCxd6C50KCwvZvn27t/2PP/5Iv3796NevH9u3b+exxx6jX79+3HTTTXa9hehV1zu658LRlJNsczByov7Bz6hJSIYNG/j2r0soKIDCQrujik6WE1Kr6+e98sor3HLLLfTt25fu3bvz17/+1Ts8JX50yBzSUPIFAwEo+/Bz7zp9UfHL/ve6te1++lOM+BOcxxcTwzLqFshXYVPEmDRpElu2bKGqqoply5YxuH4TBMx/5OfMmeN93KlTJwzDOOK2WBsmBNe+ffDqqwDsuuIGm4ORE5WeDkZSMnMOmGtEf/6rF8nONjsPouLvVIixlJD6Y/28iooKampqaNmy5VHbaP08H4RgD2l6OhQknAWAY/kyBmXXRscve3U1vPaaeb9urbsT9Tln1N1RQipim1dfhYoK6NaN8t6a7xvu6ufsD51jjjRc53qNBc+XUlEBJSU2BxeFLCWk/lg/77e//S3t27dvkNQeTuvn+SAEi5qysuC1tb1wN0uhOftZ/cpq5s4l8n/ZP/wQdu+GjAw47zy/vORS6v74+XHhdBGxwDDg2WfN++PGRd+qIREqKwu6jxkEvXoRU1XJGd+/YndIUSuoVfaPPPII8+fP54033iAhIeGo7bR+nkVVVbB3LxBaCSlAVqcYnGebyVSP3Z/Wb3Ef0cpfMNe023HeaAq+dPplK0HvkP2GDRGezYuEqM8/N+fCu1xw4412RyP+5HDAr34FQOtXnwZU3GQHSwnpiayf99hjj/HII4/wwQcf0Lt372O2dblcpKSkNLjJMdQVlHli49hDC5uDacRZ5rB9NGx/uXV9BcZbbwFw+fxr/LaV4F5acKBTd/OBhu1Fgu/pp82P11yjfUEj0fXXQ/PmJG5ay4V8aHc0UclSQurr+nl/+tOfeOihh1i4cCEDBgzwPVppXN0/CLUt22CE4tKyZ59tflyyJOKX1ah9422SKWdf65N5ZsVgv24lWH665pGK2KKoCP7xD/P+rbfaG4sERkqKORUDuJ0/2xxMdLKcvVhdP++Pf/wjDzzwALNmzaJTp04UFRVRVFTE/v37/fcuol3dElo16aGyR9NhhgyBhATYto2EzevsjiagWr49B4DykdfQP9vh160E9/fWPFIRWzz3HNTUmNey/v3tjkYC5dZbMRwOLuVdXFu+tTuaqGM5IbW6ft6zzz5LdXU1V111Fe3atfPeHnvsMf+9i2hXN8e2JiNEi78SEry9pM0/j+ChkK1bSVn6PgC7Lvf/Di4V9T2ky5ebC++LSOCVl8NTT5n3b7vN3lgksLp0oeysSwBoveApm4OJPj7tZT9p0iQmTZrU6OcOXxdv8+bNvpxCrPjhBwCq25xkcyDHcOGF8J//kLLsQyBCL+pz5uAwDBZzLimZXfz+8pWn9ILkZNi/H9asgePMxRYRP3jxRdi1i5qsU1h9ylVQYD7tj2JFCT07rr2D1E/eIf3NF2HHA9Cmjd0hRY0QnHAoltX1kFa3DdEeUjATUiA5fzFxVNscTAB4PN4dXF4kQBW4TicMGmTe1zxSkcCrroa60by7tv+G7MGxZGfjt2JFCT37Bl3AcgYSU1UJTzxhdzhRRQlpJAj1IXuAPn2gTRucFfsZyv/sjsb//vtf2LwZd7MU/slPA3eeIZpHKhI08+bB1q3UtGrL8zXXM3cu3kJFfxUrSohxOHiIB8z7M2bArl32xhNFlJCGscJCcyvOqo3mkP13lSE8ZB8TAyNHAvAT3rA5mACYOROA3RdfSyVJATnF2rWwId2cR3pg8eeRvduViN1qauD3vwfMYdwqEujRA2+hor+KFSX0/JvLqDi1jzk96s+quA8WJaRhqrDQ3IJzQLYHxzYzIR3/UGZoDyH95CcAjOJNc4g7UmzZAq+/DsDOqyf6/eXT082hweuugzPuNBPShM3rOLP7biWlIoEyezZ8/z20bh2Q32sJZQ6KbrzfvPuXv8CePfaGEyWUkIapkhJzC85/zNhJPDUYDgdvfN4+tIeQLrgAd1IyJ7GNpG9W2B2N/zz9tJlg5+RwoMtpfn/5+v2W8/Phg/x0DtQVTJ1WuVybNokEQmUlTJtm3r/vPjxJyfbGI0G39/wr4fTTobQUHn7Y7nCighLSMNcrxZw/6mjblv6D40I3GQVISKD07EsBMP4+n4ICvLew7enbvx9eeMG8f8cdATtNVtbBYcKEYeY80iFoHqlIQMyYAT/+aP7i/frXdkcjdoiJgT/+0bz/l7+AVgwKOCWkYS7+x83mnZDORA+q/fl1ALRaOJfB2TXeitUePcI0KX3hBfM/6K5dYcSI4JzzDHPY/gxUaS/idzt24HnInDu6edyDFKxxaYmnaHXxxXD++eZqC/fdZ3c0EU8JaZhzFX5n3una1d5AmqjNdcNxt86gDTtZ9+RC8vNh7lxz+kHYDT9XVBz8D/r//s/8jzoY6irtB7MssubiioSA/bfmElNWSj796TxtjJZ4imYOBzz6qHl/3jz47DN744lwSkjDnOuHDeadLv5fiD0g4uJw/vIXAHRe/CL9+5u9o2HpueeguBg6dYKxY4N33tNPx52QRBqlEb8Vq0hQLV9O8qvmesIlU57ii3ynlniKdv37ww03mPd//Wtz9QUJCCWkYc61tS4hDZMeUgBuusn8+K9/mVWs4ai83Ns7uuWX91OwOo6CgiDt3hIbS0WvgQA0+0rzSEX8oqYGJkwA4CXG0PqKM7XEUxRbu/ZgjcOX1/4Rd4tWsHq1loEKICWkYc47ZB8uPaRgdolefDEYhjlZPBw98ggUF7PRcQpdHhoT9N1byvucBZg7X4mIH/zxj1BQQG1KC37LH+2ORmxy6DJ79df1vjnpTCyvG7qfOhU2brQ3yAilhDSMJVFOfMl280E4JaQAd95pfnzxRWJ377A3Fqs2b/bOK7rLeIw5c+OCvntL2aAcAFKW/cdM7EXEd199Bb/7HQA//OYvFNPW5oDELocus1d/mzsXnqu+nqIew6Cigv1XjqHgC3d4rxATgmLtDkB815m64e6WLc1bOLnwQhgwAFasoO2sh4HpdkfUNIYBt98OVVWUDTyfN78YxQN1u7cEU3nvIVSSQOKuIvjmG+jVK7gBiESKykoYM8Ycsr/8cnaP+AX1O0dKdMrKatipYPaaOjhj7Wy+ojcpX37Ka4P+SB73kpRk7ovSunXD9priYZ16SMNYd+oKWsJp/mg9hwPy8gBIf+1Zsthic0BN9Mor5tzXuDh+uPvPgMOWMAxXAv9jqPngP/+xJQaRiHDHHfDll2YWMXOmeW0SOUR9r+nr+Z3YPe1pAP7gnMrSJ8w5/BdffHB4P6yXMbSZEtIwdjqr6+6cbm8gvsrJgQsuIKammie50+5ojm/rVrj1VvP+1KkB2ZXJiv+QU3dHCamIT+bOheefN5PQV16Bdu3sjkhCVP3mJJ0e+CWMHo3DXcsZj/6U9Yu3HzG8H5bLGIYAJaRhrDdf1d3pbW8gJ+KJJzCcsVzJG6Qt+qfd0RxdVRX87Gewd6851eC3v7U7ooMJ6eLF5sLNItJ0y5bB+PHm/QcegIsusjceCQ8OB/z1r+Y0qe3bOemOq+jfq8q7GkPYLmMYApSQhrE+fGneCeeEtHdviq6/B4B2v7+F1e//GHrbiRoGTJoEy5ZRm9KCr6e8SsFXsbbv3rKKvtS0yjC3L/3vf+0NRiQMFBaa15av/7WRmhEj4cABSs+6hILLpnivO3b/XksYSE6GN96A1FRzsfyxY7VJiR+oqClMOcv20Kl+3mW4DtnXqfnt/Xw96y1OK1vNvouv4jz+SzUukpJCZDHqKVPgr3/Fg4PLy+by3uUnez9l5+4tBjHsPfcKWr/+vHlxHD7cnkBEwkBhodl7lVKxncWMII6dFNCPcz5dQPkgZ4O22pVJjqtrV/jHP+DSS2HBArOq6S9/wa66gkigHtIw1ezrZQAcyOwSfhX2h8nq6qLl4tepbZ7GmSxl+/BxvPKy2/55OIYBv/+9eQMm8Cy/mHtJ0Jd4OpbS835i3nnrLf2HLnIMJSVmMrqu3Xl041uq2mYRt/DfLMlPbvA7HQq/1xImLrwQXn7ZHMZ/+mm4+24tw3cC1EMappp9ae6pW957CAk2x+IP7c/pAq/+HUaOpOX7f2d4bBxOXsS2H9HaWrjtNnj2WQC2TXyY52fczM02LPF0LPsGng8pKVBUZA4dnX223SGJhKT4bZtYzMWkbv8WMjNxLf4vp5/S3u6wJNz9/OewezdMnAhPPEHmlnIcPIP6+6zTVyxMNVttLjdR3vtMmyPxo4svhvnzwemk1Tsv8x4jcJbtCX4chYUwbJiZjDoc8NRTFN+QG/w4msCIi4ef1PWSzpljaywioaR+vmhBAaz72xd0vu4Mb88oixfDKafYHaJEiltugRdfBIeD1v98jjcZRcz+MrujCjtKSMNRZSXJqz4BYH/fCOsR++lP4Z//xJ3YjAv5D12vPJ0NT70XnEKnmhqYPt0sEvv0U7Pn8Z//NAuaQtTatbB+6E0AuOfNZ+s3+2yOSMR+9fNFs7MNZmY/T8cx55BYtoMvHX0p+ddSJaPifzfcAPPn44l3cTlv023cEPj2W7ujCitKSMPR4sXEVB1gKydxoHME7tBzxRXseP1TvnN0JWnPNrrcdgnbskdyY/bKwCw4XFUFs2ebxWF33gmlpTB4MKxcebD3McQcut9y95vOYh3dcFaW80i/+aGzOoGITUpKwFWxm01nXMvz3EwiBygdehktv15Ch4EappcAufpqvv3r/9hGexI3fgP9+pmbLWheaZMoIQ1H//43AO9yScTuKtLu4j4krF1F8S8mYzidjOTfrKQ/H1SchfHiLNixw6fX9Q7jfeFm/axP2fHz23BndjT/u12/3sz0nn+ewr9/SsHeU0J2KZiG+y07SL7DXE/xturHKCl22xydSHAdOjxfkG9Q+ddXWEd3On0+H2Jj4dFHSV38Fpk9m9sdqkS4il4DGcAKinqeZ66QP2ECZUOG8+NH6+wOLeQ5DCP0U/eysjJSU1MpLS0lJSXF7nDsVV0N7dvDrl1czHs8nH9xSBXZBMS338KDD2K8+ioOt5lsGQ4HFd37U3HaYCpO7UNC3x5k9G0HbduaXYcxMeZ/peXl5mL2xcWUfPYtz01eT3bt55zJZ6RwcHi7NqM9sXfdAb/6FYWlqfToYV5LDhUyy1A1Zt8+ak/qSGzZHjb9YR4n33uN3RGFnUi/zkTq+6sfnq+ogHNZzO+5n7P5FIDqrj2J/9ssc8RDJAjqfx4rKzzcxl94hHtIoIoaYqm8YSIpD+dCRobdYQbMiVxnlJCGm1dfhdGjqU5vR2LJVr7Id0Z+Qlrnhy+288JZsxlV8w/6seqYbQ2HAxwOHMdYCsndLIUt/a5g0ic/55dzcuh2ejxgJp3XXWduAXforhvp6SGajNb5ccJDtJ85haq2Wbg2fAPNmtkdUliJ9OtMpL6/gi/cTB30Li/0nE7bbz4CwONKoOy2B0j7/d0QH29zhBJtCgsPLlnoKvyOlIfuJrPgX+YTCQlw003m9LAInMushDRaeDyQnQ2rVrF9/BTavzCN/PzQWoYo0Op/0WN3bqf5iv+S9O0qnF9/SenK72lrbKcZFUcc44mNozYtnb2tuvCv9acy/Dd9yLx2KJx+OoXbnOHXG3oUqz4tJ+3sXuaGCbfeWrdIszRVpF9nIuX91V8DEjatJe3Df5Dy+mySd242PxkXB7/6Fdx7rzmSJBICCgrgnuwPmNd5CunfL/M+f2DIeSRMGAejRkHzyJhOckLXGcMHTz/9tNGxY0fD5XIZgwYNMpYtW3bM9q+++qrRrVs3w+VyGaeddprxzjvvWDpfaWmpARilpaW+hBsRtmwxjC25zxoGGLVJycarM3cZYBj5+XZHFhq2bDGM/BUeY+X/9hmr/rPT+GjeduOUhG1GM/YZ4DHM8XvDSEoy2x5xbH7D2+FtwkF+vmGM4B3D+2ZffNHukMJKsK8zuo4e1Njv4OG3lZ/sN1b8/j3jidjfGKvpdfDnHIxdtDBKf3W3YWzebPdbETnCli3m3x7wGMP4yHiXiw03joM/w3FxhpGTYxjTpxvGV18ZRm2t3SH77ESuM5YT0vnz5xvx8fHGrFmzjDVr1hjjx4830tLSjOLi4kbbf/rpp4bT6TT+9Kc/Gd98841x//33G3Fxccbq1aubfM5QvpAGw5YthjHCtcioxGUYYNzG9KMmV3JQpCSaTZWfb17bVl9xn/dCt23CQ8aW72vsDi0sBPM6E4nXUV9/3w7+sTZvMdQaJ1FonMNiYxJ/MV5knJFPP6Oa2AZJqDs2zth71iXGpmkvGYXrygPynkT85fDfj9enbzHu53dGaduuDX6uDTCM5s3NBDU31zBeftkwli83jDDJf07kOmN5yH7w4MEMHDiQp59+GgCPx0NmZia33nor99xzzxHtR48eTXl5Of+uqwwHOOOMM+jbty8zZ85s0jkjZajJspoaWLWKHY+9RMtXZxKLm73nXs7GR18HpzPk5zRKcB06mf5J7uR2zCH7DY4utMq9mRZXXwjduplzmOQIwbzOhPt19NA5cgA7d8KVVzac+hKDm/TECvL/V8FJLSvMAsOKCoo3VbD/h73E7tlJ7N4SyjaW8MV7O7nwtCJalG4ifvsWYmprGj1vVbuO1Ay9gOTLz4dLLoEWLU7ofYjY5dBivK58y2X8mxG8xxl8TnP2N3pMTYvW1LTpQE3rDsRkdaD5qe3N34G0tIO3Fi3M+oGEhIM3lwuczqC8rxO5zljal7G6upr8/Hxycw/uWhMTE0NOTg5Lly5t9JilS5cyefLkBs8NHz6cN99886jnqaqqoqqqyvu4tLQUMN9ok8yZA/PmUV1tUF3Dwf87wOwkB/Cm4fWPG2tjHNmm7qPj8Ocaa3Pc1+HIc9a1cdRUE7erCIdhkABUALtzfkbLuU/RJbHc+1ab+iWRyJeWBsuXw65dgPE7Vi/sQdun7qVN+QZ4+DeUPWy2q2mRjichGU9iEp44V93SYQ6oW0GsviAMOPi5MOTo1ZNmz09vcvv664vF/9EtC5vr6H//C1OnUnPATW21GzxuHB43nmoPZcW1JODGWXdrgZtv8NA80U2MuxY8bjOprAQGwKFnTKy71WsOdAD4GqoxbzidkJkJ3bubG1XU37KywOE4+Hq6AEqYanC9pi1wE+vX38TwX7lZ8MA39CxfTuL3a4jZsJ6aNd+SwQ7YsxP27MSxfhUGDX+vjsfjjMWId2HExWM4nRgOJzidGI4Y82OME2IcdR/rno9xEpuaRPxH7zf5PCdyHbWUkJaUlOB2u8k4bMmCjIwM1q1rfI2toqKiRtsXFRUd9Tx5eXlMmzbtiOczMzOthBt5/vMPaPsPu6OQcLenBCg5brOw9/UyWDDb8mH79u0jNTU1AAGZIvo6Wumn13G7YfNm87ZwoZ9eVCQ8XP1QAF7UXQuVtVBZfvy2h/PheujLddRSQhosubm5DXoDPB4Pu3fvplWrVjgibCH4srIyMjMz2bp1a3RNRwgwfV39L9K/poZhsG/fPtpHSHW2rqNyovR19b9I/5qeyHXUUkKanp6O0+mkuLi4wfPFxcW0bdu20WPatm1rqT2Ay+XC5XI1eC4tLc1KqGEnJSUlIn847aavq/9F8tc0kD2j9XQdDZxI/tm0k76u/hfJX1Nfr6OWtg6Nj48nOzubRYsWeZ/zeDwsWrSIIUOGNHrMkCFDGrQH+PDDD4/aXkQkkuk6KiJyJMtD9pMnT2bs2LEMGDCAQYMGMX36dMrLyxk3bhwAY8aMoUOHDuTl5QFw++23c+655/L4449z6aWXMn/+fFasWMHzzz/v33ciIhImdB0VEWnIckI6evRodu7cyZQpUygqKqJv374sXLjQO+G+sLCQmJiDHa9nnnkm8+bN4/777+fee++la9euvPnmm5x22mn+exdhzOVyMXXq1COG1uTE6Ovqf/qa+o+uo/6ln83A0NfV//Q1Pbqw2DpURERERCKXpTmkIiIiIiL+poRURERERGylhFREREREbKWEVERERERspYRURERERGylhDSEbN68mRtvvJGTTz6ZxMREOnfuzNSpU6murrY7tLAyY8YMOnXqREJCAoMHD2b58uV2hxTW8vLyGDhwIM2bN6dNmzaMGjWK9evX2x2WSKN0HfUPXUf9S9fR41NCGkLWrVuHx+PhueeeY82aNTz55JPMnDmTe++91+7QwsaCBQuYPHkyU6dOpaCggD59+jB8+HB27Nhhd2hh6+OPP2bixIl8/vnnfPjhh9TU1HDRRRdRXl5ud2giR9B19MTpOup/uo4en9YhDXGPPvoozz77LBs3brQ7lLAwePBgBg4cyNNPPw2YWzJmZmZy6623cs8999gcXWTYuXMnbdq04eOPP+acc86xOxyR49J11BpdRwNP19EjqYc0xJWWltKyZUu7wwgL1dXV5Ofnk5OT430uJiaGnJwcli5damNkkaW0tBRAP5cSNnQdbTpdR4ND19EjKSENYRs2bOCpp57i5ptvtjuUsFBSUoLb7fZuv1gvIyODoqIim6KKLB6PhzvuuIOzzjpL21ZKWNB11BpdRwNP19HGKSENgnvuuQeHw3HM27p16xocs23bNi6++GJ+9rOfMX78eJsiF2lo4sSJfP3118yfP9/uUCTK6DoqkULX0cbF2h1ANLjrrru4/vrrj9nmlFNO8d7/8ccfOe+88zjzzDN5/vnnAxxd5EhPT8fpdFJcXNzg+eLiYtq2bWtTVJFj0qRJ/Pvf/2bJkiWcdNJJdocjUUbX0eDQdTSwdB09OiWkQdC6dWtat27dpLbbtm3jvPPOIzs7m9mzZxMTo07spoqPjyc7O5tFixYxatQowBwaWbRoEZMmTbI3uDBmGAa33norb7zxBosXL+bkk0+2OySJQrqOBoeuo4Gh6+jxKSENIdu2bWPYsGF07NiRxx57jJ07d3o/p/9Mm2by5MmMHTuWAQMGMGjQIKZPn055eTnjxo2zO7SwNXHiRObNm8dbb71F8+bNvfPIUlNTSUxMtDk6kYZ0HT1xuo76n66jx6dln0LInDlzjvoLr29T0z399NM8+uijFBUV0bdvX/7yl78wePBgu8MKWw6Ho9HnZ8+efdwhVJFg03XUP3Qd9S9dR49PCamIiIiI2EoTa0RERETEVkpIRURERMRWSkhFRERExFZKSEVERETEVkpIRURERMRWSkhFRERExFZKSEVERETEVkpIRURERMRWSkhFRERExFZKSEVERETEVkpIRURERMRWSkhFRERExFZKSEVEfDBnzhwcDgebN2+2dNz1119Pp06dAhKTiEi4UkIqIhIm3n33XR588EG7wxAR8TuHYRiG3UGIiIQbt9tNTU0NLpcLh8PR5ONqamrweDy4XC7L55w0aRIzZsxAl20RiTSxdgcgIhKOnE4nTqfT8nFxcXEBiEZEJLxpyF5ExAeNzSF95pln6NWrFy6Xi/bt2zNx4kT27t3b4LjD55Bu3rwZh8PBY489xvPPP0/nzp1xuVwMHDiQL774osFxM2bMAMDhcHhvIiKRQD2kIiJ+8OCDDzJt2jRycnKYMGEC69ev59lnn+WLL77g008/PW7P6Lx589i3bx8333wzDoeDP/3pT1x55ZVs3LiRuLg4br75Zn788Uc+/PBD/va3vwXpXYmIBIcSUhGRE7Rz507y8vK46KKLeO+994iJMQefunfvzqRJk5g7dy7jxo075msUFhby3Xff0aJFCwC6devGFVdcwfvvv89ll13GkCFDOPXUU/nwww+57rrrAv6eRESCSUP2IiIn6D//+Q/V1dXccccd3mQUYPz48aSkpPDOO+8c9zVGjx7tTUYBhg4dCsDGjRv9H7CISIhRQioicoK2bNkCmL2ah4qPj+eUU07xfv5YsrKyGjyuT0737NnjpyhFREKXElIRkRBwtIp9LfEkItFACamIyAnq2LEjAOvXr2/wfHV1NZs2bfJ+/kSpql5EIpUSUhGRE5STk0N8fDx/+ctfGvRovvjii5SWlnLppZf65TzNmjUDOGIpKRGRcKcqexGRE9S6dWtyc3OZNm0aF198MZdffjnr16/nmWeeYeDAgX6ris/OzgbgtttuY/jw4TidTn7+85/75bVFROykhFRExA8efPBBWrduzdNPP82dd95Jy5Yt+dWvfsXDDz/st92ZrrzySm699Vbmz5/P3LlzMQxDCamIRATtZS8i4oMXX3yRm266ia1bt3LSSSfZHY6ISFjTHFIRER9s374dh8NBy5Yt7Q5FRCTsacheRMSC4uJiXnvtNWbOnMmQIUNISkqyOyQRkbAX8B7SZ599lt69e5OSkkJKSgpDhgzhvffeC/RpRUQCYu3atfzmN7+hS5cuzJkzx+5wREQiQsDnkL799ts4nU66du2KYRi89NJLPProo6xcuZJevXoF8tQiIiIiEgZsKWpq2bIljz76KDfeeGOwTy0iIiIiISaoc0jdbjf/+Mc/KC8vZ8iQIUdtV1VVRVVVlfexx+Nh9+7dtGrVSjuViIiIiIQgwzDYt28f7du3JybG2qzQoCSkq1evZsiQIRw4cIDk5GTeeOMNevbsedT2eXl5TJs2LRihiYiIiIgf+bIcXlCG7KurqyksLKS0tJTXXnuNv/71r3z88cdHTUoP7yEtLS0lKyuLs7mEWPyzwLSIiIiI+E8tNXzCu+zdu5fU1FRLx9oyhzQnJ4fOnTvz3HPPNal9WVkZqampDOMKYh1KSEVERERCTa1Rw2LeorS0lJSUFEvH2rIwvsfjadADKiIiIiLRK+BzSHNzcxkxYgRZWVns27ePefPmsXjxYt5///1An1pEREREwkDAE9IdO3YwZswYtm/fTmpqKr179+b999/nwgsvDPSpRURERCQMBDwhffHFFwN9ChEREREJY7bMIRURERERqaeEVERERERspYRURERERGylhFREREREbKWEVERERERspYRURERERGylhFREREREbKWEVERERERspYRURERERGylhFREREREbKWEVERERERspYRURERERGylhFREREREbKWEVERERERspYRURERERGylhFREREREbKWEVERERERspYRURERERGylhFREREREbKWEVERERERspYRURERERGwVa3cAIuHC4XRaam+43SF5Dqti4uMttfdUVwcoEhERiVTqIRURERERWykhFRERERFbKSEVEREREVspIRURERERWykhFRERERFbqcpepImCUdFu9RxWq/J9OYeq5kVEJNDUQyoiIiIitlJCKiIiIiK2UkIqIiIiIrZSQioiIiIitlJRk0iAWN1yE6wXHAWj0ErE36J1G14ROTr1kIqIiIiIrZSQioiIiIitlJCKiIiIiK0CnpDm5eUxcOBAmjdvTps2bRg1ahTr168P9GlFREREJEwEPCH9+OOPmThxIp9//jkffvghNTU1XHTRRZSXlwf61CIiIiISBgJeZb9w4cIGj+fMmUObNm3Iz8/nnHPOCfTpRY7KahW81S00fdly05mcbKm94cM5VMkvxxKM6vRQ3IY3FPmyUoe2+pVwFfRln0pLSwFo2bLlUdtUVVVRVVXlfVxWVhbwuERERETEHkEtavJ4PNxxxx2cddZZnHbaaUdtl5eXR2pqqveWmZkZxChFREREJJiCmpBOnDiRr7/+mvnz5x+zXW5uLqWlpd7b1q1bgxShiIiIiARb0IbsJ02axL///W+WLFnCSSeddMy2LpcLl8sVpMhERERExE4BT0gNw+DWW2/ljTfeYPHixZx88smBPqWIiIiIhJGAJ6QTJ05k3rx5vPXWWzRv3pyioiIAUlNTSUxMDPTpJUpYrQz2RUwQfl49lZWW2jt8qML1pTLfCl++F5FQER0p7zsUY/KF1e+H1fbBqGZXxbxEk4DPIX322WcpLS1l2LBhtGvXzntbsGBBoE8tIiIiImEgKEP2IiIiIiJHo73sRURERMRWSkhFRERExFZB36lJJBCCUYjhbNPaUntjb6nlc8SkpVk7ID7O8jkI8NahvnwvgrFdZSgKRjGe1a9VoLfUBR8KjnzZQjPABYKh+L3zRbT+7knoUQ+piIiIiNhKCamIiIiI2EoJqYiIiIjYSgmpiIiIiNhKCamIiIiI2EpV9iJNZbG61JHkw1ajLVItNTe2FVk/R4BZrdKGyKjcDcWKaF9YrZoPSkwWK+Z9YvX325f3bbWiPQjvOxJ+9yQyqIdURERERGylhFREREREbKWEVERERERspYRURERERGyloiaJCL4UGFjdCtTyNp3Wd1SEcmtFDIYP2zYGmi9FElaP8eX7HQnFG6H4HoKxVWxMog8FggH+WvmyZapP70MkSqiHVERERERspYRURERERGylhFREREREbKWEVERERERspYRURERERGylKnuJXla38Uu2ViHraZZmqT1ATGWNtfYJLsvnoGyfpeaeffsttfelAt7q1pChWG3uC6tV175soWn1+xGKW6D68jPltvi1srrlbahWzIfi91ukKdRDKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSlb1EBJ8qu3fstNb+pHRL7fedbL0KN3W9tQp4I876OWIsVtlb/trGWr+sBKPa3Cqr79uXvc0Ni8f48nPuyzGBZrWy2+HDOWLbtLbU3qgI/M+U5fftw/dOVfMSrtRDKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UlGThKRgbH/n7Jhpqb3H4usfaGn9/724zGRL7Q0f/qVMaGZti8TYvdaKPRxl5ZbaA3h2lFg7h8VtHiE0tye1eg6r21uCb8VWVvgSk1W+fC8Mi1veYrXgyIf3bblIKQgFaSqCklChHlIRERERsZUSUhERERGxlRJSEREREbGVElIRERERsZUSUhERERGxVVCq7JcsWcKjjz5Kfn4+27dv54033mDUqFHBOLVECWdamuVj3OnNLbWvbR5nqb2zylJzAA6kWfsfMfnHGsvncJZZC8xRstfyOUKR1Wpwq9XpVrc/heBU/luugg/GVqNW34cv7/vkLEvtjViL/TPbdlhrD1Bba6m5T6sLqGpewlRQekjLy8vp06cPM2bMCMbpRERERCSMBKWHdMSIEYwYMSIYpxIRERGRMBOSC+NXVVVRVXVwWLGsrMzGaEREREQkkEKyqCkvL4/U1FTvLTPT2o46IiIiIhI+QrKHNDc3l8mTJ3sfl5WVKSmNMlYn5vtU7FFprSCovJu1bT13ne1DVVOFtV9J9yrr2xfGp1s7R4vyA9ZOsKfUWnvAsFhA5Mu2jVZZ3r7Why06LW8l6QOrvxsOi6/v0/fC5bLU3KjypULQYvFetcUCwfSW1tqD5d8NY7/F7U9FwlhIJqQulwuXxQuWiIiIiISnkByyFxEREZHoEZQe0v3797Nhwwbv402bNrFq1SpatmxJVpa1teJEREREJLIEJSFdsWIF5513nvdx/fzQsWPHMmfOnGCEICIiIiIhKigJ6bBhwzAMIxinEhEREZEwE5JFTSJWOZKbWT7GnWhtK9DK1hbrj93Wp2gnZpRbal+TlGL5HM4qa++jpl2qpfZxPlTZO1tZrFh2eyyfw2qltidCKpytVvJbrZr3ZYULR6y1Pz0xHdpZPgc11qrmjZbWfpccpdZ+V4PF8uoQ2mpUQoSKmkRERETEVkpIRURERMRWSkhFRERExFZKSEVERETEVkpIRURERMRWqrKXiGDst17xWtWqvaX2Fe2tLV2WmFZpqT2Ax7BWAe9OtHwK3BXWzhG729rX1lNRYak9QEyatUp+z17rlfxWq8ctV6f7sC+9p7raUntncrLlcwS6ijomKcnyMZZjKrf+M4XV71+txZUbEnzY3tpprQ/IsXev5VOoal7ClXpIRURERMRWSkhFRERExFZKSEVERETEVkpIRURERMRWKmqSkORMS7PU3tHSWnufdLBWpNQzo9jyKTo322mp/aulAy2fw1FrbcvUqvbWtlR0xXS01B7As2GLpfZWC5R8EYwtGGOsbtNpsQjKF1YLrWIsbgMK1r+2njYWt5YFYvZZK4TyJFj7XsQc8OF7UW5tO1Nffs6NSuvFlCKhQD2kIiIiImIrJaQiIiIiYislpCIiIiJiKyWkIiIiImIrJaQiIiIiYitV2UtI8uzbZ6m9M916Fe6+TGvV5q3T9lhq/2jWG5baA8zda7Fqvsr6/5RJO6xtgVqTbK0i2ocNFYnp0M5Se/cma1X5EPjKfF+2DrVamR+Mc8Q0b26pvS8xkZhgqXltmrX2AE6Xxa1DqyyukuC2uNUoQLy1a47DZf23yWFxlQRtNSqhQj2kIiIiImIrJaQiIiIiYislpCIiIiJiKyWkIiIiImIrFTVJSLI60d5wWf9RblZs7Ryl1dYKEh4pushSe4BRrfIttf9P126Wz7F9T3tL7dt9bq0IylFda6k9AHtKLTW3urUsWC+UsyoYxSFWt/UE69uTWi5SsligBIDTWl+Is9LalpsAjlqLRUdOh6XmnmbW33dMubX2RoW17U9BRUoSvtRDKiIiIiK2UkIqIiIiIrZSQioiIiIitlJCKiIiIiK2UkIqIiIiIrZSlb2EJKuVvo79lZbPEbff2haJMTHWqs3vyvjQUnuAVVUdLLWPi7FeUeuJs1g1X2utvRFr/f9cR3KStXOUWq+Yt7olplVGVZXlYzyV1n5ufdqm0yKrqxHEWNwOEwCr25n6UGXvbm6tCr4yw9o2nYnF1r/fMQesrZIQ06a15XN4fthm+RiRUKAeUhERERGxlRJSEREREbGVElIRERERsZUSUhERERGxlRJSEREREbGVquwlJDks7r/ti+L+1qpqz8pYa6n9Lo/1va6bO61VXW/e0cryOVx7rO3ZbZW7eaLlY6zWjjsqD1g+h3vXbmvnsFrR7kMFfDCq5mPSUi21Nyosrljhy97pidZ+N4w4618nT5y1/paaJGvt4xOtx+S0uAKFYfFnFiDG4rXTU22t8l8kUNRDKiIiIiK2ClpCOmPGDDp16kRCQgKDBw9m+fLlwTq1iIiIiISwoCSkCxYsYPLkyUydOpWCggL69OnD8OHD2bFjRzBOLyIiIiIhLCgJ6RNPPMH48eMZN24cPXv2ZObMmSQlJTFr1qxG21dVVVFWVtbgJiIiIiKRKeBFTdXV1eTn55Obm+t9LiYmhpycHJYuXdroMXl5eUybNi3QoUkEqe1gvbin+Q8eS+3/8103S+1nntT4z/exzClLs9Q+Pr7W8jmqWljcOtTal4maVOsFac5yi9swui0GhfWtQ42KCsvnCDRftj/17C0NQCQHGT4UNcVYLObyZTvamBprPyPNiqwV99QmWS9qcu222LnislZ4CdZ/bq0W1vny/RZpioD3kJaUlOB2u8nIyGjwfEZGBkVFRY0ek5ubS2lpqfe2devWQIcpIiIiIjYJyWWfXC4XLh/+MxQRERGR8BPwHtL09HScTifFxcUNni8uLqZt27aBPr2IiIiIhLiAJ6Tx8fFkZ2ezaNEi73Mej4dFixYxZMiQQJ9eREREREJcUIbsJ0+ezNixYxkwYACDBg1i+vTplJeXM27cuGCcXkRERERCWFAS0tGjR7Nz506mTJlCUVERffv2ZeHChUcUOonUMyxuZ+cstV4RXZmeZKl9q7T9ltqvqbEeU7MYa9XmZ5y02fI5Fm87zVJ7d4K1gZTk1cXHb2QDz759ltpb3b7WU2lxy02CU+Ec6Kpop8WtSQGorrHU3IjxYTDPY201iZpka38Ok9ftstQegGbWrjnss3bNAevfb1XZS6gIWlHTpEmTmDRpUrBOJyIiIiJhQnvZi4iIiIitlJCKiIiIiK2UkIqIiIiIrUJyYXwRqxxWt+QDEktaWmpftDPFUvtvqtpZag/w6o6Bltp/tyvd8jlS11v7P7Qsy9rrJ21tZu0AwFFlrciFHRa3GsV68QYhWLxhtTALAl/k4pNka8U9sTusb39a28ZasVXCTms/U7UtfPg591jbzjSmJPB/oj0WC0hFAkU9pCIiIiJiKyWkIiIiImIrJaQiIiIiYislpCIiIiJiKyWkIiIiImIrVdlLRDBaWquABzBiHJba3zxgiaX26w9Yr7K/us0Xltr/dsNPLZ8j1eK/oS3XWavCjdlbbu0EAPutHWNtU0iT1a1AibV2eQzKf/e+VMBbrKK2WpXv3rXbUnsAZ1qatQNaWN+eNPaHEkvta7NaW2rv3OfDVrGlgf85FwlX6iEVEREREVspIRURERERWykhFRERERFbKSEVEREREVspIRURERERW6nKXiKC8d1my8fEndrCUvu/rj7LUvtBHbdYag/w6a5TrB1QYf1XONZicXDCJmtV1EaixWp2wNhebPkYy+eIgD27ra0LUXdMgPemj2nePKCvD8ABa/vMA7jbt7TU3lFrbXUBYnzoz7G64kGV9fdtdZUEkVChHlIRERERsZUSUhERERGxlRJSEREREbGVElIRERERsZWKmiQkWZ2Yb3krQqAqxdr/Y84NSZbar03OsNQeoKra2q9k3B7rBSuJu6x9ba0WKTlK9lpqD9a/35a3AQXL2256glAEZbXgyJFk7WcQwL13r+VjLKmttX6M1e+FD9uTxrSwVmzlKD9gqb2nufXvBdU1lppHQiGeSFOph1REREREbKWEVERERERspYRURERERGylhFREREREbKWEVERERERspSp7iQxO6/9bpX1nbQ/NfR2tVdXuW2tta1KA2hbWKpZTdlk+BYnbrVUTGxa3SHTEx1lqDxBzUntL7Y0dJZbPYbWy22oFvC9baFrdGtKoqLB8jhhfViQIMKuV/758bR3FFn85LH6dYrbvtPb6gMfi98+XbUCt/txqq1EJFeohFRERERFbKSEVEREREVspIRURERERWykhFRERERFbqahJIoJRYa1ACSB2m7Wih7TvEiy139Pd+v97KRutFQTFVhqWz1GTaq14I7Fkn6X2Rqm19sFiudjD4ut79ll/35a3TLX4Hnw5h1UxPsRktUjJkdzM8jmM/eXWzmH1fbg91tqD5cI6kWiiHlIRERERsZUSUhERERGxlRJSEREREbFVwBPSP/zhD5x55pkkJSWRlpYW6NOJiIiISJgJeEJaXV3Nz372MyZMmBDoU4mIiIhIGAp4lf20adMAmDNnTqBPJVHMpy32LFbJOqutvX6CD7tbOiwW7rpdDsvnSPjRYjV4WRCq5mutbZnqy/fbqLb2DbR6jmBs0RmMbR6tvg9PpQ8rXKSlWjvHDuvbdDoC/f3wYbtiq5X8obiqgkighOSyT1VVVVQdssdzWVmZjdGIiIiISCCFZFFTXl4eqamp3ltmZqbdIYmIiIhIgPiUkN5zzz04HI5j3tatW+dzULm5uZSWlnpvW7du9fm1RERERCS0+TRkf9ddd3H99dcfs80pp5ziy0sD4HK5cLlcPh8vIiIiIuHDp4S0devWtG79/+3da2xU1RrG8WcKdFroRQoDFGkLLVE0RmIqt0q0VaL1JGI/wDexKCFACpFADPWChBiDBBJNgCDGpPgBgkRSS8gxQpBLThREtCoITVrFShuOhTpTGDhT6Ozz4WAPKFTW1j1rD/P/JZMw7dqz3mFN6cueefYK/d21AAAAIAV5HmpqbW1VZ2enWltb1dPTo8bGRknS2LFjlZWV5fX0SBEBN2fUB2UaDU/vMkuCBwea7UsvSekXzGL2VzLMU/aBS4Zp8+7LZo/vIhl8JRw2Gp+IRLupRKSb3fzdmh5j+jzc1BQPR8zmcLPe/c1+vTmX/mM2/prg7S0jAQ/clOcN6Wuvvab333+/9/4DDzwgSdq3b5/Ky8u9nh4AAAA+53nKfvPmzXIc5w83mlEAAABIPr3sEwAAAFIHDSkAAACs8uVOTYCp+Hnz7S37GYYx+g8xC+Hd0WQeenDSzP6P2P/XqPEcipsFpwLpZuGs+PkLRuMlKS3TLGBmug2oG6Y1udlCMxH8uJVkwLCmuIv1Nt761XC93QgMHGh2gIvglOPT1yHwZzhDCgAAAKtoSAEAAGAVDSkAAACsoiEFAACAVTSkAAAAsIqUPW4LbpLEpsf0PxM2m8AwMS9Jil40Gz/IMLUrSRHzKxJ4LgFJcOPXSAKS/IngZmtPE25+9kxT84nYKtb0yg1p2dnGc7i5GgiQKjhDCgAAAKtoSAEAAGAVDSkAAACsoiEFAACAVYSacFtwE3owDRikXbliND4wbKjReEmKD8szGp92+t/Gc5hyTLcv9GnIxTR848ctN93w4/MwDVq5eQ6eh7kuGgYQ5c+1APyCM6QAAACwioYUAAAAVtGQAgAAwCoaUgAAAFhFQwoAAACrSNnDl7xOyEpSWmam0XjTrQXVGTYbLynt0n+MxrtJ7faEw0bj+2VlGY1PRCLaNJUPfzF9jbj598D0NWI6RyL+jXIzB0l+JCvOkAIAAMAqGlIAAABYRUMKAAAAq2hIAQAAYBUNKQAAAKwiZQ9fSsS+42mGCda07Gyj8W5qcgxT9m5SuKb7wMcvXTKew48SsX+6KeO14OoCvuFmLfz4GgT8gjOkAAAAsIqGFAAAAFbRkAIAAMAqGlIAAABYRagJtwVX4Z5swy0xL3of7jHdntRxMYdpGCMRwZtEbMPox4CIH2vyI1ehRQJjQFLhDCkAAACsoiEFAACAVTSkAAAAsIqGFAAAAFbRkAIAAMAqT1P2p06d0uuvv65PP/1UZ86c0ciRI/XMM8/olVdeUbphAhKpJRFb7PWc6zQaHzB8zbpJjps+Dzdz9MsyvLoA6WPPkLL3jh9T86w3cHOeNqQnT55UPB7Xpk2bNHbsWB07dkxz585VNBrV2rVrvZwaAAAAScLThrSyslKVlZW994uLi9XU1KSNGzfSkAIAAECShQvjRyIR5eXl9TkmFospFov13u/q6vK6LAAAAFiS0FBTc3Oz1q1bp3nz5vU5btWqVcrNze29FRQUJKhCAAAAJFrAcRzj3Qdra2u1evXqPsecOHFC48aN673f1tamRx55ROXl5Xrvvff6PPZGZ0gLCgpUrqfVPzDAtFzgb5GWmWk0PhFhINOglWRel2lwyo9hEqSeRAQjAVzvinNZ+9WgSCSinJwco2NdvWW/dOlSzZ49u88xxcXFvX9ub29XRUWFysrK9O677/7p4weDQQWDQTelAQAAIMm4akhDoZBCodAtjW1ra1NFRYVKS0tVV1entDQufQoAAID/8zTU1NbWpvLychUVFWnt2rXq6Ojo/d6IESO8nBoAAABJwtOGdM+ePWpublZzc7NGjRp13fdcfHQVAAAAtyFP3z+fPXu2HMe54Q0AAACQLFyHFPCCmy00TZmm092kdtMMU/NukvymdSUifZyQ9SNFnVJYbyC5kDACAACAVTSkAAAAsIqGFAAAAFbRkAIAAMAqGlIAAABYRcoetwU3idpEJLtNme4Db5rKl6RAApL8pkhEA0Bq4wwpAAAArKIhBQAAgFU0pAAAALCKhhQAAABWEWpCyvI6SOMmcGQaanL1HHy4dSgAILVxhhQAAABW0ZACAADAKhpSAAAAWEVDCgAAAKtoSAEAAGAVKXvAI6aJeSkx25mSmgcA+A1nSAEAAGAVDSkAAACsoiEFAACAVTSkAAAAsIpQE1KW1wEiN+GhRASOTJ93wHAL1PilS0bjAQDgDCkAAACsoiEFAACAVTSkAAAAsIqGFAAAAFbRkAIAAMAqUvZIWV4n2m+XbUBJzQMAvMYZUgAAAFhFQwoAAACraEgBAABgFQ0pAAAArKIhBQAAgFWk7AGPJCIB74Zf6wIApC7OkAIAAMAqzxvS6dOnq7CwUBkZGcrPz9esWbPU3t7u9bQAAABIEp43pBUVFdq+fbuampq0Y8cOtbS0aMaMGV5PCwAAgCQRcBzHSeSEO3fuVFVVlWKxmAYMGHDDMbFYTLFYrPd+JBJRYWGhpuof6q8bHwMAAAB7ruiy/qV/KhwOKzc31+jYhIaaOjs7tWXLFpWVld20GZWkVatWaeXKlX/4+r/0Ty/LAwAAwF907tw544Y0IWdIly1bpvXr1+vixYuaPHmydu3apSFDhtx0/O/PkIbDYRUVFam1tdX4CcJbXV1dKigo0M8//6ycnBzb5eAarI2/sT7+xdr4F2vjb7+9o/3rr7/qjjvuMDrWVUNaW1ur1atX9znmxIkTGjdunCTp7Nmz6uzs1E8//aSVK1cqNzdXu3btUiAQuKX5urq6lJubq0gkwgvQZ1gb/2Jt/I318S/Wxr9YG3/7K+vj6i37pUuXavbs2X2OKS4u7v3z0KFDNXToUN1111265557VFBQoEOHDmnKlClupgcAAMBtxFVDGgqFFAqFXE0Yj8cl6bq35AEAAJC6PA01HT58WEeOHNHUqVM1ePBgtbS0aPny5SopKTE6OxoMBrVixQoFg0EPq4UbrI1/sTb+xvr4F2vjX6yNv/2V9fE01PTdd9/phRde0DfffKNoNKr8/HxVVlbq1Vdf1Z133unVtAAAAEgiCb8OKQAAAHAt9rIHAACAVTSkAAAAsIqGFAAAAFbRkAIAAMCqpGtIp0+frsLCQmVkZCg/P1+zZs1Se3u77bJS3qlTpzRnzhyNGTNGmZmZKikp0YoVK9Td3W27NFz1xhtvqKysTAMHDjTe0g1/rw0bNmj06NHKyMjQpEmT9MUXX9guCZIOHjyop556SiNHjlQgENBHH31kuyRctWrVKk2YMEHZ2dkaNmyYqqqq1NTUZLssSNq4caPuv/9+5eTkKCcnR1OmTNHHH39s/DhJ15BWVFRo+/btampq0o4dO9TS0qIZM2bYLivlnTx5UvF4XJs2bdLx48f11ltv6Z133tHLL79suzRc1d3drZkzZ2rBggW2S0lpH3zwgZYsWaIVK1boq6++0vjx4/XEE0/ol19+sV1ayotGoxo/frw2bNhguxT8zoEDB1RTU6NDhw5pz549unz5sh5//HFFo1HbpaW8UaNG6c0339TRo0f15Zdf6tFHH9XTTz+t48ePGz1O0l/2aefOnaqqqlIsFtOAAQNsl4NrrFmzRhs3btQPP/xguxRcY/PmzVq8eLHC4bDtUlLSpEmTNGHCBK1fv17S/3avKygo0KJFi1RbW2u5OvwmEAiovr5eVVVVtkvBDXR0dGjYsGE6cOCAHn74Ydvl4Hfy8vK0Zs0azZkz55aPSbozpNfq7OzUli1bVFZWRjPqQ5FIRHl5ebbLAHyju7tbR48e1bRp03q/lpaWpmnTpunzzz+3WBmQXCKRiCTxO8Znenp6tG3bNkWjUaMdOaUkbUiXLVumQYMGaciQIWptbVVDQ4PtkvA7zc3NWrdunebNm2e7FMA3zp49q56eHg0fPvy6rw8fPlxnzpyxVBWQXOLxuBYvXqyHHnpI9913n+1yoP/tzJmVlaVgMKj58+ervr5e9957r9Fj+KIhra2tVSAQ6PN28uTJ3vEvvviivv76a+3evVv9+vXTs88+qyT/5IFvma6NJLW1tamyslIzZ87U3LlzLVWeGtysDwAks5qaGh07dkzbtm2zXQquuvvuu9XY2KjDhw9rwYIFqq6u1vfff2/0GL74DGlHR4fOnTvX55ji4mKlp6f/4eunT59WQUGBPvvsM+PTw/hzpmvT3t6u8vJyTZ48WZs3b1Zami/+z3PbcvOzw2dI7enu7tbAgQP14YcfXvfZxOrqaoXDYd7t8RE+Q+pPCxcuVENDgw4ePKgxY8bYLgc3MW3aNJWUlGjTpk23fEx/D+u5ZaFQSKFQyNWx8XhckhSLxf7OknCVydq0tbWpoqJCpaWlqquroxlNgL/ys4PES09PV2lpqfbu3dvb6MTjce3du1cLFy60WxzgY47jaNGiRaqvr9f+/ftpRn0uHo8b92W+aEhv1eHDh3XkyBFNnTpVgwcPVktLi5YvX66SkhLOjlrW1tam8vJyFRUVae3atero6Oj93ogRIyxWht+0traqs7NTra2t6unpUWNjoyRp7NixysrKsltcClmyZImqq6v14IMPauLEiXr77bcVjUb13HPP2S4t5V24cEHNzc2993/88Uc1NjYqLy9PhYWFFitDTU2Ntm7dqoaGBmVnZ/d+5jo3N1eZmZmWq0ttL730kp588kkVFhbq/Pnz2rp1q/bv369PPvnE7IGcJPLtt986FRUVTl5enhMMBp3Ro0c78+fPd06fPm27tJRXV1fnSLrhDf5QXV19w/XZt2+f7dJSzrp165zCwkInPT3dmThxonPo0CHbJcFxnH379t3wZ6S6utp2aSnvZr9f6urqbJeW8p5//nmnqKjISU9Pd0KhkPPYY485u3fvNn4cX3yGFAAAAKmLD/kBAADAKhpSAAAAWEVDCgAAAKtoSAEAAGAVDSkAAACsoiEFAACAVTSkAAAAsIqGFAAAAFbRkAIAAMAqGlIAAABYRUMKAAAAq/4L5+XL3cE4cU4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x800 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\" plot data \"\"\"\n",
    "\n",
    "plot_samples(data, (double_well_potential, harmonic_potential), label=\"data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define flow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "\"\"\" Define the flows. \"\"\"\n",
    "\n",
    "import distrax\n",
    "import equinox as eqx\n",
    "import lenses\n",
    "\n",
    "from flox.flow import (\n",
    "    Pipe,\n",
    "    SimpleCoupling,\n",
    "    Transform,\n",
    "    Transformed,\n",
    "    Affine,\n",
    "    DistraxWrapper\n",
    ")\n",
    "from flox.util import Lens\n",
    "\n",
    "NUM_FLOW_LAYERS = 3\n",
    "NUM_HIDDEN = 128\n",
    "NUM_BINS = 8\n",
    "\n",
    "from flox._src.nn.modules import gen_dense, KeyArray\n",
    "\n",
    "\n",
    "class MySplineTransform(eqx.Module):\n",
    "\n",
    "    num_bins: int\n",
    "    context: Callable[[State], VectorN]\n",
    "    net: eqx.nn.Sequential\n",
    "\n",
    "    def __init__(self, context: Callable[[State], VectorN], num_in: int, num_hidden: int, num_bins: int, dims: int, activation=jax.nn.silu, * , key: KeyArray):\n",
    "        chain = key_chain(key)\n",
    "\n",
    "        self.context = context\n",
    "\n",
    "        num_out = (3 * num_bins + 1) * dims\n",
    "        self.num_bins = num_bins\n",
    "        self.net = eqx.nn.Sequential(tuple(gen_dense([num_in, num_hidden, num_out], activation, next(chain))))\n",
    "\n",
    "    def __call__(self, s: State) -> DistraxWrapper:\n",
    "        params = self.net(self.context(s))\n",
    "        print(params.shape)\n",
    "        params = params.reshape(-1,  3 * self.num_bins + 1)\n",
    "        return DistraxWrapper(\n",
    "            distrax.RationalQuadraticSpline(\n",
    "                params, range_min=-5., range_max=5.\n",
    "            )\n",
    "        )\n",
    "\n",
    "class MyAffineTransform(eqx.Module):\n",
    "\n",
    "    context: Callable[[State], VectorN]\n",
    "    net: eqx.nn.Sequential\n",
    "\n",
    "    def __init__(self, context: Callable[[State], VectorN], num_in: int, num_hidden: int, dims: int, activation=jax.nn.silu, *, key: KeyArray):\n",
    "        chain = key_chain(key)\n",
    "\n",
    "        self.context = context\n",
    "\n",
    "        self.net = eqx.nn.Sequential(tuple(gen_dense([num_in, num_hidden, 2 * dims], activation, next(chain))))\n",
    "\n",
    "    def __call__(self, s: State) -> Affine:\n",
    "        shift, scale = jnp.split(self.net(self.context(s)), 2)\n",
    "        return Affine(shift, scale * 1e-1)\n",
    "\n",
    "\n",
    "def get_aux(s: State) -> VectorN:\n",
    "    return s.aux\n",
    "\n",
    "def update_aux(s: State, aux: VectorN) -> State:\n",
    "    return lenses.bind(s).aux.set(aux)\n",
    "    \n",
    "def get_pos(s: State) -> VectorN:\n",
    "    return s.pos\n",
    "\n",
    "def update_pos(s: State, pos: VectorN) -> State:\n",
    "    return lenses.bind(s).pos.set(pos)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "State -> Transformed[State] = (State, Volume)\n",
    "\n",
    "Seq: T -> \n",
    "\n",
    "Wrap: T -> V[T]\n",
    "Seq: [ ... V[T] ... ] -> V[T]\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "flow = Pipe([\n",
    "    SimpleCoupling(\n",
    "        MyAffineTransform(get_pos, num_in=1, num_hidden=NUM_HIDDEN, dims=1, key=next(chain)),\n",
    "        Lens(get_aux, update_aux)\n",
    "    ),\n",
    "    # SimpleCoupling(\n",
    "    #     MyAffineTransform(get_aux, num_hidden=NUM_HIDDEN, dims=1),\n",
    "    #     Lens(get_pos, update_pos)\n",
    "    # ),\n",
    "    SimpleCoupling(\n",
    "        MySplineTransform(get_aux, num_in=1, num_hidden=NUM_HIDDEN, num_bins=NUM_BINS, dims=1, key=next(chain)),\n",
    "        Lens(get_pos, update_pos)\n",
    "    )\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "from jaxtyping import jaxtyped, Array, Float  # type: ignore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": [
    "initializer = sample_base(next(chain))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "callable passed as argument \"fn\" has too many arguments in its declaration; expected 1 but 2 argument(s) declared",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/home/jonas/dev/PhD/so3/example.ipynb Cell 19\u001b[0m in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> <a href='vscode-notebook-cell:/home/jonas/dev/PhD/so3/example.ipynb#X25sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m flow\u001b[39m.\u001b[39;49mforward(initializer\u001b[39m.\u001b[39;49mobj)\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/jonas/dev/PhD/so3/example.ipynb#X25sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m \u001b[39m# params = jax.jit(flow.pure.init)(next(chain), initializer.obj)\u001b[39;00m\n",
      "File \u001b[0;32m~/miniconda3/envs/flox-debug/lib/python3.10/site-packages/jaxtyping/decorator.py:42\u001b[0m, in \u001b[0;36m_Jaxtyped.__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     40\u001b[0m memo_stack\u001b[39m.\u001b[39mappend(({}, {}, {}))\n\u001b[1;32m     41\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m---> 42\u001b[0m     \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mfn(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m     43\u001b[0m \u001b[39mfinally\u001b[39;00m:\n\u001b[1;32m     44\u001b[0m     memo_stack\u001b[39m.\u001b[39mpop()\n",
      "File \u001b[0;32m~/miniconda3/envs/flox-debug/lib/python3.10/site-packages/typeguard/__init__.py:1033\u001b[0m, in \u001b[0;36mtypechecked.<locals>.wrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m   1031\u001b[0m memo \u001b[39m=\u001b[39m _CallMemo(python_func, _localns, args\u001b[39m=\u001b[39margs, kwargs\u001b[39m=\u001b[39mkwargs)\n\u001b[1;32m   1032\u001b[0m check_argument_types(memo)\n\u001b[0;32m-> 1033\u001b[0m retval \u001b[39m=\u001b[39m func(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m   1034\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m   1035\u001b[0m     check_return_type(retval, memo)\n",
      "File \u001b[0;32m~/dev/PhD/so3/flox/_src/flow/api.py:147\u001b[0m, in \u001b[0;36mPipe.forward\u001b[0;34m(self, input)\u001b[0m\n\u001b[1;32m    145\u001b[0m accum \u001b[39m=\u001b[39m pure(\u001b[39minput\u001b[39m)\n\u001b[1;32m    146\u001b[0m \u001b[39mfor\u001b[39;00m t \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtransforms:\n\u001b[0;32m--> 147\u001b[0m     accum \u001b[39m=\u001b[39m bind(accum, t)\n\u001b[1;32m    148\u001b[0m \u001b[39mreturn\u001b[39;00m cast(Transformed[Output], accum)\n",
      "File \u001b[0;32m~/dev/PhD/so3/flox/_src/flow/api.py:67\u001b[0m, in \u001b[0;36mbind\u001b[0;34m(a, t)\u001b[0m\n\u001b[1;32m     65\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mbind\u001b[39m(a: Transformed[T], t: Transform[T, S]) \u001b[39m-\u001b[39m\u001b[39m>\u001b[39m Transformed[S]:\n\u001b[1;32m     66\u001b[0m     \u001b[39m\"\"\"monadic bind operator\"\"\"\u001b[39;00m\n\u001b[0;32m---> 67\u001b[0m     tb \u001b[39m=\u001b[39m t\u001b[39m.\u001b[39;49mforward(a\u001b[39m.\u001b[39;49mobj)\n\u001b[1;32m     68\u001b[0m     \u001b[39mreturn\u001b[39;00m Transformed(tb\u001b[39m.\u001b[39mobj, tb\u001b[39m.\u001b[39mldj \u001b[39m+\u001b[39m a\u001b[39m.\u001b[39mldj)\n",
      "File \u001b[0;32m~/miniconda3/envs/flox-debug/lib/python3.10/site-packages/jaxtyping/decorator.py:42\u001b[0m, in \u001b[0;36m_Jaxtyped.__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     40\u001b[0m memo_stack\u001b[39m.\u001b[39mappend(({}, {}, {}))\n\u001b[1;32m     41\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m---> 42\u001b[0m     \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mfn(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m     43\u001b[0m \u001b[39mfinally\u001b[39;00m:\n\u001b[1;32m     44\u001b[0m     memo_stack\u001b[39m.\u001b[39mpop()\n",
      "File \u001b[0;32m~/miniconda3/envs/flox-debug/lib/python3.10/site-packages/typeguard/__init__.py:1033\u001b[0m, in \u001b[0;36mtypechecked.<locals>.wrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m   1031\u001b[0m memo \u001b[39m=\u001b[39m _CallMemo(python_func, _localns, args\u001b[39m=\u001b[39margs, kwargs\u001b[39m=\u001b[39mkwargs)\n\u001b[1;32m   1032\u001b[0m check_argument_types(memo)\n\u001b[0;32m-> 1033\u001b[0m retval \u001b[39m=\u001b[39m func(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m   1034\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m   1035\u001b[0m     check_return_type(retval, memo)\n",
      "File \u001b[0;32m~/dev/PhD/so3/flox/_src/flow/api.py:131\u001b[0m, in \u001b[0;36mSimpleCoupling.forward\u001b[0;34m(self, input)\u001b[0m\n\u001b[1;32m    127\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mforward\u001b[39m(\u001b[39mself\u001b[39m, \u001b[39minput\u001b[39m: A) \u001b[39m-\u001b[39m\u001b[39m>\u001b[39m Transformed[A]:\n\u001b[1;32m    128\u001b[0m     transform \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtransform_cstr(\u001b[39minput\u001b[39m)\n\u001b[1;32m    129\u001b[0m     \u001b[39mreturn\u001b[39;00m LensedTransform(\n\u001b[1;32m    130\u001b[0m         transform, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mget_target, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mget_target\n\u001b[0;32m--> 131\u001b[0m     )\u001b[39m.\u001b[39;49mforward(\u001b[39minput\u001b[39;49m)\n",
      "File \u001b[0;32m~/dev/PhD/so3/flox/_src/flow/api.py:87\u001b[0m, in \u001b[0;36mLensedTransform.forward\u001b[0;34m(self, input)\u001b[0m\n\u001b[1;32m     86\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mforward\u001b[39m(\u001b[39mself\u001b[39m, \u001b[39minput\u001b[39m: A) \u001b[39m-\u001b[39m\u001b[39m>\u001b[39m Transformed[B]:\n\u001b[0;32m---> 87\u001b[0m     \u001b[39mreturn\u001b[39;00m unpack_volume(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mforward_lens)(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mtransform\u001b[39m.\u001b[39;49mforward)(\u001b[39minput\u001b[39m)\n",
      "File \u001b[0;32m~/miniconda3/envs/flox-debug/lib/python3.10/site-packages/jaxtyping/decorator.py:42\u001b[0m, in \u001b[0;36m_Jaxtyped.__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     40\u001b[0m memo_stack\u001b[39m.\u001b[39mappend(({}, {}, {}))\n\u001b[1;32m     41\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m---> 42\u001b[0m     \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mfn(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m     43\u001b[0m \u001b[39mfinally\u001b[39;00m:\n\u001b[1;32m     44\u001b[0m     memo_stack\u001b[39m.\u001b[39mpop()\n",
      "File \u001b[0;32m~/miniconda3/envs/flox-debug/lib/python3.10/site-packages/typeguard/__init__.py:1032\u001b[0m, in \u001b[0;36mtypechecked.<locals>.wrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m   1030\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mwrapper\u001b[39m(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs):\n\u001b[1;32m   1031\u001b[0m     memo \u001b[39m=\u001b[39m _CallMemo(python_func, _localns, args\u001b[39m=\u001b[39margs, kwargs\u001b[39m=\u001b[39mkwargs)\n\u001b[0;32m-> 1032\u001b[0m     check_argument_types(memo)\n\u001b[1;32m   1033\u001b[0m     retval \u001b[39m=\u001b[39m func(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n\u001b[1;32m   1034\u001b[0m     \u001b[39mtry\u001b[39;00m:\n",
      "File \u001b[0;32m~/miniconda3/envs/flox-debug/lib/python3.10/site-packages/typeguard/__init__.py:875\u001b[0m, in \u001b[0;36mcheck_argument_types\u001b[0;34m(memo)\u001b[0m\n\u001b[1;32m    873\u001b[0m             check_type(description, value, expected_type, memo)\n\u001b[1;32m    874\u001b[0m         \u001b[39mexcept\u001b[39;00m \u001b[39mTypeError\u001b[39;00m \u001b[39mas\u001b[39;00m exc:  \u001b[39m# suppress unnecessarily long tracebacks\u001b[39;00m\n\u001b[0;32m--> 875\u001b[0m             \u001b[39mraise\u001b[39;00m \u001b[39mTypeError\u001b[39;00m(\u001b[39m*\u001b[39mexc\u001b[39m.\u001b[39margs) \u001b[39mfrom\u001b[39;00m \u001b[39mNone\u001b[39m\n\u001b[1;32m    877\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mTrue\u001b[39;00m\n",
      "\u001b[0;31mTypeError\u001b[0m: callable passed as argument \"fn\" has too many arguments in its declaration; expected 1 but 2 argument(s) declared"
     ]
    }
   ],
   "source": [
    "flow.forward(initializer.obj)\n",
    "# params = jax.jit(flow.pure.init)(next(chain), initializer.obj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<function __main__.bar.<locals>.test(x: ~T)>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from typeguard import typechecked\n",
    "\n",
    "@typechecked\n",
    "class A:\n",
    "\n",
    "    @typechecked\n",
    "    def foo(self, a: int):\n",
    "        print(a)\n",
    "        return a\n",
    "\n",
    "T = TypeVar(\"T\")\n",
    "\n",
    "@typechecked\n",
    "def bar(fn: Callable[[T], T], x: T) -> Callable[[T], T]:\n",
    "\n",
    "    @typechecked\n",
    "    def test(x: T):\n",
    "        return fn(x)\n",
    "    return test\n",
    "\n",
    "\n",
    "# class Lens(Generic[A, B, C, D]):\n",
    "#     project: Callable[[A], C]\n",
    "#     inject: Callable[[A, D], B]\n",
    "\n",
    "#     def __call__(self, fn: Callable[[C], D]) -> Callable[[A], B]:\n",
    "#         def lifted(a: A) -> B:\n",
    "#             return self.inject(a, fn(self.project(a)))\n",
    "\n",
    "#         return lifted\n",
    "\n",
    "a = A()\n",
    "bar(a.foo, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mRunning cells with 'Python 3.10.6 ('flox-debug')' requires ipykernel package.\n",
      "\u001b[1;31mRun the following command to install 'ipykernel' into the Python environment. \n",
      "\u001b[1;31mCommand: 'conda install -n flox-debug ipykernel --update-deps --force-reinstall'"
     ]
    }
   ],
   "source": [
    "import optax\n",
    "\n",
    "NUM_ITERS = 3_000\n",
    "PLOT_INTERVAL = 50\n",
    "LEARNING_RATE = 1e-4\n",
    "\n",
    "optim = optax.adam(learning_rate=LEARNING_RATE)\n",
    "\n",
    "params, static = eqx.partition(flow, eqx.is_array)\n",
    "opt_state = optim.init(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mRunning cells with 'Python 3.10.6 ('flox-debug')' requires ipykernel package.\n",
      "\u001b[1;31mRun the following command to install 'ipykernel' into the Python environment. \n",
      "\u001b[1;31mCommand: 'conda install -n flox-debug ipykernel --update-deps --force-reinstall'"
     ]
    }
   ],
   "source": [
    "initializer = sample_base(next(chain))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mRunning cells with 'Python 3.10.6 ('flox-debug')' requires ipykernel package.\n",
      "\u001b[1;31mRun the following command to install 'ipykernel' into the Python environment. \n",
      "\u001b[1;31mCommand: 'conda install -n flox-debug ipykernel --update-deps --force-reinstall'"
     ]
    }
   ],
   "source": [
    "from flox.nn import train\n",
    "from flox.flow import DatasetSampler\n",
    "\n",
    "step = train.mle_step(joint_base_potential, optim, DatasetSampler(data, data.pos.shape[0]), 128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mRunning cells with 'Python 3.10.6 ('flox-debug')' requires ipykernel package.\n",
      "\u001b[1;31mRun the following command to install 'ipykernel' into the Python environment. \n",
      "\u001b[1;31mCommand: 'conda install -n flox-debug ipykernel --update-deps --force-reinstall'"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "from IPython import display  # type: ignore\n",
    "from IPython.display import clear_output  # type: ignore\n",
    "from tqdm import tqdm\n",
    "\n",
    "pbar = tqdm(range(NUM_ITERS))\n",
    "agg_loss = train.RunningMean(0., 0)\n",
    "losses = []\n",
    "\n",
    "\n",
    "for i in pbar:\n",
    "    \n",
    "    loss, flow, opt_state = step(next(chain), flow, opt_state)  # type: ignore\n",
    "    agg_loss = agg_loss.update(loss)\n",
    "    losses.append(agg_loss.value)\n",
    "    \n",
    "    pbar.set_description(f\"loss: {agg_loss.value:.4}\")\n",
    "    if not i % PLOT_INTERVAL:\n",
    "        \n",
    "        display.clear_output(wait=True)\n",
    "        plt.plot(losses)\n",
    "        plt.ylabel(\"NLL\")\n",
    "        plt.xlabel(\"#iters\")\n",
    "    \n",
    "        plt.show()\n",
    "        time.sleep(0.001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Check inverses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mRunning cells with 'Python 3.10.6 ('flox-debug')' requires ipykernel package.\n",
      "\u001b[1;31mRun the following command to install 'ipykernel' into the Python environment. \n",
      "\u001b[1;31mCommand: 'conda install -n flox-debug ipykernel --update-deps --force-reinstall'"
     ]
    }
   ],
   "source": [
    "a = flow.forward(\n",
    "    initializer.obj\n",
    ")\n",
    "b = flow.inverse(\n",
    "    a.obj\n",
    ")\n",
    "b, initializer, a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sample after training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "python"
    }
   },
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mRunning cells with 'Python 3.10.6 ('flox-debug')' requires ipykernel package.\n",
      "\u001b[1;31mRun the following command to install 'ipykernel' into the Python environment. \n",
      "\u001b[1;31mCommand: 'conda install -n flox-debug ipykernel --update-deps --force-reinstall'"
     ]
    }
   ],
   "source": [
    "from flox.flow import PushforwardSampler\n",
    "\n",
    "out = jax.jit(jax.vmap(PushforwardSampler(sample_base, flow)))(jax.random.split(next(chain), 10_000))\n",
    "\n",
    "weights = jax.nn.softmax(\n",
    "    -jax.vmap(joint_target_potential)(out.obj) + out.ldj\n",
    ")\n",
    "\n",
    "plot_samples(out.obj, (double_well_potential, harmonic_potential), label=\"flow_samples\", weights=weights)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "vscode": {
   "interpreter": {
    "hash": "da2293997c0764481e37fa37a118ef03db119d733d1f187fb156ba91ef764ee6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}