MODEL02_LDM_3d_main.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7184eacc-7e1f-4419-b3c2-028a9009eff5",
   "metadata": {},
   "source": [
    "# Latent diffusion model for generating ensembl members\n",
    "\n",
    "This notebook contains information for the design, hyperparameters, and training of a LDM that generates ensemble members from the bias-corrected ensemble mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "969ca223-fe95-44f5-9fed-4f91df5814d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import time\n",
    "import math\n",
    "import h5py\n",
    "import logging\n",
    "import warnings\n",
    "import numpy as np\n",
    "from glob import glob\n",
    "\n",
    "# supress regular warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' \n",
    "logging.getLogger(\"tensorflow\").setLevel(logging.ERROR)\n",
    "\n",
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers\n",
    "\n",
    "# supress tensorflow warnings\n",
    "tf.autograph.set_verbosity(0)\n",
    "tf.get_logger().setLevel('ERROR')\n",
    "\n",
    "# adjust for time step embedding layer\n",
    "tf.config.run_functions_eagerly(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "22cfec79-8be6-4215-8c3a-0dd98fcc26e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.insert(0, '/glade/u/home/ksha/GAN_proj/')\n",
    "sys.path.insert(0, '/glade/u/home/ksha/GAN_proj/libs/')\n",
    "\n",
    "from namelist import *\n",
    "import data_utils as du\n",
    "import model_utils as mu\n",
    "import verif_utils as vu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "145d2a0e-a885-4356-b54a-ed3e606c1eab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "389c7e8d-d7be-402d-9bc1-6efe81992eaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu.set_seeds(888)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "482358b2-f48d-415c-a7c9-b18b882b9d35",
   "metadata": {},
   "source": [
    "## Import VQ-VAE decoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8ea6b59e-3b7b-4f22-aeb7-789e315e02ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hyperparameters\n",
    "filter_nums = [64, 128] # number of convolution kernels per down-/upsampling layer \n",
    "latent_dim = 4 # number of latent feature channels\n",
    "activation = 'gelu' # activation function\n",
    "num_embeddings = 128 #128 # number of the VQ codes\n",
    "\n",
    "input_size = (224, 464, 1) # size of MRMS input\n",
    "latent_size = (14, 29, latent_dim) # size of compressed latent features\n",
    "\n",
    "drop_encode = False\n",
    "drop_decode = False\n",
    "\n",
    "model_name_encoder_load = model_dir+'models/VQ_VAE_encoder_stack1_tune0'\n",
    "model_name_decoder_load = model_dir+'models/VQ_VAE_decoder_stack1_tune0'\n",
    "\n",
    "encoder = mu.VQ_VAE_encoder(input_size, filter_nums, latent_dim, num_embeddings, activation, drop_encode)\n",
    "W_old = mu.dummy_loader(model_name_encoder_load)\n",
    "encoder.set_weights(W_old)\n",
    "\n",
    "decoder = mu.VQ_VAE_decoder(latent_size, filter_nums, activation, drop_decode)\n",
    "W_old = mu.dummy_loader(model_name_decoder_load)\n",
    "decoder.set_weights(W_old)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4513740d-4fd0-4bca-93a9-4e0e47da2d35",
   "metadata": {},
   "source": [
    "## LDM config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "12aef6a8-38a2-411f-9459-e4c6e669a7a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# diffusion steps\n",
    "diffusion_steps = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c7dfb9d5-d852-47a8-9f85-eb2e01941056",
   "metadata": {},
   "outputs": [],
   "source": [
    "# min-max values\n",
    "ccpa_min = -2.0606\n",
    "ccpa_max = 1.6031\n",
    "\n",
    "LEADs = np.arange(6, 144+6, 6) # forecast lead times\n",
    "N_leads = len(LEADs)\n",
    "\n",
    "# The tensor size of embedded CCPA and GEFS ensemble mean \n",
    "latent_size = (14, 29, 4)\n",
    "input_size = (N_leads,) + latent_size\n",
    "\n",
    "# model design\n",
    "widths = [32, 64, 96, 128]\n",
    "embedding_dims = 32\n",
    "block_depth = 2\n",
    "\n",
    "min_signal_rate = 0.02\n",
    "max_signal_rate = 0.95\n",
    "ema = 0.999\n",
    "\n",
    "load_weights = True\n",
    "# location of the previous weights\n",
    "model_name_load = model_dir+'baseline/LDM_step{:02d}_base/'.format(diffusion_steps)\n",
    "\n",
    "# location for saving new weights\n",
    "model_name_save = model_dir+'baseline/LDM_step{:02d}_base/'.format(diffusion_steps)\n",
    "\n",
    "epochs = 99999\n",
    "N_batch = 2\n",
    "batch_size_ = 64\n",
    "batch_size = batch_size_*32\n",
    "lr = 1e-4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "16e0f746-5211-4581-a825-22f4decb00e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rescale(x, min_val=ccpa_min, max_val=ccpa_max):\n",
    "    return ((x - min_val) / (max_val - min_val))* 2 - 1\n",
    "\n",
    "def scaleback(x, min_val=ccpa_min, max_val=ccpa_max):\n",
    "    return 0.5*(x+1)*(max_val-min_val) + min_val\n",
    "\n",
    "def to_precip(x):\n",
    "    x[x<0] = 0\n",
    "    return 10*(np.exp(x) - 1)\n",
    "\n",
    "def valid_func(valid_true_CCPA, decoder, LDM, valid_ViT, land_mask_CCPA):\n",
    "    L_valid, N_leads, EN, _, _, _ = valid_ViT.shape\n",
    "    \n",
    "    Y_pred_ViT = np.empty((L_valid, N_leads, EN, 224, 464))\n",
    "    \n",
    "    for ien in range(EN):\n",
    "        y_pred_vit = LDM.generate(L_valid, valid_ViT[:, :, ien, ...])\n",
    "        y_pred_vit = scaleback(np.array(y_pred_vit))\n",
    "        \n",
    "        for ilead in range(N_leads):\n",
    "            Y_pred_ViT[:, ilead, ien, ...] = decoder.predict(y_pred_vit[:, ilead, ...], verbose=0)[..., 0]\n",
    "\n",
    "    Y_pred_ViT = to_precip(Y_pred_ViT)\n",
    "    Y_pred_ViT[Y_pred_ViT<0] = 0\n",
    "    \n",
    "    CRPS = np.empty((L_valid, N_leads, 224, 464)); CRPS[...] = np.nan\n",
    "    \n",
    "    for ilead in range(N_leads):\n",
    "    \n",
    "        crps_ilead, _, _ = vu.CRPS_2d(valid_true_CCPA[:, ilead, ...], \n",
    "                                      Y_pred_ViT[:, ilead, ...], land_mask=land_mask_CCPA)\n",
    "        CRPS[:, ilead, ...] = crps_ilead\n",
    "\n",
    "    return np.nanmean(CRPS)\n",
    "    \n",
    "def base_func(valid_true_CCPA, decoder, valid_ViT, land_mask_CCPA):\n",
    "    L_valid, N_leads, EN, _, _, _ = valid_ViT.shape\n",
    "    \n",
    "    Y_pred_ViT = np.empty((L_valid, N_leads, EN, 224, 464))\n",
    "    \n",
    "    for ien in range(EN):        \n",
    "        for ilead in range(N_leads):\n",
    "            Y_pred_ViT[:, ilead, ien, ...] = decoder.predict(valid_ViT[:, ilead, ien, ...], verbose=0)[..., 0]\n",
    "\n",
    "    Y_pred_ViT = to_precip(Y_pred_ViT)\n",
    "    Y_pred_ViT[Y_pred_ViT<0] = 0\n",
    "    \n",
    "    CRPS = np.empty((L_valid, N_leads, 224, 464)); CRPS[...] = np.nan\n",
    "    \n",
    "    for ilead in range(N_leads):\n",
    "    \n",
    "        crps_ilead, _, _ = vu.CRPS_2d(valid_true_CCPA[:, ilead, ...], \n",
    "                                      Y_pred_ViT[:, ilead, ...], land_mask=land_mask_CCPA)\n",
    "        CRPS[:, ilead, ...] = crps_ilead\n",
    "\n",
    "    return np.nanmean(CRPS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8419e6eb-765a-4c31-87b0-27c22eeb6390",
   "metadata": {},
   "source": [
    "## Prepare validation set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "26ec52f4-4be0-496f-95e3-3ba6b83efe06",
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(save_dir+'CCPA_domain.hdf', 'r') as h5io:\n",
    "    land_mask_CCPA = h5io['land_mask_CCPA'][...]\n",
    "    \n",
    "land_mask_CCPA = land_mask_CCPA == 1.0\n",
    "ccpa_shape = land_mask_CCPA.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e3836a8a-f073-417b-8e04-cd56928474c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "BATCH_dir = camp_dir+'BATCH_LDM_member_OPT/'\n",
    "filenames = sorted(glob(BATCH_dir+'*npy'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8f8550fd-c2fc-4206-8310-1e58c776c782",
   "metadata": {},
   "outputs": [],
   "source": [
    "EN = 31\n",
    "\n",
    "BATCH_dir = camp_dir+'BATCH_LDM_member_OPT/'\n",
    "filenames = sorted(glob(BATCH_dir+'*npy'))\n",
    "\n",
    "L_valid = 20\n",
    "filenames_valid = filenames[::7][:L_valid]\n",
    "\n",
    "valid_true_CCPA = np.empty((L_valid, N_leads,)+ccpa_shape)\n",
    "valid_ViT = np.empty((L_valid, N_leads, EN)+latent_size)\n",
    "\n",
    "for i, name_ in enumerate(filenames_valid):\n",
    "    temp_data = np.load(name_, allow_pickle=True)[()]\n",
    "    valid_true_CCPA[i, ...] = temp_data['CCPA_true'][:N_leads, ...]\n",
    "    valid_ViT[i, ...] = temp_data['ViT_embed'][:N_leads, ...]\n",
    "    \n",
    "valid_ViT = rescale(valid_ViT)\n",
    "valid_ViT = valid_ViT[:, :, :5, ...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f215ec22-ba1f-47dd-839c-eecea00acfea",
   "metadata": {},
   "outputs": [],
   "source": [
    "BATCH_dir = camp_dir+'BATCH_LDM_member/'\n",
    "filename_train = sorted(glob(BATCH_dir+'*npy'))\n",
    "\n",
    "# filename_train = list(set(filenames) - set(filenames_valid))\n",
    "L_train = len(filename_train)\n",
    "\n",
    "min_del = 0.0\n",
    "max_tol = 3 # early stopping with 2-epoch patience\n",
    "tol = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8705589c-496f-42bf-8175-0493aee1277c",
   "metadata": {},
   "outputs": [],
   "source": [
    "decay_steps = batch_size/batch_size_\n",
    "cosine_decay = tf.keras.optimizers.schedules.CosineDecay(initial_learning_rate=lr, \n",
    "                                                         decay_steps=decay_steps, \n",
    "                                                         alpha=0.1)\n",
    "\n",
    "model= mu.DiffusionModel(input_size, input_size, input_size, \n",
    "                        diffusion_steps, min_signal_rate, max_signal_rate, \n",
    "                        embedding_dims, widths, block_depth, ema)\n",
    "\n",
    "model.compile(loss=keras.losses.mean_absolute_error, optimizer=keras.optimizers.AdamW(learning_rate=cosine_decay))\n",
    "\n",
    "# load previous weights\n",
    "if load_weights:\n",
    "    model.load_weights(model_name_load)\n",
    "\n",
    "batch_CCPA = np.empty((batch_size, N_leads,)+latent_size)\n",
    "batch_CCPA[...] = np.nan\n",
    "\n",
    "batch_ViT = np.empty((batch_size, N_leads,)+latent_size)\n",
    "batch_ViT[...] = np.nan\n",
    "\n",
    "for i in range(epochs):\n",
    "    print('epoch = {}'.format(i))\n",
    "    if i == 0:\n",
    "        if load_weights:\n",
    "            record = valid_func(valid_true_CCPA, decoder, model, valid_ViT, land_mask_CCPA)\n",
    "        else:\n",
    "            record = 999\n",
    "        print('Initial validation loss: {}'.format(record))\n",
    "        \n",
    "    start_time = time.time()\n",
    "    for j in range(N_batch):\n",
    "        #start_time = time.time()\n",
    "        inds_rnd = du.shuffle_ind(L_train)\n",
    "        inds_ = inds_rnd[:batch_size]\n",
    "        \n",
    "        for k, ind in enumerate(inds_):\n",
    "            # import batch data\n",
    "            name_ = filename_train[ind]\n",
    "            temp_data = np.load(name_, allow_pickle=True)[()]\n",
    "            batch_CCPA[k, ...] = temp_data['CCPA_embed'][:N_leads, ...]\n",
    "            batch_ViT[k, ...] = temp_data['ViT_embed'][:N_leads, ...]\n",
    "\n",
    "        batch_CCPA = rescale(batch_CCPA)\n",
    "        batch_ViT = rescale(batch_ViT)\n",
    "        \n",
    "        model.fit([batch_CCPA, batch_ViT], batch_size=batch_size_, epochs=1, verbose=0)\n",
    "        #print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "    # on epoch-end\n",
    "    record_temp = valid_func(valid_true_CCPA, decoder, model, valid_ViT, land_mask_CCPA)\n",
    "    \n",
    "    # print out valid loss change\n",
    "    if record - record_temp > min_del:\n",
    "        print('Validation loss improved from {} to {}'.format(record, record_temp))\n",
    "        record = record_temp\n",
    "        print(\"Save to {}\".format(model_name_save))\n",
    "        model.save_weights(model_name_save)\n",
    "    else:\n",
    "        print('Validation loss {} NOT improved'.format(record_temp))\n",
    "        \n",
    "    print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "    # mannual callbacks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4840a2cc-082f-4efc-bb97-7b7f2a1294ca",
   "metadata": {},
   "source": [
    "## Predition tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "5cd75f0b-e4f5-4296-96dd-277db63226d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = mu.DiffusionModel(input_size, input_size, input_size, \n",
    "                        diffusion_steps, min_signal_rate, max_signal_rate, \n",
    "                        embedding_dims, widths, block_depth, ema)\n",
    "# load previous weights\n",
    "model.load_weights(model_name_save);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "1ef12138-a1a0-4018-924c-7528ac6de412",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1048040653674857"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "valid_func(valid_true_CCPA, decoder, model, valid_ViT, land_mask_CCPA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0acdc8a2-aecc-4f84-a1f8-66dc6dcf2a28",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10541903031870066"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "base_func(valid_true_CCPA, decoder, valid_ViT, land_mask_CCPA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "88eade5d-a36c-49be-89ae-01feff4e8c4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def inference_func(valid_true_CCPA, decoder, LDM, valid_ViT, land_mask_CCPA):\n",
    "    L_valid, N_leads, EN, _, _, _ = valid_ViT.shape\n",
    "    \n",
    "    Y_pred_ViT = np.empty((L_valid, N_leads, EN, 224, 464))\n",
    "    \n",
    "    for ien in range(EN):\n",
    "        y_pred_vit = LDM.generate(L_valid, valid_ViT[:, :, ien, ...])\n",
    "        y_pred_vit = scaleback(np.array(y_pred_vit))\n",
    "        \n",
    "        for ilead in range(N_leads):\n",
    "            Y_pred_ViT[:, ilead, ien, ...] = decoder.predict(y_pred_vit[:, ilead, ...], verbose=0)[..., 0]\n",
    "\n",
    "    Y_pred_ViT = to_precip(Y_pred_ViT)\n",
    "    Y_pred_ViT[Y_pred_ViT<0] = 0\n",
    "    return Y_pred_ViT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "410f96b4-41cc-457a-8b54-4801dc692a7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "examples = inference_func(valid_true_CCPA, decoder, model, valid_ViT, land_mask_CCPA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "ec4b480c-46af-4bd7-a925-32366c337a29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x14d2d5badd80>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAGiCAYAAABppIV1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhZ0lEQVR4nO39f3xU9Z33/z/mR2Ymv2OAECNBqFh/FJEWLc1q3ShUQD9urex1XVq7xepH2y64Fbpby3WzRa17w0u7rbW12NvWBXtVll29Cr20Ky7+AEoFKrT5YPqDLXyhJJUgYpMwIb9m5nz/eJ8zc2YyCUnmJEzC8+7t3Cbz+8wJZl7n9X69X2+fZVkWIiIiIv3wn+kdEBERkfymYEFEREQGpGBBREREBqRgQURERAakYEFEREQGpGBBREREBqRgQURERAakYEFEREQGpGBBREREBqRgQURERAakYEFERGSMWrNmDbNmzaKsrIyysjLq6up4+eWXk/fX19fj8/nSti984QtDfh+f1oYQEREZm1588UUCgQAXXnghlmXx7LPP8vjjj/PrX/+aD33oQ9TX1/PBD36Qhx9+OPmcoqIiysrKhvQ+Qa93XEREREbHTTfdlHb9H//xH1mzZg27du3iQx/6EGCCg+rq6pzeZ0wGC4lEgnfeeYfS0lJ8Pt+Z3h0REcljlmVx8uRJampq8PtHbvS9q6uLnp6enF/Hsqw+323hcJhwODzg8+LxOM8//zwdHR3U1dUlb3/uuef48Y9/THV1NTfddBNf+9rXKCoqGvJOjTlNTU0WoE2bNm3atA16a2pqGrHvpc7OTs/2s6SkpM9tq1at6ve99+3bZxUXF1uBQMAqLy+3fvaznyXv+8EPfmBt3rzZ2rdvn/XjH//YOu+886xPfepTQ/58Y7Jmoa2tjYqKCn72s59RXFx8pndHRETyWEdHBzfeeCOtra2Ul5ePyHu0t7dTXl5OhEIgl4y3RRedNDU1pdUVDJRZ6Onp4ciRI7S1tfHCCy/wwx/+kG3btnHppZf2eezrr7/OvHnzOHDgABdccMGg92pMDkM46Zni4mJKSkrO8N6IiMhYMDrD1j58OQQLztm7M7thMEKhEDNmzABgzpw5vPXWW3znO9/hBz/4QZ/Hzp07F+DsCBZERETyUSGF+HPoSpAgQRenctqHRCJBd3d31vsaGhoAOPfcc4f0mgoWREREPFJIcc7Bwp85MejHr1y5kkWLFjF16lROnjzJ+vXr2bp1K6+88goHDx5k/fr13HDDDUyYMIF9+/axfPlyrrnmGmbNmjWk/VKwICIi4pEIRQQIDPv5ceJDevy7777LZz/7WY4ePUp5eTmzZs3ilVde4ROf+ARNTU28+uqrPPHEE3R0dFBbW8vixYt54IEHhrxfChZERETGqGeeeabf+2pra9m2bZsn76NgQURExCNFo5xZGC0KFkRERDwSoYhgDl+tMWIe7o13tJCUiIiIDEiZBREREY8UUkSQgmE/P0avh3vjHQULIiIiHimkiIIcgoXePA0WNAwhIiIiA1JmQURExCMmsxAa9vOD5L5q5UhQsCAiIuKRCMWEcggWAjkMYYwkDUOIiIjIgJRZEBER8UgRhYTIvpT0YARzaOg0khQsiIiIeCRCEWEiw36+X8GCiIjI+BYmklOwkK9UsyAiIiIDUmZBRETEI377v1yen48ULIiIiHjEZ/+Xy/PzUX6GMCIiIpI3lFkQERHxyHjNLChYEBER8ch4rVnIz70SERGRvKHMgoiIiEc0DCEiIiIDGq/BgoYhREREZEDKLIiIiHjEhy+nIsV8zSwoWBAREfHIeB2GULAgIiLikfEaLKhmQURERAakzIKIiIhHxmtTJgULIiIiHtEwhIiIiJyVhhQsrF69miuvvJLS0lKqqqq4+eab2b9/f9pjurq6WLp0KRMmTKCkpITFixdz7NixtMccOXKEG2+8kaKiIqqqqviHf/gHYrFY7p9GRETkDPJ58F8+GlKwsG3bNpYuXcquXbvYsmULvb29XH/99XR0dCQfs3z5cl588UWef/55tm3bxjvvvMMtt9ySvD8ej3PjjTfS09PDm2++ybPPPsu6dev4+te/7t2nEhEROQNi+HPe8tGQahY2b96cdn3dunVUVVWxd+9errnmGtra2njmmWdYv3491113HQBr167lkksuYdeuXXzsYx/jP//zP/ntb3/Lq6++yuTJk5k9ezbf+MY3uP/++3nwwQcJhUJ93re7u5vu7u7k9fb29uF8VhERERmGnEKYtrY2ACorKwHYu3cvvb29zJ8/P/mYiy++mKlTp7Jz504Adu7cyWWXXcbkyZOTj1mwYAHt7e385je/yfo+q1evpry8PLnV1tbmstsiIiIj4hQ+OnLYTo2HYQi3RCLBfffdx1VXXcXMmTMBaGlpIRQKUVFRkfbYyZMn09LSknyMO1Bw7nfuy2blypW0tbUlt6ampuHutoiIyIgZr8HCsKdOLl26lMbGRnbs2OHl/mQVDocJh8Mj/j4iIiLS17CChWXLlvHSSy+xfft2pkyZkry9urqanp4eWltb07ILx44do7q6OvmYX/7yl2mv58yWcB4jIiIyFp3KsUixJ08LHIe0V5ZlsWzZMjZu3Mjrr7/O9OnT0+6fM2cOBQUFvPbaa8nb9u/fz5EjR6irqwOgrq6Ot99+m3fffTf5mC1btlBWVsall16ay2cRERE5ozQMgRl6WL9+PT/96U8pLS1N1hiUl5dTWFhIeXk5d911FytWrKCyspKysjLuvfde6urq+NjHPgbA9ddfz6WXXsrf/M3f8Nhjj9HS0sIDDzzA0qVLNdQgIiJjWgc+CnL4wu8dD8HCmjVrAKivr0+7fe3atdxxxx0AfPvb38bv97N48WK6u7tZsGAB3//+95OPDQQCvPTSS3zxi1+krq6O4uJilixZwsMPP5zbJxEREZERMaRgwbKs0z4mEonw1FNP8dRTT/X7mPPPP5//+I//GMpbi4iI5L1O/PTmUHcwLpoyiYiISP868BHMYSghlqfDEPkZwoiIiEjeULAgIiLikVOQ42yIoVmzZg2zZs2irKyMsrIy6urqePnll5P3D2Zxx8FQsCAiIuKRTvycymHrHOLX8pQpU3j00UfZu3cve/bs4brrruOTn/xkcvmE0y3uOFiqWRAREckzmQsm9tfJ+Kabbkq7/o//+I+sWbOGXbt2MWXKlNMu7jhYZ2+w0Nrad3M43Sfdl5m3iYiIZOjAhz+HIsWE/dzMBRNXrVrFgw8+OOBz4/E4zz//PB0dHdTV1Z12cUcFCyIiImfAKXz4cggWLPu5TU1NlJWVJW8fqGnh22+/TV1dHV1dXZSUlLBx40YuvfRSGhoaTru442CdvcFCRUUqm+DOKjjXnfudg5ztNhERkRHgFCwOxkUXXURDQwNtbW288MILLFmyhG3btnm6P2dvsCAiIuKxU/jx5TB3wBrGc0OhEDNmzADMGk1vvfUW3/nOd/gf/+N/nHZxx8E6u2dDZGYIMjMMIiIiQ5AgQDyHLUEg931IJOju7h7U4o6DdXZnFpzCRfeQxGCeIyIiklWQ3L5ah/bclStXsmjRIqZOncrJkydZv349W7du5ZVXXhnU4o4js1ciIiKSN959910++9nPcvToUcrLy5k1axavvPIKn/jEJ4DTL+44WAoWhppZEBER6dfoZhaeeeaZAe8fzOKOg6FgQURExDMBcvtqzb1mYSSM/QLHw4cpuuKK4T/fXbcgIiIifYz5zMKkv/5rJgEtV1zBKaBr5kx44omhffm7gwX1URARkWEb3WGI0ZKfeyUiIjImjc9gYcwPQ1QAhZ8wl6VAUWMjzJ8P69bB4cNDfLGK/rMKyjaIiMhZKj9DGBERkTFpfGYW8nOvhuAPP/whc575fzlnOhQfghYgAXRt2mTqD+rrUw8eqJhxMJkDZRdERGRAmg2Rn2bPZu+XX4ASCBWaoYgIEG5uhh//GDZtMltDg9my9VPQ0IOIiEi/xnxmQUREJH9oGCJ/TZsGJUA1lB0ySZxWoBPofOkl4gAzZ8Ls2eax2bg7OWZmFJRhEBGRQVGwICIiIgMan8HC2K9ZcJilvAlMh7LpcB5QBRQCIYDGRjOV8nTTKfvLKii7ICIiZ6lxEywc/d9AlORwRGC66b0QwnzIQX1Qd0Cg4EBERIYs6MGWf/Jzr0RERMYkTZ3Ma+9885twNabRwgFz2YvpuWABiZKS0y9HrUWlRERE+lBmQURExCsRwJfD8y2gy6N98dD4CRbq6znw93AB5jj32lsMzNTJadOGnjFQcaOIiAxFIbnl7BPkZbAw5I+0fft2brrpJmpqavD5fGzatCntfp/Pl3V7/PHHk4+ZNm1an/sfffTRnD9M2549/Bemx8IJ4M+kAoY+wwuDHYrIx0ChtbX//RcREfHYkDMLHR0dXH755dx5553ccsstfe4/evRo2vWXX36Zu+66i8WLF6fd/vDDD3P33Xcnr5eWlg51V0RERPKLF5mFP3u0Lx4acrCwaNEiFi1a1O/91dXVadd/+tOfcu211/KBD3wg7fbS0tI+j/VCdM8eYldcQQAz9NMLMGVK9sJFp1ujc5aeWfyYj1kFN/d+i4jImRchtwkNca92xFsjOhvi2LFj/OxnP+Ouu+7qc9+jjz7KhAkT+PCHP8zjjz9OLBbr93W6u7tpb29P20RERGR0jGiB47PPPktpaWmf4Yq/+7u/4yMf+QiVlZW8+eabrFy5kqNHj/Ktb30r6+usXr2ahx56aNDv2/XAA4QeeYQE9pRJd3Fj5voPmWP/+b4uRLZahWzrWYiIyOgrYlxmFkY0WPiXf/kXbr/9diKRSNrtK1asSP48a9YsQqEQn//851m9ejXhcLjP66xcuTLtOe3t7dTW1vb/xvX19DzyiEmbTJuWChaGOxtiLFDAICJy5kXI7Zu1/yT7GTViwcLPf/5z9u/fz7/927+d9rFz584lFotx+PBhLrrooj73h8PhrEGEiIhIXilEwcJQPPPMM8yZM4fLL7/8tI9taGjA7/dTVVXlzZtXVMDMmSRaW7NnFU53Bq4zdBERkaQhBwvRaJQDBw4krx86dIiGhgYqKyuZOnUqYIYJnn/+ef7pn/6pz/N37tzJ7t27ufbaayktLWXnzp0sX76cz3zmM5xzzjk5fJQMs2ebFSbdfRMGGopQgCAiIrkqBApyeH6vVzvirSEHC3v27OHaa69NXndqCZYsWcK6desA2LBhA5Zlcdttt/V5fjgcZsOGDTz44IN0d3czffp0li9fnlaTICIiMiblGizkaV/lIe9WfX09lmUN+Jh77rmHe+65J+t9H/nIR9i1a9dQ33boKir6Fjfme2fGwRhoISwREZERkKcxjIiIyBgUAUI5PD8/V6gex8GCO4MwXrIKIiKS34rILVjI02/lPN0tDwwULIiIiMigjd9gQUREZLRFgFzaAo3oIgzDN36DhWyZhfEiW5HjePp8IiJjVa7Bgs+rHfFWnsYwIiIiki/Gb2bBbTyedY/HzyQiMtb5ye00PE9P4cdvsKCZDyIiMtoC5Db9UVMnRURExjk/uX3hK7MwypRZEBER8cT4DRZERERGm2oWxhhlFEREZLSN05qFPI1hPKKAQUREJGfjN7MgIiIy2pRZEBERkQH5PdiGYPXq1Vx55ZWUlpZSVVXFzTffzP79+9MeU19fj8/nS9u+8IUvDPljiYiIyBi0bds2li5dyq5du9iyZQu9vb1cf/31dHR0pD3u7rvv5ujRo8ntscceG9L7aBhCRETEKx4NQ7S3t6fdHA6HCYf7LjqxefPmtOvr1q2jqqqKvXv3cs011yRvLyoqorq6eti7pcyCiIiIVwIebEBtbS3l5eXJbfXq1YN6+7a2NgAqKyvTbn/uueeYOHEiM2fOZOXKlZw6dWpIH0uZBRERkTzT1NREWVlZ8nq2rEKmRCLBfffdx1VXXcXMmTOTt3/605/m/PPPp6amhn379nH//fezf/9+fvKTnwx6fxQsiIiIeMVHbjl7e4nqsrKytGBhMJYuXUpjYyM7duxIu/2ee+5J/nzZZZdx7rnnMm/ePA4ePMgFF1wwqNfWMISIiIhXPBqGGKply5bx0ksv8cYbbzBlypQBHzt37lwADhw4MOjXV2ZBRETEK6PcZ8GyLO699142btzI1q1bmT59+mmf09DQAMC555476PdRsCAiIjJGLV26lPXr1/PTn/6U0tJSWlpaACgvL6ewsJCDBw+yfv16brjhBiZMmMC+fftYvnw511xzDbNmzRr0+yhYEBER8Yo/YbZcnj8Ea9asAUzjJbe1a9dyxx13EAqFePXVV3niiSfo6OigtraWxYsX88ADDwzpfRQsiIiIeMUXB18st+cPgWVZA95fW1vLtm3bhr8/NhU4ioiIyICUWRAREfGKL5ZjZiGH544gBQsiIiJe8cfMlsvz85CGIURERGRAyiyIiIh4ZZwOQww5s7B9+3Zuuukmampq8Pl8bNq0Ke3+O+64o8+62QsXLkx7zPvvv8/tt99OWVkZFRUV3HXXXUSj0Zw+iIiIyBnnzIYY9ja02RCjZcjBQkdHB5dffjlPPfVUv49ZuHBh2rrZ//qv/5p2/+23385vfvMbtmzZwksvvcT27dvTeleLiIhI/hjyMMSiRYtYtGjRgI8Jh8P9rpv9u9/9js2bN/PWW29xxRVXAPDd736XG264gW9+85vU1NQMdZdERETygwocB2/r1q1UVVVx0UUX8cUvfpETJ04k79u5cycVFRXJQAFg/vz5+P1+du/enfX1uru7aW9vT9tERETyTk5DEDnWO4wgz4OFhQsX8qMf/YjXXnuN//W//hfbtm1j0aJFxONmHKalpYWqqqq05wSDQSorK5M9rTOtXr2a8vLy5FZbW+v1bouIiHgg10AhP4MFz2dD3HrrrcmfL7vsMmbNmsUFF1zA1q1bmTdv3rBec+XKlaxYsSJ5vb29XQGDiIjIKBnxqZMf+MAHmDhxIgcOHGDevHlUV1fz7rvvpj0mFovx/vvv91vnEA6HCYfDI72rIiIiuVHNwvA0Nzdz4sSJ5LrZdXV1tLa2snfv3uRjXn/9dRKJBHPnzh3p3RERERk543Tq5JAzC9FolAMHDiSvHzp0iIaGBiorK6msrOShhx5i8eLFVFdXc/DgQb7yla8wY8YMFixYAMAll1zCwoULufvuu3n66afp7e1l2bJl3HrrrZoJISIikoeGHCzs2bOHa6+9NnndqSVYsmQJa9asYd++fTz77LO0trZSU1PD9ddfzze+8Y20YYTnnnuOZcuWMW/ePPx+P4sXL+bJJ5/04OOIiIicQeO0g+OQg4X6+voB189+5ZVXTvsalZWVrF+/fqhvLSIikt9UsyAiIiJnIy0kJSIi4pm4veXy/PyjYEFERMQz4zNY0DCEiIiIDEiZBREREc9YQCLH5+cfBQsiIiKeGZ/DEAoWREREPDM+gwXVLIiIiMiAlFkQERHxTILcahZyee7IUbAgIiLiGQ1DiIiIyFlImQURERHPjM/MgoIFERERz4zPmgUNQ4iIiMiAlFkQERHxTILchhLyM7OgYEFERMQzqlkQERGRAalmQURERM5CyiyIiIh4RsMQIiIiMhDrFFg5fOFb3d7ti4c0DCEiIiIDUmZBRETEK4lTkMghs5DIz8yCggURERGvWKfAiuXw/B7v9sVDGoYQERGRASmzICIi4hWrA6zeHJ6vzIKIiMj4ZnXadQvD3KzOIb3d6tWrufLKKyktLaWqqoqbb76Z/fv3pz2mq6uLpUuXMmHCBEpKSli8eDHHjh0b0vsoWBARERmjtm3bxtKlS9m1axdbtmyht7eX66+/no6OjuRjli9fzosvvsjzzz/Ptm3beOedd7jllluG9D4ahhAREfGKdQqsHL5ah1gcuXnz5rTr69ato6qqir1793LNNdfQ1tbGM888w/r167nuuusAWLt2LZdccgm7du3iYx/72KDeR8GCiIiIV6wOT4KF9vb2tJvD4TDhcPi0T29rawOgsrISgL1799Lb28v8+fOTj7n44ouZOnUqO3fuHHSwoGEIERERryRyrFlImJqF2tpaysvLk9vq1atP/9aJBPfddx9XXXUVM2fOBKClpYVQKERFRUXaYydPnkxLS8ugP9aQg4Xt27dz0003UVNTg8/nY9OmTcn7ent7uf/++7nssssoLi6mpqaGz372s7zzzjtprzFt2jR8Pl/a9uijjw51V0RERMalpqYm2traktvKlStP+5ylS5fS2NjIhg0bPN+fIedKOjo6uPzyy7nzzjv7FEicOnWKX/3qV3zta1/j8ssv589//jNf+tKX+Ku/+iv27NmT9tiHH36Yu+++O3m9tLR0mB9BREQkT1gdYAVyeL7p/lhWVkZZWdmgn7Zs2TJeeukltm/fzpQpU5K3V1dX09PTQ2tra1p24dixY1RXVw/69YccLCxatIhFixZlva+8vJwtW7ak3fa9732Pj370oxw5coSpU6cmby8tLR3SjoqIiOQ96xRYOYzwW4mhPdyyuPfee9m4cSNbt25l+vTpaffPmTOHgoICXnvtNRYvXgzA/v37OXLkCHV1dYN+nxGvWWhra8Pn8/UZL3n00UeZMGECH/7wh3n88ceJxfqvAO3u7qa9vT1tExEROdstXbqUH//4x6xfv57S0lJaWlpoaWmhs9PUPpSXl3PXXXexYsUK3njjDfbu3cvnPvc56urqBl3cCCM8G6Krq4v777+f2267LS2d8nd/93d85CMfobKykjfffJOVK1dy9OhRvvWtb2V9ndWrV/PQQw+N5K6KiIjkLnEKEjmchyeGlllYs2YNAPX19Wm3r127ljvuuAOAb3/72/j9fhYvXkx3dzcLFizg+9///pDex2dZljWkZ7if7POxceNGbr755j739fb2snjxYpqbm9m6deuAYy//8i//wuc//3mi0WjWqSHd3d10d6dW4mpvb6e2tpatW7dSUlIy3N0XEZGzQDQapb6+nra2tiHVAQxFe3s75eXl8J0wFPqG/0KdFnype0T3dThGJLPQ29vLf//v/50//vGPvP7666f9wHPnziUWi3H48GEuuuiiPvcPdn6piIiIeM/zYMEJFP7whz/wxhtvMGHChNM+p6GhAb/fT1VVlde7IyIiMoq6T/+QMWjIwUI0GuXAgQPJ64cOHaKhoYHKykrOPfdc/vqv/5pf/epXvPTSS8Tj8WTTh8rKSkKhEDt37mT37t1ce+21lJaWsnPnTpYvX85nPvMZzjnnHO8+mYiIyGjzk9vUgTxtlTjkYGHPnj1ce+21yesrVqwAYMmSJTz44IP83//7fwGYPXt22vPeeOMN6uvrCYfDbNiwgQcffJDu7m6mT5/O8uXLk68jIiIi+WXIwUJ9fT0D1USerl7yIx/5CLt27Rrq24qIiOS/gL3l8vw8pIWkREREvKJgQURERAbkI7e6gxxmXY6kPC2lEBERkXyhzIKIiIhXNAwhIiIiAxqnwYKGIURERGRAyiyIiIh4RU2ZREREZEAahhAREZGzkTILIiIiXhmnmQUFCyIiIl4ZpzULebpbIiIiki+UWRAREfGKn9yGEvL0FF7BgoiIiFdUsyAiIiIDUs2CiIiInI2UWRCRvlpbB75eUTHwdZGzlYYhREREZEAKFkRkXHKyBv1dZuPOJFRUmMdm3iYi44aCBZGzkRMIHD5sfnY2933ZggUnCMh26f7ZHTwocJCzyTgtcFSwICIi4hUNQ4jImObOGLS2mqxCQ0P/mYVs+sskZF4/3e0iMqYoWBAREfGKOjiKyJjkzhw4NQqHD5utsTE5xOr8jfL18zIWkIhGsQCam811gJKS/jMN2bIL2bINIuOFahZEZMzIHHJwhhucIKG5mQBQQGqI1QkSsgULln0Zd/3sBAtWNEo8GsVqbk67nZISmDYte+AwbZp5EQUOImOCggURERGvqMBRRPKee8ghs5Dx8GEC0SgBzP/4fqCQVHahwH6JzL9VcfsyAfTYl3F7S2AyCVa2n6NRYo2NqUzDlCmprIIztTJzeMLJOIiMVQoWREREZECqWRCRvJZZwOjKKgSjUYKQ3CKYTEIpEAIChUCJveG6jLouoxDvTM8qOD/3ZrktAcTctzU3k2huJtHYCFdfffoiSBHJGwoWRMYqdz8EdyGja6aDM+QQwQw5RFxbYBIwG6i2N3ew4BZNbYEWCLiu08KAQUQXJpDoJRU4xICeHTvMsIRTAOkMP7g/kwIGGYs0DCEiIiIDUrBgbN++nccff5y9e/dy9OhRNm7cyM0335y837IsVq1axT//8z/T2trKVVddxZo1a7jwwguTj3n//fe59957efHFF/H7/SxevJjvfOc7lJRkO60RkeQZt7PuQrZujK2tsHUr/ubmtCGHMDABKAJ804EZrq2aVOoBsv9FiLm2qH1pZxSc7EJmxsE6brIIIUxRZC8my9CDGZJNgBmSAJNVOHw4VfjopuyCSF4YcilFR0cHl19+OU899VTW+x977DGefPJJnn76aXbv3k1xcTELFiygq6sr+Zjbb7+d3/zmN2zZsoWXXnqJ7du3c8899wz/U4iIiOQDH6R1Ohvq1l9XtDNsyJmFRYsWsWjRoqz3WZbFE088wQMPPMAnP/lJAH70ox8xefJkNm3axK233srvfvc7Nm/ezFtvvcUVV1wBwHe/+11uuOEGvvnNb1JTU9Pndbu7u+nu7k5eb29vH+pui4wdmWfX7utOBsH9s6vRUgTzP3UIkywotrfApzD1CVNI1SZE7MsgJPxgZaQ/fXHwO0UGztZqX5ZgUgVOluEAabUNvmqTbYgc79v9tsd1PWE3csqaVRAZizQMcXqHDh2ipaWF+fPnJ28rLy9n7ty57Ny5k1tvvZWdO3dSUVGRDBQA5s+fj9/vZ/fu3XzqU5/q87qrV6/moYce8nJXRfJPf0FCtt4JkNY7wRlyKMUMNxQDoemYAMHZ7CAhETJbvABikb5BAoC/xwQKvrjZAr1mI4YJEhwxzFCGu+DRUW2GJkKdZtihgNTMCffJVyLbZ3PTUITIGefpjM6WFvOXYvLkyWm3T548OXlfS0sLVVVVafcHg0EqKyuTj8m0cuVK2trakltTU5OXuy0iIuKNgAdbHhoTsyHC4TDhcPhM74bIyMg2BTIzq+D87Aw9tLYSaG6mAPM/cSEmmzChELgak0m4GKgwW6zaZBHiBWZLhKA7AF0BiPkgaC/4EEyYn8ORVHbB32OyD1YAgk4hZMx+Y+fSkdmfwebP+NkpchxweFbTKGUsGqdNmTzdrerqagCOHTuWdvuxY8eS91VXV/Puu++m3R+LxXj//feTjxEREZHT2759OzfddBM1NTX4fD42bdqUdv8dd9yBz+dL2xYuXDjk9/E0szB9+nSqq6t57bXXmD17NmCKEXfv3s0Xv/hFAOrq6mhtbWXv3r3MmTMHgNdff51EIsHcuXO93B2R/JaZUXAXL7rvdxcxtrbij0bxY2oTwpg6hTIgUAf8v8BEoAISZSabEItATxl0BE0moSsAXX6I+U1WwRG0IBI3lxU9EPFDsNfc509AIk7qL4Z7bqZT8OjOJmRcTwx0HKLR7DULzjRR93WRfDfKBY7ODMU777yTW265JetjFi5cyNq1a5PXh5OpH3KwEI1GOXDgQPL6oUOHaGhooLKykqlTp3LffffxyCOPcOGFFzJ9+nS+9rWvUVNTk+zFcMkll7Bw4ULuvvtunn76aXp7e1m2bBm33npr1pkQIuNStkDh8GFz3R00tLZCY2PazKogpliwChMo+D4BLASmQfxC6C02AUJ3kQkKugLQWpAKFP4EtANt9lZub2XAeUE7YEhALJEaVXAKHdOCBPcG/Q5B9MdHxjBEtgJHBQwy1oxysDDQDEVHOBzOOXM/5GBhz549XHvttcnrK1asAGDJkiWsW7eOr3zlK3R0dHDPPffQ2trK1VdfzebNm4lEIsnnPPfccyxbtox58+YlmzI9+eSTOX0QERGRM86jmoXMFgG51O5t3bqVqqoqzjnnHK677joeeeQRJkyYMKTXGHKwUF9fj2VZ/d7v8/l4+OGHefjhh/t9TGVlJevXrx/qW4uMD/11X3QVL9LaCs3NyUyC06vAKWYsBco+Afw1ybUdus+B9yZDNGg2J5MQBRoxWYQm+9LJLIDJKtTaWztQG4CIPbwQ9ENxIFXgCPT9qxEkPatQ4rq0Zf7tdGcVktMn3cfmdJRhkHGutrY27fqqVat48MEHh/w6Cxcu5JZbbmH69OkcPHiQ//k//yeLFi1i586dBAKDT2OMidkQIiIiY0JmF7LhPB9oamqirKwsefNwswq33npr8ufLLruMWbNmccEFF7B161bmzZs36NdRsCAy2gYxPdIfjeLD/M0JYf5HLcLUFRQ70yOXAVPgVCV0VsB7Ifh9kckeHCG9LsH52cksAKZDUsBkFNyPA5gYsAseExD3g99vujwm/ODPrFdwFzhmZhVKwN9pfnRPI3c3ZbLAFDk6n/90x86dVVCGQfKNRzULZWVlacGCVz7wgQ8wceJEDhw4oGBBJC9l68jo1tAAzc3JvzXOEEQZriBhob1Vw6lZJkhoLoT3wrAL2EAqGHCPeLZBav1o52dbk/3HqZxU0BDzpW/BbEMRmZsrQACSQYMvCgG7i6P7s1mkD0/EnVkRmcdsoIBAwYLIkDQ3N3PixAnOPffcIT1PwYKIiIhXRrkp00AzFCsrK3nooYdYvHgx1dXVHDx4kK985SvMmDGDBQsWDOl9FCyIjIbMxZ/cxY32bYHmZnykn6wXAec6GYWrgWnAFFPM2FQFLRH4mR/exC5i7CUta5Dk3OZkFpw/SPapvjNM4QxJxFx9GGK+VFbBCpjuj36nc6OzvHUXplskmILLFtdlicksxElfWM8ZirBcPyfcS1VXVKSmTrpvczj3Kbsg+WSUp04ONENxzZo17Nu3j2effZbW1lZqamq4/vrr+cY3vjHkGggFCyIjLXPYIXP2g12j4LRudoKEYsysBx7BrBZZDbGJ0FUBbSWwu8gMO7wJtHVhVmjqwQQEmWcnAVIBg/Oz67a2AjPU0YQZjogCEZ+ZTeHULvj9plW0L44JEBxd9B8w2ItXBeyhCHfdgkUqUEiuQhmNmuOSGQCcrpbhbA8YBhrico5NRYUJumRcOd0MxVdeecWT91GwICIi4hUtUS0iQ5Z5lucaduDwYfyuXgpONqEY14wHe1GoRKXJKHRWmKGH90LwANDUC3SQnlVwnG7sM4BpBRkC4tBkZwucDIPzRytoQdSCkggE7eWprYDJNCRnRjjDEk52wSlydJavLkllFgpIJTWSww+4FphyZkZkHsdszuaMQmafjmzXHU5W4b77zu5jNhq0kJSIiIicjZRZEBkp/dUp2FMknZPyEGZBqGogNIlURmEKMBG6pkF3KbxXZKZIvhU0tQpNHcApzNl8F/1nFjLTou7b45isRIF5jaYCaAqZ61cBMzNSogWRVKGj3w+BAASczIKddUhmGOzOkkTN5j+e3q/GXeDoXE9mGjKzC/0VMp6N60a4/y1t3dpngTEnU+UcTwuwmpuJNTbSM3s21NefPcfqTNAwhIgMWra0sD0EEbQDhRBm6MFpthT4IjAbs2rkRLNqZE8ZHKsyww6/CcHLwC+whx/exwxBOMGCu4+COyAAEwy4rzvNDkKuywJM4NALL5dCYwBuAP4iABdgZkVECiHiChacwCFYQqo60wkawAxH2AGDr8UMRfTab9872OPofLFlG4o4m7703MGmvYWj0eSv0vmOCmEKY50hHzC/kpNAyxNPmOEIe1Vg8V4ZGQukDZFFeo+UfKFgQURExCOl5Da+n0DBgsj4l9G2OXP4IRCNUoQZdijFzihchhl2WAhUpDIKPUVwsggOlMArmCmSvwBzitiOmZroFDc6MjMKmafxmfdHMKefCVLTKXvMz03F8IOI6d+wKABX+8zy1SU+iBRAsAASvRDvMS8VzNbV0VXg6C5yHLLB9FMY7z0Xtm41244dyURQKeZX6GyBQmAGyS6fTqFp6ACUbYbOt6O0rVsHDz44vo+VeE7BgoiIiEfKyK3sIA78yaN98ZKCBRGvZC4I1dCQzC74m5spwJwNVmJmFyYzClcAEyE+HXqLTTFjNGyKGVsL4J8wtQp0YLIKrfbP7SRrDJINl9zdGzMbMWVyBrULXD+HMKeok+z36oJfRKAxAr/ww/0RqOiFkhhE/KkMgz9hFpkKFNgFj85flgrSGjT5j/d/+LK2lRns2e94zio4malNmwg0NhImVZdQVYipc5lhb9Mw2YQpptNmImSaaAWuAOphxlfhyI4dHN+0Ce64Y9Q/ytmgjNy+WGNe7YjHFCyI5Cqze54z9LB1K367mLEAM/RQBkxwgoSPYYoZKyBeCtFqaC00AUJrAbzth/8AXu7ABAgnSS9oPE5GfwWnjSLpg6buGRKZrRedCjh3wBDBBCITSHZkbIvDyxEoD8BVYfhw0A4Y4hAJmocFu8ywhBWwhyQgVeDobIfSD12C9CDBytzdwRjPgYIzlNXQQNAOFIowx7usELPy6MWkeltUQKzErETaG4Euu213MAEltVD+TZj6Y+j+3vdonz1bhY4yaAoWREREPFJKKgYfjkHNEjoDFCyI5CrLolA0NBByDT2U2lvhJOA+ktMj46Vm2KGnDP5YbDIKf/SbYYcNYDIKf8IkDE6RGoroAKyLwVcMBUXZ98tXnLr0F6VfB7A6zJY4BYkO6D0OnR3w/mE42WWmZpZihiQmAAnYUGSmVF4VMFMqz0+YdSMASgqgoMsMRyT8EPSDP0qq50IJ+CalD0U43Rv7KClJv95f5iAzqzBesgzuf0/ONEkysgoLMUMQFWZLFJlhrFgE/lxi/i11Bewpr3GY6IfyicAMe2rf6dbbkGEpZ3wGC+rgKCIiIgNSZkFkuJwzs4aG9DPBxkaCQCFmnYcy7LUeXMtMJypNNqG7NFWn8P+F7YxCHJM9sAsMOU6qZqErAsFpEJhkLn3FqUwBmAyCw1dM3+5LAddtkDqvd1o5xs2bdm2FY29B6wnz3scx2YVp0Oi3Cx4DcJUf/sIPREwNQ4Ufiu3COisA4RLSuzmWQOC4ax0IBtHA5mwrcnT+HbkyC06BrLN2iLNmSLJWIWiOt1vQsqe6xmBCNxS9C9wBv3sbTr3wglagHCFlmP/bhqvHqx3xmIIFkaHK7KXgzHpobEybVFCFyeL76jCBwgxgInTXmiDhZBG0hkyg8F4A/g5oc4YZuoATmC/qJiBRa4KDwioTKAQmga+M/pOD7qDAmeJQAIkwWEGzJexLAF8M/DFzGeiCSAVEPgG9DfDnt6D192Zf/JigoQsaQ9BYBEcCcEMQLvOblHdXwAQN4AoWoiQDhoKMIkdHn3bPmUMR2X4P7u6O4yFQcAedzmYXyUYwwxCB6Zh/S1Ps58TMFuhNBQwVnaYYNdgFhScg8DbEb4e93/ymafcsI6YUU8w8XN1e7YjHNAwhIiIiA1JmQWQ43Gd+W7cSiEYpgOTwQylQdhkmozCTZEFjTxlEJ5m1HlpDEA2a7owvA20nMEWMrcC7pIYe/B+H8Idc2QT3Qg/9NVhwZxUKIFECiYiZVxePmC0RgV5/ei+GAKZpQqQFClqhYAIUzIHevfD+Vvj/NZl9moQ51e2AlyeYt2jzw8fCJrsAUFEEkQgEnfaCTpFjYaqTo3uZavcnSONkcjLXiMiWVRjLGYaMNUTcQ1pOkWyhu6+CMwvWWcQrZjIJQXsp8UAvpr/FV+2hh2XLlFUYBWWkT1Aeqq7TP+SMULAgIiLikTLMCcNw5VLvMJIULIgMVuYUSbteIWRnFZyzv2KgeBKmYY4zRfKc1BTJ5kKTUTgYMNMjk90Z3yeVUWgFTkWg4EMQqQdfBak/I87ooVOQaC/sYHXYRY3OQg8h+7LA1CbESlJblz81HTNORnbBD6U1UFwD4fftLEOp2ZeWr5rn9WKyC3YzpzcnmaeXA5cFXcV1Ba7MQgl91ohwL1Pt1Cw4nzC5TPVAGYVsWYWxmF1wsgoZ/76c1Umdf1vJwsaJ9F3h0/XXPGAvC86j8IdDdlZBHRtHRa7BQi7TLkeSggWRwchWeGYvDBXBFDQVY/6gh5wOjVNIDj30FkNXqSlmjAZhayC13HRq/WBMkNCOCRTCV5ovaN8E+s5gsGctxI+n90tw+ikEJtkBhp3gdwoa43Zbv1P2+x7EBCmngE775QuBD2JXaFaaL96iCIQCEP44tP7cvGzE/sAd0FYBTQVm0alav1nGustevjptYSk7WPC7ZkQ4gULA/tmJW5IBgxM0ZHICg8zrYy1g6Gc580A0SggTLFRgF8rOxgznOH+57eEHIFnoSAvwAsTfNj+2l5SoU6PkTMGCiIiIR8ows1aGS5kFkbEss6OefebnrPmQXG66DpNVsIsaM7MKXQH4dcBkFNrApPNbMdMkO+ztFCaj4PRTSHINO1jt0PNbSLxrMgpOdsFvT620OsxrECA5fGEFIR5KrS3RCuyz3/sEdmYhBgShGZNduAzTjpFp5jVCH4XYYTjVZPbVGWJwPk8Gp5tj5rSrgGtLkN5zwV346MNeL6K5OTWV0skaVFSk9wpwZxXc8jXD4OynewjClVXI7P7Zb2FjKyaF0AI0gPU2HLbvOlVSAvfdp54KoyjZC2OYTtt35AzR1EkREREZkDILIqfjHk9uaOjTfKlPVmE2fbIK0aCZUhgNmnH9dqDRySo4TZiczTchfS0H9zrU8eMmmxA/Dt1vgXXC7IgPu+CtyTRwsqaZ5wcvoM+kxASmTVwz8F9Ap/PGzqB3EP5QAZ32n4ePYNajDpVAqDq9YyQkl78ut48DmFUOwczC9GdZ/MGdWYhjsgfuro5O/YKTbfABVjSaatbkrklwsgzO9cOHs2cT8i3D4M4muDJXTrdGZ5XSCiA0nfRaBedX9XvgANAA8UPmn9Ofgbarr4abbza1Cvn2uce5YCL173+4z89HngcL06ZN449//GOf2//2b/+Wp556ivr6erZt25Z23+c//3mefvppr3dFJHdZChsza/aKgYBT1Oj0VKgwgUJvsb1MsM/0VWgJwhHgF3HSl5t2Zhi4v9ed4kVfsbmMHzdb7DD4TqRWrHEPcp7CBAxd9mpNvmIIFIC/y3RndCQwww6dTi47SlqwQAROlKRmS0Cq2yOu9/Wnfq4FpgIV9uJSkTj44q6XzcJpQJ1wXU8OPZBKyTrLWCcwQYMFEI2ScA9JuIcm3Nczix3z4cuzn0CBw4eTQajTLjxUiAkUwEQDJZhfVxTYAdZOE286I1ndV18NTzwxah9F0gUts+Xy/HzkebDw1ltvEY+n/uI1NjbyiU98gv/23/5b8ra7776bhx9+OHm9qCiXchAREREZSZ4HC5MmTUq7/uijj3LBBRfwl3/5l8nbioqKqK6uznyqSH7Jcvbnj0bTUugRIDQJM/RwMcmFfeKFZqngbjtFHw2aDEMbZhgiObLQY2/OurQBoLfDZA/8p8xtTvFivMm8YbH9xhUk+xwkq4+c5atPdsHJn5vbIkXgnwTBKAQqoCBknpc2GTxIKqtgbxMwmzPRPxiF2EGzb05DCXt/rsIkVS7BLFwUiZvsQlpmwckulKSGHNxZhV7SJ4dmFnqlZRfsLR6NmuEJpwDSnUnIHJ5wivzcwxfuy9EyQFbBPQRhr+ZhRDHDDc5yGXZBY+dOM/O1HTg1ZYopZtQ0yTNqvGYWRrTAsaenhx//+Mfceeed+Hyp//Wfe+45Jk6cyMyZM1m5ciWnTp0a8HW6u7tpb29P20RERPJNIJGqWxjOFjhbahbcNm3aRGtrK3e4Ood9+tOf5vzzz6empoZ9+/Zx//33s3//fn7yk5/0+zqrV6/moYceGsldFcku4yzQvTKis5YjM0hNabO3WMQ0JApadlYhADG/OQNsglRWoZf0AfsImGZLTdDbZG4vcL2Z62yeUlKZBYdzX7F9X6udXSg8D0Ih070xUmUeNwWYEDS1CcmMgv0iUyJwOWb6ZClQ2A7hY9DxFlR0mYyDs/52EfwFJrtQEkttwV57fYIs9QruaZNuTj9Kp1kT9J9hsOyXTmYa7ALImDvL4N6cwkd3tiHzZ/el19yrlWY0YHKyCkHSV/aIA/FOCBxwvU4LWMdNOUkzJunAAw+YdR/yoR5DxqURDRaeeeYZFi1aRE1NTfK2e+65J/nzZZddxrnnnsu8efM4ePAgF1xwQdbXWblyJStWrEheb29vp7a2duR2XCQzVdzaCtFoMlhwz4ZgBqYAzQ4UEv7UUsFgAgVHtl4EySLBiOuFnSaNzliHs7nT/6X0bezoDiycwOHEzyH+UfPYZHlQlbl+LdBsBwzvl5ihiULgOuBczFBHWTsUHQZOQMFe81kn2O9fCTMDsACo6EkFCpEes5hRchjC4bR7JpVid4YjnK7T7gDC+WiZ/fKdx/a6fnaKImOYzo+JaBSruTl5O1Om9D88kS2A8PKLNzNQyByCaG5ODs34XJ8N+7Kg015m2r7egV3j+JnPmDbOChLyxngdhhixYOGPf/wjr7766oAZA4C5c+cCcODAgX6DhXA4TDicywrhIiIiI0/BwhCtXbuWqqoqbrzxxgEf19DQAMC55547UrsiMjwZnQDdRXk+UnV/7qwCwVRWwReHWCDtJSh3X3FOr52OfGWYoQn3sEQAkx0oos8wRHlBqq8BmCGOtgL7Me7uiqXA7/8JwovN2g5FXaZYsXAiXFFmKhPd71uAySiU2ktVR/4EHIeun5lGjjX2/RWwKGCGHyZ2p2cVAvYQRFqPBWdtiGrTZDJyPJVRcBIpAy20nVlg5bSLcHWhoBczI9RZY8LJKlhArLkZq7nZfExnmAJMZiFzuqV7CmauBsoqtLZCc3Pyczp9JpzP43xu559RF9CNGYKI/T//jylolLyiPgtDkEgkWLt2LUuWLCEYTL3FwYMHWb9+PTfccAMTJkxg3759LF++nGuuuYZZs2aNxK6IDJ3zx9z9s33dGYZwtxjAGfKHZLCQsL/ZYvaAu/MHpMxvAoY2p0+BEyi4OxQ51yH1xe0ECBHTzwDMl7Q7WABojJitqYLUAlVOPXDL/4G2X0LRf4fQByDUAv4Zpo4hYY+f+BOmH0OkBcLvkexF3fE8FO81wUIFlJfCrZhahYtiJlgojoHfHSj0kD4EEcQEVlFz6Ssxu2B1puIUZ+QF+/j6JtnHdzapVSudQxSFwgNQaM8M6OlMzajIHJ5w2kc7sygS9iwKC0yvhsOHU8GBEzxkLlKVi/4ChcOHk0MPPlJBk3M8YqRqNrqA7pISM+xQX68WzjKqRiRYePXVVzly5Ah33nln2u2hUIhXX32VJ554go6ODmpra1m8eDEPPPDASOyGiIjIqArkOAwROJuGIa6//nosq+8nrq2t7dO9USRvZWQY3DMhnESAr5BUZiHj/yafnSUIJiDoM30Hyv0m6/+LCKnTaCejEHE92T1EUWqKCJ0hjJn29uG4PT7qSlsuCJkT90bsDEMEfjGJVF+GliZo+SfomgOReijqIVlZmQiDP2bvTCtwErp3QefPoLrLDD9Uw8wCk1X4GCajMLEbSk+lMgo+u3NjoDe1X8njU2FegxmkMgxRO5vgXDoZhGr7cgYmoxHJOMYxUosoHYDQAQg1mNkD7skmvaQvSpUseCTVqyHhdIN0ztad33tmoeNQMwzZCmUz+nZkdqnMnDwSB2JTppiujMom5L3xWrOghaRERERkQFpISiQb93i1q17BT3p2IXkWnCWzAKmzjEjc1C9MjJssAUBjMbT10neqpPPidn3CTPu1nIzCeZhpitVdJlsR7E1lMXojpq/DjCDUF8B7AVgEPH4etDn9FyqA43vNFlucWgrb70xm7IWurWahqpITcCHmLL8aFhXA5zGfo7oLKnqhoAtCp1IZBX8i9bM5COazJE+Zp9iXUdfmKKFvZiGCWW/Dziw49SD+BPCefd+U1OOLt5jxfbfk8hak93FwTuL8mAxDvLHR/L6duoVsBptdOM10Sb9rKm62/YwDzJxpahS0INSYoQJHkbNFxiyI/vSXlnO+MMHMDIj5TK+FiH3bXwRMYWIt0Fhgih3bMT0Yyu3b3ZdlmCChHKi2iwkreqH4fdPLwN+Ter94AZTYTaF6i6EtBFNCUBuCDRXwsjObwlnOsOn/mDdPREzAAGaxqtIumITZJpjL2mL4MmbGQ3UXlNu9FJwtGRy4OYGC+3oQ8+XuBA/uYKHC9ZyS1GX3OeazJVwNF/w9ZjFMvzM8cbF5rcAO05cgs2+Dc+msaEmW+wDizc3mB6eJ00Cp/2xf4O4gwbnMXJDM1VfB4czaiIPpCXHzzWZTkDCmaBhCREREzkrKLIj0Z5AZBiBtsSRf3Jz1WqXm5+KMirUL4lAeMG2fazGXTlbBySQ4xYxtmO6IwYTJJjjFhKFTEHkPc1belXrvQAQCJRCOmBmRoSIoLwNK4MshmFoA/1ENTU6GIWA/v6sLeuz20gFSHRonAJVQG4AbMMMfE3ugotNkE5wpkmnDDpjpo1YA4gHw+cHvZAuc5biddasgNR3SeYyzlZgMSSwCHZVmIS6nG6YztONPQEEAAs5K23Y3zcChvrNRM7tDQvbeDX7s6ZSZGYJMzqJVmQbKKNjDD+7eCZDewpqZM01GQe2bx6TxmllQsCAiIuIRf46LQflVsyAyxmScObqn3OH8HMVM3SsheXbv74FAwJx1xwvs7ILreRG/qckrsWsXppJaM6Lc/rnNvn0m5gzaqVMo7ITwSSjosN83mrE53SRLIGhv4ZMwrRbei8FtEaj1w8ul8AtnnQhnjqG7IVQFJrMQMVmFWzEZjok9qVoFJ6sQ6M1er5BsTlXgKnoMuWoMnMyCc+l0eIykMgo9RdBZCO+FU5kFp9FVJGFmVJYUQHGv/fkrSFsC2ynZdOoX3JmG5O/QJVnPEI2mlrrOlj3o73boW9DY0JBcKCpt2i2uGgX750RJibIKkpcULIicjmuuvVNB78zjL+oEXwvmS9r1f1MACLWnvvTABAzOKpROGv28jEWmjmACBqegEVwFjR0mSAifxMwCaAEa7MsW+8HOLAKnl0GFGZYoLYCCCvMFG4lAedAEARuqSQ0NuIOFIiBkZm4sIrWiZH+BgnM25Cyi5QQKiVD67AVfHAIh8IfsDo9gggX72MVcQw+xCHQEoTUELRH4HelB1SV+cyzDxXAqbhc7Yj63LwqRLEMRuH52AoUuUl/cmbMT+nRxPN3QVLaZDzt2ELA/YmZbZyfjHAfiJSVm1oMChTHNl+in2HcIzx+K7du38/jjj7N3716OHj3Kxo0bufnmm5P3W5bFqlWr+Od//mdaW1u56qqrWLNmDRdeeOGQ3kcFjiIiIh5x6ndy2Yaio6ODyy+/nKeeeirr/Y899hhPPvkkTz/9NLt376a4uJgFCxbQ1ZU5wXhgyiyIDMQ5m6yoSGYV3FPyujDrE1CCmboH5ky5BELBVFfDWMQMSYRDqZ4LsYysgnvqpKM61k9W4T1gM8nMQrzT1RCykFRmYQYwG0IRe4pjBZRUmOmVl0SgqcCevllE+lBEcSqjUIvpqzCx256m6WwZZ1DurEK8wFzGIulLdjudHYNd4Iukj88m/Ga6Z7zA9Ivo8pvhh8Nh+AGp6aXY+3QE+B8F5ngGQ/Z7VAAL7fey154oOJ5avRtSLS0gtQx0D6msUfLM3xmGcLiv91fYmGXp6RDmD617ETL3v6UYyiqMJ/5EbnUHQ33uokWLWLRoUdb7LMviiSee4IEHHuCTn/wkAD/60Y+YPHkymzZt4tZbbx38fg1tt0RERGSktbe3p23d3d1Dfo1Dhw7R0tLC/Pnzk7eVl5czd+5cdu7cOaTXUmZBJBv3GaSdXUiQ3jinB3O2mswsgDmjn0hyOqMzldHJLsQSUBGwaxfsAWunbsGdUTg/kVryubAzS61CI7AZOjvTl2gOAP5OCB2yx+wbgAPAZ8yyD0X2Wc+UUnNG/g/F8AugKWA2x1+Qqp1wGkGVxFz1Con0ZlCZWQWngVJ3kckQxOzTkmDCfPZCd9FjPPXcWAS6A6aQsSUCLUHzUd8E2pw1qIHGgMmIXBUydRjFdhMqgEA1cDOwyxwrX9ReoTKjY6RlZxxa7c+cIKNewfn9u/9NDFSz4FpJ0j1FsoC+tQqQkVWor09lFmRMG85QQubzwayl5LZq1SoefPDBIb1WS4spZpo8eXLa7ZMnT07eN1gKFkQyub8UnKChogJKSrCiUeKYP/pdmC+B0kP2l4wzG8HZYpiZCVEIRszMhHiheblgUWo+dlccIvYXZDBhFyHGzZdzRa8JEtIChV3ADmjvNPvgFFs6PQKcNai6gMhxiGyxv0DrzY0Re8bG5F4zHDIlaAoFW/2pNP95pPbB2Y+CrvTlp93p0sxAIRYxwUJrQfoMBjCfucJnLsP2H8a4HUxEg+ax74VTgcLLQNsJ+wM571lghlAemAC3RmB+EKbZi20FpmCGIyZiIgGnD4XTKdKeweJrgUALFNknWD1kKXDMJttwROay09Fo2tBDZlFjzPmdTZmSChTU0nlc8CpYaGpqoqwstQh9OBzOcc9yo2EIERGRPFNWVpa2DSdYqK6uBuDYsWNptx87dix532ApsyAyEFeBI9OmEW9sJECqSK4XOIWd8ndnFpyllasxZ7l2V8JACRTFoaDYnsJYbKYHdrkyC860yuKYPfzwZ1JZhT3AZug5BCdJnWzHSZ/61405E3CaIU54wX7+QvOkULXJEJzrTy0+5QwXxHzm/Z0MRzievgZEoJe+6ygXpIYiEiF7CCJgv66dWYj5UtkU53M6n9l5Xyer8F7AZBU2AI1dmAyBM97iiJjhiAdK4RdBuLUcLo/ARZhj7kztDHSS1mEzOVU0ChyGwgZIdKYSEH1+/wOt/+D+2ZVVcC84BqmhKyerEAMSV1+dyipo6elxY7QLHAcyffp0qquree2115g9ezZgaiF2797NF7/4xSG9loIFERERj4x2n4VoNMqBAweS1w8dOkRDQwOVlZVMnTqV++67j0ceeYQLL7yQ6dOn87WvfY2ampq0XgyDoWBBJJvMAkfntpIS4vbSwk7dgtPPqOA4hI5DwCl4dDo7uusYKkyhYbgCCorMmXrINcYP5g+NcxZf0AE0k5oquQM6j5uswinsAjnSu0s6Z7UBUifS/k4o22LG6Lka+JjJhBRjCgMjrqLEuN+0q3WmOTpn6MGujLN0GPAvSCxj8N/d8z7mqvSLBVIrc3YFTFahCVOr0NiL+bAngQ7SKzkDmIH/MtOR8uVSWBSG1eeZgkznPYOWWcbb+SzBLrO2Bq32/s8A/9tZxmQHUz+QuQ5EFk6dQtpUySlTzNLT06apTkFysmfPHq699trk9RUrVgCwZMkS1q1bx1e+8hU6Ojq45557aG1t5eqrr2bz5s1EIpEhvY+CBZFM2QocnaChoiIZLDjfhR2YgsIQ5os51GnPSDhuf3E7fQ+cbSEQBX8JhEtILpgUtxsBOMGC3/lCawR2gLUl9b3ZiRlqcAKFzBMZP+ZLyQkYWjHfs8VvQ/EBzAyJhRAImqGRRJEdKBSk90RwujM6vRWIkp6rt5ecDriGB/x+c3NxxAQFXfHUbAhIFXEGXWdQXfYMkZjPxFdH7I0u+wC3korM3AGDM0RRBtgBw5sT4Naw6RPhLtis6DGFmhO6IdFqr5ztWh7bme2SXHDKPQRxuqGILNw9OZI/u/sp2GlhGV+8KnAcrPr6eiyr/9WnfD4fDz/8MA8//PDwdwoFCyIiIp7Jp5oFLylYEOlPZobB2VpbSUSjJEhNgXNq75yTXSdL7sek7gOHoOCQyUD4oiQ7K2JnFpyplWlFeFHMWfM6iB8yHQxPYjIKzsPiAFOmpJ8B2z0hOHyYeDSa3KdkbV8nlG6B0AHgq5ihkRKzBV3DIWD/4YrRd/qhexjC/isSCNp9JewFtIJd4C9KTat0ZFvOOhiyCx7tQsdyv/m89JDKLjgVnb0RCEwyT25vMmMpxZiFr0qh7RT8oAB+MAHKC0x2oRa4NQQ32vs62c7kBDCfyT391HKOqWOIwwTO33qn4NTJWDBzZiqroIJGGWMULIiIiHhktIchRouCBZFs+uvWN21aMrPgXl/A3ZnPqWdwhumdLEOBvZVtsTsrNpBav8GZbgnmrP09TF1BA3QeStX39eDKKpSUpArkMqf4tbaa+w4fJt7aSndzc3L83DmLLjsExeswBY8Xk8xyJJeQxvVm7imHmfMLg/btTv1CxGQZiEAolJ5ZyPxD6DRyCthFll32Qat1shu99ms7xRoFc6D4UvBXmbRF7DD0HIbob+HECZNd6LUPdAe0RYAQNEXgzVKY6Ydqu9gkYL92/FD6sY1D32OaeWyd2zI7PNoSUVcxxMyZ5nfhNF5SVmFcU7AgcrbJTD871+0/9vHm5uQSw+6CR38/tzmtfxOYzooFxyHUQCpYcHqk2B0GOQA9x81IhDPzwWlxEIdUx7+BvtTsL7R4a2uy+6Qzc6IXiOx0zZCYAUzBBA3uQuls/Qlc34VAMtAAUkMTQfC7hinSOLdHzNCH0w57ol3keV7I1Cy2xUkFC5GFEKnHRAQRIAAF55kAIlIPPb+FP2+FkB00OOMvzmcpNYWOJXb/ClqBZjPc4RzfGJghCHcQNljZfgcVFeqlcJYZrzUL6uAoIiIiA1JmQeR03BmFw4dTZ4gVFSQOHyZmT6V09JdhcKYxOifLEaCgEwrehtDb9tLSgNWZyr47tX3u0YA4pFLbA03vc9Lkra1QX0/i8GESjY3JaXzOIkpFh+xFlmbb2wx754Kkhhyc+ZcHSGUWMrMLmdzZhszbnUxKxFwPxsB3juk9EQmaaZXtTiOLDuyswieACRCfALESsIJQ0AqBdvAVm+YV/iLo+N8mheNMsSwG/GYq5cRuKI/aXTGbgRdM0sIZgkhkDu0MVn9ZhWnTtOz02cbdh2S4z89DChZERES84nTeyuX5eUjBgshAsp0pOtkFuyYgcfiwud9d7GYXuDndFJ0pdE7WwVniOmQ/JoTJMjj3OVMxnR5ETkYhWdSYWdiYuZ/ZivCcplINDcn9a7dfl+NQuINUxmAGqayAc5tdR0ELadkFqzN92WW3QMb1ZJMqp07DKfDEFEY647UxZwVMp+lSZKF5Qtd50FWdqoQsqoDwe2YjDsFpZmcK7J2JkJxaOT8G53dA0bvAYWAT9Ow0Da6SJ4ND7aiYbXqt87P79yQyxilYEBkM9xeBEyg4W7aZE/b1hPMY1+JCkOoUGCN9uWtIDRE4j3Hm/ydOFyi4L/vrQum8x+HDWM3NdLv3uRMKD7iuZwsWGsA6nsrwuxfUOt0JkfP5/J2mk2TIHXhUmM0XN3HA7yDVW+F9oOI8iM6A9jIT4ThNLSgyLSgLWsEfgMSp1OpZBfbrFsFVEZjWBiUt4D8CbAVrY6qzZbKwMVvw5eivY2O24+/8nuTso2EIERERGZCCBZGznPuMPXNzZP7sWmgo0doKzc34SHX3C5Del8Hd/c+dYeiTVXCfuQ5U4OhwsiGu2xNAzN4fpw9E4LgZEkn2fXCKGFvMFj9usiBO10pnH53sh7PvDnfhp7vvBEDcyWSUYAoru8ztMb9ZDoNTmOxC+OPQMQ3+XAYnMBvAOfabxiOm2JEu6PllqnCyAqiE2mL4B6D8BATt7AgvmISFk6RIuKdLZsqWUcjM2jg/O78bDT3IOON5sPDggw/y0EMPpd120UUX8fvf/x6Arq4uvvzlL7Nhwwa6u7tZsGAB3//+95k8ebLXuyLirWx1AUMIFpwhi1hjY3JIIg5pwxOQ6oPgrFSY9kWWrQlT5n4NxB00kAoWukk1jwq4V86E5BCEdTy18KMTMCRcm7O/mXyuSz+p2SBgWmGH7ECEGdBVato9JxeR6sLUIZysNEMSTZhZDEX2C0wCQq0QeBfih8H389TwQzWUF8MjwMwOiDRjAoV18Ofj8Gf7cycy+yq4j21/v1voe/yHM4tCxp9xmlkYkT4LH/rQhzh69Ghy27FjR/K+5cuX8+KLL/L888+zbds23nnnHW655ZaR2A0REZHR5RQjDXc7m2ZDBINBqqur+9ze1tbGM888w/r167nuuusAWLt2LZdccgm7du3iYx/72Ejsjoi3MhZtSsseZHZQzCyC7Kfw0S0BUGKf1mcWNPY39HC6s9nM+50ulA0N+Ox96MIeisBextqVWYh3ptZ0cmZoOMWNlusys/mc+7O5+004mYUQrsxCl8kqHA6bBEIyBeGvSr2IPWOEQmACUPk+RI5A/LfQ9U8mo3AeMMkMPzwCXNUB5x0D9gCPwonjZgii296ydsI8ncxi0qE+X2SMGZHMwh/+8Adqamr4wAc+wO23386RI0cA2Lt3L729vcyfPz/52IsvvpipU6eyc+fOfl+vu7ub9vb2tE1ERCTv5JJVyHUIYwR5nlmYO3cu69at46KLLuLo0aM89NBDfPzjH6exsZGWlhZCoRAVGdH35MmTaWlp6fc1V69e3acOQiQvuIses41zuzMK7umM7sLHbK/pfm13NiGzVmGoPQGcfXB6RcyeTWzHjuTaFc4UzoJOU08AqamR7m6SzqJUiYwtm8xOlhap7ILTXyJkd4YM2ynYRqfRRMi+DPRAUQhmktE/4b+AJujaap44AaiAq4rhHxNw8Uk45yiE/gt6vmwSGKfs9447tQr91Rpk/s6yFZAqqyCZxmnNgufBwqJFi5I/z5o1i7lz53L++efz7//+7xQWFg7rNVeuXMmKFSuS19vb26mtrc15X0U8019TpGy3OV/+g523P9yhh/5kaS7lFF324JoZYT/cmfng9CPIFiwkexRkkXAHSNFociYImODDj72g1h6zHtS0Yqh1hiIATv07lNQDM8wylgBFCYi0ACfNFI3AXpgGVMKiCvhGD8xog/JDmKkVj8AfsQszS0pSCzsNpR9CZtCgQEGyGafBwogvJFVRUcEHP/hBDhw4QHV1NT09PbRm/JE8duxY1hoHRzgcpqysLG0TERGR0THifRai0SgHDx7kb/7mb5gzZw4FBQW89tprLF68GID9+/dz5MgR6urqRnpXREZPf2eamcMSAz0v23CDV2ew7l4Nra3JBbGc4YEC10OdOkN3N0kns3DarofQZ/povLk5uagWmAxGqBOKN0PkVpheBGurYX4NpqdCaxME/gvK34XyCGZsogc4BSf/EcqbTG+FSmAafPsUTPsThH8HvACd/9sMP7QDXH1134xCf1MlBzp27ucpsyBu4zSz4Hmw8Pd///fcdNNNnH/++bzzzjusWrWKQCDAbbfdRnl5OXfddRcrVqygsrKSsrIy7r33Xurq6jQTQkRExj4tJDU4zc3N3HbbbZw4cYJJkyZx9dVXs2vXLiZNmgTAt7/9bfx+P4sXL05ryiQy7gx0ttlf4Vx/RXQjwckutLbC7NkkduwghqkncKY2uusSnBOmtI6SAzWIGuAsPW43g/JhZkP6Af8h09GxLAJzQrCpEm6+DJPm+N2bppNjsm9kF8TfAX+TKXSsBi6EVwNw0ZvAq2CtgncwGYUugM98JntB4+nqTTLvy5ZZEBnnPA8WNmzYMOD9kUiEp556iqeeesrrtxbJTwMNSfR3fbS+gFwBQ+Lw4eSXeDepmQvuYsbkgkunCxROw3kfxyn7svAJ4A6o6IJ5M2DHZLhtDjR1/R945/+YlZ+CF4PVAVYTzAKmwecnwdfa4bwDYF2VChJOgRl6mD0b6usHt6+D6Vnh9bCQjB8ahhAREZEBKVgQkZxkdn5033amuLMD06YRa27GT+rvlXup7DjAzJmDH35wp+0zP29JCTF7GmVyDQzgz2/DOQ8Ad0DJfLjyUvj5FJj2UeBPmILHU2adGUph0Wz4bgdc8Evge/Dn/w2HgYSTTRioj8Jgjo37c2TeJnIWUbAgIiLilQ5yK1LsOv1DzgQFCyJnQr6cmbo7RNpTG2ONjcm7k2s+OAWNztl6f2fY2aaEZrtt9mwzjRKIt7bSY2c0eoBjnVCwBorWmNrF8wvB+jHELoaOqRArgIIuKD4OgSI42Al7AZYtgx/O7tvEygv58vuS/Bclt6EEBQsikpecgGH2bMAMB/Q4AYPTR8FJ59uPGVR7ZPfrZ7bDzlLMmaiooMsVhLS67qO1FY5iNvfzf56xDyIyIhQsiIiIeMVZQGW4ur3aEW8pWBARc6aeuVSzq/Cxz8/9vQb0XbI72zDEUHsUnIlppSLD0YFpBzJcPV7tiLcULIiIiHgliulIPlwKFkQkr7mzC05tgnP7cLIA/XVzVGZAZMxRsCAi6Qa7ZPNQKECQs0UHuWUHchnCGEEKFkRERLzSQW7frHnawdF/pndARERE8psyCyIiIl6JMi4zCwoWREREvNIBBHJ4fi6tokeQhiFERERkQMosiIiIeKWD3E7DE17tiLcULIiIiHglyrgMFjQMISIiIgNSZkFERMQrHYAvh+dbXu2It5RZEBER8YjVBVZnDlvX0N7vwQcfxOfzpW0XX3yx559LmQURERGPJMit7GA4z/3Qhz7Eq6++mrweDHr/1a5gQUREJM+0t7enXQ+Hw4TD4ayPDQaDVFdXj+j+aBhCRETEI3EPNoDa2lrKy8uT2+rVq/t9zz/84Q/U1NTwgQ98gNtvv50jR454/rmUWRAREfGI+wt/uM8HaGpqoqysLHl7f1mFuXPnsm7dOi666CKOHj3KQw89xMc//nEaGxspLS3NYU/SKVgQERHJM2VlZWnBQn8WLVqU/HnWrFnMnTuX888/n3//93/nrrvu8mx/FCyIiIh45EwUOLpVVFTwwQ9+kAMHDuT4SulUsyAiIuIRr2oWhisajXLw4EHOPffcHF8pnYIFERGRMerv//7v2bZtG4cPH+bNN9/kU5/6FIFAgNtuu83T99EwhIiIiEe8KnAcrObmZm677TZOnDjBpEmTuPrqq9m1axeTJk3KYS/6UrAgIiLiEYvc6g6G2u15w4YNObzb4Hk+DLF69WquvPJKSktLqaqq4uabb2b//v1pj6mvr+/TnvILX/iC17siIiIiHvA8WNi2bRtLly5l165dbNmyhd7eXq6//no6OjrSHnf33Xdz9OjR5PbYY495vSsiIiKj6kwXOI4Uz4chNm/enHZ93bp1VFVVsXfvXq655prk7UVFRSPenlJERGQ0jXbNwmgZ8dkQbW1tAFRWVqbd/txzzzFx4kRmzpzJypUrOXXqVL+v0d3dTXt7e9omIiKSbxIebPloRAscE4kE9913H1dddRUzZ85M3v7pT3+a888/n5qaGvbt28f999/P/v37+clPfpL1dVavXs1DDz00krsqIiIi/RjRYGHp0qU0NjayY8eOtNvvueee5M+XXXYZ5557LvPmzePgwYNccMEFfV5n5cqVrFixInm9vb2d2trakdtxERGRYRivwxAjFiwsW7aMl156ie3btzNlypQBHzt37lwADhw4kDVYGGhpThERkXyhYGGQLMvi3nvvZePGjWzdupXp06ef9jkNDQ0AnrenFBERkdx5HiwsXbqU9evX89Of/pTS0lJaWloAKC8vp7CwkIMHD7J+/XpuuOEGJkyYwL59+1i+fDnXXHMNs2bN8np3RERERs2ZXkhqpHgeLKxZswYwjZfc1q5dyx133EEoFOLVV1/liSeeoKOjg9raWhYvXswDDzzg9a6IiIiMKouhd2HMfH4+GpFhiIHU1taybds2r99WRERERojWhhAREfFQvmYHcqFgQURExCPjtWZhxDs4ioiIyNimzIKIiIhHVOAoIiIiA1KwICIiIgNSzYKIiIiclZRZEBER8YiGIURERGRA4zVY0DCEiIiIDEiZBREREY+M1wJHBQsiIiIe0TCEiIiInJWUWRAREfHIeM0sKFgQERHxyHitWdAwhIiIiAxImQURERGvFBeDzzf851sWdHR4tz8eUbAgIiLilfJy8OeQtE8kFCyIiIiMa+XlEAgM//nxOLzzjnf74xHVLIiIiMiAlFkQERHxSnk5BHP4ao3FvNsXDylYEBER8Up5ORQUDP/5vb3e7YuHNAwhIiIiA1JmQURExCvl5RAKDf/5PT3e7YuHFCyIiIh4paJiXAYLGoYQERGRASmzICIi4pXycgiHh//87m7v9sVDChZERES8Ul4Okcjwn9/V5d2+eEjDECIiIjIgZRZERES8Ul4OhYXDf35np3f74iEFCyIiIl4pL4eiouE/P5eZFCPojA5DPPXUU0ybNo1IJMLcuXP55S9/eSZ3R0REJDfl5Wb65HC38vJhve1If5+esWDh3/7t31ixYgWrVq3iV7/6FZdffjkLFizg3XffPVO7JCIiMuaMxvepz7Isy7NXG4K5c+dy5ZVX8r3vfQ+ARCJBbW0t9957L1/96lfTHtvd3U23azpJW1sbU6dO5Wc/+xnFxcWjut8iIjK2dHR0cOONN9La2kr5MM/cT6e9vZ3y8vKcv5ecfW1qaqKsrCx5ezgcJtzPlMyhfJ8Om3UGdHd3W4FAwNq4cWPa7Z/97Getv/qrv+rz+FWrVlmANm3atGnTNuytqalpxL7XOjs7rerqak/2s6SkpM9tq1atyvq+Q/0+Ha4zUuD43nvvEY/HmTx5ctrtkydP5ve//32fx69cuZIVK1Ykr7e2tnL++edz5MiREYsSz0bt7e3U1tb2iWhl+HRMR4aOq/fG8zG1LIuTJ09SU1MzYu8RiUQ4dOgQPR60a7YsC5/Pl3Zbf1mFoX6fDteYmA3RX/qlvLx83P2jzgdlZWU6rh7TMR0ZOq7eG6/HdDROLCORCJFcGjLlsTNS4Dhx4kQCgQDHjh1Lu/3YsWNUV1efiV0SEREZc0br+/SMBAuhUIg5c+bw2muvJW9LJBK89tpr1NXVnYldEhERGXNG6/v0jA1DrFixgiVLlnDFFVfw0Y9+lCeeeIKOjg4+97nPnfa54XCYVatW9TuGI8Oj4+o9HdORoePqPR3TsSuX79PBOmNTJwG+973v8fjjj9PS0sLs2bN58sknmTt37pnaHRERkTFppL9Pz2iwICIiIvlPq06KiIjIgBQsiIiIyIAULIiIiMiAFCyIiIjIgMZksKClrQdv+/bt3HTTTdTU1ODz+di0aVPa/ZZl8fWvf51zzz2XwsJC5s+fzx/+8Ie0x7z//vvcfvvtlJWVUVFRwV133UU0Gh3FT5FfVq9ezZVXXklpaSlVVVXcfPPN7N+/P+0xXV1dLF26lAkTJlBSUsLixYv7NE05cuQIN954I0VFRVRVVfEP//APxGKx0fwoeWXNmjXMmjUr2UGwrq6Ol19+OXm/jmnuHn30UXw+H/fdd1/yNh1XGYwxFyxoaeuh6ejo4PLLL+epp57Kev9jjz3Gk08+ydNPP83u3bspLi5mwYIFdHV1JR9z++2385vf/IYtW7bw0ksvsX37du65557R+gh5Z9u2bSxdupRdu3axZcsWent7uf766+no6Eg+Zvny5bz44os8//zzbNu2jXfeeYdbbrkleX88HufGG2+kp6eHN998k2effZZ169bx9a9//Ux8pLwwZcoUHn30Ufbu3cuePXu47rrr+OQnP8lvfvMbQMc0V2+99RY/+MEPmDVrVtrtOq4yKJ4tSTVKPvrRj1pLly5NXo/H41ZNTY21evXqM7hXYwOQtjJZIpGwqqurrccffzx5W2trqxUOh61//dd/tSzLsn77299agPXWW28lH/Pyyy9bPp/P+tOf/jRq+57P3n33XQuwtm3bZlmWOYYFBQXW888/n3zM7373Owuwdu7caVmWZf3Hf/yH5ff7rZaWluRj1qxZY5WVlVnd3d2j+wHy2DnnnGP98Ic/1DHN0cmTJ60LL7zQ2rJli/WXf/mX1pe+9CXLsvRvVQZvTGUWenp62Lt3L/Pnz0/e5vf7mT9/Pjt37jyDezY2HTp0iJaWlrTjWV5ezty5c5PHc+fOnVRUVHDFFVckHzN//nz8fj+7d+8e9X3OR21tbQBUVlYCsHfvXnp7e9OO68UXX8zUqVPTjutll12WtlLcggULaG9vT55Jn83i8TgbNmygo6ODuro6HdMcLV26lBtvvDHt+IH+rcrgjYlVJx2jtRTn2aKlpQUg6/F07mtpaaGqqirt/mAwSGVlZfIxZ7NEIsF9993HVVddxcyZMwFzzEKhEBUVFWmPzTyu2Y67c9/Z6u2336auro6uri5KSkrYuHEjl156KQ0NDTqmw7RhwwZ+9atf8dZbb/W5T/9WZbDGVLAgkm+WLl1KY2MjO3bsONO7Mi5cdNFFNDQ00NbWxgsvvMCSJUvYtm3bmd6tMaupqYkvfelLbNmyZdwunSyjY0wNQ2hpa285x2yg41ldXd2neDQWi/H++++f9cd82bJlvPTSS7zxxhtMmTIleXt1dTU9PT20tramPT7zuGY77s59Z6tQKMSMGTOYM2cOq1ev5vLLL+c73/mOjukw7d27l3fffZePfOQjBINBgsEg27Zt48knnyQYDDJ58mQdVxmUMRUsaGlrb02fPp3q6uq049ne3s7u3buTx7Ouro7W1lb27t2bfMzrr79OIpE4axf9siyLZcuWsXHjRl5//XWmT5+edv+cOXMoKChIO6779+/nyJEjacf17bffTgvEtmzZQllZGZdeeunofJAxIJFI0N3drWM6TPPmzePtt9+moaEhuV1xxRXcfvvtyZ91XGVQznSF5VBt2LDBCofD1rp166zf/va31j333GNVVFSkVepKysmTJ61f//rX1q9//WsLsL71rW9Zv/71r60//vGPlmVZ1qOPPmpVVFRYP/3pT619+/ZZn/zkJ63p06dbnZ2dyddYuHCh9eEPf9javXu3tWPHDuvCCy+0brvttjP1kc64L37xi1Z5ebm1detW6+jRo8nt1KlTycd84QtfsKZOnWq9/vrr1p49e6y6ujqrrq4ueX8sFrNmzpxpXX/99VZDQ4O1efNma9KkSdbKlSvPxEfKC1/96letbdu2WYcOHbL27dtnffWrX7V8Pp/1n//5n5Zl6Zh6xT0bwrJ0XGVwxlywYFmW9d3vfteaOnWqFQqFrI9+9KPWrl27zvQu5a033njDAvpsS5YssSzLTJ/82te+Zk2ePNkKh8PWvHnzrP3796e9xokTJ6zbbrvNKikpscrKyqzPfe5z1smTJ8/Ap8kP2Y4nYK1duzb5mM7OTutv//ZvrXPOOccqKiqyPvWpT1lHjx5Ne53Dhw9bixYtsgoLC62JEydaX/7yl63e3t5R/jT5484777TOP/98KxQKWZMmTbLmzZuXDBQsS8fUK5nBgo6rDIaWqBYREZEBjamaBRERERl9ChZERERkQAoWREREZEAKFkRERGRAChZERERkQAoWREREZEAKFkRERGRAChZERERkQAoWREREZEAKFkRERGRAChZERERkQP9/StLs7mZQ/zwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(examples[0, 0, 4, ...], vmin=0, vmax=35, cmap=plt.cm.nipy_spectral_r)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "498e80b6-3b73-4bd9-b3f9-2dd4d1fa12cf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}