Skip to content
ML_train_and_predict.ipynb 115 KiB
Newer Older
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train ML model to correct predictions of week 3-4 & 5-6\n",
    "\n",
    "This notebook create a Machine Learning `ML_model` to predict weeks 3-4 & 5-6 based on `S2S` weeks 3-4 & 5-6 forecasts and is compared to `CPC` observations for the [`s2s-ai-challenge`](https://s2s-ai-challenge.github.io/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synopsis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method: `ML-based mean bias reduction`\n",
    "\n",
    "- calculate the ML-based bias from 2000-2019 deterministic ensemble mean forecast\n",
    "- remove that the ML-based bias from 2020 forecast deterministic ensemble mean forecast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data used\n",
    "\n",
    "type: renku datasets\n",
    "\n",
    "Training-input for Machine Learning model:\n",
    "- hindcasts of models:\n",
    "    - ECMWF: `ecmwf_hindcast-input_2000-2019_biweekly_deterministic.zarr`\n",
    "\n",
    "Forecast-input for Machine Learning model:\n",
    "- real-time 2020 forecasts of models:\n",
    "    - ECMWF: `ecmwf_forecast-input_2020_biweekly_deterministic.zarr`\n",
    "\n",
    "Compare Machine Learning model forecast against against ground truth:\n",
    "- `CPC` observations:\n",
    "    - `hindcast-like-observations_biweekly_deterministic.zarr`\n",
    "    - `forecast-like-observations_2020_biweekly_deterministic.zarr`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources used\n",
