{ "cells": [ { "cell_type": "markdown", "id": "94a4daa2", "metadata": {}, "source": [ "# Example 4: Generate Synthetic Data for a Single Tracer\n", "In this example, synthetic observation data is generated for a single tracer. This is done by first generating a time series of tracer input. This data is then used as a model input, together with a pre-defined set of reference parameter values. As a result, we obtain a simulated time series of tracer concentrations at a (hypothetical) sampling location. From this time series we select a number of values as observations; those values are perturbed by noise to include measurement error." ] }, { "cell_type": "code", "execution_count": 1, "id": "b25444c9", "metadata": {}, "outputs": [], "source": [ "import PyTracerLab.model as ism\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from datetime import datetime" ] }, { "cell_type": "markdown", "id": "ac5cff6c", "metadata": {}, "source": [ "## 1. Generate Synthetic Observation Data" ] }, { "cell_type": "code", "execution_count": 2, "id": "47feadbf", "metadata": {}, "outputs": [], "source": [ "# make reproducible\n", "np.random.seed(42)\n", "\n", "### define input series\n", "# define time steps\n", "n_years = 50\n", "timesteps = np.arange(0.0, n_years * 12.0, 1.0) # 50 years of monthly values\n", "\n", "# this represents the tracer input to the aquifer\n", "# we start with a constant input of 1.0\n", "input_series = np.ones(len(timesteps))\n", "# we add a pulse of higher tracer input during some years, similar to the\n", "# increase in tritium in the atmosphere\n", "# this is just a bell-shaped pulse\n", "offset = 120.\n", "scale = - 0.001 # negative; closer to zero means wider pulse\n", "input_series += 2 * np.exp(scale * ((timesteps - offset) ** 2))\n", "# add some noise to the data\n", "input_series += np.random.normal(0.0, 0.05 * input_series, len(input_series))" ] }, { "cell_type": "code", "execution_count": 3, "id": "37eab1a7", "metadata": {}, "outputs": [], "source": [ "### define model (the true system; in practice we don't know this)\n", "# get decay constant\n", "# we assume a half life of 12.3 years\n", "t_half = 12.3 * 12.0\n", "lambda_ = np.log(2.0) / t_half\n", "\n", "# create true observations using the model\n", "# time step is 1 month\n", "m = ism.Model(\n", " dt=1.0,\n", " lambda_=lambda_,\n", " input_series=input_series,\n", " steady_state_input=1.,\n", " n_warmup_half_lives=10\n", ")\n", "\n", "# add an exponential-piston-flow unit\n", "# define the true model parameters\n", "epm_mtt_true = 12 * 40 # 40 years\n", "epm_eta_true = 1.5\n", "m.add_unit(\n", " ism.EPMUnit(mtt=epm_mtt_true, eta=epm_eta_true),\n", " fraction=.8, # 80 percent of the overall response\n", " prefix=\"epm\"\n", ")\n", "\n", "# add a piston-flow unit\n", "# define the true model parameters\n", "pm_mtt_true = 12 * 5 # 5 years (faster than EPM)\n", "m.add_unit(\n", " ism.PMUnit(mtt=pm_mtt_true),\n", " fraction=.2, # 20 percent of the overall response\n", " prefix=\"pm\"\n", ")\n", "\n", "true_params = m.params\n", "\n", "# simulate\n", "output_series = m.simulate()" ] }, { "cell_type": "code", "execution_count": 4, "id": "6adb4ffc", "metadata": {}, "outputs": [], "source": [ "# make reproducible\n", "np.random.seed(1234567)\n", "\n", "# define observations\n", "# randomly select a number of observations from the output series\n", "n_obs = 15\n", "obs_idx = np.random.choice(len(timesteps), n_obs)\n", "# select the corresponding timesteps and values\n", "obs_timesteps = timesteps[obs_idx]\n", "obs_values = output_series[obs_idx]\n", "# add some noise to the observations (observation error)\n", "obs_error = 0.01\n", "obs_values += np.random.normal(0.0, obs_error, n_obs)\n", "\n", "# make series we can later use in the model (has to be the same length as\n", "# the input series, filled with NaN-values where we do not have any\n", "# observations)\n", "obs_series = np.full(len(input_series), np.nan)\n", "obs_series[obs_idx] = obs_values" ] }, { "cell_type": "code", "execution_count": 5, "id": "a9ad8db8", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGwCAYAAACgi8/jAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAiodJREFUeJztnQecE0Ubxt+j997h6CBFQFBEivQiWAAVBKUoVqSpoAh+gooFEVERBAvYkKICiiIgXUEQkF6k9957J9/vGZww2dvkkrv0PH9+R5JNsplsdmeeedvEORwOhxBCCCGERCEpQt0AQgghhJBAQaFDCCGEkKiFQocQQgghUQuFDiGEEEKiFgodQgghhEQtFDqEEEIIiVoodAghhBAStaSSGOfatWuyb98+yZw5s8TFxYW6OYQQQgjxApQBPH36tBQoUEBSpHBvt4l5oQOREx8fH+pmEEIIISQJ7N69WwoVKuT2+ZgXOrDk6AOVJUuWUDeHEEIIIV5w6tQpZajQ47g7Yl7oaHcVRA6FDiGEEBJZJBZ2wmBkQgghhEQtFDqEEEIIiVoodAghhBAStcR8jA4hhJDAcvXqVbl8+XKom0EijNSpU0vKlCmTvR8KHUIIIQGrc3LgwAE5ceJEqJtCIpRs2bJJvnz5klXnjkKHEEJIQNAiJ0+ePJIhQwYWZSU+ieRz587JoUOH1OP8+fNLUqHQIYQQEhB3lRY5OXPmDHVzSASSPn16dQuxg/MoqW4sBiMTQgjxOzomB5YcQpKKPn+SE+NFoUMIISRg0F1FQn3+UOgQQgghJGqh0CGEEEJI1EKhQwghhJCohUKHRF1KIguTEUKSExPi6e+1116TcAbtu+WWW4L+ufPmzVPHJxxrJjG9nEQMW7dulX379kmtWrXcBqhNnDhR/v33X+nevTtXoyeE+Mz+/fud9ydMmCD9+vWTjRs3OrdlypTJZWKFNPpUqYI/lF66dEnSpEkT9M+NRGjRIRHDmDFjZM6cOS6djpV169apjmfFihVBbRshJHEgDDBAh+IPn+0NqMKr/7JmzaomVfoxJlGZM2eWadOmya233ipp06aVBQsWqElY8+bNJW/evEoIVa1aVWbNmuWy34sXL0rv3r0lPj5eva9kyZIyatQo5/Nr166Vpk2bqvdjP+3bt5cjR444n69bt6507dpVnnvuOcmVK5c0adLEq+/z6KOPSosWLWTw4MGq6B5qGnXp0sXF8l20aFEZMGCAtG3bVjJmzCgFCxaU4cOHO5/fsWOHOg4rV650boPlBttgycHz9erVU9uzZ8+utuNzwwVadEjEcezYsQTbrJ0YZzqEhB8YXN95552QfHafPn381i+8/PLLSjgUL15cDey7d++WZs2ayVtvvaVEzDfffCP33nuvmpQVLlxYvadDhw6yaNEiGTp0qFSqVEm2b9/uFDIQDfXr15cnnnhCPvjgAzl//rwSRa1bt1aTO83XX38tnTt3loULF/rU3rlz5yqRg9stW7bIQw89pNxbTz75pPM17733nvTt21def/11mTFjhvTo0UNKly4tjRo1SnT/EG+wpj/wwAPqO8Oarov9hQMUOiTiuXbtmnz22WdqATgNhQ4hJFC88cYbLgIgR44cSrxoYB2ZPHmyTJkyRVlhNm3aJN9//73MnDlTGjZsqF4DkaQZNmyYVK5cWd5++23nttGjRysBgfdCcIBSpUrJoEGDfG5v9uzZ1WegsnCZMmXk7rvvltmzZ7sInZo1ayoBB/B5EFMQXd4IHewXxwCggjHWpwonKHRIRGBndsZsCBcYzNIHDx50eS5FCnplCQk3MBmBZSVUn+0vbrvtNpfHZ86cUUHAU6dOVTE+V65cUVaZXbt2qefh8kFfVadOHdv9rVq1SllbzPgfDdxiWujAXZYUypcv77J8Aqw7a9ascXlN9erVEzz+8MMPJRqg0CERgelP1n5+7UPu1KlTgtejoyGEhBeI3YgGayviWEx69eqlrDVwZyH2Bm6bBx98UPVTIDE3DoQSXF3vvvtugufMxSytn5tUkRcXF6cs4d6iJ47mhDOSslspdEhEoDsMfbGhY9DARGyFQocQEizg5kHwbcuWLdVj9E8I0NVUqFBBCYv58+c7XVcmVapUUTEuCAoORQYXWLx4sVgfly1bVt3PnTu3uoW1Ci42YAYmAy1gkQwSbtC+TyICc/aA+6bwOX78uMfXE0JIIEHszKRJk9TgDzfUww8/7GIxgYDp2LGjsj7/9NNPKhAZ2UqI2wHIgkKSBbKeli5dqtxVCAh+7LHHgiYcFi5cqOJ/EBMEa/kPP/ygApK1ReqOO+6QgQMHyoYNG5Rg+9///ufy/iJFiihL0a+//iqHDx92mYyGGgodEhGYwgb3ExMytOgQQoLFkCFDVMBvjRo1lAsKqd+w0piMGDFCubOeffZZFRCMQOCzZ8+q5woUKKCEBkRN48aNlQUIaeQI6g1WvGHPnj1l2bJlymLz5ptvqu9kprAjOBr9KuKE0Da8xgQp6cjYQkAz0uMRhB0uxDm8LS4QpZw6dUrVSjh58iQLzIUxSN/EhQbQgZQrV07V1XEHAunQYRBCQsOFCxeU5aJYsWKSLl26UDeHeAAWJ4gX/EXSeeTt+E2LDol415UdtOgQQggBFDokIl1XiQkdxugQQggBzLoiEYEpbFB5E+ZMgAwFO+tNOEb+E0JIOLLDyBCLRmLWooOocsR5YE0SEv5YLTg7d+5Ut+78srToEEIIiWmhg3S+9evXq1Q+Ev64Ey7uhA5jdAghhMS00CGRhbuYHAodQgghnqDQIWEJim2ZlQ98FTp0XRFCCAEUOiTsgDUGK+2OHTs2UeGSIUMG2+179+6VTz/9VPbt2xewdhJCCAl/KHRI2LFnzx61rMOWLVucZdRRHh2gOJQZQO5pgcADBw6osuyEEEJiFwodEtacP39erZmybds29bh9+/aqzLq7VXmt6BLrhBDiDVivydPfa6+9JuEOJoaocoz1p9KkSaOWmMA6W7t27fJ5X/jOWJ8rUBWZP/zwQwk0FDok7DDjcSBUUP4b8Tr58+eXnDlzSqZMmVzq5WDNGHdgMTpCCPEWrNCt/zAIIw7Q3NarVy/na9EvhSrxwV3cIkQOFuCcNWuWjBw5UlnGx48fr25hDdeTxliCQoeEHRcvXnQROrDqACxwp4sEmq/FbKVSpUq2+zp9+rQzqBm369atkxMnTgT4GxBCIpV8+fI5/+Aqh0VDP/73338lc+bMMm3aNLW4Zdq0aWXBggVqtfHmzZurxSwxEYOggNAwQV/Vu3dviY+PV+8rWbKkjBo1yvn82rVrpWnTpur92A+s10eOHHE+X7duXbVQJiw1uXLlcllw0+SVV15RsYn4fOyvcOHCUrt2bbUaOizgKK3iyaJyyy23OK1WeB60bNlSHQf9GM/jdYiDxPdBrGTr1q3VmlNme61rZ7Vo0UIeffRR5/Ooh/b88887rWWBgkKHhB266rEWOnrmYsbjYIE3gKKPniohY7b1yy+/qLgfdCQ//vijfPTRRwH+BoQQOzDZOHvpbEj+/Ll+NVboHjhwoGzYsEEqVqyo3OvNmjWT2bNny4oVK+Suu+5Sq5ibrqIOHTrIuHHjZOjQoep9EAnaOo3JV/369dXK4VhBfPr06XLw4EElHky+/vpr1Q9ipXNYa6wgphHWm0ceeUQJM6t1+9lnn1WCR8c8JoauM/fll18qa5ZZdw4Wou+//171r2gvvjf27y2InyxUqJC88cYbTmtZoOASECSshc7EiROd902h065dOyWA9Gq2mCFpcB+zJ1zYsAbhAsTfTTfdFLTvQAhJyLnL5yTTOzdcz8HkTJ8zkjFNRr/sC4Nzo0aNnI9z5MjhYlUeMGCATJ48WaZMmaKsMJs2bVKiYObMmdKwYUP1muLFiztfjyxTiJy3337buW306NHKWoL3li5dWm0rVaqUDBo0yG27Dh8+rERT2bJlbZ8vW7asEnwQKbfffnui3zN37txOa7pVOKGf/uabb6RgwYLq8ccffyx33323vP/++wleaweOWcqUKZWFzJvXJwcKHRLWQsfEFDopUqRwihxtBj106JBUqVJFmVdhBoVI2r17t8saWYQQklxuu+02l8ew6MCdM3XqVGWZgCUZkyxt0Vm5cqUa1OvUqWO7v1WrVsncuXNd4g81cItpoQN3mTf403rlDrjEtMgB1atXVxYl9LOBFi6+QqFDIkbomFYbK+ggkFXgTY0dQkhoyJA6g7KshOqz/UXGjK6WIQQow1ozePBgFXsDa/KDDz7odLsnlhQBoQRX17vvvpvgOSRhuPtcOwsMrC9wjdmxYcMGNQlEG/WE0SqK/FVsNZD79hUKHRLWwcgmnmrm2EGhQ0h4gUHWX+6jcAIxMwiyRdCuFi7miuAVKlRQ1o758+c7XVcmsETDAg1rtJlskRRxgbie7777TrnXTMvK+fPn5ZNPPlFBzHAbaWFkxsacOnVKZbmaIIDZLgYS1ioEPSMZBCxevFh9vg4RsO4b+0CcZL169Vz6dHfxlf6EwcgkIl1XyRU6wTDtEkJiA8TOILgWLiq4oR5++GFnsVMAAdOxY0dldUZNGoiJefPmqbgdgEwoBAi3bdtWBfzCXYWg4ccee8xnIYA4HwgcxBAhO2z37t3yxx9/KIEDi8rw4cOdr0UA9Lfffit//vmnrFmzRrURLjYTtB1B1ijAikKuGoQO4PX4vnh/9+7dlcjS4gr7hisPf8hW69y5c4KMV+wbbUMlezPDzN9Q6JCwFTpt2rRRAXpJFTqezLzBmEUQQmKDIUOGSPbs2VUxU7igICpgpTEZMWKEcmchM6lMmTKq/pcuaAqrCKxC6JcaN26sLEBIzYYbClYSX0CtMVhXYDl5+umnpUSJEkqA4BYiygyC7tOnj4obuueee1QgMdK/8ToTBBfDLYfAaLM/hvvr/vvvV9lmaDOyz2Ax0kDUQQgh2wyfgc81rTkAVidYvvCZOvA5EMQ5YnxqC1MdaiUg/9/dApEkuCAD4ejRo8oUjMC2RYsWqe1ImdS+ZW/ATMNdRc+XXnqJxQQJCfCEBZYLlIIwEwdI5PPaa6+pvhUWrFCeR96O37TokLBDz3IgRMwT25+uK65uTgghsQGFDgkrkKWgXVdQ6KbVxZ+uKwodQgiJDSh0SFiBJRu0qEE6eXIsOp7M5RQ6hBCSdNdVMNxW/oJCh4QV8Llqaw5SUa0Vj33Bk8/W3YJ4hBBCogvW0SFhK3SAmeroq0UH9SgQdIy6PEh/XL58uUsNCDyva0AQQgiJTmjRIWEtdMziWUkppIUYH6RoIu3TBHUhPv/8c1WqnRBCSPRCiw4JS6GDhd4AVreF1QU1KuDKSiqo7mkHrD3JqURKCCEkvGEPT8LedYXCWskFwgnVSzdv3uyynUHJhBAS3dB1RcJS6KAIlD+BNQhl2W+++WaX7QxKJoSQ6IZCh4S1RcffWN1UtOgQQpIK1mr68MMPJVqYN2+emhRa16SKdCh0SNiAwOBz584FVehAWGEF3hhfCYUQYgGLYWK9JsQIIuOzSJEi0qNHD7U8TTRQt25dtZ6WCZI2sOK4vy3qoYZCh4RdsUCIkUCtjYMZmJmyjtWDkX2FTo0QEmacPCmyZ4/9c9iO5wPAtm3b5LbbblMxfePGjZMtW7bIyJEjVbZm9erV1UrjoQCLfpqrovubNGnSqNXHk5P4EY5Q6JCwYdmyZc7A4UBdaOXLl5e+ffsqwWOCReMIIWEERMxdd4nUqQPziutzeIzteD4AYqdLly5q0P/999/VytuFCxeWpk2byqxZs2Tv3r3yyiuvuEzQ2rZtq5acKViwoAwfPtz5HCzFqCKM96PgKaxD3bt3d8n67NWrl3of3l+tWjXlPtJ89dVXqjzGlClTpFy5cmofX3zxhZoIWt1LsDbVr19f3YfVCW3CfjNkyKBWQ4dg02DB5Pnz58tHH32k+lr8YRVxO9fVxIkTVb+Jz0a/idXMTbDt7bffVtYv9N34rp999plLHGTXrl0lf/78qt2wjL3zzjsSTCh0SFiAC/Ovv/5S9wNtNk2RIkWCdPNoM9USEvHAwnvoEMwr8LPcEDu4xWNsx/P/WYL9Baw1M2bMkGeffdZlrT0Aa8cjjzwiEyZMcLq733vvPalUqZKsWLFCXn75ZSU4Zs6c6RQJH3zwgXz66afKOoQVvyE6NBAAixYtkvHjx8vq1aulVatWctddd7lkh8Kd/+677yqBs27dOvX5ED/Yt2npQZvwHMB6gbfeeqtMnTpV1q5dK0899ZS0b99elixZop6HwIFlChmtcFXhLz4+PsGx+Oeff6R169bSpk0bWbNmjRJtr776qhJgJhA/sIDhGOC4de7cWTZu3KieGzp0qBJqsJ5j23fffZdgohlomF5OwgLTdWTNjAoE1irLgTQHE0KSQKFCiI69IWpw++23Iu3bX39cvPj15/E6PwKRARFTtmxZ2+ex/fjx43L48GH1uGbNmkrggNKlS8vChQuVuGnUqJGqwA5x1LBhQzW5grXj9ttvV6/Fc19++aW61RXaYd2ZPn262g4riU6Y+OSTT5SY0kB4jB07Vh5//HH1GC41WGEeeOAB9RiWHOxL061bNyXeIDbw+ZjYoQ+EtQftc8eQIUOkQYMGStzo77d+/Xol7mAV0jRr1kwJHNC7d2/1/efOnSs33XST+n4o7VGrVi1lLYJFJ9jQokPCApiDAWYZmIkEGqtFhxWSCQlDYGWAmIGogbipWdNV5NhYIfyFtwkK6LOsjzds2KDuw0Jz/vx5KV68uLKeTJ482dnXwEICSwzEQ6ZMmZx/cClt3brVuT8IkooVK7p8Biw3cDMhkQLASnL33XcrSw/AfgcMGKCsRzly5FD7hdCB6PAFfA8IORM8hhjEZ2jM9kHMQDwdgrXtPzcZFgCF6IHbDu7AYEOhQ8JK6GAmEgysFh0KHULCFIgZWHJM8DhAIqdkyZJqsNZixQq2o1J77ty5E90X3EFw18AiAzcYrB61a9dWVpozZ86oxAi4hyAE9B/2D9eSBu+zxixWrVpVSpQooVxeEFIQUNptBWBxwT5gXZk7d67ab5MmTQJWN8w6cUR7tZW8SpUqKgYSwgtthSvswQcflGBC1xUJC2AKBt50Hv6AFh1CIgS4teGuMsHjAFl0cubMqdxOECfPP/+8S5zOgQMHlPWkQ4cOTvGxePFil/fjsen2wvvvvfde9Ycg5zJlyihrTuXKlZVVBJaPO++80+d2QtigLVgmB3GHsOho4D5r3ry5tGvXTj2G6Ni0aZMKaDYne6ZVxg58D+zLBI9hhTKzVxMD5UIeeugh9QeRgzgkxELB2hQMaNEhIQcXIYLnAHzGwYAWHUIiADPwGO4qDLrajWUGKPuZYcOGqYwoWEH++OMPFUOI2BkIIFid33rrLZeBf9CgQUpIIOPqhx9+UAHJAEG7o0aNUgHBSFkfM2aMEj6IU4FYgFiBaJo0aZKyeiBYGBlJCCJODLx3+fLlqi0QD8iK0iAmBgHRSPDYsGGDPP3003Lw4EGX9yMg+O+//1bZVkeOHLGNU+zZs6eK/4E1Bt/v66+/VsfGjP9JDMT5IOPr33//VfvA8YFrS7vZggGFDgk5WuQAa5ZDoKDQISTMQZ0cU+TAglOjhmvMDp53V2cnGUAooNwFYmvgaoGbCJlL9erVU1lSpiUCYgCvhYXmzTffVAM7BBLAYI46XYhrQRwL0tN/+eUXZTUCCDqG0ME+EMPSokULWbp0qQpa9sbFhsBiZGuZbivwv//9T7mM0I66desqYYF9m0CswCoDKw8s6XbxO9gHApjhIkOSSL9+/eSNN95wCURODKScQwgiKwsuNwir3377TVmhgkWcI8ZLwqIyLiLQT548GbBqvCTx1HLMEjAj0dkLgQbmZQTnaXARmqZfQkjyJzCwUhQrVixpBUB1HR0EtVrdVNrSkyePyPTpqA/h17aTyDiPvB2/GaNDQgrMqag4GkxrjtWKBGjRISTMgHiBiEGdHGsKOUTP/PkwF1DkkESh0CEhBaZbTTCFjpXEgvIIISEAIsadkPFz/RwSvTBGh4QMFNxCwJ8mWIHIAL7iPHnyOAPiaNEhhJDohEKHhAwE9XlaWTyQYF0ZlCnHOjbg7NmzXMGcEEKiEAodEhKQyoh1WzzFzQQDLa6QcYAUT0KIf+EEgoT6/KHQISEBhbdQpdOs/YCqmcHGtCKh1gVdWIT4tygnFqUkJKno88da5NUXGIxMQrqIJ0qkQ+BgCYhgLOZpxeous1YPJYQkDdRoQQycXvMIMXjWpQwI8WTJgcjB+YPzyJdKzFYodEhI0IvRoXy5LiKFglmhFjqo2on1aIK1FAUh0YxeGVuLHUJ8BSLH0wrr3kChQ0KCdlOhaiZmeqGyotgFQKO2D4UOIckHFpz8+fOrDEcsZEmIL8BdlRxLjoZCh4QEvYqudSmGcBA67JAJ8S8YrPwxYBGSFBiMTEIqdMxg5HAROrpthBBCIh8KHRISwtGikylTJnVLiw4hhEQPFDokJISj0MHicIAWHUIIiR4odEhICEehQ4sOIYREHxQ6JCT1EcJF6JhFqHLkyKFuKXQIISR6YNYVCclK4bqsd6iFDjJBsOYV2rN161a1jUKHEEKiB1p0SNAxY2CSU9bbX6DGR968eZ1tYYwOIYREDxQ6JOhoIYH4mBQpwucU1NYlWnQIISR6CJ9RJsgMHz5cVePF8gMkuKxatSos3FZWaNEhhJDoI2aFTpcuXWT9+vWydOnSUDclpsDinfPmzQtLoaPbc/r0aXVecNVlQgiJfBiMTILCqVOnZNq0aXLt2jXntgsXLkg4WnROnDghv/32m6xbt04effTRUDeLEEJIMqDQIUEBwmHjxo0u28JV6Gh27twZsrYQQgjxDzHruiLB5ejRoxLuhJsrjRBCSPKh0CFBwW7l4nDKuHKX6o6aP4QQQiIXuq5IULCKmltuuUVuvfVWCXeLzvHjxyVXrlwhaQ8hhJDkQ6FDgi50MmTIIM2bN5dww07owOVGoUMIIZFLePkOSNRy8eJFj26scADtypIli8u2Q4cOhaw9hBBCkg+FDgkKZk2acBU6IHPmzC6PDxw4ELK2EEIIST4UOiTgoHbO+fPnnY/Tpk0r4QqFDiGERBcUOiTgoF6OXq0c1px77rlHwhVrPM6xY8dc3G6EEEIiCwodEpRlH7Qlp2/fvlKoUCEJV2rWrCmFCxeWJk2aSLZs2dS23bt3h7pZhBBCkgiFDgk4M2bMULdly5YNu9o5VtKlSyePPfaY3HHHHVK0aFG1bceOHaFuFiGEkCQS3qMOiXjgsoL7B9StW1ciiSJFiqjbXbt2hbophBBCkgiFDgkoly5dcsbnoH5OJKHjdbCaOSGEkMiEQocEFJ1tlSpVKtslFsKZ9OnTq1szY4wQQkhkQaFDAooWCVo0RBK6zci64ppXhBASmVDokKAIHQT5Rhpmm5EiTwghJPKg0CEBJZItOsgQ08UN6b4ihJDIhEKHBBRtCYlEoWO2mxYdQgiJ4tXLhw4d6vOOUYvEWk6fxFa21ejRo+XgwYMRL3ROnDhBiw4hhESz0HnuuedUNVtvF2NEJVmU+afQiV2w6rcWOZEao2MKtLFjx0rLli1l/vz50rhxY7nppptC3TRCCCH+Ejpg2bJlkidPHq9eS4FDrOtDRWrWkrkA6eTJk9Xt+PHjpX///iFsFSGEEL/G6KBTz5Qpk9c7xXpGOXLk8Pr1JPqFTsaMGSUSgduKEEJIlFt0fJ299unTJ6ntIVEmdLJmzarWuKpevbpEIlmyZJH9+/eHuhmEEEJCkXU1cOBAzniJLTpLKT4+Xq0EniZNGolEmjZtKhUqVAh1MwghhIRC6Lz99tvOBRsJsbPomDEukQgsUvfff3/ErdNFCCHED0JHL9ZISLQKHU2kpscTQkisw4KBJCBEm9CxWnQo8gkhJAaEzvr166VIkSL+aw2JOqETqfVzErPoXLlyJWRtIYQQEoA6OnYg0JQQ08qxePFiyZs3b9RZdKxCB8HWqVOnDll7CCGE+NGig5o4R44c8XKXIoULF5adO3d6/XoSHezYsUN+//13+fbbb51ZV9EidKzfw1oniBBCSARbdJBCPm3aNJWB4g1Hjx6N2Eq4JOmYGXjR5rrCSuYmXOSTEEKizHXVsWPHwLaERDxxcXHO+2fOnIkqiw6FDiGERLHQuXbtWuBbQiIeM0BXCwFadAghhIQSppcTv2Ed/GHhwRIK0YBVsF26dClkbSGEEOI9FDrEb5w/f97lMWK6UqZMKdHAbbfdJgULFnQ+ptAhhJDIgEKHBEzoZMuWTaIFxBo98cQTUqVKFfWYWVeEEBIZUOiQgAmd7NmzS7ShFyelRYcQQiIDCh3iN2JJ6Pz111+ycOHCUDeHEEJIICojIwtry5YtcujQoQQZWbVr107KLkkUCp38+fNLtGGmy8+aNUtq1qwZ0vYQQgjxs9BBif+HH35YVT62LmyILBsWCoxdrEKnRIkSEq0WHUIIIVEqdJ555hmVgTJ16lQ1YzeLxJHYBZa9c+fOqfvVq1dXf9F4bkRLAURCCIkVfBY6mzdvlh9//FFKliwZmBaRiOTs2bPKwgdx07BhwwQF9qLVogMLZrSk0BNCSDTi82hUrVo1FZ9DiMnp06fVbaZMmaJW5NgJndWrV9NdSwgh0WTR6datm/Ts2VMOHDggFSpUkNSpU7s8X7FiRX+2j0QIp06dUreZM2eWaMbqupoyZYosX75cOnXqFJWuOkIIiTmh88ADD6hbdOwadPDabcHZbWxbdKJlyQdfgpH37NkjJ0+ejKoCiYQQErNCZ/v27YFpCYka11U04y7rigUECSEkSoROkSJFAtMSEtHEikXHXdYVhQ4hhERRwcCtW7fKhx9+KBs2bFCPy5UrJz169IjKuinEO86cOaNuadEhhBASTvicHjNjxgwlbJYsWaICj/H3999/S/ny5WXmzJmBaSUJe/RAH+11ZhCHVqdOnQTbL1++HJL2EEII8bNF5+WXX5bnn39eBg4cmGB77969pVGjRr7ukkSR0ImFysF169aVVatWyYkTJ5zbaNEhhJAosejAXfX4448n2I4srPXr1/urXSRCQLbd77//rsoNxIrQAVeuXHF5PGnSJOXSJYQQEuFCJ3fu3LJy5coE27EtT548/moXiRDWrVsnixYtcj6OVaEDxowZI4cPH06w5hchhJAIcl09+eST8tRTT8m2bdukRo0aatvChQvl3XfflRdeeCEQbSRhDFawN4kVoZMrVy5VP8fKJ598ItmzZ5fu3buHpF2EEEKSKXReffVVVf32/ffflz59+qhtBQoUkNdee42de4yAwFsUhkRKuV7IM9aETsuWLWXWrFlqja9du3a5PHf8+PGQtYsQQkgyhQ6yThCMjD9dOyXay/4TVz744AO37plYETo5cuSQ1q1by9y5cxMIHUIIIRFeR0dDgRN7wJLjKQbFuvZZtONO2OklUQghhESA0KlSpYrMnj1bxR5UrlzZYweOBQ5J9JJYGnWsDe7uhA7ce7Fi3SKEkIgXOs2bN3cWgsP9WBvMyA0uXLgQ6iaEFe4sWBCEFDqEEBIhQqd///7O+wg6JrEJAo+nT5/usu32229XVbJjFbiozEKC8+bNU/dZKZkQQiK0jk7x4sXl6NGjCbajSiyeI9HLr7/+Kps2bXLZliFDBkmRwufTKCrr6WBpiIwZM6r7rJRMCCHhgc8j1I4dO1RAqpWLFy/a1hUh0cOWLVsSbEuVKlVMCx0tbDTaXUWhQwghEZZ1NWXKFJeFPbNmzep8DOGDYOVixYr5v4UkbLCLzUqZMmVMC50yZcpItWrVpHDhwuoxhQ4hhESo0GnRooVzsOvYsWOCgMyiRYuqIoIkerETNLFu0cF3v+uuu5yPKXQIISRChc61a9fULaw2S5cuVSXwSWzhzqJTqVIl+fvvv6VgwYIS61DoEEJIhBcM3L59e2BaQiISWHQaNGigRE6JEiUk1tFCh1lXhBASwZWRsb7P/PnzVel768yV611Fbxq1nZUCFh24LitUqBCSdkWiRQdLp6RLly7mqkgTQkhECJ0VK1ZIs2bNVE0VCB6s+XPkyBGVZpwnTx4KnSgFA7ddth0sOuQGWry4EzoQOUOGDJFMmTJJz549g9w6QgiJPXyOIsVinvfee69aoTl9+vSyePFi2blzp9x6660yePDgwLSShBzrKuWmRYd4b9FBeQZw5swZv382yjuMGzdOTTwIIYQkUeisXLlSzUSRbYJBDvVz4uPjZdCgQdK3b19fd0ciBFjv7KBFxzehYwZ0m1WV/cGoUaNUQcfvv//er/slhJCYEjowzet0YriqEKcDUFdn9+7d/m8hCem6VvhNv/nmG9m4caPalj9/fqlatarzNbTouILYG+BuhXczFT9QmVl2lcsJISRW8Xk6jtXLkV5eqlQpVfK+X79+ylT+7bffys033xyYVpKgA2vDBx984ByMdbYdYrFMKw4tOq5ky5ZN3cK1a4dpxcGx1YvlEkIICROLzttvv61m9eCtt96S7NmzS+fOneXw4cPy2WefBaKNJAScOnXK1uKAJQ9McUOLjiu4HtwJHVh5/v33X+djuH0JIcEH5R/skitIdOLTdByzUbirtOUG962rWZPowF1AKyw6prihRcfeogO3H4QNAvY1sHru37/f+ZhCh5DQLMSLKv6YtHXr1i3UzSHhKHRKliwp69atU64rEptCx4wzoUUnYTAyOlAEb584cUJdM9u2bVOxTqbIARQ6hAQfeB9w7eEPVh32YdGPT0IHAxwEDoIdKXRiV+iYLi1adOzdVxA6uE5+/fVX2bdvn+3ruEwEIcHHzHyE5RUTExLd+ByjM3DgQHnxxRdl7dq1gWkRCQuOHTtmu93aKXA2lJDcuXOr299//92tyAG06BASGteVKXRI9OOz0OnQoYMsWbJELeSI+ANURjb/SHTgrgOARceEFp2EFC1a1FkF2RPBFjpYmHf58uUsKEhiGnMdOk42fAexh5EmEH0epZBybLeKNYku3HUA1nRoM16HuAodTZUqVZTAsBLsThZtmDp1qrrfv39/j689efKkskaVKVPGp+v90KFDkjlzZpcgbH+COCfUKtLZbcEAcVawAoTb2mSIL1m1apVaTBexYbAg3nLLLVKkSJFQNy2sMV3GkTZg++t8jkviGI7rAMWBwauvvhox/b/PQufRRx8NTEtISMGgiyDzsmXLqkHKGj+CkgKIO8mZM6cKrtVQ9CYkS5Ys6jjpwn3169dXx9YqbIIdo6OXn/CG4cOHq5nv/fff7/WCrRA5I0aMUELnhRdeEH8D8aVLWCQm1PwJ4qwgErt27ap+13Dhjz/+UH843yB2ULUef8E8NhocH3z2gw8+qNoTLDDwHjx4UAoUKOB1X2RadGJN6KxevVoJ4tatW0vhwoWTVHZEg/4sUBMaf+OzHENMBjo0K+jUGa8RHHD8p0yZorJ6/MXff/8tv/zyiyxcuDDBIHzHHXfIE088oVIx6aryjrvvvlvd5s2bV8U12V0b/rTomIUI3S0t4UvdED0YYEkJbz978+bNXrnskopdvxMMtDVOXxvhAsSzHnw8HRtY5n788UeX2k5bt271ayV79B3YX7CXH5k4caJ88cUXqv8KtkVnwYIFMmzYsICsWxcoJk+erCasP/zwQ5Lebx67SEqm8FnouOtE0WnrdX5IYPnqq6/UKvJJPVntwKwI6JRofRJjZt6kSRNloqTI8Z5ixYrJ008/LQ8//LB6bHfs/NlRmLNUdySlQFpi+8Vir6hJgoHOnFEnJuJwjvkqiLwRc7GEOUh7+m2//vprJYp0f4EYizFjxsjo0aP9XjRv7969KhYsMfD7uVs/zxd0Ac5FixZ5ZRHEZyZm0cF5CVHrbiFjzezZs9UEHysFRBrn3Hw3TFawSLc7zKVtIim+yeuRa+jQoeoWnRkUdKZMmZzP4WKBCRX+fJLw4sKsGH5zmMCrV6+u3EOeQEeBi9Hd8gD6ZPOU0ZMYGGR/+uknJWAaNGjgdLPgAoA5WA8kXKIg6eTLl895v0WLFmrNMBN/ms3NTJLkvMbX9yBGBIMHrB61a9d2boeVQWef2fHzzz+r97Zp00ZuuummRNuxYcMGF2sBrpFgW5D96abVQs+dm+eff/5RyR0QzO4wBxpP4kILal3HyXwfAtNhdUwKWBYGQqNhw4YqSUEPnjgf4L50B9ZHnDdvnno/QiH8EVOUmGCDNQtjGFyPSKTxdA1OmjRJuXnhor/vvvvUWIf76CvhIrQK7UicAF6zOV9gmRo7dqy6j6Wd7M53U+hEktsvlS9ByPoHHjlypEsnA0sOAjCxnSRU/WvWrHE+hnk3MR/6d999py6s559/PmD+bgwc+APolHQ6OTopsyO0C8BE0Uj4eSPFPxsOYMDq3r27c8JgdhrodHA/OfU8zFkqrlHs0xooGAiLjvkZpoUmMaEDkQMwQcI5h1tYv1Bt3Q6rSyQQhd4g6oI1aM2YMUO5W9q1a+ccPE0hgEkRcNdXwHVtilBfflvzfbDk4jfGsdRL+3iLFu7W2loYACGg0I9YJ7979uyRL7/80sX9Ewyh89dff6lbTOjMAdpusNaxbOiDP/zwQ4mPj3e6+fr06aPGO9MiEop+EL/ZqFGjpGDBgnLvvff6ZZ9nDBcc+iNrhq3eHtUWHb2oY7169ZTiDWbWQ6RX4fQVHeyLWkU1atRIdhswk5s1a5aaeenOTK86Dw4cOOAc0HAB604LF7SdqscAhngdFtryDevxQqcBywZmxehwvbVueCNI8BvqldQDKXRMIWyumo4ZNGaHEHiwYnoCbi8AC+Pjjz/ulYDBYJ0UVzliUyAyMDhgANMgkBbiAcHXnhYnxvHwR/aVjimZNm2aCnI2MeNt7DJksA2uaxNv3EV2vyncFIjbALfddpvUqlVLsmbNmmhcjHkuoZ8yxRP6EC2CUHPNHDA3btzosi9/ZbIldm6bCRTmgG43WOO8MoWbGcsEAZ8rVy5lqQ+lGxXHHCIVf3ZCB+2EJcpTVtSGDRvUNaC9M+Z3xuTD/N30d4xUoeNzjM7cuXMpcrwEJ4dd4T1vXQj+8p/DN48L3ZxJmUIHfnWzE4DwAZ4GEpjV6dbyDWunDhGMAVbPKsePH+/83RGHpVPB7cBgDSsrZsR25xSsq1aRYn0NhAkGd+s5ap53iQkdc4DVcV56Bg1/Pyx/njAHcQhyq3vP3UCS1GsDsSk47tpChM4eQguCE59jPebmZ0OQYlFjuJX8hV2cSmIBn3bxFd4eD5wD5nlg9gPLli1TIsYTaA8GWW0NBta6TOZSJ7gPt6Zun9WCEiyhYx5H0/JoZ9Hx5HbTWUem0EmKSzgp4FrT38PMfrKyZcsW1Tdoq6C76wjXAMJQNKaIMc8xfC6yHTF+mNvdCR2c0xDQ4RRHlyIpJxRMZjAzw0KA1Fnzj9zAtI6Y2K1sHci4C52dpQct3Jqdk3Uw00GLFDLBj/FAJwHzPjoKDDx2vy+uwcWLF6vODq5RdCjW3xAdsRZB5vtMIKxgGdB+eY25L/35eK++j8FBDxDm+W12fN6e41bw3a0dpNkBW9tlgrbotcUSi1/T+0Rnb9Y4srohrJYKYB1AMOP/+OOPvcpQA+b3sxtozcHE7nm7bEt3A5/V0oPzwvx9re/D8U9u0DtEuCksISS1MLKKNG9chfi+2A+sce5IzKJlttsUOnai0ZOlWgucUAgdiJLBgwerc8JTID/in4DV6nfZ5rczv4d5rpkCXE9+cZ6bHgp3QgeiCBM104oWcUKnR48e6g8dH0y8COwy/4gk2tm7W17BesF6cwG98847Lh0LTm503uYAZAaO+5LtwCw6/9OzZ0/lnnEHrGvm725XxVhb3DSmCDGBUEJ8mB6ATaGD80zv23Q52VWOxTnzySefqPgiPIdBHY/xmd5kjsFdYe103YF2WQcfu4HdKtpgOXj33XeVZQarxH/++efqD8cTHTisEOb3cmdJwCBnChF3A7vZJlxvuKbHjRvn1Xe0ihfr9WiKD7uB2BycEutDrJ+FfZvnivX383TNW+P3fEH3edY+0Ruhg4wxiBy7opt24PvhteZv5E7o2PWFnkSTPvbmb+SN+MN30Nfb+vXrfQ5pwDmJcxyfhfPctCpaJwbuzu0LiQQPmxMKWHqROo++BmEPdtZ/nAs4lh999JHMnz/f2RZ9bGBZChd8jrzDLBAXdrNmzSRcwAwLAwhO0N69e6uaL+EudNxVp/Q1wBDfGbN7pIDrExSggBb87QC+Vn1xY3blKUjUhELH/0B0YjDFb29n2oX7x4yRQIdoZm8Ba/0TdH7uLAOIp9NVTM1zS1c31aBTgkhAlolZFE9bb/RABVebHnwx0/Omk9cuueLFiyeI/7AbOHGumrNqO6GD9mCQxnkOsCwNrgUd5Axg1cEsGDFCSD82A50RB2R3fcGiAfcZlroB7tb0+/TTT6VVq1bqO/mafWKtu4LvZ35fcxBFajieQ6kCfT3qwRbxdqabSGNmbFqtYXbWMXfXPI4lztOKFSuq2Tksf+6CxRMD8WJ2rny7a2DmzJnqvMfxxaCtj4e3xxnWH1hEy5cvrwoY4nc2xYt5ztnVwPF0TutjbwpEbPMUHK/rGIGHHnrIaTHHNYn34vc1jztEDLabHhJP3x3HB21G7BAw92UG2F/wQejo1d1xnnt6DazGOH9hRapTp46LMLcLZo4Yiw4OIrJuwgX8kKj1MmfOHDVrfO+99xLMUIMNTnqIQT3IWEHcgrvZn53bwFvMmYq72QyOkY6bSMw1RddVYMDgYQ0UNgNRzVkmYihgCtYBoxgYUN3UBOdZYoXa0MmZA7tVYCDTT6eXWjO4TMFuDvzo7H2pBWRnnbCbUcPNgVmvxm72i+sHFhst+jx1qrrGihnki89FRo0GA6KZdYPn0QYEC7s7nr/99putVQLXoaf4BKvbweo+Mi02+J0gDkxBo69t61IjGvP3sx5zfDdP9bf0IIn3ITgc5x32BysZzh87YeUN2IcePE3s+jgd3wWLpHl8ILbQHk8WF5yrug6MrrHjqR9Fu6znsCeho9tivgbXI8SVO0yrrBkDNWDAAGUl1dW+Ac4bTNz//PNPF8utp5o+sKigkrm+Ts3jo88lh8Phdh92gcbegN/SPLb4XLMf0s8Fy7XnV6EDywkObLgEGmEmB+WONDvMlps2bZpoAGSggbJPrOooLmScGOZxxCCiq8taFTg6XRRm09lvVrAfMy7BvLjcncDWRVgxmzehRSdwuBM66NzMwR+dPQYXdKboqHB+4LEpQk0TMWZ1WFsLtZFM8F5vA1atHb05EJvBq+jofBE6dueh3TakmmMw1gOVp3gbWDOTGpysZ/MQSdZ6NRhs3FVC1plxmFBZs7BgBRkyZIgarHU9LHefaz2+OJ5wAdhZsMxzQl/b7hZRNj/TenwTK26nv4vZRk8F5LwFg51dHBG2ow+H6MT3Mn9HXfnZFIa4DvT5bh2D8F5zkquDihOzOlrFtqfX69daz3vr8ir4LbWwMMWd3aBvttk87uZneBI6+jhANMFVbX4fWJIcDoeaDME6aIf+vr5aJq3nFtyLZv+A5/HZiCvC2JWYyzWsXFcwVSHzCjMdCAyrP9CdFcMd6NRghYG5Dh04FDuKq5lAreI1GAQQB4QYgdtvv93ZCULkaHDf9COGCph4EysNjlkAZmudO3dWr7VmPOgTDxeKnoXZZaUAnEhmIUIcS1z4SC/0JHTMGRoWBIS1AS4MQKETOBD0alpKHnjgAfX743xwF2iIa0AHy6L2Ctws5uDx7LPPOt2SepDQ4BzwVpRoM7vGnbXIV6Gjz2dvRQmEBkSFp+tZD8KJVbH1BFwH1iBkaxq0Ca4TTGTwmbA2mRYdWEEA4hogwjDoNG7cWFlS8ZvhvrXD17/h9OnTlWvQDlz/LVu2VH2unkSZ6fEm5mBqZkd5g+7PzXPQnfvOE8jMNc9vDKZ2Ax22a6sZvj+uA6ugswokPZBbzyMkyZjxiPp1iQkdfA6C+rGmG843T6/Xfamn8x6DO6w0eC1qZ3kTEwlRjGvW3fjp7fmN95vuYYyZe/bs8fgb4rpEX++rRQfXiDlhs8YTYn/405a8UJYj8dmiky1bNnXBwR+H2SMOqvnnKzgJIF4gZuyYMGGCck2hcBYCzPBaxKOEat0bb/EmDgbWGVz8MMnazXagkDFDNE3sno6jNvXqCxFBqJ5UunniYVYLkaPr7OA+K10HDnNgRR0MDGCwBKLD8+Qe0J0R3m8Voub1Z70WcX55a4W1Dkju3uer60qfi97E9QB00JjUeBoodICsXSet43cSQ5/73oLjrmNVMLi6i8XD5AXHDskC6K+0C01bqPSgDLGG9rsTORpMAt98801nDI67eBmzJpbeZ+nSpb36bhAP1kw8q6vUG6zWJrTJro8zY3ZwHE3rh3aNWCeMZvaoCa4b0yKu3VJ25xviF/UEGR4AiAwdS+bpnPYkdLAN3wFuVXwvvBZC0xuRgu9oFTk6ixDt82UtL+v1u8uwsnj6Tr4KHWCOL1ZPCvanrw1Y10JZQdrnTzZrsfgDuJrw5w6YgZ988kl57LHHnPVBUOsC67S8/PLLqjMzZ3y4r609dlj9xJ7qESQHd2Zldx2nu85fzxB9AZ0gvqM1ANOKGdegzbwoB4/jjcee6kmQ5GFWvNYDLWbBmBV5k5FhFToQSaZ11VpRO6mxFZ7w1aKDQHicm75UwvUmTR2uHutgUrVqVdWxerNMiq8zTS104K7wpVCfRrcJEwlMTtBnWYPDEwPWHHfiTPclWkTgXChVqpRX6e9om05o0HgrkHH+6c/GhNjaJt3XmgX5zAkr2gvhaILnrZ+PiSEGc2+yfCGE7dxFKHyrP9u0RFgDl63ge7jLNoRYgTXUvNZgSfGmcrJd6jzi5mA9TEwAuwNB5BCpWy37tgaxQ6yY2VKJgd8P1kl31kJc31q86+s31LX3fLboAPzQMM0iIlubOHGB+HsVV5xMcGmhXo8Grhg81rMjiBqcTOgs8Pkwg+oMJHfp2KYFyp35N7noCHhvwMWTHNO7lTvvvFN1OjBxm1YeT0LHHBghHilyAotpLdPmX186A3SeprDBY3PgM034gVp5Gx2ktcO3fq4VDAS6zocveFqeAJMfq6DB8bAOMBgYMXC4EzoQ+J4qyZodvS8TGasow4CSXIupmRmn0dl5WmzomT1ER3KXKYBLx25Gjn4C6yKhunPNmjXdxqCZMTqwptxzzz22n2MdQ6ylFLTLBDFQZlkNdyC20S5eEr+z3blqns933HGHbT/41ltv2U5GIKqsEwqIMqsb2Q67MhIgqSLHPB+2W+I6rWMeBAm+D8Zyb6wuEIPw6CT2ubgmtYUt4oQOlBp8mc2bN5cuXbo4f3Dk9vfq1cuvjdPBadZF5/BYXwD4YRCfAoWOTgzB0nadgAZrlaAD0H+JBQ0nFQQ34mRANodujzvxYAau+QOcaLfeequ67+kiw0BZrVo1dQwRO0CCB9KSrSLTl84AItW06FizjrwZsJMLTNVWF3JiQgckpbKwGYfnDRjYzWOCx4j9Q7/lzsWFW3NhUmtAtwbHPSlueqCLqMG1jT7CXbYY2oIQAbO8P6wyGjuhpYWvVeigrckVOoULF7YNoEffBdGGPs4cQK1uVbRJT4rRD3orFO2EjsZdYgasm7qvxYQcsT92wtbuNzRdaegT7777btvPsLNyJWei76m2WlKAl8TdMc5vmTRgwqIFCc5J8xqwA2Oyp9/PLIeh3Z6hFjo+u65QLBBroiB2xBQUuCgxIwoFyBayZgy5A6bzYKRN4+KvW7eus5OAoEJHAFec3SzPWgzLNAP7Cjo1bTq2C2xFO3DiIS4Es1xYyCJxBd5IBscb9Z7w++hOw1eLjjmY2A1kd911V4JO3h/gundXwiFQAYeIMbFaJz1dI3jOPCZ27cKEDdeoedxhIcYAimKo6OcQpBoIoQMhAzGK1btRfNEKnoPrwRSShQoVcjtDxuCur2FYFlCPCfVotJBOrKYJBjh3wgHg+6LftA7mpgUE+8B6bRBxVjcZficzvszb/sZTlpid2ECIA9qK8AZr3wexjLUDtSCzuteAPq9x/qAP9+V3tgbt+9KHJ1foIBFFu5LwO+E8drfPm2++WY03epKvLToAxwYWN5x3OuvRCsrL4Lvh2NgFmFsNE+6OdTDxeXRDfj/MhlbFjpoO/s52gvsHF7C5hg7AY2sRtXAGs4ty5cqp+7C0INgPF7sOsrQz5UN1JxZE5g7MvKwdO1KOtZiCuNEWH0CRExqsVgpz8MLA5MnKZ3Vd2Q1ksNZhtpYUV5EnMOC6Ezru0uaTAwShXUeJvsHdQAJXiXlM7IQOBmTrzBTHFeLDE+j7vO24ce2bA66Ol9C/vekyNmMn9CCOWCBYhTGomMfcKnRwDevzwbpelzvXFcSUjkexPo9ZPTJiNfj8xH5bCAOdem8VIbqODsB+/LXGld2Ygd/dTqBAAOh+2N3EQter0uObNdbNFzDgJ7akhia5td/MdupzGrcQ8vPnz1dtwYQbQgbnSqdOndR5gskD+ggtQnW8IK4Nq9CBQMIxQ/ybFutWoYNtOGdh2TVFcVILTfoLn+3buDDcVRT1d1wHTjYMyOasCp+Px4mtiByuwDeNiphw8+kTxg53LgCcrJg5mResFXRa5vtxgZvFxZg2Hp6Yg25i15JdjI4dplhCLIUVXwccvN7aaZkW0qS4zOzM4NgnOmesJg5RYNdOWCLxOvQFsE4iNs3cp3lMTLcPqrojmBJC0FcwCEBgeRr0TQELkQIXv87C1Fkq+jXurMtmQCy+G8Sl+Vqr0LIGo5ugL7AbsM1t5rHCb4gBUluk9Wt9EbHWgF5T6OB7BEro6O+BqsKw3phYJ3SexKppfUoqpthyVzVZ423JBXOCagKBgtR8WFtN9yzCJ3r37i1PPfWUcmdBrGj072muXae/r93vg33juOq+yXS/m5MSfNdnnnnGZXvECR34Lc10Z1z4UG5I/07KshB4LwKudNAVzKe4r60ZSC1HBVQUO4JpDjVnYAnRWViRjCfB4e45fG+Up4dZ3VuLDlS4uT8KnfDE7HitHSM6aXMAtcbouOuQcZ5gXxDVdlk6ppkZr4MAeOmll9SSA3CdaHSGC2aE5oDbvn17Vb/H3QCHz0dQp69CB8cCnSXcS/r7W8Hgj8xL9EkQFOiEcX0gGQHmdVMsmhk6OBbt2rXz6jrQFgoN2oHjaB5LDApmx272gxgctJXBm/IT+ne0G0QwWcHzmOzo4wH3JIAgdDc5QhaMjqMxMfsI8/zRge1mthL2bf7uaAMGcizu7K3Q0YMp9mP9PZMrfPCbt27d2vm74Ds0atTIRZxZPwPfyZ0A0e33tuwAzkGrCMEx0rjzQNgFx3vCXRA3vi9EDILCra6jtGnT2n4P/ZvDmqSzo/Q2u+vNug9rkU1cb/o1VgtqYkIv0Pjss0DgL04qWBRw4uJEh98YF7O3i9qZwHSGQGINhA3o2LGjWgEVa4PAf4jZKALTcGIg7sDODxhpeOpoE3MnuZtdYTtONrPTwwlodlJc2iE8MTtic5DB74fYOKQz61Lz+J3N88fd7BQDKoSLu4HENK9jRqiLTlpjgCBWkFWDtpjLM0B8meeTOcChc0MwJ2LQdAVjO+zcDNYaOHaWIrvzGIOLHmAgLtB/4Ht4EyRtBwKYUW1WL/dgusrg3kFVX8QmQqzB5YXn0Xa4BhC/ojt883rFb2V3fWMfWN8JWUV2M3fs47nnnnMZNGCVwmu1EMb3NC3gEHR68MKxQJ+qrXzm8TOFjm6rWYYDx998PWKmzCwrK9YZvLlcgF73ynoOuMs8MoHwxeci087EnZhGm7XAsl4DZkkHT0BE4jfxlNgB0Y3wAPzmOiYKvwvOfYxxCCq3thlgYmFmVrlzWSNDz5MFMiku43T/vcdMP9fb7PoL628G4QwxifEZ5701EwuTDVSxDoex2mehgx8Ugcgo5IdbWGSwGvMjjzySJDMfzKOJ1WmASsVftGEVOmZ9icRAR4JBzJrmqH8D88THfmnRiSwwKFSuXFnFc2EigMHNnBXisTngeSqT4On3RlyIXmDUOkMzOzszpdosLGYVG6bQ0ecgBnYMFu6qvpqWF7if0Gkmlvlh99l2JLfoJb4DLEBa6Jjgd0H/pWexZn0g6+9hXo9Wywomiwi6xQQSx8KTBczut9SiCZ+BQFLsS6dym5+L/uLFF1+U119/PcFzpvDSfQhEBfalzztP1hErEEKwPuAcRaV1M04Jv5vVVQNxYAoOWK9wDuCcgVBBkDNiRnDMIQbgzvMmbds8XnbiEr+FndDBOKfBccAfJveoa2YXF6aTXHCeWy09yN7DcYBbB98BpR70mGd1UeMcwPVlJhFAjEGkeiIpQie9zXjtyaJjN05b3YMmSFBCGRhPde2CRZKiUHEQIGzwR5KOtdNCJ4QaDcAaMY8L35zl4cKBGw/B4ViSw3qiWk3r5qBAoRO+wNWjrZzoWNGR6IER1oFu3bo5fz/zHEmqDxwzaT2JsHaWZmdnnjNw5yBt165WlOmGMs85T5XCTWsLBlVP8WcmwbRMaqFmxvr44tow22oVOtindb/JAX2AndDx1CZTuOj3QKzBLaktbqbQTszijOOC/sraj+FzYB2yHjdrDBFeA+GN8x2TP5wjsH5oYH2EiMWyKZ7c+O6+owbnsHW5D4gcCCsrOC/xmbCKWcuSeBJ+ekICixv+UOVYiyWrpdHO2mXGngG4Z5FYAneTDvFIyrWQzubc8BSj4+v6lhCk7ko0hL3QQcE9mKJgmjVBpWJYFxD4RLzDPDlx0pmdh3miYUYMV54VdBa4CDAr0NVMTTeAzqpAJ2oOVHRdhS8IGERwrRYMVjFhCgl0uujwtFUmKUA8uauJYe7TPH/QJqzhY/rh4bZB4U6Yr3UxT3cWAwD3gy5iaO7HUzAzLMdm5dxg1ArSYEYPgZLUtHJz9uypzpc/cOeSMq0tSP+GiIbF0GqJM99jClTTYuVtRWjrgKnPA+v5aj2uuq6X1Rqtwfth1YGb0FNQsfledwG2doLGXZkEnHO47qxCx5fJoykYrBYdXDOm6wrxmNZV6rV71rSQ+suik+a/7+GtRSdS8FnooBqyuQ6KBicdVDCFjveYF4e+sCBoMADA94mZCqw1ZtEwuwveHAzN+8j2QLoq4i5Ml1igMh5I8kFH6m0xNfjAMfh7s66aOzzVvTE7Nus5Y03NxSBoXdrBfGx2nDhnMdND3AZKRUCII/YO5yq+kzsw0zbrhQQbb9fOctd2fG/sw1x8NxCYIsJuUoN+GhYFs/8x+wd3EyFTjHgrrK3njRn3g75Ku43MYO3nn3/eq7RutCGxayUxoWN+J7ghEa+CWBtPwLUI1zICxmFRgehLav0oLeT08YdgMUWk1Z1sYrqukxOjY6J/VzuhYw3Mj2qhg4Bgu3Ls6GwDsZ5OrAkdqHet4DHb8NTx2+3HFDroBHRHgM5LrykWiFonJPigUzJjCQIpdHyxGCFjCwHLpsnd6gbD/syqs3YVi5OThhtu4Dru27dv0GtW2WW74NjbVS7GQAY3jqeYCogkvEZnwyWGVVyYgzhcUVromNY5f1qczX3ZiSctXnGcEBPkaWkDDdqqY8iScv3BI4Kac/o7myIT16O3osm8Pv3lurL73TD5hvUnHIKKk4rPVx38trA4WJUmtiVnxhOLmJ1NcuoMuBM6VsIhKIyED9YMHn+ZqhHLYU2nNQf45KSa2i3QGCmEY2FOxNFAlMKKgUBmBP56qqIMMeTLzN6TexGiAhYRq2vGnzGEiWUmYrCHBQnnZFLdv76CUghz5sxx1oLTRSVhLIDI0RltiS17Yl4LSTlmaT2II/NcRfsC7W4NND5feUilRAeJGQDqVgCkMyKFFetMEe8xlbun6PXkBDoSYlfxFq6jxOJN/OmTt7qukgrjyxLHl+OLzCjU/dGCJLGlIpLSFlj4EPJgBVYCFLID5vIe/hQc5nnnbomV5FQ/Tgo6C9EsYYAMZoQr6O9vLmTtDjP7KynHLM7yHggsO4tONIQ6+Cx0kBmEAkOIxjf9iojNwYKZkcLw4cPVXyhN4bjAYLJH5H1yFj1D54S0RKh6DgTEE4hDQDyMN0sYIOARy734A9OKk5wgYgSpYvabWBFC4j2BDuqGdU/HVrlbhT65i466w4wJSmotpUCDWB+7ApHBtG42atTIxeLvLhEhZoQODgBWKscyBjhxcYIimDDSBlgE6uIPaZhJzaTwB75WxnQHO37iDbpImjdghodqs/5YV87sOJMzY4dAQwA2cQ9io1Bc0ts4mmCAhAq4Yty1CUIIfaGnelBJwaz5FCzXVLDwh9B5/PHHVdwtXJjm8THFb0xadEx17GmtJkJIZIOOLxAZQsFMC49FYB1A3Ek4WTAwIfZUSRnnhLcB6b4AkR6qLL1IEDqFChWyDajGZAhLSuB3C/XyDSEROlhnauDAgSouB0u5W+spbNu2zZ/tI4REGRQ6gSfYcSfhig74TW6F7HAkkC6luLg4VZQxWvBZ6KCMNZZ9x2J+ySlURgiJTdhnkGABt4s3y4lEIijP8OOPPyaonEz8IHSmTZsmU6dO9WiGJIQQd9CiQ0jyQSmRZ555JtTNiAh87nHgu/O2cishhFihRYcQEtZCZ8CAAdKvXz/bpeQJISQxaNEhhIS16+r9999X64GgHDQqWlpTz7DIICGEuINChxAS1kIHVRwJISSp0HVFCAlrodO/f//AtIQQEtVgLaUFCxaoKt6EEBIs4hxJXNDmn3/+cRZiKl++vFSuXFkiEV0Z+eTJk6w9QUiAwfo80VBplRASOeO3zxYdFAls06aNzJs3z7lezokTJ6RevXoyfvx4tQIrIYTYQZFDCAk2PkcFduvWTS2qt27dOjl27Jj6W7t2rVJW3bt3D0wrCSGEEEKC4bqCmWjWrFkJ1rlasmSJWlkY1p1Igq4rQgghJPLwdvz22aKDta3szM/YZl33ihBCCCEklPgsdOrXry89evSQffv2Obft3btXrZbboEEDiRSGDx8u5cqV4wrshBBCSBTjs+tq9+7dct9996kYnfj4eOc2LOk+ZcoU2yXfwxm6rgghhJDII2BZVxA3qH6MOJ1///1XbStbtqw0bNgweS0mhBBCCAmXOjrRAi06hBBCSOQRsGBkpJAPHTo0wfZhw4bJc88953tLCSGEEEIChM9CZ+LEiVKzZs0E22vUqCE//vijv9pFCCGEEBJ8oXP06FFlKrICs9GRI0eS3yJCCCGEkFAJnZIlS8r06dMTbJ82bZoUL17cX+0ihBBCCEk2PmddvfDCC9K1a1c5fPiwqqkDZs+eLe+//758+OGHyW8RIYQQQkiohE6nTp3k4sWL8tZbb8mAAQPUtqJFi8qIESOkQ4cO/moXIYQQQkho08th1UmfPr1kypRJIhWmlxNCCCGRR8AKBprkzp07OW8nhBBCCAmvYGRCCCGEkEiBQocQQgghUQuFDiGEEEKiFp+EzuXLl6VBgwayefPmwLWIEEIIISQUQid16tSyevVqf302IYQQQkh4ua7atWsno0aNCkxrCCGEEEL8iM/p5VeuXJHRo0fLrFmz5NZbb5WMGTO6PD9kyBB/to8QQgghJHhCZ+3atVKlShV1f9OmTS7PxcXFSaQwfPhw9Xf16tVQN4UQQggh4VgZORpgZWRCCCEkesfvJKeXb9myRWbMmCHnz59Xj2NcLxFCCCEkDPFZ6Bw9elSlmJcuXVqaNWsm+/fvV9sff/xx6dmzZyDaSAghhBASHKHz/PPPqzTzXbt2SYYMGZzbH3roIZk+fXrSWkEIIYQQEg7ByL///rtyWRUqVMhle6lSpWTnzp3+bBshhBBCSHAtOmfPnnWx5GiOHTsmadOmTV5rCCGEEEJCKXTuvPNO+eabb1xSyq9duyaDBg2SevXq+bNthBBCCCHBdV1B0CAYedmyZXLp0iV56aWXZN26dcqis3DhwuS1hhBCCCEklBadm2++WRUKrFWrljRv3ly5su6//35ZsWKFlChRwp9tI4QQQghJFiwYyIKBJIDsP71f5u+cLw+UfUBSp0xt+5q/dv8lI5aNkEENB0n+zPmD3kZCCInm8dtn19WXX34pmTJlklatWrls/+GHH+TcuXPSsWPHpLWYkCikxugasuPEDnm34bvyUs2XbF9Tc3RNdXv12lUZ+8DYILeQEEKiG59dV++8847kypUrwfY8efLI22+/7a92ERIVQOSASRsmJXhu45GNMm/HPOfjbce3BbVthBASC/hs0UGhwGLFiiXYXqRIEfUcISQhDrnuIYan+MKVC5I+dXopM7yMy2uypcsWotYRQkj04rNFB5ab1atXJ9i+atUqyZkzp7/aRUhU0mRMEykxtITsO70vwXNZ02UNSZsIISSa8dmi07ZtW+nevbtkzpxZateurbbNnz9fevToIW3atAlEGwmJeGDJwd/MbTPV44JDCiZ4TYbUCQtxEkIICbLQGTBggOzYsUPV0kmV6vrbUTCwQ4cOjNEhxMBMaITrCi4rT5y/fD4IrSKEkNjCZ6GTJk0amTBhghI8cFelT59eKlSooGJ0CCE3OHv5rIvoOXnxpNevJ4QQEiKhoyldurT6I4TYc/LCDWFz+dplOXHhhMfXn7l0JgitIoSQ2CJJQmfPnj0yZcoUlWWFZSBMhgwZ4q+2ERLRnD68RwqeFNmbVWT1wdUuqeTYfjqtyKl0N15/9hItOoQQEnKhM3v2bLnvvvukePHi8u+//6olIRCzA9N8lSpV/N5AQiKSkyelYOvHZf52kbqPiuzJKtJ5amf11K1X8siErw7JoYwid7W7IXZo0SGEkDBIL+/Tp4/06tVL1qxZI+nSpZOJEyfK7t27pU6dOgmqJYczw4cPl3LlyknVqlVD3RQSjZw+LSkPH5USx0XmfSVS6D8vFm5//vSU2p7nrEjmizfewhgdQggJg7WukFa+cuVKtYBn9uzZZcGCBVK+fHkVmIxFPmHdiSS41hUJFFNmDZfyrbsqUbM1u0j7liLfThbnY23p0eRIn0OOvnQ0lE0mhJCoG799tuhkzJjRGZeTP39+2bp1q/O5I0eOJLW9hEQ8T0x5Qhp800CuXLuiHu/PlkqJGYgaiJu/Rl+/PZQ3UwKRA46dPyaP//y4LN6zODRfgBBCohCfhc4dd9yhrDigWbNm0rNnT3nrrbekU6dO6jlCYhGIm1ErRsmc7XPUauRg8d7FSszAkmMy6vm6CUSOZvTK0VLv63pBaDEhhMQGPgsdZFVVq1ZN3X/99ddV4UDU1SlatKiMGjUqEG0kJOw5eu6oS1r5wTMHZczqMSomZ84s1yrIjw2ZK22yXF+x3I7ECgsSQggJkNC5evWqSi0vXLiw0401cuRItfYVgpJZNJDEKofPHXbexzpWm45uknzHr8jCb1JJul17RYoXlxqdrrux8h06K999vE/eLnU9C8uOQ2cPOe/P3T5XpacTQggJsNBJmTKlNG7cWI4fP56EjyIkejGFyZDFQ2Tpkskq26rw0StK5Mi8eRLftLWKzTkbn1dSbN8u3V6epOrp2HHbZ7fJuwvele3Ht0v9b+pLpZGVgvdlCCEkll1XqJuzbdu2wLSGkCgQOrDmvL78A1UnZ3+e9ErkSHy8fNvyW5n9ykbJuHCpEj9xefKoooGgSYkmkitDLqmQp4J6vPvUbnl59stSdnhZ536vOa4F/4sRQkisFQx88803VR0drHV16623KveVCVO0SSxy+OwN1xVAEUAUA2xfrLEMi49X29KkTCOlc5YWySki8+dLxsyZZdiOKVIsezGpEV9Dzl0+J/3m9pM1h9Y493Px6o1CO6cvnpas6dxEMRNCCEme0HnjjTdUhhUyrQCqI8fFxTmfRzkePEYcDyGxbNExxY6joGsgspNChdRN+0rtnZsypckk8VmuiyI7Tl08RaFDCCGBEjrIsHrmmWdk7ty5vn4GITEVjFy1QFVZum+pup89fXaf9lMoy3UBZAdWP48X90KIEEJIMoSOLqCMpR4IIfYWneHNhku6VOlk6ZT/hE46/wkdWHQIIYQEMEbHdFURQm6w48QOp1BJGZfSZVkHX4jP6t5ig/o8hBBCAih0SpcunajYOXbsmI9NICSygbUTmVbgppw3yelLp53P+eq6ypcpn9vnaNEhhJAACx3E6WABLULIDfae3qtWHoclp3j24nLgzAHnc9nSZfNpX6lSpPIYo0MIISSAQqdNmzaSJ08eHz+CkOhm45GN6hYiJ3XK1JInY/KukW3dt8mWY1tk4e6F8vr8153bu03rJl+v+lpmtp8pGVJnSHa7CSEkFvC6YCDjcwixZ+PR60Lnplw3qdu0qdKqVHFQMW9Fn/eHujqNSjSS1+q+5rL90tVLasHQmVtn+qXdhBASC6TwNeuKEGJv0UF8jmbP83tk3wv7fA5GtvJwhYcTbEuZ4kawMyGEED8JnWvXrtFtRYgni44hdFDYL3/m/Mne95iWY+TNem+6bGP2FSGEBHCtK0KIZ9eVP4HL2GoVwhIR7y18T05cOOH3zyOEEIn1ta4IITdYeWCls4aOadHxJ3cWuVOlnetsrncXvqtuj5w7Iu82un6fEEKIPbToEJJELly5IJU/razuZ0mbJdnZVu64Oc/Nsr/nfulZvafLdr3MBCGEEPdQ6BCSRObtmOe837Vq14BnJlpr8pTLXS6gn0cIIdEAhQ4hSeSXjb+o26eqPCVvNXgr4J+XNa1rsc6LVy4G/DMJISTSodAhJImsOLBC3TYo3iAon2e16Jy5fCYon0sIIZFMzAqd4cOHS7ly5aRq1aqhbgqJ4KUfQOGshUMidM5eOhuUzyWEkEgmZoVOly5dZP369bJ0KQM6ie/M3zFfdp3cpe4XzFwwKJ+J2jwmZy7RokMIIYkRs0KHkKQybs04qft1XXU/TuI8rjjuT6zrW1HoEEJI4lDoEOIDB88clK7Tujofp0+dXi3kGQzis8S7PKbQIYSQxKHQIcQHes/qLcfOH3M+Pnf5XNA+O2+mvLLkiSXy3f3fqccbjmyQCWsnBO3zCSEkEqHQIcRLEPz7w/ofQtqGqgWrSpX8VZyP20xsI/tP7w9pmwghJJyh0CHES6ZunqosOMWzF5fX677uLBQYbDKmzujyuObomnL03NGgt4MQQiIBCh1CvGTVgVXqtkmJJtL3zr4yr+M8ea/xe0FvR6Y0mVwebz+xXUavGB30dhBCSCTART0J8RIdm4M1rVKlSCV1itYJSTsypnG16ICDZw+GpC2EEBLu0KJDiJccu3Bd6ORInyOk7UiTMk2CbScunAhJWwghJNyh0CHER4tOqIWOXbwOhQ4hhNhDoUNIBAudW/Ldom4pdAghxB4KHUIiUOgMbjRYmpZsKi9Uf0E9Pn7heKibRAghYQmFDiERKHR61ugpvz3ym3P5CVp0CCHEHgodQhLB4XDIU788JacungoboWNd0ZxChxBC7KHQIcQNl65ekqd/eVr6ze0nny//PIG4CAeyp8vutDblGpRL3vrjrVA3iRBCwgrW0SHEDd+u+lY+W/5Zgu2ooRMumKLr6Pmj8r+5/5NXar8S0jYRQkg4QYsOIW7Yc2qPhDvpUqVLUFfnwJkDIWsPIYSEGxQ6hLjhwpULCawn/Wr3k3AiLi4ugStt8Z7FIWsPIYSEG+FjgyckzNh5cqfL42MvHVPCItxA5tWhs4ecj+fvmC8tyrQIaZsIISRcoEWHEC+FTjiKHFC7cG2Xx9O2TAtZWwghJNyg0CHEDZuPbnbeD2cLSa3CtVwebzy6UbYe2xqy9hBCSDhBoUOIDcOWDJPD5w6r+x/d9ZF82fxLCVfuLn23FMtWTG4veLs0KNZAbRu7Zmyom0UIIWEBhQ4hNoxaMUrdPlXlKelerXtY1c6xkilNJvm367+y6PFF0rFSR7Xt61Vfh7pZhBASFlDoEGLD3lN71W2X27tIJIAU8xRxKaR5mebq8dbjW52VnAkhJJah0CHEJq1cu60KZSkkkUSWtFmc1ict1gghJJah0CHEwr7T+9Rt+lTpnUssRBJanEVCwcNwYNa2WTJvx7xQN4MQEiAodAgx+GXjL9J5amenYAjXlHJPUOgkzskLJ+XspbOq/lCjbxtJva/rOVenT/jikyJ73BxLbMfzhJCwhUKHEMNl1eqHVvL71t/V47yZ8kokUijzdaHTaUonZalo8E0DtW4Xuc7FKxelzPAycvOIm10sOXO2z5Fl+5bJP/v+ufFiiJi77hKpU0cObFgqfWf3lV0nd11/bvdutV09T7FDSNjCysiE/MfBMwfl4tWLzsfHzx+XSKRgloLO+7BU6EG8faX2IWxV+LDt+DbnemAjl410bp+4YaKMXzte3T/X95ykTplaUp0+LXLokMi2bXK1Tm35tt0Fmb5lusyq96XE1asn2ff9d47gdVmzhuYLEUI8QosOITaLYaaMSym9a/aWSCRD6gyhbkJYY7r05u6Y67yvRQ645dNbJNegXLI942WRefNEiheXgocvyLyvRNItWSEXalVXImd3tjj5Z+z7IoUKuXVpQUz1mxtea6QREkvErNAZPny4lCtXTqpWrRrqppAw4eDZg+q2aoGqcrbv2Yi1gHSo1EEq5a2UYPs1xzWJZqZumioVR1R0dT39x+vzXpfaX9aW0xdP33A9eWDT0U1y8uJJmbltpkh8vBI7W7OLlDgu8tdokQKHz8uFlDimDnnox4dcV4z/z6V1qVED6Tq2vYr5GvDHAFl/eL2/vzIhxAtiVuh06dJF1q9fL0uXLg11U0gYua50bE7aVGklUsEinyufWSn33XSfy/Yj545INHPPuHtkzaE16hb8ufNP+XDxh2qR09fmvyZ/7vpTCRdvhI7m6V+fVq4qiJ32LV2fO5pBpMhJkRmjLsmSxROvC0mInLp1lavrwPbV8tOyMc7Xn7l0xn9flhDiNYzRITHPpauX1K2elefNGJlByFbis8S7PN5/er/kyZhHoh39O9b+qrbTDan5a/dfqpiitfbQ1WtX5ezls7b7a/pdU+mS7z75drLr9sspRHZmvW7lkdZdpVmbPjJpaibJsGu/OIoXl5rNt8leI2yHQoeQ0BCzFh1CAAa4SiMrSflPysve03udFpFooHDWwi6PXdwrUUjalGltrVdXHVed999f9L789O9PLu+7fPWySwC3lUInRZ7vM0UJGrivanS6flv0pIhDboid6SNOK5FzsmAuGT2kveyxxCZHanA70+tJpEOhQyTWiwP+e+Rf2XJsiyzesziqLDpWobP/zH6JZtKnTu+8n/u93B5fi+UyNOevnHepodO4RGP5oMkH6n7Bk6ICkLXIqfuoyKLC12+12ElpCX1q2uSIPLHy9QSf+dKslyTv4LzKpab5Y+cf4b3SvJFer9xyJkyvJxEChQ6JacwMnFUHV0V0/Rwr1uUr4LqK5to4Jy6ccPv8wxUedt4f3GiwrHh6hTNgu2nJpk73JZjRboY8d8dzMuHBCXI6rcihjDdEjrbSXMifSz3ekVUk1znXz4KLC1Yga/Yb0tpRoLD+N/Xl/OXzsuHwBqnzVR0p9XEp6T6tu9w37j4lusMqaNxIr1exR1rsGLFI6nm8jpAwhUKHxDR2gakV8lSQaMDqgtNZZdGCw+GQyRsmy+ajm+WbVd94XPC0f53+zsdP3fqUVMxbUX5u87P0qdVHvrjvC/m25bcqVufXtr86X9e6fGvJmjaLdLpPpM6jIufz5XQ+1/bmtpLvtEgKh0i6q64uLVh/5n8dJ/Pq2xdpvHLtivy88WdlzVHfQxzy8ZKP5ZdNv0jZ4WWl7ld1VdXmjUc2KrcaLD63f367TFw/0adjo4PrkwXS5v9Lr3eKnb/+uiFysB3PW9PrCQkj4hy4ImKYU6dOSdasWeXkyZOSJUuWUDeHBJn3Fr6nXAqarGmzyrHex1xcG5HMi7+/KF+s+EJZO2DV+O7+7yRaGL1itDw+5XGPrxlQb4DUKlxL6hatq1xG6VKlk6oFvSwpcfKkLKmQU3KeuXrdVfXeRRXMjGVC3srcQlLeWVtSXxM5mT+HvPFGA3myxRuS8cAxKXjfI5Ji+w4lAgpZApK9pUmJJjJj6wxldcJvt/PkTrXd0d+hLEPfr/te2tzcRr5Y/oW0KtdKKuVzLSfwxvw3pP+8/jLugXHqdcnGtOBotMhB+j0hYTx+U+hQ6MQ03X7rJsOWDnOJz4DrIprA8g8dfuqg7sOC8XaDtyVSgPWj+fjmUj53eRnUaJCKpUIWVK/qveT79d+ris/u+P7B76VV+VZJ//A9e2RrxXhnfE6JVbuuD+oY9GvXFtmxQ66lSilxCxZKXLVqCUVBnjySte5iOZXuxlNI+Z+ycUqSmzSp9SQZv268EjqaGvE1ZGGnhS6vi3v9+hptsFKdfNlP8TOw5NSseePxwoUiNWr4Z9+EBHD8jo5pKyFJZNepXQlm0lHDf9kyuTPeCMx9Z8E7yi0SKdkyM7fOlN82/ybv/fWe0wIHsfPM1GcSvBZp5Ka7Lj5rMi0NhQo5g45VCrnpttmxQ6RoUUmx8C9XkaM+OF5k/nyR6dOl2W2u1pSe1XsqlxgECIiTG4vGVitYTcbeP9Zjk+7//n4XkQNgZRq7xv59sEz6ZS4L8fbII67b2re/EbMTIecTiU0odEhMY2a8tKvYTrpX6y5RgZEtU/CEa3DrkY0rIiZbxsyGQsDxhasXnI9Na85NOW+SvS/sVa4ad3WEkgKCj7XYUW4bWDR0bMoff4jcfrv9GxGzkjWriv0Z2GCgczNcUQhyhpVl2ZPLZHaH2SoYGvSo1kMq5K2QIL7IHY/d8pjz/iOTHlGL0SKQefXB1c7tcHshAw2LkSYZw4IlqVKJTJrkGrPz998Rcz6R2IRCh8QsyLTZeHSjur+jxw41KKVKESU1NI1smbKtOqssIIDb1A0bh022DMTLhLUT5PDZw7bPn750o33IUDJXG9dkS5dNxj4wVmXLFclaxLndH/WQcE5A7GwY0sfyxLdexabgfDLFS9Z0NwJ2bi1wq9QrVk/GPzhe5j86X8XSQLBpimUrJlu7b5U5HRK65xBMPbr5aJdtTcY0kRqjaqi6UCZHzx9VljwEM5vWHQQ6J5rhBUuNtmBB5Fy5ItKrl8jYsTfETq1aYXM+EWIHhQ6JWZDKixgQBCBba85EPEa2TKodu1QtmOq7rteEKXDovLJQrJnwccizZd7+821pM7GNNB7T2G2dI83fe/9OkCVXMkdJOd77uFTJX0U9Lpa9mHzd4muZ0maKpExxoyJyUoGV79Qj6+SeARPcu20S4a6Sd0nfWn1VlpcdcGPVLlJb4uLi1Irpb9R9Q5rf1FxZe1AiAGKobK6yLlaeDhWvx1z1q+26WCiOkTse/OFBafhtQyV2kNqe6Z1MKljdI5kzq1gjJWoWLLghbh5+WGTw4Bvip2hRZl+RsIXByAxGjlnGrB4j7Se3lzsL3yl/PHY91TfqsMmW0TVhBncaLw/d/FBIm4caMoi50RlFVh7/+XEZvdLVcmEyqOEgebFmIoO1v44fBnlYciByzNTqIGQdIdPqh3U/SI87eqgK0BBFAEIdcUzNxjbzel8IXP70n0+dKfnX+l1z7s8WuKNOn5Y/r26X/t92kqmfn5X0u4yaTBA5cOMx+4oEGQYjE5IIejXpm/PcLFELBh8MzgZYnBLumHBY5NNch8oOvSyHSf5M+eWLe7+Qp299WrpV6xa4xmm3jSlqkGVkrSvjbnkEP1I8e3HpXau3So83RQlcY01LNZXxD4yXyQ9NVvE/OdLnSPB+PN++Ynt1v+bomi51h2AFaj6yrkyY/r79h8MdlTmzdJrSSeZe3SINGloKT373nTrPhiwaIm/98Vai7uJOP3eSr1Z+5Qyk3n3SO8sYIUklSgISCPENzIT1kgjWCsJRBSwSsEBYKvfCooMBruOljpIpTaaQNS8x95Kd0Lkl3y3yeJXH1V9A0W4bYFpucIvH/6WQq9eFGNMyh/pCqMFjggw0uOG+XZ2wiGGT4dVl+hiRPGfnS/eBC+Qv2S2v3PmKcqel2rtfst7VXBx58six+tul0Pnr548L7dvLgV/GS8/fe6qH+Jwi2W7ESpl8/s/n8uXKL9UfCnNCdKVOkVouvXqjMjUh/oYWHRJzIAizyqdVnLPK3Bk8r4sUsVjcLoPfa+lMlUasztw/v5UHv39QuT5mbZsVcouODoz9csWXKnbk9MXTqupxyCpXZ82qUsRVqrjVLWOkkKvXhRHDmg2TUjlKuWxDBhqEix2ZL0LkXD8verz8kxz89x957OfHpNHbZeRE9crq/LlyYJ8UP3TVue7XiQI5r9fR+c+ylbHJ3c6A901HN9l+zoilI6TrtK7Ox9O3TFe3l69dtl1sF4uvnrzALC6SfCh0SEyAULRzl8+pdNul+5bKmkNrnM+ZdWaiBhu3y/mqlV3qwmDQWrt8hgoEbvRtI7X+UigtOvhtYGmDi2TwosFScWRFuXj1YoL3BPX3gohxF2D7Xwp5uIEA7U3dNsnjlW9YvPJnzq/cXnZZhajcbD0vym0+KT98ckSKHL0il4vGyxe9Gsj4iTcWN33ihRJyqdptcmXOLCV6Mu85rN6HRVBfnv2yEtDm2mO4/p797VmXzz1w5oBL9p3JsCXDpOWElqpukAbp858u+9Q/dYFITEHXFYkJOv7U0dZsH7UWHRu3S488PWTt4bUyqORpeel/09RilVi00lwLq2i2ov5tx3+BrLZiwVJkLuegG2tJgR0ndtjuEgX3SOKMuHuEOrfhRtICB5lcfWb3kcbFG0u1QtWk27RuKhg8Lj5e6j6622mx+eu/+G+ImmWfdJfX/3pVbsl4PZW/bpsTsufMEkn75vWTp9BD18WRPp/27l8uy/cvV8tudKnaRcUUIcPRilmR/PC5wy4u5JH/jHTWSoKwwT6QPg9yZsgpD5Z7UIJCYucvrrMwFLvEFWZdMesqJtAl8e3Y2HWjlM5ZWqION500ZsadhjdRg5K5PMGSJ5Z4vw6Ut5+PInKor2LNTvrPrbb0yk5p+PBVl3ZYKZC5gDPNfPfzu6M7pirIIM187aG1ck/pe1TqearFS2T+51ecz2Oh0kX/VV5omKOqTL1vvGT+tqzLau8Alhzr+QSwdEfvmr1V1thr819z247lTy2XyvkrOx9X/byqLNu3TN1f/cxqKZu7rKQekFo9frBgI7knXx1ZcG2H9KrRS0rkKHHDUuWD+Fiwa4FaD+ydBu/I7QVvT9L5qyYTYei+jBVOMeuKkOu+/o8Wf+TxNVFp0fHgdsH3hbvCOihhVh2oooVqUNB1Z4zYoRynr6oYEXese3adyq7S7abI8S8QEFgPLH3q9DK/wRiZM8v1+CLwWMfevHLvIElTpLgcfjHheaLPJ3NJC7Du8Dq1zpoWOWNajrEVFea5h7m3GZv1564/nRWys1wQefGNOVKrw/9k+rwvpMzwMvLc9OdunFceKjQjBgwrwuP1cJu9+cebymJU7YtqKh4sKedvMIskDlwwUDpM7pB4kUeSAAodEtUgu+O5Gf91hEZV2QbFGjgfwxwfS5jpx0PvGqpiOoC76sT+KFroHCz0WlHbtsnVYkVVbIhe3RuVgZuVulEPpl7R64XyXqr5koy8e6QsfXKpf9tHbrB7t6Sq31BS/rfqOgKNTxXM7YzZqZ+qlNQpUsdZ4DBj6oy2u/nnqX9kU1f7YGQsf9G2QlvbkgLmuTf538ly8uINoTJtyzSVPAAgirGavG4XRNjwpcMTFR9PTHlC4j+IV6Lmo78/kuqjqsuSvUuczy/cvVCOnjuqrFvuzl9H3bqyf8ZEdetSciAIRRIh/uByhPt99rbZAf+8aIMxOiTy8eBHX/3Pb2oWaFovsHiiaXr3WCwtCoErKG/GvCrwt3PVzirNHHEafrfoWFOx9VpRoHhx2fzDSNnzy/WKyG/Vf0teqP6CssChPbUK11LWAfw2CKJ9+rbrVh0SpHpB8fFyYtpkOVq3tpQ4dk0mjzwhcR33Oq8x1Ox5eNLDMvq+0SqgHCvMw1qj3U8fNvlQWWFqFq7pjK1B2jkWGT1+ASukXgep5ci60ucerkuIEpNfN/2aIHBaxxLhFnWhdPtPF8otg/vXkf4FC8jVq5fli+VfqKrno1aM8hj/hazDR396VMWpbeiyQS2/kTZVWnUcrs6ZLYduLy/5t22T/Hf9FxuUjGKRyPqcv3O+WnU+Q+oMXr3nzKUzzvuHzh7y+TNjHQodEtkk4kd/5dWZ0ja1yF3tboidzGkyq4DGWAXLDCArB4MOYhvyZMyTaAf6zK/PqHWnMJj5LAx10UItcrCw6KcfSIu/rhf7q16ouvS988aik/WL1ff9SxG/1wsqXKGmXPpnk1yt31Cy5M3nUi+oUYlGLi6sy6+6poijgjMwM/nuL3t/goVa4ZZEUDLcMnAf/bblNyWEIAB+eugn26VB9EKrroHT2+RSkXgpd/9u2bP9S6m6+X4ZtHCQcnt5IybeX3SjWGKL8S1k6/Gt8mOrH6V5mebyx7Xt8so955wB2r6sdeZu2RO48lqVayXft3Jdid4dEGCagExIohy6rkhkk4gfPe/BM6pGiBkHkjltZnnujudUpdiJrSdKLAL3gy4UqGOU3HWgZy+dVUsGjF0z1rlcQ3KLFl5p97Cc3bZR1XcZc/+YpHwFEoR6QWmKlpCUf/yZ5IBbxP7AnfXnY3+q6s7a2gNeqvGSU2RjIO83r58s3rNYPb7vpvvUoqfugNhRlhyDug12q+2g629dlcgBKCvhC1joF9bOcWvHqcdz/vjatkiit2uduRNVP6z/IdHXTt00VdWVMich7rIRiXsodEhkk0gcyPYcKVziQLRFB7O8b1p+45xlxjK6Lo27GB1zqQhzNu4VRuzE7lxpVBYPUpbzHTyrZuRfVO7vHABJCAlgvSAsuApXpAZxOtt7bJd3Gr4jLcu2tC1p0LpcaxVL5i5RALE5VvFhBk7vPLlT/FLjafdueaLnd876QTh/d+ZMlWBitf34dvlx/Y/OQGGIq0cmPeKy1IY50fKGKRunyD3j7lF1pbQAdPfdsNgt4o+wQj1qEqEdiOuBYDMZv3a8FP2wqEt8UixAoUMiHx0HosUOXCTbtsm5+PxSu+M15yzP144mVtCz6j2nEq7ZBHeCadbXy2b4GvtxqUghqdH+kkpVNovT1X3s9aCsFUXCC4gbuE6xzty27tuUBUcDS0+LMi3U/Q+afJDgvRAz2m2FiYwWz2aAMni4wsOqlo8nzCw+66THsXu3XK5dSxVN1Avh4vx9/PniciY+rzqvz9Ssqs7fp399Wlr90Eqen/68eu+kDZOUBfTFmS+qFH4TM5DbrO6CYGgzdnDksuu1hHSwtCeLzmvzXlPiBSvUF/2oqDT8pqE89ONDKgD7iDFRaTuxrRJKL896WWIJCh0SHaCGwvuuixI2bLRfiZzbruaTT2q+49weyrWdwhHMuAGqRZsl91Hkre7XdeXjJR87t+l6Nj7FfhQvLs/3qewUnDrGAoNHqrwFwmKtKBI6EPNVOd+NGjqPVHzEGQeG+9Mfub5UBKidorhT5BzMm0kaPJbCKZ4vFY13ih3U9WlTvo0KKrbD0d8hK59eKfMfne/c1rRkU5fX7Lp2XE5kSe0UOfr8/TvFPin74EG1fY3joGy9fEhmbpupnhu6ZKi6hv7ced1tBpdTuU/KuWRKZUyTMYGF9I35b0iewXnknrH32MblzNsxz3kf7mNrirmZyYbU+bk75irXGO6P/E8wmfFScCnGEhQ6JDoCkuvXF0erVglM2VX3iPw88oS07vWlyr7SriviOqtFijk6T8Q1jFo+SlK8nkLKDi+rKtya+CR0/ov92DzpC/nkwC9qBo/AY4BBo8NzhSXFjBkstkYkZ/obyQHWdbqalGyiyiCgDMSHrUapCswQGT16V5Ttma84z6c0fyyUi0UKOSs01ylax1boIAsLVMpXSblNx94/Vn5t+6uUyF7C5XXbHcfk3f/VkzqPirRs1E3a3tzWmQGFz8N2JDmU/MY1lggp6gt2L3DZ9uafbyprDSw4iHnTIHMNa4ChcCGuPwgmbdUxXcmmyxifD1eZiafaOhuPblS3yGbU2K1wH81Q6JDgChJ3bgrLcgA+sXGjyOrVEnflilxOIdKi9Q1T9sLRIgUOX5D0x047A5LpukoIatYAmLSf+OUJcYh9wXRPQmfNwTXOWALEKKBDXn5uqwzaOVZtQ6E400WRuXhZihyiuLPInR6tDd2qdZOjLx2VymXqKnEBkYFAaZ1RdUehO5QLO82ff8mfI/vKh61Hq4D7m3Ld5NxHh0odZHiz4WoZDBPEDN1d+m6nC1cDobH24m4V34faW1/c94XL83ZFNwFWcV9/eL3LNlhksGQG0tyPnj/qUqXcugYYihoC0+VkZeWBleq9uObsEgkwqYNVS78W6DZh+5QlYxKtx4OYI8QJgdUHVydwwXnDqgOrZPBfgxPECgUbppeT4BCocuoQSG3bily5IldTxknqqw55f6ZIz8YiP/wgkvqayNWUKWTL8Ddk7/In1Vto0UlIp8qd5PPln6tKtp4whQ5M9FnTXf+tUP/mjlF3KIGz4LEFavXrzcdcVx6H1ahg5oLOx7cVuM3v34NEJrfku0UJEG1tsQMWQfDLM/NVkUAs3fBElSdkyKIhMrTpUPUc1ux6Kf4t53ua39TceT9T6kzybFVXUWEC6076VOnlquOqsqqgvs+MrTOcz0FU4S+xLC7TcmLlyV+u90GeQGHAmvE15fyVhIvsotDp7O2zZcK6Cc6srWv9rrlkZUHkTB9zfUV6uNzWx61X3wfXro5vgtXrLmkoR968pMpNgAtXLqiFWPNlyqdcXog50nFSz8+4Hnv0VJWn1OthddMlBHSto5///Vn9DqZQveXTW5x9bihrYdGiQ4LD6dNy+cA+v5ZTx8X77qpP5Ez2jCoO5KX+NZ2WnJ8mXBc5sPCcKBUvKcuUc76PFp2EYEZ8b+l7nY8XdlooT1ZJ2CnvPrVb/tn3j/SY1kOyvZtNukztoszx209sdw4Atb6slUDkgJLZS0quDLmcj1FAjhCzfpKu0u2J2kVqy4QHJ0jBLAXV/Z/a/ORWIKHo37yO81T8Tc8aPT3uFwM01lI71CthPSmdGWiev6j07A/SpkyrRL+uCj5181TpO+dGXSmzVpEO0jZT0+GaMi06sFxD5Oh4pQInrsnHf38s57dtcsY36ZIbOpX+g0UfSJlhZaTIh0WU9QkWHI0WOeCz5Z8pkYlq8wiexjUPoXnvuHvlixVfqDIUdoHWf+35S0IJLTrEFr1i8M4TOyV7+uzKDJwcJp1eIj1a7pLl47JIbi12UHQL9Sh8KKe+99Re5a9Gtsarc16VQcsGyZCmaeXg01vkj2ktZVHLGysvg+/7t5I23UbKKbnhFqNFx56xD4yVnjN6qkrEiKWBKf+7Nd+5zGA3Hd0kt31+wxLzybJPVCwE3pMYWHyxenx1ZdXBytllcpUJ2HchRIPzE3/e4K6QqM7OQlVxpHJrd++qg6tcXlcud7kEbiuA/vPUxVMu22A5QbzO1u5bVYkHWER+2/yby2twrbxe93UplbOUKu6pXcwmny77VC2aCt5t+K4s3bdU6sqPrtWjT/ZSMYs6TV6X3IBrCZ/Ze1Zv5/7qfV1PnqjsWp3aDsQivfD7Cy5xfMgOQxFItMEaRxRKuHo5Vy9PwLTN06TlhJbyv9r/U0FyyMpJbJ2h71Z/Jx8s/kB+bP2jbV0MBLYiiwem090//5cGrnFTTh3mU8yg9MrEML1i1oEKvd1v7y6T/p3kTImG+bZa/wIybtgBdTFb930ge2rJ/35+tenkyyeTLdxiBaSXZxno+VihlP19pe+Tl2d7TlnFzBoDDroc/NOuCELCjV6/93Ipq4AsLV2HBina4Oc2P6ulL6znN6whWFICFksU+2tfqb0K/n193usycOFA5+thPYK7yLRiYSJx07CbXFx6K55ecaMdDocUGFJA9Y12HO99XE3kUg1I5ZKGr7FmkCWHPrX6yDsLbmSzapBBt+LAjTYDWN7MDDd/wdXLSZKBefLi1Yvy6txXVTT/sn3L5OCZG6mOdrSb3E7+2f+PPDvV3geuZ/zqAoMlRzyXU/9r919KmOi6FGDo30OVyAHwnZt1XyZMG+wUOVeKFlaLEppFBPMdvyydb+usBBJFjvfAzfdq7VeV4LSuTG3+VhCd7kC8AWbEutItLIUUOSScGdx4sKztvFbF7HS7/fpSJeCh8g/J+43fV0HNCFDWYCFTbTXCdqzbBoto71q9lRUI/R8KJJ7te1aeufUZGffAOGe2o0npnKWlY6WOzsfWgom4dhqXSLgshiZr2qzXCx26qR7dvqV/RA74ZdMvttutIgckqaK6H2FvQ1yAv3bu9rkJtus6EYmBNWLc+aEBZhkH73e9UC8+/JC88lUHlyyDYUuGXb9dOkzNejCT0SXZzZRJnUVQtd1LTrPstblzRGrUSFAx+ZNKfeWjph959T3IDd6o94aKW2hYvKHL9kdveVTuKnmXum9XafX5O55XHfrcjnNVUTjWLyKRRPk85ZWF5KO7PnIRGhAxCGo2A+vjs3q37hWCmUfcM0La3NzG7WuwLIq1arnJg2X/W1jUQqPijZz1h5Y8sUReKNhKfpvu6oozq0cDnTKv6VW9lxTJWsSr7+Ky0nsiwBpvptUHGwod4uTvPX+rglUnLyZM824/ub1KPbbzdJrbrH5oxPhAPGEGr02peQ+eVYLEsWCBEiJpd+6RTi98K/2/7eR8n3ZXgfRvpZcao2s4feNWUDND19a4v3N2lXaaoGIyMrpYmC7JoAM164ygQ/yy+Zdyfxn3S2gMaTJEdejI0tCZHYREEghmdreILc7pfS/sk70v7PUqRs1bzGrN1QpWS/C8adGBUIGl9Ze2v8iMdtczxEDVa/nk/Xf+kWz7jsrZ+Hwy/tNuCapHD2k8RMXloSAj3ovr+c36b8qHd32ogq+nPTJNXedWPm56o4CoN8CihZg/n5eP8SOM0WGMjgJiJPWAxAejwY0Gqwq6w5oNUzP0B75/QLYe2+oSlIfsANSgQGZAj+k95MUaL8qiRT/IVx/ucAmGK1u5kbxc5BEp0uJR5/a8SzdIphJlpObomsolYsW60jYuetSTQEolsgiyl7pZ1nS+XlvCJQUdIoc1W5IF1tFBiXkw4u4R8sxtz6jfGRVdNegcsfozZsGhTCclJFKZsHaCtJl43eKzudtm20w0hBPsP71f1f9BHJ0u8+Ds7+rUcUnyuFqwgAwc10W69p4oWfcekV05U0nWxSsla8nyHtuCdbNQGwvZVNqifuXVK5J1YFY5e/m6hQaubVjd3/vrPZeJat9afVVsToPiDSTU4zezrgJE6Y9Lh1TB+gpqR3hDr5nXFT4uvu7Vuqs1Xaw0GdNEReKjJgPABVDIkVlZXYAOhtuzbaZyiRV69EZth77znpPR8RPd+nTL5iorf+z8w/l4UutJsvf0Xik3vJycSndV6hh+cyeJZHIR78BaQAhAhIC9u9TdtqZ1uLJO9zlNCw4hSQTB/cAuhseuBpWLyDGXXgH/JXmkFJFXHhkpUvsV5caPz5NH4nIX8sqihdijaVumObchBghVpfVEFGn2yKLUQqdCngrKZQ3XX7hAoRMgIHLMCpiRWJLdU/vhonJXKXPO9jkJtu2JO60qmsLqYq4krp7LKtK3Xw2Zc2SpHNkzQ+KmPO5SAMukWPZiKg4IgcjZ02VX68YggA/FxrDgpFkgjPgXmPDfbvC2+jNBTRMsIPhhkw/VY4ocQpIO4n22dNuilrxIEv8tvaJqklknefHxIvPnS5yPFu6Xa70sG45skEcqPKIeNyzWUAkdlIhoWqqpisGEaxsT5uVPL3cJPQgHYtZ1NXz4cPV39epV2bRpk99dV0gTRLXYSAJ+ZpgpUTcFcTrIuqpaoKrM35kwLRDl/O0CUJMK/MMwf3ae2tnj635o9YOqP9Fvbj/lPkPdChIewh7C0108AyEkejh3+ZxKWsHq8HqRUr1oaDAXDPXWdRWzQkfDGB3v6upgjRZUInaXUghQxA/p6HYFswBSNe3KmoM/Hv1DmWlRI0IzrOkw6Tqtq7oPf2/X27tK/szXa+EQQgiJbU6xjg7xFzBNohDglLZTZPJDk92+DiLHrkYNKn8efvGwHOzlvhYPRA4sNSatyrdyCUKmyCGEEOIrFDrEJ/JnchUbSBvUIIDYrJWQMXVGZdpEqiIKzqH4HBaIs4LXQeSYbg8IH3NdGbpECCGEJAUKHeITplUFtXFQDG7qw1MldYrUSsSgBgMY1HCQ7Hlhj/ze7neX9z93x3Nytd9V+bXtry7xPlrITGw9UcXdIJvKrJ4LtxghhBDiK+EVGk3CHtO9hFLnSD/Eqrtn+p6RNCnTqO1n+pxRFUDdWWEgYFD/QdOgWAOXFGb8aRY/vlhV4LRbzI4QQghJDAod4hMQM09WeVJllX3b8saaVVrkAB2Fnxg9q/dUmVuox+MO1GfAHyGEEJIUmHXFrCtCCCEk4mDWFSGEEEJiHgodQgghhEQtFDqEEEIIiVoodAghhBAStVDoEEIIISRqodAhhBBCSNRCoUMIIYSQqIVChxBCCCFRC4UOIYQQQqIWCh1CCCGERC0UOoQQQgiJWih0CCGEEBK1UOgQQgghJGqh0CGEEEJI1JJKYhyHw+Fc7p0QQgghkYEet/U47o6YFzqnT59Wt/Hx8aFuCiGEEEKSMI5nzZrV7fNxjsSkUJRz7do12bdvn2TOnFni4uL8qjQhnnbv3i1ZsmTx236jER4r3+Dx8h4eK+/hsfINHq/QHyvIF4icAgUKSIoU7iNxYt6ig4NTqFChgO0fPyovAu/gsfINHi/v4bHyHh4r3+DxCu2x8mTJ0TAYmRBCCCFRC4UOIYQQQqIWCp0AkTZtWunfv7+6JZ7hsfINHi/v4bHyHh4r3+DxipxjFfPByIQQQgiJXmjRIYQQQkjUQqFDCCGEkKiFQocQQgghUQuFDiGEEEKiFgqdADF8+HApWrSopEuXTqpVqyZLliyRWOOPP/6Qe++9V1WtRNXpn376yeV5xMH369dP8ufPL+nTp5eGDRvK5s2bXV5z7NgxeeSRR1SRqWzZssnjjz8uZ86ckWjjnXfekapVq6oK3Xny5JEWLVrIxo0bXV5z4cIF6dKli+TMmVMyZcokDzzwgBw8eNDlNbt27ZK7775bMmTIoPbz4osvypUrVySaGDFihFSsWNFZfKx69eoybdo05/M8Tu4ZOHCguhafe+455zYer+u89tpr6tiYf2XKlHE+z+OUkL1790q7du3UMUEfXqFCBVm2bFn49fHIuiL+Zfz48Y40adI4Ro8e7Vi3bp3jySefdGTLls1x8OBBRyzx22+/OV555RXHpEmTkNnnmDx5ssvzAwcOdGTNmtXx008/OVatWuW47777HMWKFXOcP3/e+Zq77rrLUalSJcfixYsdf/75p6NkyZKOtm3bOqKNJk2aOL788kvH2rVrHStXrnQ0a9bMUbhwYceZM2ecr3nmmWcc8fHxjtmzZzuWLVvmuOOOOxw1atRwPn/lyhXHzTff7GjYsKFjxYoV6vjnypXL0adPH0c0MWXKFMfUqVMdmzZtcmzcuNHRt29fR+rUqdWxAzxO9ixZssRRtGhRR8WKFR09evRwbufxuk7//v0d5cuXd+zfv9/5d/jwYefzPE6uHDt2zFGkSBHHo48+6vj7778d27Ztc8yYMcOxZcuWsOvjKXQCwO233+7o0qWL8/HVq1cdBQoUcLzzzjuOWMUqdK5du+bIly+f47333nNuO3HihCNt2rSOcePGqcfr169X71u6dKnzNdOmTXPExcU59u7d64hmDh06pL77/PnznccGg/kPP/zgfM2GDRvUaxYtWqQeo2NNkSKF48CBA87XjBgxwpElSxbHxYsXHdFM9uzZHV988QWPkxtOnz7tKFWqlGPmzJmOOnXqOIUOj5er0MGAawePU0J69+7tqFWrlsMd4dTH03XlZy5duiT//POPMtGZ62nh8aJFi0LatnBi+/btcuDAAZfjhDVL4ObTxwm3MGXedtttztfg9Tief//9t0QzJ0+eVLc5cuRQtzinLl++7HK8YFYvXLiwy/GC6Thv3rzO1zRp0kQtqLdu3TqJRq5evSrjx4+Xs2fPKhcWj5M9cLnApWIeF8Dj5QrcKnC1Fy9eXLlT4IoCPE4JmTJliuqbW7Vqpdx0lStXls8//zws+3gKHT9z5MgR1fmaJzvAY/zo5Dr6WHg6TrjFBWSSKlUqNfhH87G8du2aiqGoWbOm3HzzzWobvm+aNGlUp+DpeNkdT/1cNLFmzRoVJ4FKq88884xMnjxZypUrx+NkA4Tg8uXLVRyYFR6vG2AA/uqrr2T69OkqDgwD9Z133qlWx+ZxSsi2bdvUcSpVqpTMmDFDOnfuLN27d5evv/467Pr4mF+9nJBwnH2vXbtWFixYEOqmhC033XSTrFy5Ulm+fvzxR+nYsaPMnz8/1M0KO3bv3i09evSQmTNnqsQI4p6mTZs67yPYHcKnSJEi8v3336tAWpJwQgZLzNtvv60ew6KDfmvkyJHqegwnaNHxM7ly5ZKUKVMmiMbH43z58oWsXeGGPhaejhNuDx065PI8MhgQpR+tx7Jr167y66+/yty5c6VQoULO7fi+cIueOHHC4/GyO576uWgCs+uSJUvKrbfeqiwVlSpVko8++ojHyQJcLriGqlSpombK+IMgHDp0qLqP2TWPlz2w3pQuXVq2bNnC88oGZFLBimpStmxZp7svnPp4Cp0AdMDofGfPnu2ifPEYMQTkOsWKFVMnsnmc4MuGX1YfJ9yiY0FnrZkzZ446nphtRROI14bIgQsG3xHHxwTnVOrUqV2OF9LP0amYxwsuHbPjwEweaZvWDinawDlx8eJFHicLDRo0UN8V1i/9h1k44k/0fR4ve5DivHXrVjWg87xKCFzr1hIYmzZtUlawsOvj/RbWTFzSyxFZ/tVXX6mo8qeeekqll5vR+LEAMj2QZok/nGpDhgxR93fu3OlMPcRx+fnnnx2rV692NG/e3Db1sHLlyip9ccGCBSpzJBrTyzt37qzSMOfNm+eS3nru3DmX9FaknM+ZM0elt1avXl39WdNbGzdurFLUp0+f7sidO3fUpbe+/PLLKhtt+/bt6rzBY2Rp/P777+p5HifPmFlXgMfrOj179lTXH86rhQsXqjRxpIcjAxLwOCUsV5AqVSrHW2+95di8ebPju+++c2TIkMExZswY52vCpY+n0AkQH3/8sbooUE8H6eaoERBrzJ07Vwkc61/Hjh2d6YevvvqqI2/evEoYNmjQQNVFMTl69Kg66TNlyqTSNB977DEloKINu+OEP9TW0aBzePbZZ1UqNTqUli1bKjFksmPHDkfTpk0d6dOnV500Ou/Lly87oolOnTqp+h24tjCQ4LzRIgfwOPkmdHi8rvPQQw858ufPr86rggULqsdmTRgep4T88ssvStyh/y5Tpozjs88+c3k+XPr4OPznP/sQIYQQQkj4wBgdQgghhEQtFDqEEEIIiVoodAghhBAStVDoEEIIISRqodAhhBBCSNRCoUMIIYSQqIVChxBCCCFRC4UOIYQQQqIWCh1CSJJ49NFHpUWLFkH/3K+++kri4uLU33PPPSfRenzwfv09f/rpJ7+2jZBYIlWoG0AICT8wuHqif//+arXwUBVWx0KJWFAwY8aMEmp27NihFjBcsWKF3HLLLX7bL47vwIED1aKShJCkQ6FDCEnA/v37nfcnTJgg/fr1c1mpOFOmTOovlEIMKyNHM1mzZlV/hJDkQdcVISQBEBH6D4OtFhb6DyLH6pqpW7eudOvWTbmTsmfPLnnz5pXPP/9czp49K4899phkzpxZSpYsKdOmTXP5rLVr10rTpk3VPvGe9u3by5EjR3xuc9GiReXNN9+UDh06qH0VKVJEpkyZIocPH5bmzZurbRUrVpRly5a5vG/ixIlSvnx5SZs2rdrH+++/n2C/b7/9tnTq1El9h8KFC8tnn33mfB7WHFC5cmV1nHAcTAYPHqysMjlz5pQuXbrI5cuXnc998sknUqpUKUmXLp367g8++KDP35sQ4hkKHUKI3/j6668lV65csmTJEiV6OnfuLK1atZIaNWrI8uXLpXHjxkrInDt3Tr3+xIkTUr9+fSUSIECmT58uBw8elNatWyfp8z/44AOpWbOmciPdfffd6rMgfNq1a6c+v0SJEuqxdrn9888/6rPatGkja9askddee01effVVFQdkAvFz2223qf0+++yz6ntpCxe+K5g1a5ayhE2aNMn5vrlz58rWrVvVLY4N9qv3je/bvXt3eeONN9S+8N1r166dxCNPCHGLX9dCJ4REHV9++aUja9asCbZ37NjR0bx5c+fjOnXqOGrVquV8fOXKFUfGjBkd7du3d27bv38/FIZj0aJF6vGAAQMcjRs3dtnv7t271Ws2btzoU3uKFCniaNeuXYLPevXVV53b8LnYhufAww8/7GjUqJHLfl588UVHuXLl3O732rVrjjx58jhGjBihHm/fvl3tc8WKFQmOD96L46Bp1aqV46GHHlL3J06c6MiSJYvj1KlTDk9g35MnT/b4GkKIe2jRIYT4DbiGNClTplTumgoVKji3wT0DDh06pG5XrVqlrB065gd/ZcqUUc/BEpKcz9ef5enzN2zYoCxAJni8efNmuXr1qu1+tRtP78MTcInhOGjgwtLva9SokXKvFS9eXFmevvvuO6elixDiPyh0CCF+I3Xq1C6PIQrMbTqb69q1a+r2zJkzcu+998rKlStd/iA0kuLGsfssT5+fnO/lzT48vQ/xPnCnjRs3TgkgBHxXqlRJufMIIf6DWVeEkJBRpUoVFQyMgN9UqYLfHZUtW1YWLlzosg2PS5cu7WKJ8USaNGnUrWkB8hZ854YNG6o/pOxny5ZN5syZI/fff7/P+yKE2EOLDiEkZCAL6dixY9K2bVtZunSpclfNmDFDZWklRTj4Ss+ePWX27NkyYMAA2bRpkwoYHjZsmPTq1cvrfeTJk0fSp0/vDKQ+efKkV+/79ddfZejQocqCtXPnTvnmm2+Uteemm25KxjcihFih0CGEhIwCBQooCwpEDTKyEE+D9HRYNlKkSBEUi9L3338v48ePl5tvvlm5j5AFhdR5X6wyECyffvqp+j5IZfcGfEdkaCHrDJalkSNHKjcW4noIIf4jDhHJftwfIYQEFKRnQwzFSiwL4nomT54ckuU2CIkGaNEhhEQccA8hQ6t3794SrTzzzDMhrT5NSLRAiw4hJKI4ffq0ioXR7h8UKIxGkIZ+6tQpdR9ZWeGwrhchkQiFDiGEEEKiFrquCCGEEBK1UOgQQgghJGqh0CGEEEJI1EKhQwghhJCohUKHEEIIIVELhQ4hhBBCohYKHUIIIYRELRQ6hBBCCJFo5f8MXDG2Mn8vjwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "### plot input series, output series, and observations\n", "\n", "# create figure\n", "fig, ax = plt.subplots(1, 1)\n", "# plot input series\n", "ax.plot(\n", " timesteps,\n", " input_series,\n", " label=\"Tracer Input\",\n", " c=\"grey\"\n", ")\n", "# plot output series\n", "ax.plot(\n", " timesteps,\n", " output_series,\n", " label=\"Tracer Output\",\n", " c=\"green\"\n", ")\n", "# plot observations\n", "ax.scatter(\n", " obs_timesteps,\n", " obs_values,\n", " label=\"Observations\",\n", " color=\"red\",\n", " marker=\"x\",\n", " zorder=10\n", ")\n", "ax.set_xlabel(\"Time [months]\")\n", "ax.set_ylabel(\"Tracer concentration [-]\")\n", "ax.set_yscale(\"log\")\n", "ax.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "7c66fbbd", "metadata": {}, "source": [ "## 2. Store Values" ] }, { "cell_type": "code", "execution_count": 6, "id": "80ed0487", "metadata": {}, "outputs": [], "source": [ "# generate dummy monthly time stamps\n", "start = \"1960-01\"\n", "\n", "start_date = datetime.strptime(start, \"%Y-%m\")\n", "out = []\n", "for i in range(n_years * 12):\n", " # calculate year and month offset\n", " year = start_date.year + (start_date.month - 1 + i) // 12\n", " month = (start_date.month - 1 + i) % 12 + 1\n", " out.append(f\"{year:04d}-{month:02d}\")\n", "\n", "timestamps= np.array(out, dtype=str)" ] }, { "cell_type": "code", "execution_count": 7, "id": "fd5c4fba", "metadata": {}, "outputs": [], "source": [ "# concatenate timestamps and input series\n", "input_ = np.concatenate(\n", " (timestamps.reshape(-1, 1),\n", " input_series.reshape(-1, 1)),\n", " axis=1,\n", " dtype=object\n", ")\n", "# store input series as CSV\n", "np.savetxt(\n", " \"example_input_series_1tracer.csv\",\n", " input_,\n", " delimiter=\",\",\n", " header=\"Date, Value\",\n", " fmt=[\"%s\", \"%1.3f\"]\n", ")\n", "\n", "# concatenate timestamps and observation series\n", "obs_ = np.concatenate(\n", " (timestamps.reshape(-1, 1),\n", " obs_series.reshape(-1, 1)),\n", " axis=1,\n", " dtype=object\n", ")\n", "# store observation series as CSV\n", "np.savetxt(\n", " \"example_observation_series_1tracer.csv\",\n", " obs_,\n", " delimiter=\",\",\n", " header=\"Date, Value\",\n", " fmt=[\"%s\", \"%1.3f\"]\n", ")\n", "\n", "# concatenate timestamps and full output series\n", "output_ = np.concatenate(\n", " (timestamps.reshape(-1, 1),\n", " output_series.reshape(-1, 1)),\n", " axis=1,\n", " dtype=object\n", ")\n", "# store observation series as CSV\n", "np.savetxt(\n", " \"example_output_series_1tracer.csv\",\n", " output_,\n", " delimiter=\",\",\n", " header=\"Date, Value\",\n", " fmt=[\"%s\", \"%1.3f\"]\n", ")" ] } ], "metadata": { "kernelspec": { "display_name": "PyTracerLab_app", "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.12.11" } }, "nbformat": 4, "nbformat_minor": 5 }