{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "# Gaussian Kernel Density Estimation (method 01)\n", "\n", "One of the methods used to calculate the density of particles is a 2-dimensional horizontal Gaussian Kernel Density Estimation (GKDE). This is a non-parametric (applicable to a non-Gaussian distribution) estimate of the probability density function of a field. \n", "\n", "In this example below, the GKDE is applied to a snapshot of particle positions, on all particles in the whole domain.\n", "\n", "We use the python SciPy stats gaussian\\_kde algorithm with a default kernel size (bandwidth), which is calculated following Scott's Rule [Scott, 1992]. This bandwidth selector calculates the optimal bandwidth value that avoids both over-smoothing and under-smoothing.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Author(s): Laura Gomez Navarro\n", "* Created on: 11/11/22\n", "* Last updated on:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Wed Jun 28 16:04:39 2023\n" ] } ], "source": [ "import time\n", "\n", "print(time.ctime(time.time()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 0. Imports and package versions" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from glob import glob\n", "import xarray as xr\n", "import numpy as np\n", "from scipy.stats import gaussian_kde\n", "\n", "import matplotlib.pyplot as plt\n", "import cartopy\n", "import cartopy.crs as ccrs\n", "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n", "import matplotlib.ticker as mticker\n", "import inspect\n", "import sys" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true, "tags": [] }, "outputs": [ { "data": { "application/json": { "Software versions": [ { "module": "Python", "version": "3.8.5 64bit [Clang 10.0.0 ]" }, { "module": "IPython", "version": "7.20.0" }, { "module": "OS", "version": "macOS 10.16 x86_64 i386 64bit" }, { "module": "numpy", "version": "1.19.2" }, { "module": "xarray", "version": "0.16.2" }, { "module": "scipy", "version": "1.6.0" } ] }, "text/html": [ "
SoftwareVersion
Python3.8.5 64bit [Clang 10.0.0 ]
IPython7.20.0
OSmacOS 10.16 x86_64 i386 64bit
numpy1.19.2
xarray0.16.2
scipy1.6.0
Wed Jun 28 16:04:46 2023 CEST
" ], "text/latex": [ "\\begin{tabular}{|l|l|}\\hline\n", "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n", "Python & 3.8.5 64bit [Clang 10.0.0 ] \\\\ \\hline\n", "IPython & 7.20.0 \\\\ \\hline\n", "OS & macOS 10.16 x86\\_64 i386 64bit \\\\ \\hline\n", "numpy & 1.19.2 \\\\ \\hline\n", "xarray & 0.16.2 \\\\ \\hline\n", "scipy & 1.6.0 \\\\ \\hline\n", "\\hline \\multicolumn{2}{|l|}{Wed Jun 28 16:04:46 2023 CEST} \\\\ \\hline\n", "\\end{tabular}\n" ], "text/plain": [ "Software versions\n", "Python 3.8.5 64bit [Clang 10.0.0 ]\n", "IPython 7.20.0\n", "OS macOS 10.16 x86_64 i386 64bit\n", "numpy 1.19.2\n", "xarray 0.16.2\n", "scipy 1.6.0\n", "Wed Jun 28 16:04:46 2023 CEST" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", "%load_ext version_information\n", "%version_information numpy, xarray, scipy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Loading the data" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "filedir = '../Simulations/'" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "ds = xr.open_dataset(filedir + 'toy_data_01.nc')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:     (obs: 121, traj: 144)\n",
       "Dimensions without coordinates: obs, traj\n",
       "Data variables:\n",
       "    trajectory  (traj, obs) float64 0.0 0.0 0.0 0.0 ... 143.0 143.0 143.0 143.0\n",
       "    time        (traj, obs) datetime64[ns] 2021-01-01T00:30:00 ... 2021-01-31...\n",
       "    lat         (traj, obs) float32 -39.17 -39.02 -38.9 ... -39.19 -39.21 -39.25\n",
       "    lon         (traj, obs) float32 22.5 22.55 22.61 22.66 ... 17.16 17.18 17.17\n",
       "    z           (traj, obs) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
       "    U           (traj, obs) float32 nan 2.598e-06 ... 1.845e-07 -7.776e-07\n",
       "    V           (traj, obs) float32 nan 5.796e-06 ... -1.211e-06 -2.669e-06\n",
       "Attributes:\n",
       "    feature_type:           trajectory\n",
       "    Conventions:            CF-1.6/CF-1.7\n",
       "    ncei_template_version:  NCEI_NetCDF_Trajectory_Template_v2.0\n",
       "    parcels_version:        2.3.1.dev20+g92f2fb90\n",
       "    parcels_mesh:           spherical
" ], "text/plain": [ "\n", "Dimensions: (obs: 121, traj: 144)\n", "Dimensions without coordinates: obs, traj\n", "Data variables:\n", " trajectory (traj, obs) float64 ...\n", " time (traj, obs) datetime64[ns] ...\n", " lat (traj, obs) float32 ...\n", " lon (traj, obs) float32 ...\n", " z (traj, obs) float32 ...\n", " U (traj, obs) float32 ...\n", " V (traj, obs) float32 ...\n", "Attributes:\n", " feature_type: trajectory\n", " Conventions: CF-1.6/CF-1.7\n", " ncei_template_version: NCEI_NetCDF_Trajectory_Template_v2.0\n", " parcels_version: 2.3.1.dev20+g92f2fb90\n", " parcels_mesh: spherical" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Calculating the GKDE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1. Loading the functions" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "sys.path.insert(0, filedir_root + \"Lagrangian_diags/Diagnostics/Functions/\")" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "import GKDE_01" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "import numpy as np\n", "from scipy.stats import gaussian_kde\n", "\n", "def rem_nans(ds):\n", " \"\"\"\n", " This renders lon and lat variables without nans for the last timestep.\n", " \"\"\"\n", " bad_indices = np.isnan(ds['lon'][:,-1]) | np.isnan(ds['lat'][:,-1])\n", " good_indices = ~bad_indices\n", " lon_end_nonans = ds['lon'][:,-1][good_indices]\n", " lat_end_nonans = ds['lat'][:,-1][good_indices]\n", " \n", " return lon_end_nonans, lat_end_nonans\n", "\n", "def kde_vals(lon_end_nonans, lat_end_nonans):\n", " \"\"\"\n", " Calculates the KDE values and sorts x (lon) and y (lat) in function of z (KDE)\n", " \n", " Inputs :\n", " - lon_end_nonans : Longitude values free of nans\n", " - lat_end_nonans : Latitude values free of nans\n", " \n", " Outputs: \n", " - x : longitude values\n", " - y : latitude values\n", " - z : KDE values\n", " \"\"\"\n", " x = lon_end_nonans.copy()\n", " y = lat_end_nonans.copy()\n", " \n", " xy = np.vstack([x, y])\n", " z = gaussian_kde(xy)(xy)\n", "\n", " idx = z.argsort()\n", " x, y, z = x[idx], y[idx], z[idx]\n", " \n", " return x, y, z\n", "\n", "def kde_parcels(ds, nsavedir=None):\n", " \"\"\"\n", " Inputs\n", " ds : xarray dataset of the parcels output netcdf file\n", " nsavedir: directory where to save the outputs. By default it is None, and will not save outputs. \n", " If outputs want to be saved specify directory, \n", " for example: '/data/oceanparcels/output_data/data_LauraGN/outputs_parcels/kde_calcs/'\n", " \n", " Outputs\n", " KDE for the las timestep fields of lon and lat\n", " \"\"\"\n", "\n", " # Remove nans:\n", " lon_end_nonans, lat_end_nonans = rem_nans(ds)\n", "\n", " kde_x, kde_y, kde_z = kde_vals(lon_end_nonans, lat_end_nonans)\n", " \n", " if nsavedir:\n", " savename = savedir + 'KDE_' + nfile.split('/')[-1].split('.nc')[0] + '.npz'\n", " np.savez(savename, kde_x=kde_x, kde_y=kde_y, kde_z=kde_z) \n", "\n", " return kde_x, kde_y, kde_z\n", "\n" ] } ], "source": [ "print(inspect.getsource(GKDE_01))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2. Applying the GKDE functions" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "kde_x, kde_y, kde_z = GKDE_01.kde_parcels(ds)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3. Plotting the outputs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1. Plotting cumulative distance evolution in time for the first 10 particle trajectories" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6wAAAHhCAYAAABqYjPfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3hUZd7G8e8zJT0hkAQSIPQuJYCIBaTau2vvvbP2dde2677uupbVXdeCbdFVLFiwUcQGKIIISFV6J5QkJEBImfa8fxxEUIQkTHJS7s91zUXOnJlz7slJwvzmacZai4iIiIiIiEht43E7gIiIiIiIiMi+qGAVERERERGRWkkFq4iIiIiIiNRKKlhFRERERESkVlLBKiIiIiIiIrWSClYRERERERGplXxuB6iI448/3ubn57sdo9YLBALExMS4HaPWOqxbcw7rnMYdj71N/745PPr744iP9bsdS0SiaMPmQka+PYXJM5cw6PAclmzY4XakqNLf+YZH17xh0nVvmBr6dZ89e/Yn1trjf3l/nShY8/PzmTVrltsxar3JkyczePBgt2PUasGSQibM3kxOh2YqVkXqkZXr8njslUnMmLeSs08exLczptOibRe3Y0Wd/s43PLrmDZOue8PU0K+7MSZ9X/erS7A0KP6Exgzrk83Hn34FgPHUic9sROQ3hMMRXnjnK8654zkGDB7OymWLeW70uHpZrIqIiDREKlilwTn3+j+zakMBDz73MSUlJRiP1+1IIlIFK9Zt4Zw7nmPK7OVM/2Ya9/zfEzRulo0xxu1oIiIiEiUqWKXBaZTWnHmLlrItlMgpNz3Ftu3FbkcSkUoIhyM89/YUzr3jec4+eTCTp39Pp2693I4lIiIi1UAFqzRIGRkZvPPhJxze/1Cee/srt+OISAUtW7OZs24bydfz1jJl4jv88eFR+GMT3I4lIiIi1UQFqzRo/3jiOd6YMJNN+dvcjiIi+xEKh3l2zFTO/8MLXHHl1Xw1cyHdDz8eY/TfmIiISH2mGWekQcvOzubKyy/lydFf8Pebz3A7jojsw9LVm7nz8XdolJzEd9/Npl3Hzm5HEhERkRqij6alwfvTPffx6fQfeOr1LygPBN2OIyK77NhZxqOjJnHBXS9w2XmnM3n6XBWrIiIiDYwKVqlRRUVFvPXWW4wePZrNmze7HQeAtIxMvvnqSxatLuCE655k+rwVbkcSadBC4TCjx33LsKv+ydaSCLO+ncZtD/wHrz/W7WgiIiJSw9QlWGrUK6+8QlFREQDr16/n9ttvx+dz/8ewU/dDmTB5Fm//7xmuv/kuvnjpdlKS4t2OJdKgWGuZMmspD704gYy0VD7+4F36DzzW7VgiIiLiIrWwSo2JRCK7i1WAYDBIaWmpi4n2ZozhnEtv5KQTjueFd6e5HUekQVm8ahOX3vMyDz4/nr/e9we+mrlIxaqIiIiohVVqjsfjoVOnTqxatQqA9PR0kpKSXE71a//3j3+S06MbF5zYj6yMRm7HEam3rLVMn7eSN8bPZMb8ldxx42Xcdv8/iYmJcTuaiIiI1BIqWKVGnXvuuSxatIhwOMwhhxyCMcbtSL/SunVrLj1zEA8+9zFP33uh23FE6p2ComLe+2wOb0z4jli/j0vOPYkXXhpFZvscLVMjIiIie1HBKjXK4/HQo0cPt2Mc0H0PPkGfw45i6qylHH1oJ7fjiNR51lq+nb+K18d/y5RZSzlx2JGM/PdDDD3hbPwJ6skgIiIi+6aCVWQf0lp25j9PPcXNI0YwYeTviY3xux1JpE7aum0n7302hzcnfIfX6+GiM4by/KjXyGpb+z+4EhEREfep75XIbzj1zPPo0asXz7/ztdtRRGpceSBIJBKp0nMjkQjfzl/JrY+MYcgVj7E0t4QXnn+OJas2cv8Tr6tYFRERkQpTC6vIfjw18r/07N6VC086jCaNEt2OI1ItVq7LY8K0RazdWMjqDXmsyc2ncHsJxkCztEZkZjQiKz2VzPQUMtOSyUxvRFZGIzLTG5GemoTX6yG/qJivZi1lyuylfD1nOemNk7jkvDN58c0JpGc0c/slioiISB2lglVkP1q1asWA/r2ZOmsZpw/LcTuOSFRZa3nns7n844VxnHvGCQw+fjCduvakU9futGzZkrKyMtavXc3a1StZs3IJa1cvZ93qlUyfN5dN+dvIzStie3EpjRslUloW4MicDgzu352/Pfgg3Q47Bo8v1u2XKCIiInWcClapF0pKSvD7/fj90R9revzxJ/DFxLEqWKVeKS4Ncd8z41i2bitTp82ge/fuv3pMQkICnbp0o1OXbsDJu+8Plm4jtHMrMSmZlBYXsWHNCjKbpRGf0hR/YpNaOfu3iIiI1E0VGsNqjHnNGLPRGLPdGLPUGHPVPh7zZ2OMNcYM3+O+DGPMl8aYDcaYy/e4v6Ux5l1jTL4xZpsxZoEx5rKovCKpk6y1FBcXV3rMnLWW999/n8cff5xHHnmEJUuWRD1bl559yc0rivpxRdyyKX87p454isx2vZg1Z+4+i9X98cc3Ij69Ld6YeJKaZNG59wAaNe9KTFKailURERGJqoq2sD4EXGmtLTfGdAEmG2O+t9bOBjDGtAfOAjb+4nm3AC8CHwCfGmPestaWAK8C84DWQDnQA8g86FcjdVJ5eTmjRo0iPz+fuLg4rrjiCpo0aVKh527ZsoUffviBcDgMwMcff0znzp2jli1Yup01K5awdVtp1I4p4rZZi9bSo1cfRo4c6XYUERERkf2qUAurtXaRtbb8p81dt/Z7POQp4C4g8Iuneve4+YCfPnrvB7xsrd1prQ1Za7+31k6o4muQOm7evHkUFBQQDocpKSlh8uTJFX6ux+PBWrt72+v1Ri3XO689T++e3Xj4Hw9xx+XHRu24Im4r2F5C67btD/xAEREREZdVeAyrMeYZ4DIgHvgeGL/r/rOBgLV2/D66gv0beBt4GPiLtXbnrvtnAE8bY/4DfGOtXXswL0LqNo/Hs9/t/cnIyKB///588803+P1+zjjjjKhkikQiXHrN73niD+cwrH+XSmWS+icQDLF+cyEpSfE0SorH74veByNVMe375Uyft5KiHWUU7iihcFsxCXGx5HRpSU7nbHp2bklKYtxvPj9/6zZSmjbG2gjG6GdbREREaq8KF6zW2huMMSOAI4DBQLkxJgn4O7DP5idr7UZgwD52nY3TInsf0MUYswC42lr7XeXiS32Qk5PDwoULWbNmDampqQwdOrRSzx8+fDhDhw7FGBO18XNL5kwmOSGWY47oFpXjSfQVl5QzddZS5i/PxUYiYDwQiWAMdGzdlIF9OtI0LaXKx49YmL14Pe9/OouJXy8kOSmB4pIytu/YSXJiHG1aZNA6qwmtmzehdfM02rZIp3XzNFKTE6L4Kn/B4+XVD6fz7FtTOO2EoRw6sBsZTbNIb9KYwsKtTJ8+jafHTGPB4pW0bNaYnC7Z9O7Sipyu2WRnNiE+1s+24lIWLlvPMRmp5P8wiYxDjq++vCIiIiIHyezZnbLCTzJmJPADzhjUbdbav+66fzVwlbX2s0ocKx14DDgGaGn3Eahnz572ySefrHTOhqa4uJikpCS3Y9QLq74ZzQtvf8nof/xqfjGpYcUl5Wwp2E5e4Q7yi4rZlL+dad8v57tFq+nTrS1de/bD44vDQxB/pBR/eAeLlq3lm7kryMpoxNGHduKo3h3ITEshMSGWpIQ4EuNi8Hp/3bJorWXJmi188MU8Ppg8n4SEBIYPP4Yhw44hMzNz92O2b93CpnUr2LRhNRs2rGd97mbWbtjCmtw8vB4PrZun0aZFOq2bN6FNVhqd2mTStV3mQbXUW2t5eNQkJn6zhH888k+Sk5N/8/c9Eixl87JvWfzjQub/uJK5i9exMa8Ia8Hv9zL0sK7ce82JZDRJZk1Mf9BESXWG/s43PLrmDZOue8PU0K/7kCFDZltrD/3l/VVd1saHM4Z1ENDSGHPDrvszgDHGmIettQ9X5EDW2nxjzGPApUAToOCXj4mJiWHw4MFVjNpwTJ48Wd+nKPn+6/G0aJq6130Tv17Ee5/PIXdLIYXbSthRUkYoHKFnp2zOOb4fJx/dnRi/VoqKlmDY8Mybn/HSu1+TkZZKRloq6amJpKcmcs6pw3n1lQto0bk/xvPr7rnWWkqLNvLlh/9j0mef89Qbk9laVMzO0jKKS8opKQsQF+MnMT6WxIQYkhLiiEQirMktoFFyAmedfjITPvknvXrlVKrVPlRWTO6aJSxbupTly5exfOUqpi1axdNvTaFo+06O6t2BgX06MqBPBzLTG1Xq+/H6+JlMnbWEaV+Mo2XH3gf+fT/mBOd7EQmB8YK1FG5aSXBnPsmpGcSntwMboe0+vn9Se+nvfMOja94w6bo3TLru+3bAd9fGmKbAUOBjoBQYDpwPXAD8Fdhz4cvvgNuA/U6gZIx5GGem4MU4Y2KvB5Zba39VrIq4ISu9EZvyd+zefv+Lefxr9Jc8+NcH6NI9h2aZWSQlJkK4nC8/n8RDDz/GirVbuFOTM0XFj6s2c9eTH9OqbUcWL11OixYtKvV8YwwJjZtz0qV/5MRL/kBg+xYi4XJ8sckESwoJhwLs2FZI0dY8thfmU1SYjw2Wkd0shUSvM79ci/YtK93F3BeXRKvOfWnVuS/D9rjfWsuS76cw4cN3+PLr7/j7ixNo2jiJAX06MrBvBw7r3pb4uJjfPO6yNZt5/JVJfPrB67Ts2Lty3wuP76dvCk2adwA67LFTxaqIiIjUbhVpDrI4BeVInFmF1wC3WGs/+OUDjTFhoNBaW3yAYyYAY4EsnCL4W+DUSuQWqVYtGwVYsXYzANtLQjz00idMnPQZffr0+dVjf3f+ZaRlteHOm6/d1ZIVrum49UYgGOLpN7/k9XEzefSxx7jsimsOelyyMR5iG/28apY/0VkyKTkLmh/UkSuTwdClz2C69BnMrUA4HGbmjK/54LVnePatqdz0tzfo1bklR+Z0ICHOTzAUJhAKEwyGCYbCTPrmB/78hxvoO/i0GkosIiIiUjscsGC11ubhdP09IGttmwo+bkRFHifilj7HXUVp4J+sXJfHzrIAGY0T6d37t1u2crp3YcHildhIKGoTPzU0i1Zu5s7H3qJ501Q+fnoEPYYNqbffS6/XyxFHDeKIowYRLC0ib90Svpg0jq9nzGFbOIYYnxcTKcHn85EUC9edczTnHH+Y27GlnolEIqxevRq/30/LlpXvUSAiIlITNOBOZB8SmmRz24jrePTlT7j+3EF4sITKtuOP3/e4w0ZpzQBDaXmQhP107ZRfC4Rg5Hvf8toHU3jogbu5+OJL8Cc2xuP1H/jJ9YA/PpXmnfpzUaf+XHTTr/dbaynNW0F8etuaDyf1lrWWV199ldzcXKy19O3bl+OOO87tWCIiIr+iBfhEfsPNt9zGgqXrmbVoNcaAx/Pbn++U7yzEWotHLRSVsmDZBk4b8SQr8sLMX7CIq268g9iUpg2mWK0IYwwJTTvsc3IpkaoqKipi/fr1BAIBgsEgM2fOdDuSiIjIPqlgFfkNKenNeei+m3l01CcAGO9vF6zbNy4DLLExzmO2FGxnw5YiinaUUB4IUZXlo+qz8kCIf748iSvue5kRl53BBx98QFZWltuxRBqM+Pj4vbYb8jIKIiJSu6lLcANjrWXp0qUEg0G6dOmCz6cfgf25+Kb7KSkP8cWnn1BelPub3TLjYn3E+H1c88BrbMwrIndLEXGxfkrLg5SUlpOd2YR3n7ieRsnx+3x+bWOtZe2mQpasyWPJ6k3kF+4gEAjRvmUavbu2onuH5sTG7LsVNByOsHbLDpasWE/+tp3YSBinXneK9kjE8ubE72iVlc7E/95HzgnXa+ycSA2Li4vjvPPOY+LEicTExHDaadUzoZe1lvLycmJjY/V7LiIiVaJqpYF5//33Wbx4MQDTp0/nqquu0puI/TDGw3V3PMh1dzy438fFNWnFC3+7gXJvI1p17s1RRw8jJi4RgFDpNq6+4CRueeRNThuSQ0lpgN5dW9G1Xe1rUbTW8uXMJfzn9c/ZVLCDPn0O5ZBDupLTJR0fERYsmMdfR45n+ZqNdGmXxSHtmxMX68fjMeQXFrNk9WZWrs+jWdOmdO/Zi8y0FII7txIJljjrgVqLx+vjtmvO46IrbiChaUf9/Im4pH379tx4443VdvzS0lJGjRpFQUEBKSkpXHnllWrJFRGRSlPB2sAsWLBgd/fULVu2UFRUROPGjV1OVXHWWvLy8vB4PKSnp7sdZ7eYpAzOu+Xxfe7zxTfiyRfeYMQNVzN17jri4mL51+jJ3HPtKZw26BB+anmsKdZaVqzLY3PBdvK3lbJl6w7ytm4nr2AbS1ZvImItv79wOJff9hAJTbL3eYyivA1M+/xDZs34ivLSnUQwdO2ew61/PJac/oP2+ab0p587FagiDcOsWbPYunUrkUiE7du38/XXX3P88ce7HUtEROoYFawNTGpqKkVFRc4EQR4PiYmJbkeqlHHjxjF//nystfTr149jjz3W7UgVkpyezctjJu7eXrx4MUMGDyYlwc+Qfh1rJEM4HGHS9MU88+YXbNsZoE2rlmRmZpKZlUWrzh3onRLHJakx9O2UgdfnJ3bXeqX7kprRgpPOu56Tzru+wudXoSoiIiIilaWCtYG5+OKLGT9+PIFAgGOOOYaYmLqzBEtZWRnff/89kUgEgBkzZjBs2DC83n3PnlpeXk5hYSFNmjSpda+zS5cu3HPvvXz66dgaKVgXr9rEXY+/iy82ngcf+TennnYGHo/mXBOR6tOvXz8WLlxIfn4+qampDBgwwO1IIiJSB6lgbWAaN27MhRde6HaMKvH5fHg8nt0F60/b+1JQUMCLL75IJBLB5/NxzTXX0KjRvtdQdUthwRYaJ1Zvq2MoHOa5MVMZ9f40/nTjBdz212fx/EaBLyISTXFxcVx33XUEg0H8fr96WYiISJWoiUXqDJ/Px9lnn01SUhLJycmcd955v/kGaMaMGZSVlREIBCgtLWX27Nk1nHb/goEyXnlpJCvXbuKxlz9jyepNUT/H0tWb+d0tz/LtgtV8OfZFbvvrMypWRaRGGWOIiYlRsSoiIlWmFlapUzp16sTtt99+wMclJCTg9XoJh8N4vd5frTnottKSnfTt0YnUJunENkrnoj++xBnDevP7C4eRlBB7UMcOhcM8//ZX/Hfs1/xpxKWM+MNfiE3OiFJyEREREZGao4JV6qWjjjqK9evXs27dOtq2bUu/fv0O6njbtm3D7/eTkJAQlXwpqWm8Ne7r3dt3/vFebhtxDcdd8wT3XHMSJwzsXqUWiWVrNnPn4++SFB/DN19OpHNvjRkTERERkbpLBavUSzExMVx88cUHfRxrLR9++CELFy7EWstJJ51E7969o5Bwb1kt2/DG2El8PvEDbrrhet6c+B1/ueEU2rWsWMtoKBzmpfe+4fl3pvDnP97G9TfdQkzSb8/yKyIiIiJSF6hgFdmPoqIiFi5cSCgUAmDixInVUrD+ZNjxpzFn/tE8cv8tnH3bSNJSk/B6vfi8HrxeDz6vh9TkBPp0bUWPjs1p0bQx5cEQ9zz5PsmNM5g1Zx7t2rWvtnwiIiIiIjVJBavIfni9Xqy1u7f9fn+1nzM+qTF/fvwVrrjmBlb/8B3xzboSwUsobAlHIuSuW81nH7/ByDFTyc0ronhnGTdfciz3/OsdLVUjIiIiIvWKClaR/UhJSWHYsGF8/vnnxMTEcNZZZ9XYubO79Ce7S/997jvnwsvZvm4229fMpqQsQMuex6pYFREREZF6RwWryAEcccQRHHHEETV+3sWLF7N06VLatm1Ljx499tpnPB4ate5Ho9b9mDx5Mp0zu9R4PhERERGR6qaCVaQWWrZsGe+99x7BYHD3hE89e/Z0O5aIiIiISI1SH0KRKAqUB3nyxhe4vu8feO/fH1f5OGvXriUYDAIQDAZZuXJltCKKiIiIiNQZamEViaL//fktPnl5MoHSAOuX5tKiY3P6n9in0sdp164dM2bMIBQK4ff76dSpUzWkFRERERGp3VSwikTRqoVrCZQGAAgHw6xfklulgrVt27acf/75rFixgtatW6tgFREREZEGSQWrSBSdct1xzJu8CI/HgzGGI0/rV+VjtWvXjnbt2kUxnYiIiIhI3aKCVSSKDj+5L/+e9jdWzl9Dr8GH0DQ73e1IIiIiIiJ1lgpWkSjYvCYPgGatM2jfqw3te7VxN5CIiIiISD2gWYKlztm5vYRwKOx2jN1euns0V3S9mSu63sx/73nd7TgiIiIiIvWGClapMyKRCA+e+zhnpl3O6Y0vZeG0xW5HoqyknLcf+5BAWZBAWZAxj35AeWm527FEREREROoFFaxSZ8ybvIip784gEo5QtrOc/zvncbcj4fV58Pq8e2x799oWEREREZGqU8Eqdcb6JbnYiN29XZS3zcU0Dn+Mnz+NvpnkxokkN07knjduxefX0HARERERkWjQO2upM3oP74HxGmzYKVoz2zQ96GOGQ2E+GjmJ9Us3cvzlQ+jQu22ljzHgjP4MOKP/QWcREREREZG9qYVV6oyWHZtzy7PX0DgzlXa9WvO3j/900Md8/s7/8eIfX+ODpyZw69H3sWn1ligkFRERERGRaFALq9Rqy+euYsyjH9K4aSMueeAcTrxqOCdeNTxqx//uk3mUlwQA8Hg9LJ21IiottyLiPmstxhi3Y4iIiMhBUMEqtdb2gh3cPujPlOwoxR/jY/WidTw86b6onuOw43PYsjaP8pIAkXCEzv06RPX4IlLzVq9ezVtvvUUgEGDIkCEMGDDA7UgiIiJSRSpY6yFrLeXl5cTGxtbp1oUNyzfBrvjBQIgls5ZH/RxXP3oxzTtksX5ZLsdeOphmrTOifg4RqVnvvvsuZWVlAEyZMoXu3buTmprqcioRERGpChWs9UxZWRmjRo0iPz+flJQUrrjiCpKTk92OVSVtDmlJbHwsgdIAXr+PI07tF/VzeL1eTr3huF/dP3/qD0x+axod+7Tj+CuG1unCX6ShCYVCe22Hw2GXkoiIiMjBUsFaz8yaNYuCggIikQjbt2/nq6++4sQTT3Q7VpXEJ8Xz7JxH+Py1qaSkJXPMJYNq5LzLv1/F3Sf+jfKSALEJsWzfWsy5d55WI+cWkYN34okn8sEHHwDQvXt3mjRp4nIiERERqSoVrPVcbW8ZDIVCTJ06lYKCAvr160ebNm322p+W1ZhzarhYXPTNEmzE+bq8pJzvJnyvglWkDunRowcdO3YkGAzW2R4mIiIi4tCyNvVMv379SE9PxxhDo0aNGDhwoNuR9mv8+PFMnz6dH374gddff52CggK3I9F9QBeMxyn04xJi6X9yH5cTiUhlxcXFqVgVERGpB9TCWs/ExsZy7bXXEgwG8fv9tb6Fdd26dbvHmxlj2Lx5M2lpaa5mat+rDQ9Puo+v3p1O+5y2DL/oaFfziIiIiIg0VCpY6yFjDDExMW7HqJBu3boxffp0wuEwxhiys7PdjgTAIUd25pAjO7sdQ0RERESkQVPBKq4aPHgwGRkZFBUV0a1bN3XhExERERGR3So0htUY85oxZqMxZrsxZqkx5qo99iUYY54xxuQbY7YZY6busS/DGPOlMWaDMebyPe5vaYx5d4/nLDDGXBbVVyZ1gjGG7t27M2DAAM3kKSIiIiIie6nopEsPAW2stSnAqcCDxpi+u/Y9DzQBuu7699Y9nncL8CLQGbjGGJOw6/5XgXVAayANuATYfBCvQ0SkQcvNzWXlypVac1RERETqlQp1CbbWLtpzc9etvTGmGKeAbWmt3b5r/+w9Huvd4+YDfpoBqB9wq7V2567t76sWX0REpkyZwrRp0zDG0LRpUy6//HI8Hk0CLyIiInVfhd/R7Or2WwIsBjYC44H+wBrggV3dexcYY363x9P+DVyz6zkv7lGgzgCeNsacZ4xpFY0XIiLSUE2fPp1gMEggEGDz5s3k5eW5HUlEREQkKipcsFprbwCSgYHAe0A50BLoDmwDmgM3Aa8YY7rues5Ga+0Aa22Wtfa5PQ53NvAVcB+wyhgz1xjTLxovqD4rKyvjk08+YezYsWzZssXtOCJSSyQlJe3+2lpLYmKii2lEREREosdYayv/JGNGAj/gdPV9GEiw1oZ27fsI+Mxa++8KHisdeAw4Bqdr8a8C9ezZ0z755JOVzlnfFBQUEAgEAGeyombNmu21zmpxcfFeb1ylYdB1b5j2vO6hUIiioiIikQgpKSnExcW5nE6qi37fGx5d84ZJ171haujXfciQIbOttYf+8v6qLmvjA9oDHx5UKsBam2+MeQy4FGfSpoJfPiYmJobBgwcf7KnqvIceemh3wfrT9yQzM3P3/smTJ+v71ADputc9c+bMYfbs2WRlZXHcccfh9/srfQxd94ZJ173h0TVvmHTdGyZd9307YJdgY0zTXWNNk4wxXmPMccD5wBfAVGAt8CdjjM8YcxQwGPjkAMd82BjTfddzkoHrgeXW2l8Vq/KzDh064PP58Hg8xMTEkJaW5nYkEamkNWvWMHHiRHJzc5k3bx6TJk1yO5KIiIhIrVWRFlaLU1COxClw1wC3WGs/ADDGnIazdM0fd+27xFq7+ADHTADGAllAKfAtzmzDsh9nnnkmc+fOpbS0lJycnCq1yoiIuwoKCvhp5EMoFGLzZq3oJSIiIvJbDliwWmvzgEH72b8IOKIyJ7XWjqjM48Xh9Xrp27fvgR8oIrXWTz0lfhp/3q+f5psTEakOwWCQsrIykpKS9przQ0TqlqqOYRURkSpISUnh+uuvZ+XKlWRkZNCiRQu3I4mI1Dtr165l9OjRhMNhsrOzufjii7U+tUgdpd9cEZEalpKSQk5OjopVqTBrLZ999hl5eXl89NFHhMNhtyOJ1GoTJ04kEAgQDofJzc1lxYoVbkcSkSpSwSoiIlLLzZ8/n5kzZxIKhZg/fz7Tpk1zO5JIrebz7d2J0Ov1upRERA6WClYREZFarrCwkGAwCDiTdRUUaFJ9kf05+eSTd49d7datG23btq30MVauXMmMGTPIz8+vhoQiUlEawyoiIlLLde/enenTp2OMwe/3c+ihv1pXXUT20LRpU0fLiSwAACAASURBVG6//XastVWacGnu3LmMHz+eSCTCF198wbXXXqvlBEVcohZWERGpl1bMW80Xr39Ffu5Wt6MctPT0dG688UZSU1O5/vrryc7OdjuSSJ1Q1dmB582bRzAYJBwOY63VGFgRF6mFVaQBCAaCLJu9ksaZqWS1beZ2HJFq982H3/H3C/6Fx+PB4/Xw/LzHaNoqw+1YeynYWMi8LxfSqltLOuS0JRwOs3l1Ho0zU4lPjPvV41NSUoiLi6Nx48YupJXqVlRURFFREVlZWYRCIRISErQUi4tatmzJ+vXrCYVCGGPIzMx0O5JIg6WCVaSeC5QHufnIe9iwbCPhcITbX7yOoecPdDuWVMDS9Xn879NZpKUkcs1Jh5MYF+N2pDrjg6cmUF4SACAmzs/0j2Zz2o3Hu5zqZ1vW5XNtrzsIh8JEIpbbX7yO0X97j02rtuDze3nsy7/QIafyY+6kblq6dCnvvPMO4IxR9ng8pKSkcNVVV7mcrOEaPHgwHo+HDRs20Lt3b1q1auV2JJEGSwWrSD23YOoPbFi+kdLiMgBevv8tFax1QFFxKVf+cww7ywL4fV6W5+bz9Igz3Y5VZ7Tp3oqFXy8mUBbE4/XQomOW25H28u3HswmUBQiUORMpvXz/m2xalUckHKEcePm+N3nwoz+5G1JqzJQpU3ZPqgUQDofZtm0b3377rVpZqyAYDBKJRIiNja3yMbxeL0OGDIliKhGpKhWsIvVco4wUIuEIAMZjaNIs1eVEUhFrtxTt/joYCrNo9SYX09Q9V/ztfMp2lrFk5gqOvXwwhx7by+1Ie2nRMQuP15lGIibOT0Fu0e7fU4DYhKq/0Za6JyUlhU2bNhGJ/PwzYIzB4/FgrXUxWd3zww8/MHbsWCKRCIcddhjHHXec25FE5CBp0iWReq5DTlsufeBcUpum0L5XG+763wi3I0kFtG+eRozfi8/jIS7Gx8Du6h5aGbHxsdz63HWM/P5Rzvz9SW7H+ZU+w3ty5UMX0j6nDUMuGECgtPznnQbO/9MZ7oWTGnfSSSeRnZ1NYmIiiYmJAGRkZNC/f3+Xk9UdhYWFfP3117z//vuEQiEikQizZs2isLDQ7WgicpDUwirSAJx9+6mcffupbseQSkiMi+HNuy/i429/oHFSAicf3s3tSBJlp990AqffdAIA6xZvYPmcVVhryWrXjHY9W7ucTmpSUlISl1122e7tcDiM1+t1L1Ads3PnTp5//nkCgcBerdRQ9VmCRaT2UMEqIlJLZaQmcflxh7kdQ2rAw5PuZ8KLnxMOhTnhyqF4POoA1ZCpWK2cDRs2YK3dq1j1eDwceeSRpKZqGIxIXaeCVaLqp0834+J+vSSDiIjsW1xCLGf8/kS3Y4jUSU2bNt1drPp8Prp06cLpp5+uwl+knlDBKlEze/ZsJkyYgLWWo48+mkGDBrkdSUREROq51NRULr74YmbMmEGTJk04+uijVayK1CMqWCUqrLVMmDCBcDgMwFdffcVhhx1GfHy8y8lERGpeYXEpf3j+Y1ZsLOCUw7txy5kDNZZOpBplZ2eTnZ3tdgwRqQYaJCNR88s3Y3pzJiIN1aNjvmTeilyKikt5Z+p8vl64yu1IIiIidZIKVokKY8zu8SIej4fhw4drHKuINFh5RTsJ7RpTZ7EU7ChxOZGIiEjdpC7BEjWHHHIIXbt2BdAMlyLSoF1z0uHc/Mz7eD0ekhNiGdqrg9uRRERE6iQVrBJVKlRFRKBf52w+eOByNhRso3N2U+Jj/G5HEhERqZNUsIqIiFSDjNQkMlKT3I4hIiJSp6k5TERERERERGolFawiItG0cyds2gTW7nv/N9/AqafCjTfCtm01m01ERESkjlGXYBGRaPn0Uzj9dAiH4fjj4b33YM9x3bm5cOyxTlEbEwOrVsH48e7lFZFKiUQss5etB6Bvx5Z4PFq+TUSkuqmFVUTqpo8/htatoVMnmDHD7TSOG2+EkhIoL4dPPoGpU/fev2wZeL3O14EAzJlT8xnrEWstkZ0vEym4kEjxM1gbcTuS1HP3jJrALc9+wC3PfsAfXxzndhwRkQZBBauI1D3FxXDOObB2rVMEnnyy24kce649XFYG99+/9/7evSE+3nlcYqLzGtavhyOPhKwseOihms1b15WNg+InIPgdFD+HLXnd7URSj5UGgnw2Zyml5UFKy4N8OW85xaXlbscSEan31CVYROqeHTsgskdrWlGRs+32skojR8JRR/28PX06FBZC48bOdkoKzJsHb70FmZlw1llwwgkwc6bTjfjBB2HIEDj8cHfy1zE2tBhs6a6tUggtcjWP1G+xPh/xsX6KSwMAxMX4idNyRSIi1U4Fq4jUPVlZzsRF48c7kxvddFPNFKtlZfDKK06X30svhUaNYOFCGDcOcnKc8akJCU63YHCy7dnqCtCsGfz+9z9v5+Y6xSo4r2Hz5up/HfWEiTsWu/N/e2yf4mIaqe88HsOzv/8dD47+DAvcff4wfN7a2VEtHApTsqOUpNREjNE4WxGp21SwitRTJaFyHvnhQ5bt2MgZLQ/jrNb1rNXurbdg1iyIjYWePaN33LffhgcegJYtYdQopzj+ySmnwLRpTiH67LMwZgwccYRTwMbEwFNP7V2wxsTA/PnQv/9vn+8vf4GLLwafD1q0gGOOid5rqeeMvyekvQPBWeDvhfEf4nakOisUCbM1UExabDJeUzuLsNrgkDaZvHHPRW7H2K/Vi9Zx++A/U7K9hK6Hd+LhT+/DX49aggsLC3n33XfZsWMHgwYNok+fPm5HEpFqpoJVpJ761+LxfLZpAYFIiCeXTqBdclP6NGnndqzoMQb69YvuMVescFpOS0th8WI47zyYMsXZF4nA55//vFzNmjUwdiyEQj/fXn/daUHNz//5mD91B/4tv/sd9O3rjGU99NBft8jKfhl/J/B3qvHzWhvE7nwNIrmY+HMw/o41niFaNpRs5coZIykOlZEZl8p/j7ieFH+827GkkuatzGXm4rV898IX7Ni6A2th2ZyVfPXutww9f4Db8aLmnXfeYePGjVhrGT9+PBs2bKBly5bk5OSoNVmkntLHqCL11KriLQQiIWfDwtqdBe4GqgvWrXNaOsHpprtixc/7PB7o2tXZ7/FAUhIMHPjzrL8JCc72mDHQpQukpcEjjzizGB9ImzYwYED9K1Y//BCuvRbeecftJFFnt//ZmfCp5BXs1nOw4Ty3I1XZSyu+pDCwk0AkxMayQj5aP8vtSFJJc1ds4Pp/vcvIj6Yzp2kc5a2b7N5X35be2bFjB3bXB4fhcJg5c+YwYcIEPvvsM5eTiUh1UcEqUk+d2/pI4jx+ErwxxHh9DMjo7Hak2q9/f2cypORkpwC99da993/+OVx2GZx/vtM1eMgQePVVOO00uPdeuPtuZ5zr2LFOK+tNN7nyMmqFceOc79Pzzzut1mPGuJ0ousqnAWW7NjwQWuxmmoMS6/Hh2dUy5cHg83hdTiSVNX3RGsqCISwQ8RjolonHY+g+oAsDztzPkIQ6aNCgQfh8Pnw+3+4W1WAwyNKlS11OJiLVRV2CReqp4Vk9yE5MY/3OtRyWuIQk+wXWnoQx9WcsU9TFxztro372mTN29ZdjTzMz4YUX9r7vzDOdm7XOrL+ffOJ0H/7Tn+C++2oue20zefLPY3lLSpzv6TnnuBopqmKPgtKPcYrWCPi6uJ2oyq7qMJQ5W1exeucWujVqyWkto9zVXqpd97aZxMX4KAuEiIvx8fu7fsfJb3cmNj7W7WhR17dvX1q1asWmTZv48MMPCYVC+Hw+2rZt63Y0EakmKlhF6rFOyc3oGLgaStZhMVA2EdN4pNuxarekJDj99Mo/b80aZ9bisl2tbn/7W8MuWIcPh6efdsYDJyQ4y/fUIyblAayvE4RzMfFnY7wZbkeqsrTYZN4aeAthG9GES3XUwB7t+PNFx/DFvBX065zNmQN61OvxnBkZGWRkZNCkSRNmz55Neno6/fc3uZ2I1GkqWEXqs8gmCK1jd9fF8slYa+v1GxnXpKTsvZ2a6k6O2uK445yu0Z984nSdPqV+LTljjB+TeJnbMaJKxWrddly/LhzXr+629FdFixYtaNGihdsxRKSaqWAVqc886WBiwJYDHvC2VrEKsGQJFBU5s/J6ozRer0kT+N//4OabnRbF//wHtmyBpk2jc/y66LjjnJuIy77asphPcufSI7UVZ7c+HI+KcxGROkN/sUXqMWNiMU1eh9ihEHcCpsnLbkdy3xNPQO/eTpfVY491xptGy9lnQ24uXHABnHEGtGoF//hH9I4vIpU2t3A1d899g0mb5vP00k8YvfprtyOJiEglqGAVqeeMvxOexs/iSX0c4810O477HnjAGVdZXAzffguLFh3c8T780ClKf/jB2d6xw9kuL3du990HweDB5xaRKllUtJ6IdT6YKosEmVWw0uVEIiJSGSpYRaRhSU//+etw2OnKW1VPPQXnnecsZ9OnDyxb5qzTume365/WbZX9W74cjj8eBg+G775zO43UI/3S2uE1HgwQ5/EzLLO725FERKQSNIZVpIGw5ZOxJW+Brysm6fqGu7zN2LFOkVlY6LSEVnXCjsWLnaVrSkud7fJyGDECJk6El16Ca691CtVRo6I3TrY+Gz4c1q1zumgPHw6bNjnLDIkcpE4pzRnZ/2qmbVlMl0YtGNi0q9uRRESkEir0sb8x5jVjzEZjzHZjzFJjzFV77DvHGPOjMWaHMeYHY8zpe+zLMMZ8aYzZYIy5fI/7Wxpj3jXG5BtjthljFhhjLovqKxOR3WzwB2zh76H8c9j5InbHo25Hck+PHk434NxcuOSSqh/nuOOcbsV7mj8f8vPh3HNh506ne/BZZx1c3obA2p+LVYBAwPk+ikRJt0YtubrjcBWrIiJ1UEX7qT0EtLHWpgCnAg8aY/oaY1oArwG3ASnAncDrxpifpsW8BXgR6AxcY4xJ2HX/q8A6oDWQBlwCbI7C6xGRfQn+uEc31TIIzHE1Tr2Qm7v3ts/njFVt0QKaNXOKV6kYY5wiPynJufXqVfWWbxEREalXKlSwWmsXWWvLf9rcdWsPtASKrLUTrGMcsHPXPgDvHjcf8NM75n7Ay9bandbakLX2e2vthOi8JBH5lZjDcH79fEA8xJ/scqB64IorIDHR6bbapInT4lpS4rQOFhbCbbe5nbBuee015/bSSzBlisb9NjARG+GJH8dx+pRH+ev8dwhEQjVy3uJgGbMLVpJftr1GziciIpVX4TGsxphngMuAeOB7YDxQCvxojDkVGAecApQDPzUt/Bt4G3gY+Iu1dueu+2cATxtj/gN8Y61de/AvRUR+i/FlQ9p7UPY5+Npj4oa6HanuGznSaRUsKXGWx/noI5g82dlnDPgb6BjhqvJ44LTT3E4hLhm/4XvGrptJWSRIQfkOmic05qoOw6r1nFvKtnHRtKcI2hARa3mm35UckppdreesrZauz+OliTNplBDLDaceRWqSxo+LSO1R4YLVWnuDMWYEcAQwGCi31oaNMf8DXgfigABw9k+FqbV2IzBgH4c7G7gLuA/oYoxZAFxtrdXUkCLVxPjaQVI7t2PUH8bA0D0K/9NOc7bHjXO6BP/73+5lE6ljNpYW7W5VLY+EWFdSUO3nnLRxPsWhMkI2DMCrq77iH70vqPbz1jY7Ssu58p9j2FkWwOf1sGxDPqPuPM/tWCIiuxlrbeWfZMxI4Iddt7eA44A5QF/gQ+AEa+3cCh4rHXgMOAZoafcRqGfPnvbJJ5+sdM6Gpri4mKSkJLdjSA3Tda9lrN17WZtqouveMNXX6x6IhFi9Mw8Ai6V1QgZx3urtpbA9WMrG0kIsFoOhcUwiTeMaVes5q6K6r3l5MMTqTYVEdr398hhD5+yMajufVEx9/V2X/Wvo133IkCGzrbWH/vL+qi5r48MZpxoDTLXWztp1/3fGmG+B4UCFClZrbb4x5jHgUqAJ8KuPVWNiYhg8eHAVozYckydP1vepAdJ1r2MiEXj7bdiwAbp3h8xM6Nmz0ofRdW+Y6vN1zy/fwY/bNtAhuRlZ8Y2r/XwRG+GxHz/mi00LOaRRS/7a63gSfbEVf37Ecs+oCUyavYT0RomMvPks2mYexLrOv6G6r3l5MMQZf36ZrTtK8HgMg3q249p6+jNWl9Tn33X5bbru+3bAgnXXjL9DgY9xxqwOB84HLgCKgD8aY3KstXONMb2BgcAzBzjmwzgzBS/GGRN7PbDcWlv9fYBERH5iLSxb5kyeVFOz0t5+O7zwgrN+ayTiTNp05ZXwn//UzPnloO0IlvLyyskUB8u4qO3RZCemuR2pXkiPTWZg0y41dj6P8fCHbqfyh26nVun50xat4qsFK7EW8ot28tAbn/P8rWdHOWX1i/X7GH33hYz/9keSE2I58TAt/SMitUtFpmG0OAXleqAQp/vuLdbaD6y1U4C/AO8YY3YA7wJ/t9ZOOsAxE4CxOAXvSpzlbar2P4aISFVYCxdeCDk50KEDPPVUzZx3zBhnjdaf1hwtLYVnn4Xy8v0/T2qNW2e/wpurv+H99d9xxYxnKQsH3I4kLgiGI7u/tkAwFHYvzEFqnBTPhcP6cOoRh+DzaoZuEaldDtjCaq3NAwbtZ/9TQKXe6VlrR1Tm8SIiUbdyJbz/vlMwAvzxj3DTTdV/3kMPhU8+2btAjY/XrMJ1yI/bNhDcNVFPIBJiU2kRbZKaHuBZsmnrDt6eOo+k+FjOG5JDfEzd/pkf2L0tXVo1ZeGqTfh9Xu44e7DbkURE6qWqjmEVEanbEhN/buUEqKlJDl57De69F2bNghUrnGJ11CitO1qH9E1rx/dbVxG2EZJ8cTRPiP64xfqmNBDkwn+MZltxGT6vh1lL1/H0iDPdjnVQ/D4vL9x6Nvnbd5KSEEesX2+pRKTyrLVMnDiRRYsW0a5dO8rLy4mNrfh4+oZAf11FpO6aM8dZD7VdO7jtNoiJqfhzMzPhySfhzjudYnXMmOrLuafkZC15U8c92vsixq6bSXGojNOzDyPGo/9KD2RD3jbKA856p4FQmJmL1/Lg6M+4YGhv2mXV3THAxhgyGjXcGT1F5OAtWrSI77//nmAwSHl5OZ9//jknnnii27FqFf0vKyJ107p1cPTRznjQ+HhYvhxefPHXj7PWmZU3NxfOOQeaN/953zXXODeRSoj1+jmvzVEVfvzS9Xls3VFCnw4tiGmgrXAt0hsR4/dSGgiCdWbYHTttAZNmL2Xcg1eQnBDndkQREVcUFxcT2aPH1/bt211MUzupD5qI1E3z5oHX63xdWgpffLHvx915J1xxhTNGtVcvKCysuYzS4L3++Rwue/RN7njuIy555M06PTHPwYiP9fPqXRdw3pDe+HweLM5nSdZa1m4pcjueVNCaH9dz7ykPcd9p/2D9so1uxxGpF7p3705cXByxsbEYYzjqqIp/INpQNMyPekWk7uvb13nHa4zTwnrSSft+3JtvOq2wALGxMHs2DB9eczmlQRs16TvKAiEA1ucV8cOazfRq3/wAz6qfWqQ34s6zB5NXVMxXC1YRCofx+zy0qYa1S6vL+rwiPv9+OS0zGjE0pwPGmKge/8u5y/nLq85CC3++6BiG9u4Y1eMfjHA4zO2D7md7wQ4whqXfreDNDc9H/Xsg0tAkJSUxYsQINm3axPLly1m6dCkTJkygS5cuDBw4UL9jqGAVkboqKwtmzoT//Q/atHHWMt2Xfv0gP9+ZlTcUciY7+vFHuPRSSEmp0cjS8DRNTaJwRykRawlHLGkpCW5Hct3fLj+Bd79eQNGOEk47qjuJcZUYe+6ivKJizv/7aMoDIfw+LyuOLeCakw6P2vGD4TB3/3c85UGnFf7u/07gqyfaRe34B6tkeyk7t5VgdzWPF+VtJ1AWIDZek8OIHKzY2Fhat27NwoULmTdvHsFgkPz8fBo1akSvXr3cjuc6FawiUnd16QJ///v+H/Pqq3DPPbB6tTOO9a9/dVpmR46EBQs0O69Uq0evOZl7R00kb9tOrj/5CFpmpLodyXV+n5fzBue4HaPS5q7IBSyhSIRQIMKkWUuiWrBGIpbQHmu7hn+xXZNe+csYPnr2E1p0zOT+t+8gLasxSamJdOnfkeXfrwKg2xGdVazWYmvXrmXVqlW0atWKtm3buh1HKigUChEMBgEIBoPk5eW5nKh2UMEqIu7bvv3npV2uuMJZciZakpKcWXkjEfD5nGIVnHVYN26EFi2idy6RX2ie1oj/3nGu2zEkCjq2SCcccf5+xPq9Ue/aHev3cemx/Xjts9kAXDisD/GxNb9W7bwpi3jnnx9StrOcHVuL+de1z/F/H/4RYwwPf3o/U9+ejsdjGHhW9Ip1ia7Vq1czevRoQqEQfr+fs88+m44da0/3cvlt8fHx+P1+jDFYaznkkEPcjlQrqGAVEXdZCwMHwtKlzvbrr8P06dE/j8cDnTvDsmXOOZOTISMj+ucRkXqpTWYTnrj+NMZMmUvbzCZcfWL0C7abTjuK3w3sgbXQPM2dIQtFW7bvHjMXCUcoyN26e19MrJ/hFx3tSi6puGXLlhEKOWPng8EgixcvVsFaR8TExHD11VezYcMGsrOzSUuru8t+RZMKVhFx19atsHgxBALO9syZzqy/8fHRP9cXX8C990JZGfzlL5Vbt1VEGrz+XVrRv0uraj1HVhN3x9YfdkIOTTJTKdy8jXAozKUPqIdAXZOdnY3f7ycYDOL3+2nVqnp/ZiW6MjIyyNAH6ntRwSoi7kpNdVo6N250WkGzsyGumtZkzMqCl16qnmOLiNQD8UnxPD//nyydvZKmrdJpmp3udiSppC5dunDyySezZMkS2rdvT8+ePd2OJHJQVLCKiLu8XvjmG6fF0+NxJkWK5hTu1sKqVZCQAJmZv/2YLVuc4jlWk4iISMMWExdD96O6uB1DDkLPnj1VqEq9oekxRcRdBQXw1FPOEjMPPADNoziRibXO8jXduztL34wc+evHBIMwZAi0bg1NmzrL3oiIiIhIraCCVUSqbsUKmDzZGXNaVcOGwb/+5RSthx/urJUaLatWwTvvOPnKy+HOO3/9mHHjYPZsZ//27XDLLdE7v4iIiIgcFHUJboCKi4sZM2YM+fn55OTkcMwxx+yeEVCkwt5+Gy67zOnSm5kJc+Y4S8hURjgM8+f/vNRMfr7TNTdarazx8c5yNj/Z13I5/j2WjTBGEzGJiIiI1CJqYW2AJkyYwIYNGygtLWXWrFksX77c7UhSCaFImDFrpvPvxeNZWbzZvSAPPAAlJbBjB2zaBJ99VvljeL1w2GHOuFG/35kUqVmz6GXMyoInnnAK1WbNnCL7l044AU480SlWMzPh6aejd34RkTpg6ewVLPz6RyJ7fsAnIlJLqIW1ASouLt7rP6WSkpKonyMcDmOtxefTj1i0PfbjR4zb8D3lkSBj133Hvb6B7gRp2dJZjiYcdm5ZWVU7zqefOt2By8vhxhudIjaarr/euf0WjwfeegtGjwb9vIpIA/PcHa/w0chPMR5DzpBD+Ov7d6nXlYjUKnp31gANHjyYN954A2MMiYmJdOkS3ZkA582bx0cffYS1lmHDhnHkkUdG9fgN3fS8pZRHgoDTKFgeDroT5L//hfPOg5UrnXGf/fsf+DkrVsDSpc5Y1caNnfuSk+FPf6rerBWhYlVEGhhrLe/9ezyRsPMh9pxP55O/YSsZLdNcTiYi8jO9Q2uA2rZty4gRI9i2bRvNmjXDv+cYvoNkreWjjz4iHA4D8MUXX9CnTx/iqmtdzQbosPQOTMydS3kkhLWWWLcKrebNYerUij9+0iQ44wynMIyLc8auRrP7r8gBWGv5aOQkfvhmCYPOOZIjTjnU7UgirjLGkJKWTNGWbbu3ExsluJxKRGRvKlgbqOTkZJKTk2vkXPanCXUkKu7qdhqtEzPILS3k9Jb9yJ2zxO1IFfPII86YV3CWkhk7Fq67zt1M0qC8/58J/Pfu1ykrKefrsd/yt3F302vQIW7HEnHV38ffzSOXPUV5SYCb/nMlCcnxbkcSEdmLClaJKmMMJ598Mh9//DHWWgYPHkx8vP7ziyafx8tFbX8et5pLHSlY27eHr76CQMAZN5qd7dxfXg4vvww7dzqzDjdp4mZKqce+/3wBZSXlAIQCYRZ/u1wFqzR4Hfu044X5j7sdQ0TkN6lglajLycmhe/fuWGuj2t1Y6rhHH3WWrJk71ylMTzzRuf+MM5y1XCMRZ4bexYv3XmpmX8JhZ+xsZqYzBlakAgb8rj9zPl9AeUk5vhgvOUO7ux1JREREDkAFq1QLzQ4sv5KS4nQD/qVPP4VQyPl60yZYu9Zpjf0tJSVw5JGwbJkzo/Cnn1Zswidp8I69ZDBJqYksm72Sfif0pvOh+/k5ExERkVpBVYWIuKtHD1i40Cla4+KcyZz2Z+xYWL5893hYe8c52Al3YRKvwpiYGggsddmRp/bjyFP7uR1DRESkxvw0n0xdXbJKBavUauFwmA3LNtG4WSOSGye5HUeqwyefwP33Q3Ex3HsvHGjMc9LPPwfWAyQUQPGz2NBaTOo/nPuDy7DbbofIDkj+A574E6rxBYiISFVZawkEAsTExNTZN9MitdnatWt54403CAQCHHXUUQwdOtTtSJWmgrUes9ayceNGPB4PmZmZbseptEB5kNuOvo/Vi9YD8H8f3kXvoT1cTiVRl5EBzz5b8cefcgr87nf/3969h9lRlvne/97dne6cE5IQCKcQCEk4GA4JkRGEMEhQGJURURkcRQZRGXyNiiOztzqMzrzqDLPfLVsUUWIYFRRHDm5QYNQJTETkIOEkCYQQIJBzAknspNOH5/1jVYeVTh+TdFet7u/nuupirapaVfdaz3pC/1Y9VUW6+SaYWE36p32BBtj+0I5V0msfh+aXS09e/ztS7Uyiet+9W7ckaY9s27aN+fPns2bNGkaPHs3FF1/M8OH+RGxlQwAAIABJREFUOC3tTbfddhvbtm0D4He/+x1vetOb2HffyvqbqCrvArSrpqYmtm7dusfbue2225g/fz7z5s3j7rvv3guV9a1H7l7ES8+8QkN9Aw31DVz/uR/kXZKKoKoKbrwRtr5KWngUHFwHDIbBZb8YtmwofwG0vNbXVUqSuvDII4+wbt06Ukq8/vrrLFy4MO+SpH6nqfU6IZSGBDc3N+dYze4xsBbMc889x9e//nWuvvpqbr/99t2+h2l9fT1PP/00jY2NNDY28tBDD1XcF3TIiMFvjLmvCoaPHpZzRSqSqN6XGPszGPZxYuQXiRF//8bCYR8FBkMMhUHTocaL66j/Symxbds2732tiuV3V+rc9u3bqW+9p303nXPOOdTU1FBdXc20adPYb7/9eqm63uOQ4IK58847d/wS8sc//pGTTjppt4bzDho0iKqqKlpaWgCora2lqqqyfp847vRjOOsjp/PL7/2GfQ8ey9zvXJp3SSqYqJlEjJi7y/yq4ZeR6v4c0mYYdDwRlfXdl3pq48aNzJs3j/r6evbff38uuugibyumwps5cyZPPvkk69atY+TIkbz1rW/t+kXSAPXUU09x++23A3DCCSdwduvtAbswbdo0rrjiCrZv386ICr0VoIG1YNqGyt0NmYMGDeL9738/d955J1VVVZx77rkVdzGDiODya/6Gy6/5m7xLUQWKQdPyLkHqM7/+9a/ZsmULAGvXruXxxx9n5syZOVcldW7w4MF8/OMfp6Ghgbq6uor7O0Xa2xoaGrjlllvYf//9Ofnkk6murt6x7Oc///mO0ZKPPfYYJ510EmPGjOnWduvq6qirq+uVmvuCgbUXPfzwwyxYsIChQ4fyvve9r1snOJ977rncdNNNNDY2MmPGDMaPH7/b+588eTJz5+569Enqbx5Yu4QXtqzl1PFHcvCwsXmX02+s2rCZlRs2Me2Q8Qyp9WhdkUUEEVHxty7QwBMRDB48OO8ypNytXLmSjRs3smTJEp577jkaGho488wzAVi1ahWNjY071m1ubt4pzPZ3BtY9tH79ehYsWEBtbS2nn376jqvbbdiwgXvvvZempibq6+u55ZZb+Nu//dsutzdx4kSuvPJKWlpaBtQXUeqp5tRCFcFPX3qQbz57N80tLXx36a+46eRPccDQffIur+L99unlfO47/5fq6ipGDRvMzf/zg4wYUrm/zvZ3Z5xxBi+99BKbN29mv/32Y/r06XmXJEnqgVWrVu143NTUxIsvvrjj+bp166ipqdlx2uDQoUMZNWpUn9eYFwPrHmhqauKGG25g69atVFVVsWLFCj7xiU8AsHXr1p1+4e7JCdIRsdfDakrJX9zVb/zohf/mm8/ew6CoZt/BI9nWXPrVsYZq/rBhGQcMnZFzhZXv+jt/x7bGJmiElpRY+OQLvGOWw6yLavTo0cydO5fGxkZqa2vzLkeS1EMTJ05kyZIlRAQ1NTUcc8wxO5ZNmjRpp+sSnHbaaXmUmBsD6x7YsmXLjl86WlpaWLt27Y5gOGHCBA444ABWrlxJS0tLbjfp3bBhA//+7//Opk2bOPLIIznvvPMq7uJLUrmN27fwrefupTm10JxaWFm/kbqqGhpamkgpMXlE5d1zuIjGjR5OdVXQ3JIgwZiRQ3davvKF1WzesIXDjzvU0SAFERGG1QrQ0NDAa6+9xtixY6mp8c8wSSVjxoxh3LhxHHTQQYwbN46pU6fuWDZs2DAuu+wynnvuOcaMGcPEiRNzrLTv+S/lHhg5ciQjR47k9ddfB0q/jLQexayqquJDH/oQq1atYsiQIeyzTz5DFH/5y1+yadMmUkosXbqUJUuWcOSRR+ZSi7Q3NGdXvm5VFVX81aGnsGTTSs47ZBbTRh2YU2X9y/+44M95fctWlq/ewHtOmc6sqQfvWHbn9f/Jtz89n6qqYNqsI/javV8wtErdsGbNGubNm0dKicGDB3PppZcybFjf37KtqamJBQsWsH79embNmsWkSZP6vAZJu6qpqeGUU05pd9nw4cM5/vjj+7iiYjCw7oGqqiouueQSFi1axKBBgzj22GN3WX7AAQfkVF1JU1PTjotwpJQq7l6sUlvjBo/k/RPfwk+WPwABnz/qXbzzIK+GureNHTmM7332fe0um//FH7N963YAFj/0HM8vWs6UGd7rVurKwoULaWhoAEoXTXn88cd5y1ve0ud13HXXXTz11FM0NTXx/PPPc+mllzJu3Lg+r0OSusPAuocGDx7MSSedlHcZHZozZw433ngjzc3NjB8/nmnTPAdNle//mfoOPjzpNGqqqhlW44WA+trIMcN5fe0mAFqaEyP2GZ5zRVJlGDJkCNXV1TQ3N1NVVZXbbSZWrFix45SmiGDNmjWdBtZNmzZxyy238NprrzFr1ixOPfXUvipVkgys/d2ECRO44oorqK+vZ8SIEV54SRWlsbmZQR0MNR1VO7Td+ep9X/zpZ/nK+f/Gpg1buOjL72fCYfvlXZJUEWbPns3KlStZtWoVkydP5rjjjsuljmOOOYbf/va3NDc3ExEcfPDBna7/85//nFdffZWUEgsXLmTSpEldvkaS9hYD6wBQU1PDyJEj8y5D6rbGpmY+9a3b+f3il9h/nxF89zPnc8DYgXP59qKbdMwhzHvmG3mXIVWcIUOGcPHFF+ddBqeeeir77rsvGzdu5KijjmLEiBGdrr958+ad7vG7ZcuWvihTkgDwcrE5aG5uZtWqVT261Y00kNzzyBIef/5VUoLVG7fwv2/977xLkqROtbS08Ic//IGFCxeyefPmvMvpVERw1FFHcfLJJ3fropCnn346NTU11NbWMmLECA4/3HPWJfWdHh1hjYgjgCeB/0gpfTCbdwZwLXAI8HvgopTSi9myfYFbgCnAF1JK38/mHwR8AzgNGAS8BPxbSmn+XnhPhdbQ0MB3v/vdHb9WXnjhhQPu0tQamFJK/GTBIu5/chmzph3Ch8+c2eEQ9aaWFlLr60g0Nbe0u54kFcUdd9zBM888Q3NzMw8++CCf/OQncztHdW+bNm0al19+OZs2bWLChAnejkdSn+rpEdZrgYdbn0TEOOBW4IvAGOAR4Cdl688FvgdMBS6NiNaTzn4AvAxMBMYCHwJW70b9FWfJkiVs2rSJ7du309jYyG9+85u8S5L6xM3/9Rj/6z/u58FnXuKa2xby8W/8rMN1z5o5lcP2H8ugmmpGDRvMJ89t/xLvklQUzz77LI2NjbS0tNDU1MTatWvzLmmvGjVqFAcffLBhVVKf6/a/OhHxAeA14AFgcjb7PcDTKaWfZutcBayLiGkppcVAddlUA7QeTjkR+HRK6U/Z88f28H1UjMGDB+94HBEMHeqFY7R7li56gSUPLWXQIXlX0j2PPreCprJ7qD763Ao2bKpnzMhd+8CQ2kH84MoL2Lh5KyOG1XV44SVJKooDDzyQ5cuX09zcTEqJMWPG5F2SJPUL3TrCGhEjgS8Dn22z6Gjg8dYnWQB9PpsPpWG/lwKLge+VBdQHgWsj4gMRUSF/bndfSolHHnmEG2+8kfvvv3/HhQoAjjjiCI477jhqa2vZf//9Ofvss3OsVJXq8QVPM/eUL/Dtz8znpWdWsOSR5/MuqUtzTpiy0/PqCIbUDepw/YhgzMihhlVJFeH8889n1qxZvOlNb+IjH/mIP0hL0l7S3SOsXwFuSCm93Oacs+FA2zEvrwMjAFJKK4H2xvKdD3ye0lDiaRHxJPDRlNLD7axbcRYvXsy9995LY2Mjr7zyCrW1tTvu1RoRnH322QZV7ZFf/fA+Guq3A5BaEgt/9iBTZxb7IhhnnTiNNa//iRt+8Xuqq6v4hw/N6TSwSp1pamri17/+NatXr+bEE0/kyCOPzLskDXB1dXXMmTMn7zIKZc2aNTz44IMMHz6cU045hdra2rxLklSBugysEXEc8Dbg+HYWbwHa3i9lJNDp5fFSShuBK4Ers/NgrwZuj4iDUvnhyAq1evVqGhsbAXaEVmlvOvy4SdQNfYCG+gaiKpj0psoYqPDXb5vBX79tRt5lqB+4++67efzxx2lqamLFihWMHj2aCRMm5F2WpMzWrVuZN28eDQ0NVFdXs3LlSi688MK8y5JUgaKrfBgRc4F/5o0QOpzSOanPANcBH04pnZytO4zSEdcTsnNYu1dExDGUrj48LqW0vu3y6dOnp2uuuaa7m8tdY2Mj69e/8Tb22WefPrlS4JYtWxg+fHiv70fFsGHlRv60aSsjxg9j9D7eo3SgGej9fd26dTt+GIwIRo0axZAhQ3KuqvcN9HYfiCq1zVv/Fiq/f+v++++fc1WVo1LbXXtmoLf76aef/mhKaWbb+d0ZEnw98OOy51cAhwKfyJ7/a0ScB9wFfAl4oquwGhFfp3Sl4MXAkGxbS9sLqwC1tbXMnj27G6UWx5o1a1i+fDkHHHAABx10UJ/sc8GCBRX3OWnP2e4D00Bv94ceeohf/epXpJSorq7msssuY+TItgN++p+B3u4DUaW2eUNDA9dccw1bt26lurqaqVOnVuT7yEultrv2jO3evi4Da0qpHqhvfR4RW4BtKaW12fPzgG8CP6R0H9YPdGO/Q4HbgAnA1ux17+pp8UU2fvx4xo8fn3cZktQvzZo1i3HjxrFu3TqmTJkyIMKqVEnq6uq49NJLWbRoEUOHDuWEE07IuyRJFarHN9NKKV3V5vmvgGk93MYne7pfSZLKHXbYYRx22GF5lyH1qZaWFtavX8/w4cMLPwx+1KhRnHbaaXmXIanCefdnSZKkCtDU1MT8+fNZs2YNKSXOP/98pkyZ0vULJamCdes+rOrY8uXL+c1vfsPSpUvzLkWSJPVjy5YtY+3atTQ2NtLU1MTNN9/MQw89lHdZktSrPMK6B1544QVuvvlmGhsbGTRoEO95z3uYNq1Ho6MlSZK6pa6ujrZ3d7jnnns44YQTqKnxTzpJ/ZNHWPfA888/v9P9VpcsWZJzRZIkqb865JBDvHiRpAHHwLoHDjnkEAYNGgTAoEGDOPTQQ/MtSJIk9VsRwdvf/nbOPvtsqqqqqKqq4h3veIdHVyX1a/4LtwemTJnCu9/9bp599lkOO+wwpk+fnndJkiSpnzvxxBM5/vjjiQiqq6vzLkeSepWBdQ8dffTRHH300XmXIUkVr6mpiZ/97GcsW7aMgw46iPe///3U1tbmXZZUSB5VlTRQOCRYklQIjz76KEuXLmX79u28+OKLLFy4MO+SJElSzgyskqRCqK+vp7m5GYDm5mbq6+tzrkiSJOXNwCpJKoQZM2YwZMgQ6urqGDx4MCeddFLeJeUipbTLrUskSRqoPAFCklQII0eO5FOf+hTr169nzJgx1NXV5V1Sn3vllVe46aab2Lp1K29+85s566yzuvW6pqYmNm/ezKhRo6iq8rdoSVL/YWCVJBVGbW0tEyZMyLuM3Nx+++07hkI/+uijTJ8+vcvPY8OGDdxwww00NjYybNgwPvrRjzJ06NC+KFeSpF7nz7CSJBVEU1NTp8/bc//991NfX09jYyObN29m0aJFvVWeJEl9zsAqSVJBnHPOOdTU1FBdXc3hhx/OQQcd1OVramtrdwwDjghvdyJJ6lf8v5qkQkgpsXz5crZv387kyZOprq7OuySpz02ePJnPfe5zNDQ0MHz4cCKiy9fMnj2bl19+mdWrV3PIIYdwwgkn9EGlkiT1DQOrpEL45S9/yaJFi4gIxo8fz8UXX9ytP9al/qa2tpba2tpurz906FA+9rGP9WJFkiTlxyHBkgrhD3/4A42NjWzfvp3Vq1ezcePGvEuSJElSzgyskgph5MiROz33KqeSJEkysEoqhAsvvJCJEyey//77c8EFFzB48OC8S5IkSVLOPIdVUiGMHTuWiy66KO8yJEmSVCAeYZUkSZIkFZKBVZJUcbZt28bSpUtZv3593qVIkqRe5JBgSVJFqa+v57rrrqOhoYGWlhbe+973MnXq1LzLkiRJvcAjrJKkivLcc8+xbds2tm/fTlNTEwsXLsy7JEmS1EsMrJKkilJ+C6SqqipGjx6dYzWSJKk3GVglSYXU2NjIihUr2Lx5807zJ02axKmnnsqoUaM4/PDDOfvss3OqUJIk9TbPYZUkFc62bdv4zne+Q319PSklLrjgAiZNmrRj+SmnnMIpp5ySY4WSJKkveIRVklQ4ixcv5k9/+hPbt2+nsbGR++67L++SJElSDgyskqTCGTZs2I7HVVVVDB8+PMdqJElSXgyskqTCmTx5MjNnzmTIkCEceOCBvOMd78i7JEmSlAPPYZUk7aK5uZlf/OIXLF++nClTpjBnzhwios/2HxHMmTOHOXPm9Nk+JUlS8XiEVZK0iwceeIAnnniCDRs28Oijj7Jo0aK8S5IkSQOQgVWStIuNGzfS1NQElG4vs3HjxpwrkiRJA5GBVZIKZtu2bSxevJhVq1blVsOMGTMYNGgQdXV11NbWMn369NxqkSRJA5fnsEpSgWzbto1vfetbNDQ0kFLinHPO4dhjj+3zOg488EAuu+wyVq1axQEHHMDIkSP7vAZJkiSPsEpSgSxbtoyGhoYd9x/97W9/m1sto0ePZtq0aYZVSZKUGwOrJBXI6NGjSSkBpfuPjh07NueKJEmS8mNglaQCOeCAAzjrrLMYN24ckydP5p3vfGfeJUmSJOXGc1glqWBmzJjBjBkz8i5DkiQpdx5hlSRJkiQVkoFVkiRJklRIPQqsEXFERGyLiB9mz0+KiP+MiA0RsTYifhoRE8rW3zci/isiXomIj5TNPygifhYR6yLi9Yh4MiIu2mvvSpIkSZJU8Xp6hPVa4OGy5/sA1wOHAhOBzcD3y5bPBb4HTAUujYih2fwfAC9nrxkLfAhY3cNaJEmSJEn9WLcvuhQRHwBeAx4AJgOklH7ZZp1vAveVzaoum2qAyOafCHw6pfSn7Plju1O8JEmSJKn/6tYR1ogYCXwZ+GwXq54KPF32/BvApcBi4HtlAfVB4NqI+EBEHNKzkiVJkiRJA0F3j7B+BbghpfRyRLS7QkRMB74EvLt1XkppJXBKO6ufD3we+CIwLSKeBD6aUnq4nXUlSZIkSQNQpJQ6XyHiOOBHwPEppe0RcRUwOaX0wbJ1JlMaCnxlSukHPSogYhxwNXAmcFBqp6Dp06ena665piebHZC2bNnC8OHD8y5Dfcx2H5hs94HJdh94bPOByXYfmAZ6u59++umPppRmtp3fnSOssyldVOml7OjqcKA6Io5KKZ0QEROBXwFf6WlYBUgprYuIq4EPA2OA9W3Xqa2tZfbs2T3d9ICzYMECP6cByHYfmGz3gcl2H3hs84HJdh+YbPf2decc1uuBw4Hjsuk64C7grIg4EPgNcG1K6bru7jQivh4Rx0RETUSMAD4BLE0p7RJWJUmSJEkDU5eBNaVUn1Ja1ToBW4BtKaW1wCXAYcA/RMSW1qkb+x0K3EbpqsPLKN3e5l27/S4kSZIkSf1Ot29r0yqldFXZ438E/nE3tvHJnr5GkiRJkjSwdOu2NpIkSZIk9TUDqyRJkiSpkAyskiRJkqRCMrBKkiRJkgrJwCpJkiRJKiQDqyRJkiSpkAyskiRJkqRCMrBKkiRJkgqpJu8CJEl9Y926ddx0001s3ryZWbNmceaZZ+ZdkiRJUqc8wipJA8Qdd9zBxo0baWpq4uGHH+bll1/OuyRJkqROGVglaYDYtm3bjscRQUNDQ47VSJIkdc3AKkkDxJw5cxg0aBCDBg1iv/32Y9KkSXmXJEmS1CnPYZWkAeKII45g7ty51NfXM3bsWCIi75IkSZI6ZWCVpAFk6NChDB06NO8yJEmSusUhwZIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAipZR3DV2KiEfyrkGSJEmS1GvWpZTe3nZmRQRWSZIkSdLA45BgSZIkSVIhGVglSZIkSYVkYJUkSZIkFZKBVZIkSZJUSAbWgoqIyyPikYhoiIj5ZfMPjYgUEVvKpi+WLd83Iv4rIl6JiI+UzZ8fEdvbvO7xPn5b6kRE1EXEDRHxYkRsjojHIuIdZcvPiIjFEVGftfHEsmW2e4XqrN3t7/1bRPwwIlZGxKaIeDYiLilbZn/vpzpqd/t7/xcRR0TEtoj4Ydk8+3o/17bd7es9Z2AtrleBfwLmdbB8dEppeDZ9pWz+XOB7wFTg0ogYWrbsX8peMzyldGzvlK7dVAO8DJwGjAK+CNyS/cM2Drg1mzcGeAT4SdlrbffK1WG7l61jf++fvgocmlIaCbwL+KeImGF/7/fabfey5fb3/uta4OHWJ/b1AWOndi9jX++mmrwLUPtSSrcCRMRM4KAevLS6bKoBYu9Xp96QUvoTcFXZrDsj4gVgBjAWeDql9FOAiLgKWBcR01JKi7HdK1YX7f5oFy+33StYSunp8qfZdDiltre/91OdtPv6Ll5qu1ewiPgA8BrwADA5m/0e7Ov9Wgft3hXbvQ2PsFauFyNiRUR8P/uFrtU3gEuBxcD3sj+GVYEiYj9gCvA0cDSwY9hH1q7PZ/PBdu832rR7K/t7PxUR34qIekptuBL4Bfb3fq+Ddm9lf+9nImIk8GXgs20W2df7sU7avZV9vZsMrJVnHXAiMJHSr/AjgB+1LkwprUwpnZJSmpBS+k6b114REa+VTTf2XdnqiYgYRKldb8x+ZR0OvN5mtdcptb/t3k+00+72934upXQZpXZ9K6WhgQ3Y3/u9Dtrd/t5/fQW4IaX0cpv59vX+raN2t6/3kEOCK0xKaQulcxwAVkfE5cDKiBiZUtrUxcuvTil9oXcr1J6KiCrgB8B24PJs9hZgZJtVRwKbu7FJ270CtNfu9veBIaXUDCyMiA8Cn8D+PiC0bfeU0jXY3/udiDgOeBtwfDuL7ev9VGft7v/be87AWvlS9t8BP769P4iIAG4A9gPOTik1ZoueBj5ctt4wSuc8Pb3LRlRxOmn3tuzv/VsNb/Rr+/vA0drubdnf+4fZwKHAS6V/6hkOVEfEUcB12Nf7q9l00O4ppRParGtf74JDggsqImoiYjDZSdcRMTib9+aImBoRVRExFrgGWJBSajukRJXp28CRwDtTSlvL5t8GHBMR52Xfiy8BT2TDRlX52m13+3v/FRHjI+IDETE8Iqoj4izgAuA32N/7rc7a3f7eb11PKYQel03XAXcBZ2Ff7886bHf7es8ZWIvrC8BW4Ergg9njLwCHAXdTGi7yFKXzXi7o5jb/Lna+d9O6vV+2dleU7r32MUr/sK0qa6cLU0prgfOAfwY2Am8GPtDNTdvuBdZZu2N/788SpeG/Kyj16auBuSmlO+zv/VqH7Y79vV9KKdWnlFa1TpSGAW9LKa21r/dfnbU79vUei5RS12tJkiRJktTHPMIqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSqkmrwL6A8i4nBgeN51SJIkSVI7Xk0prc27iN0RKaW8a6hYl1xySbrvvvvYuHQztQzJuxyp8Bqmjs+7BKkiVB2yLe8SpIrwF2OX512CVBGe+vEmTga+AQeklFbmXU9PGFh3Q2tQramp4bTTTuP569fnXZJUEdZ/7C15lyBVhJp3V+SP4FKfe/C4/8i7BKkyTHiOW4F/Ad5CZQVXhwQf5B6bAAASYElEQVT3QHlQPeOMM5gwYQIAz2NglSRJklRMVcB7gfcAtwIPwKtzIyoiuBpYu6GjoCpJkiRJlaISg6uBtRMGVUmSJEn9TSUFVwNrOwyqkiRJkvq7SgiuBtYyBlVJkiRJA02Rg6uBFYOqJEmSJBUxuA7owGpQlSRJkqSdFSm4DsjAalCVJEmSpM4VIbgOqMBqUJUkSZKknskzuA6IwGpQlSRJkqQ9k0dw7deB1aAqSZIkSXtXXwbXfhlYDaqSJEmS1Lv6Irj2q8BqUJUkSZKkvtWbwbVfBFaDqiRJkiTlqzeCa0UHVoOqJEmSJBXL3gyuFRlYDaqSJEmSVGx7I7hWVGA1qEqSJElSZdmT4FoRgdWgKkmSJEmVbXeCa+EDa0QcHxGce+65HHvssXmXI0mSJEnaA63B9S+Bi4AD4VUgOlq30FJKj1100UU8/vjj3HLLLaxevTrvkiRJkiRJu6kFuA14CzAKeAUO7mjdwgdWgHnz5sXzzz8fb37zm7nnnnsMrpIkSZJUYVqD6p8BvwYegoO/mVKklFZ09JqKCKytDK6SJEmSVFl2J6i2qqjA2srgKkmSJEnFtidBtVVFBtZWBldJkiRJKpa9EVRbVXRgbWVwlSRJkqR87c2g2qpfBNZWBldJkiRJ6lu9EVRb9avA2srgKkmSJEm9qzeDaqt+GVhbGVwlSZIkae/qi6Daql8H1lYGV0mSJEnaM30ZVFsNiMDayuAqSZIkST2TR1BtNaACayuDqyRJkiR1Ls+g2mpABtZWBldJkiRJ2lkRgmqrAR1YWxlcJUmSJA10RQqqrQysZQyukiRJkgaaIgbVVgbWdhhcJUmSJPV3RQ6qrQysnTC4SpIkSepvKiGotjKwdoPBVZIkSVKlq6Sg2srA2gMGV0mSJEmVphKDaqtIKeVdQ8W6+OKL03333cemZVupY3De5UiFt+3I/fIuQaoIVRO35V2CVBHePXZZ3iVIFeGJH23iROBaOLgSQmo5A+teEBFTgb+m9MOFiuMvsU2KxjYpJtuleGyTYrJdisc2KSbbpXhmp5T+Le8idoeBdS+JiEdSSjPzrkNvsE2KxzYpJtuleGyTYrJdisc2KSbbpXgquU08h1WSJEmSVEgGVkmSJElSIRlY957r8y5Au7BNisc2KSbbpXhsk2KyXYrHNikm26V4KrZNPIdVkiRJklRIHmGVJEmSJBWSgTUTEW+PiCURsTQirmxneUTENdnyJyLihK5eGxFXRcQrEbEom84uW/b32fpLIuKs3n+Hlacv2yQiDo2IrWXzr+ubd1l5eqNdsmWfzJY9HRH/UjbfvtKFvmwT+0r39NK/Xz8p+9yXR8SismX2k27oy3axr3RPL7XJcRHxYPa5PxIRs8qW2Ve6oS/bxb7SPb3UJsdGxO8i4smI+L8RMbJsWXH6SkppwE9ANfA8cBhQCzwOHNVmnbOBXwIBnAT8vqvXAlcBV7Szv6Oy9eqASdnrq/P+HIo05dAmhwJP5f2+iz71YrucDvwKqMuej8/+a18pXpvYV3Jqkzav/zfgS9lj+0kx28W+klObAPcC7yh7/YLssX2lmO1iX8mvTR4GTsseXwx8JXtcqL7iEdaSWcDSlNKylNJ24MfAu9us827g31PJg8DoiJjQzde29W7gxymlhpTSC8DSbDt6Q1+3ibqnt9rlE8DXUkoNACmlNWXbsq90rq/bRF3r1X+/IiKA9wE3l23LftK1vm4Xda232iQBrUeKRgGvlm3LvtK1vm4Xda232mQqcH/2+D+B88q2VZi+YmAtORB4uez5imxed9bp6rWXZ4fl50XEPj3Y30DX120CMCkiHouI+yLirXv8Dvqn3mqXKcBbI+L32ed/Yg/2N9D1dZuAfaUrvfnvF8BbgdUpped6sD/1fbuAfaUrvdUmc4F/jYiXgauBv+/B/tT37QL2la70Vps8Bbwre3w+cHAP9tdnDKwl0c68tpdP7midzl77beBw4DhgJaWhQt3d30DX122yEjgkpXQ88BngpvJx/Nqht9qlBtiH0hCWzwG3ZEcr7Ctd6+s2sa90rbfapNUF7HwUz37SPX3dLvaVrvVWm3wC+HRK6WDg08ANPdif+r5d7Ctd6602uRj424h4FBgBbO/B/vqMgbVkBW/8ogBwELsOU+honQ5fm1JanVJqTim1AN/ljUPp3dnfQNenbZINeVifPX6U0lj9KXvt3fQfvdIu2bJbs2EsDwEtwLhu7m+g69M2sa90S2+1CRFRA7wH+EkP96c+bhf7Srf0Vpt8GLg1e/xT/Purp/q0Xewr3dJbfxcvTinNSSnNoPSD2/M92F/fSQU4kTjvidKRhGWUTipuPRn56DbrnMPOJzI/1NVrgQllr/80pbHgAEez84nMy/Ck/7zbZN/WNqB0UvorwJi8P4eiTb3YLh8Hvpw9nkJpGErYVwrZJvaVnNokW/524L4227KfFLNd7Cs5tQnwDDA7e3wG8Gj22L5SzHaxr+TXJq0XVKwC/h24OHteqL6SewMUZaJ0Za1nKf2y8D+zeR8HPp49DuDabPmTwMzOXpvN/0G27hPAz9k5LP3PbP0lZFdMc8qvTSidZP501jn/ALwz7/df1KmX2qUW+CGlcyn+APx52TL7SoHaxL6SX5tky+a3bqPNfPtJwdrFvpJfmwCnAI9mn/3vgRlly+wrBWsX+0qubfKpbP6zwNeAKFtWmL4SWUGSJEmSJBWK57BKkiRJkgrJwCpJkiRJKiQDqyRJkiSpkAyskiRJkqRCMrBKkiRJkgrJwCpJ2qsi4qqIeCrvOvJWaZ9DRMyOiBQR4/KupTsi4qKI2JLj/udHxJ157V+SBgoDqyRViIjYLyL+v4h4LiK2RcSaiHggIj4ZEcPzrq/M1cBpvb2T9gJhRMyKiHVZmKjp7Rr2VBYQW6f6iFgWETdFxCk5lPMAMAFYn9W21wJhm/dZPn28B69/b5vZPwEO2xv1dbHvjoL8p4AP9vb+JWmgK/z/zCVJEBGHAr8FNgFfBJ6g9KPjFOBDlELGTTmVt5OU0hagz498RcTbgNuA64Er0m7caDwialNK2/d6cZ37KHAnUEcpgH0YuD8iPp9S+te+KiJ736t6cRet77Pc67u7sZTSVmDrHlW0B1JKu127JKn7PMIqSZXh20ALMDOl9OOU0h9TSk+llG5NKZ0L3Ny6YkR8JiKeiIg/RcQrEfG9iBhdtnyXI2dtjyJFxKiI+EF2FHdbduRvbtn6H4uIZ7NlayPintYjmm2PfEbEiRFxb3bkc1NELIyIP2uz/xQRl0bET7O6l0VEt49eRcT5lMLQ/5tS+mxrWI2IoyLirojYnL2XmyNi/7LXzY+IOyPi8xGxAlgREYdm9ZwXEf+ZHfn8Y0Sc2WafnW67B15LKa1KKb2YUvqvlNJFwNeAr0bE5O7ur+y9fCpr940R8f2IGFq2zqkR8WBEbImI1yPi9xFxTLZsx3cgImYD3weGlR0NvSoivtT2qHb22t9GxDXdfJ/l09bs9R1+3yJiefb6n2Z1LM/m7/Q9bv3eRcSHI2J59h6/HxG1EXFZRLwcEesj4n9FRFXZ6z4YEQ+Xfa4/jYgDs2WHAv+Vrbo22//88s+7bDt1EfG/I2J19h4ejLIj5WWf7xnZ514fEY9ExAldfG6SNKAZWCWp4CJiDHAWcG1K6U/trdPmaGILMBc4GvgrYBbwf3q4238C3gT8BTANuBh4JatnJnAt8I/AVOBtwN2dbGsE8APgrVkti4BfxK5DLL8E3AEcS2m457yImNhVoRHxMeBHwOUppa+WzZ8A3A88le33bcBw4OflgYXS8OXpwNuBM8rm/zNwTVbPw8CPIxt63YNt765/o/T/6HN7uL+3Asdky98P/CWloatE6QeFO4CF2Xt6M/ANoLmd/T9A6TtUT2mY8ARKQ73nAdMiYlbrihExFXgLcMMevN8Ov2/Aidl/P5rVceIur37DocC7s+2cB5xP6T2fCMwBLgE+SelzaVUL/AOlz+QvgHG88QPQy9l2oNSfJpB9nu34F0qf+cXA8cCTwN1Z25X7KnAlcAKlkRE/iojo5D1J0sCWUnJycnJyKvBEKVgk4C/bzF9BaejtFuC6Tl7/dqABqMqeXwRsabPO7Gwf47LnPwe+38H23kNpKOeIDpZfBTzVST0BrAQ+WDYvAV8te15DKSx9sJPtXJW9rwT8bTvLvwz8us28fbL1Z2XP5wNrgbqydQ7N1vlY2bwDs3mn9GDbnX4OZe/7vR0sWwV8q4fv5WWgpmyd7wK/yh6PydY/rYP9tf0O7PI9yebfWf59A74OPNKN97m17PvaOr2pq+9bR59T2/qyz3srMKps3n9k7VtbNm8B8M1O9jUt299B7X0uZevNB+7MHg8DtgMfKlteDTwP/FOb7ZxVts7J5ftycnJyctp18girJFWutwLHAQ8Bg1tnRsSfZ0NZV0TEZuBWSkeRejJc9dvA+yLi8Yi4OiLKL6L0n8CLwAsR8aNsCOaIjjYUEeMj4jtRGkL8OrAZGA8c0mbVJ1ofpJSaKAWN8V3UuZLS0c9PR0Tb7c0ATs2Ghm7Jho++nC07vGy9p1JKDe1s+4myx69m/22tp7vb3hNBKcz0ZH9/zD678rrHA6SUNlAKWfdkQ4s/ExEH70Zd3wU+EBFDIqIa+Gu6d3T1c5S+r+XTkmxZZ9+3nngp7Xxu6Wrg2bTzecmrKfteRcQJEXFHRLyY9ZdHskVtv0+dORwYROk8cwBSSs3A74Cj2qzb2fdKktSGgVWSim8ppeAyrXxmSumFlNJSSkciAciG0N4FPENpOOQMSkMUoRRaoTRkuO0QxEFttv1LYCKlYaDjgLsi4vvZss2UhjO+D3gJ+HtgcUQc0EH9N1IakvlpSkNHj6N0dLi2zXqNbZ4nuv7/1BZKw1/XAAvaDCGuovRZtA1JR7DzxX/aHWZdXk9KqTU4VpX9tzvb3i3ZcOl9gWU93F+nn2FK6SOUjtjfD7wLeDYizupheXdR+s6dB5wNjKbsHOpOrEopLW0zbc/q6vD71kPtvf8OP5OIGAbck72fv6b0PX17tl7b72dnWvtTexf6ajuvsZ1l/j0mSR3wH0hJKriU0nrgXuDy6Pr2NTMp/aH96ZTS71JKzwJtg+RaYGhEjCybd1w7+12XUvpBKl0E6G+AD0dEXbasKaX0m5TS31M6/3MYpfP/2nMK8H9SSnellJ6mdIS17Xl9uy2ltInSOb6vAPdFxKRs0R8onXf4YjtBafMe7rY3tw3wWUo/LNyxt/eXUno8pfT1lNJsSsNjP9zBqtspDWtt+/omSkdqL86mW1NKr/Wkhg7q6vD7Rink7VLLXjCNUkD+Hyml+1NKi9n1aGfr0dnO9r80W6/8IkvVwJ8Bf9x75UrSwGNglaTKcBmlf7MfjYgLonTF2CkRcQGli8W0XjjnuWy9uRExKVs+t822fk/pqOJXI2JyRJyXbX+HiPhyRJwbEUdExJGUzltdllJqiIi/iNKVaI/Pjmj+FaULKz3TQe3PAh/Maj4R+DFvhIC9IgttbweWUwqth1O6MNQo4CcR8eaIOCwi3hYR13c2hLmb9ua2R0fE/hFxSEScnl2F9vPAldkR9L2yv+z78LWIeEtETIyI0yn92NBRoFoODI6IM6N05eChZcu+R+liVX9B9y+21Po+y6fWi1h1+H0rq+WM7DX7dHN/3fESpfOgL88+03OAr7RZ50VKR0LPiYh92/vRKJUuhvZt4GsRcXb2Hr4N7Ad8ay/WK0kDjoFVkipASmkZpSuP3k3pD+rHKB11+wylP4jnZus9Qekqpp+hFEQuAa5os60NwIXAmZSuZHoppXu7lmugdJXcxymdlzcCeGe27DVKV6/9FbA42/4lKaX/7qD8iyld0fZRSmF1HqUAsldloeFsSgH5PkpHfU+mdKTybuBpSsGvIZv2ZF+v7sVtf5fSubjPUvps6oDZKaWr9/L+6indt/en2b5upHR15a+3t3JK6QHgOkrDfdcCf1e2bBmlz/glSkdpu6P1fZZPV2bLOvu+QemI8+mUztt9rJv761JKaS2lI8znUuov/0Cp75Sv80o2/58pnf/6zQ4293ngFkq3A1pEduXplNLKvVWvJA1E8cZpOZIkSd0TEX8EfpRS+ue8a5Ek9V81eRcgSZIqR0SMBy6gdPuf7+RbjSSpvzOwSpKknlgNrKN0n9p1eRcjSerfHBIsSZIkSSokL7okSZIkSSokA6skSZIkqZAMrJIkSZKkQjKwSpIkSZIKycAqSZIkSSokA6skSZIkqZD+fydvgigX1zKsAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(16,17)) #(17,5)\n", "\n", "cbar_ticks = [0.005, 0.006, 0.007, 0.008, 0.009] \n", "\n", "ax1 = plt.subplot(111, projection=ccrs.PlateCarree())\n", " \n", "# ax1.set_xlim([-35.1, -18])\n", "# ax1.set_ylim([29.9, 40])\n", "\n", "ss = ax1.scatter(kde_x, kde_y, c=kde_z, s=10, vmin=cbar_ticks[0], vmax=cbar_ticks[-1], cmap=plt.cm.get_cmap('viridis', len(cbar_ticks)-1), transform=ccrs.PlateCarree())\n", "ss.cmap.set_over('r')\n", "ss.cmap.set_under('gray')\n", "\n", "ax1.coastlines(resolution='50m', color='black', linewidth=1, zorder=50)\n", "ax1.add_feature(cartopy.feature.LAND, facecolor='burlywood', zorder=20) \n", " \n", "gl = ax1.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)\n", "gl.top_labels = False\n", "gl.right_labels = False\n", "gl.xformatter = LONGITUDE_FORMATTER\n", "gl.yformatter = LATITUDE_FORMATTER\n", "gl.ylocator = mticker.FixedLocator([-42., -40., -38., -36., -34.])\n", "gl.xlabel_style = {'size': 12}\n", "gl.ylabel_style = {'size': 12}\n", "\n", "cbar = plt.colorbar(ss, extend='both', orientation='horizontal')\n", "cbar.ax.get_yaxis().labelpad = 15\n", "cbar.ax.set_xlabel('Gaussian Kernel Density Estimation', size=14) #, rotation = 270\n", "\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***--> in this test case the GKDE does not show so clearly the accumulation zones as more particles should be released in the region of study***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4. Quantification" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A way to quantify then the above plot is to calculate the percentage of particles with a \"high\" GKDE value, that we choose it to be 0.008" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* first check if any particles deleted:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([26.698011, 28.164171, 27.643606, 24.718645, 37.995255, 29.234638,\n", " 38.174526, 38.089085, 36.386173, 36.925137, 37.271023, 27.692135,\n", " 34.057995, 33.911255, 36.052097, 27.156801, 37.874683, 25.212149,\n", " 44.984165, 32.263596, 29.785528, 36.0202 , 33.17997 , 30.311327,\n", " 23.818752, 25.614098, 24.09303 , 35.909134, 35.44305 , 24.450853,\n", " 38.81421 , 31.958689, 34.563725, 33.417557, 33.89991 , 32.049088,\n", " 21.763323, 14.18421 , 22.636242, 34.797607, 24.023075, 24.054976,\n", " 31.82561 , 45.362434, 45.066288, 45.34612 , 35.21189 , 35.262207,\n", " 26.059761, 15.283366, 24.001287, 22.936937, 22.500818, 24.380352,\n", " 21.790487, 21.282694, 34.865032, 45.920753, 45.098103, 34.195713,\n", " 22.949871, 13.661928, 27.963894, 23.287018, 23.10574 , 21.98923 ,\n", " 24.398281, 29.69943 , 30.03786 , 28.911987, 35.554043, 36.454082,\n", " 24.802662, 15.132946, 15.482158, 14.75911 , 22.573833, 23.885687,\n", " 21.990768, 23.868902, 23.409975, 22.211067, 27.914537, 30.84707 ,\n", " 27.61835 , 22.490934, 21.765697, 13.951972, 16.424952, 24.541779,\n", " 24.06132 , 29.723223, 34.884018, 36.20441 , 37.77424 , 40.68846 ,\n", " 18.142883, 22.070595, 24.364391, 27.339596, 27.215214, 26.326141,\n", " 23.880114, 22.573923, 22.02631 , 29.564281, 30.349543, 29.307894,\n", " 20.394741, 20.46407 , 23.062609, 23.215126, 19.929914, 15.229213,\n", " 20.331572, 24.010874, 22.541416, 26.053875, 23.698732, 28.86736 ,\n", " 16.703207, 20.428095, 14.482446, 16.92354 , 16.34547 , 15.432422,\n", " 18.356052, 18.807606, 30.947237, 13.901372, 22.827028, 19.978989,\n", " 19.663671, 17.329836, 17.753399, 17.171776, 16.870941, 17.327799,\n", " 17.529312, 18.485256, 18.942495, 17.336267, 13.677705, 17.166965],\n", " dtype=float32)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds['lon'][:,-1].values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* we calculate now the percentage of particles with a high GKDE value:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "26.39" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nondelp = sum(~np.isnan(ds['lon'][:,-1].data)) # calculating the total non-nan particles (none present in this \n", " # test case , but there could be some particles that values are\n", " # nan due to for example boundary conditions)\n", "\n", "perc_high = np.round(((np.count_nonzero(kde_z >= 0.008)) / nondelp) * 100, 2)\n", "perc_high" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.18" } }, "nbformat": 4, "nbformat_minor": 4 }