diff --git a/notebooks/tutorial_RGDR.ipynb b/notebooks/tutorial_RGDR.ipynb
index 29a7389..97ae938 100644
--- a/notebooks/tutorial_RGDR.ipynb
+++ b/notebooks/tutorial_RGDR.ipynb
@@ -51,7 +51,7 @@
"outputs": [],
"source": [
"target_timeseries = target_resampled.sel(cluster=3).ts.isel(i_interval=0)\n",
- "precursor_field = field_resampled.sst.isel(i_interval=1)\n",
+ "precursor_field = field_resampled.sst.isel(i_interval=slice(1,5)) # Multiple lags: 1 through 4\n",
"\n",
"rgdr = RGDR(eps_km=600, alpha=0.05, min_area_km2=3000**2)"
]
@@ -83,7 +83,7 @@
],
"source": [
"fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 2))\n",
- "_ = rgdr.plot_correlation(precursor_field, target_timeseries, ax1, ax2)"
+ "_ = rgdr.plot_correlation(precursor_field, target_timeseries, lag=1, ax1=ax1, ax2=ax2)"
]
},
{
@@ -116,9 +116,10 @@
"source": [
"fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 2))\n",
"\n",
- "_ = rgdr.plot_clusters(precursor_field, target_timeseries, ax=ax1)\n",
+ "_ = rgdr.plot_clusters(precursor_field, target_timeseries, lag=1, ax=ax1)\n",
"\n",
- "_ = RGDR(min_area_km2=1000**2).plot_clusters(precursor_field, target_timeseries, ax=ax2)"
+ "_ = RGDR(eps_km=600, min_area_km2=None).plot_clusters(precursor_field, target_timeseries,\n",
+ " lag=1, ax=ax2)"
]
},
{
@@ -489,17 +490,13 @@
" stroke: currentColor;\n",
" fill: currentColor;\n",
"}\n",
- "
<xarray.DataArray 'sst' (cluster_labels: 3, anchor_year: 39)>\n",
- "290.8 291.0 290.7 290.1 291.1 291.3 ... 299.0 299.5 298.9 298.9 299.2 298.2\n",
+ "<xarray.DataArray 'sst' (cluster_labels: 6, anchor_year: 39)>\n",
+ "290.8 291.0 290.7 290.1 291.1 291.3 ... 285.6 286.3 286.2 285.5 285.0 285.1\n",
"Coordinates:\n",
" * anchor_year (anchor_year) int32 1980 1981 1982 1983 ... 2016 2017 2018\n",
- " i_interval int64 1\n",
- " index (anchor_year) int64 1 13 25 37 49 61 ... 409 421 433 445 457\n",
- " interval (anchor_year) object (1980-07-02, 1980-08-01] ... (2018-0...\n",
- " target bool False\n",
- " * cluster_labels (cluster_labels) int32 -2 0 1\n",
- " latitude (cluster_labels) float64 36.05 nan 29.44\n",
- " longitude (cluster_labels) float64 223.9 nan 185.4
290.8 291.0 290.7 290.1 291.1 291.3 ... 299.5 298.9 298.9 299.2 298.2
array([[290.79588914, 290.970545 , 290.71731703, 290.0762239 ,\n",
+ " * cluster_labels (cluster_labels) <U20 'lag:1_cluster:-2' ... 'lag:4_clust...\n",
+ " latitude (cluster_labels) float64 36.05 29.44 37.33 29.58 38.14 39.78\n",
+ " longitude (cluster_labels) float64 223.9 185.4 221.8 190.2 217.8 219.3
290.8 291.0 290.7 290.1 291.1 291.3 ... 286.3 286.2 285.5 285.0 285.1
array([[290.79588914, 290.970545 , 290.71731703, 290.0762239 ,\n",
" 291.08960917, 291.31511491, 291.11538436, 290.26277142,\n",
" 290.80443321, 290.99960169, 291.53446464, 291.36075119,\n",
" 291.85483292, 291.09343404, 291.31408735, 291.41374784,\n",
@@ -509,16 +506,6 @@
" 291.44962923, 291.91882282, 290.70922506, 290.48853941,\n",
" 290.34711093, 291.99974208, 292.3259726 , 293.32741257,\n",
" 291.91085874, 291.45966391, 291.37066195],\n",
- " [291.41316602, 290.99213881, 290.7072011 , 290.28445664,\n",
- " 291.24062074, 291.15653517, 290.90280998, 290.26394997,\n",
- " 291.24222615, 291.51096688, 291.63083812, 292.04367136,\n",
- " 291.17972744, 291.09824891, 291.49867568, 291.09047242,\n",
- " 291.78362741, 291.61656609, 291.40476286, 291.22261622,\n",
- " 291.74089263, 292.15751269, 291.59638891, 291.53674801,\n",
- " 292.04091404, 292.05970401, 291.54779359, 291.4457052 ,\n",
- " 291.83593796, 291.9013184 , 291.58505835, 292.59423736,\n",
- " 291.89844604, 292.16743467, 292.15498931, 292.00892894,\n",
- " 291.98830915, 291.94098392, 291.62096043],\n",
" [298.94548042, 298.92119198, 297.71694295, 298.98291705,\n",
" 299.35080249, 298.06875363, 298.51701671, 298.495829 ,\n",
" 299.61375876, 298.35675492, 298.46956715, 299.34460285,\n",
@@ -528,64 +515,43 @@
" 298.88763028, 299.2529288 , 299.0168395 , 298.84020348,\n",
" 298.48441327, 299.30649003, 299.69018872, 299.65241405,\n",
" 299.320106 , 299.04325189, 299.48610574, 298.91044985,\n",
- " 298.89195415, 299.19568083, 298.21053747]])
anchor_year
(anchor_year)
int32
1980 1981 1982 ... 2016 2017 2018
array([1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991,\n",
" 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003,\n",
" 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015,\n",
- " 2016, 2017, 2018])
i_interval
()
int64
1
index
(anchor_year)
int64
1 13 25 37 49 ... 421 433 445 457
array([ 1, 13, 25, 37, 49, 61, 73, 85, 97, 109, 121, 133, 145,\n",
- " 157, 169, 181, 193, 205, 217, 229, 241, 253, 265, 277, 289, 301,\n",
- " 313, 325, 337, 349, 361, 373, 385, 397, 409, 421, 433, 445, 457],\n",
- " dtype=int64)
interval
(anchor_year)
object
(1980-07-02, 1980-08-01] ... (20...
array([Interval('1980-07-02', '1980-08-01', closed='right'),\n",
- " Interval('1981-07-02', '1981-08-01', closed='right'),\n",
- " Interval('1982-07-02', '1982-08-01', closed='right'),\n",
- " Interval('1983-07-02', '1983-08-01', closed='right'),\n",
- " Interval('1984-07-02', '1984-08-01', closed='right'),\n",
- " Interval('1985-07-02', '1985-08-01', closed='right'),\n",
- " Interval('1986-07-02', '1986-08-01', closed='right'),\n",
- " Interval('1987-07-02', '1987-08-01', closed='right'),\n",
- " Interval('1988-07-02', '1988-08-01', closed='right'),\n",
- " Interval('1989-07-02', '1989-08-01', closed='right'),\n",
- " Interval('1990-07-02', '1990-08-01', closed='right'),\n",
- " Interval('1991-07-02', '1991-08-01', closed='right'),\n",
- " Interval('1992-07-02', '1992-08-01', closed='right'),\n",
- " Interval('1993-07-02', '1993-08-01', closed='right'),\n",
- " Interval('1994-07-02', '1994-08-01', closed='right'),\n",
- " Interval('1995-07-02', '1995-08-01', closed='right'),\n",
- " Interval('1996-07-02', '1996-08-01', closed='right'),\n",
- " Interval('1997-07-02', '1997-08-01', closed='right'),\n",
- " Interval('1998-07-02', '1998-08-01', closed='right'),\n",
- " Interval('1999-07-02', '1999-08-01', closed='right'),\n",
- " Interval('2000-07-02', '2000-08-01', closed='right'),\n",
- " Interval('2001-07-02', '2001-08-01', closed='right'),\n",
- " Interval('2002-07-02', '2002-08-01', closed='right'),\n",
- " Interval('2003-07-02', '2003-08-01', closed='right'),\n",
- " Interval('2004-07-02', '2004-08-01', closed='right'),\n",
- " Interval('2005-07-02', '2005-08-01', closed='right'),\n",
- " Interval('2006-07-02', '2006-08-01', closed='right'),\n",
- " Interval('2007-07-02', '2007-08-01', closed='right'),\n",
- " Interval('2008-07-02', '2008-08-01', closed='right'),\n",
- " Interval('2009-07-02', '2009-08-01', closed='right'),\n",
- " Interval('2010-07-02', '2010-08-01', closed='right'),\n",
- " Interval('2011-07-02', '2011-08-01', closed='right'),\n",
- " Interval('2012-07-02', '2012-08-01', closed='right'),\n",
- " Interval('2013-07-02', '2013-08-01', closed='right'),\n",
- " Interval('2014-07-02', '2014-08-01', closed='right'),\n",
- " Interval('2015-07-02', '2015-08-01', closed='right'),\n",
- " Interval('2016-07-02', '2016-08-01', closed='right'),\n",
- " Interval('2017-07-02', '2017-08-01', closed='right'),\n",
- " Interval('2018-07-02', '2018-08-01', closed='right')], dtype=object)
target
()
bool
False
cluster_labels
(cluster_labels)
int32
-2 0 1
latitude
(cluster_labels)
float64
36.05 nan 29.44
array([36.0508552, nan, 29.4398051])
longitude
(cluster_labels)
float64
223.9 nan 185.4
array([223.86658208, nan, 185.40970765])
"
+ " 2016, 2017, 2018])
cluster_labels
(cluster_labels)
<U20
'lag:1_cluster:-2' ... 'lag:4_cl...
array(['lag:1_cluster:-2', 'lag:1_cluster:1', 'lag:2_cluster:-1',\n",
+ " 'lag:2_cluster:1', 'lag:3_cluster:-1', 'lag:4_cluster:-2'], dtype='<U20')
latitude
(cluster_labels)
float64
36.05 29.44 37.33 29.58 38.14 39.78
array([36.0508552 , 29.4398051 , 37.33257702, 29.58134561, 38.13773082,\n",
+ " 39.78162825])
longitude
(cluster_labels)
float64
223.9 185.4 221.8 190.2 217.8 219.3
array([223.86658208, 185.40970765, 221.82516648, 190.20336403,\n",
+ " 217.7810629 , 219.30300121])
"
],
"text/plain": [
- "\n",
- "290.8 291.0 290.7 290.1 291.1 291.3 ... 299.0 299.5 298.9 298.9 299.2 298.2\n",
+ "\n",
+ "290.8 291.0 290.7 290.1 291.1 291.3 ... 285.6 286.3 286.2 285.5 285.0 285.1\n",
"Coordinates:\n",
" * anchor_year (anchor_year) int32 1980 1981 1982 1983 ... 2016 2017 2018\n",
- " i_interval int64 1\n",
- " index (anchor_year) int64 1 13 25 37 49 61 ... 409 421 433 445 457\n",
- " interval (anchor_year) object (1980-07-02, 1980-08-01] ... (2018-0...\n",
- " target bool False\n",
- " * cluster_labels (cluster_labels) int32 -2 0 1\n",
- " latitude (cluster_labels) float64 36.05 nan 29.44\n",
- " longitude (cluster_labels) float64 223.9 nan 185.4"
+ " * cluster_labels (cluster_labels) "
+ ""
]
},
"execution_count": 6,
@@ -624,7 +590,7 @@
},
{
"data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEXCAYAAACzhgONAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABy0ElEQVR4nO2dd5wkZZn4v09PznlnZmdmZ3Nk2WV3ySAZEVBAQEQUBEx4eHqnnnh3htOf5+mdCiqKKIIBEck55wUWNrA5h9ndCTs555l+f3+8VTM9PR2qZzrN7vv9fPrT3dVV1U9Xd9dTTxalFAaDwWAw+MMVawEMBoPBEN8YRWEwGAyGgBhFYTAYDIaAGEVhMBgMhoAYRWEwGAyGgBhFYTAYDIaAGEURJkRkm4icHWSdfxeRP0RHoskhIq+LyOdiLYfBN7H8fkTk+yLy11i8t1NERInIXAfrzbTWTZzAe0x426mGURRhQim1RCn1epB1/lsp5ejPPRX+jKEiIueIyGsi0i4iVSFu+1kRWR0h0WL6/iJSJSK9ItLlcZseifeKF46lk+zRgFEURylx+gfsBv4IfDPabxynx8OTjyqlMj1utbEWKJ6ZAt/nUYVRFGHCuio8P8g6I1aCxxXVDSJySESaROQ/rNcuAv4duMa6utxkLc8RkXtEpE5EakTk/4lIgvXaZ0XkbRH5hYg0Az8UkTYROc7j/YusK9dpIpInIk+LSKOItFqPyyN0eABQSr2vlPoLsD+U7URkEXAXcKp1PNqs5ZeIyAci0iEih0Xk+x7b2Mf3ZhE5BLwqIgki8jPrWB8QkVs9r2r9HV9/7x9JQvl+RGSuiLxhWWpNIvKgx2sLReQlEWkRkV0i8okQZFjisW29iPy7j3XOFpFqr2Uj/wUROUlE1lnfUb2I/Nxa7U3rvs06pqda698kIjusz/yCiFR67FeJyD+JyB5gTwifw+/vxIObRKTW+u6/4bGtS0RuE5F9ItIsIv8QkXw/7/NZEdkvIp3W7+s6pzLGO0ZRxJ4zgAXAecB3RWSRUup54L+BB62ry2XWuvcBQ8Bc4ATgQsDTlXUy+iRcDPwAeBS41uP1TwBvKKUa0N/9vUAlMAPoBX7tRGAR+ZSlhPzdZoR+GPyjlNoBfAl41zoeudZL3cD1QC5wCXCLiFzutflZwCLgw8DngY8Ay4EVgPe69+Hj+AZ4/zGIyG8CHJPNIX7sUL6fHwIvAnlAOfArS54M4CXgb8A04JPAb0RkcbA3F5Es4GXgeWA6+pi8EuJnALgDuEMplQ3MAf5hLf+QdZ9rHdN3ReQy9AXSx4Ei4C3gAa/9XY7+nQf9DB44+Z2cA8xDf+ffktGLvq9Y73kW+ji0And6v4F1rH8JfEQplQWcBmwMQcb4RillbmG4AVXA+UHW+T7wV+vxTEAB5R6vvw980ntd63kx0A+keSy7FnjNevxZ4JDX+50P7PN4/jZwvR/ZlgOtHs9fR58kI3GszgeqQtzms8DqIOvcDvzC6/jO9nj9VeCLXnIoINHh8Q34/pP87XQBbdbt8VC+H+DPwN2evyVr+TXAW17Lfgd8z4FM1wIfOPgdnw1U+/g851uP3wT+Cyj0Wsf+fhI9lj0H3Ozx3AX0AJXWcwWc6/CYKmBuCL+ThR6v/xS4x3q8AzjP47VSYND6zYx8BiDD+u6u9PwNHS03Y1HEniMej3uATD/rVQJJQJ19lYr+00/zWOew1zavAekicrKIzESfbB4DEJF0EfmdiBwUkQ70HzpXLFfWVMD6XK9Z7pl29FV/oddqnsdkutdzz8dOjm8kuVwplWvdLg/x+/k3QID3RWff3WQtrwRO9rRsgOuAEgfyVAD7Jv+xuBmYD+wUkbUicmmAdSuBOzxkbUF/rjKPdbx/40GZwO/kIPq3Ysv0mIdMO4Bh9IXFCEqpbrRi/hL6N/SMiCwMVdZ4xSiK+MW7re9h9BVvoccJJVsptcTfNkqpYbSpf611e1op1Wm9/HW0y+tkpd0CtitAggkmItfJ2Awd71tYXU/2x/Gx7G/Ak0CFUioHHUfwlt9zuzq0a8amwuNxsOMbtM2yiNwV4JhsC7a9F46/H6XUEaXU55VS04Evot1Lc63P9IbH58lV2s1zi4P3PwzMdrBeN5BuP7EUWZGHbHuUUteiFe5PgIctN42v43kYbfF5ypumlHrH8+M6kMkbJ78Tz9/CDMBOJjiMdid5ypSqlKrxfhOl1AtKqQvQVsdO4PcTkDUuMYoifqkHZoqIC0ApVYf2Q/9MRLKtINscETkryH7+hr7Suc56bJOF9nu3WcG57zkVTCl1vxqboeN9O+RrO0vmVPSVu4hIqogke7z+up9AI+jjUe65vvUZWpRSfSJyEvCpIKL/A/iqiJSJSC7wLY/PFOz4+np/7+PypQDHZIm/7fzg+PsRkatlNNDdij6ZuoGngfki8hkRSbJuJ4oOztvB1yo/u30aKBWRr4lIiohkicjJPtbbDaSKDhgnAf8JpHjI9mkRKVJKudGuGSzZGq17T2V0F/BtEVlibZsjIlf7+9wh4OR38h3LilsC3AjYCQF3AT8SK6guOiHkMu+NRaRYRC6zlGA/2pXoDoPscYFRFPHLQ9Z9s4hssB5fDyQD29EnhIfRVy9+UUq9h77qm472AdvcDqQBTcAadNAy0nwIffJ7ltEA7Yser1eg4yi+eBXYBhwRkSZr2ZeBH4hIJ/BdRgOl/vi99X6bgQ8sOYbQrgQIfHx9vX8kuR3n38+JwHsi0oW+cv6qUmq/ZT1eiA5i16LdnD9h9ETu93hb214AfNTabg864Ou9Xjv6e/gDUIP+rXlmQV0EbLNkuwMdg+tVSvUAPwLettw6pyilHrPk+7vlbtuKTj6YLE5+J28Ae9EB+/9TStm/yzvQx/RFa/s16GC6Ny7gX9HHuQUd/PZpuVkW+TaP53eJyF0ez7dJnGVMiRWgMRhiinVF/A+l1GlRfM+PAHcppSqDrnwUIiIvopXKjljLYohvjKIwHDOISBr6qvhFdDDyEWCNUuprsZTLYIh3jOspzIjIc36CmeOKlQxRR9Cpmq1o19MOtCvCMMUQkTP9JQ7EWrajEWNRGAwGgyEgxqIwGAwGQ0COusZahYWFaubMmbEWw2AwGKYU69evb1JKFfl67ahTFDNnzmTdunWxFsNgMBimFCJy0N9rxvVkMBgMhoAYRWEwGAyGgBhFYTAYDIaAHHUxCl8MDg5SXV1NX19frEWZkqSmplJeXk5SUlKsRTEYDDHgmFAU1dXVZGVlMXPmTESCNkc1eKCUorm5merqambNmhVrcQwGQww4JlxPfX19FBQUGCUxAUSEgoICY40ZDMcwx4SiAIySmATm2BkMxzbHjKIwGAyGmNK0R9+mIDFTFCJSYY0n3G71X/+qj3VERH4pIntFZLOIrIiFrAaDwTBpnrgVnhp3mpsSxDKYPQR8XSm1QUSygPUi8pJSarvHOh8B5lm3k4Hf4ntoiGGCDA0NkZiY6Pe5wWAIE817ITkj1lJMiJhZFEqpOqXUButxJ7rlc5nXapcBf1aaNejh8gEnusUjVVVVLFy4kOuuu45FixZx1VVX0dPTA8Arr7zCCSecwNKlS7npppvo7+8H4LbbbmPx4sUcf/zxfOMb3wj6Hj/5yU9YunQpy5Yt47bbbgNg48aNnHLKKRx//PFcccUVtLa2AnD22Wfzta99jVWrVnHHHXeMe24wGMJMfyf0NEFXA0zBjt1xcekoIjOBE4D3vF4qQw83t6m2ltV5bf8F4AsAM2bMCPhe//XUNrbXdkxOYC8WT8/mex8NPBJ5165d3HPPPZx++uncdNNN/OY3v+HWW2/ls5/9LK+88grz58/n+uuv57e//S2f+cxneOyxx9i5cyciQltbW8B9P/fcczzxxBO89957pKen09LSAsD111/Pr371K8466yy++93v8l//9V/cfvvtAAwMDIz0xHrqqafGPDcYDGGm1WqjNNSrlUZqdmzlCZGYB7NFJBM9aexrSqkJncGVUncrpVYppVYVFflsfhhzKioqOP300wH49Kc/zerVq9m1axezZs1i/vz5ANxwww28+eab5OTkkJqays0338yjjz5Kenp6wH2//PLL3HjjjSPr5efn097eTltbG2edddaYfdtcc801Y/bh/dxgMISR1gOjj7saYifHBImpRSEiSWglcb9S6lEfq9SgB8DblFvLJkywK/9I4Z1iGijlNDExkffff59XXnmFhx9+mF//+te8+uqrYZUnIyMj4HODwRBGWqtGH3cdgcK5MRNlIsQy60mAe4AdSqmf+1ntSeB6K/vpFKBdKVXnZ9245tChQ7z77rsA/O1vf+OMM85gwYIFVFVVsXfvXgD+8pe/cNZZZ9HV1UV7ezsXX3wxv/jFL9i0aVPAfV9wwQXce++9I3GPlpYWcnJyyMvL46233hqzb4PBEAPGKIr6mIkxUWJpUZwOfAbYIiIbrWX/DswAUErdBTwLXAzsBXqAG6MvZnhYsGABd955JzfddBOLFy/mlltuITU1lXvvvZerr76aoaEhTjzxRL70pS/R0tLCZZddRl9fH0opfv5zrUeffPJJ1q1bxw9+8IMx+77ooovYuHEjq1atIjk5mYsvvpj//u//5k9/+hNf+tKX6OnpYfbs2dx7772x+OgGg6HlAOTOgLZDU9L1dNTNzF61apXyDsru2LGDRYsWxUginfV06aWXsnXr1pjJMFlifQwNhinNL1dA8RLY9Rycdiuc//1YSzQOEVmvlFrl67WYB7MNBoPhqMY9rC2J/NmQOQ06p57rySiKKDBz5swpbU0YDIZJ0FED7kHImwmZxVMyRmEUhcFgMEQSO5CdP8tSFFMvRmEUhcFgMEQSW1HkzdSuJ2NRGAwGg2EMrVUgCZBdri2KniYdt5hCGEVhMBgMkaTlAORWQEIiZBWDckN3Y6ylCgmjKAwGgyGStFZBnjVGOLNY308x95NRFAaGhoYCPjcYDJOgtUrHJ8BDUUytgLZRFFGgqqqKRYsW8fnPf54lS5Zw4YUX0tvbC8C+ffu46KKLWLlyJWeeeSY7d+4cWX7KKaewdOlS/vM//5PMzMyg77N27VpOO+00li1bxkknnURnZyd9fX3ceOONLF26lBNOOIHXXnsNgPvuu4+PfexjnHvuuZx33nnjnhsMhjDQ1w69LR6KYpq+n2IWRVy0GY8qz90GR7aEd58lS+Ej/xNwlT179vDAAw/w+9//nk984hM88sgjfPrTn+YLX/gCd911F/PmzeO9997jy1/+Mq+++ipf/epX+epXv8q1117LXXfdFVSEgYEBrrnmGh588EFOPPFEOjo6SEtL44477kBE2LJlCzt37uTCCy9k9+7dAGzYsIHNmzeTn5/PfffdN+a5wWAIA56psTBqUXQeiYk4E+XYUxQxYtasWSxfvhyAlStXUlVVRVdXF++88w5XX331yHr24KJ3332Xxx9/HIBPfepTQYcX7dq1i9LSUk488UQAsrN1v/vVq1fzla98BYCFCxdSWVk5oiguuOCCMUrB+7nBYJgknqmxAElpkJIz5VxPx56iCHLlHylSUlJGHickJNDb24vb7SY3N5eNGzfGRCbTatxgiDDeigKmZC2FiVHEkOzsbGbNmsVDDz0EgFJqpKX4KaecwiOPPALA3//+96D7WrBgAXV1daxduxaAzs5OhoaGOPPMM7n//vsB2L17N4cOHWLBggWR+DgGg8GblgOQlgepOaPLpmB1tlEUMeb+++/nnnvuYdmyZSxZsoQnnngCgNtvv52f//znHH/88ezdu5ecnNEfmu3C8iQ5OZkHH3yQr3zlKyxbtowLLriAvr4+vvzlL+N2u1m6dCnXXHMN99133xjrxmAwRBDP1FibrGI9vGgKYdqMxyk9PT2kpaUhIvz973/ngQceGFEisWAqHkODIebcsRymnwBXe8yCef7bsOHP8O+TGtYZdgK1GT/2YhRThPXr13PrrbeilCI3N5c//vGPsRbJYDCEwvAQtB+GJVeMXZ45DQa6oL8LUoKnvccDRlHEKWeeeWbQEagGgyGO6agG99BoaqyNnSLb3TBlFMUxE6M42lxs0cQcO4NhAvjKeILRorspNMDomFAUqampNDc3mxPeBFBK0dzcTGpqaqxFMRimFn4VRYm+n0IpsseE66m8vJzq6moaG6dWx8Z4ITU1lfLy8liLYTBMLVoOgCsJssvGLp+C/Z6OCUWRlJTErFmzgq9oMBgM4aK1CnJngCth7PL0fD2fYgpZFDF1PYnIH0WkQUR8DpQWkbNFpF1ENlq370ZbRoPBYJgQnl1jPXElQEbRlKqliHWM4j7goiDrvKWUWm7dfhAFmQwGg2HytB7wrSjAKrqbOq6nmCoKpdSbQEssZTAYDIaw09uqW4x7p8baZBYb11OYOVVENonIcyKyxNcKIvIFEVknIutMwNpgMMQcfxlPNpnTjEURRjYAlUqpZcCvgMd9raSUulsptUoptaqoqCia8hkMBsN4gioKy/XkHo6WRJMirhWFUqpDKdVlPX4WSBKRwhiLZTAYDIFpOaDv/SqKElDD0DM1PO9xrShEpERExHp8Elre5thKZTAYDEForYL0QkjJ8v36FBuJGtM6ChF5ADgbKBSRauB7QBKAUuou4CrgFhEZAnqBTypTXm0wGOIdf6mxNiNFd/XAcVEQaHLEVFEopa4N8vqvgV9HSRyDwWAID60HoPwk/69PMYsirl1PBoPBMOUYHoT2av+pseBlUcQ/RlEYDAZDOGk/DMod2PWUkgnJmVMmRdYoCoPBYAgnwVJjbTKnGYvCYDAYjkmCpcbaZBZPmZkURlEYDAZDOGmtgoRkyJoeeL0p1MbDKAqDwWAIJ61VkFsJriCn18yp0xjQKAqDwWAIJ4G6xnqSOQ3622GwN+IiTRajKAwGgyFcKAWtBwOnxtpMoRRZoygMBoMhXPS2Qn+HM4siy56dHf/uJ6MoDAaDIVy0Osx4gilVnW0UhcFgMIQLp6mxYFxPBoPBcEzitNgOdHdZZErUUhhFYTAYDOGitQoypkFyRvB1ExIho8hYFAaDwXBMEay9uDdTpJbCKAqDwWAIFyEriqnR78koCoPBYAgHQwPB24t7M0XaeBhFYTAYDOGg/TCgJmBRNIDbHSmpwoJRFAaDwRAOQkmNtckqAfcg9LVFQqKwYRSFwWAwhIORYrtQXE9To+jOKAqDwWAIB61VkJg6WkjnBHvdziMRESlcGEVhMBgM4cBpe3FPRqqz4ztFNqaKQkT+KCINIrLVz+siIr8Ukb0isllEVkRbRoPhmKFxN3TUxVqKqclAD1S9BaXLQttuirTxiLVFcR9wUYDXPwLMs25fAH4bBZkMhmOT+6+CF/8j1lJMTbY+An3tsPKzoW2XkgWJaXGvKBJj+eZKqTdFZGaAVS4D/qyUUsAaEckVkVKllLnsMRjCSW8rtB101nrCMJ61f4CiRVB5WmjbiUyJortYWxTBKAMOezyvtpaNQUS+ICLrRGRdY2Nj1IQzGI4aGnbo++Z9cZ/TH3fUrIe6jXDizfrEHypToOgu3hWFI5RSdyulVimlVhUVFcVaHINh6lG/Td8P91uFYwbHrP0jJGXA8ddMbPus+O/3FO+Kogao8Hhebi0zGAzhpGH76OPmvbGTY6rR2wpbH4bjr4bU7Intw1gUk+ZJ4Hor++kUoN3EJwyGCFC/HfJn68dGUThn499gqA9W3TzxfWQWa4Uz1B8+ucJMTIPZIvIAcDZQKCLVwPeAJACl1F3As8DFwF6gB7gxNpIaDEcxSukYxdKroKvRKAqnKAXr/gjlJ0Hp8RPfz0h1dgPkVgReN0bEOuvp2iCvK+CfoiSOwXBs0l4N/e1QvBgK50LTnlhLNDU48IZWqlf8bnL7ySzR93GsKOLd9WQwGCKNHZ+YtgQK5urMJ0Nw1t4Dafmw+PLJ7WcK9HsyisJgONaxM56mLYKCeTrrabA3tjLFOx21sPMZOOE6SEqd3L6mQHV2TF1PBoMhDmjYDtnlkJYLBXMABS37oXhJrCULidse2cyWmnY+NL+Is+YXsbIyj6SECF0Lb/gzqGFYddPk95VhpfQbRWEwGOKW+u06PgFQOE/fN+2ZUoqiobOPf6w7TGlOGr9/cz+/fX0fmSmJnDangLMWaMVRnpcenjcbHoL1f4I5541mik2GxGRILzCKwmAwxCnDg9C0G+ZdoJ/nz9H3Uyzz6ZnNdbgV/OmmEynOTuWdfc28sbuRN3Y18uJ2fQKeU5TBv1wwn0uPnz65N9v9HHTWwiX/FwbJLTLju+jOKAqD4VimaY+esGZbDymZkDV9yimKJzbWsrg0m7nTsgD48JISPrykBKUU+xq7eWN3Iw+8f4jvPrGNCxYXk5KYMPE3W/sHyC6DeR8Ok/TEfb8nE8w2jLC5uo3OvsFYi2GIJiMZT4tHlxXMmVKK4lBzDxsPt/Gx5eMtBRFh7rRMbj5jFt+9dDEt3QM8v3USQ4Ka98H+12HljZAQxuvszGLoNIrCEOf0DQ5z1W/f5d63q2ItSnyiFKz5LVSvj7Uk4aV+G7gSoXD+6LICq5ZCqdjJFQJPbtJdfT66LLBL6Yy5hczIT+f+9w5N/M3W/VEfrxXXT3wfvrDbeMTpMTeKwgBAbVsvA8NuDrX0xFqU+OSDv8Lzt8F7R9lIlIbtOiU2MXl0WeE86GuDnpaYieUUpRRPbKzlxJl5lOWmBVzX5RKuPWkG7x9oYW9DZ+hvNtirfwcLL9WN/MJJZrFuyNjXHt79hgmjKAwA1LTpvPn6jr4YSxKHNOyAZ7+pH7fsj60s4cYz48mmYK6+b47/Cu2dRzrZ09DFx5aPmz7gk6tXlZOUIPztvQl0yN36qFagJ06ir5M/4nwkqlEUBkBbFABH2o2iGMNADzx0ow7yLrz06Kpa7uuA9kNj4xPgoSjiP07x5KZaElzCxceVOFq/MDOFDy8p4eH1h+kbHA7tzbY8pI/NzDMnIGkQRqqzJxE/iSBGURgAqGk1isInz98GjTt0P58Zp04Zl4wj7GFF3vUSuZXgSor7nk9KKZ7cWMsZcwspyExxvN2nTp5BR98Qz2wOsRF1ezUUHzex4UTBMBaFYSpQbVkUnf1DdPcPxViaOGHLw7DhT3DGv8Dc86yqZY4eq6LBbt3hZVEkJEL+rLi3KDYcaqWmrZfLfGQ7BeLU2QXMLszgb++HGNTuaYaMwtC2cUpWfLfxMIrCAIxaFABHTJxCK4OnvgYVJ8M5/6GX2cVoLfGpKGrbeukIJb25fjskZ0HujPGvFcyNe0XxxMZaUhJdXLjEmdvJRkT41MkzWH+wlZ1HOpxt5B7WMyPSCyYgqQNScyEh2SgKQ3xT297LtCxtvtcf6+6noX54+EZwJcCV90BCkl6eVwniiruAttut+P2b+/nQT1/jJ8/tdL5hw3bdCNCXK6Vgrv6c7hD9+JGkvUbPywCGht08s7mO8xcVk5kSej3DlSvKSU508TenqbI9LYCCdP8WxV/XHOTvoVopNiJxXUthFIWBYbeirq2PVTPzAKg71hXFS9+Duk1w+W/GzgdITIGc8rhyPbV2D/C5P6/jR8/uQAEHmrqdbaiUrqHwzniyKZgLwwPQNomag3DzlyvgmX8F4O19zTR3DwStnfBHXkYylywt5bENNfQMOHC19jTp+wzfFsXQsJv/fWEXv3p1ElZY5jTojM8BnkZRGGjo7GPIrVgxQyuKY9r1tPNZXStx8pdg4SXjX8+fEzeup7VVLVz8y7dYvaeJH1y2hAsXFztPb+6s04H5aX4a/9nNAeNFKTbvg6ZdI9bckxtryUpN5OwFRRPe5adOnkFn/xBPbaoNvnK3pSj8WBTrDrbS3jtITVvvSAZhyJQshZoNMDQwse0jiFEUhpEf9pxpmWSnJh67tRRth+HxW6B0GVzwA9/rFMyB5v0xraB1uxV3vraXT969hpREF49++TSuP3UmxdmpNHQ4nLtcb7XuCGRRQPzUUux9Wd931NA3OMwL245w0ZISUpMm3rNpVWUe86ZlOnM/9TTrez/B7Fd2jLqM1h1snZhA8y+CgU44+PbEto8gRlEYqLYC2eW5aZTkpMbO9bT7BZ2CGCsev0X75K+6V7uZfJE/R48NtU8cUaaxs58b7n2f/31hFxcvLeWpr5zBcWU5ABRnpzrPWvOX8WSTUQQpOfET0N7zor7vbeWNrQfp6h/iModFdv4QEa47eQabqtvZWhOkItp2PfkJZr+8o4Ez5haSnpzAuqoJpk/POgsSU/X/IM4wisIwUpVdlpdGcXZqbCyKnhZ44JOw+hfRf2+AlgNQ9Rac9c3RNFhf2PMHYhDQXrO/mYt/+RbvH2jhxx9fyi8/uZys1KSR14uzrWQEJ99f/XbIKoX0fN+vi+jjEA+1FIO9ULUaMnRR2tsbNlOYmcKpcyafgXTFinJSk1zBU2W7rQsDH4piX2MXB5q6+fCSYlbMyGNd1QQtiuR0mPUh3cY8zno+xVRRiMhFIrJLRPaKyG0+Xv+siDSKyEbr9rlYyHm0U9PaS156EunJiZTmpMam6O7AG6Dco2M5o83u5/X9oo8GXi9GtRRKKW792wdkJCfw+D+dzrUnzUC8spWKs/VIznon7qeGbTrjKRCF8+IjRlG1Gob6YPmnADhYtZtLjy8lwTX5wrectCQuPX46T3xQQ1cgS6ynCVJzRjPgPHjZmndx3qJiVlbmsfNIR2hpyp7M/zC0VsWHgvYgZopCRBKAO4GPAIuBa0XElx38oFJquXX7Q1SFPEaobetlutVQrSQ7lcaufgaH3dEVYt+r+r5+e2yupnY9B4ULgk8sy7VTZKN7Aj3Q1E1TVz9fOmsOi0qzfa4zqiiCKPrhIWjc7d/tZFMwFzqqYcBhJlWk2POSdskcfw0AhcPNPluKT5TrTp5B98AwT2ys8b9ST7PfQPYrOxpYXJrN9Nw0TpyZj1vBB4faJiaMPePCvnCJE2JpUZwE7FVK7VdKDQB/By6LoTwxxe1W/PzFXext6Ir6e9e09Y503izOSUUp7QsPNwebu32nIioF+17X7Zv726Mfp+hr1wHEBRcFXzcxWReoRflKe4N14llRmed3Hceup5b9ulNpsFGndkA71nUje1/S/ZUsJb4go4MTKnLDtvvlFbksKs3mb+8dQvm7SOlu8ul2au0eYN3BFs5frCurl8/IJcElrJ9onCK3QrcJibM4RSwVRRng2cKx2lrmzZUisllEHhaRCh+vHxW8d6CFX766lyedpOqFEaUUNa29lOVpRVGao69Kw50iOzjs5pJfruaetw6Mf7F5n25Ot+QK/dwephMt9r4C7iGY/xFn68cgRXb9wVayUhOZW5Tpd53MlETSkxOCu56CBbJt4qE5YPM+rajmXUhjn9Cssjgpv3ec220y2EHtbbUdbKr2E9T2077jtV0NuBWcv0jHTzJTEllUmsXaicYpQLufDr2rK8HjhHgPZj8FzFRKHQ+8BPzJ10oi8gURWSci6xobG6MqYLh47AN9FR3tquj23kG6B4ZHLQrLfRHuOMXhlh66+oeobfeRY77/NX1/6j/p+2jHKXY/D2l5UH6is/UL5ujgdxRdZBsOtrJiRh6uAH55EaEkO5X6ziDfXf127T4rWhB4PTse0xRDRWGnxc47n2e31FGnCpib6rDtRghctnw66ckJ/HXNQd8r+LEoXt5RT3F2CsdNzxlZtqoynw8Ot07cfTv/IlDD+gImTnCkKERklpNlIVIDeFoI5dayEZRSzUop+/LoD8BKXztSSt2tlFqllFpVVDTxApxY0TswzLNbdHvhaBe7jWQ8ecQoIPyKYl+j9nO3dPsoJtr3qvb9ly6H7PLoWhTuYZ16Oe9C56Mt82dDf8doEVaEae8dZHdDJysDuJ1spmWnBL/YaNiuraKkwIN+SM7Qs6FjaVHseVHLmj+bZ7fU0ZU8jcy+8Le5yEpN4qqV5TyxsWbkPzGCUj4tiv6hYd7c3cS5C4vHKPATZ+bTN+hme+0EFVrZSh0PiSP3k1OL4hEfyx6e5HuvBeaJyCwRSQY+CTzpuYKIlHo8/RiwY5LvGZe8tKOerv4hCjKSo56aajcDtF1P+RnJJCe4wi7H/kYde2nt8coGGR6EA2/BnHN0Smbx4tFisGhw+H1t4s93EJ+wiXJzwI2H21AKR4qi2JFFEaB1hzcFc2NXdGenxc67AKUU22s7cOWWQUeAoPMk+OJZc1AK7n7D63vt7wD34Lhg9nv7W+jqH+KCxdPGLLdb4aydaJzClaAvXPa8qBMP4oCAikJEForIlUCOiHzc4/ZZIHUyb6yUGgJuBV5AK4B/KKW2icgPRORj1mr/LCLbRGQT8M/AZyfznvHKYxuqmZ6TyoVLSqKvKKyrJzvrSUQozkkJe9HdfsuiaPW2KKrX6WrUOefq59MWQ9NurUCiwe7ndBB97nnOt4lyiuyGg624BJY5CODqOph+/0HZgW6dfumvdYc3dhfZWGSi2Wmxcy+gvqOfzv4hkvIqtGIfCP/I3rLcNK5cUc4Daw/T4Klsu30X2728o57UJBenzRmrQIqzU6nIT5t4PQXoOEVfG1S/P/F9hJFgFsUC4FIgF/iox20F8PnJvrlS6lml1Hyl1Byl1I+sZd9VSj1pPf62UmqJUmqZUuocpVQIrTGnBo2d/by5p4nLTihjek4qrT2DoU/emgS1bb2kJrkoyBidmVySnRp2F9j+Jj8Wxf7XtL981of08+Il+uotWnnku56HytN1jrxTcmeAJETNothwqJUFJdmOuqQWZ6cyMOSmvdePom3YCajQLIq+duhu4kfPbOdP71Q5lnvS2GmxM08fyQbMnGa1RO+ITNLHLWfPYWjYzR88ky58tO9QSvHKjgbOnFfks43IiZX5rDvY6l9hB2POufoCJk7SZAMqCqXUE0qpG4FLlVI3etz+WSn1TpRkPKp5clMtw27Fx08oGwkkO+7XEwZqrBoKzyySkpy0sFs2doyirWdg7J9n36swfYUOJsNoymY04hQt+3WjuQUOs51sEpJ0y/EopI0OuxUfHGpjZWWuo/XtFFm/it5pxpON1Ryw8eA2/rD6AD96dgeHW8J/Ne+TvS/pC4ikNPY0dAJQVGZZcxFyP80szOBjy6bz1zUHR+NpPiyKHXWd1LT1jmQ7ebNyZh5NXf0cbJ7gsUrN1hcwcRKncBqjuEJEskUkSUResaqlPx1RyY4RHvugmqVlOcwrzqI4QqmpgahpHa2hsCnJ1q6nCV8NedHWM0BL9wDTslIYcis67QrY3jaoWa/jEzYF8/SVVDQyn3ZZV2uhxCds8udExfW0u76Trv4hR/EJcFCdXb8dktIhz2EuiuVm27JpHUqBAP8TysyLiWKnxc69AIA9DV3kpieRU1ypX4+QRQHwT+fMpWdgmHvftqwKHxbFKzvqEYFzFxb73MeJM3VrlAnHKUD/Lht36gy7GONUUVyolOpAu6GqgLnANyMl1LHC7vpOttZ0cMUJunykxGllbRjxLLazsd0Xbd5uogliWxN2kK+t29rvgTd12w47PgG6oK1wfnQUxe7noGihHvsZKvmz9Ykswr77DYe0n3vlDD89mbwozgryG2rYpj+zy+FfP7cS5UqioWobJ87M45az5/DMlrrJnQCdsOclfT/vfAD21ncxb1omkm1VZHdErihzXnEWHzmuhPvertIuPB8NAV/eUc+y8lyKsnw3j5xblElOWhLrJ9pJFnScAkYaIu5r7GJPfefE9zcJnCoKu8HJJcBDSqkgrRYNTnh0Qw0JLhlpRxBtRdE3OExT18A4RVGao5+Hy7KxM55WVuqTXUuPZdLvfw2SM8fXL0xbHHnXU187HHxnYtYE6CvtgS7oagivXF6sP9hKYWYyFflBUlktplmupwZ/3139dufxCQBXAv3ZleT3HuTjK8r54ofmUJqTyg+e2o7bHUElufelkbRYgL2NXcydlqVTetMLImpRgLYqOvuH+Mu7Vdr1lJim04XRx3ZTdTsXLPZtTQC4XMLKyrzJKdSCOdrC3v08w27FTfet5RsPbZr4/iaBU0XxpIjsRNcxvCIiRcAxOrQgPLjdiic21nDW/CIKM/WfOzstkdQkV9Sa8tW2jU2NtSnJsfzcYZJjX2M3SQnC8eU6YNxqK4p9r+rWDN6N1ooXQ/thfTKPFHtf1tXYocYnbKKUImsX2jmtRE5NSiA3Pcm3ku9q0FfHTjOeLA6o6cx2HeHipaWkJSfwbxctYEtNO499EJk4gWdaLEBzVz8t3QPMnWZVpWdP12NRI8hxZTmcu3Aa96w+wGBn41i30059cXCen/iEzaqZeexr7PZdO+SUBRdB1Wpe27yPg8097G/qDptLOBScKooNwIXAKuBbwP3Av0ZKqGOBNfubqWvvG3E7gZWaGoGMI394p8bajFRnh9GimJGfTpGlEFu7B7TbprVqrNvJxj6RNUSwbGbX85CW77wa25sCq3lgBOMUTV39VDX3OI5P2BRnpfqOUdjuvBAsioEhN2va85gp9eQka2V12bIylpXn8NMXdjobIxoqdlrsvNH4BMC8EUVRFnGLAuDWc+fS2jNIXV31GLfTKzvqKc9LY0FxVsDt7TjF5NxPF8HwAB+88TgAnX1DYXMJh4JTRfEdpdQh4FTgfOAO4OcRk+oY4JENNWSlJI4zX0OaUjZJar2qsm2mZaUiEj6LYn9TN7OLMslL1ym4rT2DsM9q2+EZyLaxT2SRilMMD41WY7smOCEtZ4YOukcw82mDdYIJWVHkpPp2PdnuvBAsitd3NbB9sJhEhnQ/LrRb5bsfXUx9Rz93vRGBz7/nJe3qqTxDP7UVRbGnooisRQGwYkYeZ8wtpLOlnuE0rSh6B4Z5a08T5y8qDmrlLS3LITnBNfFBRgAVJzOcnM2Mxrc4aZZWPIeilXXmgVNFYSf2XwLcrZR6BkgOsL4hAL0Dwzy/tY6Ll5aOy8GORA2DP2pae3EJlOSMrZ1MTnRRkJESlljJ0LCbg83dzCnKJCs1kQSXaIti/2uQUzHaeM6TnApIyY5cnOLwe7qYyUm3WH8kJELezIi6njYcaiMpQUYm2DmlOCvF92+ofrueXJfpvM3NoxtqaEm1Ou149HxaWZnPpceXcveb+yY+I9ofe1+CWWdCkv5d7q3vJDMlcSSGR/Z06G2JSNGdN7eeO5fs4XaqevV7v723if4hN+cv8h+fsElNSmBpec7k4hQJSWxKXcV5CR/wbxfqVOWDcawoakTkd8A1wLMikhLCtgYvXtx+hO6BYa5YMb5ZbkmOVhTR8ENWt/VSkp1KUsL4r7JkItXZ1evHdbysbu1lcFgxuygDl0vITUuivbsH9r8Js8/WbTu8EdFDdSLVymP3c+BKgjkhVGP7In+2np/tAKUU//S3Dfz53SrHu99wsJUl03NCngtdnJ1KY2c/w97B5votzusn0GnNr+ysZ/FxVos1r55Pt31kIW4FP30+jOmyXmmxoAPZc6Zljl7B55Tr+8668L2vH06elU9RQidrG1wMDLl5eUc9WSmJI1f3wVg1M48tNe0TLqKt7+jjry2LKJR2jhP9WzvUHP35IE5P9p9At9r4sFKqDcjHpMdOmEc21FCWm8ZJM8f/2KZlpYQ1NTUQNa294+ITNiXZIRbd9XXAHz8M938ChkaDd/usjKc5RTpjJC8jmayWrXruhK/4hM20xTqVMxIKc9fzMPN0XdQ0GfLnOE6R3VzdzjOb6/jFS7sdnTQGhtxsqm4L2e0EuujOrXQQeITG3VC3abQC3gFPba5jcFhx0UlLdOW6V8+n8rx0Pn/mLB7fWMsHh8LUEtsrLRZgj5UaO4KdIhuFuSUy1Eeq6uNQXzqPbKjmlZ0NfGhBEcmJzk6dqyrzGRxWbPbXvjwIf363iteGj0eJi9T9LzEtK2XiRXyTwNGnVUr1KKUeVUrtsZ7XKaVejKxoRycNHX2s3tPI5SdM99kyuiSKRXe17b3jMp5G5fDjvvDHoTW69Ub1+/DSd0YW2z2eZhfqP3peehKz2t8DRFsU/iheorOewh20bN6nT3hOZ08EomAODHZD55Ggqz647jAu0fGZRzcE969vr+ugf8g9QUXho+hu7R8gIRlW3OB4P49uqGZhSRaLp+foNE0fXWRvOXsuRVkp/PDp7eGxgr3SYtt7Bmno7PdSFJYlHoWAtl2VnZZXzH8/s4PGzn4ucOB2srG/v4m4n3oHhrn/vUOcuGguUnEy7H6eGfnpce16MoSJJzfV4lZwxQnlPl+PVi3FsFtR19Y3LpDtKUdbKH2nqt7UJ6KVN8J7d8FW3XB4f1MX+RnJ5Fm9pPLSk1nUsx6mL4f0AOa73coj3AFtu3fOZOITNvbY1CAB7d6BYZ7aWMvly8s4riybe1bvD1qDsH6CgWzwkbXW3wWbHoDFlzuOT+xv7OKDQ218fEWZdvkUzPU5lyIzJZFvXriADYfaeGrzJF1BI2mxF44s2tuoC8xGAtkwalFEIaBtV2WfcfwCOvuHSHAJZy9wHuPJz0hm7rTMCWU+PfpBNW09g9x8xixdfHdkM0uzuzkUrxaFIXw8sqGGZeU5oznhXjieezxJGjr7GHIrv66nkAcYVa3WqaYf+SlUnAxPfAUad7GvsZvZhRkjq5WkDLJweCfM9pHt5Mm0RZagYVYUu56DokU6ED1ZCpzVUjy7pY7O/iE+cWIFnztjNvsau3ljT+ABWxsOtVKWmzbyPYTCuN/Q5gd1q+yTvuB4H49/UINL4LLl1tV74VzorNVKx4srV5azZHo2P3lu5+QaWo6kxY66nfaOpMZ6pKImpenU5qgoCm1RLJs/hyXTszltTgG56aHl8ayqzGNdVUtIBYput+KPqw9wXFm2jodYFvAZagNHOvqi2jgUjKKIKjuPdLCjrmNM7YQ3dmXtkfbIpsh6z6HwJqTq7L527f+eeYZuwXH1ffrP/OBnqGtoZHbRqKI4fngLibhRvtJiPUnL0y6GcAa0e9v0iMlwWBOghyy5koLWUjy47jAzC9I5eVY+Fy8tpSQ71fdIWA82HGwNOB87EIWZybjEqs5WCt7/PZQug/JVjrZ3uxWPflDDGfOKRhXVyPzs8Z81wSV859LF1LT1cpf3LIdQ8EqLBR2fSE1yjbd8c6JTS0G3tihcmUX8/QuncNenfc5OC8iqmfl09A2NpPk64Y09jexr7ObmM2Zpi65oAeRWsqrhEabTRHVrdK0KoyiiyGMbakh0CR9dNt3vOimJCeRnJEc8RmEX25X7cz2FUp19aI3u2TTT+oNnT4er7kE17+GbA78ZY1Es7F5Lj0qhZ5qDP1y4W3nY1djhiE+AoxTZA03dvH+ghU+cWIGIkJzo4vrTKlm9t4mdR3xPQKtt66WuvY+VM3InJFZigovCzBQdozj4NjTu0NaEw+rutVUtVLf2cqVnVl6BTs30N+3ulNkFfGzZdH716l7ePzCBdFClYM8LY9JiQddQzCnKHB/Pyy6LeHU2MNrnKaOArNQkMhy0evfmRKvH2bqDzo/LH1cfYFpWCpcstc4VInDBD8jsOcwLKd9i8P17ozojxCiKKOF2Kx7fWMPZC4ooyPTdSMxGF91FR1EEdT05kaPqLR2f8Kxynn02dSu+zscS3uXs9sdGFs9oe5817kW0Djg4aRUvhsZd4RtitPt5XWHr8MraEQVzAqbI/mPdYRJcwlUrRmNSnzppBmlJCX6titH4hLMUTF+MVPi/f7e2zo670vG2j26oISM5gQsXl4wutOMxAeZn/+iK45iRn85XHthAU1eIFnH9Vl2pv/CSMYv3NnhlPNlkT49ejEISIDV3wruYkZ9OYWaK40FGu4508taeJm44bebY7Koll9Px2TfY4p7FonXfgb9cDm2HJixXKBhFESXqO/uo7+jnrAWB+8OAbvMdcYuitZfcdP9XSFmpSWSmJDqzKA68pZWE1wzmd0uv5+XhE1iw6X/0yNG2Q2R3V7HavZTWbgcn/2nWEKNwzGze9HfY8ZRuieCnGntw2B16o7sAKbJDw24eXl/NOQuKmOYRa8hNT+bqVeU8sbF27CQ1i/UHW0lLSmBhaeAWEYEozk5hqK0GdjwNJ3wm+Hxsi77BYZ7ZUjfS12mE5HQ9sGnfq3rOuA+yUpO481MraOsZ5Gt/3zi+jiMQO54GBBZcPLKou3+ImrZe5vlqlZFdpovuBsNc7OdNd5O+uHBojflCRDhxpvMGgfes3k9qkovrTp4x7rXcsnl80fVdnq74pp4O+ZtTYe09EbcujKKwcLsVv3hpN89uqWN/Y1doP3IH2Cfc6TnBg5MlOakRD2bX+mgv7k1xdkpwRdHbBkc26+Z+Xuxr6uHfhr+s/9T/uAG2PATAm+6lo40BAwoQhlYeA93w+JfhsS/qAUnnfc/naoPDbs772Rv874u7Qtt/wWwY6vVZ/PXarkYaO/v5xKqKca/dePosBt1u/rpm/BXhB4daOb48x2chpFOKs1M5q/Np7RI88WbH2724Xc9v//gKH1l5H/omHF4Dr/3I7/aLp2fzg8uWsHpvE796NYQphTufhhmnQObohdRoDY4viyJKKbI9zWMaAk6UVTPzqW7tDfp/aurq5/GNtVy5otxn0FxEKM/P5BHXhfDld7V1/My/wp8/pi2yCGEUhcWRjj5+9eoevnz/Bs792Rsc970XuOzOt7ntkc3c9/YB3tvfTPskiuDsH4iTLJbi7FSaugYYGHJP+P2C4WsOhTd2lXhAvOMTHuxv7CY3fxpyzV/0H+6VHzCUUcJeVeZMURTO12b/ROMU9dvh7nNg49/0Se6GpyDLdw7867saOdTSw1/fPRhaozu7i6yPgPaDaw9TmJnCOQvHW5GzCjM4b2Exf11zcEwGS+/AMNtqOyaUFutJaWYCV7hfYnjehSFleD26oZqy3DRO9lV5vOJ6XYfx1s+0deaHT6yq4OMryrjjlT28FSS7C9CDeeq3wsJLxyzeU+/V48mTHFtRRNj9ZFsUk2RVpbM4xV/XHGRgyM1NZ/ifkVJp11LkzoDPPA4fvQNqPoDfnKYTFyJgXRhFYTE9N43tP7iIp249g59edTzXnjSD9KQEnt92hO8/tZ1r7l7D8h++yIvbghdX+cI+4Xr3VfLFyEhUH26JcKCUCliVbeOoOrvqLUhI8dmFdX9TF7OLMnXGzSU/A2B45lmA1e8pGIkpehRnqJlPSsH6P8Hvz9EtRT7zGJz7nzr47IeH1h0mJdFFZ/8QT2wM4Sp1pJZirKJo6OjjtV0NXLmyzK9l8LkzZ9HSPcDjHu26N1e3MeRWk1YUK3reokjaaVl0veNtGjr7eHN3I1ecUOazGBSAi/8XylbCY7foam8fiAj/7/LjmDctk6/9fWNwq3Tn0/p+kZeiaOgiKUGozE8fv41tUUQ6oN3TFBaLYvH0bNKSEnhuyxF21HX4vBjpGxzmr2sOcs6CIt9WlEVlQTrVLb3a6yECKz+rrYsZJ0dsxnboIfyjGLuJ19Ly0SZsSikaOvvZXtfB5/60js3V7Vy4pCTAXnxzpKOPpAQh30EOdolHZW15no8/ySRp7x2ke2CYcj+psSNy5KTQYPUMSvB34rDrJ5LGKsBht6KqqYdz7JjMis9AQjKJ5SchG7bT4tQ6m7YYatY5WxegvxOe/hft5pp1Fnz8936tCJumrn5e3dnAjafP5K09Tfzl3YN80spSCkpOuQ7ke1kUj2yoYditfLqdbE6elc+S6dncs/oA11jvt95qhXHCjMkpisWHH+SAu5iW/FNwWh725EarGNRHD7IRElPgE3+G350FD14Hn3vFZyuU9OREfnPdCj7267f5ygMbeODzp5Doz5W242koXjrO8tnb0Mnswkzf22WV6vtIWxQ9zWGxKJISXJw2p4BnttTxzBbtpizKSqEyP53KggwqC9Jp7RmgqWuAm8+YHXBfMwrSGRh2U9/RN3qxl1sBn35UD9OaRDzFHzG1KETkIhHZJSJ7ReQ2H6+niMiD1uvvicjMGMhIcXYq5yyYRnHWBBrlWdS39zEtK9X/lZoHfovu2qv9BhJDocZPe3FvSrJTGXYr/xksI/GJ8W6n6tYeBobdY2ooWHYNCQWzyElLos2J6wl0hXbbId1LKhh1m/UJbOsj2oL4zGNBlQTo4rIht+LqVRV85tRKttd1sOFQmzP5XAl6/rRHdbZSiofWHeakmfkBrwxFhM+dOYs9DV28uUenYW442MbswgzyMybRnLluM7lN6/nr8AXUdzp3lz6xsZZlFbkBZQa0crz6Pq0cH7/Fr6tj7rQsfvzxpaytavUf++lq0N18vawJ0BlP/gpTSU63iu4iGKMYHtIWafrkLQqAO69bwZO3ns6vrj2Bb354AWfPL8LlEt7e28TPX9rNvW9Xsag0m9PnBlZMlfn6PzWu55MIpEw8ASIQMbMoRCQBuBO4AKgG1orIk0opTz/DzUCrUmquiHwS+Am6g21MKM1No659YlkWRzr6HLmdwKPfk6dSOrIV7j4Lln4CrvjthGSwsYvtgrqe7KK79j7fsZVD7+r4xKzxgWy7x5Ovk05+erLzqV92K4+GHdq09kdHLdx3iR5XecPTuumfA5RSPLy+mmUVucwvzqIsN40fP7uTv6456Nz9UzBnjKJYW9XK/qZuvnyOjxbqXlyydDo/fnYnf3hrPx+aV8iGQ62c6yOmERJrf49KTOOhvg9R5vDCpm9wmO11Hdxy1hxn7zHrTLjwh/DCv8PqX8CZvueYXba8jPcPtPC7N/ZzYmU+53uPD935DKDGxSf6Boc51NIzWhnui0jPpei14glhcD2B9lgcX57L8eW5417rGxzmcEsPhZkpQS3ZygLtZTjU0s2pcyZv7TghlhbFScBepdR+pdQA8HfgMq91LgP+ZD1+GDhPnM6EjAAlOakTHuZT39HvWFHkpSeRnOAatSiUgme/qYvFNv1Nt6GYBDV+RqB6Y7vA/FpRVat1fKJsfF2CnbEy24eiyE1Pct4d126LHayVxwv/AUP98NlnHCsJgK01Hew80snVK3WWT0ZKIleuKOOZzXVju68GIn+2VhRunXzw4NrDZKYkcvHS4C7K5EQXN5ymXV4vbKunpXtgcvGJ3lbY/BAc/wl6E7Kodxjn2nWkk2G34riyEDrqnvJlXZ/x6g9h7yt+V/vOpYtZMj2brz+0icPeDe12Pq1dTvYFgcX+xm7cyk8g2ybStRRWQ8BwuJ6CkZqUwLzirJGeaIEozUkl0SVR7SIbS0VRBhz2eF5tLfO5jlJqCGgHxn1rIvIFEVknIusaGx1kWUyQ6Tmp1Lb3htwlUynFkfa+0cErQRARpmV7DA7a8hAcegcu/j8oPg6e+ir0THwYSm1bL6lJLgqC/CiLrepsvwHtqreg4qRx8QnQU+1y05N8ulDyM0KwKHJnQHJW4ID2vtdg26P6qrbA4RWxxUPrdRDbs1r+06dUMjDs5h/rHLaxLpijexR11tLRN8gzW2r56LLppCc7M9ivO1kX4P3n41uBiTUCHGHj32CoFznp80zLcj4tcVutdu0tmR7CkCQR+NivdO+sR272m56ZmpTAb65bgdut+K+nPBR+Xzvsf0NbE17Xf3sarGaA0wK4UnIiXJ1tNQQMl0URLhITXJTlpUW1i+xRkfWklLpbKbVKKbWqqMh5Z8dQKclJo28w9FkRHX1D9A4OO1YU4DHprq8DXvxPXQOw6ma4/Df6B/z8uJCOY2radMZTMOOsMCOFRJf4TpHtbdMxAR/xCdDdRz1bd3iSm57sPEZhDzHylyI71A/PfkPHCU7/mrN9WvQNDvPExlo+vKSEnLSkkeXzirM4ZXY+97930Fk9jZ351LyPpzbV0jfo5poT/QexvclNT+bKlWU0dfWTlZrI3GAxAn+43To9csapULI0pHqcbbXtZKcmBk1wGEdyBnzyr9oF+eCn/U6dqyzI4NOnVvLarsbRrgN7XtIFlYs+Nm79vQ1dJLiEmYUBkjnsSXeRKrrriZ5FESoz8tOj2kU2loqiBvD8N5Vby3yuIyKJQA7QHBXpfGAXy4Ua0B6poXDoerLXre/ohzd/Cl312ppwuXSq6Ye+qTuC7ng6JDlsalqD11CAno1cnJ1Kva/Pe+hdQPlVFPsau/0GRfPSk2hxqihAF97V+xli9M6vdOX2xf/n07IJxEvb62nvHeTqVeOLyz5zykyqW3t5Y3dD8B3ZtRQt+/nH2sMsKM5iWXkIV+bATafrvPkTZuQ5Snjwyb5XofUAnPg5wCqYdKgottZ2sHh6trNML2/yZ8PH/6DjaC982+9qV60sZ9iteMxOB97xFGQW+0yt3tvQRWV+OimJAab7ZVvfW6QC2iOup/iyKEDHKQ5GcdJdLBXFWmCeiMwSkWTgk8CTXus8CdiTVq4CXlXRmBHqh5IRRRHaFcxIDUWIFkV6+17Umt/qFgzlHk30zvw6lCyFp7820t0yFJwU29kUZ/vJ9DrwFiSm+oxPdPYN0tjZ7zM+AXrKXd+g23mr5GlL9Ixr7+rn1oPw5v/pK1KP1tROeWh9NdNzUjltzvgTwYVLiinKSuEv7x4MvqPsMkhMpfnwDjZVt480AAyF2UWZPHrcO/zH7El0X33/d/rEa12hO3U9DQ272VnXEZrbyZv5F8Jpt8L6++Dguz5XmVOUyYoZuTy8vho10KMtigUX6wsgL/YEyniyifRcCtv1FGhuSoyozM+go29oUkXAoRAzRWHFHG5Fj1jdAfxDKbVNRH4gIrYteg9QICJ7gX8FJu5vCQN2llCoFoV9RR6KoijOSubb/BGSMuD87499MSEJLv+tdv88F9pE2r7BYZq6BhwritIcP0V3VW/5rJ8Aj6l2Rb5dT3lWLYmj6mzwaOXh5X56/jYQF1z0Y2f78aCuvZe39jRy5cpynzUiSQkurj1pBq/vbgxu4rtckDeLxqrtJCe4AraR94tSrDh4LwvWftev+yYgtRthz4vamkjUx7ckJ5Wu/iG6+gNXmu9v6qZ/yM2S6ZMcDXv2tyGnQtex+GnkeNXKCvY0dHFg7bN6OqCPtNiBITdVTd2BA9kQ+TYe3U26GWBCUtBVo80MK/PpYEt0rIqYxiiUUs8qpeYrpeYopX5kLfuuUupJ63GfUupqpdRcpdRJSilnk+wjRGGm9tlP1KKwZ004YXnXm5yRsI3Gk77pO5hWshTO+pauGdj2uOP91gbpGuuN3YV0jCHX2wpHtvjs7wS6IhtG52R7YysKxwFtX5lPO5+FXc/C2d/Sef0h8uiGGpTS7hB/XHtSBS4R7n8/uFXRnFJOYtsBLlxSPLEaiL52feLsboB194S+/ev/o09qJ39xZFFxdpBkBItttXqe86QsCtDxiov/V7c1f/dOn6tcuqyUlEQXresfgZQcmDl+hvfB5m6G3CpwIBsiPzs7TMV2kcBOkY1W5tNREcyOFgmWzz7kGEVHH3npSaQmBfC3ejLQzbJtP2W7u5Jd5Vf5X++Mr0Hpct0UrMtZtpfT1FibkpwUegaG6ejzuCo9GCQ+0dBNgkuYke9PUegrNMdJAen5uhLXtigGeuC5b0HRQp2iGSIjBXGz8qks8C0jaGvqgkXF/GPt4YBusue31vHowVRmSD3fvXRhyPIAo+6TpAxYfbvPSXJ+qVkPu5+D074CqaMn++IsZ9MSt9V0kJLo8qvYQ2LBR3QW0+v/o12DXmSnJnHxkkJmt7zF8NwLRqwfT+ypdkFdT8npuoV6pCyKMLXviAQz8u1aCqMo4pKSnFTq2kJ3PdnFa4546+ek9NTxncHPcqQzgNvAdkH1d2pl4SB8MzLZzqFFYcs95mRTtVrHJ/zMddjf1MWM/PSxvfQ9sK+4HVsUoK0Ku4vsWz+D9kNwyc8n5BZYd7CVquaekdqJQHzm1EpaewZ5dovvedD/WHuYL9+/gaHcmSQzyDT3BHMt7JPdOf+uT1Brf+9829d+rKuUPawJYKS1ebA4xbbaDhaWZPlvsREqH/mJdgk++02fv8mbKurJo5MPMnxfaOxp6ELET9dYb7LLI+h6ao7LQDboFimFmSlRC2gbRREipU46qnpxpKOPEqdup+Z98M4vGV56DevVguDpjcWLtW94x5O6liAItW29uMRZc0IYjauMKTSselPXTyT6/kz7veZke2O3T3acIgu6IKtpl67QfvsOOP6TIRXWefLQusOkJydw8dLSoOueNqeA2UUZ/GXN+KvjP7y1n397ZDOnzy3kxo+epxcGmZ/tF9t9svgymHuB/oz9ncG3O/w+7H0JTv/quPYN9ncc6DeklGJbbTuLJ+t28iSnXCu8PS+MNvzzYEnHm/STxO9qffc02tPQRXle2th5GP7Ing4dkXI9NUFGfLqewM58MhZFXFKak0ptW2hFd/VO23copV0qCSkkXPgDctKSdIpsME77Z93R85mvQ2d9wFWr23opzk51POtgnKLoadFpkH7iE8NuxYGmbr+BbNCV2QAtToYX2RQvgeEBePAzkJSu20dMgJ6BIZ7ZXMclS0sdjbUUET59ciUfHGpja4325Sul+L8XdvH/ntnBJUtL+cMNq0gtnq83CDI/2y8dtfoqPKsEzvm2jgO997vg2732I33Ve9Lnx72UmZJIRnJCwAub6tZeOvqGJh/I9ubkL+lGf899a6zCUwrXrmepzjuFV/Z1+ex0sKe+M3h8wiZ7emQsCqWsGEV8WhSg242Pq3SPEEZRhEhpThr9Q86L7vqHdJaRkzkU7H5eXx2efRtklYwW3QUjIXE0C2rDnwKu6rSGwsYOwI/IEaR+oratl/4hd0C3QVKCi6zUROdZTzAa0G7eA+d9Z8yAm1B4dssRugeGuTpAV1dvrlxZTmqSi7+uOYjbrfjOE1v59Wt7+eSJFfzy2hN0rn9WKSSmjen5FBIdNTq1NSFJK/35F+kakb52/9tUvQ37X4cz/kUHkn2gx+r6v9gYrcgOs6JISIRLf6FP4q95ZKXVbYT2w+SccAVuBY9sGGsNDLsV+5u6fY8/9UVOmT6hD4a5JX9fu26ZE6fBbNCZT3UdffQPTb5RaDCMogiRUssyqHWY+WT/SYOmxrqHdbpn0cIRX3NxKJPuihZAXiU07gy4Wm17r+NANuj2C/kZyaOKwo5PlK30uX6gHk+e5Gckh6YoihaAK1EH71fd5Hw7Lx5ad5iZBekjA++dkJOWxOXLy3h8Yw23PrCBv645xBfPms2PP750NLXW5dLHv81B3YUvOmpG0z1BuxP72mDNXf63ef3HWrkEOB5jWsH4YHttOy6BhSVhVhQAFSfqWQnv/RbqNullO54GcVG48nJOmpnPI+urx1jnh1t6GBhyM8epohhJkQ1zLUWctu/wpLIgHaXgcEuEx8FiFEXIlOaOdlR1gv0nDVqV3VGje+Wc/MWRAG2Jk1GknhTOhybfw2RAX63VtfU5To21Kc72aIZo93cKEJ8A/zUUNrnpybSGUiyUmALX3K9nIfiZeR2Mg83dvHeghatWlodcEPfpUyrpG3Tz7JYjfOuihXz7I4vG7yOnYuLD7ttrRtM9AaYv19lD796pLUVvDrypv4sz/lVn//ihJDs1YGPAbbUdzJ2W6SweMBHO/56+Kn/6X/TF0M6nofJ0yCjgqlXl7G/qZoM1gwN0fAJwblFEqpYijquybeyswkNRqKUwiiJERi0KZydwx1XZdkM1ux0E9kjUfoaGHY5ELZwPTXtHuph609DZx5BbheR6AiuA397nEZ8Yn/tus7+pi+zUxKANB/PSk5xNufNkwUX6qn2CPLy+GhF8z4MOwnFlOdx6zlx+/oll3HK2n8aDuTOg7bDv1wKhlD7RedeDnH0b9LfDmt+MX/+1/4as6fqKPQDF2boVjL+Y2rbaSVZkByMtDz783zqF98XvaIvXail+8dJS0pISeHj9qPvJbgYYNDXWJmIWhaUo4jyYDdGppTCKIkTsorsjDl1PR5xWZduKwmPKV3F2Km4FTV0OT6iF82CoF9p9n6xGUmNDbPymTzZ9cPAdAsUnQNdQzJmWGfSKPT89RNdTGHjsgxrOmFsYskVl840PLwisZHIrdJO6UGogQLuYBrvHWhSgiyoXfQze/c3YbsH7X9Oxog99PWh/q2nZqQz4iak1dfVzpKMv/PEJb5ZeracNrrGK8BZeAmC1Yi/lqU119A5oP/ve+i5Kc1LJSnWY9hypNh4j7TviV1EUZCSTkZxgFEU8MlJ057CWor6jj5RE10imj19aq7QP3sNPPZJx5DROUbhA3zft8fmy08l23pRkp9LcPcDQgTd1wLZshd919zd1Mbsw+NWg7iAbnT41oE+K1a29nDU/ct2FybEC5H4UtV9st4lnjMLm7G/r8Zbv/lo/t62JnArdAywII2N1fbif7ED24kgrChGr5iVZx5hyRxMJrlpZTlf/EC9Ys+j3Njro8eRJpIrupoDrSUSoyE+PStGdURQToDTHeXX2EWtgUVCfeGuV/vMnjKZs+px0F4hCK0XTT5xioorCdre5968OGJ/o6h+ivqM/aHwCID8jia7+IQaGHLrVJsmuI9qlEZGgrU2u5RYLNU5hz1TwpSiKF8OSK3SqbHcz7H0ZqtfCh77h93sYs7mdtebjNzTSuqM0gq4nm8K5cO0D8NHbxyw+eVY+FflpPLT+MG63Cjz+1B/ZEZhL0dOs07ADxH/igcoCoyjilpKcVMf9nur9jRH1prVq3HB5v7Oz/ZFRoCt0m3zPJ65p7SU3PclR/cAYOXJS+ZBrE8lN22DueX7XOzAy/jS4ophQ0d0k2FGnr54XlkZmpjAweqUcqqKw3SY5fpoJnvUtGOiGd+7QdRO5lbD8Oke7Lg5Qnb2ttoPyvDRyglm74WLu+TD9hDGLXC7hyhXlvLOvmbVVLfQMDDuvobCJxKS77qa4tiZsKgsyONTSg9vJ3JRJYBSFUwb74KXvQuMupuemUdfe56joTldlT0xRFGQkk+gS54oCdBqpH9dTbQjtxT2ZntzLT5PupiNrDpz0Rb/rjTYDDH5FONLGI0qKYteRTgozUyjMdN6YMWQypunRsBNRFOKCTD+jU6cthKVX6bqK2g/grH9z3LqkKMt/Y8DttR2Rj0844MoV5SgFP31BX+AE7RrrTSRmZ8d5VbbNjPx0BobcjkfeThSjKJzy3Dd1W4WnvkZpdgr9Q+6g6Z1KKa0ogqXG9ndqU9dLUbhcwrQsZ8NnBofdfOr3a3ipMYfO6u387o19vLKjnqqm7pEpbfZku1CpfO97FNDBywt/GDB4uq+hC5eMtkAOhB2zaQ2lOnsS7DzSyaJIWhOgaylyyicWo8gsGeN2HMdZ39L3ebN0+xKHpCYlkJeeNO5E0tU/xIGm7shmPDmkIj+dU2cXsP6gTpMNecJfdgSK7uK4c6wn0cp8Cs0Hcazywf2w4c9QcjwceoelFWuBTOraewO2lG7tGWRgyB3c9WR32fSR+um06G5dVSvv7Gvm7NzpXDDcym+eW0c7+g+XnOBiVmEGB5q6fQ7pCcjWR0je8Sh3qE/QoWbx8QCr7mvqpiLYVDIL+7hFI/NpaNjN7vpOrj914qm1jplIimx79fiMJ28K58EVd0P+rMAKxQe6Dmas68l2xcWDRQFw9apy3t3fTGFmMnmhtmm3XXadtaNjaSdLd/NockgcU2nXUjT3cMrsyCk2Y1EE48hW3Zl15plw84uQN5MlO25HcAfNfJpMaqxNiWexWwBe3VlPcoKLz3z0AgDevrmcR245jZ9eeTw3nj6T8rw0ZhVmcNaCELJ+Oup0/6iylTyZ9cmglk2wZoCehDy8CPj+k9t4c7ezduqeVDX30D/kjmwg2yZ3AkV3HbX+4xOeHH+13469gZiWnUqDl0WxrSZMMyjCxEXHlZCZkhh6IBs85lKE0f0Uxy3GPZmem0qCSyI+wMhYFIHoa4d/fEYPhLnqj5CUBuf8B2mPfp5LXWuo6zg+4Oa2JVCSE8QvHkBRFGen8taepqCivrKjgZNn55NWqnsgZXbuZ+WK01lZ6bxVxRiUgie/os35K35H0aNNARXWe/ub2dfQxRlznV3VjLqenCmKrv4h7nuniiPtfXwoxBTXnUeiEMi2yZ2hhw8N9urfSzCU0v71eRdETKTirBR2WcfAZlttBwUZySNZUbEmPTmROz65PHgauS/CXZ090AODPVPC9ZSY4KIsNy3iridjUfhDKXjin7Rb6Op7R5vQHXcVqngJ30h6iPqWwG2gR6qyg82iaK3SA2fSxp/Ui7ODj7Pc39jF/qZuzl9UrE9UCSkBW3k4Yv19ukHhBT+AwnmU5qT5VRT/WHeYT9/zHuX5aXz29FmOdp+SmEBGcoLjNh5VTfqKaePhNkfre7KzrpMEl0zsajVUcmboe6dT1/ra9EnJV2psmCjJSaWxs38kVgVaUSyenh1yK5NIct6iYlZWTmA+9UjRXZjajU+BYjtPKgsi30XWKAp/vHsn7HgKLvgvqDxtdLnLhZz3PSqlnspDgec/HGnvQwSmZTmwKHxYEzBqjQSKU7y6swGAcxdO032QCuZOTlG07IcX/gNmn61nMGN1Ie3sG5OGN+xW/PjZHfzbw5s5eVYBj91yekhZVbnpyY4tiv2WojjS0Rda/yu0RTGnKMNR7GTShJoiO1JDESRGMQmmjVT46zjFwJCbPQ2dceN2mjTJGdrqD5dFMdK+I/5dT6Aznw4aRREDDr6rU2EXXgqn3jr+9XkXsjNpMec13KvNVD/Ud/RRkJESfPZDAEUxUksR4OT48o56FhRnUWGNR6QocHPAgLiH4bFbdJX4ZXfqTB50g8LBYUWzdWLv7h/iS39dz+/e3M+nT5nBvTeeGHI+figdZO0aDQjdqthR1xmd+ARoiw6cKwr75DaBud9OKfZKkd1d38ngsIqbQHZYyAnjpLtu26KYGoqisiCdtp5B2nsjl0FoFIU3XY3w8I06A+ny3+j2A96I8FzJl8h3t8D7/ofL6NTYINaE261bU/uzKIK08WjvGWRtVSvnLfKYz1A4XyufIQdDj7x551dweA1c/NMxJy/Pkai1bb1cdde7vLKjnu9/dDE/vOw4x4OQPMlNT6LFoevpQFMXRVkpJCVISIqio2+Qmrbe6MQnQM+lcCU6T5G13SURtChGJ93p38N2q3XHcWVHiUUB+vg5dfcFY8pZFKOZT5EiJopCRPJF5CUR2WPd+4y4isiwiGy0bk9GXDD3MDxyk54u9ok/jxlU701f6Um87j4BtfoXen0fHGl3UGzXWacntwWzKPwMn3ljTyPDbjVeUSh36EN0jmzVlb+LPgrHXzPmJftk8+K2I1x259scbunhj589kc+ePmvCfu78jGTHldn7m7pZWJLFotJsNh72fbx9YbfuWBQti8KVoOMNoVgUgYrtwkCx18XG1tp2MlMSqcyP7/YUIRHOSXcjfZ6mTowCiGjmU6wsituAV5RS84BXrOe+6FVKLbduH4u4VK/9t+7zf8nPdefOAJTkpPKTwU8gfe3w9i99rnOkw0H7DjvjKdd3jn9GSiJZKYl+YxSv7KgnPyOZ5RUeurZwnr5v9N3KwydKweO3aOV46e3jLClb4f3y1b2kJrl49MuncfaCiU2Zs8lzGKNQSnGgsZtZhRksr8hlS3X7mMBsIHZGo3WHN6HUUjgptpskBRnJuAQarN/QttoOFpVm4XLFTyB70mSXa0sgHEV3Pc3aKgxwoRhPzMiPfNFdrBTFZYA9s/NPwOUxkmOUxt3w1s90R84TgvfRKc1JY4eqpG3O5bDmt9B5ZMzrfYPDtPUMBrco7IlofiwK0EV3vgK4Q8NuXt/VyDkLpo1OWgMosBSFn1YePmnaDUc26xkIPkzuoqwUCjOTWVWZx+NfPp35xZM/8ealJ9PRNxR03kZT1wCd/UMjiqJ7YJi9Dc5aee840klOWpKzNirhIndGCMHsamc1FJMgMcFFYaaedDfsVuyoi/AMilhgu+46w2BV9DRpayKOMsICkZGSSGFm8tHnegKKlVJ11uMjQLGf9VJFZJ2IrBGRy/3tTES+YK23rrEx9IIsQAeAP/0wXPy/jla3O6puXXAruAfhzbHbOZ5s11qlXQ85FX5X8Tc7e/3BVtp7B8e6nUB3vMyZEVpAu2q1vp99js+XE1zCa984m3988VQKwtQvKS9DB7/bggThDlgZT7MKM1hWkQvg2P20s66DhSVZ0U0DzanQLsUhB261jtqIxidsSnL0AKOq5m56BoYj31o82uSEsZaiu3nKBLJtZkS43XjEFIWIvCwiW33cLvNcT+nOev78CJVKqVXAp4DbRcTnaDGl1N1KqVVKqVVFRZOYNzD3fGdFUkBprlYAB4aLYMUNuu7AIyZgWwClThRFdjkk+m9bMDI4yItXdjaQlCCcOc/Hj7pwnt8usj6pWq0DsQFaIGSlJoXVXTFSnR3E/XTAajY4uzCTWQUZZKcmsvFwe9D9u92KXUc6WVgSRbcTWJlPKnhev11slx25jCebaVn6N2TPoDiqMp5gtA4lHNXZU6QhoCd2F9lIETFFoZQ6Xyl1nI/bE0C9iJQCWPcNfvZRY93vB14HTvC1XiwozNAZOLXtfbqbpysJXvvxyOshjUANMt6zODuFhs7+ca2EX9lRzymzC3xPAyucr11PfsaijkEpOPi2nmUcxSvv0TYegS2K/U3dJCe4KMtLw+USllXkOsp8qmnrpXtgmIWlUT4pjtRSBIlTjBTbRd6iKM5OsRRFO0kJEnor73gnnJPueqamRVHb3kv/0HBE9h8r19OTwA3W4xuAJ7xXEJE8EUmxHhcCpwPboyZhEFzWpLsj7X2QVQKnfAm2PKQzhwjR9RQgPgHabTDsVjR1j2Y+VTV1s6+xWxfZ+aJovj4JOfHZNu+DrvqAI04jge16aglmUTR2U1mQPhKHWV6Ry64jHfQM+K9WB48ZFDGxKAieImtf/UY4RgHaKm3tGeSDQ23ML84iOfEoy4wfKboLg6LobpoyGU82lQXpKAXVrc7m5IRKrH4t/wNcICJ7gPOt54jIKhH5g7XOImCdiGwCXgP+RykVN4oCtFup1poax+lf1amR23S19pH2ftKTE8gKNCRooEefoIMoitGiu1FF8fKOegDdtsMXQabdjaHqLX0fbUXhcHjRgSad8WSzvCIXt4It1YHdTzuPdCJCWALvIZFdpuNOwQLagUaghhnbsl1/sPXoczvZ5FRA3SZtIU+U4UFt6U2RGgobO0U2UgHtmCgKpVSzUuo8pdQ8y0XVYi1fp5T6nPX4HaXUUqXUMuv+nljIGojSnLTRIHNaHkxbDDUbADjS0UtJdpARqA4ynsB30d2rOxuYX5w5Wo3tja0oGh0oioNv66E7BXODrxtGbEURaHjRsFtxsLmHWUVjFQXApuq2gPvfeaSDyvz0kCf6TZqEJB3vCeZ6Gim2i7yimGY1/xt2q6Mv48lmxfV6TOy2xya+j54WfT/FLAq76O5gc2RqKY4y+zO62LOzRybdla2A2g2gFEecjEAN0DXWk9HKWq0oOvoGef9AC+cu9JcsBmQUaVM8mEWhlA5kzzwj6umAackJpCa5aAsQo6ht62Vg2D2mfXlBZgoV+WlB4xQ7o9m6wxsnKbIjxXYBvscw4flbPGotihNv1jNjXvh3PQxsIkyxqmybwsxk0pMTONRydLmejgpKc1IZGHKP+tjLVurW5C37qe/oDz7ZzqGisAumbEXxxq5GhtyK873TYj0RsQLaQRRFy36dyjnz9MDrRYi89OSAMYr9I6mxYzu/LivPZeOhNr/b9Q4Mc6C5O7qFdp7kVEB7EEXRXqMtjwgW29nYVqkILIp2cD9auBLg0l/omqbX/2di+5hinWNtRMRKkTUWRdxh9z+qs4vhpq8AwF29jnpHVdkHITkz6I8yMcFFUVbKSMrtqzsbyEtP4oQZQWZN2JlPgTj4tr6feWbg9SJEXnrgNh4HGnVq7CyvgUjLK3Kpbe8bqTb2Znd9J0oRW4uivQaGAwTcO2qikvEEuq9WcoKLWQUZ0XfFRZPyVbDyBl0EayWWhMRI+46pZVGAnlU/ODyJ+EwAjKKYBNOtWooRRVG0EJLS6Tu4jiG3clZDkTfTkcvHLrobGnbz2q6G8dXYviicB11HtJXjj6rV2k1lxzSiTF5GUkCL4kBTN1lW5aknJ8zIBfx3krWHFUV8TrY/citADWtrzR8dNVGJT4C+4pxVmMGKiQ6ymkqc9z3dfuOZrztLD/fEtiimmOsJ4NefOoE/3XRSRPZtFMUksF1Lde2WXzAhEUqXQfV6AGcxiiBuJxu76O6Dw2209Qxynr9sJ0+KrJm//qwKpaDqbT1vI0btCrRF4T9Gsb+pm1lFGeOSApZMzyHR5b+T7I66TtKTE6jIi1Hju2DtxpWyqrKjoygA/vK5k/jeRxdH7f1iRnq+Hrh1eA1seiC0bW2LIm0CA5RiTCS7DxhFMQnsors6zz5MZStJad5KIkOBYxRKhaQo7BYML++oJ9ElnDnfwRVPsBTZ1iqdeRMjtxNYjQEDuZ68UmNtUpMSWFiaFdCiWFASw8Z3OUFqKXpbdZ1LFGoobKZlpfouzjwaWX4dVJwML31nNJPJCT1NOoMxCnGjqYRRFJPALrqra/PINChbQcJwPwukOnBVdlcDDPWGZFG09w7y3JYjnDw7n2wnf/jcSkhI9t9F1o5PVMYmkA2Ql5FMW++gz26wfYPD1LT1+lQUoOMUm6vbx1WsK6XYeSSGGU8wOsvDn0UxUkMRnRjFMYfLBZf8TCvkV3/ofLue5ikXyI4GRlFMkuk5aWMtCiugvdy1b5xffQwOM55sbDfWoZYezguUFutJQiLkz/Hveqp6W/8pihY6218EyEtP0l4YH40BDzb3oNT4QLbN8oo8uvqH2Nc4tpNsfUc/bT2DsYtPACSl6rRXv4rCHoEa+T5PxywlS+HkL8G6e6FmvbNtupumZCA70hhFMUlKrFqKEfJm0p2Qw0kpVSQGmvoWoqLwtE7GdYsNROE8/66nqtU6PuGK3c8gUNGdZzNAX9iFdx94uZ92HLFbd8Q4DTSnwr/raURRGIsiopz9ba2wn/5XPZgsGD3NUzKQHWmMopgkpbm639NI0Z0Ie5PmsUz2Bd7QVhQB2ot7Yo9UnTstk8oC31fYPimcD60HdGuCMe9/UOf5V0a3bYc3eRn+23jYNRQzC30HpGcXZpCVmjguTrGzThdbLYh2jydvAhXdtdeAJOg+YYbIkZoNH/4R1G2EdX8Mvv4U7PMUDYyimCSl2akMDLtp9kjx3Oyey4zhQzAQoPiltQqypmsXhQNKctJIcIn/3k7+KJwP7qHxY1FH6idirCjS7caA411PBxq7KcpK8RuAdbmEZeW5bPJWFEc6KMtNIyctxoHb3Ao9mMhXimZHrVYSroToy3WscdyVMOsseOWHOjboD6VMjMIPRlFMktJcXXTnOYFuzUAlLty6QZk/2g46djsBZKYk8uAXTuEr54bYj6nIT+ZT1duj/aliyGircV+up+4xrTt8sbwil51HOukdGHUrxGQGhS9yZ+h56F3141/rqI5qauwxjYgObA/2wNP/4t8F1dema1+M62kcRlFMEruozu4i2zMwxHt9M/WLgQJoIaTG2qyamR96Ve3IWFQvRXFwtc52imF8AkZdT76GFx1o6mZ2UXBFMexWbK3VRYUDQ272NnTFrnWHJ4FSZKM02c5gUTgPLvwh7Hwanvpn31Zet92+wygKb4yimCSlVhsPu7PrkfY+msihJ610pJPsOAb79IkiREUxIVIy9ZWrZxfZ9mqtqGKYFmuTkZxAcoJr3PCi9p5BmrsH/GY82YyMRrX6Pu1r7GLIrWIfyAb/RXdK6RhFjsl4iiqn3AJnfQs++KtuHOjdjnykIaBxPXljqkomSUFGsp5012YpCkth9BQtJ92fRdF+GFDRURQwPvOpyo5PxF5RiAi56UnjLIoDzb6bAXpTlJVCWW4aG62W4ztHMp7iwKIYmXTnpSh6W3UNjbEoos/Z34b+Llhzp76IOvc/R1+bwn2eIo2xKCaJyyWU5KRyxGrjYXd4lbIVOg5hm7OehJgaO2kKF+haCvsK6uBq3Qun+LjovH8Q8jPGV2fbqbHBLAqA5TNGO8nurOvUze8cbBdxkjN0YNRbUURxYJHBCxGdBbXienjzf2H17aOvTdHOsdHAKIowUJqdpmdnM9ogMH2W1Zyr1of7KeqKYh4MdOr2y6DrJ2acFjcZN7npSeP6PR1o7MYlehZwMJaX51LT1ktjZz87jnQyrzgzcA1LNPFVSzFSQ2EURUwQgUtvh+Ougpe/B+//Xi+forMookGc/JumNnYtBUB9ex9ZKYmkVa4ExHecorUKEtMgM4TCuckw0vNpl76abdkf87RYT/IzkscV3O1v6qYiP93RbOflVifZTYfb2FnXER/xCZvcCh8WRfRmZRv84EqAK+6CBRfDs9+AjQ9o6z8pA5LSYi1d3GEURRjQrqc+3G7FkY4+inNSISVLd2/1FadorYK8yuh1bPXsIhtH8QmbXB8zKfw1A/TFcdNzSHAJr+5qoKGzP7atO7zJrdQjUT0Dp3axXRQm2xkCkJAEV92rayye+DLsfs4Esv1gFEUYmJ6TxsCwm5aeAY509I/OoShbqRWFd3bFBFJjJ0VmMaRk64D2wdX6ccnx0Xv/IOSlJ9HaMzhS3a6UCklRpCUnsLAkiyc3at9/XFkUORU6cG0HSsEqtiuNG9ffMU1SKnzyb1B+ora0TXzCJ0ZRhIGRuRRtfdR7zsqefoL2e3r6qENsLx4WREYzn6rehhmnxtVJKi89mWG3oqNPT4Nr6OynZ2A4aLGdJ8sqcunq19vHRQ2FjZ0i6zkWtaPaZDzFEymZ8Kl/QNkqPU/GMI6YKAoRuVpEtomIW0RWBVjvIhHZJSJ7ReS2aMoYCtOtWorq1h4au/pHG/iVrdT3nu6nnhYY6IquogAdp6jZAM174srtBB7V2VaK7P5GZ6mxntgNAgszUyjMTAmvgJNhJEXW42Kho9bEJ+KNtFy4+SUd5DaMI1YWxVbg48Cb/lYQkQTgTuAjwGLgWhGJy/FctkWxpaadYbfSMQrQ6acJyWMD2tHOeLIpnA/9usYgngLZoIPZMNrGY7+dGhukKtuTEyxFEVfxCRht+mgHtO1iO5PxFH+4XDGb9BjvxKTgTim1A4KO7jsJ2KuU2m+t+3fgMmB7xAUMkYKMZJITXHxg5fKPWBSJybon/hhFcUDfx0JRACRnQUl8mde5VmNAW1EcaOwmJdFFabBRsh7MKcqkJDuVlfE2EzotF1JyRt2PptjOMAWJ58rsMsAzAb0aONnXiiLyBeALADNmzIi8ZF64XEJxTgqbrergMZPtpq+AjX/TjchcCaMWRW5ldIW0FcWMk+NuzOOIRWF1kLUD2aGMMXW5hJf+9UOkJsVP7GUEz3bjpobCMAWJmOtJRF4Wka0+bpeF+72UUncrpVYppVYVFRWFe/eOKM1Jo9vqYFqc4+EjL1sJg92jLTRaq3QWUnLwQrKwkj9LZ9os+mh039cBuV4dZJ00A/RFVmoSSfFSaOdJbsVojMJUZRumIBG7tFRKnT/JXdQAnlN9yq1lcYmdEpvoEgozPBWFHo1KzXqYtkgrimhbE6Bzxv91R/Tf1wHZqYkkuITWngEGh90caunhI0uPooE+ORVw4C0rPlFtLTOKwjB1iMPLrxHWAvNEZJaIJAOfBJ6MsUx+sbvIFmenjnWZFMzTcQE7TtEa2hyKsCISl8E6ESEvPYmW7kGqW3sZcquQMp7intwZuoVKX5u2KEyxnWGKEav02CtEpBo4FXhGRF6wlk8XkWcBlFJDwK3AC8AO4B9KqW2xkNcJtkVRnO2VmulywfTl2qIYGtA59LFSFHFMnlWdHUozwCmDZxfZjhpTbGeYcsQq6+kx4DEfy2uBiz2ePws8G0XRJoytKOxU2TGUrYR374SWfaDcRlH4IC9dd5C1ayhCKbaLe0bmUhzWisJkPBmmGPHseppSeLqexlG2EtyDsONp/dwoinHomRSDHGjqJjc9aWTy3VFBjscAo/YaE58wTDmMoggTZXlpuAQq8nxkM9kB7a2P6HujKMZhz6QIpcfTlCE9X3clbT9sjUA1isIwtYivhPopTH5GMg9+8VSWTPfRkC67TAcvG3foSu2s0ugLGOfkerieTpt7lDVmE9FxirpNVrGdURSGqYWxKMLIiTPzSU/2oXtFdOEd6NRYlzns3uRnJDE4rNu0H1XxCZvcGaM9v0yMwjDFMGesaGE3CMyLQQ3FFMAuuoPQmgFOGXIqYKjPelweW1kMhhAxiiJalJ2g7018wif5YxTFUWpR2BiLwjDFMIoiWkxfAUnpuqOsYRx5GUkjj2cWRrm9STSwaylMsZ1hCmKC2dEiPR/+eaOZoOUHeyZFaU6q7zjPVMdOkTXFdoYpyFH4j4xjssyVpD9sRXFUup1g1PVkaigMUxDjejLEBdlpSbjkKFYUGUWQkGLiE4YpibEoDHFBgkv436uWsXxGbqxFiQwuF5z5dT1H3WCYYhhFYYgbrlx5lKeNnv2tWEtgMEwI43oyGAwGQ0CMojAYDAZDQIyiMBgMBkNAjKIwGAwGQ0CMojAYDAZDQIyiMBgMBkNAjKIwGAwGQ0CMojAYDAZDQEQpFWsZwoqINAIHJ7GLQqApTOJEAiPf5DDyTQ4j3+SIZ/kqlVJFvl446hTFZBGRdUqpVbGWwx9Gvslh5JscRr7JEe/y+cO4ngwGg8EQEKMoDAaDwRAQoyjGc3esBQiCkW9yGPkmh5FvcsS7fD4xMQqDwWAwBMRYFAaDwWAIiFEUBoPBYAjIUa8oROSPItIgIls9li0TkXdFZIuIPCUi2dbyJBH5k7V8h4h822Obi0Rkl4jsFZHb4lC+Kmv5RhFZFyP5kkXkXmv5JhE522ObldbyvSLySxGROJPvdev73WjdpoVJvgoReU1EtovINhH5qrU8X0ReEpE91n2etVys47NXRDaLyAqPfd1grb9HRG6IQ/mGPY7fkzGSb6H13feLyDe89hX2/3CY5YvIfzgsKKWO6hvwIWAFsNVj2VrgLOvxTcAPrcefAv5uPU4HqoCZQAKwD5gNJAObgMXxIp/1vAoojPHx+yfgXuvxNGA94LKevw+cAgjwHPCROJPvdWBVBI5fKbDCepwF7AYWAz8FbrOW3wb8xHp8sXV8xDpe71nL84H91n2e9TgvXuSzXuuKg+M3DTgR+BHwDY/9ROQ/HC75rNeqiMB/OBy3o96iUEq9CbR4LZ4PvGk9fgm40l4dyBCRRCANGAA6gJOAvUqp/UqpAeDvwGVxJF/ECFG+xcCr1nYNQBuwSkRKgWyl1Bql/xF/Bi6PF/nCIUcA+eqUUhusx53ADqAM/fv5k7Xanxg9HpcBf1aaNUCudfw+DLyklGpRSrVan+uiOJIvIoQqn1KqQSm1Fhj02lVE/sNhlC+uOeoVhR+2MfojuRqosB4/DHQDdcAh4P+UUi3oL/6wx/bV1rJ4kQ+0EnlRRNaLyBciKFsg+TYBHxORRBGZBay0XitDHzObWB0/f/LZ3GuZ/d8Jl2vMExGZCZwAvAcUK6XqrJeOAMXWY3+/tYj/BicpH0CqiKwTkTUicnk4ZQtBPn/Ey/ELRDT/wyFxrCqKm4Avi8h6tLk4YC0/CRgGpgOzgK+LyOwpIt8ZSqkVwEeAfxKRD8VAvj+i/4DrgNuBdyx5o81E5LtOKbUUONO6fSacAolIJvAI8DWl1Bgr0LKyYpqnHib5KpVuT/Ep4HYRmRNn8kWMMMkXzf9wSByTikIptVMpdaFSaiXwANp3CfoH/rxSatByTbyNdk3UMPbKs9xaFi/yoZSqse4bgMfQSiWq8imlhpRS/6KUWq6UugzIRftsa9DHzCYmxy+AfJ7HrxP4G2E8fiKShD6J3K+UetRaXG+7bKz7Bmu5v99axH6DYZLP8xjuR8d8ToiBfP6Il+Pnl2j+h0PlmFQUYmW0iIgL+E/gLuulQ8C51msZ6GDdTnRwdJ6IzBKRZOCTQFiyOsIhn4hkiEiWx/ILga3e+420fCKSbr0/InIBMKSU2m6Z4B0icorl0rkeeCJe5LNcUYXW8iTgUsJ0/KzPew+wQyn1c4+XngTszKUbGD0eTwLXW9lFpwDt1vF7AbhQRPKsDJoLrWVxIZ8lV4q1z0LgdGB7DOTzR0T+w+GSL9r/4ZAJZ2Q8Hm/oK8o6dPCoGrgZ+Cr6SnI38D+MVqhnAg+hfdzbgW967Odia/19wH/Ek3zoTI5N1m1bDOWbCexCB/ReRrsi7P2sQv/w9wG/treJB/mADHQG1Gbr+N0BJIRJvjPQbofNwEbrdjFQALwC7LFkybfWF+BO6zhtwSMTC+1S22vdbown+YDTrOebrPubYyRfifU76EAnK1SjEykgAv/hcMlHBP/D4biZFh4Gg8FgCMgx6XoyGAwGg3OMojAYDAZDQIyiMBgMBkNAjKIwGAwGQ0CMojAYDAZDQIyiMBgMBkNAjKIwGCaIiMwUj/bmBsPRilEUBkOMsLoAH3XvZTj6MIrCcMwiIo9bnTq32d06RaRLRH4kerDRGhEptpYXi8hj1vJNInKatZsEEfm9tY8XRSTNWn+5tf1mazt7cM3rInK76ME0X/UhU5aIHLBaiSAi2fZzEZkjIs9bMr8lIgutdT4qIu+JyAci8rKHzN8Xkb+IyNvAXyJ8OA1HMUZRGI5lblK6ceAq4J9FpADdzmONUmoZeqbF5611fwm8YS1fgW6zADAPuFMptQTdksGeffFn4FtKqePRLS2+5/G+yUqpVUqpn3kLpHRTwteBS6xFnwQeVUoNAncDX7Fk/gbwG2ud1cApSqkT0HMW/s1jl4uB85VS14Z0ZAwGD4w5ajiW+WcRucJ6XIE+6Q8AT1vL1gMXWI/PRTczRCk1DLRbVsIBpdRGj/VnikgOkKuUesNa/id0jy6bB4PI9Qf0yf5x4Ebg86LbWJ8GPCSjozJSrPty4EGrS2kycMBjX08qpXqDvJ/BEBCjKAzHJKLnZZ8PnKqU6hGR14FUYFCNNkAbJvh/pN/j8TB68mAwugO9qJR62wqUn41uTrhV9NzvNqXUch+b/Ar4uVLqSWub7zt9L4PBCcb1ZDhWyQFaLSWxEN2yPRCvALcAiEiCZTX4RCnVDrSKyJnWos8Ab/hb3w9/Rs/FuNfaZwdwQESutmQQEVnm8Vns2Qo3eO/IYJgsRlEYjlWeBxJFZAe6FfmaIOt/FThHRLagXUyLg6x/A/C/IrIZWA78IET57gfy0G3Uba4DbhYRuxW1Pe71+2iX1HqgKcT3MRiCYtqMGwxxiIhcBVymlArrSFaDYSKYGIXBEGeIyK/Qc5MvjrUsBgMYRWEwxAwR+Q/gaq/FDymlvhILeQwGfxjXk8FgMBgCYoLZBoPBYAiIURQGg8FgCIhRFAaDwWAIiFEUBoPBYAjI/wfhuxTGSUMz7QAAAABJRU5ErkJggg==",
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEXCAYAAACzhgONAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABz8UlEQVR4nO2dd3hc1Zm432/Ue7ckW7LcK8bGNh0CoYUACSFAIAUIkAZLNtndZEN2N+WXbPqmh4SQECCBEELvvRswuODei2yr99415/fHuVcajabckabJPu/zzDMzt82ZO3fud74uSikMBoPBYPCHK9YDMBgMBkN8YwSFwWAwGAJiBIXBYDAYAmIEhcFgMBgCYgSFwWAwGAJiBIXBYDAYAmIERQQQkc+KyJpYjyMQIlIpIuc53FaJyLwJfs6E953AZ31XRO6NxmdNllDO/1RmktfOp0XkhTCOJe7/l/GKERRxTDRvsobAiMj3RWSriAyJyHdjPR6IjmAUkVNE5EURaRGRRhF5UERKI/A5s6zrPdFeppS6Tyl1gcc2U+7/ICL3ikitiHSIyB4R+VysxzQRjKA4SvH8wxnCwj7gP4GnYz2QcOHwGskD7gBmARVAJ3BXBId1tPEjYJZSKhv4KPC/IrIqxmMKGSMoJoGIlIvII9ZMq1lEfudjm3EzJRF5zZ5ZiMg8EXldRNpFpElEHrCWv2FtvllEukTkKmv5JSKySUTaRORtETne47iVIvINEdkCdDsVFiJykoi8Yx2zVkR+JyLJXptdJCIHrDH+TERcHvvfICI7RaRVRJ4XkQo/n3ORiOwQkU4RqRaRrzkZ30SxZr911rl9Q0SWeqwrEJEnrZneOhH530BmCaXUPUqpZ9E3ylDG8Hnr3HRa332lj23uFpH/9Xh/tohUebz/hnW+OkVkt4icKyIXAv8FXGVdH5utbXNE5E7rd6y2vleCte6zIvKWiPxSRJqB7wYbv1LqWaXUg0qpDqVUD/A74PRQzoHH97hYRN63zvkRL83Mvt7brO9zqqepyNf/QXyYkjy1Dus3fsL6vPeAuV7bLpJRbWm3iHxiIt8rEEqp7Uqpfvut9ZgbYJe4xAiKCWL9+Z4CDqFnWzOAf0zgUN8HXkDP3MqA3wIopT5grV+ulMpUSj0gIicAfwG+CBQAfwSeEJEUj+N9ErgYyFVKDTkcwzDwb0AhcCpwLnCz1zaXAauBlcClwA0AInIp+ob1caAIeBO438/n3Al8USmVBRwHvOJrIxE5wxJa/h5nOPxezwLzgWnARuA+j3W3Ad1ACXCd9fAcw1MicqvDz/GJiFyJvhlfC9gzyuYQj7EQuAU40TpvHwIqlVLPAT8EHrCuj+XWLncDQ8A84ATgAsDT3HEycAAoBn4gIp+yJhZO+QCwPZTv4EE3+lzkoq/Rm0TkYx7HBX3dZiql3vHc0df/wcHn3Qb0AaXo6/UGe4WIZAAvAn9HXx9XA78XkSW+DiQivw9wPQY8f9a+PcAuoBZ4xsHY4wullHlM4IG+oTYCiT7WfRZYY72ehZ5FJHqsfw34nPX6r2jVvszHcRQwz+P9H4Dve22zGzjLel0J3OBw/JXAeX7WfRV41GscF3q8vxl42Xr9LHCjxzoX0ANUeH8H4DBayGVH6Df5LnCvn3W51lhygARgEFjosf5/7d8syGfcC3zX4XieB74S7Pyjb+7/67HubKDKej0PaADOA5ICfV/0zb8fSPNY9kngVY/r8vAkzu/xQAtwZgj7jLmGvdb9Cvil8v8/Gfkf+TqW93rPbTx+40Ue637I6P/yKuBNr33/CHwnQtdmAnAG8D/ev+NUeBiNYuKUA4eU81m7P/4TEOA9EdkuIjcE2LYC+A/P2Yw1juke2xwJdQAissCaQdeJSAf6D1XotZnncQ95fGYF8GuP8bRY32eGj4+6HLgIOCTa3HZqqGN1iogkiMiPRWS/9Z0qrVWFaM0nkbHfKeTz5oByYP9kDqCU2ocW3N8FGkTkHyIy3c/mFUASUOvxe/wRPWO28fs9RWSmZdbpEpEur3Xz0JOCryil3pzIdxGRk0XkVdGm2nbgS4y/zsKFr9/4kMfrCuBkr//Sp9EaZthRSg0rpdagrQY3ReIzIokRFBPnCDBTgvsBuq3ndI9lIxejUqpOKfV5pdR09Gz79+I/suMI8AOlVK7HI10p5WnqmUg54D+g1eL5Sjvd/gt9s/ek3OP1TKDGY0xf9BpTmlLqbe8PUUqtU0pdir5xPQb809dgRORMzxuWj8eZDr7Tp9AmsvPQWsQs+/BoTXAI/af19f3CxRGc2aO78XN9ACil/q6UOgN9c1PAT+xVPj6vHyj0+C2ylVJLPQ/nbxBKqcNKm3UylVKZ9nLRPqeX0Nrs3xx8H3/8HXgCKFdK5QC3M3qdTeS6HXPeRMTzvNm/sfd1a3MEeN3rus1USvm8iYvI7QGux1BMcYkYH8UxxXtoe+OPRSRDRFJFZJyTTynVCFQDn7FmuTfgcaGIyJUiYt+wWtF/GLf1vh6Y43G4PwFfsmZmYn3uxSKSNcnvkgV0AF0isgjfM56vi0ieiJQDXwFsG/HtwDfFchSLdqZe6b2ziCSLjovPUUoNWp/n9t4OQCn1pucNy8fDyYw2C33TbEbfTH7ocfxh4BHguyKSbn3nawMdTESSRCQV/Z9JtH7vhCBj+DPwNRFZZf1e88S3o38TOlgg37rZfdXjcxeKyDmi/VB9QC9jr49ZYgUWKKVq0f6un4tItoi4RGSuiJwVZJyBvvcMtC/pd0qp232s/6yIVDo8XBbQopTqE5GT0MLcphH9veb43FPj/X/YDCwVkRXWb/Nde4WP33gJY/1QTwELROQa67dNEpETRWSxrw9WSn0pwPW41Nc+IjJNRK4WkUzrv/8htCnw5QDfMS4xgmKCWBfiR9D20MNAFdru6YvPA19H37SWAp6z7ROBdy1V/wm0an/AWvdd4B5LNf6EUmq9dazfoYXKPrSddrJ8Df2n7UQLI1+OwseBDeib2tNoxzRKqUfRM9x/WCaebcCH/XzONUCltd2X0Kp+pPgr2tRQDewA1nqtvwWtadQBf0M74O3oFETkWRH5L4/t/4S+SX8S+G/r9TWBBqCUehD4AXom3YnWovJ9bPo39E2vEn2j9zz/KcCPgSZrrNOAb1rrHrSem0Vko/X6WiDZ+s6twENoZ65PLOEdaEb8OfTN+bt+zFLlwFsB9vfkZuB7ItIJfBsPjVLpiKofAG9Z1/spPvb/LmP/D3uA76G1nb2Ad9TaLUAm+rzdjUdYr1KqE+3ovxqtHdehr+MUwodCT7qq0L/F/wFfVUo9AWNMfTOt92N+C0uLGSecY4FYjhaD4ZhGRH4ClCilrgu6sWEE0ZnTX1FK7Yz1WAyRwwgKwzGJZW5KBraitbpn0JFoj8VyXAZDPGJMT0cp4hXB4vWYGfwIRz1ZaBt2N9rU83O0eS0kAjg548JkYDCEA6NRGAwGgyEgRqMwGAwGQ0COusJxhYWFatasWbEehsFgMEwpNmzY0KSUKvK17qgTFLNmzWL9+vWxHobBYDBMKUTkkL91xvRkMBgMhoAYQWEwGAyGgBhBYTAYDIaAHHU+Cl8MDg5SVVVFX19frIcyJUlNTaWsrIykpKRYD8VgMMSAY0JQVFVVkZWVxaxZsxDxLopqCIRSiubmZqqqqpg9e3ash2MwGGLAMWF66uvro6CgwAiJCSAiFBQUGG3MYDiGOSYEBWCExCQw585gOLY5ZgSFwWAwxJSmvfoxBYmZoBCRctFtEXeIbgH6FR/biIj8RkT2icgWEVkZi7EaDAbDpHn8Fnhy3G1uShBLZ/YQ8B9KqY1Wh7YNIvKiUmqHxzYfBuZbj5PRLTtPjv5Qj16GhoZITEz0+95gMISJ5n2QnBHrUUyImGkUSqlapdRG63UnsBOY4bXZpcBflWYtkCsifrt1xSuVlZUsWrSIT3/60yxevJgrrriCnp4eAF5++WVOOOEEli1bxg033EB/v26yduutt7JkyRKOP/54vva1rwX9jJ/85CcsW7aM5cuXc+uttwKwadMmTjnlFI4//nguu+wyWltbATj77LP56le/yurVq/n1r3897r3BYAgz/Z3Q0wRdDTAFK3bHxdRRRGYBJwDveq2agW6CblNlLav12v8LwBcAZs4M3Grh/z25nR01HZMbsBdLpmfznY/4bJs7wu7du7nzzjs5/fTTueGGG/j973/PLbfcwmc/+1lefvllFixYwLXXXssf/vAHrrnmGh599FF27dqFiNDW1hbw2M8++yyPP/447777Lunp6bS0tABw7bXX8tvf/pazzjqLb3/72/y///f/+NWvfgXAwMDASE2sJ598csx7g8EQZlqtMkpDvVpopGbHdjwhEnNntohkAg+je8lO6A6ulLpDKbVaKbW6qMhn8cOYU15ezumnnw7AZz7zGdasWcPu3buZPXs2CxYsAOC6667jjTfeICcnh9TUVG688UYeeeQR0tPTAx77pZde4vrrrx/ZLj8/n/b2dtra2jjrrLPGHNvmqqvGtvf2fm8wGMJI68HR110NsRvHBImpRiEiSWghcZ9S6hEfm1Sjm7fblFnLJkywmX+k8A4xDRRympiYyHvvvcfLL7/MQw89xO9+9zteeeWVsI4nIyMj4HuDwRBGWitHX3fVQeG8mA1lIsQy6kmAO4GdSqlf+NnsCeBaK/rpFKBdKVXrZ9u45vDhw7zzzjsA/P3vf+eMM85g4cKFVFZWsm/fPgD+9re/cdZZZ9HV1UV7ezsXXXQRv/zlL9m8eXPAY59//vncddddI36PlpYWcnJyyMvL48033xxzbIPBEAPGCIr6mA1josRSozgduAbYKiKbrGX/BcwEUErdjm54fxGwD+gBro/+MMPDwoULue2227jhhhtYsmQJN910E6mpqdx1111ceeWVDA0NceKJJ/KlL32JlpYWLr30Uvr6+lBK8YtfaDn6xBNPsH79er73ve+NOfaFF17Ipk2bWL16NcnJyVx00UX88Ic/5J577uFLX/oSPT09zJkzh7vuuisWX91gMLQchNyZ0HZ4Spqejrqe2atXr1beTtmdO3eyePHiGI1IRz1dcsklbNu2LWZjmCyxPocGw5TmNyuheCnsfhZOuwXO+26sRzQOEdmglFrta13MndkGg8FwVOMe1ppE/hzInAadU8/0ZARFFJg1a9aU1iYMBsMk6KgG9yDkzYLM4inpozCCwmAwGCKJ7cjOn20JiqnnozCCwmAwGCKJLSjyZmnTk9EoDAaDwTCG1kqQBMgu0xpFT5P2W0whjKAwGAyGSNJyEHLLISERsopBuaG7MdajCgkjKAwGgyGStFZCntVGOLNYP08x85MRFAaGhoYCvjcYDJOgtVL7J8BDUEwth7YRFFGgsrKSxYsX8/nPf56lS5dywQUX0NvbC8D+/fu58MILWbVqFWeeeSa7du0aWX7KKaewbNky/ud//ofMzMygn7Nu3TpOO+00li9fzkknnURnZyd9fX1cf/31LFu2jBNOOIFXX30VgLvvvpuPfvSjnHPOOZx77rnj3hsMhjDQ1w69LR6CYpp+nmIaRVyUGY8qz94KdVvDe8ySZfDhHwfcZO/evdx///386U9/4hOf+AQPP/wwn/nMZ/jCF77A7bffzvz583n33Xe5+eabeeWVV/jKV77CV77yFT75yU9y++23Bx3CwMAAV111FQ888AAnnngiHR0dpKWl8etf/xoRYevWrezatYsLLriAPXv2ALBx40a2bNlCfn4+d99995j3BoMhDHiGxsKoRtFZF5PhTJRjT1DEiNmzZ7NixQoAVq1aRWVlJV1dXbz99ttceeWVI9vZjYveeecdHnvsMQA+9alPBW1etHv3bkpLSznxxBMByM7W9e7XrFnDl7/8ZQAWLVpERUXFiKA4//zzxwgF7/cGg2GSeIbGAiSlQUrOlDM9HXuCIsjMP1KkpKSMvE5ISKC3txe3201ubi6bNm2KyZhMqXGDIcJ4CwqYkrkUxkcRQ7Kzs5k9ezYPPvggAEqpkZLip5xyCg8//DAA//jHP4Iea+HChdTW1rJu3ToAOjs7GRoa4swzz+S+++4DYM+ePRw+fJiFCxdG4usYDAZvWg5CWh6k5owum4LZ2UZQxJj77ruPO++8k+XLl7N06VIef/xxAH71q1/xi1/8guOPP559+/aRkzN6odkmLE+Sk5N54IEH+PKXv8zy5cs5//zz6evr4+abb8btdrNs2TKuuuoq7r777jHajcFgiCCeobE2WcW6edEUwpQZj1N6enpIS0tDRPjHP/7B/fffPyJEYsFUPIcGQ8z59QqYfgJc6dEL5rlvwsa/wn9Nqlln2AlUZvzY81FMETZs2MAtt9yCUorc3Fz+8pe/xHpIBoMhFIaHoP0ILL1s7PLMaTDQBf1dkBI87D0eMIIiTjnzzDODtkA1GAxxTEcVuIdGQ2Nt7BDZ7oYpIyiOGR/F0WZiiybm3BkME8BXxBOMJt1NoQZGx4SgSE1Npbm52dzwJoBSiubmZlJTU2M9FINhauFXUJTo5ykUIntMmJ7KysqoqqqisXFqVWyMF1JTUykrK4v1MAyGqUXLQXAlQfaMscunYL2nY0JQJCUlMXv27OAbGgwGQ7horYTcmeBKGLs8PV/3p5hCGkVMTU8i8hcRaRARnw2lReRsEWkXkU3W49vRHqPBYDBMCM+qsZ64EiCjaErlUsTaR3E3cGGQbd5USq2wHt+LwpgMBoNh8rQe9C0owEq6mzqmp5gKCqXUG0BLLMdgMBgMYae3VZcY9w6NtcksNqanMHOqiGwWkWdFZKmvDUTkCyKyXkTWG4e1wWCIOf4inmwypxmNIoxsBCqUUsuB3wKP+dpIKXWHUmq1Ump1UVFRNMdnMBgM4wkqKCzTk3s4WiOaFHEtKJRSHUqpLuv1M0CSiBTGeFgGg8EQmJaD+tmvoCgBNQw9U8PyHteCQkRKRESs1yehx9sc21EZDAZDEForIb0QUrJ8r59iLVFjmkchIvcDZwOFIlIFfAdIAlBK3Q5cAdwkIkNAL3C1MunVBoMh3vEXGmszknRXDxwXhQFNjpgKCqXUJ4Os/x3wuygNx2AwGMJD60EoO8n/+immUcS16clgMBimHMOD0F7lPzQWvDSK+McICoPBYAgn7UdAuQObnlIyITlzyoTIGkFhMBgM4SRYaKxN5jSjURgMBsMxSbDQWJvM4inTk8IICoPBYAgnrZWQkAxZ0wNvN4XKeBhBYTAYDOGktRJyK8AV5PaaOXUKAxpBYTAYDOEkUNVYTzKnQX87DPZGfEiTxQgKg8FgCBdKQeuhwKGxNlMoRNYICoPBYAgXva3Q3+FMo8iye2fHv/nJCAqDwWAIF60OI55gSmVnG0FhMBgM4cJpaCwY05PBYDAckzhNtgNdXRaZErkURlAYDAZDuGithIxpkJwRfNuERMgoMhqFwWAwHFMEKy/uzRTJpTCCwmAwGMJFyIJiatR7MoLCYDAYwsHQQPDy4t5MkTIeRlAYDAZDOGg/AqgJaBQN4HZHalRhwQgKg8FgCAehhMbaZJWAexD62iIxorBhBIXBYDCEg5Fku1BMT1Mj6c4ICoPBYAgHrZWQmDqaSOcEe9vOuogMKVwYQWEwGAzhwGl5cU9GsrPjO0Q2poJCRP4iIg0iss3PehGR34jIPhHZIiIroz1Gg+GYoXEPdNTGehRTk4EeqHwTSpeHtt8UKeMRa43ibuDCAOs/DMy3Hl8A/hCFMRkMxyb3XQEv/HesRzE12fYw9LXDqs+Gtl9KFiSmxb2gSIzlhyul3hCRWQE2uRT4q1JKAWtFJFdESpVSZtpjMIST3lZoO+Ss9IRhPOv+DEWLoeK00PYTmRJJd7HWKIIxAzji8b7KWjYGEfmCiKwXkfWNjY1RG5zBcNTQsFM/N++P+5j+uKN6A9RughNv1Df+UJkCSXfxLigcoZS6Qym1Wim1uqioKNbDMRimHvXb9fNwv5U4ZnDMur9AUgYcf9XE9s+K/3pP8S4oqoFyj/dl1jKDwRBOGnaMvm7eF7txTDV6W2HbQ3D8lZCaPbFjGI1i0jwBXGtFP50CtBv/hMEQAep3QP4c/doICuds+jsM9cHqGyd+jMxiLXCG+sM3rjATU2e2iNwPnA0UikgV8B0gCUApdTvwDHARsA/oAa6PzUgNhqMYpbSPYtkV0NVoBIVTlIL1f4Gyk6D0+IkfZyQ7uwFyywNvGyNiHfX0ySDrFfAvURqOwXBs0l4F/e1QvAQK50HT3liPaGpw8HUtVC/74+SOk1min+NYUMS76clgMEQa2z8xbSkUzNORT4bgrLsT0vJhyccmd5wpUO/JCAqD4VjHjniathgK5uuop8He2I4p3umogV1PwwmfhqTUyR1rCmRnx9T0ZDAY4oCGHZBdBmm5UDAXUNByAIqXxnpkIXHrw1vYWt3OBxYUcdaCIlZV5JGUEKG58Ma/ghqG1TdM/lgZVki/ERQGgyFuqd+h/RMAhfP1c9PeKSUoGjr7+Of6I5TmpPGnNw7wh9f2k5mSyGlzCzhroRYcZXnp4fmw4SHYcA/MPXc0UmwyJCZDeoERFAaDIU4ZHoSmPTD/fP0+f65+nmKRT09vqcWt4J4bTqQ4O5W39zfz+p5GXt/dyAs79A14blEG/3b+Ai45fvrkPmzPs9BZAxf/XxhGbpEZ30l3RlAYDMcyTXt1hzVbe0jJhKzpU05QPL6phiWl2cyblgXAh5aW8KGlJSil2N/Yzet7Grn/vcN8+/HtnL+kmJTEhIl/2Lo/Q/YMmP+hMI2euK/3ZJzZhhG2VLXR2TcY62EYoslIxNOS0WUFc6eUoDjc3MOmI218dMV4TUFEmDctkxvPmM23L1lCS/cAz22bRJOg5v1w4DVYdT0khHGenVkMnUZQGOKcvsFhrvjDO9z1VmWshxKfKAVr/wBVG2I9kvBSvx1ciVC4YHRZgZVLoVTsxhUCT2zWVX0+sjywSemMeYXMzE/nvncPT/zD1v9Fn6+V1078GL6wy3jE6Tk3gsIAQE1bLwPDbg639MR6KPHJ+/fCc7fCu0dZS5SGHTokNjF5dFnhfOhrg56WmA3LKUopHt9Uw4mz8piRmxZwW5dL+ORJM3nvYAv7GjpD/7DBXn0dLLpEF/ILJ5nFuiBjX3t4jxsmjKAwAFDdpuPm6zv6YjySOKRhJzzzdf265UBsxxJuPCOebArm6efm+M/Q3lXXyd6GLj66Ylz3AZ9cubqMpATh7+9OoELutke0AD1xEnWd/BHnLVGNoDAAWqMAqGs3gmIMAz3w4PXaybvokqMra7mvA9oPj/VPgIegiH8/xROba0hwCRcdV+Jo+8LMFD60tISHNhyhb3A4tA/b+qA+N7POnMBIgzCSnT0J/0kEMYLCAEB1qxEUPnnuVmjcqev5zDx1yphkHGE3K/LOl8itAFdS3Nd8UkrxxKYazphXSEFmiuP9PnXyTDr6hnh6S4iFqNuroPi4iTUnCobRKAxTgSpLo+jsH6K7fyjGo4kTtj4EG++BM/4N5p1rZS1z9GgVDXbpDi+NIiER8mfHvUax8XAr1W29XOoj2ikQp84pYE5hBn9/L0Sndk8zZBSGto9TsuK7jIcRFAZgVKMAqDN+Ci0MnvwqlJ8MH/xvvcxORmuJT0FR09ZLRyjhzfU7IDkLcmeOX1cwL+4FxeObakhJdHHBUmdmJxsR4VMnz2TDoVZ21XU428k9rHtGpBdMYKQOSM2FhGQjKAzxTU17L9OytPpef6ybn4b64aHrwZUAl98JCUl6eV4FiCvuHNput+JPbxzgAz99lZ88u8v5jg07dCFAX6aUgnn6e7pDtONHkvZq3S8DGBp28/SWWs5bXExmSuj5DJevLCM50cXfnYbK9rQACtL9axT3rj3EP0LVUmxE4jqXwggKA8NuRW1bH6tn5QFQe6wLihe/A7Wb4WO/H9sfIDEFcsriyvTU2j3A5/66nh88sxMFHGzqdrajUjqHwjviyaZgHgwPQNskcg7Czd8ug6f/HYC39jfT3D0QNHfCH3kZyVy8rJRHN1bTM+DA1NrTpJ8zfGsUQ8Nufvb8bn77yiS0sMxp0BmfDTyNoDDQ0NnHkFuxcqYWFMe06WnXMzpX4uQvwaKLx6/Pnxs3pqd1lS1c9Js3WbO3ie9dupQLlhQ7D2/urNWO+Wl+Cv/ZxQHjRSg274em3SPa3BObashKTeTshUUTPuSnTp5JZ/8QT26uCb5xtyUo/GgU6w+10t47SHVb70gEYciULIPqjTA0MLH9I4gRFIaRC3vutEyyUxOP3VyKtiPw2E1QuhzO/57vbQrmQvOBmGbQut2K217dx9V3rCUl0cUjN5/GtafOojg7lYYOh32X663SHYE0CoifXIp9L+nnjmr6Bod5fnsdFy4tITVp4jWbVlfkMX9apjPzU0+zfvbjzH5556jJaP2h1okNaMGFMNAJh96a2P4RxAgKA1WWI7ssN42SnNTYmZ72PK9DEGPFYzdpm/wVd2kzky/y5+q2ofaNI8o0dvZz3V3v8bPnd3PRslKe/PIZHDcjB4Di7FTnUWv+Ip5sMoogJSd+HNp7X9DPva28vu0QXf1DXOowyc4fIsKnT57J5qp2tlUHyYi2TU9+nNkv7WzgjHmFpCcnsL5yguHTs8+CxFT9P4gzjKAwjGRlz8hLozg7NTYaRU8L3H81rPll9D8boOUgVL4JZ319NAzWF3b/gRg4tNceaOai37zJewdb+NHHl/Gbq1eQlZo0sr442wpGcPL71e+ArFJIz/e9XkSfh3jIpRjshco1kKGT0t7auIXCzBROnTv5CKTLVpaRmuQKHirbbU0MfAiK/Y1dHGzq5kNLi1k5M4/1lRPUKJLTYfYHdBnzOKv5FFNBISIXishuEdknIrf6WP9ZEWkUkU3W43OxGOfRTnVrL3npSaQnJ1KakxqbpLuDr4Nyj7bljDZ7ntPPiz8SeLsY5VIopbjl7++TkZzAY/9yOp88aSbiFa1UnK1bctY7MT81bNcRT4EonB8fPorKNTDUBys+BcChyj1ccnwpCa7JJ77lpCVxyfHTefz9aroCaWI9TZCaMxoB58FLVr+LcxcXs6oij111HaGFKXuy4EPQWhkfAtqDmAkKEUkAbgM+DCwBPikivvTgB5RSK6zHn6M6yGOEmrZeplsF1UqyU2ns6mdw2B3dQex/RT/X74jNbGr3s1C4MHjHslw7RDa6N9CDTd00dfXzpbPmsrg02+c2o4IiiKAfHoLGPf7NTjYF86CjCgYcRlJFir0vapPM8VcBUDjc7LOk+ET59Mkz6R4Y5vFN1f436mn268h+eWcDS0qzmZ6bxomz8nEreP9w28QGY/e4sCcucUIsNYqTgH1KqQNKqQHgH8ClMRxPTHG7Fb94YTf7Grqi/tnVbb0jlTeLc1JRStvCw82h5m7foYhKwf7XdPnm/vbo+yn62rUDceGFwbdNTNYJalGeaW+0bjwrK/L8buPY9NRyQFcqDdbq1HZoxzpvZN+Lur6SJcQXZnRwQnlu2A6/ojyXxaXZ/P3dwyh/k5TuJp9mp9buAdYfauG8JTqzesXMXBJcwoaJ+ilyy3WZkDjzU8RSUMwAPEs4VlnLvLlcRLaIyEMiUu5j/VHBuwdb+M0r+3jCSaheGFFKUd3ay4w8LShKc/SsNNwhsoPDbi7+zRrufPPg+JXN+3VxuqWX6fd2M51ose9lcA/Bgg872z4GIbIbDrWSlZrIvKJMv9tkpiSSnpwQ3PQUzJFtEw/FAZv3a0E1/wIa+4RmlcVJ+b3jzG6TwXZqb6/pYHOVH6e2n/Idr+5uwK3gvMXaf5KZksji0izWTdRPAdr8dPgdnQkeJ8S7M/tJYJZS6njgReAeXxuJyBdEZL2IrG9sbIzqAMPFo+/rWXS0s6LbewfpHhge1Sgs80W4/RRHWnro6h+ipt1HjPmBV/Xzqf+in6Ptp9jzHKTlQdmJzrYvmKud31E0kW081MrKmXm4AtjlRYSS7FTqO4P8dvU7tPmsaGHg7Wx/TFMMBYUdFjv/PJ7ZWkutKmBeqsOyGyFw6YrppCcncO/aQ7438KNRvLSznuLsFI6bnjOybHVFPu8faZ24+XbBhaCG9QQmTnAkKERktpNlIVINeGoIZdayEZRSzUope3r0Z2CVrwMppe5QSq1WSq0uKpp4Ak6s6B0Y5pmturxwtJPdRiKePHwUEH5Bsb9R27lbun0kE+1/Rdv+S1dAdll0NQr3sA69nH+B89aW+XOgv2M0CSvCtPcOsqehk1UBzE4207JTgk82GnZorSgpcKMfkjN0b+hYahR7X9BjzZ/DM1tr6UqeRmZf+MtcZKUmccWqMh7fVD3ynxhBKZ8aRf/QMG/saeKcRcVjBPiJs/LpG3Szo2aCAm3GKu0PiSPzk1ON4mEfyx6a5GevA+aLyGwRSQauBp7w3EBESj3efhTYOcnPjEte3FlPV/8QBRnJUQ9NtYsB2qan/IxkkhNcYR/HgUbte2nt8YoGGR6Eg2/C3A/qkMziJaPJYNHgyHtaxV/gwD9hE+XigJuOtKEUjgRFsSONIkDpDm8K5sUu6c4Oi51/PkopdtR04MqdAR0BnM6T4ItnzUUpuON1r9+1vwPcg+Oc2e8eaKGrf4jzl0wbs9wuhbNuon4KV4KeuOx9QQcexAEBBYWILBKRy4EcEfm4x+OzQOpkPlgpNQTcAjyPFgD/VEptF5HvichHrc3+VUS2i8hm4F+Bz07mM+OVRzdWMT0nlQuWlkRfUFizJzvqSUQozkkJe9LdAUujaPXWKKrW62zUuefo99OWQNMeLUCiwZ5ntRN93rnO94lyiOzGQ624BJY7cODqPJh+/07ZgW4dfumvdIc3dhXZWESi2WGx886nvqOfzv4hkvLKtWAfCH/L3hm5aVy+soz71x2hwVPYdvtOtntpZz2pSS5OmztWgBRnp1KenzbxfArQfoq+Nqh6b+LHCCPBNIqFwCVALvARj8dK4POT/XCl1DNKqQVKqblKqR9Yy76tlHrCev1NpdRSpdRypdQHlVIhlMacGjR29vPG3iYuPWEG03NSae0ZDL3z1iSoaeslNclFQcZoz+SS7NSwm8AONPnRKA68qu3lsz+g3xcv1bO3aMWR734OKk7XMfJOyZ0JkhA1jWLj4VYWlmQ7qpJanJ3KwJCb9l4/grZhF6BC0yj62qG7iR88vYN73q50PO5JY4fFzjp9JBowc5pVEr0jMkEfN509l6FhN3/2DLrwUb5DKcXLOxs4c36RzzIiJ1bks/5Qq3+BHYy55+gJTJyEyQYUFEqpx5VS1wOXKKWu93j8q1Lq7SiN8ajmic01DLsVHz9hxogj2XG9njBQbeVQeEaRlOSkhV2zsX0UbT0DY/88+1+B6Su1MxlGQzaj4adoOaALzS10GO1kk5CkS45HIWx02K14/3AbqypyHW1vh8j6FfROI55srOKAjYe28+c1B/nBMzs50hL+2bxP9r2oJxBJaext6ASgaIalzUXI/DSrMIOPLp/OvWsPjfrTfGgUO2s7qW7rHYl28mbVrDyauvo51DzBc5WarScwceKncOqjuExEskUkSURetrKlPxPRkR0jPPp+Fctm5DC/OIviCIWmBqK6dTSHwqYkW5ueJjwb8qKtZ4CW7gGmZaUw5FZ02hmwvW1QvUH7J2wK5uuZVDQin3Zbs7VQ/BM2+XOjYnraU99JV/+QI/8EOMjOrt8BSemQ5zAWxTKzbd28HqVAgB+H0vNiothhsfPOB2BvQxe56UnkFFfo9RHSKAD+5YPz6BkY5q63LK3Ch0bx8s56ROCcRcU+j3HiLF0aZcJ+CtDXZeMuHWEXY5wKiguUUh1oM1QlMA/4eqQGdaywp76TbdUdXHaCTh8pcZpZG0Y8k+1sbPNFm7eZaILY2oTt5Gvrto578A1dtsP2T4BOaCtcEB1BsedZKFqk236GSv4cfSOLsO1+42Ft5141009NJi+Ks4JcQw3b9Xd2Ofzr51agXEk0VG7nxFl53HT2XJ7eWju5G6AT9r6on+efB8C++i7mT8tEsq2M7I7IJWXOL87iw8eVcPdbldqE56Mg4Es761lelktRlu/ikfOKMslJS2LDRCvJgvZTwEhBxP2NXeyt75z48SaBU0FhFzi5GHhQKRWk1KLBCY9srCbBJSPlCKItKPoGh2nqGhgnKEpz9PtwaTZ2xNOqCn2za+mxVPoDr0Jy5vj8hWlLIm966muHQ29PTJsAPdMe6IKuhvCOy4sNh1opzEymPD9IKKvFNMv01ODvt6vf4dw/AeBKoD+7gvzeQ3x8ZRlf/MBcSnNS+d6TO3C7Iygk9704EhYLsK+xi3nTsnRIb3pBRDUK0FpFZ/8Qf3unUpueEtN0uDD63G6uauf8Jb61CQCXS1hVkTc5gVowV2vYe55j2K244e51fO3BzRM/3iRwKiieEJFd6DyGl0WkCDhGmxaEB7db8fimas5aUERhpv5zZ6clkprkilpRvpq2saGxNiU5lp07TOPY39hNUoJwfJl2GLfagmL/K7o0g3ehteIl0H5E38wjxb6XdDZ2qP4JmyiFyNqJdk4zkVOTEshNT/It5Lsa9OzYacSTxUE1nTmuOi5aVkpacgL/eeFCtla38+j7kfETeIbFAjR39dPSPcC8aVZWevZ03RY1ghw3I4dzFk3jzjUHGexsHGt22qUnB+f68U/YrJ6Vx/7Gbt+5Q05ZeCFUruHVLfs51NzDgabusJmEQ8GpoNgIXACsBr4B3Af8e6QGdSyw9kAzte19I2YnsEJTIxBx5A/v0FibkezsMGoUM/PTKbIEYmv3gDbbtFaONTvZ2Deyhgimzex+DtLynWdje1NgFQ+MoJ+iqaufyuYex/4Jm+KsVN8+CtucF4JGMTDkZm17HrOknpxkLawuXT6D5WU5/PT5Xc7aiIaKHRY7f9Q/ATB/RFDMiLhGAXDLOfNo7RmktrZqjNnp5Z31lOWlsbA4K+D+tp9icuanC2F4gPdffwyAzr6hsJmEQ8GpoPiWUuowcCpwHvBr4BcRG9UxwMMbq8lKSRynvobUpWyS1HhlZdtMy0pFJHwaxYGmbuYUZZKXrkNwW3sGYb9VtsPTkW1j38gi5acYHhrNxnZNsENazkztdI9g5NNG6wYTsqDISfVterLNeSFoFK/tbmDHYDGJDOl6XGizyrc/soT6jn5ufz0C33/vi9rUU3GGfmsLimJPQRFZjQJg5cw8zphXSGdLPcNpWlD0Dgzz5t4mzltcHFTLWzYjh+QE18QbGQGUn8xwcjYzG9/kpNla8ByOVtSZB04FhR3YfzFwh1LqaSA5wPaGAPQODPPctlouWlY6LgY7EjkM/qhu7cUlUJIzNncyOdFFQUZKWHwlQ8NuDjV3M7cok6zURBJcojWKA69CTvlo4TlPcsohJTtyfooj7+pkJifVYv2RkAh5syJqetp4uI2kBBnpYOeU4qwU39dQ/Q7duS7TeZmbRzZW05JqVdrxqPm0qiKfS44v5Y439k+8R7Q/9r0Is8+EJH1d7qvvJDMlccSHR/Z06G2JSNKdN7ecM4/s4XYqe/Vnv7Wvif4hN+ct9u+fsElNSmBZWc7k/BQJSWxOXc25Ce/znxfoUOVDcSwoqkXkj8BVwDMikhLCvgYvXthRR/fAMJetHF8styRHC4po2CGr2nopyU4lKWH8T1kykezsqg3jKl5WtfYyOKyYU5SByyXkpiXR3t0DB96AOWfrsh3eiOimOpEq5bHnWXAlwdwQsrF9kT9H9892gFKKf/n7Rv76TqXjw2881MrS6Tkh94Uuzk6lsbOfYW9nc/1W5/kT6LDml3fVs+Q4q8SaV82nWz+8CLeCnz4XxnBZr7BY0I7sudMyR2fwOWX6ubM2fJ/rh5Nn51OU0Mm6BhcDQ25e2llPVkriyOw+GKtn5bG1un3CSbT1HX3c27KYQmnnONHX2uHm6PcHcXqz/wS61MaHlFJtQD4mPHbCPLyxmhm5aZw0a/zFNi0rJayhqYGobu0d55+wKckOMemurwP+8iG47xMwNOq8229FPM0t0hEjeRnJZLVs030nfPknbKYt0aGckRCYu5+DWafrpKbJkD/XcYjslqp2nt5Syy9f3OPopjEw5GZzVVvIZifQSXdupZ3AIzTugdrNoxnwDnhySy2Dw4oLT1qqM9e9aj6V5aXz+TNn89imGt4/HKaS2F5hsQB7rdDYEewQ2Sj0LZGhPlJVH4f70nl4YxUv72rgAwuLSE50dutcXZHP4LBii7/y5UH46zuVvDp8PEpcpB54kWlZKRNP4psEjr6tUqpHKfWIUmqv9b5WKfVCZId2dNLQ0ceavY187ITpPktGl0Qx6a6mvXdcxNPoOPyYL/xxeK0uvVH1Hrz4rZHFdo2nOYX6j56XnsTs9ncB0RqFP4qX6qincDstm/frG57T3hOBKJgLg93QWRd00wfWH8El2j/zyMbg9vUdtR30D7knKCh8JN2t+zMkJMPK6xwf55GNVSwqyWLJ9BwdpumjiuxNZ8+jKCuF7z+1IzxasFdYbHvPIA2d/V6CwtLEo+DQtrOy0/KK+eHTO2ns7Od8B2YnG/v3m4j5qXdgmPvePcyJi+ch5SfDnueYmZ8e16YnQ5h4YnMNbgWXnVDmc320cimG3Yratr5xjmzPcbSFUneq8g19I1p1Pbx7O2zTBYcPNHWRn5FMnlVLKi89mcU9G2D6CkgPoL7bpTzC7dC2a+dMxj9hY7dNDeLQ7h0Y5slNNXxsxQyOm5HNnWsOBM1B2DBBRzb4iFrr74LN98OSjzn2Txxo7OL9w218fOUMbfIpmOezL0VmSiJfv2AhGw+38eSWSZqCRsJiLxhZtK9RJ5iNOLJhVKOIgkPbzso+4/iFdPYPkeASzl7o3MeTn5HMvGmZE4p8euT9Ktp6BrnxjNk6+a5uC8uyuzkcrxqFIXw8vLGa5WU5ozHhXjjuezxJGjr7GHIrv6ankBsYVa7RoaYf/imUnwyPfxkad7O/sZs5hRkjm5WkDLJoeBfM8RHt5Mm0xdZAwywodj8LRYu1I3qyFDjLpXhmay2d/UN84sRyPnfGHPY3dvP63sANtjYebmVGbtrI7xAK466hLQ/oUtknfcHxMR57vxqXwKUrrNl74TzorNFCx4vLV5WxdHo2P3l21+QKWo6ExY6anfaNhMZ6hKImpenQ5qgICq1RLF8wl6XTszltbgG56aHF8ayuyGN9ZUtICYput+Ivaw5y3Ixs7Q+xNOAz1EbqOvqiWjgUjKCIKrvqOthZ2zEmd8IbO7O2rj2yIbLefSi8CSk7u69d279nnaFLcFx5t/4zP3ANtQ2NzCkaFRTHD28lETfKV1isJ2l52sQQTod2b5tuMRkObQJ0kyVXUtBcigfWH2FWQTonz87nomWllGSn+m4J68HGQ60B+2MHojAzGZdY2dlKwXt/gtLlULba0f5ut+KR96s5Y37RqKAa6Z89/rsmuIRvXbKE6rZebvfu5RAKXmGxoP0TqUmu8ZpvTnRyKejWGoUrs4h/fOEUbv+Mz95pAVk9K5+OvqGRMF8nvL63kf2N3dx4xmyt0RUthNwKVjc8zHSaqGqNrlZhBEUUeXRjNYku4SPLp/vdJiUxgfyM5Ij7KOxkuzJ/pqdQsrMPr9U1m2ZZf/Ds6XDFnajmvXx94PdjNIpF3evoUSn0THPwhwt3KQ87Gzsc/glwFCJ7sKmb9w628IkTyxERkhNdXHtaBWv2NbGrzncHtJq2Xmrb+1g1M3dCw0pMcFGYmaJ9FIfegsadWptwmN29rrKFqtZeLveMyivQoZn+ut2dMqeAjy6fzm9f2cd7BycQDqoU7H1+TFgs6ByKuUWZ4/152TMinp0NjNZ5yiggKzWJDAel3r050apxtv6Q8/PylzUHmZaVwsXLrHuFCJz/PTJ7jvB8yjcYfO+uqPYIMYIiSrjdisc2VXP2wiIKMn0XErPRSXfRERRBTU9OxlH5pvZPeGY5zzmb2pX/wUcT3uHs9kdHFs9se4+17sW0Dji4aRUvgcbd4WtitOc5nWHrcGbtiIK5AUNk/7n+CAku4YqVoz6pT500k7SkBL9axah/wlkIpi9GMvzfu0NrZ8dd7njfRzZWk5GcwAVLSkYX2v6YAP2zf3DZcczMT+fL92+kqStEjbh+m87UX3TxmMX7Grwinmyyp0fPRyEJkJo74UPMzE+nMDPFcSOj3XWdvLm3ietOmzU2umrpx+j47Otsdc9m8fpvwd8+Bm2HJzyuUDCCIkrUd/ZR39HPWQsD14cBXeY74hpFay+56f5nSFmpSWSmJDrTKA6+qYWEVw/md0qv5aXhE1i4+ce65WjbYbK7K1njXkZrt4Ob/zSriVE4ejZv/gfsfFKXRPCTjT047A690F2AENmhYTcPbajigwuLmObha8hNT+bK1WU8vqlmbCc1iw2HWklLSmBRaeASEYEozk5hqK0adj4FJ1wTvD+2Rd/gME9vrR2p6zRCcrpu2LT/Fd1n3AdZqUnc9qmVtPUM8tV/bBqfxxGInU8BAgsvGlnU3T9EdVsv832VysieoZPuBsOc7OdNd5OeXDjUxnwhIpw4y3mBwDvXHCA1ycWnT545bl3ujPl80fVtnir/uu4O+ftTYd2dEdcujKCwcLsVv3xxD89sreVAY1doF7kD7Bvu9JzgzsmSnNSIO7NrfJQX96Y4OyW4oOhtg7oturifF/ubevjP4Zv1n/qf18HWBwF4w71stDBgwAGEoZTHQDc8djM8+kXdIOnc7/jcbHDYzbk/f52fvbA7tOMXzIGhXp/JX6/ubqSxs59PrC4ft+7602cz6HZz79rxM8L3D7dyfFmOz0RIpxRnp3JW51PaJHjijY73e2GH7t/+8ZU+ovI+8HU4shZe/YHf/ZdMz+Z7ly5lzb4mfvtKCF0Kdz0FM0+BzNGJ1GgOji+NIkohsj3NYwoCTpTVs/Kpau0N+n9q6urnsU01XL6yzKfTXEQoy8/kYdcFcPM7Wjt++t/hrx/VGlmEMILCoq6jj9++speb79vIOT9/neO+8zyX3vYWtz68hbvfOsi7B5ppn0QSnH2BOIliKc5OpalrgIEh94Q/Lxi++lB4Y2eJB8TbP+HBgcZucvOnIVf9Tf/hXv4eQxkl7FMznAmKwgVa7Z+on6J+B9zxQdj0d32Tu+5JyPIdA//a7kYOt/Rw7zuHQit0Z1eR9eHQfmDdEQozU/jgovFa5OzCDM5dVMy9aw+NiWDpHRhme03HhMJiPSnNTOAy94sMz78gpAivRzZWMSM3jZN9ZR6vvFbnYbz5c62d+eETq8v5+MoZ/PrlvbwZJLoL0I156rfBokvGLN5b71XjyZMcW1BE2PxkaxSTZHWFMz/FvWsPMTDk5oYz/PdIqbBzKXJnwjWPwUd+DdXvw+9P04ELEdAujKCwmJ6bxo7vXciTt5zBT684nk+eNJP0pASe217Hd5/cwVV3rGXF91/ghe3Bk6t8Yd9wvesq+WKkJaoPs0Q4UEoFzMq2cZSdXfkmJKT4rMJ6oKmLOUWZOuLm4p8DMDzrLMCq9xSMxBTdijPUyCelYMM98KcP6pIi1zwK5/yPdj774cH1R0hJdNHZP8Tjm0KYpY7kUowVFA0dfby6u4HLV83wqxl87szZtHQP8JhHue4tVW0MudWkBcXKnjcpknZaFl/reJ+Gzj7e2NPIZSfM8JkMCsBFP4MZq+DRm3S2tw9EhP/92HHMn5bJV/+xKbhWuusp/bzYS1A0dJGUIFTkp4/fx9YoIu3Q7mkKi0axZHo2aUkJPLu1jp21HT4nI32Dw9y79hAfXFjkW4uyqChIp6qlV1s9RGDVZ7V2MfPkiPXYDt2FfxRjF/FaVjZahE0pRUNnPztqO/jcPevZUtXOBUtLAhzFN3UdfSQlCPkOYrBLPDJry/J8/EkmSXvvIN0Dw5T5CY0dGUdOCg1WzaAEfzcOO38iaawAHHYrKpt6+KDtk1l5DSQkk1h2ErJxBy1OtbNpS6B6vbNtAfo74al/02au2WfBx//kV4uwaerq55VdDVx/+ize3NvE3945xNVWlFJQcsq0I99Lo3h4YzXDbuXT7GRz8ux8lk7P5s41B7nK+rwNVimME2ZOTlAsOfIAB93FtOSfgtP0sCc2WcmgPmqQjZCYAp/4K/zxLHjg0/C5l32WQklPTuT3n17JR3/3Fl++fyP3f/4UEv2Z0nY+BcXLxmk++xo6mVOY6Xu/rFL9HGmNoqc5LBpFUoKL0+YW8PTWWp7eqs2URVkpVOSnU1GQQUVBOq09AzR1DXDjGXMCHmtmQToDw27qO/pGJ3u55fCZR3QzrUn4U/wRU41CRC4Ukd0isk9EbvWxPkVEHrDWvysis2IwRoqzU/ngwmkUZ02gUJ5FfXsf07JS/c/UPPCbdNde5deRGArVfsqLe1OSncqwW/mPYBnxT4w3O1W19jAw7B6TQ8Hyq0gomE1OWhJtTkxPoDO02w7rWlLBqN2ib2DbHtYaxDWPBhUSoJPLhtyKK1eXc82pFeyo7WDj4TZn43Ml6P7THtnZSikeXH+Ek2blB5wZigifO3M2exu6eGOvDsPceKiNOYUZ5GdMojhz7RZymzZw7/D51Hc6N5c+vqmG5eW5AccMaOF45d1aOD52k19Tx7xpWfzo48tYV9nq3/fT1aCr+XppE6AjnvwlppKcbiXdRdBHMTykNdL0yWsUALd9eiVP3HI6v/3kCXz9Qws5e0ERLpfw1r4mfvHiHu56q5LFpdmcPi+wYKrI1/+pcTWfRCBl4gEQgYiZRiEiCcBtwPlAFbBORJ5QSnnaGW4EWpVS80TkauAn6Aq2MaE0N43a9olFWdR19DkyO4FHvSdPoVS3De44C5Z9Ai77w4TGYGMn2wU1PdlJd+19vn0rh9/R/onZ4x3Zdo0nXzed/PRk512/7FIeDTu1au2Pjhq4+2LdrvK6p3TRPwcopXhoQxXLy3NZUJzFjNw0fvTMLu5de8i5+adg7hhBsa6ylQNN3dz8QR8l1L24eNl0fvTMLv785gE+ML+QjYdbOceHTyMk1v0JlZjGg30fYIbDiU3f4DA7aju46ay5zj5j9plwwffh+f+CNb+EM333Mbt0xQzeO9jCH18/wIkV+Zzn3T5019OAGuef6Bsc5nBLz2hmuC8i3Zei1/InhMH0BNpicXxZLseX5Y5b1zc4zJGWHgozU4JqshUF2spwuKWbU+dOXttxQiw1ipOAfUqpA0qpAeAfwKVe21wK3GO9fgg4V5z2hIwAJTmpE27mU9/R71hQ5KUnkZzgGtUolIJnvq6TxTb/XZehmATVflqgemObwPxqUZVrtH9ixvi8BDtiZY4PQZGbnuS8Oq5dFjtYKY/n/xuG+uGzTzsWEgDbqjvYVdfJlat0lE9GSiKXr5zB01tqx1ZfDUT+HC0o3Dr44IF1R8hMSeSiZcFNlMmJLq47TZu8nt9eT0v3wOT8E72tsOVBOP4T9CZkUe/Qz7W7rpNht+K4GSFU1D3lZp2f8cr3Yd/Lfjf71iVLWDo9m/94cDNHvAva7XpKm5zsCYHFgcZu3MqPI9sm0rkUVkHAcJiegpGalMD84qyRmmiBKM1JJdElUa0iG0tBMQM44vG+ylrmcxul1BDQDoz71UTkCyKyXkTWNzY6iLKYINNzUqlp7w25SqZSirr2vtHGK0EQEaZlezQO2vogHH4bLvo/KD4OnvwK9Ey8GUpNWy+pSS4KglyUxVZ2tl+HduWbUH7SOP8E6K52uelJPk0o+RkhaBS5MyE5K7BDe/+rsP0RPastcDgjtnhwg3Zie2bLf+aUCgaG3fxzvcMy1gVzdY2izho6+gZ5emsNH1k+nfRkZwr7p0/WCXj/89g2YGKFAEfY9HcY6kVO+jzTspx3S9xeo017S6eH0CRJBD76W1076+Eb/YZnpiYl8PtPr8TtVvy/Jz0Efl87HHhdaxNe87+9DVYxwGkBTCk5Ec7OtgoChkujCBeJCS5m5KVFtYrsURH1pJS6Qym1Wim1uqjIeWXHUCnJSaNvMPReER19Q/QODjsWFODR6a6vA174H50DsPpG+Njv9QX83DiXjmOq23TEUzDlrDAjhUSX+A6R7W3TPgEf/gnQ1Uc9S3d4kpue7NxHYTcx8hciO9QPz3xN+wlO/6qzY1r0DQ7z+KYaPrS0hJy0pJHl84uzOGVOPve9e8hZPo0d+dS8nyc319A36OaqE/07sb3JTU/m8lUzaOrqJys1kXnBfAT+cLt1eOTMU6FkWUj5ONtr2slOTQwa4DCO5Ay4+l5tgnzgM367zlUUZPCZUyt4dXfjaNWBvS/qhMrFHx23/b6GLhJcwqzCAMEcdqe7SCXd9URPowiVmfnpUa0iG0tBUQ14/pvKrGU+txGRRCAHaI7K6HxgJ8uF6tAeyaFwaHqyt63v6Ic3fgpd9VqbcLl0qOkHvq4rgu58KqRx2FS3Bs+hAN0buTg7lXpf3/fwO4DyKyj2N3b7dYrmpSfR4lRQgE68q/fTxOjt3+rM7Yv+z6dmE4gXd9TT3jvIlavHJ5ddc8osqlp7eX1PQ/AD2bkULQf457ojLCzOYnlZCDNz4IbTddz8CTPzHAU8+GT/K9B6EE78HGAlTDoUFNtqOlgyPdtZpJc3+XPg43/WfrTnv+l3sytWlTHsVjxqhwPvfBIyi32GVu9r6KIiP52UxADd/bKt3y1SDu0R01N8aRSg/RSHotjpLpaCYh0wX0Rmi0gycDXwhNc2TwB2p5UrgFdUNHqE+qFkRFCENoMZyaEIUaNIb9+HWvsHXYKhzKOI3pn/ASXL4KmvjlS3DAUnyXY2xdl+Ir0OvgmJqT79E519gzR29vv0T4Ductc36HZeKnnaUt3j2jv7ufUQvPF/ekbqUZraKQ9uqGJ6TiqnzR1/I7hgaTFFWSn87Z1DwQ+UPQMSU2k+spPNVe0jBQBDYU5RJo8c9zb/PWcS1Vff+6O+8VozdKemp6FhN7tqO0IzO3mz4AI47RbYcDccesfnJnOLMlk5M5eHNlShBnq0RrHwIj0B8mJvoIgnm0j3pbBNT4H6psSIivwMOvqGJpUEHAoxExSWz+EWdIvVncA/lVLbReR7ImLroncCBSKyD/h3YOL2ljBgRwmFqlHYM/JQBEVxVjLf5C+QlAHnfXfsyoQk+NgftPnn2dA60vYNDtPUNeBYUJTm+Em6q3zTZ/4EeHS1K/JtesqzckkcZWeDRykPL/PTc7eCuODCHzk7jge17b28ubeRy1eV+cwRSUpw8cmTZvLansbgKr7LBXmzaazcQXKCK2AZeb8oxcpDd7Fw3bf9mm8CUrMJ9r6gtYlEfX5LclLp6h+iqz9wpvmBpm76h9wsnT7J1rBnfxNyynUei59CjlesKmdvQxcH1z2juwP6CIsdGHJT2dQd2JENkS/j0d2kiwEmJAXdNNrMtCKfDrVER6uIqY9CKfWMUmqBUmquUuoH1rJvK6WesF73KaWuVErNU0qdpJRy1sk+QhRmapv9RDUKu9eEE1Z0vcEZCdtpPOnrvp1pJcvgrG/onIHtjzk+bk2QqrHe2FVIxyhyva1Qt9VnfSfQGdkw2ifbG1tQOHZo+4p82vUM7H4Gzv6GjusPkUc2VqOUNof445MnleMS4b73gmsVzSllJLYd5IKlxRPLgehr1zfO7gZYf2fo+7/2Y31TO/mLI4uKs4MEI1hsr9H9nCelUYD2V1z0M13W/J3bfG5yyfJSUhJdtG54GFJyYNb4Ht6HmrsZcqvAjmyIfO/sMCXbRQI7RDZakU9HhTM7WiRYNvuQfRQdfeSlJ5GaFMDe6slAN8u3/5Qd7gp2l13hf7szvgqlK3RRsC5n0V5OQ2NtSnJS6BkYpqPPY1Z6KIh/oqGbBJcwM9+foNAzNMdBAen5OhPX1igGeuDZb0DRIh2iGSIjCXGz86ko8D1G0NrU+YuL+ee6IwHNZM9tq+WRQ6nMlHq+fcmikMcDjJpPkjJgza98dpLzS/UG2PMsnPZlSB292RdnOeuWuL26g5REl1/BHhILP6yjmF77sTYNepGdmsRFSwuZ0/Imw/POH9F+PLG72gU1PSWn6xLqkdIowlS+IxLMzLdzKYygiEtKclKpbQvd9GQnrznizV+Q0lPLtwY/S11nALOBbYLq79TCwoH7ZqSznUONwh73mJtN5Rrtn/DT1+FAUxcz89PH1tL3wJ5xO9YoQGsVdhXZN38O7Yfh4l9MyCyw/lArlc09I7kTgbjm1ApaewZ5ZqvvftD/XHeEm+/byFDuLJIZZJp7grEW9s3ug/+lb1Dr/uR831d/pLOUPbQJYKS0eTA/xfaaDhaVZPkvsREqH/6JNgk+83Wf1+QN5fXk0cn7Gb4nGnsbuhDxUzXWm+yyCJqemuPSkQ26REphZkrUHNpGUIRIqZOKql7UdfRR4tTs1Lwf3v4Nw8uuYoNaGDy8sXiJtg3vfELnEgShpq0XlzgrTgijfpUxiYaVb+j8iUTf3+mAV59sb+zyyY5DZEEnZDXt1hnab/0ajr86pMQ6Tx5cf4T05AQuWlYadNvT5hYwpyiDv60dPzv+85sH+M+Ht3D6vEKu/8i5emGQ/tl+sc0nSy6Feefr79jfGXy/I+/Bvhfh9K+MK99g/8aBriGlFNtr2lkyWbOTJzllWuDtfX604J8HSzveoJ8k/ljju6bR3oYuyvLSxvbD8Ef2dOiIlOmpCTLi0/QEduST0SjiktKcVGraQku6q3davkMpbVJJSCHhgu+Rk5akQ2SDcdq/6oqeT/8HdNYH3LSqrZfi7FTHvQ7GCYqeFh0G6cc/MexWHGzq9uvIBp2ZDdDipHmRTfFSGB6AB66BpHRdPmIC9AwM8fSWWi5eVuqoraWI8JmTK3j/cBvbqrUtXynF/z2/m/99eicXLyvlz9etJrV4gd4hSP9sv3TU6Fl4Vgl88JvaD/TuH4Pv9+oP9Kz3pM+PW5WZkkhGckLAiU1Vay8dfUOTd2R7c/KXdKG/Z78xVuAphWv3M1TlncLL+7t8VjrYW98Z3D9hkz09MhqFUpaPIj41CtDlxsdlukcIIyhCpDQnjf4h50l3/UM6yshJHwr2PKdnh2ffClklo0l3wUhIHI2C2nhPwE2d5lDY2A74kXEEyZ+oaeulf8gd0GyQlOAiKzXRedQTjDq0m/fCud8a0+AmFJ7ZWkf3wDBXBqjq6s3lq8pITXJx79pDuN2Kbz2+jd+9uo+rTyznN588Qcf6Z5VCYtqYmk8h0VGtQ1sTkrTQX3ChzhHpa/e/T+VbcOA1OOPftCPZB7qtrv/JxmhGdpgFRUIiXPJLfRN/1SMqrXYTtB8h54TLcCt4eONYbWDYrTjQ1O27/akvcmboG/pgmEvy97Xrkjlx6swGHflU29FH/9DkC4UGwwiKECm1NIMah5FP9p80aGise1iHexYtGrE1F4fS6a5oIeRVQOOugJvVtPc6dmSDLr+Qn5E8Kihs/8SMVT63D1TjyZP8jOTQBEXRQnAlauf96huc7+fFg+uPMKsgfaThvRNy0pL42IoZPLapmlvu38i9aw/zxbPm8KOPLxsNrXW59Plvc5B34YuO6tFwT9DmxL42WHu7/31e+5EWLgHOx5hSMD7YUdOOS2BRSZgFBUD5ibpXwrt/gNrNetnOp0BcFK76GCfNyufhDVVjtPMjLT0MDLmZ61RQjITIhjmXIk7Ld3hSUZCOUnCkJcLtYDGCImRKc0crqjrB/pMGzcruqNa1ck7+4oiDtsRJK1JPChdAk+9mMqBna7VtfY5DY22Ksz2KIdr1nQL4J8B/DoVNbnoyraEkCyWmwFX36V4IfnpeB+NQczfvHmzhilVlISfEfeaUCvoG3TyztY5vXLiIb3548fhj5JRPvNl9e/VouCfA9BU6euid27Sm6M3BN/Rvcca/6+gfP5RkpwYsDLi9poN50zKd+QMmwnnf0bPyp/5NT4Z2PQUVp0NGAVesLuNAUzcbrR4coP0TgHONIlK5FHGclW1jRxUejkIuhREUITKqUTi7gTvOyrYLqtnlILBbovYzNOywJWrhAmjaN1LF1JuGzj6G3Cok0xNYDvz2Pg//xPjYd5sDTV1kpyYGLTiYl57krMudJwsv1LP2CfLQhipE8N0POgjHzcjhlg/O4xefWM5NZ/spPJg7E9qO+F4XCKX0jc47H+TsW6G/Hdb+fvz2r/4QsqbrGXsAirN1KRh/PrXtNZPMyA5GWh586Ic6hPeFb2mN1yopftGyUtKSEnhow6j5yS4GGDQ01iZiGoUlKOLcmQ3RyaUwgiJE7KS7OoempzqnWdm2oPDo8lWcnYpbQVOXwxtq4XwY6oV23zerkdDYEAu/6ZtNHxx6m0D+CdA5FHOnZQadseenh2h6CgOPvl/NGfMKQ9aobL72oYWBhUxuuS5SF0oOBGgT02D3WI0CdFLl4o/CO78fWy34wKvaV/SB/wha32padioDfnxqTV391HX0hd8/4c2yK3W3wbVWEt6iiwGsUuylPLm5lt4BbWffV99FaU4qWakOw54jVcZjpHxH/AqKgoxkMpITjKCIR0aS7hzmUtR39JGS6BqJ9PFLa6W2wXvYqUcijpz6KQoX6uemvT5XO+1s501JdirN3QMMHXxDO2xnrPS77YGmLuYUBp8N6gqy0alTA/qmWNXay1kLIlddmBzLQe5HUPvFNpt4+ihszv6mbm/5zu/0e1ubyCnXNcCCMNJW14f5yXZkL4m0oBCxcl6StY8pdzSQ4IpVZXT1D/G81Yt+X6ODGk+eRCrpbgqYnkSE8vz0qCTdGUExAUpznGdn11kNi4LaxFsr9Z8/YTRk02enu0AUWiGafvwUExUUtrnNfWBNQP9EV/8Q9R39Qf0TAPkZSXT1DzEw5NCsNkl212mTRkSctja5llksVD+F3VPBl6AoXgJLL9Ohst3NsO8lqFoHH/ia399hzO521JqPa2ikdEdpBE1PNoXz4JP3w0d+NWbxybPzKc9P48ENR3C7VeD2p/7IjkBfip5mHYYdwP8TD1QUGEERt5TkpDqu91Tvr42oN62V45rL++2d7Y+MAp2h2+S7P3F1ay+56UmO8gfGjCMnlQ+4NpPctB3mnet3u4Mj7U+DC4oJJd1Ngp21eva8qDQyPYWB0ZlyqILCNpvk+CkmeNY3YKAb3v61zpvIrYAVn3Z06OIA2dnbazooy0sjJ5i2Gy7mnQfTTxizyOUSLl9Zxtv7m1lX2ULPwLDzHAqbSHS6626Ka23CpqIgg8MtPbid9E2ZBEZQOGWwD178NjTuZnpuGrXtfY6S7nRW9sQERUFGMokucS4oQIeR+jE91YRQXtyT6cm9/DTpDjqy5sJJX/S73WgxwOAzwpEyHlESFLvrOinMTKEw03lhxpDJmKZbw05EUIgLMv20Tp22CJZdofMqat6Hs/7TcemSoiz/hQF31HRE3j/hgMtXlqEU/PR5PcEJWjXWm0j0zo7zrGybmfnpDAy5Hbe8nShGUDjl2a/rsgpPfpXS7BT6h9xBwzuVUlpQBAuN7e/Uqq6XoHC5hGlZzprPDA67+dSf1vJiYw6dVTv44+v7eXlnPZVN3SNd2uzOdqFS8e53KKCDlxZ9P6DzdH9DFy4ZLYEcCNtn0xpKdvYk2FXXyeJIahOgcylyyibmo8gsGWN2HMdZ39DPebN1+RKHpCYlkJeeNO5G0tU/xMGm7shGPDmkPD+dU+cUsOGQDpMNucNfdgSS7uK4cqwn0Yp8Cs0Gcazy/n2w8a9Qcjwcfptl5euATGrbewOWlG7tGWRgyB3c9GRX2fQR+uk06W59ZStv72/m7NzpnD/cyu+fXU87+g+XnOBidmEGB5u6fTbpCci2h0ne+Qi/Vp+gQ83m4wE23d/UTXmwrmQW9nmLRuTT0LCbPfWdXHvqxENrHTORENn2qvERT94UzofL7oD82YEFig90HsxY05NtiosHjQLgytVlvHOgmcLMZPJCLdNum+w6a0bb0k6W7ubR4JA4psLOpWju4ZQ5kRNsRqMIRt02XZl11plw4wuQN4ulO3+F4A4a+TSZ0FibEs9ktwC8sque5AQX13zkfADeurGMh286jZ9efjzXnz6Lsrw0ZhdmcNbCEKJ+Omp1/agZq3gi6+qgmk2wYoCehNy8CPjuE9t5Y4+zcuqeVDb30D/kjqwj2yZ3Akl3HTX+/ROeHH+l34q9gZiWnUqDl0axvTpMPSjCxIXHlZCZkhi6Ixs8+lKE0fwUxyXGPZmem0qCSyLewMhoFIHoa4d/XqMbwlzxF0hKgw/+N2mPfJ5LXGup7Tg+4O62JlCSE8QuHkBQFGen8ubepqBDfXlnAyfPySetVNdAyuw8wKqVp7OqwnmpijEoBU98Wavzl/2RokeaAgqsdw80s7+hizPmOZvVjJqenAmKrv4h7n67krr2Pj4QYojrrrooOLJtcmfq5kODvfp6CYZS2r4+//yIDak4K4Xd1jmw2V7TQUFG8khUVKxJT07k11evCB5G7otwZ2cP9MBgz5QwPSUmuJiRmxZx05PRKPyhFDz+L9osdOVdo0XojrsCVbyUryU9SH1L4DLQI1nZwXpRtFbqhjNp42/qxdnB21keaOziQFM35y0u1jeqhJSApTwcseFuXaDw/O9B4XxKc9L8Cop/rj/CZ+58l7L8ND57+mxHh09JTCAjOcFxGY/KJj1j2nSkzdH2nuyq7STBJRObrYZKzkz97LTrWl+bvin5Co0NEyU5qTR29o/4qkALiiXTs0MuZRJJzl1czKqKCfSnHkm6C1O58SmQbOdJRUHkq8gaQeGPd26DnU/C+f8PKk4bXe5yIed+hwqpp+Jw4P4Pde19iMC0LAcahQ9tAka1kUB+ild2NQBwzqJpug5SwbzJCYqWA/D8f8Ocs3UPZqwqpJ19Y8Lwht2KHz2zk/98aAsnzy7g0ZtODymqKjc92bFGccASFHUdfaHVv0JrFHOLMhz5TiZNqCGyIzkUQXwUk2DaSIa/9lMMDLnZ29AZN2anSZOcobX+cGkUI+U74t/0BDry6ZARFDHg0Ds6FHbRJXDqLePXz7+AXUlLOLfhLq2m+qG+o4+CjJTgvR8CCIqRXIoAN8eXdtazsDiLcqs9IkWBiwMGxD0Mj96ks8QvvU1H8qALFA4OK5qtG3t3/xBfuncDf3zjAJ85ZSZ3XX9iyPH4oVSQtXM0IHStYmdtZ3T8E6A1OnAuKOyb2wT6fjul2CtEdk99J4PDKm4c2WEhJ4yd7rptjWJqCIqKgnTaegZp741cBKERFN50NcJD1+sIpI/9Xpcf8EaEZ0u+RL67Bd7z31xGh8YG0Sbcbl2a2p9GEaSMR3vPIOsqWzl3sUd/hsIFWvgMOWh65M3bv4Uja+Gin465eXm2RK1p6+WK29/h5Z31fPcjS/j+pcc5boTkSW56Ei0OTU8Hm7ooykohKUFCEhQdfYNUt/VGxz8Bui+FK9F5iKxtLomgRjHa6U5fDzus0h3HzThKNArQ58+puS8YU06jGI18ihQxERQiki8iL4rIXuvZp8dVRIZFZJP1eCLiA3MPw8M36O5in/jrmEb13vSVnsRr7hNQa36pt/dBXbuDZLvOWt25LZhG4af5zOt7Gxl2q/GCQrlDb6JTt01n/i7+CBx/1ZhV9s3mhe11XHrbWxxp6eEvnz2Rz54+e8J27vyMZMeZ2QeaullUksXi0mw2HfF9vn1hl+5YHC2NwpWg/Q2haBSBku3CQLHXZGNbTTuZKYlU5Md3eYqQCGenu5E6T1PHRwFENPIpVhrFrcDLSqn5wMvWe1/0KqVWWI+PRnxUr/5Q1/m/+Be6cmcASnJS+cngJ5C+dnjrNz63qetwUL7DjnjK9R3jn5GSSFZKol8fxcs768nPSGZFuYesLZyvnxt9l/LwiVLw2E1aOF7yq3GalC3wfvPKPlKTXDxy82mcvXBiXeZs8hz6KJRSHGzsZnZhBivKc9la1T7GMRuIXdEo3eFNKLkUTpLtJklBRjIugQbrGtpe08Hi0ixcrvhxZE+a7DKtCYQj6a6nWWuFASaK8cTM/Mgn3cVKUFwK2D077wE+FqNxjNK4B978ua7IeULwOjqlOWnsVBW0zf0YrP0DdNaNWd83OExbz2BwjcLuiOZHowCddOfLgTs07Oa13Y18cOG00U5rAAWWoPBTysMnTXugbovugeBD5S7KSqEwM5nVFXk8dvPpLCie/I03Lz2Zjr6hoP02mroG6OwfGhEU3QPD7GtwVsp7Z10nOWlJzsqohIvcmSE4s6uc5VBMgsQEF4WZutPdsFuxszbCPShigW266wyDVtHTpLWJOIoIC0RGSiKFmclHn+kJKFZK1Vqv64BiP9ulish6EVkrIh/zdzAR+YK13frGxtATsgDtAP7MQ3DRzxxtbldU3bbwFnAPwhtj93Pc2a61Upsecsr9buKvd/aGQ6209w6ONTuBrniZMzM0h3blGv0854M+Vye4hFe/djb//OKpFISpXlJehnZ+twVxwh20Ip5mF2awvDwXwLH5aVdtB4tKsqIbBppTrk2KQw7Mah01EfVP2JTk6AZGlc3d9AwMR760eLTJCWMuRXfzlHFk28yMcLnxiAkKEXlJRLb5eFzquZ3SlfX82REqlFKrgU8BvxIRn63FlFJ3KKVWK6VWFxVNot/AvPOcJUkBpblaABwcLoKV1+m8Aw+fgK0BlDoRFNllkOi/bMFI4yAvXt7VQFKCcOZ8Hxd14Xy/VWR9UrlGO2IDlEDISk0Kq7liJDs7iPnpoFVscE5hJrMLMshOTWTTkfagx3e7FbvrOllUEkWzE1iRTyp4XL+dbJcduYgnm2lZ+hqye1AcVRFPMJqHEo7s7ClSENATu4pspIiYoFBKnaeUOs7H43GgXkRKAaznBj/HqLaeDwCvASf42i4WFGboCJya9j5dzdOVBK/+aGR9SC1Qg7T3LM5OoaGzf1wp4Zd31nPKnALf3cAKF2jTk5+2qGNQCg69pXsZR3HmPVrGI7BGcaCpm+QEFzPy0nC5hOXluY4in6rbeukeGGZRaZRviiO5FEH8FCPJdpHXKIqzUyxB0U5SgoReyjveCWenu56pqVHUtPfSPzQckePHyvT0BHCd9fo64HHvDUQkT0RSrNeFwOnAjqiNMAguq9NdXXsfZJXAKV+CrQ/qyCFCND0F8E+ANhsMuxVN3aORT5VN3exv7NZJdr4oWqBvQk5sts37oas+YIvTSGCbnlqCaRSN3VQUpI/4YVaU57K7roOeAf/Z6uDRgyImGgXBQ2Tt2W+EfRSgtdLWnkHeP9zGguIskhOPssj4kaS7MAiK7qYpE/FkU1GQjlJQ1eqsT06oxOpq+TFwvojsBc6z3iMiq0Xkz9Y2i4H1IrIZeBX4sVIqbgQFaLNSjdU1jtO/okMjt+ts7br2ftKTE8gK1CRooEffoIMIitGku1FB8dLOegBdtsMXQbrdjaHyTf0cbUHhsHnRwSYd8WSzojwXt4KtVYHNT7vqOhEhLI73kMieof1OwRzagVqghhlbs91wqPXoMzvZ5JRD7WatIU+U4UGt6U2RHAobO0Q2Ug7tmAgKpVSzUupcpdR8y0TVYi1fr5T6nPX6baXUMqXUcuv5zliMNRClOWmjTua0PJi2BKo3AlDX0UtJdpAWqA4insB30t0ruxpYUJw5mo3tjS0oGh0IikNv6aY7BfOCbxtGbEERqHnRsFtxqLmH2UVjBQXA5qq2gMffVddBRX56yB39Jk1Ckvb3BDM9jSTbRV5QTLOK/w271dEX8WSz8lrdJnb7oxM/Rk+Lfp5iGoWddHeoOTK5FEeZ/hld7N7ZI53uZqyEmo2gFHVOWqAGqBrryWhmrRYUHX2DvHewhXMW+QsWAzKKtCoeTKNQSjuyZ50R9XDAtOQEUpNctAXwUdS09TIw7B5TvrwgM4Xy/LSgfopd0Szd4Y2TENmRZLsAv2OY8LwWj1qN4sQbdc+Y5/9LNwObCFMsK9umMDOZ9OQEDrccXaano4LSnFQGhtyjNvYZq3Rp8pYD1Hf0B+9s51BQ2AlTtqB4fXcjQ27Fed5hsZ6IWA7tIIKi5YAO5Zx1euDtIkReenJAH8WBkdDYsZVfl5flsulwm9/9egeGOdjcHd1EO09yyqE9iKBor9aaRwST7WxsrVQEFkfbuR8tXAlwyS91TtNrP57YMaZY5VgbEbFCZI1GEXfY9Y9q7WS46SsBcFetp95RVvYhSM4MelEmJrgoykoZCbl9ZVcDeelJnDAzSK8JO/IpEIfe0s+zzgy8XYTISw9cxuNgow6Nne3VEGlFeS417X0j2cbe7KnvRCliq1G0V8NwAId7R3VUIp5A19VKTnAxuyAj+qa4aFK2GlZdp5NgrcCSkBgp3zG1NArQveoHhyfhnwmAERSTYLqVSzEiKIoWQVI6fYfWM+RWznIo8mY5MvnYSXdDw25e3d0wPhvbF4XzoatOazn+qFyjzVS2TyPK5GUkBdQoDjZ1k2VlnnpywsxcwH8lWbtZUcT7ZPsjtxzUsNbW/NFRHRX/BOgZ5+zCDFZOtJHVVOLc7+jyG0//h7PwcE9sjWKKmZ4AfvepE7jnhpMicmwjKCaBbVqqbbfsggmJULocqjYAOPNRBDE72dhJd+8faaOtZ5Bz/UU7eVJk9fz1p1UoBZVv6X4bMSpXoDUK/z6KA03dzC7KGBcUsHR6Doku/5Vkd9Z2kp6cQHlejArfBSs3rpSVlR0dQQHwt8+dxHc+siRqnxcz0vN1w60ja2Hz/aHta2sUaRNooBRjIll9wAiKSWAn3dV61mGasYqU5m0kMhTYR6FUSILCLsHw0s56El3CmQsczHiChci2VurImxiZncAqDBjI9OQVGmuTmpTAotKsgBrFwpIYFr7LCZJL0duq81yikENhMy0r1Xdy5tHIik9D+cnw4rdGI5mc0NOkIxij4DeaShhBMQnspLvaNo9IgxkrSRjuZ6FUBc7K7mqAod6QNIr23kGe3VrHyXPyyXbyh8+tgIRk/1Vkbf9ERWwc2QB5Gcm09Q76rAbbNzhMdVuvT0EB2k+xpap9XMa6UopddTGMeILRXh7+NIqRHIro+CiOOVwuuPjnWiC/8n3n+/U0TzlHdjQwgmKSTM9JG6tRWA7tFa794+zqY3AY8WRjm7EOt/RwbqCwWE8SEiF/rn/TU+Vb+k9RtMjZ8SJAXnqStsL4KAx4qLkHpcY7sm1WlOfR1T/E/saxlWTrO/pp6xmMnX8CIClVh736FRR2C9TI13k6ZilZBid/CdbfBdUbnO3T3TQlHdmRxgiKSVJi5VKMkDeL7oQcTkqpJDFQ17cQBYWndjKuWmwgCuf7Nz1VrtH+CVfsLoNASXeexQB9YSfeve9lftpZZ5fuiHEYaE65f9PTiKAwGkVEOfubWmA/9e+6MVkwepqnpCM70hhBMUlKc3W9p5GkOxH2Jc1nuewPvKMtKAKUF/fEbqk6b1omFQW+Z9g+KVwArQd1aYIxn39Ix/lXRLdshzd5Gf7LeNg5FLMKfTuk5xRmkJWaOM5PsatWJ1stjHaNJ28CJd21V4Mk6DphhsiRmg0f+gHUboL1fwm+/RSs8xQNjKCYJKXZqQwMu2n2CPHc4p7HzOHDMBAg+aW1ErKmaxOFA0py0khwif/aTv4oXADuofFtUUfyJ2IsKNLtwoDjTU8HG7spykrx64B1uYTlZbls9hYUdR3MyE0jJy3Gjtvcct2YyFeIZkeNFhKuhOiP61jjuMth9lnw8ve1b9AfShkfhR+MoJgkpbk66c6zA93agQpcuHWBMn+0HXJsdgLITEnkgS+cwpfPCbEeU5GfyKfKt0brU8WQ0VLjvkxP3WNKd/hiRXkuu+o66R0YNSvEpAeFL3Jn6n7oXfXj13VURTU09phGRDu2B3vgqX/zb4Lqa9O5L8b0NA4jKCaJnVRnV5HtGRji3b5ZemUgB1oIobE2q2flh55VO9IW1UtQHFqjo51i6J+AUdOTr+ZFB5u6mVMUXFAMuxXbanRS4cCQm30NXbEr3eFJoBDZKHW2M1gUzocLvg+7noIn/9W3ltdtl+8wgsIbIygmSalVxsOu7FrX3kcTOfSklY5Ukh3HYJ++UYQoKCZESqaeuXpWkW2v0oIqhmGxNhnJCSQnuMY1L2rvGaS5e8BvxJPNSGtUq+7T/sYuhtwq9o5s8J90p5T2UeSYiKeocspNcNY34P17deFA73LkIwUBjenJG5NVMkkKMpJ1p7s2S1BYAqOnaAXp/jSK9iOAio6ggPGRT5W2fyL2gkJEyE1PGqdRHGz2XQzQm6KsFGbkprHJKjm+ayTiKQ40ipFOd16CordV59AYjSL6nP1N6O+CtbfpSdQ5/zO6bgrXeYo0RqOYJC6XUJKTSp1VxsOu8CozVmo/hK3OehJiaOykKVyocynsGdShNboWTvFx0fn8IORnjM/OtkNjg2kUACtmjlaS3VXbqYvfOdgv4iRnaMeot6CIYsMigxciOgpq5bXwxs9gza9G103RyrHRwAiKMFCanaZ7ZzNaIDB9tlWcq8aH+SnqgmI+DHTq8sug8ydmnhY3ETe56Unj6j0dbOzGJboXcDBWlOVS3dZLY2c/O+s6mV+cGTiHJZr4yqUYyaEwgiImiMAlv4LjroCXvgPv/Ukvn6K9KKJBnPybpjZ2LgVAfXsfWSmJpFWsAsS3n6K1EhLTIDOExLnJMFLzabeezbYciHlYrCf5GcnjEu4ONHVTnp/uqLfzCquS7OYjbeyq7YgP/4RNbrkPjSJ6vbINfnAlwGW3w8KL4Jmvwab7tfaflAFJabEeXdxhBEUY0KanPtxuRV1HH8U5qZCSpau3+vJTtFZCXkX0KrZ6VpGNI/+ETa6PnhT+igH64rjpOSS4hFd2N9DQ2R/b0h3e5FbolqiejlM72S4Kne0MAUhIgivu0jkWj98Me541jmw/GEERBqbnpDEw7KalZ4C6jv7RPhQzVmlB4R1dMYHQ2EmRWQwp2dqhfWiNfl1yfPQ+Pwh56Um09gyOZLcrpUISFGnJCSwqyeKJTdr2H1caRU65dlzbjlKwku1K48b0d0yTlApX/x3KTtSatvFP+MQIijAw0peirY96z17Z00/Qdk9PG3WI5cXDgsho5FPlWzDz1Li6SeWlJzPsVnT06W5wDZ399AwMB02282R5eS5d/Xr/uMihsLFDZD3bonZUmYineCIlEz71T5ixWveTMYwjJoJCRK4Uke0i4haR1QG2u1BEdovIPhG5NZpjDIXpVi5FVWsPjV39owX8ZqzSz57mp54WGOiKrqAA7aeo3gjNe+PK7AQe2dlWiOyBRmehsZ7YBQILM1MozEwJ7wAnw0iIrMdkoaPG+CfijbRcuPFF7eQ2jCNWGsU24OPAG/42EJEE4Dbgw8AS4JMiEpftuWyNYmt1O8NupX0UoMNPE5LHOrSjHfFkU7gA+nWOQTw5skE7s2G0jMcBOzQ2SFa2JydYgiKu/BMwWvTRdmjbyXYm4in+cLli1ukx3olJwp1SaicEbd13ErBPKXXA2vYfwKXAjogPMEQKMpJJTnDxvhXLP6JRJCbrmvhjBMVB/RwLQQGQnAUl8aVe51qFAW1BcbCxm5REF6XBWsl6MLcok5LsVFbFW0/otFxIyRk1P5pkO8MUJJ4zs2cAngHoVcDJvjYUkS8AXwCYOXNm5EfmhcslFOeksMXKDh7T2W76Stj0d12IzJUwqlHkVkR3kLagmHly3LV5HNEorAqytiM7lDamLpfw4r9/gNSk+PG9jOBZbtzkUBimIBEzPYnISyKyzcfj0nB/llLqDqXUaqXU6qKionAf3hGlOWl0WxVMi3M8bOQzVsFg92gJjdZKHYWUHDyRLKzkz9aRNos/Et3PdUCuVwVZJ8UAfZGVmkRSvCTaeZJbPuqjMFnZhilIxKaWSqnzJnmIasCzq0+ZtSwusUNiE11CYYanoNCtUaneANMWa0ERbW0CdMz4v++M/uc6IDs1kQSX0NozwOCwm8MtPXx42VHU0CenHA6+afknqqxlRlAYpg5xOP0aYR0wX0Rmi0gycDXwRIzH5Be7imxxdupYk0nBfO0XsP0UraH1oQgrInHprBMR8tKTaOkepKq1lyG3CiniKe7JnalLqPS1aY3CJNsZphixCo+9TESqgFOBp0XkeWv5dBF5BkApNQTcAjwP7AT+qZTaHovxOsHWKIqzvUIzXS6YvkJrFEMDOoY+VoIijsmzsrNDKQY4ZfCsIttRbZLtDFOOWEU9PQo86mN5DXCRx/tngGeiOLQJYwsKO1R2DDNWwTu3Qct+UG4jKHyQl64ryNo5FKEk28U9I30pjmhBYSKeDFOMeDY9TSk8TU/jmLEK3IOw8yn93giKceieFIMcbOomNz1ppPPdUUGORwOj9mrjnzBMOYygCBMz8tJwCZTn+Yhmsh3a2x7Wz0ZQjMPuSRFKjacpQ3q+rkrafsRqgWoEhWFqEV8B9VOY/IxkHvjiqSyd7qMgXfYM7bxs3KkztbNKoz/AOCfXw/R02ryjrDCbiPZT1G62ku2MoDBMLYxGEUZOnJVPerIP2SuiE+9Ah8a6zGn3Jj8jicFhXab9qPJP2OTOHK35ZXwUhimGuWNFC7tAYF4MciimAHbSHYRWDHDKkFMOQ33W67LYjsVgCBEjKKLFjBP0s/FP+CR/jKA4SjUKG6NRGKYYRlBEi+krISldV5Q1jCMvI2nk9azCKJc3iQZ2LoVJtjNMQYwzO1qk58O/bjIdtPxg96QozUn17eeZ6tghsibZzjAFOQr/kXFMlplJ+sMWFEel2QlGTU8mh8IwBTGmJ0NckJ2WhEuOYkGRUQQJKcY/YZiSGI3CEBckuISfXbGcFTNzYz2UyOBywZn/ofuoGwxTDCMoDHHD5auO8rDRs78R6xEYDBPCmJ4MBoPBEBAjKAwGg8EQECMoDAaDwRAQIygMBoPBEBAjKAwGg8EQECMoDAaDwRAQIygMBoPBEBAjKAwGg8EQEFFKxXoMYUVEGoFDkzhEIdAUpuFEAjO+yWHGNznM+CZHPI+vQilV5GvFUScoJouIrFdKrY71OPxhxjc5zPgmhxnf5Ij38fnDmJ4MBoPBEBAjKAwGg8EQECMoxnNHrAcQBDO+yWHGNznM+CZHvI/PJ8ZHYTAYDIaAGI3CYDAYDAExgsJgMBgMATnqBYWI/EVEGkRkm8ey5SLyjohsFZEnRSTbWp4kIvdYy3eKyDc99rlQRHaLyD4RuTUOx1dpLd8kIutjNL5kEbnLWr5ZRM722GeVtXyfiPxGRCTOxvea9ftush7TwjS+chF5VUR2iMh2EfmKtTxfRF4Ukb3Wc561XKzzs09EtojISo9jXWdtv1dErovD8Q17nL8nYjS+RdZv3y8iX/M6Vtj/w2EeX0T+w2FBKXVUP4APACuBbR7L1gFnWa9vAL5vvf4U8A/rdTpQCcwCEoD9wBwgGdgMLImX8VnvK4HCGJ+/fwHusl5PAzYALuv9e8ApgADPAh+Os/G9BqyOwPkrBVZar7OAPcAS4KfArdbyW4GfWK8vss6PWOfrXWt5PnDAes6zXufFy/isdV1xcP6mAScCPwC+5nGciPyHwzU+a10lEfgPh+Nx1GsUSqk3gBavxQuAN6zXLwKX25sDGSKSCKQBA0AHcBKwTyl1QCk1APwDuDSOxhcxQhzfEuAVa78GoA1YLSKlQLZSaq3S/4i/Ah+Ll/GFYxwBxlerlNpove4EdgIz0NfPPdZm9zB6Pi4F/qo0a4Fc6/x9CHhRKdWilGq1vteFcTS+iBDq+JRSDUqpdcCg16Ei8h8O4/jimqNeUPhhO6MXyZVAufX6IaAbqAUOA/+nlGpB//BHPPavspbFy/hAC5EXRGSDiHwhgmMLNL7NwEdFJFFEZgOrrHUz0OfMJlbnz9/4bO6y1P5vhcs05omIzAJOAN4FipVStdaqOqDYeu3vWov4NTjJ8QGkish6EVkrIh8L59hCGJ8/4uX8BSKa/+GQOFYFxQ3AzSKyAa0uDljLTwKGgenAbOA/RGTOFBnfGUqplcCHgX8RkQ/EYHx/Qf8B1wO/At62xhttJjK+TyullgFnWo9rwjkgEckEHga+qpQaowVaWlZM49TDNL4KpctTfAr4lYjMjbPxRYwwjS+a/+GQOCYFhVJql1LqAqXUKuB+tO0S9AX+nFJq0DJNvIU2TVQzduZZZi2Ll/GhlKq2nhuAR9FCJarjU0oNKaX+TSm1Qil1KZCLttlWo8+ZTUzOX4DxeZ6/TuDvhPH8iUgS+iZyn1LqEWtxvW2ysZ4brOX+rrWIXYNhGp/nOTyA9vmcEIPx+SNezp9fovkfDpVjUlCIFdEiIi7gf4DbrVWHgXOsdRloZ90utHN0vojMFpFk4GogLFEd4RifiGSISJbH8guAbd7HjfT4RCTd+nxE5HxgSCm1w1LBO0TkFMukcy3weLyMzzJFFVrLk4BLCNP5s77vncBOpdQvPFY9AdiRS9cxej6eAK61ootOAdqt8/c8cIGI5FkRNBdYy+JifNa4UqxjFgKnAztiMD5/ROQ/HK7xRfs/HDLh9IzH4wM9o6xFO4+qgBuBr6BnknuAHzOaoZ4JPIi2ce8Avu5xnIus7fcD/x1P40NHcmy2HttjOL5ZwG60Q+8ltCnCPs5q9IW/H/idvU88jA/IQEdAbbHO36+BhDCN7wy02WELsMl6XAQUAC8De62x5FvbC3CbdZ624hGJhTap7bMe18fT+IDTrPebrecbYzS+Eus66EAHK1ShAykgAv/hcI2PCP6Hw/EwJTwMBoPBEJBj0vRkMBgMBucYQWEwGAyGgBhBYTAYDIaAGEFhMBgMhoAYQWEwGAyGgBhBYTAYDIaAGEFhMEwQEZklHuXNDYajFSMoDIYYYVUBPuo+y3D0YQSF4ZhFRB6zKnVut6t1ikiXiPxAdGOjtSJSbC0vFpFHreWbReQ06zAJIvIn6xgviEiatf0Ka/8t1n5245rXRORXohvTfMXHmLJE5KBVSgQRybbfi8hcEXnOGvObIrLI2uYjIvKuiLwvIi95jPm7IvI3EXkL+FuET6fhKMYICsOxzA1KFw5cDfyriBSgy3msVUotR/e0+Ly17W+A163lK9FlFgDmA7cppZaiSzLYvS/+CnxDKXU8uqTFdzw+N1kptVop9XPvASldlPA14GJr0dXAI0qpQeAO4MvWmL8G/N7aZg1wilLqBHSfhf/0OOQS4Dyl1CdDOjMGgwdGHTUcy/yriFxmvS5H3/QHgKesZRuA863X56CLGaKUGgbaLS3hoFJqk8f2s0QkB8hVSr1uLb8HXaPL5oEg4/oz+mb/GHA98HnRZaxPAx6U0VYZKdZzGfCAVaU0GTjocawnlFK9QT7PYAiIERSGYxLR/bLPA05VSvWIyGtAKjCoRgugDRP8P9Lv8XoY3XkwGN2BViql3rIc5WejixNuE933u00ptcLHLr8FfqGUesLa57tOP8tgcIIxPRmOVXKAVktILEKXbA/Ey8BNACKSYGkNPlFKtQOtInKmtega4HV/2/vhr+i+GHdZx+wADorIldYYRESWe3wXu7fCdd4HMhgmixEUhmOV54BEEdmJLkW+Nsj2XwE+KCJb0SamJUG2vw74mYhsAVYA3wtxfPcBeegy6jafBm4UEbsUtd3u9btok9QGoCnEzzEYgmLKjBsMcYiIXAFcqpQKa0tWg2EiGB+FwRBniMhv0X2TL4r1WAwGMILCYIgZIvLfwJVeix9USn05FuMxGPxhTE8Gg8FgCIhxZhsMBoMhIEZQGAwGgyEgRlAYDAaDISBGUBgMBoMhIP8fIlosi9bO2wQAAAAASUVORK5CYII=",
"text/plain": [
""
]
@@ -639,8 +605,8 @@
"import matplotlib.pyplot as plt\n",
"\n",
"clustered_data -= clustered_data.mean(dim='anchor_year')\n",
- "clustered_data.sel(cluster_labels=1).plot.line(x='anchor_year', label='pos. corr')\n",
- "clustered_data.sel(cluster_labels=-2).plot.line(x='anchor_year', label='neg. corr')\n",
+ "clustered_data.sel(cluster_labels=\"lag:1_cluster:1\").plot.line(x='anchor_year', label='pos. corr')\n",
+ "clustered_data.sel(cluster_labels=\"lag:1_cluster:-2\").plot.line(x='anchor_year', label='neg. corr')\n",
"plt.legend()\n"
]
},
diff --git a/s2spy/rgdr/rgdr.py b/s2spy/rgdr/rgdr.py
index 0e505ca..00edd9b 100644
--- a/s2spy/rgdr/rgdr.py
+++ b/s2spy/rgdr/rgdr.py
@@ -79,24 +79,96 @@ def remove_small_area_clusters(ds: XrType, min_area_km2: float) -> XrType:
valid_clusters = np.array([c for c, a in zip(clusters, areas) if a > min_area_km2])
ds["cluster_labels"] = ds["cluster_labels"].where(
- np.isin(ds["cluster_labels"], valid_clusters), 0
+ np.isin(ds["cluster_labels"], valid_clusters), "0"
)
return ds
-def masked_spherical_dbscan(
+def add_gridcell_area(data: xr.DataArray):
+ """Adds the area of each gridcell (latitude) in km2.
+
+ Note: Assumes an even grid (in degrees)
+
+ Args:
+ data: Data containing lat, lon coordinates in degrees.
+
+ Returns:
+ Input data with an added coordinate "area".
+ """
+ dlat = np.abs(data.latitude.values[1] - data.latitude.values[0])
+ dlon = np.abs(data.longitude.values[1] - data.longitude.values[0])
+ data["area"] = spherical_area(data.latitude, dlat, dlon)
+ return data
+
+
+def assert_clusters_present(data: xr.DataArray) -> None:
+ """Asserts that any (non-'0') clusters are present in the data."""
+
+ if "i_interval" in data.dims:
+ n_clusters = np.zeros(data["i_interval"].size)
+ for i, _ in enumerate(n_clusters):
+ n_clusters[i] = np.unique(data.isel(i_interval=i).cluster_labels).size
+
+ if np.any(n_clusters == 1): # A single cluster is the '0' (leftovers) cluster.
+ empty_lags = data["i_interval"].values[n_clusters == 1]
+ raise ValueError(
+ f"No significant clusters found in lag(s): i_interval={empty_lags}."
+ "Please remove these intervals from the model before continuing."
+ )
+
+ elif np.unique(data.cluster_labels).size == 1:
+ raise ValueError("No significant clusters found in the input DataArray")
+
+
+def _get_dbscan_clusters(
+ data: xr.Dataset, coords: np.ndarray, lag: int, dbscan_params: dict
+) -> np.ndarray:
+ """Generates the DBSCAN cluster labels based on the correlation and p-value.
+
+ Args:
+ data (xr.DataArray): DataArray of the precursor field, of only a single
+ i_interval. Requires the 'latitude' and 'longitude' dimensions to be stacked
+ into a "coords" dimension.
+ coords (np.ndarray): 2-D array containing the coordinates of each (lat, lon) grid
+ point, in radians.
+ lag (int): The i_interval value of the input data.
+ dbscan_params (dict): Dictionary containing the elements 'alpha', 'eps',
+ 'min_area_km2'. See the documentation of RGDR for more information.
+
+ Returns:
+ np.ndarray: 1-D array of the same length as `coords`, containing cluster labels
+ for every coordinate."""
+
+ labels = np.zeros(len(coords), dtype="= 0, data["corr"] < 0]):
+ mask = np.logical_and(data["p_val"] < dbscan_params["alpha"], sign_mask)
+
+ if np.sum(mask) > 0: # Check if the mask contains any points to cluster
+ db = DBSCAN(
+ eps=dbscan_params["eps"] / RADIUS_EARTH_KM,
+ min_samples=1,
+ algorithm="auto",
+ metric="haversine",
+ ).fit(coords[mask])
+
+ cluster_labels = sign * (db.labels_ + 1)
+ labels[mask] = [f"lag:{lag}_cluster:{int(lbl)}" for lbl in cluster_labels]
+
+ return labels
+
+
+def _find_clusters(
precursor: xr.DataArray,
corr: xr.DataArray,
p_val: xr.DataArray,
dbscan_params: dict,
) -> xr.DataArray:
+ """Computes clusters and adds their labels to the precursor dataset.
- """Determines the clusters based on sklearn's DBSCAN implementation. Alpha determines
- the mask based on the minimum p_value. Grouping can be adjusted using the `eps_km`
- parameter. Cluster labels are negative for areas with a negative correlation coefficient
- and positive for areas with a positive correlation coefficient. Areas without any
- significant correlation are put in the cluster labelled '0'.
+ For clustering the DBSCAN algorithm is used, with a Haversine distance metric.
Args:
precursor (xr.DataArray): DataArray of the precursor field, containing
@@ -109,43 +181,74 @@ def masked_spherical_dbscan(
'min_area_km2'. See the documentation of RGDR for more information.
Returns:
- xr.DataArray: Precursor data grouped by the DBSCAN clusters.
+ xr.DataArray: The input precursor data, with as extra coordinate labelled
+ clusters.
"""
- orig_name = precursor.name
data = precursor.to_dataset()
data["corr"], data["p_val"] = corr, p_val # Will require less tracking of indices
+ if "i_interval" not in data.dims:
+ data = data.expand_dims("i_interval")
+ lags = data["i_interval"].values
+
data = data.stack(coord=["latitude", "longitude"])
coords = np.asarray(data["coord"].values.tolist())
coords = np.radians(coords)
# Prepare labels, default value is 0 (not in cluster)
- labels = np.zeros(len(coords), dtype=int)
+ labels = np.zeros((len(lags), len(coords)), dtype="= 0, data["corr"] < 0]):
- mask = np.logical_and(data["p_val"] < dbscan_params["alpha"], sign_mask)
- if np.sum(mask) > 0: # Check if the mask contains any points to cluster
- db = DBSCAN(
- eps=dbscan_params["eps"] / RADIUS_EARTH_KM,
- min_samples=1,
- algorithm="auto",
- metric="haversine",
- ).fit(coords[mask])
-
- labels[mask] = sign * (db.labels_ + 1)
+ for i, lag in enumerate(lags):
+ labels[i] = _get_dbscan_clusters(
+ data.isel(i_interval=i), coords, lag, dbscan_params
+ )
precursor = precursor.stack(coord=["latitude", "longitude"])
- precursor["cluster_labels"] = ("coord", labels)
+ if "i_interval" not in precursor.dims:
+ precursor["cluster_labels"] = ("coord", labels[0])
+ else:
+ precursor["cluster_labels"] = (("i_interval", "coord"), labels)
precursor = precursor.unstack(("coord"))
- dlat = np.abs(precursor.latitude.values[1] - precursor.latitude.values[0])
- dlon = np.abs(precursor.longitude.values[1] - precursor.longitude.values[0])
- precursor["area"] = spherical_area(precursor.latitude, dlat, dlon)
+ return precursor
+
+
+def masked_spherical_dbscan(
+ precursor: xr.DataArray,
+ corr: xr.DataArray,
+ p_val: xr.DataArray,
+ dbscan_params: dict,
+) -> xr.DataArray:
+
+ """Determines the clusters based on sklearn's DBSCAN implementation. Alpha determines
+ the mask based on the minimum p_value. Grouping can be adjusted using the `eps_km`
+ parameter. Cluster labels are negative for areas with a negative correlation coefficient
+ and positive for areas with a positive correlation coefficient. Areas without any
+ significant correlation are put in the cluster labelled '0'.
+
+ Args:
+ precursor (xr.DataArray): DataArray of the precursor field, containing
+ 'latitude' and 'longitude' dimensions in degrees.
+ corr (xr.DataArray): DataArray with the correlation values, generated by
+ correlation_map()
+ p_val (xr.DataArray): DataArray with the p-values, generated by
+ correlation_map()
+ dbscan_params (dict): Dictionary containing the elements 'alpha', 'eps',
+ 'min_area_km2'. See the documentation of RGDR for more information.
+
+ Returns:
+ xr.DataArray: Precursor data grouped by the DBSCAN clusters.
+ """
+ precursor = add_gridcell_area(precursor)
+
+ precursor = _find_clusters(precursor, corr, p_val, dbscan_params)
if dbscan_params["min_area"]:
precursor = remove_small_area_clusters(precursor, dbscan_params["min_area"])
- precursor.name = orig_name
+ # Make sure a cluster is present in each lag
+ assert_clusters_present(precursor)
+
return precursor
@@ -282,10 +385,11 @@ def get_clusters(
corr, p_val = self.get_correlation(precursor, timeseries)
return masked_spherical_dbscan(precursor, corr, p_val, self._dbscan_params)
- def plot_correlation(
+ def plot_correlation( # pylint: disable=too-many-arguments
self,
precursor: xr.DataArray,
timeseries: xr.DataArray,
+ lag: Optional[int] = None,
ax1: Optional[plt.Axes] = None,
ax2: Optional[plt.Axes] = None,
) -> List[Type[mpl.collections.QuadMesh]]:
@@ -296,6 +400,8 @@ def plot_correlation(
precursor: Precursor field data with the dimensions
'latitude', 'longitude', and 'anchor_year'
timeseries: Timeseries data with only the dimension 'anchor_year'
+ lag: The i_interval which should be plotted. Required if the precursor
+ has the dimension "i_interval".
ax1: a matplotlib axis handle to plot
the correlation values into. If None, an axis handle will be created
instead.
@@ -305,6 +411,15 @@ def plot_correlation(
Returns:
List[mpl.collections.QuadMesh]: List of matplotlib artists.
"""
+
+ if "i_interval" in precursor.dims:
+ if lag is None:
+ raise ValueError(
+ "Precursor contains multiple intervals, please provide"
+ " the lag which should be plotted."
+ )
+ precursor = precursor.sel(i_interval=lag)
+
corr, p_val = self.get_correlation(precursor, timeseries)
if (ax1 is None) and (ax2 is None):
@@ -326,6 +441,7 @@ def plot_clusters(
self,
precursor: xr.DataArray,
timeseries: xr.DataArray,
+ lag: Optional[int] = None,
ax: Optional[plt.Axes] = None,
) -> Type[mpl.collections.QuadMesh]:
"""Generates a figure showing the clusters resulting from the initiated RGDR
@@ -335,18 +451,30 @@ class and input precursor field.
precursor: Precursor field data with the dimensions
'latitude', 'longitude', and 'anchor_year'
timeseries: Timeseries data with only the dimension 'anchor_year'
+ lag: The i_interval which should be plotted. Required if the precursor
+ has the dimension "i_interval".
ax (plt.Axes, optional): a matplotlib axis handle to plot the clusters
into. If None, an axis handle will be created instead.
Returns:
matplotlib.collections.QuadMesh: Matplotlib artist.
"""
- clusters = self.get_clusters(precursor, timeseries)
-
if ax is None:
_, ax = plt.subplots()
- return clusters.cluster_labels.plot(cmap="viridis", ax=ax)
+ clusters = self.get_clusters(precursor, timeseries)
+
+ if "i_interval" in precursor.dims:
+ if lag is None:
+ raise ValueError(
+ "Precursor contains multiple intervals, please provide"
+ " the lag which should be plotted."
+ )
+ clusters = clusters.sel(i_interval=lag)
+
+ clusters = utils.cluster_labels_to_ints(clusters)
+
+ return clusters["cluster_labels"].plot(cmap="viridis", ax=ax)
def fit(self, precursor: xr.DataArray, timeseries: xr.DataArray):
"""Fits RGDR clusters to precursor data.
@@ -400,12 +528,18 @@ def transform(self, data: xr.DataArray) -> xr.DataArray:
data["cluster_labels"] = self._clusters
data["area"] = self._area
- # Add the geographical centers for later alignment between, e.g., splits
reduced_data = utils.weighted_groupby(
data, groupby="cluster_labels", weight="area"
)
- return utils.geographical_cluster_center(data, reduced_data)
+ # Add the geographical centers for later alignment between, e.g., splits
+ reduced_data = utils.geographical_cluster_center(data, reduced_data)
+
+ # Remove the '0' cluster
+ reduced_data = reduced_data.where(reduced_data["cluster_labels"] != "0").dropna(
+ dim="cluster_labels"
+ )
+ return reduced_data
def fit_transform(self, precursor: xr.DataArray, timeseries: xr.DataArray):
"""Fits RGDR clusters to precursor data, and applies RGDR on the input data.
diff --git a/s2spy/rgdr/utils.py b/s2spy/rgdr/utils.py
index 92d45ae..b06dd61 100644
--- a/s2spy/rgdr/utils.py
+++ b/s2spy/rgdr/utils.py
@@ -68,25 +68,45 @@ def geographical_cluster_center(
for i, cluster in enumerate(clusters):
# Select only the grid cells within the cluster
- cluster_data = stacked_data.where(
+ cluster_area = stacked_data["area"].where(
stacked_data["cluster_labels"] == cluster
- ).dropna(dim="coords")
+ )
+
+ if "i_interval" in cluster_area.dims:
+ cluster_area = cluster_area.dropna('i_interval', how='all')
+ cluster_area = cluster_area.dropna("coords")
# Area weighted mean to get the geographical center of the cluster
- # for the 0 clusters (leftovers), set to nan as this will avoid them in e.g.
- # plots
- if cluster == 0:
- cluster_lats[i] = np.nan
- cluster_lons[i] = np.nan
- else:
- cluster_lats[i] = (
- cluster_data["latitude"].weighted(cluster_data["area"]).mean().item()
- )
- cluster_lons[i] = (
- cluster_data["longitude"].weighted(cluster_data["area"]).mean().item()
- )
+ cluster_lats[i] = (
+ cluster_area["latitude"].weighted(cluster_area).mean().item()
+ )
+ cluster_lons[i] = (
+ cluster_area["longitude"].weighted(cluster_area).mean().item()
+ )
reduced_data["latitude"] = ("cluster_labels", cluster_lats)
reduced_data["longitude"] = ("cluster_labels", cluster_lons)
return reduced_data
+
+
+def cluster_labels_to_ints(clustered_data: xr.DataArray) -> xr.DataArray:
+ """Converts the labels of already clustered data to integers.
+
+ Args:
+ clustered_data: Data already clustered and grouped by cluster.
+
+ Returns:
+ Same as input, but with the labels converted to integers
+ """
+ un_labels = np.unique(clustered_data.cluster_labels)
+ label_vals = [int(lb[-2:].replace(":","")) for lb in un_labels]
+ label_lookup = dict(zip(un_labels, label_vals))
+
+ clustered_data['cluster_labels'] = xr.apply_ufunc(
+ lambda val: label_lookup[val],
+ clustered_data['cluster_labels'],
+ vectorize=True
+ )
+
+ return clustered_data
diff --git a/tests/test_rgdr/test_rgdr.py b/tests/test_rgdr/test_rgdr.py
index 3cfb548..4c719eb 100644
--- a/tests/test_rgdr/test_rgdr.py
+++ b/tests/test_rgdr/test_rgdr.py
@@ -6,6 +6,7 @@
import xarray as xr
from s2spy import RGDR
from s2spy.rgdr import rgdr
+from s2spy.rgdr import utils
from s2spy.time import AdventCalendar
from s2spy.time import resample
@@ -38,6 +39,12 @@ def example_field(raw_field, dummy_calendar):
return resample(cal, raw_field).sst.isel(i_interval=1)
+@pytest.fixture(autouse=True, scope="class")
+def example_field_multiple_lags(raw_field, dummy_calendar):
+ cal = dummy_calendar.map_to_data(raw_field)
+ return resample(cal, raw_field).sst.isel(i_interval=slice(1, 4))
+
+
@pytest.fixture(autouse=True, scope="class")
def example_target(raw_target, raw_field, dummy_calendar):
cal = dummy_calendar.map_to_data(raw_field)
@@ -146,6 +153,7 @@ def test_dbscan(
clusters = rgdr.masked_spherical_dbscan(
example_field, corr, p_val, dummy_dbscan_params
)
+ clusters = utils.cluster_labels_to_ints(clusters)
np.testing.assert_array_equal(clusters["cluster_labels"], expected_labels)
@@ -158,6 +166,8 @@ def test_dbscan_min_area(
clusters = rgdr.masked_spherical_dbscan(
example_field, corr, p_val, dbscan_params
)
+ clusters = utils.cluster_labels_to_ints(clusters)
+
expected_labels[expected_labels == -1] = 0 # Small -1 cluster is missing
np.testing.assert_array_equal(clusters["cluster_labels"], expected_labels)
@@ -186,7 +196,8 @@ def test_fit(self, dummy_rgdr, example_field, example_target):
def test_transform(self, dummy_rgdr, example_field, example_target):
dummy_rgdr.fit(example_field, example_target)
clustered_data = dummy_rgdr.transform(example_field)
- cluster_labels = np.array([-2.0, -1.0, 0.0, 1.0])
+ clustered_data = utils.cluster_labels_to_ints(clustered_data)
+ cluster_labels = np.array([-1, -2, 1])
np.testing.assert_array_equal(clustered_data["cluster_labels"], cluster_labels)
def test_fit_transform_fits(self, example_field, example_target):
@@ -198,12 +209,32 @@ def test_fit_transform_fits(self, example_field, example_target):
def test_fit_transform(self, example_field, example_target):
rgdr = RGDR(min_area_km2=1000**2)
clustered_data = rgdr.fit_transform(example_field, example_target)
- cluster_labels = np.array([-2.0, -1.0, 0.0, 1.0])
+ cluster_labels = np.array(
+ ["lag:1_cluster:-1", "lag:1_cluster:-2", "lag:1_cluster:1"])
+ np.testing.assert_array_equal(clustered_data["cluster_labels"], cluster_labels)
+
+ def test_fit_transform_multiple_lags(self, example_field_multiple_lags, example_target):
+ rgdr = RGDR()
+ clustered_data = rgdr.fit_transform(example_field_multiple_lags, example_target)
+ cluster_labels = np.array(
+ ["lag:1_cluster:-2", "lag:1_cluster:1", "lag:2_cluster:-1",
+ "lag:2_cluster:1", "lag:3_cluster:-1"])
np.testing.assert_array_equal(clustered_data["cluster_labels"], cluster_labels)
def test_corr_plot(self, dummy_rgdr, example_field, example_target):
dummy_rgdr.plot_correlation(example_field, example_target)
+ def test_corr_plot_multiple_lags(
+ self, dummy_rgdr, example_field_multiple_lags, example_target
+ ):
+ dummy_rgdr.plot_correlation(example_field_multiple_lags, example_target, lag=1)
+
+ def test_corr_plot_multiple_lags_fail(
+ self, dummy_rgdr, example_field_multiple_lags, example_target
+ ):
+ with pytest.raises(ValueError):
+ dummy_rgdr.plot_correlation(example_field_multiple_lags, example_target)
+
def test_corr_plot_ax(self, dummy_rgdr, example_field, example_target):
_, (ax1, ax2) = plt.subplots(ncols=2)
dummy_rgdr.plot_correlation(example_field, example_target, ax1=ax1, ax2=ax2)
@@ -211,6 +242,17 @@ def test_corr_plot_ax(self, dummy_rgdr, example_field, example_target):
def test_cluster_plot(self, dummy_rgdr, example_field, example_target):
dummy_rgdr.plot_clusters(example_field, example_target)
+ def test_cluster_plot_multiple_lags(
+ self, dummy_rgdr, example_field_multiple_lags, example_target
+ ):
+ dummy_rgdr.plot_clusters(example_field_multiple_lags, example_target, lag=1)
+
+ def test_cluster_plot_multiple_lags_fail(
+ self, dummy_rgdr, example_field_multiple_lags, example_target
+ ):
+ with pytest.raises(ValueError):
+ dummy_rgdr.plot_clusters(example_field_multiple_lags, example_target)
+
def test_cluster_plot_ax(self, dummy_rgdr, example_field, example_target):
_, ax = plt.subplots()
dummy_rgdr.plot_clusters(example_field, example_target, ax=ax)