{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Reconstruction over a smaller domain using GraphEM" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Expected time to run through: 10 mins**\n", "\n", "This tutorial demonstrates how to get a reconstruction over a smaller domain (30N-70N, 130W-70W) using GraphEM, leveraging HadCRUT4 and the PAGES2k proxy database." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test data preparation\n", "\n", "To go through this tutorial, please prepare test data following the steps:\n", "\n", "Download the test case named “PAGES2k_HadCRUT” with this [link](https://drive.google.com/drive/folders/12V0Ny4eO2HKU40-Iy5Jiw2SAfVNzXXcd?usp=sharing).\n", "Create a directory named “testcases” in the same directory where this notebook sits.\n", "Put the unzipped direcotry “PAGES2k_HadCRUT_cropped_domain” into “testcases”.\n", "\n", "Below, we first load some useful packages, including our `GraphEM`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", "\n", "import LMRt\n", "import GraphEM\n", "import os\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Low-level workflow" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "job = GraphEM.ReconJob()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mGraphEM: job.load_configs() >>> loading reconstruction configurations from: ./testcases/PAGES2k_HadCRUT_cropped_domain/configs.yml\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_configs() >>> job.configs created\u001b[0m\n", "\u001b[1m\u001b[36mGraphEM: job.load_configs() >>> job.configs[\"job_dirpath\"] = /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/recon\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_configs() >>> /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/recon created\u001b[0m\n", "{'anom_period': [1951, 1980],\n", " 'calib_period': [1930, 2000],\n", " 'job_dirpath': '/Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/recon',\n", " 'job_id': 'GraphEM_tutorial',\n", " 'obs_path': {'tas': './data/obs/HadCRUT.5.0.1.0.analysis.anomalies.ensemble_mean.nc'},\n", " 'obs_varname': {'lat': 'latitude', 'lon': 'longitude', 'tas': 'tas_mean'},\n", " 'proxydb_path': './data/proxy/pages2k_dataset.pkl',\n", " 'ptype_list': ['coral.d18O',\n", " 'coral.SrCa',\n", " 'coral.calc',\n", " 'tree.TRW',\n", " 'tree.MXD'],\n", " 'recon_period': [1500, 2000]}\n" ] } ], "source": [ "job.load_configs('./testcases/PAGES2k_HadCRUT_cropped_domain/configs.yml', verbose=True)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mGraphEM: job.load_proxydb() >>> job.configs[\"proxydb_path\"] = /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/data/proxy/pages2k_dataset.pkl\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_proxydb() >>> 692 records loaded\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_proxydb() >>> job.proxydb created\u001b[0m\n" ] } ], "source": [ "job.load_proxydb(verbose=True)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Proxy Database Overview\n", "-----------------------\n", " Source: /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/data/proxy/pages2k_dataset.pkl\n", " Size: 111\n", "Proxy types: {'tree.TRW': 70, 'tree.MXD': 41}\n" ] } ], "source": [ "# filter the proxy database by spatial region larger than the target (30N-70N, 130W-70W); using (10N-90N, 150W-60W)\n", "assim_pids = []\n", "for pid, pobj in job.proxydb.records.items():\n", " if pobj.lat >= 10 and pobj.lat <= 90 and pobj.lon >= np.mod(-150, 360) and pobj.lon <= np.mod(-60, 360):\n", " assim_pids.append(pid)\n", " \n", "job.proxydb.filter_pids(assim_pids, inplace=True)\n", "print(job.proxydb)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mGraphEM: job.seasonalize_proxydb() >>> job.configs[\"ptype_season\"] = {'tree.TRW': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], 'tree.MXD': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]}\u001b[0m\n", "\u001b[1m\u001b[36mGraphEM: job.seasonalize_proxydb() >>> seasonalizing proxy records according to: {'tree.TRW': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], 'tree.MXD': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]}\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.seasonalize_proxydb() >>> 111 records remaining\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.seasonalize_proxydb() >>> job.proxydb updated\u001b[0m\n", "\u001b[1m\u001b[36mGraphEM: job.filter_proxydb() >>> filtering proxy records according to: ['coral.d18O', 'coral.SrCa', 'coral.calc', 'tree.TRW', 'tree.MXD']\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.filter_proxydb() >>> 111 records remaining\u001b[0m\n" ] } ], "source": [ "ptype_season = {}\n", "for k, v in job.proxydb.type_dict.items():\n", " ptype_season[k] = list(range(1, 13)) # annual\n", " \n", "job.seasonalize_proxydb(ptype_season, verbose=True)\n", "job.filter_proxydb(verbose=True)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mGraphEM: job.load_obs() >>> loading instrumental observation fields from: {'tas': '/Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/data/obs/HadCRUT.5.0.1.0.analysis.anomalies.ensemble_mean.nc'}\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_obs() >>> job.obs created\u001b[0m\n" ] } ], "source": [ "job.load_obs(varname_dict={'lat': 'latitude', 'lon': 'longitude', 'tas': 'tas_mean'}, verbose=True)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[30mGraphEM: job.seasonalize_obs() >>> seasonalized obs w/ season [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]\u001b[0m\n", "Dataset Overview\n", "-----------------------\n", "\n", " Name: tas\n", " Source: /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/data/obs/HadCRUT.5.0.1.0.analysis.anomalies.ensemble_mean.nc\n", " Shape: time:171, lat:36, lon:72\n", "\u001b[1m\u001b[32mGraphEM: job.seasonalize_obs() >>> job.obs updated\u001b[0m\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/fzhu/Github/LMRt/LMRt/utils.py:261: RuntimeWarning: Mean of empty slice\n", " tmp = np.nanmean(var[inds, ...], axis=0)\n" ] } ], "source": [ "job.seasonalize_obs(verbose=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cropping the domain" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dataset Overview\n", "-----------------------\n", "\n", " Name: tas\n", " Source: /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/data/obs/HadCRUT.5.0.1.0.analysis.anomalies.ensemble_mean.nc\n", " Shape: time:171, lat:36, lon:72\n" ] } ], "source": [ "print(job.obs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The current domain is global.\n", "We may call the `job.obs.crop()` method with the `domain_range` parameter to crop the domain.\n", "\n", "There are two ways of cropping the domain:\n", "+ crop only the latitudes with the input argument `domain_range = [lat_min, lat_max]`\n", "+ crop both the latitudes and longitudes the input argument `domain_range = [lat_min, lat_max, lon_min, lon_max]`" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dataset Overview\n", "-----------------------\n", "\n", " Name: tas\n", " Source: /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/data/obs/HadCRUT.5.0.1.0.analysis.anomalies.ensemble_mean.nc\n", " Shape: time:171, lat:8, lon:12\n" ] } ], "source": [ "# crop the 30N-70N, 70W-130W region\n", "\n", "job.obs.crop([30, 70, np.mod(-130, 360), np.mod(-70, 360)], inplace=True)\n", "print(job.obs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data preparation" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.recon_time created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.calib_time created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.calib_idx created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.temp created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.df_proxy created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.proxy created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.lonlat created\u001b[0m\n" ] } ], "source": [ "job.prep_data(verbose=True)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
NAm_153NAm_165NAm_046NAm_193NAm_194NAm_188NAm_156NAm_154NAm_122NAm_186...NAm_172NAm_201NAm_202NAm_136NAm_137NAm_203NAm_204NAm_108NAm_018NAm_143
1500NaNNaN1.026NaNNaN0.613NaNNaNNaNNaN...NaNNaNNaN0.764NaN0.8140.901NaN0.769NaN
1501NaNNaN1.058NaNNaN0.592NaNNaNNaNNaN...NaNNaNNaN0.928NaN0.9460.863NaN1.060NaN
1502NaNNaN1.088NaNNaN0.633NaNNaNNaNNaN...NaNNaNNaN0.937NaN0.9690.914NaN0.953NaN
1503NaNNaN0.875NaNNaN0.675NaNNaNNaNNaN...NaNNaNNaN0.775NaN1.1090.927NaN1.230NaN
1504NaNNaN1.139NaNNaN0.882NaNNaNNaNNaN...NaNNaNNaN0.834NaN1.0900.838NaN1.203NaN
..................................................................
19961.346NaN1.373NaNNaN0.823NaN1.159NaN1.354...NaNNaNNaN0.7730.9780.6481.018NaNNaNNaN
1997NaNNaN1.153NaNNaN0.564NaN1.358NaN1.134...NaNNaNNaN0.6021.0170.7520.976NaNNaNNaN
1998NaNNaN1.369NaNNaN0.886NaN1.575NaN1.368...NaNNaNNaNNaNNaN0.9741.029NaNNaNNaN
1999NaNNaN1.502NaNNaN0.798NaN1.050NaN0.830...NaNNaNNaNNaNNaN0.7021.031NaNNaNNaN
2000NaNNaN1.674NaNNaN0.720NaN1.237NaN1.254...NaNNaNNaNNaNNaNNaN0.998NaNNaNNaN
\n", "

501 rows × 111 columns

\n", "
" ], "text/plain": [ " NAm_153 NAm_165 NAm_046 NAm_193 NAm_194 NAm_188 NAm_156 NAm_154 \\\n", "1500 NaN NaN 1.026 NaN NaN 0.613 NaN NaN \n", "1501 NaN NaN 1.058 NaN NaN 0.592 NaN NaN \n", "1502 NaN NaN 1.088 NaN NaN 0.633 NaN NaN \n", "1503 NaN NaN 0.875 NaN NaN 0.675 NaN NaN \n", "1504 NaN NaN 1.139 NaN NaN 0.882 NaN NaN \n", "... ... ... ... ... ... ... ... ... \n", "1996 1.346 NaN 1.373 NaN NaN 0.823 NaN 1.159 \n", "1997 NaN NaN 1.153 NaN NaN 0.564 NaN 1.358 \n", "1998 NaN NaN 1.369 NaN NaN 0.886 NaN 1.575 \n", "1999 NaN NaN 1.502 NaN NaN 0.798 NaN 1.050 \n", "2000 NaN NaN 1.674 NaN NaN 0.720 NaN 1.237 \n", "\n", " NAm_122 NAm_186 ... NAm_172 NAm_201 NAm_202 NAm_136 NAm_137 \\\n", "1500 NaN NaN ... NaN NaN NaN 0.764 NaN \n", "1501 NaN NaN ... NaN NaN NaN 0.928 NaN \n", "1502 NaN NaN ... NaN NaN NaN 0.937 NaN \n", "1503 NaN NaN ... NaN NaN NaN 0.775 NaN \n", "1504 NaN NaN ... NaN NaN NaN 0.834 NaN \n", "... ... ... ... ... ... ... ... ... \n", "1996 NaN 1.354 ... NaN NaN NaN 0.773 0.978 \n", "1997 NaN 1.134 ... NaN NaN NaN 0.602 1.017 \n", "1998 NaN 1.368 ... NaN NaN NaN NaN NaN \n", "1999 NaN 0.830 ... NaN NaN NaN NaN NaN \n", "2000 NaN 1.254 ... NaN NaN NaN NaN NaN \n", "\n", " NAm_203 NAm_204 NAm_108 NAm_018 NAm_143 \n", "1500 0.814 0.901 NaN 0.769 NaN \n", "1501 0.946 0.863 NaN 1.060 NaN \n", "1502 0.969 0.914 NaN 0.953 NaN \n", "1503 1.109 0.927 NaN 1.230 NaN \n", "1504 1.090 0.838 NaN 1.203 NaN \n", "... ... ... ... ... ... \n", "1996 0.648 1.018 NaN NaN NaN \n", "1997 0.752 0.976 NaN NaN NaN \n", "1998 0.974 1.029 NaN NaN NaN \n", "1999 0.702 1.031 NaN NaN NaN \n", "2000 NaN 0.998 NaN NaN NaN \n", "\n", "[501 rows x 111 columns]" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "job.df_proxy" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(501, 96)\n", "(501, 111)\n", "(207, 2)\n" ] } ], "source": [ "print(np.shape(job.temp))\n", "print(np.shape(job.proxy))\n", "print(np.shape(job.lonlat))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mLMRt: job.save_job() >>> Prepration data saved to: /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/recon/job.pkl\u001b[0m\n", "\u001b[1m\u001b[36mLMRt: job.save_job() >>> job.configs[\"prep_savepath\"] = /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PAGES2k_HadCRUT_cropped_domain/recon/job.pkl\u001b[0m\n" ] } ], "source": [ "job.save(verbose=True)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimating graph using neighborhood method\n", "Running GraphEM:\n", "\n", "Iter dXmis rdXmis\n", "\n", "001 0.0386 0.0930\n", "002 0.0959 0.2285\n", "003 0.1653 0.3785\n", "004 0.1087 0.2143\n", "005 0.0851 0.1510\n", "006 0.0792 0.1310\n", "007 0.0694 0.1079\n", "008 0.0585 0.0861\n", "009 0.0493 0.0694\n", "010 0.0429 0.0581\n", "011 0.0385 0.0505\n", "012 0.0356 0.0455\n", "013 0.0351 0.0437\n", "014 0.0330 0.0401\n", "015 0.0317 0.0377\n", "016 0.0306 0.0356\n", "017 0.0294 0.0336\n", "018 0.0283 0.0318\n", "019 0.0273 0.0300\n", "020 0.0262 0.0283\n", "021 0.0252 0.0267\n", "022 0.0242 0.0253\n", "023 0.0232 0.0239\n", "024 0.0223 0.0227\n", "025 0.0215 0.0215\n", "026 0.0207 0.0204\n", "027 0.0199 0.0194\n", "028 0.0192 0.0185\n", "029 0.0186 0.0177\n", "030 0.0212 0.0199\n", "031 0.0190 0.0177\n", "032 0.0196 0.0180\n", "033 0.0183 0.0167\n", "034 0.0186 0.0167\n", "035 0.0163 0.0145\n", "036 0.0156 0.0138\n", "037 0.0151 0.0132\n", "038 0.0147 0.0127\n", "039 0.0143 0.0123\n", "040 0.0140 0.0119\n", "041 0.0137 0.0115\n", "042 0.0134 0.0112\n", "043 0.0131 0.0109\n", "044 0.0129 0.0106\n", "045 0.0126 0.0104\n", "046 0.0124 0.0101\n", "047 0.0122 0.0099\n", "048 0.0120 0.0097\n", "049 0.0118 0.0094\n", "050 0.0117 0.0092\n", "051 0.0115 0.0091\n", "052 0.0113 0.0089\n", "053 0.0112 0.0087\n", "054 0.0110 0.0085\n", "055 0.0109 0.0084\n", "056 0.0108 0.0082\n", "057 0.0106 0.0081\n", "058 0.0105 0.0079\n", "059 0.0104 0.0078\n", "060 0.0102 0.0076\n", "061 0.0101 0.0075\n", "062 0.0100 0.0074\n", "063 0.0099 0.0073\n", "064 0.0097 0.0071\n", "065 0.0096 0.0070\n", "066 0.0095 0.0069\n", "067 0.0094 0.0068\n", "068 0.0093 0.0067\n", "069 0.0092 0.0066\n", "070 0.0091 0.0065\n", "071 0.0089 0.0064\n", "072 0.0088 0.0063\n", "073 0.0087 0.0062\n", "074 0.0086 0.0061\n", "075 0.0085 0.0060\n", "076 0.0084 0.0059\n", "077 0.0083 0.0058\n", "078 0.0082 0.0057\n", "079 0.0081 0.0056\n", "080 0.0080 0.0055\n", "081 0.0080 0.0054\n", "082 0.0079 0.0054\n", "083 0.0078 0.0053\n", "084 0.0077 0.0052\n", "085 0.0076 0.0051\n", "086 0.0075 0.0050\n", "087 0.0074 0.0050\n", "\u001b[1m\u001b[32mGraphEM: job.run_solver() >>> job.G created and saved to: ./testcases/PAGES2k_HadCRUT_cropped_domain/recon/G.pkl\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.run_solver() >>> job.recon created\u001b[0m\n", "CPU times: user 1min 26s, sys: 24 s, total: 1min 50s\n", "Wall time: 15.2 s\n" ] } ], "source": [ "%%time\n", "\n", "save_path = './testcases/PAGES2k_HadCRUT_cropped_domain/recon/G.pkl'\n", "job.run_solver(save_path=save_path, verbose=True, distance=1e3, maxit=100)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(501, 8, 12)" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.shape(job.recon)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mLMRt: job.save_recon() >>> Reconstruction saved to: ./testcases/PAGES2k_HadCRUT_cropped_domain/recon/recon.nc\u001b[0m\n" ] } ], "source": [ "job.save_recon(f'./testcases/PAGES2k_HadCRUT_cropped_domain/recon/recon.nc', verbose=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Validation" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Dimensions: (lat: 8, lon: 12, year: 501)\n", "Coordinates:\n", " * year (year) int64 1500 1501 1502 1503 1504 ... 1996 1997 1998 1999 2000\n", " * lat (lat) float64 32.5 37.5 42.5 47.5 52.5 57.5 62.5 67.5\n", " * lon (lon) float64 232.5 237.5 242.5 247.5 ... 272.5 277.5 282.5 287.5\n", "Data variables:\n", " recon (year, lat, lon) float64 ...\n" ] } ], "source": [ "with xr.open_dataset('./testcases/PAGES2k_HadCRUT_cropped_domain/recon/recon.nc') as ds:\n", " print(ds)\n", " recon = ds['recon']\n", " lat = ds['lat']\n", " lon = ds['lon']\n", " year = ds['year']" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(81, 8, 12)\n" ] } ], "source": [ "mask = (ds['recon'].year>=1850) & (ds['recon'].year<=1930)\n", "tas_recon = ds['recon'].values[mask]\n", "print(np.shape(tas_recon))" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(81, 8, 12)\n" ] } ], "source": [ "tas_obs = job.obs.fields['tas']\n", "mask = (tas_obs.time>=1850) & (tas_obs.time<=1930)\n", "print(np.shape(tas_obs.value[mask]))" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "ce = LMRt.utils.coefficient_efficiency(tas_recon, tas_obs.value[mask])\n" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAH9CAYAAAAeU6YWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAADN0ElEQVR4nOzdd3QU1cPG8e/M7qZseq+kQQiEAKH3KkgVQRBRULB3aXZR7L37U6yvCopiQ4qIINKk906AEAKk955su+8fa1ZCAqQn4P2ck6M7e2fm7iZkn9yqCCEEkiRJkiRJVzC1qSsgSZIkSZLU0GTgkSRJkiTpiicDjyRJkiRJVzwZeCRJkiRJuuLJwCNJkiRJ0hVPBh5JkiRJkq54MvBIkiRJknTFk4FHkiRJkqQrngw8kiQ1e3J9VEmS6koGHklqYjt37uTmm28mJCQER0dHWrZsyV133UVCQkKD3fPdd9/F398fR0dHXnzxRQ4cOECnTp2wt7cnOjqar776CkVROHXqVLWuV9PyNfHiiy/y5ptv1vk6JpOJadOm4eLigqurK2vXrq1U5tSpUyiKctGvjz/+2Fb++++/JzQ0FHt7e+6++27Onj1L//79cXBwwNfXlxUrVqAoCuvWratWHdetW1ej8pIkVZ8it5aQpKbz4YcfMmPGDAYNGsS0adMIDAzk+PHjvPHGG2RlZfHXX3/RsWPHer1nfn4+Hh4ejB49mtmzZxMeHs7MmTP566+/+Oabb/D19SU0NJT4+HhbCLqUjIyMGpWvCUVRmDt3Ls8++2ydrrN8+XKuueYann76aYYMGULnzp1xdnauUObUqVOEh4czZ84cRo0aVeV1IiIi8PX1BcDb25vIyEheeeUVgoKC+OSTT3jvvff49ttvCQoKon379hw+fJjo6GhcXV0vWcf8/PwalZckqfq0TV0BSfqv2rRpE9OnT+eBBx7g3XfftR0fOHAgY8eOpVOnTtx2223s2rWrXu+bk5ODxWJh7Nix9O/fH4CsrCzat2/PyJEjbeV8fHyqfU0fH58alW8KWVlZANx6662Eh4dftGzLli3p2bNnta55//33M3DgQNvjwMBAJk6caCtTneuUc3V1rVF5SZJqQEiS1CSuvfZa4eXlJYqKiqp8ftGiReL5558XhYWFQgghTCaT+PDDD0VMTIxwcHAQLVq0EI899pgoKSmpcN6GDRtE//79haOjo/Dw8BC33HKLSE9PF0II8eWXXwrgol9ffvmlrVxCQoLtur/99pvo3bu30Ov1IiAgQNx9990iJyenwnXPLX+xepSfo9FoxNatW0XPnj2Fvb29CAkJEW+88YatzPl1u5BLvTdTp06tcJ0BAwZUeZ2EhATbe3Axa9eurVS30NDQCo/nzp1rK7d27VrbuVu2bBFDhw4VLi4uwtvbW0yaNEmcPXu2wnXPLX/gwAExatQo4eLiIlxcXMTYsWNFfHx8pbr8+eefYujQocLR0VH4+fmJRx99VJhMJlu5srIyMWfOHBEeHi4cHBxEu3btxFdffSWEEOJ///ufAERcXFyF17lgwQKhqqo4ffr0Rd8PSbocyMAjSU3AYrEIBwcHMXHixGqfc/vttwudTieefvppsWrVKvHaa68JvV4vrr76amGxWIQQQqxfv17odDoxfPhwsWzZMvH111+LkJAQ0a5dO1FcXCzS09PFL7/8IgAxZ84csWXLFrFlyxbRqVMn0alTJ7FlyxaRnp5eKcAsW7ZMKIoixo4dK5YvXy7mz58vfH19xdVXXy2EqBx4LlWP8nMURREhISHi3XffFWvWrBE33XSTAMTKlSuFENZwAIjbb79dbNmypdbvzYkTJ8ScOXMEIH755Rdx6NChKq9THni++OILYTQaK32VB4i8vLxKddu6dasYOXKk8Pf3F1u2bBFnzpypFGB2794t7OzsRL9+/cTixYvFTz/9JFq1aiXatWsnjEZjpfJxcXHCxcVFdOvWTfzyyy/ihx9+EB06dBD+/v4iLS1NCPFv4PHz8xPPP/+8WLNmjZg5c6YAxMcff2x7bRMmTBCOjo7ipZdeEn/++aeYNWuWAMTChQtFTk6OcHBwEHPmzKnwfgwZMkQMHTq0Oj+ektTsycAjSU0gPT1dAOKxxx6rVvlDhw4JQLzyyisVji9YsEAA4rfffhNCCNG7d28RExNT4S/7uLg4odFoxP/+9z8hRNWtGAMGDKjQ6nF+gOnSpYvo1KmTLVgJIcT3338vWrduLVJTUyuVr049ys/5/PPPbWVKS0uFg4ODeOCBB2zHyltL6vreVNUKdb7y9+ZCX05OThXKn1+3qVOnitDQUNvj8wPM+PHjRUBAQIVWuc2bN4uwsDCxZ8+eSuVvuukm4efnJ/Ly8mzls7KyhJubm3j44Ycr3OP8sBIeHi5Gjx4thLC2EgHi3XffrVDmuuuuE3feeacQQogbb7xRhIWF2b7HZ86cEaqqioULF17w/ZKky4mcpSVJTUCrtQ6fM5vN1Sq/fv16AG688cYKxydNmoRGo2HdunUUFxezdetWRo0ahRACk8mEyWQiIiKCtm3bsnr16lrVtaSkhN27dzNu3DgURbEdv+GGG4iLi8PPz69C+ZrWo1evXrb/t7e3x8fHh6KiomrXrzrvTU3NnTuXHTt2VPrauHFjja91rr///psRI0bg4OBgO9arVy8SEhKIjY2tVH7NmjUMHDgQvV5vex9dXV3p16/fRd9HgODgYNv7+PfffwNw3XXXVSjz888/8+mnnwJw++23c+rUKdtrnD9/Pi4uLowbN65Or1mSmgs5aFmSmoCHhwcuLi4kJiZesExRUREGgwEPDw+ys7MB8Pf3r1BGq9Xi7e1Nbm6ubTDya6+9xmuvvVbpeo6OjrWqa3Z2NkII28ykS6lpPfR6fYXHqqpisVhqVD+4+HtTU2FhYXTt2rXG511KVlZWtd/H8vKLFi1i0aJFlZ47f5D4xd7H8gHbF7v34MGDCQ8PZ/78+fTv35/58+czadKkCuFMki5nMvBIUhMZNmwYa9eupbS0tMoPlc8++4zZs2ezY8cOPD09AUhNTSU0NNRWxmg0kpmZibe3N66uriiKwsyZMyu1dkDlD8TqcnNzQ1EUMjIyKhwvLS1l7dq19OjRo8LxhqrHhVTnvWku3N3dK72PACtWrKBTp05Vlh8yZAizZ8+u9Fx5K2F17wvW5QOCg4Ntx48ePUpWVhZ9+vRBURSmTZvGe++9x7333ktcXBxff/11te8hSc2d7NKSpCYye/ZssrKymDNnTqXnUlNTefPNN4mOjqZz584MGDAAgO+++65Cue+//x6z2Uzfvn1xcXGhc+fOHD16lK5du9q+2rVrx9y5c2u9mJ2zszOxsbEsW7aswvHff/+dkSNHkpycXOF4fddDVS/+a6o6701z0a9fP1atWoXBYLAd27NnD6NGjapy+YEBAwZw+PBhYmNjbe9jly5dePvtt1m8eHG171v+Hpz/PXzssceYPn267fG0adPIzc1l9uzZtG3btlKYlaTLmWzhkaQm0rNnT1544QXmzJnDkSNHmDp1Kt7e3hw8eJA33niDkpIS2ziN6Ohopk6dyjPPPENxcTH9+/dn7969PPvsswwaNIjhw4cD8PLLLzNy5EgmT57M5MmTMZvNvPnmm2zbto2nn3661nV9/vnnGTNmDDfeeCNTp04lNTWVJ554grFjxxITE8POnTsrlK/Peri7u7Np0yY2bNhAv379Kowjqsl7UxPx8fFs3bq1yuc8PT1p3bp1ja8J8PTTT9OrVy9GjRrF9OnTKSkpYc6cOXTv3p2rr76azZs3Vyj/zDPP0KtXL0aPHs29996Lg4MDn3zyCb/++is//fRTte/bsWNHrr/+eh555BGKi4uJjY3l999/Z9myZfzyyy+2ciEhIQwdOpQ//vijyu5ISbqsNfGgaUn6z1uxYoUYOXKkCAgIEPb29qJVq1binnvuqbT2iclkEi+++KKIiIgQOp1OhIWFiSeffLLSOjx//vmn6Nevn3B0dBRubm5i8ODBYuPGjbbnazNLSwghli9fLrp16ybs7e1FcHCwmD17tm2NoKrKX6oeF5o1FRoaKqZOnWp7/NZbbwl3d3eh1+tFYmJile9hdd6b+pilBYhrr73WVp4aztISwjrVfuDAgbb1cm677TaRkZFxwfK7du0Sw4cPFy4uLsLZ2Vn07NlTLFmy5KL3EKLy97SsrEw88cQTIjg4WDg4OIhOnTqJn3/+udJ78P777wuNRiOSk5Mv+D5J0uVIbi0hSZIk2ZTPIqtJl5kkXQ5kl5YkSZLECy+8QFxcHH/88YdtGrskXUlk4JEkSZJYunQpJ06c4I033qB3795NXR1JqneyS0uSJEmSpCuenJYuSZIkSdIVTwYeSZIkSZKueDLwSJIkSZJ0xbtiBy2fPn2azMzMpq6GJEmSJEkNyNvbm5CQkEuWuyIDz+nTpyvsqSNJkiRJ0pVJr9dz5MiRS4aeKzLwnN+y4+DggLe3N15eXuzbt69S+bCwMCZcP55uPbrVaEM+gJPxCTz1xFO0bt2a2bNn06FDh0ue89577zF//nwAbph0A+PGX4tOZwfARx/MY8eOHUyYMIGYmBj69et3wb2Etm7dygMPPEBtJtoFBQXRunVrYmJiaN++PdHR0bXeTbs+zJgxg507d3L3PXfRu28vNJor8kfzP6OsrJTUlDT8/P0qbYx65PARXn35NQoLC3FwcMDHx4czZ87g6OhI3759aR3VmsW/LCY5OZkuXbpw3YRxeHp5snLFStb+tY61a9fWa10LCgqYMmUKZ8+exd/fn+DgYBwdHYmLiyM9PR2Ap556ir59+150t3EhBPPnz+ezzz6jpKQEAI1GQ2RkJP369yW6XTuCggPq7WdbCAvHj8Wz6e9NLF+2HLBuw3HXXXcxduxY7O3t6+U+zU1mZibDhg0DrK83NzfX9pyvry+eXp7Y6ewwmU2YjCZSU1PJz89Hr9cTGxtLVFQUQUFBCCGwt7fHzs4Oe3t72/tVWFjImTNnSE1NJTMzk6SkJMrKynByciI0NJRWrVrRunVrWrduTUlJCWvWrOHPP//kyJEjtnp4e3tjZ2eHm5tbheM+Pj6Eh4eTmJhIWlqa7bhOp+Oxxx6jTZs2hIWFVfpdbLFY2LZtG0uWLMHBwYHhw4fTs2dPysrKePbZZzly5Ah6vR5XV1eys7OJj48HYPKUm7huwnVkZWaSlZWDm5srnt6e6LS6ev++NJWT8fE8MvsxMjMzLxl4rshp6bt376ZLly48/uRjDBk6hMCgAFtoOHjgEPff8wBGo5GrrrqKsdddS6cusZfcoPBiDh86wjNPzeXgwYOMGDGC559/nq5du1ZZduHChUyePJmYmBgW/vANdnZ2FZ7PSM/ktZdfZ/ly6y+wzZs306tXr0rXEUJw9dVXk5WVxedffYbJaESn1SEUMBqMmC1mtFotWq0Wi8VMUWExRUVF5GTnkJaaRnx8PAcPHGTXrt0UFRWh1WoZMGAAN9xwAxMnTsTNza3W70dtFBYW8sADD/D1118z+prRjBg5nMFDBlXaN0lq/oxGIyOGjuLs2bMAjBw1koEDBxLVJhJPLy/efes9fv75ZxYtWsTYsWOxs7Pj5MmTfPPNNyxevJh9+/bRpk0bpk+fzlNPPUVWVpbt2vfffz//+9//6qWe5XuVrVq1ig8//JBPPvuEAYP62Z43mUzM/3IBr7/2hu3Y1VdfzfPPP3/RTTULCwuZO3cuP/zwAy4uLqiqyokTJygrK8PBwYGYmBjaxbQjomUE3y74lry8PB586AHGTRiLRqOpdv1/W7qC2bMervK5wYMH8+abb9KhQwdMJhPr16/n+++/57fffrOFuAULFjBlypRq3685iY6OxmAwsHzlUjLSMzGZjHh4elT5R5vFYuH4sRNs3byVXTt3c/jwYdvPZnVMmDCBH3/88ZLlEhIS+Pnnn1m3bh1r1qyhtLTU9lxQUBBdu3ZhyZKltmNarRZfX1/b5rvr1q2zbYQL1p/PDz74gBUrVrBnzx7y8/Mr3O/kyZOEh4dXOCaEoGXLliQkJAAwfvx4FEWpsO9a+/bt+XHxomq//ubu0MHDjB87gV27dtG5c+eLlr2iA8/Pv/5Eu5joSs8LIer9g9RisbDur/W8+/Z7HDt2jEGDBnH33XczcuRIXFxcbOW2bNnC0KFDcXd359vvF+Af4F/hOgcPHOKJR58kISGBhx56iNdff73SL8EFCxYwb948tmzZwmdffEq/AbXfDdpsNpN46jR7du9lzeo1ZGdnExcXx3XXXce0adMYPHhwncJgTe3Zs4dnnnmG5cuXM3jwYN5+/81KLQRSzQkhSElOocxgICgosFLQrm8fvj+PD97/4KJlZs+ejZ2dHaWlpbRq1YrJkyfj5uaG2WxGVVWOHz9OVFSUrfygQYNYvXp1jULBxbz00kuVdqq/8847mPXoTNvvB6PRyIKvvkFVVfz8/fjis//j4MGDdOvWjQEDBjB48GAGDx5coTVl7ty5PP/88xWuq9VqufPOO4mMjGTPnj3s2bOHI0eO0LFjRzp06MDXX39Nt27duPX2afj6+iIQnD1zlo3r/yYhIYE2bdvw4PT78fD0sF2zqKiYH7//EXd3d9w93AkMDMRgKGP2zEdITEwErC0HRqMRsDb7FxcXV6hXhw4dcHZ2pqSkhOPHjxMWFsZXX31FREQE7u7uzfYPjnHjxvHrr7+yc+92nJ2da3x++c+YyWDCYDRgMFi/hIC9e/Yy46GZtrJTpkxhwYIFNbp+QkICcXFxZGRkcOzYMV599VVMJlOFMsOGDcPDw4Pu3bszYsQI2rRpY3tu165dF/yj2cfHh2uuuYb3338fJyenCs/l5uby2Wef8cUXX+Do6MjBgwdxd3e39XrExsYyY9Z0eva+cGC/3MjAc4nA05BMJhMb1m1kwdffsGXLFnQ6HbGxsXTq1ImWLVvi7OzMokWL2LBhAy++/CITJl5X4dzunXtSXFzM2LFj6dChA6mpqZjNZttXfHw8mzdvZujQoYwbP5bBQwbVa/0z0jP44/dVLPr+B44fP46vry/9+/cnKiqKiIgIAgMD8fX1xc/PDx8fnwb74Pzjjz+47rrrUBSFzp07ExvbkXbt2xEaHkpAgL8MQTWwZ9denpv7PEePHgXAycmJp+Y8xXXXj23Q+5aUlPD78j947dXXMBgM3HPPPcybN8/2l6+npydubm6oqkJ8/El69uzJli1bbOebTCZ0uopN74WFhZV+ydfWM888wwsvvABAnz592LRpEwAbNq3H18+nynPMZjMbN2xi5YqV7Nyxk7Nnz+Lq6sro0aOZPHkyI0eOtH0Ylxsz5hqWLl3GtddeW+G40WhEq9WiKAp///03M2bMYNeuXRet8/bdW3F1db3kazMYDMQdiSMu7jjOzk6ER4QTHhHGwgXfYTKZ8PT0Ii4uDkVR+H3F7xW6V8pNnz6dd99995L3agrz589n6tSpLPttKZFRrer12iuWr2TWjFm2x3q9nmeeeYbHHnus1tcUQvDXX3+xZs0aMjMzeeaZZwgODr5g+QULFnDLLbfg6emJVqtFVVUCAwN55JFHuOGGG6oMorfcckuFYObn54dWq6V79+507NiR5557DiEEHh4etGrVis5dOtOjZ3diO3dEr9fX+rU1NRl4mjDwnCs1JZUtm7axd+9e4o7GcfLkSYqKimjTpg3XjR/H2PFjcXL69wdNCMHnn3zBHytXkZ+fj9FoxM/PD51Oh0ajQaPR4OnpwZixY+rUqlMdQggOHzrKhnXr2b1rD6dOneLs2bOVxgu5u7vj6+tr+woJCaFHjx6MHj26Vn95nSsuLo4lS5awfft2tm/fzpkzZ2zPubm54ePjjYuLK46Ojjg6OmLvYI/FYsFsMmM0GrFYLCiqgqqoaDQaXFxdcHd3Jyg4iNDQEFqEBBPcIrjG47YuJ6tXrmbG9Fn07duXhx9+GDc3N5588kmOHz/OXxv+rLfWkgsRQtC+bUfbX7c9evSgb78+dOvejQ6x7VFVlf179zNxwiQCAgJ47733mDBhgu0X+q5du3juuee46qqruOOOO+ol7BgMBjZs2MD//vc/lixZwsSJE3n+5WcRQmAwGKo99kUIwalTiaxbs44Vv/3OgQMHmDZtGlOnTmXSpEk4ODjYWloiIyP56aefLjrGLz4+nlatLv7hvWHTOnz9LjyOqDr+Xr+JO26/0/ZYVVUsFkulcm+88QYPP1y5yywxMZGPPvqIgwcPYmdnx7XXXsuUKVMa9d/RkiVLGDt2LEuXL6F1m8h6v77RaCQ1JY2zZ87y3bffs2rVKmJiYvD19UWj0eDr60ubNm2IjY2lS5cuBAQE1Ov9S0tLueOOO/j2228Baxfe4cOHcXNzIycnp1LgKSwsxNPTk6FXD2XqtFtIS0tn8c+LbePdjh49ikajYfPmzSQlJXH48GE2btxIYmIi3t7efPTx/+gQe+nxp82RDDzNJPBcaYxlRnJyc8nNySU3N4fs7Byys7LJzsomIzOTrMwsTp1K4MSJeBwdHRk1ahTXXHMNQ4YMITAwsM73T09P5/jx48THx5OWlkZ6ejqFhYUUFxdTXFxMSUkJGo0GrVZrC4lCCMxmMyaTiby8PDIzMzl58iQFBQUAXD3sat7/8N061625ObDvAJ9+/Dl//vknU6dO5bPPPrN9IH388cfce++9uLu7ExMTg5OzE0aDEUe94z8B1ofQsFCCggJxcnbGy8sTN/fajenKzcmlZzfrvkxLli0mqm1UpTLFhcVMf3AmR44cITMzky+//JJp06ZV+x55eXmsWbMGIYStVRAgJSWFjz76iPT0dAoKCsjPzyc1NZWjR49SVFREZGQkU26ezOhrR1f4w6O2Vq74g9kzH+bmm2/miy++QFVV0tPTURQFb2/vS3YP/fXXX1x11VWEhoby9DNz8A/wY/TIMbi4uDBz1gyuHTcGJ+e6B76ks8ncesttnD592nasX79+bNy4EY1GQ79+/WjRogX/+9//KrQmCSH49ttvue+++7C3t6db924U5Bfw999/Exsby/PPP8+oUaMapQu8uLiYdu3akZKSQnh4OEFBQVw9/GoGDxlYYQhBfTh7JokbJ96Eu7s77dpFI4QgPT2Do0ePkpOTA8CLL77IU089Ved7HT16lH379vHFF1+wdu1a/P39CQwMZPv27fj6+vL8889z9913VzqvtLSUmJgYioqKGDlqJHm5eSxdupSePXvy1ltv0a1bt0rnCCHYsmULffr0oW3btixe9nOd698UahJ4rtw/baV6p7PX4evnc8Hm/nKpKWmsWb2GFSt+Z+rUqQCEh4cTGxtLdHQ0kZGRhIWFERYWRlBQULX+Mty3bx8mkwkfHx8yMzPx8vJi4MCBtG7dusYDrK2/sNL59ttvmT17NnFH4qr8IL4cGAwGiotKyMrMJC0tndSUNP5c/Sd//fUXsbGxfPbZZ9x6660VPoTuvvtuunTpwsqVK20f/q4urhQVFRF/Ip51a9eRlJRkK+/o6Mi6jX/VKvS4e7iz/u+1eHl7XfD7bBZmHprxANu37eCN19+sVgvLd999x4IFC3B3d+fUqVO2rrA333yT2bNnI4Sgf//+nDhxgs5dOuPs7IyzszMxMe0YNXok3Xt0o2Vky3odozJ85DA0Gg0PPTCd0tJSPv/884vO6jpXVlYWe/bsAawtKDt27KRvvz7M++Qj+vbvU6lrry6CggNZ9ddKcnPy+HP1Gn75+Rc2btzI3Llzeeihh/D09Kx0Tn5+PpMnT2b58uXcMOkGHn38YVv4ijsSx8svvsqYMWMIDQ1l/PjxjB8/np49ezZY+NHr9ezfv59PP/2UM2fOcODAAR575DH6D+jPp198XK/3Cm4RxMYt6yscS01JY+2atbz6ymuUlZURFxdXp3ukpKRw/fXX27pVg4ODefb5Z1myeAk7duzg3Xff5Z577rngvw0HBwdWrVrFU089xc4dO3FxceGxxx7jqaeewmKx8Oijj7JmzRqKiorw8/PDy8uLgoICNm7ciIuLC/c/eB8FBQW2WWtGo5FfflzM3GeeBWD2w7O485476vQamwPZwiM1qJzsHPbtO8De3Xs4Fneco0eP2mYlgHXabkREBDExMRgMBhwdHfH29rYtI6DVatmwYcNFZ0mEhYXRr18/unXrRkxMDJ07d65WCLJYLLi6ujJz1gymTJ1cL6+3IQghOHzwCCtXrGTv3n2cOXPG1qp1/kBIrVZL165dmTVrFuPHj6/wgVNaWsrevXvZsWMH+/btIzExEUdHR1q2bElsbCyDBg2yTessKiriySef5P333wdg8pTJ+Pr64OXlRW5uLtu2bgcgtlMs7dq3w9vbC1cXF7y8vS7YCmE2mzly+Chnz5wlIz2DvXv3sXfPXlu40uv13HTTTXzyySeX/KBs2bIlycnJREVFVVhq4sEHH7TVedCgQWzbto377r+XmPYx/PrLr2zfvoPIyEjeeu8NW0tAcWExikapt2UZ1q3dwKzps/Dx8eGxxx6jZ8+ebN26FSEEt912W6XxZ19//TX33XcfBoOBESOGM/nmybh7uDF86EgABgwcwIcff9BgXUZCCD77+AvefuttbrnlFubNm1dpTMesWbP47LPPeOe9ty/YnX740BGWLVnOit9WkJaWRqtWrXj66aeZMmVKg7b6HDx4kGuuuYZTp07x4EMPcv9D9zbYvb7+v/m88vKrtsd2dna8/vrrPPDAA7XuHj527BjdunUjPz+fu+66E52djuPHTvD333/TokULPv30U/r371/rOn/xxRfccccd3DT5JtzcXMnMyCInJwd7B3tatWpJWHgYa9esZenSZSiKQmxsLBkZGRVmsk2degtPPP14revQkGSXlgw8zVpZWRnpaRmkpqSSnJzCiWPHSUg4hZ29HWWlZWRnZ5OVlUVmZiZms5mysjLbud26deOV11/GwdGBnKwcEhMTOXjgINu37eDgwYMYDAY0Gg09e/Zk2LBhXHvttRccN2EwGHBwcOCBB+7n/un3NdbLvySTyURKciqnT51mz+49rFy5khMn4gkNDaV///5ERkbi5OSEXq9Hr9fj7OxMcHAwLVq0wNfXF1VVSUxMZNeuXbavEydOcPr0aUwmE3Z2drRr146QkBaUlJZyMj6eEyes63aMGzeOL774Ag8PD7KysnjxxRc5ceIEaWlppKWlkZqair29PcOGDUMIwbZt2ypN8XVxccHf35/S0lJycnIoKSlBp9NhsVgwGAyAdexXhw4d6NWrFw4ODhw4cIBt27bRrl077rnnHsaNG3fB96e0tBRHR0d8fX154qnHeemFlyksLOSbb75h8ODBeHhYZzLl5eXx6KOPsmjRIvLy8gDw8PCwjXeIjIy0hUBnZ2fmf/M1kVGt6qXV5+zZJD5454MK05ABHnjgAT74oOLstRYtWlBQUMCTc55Eq9GQlJTEr4t/BRTeffddxowZwzcLF9C1e5c61+tChBCMHX0dcXFxrF69miFDhties1gshIWFMWLkCGY9MuOS1zKbzRzcf4hvF3zL0qXL6NevH3/88Ue9BUqDwcC8efPIz8+na9euPPDAA5w8eRKA5194jok3Xl8v96nKQ/fNYNWqVZWOx8bG8vvvv+Pv71/FWRd3+vRpWrVqZZtN5+bmRkxMDGPHjuWOO+7A3d29TnV+8803eeSRRzhwyPqHwdGjR/l74ybW/rWO/fv3AxASEsITTzyBqqps2rQJV1dXunfvzi233AJYg93u/TsrhO7ymW5NPZNPBh4ZeK4oJSUlxJ84SV5OHu3at8Pdo+rWG5PJRPLZZPbs2cffGzaybt16CgoK6NOnD4899hijR4+u9I/z8ccf57XXXmPqtKk8/MgsdPb1vyCXxWKhoKCQgoICykrLKC0ppbS0lJLiErKzc8jKyiQ9PYO0tDROxp8kPj7e9svPz8+PMWPGMHHiRNsSAadPnyYlJYWcnByys7PJzs4mLS2NxMREEhMTOXTokG3tmsDAQDrGdqRVq1YEBQfRNroNEeHh6Ox11iUJEs9wJvEMWq2G1NQ03nrjLWJjY/nzzz+r/EUmhEAIUeEv9oyMDDIzM22LtJUv2qbX6/HwsK6NYjKZUBSFrl27Ehsbi4ODA2lpaRU+IIYOHcrq1atRVZXS0tILduMIIRg9ejTbt28nMzMTvV7Pyy+/zPTp06ssbzKZSEtLQ6vV4ufnx+7du/n2229JT09Ho9EQHR3NK6+8Qm5uLm5ubnTo0IEuXTvTo2cP2neMqXHLSm5OLvv2HuCdt96xzYwD6NSpEx9++GGldbUmTJjAzz//O34iICCAHj168Oqrr+Lg4EBYWBgzZ81g5KgRtAhtUaO61MSH789j3kfzMBgMFb73r776Kk888QQ/Lf6RmPbtanTNPbv2cuvU25g1axYvvfRSneu4fPly7rrrLlJTU3F1dbUF2XJ33nkHsx+bdYGz65fZbKawsIiVv61k7jPP8sknn3DXXXfV6lqpqans3r2b0NBQ2rZtW68tYjfddBPfffcdrVq15PTpMxgMBnx8fBg2bBjDhw+3jdkCeP/99/m///s/LBYLHTp0wGg02lrXb7zpRh6a8QAenh5s2riZ22+1dnGFhoby5tuv075j+3qrc03IwCMDj4R1psXmv7ey4OsF/P333/Tv359PPvmkwnoXAJ9//jkPPPAALVu2ZNbsmbSPbY+jowNCCEpLSsnPyycnJ4ecnDwMZaUYjCZMpn++jEaMJhNGg4Hi4hJOJyaydes2srKysLOzw9nZmaysLMxmc5V1VFUVX19f/P39CQgIoGXLlrRt25Y2bdrQtm1b/P39bR8+BQUF3HvvvbaZG+Xs7Ozw8fGhRYsWBAcH07JVBG2j29K6TWSl8Rh5efn8seIP/t64iS1bttgGb4P1r9TZj8zi5sm3sGTJEsaMGcPGjRuZMmUKxcXFODs74+npydixY7n11lsvOq22OnJycujbty+HDx+ucLx79+5s3br1kn85lq8+26ZNG1urTm1lZ2ezdetWdu3axdKlSzl69CiFhYX4+vrSp08fwsLDcHGxrlejqiqtIlsR2ToSP39fFEWhtLSUndt28cOiH1m9ejVCCAYNGsRTTz1FVFQULi4uF+1mzcrKQgiBu7u7LWAVFBSwatUqJkyYYCu39+DuBluS4cP357Ho+0WkpqZWGPgdGBjIffffy0MzH6zVdd9+/R2WLl1GYmJinZex8PT0JDo6mjnPPkWLFsFkZmTh6OiAs4tzg8w4NJlM5Ofnk5KUiqOTnrCwEMrKyog/Hs+B/QdZ+fsfbNu2jZCQENatW1dpIcDm4MSJEyxYsICioiJCQ0Pp1asXnTt3rhSqVq1axbBhwxg79lp8fHxYtWo1ZWVl/PHHH2zcuJGnn37auqjoiBGYzCYW/7IYAGdnZwYOHMib777eFC9PBh4ZeKRzCSHYsW0nLzz3IklJScyZM4cZM2ZUaGLft28fDz74IBs3bqzWNVVVRafTYWdnZ/uvs7MzgYGB9OvXj6CgINs4Gx8fH+uS956eODk5VeiO8vLyqtYv6r1793LDDTeQnJzM3Gfn0q59NC4uLri4OmNvb3/BcFBWVkZiwmkMhjIcHR25bdod5ObmEhYWRmhoKL6+vsTFxbFjxw6m3TqNx596lHHXjCcgIIDvv/+ezZs3M3bs2ErXbdGiBYmJifXSnG2xWIiPj2fbtm2oqsqNN97YKM3k+fn5xMfHExQUZBtcvG3bNnr27FlleZ1Oh6enJyaTydaC5ubmhp2dHTk5OZhMJrp168Zdd91l63qszevYtWsXr732GosXL8ZkMuHh4UH79jHcdsdtDbpg3LhrxpOWlsamTZto06YNZrOZuXPn8tJLL7Hw+2/p3LVTra6bmJDI8KtH8thjj/HKK6/UqY7+/v5MvGFig47TKTdtym1s3bq1wjFvb29rS22Z9d/ToEGDuP322xkzZsxlv7xF+SzBoKAghgwdQq8+vbjnzntsi4Pm5OTw7rvvsnr1alJTUyktLaW4uJj8/Hxef+N1xowbDWBrzW3oJS/KycAjA89/XkZ6JmlpaTg6OODg6IiDowO5ObmMGm79R/nnn39y1VVXVTrv8OHDHDlyhKKiIjQaDXq9Hm9vb3x9ffHx8UGv19umvDeWHTt20L17d1q3bs1Hn35IcHBQhecNBgPff7uIv//ehJOTEy2Cg3F1cyXh5ClWrFhh29cJrK1BvXv3ZsOGDdZ1ihSF7t27M+mmGxg24mpUVbUNurVYLEyaNAmNRkNOTg5paWns37+f0tJSbrnlFj7//PNGew/q2/79++nYsaPtcfv27dHr9ezdu9c2Zmzus8/g4emBv78/J+MTmPfRPDIzM/nf//5HTEwMR48eZf/+/Rw5cgR/f38efvhh2rZtaxtHVtOfkby8PJ5//nneffddIiIimDrtFnr17UVQUGCjBMA/VqzipRdfJjc3l7vvvpu1a9dy8OBBHnn0YW6789Ya1WHF8pV8v/B7rp94PUOGDea5Z15g3959trE2VUlJSWHt2rWkp6czYcKESi2I69atY9CgQXzx5ef06de71q+zul5+/hXmz1+AVqtl9+7dpKens3HjRtzd3Rk0aBDR0dH1OnsuLy+PF154gYULF6LT6fjwww8ZPXp0na4phODQoUNs3ryZgwcPkpGRgbOzM4WFhZw8eRKDwYCnpycRERG0bt0aR0dH9u7dy7p162z7cYWEhNjWkzrfihUrGDVqFAMHDqRzl84c2H+A1atX07t3b/5vfuP8fpCBRwae/6SSkhKOHDrK7l27+ejDeZWW0S83cuRIli1b1qhbZtTFp59+alt7Y9ToUQwffjV6vRN5+fmoqsLG9X/bxoCMGjWKU6dOkZubi4+PDzfccAMDBgxg/vz5fPzxv9N1p894iGvGXoOnp0eVq6zm5eWz+KfFLP7l1wpTbl1cXFi/fj2dOtXur/3morCwkFatWpGWlkZkZCTde3THaDDSLiaaFiEt/vmq+IFbUlLCKy++xg+LfqjymnfffTcbNmzgyJEjTJo0ie+++65adTl58iTvv/8+8+fPp7S0lEcefZiJN17fJC0GBoOBD9+fxycff0JwcDAfzvugxks25Obk0bNb5f3/pk2bxpdfflnlOb/++is333yzbUNZIQQ//fST7QO/rKyM2NhYXF1d+eqb/2u0f7vjr72eQ4cO8eabb2JnZ0dhYSFPPPFEg9zr119/rTBYv23btpW6e2siISGB2267jXXr1qHRaGjZsiX+/v4UFRfh6OBIaFgodnZ2ZGRkkJKSwtEj1iUqwLp9Rfk0+IkTJxIZeeHFHX/++Wfee8+6pVLHjh1tg7o/+PB9hlx9VYOHdRl4ZOD5z0lMSGTihEnk5eXh4eHBuHHjmDlzJiUlJRQVFVFYWIhOpyMkJIQ2bdo0+cyC85nNZubPn8+WLVt46aWX8PGpuNbR8ePH+fXXX/nll18qNbOXGzduHP3796ewsJDNmzdz9OhRzGYzaWlplJWVMXr0aLQ6LWVlZTz34txqbVFQXrdtW7ZzMv4kL77wEv7+/jz55JOEh4fj4uKCwWCgpKSE9PR0Tp8+zc6dO9m/fz/9+vVj6tSpXH311c0uXP71118sXrzYthHpd4sW0qlL7CXPE0Jw5PBRiotLeO3l1zhw4IDtubZt21JcXExiYiKLFi1i4sSJtnN+//13Dh8+THp6OhkZGbRv3x53d3dWrlzJzz//jKenJzdMmsiNU27E29urQV5zTSSdTcbZxRk3t+r9jJxLCMHwIdYtNsaOHcvZs2dxc3Pj6quvrnIMT0JCApGRkbRu3Zr5C79Co2i4ZtS1DB06lC+//NK2jsz777/PshVLCAsPq4dXWD0mk4mPPpjHRx/OA6yzC8sXG6xvpaWlzJo1i7Vr19KvXz+eeOKJWo8J+vbbb7nnnnvw9vZmztNP0a1n10uub2WxWEhPSyfh5Ck+/ugTiouL2b9/P7t27WL+/Pm2yQn29vaEhoYSERFBdHQ04eHh2NnZ4ejoiI+PDykpKcycOZOffvqJYcOH8fpbr1Z79fLakIFHBp7/DGOZkUWLfuS9d96joKCAP/74gyFDhjSrD1ghBEePHmXlypXEx8eTkZFhXR/mvvt45JFHUBSlwoylESNGsHbtWlxcXHj00UeZPn16habz5OTkf7YZsQ5IHjp0KOvXWxdGc3BwwNHRkS5dutA2ui06nRZPT086dOxQL0vwb1z/N/O/WsC2bdtsU8zLKYpCQEAAbdq0oU3bKNav28ChQ4eIiIjgnnvu4e677652yGpoYWFhlZrpP/50HlqNtfvixIl4SopL8PbxJja2Iz179SAkLMQWlIUQXDVgKMnJybRu3Zp33n+blq0imDX9YZLOJrFnzx5b2Xnz5nHffffh6uqKh4cHXt5eHDl8hJKSEtq1a8fNt0xh+KhhV9T+cKOHj2H06NG88847lyyblZWFt7c3rq6ubNu1BUVReH7ui+zds5dPP/2UGTNmsHPnTp6Z+zQ33XxjI9S+si2btnLr1Nt45513mDFjRpPUoTrMZjOzZ8/mvffe44ZJN/DY44+gd67ZKuJGo5HRw8fQpk0bVq5cSWxsLPv27aNfv354eXlSWlbGmdNnOHnyZIXuclVVcXBw4MMPP2TatGksX76ciRMn0qdPH15/+9UG269LBh4ZeP4zflu6gtmzrPv9BAYG8sorr6AoCqqq4u/vT1RUFEFBQQ3eolNWVsauXbvYs2cP3t7etG7dmkOHDrFhwwb++OMPTp8+jb29Pa1bt0av1xMcHMzPP/9MREQEM2fO5P777yc5OZlvv/0WT09P7rzz372OVq1axdChQyvdUwjBCy+8wNy5cwHYumPLBafs1zeTyURubh6lJaXodFrs7e1xcnaqEMzK92P74fsfWPzLYvR6PZMnT2bgwIFER0fTqlWrGs3ayc/PZ//+/SQnJ2M0Gjl9+rRtnIGdnR2urq62L29vb4KCgmjRogWhoaGVvv/PPvsszz33XJX3CQsLo3379jg7O5OUlMT27dttS/c/9fSTtpag8rE+5X+9pqWmMaj/VcybN882PTkhIYH27dsz4foJPDHn380ny9ckupJCzrmeeepZ1q1dx5EjR6pcufl8rVq1wsvLiy++tm6B8vQTc1m6dCkGg4HOnTvz5JzHiekQ0wg1r5rJZGLKpFs4deoUr7zyCrfffnu9jt+pL7feeivz58/nxZdeZNyEa2v1e89gMNClYzfeeOMNpk+fznPPPcezzz7L7XfcxqCrBuPp5YmdToeqKqSmWNfnMhqsW/d89ulnODk52cZqrVmzhrFjxxIQEMD7H71HeAO0zsnAIwPPf4bBYGDLpm3s2rmLFb+tqLQIHlhnVkRGRuLu7o6dnR0BAQEEBgbaBiCXry0jhMDZ2RlfX18cHR1RVdV2vKCggA0bNnDq1Kl/PuxzycnJITc31zZrx2KxVNqIMSoqigEDB9C3X286xHao8AG3Y9tOli1dzg+LfqBXr154eXmhqiqqquLl5WX7wJ4yZUqVH4zlM4rGjLmGWY/MxD+g5oueNZbMjEy++/Z7Vvz2OwkJCYB1VeiePXui0Wg4e/YsdnZ2BAcH07FjR8LDw8nKyiIuLo5jx45x4sSJSl0Jer2eNm3boKoqhjIDhYXWtY7y8vIqtD55e3vToUMHQkJCUFWVM2fOsGnTJtsYLxcXF1586QWmPzQDsH5oTJgwgb59++Lq6kpJSQlLlizhxhtvpHv37sxf+FWVr/GH737kxRdeIisry7Zx7rXXXsvKlSvZsmNzvezXdbnIyspm5NWj6NixIwsXLiQoKOii5R977DFef/116+QAXx8OHTyEk5MTL7z0AsNHNo8u0cLCQt587W2+/+57wsPDmTVrFpMnT67zkgj1ZePGjfTv35+33n6TUWNG1ulaN4y/kW7duvHZZ58hhOCdd97hhRdeIDc395LnRkRE2P4QAetK0hMmTCApKYkFC+fTslVEnep2Phl4ZOD5zyorK0NFxWQxkZ2Vw6mEUxyLO8bp02coyC/AYDCQmppqm1ZpNBpRFMX2VVRUVGlH+HJhYWG0jmqNTqezravi6uqCVqvFw8OTqDaRtIxsSUZ6Bvl5+fgHBlRrDMRff67lj9//wGgy2XZ7T0lJ4eDBgyxatIiOHTvy3nvv0bdvX2666SYAFixYwLRp01BVFR8fHx546AHGXjemUWeP1VZuTi6Jp05zLO4Yu3fvQaNqCAjwx2gycvZMEocOHeL06dN4eHgQGRlJRMsIwsJC8fPzpVXrSPz8fFE11ubzCy2OWFJUQlZWFsnJyRw8cIjjx46TnJyCQODn60ds51hiYtrh6u6Gn68PTs5OFBYW8suPi/n6q/kkJSWhqipt2rTB2dmZPXv2YDQaL7iS7/atO7jvnvsZNmwYP/30k+34lClT+Pbbbxk9ejQvvvr8FduiU5V9e/Yzc/osSkpKePHFF7nzzjsv2iqyb98+Fi5cyMaNG9m6dSvLfltCq9YX3z2+KZw8eYrPPv6MZUuXodVqufvuu3nllVfqbSXp2nrttdd4/PHH2bV3R503mZ025TZatWrFggULbMfMZjMJCQlkZWVRVlZm+90J2FrVnZycCAsLw9vbu8L1cnJyGDJkCElJSSxZvhh3D/c61e9cMvA0UOBpP+GNeruW1DwJYQFTGQgzlP/TUBRQNKCt+gO2oZiOLQMUhKkUjMUgzKj+nVB922NOWIMo+GeDT+dAKExGcY9AGzag0erXkIQQjT6w3N7F03ZvUZqLOe8MlsIMhNmA6uKPe9ur0Dp7VzpPCEHqzzMQ5jLa3Pc9qp3jOc9ZyD/2N8mr3setzUAChz4AgLtv9T6Qgv2dq1WuTUDzGBt1vpKCXFZ9+jr7/vwVz8AQBk+bSdu+V9u+t618Kr6+ovxc3n/kbgylJcz9emlVl2wyJzIKKzwuzM7gjx8WcHTFfNyCWzHw0Y9Izam4t11uelGl6+SmpFR4nJ90rMLjsoLsWtVPGIowHf4RNaAzGr+qt9OpDkt+EuaTq1GDeqDxaVvr61Sqn7EYU9wSRlw9iLffe7Pe/n3XJPA0fTuhJDUjiqKi6BxR7JxR7F2sX3bO1mON/AGs+nZAFGeCsRht1BhUvw5YUvdiOvoLoiz/3zpjQfWKQhScvWDr1OWmKWfRKYqC6uiBzr8D9q2uwiFqBHaBnaoMO+XlnSIHgBAkr7JuXCosZoyFWZhLC3GL6o9f/9vJPbiK4uQjjflSmpyjizvXzn6Zuz9cjGdgCD+9NIP/mzmJtIR/P+RLi4tY+sX7PHxtX+67qiMn9u9i5NSGX1iwrpw9fWh37R10vvlRsk8eIj8loVKZ6gRb16DWFR6XB++aUuycUL3bYEndi6U4s9rnCYsJS1Ea5vSDmE6sxHxyFYpLAKpX3Sc5VKifTo+mRR9+X/E7P/+4uF6vXV2X99KQknQFU91DUdqMBUVFsXdFE9AF1T0Cc/oBRM6/v1wV1xAUvRdkxWFJ3YPq36nZTbu/HJQVZNf6w8a143WoDq7k7/mR+PkPYshLRZisg5qDRz+Oa+s+5OxdTvae5egD25KbXlStD8OzqYXVbuVpzvwiorjphU9J2LeN3z98gR9ffIjb3/2eBT9/xfZfvqS4II/eI8Zxza0P0LHvYNy9fZu6ypW08nGu1MoDkHF0Nw5uXrgHR+Kp1XI2tWIZd1+nCi097gEBlVp5zmfv4lmrlh41sCuiOBPzid+xuEegaB1AAcxGhNkAFiMIi7X12mz4p/W4yHpM0aA4+aAJHYDiHt4gv0NUtxDuuf9+5jw5hz279/DSqy/U+z0uRgYeSWrGFAd32/8LYwnCbEAT0heCumNJ2Y0lKw5hKkPj7I8I6Gw9ln0CTVB3VPewJqv35epioSc/6Vilv8bLWUxl6Fv2x1JagMVYjE/0YNI3WhfZO7v8VVs5fYt/N1isbuipjqMp+c22W+tc4R17MO6R1/h8+kTeuL4XqlZHnxFjGXf3LLwD6rY3W1PwLEsjcfPvdLxxBmoNFoo8P/S4BrWu1LVVG4qqRdNyGJbUPVgKU6whRwjQ2KFo7ECjA1RQFdDpUbX2YOeMqvcBR08UpeE6fYShCFGWh15v7W4z/bNBcmOSgUeSLhOW1D1YsuLOO6qAoQBzxiFQNGjCr8KcsAbzqbUoHafJlp5aqGnoKU3aR/bGj0DV4tJuFKI0D1NBJr59bgZVg6pzIPUv6yrXHh2G17g+1W3lOZrybzdncw4/AZHtuO2d70lPiKNl1364evvh7VO3VqzMlLP8vfwnTuzfhZOrO5Edu9Jz2Bic3ep3BlV+Zhond2/GYjaSlnCM/WuW4uwXTMuB/66QHOzvfMlWnuqobSuPotGhCepOc5m+YCnOxHz6byi1zrL8/PPtPDT9Qe66985LnFn/ZOCRpMuEGtgVUZaHKEw956hA5J5E5JbvUfRvwDGfXI225dWNWscrRU1Cj9Y10Po/FhMFB5YAUFLViYC9Z4sKjxuqa6u5h5+gqPYERf3b2nUio7DSAObq+nv5T3w6d2aFY1tW/sqP/3uNPqPHM3rafXj6BtSpvgC71v3BvKenYyi2BhdXnwC6jLoBr17j0OguvZ7Upbq26quVpzmxFKZZxwQ5uKOGDUJx9OTvX59tstmkMvBI0mVC0dihbTUCYSzBfHYzIu80KCpo7MBUijZqLJaCJDCXYUnbD+ayJpntdKWobujRuvgQMOEDik6sJ3//r2D5d6aOc3hXFI0OR79I3NoOQFEb/xd9cw8/5Wobegyl1ngZFNGaAeNuJCn+GOt//Y6SogL+XPQVfy76isfnfY/exZWwtu0vcbXK0s8msm7xd/w2fx5Rva5izMwXsXdysf27Ovf9LVdVK091nB96atvK0xwIUxnmM3+jOHigaTkMRWNdkqApl86QgUeSLjOKzhFt+FWI4kxMpzdBqfUXoilhDdrIUSg6RzQBXZq4lleG6g5kVrR2OLcZir5lXwoOLKPo2BoALBonwq6ZfdFzG3MAc3MPP7UJPYMn3MzgCTdXODZg7CS+evkJTh+zbr756r2TAOh21SjueOYNHJ1dLnpNi8XCwa0bWPZ/HxC3Zzs6e3vG3TmTdtfehnLeIohtAlyrDD1VacwBzE1BCAuWjCNYso7CPzNJ1fAhtrDT1GTgkaTLlKL3Rhs1BpF/BnPCGjAUYDr0PZqwQXLAciOoajyPqnPErfNE9OE9yfjjJXRu1etKqc8BzNXV3MKPsayUncu/Y4+pCI1Wh0artf3XYjZTWlxEfnYmp48fwd7BkbZdezHi5ruxs6+8mGOr9p15YeFKjuzcwhcvPEJG0mlueOhJln7xAc9OHcPdL7xLRHRHW3mzycS2VcvYsvJXMpLPkJWaRFlJMeHRHXjg1XnE9OyP3sW1yllaF1LdVp4roWtLCAsi7wzmtL1QkoPiEYHqH4ui90axb5ztbqpDBh5JuowpioLiFoIScyPmxA2IgiTMp9ZidvBAG1W7vXSkimozc0vnEYL/de+i2jmSm5KCe0Ddx5BAw01TLw8/TRl8ti/9hjVfvAWAp18AZpMJk8mI2WRCo9Fg7+iEs5s7LSLbUlpcxJLP32f76uXcNud1WrbvVOl6iqIQ3a03T3zyAy/dMZ6V337GLY+/yIqv5/Hc1DHc9ezbBIS3Im73Ntb88DXpSadp27U3MT374+UfSMuYTkR27Frv/4aqM4D5cunaEqYyLNnHsGQeBUMhirM/auQoVCefpq5alWTgkaQrgKJ1QNvyaiyFqZhP/A6ludY1NzTV35xTurDahJ5zV1yuTuhpDmvzNGWrj3eLlgA4OLsR1rkfMR3aE962I6FtYtBWsSXF6eNH+GzuLJ6bNobobn2I7T+E6G59CI5ojarRUFyQz+bff0FVNdz61Cv88vFbfP3Kk8x4+/9Y/MlbfPLMDMC62Kizmzvj732YyA5die7ep95e04VaeS7nri1hNiLyz2DJSfhntXeB4h6OGjqw2QadcjLwSNIVRHX2R429tamrcUWqy8KE1dUUXVsX0tjhJ6rnIKa+MZ+9qxaTdHQfe1cvxmIyYu+op1P/Idw0a26FBQlDItvy7PxlbF21lE2//cwP77+CyWjAwcmZiOiOZCSdJis1GUVRMJv/HUj+44ev8dRnP5IYdwhVo2HDkkWs+XE+P897E4APVu3GzatmH9w1GcdzIc29a0sIgSVlF5aMwyDMKHpv1IDOqB4tUXRNu49YdcnAI0mSVE0XCj0XW5Sw3OXQtXUhjRV+Qtt3I7R9NwDMRgOarFPsWvsHv339Ef4hEVx3T8UB4Bqtlk79huAbHIqLuye5Gakc37+bY3u3k52egsVixtMvEJ+gFrh6euPi7kn/a29Aq7OjZYy1GyysTXuiOvcg8ehBwqM71jjsXEp1W3mq0qy6toxFWNIPAKBtcx2KQ/MZm1NdMvBIkiTVg/oKPQ3VtSWEYOu8OZzZ8Sc+UZ0pykyhOCsFRw8fPMPb4RkRjaO7D2aTAQdXTzwj2mEozCNu5Tc4eQWi9/bHM6wtR0SEbVxLQ7f8bN+2nU1LF6F3dqXnsGsRQpB08hiZyWeJ6dmPkqIinpt6DelnEwG47p7ZXHvHdAylt/H8tDGcOXGU7LRkstOSUf+ZDp2bkcakGXPwDwkHrGN9el49hp5Xj2nQ11KVy6lrS7FzRg3ojCVlD6aENWgCu6K4trisxgnKwCNJklQDtd1+oiYaomvLYjKQc9q6UndG3G7b8ZKcDJJy1pG0e121ruPo4UOL7kMJ6z2Co/z7Wusr/AiLhQNrl7Nu/vvkZaTQcchYbnnoUc4cP8x7s+8gJTEegLC27ek6aATpZxPpN2Yibp4+/PLxW6QmnuTOZ9/muW9+4+jubfz04eucOXGEN37dyIr5n7Dquy/ITE3ixYUr66W+5S7WrVWTdXmac9eWxq8jqmsLzEk7rDNDHTxQPSJQvVpb9+1q5hRxpWyvfI7du3fTpUsXfv71J9rFRNfbddtPeKPeriVJ0uXtQqGnOoGnOl1b1Q08NW3lyU2MI+XgVhI2LKHdtXfiGdGOsvxsSnIzyDsbT2HaWVIPbsVY8u8HtMbOng4TH8TZO5CUg1s5s201ZQU5uAZF4BfdDSfvALxaxuDV0rqwX23CjxCCEzs28NeX75CWEEdUr6sYNHU6vmGRbF/yDX98/DIxPfsz7KY7AFj25f+IP7AHi8WMk6s7I2++G72LK1+/+hRj75jB2LtmAHBs7w5evP06dPb2GMvKKtwzsmNXHPRODLvpDjr0Hlitel5savrFxvFcLPCc37VVVSvP+aGnKQcwCyEQBclYso//swCqBm2r4dZNjC/hwE+P1GtdDh08zPixE9i1axedO3e+aFkZeGpABh5Jks51OYaemigryCUjbjenNq0gee9GnHyC6DjxAQJj+5N6cCtnd/5F5vF9FGenYzEZCO93DaG9huMT1dm2QF91wk/ayThWffoqCXu3EtK+K1fdOosW0dYxNid2bmThnLvoOf5W7nn8WdRzFv7LSk1m5qgeePj6U5iXg3eAtYsl7cwpPt14BJ2dPRaLhf2b1pJ6JoF1vywkOeF4pftfe8d0xt/7cLXek9oGHqhb6KmqladZzNoylWKOX4UwlaBtPeaSA5ibMvDILi1JkqRaqssg5suBvYs7wV0HE9x1MLlnjrP/xw/Z/OETeEd2pMPEB+l+xzOAtRvq6Ir5nNywlISNy3DyCaTVVROJ6H8NR8/53D4//BjLStnw7Uds/ukLPANCuGHuh7TuOajCuJB18z8gtH03ht7xSIWwA9YZVwBB4ZHcNHsub8+4lazUJFrHdmP7n7/RZeAwHPROxPa7CoChE6eRdPIYJqMBs9nMiX07cffxo8fQaxri7auTmnZtWYozEIVpKHYu1gHFGjtQNNad0utpnI0QFkT2CSw5JxFlBWAuA1VrvZexGEv6fjRBPerlXg1BBh5JkqQGcKnQc7mszVPOvUUk/We9S+rBrexb9AF/vXQH7iGt8Y/pid7LDxf/ELrf8TSF6Ukk7VrH/h8/4NCvnxLQsQ9ugRFoHRw5obNHmM142lvITUsibvMaigvyGHjLQ/Qef2ulTTjz0pNJPnaACU++g6IolbaeGHTdTaQmxnNw20b2b1rLU5/9xBsPTuHY3h3E7d6Gg96J3iOvY/w9D+Pi4YlGqyWk9b+t/pEd6ncLlrpMT6/LjuqWrOOYz2yy7q0nzBULqToUe1dwcLeufKz3QnH0QlG1CIvZumaXqgWd3rYFhCXvDJbMIwhjEYqiATsnQEEUZ4GxEMUlCNUjHLT2YDYhSnOt5zv51eq1NxYZeCRJkuqgLuvz1HfogYbr3irnH9MT3+hupOzbxOltq0jc/Dul+dnWD8/zmMzFnNm2mjMXuV6/m+6l13XTqtxxPOWEdS+sFu3+7ao4N/REderBs/OXs+j9V1j0wSss/eIDHJ1dePmHP7F30LNx2Q+s+u4LTuzfxQ0PPUmLVm1wdvesciHDxnCpwcu13VHdkn0cEGjbTQRh+bf1xWJGGAoRZXlQmoMlNwGEBVBQHL0QJZkVL6TqrKHJXIbi5IfqHGD9vhqLAIHq1gLFo2WzX2DwQmTgkSRJqqPm1LXVGMFHVTUEdepPUKf+gLVLy1BcQFl+NqX52ZTl51Can43ZWGZ73mwow2woxWwow2QoQS3OJTftLBsXzsNsMjLktsqbrKadjMPR1R1nz4ofsOe39Ey47xHMJiOlxUWMv/dh21o64+6aSWy/q3hnxq288cAUW3kPX38CwyNp170vV0+6FTuHy2PhPLCGnryzR7EUpGIuSEUUpCKK0gAQhamo7mEoOn2V55a36FiKM6znnBd4FPdQFHtXFHs3FLfQy2rKeXXIwCNJklQPaht66rOV51zntiQ0dKuPoqrYO7th7+yGa2B4tc9zyT/FVw/fTGlBXpXPp508il9Emyo/eM8NPRqtlptmPVPlNcLbduDd33eQfjaR5ITjFORmk5F0mjPHj7D4k7fZtOJn7njmTSLaxTbKB3xdW3lKzuymeOd3iLJ8UDSozr7oAmKxOHihuLa46L0VVQN6LzR6L/Bug2jRF5FzAvOZzQCI7BNo2k5Asb/4bvKXKxl4JEmSGlhThZ5yjRl+qiv9yE5+fnc2AZHtGHzrzErPmwxlJOzbRvdrp1Rxds2oqop/SLhtscFyZ08c5cMn7+e5qWNw8/Ihtu9V9L/2Blp16FLr8FMf20ycrzz0GHNOk7PpUxyCOqB4t0V19reGmBo4d2aXompQvKIwp+4Ho/VnRBiLKgQeISxgKgWt42Xf4iMDjyRJUj1pjEUJ66qxxvpciBCC+LU/s+/79/Fu3ZHuD72B3tWjUrndK3/CUFJExyHXXvBa53dt1VRwqza88O3vHN29jUPb/mbb6mWsX/I9AaEt6dBnIH4twvEJCiG4ZRRe/oG1vk+l+9awladc4ZFVqI5uePS5G0XV1GpBwqp+PoVXBKbU/Wg8wnAMsA7qtpTmUXLwF0SJNSBp/WKwbz2sWUyFry0ZeCRJkupRbQcxN3Qrz/maotWnNC+LHf/3Iin7N9Ny8Hg63vAQWjsHjqbkV5iybjYa2PTDZ8QMHIVnYOhFr1nX0KPV2RHTox8xPfpx/QOPcXTnFjYu/4k9G/4kKyXJtvFoSFQ7rr39ITw69K/1vWri/NDj7O5E8ukduMaOr3GrzqU4RA5FhPUFzTmDuS1mW9hRHNwxpR2kTKNDcXBH0dih9YwAnSMgQAhEaa61ZUi1XqM5tgbJwCNJklTPmtt4nktpjFaf5L0b2fF/L4ECfWe8RWDHvhWePzf0nD60m4LMNHqOm1qta9c19JRTVZXo7n2I7t4HALPJRE5GKicP7WP9r9/xwaN303nEREY+8DSq5sIfn9Xp1qrJdhMAGnsndK5+FMVvxCG4E1pnnzpvOyGEwJi0C0tJDpiNaFwD0biHoDi4oTh6oO9xDyX7FyFKctC4hWDKPIYwlYHFhNE1GNXeBVNmHKgaMBtt11XsXdF3u6PZhR4ZeCRJkhpAc5q5VV0N0epjMZvY8+1bxK/9hcDYfnS99UkcXC/eAlZWVACAi3fTruui0WrxDgjGOyCYbleNZMPSRfzfi49jLC1h3GOvN/j9z23lUVQNoeNfIOHHp8j4/Xn0rQagD++F3tWdorxsFEW9xNWsLKX5WArTEGYDluIsjGd3WKeZq1pMGUeBfzZfUNR/prBb2bcajPrP1hGGpN0YTq7FAmg8wtG4BaM6+1F2fBWiLB+Ni3+zCzsgA48kSVKju1joaapWnvPVV6vP/h8/5OSGJXS+5TFaDhx30Q/C8lYe33Dre3N8+3o6DRtfrfvUVyvPhSiKwoBrJ5FVbObXNx8ndvh4wjvWbVXh6rTynBt67Nz9aXXLBySt+5qi42spilttLaSxQ+vZElXvASgodnpUJx8UnRPmvDOY886CMGPOT0aUVpwRp7oG4thhEoqiIIylmAvTrOv4WEygtUN1cEfVe6Fo7W3n6PyisRSkoNg7Yxfa19bFZnTywVyWjzk/GWPqQbR+7ZpV8JGBR5IkqYHUdTyPMBspTolDH9i20riNhg49ptIico5s5nRWMoqqIXjgTbQI/HecjbBYbPtlnc9iNpF+ZCcJG5dxZvufdJj4IK0GXVet+x5NyadNYCjtBoxg3fz36Tjk2ot2H52rJqFHCEFmyllOHTlAbmY6jk7O9Bw2Bm0VCyCea9ykm9jyy1ds+PajOgee2tDY6wkZdi/ZbYZhzD0LFhNlGccpTtyJOTcREAhjSYVzyltwtJ4RaNxaoLoGWnc3FxZQtbZQougc0HpcfMwUgKJ1wKHNqErHHduNxVyUgfH0NsqO/4Ex7SAOkUNtLUNNTQYeSZKkBlSXrq2S1OMk/vgkAH79b8UptDP2nsH1Pmi1Kqd++4j0nSttj8tyUsjrMIiyvAzMpUWc+m0egZ36Ez3mNpx9gjCWFpN98hAp+zeRsm8zZQU5OPsG02Xq40T0v/BMq3MJIShIPc3S1VtJPnaQkoI8TIYy7Byr/1FVndCz86/f+f69l0g/mwiAqtFgMZtJOnmc6+9/FFVT9ftbUljA588/QtrJowD89dW7ZCclggI+Ia3Q2Tuw87fv0bu4033szeja9L9gKCxX01YeAGNhFoaseCylBVjKCnFq2Q/X9mNszwuzEWNuEhZDMVoXH7TOjbMycn7SMTROPmjajsaU056yE39SvHs+Wp+26AI7onG59Ia5DUkGHkmSpAZWm9CTm5KCW0CU7XHahi+BLwHQOnniGtmb/OAYAmK7o3Nyq9f6Zh/dWiHsOHgGkn2k4jGA5D0bSN6zocIx16AIwvtdQ4tuV+EeGlWtLo20wzs5s30VqQe3UZyViqrV4RfdjTGzX8bOseatWBcLPYayUj6e8xCRsd2YPHsuLWM64aB34r7BHfjt649Y+c2nuPv44eUfSJeBwxg+5S7bazi4bSM71vxmu9b+NUvxDAxBCMGpfQspzssmssdAhMXCr288hp2zG6PfWorWzqHGr+F85aHHYizlxFf3IoyltufMJTm4dZpoe6xodNh5hdX5njVV4Wc5qDWizQCK4jdQeHQ1JXsP4d5jGgAGg4GUpBRCwkIatctLBh5JkqQmdKnQY+8VQlnW6QrHTUXZZO9dTvbe5Zxdbj3m23U4rmHtsfcMBGEhc/96HL2DCexbvTEw59L7hdn+3ykwkvDR9+ES0o6SrLPYu/lQlJpA7rEd6H1DsPfwpyw3HY29noh2bXHyrv5f8RaTid3fvMHJ9b/i4h9CUKf++Mf0xKdNZ7T2joSet7t6XZWVlPDzvDcwlJUyedYzBLdqY3vu7d+2knBoH1mpSWSlJZN88jjfvfsioW3aE92tNwBdB4/graWbOJVdjLOHNzr7ikHGbDKi0VqnZX9w69XkpJzhl7sHENR5IG1G3oxXy5gq61WTGVuK1h6f7hPJP74ZB59QhKM/+vBetXk7GpyitcM5aghOkYPI3foleTu/ZeZDafy+4ncAPvnsEwYM6td49RFCiEa7WyPZvXs3Xbp04edff6JdTPSlT6im9hPeqLdrSZL033Oh8TwX69py9fYg58AfFCcdoiT1BL69J6NxdKUk5Sj5xzZhyE2u1r2dW7TFpUVbXCM64hreAZ2+YpgoTk8EFPS+IZRmJXPoi0cx5GfS5bGF2NVwHFJ1BjpbLGa2fPQkyXs30nnKI0QMGAtCkJcUj9loQNVoMZUV40kRBdkZFGSloygKPa+bhrOHd7XqUd7KYygrZcOSRSz94n0KcrK5/oHHGHnLPZeon4Xpw7vh4etH+579cfX0wT80gtDW0WRy6Vangqx0Vv3wDYeWfAZAcNdB9L7/1QuWr27gOX9BwnO3nbAYSylLPYzWxQ+tW2CzGjAsTAbydn9PtGcZO3fuBGDw4MEoioLJbKJVq5ZcN34cLSNb1ui6hw4eZvzYCezatYvOnTtftKwMPDUgA48kSXVVm9BzsVlbwmKmLOs0ObsWkXNkc7Xr4egXhpN/S3RObpRmJ5NzdCsAqp0DFkMpdm4+RN30DC4hdfsdeqHwk3l8P3+9fCfuIVGE9hpGUUYySXvWU5KTUamszt4RFy8finKzCYvtwQ3P/O+S97WYTWxf+i1n9/xN4tEDFBfk03vkdYy7aya+wZcemAuwZeWvrJj/CUX5ueRlZWA0WDdDdfUJICiqPUFRHQhq04GgqA5o7ext9z20/ne2LfmG5Lj92Lt40HLQdbS+ehJ2ThdvsapL6Mnb+zPFJ9Zb18kBNHov7ANjcIochM6tacfOnOvswjt5/PHH2bp1K05OTuj1elRVZdu2bZw+fZoxY66hR6+ejL9+XLWuJwOPDDySJDVj9R16AE4unEVp2gkAvNoPwL/HNbiEtQdhIfvIFs7+tYDi1IRK59m5+eLgFYBndB8cvIIoST+NvZs3Hm17o6mHsSflzg8+xpIitsx7ipzEo5gNZTi4ehLQoTeBnQZg7+KGxWxCa+eAo4cPWgcn2ga6sWPZQlbOe5mZ36yttIP6uYQQLH79UQ6tX0HrHoNo0y6G3iOvIyA0otb1t1gsZKUmkXB4P7t2bCcpbj8pxw5hLCtBZ+9ISPuueAaGcGLHRnJSTtOyS186DZ+ACO2C5hIzv8pVFXgKTh8hc/9aynJSMZUUoGh0mC0a7NwD0AdE4RzWhfysXDLXvosh7Qiqgyt2XuGgajFkxmMpySVgwgco2urVoaElfXeX7f8NBgMbN26kpKQEIQQ33HADJSUl2NnZse/Qnmq1UNUk8MgxPJIkSY2sIRYlDL/hdYTFjKK1Q1EU3M6Zsu7dfgBeMf0pPHOU3GM7KM1KoigtgeLUk5hLC3EOak1Ar7HW2V9tetb6dV3M+Ysa6hyd6D/r3WqffzQln5iBo1j16Wvs/2sZvSfcVmW5otxs1n79LgfXLmf8k2/Trv+IeluF2SewBT6BLfDsOACwtuakJRzj5O5NJB7YScKeLQS16cCEp94hoFW0rd51cfzHVynNPIujbwhOgZEIswnVUErBia1k716Czs2f0PHP0/LGVylJPkLWrl8pOLkdhAVV5wiAHXk4BXSo2xvQAL755htuv/1222MnJyduve1Wbr19aoN0x8nAI0mSdBm41IKEikaLcs56NeXdHuVr9SiKgktIW1xC2trKlOWmkbL5V5I3/oBrRCyeDRR26ktioUKb3kPY8/uP9Bw3FVWjQQhBbupZEg/sJG7rGo5v34BGq2P09Odp138E0HCLEqoaLQGtogloFU2fiXdWWaYmu6dXNXi55dgZHFv0CqVZKbhHdiOg9zgcPAPISSugNO0Ep354ktOLn6PVtHnog6LRB0VjKs6l6PQ+DLkpaPXu6IPa1fm1NoRJkyYRHx/Pyy+/DMC8Tz6ke8/uDXY/GXgkSZKaQG2nql+qa6vSOeeM9zh/oUJ7dz9CR9xFypbFlGWnnH9qgzmbWljrFZx7XjeNL6ZP5Jun7kBYLKSfOkZJfi4oCkFR7Rl6x8O0H3xNpR3YG3ol5obi1rITnWd/TfLGH0nZ/AspmxfjHtkFjZ0j+WeOIcwG7Nwr7uSu1bvj1mZAE9W4+vR6PS+99BL33nsvU6dO5ZYp0xg2fBiTp9xIVJso9uzey3vvvM+sh2fSb0DfS1/wEmTgkSRJaiKNvd/Wua0+FqMBQ0EmJelnQAigcWf01Db0FLiGMuL+pzm2dS32rs50HzOFwNYxBEV1QO/mcdFzqxN6LBYLuZlpWMwWFEUhOeE4aWcSAIWwtu1p1f7i40TqqqpWHo29Iy2G3EJg/+vJ2PMnWQc3YCopxKttd3S+HXAO79KgdWpowcHBrF69mu+++47XXnuNW6ZMq/D83xv+loFHkiTpclfT0HOxVh6zoRhzcS4lqccxFedi5xGMzsULjYMLGgdXFFWl6Mx+Ti9eSmHiHtvmkE5BUfh0Hlq/L6waaht6XLqOZvI1N9XqnuWhpyAnm0+fncXJQ3sJDGtJi8i2FORmc3j7JgrzcxGWfzfO1Gh1IARms4n7Xv4f3p2vqtW960pj54h/j2vw73GN7dj5M7YuV6qqMnnyZG666Sb27dvHe++9x6+//kpubi5ffz2fdevW06VrFzrGdmTQ4AH4+vnW+B4y8EiSJF1mqgo9ZTlJxM9/ACxmABSNHcJsqHjiPztgO/hG4D/wTuw8gtA6uWPvFYLWoWE3I72Q2oae8o1GayM7LYXX759MQW42V024meSEExzdtRV7vRODxk8mNCoGe0c9FrOZgLAIfIJCyU5LZtboXuRnZ9HTx5kTGdWbQg41G8cDNVuI8PxtJxqCsSCDorMHcfCJqPbWJrXd501RFGJjY1m3bh25ubm244mJiZw9e5ZfF//K3KctdO3alXbtosnKzq72tWXgkSRJamL10bWldXRDY6fHXFqAaqfHqUUH9EHROAa0wVxagLm0AGEyYO8ThqN/5S0fzh/kfDmoTejJOpvAe0/ejlZVmfPFL5ecqi6EYO/GNXzz5ly8/IPoM7rmK1c3tJqGnpp+jw8tfZa8+N2VjrcYMpUWV91So2tV19GjR/n555/566+/+OKLL4iKimLU6FGsXrWKAwcOkp2dzYEDBykqqv7rloFHkiSpGahL11bOwVVkbF6I5Z9F5yyGYgrit1KakUBkl7E1qsfFBjk3lLoMYq5u6BEWC7tX/sRfX72Dk7sXU176/JJh5+C2jXz/3kucjjtE69hu3PPi+zi51O++ZRdyqVYei9lE3vFdZB3cQGHycUozzqB1csPezQet3hWtkxt631BcQmNwDoxEreZaQFWJnPgYKZt/JW37ckwlBbbjucd3EjTwJgx56di5+VZ7V/vqsLe356abbuKmm24iNDSUZ555hlOnThEWFgbAyZMn+frrr4mJiaFLl+qNYZKBR5IkqZm4VOgRQmDMTkTROaB18SX71AmKDy0m/9jfOAa0xj30auw9gjDkJGPITcG9Xd3GmjRmq09dQk917PnjZ357fy7R/YYz8sFn0Lt6XHAQc1LCcVZ8PY+Ny34kqlN3nvjkB9p06VnntWGEECTv2UDGsb1odHYYivLR6OyJGj4ZR4/q7Whemp1MyuZfydj7J6aiPBx9Q3ANa49vp6EYi/MxFGRhKs6nJP00WfvXYzGWgqri6BOCU0BLnAJa4RRo/W91N521c/UmdPgdhA6/AwCLyUBu/B6yDqxnx4vjMZcWonFwxrfrcPy7j8bRp0Wt36OqPP3001xzzTVs2LCBhIQE+vTpg5+fH6NHj+bUqVPVvo4MPJIkSc3IxUKPk5c/matfqfK80PEvoOrqb2XkczVWq09DjufxDApFo9ORkXgci9mMEIKCzDS2JxvoHB2JqtGwf9NaVn33BQe3bcTNy4epj7/EoPFTUFW1wrVqMn6nnLGkkK2fPEPKvk04+QQhLGa0DnqKMpKJX/cLnSY/TET/MbbyVbXu5J86yKEvHkajc8C363B8YoegD2h5wSAmzGaKUuMpTDpGUXI8xSknyD60yRqCsAYZ98iuuEd1J//kPopS4tE5u6P3j8DBww/XiFgcPPwRFgtFqfHkHd9FwenD5CcexFSUh72HPwG9x+ISEk3eyb2k7/idlL9/wjm4DUEDbsCzXb96W0AwNjaW2NjYSsdrEnjk1hI1ILeWkCSpsVxo+wk7pYzsDZX3krLzjaLV5Mb9HdVQ4acuLT0XCz7pp47zxfQb0Oh02Dk6kZ9hXXtIZ++Io5MT+dmZhEd34Oobb6f7kFHo/tkf61y1CTv745NY9+q9FGWm0OOu5wiM/XeKtamslJ1fvczpbasY/dZS9B6+F+zKOvHLW+Sf3EfHhz6t9bYfwmKmNCuZopR4jn33gu24nbsvrmHtMRbkUJx2EmNhLgD2ngG2NZpUOwdcWrTFJSQaj6geOLdoi3JOGDQby8iN20bq1mXkxe/GKTCSgN7jsPfwx97dBwfPQDY9PrhW9b6Q8s97ubWEJEnSZaqswDr7pDz4CGMJxvQjGMS/06W9r34CRdGQu/NbdO5Btp2za7o4YW01VMtPQ43p8Q2LpMe4W0jYu4WQdl0JiemMzt6RtIRjlBbmMXDoMFp16HLBVonahJ2jKfkcXbGAosxkBj/5GW7BFXcD19o74BHWlrM7/0Kjs7/ouB1Vo0VRNXXa40xRNTj6tMDRpwVmQwln/pxPqwkP49ayc4XXbSotJPvwFhJXfgqAc4u2xNz1DqpWd8Fra3T2eMX0xyumP7kndpG88UdO/PT6P5VX6XDfh0D9Bp6akIFHkiSpGSvNPIVI3Y0559+NPxWdHl1IT+w8wwDwGfp4hXPKgw9cvuGn/IO/vru4Bk+bAcyocCyic2/b/9d32AHIO3sCn6jOlcIOQO6ZExxb+S1+7bqTUXjx6d6uEbGkbl1KaVYyDl6BFy1bHX5dR+DXdUSVz2kdnPHtPBTfzkMpyTyLztnjomHnfO6tuuDeqgtluWkUnDnKsYXPE7/4Hb7t6cKQIUPw8/Orc/1rSr10EUmSJKnJWEwVwg6AMBZjPL2F/KRj5Ccdu+jpuSkpFQJQY8hNL7J91VV116M5X2037awq2NQl7AAoikpJbiaG4oIKZfKSTrLutXuxd3EnaORDl7ymR1R3tHpXTq38FPHPekuNwdE7uNbrNOmcPDi28HkAipKOMWXKFIKDg5k6dSo5OTn1Wc1LkoFHkiSpOXNwR/X+d8NP+7bXogvqin3USNuxS4Ue+Df4XI7hpylDT13DDkBwt8HkJsax5MFh7Pi/lygrzMNYWsyWeU/h4O5N61vfxM7V+5LX1dg50nLcTLIPb2L/h/eRfWQzzX0YrqLRENj/BoIG3Ej3uUvYumMzT899mmXLlhETE8PHH3/Mu+++S1xcXIPXRXZpSZIkNUPCYsaScQhLxiEwlYK9GxqfdpjtPFB8PNGeN6i5PPRUZ6HCxh7rY7tvHaa5N/aKzLUJOuX3O19Y75H4RHXm7I41HF72JUm716HTu1BWkEPMPR+gdaz+6/KK6U/MXe9yetUXHJ3/NA7ewQT1n4hv15H1NiOqPimqhrARd9keu3u44+zkRE5ODjk5Odx7770AzJw5k169evHXX3/h4NAwsw1lC48kSVIzIwxFmOKWYEnZjeoWhrb1GLRtxqF6/7tCcvmg5vNVp5urXFO2+tRGY7f01Od9nLz8iRo+meEvfkdY39H4RHWi3V3vovcLq/F9XEPb0e6Ot3CP7Epp5lnif3kbYTbWoeaNy9PLGtZbt27N62+8jqurNZAmJiY2aGiTLTySJEnNiDCVYk7eAWV5aFr0RvWKumDZC63ZAzXfcb2xBzrnphddFi09Nbl+dTh6+BA7aXqtw1u5gtOHyT2+E9fIPvgNvANVW/uVlBtbrz49OXriMABxR+LIz7e+d+Hh4djbV14KoL7IFh5JkqRmxJJ+AJGbAIoGxfHS4zrKCrLrpbXnXI3V6nOltPTU9Lp1DTtgnQIO4Nl5DDpnr8t21/SotlF88eXnAGzatIlly5ZhOWen+vokW3gkSZKaE8X6d6g25kYUTfWnAV+qtQeqN77nXI3R6lPbcT11aemBiy9QWJvrVVd9hB0hBEnb1wCg2jnajte21ayp9enXm1V/ruS5uS8wZswYHBwccHR0xN3dnbZt29KzZ0969+5NVFQUer2eM2fOkJeXR4cOHWp0Hxl4JEmSmglLQTKW9IOoPtE1CjvlLhZ6oPbBBxp+oHNtPqwbY9PRS12jJuoj7JhKCon78V3yjqzFp88UHLzD6nzN5iAkLITPv/qUuCPHOHjgIGVlZeTm5HL0aBxvv/02zzzzTKVz7O3tef/996t9Dxl4JEmSmgFLQQrmk3+iOAegBnSt9XXOX6G5KjUd33Ouhmz1qW3ogfpfoLA659ZEXcKOxWzCVJRHwZkjxC18DiwWAofPxLVVLwy5qejc/GyDfS/XVh6wLvrYJjqKNtEVx61ZLBbOJJ4hLS2N4pJSfH190OsduWXyNF55peq95aoiA48kSVITE6W5mBP+QnHyQxM+GEW9+Iq71dGQrT3lGqLV53IYzNxYYSdp4w8krvjE9ljRaOGf8S1p6z4neeU7AASPehTX1n2rvMaVQFVVQsNDCQ0PrXD8ldde5rVXXq/+deq7YpIkSVL1CWHBlLgBdI5owgehqPX3d+iFBjOfq7YDm89V39Pbm/Ng5sYKOxaTkTOrv/r3gKLiN+BOWox5iqCRj6LVuwOgD47BOaJbhXMv1wHMNdWnX29eef3lapeXLTySJElNQJgNmE+tQ5RkgakMTeQoFE39Ty2uThcX1E+LD9Rfl1dzbOlpzG4sVauj21M/g6KgsXOoFGIsxhJSVn9A8dmDxH18M63v/BKN/b/v1+XctdVQZAuPJElSExAlOYiCJDCVookYgurk06D3q05rD1Rvm4rqqmurT3Nq6WmSqef2jlWGHQCPmKE4hXUGwCWie60Guf/XyMAjSZLUBFRnP1TvNqDTo7oGN8o9axJ6GiL41Cb8NIfQ0xRhp9zFXr9TcHsA3NoMqHLhwf9K11Z1ycAjSZLUBIQQWPKTUFwCG/W+F1uo8Hz1HXyAyy70NNewA+DVZSzOYV04s+QFSlKr/j7J0PMvOYZHkiSpCYjiTDAUoHr0bpL7X2oW17nqa3xPudyUlBqP72mKMT21uVd9qU5QUVQN5jLrPRO+e5jAYTNxjx5Ub3VoTiwWCytXrMJkNOHs4oQQAiEEpxMTq30NGXgkSZKagCg4Cxo7FGf/JqtDTUIP1G39nvPVNvRA463KXNN7XIrFZKAsNwMHD38UzYWXHqhJq4xX57Fkbv+B0owEMjZ/U2XguRIGMOfn5TNrxqw6XUMGHkmSpKZgKgOdE4rStCMLqjuLq1x9tvbUJvRA4y9QWB8KTh/mwLwHAXCP6kH0tKqnU9e0C8q1dR+cwzpz9MMbcA7rcsFyl3vocfdwZ8Ufv7Fvz34KCvIRQlBSUsLpxNP88svial1DBh5JkqSmoKggGmaTxNqoTWsP1N809ubexXWpa15MUWoCBz/9t3XCr9vIKsvVJOwIYSH34GpMxbnk7F0BgM696VoLG0NEy3AiWoZXOHbo4GEZeCRJkpo1RWlWgQdqHnqgfoPP5Rh6qtOVte+9OwAIG3UvAb3GVdmdVdOWneKkw6T8+WGFY04tOl70nMu9laeu5CwtSZKkpqCoIERT16KSmsziOld9zOa6HGZw1eYa7e54k9gZXxDYd0K9hB1TST5Jv78FQOSdX+HVZRwAwlRWo+v818jAI0mS1BQsFlCb76/g2oae+timosbnNEHoqcm5bi07ofcLq/K5WtVdCEyFWQDonD3x638r0TOXog+KvuSp/+Vp6s33X5skSdIVTJjLUNT630qiPtUm9EDdg09zDz31Nf28pnUWQpC973dOLngIAJ9eNzXKfa8UcgyPJElSI7PkJSKyj6N6RTV1VS6pprO4zlWX8T3NdUxP04UdCyl/fkjuwdW4RQ/Gq8tYHLzD6qUu/xWyhUeSJKmRmRP+AsCSfwZhMTdxbaqntq09UPvxPc2tpaepwo6pJJ+zy14l9+CfBA6bSdCwGXUOO//FVh4ZeCRJkhqZtv1kNGGDwViM+cTvWHLim7pK1VLX0FOb4NNcQk99rqJcE8b8DE5+M4PipEO0GPNkva6kfLmHHiEEBQXV/77ILi1JkqRGpmjsUNxDEUHdETkJmBM3YMlPQvVshSjNQ+SfAYsJYTGhaOxA54iisbc+1jmi+rRD0To0Sd3r0sUFtVutubFXZYaKCxQ29pYR58ravQRhMRIx5T10Lt71Vo8rwXVjJnDkyJFql5ctPJIkSU1E49MObevRqMG9EDnxmOP/wJK0FRBg54yq9watPRiKsBSmQGkOlvRDmJN3NnXVG72Lq6l2W2/KsANQlpmIztkbYTbVWz3OdTm38gwY2L9G5WXgkSRJamIa7zaoQT0AUNxC0LYchja0P5oWvdGGDUIbORJdm3FoW1+D6tUaUZiCaAZr+MhxPQ1fB9c2/SnNSODEl3eT9PtbGPJS661O5S7X0DNj9nR+/vWnapeXgUeSJKkZEEVpAKhuoRctp7iHgaEQ0UzG/dR2oUK4PMb11Ie63NsjZihR93xDwFX3UnTmACe+uo+UPz+iJPUYopmt1N3cyTE8kiRJzYDIPQWA4hJ40XKqsz8WjwjMSdtQXAJRdPpGqN2l1WZbinKNOa6nsbdWqI+gpXFwxqPDcNzaDiJ7z1Ky9/5GzoGVaPUeuET2wrPjKOy9WtS5nlf6thOyhUeSJKkZ0IRfhSZiaLUCjCaoBygq5qTtjVCz6muKcT01PqcRW3rq+16qzh7v7tcTeccXhE54Cdc2/Sk4voWT304ndcP/YSrJr9P1L9eureqSgUeSJKkZUN1CUF2Dq1VW0TqgCeiKyE3AUpzZwDWrmbp2cdVUcw09DXkPRdXg1KI9/gNup9Xtn+PdfSK5+//gxP/dRcb2Hy6btZ0amww8kiRJlyHFsyXYuWBJP9jUVanSfzn0NGZLiarV4dNzEq1u+xT3dkPI2LyQpJXv1Dr0XMmtPDLwSJIkXYYURUX1iUbknkKUFTR1dar0Xww9TRUYtHo3/AfeQfDIR8g/vomkFW822FT2y5UMPJIkSZcp1TMSdHrMZzY1i2nqValtF1djhp76CinNoXXEtXUfgkc9Rn78NlLXflqrazSH19EQZOCRJEm6TCkaHZqQPojCFCyZ1V9xtik059ADdf+Qr8n5eUfXE/fJLZSk1n5H+YtxbdUT/4F3kHNgJSWpx2t1jSsx9MjAI0mSdBlTXYJQvdtgSd7ZbLu2yjVm6GnMLq6anpf0+1uYi3MpzThVq/tVh0f7Ydh7hZK67jO5Xs8/5Do8kiRJlzk1oCuW3EQsGYfQBPds6upc1Lmhp7rr9pwbemqyXk956KnJej0X2oOrPls8vLqMozTjJO7Rg+vtmudTVA3+g+8m8ccnyTu6Hve2Nd909NzXfCWs0SMDjyRJ0uVO1YLWEWEsrvJpIYR1deaSLOtXcRaiLB9MpaAooNOj6JxQnP1R3ULA3g1FURq82lW1+FwqBFXV4nOpEFRVa8+lQlBDdun49b+1wa59LqfgGJwjupO145daBZ5zXer9uBwCkQw8kiRJlzlL2j4ozUb1boMwlYLFjChK+zfclGSB2WAtrHVE0XuhuoeC1gGEAGMxwlCAJW0flpRd4OCO6tES1T0M7FwaJfyUa84hqL7UdpxRTZS/Fve2gzj722sYC7PROdduJezquFggai5hSAYeSZKky5QwlmCK/wNKcwAwn90MZzf/W0DnjKL3RPVph6L3QnH0uuhKzsJiQhSkYMk5iSV1rzX86PQoTn4oDu4o9q4o9q7WFiCNrqFfnk1ThqDLVflrMWENG9nH92Lv39b2fGOFO2g+A6Bl4JEkSbpMmc9usYUdNaCLNaBoHdEE97SGFJ1jja6nqFoUtxaobi0QZiOiMNX6VZyOpTDF2gVmLQkOHij2Lih2LiguASguQVdkS1B9K0uLoyRxG87RI9DorfVVVE2D3U/j5A2qFmN+SoXAc7Fw15hhqDHJwCNJknSZ0gT3BL8O4OiFoiho/DrU27UVjQ7FrQW4/bsppTCVWcf+lGZbt7QwFGLJTYCMgyh6b9SgHqhOvvVWh5q6HEJQYdxqypIPoPMIoSwtDkN6HN5Dn0Dr4tMg91NUDaqdE8JQUu1zLtXSdbkGIhl4JEmSLlOKTg+NuFu6orVH0fqAkw+qVxRgHRAtitIwn92K+fhvWFyCUH3bozj7N2qLz4XUVwiqL4prGDqNEyYHX0rPfg8IcuM2oOoc0HiGo6i1/1i+UFBTVBUh6m9/rcu1dUgGHkmSJKnWFEWxhpuoaxG5pzCn7cUcvxLs3VA9IlA9I1Hsmseg1XK1CUH1RevTGq2PNZjYtRyEMfUAxtR9iMI0ULU4tB2DxiPMWliYaxSALhTUzIYSFK1dXateLc05DMnAI0mSJNWZoigoHuEo7mGIwlQs2cewpB/Akn4ANaALqnfbZtHicyGNFYIshiLMmccxZcdjzjlV4TlV70XpoV/QeLZEGIuxFGWg7zwV1dG91vcTFjOYyjAUFV4wEDXWOCZrGGrVKPeqigw8kiRJUr1RFAXFJQDVJQBhNmBJ2YUlaRuiKA1NSP8GHaBb32q7+en5yoNT6Yk/MaXsq7KMLqgrdsHdKItfgynz32BizjtTp8BjKcoABKqjxwXLXKwLr7EHdTckGXgkSZKkBqFo7NAE90JxDsScuB7zyVVowq9C0TRO90pzUVaQjbCYK4YdO2cwFNoemgpSMW77+LwzFcyqHkstVqcuZ0w9iKJzQnULqk3Vr6gwJAOPJEmS1KCsixwOw5zwJ6Zjy1F9olGd/MDBDUX5b2zpqKgaNK2vQRSlYUk7WCHsKE5+iPyz/xbWOoKqQdOiD8p5LTPVbXWyd/HEmBGHKXU/dmF9GuR9vtzCkAw8kiRJUoNTnf1QIkdhPrsNy9mtWBCgaMDe1bqooZMPqnMgOLg367E+daEoKuak7ZWOi6I02/+rvu3RBHat032EsFAc9weW9IMoHhFY3FrWag+zurhwGOrX4Pe+EBl4JEmSpEahOLijbTXMuqhhcSaiNAdKcxGluViSd2IRFuvWF84BqC7+KE5+jbavV2MQZmOVxxXnAERZHpgNqL4xdbuHqRRzwlpEURpqYFdUn5hK79/FWokaa7ZaU5CBR5IkSWpUikaH4hIALv9OUxYWE6IoHVGQjChMwXwmARCgsUOxd0PR+6C4BaM4B1y23WCqsx9Kh5tBWLBkHEZx8EBx9LSO50GAEHUa1C0MRZhOrgJTKZpWI1Cd/Wp8jQuFoSshCF3RgeePFX/QIiQYV1fXpq6KJEmSdBGKqkVxCQSXQACE2YAozkAUZSLK8rDknYbMw6BzQuPfCcWz1WXZ8lO+ro7GP/b8Z6CWL0cICyL7OOaknaDRom01AsXBvS7VrORKaBW6ogPPp59+xqeffkZoWCgLF32Ll9fl8U2RJEn6r1M0diguQeBinV0khECUZGJJP4T5zN8o+WfQ+HVA0Xs3cU2blqUg2TouqDQHxaMlmqAeKFr7Rq3D5RKGrujAUy7xVCKO9g5NXQ1JkiSplhRFQdH7oIYNxJITijl5O6Zjy1A8WqJ6RgICResIOkdQVDAbEGUFYC4DIRCGAut1HL2srR9a+zpt49DUhKEI89mtiPzTKE6+qJGjUZ0aZj+uuqivtYzqw+X73a6mESNH8OIrL6B3arz9ZiRJkqSGo3qEo7i1wLR/ASInHnNOfDVO0gIqWAz/HtM6oti7gKoDswG0Dqg+7UBRUOycUeycG+w11IYQFkR+EpaceEReImjs0IQNRHELuyy79xrbFR14goKDePOd19FoLp+VPSVJkqRLM59Yaf0feze0YYNA1SJMxWAsASFAo7OGGY0DKIo11AAY8hGl+WAuRZQVWlt+zEZwcLcOlo5fabuH4uSH6h+L+s+4oqYijCVYchOwpB8CYyE4uKMGdEb1ivrPLeJYF1d04HnzbRl2JEmSrkSKsz/Yu6IJ6Wdr3VDsXS59or0bir1blU8JiwlKc63hqSQbS8ZhzPF/QGh/VI+W9Vj7ixPCgijOQhQkWb+KMgBQ3MPQ+A4CRy/ZolMLV3TgsbNr3IFbkiRJUuOo6+J8VVFULfwzCFpxcEdxD8ecuA5z0nYU1xAUja7e7wnWAdkYCq2rMBemIgqSwVgEqh2KSwCaFr1R3EJQtHIsal1c0YFHkiRJkmpLURRU93DMuafAUACO9TfjSJhKEXmJWArTEIWp1oAD4OCB6haC4h6G4uR72a451BzJwCNJkiRJF2DJOgYOnuBw4d3Ga8MU/weUZFsHSnu0RHG2rizd2FPK/0tkdJQkSZKkKoiSbERBEhrfytsz1JXGP9Y6S8zO2boFhFuIDDsNTAYeSZIkSaqCJTMOdHoUj/B6v7bqFoomfJB1T7Gck/V+famyKzrwbN60mQVffUNaatqlC0uSJEnSOSyFKdaWlwYaR6M6+aG4h2E+uwVL3hnr4GWpwVzRY3jeeuNtADIzs5j58PQmro0kSZJ0OVE0OixF6aim0gabIaUJ6Yf51FrMCX+CokHbZiyKvdz/sSFc0S08Q4cOZd4nH/HgjPubuiqSJEnSZUYN6gnGYkzHf7Ou0dMAFFWLJnyI9YEwY4pbhijNbZB7/ddd0YHnnvvvZtBVA9FqKzdkmc1mziSe4bWX3+CZJ58lMSGxCWooSZIkNVeqkw/aViOgrMA6W6uBKIqCtt0NKE6+YDFgyTvTYPf6L7uiu7QuZMXylcx9ei4FBQW2Y0eOHOHHxYuasFaSJElSc2NdgDAUS9ZxND7RDXcfnR7Vuy3monQs6QfQ+LVvsHv9V/2nAo8Qgrdff4fPPvu8wvE5z8xhwMD+TVQrSZIkqVlT7QBLg9/GkpsAWMf1SPXviu7SOrcFB+DQwcO2sOPu7m47npaSip+/LyUlJZSVlcmR8pIkSZKNKExBdfZv+BupOtDYW/cJk+rdFR14Ppn3KUajEbC27uzasQuAzp078b95/yMw0LoD7meffU6H6Fg6te9Cx3adaBvZjkULf5DBR5IkSQKNHcJibvjb+MeCuQyRL8fwNIQrOvBs2byFCWMnsmjhD1w1YCivvPwqALt372HW9Fl06hTLvffdQ9u2bSudO/eZZxky8OrGrrIkSZLUzCh6L0RxVp2vYz6zGePeL7EUpVd9H3tXcPTCIgNPg7jix/AUFRUx95ln6dChA1qtltOnT+Ph4UF6ejqrVq3GaDQSGBjI4qW/0Da6DQDpaen07zOQpKQkDAYDdnZ2TfwqJEmSpKYgTKWI/GQUexeEqQxL9gkwFiJMZWAqBXMZoKDovVHcQlBdAi98LUOh9b8FyeDkW2UZ1ckXS/7Zhngp/3lXdODx9PRk0aJFPPnkk6xfv952PCcnh4emP0hhYRELv11IZmYmkydNYePm9Tg5O+Hr58uR44coKSqRYUeSJOk/zJK2HywGVN/2mOKWWEOOnbN1IUKN/T+bigosOSch8wjCpx1qQBcUVVPpWpqIoWAsBp3jBe+nOAdYr1NWgGLv0oCv7L/nig48v//+O927d+ett96ia9euALi5uaHX67nzrjvQ2eu4ZdrNjLh6JMXFxezetYd+A/pSVFiE0WjC3cOtiV+BJEmS1JREaQ6KnQuWwhQwFqE4eqF4tkJx8kVx8LAFGyEElozDWFJ2Iooz0IRfVWl1ZkVRwM7povdTXAJB1WJO2YUmpC+KekV/TDeqK3oMT/mCg126dGH37t08+OCDxMbGkp6eTnpGBgD+AX5cf/0EdDodH3/0CXfeejddYrvRs1uviw5athRlNOhCVJIkSVLTU73bIEpzEFnHUDxagarBkrwD87FlmA58i+n4b5iTdyGK0lB92qJpNQJRlo85cUOt7qdodGiCeyNyE2p9Dalq/5no2KlTJ7788kvWr1+PqqrMnvEwg68azF9r/mLU6JGYzWb27NmDxWJda+Gh6Q9a0/gFWDKPInJOgL0bqrNfY70MSZIkqRGpbqEoba7DnLQNkXMCTUh/lJahiJJsRFEGojgDS/YxSN8PqhbsXUEIREEyQlhqtfGo6tkSc+oeuMhnkFRz/5nAA9C7d28++OADLBYL+/btY+/evQDcfudttqADoNPpmDJ1MsXFxej1+grXEEJgST9gDTsADbSLriRJktQ8KPYuaMKvwnxmE+bTG1CNXVB926P+M/BYCIEoyUQUpkFZPkLnhOraok67rCs6RzBYB0crWvv6ein/af+pT+tJkybx1FNPAdYf0Kuvtk47z8vNs5VRFAWj0Uj3zj0ZO/q6CucLQyHm479hSbGu56MG9UB18mmk2kuSJElNRVEUNC36oPp2wJKyC/PpjQhhsT2n6n3Q+MagadEbbcQQVO+oOt1P9e+MKMvHdGwZwlRaHy/hP+8/FXgAXnzxRX744Qf8/f2ZNGkSAHOeepqAgACACuN27n/wfoxGI3FHj/HpvM8xHf4RUWwd+6M4+6M6+crFCSVJkv4jFEVBE9gFTegARE48ltR9DXYv1SUAbevRYCi0bTkh1c1/LvAAXH/99SQnJzNx4kTbsTlz5tgCUDm93pH2bTty7eixvP3W29aDjp4AiMJUTMeWYU5cjyRJkvTfYCnKsP3et2QeQZiNDXYvxd4NEFjObkVYTA12n/+K/9QYnnMpioKTkxM7duygsLCQgQMHsmnTJgBiY2O55/57uOfOe2zl73/gPj5Zm42ic0SU5mIpTLH+EOYmYEpU0YbKzUclSZKudJb0g/8+MJdhOrYMTXAPFOeAWo3ZEUJAaS4gUP75g9p2r8I06/9o7EH2JtTZfzbwlCtfnwesCxUCTLttaoVxPeWrMH/69xsAKA7uaBzcUV2CMB352dq06dW6cTaXkyRJkpqMJmwAWPqiaHSI0lxMiesxx68CjR2qRwSqT0y1Fww0nd6IyEkAYd2nS9t2PIq9q3UQdE485rNbUZz80LQcVuVChlLN/OcDz7nee+898vPzeeqJOXh7e9uOf/n5l7z+9muVyiv2rmg73IzITaiUzCVJkqQrj6KooLG25CgO7mhbj0EUZyDyz2DJjMOSdRw1oBOqT7uLtvgIIRDZJ/69rt4H7KxByZK8A0vGIRT3MDQt+sqwU09k4DmHqqq8+uqrfPfdd6SlpdGlSxdSUlJ4Ys7jFcoJIf5J5AooGlTPyKapsCRJktSkFEVBcfIFJ19Uv45YUnZjSd6JyE1EEz4YRae/4HnadpMQhSmAQHELRVEULPlnsWQeRfXrgCagS+O+mCucDDzn8ff3p1+/fmzfvp2ioiK6dOmMh6cHZ06fxXj4JzAUVHmeJqQ/qmfLRq6tJEmS1FwoqhZNUHcU9zDMp9ZiTliDJnLUBVt6FJ0jikcEoiQHDEXg4IYlJx6EGcUlqJFrf+WTgec86enp/PXXXwAcP36c+Ph4OsZ25MUXXrrwSTq9NeFLkiRJ/3mqky+ED8Z87Dcs6QfQ+HWsspwwFmM+uQZRkgkoqAGd0fi0w1yWbx0XFD4Y1TW4cSt/BZOB5zx+fn5cc801LFu2DLPZTI8ePWxhR3ENRhM6EEWja+JaSpIkSc2ZqvdB+LTDkroP1bO1deXkcwhTGeZT6xCGAjShA7BkHcOSutcakAK7YT6zCUvGIRl46pEMPFVYuHAhn3/+OTqdjs6dO9O7d28mT5nMogO6i+6vJUmSJEnwzwbTGYcAgenoL6ju4Sh6HxS3FihaB8zJ2xGlOWjCr0J19kf1iECYyjCdWIH5zCawd0MT3LupX8YVRQaeKjg7OzNjxgwAXn/9dZycnHhizmP8MOmdpq2YJEmS1KwJYcF8aj0i75TtmOIcgKUoDbLi4IyC4haCyEtEDehSYTkTRWuP6tUGS9JWVM+W1Z7eLlWPDDyXEBQURFFREU8+OgdhcECxc6qynDAbEAXJqO5hjVtBSZIkqdmwZB61hR1Niz4onpG2ngFhKsWScxJL8k7r47zTCN+YCoOaVa/WUJaL6iEnwdQ3GXgu4cYbbwTg4YcfxpSejerbDtUtFBzcwWIEUynCXIbl7DZEcQZK1FgUR4+mrbQkSZLUJFSPCOvsK9cWKGrFj1hF64DGJxrVPQzToUWI4gxM+75G23aCrTVHUTVogns1RdWveDLwnKO0tJSPPvqItWvXkp6ejqqqREdH06tXLzZu3EhUr2uwpB/AkroHVB1YTEDF5b7NpzeiaTWi2gObzan7sOSeQhs1Ro4PkiRJuswpWgcU9/CLl9Hp0bafjOnYMijLx3TkJ7Ttbrjgmj1S/ZCB5xy33norP/zwA4MGDSIqKgqjyciePXv4+uuvURQFofdHDewBOkcoyQaNDsXBA0VjZ51emPAnoiQLkX8ai9mIKEhG0elRA7tWSvpg7eu1pO7+54EZFPntkCRJ+i9QNHbo2o7HnH4AS/JOLGkH0AT3aOpqXdHkJ+w5UlNTcXNzY/ajs4hoaU3oK5av5Pfffmf16tWQ///t3Xd8VFXeP/DPuXcmk0ky6T0EQkgIVZoKgnQBQRALVlTclVVZVtd1XX+KXde26vI8j7vrCvZ1FQugIggqRUIRlN5bQkJIgPRCyszce35/DBkYEiAJmZLJ5/168WLmljPfCZnkw7nnnpMHWZEHCBUiJA4iJAFQVMAcAWi1zna0nNXOxxKOa7JSqLAfWAy1y1goQaeWrbBWnX5xzQo0EoqIiMh/qbG9IYLjIIyNjw+l1sPfsGf46KOPMGHCBFx7zWQMuXIIflr101lHCBi63wi9/DBk1THox7cDBZsAxQBhCgeE6lwEDgCgBgDGINhz1wCQgFYHLWc19KAYQAjIkgOOVkMSAIMZRETU/iicuNYjGHjOkJycjF9//RVvvfUWnnvuuYYHqAEQJgvU2N5AbG9IqUNWF0OePAZZW+EIPNZywF4HR8CxArp2aryP1RFqVBNk6SHgjPE6aspIjt8hIiJyIwaes5hMJjz00EOw2Wx49NFHMWTIELz74VxkHcrGpJn/cDlWCAUiOAYIjnHZLqUE6ioc43mqiyDrKiBt1UBdOVBTDEACpjCgrhKQGvSybKjR3RqtR0oJ6HbO7kxERHQRGHjO4ZFHHsGCBQuwdu1aTLx6EmpqaqGX1kFGpgPWKijxfRsdiAw4VsFFYBhEYBgQkercLu21jlH51iqgtsy5Xc//BYolqdFJpvTivdDzfoYITYbaeTR7goiIiFqg8SVcCUIITJo0CUlJSRg5chTuueceSOtJ6HnroZ/YAVma3fw2DYEwpIwCcCq0mMIcf0tAy10NKfWG5ygBjkMqjgD2mpa+HSIionaNPTznMWvWLMyaNcv5/KUv90IvPQS9YBP0ijxoRXuAugogKAqGxMsg6u++Og8RFAWl41Douasdq6yHdoAwhkDP3wjt4FKoHYc6e3qkvQ5aruOOLyWuD+doICIiaiEGnmYQAcFQ4y4BpA69cLdjpmWpA1XHYD+wxDHhYFCUyzThjVEju0AI4bh9XTVCSegPYQyEdmQt7HsXQIQmA4oRsirfeY4S1tHdb4+IiMhvMfC0gBrfF0pcH9i3/8exwWAG7DXQDnx7+qDAcBi6jDtnr4wSkQooRmjZP0Iv2AwRGH5q5mZAlucAxiCXiQjt+xc5zovuDrXDIOinbotXO1zBBeaIiIgugGN4WkgIATXtasAc2fjYmtoy4AI9PUpYMoQ5CnrhLsAcDbXLOMfcPQBgqwbs1Q1Pkjr06kJoB7+DrDyKs5e2ICIiooYYeC6CEhwLQ9drYciYDIQkODaawgBzFBDWEVrxQeiV+edtQ+08CjCFQtv/DfSSA1A6XgmlfuE43Q4R0+P0waYwQDFAO7DE8TwgBMIUCllXCf3kCTe8QyIiIv/AS1oXSQgBmCOhRqRCqy11jOs5Nd+OLM+FBgBp46GExDd+fkAIDBnXQS/eB71oD2RpFiAUx1TjoR2gRKYBURnQK/IgK/OhlxyACImHrMwHrFWwbX2/viUofe/20LsmIiJqWxh4WokS1dWxZpaUjkVET625BalBO/gdtMAIGDqPgjCFupwnpYSWuwZKVDoM3W4A6sqhVxZAVuVDP7YVesFmKAn9ocT2gojt5TzPfnglZNlh53M1ZYSH3ikREVHbw8DTyoQQMKSOgdRsjqUnivZAP7YFqC2FtNU4LkFJHUIokLodsiIPsiwLum6HwZIIBIZDDQwHYrpD2uugn9juuA2+Mh+oLoLaeTQUSwLUuL6wnwo8wpIIaa91tktERESuGHjcRKhGCAAyPAV6yQHAWgUtbz10xQBZXegY0HzGRIPKGb03zjYMJqiJl0GYI50rsGuHlkILioUhdTRERCpkaRZkbZnjEldtOdQOAz31FomIiNoMdge4mQgMh7HHTTB0vRbCYAIAKB2ugAjr5DxG7TQcSkjcudsIT4UISYRzhubqE9DyN0FNHgIl6XIYMiZDie/nWIbiZGGjMzYTERG1Zww8HiKComBIGw9D14mOhUIVx2KgSnQ3iHMMaHaeK4Tjbq7A8NPbQuIhFAOUqAzopVmQdiugmqAd+BZa7hp3vhUiIqI2h5e0vERN6AcNgF58AHrRXojgWAhLBwhLYqOzNQvVCLXDFdAOOm5Jl1X50OrKIKtLICvzHAeFJEKISMe6W0REROTEwOMlwhgEQ8chkEmXQZblQC/PhX5iB3BsM6AEQIQmQQlNdgQgoxkAoITEQUZlQC/eB1ly8KwpBwVEcCzk8a0AAKnbz7maOxERka+z1dnw348/QWxcHPr0vQSxcTEwGo0tbo+/Eb1MqAEQUenQK/MhgqKA4HhAakBlgXPhUJjCIAIjIALDIGvLnOcqcX2hF+4EpASk5liSop6tBuCSE0RE1EaVlJbglZdfbXTfvff+DlNuvrFZ7THw+AApJWRZluNJ1TEAAsIcCRHdw7FOVl0FZG0p9OJjgG6HknQ5RIAFIjQZakI/x+3t5Ueg5axytmk/vALGjMleeT9EREQXKy4+Dhs3/4y/vfQ6vvzyS5d9c+bMxZw5c5vVHgOPDxBCAKoJSkwPKBFdHDMqnzwOWbIfUrcDigFQTYAhCMJggqw6AamWQh7bChEQAiWmB/SKXNdGa0ogrSchAoK986aIiIguUmhoKP76yvP46yvPO7eVlZZj5fJVmPX4LEjZ9PUkGXh8hcEE2GshTBYIUwaU6AxIex1k1TFIaxWg1QH2OshTf8NaBWGyQFYdh1Z/KevUzM717Pu+hqH7DRCGQC+9KSIiotYVHhGG66dMxvVTJmPnzl2Yct1NTTqPgcdHKCHx0MuyHbepn7r9XBhMEOGdznue1GyAtRJSt51eVNTlAK6mTkRE/knUz0/XBJyHx0cocX0B1QT7/kXQK442+TyhGgFDILSczIZtdhrmvMOLiIjI3xTkFzT5WAYeHyECgmHoOgkiOBZabiakZm3SedJeC/vuLwBrZYN951qhnYiIqK1bvWoNHn3k/zX5eAYeH+KYXHAwYK9xrI11AVLXYN/5qcuaXC5s1a1cIRERkXdZrVb843/+iXun34uePXs2+TwGHh8jTBbAFNqkwKMf2+x8rKZPhKH7FJf9Uqtr9fqIiIi8ZfvW7bhm3CT8+99v49VXX0VpaWmTz2Xg8UFKWEfH+li2mvMeJ2tKnI+1A9/Cvu9r1wM0uzvKIyIi8riD+w/int/8Dp07d8bHH3+MDz/8EFlZWU0+n4HHBymxvQEhoB/fdt7j1OQroUR3P71Bt51uI77fBRclJSIiaguO5Obht3dPR7du3fDCCy/g/vvvh8lkwt//540mt8HA44OEIRBKRCr0oj3Qz3NpSwQEQ549TkcJgKH7jVDj+zomNCQiImqjThwvxMsvvIIJ465BQEAAxo8fj2uuuQbp6el49/25SO6Y3OS2OA+Przo1F492aBk0AGryYChRGQ0OUzsNA6xVgFAha0ogQuI40SAREbVpNpsN//PG/+G9d99zzqacl5eH559/HtOm3YU//ulBBIUENatNBh4fpYSnQs/7+fQGY+P/sEIxOMOR4GKhRETUhkkpkZ2dgynXTUF19ekrGAkJCbjzrjswYeIExCfEtahtBh4fJQwmiOBYyNoyx3gcSwdvl0RERNSq9u3dj61btuH4sWP41z/farB/zJgxuHPaVFx6+aVQlIsbhcPA48OUmB7QDq+CfnSDY/V0DkImIqI2zmq14j8ffIzX/vZ6o/v79OmD+2bcixGjhl90yDkTA48P08uPOB+LgBAvVkJERHTxcrJzcPXYCS6rnL8x+w3ExESha0ZXhEeEu+21GXh8mDBHQZYeAkxhgBoAaa/lgGQiImpTpJSQ1YXQDizGuDHvO7dv3PwzQkNDPVYHA48vs510/F1XDvuO/7rsUtMnQgRF89ZzIiLySXrlUegndkFWnl4QW1EUfPrZJ+jdp1erXq5qCgYeH6bE9IBeuKvRfdqBbwFjMIw9b/ZwVUREROcmpYQsOwwtZ5XLdjVtPHYtbfpEga2NgceHiYAQiIgukGWHoaaNh1BUIDAcsjQbeuFuKHG9vV0iERGRCy37R8iKPACACEmAEtkFIqwThBrg1boYeHyd1CCC46AExzg3icguUCK7eLEoIiKihrTCXc6wY+g+xafmh2Pg8XVSBzhOh4iIfJSsq4BeUwZZlgVZlg0AMGRc51NhB2Dg8XnCGAy9qsDbZRAREbmQUoeWvQKy4ojLdrXrRAhzhJeqOjcGHl9nMAG2Gm9XQURE7ZzU7dCyl0PaqqGEdYJelg3UVbgco3QYDCUo5hwteFeL7wnbv38/QkJC8MEHHzi3/e53v4MQwuVPSkqKc/+aNWvQu3dvdOnSBQsXLnRuHzFiBMLDw5GXl9fgdZ599lmXNtobWVMGYfLcPAVERESNkhKyMh+oLYN+fFuDsAMInx5f2qLAY7PZMHXqVJw8edJl+/bt2zFr1iwUFBQ4//zyyy/O/dOnT8cLL7yAuXPn4oEHHkBdXZ1zX3l5OX73u9+18G34Mwko7IgjIiLvEqoRhm7XQ5zdg2MIBMyRUDoOdSxo7aNaVNkzzzzTYHZEKSV27dqFxx57DPHxja/5VFdXh379+iEkJAQGgwFWqxUmkwkAkJqaiqVLl+Ldd9/FPffc05Ky/JKsLYViSfJ2GURE1I5JaxXs2cuBmhKX7WrKSMct523g5ppmB57Vq1fj7bffxtatW9GxY0fn9kOHDuHkyZPo3r37Oc997rnn0K1bN9jtdjz//POwWE6P4B46dChGjBiBhx9+GGPHjkVycnJzS/M7UrcDdRUQsZxvh4iIPE/qdtgPLAFqip3blNheEJZEiJAECOHZ2ZIvRrMCT1lZGe688068+eabDQLJjh07AAD/93//h++++w6KomD8+PF48cUXERYWBgC46667cOONN8Jutzu3nWn27Nn44YcfMH36dCxbtqzJdaWmpro8P/NSWVsma8scDwJ9b7Q7ERH5P73koDPsKImXQonu4ZgEtw1qVuCZMWMGBg8ejNtvv73Bvp07d0JRFCQmJmLRokU4dOgQHnnkEezcuRMrVqxwrpkRHBx8zvZDQ0Mxd+5cXH311Zg7dy7H9JzqOhSB4d6tg4iI2iUlqiuU0A4QASHeLuWiNTnw/Oc//0FmZqazJ+dsTzzxBH7/+98jKioKANCrVy/Ex8dj0KBB+OWXXzBw4MAmvc64ceMwffp0PPLIIxg3blyTzsnKynJ5vnnzZgwYMKBJ5/oyaT0JGMwQqtHbpRARUTskhAL4QdgBmnGX1nvvvYfjx48jOTkZISEhCAlxfAHuv/9+jB8/HoqiOMNOvV69egFAo7ebn88bb7yB8PBwTJ8+vVnn+R3VCOg2b1dBRETU5jW5h+fjjz9GTY3rBHjp6el4/vnnMXXqVNx1113Iz8/Hjz/+6Nxff0t6z549m1VUaGgo3nnnHYwdO7ZB7027YqsGDGZvV0FERNTmNTnwJCU1fmt0bGwskpKSMGXKFEyePNkZgPbv34+ZM2fi9ttvR7du3Zpd2JgxY3Dvvfdizpw56NSpU7PP9wdSswJS83YZREREbV6r3U927bXX4vPPP8dXX32F3r1745577sENN9yAd999t8Vtvv766+027ACArDgCJSLN22UQERG1eRc1JaKU0uX5TTfdhJtuuqnZ7axatarR7RaLBYcPH25BZW2f1O2AvZbLShAREbWCtjNjUHtTv2BoQJB36yAiIvIDDDy+SrM6/lYDvFsHERGRD9J1Hbt27Gzy8b67ylc7J3U7AEAonIOHiIjoTEWFRXjwDw9h86bNTT6HPTw+SlYXOh6oJu8WQkRE5GPumz4Dmzdtxt/+9rcmn8PA44OkrkE/sRMioguEkfPwEBER1cs9nItdu3ahY8eOGDVqVJPPY+DxQXrRbsBeA8HeHSIiIqdfNvyKm268Bb169cL69eshhGjyuQw8Pkbqduj5vwIARGiil6shIiLyDYezD+M3036LIUOGYM2aNUhMbN7vSA5a9jGyMh8AYMiYDGGO9HI1RERE3mez2fDyi68iJSUFCxcuREBA8+9gZuDxMXrZYSAwnGGHiIjolCcfewo/rfoJixYtalHYAXhJy6dIWzVkeQ6U8M7eLoWIiMgnHD92HF9//Q3+/e9/Y+LEiS1uh4HHh+ilWYCUUKK7e7sUIiIin/DNV4tgNBoxderUi2qHgceHCNUISA2yttTbpRAREfmE4OBg2O12/POf/0RBQUGL22Hg8SEiMh0iKBp6QdNnjiQiIvJnvXr3hJQSjz32GP7xj3+0uB0OWvYhQigQliToJQe8XQoREZHXHcnNw4z7ZiImJgYzZszAww8/3OK2GHh8jWblgqFERNTuHc4+jKvHTIDZbEZubi6io6Mvqj1e0vI1xiDAehJS6t6uhIiIyGuOHzsOAKipqcGll16Kn3766aLaY+DxMYolCdBtkGWHvV0KERGR1wy8YiB+/mU9PvzP+4iPj8cdd9yByspKAMCaNWswdOhQrFu3rsnt8ZKWrzFHQoQmQ8vNBAxmKJYEb1dERETkFeERYRh4xUBERUVh4oRrcf/99+O///0vPvvsM6xZswZr1qxpclvs4fExQgioKSMB1QRZftjb5RAREXldXEIcLrnkEnzyySfIzc3Fa6+9hr/85S/NaoOBxydJwF4LSG/XQURE5DmbftmMbmk98MyTz+HzT79wbq+tqcPDj/wJBoMBN998M5599lls2bKlWW3zkpYvEiqUqK7Qi/dCie0BYQrzdkVERERut3aNY0zOZ/M+g9FoxPU3XIdPP5mHl158GQDw/x57FK++8jds2LCh2W0z8PggIQSUpMuhl2ZBL82GGt/X2yURERG53XdLvnM+ttlsGD50JMrLyzFgwABs2rQJr77yN0RFRWHYsGFIT0+H0WjEiy++2KS2GXh8lFAMECHxkFXHvV0KERFRq6uurkZJSSmMBgOqqk5iz649OHbsmHN/RkYGbr31Vlx55ZV44403nNt/+OEH9OvXDwAwc+bMJr8eA48v02yAavR2FURERC6k1CErj0JYkiBEw+HAUkrIilzo5Ucg1ABA6pBSw+03b0N1dTWqq6uRm5vb4LzRo0dj7ty5CAsLQ2RkJADg4YcfxtKlS53HVFZWYvHixc1eOZ2Bx0dJey3kyeNQk6/wdilERESuasuhZf0IGINg7HkLgFMhp6YIsvgA9PJcwF4DGIMgFSMgFAhFRVRUJ1RWVmLcuHEYOnQokpKSYLPZEBwcjD59+sBisTR4qRdffBFjxoxBamoq0tLSoKqqSwBqKgYeX2WtAiAhzFHeroSIiMhVYLjjb1s1pNShF+2FXrwPqC0DjEFQIlIBSxJEcByEboO0VUE/+gt++OEHAEBRURHmzp0LAPjyyy8xbty4c76U2WzG+PHjXbZdffXVkFJi8+bNGDBgQJNKZuDxUVJK59/Cy7UQERGdSYjTv5ns+74BaksBYzBgCgNMYY4ensJd5zy/tLTU+Tg2NtattdZj4PFRetEeAIAwR3q5EiIiIld6adbpJ7WO8CKMQUBACGCvhbAkAeYIQLNCVhdBlue4nP/AAw/gN7/5jXPwsScw8PggvSwHsvQQoJogFNXb5RARETlphbugH914ekN4KtTobhABIZC1pY7LXHWV0Iv2AbUljmMCLFBCk7B20fu49NJLYTB4Pn4w8Pig+t4dNXW0lyshIiI6TTu+HXrBJteNZVnQKvMAzXp6myEQwpIEJa43hCURwhAIABg0aJAHq3XFwONjpK0asqoAavIQKMFx3i6HiIjIyXkpSxgASEBqEJHpECYLRGA4RGAkYDQDQnUZ5+MLGHh8zqlvEDXAu2UQERGdRUnoD/3YFqDGcanK0ONmiIBgL1fVNFw81NcYAoGAEOhFeyDtdd6uhoiIyEkN6whjxmSoqWMcG3SbdwtqBgYeHyOEgNphMGRNKez7voKsLvJ2SURERC5keQ4QYHHcht5GMPD4ICU0CYaMyRDGINgPLYM81XVIRETkbtqJXdBLDp1zv6yrhF683zFmx8fG6ZwPA4+PEgHBUFPHAgEWR+ipq/R2SURE5KekvRa2vV/DtvV96PkboZ1n0kDtyLpTJ+keqq51MPD4MGEwwdBlLCAM0PI3XvgEIiKiZpKaDfadn56eMwcChrTxjR8rJWRVPgBAiUzzUIWtg3dp+ThhCIQS2xN6/i+QmtWx6iwREVELSc0K1FUA5igIIWDf9xUAQESmw9DxyvOeq5/Y4Xwswju7s8xWx8DTBigh8dClDllbBhHsmTVHiIjI/0h7raM3B4CwJDmuIlirAABqhyvOe659/yLnjTRqyog2NX4HYOBpE2RdBQBAGMxeroSIiNoiWVcB+575rtsqj7o813JWQ+00rNEljfSTxx1hxxQGtctYKAEhbq3XHTiGpw2Q9hoAgpMREhFRs0h7LbSjGxqEHQAQQTEAAEPGZMex5Ydh3/0F9JIDLvPAaSd2QTuwBDAEwpB2dZsMOwB7eNoEJSwFet4G6KWHoMb08HY5RETUVthqoBfuPvVEAJBQkga6/C4R5kgoCQOgF2yCCI6FlrvGsT00GWp8X+jHNjsedxoOoRo9/x5aCXt42gBhNEOEdoAsz/V2KURE1JYYg854IgG4DjyuJ05NIKjG9YGaNgEiMg2y4gjsBxYDBnObDzsAA0+bIUwWSGsVpJTeLoWIiHyUlBJ6xdHTvyuUhhdylIguDbaJ0A5AYDjsuZkQ5kgI4+n1sdSUEW0+7AAMPG2GCEkArJXQDq+C1O3eLoeIiHyQLM+BlvU9ZEUeADgGIJtCIUKTTx9kbLjYp1BUGDqNAKxV0A4tc17mMnSfAiUo2kPVuxcDTxuhhHWEmjIKsjIP2oElkLYab5dEREQ+RNaWQTu80vHk1E0uUkrAetJ1VuS6ikb/4yzMEVBTr4KsLoRengM17pI2sxJ6UzDwtCFKeCcY0iZA2qody01YT3q7JCIi8gHSVg373oUAHFcElJA4SGsV7Ns+AKTmGMsjHL/y9aLdjruuGnNqGSNZdcITZXsUA08bI4KiYEi7GrDXwr5nPrTjDQefERFR+yIrjpx+EhAC/eQJOO7KOrW/5AAgJRAYCSWmB5TEyxptRwQ7blWHrQpSs7mxYs/jbeltkAgMh6H7DdALNjlvI1RC4rxdFhEReYm0Vp9+XHIAWskBx4DlwAigtrR+D1BXDmkOh35oKTRTKAzdrocQZ/R9nDnfm/CvPhH/ejftiFADoCReDgDQDi6BXprl3Cel9LtkTkRE56bE9IAS3c1xg0s93X5G2AFEWEco0RmQ9b8v6ipcxvZIXYP94DLHk5D4RmdcbsvYw9OGCUWFmjwY2pG10HJ+ghKR6lgn5eB3QG0ZYAwBIAHbSSAwAsJohqyrgBLaAUrSQNdUT0REbZYwmKB2uAL2/YvOeYwsz4Xa+w4o8f2hZf8IWXUMsvIoRFgnAIB+bAtQVwYAUOP7e6Jsj+JvvDZO1hQ7H+ulWZC15Y6wAwC2KkfYAYDaUsjKfMBaBb1oL1BT2rAxIiJq09ROw08/ESrUjlcCAZYztikQqhFq8hAAgJa9AnplPvSTJ05PSCgUvxwmwR6eNk5EpkFRTZB15dByfnLdZ0mE2mGw44luh9TtEEYzoJr8YhIpIiJyJUyhMPSZBvu2DwGpQTuyHoaMa6EX7oJevB96wSaoSZcDhkDHDP4VedALd0PaTo8BMvS6zYvvwH3Yw9PGKUExUBP6w5AyEmraeJdBZrKmBJAahMkCYY6AEhwDERDCsENE5MeEUGC45E7HE6nBvnchlKhujt8Pp2ZgFmqAczkJGAKBmmKIsBSo6RMh/HShagYeP6KExDu7KQEA9lpouWu4HAURUTsjFAMMfaY5x+fY938DSB0iIhUAoFfkQS/cBSXxcsiqAojQZBg6j4RSf1u6H2Lg8Tf1vTfmKACArC6ELMv2YkFEROQNeuFuKJHpLnPuyKoCx77jOwBjiOO5tQpKWEdvlekxDDx+RoQkAKoJSlC0c/CarB/ETEREfktaT0LWlkPaqh2LiOb/Ai37R0CrgxLd3XFMZT5kXQXkyWOOyQVrih137TayoKi/4aBlPyPUACiRadDLcmBMHgzlVPclERH5Lyl12Hd/7rJNWBIhK/OhH99+eltQNOTJ08tGqB0GQwlLRnvAHh5/ZAoFbCehVxz1diVEROQBQihQU0ZCnLGyuazMdz1IMUCJu+T0zS2GQIjQJA9W6V3s4fFDSmQaZEWeoyuzwyCI0GTAYIYQ4sInExFRm6SEp0AJTwEAyLpKx8Dkoz87dgoF0O3QSw5Cie4GpaYUIiy5XU1Ay8Djh4RigNp5FLTcTGhH1jk2GsxQOw2FYmk/aZ6IqL0SJgvUmO4AJPSjG5xLSOhHN0CWHYYISYAI8t87shrTfqJdOyOEArXjMBi6T4GaMgoiMAxazmreok5E1I6oMT2gJA0EjEEQliQocX0BxQD9+FaXsTztAXt4/JgQAjBZIEwWyLpyyKpjjsXkOPEgEVG7ocb0gBrTw/lcSh32bR9CP7EdSsgYL1bmWezhaQdkdRH0Y1sgIrtylmUiovaupgQAICvyIHXNy8V4DgNPOyBrywCpO9ZPISKidk3WVZ5+XHHEi5V4FgNPe2ByrJTb4BZFIiJqV6StBjCYnM+1o794sRrP4hiedkAExUBYkqAdXgl0uAJKdIa3SyIiIjfTjm0DtDqoSZdDSglZdhhazirXg2xVkLodQvH/OMAennZACAVq6lVQortBy1sHrWiPt0siIiI3049thl64C/rJQmhZ37uGHXOk86GW97Pni/MCBp52QggFStJAKNHdoedtgF6e6+2SiIjIjURIAgBAO/Cty5AGERIPNb4foAYAAGTJAUhbtVdq9CQGnnZECAEl6XKI0GRouZlcVJSIyI8pkWmAemq8jlCBwAiI8BTIqmPQspcDmtV5rH33l16q0nP8/6IduRBCgZp8BewHl8J+aBkM3W9sF9duiYjaGyUyzRF6ziBt1bCXHXbs73AF9Lz1gMEMQEDWlUOYwjxfqIewh6cdEsYgGFKvAuy1sO/9CnpFnrdLIiIiTzJHOm5JNwRC6TgMsFdDP7HL21W5FQNPOyVMoTCkT4QIsEDLXg699BDkqbVWiIiobZJ1lbBtfR/2/d9CrzoOe+4a2HZ9BlldBODUf3i7TgI0G2TVcagdroAS6OjVqT/GX/FaRjsmgqKgpl4FLWcVtJzVEBV5MHQa7u2yiIiopYxmAICsLoR2cIlzs33/IghLIkRgJPSSA4DRDEPGZAiTBXplgeOgwAhvVOwx7OFp54SiwtB5NER4Z8iKo1xclIioDROKAYY+06CmjoEwR53eEWABIKCXZUGExMOQNgEICIasKYF2aCkAQInu5p2iPYQ9PAQAEMGxkGXZsO9dADV5CJSQeG+XREREzaRX5kPP3wS1yxgYMq6FdmIn9PxfHJeuQpMAAFKzOu7K0uqct64DgAiK9lbZHsHAQwAAJbo7RGAYtLwN0I9uhJJxrbdLIiKi5tLtkDVFsO/8FEpUBvTifafm29EhpQ4hFMB6EtDqAMAxVUmApV0sLM3AQwAcc/QISxL0gBBvl0JERC2khHWEjO4OvWiPI+wAgGaFlvWjY04erQ6G3lNh7Psb7xbqBQw85EKYQqG3o9VziYgaY7I4ll6oqyzxciUXJqUENCvEqUVB1Q6DHLPqVxVAlh6C1DWgptjZqwPF/3tzGsPAQ65UI8Db04mI2gy9eB/0vPVQEi+FXrjb8TNc6i4zKdcz9LkbQggvVOl9DDzkStcAWzX0quNQQuKcd2211w8IEbU/9b079Y99vpfHWgkA0PN/hQhJhAiJA4QCCAX6iR2AvRYAYOgzrV3/LGfgIRdKZBrkyRPQDi6BZggE7HUAJBAQAsWSBASGA7odwmSBCOvkGABHROQnzgw7Z27z1dCjlx6CfmInAEBNGXnq5/LpUKPG9vJWaT6HgYdcCHMk1LSrHdd9bTWAwQQhFMiaUuiVR4GSA47/Oeh2wGCGEt8XSlRXBh8iavMaCzu+Sup2aDmrIctzIMI7Q026HMIY5O2yfBoDDzUgFANEVEaD7SpwemLC2jJoJ3Y4Fp7TrFDjLvFskUREbhSa1BUAUHF0PwDf6uXRK49CL9wNWXUMaqcRUCI6e7ukNoGBhxqlVx2HfmInlMQBEFKHBKAf3wERFAWhBgCmUBg6DYNdUaEXbIKsOAJhSYKon7BQ6hDmKOddA0REvuzM3p36sFP/uD70+AKp26Ed+h4AoKaMgBLOsNNUDDzUKD3/F8daLBW5LttlWZbzsegzDWqHKyCD46GX5ziuIx/bcvpggxmGjGvZzUpEPu1cYaex4zzdyyOlDll96pZyISBrK5z7RFiKR2tp6xh4qFFql3GQtWXQcjOBunIACpTYXpC2KghjkMuAZRHZBUpkF8dq67XljjE+UoP90PfQjm6AIWWkd98MEdE5NGXcjjd6eaRmg35iO/TCPYBua7Bfie3dru+4agkGHmqUUI0QwTEQ3a53BJ6AEAjl/N8uQiiA+fRqu8KSCFl6CLadn0IEWKDE93Ou5UJE5G1nh50ze3fCExxrTJUVFDQ4x929PHpZDrRT4yOVmO4QYSmnesqlY34dg7ldLAXR2hh46LyEEI5b0VtAiUyDbq8BhAroNmhZ30PG9HD8z4SXuYjIh5zvUlb9fk/08uilWdByfoII7Qi1w0AILvfTahh4yC30koOA1KF2Hg2hGCClhH5iB/Tj26EXH4DacQgH2xGR15xv3E5970794/penvrQ465eHr3qGLTcTIiILlA7DuUlq1bGyVOo1UkpoeWugXZkLezb/wNZVw4hBNS4S2DocRNEaAdoh1dBy/vZscYLEZEHnW/czplh53zbWnPOHiklbHvmQzv4HURwLNTkIQw7bsDAQ82ml+fAtvV9aEc3NLpfCAFDxrWnjy/cc3qfwQS103AoHQZBL94HLftHx2BnIiIPaMq4nfO50KWvFqktA+ocd1+pyVdCKGrrvwYx8FALqI65dfTyc6+qLsyRMPS6DYaMyVAS+rvuEwJqdHeoqVdBVuZDP77NreUSETXmfOElPDYY4bHBp5+7sZdHytM93cJkaZU2qSGO4aFmU0LiIfpMO+d+qduh5WZCCesIJaLLuduxJEHG9YF+bBuU6O4QhkB3lEtEBKAZ43bODDqxwSg7cdLl2NYawKyXHYZetBey6hgQGAFD6lUX3SadG3t4qEWEUM69fpauQZYdhpazGnpF3nnbUaK7A5CQ5+ktIiK6WM0dt9PU41ray1M/z5m0VUNJuhyGrpN4R5abMfBQqxMGE9TkIQAALesHSOvJcx9sCASECqlZPVQdEbU3zRm3c2bvzvm2tXQsj1awGbat78O+72vAGAxD2tVQY3pw3I4H8JIWuYUS1RUICAbsVuCCc+5I6PkboYQmQbRwzh8ioqa40Lideh3iHb0teceqXI854zb1ek25TV3qGmCrgl68DwCgxPWFEtvzghO6UuvhV5rcRrFceFZlIQSUxMugH90AvTQL6lkDnImILkZLxu005kJjec4OPbK2DFr+rxDBcQB06Me3A7odUAOgdh4NJaxjS98StRADD3mdGtMDsvwI5MlCb5dCRO3E+cbt1Pfu1D9uSi9PPanZHNNtVB1zPK84AkBARKZBCesEERLPZSG8hIGnGaxb3vN2CX5r7ty5uP/++7HgxeuRnpHm7XKIiM7Q2K3i9T+nhrpsLS4uwZCBHzufBwcH4/jx4wgOPn8PErkfBy2TT5g6dSr69euH22+dikMHs7xdDhFRi0RFRSImJgYAEBsbi3feeYdhx0cw8JBPCAoKwqpVqxAREYEvP5/v7XKIiFpkzlvvoLCwEEOGDEFBQQFuvfVWb5dEpzDwkM8ICQnBLbfcgi+/+BI/fr8cUsoWtXPieGGLzyUiagmbzYad23fi/ffex+jRo/HDDz9AUfgr1pfwX4N8yuOPP47x48fjD79/AK+98kazzz9WcBzDhgxH9/SeWPHjSthsNjdUSUTk6vmn/4opN9yMiIgIzJkzB2az2dsl0VkYeMinRERE4NlnnwUALFywsNnnx8XHomNHx+2ev79/JoYNHoEvP5sPTeOq7ETkHjnZOfjiiy8AAPv27UNqaqqXK6LGMPCQz4mIiAAAXHb5Zc0+VwiBf899CwAQHx+PiRMn4qknn8ZdU+/Gvr37Yauz4f7pv8dDDzyMtavX8dIXEV20t/71NgDg22+/hcHAm599FQMP+Zz4+HhceumlOJJ7BOVl5c0+f83qNQgJCcHu3bvx0UcfYe3atairrcPkidehX58BWLVqFfbu2Yt7fjsdf/zDn2C3293wLoiovfhq4VcAgKuu4uKfvoyBh3zS+++/j+LiYtx5+zSUlJx/yvazHT2aj4yMDGdP0RVXXIEtW7Zg/vz5eOqpp/Djjz8iOzsb77zzDr5f9j3WZK6DruvueBtE5OfmvDUXAHDHHXfAZDJ5uRo6HyH9sE9/8+bNGDBgADZt2oT+/blUQVt18OBBDBs2DJ07d8bb774Fo7Fps5P++aG/QOoSixcvPu9xZWVl6NOnD3Jzc2E2m5Geno7BQ67ANROv4eSHRHRB1VXV6N/3UgBAaWkpwsPDvVtQO9Sc3/fs4SGflZaWhi+++ALr1q3DquWrmnxeSXFJkyb6Cg8PR1ZWFlasWIEXX3wRl112Gf791tuYdM21KCtt/qU0Impf/u9//wEAuOyyy2CxNDYbM/kSBh7yaUOGDEFSUhIOHjzU5HMGDrocixcvxsGDBy94rKqqGDlyJP70pz/hnXfewSWXXAIAkJKXuIjo/Pbt3Qej0YgffvgBqqp6uxy6AAYe8nl9+/bFnLfn4vvvvm/SXVXTfnsX4uLi8PTTTzf7teLj4xEdHY2IyIiWlEpE7UR1dTXWr18Pm82G0NBQb5dDTcDAQz7vv//9L6ZMmYIHH3gIn34874LHm81mTLnpRixZsqTZEw/a7XYMHnxFS0slonbiSG6e87EQwouVUFMx8JDPCwsLw0cffYRbbrkFixcvadI5Ay4dgPLy8iZd1jpTamoqtm3bzhmaiVrBlk1bMXLoaGz6ZbO3S2l1CYnxEEJg+PDh3i6FmoiBh9qMzp07o7iouEnHdkrpCEVRsGrVqma9xsyZM3H06FHM+dfcFlRIRGf68IOPUFBQgBee+2uDfXa7Ha+8+Dd0S+uBZ598zrld0zTs2rEbO3fsws4du2C1Wj1ZcpOFhobiuusmo66uztulUBNxSkhqMzIyMpCTk4NNv2xG/0v7nbcbOSIyAoMHD8by5csxY8aMJr9G37598eKLL+Kxxx7DVWNHI6N7RmuUTtSmHNx/ENXV1UhLS0NQSFCD/Xa7HVKT2H/gINaszkRlZSVGjBqJmOgoVNfUoLamFgEmE5Z+txTR0dHYu3cvcrJz0KlzJxw6mIXS0jIUFxbhg/c/AADMm/cZ+g8YgGPHjmHxt4uxb98+l9fbtmuLT85x07dfXyxd+io0TeOg5TaAgYfajGuuuQa9evXC1NvuwPMvPIebb7vpvMfb7fYW3Sr68MMP49NPP8Vzz76A/877D6/PU7vyzcJv8ehfHgXgmLrhlddewRWDBzoDx5efL8CzTz/rnKE8KioK4eHheGfuu422d+edd2L27NkYN2Y8nnv+WTzz9LONHvf/Hv1/iI+Px+DBgxsEns8//QKTJk9ERXkFFn2zGOOvuRqpXTq3zhu+CF0zuqKmpgb79u1Djx49vF0OXQADD7UZMTEx2LZtGyZOnIgff/jxgoGnpqamyZMVnslgMOC1117DmDFjsGvHLvS6pFdLSyZqU47k5uGJWU9gxowZmDFjBh5++GHc/7v7nftHjBiBVatW4cYbb0SvXr0wYcIE9O/fHzabDUFBjp6gZ555Bhs2bMDSpUsBALNnz3aeXx92unfvjltuuQURERHo3bs3unTpgoSEBJfPa25uLlavXo0FCxbgb6++hldeftW5CLDRaMS9M6a7+8txQV0z0mEwGPDjjz8y8LQBDDzUpggh0LNnT8yfP/+Cx/Yf0B8LFy7Eq6++6lxmoqlGjRqF5ORk/PD9cgYeajd+XvczbDYbioqKsGLFigYzB69atQqxsbGYP38+5s+fj1mzZiE3NxeLFy/G4sWLsWzZMtx555249957sWDBAqSlpcFut0PTNGRkZCAuLq7Jn8WOHTvijjvuwB133IGioiK8++67eOyxxwAAa9esvWDg0TQN7819H+vX/YyDBw8iMjISMx/4PUZdNbLVLj8FhwRj7Lix+Pzzz/Hggw+2SpvkPn4deDZv9r87AwjIzs5GVlYWfv11M8yBgec87sphQ/D5Z59jyJAhmD17NmJiYpr1OmlpaVj902qMufoqCPCyFnnW4kVLsGTxEsx6ahaSkhI98prdemTggQf/gIULvsLChQuRkpKCjIyu2Ldvv/OYEydOAHB8Pnbu3ImPP/4Ys2fPxl//+ldMmzYN5eXlKC8vx+DBg13arq6uRnZ2NrKzs1tU25gxYzBs2DDceuut2LBhA5596jlcd+N1MBoa9uKeOH4CXy38GsuWLsPIkSNx3XXXYe/evXhg5oPo1KkT7pn+W/Ts3bNFdZytc+fOePvfb2Pt2rUwm82t0iY13Z49e5p8rF+upVVQUIDERM/8gCAiIiLvCQoKwp49e9CxY8fzHueXgQdwhJ6CggJvl+E3Jk2aBABYtGiRlyshIiI6LTo6+oJhB/DjwEOtKzU1FQCQlZXl5UqIiIiajxMPEhERkd9j4CEiIiK/x0taRERE5PfYw0NERER+j4GHiIiI/B4DDxEREfk9Bh4iIiLyeww8fu7ll1/GiBEjXLYtWrQIl112GUJCQpCSkoK//OUvqKmpce6vra3FzJkzERsbC4vFgttvvx1FRUXO/VVVVZgyZQoSEhJw7733QtM0bN++HUIIfPfddy6vtXTpUgghMHLkyAa1paenY+bMma37homIiBrBwOPH/vWvf+HJJ5902ZaZmYnrr78e119/PbZs2YK33noL8+bNcwkeM2bMwLJlyzB//nwsX74ce/fuxY033ujcP3v2bERGRmLZsmXIz8/Hp59+it69eyMmJgbr1q1zeb2lS5ciOTkZa9euRVVVlXN7QUEBDh48iLFjx7rp3RMREZ3GwOOH8vPzMWnSJDz66KPo2rWry763334bI0eOxKxZs5Ceno7x48fjpZdewscff4y6ujocPXoUH330Ed58800MHToUl19+OebNm4fVq1dj/fr1AIDKykp07doVvXv3RkpKCioqKiCEwKhRo7B27VqX11u2bBn+/Oc/Q1VVLF++3Lk9MzMTBoOh0Z4fIiKi1sbA44c2bdqEgIAAbN++HQMHDnTZ9+c//xmvv/66yzZFUWCz2VBZWekMLGcGka5duyIpKQmrV68GADz44IOYM2cOjEYjNm7ciLvuugsAMHr0aGzcuBF2ux0AkJOTg71792Ly5MkYNmwYli1b5mxz9erVGDhwIEJDQ1v/C0BERHQWg7cLoNY3adIk52KfZ+vXr5/Lc5vNhr///e+49NJLER0djby8PERHRyMwMNDluMTERBw5cgQA0KFDB+zduxeFhYWIi4tzHnPVVVfh5MmT2LZtGwYMGIBly5YhPT0dKSkpGDt2LP71r385j83MzMQNN9zQWm+ZiIjovNjD047Z7Xbceeed2LVrlzOMVFdXw2QyNTg2MDAQtbW1zueKoriEHQDo3LkzOnfu7OwlWrp0KcaNGwcAGDduHLKysnDgwAGUlZVh586dHL9DREQew8DTTlVWVmLSpEn4+uuvsWDBAlx22WUAALPZjLq6ugbH19bWIjg4+ILtjh49GmvXroXdbsfy5cudoaZXr15ITEzEypUrkZmZCYvFgssvv7x13xQREdE5MPC0QwUFBRg6dCjWr1+PZcuWYcKECc59ycnJKC4uhtVqdTknPz8fSUlJF2x79OjRWL9+PTZs2ICamhqXsUBjx47F6tWrkZmZiVGjRkFV1dZ7U0REROfBwNPOlJaWYtSoUSgsLERmZiaGDRvmsv/KK6+EruvIzMx0btu/fz+OHj3a4NjGjB49Gnl5efjss88wePBghISEOPeNHTsWmzdvxpo1azBmzJjWe1NEREQXwMDTzvzpT39CVlYWPv74Y8TExODYsWPOP5qmITExEbfddht+97vfYdWqVfjll19w6623YsSIERg0aNAF24+JiUHv3r3x/vvvO8fv1BszZgwOHjyIzZs3c/wOERF5FO/Sakc0TcNnn30Gq9WKUaNGNdifnZ2NlJQUzJkzBw899BCuv/56AMD48ePx5ptvNvl1Ro8eje3btzcINdHR0ejTpw+Ki4vRpUuXi3szREREzSCklNLbRRARERG5Ey9pERERkd9j4CEiIiK/x8BDREREfo+Bh4iIiPweAw8RERH5PQYeIiIi8nsMPEREROT3GHiIiIjI7zHwEBERkd9j4CEiIiK/x8BDREREfo+Bh4iIiPweAw8RERH5PQYeIiIi8nsMPEREROT3GHiIiIjI7zHwEBERkd9j4CEiIiK/x8BDREREfo+Bh4iIiPweAw8RERH5PQYeIiIi8nsMPEREROT3GHiIiIjI7zHwEBERkd9j4CEiIiK/x8BDREREfo+Bh4iIiPweAw8RERH5PQYeIiIi8nsMPEREROT3GHiIiIjI7zHwEBERkd9j4CEiIiK/x8BDREREfo+Bh4iIiPweAw8RERH5PQYeIiIi8nsGbxfgD6qqqjB79mzYbDZvl0JEROQ1ISEh+NOf/gSj0ejtUhpg4LlIVVVVCI1Pg6wugggIcstrKKr7/pmEqrqlXUVxX+ehQRVuateNNStuqtlN7QKA6qa23VgyFEj3NKxr7mkXADT3tS01e5tqFwB0u3u+HrrNfV9n3a67pV1Nd9P3MwC7m9qugoaff/4Zn332mc+FHiGldN9X1M85w05NCcy9p0ANTXTL64QmdXVLuwAQnpDgnnZjg93SLgB0iA9xS7vdEkLd0i4ApMW4p+aUcLNb2gWA+JAAt7QbZXZPyAaAEL3aLe2qFQVuaRcA5Ilct7VtKzjslnarD7unXQAoP3TULe2WHix0S7sAUHKwxC3tHj3hnu9nANhfZXVLu1uslfgehZh8/fU+F3o4hqeFzgw7apexbgs7REREbUUnBGEsYvD1woW45ZZbfGqoBwNPC5wddpTgWG+XRERE5BN8NfQw8DQTww4REdH5+WLoYeBpBoYdIiKipvG10MPA00QMO0RERM3jS6GHgacJGHaIiIhaxldCDwPPBTDsEBERXRxfCD0MPOfBsENERNQ6vB16GHjOgWGHiIiodXkz9DDwNIJhh4iIyD28FXoYeM7CsENERORe3gg9DDxnYNghIiLyDE+HHgaeUxh2iIiIPMuToYeBBww7RERE3uKp0NPuAw/DDhERkXd5IvS068DDsENEROQb3B162m3gYdghIiLyLe4MPe0y8DDsEBER+SZ3hZ52F3gYdoiIiHybO0KPRwPPr7/+ijvvvBMdO3aE2WxGly5dcO+99yI7O9t5zN133w0hxDn/xMfHt/j1q6qqMH78eIYdIiIiH9faocfQSnVd0D//+U889NBDGDlyJF555RUkJibiwIEDeO211zB//nysWLECffr0AQDEx8dj4cKFjbYTEBDQ4hpeeeUVrFmzBmraeIYdIiIiH9cJQRiJaCxcuBDvvPMOZsyY0eK2PBJ41q5diz/+8Y/4wx/+gP/5n/9xbh8xYgSuu+469OvXD7/97W+xadMmAIDJZMKgQYNavY7hw4fjxZdehl64CyIoBkJRW/01iIiIqHXYoGMXKmGEwMCBAy+qLY9c0nrttdcQHh6Ol156qcG+mJgY/P3vf8d1112HkydPurWOMWPGYNE3X0NW5EHLWQWpa259PSIiImoZG3QswQkUw4rV69ehf//+F9We23t4pJRYtmwZrr32WgQFBTV6zM0339xgm91ub/RYVVUhhGhxPRMnTsSib77GpGsnQ8tZBbXTCPb0EBER+ZCzw05rXPVxew9PUVERamtr0blz5yafk5OTA6PR2OifN95446Jrqg897OkhIiLyLe4IO4AHengMBsdLaFrTQ0VCQgK++eabRvclJye3Sl3s6SEiIvIt7go7gAcCT0REBCwWC3Jycs55zMmTJ2G1WhEREQHAcSfWpZde6u7SGHqIiIh8hDvDDuChQcvjxo3DypUrUVtb2+j+uXPnIjo6Gps3b/ZEOS54eYuIiMi73B12AA8Fnj//+c8oLi7Gk08+2WDfsWPH8Prrr6NHjx4XPQK7pRh6iIiIvMMTYQfw0Dw8gwYNwgsvvIAnn3wSe/bswbRp0xAdHY2dO3fitddeQ01NDX744Qfn8XV1dfj555/P2d4ll1xyzju+WoqXt4iIiDzLU2EH8OBMy0888QT69++Pf/zjH3jooYdQUlKC5ORkTJw4EbNmzXIZjHzs2DFcccUV52xry5Yt6Nu3b6vXyNBDRETkGZ4MOwAgpJTSra/QBn377beYdO1kiNAOTQ49Jkuk2+oJTerqtrbDExLc025ssFvaBYAO8SFuabdbQqhb2gWAtBj31JwSbnZLuwAQH9LyZVzOJ8rsvv9EhOjVbmlXrShwS7sAIE/kuq1tW8Fht7Rbfdg97QJA+aGjbmm39GChW9oFgJKDJW5p9+gJ93w/A8D+Kqtb2i2xNm1IiKfDDtAOV0tvCo7pISIicg9vhB2AgeecGHqIiIhal7fCDsDAc14MPURERK3Dm2EHYOC5IIYeIiKii+PtsAMw8DQJQw8REVHL+ELYARh4moyhh4iIqHl8JewADDzNwtBDRETUNL4UdgAGnmZj6CEiIjo/Xws7AANPizD0EBERNc4Xww7AwNNiDD1ERESu6sPOSUuAT4UdgEtLXLT6ZSjU8E5QIzu75TXM4bFuaRcAzKFh7mnX4p5lCQAgPNTklnbjQwPd0i4AxFjcU3Ok2X1f51CTe5aACDa67/9ZJume6fKV2gq3tAsAstI9yxIAgFZR7JZ2bcXuq7mmqMwt7VYXVrmlXUfb7lkCoqzSPd/PAHCizu6WdrdrVThpCcD333/vU2EHYOBpFfHx8Thx4gSMRqO3S2kWq9XxYQoIcN8vzdbGmj2DNXsGa/YM1uwZVqsViqJg7dq1Phd2AAaeVpGamgoAyMrK8nIlzdMW62bNnsGaPYM1ewZr9gxfr5ljeIiIiMjvMfAQERGR32PgISIiIr/HMTxERETk99jDQ0RERH6PgYeIiIj8HgMPERER+T0GHiIiIvJ7DDxERETk9xh4LtJ9992Hu++++4LHHT58GBMnTkRoaCgSEhLw1FNPQdM8t+BobW0tZs6cidjYWFgsFtx+++0oKio67zmHDh3CpEmTEB4ejoSEBNx3330oLy/3UMUtq7myshIzZsxATEwMwsLCMGnSJGRnZ3uoYoeW1H2mF198EUIIN1bYUEtqXrduHUaMGIGwsDAkJSVh+vTpKClx3xpLuq7jmWeeQVJSEoKDgzFhwoTz/tsWFxdj6tSpiIiIQGRkJGbOnInqaveseXQuza15165duOaaaxAVFYXY2FjcdNNNyM3N9WDFza/5TP/9738hhMDhw4fdW+RZmluzzWbD448/7jx++PDh2Lp1q+cKRvNrPnHiBKZOnYqYmBhER0fj1ltvRX5+vgcrdvXyyy9jxIgR5z3GFz6DLiS1iKZp8vHHH5cA5LRp0857rNVqlV27dpXXXHON3LFjh1y4cKGMjIyUTz/9tGeKlVLefffdskuXLnL16tVyw4YNsl+/fnLYsGHnrTk9PV1ef/31cvfu3TIzM1N27dpV3nTTTT5bs5RSXnXVVTIjI0OuWbNGbtu2TQ4dOlT27NlTaprmoapbVne9jRs3SoPBID390Wxuzfv27ZPBwcHygQcekHv27JGrV6+WvXr1kqNGjXJbjc8++6yMjo6W3377rdy6dascO3asTE9Pl3V1dY0eP2LECHnZZZfJTZs2yeXLl8tOnTrJu+66y231XWzNRUVFMj4+Xk6ZMkXu2LFD/vrrr3LYsGGye/fusqamxidrPtPhw4dlWFiYBCCzs7M9U+wpza35nnvukXFxcXLp0qVyz5498sYbb5Tx8fGyrKzMZ2sePny4HDJkiNyyZYvcvHmzHDRokLzssss8Vu+Z/vnPf0pFUeTw4cPPe5wvfAbPxMDTArt375aDBw+WMTExsmPHjhcMPJ988ok0mUyypKTEue3tt9+WoaGhsra21s3VSpmXlycVRZFLlixxbtu3b58EINetW9foOVu2bJEA5Pbt253b/vd//1daLBa31ytly2peuXKlFEK41Lxr1y7ZsWNHuXfvXrfXLGXL6q5XVVUl09PT5ahRozwaeFpS8xNPPCG7du0qdV13blu9erUEIA8dOtTqNdbV1UmLxSL/9a9/ObeVlpZKs9ksP/nkkwbHr1u3TgKQu3fvdm5btmyZFELIvLy8Vq+vNWqeO3eutFgssrq62rktNzdXApDLly/3yZrraZomr7zySuf3ricDT3NrzsrKkkII+e2337ocn5KS4rNf59LSUglAfvPNN85tX3/9tQQgi4uLPVKzlFIePXpUTpw4UQYHB8tu3bqdN/D4wmfwbLyk1QIrV65E9+7dsXPnTnTu3PmCx2dmZqJ///6IiIhwbhs1ahQqKio80o26du1aAMDIkSOd27p27YqkpCSsXr260XOio6OhKArmzJmDuro6FBYW4osvvsDAgQPdXm9La162bBl69+6N3r17O7f16NEDOTk5yMjIcG/Bp7Sk7np//OMf0bt3b9x5551urfFsLan5jjvuwIcffuhy6U1RHD9OSktLW73GrVu3orKyEqNHj3ZuCw8PR//+/RutMTMzEwkJCejevbtz24gRIyCEwJo1a1q9vtao+aqrrsLXX38Ns9ns3ObOr2ljmltzvZdeeglWqxWPP/64J8p00dyav//+e4SFhWH8+PEux2dnZ2PUqFE+WbPZbIbFYsGHH36IiooKVFZW4j//+Q8yMjIQHh7ukZoBYNOmTQgICMD27dsv+LvAFz6DZ2PgaYHf//73eOeddxAbG9uk4/Py8pCcnOyyLTExEQBw5MiRVq+vsdePjo5GYGBggxrO9fodOnTAm2++iQ8++ABBQUGIjY1FeXk55s2b5/Z6gZbVvG/fPqSlpeGtt95Cz549kZSUhFtuuQVHjx71RMkAWlY3ACxYsABLlizBnDlz3F1iAy2puVu3bhg0aJDLtldffRUJCQm45JJL3FIjgEY/R43V2NhnLiAgAFFRUR75zNXXADS95pSUFJfQCQCvvPIKzGYzhg0b5r5Cz9DcmgFg48aNeP311/Hxxx/DYDC4vcazNbfmffv2ITU1FQsWLMCAAQMQHx+PCRMmYM+ePR6pF2h+zSaTCR988AF++uknhIeHIzw8HD///DOWLFniDMWeMGnSJMyfP9+5Kvr5+MJn8GwMPGc5fPgwhBDn/NOcwaf1qqurYTKZXLbV/3Kpra11e82NvX59Ded6favViu3bt+OGG27A+vXrsWTJEmiahptvvrlVBlu7o+aKigqsWLEC8+bNw7///W/MmzcPeXl5GDlyZKt8nd1Vd35+Pu677z689957iIqKapU63V3z2R555BF8++23eOutt2A0Glv7LTgHOjb2OWqsxtZ4TxeruTWf7c0338Q//vEPvPrqq4iJiXFLjWdrbs0nT57E1KlT8eqrryI9Pd0jNZ6tuTVXVFTg4MGDeOGFF/Dyyy/jm2++gdFoxNChQ3HixAmfrFlKia1bt2Lw4MHIzMzEihUr0KlTJ0yePBkVFRUeqbm5fOEzeDbPx3Efl5SUdN6kf+ZlqaYym82oq6tz2Vb/Dx4cHNzs9s52oZoXL17c4PXrazjX68+ePRsrVqzAnj17oKoqACA9PR3p6elYtGgRrrvuOp+r2Wg0ora2Fl999ZXz32nBggVITEzEokWLcNNNN11Uze6oW0qJadOm4eabb8bVV1990fU1xh1f63o2mw333nsvPvroI8yZMweTJ0++6HobU3+Zp66uzuWSz7lqbOwzd77j3aG5NdeTUuLpp5/GX//6Vzz55JN44IEH3F5rvebW/OCDDyIjIwP33Xefx2o8W3NrNhqNqKiowLx585yXW+bNm4fk5GR8+OGH+Mtf/uJzNX/++ed48803kZubC4vFAgBYtGgROnXqhPfeew8PPfSQ22tuLl/4DJ6NgecsRqMR3bp1a9U2k5OTsWPHDpdt9bcTJiUlXXT7F6p5+/btKC4uhtVqRUBAgEsN53r9+nFH9WEHANLS0hAdHY0DBw74ZM0dOnRAUlKSSyiNi4tDVFRUq92a3tp15+bm4scff8TatWvx4YcfAgDsdjsAICQkBG+//TamTp3qUzXXq6ysxPXXX4/MzEzMmzevVQLludR3jefn56NLly4uNTZ2CS05ORlfffWVyzar1Yri4uJW+cw1RXNrBhwB8je/+Q0++eQTzJ492+O/yJpb83vvvQeTyYSQkBAAcPb+9uzZE0888QRmzZrlczV36NABBoPBZWyJ2WxGamqqx6awaG7NmZmZyMjIcIYdwPGf74yMjFb5eewOvvAZPBsvaXnAsGHDsHnzZpeuxxUrVsBisaBv375uf/0rr7wSuq4jMzPTuW3//v04evToOccGdOjQATt37oSU0rktPz8fxcXF6Nq1q0/WPHz4cOTk5KCgoMC5raCgAEVFRUhLS3N7zUDz605KSsKBAwewfft2bN26FVu3bsXzzz8PwDGw8dprr/W5mgHHD65rrrkGGzduxLJly9wadgCgT58+CA0NxapVq5zbysrKsHnz5kZrHDZsGPLy8nDw4EHntvpzhwwZ4tZa6zW3ZgC488478fnnn+OTTz7xyv/am1vzgQMHsHPnTuf37jvvvAMAWLJkCe6//36frHn48OGw2+349ddfndtqampw6NAhj/2caG7NHTp0wIEDB1wuBZ08eRJZWVke+XncEr7wGWzAK/eG+ZHhw4c3uC29rq5OFhQUOOdTqKmpkV26dJFXX3213LZtm/zqq69kZGSkfPbZZz1W5+233y47d+4sV65cKTdu3Cj79esnR4wYcc6at2/fLk0mk5w+fbrcs2ePXL9+vRw0aJDs27evtFqtPllzbW2t7NGjhxwyZIj89ddf5ebNm+WwYcNkt27dPHL7f0vrPtv777/v8Xl4mlvzM888I4UQ8tNPP5UFBQUufy40X0tLzZo1S0ZFRcmvv/5abtu2TY4dO1Z27dpVWq1WabfbZUFBgfOWbl3X5ZAhQ+SAAQPkxo0b5YoVK2RKSoq8++673VJba9Rc/+/+2muvNfiannmrui/VfLaVK1d6ZR6e5tZ81VVXye7du8vVq1fLXbt2ySlTpsjY2FhZWFjokzXn5+fLqKgoee2118pt27bJrVu3yokTJ8qkpCSPzh10pmnTprnclu6rn8EzMfBcpMYCT/2HfuXKlc5tBw4ckGPGjJGBgYEyMTFRPvXUUx6dDK+qqkpOnz5dhoeHy/DwcHnbbbfJoqKi89a8fv16OXLkSBkWFiYTExPl3XffLU+cOOHTNRcUFMjbbrtNhoaGSovFIm+44QaPz/nQkrrP5I3A09ya09PTJYBG/5zrfV0su90uH330URkTEyNDQkLkhAkTnL9Ys7OzJQD5/vvvO48/fvy4nDJligwODpbR0dFyxowZHp3Ar7k1jxkz5pxf0zPfly/VfDZvBZ7m1lxRUSFnzJgho6OjZVBQkBwzZozctWuXT9e8e/duOWnSJBkVFSVjYmLkDTfc4PGv85nODjy++hk8k5DyjGsWRERERH6IY3iIiIjI7zHwEBERkd9j4CEiIiK/x8BDREREfo+Bh4iIiPweAw8RERH5PQYeIiIi8nsMPEREROT3GHiIiIjI7zHwEBERkd9j4CEiIiK/9/8BvX2nZQ17H/QAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import seaborn as sns\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n", "\n", "fig = plt.figure(figsize=[8, 8])\n", "ax = plt.subplot(projection=ccrs.PlateCarree(central_longitude=180))\n", "ax.set_title('Coefficient of Efficiency')\n", "# latlon_range = [0, 360, -90, 90]\n", "latlon_range = [np.mod(-130, 360), np.mod(-70, 360), 30, 70]\n", "transform=ccrs.PlateCarree()\n", "\n", "ax.set_extent(latlon_range, crs=transform)\n", "lon_formatter = LongitudeFormatter(zero_direction_label=False)\n", "lat_formatter = LatitudeFormatter()\n", "ax.xaxis.set_major_formatter(lon_formatter)\n", "ax.yaxis.set_major_formatter(lat_formatter)\n", "\n", "lon_ticks=[60, 120, 180, 240, 300]\n", "lat_ticks=[-90, -45, 0, 45, 90]\n", "lon_ticks = np.array(lon_ticks)\n", "lat_ticks = np.array(lat_ticks)\n", "lon_min, lon_max, lat_min, lat_max = latlon_range\n", "mask_lon = (lon_ticks >= lon_min) & (lon_ticks <= lon_max)\n", "mask_lat = (lat_ticks >= lat_min) & (lat_ticks <= lat_max)\n", "ax.set_xticks(lon_ticks[mask_lon], crs=ccrs.PlateCarree())\n", "ax.set_yticks(lat_ticks[mask_lat], crs=ccrs.PlateCarree())\n", "\n", "levels = np.linspace(-1, 1, 21)\n", "cbar_labels = np.linspace(-1, 1, 11)\n", "cbar_title = 'CE'\n", "extend = 'both'\n", "cmap = 'RdBu_r'\n", "cbar_pad=0.1\n", "cbar_orientation='horizontal'\n", "cbar_aspect=10\n", "cbar_fraction=0.35\n", "cbar_shrink=0.8\n", "font_scale=1.5\n", "land_color=sns.xkcd_rgb['light grey']\n", "ocean_color=sns.xkcd_rgb['white']\n", "\n", "ax.add_feature(cfeature.LAND, facecolor=land_color, edgecolor=land_color)\n", "ax.add_feature(cfeature.OCEAN, facecolor=ocean_color, edgecolor=ocean_color)\n", "ax.coastlines()\n", "im = ax.contourf(ds['lon'].values, ds['lat'].values, ce, levels, transform=transform, cmap=cmap, extend=extend)\n", "cbar = fig.colorbar(\n", " im, ax=ax, orientation=cbar_orientation, pad=cbar_pad, aspect=cbar_aspect,\n", " fraction=cbar_fraction, shrink=cbar_shrink)\n", "cbar.set_ticks(cbar_labels)\n", "cbar.ax.set_title(cbar_title, x=-0.05, y=0.1)\n", "\n", "LMRt.showfig(fig)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 4 }