diff --git a/PatchTST_self_supervised/src/callback/tracking.py b/PatchTST_self_supervised/src/callback/tracking.py index 0f25a96b..27e7b7f3 100644 --- a/PatchTST_self_supervised/src/callback/tracking.py +++ b/PatchTST_self_supervised/src/callback/tracking.py @@ -7,6 +7,7 @@ import time import numpy as np from pathlib import Path +import logging class TrackTimerCB(Callback): @@ -121,7 +122,10 @@ def after_epoch_valid(self): self.recorder['valid_'+name].append( values[name] ) - def after_batch_train(self): self.accumulate() # save batch recorder + def after_batch_train(self): + self.accumulate() # save batch recorder + logging.info(f"Batch loss: {self.batch_recorder['batch_losses'][-1]}") + def after_batch_valid(self): self.accumulate() def accumulate(self ): @@ -129,7 +133,7 @@ def accumulate(self ): bs = len(xb) self.batch_recorder['n_samples'].append(bs) # get batch loss - loss = self.loss.detach()*bs if self.mean_reduction_ else self.loss.detach() + loss = self.loss.detach()*bs if self.mean_reduction_ else self.loss.detach() self.batch_recorder['batch_losses'].append(loss) if yb is None: self.batch_recorder['with_metrics'] = False diff --git a/PatchTST_self_supervised/src/learner.py b/PatchTST_self_supervised/src/learner.py index 276cea1e..74a4fbe3 100644 --- a/PatchTST_self_supervised/src/learner.py +++ b/PatchTST_self_supervised/src/learner.py @@ -122,7 +122,7 @@ def fit(self, n_epochs, lr=None, cbs=None, do_valid=True): def fit_one_cycle(self, n_epochs, lr_max=None, pct_start=0.3): self.n_epochs = n_epochs self.lr_max = lr_max if lr_max else self.lr - cb = OneCycleLR(lr_max=self.lr_max, pct_start=pct_start) + cb = OneCycleLR(lr_max=self.lr_max, pct_start=pct_start, verbose=True) self.fit(self.n_epochs, cbs=cb) def one_epoch(self, train): diff --git a/SeasonTST/SeasonTST_finetune.py b/SeasonTST/SeasonTST_finetune.py index 233db0f7..f6a5e9ac 100644 --- a/SeasonTST/SeasonTST_finetune.py +++ b/SeasonTST/SeasonTST_finetune.py @@ -39,7 +39,7 @@ datefmt="%m/%d/%Y %I:%M:%S %p", filename=f'logs/{datetime.datetime.now().strftime("%Y_%m_%d_%I:%M")}_finetune.log', encoding="utf-8", - level=logging.DEBUG, + level=logging.INFO, ) @@ -51,10 +51,6 @@ def finetune_func(learner, save_path, args, lr=0.001): print("end-to-end finetuning") - if not os.path.exists(save_path): - os.makedirs(save_path) - - print(save_path) # fit the data to the model and save learner.fine_tune( n_epochs=args.n_epochs_finetune, base_lr=lr, freeze_epochs=args.freeze_epochs @@ -107,20 +103,6 @@ def save_recorders(learner, args): ) -def test_func(weight_path, learner, args, dls): - - out = learner.test( - dls.test, weight_path=weight_path, scores=[mse, mae] - ) # out: a list of [pred, targ, score] - print("score:", out[2]) - # save results - pd.DataFrame(np.array(out[2]).reshape(1, -1), columns=["mse", "mae"]).to_csv( - args.save_path + args.save_finetuned_model + "_acc.csv", - float_format="%.6f", - index=False, - ) - return out - def load_config(): @@ -135,13 +117,14 @@ def load_config(): "revin": 0, # reversible instance normalization "mask_ratio": 0.4, # masking ratio for the input "lr": 1e-3, - "batch_size": 128, + "batch_size": 64, + "drop_last": False, "num_workers": 6, "prefetch_factor": 3, - "n_epochs_pretrain": 1, # number of pre-training epochs, + "n_epochs_pretrain": 20, # number of pre-training epochs, "freeze_epochs": 0, - "n_epochs_finetune": 250, - "pretrained_model_id": 2500, # id of the saved pretrained model + "n_epochs_finetune": 10, + "pretrained_model_id": 2, # id of the saved pretrained model "save_finetuned_model": "./finetuned_d128", "save_path": "saved_models" + "/masked_patchtst/", } @@ -186,17 +169,18 @@ def main(): # Create dataloader dls = get_dls(config_obj, SeasonTST_Dataset, data, mask) - # suggested_lr = find_lr(config_obj, dls) # This is what I got on a small dataset. In case one wants to skip this for testing. - suggested_lr = 0.00017073526474706903 + suggested_lr = 0.0002 # 0.000298364724028334 + learner = get_learner(config_obj, dls, suggested_lr, model) + suggested_lr = learner.lr_finder() print(suggested_lr) - learner = get_learner(config_obj, dls, suggested_lr, model) # This function will save the model weights to config_obj.save_finetuned_model. ie will not overwrite the pretrained model. - # However, there is currently no set-up to do finetuning from the result of a previous finetuning. + # To continue training from a previous fine-tuning checkpoint, the path needs to be explicity fed to the get_model function finetune_func(learner, pretrained_model_path, config_obj, suggested_lr) if __name__ == "__main__": + # PYTHONPATH=$(pwd) python SeasonTST/SeasonTST_finetune.py main() diff --git a/SeasonTST/SeasonTST_pretrain.py b/SeasonTST/SeasonTST_pretrain.py index a84c9190..6032c448 100644 --- a/SeasonTST/SeasonTST_pretrain.py +++ b/SeasonTST/SeasonTST_pretrain.py @@ -29,7 +29,7 @@ datefmt="%m/%d/%Y %I:%M:%S %p", filename=f'logs/{datetime.datetime.now().strftime("%Y_%m_%d_%I_%M")}_train.log', encoding="utf-8", - level=logging.DEBUG, + level=logging.INFO, ) @@ -95,10 +95,11 @@ def load_config(): "mask_value": -99, # Value to assign to masked elements of data input "lr": 1e-3, "batch_size": 128, + "drop_last":True, "prefetch_factor": 3, "num_workers": 6, - "n_epochs_pretrain": 1, # number of pre-training epochs - "pretrained_model_id": 2500, # id of the saved pretrained model + "n_epochs_pretrain": 20, # number of pre-training epochs + "pretrained_model_id": 2, # id of the saved pretrained model } config_obj = SimpleNamespace(**config) @@ -109,37 +110,42 @@ def main(): data, mask = load_data() config_obj = load_config() + save_path = "saved_models" + "/masked_patchtst/" + pretrained_model = ( + "patchtst_pretrained_cw" + + str(config_obj.sequence_length) + + "_patch" + + str(config_obj.patch_len) + + "_stride" + + str(config_obj.stride) + + "_epochs-pretrain" + + str(config_obj.n_epochs_pretrain) + + "_mask" + + str(config_obj.mask_ratio) + + "_model" + + str(config_obj.pretrained_model_id) + ) + pretrained_model_path = save_path + pretrained_model + ".pth" + # Creates train valid and test datasets for one epoch. Notice that they are in different locations! dls = get_dls(config_obj, SeasonTST_Dataset, data, mask) - model = get_model(config_obj) + + model = get_model( + config_obj, headtype="pretrain", weights_path=pretrained_model_path, exclude_head=False + ) # suggested_lr = find_lr(config_obj, dls) # This is what I got on a small dataset. In case one wants to skip this for testing. suggested_lr = 0.00020565123083486514 - save_pretrained_model = ( - "patchtst_pretrained_cw" - + str(config_obj.sequence_length) - + "_patch" - + str(config_obj.patch_len) - + "_stride" - + str(config_obj.stride) - + "_epochs-pretrain" - + str(config_obj.n_epochs_pretrain) - + "_mask" - + str(config_obj.mask_ratio) - + "_model" - + str(config_obj.pretrained_model_id) - ) - save_path = "saved_models" + "/masked_patchtst/" + + pretrain_func( - save_pretrained_model, save_path, config_obj, model, dls, suggested_lr + pretrained_model, save_path, config_obj, model, dls, suggested_lr ) - pretrained_model_name = save_path + save_pretrained_model + ".pth" - - model = transfer_weights(pretrained_model_name, model) + model = transfer_weights(pretrained_model_path, model) if __name__ == "__main__": diff --git a/SeasonTST/dataset.py b/SeasonTST/dataset.py index ac41d68c..1776b9d0 100644 --- a/SeasonTST/dataset.py +++ b/SeasonTST/dataset.py @@ -125,7 +125,7 @@ def scale(self, batch): for var, data_var in batch.data_vars.items(): batch[var] = ( data_var - self.scaling_factors["mean"][var] - ) / self.scaling_factors["mean"][var] + ) / self.scaling_factors["std"][var] return batch def __len__(self): diff --git a/SeasonTST/utils.py b/SeasonTST/utils.py index 66aa7fae..57343c7b 100644 --- a/SeasonTST/utils.py +++ b/SeasonTST/utils.py @@ -30,6 +30,7 @@ def get_dls( batch_size=config_obj.batch_size, workers=config_obj.num_workers, prefetch_factor=config_obj.prefetch_factor, + drop_last=config_obj.drop_last ) dls.vars, dls.len = dls.train.dataset[0][0].shape[1], config_obj.sequence_length @@ -74,17 +75,14 @@ def get_model(config, headtype="pretrain", weights_path=None, exclude_head=True) return model -def find_lr(config_obj, dls): +def find_lr(model, config_obj, dls): """ # This method typically involves training the model for a few epochs with a range of learning rates and recording the loss at each step. The learning rate that gives the fastest decrease in loss is considered optimal or near-optimal for the training process. - :param config_obj: - :return: """ - model = get_model(config_obj) # get loss loss_func = torch.nn.MSELoss(reduction="mean") # get callbacks diff --git a/SeasonTST_evaluation.ipynb b/SeasonTST_evaluation.ipynb new file mode 100644 index 00000000..93d1d64d --- /dev/null +++ b/SeasonTST_evaluation.ipynb @@ -0,0 +1,655 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "fa27dbbf-e087-4614-955f-a84410742504", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "f0eba995-b17e-4d62-a2e5-60641b2e574d", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "dd1bb57f-b345-4b2d-8957-393d609a7fbd", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/paolo/miniforge3/envs/PatchTST/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n", + " from .autonotebook import tqdm as notebook_tqdm\n" + ] + } + ], + "source": [ + "import os\n", + "from types import SimpleNamespace\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import torch\n", + "import xarray as xr\n", + "from dask.cache import Cache\n", + "\n", + "from PatchTST_self_supervised.src.callback.patch_mask import PatchCB, ObservationMaskCB\n", + "from PatchTST_self_supervised.src.callback.tracking import SaveModelCB\n", + "from PatchTST_self_supervised.src.callback.transforms import RevInCB\n", + "from PatchTST_self_supervised.src.learner import Learner, transfer_weights\n", + "from SeasonTST.dataset import SeasonTST_Dataset\n", + "from SeasonTST.utils import find_lr, get_dls, get_model, load_data\n", + "from PatchTST_self_supervised.src.metrics import mse, mae\n", + "\n", + "\n", + "#\n", + "# SETUP\n", + "#\n", + "\n", + "# Set up Dask's cache. Will reduce repeat reads from zarr and speed up data loading\n", + "cache = Cache(1e10) # 10gb cache\n", + "cache.register()\n", + "\n", + "import logging\n", + "import datetime\n", + "\n", + "logger = logging.getLogger(__name__)\n", + "logging.basicConfig(\n", + " format=\"%(asctime)s %(levelname)s %(module)s - %(funcName)s: %(message)s\",\n", + " datefmt=\"%m/%d/%Y %I:%M:%S %p\",\n", + " filename=f'logs/{datetime.datetime.now().strftime(\"%Y_%m_%d_%I_%M\")}_evaluation.log',\n", + " encoding=\"utf-8\",\n", + " level=logging.DEBUG,\n", + ")\n", + "\n", + "\n", + "#\n", + "# FUNCTIONS\n", + "#\n", + "\n", + "\n", + "\n", + "def get_learner(args, dls, lr, model):\n", + " \"\"\"\n", + " Learner set-up\n", + "\n", + " TRAINING\n", + " - Input is [bs, seq_len, n_vars]\n", + " - Before forward pass:\n", + " - RevInCB normalized inputs\n", + " - ObservationMaskCB masks random observations with fill value\n", + " - PatchCB reshaped to [bs, num_patches, n_vars, patch_len]\n", + " - Forward pass in: [bs, num_patches, n_vars, patch_len]; out: [bs, pred_len, n_vars]\n", + " - After forward pass\n", + " - RevInCB denormalized outputs\n", + " - ObservationMaskCB custom loss function on outputs for just the masked values\n", + " - Loss is therefore mean squared difference on denormalized masked values.\n", + " - Will give more weight to variables with larger numerical range\n", + " \"\"\"\n", + "\n", + " # get loss\n", + " loss_func = torch.nn.MSELoss(reduction=\"mean\")\n", + " # get callbacks\n", + " cbs = [RevInCB(dls.vars, denorm=True)] if args.revin else []\n", + " cbs += [\n", + " # ObservationMaskCB(mask_ratio=0.2, mask_value=-99),\n", + " PatchCB(patch_len=args.patch_len, stride=args.stride),\n", + " SaveModelCB(\n", + " monitor=\"valid_loss\", fname=args.save_finetuned_model, path=args.save_path\n", + " ),\n", + " ]\n", + " # define learner\n", + " learner = Learner(dls, model, loss_func, lr=lr, cbs=cbs, metrics=[mse])\n", + " return learner\n", + "\n", + "\n", + "def save_recorders(learner, args):\n", + " train_loss = learner.recorder[\"train_loss\"]\n", + " valid_loss = learner.recorder[\"valid_loss\"]\n", + " df = pd.DataFrame(data={\"train_loss\": train_loss, \"valid_loss\": valid_loss})\n", + " df.to_csv(\n", + " args.save_path + args.save_finetuned_model + \"_losses.csv\",\n", + " float_format=\"%.6f\",\n", + " index=False,\n", + " )\n", + "\n", + "\n", + "def test_func(weight_path, learner, args, dls):\n", + "\n", + " out = learner.test(\n", + " dls.test, weight_path=weight_path, scores=[mse, mae]\n", + " ) # out: a list of [pred, targ, score]\n", + " print(\"score:\", out[2])\n", + " # save results\n", + " pd.DataFrame(np.array(out[2]).reshape(1, -1), columns=[\"mse\", \"mae\"]).to_csv(\n", + " args.save_path + args.save_finetuned_model + \"_acc.csv\",\n", + " float_format=\"%.6f\",\n", + " index=False,\n", + " )\n", + " return out\n", + "\n", + "\n", + "def load_config():\n", + "\n", + " # Config parameters\n", + " # TODO maybe load from a JSON with a model key?\n", + " config = {\n", + " \"c_in\": 5, # number of variables\n", + " \"sequence_length\": 36,\n", + " \"prediction_length\": 2, # Sets both the dimension of y from the dataloader as well as the prediction head size\n", + " \"patch_len\": 4, # Length of the patch\n", + " \"stride\": 4, # Minimum non-overlap between patchs. If equal to patch_len , patches will not overlap\n", + " \"revin\": 0, # reversible instance normalization\n", + " \"mask_ratio\": 0.4, # masking ratio for the input\n", + " \"lr\": 1e-3,\n", + " \"batch_size\": 128,\n", + " \"drop_last\": False, # Whether to drop the last observation that don't make a full batch\n", + " \"num_workers\": 0,\n", + " \"prefetch_factor\": 2,\n", + " \"n_epochs_pretrain\": 1, # number of pre-training epochs,\n", + " \"freeze_epochs\": 0,\n", + " \"n_epochs_finetune\": 250,\n", + " \"pretrained_model_id\": 2, # id of the saved pretrained model\n", + " \"save_finetuned_model\": \"./finetuned_d128\",\n", + " \"save_path\": \"saved_models\" + \"/masked_patchtst/\",\n", + " }\n", + " config_obj = SimpleNamespace(**config)\n", + "\n", + " save_pretrained_model = (\n", + " \"patchtst_pretrained_cw\"\n", + " + str(config_obj.sequence_length)\n", + " + \"_patch\"\n", + " + str(config_obj.patch_len)\n", + " + \"_stride\"\n", + " + str(config_obj.stride)\n", + " + \"_epochs-pretrain\"\n", + " + str(config_obj.n_epochs_pretrain)\n", + " + \"_mask\"\n", + " + str(config_obj.mask_ratio)\n", + " + \"_model\"\n", + " + str(config_obj.pretrained_model_id)\n", + " )\n", + " save_path = \"saved_models\" + \"/masked_patchtst/\"\n", + " pretrained_model_path = save_path + save_pretrained_model + \".pth\"\n", + "\n", + " return config_obj, save_path, pretrained_model_path" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "aed4fece-e91e-4a79-84d0-e7f08e970f58", + "metadata": {}, + "outputs": [], + "source": [ + "#\n", + "# EVALUATION STEPS\n", + "#" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "36401829-11bf-4da5-9668-ce278d1a363a", + "metadata": {}, + "outputs": [], + "source": [ + "data, mask = load_data()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "1340dfdf-6630-47db-b5f8-6435268ed0f8", + "metadata": {}, + "outputs": [], + "source": [ + "data = data.compute()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "dbc6b92e-0310-4306-87db-e48d83cd57c4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAHFCAYAAAAdTZjVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABC30lEQVR4nO3deXgUVdr38V+TkIQACUIgAQ0kKBJQRxAEg8MkPiOLOy4oMqK4MCLDQGAYlEElgICgD6IDiAsQRlF5lEXeGWRRgREJmxAXQHA0ENSETUhAJCHd5/0D00OThaS6sxT1/VxXXZd9+pyqu48Lt+c+VeUyxhgBAAA4RK3qDgAAAKAqkfwAAABHIfkBAACOQvIDAAAcheQHAAA4CskPAABwFJIfAADgKCQ/AADAUUh+AACAo5D8wPHWr1+v1NRUHT16tNh3ycnJSk5OrvKYqsLu3bs1YsQIdejQQQ0aNFDDhg117bXX6r333iux/4EDB9S/f39FRUUpPDxciYmJ+uijj3z65OXlacKECUpOTlZMTIzq1aunK664QpMnT9bJkyfLjOfDDz+Uy+WSy+XSoUOHyv07KjuuIhWZr++//14pKSlKSkpSgwYN5HK5lJaWVu7fBKBykfzA8davX6+xY8eWmPzMnDlTM2fOrPqgqsDKlSv1r3/9S3feeafeffddzZ8/X61atVLv3r01btw4n775+fn6/e9/r48++kgvvvii3n//fUVHR6tnz55au3att19WVpamTZumq666Sq+++qqWLl2qu+66S6mpqbr55ptV2tt0jh8/rgEDBqhZs2YV+g2VHZfV+frPf/6j+fPnKyQkRDfeeGOFfhOAKmAAh3vuueeMJJOZmVndoVSpgwcPGo/HU6z9pptuMuHh4ebkyZPethkzZhhJZv369d62U6dOmbZt25pOnTp5244fP26OHz9e7JxFc/zJJ5+UGMuf/vQn0759e/Pkk08aSebgwYPl+g2VHdeZKjJfbrfb+9ebN282kszcuXPL9ZsAVD5WfuBoqamp+utf/ypJio+P95Zd1qxZI6l42WvPnj1yuVx67rnnNHnyZMXFxalOnTpKTk7W7t27derUKT3xxBNq1qyZIiMjdfvtt+vAgQPFrrtgwQIlJiaqbt26qlevnnr06KFt27ZVxU/2ioqKksvlKtbeqVMnnThxQj/99JO3bfHixWrdurUSExO9bcHBwbrvvvu0adMm/fDDD5KkunXrqm7duiWeU5L27dtX7LtPPvlEr776ql5//XUFBQVV6DdUZlxnq8h81arFf1qBmox/Q+FojzzyiP785z9LkhYtWqT09HSlp6frqquuKnPcjBkz9Omnn2rGjBl6/fXX9fXXX+uWW27Rww8/rIMHD2rOnDmaMmWKPvzwQz3yyCM+YydOnKh7771Xbdu21f/93//pjTfe0LFjx9S1a1ft2LHjnDEXFhaW6zDlKOWUZPXq1WrcuLGaNGnibfvqq6/0m9/8pljforbt27eXec6PP/5YknTZZZf5tP/yyy96+OGHlZKScs45L0llxVURJc0XgJotuLoDAKrTRRddpObNm0uS2rdvr7i4uHKNa9CggZYsWeL9P/xDhw4pJSVFCQkJev/99739vv76a02bNk15eXmKiIjQvn37NGbMGA0ePFgvvfSSt1+3bt3UqlUrjR07VgsWLCj1unv27FF8fHy5Yly9enWFN2u//vrrWrNmjV588UWfVZjDhw+rYcOGxfoXtR0+fLjUc37xxReaMmWKbr/99mKJylNPPSW3262xY8dWKM7Kjqu8SpsvADUbyQ9gwY033uhT2mjTpo0k6aabbvLpV9SelZWlyy+/XCtWrFBhYaHuv/9+FRYWevuFhYUpKSlJq1evLvO6zZo10+bNm8sVY+vWrcvVr8gHH3ygP/3pT7rrrru8q2FnKqnkc67v9uzZo5tvvlmxsbF6/fXXfb7btGmTpk2bpuXLl6tOnTqlntsYI7fb7dMWHPzf/3QFOq4z/75IUlBQUInnOdd8Aai5SH4AC85ebQgJCSmzveh26v3790uSrr766hLPe669IiEhIWrXrl25YqzISsSKFSt0xx13qFu3bpo/f36xP+wbNWpU4ipK0T6XklZf9u7dq+uuu07BwcH66KOPivV56KGHdMcdd6hjx47eO+2K5ikvL0+hoaGqX7++5s2bpwcffNBnbFFJrzLiql27ts/nuXPnqn///j5t55ovADUbyQ9QhaKioiRJ7733nlq0aFHh8ZVR9lqxYoV69eqlpKQkLVy40JuwnemKK67Ql19+Way9qO3yyy/3ad+7d6+Sk5NljNGaNWt00UUXFRu7fft2bd++Xe+++26x7y6++GJdeeWVysjI0C233FLqaldlxHX2tc6e7/LMF4CajeQHjhcaGirp9ObbytajRw8FBwfr22+/1Z133lnh8YEue61cuVK9evXSb3/7Wy1ZssQ7F2e7/fbbNWjQIG3cuFGdO3eWdLo89Oabb6pz584+z+fJyspScnKy3G631qxZU2qSV1KJLy0tTfPmzdOSJUt04YUXSjq9utOoUaMqi6tjx44ltkvlny8ANRvJDxzviiuukCS9+OKLeuCBB1S7dm21bt1a9evXD/i14uLiNG7cOI0ePVrfffedevbsqQsuuED79+/Xpk2bVLdu3TI3/4aEhJT5h3NFrFu3Tr169VJMTIz+9re/KSMjw+f7tm3bKiIiQtLpEtWMGTPUu3dvPfvss2rSpIlmzpypXbt26cMPP/SOOXDggK677jplZ2dr9uzZOnDggM+t/hdddJF3taWkVamiRwxce+213lWyslRGXIGYL0neJz9/9913kqQtW7aoXr16kqS77rrrnL8NQCWq1qcMATXEqFGjTLNmzUytWrWMJLN69WpjjDFJSUkmKSnJ2y8zM9NIMs8995zP+NWrVxtJ5t133/Vpnzt3rpFkNm/e7NO+ZMkSc91115mIiAgTGhpqWrRoYe666y7z4YcfVsrvK8mYMWOMpFKPojkokpOTY+6//37TsGFDExYWZq655hqzatUqnz5F81DaMWbMmHLFVN6HHFZVXGfGVt75KqsvgOrlMsbiw0AAAABsiIccAgAARyH5AQAAjkLyAwAAHMU2yc+RI0fUr18/RUZGKjIyUv369fM+GK00+/fvV//+/dWsWTOFh4erZ8+e+uabb6omYAAAUCPZJvnp27evMjIytHz5ci1fvlwZGRnq169fqf2NMerVq5e+++47vf/++9q2bZtatGih66+/Xj///HMVRg4AAGoSW9zttXPnTrVt21YbNmzwPshsw4YNSkxM1Ndff13iw9x2796t1q1b66uvvvK+sdntdqtJkyaaPHlysTdtAwAAZ7DFQw7T09MVGRnpTXwk6ZprrlFkZKTWr19fYvKTn58v6fQLI4sEBQUpJCRE69atKzX5yc/P946VJI/Ho59++kmNGjXi/T0AgDIZY3Ts2DE1a9bsnO/q88fJkydVUFDg93lCQkJ8/px0ClskPzk5OWrSpEmx9iZNmignJ6fEMQkJCWrRooVGjRqlV155RXXr1tXUqVOVk5Oj7OzsUq81adKkMp+wCwDAuezbt++cTw236uTJk4pvUU85B9x+nysmJkaZmZmOS4CqNflJTU09Z6JR9B6jklZdjDGlrsbUrl1bCxcu1MMPP6yGDRsqKChI119/vW644YYyrzdq1CgNHz7c+zk3N1fNmzdX7NgnVcth/3AAACrGc/Kk9o15plJej1OkoKBAOQfc2vtZnCLqW19dyjvmUYsOe1RQUEDyU5UGDx6sPn36lNknLi5OX3zxhfbv31/su4MHDyo6OrrUsR06dFBGRoZyc3NVUFCgxo0bq3PnzmW+Gyk0NLTElxXWCgsj+QEAlEtVbJOoV9+levWtX8cj527lqNbkJyoqqlwvL0xMTFRubq42bdqkTp06SZI2btyo3NxcdenS5ZzjIyMjJUnffPONtmzZovHjx/sXOAAA1cxtPHL7ccuS23gCF4zN2OJW9zZt2qhnz54aMGCANmzYoA0bNmjAgAG6+eabfTY7JyQkaPHixd7P7777rtasWeO93b1bt27q1auXunfvXh0/AwCAgPHI+H04lS02PEvS/PnzNWTIEG/icuutt2r69Ok+fXbt2qXc3Fzv5+zsbA0fPlz79+9X06ZNdf/99+upp56q0rgBAEDNYpvkp2HDhnrzzTfL7HP2I4uGDBmiIUOGVGZYAABUC4888qdw5d9oe7NN8gMAAP7LbYzcfjyn2J+xdmeLPT8AAACBwsoPAAA25O+mZTY8AwAAW/HIyE3yYwllLwAA4Cis/AAAYEOUvawj+QEAwIa428s6yl4AAMBRWPkBAMCGPL8e/ox3KpIfAABsyO3n3V7+jLU7kh8AAGzIbeTnW90DF4vdsOcHAAA4Cis/AADYEHt+rCP5AQDAhjxyyS2XX+OdirIXAABwFFZ+AACwIY85ffgz3qlIfgAAsCG3n2Uvf8baHWUvAADgKKz8AABgQ6z8WEfyAwCADXmMSx7jx91efoy1O8peAADAUVj5AQDAhih7WUfyAwCADblVS24/CjjuAMZiNyQ/AADYkPFzz49hzw8AAIAzsPIDAIANsefHOlZ+AACwIbep5fdhxcyZMxUfH6+wsDB16NBBn3zySZn9165dqw4dOigsLEwtW7bUrFmzSu37zjvvyOVyqVevXpZiKy+SHwAAUC4LFixQSkqKRo8erW3btqlr16664YYblJWVVWL/zMxM3Xjjjeratau2bdumv/3tbxoyZIgWLlxYrO/evXs1YsQIde3atbJ/BskPAAB25JFLHtXy46h42Wvq1Kl6+OGH9cgjj6hNmzaaNm2aYmNj9fLLL5fYf9asWWrevLmmTZumNm3a6JFHHtFDDz2k559/3qef2+3WH/7wB40dO1YtW7a0NB8VQfIDAIANFe358eeQpLy8PJ8jPz+/xOsVFBTos88+U/fu3X3au3fvrvXr15c4Jj09vVj/Hj16aMuWLTp16pS3bdy4cWrcuLEefvhhf6ak3Eh+AABwsNjYWEVGRnqPSZMmldjv0KFDcrvdio6O9mmPjo5WTk5OiWNycnJK7F9YWKhDhw5Jkj799FPNnj1br732WgB+TflwtxcAADbkz6bl0+ONJGnfvn2KiIjwtoeGhpY5zuXyLZcZY4q1nat/UfuxY8d033336bXXXlNUVFSF4vcHyQ8AADZ0es+PHy82/XVsRESET/JTmqioKAUFBRVb5Tlw4ECx1Z0iMTExJfYPDg5Wo0aNtH37du3Zs0e33HLLf+PyeCRJwcHB2rVrly6++OIK/a7yoOwFAADOKSQkRB06dNCqVat82letWqUuXbqUOCYxMbFY/5UrV6pjx46qXbu2EhIS9OWXXyojI8N73HrrrbruuuuUkZGh2NjYSvktrPwAAGBDHj/f7eWRqfCY4cOHq1+/furYsaMSExP16quvKisrSwMHDpQkjRo1Sj/88IP+8Y9/SJIGDhyo6dOna/jw4RowYIDS09M1e/Zsvf3225KksLAwXX755T7XaNCggSQVaw8kkh8AAGwoUHt+KuKee+7R4cOHNW7cOGVnZ+vyyy/XsmXL1KJFC0lSdna2zzN/4uPjtWzZMg0bNkwzZsxQs2bN9NJLL+nOO++0HHcguIyx8OsdJC8vT5GRkWox+RnVCgur7nAAADWY5+RJ7X38SeXm5pZrH40VRX8uvZVxucLrB1k+z4ljbvVt91WlxlpTsecHAAA4CmUvAABsyG1cchs/Xmzqx1i7I/kBAMCG3H5ueHZb2PB8vqDsBQAAHIWVHwAAbMhjasnjx91eHgff70TyAwCADVH2so6yFwAAcBRWfgAAsCGP/LtjyxO4UGyH5AcAABvyqJY8fr3ewrnFH+f+cgAA4Eis/AAAYEP+v9vLuesfJD8AANiQRy555M+eH57wDAAAbISVH+uc+8sBAIAjsfIDAIAN+f+QQ+euf5D8AABgQx7jksef5/w4+K3uzk37AACAI7HyAwCADXn8LHs5+SGHJD8AANiQ/291d27y49xfDgAAHImVHwAAbMgtl9x+PKjQn7F2R/IDAIANUfayzrm/HAAAOBIrPwAA2JBb/pWu3IELxXZIfgAAsCHKXtaR/AAAYEO82NQ62/zyCRMmqEuXLgoPD1eDBg3KNcYYo9TUVDVr1kx16tRRcnKytm/fXrmBAgCAGs02yU9BQYF69+6txx57rNxjpkyZoqlTp2r69OnavHmzYmJi1K1bNx07dqwSIwUAoPIZueTx4zDc6l7zjR07VpKUlpZWrv7GGE2bNk2jR4/WHXfcIUmaN2+eoqOj9dZbb+nRRx+trFABAKh0lL2sO29/eWZmpnJyctS9e3dvW2hoqJKSkrR+/fpSx+Xn5ysvL8/nAAAA54/zNvnJycmRJEVHR/u0R0dHe78ryaRJkxQZGek9YmNjKzVOAACs8BiX34dTVWvyk5qaKpfLVeaxZcsWv67hcvn+zTXGFGs706hRo5Sbm+s99u3b59f1AQCoDO5f3+ruz+FU1brnZ/DgwerTp0+ZfeLi4iydOyYmRtLpFaCmTZt62w8cOFBsNehMoaGhCg0NtXRNAABQ81Vr8hMVFaWoqKhKOXd8fLxiYmK0atUqtW/fXtLpO8bWrl2ryZMnV8o1AQCoKv6Wrih72UBWVpYyMjKUlZUlt9utjIwMZWRk6Pjx494+CQkJWrx4saTT5a6UlBRNnDhRixcv1ldffaX+/fsrPDxcffv2ra6fAQBAQHhUy+/DqWxzq/vTTz+tefPmeT8XreasXr1aycnJkqRdu3YpNzfX22fkyJH65ZdfNGjQIB05ckSdO3fWypUrVb9+/SqNHQAA1By2SX7S0tLO+YwfY4zPZ5fLpdTUVKWmplZeYAAAVAO3ccntR+nKn7F2Z5vkBwAA/Bd7fqwj+QEAwIaMn291NzzhGQAAwBlY+QEAwIbccsntx8tJ/RlrdyQ/AADYkMf4t2/HY87d53xF2QsAADgKKz8AANiQx88Nz/6MtTuSHwAAbMgjlzx+7NvxZ6zdOTftAwAAjsTKDwAANsQTnq0j+QEAwIbY82Odc385AABwJFZ+AACwIY/8fLeXgzc8k/wAAGBDxs+7vQzJDwAAsBPe6m4de34AAICjsPIDAIANcbeXdSQ/AADYEGUv65yb9gEAAEdi5QcAABvi3V7WkfwAAGBDlL2so+wFAAAchZUfAABsiJUf60h+AACwIZIf6yh7AQAAR2HlBwAAG2LlxzpWfgAAsCGj/97ubuUwFq87c+ZMxcfHKywsTB06dNAnn3xSZv+1a9eqQ4cOCgsLU8uWLTVr1iyf71977TV17dpVF1xwgS644AJdf/312rRpk8XoyofkBwAAGypa+fHnqKgFCxYoJSVFo0eP1rZt29S1a1fdcMMNysrKKrF/ZmambrzxRnXt2lXbtm3T3/72Nw0ZMkQLFy709lmzZo3uvfderV69Wunp6WrevLm6d++uH374wfLcnIvLGGM1+XOEvLw8RUZGqsXkZ1QrLKy6wwEA1GCekye19/EnlZubq4iIiEq5RtGfS//zr4EKrhtq+TyFP+fr45tmVSjWzp0766qrrtLLL7/sbWvTpo169eqlSZMmFev/+OOPa+nSpdq5c6e3beDAgfr888+Vnp5e4jXcbrcuuOACTZ8+Xffff38Ff1X5sPIDAIANBWrlJy8vz+fIz88v8XoFBQX67LPP1L17d5/27t27a/369SWOSU9PL9a/R48e2rJli06dOlXimBMnTujUqVNq2LBhRaek3Eh+AACwoUAlP7GxsYqMjPQeJa3gSNKhQ4fkdrsVHR3t0x4dHa2cnJwSx+Tk5JTYv7CwUIcOHSpxzBNPPKELL7xQ119/fUWnpNy42wsAAAfbt2+fT9krNLTsUprL5btXyBhTrO1c/Utql6QpU6bo7bff1po1axRWiVtNSH4AALChQN3qHhERUa49P1FRUQoKCiq2ynPgwIFiqztFYmJiSuwfHBysRo0a+bQ///zzmjhxoj788EP95je/qchPqTDKXgAA2JAxLr+PiggJCVGHDh20atUqn/ZVq1apS5cuJY5JTEws1n/lypXq2LGjateu7W177rnnNH78eC1fvlwdO3asUFxWkPwAAIByGT58uF5//XXNmTNHO3fu1LBhw5SVlaWBAwdKkkaNGuVzh9bAgQO1d+9eDR8+XDt37tScOXM0e/ZsjRgxwttnypQpevLJJzVnzhzFxcUpJydHOTk5On78eKX9DspeAADYUNHDCv0ZX1H33HOPDh8+rHHjxik7O1uXX365li1bphYtWkiSsrOzfZ75Ex8fr2XLlmnYsGGaMWOGmjVrppdeekl33nmnt8/MmTNVUFCgu+66y+daY8aMUWpqqrUfdw4kPwAA2FB1vd5i0KBBGjRoUInfpaWlFWtLSkrS1q1bSz3fnj17LMXhD8peAADAUVj5AQDAhqxsWj57vFOR/AAAYEO81d06kh8AAGyIlR/r2PMDAAAchZUfAABsyPhZ9nLyyg/JDwAANmQk/fqaLMvjnYqyFwAAcBRWfgAAsCGPXHJV8ROezxckPwAA2BB3e1lH2QsAADgKKz8AANiQx7jk4iGHlpD8AABgQ8b4ebeXg2/3ouwFAAAchZUfAABsiA3P1pH8AABgQyQ/1pH8AABgQ2x4to49PwAAwFFY+QEAwIa428s6kh8AAGzodPLjz56fAAZjM5S9AACAo7DyAwCADXG3l3UkPwAA2JD59fBnvFNR9gIAAI7Cyg8AADZE2cs6kh8AAOyIupdltil7TZgwQV26dFF4eLgaNGhQrjGLFi1Sjx49FBUVJZfLpYyMjEqNEQCAKvPryo/VQw5e+bFN8lNQUKDevXvrscceK/eYn3/+Wddee62effbZSowMAADYiW3KXmPHjpUkpaWllXtMv379JEl79uyphIgAAKg+POHZOtskP1UlPz9f+fn53s95eXnVGA0AACVjw7N1til7VZVJkyYpMjLSe8TGxlZ3SAAAIICqNflJTU2Vy+Uq89iyZUuVxjRq1Cjl5uZ6j3379lXp9QEAKJeiTcv+HA5VrWWvwYMHq0+fPmX2iYuLq5pgfhUaGqrQ0NAqvSYAABXFnh/rqjX5iYqKUlRUVHWGAAAAHMY2e36ysrKUkZGhrKwsud1uZWRkKCMjQ8ePH/f2SUhI0OLFi72ff/rpJ2VkZGjHjh2SpF27dikjI0M5OTlVHj8AAAFlAnA4lG3u9nr66ac1b9487+f27dtLklavXq3k5GRJp5Ob3Nxcb5+lS5fqwQcf9H4uKrGNGTNGqamplR80AACVhLu9rLNN8pOWlnbOZ/yYswqY/fv3V//+/SsvKAAAYDu2SX4AAMBZHFy68gfJDwAANkTZyzqSHwAA7Ii3ultmm7u9AAAAAoGVHwAAbMn16+HPeGci+QEAwI4oe1lG2QsAADiK5eTnk08+0X333afExET98MMPkqQ33nhD69atC1hwAACgFDzh2TJLyc/ChQvVo0cP1alTR9u2bVN+fr4k6dixY5o4cWJAAwQAACXgre6WWUp+nnnmGc2aNUuvvfaaateu7W3v0qWLtm7dGrDgAAAAAs3Shuddu3bpd7/7XbH2iIgIHT161N+YAADAORhz+vBnvFNZWvlp2rSp/vOf/xRrX7dunVq2bOl3UAAA4BzY82OZpeTn0Ucf1dChQ7Vx40a5XC79+OOPmj9/vkaMGKFBgwYFOkYAAICAsVT2GjlypHJzc3Xdddfp5MmT+t3vfqfQ0FCNGDFCgwcPDnSMAADgbP5uWnbwhmfLDzmcMGGCRo8erR07dsjj8aht27aqV69eIGMDAAClcJnThz/jncqvJzyHh4erY8eOgYoFAACUF094tqzcyc8dd9xR7pMuWrTIUjAAAACVrdzJT2RkpPevjTFavHixIiMjvSs/n332mY4ePVqhJAkAAFjEnh/Lyp38zJ071/vXjz/+uO6++27NmjVLQUFBkiS3261BgwYpIiIi8FECAABflL0ss3Sr+5w5czRixAhv4iNJQUFBGj58uObMmROw4AAAAALNUvJTWFionTt3FmvfuXOnPB6P30EBAIBz4CGHllm62+vBBx/UQw89pP/85z+65pprJEkbNmzQs88+qwcffDCgAQIAgBJQ9rLMUvLz/PPPKyYmRi+88IKys7MlnX7lxciRI/WXv/wloAECAAAEkqXkp1atWho5cqRGjhypvLw8SWKjMwAAVYm7vSzz6yGHEkkPAADVgSc8W2cp+YmPj5fLVXrG+N1331kOCAAAoDJZutsrJSVFQ4cO9R6DBg1SYmKicnNz9cc//jHQMQIAgLNV091eM2fOVHx8vMLCwtShQwd98sknZfZfu3atOnTooLCwMLVs2VKzZs0q1mfhwoVq27atQkND1bZtWy1evNhacOVkaeVn6NChJbbPmDFDW7Zs8SsgAABQMy1YsEApKSmaOXOmrr32Wr3yyiu64YYbtGPHDjVv3rxY/8zMTN14440aMGCA3nzzTX366acaNGiQGjdurDvvvFOSlJ6ernvuuUfjx4/X7bffrsWLF+vuu+/WunXr1Llz50r5HS5jTMCqft99953atWvn3QR9PsjLy1NkZKRaTH5GtcLCqjscAEAN5jl5Unsff1K5ubmVtic2UH8uWYm1c+fOuuqqq/Tyyy9729q0aaNevXpp0qRJxfo//vjjWrp0qc+zAQcOHKjPP/9c6enpkqR77rlHeXl5+uCDD7x9evbsqQsuuEBvv/221Z9XJktlr9K89957atiwYSBPCQAAKlFeXp7PkZ+fX2K/goICffbZZ+revbtPe/fu3bV+/foSx6Snpxfr36NHD23ZskWnTp0qs09p5wwES2Wv9u3b+2x4NsYoJydHBw8e1MyZMwMWHAAAKEWAbnWPjY31aR4zZoxSU1OLdT906JDcbreio6N92qOjo5WTk1PiJXJyckrsX1hYqEOHDqlp06al9intnIFgKfm57bbbfJKfWrVqqXHjxkpOTlZCQkLAggMAAKUI0BOe9+3b51P2Cg0NLXPY2Xd7G2PKvAO8pP5nt1f0nP6ylPyUlBECAAD7iYiIKNeen6ioKAUFBRVbkTlw4ECxlZsiMTExJfYPDg5Wo0aNyuxT2jkDwdKen6CgIB04cKBY++HDh33e9A4AACpJFd/qHhISog4dOmjVqlU+7atWrVKXLl1KHJOYmFis/8qVK9WxY0fVrl27zD6lnTMQLK38lHaDWH5+vkJCQvwKCAAAnFt1POF5+PDh6tevnzp27KjExES9+uqrysrK0sCBAyVJo0aN0g8//KB//OMfkk7f2TV9+nQNHz5cAwYMUHp6umbPnu1zF9fQoUP1u9/9TpMnT9Ztt92m999/Xx9++KHWrVtn/cedQ4WSn5deeknS6drc66+/rnr16nm/c7vd+ve//82eHwAAzlP33HOPDh8+rHHjxik7O1uXX365li1bphYtWkiSsrOzlZWV5e0fHx+vZcuWadiwYZoxY4aaNWuml156yfuMH0nq0qWL3nnnHT355JN66qmndPHFF2vBggWV9owfqYLP+YmPj5ck7d27VxdddJFPiSskJERxcXEaN25cpQZc1XjODwCgvKryOT9xz0zw+zk/e54cXamx1lQVWvnJzMyUJF133XVatGiRLrjggkoJCgAAnEOA7vZyIkt7flavXh3oOAAAAKpEuZOf4cOHa/z48apbt66GDx9eZt+pU6f6HRgAAChddWx4Pl+UO/nZtm2b91HUW7durdSHDwEAgHMI0BOenajcyc+Zpa41a9ZURiwAAKC82PNjmaWHHD700EM6duxYsfaff/5ZDz30kN9BAQAAVBZLyc+8efP0yy+/FGv/5ZdfvA82AgAAladoz48/h1NV6G6vvLw8GWNkjNGxY8cUdsbzBdxut5YtW6YmTZoEPEgAAHAWyl6WVSj5adCggVwul1wuly699NJi37tcLo0dOzZgwQEAAARahZKf1atXyxij//mf/9HChQvVsGFD73chISFq0aKFmjVrFvAgAQDAWfwtXbHyUz5JSUmSTj/pOTY2VrVqWdoyBAAA/EXZyzJLT3gueoHZiRMnlJWVpYKCAp/vf/Ob3/gfGQAAQCWwlPwcPHhQDz74oD744IMSv3e73X4FBQAAzoGVH8ss1a1SUlJ05MgRbdiwQXXq1NHy5cs1b948tWrVSkuXLg10jAAA4Czc6m6dpZWfjz/+WO+//76uvvpq1apVSy1atFC3bt0UERGhSZMm6aabbgp0nAAAAAFhaeXn559/9j7Pp2HDhjp48KAk6YorrtDWrVsDFx0AAECAWUp+WrdurV27dkmS2rVrp1deeUU//PCDZs2apaZNmwY0QAAAUAITgMOhLJW9UlJSlJ2dLUkaM2aMevTooTfffFMhISGaN29eQAMEAADF+btvhz0/FfSHP/zB+9ft27fXnj179PXXX6t58+aKiooKWHAAAACBVu7kZ/jw4eU+6dSpUy0FAwAAKsDBqzf+KHfys23btnL1c7lcloMBAADlxHN+LCt38rN69erKjAMAAKBKWNrzAwAAqhcbnq0j+QEAwI4oe1nGa9kBAICjsPIDAIANUfayzjYrPxMmTFCXLl0UHh6uBg0anLP/qVOn9Pjjj+uKK65Q3bp11axZM91///368ccfKz9YAAAqG094tsw2yU9BQYF69+6txx57rFz9T5w4oa1bt+qpp57S1q1btWjRIu3evVu33nprJUcKAABqMtuUvcaOHStJSktLK1f/yMhIrVq1yqft73//uzp16qSsrCw1b9480CECAFB12PBsmW2Sn0DIzc2Vy+Uqs2yWn5+v/Px87+e8vLwqiAwAgIphz491til7+evkyZN64okn1LdvX0VERJTab9KkSYqMjPQesbGxVRglAADlxJ4fy6o1+UlNTZXL5Srz2LJli9/XOXXqlPr06SOPx6OZM2eW2XfUqFHKzc31Hvv27fP7+gAAoOao1rLX4MGD1adPnzL7xMXF+XWNU6dO6e6771ZmZqY+/vjjMld9JCk0NFShoaF+XRMAgErHnh/LqjX5iYqKUlRUVKWdvyjx+eabb7R69Wo1atSo0q4FAEBVYs+PdbbZ85OVlaWMjAxlZWXJ7XYrIyNDGRkZOn78uLdPQkKCFi9eLEkqLCzUXXfdpS1btmj+/Plyu93KyclRTk6OCgoKqutnAACAamabu72efvppzZs3z/u5ffv2kk6/bT45OVmStGvXLuXm5kqSvv/+ey1dulSS1K5dO59znTkGAABbouxlmW2Sn7S0tHM+48eY//6djIuL8/kMAMD5hLKXdbYpewEAAASCbVZ+AADAGSh7WUbyAwCAHZH8WEbZCwAAOAorPwAA2JDr18Of8U5F8gMAgB1R9rKM5AcAABviVnfr2PMDAAAchZUfAADsiLKXZSQ/AADYlYMTGH9Q9gIAAI7Cyg8AADbEhmfrSH4AALAj9vxYRtkLAAA4Cis/AADYEGUv60h+AACwI8pellH2AgAAjsLKDwAANkTZyzqSHwAA7Iiyl2UkPwAA2BHJj2Xs+QEAAI7Cyg8AADbEnh/rSH4AALAjyl6WUfYCAACOwsoPAAA25DJGLmN9+cafsXZH8gMAgB1R9rKMshcAAAioI0eOqF+/foqMjFRkZKT69euno0ePljnGGKPU1FQ1a9ZMderUUXJysrZv3+79/qefftKf//xntW7dWuHh4WrevLmGDBmi3NzcCsdH8gMAgA0V3e3lz1FZ+vbtq4yMDC1fvlzLly9XRkaG+vXrV+aYKVOmaOrUqZo+fbo2b96smJgYdevWTceOHZMk/fjjj/rxxx/1/PPP68svv1RaWpqWL1+uhx9+uMLxUfYCAMCOamjZa+fOnVq+fLk2bNigzp07S5Jee+01JSYmateuXWrdunXxUIzRtGnTNHr0aN1xxx2SpHnz5ik6OlpvvfWWHn30UV1++eVauHChd8zFF1+sCRMm6L777lNhYaGCg8uf0rDyAwCAg+Xl5fkc+fn5fp0vPT1dkZGR3sRHkq655hpFRkZq/fr1JY7JzMxUTk6Ounfv7m0LDQ1VUlJSqWMkKTc3VxERERVKfCSSHwAAbClQZa/Y2Fjv3pzIyEhNmjTJr7hycnLUpEmTYu1NmjRRTk5OqWMkKTo62qc9Ojq61DGHDx/W+PHj9eijj1Y4RspeAADYUYDKXvv27VNERIS3OTQ0tMTuqampGjt2bJmn3Lx5syTJ5XIVv5wxJbaf6ezvSxuTl5enm266SW3bttWYMWPKPGdJSH4AALChQL3eIiIiwif5Kc3gwYPVp0+fMvvExcXpiy++0P79+4t9d/DgwWIrO0ViYmIknV4Batq0qbf9wIEDxcYcO3ZMPXv2VL169bR48WLVrl37nLGfjeQHAACcU1RUlKKios7ZLzExUbm5udq0aZM6deokSdq4caNyc3PVpUuXEsfEx8crJiZGq1atUvv27SVJBQUFWrt2rSZPnuztl5eXpx49eig0NFRLly5VWFiYpd/Cnh8AAOzIBOCoBG3atFHPnj01YMAAbdiwQRs2bNCAAQN08803+9zplZCQoMWLF0s6Xe5KSUnRxIkTtXjxYn311Vfq37+/wsPD1bdvX0mnV3y6d++un3/+WbNnz1ZeXp5ycnKUk5Mjt9tdoRhZ+QEAwKZq6pvZ58+fryFDhnjv3rr11ls1ffp0nz67du3yeUDhyJEj9csvv2jQoEE6cuSIOnfurJUrV6p+/fqSpM8++0wbN26UJF1yySU+58rMzFRcXFy54yP5AQAAAdWwYUO9+eabZfYxZ71bzOVyKTU1VampqSX2T05OLjbGKpIfAADsyJjThz/jHYrkBwAAGwrU3V5OxIZnAADgKKz8AABgRzX03V52QPIDAIANuTynD3/GOxVlLwAA4Cis/AAAYEeUvSwj+QEAwIa428s6kh8AAOyI5/xYxp4fAADgKKz8AABgQ5S9rCP5AQDAjtjwbBllLwAA4Cis/AAAYEOUvawj+QEAwI6428syyl4AAMBRWPkBAMCGKHtZR/IDAIAdcbeXZZS9AACAo7DyAwCADVH2so7kBwAAO/KY04c/4x2K5AcAADtiz49l7PkBAACOwsoPAAA25JKfe34CFon9kPwAAGBHPOHZMspeAADAUWyT/EyYMEFdunRReHi4GjRoUK4xqampSkhIUN26dXXBBRfo+uuv18aNGys3UAAAqkDRre7+HE5lm+SnoKBAvXv31mOPPVbuMZdeeqmmT5+uL7/8UuvWrVNcXJy6d++ugwcPVmKkAABUAROAw6Fss+dn7NixkqS0tLRyj+nbt6/P56lTp2r27Nn64osv9Pvf/z6Q4QEAAJuwTfLjr4KCAr366quKjIzUlVdeWWq//Px85efnez/n5eVVRXgAAFSIyxi5/Ni07M9Yu7NN2cuqf/7zn6pXr57CwsL0wgsvaNWqVYqKiiq1/6RJkxQZGek9YmNjqzBaAADKyROAw6GqNflJTU2Vy+Uq89iyZYtf17juuuuUkZGh9evXq2fPnrr77rt14MCBUvuPGjVKubm53mPfvn1+XR8AANQs1Vr2Gjx4sPr06VNmn7i4OL+uUbduXV1yySW65JJLdM0116hVq1aaPXu2Ro0aVWL/0NBQhYaG+nVNAAAqG2Uv66o1+YmKiiqzBFUZjDE+e3oAALAl3u1lmW32/GRlZSkjI0NZWVlyu93KyMhQRkaGjh8/7u2TkJCgxYsXS5J+/vln/e1vf9OGDRu0d+9ebd26VY888oi+//579e7du7p+BgAAgVH0hGd/Doeyzd1eTz/9tObNm+f93L59e0nS6tWrlZycLEnatWuXcnNzJUlBQUH6+uuvNW/ePB06dEiNGjXS1VdfrU8++USXXXZZlccPAABqBtskP2lpaed8xo85I4sNCwvTokWLKjkqAACqh79PaXbyE55tk/wAAIAz8GJTy2yz5wcAACAQWPkBAMCGXJ7Thz/jnYrkBwAAO6LsZRllLwAA4Cis/AAAYEc85NAykh8AAGyI11tYR9kLAAA4Cis/AADYERueLSP5AQDAjowkf25Xd27uQ/IDAIAdsefHOvb8AAAAR2HlBwAAOzLyc89PwCKxHZIfAADsiA3PllH2AgAAjsLKDwAAduSR5PJzvEOR/AAAYEPc7WUdZS8AAOAorPwAAGBHbHi2jOQHAAA7IvmxjLIXAABwFFZ+AACwI1Z+LCP5AQDAjrjV3TKSHwAAbIhb3a1jzw8AAHAUkh8AAOyoaM+PP0clOXLkiPr166fIyEhFRkaqX79+Onr06Dl+jlFqaqqaNWumOnXqKDk5Wdu3by+17w033CCXy6UlS5ZUOD6SHwAA7Mhj/D8qSd++fZWRkaHly5dr+fLlysjIUL9+/cocM2XKFE2dOlXTp0/X5s2bFRMTo27duunYsWPF+k6bNk0ul/UNT+z5AQAAAbNz504tX75cGzZsUOfOnSVJr732mhITE7Vr1y61bt262BhjjKZNm6bRo0frjjvukCTNmzdP0dHReuutt/Too496+37++eeaOnWqNm/erKZNm1qKkZUfAADsqIaWvdLT0xUZGelNfCTpmmuuUWRkpNavX1/imMzMTOXk5Kh79+7ettDQUCUlJfmMOXHihO69915Nnz5dMTExlmNk5QcAAFvyN4E5PTYvL8+nNTQ0VKGhoZbPmpOToyZNmhRrb9KkiXJyckodI0nR0dE+7dHR0dq7d6/387Bhw9SlSxfddtttluOTWPkBAMDRYmNjvRuTIyMjNWnSpBL7paamyuVylXls2bJFkkrcj2OMOec+nbO/P3PM0qVL9fHHH2vatGkWfqUvVn4AALCjAD3hed++fYqIiPA2l7bqM3jwYPXp06fMU8bFxemLL77Q/v37i3138ODBYis7RYpKWDk5OT77eA4cOOAd8/HHH+vbb79VgwYNfMbeeeed6tq1q9asWVNmbGci+QEAwI48RkWlK+vjpYiICJ/kpzRRUVGKioo6Z7/ExETl5uZq06ZN6tSpkyRp48aNys3NVZcuXUocEx8fr5iYGK1atUrt27eXJBUUFGjt2rWaPHmyJOmJJ57QI4884jPuiiuu0AsvvKBbbrnlnHGdieQHAAAETJs2bdSzZ08NGDBAr7zyiiTpj3/8o26++WafO70SEhI0adIk3X777XK5XEpJSdHEiRPVqlUrtWrVShMnTlR4eLj69u0r6fTqUEmbnJs3b674+PgKxUjyAwCAHRnP6cOf8ZVk/vz5GjJkiPfurVtvvVXTp0/36bNr1y7l5uZ6P48cOVK//PKLBg0apCNHjqhz585auXKl6tevH/D4SH4AALCjGvxW94YNG+rNN988x+V9r+9yuZSamqrU1NRyX+fsc5QXyQ8AAHYUoD0/TsSt7gAAwFFY+QEAwI5qcNmrpiP5AQDAjoz8TH4CFontUPYCAACOwsoPAAB2RNnLMpIfAADsyOOR5MezejyV95yfmo6yFwAAcBRWfgAAsCPKXpaR/AAAYEckP5ZR9gIAAI7Cyg8AAHbE6y0sI/kBAMCGjPHI+PFmdn/G2h3JDwAAdmSMf6s37PkBAABwBlZ+AACwI+Pnnh8Hr/yQ/AAAYEcej+TyY9+Og/f8UPYCAACOwsoPAAB2RNnLMpIfAABsyHg8Mn6UvZx8qztlLwAA4Cis/AAAYEeUvSwj+QEAwI48RnKR/FhB2QsAADgKKz8AANiRMZL8ec6Pc1d+SH4AALAh4zEyfpS9jIOTH9uUvSZMmKAuXbooPDxcDRo0qPD4Rx99VC6XS9OmTQt4bAAAVDnj8f9wKNskPwUFBerdu7cee+yxCo9dsmSJNm7cqGbNmlVCZAAAwE5sU/YaO3asJCktLa1C43744QcNHjxYK1as0E033VQJkQEAUPUoe1lnm+THCo/Ho379+umvf/2rLrvssuoOBwCAwDEe+bfh2bllr/M6+Zk8ebKCg4M1ZMiQco/Jz89Xfn6+93Nubq4kyXPyZMDjAwCcX4r+rKiKVZVCnfLrGYeFOhW4YGymWpOf1NRUbzmrNJs3b1bHjh0rfO7PPvtML774orZu3SqXy1XucZMmTSoxpn1jnqlwDAAAZzp8+LAiIyMr5dwhISGKiYnRupxlfp8rJiZGISEhAYjKXlymGot+hw4d0qFDh8rsExcXp7CwMO/ntLQ0paSk6OjRo2WOmzZtmoYPH65atf67p9vtdqtWrVqKjY3Vnj17Shx39srP0aNH1aJFC2VlZVXaP8jnq7y8PMXGxmrfvn2KiIio7nBsh/mzjrmzjrnzT25urpo3b64jR45YujO5vE6ePKmCggK/zxMSEuLzZ6xTVOvKT1RUlKKioirl3P369dP111/v09ajRw/169dPDz74YKnjQkNDFRoaWqw9MjKS/xBYFBERwdz5gfmzjrmzjrnzz5n/410ZwsLCHJm0BIpt9vxkZWXpp59+UlZWltxutzIyMiRJl1xyierVqydJSkhI0KRJk3T77berUaNGatSokc85ateurZiYGLVu3bqqwwcAADWEbZKfp59+WvPmzfN+bt++vSRp9erVSk5OliTt2rXLu0EZAACgJLZJftLS0s75jJ9zbV8qbZ9PWUJDQzVmzJgSS2EoG3PnH+bPOubOOubOP8yfPVTrhmcAAICqZpvXWwAAAAQCyQ8AAHAUkh8AAOAoJD8AAMBRHJ38FBYW6sknn1R8fLzq1Kmjli1baty4cfJ4yn7Z24wZM9SmTRvVqVNHrVu31j/+8Y8qirhmOXbsmFJSUtSiRQvVqVNHXbp00ebNm8scs3btWnXo0EFhYWFq2bKlZs2aVUXR1jwVnb/s7Gz17dtXrVu3Vq1atZSSklJ1wdYwFZ27RYsWqVu3bmrcuLEiIiKUmJioFStWVGHENUdF527dunW69tpr1ahRI9WpU0cJCQl64YUXqjDimsXKf/eKfPrppwoODla7du0qN0icm3GwZ555xjRq1Mj885//NJmZmebdd9819erVM9OmTSt1zMyZM039+vXNO++8Y7799lvz9ttvm3r16pmlS5dWYeQ1w913323atm1r1q5da7755hszZswYExERYb7//vsS+3/33XcmPDzcDB061OzYscO89tprpnbt2ua9996r4shrhorOX2ZmphkyZIiZN2+eadeunRk6dGjVBlyDVHTuhg4daiZPnmw2bdpkdu/ebUaNGmVq165ttm7dWsWRV7+Kzt3WrVvNW2+9Zb766iuTmZlp3njjDRMeHm5eeeWVKo68Zqjo/BU5evSoadmypenevbu58sorqyZYlMrRyc9NN91kHnroIZ+2O+64w9x3332ljklMTDQjRozwaRs6dKi59tprKyXGmurEiRMmKCjI/POf//Rpv/LKK83o0aNLHDNy5EiTkJDg0/boo4+aa665ptLirKmszN+ZkpKSHJv8+Dt3Rdq2bWvGjh0b6PBqtEDN3e23317mfyfPV/7M3z333GOefPJJM2bMGJKfGsDRZa/f/va3+uijj7R7925J0ueff65169bpxhtvLHVMfn5+sfep1KlTR5s2bdKpU6cqNd6apLCwUG63u8S5WLduXYlj0tPT1b17d5+2Hj16aMuWLY6aO8na/OG0QMydx+PRsWPH1LBhw8oIscYKxNxt27ZN69evV1JSUmWEWKNZnb+5c+fq22+/1ZgxYyo7RJRXdWdf1cnj8ZgnnnjCuFwuExwcbFwul5k4cWKZY0aNGmViYmLMli1bjMfjMZs3bzZNmjQxksyPP/5YRZHXDImJiSYpKcn88MMPprCw0LzxxhvG5XKZSy+9tMT+rVq1MhMmTPBp+/TTTx05d8ZUfP7O5OSVH2P8mztjjJkyZYpp2LCh2b9/fyVHWvNYnbsLL7zQhISEmFq1aplx48ZVUbQ1T0Xnb/fu3aZJkyZm165dxhjDyk8N4eiVnwULFujNN9/UW2+9pa1bt2revHl6/vnnfd4hdrannnpKN9xwg6655hrVrl1bt912m/r37y9JCgoKqqLIa4Y33nhDxhhdeOGFCg0N1UsvvaS+ffuWOQ8ul8vns/n1AeNntzuBlfnDaf7M3dtvv63U1FQtWLBATZo0qYJoaxarc/fJJ59oy5YtmjVrlqZNm6a33367iiKuWSoyf263W3379tXYsWN16aWXVkO0KFV1Zl7V7aKLLjLTp0/3aRs/frxp3br1OccWFBSYffv2mcLCQu8maLfbXVmh1mjHjx/3rtzcfffd5sYbbyyxX9euXc2QIUN82hYtWmSCg4NNQUFBpcdZU5V3/s7k9JWfIhWdu3feecfUqVOn2J4NJ7Lyz12R8ePHl3uV7XxVnvk7cuSIkWSCgoK8h8vl8rZ99NFHVR02fuXolZ8TJ06oVi3fKQgKCjrnre6SVLt2bV100UUKCgrSO++8o5tvvrnYuZyibt26atq0qY4cOaIVK1botttuK7FfYmKiVq1a5dO2cuVKdezYUbVr166KUGuk8s4fiqvI3L399tvq37+/3nrrLd10001VGGXN5M8/d8YY5efnV2J0NV955i8iIkJffvmlMjIyvMfAgQPVunVrZWRkqHPnztUQOSQ5e+XngQceMBdeeKH3VvdFixaZqKgoM3LkSG+fJ554wvTr18/7edeuXeaNN94wu3fvNhs3bjT33HOPadiwocnMzKyGX1C9li9fbj744APz3XffmZUrV5orr7zSdOrUybuKc/bcFd3qPmzYMLNjxw4ze/ZsR9/qXtH5M8aYbdu2mW3btpkOHTqYvn37mm3btpnt27dXR/jVqqJz99Zbb5ng4GAzY8YMk52d7T2OHj1aXT+h2lR07qZPn26WLl1qdu/ebXbv3m3mzJljIiIiKnR32PnEyr+3Z2LPT83g6OQnLy/PDB061DRv3tyEhYWZli1bmtGjR5v8/HxvnwceeMAkJSV5P+/YscO0a9fO1KlTx0RERJjbbrvNfP3119UQffVbsGCBadmypQkJCTExMTHmT3/6k88fJmfPnTHGrFmzxrRv396EhISYuLg48/LLL1dx1DWHlfmTVOxo0aJF1QZeA1R07pKSkkqcuwceeKDqg69mFZ27l156yVx22WUmPDzcREREmPbt25uZM2c6tsxv5d/bM5H81AwuY37dcQoAAOAAztykAgAAHIvkBwAAOArJDwAAcBSSHwAA4CgkPwAAwFFIfgAAgKOQ/AAAAEch+QFsLDk5WSkpKefNNfv3769evXpVyrkBoEhwdQcAwF4WLVrk8y62uLg4paSkVHkSBgBWkfwAqJCGDRtWdwgA4BfKXsB54siRI7r//vt1wQUXKDw8XDfccIO++eYb7/dpaWlq0KCBVqxYoTZt2qhevXrq2bOnsrOzvX0KCws1ZMgQNWjQQI0aNdLjjz+uBx54wKcUdWbZKzk5WXv37tWwYcPkcrnkcrkkSampqWrXrp1PfNOmTVNcXJz3s9vt1vDhw73XGjlypM5+244xRlOmTFHLli1Vp04dXXnllXrvvfcCM2EAHIvkBzhP9O/fX1u2bNHSpUuVnp4uY4xuvPFGnTp1ytvnxIkTev755/XGG2/o3//+t7KysjRixAjv95MnT9b8+fM1d+5cffrpp8rLy9OSJUtKveaiRYt00UUXady4ccrOzvZJpM7lf//3fzVnzhzNnj1b69at008//aTFixf79HnyySc1d+5cvfzyy9q+fbuGDRum++67T2vXri3/xADAWSh7AeeBb775RkuXLtWnn36qLl26SJLmz5+v2NhYLVmyRL1795YknTp1SrNmzdLFF18sSRo8eLDGjRvnPc/f//53jRo1Srfffrskafr06Vq2bFmp123YsKGCgoJUv359xcTEVCjmadOmadSoUbrzzjslSbNmzdKKFSu83//888+aOnWqPv74YyUmJkqSWrZsqXXr1umVV15RUlJSha4HAEVIfoDzwM6dOxUcHKzOnTt72xo1aqTWrVtr586d3rbw8HBv4iNJTZs21YEDByRJubm52r9/vzp16uT9PigoSB06dJDH4wlovLm5ucrOzvYmNZIUHBysjh07ektfO3bs0MmTJ9WtWzefsQUFBWrfvn1A4wHgLCQ/wHng7L0yZ7YX7cOR5HOXliS5XK5iY8/sX9a5y1KrVq1i484sv5VHUcL1r3/9SxdeeKHPd6GhoRWOCQCKsOcHOA+0bdtWhYWF2rhxo7ft8OHD2r17t9q0aVOuc0RGRio6OlqbNm3ytrndbm3btq3McSEhIXK73T5tjRs3Vk5Ojk8ClJGR4XOtpk2basOGDd62wsJCffbZZz6/KTQ0VFlZWbrkkkt8jtjY2HL9JgAoCSs/wHmgVatWuu222zRgwAC98sorql+/vp544gldeOGFuu2228p9nj//+c+aNGmSLrnkEiUkJOjvf/+7jhw5Umw16ExxcXH697//rT59+ig0NFRRUVFKTk7WwYMHNWXKFN11111avny5PvjgA0VERHjHDR06VM8++6xatWqlNm3aaOrUqTp69Kj3+/r162vEiBEaNmyYPB6Pfvvb3yovL0/r169XvXr19MADD1iaKwBg5Qc4T8ydO1cdOnTQzTffrMTERBljtGzZsmKlrrI8/vjjuvfee3X//fcrMTFR9erVU48ePRQWFlbqmHHjxmnPnj26+OKL1bhxY0lSmzZtNHPmTM2YMUNXXnmlNm3a5HNXmST95S9/0f3336/+/fsrMTFR9evX9260LjJ+/Hg9/fTTmjRpktq0aaMePXro//2//6f4+PgKzAwA+HIZKwV9AI7g8XjUpk0b3X333Ro/fnx1hwMAAUHZC4DX3r17tXLlSiUlJSk/P1/Tp09XZmam+vbtW92hAUDAUPYC4FWrVi2lpaXp6quv1rXXXqsvv/xSH374Ybk3TQOAHVD2AgAAjsLKDwAAcBSSHwAA4CgkPwAAwFFIfgAAgKOQ/AAAAEch+QEAAI5C8gMAAByF5AcAADgKyQ8AAHCU/w+SDrSA1oHLFQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "mask.plot()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "f6e0e2b9-a646-41a2-9b45-f1ccf663224c", + "metadata": {}, + "outputs": [], + "source": [ + "config_obj, save_path, pretrained_model_path = load_config()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "eb97ea83-0742-493e-b41a-3aa51afd5c4c", + "metadata": {}, + "outputs": [], + "source": [ + "# Create dataloader\n", + "dls = get_dls(config_obj, SeasonTST_Dataset, data, mask)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "07be0a63-d531-4d5a-bdf8-d117f1e93a21", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "number of model params 3176450\n", + "weights from saved_models/masked_patchtst/finetuned_d128.pth successfully transferred!\n", + "\n" + ] + } + ], + "source": [ + "# Use the finetuned checkpoint\n", + "path = save_path + config_obj.save_finetuned_model[2:] + \".pth\"\n", + "model = get_model(config_obj, headtype=\"prediction\", weights_path=path, exclude_head=False)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "0d939597-9faf-41af-8126-2d0fb4b247e4", + "metadata": {}, + "outputs": [], + "source": [ + "suggested_lr = 0.00020565123083486514 # Irrelevant as no learning in this notebook\n", + "learner = get_learner(config_obj, dls, suggested_lr, model)" + ] + }, + { + "cell_type": "raw", + "id": "f333bb63-5ad8-4a02-84a5-f91e2e149d6e", + "metadata": {}, + "source": [ + "def inference(learner, data):\n", + " learner.model.eval()\n", + " with torch.no_grad():\n", + " pred = learner.model_forward()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "fab729dd-6811-4eab-96ca-97e3140b8a7b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "score: [array(0.37956834, dtype=float32), array(0.28010777, dtype=float32)]\n" + ] + } + ], + "source": [ + "# Evaluate on test data\n", + "# Pass None as weight_path as weights are already in learner.model\n", + "pred, targ, score = test_func(None, learner, config_obj, dls)" + ] + }, + { + "cell_type": "raw", + "id": "1ea14a02-07f1-4f1e-9f3e-f37eb40d7370", + "metadata": {}, + "source": [ + "score: [array(0.32490477, dtype=float32), array(0.27878124, dtype=float32)]\n", + "score: [array(0.7188467, dtype=float32), array(0.3590506, dtype=float32)]\n" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "6287ffcc-2b08-48af-94cb-fd1fe79453d2", + "metadata": {}, + "outputs": [], + "source": [ + "# Var name sorted as in model output\n", + "var_names = list(dls.valid.dataset.dataset.data_vars.keys())" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "2c06737a-dcfd-4594-96c5-30cf7bdeb112", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'ET0': 0.15539187,\n", + " 'LST_SMOOTHED_5KM': 0.121342644,\n", + " 'NDVI_SMOOTHED_5KM': 0.13158819,\n", + " 'RFH_DEKAD': 1.2999749,\n", + " 'SOIL_MOIST': 0.38951224}" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Scaled metrics\n", + "rmse = np.sqrt(((pred-targ)**2).mean(axis=(0,1)))\n", + "{k: v for k, v in zip(var_names, rmse)}" + ] + }, + { + "cell_type": "raw", + "id": "5c08001e-8eed-4a6d-923e-82e079e9ca17", + "metadata": {}, + "source": [ + "# Scaled metrics\n", + "\n", + "{'ET0': 0.18551289,\n", + " 'LST_SMOOTHED_5KM': 0.15556636,\n", + " 'NDVI_SMOOTHED_5KM': 0.08681418,\n", + " 'RFH_DEKAD': 1.3670056,\n", + " 'SOIL_MOIST': 0.3804626}\n", + "\n", + "{'ET0': 0.14931862,\n", + " 'LST_SMOOTHED_5KM': 0.16489653,\n", + " 'NDVI_SMOOTHED_5KM': 0.15156893,\n", + " 'RFH_DEKAD': 1.1753004,\n", + " 'SOIL_MOIST': 0.41319743}\n", + "\n", + "{'ET0': 0.077757694,\n", + " 'LST_SMOOTHED_5KM': 0.09408079,\n", + " 'NDVI_SMOOTHED_5KM': 0.1081092,\n", + " 'RFH_DEKAD': 2.0208912,\n", + " 'SOIL_MOIST': 0.57550615}\n", + "\n", + "{'ET0': 0.15021749,\n", + " 'LST_SMOOTHED_5KM': 0.15026596,\n", + " 'NDVI_SMOOTHED_5KM': 0.08280951,\n", + " 'RFH_DEKAD': 1.7938881,\n", + " 'SOIL_MOIST': 0.56937885}\n", + "\n", + "{'ET0': 0.10630354,\n", + " 'LST_SMOOTHED_5KM': 0.08450619,\n", + " 'NDVI_SMOOTHED_5KM': 0.1100593,\n", + " 'RFH_DEKAD': 2.038002,\n", + " 'SOIL_MOIST': 0.6122985}\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "df3355a5-cc30-4d3f-8114-8a7cbcfc481b", + "metadata": {}, + "outputs": [], + "source": [ + "# Unscale predictions and targets\n", + "u_pred = pred.copy()\n", + "u_targ = targ.copy()\n", + "for v in range(u_pred.shape[-1]):\n", + " u_pred[:,:,v] = u_pred[:,:,v]*dls.valid.dataset.scaling_factors['std'][var_names[v]] + dls.valid.dataset.scaling_factors['mean'][var_names[v]] \n", + " u_targ[:,:,v] = u_targ[:,:,v]*dls.valid.dataset.scaling_factors['std'][var_names[v]] + dls.valid.dataset.scaling_factors['mean'][var_names[v]] " + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "49ad79d5-cb1e-4a8a-96e1-512d798b742e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'ET0': 0.35864443,\n", + " 'LST_SMOOTHED_5KM': 1.1727766,\n", + " 'NDVI_SMOOTHED_5KM': 0.03459453,\n", + " 'RFH_DEKAD': 38.32326,\n", + " 'SOIL_MOIST': 0.052038837}" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Unscaled metrics\n", + "rmse = np.sqrt(((u_pred-u_targ)**2).mean(axis=(0,1)))\n", + "{k: v for k, v in zip(var_names, rmse)}" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "081f5bc4-1047-4491-9b25-365d18aa0cf7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'ET0': 2.308,\n", + " 'LST_SMOOTHED_5KM': 9.665,\n", + " 'NDVI_SMOOTHED_5KM': 0.2629,\n", + " 'RFH_DEKAD': 29.48,\n", + " 'SOIL_MOIST': 0.1336}" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Compare with std\n", + "dls.valid.dataset.scaling_factors['std']" + ] + }, + { + "cell_type": "raw", + "id": "6498d07f-40cc-4885-bc30-6f4ae7964e5a", + "metadata": {}, + "source": [ + "# Unscaled metrics\n", + "\n", + "{'ET0': 0.4281636,\n", + " 'LST_SMOOTHED_5KM': 1.5035485,\n", + " 'NDVI_SMOOTHED_5KM': 0.022823464,\n", + " 'RFH_DEKAD': 40.299324,\n", + " 'SOIL_MOIST': 0.050829805}\n", + "\n", + "{'ET0': 0.3446275,\n", + " 'LST_SMOOTHED_5KM': 1.593725,\n", + " 'NDVI_SMOOTHED_5KM': 0.039847463,\n", + " 'RFH_DEKAD': 34.647865,\n", + " 'SOIL_MOIST': 0.055203173}\n", + "\n", + "{'ET0': 0.17946891,\n", + " 'LST_SMOOTHED_5KM': 0.90929276,\n", + " 'NDVI_SMOOTHED_5KM': 0.028422236,\n", + " 'RFH_DEKAD': 59.57462,\n", + " 'SOIL_MOIST': 0.076887764}\n", + "\n", + "{'ET0': 0.34670168,\n", + " 'LST_SMOOTHED_5KM': 1.4523225,\n", + " 'NDVI_SMOOTHED_5KM': 0.02177062,\n", + " 'RFH_DEKAD': 52.88384,\n", + " 'SOIL_MOIST': 0.076069005}\n", + "\n", + "{'ET0': 0.24534859,\n", + " 'LST_SMOOTHED_5KM': 0.81675214,\n", + " 'NDVI_SMOOTHED_5KM': 0.028934589,\n", + " 'RFH_DEKAD': 60.080257,\n", + " 'SOIL_MOIST': 0.08180307}\n", + "\n", + "{'ET0': 0.28797588,\n", + " 'LST_SMOOTHED_5KM': 1.0593605,\n", + " 'NDVI_SMOOTHED_5KM': 0.041425053,\n", + " 'RFH_DEKAD': 60.02326,\n", + " 'SOIL_MOIST': 0.08142807}\n", + "\n", + "{'ET0': 0.24534859,\n", + " 'LST_SMOOTHED_5KM': 0.81675214,\n", + " 'NDVI_SMOOTHED_5KM': 0.028934589,\n", + " 'RFH_DEKAD': 60.080257,\n", + " 'SOIL_MOIST': 0.08180307}\n" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "1fe166da-078b-4465-8ddc-dec67b3ce9d9", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(312, 2, 5)" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "u_targ.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "a1b73b45-5934-4c0c-8e85-ddddaa2479d8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "312" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(dls.test.dataset)" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "id": "4ecd22a8-38a2-4f7b-afda-3f2f3c273f8d", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGxCAYAAABBZ+3pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABtbUlEQVR4nO3deVxU9foH8M/MwMywDiA7IuKCorjiiuKSiJlatphWahrert26pXZvZdatbLHVn95KbylkWqaVVuZSYikuuOKuuCAqKIMswoAsA8yc3x/DTCKgDAJnls/79Tqv9Mx3zjzHQ87jd3m+EkEQBBARERFZMKnYARARERHdCRMWIiIisnhMWIiIiMjiMWEhIiIii8eEhYiIiCweExYiIiKyeExYiIiIyOIxYSEiIiKLx4SFiIiILB4TFrIrK1asgEQigVKpxOXLl2u9PmzYMERERNQ417ZtW0gkEkgkEkilUqhUKoSHh2Pq1KnYunVrjbazZ8+GRCLBmTNn6o1h3rx5kEgkOHz4sOn6Y8eONes+BEHAmjVrEB0dDV9fXyiVSrRu3RqjRo3C8uXLa7Q1xj5t2rQ6rzV//nxTm0uXLtX6nNWrV+Oee+6Bp6cnFAoF2rVrh2effRaZmZn1xtaQ9+zYscP0uXc6gL+e3aFDh+r83LFjx6Jt27Y1zt387G49hg0bZmpnvLbxUCqV8Pf3x/Dhw7FgwQLk5OTU+Zl3Ut9nv//++zXavfnmm5BIJMjLy6tx/sKFC2jXrh38/Pxw9OhRAMC0adMgkUjg5uaGGzdu1PrMy5cvQyqVQiKR4M0332xU3MZ4bj2USmWttsXFxXj++ecRFBQEhUKBsLAwfPjhh9DpdI36bKL6OIgdAJEYtFotXnvtNaxatapB7QcNGoSPP/4YAHDjxg2cPXsWa9aswahRo/Dwww/ju+++g6OjI+Li4rBo0SIkJCTgww8/rHUdvV6PlStXomfPnujdu3ej4587dy4++OAD/O1vf8O///1vuLm54fLly/jzzz/xyy+/YMaMGTXau7m54YcffsCnn34KNzc303lBELBixQq4u7ujqKioVqyPP/441q5di8ceewwrVqyASqXC8ePH8dFHH2H16tXYuHEjBg0a1Kj39O7dG3v37q3xmQ8++CDat29v+rNuCjc/u5u5u7vXOvfVV1+hc+fOqKysRE5ODnbv3o0PPvgAH3/8MdauXYuYmBizP/+RRx7Biy++WONcmzZt7vi+EydOYNSoUXB0dMTu3bvRsWNH02uOjo6oqqrC2rVrERcXV+se3Nzcaj3Pxvjtt9+gUqlMv5dKa/4bt6qqCiNHjsS5c+fw9ttvIywsDL/99hteeeUVXLlyBf/973/vOgYiE4HIjnz11VcCAOHee+8VpFKpcPTo0RqvDx06VOjatWuNcyEhIcKYMWPqvN4bb7whABBeeukl07l+/foJ/v7+QmVlZa32W7ZsEQAIn376aYOuX5fS0lJBoVAIU6dOrfN1nU5X4/cAhMmTJwtOTk7Cl19+WeO1bdu2CQCEv/3tbwIA4eLFi6bX3nvvPQGA8P7779f6jOzsbCEkJETw8/MTCgoK7uo9N7vdn4Xx2R08eLDO18eMGSOEhIQ0+HoNvfbly5eF4OBgwc3NTcjOzr7jtW4GQHj22Wfv2M74c5SbmysIgiDs3btX8PT0FDp37ixkZmbWaPvkk08KLi4uwqRJk4SoqKgar+n1eiEkJMT0PN944w2z4q0vnvp89913AgBh3bp1Nc4//fTTglQqFc6cOdOozyeqC4eEyC699NJLaNWqFV5++eW7us6bb76Jrl274rPPPkN5eTkAIC4uDtnZ2diyZUut9l999RUUCgWeeOKJRn9mSUkJtFotAgIC6nz91n8FA4BKpcKDDz6IhISEGucTEhIwaNAghIWF1ThfUVGBjz76COHh4XjppZdqXc/Pzw8LFizAtWvXEB8f3+j3WIM2bdrgk08+QXFxMb744otm/7zExETExMSgffv22LVrF1q3bl1nu6eeegrJyck4e/as6dy2bdtw+fJlTJ8+vdnjBIA9e/ZAIpFg9OjRNc6PHTsWer0eP/30U4vEQfaBCQvZJTc3N7z22mv4/fff8eeff97VtcaNG4fS0lLT3IrHHnsMzs7OtZKDgoIC/PLLL3jwwQfh6enZ6M/z9vZGhw4dsGTJEixcuBBnzpyB0IBN1+Pi4rBv3z6kpqYCAAoLC7F+/fpaQwoAkJKSgoKCAtx///2mOSS3GjduHKRSKRITExv9nsbQ6XSoqqqqddT3ZyAIglnt63LfffdBJpNh586dZse7evVqODk5QaFQIDIyEl999VW9bdetW4exY8eib9+++PPPP+Ht7V1v25iYGISEhNT4OYuPj8eQIUNqDB/djW7dukEmk8HPzw9Tp05FRkZGjdcrKioglUrh6OhY47xCoQAAHD9+vEniIAKYsJAdmzlzJtq1a4eXX37ZrC+vW4WEhAAAsrKyABiSoQkTJmDTpk01Jmt+++230Gq1dSYI5lq9ejU8PT3x4osvIjw8HCqVCuPGjcOqVavqvZfhw4cjNDTU9AW3evVqODg4YMKECbXaGr+YQkND643B1dUVPj4+praNeU9jDBgwAI6OjrWOzZs319l+8+bNdbZ/9913G/yZLi4u8Pb2Nj3jhnr88cfx2WefYevWrVi9ejX8/Pzw1FNP4fXXX6+z/cyZM9G6dWts2bKlxlyjuhgnUq9cuRJVVVW4fv06fv75Zzz11FNmxViX9u3b491330VCQgK2bduGOXPmYNOmTejXrx+uXr1qatelSxfodDrs27evxvt3794NAMjPz7/rWIiMmLCQ3ZLL5XjnnXdw6NAhfP/9942+Tl0JQlxcHCorK2tM6v3qq68QEhKCESNGNPqzjPr27Yu0tDT89ttvePXVVzFw4ED88ccfmDp1Ku6///46YzJ+wa1atQpVVVWIj4/Ho48+CldX10bHIQhCvb0pTfmem61cuRIHDx6sdQwePLjO9oMHD66zvbmJY2OS2m+//RaPP/44oqOj8fDDD2Pz5s0YO3Ys3n//feTm5tZqf//99yM9Pb3Bq3umT5+Oa9euYcuWLfj2228hl8vrTEDNNWXKFLz66qsYPXo0hg8fjpdffhlbtmxBbm5ujcnkTzzxBLy8vPD0009j//79KCwsxHfffWeabFvX8CRRY/GniezapEmT0Lt3b8ybNw+VlZWNuoZxeXRgYKDpXHR0NMLCwkzd/8ePH8fhw4cxffr0u/qyvpmjoyNGjRqFd999F7///jsyMzMxbNgwbNy4sc75M4DhCy43NxfvvfceDh8+XO+XtnEVy8WLF+v9/JKSEuTl5SE4OLjR72mM8PBw9OnTp9Zx82qWm6lUqjrb1zcHqL648/Pzazzjxpo8eTKqqqrqXJ69bNkyTJs2DR988EGd84BuZUyAExISkJCQgEmTJsHZ2fmuY6xLv379EBYWVqM3xdvbG7/99hsAQ8+Xp6cn/vnPf2LhwoUAgKCgoGaJhewTExayaxKJBB988AEuXLiAL7/80uz3C4KAX3/9FS4uLujTp0+N15566imcOnUKBw4cQEJCAqRSab21UJpCq1atMGvWLADAyZMn62wTHByMmJgYvPXWW+jUqROioqLqbBcZGQlPT09s2LCh3p6FDRs2QK/XY+TIkY1+j7XYtGkTdDpdjdotjWX8s6mr90EqlSI+Ph7Tp0/HRx99hH/96193vN5TTz2FDRs24OjRo00yHHQ7giDUirtv3744ffo0Ll68iJMnTyIrKwvh4eEAgCFDhjRrPGRfmLCQ3YuJicHIkSMxf/78Ogtx3c5bb72F06dP44UXXqhVVOvJJ5+Eg4MDvvjiC3z77bcYMWKEab7L3aisrKx3boBxQu3tegJefPFFjBs3rt55FIBhuOzf//43UlNT8dFHH9V6PScnB3PnzoWfn5+p5ktj3mMNMjIy8K9//QsqlQp///vf7/p6q1atgqOjIyIjI+t83Zi0zJgxA5988gnmzJlz2+s9+OCDePDBB/HUU09hwIABdx1fffbt24fz58/X+xlt27ZF165d4ejoiE8++QSBgYFNMjxFZMTCcUQAPvjgA0RGRiInJwddu3at9XphYaGpK7ykpMRUOG7Xrl149NFH8dZbb9V6j7+/P+677z589dVXEAShSSbbAoBGo0Hbtm0xYcIExMTEIDg4GDdu3MCOHTuwePFihIeH46GHHqr3/bGxsYiNjb3j57z88ss4duyY6b8TJ06sUQSuuLgYGzdurDEU05j3NLebn93NFAoFevXqVePcyZMnTauIcnJysGvXLnz11VeQyWT46aef4OPj0+DP/eijj3D69GmMGDECrVu3Rk5ODuLj47F161a8+eabt10BJJFI8OWXX0IikeD//u//IAgC/u///q/OtkqlEj/++GODYjJWAr61ovGtevTogcmTJyM8PBxKpRIHDhzARx99BH9//1pDVfPmzUO3bt0QEBCAjIwMJCQkYP/+/di0aROcnJwaFBdRQzBhIQLQq1cvPPbYY1i9enWdr+/ZswcDBw6ERCKBi4sLgoKC0K9fP7z22mu3/fKPi4vDhg0b4OXlhfHjxzdJrO7u7njrrbfwxx9/4NVXX8W1a9cgkUgQGhqKWbNm4eWXX26SeQxSqRTfffcd7r//fixbtgxPPvkkSktLERQUhLFjx+KVV16pVbG1Me9pbsZnd6ugoCBcuXKlxjlj/RK5XA4PDw+Eh4fj5ZdfxowZM8xKVgCgc+fO2LBhAzZt2oSCggI4OTmhZ8+e+O677zBp0qQ7vl8ikeCLL76ATCbDokWLoNfrsXjxYrNiuFVJSQk6dOhwx3ZdunTBl19+CbVajYqKCgQGBmLSpEn4z3/+U2vuT0FBAV5++WVkZ2fD3d0dQ4cOxf79+9GtW7e7ipXoVhLhbtZzEhGRVTh9+jS6du2KjRs3YsyYMWKHQ2Q2zmEhIrID27dvx8CBA5mskNViDwuRBdHpdLet9yGRSCCTyVowIrqVIAh33IlYJpM12fJ1IjJgDwuRBRkxYkSdVVmNR/v27cUO0e4lJSXd9hk5Ojri66+/FjtMIpvDHhYiC3L27FkUFxfX+7pCoeBkRpEVFxfX2HCwLqGhoWjVqlULRURkH5iwEBERkcXjkBARERFZPJupw6LX65GVlQU3NzdOdiMiIrISgiCguLgYgYGBt90w02YSlqysrLvaUI2IiIjEk5mZidatW9f7us0kLG5ubgAMN+zu7i5yNERERNQQRUVFCA4ONn2P18dmEhbjMJC7uzsTFiIiIitzp+kcnHRLREREFo8JCxEREVk8JixERERk8ZiwEBERkcVjwkJEREQWjwkLERERWTwmLERERGTxmLAQERGRxWPCQkRERBaPCQsRERFZPCYsREREZPGYsBAREZHFY8LSxM5mF+PTP84jW1MudihEREQ2w2Z2axabtkqHz7dfwNIdaajUCfhm/2Usn9oX3VqrxA6NiIjI6rGHpQkczijA2P/uxn//OI9KnQB3pQOuFWnx6Bd7sfVUttjhERERWT0mLHehRFuFt349hYeXJuN8zg14u8rx2eO9sPuVexDd0RtllTr8/ZsULNuZDkEQxA6XiIjIajFhaaSd53IR+3878dWeSxAE4KHeQUicPRRjuwfCXemIr6b1xRP920AQgHc3p+LVn06iUqcXO2wiIiKr1KiEZcmSJQgNDYVSqURkZCR27dp12/bffvstevToAWdnZwQEBGD69OnIz8+v0WbdunXo0qULFAoFunTpgp9++qkxoTW7wtIK/OuHY5iacABXC8sQ5OGEr5/qh4WP9oSni9zUzkEmxTvjI/D62C6QSIDvDmRg+lcHoSmrFDF6IiIi62R2wrJ27VrMmjUL8+bNw5EjRxAdHY3Ro0cjIyOjzva7d+/G1KlTERcXh1OnTuGHH37AwYMHMWPGDFObvXv3YuLEiZgyZQqOHTuGKVOm4NFHH8X+/fsbf2dNTBAEbD6hRszCnfgx5QokEmBaVFtsnT0EQ8N86nyPRCJB3OBQLJvSB85yGXan5eHhpcnIyC9t4eiJiIism0Qwc3JF//790bt3byxdutR0Ljw8HOPHj8eCBQtqtf/444+xdOlSXLhwwXTu008/xYcffojMzEwAwMSJE1FUVIQtW7aY2tx7773w9PTEd99916C4ioqKoFKpoNFo4O7ubs4t3dG1onK8/vNJbD19DQDQwdcVHzzcHZEhng2+xsmrGsz4+hCyi8rh5SLHsqmRiAzxatI4iYiIrE1Dv7/N6mGpqKhASkoKYmNja5yPjY1FcnJyne+JiorClStXsHnzZgiCgGvXruHHH3/EmDFjTG327t1b65qjRo2q95oAoNVqUVRUVONoaoIgYO3BDMQsTMLW09fgIJXg+Xs6YNPzg81KVgAgIkiFn58dhK6B7rheUoHHlu3HL0evNnnMREREtsishCUvLw86nQ5+fn41zvv5+SE7u+7lu1FRUfj2228xceJEyOVy+Pv7w8PDA59++qmpTXZ2tlnXBIAFCxZApVKZjuDgYHNupUG0VXr8LykdxeVV6N5ahV//ORhzYjtB4SBr1PX8VUr8MHMgRnbxQ0WVHi+sOYrF285zBREREdEdNGrSrUQiqfF7QRBqnTM6ffo0nn/+efznP/9BSkoKfvvtN1y8eBEzZ85s9DUBYO7cudBoNKbDOLzUlJSOMrz/UDfMuy8c65+JQnjA3Q81Ocsd8L/JkfhbdCgA4P+2ncOc749BW6W762sTERHZKrMq3Xp7e0Mmk9Xq+cjJyanVQ2K0YMECDBo0CP/+978BAN27d4eLiwuio6PxzjvvICAgAP7+/mZdEwAUCgUUCoU54TdK/3at0L9dqya9pkwqwbwxXRDq7YrXfzmJn45chcJBivcf7t6kn0NERGQrzOphkcvliIyMRGJiYo3ziYmJiIqKqvM9paWlkEprfoxMZhhSMQ6FDBw4sNY1t27dWu81bcXj/dtg2dRIAMAPKVeQeZ2rh4iIiOpi9pDQnDlzsHz5ciQkJCA1NRWzZ89GRkaGaYhn7ty5mDp1qqn9uHHjsH79eixduhTp6enYs2cPnn/+efTr1w+BgYEAgBdeeAFbt27FBx98gDNnzuCDDz7Atm3bMGvWrKa5Swt2T2c/RHf0hk4v4H9JF+78BiIiIjtk9uaHEydORH5+PubPnw+1Wo2IiAhs3rwZISEhAAC1Wl2jJsu0adNQXFyMzz77DC+++CI8PDxwzz334IMPPjC1iYqKwpo1a/Daa6/h9ddfR/v27bF27Vr079+/CW7R8j03vAN2nc/DD4eu4J/3dIS/Sil2SERERBbF7Doslqo567C0hAn/S8bBSwWIGxyK18d2ETscIiKiFtEsdVio+Tw7vAMAYPX+DOTf0IocDRERkWVhwmIhhob5oFuQCmWVOiTsuSh2OERERBaFCYuFkEgkpl6WlcmXuUkiERHRTZiwWJDYLn4I83NFsbYKq/ZeEjscIiIii8GExYJIpX/1ssTvvogSbZXIEREREVkGJiwWZky3AIS0ckZBaSW+O5Bx5zcQERHZASYsFsZBJsUzQ9sDAL7cmY7ySu4xRERExITFAj3UuzUCVErkFGvxY8oVscMhIiJbo9UCt+zhZ+mYsFgguYMUfx/SDgCwdMcFVOr0IkdEREQ24/JlYPBgYOxYQ+JiJZiwWKhJ/drA21WOq4Vl+OVoltjhEBGRrZBKgYsXgZQU4MUXxY6mwZiwWCilowxxgw29LEt2pEGnt4kdFIiISGzBwcCqVYZff/458P334sbTQExYLNjkAW2gcnJEem4JtpxUix0OERHZitGjgVdeMfx6xgzg/Hlx42kAJiwWzE3piGlRbQEAn2+/ABvZp5KIiJpYVmGZ+W96+20gOhooLgYefRQoL2/6wJoQExYLN31QW7jIZUhVF+HPMzlih0NERBYmp6gcwz/egakJB3DDnIKjDg7Ad98BPj7A0aPArFnNFWKTYMJi4Tyc5Zg8MAQA8OmfaexlISKiGpYmXYC2So9SbRVc5DLz3hwUBHzzDSCRAF98YUhgLBQTFiswY3A7KBykOJpZiOQL+WKHQ0REFiKnqByr9xuqor8Q0xESicT8i8TGAvPmGX799NPA2bNNGGHTYcJiBXzcFJjUNxgA8NmfaSJHQ0REluJ/SenQVukRGeKJwR28G3+hN98Ehg0DbtwwzGcpa8ScmGbGhMVKPD20PRxlEuxNz0fK5etih0NERCLLKS7Ht/svAwBeGNHI3hUjmQxYvRrw9QWOHweef76Jomw6TFisRJCHEx7q1RoAe1mIiAj4orp3pXcbD0R3vIveFaOAAEPSIpEAy5cb5rZYECYsVuSZYe0hlQDbz+bi5FWN2OEQEZFIavSuxITdXe/KzUaMAP7zH8Ov//53IDW1aa7bBJiwWJG23i4Y2z0QgGGPISIisk9fJqWjvFKPXm08MKQpeldu9vrrhsSltBSYMAEoKWna6zcSExYr8/ehhnL9iaevocSc9fZERGQTcou1+Kap5q7URSYDvv0W8PcHTp0Cnnuuaa/fSExYrEyXAHe08XJGhU6P3Wl5YodDREQt7MudF1BeqUfPYA8MDfNpng/x8zPUZJFKgRUrDIfImLBYGYlEgns6+wIA/kxl5VsiInuSd0OLVfuMc1eaoXflZsOGAW+9Zfj1P/5h6G0RERMWKzQivDphOZsDPXdxJiKyG1/uNMxd6dFahWHN1btys1dfBUaONNRlmTDBUKdFJExYrFC/UC+4yGXILdbiZBZXCxER2YO8G1qs2mvoXZnVlCuDbkcqNSxvDgw0rBj64Yfm/8z6QhHtk6nRFA4yRHc0ZNZ/cFiIiMguLNuZjrJKHbq3VmFYpxboXTHy9QXWrAFWrQKmT2+5z70FExYrdY9xWIg7OBMR2bz8G1qsNPWuNPPclbpERwOTJ7fsZ96CCYuVGt7JkLCcuKrBtaJykaMhIqLm9OWuv3pXjH//2xsmLFbKx02BHsEeAIDt7GUhIrJZ+TfNXWmWuitWolEJy5IlSxAaGgqlUonIyEjs2rWr3rbTpk2DRCKpdXTt2tXUprKyEvPnz0f79u2hVCrRo0cP/Pbbb40Jza6MqF7e/AcTFiIim7Vs10WUVujQLUhlKmthj8xOWNauXYtZs2Zh3rx5OHLkCKKjozF69GhkZGTU2X7x4sVQq9WmIzMzE15eXpgwYYKpzWuvvYYvvvgCn376KU6fPo2ZM2fiwQcfxJEjRxp/Z3bA+IO7+3weyit1IkdDRERN7XpJBVbuvQTAvntXgEYkLAsXLkRcXBxmzJiB8PBwLFq0CMHBwVi6dGmd7VUqFfz9/U3HoUOHUFBQgOk3zTRetWoVXn31Vdx3331o164dnnnmGYwaNQqffPJJ4+/MDnQNdIefuwJllTrsS88XOxwiImpiy3alo7RCh4ggd1MNLntlVsJSUVGBlJQUxMbG1jgfGxuL5OTkBl0jPj4eMTExCAkJMZ3TarVQKpU12jk5OWH37t31Xker1aKoqKjGYW8MVW/9AHC1EBGRrbleUoGVyZcAAC+MaKG6KxbMrIQlLy8POp0Ofn5+Nc77+fkhOzv7ju9Xq9XYsmULZsyYUeP8qFGjsHDhQpw/fx56vR6JiYn45ZdfoFar673WggULoFKpTEdwcLA5t2IzTPNYUnMgCKx6S0RkK5bvSkdJhQ5dA90RY+e9K0AjJ93emuUJgtCgzG/FihXw8PDA+PHja5xfvHgxOnbsiM6dO0Mul+O5557D9OnTIZPJ6r3W3LlzodFoTEdmZmZjbsXqDergDYWDFFcLy3Dumnglk4mIqOkUlFTga1Pvin3PXTEyK2Hx9vaGTCar1ZuSk5NTq9flVoIgICEhAVOmTIFcLq/xmo+PD37++WeUlJTg8uXLOHPmDFxdXREaGlrv9RQKBdzd3Wsc9shJLkNU+1YAgD/OXBM5GiIiagrLdxt6V7oEuGNkl9t/v9oLsxIWuVyOyMhIJCYm1jifmJiIqKio2743KSkJaWlpiIuLq7eNUqlEUFAQqqqqsG7dOjzwwAPmhGe37gmvnsfCMv1ERFbP0LvSQjsyWxGzh4TmzJmD5cuXIyEhAampqZg9ezYyMjIwc+ZMAIahmqlTp9Z6X3x8PPr374+IiIhar+3fvx/r169Heno6du3ahXvvvRd6vR4vvfRSI27J/hiXNx/OKMD1kgqRoyEiorux7vAV3NBWITzAHbHsXTFxMPcNEydORH5+PubPnw+1Wo2IiAhs3rzZtOpHrVbXqsmi0Wiwbt06LF68uM5rlpeX47XXXkN6ejpcXV1x3333YdWqVfDw8DD/juxQkIcTOvu74Ux2MZLO5eDBXq3FDomIiBop6VwuAOCRyNbsXbmJRLCRpSVFRUVQqVTQaDR2OZ/lo9/P4PPtFzC2ewA+e7y32OEQEVEjlFXo0GP+VlRU6bFtzhB08HUTO6Rm19Dvb+4lZCOM9ViSzuWiUqcXORoiImqMfRfzUVGlR5CHE9r7uIodjkVhwmIjegZ7wMtFjuLyKhy6VCB2OERE1AhJZw3DQUPCfDgcdAsmLDZCJpVgWCcfAMCfXN5MRGSVdlbPXxka5i1yJJaHCYsNGVE9LMTdm4mIrE/m9VKk55VAJpUgqgMTllsxYbEh0WHecJBKkJ5bgkt5JWKHQ0REZjCuDops4wl3paPI0VgeJiw2xF3piH6hXgC4GSIRkbUxJixDOBxUJyYsNsZYRI4JCxGR9aio0mPvhXwAwNAwbnRYFyYsNmZEdZn+/RfzUVxeKXI0RETUEIczCnBDW4VWLnJ0DbS/WmINwYTFxoR6u6CdtwsqdQJ2n88TOxwiImoA4+qg6I7ekEq5nLkuTFhskHFYiKuFiIisg3H+ytDq8hRUGxMWG3RPuCFh2X4mB3q9Tey8QERks3KLtTiVVQQAiO7IhKU+TFhsUN+2XnBTOCC/pALHrhSKHQ4REd3GrvOG3pWIIHd4uypEjsZyMWGxQY4yKYaYqt5yWIiIyJKZhoPC2LtyO0xYbNQI4zyWVCYsRESWSq8XsKt6gcQQDgfdFhMWGzWsky8kEuC0ughqTZnY4RARUR1OZmlwvaQCrgoH9A7xFDsci8aExUZ5ucjRu43hh5/DQkRElsm4O/OgDq3gKONX8u3wT8eGmarecliIiMgi7TxvLMfP4aA7YcJiw0ZUL2/enZaHsgqdyNEQEdHNNGWVOJxRCIDzVxqCCYsN6+TnhiAPJ2ir9Nibzqq3RESWJDktDzq9gPY+Lgj2chY7HIvHhMWGSSSSv6recliIiMiicDjIPExYbJyx6u2fZ3IgCKx6S0RkCQRBME24Zf2VhmHCYuMGtmsFpaMUak05zl4rFjscIiICkJZzA1macsgdpOgf2krscKwCExYbp3SUYWA7w/8MxmyeiIjEZaxu2z/UC05ymcjRWAcmLHbA2N24gwkLEZFFYDl+8zFhsQPDOhnmsRy6fB03tFUiR0NEZN/KK3U4cPE6ACYs5mDCYgfaersgpJUzKnUC9l7IFzscIiK7ti89H9oqPQJVSnTwdRU7HKvBhMVOGLP4pHNc3kxEJCbTcFAnH0gkEpGjsR5MWOzEsE5/zWPh8mYiIvHsrE5YWN3WPExY7MSAdq0gl0lxpaAM6XklYodDRGSXMq+X4kJuCWRSCaI6eIsdjlVhwmInnOUO6BfqBYDLm4mIxGKsbtu7jQdUTo4iR2NdmLDYEdOw0DkmLEREYuBwUOM1KmFZsmQJQkNDoVQqERkZiV27dtXbdtq0aZBIJLWOrl271mi3aNEidOrUCU5OTggODsbs2bNRXl7emPCoHsaJt/vT81Feyd2biYhaUqVOjz1phpWaQzsxYTGX2QnL2rVrMWvWLMybNw9HjhxBdHQ0Ro8ejYyMjDrbL168GGq12nRkZmbCy8sLEyZMMLX59ttv8corr+CNN95Aamoq4uPjsXbtWsydO7fxd0a1dPB1RaBKCW2VHvvSubyZiKglHckoxA1tFbxc5IgIVIkdjtUxO2FZuHAh4uLiMGPGDISHh2PRokUIDg7G0qVL62yvUqng7+9vOg4dOoSCggJMnz7d1Gbv3r0YNGgQHn/8cbRt2xaxsbF47LHHcOjQoXrj0Gq1KCoqqnHQ7UkkEgytLiLHqrdERC3LWFYiuqM3pFIuZzaXWQlLRUUFUlJSEBsbW+N8bGwskpOTG3SN+Ph4xMTEICQkxHRu8ODBSElJwYEDBwAA6enp2Lx5M8aMGVPvdRYsWACVSmU6goODzbkVu2UcFtrJeSxERC2K5fjvjoM5jfPy8qDT6eDn51fjvJ+fH7Kzs+/4frVajS1btmD16tU1zk+aNAm5ubkYPHgwBEFAVVUVnnnmGbzyyiv1Xmvu3LmYM2eO6fdFRUVMWhogqkMrOEglSM8rQUZ+Kdq0chY7JCIim5d3Q4uTVw0jAdGccNsojZp0e2tlPkEQGlStb8WKFfDw8MD48eNrnN+xYwfeffddLFmyBIcPH8b69euxceNGvP322/VeS6FQwN3dvcZBd+audETvEE8ArHpLRNRSdlUvZ+4a6A4fN4XI0VgnsxIWb29vyGSyWr0pOTk5tXpdbiUIAhISEjBlyhTI5fIar73++uuYMmUKZsyYgW7duuHBBx/Ee++9hwULFkCv15sTIjWAcXlzEoeFiIhahLH+FYeDGs+shEUulyMyMhKJiYk1zicmJiIqKuq2701KSkJaWhri4uJqvVZaWgqptGYoMpkMgiCwjHwzMP4Pk3whH9oqLm8mImpOer2AXefzAABDmLA0mllzWABgzpw5mDJlCvr06YOBAwfiyy+/REZGBmbOnAnAMLfk6tWrWLlyZY33xcfHo3///oiIiKh1zXHjxmHhwoXo1asX+vfvj7S0NLz++uu4//77IZPJGnlrVJ8uAYYuydxiLVIuFbA8NBFRMzqVVYT8kgq4KhzQu42n2OFYLbMTlokTJyI/Px/z58+HWq1GREQENm/ebFr1o1ara9Vk0Wg0WLduHRYvXlznNV977TVIJBK89tpruHr1Knx8fDBu3Di8++67jbgluhOJRIKhYT74MeUKdpzLZcJCRNSMjOX4o9q3gtyBBeYbSyLYyJhLUVERVCoVNBoNJ+A2wK/HsvDP746gk58bfp89ROxwiIhs1uTl+7E7LQ9vP9AVUwa2FTsci9PQ72+menYquqM3pBLg7LViqDVlYodDRGSTqnR6HM4oAAD0rd6AlhqHCYud8nCWo2ewBwDu3kxE1FxOZRWhtEIHlZMjwnzdxA7HqjFhsWNDwwxl+rm8mYioeRy8dB0A0CfEk+X47xITFjtm3C109/k8VOpY74aIqKmZEpa2HA66W0xY7Fj3IBW8XOQo1lbhSEah2OEQEdkUQRBw6JJh/kq/UC5nvltMWOyYVCpBdEfDkmaW6ScialrpeSXIL6mAwkGKiCCV2OFYPSYsds5Y9ZbzWIiImtbBi4bhoB7BHlA4sAjq3WLCYueMZaJPXi1CTnG5yNEQEdmOA9XzV/px/kqTYMJi57xdFehW3VW561yeyNEQEdkO4/wV1l9pGkxYiMNCRERN7FpROTKul0IqAXq38RA7HJvAhIVMy5t3nc+FTm8TOzUQEYnKuJw5PMAdbkpHkaOxDUxYCL2CPeCmdEBBaSWOXykUOxwiIqtnnHDbl/NXmgwTFoKDTHrT8mYOCxER3a0DxvkrTFiaDBMWAsB5LERETaWovBJnsosAAH1ZMK7JMGEhAH/tK3Q0sxAFJRUiR0NEZL1SLhdAEIC2rZzh66YUOxybwYSFAAD+KiU6+7tBEIBdaVzeTETUWMb5K9w/qGkxYSET07DQWQ4LERE1lmn/ICYsTYoJC5kYlzcnncuFnsubiYjMpq3S4Wj1aksWjGtaTFjIpE+IF5zlMuTd0OK0ukjscIiIrM7xKxpUVOnh7SpH21bOYodjU5iwkIncQYqo9lzeTETUWMaCcX3bekEikYgcjW1hwkI1DOvEeSxERI3FCbfNhwkL1WCceJuSUYCi8kqRoyEish46vYBDlznhtrkwYaEagr2c0d7HBTq9gO1ncsQOh4jIapy7Vozi8iq4yGUID3ATOxybw4SFarmvWwAA4NdjapEjISKyHsb5K71DPOEg49drU+OfKNUytnsgAGDnuVxoyjgsRETUEAe44WGzYsJCtXTyd0OYnysqdHoknr4mdjhERBZPEIQaK4So6TFhoToZe1l+PZYlciRERJbvSkEZrhVp4SiToGewh9jh2CQmLFSnsd0N81j2pOXhOjdDJCK6LeNwUESQCk5ymcjR2CYmLFSndj6u6Brojiq9gN9OZosdDhGRRTt02ZCwcDlz82HCQvUa18MwLLTxOIeFiIhu5wALxjW7RiUsS5YsQWhoKJRKJSIjI7Fr1656206bNg0SiaTW0bVrV1ObYcOG1dlmzJgxjQmPmsiY6uXN+9LzkVNcLnI0RESWKf+GFhdySwAAfUI8RY7GdpmdsKxduxazZs3CvHnzcOTIEURHR2P06NHIyMios/3ixYuhVqtNR2ZmJry8vDBhwgRTm/Xr19doc/LkSchkshptqOUFezmjVxsP6AVgywkOCxER1cVY3TbMzxWeLnKRo7FdZicsCxcuRFxcHGbMmIHw8HAsWrQIwcHBWLp0aZ3tVSoV/P39TcehQ4dQUFCA6dOnm9p4eXnVaJOYmAhnZ2cmLBaAq4WIiG6P+we1DLMSloqKCqSkpCA2NrbG+djYWCQnJzfoGvHx8YiJiUFISMht20yaNAkuLi71ttFqtSgqKqpxUNMb0y0AEonhXxBZhWVih0NEZHGM9Vc44bZ5mZWw5OXlQafTwc/Pr8Z5Pz8/ZGffechArVZjy5YtmDFjRr1tDhw4gJMnT962DQAsWLAAKpXKdAQHBzfsJsgs/iqlqQjS5hMs1U9EdLPSiiqczDL8g7lvKBOW5tSoSbcSiaTG7wVBqHWuLitWrICHhwfGjx9fb5v4+HhERESgX79+t73W3LlzodFoTEdmZmaDYifzjetu3FuIw0JERDc7klEInV5AoEqJIA8nscOxaWYlLN7e3pDJZLV6U3Jycmr1utxKEAQkJCRgypQpkMvrnpRUWlqKNWvW3LF3BQAUCgXc3d1rHNQ8RncLgFQCHLuiQUZ+qdjhEBFZDNP+QexdaXZmJSxyuRyRkZFITEyscT4xMRFRUVG3fW9SUhLS0tIQFxdXb5vvv/8eWq0WkydPNicsambergpEtfcGAPzKmixERCbGgnHcP6j5mT0kNGfOHCxfvhwJCQlITU3F7NmzkZGRgZkzZwIwDNVMnTq11vvi4+PRv39/RERE1Hvt+Ph4jB8/Hq1atTI3LGpm43oYhoU2Huc8FiIiAKjU6XH4ciEAJiwtwcHcN0ycOBH5+fmYP38+1Go1IiIisHnzZtOqH7VaXasmi0ajwbp167B48eJ6r3vu3Dns3r0bW7duNTckagGjuvpj3k8nkaouQlrODXTwdRU7JCIiUZ3KKkJZpQ4qJ0d05N+Jzc7shAUA/vGPf+Af//hHna+tWLGi1jmVSoXS0tvPfQgLC4MgCI0Jh1qAh7Mc0R29sf1sLjYez8KsmDCxQyIiEtWh6uXMfUI8IZXeeeEJ3R3uJUQNZtxb6NdjWUwuicjuccJty2LCQg02sosf5A5SXMgtwZnsYrHDISISjSAIppL8nL/SMpiwUIO5KR0xvJMPAO7gTET27UJuCa6XVEDhIEW3IJXY4dgFJixklr/2FlJzWIiI7JaxHH/PYA/IHfhV2hL4p0xmGRHuCydHGTKul+LEVY3Y4RARicK44WE/zl9pMUxYyCzOcgeMCPcFwFL9RGS/Dl7mDs0tjQkLmc24WmjTcTX0eg4LEZF9ydaUI/N6GaQSoHcbD7HDsRtMWMhsQ8N84KZwQJamHIczCsQOh4ioRR2onr/SJdAdbkpHkaOxH0xYyGxKRxlGdjFsdslS/URkb/4qGMfhoJbEhIUaxTQsdEINHYeFiMiOHOCEW1EwYaFGGdTBGyonR+QWa7H/Yr7Y4RARtYi8G1pT4cz+TFhaFBMWahS5gxSjI/wBGGqyEBHZg70XDP9A6+zvhlauCpGjsS9MWKjRjEXkfjupRqVOL3I0RETNL/lCHgBDLzO1LCYs1GgD2nnB21WOgtJKJF/gsBAR2T7j33VR7VuJHIn9YcJCjeYgk2J0RAAAFpEjItt3paAUl/NLIZNKOOFWBExY6K4YVwv9fiob2iqdyNEQETUfY+9K99Yq1l8RARMWuit9Qjzh765EcXkVks7mih0OEVGzSU6rnr/SnvNXxMCEhe6KVCrBuB6GYaHvD2WKHA0RUfMQBOGv+SsdOH9FDExY6K5N6tcGAPDnmRxkFZaJHA0RUdO7kHsDOcVaKByk6N3GU+xw7BITFrpr7X1cMaCdF/QCsPYge1mIyPbsSTP0rvRp6wmlo0zkaOwTExZqEo/3DwFgSFiqWJOFiGyMsf5KFOeviIYJCzWJUV394OUiR3ZRObZz8i0R2RCdXjBVuGX9FfEwYaEmoXCQYUJkawDA6v2XRY6GiKjpnM4qQlF5FdwUDugWpBI7HLvFhIWazGPVk293nMtF5vVSkaMhImoae6qHg/q384KDjF+bYuGfPDWZtt4uGNShFQROviUiG2JczjyQ81dExYSFmtTj/aon3x7K5IaIRGT1Kqr0OHjxOgBgEOuviIoJCzWpkV384O0qR26xFn+kXhM7HCKiu3I0sxBllTq0cpEjzNdN7HDsGhMWalJyBykm9AkGAHy7P0PkaIiI7s6e6nL8A9u3glQqETka+8aEhZrcY30Nk293nc9DRj4n3xKR9WL9FcvBhIWaXJtWzojuaPif+7uD7GUhIutUWlGFIxmFADh/xRIwYaFm8UR/Qy/LD4cyUVHFybdEZH0OXLyOKr2AIA8ntPFyFjscu9eohGXJkiUIDQ2FUqlEZGQkdu3aVW/badOmQSKR1Dq6du1ao11hYSGeffZZBAQEQKlUIjw8HJs3b25MeGQBRoT7wcdNgbwbFUg8zcm3RGR9bq5uK5Fw/orYzE5Y1q5di1mzZmHevHk4cuQIoqOjMXr0aGRk1N31v3jxYqjVatORmZkJLy8vTJgwwdSmoqICI0eOxKVLl/Djjz/i7NmzWLZsGYKCghp/ZyQqR5kUE6sn364+wMq3RGR9jAXjBnXg/BVL4GDuGxYuXIi4uDjMmDEDALBo0SL8/vvvWLp0KRYsWFCrvUqlgkr1Vynjn3/+GQUFBZg+fbrpXEJCAq5fv47k5GQ4OjoCAEJCQm4bh1arhVarNf2+qKjI3FuhZjapXzA+35GGPWn5uJRXgrbeLmKHRETUIIWlFTiVZfheGcj9gyyCWT0sFRUVSElJQWxsbI3zsbGxSE5ObtA14uPjERMTUyMh2bBhAwYOHIhnn30Wfn5+iIiIwHvvvQedTlfvdRYsWGBKhlQqFYKDg825FWoBrT2dMTTMBwDw3QFOviUi67EvPR+CAHTwdYWfu1LscAhmJix5eXnQ6XTw8/Orcd7Pzw/Z2dl3fL9arcaWLVtMvTNG6enp+PHHH6HT6bB582a89tpr+OSTT/Duu+/We625c+dCo9GYjsxMloK3RI9X7y/0Q8oVaKvqT0CJiCzJnjTuzmxpzB4SAlBr8pEgCA2akLRixQp4eHhg/PjxNc7r9Xr4+vriyy+/hEwmQ2RkJLKysvDRRx/hP//5T53XUigUUCgUjQmfWtA9nX3h765EdlE5fj91Dff3CBQ7JCKiO2L9FctjVg+Lt7c3ZDJZrd6UnJycWr0utxIEAQkJCZgyZQrkcnmN1wICAhAWFgaZTGY6Fx4ejuzsbFRUVJgTIlkYB5kUE/tWT77dz8m3RGT5sjXluJBbAokEGNDOS+xwqJpZCYtcLkdkZCQSExNrnE9MTERUVNRt35uUlIS0tDTExcXVem3QoEFIS0uDXv9XvY5z584hICCgVnJD1mdSv2BIJcC+9Ou4kHtD7HCIiG5rb7qhdyUiUAUPZ9v8DiqrLIMgCGKHYRazlzXPmTMHy5cvR0JCAlJTUzF79mxkZGRg5syZAAxzS6ZOnVrrffHx8ejfvz8iIiJqvfbMM88gPz8fL7zwAs6dO4dNmzbhvffew7PPPtuIWyJLE6Bywj2dfQEA33F/ISKycLY+f0Wn1+Hh7x/G5J8mo6yyTOxwGszsOSwTJ05Efn4+5s+fD7VajYiICGzevNm06ketVteqyaLRaLBu3TosXry4zmsGBwdj69atmD17Nrp3746goCC88MILePnllxtxS2SJHu/fBttSc/Dj4Sv416hOUDrK7vwmIqIWJgjCXwXjbLT+yoGrB5CYnogqfRXO5J3BzxN/RrDK8lfaSgRr6xOqR1FREVQqFTQaDdzd3cUOh26h0wsY8uF2XC0sw6KJPTG+F4sCEpHluZRXgmEf74CjTIJjb8TCWd6otSkWb8elHXjk+0eQX5YPXxdfrH90PQa1GSRKLA39/uZeQtQiZFLJTZNvOSxERJYpubp3pVewp80mKwAwrO0wHHr6ELr7dUdOSQ6Gfz0cy1KWiR3WbTFhoRYzsW8wZFIJDly6jvPXisUOh4ioFmM5/ig72J25rUdbJD+VjEe6PIJKfSWe3vg0ntv8HCp1lWKHVicmLNRi/NyVGFE9+XY1K98SkYXR62+av2In9Vdc5C74/pHv8fbwtwEAnx/8HLHfxCK3JFfkyGqz3f4uskiP92+DraevYV3KFbx8b2dOvqVG+yP1Gt7ZlIrrJRVQOTnCw9kRKidHuDs5wsPJscY5wyGHyskR7Xxc+HNHdTp7rRjXSyrg5ChDz2APscNpMRKJBK8NeQ3dfLth8k+TsePSDvRd1he/TPoFPfx7iB2eCRMWalFDOvqgtacTrhSUYdNxNR6ObC12SGRlisor8c7G0/j+0BXTOU1ZJTKuN+z9QR5OWPv3AWjt6dxMEZK12pNmGA7qG+oFuYP9DUA80PkB7IvbhwfWPIALBRcQlRCFr8d/jUe6PCJ2aACYsFALk0oleKxfG3z0+1msOZjBhIXMsictDy/9eBxXC8sgkQB/i26HCZGtUVReicLSSmjK/vqv8SgsrTD9+lqRFlcLyzDj60P48ZkouCr4VyD9xTgcNMhG6680RFffrjjwtwOY9OMkJKYnYsIPE/Ba9Gt4a/hbkErETeL4fyu1uEciW+OTrWdx8FIBLuTeQHsfV7FDIgtXWlGFD7acwdd7Dds7tPFyxieP9kDftuaVTb9aWIYHPtuDM9nFeOG7I/hyah/IpHfeB41sX5VOj/0XDd109jJ/pT5eTl7Y/MRmvJz4MhbuW4h3dr2D4znHserBVXBXiFc2xP76vEh0fu5KDO9kmHz7/SHusk23l3L5Ou5bvMuUrEwZEIItL0SbnawAhuGgZVMjIXeQ4o8zOXh/S2pTh0tW6vhVDW5oq6ByckSXQNbycpA64JNRn+Dr8V9DIVNgw9kNGBg/EGnX00SLiQkLieLR6pos61KuolKnv0NrskfaKh3e33IGE/63F5fySxGgUmLlU/3w9vgIuNzFUE6vNp74eIJhIuGyXRex9iBXrBGQXD1/ZUA7L/a63WRqj6nYOX0nAt0CkZqbKmrCwiEhEsU9nX3h7SpH3g0ttp/JQWxXf7FDIgty8qoGL35/DGer6/U83Ls1/jOuC1ROjk1y/ft7BCIt5wb++8d5zPvpJEJauWBAO/udt0B/FYwbZKPl+O9Gv6B+OPS3Q9h+aTvu7XCvaHGwh4VE4SiT4qHehgm3HBYio0qdHou3ncf4z/fg7LVieLvK8cWUSHzyaI8mS1aMZo3oiDHdA1ClFzDzmxRcyitp0uuT9Siv1OHQ5QIAnL9SnwC3ADze7XFRY2DCQqJ5tI9hWGj72VzkFJWLHA2J7WphGR5emoz/23YOVXoBoyP88fusIRjVTL1vUqkEn0zogR6tVSgsrUTc1wehKbPMCp/UvFIuF6CiSg9fNwXa+7iIHQ7VgwkLiaaDrysiQzyh0wv48fCVO7+BbNbFvBJMWJqM41c0cFc6YPGknljyRG+0clU06+cqHWVYNrUPAlRKXMgtwXOrD6OKc6rsTnJ1Of5BHbwhkXD+iqViwkKimljdy/LDoSuwkY3DyUxns4sx4X97kaUpRzsfF/w2awge6BnUYl8cvu5KLJvaB06OMuw6n4e3fj3dIp9LlmNPmmH+ykA7rr9iDZiwkKjGdA+Ai1yGi3klOHCxgaVKyWacuKLBxC/3Iu+GFp393fD93wci0MOpxeOICFJh0aSekEiAVfsu4+vkSy0eA4kjp7gcx68UAuCEW0vHhIVE5aJwwNjugQCAtZx8a1cOXrqOx5ftQ2FpJXoEe2DN0wPg3cxDQLczqqs/XhrVGQDw1q+nkHTO8jZ/o6a34WgW9ALQM9gDQSIky9RwTFhIdMaaLJtPqFFczkmP9mD3+TxMjT+AYm0V+od64dsZ/eHhLBc7LMwc2g4P924NvQA89+1hpOUUix0SNbOfjlwFADzUO0jkSOhOmLCQ6Hq38UAHX1eUV+rx6zG12OFQM9t2+hqeWnEQZZU6DA3zwYrp/SxmTx+JRIL3HopA37aeKNZW4akVh3C9pELssKiZnLtWjFNZRXCQSkw9vWS5mLCQ6CQSiWnyLYeFbNuvx7Iw85sUVOj0GNXVD19OjYSTXCZ2WDUoHGT43+RIBHs5IeN6KWZ+k8JqzDbK2LsyrJMvvFzE7+Gj22PCQhbhwd5BcJBKcCyzEGez2Q1vi74/mInn1xxBlV7Ag72C8PnjvaFwsKxkxaiVqwLxT/aFm8IBBy5ex6rqfYzIduj1An7hcJBVYcJCFsHbVYGYcD8AwNqD7GWxNSv2XMRL645DEIDH+rXBJxN6wEFm2X/9hPm54ZX7DJNwF/9xHoWlHBqyJfsu5iNLUw43pQPu6ewrdjjUAJb9NwbZlYnVk29/OnIF2iqdyNFQU/l8exrerK5tMmNwKN57MAJSK9lcbmKfYHTyc4OmrBKL/zgvdjjUhH46bOhdGds9AEpHy+zpo5qYsJDFGBLmA393JQpKK7HtdI7Y4dBdEgQBH/1+Bh/9fhYA8PyIjpg3JtyqKok6yKR4bWw4AGDV3stIz70hckTUFMoqdNhyMhsAML4nh4OsBRMWshgyqQSPRBo2ROTkW+sXv/siPt9+AQAwd3RnzBkZZlXJilF0Rx8M7+SDKr2A9zafETscagKJqddwQ1uFIA8n9G3rJXY41EBMWMiiGDdE3HU+F1cLy0SOhhpLU1aJ/1YPocy7Lxx/H9pe5Ijuzrwx4ZBJJdiWeg3JaXlih0N36efqybYP9gqymuFJYsJCFqZNK2cMbNcKggD8eIgbIlqrhN0XUVRehY6+rnhqcKjY4dy1Dr5umNy/DQDg7U2p0Om575W1yruhNVUxfpCrg6wKExayOI/2NQwL/ZCSCT2/GKyOprQSCbsvAgBmxYRBZiP/gn0hJgxuSgekqovwYwqHLK3Vr8eyoNML6NFahfY+rmKHQ2ZgwkIWZ3REANyUDrhSUIbkC/lih0NmWr47HcXaKnT2d8PoCH+xw2kyXi5yvDCiIwDg463ncENbJXJE1Bg/3TQcRNaFCQtZHKWjDA/05IaI1qigpOKm3pWONjc/YOrAtmjbyhm5xVr8b8cFscMhM6Xl3MDxKxrIpBKM68FS/NaGCQtZpIl9DPMFfj+VzYJdVmTZrnSUVOjQJcAdsV1sp3fFSO4gxSujDcucl+1K58RwK/PTEcO8uKFhPmgl4s7g1DiNSliWLFmC0NBQKJVKREZGYteuXfW2nTZtGiQSSa2ja9eupjYrVqyos015eXljwiMbEBHkjvAAd1RU6U0z+smy5d/QYkXyJQDA7JFhNte7YjSqqx/6h3pBW6XHh79xmbO10OsF/HwkCwCHg6yV2QnL2rVrMWvWLMybNw9HjhxBdHQ0Ro8ejYyMjDrbL168GGq12nRkZmbCy8sLEyZMqNHO3d29Rju1Wg2lUtm4uyKrZ9gQ0ViT5QoEgZNvLd2XO9NRWqFDtyAVYsJtt9S5RCLB62O7QCIBfjmahSMZBWKHRA1w8NJ1XC0sg5vCASO7+IkdDjWC2QnLwoULERcXhxkzZiA8PByLFi1CcHAwli5dWmd7lUoFf39/03Ho0CEUFBRg+vTpNdpJJJIa7fz9ba87mcwzvlcQ5A5SpKqLcPJqkdjh0G3kFmuxsnqDwNkjO1plgThzRASp8HBvQ0L99sbTTKitgHGy7ehu/izFb6XMSlgqKiqQkpKC2NjYGudjY2ORnJzcoGvEx8cjJiYGISEhNc7fuHEDISEhaN26NcaOHYsjR47c9jparRZFRUU1DrItHs5yjOpqSFzXHqq7B48swxdJF1BWqUOPYA8M72S7vSs3+/eoTnBylOFwRiE2HleLHQ7dRnmlDptOGJ7Rg71aixwNNZZZCUteXh50Oh38/Gp2p/n5+SE7O/uO71er1diyZQtmzJhR43znzp2xYsUKbNiwAd999x2USiUGDRqE8+fr32xswYIFUKlUpiM4ONicWyErMbG68u0vR7NQXskNES1RTlE5Vu2r7l2Jsf3eFSM/dyVmVlfwfX/LGf58WrA/UnNQXG4oxd8/lKX4rVWjJt3e+heSIAgN+ktqxYoV8PDwwPjx42ucHzBgACZPnowePXogOjoa33//PcLCwvDpp5/We625c+dCo9GYjsxMLn+1RVHtW6G1pxOKy6uw5ST/FWuJliZdgLZKj95tPDA0zEfscFrU00PaIUClxNXCMiTsuSh2OFQP4+qgB3oG2uxkcHtgVsLi7e0NmUxWqzclJyenVq/LrQRBQEJCAqZMmQK5XH77oKRS9O3b97Y9LAqFAu7u7jUOsj1SqQQTIg29LD+wVL/FydaU49v9huG6OSM72U3vipGTXIaX7u0EAFiy/QJyi7UiR0S3yr+hxY6z1aX4uTrIqpmVsMjlckRGRiIxMbHG+cTERERFRd32vUlJSUhLS0NcXNwdP0cQBBw9ehQBAQHmhEc26qHq/T72pufjWhGXuluSJTvSUFGlR9+2nhjUoZXY4YjigR5B6N5ahRvaKixMPCd2OHSLTSfUqNILiAhyR0c/N7HDobtg9pDQnDlzsHz5ciQkJCA1NRWzZ89GRkYGZs6cCcAwVDN16tRa74uPj0f//v0RERFR67W33noLv//+O9LT03H06FHExcXh6NGjpmuSfQv2ckbvNh4QBGATJzdajKzCMqw5YBiKnT0yzO56V4ykUsMyZwBYezADZ7K5AMCSrD9sLMXPybbWzuyEZeLEiVi0aBHmz5+Pnj17YufOndi8ebNp1Y9ara5Vk0Wj0WDdunX19q4UFhbi6aefRnh4OGJjY3H16lXs3LkT/fr1a8QtkS26v7qM9oZjWSJHQkafb09DhU6PAe28ENXeW+xwRNW3rRfGdAuAXgDe2ZjKZc4WIj33Bo5mFkImlZj+DiHrJRFs5P+soqIiqFQqaDQazmexQTnF5Rjw3h/QC8Cul4Yj2MtZ7JDs2pWCUgz/eAcqdQLWPj0A/dvZ53DQzTLySxGzMAkVOj2+mt7XbpZ3W7KFW8/iv3+mYVgnH6yYzn8AW6qGfn9zLyGyCr5uSgxsb/hSZC+L+D7fnoZKnYBBHVoxWanWppUzpg1qCwBYtO08e1lEJggCfjrKnZltCRMWshrjuhu6dH9lwiKqjPxS04qt2TFhIkdjWZ4e0g5KRymOZRZi1/k8scOxa4cuFyDzehlcFQ42uRGnPWLCQlbj3gh/OMokOJNdjPPXisUOx259+ud5VOkFDAnzQZ+2LMJ1M29XBR7vZ5jP998/2MsiJmMp/nsj/OEkZyl+W8CEhayGh7McQzoaCpOxl0Ucl/JKsL76i2B2TEeRo7FMfx/aDnIHKQ5dLsDe9Hyxw7FL2iqdaUUhh4NsBxMWsir39/xrtRD/9dry/vvneej0AoZ38kGvNp5ih2OR/NyVmNTXUOzw0z/SRI7GPm0/kwNNWSX83ZUYwDlWNoMJC1mVmHA/KB2luJRfyh2cW9iF3Bv42di7MpJzV25n5tD2cJRJsDc9HwcvXRc7HLtjrL3yQK9AyFiK32YwYSGr4qJwwIhwwzYQG45dFTka+/L59jToBSAm3BfdW3uIHY5FC/RwwiPVW0r894/6txihppd3Q4vtZ3MAAA+xWJxNYcJCVse4WmjjcTX0eg4LtYRKnR6Jp64BgGmHYrq9fwxrD5lUgl3n83Ako0DscOzGf/84j0qdgB7BHujkz1L8toQJC1mdYZ184KZwgFpTjkOX+UXQEo5mFqJYWwVPZ0fOXWmgYC9nPFQ94fPTPzmXpSVcyL1h2ozz5epNKcl2MGEhq6N0lCG2q6GuAlcLtYyd5wy73Q7u6MM5AWZ4dngHSCXAn2dycPKqRuxwbN77W85ApxcQE+5r99tF2CImLGSVjKuFNp9Qo0qnFzka22dMWIZ05JeAOdp6u5j2sPn0T85laU770vORePoaZFIJXhndWexwqBkwYSGrFNW+Fbxc5MgvqUDyBda6aE7XSypwvLp3YEiYj8jRWJ/n7ukAiQT4/dQ1pKq5sq056PUC3tucCgCY1DcYHXw5d8UWMWEhq+Qok+K+boZhIe4t1Lx2p+VBEIDO/m7wc1eKHY7V6eDrhvu6BQAAPtvOuSzN4dfjWTh+RQMXuQyzuF2EzWLCQlbLuFro95PZ0FbpRI7GdpmGg9i70mj/vKcDAMMQZloOt5VoSuWVOnz421kAwDPD2sPHTSFyRNRcmLCQ1erb1gv+7koUa6uw42yu2OHYJEEQsOu84c82mvNXGq2zvztGdfWDIACfccVQk/o6+RKuFpbB312JuMHtxA6HmhETFrJaUqkEY7sbutq5Wqh5nL1WjGtFWigdpejLjQ7vyj/vMey9tOFYFi7mlYgcjW24XlJhGmb716hO3OTQxjFhIatmXC20LfUaSrRVIkdje4zDQf1DW0HpyC+DuxERpMI9nX2hFwxVg+nu/feP8ygur0KXAHducmgHmLCQVesWpEJIK2eUV+qxLfWa2OHYnJ3n8gBw/kpTMc5l+enIVWReLxU5Gut2Ma8E3+y7DACYNyac9YHsABMWsmoSicRU54LDQk2rrEKHA9Ub9w0N4/yVptCrjSeiO3pDpxewZMcFscOxah9sOYOq6p3DB3Xgz6c9YMJCVm9cdcKSdC4XmtJKkaOxHfsu5qOiSo9AlRLtfVzFDsdmvDDCMJflx5RMXC0sEzka63Tw0nX8diobUgkw975wscOhFsKEhaxemJ8bOvm5oVIn4LdTarHDsRk3L2eWSNjd3lT6tPXCwHatUKkT8EUSe1nMJQgC3tlkKBI3sW8bhPmxSJy9YMJCNsE4+fbXY0xYmgrrrzSff44wzGVZczAT14rKRY7Gumw8rsaxzEI4y2WYPbKj2OFQC2LCQjbBuLw5+UIecor5BXC3rhaW4UJuCaQSYBA3kWtyA9u1Qt+2nqio0uOLpHSxw7Ea2iodPvjtDABg5tD28HVj5WV7woSFbEJIKxf0CPaAXgC2nMgWOxyrt6u6d6VnsAdUzo4iR2N7JBKJqS7L6gOXkXdDK3JE1mFl8mVcKSiDn7sCM6JDxQ6HWhgTFrIZ46p7Wbi30N3beZ7DQc0tuqM3egR7oLxSz7osDVBYWmHa8frFkZ3gLHcQOSJqaUxYyGaM7R4IiQRIuVyAKwWscdFYVTo9dp9n/ZXmJpFI8K9Yw0Z9K/de5k7Od/Dpn2koKq9CZ383PBzZWuxwSARMWMhm+KuU6FddPn7jcU6+baxjVzQoKq+CyskRPVp7iB2OTYvu6IP7uvlDpxcw76cT0OsFsUOySJfzS7By7yUAwKv3sUicvWLCQjblr9VCHBZqLOPqoMEdvPnF0AL+M7YrXOQyHM4oxA8pmWKHY5E+/O0sKnUChoT5sNfPjjFhIZsyOiIADlIJTmUV4ULuDbHDsUp/zV/h6qCW4K9SYvZIw9DQgi1ncL2kQuSILMveC/nYdEINqQR49b7OYodDImLCQjbFy0WOwR0NX7TsZTGfprQSxzILAXD+SkuaFtUWnf3dUFhaife3pIodjsUoKKnA7LVHARiKxHX2dxc3IBIVExayOeO6G4aFNhzLgiBwToA5dqflQS8AHX1dEaByEjscu+Egk+LdByMAAN8fuoJD1Xs42TNBEPDvH48hu6gc7bxd8NoYluC3d41KWJYsWYLQ0FAolUpERkZi165d9badNm0aJBJJraNr1651tl+zZg0kEgnGjx/fmNCIENvVD3IHKdJzS5CqLhY7HKvC6rbiiQzxwsQ+wQCAeT+dRKVOL3JE4lqRfAnbUnMgl0nx6eO94KLgMmZ7Z3bCsnbtWsyaNQvz5s3DkSNHEB0djdGjRyMjI6PO9osXL4ZarTYdmZmZ8PLywoQJE2q1vXz5Mv71r38hOjra/DshquamdMSw6i/cTSc4LNRQgiCY5q9Ed+T8FTG8MrozPJ0dcfZaMVbsuSR2OKI5eVWDBZsNFW3njQlH10CVyBGRJTA7YVm4cCHi4uIwY8YMhIeHY9GiRQgODsbSpUvrbK9SqeDv7286Dh06hIKCAkyfPr1GO51OhyeeeAJvvfUW2rVrd8c4tFotioqKahxERmOrd3DedFzNYaEGSsu5AbWmHHIHKfqHthI7HLvk6SLH3NGGoY//23YOWXa4m/MNbRX++d0RVOj0GNnFD1MHhogdElkIsxKWiooKpKSkIDY2tsb52NhYJCcnN+ga8fHxiImJQUhIzR/C+fPnw8fHB3FxcQ26zoIFC6BSqUxHcHBww26C7MKIzr5QOkpxKb8Up7KYzDZEUvVwUP9QLzjJZSJHY78eiWyNyBBPlFbo8PbG02KH0+L+8/NJXMwrQaBKiY8e6c6dwsnErIQlLy8POp0Ofn5+Nc77+fkhO/vO+7eo1Wps2bIFM2bMqHF+z549iI+Px7Jlyxocy9y5c6HRaExHZibrF9BfXBQOuKezLwAWkWuoncbqth05f0VMUqkE74yPgEwqwZaT2dh+NkfskFrMupQrWH/kKqQSYPFjveDhLBc7JLIgjZp0e2vGKwhCg7LgFStWwMPDo8aE2uLiYkyePBnLli2Dt3fDx80VCgXc3d1rHEQ3G9PNMCy08ThXC91JeaUO+9PzAXDCrSUID3DHU4PaAgDe+OUUyit14gbUAi7k3sDrv5wEAMyOCUPf6qrVREZmJSze3t6QyWS1elNycnJq9brcShAEJCQkYMqUKZDL/8qaL1y4gEuXLmHcuHFwcHCAg4MDVq5ciQ0bNsDBwQEXLlwwJ0Qik+GdfeDkKMOVgjIcv6IROxyLduDidWir9PB3VyLMz1XscAjACzFh8HdXIuN6KZbY+OaI5ZU6/HP1EZRW6DCwXSv8Y3gHsUMiC2RWwiKXyxEZGYnExMQa5xMTExEVFXXb9yYlJSEtLa3WHJXOnTvjxIkTOHr0qOm4//77MXz4cBw9epRzU6jRnOUOGBFuHBbiaqHbMS5nju7ozTkDFsJV4YA3xnUBAPwvKR3pNly5+f0tZ3BaXQQvFzkWTerJLSGoTmYPCc2ZMwfLly9HQkICUlNTMXv2bGRkZGDmzJkADHNLpk6dWut98fHx6N+/PyIiImqcVyqViIiIqHF4eHjAzc0NERERNXpjiMw1tnsAAK4WupO/yvFzOMiS3Bvhj2GdfFCh0+P1X07a5M/w1lPZWJF8CQDwyYQe8HNXihsQWSyzE5aJEydi0aJFmD9/Pnr27ImdO3di8+bNplU/arW6Vk0WjUaDdevWNXgFEFFTGdbJFy5yGbI05ThSXXKeasrWlOPctRuQSAwbHpLlkEgkeOv+rlA4SLEnLR+/2tgE8qzCMvz7x+MAgL9Fh2J49UR5orpIBBtJ2YuKiqBSqaDRaDgBl2p4Yc0R/HI0C08NCsV/qrvY6S/fH8rESz8eR49gD/zy7CCxw6E6/PeP81iYeA4+bgr88eJQuCsdxQ7prlXp9Hhs2T4cvFSA7q1V+HFmFOQO3C3GHjX0+5s/HWTzxlbvLbT5hBp6vU3k503KOH9lKKvbWqy/D22HUG8X5BZrsXDrObHDaRL//eM8Dl4qgKvCAZ8+1ovJCt0Rf0LI5g0J84abwgHZReVIySgQOxyLotML2J1WXX+F81cslsJBhrcfMMz/W7n3Ek5ete5Vb8kX8vBp9cqn9x7qhpBWLiJHRNaACQvZPIWDDCO7Gpbdb7KxOQB368RVDQpLK+GmdEDPYA+xw6HbGNzRG+N6BEIvAM+tPozM66Vih9Qo+Te0mLXmKAQBmNgnGPdXb6NBdCdMWMgumFYLnVBDx2EhE+Nw0KD23nCQ8a8DS/f62HAEeTjhUn4pHvlfMs5mW9du5IWlFXjyqwPIKdaig68r3rifc8qo4fg3FNmFwR184K50QG6xFgcvXRc7HIthTFg4HGQdfN2UWPdMFML8XHGtSItHv9iLlMvWMcxZUFKBx5ftx8mrRWjlIsf/JveGs9xB7LDIijBhIbsgd5BiVFd/ABwWMioqrzQt9R4Sxgm31sJfpcT3fx+I3m08oCmrxOTl+7HDwvcbyr+hxWPL9uG0ugjergqseXoAOvi6iR0WWRkmLGQ3xlaPlW85qUaVTi9yNOJLTsuDTi+gnY8LWns6ix0OmcHDWY5vZvTH0DAflFXqMOPrQ/jl6FWxw6pT3g0tHl+2H2eyi+HjZkhWOvoxWSHzMWEhuxHVvhU8nR2Rd6MCBy5yWCjpHHdntmbOcgcsm9oH9/cIRJVewKy1R/F1dcVYS5FTXI7HvtyHs9eK4etm7FnhXlXUOExYyG44yqS4N8IwLGRrFUPNJQjCTfNXOBxkreQOUiya2BPTotpCEIA3NpzCwsRzFlHCP6fIkKycz7kBf3cl1v59INr7MFmhxmPCQnZlTDfDsNBvdj4sdCH3Bq4WlkHuIMWAdq3EDofuglQqwRvjumDOyDAAhoJsr/9yUtTVcNmackz6ch8u5JYgUKXE2r8PQKg3a63Q3WHCQnZlQDsvtHKRo6C0EskX8sUORzR/njFM0hzQrhVXatgAiUSC50d0xNvjIyCRAN/sy8ALa46goqrlk3K1pgyTvtyL9LwSBHk4Yc3TA1kYjpoEExayKw43DQvZ82qh7WcMw0HDO3H+ii2ZMiAEnz7WC44yCTYeVyPu64Mo0Va12OdfLSzDxC/24VJ+KVp7OmHN0wPQphUndFPTYMJCdmdMdRG5305lo9IOh4WKyitNtWiGd+LuuLZmbPdAxD/ZF85yGXadz8MTy/ejoKSi2T8383opJn6xFxnXS9HGyxlrnh6AYC8mK9R0mLCQ3ekf2grergpoyipN++jYkz3n81ClF9DO2wVtOa/AJg0J88G3M/rDw9kRRzMLMfbT3Vh7MKPZEvTM66WY9OU+XCkoQ0grQ7LCpfLU1JiwkN2RSSW4r5v9DgsZ568M78zeFVvWq40nfvj7QAR5OOFqYRleXncCMQuTsP7wlSadkHsqS4OJX+zF1cIyhHq7YO3TAxHo4dRk1ycyYsJCdmlsd8Nqod9PZUNbpRM5mpaj1wvYcc44f4UJi63r6OeGbXOGYt594fBykeNyfinmfH8MI/8vCRuOZUHfyMSlsLQCK/dewgOf7caY/+5GlqYc7XxcsObpAfBXKZv4LogMuDyA7FKfEE/4uStwrUiL3efzMCLcT+yQWsSprCLkFmvhIpehb6in2OFQC3CSy/C3Ie3weP82+HrvJXyRlI703BI8/90RfP5nGmaP7IhRXf0hkUhue50qnR47z+fix5Qr2HY6BxXVw0sOUglGhPvi7Qci4OvOZIWaDxMWsktSqQT3dQvAV3suYeNxtd0kLNur95wZ1MEbCgeZyNFQS3JROOAfwzpgyoAQJOy+hOW703H2WjFmfnMYXQPdMTsmDCPCfWslLmezi7Hu8BWsP3wVeTe0pvPhAe54JLI1HugZCG9XRUvfDtkhJixkt8Z2NyQsiaevobxSB6Wj7X+BG+ev3MP5K3bLTemIF2I6YlpUWyzfnY6E3RdxKqsIM1YeQo9gD8wZGYbuQSpsOJaFdYev4PgVjem9Xi5yjO8ZhIcjg9A1UCXiXZA9YsJCdqtXsCcCVUpkacqRdC7XtJuzrcq/ocWxK4UAgGGcv2L3VM6OeDG2E6YPCsWXO9PxdfIlHMssxJMJByCVAMbpLQ5SCe7p7ItHIltjWCdfyB049ZHEwYSF7JZxWGj57ovYdFxt8wlL0rlcCIKhK58TI8nIy0WOV0Z3RtzgUPwv6QK+2XcZ2io9utw05NOKQz5kAZiwkF0b2yMQy3dfxLZU2x8W2n7WsDrons6sbku1+bgp8PrYLnhueAcUl1exQi1ZHPbtkV3r0VqF1p5OKK3QYXv1/A5bVKXTm3Zn5nJmuh1PFzmTFbJITFjIrkkkElOp/o0nbLeI3JHMQmjKKuHh7IhebbicmYisDxMWsntjuxmKyP2ZmtOiG8W1JGPv0ZCOPpBJb19vg4jIEjFhIbsXEeSOUG8XlFXqsPF4ltjhNAsuZyYia8eEheyeRCLB4/3aAABW7r0MQWi6fVYsgVpThjPZxZBIDJviERFZIyYsRAAeiWwNhYMUp7KKcDSzUOxwmtT2M4bJtr2CPeDlIhc5GiKixmHCQgTDygjjhoir9l0WOZqmZSzHz9VBRGTNmLAQVZsyMAQAsPG4GtdLKkSOpmloq3TYk5YHABjO+StEZMUalbAsWbIEoaGhUCqViIyMxK5du+ptO23aNEgkklpH165dTW3Wr1+PPn36wMPDAy4uLujZsydWrVrVmNCIGq1HaxW6BalQUaXHD4cyxQ6nSRy4eB2lFTr4uinQNdBd7HCIiBrN7IRl7dq1mDVrFubNm4cjR44gOjoao0ePRkZGRp3tFy9eDLVabToyMzPh5eWFCRMmmNp4eXlh3rx52Lt3L44fP47p06dj+vTp+P333xt/Z0RmkkgkmDLA0Mvyzf7L0Outf/KtcXXQsE4+tXbhJSKyJmYnLAsXLkRcXBxmzJiB8PBwLFq0CMHBwVi6dGmd7VUqFfz9/U3HoUOHUFBQgOnTp5vaDBs2DA8++CDCw8PRvn17vPDCC+jevTt2797d+DsjaoRxPQLhrnRA5vUyJJ3PFTucu7bDVI6fw0FEZN3MSlgqKiqQkpKC2NjYGudjY2ORnJzcoGvEx8cjJiYGISEhdb4uCAL++OMPnD17FkOGDKn3OlqtFkVFRTUOorvlJJdhQp9gAMA3e6178u3FvBJczCuBo0yCQR28xQ6HiOiumJWw5OXlQafTwc/Pr8Z5Pz8/ZGdn3/H9arUaW7ZswYwZM2q9ptFo4OrqCrlcjjFjxuDTTz/FyJEj673WggULoFKpTEdwcLA5t0JUryf6G2qy/Hk2B5nXS0WOpvGM1W37tvWCm9JR5GiIiO5Ooybd3joWLghCg8bHV6xYAQ8PD4wfP77Wa25ubjh69CgOHjyId999F3PmzMGOHTvqvdbcuXOh0WhMR2ambUySJPG183FFdEdvCAKw+kDdc7OsAZczE5EtcTCnsbe3N2QyWa3elJycnFq9LrcSBAEJCQmYMmUK5PLaxaukUik6dOgAAOjZsydSU1OxYMECDBs2rM7rKRQKKBQKc8InarDJA0Kw63we1h7MxKyYjlA4yMQOySwl2irsT78OgMuZicg2mNXDIpfLERkZicTExBrnExMTERUVddv3JiUlIS0tDXFxcQ36LEEQoNVqzQmPqMmM6OyLAJUS10sqsOXEnYc7Lc2etDxU6PQI9nJCex8XscMhIrprZg8JzZkzB8uXL0dCQgJSU1Mxe/ZsZGRkYObMmQAMQzVTp06t9b74+Hj0798fERERtV5bsGABEhMTkZ6ejjNnzmDhwoVYuXIlJk+e3IhbIrp7DjKpaX8ha6x8u924OqiTL5czE5FNMGtICAAmTpyI/Px8zJ8/H2q1GhEREdi8ebNp1Y9ara5Vk0Wj0WDdunVYvHhxndcsKSnBP/7xD1y5cgVOTk7o3LkzvvnmG0ycOLERt0TUNCb2C8biP84j5XIBTmVp0DVQJXZIDSIIAnZUz18ZxuEgIrIREsFGtqYtKiqCSqWCRqOBuzsrelLTeHb1YWw6rsZj/dpgwUPdxA6nQVLVRRi9eBeUjlIc/U8slI7WNf+GiOxLQ7+/uZcQ0W0YK9/+fOQqisorRY6mYYyrg6LaezNZISKbwYSF6Db6h3qho68ryip1+OnwVbHDaRBj/RWuDiIiW8KEheg2JBKJaRfnVfsuw9JHUDWllUi5XAAAGBbmI3I0RERNhwkL0R082CsIznIZ0nJuYF91bRNLlXQ+F3oB6OjrimAvZ7HDISJqMkxYiO7ATemIB3sFAQC+sfAlzjuqh4O42SER2RomLEQNMLl68u3vp7Jxrahc5GjqptML2HHOUH9lGMvxE5GNYcJC1ADhAe7o29YTVXoBaw5Y5r5Vx68U4npJBdwUDujT1lPscIiImhQTFqIGMvayrD5wGZU6vcjR1GZcHRQd5g1HGf/XJiLbwr/ViBro3gh/eLvKca1Iiz9Sr4kdTi3GcvzcnZmIbBETFqIGUjjIMLFvMADL218o83opTlzVAACGduJyZiKyPUxYiMzwWL82kEqAPWn5SMu5IXY4Jh/+fhYAMLiDN3zdlCJHQ0TU9JiwEJmhtacz7unsBwD4dr9l9LKkXC7Ar8eyIJEAc+/rLHY4RETNggkLkZkmD2gDAPgx5QpKK6pEjUWvF/D2xtMAgAmRra1mR2kiInMxYSEy05COPmjj5Yzi8ip8kZQuaiy/Hs/C0cxCOMtl+FdsJ1FjISJqTkxYiMwklUrwYmwYAOCz7Wk4klEgShzllTp8sOUMAOAfw9rD151zV4jIdjFhIWqEB3oGYVyPQOj0AuZ8f0yUoaHlu9KRpSlHoEqJGdHtWvzziYhaEhMWokZ654EIBKiUuJhXgnc2pbboZ+cUl2PJjgsAgJdHd4bSUdain09E1NKYsBA1ksrZER9P6AEAWL0/o0WLyX3y+zmUVujQM9gD9/cIbLHPJSISCxMWorswqIM34gaHAgBeXncceTe0zf6Zp7I0+D7FsJ/R62PDIZFImv0ziYjExoSF6C79e1QndPJzQ96NCryy7gQEQWi2zxIEAe9uSoUgAGO7ByAyxKvZPouIyJIwYSG6S0pHGf5vYk/IZVJsS72GtQebbzfnbak5SL6QD7mDFC/fyyJxRGQ/mLAQNYEuge741yjDUuf5G0/jUl5Jk39GRZUe7202TO6NGxyKYC/nJv8MIiJLxYSFqInMGNwOA9p5obRCh9nfH0WVTt+k1/9m32VczCuBt6sc/xjWvkmvTURk6ZiwEDURqVSCTx7tCTelA45kFJqWHTeFwtIKLP7jPABgzshOcFM6Ntm1iYisARMWoiYU5OGEtx+IAAAs/uM8jmUWNsl1F/9xHpqySnT2d8PEvsFNck0iImvChIWoiT3QMxBjuwdApxcwe+3Ru66CeyH3BlbtNewM/dqYLpBJuYyZiOwPExaiJiaRSPDu+G7wd1ciPa/ENFG2sRZsTkWVXsA9nX0xuKN3E0VJRGRdmLAQNQOVsyM+edRQBfebfRnYfianUdfZk5aHbak5cJBK8Op94U0ZIhGRVWHCQtRMBnXwxlODDFVw//3jceSbWQVXpxfw9sbTAIDJA0LQwde1yWMkIrIWTFiImtFL93ZCmJ8r8m5oMXe9eVVwfziUiTPZxXBXOuCFER2bMUoiIsvnIHYARLZM6SjDoom98MDnu7H19DVEvPE7lI4yKBykhv+afm34vdJBBoWjFEoHGbZVb6b4/IiO8HSRi3wnRETialQPy5IlSxAaGgqlUonIyEjs2rWr3rbTpk2DRCKpdXTt2tXUZtmyZYiOjoanpyc8PT0RExODAwcONCY0IovTJdAd/xnbBQ5SCUoqdMgvqUCWphzpeSVIVRfhaGYh9qVfx46zufjtVDZ+OZqFtYcykV9SgbatnDF1YFuxb4GISHQSwcyd2tauXYspU6ZgyZIlGDRoEL744gssX74cp0+fRps2bWq112g0KCsrM/2+qqoKPXr0wD//+U+8+eabAIAnnngCgwYNQlRUFJRKJT788EOsX78ep06dQlBQUIPiKioqgkqlgkajgbu7uzm3RNQiCksrUFhaCW2VHuWVOsNRpYe2+r/llTpojb+v1KFSJ2BcjwB08HUTO3QiombT0O9vsxOW/v37o3fv3li6dKnpXHh4OMaPH48FCxbc8f0///wzHnroIVy8eBEhISF1ttHpdPD09MRnn32GqVOn1tlGq9VCq/1rEmNRURGCg4OZsBAREVmRhiYsZg0JVVRUICUlBbGxsTXOx8bGIjk5uUHXiI+PR0xMTL3JCgCUlpaisrISXl5e9bZZsGABVCqV6QgOZvVPIiIiW2VWwpKXlwedTgc/P78a5/38/JCdnX3H96vVamzZsgUzZsy4bbtXXnkFQUFBiImJqbfN3LlzodFoTEdmZmbDboKIiIisTqNWCUkkNUuDC4JQ61xdVqxYAQ8PD4wfP77eNh9++CG+++477NixA0qlst52CoUCCoWiwTETERGR9TIrYfH29oZMJqvVm5KTk1Or1+VWgiAgISEBU6ZMgVxe9xLNjz/+GO+99x62bduG7t27mxMaERER2TCzhoTkcjkiIyORmJhY43xiYiKioqJu+96kpCSkpaUhLi6uztc/+ugjvP322/jtt9/Qp08fc8IiIiIiG2f2kNCcOXMwZcoU9OnTBwMHDsSXX36JjIwMzJw5E4BhbsnVq1excuXKGu+Lj49H//79ERERUeuaH374IV5//XWsXr0abdu2NfXguLq6wtWV5ciJiIjsndkJy8SJE5Gfn4/58+dDrVYjIiICmzdvNq36UavVyMjIqPEejUaDdevWYfHixXVec8mSJaioqMAjjzxS4/wbb7xhqtVCRERE9svsOiyWioXjiIiIrE+z1GEhIiIiEgMTFiIiIrJ4TFiIiIjI4jFhISIiIovHhIWIiIgsXqNK81si42KnoqIikSMhIiKihjJ+b99p0bLNJCzFxcUAwF2biYiIrFBxcTFUKlW9r9tMHRa9Xo+srCy4ubk1aCPGhioqKkJwcDAyMzPtpr4L75n3bKt4z7xnW2XN9ywIAoqLixEYGAiptP6ZKjbTwyKVStG6detmu767u7vV/RDcLd6zfeA92wfes32w1nu+Xc+KESfdEhERkcVjwkJEREQWjwnLHSgUCrzxxhtQKBRih9JieM/2gfdsH3jP9sEe7tlmJt0SERGR7WIPCxEREVk8JixERERk8ZiwEBERkcVjwkJEREQWjwkLERERWTwmLHewZMkShIaGQqlUIjIyErt27RI7pGbz5ptvQiKR1Dj8/f3FDqtJ7dy5E+PGjUNgYCAkEgl+/vnnGq8LgoA333wTgYGBcHJywrBhw3Dq1Clxgm0id7rnadOm1XruAwYMECfYJrBgwQL07dsXbm5u8PX1xfjx43H27NkabWztOTfknm3tOS9duhTdu3c3VXYdOHAgtmzZYnrd1p4xcOd7trVnfCsmLLexdu1azJo1C/PmzcORI0cQHR2N0aNHIyMjQ+zQmk3Xrl2hVqtNx4kTJ8QOqUmVlJSgR48e+Oyzz+p8/cMPP8TChQvx2Wef4eDBg/D398fIkSNNm2taozvdMwDce++9NZ775s2bWzDCppWUlIRnn30W+/btQ2JiIqqqqhAbG4uSkhJTG1t7zg25Z8C2nnPr1q3x/vvv49ChQzh06BDuuecePPDAA6akxNaeMXDnewZs6xnXIlC9+vXrJ8ycObPGuc6dOwuvvPKKSBE1rzfeeEPo0aOH2GG0GADCTz/9ZPq9Xq8X/P39hffff990rry8XFCpVML//vc/ESJserfesyAIwpNPPik88MADosTTEnJycgQAQlJSkiAI9vGcb71nQbD95ywIguDp6SksX77cLp6xkfGeBcH2nzF7WOpRUVGBlJQUxMbG1jgfGxuL5ORkkaJqfufPn0dgYCBCQ0MxadIkpKenix1Si7l48SKys7NrPHOFQoGhQ4fa9DMHgB07dsDX1xdhYWH429/+hpycHLFDajIajQYA4OXlBcA+nvOt92xkq89Zp9NhzZo1KCkpwcCBA+3iGd96z0a2+owBG9qtuanl5eVBp9PBz8+vxnk/Pz9kZ2eLFFXz6t+/P1auXImwsDBcu3YN77zzDqKionDq1Cm0atVK7PCanfG51vXML1++LEZILWL06NGYMGECQkJCcPHiRbz++uu45557kJKSYvVlvgVBwJw5czB48GBEREQAsP3nXNc9A7b5nE+cOIGBAweivLwcrq6u+Omnn9ClSxdTUmKLz7i+ewZs8xnfjAnLHUgkkhq/FwSh1jlbMXr0aNOvu3XrhoEDB6J9+/b4+uuvMWfOHBEja1n29MwBYOLEiaZfR0REoE+fPggJCcGmTZvw0EMPiRjZ3Xvuuedw/Phx7N69u9Zrtvqc67tnW3zOnTp1wtGjR1FYWIh169bhySefRFJSkul1W3zG9d1zly5dbPIZ34xDQvXw9vaGTCar1ZuSk5NTK2u3VS4uLujWrRvOnz8vdigtwrgiyp6fOQAEBAQgJCTE6p/7P//5T2zYsAHbt29H69atTedt+TnXd891sYXnLJfL0aFDB/Tp0wcLFixAjx49sHjxYpt+xvXdc11s4RnfjAlLPeRyOSIjI5GYmFjjfGJiIqKiokSKqmVptVqkpqYiICBA7FBaRGhoKPz9/Ws884qKCiQlJdnNMweA/Px8ZGZmWu1zFwQBzz33HNavX48///wToaGhNV63xed8p3uui7U/57oIggCtVmuTz7g+xnuui809Y7Fm+1qDNWvWCI6OjkJ8fLxw+vRpYdasWYKLi4tw6dIlsUNrFi+++KKwY8cOIT09Xdi3b58wduxYwc3Nzabut7i4WDhy5Ihw5MgRAYCwcOFC4ciRI8Lly5cFQRCE999/X1CpVML69euFEydOCI899pgQEBAgFBUViRx5493unouLi4UXX3xRSE5OFi5evChs375dGDhwoBAUFGS19/zMM88IKpVK2LFjh6BWq01HaWmpqY2tPec73bMtPue5c+cKO3fuFC5evCgcP35cePXVVwWpVCps3bpVEATbe8aCcPt7tsVnfCsmLHfw+eefCyEhIYJcLhd69+5dY5mgrZk4caIQEBAgODo6CoGBgcJDDz0knDp1SuywmtT27dsFALWOJ598UhAEw5LXN954Q/D39xcUCoUwZMgQ4cSJE+IGfZdud8+lpaVCbGys4OPjIzg6Ogpt2rQRnnzySSEjI0PssButrnsFIHz11VemNrb2nO90z7b4nJ966inT380+Pj7CiBEjTMmKINjeMxaE29+zLT7jW0kEQRBarj+HiIiIyHycw0JEREQWjwkLERERWTwmLERERGTxmLAQERGRxWPCQkRERBaPCQsRERFZPCYsREREZPGYsBAREZHFY8JCREREFo8JCxEREVk8JixERERk8f4f7K026k0PQ5kAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot\n", + "v = 2 # index of variable\n", + "\n", + "batch_id = np.random.choice(len(dls.test.dataset))\n", + "#batch_id =0\n", + "var_names = list(dls.valid.dataset.dataset.data_vars.keys())\n", + "\n", + "gt = dls.test.dataset.batch_gen[batch_id].get(var_names[v]).isel(latitude=0, longitude=0).values[:36]\n", + "gtp = np.hstack([\n", + " np.ones(36)*np.nan, \n", + " dls.test.dataset.batch_gen[batch_id].get(var_names[v]).isel(latitude=0, longitude=0).values[36:]\n", + "])\n", + "targ_2 = np.hstack([np.ones(36)*np.nan, u_targ[batch_id,:,v]])\n", + "p = np.hstack([np.ones(36)*np.nan, u_pred[batch_id,:,v]])\n", + "\n", + "plt.plot(gt)\n", + "plt.plot(gtp, color='red', label='gt')\n", + "#plt.plot(targ_2, color='yellow', label='gt')\n", + "plt.plot(p, color='green')\n", + "plt.title(f\"{var_names[v]}, {batch_id}\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5447e15b-dd13-4ac3-875e-eb470a6f5ec8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "PatchTST", + "language": "python", + "name": "patchtst" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}