{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0474d87d",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 0b - Extend the timeseries to 5000 years for every region\n",
    "\n",
    "Uses the scaled dataset to extend it to 5000 years everywhere. \n",
    "\n",
    "*this step has to happen after the notebook `0a_analysis_regional_model_dataset_merging_and_initial_state_comparison.ipynb`*. \n",
    "\n",
    "applied option: use the 101 last years and add them until the end of the timeseries. We use 101 years ecause that is more consistent with the rolling average of 101 years.\n",
    "\n",
    "\n",
    "creates `../data/GMIP3_reg_glacier_model_data/glacierMIP3_Feb12_2024_models_all_rgi_regions_sum_scaled_extended_repeat_last_101yrs.nc`\n",
    "- dataset useful for plotting timeseries or steady-state estimates\n",
    "\n",
    "---\n",
    "other options that we thought of, but did not apply at the end: do a fit to interpolate until 5000 years, but this is not so straightforward and most likely not necessary because the regions that run only until simulation year 2000 mostly reached steady-state until then. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "304d2697",
   "metadata": {},
   "outputs": [],
   "source": [
    "DATE = 'Feb12_2024'\n",
    "# download it here https://cluster.klima.uni-bremen.de/~lschuster/glacierMIP3_analysis/glacierMIP3_Feb12_2024_models_all_rgi_regions_sum_scaled.nc\n",
    "# and change the path to your local path\n",
    "\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "import glob\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from help_functions import pal_models, model_order, d_reg_num_name\n",
    "\n",
    "try:\n",
    "    #path_merged_runs_scaled = f'/home/www/lschuster/glacierMIP3_analysis/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled.nc'\n",
    "    path_merged_runs_scaled = f'../data/GMIP3_reg_glacier_model_data/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled.nc'\n",
    "    ds_reg_models = xr.open_dataset(path_merged_runs_scaled)\n",
    "\n",
    "except:\n",
    "    # add your local path... \n",
    "    path_merged_runs_scaled = f'/home/lilianschuster/Downloads/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled.nc'\n",
    "    ds_reg_models = xr.open_dataset(path_merged_runs_scaled)\n",
    "\n",
    "fill_option = 'repeat_last_101yrs'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85bb5f8f",
   "metadata": {},
   "source": [
    "## Extend by using the last 101 years (because of long-term decadal variability ...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "74e19f34",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "repeat_last_101yrs\n"
     ]
    }
   ],
   "source": [
    "import pymannkendall as mk\n",
    "\n",
    "# these regions should run until 5000 years:\n",
    "rgi_regs_5000 = ['01', '03', '04', '05', '06','07', '09', '17','19']        \n",
    "\n",
    "\n",
    "n_increasing = 0\n",
    "n_decreasing = 0\n",
    "n_no_trend = 0\n",
    "deltaV_l = []\n",
    "deltaV_l50 = []\n",
    "ds_reg_models_extend = ds_reg_models.copy()\n",
    "for gcm in np.arange(0,len(ds_reg_models_extend.gcm.values),1):\n",
    "    for period_scenario in np.arange(0, len(ds_reg_models_extend.period_scenario),1):\n",
    "        for m in np.arange(0, len(ds_reg_models_extend.model_author),1):\n",
    "\n",
    "            ds = ds_reg_models_extend.isel(model_author=m).isel(gcm=gcm).isel(period_scenario=period_scenario)\n",
    "            #print(ds)\n",
    "            for rgi_reg_id,rgi_reg in enumerate(ds.rgi_reg.values): \n",
    "                if np.all(np.isnan(ds.sel(rgi_reg =rgi_reg).volume_m3.values)) and np.all(np.isnan(ds.sel(rgi_reg =rgi_reg).area_m2.values)):\n",
    "                    # ok we do not have any regional data for that region, model_author, gcm, period_scenarios ... just keep the values np.NaN...\n",
    "                    pass\n",
    "                #elif rgi_reg in rgi_regs_5000:\n",
    "                #    assert not np.any(np.isnan(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).volume_m3.values))\n",
    "                else:\n",
    "                    # check that all are not nan-values! -> then do not need to extend\n",
    "                    if not np.any(np.isnan(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).volume_m3.values)):\n",
    "                        pass\n",
    "                    else:\n",
    "                        #try:\n",
    "                        # check that it is really always np.NaN values after simulation year 2000 for that region\n",
    "                        assert np.all(np.isnan(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).volume_m3.values))\n",
    "                        # Huss has area 0 for some regions \n",
    "                        assert np.all(np.isnan(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).area_m2.values)) or np.all(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).area_m2.values==0)\n",
    "\n",
    "                        # fill them up with the last simulation year values \n",
    "                        # we fill that up that later \n",
    "                        #for y in np.arange(2001,5001):\n",
    "                        if fill_option == 'last_value':\n",
    "                            ds['volume_m3'].data[..., rgi_reg_id, 2001:] = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=2000).volume_m3.values\n",
    "                            ds['area_m2'].data[..., rgi_reg_id, 2001:] = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=2000).area_m2.values\n",
    "                        elif fill_option == 'repeat_last_21yrs':  # 21-yr period is repeated ~142.9 times to fill up the additional 3000yrs\n",
    "                            ds['volume_m3'].data[..., rgi_reg_id, 2001:] = np.tile(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1980,2000)).volume_m3.values, 143)[:3000]\n",
    "                            ds['area_m2'].data[..., rgi_reg_id, 2001:] = np.tile(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1980,2000)).area_m2.values, 143)[:3000]\n",
    "                        elif fill_option == 'repeat_last_101yrs':  # 101-yr period is repeated ~29.7 times to fill up the additional 3000yrs\n",
    "                            ds['volume_m3'].data[..., rgi_reg_id, 2001:] = np.tile(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1900,2000)).volume_m3.values, 30)[:3000]\n",
    "                            ds['area_m2'].data[..., rgi_reg_id, 2001:] = np.tile(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1900,2000)).area_m2.values, 30)[:3000]\n",
    "                        #print(gcm, period_scenario, ds_reg_models.isel(model_author=m).model_author.values, rgi_reg)\n",
    "\n",
    "                        #except:\n",
    "                        #    # ok some models did run over all regions for 5000 years\n",
    "                        #    print('runs for 5000 years: ' , ds_reg_models.isel(model_author=m).model_author.values, rgi_reg)\n",
    "                        #    pass\n",
    "\n",
    "                        dend = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1901,2000)).volume_m3\n",
    "                        dendm = dend.mean(dim='simulation_year')\n",
    "                        dendm_e = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1801,1900)).volume_m3.mean(dim='simulation_year')\n",
    "                        deltaV_l.append((dendm.values-dendm_e.values) / (ds.sel(rgi_reg =rgi_reg).volume_m3.sel(simulation_year=0).values))\n",
    "                        mk_output = mk.original_test(dend, alpha=0.01)\n",
    "                        if mk_output.trend =='no trend':\n",
    "                            n_no_trend +=1\n",
    "                        elif mk_output.trend =='decreasing':\n",
    "                            n_decreasing +=1\n",
    "                        elif mk_output.trend =='increasing':\n",
    "                            n_increasing +=1\n",
    "                        #    print(mk_output, gcm, period_scenario, ds_reg_models.isel(model_author=m).model_author.values, rgi_reg)\n",
    "\n",
    "                        dend50 = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1951,2000)).volume_m3\n",
    "                        dendm50 = dend50.mean(dim='simulation_year')\n",
    "                        dendm_e50 = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1901,1950)).volume_m3.mean(dim='simulation_year')\n",
    "                        deltaV_l50.append((dendm50.values-dendm_e50.values) / (ds.sel(rgi_reg =rgi_reg).volume_m3.sel(simulation_year=0).values))\n",
    "\n",
    "            assert np.shape(ds_reg_models_extend['volume_m3'][m,gcm,:,:,period_scenario])== (19,5001)\n",
    "            # add it to the big file \n",
    "            ds_reg_models_extend['volume_m3'].data[m,gcm,:,:,period_scenario] = ds['volume_m3'].values\n",
    "            ds_reg_models_extend['area_m2'].data[m,gcm,:,:,period_scenario] = ds['area_m2'].values\n",
    "\n",
    "ds_reg_models_extend.coords['extend_option'] = fill_option\n",
    "\n",
    "# OGGM_v16 gave global estimates -> after extending the timeseries, there should be no NaN values left \n",
    "assert not np.any(np.isnan(ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').values))\n",
    "\n",
    "if fill_option == 'repeat_last_21yrs': \n",
    "    # check if values are all equal to the last 20 yr timeseries for a RGI region where there are no \n",
    "    for j in np.arange(0,3000,21):\n",
    "        np.testing.assert_allclose(ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='02').sel(simulation_year = slice(1980,2000)).isel(gcm=0).isel(period_scenario=0),\n",
    "                                   ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='02').sel(simulation_year = slice(2001+j,2021+j)).isel(gcm=0).isel(period_scenario=0))\n",
    "\n",
    "    for j in np.arange(0,3000,21):\n",
    "        np.testing.assert_allclose(ds_reg_models_extend.volume_m3.sel(model_author='PyGEM-OGGM_v13').sel(rgi_reg='02').sel(simulation_year = slice(1980,2000)).isel(gcm=0).isel(period_scenario=0),\n",
    "                                   ds_reg_models_extend.volume_m3.sel(model_author='PyGEM-OGGM_v13').sel(rgi_reg='02').sel(simulation_year = slice(2001+j,2021+j)).isel(gcm=0).isel(period_scenario=0))\n",
    "elif fill_option == 'repeat_last_101yrs':\n",
    "    print(fill_option)\n",
    "     # check if values are all equal to the last 101 yr timeseries for a RGI region where there are no \n",
    "    for j in np.arange(0,2902,101):\n",
    "        np.testing.assert_allclose(ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='02').sel(simulation_year = slice(1900,2000)).isel(gcm=0).isel(period_scenario=0),\n",
    "                                   ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='02').sel(simulation_year = slice(2001+j,2101+j)).isel(gcm=0).isel(period_scenario=0)\n",
    "                                  )\n",
    "    \n",
    "    np.testing.assert_allclose(ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='02').sel(simulation_year = slice(1900,2000)).isel(gcm=0).isel(period_scenario=0).mean(dim='simulation_year'),\n",
    "                                   ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='02').sel(simulation_year = slice(4900,5000)).isel(gcm=0).isel(period_scenario=0).mean(dim='simulation_year'),\n",
    "                               rtol=1e-3\n",
    "                                  )\n",
    "    for j in np.arange(0,2902,101):\n",
    "        np.testing.assert_allclose(ds_reg_models_extend.volume_m3.sel(model_author='PyGEM-OGGM_v13').sel(rgi_reg='02').sel(simulation_year = slice(1900,2000)).isel(gcm=0).isel(period_scenario=0),\n",
    "                                   ds_reg_models_extend.volume_m3.sel(model_author='PyGEM-OGGM_v13').sel(rgi_reg='02').sel(simulation_year = slice(2001+j,2101+j)).isel(gcm=0).isel(period_scenario=0)\n",
    "                                  )\n",
    "    np.testing.assert_allclose(ds_reg_models_extend.volume_m3.sel(model_author='PyGEM-OGGM_v13').sel(rgi_reg='02').sel(simulation_year = slice(1900,2000)).isel(gcm=0).isel(period_scenario=0).mean(dim='simulation_year'),\n",
    "                                   ds_reg_models_extend.volume_m3.sel(model_author='PyGEM-OGGM_v13').sel(rgi_reg='02').sel(simulation_year = slice(4900,5000)).isel(gcm=0).isel(period_scenario=0).mean(dim='simulation_year'),\n",
    "                               rtol=1e-3\n",
    "                                  )\n",
    "    \n",
    "ds_reg_models_extend.attrs['description'] = 'Scaled and extended regionally aggregated glacier model projections from GlacierMIP3 (volume and area).'\n",
    "ds_reg_models_extend.attrs['postprocessing_phase'] = ('Volume scaled to match regional Farinotti et al. (2019) multi-model glacier volume estimate at the beginning, '\n",
    "                                                'area scaled to match the RGI6.0 area at the beginning. Scaling done for each climate scenario and model time '\n",
    "                                               'series (experiment) individually. For most experiments, scaling was not necessary as models aimed to match regional '\n",
    "                                               'volume/area, and upscaled already internally. In addition, where necessary, the experiments are extended from the year '\n",
    "                                               '2000 to the year 5000. More in Methods and Supplementary Data Table 2 of Zekollari and Schuster et al., submitted') \n",
    "ds_reg_models_extend.attrs['dataset_version'] = 'v1.0'\n",
    "ds_reg_models_extend.attrs['source'] = 'Data from Glacier Model Intercomparison Project Phase 3 (GlacierMIP3), https://doi.org/10.5281/zenodo.14045269'\n",
    "ds_reg_models_extend.attrs['contact'] = 'lilian.schuster@uibk.ac.at, harry.zekollari@vub.be' \n",
    "\n",
    "try:\n",
    "    ds_reg_models_extend.to_netcdf(f'/home/www/lschuster/glacierMIP3_analysis/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled_extended_{fill_option}.nc')\n",
    "    ds_reg_models_extend.to_netcdf(f'../data/GMIP3_reg_glacier_model_data/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled_extended_{fill_option}.nc')\n",
    "except:\n",
    "    ds_reg_models_extend.to_netcdf(f'/home/lilianschuster/Downloads/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled_extended_{fill_option}.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "74b23e85-e257-47cd-a96b-1c5699b22564",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2239, 389, 252)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# there are still a few simulations with signigicant trends, but it doesn't make sense to correct here\n",
    "# as anyways climate is not stable over such long periods even without antropogenic changes...\n",
    "n_no_trend, n_decreasing, n_increasing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "83f56e4b-1ccc-4f9d-b612-0aa007d6a9f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01875"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# ratio of experiments and models going for 2000 years which are not stabilising after ~2000 years \n",
    "len(np.array(deltaV_l)[list(np.abs(deltaV_l)>0.01)])/len(deltaV_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c605967-acba-43c4-b428-6c66346768f5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "bc2389c4",
   "metadata": {},
   "source": [
    "### Using a fit instead as extension?\n",
    "- this is not so straightforward, as most fits are not working very well (and we would need to find one that works on every simulation). Therefore, we decided to not use that approach. If you are interested, you can still have a look into the code below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "45c1da68",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_671031/2453980801.py:2: RuntimeWarning: divide by zero encountered in log\n",
      "  plt.plot(np.log(x))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5b34279210>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0xklEQVR4nO3deXxb1Z3//7dkW7Id2/K+xUucnaxkJ4SdlBAohZYftDRtA/RLpxAGaKaFpv21lMf8IOnMd5jptEwKnRJmCjSUKVCgBSaEJdCG7AlxQjayOfGWxLHlVbal8/tDthKThciRdO3r1/Px0MOWdGV9dADdN+ece47DGGMEAAAQAU6rCwAAAPZBsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABETH+s3DAQCqqysVGpqqhwOR6zfHgAA9IIxRo2NjSosLJTTeeZ+iZgHi8rKShUXF8f6bQEAQARUVFSoqKjojM/HPFikpqZKChaWlpYW67cHAAC94PV6VVxcHDqPn0nMg0X38EdaWhrBAgCAfubzpjEweRMAAEQMwQIAAEQMwQIAAEQMwQIAAEQMwQIAAERMWMFiyJAhcjgcp9wWLFgQrfoAAEA/EtblpuvWrZPf7w/dLy8v1xe+8AXdcsstES8MAAD0P2EFi5ycnB73lyxZomHDhunyyy+PaFEAAKB/6vUCWe3t7Xr22We1cOHCsy6W4fP55PP5Qve9Xm9v3xIAAPRxvZ68+corr6i+vl633377WY9bvHixPB5P6MY+IQAA2JfDGGN688I5c+bI5XLptddeO+txp+uxKC4uVkNDA0t6AwDQT3i9Xnk8ns89f/dqKOTAgQN6++239dJLL33usW63W263uzdvAwAA+pleBYtly5YpNzdX119/faTrAQAA58AYo/qWDlU1tKna2xr82XVbcvMExTnPvllYtIQdLAKBgJYtW6b58+crPj7mm6MCAGB7xhjVNbersr5NVQ2tqva2hYJDVUNr1882+ToDp339D+aMUm5aYoyrDgo7Gbz99ts6ePCg7rzzzmjUAwCA7TW2BXsaKutbQ+Ghsj54v6qh9ayh4bOyU1zK9yQqPy1JBZ5E5XsSlRBn3cLaYQeLa665Rr2c7wkAgO35Ov2qbmjT4fpWVXWFhsNdP6vq21TZ0KrGts5z+ls5qW4VehJV4ElSvicxFBwKPMEQkZvmljs+LsqfKDyMZQAAcI78AaPaxraTehl69jhUNbTqaFP7Of0tT1KCCjyJKkxPUmF6MCx0/xycnqS8tES54vvfll4ECwAAurR3BoI9DMdbdai+VYeOB38/XN+iQ8eDcxs6A5/fa5+Y4AwGhq6ehYL0JA3+THgY5LbnKdienwoAgNNobffrcH2rDh1v0eH6rgBxvDX0e01jmz5vtD/e6VBeWqIGpyep4KSwUOgJ3i/0JCk9OeGsq1LbGcECAGAb3raOYA9DV1j4bIA41vz5wxTueKeKMpI0OCNZg9OTVJQRvA1OT9LgjCTlpiZadilnf0CwAAD0G20dflXUtajieIsq6lp1sK6l636rDh9vkfccJkWmuONDQSEYIJI0OD059HvWINeA7W2IBIIFAKDP8AeMqr1tqqhr0cG6Fh3q+llxPBgijjT6PvdvZCQnaHBGkorSk7tCw4kAUZSerLSkeIJDFBEsAAAxY4xRQ2vHid6G4y0neh3qgsMWHf6zT3JIdcerODNZxZlJKslMDv6ecSJE2HVSZH9B6wMAIqrTH1BVQ5v2H2vW/mMtOnisuUeQ+Lw1HOKdDhVlJHWFh2BoKDkpSHiSBu7EyP6AYAEACFt7Z0CHjrfowLEW7T/WHPp58FgwPHxer0NOqlvFGSf1OHQHiKxk5acxObI/I1gAAE6rtd2vg3UtOnBScOj+WVnfqrMt5+CKc6okK1lDsoKhobQrPJRkJqsoI1lJrr61WiQih2ABAANYS3un9h1t1v6jLTpQ16wDR08EiGpv21lfm5QQp9KsZA3JGqTS7K6fmckqzR5Er8MARrAAAJvr9Ad06Hir9h5t0t4jzdp3NHjbe6T5c8NDamJ8MDB0B4isZA3JDv7MSXEz1wGnIFgAgA0YY3SkyXdKcNh7tEkHj7WcdRnqjOQElWUP0pCsQV3DFyeCxEBeQRK9Q7AAgH6kydep/Ueb9emRph4BYt/RZjX5zny1RWKCU0OyBmloziANzU5RWfYgleUM0tDsQUpPdsXwE8DuCBYA0Md09z7sqW3Sp7VN2lPbpD1Hgj9rvGdeIMrpkIoyklWW3R0gBqksO0VlOYNUkJYoJ3MeEAMECwCwiD9gdPh4q/YcaQyGh5NuZ1uaOjvFFexxyB6koTnB3oeh2cFhDHc8V1vAWgQLAIgyX6df+4429wgOnx5p1t4jTfJ1Bk77GqdDKs5M1vCcFA3PTdGw3K6fOSnyJCXE+BMA545gAQAR4uv0a++RZu2qadSumkbtrG7SntpGHaxrOeOaD654p4ZmDwoGh64QMTw32AuRmEDvA/ofggUAhKnTH9CBuhbtqm7UzppG7a5p0s6aRu072iz/GRJEamJ8MDScFB6G56aoKCOZ9R5gKwQLADiDQMDocH1rsPehprErSDTp0yNNaj/DEEZqYrxG5aVqZH6qRuamaGReqobnpignlTUfMDAQLABAUl1zu3ZUebW9ytsVJJq0u6ZRLe3+0x6flBCnkXnB4DCyK0iMyktVXhoBAgMbwQLAgNLpD2jf0WZtr/JqR3WjPqny6pMq7xkv43TFOTU0Z5BG5aeGQsSovFQVZSRx+SZwGgQLALbV0NKh7V3BYUe1V59UBSdVnulKjNKsZI3OT9Wo/DSN7goSQ7KSFR/njHHlQP9FsADQ7wUCRvuONYd6Hz6patSOKq8qG06/D0ayK06j81N1QUGaRhekaUxBMEykuPlKBM4X/xUB6Fc6/QF9eqRZ5YcbVF7ZoPLDDdpe6VXzGeZCFGUk6YKCNF3QFSBG56epJDOZYQwgSggWAPosX6dfu2uaTgoRwR6J0w1lJCY4NSo/GB66g8So/FSlJbKYFBBLBAsAfUJbh1+fVHlVXulV+aFgkNhV06gO/6nrQqS44zWmME3jCj0aNzhN4wZ7NDR7EHMhgD6AYAEg5jr8Ae2sbtSWQ/X6uKJBWw7Va3dt02kXl/IkJQTDQ6FHYwd7NH6wR6UMZQB9FsECQFQZY7T/WIu2VNRry6F6bamo17bK0w9nZKe4NG6wJ9QTMbbQo6KMJNaFAPoRggWAiKr1tmnLoYZQkPj4UIMaWjtOOS41MV4Ti9I1sdijCUXpmliUzuJSgA0QLAD0WluHXx8fatDGg8e1+WAwSFSd5hJPV7xTYwvTNLEoXRcWp2tCkUdDsgYxnAHYEMECwDkxxqiyoU0bDxzXhgPHtengcW2r9KrzM/MinA5pRG5qqCfiwuJ0jcxLlSueiZXAQECwAHBavk6/tlV6tfHAcW08GAwTp1v2OifVrSklGZpUEgwR4wZ7NIiFpoABi//6AUiSahvbtGH/iRBRftirdn/PCZZxTofGFKRpSmkwSEwuyWByJYAeCBbAAGSMUUVdq9bur9Pafce0bv9x7TvafMpxmYNcmlySocmlwRAxocijZBdfGwDOjG8IYAAIBIx21zZp7b5jWrv/uNbtq1O1t+ckS4dDGpWXqimlGZpckqEppRkqzUqmNwJAWAgWgA11+AMqP9ygdfvrtHbfca0/UKf6lp6XfMY7HZpQ5NG0skxNH5KpqaWZ8iSz/DWA8xN2sDh8+LAeeughvfHGG2ppadHw4cO1bNkyTZ06NRr1ATgH/oDRtsoGrf70mP726TGt21+nls9sypWUEKfJpemaPiRL08oyNKk4Q0muOIsqBmBXYQWL48ePa9asWbryyiv1xhtvKCcnR7t371ZGRka06gNwGoGA0a7axlCQ+GjvMTW2dfY4xpOUoGlDMjW9LEPTy7I0tjBNCeylASDKwgoWP//5z1VcXKxly5aFHisrK4t4UQB66l4W+2+fHg0GiU+P6Vhze49jUt3xmjE0UzOHZeviYVkalZfKAlQAYi6sYPHqq69qzpw5uuWWW/T+++9r8ODBuueee3TXXXdFqz5gwDrS6NOHe47og13BMPHZyZZJCXGaOiRDF3cFibGFaezuCcByYQWLvXv3aunSpVq4cKF+9KMfad26dbrvvvvkcrk0f/78077G5/PJ5zuxqI7X6z2/igGb8nX6tWH/cb2/Oxgmtlf1/G/FFefU5NJ0zRyarYuHZ2liUTqrWQLocxzGmFP3KT4Dl8ulqVOn6m9/+1vosfvuu0/r1q3T6tWrT/uan/3sZ3rkkUdOebyhoUFpaWm9KBmwB2OMPj3SpFW7jmrV7iNas7dOrR09J1yOG5ymS0fk6JLh2ZpSmqHEBCZbArCG1+uVx+P53PN3WD0WBQUFGjNmTI/HLrjgAv3xj38842sWLVqkhQsX9iisuLg4nLcFbKOhtUMf7j6qVbuO6IPdR1T5mQ27clLdumxEji4bma1Zw7OVneK2qFIA6J2wgsWsWbO0c+fOHo/t2rVLpaWlZ3yN2+2W282XIwam7l6JlZ/U6p0dtVp/4Lj8J23a5Yp3akZZpi4dka3LRuZoVF4qC1IB6NfCChbf+973dPHFF+uxxx7TrbfeqrVr1+qpp57SU089Fa36gH6nrcOvNfvq9O6OWq3cUaOKutYezw/PTdHlI3N02cgcTR+SyVoSAGwlrDkWkvT6669r0aJF2r17t8rKyrRw4cKwrgo51zEaoD+p8bZ1BYla/XXP0R6LU7ninLpoWJauHp2rq0bnqjgz2cJKAaB3zvX8HXawOF8EC9iBMUZ7apv01rZq/e/2Gn18qKHH83lpbl01OldXjsrVrOHZbCMOoN+LyuRNYCALBIw2H6rXW9uqtWJbjfaetBuowyFNLErXVV29EmML05grAWBAIlgAZ9HeGdCafceCYWJ7jWq8J9ZkccU5NWt4luaMzdfVF+QpJ5VJygBAsAA+w9fp1we7jurPW6v09ic1PfbgSHHH64pROZozNl9XjMpRaiK7gQLAyQgWgII9E3/dc1SvfVypFdtq1Og7ESayU1z6wpg8XTM2XxcPy5I7nqs4AOBMCBYYsDr8Af3t02P688eVemtbjRpaO0LP5acl6rrxBZo7Pl+TSzIUx2ZeAHBOCBYYUPwBo4/2HtPrH1fqzfJqHW85ESZyUt26bly+vjixUFNKMtgZFAB6gWCBAeGTKq9e2XRYr2w+3GMCZtYgl+aOz9cXJxRq2pBMeiYA4DwRLGBb1Q1tenXLYb208bB2VDeGHvckJei6rjAxoyyTrcYBIIIIFrCVJl+n3iqv1subDuuvnx5V9/Jvrjinrhqdqy9PHqwrRuUwARMAooRggX7PGKM1++r0h3UVeqO8usfW49OGZOjLk4p0/fgCeZK5NBQAoo1ggX6r1tum/9l4SC+uP6R9J62COTR7kL48abBumjSYfTkAIMYIFuhXOv0BvbvziF5YV6F3d9aGtiAf5IrTDRMLdeu0Yk0qTmc5bQCwCMEC/UJFXYueX3tQ/7PhkI40nriqY0pphr46rVjXjy9goy8A6AP4JkafFQgYrdp9RL9bfUDv7KwNTcTMGuTSzVOKdOvUIg3PTbW2SABADwQL9Dn1Le16cf0hPbvmgA4cawk9fumIbH19eomuviBPrnguEQWAvohggT7j40P1+t3qA3p1S6V8nQFJUmpivG6ZUqx5F5VoWE6KxRUCAD4PwQKW8geMVmyv0W8/3Kt1+4+HHr+gIE3fmlmqGy8sVLKLf00BoL/gGxuWaPZ16sX1FXr6r/t1sC443JEQ59B14wv0rZmlmlySwZUdANAPESwQU1UNrfqvvx3Q82sOyNsW3Jrck5Sgb1xUom/NHKK8tESLKwQAnA+CBWJiV02jlr73qV7bUqnOrrUnyrIH6c5LynTz5MEMdwCATfBtjqj6+FC9nnh3j97aVhN6bEZZpv7PpUN19ehctiYHAJshWCDiuvfueOLdPfpg91FJksMhXTs2X3dfMUwTitKtLRAAEDUEC0SMMUbv7TyiX727RxsOBK/wiHM6dOOFhbrnimEsZgUAAwDBAufNGKMP9xzVv/zvLm2uqJckueKdunVqkf7usmFsBAYAAwjBAudlzd5j+pcVu7R2X50kKTHBqW9eVKq7Lh2qXK7wAIABh2CBXtl08LgeX7ErNIfCFe/UvBkluvuKYcpNJVAAwEBFsEBYdtc06udv7tDbn9RKkuKdDn11WrHuvWq4CjxJFlcHALAawQLnpLaxTf+6YrdeWHdQAROclPmVSYN139UjmEMBAAghWOCsWto79ZtV+/Tkqk/V0u6XJF0zJk8PzR3NpmAAgFMQLHBa/oDRHzcc0v/9352qbfRJkiYWp+vH112g6WWZFlcHAOirCBY4xYYDx/Xwq+UqP+yVJBVnJunBOaP1xQkFbAwGADgrggVCjjT69PM3d+h/NhySJKW643Xf1SP0rYtL5Y6Ps7g6AEB/QLCAOvwB/W71Af3ril1q9AV3HL1lSpEevHa0clLdFlcHAOhPCBYD3IYDdfrRS+XaWdMoSRo/2KNHbhyrySUZFlcGAOiPCBYDVGNbh/7pzZ16ds0BGSOlJyfowTmj9dVpxYpjx1EAQC8RLAagFdtr9JNXylXtbZMk/T9TivTj6y5QxiCXxZUBAPo7gsUAUtvYpkde3a4/b62SJJVmJeuxL4/XrOHZFlcGALALgsUA8ZetVfrxy1t1vKVDcU6H7rp0qO6/eoSSXFztAQCIHGc4B//sZz+Tw+HocRs9enS0akMENLR06P7lm3TPcxt1vKVDYwrS9Oq9s/TDuaMJFQCAiAu7x2Ls2LF6++23T/yBeDo9+qr3dx3Rg/+zRTVen+KcDt1zxTD9/VUj5IoPK08CAHDOwk4F8fHxys/Pj0YtiJC2Dr/+vz9v17MfHZQkDc0epH+5daImcQkpACDKwg4Wu3fvVmFhoRITEzVz5kwtXrxYJSUlZzze5/PJ5/OF7nu93t5VinOyu6ZR9z6/KbQuxe0XD9FD1zLsAQCIjbD6xGfMmKFnnnlGb775ppYuXap9+/bp0ksvVWNj4xlfs3jxYnk8ntCtuLj4vIvGqYwxenF9hb70q79qZ02jslNc+u87p+tnXxpLqAAAxIzDGGN6++L6+nqVlpbq8ccf17e//e3THnO6Hovi4mI1NDQoLS2tt2+NkzT5OvWTV8r18qbDkqRZw7P0r1+9ULmpiRZXBgCwC6/XK4/H87nn7/OaeZmenq6RI0dqz549ZzzG7XbL7Wa/iWjZXdOov3t2g/YeaZbTIS38wkjdfcVwVs8EAFjivC4PaGpq0qeffqqCgoJI1YMwvFlerZue+Kv2HmlWflqiln9npu69agShAgBgmbB6LL7//e/rhhtuUGlpqSorK/Xwww8rLi5Ot912W7Tqw2n4A0b/9vYu/fKdYE/RRUMz9cTXJysrhZ4hAIC1wgoWhw4d0m233aZjx44pJydHl1xyiT766CPl5OREqz58RkNrhx5Yvknv7jwiSbpzVpl+dN1oxcexNgUAwHphBYvly5dHqw6cg4q6Ft3xzDrtqW2SO96pJTeP15cnFVldFgAAISyb2U9sOnhc/+e/1utYc7vy0xL1n/Onatxgj9VlAQDQA8GiH3hja5UeeGGzfJ0BjSlI09O3T1O+h0tJAQB9D8GiDzPG6D8/2KfH3vhExkhXjsrRL78+WSlu/rEBAPomzlB9lDFGS97YoSdX7ZUkfWtmqX76xTFM0gQA9GkEiz7IHzD68ctbtXxdhSTpR9eN1l2XDpXDwfoUAIC+jWDRx7R3BvS9Fzbrz1ur5HRIS74yQbdOY38VAED/QLDoQ9o6/Pq7323Q+7uOKCHOoV98bZKuG8+qpgCA/oNg0Ue0dfh113+v1we7jyopIU5PfnOKLhvJwmMAgP6FYNEHtHX49Z3fbdAHu48q2RWnZ+6YrullmVaXBQBA2LjEwGK+Tr/ufnaDVu06oqSEOC27fRqhAgDQbxEsLNTeGdA9z27UuzuPKDHBqadvn6YZQ7OsLgsAgF4jWFgkEDD6/otbtHJHrdzxTv12/jTNHEaoAAD0bwQLCxhj9I9/3q5Xt1Qq3unQk9+colnDs60uCwCA80awsMB/vPeplv11vyTp/94yUVeMyrW2IAAAIoRgEWN/WFehf35rpyTpJ18co5smDba4IgAAIodgEUN/23NUi17eKkm6+4ph+vYlZRZXBABAZBEsYmTf0Wbd/dxG+QNGN15YqAfnjLK6JAAAIo5gEQMNLR369jPr1NDaoQuL0/XzmyewoRgAwJYIFlHW6Q9owfMbtfdoswo9iXrqW1OUmBBndVkAAEQFwSLKlryxQx/uCS7V/Z/zpyk3NdHqkgAAiBqCRRS9WV6l//xwnyTp8VsnakxhmsUVAQAQXQSLKNl3tFk/ePFjSdLfXTZU145j+3MAgP0RLKKgrSO4sVijr1PTh2Tq+1wBAgAYIAgWUfDIa9u0o7pR2Sku/fLrk5QQRzMDAAYGzngRtmJ7jX6/tkIOh/SLr01SXhqTNQEAAwfBIoKONPr0wz8G51XcdelQNhYDAAw4BIsIMcbooT9+rGPN7Rqdn6p/uGak1SUBABBzBIsIeX7tQb2zo1aueKd+8bVJcsezCBYAYOAhWETA4fpWPfbnTyRJD84ZpVH5qRZXBACANQgW58kYo//35a1qbvdrammG7pzFjqUAgIGLYHGeXvu4Su/uPCJXnFNLbh4vp5PNxQAAAxfB4jwcb27XI69ukyTde9VwDc9lCAQAMLARLM7Dz9/coWPN7RqZl6LvXj7M6nIAALAcwaKXth5q0AvrKyRJj315vFzxNCUAAJwNe8EYo4dfLZcx0k0XFmrqkEyrSwIAoE8gWPTCK5sPa+PBeiW74vTDuRdYXQ4AAH0GwSJMzb5OLXljh6TghM18D3uBAADQjWARpt9+uE81Xp9KMpP17UtYswIAgJOdV7BYsmSJHA6HHnjggQiV07fVNbfrqVV7JUk/mDOKZbsBAPiMXgeLdevW6cknn9SECRMiWU+ftvS9PWrydWpMQZquH19gdTkAAPQ5vQoWTU1Nmjdvnn7zm98oIyMj0jX1SZX1rfqv1QckSQ9eO4oVNgEAOI1eBYsFCxbo+uuv1+zZsz/3WJ/PJ6/X2+PWH/37yt1q7wxoelmmLh+ZY3U5AAD0SfHhvmD58uXauHGj1q1bd07HL168WI888kjYhfUlFXUtenHDIUnSQ9eOksNBbwUAAKcTVo9FRUWF7r//fj333HNKTDy3yywXLVqkhoaG0K2ioqJXhVrpqVV75Q8YXToiW1NKWQwLAIAzCavHYsOGDaqtrdXkyZNDj/n9fq1atUq/+tWv5PP5FBfX80oJt9stt9sdmWotUOttCy3dfc8Vwy2uBgCAvi2sYHH11Vdr69atPR674447NHr0aD300EOnhAo7+O2H+9TeGdCU0gxdNJTeCgAAziasYJGamqpx48b1eGzQoEHKyso65XE7aGjp0LMfBa8EWXDlMOZWAADwOVh58yyWrzuo5na/Ruen6spRuVaXAwBAnxf2VSGf9d5770WgjL6n0x/Qf3etW3HnJWX0VgAAcA7osTiDtz+p0eH6VmUOculLEwutLgcAgH6BYHEGT/91vyTp69NLlJhgv0mpAABEA8HiNLZVNmjtvjrFOx36xkWlVpcDAEC/QbA4jefWHJQkXTsuX/mec1sIDAAAECxO0dru12ubKyUFh0EAAMC5I1h8xhvlVWr0dao4M0kXDc2yuhwAAPoVgsVnvLAuuHz3rVOK2RodAIAwESxOsv9os9bsq5PDId08pcjqcgAA6HcIFid5cUOwt+KyETkqTE+yuBoAAPofgkUXY4xe2RSctHnLVHorAADoDYJFl40H63W4vlWDXHGafUGe1eUAANAvESy6vLYl2FvxhTF5rLQJAEAvESwk+QNGf95aJUm6gX1BAADoNYKFpLX76nSk0ae0xHhdOiLH6nIAAOi3CBaS/rw1OAxy7bh8ueJpEgAAemvAn0WNMXp7e62kYLAAAAC9N+CDxbZKr6q9bUpKiNPFw7KtLgcAgH5twAeLFdtrJEmXjczmahAAAM4TwaIrWLB2BQAA529AB4vD9a3aXuWV0yFdNTrX6nIAAOj3BnSweHdHcNLm5JIMZaW4La4GAID+b0AHiw93H5UkXT6StSsAAIiEARss/AGjv30aDBazRnA1CAAAkTBgg8XHh+rlbetUamK8Jgz2WF0OAAC2MGCDxV/3BHsrLh6Wpfi4AdsMAABE1IA9o37QNb/iEvYGAQAgYgZksGhp79TGg8clSZcOZ34FAACRMiCDxcYD9erwGxV6ElWalWx1OQAA2MaADBbrD9RJkqaVZcrhcFhcDQAA9jEwg8X+4DDI1NIMiysBAMBeBlyw6PQHtKlrfsXUIZkWVwMAgL0MuGCxo7pRze1+pbrjNTIv1epyAACwlQEXLNbvD86vmFSaoTgn8ysAAIikgRcsDgSHQaYxvwIAgIgbcMFiQ1ewmDKEYAEAQKQNqGBR621TVUObHA5pYlG61eUAAGA7AypYbD3cIEkanpOiQe54i6sBAMB+wgoWS5cu1YQJE5SWlqa0tDTNnDlTb7zxRrRqi7juYDGe3UwBAIiKsIJFUVGRlixZog0bNmj9+vW66qqrdOONN2rbtm3Rqi+iyg97JUnjCBYAAERFWOMBN9xwQ4/7jz76qJYuXaqPPvpIY8eOjWhh0VDe3WNRRLAAACAaej3RwO/368UXX1Rzc7NmzpwZyZqi4kijT9Xe4MTNMQVpVpcDAIAthR0stm7dqpkzZ6qtrU0pKSl6+eWXNWbMmDMe7/P55PP5Qve9Xm/vKj1P5ZXB3ophTNwEACBqwr4qZNSoUdq8ebPWrFmju+++W/Pnz9f27dvPePzixYvl8XhCt+Li4vMquLd2VjdKki6gtwIAgKgJO1i4XC4NHz5cU6ZM0eLFizVx4kT94he/OOPxixYtUkNDQ+hWUVFxXgX31q6uYDEqL8WS9wcAYCA47zGBQCDQY6jjs9xut9xu9/m+zXnbWRMMFmw8BgBA9IQVLBYtWqS5c+eqpKREjY2Nev755/Xee+/prbfeilZ9EeEPGO2ubZIkjconWAAAEC1hBYva2lp961vfUlVVlTwejyZMmKC33npLX/jCF6JVX0QcONas9s6AEhOcKs5ItrocAABsK6xg8dvf/jZadUTVrq5hkBG5qXKyVToAAFEzIPYK2VUTHAZhfgUAANE1QIJF98RNrggBACCaBkSw2He0WVJwcSwAABA9tg8Wxhjt7woWQ7KZuAkAQDTZPlgcafKpud0vp0MqziRYAAAQTbYPFgeOtUiSCtOT5I6Ps7gaAADszfbBont+xZCsQRZXAgCA/dk+WDC/AgCA2LF9sOgeCqHHAgCA6LN9sOgeCinLJlgAABBttg4WxhgdOBYMFqX0WAAAEHW2DhZ1ze1qbvdLkoozkyyuBgAA+7N1sKhqaJMk5aS6udQUAIAYsHWwOFzfKim4hgUAAIg+WweLqu5g4Um0uBIAAAYGWweLyq6hEHosAACIDVsHC4ZCAACILVsHC4ZCAACILVsHi8p6hkIAAIgl2waLDn9AtY0ECwAAYsm2waLG26aAkVxxTmUNclldDgAAA4Jtg0X3MEhBeqKcTofF1QAAMDDYNlhUNQQnbhYwcRMAgJixbbCo9fokSflpBAsAAGLFtsHiSFMwWGSnuC2uBACAgcO+waIxGCxyUgkWAADEim2DxdEmggUAALFm22DR3WPBUAgAALFj+2BBjwUAALFjy2DR6Q+orqVdEj0WAADEki2DRV1zu4yRnA4pk1U3AQCIGVsGi9quYZCsFLfiWHUTAICYsWWwOMoaFgAAWMKWweJ41/yKzEEJFlcCAMDAYs9g0dwhSUpPZn4FAACxZMtgUd/VY5GRTI8FAACxZMtgcbylq8ciiR4LAABiyabBIthjkU6PBQAAMRVWsFi8eLGmTZum1NRU5ebm6qabbtLOnTujVVuv1Xf1WGQwxwIAgJgKK1i8//77WrBggT766COtWLFCHR0duuaaa9Tc3Byt+nqlu8cig6tCAACIqfhwDn7zzTd73H/mmWeUm5urDRs26LLLLotoYeeju8eCq0IAAIitsILFZzU0NEiSMjMzz3iMz+eTz+cL3fd6vefzlufkxFUhBAsAAGKp15M3A4GAHnjgAc2aNUvjxo0743GLFy+Wx+MJ3YqLi3v7luekvTOg5na/JC43BQAg1nodLBYsWKDy8nItX778rMctWrRIDQ0NoVtFRUVv3/KcdPdWOB1SWiLBAgCAWOrVUMi9996r119/XatWrVJRUdFZj3W73XK7Y7dnR/caFp6kBDnZgAwAgJgKK1gYY/T3f//3evnll/Xee++prKwsWnX1WmNbMFikJdFbAQBArIUVLBYsWKDnn39ef/rTn5Samqrq6mpJksfjUVJSUlQKDFejr1OSlOI+r3mpAACgF8KaY7F06VI1NDToiiuuUEFBQej2wgsvRKu+sDW1ESwAALBK2EMhfV1TV49FKhM3AQCIOdvtFdI9xyI1kR4LAABizXbBgqEQAACsY7tgEZq8SY8FAAAxZ7tgQY8FAADWsV+wCE3eJFgAABBrtg0W9FgAABB7tgsWjQyFAABgGdsFC9axAADAOvYLFm3MsQAAwCr2CxbMsQAAwDK2Chb+gDkRLOixAAAg5mwVLJrbO0O/02MBAEDs2SpYtPj8kqQ4p0PueFt9NAAA+gVbnX1bO4LBIikhTg6Hw+JqAAAYeGwVLNq6gkVigq0+FgAA/YatzsAngkWcxZUAADAw2SpYtBIsAACwlK2Cha8jIImhEAAArGKrM3DbSZM3AQBA7NkqWDAUAgCAtWwVLNq6hkLc8QQLAACsYLNg0TUU4iJYAABgBVsFi9BQCKtuAgBgCVudgX3MsQAAwFK2ChZtncE5FgyFAABgDVsFi9Z2hkIAALCSrc7A3ZM33QyFAABgCXsFi+6hEIIFAACWsFWwCA2FECwAALCErYKFr7N7HQtbfSwAAPoNW52BQ9ums/ImAACWsFWwYK8QAACsZatg0RbaNp1gAQCAFWwWLLp7LGz1sQAA6DdsdQZuYygEAABL2SxYsI4FAABWslmwoMcCAAAr2SZY+ANGnQEjSXKxVwgAAJYI+wy8atUq3XDDDSosLJTD4dArr7wShbLC1+EPhH4nWAAAYI2wz8DNzc2aOHGinnjiiWjU02snB4uEOIeFlQAAMHDFh/uCuXPnau7cudGo5bx0+E3o9wQnPRYAAFgh7GARLp/PJ5/PF7rv9Xqj8j7dPRZxToecTnosAACwQtT/137x4sXyeDyhW3FxcVTep71ry3SGQQAAsE7Ug8WiRYvU0NAQulVUVETlfbqvCEmIYxgEAACrRH0oxO12y+12R/ttQkMhLoIFAACWsc1Z+MRQiG0+EgAA/U7YPRZNTU3as2dP6P6+ffu0efNmZWZmqqSkJKLFhaO7xyKeORYAAFgm7GCxfv16XXnllaH7CxculCTNnz9fzzzzTMQKC1f35aYMhQAAYJ2wg8UVV1whY8znHxhj3T0WDIUAAGAd25yFQ8EinqEQAACsYqNgweWmAABYzTZnYYZCAACwnm3OwieCBUMhAABYxTbBgnUsAACwnm3OwsyxAADAerY5C3cGWNIbAACr2eYs3D0UwsqbAABYxzbBwt+1u2m80zYfCQCAfsc2Z+HubdMZCQEAwDq2OQ0HQsHCNh8JAIB+xzZnYXosAACwnm1Ow4GujdHiHEzeBADAKrYJFp0MhQAAYDnbnIUDDIUAAGA525yG6bEAAMB6tjkL++mxAADAcrY5DTN5EwAA69kmWDAUAgCA9WxzFmbyJgAA1rPNaZgeCwAArGebszA9FgAAWM82p2F6LAAAsJ5tzsL+0FUhFhcCAMAAZp9g4e/usSBZAABgFfsEC8NQCAAAVrPNWZiVNwEAsJ5tTsN+Jm8CAGA525yFQ0t62+YTAQDQ/9jmNNzpp8cCAACr2eYs7GcTMgAALGefYMHkTQAALGeb0zCTNwEAsJ5tzsL0WAAAYD3bnIbpsQAAwHq2OQuHggWTNwEAsIx9goVhrxAAAKzWq2DxxBNPaMiQIUpMTNSMGTO0du3aSNcVthNDIQQLAACsEnaweOGFF7Rw4UI9/PDD2rhxoyZOnKg5c+aotrY2GvWdMyZvAgBgvbBPw48//rjuuusu3XHHHRozZox+/etfKzk5WU8//XQ06jtnTN4EAMB6YZ2F29vbtWHDBs2ePfvEH3A6NXv2bK1evfq0r/H5fPJ6vT1u0cDkTQAArBdWsDh69Kj8fr/y8vJ6PJ6Xl6fq6urTvmbx4sXyeDyhW3Fxce+rPQsmbwIAYL2ojxssWrRIDQ0NoVtFRUVU3ueOWUN0zxXDlJ3qisrfBwAAny8+nIOzs7MVFxenmpqaHo/X1NQoPz//tK9xu91yu929r/Ac3XPF8Ki/BwAAOLuweixcLpemTJmilStXhh4LBAJauXKlZs6cGfHiAABA/xJWj4UkLVy4UPPnz9fUqVM1ffp0/du//Zuam5t1xx13RKM+AADQj4QdLL761a/qyJEj+ulPf6rq6mpdeOGFevPNN0+Z0AkAAAYehzFdl1PEiNfrlcfjUUNDg9LS0mL51gAAoJfO9fzNalIAACBiCBYAACBiCBYAACBiCBYAACBiCBYAACBiCBYAACBiCBYAACBiCBYAACBiCBYAACBiwl7S+3x1L/Tp9Xpj/dYAAKCXus/bn7dgd8yDRWNjoySpuLg41m8NAADOU2Njozwezxmfj/leIYFAQJWVlUpNTZXD4YjY3/V6vSouLlZFRQV7kEQR7RwbtHPs0NaxQTvHTrTa2hijxsZGFRYWyuk880yKmPdYOJ1OFRUVRe3vp6Wl8S9tDNDOsUE7xw5tHRu0c+xEo63P1lPRjcmbAAAgYggWAAAgYmwTLNxutx5++GG53W6rS7E12jk2aOfYoa1jg3aOHavbOuaTNwEAgH3ZpscCAABYj2ABAAAihmABAAAihmABAAAixhbB4oknntCQIUOUmJioGTNmaO3atVaX1K8sXrxY06ZNU2pqqnJzc3XTTTdp586dPY5pa2vTggULlJWVpZSUFN18882qqanpcczBgwd1/fXXKzk5Wbm5ufrBD36gzs7OWH6UfmXJkiVyOBx64IEHQo/RzpFz+PBhfeMb31BWVpaSkpI0fvx4rV+/PvS8MUY//elPVVBQoKSkJM2ePVu7d+/u8Tfq6uo0b948paWlKT09Xd/+9rfV1NQU64/SZ/n9fv3kJz9RWVmZkpKSNGzYMP3jP/5jj70kaOfeWbVqlW644QYVFhbK4XDolVde6fF8pNr1448/1qWXXqrExEQVFxfrn/7pn86/eNPPLV++3LhcLvP000+bbdu2mbvuusukp6ebmpoaq0vrN+bMmWOWLVtmysvLzebNm811111nSkpKTFNTU+iY7373u6a4uNisXLnSrF+/3lx00UXm4osvDj3f2dlpxo0bZ2bPnm02bdpk/vKXv5js7GyzaNEiKz5Sn7d27VozZMgQM2HCBHP//feHHqedI6Ours6Ulpaa22+/3axZs8bs3bvXvPXWW2bPnj2hY5YsWWI8Ho955ZVXzJYtW8yXvvQlU1ZWZlpbW0PHXHvttWbixInmo48+Mh988IEZPny4ue2226z4SH3So48+arKysszrr79u9u3bZ1588UWTkpJifvGLX4SOoZ175y9/+Yv58Y9/bF566SUjybz88ss9no9EuzY0NJi8vDwzb948U15ebn7/+9+bpKQk8+STT55X7f0+WEyfPt0sWLAgdN/v95vCwkKzePFiC6vq32pra40k8/777xtjjKmvrzcJCQnmxRdfDB3zySefGElm9erVxpjgfwROp9NUV1eHjlm6dKlJS0szPp8vth+gj2tsbDQjRowwK1asMJdffnkoWNDOkfPQQw+ZSy655IzPBwIBk5+fb/75n/859Fh9fb1xu93m97//vTHGmO3btxtJZt26daFj3njjDeNwOMzhw4ejV3w/cv3115s777yzx2Nf+cpXzLx584wxtHOkfDZYRKpd/+M//sNkZGT0+O546KGHzKhRo86r3n49FNLe3q4NGzZo9uzZocecTqdmz56t1atXW1hZ/9bQ0CBJyszMlCRt2LBBHR0dPdp59OjRKikpCbXz6tWrNX78eOXl5YWOmTNnjrxer7Zt2xbD6vu+BQsW6Prrr+/RnhLtHEmvvvqqpk6dqltuuUW5ubmaNGmSfvOb34Se37dvn6qrq3u0tcfj0YwZM3q0dXp6uqZOnRo6Zvbs2XI6nVqzZk3sPkwfdvHFF2vlypXatWuXJGnLli368MMPNXfuXEm0c7REql1Xr16tyy67TC6XK3TMnDlztHPnTh0/frzX9cV8E7JIOnr0qPx+f48vWUnKy8vTjh07LKqqfwsEAnrggQc0a9YsjRs3TpJUXV0tl8ul9PT0Hsfm5eWpuro6dMzp/jl0P4eg5cuXa+PGjVq3bt0pz9HOkbN3714tXbpUCxcu1I9+9COtW7dO9913n1wul+bPnx9qq9O15cltnZub2+P5+Ph4ZWZm0tZdfvjDH8rr9Wr06NGKi4uT3+/Xo48+qnnz5kkS7RwlkWrX6upqlZWVnfI3up/LyMjoVX39Olgg8hYsWKDy8nJ9+OGHVpdiOxUVFbr//vu1YsUKJSYmWl2OrQUCAU2dOlWPPfaYJGnSpEkqLy/Xr3/9a82fP9/i6uzjD3/4g5577jk9//zzGjt2rDZv3qwHHnhAhYWFtPMA1q+HQrKzsxUXF3fKrPmamhrl5+dbVFX/de+99+r111/Xu+++22Nr+/z8fLW3t6u+vr7H8Se3c35+/mn/OXQ/h+BQR21trSZPnqz4+HjFx8fr/fff17//+78rPj5eeXl5tHOEFBQUaMyYMT0eu+CCC3Tw4EFJJ9rqbN8d+fn5qq2t7fF8Z2en6urqaOsuP/jBD/TDH/5QX/va1zR+/Hh985vf1Pe+9z0tXrxYEu0cLZFq12h9n/TrYOFyuTRlyhStXLky9FggENDKlSs1c+ZMCyvrX4wxuvfee/Xyyy/rnXfeOaVrbMqUKUpISOjRzjt37tTBgwdD7Txz5kxt3bq1x7/IK1asUFpa2ilf8APV1Vdfra1bt2rz5s2h29SpUzVv3rzQ77RzZMyaNeuUS6Z37dql0tJSSVJZWZny8/N7tLXX69WaNWt6tHV9fb02bNgQOuadd95RIBDQjBkzYvAp+r6WlhY5nT1PI3FxcQoEApJo52iJVLvOnDlTq1atUkdHR+iYFStWaNSoUb0eBpFkj8tN3W63eeaZZ8z27dvNd77zHZOent5j1jzO7u677zYej8e89957pqqqKnRraWkJHfPd737XlJSUmHfeecesX7/ezJw508ycOTP0fPdlkNdcc43ZvHmzefPNN01OTg6XQX6Ok68KMYZ2jpS1a9ea+Ph48+ijj5rdu3eb5557ziQnJ5tnn302dMySJUtMenq6+dOf/mQ+/vhjc+ONN572cr1JkyaZNWvWmA8//NCMGDFiwF8GebL58+ebwYMHhy43femll0x2drZ58MEHQ8fQzr3T2NhoNm3aZDZt2mQkmccff9xs2rTJHDhwwBgTmXatr683eXl55pvf/KYpLy83y5cvN8nJyVxuaowxv/zlL01JSYlxuVxm+vTp5qOPPrK6pH5F0mlvy5YtCx3T2tpq7rnnHpORkWGSk5PNl7/8ZVNVVdXj7+zfv9/MnTvXJCUlmezsbPMP//APpqOjI8afpn/5bLCgnSPntddeM+PGjTNut9uMHj3aPPXUUz2eDwQC5ic/+YnJy8szbrfbXH311Wbnzp09jjl27Ji57bbbTEpKiklLSzN33HGHaWxsjOXH6NO8Xq+5//77TUlJiUlMTDRDhw41P/7xj3tcvkg7986777572u/l+fPnG2Mi165btmwxl1xyiXG73Wbw4MFmyZIl510726YDAICI6ddzLAAAQN9CsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABFDsAAAABHz/wMWbqFF3ne4DwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=np.arange(0,1000,1)\n",
    "plt.plot(np.log(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bf40ece9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.12705882352941175"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_total = n_increasing+ n_increasing+n_no_trend\n",
    "n_increasing/n_total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c80ee121",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.15899159663865547"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_decreasing/n_total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d41ff02",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "12ef65c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Mann_Kendall_Test(trend='increasing', h=True, p=2.08909639254351e-06, z=4.7446088308081515, Tau=0.3203960396039604, s=1618.0, var_s=116150.0, slope=300303703.2010944, intercept=4719902094583.945)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dend = ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='13').sel(simulation_year = slice(1900,2000)).isel(gcm=0).isel(period_scenario=0)\n",
    "dall = ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='13').sel(simulation_year = slice(0,2000)).isel(gcm=0).isel(period_scenario=0)\n",
    "mk_output = mk.original_test(dend, alpha=0.01)\n",
    "mk_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "64cc6e1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebfa9d8b-f72c-47cc-bbe3-eaef7fe9c4e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymannkendall as mk\n",
    "\n",
    "def func(x, a, b,c,d,e): # x-shifted log\n",
    "    return (a*np.exp(-b*x**d)+c)/e\n",
    "\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "# these regions should run until 5000 years:\n",
    "rgi_regs_5000 = ['01', '03', '04', '05', '06','07', '09', '17','19']        \n",
    "# expand the rgi regions that only run until 2000 by using the last simulated year\n",
    "n_increasing = 0\n",
    "n_decreasing = 0\n",
    "n_no_trend = 0\n",
    "fill_option = 'fit'\n",
    "\n",
    "ds_reg_models = xr.open_dataset(path_merged_runs_scaled)\n",
    "ds_reg_models_extend = ds_reg_models.copy()\n",
    "for gcm in np.arange(0,len(ds_reg_models_extend.gcm.values),1):\n",
    "    for period_scenario in np.arange(0, len(ds_reg_models_extend.period_scenario),1):\n",
    "        for m in np.arange(0, len(ds_reg_models_extend.model_author),1):\n",
    "        \n",
    "            ds = ds_reg_models_extend.isel(model_author=m).isel(gcm=gcm).isel(period_scenario=period_scenario)\n",
    "            #print(ds)\n",
    "            for rgi_reg_id,rgi_reg in enumerate(ds.rgi_reg.values): \n",
    "                if np.all(np.isnan(ds.sel(rgi_reg =rgi_reg).volume_m3.values)) and np.all(np.isnan(ds.sel(rgi_reg =rgi_reg).area_m2.values)):\n",
    "                    # ok we do not have any regional data for that region, model_author, gcm, period_scenarios ... just keep the values np.NaN...\n",
    "                    pass\n",
    "                #elif rgi_reg in rgi_regs_5000:\n",
    "                #    assert not np.any(np.isnan(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).volume_m3.values))\n",
    "                else:\n",
    "                    # check that all are not nan-values! -> then do not need to extend\n",
    "                    if not np.any(np.isnan(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).volume_m3.values)):\n",
    "                        pass\n",
    "                    else:\n",
    "                        #try:\n",
    "                        # check that it is really always np.NaN values after simulation year 2000 for that region\n",
    "                        assert np.all(np.isnan(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).volume_m3.values))\n",
    "                        # Huss has area 0 for some regions \n",
    "                        assert np.all(np.isnan(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).area_m2.values)) or np.all(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(2001,5001)).area_m2.values==0)\n",
    "\n",
    "                        # fill them up with the last simulation year values \n",
    "                        # we fill that up that later \n",
    "                        #for y in np.arange(2001,5001):\n",
    "                        if fill_option == 'last_value':\n",
    "                            ds['volume_m3'].data[..., rgi_reg_id, 2001:] = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=2000).volume_m3.values\n",
    "                            ds['area_m2'].data[..., rgi_reg_id, 2001:] = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=2000).area_m2.values\n",
    "                        elif fill_option == 'repeat_last_20yrs':  # 20-yr period is repeated 150 times to fill up the additional 3000yrs\n",
    "                            ds['volume_m3'].data[..., rgi_reg_id, 2001:] = np.tile(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1981,2000)).volume_m3.values, 150)\n",
    "                            ds['area_m2'].data[..., rgi_reg_id, 2001:] = np.tile(ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1981,2000)).area_m2.values, 150)\n",
    "                        elif fill_option =='fit':\n",
    "                            # let's only take the years 1000 to 2000\n",
    "                            skip_yrs = 1000\n",
    "                            dall = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(skip_yrs,2000))\n",
    "                            #dall = ds_reg_models_extend.sel(model_author='OGGM_v16').sel(rgi_reg='13').sel(simulation_year = slice(0,2000)).isel(gcm=0).isel(period_scenario=0)\n",
    "\n",
    "                            ### volume\n",
    "                            p0 = [-dall.volume_m3.max(), 1, dall.volume_m3.max(),1,1]\n",
    "                            fittedParameters, pcov = curve_fit(func, dall.simulation_year.values, dall.volume_m3.values, p0 = p0)\n",
    "                            t = np.arange(skip_yrs,5001,1)\n",
    "                            d_fit_volume_m3 = func(t, *fittedParameters) \n",
    "                            np.testing.assert_allclose(dall.volume_m3.values[:2000-skip_yrs], \n",
    "                                                       d_fit_volume_m3[0:2000-skip_yrs], rtol=0.1)\n",
    "                            # area\n",
    "                            p0 = [-dall.area_m2.max(), 1, dall.area_m2.max(),1,1]\n",
    "                            fittedParameters, pcov = curve_fit(func, dall.simulation_year.values, dall.area_m2.values, p0 = p0)\n",
    "                            t = np.arange(skip_yrs,5001,1)\n",
    "                            d_fit_area_m2 = func(t, *fittedParameters) \n",
    "                            np.testing.assert_allclose(dall.area_m2.values[:2000-skip_yrs], \n",
    "                                                       d_fit_area_m2[:2000-skip_yrs], rtol=0.1)\n",
    "                            \n",
    "                            ds['volume_m3'].data[..., rgi_reg_id, 2001:] = d_fit_volume_m3[2000-skip_yrs+1:]\n",
    "                            ds['area_m2'].data[..., rgi_reg_id, 2001:] = d_fit_area_m2[2000-skip_yrs+1:]\n",
    "                        #print(gcm, period_scenario, ds_reg_models.isel(model_author=m).model_author.values, rgi_reg)\n",
    "\n",
    "                        #except:\n",
    "                        #    # ok some models did run over all regions for 5000 years\n",
    "                        #    print('runs for 5000 years: ' , ds_reg_models.isel(model_author=m).model_author.values, rgi_reg)\n",
    "                        #    pass\n",
    "                        \n",
    "                        dend = ds.sel(rgi_reg =rgi_reg).sel(simulation_year=slice(1900,2000)).volume_m3\n",
    "                        mk_output = mk.original_test(dend, alpha=0.01)\n",
    "                        if mk_output.trend =='no trend':\n",
    "                            n_no_trend +=1\n",
    "                        elif mk_output.trend =='decreasing':\n",
    "                            n_decreasing +=1\n",
    "                        elif mk_output.trend =='increasing':\n",
    "                            n_increasing +=1\n",
    "                        #    print(mk_output, gcm, period_scenario, ds_reg_models.isel(model_author=m).model_author.values, rgi_reg)\n",
    "\n",
    "            assert np.shape(ds_reg_models_extend['volume_m3'][m,gcm,:,:,period_scenario])== (19,5001)\n",
    "            # add it to the big file \n",
    "            ds_reg_models_extend['volume_m3'].data[m,gcm,:,:,period_scenario] = ds['volume_m3'].values\n",
    "            ds_reg_models_extend['area_m2'].data[m,gcm,:,:,period_scenario] = ds['area_m2'].values\n",
    "            \n",
    "ds_reg_models_extend.coords['extend_option'] = fill_option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62d5d8b2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0a52a140",
   "metadata": {},
   "source": [
    "**below you can see one of the issues**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 313,
   "id": "f029d37a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3102955/2645228258.py:4: RuntimeWarning: overflow encountered in power\n",
      "  return (a*np.exp(-b*x**d)+c)/e\n",
      "/tmp/ipykernel_3102955/2645228258.py:4: RuntimeWarning: overflow encountered in exp\n",
      "  return (a*np.exp(-b*x**d)+c)/e\n",
      "/home/users/lschuster/.local/lib/python3.10/site-packages/scipy/optimize/_minpack_py.py:881: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fce5279b190>]"
      ]
     },
     "execution_count": 313,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGsCAYAAAD+L/ysAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8GElEQVR4nO3de3xU1aH3/+/MJJkkQhIwISEQCIiAIASFkkbBymM0XH4c9fTnQaSCacUHKqdKqpYoF62tse0xxXpQbAuiPj8FVKR9DjSWRkE5BjgGolIugqDhlnAzmSRCEmbW7w9kYEhCMmFuST7vV/erM3uvvWbtVc7J97X22mtbjDFGAAAAIcwa7AYAAAA0h8ACAABCHoEFAACEPAILAAAIeQQWAAAQ8ggsAAAg5BFYAABAyCOwAACAkEdgAQAAIY/AAgAAQl67CywffvihJk6cqOTkZFksFq1evdqr80+fPq377rtPQ4YMUVhYmO64444GZVatWqVbb71VCQkJiomJUUZGht577z3fXAAAAGig3QWWmpoapaWladGiRa063+l0KioqSj/72c+UmZnZaJkPP/xQt956q9auXavi4mKNGTNGEydO1LZt2y6n6QAAoAmW9vzyQ4vFonfffddjlKS2tlZPPPGE3nzzTVVUVOjaa6/Vb37zG918880Nzr/vvvtUUVHRolGawYMHa9KkSZo/f77vLgAAAEhqhyMszZk1a5aKioq0fPlyffbZZ7rrrrs0duxY7dmzp9V1ulwuVVVVqWvXrj5sKQAAOKdDBZbS0lK98soreuuttzR69GhdddVVeuSRRzRq1Ci98sorra73P/7jP1RdXa1/+7d/82FrAQDAOWHBbkAgff7553I6nerfv7/H/traWl155ZWtqvONN97QU089pb/85S/q1q2bL5oJAAAu0qECS3V1tWw2m4qLi2Wz2TyOderUyev6li9frvvvv19vvfVWkxN0AQDA5etQgeW6666T0+nU0aNHNXr06Muq680339SPf/xjLV++XBMmTPBRCwEAQGPaXWCprq7W3r173d/379+vkpISde3aVf3799eUKVM0depUPffcc7ruuut07NgxFRYWaujQoe7gsWPHDtXV1enkyZOqqqpSSUmJJGnYsGGSzt4GmjZtmp5//nmlp6errKxMkhQVFaXY2NiAXi8AAB1Bu3usef369RozZkyD/dOmTdOyZctUX1+vX/3qV3rttdd06NAhxcfH6/vf/76eeuopDRkyRJKUmpqqr7/+ukEd57rq5ptv1oYNG5r8DQAA4FvtLrAAAID2p0M91gwAANomAgsAAAh57WLSrcvl0uHDh9W5c2dZLJZgNwcAALSAMUZVVVVKTk6W1XrpMZR2EVgOHz6slJSUYDcDAAC0woEDB9SzZ89LlmkXgaVz586Szl5wTExMkFsDAABawuFwKCUlxf13/FLaRWA5dxsoJiaGwAIAQBvTkukcTLoFAAAhj8ACAABCHoEFAACEPAILAAAIeQQWAAAQ8ggsAAAg5BFYAABAyCOwAACAkEdgAQAAIY/AAgAAQh6BBQAAhDwCCwAACHkElmaccbq0ZON+/fNwZbCbAgBAh9Uu3tbsT29uKdXT/7VDkvTVsxOC3BoAADomRliasf2QI9hNAACgwyOwNMPIBLsJAAB0eAQWAAAQ8ggszTAMsAAAEHQElmaQVwAACD4CCwAACHkElma8XXww2E0AAKDDI7AAAICQ53Vg+fDDDzVx4kQlJyfLYrFo9erVzZ6zfv16XX/99bLb7erXr5+WLVvmcfzJJ5+UxWLx2AYOHOht0wAAQDvldWCpqalRWlqaFi1a1KLy+/fv14QJEzRmzBiVlJTo4Ycf1v3336/33nvPo9zgwYN15MgR97Zx40ZvmwYAANopr5fmHzdunMaNG9fi8osXL1afPn303HPPSZKuueYabdy4Ub///e+VlZV1viFhYUpKSvK2OQAAoAPw+xyWoqIiZWZmeuzLyspSUVGRx749e/YoOTlZffv21ZQpU1RaWtpknbW1tXI4HB4bAABov/weWMrKypSYmOixLzExUQ6HQ6dOnZIkpaena9myZSooKNBLL72k/fv3a/To0aqqqmq0zry8PMXGxrq3lJQUv7TdsGocAAAhISSeEho3bpzuuusuDR06VFlZWVq7dq0qKiq0cuXKRsvn5uaqsrLSvR04cMAv7XKRVwAACAlez2HxVlJSksrLyz32lZeXKyYmRlFRUY2eExcXp/79+2vv3r2NHrfb7bLb7T5v68WcJBYAAEKC30dYMjIyVFhY6LFv3bp1ysjIaPKc6upqffnll+revbu/mwcAANoArwNLdXW1SkpKVFJSIunsY8slJSXuSbK5ubmaOnWqu/yMGTO0b98+PfbYY9q1a5defPFFrVy5UrNnz3aXeeSRR7RhwwZ99dVX+vjjj3XnnXfKZrNp8uTJl3l5lyciLCTumAEA0OF5fUvok08+0ZgxY9zfc3JyJEnTpk3TsmXLdOTIEY8nfPr06aM1a9Zo9uzZev7559WzZ0/9+c9/9nik+eDBg5o8ebJOnDihhIQEjRo1Sps2bVJCQsLlXBsAAGgnLKYdPArjcDgUGxuryspKxcTE+LTu1Dlr3J+/enaCT+sGAKAj8+bvN/c8AABAyCOwAACAkEdgAQAAIY/AAgAAQh6BBQAAhDwCixdO1tQFuwkAAHRIBBYv1NSeCXYTAADokAgsXqh3uoLdBAAAOiQCixdcbX+NPQAA2iQCixcYYAEAIDgILF5ghAUAgOAgsHjh/356ONhNAACgQyKweOHF9V8GuwkAAHRIBBYAABDyCCwAACDkEVgAAEDII7A0o7M9LNhNAACgwyOwNKPw5z8IdhMAAOjwCCzN6BYTGewmAADQ4RFYAABAyCOwAACAkEdgAQAAIY/A0pzSzfqRbV2wWwEAQIdGYLmUI59KS2/T/LDXFKNqDe/dJdgtAgCgQyKwXEr3NCn8CkVYnLrSUqXMaxKD3SIAADokAktzbOcXjjMyQWwIAAAdF4HFC4a8AgBAUBBYAABAyCOweMEwxAIAQFAQWFrIIsMtIQAAgoTA0iyL+5OLwAIAQFB4HVg+/PBDTZw4UcnJybJYLFq9enWz56xfv17XX3+97Ha7+vXrp2XLljUos2jRIqWmpioyMlLp6enasmWLt03zO54SAgAgOLwOLDU1NUpLS9OiRYtaVH7//v2aMGGCxowZo5KSEj388MO6//779d5777nLrFixQjk5OVqwYIG2bt2qtLQ0ZWVl6ejRo942z6+4JQQAQHCENV/E07hx4zRu3LgWl1+8eLH69Omj5557TpJ0zTXXaOPGjfr973+vrKwsSVJ+fr6mT5+u7Oxs9zlr1qzR0qVLNWfOHG+b6DfkFQAAgsPvc1iKioqUmZnpsS8rK0tFRUWSpLq6OhUXF3uUsVqtyszMdJe5WG1trRwOh8fmbxYZhlgAAAgSvweWsrIyJSZ6LmmfmJgoh8OhU6dO6fjx43I6nY2WKSsra7TOvLw8xcbGureUlBS/tV+W85NuiSsAAARHm3xKKDc3V5WVle7twIEDAfldBlgAAAgOr+eweCspKUnl5eUe+8rLyxUTE6OoqCjZbDbZbLZGyyQlJTVap91ul91u91ubm8JTQgAABIffR1gyMjJUWFjosW/dunXKyMiQJEVERGj48OEeZVwulwoLC91lQgUjLAAABIfXgaW6ulolJSUqKSmRdPax5ZKSEpWWlko6e7tm6tSp7vIzZszQvn379Nhjj2nXrl168cUXtXLlSs2ePdtdJicnR3/605/06quvaufOnZo5c6ZqamrcTw2FChaOAwAgOLy+JfTJJ59ozJgx7u85OTmSpGnTpmnZsmU6cuSIO7xIUp8+fbRmzRrNnj1bzz//vHr27Kk///nP7keaJWnSpEk6duyY5s+fr7KyMg0bNkwFBQUNJuIGG7eEAAAIDotpB2/0czgcio2NVWVlpWJiYnxb+W/6SKdO6pba3ylz9Gjljr/Gt/UDANBBefP3u00+JRQsbT7ZAQDQRhFYvNAOBqMAAGiTCCxeIK8AABAcBJYWsjDlFgCAoCGwNOfCpflJLAAABAWBxQuMsQAAEBwEFi8wwgIAQHAQWLzAU0IAAAQHgaWFLGIdFgAAgoXA0iwm3QIAEGwEFi8w6RYAgOAgsHiBERYAAIKDwOIF8goAAMFBYPECIywAAAQHgaWFLDJijAUAgOAgsDSHpfkBAAg6AosXXCQWAACCgsDihZWfHFTV6fpgNwMAgA6HwOKlF9d/GewmAADQ4RBYWsjy3YTbE9W1QW4JAAAdD4GlWZaLvlmaKAcAAPyFwOIlKz0GAEDA8efXSxYLIywAAAQagcVLxBUAAAKPwOIlKyMsAAAEHIGlhc7FlIpTrMMCAECgEViaY7n4KSEAABBoBBYvWUksAAAEHIHFS8xhAQAg8AgsAAAg5BFYWujc0vy8rxkAgMAjsDTL8xaQyxBZAAAItFYFlkWLFik1NVWRkZFKT0/Xli1bmixbX1+vX/7yl7rqqqsUGRmptLQ0FRQUeJR58sknZbFYPLaBAwe2pml+R14BACDwvA4sK1asUE5OjhYsWKCtW7cqLS1NWVlZOnr0aKPl586dq5dfflkvvPCCduzYoRkzZujOO+/Utm3bPMoNHjxYR44ccW8bN25s3RUBAIB2x+vAkp+fr+nTpys7O1uDBg3S4sWLFR0draVLlzZa/vXXX9fjjz+u8ePHq2/fvpo5c6bGjx+v5557zqNcWFiYkpKS3Ft8fHzrrsjPhvaMDXYTAADocLwKLHV1dSouLlZmZub5CqxWZWZmqqioqNFzamtrFRkZ6bEvKiqqwQjKnj17lJycrL59+2rKlCkqLS1tsh21tbVyOBwem7+dm3R7hT3M778FAAA8eRVYjh8/LqfTqcTERI/9iYmJKisra/ScrKws5efna8+ePXK5XFq3bp1WrVqlI0eOuMukp6dr2bJlKigo0EsvvaT9+/dr9OjRqqqqarTOvLw8xcbGureUlBRvLuOyMOkWAIDA8/tTQs8//7yuvvpqDRw4UBEREZo1a5ays7NltZ7/6XHjxumuu+7S0KFDlZWVpbVr16qiokIrV65stM7c3FxVVla6twMHDvjvAi5aKI68AgBA4HkVWOLj42Wz2VReXu6xv7y8XElJSY2ek5CQoNWrV6umpkZff/21du3apU6dOqlv375N/k5cXJz69++vvXv3NnrcbrcrJibGYwsU8goAAIHnVWCJiIjQ8OHDVVhY6N7ncrlUWFiojIyMS54bGRmpHj166MyZM3rnnXd0++23N1m2urpaX375pbp37+5N8wKDIRYAAALO61tCOTk5+tOf/qRXX31VO3fu1MyZM1VTU6Ps7GxJ0tSpU5Wbm+suv3nzZq1atUr79u3TRx99pLFjx8rlcumxxx5zl3nkkUe0YcMGffXVV/r444915513ymazafLkyT64RN9ykVcAAAg4rx95mTRpko4dO6b58+errKxMw4YNU0FBgXsibmlpqcf8lNOnT2vu3Lnat2+fOnXqpPHjx+v1119XXFycu8zBgwc1efJknThxQgkJCRo1apQ2bdqkhISEy79CHzk3k8UwwgIAQMBZTDv4C+xwOBQbG6vKykrfz2d57hqp6rAm1D6jf5pULZg4SNk39vHtbwAA0AF58/ebdwl5qe3HOwAA2h4CSwsN6t5ZEk8JAQAQDASWFoqMODvd50jFqSC3BACAjofA0kJbvz4pSfrzxv1BbgkAAB0PgQUAAIQ8AktzLlqaHwAABB6BBQAAhDwCCwAACHkEFgAAEPIILC1kYQUWAACChsDSLCbdAgAQbAQWAAAQ8ggsAAAg5BFYWqhHXFSwmwAAQIdFYGmh24clS5KyBicGuSUAAHQ8BJbmfLfSre27/3a6gtkYAAA6JgJLC51bod8YHm8GACDQCCwtdC6wuAgsAAAEHIGlhazfJRYXeQUAgIAjsLSQ1XI2qTDCAgBA4BFYWuj8CAuBBQCAQCOwNOtsULGcCyw8JQQAQMARWFrIKkZYAAAIFgJLC1l5SggAgKAhsLSUO7AEtxkAAHREBJYWYoQFAIDgIbA057ugYmOEBQCAoCGwtND5p4RILAAABBqBpYVYmh8AgOAhsLQQS/MDABA8BJYWcn23YtzOI44gtwQAgI6nVYFl0aJFSk1NVWRkpNLT07Vly5Ymy9bX1+uXv/ylrrrqKkVGRiotLU0FBQWXVWcw/PMwQQUAgGDxOrCsWLFCOTk5WrBggbZu3aq0tDRlZWXp6NGjjZafO3euXn75Zb3wwgvasWOHZsyYoTvvvFPbtm1rdZ2BdfZWUL2Te0EAAASL14ElPz9f06dPV3Z2tgYNGqTFixcrOjpaS5cubbT866+/rscff1zjx49X3759NXPmTI0fP17PPfdcq+sMBibbAgAQPF4Flrq6OhUXFyszM/N8BVarMjMzVVRU1Og5tbW1ioyM9NgXFRWljRs3XladDofDY/O363vF+f03AABA47wKLMePH5fT6VRiYqLH/sTERJWVlTV6TlZWlvLz87Vnzx65XC6tW7dOq1at0pEjR1pdZ15enmJjY91bSkqKN5fRKr2vvMLvvwEAABrn96eEnn/+eV199dUaOHCgIiIiNGvWLGVnZ8tqbf1P5+bmqrKy0r0dOHDAhy1unEVnbwmFnVujHwAABIxXqSE+Pl42m03l5eUe+8vLy5WUlNToOQkJCVq9erVqamr09ddfa9euXerUqZP69u3b6jrtdrtiYmI8Nr/5bv2VczGFmSwAAASeV4ElIiJCw4cPV2FhoXufy+VSYWGhMjIyLnluZGSkevTooTNnzuidd97R7bffftl1BhQr3QIAEDRh3p6Qk5OjadOmacSIERo5cqQWLlyompoaZWdnS5KmTp2qHj16KC8vT5K0efNmHTp0SMOGDdOhQ4f05JNPyuVy6bHHHmtxnaHg3Eq35BUAAALP68AyadIkHTt2TPPnz1dZWZmGDRumgoIC96TZ0tJSj/kpp0+f1ty5c7Vv3z516tRJ48eP1+uvv664uLgW1xkKzr38UJKMMR7fAQCAf1mMaftjBg6HQ7GxsaqsrPT9fJbn06RvvpLjnrUaurRCkvTlM+NlY/ItAACXxZu/37xLqIWsF42wAACAwCGwNOu7p4Qu6Cne2AwAQGARWFqBJ4UAAAgsAksLWcWcFQAAgoXA0kJWy/lRla9O1ASxJQAAdDwElpa6YIClcOfR4LUDAIAOiMDSnIuW5pd4SggAgEAjsLSQ5YLI4nQFsSEAAHRABJYWunBhW54SAgAgsAgsLeSxcFwQ2wEAQEdEYGkhy4UxhREWAAACisDSrO8m3TLCAgBA0BBYWoE5LAAABBaBpRV4lxAAAIFFYGkFRlgAAAgsAktLXRBSyCsAAAQWgaUVXNwTAgAgoAgszbE0fEvzwO4xQWgIAAAdF4HFC4OTzwaVTvawILcEAICOhcDihSsizgYVXn4IAEBgEVi8cO7uEFNYAAAILAJLi5kLAguJBQCAQCKwNOv8pNtzL0AkrgAAEFgEFi+4AwsjLAAABBSBxQvcEgIAIDgILF44N8LicgW5IQAAdDAElpYyRlZGWAAACAoCixfOz2EJckMAAOhgCCzNuWBpfsu5W0IkFgAAAorA4gUrC8cBABAUBBYvWBlhAQAgKFoVWBYtWqTU1FRFRkYqPT1dW7ZsuWT5hQsXasCAAYqKilJKSopmz56t06dPu48/+eSTslgsHtvAgQNb0zS/sn7XW6zDAgBAYHn92uEVK1YoJydHixcvVnp6uhYuXKisrCzt3r1b3bp1a1D+jTfe0Jw5c7R06VLdcMMN+uKLL3TffffJYrEoPz/fXW7w4MH6xz/+cb5hYaH2RmRzwRyWIDcFAIAOxusRlvz8fE2fPl3Z2dkaNGiQFi9erOjoaC1durTR8h9//LFuvPFG3XPPPUpNTdVtt92myZMnNxiVCQsLU1JSknuLj49v3RX5nKXBJ24JAQAQWF4Flrq6OhUXFyszM/N8BVarMjMzVVRU1Og5N9xwg4qLi90BZd++fVq7dq3Gjx/vUW7Pnj1KTk5W3759NWXKFJWWljbZjtraWjkcDo8tEHisGQCA4PDqvsvx48fldDqVmJjosT8xMVG7du1q9Jx77rlHx48f16hRo2SM0ZkzZzRjxgw9/vjj7jLp6elatmyZBgwYoCNHjuipp57S6NGjtX37dnXu3LlBnXl5eXrqqae8abpPfPNtnSRp24GKgP82AAAdmd+fElq/fr2eeeYZvfjii9q6datWrVqlNWvW6Omnn3aXGTdunO666y4NHTpUWVlZWrt2rSoqKrRy5cpG68zNzVVlZaV7O3DggL8vQ5L00Z7jkqT/++nhgPweAAA4y6sRlvj4eNlsNpWXl3vsLy8vV1JSUqPnzJs3T/fee6/uv/9+SdKQIUNUU1OjBx54QE888YSs1oaZKS4uTv3799fevXsbrdNut8tut3vT9MvHfSAAAILGqxGWiIgIDR8+XIWFhe59LpdLhYWFysjIaPScb7/9tkEosdlskpp+PLi6ulpffvmlunfv7k3zAmrv0epgNwEAgA7D62eHc3JyNG3aNI0YMUIjR47UwoULVVNTo+zsbEnS1KlT1aNHD+Xl5UmSJk6cqPz8fF133XVKT0/X3r17NW/ePE2cONEdXB555BFNnDhRvXv31uHDh7VgwQLZbDZNnjzZh5faShcszX+h2jPOADcEAICOy+vAMmnSJB07dkzz589XWVmZhg0bpoKCAvdE3NLSUo8Rlblz58pisWju3Lk6dOiQEhISNHHiRP361792lzl48KAmT56sEydOKCEhQaNGjdKmTZuUkJDgg0v0j8hwW7CbAABAh2Ex7WDZVofDodjYWFVWViomJsa3lS9Kl47tkqb9l376cbTWfl4mSdrw6M3qfeUVvv0tAAA6EG/+fvMuoRYzur5XF/c3axO3igAAgO8RWAAAQMgjsDTr/EjKhTfPXMbodD0TbwEACAQCSytN/uMmDZxXoJM1dcFuCgAA7R6BxQtG54dYDleeliT9Y2d5U8UBAICPEFi80PafpwIAoG0isLQUaQUAgKAhsDTngseXG4ssPNwMAID/EVi8wCALAADBQWDxgml0jAUAAPgbgQUAAIQ8AkuLmUZvCVlYoh8AAL8jsHhh0vdSGuxrB++OBAAg5BFYmnV+BCW+k12d7WFBbAsAAB0TgcVL4WGeXcYtIQAA/I/AAgAAQh6BxUsu5qwAABBwBJaW+i6okFcAAAg8AktzLpqjcvEICzNYAADwPwKLlxhhAQAg8AgsXmIOCwAAgUdg8RJ5BQCAwCOwtNjZpHLxCAsjLgAA+B+BxUu1Z1we38krAAD4H4GlWZd+DshJYgEAwO8ILJfJ6SKwAADgbwSWy8TbmgEA8D8CS0s1EUwYYQEAwP8ILJeJvAIAgP8RWJrTzNr7PNYMAID/EVi8lBwb6fGdwAIAgP+1KrAsWrRIqampioyMVHp6urZs2XLJ8gsXLtSAAQMUFRWllJQUzZ49W6dPn76sOoPFavUccnG6migIAAB8xuvAsmLFCuXk5GjBggXaunWr0tLSlJWVpaNHjzZa/o033tCcOXO0YMEC7dy5U0uWLNGKFSv0+OOPt7rOUGLECAsAAP7mdWDJz8/X9OnTlZ2drUGDBmnx4sWKjo7W0qVLGy3/8ccf68Ybb9Q999yj1NRU3XbbbZo8ebLHCIq3dQZH48GEO0IAAPifV4Glrq5OxcXFyszMPF+B1arMzEwVFRU1es4NN9yg4uJid0DZt2+f1q5dq/Hjx7e6ztraWjkcDo8tWFw8JgQAgN+FeVP4+PHjcjqdSkxM9NifmJioXbt2NXrOPffco+PHj2vUqFEyxujMmTOaMWOG+5ZQa+rMy8vTU0895U3TL4PnnBXLRU8NEVcAAPA/vz8ltH79ej3zzDN68cUXtXXrVq1atUpr1qzR008/3eo6c3NzVVlZ6d4OHDjgwxZ7h6eEAADwP69GWOLj42Wz2VReXu6xv7y8XElJSY2eM2/ePN177726//77JUlDhgxRTU2NHnjgAT3xxBOtqtNut8tut3vTdL/hjhAAAP7n1QhLRESEhg8frsLCQvc+l8ulwsJCZWRkNHrOt99+K6vV82dsNpuks+/haU2dQdFEMOFdQgAA+J9XIyySlJOTo2nTpmnEiBEaOXKkFi5cqJqaGmVnZ0uSpk6dqh49eigvL0+SNHHiROXn5+u6665Tenq69u7dq3nz5mnixInu4NJcnaGkvLLW4zu3hAAA8D+vA8ukSZN07NgxzZ8/X2VlZRo2bJgKCgrck2ZLS0s9RlTmzp0ri8WiuXPn6tChQ0pISNDEiRP161//usV1BtVFs2zrLlopjrwCAID/WUw7uKfhcDgUGxuryspKxcTE+Lbyl2+SjnwqTXlHujpTqXPWeBye8YOrNGfcQN/+JgAAHYA3f795l9Blagd5DwCAkEdguUzMYQEAwP8ILC3G0vwAAAQLgaVZlkseJa8AAOB/BJbLxAgLAAD+R2C5TIYxFgAA/I7AcpkYYQEAwP8ILC1FMgEAIGgILAAAIOQRWJpz0dL8/3p9D4/vLBwHAID/ef0uoY7u13cM0W2DkrRl/0kt/e/9TLkFACAAGGHxUlSETWOvTVIn+9k3TTPAAgCA/xFYWuyiZGK59IJyAADAdwgsrXQurrAOCwAA/kdgaVbjIynnBli4JQQAgP8RWFrJ8l2QIa8AAOB/BJbLxAgLAAD+R2BpJebcAgAQOASWlrpoKOV8XmGIBQAAfyOwtBKTbgEACBwCS3OaufdDYAEAwP8ILK1ksZx7SojEAgCAvxFYLtPKTw7qLyWHgt0MAADaNQJLi1006faCO0UPLS8JbFMAAOhgCCytZGliBVwAAOB7BJZmNR5MHKfrPb4bZt8CAOA3BJZWemn9l8FuAgAAHQaBxUcYYAEAwH8ILD5CXgEAwH8ILC3VzBAKc1gAAPAfAouPEFcAAPAfAktzWvhaZgZYAADwn1YFlkWLFik1NVWRkZFKT0/Xli1bmix78803y2KxNNgmTJjgLnPfffc1OD527NjWNC1gMvpeGewmAADQYXgdWFasWKGcnBwtWLBAW7duVVpamrKysnT06NFGy69atUpHjhxxb9u3b5fNZtNdd93lUW7s2LEe5d58883WXVGATPpeisd33ikEAID/eB1Y8vPzNX36dGVnZ2vQoEFavHixoqOjtXTp0kbLd+3aVUlJSe5t3bp1io6ObhBY7Ha7R7kuXbq07or8xjOQpKXEeR4lrwAA4DdeBZa6ujoVFxcrMzPzfAVWqzIzM1VUVNSiOpYsWaK7775bV1xxhcf+9evXq1u3bhowYIBmzpypEydONFlHbW2tHA6HxxZoLMwPAEDgeBVYjh8/LqfTqcTERI/9iYmJKisra/b8LVu2aPv27br//vs99o8dO1avvfaaCgsL9Zvf/EYbNmzQuHHj5HQ6G60nLy9PsbGx7i0lJaXRcr7ReDS5eC4uIywAAPhPWCB/bMmSJRoyZIhGjhzpsf/uu+92fx4yZIiGDh2qq666SuvXr9ctt9zSoJ7c3Fzl5OS4vzscDj+HluYxhwUAAP/xaoQlPj5eNptN5eXlHvvLy8uVlJR0yXNramq0fPly/eQnP2n2d/r27av4+Hjt3bu30eN2u10xMTEeW6BdPKLy4RfH9U1NXcDbAQBAR+BVYImIiNDw4cNVWFjo3udyuVRYWKiMjIxLnvvWW2+ptrZWP/rRj5r9nYMHD+rEiRPq3r27N83zr2bu+cz4P8W67ul1AWoMAAAdi9dPCeXk5OhPf/qTXn31Ve3cuVMzZ85UTU2NsrOzJUlTp05Vbm5ug/OWLFmiO+64Q1de6bl+SXV1tR599FFt2rRJX331lQoLC3X77berX79+ysrKauVl+R83gAAACByv57BMmjRJx44d0/z581VWVqZhw4apoKDAPRG3tLRUVqtnDtq9e7c2btyov//97w3qs9ls+uyzz/Tqq6+qoqJCycnJuu222/T000/Lbre38rJ8qImVbpt6d5AxRpYWro4LAABaplWTbmfNmqVZs2Y1emz9+vUN9g0YMKDJP/BRUVF67733WtOMkOR0GYXZCCwAAPgS7xJqpaZuCTl5vhkAAJ8jsPhYSWlFsJsAAEC7Q2BpMc+Rk6YGUgp3Nf5OJQAA0HoEFh9zurglBACArxFYmtXUBNrGg4mLOSwAAPgcgaWVesRFN7qfvAIAgO8RWFopKsKmrfNubbCfW0IAAPgegaWlGhk66XpFRIN93BICAMD3CCw+xgALAAC+R2BpjpfL7De1oi8AAGg9AouPMYcFAADfI7D42E39E4LdBAAA2h0Ci49FhNGlAAD4Gn9dW6xlt3qYwwIAgO8RWC7TG/ene3wnrwAA4HsElmZd+imhgd1jPL4z5xYAAN8jsFymi+OMaeGtIwAA0HIElstkvWidljNOo3qnK0itAQCgfQoLdgPajC/fl7492WB3eL1Tk2073N83v12o/3nHoqdvv1ZW79acAwAgdFlt0vVTg/bzBJbmhH33vqDiZY0ejpaUF97IgTX+ahAAAEFgsxNYQtpNj0mRsZLL2ejhMy6XCncebbD/5gEJsofZ/N06AAACwxrcyEBgaU7qjWe3JtTVndH/nv9eg/1b77xV9kbe5gwAALzHpNvLdPGkWwAA4HsEFj9hxVsAAHyHwHKZmhphYQE5AAB8h8BymZq6I8QCcgAA+A6Tbi9TUzNYuCPkqab2jOxhVoXZrPqmpk7ffFunOqdLX5/4VjW1Z2SzWuQyRjarVRZJtWdccpyql8sY1dQ6Ve906VS9U06XkdVikZGR02VUd8alcJtVlafqdareqTCrRU7X2bh4qs6pyHCr7OE2VZ0+o+NVtepkD5PVKkWG2+Qy0smaWnW2hys6wiaXMTpV71RkuE01tWckSWFWq8JsFtXWu2SzWhQedrZ9Rudv+xlzNqBaZFGYzSKXOXvMYrEo3GrRGS+G21o6Jaq23iWL5Wz50/UuWb471yKLvvuP+/u5cuc+S9LUjFTdOiixxe0CgGAjsFympm8Jtd/EUnvGqaOOWvXsEqWTNXWqOFWvVVsPqvTkKe07Vq2ocJuOVtWqV9dofVFepaNVte5zo8JtOlXf+CPiCJzMawgrANoWAstlsjaxnG17yCtnnC6VHKjQH97fq6Ivj6ve6d1FlZ78tsG+C8NKVLhN0RE2JXS2yxipvOq0BiR2lpEUZrUoJjJcFosUFx2uCJtVEWFW2axWWS1SvdMlm9Wq6AibjlfXqltnu6prneoeGymLRbJZLTpd75IxRjarRVHhNl1hD5PLGLmMUb3TyGaxKDrCppo6p5wul3sEwujsCEV0RJjOuFyqrXfJZYwiwqweoycWSVbr+ZELY86269wx810fhtmsTY/EXaL/mpq4bSTZw6zflZEiw89/Njo/4nP+u+dIkDHS0J5xl/hlAAg9BBY/aUsjLEerTstus8lxul4HTn6r6a99opo670dBesRF6apunRQTGaYu0RGyWS1KS4lV99goWSRdYQ9TZLhNV14RoZiocNl4dwEAoIUILH4S6nml5ECF7lj0316f1z+xk864jH4xdqD6deukssrTGt67i6wWiyLCmMMNAPCPVv2FWbRokVJTUxUZGan09HRt2bKlybI333zz2SHyi7YJEya4yxhjNH/+fHXv3l1RUVHKzMzUnj17WtO0oEjrGdtgX6gGlm9q6pQ6Z02LworFIq383xnanzdeXz07QV89O0F/n/0Dvf/zm5U1OElXJXTSjf3iFRluI6wAAPzK6xGWFStWKCcnR4sXL1Z6eroWLlyorKws7d69W926dWtQftWqVaqrq3N/P3HihNLS0nTXXXe59/32t7/VH/7wB7366qvq06eP5s2bp6ysLO3YsUORkZGtvLTAWTkjQ+Oe/0j7jtW494XCLaHDFac0b/V2dY4M0+qSw5csm3lNoob37qKpGb0VFW5rcm4OAADB4HVgyc/P1/Tp05WdnS1JWrx4sdasWaOlS5dqzpw5Dcp37drV4/vy5csVHR3tDizGGC1cuFBz587V7bffLkl67bXXlJiYqNWrV+vuu+/2+qICzR5mU3JslEdgCXZcMcbohmffv2SZdbNv0tWJnQPUIgAAWs+rwFJXV6fi4mLl5ua691mtVmVmZqqoqKhFdSxZskR33323rrjiCknS/v37VVZWpszMTHeZ2NhYpaenq6ioqNHAUltbq9ra84/KOhwOby4jIAI9wnLg5Lca/dsPdO/3e+v1TV9fsuynC25TbFR4gFoGAMDl8yqwHD9+XE6nU4mJnms4JCYmateuXc2ev2XLFm3fvl1Llixx7ysrK3PXcXGd545dLC8vT0899ZQ3Tfe7i5djCfS7hEb/9gNJajSsvDMzQ8ZIaSlxCrcx1wQA0PYE9CmhJUuWaMiQIRo5cuRl1ZObm6ucnBz3d4fDoZSUlMttnk8ZIxV/fVK7yqp0z8hesvj4rc4f7Tmme5c0Pdn5nC1P3KJunUN/HhAAAJfiVWCJj4+XzWZTeXm5x/7y8nIlJSVd8tyamhotX75cv/zlLz32nzuvvLxc3bt396hz2LBhjdZlt9tlt9u9aXrAuYz0w5fO3iZL6RKtm/on+KTeXWUO/euLH+vbJtZJuft7KfpfA7tpWK84nXEawgoAoF3w6v5ARESEhg8frsLCQvc+l8ulwsJCZWRkXPLct956S7W1tfrRj37ksb9Pnz5KSkryqNPhcGjz5s3N1hnKLpzDsrusqlV1fPLVSc38P8U6VHFKkvTZwQqNXfhRk2Fl67xb9ewPh+q2wUnq1jlSyXFRrfpdAABCjde3hHJycjRt2jSNGDFCI0eO1MKFC1VTU+N+amjq1Knq0aOH8vLyPM5bsmSJ7rjjDl155ZUe+y0Wix5++GH96le/0tVXX+1+rDk5OVl33HFH668syP79zW3uz3VOV6vq+H8Xnx2h+dv2Mv1h8nX62QV1nvPVsxPk/O7leqwcCwBor7wOLJMmTdKxY8c0f/58lZWVadiwYSooKHBPmi0tLZXV6jlws3v3bm3cuFF///vfG63zscceU01NjR544AFVVFRo1KhRKigoaBNrsDRl79Fqn9bXWFh57+GbJBFUAADtn8UE+nEWP3A4HIqNjVVlZaViYmKC0oZ7l2zWR3uON3rssbED9NOb+7W4rh2HHRr/h4+aPN6zS5R+88OhurFfvNftBAAgVHjz95tnXAPgjc2lHt+NMao70/htos8PVl4yrPzqjmu18Rf/i7ACAOhQCCwBcPCbU+7PZZWn1Sd3rYY+9Z4qT9Urd9Vneqf4oCTJcbpeE/9zY5P1zPt/BulH3+/t9/YCABBqeFuzj7R0nZV/ffHsSwdP17s0553P9LftZXpzywEt2bhfO440XLF34y/G6Juaeu084tBdI3r6tM0AALQVBBYfaW4q0KGKU+oRF6XDlafd+y78fHFYiY6w6YNHblZiTKR6dpGGNPJGaAAAOgpuCQXIT5b9T4N9nx6oaLL8jl+OVWJM231KCgAAXyKw+Ehzt4R2ebF43OzM/pfbHAAA2hUCSwA9vLzhWiqN+dktLX8EGgCAjoDAEkCrSw43Wyb/39J8/qJEAADaOibd+ojLdfnr77338E0akNTZB60BAKB9YYTFR+pb+b6gc6ak9yKsAADQBEZYfKRvwhXavP9kq879473DNWZgNx+3CACA9oPA4iNzxl4jq8Wi/++iZfhb4rbBSX5oEQAA7Qe3hHwkNjpcv75zSIvLx0WH+7E1AAC0LwQWH3v53uEe3x+4qW+j5dY/crPGDEjQmp+NCkSzAABo07gl5GNZF93e6Wxv2MVhVovioiP0SvbIQDULAIA2jREWP4sIa9jFUeG2ILQEAIC2i8DiZ40FlmG94gLfEAAA2jACi5+F2xp28XN3pQWhJQAAtF0EFj/4Qf8E9+eLnwb66c1XqRtvYQYAwCsEFj945LYB7s9dr4hwf340a4B+dsvVwWgSAABtGk8J+UGfhCvcn69L6aIFEwcpNf4KjRnAarYAALQGgcUPOtnD9Er292SRFBVhU/aNfYLdJAAA2jQCi58wmgIAgO8whwUAAIQ8AgsAAAh5BBYAABDyCCwAACDkEVgAAEDII7AAAICQR2ABAAAhj8ACAABCXqsCy6JFi5SamqrIyEilp6dry5YtlyxfUVGhBx98UN27d5fdblf//v21du1a9/Enn3xSFovFYxs4cGBrmgYAANohr1e6XbFihXJycrR48WKlp6dr4cKFysrK0u7du9WtW8PVXevq6nTrrbeqW7duevvtt9WjRw99/fXXiouL8yg3ePBg/eMf/zjfsDAW4QUAAGd5nQry8/M1ffp0ZWdnS5IWL16sNWvWaOnSpZozZ06D8kuXLtXJkyf18ccfKzw8XJKUmprasCFhYUpKSvK2OQAAoAPw6pZQXV2diouLlZmZeb4Cq1WZmZkqKipq9Jy//vWvysjI0IMPPqjExERde+21euaZZ+R0Oj3K7dmzR8nJyerbt6+mTJmi0tLSJttRW1srh8PhsQEAgPbLq8By/PhxOZ1OJSYmeuxPTExUWVlZo+fs27dPb7/9tpxOp9auXat58+bpueee069+9St3mfT0dC1btkwFBQV66aWXtH//fo0ePVpVVVWN1pmXl6fY2Fj3lpKS4s1lAACANsbvE0VcLpe6deumP/7xj7LZbBo+fLgOHTqk3/3ud1qwYIEkady4ce7yQ4cOVXp6unr37q2VK1fqJz/5SYM6c3NzlZOT4/5eWVmpXr16MdICAEAbcu7vtjGm2bJeBZb4+HjZbDaVl5d77C8vL29y/kn37t0VHh4um83m3nfNNdeorKxMdXV1ioiIaHBOXFyc+vfvr7179zZap91ul91ud38/d8GMtAAA0PZUVVUpNjb2kmW8CiwREREaPny4CgsLdccdd0g6O4JSWFioWbNmNXrOjTfeqDfeeEMul0tW69k7UF988YW6d+/eaFiRpOrqan355Ze69957W9Su5ORkHThwQJ07d5bFYvHmkprlcDiUkpKiAwcOKCYmxqd14zz6OTDo58ChrwODfg4Mf/WzMUZVVVVKTk5utqzXt4RycnI0bdo0jRgxQiNHjtTChQtVU1Pjfmpo6tSp6tGjh/Ly8iRJM2fO1H/+53/qoYce0r//+79rz549euaZZ/Szn/3MXecjjzyiiRMnqnfv3jp8+LAWLFggm82myZMnt6hNVqtVPXv29PZSvBITE8P/MQQA/RwY9HPg0NeBQT8Hhj/6ubmRlXO8DiyTJk3SsWPHNH/+fJWVlWnYsGEqKChwT8QtLS11j6RIZ2/TvPfee5o9e7aGDh2qHj166KGHHtIvfvELd5mDBw9q8uTJOnHihBISEjRq1Cht2rRJCQkJ3jYPAAC0QxbTkpkuHZjD4VBsbKwqKytJ735EPwcG/Rw49HVg0M+BEQr9zLuEmmG327VgwQKPSb7wPfo5MOjnwKGvA4N+DoxQ6GdGWAAAQMhjhAUAAIQ8AgsAAAh5BBYAABDyCCwAACDkEViasWjRIqWmpioyMlLp6enasmVLsJsUsj788ENNnDhRycnJslgsWr16tcdxY4zmz5+v7t27KyoqSpmZmdqzZ49HmZMnT2rKlCmKiYlRXFycfvKTn6i6utqjzGeffabRo0crMjJSKSkp+u1vf+vvSwspeXl5+t73vqfOnTurW7duuuOOO7R7926PMqdPn9aDDz6oK6+8Up06ddIPf/jDBq/UKC0t1YQJExQdHa1u3brp0Ucf1ZkzZzzKrF+/Xtdff73sdrv69eunZcuW+fvyQsZLL72koUOHuhfKysjI0N/+9jf3cfrYP5599llZLBY9/PDD7n309eV78sknZbFYPLaBAwe6j7eJPjZo0vLly01ERIRZunSp+ec//2mmT59u4uLiTHl5ebCbFpLWrl1rnnjiCbNq1Sojybz77rsex5999lkTGxtrVq9ebT799FPzL//yL6ZPnz7m1KlT7jJjx441aWlpZtOmTeajjz4y/fr1M5MnT3Yfr6ysNImJiWbKlClm+/bt5s033zRRUVHm5ZdfDtRlBl1WVpZ55ZVXzPbt201JSYkZP3686dWrl6murnaXmTFjhklJSTGFhYXmk08+Md///vfNDTfc4D5+5swZc+2115rMzEyzbds2s3btWhMfH29yc3PdZfbt22eio6NNTk6O2bFjh3nhhReMzWYzBQUFAb3eYPnrX/9q1qxZY7744guze/du8/jjj5vw8HCzfft2Ywx97A9btmwxqampZujQoeahhx5y76evL9+CBQvM4MGDzZEjR9zbsWPH3MfbQh8TWC5h5MiR5sEHH3R/dzqdJjk52eTl5QWxVW3DxYHF5XKZpKQk87vf/c69r6KiwtjtdvPmm28aY4zZsWOHkWT+53/+x13mb3/7m7FYLObQoUPGGGNefPFF06VLF1NbW+su84tf/MIMGDDAz1cUuo4ePWokmQ0bNhhjzvZreHi4eeutt9xldu7caSSZoqIiY8zZcGm1Wk1ZWZm7zEsvvWRiYmLcffvYY4+ZwYMHe/zWpEmTTFZWlr8vKWR16dLF/PnPf6aP/aCqqspcffXVZt26deYHP/iBO7DQ176xYMECk5aW1uixttLH3BJqQl1dnYqLi5WZmeneZ7ValZmZqaKioiC2rG3av3+/ysrKPPozNjZW6enp7v4sKipSXFycRowY4S6TmZkpq9WqzZs3u8vcdNNNHi/OzMrK0u7du/XNN98E6GpCS2VlpSSpa9eukqTi4mLV19d79PXAgQPVq1cvj74eMmSI+5Ua0tl+dDgc+uc//+kuc2Ed58p0xH//TqdTy5cvV01NjTIyMuhjP3jwwQc1YcKEBv1BX/vOnj17lJycrL59+2rKlCkqLS2V1Hb6mMDShOPHj8vpdHr8jyNJiYmJKisrC1Kr2q5zfXap/iwrK1O3bt08joeFhalr164eZRqr48Lf6EhcLpcefvhh3Xjjjbr22mslne2HiIgIxcXFeZS9uK+b68emyjgcDp06dcoflxNyPv/8c3Xq1El2u10zZszQu+++q0GDBtHHPrZ8+XJt3brV/dLcC9HXvpGenq5ly5apoKBAL730kvbv36/Ro0erqqqqzfSx1y8/BBA6HnzwQW3fvl0bN24MdlPapQEDBqikpESVlZV6++23NW3aNG3YsCHYzWpXDhw4oIceekjr1q1TZGRksJvTbo0bN879eejQoUpPT1fv3r21cuVKRUVFBbFlLccISxPi4+Nls9kazJIuLy9XUlJSkFrVdp3rs0v1Z1JSko4ePepx/MyZMzp58qRHmcbquPA3OopZs2bpv/7rv/TBBx+oZ8+e7v1JSUmqq6tTRUWFR/mL+7q5fmyqTExMTJv5f3CXKyIiQv369dPw4cOVl5entLQ0Pf/88/SxDxUXF+vo0aO6/vrrFRYWprCwMG3YsEF/+MMfFBYWpsTERPraD+Li4tS/f3/t3bu3zfx7JrA0ISIiQsOHD1dhYaF7n8vlUmFhoTIyMoLYsrapT58+SkpK8uhPh8OhzZs3u/szIyNDFRUVKi4udpd5//335XK5lJ6e7i7z4Ycfqr6+3l1m3bp1GjBggLp06RKgqwkuY4xmzZqld999V++//7769OnjcXz48OEKDw/36Ovdu3ertLTUo68///xzj4C4bt06xcTEaNCgQe4yF9ZxrkxH/vfvcrlUW1tLH/vQLbfcos8//1wlJSXubcSIEZoyZYr7M33te9XV1fryyy/VvXv3tvPv2SdTd9up5cuXG7vdbpYtW2Z27NhhHnjgARMXF+cxSxrnVVVVmW3btplt27YZSSY/P99s27bNfP3118aYs481x8XFmb/85S/ms88+M7fffnujjzVfd911ZvPmzWbjxo3m6quv9nisuaKiwiQmJpp7773XbN++3SxfvtxER0d3qMeaZ86caWJjY8369es9HlH89ttv3WVmzJhhevXqZd5//33zySefmIyMDJORkeE+fu4Rxdtuu82UlJSYgoICk5CQ0Ogjio8++qjZuXOnWbRoUYd6DHTOnDlmw4YNZv/+/eazzz4zc+bMMRaLxfz97383xtDH/nThU0LG0Ne+8POf/9ysX7/e7N+/3/z3f/+3yczMNPHx8ebo0aPGmLbRxwSWZrzwwgumV69eJiIiwowcOdJs2rQp2E0KWR988IGR1GCbNm2aMebso83z5s0ziYmJxm63m1tuucXs3r3bo44TJ06YyZMnm06dOpmYmBiTnZ1tqqqqPMp8+umnZtSoUcZut5sePXqYZ599NlCXGBIa62NJ5pVXXnGXOXXqlPnpT39qunTpYqKjo82dd95pjhw54lHPV199ZcaNG2eioqJMfHy8+fnPf27q6+s9ynzwwQdm2LBhJiIiwvTt29fjN9q7H//4x6Z3794mIiLCJCQkmFtuucUdVoyhj/3p4sBCX1++SZMmme7du5uIiAjTo0cPM2nSJLN371738bbQxxZjjPHNWA0AAIB/MIcFAACEPAILAAAIeQQWAAAQ8ggsAAAg5BFYAABAyCOwAACAkEdgAQAAIY/AAgAAQh6BBQAAhDwCCwAACHkEFgAAEPIILAAAIOT9//Ynkd8e19OBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p0 = [-dall.volume_m3.max(), 1, dall.volume_m3.max(),1,1]\n",
    "fittedParameters, pcov = curve_fit(func, dall.simulation_year.values, dall.volume_m3.values, p0 = p0)\n",
    "t = np.arange(0,5001,1)\n",
    "d_fit_volume_m3 = func(t, *fittedParameters) \n",
    "\n",
    "plt.plot(dall.volume_m3)\n",
    "\n",
    "plt.plot(t,d_fit_volume_m3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 294,
   "id": "70f27c47",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2.25290996e+12  1.32947258e-02  6.59375456e+12  7.45563511e-01\n",
      "  1.38590735e+00]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/users/lschuster/.local/lib/python3.10/site-packages/scipy/optimize/_minpack_py.py:881: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fce52025bd0>]"
      ]
     },
     "execution_count": 294,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGsCAYAAAD+L/ysAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRxUlEQVR4nO3deXxU9b3/8deZmWSy72QhBNlBloCAYqxWK5FFqlB3f1jUulyt9uptr1psa0u1N1S996rVolUreltMqxVs3dCiQa2IgCCLiuxhC4FA9mSSmfn+/ggZMiQhe2aSvJ+Px3nkLN8585kT47z5nnO+xzLGGERERESCmC3QBYiIiIi0RIFFREREgp4Ci4iIiAQ9BRYREREJegosIiIiEvQUWERERCToKbCIiIhI0FNgERERkaCnwCIiIiJBT4FFREREgl6vCywffvghl1xyCf3798eyLJYtW9am11dXV3PDDTcwbtw4HA4Hc+bMadTmtdde46KLLqJfv37ExMSQlZXF8uXLO+cDiIiISCO9LrBUVFQwfvx4nnrqqXa93uPxEB4ezr//+7+TnZ3dZJsPP/yQiy66iLfeeot169bxne98h0suuYT169d3pHQRERFphtWbH35oWRZLly716yVxuVz87Gc/4+WXX6a4uJixY8fy29/+lgsuuKDR62+44QaKi4tb1UszZswYrr76ah544IHO+wAiIiIC9MIelpbceeedrFq1itzcXDZu3MiVV17JjBkz2LZtW7v36fV6KSsrIyEhoRMrFRERkXp9KrDk5+fzwgsv8Morr3DeeecxdOhQ/vM//5Nzzz2XF154od37ffTRRykvL+eqq67qxGpFRESkniPQBXSnTZs24fF4GDFihN96l8tFYmJiu/a5ZMkSFixYwOuvv05ycnJnlCkiIiIn6VOBpby8HLvdzrp167Db7X7boqKi2ry/3Nxcbr75Zl555ZVmL9AVERGRjutTgeWMM87A4/FQWFjIeeed16F9vfzyy/zgBz8gNzeXWbNmdVKFIiIi0pReF1jKy8vZvn27b3nXrl1s2LCBhIQERowYwdy5c5k3bx7//d//zRlnnMHhw4dZsWIFmZmZvuDx5ZdfUlNTw9GjRykrK2PDhg0ATJgwAag7DXT99dfz+OOPM2XKFAoKCgAIDw8nNja2Wz+viIhIX9DrbmvOy8vjO9/5TqP1119/PYsXL6a2tpaHHnqIl156if3795OUlMTZZ5/NggULGDduHACDBg1iz549jfZRf6guuOACVq5c2ex7iIiISOfqdYFFREREep8+dVuziIiI9EwKLCIiIhL0esVFt16vlwMHDhAdHY1lWYEuR0RERFrBGENZWRn9+/fHZjt1H0qvCCwHDhwgIyMj0GWIiIhIO+zdu5cBAwacsk2vCCzR0dFA3QeOiYkJcDUiIiLSGqWlpWRkZPi+x0+lVwSW+tNAMTExCiwiIiI9TGsu59BFtyIiIhL0FFhEREQk6CmwiIiISNBTYBEREZGgp8AiIiIiQU+BRURERIKeAouIiIgEPQUWERERCXoKLCIiIhL0FFhEREQk6CmwiIiISNBTYBEREZGg1ysefigiLXv6X2t5fdMW5s8cR2yEHWMMofZQKmoriHfGU1RdREJYAqU1pYQ5wnB73QCE2EKorK0k1hlLUXURiWGJlLhKCA8J97Vx2BxUu6uJCY2hqLqIpLAkil3FRIREUOOpwbIsHJYDl8dFVGgUx6qPkRiWyDHXMaJConB5XNgsG3bL3mIbm2Wj1ltLhCOCYlexr010aDTV7mrslh3LsnxtSlwlJIQlcLT6KLHOWCrdlTisuv/1eYyHcEd4k21CbCEYY/AYD2GOMMpqynzHKc4Z52vjNV68xtuoTXxYPOU15TjtTjzGg9d4cdqdlNeWN9nGbdy+30llbSVxzrg2/04iQiKo9dZ2y++k/nj39t9Jw+N9qt+Jw+agqraq1b+To9VHSQxL9PudtOVvoMZT0+zv5OTj3Rm/k4SwBAbGDOyO/1U1yzLGmIBW0AlKS0uJjY2lpKRET2sWacJf129m/ptvA5AaG8Zv5owLcEUi0tOclXoWESERnbrPtnx/65SQSA9mjOHZVZ/z/je7Kat2NduuPqwAFJRUd0dpItLL1PcUBUqHAsvChQuxLIu777672TYXXHABlmU1mmbNmuVrc8MNNzTaPmPGjI6UJtKjfVN4lCXrNtFSB+g/Nn/DwhUruOWvrzDlf3+Px+Nt1KYXdKKKSBCwLCug79/ua1jWrFnDM888Q2Zm5inbvfbaa9TU1PiWi4qKGD9+PFdeeaVfuxkzZvDCCy/4lp1OZ3tLE+nxZv7heQA8XsP3z2z+b2zTwUO+eZfHTYmrmoQI/y7bare7a4oUaS2vFzBgDJj6eW+D+QbbjLdu3tfGAF7/NpzU1q99U6+vb8OJ5bqFBtupe5+GbQz+r23UvsG+GtZBM69v8j0a7LvhsfBbRxP7O0WNPg0/M43naWq9aaa9oTLmfWKyf02gtCuwlJeXM3fuXJ599lkeeuihU7ZNSEjwW87NzSUiIqJRYHE6naSmpranHJEe6cY/L6XM5eKVG69u9l8un+zc02xgufPVN3n76y/91n2+t4DskUP81lXWNO7Grar1EB5ib2fl0i5eLxgPeN3Hpybm67d7POCtrfvC9dae1N7TxGuO79t4Gyx769r41nv923q9QH2bhu0btvXUfWGdaj9NBYaGYUE9fL3G1zFppPa0wHLHHXcwa9YssrOzWwwsJ3v++ee55ppriIyM9Fufl5dHcnIy8fHxXHjhhTz00EMkJiY2uQ+Xy4XLdeJ8fWlpads/hEgAVdbU8uGu7QBsO3yMEckJTbar9XpYt/cg//VeHvdN/TZnnZbu23ZyWAH496Wvs2DGdK6cMBqATQcO8aO//aNRu22HysgcENcJn6SH8XrAXQ2emrqfblfdfG318aBQA57aE5O34fLxee/xqan1vmV3g7b1waLx6TppgmUDmw2w6uYt6/j88WWOr6tfPmUbG1jUbcM6Me/7B0LD1x2f56Q2vp+2ZrY1fM1J71c/S8P9N/F6vxobrjt5fQt1NPqHT4Nl6+R5q3G7hsfl5PcHcEYTSG0OLLm5uXz++eesWbOmzW/22WefsXnzZp5//nm/9TNmzOCyyy5j8ODB7Nixg/vvv5+ZM2eyatUq7PbG/wrMyclhwYIFbX5/kWBRXHXiwtdLn1vMknnXMnFAGgAFpeW+bR/s2MYHO7YBcO3/LeHl7/8//r75az7euavJ/brcbn76xpv89I03yUxLZ+PB/U22C7H3gOvtPW5wV0FtFdRWHv95fHI3WHa7/CdP/Xx9MGkw7wmy02N2B9gcYNnrftocYGtpvsE6y1735W5zALa6bfVf+Jb9+NRwuUEby3bi9U3NW3awWSftp8Frbfbj79kwMNTPNwgM9ettLbVp6gtX5IQ2BZa9e/dy11138d577xEWFtbmN3v++ecZN24cZ511lt/6a665xjc/btw4MjMzGTp0KHl5eUydOrXRfubPn8+Pf/xj33JpaSkZGRltrkekq5RUVRMb3vzfSFFFpW++1uvhP5a+wR+u+h5HK6v45Tv/bPZ11/7fklbX0FxYAXC5u+Ff+14P1FRCbTm4KqCmvG6qn6+tAFc51FT4B5L6kOLpwjsSLAscYWAPBYezbrKHHp9C6iabo27ZFnJiXf1220nLTa2zhZwUSEKOh4uQBl/8+oIWaa02BZZ169ZRWFjIxIkTfes8Hg8ffvghTz75JC6Xq8keEYCKigpyc3P59a9bPv81ZMgQkpKS2L59e5OBxel06qJcCUpl1S7+uHo9T3z0Ef92dhb3Zp/bZLvHV67yWz5aWcnFz77QZNvONigxktMS2jGWgrsaqkvrJlcJVJccny89Pl9yPJRU1IWS2sqW99kaDieEhPtPjojjP8MgJOx44HBCyPGfvhBSP98gnNiddYFCYUGkTUbGjwzo+7cpsEydOpVNmzb5rbvxxhsZNWoU9913X7NhBeCVV17B5XJx3XXXtfg++/bto6ioiLS0tLaUJxJQxhhmPLOYgrK6a6qe+XQVI1KSmDNuVKO2n+3N91uurK1p1Kar/OK7o08sGFPXy1FVBFXFUHkUqo5C1bG6qbrkeCAprTu10h4h4RAaWTc5oyA0CkIazIdGQEj9dHIwCT9+6kFEAmlo3FDSogL7ndymwBIdHc3YsWP91kVGRpKYmOhbP2/ePNLT08nJyfFr9/zzzzNnzpxGF9KWl5ezYMECLr/8clJTU9mxYwf33nsvw4YNY/r06e35TCLdoqi8koTIcN8dPhWuWl9YqfeT1/9BfHgY5w8b5Lf+snHj+L91a7u8xghcJFqlJFFColXCz0amEPLJeqLLD1NeXkBk2SEqjRunMXgBjwVOY6iwbER7vZTZbER5vVTabIQZg8fhxBORhDMigYqIBKIjkykLiyEqKo3KsGjCIvrhDovGOGMIjehHpcNBVFg8ZTVlRIdGU15bTrgj3DcAVf0w9FEhUc22CbGFUOWuaraNhVU3LHoTbU4e8ry5Ni63C7vNjt2yU+2uJiIkgvKacr821e5qQmwhWFi4PK5TtgGo8dQ0ahMZEkllbSVOuxODodZbS7gj3NemrKaMqNAoKmsrCXOE4TEePF4PToeTipqKFts0PJYN27i9bgymU453d/9ObNhaPN6t/Z2cfLxP/p148eL2ulv1OwlzhDXbprnjHcy/k4Z/A1XuKiJDIv2OU3JEcpf//6olHR6a/4ILLmDChAk89thjvuVBgwaxePFiX5utW7cyatQo3n33XS666CK/11dVVTFnzhzWr19PcXEx/fv3Z9q0aTz44IOkpKS0qgYNzS/dbdWuvVz351zOGzyMxXO/B8Cf127igXfeadT2tqxz+PEFWdjtNmrdHt7dupONBwp4bvWnHa4jBA/9rGKSKSbRKiHpeDhJskpItEqJxL9X5IKqqqZ3FJEE0WkQk1b3MzoNolMhKhki+0FkUt3P0CidShGRTtOW7289S0ikHcY//ATlNXVhYMfP7wFgzMLHqHY3faHo2JT+vH7LXM594g8cLC1p03vZMCRYZaRwlBTrGCnWMVI5Rop1lESrtMXhqssI54iJITttJBmpwyEm/Xgw6V8XSqJT667tEBHpZm35/tbTmkXaoT6sNBQXFk5BedOBZfOhAwAthpVESkm3jjDAOkL/4z9TraOE4Gn2NVWEUmjiKDSxFBHLERNLkYnhCHU/XYSw6IrLyBg1tA2fUEQkuCiwiAD/t2Yj6bExXDhi0CnbeTxeymqavkB2ZHIyBeXND2LYsDMzFDcZViGnWYcYYB0m3TpCulVEOE3vuxYbhSaOQyaBQ8RTYOLr5k08pbR8x094SEiLbUREgpkCi/R5//PBKp7618dA3ekdt8fL0coqth0+yuubvuJn077tG1PlysW5fNHM+CZVtU33roTgZqB1mFtyruNGewGnWYfobxVhp/HZWDc2CkwCOAewuiqKAyaJ/SRy1dkzmDpqOJMdDpZu/Jqt+XvYdvCA73UPzpzBL95ufP0MQL/IKM5I12MvRKRnU2CRPq8+rAB8uD2fG3P/4rfdGMMjc6bj9nibDCuX//FlLh8/lu1HjgAQTxnDrAMMs+1nuLWfAdbhJsNJCeHs8aaSTzL7vUnsox+FJh43Nn767Qu5MDKSn7/9Dr+cns1VZ5y4O++nF50LnBjfpai8ksSoCL/A8vvLL+O1jVtYeMlFxIWHBfwpqyIiHaXAIn1erDOcElfd3TMnhxWArwoPA7DjyNFG2ywMRw5u4P2CN7jC2s/QkAMk0fi0UCnh7PamsIcU9nhT2W1SOEY0I5KS+eZIYaP2kc5Qvjd+FJeOHYG9hWH0E6PqTgmd0T+D9Qf2ct7gYUw/fSjTT9c1KyLSeyiwSJ/37aFD+ceXm5vdXuNx89me/fzob38HIMkq5XRrT91kyycG/1uFvUC+SWa7SWe7N52dpj9FNP3QsIcvncmcP77YaH1EaN01Jy2FlYYe+95M/rbxK34w5YxWv0ZEpKdQYJFexxjT7CmQyppath0uYnyDazpqW3ggXm11KW+89RQXu9ZyumMPyZb/nT4uHGz3prPNpLPdpBOfdiYRkXHcdPakZp/9c9nYTO66IIsBcTGsvvuHTHns937bwx1tv0h2QHwsd51/dptfJyLSEyiwSK9S6/Zw8R9eIi0mhpeuu7zR9nEPPwbAHd86l0kZafwg95Um95NklTDe2sl4awcjavYSUuOlfsATDxa7TBpfmtP4yjuQnaY/7gajobzz3e8yPDmB5oY4GhiXwIOzphIWUvfnlxQVyZSM01i9d4+vzZDE+PZ8fBGRXkuBRXosYwzvfr2TiNAQbnj5L9zzne8wOiWZnUePsPPokVP2tDz1r49JiIj0WzfQKmSybSsTrB2kW0V+2wpNHBvNYLZ4B/GNGUA1oc3W5bDXvWdz7/3tIYN9YaXekuuvAuDJD1djt9kZnpxw6g8vItLHKLBIj/Xo+5/w9KpPfMuPfPABL1xztW/5WFUVCRHNj1FytLKCAdYRzrR9zWRrK6lWsW+bB4vtJp0vvEPYaIZy0LQ+QNS4mx/kDSDU0fyf3Z3fntLq9xER6UsUWKTHahhW6jkdJ57su/3wMc46rXFg6WcVc47tS860vibNOuZbX4OdTd4hrPMOZ7MZTAVhp3z/hPBIjlZVNFo/IO7E8NKPXnIJ8998i1pvgxDT85+GISLS7RRYpMe6cvwEXvlig986V4PejZKqBsPnu8rhy9e51/EXRlr7fKtrj4eUNd4RbDRDTnmq52Q/m3YhReVV/NeKf/rWDU9KJtJ5Yh/fGz+KmaOHMea3/+tbd6i8vNXvISIidRRYJGg1vAbl8ZWfsmr3Hp664hK+KjjCl4cKCWvi1EpFg2Hzq2trYO8aWPcCbFnGYW8VI0OdeIGvvKfxiXcMX5ghVOH/4L/IUCcVTTwr6GRbDhbys2nf9gssdlvj61ZOvl6lf2xsi/sWERF/CiwSlPYdK2HWsy8yY9Tp/PbSi3jio48AeOrDNby+ZTPF1ZUkhkc2et2OI0WE4GaK7WvGvP83OPYVNUCpzcb7oSl87BnLp97RHG1mXBSAu847zy+ENOR0OHC5626D/tbggQB89h93cNb/PgWArcVnJ8Nt35rcYhsREfHX+lGpRLrJPzZv5ft/fpXyGhevbtzgt624uori6koAik66fiSRUjK2PsOjIc9wo305kaXfgN3JPxIv4Dr79dzvvom3vFOaDSuPzbmUlXf+m981KNdNnOTXZv7UC33zkzP6171v5InrZE5LaPp25InpGQCE2OzEhZ/62hgREWlMPSwSVApKy7l72d/91jUcz8TbxAWrA6wjzLB9xlm2r7Efrdt+hBi+HDaH15Mv4rerNrT4vrecncUlY0cCkB4bzZXjJzCiXxKpMVH86fN1vnbpDU7nRDpPDO72+Pfm8MqGjfxqxgVN7n/J969k+5FjnJ6a1GItIiLSmAKLBJVjldWN1lXVnhiJ1u31+uaHW/uZafuM8badvnVfegeywjuRjWYI3h0W7NhwyvezYfGDKVO4u8EIsZZlsfCSiwD4ZOde3/opA08jNTrSr129744ZznfHDG/2fUIcdoUVEZEOUGCRoOKwNT5LWe2u9c0XlJYx0trHbPvHjLTqnpzsBdZ5R/C29yz2mJQ2vV/WoMHMv+i8ZrdPOS3dN3/O4EGMTuvHr6ZPIyUqqk3vIyIiHaPAIkGlqRFK1u8tAGCwdZDvHHqVMY66IexrsfGJdyzLvZM5ZFo/lP3W+T9h8Wcb+HR3Pr+a8Z1TtvV7+ODx01HfP3N8q99LREQ6hwKLdIvCsgoeX/kpcydlMjqtX7Pt3J7Go8Q+8Y8X+ZH9n0yw7ahrg42PvON4w3M2xbS9p8Nht3Fz1kRuzprY5teKiEhgKLBIt/jJ62/zye5dvL55E2/deiMPv/8xN589iQkDUv3afVV4xDcfTRWz7f/ifM8X2Gx1p34+8Y7hH94sjpi6i19/c/FMXt/0JZ81eHDgqQxNbD4sNcdu2fAYLwPiNH6KiEigKLBIt9hSUHdap8pdy91L3+SLg/v55zdb+dddt/F/a77g2kmZpERHcrSiEgdeLrSt5xL7KiKoG8BtnXc4S73n+j3T51fTp3PNxLFcM3Es/9q5l3lLclusIzwkpMU2J3v00u+y7XARl4wZ0ebXiohI51BgkS6Xt203JdVVvuUvDtZdLFvr9XDP68tZuXM7f9nwBavuvo3Qg5/yK8di3zN+8k0/cj0XstUM8NtnZlo63z8z07d8emqi3/ahif3YUXS4US1VtTWN1rXk0nEj2/waERHpXAos0uVu+ssrzW5buXM7AFXlhXj+dgunbf8HWFBCOK95zuMT71i8NB7u/uS7iRr2nESFOpsMKwCVtbVNrhcRkeCmwCIBd461hasceWz4phKvzcYH3gks9Zzb6Bk/Z/TPYP2BunFRTg4sDZ8rZLMaB5x61QosIiI9kgKLBEwc5dxoX85Y224AvrRSeMl9ETtNWpPtNxcc8M2Xu/xP7TQcxM2yLEYkJfPNkcJG+7j57CmdULmIiHQ3PUtIusy+YyVc8+Jfmtx2prWVX4csZqxtNzXYecXzbR50X+cXVv5w5eV+rxnRL9k3v62JMFLPbtm4cPhQv3U2LN64+Ub+7ZxJzbxKRESCmXpYpMvk/PNj1uzN91sXQTXX2VcwxfY1ALtNCs95Lva7+wfg7IGDyRqc4bdu0ZWz+faTTwN1F+w2p6LGxaVjR/H0qlUMjEvgD1d9j/gIJ0lRjZ/uLCIiPYMCi3TYG1u28cmuPfx65oU4GowMu/XwIb92Q6yD/JvjDZIoxYPFm96zecNzNp4mOvpemnsZNpv/tSjpcU0/ZflkLo+bkSlJ5N1xK4mREUSEtv1WZhERCS4KLNJhdy1dBsDgxARuaTB67K6jRb75bNs6rrR/iAMvhSaOZz0XN3utyoDYOP8h8dspI14DvYmI9Ba6hkU6pKrBXTf5R4sbbY/AxQ/tf+daex4OvKzxjuDXJ12r0mifNc3fyfO/sy/FaXfwq+nTGm373WVzAMidN7fV9YuISM+gHhbpkI37T1z8WnvSc4D6W8e40/4aKVYxbmz8xXMB73vPaHGfRVUVjdZNHjAQqBvE7eLRw/1OPdW7ePRwLh59T1s/goiI9ADqYZEO+erQicDyr927mLboBVbt2svutUuZ7/gzKVYxR4hhofvaRmHlxsln+eb/ev11Te4/JSoGgO+OGeVb11RYERGR3k09LNIhj36w0jd/oLQEgMUv38PV9jwigG9MOr93z8YZ0Y8bR4/mhbWf+dqP7Z/Cwu/Ool9kBPERYU3uf9lN17FhXwEXjRrSlR9DRESCnAKLdEh4SChV7rprTux4+X/2FVxg2wjAx96x/MmTTS125o7P5D8vPIf/d2YmN/z5VfaXFnPe0IEkRkYAkH+0xLfP/7p4pm8+OTqSaaf7j6kiIiJ9jwKLdMighASO7q8gBDe32t9kom07XuAVz/m8653sa+cxBsuyGJIYz4of/oCK2lriwk/0qjQcV2XaKAUUERHxp8AiHVJV6yYCF3c6ljLS2k8tNv7g/i6fm+F+7YYknhgYLsRhJ85h99ue3GBQt2in/zOEREREFFik3TweLwcKt3Ov429kWIepJJSn3HP42viPUDtj5GiunHD6KfcVHebk7VtvwmG3dFGtiIg0osAi7fJ/a77g8eV/4z7HX0mxiikmksfcl7PX9GvU9qkrZ7VqnyOSE1puJCIifZICi7TK1kNHsCybL1Q8sfxv3Ov4K8lWMUdMDI96ruSwiQtskSIi0mt1qO994cKFWJbF3Xff3WybxYsXY1mW3xQW5n8LqzGGBx54gLS0NMLDw8nOzmbbtm0dKU06UXWtm6teXMLMPzxPUXklFO/lHsdfSLaKOWxieNh9NcVWYpOvXXnHrd1crYiI9EbtDixr1qzhmWeeITMzs8W2MTExHDx40Dft2bPHb/vDDz/ME088wdNPP83q1auJjIxk+vTpVFdXt7c86USl1S7Ka1wA7Nq9hdI/XkyyVUKhieVh9zUUEUOIrfF/Sgu/O4sBep6PiIh0gnYFlvLycubOncuzzz5LfHx8i+0tyyI1NdU3paSk+LYZY3jsscf4+c9/zuzZs8nMzOSll17iwIEDLFu2rD3lSSdzud0ARFFF0Zs38nntYQpNHA+7r+YodU9QtjURWMpdrm6tU0REeq92BZY77riDWbNmkZ2d3ar25eXlnHbaaWRkZDB79my2bNni27Zr1y4KCgr89hUbG8uUKVNYtWpVk/tzuVyUlpb6TdJ1qmo9hFHD3Y6/4eQIR4niUfdVHDseVgBfD0xDFad4iKGIiEhbtDmw5Obm8vnnn5OTk9Oq9iNHjuSPf/wjr7/+On/605/wer2cc8457Nu3D4CCggIAv16X+uX6bSfLyckhNjbWN2VkZDTZTjpHTXU5d9pfZ7B1iDLC+R/3lRQ1CCv1DyY82cWjhze5XkREpK3adJfQ3r17ueuuu3jvvfcaXTjbnKysLLKysnzL55xzDqeffjrPPPMMDz74YNuqPW7+/Pn8+Mc/9i2XlpYqtHQVr5eED37C6bZ8qgjhMfflHDQnbj9+59abGJ6cwNCHHvGte/PmG4mLCCM1JioQFYuISC/UpsCybt06CgsLmThxom+dx+Phww8/5Mknn8TlcmG320+xBwgJCeGMM85g+/btAKSmpgJw6NAh0tLSfO0OHTrEhAkTmtyH0+nEqdFQu5wxhqI35vNNwUfUYuNJ9xx2mxM9YU9fcTnDTxo7JSIklFGpSd1dqoiI9HJtOiU0depUNm3axIYNG3zT5MmTmTt3Lhs2bGgxrEBdwNm0aZMvnAwePJjU1FRWrFjha1NaWsrq1av9emaka9S6Pew7VkJ1rdtv/b927uWGnOvZ/NWLACx2z+Br43/qx+k4kXcfmjmTgXEJLPn+NV1ftIiI9Dlt6mGJjo5m7NixfusiIyNJTEz0rZ83bx7p6em+a1x+/etfc/bZZzNs2DCKi4t55JFH2LNnDzfffDOAbxyXhx56iOHDhzN48GB+8Ytf0L9/f+bMmdMJH1FO5eoX/8IXB/cDsHX+T3DYbfzsjRVs+GIZ/+H4JwB/92bxqWk8tP6YtBM9KddOGsu1k8Y2aiMiItIZOn2k2/z8fL9bXI8dO8Ytt9xCQUEB8fHxTJo0iU8++YTRo0f72tx7771UVFRw6623UlxczLnnnss777zT6utkpP3qwwrAO1/v4KuCQ3zwxQp+5vg7Drys9o7idc85jV53+bjxJEZGdGepIiLSh1nGGBPoIjqqtLSU2NhYSkpKiImJCXQ5PUrDi2UBnNRyv2MJA6wjbDdpPOq+itqTcm1SRBSrf3x7d5YpIiK9UFu+v/VYXPEzz/4uA6wjFBPJ792zG4UVgCOV5QGoTERE+jIFFvGZalvP2bav8WDxtPsSSogMdEkiIiKAAkufdLSykvozgTar7j+BYdZ+rrLnAfCK53y2mfQAVSciItKYAksfk7dtN2f971Pc9/f3+KrgCF7jJYJqbnW8iQMvn3lH8p53UqDLFBER8aPA0sc8/uEnGAN/2/QF333uBQDm2f9JImUUmjhe9Exr9rX/O+dSwh0h/Gp6821ERES6Qqff1izBLT483G/5XNtmzrRtxY2NP3gupprQJl83Pi2dS8eOZNbpw7HblXNFRKR7KbD0ITe/vIyVO7f7llOsY1xrrxth+HXPt9hl0hq9Zu2P7+S9rTu46oy6QeEUVkREJBD07dMLVde6+dfOvbg9Xt+6p/+1lg92bPMt2/Fys/1NwnDzlTeDd7xnNtrPzFGjiY8I94UVERGRQFEPSy9021//wUe7tnP95LN4YMb5ADy7arVfm2m2tQyxDlGBk+c9M/FiMX3EKP7nezMJC3HwTeFRhibGBaB6ERGRxhRYeqGPdtWd9nlx7WecMziD84YOpLi60rc9zTrKbPu/AMh1f4dRp2Xyo/OymDgglRBH3QMsR5z0FGYREZFAUmDp5f7tlb8xZeBpvmUbhhvt7xCCl43eQXxixvDuzAsZmqSAIiIiwUvXsPQBq/P3+OazbesYah2kklBeOn4Ls8KKiIgEOwWWPiTJKuF7x08F/dVzAceIJi1aD4sUEZHgp8DSy/z7395qdtu1tvcJPX5X0EfecQC4vd5m24uIiAQLBZZepKi8kje/2tLktgnWDibYdnJ6rZc/e7N962PCwrqrPBERkXZTYOlF3McfaHiyUNxc63gfgGMjr6N/+gTftsfmfLc7ShMREekQ3SXUizQcKK6hWfZPSaKUI8QQknkrLw0ZxpaCw0wckIplWd1cpYiISNuph6UXOVha3mhdinWMGbY1AOS6L8DriCAsxMGkjDSFFRER6THUw9KLPPXRp43WXWH7CAdeSkNG4U78NlmDBgSgMhERkY5RYOlF1u7L91seYe1jom0bIQYuvW4RlyafHqDKREREOkanhHqJbYVHqayt8S1bGK605wGQNvC7oLAiIiI9mHpYeon8Y8V+y2dZXzPEOsRgj53TLv2vwBQlIiLSSRRYeomwkBDf/MDoSB6rXUVVdS3p3/oPiEoOYGUiIiIdp1NCvYTTcSJ7/i2zmuTyA5wWkYLjWz8KYFUiIiKdQ4Gll/AcH2I/NTyEhLW/r1t5/r0QGhHAqkRERDqHTgn1Aj9YspSVO7cDcL53DVQchvhBMGFuYAsTERHpJOph6eEqXDW+sBKOi3O8dU9j5vyfgj3kFK8UERHpOdTD0oM9suJfmAbPD8q2fU4U1ZA4HDKvCmBlIiIinUuBpYf6quAIT6/6xLccQTXT7OsYWlsDF/wUbPYAViciItK5dEqoh/rhq6/7LV9o20AELjyRg2HMZQGqSkREpGsosPRQ+cVHffOhuMm2fw7AvhHzwKZfq4iI9C76ZuthdhYd42hlpd+682ybiKaKwyaGL6LODFBlIiIiXUfXsPQgG/YVcOXiP+HlxIW2drxMs68FYLn3TBZOGh+o8kRERLqMAksPsuKbnX5hBWCK9RVJlFJuInj2P5/F0kBxIiLSC+mUUA/xyoYv+eLAQb91FoYZ9s8AWO6dpLAiIiK9lnpYeoBNBw7x0zfebLR+rLWbdOsolYSS59WpIBER6b3Uw9ID7Cg61uT6qba6O4M+9o6jEmd3liQiItKtFFh6AJtlNVqXZh1lnG03XuB97wT+ev113V+YiIhIN1Fg6QHsVuNf04W29QBs9A7lsIkj2qkeFhER6b10DUsPcHIPSwQuzrFtBmCD81y+lTyE4f3iA1GaiIhIt+hQD8vChQuxLIu777672TbPPvss5513HvHx8cTHx5Odnc1nn33m1+aGG27Asiy/acaMGR0prVc5ObCca9tEGG72mSR+dMVdvHTd5VhNnDYSERHpLdodWNasWcMzzzxDZmbmKdvl5eVx7bXX8sEHH7Bq1SoyMjKYNm0a+/fv92s3Y8YMDh486Jtefvnl9pbW69R43L75y8eOY7p9AwArvGcQYtdZPRER6f3a9W1XXl7O3LlzefbZZ4mPP/WpiD//+c/88Ic/ZMKECYwaNYrnnnsOr9fLihUr/No5nU5SU1N9U0v77UvcnhODxV0RV8ScqgJCjJP8qCxOT+kXwMpERES6R7sCyx133MGsWbPIzs5u82srKyupra0lISHBb31eXh7JycmMHDmS22+/naKiomb34XK5KC0t9Zt6sw37D/jmzzr0NgBZ465k+R23EeKwB6osERGRbtPmi25zc3P5/PPPWbNmTbve8L777qN///5+YWfGjBlcdtllDB48mB07dnD//fczc+ZMVq1ahd3e+As5JyeHBQsWtOv9expjDH/6fB0AsVTA1rrAYpt8IzadDhIRkT6iTYFl79693HXXXbz33nuEhYW1+c0WLlxIbm4ueXl5fq+/5pprfPPjxo0jMzOToUOHkpeXx9SpUxvtZ/78+fz4xz/2LZeWlpKRkdHmenqCFd/s8s1/y7YF3G4YcCakjAlgVSIiIt2rTf9EX7duHYWFhUycOBGHw4HD4WDlypU88cQTOBwOPB5Ps6999NFHWbhwIe+++26LF+oOGTKEpKQktm/f3uR2p9NJTEyM39RbvbSmbrwVC8N5tk11KydeH8CKREREul+belimTp3Kpk2b/NbdeOONjBo1ivvuu6/J0zcADz/8ML/5zW9Yvnw5kydPbvF99u3bR1FREWlpaW0pr1fqFxkFwEhrL6O9RRAaDWMvC3BVIiIi3atNgSU6OpqxY8f6rYuMjCQxMdG3ft68eaSnp5OTkwPAb3/7Wx544AGWLFnCoEGDKCgoACAqKoqoqCjKy8tZsGABl19+OampqezYsYN7772XYcOGMX369M74jD1aZU0NAOfbNpLqdsMZ10FoZICrEhER6V6dftVmfn4+Bw8e9C0vWrSImpoarrjiCtLS0nzTo48+CoDdbmfjxo1ceumljBgxgptuuolJkybx0Ucf4dRw85TXuojAxXfNNyR5vDBxXqBLEhER6XYdHpo/Ly/vlMu7d+8+5evDw8NZvnx5R8votSpcNUyyfUOCxwX9RkHahECXJCIi0u10X2yQ++LgfrJsW3BgIPNq0BD8IiLSBymwBLGq2lqSrFJGWvuptmyQeVWgSxIREQkIBZYgVev28Ke1m5hifQVAYv8pEDsgwFWJiIgERoevYZGucffSd3hn65c85NgCQNiEawNckYiISOCohyUIGWN4Z+uXDLYKSLOOUYMda/TsQJclIiISMAosQWjuS68AMMX2NQDrvcMgrPeO5isiItISBZYgtHrvHiwMk2xbAVjjHRXgikRERAJLgSVIDbEOkkA5VYSy2QwKdDkiIiIBpcASpCbbvgHgC+8QanVttIiI9HEKLEHG4/EePx1UF1jWekcEuCIREZHAU2AJMv/48huGWAdJpIwqQthsBge6JBERkYBTYAki3xQe5Sev/6PB6aCh1OIgIUJPZxYRkb5NF0cEkZl/eL7R6aCZo0bz2PdmBrgyERGRwFIPS5AZZBX4nQ760XlTcNj1axIRkb5N34RBZoJtBwCbvYOpxUGkMzTAFYmIiASeAkuQGW/VBZYN3mEApEVHBbIcERGRoKBrWIKE2+MlySolwzqCB4sF3/8ZQzOGYFlWoEsTEREJOPWwBInymhpf78p205+w6CSFFRERkeMUWILEodJyxlvbgbrTQWEh6vwSERGpp8ASJIpLChlp2wfAF2YoYY6QAFckIiISPBRYgkTY3o9w4OWgSeCQiVcPi4iISAMKLEEict9KADaYoQAae0VERKQBfSsGA6+HmIJVQN1w/E9fcXmACxIREQkuCixBwBzYgM1VTCWhWAnjuGjUkECXJCIiElQUWILAvX/KYbPTydfe0/im6GigyxEREQk6CiwB9v43uxnOLgA2m0GBLUZERCRIKbAE2NN5/2SodQCALQosIiIiTVJgCRBjDL946328hz/DjuGgieeIieGcQYMDXZqIiEjQUWAJkLe+3M6Sz9cx1rYbONG7csOZEwNXlIiISJBSYAmQf1+6DIAx1m4AtngHARAd5gxMQSIiIkFMw6l2s6raWjYfPAxAinWMflYpbmxsNRkAjE9PCWR5IiIiQUmBpZtd/+fXWLcvH4DR1h4AtnnTcVH37CCnQ78SERGRk+mUUDerDysAI629AHxtBgKQHBUdkJpERESCnQJLAI201QeWutNBi66YE8BqREREgpcCS4D0t4qIoQoXDnabNEb2S2HCgNRAlyUiIhKUFFgCZJRVd2pouzcdNzb+4/xzA1yRiIhI8FJgCZBRlv/poLiIsECWIyIiEtQUWALAwviuX6m/ndluswJZkoiISFBTYOlG1bVuANKtIqKophoHu03ddStnpOv6FRERkeZo0I9uVO5yATDKN/7KALJHjOaeqediWephERERaU6HelgWLlyIZVncfffdp2z3yiuvMGrUKMLCwhg3bhxvvfWW33ZjDA888ABpaWmEh4eTnZ3Ntm3bOlJa0DHG8OvlKwEYae0D6k4HTRs1nMGJ8YEsTUREJOi1O7CsWbOGZ555hszMzFO2++STT7j22mu56aabWL9+PXPmzGHOnDls3rzZ1+bhhx/miSee4Omnn2b16tVERkYyffp0qqur21te0Hl901be/GoLAMNt+wH4xmQwKjkpkGWJiIj0CO0KLOXl5cydO5dnn32W+PhT9w48/vjjzJgxg3vuuYfTTz+dBx98kIkTJ/Lkk08CdT0Pjz32GD//+c+ZPXs2mZmZvPTSSxw4cIBly5a1p7ygtPvoMQBSraNEU0UNdvaYFEalKrCIiIi0pF2B5Y477mDWrFlkZ2e32HbVqlWN2k2fPp1Vq1YBsGvXLgoKCvzaxMbGMmXKFF+bk7lcLkpLS/2mYFfr8QIwzDpAKIZh0SP45O4fBbgqERGRnqHNF93m5uby+eefs2bNmla1LygoICXF/wnEKSkpFBQU+LbXr2uuzclycnJYsGBBW0sPqI937QZgqLUfNzBkxPkQFRHQmkRERHqKNvWw7N27l7vuuos///nPhIUFbqCz+fPnU1JS4pv27t0bsFpaa3PBAQCGWwfwYkHG2QGuSEREpOdoUw/LunXrKCwsZOLEib51Ho+HDz/8kCeffBKXy4Xdbvd7TWpqKocOHfJbd+jQIVJTU33b69elpaX5tZkwYUKTdTidTpxOZ1tKDwpRVJFmHa1byDgrsMWIiIj0IG3qYZk6dSqbNm1iw4YNvmny5MnMnTuXDRs2NAorAFlZWaxYscJv3XvvvUdWVhYAgwcPJjU11a9NaWkpq1ev9rXpLYZadb0sB00CRCQEuBoREZGeo009LNHR0YwdO9ZvXWRkJImJib718+bNIz09nZycHADuuusuzj//fP77v/+bWbNmkZuby9q1a/nDH/4A4BvH5aGHHmL48OEMHjyYX/ziF/Tv3585c+Z0wkcMHvW3M2fpdJCIiEibdPpIt/n5+dhsJzpuzjnnHJYsWcLPf/5z7r//foYPH86yZcv8gs+9995LRUUFt956K8XFxZx77rm88847Ab1OpjPVD8k/xDrAOJeLxNEXBrgiERGRnsUyxphAF9FRpaWlxMbGUlJSQkxMTKDLaaSs2sXkRx/jdyFPMMVVQewP10DS8ECXJSIiElBt+f7Www+7Qa3Xw0DrEKF4iHLGQ+KwQJckIiLSoyiwdAO3xzDYqhtTxuo/CfSgQxERkTZRYOkGtR4vg6yDAFgDJrbQWkRERE6mwNINPF4vg61DWBis9EmBLkdERKTHUWDpBu6qElKso1gA/dXDIiIi0lYKLN3AOvQFNuCoiYWofoEuR0REpMdRYOkGjkMbAci30lpoKSIiIk1RYOkGtkMbANhlFFhERETaQ4GlG6w69AUA33h0OkhERKQ9FFi6WvlhkijFC+wxKYGuRkREpEdSYOlie79aCUCBSaAKZ4CrERER6ZkUWLrYli/zANhtUvlZdnZgixEREemhFFi6WEHhJgD2mGR+cPYZAa5GRESkZ1Jg6WJR3roh+fN1/YqIiEi7OQJdQG9VVVvLx5s3kkQpAPuM7hASERFpL/WwdJEH3nyfx95aDMBhE4MtNCawBYmIiPRgCixd5LXNGxloHQZgr0nhjPSMAFckIiLScymwdJEwRwgZViEA+aYfqdFRAa5IRESk51Jg6SJe4yXD18OSzI++fXaAKxIREem5FFi6iPHUkGYVAbCXZNLjogNckYiISM+lwNJFxkVU48BLBU6S+w0PdDkiIiI9mgJLFzkrogyAfG8yl2WODXA1IiIiPZsCSxdJqskH6k4H2awAFyMiItLDKbB0kX41+wDY6+2HzdJhFhER6Qh9k3YFY0hyHwBgH/2YNVrXsIiIiHSEAktXKD9EmKnEC3z//MtJjIoIdEUiIiI9mgJLVyj8CgMUmjgiIjRgnIiISEcpsHSFw1/jxeKAScJh0yEWERHpKH2bdoXCryiz2dhPIiVVrkBXIyIi0uMpsHQBb+FXABzwJhHlDAlwNSIiIj2fAktnM4aSwq8BOEASl2WeHuCCREREej4Fls5WeoCdlgsPFodMPCEOe6ArEhER6fEUWDrb4a+wgEITTy0KKyIiIp1BgaWzFX5Nqc3GfpMY6EpERER6DQWWTuY59CVQd/2KiIiIdA4Flk62Ztc6APZ7E7lx8lkBrkZERKR3UGDpTMZwrLbuGUIHSKJfdGSACxIREekdFFg6U9lBwqnFg0WhiSc2PCzQFYmIiPQKCiydqWg7AEdMLG5snHVaeoALEhER6R0UWDrTkW0AFJgEAKJCQwNZjYiISK/RpsCyaNEiMjMziYmJISYmhqysLN5+++1m219wwQVYltVomjVrlq/NDTfc0Gj7jBkz2v+JAul4D8sh4gGIciqwiIiIdAZHWxoPGDCAhQsXMnz4cIwxvPjii8yePZv169czZsyYRu1fe+01ampqfMtFRUWMHz+eK6+80q/djBkzeOGFF3zLTqezrZ8jOBwPLPU9LBGheo6QiIhIZ2hTYLnkkkv8ln/zm9+waNEiPv300yYDS0JCgt9ybm4uERERjQKL0+kkNTW1LaUEp+OnhA6Z+AAXIiIi0ru0+xoWj8dDbm4uFRUVZGVlteo1zz//PNdccw2Rkf63++bl5ZGcnMzIkSO5/fbbKSoqOuV+XC4XpaWlflPAuV1UlOQDdT0si664LMAFiYiI9B5t6mEB2LRpE1lZWVRXVxMVFcXSpUsZPXp0i6/77LPP2Lx5M88//7zf+hkzZnDZZZcxePBgduzYwf3338/MmTNZtWoVdnvTz+LJyclhwYIFbS29S5Uf+obPnSFUEUIxkUwemBbokkRERHoNyxhj2vKCmpoa8vPzKSkp4dVXX+W5555j5cqVLYaWf/u3f2PVqlVs3LjxlO127tzJ0KFD+ec//8nUqVObbONyuXC5XL7l0tJSMjIyKCkpISYmpi0fp9Ps/OTP5H/0E3abFB50X8eW+/6DsJA250EREZE+o7S0lNjY2FZ9f7f5lFBoaCjDhg1j0qRJ5OTkMH78eB5//PFTvqaiooLc3FxuuummFvc/ZMgQkpKS2L59e7NtnE6n706l+inQak66fkVhRUREpPN0eBwWr9fr19vRlFdeeQWXy8V1113X4v727dtHUVERaWk965RK6WH/MVhERESk87SpG2D+/PnMnDmTgQMHUlZWxpIlS8jLy2P58uUAzJs3j/T0dHJycvxe9/zzzzNnzhwSExP91peXl7NgwQIuv/xyUlNT2bFjB/feey/Dhg1j+vTpHfxo3Wv/0W3EAwW6Q0hERKTTtSmwFBYWMm/ePA4ePEhsbCyZmZksX76ciy66CID8/HxsNv9Om61bt/Lxxx/z7rvvNtqf3W5n48aNvPjiixQXF9O/f3+mTZvGgw8+2OPGYjndc5gCCwpQD4uIiEhna1NgOfkOn5Pl5eU1Wjdy5Eiau643PDzc1zvTo1UVE+EuhZAQjcEiIiLSBfQsoc5wbBdeoIQIXISQHBUd6IpERER6Fd3K0hmO7cZYcNjEMSg+kb/ecE2gKxIREelV1MPSGY7uwgCHTSwT+qeTGBkR6IpERER6FQWWznBsd11gIY4Quw6piIhIZ9O3a2c4thsvFoe9sTgUWERERDqdvl07w/GLbg8TR4it6ecfiYiISPspsHSUpxZK9mEsKDRxhDoUWERERDqbAktHFeeD8VJlhVBCJNFhPWvAOxERkZ5AgaWjju0G4Ig9CYBYBRYREZFOp3FYOshTtItdISFs8dQNFhfdwx4pICIi0hMosHTQmm/WUu1wcNgbB0BMWFhgCxIREemFdEqog6pL84G6UW4BYnRKSEREpNMpsHRQaHUBUDfKLegaFhERka6gwNIRxlDsPgLUjcECEBuuwCIiItLZFFg6orqECGoAKDIxACREhAeyIhERkV5JgaUjSvbiwFBGODU4mHX6aA3NLyIi0gV0l1BHlOzDYQxF1PWu5Hz3ogAXJCIi0jupO6AjSvbhxeKoieHl7/8/Ip2hga5IRESkV1Jg6YjifLwWFBFNeIg6q0RERLqKAktHlOzDCxw1MYQ5QgJdjYiISK+lwNIRJXsx1N0hFBaipzSLiIh0FQWWDqgq3ovBoogYnRISERHpQvqWbaedhYXs9pRgA46aaJ0SEhER6ULqYWmnd9d+jA2owU4pEYSph0VERKTLKLC0k6Oy7hlCRSYGu2XTgHEiIiJdSN+y7RTpKgTqTgd5jDfA1YiIiPRuCiztFFpdF1jqR7kVERGRrqPA0k6p3mNA3RgsIiIi0rUUWNrJVB4C4CjRAa5ERESk91Ngaaf8qrqLbtXDIiIi0vUUWNopzqoAoJgoLhw2MsDViIiI9G4KLO1RU0EELgCKTSS/mvGdABckIiLSuymwtEPtsX0AVBFCjRVOepyuYxEREelKCiztsGfvVgCKTRTj+6cHuBoREZHeT4GlHdxlB4G6wPLM1ZcGuBoREZHeT4GlHVwl+wEwjngSIiICXI2IiEjvp8DSDqa87pZmlz0usIWIiIj0EQos7bD/aD4AR4gKcCUiIiJ9gwJLOxw8PsrttipHgCsRERHpGxRY2qjW7WkwaFxkgKsRERHpG9oUWBYtWkRmZiYxMTHExMSQlZXF22+/3Wz7xYsXY1mW3xQWFubXxhjDAw88QFpaGuHh4WRnZ7Nt27b2fZouZozh3MefJo5yAOadMz3AFYmIiPQNbQosAwYMYOHChaxbt461a9dy4YUXMnv2bLZs2dLsa2JiYjh48KBv2rNnj9/2hx9+mCeeeIKnn36a1atXExkZyfTp06murm7fJ+pCbo+XmqojOPACkJw8KLAFiYiI9BFtugjjkksu8Vv+zW9+w6JFi/j0008ZM2ZMk6+xLIvU1NQmtxljeOyxx/j5z3/O7NmzAXjppZdISUlh2bJlXHPNNW0pr8u53B7irLrelRLCGRQXF9iCRERE+oh2X8Pi8XjIzc2loqKCrKysZtuVl5dz2mmnkZGR0ag3ZteuXRQUFJCdne1bFxsby5QpU1i1alWz+3S5XJSWlvpN3cHlcftOBxWbaFJjNCS/iIhId2hzYNm0aRNRUVE4nU5uu+02li5dyujRo5tsO3LkSP74xz/y+uuv86c//Qmv18s555zDvn11z+IpKKgbzyQlJcXvdSkpKb5tTcnJySE2NtY3ZWRktPVjtEt17YkelmITRVJEeLe8r4iISF/X5sAycuRINmzYwOrVq7n99tu5/vrr+fLLL5tsm5WVxbx585gwYQLnn38+r732Gv369eOZZ57pUNHz58+npKTEN+3du7dD+2utyppaYqgEoIQI7HbdZCUiItId2jyQSGhoKMOGDQNg0qRJrFmzhscff7xVISQkJIQzzjiD7du3A/iubTl06BBpaWm+docOHWLChAnN7sfpdOJ0OttaeoeV1biIOX5Lc2hoQre/v4iISF/V4S4Cr9eLy+VqVVuPx8OmTZt84WTw4MGkpqayYsUKX5vS0lJWr159yutiAqW8uq6HxcLwH2efH+hyRERE+ow29bDMnz+fmTNnMnDgQMrKyliyZAl5eXksX74cgHnz5pGenk5OTg4Av/71rzn77LMZNmwYxcXFPPLII+zZs4ebb74ZqLuD6O677+ahhx5i+PDhDB48mF/84hf079+fOXPmdO4n7QQVNTXEWBVEGIMztuk7n0RERKTztSmwFBYWMm/ePA4ePEhsbCyZmZksX76ciy66CID8/HxsthOdNseOHeOWW26hoKCA+Ph4Jk2axCeffOJ3ke69995LRUUFt956K8XFxZx77rm88847jQaYCwZvbPmKc6mkwrJBZHKgyxEREekzLGOMCXQRHVVaWkpsbCwlJSXExMR02fsMfegRfhfyJBG4uOAHK6HfyC57LxERkd6uLd/fus2lDa4bP44IXER7vRClHhYREZHuosDSBtHeMgDCjQ3C4gJbjIiISB+iwNIGITXHAHDZY8CyAlyNiIhI36HA0gahtcUAVDliA1uIiIhIH6PA0koej5dPd28GwKXAIiIi0q0UWFpp08FCYo8Py+8KiQtsMSIiIn2MAksrlblqfMPyVyuwiIiIdCsFllaqrKn19bBUh8YFthgREZE+RoGlld7fttPXw1KjHhYREZFupcDSSq9u3ED08R6WSnvXjaYrIiIijSmwtEGkVQXAeadnBrgSERGRvkWBpZUyU1OJoppBtbWMGaRnCImIiHQnBZbWqi3HjiHKeCEiIdDViIiI9CkKLK105NgBALz2CHA4A1yNiIhI36LA0krh5vgYLKEa5VZERKS7KbC0Qq3bQxR1F9x6FVhERES6nQJLK1S73UQeDyz9Y1MCXI2IiEjfo8DSCi63m6jjtzTbdcGtiIhIt1NgaYWqWjfRVhUWBntUv0CXIyIi0ucosLRCVa2HSKqwACISA12OiIhIn6PA0goHSkqJphovFoTHB7ocERGRPkeBpRVe/nwjkVbdc4TUwyIiItL9FFhaYdfRo77bmjXKrYiISPdTYGmFHUWHfXcJqYdFRESk+ymwtMK04SOJopo0jxvC1cMiIiLS3RRYWsFpqrBjCPManRISEREJAAWWVrDVlAFgbKEQEh7gakRERPoeBZZWsGpLAXA7ogJciYiISN+kwNIKdnfdk5o9jsgAVyIiItI3KbC0gsNTDoBxRAe4EhERkb5JgaUVbMd7WEyITgmJiIgEggJLK4R46gILTvWwiIiIBIICSyu4a+vuEorQLc0iIiIBocDSgqraWszxHpb42H4BrkZERKRvUmBpwS25rxNBNQDR0UkBrkZERKRvUmBpwao9u4jABYAVHhvgakRERPomBZZWCLfqAgthCiwiIiKBoMByCsYYAF8PiwKLiIhIYCiwnEJpdV1QiVAPi4iISEApsJzC0coqACKoJsp4ISwusAWJiIj0UW0KLIsWLSIzM5OYmBhiYmLIysri7bffbrb9s88+y3nnnUd8fDzx8fFkZ2fz2Wef+bW54YYbsCzLb5oxY0b7Pk0nO1pZjR0v4dQyzuVSYBEREQmQNgWWAQMGsHDhQtatW8fatWu58MILmT17Nlu2bGmyfV5eHtdeey0ffPABq1atIiMjg2nTprF//36/djNmzODgwYO+6eWXX27/J+pEo1OTWDZ3DpOrq3EaICwm0CWJiIj0SZapv7K0nRISEnjkkUe46aabWmzr8XiIj4/nySefZN68eUBdD0txcTHLli1rdw2lpaXExsZSUlJCTEwnh4rifHhsHNid8IvCzt23iIhIH9aW7+92X8Pi8XjIzc2loqKCrKysVr2msrKS2tpaEhL8h7jPy8sjOTmZkSNHcvvtt1NUVHTK/bhcLkpLS/2mLuOqe1IzTj34UEREJFAcbX3Bpk2byMrKorq6mqioKJYuXcro0aNb9dr77ruP/v37k52d7Vs3Y8YMLrvsMgYPHsyOHTu4//77mTlzJqtWrcJutze5n5ycHBYsWNDW0tun5viDD0Mju+f9REREpJE2nxKqqakhPz+fkpISXn31VZ577jlWrlzZYmhZuHAhDz/8MHl5eWRmZjbbbufOnQwdOpR//vOfTJ06tck2LpcLl8vlWy4tLSUjI6NrTgnteB/+73uQPAZ++Enn7ltERKQP69JTQqGhoQwbNoxJkyaRk5PD+PHjefzxx0/5mkcffZSFCxfy7rvvnjKsAAwZMoSkpCS2b9/ebBun0+m7U6l+6jLqYREREQm4Np8SOpnX6/Xr7TjZww8/zG9+8xuWL1/O5MmTW9zfvn37KCoqIi0traOldQ5dwyIiIhJwbQos8+fPZ+bMmQwcOJCysjKWLFlCXl4ey5cvB2DevHmkp6eTk5MDwG9/+1seeOABlixZwqBBgygoKAAgKiqKqKgoysvLWbBgAZdffjmpqans2LGDe++9l2HDhjF9+vRO/qjtVHM8sKiHRUREJGDaFFgKCwuZN28eBw8eJDY2lszMTJYvX85FF10EQH5+PjbbibNMixYtoqamhiuuuMJvP7/85S/51a9+hd1uZ+PGjbz44osUFxfTv39/pk2bxoMPPojT6eyEj9cJfIElOrB1iIiI9GFtCizPP//8Kbfn5eX5Le/evfuU7cPDw329M0FL17CIiIgEnJ4l1JL6wKJrWERERAJGgaUlrrK6n+phERERCRgFlpb4TgnpGhYREZFAUWBpie4SEhERCTgFlpboGhYREZGAU2BpSX1gCVEPi4iISKAosLTEfXwUX0eQjAsjIiLSBymwtMRdXffTERbYOkRERPowBZaWeGrqfjpCA1uHiIhIH6bA0hL1sIiIiAScAktL3Md7WOzqYREREQkUBZaWqIdFREQk4BRYTsXrBW9t3bwCi4iISMAosJyKx3ViXhfdioiIBIwCy6nUnw4C9bCIiIgEkALLqdRfcIsFNkdASxEREenLFFhOpeEFt5YV2FpERET6MAWWU9GgcSIiIkFB5zlOJSwOvn2vTgeJiIgEmL6JTyWqH1z4s0BXISIi0ufplJCIiIgEPQUWERERCXoKLCIiIhL0FFhEREQk6CmwiIiISNBTYBEREZGgp8AiIiIiQU+BRURERIKeAouIiIgEPQUWERERCXoKLCIiIhL0FFhEREQk6CmwiIiISNDrFU9rNsYAUFpaGuBKREREpLXqv7frv8dPpVcElrKyMgAyMjICXImIiIi0VVlZGbGxsadsY5nWxJog5/V6OXDgANHR0ViW1an7Li0tJSMjg7179xITE9Op+5YTdJy7h45z99Bx7j461t2jq46zMYaysjL69++PzXbqq1R6RQ+LzWZjwIABXfoeMTEx+mPoBjrO3UPHuXvoOHcfHevu0RXHuaWelXq66FZERESCngKLiIiIBD0FlhY4nU5++ctf4nQ6A11Kr6bj3D10nLuHjnP30bHuHsFwnHvFRbciIiLSu6mHRURERIKeAouIiIgEPQUWERERCXoKLCIiIhL0FFhO4amnnmLQoEGEhYUxZcoUPvvss0CXFNQ+/PBDLrnkEvr3749lWSxbtsxvuzGGBx54gLS0NMLDw8nOzmbbtm1+bY4ePcrcuXOJiYkhLi6Om266ifLycr82Gzdu5LzzziMsLIyMjAwefvjhrv5oQSUnJ4czzzyT6OhokpOTmTNnDlu3bvVrU11dzR133EFiYiJRUVFcfvnlHDp0yK9Nfn4+s2bNIiIiguTkZO655x7cbrdfm7y8PCZOnIjT6WTYsGEsXry4qz9e0Fi0aBGZmZm+gbKysrJ4++23fdt1jLvGwoULsSyLu+++27dOx7rjfvWrX2FZlt80atQo3/YecYyNNCk3N9eEhoaaP/7xj2bLli3mlltuMXFxcebQoUOBLi1ovfXWW+ZnP/uZee211wxgli5d6rd94cKFJjY21ixbtsx88cUX5tJLLzWDBw82VVVVvjYzZsww48ePN59++qn56KOPzLBhw8y1117r215SUmJSUlLM3LlzzebNm83LL79swsPDzTPPPNNdHzPgpk+fbl544QWzefNms2HDBnPxxRebgQMHmvLycl+b2267zWRkZJgVK1aYtWvXmrPPPtucc845vu1ut9uMHTvWZGdnm/Xr15u33nrLJCUlmfnz5/va7Ny500RERJgf//jH5ssvvzS/+93vjN1uN++88063ft5A+fvf/27efPNN880335itW7ea+++/34SEhJjNmzcbY3SMu8Jnn31mBg0aZDIzM81dd93lW69j3XG//OUvzZgxY8zBgwd90+HDh33be8IxVmBpxllnnWXuuOMO37LH4zH9+/c3OTk5Aayq5zg5sHi9XpOammoeeeQR37ri4mLjdDrNyy+/bIwx5ssvvzSAWbNmja/N22+/bSzLMvv37zfGGPP73//exMfHG5fL5Wtz3333mZEjR3bxJwpehYWFBjArV640xtQd15CQEPPKK6/42nz11VcGMKtWrTLG1IVLm81mCgoKfG0WLVpkYmJifMf23nvvNWPGjPF7r6uvvtpMnz69qz9S0IqPjzfPPfecjnEXKCsrM8OHDzfvvfeeOf/8832BRce6c/zyl78048ePb3JbTznGOiXUhJqaGtatW0d2drZvnc1mIzs7m1WrVgWwsp5r165dFBQU+B3T2NhYpkyZ4jumq1atIi4ujsmTJ/vaZGdnY7PZWL16ta/Nt7/9bUJDQ31tpk+fztatWzl27Fg3fZrgUlJSAkBCQgIA69ato7a21u9Yjxo1ioEDB/od63HjxpGSkuJrM336dEpLS9myZYuvTcN91Lfpi38DHo+H3NxcKioqyMrK0jHuAnfccQezZs1qdDx0rDvPtm3b6N+/P0OGDGHu3Lnk5+cDPecYK7A04ciRI3g8Hr9fDEBKSgoFBQUBqqpnqz9upzqmBQUFJCcn+213OBwkJCT4tWlqHw3foy/xer3cfffdfOtb32Ls2LFA3XEIDQ0lLi7Or+3Jx7ql49hcm9LSUqqqqrri4wSdTZs2ERUVhdPp5LbbbmPp0qWMHj1ax7iT5ebm8vnnn5OTk9Nom45155gyZQqLFy/mnXfeYdGiRezatYvzzjuPsrKyHnOMe8XTmkX6qjvuuIPNmzfz8ccfB7qUXmnkyJFs2LCBkpISXn31Va6//npWrlwZ6LJ6lb1793LXXXfx3nvvERYWFuhyeq2ZM2f65jMzM5kyZQqnnXYaf/3rXwkPDw9gZa2nHpYmJCUlYbfbG10hfejQIVJTUwNUVc9Wf9xOdUxTU1MpLCz02+52uzl69Khfm6b20fA9+oo777yTN954gw8++IABAwb41qemplJTU0NxcbFf+5OPdUvHsbk2MTExPeZ/cB0VGhrKsGHDmDRpEjk5OYwfP57HH39cx7gTrVu3jsLCQiZOnIjD4cDhcLBy5UqeeOIJHA4HKSkpOtZdIC4ujhEjRrB9+/Ye89+zAksTQkNDmTRpEitWrPCt83q9rFixgqysrABW1nMNHjyY1NRUv2NaWlrK6tWrfcc0KyuL4uJi1q1b52vz/vvv4/V6mTJliq/Nhx9+SG1tra/Ne++9x8iRI4mPj++mTxNYxhjuvPNOli5dyvvvv8/gwYP9tk+aNImQkBC/Y71161by8/P9jvWmTZv8AuJ7771HTEwMo0eP9rVpuI/6Nn35b8Dr9eJyuXSMO9HUqVPZtGkTGzZs8E2TJ09m7ty5vnkd685XXl7Ojh07SEtL6zn/PXfKpbu9UG5urnE6nWbx4sXmyy+/NLfeequJi4vzu0Ja/JWVlZn169eb9evXG8D8z//8j1m/fr3Zs2ePMabutua4uDjz+uuvm40bN5rZs2c3eVvzGWecYVavXm0+/vhjM3z4cL/bmouLi01KSor5/ve/bzZv3mxyc3NNREREn7qt+fbbbzexsbEmLy/P7xbFyspKX5vbbrvNDBw40Lz//vtm7dq1Jisry2RlZfm219+iOG3aNLNhwwbzzjvvmH79+jV5i+I999xjvvrqK/PUU0/1qdtAf/rTn5qVK1eaXbt2mY0bN5qf/vSnxrIs8+677xpjdIy7UsO7hIzRse4MP/nJT0xeXp7ZtWuX+de//mWys7NNUlKSKSwsNMb0jGOswHIKv/vd78zAgQNNaGioOeuss8ynn34a6JKC2gcffGCARtP1119vjKm7tfkXv/iFSUlJMU6n00ydOtVs3brVbx9FRUXm2muvNVFRUSYmJsbceOONpqyszK/NF198Yc4991zjdDpNenq6WbhwYXd9xKDQ1DEGzAsvvOBrU1VVZX74wx+a+Ph4ExERYb73ve+ZgwcP+u1n9+7dZubMmSY8PNwkJSWZn/zkJ6a2ttavzQcffGAmTJhgQkNDzZAhQ/zeo7f7wQ9+YE477TQTGhpq+vXrZ6ZOneoLK8boGHelkwOLjnXHXX311SYtLc2Ehoaa9PR0c/XVV5vt27f7tveEY2wZY0zn9NWIiIiIdA1dwyIiIiJBT4FFREREgp4Ci4iIiAQ9BRYREREJegosIiIiEvQUWERERCToKbCIiIhI0FNgERERkaCnwCIiIiJBT4FFREREgp4Ci4iIiAQ9BRYREREJev8fK/TwsQNbYTIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "print(fittedParameters)\n",
    "#absError = modelPredictions - yData\n",
    "\n",
    "#SE = numpy.square(absError) # squared errors\n",
    "#MSE = numpy.mean(SE) # mean squared errors\n",
    "#RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE\n",
    "#Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))\n",
    "\n",
    "t = np.arange(0,5001,1)\n",
    "modelPredictions = func(t, *fittedParameters) \n",
    "plt.plot(dall.simulation_year, )\n",
    "plt.plot(t,modelPredictions)\n",
    "\n",
    "plt.plot(t, ds_reg_models_extend.volume_m3.sel(model_author='OGGM_v16').sel(rgi_reg='13').isel(gcm=0).isel(period_scenario=0), alpha = 0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "id": "087b6c00",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fce54bec940>]"
      ]
     },
     "execution_count": 259,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCCUlEQVR4nO3deXxU9b3/8fdkmyRkZUmAEAw7skQFFAPiRnDjcqWL5WKuAcV61fgr1KUYN6SKiUu95ValiAtdxLR6Rb0VpAgEqiKEJRhAQZSYiEBEIAshk2Tm/P4IGTJknWQmJ8m8no/HPHrmzPfM+eQ8rPP2e77n+7UYhmEIAADAJH5mFwAAAHwbYQQAAJiKMAIAAExFGAEAAKYijAAAAFMRRgAAgKkIIwAAwFSEEQAAYCrCCAAAMBVhBAAAmKpThZFNmzZp2rRp6tu3rywWi9599123jq+oqNDs2bM1evRoBQQEaPr06fXavPPOO5oyZYp69eqliIgIJSUlac2aNZ75AwAAQD2dKoycOnVKF1xwgV588cVWHW+32xUSEqJf/epXSk5ObrDNpk2bNGXKFK1atUrbt2/XVVddpWnTpmnnzp1tKR0AADTC0lkXyrNYLFq5cqVL74bNZtPDDz+sN998UydPntSoUaP09NNP68orr6x3/OzZs3Xy5MkW9a6MHDlSM2bM0GOPPea5PwAAAEjqZD0jzbnnnnu0efNmZWVl6fPPP9dNN92k6667Tl999VWrv9PhcKi0tFTdu3f3YKUAAKBWlwkjBQUFev311/XWW29p0qRJGjRokO6//35ddtllev3111v9vc8995zKysr0i1/8woPVAgCAWgFmF+ApeXl5stvtGjp0qMt+m82mHj16tOo7V6xYoYULF+q9995TTEyMJ8oEAADn6DJhpKysTP7+/tq+fbv8/f1dPgsLC3P7+7KysnT77bfrrbfeanSwKwAAaLsuE0Yuuugi2e12FRUVadKkSW36rjfffFO33XabsrKyNHXqVA9VCAAAGtKpwkhZWZkOHDjgfH/w4EHl5uaqe/fuGjp0qFJSUpSamqrf/e53uuiii/TDDz9o3bp1SkxMdIaKvXv3qrKyUsePH1dpaalyc3MlSRdeeKGkmlszs2bN0uLFizV+/HgdOXJEkhQSEqLIyMh2/XsBAPAFnerR3uzsbF111VX19s+aNUvLly9XVVWVnnzySf35z3/WoUOH1LNnT1166aVauHChRo8eLUlKSEjQt99+W+87ai/DlVdeqY0bNzZ6DgAA4FmdKowAAICup8s82gsAADonwggAADBVpxjA6nA49P333ys8PFwWi8XscgAAQAsYhqHS0lL17dtXfn6N9390ijDy/fffKz4+3uwyAABAKxQWFqpfv36Nft4pwkh4eLikmj8mIiLC5GoAAEBLlJSUKD4+3vk73phOEUZqb81EREQQRgAA6GSaG2LBAFYAAGAqwggAADAVYQQAAJiKMAIAAEzVpjCSmZkpi8WiefPmNdnu97//vYYNG6aQkBDFx8fr17/+tSoqKtpyagAA0EW0+mmanJwcLV26VImJiU22W7FihR588EG99tprmjBhgvbv36/Zs2fLYrHo+eefb+3pAQBAF9GqnpGysjKlpKRo2bJlio6ObrLtp59+qokTJ+rmm29WQkKCrrnmGs2cOVNbt25tVcEAAKBraVUYSUtL09SpU5WcnNxs2wkTJmj79u3O8PHNN99o1apVuuGGGxo9xmazqaSkxOUFAAC6Jrdv02RlZWnHjh3KyclpUfubb75Zx44d02WXXSbDMFRdXa0777xTDz30UKPHZGRkaOHChe6WBgAAOiG3ekYKCws1d+5cvfHGGwoODm7RMdnZ2Xrqqaf00ksvaceOHXrnnXf0wQcf6Iknnmj0mPT0dBUXFztfhYWF7pQJAAA6EYthGEZLG7/77rv6yU9+In9/f+c+u90ui8UiPz8/2Ww2l88kadKkSbr00kv17LPPOvf99a9/1R133KGysrImV/GrVVJSosjISBUXFzMdPAAAnURLf7/duk0zefJk5eXluey79dZbNXz4cM2fP79eEJGk8vLyeoGjtp0bOQgAAHRRboWR8PBwjRo1ymVft27d1KNHD+f+1NRUxcXFKSMjQ5I0bdo0Pf/887rooos0fvx4HThwQI8++qimTZvWYHgB4L5vfzylO/68XfdeM1TXjuxtdjkA4BaPr9pbUFDg0hPyyCOPyGKx6JFHHtGhQ4fUq1cvTZs2TYsWLfL0qQGfVFJRpSuezZYk/ddftis/c6q5BQGAm9waM2IWxozAV/1YZtORkgqN7BvZaJvH3tutP2/+1vn+YMYNzS7XDQDtwStjRgB4xolTlYoKDWw2NIx/ap2qHTX/vfDRvVdocExYvTb/t+t7l/e2aoeCA7kFCqDzYKE8oJ1t2Feki55Yq9/+Y2+zbWuDiCQ9+u7uBtuc27VZZXe0pTwAaHeEEcALDhef1vJPDqrMVl3vs8xVX0qSXv8kv9Hjj5XZ9OyaL1325eQfb7Bttd01jry5tcDNagHAXIQRwAt+vmSzHv+/vVr0wRf1Pgvwd701892Jcp04Vemy7/63dunFDV+77KvbS1Kr8Hh5vcDz1Kov67UDgI6MMAJ4waGTpyXV9FKcGzQKfix3bm/4skiXPb1BFz2xVnaHoSq7Q58cOKaN+39o8HvHP/WR3ss9pK+Olqq4vEqTntngvT8CANoJA1gBL7vi2Q36/PFrne9L6/Rk3Lr87BpPgx5apZ+P7ae3t3/X6HcdLbFpblZuk+fjQRoAnQ1hBHDTl0dKFBLor/N6dGtR+5KKahUeL9eKrQUa3ju8ybZNBZGWiosKafN3AEB7IowALVRRZdeaPUecPRPuTC5231u7tPVgwwNQPWlgr25af9+VXj8PAHgSY0aAFnr8/T0ut0je2PKtfii11WvnaGCgaXsEEUl6+ZZx7XIeAPAkwgjQALvDkK3a7rIvK6fQ5f3DK3frxhc+rndsQ0+9tIf064c3OCkaAHR0hBGgAT996RONe/Ijna60N9nu++IKl/fHymw6Vla/t6Q9TE3sY8p5AaCtCCPAOarsDu36rlilFdXaWXBCknTw2KlG21/z3xtVVFKhw8WnNe7JjzQhc317lep07chY9YsObffzAoAnEEaAc/x929nbMbV3XE6WVzbSWtp/tEy/X/eVZiz9zGs1PXHjSOf2tAv61vt8KWNFAHRihBF0eYZhqLyy/rTsjTlWejZ4GGdWfjld1fTtmmq7QwXHy5ts0xb+fmf/r1p4znkC/JhYBEDnRhhBlzcgfZVGPLZG+46Uuuw/VmZzefLluxPluvdvufr2+NlbMsaZjyuaCSN/3+be/CB3XTnIuT1xcA+ljO+vZ3+eqImDezTY3mGcrfPIOeNU5lw2wK1zA0BHwzwj6NLqjvV4ccMBXT60l37/0X79ctJALXh/j24Y3VsvpYyVJKW+ulXfnDM25ER5pQqPl+t0pWdXwp1+YZyWZNesPfPG7Zc69980Ll6GYej4qUqdPF2lXYUnNax3uHZ8e8LZZlnqOE174WNNGtJTt05M0MTBPT1aGwC0N8IIurTK6rMh4v1d3+v9Xd9Lkha8v0eStCrviPPzc4OIJOe8IlNHe/ZJlWG9w7UsdZxiwq31PrNYLOoRZlWPMKsG9ap5VPfz74qdn4/uF6m9v71WoUH83xdA18C/zdCl1b290WS7ZuYG+SDvcKvOH+Tvp0p7w70qU0bEtvh7Jp8fo6jQQCWfX3MMQQRAV8K/0dCl1e0ZaUyZrVrXPL/Ro+cN8vfTq7PHKSokSNMamBjNXTHhwdr2cLIC/BnmBaDr4d9s6HQaG0xqdxh6c2uBc4r2KrtDthaEkV2FJ+tNXtaUN24fr72/vbbRz5PPj9X+Rddr0pBeGt0vssXf2xyCCICuin+7oVP57f/t1fBHP9QXh0vqfTY3a6fS38nTxYs+UpXdoSEPr9Yvlm5u9jtPNDGHyLkiggM0cXDPRm+TdAvy17LUsS771t93heK7n11J9+bx/Vt8PgDwBdymQYf1Y5lNFdUO3fu3XFU7DL31X0l67ZODkqTn1+7XslTXib7+8fnZcR3/+uqHFp8nY9WXLW7b3LozkSGBslhc5/0Y2CtM//rN1fr+5Gl99MVRzbg4vsXnAwBfQBhBh1RZ7dDYJz9y2Xfwx7NPu+w5VHzuIUrsF+l86uS25dtafK5DJ0+3uG21/WwYmT0hQcs/zXf5PDCg8c7GvlEhSk1KaPG5AMBXcJsGHVLx6ap6++o+8dLNWj9Ht8c8pHUnK/v1lKGynhM+vv3Re7OwAkBXRRhBh9TQI7mn6qyg2+DAVIvn4khEcIBeSa2/3sudV5wNI5Ehgdrx6BRdOayXx84LAL6IMALTlVdW6285Bc7ekFO2alU1MDfHNz+UObdr14ypa1fhSY/VVFJRreQRsRrUq5vL/qBzekK6WQPUo1v9icsAAC1HGIEpPjlwTH/ZnC9JWrGlQPP/N083vvCx1uw5otGPr9HLm76pd8y9f9/V6Pdt3N+yAau/uW6YW3W++ctLXd77N7AonVGnF4c16wDAfQxghSlSXtkiSRrQM0x/OhNK8n8s1yv/+kYOQ/rz5m+bPL729/9YmU2zXtuqPd/Xf9S3rquHx6hPZLDuumKQnvlwnyRpaGyY9h8ta7B9/+6hkqSYiGAlnx+jj74oavS7+/cIdW7/94wLm6wDAFAfPSNoV+/uPKRpfzg7I+lj7+1WgN/Zfwxz8k80dJgk6cYL+zq3a9d0Wf9FUbNBRJJem32xFv1ktCwWi8YP6C5Jev4XF7q0eeDas70mQ2PDnNu1T8BMGtLwgnR3XD5Qk4b01E1j+2laYt8G2wAAGkfPCNrVvL/lurw/XFyhPlHBzvcBfpZG5/LY8s1x5/aOgpNauvFrZaxufo6QGeNc5/X46+3j9WNZpXpHBrvsv/vKQXp2TU2vSWCd2U4vH9pLH917hfpFh6ghoUEB+suc8c3WAQBoGD0jMNXtkwbIv85TME1NKnbuoNWWBJGdj05R5s9Gu+wL9PerF0SkmtVyLx9a82TMrAkJLp8NjglTcKB/s+cDALiPnhGYKswaoK+KGh63ca6jJTa3vz+6W5Bb7V+bNU7HGug1AQB4D2EEpmpmdnW3zbykv576ySj9LadQF/WPdvv4gEZ6TQAA3kMYgUe8s+M7hQb567pRfdw67rsTnp2xdMbF8bJYLPqPS1iMDgA6izaNGcnMzJTFYtG8efOabHfy5EmlpaWpT58+slqtGjp0qFatWtWWU6MDKSqt0L1/36U7/7pDdoehd3Z8p/VfHm3RsW9sKWj0s/umDNXdVw7S/7t6cItruTA+qsVtpbMr6MZFNTw4FQDgfa3uGcnJydHSpUuVmJjYZLvKykpNmTJFMTExevvttxUXF6dvv/1WUVFRrT01OpjSimrn9sFjZc7JyQ5m3KDC46cVHOSnmHD3bn3ERlh1z9WDnSvg/mH9Ac8VXMddVwxSkL+fbps4wCvfDwBoXqvCSFlZmVJSUrRs2TI9+eSTTbZ97bXXdPz4cX366acKDAyUJCUkJLTmtOiADMNQ3ndnV9CtO0/I4eIKXbd4k8or7Tqw6HoF+Le8I+7dtInOICJJyefH6qMvanpbLh/aS5taOONqc+K7h+rxfx/pke8CALROq27TpKWlaerUqUpOTm627fvvv6+kpCSlpaUpNjZWo0aN0lNPPSW73d7oMTabTSUlJS4vdExvbi10mTsk/Z085/aEzPUqP7O43Y+nKvWPz79v8ffWnQhNkqyBZ9/XTngGAOga3O4ZycrK0o4dO5STk9Oi9t98843Wr1+vlJQUrVq1SgcOHNDdd9+tqqoqLViwoMFjMjIytHDhQndLgwme++e+FrV75N3dWru38XEko+IitPvQ2dAZ6O+6yIu1Tq9K7YJ6AICuwa2ekcLCQs2dO1dvvPGGgoNbNgbA4XAoJiZGL7/8ssaOHasZM2bo4Ycf1h//+MdGj0lPT1dxcbHzVVhY6E6ZaEf/ltiyp2eaCiKSXIKIpHq3dOqultvYtOwAgM7JrZ6R7du3q6ioSGPGjHHus9vt2rRpk1544QXZbDb5+7vOUtmnTx8FBga67D///PN15MgRVVZWKiio/qRUVqtVVitd8R2d3WEoOtS9ScVaKuCc5W/rhpErh8ZI2lPvmD/+55h6+wAAHZ9bPSOTJ09WXl6ecnNzna9x48YpJSVFubm59YKIJE2cOFEHDhyQw+Fw7tu/f7/69OnTYBBB52B3GLrmvzdq8bqvWnX8xgeu1LUjY53vZ58z/Xq9MFKnpyQ46Oz2BfFRmjSkp3YtuMbtOU4AAB2DW2EkPDxco0aNcnl169ZNPXr00KhRoyRJqampSk9Pdx5z11136fjx45o7d67279+vDz74QE899ZTS0tI8+5egXR0/VamvfzjV6uPP69FNVfaz06/OuNh1MTv/c8JIpf1smA0J9Ff69cOVNLCHsn55qf4yZ7wiQwJbXQsAwFweXyivoKBAhw8fdr6Pj4/XmjVrlJOTo8TERP3qV7/S3Llz9eCDD3r61PCCb388pduW5ygn/7jL/lO26kaOaN7/3jVBklRw/Ozsq2HWAD32byOc7+s+1itJ7+WefRInONBf/3XFIL15x6UKCWLxOgDo7No8HXx2dnaT7yUpKSlJn332WVtPBRP8KitXuwpPav2XRcrPnOrcv6PgRBNHNa32SZljZWcXvrM7DA2NDW/0mJ9cFKfln+afOZ7FpgGgK2FtGjRpV+HJBvdv+7YtYaQmTJwsP/uIbnz3UMVFh+iC+CgN7hVW75hHpp6v706c1jUjYut9BgDo3AgjaJVx50VrRZ11ZX46Jk7v7DjUomNre0Z6hgXpWFmlpJoxIv6y6L20iQ0eE+Dvp1dmjWtj1QCAjoj+brTYMx9+qQ93H5HUeI9JQx66YbjL+9rZVedcNlCSFBzIP4YA4MvoGUGjDhSVubx/KftrSVJ4cIDL4niS9J+XnlevZ+TtO5O09oujmjUhQU+t+tK5v/ZJmTmXDVBkSCCTmAGAj+M/SdGo706UN7j/3CAiSWP6R2vFL8croUeoJGlqYh+NS+iu9OvPlzXA9YmXuKgQSTUTmd08vr/iu4d6uHIAQGdCzwgaVVHV+GKGDZkwqKfe/3+Xaf0XRUpuZKDpiD4R8jtnDhEAgG8jjKBRFVWO5hudIyI4UNMvimv087jokLaUBADogrhNg0b9eKqyRe2emD6q2Tav33qxrhrWS0+2oC0AwLfQM4IGVdkdeuIfe5tt9/ItY3XNyN7NtrtqWIyuGhbjidIAAF0MPSOoZ1v+cQ15eHWL2rYkiAAA0BTCCPRjmU0f7j6s6jOL0d22PMfl83MXrQMAwJMII9Bz/9ynO/+6Q09/WDMXSNk5i+AN6tWtweOW33qx12sDAHR9hBHoza2FkqQ3zkzv7jBcP68+d8cZVzIGBADgAYQROJVX2pXw4Af19lfbGw4jAAB4AmEEzXIYhBEAgPcQRtCs2unb67p6OLdoAACeQRiBIoIbnm4mOjRQnzx4tctaNMN7h2tJyhgt/o8L26k6AEBXx6RnPq6kokolDSx8J0k7Hp0ii8WiMedFae/hEknSn267RLERwe1ZIgCgi6NnxEcVHi/XvX/PVeLj/2y0jcVSM7/I5OE1i94FBfgRRAAAHkfPSBdWUlGlJdlfK8waoDsuH6hA/5rs+d2Jck16ZkOLv+fKYb30lzmXaGhsuLdKBQD4MMJIF7ZiS4GWZH8tSbJVO3TvlKFakv21c3KzpiyYNsK5bbFYNGlIL6/VCQDwbYSRLuzVjw86t/+Y/bVKTldp+af5zR4385J4zZ6Q4L3CAACogzEjXdjpSrtzu9LuaFEQkaTZEwY4x4sAAOBthJEuLDTIv1XHWQP4xwIA0H58+zaNvVoq+FSqOm12JV5xd1y+Nu7/we3jIgrt0nGrFyoCAHRY8ZdIIdGmnNq3w8gn/y2tf9LsKrxmtqTZQa048D0PFwIA6PjmfCTFm7Mau2+HkeLvav43vE/Nqwv49ni5HA5DCT27add3J1t8XFRIkE6erpQkjY6LlD9jRgDAtwSFmnZq3w4jtcbNka54wOwq2ux0pV1XPPahJGlSbE/9q/JYi4/NfXCKrvztWgX5+2nfHddJhBEAQDshjHQhpyrPTuv+r69aHkQkKSo0SJ88eLVCA/15kgYA0K54bKKLKK+s1rs7D7l1zOePX6Oo0EAt/PeRkmpW543u1ppBJgAAtB49I12AYRga8dgat4+LCA5U7mPXeKEiAABazrfDiGGYXYHbHA5D//3Rfo05L1pXDYuRJO35vqRFx/5y0gBFhgRq13fFuiShuzfLBACgxXw7jHRC/8g7rD+sPyBJys+cKkn68VRlk8dMGtJTw3uH6+GpI5psBwCAGdo0ZiQzM1MWi0Xz5s1rUfusrCxZLBZNnz69Laf1vE40XvO7E+XO7YX/t0clFVX6cPeRJo/5y5zxBBEAQIfV6p6RnJwcLV26VImJiS1qn5+fr/vvv1+TJk1q7Skh1ztLr3+Sr9c/yW+y/Y0X9vVuQQAAtFGrekbKysqUkpKiZcuWKTq6+alj7Xa7UlJStHDhQg0cOLA1p0Qr/e6mC8wuAQCAJrUqjKSlpWnq1KlKTk5uUfvf/va3iomJ0Zw5c1rU3mazqaSkxOWFGoabg24D/Hl6GwDQsbl9myYrK0s7duxQTk5Oi9p//PHHevXVV5Wbm9vic2RkZGjhwoXultYKnetpmoIfy/XcP/ebXQYAAB7l1n82FxYWau7cuXrjjTcUHBzcbPvS0lLdcsstWrZsmXr27Nni86Snp6u4uNj5KiwsdKfMLuvyZzc02+bfErvGGjsAAN/hVs/I9u3bVVRUpDFjxjj32e12bdq0SS+88IJsNpv8/f2dn3399dfKz8/XtGnTnPscDkfNiQMCtG/fPg0aNKjeeaxWq6zW9lzCvhM9TtOEhB6hmji4p/7x+WFJ0sxL4k2uCACA5rkVRiZPnqy8vDyXfbfeequGDx+u+fPnuwQRSRo+fHi99o888ohKS0u1ePFixcfzY+lJ+T+W6xfj4lVld2hM/2id3yfC7JIAAGiWW2EkPDxco0aNctnXrVs39ejRw7k/NTVVcXFxysjIUHBwcL32UVFRklRvPzzD38+i1KQEs8sAAKDFPD4Da0FBgfz8eIKjPfXvHqqC4+WaPSHB7FIAAHBbm8NIdnZ2k+/PtXz58rae0nM64do0DVl59wRt3P+DrhvV2+xSAABwG2vTSJKl4w9g/eTAsUY/6xFm1U/H9GvHagAA8Bzup3QSKa9saXD/a7PHtXMlAAB4FmGkE5p/3XBJ0uCYMF09PNbkagAAaBtu03RCt05MUP/uobpqeC+zSwEAoM18PIx0vgGsPx0Tp+BAf01lplUAQBfBbZpOYnjvcEnSzxioCgDoYny8Z6RWx32apqi0Qv+z7it9d+K0pJpJzQAA6EoIIx3cJYvWubwP9KczCwDQtRBGOiiHw9CxU7Z6+wf3CjOhGgAAvIcw0kHd+dft+ufeo/X2R4YGmlANAADe49t9/h30YZoqu6PBIMJ4EQBAV+TbYaSD+vu2wgb32x0dND0BANAGhBGpw61NU/vkDAAAvoAxIx3I0ZIKHTx2Sv4dLBwBAOBNhJEO5MYXPtGRkopGP//7fyW1YzUAALQPbtN0IE0FkT/ddokuGdC9HasBAKB9+HjPSMcYEHr8VKW+OlraZJtAnqQBAHRRPh5GOoarnstW8emqJttYA+nEAgB0TfzCSTJ7bZrmgogkXdAvyvuFAABgAsJIJxHAmjQAgC6KX7hOICbcanYJAAB4DWGkE5g9McHsEgAA8BrfHsBqdIynaRoT5O+nv8y5ROMSeKQXANB1+XYYqdVBZzztZvXX+IE9zC4DAACv4jZNB/PlE9c5txN6djOxEgAA2gdhpAOZPSFBwYH++tsdl+qqYb30+xkXml0SAABex20ak9kdZ8etzJ6QIEkaP7AHt2cAAD7Dx3tGzB/AunTT185tZlkFAPgifv1M9syH+5zbwQH+JlYCAIA5CCOSzJoOvqLK7vI+MiTQlDoAADATYcRE350od3nvx8q8AAAfRBgxUWlFtdklAABgOsKIicpsZ8PIhvuvNK8QAABM5NthxOTp4MvO9IxcnBCtAUxwBgDwUW0KI5mZmbJYLJo3b16jbZYtW6ZJkyYpOjpa0dHRSk5O1tatW9ty2i6j9jZNNyvTvQAAfFerw0hOTo6WLl2qxMTEJttlZ2dr5syZ2rBhgzZv3qz4+Hhdc801OnToUGtP7XkmrU3zu7U1j/Vm7/vBlPMDANARtCqMlJWVKSUlRcuWLVN0dHSTbd944w3dfffduvDCCzV8+HC98sorcjgcWrduXasK7kqOltjMLgEAANO1KoykpaVp6tSpSk5OdvvY8vJyVVVVqXv37o22sdlsKikpcXl1Ne/lnu0Z+q/LB5pYCQAA5nJ7sEJWVpZ27NihnJycVp1w/vz56tu3b5NBJiMjQwsXLmzV93cGXx4p0dysXOf7wTFh5hUDAIDJ3OoZKSws1Ny5c/XGG28oODjY7ZNlZmYqKytLK1eubPL49PR0FRcXO1+FhYVun6tlzHma5tzbMwH+THYGAPBdbvWMbN++XUVFRRozZoxzn91u16ZNm/TCCy/IZrPJ37/h9VWee+45ZWZm6qOPPmp20KvVapXVanWntE5jyzc/atZrrk8TBfj59hPWAADf5lYYmTx5svLy8lz23XrrrRo+fLjmz5/faBB55plntGjRIq1Zs0bjxo1rfbVe0349E2krdtTbF8A08AAAH+ZWGAkPD9eoUaNc9nXr1k09evRw7k9NTVVcXJwyMjIkSU8//bQee+wxrVixQgkJCTpy5IgkKSwsTGFhvjdWotpR/9ZQgD89IwAA3+XxX8GCggIdPnzY+X7JkiWqrKzUz3/+c/Xp08f5eu655zx96k6hoUlfGTMCAPBlbZ76Mzs7u8n3+fn5bT1Fl2JvoGfkh1LmGwEA+C7fvj9gwto0dRfHc+5j9V4AgA/z7TBSy6Tp4GuNios09fwAAJiJMNJOjhRX6KcvfVJv//De4bpkQOOz0QIA0NURRtrJb/73c+0oOOmy76dj4vThvMvNKQgAgA6CMNJO9h2pv77O2j1HTagEAICOxcfDSPsNYA0KqH+pG9oHAICv4dewnUSFBNXb95c5402oBACAjqXN84x0Dd5/mqZX+Nm1du65arDuv3aY188JAEBnQM9IO6myO5zbLEUDAMBZhJF2Ul5pd277kUYAAHAijLSTE+WVzu1qe/vP/AoAQEfl22GkPaeDr3Oqj77gkV4AAGr5dhhpR9V1FsjrFx1qYiUAAHQshBGpXdamqTuA9cHrh3v9fAAAdBaEkXZSdWacyKpfTdLgmDCTqwEAoONgnhEvMwxDP5TZdKzMJolZVwEAOBdhxMvuWbFTH+Qddr4P9OexXgAA6vLx/0z3/tM0dYOIJAX6+/glBwDgHPwyepHRwKPDwYH+JlQCAEDHxW0aSd5Ym2bv9yX606f59fZHhwZ6/FwAAHRmhBEvueF//tXgfks7PEYMAEBnwm0aLzhSXGF2CQAAdBqEES84UFTW4P5fjOvXzpUAANDx+fZtGg+vTVNaUaV5WbnqGxXS4OdP/yzRo+cDAKAr8O0wUstD4zhe3vSN1n1Z1MRpGC8CAMC5uE3jQcdPVZpdAgAAnQ5hpJ2su+8Ks0sAAKBDIox4ULW98TEog3qxOB4AAA3x8TDi2QGs1Y6Gv29473CPngcAgK7Ex8OIZ1U7HA3u//JIaTtXAgBA50EYkeSp6eCPldk88j0AAPgSwogHfXLgxwb3L/6PC9u3EAAAOhHmGfEgP4t07rCRP8y8SNMu6GtOQQAAdAL0jHjIKVu1+ncPrbefIAIAQNPaFEYyMzNlsVg0b968Jtu99dZbGj58uIKDgzV69GitWrWqLaf1HA9OB3/94n8p/8dyj30fAAC+otVhJCcnR0uXLlViYtPrrXz66aeaOXOm5syZo507d2r69OmaPn26du/e3dpTd0gFx+sHkf+9K8mESgAA6FxaFUbKysqUkpKiZcuWKTo6usm2ixcv1nXXXacHHnhA559/vp544gmNGTNGL7zwQqsK9goPrxnTv3uobrn0PI09r7tHvxcAgK6oVWEkLS1NU6dOVXJycrNtN2/eXK/dtddeq82bNzd6jM1mU0lJicurIzPOud2T+dPRemL6KJOqAQCgc3H7aZqsrCzt2LFDOTk5LWp/5MgRxcbGuuyLjY3VkSNHGj0mIyNDCxcudLc003xx2HVSs6AAxgUDANBSbv1qFhYWau7cuXrjjTcUHBzsrZqUnp6u4uJi56uwsNBr5/KERav2urwnjAAA0HJu9Yxs375dRUVFGjNmjHOf3W7Xpk2b9MILL8hms8nf39/lmN69e+vo0aMu+44eParevXs3eh6r1Sqr1epOaa3kmadpjpVWurwnjAAA0HJu/WpOnjxZeXl5ys3Ndb7GjRunlJQU5ebm1gsikpSUlKR169a57Fu7dq2SkrrOkyaDY11X5I0KCTKpEgAAOh+3ekbCw8M1apTrwMxu3bqpR48ezv2pqamKi4tTRkaGJGnu3Lm64oor9Lvf/U5Tp05VVlaWtm3bppdfftlDf4L5ugWdDWH3TRmq3pHeu4UFAEBX4/H7CQUFBTp8+LDz/YQJE7RixQq9/PLLuuCCC/T222/r3XffrRdqOrO/b/vOuf3/Jg8xsRIAADqfNq9Nk52d3eR7Sbrpppt00003tfVUAACgC2KkpQf0iw6RJM1KOs/kSgAA6Hx8O4x4aG2a4b3DJUnn94nwyPcBAOBLfDuM1GrjdPB2R02o8fPz7LTyAAD4AsKIB1SfCSMBhBEAANxGGPGA2p4Rf8IIAABuI4x4wNmeES4nAADu8vFfT88MYHU4e0Y88nUAAPgUfj49oNoZRricAAC4i19PSVLbxnocKCqTxABWAABagzDSRp8cOKYyW7UkyVbtMLkaAAA6H8JIG93/1i7nduHxchMrAQCgcyKMtNHh4grn9k/HxJlYCQAAnZNvhxEPTQdfq0eY1aPfBwCAL/DtMOJBsREEEQAAWoMwIrV6bZqlG792bi/5z7GeqgYAAJ9CGGmDjNVfOrcjggNNrAQAgM4rwOwCOiPDMLR43Vcu++KiQkyqBgCAzo2ekVb4/Lti/f4j1zASEuRvUjUAAHRuhJFWYHIzAAA8hzAiyd3p4IMCuGwAAHgKv6qtYHe4zk/ym+uGmVQJAACdH2GkFartrrdpbr9soEmVAADQ+RFGWqH6nJ4RbtsAANB6/Iq2wrlhBAAAtJ5vh5FWrk1z7m0aAADQer4dRmq5OR18lZ2eEQAAPIUw0gp1n6ZJ7BdpYiUAAHR+hJFWqHacvU3zSuo4EysBAKDzI4y0QvWZ2zRXDO2lmIhgk6sBAKBz8/Ew0soBrGd6RgL83BtrAgAA6vPxMNI6tQNY/QgjAAC0GWFEkrtr0/x5c74kae3eo16oBQAA30IYaYX9R8vMLgEAgC7DrTCyZMkSJSYmKiIiQhEREUpKStLq1aubPOb3v/+9hg0bppCQEMXHx+vXv/61Kioq2lR0R/Hk9FFmlwAAQKcX4E7jfv36KTMzU0OGDJFhGPrTn/6kG2+8UTt37tTIkSPrtV+xYoUefPBBvfbaa5owYYL279+v2bNny2Kx6Pnnn/fYH9GeSiuqnNv+jBkBAKDN3Aoj06ZNc3m/aNEiLVmyRJ999lmDYeTTTz/VxIkTdfPNN0uSEhISNHPmTG3ZsqUNJXtQK6aDL6modm5fOayXJ6sBAMAntXrMiN1uV1ZWlk6dOqWkpKQG20yYMEHbt2/X1q1bJUnffPONVq1apRtuuKHJ77bZbCopKXF5dRSr8w47t/tEhphYCQAAXYNbPSOSlJeXp6SkJFVUVCgsLEwrV67UiBEjGmx7880369ixY7rssstkGIaqq6t155136qGHHmryHBkZGVq4cKG7pbWeG2vTLP8033t1AADgg9zuGRk2bJhyc3O1ZcsW3XXXXZo1a5b27t3bYNvs7Gw99dRTeumll7Rjxw698847+uCDD/TEE080eY709HQVFxc7X4WFhe6W6RWGYei7E6fNLgMAgC7F7Z6RoKAgDR48WJI0duxY5eTkaPHixVq6dGm9to8++qhuueUW3X777ZKk0aNH69SpU7rjjjv08MMPy8+v4SxktVpltVrdLc3rcvJPOLeHxISZWAkAAF1Hm+cZcTgcstlsDX5WXl5eL3D4+/tLqull6GwqquzO7a+KmGsEAABPcKtnJD09Xddff7369++v0tJSrVixQtnZ2VqzZo0kKTU1VXFxccrIyJBU8/TN888/r4suukjjx4/XgQMH9Oijj2ratGnOUGIu9wJR3bVowqxudyoBAIAGuPWLWlRUpNTUVB0+fFiRkZFKTEzUmjVrNGXKFElSQUGBS0/II488IovFokceeUSHDh1Sr169NG3aNC1atMizf0WbtWwAa0X12Z6Rf0vs461iAADwKW6FkVdffbXJz7Ozs12/PCBACxYs0IIFC9wurCMqrTPHSIA/E54BAOAJrE3jhsPFZ6exD2hk8C0AAHAPv6huOFl+dir4iGDGjAAA4AmEETfUfZpmZFykiZUAANB1+HYYcfPxYtuZAayhQf66ZkSsNyoCAMDn+HYYqdWC6eD/llOgN7fWzAT76+ShsrgxhTwAAGgcYaSF5v9vnnP7dJ3bNQAAoG0II62wqs7KvQAAoG0II61wfp8Is0sAAKDL8PEw0rIBrA6Ha7v/uDjeG8UAAOCTfDyMtExZZbXLe38/Bq8CAOAphBFJza1NY7e79ox0vvWGAQDouAgjLVDlcLi8d3N6EgAA0ATCSAtUn9MzMqBnN5MqAQCg6yGMtIC9zgDWj+69Qr3CrSZWAwBA1+LbYaSF91uq7DW3acKDAzQ4JsybFQEA4HN8O4y0UPWZnpFAfy4XAACexq+r1OzaNLU9IwE80gsAgMcRRlqgdgArYQQAAM8jjLRA9ZlHewO4TQMAgMfx69oC5ZU1q/SGBvmbXAkAAF2Pj4eRlj1NU1ZRMx18mDXAm8UAAOCTfDyM1Gp6LEiZrSaMdCOMAADgcYSRZqzYUqAH3v5cEj0jAAB4A2GkGQ+tzHNud7MyZgQAAE8jjDTB4XAdU8JtGgAAPI8w0oTj5ZUu73cfKjapEgAAui7fDiPNrE2T951r+Cg981QNAADwHN8OI7UamQ7+z5vzXd5PvyiuHYoBAMC3EEaacPGA7i7v51w2wKRKAADouggjTah7F+eaEbGs2gsAgBfw69oEW3XNmjS/GNdPL6eOM7kaAAC6Jh8PI00PYK08E0bCgwPboxgAAHySj4eRptWGkaAALhMAAN7Cr6ykxtamqbTXrNYbxFgRAAC8xq1f2SVLligxMVERERGKiIhQUlKSVq9e3eQxJ0+eVFpamvr06SOr1aqhQ4dq1apVbSq6vdAzAgCA97k1v3m/fv2UmZmpIUOGyDAM/elPf9KNN96onTt3auTIkfXaV1ZWasqUKYqJidHbb7+tuLg4ffvtt4qKivJU/V5jGIZOlFdJkqyEEQAAvMatMDJt2jSX94sWLdKSJUv02WefNRhGXnvtNR0/flyffvqpAgNrBoEmJCS0vtp2NOdP27T+yyJJ9IwAAOBNrf6VtdvtysrK0qlTp5SUlNRgm/fff19JSUlKS0tTbGysRo0apaeeekr2M2MxGmOz2VRSUuLy8oompoOvDSISY0YAAPAmt5ehzcvLU1JSkioqKhQWFqaVK1dqxIgRDbb95ptvtH79eqWkpGjVqlU6cOCA7r77blVVVWnBggWNniMjI0MLFy50t7TWa2Q6+Fr0jAAA4D1u/8oOGzZMubm52rJli+666y7NmjVLe/fubbCtw+FQTEyMXn75ZY0dO1YzZszQww8/rD/+8Y9NniM9PV3FxcXOV2FhobtlehRhBAAA73G7ZyQoKEiDBw+WJI0dO1Y5OTlavHixli5dWq9tnz59FBgYKH9/f+e+888/X0eOHFFlZaWCgoIaPIfVapXVanW3NK/hNg0AAN7T5l9Zh8Mhm83W4GcTJ07UgQMH5HA4nPv279+vPn36NBpEOgKHw3UsCT0jAAB4j1u/sunp6dq0aZPy8/OVl5en9PR0ZWdnKyUlRZKUmpqq9PR0Z/u77rpLx48f19y5c7V//3598MEHeuqpp5SWlubZv8LDquqEJwAA4F1u3aYpKipSamqqDh8+rMjISCUmJmrNmjWaMmWKJKmgoEB+fmfzTXx8vNasWaNf//rXSkxMVFxcnObOnav58+d79q9otYafpjllc33ax9LMAFcAANB6boWRV199tcnPs7Oz6+1LSkrSZ5995lZR7c81bKzefdjl/WWDe7ZnMQAA+BQGQzQgMsR1lV5/P3pGAADwFsJIA+qOX1133xXmFQIAgA8gjDTgV2/udG4P6hVmYiUAAHR9hBEAAGAq3w4jTaxNI0kTB/dop0IAAPBdvh1Gap3z6G589xBJ0r1ThplRDQAAPoUw0oDaeUbCg92eLR8AALiJMNKA46cqJUndrIQRAAC8jTByjtzCk87tMMIIAABe5+NhpP4A1jV7jji3uwX51/scAAB4lo+HEVelFVU6dOK0832AP5cHAABv4z5EHXOzcrX+yyJJUvL5MSZXAwCAb+A//euoDSKSFBzILRoAANoDYaQR1gDCCAAA7YEwUsfAXt2c25V2h4mVAADgO3w7jJwzHXz/7qHO7a+OlrZ3NQAA+CTfDiO1zkwHf7rS7tyVcul5ZlUDAIBPIYzUUVF1NoyM6R9lXiEAAPgQwkgd5Wd6Rv7rioEa2TfS5GoAAPANhJE6asPI9aP6mFwJAAC+gzBSx6GTNbOvhjDHCAAA7YYwcsaPZTbnNrPAAwDQfvjZlSRZZHecfcw3Liq0ibYAAMCTCCNn2OvMORLCar0AALQbwsgZ1faaMGIN4JIAANCe+OU9o/Y2TYCfxeRKAADwLYSRM2pv0/gTRgAAaFe+HUbqjBP5556jkqSSimqzqgEAwCf5dhipZbHo6Q+/NLsKAAB8EmEEAACYijACAABMRRgBAACm8vEwUjOA9WiprZl2AADAW3w8jNTY/V2x2SUAAOCz3AojS5YsUWJioiIiIhQREaGkpCStXr26RcdmZWXJYrFo+vTpranTq6K7BTm3/3zbJSZWAgCA73ErjPTr10+ZmZnavn27tm3bpquvvlo33nij9uzZ0+Rx+fn5uv/++zVp0qQ2FestldUOSdLlQ3vp8qG9TK4GAADf4lYYmTZtmm644QYNGTJEQ4cO1aJFixQWFqbPPvus0WPsdrtSUlK0cOFCDRw4sM0Fe0OlvSaMhARy1woAgPbW6l9fu92urKwsnTp1SklJSY22++1vf6uYmBjNmTOnxd9ts9lUUlLi8vKm0jOzrgb4EUYAAGhvAe4ekJeXp6SkJFVUVCgsLEwrV67UiBEjGmz78ccf69VXX1Vubq5b58jIyNDChQvdLc19Z6aD/8fnhyX11wd5h/Wi988KAADqcLsrYNiwYcrNzdWWLVt01113adasWdq7d2+9dqWlpbrlllu0bNky9ezZ061zpKenq7i42PkqLCx0t0y3GGJxPAAAzOJ2z0hQUJAGDx4sSRo7dqxycnK0ePFiLV261KXd119/rfz8fE2bNs25z+GoGZsREBCgffv2adCgQQ2ew2q1ymq1ultaqw3sGSoVSbdfNqDdzgkAAGq4HUbO5XA4ZLPVnzRs+PDhysvLc9n3yCOPqLS0VIsXL1Z8fHxbT+0xYcGBkqSBvcJMrgQAAN/jVhhJT0/X9ddfr/79+6u0tFQrVqxQdna21qxZI0lKTU1VXFycMjIyFBwcrFGjRrkcHxUVJUn19put2lEzdiQogAGsAAC0N7fCSFFRkVJTU3X48GFFRkYqMTFRa9as0ZQpUyRJBQUF8uuET6RUn3m0lzACAED7cyuMvPrqq01+np2d3eTny5cvd+d07aCmR6Siyi5JCvInjAAA0N749ZX07fHTkiT7mds1AACg/RBG6ojuFmh2CQAA+BzCSB2J/aLMLgEAAJ9DGKkjJNDf7BIAAPA5hJEz/P0s8vdjJlYAANqbb4cR4+yA1f7dQ00sBAAA3+XbYcTJwmO9AACYhF/gM5jwDAAAc/ALfAZhBAAAc/ALfIaVMAIAgCn4BT6DnhEAAMzh47/AZ5+mYQArAADm4BdYNZGkm9WtNQMBAICHEEbOiI0INrsEAAB8EmHkjN4RVrNLAADAJ/l0GDlZXuXcjgoNMrESAAB8l0+Hke9Plju3GTMCAIA5fDqM1DJkEWvkAQBgDp8OI0adhfLC6BkBAMAUPh1Gqh1nw8iF/aPMKwQAAB/m02HEXqdnxBrgb2IlAAD4Lp8OI6GBNQGkX3SoyZUAAOC7fHqgRMmwm/TJiQLdesVks0sBAMBn+XQYueSm+80uAQAAn+fTt2kAAID5CCMAAMBUhBEAAGAqwggAADAVYQQAAJiKMAIAAExFGAEAAKYijAAAAFMRRgAAgKncCiNLlixRYmKiIiIiFBERoaSkJK1evbrR9suWLdOkSZMUHR2t6OhoJScna+vWrW0uGgAAdB1uhZF+/fopMzNT27dv17Zt23T11Vfrxhtv1J49expsn52drZkzZ2rDhg3avHmz4uPjdc011+jQoUMeKR4AAHR+FsMwjLZ8Qffu3fXss89qzpw5zba12+2Kjo7WCy+8oNTU1Bafo6SkRJGRkSouLlZERERbygUAAO2kpb/frV4oz26366233tKpU6eUlJTUomPKy8tVVVWl7t27N9nOZrPJZrM535eUlLS2TAAA0MG5HUby8vKUlJSkiooKhYWFaeXKlRoxYkSLjp0/f7769u2r5OTkJttlZGRo4cKF9fYTSgAA6Dxqf7ebvQljuMlmsxlfffWVsW3bNuPBBx80evbsaezZs6fZ4zIyMozo6Ghj165dzbatqKgwiouLna+9e/caknjx4sWLFy9enfBVWFjY5O9+m8eMJCcna9CgQVq6dGmjbZ577jk9+eST+uijjzRu3Di3z+FwOPT9998rPDxcFoulLeW6KCkpUXx8vAoLCxmL4kVc5/bDtW4fXOf2wXVuH968zoZhqLS0VH379pWfX+PPzLR6zEgth8PhMr7jXM8884wWLVqkNWvWtCqISJKfn5/69evX2hKbVfuoMryL69x+uNbtg+vcPrjO7cNb1zkyMrLZNm6FkfT0dF1//fXq37+/SktLtWLFCmVnZ2vNmjWSpNTUVMXFxSkjI0OS9PTTT+uxxx7TihUrlJCQoCNHjkiSwsLCFBYW5u7fAwAAuiC3wkhRUZFSU1N1+PBhRUZGKjExUWvWrNGUKVMkSQUFBS7dMEuWLFFlZaV+/vOfu3zPggUL9Pjjj7e9egAA0Om5FUZeffXVJj/Pzs52eZ+fn+9uPe3KarVqwYIFslqtZpfSpXGd2w/Xun1wndsH17l9dITr3OYBrAAAAG3BQnkAAMBUhBEAAGAqwggAADAVYQQAAJjKp8PIiy++qISEBAUHB2v8+PHaunWr2SV1WJs2bdK0adPUt29fWSwWvfvuuy6fG4ahxx57TH369FFISIiSk5P11VdfubQ5fvy4UlJSFBERoaioKM2ZM0dlZWUubT7//HNNmjRJwcHBio+P1zPPPOPtP61DycjI0MUXX6zw8HDFxMRo+vTp2rdvn0ubiooKpaWlqUePHgoLC9PPfvYzHT161KVNQUGBpk6dqtDQUMXExOiBBx5QdXW1S5vs7GyNGTNGVqtVgwcP1vLly73953UYS5YsUWJionOSp6SkJK1evdr5OdfYOzIzM2WxWDRv3jznPq61Zzz++OOyWCwur+HDhzs/7/DXuQXL0XRJWVlZRlBQkPHaa68Ze/bsMX75y18aUVFRxtGjR80urUNatWqV8fDDDxvvvPOOIclYuXKly+eZmZlGZGSk8e677xq7du0y/v3f/90YMGCAcfr0aWeb6667zrjggguMzz77zPjXv/5lDB482Jg5c6bz8+LiYiM2NtZISUkxdu/ebbz55ptGSEiIsXTp0vb6M0137bXXGq+//rqxe/duIzc317jhhhuM/v37G2VlZc42d955pxEfH2+sW7fO2LZtm3HppZcaEyZMcH5eXV1tjBo1ykhOTjZ27txprFq1yujZs6eRnp7ubPPNN98YoaGhxr333mvs3bvX+MMf/mD4+/sbH374Ybv+vWZ5//33jQ8++MDYv3+/sW/fPuOhhx4yAgMDjd27dxuGwTX2hq1btxoJCQlGYmKiMXfuXOd+rrVnLFiwwBg5cqRx+PBh5+uHH35wft7Rr7PPhpFLLrnESEtLc7632+1G3759jYyMDBOr6hzODSMOh8Po3bu38eyzzzr3nTx50rBarcabb75pGIbhXOwwJyfH2Wb16tWGxWIxDh06ZBiGYbz00ktGdHS0YbPZnG3mz59vDBs2zMt/UcdVVFRkSDI2btxoGEbNdQ0MDDTeeustZ5svvvjCkGRs3rzZMIya4Ojn52ccOXLE2WbJkiVGRESE89r+5je/MUaOHOlyrhkzZhjXXnutt/+kDis6Otp45ZVXuMZeUFpaagwZMsRYu3atccUVVzjDCNfacxYsWGBccMEFDX7WGa6zT96mqays1Pbt25WcnOzc5+fnp+TkZG3evNnEyjqngwcP6siRIy7XMzIyUuPHj3dez82bNysqKsplfaLk5GT5+flpy5YtzjaXX365goKCnG2uvfZa7du3TydOnGinv6ZjKS4uliR1795dkrR9+3ZVVVW5XOvhw4erf//+Ltd69OjRio2Ndba59tprVVJSoj179jjb1P2O2ja++M+/3W5XVlaWTp06paSkJK6xF6SlpWnq1Kn1rgfX2rO++uor9e3bVwMHDlRKSooKCgokdY7r7JNh5NixY7Lb7S4XXZJiY2Od6+eg5WqvWVPX88iRI4qJiXH5PCAgQN27d3dp09B31D2HL3E4HJo3b54mTpyoUaNGSaq5DkFBQYqKinJpe+61bu46NtampKREp0+f9saf0+Hk5eUpLCxMVqtVd955p1auXKkRI0ZwjT0sKytLO3bscK5ZVhfX2nPGjx+v5cuX68MPP9SSJUt08OBBTZo0SaWlpZ3iOrd51V4A3pGWlqbdu3fr448/NruULmnYsGHKzc1VcXGx3n77bc2aNUsbN240u6wupbCwUHPnztXatWsVHBxsdjld2vXXX+/cTkxM1Pjx43Xeeefp73//u0JCQkysrGV8smekZ8+e8vf3rzeS+OjRo+rdu7dJVXVetdesqevZu3dvFRUVuXxeXV2t48ePu7Rp6DvqnsNX3HPPPfrHP/6hDRs2qF+/fs79vXv3VmVlpU6ePOnS/txr3dx1bKxNREREp/gXlycEBQVp8ODBGjt2rDIyMnTBBRdo8eLFXGMP2r59u4qKijRmzBgFBAQoICBAGzdu1P/8z/8oICBAsbGxXGsviYqK0tChQ3XgwIFO8c+0T4aRoKAgjR07VuvWrXPuczgcWrdunZKSkkysrHMaMGCAevfu7XI9S0pKtGXLFuf1TEpK0smTJ7V9+3Znm/Xr18vhcGj8+PHONps2bVJVVZWzzdq1azVs2DBFR0e3019jLsMwdM8992jlypVav369BgwY4PL52LFjFRgY6HKt9+3bp4KCApdrnZeX5xL+1q5dq4iICI0YMcLZpu531Lbx5X/+HQ6HbDYb19iDJk+erLy8POXm5jpf48aNU0pKinOba+0dZWVl+vrrr9WnT5/O8c90m4fAdlJZWVmG1Wo1li9fbuzdu9e44447jKioKJeRxDirtLTU2Llzp7Fz505DkvH8888bO3fuNL799lvDMGoe7Y2KijLee+894/PPPzduvPHGBh/tveiii4wtW7YYH3/8sTFkyBCXR3tPnjxpxMbGGrfccouxe/duIysrywgNDfWpR3vvuusuIzIy0sjOznZ5RK+8vNzZ5s477zT69+9vrF+/3ti2bZuRlJRkJCUlOT+vfUTvmmuuMXJzc40PP/zQ6NWrV4OP6D3wwAPGF198Ybz44os+9Sjkgw8+aGzcuNE4ePCg8fnnnxsPPvigYbFYjH/+85+GYXCNvanu0zSGwbX2lPvuu8/Izs42Dh48aHzyySdGcnKy0bNnT6OoqMgwjI5/nX02jBiGYfzhD38w+vfvbwQFBRmXXHKJ8dlnn5ldUoe1YcMGQ1K916xZswzDqHm899FHHzViY2MNq9VqTJ482di3b5/Ld/z444/GzJkzjbCwMCMiIsK49dZbjdLSUpc2u3btMi677DLDarUacXFxRmZmZnv9iR1CQ9dYkvH6668725w+fdq4++67jejoaCM0NNT4yU9+Yhw+fNjle/Lz843rr7/eCAkJMXr27Gncd999RlVVlUubDRs2GBdeeKERFBRkDBw40OUcXd1tt91mnHfeeUZQUJDRq1cvY/Lkyc4gYhhcY286N4xwrT1jxowZRp8+fYygoCAjLi7OmDFjhnHgwAHn5x39OlsMwzDa3r8CAADQOj45ZgQAAHQchBEAAGAqwggAADAVYQQAAJiKMAIAAExFGAEAAKYijAAAAFMRRgAAgKkIIwAAwFSEEQAAYCrCCAAAMBVhBAAAmOr/AxWoqnPHh2yhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afe9c0ce-47d5-489c-a434-ef27772698f3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f608b8b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:oggm_gmip3_working]",
   "language": "python",
   "name": "conda-env-oggm_gmip3_working-py"
  },
  "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
}