From 1d0f412e02bdf4177d22d5ce14ee830cd38bfa72 Mon Sep 17 00:00:00 2001
From: "Daniel.Frisinghelli" <daniel.frisinghelli@eurac.edu>
Date: Wed, 6 Oct 2021 13:37:31 +0000
Subject: [PATCH] Refactor.
---
Notebooks/eval_precipitation.ipynb | 36 +++++---
Notebooks/pr_distribution.ipynb | 135 +++++++++++++++++++++++++----
Notebooks/pr_sampling.ipynb | 31 +++++--
3 files changed, 163 insertions(+), 39 deletions(-)
diff --git a/Notebooks/eval_precipitation.ipynb b/Notebooks/eval_precipitation.ipynb
index a0c652d..40f103a 100644
--- a/Notebooks/eval_precipitation.ipynb
+++ b/Notebooks/eval_precipitation.ipynb
@@ -59,12 +59,13 @@
"PLEVELS = ['500', '850']\n",
"# PLEVELS = []\n",
"SPREDICTORS = 'p'\n",
- "DEM = ''\n",
- "DEM_FEATURES = 'dem'\n",
+ "DEM = 'dem'\n",
+ "DEM_FEATURES = ''\n",
"DOY = ''\n",
- "WET_DAY_THRESHOLD = '2'\n",
+ "WET_DAY_THRESHOLD = '1'\n",
"# LOSS = 'MSELoss'\n",
- "LOSS = 'BernoulliGammaLoss'"
+ "LOSS = 'BernoulliGammaLoss'\n",
+ "SEASON = 'season'"
]
},
{
@@ -202,6 +203,7 @@
"PATTERN = '_'.join([PATTERN, DOY]) if DOY else PATTERN\n",
"PATTERN = '_'.join([PATTERN, '{}mm'.format(str(WET_DAY_THRESHOLD).replace('.', ''))]) if WET_DAY_THRESHOLD else PATTERN\n",
"PATTERN = '_'.join([PATTERN, LOSS])\n",
+ "PATTERN = '_'.join([PATTERN, SEASON]) if SEASON else PATTERN\n",
"PATTERN"
]
},
@@ -266,14 +268,26 @@
"metadata": {},
"outputs": [],
"source": [
- "# align datasets and mask missing values in model predictions\n",
- "if LOSS == 'BernoulliGammaLoss':\n",
+ "# align datasets\n",
+ "if len(y_pred.data_vars) > 1:\n",
" y_true, y_refe, y_pred_pr, y_pred_prob = xr.align(y_true.precipitation, y_refe.precipitation, y_pred.precipitation, y_pred.prob, join='override')\n",
- " y_pred_prob = y_pred_prob.where(~np.isnan(y_true), other=np.nan)\n",
"else:\n",
- " y_true, y_refe, y_pred_pr = xr.align(y_true.precipitation, y_refe.precipitation, y_pred.precipitation, join='override')\n",
- "y_pred_pr = y_pred_pr.where(~np.isnan(y_true), other=np.nan) \n",
- "y_refe = y_refe.where(~np.isnan(y_true), other=np.nan)"
+ " y_true, y_refe, y_pred_pr = xr.align(y_true.precipitation, y_refe.precipitation, y_pred.precipitation, join='override')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "aa3f8e67-6b49-46d4-a956-a0cee7b3923a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# mask missing values\n",
+ "mask = ~np.isnan(y_true)\n",
+ "y_pred_pr = y_pred_pr.where(mask, other=np.nan) \n",
+ "y_refe = y_refe.where(mask, other=np.nan)\n",
+ "if len(y_pred.data_vars) > 1:\n",
+ " y_pred_prob = y_pred_prob.where(mask, other=np.nan)"
]
},
{
@@ -523,7 +537,7 @@
"# fig.suptitle('Average yearly mean absolute error: 1991 - 2010', fontsize=20);\n",
"fig.subplots_adjust(hspace=0, wspace=0, top=0.85)\n",
"\n",
- "# add colorbar for bias\n",
+ "# add colorbar for dem\n",
"axes = axes.flatten()\n",
"cbar_ax_bias = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,\n",
" 0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])\n",
diff --git a/Notebooks/pr_distribution.ipynb b/Notebooks/pr_distribution.ipynb
index 51d2351..1a7f100 100644
--- a/Notebooks/pr_distribution.ipynb
+++ b/Notebooks/pr_distribution.ipynb
@@ -24,7 +24,7 @@
"import xarray as xr\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
- "import seaborn as sns\n",
+ "import seequantilesn as sns\n",
"import pandas as pd\n",
"import scipy.stats as stats\n",
"from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
@@ -159,6 +159,18 @@
"gamma = stats.gamma(alpha, loc=loc, scale=beta)"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "dcd9bfeb-67dc-4b63-98fd-c86c3a07c2b0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# fit lognormal distribution\n",
+ "alpha, loc, beta = stats.lognorm.fit(y_valid, floc=0)\n",
+ "lognorm = stats.lognorm(alpha, loc=loc, scale=beta)"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -171,6 +183,30 @@
"genpareto = stats.genpareto(alpha, loc=loc, scale=beta)"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d489a3e7-7ece-440e-bbd9-1cfd739d822c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# fit exponential distribution to data\n",
+ "loc, beta = stats.expon.fit(y_valid, floc=0)\n",
+ "expon = stats.expon(loc=loc, scale=beta)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "01d8c7d9-541e-481d-b0de-e8590c571ca5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# fit weibull distribution to data\n",
+ "alpha, loc, beta = stats.weibull_min.fit(y_valid, floc=0)\n",
+ "weibull = stats.weibull_min(alpha, loc=loc, scale=beta)"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -182,12 +218,18 @@
"eq = np.quantile(y_valid, quantiles)\n",
"tq_gamma = gamma.ppf(quantiles)\n",
"tq_genpareto = genpareto.ppf(quantiles)\n",
+ "tq_expon = expon.ppf(quantiles)\n",
+ "tq_lognorm = lognorm.ppf(quantiles)\n",
+ "tq_weibull = weibull.ppf(quantiles)\n",
"\n",
"# Q-Q plot \n",
"RANGE = 40\n",
"fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n",
- "ax.scatter(eq, tq_gamma, color='grey', label='Gamma')\n",
- "ax.scatter(eq, tq_genpareto, color='k', label='GenPareto')\n",
+ "ax.scatter(eq, tq_gamma, marker='*', color='k', label='Gamma')\n",
+ "ax.scatter(eq, tq_genpareto, marker='x', color='k', label='GenPareto')\n",
+ "ax.scatter(eq, tq_expon, marker='o', color='k', label='Expon')\n",
+ "ax.scatter(eq, tq_lognorm, marker='+', color='k', label='LogNorm')\n",
+ "ax.scatter(eq, tq_weibull, marker='^', color='k', label='Weibull')\n",
"ax.plot(np.arange(0, RANGE), np.arange(0, RANGE), '--k')\n",
"ax.set_xlim(0, RANGE)\n",
"ax.set_ylim(0, RANGE)\n",
@@ -233,6 +275,9 @@
"# fit distribution to observations for each month\n",
"month_gamma = {}\n",
"month_genpareto = {}\n",
+ "month_expon = {}\n",
+ "month_lognorm = {}\n",
+ "month_weibull = {}\n",
"for month, idx in month_idx.items():\n",
" print('Month: {}'.format(calendar.month_name[month]))\n",
" # select the data of the current month\n",
@@ -249,7 +294,22 @@
" # genpareto\n",
" alpha, loc, beta = stats.genpareto.fit(data, floc=0)\n",
" genpareto = stats.genpareto(alpha, loc=loc, scale=beta)\n",
- " month_genpareto[month] = genpareto "
+ " month_genpareto[month] = genpareto \n",
+ " \n",
+ " # exponential\n",
+ " loc, beta = stats.expon.fit(data, floc=0)\n",
+ " expon = stats.expon(loc=loc, scale=beta)\n",
+ " month_expon[month] = expon\n",
+ " \n",
+ " # lognormal\n",
+ " alpha, loc, beta = stats.lognorm.fit(data, floc=0)\n",
+ " lognorm = stats.lognorm(alpha, loc=loc, scale=beta)\n",
+ " month_lognorm[month] = lognorm\n",
+ " \n",
+ " # weibull\n",
+ " alpha, loc, beta = stats.weibull_min.fit(data, floc=0)\n",
+ " weibull = stats.weibull_min(alpha, loc=loc, scale=beta)\n",
+ " month_weibull[month] = weibull"
]
},
{
@@ -276,10 +336,16 @@
" # compute theoretical quantiles\n",
" tq_gamma = month_gamma[month].ppf(quantiles)\n",
" tq_gpare = month_genpareto[month].ppf(quantiles)\n",
+ " tq_expon = month_expon[month].ppf(quantiles)\n",
+ " tq_lognr = month_lognorm[month].ppf(quantiles)\n",
+ " tq_weibu = month_weibull[month].ppf(quantiles)\n",
" \n",
" # plot empirical vs. theoretical quantiles\n",
- " ax.scatter(eq, tq_gamma, color='grey', label='Gamma')\n",
- " ax.scatter(eq, tq_gpare, color='k', label='GenPareto')\n",
+ " ax.scatter(eq, tq_gamma, marker='*', color='k', label='Gamma')\n",
+ " ax.scatter(eq, tq_gpare, marker='x', color='k', label='GenPareto')\n",
+ " ax.scatter(eq, tq_expon, marker='o', color='k', label='Expon')\n",
+ " ax.scatter(eq, tq_lognr, marker='+', color='k', label='LogNorm')\n",
+ " ax.scatter(eq, tq_weibu, marker='^', color='k', label='Weibull')\n",
" ax.plot(np.arange(0, RANGE), np.arange(0, RANGE), '-k')\n",
" ax.set_title(calendar.month_name[month], fontsize=14)\n",
" ax.set_xlim(0, RANGE)\n",
@@ -290,7 +356,7 @@
" ax.set_yticklabels([str(t) for t in np.arange(0, RANGE + 5, 5)], fontsize=12)\n",
"\n",
"# add legend\n",
- "axes[0].legend(frameon=False, fontsize=12)\n",
+ "axes[0].legend(frameon=False, fontsize=12, loc='upper left')\n",
"\n",
"# add figure title\n",
"fig.suptitle('Reference period: {} - {}'.format(str(PERIOD[0]), str(PERIOD[-1])), fontsize=14)\n",
@@ -330,6 +396,9 @@
"# iterate over the grid points\n",
"gammaq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan\n",
"genpaq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan\n",
+ "exponq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan\n",
+ "lognrq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan\n",
+ "weibuq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan\n",
"for i, _ in enumerate(y.x):\n",
" print('Rows: {}/{}'.format(i + 1, len(y.x)))\n",
" for j, _ in enumerate(y.y):\n",
@@ -351,18 +420,42 @@
" alpha, loc, beta = stats.genpareto.fit(point, floc=0)\n",
" genpa = stats.genpareto(alpha, loc=loc, scale=beta)\n",
" \n",
+ " # fit Exponential distribution to grid point\n",
+ " loc, beta = stats.expon.fit(point, floc=0)\n",
+ " expon = stats.expon(loc=loc, scale=beta)\n",
+ " \n",
+ " # fit LogNormal distribution\n",
+ " alpha, loc, beta = stats.lognorm.fit(point, floc=0)\n",
+ " lognr = stats.lognorm(alpha, loc=loc, scale=beta)\n",
+ " \n",
+ " # fit Weibull distribution\n",
+ " alpha, loc, beta = stats.weibull_min.fit(point, floc=0)\n",
+ " weibu = stats.weibull_min(alpha, loc=loc, scale=beta)\n",
+ " \n",
" # compute theoretical quantiles of fitted distributions\n",
" tq_gamma = gamma.ppf(quantiles)\n",
" tq_genpa = genpa.ppf(quantiles)\n",
+ " tq_expon = expon.ppf(quantiles)\n",
+ " tq_lognr = lognr.ppf(quantiles)\n",
+ " tq_weibu = weibu.ppf(quantiles)\n",
" \n",
" # store theoretical quantiles for current grid point\n",
" gammaq[:, j, i] = tq_gamma\n",
" genpaq[:, j, i] = tq_genpa\n",
+ " exponq[:, j, i] = tq_expon\n",
+ " lognrq[:, j, i] = tq_lognr\n",
+ " weibuq[:, j, i] = tq_weibu\n",
"\n",
"# store theoretical quantiles in xarray.DataArray\n",
"gammaq = xr.DataArray(data=gammaq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),\n",
" name='precipitation')\n",
"genpaq = xr.DataArray(data=genpaq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),\n",
+ " name='precipitation')\n",
+ "exponq = xr.DataArray(data=exponq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),\n",
+ " name='precipitation')\n",
+ "lognrq = xr.DataArray(data=lognrq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),\n",
+ " name='precipitation')\n",
+ "weibuq = xr.DataArray(data=weibuq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),\n",
" name='precipitation')"
]
},
@@ -375,7 +468,21 @@
"source": [
"# compute bias in theoretical quantiles\n",
"bias_gamma = gammaq - equantiles # predicted - observed\n",
- "bias_genpa = genpaq - equantiles"
+ "bias_genpa = genpaq - equantiles\n",
+ "bias_expon = exponq - equantiles\n",
+ "bias_lognr = lognrq - equantiles\n",
+ "bias_weibu = weibuq - equantiles"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "23abd0d1-7c27-4f02-b7ae-9165c2dde0b6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# distributions\n",
+ "dists = {k: v for k, v in zip(['gamma', 'genpareto', 'expon', 'lognr', 'weibu'], [bias_gamma, bias_genpa, bias_expon, bias_lognr, bias_weibu])}"
]
},
{
@@ -390,8 +497,7 @@
"fig, axes = plt.subplots(3, 3, sharex=True, sharey=True, figsize=(12, 12))\n",
"axes = axes.flatten()\n",
"\n",
- "for dist in ['gamma', 'genpareto']:\n",
- " biasq = bias_gamma if dist == 'gamma' else bias_genpa\n",
+ "for dist, biasq in dists.items():\n",
"\n",
" # iterate over quantiles to plot\n",
" for ax, q in zip(axes, plot_quantiles):\n",
@@ -410,17 +516,8 @@
" cbar_bias.ax.tick_params(labelsize=14, pad=10)\n",
"\n",
" # save figure\n",
- " fig\n",
" fig.savefig('./Figures/pr_distribution_{}_grid.png'.format(dist), bbox_inches='tight', dpi=300)"
]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a5ee4d3f-608b-4598-b235-3cd20a184aff",
- "metadata": {},
- "outputs": [],
- "source": []
}
],
"metadata": {
diff --git a/Notebooks/pr_sampling.ipynb b/Notebooks/pr_sampling.ipynb
index 5d23a94..9d2088e 100644
--- a/Notebooks/pr_sampling.ipynb
+++ b/Notebooks/pr_sampling.ipynb
@@ -10,7 +10,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 1,
"id": "eeba7f9b-066a-4843-bd64-5b6326c0b536",
"metadata": {},
"outputs": [],
@@ -41,7 +41,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 2,
"id": "e75b3217-26f7-4a4a-ae2a-4fbb92a9f2a2",
"metadata": {
"tags": []
@@ -54,7 +54,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 3,
"id": "3aa8466e-84a9-4c2e-ae19-403b6246e27f",
"metadata": {},
"outputs": [],
@@ -65,7 +65,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 4,
"id": "4f1a58a2-8c4c-4d73-a116-e64e68fdd507",
"metadata": {},
"outputs": [],
@@ -76,7 +76,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 5,
"id": "5e6696df-8660-4083-9a32-0dd282112948",
"metadata": {},
"outputs": [],
@@ -88,7 +88,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 6,
"id": "b87accd6-d5e4-4dc6-9532-3ef8aa162d24",
"metadata": {},
"outputs": [],
@@ -99,7 +99,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 7,
"id": "559d1450-09db-4b2f-844a-d572485973e0",
"metadata": {},
"outputs": [],
@@ -111,10 +111,23 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 8,
"id": "7fd013f9-77d0-48de-8d5f-2c6a1cb3ed17",
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnsAAAKQCAYAAADnt7zWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABPJklEQVR4nO3dedxcZX3//9eHRCCQhDUaBExQUCpiRYJgQRaxWBWpyKagLLYuLF+g1gV3oP2JuFtZLCJE0LSCC1ZpVaAmYEEMiBShKAgh7FsIJAEMJJ/fH+dM7slkZu79nnDl9Xw85nHuOee6zrlm5sx1v+eskZlIkiSpTGv1ugGSJEkaPYY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9rSKiJgdERkRR3aYvlNE/CQiHomI5XXZkwdSd7RExMzmdmhwOn1uETG9Hu81mp6jImLP+jOc1+u2jJRe9TOjpfEdi4jpz5Xl2uc+t4zvdQOeSyJiPPAu4B3AXwKbAEuAB4A7gCuB/87MuS31XgW8DZiXmTPHrsUrq7/QRwILM/OrQ5zHNsBsYD1gOfBIPVw8Em2U1hQRsSFwIkBmntzLtqwufE/KFxFrAXsAOwEz6uH0evLRmfmNHjWtaIa9AYqIKcB/Uq2cDU8DAbwM2BZ4M/A4sGFL9VcBnwHmADNHt6VdTa/bcRfw1S7l5gN/oHotrd5HFfSuAvbLzIWDqKvnnmeoPk+NvA2pvo8AJ/euGauVDfE9Kd1k4L973Yg1jWFv4L5DFfQWAf8EXJiZDwBExCRgZ2B/4C09a+EIyczDu0zerh5e1Cbo9VdXzzGZeS/VDxlJGilLgN8C1wFzgS8DU3vaosIZ9gYgIrYF9qmfviczv988PTMXAZcDl0fEh8a6fWNsQj10t60kabAeByZn5vLGiIg4rYftWSN4gsbAbN/090+7FczMp5qf1we2n18/3aPpgNjGY8/mso2DZSPiLyLi2xFxd0Q8ExGXNJV7cUT8Y0RcERF3RsTTEbEwIn5dj59Ai/rg7F/WT6e1aceRTWVXOfg5IubVr6XR3vOb6s7rVrelHWtFxLsj4rKIeDgilkbEfRHxvYjYudt7GxE71yeGLIiIxRHxu4g4oT4GZMgi4m8j4j8j4sH6vV4QEX+IiH+LiEPalB/0+1/XO7J+b2bXz98ZEVdHxBP1e/GjiPiLpvKbRcTX6/f+6Yi4PSJOiohxbea90okUEbFrRPy0nu+T9Xt13GDfq9b5tkxbcYB2RIyLiBMj4sZ6eQvq5c9oN9+meewaEZfW5ZfU9U+s15MRPwA8Ir5Vz/PzbaYd2LROf6TN9A80f35tpu8WEf8eEfdExJ8j4tGIuLz+nKOl7Gzgzqbnrd/Hkwf5ujaIiC82rY93R8Q3I2KLfuptHBFHRMQPIuLWiFhUfw63RMSXI+KFbeocXrfx/qiOY+40773qck9GxAYDeA2zGcJ7EhET6nXwDxHxVEQ8VH8O23Qo37zerhMRn4iI/61fe0Z13GCj7JD6q4jYIyK+X68LSyPi8Yi4LSIuiYj3d/seRsSL6s+usR7dWX+2k/t5/94eET+r2/nnuv53I+LV3er1M88R73Ozsrz/koMXTf9/ImJyRHw+Iv5Urxd3RMSpEbFuU/m9I+LnUZ1suCQiroyI13WY98n1vGdG5diIuKGud39U/6+3aCq/TT3unvo7+fuIeO9ovO4ByUwf/TyAg4CsHy8ZZN0HqH7JJLC0ft78+Kumso1lvJtqM3cCTwBPAZc0lbuuqexy4LF62Bg3F5jU0o65wIJ6+rI27TikqezsutyRLfUfqF9D1q+pUXdut7pN0yYBl7W0/fGm58uA4zq8j+8Anm0q+xjV8WQJfB/4dv33yYP8fP6/pnk2v9+N5w+0qTPo97+ud2Q9fTZwev33M/UyG3UfBV4KbAPc3dSm5td+Zpt5T2+afkDTe9P8PiXwI2B8m/ptP7fm+bapM7Oe9s/Af9G3ji9qWt5TwGs7vPeH1595u8/0h03zX+UzbZo2b5Cf9+F1vWvbTPuXprZc2mb6rHraKW2mnd5Ut/GZNb+2fwPWair/Q+Dh5vWs5fGhQbymzYDbWt7zxmfwEPB3nd4r4Ist7X68ZV17CHhlS50JwMJ6+lu7tOvCusx3Bvg6BvyeNK2vx1PtDkyqY6ifbKr/KG3666Z153PAtU3rbeM1bTic/orquObm93QJ1Z6Q5nHrttRpjP/but2Ndaj5uzsXeF6b5a1FX/+X9ef3WEs7j+7wnjfKTB+rPrdDO+bV8/rAMOfTWC/+Afi/+u/F9P3fSuA/6rLH1J/pspbP9c/Arm3mfXI9fSbw701lmz/bO4ApwC5Nn8FCVv7/8OHhvl9Dem96sdDn2gN4cdMH9TNgyiDrH1nXnd1PucYyFtUr7Svq8UFTpwV8EzgBeAmwdj1uHeCtVAfTJ+0DwZ4M4B8k3QNbx2kDqPujetqNVCezTKjHbwh8rP7iLGv9otWvsxHAfg68uB6/HvBBqg5p4WA7Hqog0/iH/Flg06Zpz6cKTd9qU2+o739jPVhI1fmcAKxXT9seuJW+oHMtcDXwl02v9RP0/dN5RZvXkk3z/y9gq3ra+sCHm17rxwf6uTGwsPcY1T+og5vej1cCN9XTf9Om7rb1553ApdT/bKiCxP+j+qfS6CxX+UwZethrvJ5ngIkt026k75/s48C4lun31tNf3zL+BPqC0dH0hYV1qX4o3ldP/9hA39tBvqZf1PN5GNiPOlQCrwNup++7scp7RfVP8TRgh8b7AYwDdqTq6xL4PRAt9c5qrKsd2jSZvh+srx/EaxnQe9K0vj5GtTXwjXW716pfd+OH0kVd1p1Fdf1DmtbbadSBiiH0V1Tf00bQ/hawZdO0jYG/ofrRsHZLmxrf3ceAK+jr+9cB3kMVZBM4ps3rOYm+fuGT1D80gc2Bi+gLfLu3qdtY7vSW8aPS53b5POcxsmFvIVV/uls9fm3g7+kLq5+i6oM/S9/3dRpVn9upzzq5ad6LgMPq+Ua9zt1fT/9G/Xp+0vS+TQbOpu/H2CbDfc8G/d6M9QKfqw9W/uX0Z6pj9P6Z6pdY1/DH4MPen6g7liG088X1Cr2EOkg0TduTHoU94A31+DuBjTvU/Uhd5qct479Vj7+Vll/E9fRPNr13Jw/ivTq4rvN/I7iedHv/j2xq52fa1H1d0/QFjU6opcwV9fRPt4yf3lT398A6beo2OqvHgfUH+LmtmG+b+c1sWuZubabv2DR9Wofv0020/ONrWRfafqYMMezVde+q676xadzGVP8sb6Hvn/yOTdO3pu+7P6Fp/IZUHf8zwGs6LG+Xet4Lml9rt/d2EK+leZ3Zq830rekLCoN6r6iCxs113T1apu1Qj19Km/6Pvq1bd9ASFPtZ5oDek6b19Ulg6zbTD6inP926frWst/t0mP+Q+ivgNfW4xbT8WOjn9fT33f16Pf2/W8avT99WqdPa1BtHdeWEBK7sstzpLeNHpc/t8vrnMbJh75kO68W3mtp9Xpvp0+jbCveilmknN9U9ok3ddzdN/yMte1Cofog0tsAfPtz3bLAPj9kbuPdSnTG0lCrN7021peUS4KGI+E1EHBax8rE5Q3RGthz7N1CZeQdVB70e1SVfVhdH1MOZmbmgQ5lZ9XCvqI9Lq9/Pt9fjv5KZT7ep91WqTn+wnqiHG0TEekOov4oBvv9LqdalVv9D9c8J4Oxsc7YzVdgDeEWXZnwpM//cZvyX6/lPBv66S/3Buiozf9U6MjOvB+6pnzbO4m5cZ+tt9dOvZubSNvM8gyowt5WZR2ZmZOb0IbT3ynq4R9O411H9Qp/dYXrj77kt380DgInArzLzNx3a+muq0LMRVQAeSQfWw19n5i/bLPt24HtDmXG9Dl1WP921ZdoNVLtPn0d17dFWR9XDmVn/pxsl369fY6v/oPqnug5V4G3nfzPzFx2mDam/oq9PeR7VdVgH68sdvruX1MPW7/0+VN/npcAqx6Fm5jKqq0cAvC4i+j3jdZT73LFycYf14vKmv1c5KSQz76LaGg6d+9h7qA5R6DbvL2Tmsy3zXk7fcfPd+u9RYdgboMxcmpn/CGwJfIDqGJxGSofqwpDfAb43nINXa9f0VyAi/jqqEwj+VB8AveJAZqoLPgOscnB1D/1VPfyHiHig3YPqWDioglKjo3wxfdctnNNuxpm5GLh+CG26lmpry2bANRHxvojYaiAVh/n+z8vqDO7W19G4SDVUv/DbebAebtSlebPbjczMJ4Ab6qdDPmi7jbldpt1bD5vb+2Kqf1AAq4REgMx8kqF9pgPRWI/ahbk5A5jerLFe79xpva7X7RfV5bYcgfY3a3yObb8bA5hGRGwbEWdEdZLCE9F3V5yk2kUN7dflc+vhUc0jozrJqLE1c2Z/L2CY2q57mfkM1W516Pxd6dbPDrW/uq1+rE3Vp/xD/f4OdCNAp+9Su+8R9H3+N2bmYx3qXkm127W5fDej2eeOlZs6jG+sE0/TF+pa9dfH3pLtTzB5qOnv4fTfo8JLrwxSZj4E/Gv9ICJeQHWs1qepOvKDqLbQfG0Yi3m428SI+Beq45oanqEKLc/Uzzem+mW5/jDaMNI2q4cb1I/+NLa0TWkad1+X8vd2mdZWZj4WEe8Gvkt1jFnjM32A6jio8zJzlc5uBN7/+7s0a1k/ZRrTn9dlHt3ei8a0KV3KDNYqwbVJY6tAc3s3bfq723vR7fMejsaWu50iYr06WDaHuYeodo29LiLWqjv2TmGvsV5PoO+yRN2MyBbkJo3PcUjfjYh4B3ABfZ9P4ySExtaliVTrcbt1+btUJ3hsHxE71ltyoTrGDODyzJzf7ysYnsGue8269bND6q8yc1lEHEq1Je7FVFvTvwwsiIj/ptoi9JMuWzs7vZ7Ga2n9n934/Dt+xpn5dEQ8CryAgX3vR63PHUP99Z8PdvkM+utj2867/uwHuvxu/feocMveMGXmg5l5LtUvpkZqf0+XKgOxrNOEiHgTVdBYRnUMwdZUx3hskplTM3Mq1RYrqHZLrS4a69rf1rvf+nvMG+T8h/RaM/M/qY4Teh/Vwcz3UV3c83BgdkScs9JCnrvvf8Pq0KaetiEz/0jVGT8PeG1UlwV5FfCHzHygDne/ovr1vX1Utxl8EdXWkatbZtdYr78ywPV65ui/wlW0fb+juivQN6neh+9RXTR+3czcqGld/kqnedRbii+unx5Vz7NxS0mA80bsFYyOjv0sw+ivMvM6qrPp30UVpO+g+gF4IPBj4NJoc/mkYVpnhOc3EKtDX6IBMuyNkMx8hOqLDNWlM0bLQfXw3Mw8JTP/1OYXygtGcflD1QjCLx9kveZf3912S2/WZVpXmfl4Zn4zMw/JzM2pji/7Zj35vRHRfFeU58L7P5D3qevW41HWvOxun9uQP9MBaD4ubzeqvnB20/TmXbmNrXq/rXdfNRvqej1SGu/lUL4bb6LacncLcGhmXl/v/mzW37rc2JV7aESsQ3XW6lSqLd2X9FN3dTaszzUzn8rM72bmEZn5EqqtfKdRHfbzJqpDgUZC4/Of1qlAfV25xm7mgXzvR73P1dgz7I2sxgHlrQecN/bvj8QvocZFG29oNzEiptH5gOSRbMdgNY6POWCQ9e6gOtUdYPd2BSJifVa+Z/GwZOYtmfk+4Nf1qOZjt4bz/o+VPdqNjOq2fo1jdn47ds1ZxR30Hci+W7sCUV2YeqRPZmjWHPba7aLtb3pDY73eIyIGe0B+8x0EhvqdbHyObb8btbbrA33r8v+2OwapbtPruy28PjHnVqqtoG+j7/i9WR1ONOjPSLwnI2Go/VVbmXlnZn6cvpNlOn0mg9X4/LeJiM07lNmdvt2/A/nej3mfq9Fn2BuAiNgqIl7ST5n16DvD8Hctkxv/2DYcgeY8Xg+37zD9s3QOcyvOPh2BdgzWzHo4IyK63j83IlYcvFpvNftB/fTEeutBq+MZwrFQEbF2P0UaZ102L3M47/9Y+ccOr+1Eqmu/PUF1TGJP1MGisRX8hIhod/zKMVRbnUZLI7jtTHWNtuZxUB18vpjqn92ebaY3XEz1I29d4AvdFti8XteeaPp7w/4a3EFjN+prI2KVf8wR8WKq68i101iXX9EhWL2X6npr/Wls3fsgffcGH+ou3JF4T0bCzHo4qP5qiH3KcPyC6j17HtW1NFvbNo7qmnJQnTX/QH8zHM0+V71j2BuY7YA/RMQPI+LgiFix+Toi1o+It1Jdy6hxJmfryRk318OXRz+3BBuAxqUQ3h8R72l0LlHdYufbwDupLszZzm1UJxFsEBEj8ot1oDLzZ1QXCwY4LyJOaXkfN4rqtmU/ZtXLkpxGdYDyXwCXNM6Yjeo2SSdSXVrgcQbv6KhulXNoS1s2jIiP0/dP/udNdYbz/o+VFwE/qo81IyLWi4gPAp+pp59en5TQS6dRbQHfHvhBvUWUiFg3Io6tpy/sVDn6bnk1b4jLv4XqzOd1qE7OuS0zVxyMXl824Wqqk0m2otri1O7yMo9SXWAX4KiIuCgiVlxWoX49u0XEmVQnbjXXXUjfAfArndE6UPWWtcY6+f2I2DfqqwFExK5UF0butIXtcqrdiq8A/iXq24RFdZupDwNnUl0suz8XUH2Wr6EKHTdkdWmWobyehQzzPRkJw+iv3hwR10TEexvrdF1+vahulXVYPaq5TxlOO5dQ/cAEOD6qW79NrJe5OdVVI3aj74LLAzVafS71fDaIiE0bD/qyyPrN4zsETQ1FjvGF/Z6LD6pf/tnyeJK+K4g3Hs/S5u4E9TzmNJV7lOoikvOAXZrKtL3AZct81qbaxdC8zMeann+K7hdFbr449MKmdhzYVKZb/Y7TBlB3ffouWNvchsdbxp3fpm63W/f8gCHcuodqS1fzche3vJcJ/OtIvf8M4OLa9F1cdM8O09vOg4HfLu0SRv52aR3f837Wh6NY+TZCC+i7rdFFTZ/px7ose94wvtc/bFr2OW2mf7xp+vX9zOuTLa9lSf16mm+Zdmebeqe0rH/z6seJg3gdrbdLe5KB3y7ty031Gp9B43v2M6oLxyfV9ea6teHipnm0veXhIF5Pv+9Jt/Wqv+/SQNbbutyg+yuqvTvN056s39PmdeNSVr3gbte+n+7fw3Gseru05mUuo82dN/pbLqPQ57b5/Pp7dPx8+5lv23oM4MYCneZB0+3SutTt73Psdx6j9XDL3gBk5s+BlwEfovpn2bg+z0SqL/9vqS4y+ZeZ+dlV5wBUF6k8i+qK7BOpDqidRrX7ZzBtWUp1dffPUR1bsZzqC3kZ1X0q/6lLdagODD6N6rZe6zS1YzR3mQGQmUsyc39gX6p/tPdSXa5ibar3dBbVGWvHtKn771QXdr2U6j1fm2rrzIn03bt4sGZR7ar6HtV9FJ+heh/up7oo699m5vtb2jHc93/UZeYPgL2o3qtlVO27keos4rdny8U+eyUzz6faTfozqn+g61B9pidQ/aPZoC66cJSaMKfD3+3GXdlm+gqZ+c9U11c8hyp4BVVYuJ/q1nVHU+0ybnUq8FHgf+s6je/jhgN5AfWy76e6zueXqe4OMo7q/fwW1TGaf+pS94NUZ6LfQLUFcDzVYSgnUu2SHei60tgK9meqS7IMx7Dfk5EwxP7qv6nupPBtqmu9PUl1j91HqbakHkHVT4zYdzAzl2XmEXVbfkH1fWn0Y/9GdWeXs4Yw39Hoc9UjUadNSc9h9S7bOwEys9fHDA5bfQzZXVTXrtwrM2f3tkXqJiK+SXXv0e9l5jt63R5JK3PLnqTV0Tuogt4TQNvbkGn1UF+nsBHwzulWVlJveAcNST1RnwSziOrQiHszc3l9ZuPh9N238qzs/ckk6qA+QenLVLsN/5e+e39KWo24G1cqwHNxN25EfIe+sxOXUp3UsCF9l665nOr4pnY3YlcPRcSBVLdK25Tq2MQE9snMy7tWlNQTbtmT1CtnUe2m3Y3qjNINqc4i/F/gO8AFq8vJJFpF4ySzP1Od3HGqQU9afbllT5IkqWCeoCFJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe9JqLCJmR0RGxJEdpu8UET+JiEciYnld9uSB1B0tETGzuR2Sxo59htoZ3+sGqDciYiZwRP30emCnzMwOZb8DHAZ8OzOPHME2nFz/+dXMXDhS8+1nma8C3gbMy8yZY7HMDu2YDhwJLMzMrw5xHtsAs4H1gOXAI/Vw8Ui0USpNRIwH3gW8A/hLYBNgCfAAcAdwJfDfmTm3qc6rsM/Qc5xhTwA7AvsDPxzj5X6mHs4EFo7RMl9VL3dOvdxemV634y7gq13KzQf+ADzeZtr7qDrtq4D92gTmbnWlNUpETAH+E5jRNPppIICXAdsCb6b6vmzYVOZV2GfoOc6wp4ZTI+KSzFze64aoT2Ye3mXydvXwonZbRvupK61pvkMV9BYB/wRcmJkPAETEJGBnqh+9b+lZC0eAfYba8Zg9zQGepOoEDu1xWzQ4E+qhu2CkLiJiW2Cf+ul7MvMLjaAHkJmLMvPyzDwW+IueNHJs2GesoQx7egA4o/775PqYlkGLiLdHxM8i4uGI+HNE3BMR342IV7cpOzMimo8PvLM+OLfxmDmI5U6KiE9FxPURsSgilkbEfRFxXUR8ISJe0VQ2gfPrp3u0LDMjYs/msvVjekT8RUR8OyLujohnIuKSpnIvjoh/jIgrIuLOiHg6IhZGxK/r8Y3OtbnN84Bf1k+ntWnHkU1lVzlgOiLm1a+l0d7zm+rO61a3pR1rRcS7I+Ky+nNrvHffi4id+3nfd64P8l4QEYsj4ncRcUJE2KdodbR9098/7VYwM59q/G2fsUo77DOeo9yNK4DTgQ8ALwGOAr450Ir1F/V8oLH5fxnVbpLNqbYUviMijsvMs5uqPQ48CLygfv5IXa95+kCWvQFwNfDyetTyuu4LgM2ojkVcBpxUT3+Q6pftZOAZYEHLLJe2WczrgG9QHeeyCHi2ZfpF9XIAsl7+BlS7hHamev2vz8xFTXUertuwUd3mh1vm+RTdPQysC2wMPA94oqlO67zaimq31Q+BNzS1fRHV+3YwcGBEnJCZZ7Sp+w6qXWLj6lELqbYMf5Xq/VoykDZIPbI58KcBlrXPqNlnPMdlpo818EF1oHEC/14/P7l+Ph9Yp6Xsd+ppM9vM56R62nLgk8CkevzmVJ1aUgWu3dvUzfoxfYiv4dN1/YeojrMZX49/HrAN8FHgvS11jqzrzO5n3o22LaI6e+0V9fgAXtJU7pvACVRBee163DrAW6kOdE7gzDbz37OeNq+fdsyuyx05mGkDqPujetqNVAelT6jHbwh8DPhz/bnt2lLvJVT/JBL4OfDievx6wAep/rEtrKef3Ov13IePzAR4cdN3+mfAlEHUtc9I+4zn+sMte2r4MvD/gC2ptvJ9rb8KEbE+1Zcc4PTM/OfGtMy8NyLeSfWrbzfgn4HdR7jNu9TDL2XmpU3Lfga4jWqL5XA9BLwp6107WfVSK7YKZOZ7Wytk5p+Bn0TEzVSd95ER8eHMfHIE2jNsEfEG6ktJAHtl5oqtFVkdtH1aRCyjev8+BuzbVP3jVFsI/gD8bWY+Xdd7EvhyRKxHdfC7tNrIzDsi4gKqPRBvBO6JiKuAXwNzgaszc0BbuAbAPsM+Y7XjvnIBkJlPAJ+vn368DnL92Ydq18LSprrN81xG35f4dRExdSTa2uSJerjZCM+32RnZdAzPYGTmHcDNVL9gXzWSjRqmxvUVZzZ32i1m1cO9ImIcQEQE8PZ6/FcanXaLr1Kd8COtbt5L9aN2KbA2sDfwCeAS4KGI+E1EHFav58Nhn2Gfsdox7KnZ16mOUXk+cPwAyjdOvrgxMx/rUOZK+o5ZWeVkjWH6z3p4fERcGBFvqo8rGUnX9FcgIv46Iv4tIv4UEU82HzhNdeFWgBeOcLuG46/q4T9ExAPtHsB1dZn1qC48C9WusA3rv+e0m3FmLqa6SLe0WsnMpZn5j/Ttvfg3qj0AjZPFdqI6ZOV7wzxpwD7DPmO1425crZCZT0bEZ6l24X44Is7KzG4nS0yph/d2mefTEfEo1UkTUzqVG4rMvCAidqW6UOi76sfyiPhf4CfA2Zl5/zAX03XXTkT8C9Xu74bGQdzP1M8bB0QPZEvpWGlsCd2gfvRnvXrY/Pnd16V8x/VB6rXMfAj41/pBRLyA6ni5T1MFwYOA/2EAh7J0YJ9hn7HaccueWv0rcDfVWV//OMA664xec7rLzPcDrwBOpTq4+M9Uuz8+BdwWEX89zEUs6zQhIt5E1WkvozrBZWuqk1s2ycypmTkVuLZRfJjtGEmN7/3fZmYM4DFvkPNfnV6r1FVmPpiZ51LteXiwHv2eYczSPsM+Y7Vj2NNK6gOFG8fZnRgRm3Yp3vgFO61TgYhYl75N+iN1APRKMvPmzPxMZu5FtcvgrcBNVL+Mvx0RzxuN5VJtAQA4NzNPycw/1QdjN3tBa6XVQOMf2su7llpV8+fXbRfTaB5DKY2KzHwE+HH99KWjtBj7jPbsM0aZYU/tnE919tgk+q5R185v6+E2EbF5hzK703e4wG9bpjU6uRH7VVcfl/NT+jrVzaguw9LQuB3cSCxzi3p4Q7uJETGN6pd7OyPZjsFqHFN0wCDr3UHfPYzbnlldn9gzo9006Tmgcb235uvn2WfYZzznGfa0isx8lmoXA8AxdP5F9guqM2KfB3y4dWJ9Rtan6qdXZdPtiWqNs2k3HEo7I2LtLpObz4Zr3s08rGW2aBzPuH2H6Z+lc8fcaMdAjn8ZaTPr4YyI6HovzIjYqPF3vQXiB/XTEyOi3e774+k7XkdaLUTEVhHxkn7KrEd1eRGA3zVNss+wz3jOM+ypk1nALVRXj9+rXYHMXELVOUF1RuwnImIiQL2l79+orrHXuOByq5vr4eGNU/UH6fKI+JeI2D2abjEUEdvR1zndT7VLt3WZL+/v9j4DcFk9fH9EvKcRPiPiRRHxbeCdQKezlG+jOiB7g4gY7K/lYcnMn1FdCR/gvIg4JSJW7EaJiI0i4m8j4sdUl6podhrwNNX9Qy+JiK3qOhMi4kSqQwAGdAcUaQxtB/whIn4YEQe3rO/rR8RbgauArerRzSdn2GfYZzz3jcaVmn2s/g9a7qDRocwB9F0VvtMdNMYB324q8yzVmWXL6bt7xjEd5n9UU72ngLuoLtr5xQG+ht811V9WL/eppnFLgL3b1JvTVObRepnzgF2ayvR7dw+qa3Vd0/LaH2t6/im6X5G++X1b2NSOA5vKdKvfcdoA6q5P3xXxm9vweMu489vUfUf9WhtlHqP6J9T4Fd94XSf3ej334SMzobqQcrY8nqTvzg3N3+GPt6lvn2Gf8Zx+uGVP3fyQVY+zW0lmLsvMI4ADqXbrLgQmUm1R+zfgNZl5Voe651Nd6PQ3VB3BllQne3Q7KaTZ3wOfobpB+HyqrZAAtwJnUN2u6Io29d4OnAXcWbd1Wv1Yd4DLbbR/KdV9Ij9HdWzK8vp1XAa8NTP7uyr8B6h+9f6Baldzox0TB9OOocjMJZm5P9WV7n9IdemDCVT/jG6n2rJ7INVu/Na6/w7sClxK9XmvTbUV+ESqYyVbDziXeiozfw68DPgQ1UWUb68nTaRah39LdXHfv8zMz646B/sM+4zntqhTtyRJkgrklj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYOP7L7Jm2nTTTXP69Om9boak1dj111//SGZO6b/k8NgfSepPt/7IsNfB9OnTue6663rdDEmrsYi4ayyWY38kqT/d+iN340rqmTe96U18+9vfHvGykjQUpfZJXmevgxkzZqS/pKVVTZzYd/3WJ598knXWWYdx46q73f3rv/4rhx12WK+aNuYi4vrMHPWbuNsfSZ3ZJ1W69UfuxpU0KIsXL17x9/Tp0zn33HN5wxvesEq5Z599lvHj7WIkjS77pP65G1fSiJg9ezZbbLEFp59+OlOnTuWoo47iscceY99992XKlClstNFG7Lvvvtxzzz0r6uy5556ce+65AMycOZPddtuND33oQ2y00UZstdVW/Nd//deQyt55553svvvuTJo0iTe84Q0ce+yxvOtd7xqjd0LS6sA+qY9hT9KIeeCBB1iwYAF33XUX55xzDsuXL+eoo47irrvuYv78+UyYMIHjjjuuY/1rr72Wl73sZTzyyCN85CMf4e/+7u/odKhJt7KHHnoor3nNa3j00Uc5+eSTufDCC0fl9UpavdknVQx7kkbMWmutxSmnnMI666zDhAkT2GSTTTjggANYb731mDRpEp/4xCeYM2dOx/rTpk3jve99L+PGjeOII47g/vvv58EHHxxU2fnz5zN37lxOPfVU1l57bXbbbTf222+/0XrJklZj9kkVw56kETNlyhTWXXfdFc+ffPJJ3v/+9zNt2jQmT57M7rvvzsKFC1m2bFnb+lOnTl3x93rrrQesfDzOQMred999bLzxxivGAWy55ZZDf1GSnrPskyqGPUkjJiJWev6lL32JP/zhD1x77bU88cQTXHnllQAdd4OMhM0224wFCxbw5JNPrhh39913j9ryJK2+7JMqhj1Jo2bRokVMmDCBDTfckAULFnDKKaeM+jKnTZvGjBkzOPnkk1m6dCnXXHMNP/nJT0Z9uZJWf2tqn2TYkzRqTjzxRJ566ik23XRTdtllF/7mb/5mTJb73e9+l2uuuYZNNtmET37ykxxyyCGss846Y7JsSauvNbVP8qLKHXgRU6kchxxyCNtuu+2I/4r3osqShmI0+iQvqjxGdvzwBb1uggbo+i8c3usmaBTNnTuXjTfemK222opf/OIX/PjHP+akk07qdbMkraF63ScZ9iQV54EHHuDtb387jz76KFtssQVnn302O+ywQ6+bJWkN1es+ybAnqThvfetbeetb39rrZkgS0Ps+yRM0JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKtiAw15EvCwiToiI70TErRGxPCIyIg7sUmdmXabT49Z+lnloRFwVEY9HxOKIuC4ijo2Iru0eaj1JkqTSjB9E2aOBE4a4nP8Bbm8z/v5OFSLiTOAY4GngCuAZYG/gDGDviDgoM5eNVD1JkqQSDSbs/R74AnAdcD3wLWCPAdY9NzNnDnRBEXEAVWB7ANg9M2+rx78A+CWwP3Ac8LWRqCdJklSqAe/WzMxzM/MjmXlRZv5pNBsFfKwefrQR2Oo2PEi1hRHgpDa7ZYdaT5IkqUirXeiJiC2AHYGlwMWt0zNzDnAvMBXYZbj1JEmSSjaY3bjDsVdEvBKYCDwI/Aq4LDOXtym7Qz28OTOf6jC/ucDmddmrh1lPkiSpWGMV9g5vM+6WiHhHZt7UMn6renhXl/nNbyk7nHqSJEnFGu3duL8Djge2o9qq90JgX+BG4OXA5RGxeUudifVwSZf5Lq6Hk0ag3goR8b76Mi3XPfzww11mI0mjy/5I0kgZ1bCXmV/NzK9n5i2ZuSQz78/MS4HXAL8Gnk/fSRUN0ag+yMUNtV5ze8/JzBmZOWPKlClDnY0kDZv9kaSR0pMTNDJzKXBa/fTNLZMX1cOJdNaYtqhp3FDrSZIkFauXZ+M27p7Ruht3Xj2c1qXuli1lh1NPkiSpWL0Me5vUw8Ut42+oh9tFxIQOdXdqKTucepIkScXqZdg7uB7ObR6ZmXcDvwXWBg5qrRQRewBbUN0l45rh1pMkSSrZqIW9iHhVROwbEeNaxo+PiA9SnaUL8JU21RvH850eEVs31X0+cFb99HNtrtM31HqSJElFGvB19iLi1fQFJqgunQLw2Yj4UGNkZjbuTjEd+BGwICL+CNxDdcmT7akuwbKc6rZmP29dVmZ+PyLOprrF2U0RcTnwDLA3MBm4BDhjpOpJkiSVajAXVZ4M7Nxm/DYdyt8IfI3qMivTqO5akVSh73zgzMy8vtPCMvOYiPgVcCywBzCO6qSO84CzO22dG2o9SZKkEg047GXmbPquZTeQ8ncCJw6+SSvNYxYwa6zqSZIklaaXJ2hIkiRplBn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKNr7XDZAkjY4dP3xBr5ugAbj+C4f3ugkqnFv2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSrYgMNeRLwsIk6IiO9ExK0RsTwiMiIOHEDdQyPiqoh4PCIWR8R1EXFsRHRd/ljXkyRJKs34QZQ9GjhhsAuIiDOBY4CngSuAZ4C9gTOAvSPioMxc1ut6kiRJJRrMlq7fA18ADgG2Bub0VyEiDqAKXg8Ar8zMfTNzf2Ab4P+A/YHjel1PkiSpVAMOe5l5bmZ+JDMvysw/DbDax+rhRzPztqZ5PUi1pRDgpDa7V8e6niRJUpFGLfRExBbAjsBS4OLW6Zk5B7gXmArs0qt6kiRJJRvNLVw71MObM/OpDmXmtpTtRT1JkqRijWbY26oe3tWlzPyWsr2oJ0mSVKzRDHsT6+GSLmUW18NJPay3QkS8r75My3UPP/xwl9lI0uiyP5I0UkYz7EU9zNW83gqZeU5mzsjMGVOmTBnqbCRp2OyPJI2U0Qx7i+rhxC5lGtMWNY0b63qSJEnFGsxFlQdrXj2c1qXMli1le1FPGjXzT92+103QAL3o0zf1ugmSNCpGc8veDfVwu4iY0KHMTi1le1FPkiSpWKMW9jLzbuC3wNrAQa3TI2IPYAuqu11c06t6kiRJJRvtO0mcVg9Pj4itGyMj4vnAWfXTz2Xm8h7XkyRJKtKAj9mLiFfTF5gAXl4PPxsRH2qMzMxdmv7+fkScTXWrspsi4nLgGWBvYDJwCXBG67LGup4kSVKpBnOCxmRg5zbjt+lWKTOPiYhfAccCewDjgFuB84CzO21lG+t6kiRJJRpw2MvM2fRdy25QMnMWMGt1rydJklSa0T5mT5IkST1k2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCje91AyRJ0tiYf+r2vW6CBuBFn75pROfnlj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlgox72ImJmRGSXx61d6h4aEVdFxOMRsTgirouIYyOia7uHWk+SJKk048dwWf8D3N5m/P3tCkfEmcAxwNPAFcAzwN7AGcDeEXFQZi4bqXqSJEklGsuwd25mzhxIwYg4gCqwPQDsnpm31eNfAPwS2B84DvjaSNSTJEkq1eq6W/Nj9fCjjcAGkJkPAkfXT09qs1t2qPUkSZKKtNqFnojYAtgRWApc3Do9M+cA9wJTgV2GW0+SJKlkY7kbd6+IeCUwEXgQ+BVwWWYubym3Qz28OTOf6jCvucDmddmrh1lPkiSpWGMZ9g5vM+6WiHhHZt7UNG6renhXl3nNbyk7nHqSJEnFGovduL8Djge2o9qq90JgX+BG4OXA5RGxeVP5ifVwSZd5Lq6Hk0ag3goR8b76Mi3XPfzww11mI0mjy/5I0kgZ9bCXmV/NzK9n5i2ZuSQz78/MS4HXAL8Gnk/fiRUA0ag6yEUNtV5zW8/JzBmZOWPKlClDnY0kDZv9kaSR0rMTNDJzKXBa/fTNTZMW1cOJdNaYtqhp3FDrSZIkFavXZ+M27p7RvBt3Xj2c1qXeli1lh1NPkiSpWL0Oe5vUw8VN426oh9tFxIQO9XZqKTucepIkScXqddg7uB7ObYzIzLuB3wJrAwe1VoiIPYAtqO6Scc1w60mSJJVsVMNeRLwqIvaNiHEt48dHxAepztIF+EpL1caxfKdHxNZN9Z4PnFU//Vyba/QNtZ4kSVKRRvs6e9OBHwELIuKPwD1Ulz3ZnuoSLMupbm328+ZKmfn9iDib6hZnN0XE5cAzwN7AZOAS4IzWhQ21niRJUqlGO+zdCHyN6jIr06juXJFUoe984MzMvL5dxcw8JiJ+BRwL7AGMozqh4zzg7E5b54ZaT5IkqUSjGvYy807gxGHUnwXMGqt6kiRJpen1CRqSJEkaRYY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWDFhr2IODQiroqIxyNicURcFxHHRkSxr1mSJKlVkcEnIs4EvgvMAK4CLgNeCpwBfD8ixvWweZIkSWOmuLAXEQcAxwAPAK/MzH0zc39gG+D/gP2B43rYREmSpDFTXNgDPlYPP5qZtzVGZuaDwNH105PcnStJktYERQWeiNgC2BFYClzcOj0z5wD3AlOBXca2dZIkSWOvqLAH7FAPb87MpzqUmdtSVpIkqVilhb2t6uFdXcrMbykrSZJUrPG9bsAIm1gPl3Qps7geTmqdEBHvA97XKBcRfxjBtj2XbQo80utGjKT44hG9bsJzWXHrAwCfiaHUmjbSzWiwP+qouPXP/mjYilsnRro/Ki3sNd6dHErlzDwHOGfkmlOGiLguM2f0uh1aPbg+jA37o/Zc/9TKdaJ/pe3GXVQPJ3Yp05i2qEsZSZKkIpQW9ubVw267VrZsKStJklSs0sLeDfVwu4iY0KHMTi1l1T93JamZ64N6yfVPrVwn+hGZQzq8bbUVEdcDrwaOyMwLWqbtAcymurvG5pm5fOxbKEmSNHZK27IHcFo9PD0itm6MjIjnA2fVTz9n0JMkSWuC4rbsAUTEWVS3RnsauBx4BtgbmAxcAhyYmct61kBJkqQxUmTYA4iIQ4Fjge2BccCtwHnA2W7VkyRJa4piw54kSZLKPGZPkiRJNcOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFWx8rxuwutp0001z+vTpvW6GpNXY9ddf/0hmThnt5dgfSepPt/7IsNfB9OnTue6663rdDEmrsYi4ayyWY38kqT/d+iN340oaUxHB7bffDsAHPvAB/umf/mlAZQfru9/9Lvvss8+Q6kpaM6wp/ZFhT9KgvfGNb+TTn/70KuN//OMfM3XqVJ599tkBzecb3/gGn/rUp4bdnnnz5hERKy33sMMO4xe/+MWw5y1p9WZ/1D/DnqRBO/LII7nwwgvJzJXGX3jhhRx22GGMH+8RIpLGhv1R/wx7kgbtbW97GwsWLOCqq65aMe6xxx7jpz/9Kfvttx+vfe1r2XDDDdlss8047rjjWLp0adv5HHnkkXzyk59c8fwLX/gCm222GS984Qs577zzVip76aWXssMOOzB58mS23HJLTj755BXTdt99dwA23HBDJk6cyDXXXMPMmTPZbbfdVpS5+uqr2Wmnndhggw3YaaeduPrqq1dM23PPPfnUpz7FrrvuyqRJk9hnn3145JFHhvUeSRob9kf9M+xJGrQJEyZw8MEHc8EFF6wYd9FFF7HtttsyceJEvvKVr/DII49wzTXXcMUVV3DWWWf1O8+f/exnfPGLX+Syyy7jtttu4/LLL19p+vrrr88FF1zAwoULufTSSzn77LO55JJLALjyyisBWLhwIYsXL+a1r33tSnUXLFjAW97yFo4//ngeffRRPvjBD/KWt7yFRx99dEWZWbNmcf755/PQQw+xdOlSvvjFLw717ZE0huyP+mfYkzQkRxxxBBdffDFPPfUUABdccAFHHHEEO+64I7vssgvjx49n+vTpvP/972fOnDn9zu+iiy7iqKOO4hWveAXrr7/+Sr+Uofq1u/3227PWWmvxyle+kne+850Dmi9Uv8K32WYb3v3udzN+/Hje+c53su222/KTn/xkRZmjjjqKl770pSv+cfzud78b8Hshqbfsj7oz7Ekakt12240pU6bw4x//mDvuuIO5c+dy6KGH8sc//pF9992XqVOnMnnyZD7+8Y8PaBfEfffdx5Zbbrni+bRp01aafu2117LXXnsxZcoUNthgA77xjW8MeNfGfffdt8r8pk2bxr333rvi+dSpU1f8vd5667F48eIBzVtS79kfdWfYkzRkhx9+OBdccAEXXngh++yzDy94wQs4+uij2Xbbbbntttt44okn+OxnP7vKgdPtbLbZZtx9990rns+fP3+l6Yceeij77bcfd999N48//jgf+MAHVsw3IrrO+4UvfCF33bXyJajmz5/P5ptvPtCXKmk1Z3/UmWFP0pAdfvjhXH755Xzzm9/kiCOOAGDRokVMnjyZiRMncuutt3L22WcPaF4HH3wwM2fO5JZbbuHJJ5/klFNOWWn6okWL2HjjjVl33XX5zW9+w6xZs1ZMmzJlCmuttRZ33HFH23m/+c1v5o9//COzZs3i2Wef5Xvf+x633HIL++677xBfuaTVjf1RZ4Y9SUM2ffp0/uqv/oolS5aw3377AfDFL36RWbNmMWnSJN773vdyyCGHDGheb3rTmzjxxBN5/etfz9Zbb83rX//6laafddZZfPrTn2bSpEmceuqpHHzwwSumrbfeenziE59g1113ZcMNN+TXv/71SnU32WQTfvrTn/KlL32JTTbZhM9//vP89Kc/ZdNNNx3mOyBpdWF/1FkMZHPmmmjGjBnp7YkkdRMR12fmjNFejv2RpP5064/csidJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUsAGHvYh4WUScEBHfiYhbI2J5RGREHNilzsy6TKfHrf0s89CIuCoiHo+IxRFxXUQcGxFd2z3UepIkSaUZP4iyRwMnDHE5/wPc3mb8/Z0qRMSZwDHA08AVwDPA3sAZwN4RcVBmLhupepIkSSUaTNj7PfAF4DrgeuBbwB4DrHtuZs4c6IIi4gCqwPYAsHtm3laPfwHwS2B/4DjgayNRT5IkqVQD3q2Zmedm5kcy86LM/NNoNgr4WD38aCOw1W14kGoLI8BJbXbLDrWeJElSkVa70BMRWwA7AkuBi1unZ+Yc4F5gKrDLcOtJkiSVbDC7cYdjr4h4JTAReBD4FXBZZi5vU3aHenhzZj7VYX5zgc3rslcPs54kSVKxxirsHd5m3C0R8Y7MvKll/Fb18K4u85vfUnY49SRJkoo12rtxfwccD2xHtVXvhcC+wI3Ay4HLI2LzljoT6+GSLvNdXA8njUC9FSLiffVlWq57+OGHu8xGkkaX/ZGkkTKqYS8zv5qZX8/MWzJzSWben5mXAq8Bfg08n76TKhqiUX2Qixtqveb2npOZMzJzxpQpU4Y6G0kaNvsjSSOlJydoZOZS4LT66ZtbJi+qhxPprDFtUdO4odaTJEkqVi/Pxm3cPaN1N+68ejitS90tW8oOp54kSVKxehn2NqmHi1vG31APt4uICR3q7tRSdjj1JEmSitXLsHdwPZzbPDIz7wZ+C6wNHNRaKSL2ALagukvGNcOtJ0mSVLJRC3sR8aqI2DcixrWMHx8RH6Q6SxfgK22qN47nOz0itm6q+3zgrPrp59pcp2+o9SRJkoo04OvsRcSr6QtMUF06BeCzEfGhxsjMbNydYjrwI2BBRPwRuIfqkifbU12CZTnVbc1+3rqszPx+RJxNdYuzmyLicuAZYG9gMnAJcMZI1ZMkSSrVYC6qPBnYuc34bTqUvxH4GtVlVqZR3bUiqULf+cCZmXl9p4Vl5jER8SvgWGAPYBzVSR3nAWd32jo31HqSJEklGnDYy8zZ9F3LbiDl7wROHHyTVprHLGDWWNWTJEkqTS9P0JAkSdIoM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBVsfK8bIEkaHTt++IJeN0EDcP0XDu91E1Q4t+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVLABh72IeFlEnBAR34mIWyNieURkRBw4gLqHRsRVEfF4RCyOiOsi4tiI6Lr8sa4nSZJUmvGDKHs0cMJgFxARZwLHAE8DVwDPAHsDZwB7R8RBmbms1/UkSZJKNJgtXb8HvgAcAmwNzOmvQkQcQBW8HgBemZn7Zub+wDbA/wH7A8f1up4kSVKpBhz2MvPczPxIZl6UmX8aYLWP1cOPZuZtTfN6kGpLIcBJbXavjnU9SZKkIo1a6ImILYAdgaXAxa3TM3MOcC8wFdilV/UkSZJKNppbuHaohzdn5lMdysxtKduLepIkScUazbC3VT28q0uZ+S1le1FPkiSpWKMZ9ibWwyVdyiyuh5N6WG+FiHhffZmW6x5++OEus5Gk0WV/JGmkjGbYi3qYq3m9FTLznMyckZkzpkyZMtTZSNKw2R9JGimjGfYW1cOJXco0pi1qGjfW9SRJkoo1mmFvXj2c1qXMli1le1FPkiSpWKMZ9m6oh9tFxIQOZXZqKduLepIkScUatbCXmXcDvwXWBg5qnR4RewBbUN3t4ppe1ZMkSSrZaN9J4rR6eHpEbN0YGRHPB86qn34uM5f3uJ4kSVKRxg+0YES8mr7ABPDyevjZiPhQY2Rm7tL09/cj4myqW5XdFBGXA88AewOTgUuAM1qXNdb1JEmSSjXgsEcVlnZuM36bbpUy85iI+BVwLLAHMA64FTgPOLvTVraxridJklSiAYe9zJxN37XsBiUzZwGzVvd6kiRJpRntY/YkSZLUQ4Y9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKtj4XjdAkiSNjfmnbt/rJmgAXvTpm0Z0fm7ZkyRJKphhT5IkqWCGPUmSpIJ5zN4I2vHDF/S6CRqg679weK+bIEnSmHDLniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUsFEPexExMyKyy+PWLnUPjYirIuLxiFgcEddFxLER0bXdQ60nSZJUmrG8g8b/ALe3GX9/u8IRcSZwDPA0cAXwDLA3cAawd0QclJnLRqqeJElSicYy7J2bmTMHUjAiDqAKbA8Au2fmbfX4FwC/BPYHjgO+NhL1JEmSSrW67tb8WD38aCOwAWTmg8DR9dOT2uyWHWo9SZKkIq12oScitgB2BJYCF7dOz8w5wL3AVGCX4daTJEkq2Vjuxt0rIl4JTAQeBH4FXJaZy1vK7VAPb87MpzrMay6weV326mHWkyRJKtZYhr3D24y7JSLekZk3NY3bqh7e1WVe81vKDqeeJElSscZiN+7vgOOB7ai26r0Q2Be4EXg5cHlEbN5UfmI9XNJlnovr4aQRqLdCRLyvvkzLdQ8//HCX2UjS6LI/kjRSRj3sZeZXM/PrmXlLZi7JzPsz81LgNcCvgefTd2IFQDSqDnJRQ63X3NZzMnNGZs6YMmXKUGcjScNmfyRppPTsBI3MXAqcVj99c9OkRfVwIp01pi1qGjfUepIkScXq9dm4jbtnNO/GnVcPp3Wpt2VL2eHUkyRJKlavw94m9XBx07gb6uF2ETGhQ72dWsoOp54kSVKxeh32Dq6HcxsjMvNu4LfA2sBBrRUiYg9gC6q7ZFwz3HqSJEklG9WwFxGvioh9I2Jcy/jxEfFBqrN0Ab7SUrVxLN/pEbF1U73nA2fVTz/X5hp9Q60nSZJUpNG+zt504EfAgoj4I3AP1WVPtqe6BMtyqlub/by5UmZ+PyLOprrF2U0RcTnwDLA3MBm4BDijdWFDrSdJklSq0Q57NwJfo7rMyjSqO1ckVeg7HzgzM69vVzEzj4mIXwHHAnsA46hO6DgPOLvT1rmh1pMkSSrRqIa9zLwTOHEY9WcBs8aqniRJUml6fYKGJEmSRpFhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSqYYU+SJKlghj1JkqSCGfYkSZIKZtiTJEkqmGFPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWCGPUmSpIIZ9iRJkgpm2JMkSSrY+F43QCrZ/FO373UTNEAv+vRNvW6CJI0Kt+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBDHuSJEkFM+xJkiQVzLAnSZJUMMOeJElSwQx7kiRJBTPsSZIkFcywJ0mSVDDDniRJUsEMe5IkSQUz7EmSJBXMsCdJklQww54kSVLBig17EXFoRFwVEY9HxOKIuC4ijo2IYl+zJElSqyKDT0ScCXwXmAFcBVwGvBQ4A/h+RIzrYfMkSZLGTHFhLyIOAI4BHgBemZn7Zub+wDbA/wH7A8f1sImSJEljpriwB3ysHn40M29rjMzMB4Gj66cnuTtXkiStCYoKPBGxBbAjsBS4uHV6Zs4B7gWmAruMbeskSZLGXlFhD9ihHt6cmU91KDO3pawkSVKxSgt7W9XDu7qUmd9SVpIkqVjje92AETaxHi7pUmZxPZzUOiEi3ge8r1EuIv4wgm17LtsUeKTXjRhJ8cUjet2E57Li1gcAPhNDqTVtpJvRYH/UUXHrn/3RsBW3Tox0f1Ra2Gu8OzmUypl5DnDOyDWnDBFxXWbO6HU7tHpwfRgb9kftuf6pletE/0rbjbuoHk7sUqYxbVGXMpIkSUUoLezNq4fddq1s2VJWkiSpWKWFvRvq4XYRMaFDmZ1ayqp/7kpSM9cH9ZLrn1q5TvQjMod0eNtqKyKuB14NHJGZF7RM2wOYTXV3jc0zc/nYt1CSJGnslLZlD+C0enh6RGzdGBkRzwfOqp9+zqAnSZLWBMVt2QOIiLOobo32NHA58AywNzAZuAQ4MDOX9ayBkiRJY6TELXtk5jHAYcBvgT2ANwK3A8cBB6wpQS8i5kVENj2WR8QTEXF3RFweEadFxCt73U6NrojYuv78n42IyR3KnNS0nuzcocy+9fT7R7fFKpV9kuyPeqPIsAeQmbMyc9fMnJyZ62fmjpl55hq6+/bnwLeBC4CfUQXfGcBJwI0R8R8RMbWH7dMoyszbgXuAccDrOhTbs8Pf7cr8cqhtafpnP32o81AR7JPWUPZHvVFs2NNKPpeZR9aPgzNzL2Bj4G3AHcBbgTkRsUkvG6lRNace7tk6ISLGA7sCN1Md8rBKmdoe9XD2yDZNayD7pDWb/dEYM+ytoTJzeWb+mOrX9O3AS4Ev9bZVGkWz6+GebabtSHWx8Z8Bc4Fd6w53hXp3yw710yH/kpY6sU9ao8yuh3u2mWZ/NAoMe2u4zHwMOLF++q7GrpOImF5v3p4XEeMj4kMRcWNELImIhfX0ZRGxoNM1DSPieRFxfz2fl4/RS1J7jQ5xhzbHyTR+Ic8BrqS6b/SOLWV2o9rtcl9m3tYYGRE7R8S/R8Q9EbE0Ih6ud8Ht1lw5Io6MiKTvgud3thy7NX34L1ElsE9aI9gfjTHDngD+E1hA9eXZq2VaAD8A/j/gIeA/gJszcx7wE2Aj4J0d5nsAMBWYnZm3jHyzNVCZ+SfgbtofJ7MnsBz4FZ13rzSer/gVHRH/CFwDHEx17cofU22ReQvVLrj3NtW/neoYrSX18x/UzxuPxUN5XSqWfVLB7I96IDN9FPqguiVcAnsOoOzlddl/rp9Pr58ncBewdZs6e9fTr+swzyvr6Qf2+r3wkQAX1p/H55vGjQMeB26on08CngX+s6XutXXdv6+f/039/F5g55ayu9bzXAq8tGVaY52c3uv3w8fYP+yTfDR9FvZHY/hwy54aHq6H7Q6I/lhWZ1CtJDOvAG4BdoyI1zRPi4hXUP1iu4/q2obqvdn1cM+mca+muv7klQCZuYjqVoK7RcQ4gIiYWJeDvl/Sp9TDv8/Ma5sXkpn/A/wT8Dzg/SP6CrQmsU8q2+x6uGfTOPujUWLYU0NjXWh3aZofdal3Rj08pmX8sfXwnMx8djgN04hpdIyvbjpOZs96OKepXOtxMq8DxgP3ZOafImJTqntMPwH8osOyGvN77Qi0W2sm+6Sy2R+NIcOeGjathwtaxj+UmU91qXcB1SbyQyJiY4CImAS8i2rzuzeoXk1k5h3AfKpdJY0Dlveg2o1xZVPR1uNkGgdMNzrnraiOm5oMPNtyYHPWBz7/pi47ZcRfiNYU9kkFsz8aW+P7L6LSRUTQdxr7TS2Tu3WqZOaSiDgP+AfgPcAXgSOoTp2/ODO9uvnqZQ7wbmDPiPg5VSd7S2Y+0lTmKurjqoDP09fJzq6H4+rh4/S/O+yRfqZLq7BPWmPYH40Rw56gOltpI6pfvbOHUP9M4ATgAxHxZar7EjfGa/Uym7pzBV4FbADMai6QmY9FxE1Ux8lsQN/uk8Yv6bvr4TOZeeToNldrKPukNcNs7I/GhLtx13ARsRHwlfrptzPzocHOI6vT6P8LeAnwWeDlVJdCmNO1onphxXEyVHcpgJWPj2loHCdzPNWPwvmZeSdAZt5LtbVl04jYc5DLX1oP/aGptuyT1ij2R2PEsLeGioi1ImI/qiuUbw3cCnx4GLP8ej38aD08axjz0iipO8i7qHZ9HF+PvrJN0UaH+w/1sPUq9Z+qh9+JiH1aK0fE2hGxX0S0HhB9bz38i0E1XMWzT1rz2B+NnaivNaMCRcQ8qiuE/5zqIpMA61IdpPpqYMN63CXA+5t/QddXEL8TuCszpw9gWQH8H/AyYBGweX3avFYzEfFt4PD66R8z82VtyjwfeLBp1FGZObOlzAepjqEZB/wR+APVL+UtqdaDDYCjM/MbTXWOB75GtY78AlhYT/poZj463Nem1Zt9klrZH42N4jddCoA31sOkumL4QuB6qjOUZmXm74e7gMzMiLic6kt1gZ3qau2X9HWu7X5Fk5kPRcStwLb1qNltynw5Iq4A/h/VMTd/TXWM1f1Uv8R/AvywpdoZVGfNHQbsC6xTj/9noKjOVV3ZJ6nB/mgMuGVPIyIi1qY6jf4FwHbprYgk9ZB9ktTHY/Y0Uo6l6lR/ZqcqaTVgnyTV3LKnIYuIl1EdQP1Cqt0yy4CdMvPGnjZM0hrJPklqz2P2NBybAX8H/Bm4EfiknaqkHrJPktpwy54kSVLBPGZPkiSpYIY9SZKkghn2JEmSCmbYkyRJKphhT5IkqWD/P3UfZ5d+1YMaAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<Figure size 720x720 with 4 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
"source": [
"# plot distribution of wet days in calibration period\n",
"fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(10, 10))\n",
--
GitLab