Aaron Spring's avatar
Aaron Spring committed
    "for training, details in reproducibility\n",
    "\n",
    "- platform: MPI-M supercompute 1 Node\n",
    "- memory: 64 GB\n",
    "- processors: 36 CPU\n",
    "- storage required: 10 GB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Safeguards\n",
    "\n",
    "All points have to be [x] checked. If not, your submission is invalid.\n",
    "\n",
    "Changes to the code after submissions are not possible, as the `commit` before the `tag` will be reviewed.\n",
    "(Only in exceptions and if previous effort in reproducibility can be found, it may be allowed to improve readability and reproducibility after November 1st 2021.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Safeguards to prevent [overfitting](https://en.wikipedia.org/wiki/Overfitting?wprov=sfti1) \n",
    "\n",
    "If the organizers suspect overfitting, your contribution can be disqualified.\n",
    "\n",
    "  - [x] We didnt use 2020 observations in training (explicit overfitting and cheating)\n",
    "  - [x] We didnt repeatedly verify my model on 2020 observations and incrementally improved my RPSS (implicit overfitting)\n",
    "  - [x] We provide RPSS scores for the training period with script `print_RPS_per_year`, see in section 6.3 `predict`.\n",
    "  - [x] We tried our best to prevent [data leakage](https://en.wikipedia.org/wiki/Leakage_(machine_learning)?wprov=sfti1).\n",
    "  - [x] We honor the `train-validate-test` [split principle](https://en.wikipedia.org/wiki/Training,_validation,_and_test_sets). This means that the hindcast data is split into `train` and `validate`, whereas `test` is withheld.\n",
    "  - [x] We did use `test` explicitly in training or implicitly in incrementally adjusting parameters.\n",
    "  - [x] We considered [cross-validation](https://en.wikipedia.org/wiki/Cross-validation_(statistics))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Safeguards for Reproducibility\n",
    "Notebook/code must be independently reproducible from scratch by the organizers (after the competition), if not possible: no prize\n",
    "  - [x] All training data is publicly available (no pre-trained private neural networks, as they are not reproducible for us)\n",
    "  - [x] Code is well documented, readable and reproducible.\n",
    "  - [x] Code to reproduce training and predictions is preferred to run within a day on the described architecture. If the training takes longer than a day, please justify why this is needed. Please do not submit training piplelines, which take weeks to train."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Todos to improve template\n",
    "\n",
    "This is just a demo.\n",
    "\n",
    "- [ ] use multiple predictor variables and two predicted variables\n",
    "- [ ] for both `lead_time`s in one go\n",
    "- [ ] consider seasonality, for now all `forecast_time` months are mixed\n",
    "- [ ] make probabilistic predictions with `category` dim, for now works deterministic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.keras.layers import Input, Dense, Flatten\n",
    "from tensorflow.keras.models import Sequential\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import xarray as xr\n",
    "xr.set_options(display_style='text')\n",
    "import numpy as np\n",
    "\n",
    "from dask.utils import format_bytes\n",
    "import xskillscore as xs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get training data\n",
    "\n",
    "preprocessing of input data may be done in separate notebook/script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hindcast\n",
    "\n",
    "get weekly initialized hindcasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "v='t2m'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocessed as renku dataset\n",
    "!renku storage pull ../data/ecmwf_hindcast-input_2000-2019_biweekly_deterministic.zarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "hind_2000_2019 = xr.open_zarr(\"../data/ecmwf_hindcast-input_2000-2019_biweekly_deterministic.zarr\", consolidated=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocessed as renku dataset\n",
    "!renku storage pull ../data/ecmwf_forecast-input_2020_biweekly_deterministic.zarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "fct_2020 = xr.open_zarr(\"../data/ecmwf_forecast-input_2020_biweekly_deterministic.zarr\", consolidated=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Observations\n",
    "corresponding to hindcasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocessed as renku dataset\n",
    "!renku storage pull ../data/hindcast-like-observations_2000-2019_biweekly_deterministic.zarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_2000_2019 = xr.open_zarr(\"../data/hindcast-like-observations_2000-2019_biweekly_deterministic.zarr\", consolidated=True)#[v]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocessed as renku dataset\n",
    "!renku storage pull ../data/forecast-like-observations_2020_biweekly_deterministic.zarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_2020 = xr.open_zarr(\"../data/forecast-like-observations_2020_biweekly_deterministic.zarr\", consolidated=True)#[v]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ML model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "based on [Weatherbench](https://github.com/pangeo-data/WeatherBench/blob/master/quickstart.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run once only and dont commit\n",
    "!git clone https://github.com/pangeo-data/WeatherBench/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(1, 'WeatherBench')\n",
    "from WeatherBench.src.train_nn import DataGenerator, PeriodicConv2D, create_predictions\n",
    "import tensorflow.keras as keras"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs=32\n",
    "\n",
    "import numpy as np\n",
    "class DataGenerator(keras.utils.Sequence):\n",
    "    def __init__(self, fct, verif, lead_time, batch_size=bs, shuffle=True, load=True,\n",
    "                 mean=None, std=None):\n",
    "        \"\"\"\n",
    "        Data generator for WeatherBench data.\n",
    "        Template from https://stanford.edu/~shervine/blog/keras-how-to-generate-data-on-the-fly\n",
    "\n",
    "        Args:\n",
    "            fct: forecasts from S2S models: xr.DataArray (xr.Dataset doesnt work properly)\n",
    "            verif: observations with same dimensionality (xr.Dataset doesnt work properly)\n",
    "            lead_time: Lead_time as in model\n",
    "            batch_size: Batch size\n",
    "            shuffle: bool. If True, data is shuffled.\n",
    "            load: bool. If True, datadet is loaded into RAM.\n",
    "            mean: If None, compute mean from data.\n",
    "            std: If None, compute standard deviation from data.\n",
    "            \n",
    "        Todo:\n",
    "        - use number in a better way, now uses only ensemble mean forecast\n",
    "        - dont use .sel(lead_time=lead_time) to train over all lead_time at once\n",
    "        - be sensitive with forecast_time, pool a few around the weekofyear given\n",
    "        - use more variables as predictors\n",
    "        - predict more variables\n",
    "        \"\"\"\n",
    "\n",
    "        if isinstance(fct, xr.Dataset):\n",
    "            print('convert fct to array')\n",
    "            fct = fct.to_array().transpose(...,'variable')\n",
    "            self.fct_dataset=True\n",
    "        else:\n",
    "            self.fct_dataset=False\n",
    "            \n",
    "        if isinstance(verif, xr.Dataset):\n",
    "            print('convert verif to array')\n",
    "            verif = verif.to_array().transpose(...,'variable')\n",
    "            self.verif_dataset=True\n",
    "        else:\n",
    "            self.verif_dataset=False\n",
    "        \n",
    "        #self.fct = fct\n",
    "        self.batch_size = batch_size\n",
    "        self.shuffle = shuffle\n",
    "        self.lead_time = lead_time\n",
    "\n",
    "        self.fct_data = fct.transpose('forecast_time', ...).sel(lead_time=lead_time)\n",
    "        self.fct_mean = self.fct_data.mean('forecast_time').compute() if mean is None else mean\n",
    "        self.fct_std = self.fct_data.std('forecast_time').compute() if std is None else std\n",
    "        \n",
    "        self.verif_data = verif.transpose('forecast_time', ...).sel(lead_time=lead_time)\n",
    "        self.verif_mean = self.verif_data.mean('forecast_time').compute() if mean is None else mean\n",
    "        self.verif_std = self.verif_data.std('forecast_time').compute() if std is None else std\n",
    "\n",
    "        # Normalize\n",
    "        self.fct_data = (self.fct_data - self.fct_mean) / self.fct_std\n",
    "        self.verif_data = (self.verif_data - self.verif_mean) / self.verif_std\n",
    "        \n",
    "        self.n_samples = self.fct_data.forecast_time.size\n",
    "        self.forecast_time = self.fct_data.forecast_time\n",
    "\n",
    "        self.on_epoch_end()\n",
    "\n",
    "        # For some weird reason calling .load() earlier messes up the mean and std computations\n",
    "        if load:\n",
    "            # print('Loading data into RAM')\n",
    "            self.fct_data.load()\n",
    "\n",
    "    def __len__(self):\n",
    "        'Denotes the number of batches per epoch'\n",
    "        return int(np.ceil(self.n_samples / self.batch_size))\n",
    "\n",
    "    def __getitem__(self, i):\n",
    "        'Generate one batch of data'\n",
    "        idxs = self.idxs[i * self.batch_size:(i + 1) * self.batch_size]\n",
    "        # got all nan if nans not masked\n",
    "        X = self.fct_data.isel(forecast_time=idxs).fillna(0.).values\n",
    "        y = self.verif_data.isel(forecast_time=idxs).fillna(0.).values\n",
    "        return X, y\n",
    "\n",
    "    def on_epoch_end(self):\n",
    "        'Updates indexes after each epoch'\n",
    "        self.idxs = np.arange(self.n_samples)\n",
    "        if self.shuffle == True:\n",
    "            np.random.shuffle(self.idxs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.DataArray &#x27;lead_time&#x27; ()&gt;\n",
       "array(1209600000000000, dtype=&#x27;timedelta64[ns]&#x27;)\n",
       "Coordinates:\n",
       "    lead_time  timedelta64[ns] 14 days\n",
       "Attributes:\n",
       "    comment:  lead_time describes bi-weekly aggregates. The pd.Timedelta corr...</pre>"
      ],
      "text/plain": [
       "<xarray.DataArray 'lead_time' ()>\n",
       "array(1209600000000000, dtype='timedelta64[ns]')\n",
       "Coordinates:\n",
       "    lead_time  timedelta64[ns] 14 days\n",
       "Attributes:\n",
       "    comment:  lead_time describes bi-weekly aggregates. The pd.Timedelta corr..."
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 2 bi-weekly `lead_time`: week 3-4\n",
    "lead = hind_2000_2019.isel(lead_time=0).lead_time\n",
    "\n",
    "lead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mask, needed?\n",
    "hind_2000_2019 = hind_2000_2019.where(obs_2000_2019.isel(forecast_time=0, lead_time=0,drop=True).notnull())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## data prep: train, valid, test\n",
    "\n",
    "[Use the hindcast period to split train and valid.](https://en.wikipedia.org/wiki/Training,_validation,_and_test_sets) Do not use the 2020 data for testing!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# time is the forecast_time\n",
    "time_train_start,time_train_end='2000','2017' # train\n",
    "time_valid_start,time_valid_end='2018','2019' # valid\n",
    "time_test = '2020'                            # test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n",
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n",
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n",
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n",
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n"
     ]
    }
   ],
   "source": [
    "dg_train = DataGenerator(\n",
    "    hind_2000_2019.mean('realization').sel(forecast_time=slice(time_train_start,time_train_end))[v],\n",
    "    obs_2000_2019.sel(forecast_time=slice(time_train_start,time_train_end))[v],\n",
    "    lead_time=lead, batch_size=bs, load=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n",
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n",
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n",
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n",
      "/work/mh0727/m300524/conda-envs/s2s-ai/lib/python3.7/site-packages/dask/array/numpy_compat.py:40: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.divide(x1, x2, out)\n"
     ]
    }
   ],
   "source": [
    "dg_valid = DataGenerator(\n",
    "    hind_2000_2019.mean('realization').sel(forecast_time=slice(time_valid_start,time_valid_end))[v],\n",
    "    obs_2000_2019.sel(forecast_time=slice(time_valid_start,time_valid_end))[v],\n",
    "    lead_time=lead, batch_size=bs, shuffle=False, load=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# do not use, delete?\n",
    "dg_test = DataGenerator(\n",
    "    fct_2020.mean('realization').sel(forecast_time=time_test)[v],\n",
    "    obs_2020.sel(forecast_time=time_test)[v],\n",
    "    lead_time=lead, batch_size=bs, load=True, mean=dg_train.fct_mean, std=dg_train.fct_std, shuffle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((32, 121, 240), (32, 121, 240))"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X, y = dg_valid[0]\n",
    "X.shape, y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x2b55b159ba90>"
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAELCAYAAAA1AlaNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACgJ0lEQVR4nOydd5gkZbX/P2/Fru6e6Qk7szs7myNRckYBBQFRETFwvQaucBUvXvNVL8rP7MXINYKYMFzFiKCIZERAcl7C5hxmdnLHSu/vj7eqpmd2dmZ2d3Z3dre/zzPPdKiufqu66vue95zvOUdIKamhhhpqqOHAgLa3B1BDDTXUUMOeQ430a6ihhhoOINRIv4YaaqjhAEKN9GuooYYaDiDUSL+GGmqo4QBCjfRrqKGGGg4gTDrSF0KcI4R4SQixXAjxyb09nhpqqKGG/QliMun0hRA6sBQ4C1gPPAr8i5Ty+b06sBpqqKGG/QSTzdI/HlgupVwppXSBG4Dz9/KYaqihhhr2Gxh7ewDD0A6sq3q+HjhhtA9MmTJFzp41a7cOqoYaatg/8MSTT26VUrZM5D4PE/VyERn+KDeJidzv7sJkI/2RTto2/ichxHuA9wDMnDmTBx54YHePq4YaatgP4KTTayZyf28QbbKDCv14E7nb3YrJ5t5ZD8ysej4D2Dh8IynldVLKY6WUx7ZMmbLHBldDDTXUUI0n6OU0mjHQeINomzwB0lEw2Uj/UWChEGKuEMICLgJu3stjqqGGGmrYBm8QbdJAoxmLY8jxBL17e0jjwqRy70gpfSHE+4HbAB34iZRyyV4eVg017Jc45lN3Ylg6AJou0A0NXdcwLB0hBEEQEvohQRACIENJEIQ8dOUZ4/6OQz/8VzRdQzc0DFNHN7Tk+wxLx7INDFMj5ZhkUwZpS8cydHRtqKc3CJUR/b3zD5qIQ58QPEEvJ9AIQDNWYu3/aZL79ieVZHNncMzRR8vJ7tM/67sPUS56WLZBtt4GQNcErh/S21ng/k++YsK+K3fqB+i7/9u0nX8Vuu2w/rcfnLB917Bv4Kgr7sCr+AAEfkjouxiWjWHqpDIWhqmINwhCfHeQ0A1LR9MFYSCplDzKhTIyDNAMi1TGAkDXBz8rQ4lpG+iGQIihPKcbGkIT6LqWbOt7AYEvCX31nWHEPZoQCE1gWNHEEE08AKatY9sGjqXTkLZwLB09+i5dE9jRJGIZGp8+fc6Y58ZJpx+XUh670yc3whtEm3yGfl7NYEy4C5eH6WGtLE1q0p9Ulv7uwodvWYqla9iGRsUPWdNVoHNrkVTaJJ1SpyBt6XTlXf5yyTET/v13vP/ECd/n9tB3/7epP/ly0s3TAZh50XfRDXXDrv7le/bYOGrYeygNlDBt9ZunMhb1jTl8LyDfW6Zr3SZkGCA0HbuugWyDAyhy9b2AMJAUByo8+/Vzhuxz0ft+D0A6V08qY2LrBkEQEvgS31WkHiOUUq0KInKvhozeiyE0NWFomsCrqOeGqWNYGqZt4FZ8PDvATRnomsAyNHRd0Jy1+O9XzJ7wczdeVFv5MfYVa/+AIP2rz1sEwGU3voCuCXRNcNcHT+K1P34cK7IUfvbWw/fmECcU/Q9+b8TXW8/9LB23fpa2868CYNNNtYTnyY4jP3k7gR9uQ8IxFr7nN4S+y4qfvIM5b7+OwHfRDYsXR5jgZ73tGjRNR7ccDCeLaRuEEQEPdJcoDQyw9Jo3JftNNzRgmDpBEJLO1Sf7qZT8QWs9+nwYhIgql4wMJWEoCX236rUAGQTJc6HraIaF0HQ0TVn7OmplIHyBrqvvCPQQKSVBKCm5PnrK3KuEX+3LH45jyPEwPXthVOPHAefe+c4jG/jP49t344h2HNMv+Bq67bDuhvdP6H7rT758mwmg9dzPYqQyCE1PbkjNsJBhkDzf/OdPTeg4ahgbR3z8tsR3DvDcN87lsI/eynPfOBeAgz94Cy9867wRPzvn7dcBaiU34y3fIvRcZKjIVWh68t/K5DBSWUW+YYBuO8k+YrIfjrbzr6J+xmLMVApQ7qLALaEZFrqhJo3ALSXfAxBUSoS+i9D05C8eT/w/RrJNNAHohoEWuXcMU+0zlTFpasnQ1uDgWDp1tsHnz5w3zjM7FBPh3pklHHkCjSOSPsDtdPIy6pms1v4BR/p7C81nXYEMAgK3RN/9356w/daffDmAst5sB91KJTeZDEOEpg258Q1H3fSh56LbDlY6h247hJ6LV+wDJsYN1HL2lQAEnjvE2uu972rqT74cLXI59d539S5/1/6IBZf+KiFOzbDQTCuZqEci6Flvuyb5XUPfTVZxred+FjOtrHQZBliZHGY6h+FkCX2Xl753wS6Pte38q6ibvoDQU79z4JYIfHcI8ccuJcNyhkwAsaUvNB3ddjAsOwn8xvEDoQlsxyDbkOLIec0cMTOHqWm8+eDmHR7rrpK+EOKkuaQfrPblD8dk9+3XSH83IyZl2L7bZUfQfNYVyc1VbTXptpP47jXDRI9uruptNMNC6IM3YWwRxqQSf756khCazoqfvCPZx/x3/2LI8/Fg2uu+RGWgh557v75zB72fYvoFX8PK5NT5NwatRs20kEHAsuveOmT7he/5DUAyica/w6y3XYNuWOiWw/IfvY1Zb7uGoFICwCv203HrZwE1EXulPL33XU3ruZ9FhgFGKpPsf1fcfe0XfhOncRoAbrGPoFJKrr8h11N0/QHosWvHsDAzOYSmE7gldMuhvinNlOl10ecE2ZRBW4PDiXObeOshO5+bMwGk/53X0vr+dpxRt7uZzWyi0iKl3Lqz37W7cED49PcWml75CTTDIvTdcRH+SBNE7KKpfi9GbCUBahWBi6bphD7IULkKqi396glABfKaktfDqveGP499wTE5zX3njwHwynnym1ePaq23nX8VZiaHVdc0hIwCt3TAu5ECt4wXuV2W/+htLLj0VwD4pTwA8y6+Hhj8nYUekWeob7MvGQaU+zsB6F+/dMhvUr2yCn2Xqed9ASOVwSv245cLABipDO0XfhMZBmy88b92+FgGNq1Ixpeqb6FYyuOXCwS+uiYVuZuJgRFDtxzMTA7D1LEdE0gRBCFO1sKI4m2+H+JGMYRdIfwJgmGOI70p2mZS8usBb+lPPe8LwCAhBpEVpRtWYiHtKnKnfgAZBqMS/3BSryb9kRATsG5ayXOrrhErnUvcAL5bwrCcITd89X8AI5VNXD5eoW+b9wPfTaw1M5VNJjEYXBEAVPLdaIaFYTmEvpt8Z3XsoHoy0aqsv+rjiYkjqJR2eEUxUWh4xYeHWN5h5KqYiJXaeNB81hUYloPQNKy6puQ3He4Pl2HAip+8g/YLv8mGP3wEUKsHzbSSyTW23utPvjyxquPfUzOtZNWomRZmZPWHvodupQh9D98t0XXHl8c17vhazU6dQ+i7uIV+/HI++Y2tyJrXLRUfiH9vK50bMhnolkOmPkVDixqPYWnUZSzSlr7LgosJsPSvebNou2yasEfd7uZwC2sotUkpN+/sd+0uHFCkP+11X1JEaTsElcHgU+zaiJ8HVSSllAXqda+UJ4xIUIYhW25RfuvYGm88/WPAIGnqlpOQchi5UwK3NK6xjkX68Q0MYNgOqcapGJaDW+gbErCrthC9Qh9eOR8FzKwhJBIOI5QY1RZ+fOwAmZbpZOpTGKZGpeRT6FeabtO2cEsl6pqUH9l2DPq71TFXCnmCSim5wauDjclYNR3dMBIpn5Qy8e0apj5EIeJ7Ab6n9qHrGlpkGZYLLn70W5X7OhnYtGJMwp7xlm9hpXPJmMxMTo0/k8WruImrwi/ltzuhGU5WPQ6CZPvhqytQ10dQNbHGbrV42zAMkmsuDIPk/cpAd7I60y0n2Zfvlob480Pf28aIGX7tCU3fJrZUHfjf3nUHY7spG0//GE7jVPJbVo/4fnxugcTQSDdP3yMrv4kg/Yu0sUn/T8EWVtdIf/dgvKTfcvaVZKfOAZRbQgZBEhwDdcPphpXccNXEX+lTbrnO29SqICZ3w3aipXcccIpvam+ItRv6Ln6kaBhure0Kqi3t2IqK/bSaYWE62SHbVJOOFvmAYx9rUCklpFBNwGKYRR76Lk5jKwBO1kY3NHw3wI2SgeIknjgBKJWxsGydvq4iXjlIEoDKRQ+35OG7lWT/pm0l0j2hCXw3wLINLEetkqWUBL66XnVDRGNVipKY/AG8ik+5r4tC51qEpuOXC1QGupPzU01auVM/kBCp0HSsukbMVDYh1dj6HK5AiV0u1asiy3HQNKGkikFIGEq0qkkqljDGVrhfzuOV8uS3rI5+NzVppKdMVyum6FqMf6tqt0gMGQaJ/zyeANRvWU62qQyoVVg16e/IqnMkjPb5xtM/hmE7VAa6R92HXddEffsiulc+Pa54z9TzvsCWW66k5ewr0a0URtX1PV7l20SQ/tv06Ze1jUH6fww2s0rWSH+3YDykHysYrCoLA8At9CXLYhgMiA0nyvhmGb4CqLZQh1t18Y2otldZkWF0g8JQF8pEICauMNJC66aFbqXUaiOS5sXjiV0+w10Gw8dUfR6AbfyxmcYctmMkxKsbWqLbjkk/Tr8XmkiyMoUmyPeWt0nxV+OXiVQvJkzTNqJxR9mg/iCZuhUf31XH4JbUuS12bdiGBBpe8eHtumhmXvRd7GxT4vKodu3E44hXGHHCkWGqEgKD50ZLjqP6mHRdI6xarcTvV//33Upipattgm3+69FvVm2oxL9fUCklyqv4Gov35ZcLicERJ+7FqD4X4yH64Z8ZDW3nX0Whc92I78Xj3xkVW9v5VxH6HpWB7kQdtCP7mgjSf5cx/bLpWmrU7X7rbWLFJCX9SRlo2FW0nH1lIlPTI9dHbDHFF8nyH70NGNQ4x3KyarcJMMQCM6Ln1W6JkTTIw6380HPRrRSBW058qHrk+95Ry3/4RBPDj9xVummhpzKRXFMdaxwYjBU+Xlk9j1c78TbxmKuJffixxioRUJmfvqusc83QklothqmTakhFn1cZmJWSRxjIxB1k2QbYUBxQlr5h6xT7+jFTaaSU2CkTzdDQNJFY+qEfRqn9kkrJx/eCwckiCJNzm2nZtr/C8MAmDBLYuhvez+LLb0yOP56ogCGknZuSVjr1quOMUSl7yBCEIRC+muCkkLgVHyklWsqMzqPANAy8io8IBbqpoekOoWXjFgaS/VVb5YHv4pXy2HVNhL47RF/vl5TP3EhllfJl2HWhGR490Qo1xvaIe/jr450IRsJohL8rkuXJkFCoCdDHEGMKwQhF4ScH9jvSb7/wm9i5KckNHAfBqqWK1Trn2JetV7kzQt9NtMbp5vbEOq9OaNGG+WGrH8dZh3EAMDAiaaRpDbGmw0iaF283HgyXYCb+ZV0RcqyQiL+/OnEmqJQIPXdQBaLrQ2Sb20xenruNVelH5yXev1+K5J6WE6XPD10dCE0krhdNF/heiFfx0XUtcQmBcsnIQK0C4oSkY6+8S2VhDqjvMywNAvDdwezPcqG4w1rz4eS26H2/V+cy2mc84YCy8pumZZFS4lUCdYymTqXsYdpGYtXbKVPVlglCDFMjXkFrhpasegDCQOK7AbZjUil5Q0oSxNdWUClR6tlMqWcLqdwUrLom9Vv5rjIWvMFAeuyik0GQyHT1MNiuu2M0wt+emGBnsL3YwETmqOwt6EIk9X+2BzFia5DJgf2O9IWu4zRMG2INDQ+qLbj0V+i2w0vfuwAzsuIDVy2BV/38EuZdfD1GKksq8hka1lD/XUw46oZVWYnDA3YyDBDh0NWFH2UyJjEA3yPARded7Vrw1b70aj97TNzVKhjdSiXHHXruNn7n4RNU9T6BIZ/Vo3hFnEwT72e4n79ayeOFOmFgqJopJdVUIl1vAwZCC5LU/Xgi0DSRkKZpG7gFEsIHeOwLr0oyVQ1T1YYxTJ3HvvCqMa6C8aPt/KvYdNMnmf/uX6DnmpkyvR7TNujarKxuTQh6OtTKKPAltmNQ6K/w5JfPAlQVyfiYsg0phBBYtoHvBbiVAMPUkuJjQOISit1gvhsQSok6jeo93aijfsZislPnqu91S2ANXiPVmdTxRBBDM6wJUz2NlNE90mujYTyB4X0NgslXk35HsN/59NvOv4rcrEMS63zZdW9NdOWglAPVlv6Mt3wLGHRfGJaTBM1i60/Th/7EcbCuunjUSHVGgCEp8YFbwi30DbGoZRgOeT5cRTO8Vkn8XDe3jT3Ek0lMAoFbSvT6oJK2qo+1msyBRH0SJ9kITccv5ZPxxS6patWTDNTqJ14t6IaGEGKwKqMhsFMmf/+vU9lRvPyr/0iqRe5IOd9dxeLLb8SpU3LBuiaHdNaiv7vEPz99+rj3ccLn7qFccKlrchCaQBODMYjhAenAl3gVP4kXxNdVXOIg/s1jiWX1JB6f9xjxRLCjxD+WWmciVwF7Es1nXTFEcjoRPv3/sGdcNlMf3ad/fXkjL4XFSenT3+9IP0ac2LLy+osBlUkap3pXuwPaL/wmoAg1DnDGss5qEhuO6hsTGKLKqLa+RvL3x+6R6qxFZcF52yV+oQ9KR0GVN9BHUHPAoOQvJveY7JP9RYHBeDyGk8XONg1ZHaUyaQAsxyT0Q0IpKRdcdENLAquaEIkaB6BcVO6K6uqK2ysUNhpO+9r9quiWpqzmuz540g7vY1dx4hfuAdRkc+IX7tmpSeeIj99GKmORrrcT48Cr+OrYhCphHIYymdhi+F4w5JqTUiKEWh1Vx4HigG016Q8Pvlfy3eNWtuyoNT7ZiX/6BV8j9F3MTC45BxNB+u9Pzbxs1hik/5PyBl4MJifp73funRhWpGZZ+J7fqJoeTnYb3++st12DGVm3oAJmQWTJmpEFFQ6zqGMkZB1UWenV/np/aNGr4ZI5UBNNbL0PWt0e+C5B7LfV9SHkXq3hhsHYglvow7CVCyCslPDcEmYmp2SoEaeohBs3eT0O/vWtX0quff6QWieBH2JYekJIpm2QbXDQDZEEMIUmKBeUG0fTFUFXSh6l3u4kUL6jOOZTd6IZGqmMuccJ/9AP/xXTNnjqqlcPIfmdXWU8/dWzOfiDtyQWPygJq6pV76lyxm6QuK6qJ0vfU3LV+LeASM1UlRxnONnBmEy0CqguewAqpjX3nT/GK+cJPXensm33ZaRyLYkQY6KgjyeQO6HfOLHYLyz9zYsvY8VP3sH8d/8CM5MbsuTVq6R0MQLfV6QcVLtZBok3DmpW69ZjVOv6w6rPDJd6JtZ7pP8fTvrx5DJc6umX8gRuOSH95HujyUHlBphDYgVWJpe4Xgwni53JokV+c8PUqZSUQibbkKFnw8YoW7KPQuc6GmYdnByXH43PsBzsXAu6YSTnKx5jnHQFqvqhaRuUCy5aVAfdrShlThykjScOIcS43SOnfe3+nXIHTVbEfv8lV7+GIz95O4apk6m3I19+gO+FQxqReJXBpDOv4icxo+HKsdj96JWLALgD3ZhRgtnw0hvVRkkQrUo3/OEjE+Zrn0xW/3GfuRvfC3jyy2cx623XUN82J3lvxTUX7rKl/9HMzMvm6KPX3rmmuIHn/ULN0t9dWPGTd7D48huTyoFBpTQYlPS3lVUmgciq4k/Vipdq8h6+fK5+P5bHhZGqZXjBsnhSiEsZhFU3bfVKQeh6MmY9iinEE0Q8AfnlApo16Duvn7EYoenkpig3TCptUeivUC64WLZBKCWWoZFKm4llXi66tM6biVvyKeYrKjM5qr8ejwnUqsF3S1jpXFJGN9aIF/rLyURaLrqqc5KuDWnCEUqZbAPKUj7hc/eM+/fcnwgfFNnHeOqqV3PcZ+5ONPq2Y2LZMpksAbINKcpFtRpQVn4ImoGmK3djvPry3QohYKbS+G4FM51LVo9DjJMqtQ+o4L80Ldov/CapXEuy3c4kEG6vLtTexKOfe2XyWDcsfNeb0P3rjEe9M3mxX5D+4stvBAabNFQrTKotI1Ak6kYZikLfVmOfKCSigJkZVUEERsxYVSqdod+ZELamE0bb+5FeP8YQ5UsYSUUjiaTQ9SRVPf5MffsiQt/FK+cxLAfTVi3snKxSFg10lwgClSClGRp+yaOuyaHQX0ncNqm0he8NuhAMy8Gqa8KwzCErocpA72DNdNuh1LNZrWRcHa/QN7iKsh0C34kKZSl3j1fxhxBYDLfic8yn7uTxL525Yz/ufohHP/dKTvny3wF1znRDw3KM5JxVSkrSatk6YKnMZS+gUshT19SAiM635ZhJwlu1fNcv5beJKRU61iauzGoFWFxlM3DLGAzKiKtjUMOzmHdUwbM3sernlwCDFUonAuPS6U/Yt0089gvSj7MSq+uPFzrXAlEhsnRuiOtDb5yStGiLpXR9G9ehWw5OQxNhEFIu5RNCHAyqOYnlHU8M1fr2GHokr4PB2ijVKphqRY0aY1Xyk+0M8dnHyVBmJodX6CM7dS5C03FLJTRN0B/5e9M5W1nwAxVMe1AWmc5aVEpRv9QgJKyE0es2pS6XSl8nFaBp5uxEPhmGuUSZk9+8CruuCRkG2LmWpKYNKHLxCn0YZjvlQjk5/5bjJJNI7I/uXb82uQFrgAeuOI1Tvvx3ZCgxbJ2UYybNv23HwLTjHrFGIhM17fokPwBUVrCmiUQiO5gNnUvyBHxXFdKLSz2AKoim207U30GVa0jclJqOsAYNkjjEHJce6bn360MIPy5OF5dGgMnVhCcm+1jiPBGokf4kQGyJaoZF79oXaJxz2GAxssgl4pXzKnW7cVrSrCEMlFULqoBYsWtLksofV54MfD1J2ok7Fx38wVsSN0ucWFUdxI0LYMGg7DFwy8kSujp4DKAZGey6pqRkQrXvNvSi1YYmSLVOpVwok87ahNJC1zXyvao0rlNnkUqblAsubmTlx0XI3EpUfkHTCYWgNDBAZaCX7NSZyYQWSwkBKn2deFF1xHgcVpQNGh8TQH7Laqy6RgY2r8VIZZOJ1fcCvIra9sXvvD46xpGVRgcyHrjiNM79waPohkY2ZSSkP7s5Q1/Jo9Tg01v0cN0AXdfwvYBi3h1sVSgEhmOoXs9ZC98LKBe8qBDdoGFhZnK0HnIKA5uWUxnooZLvwYiux+EGSyz9jY0ZK5NLSoiAIvnq7ObJ3gQn7kkw/92/mLB9mkJgaaPT+hhv71XsF6Qfei6FjnWkclNI5aYkcksYtGDW3fB+Zl70XdxiH2aVPz0uOeCVy4krR4ZBQlYw2BQ6VgINIfwqxUSpZzBmIzRVr96PVDSqhk8ZJz0YCNVMKyn9AINuoWo9fjz5oBk4WYuGlgx9XUXSWYu6JicpP1voL+N7AQ0tGQxLTwKppXxlSB0cgGzDNIIgjMoh6LgVlSgVf6Z+2nQqUXKVl+5L4iSxS8qN6rz4bgkzrFdxgSiuITQdd6B7m+5be6tM8mSHZem0NThkbYN0lM1c8UMsQ8MyrOi/Tm+9zZbOIpqhYUclHYIgVKtVIRCayv61HBOv4jMQVTYNQ6km72BQEmxYTnKtVrvq4mtO+G4kXza3ERWYTpaWs69Mig/uK4ivP+eG9+7yvsaXkTt5sV+QvgwDMq0zEwKtTr6Km0XPvOi7CfHHvnG9yg862D1KwysPVc6M1J4uTt0Xmo7T2Erg+4n6xYysXquuCTtW6OSCIdLMGHEBLcMyCfyQSl8nMgxwk5Kz7aTqc1i2kWSwZhtSyY2flBY2NOUqsJQM0HcD0nW2qvUSuQq8iioTEMrBomiaoWFHssB0Ng5EC5yshVcJGIis92p1k4yqJ5qpTFL7JX6vku9m/W8/uM35ihuolHq20H33V3b0J95vMa81S2udPcTSL7kBQV7SkDaZ1pDjyTW9pC2due11DJT9pKGI64fomsCPnutGSGnAxasEOHVxtVOdSlkjDmWmclOQYZgYN8MTveLALig//75G7jHazr8Kv1xIxh+375wIjLv2ziTFfkH6ViZHEJVQiBFn2q7/7QeZedF3k6JjcZLGgkt/NaQuTbm/k0zLLGQox1XLZek1b2LR+36vLCdTR9M10s2q4XocSPMKfYn0sW/DMuraFiBDVXM+RpzkJUOJZRtMOWQ+5YJHpl4FaEt5F8NUCVHVMr44OBtb707WSsr5OnUWTp2F7w6NHdiOkWjCk6SyQCVTGaae+PR9N0jGZVh2Yi2Cqn8er46cxmmUejYPqVQ6HNMv+BoyDCZFoazJiNnNURKcodHkRHXz/YCmrIWlawRSomuCQ9tVYH/Jhj668oqosykD29DwQ0nJDSiWVG6FGQyWmlYF6jQCw8JM10eB2XBIzwe13WAviVg8sK+2t5z2ui9RP2MxuuWw+PIbKfVspmH2YQDk/7Hr+ze1sd07tdo7uxkjNfI2quqRxEQ/4y3fws42JcXUAOzIsklIrVxkwaW/GlJHPZ1TPnjfDSj390VuGRvdcrAzWQY6NpFunpro2kPTwiv0Je4QPZWiec5BUdDYTILHWpRlGRBSGegndLJ0rushlUmhG4PjDyOFhu2o1YCdU0lSmq4l8QahCcKoRk28IjAsDaGpmvOgKkACaLqG0GTiJzZtXZVGDrRkXPnecuLuCdwSS695Ewd/8JYhvXQ108LN94z622y88b9oeuUnRv8BD2AEoZLW6prAjibwtKlhGzqhlAy4yrJ3LJ0glOQci42R60aPy0v7YbJKEELlSySBdG1QlZXKteCV8/ilPJmWmUOarVSrybr2UbKPsfnPn+LQD/9VlU8xHVKZeeS3dkzY/jUY270zeTl//yD9kTB8Iph50Xe3W4xq3sXXs/L6i5nxlm+x/rcfZPHlNybVKC3HYaBD9R51Gpqob2lOZHaGqYKklj2DVGZQD79uSQeB72Klc1iOg25oONnBEsSglt1xvRVNE4SNDuWiR645TbYhlWix6xodZrVlWbW+HyFUUa9YppeL3DFHzGzAMjRyaRPXD9ncW8aPbvrF0+qSpagXSh5e2c2KDX3IEKShasD4no8MJZWycgL0bS3iDnQny//4XMaB7Bhz3n5dUsJ6e2g8/WP7rMW4J/DRk2fys2c62JqvJG6betvANnQCKUmbOmcsVobJAyu6WDg1S11K3bYdAxW68hV0TahJwzYAVeYhvs70yBAIQ4nR1EIYNOOVi0Oa5sToW790v/ut3FIpElKMr2PdeKCJwVpK28Mk5vz9g/TbL/zmYEA28oVXV4LceON/jVp/JK7PE/uiq9078YQASrWjGxqmYaDpSu6ZyphJE49YStd+yHy8io/vKt18Kq3871rkd4eq2vCGqi/vVXycrEW2IUVz1mJrtF19nYEeNSCZ0urQnLXJOSaGJpjfmmVrXmXbHtSapaPgUm8bBPU2rh+ysrPA6q0FOqKa9TMaHVrrbJYUPXRdQ0otGbtKFiIa22DNdhkEQ3qwDsfwxjQwtIbLSCUsahiKXEol08VEommCih8QSPDCkLSpU/FDjprVgCYEnf3q98yXPXRNEIQyIf5KZei+4980jvlIGTWmyeQwI2UOqGxeu65xzxzwHsCSq19D7tQPbFPK2Ul/aJf3rZkaujV6nU0xhvtnb2K/IH0ZKi2ymcmRbm7HLfYlxc/KUavDnUVM+KAmFMtuxjA1bEe5aWzHIPBDvEqQWOdxFmzse4+bgMTVOePH5YKHH4ZR+V2TctFj9fMdbE5byUXT3KYsu4VzGli9SZX77S26BKGkEgXyAHrqPA5pybK16JK1DZ5Y24sVWXvx0n9aLsXariKXnrOYn96xDN9VEkDdEKopSTQBaIalrMMoRrA9wu9e+fR2S+fuK8k7kwFvWNTEb57fShAF2ItekLgPUrqGF0gCKekreliGRiVaEeTSFiXXj1xEOgNlL4kPJbLO2O0TGQ5CCMyUiiMEvp/47wPfnVT6+olA07wjWPS+3zOwcfmQYoK7CqEJhD5WceUa6e9WVNfLSZp4e4ONPiYKuu1gVWWf6rqG74ZJgLQUBdjKRY/6RofAUSUKBoOmIWEgsaLleRycjRU38X5kKBMlzZa1fUNqya9Z1UMqo1xFvQOD2bZrugoUFrfQEKl6dE3Q0V8mCCW9RTWu3z60ltXPdzDnkFZ8T2m/44YfccVMIKn/7ruKDA7+4C3buHZgkNjbzr9K5UhMcs32ZMZbD5nC9x/bCMCMnMOmgTK6JtCkoOwHbO4tYxmaaioTBWmVekdD1yBrG5RcZXTohkYQx2e8EAwNEU0EcS5GXH/KLw8VOOwvWHDprzBSWZXDYzuJETgRELpA7MPZWfsF6eu2k9TICdySao0YuXtMJ0vzWVeo7FnTGlJbO1b4wGC5g5GqEMb1+O1ci/LLR1aT4QyWXXZSRuJjL7kBxYJLOmPh+yG5rIXrh1TcQHVbiqywTJ2F7ajiZEU/BAT53iKptEW5qPzrqbTJ1o39OFmbcsGlv6uI74WJe8iOVhFuxWfJhn4uOrqdm57dTOdAmc0dBT5w2nw+f8sLAGxe05vU6omRqVdZlF7KpxC5DXRbEUbohxiWiVsssOh9vx9Ruqq282qEPwH4j2OnJ49/+vQWnlnfx+zmdLJSGyj7lLwgee76AY5lqMlAqpWfKnpX1UYysvxjxZYQAonq/bC933Nfx2EfvZXALdG1/Aka5ijVzkRb+toYpF8L5O4BqNr3VlLcLEk00XSISs9qmk7DKz6Mlclt0yhdpJV6Yd7F1w+qbiwHu66BbKuSYlq2gWFqOFkb09aHKGeCSDYH0N1ZYFpbHSU3oC5tUnIDLEMjnTIouYMB2LSl09FfQYsahcdWfqXkJfVsAh9MW6fYX6F7Qwe67TDQnSeVSUFAIrMUmsD1Q0IJ3YUKq9b3I0PJ8535qvR8HTtl0tORV+4bqZpzx60LE1WRrmq9O3UWQgjKtkH/5u0XC9xyy8RpoGtQsA2Nq89bxIdvWUpbTk3Ma7oKNKStxG0HRuLX16Pf362o+FBcegMg31tIkg3nvP06vFJ+iB5/f0Ndo4MQ7di5FoJKCb88WAJlIqAbGro5eqxqPD59IcRPgNcCHVLKw0Z4XwDfAl4DFIGLpZRP7MyYq7FfkH7cjCSoymyN69tohomRyuCXVbkCu64JI5VJKkzG7RKrG6InbiLToti1KanN3zS1JSknHN9osU89vukA/vHxl3PSF+9l3nz1uSCUFMuqEJntmMl2+fLQG9St+Gi6RiY3tD1jKa+qWda3NBMEIYEfkt+6FRkG1LVOBVQxtbtvX8KLK7uVZt8P6esq8qXv3EUqpwJ0LTNyvOqo6Ty9rpd1q3uYMr0er+JHfn0tce9k6m1KAy6WYyQrm4FaQHaP4m2HtvCjpzZz9XmLuOL2FQShJF9W/vvmqMheLPccKPvYhkah4HLH+08cdb+xEqvl7CuZ8/brRpQ77wvYXtG3Y6+8K3ksohW/QXZC3bxoE+beuR74LvDz7bx/LrAw+jsBuCb6v0vYL0g/bhrtlwtJdmisw49r7xipTFQjR7lN9CgGUF2N04h64kJ0wWiCTGMu6Yd61BV3kMqYOJaeEL4R9z6tKiV81ncfSmrHn/nth+jtLOBVfFpm5CKSH5oYle8tE4aSVMZKmpTEloJh6jhZO6pn46P5GroeYrQqSy0Ovg70lMi1NkX7V8eYqbfJNc9K9PZTpqRZ2JpViT+HT6Ov6LG+p8Sm3hIlN0gyO9Mpg5ktGYpuwNaeEltWbtjppig17DzyZZ//uW8NX371fM774WMYhkbvQIWGdFQ6QRMMlAPyZZ880LG+f9z7rgx0JwXS9kUMJ/zms65g+hGvSMqK64bKYakMq6I7ERC6NiGBXCnlfUKIOaNscj7wc6kCMQ8JIRqEEG1Syk3jH+222C9IX5VB9vCjYI0fFQvLTp1DqWdL0ngEwKprTBqGVAd+DSeLYdlJgDbWOduOwUlfvBfLNhLyf+PPnqS1PoVTtcQLpHLvlFyfujqbc659JAnCZhtS3PvRU0Y9hkM//FdSGTXGWO45a4FaKbh+GAVy00Na64WhVBd2SRG9aeuUBlyyDSnS9TZnHKpWAR39FWZPSVOfMjliWh1ZS2fTgMtD63porbNprbMZKPsYmmAg2v+m3hKrlnfz8GfOAE5gztuvA0ZOhKth9+BDJ84A4N//sIRb/n2w78fZ33+Y3s4i6XqbUr5C4Esy9XZyfY4Hw0mz6ZWfoPvur+xTZZNjnHrVfcw85nRkKJPCgV65vE3wdqKIXzMEujmWZBOAtwshLqp6+Top5XU78FXtwLqq5+uj12qkryx2LSF3K6MShgJXVbRUXbDKQ2rgm6nsYCVOJ6uyYCOyVyWXVXajZRuJ/j6uB//sA0tZ+MZjACh5iuiLbkBfcWizhinT6/A9Rdgv/+o/CP2QB644bcRjWHL1a5h38fW0HzIfJ6t6qnZvLXLCoVNZ31MkCA1a61P0FV38UHLc3CYsQ6M7UgwNlH1e3NTPyq0dNLdmOHF+M/cs2cLC6fW01tn0FT3yZZ9/SompCf74xAZOWTiFgbJPXZTOH+v5n3h28zZdrmpkv/fwwwsPHfJcaILGqRkK/RVKAy5zDm7BdcdPaNNe9yWKXRuHvBYT/b5G+Ad/8BZgsLlRdTe6+K+66u1EQNM0le07CiKBxy+llLuS7TbScmGXWx3uF6TvFvuTqoAwNCFIPTaTWvVxCQFVK99OJJixXFJ9RmXcxfr6ujrlQ9WjMgbLrnsr33hwHY6ls2RDP71FL1FUJN8bhe9DP0xa4nWv2zDqcay8/mJO/MI9/OPjL+e8Hz4GqGBv2tLpyrts6i0xryVL2tKZ15xhWWeeT58+Z5v9vPt3z7FkQx//evJsvFBldQJkLZ3WjE0uZbC6q8ijq7q54V+P2HYg5x806jhr2LvI96oa+Km0xYyFzTRnLXqLHkd+8naeuurV22zfdv5VuMV+gkqJ3vuuHlGPH+dZxJ2w9gXyn3fx9UPIXIYBVl1T1CcgqjIalT2PhR0TATFxPv2xsB6YWfV8BrBxO9uOG3uN9IUQq4EBIAB8KeWxQogm4DfAHGA18BYp5ejFXYgbjaSSCyC26mHQnRM3Arfr6hP5o2UbpKOgqWWrLNukTo0G6YyFrokk7T0OwIJKn//ivatZODVLyfWThJl4OyllYuUDlAvuuMoLx024D5vZAEDRDSi6AUfPbuSFTf0cMr2eUErW9BT56MkzR9zHT958GPes6efw1gxLu0vMitQfs8wSUhMI6TGtIZVMGG/82ZP88V1HjTm2GiYH4rLJueY0QhPMbs5w9013bLcBeKlnC6Hvjkrk1e9Vt0CczOS/demjpJsHZa5mJoc70E0YBmhVnezirnQTBZWRO4Z6Z2I0mzcD7xdC3IAK4Pbtqj8f9r6lf4aUsjpl9pPAXVLKq4QQn4yej1mtK275Vu2uSd6znCE9XmPNsu2Y1DU5TJ+mto0JPQglblTPvDlr0ZQZVNJkUyZHXXFH4jv9/d+W8vlLjuPFTQMU3SCZHAbKflI6QYaS1U8u2WH3yCsWNAOwsrtIW0OKTb1lmrM2XhjSnXf5+KmzRv38F/+0hD9ffiKnmRuRniL9cOnT6POOws9N55KjshAl6oyH8PcV6+9AwNNfPTt5fMTHb+Mr5x/GsZ99M1+76XmOvfIuutasTtpt9q19YVzKleEkP/y3noyTQP+D30tybULPpdi1kVT9FEBZknGuzkS6dmDiLH0hxK+B04EpQoj1wGcAE0BKeS3wV5RcczlKsvlvOz/qQext0h+O81EnAeBnwL2Mg/TjmvRxgwirrmlIEEdoOoaTTfzzpq2TrrM5dnELM5pU0saarUV0TVDyAvJlj1zaojlj0d7oJDVR/rmii39/46EJ8a+87yZWvuEQTpjXxD+WdpKLVBVbekpsXt3LxqcfpOPWzwKjB3FHgm0oS+KV86Zw98qttNbbrOkq8qfHN/CXS44Z8/OvO2EmnSWfTdZcFq6+Q52neUexSp+GUQzJ2YJ7Vvdy3oLx1VvRrYlLbqlh4hBPAJdc8yDN0+qYPaMe/WVTyUaZ2bf+ORjS3Gd7GIvMJxPZx5h+wdeSx0LTSTdPJ6iUkgzcWJ03kdm4QCJjHg3jsfSllP8yxvsSmPCu83uT9CVwuxBCAj+IotpT4+WLlHKTEKJ1PDtqnnMQvhck/vlywR2iuTdTKUxbNZ52shaZeptj5jZx0twm+iO1yswGRe4dhQrdeZcZjQ5tWZuXTc3SnFJBm7cd3kp67WP8x5deBagb4XN3r6IxIvvuqPjZ8seWRxLHbf2r48WZc1QwOkRQ9Bp5aWueB5ZtHRfhw6DyA0CbdYj6X+xh1rRp5EOdho7neINe5DcvLuCNBzWPuq/hdXVqmHyI3YLD0XT1T8hOnYPQ9CG/42Qk8R3FBRe/jgcfVuIWw9RJZUwy9TZbNwxQiZrJh76LH/W7BuidiC8eT+2dWkbuiDhFSrkxIvY7hBAvjveDQoj3AO8BmDlzJq1NDqm0RfvUDFv7KxT6K0lzEbfiE/oh02fmSFs6s5szHN6e48QZOeZ66wnq2wDQVz7Cimkn4gaRqwjBvEYLect32fTPpwHoX7uVJ+9azTsf/iln/VUVR4uVNsvuuXHwRnrr4RN1jrDXPArWIZSDkC0bxq/DHoKt6sY4/GcVHvri4TR0L6Xy7AOEpQJHnH70mB+Pfbw1F8++h7hLmdGzliufDPjGp/ef0slHzW7gmjNVvkqX0UjJD9nQ77IpX+GBFV088NRGTjlyOks29FEuKGXd6gn4Xn08VTYncR2GvUb6UsqN0f8OIcSNwPHAljj5QAjRBozY+SBaFVwHcMzRR8ufXnYi3UUPTQhaMib3r+nhqXWqxHJ3ocK0nENbQ4pTZjdSZxssNgfgpVvZdOttDKxVX9F82BzmXbwIzR+g8shtSLfM5sdfYtXtL7EqSnpZmnc5IpfiiQ98iiPe+y26t+SrpI1n7pbz5M0+ltYBjz//bSMdq8depo+IKDv58S+pMd5SmMrxZ7yHwv+7hPvn9TC/YfQFVUz2NeLfd+E3zuIzr4RvjHP76t98sv7elx45jfIN/wNAy9GnEsw9lvZWA6YIXj1vPsbZC5DAks4SmyI58gMT0cBNH4+lXyP9IRBCZABNSjkQPX418HlUtPpdwFXR/5vGs78tZ57Joter0hUbHl5JzyObmA2ceHgLL7v0DJxTLsPoXkP54d/hF0vc9/Vb2bQlzwNd1b6+Z+ETfx7zu+7bWuQ1bziBTy77IcUv3wSnfIK/z7+QB1ervrE//8Nz9K9/adQWgjuK6VkTy9BYdt1bd+rz/kGnRTkG6vmr5+ZoO/8qNt30U941js/HN/0FH7wMgFO+/HdWPvzPWgvEfRDVpbCHl8Xuf/B7yLt+wuF/VavfU6+6j3c/dTd/nKVcihc9dD1+8xzecMNy7vzBjybFZGBd9N8A+MB3Wo5gRcHl3WfPY9pRszEyKdb9/UV6V/biANOOHJe3eExo+3jBNRGXWt2jXyrEPODG6KkB/EpK+SUhRDPwW2AWsBZ4s5Sye7R9HXP00fI/15R4qHtosKYlagZ+zrHT2bSsG8fSSU9x6NgwwD2dxQk5jtdMy7I87/KmNU9wzT/Xcv31dwPQtfyJ3XZDtJx95Q43q/7TUnUK37CoaYc+N9zK6ygF/OJJJRP+3MeumhQ3fQ07DykEIrr/D/3wX1n38C1865r/x3V/foETvvy+IduePVW5PU/84Ol8+2WX88e/LWXV/TdPqmvgptnHbPfePizqOf2f/Usfl1IeO+JG44AQ4ppfn3vKZS+bMroA4n13P8J9GzrapJQ7uTzffdgrlr6UciWwTVaQlLILeNWO7u+i9U8Q5zp/KH0wAJlo+XXzIxvo8yIN/Q7UJhkNc9Imq4seT/SWefslRyMl3HDT83QtVwXwdueN8C/vf/sOfyYm+8/dvYpvfPrrOz2+VkfnY8dNQZoOH51EN3sNOwdRZfAtufo11J98Cx983+eBbat6vTSgZJ8n6xqfn7KUa196kqZ5RwzpLLe3cf6axzmfQQ6IMT1l8Fx/ZeQP7QT2YHLWbsFkk2zuMt7/tkNJNWZINSud8pJfP86Ny5Sl22LrNJo601IG920tcmKTk1wMugBH19hc9re77xiro3ILm8s+v/zxE3zqHXfy8Xcfy7+N0GhkovHt1y0e97a3rujlh/evSnT4a7oK4/pctcpjuLUvzZp0c3/GtT/8LGtf/kqGm6fxNX/lnEv4r7nz2fIZl3Cgl1P/OZuzvvsQm1f3JomImiZYcvVr9vDIB/G/xRdYfumbWPCj32/z3g/S6RE+sWOYqIJrewv7HekP/6EX9l9CyzUP01kJKAWSMxY1UNxa4pRmRV6vaq8DoL/o8cLAjmftHVxn4W9axYtyPhwxddcPYAJx7vwG/rkqx6L3qXOy9Jo38fsd7HUymZbvNex+pLZj+Hzy46pmlLX0B2TO+yZLv/8zbr1pGV996gEeWdfLrZpIalRVFwXcWxiJ8CcKQgiENpZ6Z7d9/S5jvyP94Wh71Slc1lBHUHZZc+9Sbnh4A6+b28iMqWkaZueSZVpqZS/3PLVlzP29ZlqWh7pLzM+YTK+3mX78dMTZl/GV3X0gO4nPnzmPVx+kAlhz3/lj5r3i/G2s9+kXfI38ltVJAtbisy5k2T1/HjOT86D/vDlpzlHDvo36ky/n6SO38q2jTx3y+kXHtHHwW44ld+SRAJiLjqHyt+/RetQCsresYHFzmme3DPCql7WxvkfF1bryFc659hH+dtnxe/ow9gg000C3xqDOWmP0vQehabS8+tXge4RhyEXAlIOaKWwpEAYSgqj1XMHjjJY03W7I033l7e4v74d0uwGNps6xx09n8c/+uIeOZOexMMo6TjVOY9OzD2/zftO8w0jlWih0Kj3/9y87kYv7y6x7WFUw3J5s78P/dtxuHHUNEwHtiT8THv267b5/+yolbe799CFsvOcRdJHcEhzdkGLmKbNpPus8CouUpf9ib4XF572fnFfm0nd+mId6yizZ0E9rnZ2UGk9bOltDyUlfvHebaq37A8bTGL2m098L6PySavTc/m/vpevmX5Oe1ozdkOWYj5zP1ideQDMNQs8niErSzn/1AjRTnY6LTzuO1Tfdw7euf2ab/d63VakDFmTNfYLwAdb2qbjFP792Lq/6nL1NK7fnvnEu7/7dc5x96GsB+Mbdy2iZkWPJ1YNEb/SspftX36f+8quS117cNAC0YfSu571/z5NzTL786vm7/4BqGDdGI3yAsztUp6n8pi3Uz2njXa+ex0uPbaLJMWmc18C0Ew8jnP0yrGhF7AWST966jNZ6m9uf2ohh6hTzLpl6m3mtgzWvMhmLMJS89sePE/ghL/zjyf2mPLfQtDHdOzWf/l5Ay6e+C0Dndz5O3ezpuH0DmJkUwk5RP7cNI5Oi8cjDkIUqRY+mozdPQz/oBMxTLuXqC39O9+PP8IXP3pZsMtMxWVfyOHfVY3v6kHYax7UNBq9OPH4Ga7YWePlX/8HTf7oBzbDove9qlq3qwYrq/QwUPb7yL0cmn8nf+1VE99qE8Au+MgWvPm8RwivhN8zguhNf4rrN6qZfsrXCoVOURK7w4/9H7g0X4zfPG3FsT2xRk+jRU3c9wFbDjqPls8qwWfUmKHb24pd8mhyTGSe1Y9fbBGUXK/ATpY8XJfrd91Jn0vFt3TNPM/NlR5Ctqkbr+iG2Y+J7qkl720GLOPKTt7PhyX9QGVDCin01XiRMA80yR99mErt39opOfyJxzNFHywceeGDM7fI//DS9S9cx54ovELz4MMK06H/yUUDN3KnmHJXXfoQVPRWmpA2mZkz6KmoV8MIrVF2TzZvydLsBOVPjzet2uT/xXseeyLbU+zYiQh/viTvRTzwfseIxxKxD8RtmsP6/VGrYjK/9bLeOoYaR8c+TTwfgqn//JreeY1N6+G8YU2fhrl9JuasPM+NQ/8rX409RE/YqN8XTm/Pc9vwWXlzdQ/eWPAC6ruHUWdQ3KjdiOmXQP1Bh64YBzJTq5hYEIZWB/iROZDhZ/FKewC3xx2+8fY9N+k46vcs6/Zveed5lR04fvbH8xb+7k3tWrK/p9PcmSh29yCDEb5iBNasHd/rh5A46mTCtkiw2lwV1gKkLthZ92l+6ldYZiym2LOLU76ns2gc+cDX3PbWFaan947Rt/cAM3Bv+B3vuYg79rUN+y2Bntolaipfrp+P0rMacuZCeX3yDusWLyD/8D7L//sURyd7oXQ+bV+AfNHKHsRomDic9eC93LTqBD/zP+/APvwLzNZeBDEkdXsbxK4ieDQTZFrpFhpIXYmkQSsnqjjy9HQVaZ9TT313CsHRSaTOxbl0/RGgCp85KCiGWCy6aaVHq3cy5bz0z2e7sw6Zx7Ia7CKeO7oaaTBDa2JLNyWvnH0CkP7C2g1RzHdrjfyZYeBy9Hvh6E0Q5G+Ug5IHlvbxyXiP1d17DptMvY0ZxFfbTf2X9H28GoK4ti/NsB1Pt/eO0xSnsEnhueDYOYG14Grf9CMyOpXgvPQYvfxuPbipy2G2qpK1z8WfG/A5TgOjbQlAcoO6yL6vXNn1pu9v7DTOgQVUIXZ/3mZHdP871ZMWrlqrAfvbkyyn/fh6PGAs4qHkK6aBIb3o6Oc3D9yRLo5IlfWWfI2Y3MnNKhiAMWemFGKZqH5g0DHJDCv1ldENDhpDOWuiGwKu4GJZDzlFVaT9y2lxmdzzO6l/8kfYxYg+TCcLQk/jfdjGJ3TsHzB3Vs7KXhS+bh97YQgg0WBr9bsiTm1XC0pS0xUkzczzfWeSIV/8HMzuew2tdhGamEfqfAGh/+UFc0phi81Mj1oHb7+C2q6Rpr3URtC4CovjAOMi+Gt7coTOK/TbVru/Ft7+Bg375p+1+rkb4ew73ZJZReuR2jjvRoKgdTF+YpuAGOGkbIwxpi1qGzqhPUWfpPL95gIdWdDFjWpaNnQXcioeMWobqhoYRdZZKpU0sx8D3Ahpa6zjpnIP43xOVP9x/5Jcs/9MdzP7W/+2dg95JiFrBtX0Dx9xxO9amJXQ0H4wfSEQ5pK8SYEaqhPmNFilN0pZOo5V7CZ0c4pnbYeHxFDapwNO0U4+m6bhjmdfXtTcPZdJBf+Y2gp4OjMNO2W7Alr//Ak4b2i5yNMKvYc/imDtup/7kyylevxUntR69YSa5oB+RL2MbFo05lcQohcZCW+PXYZaBsk/a0mnK2vQVXbZEvXvdio+ua8hQYtoGuqGxeGEzXXmXb5+axXvsbwAUV6/e5wgfxunembycf+CQvvH83aybcxp6ILlvTS+nzWnEMQSHt6pCUiVf4pghWn4rUjcJXnwYjns9+rqnaD1qAQBuVxd0KcJ//rxzOOqWv434XZs/exnTPnvtnjmwSYCug85CA5rya9Ge+ivrf/sHpn/5x0O20TL1hCN/vIZJgv4HvwddKwleegSrcRXu0qcQpok5+2C0+sFifdJyuGhhO2cvaGJZVxldg3tWdBFElv5AWcd1Qiolj6MXTKEuZXDY9Hq25Cto5f5E7tj51DIye+VIdw1CV+6s0TeavKx/wJD+0vaX88x6lYjSmrUZcAPqLD1xvVm6gMBDpurQBrbAwSchn78Xv1RIZnVrahuy0E9+3WamLN42el+MpIwv/v4Zpn12jxzWXoEIVZq93rcR/7n7aTz5TQivhLQy4HvbED6AtuAYKr/7Cuabx+x+OSoePOE0Tn7477u0jxq2jx8ediGX3nYV7vJnQNPQsg2Unn0IPTOowddzzZhzS0zJ99LS1EbgNJI7ZBqHt6lub3e91EnJC5jR6HDavGZm5mwMTfD31b3Ins1suE39fnO+8+u9coy7Ck0f26dfS87ay3jhX87n8O9eTW7WHNb1V0gZGkFkdsbF8oQQqpuyDAnq29CKPWjT5qIVekg1LwMgv2w5dkMdZsZh5jd/key/4EvqSx2Ur1PFGEqFHa/hsy/A+91XsF/7HkXugAhczIVHIVc8TNh+CIQB+pxD4bGbCI89H4Bl776QhT/5A6JnA8XNXeR2cQwnP/x3fj79KN658cld3FMNI+HfNj1JAFgzDyZc9TR6+yLyLz5P/zMrkm2aD50LhoV2yClIKdEqA0yrbyJnq4mhLWvzfGeeQ1qyHOkMEKZMVhUEF80MEV0Z6ue07aWjmxjUfPr7CMI1S5iW20z9vJMBsEIXqYHSrkT/dRNR7EGLLFn32fvpfOAxMu1TALAb6hLFSzXqwiLBU3ex4haV6LIvJW7tCKwLPkRJs0mFkeSpr4NgoBdv9QuETz+IOWUqsq6RcKAH/49K4TPv9Scj7/oJlUqJzqdXsvSVZ3Lc3Xfu0jjGIvwPpQ/mf4sv7NJ3HOjwG2fxmTnn8qWHv0/9CS8neIeqeC5/8EkqvQNkZy7EfeJ2jJZ2AOrsOnKeUvi0ySJHzG2hU1hQ8dD7N7NAhrgP/hmvv5/swgV77bgmAkKMJyN38uKAIP2Df30T9xxxCoe85SgaFv0D3baRja0IKwWAMWUa0vPw1i3FmHcoQU8n/ua1+Pk8Zn2a7iWrAIZY99WQmoG7YQ1HvOdVPPzVW7Z5/56DT+SMFx7afQe4m5D/4afJLliASNeD7yLDEL1zAx7g9RcpbO5iYO0Wyj1lAjeg1FPGK3gEXoAe1WGx623MrEm62aHS7yJDyb2HnYxuabz8ift3aDz/PPl0Tnrw3jG3qxH+xOBzfUuQ+U7CZ+8l9/gf0A8+iaXPrUZPmXiFn1LpzVPq6EVPmRiOhWYqVU7D/HbCIKBpzhzCVJqV/3cj5Z4yXct6MBxFOSf+4+69eGS7BlVwbayM3Mk7KRwQpA9w2tMP0P2Nj9C/ahOGY2F1K/++0DVYugSrsQEt14y/eS0Dzz9HpSeP4VgEZRc9ZW13vxuvuIS2M05kxU0PcfCvb+LlI6wE9jXC/37rEbz9M+eSaW8hLBUItm6h3NWH0DXKXapsRWFTF6WuAsWtJSr9FQI3xCu4eFEWsxYtb8NAUoomBacxRaY1Q6oxRaWvwt+POIXTnh47mzrGeAi/honFhq9fSeuxh6Jl6vCfvJO+9f04jSkG1i+lZ1UvRsrArrfRLZ0w8pn2Lt9Cw4KpdDz+EsWtJQIvoGdlL5u6SnRWAv6z8+m9fFS7hvEVXNtDg9kJHDCkD1Du6sMrlJFhiFWXxsw4yOhCtXsHsBt6AUhPnUJ+w1Y6nlyF0AWlHiVFmz7CPqd/+cdsvOIS6mfvWCvCyYycqbHpkaVMOzZMLBZ3oIjbXyAou5R7Vb0ct+BhOAZuwUPTJZqlY+kagRugRzrtwAswqzKYGxe0EJRdil2FxOqrYXLip21H8boPvpye51fQdNShAGi6RnFria1Lu3m+v4IuBAdHq7l4dVfaWmLrS0rl9tLGPFsqPp/qeW6vHcdEY1wF1yYx6x9Qd507UMQruISBpNJfwa4vYtWpmh9hlyJ/qy6Dly9gN2QJ3ICnH9rAv24Y3Ye8/sHVHH/vXbt9/HsKj/aUmftiF2bKwMzYOK2NpJqVMqMShpgZtfIRmqBvTT+hG+CVfUI3QOjaEDLX9MHWckLXCMoudkOk+Q5Cnn696rDUs7KX0597cE8eZg1jwNJgYO0WWo5ayPM//DMAgRvwwrMdlAJJIJUQYmPJoykIcaKJXtMFgRuyruhx6ean9uIR7CaMQ6c/mUl/8jqedgPmfOfXVPoq9K3pwy/5lLYW8QplvMJg/XzNMpBBSOj6VPrdMQkf2K8IH9SNXNxaomdlL31reih39UW+W+W3NTMpzEyK+lmtNC1sxI6aTodRIXatqn9odSKLpgsCz8crlvDLLrqlY6QMjJSBbh1Ql+I+gXmzldbKzAxtkenoGroQFIKQUhDS50X/yz59ZZ9y2WdrxacvKsuwv0GLyjCM9jceyaYQ4hwhxEtCiOVCiE+O8P7pQog+IcRT0d//m4jxH1CWPsDSRzcxrS1L35o+7Ho7ISqvUMEvlBCaRmFzN7ppJD7KAw3fKLzAT9uOQl/fjzMlTeCFNC+GullTqZ/ThozK6xY2dOIOFMhOz1LqKeOXfFXUrhQm1n48AchAErgBQvMSl1oYyGS7aUdOrlaTNahmQ0bK4sUbHmBgo6qoWemv0BnFbdxQktE1HF1QCiQ5U/3WpeieKof7dgXf7UG5d/QxNhqd9IUQOvA94CxgPfCoEOJmKeXzwzb9h5TytTs/2m1xwJH+G9c+zu3zjiOQkgZdoEUWZkxETksZM5NiYO2WXZYW7svodgOs/gpT3ZDQDQjcAKsuQ6q5HrtRuWeclkZKXf1YmTJ107NousAteIk7J0boBsiUgV/y0S2dck8ZoQsq/S6pRqWgcvPeHj/GGkZH4AaUe4sUthTwo7653SUPRxdsLvu4oWSr6+NJnUBKVqtQD24oyRoan+x+di+OfjdC09XfqBjT0j8eWC6lXAkghLgBOB8YTvoTjgOO9EEVhOoreaQLXuKaMLMWmq5R6uwhv6nvgK8Ls66kSLizEtDu+rSWfTY9thorYyaro2xbHe6AUuakGlPK2t9aQugCGW1jOEb0mZBSTxnd0vHKPjKQCF2gm2rSXfyzm/fOgdawXVT6XXRTxWiCyFXTaBlYfoAb6oSo9qFuKOl2B105uoAv9i/ZS6PeA9B1hDm6ZHMcVTbbgXVVz9cDI9S65SQhxNPARuBjUspdPrEHJOn3lX10IfBLPj0rewEodBRw8x7dJW+/aJAyEdhSUdZdvx/S7YZM6y1zT2dxVB38P44+lb7uMm60tHfCkFLBxclYGJGKR7d0AjcgPSXN1pdUMbt9O11n/0TgBlT6KxS3lpKYS+O8BjIFl6ayT66zSGclSH7rvK9Wy1/N73Zjde9i/Jb+24UQF1W9eJ2U8rrqDYZhuD/sCWC2lDIvhHgN8Cdg4U6MeAgOSNJ/49rH9/YQJj22R+znj/G5lz9x/5CM2E/XH4obSpzuElMsA11Au2PQ2l7H2kc3AruewVzLwJ143DT7GHQhmH5Qc+KCk0FIsUv5cEq9Kis7Z2rjEjvsV9C0sUlf+fR/KaX8+na2WA/MrHo+A2XNJ5BS9lc9/qsQ4vtCiClSyq07M+wYNclEDROOagL+Yv8S3FDS54WsKLiJ2+jJ5T2cu+qxCSlZ8b/FF/hQ+uBd3k8NCt9ofhl9XkiuKUWlv4KZMpKYl5EyEJqGbmhYmsCaxM1CdheEbiAMc/S/sdU7jwILhRBzhRAWcBEwxMcphJgmoh0JIY5H8fUu13U/IC39GvYs9oQVnjU0fjTtSID9Uxu+B/HRrmf4fusR+CU/eS2W3foln8ALkIHk/DUH6IpZjNvS3y6klL4Q4v3AbYAO/ERKuUQIcVn0/rXAm4D3CSF8oARcJCegqXmN9GvYL5D3wxrZTxA+nj2ERVmLMJBJnoWski/LQHLmikf24gj3LsaVkTuOLrlSyr8Cfx322rVVj78LfHdnxjgaaqRfwx5F7IaZKOv/49lDmJ4yaj79CYYbyiTHQtMFYeQJPvWxe/fiqCYJxu/Tn5So+fRr2KOYSHL+atPhTE8ZfGjrMxO2zxoU4UeKW0X4gZoANH3yEtmehNBNhGmN+jeZSb9m6dewx7GrxP+h9MFko0Dix/fXBKC9iBZbx9IEhmOoInoZky0v1fpCJ9A09TcaaqRfQw0Ti7wf1lw6uwHfaTmCmY6prH03RGj+Pl37frdgPO6dcfj09xZqpF/DPoca2e9edLsBs9Im6SkOuqnzxDlnU47Ki9f6Ew9KNkfdZhJLWWukX0MNNSRwdMHc+pTqeJYxOeqWv+3tIU0+iHEUXJvEln4tkFtDDTUkuHTzU6oiqhdQ2lrirkUjlYM5wBGXYRjtbxL79GukX0MNNQzBmSseQQaSfGeR9BRn7A8caBBiMJi7vb8D1dIXQvxECNEhhHiu6rUmIcQdQohl0f/Gqvf+O2oo8JIQ4uzdObYaaqhh+zj9uQd59cpH6d2U5875x+/t4UwqCGMcks1J7NPf3Zb+9cA5w177JHCXlHIhcFf0HCHEIaj6E4dGn/l+1Gighhpq2Es4d9VjB3T27YiI1Tuj/R2olr6U8j6ge9jL5wM/ix7/DHhD1es3SCkrUspVwHJUo4EaaqihhkmDuAzDaH+T2ae/N9Q7U6WUmwCklJuEEK3R6+3AQ1XbrY9eq6GGGmqYPNAM0MdoolIj/XFhPE0F1IZCvAd4D8DMmTNH2qSGGmqoYfdAaOpv9I32yFB2BntDvbNFCNEGEP3viF4fs6lADCnldVLKY6WUx7ZMmbJbB1tDDTXUUA0pBFJoo/5NZuyN0d0MvCt6/C7gpqrXLxJC2EKIuai2YLUIUg011DC5IMYRyD1Q3TtCiF8DpwNThBDrgc8AVwG/FUJcAqwF3gwQNRD4LaobvA9cLqUMRtxxDTXUUMPeQq3g2vYhpfyX7bz1qu1s/yXgS7tvRDXUUEMNuwbJ5HfhjIbJFMitoYYaapj8EGKfDuSOi/SFEDlUwlQ7SlGzEbhNStm7+4ZWQw011DAJoelKtjkaJrF7Z8w1ihDincATKN98GsgAZwCPR+/VUEMNNRw4GEO5o1w/k5f0x2Ppfwo4ZrhVH9XMeRj4+W4YVw011FDDJMV4dPqTF+MhfcHISVIhk3k6q6GGGmrYHdjHG6OPh/S/BDwhhLgdWBe9Ngs4C/jC7hpYDTXUUMNkRJycNeo2k9geHnONIqX8GXAs8HegArjAvcCxUsrrd+fgaqihhhomHcQYtfQ1bVL7QMal3pFS9gA3jLaNEOKfUsqTJmRUNdRQQw2TFeNR70zi/lQTqdNPTeC+aqihhhomJ8ZTcG0SW/oTOR2NWBGzhhpqqGF/wvgKro3N+kKIc6IugcuFEJ8c4X0hhPh29P4zQoijJ2L8k3cNUkMNNdQwKaENWvvb+xsDUVfA7wHnAocA/xJ1D6zGuajCkwtRpeSvmaDRTxgm8YKmhhpqqGGCMJ52iWNLNo8HlkspV0opXVTM9Pxh25wP/FwqPAQ0xGXpdwU77NMXQtRXf05KGbdDfMeuDqaGGmqoYbJjfDXzxyT9dgYl8KD6iZwwjm3agU3jGOZ2MW7SF0K8F/g8UGLQfy+BeQBSyud2ZSA11FBDDfsGxlNwDYC3CyEuqnp+nZTyusGdbIPhcdFxdxPcEeyIpf8x4FAp5dZd/dIaaqihhn0VEkE4hiUfJWf9Ukr59e1sMp5OgePuJrgj2BGf/gqguKtfWEMNNdSwLyNEEsrR/+TYBvmjwEIhxFwhhAVchOoeWI2bgXdGKp4TgT4p5S65dmDHLP3/Bh4UQjyMyswFQEr5gV0dRA011FDDvgLJrvtYpJS+EOL9wG2ADvwk6h54WfT+tcBfgdcAy1EG97/t4tcCO0b6PwDuBp5FFVuroYYa9mFMv+BrGE4WAMNy0C0HK1OH5Zik0iaZ+hT1OZu2BofWOhuA9iaHS4+ctkPfoz93B9L3wPeQngvhsC6okeJFmCbCSiESFcygI0LoOt78k3ftgCcIUkI4AVlJUsq/ooi9+rVrqx5L4PJd/6ah2BHS96WUH5noARwo+H93rgTg82fO28sjqeFAQP3Jiiv6H/xe8lyReg7NMBGajgwD/FIeAK/QB4DQdHTbwUxlMTM5nLoMf7zq1SN+x6+WdNKWtWmrs2lyDJpSqvKkVelDK/aglQeQpQFlIYYBMgiGVqcMg0HC1/VtK1eGkW0Zkb+54sEhb0vPI8z3qk2PHa523H0IQgjGYH05iVNVd4T07xFCvAf4M0PdO93b/8jkwM3Lenj9wkYAvvHgOnRN8KETZwzZ5pfPdvD2w1t32xj2BtlfduMLQ55fe8HBe3wMNewd5Gap33ruO3+Mmckx++TX8ezXz+HYK+8iCEK8io9XLhNUSgCEvkvouwCs/uV7trvfIz95O9mGFNmGFM25FPNbshzUVsfhrXVkTEXOtltC9G3B7+lAVspDLHsZhsMsfQ9hmEiswRruWgja0NVAwqFhiAwDwiNfswtnZ9cxFqdPYs7fIdJ/W/T/v6teSySbkxkx4QM4lk69PXjYf3yxC2C3Ev7eQo3k9w90lAJanZHrt5/7g0e59b3HJc+P+dSdPP6lM1l3w/u32TZ+76gr7sAwdQwzg8ykAQhDybNfP2fI9kd8/DYMU1eByciylaEk31um0F+hxy6wfn0/j66wOWhGjuPmqPvsuPYmGqe1EExV+9EF2IZGnSHR+zejdSvpedC1GRlNOvguMtQhCECPjlUbqjORnoc84Y07cOZ2D0Imxr2ztyDkZF6HjAPHHH20fOCBB/b2MHYZZ3//YW77j+G5GROPi/7vaXRNyc2CUHLDvx7B+T99giCU/OWSY3b799cwMszNL+BNGzpJ6wNbCOqm7vC+jrriDp788lnj2var968F4OOnztruNpff9CIrNvUD0LOlwMOfOYO286+i0LluyHZ2XRMATvN0nMZpOHUZnKyF7ZiYto7tmDRnLWY3ZwA4qK2OhU0ZZuVspqZ19P7N6PlOwnwv0i2rOAAo15DnIcsFOG3XckCddPpxKeWxO/t5IcQ1f73r3suOOmb0Xbz9LRdy1+23tUkpN+/sd+0ujGnpCyFeKaW8Wwgx4hQrpfzjxA9r/8SpV92Hrmv8/b9O3ea9iSD8Yz51J5qhoWkja4iFBrquoRkauq5hOwbn/uBRZJUlt6ewpya5fQXDCR8gqJuKVujCz0xB2wGHQe/6tcnjE79wD3PmNXHGwa0c155jdi7FFCtEG+hAq+T5+KmLxtzf984/aJvXNt00WB8sjh9UBrqT/72rnyOVayE7dQ6ZKdNI19lYjke+t8yGLQUAHl3RRS5rMaMxzcJpWRY1Z1g85WDaWwOMrtVQ7AVQLiKAugbkc3cQHDa+CW13Qcpx+OwnsS09HvfOK1CqndehDkUM+18j/VFw3Gfu5tHPvRKA+z/5ilG3PemL9wIQBuqKefgzZyTvHfSfNyPDgJe+d0HyvBqh7xJUSnjlPF6hD82wVHAMFZyL/2tVwTLNsHAaW8k2OGi64KQv3ouUksCX6IbATplYjoFlG4hoIvEqPgD53jKVkk8QhGNala/61j8RQkS9J2olmqrxm+e3ognBmw9uTl772TMdpAyNAdfnhze9wEt33Uzou/Q/+D0aT/8YAD33Ds35Oe1r92M5Bne8/0Q++aFBf/dDV6prSCv1oPetwUsdDA//CWmayHQ9RvcmRFMbXusiHt1U5EM/fZQHrjhth46hOlhcjXJfJ+W+TrYuVc/TzdNJN7cDoFuO+kxLM+vrLJ5cZuJkLea1ZjlxXhNnz18IOWj1OtG61yErJfyDdmxcuwvhONQ7k5jzx3bvCCE+yrZkT/QYKeU3d+cAx8JEuXeq1Q7zLr6elddfvMv73BEc+uG/Evg+7kA3lYFuNvxBCaVazr4SACOVwXCyGJaDZlgABG4psa7cYj9+KY+MgmS65WDXNaJbTkL+oIg/ngTMVJb6adNIpU00TSTEDsp3K6LXDFMn9EN8L8D3lKKiUvII/BApJUIILNvAcgx0Q0NoAl1X/ljd0ND0kYm++jsNU8cyNCxDS5QRv3vHkbt+YncCtyzvoRgdZ9ELKLg+rh9uE/zfnVh8+Y2UejbTu1pVN9Eth1RuCnZdE1bkRll6zZu2+Vz3Nz6CXyhjZFR7i6aPDt6et6/qY3lXgSUb+lndkaevq0hvp7K61z96exLIrUZM6ONB9T00fAKIoVsOjXMOA8BpbMVyTCzbwLR1MvU2Jyxq4d3Hz2SBp1xHcsMygpedPe4xjIWJcO/cdPs9lx05hnvn4osu5J47bt833TtANvq/GDgOuAlF/K8D7ttN49qt2N4F6TROY97F1wNK9WCksui2k5DkC986b8LGsPA9vyFwS/huibW/eh9LrlbW2dTzvpAQNyhrXEQBraBSIvRcZBgQuGV8tzRUaldF8JqmE/oegVtOJgnNMNEMC9120KP/bslDhhLD1BMC1jSBpgt0TSBDSbng4lUCvIpPGEpiQ0EINSHELiO1QgkxrKp96SIhdxlKwojQ9cjFFG9nGRq6JghCia6p737tjx8f+vtYOq31KXRN8O3XLQbgw7csxdLVZ42qfX369Dk7/Js8uCFPxQ/RhMDU45VNyObeMr0lj8tufAHL0HBMnbSlE4QSPzoeIxpzEMpEqfVsZ5nfPb2RjoFKEkc5fm4TZT9gyYZ+VmwewK34GKaOVwm48wMnJmMp9Wxm7a/ex/QLvoZfLqBbqSEulZHQ/71PDpL833+hxtW1Er9ZjefVc3O8em6O+g98SW0/hNDPBba9N+pPvnzcxF+93faIP3BLbF36KKBiAJmWmZjpXHLNf/MtR7KosBS/Ybra/ohzEJMs7ijZz907UsrPAUSN0Y+WUg5Ezz8L/G63jm43YPiFGAefrKxSHrjFPoSmk2mZhWlbhEFI4CuXxqL3/Z7QcwlcpThY9fNLAJj2ui+x+c+fSvbZ8IoPI8MAGQYITcdwspipDDIM8aPPAuiGhdB0pr3uS8gwVC4azx1imftuCcManHgGLfkURiqD06gCfaHnRrI7L9kmiCy3INJi66aFmdYJfTcZnzvQrfTbdU3Yjqm2G0LiStkRBCFCE2iQTELxdpom1DZ+SBgKykUviRGYtkEqo6y52MWjGxpCqM/E6+Sip8YshEA3NGT0nYEfIiO5duCH1KUMcmmLj/11GZahUWcb6JrAsfSEWAG+88gGHEvH1DTe9bKxlVm3LO/BCyVeIPHCkI4BpUre1Fum5AWUXB/HMhLCj1ckxaIKNnaWXIpugOuHnPuDRyn2V1jx6ONU+rYy5aATaZ8/6L45qK2OeS0Z5rVkkkkE4IrbVwDwxMou2g9RPv6NN/4XZ1z9IJl6lRz1vw+t56jpOQDqbJ3WjMVs2YVYt4SmV52Lj1L7/G9FJTLd/f01rH/2z5R6ttB1x5e3a/DAUIKv3q768XgmgNG+A6Bp3hEAtB20iIaWDE31Ng1pK4kduByabDvZCB9ARqUWRt1mD41lZ7Ajks1ZqKboMVxgzoSOZjei8fSPIcMg8SXKMEAzLDpv+0KyTe7UD2A4WVL1UxjYuBxQmYDrf/vBEffZcvaVdN72BYpdG7d7ocswwCv0JRb5EGRyGNF4qq32eCygLH0vCBC6nkwSoEhfM62h/nrTIqiUBieGEcYTE35s/VfvM4bQBKEf4gVDE69jIg/8MHluagKq/P1uxSf03ap9prBsHV8EibsIwLINjCpXTlDlPopdS5om1GrCUpNMU71NS12KupRBNmXgWDoZU8fUNbwgxIv3JSW6EPzbEVP544td3Lysh5ShkbV0Kn7IGbPrhxzXLct7OG+BmvRvXtaDjYZer8apa4LuvEvR1GlIm8nkomsC1w+HTDS9xQLd3SV6Owp0rHiR/vVLSeVamH1wC68/VrmGuvMuKzsL5Mtqssg5Fs1Zi2zK4PhI8jivValbfvlsBwAXv2oBq7cW+MCfX0q+F2Buc4YghJLTSLnlZDQhKG4q8o27l/HQrQ8D0Ld+qdK1++6YVvv2fPMjvbcjbp/hcKPrPFOf4rj5zePOYfnE35bzlXMW7PT3ThQmogzD3sSOkP4vgEeEEDeijvkC4Ge7ZVQTjJazryT0XXTLwXSyyDCg3Ne5jQ9zewQ90oUevzaWVRPDSGUxbAc7NyUhRL+Uj0g4TMg/qJSGjCt+LDSd0FNkWh2g1W31uZi87WxTkm2ZWPPR8XjF/ui9MCECva4JzbQIfZfSgPouw7JJpZXUzq34eBUfUdUUIowmA7dUodjrJudOZXKmMEwn2VZogkrJp1zw0AwNw9SwIz+u+pxM9qkbGrqhCF4IMRgvsHTmtSov44wmB8fSSekapq5hagJD15RLJrK+3nrIlOT733jQoIW9PcSED2BGhJ6OrO962yBfl6LiB2iaQBcimWC68y4DZbUK7C5UKJbVeUrnbFrnH0TTrIU0tGQ4Yk4TDfEqKiLtfNljoOzTlXdZ2cmQyUPXBNe++XCypsamvM/dq7p5bkMfy5Z0sOWFJ3CLSj6ZaZlJy9yFNLRmKA5U2PD8C0kMYCQMd7nsrL9+LIzm0we1etkZfOWcBbz7d8/xhiPbydkGy7oLO1wSYiKgMnJH32YyTwo7pNOPejS+PHp6n5Tyyd0yqh3A9gK5c9/5Y2DQgi73bSX0XTTDove+q5PtxkvauwrNsDCdLGZ60MrUTEv52A0L3y1t10UTB2jjIGxs1ZtRcBdUUFaGAb5bQmg6ZiqL4WSRQYBfVu4dt9BH4JYS0ldxgFTi44/dQWYqi2ZYyfji86Ybiqhjv3zglpI0/ngVFcdAbEe5I7INKYJg0EWTrreZ2pqhtT6VEB9AyQ3QNYFlaKQtHccykuezm9NJ7ZcYuhCYmsA2FPkPOddCcOacodb8eHDrit7kcbxPXUAgYUN/mfuXdfLUM1vQDcGU6fXIUDInmowWTs1ScgPWdBVZtamffG+ZIAixUyZCE5i2mkQacmqlAtBb9CiWPHw3pFL2cEt+9JuHyWpLMzQy9XZ0DiU9W/L0bVgGQK59IdPnNzG9JcPZh01joOLzx4fWsn5ZF1uXqnjIlluu3OHzEGPa675EuW/rkPtlRzDaimFvYSICub+99e7LXnb06Dkt7/3XN/P3O/fdQG4CKeUTwBO7aSwTgoP+82byW1YlZKRHhaR02xmSar6nEfoulYFu5ZaJXCuxhDImaWkp0pbhsIJUjVORQZBMCF4pj5nKIDQ9SaOP95fs17AwU2kC309WBmYmR7FrQxIMjv/ifejRuCAi9HJekX+0igh8hrhu9KpYQ+xWciM1Uegrv7PlmKSzVpIb0NioSK/k+gyUffKRT1xKSSBUIDRW7Xzsr8v4+msWAsrtAoqEK35I0Qso+9BX8fm3I3Y8gakav3y2g6xtYOuDweSOgvLpD0QS1WfW9fHUM1so5SvIUFIacDFMne4t6jrL2jN4xYJmzlzYwqMbevnz4xvo6YgnRJmonorRJNfW4NAbxT5sx0A3xBD3WV9HN8WuDZipLK0vP4LXH9NOS9piQ3+Zvz7ZAsBz9zxM1/Ined5yWLHqWE45cjrvfeUC+k6Zw78dseNlClrP/SxzTzyNY49s43vnHzQkTrUzGMvi31cxHp3+JAxFJNjvMnKnnvcFQt/Fj4gsdnUEVQFUGNvi2B0Xa2ztD4eRygzR1WuGlejpYwKPH+u2kmwGbgkZBJT7OwGS1YCSajYllroMgmQSif35MYavDIb79uOxxC6kOBcg3k4zLXTLQdMEvqtIMvTcaOWg7AnTNsjU21iOmfjqFQkGmLaRSDullNz1wZN28QzvOm5e1kNX0WVDj7peugouJTdISlq86lv/HHWcJ37hHkzboH9rkb6Na7BzLaRzWbTIPRYEIZ3LnmNg0wrmn/4GZi1sJggllZJHOZoAfTeg2F8hv7UjcUMKTcdpnIYMA0o9ynisDPQMigVsJymkNpbKZ3uoP/lypiw6jnlHH0SmzsL3Q2759502iiclJsLS//Utd192+FGjW/rve8eb+cddk9PS3+9If3sYzScZB7j2hFWSaVGNcGLyjn35oe8ShkFS4jYm5zjoWm1hV5NtDK9cTCx4oQ0qdHTLQTPVpCEDNfnFRB5b+fG+Y1dRHOyOYTsmtmNSLrhUSorcq1cY8ViVy0mtLuK4yIqf7Hja/HGfuZv6JodUWvnB91XiGX5dmZkcVroezbCom76A3JQ6QMlabcdMEt98L8St+NQ3Opx21HQOb89R8QNePqeRnsgF9Ok/P8/Sx1Ylk308gVe7CAE6bv3sDo15+gVfo9i1ccwV8d5204wXC9/zG/JbVgEqi3giSP//brlrTNK//B1vmbSkv8ON0fcVDE8U2d5FuqMB2V1BvOqIH8eJVsNvWFEeTLIKPTex2CEifMNKcgji7fxSnlLP5kgeOmi1q4QuRZ6+6yX7MCybwPfxS3nKfZ0YTnZQ0TOM9AM/pFLycCNSAhJ3Wby/oFKi1LOZ/kLfLrkF2i/8Jk7jNAa6s+Sm1A3JSt7XEF9z8XmtFglopoXtLGTTs/+kacHRtMxsxIyC24YpMUwNt+KTTRm8tGWAe57exNUdBYQmKBfVeY+FAGEUlI/3G+dsABz+sb8xdVYDnev76Nu4ZtQKmsd95m7qZyxGtx28Qh+BW07Gb6SUokiGAfktq3dIv7834UUS7J0NHo+EA6EMwz6N4Rfn3vIxaoZFKjdl0JqOrPfQd3ELfXjF/kSjH3pD1TsxaZhOFs0wCaKgq+mrmuegshvTzVOpDPSTqs9hWDrlQrT/YhygTSNDS93QFeUicgt9g5ZiFPANfRcvelydIxD6LmY6N0Q9tE38ASWPHV4mYLzY8IeP0H7hN9GLDi9+5307tY/JhpGEA31rX6BvrSp9vfHx2yh2HZbEoey6Jppnz8Gr+Pzm1qVJjkToh5QGBpLtfLeEmcqiA5V8NzIIklVdql75/Vtn5pSrZmqW0kAT7Rd+E79c2CbA23L2lehWKnke526AulYzLSpZavZBLRTzLs/97dbkeCYb+bee+1k0w2Tznz+VTHJxZvtEYHw6/cnL+vsd6Y9E6nszmCQ0HbuuKbmhYp+4G/m+AZVdWykNcb3ExB9LK61MbojP3bAc5S/OKlWLZqgCasb0+kT2GIaS0FfKEIDAl3gVg1QmTaG7i8BXq4hM66ykbg+AX+xLCCR2EVW7mYgM/nhsYRgkAWmh6ZhOluazrkjGGU9sffd/e1znLHBLSRmK/Q3bcyNWSy3zW1aT37Iap3k6VjqHZlq4A93q2qhy/WmajlvsS9RW6ant2I5Nsbc3mZj7u0vke8sUo2QzK5PDyuRY9L7fDynj4JXyiXvRqmvEaZyGnVEGipQqYxugt7NAob+CbqWoP/ly5r3ifL7/2Eb+49jpu+eE7QSMSOQQo/msK4a4I3cVfpTENxoms9d8t5K+EOInwGuBDinlYdFrnwX+HeiMNrsiahuGEOK/gUuAAPiAlPK2Hfm+yagUMCLrfNBHH1Lu2xg9Dsi0zFTZtf6gPx4YctHKMFDEirLAMq2zcOrqMEydIBIMp3M26axFX1cx0dsLTSRlbQFMW+CWPPo3r0eGAalcy5BAsVtUrgflJw4RbonAsPBKefxyIZoAUslkFVv6cfA2KQ9RlSDmei499359h36bXVWN7A+oDHQThgG9hcHJwMzksLODOQWGk0WP1FWmpifBc8PJJvEer+KrWkiGRjpXj9naRCpjousar/rWP3Gy6rd/5SXvxPVDejsLlAZcgiBE1zVCKRno6Ex+z/xWncAtYWZyZDSdp7bTVWtvYrgrZ8TEyF1AICXBAZKRuzO4Hvgu8PNhr18tpRyy/hdCHILqCH8oMB24UwixSEq5rf9gGCYj2cdQwdPykPIK1S6RwC2TapxK4JYTV46ZyiT+VD+qKusV+kg3T6e+bR6arpHKWHRv6CDTpEggTnYKfEnfxjWAmiD8xinUNSqSTtfZZBtSGOZs+ju78Ap9VAa6lR/XcrDSylUUxxZkGFDJ9+CX8snEUO7rTMZuZnKJvt+Mksyq3UQ1jIztZb7G5zj2zw8nqzgmEPvrHQZzPeyccucUurtU3CedST5nOwZeOcByDAxLTwqc+a6K1QDMaM8xUPYoOWZSMC9eMY5VRbX5rCvouuPLO3Mq9giq3U9O+qe7vD/VI7fm3hkRUsr7hBBzxrn5+cANUsoKsEoIsRw4Hvjn7hofkPjEt2cNJC6NMRBr2WEwO1WLrOjq5bg/bJlZ7utMXDsweMPXt6s653WtU5MbMy5xnMqYFPoqWJm6ZD/dW/Jk6m36OroHdfqZHLquJRrx/u4SoR+S7x1ISN2okpDGPlzNsPDLeQqd64acF9PJoplW8ppX6MMDUrkWArecVPysRnzDTTa/72RANfmnci1DJtTREF8r+S2rybTMRGSbIiI3MMyhGchx2YtK2cOps9CEwDA10o6JnhEMFNT1tm6rChKn0iaaLhJ3zq3vPY5zrn2EDcvVbzu8u1bruZ/FK/RRf/LltBykCsbtjGJrd2KijcKae2fn8H4hxDuBx4CPSil7gHbgoapt1kev7TZUq1ScxmmU+zqHWOHbS1uP6+JApGywHYxUJvlssWtjFHS11DI9IvLYah+O0HdJRZZaua+TKYuOxk4pxU2+dwBQpRHCUBJG+vZUxsIwtUTbXRoo4bsBU+dOxXen4FUCDEsbUj7BsHSKfRXMVAozlcJ3vSSIG/puEkeIJZy6lUpcOKaTJdU4Va1GosktJnk/kn1WT5A1kh8fms+6glSuhbq2+clr4yV/gHRzO6Hvkt/aQSrXjNAETtbGilx6biVACEHj1CyptEldnU2h4NITGQCxxdoyNcvsKRmCULKpt8RAwSXtmLzxZ09imDr5jg0jfv+OSkL3Biba0g/H496pkf4QXAN8AeX2+gLwDeDdDNbpr8aIpy5q0P4egJkzZ0KVcbM9RclIAbTAHbSwq0sUA0nQcSTyqpuublCnYVpSHrm6OcmURUcDUBnoZ8tzY1eftjI5nMap1E2bRbYhxda1W/DKgwlZMdxSKfqvY5g60hisW29YqiZ+6KvSxnGlSi36D6o8gVNn4dRZUYEzG98L0ISIip1FKxJdw6soF08cVI6DujDohgA1cRq2gw+JW6dG+OOHV+jDbJ5O54sPjb3xCOh88SHq2uYnMRU714Jb8pLfMpU21fUQXSdBKDlkdiPPrepm/eoeBjapwoKrgGcbpnHyqw6hIW3y5N+f59gzDiVt6RTdIKkou/jyG+lf/9JOJ4DtD9jXm6jscdKXUm6JHwshfgj8JXq6HphZtekMYON29nEdcB2o5Kxi+6JtKmEOT8aa9bZrtiGjWGGSbp6ODAOKWzdi1zWOeUHH33XMp+4klTHp7y5R7u9LJpByX09SxKz1kFOGZAXH9W3sbBNWpg6vXCSI3Du2Y1Iueth19Ukly9B3E029DAOsdIZ0vU0YyKi5tY5pGzS0ZLAd5f4J/BCvEiC0bceuaYLAV2WLw6jGC6hSAcV8bOmHUQ38dso9m3GDvkHZpueiWykVi8i1UBno3iesvcmMYte2l3m1xNcd6EFoepKJG8NpnJaUuW6aNZ9KyUM3tCQLGlStI1AVStOWTjZlYukaKx57kXdffDrZlCpzXHIDSm5ALm3S0V/hfe9+OR89eSbDEXdum+w4/RsPUC665HvLlPtUCY9Yzrqr8MNtK9AOx2ROet3jpC+EaJNSboqeXgDE8oSbgV8JIb6JCuQuBB4Zzz63V/q4Gmt/pTTfcTlkgMbZh7H8R28bst2ct18HKIsmnctS3+hgWBHhhZLuTQOkMlaSMenUWVRKPqXeIAmmyTBANwycurqk6UjsZomtdVWNUsd3B2V3pXyFxqlZMvX1SVGzrRsHCP2QdH2KVMZU1STtwZ8tvrikVEROQNLRilBZ8LF/Np4EhAaaFBAF6mSodMexSyAwRHK8qcZpCE1n5fUXM/Oi7ybfq1spKvmecckwJ6OWe7JheOwolmfGmnwjlaF+xiL61y9N4lAqb6I+EQiks3ZSydSpUwTX2pDCsQzacinmTMmgCzXxn/O5N2Lqqtw0wL0ru3AsHdcPJ0X54l2F7wVsXaNoxq5rSHpiTARCqYrwjYbJS/m7X7L5a+B0YIoQYj3wGeB0IcSRqPOyGngvgJRyiRDit8DzKCX45eNR7mwP4yEZw8my+PIbgUELJjt1Jsd86k5aZjYntWHieu5eJaBlRg6hqWqJrfUp1mwt0NCSoVxsxC2q1nO65eC7FTTdiRqMkDQaMVNKr+9kbTRN4GQbaZnZjJO1aMil6B+o4FYCXn64Khn7XJ3NQKSxBpDhINH7bjBERWCYGpquIUNF2DIMlXY/In3fC6LXo4khalSiR6WJZWz1e4MBwOpuYetueD9A0l1sIrMcD2ToUS5DtRBA6DoDm1agGRZTFh1HXWsLhZ4+ils3JqU8Mi2zcAe6seqayDakkk5kpq0zvTkNQDZlMqPRYXZTmkbHJBs1lkmbOqGUtGVV7OiQqXW81Jnn46fO2jsnYQJx+Mf+Rn2TQ11rC2G0Yq5uBbqrCMel3pm8OGBq74BSGgx3RRz2UZVZ+Nw3zk1eO+qKOxLJWn3zYI2bYt7FLflouuoJG3eNitU1sR91zWMPJ8XV6qbNIjclTabepntLARlKUhl1oxmmju0o1YUMJbal099XIQhCjlmsVg26Jnjspc4hbQdBTUCVklfl09dxslZSzMz3wsjN4yfqHbfiRwogkpK/cbOSOJErPg7fCwlDiVfxt6vFnv/uX4yp1KhZ+WOj/uTLMTO5IUopI5XFyqg6PfUzFmPaFgOb1zL3mMMo5aMKoN0lNF2jaWoWt+KjaWoVOHtGPe2NivQb0or0FzZlyEUlnZsdAzcICSSs7VPqgpn1KWbW7T+5mnHPacMy8V2PF7/zemBiCq595Ve3XLbw8KNG3e6L//EOnvjHXTtde0cI0QT8BtWoajXwlkjwMny71cAAKrfJH8+x7Re/8qL3/T6p9+4WBuWIwJAG5yP5nqvJPsZouuRTr1KB2VTGpNBfoVLy0Q1BrjmduGQWnHIK5YKHW/JonJrFdgwWtOegPUdv0aW7P8qOjLowxURfJLLWpeCxF1TXJNsxKA5U0A0t8t/r0c2t47vBkF60XkU1K4ldOkITamKKxmXZBuWCi2kb2KnBpuWGqSOlxKuoc2baBr6neuJ6le0vi4Wmb5PZCYrEqssv1zA6+h/8Ho2nfywJkMctMP1yIVHnYFs0z56HF/XUBbVa1A2l4DJtnVTaIlNnkUtbSb/gpqyFqWk0p02ylkbG1HADiWNqDFRCpqTVdwZS8uimIse1pffCGZh4xD2ndwf8UCZd2raHCdDpfxK4S0p5lRDik9HzT2xn2zOklFvHu+P9gvTrmhpUT1c/xDAbkuYVuwOBH1IuevR2Fgj8EMs2cOoswlDSHCVB0ejQ01dGaELpoaNuTA2Osrqs2co639Rbxg1ClpY8hK9IuNBfIZtLkY8ssHLBSwg/7iyl6xqGodG1KZ9Y506dlRC90ETiysk2pPAjkujZkqdSqgxa+LrK1HRSBm4U/AUo5StJLELXR4gGRxgeDwFF+MMb1dQwNqqzlgud6zBSWfyyKo2gGxbNi46jaWo28dkDZOptMhmL/r6KitNEvQAsXSMXVShtr0sxJW2hCciYGoYmcAyNvBvQktbJu2pyEEIwp3733DN7Cwsu/RWN7dOQoaTQvx299E4gHE/tnV13oJyPco2D6lB4L9sn/R3CfkH6lmMoF4WpYTkGvhdSjpJOFl9+44QpDg776K3ouka5UIxq3g8mscQ3HChZ3OKZDclzAMfUyaYMLENLWuO1NaTY1FvmsFmNPLW8i851/TS1RZ2Y5qpM2y29ZSWr1JUryTC0pLvUen9QQVCJSu7Gah7D1JLJIgzUuWhoySCb04m2X4bKTVRk0F0EJA3Q3Yoq8bvg0l+NSPAxml75Cbrv/gpQk2vuCuJzN/2CryW1luLkvkLnWtyWg8nYqSRTtiHqJmZYGr4bYjqqb/DRsxtojIyeY6bX0ZTSkYApwI/80elo4pjijNRJed/H9Au+RuOcQ7FsHc3QEpfq6gnYt5QyWdXvRkyNBS9Syk1CiNbtDQe4XQghgR9EysZRsV+QfiqjlqiJzlwTyFD9yG6pxLyLr1eli313iNLn1KvuI9+rLIBKyRsStBwJajmtiN5MqTR2GQ5aw3HDasfSydrKeq6rInrL0JJtAHJpk0Pb61myoZ9MnYVXSSWrhUq0nW3ppCNfbPzZ1voUx85p5PnnOxP3SxCooG01hBA0Zy2I6qsMlH10TWVhlgsepXyFcsEjlJJMvZ1Y+oapEYZSNUfxVPJW2/lXbVfKGhN+DRODjTf+F+0XfpMNf/gIB3/wFip9ncw87BDqo+BsW1NUqC8invlt9fQWXXrzLjOa0hzSkmVqRHJNKR3LHQDfpZJuRqDyNbRJHWrcNcx9549pWfiyyPhRpD+R8EPGdO9Eb79dCHFR1cvXVZOyEOJOYKQmvztSfOoUKeXGaFK4QwjxopRy1OSg/YL0DVNLerDGbo3YGkrnsvREmbZOwzQWvuc3yDCgsX0atmMkBacMU+fYK+/CdgwV1Iw0WXHwE5QfX4YyCYKWQpdCf5lUxsQwdSoRAd/0b0fzib8tJ5CSohtgGWq5bekarRkrIXQvlJjRZHBQWz2dWUt1UfLD5IZuzlps6i0T+CEzpqTJptREsbmvzMx5jXRvLQ45F0kJ3oqP7Zj0Fr1kX6VYhx89970Q3wsoF4pJPZbknBoa+d4y9U1pOvs6D+hknL2BuMroC986j3kXX4/lmJy4cApdBTeZ/LvyFWY0pZnR6JBN5ViztchZi1uYWW9Tb6vr3wxdCH3CTDODzhtJJQR7Yrlw0sBIxeXLtaT6rJxAyzxgPAXXJMAvh9cYG7KNlGdu7z0hxJZY3i6EaAM6trOPjdH/DiHEjajSNfs/6ed7y0pt4it3hJRyyI/cOGsulm1UqVmCIVmLmhDJxRG7NizbQAhB1+aBIY08Tr3qPizHYHqLKmjlh5IglOTLPkFl0IqPtc7nXPsIMpTc9h8n8Im/LaevYbBmOZC4anqLHi11KQ5qqyNnG/RFE0i+4gO95ByLfMWnr+hy55ItpKMgcEtrhq1bi2TqbRxL57i5Tfz14XXoaLzhpFm859h2/uuWlwDoSldorU/R4Jis6SoQhGpSKrnBED2/V/EJfIk11aDQXx4SDK9hzyJWP827+Hr+0lvm+KOnJ03izzl0KrausTlfwQ1Cvvzq+dt8XuoW0tk2KcnWUP69kTL49lFMv+Br1M9YTPviGYnRB+C76r6cKKWiDMfh3tn1r7oZeBdwVfT/puEbCCEygCalHIgevxr4/Fg73i9IPwwG68bHmaxxdqyqo92MrmuYdlSK2It7xoZJopMWVRZMZUxSjkm55OFWArLDSFrXNcJAsrmryLTmNGlLZyBqdh2vCM79waPc+t7jAPjbZccnn61Oenn3755j2aoe/vHxlwPwnUc20OiYNKctlncVWBNZ8JahMa8ly6a+MjManWQl0Ffy6OhXll5z1qbk+uTSFnc8tZE1T7/ImeefyCmzG7nusQ28tK4XgPZoosqmDBzLwDF1OgfKlNxg0BUEVFIGA0WPYl65gWrYe4j9/Cuvv5g5b7+O5VMzXPAa1a/30JY0lq7hBiGPbhzY8Z3vR4R/xMdvI5VrIdugrvGhpCwn1tLfM8lZVwG/FUJcAqwF3gwghJgO/EhK+RpgKnBjlPhpAL+SUv5trB3vF6Tf311MmnbHzT7ilGvNsPDdCpWSNsR9IaIgbJx4Fb8G4EflCQxTBU+P+8zdADz6uVfy9/86lbO//zAARTcgmzJJW4MBXFC+99f++HH+cskx/Ouvn2F+S5bPnzlvyJh/8ubDhjz/z+PbuWW5asrtWDpvPlI1pegueTy/ZYBjZqvA7kubB2ittwlCyRkHt+JEGZXLtuSpsw2OXjCFRTNP4aR5TTSnTYJQJjr9jv4KXXmXOtvA9QM6+svR+wG2bdARxTfKReXv71rXwbLr3rpLv00NE4fVv3wPx33mbm58WpVt6FjQzMtnNTI//xL6jIP28uj2DuI8G9/1qAx043vT8b0gSVSrljFPFHwp8cLRyzCMpe4ZC1LKLuBVI7y+EXhN9HglcMSO7nu/IP24V2hcMbIauu3g1Dk0tGSYMiVN/0CFvq4igR9ntaofz7SVJDLwJYapLPq0o7ygLVOGapdv+48TePfvnkvqmFSrdkAlxOiaxuU3vUhz1qar4HL5TS8ShJJrLzh4u8dx3oJG/rS0myywuletVL5/y4u89Yx53L+sk6aMzV0PrOHt5y1O3ED1kUqjwTHpKrj85/Ht/OyZDs6ar6rQHTU9xz31yh24ddMAhyyeworOPOvX9yf9Vr1ygGHp9Hb0AvtOfZUDEY9+7pX89GlVvurotnp0DZ5LLSS3Exzz7t89t43xsa8hzoJ3B7px8z0MdGzBqzSTrk+hG4K+rcVk1T9RCCOX7r6K/YL0A7eUdJ0yUlkMJ5u0EWxozZCN/N0A9XV20nHKqwT4rpokgiBECKHKEWsCGUpcP0TXBHbkHzzruw9xx/tVzfB5LVmyKYNS9PmSJnD9wQlnU2+JwoBLfc5m9pQMMxodpkRjGg1vWKRq2t++SmVnvvc1i5mVc5jbmObxdb1ceJZyEZW8gKPa6jmyddvkp/b6FOUgxA+gr+LzufPURPPRXz/Fo49uwDB1ejZsxC32JaUVath3sHiKcmEsaLTRhGAzHnXW6K6aP77YxRsPGlprfzjhT7/ga/tcaQ03Ku8dRl3orLomUmmLVMbEsg0Gukt4nruNMbgrGA/pT+ZCB/sF6Ye+S+C7TJm3mPb5TbQ1pNC1WEYZJLr4+IeKg6AVPyQf+ePLJS+y8nUsQ1nvcZA1/nx1cOjTp8/hl8920KMpdUxQkJRcta+iG7B5bR+BH7JlbS8vBSFO1uaxL2yzWhsTj67uoe0wm7wb8LL2HBU/RBdw10ud2+1Leuaceh7ckOfk9ixLOuH2pao+++zpdRT6y8hQNeCIm2Mf95m7efRzr9zhsdWwd3Byu1Kn9Ear1DpLZ1NeXXvzG7ZNsHpiSxFjlCS7GNWEv6+Vz5BhQBAqMYLlGNQ3OsyfVsfs6XX8835JZQd6FIwFLwiHSK9Hwq66d3Yn9gvSr5s2i9kHt/C2E2dR9AI2dJfoiBtBG4N+fCuqHplzTOyI2AcilUxf0aUr79KQNnEsdVqyttLY90VJS7Zt8O9/WMIPLzwUgLcf3sr/PrQ+mUziz+maQGgw0N1PXVP9TpH9wZE1N+fkOXQUXJ7f1M85B7Wytehy2wsdfO/80X24J7dn6fnfjzH1LZ/Dj/YVB4LXbx7aznA8hP+Jvy3fL6ov7k9oiKx7o3s1uWwLv1xW5IN/WMtfLjmGXy1RJHfLM5uS+M+OYF8h/IYZcwDo37wRw8nS0JqhvtHh4PYczVkLpyXDvAuy3PqgaiG6aZR9jRehpGbp722876KX0Zy2cEwdL1AWei7yx9eljGQCcEydhrTJ9JxS5HihTAi921aKlqxtMK81gx5l2HbnXR58VtVM6to0wPIghIj0AVZ2FjhmdiNb8xVe3NQPwLrOAr4b7pJvPC5+5Uu4b003lqHx9buWsWlDP/d/8hXj2kfjh77OCcDS6Fw835nnX4+fxaaBCr/UxQ4t50+Z3zz2RjXsFfhNcwD41Fdu4BXnHc+mgp8YM3f+8ibSzdO54V8v2Ysj3H3o27gOgHRzGzMXT2HO1DqytsHRsxswNQ3b0Ogve1x0pjJYnvjirn9nMB73zq5/zW7DfkH6R7bVkzI0crZJd8ljStqk6IUcNCXDsq4CTVmLUpQk9fLZTRza4tBd8in4IR15jxMaKvxmtcPmgTL1KZNDWrJsyld4bE0PlqFx5jHtvOuYGfzXzUtoztr87oUu3nywIsFlG/u54GVtPL6mh47+Cjf929Gcc+0jO2Xdj4SOos+SDf1YusrqHS/hV2NxnST8+/+xdv6byNkGh7dmOG1OI//O+MtUvH5h406MvoY9ibi7FcB7j25T/2+5kvqTL2fBpQ52XQPFrk1cfPEr+cwr5+6tYU4Y6k++nPbjXsNdXz6b3z67mVkNDq1Zm76yx5yGNBlLI21q5N2Q3rLHQc0OH5+A7625dyYBTll/O0GXWrjNKuRxzvwXOtMzADi8JYXwSuRFClsXOFteoHLzbfD8SrSOHqYBlZMO421Hvxy5YDpauRuvJceqTJaj2+rImBrpP/wPbIRfOhb3fuR3NB/UxIr/+zP3relO9PhnzK5PxlOtzd9VTM8YvGxGjtue28xLS0ZMyhsbQqCd/nbO1JWM9cxvP8QjN/yCE972DhYcuW13pBr2L1S7avrdkPoxgr77Cvof/B5PnHM2my/8Nh/86JvQjnsXet9G8vf8gi2PvUi6tRGntYH2V0YVN/Oj72+8GJd6Z/Jy/v5B+k996UcYKYOpR8+h0jvAYx+7gKMvPASA+lccib3oKBqcDOFAL+v+eCMP/uRRVhU8miydY06cTs+q+5ixdgvp1kZ6lq5j1uUfZoFpI60Mod5AOFVZTU9cdQNP9JRY/Fwnb9S6mPWyadyyvIfzFjTyib8tZ33PYEmE//uXl03Y8b378Ga+/8clo5Z8Hg1SH5qReecHTqT11r/xwn0Pj6vrWIxTvvx3HrjitJ0aQw2TA/nP/juNF72F4GVnj/h+/cmXM/fU1/P0V8/epuXoZMTRf7sNAGvD08iu1UgrTbGjh1t+/RxL8y5H5FKct6mbNfcunbDvHJdPfxKz/n5B+jfcs5pAwpfefAq//uANbC773HPtYwCc8rslnHTRMwhdw0w7rLr9Re6Lsl3XlTyevm0lAOe/2MWMk9pZde9amg7+I2Z9PX4+j9A1Vv5FJWPd+8RmNpZ9LE0QPP43Xvj273jtxy+GuefzpiOm8427lwGMufTbUUjN2CV1zefuXgUwZEm/M31ta4S/72P6l39MAMjbrkWcfdk27+e/ew7Z999M/ck3M/fU1/Pcpw7nf+5fyy/+9DwAax78c7LtZJoM3PYj6Pj8f6BZBmvvfYmlUZ2p5/rLbLz6PqbaE0d1QSjxa4HcvYuZjkkpkDz/i7+zNO+S90OsSGb5aE+ZVT9+grNPncnTj21ic3nkpiA3relj2pYCm8s+f730lwBceu58Ghe2UImanqyLgr6rix7e1i2kmx16H3uUnFvmhBPfwGHtqnfpp0+fs9uO9av3r93hlnYT6b+1NjzN74vq+2t+/n0XIxH+lmJA26Gn8+Mfn8wll/w/lvz3wYjudVx+4vFko0qvn13yEOVI/jjZZJ2t/+/7vPAv57P6+cF+Ii22weayT2dl4nT64/PpT9jXTTj2C9L/j46nk8c/TR+MLuCMFpVF27qgiV88sI6f3rkq2ebdZ8+je1kPfQMVDj1LEeLApjxNC5q5+odPJNv96NYVcOsKFmW3LVj1P//+C+akTeY+sZmjL0/RcPAm3nHM7g+O7ckepiPd1G77Ebx+j42ghj2JhWd+gPbjXsPfm2+h7zMvY+Cmn5I57Ehy9au45UmVBFXu6+TM917KnT/4ETD5iP/gX9/EijmqY+Dbjp9OutnhhUc2Mnu+MlB+8MiaXf6OQI6jyuYkNvX3C9Kvxv8WXxjyfNm7L+TfX7uQH/5lWfJay2HtvOwDb8ZadBSV5x8BICz0U+rsAZ5gOOKlYjW63YCptkG2LUtQdhF+hdl1k6/zUMvZV9J52xeS5y//6j+SIm+jYdbbriE7dc5uHFkNkw1f+t8r+OiMbrb8Ns2z1/6ZwA04av5CdK/Cxacqg6av66384eKjsc77MNdtqOOD7/t84vuvxt6cCF67+jFeW/V8SN5xetfbQY7Ppz95sd+R/nAs/sAlsPAEvm2m+Fb7yawouATlCt7pFxMKgdem6hUt6SzSkDJ5zf/cxl83jy/M/8JAhaZnOjjivdNwpx++Ow9jp9F52xe47EY1EV57wcF8+PWHjPmZ07/xAGecfyo3fuva3T28GiYRlm3Jc0P9dN6waC6ljh7shjpKK5ax6fZ7ecs5Kqb0hk9chJQqzvTKeU2c/M538fiNfyY7dQ5dy7c1mPZH7OuSzf1DuzUKruNoVoX1rKpYvO8Hb+fSc+ez9r7V+D/7HLpf5qp7V3HVvat483//ge/cv4qev9zKJz9+Gjlz9FPztuOnc+FBzUxrSaO9/kN75mB2EtdecDDXXnAwv3y2Y5v6K6BWAwD3r89z//o8l756IUU3oGneDhfwq2EfxeU3vchFR7fjmDr/OOydtJ12HJplMLB2C+WuPvqeeoq+p57CeeFuQinJ52ZT8SVTGx3mnHhG0rgkxkjW//4CJdkMR/2bxJy//1v6bzmklYZwAOGWCF/3ft50z2I2brqNC5w3wh+W8utXaNzQN43OFx/ix196iDf/7ivY//19Pjbt01z5kRsBOK4xxfR6m5vWqCJopzQ7HPOh16C/4aN789B2GG8/vJX6ky9n3ivO56mrXp28Hrt/FvzoY/Qu38zMa37HU+v6+N0v3wOo4G3hwb9hvnn0vsxL3vJ6Dv3tzbvvAGrYLQhv/l/+9POQb77uIyx5zdnMfuXByI9/i0+XX+CbrzuIHJBd8zDCSvFc+mA+8INHKPRX0HWNrk0DbF321Ij7rT/5coSmc/olF3PTvx29R49pd6JWZXOSo97S8P/yc7z+ftKvu5SNjytd743fupa6tvm0/BZOeP0ZnPzOd3HEnCZ6Sh5/X9PLoRd9ls8HIT1L17Hy9hUAvOXwVh5f2cvSvLvPEX41Vt53E/Un37SN3/XH37ofgPcaGv9xypzk9Y2/+DH1c9rG3G//elWGQiv1oBe6APAbZrLVN2lO7feLyn0S/W4I53yAr83v4vXXPswdV30QFp/EL5d08IuvfIc//Gg6W6//F2R9C6WmeXzuhmfI95bJNafZsraPvo2rKHSu22a/9TMWcfcP3sfn/vYif/7OD6j/4Y/V902ioO/Owg3CpOXp9jCZJ4UD4k7UXvufCE0jfP4BXrrjW0w7Qi1HU7kWMq2zeO6hFXzh/MN4Zk0P7fUp3myvZNFLf6HSm6d7eRf/v707D46rvhI9/j33tlrdrc2yZMu2vBvb2CxeABMgmTgJCYSBR0JBHklmhsyDCnmPDI96vJqQ8CpJVSpTTGVIyOTlMcMECGFICBMgMCQFAYPZTGxsA8YLi/dNthZrbam3e8/747aMMNZiW1Jv51OlUvfV7avfT1d9+tfn/vr83mrp5a2WXnpaeom6wpKayPC/NE8N9aS7vX0zt7dvpj7qMrfmg4vS31v0Pwhf++1hj33BmtUA+NFakuueIbnuGdyd6/DvvJmyw+8eXfDC5IeHpy+nJ+3Tk/a5rq6Vve+2ctavXZq0mgdf2Q1Ab9tBbngpyb0HKnno7WYOt/Sy6cnf8vL99/PeqseOG/ABZi1dyjutcZbMmPDBtguvoObjN/PV32ziq7/ZRN1nvzMOvRx9vurR+juDfeUzyeepRSNxzvLl+uqrrw67X9mutXTMWEEiozT07Uc6mnil/Ey++t1g6clPfX4J295v47WbFtJ6748AOLT2fQ69eZiu3mB+fp+nNFSF+dS2P49dh4pM6x03s/PpLYQiId76x18xv66Cj08P8r+y9jH0/Kty3MLStvPGLwEw7+5fM/0r/0Ln3m2cdfmXePupRz603xmXXQNAR0ucA6//cdjjVk2dRyYRx8+kcMMR3HCw7oO4LgsvCsqUeBmfd1av5v67buLz8yaMYq8GF43FNqjquSf7eBG5+6LbfvGN2jlDT4hY97Nbad782lRVPXSyv2usFH16p583+TR+t7WF1p4k1587nUnpBItrY9z73SsAaE+kueGC2TzVnObyG75F31O/wAnvZMrSBuZODoJU9/5OuvZ3s/rMC1m5eU0uuzOq3M6DTLvxMb79v4MAPFid/pOR7OimenoV657fw/uHe/hvZ9Yenc5mAT/35v5rENx/sHo3c89Zyht7t30k4EdqJlFRHby7bTvQ+pFjXPg31+FlfGY0BM+T3/3kbrqbgpRoWUUNNQ1z6G7aQaxuGtVTpnN4bwcA6USwPOfl8bV4HL8sRD7qH+kPJZ+H0iUT9GXHeq5e/Dme3dlOczxDpmoe7zXHmZ79Zz5jUoz6qEu4aQuZ7XsIN85i8rI2KqbUEVkQzPTte2cTu5/ZQPvOjhz2ZPR5NdPGbAWtxjvuB6D6O9fzN5+bl9dPhlL2f3Q1D7wfpe605VRMmkZVbZTpsyZw+ZJpxMpc7n0xCOKH3nrhQ49b/fiPmFsboSPhEazPDY/+1D26UtX0ZZ8kUhFm4fkLAPj0GQ3817OD60PrD3axdMrlxJ+4g8ggtYDyUTozgimbeZziKZmg7519CdVdh/jc3AYqJM2BhNKZyDAzm58/HM/gOkID4NZOJrP3PWKTa4mcvgS3JpjmWHnBxSw+ewVbf3R3DntSmKb9w725boIZyqe+xvPnZpgZSbM7EeKNpm6WTqkinvLZ2tLDfzk3qFqbTv4VGx59hKpp86hpnMdbh7pZXB9lUizE7o7gQ4wNZ/0Fk2ZMon5qFWc01uD5yt+umEEk5DCrJkx4Z5AenbT+OfxUhsjXvpezbp8UVTSPg/pwSiboA/ixWmKuEmrbS3nlXM6ZVkVNebCyVmNVGaHOJpKbXkaTCZJHOglVRBDHpfmJ/wDAcR1qL/oEsy5elstuFKRtX76SRb95ItfNMEP53t9y8Af3MyfcR+P8iYReepDw3DNomL2Et5uDBch/cPVZVH1lGbf9fjNHDvfw29f2cMlpdVSWOdzy8BsA1DXWI46wfFYtjROjlDkO1eUuIUdw3/gDXW++DkBfWycTb/1xzrp7snx/+JF8Pl8qLZmg72z8T3TJJYTa9yKpPianmvEr63Hih4MdVJFMArduKup5lIea6D3QRO+Lqzi4JijLmmhPcN4nVlKxZAVDlW8qtpz/cMpatwMgyTiHHn7guE9kLzV6Ba/M2Jjx4wcJbX2eA7M/yfS27fjLPgvdh5mybw2TZwQpzmSkGleEr318Dmt3HeH8OROJhhzSPtzw6WB1qp2tcaJhl4WTKplTG8V1BFeEMkfwuzs48s5u4IPUX6HxfR/fG6aSbh5H/ZIJ+qGpc9HWXZDoxutsI95wJmHHwauaQnnzu8FJUh8/3oWEI5BJk2jrItHWSc2sWpyyEMmOHg48+nvath7kzEc/yEHqqvuQ7ELsLWvW07hi+DnthSqd/V+OxFvQUBi8DGTSEG8nsWUdqe7eozNC+i8SAsy5bDk7b/zSh7aZ/JNZ/GkagOjFPyHxH1+n7/XngufD9k2UzVxA51P/SaSummsv+CRfWVhB+uDrSEcZrLiSaxdWczDh8OXeV+Hsi8FLEerYTmbCNCTVh9vZTiZawYyrv8D+xwr3Q3zq6/Aj/XFqy8komaCf3PQyzoTJONEK0ru3UTF7CX5lffBDVbSzGZ12OuHTzkYTcXq2vk1fczteOkMoEsZPZ6ia2UD9bf/MsSE9vnMnPQda6dwVrGy18IHHxrdz4yS8byNuTTCzR7wUTk8rJOOQSdG9ZhUdOw7Qs7+Vrv3dtL53hO1zz+NIymPxeVOpnl6NuMKmL/4lZz/+hxz3xAyna83PSQFrr7qZC//pRsKzTwdg15+2Un96HRWzZ+FOakRTCfp2vE/N7MW8E1sAKDL/PLyyKKHediD43Iab6qP7hd8jjsOWB18EYPTmiI2v7PhwmJ3GpSknpWSC/pG33qF6TideIkXrpu1MONzKhGu+AV4KDZWhUxegZRE0E8zJL6uIIq5DbEItXjqowV9/2z8f99heIoV6Pn6RpzBS2zcRaoyjiV4U8FMJNNlHYs8Odv1xPb6nqOeze2vr0frlKV/Z8OeDLDqtjwmzqnGy6xyYwnCoqQft7SK1fRMAU89pRD0fJ1aFM3Mx7rTTcGu20TV5MX3twYXceO1URIU9bgOzd79G8oXfE29qIxQJ0/zGdjJ9x1/TolCoDl86OY9jfukE/bq/v4vWO26mfEIlTlmI3qYj1OzfhjQuBEDDMfxwBU5VPd77G0l1x3Ej5cSmTiTZPnTVTS+V4ci7h+huGqVFOPPEkTv/F6nuXspiEeouWEF4wTJwQqRbDgDgtbfgpxJs/L/PkmhP4IRd/JTH3t40qezbX5+g/niiPUGXK4QiIdat/Ax7drRzzb7SqMpYyP7yjqsITZ2DU1EFQOOlPhrvwjlrJX55BZLoRk6/kLDrMK82WHcilupAy2LM2/cqRzZuwAmHSHUHq9XN/cIncS7/u5z1ZzT4no833Op4p5jTF5FrgO8Di4AVqrp+kP0uBX4KuMAvVPWO4Y5dMkEfoOdAC5lECnFdeg60En/jz5Q37QYgfNrZuFV10N1GpmkX4jhUzZyMOA5e4qP19AdSzydaGymqYmP+k3cRm1xL+YRKkh09+N3tZFoOIKEBawY4Dn4qQ6w+xp69XXiqHEp4tKc9asuCWVE1IQdXoKUzQbQnRVVNOeXV5dSGS+pfr2Bl+pLgOGgy+CBV6nATibZO6tv3I+EoqI+6YcKZJGEveJ54E4LpnYl3NtDb0sHUz67kvUfvJ1YfLcjZOsdSf/gpm6NwHXczcBXwr4PtICIu8HPgs8B+4HUReVJVtw514DF95onIDOBXwBSCQd89qvpTEZkI/BaYDewGvqSq7dnHfBu4HvCAm1X1mdFqz75X95Hs3EHa96msjeBGyiivDUatk5MJ3NpJEAry9xCkeLY+9DJuOAhgkwc57p7ntxGuyL8FVE7W9huupmLqRCob6+neexgnXEb3jj1UZdJIeeRoAEj3xOltbidWH2X6lErW7uog7vmkfCU94L/eU2hJeoQdn4gbXPAORUM8PH051+630X4+23TvS5w/teHofT+Vwfd8klvWQiaNWzcFiVYgoTBkJzOUJeMgDmUfu5Rps3fjL7+C5U9fm6sujDpfdczr5avqNgCRIdOhK4Dtqrozu+/DwJVA7oI+kAFuVdWNIlIFbBCRZ4GvAatU9Q4RuQ24DfiWiCwGrgXOILjO85yILFDVUUmWf2LjK9xZdzaN0RD1EyKUVUSomDIRCPLTTiaNOC6ZeAIvkWLHUxtY9oenhz3u8qdH7XUpL5z2i9+x+swLmTi/llhdlPLqKI7jIM4hInU1qB+8tU3H+8jEE2T6MkRrIyxsDXOgL0ND+QfHirrCkZSHK0LYEdK+j/RlqGioYFLK46ns0naX7z7uu1eTY1t2dnBeopdMPDvS747jloXw4j14iRShTJpQ3RQkVIafCFI4mknjzFyMk+jGX35FLps/JtRT/GHSO+NU06wRGFjxbj9w/nAPGtOgr6pNQFP2dreIbCNo6JXAyuxuDwCrgW9ltz+sqklgl4hsJ3g1e2202jSnooy62gix+hjJjp6jqZtEWxcNlyzAa28mHe8jHU8w69NnjNavLTh98RSJ9gRlkRDl1cFFbQgCfT/1/OBvlQhGf9GwSyP9a4gG+0ysLmcmcKgjQdSVINef9oimPBzXLurms1tii/jv1yziyOYP1pfubWknHU9RNXMSjuPgpTNE6g7juA5uJMjpV0Vi6MH3STbtRvJ0RblT4evwUzazV3L/SkQGvsW5R1Xv6b8jIs8RZEGOdbuqjuSTjMd7Ag37ajNuiVURmQ0sA9YCDdkXBFS1SUT6MyeNwMASlvuz20bNVXs3sPrMC0l2JQlFQ7S+E9R9nzC7l5rtW4kfaiPZ0UMm3oeXylAzmr+8gNRMrsBLeaQTmex1kOBJ7bguiY5uANJdvYgbpL6SXSkcV4i4Dr4q0QnBcD8yIShzMT0SOvoOAcD3FC/lU1UT7Hdn3dnc2rZpPLtoRqB9ZwdlkRBuJEhf9rX2kuhK0n2wh1h9lFAkRCbex4QFMwhFgnPptbeQ2LMDP5UhmsvGjxXVEczeUYB/V9V/GvwwevEptmQ/MGPA/enAweEeNC5BX0QqgUeBW1S1a4g81YheuUTk68DXAWbMmPGRBwxn5eY1PNS4jMb9wb9kxeQY4gjNG96hp6mTd1/ex5V7NpzwcYvJhWtf5GeTlrDjzcNEXeG6K5NUTe/GT2Wvd1RXoJ7P1B/cw1Rg3crP4LhCOhGkelI9wdTXeHMvXspDPWXK0smUVYZJdiVp3tSCG3bwUx6+pxbw89B5tRH8lEfru8HAKFYXJVYfw/eUA+ubOLQpw5wpFdTOmUCoIkrH9uDT7aFoiPn3PZrLpo8pz1O8TF6UYXgdmC8ic4ADBKnxrwz3oDEP+iJSRhDwH1LV/k8tHRaRqdlR/lSgObt9RK9c2bdI90BQT/9k2nUk5dPX2suc6nK69ncTiobobetDPS35gN/vrIYKdu9KsaCynL1rDlA3r5dznv3TcfddsXoVECyZ6KW9o/dfWPQxANK+T29bHxWug1vmEqmN4LjChWtfHJ/OmBMWdgRxBd9TyiIhkl0pWra2Eq2PEY6VUZnxSfWk6W3ro6orTixbgnzWTx/KccvHlvoMX3DtFIO+iHwR+BkwCfiDiLypqpeIyDSCqZmXqWpGRL4JPEMwZfM+Vd0y3LHHevaOAPcC21R14FytJ4HrgDuy358YsP3XIvJjggu584F1Y9G2v2t5aywOW1QOdSS4M77thB5z7LTV6unVtO1oxxUh1ZOmvDqDuEJ5dZieIvtcQ7E5kvKpPdyLG3ZIZmeweSmfRHsCcYXaijK8lM+K1avYeOklRTehYVDjUGVTVR8HHj/O9oPAZQPu/xEYflWbAcZ6pH8R8NfA2yLyZnbbdwiC/SMicj2wF7gGQFW3iMgjBFOOMsBNozVzx5y40ZhOec6zf+Kp2edSWRmmvDqMl/Jwwy6x+hjnPf/cKLTSjIW/r1zMvIowad/HTyqHOpMA9GR8vK4kM2MhWpLBrCwovhlsQxnJlM18XpFwrGfvvMLx8/QAnxnkMT8EfjhmjTLjrs/zifRl6NrfjZfxmbasgd7W3lH9HbfEFnFX74m9KzGDS/nKvr7gk9WeKpWhYPaWK8LBRJqDiTSeQkN56X3Izvf8Yads5nMdhtI7Y2bcWbmFwtST8XElQ9gRdmfXie5Xyi+wOqJ6+vkb9S3oG2M+5JbYoqO3O9PDlZMsPTqCKZv5zMl1A0xpWX3mhRz6/jdG7Xj9AaqUR56jbai/5V2920r+b91fcG2or3x+UbCRvhlXo7mi2A9rz6SmzMYtY2FBZZj3eoYuNFiqVP2jC78PvtP4tOVkWNA3Bev29s25bkJRGpjeAXsXdSz1hw/6msdR34ZJxpghHfsiUOrU94b9sjVyjTEF4XgB3kb6H6aZNH56mDU2hl1PMXdspG+MOerYAG8B/6P6c/qFOtK3oG+MGZSldj6qP6c/ZNC3nL4xplDc1buNj02M2ih/EOql8TOpIb/GujbPqbCcvjHmI/rrLlng/yhfffzhpmzm8Ujfgr4xxpyAD1I4Q+yTxzl9C/rGGHMCNJMZdvYOeTx7x4K+McacABvpG2NMCRlRGQbL6RtjTHEYyUg/n+fpW9A3xpgToJ6Hn0kPvY8FfWOMKQ6qNtI3xpiS4fs+UsBVNi3oG2PMCRhJwTV8m7JpjDFFYST19POZBX1jjDkBI8np53NpZQv6xhgzchv8jr1ItG7QHTSTRFO9AG3j1qoTYFU2jTFm5H7ld+5BM8lBd/BbtuJOOh1VHXpeZ45Y0DfGmBFS1ZQ7aRF+y9bj/zyTxO/ci7dvTfk4N23ELOgbY8wJ8PatKR9stD9glD/M9J7csaBvjDEnYLDRfiGM8sGCvjHGnLDjjfYLYZQPFvSNMeaEHTvaL5RRPljQN8aYkzJwtF8oo3ywefrGGHNSVDUVmnkR/qE38HsOQ+JI3o/ywUb6xhhz0oLR/r6CGeUDSD7XfR4JEWkB4kBrrtsyhuop3v4Vc9/A+pdvZqnqpNE8oIhUAn2qWhAFeQo+6AOIyHpVPTfX7Rgrxdy/Yu4bWP9M/rH0jjHGlBAL+sYYU0KKJejfk+sGjLFi7l8x9w2sfybPFEVO3xhjzMgUy0jfGGPMCBR00BeRS0XkXRHZLiK35bo9o0FEdovI2yLypoisz26bKCLPisj72e+1uW7nSInIfSLSLCKbB2wbtD8i8u3s+XxXRC7JTatHbpD+fV9EDmTP4ZsictmAnxVM/0Rkhoi8ICLbRGSLiPzP7PaiOX+lqGCDvoi4wM+BzwOLgS+LyOLctmrUfEpVlw6YCncbsEpV5wOrsvcLxS+BS4/Zdtz+ZM/ftcAZ2cf8v+x5zme/5KP9A/hJ9hwuVdU/QkH2LwPcqqqLgI8BN2X7UEznr+QUbNAHVgDbVXVn9pNwDwNX5rhNY+VK4IHs7QeAL+SuKSdGVV8CjhyzebD+XAk8rKpJVd0FbCc4z3lrkP4NpqD6p6pNqroxe7sb2AY0UkTnrxQVctBvBPYNuL8/u63QKfAnEdkgIl/PbmtQ1SYInojA5Jy1bnQM1p9iOqffFJFN2fRPf/qjYPsnIrOBZcBaSuP8Fa1CDvpynG3FMBXpIlVdTpC2uklE/iLXDRpHxXJO7wbmAUuBJuDO7PaC7F+2zMCjwC2q2jXUrsfZlvf9KzWFHPT3AzMG3J8OHMxRW0aNqh7Mfm8GHid4e3xYRKYCZL83566Fo2Kw/hTFOVXVw6rqqaoP/BsfpDgKrn8iUkYQ8B9S1ceym4v6/BW7Qg76rwPzRWSOiIQJLiA9meM2nRIRqRCRqv7bwOeAzQT9ui6723XAE7lp4agZrD9PAteKSLmIzAHmA+ty0L5T0h8Qs75IcA6hwPonIgLcC2xT1R8P+FFRn79iV7D19FU1IyLfBJ4BXOA+Vd2S42adqgbg8eC5Rgj4tao+LSKvA4+IyPXAXuCaHLbxhIjIb4CVQL2I7Ae+B9zBcfqjqltE5BFgK8HMkZvyvXLhIP1bKSJLCVIbu4EboSD7dxHw18DbIvJmdtt3KKLzV4rsE7nGGFNCCjm9Y4wx5gRZ0DfGmBJiQd8YY0qIBX1jjCkhFvSNMaaEWNA3xpgSUrDz9E1pEpHvAz1ANfCSqj53isd7mqCC5Cuqevmpt9CY/GZB3xQkVf3uKB3qR0CM7AeojCl2lt4xeU9Ebs8uyvEcsDC77ZcicnX29m4R+QcReU1E1ovIchF5RkR2iMg3hjq2qq4Cuse+F8bkBxvpm7wmIucQ1FVaRvD/uhHYcJxd96nqBSLyE4KFTS4CIsAW4F/Gp7XG5D8L+ibffQJ4XFV7AURksKJ6/dvfBiqzi350i0hCRCaoasfYN9WY/GfpHVMIRlIgKpn97g+43X/fBjfGZFnQN/nuJeCLIhLNlp2+ItcNMqaQ2QjI5DVV3SgivwXeBPYAL4/m8UXkZeB0oDJbGvl6VX1mNH+HMfnESisbY0wJsfSOMcaUEEvvmKInImcBDx6zOamq5+eiPcbkkqV3jDGmhFh6xxhjSogFfWOMKSEW9I0xpoRY0DfGmBJiQd8YY0rI/wcBN+uolQel6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# short look into training data: large biases\n",
    "# any problem from normalizing?\n",
    "i=4\n",
    "xr.DataArray(np.vstack([X[i],y[i]])).plot(yincrease=False, robust=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `fit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:AutoGraph could not transform <bound method PeriodicPadding2D.call of <WeatherBench.src.train_nn.PeriodicPadding2D object at 0x2b5766ccf450>> and will run it as-is.\n",
      "Please report this to the TensorFlow team. When filing the bug, set the verbosity to 10 (on Linux, `export AUTOGRAPH_VERBOSITY=10`) and attach the full output.\n",
      "Cause: module 'gast' has no attribute 'Index'\n",
      "To silence this warning, decorate the function with @tf.autograph.experimental.do_not_convert\n",
      "WARNING: AutoGraph could not transform <bound method PeriodicPadding2D.call of <WeatherBench.src.train_nn.PeriodicPadding2D object at 0x2b5766ccf450>> and will run it as-is.\n",
      "Please report this to the TensorFlow team. When filing the bug, set the verbosity to 10 (on Linux, `export AUTOGRAPH_VERBOSITY=10`) and attach the full output.\n",
      "Cause: module 'gast' has no attribute 'Index'\n",
      "To silence this warning, decorate the function with @tf.autograph.experimental.do_not_convert\n"
     ]
    }
   ],
   "source": [
    "cnn = keras.models.Sequential([\n",
    "    PeriodicConv2D(filters=32, kernel_size=5, conv_kwargs={'activation':'relu'}, input_shape=(32, 64, 1)),\n",
    "    PeriodicConv2D(filters=1, kernel_size=5)\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"sequential\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "periodic_conv2d (PeriodicCon (None, 32, 64, 32)        832       \n",
      "_________________________________________________________________\n",
      "periodic_conv2d_1 (PeriodicC (None, 32, 64, 1)         801       \n",
      "=================================================================\n",
      "Total params: 1,633\n",
      "Trainable params: 1,633\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "cnn.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnn.compile(keras.optimizers.Adam(1e-4), 'mse')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/3\n",
      "30/30 [==============================] - 16s 505ms/step - loss: 0.1810 - val_loss: 0.0963\n",
      "Epoch 2/3\n",
      "30/30 [==============================] - 15s 495ms/step - loss: 0.0917 - val_loss: 0.0657\n",
      "Epoch 3/3\n",
      "30/30 [==============================] - 14s 462ms/step - loss: 0.0663 - val_loss: 0.0601\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<tensorflow.python.keras.callbacks.History at 0x2b55b1642e90>"
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cnn.fit(dg_train, epochs=3, validation_data=dg_valid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `predict`\n",
    "\n",
    "Create predictions and print `mean(variable, lead_time, longitude, weighted latitude)` RPSS for all years as calculated by `skill_by_year`. For now RPS, todo: change to RPSS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scripts import add_valid_time_from_forecast_reference_time_and_lead_time\n",
    "\n",
    "def _create_predictions(model, dg, lead):\n",
    "    \"\"\"Create non-iterative predictions\"\"\"\n",
    "    preds = model.predict(dg).squeeze()\n",
    "    # Unnormalize\n",
    "    preds = preds * dg.fct_std.values + dg.fct_mean.values\n",
    "    if dg.verif_dataset:\n",
    "        da = xr.DataArray(\n",
    "                    preds,\n",
    "                    dims=['forecast_time', 'latitude', 'longitude','variable'],\n",
    "                    coords={'forecast_time': dg.fct_data.forecast_time, 'latitude': dg.fct_data.latitude,\n",
    "                            'longitude': dg.fct_data.longitude},\n",
    "                ).to_dataset() # doesnt work yet\n",
    "    else:\n",
    "        da = xr.DataArray(\n",
    "                    preds,\n",
    "                    dims=['forecast_time', 'latitude', 'longitude'],\n",
    "                    coords={'forecast_time': dg.fct_data.forecast_time, 'latitude': dg.fct_data.latitude,\n",
    "                            'longitude': dg.fct_data.longitude},\n",
    "                )\n",
    "    da = da.assign_coords(lead_time=lead)\n",
    "    # da = add_valid_time_from_forecast_reference_time_and_lead_time(da)\n",
    "    return da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optionally masking the ocean when making probabilistic\n",
    "mask = obs_2020.std(['lead_time','forecast_time']).notnull()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scripts import make_probabilistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!renku storage pull ../data/hindcast-like-observations_2000-2019_biweekly_tercile-edges.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "cache_path='../data'\n",
    "tercile_file = f'{cache_path}/hindcast-like-observations_2000-2019_biweekly_tercile-edges.nc'\n",
    "tercile_edges = xr.open_dataset(tercile_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is not useful but results have expected dimensions\n",
    "# actually train for each lead_time\n",
    "\n",
    "def create_predictions(cnn, fct, obs, time):\n",
    "    preds_test=[]\n",
    "    for lead in fct.lead_time:\n",
    "        dg = DataGenerator(fct.mean('realization').sel(forecast_time=time)[v],\n",
    "                           obs.sel(forecast_time=time)[v],\n",
    "                           lead_time=lead, batch_size=bs, mean=dg_train.fct_mean, std=dg_train.fct_std, shuffle=False)\n",
    "        preds_test.append(_create_predictions(cnn, dg, lead))\n",
    "    preds_test = xr.concat(preds_test, 'lead_time')\n",
    "    preds_test['lead_time'] = fct.lead_time\n",
    "    # add valid_time coord\n",
    "    preds_test = add_valid_time_from_forecast_reference_time_and_lead_time(preds_test)\n",
    "    preds_test = preds_test.to_dataset(name=v)\n",
    "    # add fake var\n",
    "    preds_test['tp'] = preds_test['t2m']\n",
    "    # make probabilistic\n",
    "    preds_test = make_probabilistic(preds_test.expand_dims('realization'), tercile_edges, mask=mask)\n",
    "    return preds_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `predict` training period in-sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!renku storage pull ../data/forecast-like-observations_2020_biweekly_terciled.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!renku storage pull ../data/hindcast-like-observations_2000-2019_biweekly_terciled.zarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scripts import skill_by_year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#on renku with small memory\n",
    "#step = 3\n",
    "#for year in np.arange(int(time_train_start), int(time_train_end) -1, step): # loop over years to consume less memory on renku\n",
    "#    preds_is = create_predictions(cnn, hind_2000_2019, obs_2000_2019, time=slice(str(year), str(year+step-1))).compute()\n",
    "#    print(skill_by_year(preds_is))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>RPSS</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>year</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2000</th>\n",
       "      <td>-0.634723</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2001</th>\n",
       "      <td>-0.785985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2002</th>\n",
       "      <td>-0.879244</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2003</th>\n",
       "      <td>-0.817812</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2004</th>\n",
       "      <td>-0.847618</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2005</th>\n",
       "      <td>-0.922555</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2006</th>\n",
       "      <td>-0.835650</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2007</th>\n",
       "      <td>-0.930159</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2008</th>\n",
       "      <td>-0.816114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2009</th>\n",
       "      <td>-0.932632</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2010</th>\n",
       "      <td>-0.908160</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2011</th>\n",
       "      <td>-0.788804</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012</th>\n",
       "      <td>-0.888000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013</th>\n",
       "      <td>-0.891895</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014</th>\n",
       "      <td>-0.886370</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2015</th>\n",
       "      <td>-0.922456</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016</th>\n",
       "      <td>-1.004312</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017</th>\n",
       "      <td>-0.965584</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          RPSS\n",
       "year          \n",
       "2000 -0.634723\n",
       "2001 -0.785985\n",
       "2002 -0.879244\n",
       "2003 -0.817812\n",
       "2004 -0.847618\n",
       "2005 -0.922555\n",
       "2006 -0.835650\n",
       "2007 -0.930159\n",
       "2008 -0.816114\n",
       "2009 -0.932632\n",
       "2010 -0.908160\n",
       "2011 -0.788804\n",
       "2012 -0.888000\n",
       "2013 -0.891895\n",
       "2014 -0.886370\n",
       "2015 -0.922456\n",
       "2016 -1.004312\n",
       "2017 -0.965584"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# with larger memory, simply do\n",
    "preds_is = create_predictions(cnn, hind_2000_2019, obs_2000_2019, time=slice(time_train_start, time_train_end))\n",
    "skill_by_year(preds_is)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `predict` validation period out-of-sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,