notebooks/iclr24_inverse.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inverse problem [![Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/tumaer/jax-sph/blob/main/notebooks/iclr24_inverse.ipynb)\n",
    "\n",
    "The inverse problem of finding the initial state of a 100-step-long SPH simulation given the final state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"0\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-05-25 04:58:34.723939: W external/xla/xla/service/gpu/nvptx_compiler.cc:698] The NVIDIA driver's CUDA version is 12.2 which is older than the ptxas CUDA version (12.4.99). Because the driver is older than the ptxas version, XLA is disabling parallel compilation, which may slow down compilation. You should update your NVIDIA driver or use the NVIDIA-provided CUDA forward compatibility packages.\n"
     ]
    }
   ],
   "source": [
    "import contextlib\n",
    "import copy\n",
    "\n",
    "import jax\n",
    "import jax.numpy as jnp\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from jax import jit\n",
    "from omegaconf import OmegaConf\n",
    "\n",
    "from jax_sph import partition\n",
    "from jax_sph.case_setup import load_case\n",
    "from jax_sph.defaults import defaults\n",
    "from jax_sph.integrator import si_euler\n",
    "from jax_sph.io_state import read_h5, write_h5\n",
    "from jax_sph.jax_md import space\n",
    "from jax_sph.jax_md.partition import Sparse\n",
    "from jax_sph.simulate import simulate\n",
    "from jax_sph.solver import WCSPH\n",
    "from jax_sph.utils import Tag\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Case A\n",
      "Loading from data/debug/2D_CW_SPH_123_20240525-045900\n",
      "Case B\n",
      "Loading from data/debug/2D_CW_SPH_123_20240525-045903\n"
     ]
    }
   ],
   "source": [
    "# generate reference simulations\n",
    "cfg = OmegaConf.create({\"config\": \"../cases/cw.yaml\"})\n",
    "cfg = OmegaConf.merge(defaults, OmegaConf.load(cfg.config), cfg)\n",
    "cfg.case.viscosity = 0.02\n",
    "cfg.case.r0_noise_factor = 0.02\n",
    "cfg.solver.t_end = 0.5\n",
    "cfg.solver.dt = 0.005\n",
    "cfg.io.write_every = 1\n",
    "\n",
    "# suppress prints\n",
    "print(\"Case A\")\n",
    "with open(os.devnull, \"w\") as f, contextlib.redirect_stdout(f):\n",
    "    simulate(cfg)\n",
    "\n",
    "# newest dir in cfg.io.data_path\n",
    "target_dir_A = max(\n",
    "    [os.path.join(cfg.io.data_path, d) for d in os.listdir(cfg.io.data_path)],\n",
    "    key=os.path.getmtime,\n",
    ")\n",
    "\n",
    "print(f\"Loading from {target_dir_A}\")\n",
    "target_state_A = read_h5(os.path.join(target_dir_A, \"traj_000.h5\"))\n",
    "final_state_A = read_h5(os.path.join(target_dir_A, \"traj_100.h5\"))\n",
    "\n",
    "print(\"Case B\")\n",
    "cfg.case.special.cube_offset = [0.5, 0.3]\n",
    "\n",
    "with open(os.devnull, \"w\") as f, contextlib.redirect_stdout(f):\n",
    "    simulate(cfg)\n",
    "\n",
    "# newest dir in cfg.io.data_path\n",
    "target_dir_B = max(\n",
    "    [os.path.join(cfg.io.data_path, d) for d in os.listdir(cfg.io.data_path)],\n",
    "    key=os.path.getmtime,\n",
    ")\n",
    "print(f\"Loading from {target_dir_B}\")\n",
    "target_state_B = read_h5(os.path.join(target_dir_B, \"traj_000.h5\"))\n",
    "final_state_B = read_h5(os.path.join(target_dir_B, \"traj_100.h5\"))\n",
    "\n",
    "# del cfg.case.special.cube_offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimization Algorithm:\n",
    "# Step 1: Define Loss Function which takes in the random initial position and other \n",
    "#   solver parameters and runs the forward simulation\n",
    "# Step 2: Load the target ground truth position from the paraview file\n",
    "# Step 3: Inside the loss function, compute the mse between the target and predicted r\n",
    "# Step 4: Return the loss and the state computed from forward pass.\n",
    "# Step 5: Define the optimization loop, where the initial state is updated using the \n",
    "#   gradients of the loss: init_state[r] = init_state[r] - learning_rate * grads[-1]\n",
    "\n",
    "def optimization_case_setup(cfg):\n",
    "\n",
    "    Case = load_case(os.path.dirname(cfg.config), cfg.case.source)\n",
    "    case = Case(cfg)\n",
    "    cfg, box_size, state, g_ext_fn, bc_fn, nw_fn, eos_fn, _, displ_fn, shift_fn = case.initialize()\n",
    "\n",
    "    optim_init_state = state\n",
    "    # NOTE make it easier on us\n",
    "    mask = optim_init_state[\"tag\"] == Tag.FLUID\n",
    "    fluid_min = jnp.min(optim_init_state[\"r\"][mask], axis=0)\n",
    "    optim_init_state[\"r\"] = (\n",
    "        optim_init_state[\"r\"]\n",
    "        .at[mask]\n",
    "        .set(optim_init_state[\"r\"][mask] - fluid_min + jnp.array([0.1, 0.5]) + 0.15)\n",
    "    )\n",
    "\n",
    "    displacement_fn, shift_fn = space.periodic(side=box_size)\n",
    "\n",
    "    # Solver setup\n",
    "    solver = WCSPH(\n",
    "        displacement_fn,\n",
    "        eos_fn,\n",
    "        g_ext_fn,\n",
    "        cfg.case.dx,\n",
    "        cfg.case.dim,\n",
    "        cfg.solver.dt,\n",
    "        cfg.case.c_ref,\n",
    "        cfg.solver.eta_limiter,\n",
    "        cfg.solver.name,\n",
    "        cfg.kernel.name,\n",
    "        cfg.solver.is_bc_trick,\n",
    "        cfg.solver.density_evolution,\n",
    "        cfg.solver.artificial_alpha,\n",
    "        cfg.solver.free_slip,\n",
    "        cfg.solver.density_renormalize,\n",
    "        cfg.solver.heat_conduction,\n",
    "    )\n",
    "    forward = solver.forward_wrapper()\n",
    "\n",
    "    # Initialize a neighbors list for looping through the local neighborhood\n",
    "    # cell_size = r_cutoff + dr_threshold\n",
    "    # capacity_multiplier is used for preallocating the (2, NN) neighbors.idx\n",
    "    neighbor_fn = partition.neighbor_list(\n",
    "        displacement_fn,\n",
    "        box_size,\n",
    "        r_cutoff=solver._kernel_fn.cutoff,\n",
    "        backend=cfg.nl.backend,\n",
    "        capacity_multiplier=2.0,\n",
    "        mask_self=False,\n",
    "        format=Sparse,\n",
    "        num_particles_max=state[\"r\"].shape[0],\n",
    "        num_partitions=cfg.nl.num_partitions,\n",
    "        pbc=np.array(cfg.case.pbc),\n",
    "    )\n",
    "    num_particles = (state[\"tag\"] != Tag.PAD_VALUE).sum()\n",
    "    neighbors = neighbor_fn.allocate(state[\"r\"], num_particles=num_particles)\n",
    "\n",
    "    # Instantiate advance function for our use case\n",
    "    advance = si_euler(cfg.solver.tvf, forward, shift_fn, bc_fn, nw_fn)\n",
    "    advance = advance if cfg.no_jit else jit(advance)\n",
    "\n",
    "    # compile kernel and initialize accelerations\n",
    "    _state, _neighbors = advance(0.0, optim_init_state, neighbors)\n",
    "    _state[\"v\"].block_until_ready()\n",
    "\n",
    "    return advance, optim_init_state, neighbors, neighbor_fn, num_particles\n",
    "\n",
    "\n",
    "def forward_simulation(advance, state, neighbors, neighbor_fn, num_particles, cfg):\n",
    "\n",
    "    for step in range(cfg.solver.sequence_length):\n",
    "\n",
    "        state_, neighbors_ = advance(cfg.solver.dt, state, neighbors)\n",
    "\n",
    "        # Check whether the edge list is too small and if so, create longer one\n",
    "        if neighbors_.did_buffer_overflow:\n",
    "            edges_ = neighbors.idx.shape\n",
    "            print(f\"Reallocate neighbors list {edges_} at step {step}\")\n",
    "            neighbors = neighbor_fn.allocate(state[\"r\"], num_particles=num_particles)\n",
    "            print(f\"To list {neighbors.idx.shape}\")\n",
    "\n",
    "            # To run the loop N times even if sometimes did_buffer_overflow > 0\n",
    "            # we directly rerun the advance step here\n",
    "            state, neighbors = advance(cfg.solver.dt, state, neighbors)\n",
    "        else:\n",
    "            state, neighbors = state_, neighbors_\n",
    "\n",
    "    return state\n",
    "\n",
    "\n",
    "# Loss function:\n",
    "def loss_fn_wrapper(advance, cfg):\n",
    "    def loss_fn(state, neighbors, neighbor_fn, num_particles, target_position):\n",
    "        # Forward Simulation\n",
    "        state = forward_simulation(\n",
    "            advance, state, neighbors, neighbor_fn, num_particles, cfg\n",
    "        )\n",
    "        mask = state[\"tag\"] == Tag.FLUID\n",
    "        fluid_position = state[\"r\"][mask]\n",
    "\n",
    "        loss = jnp.sqrt(jnp.mean((fluid_position - target_position) ** 2))\n",
    "        return loss\n",
    "\n",
    "    return loss_fn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimizing case A\n",
      "dt_convective : 0.0011363636363636365\n",
      "dt_viscous    : 0.03124999999999966\n",
      "dt_body_force : 0.055901699437494734\n",
      "dt_max_allowed: 0.0011363636363636365\n",
      "dt_final      : 0.005\n",
      "Total number of particles =  312\n",
      "Step 0, Loss: 0.300\n",
      "Step 1, Loss: 0.279\n",
      "Step 2, Loss: 0.258\n",
      "Step 3, Loss: 0.237\n",
      "Step 4, Loss: 0.216\n",
      "Step 5, Loss: 0.196\n",
      "Step 6, Loss: 0.175\n",
      "Step 7, Loss: 0.154\n",
      "Step 8, Loss: 0.133\n",
      "Step 9, Loss: 0.112\n",
      "Step 10, Loss: 0.091\n",
      "Step 11, Loss: 0.071\n",
      "Step 12, Loss: 0.050\n",
      "Step 13, Loss: 0.029\n",
      "Step 14, Loss: 0.008\n",
      "Best Loss: 0.008\n",
      "Optimizing case B\n",
      "dt_convective : 0.0011363636363636365\n",
      "dt_viscous    : 0.03124999999999966\n",
      "dt_body_force : 0.055901699437494734\n",
      "dt_max_allowed: 0.0011363636363636365\n",
      "dt_final      : 0.005\n",
      "Total number of particles =  312\n",
      "Step 0, Loss: 0.315\n",
      "Step 1, Loss: 0.294\n",
      "Step 2, Loss: 0.273\n",
      "Step 3, Loss: 0.252\n",
      "Step 4, Loss: 0.232\n",
      "Step 5, Loss: 0.211\n",
      "Step 6, Loss: 0.190\n",
      "Step 7, Loss: 0.169\n",
      "Step 8, Loss: 0.148\n",
      "Step 9, Loss: 0.127\n",
      "Step 10, Loss: 0.107\n",
      "Step 11, Loss: 0.086\n",
      "Step 12, Loss: 0.065\n",
      "Step 13, Loss: 0.044\n",
      "Step 14, Loss: 0.023\n",
      "Best Loss: 0.023\n"
     ]
    }
   ],
   "source": [
    "STEPS = 15\n",
    "LR = 1.5\n",
    "\n",
    "# Load the target state\n",
    "for case, states in zip(\n",
    "    (\"A\", \"B\"), [(target_state_A, final_state_A), (target_state_B, final_state_B)]\n",
    "):\n",
    "    print(f\"Optimizing case {case}\")\n",
    "\n",
    "    _, final_state = states\n",
    "\n",
    "\n",
    "    advance, optim_init_state, neighbors, neighbor_fn, num_particles = (\n",
    "        optimization_case_setup(cfg)\n",
    "    )\n",
    "\n",
    "    state = copy.deepcopy(optim_init_state)  # init_state required for plotting\n",
    "    mask = state[\"tag\"] == Tag.FLUID\n",
    "    target_position = final_state[\"r\"][mask]  # Extract only the fluid particles\n",
    "\n",
    "    # Main Optimization Loop:\n",
    "    loss_fn = loss_fn_wrapper(advance, cfg)\n",
    "\n",
    "    write_h5(optim_init_state, f\"data/optim_{case}_init.h5\")\n",
    "    \n",
    "    best_loss = 1e6\n",
    "    best_state = state\n",
    "    mask = state[\"tag\"] == Tag.FLUID\n",
    "\n",
    "    for optimization_step in range(STEPS):\n",
    "\n",
    "        loss, grad = jax.value_and_grad(fun=loss_fn, allow_int=True)(\n",
    "            state, neighbors, neighbor_fn, num_particles, target_position\n",
    "        )\n",
    "\n",
    "        state[\"r\"] = state[\"r\"].at[mask].set(state[\"r\"][mask] - LR * grad[\"r\"][mask])\n",
    "\n",
    "        if loss < best_loss:\n",
    "            best_loss = loss\n",
    "            best_state = copy.deepcopy(state)\n",
    "\n",
    "        print(f\"Step {optimization_step}, Loss: {loss:.3f}\")\n",
    "    \n",
    "    print(f\"Best Loss: {best_loss:.3f}\")\n",
    "\n",
    "    write_h5(best_state, f\"data/optim_{case}_final.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "init_optim_A = read_h5(\"data/optim_A_init.h5\")\n",
    "final_optim_A = read_h5(\"data/optim_A_final.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAHqCAYAAADLbQ06AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABXIElEQVR4nO3de3xMZ/4H8M+ZkAsRQZSwCJGWoi6lokGwSLZ1a0uV3W2iTbSLqrrs0huWtqro/qo3l7oUtW1Xq3SLat1VQkrc6hIkqUXcb4m4ZZ7fHzHDmJmYSU5OvjP5vF+vvLI9c+Y7zzMbz2fmzJnz1ZRSCkRERCSSqaQHQERERM4xqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREgjGoiYiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREgjGoiYiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMVg4yMDGiahvj4+GJ7jA4dOkDTtGKrT0QyMKiJCiksLAxhYWElPQwi8nJljH7A1ORdmDVlPtL2HUZEw3AkjoxDs9YPeWQNImdq1qyJffv2oWLFisX2GJ9//jmuXLlSbPWJSpqUtb6k80JTSimjHiw1eRcSur8EBQVznhkmHxM0aJi9fLrLk5ZSg8jybjojI6NEx0HkjaSs9RLywtBD37OmzLdOFgDMeWYoKMyaMt/japDnmDt3Llq3bo3AwEAEBgaidevWmDdvns0+69atg6ZpGDduHDZt2oQOHTqgQoUKCA4OxlNPPYVDhw5Z97V8/pyZmYnMzExommb9GTdunM0+d39Gbflc+dq1a3j11VdRu3ZtBAQE4OGHH8ZPP/0EALh48SIGDx6MGjVqwN/fH23atMHWrVvt5uXoM+o7x+Lo5+55p6enIyEhAbVr14afnx9CQ0MRHx+PzMxMu8fTNA0dOnTAsWPH8Oyzz6J69eowmUxYt26da/9HELlBylovIS8MPfSdtu+wdbIW5jwz0vYd9rga5BmGDh2K6dOno2bNmnj++ecBAEuWLMGAAQOwY8cO/N///Z/N/klJSXjnnXcQGxuLl156CXv37sW3336LjRs3IikpCfXq1UNwcDDGjh2Lf/3rXwCAYcOGWe/foUMHl8bVt29f7N69Gz169EBubi4WLVqEbt26YfPmzRg4cCCuX7+OPn364PTp0/jyyy8RGxuL9PT0ex5KHzt2rMPtn3zyCU6dOoVy5cpZtyUnJyMmJgY5OTno1q0bIiIikJGRgUWLFmHFihXYsmUL6tWrZ1Pn7NmzaNOmDSpXroxnnnkGV69eRVBQkEtzJnKHlLVeQl4YGtQRDcNx9tQ5m0mbfEyIaBjucTVIvg0bNmD69Olo2LAhtmzZYg25cePGITIyEh988AF69+6Ndu3aWe+zatUqfPrpp3jhhRes22bMmIEXX3wRL7/8MpYvX47g4GCMGzfO+u7U8i7aHWfPnsWuXbtQvnx5AEBMTAz69u2Lzp07o0uXLvjiiy9Qpkz+P89mzZrhH//4Bz777DMMHz68wLqOxvLuu+/i1KlT6NmzJ3r37g0AuHHjBp555hmYzWZs3boVzZs3t+5vOaJgme+d9uzZgwEDBmDWrFnw8fFxe95ErpKy1kvIC0MPfSeOjIMGDSaf/Ie1HOtPHBXvcTVIvvnz8w9NjRs3zuadaKVKlazvPO8+FHz//fcjMTHRZltiYiIiIiLw3//+F6dPn9ZlbG+99ZY1pAGgd+/eKFu2LC5cuIApU6ZYQxoA+vXrBwDYuXOn24/zzTffYMyYMWjRogUWLVoEkyn/b/77779HRkYGRo0aZRPSANC2bVv07NkTP/zwAy5dumRzm6+vLyZPnsyQpmInZa2XkBeGvqNu1vohzF4+3fbsuVHxaPZIE4+rQfLt2LEDgOPD0R07dgQApKam2myPioqyhpmFyWRCVFQU0tLSsHPnTnTu3LnIY2vWrJndY9x33324cuUKateubXNbaGgoAOD48eNuPUZKSgr++te/okaNGli+fLnNC4OkpCQAwIEDBxy+C8/KyoLZbMbBgwfRsmVL6/a6desiJCTErXEQFYaUtV5CXhj+9axmrR/CR19P9YoaJNulS5dgMplQtWpVu9uqVasGTdPs3jFWq1bNYS3L9osXL+oyNkef65YpU8bpdiD/cLWrjh49iu7du0PTNCxfvhw1atSwuf3cuXMAgEWLFhVYJycnx+a/nT0/RMVBylpf0nnBC56Q1woKCoLZbHZ4uPrUqVNQStkF48mTJx3Wsmwvzu9F6+Xy5cvo1q0bTp06hS+++MLu0DZw+4XC8uXLoZRy+hMdHW1zP14Jjch4DGryWpaAcvT1Icu2uw9Bb968GWbzXWd4ms345ZdfoGkamjZtat3u4+ODvLw8XcdcVHl5eXjmmWewa9cuvPfee+jRo4fD/Vq3bg0A2LJli5HDI6JCYFCT14qLiwMAjB8/3uYQ98WLFzF+/HibfSwOHjyIWbNm2WybNWsWDh48iMcff9zmMHrlypVx5swZXL16tbim4LZhw4bhhx9+wMCBAws8Q7xnz56oXbs2pk2bhg0bNtjdfuPGDWzatKk4h0pELjL8M2oio7Rv3x4vvfQSpk+fjsaNG+Opp56CUgpLlizB//73PwwdOhTt27e3uU9MTAyGDh2KH374AY0aNcLevXuxfPlyhISE2H3nulOnTkhJScGf/vQntGvXDr6+vmjfvr1dTaNs3boVH374IQICAlC1alWHJ4n16tULzZo1g5+fH/7zn//gT3/6E6Kjo9GpUyc0adLEeiGXjRs3okqVKti/f7/xEyEiGwxq8moffPABmjdvjk8++QQzZ84EADRq1Aj//Oc/MWDAALv9IyMj8frrr+P111/HBx98AB8fH/Tq1QuTJ0+2u/jHG2+8gfPnz+P777/Hxo0bkZeXh7Fjx5ZYUFuu+52bm4u33nrL4T5hYWHWw/2tWrXCzp078d577+GHH37A5s2b4efnh5o1a6JXr17Wr4URUcky9FrfRFKtW7cOHTt2xNixYwt1ARMiouLCz6iJiIgEY1ATEREJxqAmIiISzPDPqKU08S7pRuBERFQwKWt9SeeFoUEtpYm3hEbgRETknJS1XkJeGHroW0oTbwmNwImIyDkpa72EvDA0qKU08ZbQCJyIiJyTstZLyAtDgzqiYbi1p6d1AIVo4i2hBnmmhQsX4oUXXkDLli3h5+cHTdPselLfady4cdA0zelPRkaG7mPUNM1ha06i0kTKWi8hLwy9MlniyDgkr0+Bycdkc6zf3SbeEmqQZ3r99deRmZmJkJAQhIaGIjMz06X7xcXFISwszG57cHCwvgMkIgBy1noJeVHyZ30XogG3lBrkeX766SdERESgTp06mDRpEsaMGYO5c+ciPj7e4f7jxo3D+PHjsXbtWsPe5WqahujoaIddv4hKEylrfYnnhSIqpd555x0FQM2dO9fpPmPHjlUA1Nq1a4v8eGvWrFGxsbEqNDRU+fr6qvvuu0+1bdtWzZgxQyml1Nq1axUAhz93j3Hp0qWqU6dOKjg4WPn5+alGjRqp9957T928edNmv7lz51rvv3TpUtWqVSsVEBCgQkJC1IABA1RWVpbdOH/99Vf11FNPqVq1ailfX18VEhKiWrZsqSZOnFjk54CI3MemHEQu2LBhA5KTk2EymRAREYHOnTsjMDDQ5fv/97//Rffu3REcHIyePXsiNDQUp0+fxs6dO7FgwQIMHDgQYWFhGDt2LMaPH486derYvMu/s2/2mDFjMGnSJNSsWRNPPvkkKlasiI0bN2LUqFFITk7G119/bff4S5YswapVq9C7d2907twZSUlJmDt3LjZu3IitW7eiUqVKAIDU1FQ8+uij8PHxQc+ePVGnTh1cuHABv/32G2bOnInXXnut0M8hERVSSb9SICop7ryjvvsnODhYzZ8/3+XHevLJJxUAlZqaanfbmTNnbP4bgIqOjnZY58cff1QAVExMjMrOzrZuN5vN6sUXX1QA1H/+8x/rdss7agBq5cqVNrVGjx6tAKghQ4ZYtw0fPlwBUEuXLr3nOInIGLyEKFEBmjZtijlz5uDIkSPIzc1Feno6pk+fDk3TEB8fj2XLlrlVLyAgwG5blSpVXL7/hx9+CACYOXMmypcvb92uaRomTZoETdOwePFiu/t17twZMTExNttee+01BAcH4/PPP4fZbPv1k6KOk4j0w0PfRAV44oknbP47LCwMQ4YMQcOGDdGlSxe8/vrr6NGjxz3rPPPMM/jmm28QGRmJ/v37449//CPatWuHkJAQt8aTlJSE8uXLY86cOQ5vDwgIwP79++22t2vXzm5bYGAgmjVrhnXr1uHIkSOoX78+nn76afzrX//CE088gb59+6JLly5o3749atas6dY4iUg/DGqiQvjjH/+I8PBw7N69G5cuXUJQUFCB+/fp0wdLly7FtGnT8Omnn+Kjjz6Cpmno2LEjpk6davMZdEHOnTuHmzdvYvz48U73ycnJsdtWrVo1h/tatl+8eBEA0Lp1a6xbtw5vv/02vvjiC8ydOxcA0KpVK7z77rvo2LGjS+MkIv3w0DdRIVneDV+5csWl/Xv27In169fj/PnzWLFiBRISErBu3TrExsbiwoULLtUICgpClSpVoJRy+pOenm53v5MnTzqsZ9lesWJF67Z27dphxYoVOH/+PNauXYvhw4dj9+7dePzxx3HkyBGXxklE+mFQExVCTk4O9u7di/Lly7t9+LpChQqIjY3FzJkzER8fj5MnTyI5Odl6u8lkQl5ensP7tm7dGmfPnkVaWppbj7lx40a7bdnZ2UhNTUVQUBDq1atnd3tAQAA6dOiAqVOn4tVXX0Vubi5Wr17t1uMSUdExqImcuHz5Mg4ePGi3PTc3F4mJibh8+TKefvpplClz70+QNmzY4DB8T506BQDw9/e3bqtcuTL+97//OawzdOhQAMBzzz2Hs2fP2t2elZWFffv22W3/6aefsGrVKpttb731Fi5cuIBnn30WJlP+UrBlyxZcvXrV7v6Wd953jpOIjGH4lcmIStLs2bOxadMmAMDu3buxfft2REVFoX79+gCAtm3bIiEhAQCQkZGBevXqoVWrVmjYsCGqV6+OkydP4qeffsL//vc/NGnSBGvXrnXpbOhmzZrh+PHjaNu2LcLCwqBpGjZt2oStW7ciMjISmzZtgo+PDwCgb9+++Oqrr9CzZ080b94cPj4+6NGjBx56KL+l3ptvvokJEyYgODgYsbGxqFOnDs6ePYtDhw5h48aNmDhxIkaPHg0AmDdvHgYMGIBu3bpZv0cdFhaGpKQkrF27FuHh4di2bZv1e9S9evXC2rVr0b59e9StWxf+/v7Yvn07fv75Z9SrVw87duy45+fxRKQzo78PtiNppxrUe7jq0qinGtR7uNqRtNNja5DniYuLc3r1LwAqLi7Ouu/FixfV4MGDVatWrVTVqlVVmTJlVIUKFdQjjzyiJk+erK5cueLy4/773/9WTz/9tAoPD1flypVTFStWVE2bNlXvvvuuunz5ss2+J06cUE8//bQKCQlRJpPJ4Xe9V69erbp3766qVq2qypYtq6pXr67atGmjJkyYoH7//Xfrfs6uTFalShUVHx+vTpw4YVN35cqV6tlnn1UPPPCAqlChggoMDFQPPvigevXVV9Xp06ddf6KJdCBlrS/pvDD0HbWUJt4SGoETGcHyjrqg65kTSSRlrZeQF4Z+Ri2libeERuBEROSclLVeQl4YGtRSmnhLaARORETOSVnrJeSFoUEtpYm3hEbgRETknJS1XkJeGBrUiSPjoEGzTrqwTbwl1CDyBPHx8VBK8fNp8jhS1noJeWH417OkNPEu8UbgRERUIClrfUnnBb9HTUREJBivTEZERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRE5FBGRgY0TUNsbGxJD6XYXLt2Df/85z8REREBf39/1KhRAwMHDrR2NXPHtm3b8NhjjyE4OBjly5dHZGQkvvrqK4f7btiwASNHjkTHjh1RsWJFaJrGr9CRU/fuz0dE5IXMZjN69uyJVatWITIyEk899RTS0tIwe/Zs/Pzzz0hKSkLVqlVdqrV27VrExMTA398fzzzzDCpUqIAlS5agb9++OHr0KEaMGGGz/5w5czB//nyUK1cOtWvXxqVLl4pjiuQtDG0BQkQeIz09XQFQMTExJT2UYjFnzhwFQPXr10+ZzWbr9k8++UQBUAMHDnSpzo0bN1R4eLjy8/NTO3bssG6/cOGCuv/++5Wvr6/KyMiwuc+2bdvUnj171M2bN9WWLVvsOrcR3YmHvolIF5mZmXj++edRs2ZN+Pr64g9/+AOef/55/P7773b7njhxAi+//DIiIiIQEBCA4OBgNGzYEC+++CIuXrxo3e/ixYt488038eCDDyIwMBBBQUGoX78+4uLikJmZWaTxzpo1CwDwzjvvQNM06/YXXngB9erVw6JFi5Cbm3vPOmvWrMHhw4fRv39/NGvWzLq9YsWKePXVV3H9+nXMn2/bwKFly5Zo1KiRtQc5UUF46JtIoI8++gi//fYbmjZtiqZNm6Jx48YoX758SQ/LqYMHD6Jt27Y4ffo0unfvjkaNGmHPnj2YM2cOli9fjk2bNuH+++8HAFy5cgVRUVHIyMhA165d8cQTT+D69etIT0/HggULMHLkSFSsWBFKKcTExCA5ORlRUVGIjY2FyWRCZmYmli1bhr/+9a+oU6eOdQyWsFUuXMPp6tWrSE5OxgMPPGBTw1KnS5cumDFjBlJSUtCuXbsCa61btw4A0LVrV7vbYmJiAADr16+/55iInDE8qO0uxTYyzu2enlJqEBUHpRReeeUV3Lx50xo6mqYhLCwMLVq0QLNmzawBXqtWLZt3gyXlxRdfxOnTpzFjxgwMHDjQuv3jjz/G4MGD8be//Q0///wzAODnn39Geno6hg0bhvfff9+mTnZ2NsqWLQsA2LNnD5KTk9GrVy98++23Nvtdu3YNN27cKPR4Dx8+DLPZjIiICIe3W7anpaXdM6jT0tJs7nOn6tWrIzAw0LoPuUfKWl/SeWFoUN/dgPvsqXNIXp9SpCbeJVVDbxkZGahbty5iYmKwcuXKEhlDcbt27RreffddLFiwAEePHkXlypXRrVs3TJw4Effdd59LNa5evYpXX30VKSkpOHToEM6dO4fg4GCEh4cjISEBf/nLX6wL/Z0uXbqEcePGYcmSJcjKykJoaCj69OmDsWPHIjAwUO+pFommaahevTqOHj1q3aaUQnp6OjIzM7F06VLk5eUBAAIDA9G0aVO0aNHCGt6NGjVCQECAYeP9/fffsXbtWjz44INITEy0ue3FF1/E9OnTsWbNGhw9ehS1atWy3uZojI7+v3C0n5+fH/z8/Gy27du3z+UxWw6vV6xY0eHtQUFBNvsVtZYrdciWlLVeQl4Y+hm1lCbeEhqBlzaWM2zHjh2LkJAQDBs2DG3atMHs2bPRpk0bnD592qU62dnZ+OSTT6BpGh5//HEMHz4cTzzxBI4dO4bnnnsO3bp1g9ls2zs2JycH0dHReP/999GgQQO88soreOCBBzBlyhR06tQJV69eLY4pF0nt2rUdbjebzdaQBvKfj82bN+PTTz9FQkICWrVqhcDAQNSvXx99+/bF2rVri32sqampAIDo6Gi7d/cmkwnt27e32a99+/YIDQ3FpEmT8Pjjj+OTTz7Bb7/9ZnfIumHDhnjooYewePFitG/fHtOmTcP27dvt/v+1aNCgARo0aKDv5KjESFnrJeSFoe+opTTxltAIvLSZP38+Vq1ahX79+mHRokXWBf3TTz/F3/72N7z++uuYMWPGPetUrlwZFy9ehK+vr832mzdvokuXLvjxxx+xYsUKPP7449bbJk+ejNTUVPzjH//ApEmTrNtHjx6Nd999F++//z7GjBmj00z1ERYWhqSkJJtQLsidh4HNZjMOHz6Mw4cPo2bNmujYsWNxDRMArF8tqlatmsPbQ0NDbfarWLEikpKS8Oabb2L58uX44YcfAAC1atXC6NGjMWjQIABAmTJlsGbNGuuREMtXnKpWrYohQ4bgtddeK/TJWJZ3v87e6d45Vj1qVapUqTDDLNWkrPUS8sLQd9RSmnhLaAReVKX1DFuTyWQX0kD+ov7EE08AAA4dOmTdrpTC7NmzERgYiDfeeMPmPm+88QYCAwMxe/bsQs1JD9euXcORI0ewYcMGzJ8/H4MHD0aLFi2waNEil0P6Tpbgio6Oxrfffov33ntP7yHbsRwmPnnypMPbs7KybPYD8o8YzJs3D6dPn8aOHTvw7rvvwmw2Y/DgwVi8eLF1vypVqmD69Ok4duwYfvvtN3z44YeoXLkyxo4di8mTJxd6zPXq1YPJZHL62XFBnzvf7c7Ps++WlZWF7Oxsl+qQLSlrvYS8MDSopTTxltAIvCgOHjyIVq1aYc6cOXj44YcxYsQING/eHHPmzEHLli1x8OBB676WM2ynT5+O8PBwvPTSS4iPj8f999+PBQsWWA85W86wnTBhAipXroyBAwdi4MCBaN68OZYtW2a3CGma5vJJTK6cYZuTk4OUlJRCPydms9n62X7jxo2t29PS0nD8+HFERUXZnTVdvnx5REVF4ciRIzafB+vFWQhbnjtN0+Dv74/w8HBER0cjPj4eH3/8MXbs2OH2Y1lqJSQkYM+ePVi3bh169eplyNd/LF9J2rBhg93ha6UUNmzYYLPfnUwmE5o1a4a///3v1oBetmyZ3X6apqFhw4YYPHgwVq9e7XQ/VwUEBOCRRx7BgQMH7F6EKqWwevVqlC9fHi1btrxnrejoaADAjz/+aHfbqlWrbPYh10lZ60XkhdFf3N6RtFMN6j1cdWnUUw3qPVztSN7lsTX05M7FJTp27KgAqBkzZths/+ijjxQA1alTJ+u2ZcuWKQBq2LBhdnUuX76srl69qpRSateuXQqA6tWrl91+V69eVZcvX7bZBkC5+uezZ88eBUB169bN4e1TpkxRANRnn33mUj2llLp27ZoaO3asevPNN9XgwYNVgwYNFAA1YMAAm/2+//57BUANGTLEYZ0hQ4YoAOrnn392+bGVyn9ODh8+rNavX6/mzZunBg0apJo3b259Xtz5qVWrlnr66afV+++/r1asWKH27NmjFi5ceM/7+fj4KACqTp06atq0aer8+fNuzeFeCvM3OXv2bJvtn376qd3f5J49e1RWVpZdja+//loBUPHx8dbHT09Pt9tv27ZtCoDq0KGDzfZ9+/apffv2uTI1pZT7Fzy5fv262rdvnzp06JDN9hs3bqh69eoVeMETR/Ow4AVPnJOy1pd0Xhj+9axmrR/CR19P9YoaJaG0n2Frcf36dYwfP97635qmYeTIkXjnnXd0e2ylFL744gv88ssv2LJlS6He6daqVQtt2rRBmzZt0KBBA9SqVQt/+MMf7vnZ57Vr15ze5uPjg7y8PHTs2BHDhg3Dn/70J5hMxXdwbPfu3U6vQ92gQQOMHj0an3zyCdq2bYvExEQsX74cDz74IPbu3Ytly5ahatWq+OSTT6z3Wb16NUaNGoWoqCjcf//9qFKlCo4cOYJly5bB398fgwcPBpB/8tmTTz6JRx55BA8++CCqV6+OY8eOYenSpTCZTHjllVdsxtKwYUMArn2PGgDi4uLw5ZdfYvHixUhPT0d0dDQOHTqEb775BnXr1sXEiRNt9j927BgaNmyIOnXqICMjw7q9TJkymD17NmJiYtC+fXubS4hmZmZiypQpCAsLs6m1adMm60culqNamzZtsj7PISEhmDJlikvz8GZS1voSzwtDXxaQU66+e/nuu+8UAPW3v/3N4e0DBw5UANSyZcuUUvmv6kNDQ5Wmaeqxxx5TH3/8sdq7d6/NOwil8t8VPPTQQwqAateunZo6dar69ddfVV5eXpHntnnzZgVA/fnPf3Z4+8yZMxUANW3aNLdr5+XlqaNHj6qPP/5YBQcHq6ioKHXx4kXr7YsWLVIA1Guvvebw/q+++qoCoL755hu727799luX3glPmzbN+k74woULbs/BkZMnT9o9nqZpqly5cmrIkCFq//79ujxOQSx/kwX9REdHW/fPyMhQAwYMUKGhoapMmTIqNDRUDRgwwO7ymb/99pt6+eWXVfPmzVWVKlWUn5+fqlevnoqLi1N79+617nf06FE1evRoFRkZqe677z7l6+urateurZ588km1ZcsWu/FaxuSOq1evqnHjxqnw8HDl6+urqlevrhISEhy+47c8H3Xq1HFYKzk5WcXGxqqgoCAVEBCgHnnkEfXvf//b4b5z584t8Hl19hhUOjGohXA1qBcsWKAAqHHjxjm8fezYsQqAWrhwoXVbZmamiouLU5UrV7YJmY8++sjmvmfOnFFDhgxRoaGh1v2qVq2qxo8fr27evFnouRXHoW9HvvrqKwVA/f3vf7duK8qh77y8PPXFF1/oHsKuyMvLU2XLllWapikAKjw8XE2fPl1dunTJsDEQkQy8hKiHKcoZtmazGbt27cKPP/6IDz74AIMHD0alSpXQr18/ALfPsP3ggw+wf/9+rFmzBtOnT8fYsWNRtmzZQn+FSc8zbAtiuYSj5ZKOd9YszGObTCbrc2M0k8mE4cOHY8+ePXj55ZfRuXNnEVcgI6ISUNKvFCifq++oMzMzFQDVqFEju8PXZrNZNWzYUAFQv//+e4F1NmzYoACoZ555psD9fv/9dwVARUZGujYRJyIjIxUAu8OgZrNZ1atXT5UvX15duXKlSI/x22+/KQAqKirKpn6NGjVUYGCgys7Ottk/OztbBQYGqrp16xbpcYmIihO7Z3mY2rVro2PHjti7dy/mzJljc9vMmTOxb98+dOrUyXoi2d69ex2++7Zs8/f3B5B/CdM7T5Bxtp/F/v37sX//fpfHbbn+85gxY2xO9pkxYwaOHDmCP//5zzYnst24cQP79+/H4cO2FxX47bffcOXKFbv6V65cwfDhwwEAjz32mHW7pmlISEhAdnY2JkyYYHOfCRMmIDs72+6kPCIiSTSlXDxFkoqV5VrfNWrUQJcuXRzuYznD9sCBA2jbti3Onj2LHj162J1he2enon/9619Oz7AFgI0bN6Jly5ZYunSp0zNss7Oz8e2336JHjx7WsbjTqQjI/57zY489hlWrViEyMtLmDNuwsDAkJyejatWqds/H3WfYjhs3DtOmTUPbtm0RFhaGoKAgHDt2DCtWrMDZs2fRrl07rFq1yib0c3JyEBUVhZ07d6Jr165o0aIFtm/fjh9//BGtWrXC+vXrDb02NhGRW0r2DT1Z8Axbx8/H3We/btu2TSUmJqpGjRqp4OBgVaZMGVWlShXVsWNHNWPGDHXjxg2Hj33hwgU1bNgwVatWLVW2bFlVu3ZtNWLECJ6cRUTi8R01ERGRYPyMmoiISDAGNRERkWCGf486NXkXZk2Zj7R9hxHRMByJI+Pcbr4tpQaRJzl32Yy043k4d0mhcpCGiBo+qFwh/7X61d2/4Pz8Sbh+ZC986zVCpbjR8G/yqFt13KlR3PMxsoaedciWlLW+pPPC0M+oU5N3IaH7S9Ym3JYuJLOXT3d50lJqEHmSc5fNWL39Bs5cUgjwBXKvAyFBGrq0KItyGUk4PjQGUAow5wEmH0DTUOODVXZB66xOR/8UZI+OdamGEfNx5QWDHjX0rEO2pKz1EvLC0Jd8s6bMt04WyG++raAwa8p8j6tB5EnSjufhzCWFsGoaqlc2IayahjOXFNKO5+H8/Em3AxbI/61U/nYX65yd53oNID/ckg/cwIpt15F84AbOXc7/t3h19y84MbIHMp8Mx4mRPXB19y9uzSd9w2YcHxqD3JQ1yDtzHLkpa3B8aIzDOnrU0LMO2ZKy1kvIC0ODOm3fYetkLcx5ZqTtO+zkHnJrEHmSc7fe7Vm+/65pGgJ887dfP7L3dsBamPPyt7tYx3TUjRq33oGmHMzDiXNmpBzMw+rtN3Aq2fVgczaOsksnu/yCQY8aetYhW1LWegl5YWhQRzQMtzbftg7Ax4SIhuEeV4PIk1QO0pB7/fYFapRSyL2ev923XqP8Q9V3Mvnkb3exjrmW6zX0eFfubBwBJ11/waBHDT3rkC0pa72EvDA0qBNHxkGDZp205Vh/4qh4j6tB5EkiavggJEhDxkmFrHNmZJxUCAnScH9NH1SKGw1o2u2gvfX5cqU4+yYszuqExLteQ4935c7G4Rfe2OUXDHrU0LMO2ZKy1kvIC8MveGJ39tyoeDR7pIlH1iDyJHefmXx/TR9UCnR21vcY+Ddp41YdV2skH8g/7B1WTYOmaVBKIeOkQqcf+sBv31rbsDb5IKBlJ4ROWebSOALSnZ0Y96PDsehRQ886ZEvKWl/SecErk3kIdSoV2DMbuJAGBEcAjROg3desxOoQFZazs6Q7+afgssMzx90LNndedBRnDT3rUOnGoPYA6lQq8FMiAAUoM6CZAGhA51luhaxedfSgxwsGKTXIfUV9V05UmjCoPYBaMwTISsoPVwvNBFSPhNbpQ0Pr6BaORXzBIKUGEVFx46VzPMGFNNtwBfL/+0KaoXWswZaVBOSezv/9U2L+dnfsmQ1rOFrGAHVru4fVICIqZgxqTxAccevd3h00U/52I+voFWx6vPCQUoOIqJgxqD1B4wQA2u2QtRyibZJobB29gk2PFx5SahARFTMGtQfQ7msGdJ4FVI8EAqrm/+4yG1rVpsbW0SvY9HjhIaUGEVEx48lk5DKnJ18V4kWD3UlpTRI9tgYRUXFiUJNbGGxERMZiUBMREQlm+GfUqcm7MLjPCHRt3AuD+4xAavIuj61BRETFR8paX9J5Yeg7ailNvCU0AicqCXdfESyihg8qV3Dv9bqzGvZXFRsN/yaPFlsNSWPR43klW1LWegl5YWhQD+4zAknrt9n09jT5mBAZ3QoffT3Vo2oQeRpn19ju0qIsymUkuRRKzmp09E9BtsPrdK+yq6NHDUlj0eN5JXtS1noJeWHoSz4pTbwlNAInMpqzPtDpGzbj+NAY5KasQd6Z48hNWYPjQ2NwdfcvLtdwp5e0HjUkjUWP55XsSVnrJeSFoUEtpYm3hEbgREZz1ge67NLJLoeSHr2k9aghaSx6PK9kT8paLyEvDA1qKU28JTQCJzJa5SANudcBy6ddSinkXgcCTroeSs5qmGs1yj88fCeTD3zrNSqWGpLGosfzSvakrPUS8sLQoG7W+iHMXj4dkdGtUDU0BJHRrTD7+w/dasAtpQaRp4mo4YOQIA0ZJxWyzpmRcVIhJEiDX3hjl0PJWY2Q+NGApt2uc+sz3UpxY4qlhqSx6PG8kj0pa72EvOD3qIlKEUd9oAPSk3B8aIyDE6d+dNgLWo9e0nr1o5YyFj2eVyJnGNSliG69pItYQ9pYCG4HJLmGzyvpgUHtIYoaSk6v0915lst19KghciwMeyISjN/I9wDWUMpKAnJP5//+KTF/u6v06CWtVz9qIWPR5XklIipmDGpPoEew6dFLWq9+1FLGotcLDyKiYsSg9gR6hJIevaT16kctZSx6vfAgIipGDGpPoEcoNU4AoN2uY/lMt0misTUkjUWvFx5ERMWIJ5N5AKcnTnWZ7VYvaD16SevVj1rCWPR6XomIihOD2kPoFZBki88rEUnHoCYiIhKsjNEPmJq8C7OmzEfavsOIaBiOxJFxbvf0lFKDqDTytt7L3jYfbyJlrS/pvDD0HbWUJt4SGoETGc1ZINlfPct5z2S9ei/rMZaC6rjznEiaD90mZa2XkBeGBrWUJt4SGoETGclZIHX0T0H26FgH16Ne5TBQkg/cQMrBPIRV06BpGpRSyDipEGnahkpTHnOpjl5j0SNkJc2HbElZ6yXkhaHHd6Q08ZbQCJzISGnH83DmkkJYNQ3VK5sQVk3DmUsKZ+dNcqtnsh69l/Uai7M66Rs24/jQGOSmrEHemePITVmD40NjcHX3L6LnQ7akrPUS8sLQoJbSxFtCI3AiIzkLJNNR93om69F7Wa+x6BGykuZDtqSs9RLywtCgltLEW0IjcCIjOQskc61GbvVM1qP3sl5j0SNkJc2HbElZ6yXkheFfz7I7e25UvNsNuKXUIPIUzj5H7eSfgssOP0d13jO5qL2X9RqLszptlvbBzR1rbMPa5IOAlp0QOmWZ2PmQPSlrfUnnBb9HTVRKOAqkSoGOzkwuXM9kd+roNZaihqy0+RA5wqD2AHr1TNajjqT+zd42H9IPA5K8CYNaOKfXo+48y61A0aOOXmOx1ipCQEqbDxFRceHld6TTq2eyHnV0Gos1ILOSgNzT+b9/SszfbuRY2I+aiDwAg1o6vXom61FHr7HoEZCS5kNEVIwY1NLp1TNZjzp6jUWPgJQ0HyKiYsSglq5xAgDtdqBYPkdtkmh8Hb3GokdASpoPEVEx4slkHkCvnsl61NGthqOTuLrMdquWlPkQERUnBjWVCAYkEZFrGNRERESCGf4ZdWryLgzuMwJdG/fC4D4jkJq8y2NrEBFR8ZGy1pd0Xhj6jlpKE28JjcCJSrO7L7kZUcMHlSs4uuSm4z7S96pT1HHoNR8qPClrvYS8MPQvadaU+dbJAvk9PRUUZk2Z73E1iDzRuctmJB+4gRXbriP5wA2cu5z/b+Dq7l9wYmQPZD4ZjhMjezjs3azH/S01Vm+/gZSDeThxzoyUg3lYvf0GTiW73ke6oDrnLptdGk9R7693HbIlZa2XkBdlDHskyGniLaEROJHR7u7ydCQLyMgyo6N/CrLv6PKUe+4kcn9dixofrLJ5N1vU+1ukHc/DmUsKYdU0aJoGpRQyTiqc/fck+N3dR9rkg/PzJznseuWsTvqGzag05bF7jqeo99e7DtmSstZLyAtD31FLaeItoRE4kdHuDJTqlU0Iq6bhzCWFs/Mm3e40BeT/Vgrn50/S9f4W524FvaZpAPJ/B/gCpqOu95EuqE7ZpZNdGk9R7693HbIlZa2XkBeGBrWUJt4SGoETGa2oAalXwFYO0pB7HbCcHqOUQu51wFyrUX5LyjuZfOBbr5FbdQJOujaeot5f7zpkS8paLyEvDA3qZq0fwuzl0xEZ3QpVQ0MQGd0Ks7//0K0G3FJqEHmaogakXgEbUcMHIUEaMk4qZJ0zI+OkQkiQhpD40YCm3a51q490pbgxbtXxC2/s0niKen+965AtKWu9hLzg96iJSom7P2POvQ6EBGno5J+Cy3d8xmwJyBof/GjTw7mo9797LHeeJX1/TR9UCnR01nfBfaQd1QlIT8LxoTEujaeo99e7DpEjDGpyS1H7SFPJKmpA6hWwxa2o49FrPtKeF/JMDGoPUdSA1CNgnV6ju/OswtUq4fnoWYeIqLgwqD1AUQNSr4BVa4YAWUm2LSo1E1A9ElqnDz1vPjq+8CAiKi68dI4n2DMb1jABbv1Wt7YbcH8LPfpI6zEeveajVx0iomLEoPYERQ1IvQJWjz7SeoxHr/noVYeIqBgxqD1BUQNSr4BtnABAu13Lcqi4SaJ7daTMR686RETFiEHtCYoakDoFrHZfM6DzLKB6JBBQNf93l9nu95EWMh/d6hARFSOeTOYh7M5ObpLoVkAW9f56kzIfac8LEdHdGNRERESCGX7oW0oT75JuBE5ERAWTstaXdF4Y+o5aShNvCY3AiTzV3Vcni6jhg8oVSu/pLnw+ioeUtV5CXhj61ySlibeERuBEJeHcZTOSD9zAim3XkXzgBs5dzv83cHX3LzgxsgcynwzHiZE9cHX3L07vv3r7DaQczMOJc2akHMzD6u35dVytocc4pNTQ6/kge1LWegl5UcawR4KcJt4SGoETGe3uphpHsoCMLDM6+qcg+46mGrnnTiL317Wo8cEq+Dd51KbGnT2pNU2DUgoZJxXSN2xGpSmPuVRDj3FIqaHH80GOSVnrJeSFoe+opTTxltAInMhod4ZK9comhFXTcOaSwtl5k253eALyfyuF8/Mn2dVw1pO67NLJLtfQYxxSaujxfJBjUtZ6CXlhaFBLaeItoRE4kdGchYrp6N7bgWJhzsP1I3vtajjrSR1w0vUaeoxDSg09ng9yTMpaLyEvDA1qKU28JTQCJzKas1Ax12qU3yv5TiYf+NZrZFcjooYPQoI0ZJxUyDpnRsZJhZAgDX7hjV2uocc4pNTQ4/kgx6Ss9RLygt+jJiol7v5MNvc6EBKkoZN/Ci7f8ZksTD6ApqHGBz+63JM6ID0Jx4fGuFRDj3FIqaHH80F0LwxqolLEUahUCjTh6u5fcH7+JFw/she+9RqhUtwYtwPFnRp6jENKDT2eD6KCMKg9hN2lLhsnFKpnclHrSBmHHjX0mgsRUXFiUHsAdSoV+CkR1t7JluYRnWe5H0xFqCNlHJLmQkRU3Hj5HE+wZzasgQLc+q1ubTewjpRx6FFDr7mUgK0HbuA/m66V9DCIyCCGXvCECulC2u1AsVDm/O1G1pEyDj1q6DWXEvDmgqsAgN5t/Up4JERkBL6j9gTBEbd7JltopvztRtaRMg49aug1FyKiYsag9gSNEwBot4PF8nlqk0Rj60gZhx419JoLEVEx48lkHsLuDOUmidCqNjW8jpRx6FFDr7kYLfb1ywCAlRMrlPBIiMgIDGoiD8OgJipdeOibiIhIMMPP+k5N3oVZU+Yjbd9hRDQMR+LIOLebb0upQVTa3H0lr4gaPqhcwf3X+3rVkcLb5iOFlLW+pPPC0EPfqcm7kND9JWsTbksXktnLp7s8aSk1iEpKUQ59OwsU+8tdjr5n/2bLtbG7tCjrcg296xR2LnrXcDafchlJLtchW1LWegl5YWhQD+4zAknrt9k04Tb5mBAZ3QoffT3Vo2oQlZTCBrWzQOnon4Jsh00oVtmESvKBG0g5mIewaho0TYNSChknFVre74Om17c5aUKxyi6Y9KhT1LnoVaOg+USatqHSlMdcrkO2pKz1EvLC0GMzafsO20wWAMx5ZqTtO+xxNYg8TdrxPJy5pBBWTUP1yiaEVdNw5pLC2XmTbocJkP9bKZyfP8nm/s76N5+7dGtfF2roVaeoc9GrRkHzKbt0slt1yJaUtV5CXhga1BENw63Nt60D8DEhomG4x9Ug8jTOAsV0dO/tMLEw5+H6kb02m5z1b64cpOXv60INveoUdS561ShoPgEn3atDtqSs9RLywtCgThwZBw2addKWY/2Jo+I9rgaRp3EWKOZajfIPy97J5APfeo1sNkXU8EFIkIaMkwpZ58zIOKkQcqs1pG8912roVaeoc9GrRkHz8Qtv7FYdsiVlrZeQF4Z/j9ru7LlR8Wj2SBOPrEFUEvT+jLqTfwouO/xM9ke7/skF9W92/NmyfQ096ug1l6LWKGg+AelJbj0nZE/KWl/SecELnhB5GD3P+r4zIG3PTh7jdpjoUcOdOnrMpTifDz2fEyrdGNQewu5yl40T3O6brEcNPepIGYeeYzESr0xGVLowqD2AOpUK/JQIa/9kSwOJzrNcDhU9auhRR8o49ByL0RjURKULL53jCfbMhjVMgFu/1a3tBtbQo46Uceg5FiKiYsSg9gQX0m6HiYUy5283soYedaSMQ8+xEBEVIwa1JwiOuN032UIz5W83soYedaSMQ8+xEBEVIwa1J2icAEC7HSqWz1KbJBpbQ486Usah51iIiIoRTybzEHZnJzdJhFa1qeE19KgjZRx6jsVIPJmMqHRhUBN5GAY1UenCQ99ERESCGR7Uqcm7MLjPCHRt3AuD+4xAavIuj61BRETFR8paX9J5YeihbylNvCU0AicqLB76ptJAylovIS8MfUc9a8p862SB/J6eCgqzpsz3uBpERFR8pKz1EvLC0KCW0sRbQiNwIiJyTspaLyEvDA1qKU28JTQCJyIi56Ss9RLywtCgltLEW0IjcCIick7KWi8hLwz/HrWUJt4l3QicqLB4MhmVFlLW+pLOC17whMjDMKiJShde8ISIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMMODWkoT75JuBE5ERAWTstaXdF4YeglRKU28JTQCJyosXkKUSgMpa72EvDD0HbWUJt4SGoETEZFzUtZ6CXlhaFBLaeItoRE4ERE5J2Wtl5AXhga1lCbeEhqBExGRc1LWegl5YWhQS2niLaEROBEROSdlrZeQF4b3o5bSxLukG4ETFRZPJqPSQspaX9J5YXhQE1HRMKiJShde8ISIiEgwBjV5pXXr1kHTNIwbN66khyLG3sybeH3+Ffx5cjZen38FezNvlmgdInINg9rL5OTk4O2330aLFi0QGBgIPz8//OEPf0C7du0wZswYHD5s+5WCDh06QNM064/JZEKlSpXQrl07zJs3D3d+MpKRkQFN0xAbG+v08S0B+eKLLxbbHEsTs1nh3GWzzY/F3duv33D+KdbezJsY9Vkuth/Kw9lLCtsP5WHUZ7luh6xedfSgxwsGvnghT1CmpAdA+rl8+TLatm2LXbt2oX79+vjLX/6CKlWq4MyZM9i6dSsmTZqE8PBwhIfbf61gxIgRCAwMRF5eHo4cOYJvvvkGmzZtwq+//orp06eXwGwIAL7eeB1zV193eFv/d3Ns/rt5uA/eGVDO4b6L110HFGC+leVmlf8qffG665gY5/oyoFedvZk3sXjddaSfNKNuNRP6dfBFozru3X/UZ7nWsZy/nIfth3Px3vMBLtfRo4aedYic4V+RF/nXv/6FXbt2ISEhATNnzoSmaTa3p6en49q1aw7vO3LkSFSvXt3637t370br1q3x0UcfYfjw4ahbt26xjp0cu+7GG7O7rslgI/2k2RquFmaVv90detTRI9j0eMEg7cULkTM89O1FtmzZAgAYPHiwXUgDQN26ddGgQQOXajVp0gTR0dFQSiElJUXXcd7p6tWrmDp1Kpo2bYqKFSuifPnyCAsLw9NPP42dO3fa7Hvz5k1MmzYNTZs2RUBAACpWrIiOHTti+fLlLj1W/fr1UaFCBVy5csXh7T169ICmaTh48KDN9u+++w5//OMfUalSJfj7+6Nx48aYMmUK8vLybPabN28eNE3DvHnzsHz5ckRFRaFChQoICwtz/Qm5S2QD1xZ6TQMebeh837rVTDDd9Sdh0vK3u0OPOo6CDerWdhfp8YJB0osXooIwqL1IlSpVAMAuaIrKUejrJS4uDiNHjgQADBgwAEOGDMGjjz6KjRs3Ytu2bdb9lFLo3bs3RowYgatXr2Lw4MHo378/du7ciR49euD999+/52P95S9/QXZ2NpYuXWp325kzZ7By5Uq0bt0a999/v3X7mDFj0KtXLxw4cABPPvkkBg0ahICAAIwaNQrPPPOMw8f5+uuv8eSTT+K+++7DoEGD8Kc//cnNZ+W28FATKpa79/OvFPDIA86Dul8HX0CDNWRNGgAN6N/Rz63x6FFHj2DT4wWDpBcvRAXhcRkv0qdPHyxcuBAJCQnYunUrunbtiocfftga4O7Yu3cv1q9fD03T0LJlS5vbDh065PRs6oyMDJcf4+LFi/j666/x8MMPIzk5GT4+Ptbb8vLycPnyZet/L1iwAN999x2io6Px448/wtfXF0B+kD788MP4+9//jp49e6JevXpOH+8vf/kLxo8fj4ULF6J///42t/373//GjRs38Ne//tW6bfXq1Zg0aRJiYmKwZMkSlC9fHkD+i4ZBgwbh008/xZIlS/DUU0/Z1Fq5ciVWrVqFzp07u/xcOGMyaWjzoA9Wb79Z4KHtGpU11KjiPBga1SmD954PsPlcuH9HPzxY28fpfYqrTt1qJpy/nGcT1u4GW78Ovth+OBcm3DrUXIgXDHrU0LMOkTMMai/So0cPTJ06FWPHjsXUqVMxdepUAEB4eDhiY2Px8ssvIyIiwuF9p0yZYj2ZLD09Hd988w1yc3MxdOhQu0O3hw8fxvjx44s8Xk3ToJSCv78/TCbbRdrHxwfBwcHW/54/P79TzeTJk60hDQC1a9fGK6+8gtdeew2LFi3CG2+84fTx6tevjzZt2mD16tU4deoU7rvvPuttCxYsQNmyZdG3b1/rtg8//BAAMHPmTGtIW8Y9adIkzJgxA4sXL7YL6p49e+oS0haRDcpiZYrzD6t9TEBUo7L3rNOoThldPjMtah09gk2PFwySXrwQFcTwoLa7FNvIOLd7ekqpIdHw4cORmJiIlStX4pdffkFKSgqSk5Px0Ucf4bPPPsOXX36JHj162N3PEuqapiEoKAgtW7bE888/j2effdZu35iYGKxcudLh469btw4dO3Z0aaxBQUF47LHH8MMPP6BFixbo06cPOnTogFatWqFsWdvg2bFjB8qVK4dHHnnEro7l8VJTU+/5mH/961+xZcsWLF68GC+//DIAIC0tDVu3bkX37t0REhJi3TcpKQnly5fHnDlzHNYKCAjA/v377bY7GmNRNKvngzI+wM08x7fnmYHWD3hOKOgZkEV94SHlxQs5JmWtL+m8MPQv6+4G3GdPnUPy+pQiNfEuqRqSVahQAX369EGfPn0A5B9ifvXVV/Hxxx/j+eefx7Fjx2zelQLAiRMnbM76NsrXX3+Nt99+G1988QVee+01APkBPmDAALz99tsoVy7/60aXLl1CrVq1HNYIDQ217nMvffv2xbBhw7Bw4UJrUC9YsAAAbA57A8C5c+dw8+bNAo8e5OTk2G2rVq3aPcfhDn9fDc3DffBrWp7dZ7sAUM4PaOhh794YbHQvUtZ6CXlh6NkOUpp4S2gEbqSKFSviww8/RJ06dXDmzBns3r27pIdkVa5cOUycOBFHjhzBkSNH8Nlnn+GBBx7A//3f/+GVV16x7hcUFIRTp045rJGVlWXd514qV66Mxx57DCkpKThw4AAAYOHChahYsSK6d+9us29QUBCqVKkCpZTTn/T0dLvHKI6T7yIblHEY0j4moHWDMvC5+2wmIg8nZa2XkBeGBrWUJt4SGoEbTdM0m89ZJapbty6ee+45rF+/HoGBgVi2bJn1tubNm+PKlSvYunWr3f3WrVsHAGjWrJlLj2N557xw4UJs3rwZ6enp6N27N/z9/W32a926Nc6ePYu0tLTCTUhHrZ18TSvPDLRx8StcRJ5EylovIS8MDWopTbwlNAIvDjNmzLD5StOdli5din379iE4OBiNGzc2eGSOnT59Gnv27LHbfv78eVy7ds0mOOPi4gDkn+V948YN6/ajR49i2rRpKFOmDP785z+79LiPP/44KlWqhEWLFuHzzz8HYH/YGwCGDh0KAHjuuedw9uxZu9uzsrKwb98+lx6zqEKCTA7PijZpwMMRDGryPlLWegl5Yei/8MSRcUhenwKTjwnmPHOhm3hLqCHRihUr8OKLL6J+/fqIiopCjRo1kJOTgx07dmDjxo0wmUz4+OOP4ecn42sjx44dQ/PmzdG0aVM89NBDqFmzJs6ePYvvvvsON27csH6/GsgP0m+++QbfffcdHnroIXTr1g05OTn48ssvce7cOUydOrXAr2bdyc/PD08//TRmzJiBuXPnok6dOmjfvr3dfrGxsXjjjTcwYcIE1K9fH7GxsahTpw7Onj2LQ4cOYePGjZg4cSIaNmyo23NSkEcfLIPMU9eth8A1DWgc5oPy/jzsTd5HylovIi+UwXYk7VSDeg9XXRr1VIN6D1c7knd5bA1p9u/fryZPnqy6dOmi6tatq/z9/ZW/v78KDw9XcXFxKiUlxe4+0dHRCoA6ceLEPeunp6crAComJsbpPmvXrlUA1AsvvHDPeufPn1fjxo1T7du3V6GhocrX11fVqFFDxcbGqhUrVtjtf+PGDTVlyhTVpEkT5efnpypUqKCio6PVd99953QcY8eOdfjYmzZtunU9LKgxY8YUOM7Vq1er7t27q6pVq6qyZcuq6tWrqzZt2qgJEyao33//3brf3LlzFQA1d+7ce869MA7876aKee2S9Sf29Uvqm83XiuWxiCSQstaXdF5oSikHp6gQkTRms0K/STm4eOX2P9k5r5Qv8EInROT5+C+cyENYrlJm+bjsXlcjIyLvwH/lRB4kskFZ5JnzP5925WpkROT5eLookQe58yplJXU1sqL2ktarhqSx6DUfIkf4GTWRh3l9/hXsO5qHr8YEwsfHvTO+ixood/eStlyn251e0nrUkDQWveZD5AwPfRN5mAnPBuA/rxUupEd9lovth/Jw9pLC9kN5GPVZLvZmOm/4cTc9eknrUUPSWPSaD5EzDGoiD6NpWqEuU6pHoOjRS1qPGpLGotd8iJxhUBOVEnoESt1qJtx9WXF3e0nrUUPSWPSaD5Ez/EvycgsXLsQLL7yAli1bws/PD5qmYd68eQXe59KlSxg+fDjq1KkDPz8/hIWFYdSoUcjOzna4v9lsxvTp09GkSRMEBASgatWq6NevH44cOVIMM8p/R9mhQ4diqe3N9AiUfh18AQ3WOoXpJa1HDUlj0Ws+RM7wZDIvFxYWhszMTISEhKB8+fLIzMzE3LlzER8f73D/nJwctG3bFqmpqejatSuaN2+OHTt24Mcff0SrVq2wYcMGu+YViYmJmD17Nho1aoTHH38cx48fx1dffYXAwEAkJSUhIiJC1zlpmobo6GhrMw5yjbOTnqYklHOrF/TdJ6QVppe0HjUkjUWv+RA5ZOh10JSDS7El7fTYGp5g9erVKiMjQyml1DvvvHPPS1y++eabCoD6xz/+YbP9H//4hwKg3n77bZvta9asUQBU+/bt1bVrty9n+cMPPygAqmvXrvpN5hYAKjo6Wve6pcGejBvqtXk5qv+7l9Vr83LU3sybJT0kIqekrPUlnReGBvWOpJ3q4artVYuq7VSzylGqRdV26uGq7d2atJQanuheQW02m1WNGjVUYGCgys7OtrktOztbBQYGqnr16tls79evnwKg1q9fb1evQ4cOCoDKzMx0aXxr1qxRsbGx1ut+33fffapt27ZqxowZSqnb1+929HP3nJYuXao6deqkgoODlZ+fn2rUqJF677331M2btsF05/W5ly5dqlq1aqUCAgJUSEiIGjBggMrKyrIb56+//qqeeuopVatWLeXr66tCQkJUy5Yt1cSJE12aJxHdm5S1XkJeGPoZtZQm3hIagUuUlpaG48ePIyoqyq53dfny5REVFYUjR47g6NGj1u3r1q2z3na3mJgYAMD69evv+dj//e9/8cc//hHJycmIiYnBiBEj0KNHD1y7dg0LFiwAkH8Yf+zYsQCAOnXqYOzYsdafO3tRjxkzBr169cKBAwfw5JNPYtCgQQgICMCoUaPwzDPPOHz8JUuWoE+fPqhfvz6GDRuGJk2aYO7cuWjbti3Onz9v3S81NRWPPvooVqxYgbZt22L48OHo3bs3ypUrh5kzZ95znkTkGilrvYS8MPTb+FKaeEtoBC5RWloaADj9TDkiIgKrVq1CWloaatWqhZycHJw4cQKNGzeGj4/953GWOpa6BZkzZw6UUli7di2aNm1qc5ulF3RYWBjGjRuH8ePHW//33VavXo1JkyYhJiYGS5Yssb7gUEph0KBB+PTTT7FkyRI89dRTNvf7/vvvsXLlSuuLCyA/8CdNmoQ333wT06dPBwAsWLAA165dw9KlS9GzZ0+H4ySiopOy1kvIC0PfUUtp4i2hEbhEFy9eBABUrFjR4e1BQUE2+7m7vysCAgLstlWpUsXl+3/44YcAgJkzZ9ocFdA0DZMmTYKmaVi8eLHd/Tp37mwT0gDw2muvITg4GJ9//jnMZtt/qEUdJxEVTMpaLyEvDA3qxJFx0KBZJ13YJt4SapC+LIekIyMjMWTIEHz77bc4c+aM23WSkpJQvnx5zJkzB+PGjbP5ef/99xEQEID9+/fb3a9du3Z22wIDA9GsWTNcunTJ+lWzp59+GiaTCU888QSee+45LF68GMeOHXN7nERUMClrvYi8MOzT8FukNPEu6UbgJeFeJ5N9//33CoAaMmSIw9uHDBmiAKiff/5ZKZV/ghkA1bhxY4f7/+c//1EA1BtvvOHS+JYuXarat2+vfHx8FAClaZrq1KmT2rFjh81+KOCs7zJlyjg94czyExYWZt3fcjLZp59+6rBe3759FQCVkpJi3bZhwwYVGxur/Pz8rDVbtWql1qxZ49I8icg1Utb6ks4Lw4OaSs69gvrAgQMKgIqJiXF4e0xMjAKgfv/9d+u20NBQVb58ebuzqe98vM8//9ytcV66dEmtWLFCJSYmKpPJpKpVq6bOnz9vvb2goK5cubKqUqWKy49lCerx48c7vN1y5npaWprdbVeuXFFr165Vw4cPV/7+/iogIEAdPnzY5ccmInIFr0xGVhEREahRowY2b96MnJwcm9tycnKwefNm1K1bF7Vq1bJuj46Ott52t1WrVgEA2rdv79Y4KlSogNjYWMycORPx8fE4efIkkpOTrbebTCbk5eU5vG/r1q1x9uxZl05gu9PGjRvttmVnZyM1NRVBQUGoV6+e3e0BAQHo0KEDpk6dildffRW5ublYvXq1W49LRHQvDGqy0jQNCQkJyM7OxoQJE2xumzBhArKzs5GYmGizfeDAgQCAN954A9ev327usGLFCqxbtw5du3ZFnTp17vnYGzZscBi+p06dAgCbq6FVrlwZ//vf/xzWGTp0KADgueeec3gWdlZWFvbt22e3/aeffrK+sLB46623cOHCBTz77LMwmfL/qWzZsgVXr161u//JkyftxklEpAdeQtTLzZ49G5s2bQIA7N69G9u3b0dUVBTq168PAGjbti0SEhKs++fk5CAqKgo7d+5E165d0aJFC2zfvt16CdH169fbnfF89yVET5w4gS+//BKBgYHYsmUL7r///nuOs1mzZjh+/Djatm2LsLAwaJqGTZs2YevWrYiMjMSmTZusXwHr27cvvvrqK/Ts2RPNmzeHj48PevTogYceeggA8Oabb2LChAkIDg5GbGws6tSpg7Nnz+LQoUPYuHEjJk6ciNGjRwMA5s2bhwEDBqBbt25YtWoVevfujbCwMCQlJWHt2rUIDw/Htm3bUKlSJQBAr169sHbtWrRv3x5169aFv78/tm/fjp9//hn16tXDjh07rGe7S1PUXtR61vG2seg1HyJHGNReLj4+HvPnO/9iflxcnF2TjosXL2LcuHFYsmQJsrKyEBoaij59+mDs2LGoUKGCXQ2z2YwPP/wQM2fOxKFDhxAYGIjOnTvjrbfeQni4a19h+PLLL/HNN9/g119/xYkTJ1C2bFmEhYWhf//+GDRoEAIDA637ZmVl4eWXX8aaNWtw7tw5mM1mu+uX//TTT/jggw+QlJSECxcuoEqVKqhbty4ee+wxxMXFWQ/fW4J67ty5qFSpEt566y3s2bMH5cqVQ/fu3fHOO++gevXq1rqrVq3CF198geTkZBw/fhxKKdSuXRu9evXCK6+8gpCQEJfmazRn1/l+7/kAtwJFjzreNha95kPkDIOaSrU7g9pZoxJv8Pr8K9h+KM+mzaVJA1rU98HEuHKG1vG2seg1HyJn+Bk1USmgRy9qvep421j0mg+RMwxqolJAj17UetXxtrHoNR8iZ/iXRFQK9OvgC2iwBorlc9T+Hf0Mr+NtY9FrPkTO8DNqolLi7jOT+3f0w4O17ZupGFHH28ai13yIHGFQExERCWb4oe/U5F0Y3GcEujbuhcF9RiA1eZfH1iAiouIjZa0v6bww9B11avIuJHR/ydqE29KFZPby6WjW+iGPqkFERMVHylovIS8MfUc9a8p862SB/ObbCgqzpji/IIfUGkREVHykrPUS8sLQoE7bd9g6WQtznhlp+w57XA0iIio+UtZ6CXlhaFBHNAy3Nt+2DsDHhIiGrl1mUlINIiIqPlLWegl5YWhQJ46MgwbNOmnLsf7EUfEeV4OIiIqPlLVeQl4Y/vWs1ORdmDVlPtL2HUZEw3AkjopHs0eaeGQNIiIqPlLW+pLOC36PmoiISDBeQpSIiEgwNkslIsPdfcnNfh18Pbp3s7fNh2ThoW+iUqSogaJHIO3NvIlRn+UCKr8dpKWJxXvPBxSqljfNh8gRBjVRKVHUQNErkF6ffwXbD+XZ9HA2aUCL+j6YGFeu1M6HyBl+Rk1USixed90aSsCt3+rWdgPub5F+0mwTapZa6SfNju9QTOORNh8iZxjURKVEUQNFr0CqW81k7d1sYdLyt7vD2+ZD5Az/kohKiaIGil6B1K+DL6DBWstyyLl/Rz+36njbfIicYVATlRJFDRS9AqlRnTJ47/kAtKjvgypBGlrU98GUhHJ4sLaPW3W8bT5EzvBkMqJS5O6znPt39HMrUIp6f71523yIHCn5S4iOjHO7p6eUGkREVHykrPUlnReGBrWUJt4SGoETEZFzUtZ6CXlh6GfUUpp4S2gETkREzklZ6yXkhaFBLaWJt4RG4ERE5JyUtV5CXhga1FKaeEtoBE5ERM5JWesl5IWhQS2libeERuBEROSclLVeQl6U/FnfpbQROBERFUzKWl/SecHvURMREQnGK5MREREJxqAmIiISjF3NiUqRuy+Z2a+Dr1u9l1mjeGoQFYSfUROVEnszb2LUZ7nWHsyWJhTvPR/gcrCwhv41iO6Fh76JSonF665bAwW49Vvd2s4aJVaD6F4Y1ESlRPpJszVQLMwqfztrlFwNonthUBOVEnWrmay9ly1MWv521ii5GkT3wr8molKiXwdfQIM1WCyfp/bv6McaJViD6F54MhlRKXL3Gcr9O/rhwdo+rFHCNYgKUvKXEC2ljcCJiKhgUtb6ks4LQ4NaShNvCY3AiYjIOSlrvYS8MPQzailNvCU0AiciIuekrPUS8sLQoJbSxFtCI3AiInJOylovIS8MDWopTbwlNAInIiLnpKz1EvLC0KCW0sRbQiNwIiJyTspaLyEvSv6s71LaCJyIiAomZa0v6bzg96iJiIgE45XJiIiIBGMfNiJyGXsvExmPh76JSpGiBK2evZf1CHwpNfSsQ+QIg5qolChq0L4+/wq2H8qzaeto0oAW9X0wMa6cYeOQVEPPOkTO8DNqolJi8brr1jABbv1Wt7a7QK/ey0Udh6QaetYhcoZBTVRKFDVo9eq9rEfgS6mhZx0iZxjURKVEUYNWr97LegS+lBp61iFyhn9JRKVEUYO2UZ0yeO/5ALSo74MqQRpa1PfBlIRybvde1iPwpdTQsw6RMzyZjKgUufvs5P4d/dwOWinjkFJDzzpEjpT8JURLaSNwIiIqmJS1vqTzwtCgltLEW0IjcCIick7KWi8hLwz9jFpKE28JjcCJiMg5KWu9hLwwNKilNPGW0AiciIick7LWS8gLQ4NaShNvCY3AiYjIOSlrvYS8MDSopTTxltAInIiInJOy1kvIi5I/67uUNgInIqKCSVnrSzov+D1qIiIiwXhlMiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBDA/q1ORdGNxnBLo27oXBfUYgNXmXx9YgIqLiI2WtL+m8MPQSolKaeEtoBE5ERM5JWesl5IWh76ilNPGW0AiciIick7LWS8gLQ4NaShNvCY3AiYjIOSlrvYS8MDSopTTxltAInIiInJOy1kvIC0ODWkoTbwmNwImIyDkpa72EvDC8H7WUJt4l3QiciIgKJmWtL+m8MDyoiYiIyHW84AkREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwQwP6tTkXRjcZwS6Nu6FwX1GIDV5l8fWICKi4iNlrS/pvDD0Wt+pybuQ0P0laxNuSxeS2cuno1nrhzyqBhERFR8pa72EvDD0HfWsKfOtkwXym28rKMyaMt/jahARUfGRstZLyAtDgzpt32HrZC3MeWak7TvscTWIiKj4SFnrJeSFoUEd0TDc2nzbOgAfEyIahntcDSIiKj5S1noJeWFoUCeOjIMGzTppy7H+xFHxHleDiIiKj5S1XkJeGHoyGZD/wfysKfORtu8wIhqGI3FUPJo90sQjaxARUfGRstaXdF4YHtRERETkOl7whIiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgww4NaShNvvRqBSxmLN9WQNBZvqiFpLN5UQ9JYvKmGpLHoNZ/CMvQSolKaeOvVCFzKWLyphqSxeFMNSWPxphqSxuJNNSSNRa/5FIWh76ilNPHWqxG4lLF4Uw1JY/GmGpLG4k01JI3Fm2pIGote8ykKQ4NaShNvvRqBSxmLN9WQNBZvqiFpLN5UQ9JYvKmGpLHoNZ+iMDSopTTx1qsRuJSxeFMNSWPxphqSxuJNNSSNxZtqSBqLXvMpCkODWkoTb70agUsZizfVkDQWb6ohaSzeVEPSWLyphqSx6DWfojC8H7WUJt56NQKXMhZvqiFpLN5UQ9JYvKmGpLF4Uw1JY9FrPoVleFATERGR63jBEyIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBDA9qKU282dhcbg1JY/GmGpLG4k01JI3Fm2pIGote8yksQy8hKqWJNxuby60haSzeVEPSWLyphqSxeFMNSWPRaz5FYeg7ailNvNnYXG4NSWPxphqSxuJNNSSNxZtqSBqLXvMpCkODWkoTbzY2l1tD0li8qYaksXhTDUlj8aYaksai13yKwtCgltLEm43N5daQNBZvqiFpLN5UQ9JYvKmGpLHoNZ+iMDSopTTxZmNzuTUkjcWbakgaizfVkDQWb6ohaSx6zacoDO9HLaWJNxuby60haSzeVEPSWLyphqSxeFMNSWPRaz6FZXhQExERket4wRMiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwQwPailNvNnYXG4NSWPxphqSxuJNNSSNxZtqSBqLXvMpLEMvISqliTcbm8utIWks3lRD0li8qYaksXhTDUlj0Ws+RWHoO2opTbzZ2FxuDUlj8aYaksbiTTUkjcWbakgai17zKQpDg1pKE282NpdbQ9JYvKmGpLF4Uw1JY/GmGpLGotd8isLQoJbSxJuNzeXWkDQWb6ohaSzeVEPSWLyphqSx6DWfojA0qKU08WZjc7k1JI3Fm2pIGos31ZA0Fm+qIWkses2nKAzvRy2liTcbm8utIWks3lRD0li8qYaksXhTDUlj0Ws+hWV4UBMREZHreMETIiIiwRjUREREgjGoiYiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREgjGoiYiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREgjGoiYiIBPt/NHpvoP7ywO8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "def ploting_A(\n",
    "    optim_final_state, optim_init_state, target_state_final, target_state_init, mask\n",
    "):\n",
    "    # Add two plots side by side\n",
    "    fig = plt.figure(figsize=(5, 5))\n",
    "    turbo_map = cm.get_cmap(\"turbo\")\n",
    "    global_size = 15\n",
    "    plt.scatter(\n",
    "        target_state_init[\"r\"][mask][:, 0],\n",
    "        target_state_init[\"r\"][mask][:, 1],\n",
    "        label=\"target 0\",\n",
    "        s=global_size,\n",
    "        c=turbo_map(0.15),\n",
    "        alpha=0.4,\n",
    "    )\n",
    "    plt.scatter(\n",
    "        target_state_final[\"r\"][mask][:, 0],\n",
    "        target_state_final[\"r\"][mask][:, 1],\n",
    "        label=\"target 100\",\n",
    "        s=global_size,\n",
    "        c=turbo_map(0.15),\n",
    "    )\n",
    "    plt.scatter(\n",
    "        optim_init_state[\"r\"][mask][:, 0],\n",
    "        optim_init_state[\"r\"][mask][:, 1],\n",
    "        label=\"init\",\n",
    "        s=global_size,\n",
    "        c=turbo_map(0.7),\n",
    "    )\n",
    "    plt.scatter(\n",
    "        optim_final_state[\"r\"][mask][:, 0],\n",
    "        optim_final_state[\"r\"][mask][:, 1],\n",
    "        label=\"updated\",\n",
    "        s=global_size,\n",
    "        c=turbo_map(0.85),\n",
    "    )\n",
    "    wall = target_state_init[\"r\"][~mask]\n",
    "    top_lid = wall[:, 1] < 1.17\n",
    "    plt.scatter(\n",
    "        wall[top_lid][:, 0],\n",
    "        wall[top_lid][:, 1],\n",
    "        label=\"wall\",\n",
    "        s=global_size,\n",
    "        color=turbo_map(0),\n",
    "    )\n",
    "    plt.arrow(0.52, 0.96, 0.09, 0.01, head_width=0.03, head_length=0.05, fc=\"k\", ec=\"k\")\n",
    "    plt.text(0.6, 1.03, \"optimizer\\n15 steps\", fontsize=14, ha=\"center\", rotation=0)\n",
    "    plt.text(0.22, 0.94, \"Loss: 0.30\", fontsize=14)\n",
    "    plt.text(0.7, 0.96, \"Loss: 0.01\", fontsize=14)\n",
    "\n",
    "    plt.arrow(\n",
    "        0.63,\n",
    "        0.7,\n",
    "        0.0,\n",
    "        -0.15,\n",
    "        head_width=0.03,\n",
    "        head_length=0.05,\n",
    "        fc=turbo_map(0.15),\n",
    "        ec=turbo_map(0.15),\n",
    "    )\n",
    "    plt.text(0.42, 0.45, \"SPH solver\\n100 steps\", fontsize=14, ha=\"center\", rotation=0)\n",
    "\n",
    "    plt.axis(\"equal\")\n",
    "    plt.axis(\"off\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    fig.savefig(\"optim_A.pdf\", bbox_inches=\"tight\")\n",
    "    print(\"done\")\n",
    "\n",
    "\n",
    "mask = target_state_A[\"tag\"] == Tag.FLUID\n",
    "\n",
    "ploting_A(final_optim_A, init_optim_A, final_state_A, target_state_A, mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_optim_B = read_h5(\"data/optim_B_init.h5\")\n",
    "final_optim_B = read_h5(\"data/optim_B_final.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAHqCAYAAADLbQ06AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABdBklEQVR4nO3deVxU9foH8M8ZlEUQUXBBb4IiJi6Fu+WCmguVqdettBL0qnXVrFzKFlPLykq73mxz+bmUtqrXtDIzRVMTyhS3LFesVNx3wYX5/v6gmRhmBgY4c+aZ4fN+vXxZZ2ae+Z4Bz2fO9n00pZQCERERiWTy9ACIiIjIOQY1ERGRYAxqIiIiwRjUREREgjGoiYiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREgjGoiYiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREghke1OlpOzGi7xh0adgTI/qOQXraTq+tQURE7iNlW+/pvNCUUsqoN0tP24kh9z0GBQVzjhkmPxM0aJi7cibiW97mVTWIiMh9pGzrJeSFoXvUc6YttK4sAJhzzFBQmDNtodfVICIi95GyrZeQF4YG9f69B60ra2HOMWP/3oNeV4OIiNxHyrZeQl4YGtSxcTEw+dm+pcnPhNi4GK+rQURE7iNlWy8hLwwN6qFjk6BBs6605Vj/0HHJXleDiIjcR8q2XkJeGHoxGZB7Yn7OtIXYv/cgYuNiMHRcMuJbNPLKGkRE5D5StvWezgvDg5qIiIhcxwlPiIiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiAQzPKilNPH2dCNwIiIqmJRtvafzwtApRKU08ZbQCJyIiJyTsq2XkBeG7lFLaeItoRE4ERE5J2VbLyEvDA1qKU28JTQCJyIi56Rs6yXkhaFBLaWJt4RG4ERE5JyUbb2EvDA0qKU08ZbQCJyIiJyTsq2XkBeG96OW0sTb043AiYioYFK29Z7OC8ODmoiIiFzHCU+IiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREgjGoiYiIBDM8qKU08fZ0I3AiIiqYlG29p/PC0ClEpTTxltAInIiInJOyrZeQF4buUUtp4i2hETgRETknZVsvIS8MDWopTbwlNAInIiLnpGzrJeSFoUEtpYm3hEbgRETknJRtvYS8MDSopTTxltAInIiInJOyrZeQF4b3o5bSxNvTjcCJPC0jIwO1atVCUlISFixY4Jb3aN++PTZs2AC2vafikLKt93ReGB7URGSc6OhoALmhnB+Dmsg7lPH0AIjIM2rUqIG9e/eiQoUKbnuPDz74AFevXnVbfaLSgHvURD6soD1qIvIOnEKUSID58+ejZcuWCAkJQUhICFq2bGl3OHr9+vXQNA2TJk3Cpk2b0L59e5QvXx5hYWHo3bs3Dhw4YH1uRkYGNE3DkSNHcOTIEWiaZv0zadIkm+ckJyfbvE/79u2haRquXbuGZ599FjVr1kRQUBCaNm2K7777DgBw4cIFjBgxAtWrV0dgYCDuuOMO/Pjjj3brZamVV96xOPqTf70PHz6MIUOGoGbNmggICEBkZCSSk5Nx5MgRu/fTNA3t27fH0aNHMXDgQFSrVg0mkwnr16937QdBJBAPfRN52KhRozBz5kzUqFED//rXvwAAS5cuxaBBg7B9+3b897//tXl+amoqXn31VSQmJuKxxx7Dnj178L///Q8bN25EamoqateujbCwMEycOBEzZswAADzxxBPW17dv396lcd1///3YtWsXunfvjqysLCxevBjdunXD5s2bMWzYMFy/fh19+/bFqVOn8OmnnyIxMRGHDx8u9FD6xIkTHS5/7733cPLkSZQrV866LC0tDV27dsWVK1fQrVs3xMbGIiMjA4sXL8aqVauwZcsW1K5d26bOmTNncMcdd6BSpUp44IEHkJ2djdDQUJfWmUgkRUQes2HDBgVAxcXFqfPnz1uXnz17VtWtW1cBUN9//71SSqmUlBQFQAFQ77//vk2d999/XwFQ3bp1s1keFRWloqKiHL734cOHFQCVlJRkszwhIUEBUG3atFGXL1+2Lv/0008VABUWFqb69u2rbty4YX3stddeUwDU9OnTHdYqzNSpUxUA1aNHD5WTk6OUUur69esqOjpalS9fXm3bts3m+Rs3blR+fn5262v5fAYNGqRu3rxZ6PsSeQMGNZEHDR48WAFQn376qd1jixcvVgDU4MGDlVJ/B3XdunWtYWaRk5OjYmNjlaZp6uTJk9blJQnqDRs22L1H2bJlFQB15MgRm8d+//13BUANHDjQYa2CLF26VGmappo0aWLzxWDZsmUKgHrxxRcdvq5Xr17KZDKpCxcuWJcBUP7+/urUqVMFvieRN+GhbyIP2r59OwDHh6M7dOgAAEhPT7dZ3rp1a5hM+WZKMpnQunVr7N+/Hzt27ECnTp1KPLb4+Hi796hSpQquXr2KmjVr2jwWGRkJADh27FiR3mPr1q14+OGHUb16daxcuRLBwcHWx1JTUwEAv/32m/W8el6ZmZkwm83Yt28fmjVrZl1eq1YtREREFGkcRJIxqIk86OLFizCZTKhcubLdY1WrVoWmabh48aLdckcsyy9cuKDL2Byd1y1TpozT5QBw48YNl+v/8ccfuO+++6BpGlauXInq1avbPH727FkAwOLFiwusc+XKFZv/d/b5EHkrXvVN5EGhoaEwm804deqU3WMnT56EUsouGE+cOOGwlmW5O++L1sulS5fQrVs3nDx5Eh999BEaN25s9xzLeq9cuRIq9zSdwz8JCQk2r8t/lTmRtzM8qKU08fZ0I3AiANaAcnT7kGVZ/kPQmzdvhtmcr5uP2YwffvgBmqbh9ttvty738/NDTk6OrmMuqZycHDzwwAPYuXMn3njjDXTv3t3h81q2bAkA2LJli5HDI0GkbOs9nhdGnhDfnrpDNa3cTjWp3FbFV2qtmlRuq5pWbqe2p+7wuhp6s1zY07VrV4+Nwd2ys7PV5MmTVZ06dVRAQICKjIxUQ4cOVSdOnHC5RlZWlnryySdV27ZtVWRkpAoICFBVq1ZVd955p5o3b566fv26zfOvX7+ulixZogYOHKjq1aungoODVUhIiGrRooV69913PX5lsOWq7/r169tcFHX+/HlVr149m4u6inPVd7NmzVRgYKDKysqye+/CLiZzpKCL0wCohISEQmuNHDlSAVDDhg1zWMciOztb1axZUwUGBtpd2KZU7s9248aNhY6BvJOUbb2EvDD0HLWjBtwmPxPmTFuIdz6f7lU1qGjMZjN69OiB1atXo1WrVujduzf279+PuXPnYu3atUhNTXV4nja/y5cv47333kOLFi1w7733onLlyjh37hxWrVqFwYMH45NPPsGqVausF1sdPHgQffr0QUhICO666y50794dFy5cwMqVKzF8+HB8/fXXWLFihccOl7Zr1w6PPfYYZs6ciYYNG6J3795QSmHp0qX4888/MWrUKLRr187mNV27dsWoUaPw9ddfo0GDBtizZw9WrlyJiIgIu3uuO3bsiK1bt+Luu+9G27Zt4e/vj3bt2tnVNMqPP/6It99+G0FBQahcubLDi8R69uyJ+Ph4BAQEYMmSJbj77ruRkJCAjh07olGjRtaJXDZu3Ijw8HD8+uuvxq8IuZ2Ubb2EvDA0qKU08ZbQCLy0WbhwIVavXo3+/ftj8eLF1mB8//338e9//xvPP/88Zs2aVWidSpUq4cKFC/D397dZfvPmTXTu3BnffvstVq1ahXvvvRcAUL58ebzzzjtISkqyuaJ4+vTpaN++Pb788kssWbIEffv21XFti+att95C48aN8d5772H27NkAgAYNGuDFF1/EoEGD7J7fqlUrPP/883j++efx1ltvwc/PDz179sTrr79uN/nHhAkTcO7cOXz55ZfYuHEjcnJyMHHiRI8FtWXe76ysLLz88ssOnxMdHW093N+8eXPs2LEDb7zxBr7++mts3rwZAQEBqFGjBnr27In+/fsbNXQymJRtvYS8MDSoY+NicObkWZuVLk4Tbwk1PO3IkSN48cUX8c033+DUqVOoUqUKunbtiokTJ9rdOnP8+HFMnToVX3/9Nf7880/rNIwJCQl47bXXrBcfXbhwAdOnT8eSJUvw+++/W2/Had26NV588UVERUUVe7xz5swBALz66qs2e6+PPPII3njjDSxevBgzZsxAUFBQgXVMJpNdSAO5Vx3/85//xPr1622m0qxRowaGDx9u9/zg4GCMHj0aAwYMwIYNGzwa1AAwaNAgh6HsTJs2bVyaFjMkJMQa/vlFR0c77GpVUN2C5gx3pVb79u2L3EmrRo0amDFjhnWWtYIUtTbJJWVbLyEvDL2YTEoTbwmNwEti3759aN68OebNm4emTZtizJgxaNy4MebNm4dmzZph37591udevXoVrVu3xsyZMxETE4PHHnsMycnJqFu3Lj788EPr1cZKKXTt2hUvvfQSKlWqhGHDhmHYsGFo3LgxVqxYgf3799uMwTIvsyuys7ORlpaGW2+91S7sNU1D586dceXKFWzdurXYn4nZbMY333wDAGjYsKFLrylbtiyAv28tIiI5pGzrReSFYWfD/7I9dYca3me06tyghxreZ7TanrbTa2voqSgXk3Xo0EEBULNmzbJZ/s477ygAqmPHjtZlK1asUADUE088YVfn0qVLKjs7Wyml1M6dOxUA1bNnT7vnZWdnq0uXLtksw18XNbli9+7dDi90spg2bZoCoP7v//7PpXpKKXXt2jU1ceJE9cILL6gRI0ZYL7waNGiQyzXuvvtuBUB99dVXLr/GkywXk02cONHTQyEyhJRtvafzglOICuFqUB85csR6lbDZbLZ5LCcnxxpYv//+u1Lq76B+5plnCqxrCer+/fu7NN69e/eqvXv3uvTczZs3KwDqwQcfdPj47NmzFQD15ptvulRPqdwvGZYvCwCUpmlq7NixNvNPF2TWrFl2X2qkY1ATlU485udlLNNJJiQk2B16NplMaNeuHX799Vekp6fjlltuQbt27RAZGYmpU6dix44d6NatGxISEhAXF2fz+ri4ONx22234+OOP8eeff6Jnz55o37494uPj7aarBIB69eq5dT0LExISAqUUzGYzjh07hpUrV+LZZ5/Fli1b8PXXXxfYLenLL7/EyJEjERUVhUWLFhk46pIpzvldIvJ+nJnMy1imk3Q2TaJlzmXL8ypUqIDU1FQMHDgQqampGD58OBo0aICoqCi8++671teVKVMG69atw8iRI3HgwAGMGTMGTZs2RbVq1fDiiy+WaNKMvBerFbROxZlRy2Qy4R//+Af+/e9/Y/bs2di8ebPTq4kB4Ouvv0afPn1QtWpVrFu3zvp5ERFJxaD2MpY9RWfTSGZmZto8DwBq1qyJBQsW4NSpU9i+fTtee+01mM1mjBgxAh9//LH1eeHh4Zg5cyaOHj2KX375BW+//TYqVaqEiRMn4vXXXy/2mGvXrg2TyWR3QZqFZXlsbGyx3wMAunTpAsD5VctfffUVevXqhYiICKSkpNjdykREJBGD2stY7i/9/vvv7Q6DKqXw/fff2zwvL5PJhPj4eDz11FPWgF6xYoXd8zRNQ1xcHEaMGIE1a9Y4fZ6rgoKC0KJFC/z22284cuSI3ZjXrFmD4OBgmw5IxWHp3GS5mjuvr776Cr1790alSpWQkpKCOnXqlOi9iIiMwqD2MjVr1kSHDh2wZ88ezJs3z+ax2bNnY+/evejYsSNuueUWAMCePXsc7n1blgUGBgLIvT/W0T2y+Z9n8euvvxZpRqhhw4YBAJ555hmbLxizZs3CoUOH8OCDD9rcQ33jxg38+uuvOHjQdlKBX375xTppRl5Xr17F6NGjAQD33HOPzWOrVq1C7969UbFiRaSkpJR4z93dFi1ahEceeQTNmjVDQEAANE3DggULnD5/0qRJ1tvlHP0p6N7n4tI0zWFrTiLSHy8mE2bXrl1ITk52+Fi9evUwfvx4vPfee2jTpg2GDh2KlStXon79+tizZw9WrFiBypUr47333rO+Zs2aNRg3bhxat26NunXrIjw8HIcOHcKKFSsQGBiIESNGAMi9SK1Xr15o0aIF6tevj2rVquHo0aNYvnw5TCYTnnzySZuxxMXFAXB9gomkpCR8+umn+Pjjj3H48GEkJCTgwIEDWLZsGWrVqoUpU6bYPP/o0aOIi4tDVFSUTdB89tlnePPNN9GmTRtER0cjNDQUR48exapVq3DmzBm0bdvWZqy//vor/vnPf+LatWto3769zaF+i+joaKefuSc8//zzOHLkCCIiIhAZGWl3FMKZpKQkREdH2y0PCwvTd4BEZCwPXnFOeVhuzyroT95mAxkZGWrQoEEqMjJSlSlTRkVGRqpBgwapjIwMm7q//PKLevzxx1Xjxo1VeHi4CggIULVr11ZJSUlqz5491uf98ccfavz48apVq1aqSpUqyt/fX9WsWVP16tVLbdmyxW68KMJ91BbZ2dlq0qRJKiYmRvn7+6tq1aqpIUOGqMzMTKefR/4GED/99JMaOnSoatCggQoLC1NlypRR4eHhqkOHDmrWrFl2t2flbWThyufqqnPnzqn9+/cX+XWuWLNmjfXn+OqrryoAav78+U6fP3HiRAVApaSkuGU8jhT3cyOiomNQExXDhAkTFAA1b948t76P0UG9bt06lZiYqCIjI5W/v7+qUqWKatOmjXVynYK++OQf4/Lly1XHjh1VWFiYCggIUA0aNFBvvPGGXcey+fPnW1+/fPly1bx5cxUUFKQiIiLUoEGDHH6R+/nnn1Xv3r3VLbfcovz9/VVERIRq1qyZmjJlSok/AyJpeOibqBgGDBiAl156CYMHD0azZs3QqFEjTw8J33//PdLS0mAymRAbG4tOnTohJCTE5dd/9dVXuO+++xAWFoYePXogMjISp06dwo4dO/Dhhx9i2LBhiI6OxsSJEzF58mRERUXZnDLIewHjM888g6lTp6JGjRro1asXKlSogI0bN2LcuHFIS0vD559/bvf+S5cuxerVq9GnTx906tQJqampmD9/PjZu3Igff/wRFStWBJB7mubOO++En58fevTogaioKJw/fx6//PILZs+ejeeee67YnyGRSJ7+pkDkrT7++GPr3uTZs2fd8h5F2aPO/ycsLEwtXLjQ5ffq1auXAqDS09PtHjt9+rTN/6OAQ9/ffvutdZa9y5cvW5ebzWb16KOPKgBqyZIl1uWWPWoA6ptvvrGpNX78eAVAjRw50rps9OjRCoBavnx5oeMk8gWGX/WdnrYTI/qOQZeGPTGi7xikp+302hpUuj3wwAMYOHAggNz2myWZFKYkbr/9dsybNw+HDh1CVlYWDh8+jJkzZ0LTNCQnJxf51jpHHczCw8Ndfv3bb78NIPcuhLytRTVNw9SpU6FpmsOL+jp16oSuXbvaLHvuuecQFhaGDz74AGazbavBko6T5JOyrfd0XmhKGTcnYXraTgy57zFrE25LF5K5K2civuVtXlWDCMjt2hUcHIzs7Gz06dPH4SHdkpg6dSqeeeYZzJ8/v8hXpq9duxadO3dGw4YNsXNn4RuWzz//HP369UPFihUxYMAA3HXXXWjbti0iIiLsnqtpGhISEhxOLlO1alVcuXIFY8eOdfg+b7zxBmrVqoXdu3cDABYsWIBBgwZh8uTJeOGFF+ye36FDB6xfvx779+9HnTp1kJaWhjvvvBOBgYG4//770blzZ7Rr1w41atQodB3Je0jZ1kvIC0P3qOdMW2hdWSC3+baCwpxpC72uBhGQO4mM5V7zJUuWYP78+R4e0d/uuusuxMTEYNeuXdZpWgvSt29fLF++HI0aNcL777+PXr16oUqVKrjrrrusc8y74uzZs7hy5QomT57s8M/Vq1dx5coVu9c5mxbXstwyBW3Lli2xfv16tGvXDh999BEGDBiAf/zjH2jRogVSUlJcHifJJmVbLyEvDA3q/XsP2jTfBnJXev/eg05eIbcGkUVoaCj27t0LABg8eDB27drl4RH9zbI37GiSGEd69OiBDRs24Ny5c1i1ahWGDBmC9evXIzExEefPn3epRmhoKMLDw6Fy7ypx+Ofw4cN2r3M2La5led654Nu2bYtVq1bh3LlzSElJwejRo7Fr1y7ce++9OHTokEvjJNmkbOsl5IWhQR0bF2Ntvm0dgJ8JsXExXleDKK969epZz7vedtttOHfunIdHBFy5cgV79uxBcHCww8PXBSlfvjwSExMxe/ZsJCcn48SJE0hLS7M+bjKZnJ6Tb9myJc6cOeN0bndnNm7caLfs8uXLSE9PR2hoqMO52YOCgtC+fXtMnz4dzz77LLKysqzT3pJ3k7Ktl5AXhgb10LFJ0KBZV9pyrH/ouGSvq0GUnycuLrt06RL27dtntzwrKwtDhw7FpUuX0K9fP5QpU/idmN9//73DMZ88eRKA7TSylSpVwp9//umwzqhRowDkHl04c+aM3eOZmZnWIxB5fffdd1i9erXNspdffhnnz5/HwIEDre1Wt2zZguzsbLvXO5vulryTlG29hLww9GIyIPfE/JxpC7F/70HExsVg6LhkxLco2j2oUmoQ5afHxWVz587Fpk2bAOROKbtt2za0bt3a2kikTZs2GDJkCIDcOdpr166N5s2bIy4uDtWqVcOJEyfw3Xff4c8//0SjRo2QkpLi0tXQ8fHxOHbsmHV6Vk3TsGnTJvz4449o1aoVNm3aBD8/PwDA/fffj88++ww9evRA48aN4efnh+7du+O223IvrnnhhRfw0ksvISwsDImJiYiKisKZM2dw4MABbNy4EVOmTMH48eMB/H0xWbdu3az3UUdHRyM1NRUpKSmIiYnBTz/9ZL2PumfPnkhJSUG7du1Qq1YtBAYGYtu2bVi7di1q166N7du3F9iPnLyHlG29x/PCQ7eFURGZT2xX5rUjlHlpl9y/T2z3aB1y7sKFC9b7goszc1lSUlKBU54mJSXZvNeIESNU8+bNVeXKlVWZMmVU+fLlVYsWLdTrr7+url696vL7fvLJJ6pfv34qJiZGlStXTlWoUEHdfvvt6rXXXlOXLl2yee7x48dVv379VEREhDKZTA7v9V6zZo267777VOXKlVXZsmVVtWrV1B133KFeeukl9fvvv1uf52xmsvDwcJWcnKyOHz9uU/ebb75RAwcOVLfeeqsqX768CgkJUfXr11fPPvusOnXqlOsfNJGXMHyPmopOnUwHvhsKQAHKDGgmABrQaQ60KvGG19GDOpkO7J4LnN8PhMUCDYcUeQxSajjy66+/WhuX7Ny5U8TMZVJZ9qiLcwsaUWnANpfeYPdcWMMV+Otv9ddyY+uok+lQ60ZCLeua+/fJ9KKNAXm+MGSmAlmncv/+bmiRakmp4YzEi8uIyDsxqL3B+f1/h6uFMucuN7CObsGmxxcPKTUK8MADDyApKQmAZ2cuIyLvxqD2BmGxfx2mzkMz5S43so5ewabHFw8pNQoxb94861XIDzzwgG51iaj0YFB7g4ZDAGh/h6zl3HKjocbW0SvY9PjiIaVGISTPXCZFcnIylFI8P03kBIPaC2hV4oFOc4BqrYCgyrl/d54LrfLtxtbRK9j0+OIhpYYLJM9cRkTy8apvcpnTq8aL8aXB7mrrRkO9toarPvnkE/Tv3x9A7nzYlvuCiYgKwqCmIjEy2HxRcnIyFi7Mncz/5s2b1glEiIicYVATGchsNiMkJARZWVluaYtJRL7H8HPUUpp4e7oROJVOJpMJmZmZAHhxGVFhpGzrPZ0Xhu5RS2niLaEROJVueWcu27Fjh3WObCLKJWVbLyEvDN2jltLEW0IjcCrd8s5cdvvtt3PmMqJ8pGzrJeSFoUEtpYm3hEbgRJy5jMg5Kdt6CXlhaFBLaeItoRE4EZA7c1m5cuUAcOYyorykbOsl5IWhQS2libeERuBEQO7FZcePHwfAi8uI8pKyrZeQF4bfniWlibfHG4ET5cGLy4jsSdnWezoveB81kRCcuYyIHGFQlyJ2s4o1HJI7/7fBNaSNRRLOXEZE+TGovURJQ8npPN2d5rhcR48aIsciKOzNZjPKly+Pq1evcuYyIgLA7llewRpKmalA1qncv78bmrvcVXr0ktarH7WQsejyueqMF5cRUX4Mam+gR7Dp0Utar37UUsai1xcPneVvi7lzp+10hUuXLkWdOnW4t01USjCovYEeoaRHL2m9+lFLGYteXzzcwNHMZTk5OXjmmWfQp08fHDx4EJMmTQLPXBH5Pga1N9AjlBoOAaD9XcdyTrfRUGNrSBqLXl883CT/zGX33HMPXnvtNevjv/zyC9LS0jw1PCIyCC8m8wJOL5zqPLdIvaD16CWtVz9qCWPR63N1J7PZjLJly8JsNkPTNJs96DJlyuDBBx/EggULPDdAInI7BrWX0CsgyZb0z3Xp0qV46KGHkJ2d7fBxf39/ZGZm8p5rIh/GoCYSKCcnB88//zymTp1qtyedl6ZpmDFjBkaNGmXwCInIKIafo5bSxNvTjcCJCpKcnGw9H13Qd2mlFN555x1eVEY+Scq23tN5YegetZQm3hIagRMVpHfv3vjf//7ncgBv3LgRbdq0cfOoiIwjZVsvIS8M3aOW0sRbQiNwooJ8+umnmDx5Mvz8/AqdRrRMmTJ47733DBoZkTGkbOsl5IWhQS2libeERuBEBSlTpgwmTJiAtLQ01K5dGyaT83+qN2/exOeff44zZ84YOEIi95KyrZeQF4YGtZQm3hIagRO5omnTptixYwcee+wxAHAa2Dk5OdZmHkWRkZEBTdOQmJhYonFKdu3aNbz44ouIjY1FYGAgqlevjmHDhuHkyZNFrvXTTz/hnnvuQVhYGIKDg9GqVSt89tlnds+7cuUKFi1ahH79+qFu3boICgpCWFgYEhISrBPZUMGkbOsl5IWhQS2libeERuBErgoKCsKMGTOwdu1aVK1a1eGhcLPZzIvKHDCbzejRowcmTpyIiIgIPPHEE7jjjjswd+5c3HHHHTh16pTLtVJSUtC6dWts2rQJ/fr1w6OPPorMzEzcf//9mD59us1zN27ciIcffhjr1q1D48aN8cQTT6B3797YuXMnBgwYgJEjR+q9qj5HyrZeRF4og21P3aGG9xmtOjfooYb3Ga22p+302hpERjt37px66KGHFAClaZoCYPNn3bp1Rap3+PBhBUB17drVTSP2rHnz5ikAqn///spsNluXv/feewqAGjZsmEt1bty4oWJiYlRAQIDavn27dfn58+dV3bp1lb+/v8rIyLAu3759u/rwww/VtWvXbOpkZmaqqKgoBUClpaWVbOVKASnbek/nheFBTUQlt2TJEhUWFqb8/PysIV2mTBnVr1+/ItUpalBnZGSowYMHq+rVq6uyZcuqGjVqqMGDB6sjR47YPffYsWNq1KhRqk6dOiowMFBVqFBB1atXTz3yyCPq/Pnz1uedP39eTZgwQcXFxang4GBVvnx5FRMTowYOHGgTfsVxxx13KAB2dcxms6pdu7YKDg5WV69eLbTO6tWrFQA1aNAgu8cWLFigAKjJkye7NKZXXnlFAVBvvPGGaytBpV4Z4/bdqbj06pmsRx1J/Zt9bX2Konfv3mjdujUGDx6MVatWAci9qGzZsmU4efIkqlSpovt77tu3D23atMGpU6dw3333oUGDBti9ezfmzZuHlStXYtOmTahbty4A4OrVq2jdujUyMjLQpUsX/POf/8T169dx+PBhfPjhhxg7diwqVKgApRS6du2KtLQ0tG7dGomJiTCZTDhy5AhWrFiBhx9+GFFRUdYxaJoGoOB7yy2ys7ORlpaGW2+91aaGpU7nzp0xa9YsbN26FW3bti2w1vr16wEAXbp0sXusa9euAIANGzYUOiYAKFu2LIDcCwaJXMHfFOHs5qPOPgNkpkF1mlOkQNGjjl5jsdYqQUBKWx9PqFatGr766ivMnTsXo0aNwvXr13Hz5k3Mnz8fTz/9tO7v9+ijj+LUqVOYNWsWhg0bZl3+7rvvYsSIEfj3v/+NtWvXAgDWrl2Lw4cP44knnsB//vMfmzqXL1+2htXu3buRlpaGnj174n//+5/N865du4YbN24Ue7wHDx6E2WxGbKzjJiuW5fv37y80qPfv32/zmryqVauGkJAQ63MKkpOTgw8++ACapqFTp06FPp8IYPcs+fTqmaxHHZ3GYg3IzFQg61Tu398NzV1u5FiE9qMuCk3TMHToUOzevRstWrQAkBucZrO5kFcWze+//46UlBTUr18fQ4fadih79NFHUa9ePaxbtw5//PGHzWNBQUF2tUJCQhAQEFDo8wICAhASEmKzbO/evdZe3YW5cOECAKBChQoOHw8NDbV5XklruVJnwoQJ2LVrFwYNGoSGDRsW+nwigHvU8unVM1mPOnqNxVFAaqbc5R3fNm4sgvtRF1VMTAw2bdqEadOm4ebNm9ZDxHpJT08HACQkJNjVNplMaNeuHX799Vekp6fjlltuQbt27RAZGYmpU6dix44d6NatGxISEhAXF2fz+ri4ONx22234+OOP8eeff6Jnz55o37494uPjHd6KVq9ePV3Xy0jvv/8+Xn31VTRu3Bj//e9/PT0c8iLco5ZOr57JetTRayx6BKSk9RHCz88PTz/9NJ577jndg/rixYsAgKpVqzp8PDIy0uZ5FSpUQGpqKgYOHIjU1FQMHz4cDRo0QFRUFN59913r68qUKYN169Zh5MiROHDgAMaMGYOmTZuiWrVqePHFF5GTk1PsMVv2fp3t6eYdqx61Cqozd+5cDB8+HI0aNcKaNWvsjhQQFYRBLV3DIQC0vwPF0jO50dCCXuWeOnqNRY+AlLQ+pYDlMPGJEyccPp6ZmWnzPACoWbMmFixYgFOnTmH79u147bXXYDabMWLECJtJP8LDwzFz5kwcPXoUv/zyC95++21UqlQJEydOxOuvv17sMVtmdHN27rig88755T2fnV9mZiYuX77stM6cOXMwbNgw1K9fH2vXrkV4eLirq0AEgEEtnlYlHug0B6jWCgiqnPt357lF7pmsRx29xqJHQIpan1IgPj4eAPD999/bXXGtlML3339v87y8TCYT4uPj8dRTT1kDesWKFXbP0zQNcXFxGDFiBNasWeP0ea4KCgpCixYt8Ntvv+HIkSN2Y16zZg2Cg4PRrFmzQmslJCQAAL799lu7x1avXm3znLzmzJmDRx55BHFxcVi3bh0qV65cnFWh0s6jN4dRqWU+sV2Z145Q5qVdcv8+me7pIZVKRbmPukOHDgqAmjt3rs3y999/XwFQHTt2tC7bvXu3yszMtKvx+eefKwAqOTnZ+v6HDx+2e95PP/2kAKj27dvbLN+7d6/au3evK6umlCr6hCfXr19Xe/fuVQcOHLBZfuPGDVW7du0CJzzJvx5z5sxRmqapuLg4h58FkasMbXNJRLJkZGSgVq1aqF69Ojp37uzwOfXq1cP48ePx22+/oU2bNjhz5gy6d++O+vXrY8+ePVixYgUqV65scx/1jBkzMG7cOLRu3Rp169ZFeHg4Dh06ZN1D3rhxI5o1a4bly5ejV69eaNGiBerXr49q1arh6NGjWL58OS5fvoz//e9/6N69u3UsRbmPGsidQvSee+7B6tWr0apVKyQkJODAgQNYtmwZoqOjkZaWZrOXa/k8oqKikJGRYVMrJSUFXbt2RWBgIB544AGUL18eS5cuxZEjRzBt2jSMGTPG+tx169ahU6dOUErhkUceQbVq1ezGFh8fj549e7q0HlTKGf3NwG4qttQdXluDyNtZ9qgL+pOQkGB9fkZGhho0aJCKjIxUZcqUUZGRkWrQoEF2M3/98ssv6vHHH1eNGzdW4eHhKiAgQNWuXVslJSWpPXv2WJ/3xx9/qPHjx6tWrVqpKlWqKH9/f1WzZk3Vq1cvtWXLFrvxWsZUFNnZ2WrSpEkqJiZG+fv7q2rVqqkhQ4Y43Mu1fB5RUVEOa6WlpanExEQVGhqqgoKCVIsWLdQnn3xi97z58+cX+rkmJSUVaT1KIynbek/nhaF71FKaeEtoBE5ERM5J2dZLyAtDLyaT0sRbQiNwIiJyTsq2XkJeGBrUUpp4S2gETkREzknZ1kvIC0ODWkoTbwmNwImIyDkp23oJeWFoUEtp4i2iETgRETklZVsvIS8Mvz0rPW0n5kxbiP17DyI2LgZDxyUjvkUjr6xBRETuI2Vb7+m84H3UREREgrF7FhVJSftIE+khe9cPOLdwKq4f2gP/2g1QMWk8Ahvd6elhEbkF96i9REkDUo+AtfaRtrSotMzR3WlO8Wp5eH30rFNa6BGQJa2RvesHHBvVFVAKMOcAJj9A01D9rdUMa/JJDGovUNKA1Ctg1bqRQGaqbYtKzQRUawXN1T7SktZHxy8epYEeAalHjeNjuyNr67rc11uY/BDUrCMip7nexIN75eQt2D3LG+yeC2uYAH/9rf5absDrLfToI63HePRaH73qCJW96wccH9sdR3rF4PjY7sje9UOJ6p1bOPXvgAVy/1Yqd7mBNa4f2mMb0n/VuX5oj8s1LF8YsrauQ87pY8jaug7HRnUt8Wek92dOBDCovUNJA1KvgNWjj7Qe49FrffSqI5A7gkiPgNSjhn/tBrl74nmZ/HKXu0iPLwz5uSv8iRjU3qCkAalXwOrQR1qX8ei1PnrVEcgdQaRHQOpRo2LSeEDT/q7z1+HziknPuFxDjy8M+bnjMycCGNTeoaQBqVPAalXigU5zgGqtgKDKuX93ngut8u1FqiNlfXSrI5A7gkiPgNSjRmCjO1H9rdUIatYRfhHVEdSsI6q/9S0CG93hcg09vjDk547PnAjgxWRew+7q5EZDixSQJX293qSsj7TPRS96XXCVn/0FWM8UKSD1qlFSzi9qK1rg5+Wuz5yIQU3kg9wRRL5G7y8M/MzJXTw/hejYpCL39JRSg0gyCXuupQ0/c31J2dZ7Oi8MDWopTbwlNAInMtrZS2bsP5aDsxcVKoVqiK3uh0rli36Zil51SkqPcfjaZ+JLpGzrJeSFob9JUpp4S2gETmSks5fMWLPtBrbuy8Hxs2Zs3ZeDNdtu4Owlc+EvdlOdtN9uYNVP15H2m2fGIe0zIVtStvUS8sLQoJbSxFtCI3AiI+0/loPTFxWiq2qoVsmE6KoaTl9U2H8sp/AX61xHj2DTYxySPhOyJ2VbLyEvDA1qKU28JTQCJzLS2YsKQf6ApmkAcv8O8s9dbnQdXcJeh3FI+kzInpRtvYS8MDSopTTxltAInMhIlUI1ZF0HLJekKKWQdT13udF19Ag2PcYh6TMhe1K29RLywvNXfZfSRuBERrIcbj79V0hmXQciQjV0aVoWFUNc/76uR52033IPe0dX1aBpGpRSyDih0KyuH1reWtawcUj6TMgxKdt6T+cF76MmKiXyX5lct4ZfsYKkOHXy3rakRTXArsZjkFGxZYkDcv+xHGTt2IKqKW8gMHMPAmMaFqkLlic/EyJXMaiJyK0cTgQC4Ny4VThZpWWJgo29qak0KOPpAZBr7Ka6bDikWD2TS1pHyjj0qKHXuvgavfs0O2xWYfJD9ZQ30LSEU2s6q31u4VSdp0plr2ryHO5RewF1Mh34biisvZMtzSM6zSl6MJWgjpRxSFoXX+OOPdQjvWKQc/qY3XK/iOqIWlayW1zcUZt76SQNT6J4g91zYQ0U4K+/1V/LDawjZRx61NBrXXyM1PaYRtZmu0qShkHtDc7v/ztQLJQ5d7mRdaSMQ48aeq2LANm7fsDxsd1xpFcMjo/tjuxdPxS7ltT2mEbWdsdnoOfPiEofBrU3CIv9u2eyhWbKXW5kHSnj0KOGXuviYZbDtFlb1yHn9DFkbV2HY6O6FjsI3LGHqkf/aCNr6/0Z6P0zotKHQe0NGg4BoP0dLJbzqY2GGltHyjj0qKHXuniY3odp3bX3G9joTkROW4GoZQcROW2Frh2l9K6t92fAQ+lUUryYzEvYXaHcaCi0yrcbXkfKOPSoode6eJK7LqYq7a0a9fwM3HkxHZUODGoiL3Z8bHdkbV1ne07V5IegZh1LdHsS6Yc/IyopHvom8mLuvFCL9MGfEZWU5+f6HptU5ObbUmoQSeCNh6rzT7kZW90PlcqXfOrO4tTRaywF8cafkQRStvWezgtDgzo9bSeG3PeYtQm3pQvJ3JUzXV5pKTWIvJEngi1/SPn1fRop2c3smlh0blK2SGPJ2wyj6olU3LJxGiqc/gXl6jZElcGuzSTmrKFGccbi7rAvbaRs6yXkhaG/SXOmLbSuLJDbfFtBYc60hV5Xg8jbWEJp674cHD9rxtZ9OViz7QbOXsr9d+DKvb6F1cjP0a1Jl57uipu/bClRP2rg777WcVfT0HhRN0RkrIf/peO4sd31258K6o3t6r3PRf1MyDVStvUS8sLQoN6/96B1ZS3MOWbs3+v6lY9SahB5m8JCyZV7fQuq4YjjW5OAumnTStSPGvi7r3XEd9MABWgq9z20Itz+5Kw3dtaOLS7f+1zUz4RcI2VbLyEvDA3q2LgYa/Nt6wD8TIiNi/G6GkTexlkonb2oXL7Xt6Aajjia5UtTOQg++QssZ92UUsi6DlQK1Yq0PpVCNWRdBwKO77GGtJWLM4lZauQfS9WUN1y+97monwm5Rsq2XkJeGBrUQ8cmQYNmXWnLsf6h45K9rgaRt3EWSpVCNZenzSyohiPOZvnKrlYfGScUMs+akXFCIeKvVpdFEVvdDxGhGs6H14fSijeTmKVG/rEEZro+jWhRPxNyjZRtvYS8MDSo41vehrkrZ6JVQnNUjoxAq4TmmPvl24hv0cjrahB5G2ehVLeGn8vTZhZUwxFntyb945Fn0ayuHyIrmdCsrh+6NC1b5H7Ulcqb0LlJWZTp9zSgAaoYtz9ZauQfS2BMQ5enES3qZ0KukbKtl5AXnPCEyAvo1R85/9XJdWv4oWKIqYDWjvbzZjur4frY9b81Se/3KMrnART9M3E+bva9JnsMai9hN91lwyFF7pusRw096kgZh55jcSej+iPzXl9b7v482PeaXMWg9gLqZDrw3VBY+ydbGkh0muNyqOhRQ486Usah51jcjVNQ+ib+XMlVvCPfG+yeC2uYAH/9rf5abmANPepIGYeeY3Ezd/RHJs/jz5VcxaD2Buf3/x0mFsqcu9zIGnrUkTIOPcfiZu7oEU2ex58ruYpB7Q3CYv/um2yhmXKXG1lDjzpSxqHnWNyMTR18E3+u5CoGtTdoOASA9neoWM6lNhpqbA096kgZh55jcbPARnei+lurEdSsI/wiqiOoWUenVx+T9+DPlVzFi8m8hN3VyY2GQqt8u+E19KgjZRx6joWIyF0Y1ETkNrxPmKjkGNRE5BZG3yfMLwXkqwwPailNvD3dCJzIW7nae9nI+4RL8qXAEz26yTVStvWezgtDg1pKE28JjcCJjKZXIK3ZdgOn/+oYlXUdiAjV0LlJWbtaR3rFIOf0MbsafhHVUX7hfl2CzbJOQdN6IXh/Sm6LSwsXvhQUZX3cWYPsSdnWS8gLQ3+LpDTxltAInMhIljDZui8Hx8+asXVfDtZsu4Gzl8yFvziPovRednafsBbVQJex5F2nskf32IY04NLkIXr0kmY/aveQsq2XkBeGBrWUJt4SGoETGUmvMClK72Vn9wmf6vSULmPJu043ajQoVqvLkvSSzt71A46P7Y6wMbeixdI+KJeRVuQa5JyUbb2EvDA0qKU08ZbQCJwoL8tG/0ivGBwf2x3Zu37Qtf7ZiwpVT6Tiljl9EDO5Hm6Z0wdVT6QWOUyK0nvZ2X3CJ6u0KHY45l8nS50zncfltrrUijZ5SHF7SVvOiWdtXYeyF4+j4uH1qPnuPQg6nOq2ftTu/h2RRsq2XkJeGBrUUpp4S2gETmSRd6Ofc/oYsrauw7FRXXXdEFc5mYb4D7sheF8Kyl48juB9KYj/sBuqnEwrUp2i9l4ObHQnIqetQNSyg4ictgKBje4odjjml7dOVq1WOPLvr3A6uj1ywiJdnjykuL2kzy2c+veFawBMKgdKASFfv+6WftRG/I5II2VbLyEvPH/V97jkIjfgllKDSA9GXB39+5PdcWPbOmjq7/dQmh/KNumImv8p2nsUt/dy3tc7uviqS9OyHqtT1PVxdqHcjdBInJ/+W5E/k8KU1k5bUrb1ns4L3kdN5GEFXR0dtUyf82BGvEdRlDTsXanjzvuqjQ5OaT8/MlYZTw+AqLTzr90AWWdP2G309eyi5O73KGooVipvQstbS77H6axO/vuqs86eQNbPKbpNtlIxaTyyfk7JPRee575tdzXUMOJ3hOTiTX5EHmZEFyV3vofE86f5zyHDnAMolbu8mPJezHVu4VSEPz7dsIYa7LRVuvHQN5EA9nukz+i+0Xf0HoAq8eFhiedP9T5UbPR0qM7G4O7fEZKJQU1USukVPhLPn+r95UHilxEqPXjom6iU0uvwsLMZyDx5/lTvQ8XXD+2xDWnApZnPiPTAoCYqpfQKH4nnT51NtlLcQ8USv4xQ6cFD30SllJ6Hc339/Knz0wTuu4CMyIJBTVRKMXyKRsKXEfbcLp0Y1EQeImGjKyF8yDUSrjwnz/D8FKKltBE4lW6lfaOrR2/s0qY0XnkuZVvv6bwwNKilNPGW0AicSjdPbXTdHZCuHCVwNkd35yZlizwWZ+vj6tEKvT4PI754SLwNzp2kbOsl5IWhX2GlNPGW0AicSjdP3O5jCcit+3Jw/KwZW/flYM22Gzh7yVz4i/PUSPvtBlb9dB1pv9m+1tUZynTrje1kfU6mbXZpHHp8HoWNQ8+2lKXtynMp23oJeWFoUEtp4i2hETiVbu7a6BbUs7ikAVlYsLl6X3bePtJA8ftRO1ufMwtcG4deXxgc1bn5yxZcelrfaVVLchucN/aylrKtl5AXhga1lCbeEhqBk28p6obQHfceF7ZHW9KALCzYXD1KoFc/amfrY/rDtXHo9YXBUZ26adMABV3nGi/uveES52J3hZRtvYS8MDSopTTxltAInHxHcTaEek/IARS+R1vSgCws2Fw9ShBb3Q8RoRoyTihknjUj44RCxF8tKovC2fqYb3FtHHp9YXBUJ/jkLza9vwHocmojsNGdiJy2AlHLDiJy2gqXfl/c0aDECFK29RLywtCgjm95G+aunIlWCc1ROTICrRKaY+6XbxepAbeUGkQWxd0QFmejW5DC9mhLGpCFBZurRwkqlTehc5OyaFbXD5GVTGhW1w9dmpYtcj9qZ+sTkezaOPT6wuCoTna1+mLOJ3vr9KdStvUS8oL3UROVkJSrcV25kjz/1cl1a/i5HJDOrtbOG7JG35ftbH1cHUdJPo+C6tQ+/yMujU8UMZlMabyty9cwqIlKSMqG0IiZxvQKttJAymQynIHO+zGoiUpI0oZQSjiQLPy98G4MaiIdcENIRO7CoCYiIhKMJ5eIiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMMODOj1tJ0b0HYMuDXtiRN8xSE/b6bU1iIjIfaRs6z2dF4beniWlibeERuBEnpB/ZrHY6n6oVL5o39dZQ/8aZE/Ktl5CXhj62ySlibeERuBERiusnzRreKYGOSZlWy8hLwwNailNvCU0AicyWmH9pFnDMzXIMSnbegl5YWhQS2niLaEROJHRCusnzRqeqUGOSdnWS8gLQ4NaShNvCY3AiYxWWD9p1vBMDXJMyrZeQl4YPtd3etpOzJm2EPv3HkRsXAyGjksucgNuKTWIvIkr/aRZw/ga5JyUbb2n84JNOYhKET36SbOG/jWICsKgJiIiEoxf+4iIiARjUBMREQnGoCYSICdHYfy8q1iz7Yanh0JEwjCoiQS4eh1IP5SDGcuzPT0UIhKGQU0kCK/sJKL8GNRERESCMaiJiIgEY1CT11q/fj00TcOkSZM8PRQiIrcxPKilNPH2dCPwkrhy5QpeeeUVNGnSBCEhIQgICMA//vEPtG3bFs888wwOHrTt6tK+fXtommb9YzKZULFiRbRt2xYLFixA3jlvMjIyoGkaEhMTnb6/JSAfffRRt60jEZGUbb2n86KMkW+WvwH3mZNnkbZha4maeHuqhqdcunQJbdq0wc6dO1GnTh089NBDCA8Px+nTp/Hjjz9i6tSpiImJQUyMfWeXMWPGICQkBDk5OTh06BCWLVuGTZs24eeff8bMmTM9sDZktPzTXcZW90Ol8q5/Xy/p6/WsI6WGnnXob1K29RLywtCgdtSA2+RnwpxpC/HO59O9qoanzJgxAzt37sSQIUMwe/Zsa3s9i8OHD+PatWsOXzt27FhUq1bN+v+7du1Cy5Yt8c4772D06NGoVauWW8dOnpW/gcShTCAj04zOTcq6FColfb2edaTU0LMO2ZKyrZeQF4b+Fklp4i2hEXhxbdmyBQAwYsQIu5AGgFq1aqFevXou1WrUqBESEhKglMLWrVt1HWde2dnZmD59Om6//XZUqFABwcHBiI6ORr9+/bBjxw6b5968eRNvvvkmbr/9dgQFBaFChQro0KEDVq5c6dJ71alTB+XLl8fVq1cdPt69e3domoZ9+/bZLP/iiy9w1113oWLFiggMDETDhg0xbdo05OTk2DxvwYIF0DQNCxYswMqVK9G6dWuUL18e0dHRrn8gHrL/WA5OX1SIrqqhWiUToqtqOH1RYf+xnMJfrMPr9awjpYaedciWlG29hLwwNKilNPGW0Ai8uMLDwwHALmhKylHo6yUpKQljx44FAAwaNAgjR47EnXfeiY0bN+Knn36yPk8phT59+mDMmDHIzs7GiBEjMGDAAOzYsQPdu3fHf/7zn0Lf66GHHsLly5exfPlyu8dOnz6Nb775Bi1btkTdunWty5955hn07NkTv/32G3r16oXhw4cjKCgI48aNwwMPPODwfT7//HP06tULVapUwfDhw3H33XcX8VMx3tm/9vgsP2tN0xDkn7vciNfrWUdKDT3rkC0p23oJeWHooe+hY5OQtmErTH4m6+GD4jTxllDDU/r27YtFixZhyJAh+PHHH9GlSxc0bdrUGuBFsWfPHmzYsAGapqFZs2Y2jx04cMDp1dQZGRkuv8eFCxfw+eefo2nTpkhLS4Ofn5/1sZycHFy6dMn6/x9++CG++OILJCQk4Ntvv4W/vz+A3CBt2rQpnnrqKfTo0QO1a9d2+n4PPfQQJk+ejEWLFmHAgAE2j33yySe4ceMGHn74YeuyNWvWYOrUqejatSuWLl2K4OBgALlfGoYPH473338fS5cuRe/evW1qffPNN1i9ejU6derk8mfhaZVCNRzKzF03TdOglELW9dzlRrxezzpSauhZh2xJ2dZLyAtDgzq+5W2Yu3JmiRpwS6nhKd27d8f06dMxceJETJ8+HdOn554jiYmJQWJiIh5//HHExsY6fO20adOsF5MdPnwYy5YtQ1ZWFkaNGmV36PbgwYOYPHlyicdr2XAFBgbCZLL9Vurn54ewsDDr/y9cuBAA8Prrr1tDGgBq1qyJJ598Es899xwWL16MCRMmOH2/OnXq4I477sCaNWtw8uRJVKlSxfrYhx9+iLJly+L++++3Lnv77bcBALNnz7aGtGXcU6dOxaxZs/Dxxx/bBXWPHj28KqQBILa6HzIyzcg4oRDknxsmEX/1Tzbi9XrWkVJDzzpkS8q2XkReKPJKFy9eVJ999pl64oknVJs2bVTZsmUVABUYGKi++OILm+cmJCQo5M5OqQAoTdNUhQoVVNu2bdWCBQuU2Wy2Pvfw4cMKgOratavT905JSVEA1COPPOLSWO+55x4FQMXHx6uXX35Zbd68WV2/ft3ueRUrVlTlypVzWOOHH35QAFSvXr3sxjFx4kSb57777rsKgJoxY4Z12b59+xQAdd9999k8t0qVKio4OFhNnDjR4Z9y5cqpBg0aWJ8/f/58BUC9+uqrLq27qy5eNauuz11Ud0+4qGvd/M5czFGpv15XX/94TaX+el2dvZRj6Ov1rCOlhp51iBwxdI+a9FO+fHn07dsXffv2BZB7iPnZZ5/Fu+++i3/96184evSozV4pABw/ftzmqm+jfP7553jllVfw0Ucf4bnnngMAhIaGYtCgQXjllVdQrlw5AMDFixdxyy23OKwRGRlpfU5h7r//fjzxxBNYtGgRHn/8cQC5e9MAbA57A8DZs2dx8+bNAo8eXLlyxW5Z1apVCx2HRJXKm9Dy1uJfmlLS1+tZR0oNPesQOcLfLB9RoUIFvP3224iKisLp06exa9cuTw/Jqly5cpgyZQoOHTqEQ4cO4f/+7/9w66234r///S+efPJJ6/NCQ0Nx8uRJhzUyMzOtzylMpUqVcM8992Dr1q347bffAACLFi1ChQoVcN9999k8NzQ0FOHh4VBKOf1z+PBhu/dw58V3RER5Mah9iKZpNudZJapVqxYGDx6MDRs2ICQkBCtWrLA+1rhxY1y9ehU//vij3evWr18PAIiPj3fpfSx7zosWLcLmzZtx+PBh9OnTB4GBgTbPa9myJc6cOYP9+/cXb4WIiNyMQe1lZs2aZXNLU17Lly/H3r17ERYWhoYNGxo8MsdOnTqF3bt32y0/d+4crl27ZhOcSUlJAHKv8r5x44Z1+R9//IE333wTZcqUwYMPPujS+957772oWLEiFi9ejA8++ACA/WFvABg1ahQAYPDgwThz5ozd45mZmdi7d69L70lE5A48R+1lVq1ahUcffRR16tRB69atUb16dVy5cgXbt2/Hxo0bYTKZ8O677yIgIMDTQwUAHD16FI0bN8btt9+O2267DTVq1MCZM2fwxRdf4MaNG9b7q4HcIF22bBm++OIL3HbbbejWrRuuXLmCTz/9FGfPnsX06dMLvDUrr4CAAPTr1w+zZs3C/PnzERUVhXbt2tk9LzExERMmTMBLL72EOnXqIDExEVFRUThz5gwOHDiAjRs3YsqUKYiLi9PtMyEiKgoGtZd57bXX0Lp1a6xZswbff/89jh8/DgCoUaMGkpKS8Nhjj6Fp06YeHuXfoqOjMWnSJKxbtw7fffcdzpw5g4iICDRp0gSPP/64TfMPTdOwZMkS/Pe//8XChQsxc+ZM+Pv7o0mTJhg9ejS6d+9epPd++OGHMWvWLNy4cQMDBgxwel75xRdfRLt27fDWW29h7dq1OH/+PMLDw1GrVi1MmjTJ5b14IiJ30JRSnD6HyMMuZSn0ffkyTCbg6xfLe3o4RCQIz1ETEREJxqAmIiISjEFNREQkmOFBnZ62EyP6jkGXhj0xou8YpKft9NoaRETkPlK29Z7OC0MvJktP24kh9z1mbcJt6UIyd+VMxLe8zatqEOmJF5MR2ZKyrZeQF4buUc+ZttC6skBu820FhTnTFnpdDSIich8p23oJeWFoUO/fe9C6shbmHDP27z3odTWIiMh9pGzrJeSFoUEdGxcDk5/tW5r8TIiNi/G6GpIsWrQIjzzyCJo1a4aAgABomoYFCxYU+JqLFy9i9OjRiIqKQkBAAKKjozFu3DhcvnzZ4fPNZjNmzpyJRo0aISgoCJUrV0b//v1x6NAhN6xR7uQn7du3d0ttIpJPyrZeQl4YGtRDxyZBg2Zdacux/qHjkr2uhiTPP/88Zs+ejSNHjljbQRbkypUrSEhIwH/+8x/Uq1cPTz75JG699VZMmzYNHTt2RHZ2tt1rHnnkEYwaNQpKKYwaNQqJiYlYtmwZmjdvzoYWRKQ7Kdt6EXlhdAPs7ak71PA+o1XnBj3U8D6j1fa0nV5bQ4o1a9aojIwMpZRSr776qgKg5s+f7/T5L7zwggKgnn76aZvlTz/9tAKgXnnlFZvl69atUwBUu3bt1LVr16zLv/76awVAdenSRb+V+QsAlZCQoHtdqS5eNauuz11Ud0+46OmhEIkhZVvv6bwwPKjJvQoLarPZrKpXr65CQkLU5cuXbR67fPmyCgkJUbVr17ZZ3r9/fwVAbdiwwa5e+/btFQB15MgRl8a3bt06lZiYqCIjI5W/v7+qUqWKatOmjZo1a5ZSSqmUlBQFwOGf/Ou0fPly1bFjRxUWFqYCAgJUgwYN1BtvvKFu3rxp87z58+dbX798+XLVvHlzFRQUpCIiItSgQYNUZmam3Th//vln1bt3b3XLLbcof39/FRERoZo1a6amTJni0noWFYOaiJzhhCelzP79+3Hs2DG0bt3arnd1cHAwWrdujUOHDuGPP/6wLl+/fr31sfy6du0KANiwYUOh7/3VV1/hrrvuQlpaGrp27YoxY8age/fuuHbtGj788EMAuU08Jk6cCACIiorCxIkTrX/y9qJ+5pln0LNnT/z222/o1asXhg8fjqCgIIwbNw4PPPCAw/dfunQp+vbtizp16uCJJ55Ao0aNMH/+fLRp0wbnzp2zPi89PR133nknVq1ahTZt2mD06NHo06cPypUrh9mzZxe6nkREuvL0NwXSV2F71F9++aUCoEaOHOnw8ZEjRyoAau3atUqp3L1sAKphw4YOn79kyRIFQE2YMKHQsfXq1UsBUOnp6XaPnT592ub/UcCh72+//VYBUF27drU5KmA2m9Wjjz6qAKglS5ZYl1v2qAGob775xqbW+PHj7T6P0aNHKwBq+fLlhY5TL9yjJiJnuEddyly4cAEAUKFCBYePh4aG2jyvqM93RVBQkN2y8PBwl1//9ttvAwBmz55tc1RA0zRMnToVmqbh448/tntdp06drEcALJ577jmEhYXhgw8+gNlsewtGScdZkG0HbmJxyjXrnyUbrwMAzGbYLF+ccg2XrrLBHVFpxn7UZJgHHngAy5YtQ6tWrTBgwADcddddaNu2LSIiIopUJzU1FcHBwZg3b57Dx4OCgvDrr7/aLW/btq3dspCQEMTHx2P9+vU4dOgQ6tSpg379+mHGjBn45z//ifvvvx+dO3dGu3btUKNGjSKN0xmlFF75NAtXsgGTgxbZH6Vct/53jhm4mQMkdQrQ5b2JyPswqEsZy56xsz3gixcv2jyvqM8vSN++fbF8+XK8+eabeP/99/HOO+9A0zR06NAB06dPtzkHXZCzZ8/i5s2bmDx5stPnXLlyxW5Z1apVHT7Xstyyji1btsT69evxyiuv4KOPPsL8+fMBAM2bN8drr72GDh06uDROZzRNQ71/+OHnAznIN48CANgtq/cPvxK9HxF5Nx76LmViY2MBwOm9z5bllucFBwcjMjIShw8fRk5OTqHPL0yPHj2wYcMGnDt3DqtWrcKQIUOwfv16JCYm4vz58y7VCA0NRXh4OFTuXQsO/xw+fNjudSdOnHBYz7I875eNtm3bYtWqVTh37hxSUlIwevRo7Nq1C/fee68uk7zcEVcGrsyyX9YPiI9hUBOVZgzqUiY2NhbVq1fH5s2b7fY6r1y5gs2bN6NWrVq45ZZbrMsTEhKsj+W3evVqAEC7du2KNI7y5csjMTERs2fPRnJyMk6cOIG0tDTr4yaTyeEXAyB3j/fMmTNFnmhl48aNdssuX76M9PR0hIaGonbt2naPBwUFoX379pg+fTqeffZZZGVlYc2aNUV6X0da1iv8YJZJAxrH+CGgrIPj40RUajCoSxlN0zBkyBBcvnwZL730ks1jL730Ei5fvoyhQ4faLB82bBgAYMKECbh+/e/zp6tWrcL69evRpUsXREVFFfre33//vcPwPXnyJAAgMDDQuqxSpUr4888/HdYZNWoUAGDw4ME4c+aM3eOZmZnYu3ev3fLvvvvO+sXC4uWXX8b58+cxcOBAmEy5/xy2bNnicHY2y5533nEWV0SoCdFVC/7nZ1ZAqzienSIq7Qxtc0nuMXfuXGzatAkAsGvXLmzbtg2tW7dGnTp1AABt2rTBkCFDrM+/cuUKWrdujR07dqBLly5o0qQJtm3bhm+//RbNmzfHhg0b7K54Hjp0KObOnYsGDRrg3nvvxfHjx/Hpp58iJCQEW7ZsQd26dQsdZ3x8PI4dO4Y2bdogOjoamqZh06ZN+PHHH9GqVSts2rQJfn65h3nvv/9+fPbZZ+jRowcaN24MPz8/dO/eHbfdlttW7oUXXsBLL72EsLAwJCYmIioqCmfOnMGBAwewceNGTJkyBePHjwcALFiwAIMGDUK3bt2wevVq9OnTB9HR0UhNTUVKSgpiYmLw008/oWLFigCAnj17IiUlBe3atUOtWrUQGBiIbdu2Ye3atahduza2b99uvdq9JD5cew0fr78OcwH/Ahc/FYzwUH6fdmTPkZv4eP11HD5hRq2qJvRv748GUfxiQ77H8KBOT9uJOdMWYv/eg4iNi8HQsUlF7ukppYYUycnJWLjQecu1pKQkuyYdFy5cwKRJk7B06VJkZmYiMjISffv2xcSJE1G+vH0/ZLPZjLfffhuzZ8/GgQMHEBISgk6dOuHll19GTIxrk9N/+umnWLZsGX7++WccP34cZcuWRXR0NAYMGIDhw4cjJCTE+tzMzEw8/vjjWLduHc6ePQuz2Yz58+cjOTnZ+pzvvvsOb731FlJTU3H+/HmEh4ejVq1auOeee5CUlGQ9fG8J6vnz56NixYp4+eWXsXv3bpQrVw733XcfXn31VVSrVs1ad/Xq1fjoo4+QlpaGY8eOQSmFmjVromfPnnjyySeLfJW6M/uO5mDUe1edPl67mgnvjgx2+rgnlTQk9Xj9uP/LAlTukQeTBkAD3vhXEMPah0jZ1ns6LwwNailNvCU0Aifj5A3qvEHvaWazQv+pV3DBwX3SJg0Y0MEfD3WUd1tWSUNSj5B9fuFVbDuQY3M0wqQBTer4YUpSuSKtC/fKZZKyrZeQF4YeU5PSxFtCI3Aik0nDHfX94OfgX6FZuXbBWWH2HLmJ5xdexYOvX8bzC69iz5GbJa758frr1pAF/vpb/bXcgNcDwOETZrtTBmaVu9xVli8M2w7k4MxFhW0HcjDu/7JK/Bm54zMvjaRs6yXkhaFBLaWJt4RG4EQA0PLWsg7vpa4QrKFOZMn+eboriEoaknqEbK2qJrvJYkxa7nJX6fGFIT93fealkZRtvYS8MDSopTTxltAInAjIvf2qTL7bpP1MwJ1xftC0kt2W5Y4gAkoeknqEbP/2/oD298xulsPnAzq4fqpAjy8M+bnrMy+NpGzrJeSFoUEtpYm3iEbgZJjk5GQopUSdn7YI9NcQX9sPeTM5xwy0qle2xLXdEURAyUNSj5BtEFUGb/wrCE3q+CE8VEOTOn6YNqQc6td0fXIYPb4w5Oeuz7w0krKtl5AXnr/qe1wy4ls08soaRHr46sfrmLnimvX/y/oBS54PKfFEJ3pdcOVI/ouwBnQIKFJIlvT1enB2UVtRAz8vd37mpZGUbb2n84L3URN52OmLZjz0eu4scSYNaBbrhxcHlnyj7o4g8jV6f2HgZ07uwKAmEuDRmVeQ8dfh0VE9AnBPc39d6krYcy1t+JmT3hjURAJ8uPYaFv/V3pKzkRFRXtwaEAlguWe6djUTQ5qIbHAKHiIB6kSakNisLBqzpSUR5cND30RERILxGBsREZFgPPRNRIZgAwyi4uGhbyJyO3e3peSXAPJlDGoicjt3z5Lmri8B/AJAEhh+jjo9bSdG9B2DLg17YkTfMUhP2+m1NYh8ld6tGt05B7a7GmGwE5bnSdnWezovDN2jltLEW0IjcCKp3LGH6s496gdfv4wzF+03Y+GhGhY/FVLsupy327OkbOsl5IWhe9RSmnhLaAROJJU79lD16JjljDu6YAHshOVpUrb1EvLC0JMtUpp4S2gETqQnPc+luiOgLG0p3TEHdv/2/th2MAsm2B4BKOmXgFpVTTh3yX6PuiRfAHjO23VStvUS8sLQ35DYuBicOXnWZqWL08RbQg0iKfIfqj53KQfbDmYV+1C1OwIKyA3rKUn6b3Lc9SVA7y8Aev+cfJ2Ubb2EvDD00LeUJt4SGoET6UXvQ9XuPEztLrlfAsph8VMhmJKkT0tJyxeAJnX8EB6qoUkdvxK1q3TXRW++Ssq2XkJeGH57lpQm3p5uBE6kF3dcTMVWjfpz10VvvkzKtt7TecH7qIm8HK9O9g78OVFxca5vIi/njYeqSyP+nKi4uEdN5AN4qNo78OdExcGgJiK3M+q2JN7+RL6IQU1ENvQOO3c35HDn+zD4SQKeoyYiK3fMb23UbUl6vw/n+iYpGNREZOWOUDVqKk6934f3PZMUDGoisnJHqLprLm53vw/n+iYpGNREZOWOUDXqtiS938eoLxhEheHFZERewogLm5xdkFWSqTMtdY24LUnP93HXZ1HYmHnBGuXn+SlExyYVuaenlBpERjHqymnLe/Fe31zu/iyM/Ll6Iynbek/nhaFBLaWJt4RG4ERFweknfRN/rs5J2dZLyAtDT7ZIaeItoRE4UVHwwibfxJ+rc1K29RLywtCgltLEW0IjcKKi4IVNvok/V+ekbOsl5IWhvw2xcTHWnp7WARSjibeEGkRGYkMH38Sfq3NStvUS8sLQoJbSxFtCI3CiomgQVQZv/CsITer4ITxUQ5M6frpffUzG48/VOSnbegl54fmrvktpI3AiIiqYlG29p/OC91ETkVvxPmGikmFQE5HbGH2fML8UkC/ipYVE5DZGNrZgtyvyVQxqInIbI+8TZrcr8lUMaiJyGyPvE+bkIeSrePKGqJQx8jxu//b+2HYwCybYnqN2x33CtaqacO6S/XSc7pw8xOhz4jwHXzrxYjIiITzZHcudTSCM7JxlRLerwt7PnRfKsYFH6cSgJhLAqI2wrzeBMLLzl9Gfpa//7Mg5fg0jEsDRhVCmv5ZPSdLvn6mvn8dtEFVG18+rIEZ/lr7+syPnGNREAhi1EXbneVyJ50/dOSajz4l74hw8yWD4Tzg9bSdG9B2DLg17YkTfMUhP2+m1NYj0YtTV0e5qAiHxHmZ3j8nohhqlsYGHlG29p/PC0HPUUpp4S2gETpSXkRdC5d/LbBXnh9S9OSXa65R4/tQdY7L77Or5IfXXHEPOiTt6f3e/nydJ2dZLyAtD96ilNPGW0AicKC8juyjlnscth8VPhaB/e3+8++X1Eu91Sjx/qveYHO2hv/vVdfRv74/FT4VgSpL7u17l/dkZ8X6eJGVbLyEvDD2BJKWJt4RG4ET5GXkhlIVeF7FJPH+q95iMuuCPcknZ1kvIC0P/FUlp4i2hETiRBHrtdUo8f6r3mCQeNfBlUrb1EvLC0KCW0sRbQiNwIgn0uojNyEP3nhqTkdOhkpxtvYS8MHzCEylNvD3dCJxIAqNn8/JmUj4ribfBuYuUbb2n84IzkxF5mKc3vKXpSuKS8vRnxWlESycGNZEHccNLRSHxNjhyP55cIfIg9lCmouAFbaUTv7ITeRA3vO7j6VMK7iDxNjhyP+/+rSXycr644ZUQkPlPKZy7lINtB7O8/pSCkf29SQ6eoybyIClXEutFyjl3Sedy9f7i4ukL2sh4DGoiD9N7w+vJPVopAfng65dx5qL9pi08VMPip0IMG4eULy7k3fibQlRCJQ1GPacO9fQhXynn3KWcUpA47aiEUxNUNN57IoxIAGntHT19FbmU2bukTGkq5YuLhbTfV3INg5qoBDwdjPl5OhikBKSUKU2lfHGxkPb7Sq7h8Q6iEvB0MObn6UO+loCUcLGTJ7qR5SftKm1pv6/kGsN/i+3mTB2bVOTm21JqEHk6GPOTEAwSAlIKSV9cAHm/r4WRsq33dF4YetV3etpODLnvMWsTbksXkrkrZ7q80lJqEAEyb6/i7TvkjMTfV2ekbOsl5IWhX6PmTFtoXVkgt/m2gsKcaQu9rgYRIOdcaP4xTUkqh8VPhWBKkrwNMHmOxN9XZ6Rs6yXkhaHHp/bvPWhdWQtzjhn79x70uhpEFjzUS97EW35fpWzrJeSFoXvUsXEx1ubb1gH4mRAbF+N1NYiIyH2kbOsl5IWhQT10bBI0aNaVthzrHzou2etqEBGR+0jZ1kvIC8OnELW7em5cMuJbNPLKGkRE5D5StvWezgvO9U1ERCSYzJvniIiICACDmoiISDQGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCSY4UGdnrYTI/qOQZeGPTGi7xikp+302hqSxuJLNSSNxZdqSBqLL9WQNBZfqiFpLHqtT3EZOoWolCbeejUClzIWX6ohaSy+VEPSWHyphqSx+FINSWPRa31KwtA9ailNvPVqBC5lLL5UQ9JYfKmGpLH4Ug1JY/GlGpLGotf6lIShQS2libdejcCljMWXakgaiy/VkDQWX6ohaSy+VEPSWPRan5IwNKilNPHWqxG4lLH4Ug1JY/GlGpLG4ks1JI3Fl2pIGote61MShga1lCbeejUClzIWX6ohaSy+VEPSWHyphqSx+FINSWPRa31KwvB+1FKaeOvVCFzKWHyphqSx+FINSWPxpRqSxuJLNSSNRa/1KS7Dg5qIiIhcxwlPiIiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiAQzPKilNPFmY3O5NSSNxZdqSBqLL9WQNBZfqiFpLHqtT3EZOoWolCbebGwut4aksfhSDUlj8aUaksbiSzUkjUWv9SkJQ/eopTTxZmNzuTUkjcWXakgaiy/VkDQWX6ohaSx6rU9JGBrUUpp4s7G53BqSxuJLNSSNxZdqSBqLL9WQNBa91qckDA1qKU282dhcbg1JY/GlGpLG4ks1JI3Fl2pIGote61MShga1lCbebGwut4aksfhSDUlj8aUaksbiSzUkjUWv9SkJw/tRS2nizcbmcmtIGosv1ZA0Fl+qIWksvlRD0lj0Wp/iMjyoiYiIyHWc8ISIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMMODWkoTbzY2l1tD0lh8qYaksfhSDUlj8aUaksai1/oUl6FTiEpp4s3G5nJrSBqLL9WQNBZfqiFpLL5UQ9JY9FqfkjB0j1pKE282NpdbQ9JYfKmGpLH4Ug1JY/GlGpLGotf6lIShQS2liTcbm8utIWksvlRD0lh8qYaksfhSDUlj0Wt9SsLQoJbSxJuNzeXWkDQWX6ohaSy+VEPSWHyphqSx6LU+JWFoUEtp4s3G5nJrSBqLL9WQNBZfqiFpLL5UQ9JY9FqfkjC8H7WUJt5sbC63hqSx+FINSWPxpRqSxuJLNSSNRa/1KS7Dg5qIiIhcxwlPiIiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREgjGoiYiIBGNQExERCcagJiIiEoxBTUREJBiDmoiISDAGNRERkWAMaiIiIsEY1ERERIIxqImIiARjUBMREQnGoCYiIhKMQU1ERCQYg5qIiEgwBjUREZFgDGoiIiLBGNRERESCMaiJiIgEY1ATEREJxqAmIiISjEFNREQkGIOaiIhIMAY1ERGRYAxqIiIiwRjUREREgjGoiYiIBGNQExERCcagJiIiEuz/AVlprPuxIPQ4AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "def ploting_B(\n",
    "    optim_final_state, optim_init_state, target_state_final, target_state_init, mask\n",
    "):\n",
    "    # Add two plots side by side\n",
    "    fig = plt.figure(figsize=(5, 5))\n",
    "    turbo_map = cm.get_cmap(\"turbo\")\n",
    "    global_size = 15\n",
    "    plt.scatter(\n",
    "        target_state_init[\"r\"][mask][:, 0],\n",
    "        target_state_init[\"r\"][mask][:, 1],\n",
    "        label=\"target 0\",\n",
    "        s=global_size,\n",
    "        c=turbo_map(0.15),\n",
    "        alpha=0.4,\n",
    "    )\n",
    "    plt.scatter(\n",
    "        target_state_final[\"r\"][mask][:, 0],\n",
    "        target_state_final[\"r\"][mask][:, 1],\n",
    "        label=\"target 100\",\n",
    "        s=global_size,\n",
    "        c=turbo_map(0.15),\n",
    "    )\n",
    "    plt.scatter(\n",
    "        optim_init_state[\"r\"][mask][:, 0],\n",
    "        optim_init_state[\"r\"][mask][:, 1],\n",
    "        label=\"init\",\n",
    "        s=global_size,\n",
    "        c=turbo_map(0.7),\n",
    "    )\n",
    "    plt.scatter(\n",
    "        optim_final_state[\"r\"][mask][:, 0],\n",
    "        optim_final_state[\"r\"][mask][:, 1],\n",
    "        label=\"updated\",\n",
    "        s=global_size,\n",
    "        c=turbo_map(0.85),\n",
    "    )\n",
    "    wall = optim_init_state[\"r\"][~mask]\n",
    "    top_lid = wall[:, 1] < 1.17\n",
    "    plt.scatter(\n",
    "        wall[top_lid][:, 0],\n",
    "        wall[top_lid][:, 1],\n",
    "        label=\"wall\",\n",
    "        s=global_size,\n",
    "        color=turbo_map(0),\n",
    "    )\n",
    "    plt.text(0.72, 0.91, \"optimizer\\n15 steps\", fontsize=14, ha=\"center\", rotation=0)\n",
    "    plt.text(0.22, 0.94, \"Loss: 0.32\", fontsize=14)\n",
    "    plt.arrow(\n",
    "        0.6,\n",
    "        0.5,\n",
    "        0.0,\n",
    "        -0.10,\n",
    "        head_width=0.03,\n",
    "        head_length=0.05,\n",
    "        fc=turbo_map(0.15),\n",
    "        ec=turbo_map(0.15),\n",
    "    )\n",
    "    plt.text(0.40, 0.38, \"SPH solver\\n100 steps\", fontsize=14, ha=\"center\", rotation=0)\n",
    "    plt.text(0.68, 0.78, \"Loss: 0.02\", fontsize=14)\n",
    "    plt.axis(\"equal\")\n",
    "    plt.axis(\"off\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    fig.savefig(\"optim_B.pdf\", bbox_inches=\"tight\")\n",
    "    print(\"done\")\n",
    "\n",
    "\n",
    "mask = target_state_B[\"tag\"] == Tag.FLUID\n",
    "\n",
    "ploting_B(final_optim_B, init_optim_B, final_state_B, target_state_B, mask)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}