LSTM_for_rainfall_runoff_modelling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Long Short-Term Memory (LSTM) for rainfall-runoff modelling\n",
    "\n",
    "Recently, Kratzert et al. (2018a, 2018b) have shown the potential of LSTMs for rainfall-runoff modelling. Here, I'll show some example for setting up and training such a model. For single basin calibration (as done in this notebook) no GPU is required, and everything should run fine on a standard CPU (also rather slow, compared to a GPU).\n",
    "\n",
    "In this example, we use the CAMELS data set (Newman et al. 2014) that provides us with approx. 35 years of daily meteorological forcings and discharge observations from 671 basins across the contiguos USA. As input to our model we use daily precipitation sum, daily min/max temperature as well as average solar radiation and vapor pressure. \n",
    "\n",
    "\n",
    "- Kratzert, F., Klotz, D., Brenner, C., Schulz, K., and Herrnegger, M.: Rainfall–runoff modelling using Long Short-Term Memory (LSTM) networks, Hydrol. Earth Syst. Sci., 22, 6005-6022, https://doi.org/10.5194/hess-22-6005-2018, 2018a. \n",
    "\n",
    "- Kratzert F., Klotz D., Herrnegger M., Hochreiter S.: A glimpse into the Unobserved: Runoff simulation for ungauged catchments with LSTMs, Workshop on Modeling and Decision-Making in the Spatiotemporal Domain, 32nd Conference on Neural Information Processing Systems (NeuRIPS 2018), Montréal, Canada. [https://openreview.net/forum?id=Bylhm72oKX](https://openreview.net/forum?id=Bylhm72oKX), 2018b.\n",
    "\n",
    "- A. Newman; K. Sampson; M. P. Clark; A. Bock; R. J. Viger; D. Blodgett, 2014. A large-sample watershed-scale hydrometeorological dataset for the contiguous USA. Boulder, CO: UCAR/NCAR. https://dx.doi.org/10.5065/D6MW2F4D\n",
    "\n",
    "Date: 29.05.2019<br/>\n",
    "Created by: Frederik Kratzert (kratzert@ml.jku.at)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "from pathlib import Path\n",
    "from typing import Tuple, List\n",
    "\n",
    "import gcsfs\n",
    "import matplotlib.pyplot as plt\n",
    "from numba import njit\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from torch.utils.data import DataLoader, Dataset\n",
    "import tqdm\n",
    "\n",
    "# Globals\n",
    "FILE_SYSTEM = gcsfs.core.GCSFileSystem(requester_pays=True)\n",
    "CAMELS_ROOT = Path('pangeo-ncar-camels/basin_dataset_public_v1p2')\n",
    "DEVICE = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\") # This line checks if GPU is available"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data loading\n",
    "### CAMELS specific data loader functions\n",
    "\n",
    "Next we define two functions to load the meteorological forcings and the discharge for any specific basin from the camels data set. From the header of the forcing file we also extract the catchment area, to normalize the discharge (to mm/d)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_forcing(basin: str) -> Tuple[pd.DataFrame, int]:\n",
    "    \"\"\"Load the meteorological forcing data of a specific basin.\n",
    "\n",
    "    :param basin: 8-digit code of basin as string.\n",
    "    \n",
    "    :return: pd.DataFrame containing the meteorological forcing data and the\n",
    "        area of the basin as integer.\n",
    "    \"\"\"\n",
    "    # root directory of meteorological forcings\n",
    "    forcing_path = CAMELS_ROOT / 'basin_mean_forcing' / 'daymet'\n",
    "    \n",
    "    # get path of forcing file\n",
    "    files = list(FILE_SYSTEM.glob(f\"{str(forcing_path)}/**/{basin}_*.txt\"))\n",
    "    if len(files) == 0:\n",
    "        raise RuntimeError(f'No forcing file file found for Basin {basin}')\n",
    "    else:\n",
    "        file_path = files[0]\n",
    "    \n",
    "    # read-in data and convert date to datetime index\n",
    "    with FILE_SYSTEM.open(file_path) as fp:\n",
    "        df = pd.read_csv(fp, sep='\\s+', header=3)\n",
    "    dates = (df.Year.map(str) + \"/\" + df.Mnth.map(str) + \"/\"\n",
    "             + df.Day.map(str))\n",
    "    df.index = pd.to_datetime(dates, format=\"%Y/%m/%d\")\n",
    "\n",
    "    # load area from header\n",
    "    with FILE_SYSTEM.open(file_path) as fp:\n",
    "        content = fp.readlines()\n",
    "        area = int(content[2])\n",
    "\n",
    "    return df, area\n",
    "\n",
    "\n",
    "def load_discharge(basin: str, area: int) ->  pd.Series:\n",
    "    \"\"\"Load the discharge time series for a specific basin.\n",
    "\n",
    "    :param basin: 8-digit code of basin as string.\n",
    "    :param area: int, area of the catchment in square meters\n",
    "    \n",
    "    :return: A pd.Series containng the catchment normalized discharge.\n",
    "    \"\"\"\n",
    "    # root directory of the streamflow data\n",
    "    discharge_path = CAMELS_ROOT / 'usgs_streamflow'\n",
    "    \n",
    "    # get path of streamflow file file\n",
    "    files = list(FILE_SYSTEM.glob(f\"{str(discharge_path)}/**/{basin}_*.txt\"))\n",
    "    if len(files) == 0:\n",
    "        raise RuntimeError(f'No discharge file found for Basin {basin}')\n",
    "    else:\n",
    "        file_path = files[0]\n",
    "\n",
    "    # read-in data and convert date to datetime index\n",
    "    col_names = ['basin', 'Year', 'Mnth', 'Day', 'QObs', 'flag']\n",
    "    with FILE_SYSTEM.open(file_path) as fp:\n",
    "        df = pd.read_csv(fp, sep='\\s+', header=None, names=col_names)\n",
    "    dates = (df.Year.map(str) + \"/\" + df.Mnth.map(str) + \"/\"\n",
    "             + df.Day.map(str))\n",
    "    df.index = pd.to_datetime(dates, format=\"%Y/%m/%d\")\n",
    "\n",
    "    # normalize discharge from cubic feed per second to mm per day\n",
    "    df.QObs = 28316846.592 * df.QObs * 86400 / (area * 10 ** 6)\n",
    "\n",
    "    return df.QObs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data reshaping for LSTM training\n",
    "Next we need another utility function to reshape the data into an appropriate format for training LSTMs. These recurrent neural networks expect sequential input of the shape `(sequence length, number of features)`. We train our network to predict a single day of discharge from *n* days of precendent meteorological observations. For example, lets assume that _n_ = 365, then a single training sample should be of shape `(365, number of features)`, and since we use 5 input features the shape is `(365, 5)`.\n",
    "\n",
    "However, when loaded from the files the entire data is stored in a matrix, where the number of rows correspond to the total number of days in the training set and the number of columns to the features. Thus we need to slide over this matrix and cut out small samples appropriate to our LSTM setting. To speed things up, we make use of the awesome Numba library here (the little @njit decorator JIT-compiles this function and dramatically increases the speed)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "@njit\n",
    "def reshape_data(x: np.ndarray, y: np.ndarray, seq_length: int) -> Tuple[np.ndarray, np.ndarray]:\n",
    "    \"\"\"\n",
    "    Reshape matrix data into sample shape for LSTM training.\n",
    "\n",
    "    :param x: Matrix containing input features column wise and time steps row wise\n",
    "    :param y: Matrix containing the output feature.\n",
    "    :param seq_length: Length of look back days for one day of prediction\n",
    "    \n",
    "    :return: Two np.ndarrays, the first of shape (samples, length of sequence,\n",
    "        number of features), containing the input data for the LSTM. The second\n",
    "        of shape (samples, 1) containing the expected output for each input\n",
    "        sample.\n",
    "    \"\"\"\n",
    "    num_samples, num_features = x.shape\n",
    "\n",
    "    x_new = np.zeros((num_samples - seq_length + 1, seq_length, num_features))\n",
    "    y_new = np.zeros((num_samples - seq_length + 1, 1))\n",
    "\n",
    "    for i in range(0, x_new.shape[0]):\n",
    "        x_new[i, :, :num_features] = x[i:i + seq_length, :]\n",
    "        y_new[i, :] = y[i + seq_length - 1, 0]\n",
    "\n",
    "    return x_new, y_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyTorch data set\n",
    "Now we wrap everything into a PyTorch Dataset. These are specific classes that can be used by PyTorchs DataLoader class for generating mini-batches (and do this in parallel in multiple threads). Such a data set class has to inherit from the PyTorch Dataset class and three functions have to be implemented\n",
    "\n",
    "1. __init__(): The object initializing function.\n",
    "2. __len__(): This function has to return the number of samples in the data set.\n",
    "3. __getitem__(i): A function that returns sample `i` of the data set (sample + target value)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CamelsTXT(Dataset):\n",
    "    \"\"\"Torch Dataset for basic use of data from the CAMELS data set.\n",
    "\n",
    "    This data set provides meteorological observations and discharge of a given\n",
    "    basin from the CAMELS data set.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, basin: str, seq_length: int=365,period: str=None,\n",
    "                 dates: List=None, means: pd.Series=None, stds: pd.Series=None):\n",
    "        \"\"\"Initialize Dataset containing the data of a single basin.\n",
    "\n",
    "        :param basin: 8-digit code of basin as string.\n",
    "        :param seq_length: (optional) Length of the time window of\n",
    "            meteorological input provided for one time step of prediction.\n",
    "        :param period: (optional) One of ['train', 'eval']. None loads the \n",
    "            entire time series.\n",
    "        :param dates: (optional) List of pd.DateTimes of the start and end date \n",
    "            of the discharge period that is used.\n",
    "        :param means: (optional) Means of input and output features derived from\n",
    "            the training period. Has to be provided for 'eval' period. Can be\n",
    "            retrieved if calling .get_means() on the data set.\n",
    "        :param stds: (optional) Stds of input and output features derived from\n",
    "            the training period. Has to be provided for 'eval' period. Can be\n",
    "            retrieved if calling .get_stds() on the data set.\n",
    "        \"\"\"\n",
    "        self.basin = basin\n",
    "        self.seq_length = seq_length\n",
    "        self.period = period\n",
    "        self.dates = dates\n",
    "        self.means = means\n",
    "        self.stds = stds\n",
    "\n",
    "        # load data into memory\n",
    "        self.x, self.y = self._load_data()\n",
    "\n",
    "        # store number of samples as class attribute\n",
    "        self.num_samples = self.x.shape[0]\n",
    "\n",
    "    def __len__(self):\n",
    "        return self.num_samples\n",
    "\n",
    "    def __getitem__(self, idx: int):\n",
    "        return self.x[idx], self.y[idx]\n",
    "\n",
    "    def _load_data(self):\n",
    "        \"\"\"Load input and output data from text files.\"\"\"\n",
    "        df, area = load_forcing(self.basin)\n",
    "        df['QObs(mm/d)'] = load_discharge(self.basin, area)\n",
    "        \n",
    "        if self.dates is not None:\n",
    "            # If meteorological observations exist before start date\n",
    "            # use these as well. Similiar to hydrological warmup period.\n",
    "            if self.dates[0] - pd.DateOffset(days=self.seq_length) > df.index[0]:\n",
    "                start_date = self.dates[0] - pd.DateOffset(days=self.seq_length)\n",
    "            else:\n",
    "                start_date = self.dates[0]\n",
    "            df = df[start_date:self.dates[1]]\n",
    "\n",
    "        # if training period store means and stds\n",
    "        if self.period == 'train':\n",
    "            self.means = df.mean()\n",
    "            self.stds = df.std()\n",
    "\n",
    "        # extract input and output features from DataFrame\n",
    "        x = np.array([df['prcp(mm/day)'].values,\n",
    "                      df['srad(W/m2)'].values,\n",
    "                      df['tmax(C)'].values,\n",
    "                      df['tmin(C)'].values,\n",
    "                      df['vp(Pa)'].values]).T\n",
    "        y = np.array([df['QObs(mm/d)'].values]).T\n",
    "\n",
    "        # normalize data, reshape for LSTM training and remove invalid samples\n",
    "        x = self._local_normalization(x, variable='inputs')\n",
    "        x, y = reshape_data(x, y, self.seq_length)\n",
    "\n",
    "        if self.period == \"train\":\n",
    "            # Delete all samples, where discharge is NaN\n",
    "            if np.sum(np.isnan(y)) > 0:\n",
    "                print(f\"Deleted some records because of NaNs {self.basin}\")\n",
    "                x = np.delete(x, np.argwhere(np.isnan(y)), axis=0)\n",
    "                y = np.delete(y, np.argwhere(np.isnan(y)), axis=0)\n",
    "            \n",
    "            # Deletes all records, where no discharge was measured (-999)\n",
    "            x = np.delete(x, np.argwhere(y < 0)[:, 0], axis=0)\n",
    "            y = np.delete(y, np.argwhere(y < 0)[:, 0], axis=0)\n",
    "            \n",
    "            # normalize discharge\n",
    "            y = self._local_normalization(y, variable='output')\n",
    "\n",
    "        # convert arrays to torch tensors\n",
    "        x = torch.from_numpy(x.astype(np.float32))\n",
    "        y = torch.from_numpy(y.astype(np.float32))\n",
    "\n",
    "        return x, y\n",
    "\n",
    "    def _local_normalization(self, feature: np.ndarray, variable: str) -> \\\n",
    "            np.ndarray:\n",
    "        \"\"\"Normalize input/output features with local mean/std.\n",
    "\n",
    "        :param feature: Numpy array containing the feature(s) as matrix.\n",
    "        :param variable: Either 'inputs' or 'output' showing which feature will\n",
    "            be normalized\n",
    "        :return: array containing the normalized feature\n",
    "        \"\"\"\n",
    "        if variable == 'inputs':\n",
    "            means = np.array([self.means['prcp(mm/day)'],\n",
    "                              self.means['srad(W/m2)'],\n",
    "                              self.means['tmax(C)'],\n",
    "                              self.means['tmin(C)'],\n",
    "                              self.means['vp(Pa)']])\n",
    "            stds = np.array([self.stds['prcp(mm/day)'],\n",
    "                             self.stds['srad(W/m2)'],\n",
    "                             self.stds['tmax(C)'],\n",
    "                             self.stds['tmin(C)'],\n",
    "                             self.stds['vp(Pa)']])\n",
    "            feature = (feature - means) / stds\n",
    "        elif variable == 'output':\n",
    "            feature = ((feature - self.means[\"QObs(mm/d)\"]) /\n",
    "                       self.stds[\"QObs(mm/d)\"])\n",
    "        else:\n",
    "            raise RuntimeError(f\"Unknown variable type {variable}\")\n",
    "\n",
    "        return feature\n",
    "\n",
    "    def local_rescale(self, feature: np.ndarray, variable: str) -> \\\n",
    "            np.ndarray:\n",
    "        \"\"\"Rescale input/output features with local mean/std.\n",
    "\n",
    "        :param feature: Numpy array containing the feature(s) as matrix.\n",
    "        :param variable: Either 'inputs' or 'output' showing which feature will\n",
    "            be normalized\n",
    "        :return: array containing the normalized feature\n",
    "        \"\"\"\n",
    "        if variable == 'inputs':\n",
    "            means = np.array([self.means['prcp(mm/day)'],\n",
    "                              self.means['srad(W/m2)'],\n",
    "                              self.means['tmax(C)'],\n",
    "                              self.means['tmin(C)'],\n",
    "                              self.means['vp(Pa)']])\n",
    "            stds = np.array([self.stds['prcp(mm/day)'],\n",
    "                             self.stds['srad(W/m2)'],\n",
    "                             self.stds['tmax(C)'],\n",
    "                             self.stds['tmin(C)'],\n",
    "                             self.stds['vp(Pa)']])\n",
    "            feature = feature * stds + means\n",
    "        elif variable == 'output':\n",
    "            feature = (feature * self.stds[\"QObs(mm/d)\"] +\n",
    "                       self.means[\"QObs(mm/d)\"])\n",
    "        else:\n",
    "            raise RuntimeError(f\"Unknown variable type {variable}\")\n",
    "\n",
    "        return feature\n",
    "\n",
    "    def get_means(self):\n",
    "        return self.means\n",
    "\n",
    "    def get_stds(self):\n",
    "        return self.stds\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build LSTM model\n",
    "\n",
    "Here we implement a single layer LSTM with optional dropout in the final fully connected layer. Using PyTorch, we need to inherit from the `nn.Module` class and implement the `__init__()`, as well as a `forward()` function. To make things easy, we make use of the standard LSTM layer included in the PyTorch library.\n",
    "\n",
    "The forward function implements the entire forward pass through the network, while the backward pass, used for updating the weights during training, is automatically derived from PyTorch autograd functionality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Model(nn.Module):\n",
    "    \"\"\"Implementation of a single layer LSTM network\"\"\"\n",
    "    \n",
    "    def __init__(self, hidden_size: int, dropout_rate: float=0.0):\n",
    "        \"\"\"Initialize model\n",
    "        \n",
    "        :param hidden_size: Number of hidden units/LSTM cells\n",
    "        :param dropout_rate: Dropout rate of the last fully connected\n",
    "            layer. Default 0.0\n",
    "        \"\"\"\n",
    "        super(Model, self).__init__()\n",
    "        self.hidden_size = hidden_size\n",
    "        self.dropout_rate = dropout_rate\n",
    "        \n",
    "        # create required layer\n",
    "        self.lstm = nn.LSTM(input_size=5, hidden_size=self.hidden_size, \n",
    "                            num_layers=1, bias=True, batch_first=True)\n",
    "        self.dropout = nn.Dropout(p=self.dropout_rate)\n",
    "        self.fc = nn.Linear(in_features=self.hidden_size, out_features=1)\n",
    "        \n",
    "    def forward(self, x: torch.Tensor) -> torch.Tensor:\n",
    "        \"\"\"Forward pass through the Network.\n",
    "        \n",
    "        :param x: Tensor of shape [batch size, seq length, num features]\n",
    "            containing the input data for the LSTM network.\n",
    "        \n",
    "        :return: Tensor containing the network predictions\n",
    "        \"\"\"\n",
    "        output, (h_n, c_n) = self.lstm(x)\n",
    "        \n",
    "        # perform prediction only at the end of the input sequence\n",
    "        pred = self.fc(self.dropout(h_n[-1,:,:]))\n",
    "        return pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train utilities\n",
    "\n",
    "Next we are implementing two functions for training and evaluating the model.\n",
    "\n",
    "- `train_epoch()`: This function iterates one time over the entire training data set (called one epoch) and updates the weights of the network to minimize the loss function (we use the mean squared error here).\n",
    "- `eval_model()`: To this function evaluates a data set and returns the predictions, as well as observations.\n",
    "\n",
    "Furthermore we implement a function `calc_nse()` to calculate the Nash-Sutcliffe-Efficiency for our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_epoch(model, optimizer, loader, loss_func, epoch):\n",
    "    \"\"\"Train model for a single epoch.\n",
    "\n",
    "    :param model: A torch.nn.Module implementing the LSTM model\n",
    "    :param optimizer: One of PyTorchs optimizer classes.\n",
    "    :param loader: A PyTorch DataLoader, providing the trainings\n",
    "        data in mini batches.\n",
    "    :param loss_func: The loss function to minimize.\n",
    "    :param epoch: The current epoch (int) used for the progress bar\n",
    "    \"\"\"\n",
    "    # set model to train mode (important for dropout)\n",
    "    model.train()\n",
    "    pbar = tqdm.tqdm_notebook(loader)\n",
    "    pbar.set_description(f\"Epoch {epoch}\")\n",
    "    # request mini-batch of data from the loader\n",
    "    for xs, ys in pbar:\n",
    "        # delete previously stored gradients from the model\n",
    "        optimizer.zero_grad()\n",
    "        # push data to GPU (if available)\n",
    "        xs, ys = xs.to(DEVICE), ys.to(DEVICE)\n",
    "        # get model predictions\n",
    "        y_hat = model(xs)\n",
    "        # calculate loss\n",
    "        loss = loss_func(y_hat, ys)\n",
    "        # calculate gradients\n",
    "        loss.backward()\n",
    "        # update the weights\n",
    "        optimizer.step()\n",
    "        # write current loss in the progress bar\n",
    "        pbar.set_postfix_str(f\"Loss: {loss.item():.4f}\")\n",
    "\n",
    "        \n",
    "def eval_model(model, loader) -> Tuple[torch.Tensor, torch.Tensor]:\n",
    "    \"\"\"Evaluate the model.\n",
    "\n",
    "    :param model: A torch.nn.Module implementing the LSTM model\n",
    "    :param loader: A PyTorch DataLoader, providing the data.\n",
    "    \n",
    "    :return: Two torch Tensors, containing the observations and \n",
    "        model predictions\n",
    "    \"\"\"\n",
    "    # set model to eval mode (important for dropout)\n",
    "    model.eval()\n",
    "    obs = []\n",
    "    preds = []\n",
    "    # in inference mode, we don't need to store intermediate steps for\n",
    "    # backprob\n",
    "    with torch.no_grad():\n",
    "        # request mini-batch of data from the loader\n",
    "        for xs, ys in loader:\n",
    "            # push data to GPU (if available)\n",
    "            xs = xs.to(DEVICE)\n",
    "            # get model predictions\n",
    "            y_hat = model(xs)\n",
    "            obs.append(ys)\n",
    "            preds.append(y_hat)\n",
    "            \n",
    "    return torch.cat(obs), torch.cat(preds)\n",
    "        \n",
    "def calc_nse(obs: np.array, sim: np.array) -> float:\n",
    "    \"\"\"Calculate Nash-Sutcliff-Efficiency.\n",
    "\n",
    "    :param obs: Array containing the observations\n",
    "    :param sim: Array containing the simulations\n",
    "    :return: NSE value.\n",
    "    \"\"\"\n",
    "    # only consider time steps, where observations are available\n",
    "    sim = np.delete(sim, np.argwhere(obs < 0), axis=0)\n",
    "    obs = np.delete(obs, np.argwhere(obs < 0), axis=0)\n",
    "\n",
    "    # check for NaNs in observations\n",
    "    sim = np.delete(sim, np.argwhere(np.isnan(obs)), axis=0)\n",
    "    obs = np.delete(obs, np.argwhere(np.isnan(obs)), axis=0)\n",
    "\n",
    "    denominator = np.sum((obs - np.mean(obs)) ** 2)\n",
    "    numerator = np.sum((sim - obs) ** 2)\n",
    "    nse_val = 1 - numerator / denominator\n",
    "\n",
    "    return nse_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare everything for training\n",
    "\n",
    "Now that we have everything needed, we have to do the following steps:\n",
    "\n",
    "- Specify a training, as well independent validation and test periods. For each of the three create a PyTorch data set and with these data sets PyTorch DataLoader. These DataLoader put single samples together to mini-batches that we use to train the network.\n",
    "- Initialize the model\n",
    "- Create a optimizer\n",
    "- Define the loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "basin = '01022500' # can be changed to any 8-digit basin id contained in the CAMELS data set\n",
    "hidden_size = 10 # Number of LSTM cells\n",
    "dropout_rate = 0.0 # Dropout rate of the final fully connected Layer [0.0, 1.0]\n",
    "learning_rate = 1e-3 # Learning rate used to update the weights\n",
    "sequence_length = 365 # Length of the meteorological record provided to the network\n",
    "\n",
    "##############\n",
    "# Data set up#\n",
    "##############\n",
    "\n",
    "# Training data\n",
    "start_date = pd.to_datetime(\"1980-10-01\", format=\"%Y-%m-%d\")\n",
    "end_date = pd.to_datetime(\"1995-09-30\", format=\"%Y-%m-%d\")\n",
    "ds_train = CamelsTXT(basin, seq_length=sequence_length, period=\"train\", dates=[start_date, end_date])\n",
    "tr_loader = DataLoader(ds_train, batch_size=256, shuffle=True)\n",
    "\n",
    "# Validation data. We use the feature means/stds of the training period for normalization\n",
    "means = ds_train.get_means()\n",
    "stds = ds_train.get_stds()\n",
    "start_date = pd.to_datetime(\"1995-10-01\", format=\"%Y-%m-%d\")\n",
    "end_date = pd.to_datetime(\"2000-09-30\", format=\"%Y-%m-%d\")\n",
    "ds_val = CamelsTXT(basin, seq_length=sequence_length, period=\"eval\", dates=[start_date, end_date],\n",
    "                     means=means, stds=stds)\n",
    "val_loader = DataLoader(ds_val, batch_size=2048, shuffle=False)\n",
    "\n",
    "# Test data. We use the feature means/stds of the training period for normalization\n",
    "start_date = pd.to_datetime(\"2000-10-01\", format=\"%Y-%m-%d\")\n",
    "end_date = pd.to_datetime(\"2010-09-30\", format=\"%Y-%m-%d\")\n",
    "ds_test = CamelsTXT(basin, seq_length=sequence_length, period=\"eval\", dates=[start_date, end_date],\n",
    "                     means=means, stds=stds)\n",
    "test_loader = DataLoader(ds_test, batch_size=2048, shuffle=False)\n",
    "\n",
    "\n",
    "#########################\n",
    "# Model, Optimizer, Loss#\n",
    "#########################\n",
    "\n",
    "# Here we create our model, feel free \n",
    "model = Model(hidden_size=hidden_size, dropout_rate=dropout_rate).to(DEVICE)\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)\n",
    "loss_func = nn.MSELoss()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model\n",
    "Now we gonna train the model for some number of epochs. After each epoch we evaluate the model on the validation period and print out the NSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "58058f6fbb5e46fda735ebf718502446",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.14\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7bdc746ef24a42259d08c8daca72367d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.20\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "820c931c02344161b16ee34b267c6214",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.26\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "04db65eeba9b405c98c74e7fa8cdd5eb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.33\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e3fee3b4a02e48ba96b037deff75c0e4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.41\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c9cbf34adf2049399e9c1de000e2ebe6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.45\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "49168d02c91f4bb2aa507340c3c51b96",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.49\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0f8e8463db944788bf9ded22c970c505",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.53\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "faf367a529f94ead98274830b277e3e5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.55\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7fdd7d1266534e80ba38aaa191dba9f3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.58\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cb93d6a94c3f4fe591be423b0285d5d7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.58\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "823fb4a444a44115934eadc02c4ec3bd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.60\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1ef8ee735a0f4a1194f050896a22dba7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.61\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0ef4e529fa2c45ffa315fdb30635d51d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.62\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "85f07691e82041beadb1c4d5a70eedf8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.63\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4520be453e0e4cc793c75b1a6e9c4b35",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.64\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7a1273d4a7c34999a5f7c00ac9135a7e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.64\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "06fc117bc823436b9896352ffffcc692",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.65\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "40220ac701da46aab10b179cb6e2b5f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.66\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3fb0c99ea4e84589a1945217d7f06e2e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validation NSE: 0.66\n"
     ]
    }
   ],
   "source": [
    "n_epochs = 20 # Number of training epochs\n",
    "\n",
    "for i in range(n_epochs):\n",
    "    train_epoch(model, optimizer, tr_loader, loss_func, i+1)\n",
    "    obs, preds = eval_model(model, val_loader)\n",
    "    preds = ds_val.local_rescale(preds.numpy(), variable='output')\n",
    "    nse = calc_nse(obs.numpy(), preds)\n",
    "    tqdm.tqdm.write(f\"Validation NSE: {nse:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate independent test set\n",
    "Finally, we can can evaluate our model on the unseen test data, calculate the NSE and plot some observations vs predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAEmCAYAAACZJK1dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydebwcZZX3v6f73uRmYd9EQIOKymIIGEAJCAwOAiKCAm5sLyKK8qq4zIuMOojLIDoQEBBh2EZZFGRT0FF2wp5AIISwJ4SQEEL27d7b3XXeP2rp6urq5fat6u7bfb6fT3+6u6r6qdO1nufU75xHVBXDMAzDMAzDMCqTabUBhmEYhmEYhtHumNNsGIZhGIZhGDUwp9kwDMMwDMMwamBOs2EYhmEYhmHUwJxmwzAMwzAMw6iBOc2GYRiGYRiGUQNzmg3DaCoi8jcROaHVdhiGYRjGUDCn2TCMWERknoisF5E1IrJcRO4Qke2G266qHqKq1zRo06YicouIrBWR10Tki6F5o0TkJs9uFZH9I78VEfmliCz1XueKiHjzthSR60VkoYisFJGHRGSv0G8/KSLTRGSFiLwpIpeLyAah+VeLyKC3rfxXNjR/kojMEJF13vukiG2ne+2uFJErRWR0g9tndmj9BRHpD30/s5E2vXZvEJEfNvr7Ku0eLCIv17FuR0R2DU3bRUT6Q993FZG7veN0uYg8ISIfD63DieybNSKyW502biEif/GOubkiclSVZa+OrGNARJbU25aIHCIiL3rz7xKRbeux0TCM5mBOs2EY1fiUqo4HtgYWA79psT0XA4PAVsCXgN+KyM6h+dOAY4E3Y357CnAEsCswETgM+Ko3bzzwBPBhYFPgGuAOERnvzd8I+BnwTmBHYFvgV5H2z1XV8aFXAVxnHrgN+AOwidf2bd50ROQTwBnAgcAE4D3AT4a0VTxUdWd//cCDwGkhe37RSJttwnLgp3EzRCQD3IG7jbcE3gF8F1gbWuzVyL4Zr6pP1bnuy7z1bwl8GbhSRHaIW1BVTwyvA7gFuLGetkRka+CPwPeBzYHncI8ZwzDaBVW1l73sZa+yFzAP+Hjo+6HAi6HvnwSeAlYBrwNnheb14d7wlwIrcB3Srbx59wEne59PxHV0f43rTMwFDqlgzzhch/n9oWm/B86JWXYBsH9k2sPAKaHvXwYerfL/VwEfrjDvM8Cs0PergZ9VWPYg4A1AQtPmAwd7n68DfhGadyDwZgL7L9jOkelfBV4AluE6m9t407PARcASYCXwNPAB4JtADhgA1gA3xrQZ+1tv3hhgqneMvInb8RoNbAasBxyv3TXAZjFt3wCcC7wN7OVN2wXo9z5vCygwpsJ2OBh4ucFtuAmQB94dmnZj+Fiv8tuNvP+3Vz1tedv5ntC8jb3jfcJwjwV72cteybws0mwYRk1EZCzwOeDR0OS1wPG4N/dPAqeKyBHevBNwnYbtcJ2jr+E6EHHshevEbY7rHF3hyyYivB8oqOqLoWlPAzvHLBvHzt7yNX/rySdGAZWkAx8DZkemfV1Elnnyi89G1vuMqmpo2jOhdcfZtZWIbFbtzzSCiHwe+DbwKdxo/VMUo5mH4Uba34vr4H0RWK6qFwJ/Bn6qbgT16JimY3/rzTsf17H9EK4T/n7gDFVdChxJaRR4aQXTV+EeGz+Pmfcm8BpwnYh8WkS2rGtjeIjIf4jITRVmfxBYraqvhabVe8x9Dpinqo/V2VbJcaCqK3A7V/Ue34ZhpIw5zYZhVONWEVmB67T8KyFJgqrep6qzVNVR1WeA64H9vNk5XGf5fapaUNUZqrqqwjpeU9XL1ZUzXIMrBdkqZrnxuFHMMCuBDWKWjSP6+5XA+KiDLiIb4kawf6Kq0fUhIv+K2yn4cWjyhcAOuI/dfwRcLSJT6rQ7zi6o/38Nha/iRsRfVNUcrgxkHxHZCnefbYjr3KmqzlbVt+psN/a3ItIDnAR8S1VXeNvzHODzDdh+EbBzVKuuqnnc424xbkR7kadv3j602PaeHj38ynq//4mqVtIpD+eYOwH3eK63reEe34ZhpIw5zYZhVOMIVd0Y93H6acD9IvIOABHZS0TuFZElIrISN5q8ufe73wP/C9zgJdedKyK9FdYR6I9VdZ33cXzMcmtwHbMwGwKr6/wv0d9vCKwJR4BFZAzwF1zZxn9GGxCRj+DKKY4KR7xV9UlVXaqqeVW9E7gWV8JRj91xdhH3v0TkzFCS2aU1/3E57wYu9R1HXDlFHjcS/DfgCuB3wGIRuSSk6a5Fpd++E+gFZofWeStu52JIeMfGL4iJNqvqa6r6NVXdHlcTDnBlaJG5qrpx5FWoY7UNHXMi8j7cJyhhTfJQj4O61mUYRvMwp9kwjJp40eKbgQKwjzf5OuB2YDtV3Qi4FBBv+ZwXwdsJ2Bv38f3xwzTjRaAnkoS1K+UyiUrM9paP/a1XseJWXP3xV4ngVVu4HThJVe+usS7F2xbeOiZGItoTQ+uOs2txnFRBVX8RkjJ8rYYNcbwOnBhxHsd4TwJUVc9T1d08+3YFvhX6PxWp8ttFuE75e0Pr20hVfelJ1XZj+B2ug39oFVteA36Lq3seLs8DG4rIu0LT6jnmTsDVJy8YQlslx4GIbITbyan3+DYMI2XMaTYMoybi8mlcveocb/IGwDJV7ReRPXF1rP7yB4jIh7xH4KtwH9/XE9mriKquBW4GzhaRcZ784dO4UW1/vaNFpM/7OkpE+kLO6v8A3xGRbUTknbgVFq72ftcL3ISruz5eVZ3I/98F+Dvwf1X1L1HbROQoERkvIhkROQi3gsft3uz7vP/+Tc++07zp94Ts+rKI7CQimwA/9O1KgUuBH4rIBzy7N/H11yLyERGZ7Ekq1uImofn7bDHFCG4ZlX7rSUCuBC4Qkc2942g7T+Lit7tlvRFtVR0Ezgb+X2jdW4nIj0XkPV77W+ImmD5aoZm6UdXlwF9xj7mxnjTkYNwnCbF4x9txRPZhHW3dBOwhIp/yjuGfAA+r6rzh/g/DMBKimVmH9rKXvUbOC7d6xnrcx8argWeBL4XmH4WbgLUa1xm4CPiDN+8LuMl9a3EdowuBHm/efUSqZ0TWq7ha6DibNsWNBq/FTZL6YozNGnlN8OYJbjLZMu91Ll5FC1xNrALrKFZyWAPs682/itIqD2uA2aH1PoirP12Fm8z1+YhduwEzvO35JLBbZP53vO20ylvX6AT2X7CdI9O/jBu9XOXtv0u96Qd7+3gNrmzjamCsN28nYBZuJZQbYtqs9tsx3rae561zNnBqaJ+Eq6xsGtP2DcAPQ997cJ86+NUzNsLtOM331r/I+75VyLbovlsDHObN/wlwS5XtuAXu8b3O+w9Hh+a932try9C0A7xjoayaR7W2vPmHAC958+8Ctm31dcBe9rJX8eXfMAzDMAzDMAzDqIDJMwzDMAzDMAyjBuY0G4ZhGIZhGEYNzGk2DMMwDMMwjBqY02wYhmEYhmEYNTCn2TAMwzAMwzBq0NNqA+ph88031wkTJrTaDMMwDMMwDKPDmTFjxtuqukV0+ohwmidMmMD06dNbbYZhGIZhGIbR4YjIa3HTTZ5hGIZhGIZhGDUwp9kwDMMwDMMwamBOs2EYhmEYhmHUYERomg3DMAzDMDqVXC7HggUL6O/vb7UpXUVfXx/bbrstvb29dS1vTrNhGIZhGEYLWbBgARtssAETJkxARFptTlegqixdupQFCxaw/fbb1/Ubk2cYhmEYhmG0kP7+fjbbbDNzmJuIiLDZZpsNKbpvTrNhGIZhGEaLMYe5+Qx1m5vTbBiGYRgJ8fHz7ucb1z7ZajMMY9jMmzePXXbZpdVmMHPmTO68887g++23384555zTElvMaTYMwzCMhHj5rTXcMWtRq80wjLYkn88P+TdRp/nwww/njDPOSNKsujGn2TAMwzAMo8s577zz2GWXXdhll12YOnUq4Dq5J5xwAhMnTuSoo45i3bp1AJxxxhnstNNOTJw4ke9973sALFmyhM9+9rPsscce7LHHHjz00EMAnHXWWZxyyikcdNBBHH/88ey1117Mnj07WO/+++/PjBkzePzxx9l7773Zbbfd2HvvvXnhhRcYHBzkxz/+MX/84x+ZNGkSf/zjH7n66qs57bTTAHjttdc48MADmThxIgceeCDz588H4MQTT+Sb3/wme++9N+95z3u46aabEtlGVj3DMAzDMAyjTfjJX2bz3MJViba50zs35D8+tXPF+TNmzOCqq67iscceQ1XZa6+92G+//XjhhRe44oormDJlCieddBKXXHIJJ510ErfccgvPP/88IsKKFSsA+Na3vsXpp5/OPvvsw/z58/nEJz7BnDlzgvanTZvGmDFjOP/88/nTn/7ET37yExYtWsTChQv58Ic/zKpVq3jggQfo6enhrrvu4swzz+TPf/4zZ599NtOnT+eiiy4C4Oqrrw7sPu200zj++OM54YQTuPLKK/nmN7/JrbfeCsCiRYuYNm0azz//PIcffjhHHXXUsLejRZoNwzAMwzC6mGnTpnHkkUcybtw4xo8fz2c+8xkefPBBtttuO6ZMmQLAsccey7Rp09hwww3p6+vj5JNP5uabb2bs2LEA3HXXXZx22mlMmjSJww8/nFWrVrF69WrAlVSMGTMGgGOOOYYbb7wRgD/96U8cffTRAKxcuZKjjz6aXXbZhdNPP70kGl2JRx55hC9+8YsAHHfccUybNi2Yd8QRR5DJZNhpp51YvHhxItvJIs2GYRiGYRhtQrWIcFqoauz0aHUJEaGnp4fHH3+cu+++mxtuuIGLLrqIe+65B8dxeOSRRwLnOMy4ceOCz9tssw2bbbYZzzzzDH/84x/53e9+B8CPfvQjDjjgAG655RbmzZvH/vvvP+T/EbZ39OjRNf/fULFIs2EYhmEYRhfzsY99jFtvvZV169axdu1abrnlFvbdd1/mz5/PI488AsD111/PPvvsw5o1a1i5ciWHHnooU6dOZebMmQAcdNBBgYQCCKbH8fnPf55zzz2XlStX8qEPfQhwI83bbLMNUCrB2GCDDYKIdZS9996bG264AYBrr72WffbZp/GNUAfmNBuGYRiGYXQxu+++OyeeeCJ77rkne+21FyeffDKbbLIJO+64I9dccw0TJ05k2bJlnHrqqaxevZrDDjuMiRMnst9++3H++ecDcOGFFzJ9+nQmTpzITjvtxKWXXlpxfUcddRQ33HADxxxzTDDt3/7t3/jBD37AlClTKBQKwfQDDjiA5557LkgEDHPhhRdy1VVXMXHiRH7/+99zwQUXJLxlSpGkQtZlDYtsB/wP8A7AAS5T1QtE5CzgK8ASb9EzVfXO+FZcJk+erNOnT0/FTsMwDMNIigln3AHAvHM+2WJLjJHEnDlz2HHHHVttRlcSt+1FZIaqTo4um6amOQ98V1WfFJENgBki8k9v3vmq+usU120YhmEYhmEYiZGa06yqi4BF3ufVIjIH2Cat9RmGYRiGYRhGWjRF0ywiE4DdgMe8SaeJyDMicqWIbNIMGwzDMAzDMAyjUVJ3mkVkPPBn4Nuqugr4LfBeYBJuJPq/KvzuFBGZLiLTlyxZEreIYRiGYRiGYTSFVJ1mEenFdZivVdWbAVR1saoWVNUBLgf2jPutql6mqpNVdfIWW2yRppmGYRiGYRiGUZXUnGZxK0xfAcxR1fNC07cOLXYk8GxaNhiGYRiGYRhGEqQZaZ4CHAf8i4jM9F6HAueKyCwReQY4ADg9RRsMwzAMwzCMJjN+/HgAFi5cyFFHHVV12alTp7Ju3brg+6GHHsqKFStSta8R0qyeMQ2QmFlVazIbhmEYhmEY7UehUCCbzQ7pN+985zu56aabqi4zdepUjj32WMaOHQvAnXe2p6toIwIahmEYhmF0OfPmzeODH/wgJ5xwAhMnTuSoo45i3bp1TJgwgbPPPpt99tmHG2+8kVdeeYWDDz6YD3/4w+y77748//zzAMydO5ePfvSj7LHHHvzoRz8qaXeXXXYBXKf7e9/7Hh/60IeYOHEiv/nNb7jwwgtZuHAhBxxwAAcccAAAEyZM4O233wbgvPPOY5dddmGXXXZh6tSpQZs77rgjX/nKV9h555056KCDWL9+ferbKM3BTQzDMAzDMIyh8Lcz4M1Zybb5jg/BIefUXOyFF17giiuuYMqUKZx00klccsklAPT19TFt2jQADjzwQC699FJ22GEHHnvsMb7+9a9zzz338K1vfYtTTz2V448/nosvvji2/csuu4y5c+fy1FNP0dPTw7Jly9h0000577zzuPfee9l8881Llp8xYwZXXXUVjz32GKrKXnvtxX777ccmm2zCSy+9xPXXX8/ll1/OMcccw5///GeOPfbYYW6o6lik2TAMwzAMw2C77bZjypQpABx77LGBo/y5z30OgDVr1vDwww9z9NFHM2nSJL761a+yaNEiAB566CG+8IUvAHDcccfFtn/XXXfxta99jZ4eN2a76aabVrVn2rRpHHnkkYwbN47x48fzmc98hgcffBCA7bffnkmTJgHw4Q9/mHnz5g3jn9eHRZoNwzAMwzDahToiwmnhFj4r/z5u3DgAHMdh4403ZubMmXX9Poqq1lwmunwlRo8eHXzOZrNNkWdYpNkwDMMwDMNg/vz5PPLIIwBcf/317LPPPiXzN9xwQ7bffntuvPFGwHVqn376aQCmTJnCDTfcAMC1114b2/5BBx3EpZdeSj6fB2DZsmUAbLDBBqxevbps+Y997GPceuutrFu3jrVr13LLLbew7777JvBPG8OcZsMwDMMwDIMdd9yRa665hokTJ7Js2TJOPfXUsmWuvfZarrjiCnbddVd23nlnbrvtNgAuuOACLr74YvbYYw9WrlwZ2/7JJ5/Mu971LiZOnMiuu+7KddddB8App5zCIYccEiQC+uy+++6ceOKJ7Lnnnuy1116cfPLJ7Lbbbgn/6/qRaqHvdmHy5Mk6ffr0VpthGIZhGFWZcMYdAMw755MttsQYScyZM4cdd9yxpTbMmzePww47jGef7a4x5+K2vYjMUNXJ0WUt0mwYhmEYhmEYNTCn2TAMwzAMo8uZMGFC10WZh4o5zYZhGIZhGIZRA3OaDcMwDMMwWsxIyDHrNIa6zc1pNgzDMAzDaCF9fX0sXbrUHOcmoqosXbqUvr6+un9jg5sYhmEYhmG0kG233ZYFCxawZMmSVpvSVfT19bHtttvWvbw5zYZhGIZhGC2kt7eX7bffvtVmGDUweYZhGIZhGIZh1MCcZsMwjA7FcZRn34gfmcswDMMYGuY0G4ZhdCi/vf8VDvvNNGa8trzVphiGYYx4zGk2DMPoUJ5btAqAN1asb7ElhmEYIx9zmg3DMFrM359dxJsr+xNvNysCWP1XwzCMJDCn2TAMo4U4jvK1PzzJ0b97OPG2M67PTMExp9kwDGO4mNNsGIbRQnx39vVlyUsoMp7XbD6zYRjG8DGn2TAMo0PJiO80m9dsGIYxXMxpNgzD6FB8eYZjoWbDMIxhY06zYRhGC0kzSS9r8gzDMIzEMKfZMAyjQxFPnlEweYZhGMawMafZMAyjhaTpzlrJueZi29kwOhtzmg3DMDoUKzlnGIaRHOY0G4ZhtJA0g5Mipmk2DMNICnOaDcMwOpSMyTOaim1mw+hsUnOaRWQ7EblXROaIyGwR+ZY3fVMR+aeIvOS9b5KWDYZhGO2OpqhqznpXeKvTbBiGMXzSjDTnge+q6o7AR4BviMhOwBnA3aq6A3C3990wDMNIGD/SXHBabEiXYF0Tw+hsUnOaVXWRqj7pfV4NzAG2AT4NXOMtdg1wRFo2GIZhtDupBoG9RMA0o9mGYRjdQlM0zSIyAdgNeAzYSlUXgetYA1tW+M0pIjJdRKYvWbKkGWYahmF0FIKvaW6xIYZhGB1A6k6ziIwH/gx8W1VX1fs7Vb1MVSer6uQtttgiPQMNwzA6FE+dYTQJS7g0jM4mVadZRHpxHeZrVfVmb/JiEdnam7818FaaNhiGYbQzqZacC9ZhzpxhDJe3VvWzw7/fyYzXlrXaFKNFpFk9Q4ArgDmqel5o1u3ACd7nE4Db0rLBMAyjmymWnGuxIYbRAfx99pvkCspNMxa02hSjRfSk2PYU4DhglojM9KadCZwD/ElEvgzMB45O0QbDMIy2Js0kPQkSAY1mYNu5s+nPFQAYPzpN18loZ1Lb86o6jeLTwSgHprVewzAMw6Uoz2ipGYbREfiJtTbCZvdiIwIahmG0kHRLznnyDIuBNgXrnHQ2wZMb289dS81Is4hMBvYF3gmsB54F7lJVU8IbhmGMAOwmbxjJceVDc/nxp3ZqtRlGC6gYaRaRE0XkSeAHwBjgBdxKF/sA/xSRa0TkXc0x0zAMozNpwtgmFmc2jASQLqzhWDAtSgnVIs3jgCmquj5upohMAnbATeYzDMMw2gwxUXNTMRlMZ9NtLvOC5evY55f38qujJnL05O1abU5bUDHSrKoXV3KYvfkzVfXudMwyDMPoDtKsoRyMCJjaGgyje+i2QPMrS9YCcPvTC1tsSftQMdIsIhdW+6GqfjN5cwzDMIyksMSl5mLbubPpMp+ZjPeHHTuwA6pVz5jhvfqA3YGXvNckoJC+aYZhGCOfl99aw5sr+yvOb46mufJalq4Z4Kd/fY58wUnRktbzwpurueiel1pthjGCyWS6y23Oer1up7MvDUOiYqRZVa8BNyEQOEBVc973S4F/NMU6wzCMEUzBUT5+3v0AzDvnk01ffz2R5p/85Tluf3ohe0zYhIN32bo5hrWAIy95iHWDBb6633vpzVq1VWPodJfLXEx8tEhzkXquHO8ENgh9H+9NMwzDMKrwyQsfrLlMmvcjkdqa5oF8dzw4HMy74TK7/xsN02Wi5ozJu8qoZ0TAc4CnRORe7/t+wFmpWWQYhtEhPP/m6labAFS/6fkVpTq9nFZGBNCmRc1UteO3abfRbXvTl6MUzGsOqJYI2KOqeVW9SkT+BuzlzTpDVd9sjnmGYRgdTqqRZn8VlVfiV+/IdLiD1+ykSEch29mbtOvo8FOkjIzJM8qoFml+VEQWAH8H/q6qtzXJJsMwDCMB/JJz1RxzP9Lc6TlOvgOQZtQs3LTbGenwjdplSJftz2L1jNba0U5USwScLCLvBg4BporINsA04G/A/ao60CQbDcMwOpY0B8SQ2j5zMOJXt0SamybPaMpajGbS4adIGf41Ic1a8iONqomAqvqaql6qqkcAewN/AT4OPCgidzTDQMMwDKMxigMCVr7p+U5kpzsEwaPmFMNm4Q6QPdLuPDr8FCkjeDpjoeaAehIBAfBKzt3jvfAiz4ZhGMYwSLd6Ru11aCDP6GyXQJr8qNl85s6jw0+RMjJeWNV85iI1S86JyGEi8pSILBeRVSKyWkRWqeobzTDQMAzDaIx6htF2aiQCzlm0ijUD+aRNazoWNTOGS7dpmoPrh/UAA+qp0zwVOAHYVFU3VNUNVHXDlO0yDMPoClIdEbCOSHPRaS6fp6occsGDnHT1EylY11yK26I5iYAmz+hAustnDrBjuUg9TvPrwLNqXQ3DMIwRSbVkw2p1mv15j89dloZZTaUZ1TPC2B2z8+hSn9mezoSoR9P8b8CdInI/EFTMUNXzUrPKMAyjS0gzHhGMCFhV01w50txJEaZml8/qnC1n+HTbYDV+Z9uO5SL1RJp/DqwD+nCH0/ZfhmEYRhtTzy1eq0SaOyvClHz1jLlvr2XdYFHvHW7ZHs52Ht3lMtvTkjjqiTRvqqoHpW6JYRhGF9IMTXM9y8RFlYsOdYJGtYhMlf/ZKAf8+j6mvG8zrj35I2XzOqq/YQCdcR4Yw6OeSPNdImJOs2EYxgilWtSzmCFfPq+z5Bn+kMDJtvvQy0uDzyXbuXM2ndHt2LEcUI/T/A3g7yKyPlxyLm3DDMMwuoFm+KVVV1GlqkRnOc3ue7MkJ2mO9Gi0hg46HerC/79d9rerUlOeoaqmXzYMwxiBBDe9Kne9YNTAmHmOk7RFrUOaPCSwyTOMkY51/Mqpa0RAEZkITAgvr6o3p2STYRhG19CMG1O1dVTTNAdDbKdiVXPx/2dSJefinG9LBOxsum2PFjvd3fbPK1PTaRaRK4GJwGzAjzsoYE6zYRhGG+Pf6qpFPbtO05xQ9LzWpumcLWd0O3YsF6kn0vwRVd0pdUsMwzC6gGN+9wgz56/gxZ8f4k5ohqa5mtNcNdKckkEtII3qGdXopA6H4RKOuKpqx9dttiO4nHoSAR8REXOaDcMwEuDxucsYLDRHLFy8yde+/cVqmjvI8ZOgekZC8oy4aVpjAaNj6KBToyL+9aMb/mu91OM0X4PrOL8gIs+IyCwReabWj0TkShF5S0SeDU07S0TeEJGZ3uvQ4RhvGIYx0mn1/UiqZAJ2ltPsvjeveobRaYT3aSedG5Xo/H84dOqRZ1wJHAfMoqhproergYuA/4lMP19Vfz2EdgzDMIxhUL16RuUIbGfJM5Kt0xybHBWa1A1O1UjhxcWr2XqjPjbo602szW7au1ZFo0g9TvN8Vb19qA2r6gMiMmHIFhmGYXQRafpWdZWck8rLJDnkdKtptqbZfOb24aDzH2DX7Tbmtm9MGV5DXdYpquf60W3U4zQ/LyLXAX8BBvyJwyg5d5qIHA9MB76rqssbbMcwDMOog3oiRdVKzqXJYN7BUaWvN5vqeoKIekIdgVqtmJ/RXjz9+opE2+sOR7Ir/uSQqEfTPAbXWT4I+JT3OqzB9f0WeC8wCVgE/FelBUXkFBGZLiLTlyxZ0uDqDMMw2ps0H336bVePNHsl52Lm+f5lmlUCDrngAT74o7+n1r5P0nWa4wjvy06K0hsu4f3bDU5zN/zHoVLPiID/J6mVqepi/7OIXA78tcqylwGXAUyePNl2nWEYRoNUu4AGeYAtijS/smRt6uuAoqY5qb9kDkVrUFVeWLyaD75jw9ba0UVRWDvWi1SMNIvID0Vk0yrz/0VEhhRxFpGtQ1+PBJ6ttKxhGEY30M6a5k4aCSzj3e2Sqp5Ry2nqoE3XVvz3g3M5eOqDPDW/PmVnksdwuKlueJDQBX9xyFSLNM8C/iIi/cCTwBKgD9gBV15xF/CLSj8WkeuB/YHNRWQB8Jipz2wAACAASURBVB/A/iIyCXdfzAO+Ovy/YBiGYVSj6jDa3nucE+AETvfIv31mEq7THEepUzXyt1k7Mu3ltwFYvm6wruWT3A3dNkx6F/zFIVPRaVbV24DbRGQHYAqwNbAK+ANwiqqur9awqn4hZvIVw7DVMDqG9YMF/j57EUdM2qbjR5UyqpPmfUnLPpRT1DRXXqgT7p3FzkFCkWYbRrslrOrPAbDRmPrKx6W1H7oi0mxecxn1aJpfAl5qgi2G0TX8/M7n+MOj89lqgz72ft/mrTbH6FACeUaVZapFmjup5FQwImCKgzF2WySyFfgJltlMPXUM0pNndFOvyI7lIvUddYZhJMrbq91HiyvW51psidFqmnFDqraOoqa5fJlOSnaSJtdp7oZI5EggvUhz5+9gjbwb5jQbRkvo7XFPvVwhxbCX0fUEJeeqLpVsVYl2pRma5lI6fIOOEJLVNIdKziXXbNvS6deERjCn2TBaQK83PFmuYFelbqFStLcZN6a6qmfEuAGddNMsjgiYTHu1qo100rZrR2o9oXlt6Vqefn1Fak9LBvIFvvOnmSxYvi6V9tuBeuq8dxs1nWYReb+I3C0iz3rfJ4rID9M3zTA6l56sewfPW6S5a4jeeJojy/DeqywTaJo7/FD0Nc1JlZyrhckzUqLKYDxh9vvVfXz64oeSjTSH2nrwxbe5+ck3+NGtnV85t5NkWsOlnkjz5cAPgByAqj4DfD5Nowyj0+nJevIMu7N2DdE93SznDerUNMf+Lh17WkEmYU1zzTrN5mikQnEwnvqWT+sYTnNkybahC/7iUKnHaR6rqo9HpuXTMMYwugVfnmGR5u4h6qw186ZbPdLcbK1va2hKnebw587enG1AfRs4yc5LuCX/OOrkkqFBIqAdywH1OM1vi8h78bafiBwFLErVKsPocCThIX2N9ie6r305RFOOgTo0zXHLNDNa6qQceQ+qZyTUT6213zq9E9JqWh1p7obd2w3/cajUrNMMfAO4DPigiLwBzAWOTdUqw+gS7JrUPbQi0uzLMqqOCNjkUmyVcFTJkF7Uzo80J7Xda8lZzOFIl3o3b6K7oSTRs3t2cPf809rUM7jJq8DHRWQckFHV1embZRidTbXauEZ34Gua04zm1hqcZMnqAfKFymXpmnl4pi3xLj7dsXOuE6g/0pzO/u6GdBTT5ZdT02kWke9EvgOsBGao6syU7DKMjkZSjKgZ7Um5PKN55ZwqlUfb4+d3Fe1pg0hzmviJgEmlEdRyxlq9PZNGVXlx8Ro+8I4NWmrHUJ+MJLkXYjXNtX6jyktvreH9W7V2uzVCJ40ImhT1aJonA18DtvFepwD7A5eLyL+lZ5phdC7FSHNr7TCaR0vkGcF77RrMrT4W03Yyi8OFp5kI2Ll1mq+YNpdPTH2Ap+Yvb7UpQOs1zfVGmq95eB4Hnf8Aj89dlo4hKdJhh3Ai1OM0bwbsrqrfVdXv4jrRWwAfA05M0TbD6FiKuVftcVkqOJp6Ila3E926QaS5GetWWL52sOTGHV1vq2ULaZfgS7p6Rq1WOu1sevTVpQC8tXqgxZa41HvtTPK4LtWs19fuswtXATBv6drE7Gg+nXY0N049TvO7gMHQ9xzwblVdD7TH2WMYI4x2izS/98w7+eylD7fajI6mNYmA3jtw7BWPcczvHqnYOeoWTXOqncNQ050gz3j5rTU8Mc/taK0dKAAwblQ99QOaQJtUz6hVca5YV3rkHQ9BIvHIMz016jn6rwMeFZHbvO+fAq73EgOfS80yw+hg2rG251PzV7TahI4meuMJEgGb5DzPWeRGvPKOMiojZett9ZOGtLdDoGlOaDW1zO0ER+Pj590PwLxzPsmgJwYf3VtPrC19WlE9I3yM1tspyozg8qIj0OTUqXn0q+pPga8AK3ATAL+mqmer6lpV/VLaBhpGO/LGivUsXztYe8Ea2EWpi6hQpzndVRZVzcEolJ7zUybPqPr79ElbntH8ijWddXZn2uTp2FC16a2unpHJDG35elBVZi9cmVyDFVeU/ipGGlWdZhHJiMizqjpDVS9Q1amqOr1ZxhlGuzLlnHvY55f3NPz7oQ4Fa4x8KskzmqVp7vG8nsBpbrtEwHTbD+o0J7WiuIokoc+dliIgTRhRcSjUa0aSMqi46hm162ckv92ufngen7xwGo+8sjSxNqvRHnu8PajqNKuqAzwtIu9qkj2GMWJYO1ho/Md+1MYuR12DAgP54jGTdmQVSjXNWc9pHgwizRF5RsxNvZn+UfryDN95Saa9Wudum/iWiZAvOGSboQmvg6Dedp3L55PS4zRIEKFPsE1fajV/WbrJhf4xPhL12GlRj6Z5a2C2iDwOBHtIVQ9PzSrD6HD8Os12LeoeHFX++dziku+Q7jEQiDNUQ5Hm9jzoUk+MbPLIh53kaBRUg05XMxJY66He7Ztk5zS8yno7D5kUBtURku0AVqJNdnVbUY/T/JPUrTCMLqMN8wCNlFEt3kChOZFmH0ch64kr8xXkGbGR5tQtC68/3fYzCUdK4weMCX1OZC3tgeMUr1nNPG6rUa8VuaRGs4mss97NEAzGkuB283XSaTu14SdVhks9w2jf3wxDDKObGMlliIzGULSC05xmqLmom/YjzYP5eCei1Ydi2o/9M0GkOdXVBLSL9jcJCqqJ17lulKFeO/Mp7fB6I+6ZIcpJ6qHZ+vIOOpSHTc3qGSLyERF5QkTWiMigiBREZFUzjDOMTsUizd2HG2kufm+m86GqZU5jeaQ5/nfNolkjAiYlL6hZbaSDHI2CE5JnNKHqSz3UuxuT1DSHz4diDeP62k/Sd29W0KWDDuHEqKfg4kXAF4CXgDHAyd40w+hKkng8aZrm7qOSPKM5mubitEBLHbkltjopNXV9pv/erOhcU9bSHBxHQ9VH2sNrrttpTsle//yt1dlLQ9OcdFJrJexJaDl1VSlX1ZeBrKoWVPUqYP9UrTKMNiaJiFhQM3bYLRkjBUc10CL635u57lqPdGM1upHvuYLD22vSGQg2ba2s//+SWk8th6LVMoYkKYSeVLRNpLnO5dKSZ/jN1mo+jfrWxadGzYk0m/NcpB6neZ2IjAJmisi5InI6MC5luwyjbUny+mHXou5BiUaai9NTW2eMFCNumvu9tiXf+dPTTP7ZXancRJt1Y051FO2Y7dwJlESa2+SPtaJ6Rtz6azmukoKDK02KNPu0xx5vD+pxmo/zljsNt+TcdsBn0zTKMNqZJB5jB5o0uxx1DaqVEgHTJ3zDLsozosuU/y56n//rMwsrLjtcmiXPSMp5qdVKJ53ZhdBTktbXafbsqNOMRKtnxMicapFGImAako842qR/1FbUUz3jNe9jP1Z+zjCSuZCIaZq7jaimuTl1mssd5GIiYETTPAQ7wolhSZFGJyJXcJj20tsc8MEti5HBJjl9HSXPCEWa05I7DJ06q2ekVJe8KM+oFWlOvtKFpCD5iEdL3oz6qmdMEZF/isiLIvKq/2qGcUZncsPj8/nN3S+12oy2wK5F3UO0ekaQCJjiUaAxDnLlSHN8PYgwvvlpOIRptHnh3S/xf65+gmkvvR1MS6x6Ri0NeAed3KrFESVbHWn2qXsY7SQHN6H8PKqVZzhUB/eWpxbwq/99vuoyTdM0t8eubivqGdzkCuB0YAYwjHGDDcPljJtnAfB/D9yhxZa0kCZHvYzW4yYChuQZTYg0+6iGb97x6x2KGSPFaX592ToAFq/qj60kMhzCDpSGEi3j5o90Cu2oaa5zuSTlGWGcOju9mch5V4vT//g0AN//xAertGma5lZRj6Z5par+TVXfUtWl/qvWj0TkShF5S0SeDU3b1Itav+S9bzIs642O4YuXP8ov/169d90uJHHPqPfRXitwHGWvX9zFn2csaLUpHUU0ETC46aYqz/DWVRJprrBs3IiAkUlpJiCl0WZv1hsF0XGCjZGGDKRYPjC0ndukykQSFFTbZkTAoZbrTG0YbS19r0QaDq4EOumUI83+exvep1pFRadZRHYXkd2Be0XkVyLyUX+aN70WVwMHR6adAdytqjsAd3vfDYOHX1nKb+97pdVmNI3g0V4bXosG8g6LVw1w5i2zWm1KR6EhxwOaI88I1k3McL7RSPMQNc1JcPofZwaf0+hA9nhO82BBg+2c2HpCzcRFX9vw1G4Yx1GybTIiIEOUJuQSvMjOXlgc162Yk1BD0xxZPgmapWlu9a5uR6rJM/4r8n1y6LMC/1KtYVV9QEQmRCZ/mmKN52uA+4D/V8NGw2grknBy6i2M3wpyXogs6USvbmTsqCzrBl1Vm6MVBhlJM9IcEw0r+sylK447FqNTkh6J7Jan3iiuP4UeZG/WS14rOKFtkYK0JCaqPNKjc8vXDgafC1pM/KyWWFdwFEc1iPCnSb1bN5+gPOP2pxcGnyuNrBlFUkj6zkQ7wCnRSRKjpKjoNKvqASmsbytVXeS1v0hEtqy0oIicApwC8K53vSsFU4xWsG4w32oThk0SF79coX01zYN5z2m2sb6HzdhRWXoywqr+PKpa4kgFdZqbomnWss/1aJor2ZbGI/pU5RkFrTtxq17C5ob16e+RhSzSTdvyKdJQOPXaGcHnglN7cByAz1zyEE8vWMm8cz6Zun31dkrSqp5Rb53mNEcEbNYhNsIP5USpp3rGL0Rk49D3TUTkZ+maBap6mapOVtXJW2yxRdqrM5rET/86p9UmDJskLiD+ULTteGP1E2cyFmlOBF8iUFAt2d+Bo5Vm9QyKzpyvBa2gzqgrAjvUGrlDIR15hmvwYCjSnEYiW7gTcc/o73HNqF8y0l2NN1f2B58dh9CIgJX/19MLVqZt1pArR+RSEpfXK7FL45xJs4pNGFUYx/oR/9QkSep5hnKIqq7wv6jqcuDQBte3WES2BvDe32qwHWOEsjSlIXhHGvk2lmcEkWZzmoeNKvSEHmuXlKxqQiJgsK6YknNlxEyO3iyLSVhpSBxSkGdkipFmv/XEBjfR8Ge/J+KeO3tmXhjxetBxo4sPosPyjFZXz/CfHixc0V9jSZdcPhmnOXrMF58U1Yo0e8sl2Ilq1oiAY9a+zuy+L3O03J3uikYQ9TjNWREZ7X8RkTHA6CrLV+N24ATv8wnAbQ22Y4xQRvWkr3VLmyQchnbTNO+fmcnOMhcIRZpNnpEIxQoOWpo8FiQCpkiM7rJSAtNQjsU0HKc0HAD/aYmjWnH48CQIoq9OUX7Wjk+RhkJPuDxiuORcSnKHehnt3UN+9b8v1LV8UoOxRJvRuiPNyTu4w5V8PPDiEn7/6Gs1lxu3yk3O/7g80dB6OpF66jT/AbhbRK7CvQSfhJvEVxURuR436W9zEVkA/AdwDvAnEfkyMB84ukG7jRFKJaf589l7+Ez2QSB9LdxwSeLaN/215UB7OM2qytWjzvW+nRborc1nHj5KJBkN2JRVnNpzO05hRyDdYyAcXY1m3EfXWnOwDghVLkjGvjBpOOLZkHOxbf9LnNHzV+YUvpNI2+HIoW+7aCF2/kgkLM9ytH3qNA818DKYUCKgf56OoR+HTGiQoOrbQ4YoJ6mH4bZ5/JWPA3DcR95dfT1+xRnsZuBTzzDa54rIM8DHcS+ZP1XV/63jd1+oMOvAoZlodBKjYrKqHUc5p/e//S+QGfnR6Fq8/NYaoD2iUY5CtuR7fCTSqI87Zy1i93dtwjs26gPgnfI2q8iT9yoL/Kj39xyZfYhpb/4r8B5U4YjMNJ5wPpCqXT8YuIAPjHqOV/VBoOgk95JnY9YM6VhMQ0qRRuchE3Lyz3jzdPp61vPQmwpcN+y2S6L3/qP6QtFpbodzezhkIzXFM5pnU1a1PHm53socGXH3QS6frBxnTt9JLNJN+U+9GaidWBo8sUs00uy+D3dXxA3KE8bvBGa0g4qOD5N6EgHHAf9Q1e8BlwGjRaQ3dcuMjiQuSlCSqJFb10RrGiPJe3urb0BQ7hzHlSm7+ckFfOL8B5po1cjEcZSvX/skR136MOBu2+vWnMzjo7/u6mrVdVLdhXPee56poy7h5tH/kbg94Qz/g/P3sn1mcVnJufN7L+GJvq+jMTfGssFNvPc0HNw0OmmZkA63T9cDMGXVHYmvJ4i+hiLN7XBuD4fo6JWHvfafPNn3NbSQa7otB099gBunvw5AT42gynMLVzHp7H8Ex3k+oUTA8DG/tSwLJQLWp2lONtKcTBm7gVp6b2/bjRHLRfKpp8v2ANAnItsAdwH/B3fgEsMYMnHJZSXZ2CPAaU4yYtAO8oxK9/bwfvnOn57mhcWrR7wjkDb+/lywfH3J9B5xyDmuPEP9ZDo/Kum471vJCtKiNCpaKJl3WPZRoFRaUIt0RtVLvMnQaGwpOPnALaN+zOk9NwXnhYQ0zWkN39wsSiPNsPMyNxlMC4OVfpIaz7+5mu/f9Iy7fu8CvM3GY2KXvWLaXFasKzr2SckzoodQvXXWo1VrkiCpMnYlCcKOMn3essgS7jVhLPUlXXYD9TjNoqrrgM8Av1HVI4Gd0jXL6FSij/ygWLMYgOWvwdwHm23WkEhSq9gOPmjUoag2CMT3b3qG/lz9zlW3UaYTDn32awU73mXXj+zGRXgTsydmX0phoMS4nLrinEyMMxQ91kdayblMSMe9JrMBAMuymyfTuCq7ZV7mWz03hxyo4rmR5KAarSAbiTQ74om4mhxpruQYbrtJvNMcVW8kVac5enwGkpw6Nc2pDG4yzEbDSZIX3vMSR136CDO8fBsACXTcFmn2qctpFpGPAl8C/Oda9SQQGkYZ0QsxwA2Pzy8ucMXH4ZrDYN5DzTatqWzV53B45iGcNrixVrruhqPKfjLbn59cwM1PvhH/A6PqjbHgOCVetHoRXy2kPeCPMlbXFL/mB72pLr4Tn9HqdoSPhxEjz/CT1xzlobFuOs0zY/ZMpG0NRewLMZFmBtcmsp5WUZII6FR2mv/y9EIuvPul1Owo66BV6dRDMbLrk1TnpcxpHnL1jOSO70xC8oxwJZRZXo3tt8NlYb0O/WiaL8lpV+pxmr8N/AC4RVVni8h7gHvTNcvoVDKRMkYAv/rbs+ULzry2WSYNmSSufd/h91w46mLes/ap4Tc2TKKREv97OEs+nHyzfF0Cj2cLOUjdWRw6awfy/P7R1xp24CpF7cF9oqIom7AagGzBe+Q5BFnEUFHgiMxD/GXNl4Jp4kWUg8E+/NuAxtwYQ/bnHQ09ah4Z8oy4KgNZJyF5QegJQdB+aF9Kbk30FyOKaMk532nOOKXHyf+9/inO++eLqdkRlQJphek+UcnzYEKR5mgrTg3nvez3qVTPGF474Wv8xmNHAbAifH33wulBHoZRV/WM+4H7Q99fBb6ZplFG5xKWZ/gXvdgT8l0fbZZJQyaJS98WuhQE+gqtv7GGI4gaqmcbzp8J30ATufj/dHPYehJ89f7ayzaRn90xh+sfn892m4xh/w9smWjbecdBNRuUccrmPf2+k67c5aOZ50q++/KMoHPkOc1Zp/qNMVdwgpt1kprmDK7WO6mErZK2QxG5jOfQZmtE1OtFQ/YGDn9oWmZgdSLraRXhOu0FVQpejR0txD+qH0xoEJEo5R1Rv1Mfv3y0vnxSkeaoiqrScPRlvwtK0yVHUtHr8Hm86Ti3vsPykB7c/9PmNBep6DSLyFRV/baI/IWY/a2qh6dqmdGRhCPN+SpOcyE/WFIGrdPIB9G91uuDndDj1vD4G06FSHNiAZNFMxNqKDn8KMvagcb2S6UbPPiaZhjAvTll816yYHhADEcTHb5c1S15V4LvNPudo6D4cvl5GNVk+5YlGWi+adRZbCdLeKDwcHKNeoS1nxI4zck8ahbCTrO/MYvbUHIjW54R1gY7jqJepFkqPCFKo9MDMeeU916osL6o05xUQmZleUZ9J0OiHc2ELhFhm8aOct3BleuL54d/jI8yeUZAtUjz7733XzfDEKM7iI80lzsohcH1bes0JxFpzWsGZGgVC9JCC6WPmWVgJX8ddSbfyX89mF460EFTzWsqwx0qOPqzTMixyjuKqgaauJ4g0ly6zKgEneat185m32yp/EkimtQg0lzDmQxXIUjSAdg983LibfpkQyMC+vuix0nIAXCqyzMyI95pjk8EjEsYhdKE7lr1f4dCmTzDlxVV8IWjFZpyKSUCBiN61tl8kod3UqUfw9s2VvbiPQUbZZHmgIpOs6rO8N7vF5EtvM9LmmWY0ZmEoxf5KnqpwuD6smktYdVCeHMWvP8TwaQkrn0F9UNg7ZAIWNz+BUcZv/Bh3puZx/ezNwBfBUojG+1QJi8tAierwTtc9FfZ0LHtjwjoO2/Zguc0a3lCWVJsMrCwbFpRnuGt008EjIkghnd13nFSSWryyaVwLoSHMM4kHGnWOE2zE3aaR0D5zCpkQ+LgglOs+iJayWkObw/IJtT3Kxu+2p9e4VyJ+uphu4bjzEfX1uNth1rnQnF2gpHmTDLnYT4izYNotRHvHi2tD+60CxUTAcXlLBF5G3geeFFElojIj5tnntFpxCUC9kr5zdppF6f5ioPgumNKJg3XX1DVQJ7RDpFmJ+SsaD5Hrm8zADaTVcH08CPPTh4pMFxtoRGiN7FwhzCXL6CqZP0bkec0l1RhSHjbaszwt8VEQO/xcnAsVncmwyOrpeE0J1UaLExQp9kpRpqT0jSHRa5B9YxwpLnQJtewBgk7vaoEkebokwqf8P5LsvMXdY4D565OeUbYMRyOWdFjPuv0e/ZU/52fO5BknzDcGRwOJZFm72N4u0oHX+sbpVr1jG8DU4A9VHUzVd0E2AuYIiKnN8U6o+MIyzP8i2zcox8n1ybF1Fe6o1AlmawVjtq0g6a5xGnLrUczruY2XGYofBsa1oV6+Tx45sZhNJAugdOckDyjN+yg5fvdhLQg0uwd4yEdbCEFx7GMUPWMIzLTgkFVolURoLSySs5xQo+FkzcrjcFAwppmP9LcUyFSOmTCo//5yV6haH12hMszogGOgrgPpqWCvCVXiIm8J0AlTXOlYzCqbsolJCtSLdWx9xR8p7m+SHOS9f2LuQXJyTOcoDNSnCY2fHYZ1Zzm44EvqOpcf4JXOeNYb55hDJlwQlk1TXPbRJp9Qjq+4V78cgWloH50r/UXJSd0o3dy/cFVPhu6QYQfaQ7rhnjlwXDzyY3/PmV6hivPqBJp1vx6FMh6x0/Gd5Y1rGlO9niIO74yweAmytRRlxSnV+jAbcMSQN3qGTiMZjCxSKJENN9p4YQ6Kz2pRJrd95JIc754Dbv/xSWl9ehHAO9b9zT/3fsrBMfVNPtPJCqU7EvKOY1S1oENNM3x66iWCDica5ejyhiK/913mmv91VpOfiMUn6AMr51CTBTe77jnCg5/mxWqyW9RZ6C609yrqm9HJ3q65t70TDI6mXAUwL9J9sREmjXfJpFmn3y44Pvwmpp614uh2rhtEGkO3fw15DSHHb7tdCE7i9t/HtbfX71oOL9OnUzCiYDH8I/gs+TWu5Fmcbf3yrXrWLpmIBJpTjZLXWL2VuD0RBzquEjz2GXP81Dftzgp+3fyBeUMuZoX+k4skfQMh3CH2UlxpDnHCUWaE0pqKi055yeFFf9PT6GoaT7hysc54+ZZiay3WXzhjV/w8exTHJZ5tGRwk7jjBEo7PUnKjKKHWlAqsWKd5sqJgMMxy1EYFxpOOlvoR3Dq1jQnO7iJ1/Ywb0aliYClkeZ/zF5MLh+6P+VtVECo7jRXe4bV/MHnu5B1g3lW9XdWqZfwKe6fsHHyDM21W6Q5uf1wxbS5IR1p651mCiHHJdePn/wRjjTfMPAN7hj97+4ySUYc2iARMkzwOD+hRMCvyc3FL7n1OCFNc39/P8df+XjpIBn9Cdf2jTm+Ak1zZNvHaX1Hr3M7OftmnmGw4HA0/3R/m0/mFtATcpo1BUlWMcoXqp6RVMm58H5zHFjyAr2rXg+mZdut4z9E1mfdYcd/M+qiSPWM+O03mHfYnJVsyJpEZUaVBgyqHGku/V4SAR/GtUtV2ViKdfXfkXuduX3HcrBTvdZ84NgOcdX1SC+GG70OP9kqbld32rK1AyXVf0iq6swIp1rJuV1FQplARQToS8keI8QBv76PxasGmHfOJ1ttSmJEs/EhXp6R1E152EjWdTxCBf3dR+wFCmQaqqvbm820lzwjEmn2nameShnTST6ly62D0eMTbHB4+PKMpBIBS8j1l1TP+GjmOeYuWQXOuOLv+1cDWze07jjiji9/qGeNnHfhRMC75yxm03GjkJ6xAIyVAfIFJU+Pe74m5BCGo76aSz6SVYzyFTuoPQnVnC15QlMYhIv3ZKvQ/ExSIw+2CJViTM0J5WFU+l95R5nedyp5zbBCj0jMjkrnYiUHuFyekUwCqypsL8UnZe/OvQrAoc59NX/XyLpVyyuBBPNorM0o4d/7gYKc957JSKnTnOKToJFExUizqmZVdcOY1waqavKMJrB4lXsT+cOjr/HwK2VKmRFJ+BT3O7nvy7xRvlyN0cmahn/jCGuaCwVe6TuOf++5tqGLVkaKZb4kHN372Vbwh88Oy9xGCCcvaW59zVGgEo00D7ZXstQGheWc1XN1w1KBaptG8uu8Os1eJF+U07I3l0bbE94ecfKMQNPsROUZxf395Wumc+QlD1PwnWb6yRUccn719ISc5nCHWQvpRWYLoUTA3sQ0zSHbYzr54YS5yfI8h2YeHVmVZyRScs6PNIc6V+H/40d0e8Rp+ElNHJUizZVLzlXWNA8nRuE+rQhLPerTyIvm2Scza8hR4Xqus8M9nMIVT/xP/lOCnoyU/F/m3D68lXUI1eQZRpvww1uf5ctXT2+1GYlQMkKad9P+ee+V5cu1S6824zkJIXv6Zl8HwFd67mzocd9m40fHD12c74eX72rc1gYJR8zI9wd3lrA8I0yi+VqDrR9GPMzBCy7gxJ5/sN1b9zX0+6hT9ITuGJR9E8/RDG/XibxSqmlO+glLXKTZO5ajj7bqkgAAIABJREFUtmZiZAu+7T045AoOec9plnwy8qmwPCOTRqQZX2tc7Kz0JjW6WfjRdszQ0mGn+abRZ3PJqAtTTXZMmnCkuaDgUK5pDv+fpGQQUcrrNJdXeQgTffAXHkZ72stvM+GMO5j39tA7p+FkUig+uah1PO36+rX8YdR/stOaoY14WW0LVote3/fCW/zPI/Oqtu1vo/B+clQZz7pge/VS4PPZe4s/+su32i7I0QrMaR4hrM+1gfY1Yao+Am8Hp9lxihE1PwmikGPju75XsshQ+dj7Nw+c5sTKXw2DkhEBQ/KMpkSa334pubYSwL/fVqoQUIvolhlFjvU9G7pt5n1Nc3GpMQyUapoTdprjEwG9cytadzY2YqZBOzlPngEklhRUMmhCjOM5XIpRSVdSBcklAob3m+bLr1fZmGMojVrUcbyxYj3fuO5J1g82ft8I1/h2KkSal68t/sfU6jQPMdKcjUSaxwwuY17fF9knM4tbnloAwCOvLh2yHRqJNI/26jRnaujVxg66T4m3HBha9ZRq19lqyZAnXvUEP75tdtW24+rRv3vVDJ7tO5nJa+4GYLu37ua9mUjidjvcl1uMOc1GUynVNFe5KLTDyfnCHcXPvj2RyF0jDmTBgbG97g0oqdHJhkM4459QIuBGsg5e+HvM8gmuvH9lgo0NH8erUS0VhgquRVmdZvIMZF3NtuRK6zQD9DEAGi75l3D1jLhEQMePNJfOy8Q5zd4fyuCQLzjFqi8JnZ/HZ4vVRZKSfIQpSQSsITsactslTnP58RJXZSJfaE7w48ybZ3HHM4t45NXGZX1KONKsFChNBFRV9vzF3cEyJaXdEkzVqFBxrrKmORJqfnfOHab9lOxf6fFGOWykJrgbaS6uc6zjRl2dmAGEwqzv2RiAjfJDG1C5nutsowEM32kO34N3Wn4PAFsPznOX8Z+ylhjV+hycVmNOs9FUwiVyor3kJbph8Us7aJoHQ8PgzrnNfY9cpBp5DOk4GkTYAqd5TetGqA/rxzXfX3rHu/5z5csn6TW32VDDmvEGcIgZUroeojexHvIMZN1Ev0zBLTkXlmdkKJRUsdCEM9RjNc1OfAcwfqQ3d5n3Z95gMF/Aj8Vrg52KKKf0FDumkmJJK0eLNaFHkU+m5xcaFOje2a+XzY7rhORzzXmytGK9uy83Hjuq4TY0XJvdcYJIs/9E4sn5K0qWDyfcJTuyZXmk+ec9V/BlvSV26Wjy3Kq8e073ySBZz6EezA/d+VPVkrri43RtYE81Vva9E4AtBstzd6qvr/Y8/9Jx28w32Pfce+rWkvsjpJcs7zW6nE0qrz+h834kY06z0VRUYTSDjKU/6OXO7J3EG2zBJflPB8tVGqq1qYR72g9d4H2IXMAb6HjnHaVHvGQL32n++xkNGJgMTjjqku8vuTHELp9o9Yz2Ki3oiJfjnJA8owcncJqzhf6S0mcAGXVKosFOzGP+YRFbPSO+5JzEdFTDw+gW8oPFqFqDnYqqpFFtwnHYTha7w5eHI+tJdMpD2+aJl0odoiW6YXykOQXddhxBHd86nddPXzSNM28prSMdjjSrky9Wz/CuWb+975WS5dMa3KT8Lyhf6rmb72ZviF0+Wj1jTcE9p/soOs0DDTjN0UjzaNx9WUue4f+BoZY6rCeK7HdOvn/jM7y+bH3d/ysu0uyf2z2OP/hReVtr1rVXkKMVmNPcjix/DdYta7UVqaDAg6O/zXN9JwX1IEFZwmYMhsfMaYdIc0y9H40Mp91QpFk1KOcWZF638LGXEJFnRGz59f++EPqmw9Q0R4uotldiiWbdY7DSAA61CEduCo4iKPmMW6EzWygtOQde9NMJV5BIWNPs7av+UMGjipHmuHMunCg02F8c/r3B81NVSxKzwmRT0DS/b8GfeXD06WzfPztSc3b415fweXPgVkVnYoWOY5FuFiu9KjRpgAhf1/tUJBpciacXrOS6x0o1t04oEZB8cXv5zt9+79+8ZPlwzd9GrxH7/PIevvTfj5ZMi7ZUq+loIqB/X+ljsGIJt3qIVs8Yre6+zNaQ+/hlH4cqxauaCOi/extj0Dun6h1RtDiiYLh6huc0qyeTirknLVoWV4W4yHWPzedvs9p7AKvhYk5zO3LBRPTyf2m1FamgCluKeyH3E0cy6oBkGAg5zdIOhdSlXNMVvSg1pmlWerzLXpAs1NO60udOuCOQH0AjUaKL7n05+JytYwSsqmQjj4vbbJQpFS/RLYFBZ9zImzKYHQO4TjMReUaW0khzXELZcPCfGrypmxan+edWtORc3E29xLb+YnJYg+fn2X99jvf9+99iI6AS4zQ//MrbfPuGpxqWBG28xo2Gvqf/uSAREKisyVatOzm1ENInb90XKsOGMEgv2ZhtVBhszuNt3yn62R1zGm4jHGnGGSw6f97/6s268zcb557TufzwEwEXLF/PQy+XJulV0jRDfDKgVNAY91AInpI1cjyFR/MEGOWN8dZT48mcL3HqHeKTlGrX2fFrX+fZ0Sfx7vXPlUyvd7v7nYdwpNl36nu9SHNcPsRgjeP3zFtmceq1T9Zlw0jFnOZ2w3MiZPncFhuSDnGaZkFRyTAYLv/ttEG1kJhEiOgj7UbqkRYcLWby+5HmsZtW+UW6lGhq8/0sXVNZMtFDYXgVANrcaQ72eYPHXzTRVQBHenDIkCkMoigZKS6UwSEfcpQTjzR7N+xwh7QYRS/dj/GSqHCkeX0o0tyY03z1w/OA+CTgTIzT/KX/foxbZy5s2Anr790IgLGFVRRC0dKKkebHfgcXTYYFtUt8hp3mcOnEHgrktCe2GkmzIs1DeQJWyYEMV8+QwmBwLPmdKz+6efPBgxydvY9czLDiSRAdKnrDXDG5Me5/Rp3NYCRICsG8RvpgjmpJjkBfEGmufq0Qx480D+3crvbwcdsl9zNe+tl7zT9LpufqrGASVz3D7wz58gyNMaBtBh1rIeY0txvt4CymyJj+YsKbf+OUuEhzG1SViIs0R53kRu4NhRJN8yDPvrESHb9VjV+lh4Zu/jNfXcQfHqncYRtFfng3RIlcctossSSIrjUYaXZU+dfMdHrJk8s7ZHBQMuSkl4wziBNNBFSHe+e8WVx/wlp+Pzro11eGotMcvSnGVs8I7+twpLlBTXPczTqYF+M0x2kvh0Ih43bSFi9fVbJPK2rH3/Cc5aWvxM8Pkc+HZU1FecYocuSlJ7acpDbLaR7C9qq0bEmkuTAYSH38jr6fTPfuO77Ir3ovK+lMJ1qWMsL33irmf8TZHp3kO7rhEU4bsU4jmuY+8ZzMGk6zH2ke6uiv0c5CmIKfexFpM7w9qlUICeo0h5b3z/9e/7iNkXo4TdLktzPmNLcd6V1s2oGPvfzL4HM40kwmw0BoVPdMW2ia406P4VfPyDvKv+bc8j4DAwMc9ptpPPHKYndmptrI9ukQdp6WLF8VW3HBp4f88DLjozeONos0F7yOkjTYeR29YBqXjzqP7/bcSM4pFqMq0EPWGWTs2tfZWor5CllxePcmRWlOWk5zLuQ0ByXnIv9RYpzmcBRS8wOhRMDG7PS3R5wTHFfX2L+5N1IiDIpO8yjyJZ2VQqXOmt9RrqPTVAhtv558UZs/WvLkJT7SnE+4pGBl24bgNFeKNEs40pwrjmTpBTSiSWcliYAN7q9YOyLmvSO3IPgc55xHnc3AacYJLt+NOPVlkWZPnlEz0qyNOc3VdqFfGjPaMQvvg8GqTnNM59XbJr1VEgHzFmk2p7nt6KI6iP4Jm6GASKYkEbAtnObovijkE5FnjM0Xk3M2ltVc2ns+ube96G4rnjSE1jmawbJs8PBgED0UhjeqWdQZaddIc4P7Qby60xPkTe9RqSs9ymd6yTg5PvTyb0uWH02OsSFVUuJOsz8IQkykuVA2jHb5uksqqeT7Q9unMTszImzBCh74ry+VdZjiEgH9m3uuEUnQsrlsvOZVADaTVSUDNTj5+OvLek+Xm6swP0whFGkelS8d2dKR3thqCU6TSs4NxWmulDvmRCPNvszB6wwUHPepik94HzkJ5qRE/dvwMRl3LYou71/PsiF5RkN2hNrKkwmc5pqRZt9prql9jvysiq15rzMY7ZiVjNBYpZKGX8u6tMPklWRU32ku/18mz4Dmh7WMyqxfDtd9vtVWpErwWIlQpFkVlSyDWjwc46JeTSd60cj3lz3SbuQi/K3F/x583iPzovshSEpW12GLKyyfEuFBLuKc5nCEbhR5CsPRNJdFmpMf0GI4BJHmBo+/IALkyTN6UUDIyygymmN9b6kMZzQ5RueKA7wkPqhPjDyjYiJgbPWMkN491x9EmhutYy0CZ/Zey6GDD8Hzf42sPy7S7MkzGolcXjiJD3gfD88+UjIrnx9kdMxPnn5jNR8BZi1Yzu6TqzdfCG2D0YWI05zpjY00O+0oz6hH0+wMBsdCxsmxfO0gjiqXjzovWCbc0dA6Oh31Uh45LhIXtCgbHt77fW+CmuYBRjNO3NyPmiNMBpHmoXXEq+1C9a5TPZHOSbgjXD3S7L6HnWw/Eh7IM+L04ibPaE2kWUTmicgsEZkpIrUzLrqFmdfD649WnJ3ooBItIi/FRLBA04yDSIZcqA8XP6Rvk4k6EfmBMqe5EX3vewefr75Ak6Ov4ej5KHKlpbkoLZHWI8OMNEciuElqPFWVuW8Pr4TdE/M9B7bRREDPaR5FnrzjuDdZcTuLWWcwkAv4jGaQcYOhUdsS3vdBZExDTrPnmDuRm3jcORe+5kh+IIg0a4NPgkRgv8zT7pdCnjmZHVxbkOryjESLg1fRNHuSLKeOkfvCiYBjIk6zZnp5ty6EmddF1tucc3vrjfuK5/Ib1asZVNQ0V5BnSG4tB/30RgZzpceAhmquJ/k/4+o0+wxF05wNV89oyI5ijfUBKXa5alXPKCZQJhdp9mdFy9iFo/3haiZR4krO+Zp1P9Ic5+QX2mH8hBbTSnnGAao6SVVr9Oe7iMHSG/575Q1+ukeOLTdwT9BqPceRQiFTjDSfe+ezFBzvQpTJsJJxwbzYpKRmE3Wc8uvL5Rlp9GOaLVkIOUA7Oy/ym1EXlcwOR5p7yZc91h8S0Uh9gk7zFdPmcsCv7+PZNxofmnvtoH9Xbew/Op4mvZcCg3k3MqVkKGRG0aO5IELkM1ryJeUVky45FxdpzniRpIXLSgcqiDvnSnTPhVCkucHzMyPCpuI5mE6OAhmeHzuZvPSSdarIMxoYjKIalZxmyQzFaS7aNFZLt6Vf75tbT42stznn9gffsSEv9p3Aq33HwuUHwPzKwZhKEjOn5OlEMRHwX7IzeaLv6wwOlD4lGly/uvglhZrbPuEnYfVUzwhrmn0ddiOyuvDgJgOEneY6Nc1DdNWrLq1+RRD3OPYHbQknY9ajaS6JNAcjZvqRZqueEYdpmtuJ3tJavXeP/j7HzTqBr+z7HqCxoT/bjbA8o3+gn1eXrEFQRLIltWTb02keKHtUmEaW+O5n3cllD9TO3k+KcF3m7XRh2fxRocePvQlrmpOMNPsDOcxb2ni0efstNnA/NDp4h/ckpUfy5Ap+pFkoZHrp0cFSnahP6DhLekRA39HZYGzxJu9rl8/+6+ySZeMGXwgnFkl+oBh9rBVxeumfcNZG7kBNIUpGa3PybnURyZCTUbHVJnztZb2DNtRL5URAL5Jex6P0fOgY2YDSY06jpRU9miXPKCsXtvrNsmUmnHEHF9z1Ul3yjIzmyzS5OrC65PvA2mK0PcnjuNolthCznnJNs5fASIH+QXe/NnINU4Wzev8HgMFQpLm3ijxj0cr1LFnlRuClZpWNUqrdW4KOgHfOZqX8PKmWPDvJmc2UzKxIAMRt0x+0xUrOxdMqp1mBf4jIDBE5pUU2pMOd36+rxmcsFQa4GIcbxegEpzkszxhFjoG8UxzcREazZ//F/C7/SXfI22duhId/0zpjy+QZ/TGR5qFffHM1Ugl6yXPeP18ccrsNU6OD4ie8gB9pbtBpjmyrnGYTrZ4RF20ZKoEUpUGn2ddE94bkGYrgeJFmJ6YiS0bzrFf3vHASTwQsuNFhKU8EjGrX4zqq4SHWpTAQOP2xoweG8WUJb5ReC0tGZHMKbj4DGfIyip6q1TOS7ZxW0twGTwLq2A+Dg+4y+cxoNpDS2uaaiXeak67DXcK6ZcE5VvZovXdMqR3ecuff9WKVSHPxWJXCYFn1B2egtKMwEIo0b7zwgaHZ3iBxGtuorME/5HrEYX3Oc5obeGqbWV8cdCWXKTrN2SryjCvO/S5fWH6pu9wQNc3Vbi3+f/QlVd4DktJEwCr/8eLBH3LtqP8kvIi/f/1Is4Tudb/IfcFdrznNLXOap6jq7sAhwDdE5GPRBUTkFBGZLiLTlyxZUt5CO/LKPfD4ZfDfBzb2+2xv7OQjHzoCKC/xMxLJZ8JOc4GBfMEb3ETozWZ4i03oZ7Rbxufmk+EfP2ydsf5Fbv8fuO/5/jI9RiMO5COjpwDw0Hu/Ezu/V/KlEbmU8TsCvuMWZYwUb0zDcpojjtZ6RqMJJpb0ZP1KC8M4T5zGknaiFOUZgCfP6NVcySPv4joLQeWYpB/fC4qDICFnveg0l26nsKb59J4b+Vz23pKyalIYCKKPNUfsDAaJKV3H9oSeZIQizfnMqGICUrgZSWCfxlCosJ19zXZc1DvKgKfpLfSOZxNKo65lg/j47Sctv/F5+2U4d3uYfiVQu7xZ+BSuq+SckyuPNA+W6rhz64tDLO8w/ax6rK5ISanDKkKFuIh2eOl/yTzJ1lJ0dvtzjUeaM/3Lg8+5TDHA1SuVrxU/7L02+DxUeUY9kebeSKQ5fJ7Uc86URJq99fX9f/bOO06uqn7/73PvnbZ9N71AOgESIBBa6EU6qPQiKCAiWBEVAUWQjuIXlV5EEAELoCBBOqH3lhAI6aQn23dnp97y++PcPndmSxYEf/m8Xnll5ra9M3PuOc95zvN5Pm7b9/b929hNbvmcuR39N+K/ApotS64BW5a1AfgnsHPEMbdZlrWjZVk7Dhs27LO+xYFFT2vvx1SKMg9JKi8nDf8LTLMpPJb1DG02uaIsAIFQiNm0UtEKAYv/VsEXB+Q1TpD/63ms8MAxsIwS1mpjWTts98jdcT5r0Cy/30ykn0CQaa4V2YGD5hB7lyExqNrHmE23bIx8xF2SHHBFQDsD3WaaZduWiWGSaQ627Y/MzRBWEUOoFCz1U/FpNlE8KgpPhhEexP1M2A+1f3JN7PbAoKr6ipv0yjQ7ID0E3u7jAve1ZRRRMCRo7kWeMdhMsxly/zjwuuf5x1srMe3nLor1DofjFpEZOp1qEWzHQvNAc+B5+bTkGR3L5f8LZsv/w+03BHb899SX51kxi0zIBUtyi1AOTjHjgWY9Xt/rNStFhH1w9HHFUvcdP9i8M34tf4jf6L53mOYBTcJ8bVlXo1eFK4WCCbmuXgsD/Vz7Cw/Ff1l5bHELzchneUvlE/ZR3gv8loUKiYBO+J2QnElRIkLT7K6ObgLNnz1oFkJUCyFqndfAgcAHn/V9fBpRHMCDFIhe2K3/iURAn6b5LO1R8kWplTOFSkyzQU+Yjftv2ZI5A0+82r0PKzTAWLn+J50plo4hYijxVOT+mcpCPjvI7AG9bB9A813xXzOkWKp77ks88m5Q31qwNIgY9AYaDtM8IHsyO6wB2kN55zv2Vro7MFsILDUu7fp8zedxYycEspCKKVR0tMEvo23JioR+eYbjUlFqLVgK2P0OEYrpyTMifZr9YKAMaK7BkzEUCzkUO1FSVxKRvsauNdYg933+yYllWSxcn+anD8xFtScDUUmJ4SjaPs3myBkl+xQfaHbYTZDymw3duY12eSkJze5L7L6ypP2GwLofWPrJxqhiFwBKVLsMgeZ8j9cXrht9QF/uumz4tblh6OeXjUTJBZyPEF5JAcjamua+TBTueHEpH631JgL+pFhDje67K4VimXD1ZvDvH1Q87lvaY+ygLO6FkHESAWV7fZDzuCv+64A0rS8TA6F7z6OT/5Bw5Rne5y04oHmTPOO/wjSPAF4SQrwPvAHMtizr8f/CfQx6vL16IwFAL7RlvvjFB83FkNavkMugWNLLNqbK5lgMg+ZiUC/4mYXdaSzusL/3CMu5Ua9f0e/LalYRQ2iIWFXk/t/Ebuv3NTcqHKbZigbNKRHsKIfqpUlFfYlL//lu4L2Bgihmyhzd/3Daz8YwzSN1WW1soBUBHQQSE557BkJBaAniFMnqEgX+lHPQUYmho1gGBqpkcz6F4iYmwmWaC5ZKrBxojpgomD7w4i9zXeLpPP+fcNkQaFlk/+Hey5EXs90olpRnGEqsojxjo5JPfTFniPTB98tg/N2u4ix392UFxG4jWvWQkl3CJ8/IFnzflV5g5yueYd9r5/TjrvsQTj6M01eG5RmhPtQPmv3yjEDip2W5YEkzS8c2pRgEzSkrg271PZGyUvjBX1ij7CdVohIrncOjEvQ8prn39nT57I845Pcveht8z4IRIsj64sZRL+T3ZYVsCMtFxXwZZ3JeYjlXKs+4ec4Sxp8/O3LiafnaucM0JynYkj0fa223g0H3kf8CxmcOmi3LWmpZ1nb2v2mWZfUfdXxOQ030f/YZiF50aIU+2CB97iPcAebSCFvX6Mgz9HCi3H8NNMtO97xHpJOFWSwtbqL6dG59Dc3SMRUNNRENmgFmWvO8N/n0oDKygehe7w4G5eQZKYIDU8EcWOGVt5Ih+y2UANOxsTEYiYD75p4F+mmptuRZ+MepQJBplvIMOSEUWkKyz/aS/hvK9hTQiKEjLB0TFV2oWAMsGlI2LBMLgWKvXWRJELdkWwozcVrEZ/ZrmmNmzk16KpFnOLKANfbEqAzT7A+zp9V99g0lQbwCaB6MVbacFSOx/bEA6D4trB+cOMmQSh9As+lIcWpKQbMS84HmnPfsfmqJgG4mmLzvuvz6wO4N7cEVsXLyDH/ejMCigF2G3Ch9Tv2lwwGqyaEJu30MAFz5wbF/khR+mt9PzPTuvUIi4PbK4pJ9biJgL24sUR7J/omAEasO7OtPQnhftc2VjvJs9ILf89WPezUAHDnnNfa2XIS8U/jbuZNcKCyKhdxnL89Y9FRFa8TPS2yynBvEqI1v5NdZgd2qJfM/kQgYlYEtdZ+ivDzjvwWa7U6yB8kqFHOZ0s50AOViVYoYIkYsWV32mGN52ntz1Ri4JVr/vFGx9n347RaMWPxXALJWtLzIL88A0IyNY4cftpNK3jKnogwi0+wmAg6CPVm/5Bn3HCmZVqPoJi058gwnyVWNJYkLnXzRTsJTFHRLRRMGmAaGUDHQBlyeuuznwJLL2bbMQYLmPDfPWVIiAVIjLLH87hmakfOAdhhgO8yqw/w5z2wF8GRm2t1EQCdRMhwuFhzIRKgqCGZNFOJxOTEsFrw27WdaHXmGYvRhkuqsKtQ2lexq8LlpFAq+a31aoMPV4sv7P2XR9wO7VzcHJ/f+RyTjY8KDS/qmm6AaM0v74MML/wm890tvBjI5CEi/ffcR7nItH2wxIsgE5/Bvqo+V7EvaE8be2lOUfMOR5i390h1YWpAg67d8qKt3iVvFREDHpzn0HC5t9iYy4YlmlNe54pM++hM9i+kWd7Xt+PxFGCiYlujTZHLAce8xcOdBn971Byk2geZBjCgNVb+iAiuzrbLkfwI0h+fPRj5tAwuVpCbBchg0W+n10k7psw57AHIY2Hy+p+Q3GgijoloGpohRk4pmdoFSfWdrKWuy0WEvpTeulUuQZRMBQ/KMjQW6zxvbMj53H2usIZLRHaRET81mmjeqzLcdkSWly4bji5b1luwxKBpOcROBEpPyDKeKmqqqtjzDAEPHEBoG6qcgzzCle4b9PmvFiVt5rnn8oxI3BI1iyTKzn2nWTA80l5TRdpx/HLCUtZ/XCvkIViFjg2YVU0l4RRV8sVFltCMcLBIJOTEsFLzB349NnCTJSA1vKByNa6zaA82rE5O5qngiU6s9N41c1g+aP6XlbQcF25Ouar0juDtEPPgnCumc91sed+urmKbFuyvaEZZJ3gbNqi6f+U/M4e6x21hBW8wG4fdp7j9o7mtyor/d6vbkp2iYtKTlb+q04QPU0kqI48QGeV4vE+soRxGXaRZKSU2FfKGfv2t6fa+H9KUiYDwiD8GJsAQlarXGzzT7Ca18x3p3rHvd2hIQUqIxyJP6L2JsAs2DGRup46oEmofR+T/hnhGWN5j5jJsMlIzZTHPIPUPcfbi0U/osY+ET8MylAGRsBlbPZ0sSAUUfsuzDoVlFDEWjJlHer7lSZzhoYVevc5b6yiUChuUZ6kYyzU44LNZgAQlXEz9A/WtnxruPqEIfZcOxV9NzLjsU9zHNINDiSbnNBs0Th9dRREXDwDKlPMMQWu9Wbv0MYcszcH/jJComMQxX03x24Yc8Y2wvbfJCA6ufaa4yun0+1qH7VO224/yWznNeYZXIKmbtZ19gqPHINr9R8gwlaOFpCoWEvbpj5L027Gf0XNDch0RANzEsUetue3TSxdyjfpXqmDe05nK+7+AzYpoNEexbrBAj64DSiWIN6uo33O1Lm3u446WlHHnTK3Rk8hRsX33FllG1UVf6t4+5E716BI3CZ7s3ANDs/x2Cz3DYd9kHmm15xgUPzWPHy5+maJgVq7Q69zgQptmZmCiKggjVVMjl+7ka2ofvJ1wTwB9Onx1+ZvzEXTgRMAo/BJlji7wl202xewNgYVoChxQofAo5F1/E2ASaBzEqNfI+XqDsrgaR/p8AzeHP2ND5IcNpA6GQiEnwUeil+MdnEu/82X3pgMliLlealKLV9O16C5+EtLQO1NAxhQTNS8xRkYcnnM7wmvH9u+/+hBL8nsslAoblGepG6pA129e04KwoDNKSn5sIOED966WPfui+VuwJ8IbuHM8v7MUn3vkeixnwyzN0t+i0DZqLrqb5iqO2I5FIeppm8emAZsW2nPMmRhKym4/NAAAgAElEQVQEpcjxSOIiAI7eZQvGDK1HwyBXNFjf5QEs0wZhy8UYRhlrpAMARIBmZwJk/5YOoKzENOs5V55hqgniVrHk+fLcM/o/EUoXgmy4iSCeskFzwQ+afR/DXu5W+yDPcB02fAAqQ4K4piB84Dif8/6W+NRAs/192997UQSfZcXIBTQZDkB9NvETdnnuxMCxC9dLxjhX0F1N815ZKRfrVhvc43qsBAYKTDsKUk0M8XlVD0SeEWCaA4mAweOE771jOffA2zKBN1MwKvo619mVG3ubWEfKM3xMswg5H+Xz/cw56UOfV7mQiOOeESTq/CA6DJqjVqoDTDMWG6xGefWOVWDJVSr3emgoAyCJ/tdiE2gexLD8S7oDSeiJYKqLI7cHoJ6e/w15Rmhisf+SqwGwhELSBs3OkuB/NXwVtJz7MYrZkqzw7qZpvV9LL8B9x8I9skiNTASMUZPUOLLwK67X5fbXzK3cU+LCbj/Z/ica9jlCxXTKJwIGO8qEmd2oCdzDhtRnRzHNlmWx1UWP8+dXl/frmk/OX8df31wBDNxpIa8brLFLuTuM48l3vM437nyjskWVY+dW9ICJKixPniEEsUSKhNDdbP/66iTbjRsq5RmmYTPNsd79j/sZjjzECaeATbVv9WBYfYpksooERXJFk8tne368DtPcLhqpMbs8Jivcv2kO02y3FQcsV2CaRTHnyjNQEyREsaSPUyLKA/c1OjNhYCJIVclJbiHb49q+mZbFOLGOJHm38IuoBGrWfwgv/MYD1r6+oseME1cV2Ps8d1s+wDQH2/qghTNJsa/vL74BsMPqe+HSRu82KrRnZ/KpYEnJEFBjye+qW/NAc7XI00ktCIFS1RRkmgfASBqBRMDylnN+ptkoBvumbMGoaEJVjaNp7kWeEalptqVJiooSqrBYyOcj2eOSmgN2zP1kQ8W/L2+iQhss8yH9oDncR0dZ0KkhecYn1nAyVgKlZSHCMkpA82BP6r+IsQk0D2IEmOaBsGcRTLMQKmai7ovNNBtFePIiyHaUZ9OFQm1SMnZ5oqtpsbpUo/aphY89MlApWip6PlvSWUUyR61L4L7jwWGznM/c/DEAMYoYSpzqhEYXNdymH84cYzsuKJ5BlyUdNT5LeYYTYXnGHGM7oNRyrkZkXb/T/kSPleB2/VD393Uzsn22UQVDlrr95cPz+3XtM+95m5VtEpwM1D0jpiquZMFhHJe3yN+wp1ABzLryjGwACD307moULIRQ0GyHhZQufV+FoiDUGBo6wpRMs6monwJoNjGF4mma7d+42peopptALEVSFMgVDdfFBjwJQlqtp8ZKe5rm8ODpSCEc4ODIAT56RDq0RN2bm1iogJYgQSHgaQweaH5r+QBcakIsnIkgWSWZ5pc+WsW+185BN0wsE55PnMsdsWsru2cseU72ZfcfD89eTnXBLmbl6yvSZpy4psDQKaw/+FYACllP6+uXcw1qwZaQPMNf5jkqKiWZxVSnjmWQaQTIxoJJj87EV6luYmgANA9AnmH6QXMFphnLJTLCZbQzBb3iZMQhI3qbWPtB8yl/fJ3FG7rRbfcqTVNRQ0xzcuHDcPkwOaHyX6cMxGrtTEdu90cJ02xZ8PbdkO8OJtT7Pm8Cf1Kn7bBhP8+9yTOEPUlqseoQmRabaVb45h5SGlmwtF61/seoz/Ml5e1eP9sXOTaB5kEMv/n5gCo/RTzsSrYFK9lAveihoH9BLefm/xNe+QM89cuyoNlCMKpODj5lmebb9/207rA0QkxCjnik5ZyIWn5+4kJY+Dgse0G+d9lp+ftWkaOoVlMTl6CxmypOLf6MZdYot5OdKNaWsPKDXeABEXz8w/KMX+vHA3CW9u/A9np6KoPIMqFhBJI8i1apjdFgTAwHwkqCHFxc0GwPPglbm9qTL/N5TdMtydzd3RVoHx+t7UJgoakK8dqhANQaEvwJBMSSqMJCs/KYQsUSscFncuyB75lRZ/KOOZlnhSy+WoPXbg0TRLyKFHm6c7rLNAIYNqOc1RqIo1NlOySUgHvnvZ0k6sofOlbAHftH3pqi28VNhNSIOkx34Bh7wP/H26v63f7VksRsQVUiTt7S3NLwPXnDTXbcQ53vFX4Jg+YVr8mVomcvd0HykPxKuS8AmmMkbBcgx6mj4NNP+xnYQS1WFZJnLKndMfKwdE5+vkqPiJ9pJgSa8/GgI0nBKVaVaqRJeIVABtKO/UC1kk8zWBTtibezcnO08gJPxM8jky8yvutN/h7/FT0RcjOHjOiVafb9zRcXtXDxI/MxbGmVpmrEk0G70LGvXixfrAl60Zc4QdkRVckwHCWgecVrsjDK7J8EJSi+NhUlz3BchaLaW5hpRii0Uo+SaXGriZ57wBacuPPmtjyj8u96bexW7oj/ttfP9kWOTaB5ECMAqAYCmn2g+/2xJ8nL7HsRJGqpJveFlWc430tHZ2cJ6HSPEQoj623QbH0O5BmhRI88McxitsRdoLJG0fJYdnAnRSkrR1FNuYDAHw6zUyNy8PrNgX2FYmFwy/CGfouwPCNXhvGvE5mATVVfQ8MIFK5xtesbPnRXEQaDfRvoNTRVcZd+HabZAUAOaP54XXcQQD/2E3eZvqW9q2SAF1ik8wZaSpYVrrV6KFoqiqKgxeXAuwMfYwqVJDlmFN+DdYNXINVJtGtNjeOowqWkNQne/UzzNps1kkxVk6TI+q4cw4qr3H2OD66ekkv7DbZutQQUOc+BXSWuO+1jHTtXRt6bauZdTbOIJ0mKIrlQu/I/Ipli/0iDMNOsWTqaqpAjwSzlQw5Q3qI7X8T09btxQ95/CWh23HuaP3YT/6p1e5uW5Ibhl3JH6jS6ClCTlP1XPC77ECPn2YAJo4CGjooRaQE24HBQcL4T8t1k1ehci85O6aoRdoc4X7uPHYR0w3BAs8Byy4o7oaeCoLnoEByppkDCcABcFTLw5h2VkTr9kGdYpgvWnQnrb+O3MFVZRS6b4bBV17Gz8nFJaXOQCbry+pX7CG3JUwzHW93IFAy3OqamKVRV15Y5M+Q+UwZiWX3ox/UwsHZySZoXBKcyPjejhCgyjHYmidX8693VAMScwkZRTLMZZJpVRaHFqkPLOUyzQAi46qhtUGOJTe4ZbALNgxsBpnkAxSh8IGZ1jbTlUqYfhRKvoVrkyrNdn/NwdF2vLFpXvuqhUBnRG9P8WUaIaZagOe8mVJ1d+CHrrMboZchuX8W8Dx+Gt/4oX1smmCZVIl+2DKu/RCyLngzsS9y+J1w+nEGLEGjOEsoIt6JBcz09/PbJhZH7yoZpogoLw4oAzX89yV1F6Evp195ioIx8TIEmGxQ6nsGue4NuYZoWB/3uBc64+y3vpHf/4r7M5yIsCYHJI+rcUuyNoltqAwUBn+6CqbBZYal88+A34fVBqgppD3wOi6fZBZj8THMqHqO6ppaEKLKus4cjWv/k7jNtECxSclneSVAtKTPtsF2FHrAsGooRms3uYCVJ1cjbFQFVFHuSGnYhUHygLZPvH2hOKMG+xqkQVxBxZihLuD3+f6TzOqZPn520QbPQc/z9TR/Yd6RMpu7+lrV6u+yrFIWlQ/fmTvMI0rkiNQnZxmMJOQn1O3VYeoHFya/zaPzng8w0+671l6Pdksjh6OqQQDCs2T1Le5SHEpfI+1ad9FUr4IkMgD35c6LoMM1VQdlGADQ/dwXM/jEseLTiR/Bjaj+oFUaBa2O3MJoW975029UjDD5z2R7ySvmiUU777W1i3fSvk3kk8QvA4mjlBcx8j7sioaoa1XWNkee19QTHg3JMc18S6sxCCEM4uRNr3wv83pYvbyCOzkuJH/JM4qcsbelBN0xUVaBgUr3qhdL7CDHNsZhKq1VHLNvqMs3OM2gqsU8vkfULFJtA88ZG6xKY94B8vTHyjPn/wnzLG6xWbZB6OVURiEQ1tUqe9swXc5Zn2A+7sAywLDn7/vHHfJTY1j3GEoIh1bIjLAeaO2smffo360SYabZiUMy6uvVDtxlNwdJKNHWA7NRAThACHsSWO1M37HLiM8c1umwmBJkJMx/Uvamt/QSqvUVoid1US9n1qKgXPfzngz6U0l4wGy6ph641rk+s/5qFiOsPVJ7RSBcjkc/MQBMBd2qfjWqn5m/OWlj3QcDyzAE5ry5t9U5KeIxeMddTMimMKTC8Luke10APOipCCOI+0Nye833u5gXwn5+yfNmiAX0OfzjAx1mlStqJcFU+0IxQqLK3t3Z0ocW9dmDZgC9WHQQJWnjQd5nmtGQVw7HuA/jt1OA1rLwrz1BsjajfaQKCTHO6n6SBEmKaE6I0Sa4nrwdyUaot+czFrALnPTjX01g7unVTB7saXI3e7soEhtUmaEkX6MwWqbZlV7G4Y1XpMc2Oq8RWyorBzVHxJyivfJ1y9eS6uyXTXEnTrCmePCOsaa5Khop64Mkz/BHzl912WPpXbyj7NyHINPtBfd2alzhGfYErYn/EsmRybVEJyjOcKOSy7ngDkFODBaTiog/yDPs+Rop2dhEL+G38Fr6Zvs1jmlWN2oahkae+vXQ9G7p97jNlIFZffMBLmGZ/Dorvuyr4VjLiFEnYum0Fk1eXtqIpCmepjzBjzmlSl+8LLcw0qyrFZBOJQnuJe4apxjclArIJNG983Ly7ZIaMYkB6YPXXlusf30DJeNZWa9dLUCKEgHg1dUq+ZBb7uY3Vb0PbMvetac+2NUzSuYLM5q8dSXfCZ7cmFBpd0Cz/L4Qyj5/pjLZn+1QigmkWuqdpTsVVCsQCs/ySsIwS3bALVO1B+G9n7sq8Sw5imzH17DZpSAA0t7W1MiiR7Yh2cwkVFWloDA58h24f8sbe41zAs23qNV6xB8nmBRgF+T35JR/FkLWg7gOm/Y23EmfzWvL7jBXNAwbNU7q8Eq4qJtyyO0eaT3KO9gBFw4y+btwblIu5nkh5BgiIS1BaL9Iuw55IeYyYgVpSPn72q3MH9DkCf9/2aXZAc1W1vI9qEQTNii0V6e7uoii838i0tcmp2jBoDjPNPtC8ujQRqHXZ+4H3q8QoNEsnQQGQFRMBCrnyK3ThJMHeQi3jm6/7kuS6c3ogF8WxV2y0C3V0O4U//EyzJr+fKqPb9TGeNrqegmGypLmHITVym7C9q3MZb/LrrzQ4qExz6FkuxzRnbdAsQWnpMd9XH8IydWrw/PP9UVdTxbcL5/j+rg2iUh7TnLESDDFbPWDn/A4rX6/4Efw2c/4VJ8vuK1Vk/yAsC92x1AvpfvP5LKbPozqnBn2lnUS5iitavu/yqtjtAFQVWnjuIzkma5pKcng0gTNtyW2cceNs972faW6zvAm2Yhagay1sWBA4/we3e1UWw84gq7v8773vKt++1n3t1zTvKD7mlD++gaYItlPsVayMHFOcPsjPNKsYmEIlUTcCDZ1EsRMT4a32KLE+gf3/9dgEmjc2HHDcvTbwsOmFjdOepuqHM2W4/ZDFa6gm98UBzbfvB3+YAWveg3zatS1SMVjTnsaw7OWeZIPvJIXGKhs02wliVojlCPsFfxbxsTlW3hNxmfRn/8aJRJw8sWhtmsN06HmvDrAdlrOMbbsNaKpCXFN4+Lu7c/fpOweYCVHoPcO61zANuGYcPPrDkl3pbPDeq0M6vXMO3iZ4whYHY+xwKvVCguaKNmzgggvy3W5imB80F6wgSFyyIegQ0x9LLochfinxQ0ZkBsbQWkop8/0z/RbO0R6ioJvRGtS4951JdxXvmMaqGGDJiZMNrutElqF2wpTD+iKPCoBVKHUtGUg4ZbS/MWscADtvIdtzdYhpdiqcpdPd6IrvPmy9ZLwmBJqt0L05k8FlL8D795fcRzbvtbV3krsyu+orAMSEgaWoqDYru769k5//cx55O+nZzz7m+5kIHWaanTB8KyrpvO4yiOGYpcynO+c8rw7TbLiVBjV016Vip/He9zOq3p5w25aOeZ97Rlz3Xg+GFMmNEllQ9LVnvXw6WBamZXle8L74cewBdl98LR8kzyAuiqV9cCLJE+bO7vvpLJEvfExzs1UvwZvTfxX6NskuxzQ7IDjmVtk03WclvNJXyGcDk8+8GpRqJOwxpOKqhW+yNVGRQNlCYUOXHOs1VQ0UtPHHaNHGLzLXuO/9hbqOLVzMdcWjARjV/jb835Zw0y6B83+3yvPMNkLyjJa093v53TPMVq+/c61KgW9rUg6TiCnu2Gll5aTJISz8CYmKZWCiodZJCeDw7g8wUVy5DmocpY9FnwbVTvFzFptA82BF5+qAh69eGFgBiN/pR/G9mut4jN3ZapQ9S45XU0WO5vSnWPf904jb9oa/nuQCQQ2DRrrdZXmR8kCzpaguQ3PIjPH21mCH3RD7DN1D7E7p+IJM4ssRRxh514tb02IUbfaZfHfw3CE2C1HMlDDNphFa7rVDUQQxVUH4jq+yNq6IiPyDdif63n0lu9rSwaXw+vqGwPt4PKRpFgpqVSNNiizi0Svz50hcsu2YDtPs00mH0yA/XrWeQ37/ovu+v8vxTriuBv0NJVp/CJIVLEYlMvnkGWYhyDS/+50JUpJkrxaFI5nyzt2zbn0JaFYGkhcRCgXJNO84vonlVx/G2OEykcufCCgF1hJc5HqCTLOT3BUGzbEwaO6FgfID04xaSyHmMYAWiisJufGpD7j39RU8MV/a1PldFMLOGhVj7t9dJ5SSiHlAqieTDboeAWm7AuhY0VzKNFuGV/0QT+oxqj7FuCHyuluMsH9XG1zrvuXzlOGB5k9NngHlfXz1NBTSmGZo4uSLrVplIZOZyqISeYYak5/pJv3LAHxoysmYX9PciWzrxR47kW7MTPl/MqiHDkdZ9wyHYBBOwSAL3dZSF0PAcm1rB2u7vX5jbWKi+zovEsTRqY6rdGZLC+m4EWH7aCJcu0VNq1x8a5zwLBZ1H8RaYo3hD8aRAEzpeLHkPADFV7klLD2p6/zI9873XfkkTQmKbkGa/VXp5FGXjNFj56ssX77E/TwQrBSpYGAIlXjT5gA05lbZlpkeaFbNUku/na94mqNuejmwbYCLfV+I2ASaNzaqh8n/ezYgfNYven/LasZrSVeN5Qb9q3zIRLryJnUp++GM15CyMqxuz3zxZnDLnke3mde91HmcoM2h1h6wNZ9OUjMLxFSF5VcfxjkHTwegleBsvlb7DEGz6VStk51zzoqjGjkXNFtCxVDijO96E64aG8xicRjLYq4ENHf12B2cGt3xFhRPFqJZ0ZMkQ++HrswBBBGuJWGw0BACzYlY6B4tA5INKJZOijyZ3ryaneIp2Q4Mu3P2a5obRJBJ/2DZmsD7joFq+AfodWwq5QfDgm5GJw+pfilDNugOc4PP9ite6mbgSCIAUno3meSI4P5BAM0ipEusr5PA5Qfav3wHKe5EL7/2I1p8XZfDdqupYNuIhTTNVkh//8cJ/8e/jN28/f4KfKiYCQ9AWUIlUSf70a+YzwCeFMOvve2zPMOy4KFvld0dS/gKF2XSgURAgDmm9CfXMDzQ7HyHph4oCpTXvD7qZwdvCcDUkfaEwJ5QKflO+7MIUuanxDSbYdBc4dq5LgzLokpEt6+M4pvgCYFRPdJ9q8VT3HnqjvTs+XOOy1/EmUUp2aLGa7s9tiRi+YpP5AZnXExElOB2Yv2HNLx5HYcpr7Gb8kFACuXI+5wy7wILQ2gU0dDzwYn/E+99QsHy+ly/7VxRSRGnSH0qRtGwyk/Cwt8lsMwa6U3CwpK7UMTReXL+Ol5b2lqiaS5JrKwQRohFn/D6L93XooKmWQ+VUE/FVffe82nJNDs9gn+lVLMrkyY338Hd1uTz3hZaAg2dbOg53NCd550VHQFs0usq5Bc4NoHmjQxdkx2MlesKPGxG6GHuNcwiS4fuj45GXjdJ53RqEnbnXDMC1dKJF7u+kMmAYZs2J2K1w9zX9XlPl9VQV88DI37A+qMeZDXeMUofStsOWtiDjgM42qijobDO1QZbiupljgPkOrzXzqCqZ0sqY7V22MeJaFbz9ep93NexMkth6e7uyO2RUQFAWiGw0FDvAZkZuVu9ZTknjKLLFtXT0zuIcTrRbDuWT9N82DZSm/6iGZJ/hJZx2zMDkyeYA/RpNstkugOkc3p08pAPQMlkx4i2nt4QzWLHfImXRoEhJwUdM8RGlisHp7iJ97dHDGmIOEhA3RgAhogu3lvj9V1VNtM8ddyYwCmKmeeeV5e77zds8BUw2e8i5iZ24KfFs/h18Th5fLbN3W0KFcsvzRKC6s1kWzhLe5QpYpWrwfYPvn1mmnupRhf3yWJueHIeR974UmB/2pKgOkHRk2c4v6upB39Lnyzh0G1GsfDyQ5gw1AadNkh0bPo6qGaE6bmKDKqFaIhM2an9sfLH5joxTKss09wjPNBsoWB98yn3vVk/jv22HMFPD9qSEdvuz/knHih3+L6HcdNmAZBfZSdEu9UhK4yJfzqEoW9ey43xP3Bf/EoufXS+y8Q7k1kVg4Ju2tIEQVqpI5YPFr1JiGLQC96fm69VExe6W0K7M1umnUT0mQmK1GHfvz0Bbpn41cjT4xQ58563OeG217wKmsDpu8sckdtjXyvzJQSjsi2d93t3d3e7uTDnaX+jyvKVbcdkbH4pB6tvAqAUZVt07ksxvD5GsTXNo4YNYavcnSV/MZZIEKfI+q7o+/IP8xul1x9g//1ZxSbQvJGxOiOBRXNLS8DD0Mj1A9gUemTHYmvAurJFCobpVsijVs70R4h2VrcPwpL9ZxxmGc2gNmJL93V93mMZFUVwzNmXsf12O3Bh1cXu9pjxGX52FzQr/PM7u7GqfiZ1ZgfxjsX2TWp0677B02+n5bCPxVzJsnVHp2SdRBlWc9bXr+Ci8fcyv3qXyP0A6Z5+tK0yyVBQyjQP8SUC3nTGl7xlOSfMItiMo/Rq7gU0O0mSL/9OFrkATtptCrtOsr1e49XcaC/zAqS7OgKnf9ZMs6O9tyImNOu7MgFm0F1a9+ugixl34qwnGmHHb8rtHZ8ErtW931XyRf1mgXuOjZ3Bwvrd3U16oZ8T74hwbKOciMWibASFC3oa6Q74Gydt0BxL1tCmeh699aS55uE3ve8h5wMvQyaTKRgU0VhkSQ21kvOSWi1FRfPJPQwlTqpprPt+gljrXjcImvvINPciFanygeYqkQsAG/BcNiRottuSW3HPCJAj4dWBuM8Jh3gNuhJnO0UuiTeFVlYGtyJg8LuJWxUIhnwXplUeNDcXvcmcioHW6LXTWMoD1NefuD2HbztavhECJu4DQO2oyRiWwGiTz7wDmsOrEYEIAcSVbVkenRtcedJcplm6LfXEGqkqtgWOiRMEzRnfxCSeqiFBkZ0nSClJ1KR8ZVuG91e2lWw/TXuCH8UeQLcU97kd+rU/wk8Wlxwb91Xl03xyi18cthULLjuYbKhATLkVh0hnJjv8rG463eWC5knK2sBx1eS4uP1877wQLhF6wSUDVMvAFBqTh9cEVjzdz5WS1rfrOqPbjf9ZHYiPv3ehz3fu1ibQvJHRY9oZ0un2YCJgtqvcKaXx2k0AbNH8OADdtpbTA82SmRsh2lnVvvEDadloXVLeR7lSWBYseipyWQug6aWLI7dXjZkuOyGgqJZqPgEyKc8xo8HqLMtaDzjSG6KZKXsQMhFsN7aBbEJaDKmOw4miofu9ltN+Czb7HvVsBGjucM+PivHDarjs1MNRE9HFCQDS3f1oWxVm7Wbo9xrW5AGZ3SZHWCoZRbAZwnp6SpbpSsLHLNU9KpfLLS3FAVuNIBlTOP+QLWm1PHY7kw6CZv+glinodOX6mITSD6biwzVdvLrEtqpzwLKW5C79wMBxbR1dAZDTkbXvTfWD5pz7+Cw99H7Y/5cwbg846KrAtWpr7c9c1QRnvyJfH3QlAPnq0e5xZn9XqyLC0TRXDMsENUZeraZJdBPzJRO5yYiKyop6LwEsLgw+SJ7hMrHVRjc5pyjRqG3dMutdlnyujbQHRCyhkqzxgIOpxEEIFloSjGwmNrigeZK+iH0UyVgubeljYqzvmWtTS9uxP+G1znaJ8Mcem0nQmKTgtTl/mWqf9CGXiLYeA0BR6K6fylgh/YXXiCDALugmXbnixgEMJ8r0vZGR68QwTE7XHo/cvaTo/TZ1pEEI3jMn8aviKSS18qsxVNl6+VSSDTSidMkiOY5uVhj58qsAESsx7qTcfp5jDtOMiYVCLj6EOiPENFMMuPKk895vVZVKsd/ker43Q+5f3lKaoHjhP+dxxp9eK9kOMFa0sNoa6iU4qxrUDGPduC8HjosJg/til3OU8kKgKqWiCJIxFVEbbAePz11FVKjF8uRIvii/xxxxMj3d7oQ/HBPEOjTfiqUo2Eyz3YYTosDiZvlcOe4ZyZjKvWeUkjbJmibqyLCiLTqx07/Cl+2NUKkUg5EA/ynGJtDc31jyXACIuHXdezoQltf5lQXNhZ5SD2c7o/WOCb8LbI5kmjuCbOvazmz5Zab+xMo34Pod4K3SZZleY+lzcO8x8MK1kbvj7dFuBkPqq5mc/wvfKpzLX8f9KvIYS6ti99zvuUU/gqF00pWRn78zU3QfzKJhct1TC/sMqtwwdLh2Cjzy/ZJdRTtTf1RDNYoiMOJyqVXJ2oyZUBk+xBtcil2+Yg5uSeFsiSVSuksyzZWSzgDUZHnQnO1uL7uvJCowzWEbuqaGiKV7f9SPdeUZe6lzex/so7LmtRQj65MsuOwQDttmFGusIb7DOwOH+pnm/a59nm0vCRZ7KRdhBj0c76/scMHeoX94kRNvl4OkO/BYJjfFTw+c8+Dri7jhWY9Vcu/Np29U9JzHAAlFsvKnzYYRWwdvwK9VHjENLumEXc4EwKzzZBAV7Qz7GIpllq1K5ob9d/RkEw0iHWCaHXkGQkHUlVo+dmbyYFmkjG7+ZBzMdTvNgaaJ7oTKSQrT0y3uOZZQSdV7YNO02fqD8nJycVHsXnIF+dz8Mf9T7or/GhV97bgAACAASURBVIAbn1vSt5wOHzB7rOFE+N5bgd3CpyUfJjpKym23zLqAgqWSEEW6slGg2SfDi5V/TgGU2pFuDse9DWcF9unFPNte8iR7/XpO75+pt6ikYQ5HrpOa9a9zuBoNDv0ON5MMaR361cJl/Mk4xMuziQobNCuFNM3qcFI9siKd4W/H4aRpJyJ0ws7EyXIS8DDoyBRt73FBMTmEBrMzkFCZoBhwrKhN+VZWmhcQX/Y0W/19T6Yoq/h4fem9vLiohf3U98p+xE+sESXbhm53UMm23dQP+b/4La6LyTmF73j32DAycGz78nmRf0vN+zBELtg3vrNcTkLzSopcpqcsaP5L/MpA61YL3ZBpI+lbTVneIifnqu2eAbDrxCEwaf9ALkZd41BSosBbS3xSLOBYdQ4HKm+i+3JterIbIaVcOqfEiu/zFJtAc39i6Ry456vwkldbPWZXnTKynSRyns+ykSsDmq8cDbfuHdxWzELVEBYokwObXU1z7SgsNc7O2hJW+eQZS5vTzLrqWbb7Vd/ARMVot5eQP3m58nFRkbVBXPNHlY+zo7jPLwCoTcaYNXEIT5k7sqA7GXnsybuOYzXDWGUNRREWPcveAr3Adpc+yUG/kxWOnpi/jt8/s4jfPP5x/+7bkdPM/VvJrteWyEF+ZYd8+B17PDVnM2aKyvZTvGXLBUuXe5e1QVsu21PCNI9f9TAAQi21N/NHLOWxYR/MvJx5dV6bqV/4QMVzA1FBqmCG9mnxRJkj7RgxzZVnfF/7V+9sQq4DtjkWtjzc+5s+54Km6jjjJmzhvleLkmHYTKwHLC5+ZL67b11XPzrhChMFo5DjjVvP5oKbS39zZ+ARI7fhmfO+FCjnnqLA7Hne0me7Y/9oFFmqTWKBuZnU3Nt/u0TaAnDIb+T/I7cp3WeH0rC596aSBrSPITCweklccibxVfXDGJvIBpaXx4oN0i9dKFQ1SRZ8tW+i8/GqFihmUK0iHVYNbUX5nTksYafNNLvPDVKa1FTrrSwZtm1bIEmqO7jM7ERXrg+srE8mZwnNS9Z2wle4qIEeFBEEnLVjpxEXBt/RHmGtQ1I4EzFTD6zGmbHoFTInkpP39P5sdRP36fu57/UuCT5aBsMVqdLkOBzv/JnpT55UdnccPeApDHDpV6YxpiHFmIboSqaAC5rpaaYzPpLJ2ffh6UuCSfH3nwAPf6/03IjnxdXF2iSVJgxa0rIqqyUESu1whopO1vhWX7+ivsIw4a1Y7bfVKGi0/eZ9k9W9a9ewaH0a/vwV+Kc3mUloCtfY3sxuTNrffbntuFBbArQtDvI+eygUTO7X9+Vf5h7utmHjgpPo/7w5P7IvjRV8K2837RbYt72yGNMS6FoNxVxP2cqD9SKDZXnfbUxPwy3evSQpuKvXKnpwXDr5QbhwtfdZbMLk3UWfBCavv4ndxm3x6zCK3nNX6AmuGvYrHvxmiRXf5yk2geb+hKPJWvU2dKyElsWuP6yV62Lykj8D0GFVY+YqLDE0fwRXjvVmj3oOYlX05HWmj6lzdXH1KbsBx5KICXtznHiada3ejHNJcx+LTPQlnGIeA9FROkCoHIvgi6wVR9nrJ+77S748DSAwGfDHnlMkIzXfHA/AmAePgMuHcbL6FCva5L06mrC2/iSO9bS4DH8US5MrFDAtgZNn7CQuKfYEwRIais+rM9vhzb4zOXkfc5etxQwxzTPanwAILJlFRaJaMtsZK0HH1BPoqvOqqVW1fdjrx3OjAuvqZ2QNrcJg6A9fAldF3/Bcl9QxD9syIEUp+DLxhRCcd9pxPG/IypA1Iss0sZwXEz/i66qcCPZ79QBo0r3J639en8ufH/fsnfQNH/Mt7TGOaQsm3a3rzHmyjq/9neq45to0gSx/7Y9W57PrOQrEyBFDMXLSYg6i3VF2OVOyypvtXLrPDmXcru7rBt/nGGj0iWm2Ew5FVROTa4poGK5sqk5kWWKNBiGYPF7aiw0fPd49dUNrmztp7qDGXfXKFnR2Ht/EhUfJwS/pAwBFS7D5kCoytrNBVZXsPw6a5rF4sXRQz7qDkNUwm7v7MHnyTVRjsVip1VnSc3H4Wfzv/Eh7MLDbDwxbO+z+1vIBOF9/0RvTnJj1bfd1TU01fzP2cd/rndETg7c/aeu3+8Da9n6MBcuj7c7Y/hRA2hGaIjip//qs8bx8/n5oaoW2tJkNdKqGoqdscPnSdSjrfIVtVr4O795Tem4Fphkf09zcnUdYJkIoJBtGkhIFVjd7qxgHqG+zr+r9vURMkxKony0PXHtGqoUP1nRKIsznKz68KgJ8jpkJX5NERUNTKWimZhict5Q7J5SutCqYJQ4aW44fyzx7PAP4ljqblRGyS7XgI966ghKO/dV3KaIi4lXEjCy6VV6C5W9JSTMDXR4QrlJ0VrZloHsdSQoYsaBzSiDs50jPdDDhgsdYFGLq01lvDPevLP2vxSbQ3J9wOttMK/xuOtww0+1AhQ8wdltVEkBaVnmNcKEblttZ28UMxFKk8zo1CY3NGmWnHVgKsxvwbmvvZnlLDwvWdZGKVV7i71c4oNlhtzpXQUtpkkNUOF6ludXRy0z+SIkCqq827pThNRy1wxiuPGp65PFDahLccNL2fO+UE2wQK+Py2J+4O3Y1qzuyOIREui8slBO/mSTlKGWiSiNg1SXsRCnNSWhS1EDVwEKXB3AcC6tErpmudPRg1jx898jtTtSOkgzsQmss1UmNoZp3na5YBR1lOPooz1h02vzyx/nDB0BWlktKbV0CN+wkX4/c1gUx75sTSyotqlqMkWdI1reWDKOE/H73tXWsK9uCg8k9rwWT6qKiXveSzg58bC++/prHdDu+rpNEEJR9tK4LYRbppgpSjSiKIKd4rPhksTpw/LzVzoQ3T544eeKM0Vew7etyQmgmo5mn3mL42MlsnbuT98xJjDTWbbQtWdSgzeY+1irVCJMPsF83UaV3EEOnnVqcCaMj2xET94WGzYkdejV8WVZ7bO3ocCu9tVs1tPZI1jRTMJgwtJojdtoCwxI0Wp6kKG8qTBhaTasl+9O6Kgmebz1lR8xTpc421h0ECQ8lLgGkvVWv4ZNnDKmtKh38fRKY4bRxtBoEkf5VgnSHPfg7/biRD0xErQgrwUDEkrwSl8/6hNHDWaBuwRH5y+W5XaWgeeH6bo6++VWufKxvq3ZOZPOlk8uCT9u7Ze5PvV9klLTaqyHbqxdxZEzaF854Bnb+FoUR27mb453LKpwkIyr51i1mY3qgeUVbRlbZFILG4fJ3vOCeZ8tfWAiIV5WU+Z5SW+STVl/f0roETJNbTc/W7UHDZmQ32wUm7gt7nQcHXlH2T5100unc1vCjwDaV0knr+CHVgWdyP/U9Vq4NSh4AtHxltraIhkg1Ui/S5RzJgWCRsBqCfXZdzGDOwmaKT15MnhgfD/lS+QvZ32GT7QZz58vLA7u7M961jfQgVbP9HMYm0NxbFLNSm3nFaJhnL4tnvAah2kvcqYK3LU0Ske+EK0bBvcd61wqD0L+eBJfUy5LT8WrSeYOaRIxbT9mRY2eOZeJQX4d8xB8AyGe62efaORz8uxfJFnQOVN5ELVP5ql/haGydJazrpslJQR9iTZt8iLRsS7/tYhRF8H/HzWDmuKayxxy+7Wh2nDCEe439A9v3Vufy+gtPuIxEvysmhpe/3/2L1KxnO6jSTArE3PG2OpWkx0q4mmahxmCJ11nr3c3ekpXdJoYV1wSqgQX+dKoy8K2ZeTy364dyWfFkahIanTM93bXS0w8GsgLT7E8EHF5fZvD/6i3yf4eJVlTYUwLDVS1lOvWXrnMTI59oHYb1pV+RHncAXytciBqxDDt181HoWhWjRBt7KbJ09FZNsmtaFkrWuehfH/QqC4nbBSSe/nC9uxLk6K91uzSt45BwlvYohyuvsmRDGmHpAW3gesNjJL8zazhgcbF2N9uLRbz9iQ0C9Rx5YmStBGPwfhcrVbmQQ7kYWhPnjP2mYw2dwlRlFes6+qhrbl3iObg8cxmskYUNhGWUWulpPhnOsXd71ndVQ4jnmhlerZBKVeFwVCOG2ICjZhicM08y5bYu+NCFv4BHZLXJ1dZQmm1Qmy0apOIqCIEqLBK+5MJs0aIqrnFK8XxmGzvTPe4Ad58yWq46KF2reOrDUiDR3CfQ7PUDqpO0dcQf4MS/uZ+zbGx7fOBtT1crpml5havS62He3939VkTRmnC8vPUlfLvwI2rGTGPnCU0ss0ZiWIKGDm+ims7r3P7CUtch5NkFG8pdLjLiailsUjDJHfhrur9yF+9c+uWIs2DVHr4kVXvFsFZkEarGv8ddwPMzr+/XfTB2R1BjiOnHcHz+InfzXHMCrZbPdz+U12Nlgo4Vw2gvcS5JUmD+mk67YJBK3RApFxpOhRyPKJkUMK66GHQ6uX4HePk6ttblZOUNcyo/KZ7Fx0c+DlO+JFeO9vu5fAbKRDKmssesIBkSwygpEJOKq6y3giB+/YbS37suXzqpao+NcL2ni2jE64fTRHdQl7/VEXDmHPdtIEdBBL/3zWpVUh0Lic37K/fp+9FTO6Hs53Nyq84ZIQmNJc3BsS2d8VaBjE1M8/+nsf5DuGIk/OM0KPbAO3cDYPpAc60hH/bNi0vdbRlSjGt5QS57LvY8LsuC0LXvQf1mpPNFahIqk4fX8Jtjtws+1HWj6G7Yki2Er9rZoie5LX4d31UfjvaQ7U+4yWt9lGf4mHUHiGgYkkH3x8xTA2/XW70km5WJ+lSMWxOnlmwvbFhEwWYkNrrM+MPflZr1a8Yxbc1DFNCY/X2pSRzTkKKTamK2e0YsFocdvu6emtI7WG4zF8IetEdaLcRbS6UUG6wGKq1yAhCv4vfqqbxtTWV0Q4rRYzbnPn1fVllDGdXzoacj7y3KgOY3l7fR2i0B2bLtzqWppoyeecaJ8MO58KMPvG2Ncpl+4eLFnpuJaUD7clj6PHRJFrc7NoxvP7ya59saWHrAHaSpIqFFfHAhYPjWnKY9wdc1+bwMi8vv0J9858TazspAMmn00Nyd54w/ewlgjk2SU94+SYGiYXK+dh83xK9nSXMPiqkHtIHDNvf01lPfvoRj1ec5TXuCG6rv4P2VHWSzOVjzDnohS8FXrvcm/cuIMj7cvYUQgnMPnErNlD0ZJjpp+aS0/fzg/nfZ9cpnghuv3wF+u6VcPXjxWrhtX7ndMkuBg0/T6y/OQqoRoefYP/8MNVW9yHWqZandLYsfYQgVE8F8awLN3XkM06Inr1MVj/4ORiblhGW5NYrvFs8hXuWbYMSryWgN7JN7lgP+vkXJuRu68qxozfCBzfSv78rxswfmBidSftAcs2UGM78BUw+Wr8fsAE2ToGFc8OLH3QNH2bKdJllJLmV0s6ItQ75c0mulgh12/PDwHfnyCd9ml4lDmDVpCGmqSJNivw13U4ecFF4x+yOueOwj/v2+fHZ6+lsNM2JFScMkudu3qd3+SKriGtYZz5QcU2iYDBP2gmPudIux7KIsoD6/liNOO5+9j/h6yTl9ieljG3jd2oonDTnmjRKtPG967DM9PlDVsrjE9u/N5HdJ28ylk0xXI3Isa8mgoWMocaiWxMNoUWoR54Zf9nHyQ+7LpJHmJweG2tdijwTZYDVgoZBt2qrXz+qP+i2CoLlWZEtXeoC3trmYO/WD3fexlS+VHDO8uKpkpbpbbaTaBr6NIk1VwwiGqd1BB5hUI4zeHg6RCbSJiAmVEzVkeDx2HgB3Gwe6dnyRMUxKBPfq+Bd/b7qF95at54rZXv/U3u0RHEZHP6uyfoGKtm0CzZVig90gFj0R2Kz4ClnUmcGs1vOLZxCrGgDLVD2UnrxBTbL8slh87Pbsrc7lfO1+BCYbNkhmaZKyZuMBo7NUH/bSjGrMi56SVfCWyWVN028j9M6fvdeKBod6Oq/Vow8i/a3XB3yL208eyxyCE4/uzjbXCmyjvwNfxI0eCsTY3C6NO7YpxTqrye3AlVQdbHOMe3wT3TzwtuwoYnqaZeYIFGFRv+Ylmq06Lih+k9vjUjP4rDEDVen90fvXd3fn1lNmUp3QGN1QxZvbXMLvk2dL8/oFs/v2QaLK687+MVfdejevLpJMnrGrzz3kh+/DT5cEz2kc5w5QgAuY1Ewzzy+y2dWHzoTfbwd//jIseQYm7MWvp94PCFa0ZVxQkyoDpLSpwQx0x25pfUQC4NoyPqFO1Ike3l60iptj17nb1rdJfaBue58mKLLGx+Iu2dCNMHV035L2kM08H3GQCS8AtTFZmGLhs7KtzzUn8X7tXu5xv9ZPKEdw9TnqpsrrZRe/ULLvkffXsK4rF6F7tXzJTnKfsMzSpe/tT/Ze+xN//I4uWtJL4Axl7gOBZEY138GN+leoSWi0Z4pc99RCTAu2HGUDyp2/HTh1V1vWPqretnYLycy6GrdmihKUwwDUJjRWtme48fdXcM2NN2FZFjfPWcLf3lrpPntA0BUmiglONcIP3oFz5sJuP/C2+3IUOPoOAOpFD++ubCdbKKOtryl1UwhHXFM4bNtRCCHYxQYlb5iybe2oyORlh7Vz2navHuihcOy+/FrZcIixO/JxTVBPr8WT8I1/w/SjYexO/fqblWJ0Q4qbv7YDfzHkcn+BGJ+Yvu8q4wPNXaW/NYDeIbf7LSRb0nlZylnR3H5o29oKuTR+0Dx5f5lTMHEf6FrL6dsHx2ljpTc23aR9gy1H1jJleC/ym1CMbkhhhPTFQ2tT3HJyUAp45qG7smzHi8h9T66sFbtKV1WqrCyEVhWTenB1T1QPpc5Kk8DXPp28E3tCFzNtdvnIYB7Hx417Q7uUzjQ3zOCG7x3L/ltVaM9aAo7+IwA7Z17gWPV5bn/Rk96k13kER23bByWnV4wBeuv/N2ITaK4UvenV7OiwvI55qTkKdWxIK5ttl6x1pdjyMNI5nepEedCcmH4EAGdp/+Z87X5eWy4HswRFmjc2A9tptJmQFinfBaveCtqHLbQnER8/BqbJ+Jd/5u17Ujpj8KMP+eC4V/ikwwOy7Xv+ikljg3Y7/Ykjtx/Nb/JHBbYV0m1ulnW2aPRNA2r07QHNWzHXUnCrkXW0+5YXtaT9+svXQ9NERmnd3P3KcjZ09JDUu1zPWQWTF81tud/Yn/fGnc6h+Su5TD8lUqYQjsnDazhomvy+FEVw3fEzmDzry+SsGPk1feyU/EyzZckVgjfv4KHEJdQJyYzXpHzMY+P4IECOCnv/MKWTq+96UDK/H4QcPUZtR22N/I7uenk5Ly2Wg2QyVqbL8U1AAOhYwetDfkVXTidXNEiSp8auyOUHzfNWSe9uv25whrKUVx/8PYfYVbAAWtvk3298+TJAyjPeX+kNQCM3vIhqFQIlaFNbBb2anajLrmRaVSfbvflTiiLGNfrxzKvalcWxqVyvyyphFf1s+xDDx0uN/24fXgpv3+XtyKepRwKsDU5SXI/vmQ25tThL2YHY6nAYYecQ+CfFfoApFDjgUvl6j3NLbzDVEEiues+cxNd2lc4fNzwnB89dHdZqT+98q3ECib2l7vO642cwc1wj44Z4LD2AtV10xbR54jj+/OonXKPcwD3xq1nSnOaJ+ZI4WO7Tp1qGrzSwFlXMxRcHXgb1tmOJz9nF0T3fGL8e46GzMeYG2/dPi2fyhjkVtbp/2vWZ45r4x1mzWLaL1DXPVBYCFm8sk2yp402ersQ0p5vhsZ/K4kl2OHkUPy9+s+Lf33zKtoH3mw31MeV1o73vYJez+/JxKsYh24xixr7HsFf+Or6av5S/G/uQjtt9ix8MlpG4FFpsQOZrowkKaBRloqLdD31pdIWxL8oXOtkA6+eh/XZSYLNqu668o83gsYu/xuPn7FVxPI4KIQSHFa7iYV8J+WH1VRw8PWjXOLQmwWVfnU5yyOboaJzUeUfkqmBx9fuB9zVGpyeVAxi5DQKLlCjQ43y3Uw+V/6d8q7oHXw1TPBnULfrhxCb4HD2++zjTx/SB7PP1019S3uZAxetjt1l0s/t6Rutj0LqEtS2t/baJ9LZ9PoH0JtBcKfpYxta/tNROLRO/cgGrLB/wuGY83DzLfXtQ/mom5P7Cn3XZiNtG7cXjuWkUDJO6ZAUrsslfkjNs4NvabK6Py2ScJAUeuPlirrny59GsUF/ib/ZAlQ9Z5TV/DHfsL63yLqmHa7eAN21LntdugksbiRVCf3Prr2LVjebwuxaz92/m8OLhL/DV/KVYtaPZmNh36nDGbr0rB+R/7W6rKbYGluujmMmSiPpdO0vZjs2UZmK2jqKxOs5zQz2bpljC7rh2+DpsfzI1Zhdvcgo1v5PLuovqvU5zrdXEybtuzi8O34oPrfH0kCotUd3H2HpMI3OtiegfP9m3ggZ+7WCmLVDu+0xNstUN1b1YzYWjSoKhaxoe4YnE+Xx80wmlx0w7kltfkJKlpS09XG/LLMKsohtNEyXL7YsRPR9TbXZj3LYfC5Kn8UHyDACu/s8CTNNi0fpujrjhJa5+fAFdojZw7q9idwfed7RKJifR6iVYfflfnvXTH8wrmdn1DGMtX5EaOzGKVCMMCy7TzjYlqHhen46FwlUnzmL4j1/mt7osGz20thew1ksIRWFJlfz7y1+8n5vn2Oz/rXvyfvJMpooVXnVQP1Pn94Be+ITNNEd089PtyacfsOzwDRhrM5FNE2DIJPjFBpgWXS7Yn1z1vjmZfacO57mf7MMh00dyyq7jGF5nT8ZqR2JO+hJ6w0TED96V1oVIL9gHz96N2lCfN2Kno+i2ouUhMZ8d3s8enOdOoPwJo/mc1+bLSUQCse8F8v9Gn1yjdiSocZIUOEZ9gaHL/w3AR+ZmHJK/in8Y+3Bc4WKS8f4nzO00voktJkt70e9qj/ATzdNIr2jrgzzuuSvgjdvgw3+5myy76upKq7zmFiBVG1x+F36LQ4DRNuEz9WAGIw7dZiQrrBE004jWNJZDu38ud7zkq0lQxmO6NreGFa2ZwP6txAo0DOntrSUgWc/4ZX8FwJxS6pfMmAhJ5BaVP9ukCRMrf6he4sqzT+CtHX5NWpPPR5Q8ww0h+HDkVwDIvu25eLxpboFuKXR//DwAPcjnocdKwpnPSaeTH74PE/aGuOz7qnY8Cc5fCeNsrOG4mQDseLrbbwOc+fObmHj4j2W/O2pGSYJ2xdhfFivbV32f2+Leat7kjMyjeMeynZ6u34Glvz+CW55fWnKJkoiqAtjRe9L3fyM2geZysXYu/OPUyF2/14/iBWMbbtC/wrvmZC7XT2aBKZnFFdZwUrWN9Jz9HnvmrwucZ808jReOeJmPrc3ZYkQ9v9RPY0LuL/x22BWc9Zd3ADhkegUmVksgLlxTsnmsaOZi9U/8rHAD1s27RZzYz8h1eiz7Hw8I7kvby0g+EJFPDudbhXN5xLAf1gMuDTAlT6xSeM+aHHDNGEgIIThr70ksssby9HEL6Rm6LVOVVcx9/TlOUZ9kGB2uUXvFiCrnesf+pYdZscA9H3DIkbxubskjxixpZeTEzNMA6QxSZZemjc/wEkBfsbbjkiOmMao+xUm7bM7EYdXsM3V4Hz91MHYa38T9HER11xK463B47krv85hGaeEcv0Z9/bxILXRZIFsuquSEcEiPLFqzT34OACcUfsEfa8/Cmnk6jJkZaZvVWFUBTDaOh31/Hlj6fzXxfaqbvWID24ilJNKr+Mvrn7gA4/EP1tFF5YSs3NqPvLLeocgM91aGXjB8HsqxlEx8O3cBfPc1yazOCnrMnl38EcNrE4xpSFGXjHHXaTtx9A5jqRoAmApH9sSHyFsa4zteo/2p3/Dde9+BNjkAPZE432NXfUA57/fEve84JuqLI50J2ONcWfTDX3glXgWnPyFZKTvxOJA0GBUn/pX1O1/IZV/bl10mNDFhaDU3nzyTy74adMNRTnkQ7Zx3yyZmBY5NVPHwfk/zqCEH/fWzLnKXhb9f85x7XOPKp/l3/EKa6HLlDR88eBXJv0k27O/63miT96PXmHES/LLNTXRyY5ezSg69a9qf+MjywHVigCsKe0wZxi+Kst/4nvYwI2llslhFtm0dl2h3kSRPOq/zybN30L1mYfRFcp3QvBAWP+Mm9aap4tLiKRX+sM/d4dTHAkAKgK/cIPXNE0K1BAYYU0d4k9mdxw/xQP3yF+W9Q6DP6rSqWH6CBIoHKG/zwFufBCRmR2ovEUPHcsrXb/0Vd58yYS+4cC384F35zJ7zAUw7svSmtj0O6saWbm+ahKUlqT/qutJ9/YgdNm/ksiO3YcNkOYFer1Tu6/WDf81ccwL5py6jiMZH5uaYx9zNB9YEOuY/jW6YCAFtVg03jvsDDN9K/k6N46Wj11ZyBVqosYCdIlVN8vfe5wLvOf7JYvjZJyipeinHOvN5+cz3J/Y81x3zouLF7TxSa3d1Pi++7EnMCroZnchtr3T7nT6Ka94vPe5zEBvfs/+vRqdnebSsejuGphdSK7JkrTidO/+Ir78iNXTXAqoieK24NQ2kaaiTncTUkbXc9N2jOeTeoYzMLOG4reI80Lozz7wsZ083fm17HnxnNTfPWcK9b8hr/fSgqUwc1oskREvIDuH/PL3lZMUD0sJ330ub07y4qIWvzxoXXWyhfbnULoYHi6tD7MMWh8ilu8n7S8ePQ37Di01HMn1YnMZYkTnLdZ66522eMnek6rjbmaYMZeEKj9H8y2sr5K0PkF31x/abNzL3kgOpS8awVu7Nri3X80hCZml/U/0PT6+bxR5ThtKT19FUET2oPXB66Ta7mMKFk/7JffOzzBCLaYmNwp+esfcWw3j2pEdYsqKDw4f6QFpVE5z4V2ncD6w9+mFO33o6cC7pdYu5YJ8zXH/TK48sX9iiL5GMqXRMOIJnVr7N/itehhWvSKeERU/KZebOFfDz9Z4jgp95fO1mb+luYyJRK5f8fYPZbGNnXjO35rVmeFCrQ71efnPf3msiJ+y8OSfc9ipH2TiKswAAFCVJREFUbDua0ZWKIwDsfZ4EdLfvC+vmukkvTvw7IeU/Ux++i1lT5RJ6azpP3ucFXRy3D7FP5gTOO2P1RXCFbCfLzRFc1HgNZ7ZdSzdVHHrW4zz99nzOeHAFp+8+kb38J/qZuFQjHHSFfG4WPIqx53ncNGpXV5sLsM/U4QOeEIVj+mZDeX3K99ll8XVcGLuf8fMO50afkubjVRuwdhhDMZdxa7itXvIhYZ7McRAJhBAwdErpdkWBXfuxND/1EEZMPYRD+n5Gn+LkvaeT3/0/NH/0HCO23lvmRzx9CT/o9KzT7ojLIlPvJM/i1Nbz+OiX32O64mmbW6adxjab9THxOKpC567fga7VWNsez/XPLOSmNZN56bBtmbVFC4oieGVxC1NG9E/z6oSmKhx/9q+Y9+BKtml/mteSMq/gJv3LnKo9yWprKG8s2on9Xvgx6eeTtJ+3ksZq+StbsSoJLf5znnu9iTbYMBHcaRzCv41ZHLV1HReE/3C8Ci7ukMvhUdKVpgny3yCFEIJ/nDWL1nQB3TR58J1V/Kp4ChfH7oEbd5IaY5tlfNGYzreL5/LC2K1g4r4cuPQ5tDfPg52lo4klVL5h+7e7pbKP+IOUK7xxK2x3ovx8Tb0wxYoK586Hp34JL/9egspELez6nehxcoCx2dGX85/Z27Hr3hHA3RczNh/C1UNO58J22T+tGLIHB227NbfN3pMz83fRedUU6q0sz6cO42dfO6z0AodcI5P0bPAciC9dEnwfdv9I9p7IGhn7/QLqx2K+/1eU1kUs2+5ccu8+8P/au/PoqOurj+Pvmz0ESFiCaNgFBDegIkLRB21V1FqtFau1YtUKxwqttS5PW1sf9EHrOT21atXaHtRW+2gVqRUrLS5FRUEFEdlRVmUnECB7ZrnPHzOELWFgmGTIzOd1jufM/H6/SW64TnLn+7vf75dXQ8MY/fUz+HB+P87IiOzqd0nda2wvv5LFmyp5fMZKZq/axpoH9vo5Ni2M5AGYELiWd8OnMiP3NrKnXA9Z2Q3/XElkh9RvkuhvanYB8DCQCUxy9wcOdv3gwYN97ty5B7sk8cIhdv31Wt77fAs3B34COO0op61V8cztV/HTFz9laK/2PDZjJcVtcuuXQnr3jnPqJ48BfLGtirunLuLt5Xt6uC4bVMLvrhwIRG4r3vPqEopaZXP/Zafsu2LGwcx6FN57MNJ7NmPigeez8vl9wTimbDmWSeeE6f3VyyKtF+FQ5JbsF7MjO+8A9DkfPn+d+QXD6Zu1mVY7VwAGF/0mUhCefPmef5awU1pRy5D736JTm1w+/MXXmbZwE+Oei4yU9+vchnVl1fUjzaOHdq9fW/et20ZwfKwPBYejYitfPn0dXbftO/O4IqcTTwQupP2xPbnh3EGRHaHWfxxdUu7AGeS7nVTzZP1tsN32eXPHEg5FRnJj9QQfodcWbOTnz81kQd6YBs/PbPV12n3vaU7O3YK//wj2yTPsOvartN04q/6aR4vu4Lhts5kX7sPE+x5q8OscVPlm+MdNkYlUhV1Yn1lCfk4W9766mH/M3/Mh7v9uPIPhveP899iylLonL6K8JsCvAtfzeM4j9adKvS0TAt/n2qzXGVd3C8/l3BeZPHbipfDtSTCx8dvUl9bey903jeaqP33AdwZ35b7D/SATrIvcMs7Oi31tAlQ/Moz8Rja0qSaX8qL+dNrR+Na/b2aP4Ny7pjZVeM1n7Wx4+tDaBuaFe5Nz3cuc3Ktb7IsPQTjsVNQFD94+F4+ty+GxAze7WerdeSo4sn7yaa1nEz73HvLXz4Jl/2zwS70ROo0xgdsA+OU3+nPjWUfWZpBoFbVBbnh6DgVUMmnjKDLNWT3wTnqUzsDWzeGbtRNZ6L1Ycd+FZG1bDo8P3ef12/tfQ/ulfwXg1W538s0b7jryoKq2HzjangSVNQHe+e3V/Ffdu/ylz8OMu+YqFq3ZxLInxzAqMzJS+3LJ7Vw25lcxvlIzc4/8zcsr5MG3VnJ8cQGXDiyhtLyGt196jFFr762/9OXQcE61VdSSwzFfvZrWbYvILTxmnzv6dwTGMjl0NhdkfMQTOQ8RLOxB1o/n7jthuZmY2cfuPviA481dNFtkLabPgPOAdcAc4Lvu3uhMuWQUzevKqrh/2lKmLdzT4zigaxG1gRDTfnwWGdHb9pNmrqJf57a0ycti1spt/PDs4xv8eovW7+S9FaVcfOqxdGnXqsFr4jbvWZg6nsrcThTUHt76nrsFc4s4f9cv2eztWHxzV17Z3oXnP/qCn13Yn3lry2idm0VtKMx9ry2hJrCnxyw/OxPHqQmEufwrXZgyb89It1nkQ8Sq0kqenb2WR68edPitAIdgw9rPaJWdSf4zI/fZyjyW8dzJstqOrPASrhzcjRfm7rtMztBe7fnb2GGNvDp5QmFn7DNzmblsPXdl/ZXvZ71xwDWr/Th67rWBx5m1D/Fe7k/qn7900XxWl9XSq2NrLj+tgVuVcXKPTGoqaZfPjqoAJx3X9ohHbz5ctY0TOrehKC8L3vwfmPVI4xdPiPbXV5fhW5djFVsIlwzm7Rd+R2DDQm6qGUfPjm1467YR1IXCZGdk1L+Xj1qBGrhvz6z2zZ3PhrbHccxnz+1z2XZvTXuLtCk8FzyHd8MD2EFrTh80iNuuOMimBS3J52/A9LtgyBiCGxZiGz8hs3RZ/WjlxtN/TruzxjBrQ5BzTuiU0JHDJrNrY2Qi7eJ/wPpD/zv3cPDb3JIVWUbt84zj+Un1DVx1ycWMOq1ro6vUHA3cnUnT5/K12dfWbzIUduPlEdO4ZMTQ+nkkwboann3wdq6viewguPLb/+b4vicxf3MtPYuLKGzV/IVUU1pTWsnEqfO5Znif+rtVdcEwj03+FxsWvc2Ai8dxzbAeyQ3yMPn21az+8xh67ZoT+2JgbN2tvB6OrOIyNGMJPVoHmXD7HXHNHThSR1PRPAyY4O4jo89/DuDuv27sNckomifP/ZI7XlpAm7wsPvrFuWzaVUPPjgW4+9H5izhQA9l5vPOvF+n64T2UZnRgSPhTduR0ZktNJmZGdU57gnkd+Ur5nvUo/148jozuQ/nZh9n1xfBJx7Vl8YZdjX0nAI4vLmBIzw48/1Gk9aJdq2ze/9nXWLBuJ2WVdZzVt5jtFXX7jLo3uXCYqvJtzHryDrzsC07PWEY+dSzK7E8320xxaDO/DFxPr7wKXqoaxBLvAUCXdvms31GNO9x98YlcM7Q7lbVB8nMym6TIT4Rw2Jm5opTpizYybMVvyQgFqKyupluHAoaWvbrPtbfU3cwr4TMpbpNLTXkZFeQx+abhDO6R/BGWuNSWwx+GRyaKZOVRWTyIgo2zYcjYyN2RRtQEQny8toxB3YoS0nPc7DYviaw60CvSb1pTXcV7i9cwdNpIanKKKLvmLR54fjpbMjvz6yuHUFEbZFC3IrIyMo54PsFRb/3HkN8+oa0FSbFsWqRt4rPp7Nq1k/zug8g+ZRQfzXqTjfP+RZ+q+UwJncnM8KlszuvFsNr32eZtmeP9yMvOYOGEkfVF59FuV3UdL/7nA96Z8wk7Mzow6dZRdGqz792bnVUB/v7GDEpWvcCZN/2eVvnN+PdEEmbnzjJCr93J+u3lbKEjqzeV0sfWMzhjOQVWyx+D32CrF/FU6EJ6H9OWkSd1ZuqnG7jt/BO4ZMCRLSAQr6OpaB4FXODuN0afjwbOcPfxjb0mGUVzeU2AdWXVtMnLSvzIcDOb90UZ0xZsZPnmclZuqaCyLkRO7XZy8ltTa3mURper69mxgO+d0Y03l25m084aendqg1mkqLx0YAk1gRCtcjLp2bGAgpysyFbDgRBhdwIhpzD/6PnkHwiFqQ2GmfLxOt5bUcrKrRW0zs2iojZIZW2QwvxshvXqwJl9ijm3fyeWby5n3fZqvtav09E/8tiIYCiMh4Nkv/JD6gaM5sVtPVmwbge5WZl0a9+KRRt2Egw5J5cUctOIXkfnhz85Iu5OKOz1PfSSWnZU1fHJlzs4paSQwvxs5q0t49jCfLZX1dE6N4veh7mu8NGgNhjCsENvTZQWz90Jhp3ZK7exaks5rfNzWFNayfXDe9ChsY22mtnRVDRfAYzcr2ge4u4/2u+6scBYgG7dup22du3RufxIKthZFWBndYCu7fNVSImIiEhaa6xoTsZHu3VA172edwEOWEfN3f/k7oPdfXBx8cHXnpQjU9gqsvOdCmYRERGRhiWjaJ4D9DGznmaWA1wFpMDUbhERERFJVc0+G8bdg2Y2HphOZMm5p9x9cXPHISIiIiJyqJIyhdzdpwHTkvG9RUREREQOl6arioiIiIjEoKJZRERERCQGFc0iIiIiIjGoaBYRERERiaHZNzeJh5ltBbS7ScvQEShNdhDSpJTj9KA8pwflOT0oz4enu7sfsElIiyiapeUws7kN7aIjqUM5Tg/Kc3pQntOD8pwYas8QEREREYlBRbOIiIiISAwqmiXR/pTsAKTJKcfpQXlOD8pzelCeE0A9zSIiIiIiMWikWUREREQkBhXNIiIiIiIxqGgWEREREYlBRbOIiIiISAwqmiVuZjbSzH5gZj32O35DciKS5mBm/0l2DJJYZnaZmbWPPi42s2fMbKGZvWBmXZIdnySGmbU3s7vN7EaLuMvM/mlmvzGzdsmOTxIj+rf5D2Y21cxeiT6+INlxpQKtniFxMbP7gTOBecA3gYfc/ffRc/Pc/SvJjE8Sw8wW7H8I6AssB3D3U5s9KEk4M1vi7idGH78AfABMBs4Fvufu5yUzPkkMM5sGLATaAv2jj18EzgMGuPulSQxPEsDMHiLyO/oZYF30cBfgWuBzd78lWbGlAhXNEhczWwgMcvegmRUBzwHL3f1WM/vE3QclOURJADObCuwCJgLVRIrmmUQ+MOHua5MXnSSKmS139xOijz9299P2Ojff3QcmLzpJlN25NDMD1rl7yf7nkhieJICZfebufRs4bsBn7t4nCWGlDLVnSLyy3D0I4O47iIw2tzWzyUBOUiOThHH3S4ApRBbGH+Dua4CAu69VwZxS3jaze80sP/r4WwBmdg6wM7mhSQJlRNswugKtd7fWmVkH9Hs7VdSY2ZAGjp8O1DR3MKkmK9kBSIu10sxGuPs7AO4eAn5gZhOBy5MbmiSSu79sZq8D/2tmN6I/rqloPHAX0bYb4FYzqwReBUYnLSpJtF8Dy6KPbwAmmZkDJwL3JC0qSaTrgD+YWRv2tGd0JXLH8LokxZQy1J4hcYmOSOHu1Q2cK3H39c0flTQ1MxsADHP3J5IdizQNMyskcidpW7JjkcQzs0wif/uDZpYFDATWu/vGJIcmCWRmnYESIi1169x9U5JDSgkqmiVu0R6pIUTemA5sAD5y/U+VUpTn9KA8pwflOX2ZWT93Xxb7SmmMimaJi5mdDzwOfA7sHlXuAvQGbnb315MVmySO8pwelOf0oDynNzP7wt27JTuOlkxFs8TFzJYCF0Ynhu19vCcwzd37JyUwSSjlOT0oz+lBeU59ZvZIY6eA77t72+aMJ9VoIqDEK4s9kwz2th7IbuZYpOkoz+lBeU4PynPqux64Daht4Nx3mzmWlKOiWeL1FDDHzP4GfBk91hW4CngyaVFJoinP6UF5Tg/Kc+qbAyxy91n7nzCzCc0fTmpRe4bEzcxOBC5hrxm6wFR3X5LUwCShlOf0oDynB+U5tZlZe6DG3auSHUsqUtEsIiIiIhKDdgSUuJhZoZk9YGbLzGxb9L+l0WNFyY5PEkN5Tg/Kc3pQnlOfcty0VDRLvF4EyoCz3b2Du3cAzgF2AJOTGpkkkvKcHpTn9KA8p77GclyGcnzE1J4hcTGz5e5+wuGek5ZFeU4PynN6UJ5Tn3LctDTSLPFaa2Z3mtkxuw+Y2TFm9t/smZUtLZ/ynB6U5/SgPKc+5bgJqWiWeF0JdADeMbMyM9sOvA20B76TzMAkoZTn9KA8pwflOfUpx01I7RkSNzPrR2QL1g/cvWKv4xe4+7+TF5kkkvKcHpTn9KA8pz7luOlopFniYmY/Bl4BxgOLzOzSvU7fn5yoJNGU5/SgPKcH5Tn1KcdNSzsCSrzGAKe5e4WZ9QBeMrMe7v4wkQXzJTUoz+lBeU4PynPqU46bkIpmiVfm7ts+7r7GzM4m8ubsjt6YqUR5Tg/Kc3pQnlOfctyE1J4h8dpkZgN3P4m+SS8GOgKnJC0qSTTlOT0oz+lBeU59ynET0kRAiYuZdQGC7r6pgXPD3f39JIQlCaY8pwflOT0oz6lPOW5aKppFRERERGJQe4aIiIiISAwqmkVEREREYlDRLCLSgphZyMzmm9liM/vUzH5qZgf9XW5mPczs6uaKUUQkFaloFhFpWardfaC7nwScB1wE/E+M1/QAVDSLiBwBTQQUEWlBzKzC3Vvv9bwXMIfIklLdgWeBgujp8e4+y8w+APoDq4G/AI8ADwBnA7nAY+7+x2b7IUREWiAVzSIiLcj+RXP0WBnQDygHwu5eY2Z9gOfdfXB0g4Pb3f3i6PVjgU7uPtHMcoH3gSvcfXWz/jAiIi2IdgQUEWn5du/0lQ08Gt3cIAT0beT684FTzWxU9Hkh0IfISLSIiDRARbOISAsWbc8IAVuI9DZvBgYQmbNS09jLgB+5+/RmCVJEJAVoIqCISAtlZsXAE8CjHum1KwQ2unsYGA1kRi8tB9rs9dLpwA/NLDv6dfqaWQEiItIojTSLiLQs+WY2n0grRpDIxL8Ho+ceB6aY2RXADKAyenwBEDSzT4E/Aw8TWVFjnpkZsBX4VnP9ACIiLZEmAoqIiIiIxKD2DBERERGRGFQ0i4iIiIjEoKJZRERERCQGFc0iIiIiIjGoaBYRERERiUFFs4iIiIhIDCqaRURERERiUNEsIiIiIhLD/wOcCQm0uZ98uwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Evaluate on test set\n",
    "obs, preds = eval_model(model, test_loader)\n",
    "preds = ds_val.local_rescale(preds.numpy(), variable='output')\n",
    "obs = obs.numpy()\n",
    "nse = calc_nse(obs, preds)\n",
    "\n",
    "# Plot results\n",
    "start_date = ds_test.dates[0]\n",
    "end_date = ds_test.dates[1] + pd.DateOffset(days=1)\n",
    "date_range = pd.date_range(start_date, end_date)\n",
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "ax.plot(date_range, obs, label=\"observation\")\n",
    "ax.plot(date_range, preds, label=\"prediction\")\n",
    "ax.legend()\n",
    "ax.set_title(f\"Basin {basin} - Test set NSE: {nse:.3f}\")\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "ax.set_xlabel(\"Date\")\n",
    "_ = ax.set_ylabel(\"Discharge (mm/d)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congratulations, you have trained a LSTM. Although results for single basin calibration might not look superior to traditional hydrological models don't forget that this model had to learn every hydrological process from scratch. Also note that, due to computational contraints of the Pangeo servers, we trained very small LSTM networks. However, while getting the overall dynamics farely well, we can see that the LSTM (as most hydrological models) has problems with peak flows and is mostly underestimating them. \n",
    "\n",
    "I hope you are able with the code provided in this notebook to try LSTMs on your own data/problems. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}