Refactor.

1d0f412e · Frisinghelli Daniel · 432286c1 · 1d0f412e · 1d0f412e · 1d0f412e
Commit 1d0f412e authored 3 years ago by Frisinghelli Daniel
--- 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 = 'dem'\n",
-    "DEM_FEATURES = '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",
+    "# align datasets\n",
-    "if LOSS == 'BernoulliGammaLoss':\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_true, y_refe, y_pred_pr = xr.align(y_true.precipitation, y_refe.precipitation, y_pred.precipitation, join='override')"
-    "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)"
+  },
+  {
+   "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",

 %% Cell type:markdown id:4735431f-6741-437e-b0fc-dd6d8eaa22ca tags:
 # Evaluate ERA-5 downscaling: precipitation
 %% Cell type:markdown id:a87da113-4b0f-4ac8-9721-19c85848acec tags:
 We used **1981-1991 as training** period and **1991-2010 as reference** period. The results shown in this notebook are based on the model predictions on the reference period.
 %% Cell type:markdown id:ad1769d4-9c0c-4e3f-9adf-ef02bb43c047 tags:
 **Predictors on pressure levels (500, 850)**:
 - Geopotential (z)
 - Temperature (t)
 - Zonal wind (u)
 - Meridional wind (v)
 - Specific humidity (q)
 **Predictors on surface**:
 - Mean sea level pressure (msl)
 **Auxiliary predictors**:
 - Elevation from Copernicus EU-DEM v1.1 (dem)
 - Day of the year (doy)
 %% Cell type:markdown id:f9334da7-17d1-45ef-9ed9-5c2bee9fcdcc tags:
 Define the predictand and the model to evaluate:
 %% Cell type:code id:a81acde7-16a2-4087-bc08-95b084adbd06 tags:
 ``` python
 # define the model parameters
 PREDICTAND = 'pr'
 MODEL = 'USegNet'
 PPREDICTORS = 'ztuvq'
 # PPREDICTORS = ''
 PLEVELS = ['500', '850']
 # PLEVELS = []
 SPREDICTORS = 'p'
-DEM = ''
+DEM = 'dem'
-DEM_FEATURES = 'dem'
+DEM_FEATURES = ''
 DOY = ''
-WET_DAY_THRESHOLD = '2'
+WET_DAY_THRESHOLD = '1'
 # LOSS = 'MSELoss'
 LOSS = 'BernoulliGammaLoss'
+SEASON = 'season'
 ```
 %% Cell type:markdown id:dd188df0-69ee-44b0-82b2-d212994dc271 tags:
 ### Imports
 %% Cell type:code id:06792bf2-b33b-4728-ba60-d60fab46779d tags:
 ``` python
 # builtins
 import datetime
 import warnings
 import calendar
 # externals
 import xarray as xr
 import numpy as np
 import matplotlib.pyplot as plt
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 import scipy.stats as stats
 from IPython.display import Image
 from sklearn.metrics import r2_score, roc_curve, auc, classification_report
 # locals
 from climax.main.io import ERA5_PATH, OBS_PATH, TARGET_PATH, DEM_PATH
 from pysegcnn.core.utils import search_files
 from pysegcnn.core.graphics import plot_classification_report
 ```
 %% Cell type:code id:cd14134e-f9be-4935-877b-ef6d34e03d2e tags:
 ``` python
 # mapping from predictands to variable names
 NAMES = {'tasmin': 'minimum temperature', 'tasmax': 'maximum temperature', 'pr': 'precipitation'}
 ```
 %% Cell type:markdown id:a3d92474-2bc1-4035-8938-c5fbb07ae891 tags:
 ### Model architecture
 %% Cell type:code id:c68e022b-41c2-438b-bfb0-e27eddee89bf tags:
 ``` python
 Image("./Figures/architecture.png", width=900, height=400)
 ```
 %% Cell type:markdown id:c8833efe-c715-4872-9aee-a0b5766f5c67 tags:
 ### Loss function
 %% Cell type:markdown id:bb801367-6872-4a0b-bff5-70cb6746e057 tags:
 For precipitation, the network is optimizing the negative log-likelihood of a Bernoulli-Gamma distribution after [Cannon (2008)](http://journals.ametsoc.org/doi/10.1175/2008JHM960.1).
 %% Cell type:markdown id:a8775d5e-5ad4-47e4-8fef-6dc230e15dee tags:
 Bernoulli-Gamma distribution:
 %% Cell type:markdown id:ab10f8de-d8d2-4427-b9c8-5d68803543c3 tags:
 $$P(y \mid, p, \alpha, \beta) = \begin{cases} 1 - p, & \text{for } y = 0\\ p \cdot \frac{y^{\alpha -1} \exp(-y/\beta)}{\beta^{\alpha} \tau(\alpha)}, & \text{for } y > 0\end{cases}$$
 %% Cell type:markdown id:6b6dbd06-1f0e-4c52-84f2-b7ff31c75726 tags:
 Log-likelihood function:
 %% Cell type:markdown id:e41c7b39-f352-4a98-820f-9a7345b3283c tags:
 $$\mathcal{J}(p, \alpha, \beta \mid y) = \underbrace{(1 - P(y > 0)) \log(1 - p)}_{\text{Bernoulli}} + \underbrace{P(y > 0) \cdot \left(\log(p) + (\alpha - 1) \log(y) - \frac{y}{\beta} - \alpha \log(\beta) - \log(\tau(\alpha))\right)}_{\text{Gamma}}$$
 %% Cell type:markdown id:5a0c55f0-79fb-4501-b3cf-b5414399a3d9 tags:
 ### Load datasets
 %% Cell type:code id:efa76f8e-c089-47ff-a001-d4c2a11c4d6d tags:
 ``` python
 # construct file pattern to match
 PATTERN = '_'.join([MODEL, PREDICTAND, PPREDICTORS, *PLEVELS])
 PATTERN = '_'.join([PATTERN, SPREDICTORS]) if SPREDICTORS else PATTERN
 PATTERN = '_'.join([PATTERN, DEM]) if DEM else PATTERN
 PATTERN = '_'.join([PATTERN, DEM_FEATURES]) if DEM_FEATURES else PATTERN
 PATTERN = '_'.join([PATTERN, DOY]) if DOY else PATTERN
 PATTERN = '_'.join([PATTERN, '{}mm'.format(str(WET_DAY_THRESHOLD).replace('.', ''))]) if WET_DAY_THRESHOLD else PATTERN
 PATTERN = '_'.join([PATTERN, LOSS])
+PATTERN = '_'.join([PATTERN, SEASON]) if SEASON else PATTERN
 PATTERN
 ```
 %% Cell type:code id:ecaba394-f802-481f-8274-c44e2f5fdf1a tags:
 ``` python
 # digital elevation model
 dem = xr.open_dataset(search_files(DEM_PATH, 'eu_dem_v11_stt.nc').pop())
 dem = dem.Band1.to_dataset().rename({'Band1': 'elevation'})
 ```
 %% Cell type:code id:a5db133d-2c36-4e84-879e-20e617e821f1 tags:
 ``` python
 # model predictions and observations NetCDF
 y_pred = xr.open_dataset(search_files(TARGET_PATH.joinpath(PREDICTAND), '.'.join([PATTERN, 'nc$'])).pop())
 y_true = xr.open_dataset(search_files(OBS_PATH.joinpath(PREDICTAND), 'OBS_pr(.*).nc$').pop())
 try:
    y_pred = y_pred.rename({'pr': 'precipitation'})
 except ValueError:
    pass
 ```
 %% Cell type:code id:325f3086-85f0-4c28-b37d-7370d2d92405 tags:
 ``` python
 # reference dataset: ERA-5 precipitation
 y_refe = xr.open_dataset(search_files(ERA5_PATH.joinpath('ERA5', 'total_precipitation'), '.nc$').pop())
 ```
 %% Cell type:code id:528c3116-7707-45ca-b811-0adad7bc20f3 tags:
 ``` python
 # subset to time period covered by predictions
 y_true = y_true.sel(time=y_pred.time)
 y_refe = y_refe.sel(time=y_pred.time).drop_vars('lambert_azimuthal_equal_area')
 y_refe = y_refe.rename({'tp': 'precipitation'})
 ```
 %% Cell type:code id:70b903cb-e597-45d3-b575-0ebaf7a45649 tags:
 ``` python
-# align datasets and mask missing values in model predictions
+# align datasets
-if LOSS == 'BernoulliGammaLoss':
+if len(y_pred.data_vars) > 1:
    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')
-    y_pred_prob = y_pred_prob.where(~np.isnan(y_true), other=np.nan)
 else:
    y_true, y_refe, y_pred_pr = xr.align(y_true.precipitation, y_refe.precipitation, y_pred.precipitation, join='override')
-y_pred_pr = y_pred_pr.where(~np.isnan(y_true), other=np.nan)
+```
-y_refe = y_refe.where(~np.isnan(y_true), other=np.nan)
+%% Cell type:code id:aa3f8e67-6b49-46d4-a956-a0cee7b3923a tags:
+``` python
+# mask missing values
+mask = ~np.isnan(y_true)
+y_pred_pr = y_pred_pr.where(mask, other=np.nan)
+y_refe = y_refe.where(mask, other=np.nan)
+if len(y_pred.data_vars) > 1:
+    y_pred_prob = y_pred_prob.where(mask, other=np.nan)
 ```
 %% Cell type:code id:4e71997f-7808-463a-8b45-dcac639ebe88 tags:
 ``` python
 # align digital elevation model
 _, dem = xr.align(y_true.isel(time=0), dem, join='override')
 dem = dem.where(~np.isnan(y_true.isel(time=0)), other=np.nan)
 ```
 %% Cell type:markdown id:b269a131-cf5b-4c6c-9f8e-a5408250aa83 tags:
 ## Model validation: precipitation amount
 %% Cell type:markdown id:0fa1fe82-0d6e-4676-b5b8-9eac2fd28ffb tags:
 ### Coefficient of determination: monthly mean
 %% Cell type:code id:4d11cee6-1ebd-4424-8ec6-92e5a196bac4 tags:
 ``` python
 # calculate monthly mean precipitation (mm / month)
 y_pred_values = y_pred_pr.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim='time').values
 y_true_values = y_true.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim='time').values
 ```
 %% Cell type:code id:6e55cdd5-8adb-40c5-9742-46f4fc3d4be9 tags:
 ``` python
 # apply mask of valid pixels
 mask = (~np.isnan(y_pred_values) & ~np.isnan(y_true_values))
 y_pred_values = y_pred_values[mask]
 y_true_values = y_true_values[mask]
 ```
 %% Cell type:code id:13a9ff21-34ea-4db1-9c0a-bcb7ea1001f7 tags:
 ``` python
 # calculate coefficient of determination
 r2 = r2_score(y_true_values, y_pred_values)
 r2
 ```
 %% Cell type:code id:a1b33431-9f2b-4e42-bbb0-f4fa258edd98 tags:
 ``` python
 # group timeseries by month and calculate mean over time and space
 y_pred_ac = y_pred_pr.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim=('y', 'x', 'time'), skipna=True).values.squeeze()
 y_true_ac = y_true.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim=('y', 'x', 'time'), skipna=True).values.squeeze()
 ```
 %% Cell type:code id:c683858e-8b7f-4c76-a40f-6f68397b3479 tags:
 ``` python
 # scatter plot of observations vs. predictions
 fig, ax = plt.subplots(1, 1, figsize=(10, 10))
 # plot only a subset of data: otherwise plot is overloaded ...
 # subset = np.random.choice(np.arange(0, len(y_pred_values)), size=int(1e3), replace=False)
 # ax.plot(y_true_values[subset], y_pred_values[subset], 'o', alpha=.5, markeredgecolor='grey', markerfacecolor='none', markersize=3);
 # plot entire dataset
 ax.plot(y_true_values, y_pred_values, 'o', alpha=.5, markeredgecolor='grey', markerfacecolor='none', markersize=3);
 # plot 1:1 mapping line
 interval = np.arange(0, 300, 50)
 ax.plot(interval, interval, color='k', lw=2, ls='--')
 # add coefficient of determination: calculated on entire dataset!
 ax.text(interval[-1] - 2, interval[0] + 2, s='Coefficient of determination R$^2$ = {:.2f}'.format(r2), ha='right', fontsize=18)
 # format axes
 ax.set_ylim(interval[0], interval[-1])
 ax.set_xlim(interval[0], interval[-1])
 ax.set_xticks(interval)
 ax.set_xticklabels(interval, fontsize=16)
 ax.set_yticks(interval)
 ax.set_yticklabels(interval, fontsize=16)
 ax.set_xlabel('Observed', fontsize=18)
 ax.set_ylabel('Predicted', fontsize=18)
 ax.set_title('Monthly mean {} (mm / month)'.format(NAMES[PREDICTAND]), fontsize=20, pad=10);
 # add axis for annual cycle
 axins = inset_axes(ax, width="30%", height="40%", loc=2, borderpad=0.25)
 axins.plot(y_pred_ac, ls='--', color='k', label='Predicted')
 axins.plot(y_true_ac, ls='-', color='k', label='Observed')
 axins.legend(frameon=False, fontsize=12, loc='lower center');
 axins.set_yticks(np.arange(0, 200, 50))
 axins.set_yticklabels(np.arange(0, 200, 50), fontsize=12)
 axins.yaxis.tick_right()
 axins.set_xticks(np.arange(0, 12))
 axins.set_xticklabels([calendar.month_name[i + 1] for i in np.arange(0, 12)], rotation=90, fontsize=12)
 # save figure
 fig.savefig('../Notebooks/Figures/{}_r2.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')
 ```
 %% Cell type:markdown id:3538a109-4d3a-418e-9d4d-c3c120256ab2 tags:
 ### Bias
 %% Cell type:markdown id:43c922bd-9e9b-4812-b096-f3bde06fb249 tags:
 Calculate yearly average bias over entire reference period:
 %% Cell type:code id:52efc445-054e-4184-9667-4e761bc12a84 tags:
 ``` python
 # yearly average bias over reference period
 y_pred_yearly_avg = y_pred_pr.groupby('time.year').mean(dim='time')
 y_refe_yearly_avg = y_refe.groupby('time.year').mean(dim='time')
 y_true_yearly_avg = y_true.groupby('time.year').mean(dim='time')
 bias_yearly_avg = ((y_pred_yearly_avg - y_true_yearly_avg) / y_true_yearly_avg) * 100
 bias_yearly_avg_ref = ((y_refe_yearly_avg - y_true_yearly_avg) / y_true_yearly_avg) * 100
 print('(Model) Yearly average relative bias: {:.2f}%'.format(bias_yearly_avg.mean().item()))
 print('(ERA-5) Yearly average relative bias: {:.2f}%'.format(bias_yearly_avg_ref.mean().item()))
 ```
 %% Cell type:code id:92a094cb-2a88-4f4a-8c87-1025672d6fe7 tags:
 ``` python
 # mean absolute error over reference period
 mae_avg = np.abs(y_pred_yearly_avg - y_true_yearly_avg)
 mae_avg_ref = np.abs(y_refe_yearly_avg - y_true_yearly_avg)
 print('(Model) Yearly average MAE: {:.2f} mm'.format(mae_avg.mean().item()))
 print('(ERA-5) Yearly average MAE: {:.2f} mm'.format(mae_avg_ref.mean().item()))
 ```
 %% Cell type:code id:25d397e2-1b46-4f48-b39d-8632a7e56288 tags:
 ``` python
 # root mean squared error over reference period
 rmse_avg = ((y_pred_yearly_avg - y_true_yearly_avg) ** 2).mean()
 rmse_avg_ref = ((y_refe_yearly_avg - y_true_yearly_avg) **2).mean()
 print('(Model) Yearly average RMSE: {:.2f} mm / day'.format(rmse_avg.item()))
 print('(ERA-5) Yearly average RMSE: {:.2f} mm / day'.format(rmse_avg_ref.item()))
 ```
 %% Cell type:code id:d6bbdcd6-3920-4856-b90a-d56f0e5ab2df tags:
 ``` python
 # Pearson's correlation coefficient over reference period
 for year in y_pred_yearly_avg.year:
    y_p = y_pred_yearly_avg.sel(year=year).values
    y_t = y_true_yearly_avg.sel(year=year).values
    r, _ = stats.pearsonr(y_p[~np.isnan(y_p)], y_t[~np.isnan(y_t)])
    print('({:0d}) Pearson correlation: {:.2f}'.format(year.item(), np.asarray(r).mean()))
 ```
 %% Cell type:code id:ca683581-6b3b-4abd-ad46-157ebded19c6 tags:
 ``` python
 # plot yearly average MAE of reference vs. prediction
 vmin, vmax = 0, 5
 fig, axes = plt.subplots(1, 3, figsize=(24, 8), sharex=True, sharey=True)
 # plot bias of ERA-5 reference
 reference = bias_yearly_avg_ref.mean(dim='year')
 im1 = axes[0].imshow(reference.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
 axes[0].text(x=reference.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(reference.mean().item()), fontsize=14, ha='right')
 # plot MAE of model
 prediction = bias_yearly_avg.mean(dim='year')
 im2 = axes[1].imshow(prediction.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
 axes[1].text(x=reference.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(prediction.mean().item()), fontsize=14, ha='right')
 # plot topography
 im_dem = axes[2].imshow(dem['elevation'].values, origin='lower', cmap='terrain', vmin=0, vmax=4000)
 # set titles
 axes[0].set_title('ERA-5', fontsize=14, pad=10);
 axes[1].set_title('DCEDN', fontsize=14, pad=10);
 axes[2].set_title('Copernicus EU-DEM v1.1', fontsize=14, pad=10)
 # adjust axes
 for ax in axes.flat:
    ax.axes.get_xaxis().set_ticklabels([])
    ax.axes.get_xaxis().set_ticks([])
    ax.axes.get_yaxis().set_ticklabels([])
    ax.axes.get_yaxis().set_ticks([])
    ax.axes.axis('tight')
    ax.set_xlabel('')
    ax.set_ylabel('')
    ax.set_axis_off()
 # adjust figure
 # fig.suptitle('Average yearly mean absolute error: 1991 - 2010', fontsize=20);
 fig.subplots_adjust(hspace=0, wspace=0, top=0.85)
-# add colorbar for bias
+# add colorbar for dem
 axes = axes.flatten()
 cbar_ax_bias = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,
                             0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])
 cbar_bias = fig.colorbar(im_dem, cax=cbar_ax_bias)
 cbar_bias.set_label(label='Elevation (m)', fontsize=14)
 cbar_bias.ax.tick_params(labelsize=14, pad=10)
 # add colorbar for predictand
 cbar_ax_predictand = fig.add_axes([axes[0].get_position().x0, axes[0].get_position().y0 - 0.1,
                                   axes[-1].get_position().x0 - axes[0].get_position().x0,
                                   0.03])
 cbar_predictand = fig.colorbar(im1, cax=cbar_ax_predictand, orientation='horizontal')
 cbar_predictand.set_label(label='Relative mean error (%)', fontsize=14)
 cbar_predictand.ax.tick_params(labelsize=14, pad=10)
 # add metrics: MAE and RMSE
 #axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_avg[NAMES[PREDICTAND]].item()) + 'mm day$^{-1}$', fontsize=14, ha='right')
 #axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_avg[NAMES[PREDICTAND]].item()) + 'mm$^2$ day$^{-2}$', fontsize=14, ha='right')
 # save figure
 fig.savefig('../Notebooks/Figures/{}_rbias_ERA_vs_model.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')
 ```
 %% Cell type:code id:f022fcbc-077f-4f76-8330-47fa72a1fcec tags:
 ``` python
 # plot average of observation, prediction, and bias
 vmin, vmax = 0, 5
 fig, axes = plt.subplots(1, 3, figsize=(24, 8), sharex=True, sharey=True)
 axes = axes.flatten()
 for ds, ax in zip([y_true_yearly_avg, y_pred_yearly_avg, bias_yearly_avg], axes):
    if ds is bias_yearly_avg:
        ds = ds.mean(dim='year')
        im2 = ax.imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
        ax.text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
    else:
        im1 = ax.imshow(ds.mean(dim='year').values, origin='lower', cmap='BuPu', vmin=vmin, vmax=vmax)
 # set titles
 axes[0].set_title('Observed', fontsize=16, pad=10);
 axes[1].set_title('Predicted', fontsize=16, pad=10);
 axes[2].set_title('Bias', fontsize=16, pad=10);
 # adjust axes
 for ax in axes.flat:
    ax.axes.get_xaxis().set_ticklabels([])
    ax.axes.get_xaxis().set_ticks([])
    ax.axes.get_yaxis().set_ticklabels([])
    ax.axes.get_yaxis().set_ticks([])
    ax.axes.axis('tight')
    ax.set_xlabel('')
    ax.set_ylabel('')
 # adjust figure
 fig.suptitle('Average {}: 1991 - 2010'.format(NAMES[PREDICTAND]), fontsize=20);
 fig.subplots_adjust(hspace=0, wspace=0, top=0.85)
 # add colorbar for bias
 axes = axes.flatten()
 cbar_ax_bias = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,
                             0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])
 cbar_bias = fig.colorbar(im2, cax=cbar_ax_bias)
 cbar_bias.set_label(label='Relative bias / (%)', fontsize=16)
 cbar_bias.ax.tick_params(labelsize=14)
 # add colorbar for predictand
 cbar_ax_predictand = fig.add_axes([axes[0].get_position().x0, axes[0].get_position().y0 - 0.1,
                                   axes[-1].get_position().x0 - axes[0].get_position().x0,
                                   0.05])
 cbar_predictand = fig.colorbar(im1, cax=cbar_ax_predictand, orientation='horizontal')
 cbar_predictand.set_label(label='{} / '.format(NAMES[PREDICTAND].capitalize()) + '(mm day$^{-1}$)', fontsize=16)
 cbar_predictand.ax.tick_params(labelsize=14)
 # add metrics: MAE and RMSE
 axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_avg.mean().item()) + 'mm day$^{-1}$', fontsize=14, ha='right')
 axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_avg.mean().item()) + 'mm day$^{-1}$', fontsize=14, ha='right')
 # save figure
 fig.savefig('../Notebooks/Figures/{}_average_bias.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')
 ```
 %% Cell type:markdown id:f9eccfb4-c8f6-41f1-9636-9a46221cf53e tags:
 ### Seasonal bias
 %% Cell type:markdown id:fe07f855-a68a-4bdf-a923-ef334fb82ad3 tags:
 Calculate seasonal bias:
 %% Cell type:code id:65813f9f-094d-43c0-8261-3364294631f2 tags:
 ``` python
 # group data by season: (DJF, MAM, JJA, SON)
 y_true_snl = y_true.groupby('time.season').mean(dim='time')
 y_pred_snl = y_pred_pr.groupby('time.season').mean(dim='time')
 y_refe_snl = y_refe.groupby('time.season').mean(dim='time')
 bias_snl = ((y_pred_snl - y_true_snl) / y_true_snl) * 100
 bias_snl_ref = ((y_refe_snl - y_true_snl) / y_true_snl) * 100
 ```
 %% Cell type:code id:2a70e6ac-54fb-47b8-ae76-a1b51d9d1f17 tags:
 ``` python
 # print average bias per season: ERA-5
 for season in bias_snl_ref.season:
    print('(ERA-5) Average bias for season {}: {:.1f}%'.format(season.values.item(), bias_snl_ref.sel(season=season).mean().item()))
 ```
 %% Cell type:code id:40fa07e4-2df1-4ac9-b27f-cd902e16ed30 tags:
 ``` python
 # print average bias per season: model
 for season in bias_snl.season:
    print('(Model) Average bias for season {}: {:.1f}%'.format(season.values.item(), bias_snl.sel(season=season).mean().item()))
 ```
 %% Cell type:markdown id:a73b3ed4-eb44-4f3e-b240-7f6fe44bedd3 tags:
 Plot seasonal differences, taken from the [xarray documentation](xarray.pydata.org/en/stable/examples/monthly-means.html).
 %% Cell type:code id:2b2390e4-8bf1-41bd-b7d1-33e92fe8bc65 tags:
 ``` python
 # plot seasonal differences
 seasons = ('DJF', 'JJA')
 fig, axes = plt.subplots(nrows=1, ncols=len(seasons) + 1, figsize=(24,8), sharex=True, sharey=True)
 axes = axes.flatten()
 # plot annual average bias
 ds = bias_yearly_avg.mean(dim='year')
 axes[0].imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
 axes[0].set_title('Annual', fontsize=16);
 axes[0].text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
 # plot seasonal average bias
 for ax, season in zip(axes[1:], seasons):
    ds = bias_snl.sel(season=season)
    ax.imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
    ax.set_title(season, fontsize=16);
    ax.text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
 # adjust axes
 for ax in axes.flat:
    ax.axes.get_xaxis().set_ticklabels([])
    ax.axes.get_xaxis().set_ticks([])
    ax.axes.get_yaxis().set_ticklabels([])
    ax.axes.get_yaxis().set_ticks([])
    ax.axes.axis('tight')
    ax.set_xlabel('')
    ax.set_ylabel('')
 # adjust figure
 fig.suptitle('Average bias of {}: 1991 - 2010'.format(NAMES[PREDICTAND]), fontsize=20);
 fig.subplots_adjust(hspace=0, wspace=0, top=0.85)
 # add colorbar for bias
 axes = axes.flatten()
 cbar_ax = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,
                        0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])
 cbar = fig.colorbar(im2, cax=cbar_ax)
 cbar.set_label(label='Relative bias / (%)', fontsize=16)
 cbar.ax.tick_params(labelsize=14)
 # save figure
 fig.savefig('../Notebooks/Figures/{}_average_bias_seasonal.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')
 ```
 %% Cell type:markdown id:ed9dee0e-9319-4ad5-a619-4ea30b116917 tags:
 ### Bias of extreme values
 %% Cell type:code id:92f0ce1d-2cd3-46f6-be1b-86cfe5c8e9d7 tags:
 ``` python
 # extreme quantile of interest
 quantile = 0.98
 ```
 %% Cell type:code id:bf98b3e0-7886-485f-a3bc-4bd64b7b2814 tags:
 ``` python
 # calculate extreme quantile for each year
 with warnings.catch_warnings():
    warnings.simplefilter('ignore', category=RuntimeWarning)
    y_pred_ex = y_pred_pr.groupby('time.year').quantile(quantile, dim='time')
    y_true_ex = y_true.groupby('time.year').quantile(quantile, dim='time')
    y_refe_ex = y_refe.groupby('time.year').quantile(quantile, dim='time')
 ```
 %% Cell type:code id:89893df3-987a-4794-a6e9-0893c925218c tags:
 ``` python
 # calculate bias in extreme quantile for each year
 bias_ex = ((y_pred_ex - y_true_ex) / y_true_ex) * 100
 bias_ex_ref = ((y_refe_ex - y_true_ex) / y_true_ex) * 100
 ```
 %% Cell type:code id:9a610a70-bb94-47e6-b470-e01bf4a295c0 tags:
 ``` python
 # bias of extreme quantile: ERA-5
 print('(ERA-5) Yearly average bias for P{:.0f}: {:.1f}%'.format(quantile * 100, bias_ex_ref.mean().item()))
 ```
 %% Cell type:code id:423962d7-453d-4c86-b93e-a4571bfec1c4 tags:
 ``` python
 # bias of extreme quantile: Model
 print('(Model) Yearly average bias for P{:.0f}: {:.1f}%'.format(quantile * 100, bias_ex.mean().item()))
 ```
 %% Cell type:code id:180b723a-49b0-4692-8249-b457eda48d44 tags:
 ``` python
 # mean absolute error in extreme quantile
 mae_ex = np.abs(y_pred_ex - y_true_ex).mean()
 mae_ex_ref = np.abs(y_refe_ex - y_true_ex).mean()
 ```
 %% Cell type:code id:7331acdd-60bb-414e-a0e7-e81956c3a1bf tags:
 ``` python
 # mae of extreme quantile: ERA-5
 print('(ERA-5) Yearly average MAE for P{:.0f}: {:.1f} mm / day'.format(quantile * 100, mae_ex_ref.item()))
 ```
 %% Cell type:code id:d7bf1a6c-e111-4c72-8774-64a3939cbe50 tags:
 ``` python
 # mae of extreme quantile: Model
 print('(Model) Yearly average MAE for P{:.0f}: {:.1f} mm / day'.format(quantile * 100, mae_ex.item()))
 ```
 %% Cell type:code id:e6309621-7d16-43b7-a9d4-c2857dbe8270 tags:
 ``` python
 # root mean squared error in extreme quantile
 rmse_ex = ((y_pred_ex - y_true_ex) ** 2).mean()
 rmse_ex_ref = ((y_refe_ex - y_true_ex) ** 2).mean()
 ```
 %% Cell type:code id:712a2109-14b9-45a3-9194-e2b917c5ba3f tags:
 ``` python
 # rmse of extreme quantile: ERA-5
 print('(ERA-5) Yearly average RMSE for P{:.0f}: {:.1f} mm / day'.format(quantile * 100, rmse_ex_ref.item()))
 ```
 %% Cell type:code id:1dfad687-c16f-4767-8ddd-96c3e56bd96a tags:
 ``` python
 # rmse of extreme quantile: Model
 print('(Model) Yearly average RMSE for P{:.0f}: {:.1f} mm / day'.format(quantile * 100, rmse_ex.item()))
 ```
 %% Cell type:code id:95338821-a78e-4a1c-ae73-25031ba7f6ac tags:
 ``` python
 # plot extremes of observation, prediction, and bias
 vmin, vmax = 10, 40
 fig, axes = plt.subplots(1, 3, figsize=(24, 6), sharex=True, sharey=True)
 axes = axes.reshape(1, -1)
 for ds, ax in zip([y_true_ex, y_pred_ex, bias_ex], axes[i, ...]):
    if ds is bias_ex:
        ds = ds.mean(dim='year')
        im2 = ax.imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
        ax.text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
    else:
        im1 = ax.imshow(ds.mean(dim='year').values, origin='lower', cmap='BuPu', vmin=vmin, vmax=vmax)
 # set titles
 axes[0, 0].set_title('Observed', fontsize=16, pad=10);
 axes[0, 1].set_title('Predicted', fontsize=16, pad=10);
 axes[0, 2].set_title('Bias', fontsize=16, pad=10);
 # adjust axes
 for ax in axes.flat:
    ax.axes.get_xaxis().set_ticklabels([])
    ax.axes.get_xaxis().set_ticks([])
    ax.axes.get_yaxis().set_ticklabels([])
    ax.axes.get_yaxis().set_ticks([])
    ax.axes.axis('tight')
    ax.set_xlabel('')
    ax.set_ylabel('')
 # adjust figure
 fig.suptitle('Average P{:.0f} of {}: 1991 - 2010'.format(quantile * 100, NAMES[PREDICTAND]), fontsize=20);
 fig.subplots_adjust(hspace=0, wspace=0, top=0.85)
 # add colorbar for bias
 axes = axes.flatten()
 cbar_ax = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,
                        0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])
 cbar = fig.colorbar(im2, cax=cbar_ax)
 cbar.set_label(label='Relative bias / (%)', fontsize=16)
 cbar.ax.tick_params(labelsize=14)
 # add colorbar for predictand
 cbar_ax_predictand = fig.add_axes([axes[0].get_position().x0, axes[0].get_position().y0 - 0.1,
                                   axes[-1].get_position().x0 - axes[0].get_position().x0,
                                   0.05])
 cbar_predictand = fig.colorbar(im1, cax=cbar_ax_predictand, orientation='horizontal')
 cbar_predictand.set_label(label='{} / '.format(NAMES[PREDICTAND].capitalize()) + '(mm day$^{-1}$)', fontsize=16)
 cbar_predictand.ax.tick_params(labelsize=14)
 # add metrics: MAE and RMSE
 axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_ex.item())  + 'mm day$^{-1}$', fontsize=14, ha='right')
 axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_ex.item())  + 'mm$^2$ day$^{-2}$', fontsize=14, ha='right')
 # save figure
 fig.savefig('../Notebooks/Figures/{}_average_bias_p{:.0f}.png'.format(PREDICTAND, quantile * 100), dpi=300, bbox_inches='tight')
 ```
 %% Cell type:markdown id:f7758b71-e844-4a47-b741-25cc2d277814 tags:
 ### Bias of extremes: winter vs. summer
 %% Cell type:code id:de16b327-dc4b-4645-b3a0-abc8f72c4225 tags:
 ``` python
 # group data by season and compute extreme percentile
 with warnings.catch_warnings():
    warnings.simplefilter('ignore', category=RuntimeWarning)
    y_true_ex_snl = y_true.groupby('time.season').quantile(quantile, dim='time')
    y_pred_ex_snl = y_pred_pr.groupby('time.season').quantile(quantile, dim='time')
    y_refe_ex_snl = y_refe.groupby('time.season').quantile(quantile, dim='time')
 ```
 %% Cell type:code id:185a375c-b7bc-42cf-a57d-08268e824f21 tags:
 ``` python
 # compute relative bias in seasonal extremes
 bias_ex_snl = ((y_pred_ex_snl - y_true_ex_snl) / y_true_ex_snl) * 100
 bias_ex_snl_ref = ((y_refe_ex_snl - y_true_ex_snl) / y_true_ex_snl) * 100
 ```
 %% Cell type:code id:e5f82b16-7ebf-4b73-80dd-0bc7e6a305a1 tags:
 ``` python
 # print average bias in extreme per season: ERA-5
 for season in bias_ex_snl_ref.season:
    print('(ERA-5) Average bias of P{:.0f} for season {}: {:.1f}%'.format(quantile * 100, season.values.item(), bias_ex_snl_ref.sel(season=season).mean().item()))
 ```
 %% Cell type:code id:0671b4f3-52c6-4167-b35b-03902dbe11a3 tags:
 ``` python
 # print average bias in extreme per season: Model
 for season in bias_ex_snl.season:
    print('(Model) Average bias of P{:.0f} for season {}: {:.1f}%'.format(quantile * 100, season.values.item(), bias_ex_snl.sel(season=season).mean().item()))
 ```
 %% Cell type:code id:322ee9c4-92db-4041-8a36-45aa8d2021b5 tags:
 ``` python
 # plot seasonal differences
 seasons = ('DJF', 'JJA')
 fig, axes = plt.subplots(nrows=1, ncols=len(seasons) + 1, figsize=(24,8), sharex=True, sharey=True)
 axes = axes.flatten()
 # plot annual average bias of extreme
 ds = bias_ex.mean(dim='year')
 axes[0].imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
 axes[0].set_title('Annual', fontsize=16);
 axes[0].text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
 # plot seasonal average bias of extreme
 for ax, season in zip(axes[1:], seasons):
    ds = bias_ex_snl.sel(season=season)
    ax.imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
    ax.set_title(season, fontsize=16);
    ax.text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
 # adjust axes
 for ax in axes.flat:
    ax.axes.get_xaxis().set_ticklabels([])
    ax.axes.get_xaxis().set_ticks([])
    ax.axes.get_yaxis().set_ticklabels([])
    ax.axes.get_yaxis().set_ticks([])
    ax.axes.axis('tight')
    ax.set_xlabel('')
    ax.set_ylabel('')
 # adjust figure
 fig.suptitle('Average bias of P{:.0f} of {}: 1991 - 2010'.format(quantile * 100, NAMES[PREDICTAND]), fontsize=20);
 fig.subplots_adjust(hspace=0, wspace=0, top=0.85)
 # add colorbar for bias
 axes = axes.flatten()
 cbar_ax = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,
                        0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])
 cbar = fig.colorbar(im2, cax=cbar_ax)
 cbar.set_label(label='Relative bias / (%)', fontsize=16)
 cbar.ax.tick_params(labelsize=14)
 # save figure
 fig.savefig('../Notebooks/Figures/{}_average_bias_seasonal_ex.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')
 ```
 %% Cell type:markdown id:1f1a76c5-0bcb-4f79-832c-caf958e5703a tags:
 ### Frequency of wet days
 %% Cell type:code id:4491395e-f69a-4749-ac73-c4cbce60c3bc tags:
 ``` python
 # minimum precipitation (mm / day) defining a wet day
 WET_DAY_THRESHOLD = float(WET_DAY_THRESHOLD)
 ```
 %% Cell type:code id:40676a79-26a2-4714-900e-d530fda32f8f tags:
 ``` python
 # true and predicted frequency of wet days
 mask = (~np.isnan(y_true)) & (~np.isnan(y_pred_pr))
 wet_days_true = (y_true >= WET_DAY_THRESHOLD).where(mask, other=np.nan).astype(np.float32)
 wet_days_pred = (y_pred_pr >= WET_DAY_THRESHOLD).where(mask, other=np.nan).astype(np.float32)
 ```
 %% Cell type:code id:1ad9ad08-ad1f-4a76-9020-c09a442385c9 tags:
 ``` python
 # number of wet days in reference period: annual
 n_wet_days_true = wet_days_true.sum(dim='time', skipna=False)
 n_wet_days_pred = wet_days_pred.sum(dim='time', skipna=False)
 ```
 %% Cell type:code id:54088572-1873-4949-9937-0ea099e9c2b6 tags:
 ``` python
 # frequency of wet days in reference period: annual
 f_wet_days_true = (n_wet_days_true / len(wet_days_true.time)) * 100
 f_wet_days_pred = (n_wet_days_pred / len(wet_days_pred.time)) * 100
 ```
 %% Cell type:code id:3534d1d9-f50d-40f3-b1f0-d34c4152a37c tags:
 ``` python
 # frequency of wet days in reference period: seasonal
 f_wet_days_true_snl = wet_days_true.groupby('time.season').mean(dim='time', skipna=False)
 f_wet_days_pred_snl = wet_days_pred.groupby('time.season').mean(dim='time', skipna=False)
 ```
 %% Cell type:code id:60ac0f93-f9ca-4f75-8786-e692efe3a556 tags:
 ``` python
 # relative bias of frequency of wet vs. dry days: annual
 bias_wet = ((f_wet_days_pred - f_wet_days_true) / f_wet_days_true) * 100
 # relative bias of frequency of wet vs. dry days: seasonal
 bias_wet_snl = ((f_wet_days_pred_snl - f_wet_days_true_snl) / f_wet_days_true_snl) * 100
 ```
 %% Cell type:code id:e8e9d30a-2fbc-4f29-bc82-bfbc155430a0 tags:
 ``` python
 # plot average of observation, prediction, and bias
 fig, axes = plt.subplots(2, 3, figsize=(24, 12), sharex=True, sharey=True)
 axes = axes.flatten()
 # plot annual average bias of extreme
 ds = bias_wet
 im = axes[0].imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
 axes[0].set_title('Annual', fontsize=16);
 axes[0].text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
 # plot seasonal average bias of extreme
 for ax, season in zip(axes[1:], bias_wet_snl.season):
    ds = bias_wet_snl.sel(season=season)
    ax.imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
    ax.set_title(season.item(), fontsize=16);
    ax.text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
 # adjust axes
 for ax in axes.flat:
    ax.axes.get_xaxis().set_ticklabels([])
    ax.axes.get_xaxis().set_ticks([])
    ax.axes.get_yaxis().set_ticklabels([])
    ax.axes.get_yaxis().set_ticks([])
    ax.axes.axis('tight')
    ax.set_xlabel('')
    ax.set_ylabel('')
 # turn off last axis
 axes[-1].set_visible(False)
 # adjust figure
 fig.suptitle('Frequency of wet days (>= {:.1f} mm): 1991 - 2010'.format(WET_DAY_THRESHOLD), fontsize=20);
 fig.subplots_adjust(hspace=0.1, wspace=0, top=0.925)
 # add colorbar
 cbar_ax_predictand = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,
                                   0.01, axes[0].get_position().y1 - axes[-1].get_position().y0])
 cbar_predictand = fig.colorbar(im, cax=cbar_ax_predictand)
 cbar_predictand.set_label(label='Relative bias / (%)', fontsize=16)
 cbar_predictand.ax.tick_params(labelsize=14)
 # save figure
 fig.savefig('../Notebooks/Figures/{}_bias_wet_days.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')
 ```
 %% Cell type:markdown id:fde67674-8bed-45b2-b4ca-4a507084b5a8 tags:
 ### Mean wet day precipitation
 %% Cell type:code id:aa2d2b5a-e192-4c34-975c-d4e61c969a15 tags:
 ``` python
 # calculate mean wet day precipitation
 dii_true = (y_true * wet_days_true).sum(dim='time', skipna=False) / n_wet_days_true
 dii_pred = (y_pred_pr * wet_days_pred).sum(dim='time', skipna=False) / n_wet_days_pred
 ```
 %% Cell type:code id:d6331bb6-afe8-4e62-9fa8-60bce9af42ec tags:
 ``` python
 # calculate relative bias of mean wet day precipitation
 bias_dii = ((dii_pred - dii_true) / dii_true) * 100
 ```
 %% Cell type:code id:36e3cbac-233d-4892-9756-941d91981534 tags:
 ``` python
 # plot average of observation, prediction, and bias
 fig, axes = plt.subplots(1, 3, figsize=(24, 6), sharex=True, sharey=True)
 for ds, ax in zip([dii_true, dii_pred, bias_dii], axes):
    if ds is bias_dii:
        im2 = ax.imshow(ds.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)
        ax.text(x=ds.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(ds.mean().item()), fontsize=14, ha='right')
    else:
        im1 = ax.imshow(ds.values, origin='lower', cmap='BuPu', vmin=0, vmax=15)
 # set titles
 axes[0].set_title('Observed', fontsize=16, pad=10);
 axes[1].set_title('Predicted', fontsize=16, pad=10);
 axes[2].set_title('Bias', fontsize=16, pad=10);
 # adjust axes
 for ax in axes.flat:
    ax.axes.get_xaxis().set_ticklabels([])
    ax.axes.get_xaxis().set_ticks([])
    ax.axes.get_yaxis().set_ticklabels([])
    ax.axes.get_yaxis().set_ticks([])
    ax.axes.axis('tight')
    ax.set_xlabel('')
    ax.set_ylabel('')
 # adjust figure
 fig.suptitle('Mean wet day (>= {:.1f} mm) precipitation: 1991 - 2010'.format(WET_DAY_THRESHOLD), fontsize=20);
 fig.subplots_adjust(hspace=0, wspace=0, top=0.85)
 # add colorbar for bias
 axes = axes.flatten()
 cbar_ax_bias = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,
                             0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])
 cbar_bias = fig.colorbar(im2, cax=cbar_ax_bias)
 cbar_bias.set_label(label='Relative bias / (%)', fontsize=16)
 cbar_bias.ax.tick_params(labelsize=14)
 # add colorbar for predictand
 cbar_ax_predictand = fig.add_axes([axes[0].get_position().x0, axes[0].get_position().y0 - 0.1,
                                   axes[-1].get_position().x0 - axes[0].get_position().x0,
                                   0.05])
 cbar_predictand = fig.colorbar(im1, cax=cbar_ax_predictand, orientation='horizontal')
 cbar_predictand.set_label(label='Mean wet day precipitation / (mm day$^{-1}$)', fontsize=16)
 cbar_predictand.ax.tick_params(labelsize=14)
 # save figure
 fig.savefig('../Notebooks/Figures/{}_bias_wet_days_p.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')
 ```
 %% Cell type:markdown id:359e6d22-5f7f-4203-bb2a-fa51e731d867 tags:
 ## Model validation: precipitation probability
 %% Cell type:markdown id:5bd5c551-d426-45d7-9338-5b5003ca676e tags:
 ### ROC: Receiver operating characteristics
 %% Cell type:code id:1bff2faf-c894-4917-8dcf-310106c153ad tags:
 ``` python
 # true and predicted probability of precipitation
 p_true = (y_true > float(WET_DAY_THRESHOLD)).values.flatten()
 p_pred = y_pred_prob.values.flatten()
 ```
 %% Cell type:code id:6c4a7097-26af-4497-bf52-43601038be41 tags:
 ``` python
 # apply mask of valid pixels
 mask = (~np.isnan(p_true) & ~np.isnan(p_pred))
 p_pred = p_pred[mask]
 p_true = p_true[mask].astype(float)
 ```
 %% Cell type:code id:e4c35894-e31f-4879-806e-0b3935412730 tags:
 ``` python
 # calculate ROC: false positive rate vs. true positive rate
 fpr, tpr, _ = roc_curve(p_true, p_pred)
 area = auc(fpr, tpr) # area under ROC curve
 rocss = 2 * area - 1 # ROC skill score (cf. https://journals.ametsoc.org/view/journals/clim/16/24/1520-0442_2003_016_4145_otrsop_2.0.co_2.xml)
 ```
 %% Cell type:code id:ac418718-44bf-4378-a1da-7fb75e45b3a8 tags:
 ``` python
 # plot ROC curve
 fig, ax = plt.subplots(1, 1, figsize=(10, 10))
 ax.plot(fpr, tpr, lw=2, label='Area={:.2f}, ROCSS={:.2f}'.format(area, rocss), color='k')
 # plot classifier with no skill
 interval = np.arange(-0.05, 1.1, 0.05)
 ax.plot([0, 1], [0, 1], lw=2, linestyle='--', color='k')
 ax.text(0.95, 0.975, 'Random Classifier', ha='right', va='top', rotation=45, fontsize=12)
 # plot perfect classifier
 ax.plot(0, 1, '-o', markersize=5, markerfacecolor='k', markeredgecolor='none')
 ax.text(0.02, 1, 'Perfect classifier', va='center', fontsize=12)
 # plot direction of increase / decrease
 ax.arrow(np.median(interval), np.median(interval), 0.1, -0.1, head_width=0.01, facecolor='k')
 ax.arrow(np.median(interval), np.median(interval), -0.1, 0.1, head_width=0.01, facecolor='k')
 ax.text(np.median(interval) + 0.05, np.median(interval) - 0.05, s='Worse', rotation=45, ha='left', fontsize=12)
 ax.text(np.median(interval) - 0.05, np.median(interval) + 0.05, s='Better', rotation=45, ha='left', fontsize=12)
 # adjust axes
 ax.set_xticks(np.arange(0, 1.1, 0.1))
 ax.set_xticklabels(['{:.2f}'.format(i) for i in np.arange(0, 1.1, 0.1)], fontsize=12)
 ax.set_yticks(np.arange(0, 1.1, 0.1))
 ax.set_yticklabels(['{:.2f}'.format(i) for i in np.arange(0, 1.1, 0.1)], fontsize=12)
 ax.set_xlim(interval[0], interval[-1])
 ax.set_ylim(interval[0], interval[-1])
 ax.set_xlabel('False Positive Rate', fontsize=14)
 ax.set_ylabel('True Positive Rate', fontsize=14)
 ax.set_title('ROC of precipitation probability: 1991 - 2010', fontsize=14, pad=10)
 ax.legend(frameon=False, loc='lower right', fontsize=14);
 # save figure
 fig.savefig('../Notebooks/Figures/{}_ROC.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')
 ```

--- 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_gamma, marker='*', color='k', label='Gamma')\n",
-    "ax.scatter(eq, tq_genpareto, color='k', label='GenPareto')\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_gamma, marker='*', color='k', label='Gamma')\n",
-    "    ax.scatter(eq, tq_gpare, color='k', label='GenPareto')\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",
+    "for dist, biasq in dists.items():\n",
-    "    biasq = bias_gamma if dist == 'gamma' else bias_genpa\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": {

 %% Cell type:markdown id:63805b4a-b30e-4c10-a948-bc59651ca7a6 tags:
 ### Imports
 %% Cell type:code id:28982ce9-bf0c-4eb1-8b9e-bec118359966 tags:
 ``` python
 # builtins
 import datetime
 import warnings
 import calendar
 # externals
 import xarray as xr
 import numpy as np
 import matplotlib.pyplot as plt
-import seaborn as sns
+import seequantilesn as sns
 import pandas as pd
 import scipy.stats as stats
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 import scipy.stats as stats
 from IPython.display import Image
 from sklearn.metrics import r2_score, roc_curve, auc, classification_report
 from sklearn.model_selection import train_test_split
 # locals
 from climax.main.io import ERA5_PATH, OBS_PATH, TARGET_PATH, DEM_PATH
 from climax.main.config import CALIB_PERIOD, VALID_PERIOD
 from pysegcnn.core.utils import search_files
 from pysegcnn.core.graphics import plot_classification_report
 ```
 %% Cell type:code id:de6ae734-3a6a-477e-a5a0-8b9ec5911369 tags:
 ``` python
 # entire reference period
 REFERENCE_PERIOD = np.concatenate([CALIB_PERIOD, VALID_PERIOD], axis=0)
 ```
 %% Cell type:code id:534d9565-4b58-4959-bef3-edde969e2364 tags:
 ``` python
 # empirical quantiles
 quantiles = np.arange(0.01, 1, 0.005)
 ```
 %% Cell type:markdown id:12382efb-1a3a-4ede-a904-7f762bfe56c7 tags:
 ### Load observations
 %% Cell type:code id:2373d894-e252-4f16-826b-88731e195259 tags:
 ``` python
 # model predictions and observations NetCDF
 y_true = xr.open_dataset(search_files(OBS_PATH.joinpath('pr'), 'OBS_pr(.*).nc$').pop())
 ```
 %% Cell type:markdown id:5d30b543-aa3b-45f3-b8e8-90d72f4f6896 tags:
 ### Select time period
 %% Cell type:code id:f902683a-a560-48f9-b2d1-ef9c341ca69a tags:
 ``` python
 # time period
 PERIOD = REFERENCE_PERIOD
 ```
 %% Cell type:code id:0c2c1912-a947-4afe-84a7-895726be5cfd tags:
 ``` python
 # subset to time period
 y = y_true.sel(time=PERIOD)
 ```
 %% Cell type:markdown id:f6d01e1e-9dc2-4c31-a31a-a6c91abc7fb4 tags:
 ### Fit distributions: annually
 %% Cell type:code id:0ffce851-50fc-4795-84b9-972e4f1a5169 tags:
 ``` python
 # helper function retrieving only valid observations
 def valid(ds):
    valid = ds.precipitation.values
    valid = valid[~np.isnan(valid)]  # mask missing values
    valid = valid[valid > 0]  # only consider pr > 0
    return valid
 ```
 %% Cell type:code id:6f68803b-4dbc-4d43-99c0-a32e482b647a tags:
 ``` python
 # valid observations
 y_valid = valid(y)
 ```
 %% Cell type:code id:5de4933a-ef9d-4afe-8af6-ff68d91860ce tags:
 ``` python
 # fit gamma distribution to data
 alpha, loc, beta = stats.gamma.fit(y_valid, floc=0)
 gamma = stats.gamma(alpha, loc=loc, scale=beta)
 ```
+%% Cell type:code id:dcd9bfeb-67dc-4b63-98fd-c86c3a07c2b0 tags:
+``` python
+# fit lognormal distribution
+alpha, loc, beta = stats.lognorm.fit(y_valid, floc=0)
+lognorm = stats.lognorm(alpha, loc=loc, scale=beta)
+```
 %% Cell type:code id:75b74f7c-c9d7-4d52-b140-e0ad9de17b69 tags:
 ``` python
 # fit generalized pareto distribution to data
 alpha, loc, beta = stats.genpareto.fit(y_valid, floc=0)
 genpareto = stats.genpareto(alpha, loc=loc, scale=beta)
 ```
+%% Cell type:code id:d489a3e7-7ece-440e-bbd9-1cfd739d822c tags:
+``` python
+# fit exponential distribution to data
+loc, beta = stats.expon.fit(y_valid, floc=0)
+expon = stats.expon(loc=loc, scale=beta)
+```
+%% Cell type:code id:01d8c7d9-541e-481d-b0de-e8590c571ca5 tags:
+``` python
+# fit weibull distribution to data
+alpha, loc, beta = stats.weibull_min.fit(y_valid, floc=0)
+weibull = stats.weibull_min(alpha, loc=loc, scale=beta)
+```
 %% Cell type:code id:14ade547-443a-457a-bccd-d88d049b9d81 tags:
 ``` python
 # empirical quantiles and theoretical quantiles
 eq = np.quantile(y_valid, quantiles)
 tq_gamma = gamma.ppf(quantiles)
 tq_genpareto = genpareto.ppf(quantiles)
+tq_expon = expon.ppf(quantiles)
+tq_lognorm = lognorm.ppf(quantiles)
+tq_weibull = weibull.ppf(quantiles)
 # Q-Q plot
 RANGE = 40
 fig, ax = plt.subplots(1, 1, figsize=(6, 6))
-ax.scatter(eq, tq_gamma, color='grey', label='Gamma')
+ax.scatter(eq, tq_gamma, marker='*', color='k', label='Gamma')
-ax.scatter(eq, tq_genpareto, color='k', label='GenPareto')
+ax.scatter(eq, tq_genpareto, marker='x', color='k', label='GenPareto')
+ax.scatter(eq, tq_expon, marker='o', color='k', label='Expon')
+ax.scatter(eq, tq_lognorm, marker='+', color='k', label='LogNorm')
+ax.scatter(eq, tq_weibull, marker='^', color='k', label='Weibull')
 ax.plot(np.arange(0, RANGE), np.arange(0, RANGE), '--k')
 ax.set_xlim(0, RANGE)
 ax.set_ylim(0, RANGE)
 ax.set_xticks(np.arange(0, RANGE + 5, 5))
 ax.set_yticks(np.arange(0, RANGE + 5, 5))
 ax.set_xticklabels([str(t) for t in np.arange(0, RANGE + 5, 5)], fontsize=12)
 ax.set_yticklabels([str(t) for t in np.arange(0, RANGE + 5, 5)], fontsize=12)
 ax.set_ylabel('Theoretical quantiles', fontsize=14);
 ax.set_xlabel('Empirical quantiles', fontsize=14);
 ax.legend(frameon=False, fontsize=14);
 ax.set_title('Reference period: {} - {}'.format(str(PERIOD[0]), str(PERIOD[-1])), fontsize=14)
 # save figure
 fig.savefig('./Figures/pr_distribution.png', bbox_inches='tight', dpi=300)
 ```
 %% Cell type:markdown id:5fd0e9d8-759d-45ee-bb1f-9c749ac23e8e tags:
 ### Fit distributions: monthly
 %% Cell type:code id:156e5415-4065-4887-b759-0e665d671b38 tags:
 ``` python
 # get the indices of the observations for each month
 month_idx = y.groupby('time.month').groups
 ```
 %% Cell type:code id:092e865d-f033-4f60-8098-86ae5068e045 tags:
 ``` python
 # fit distribution to observations for each month
 month_gamma = {}
 month_genpareto = {}
+month_expon = {}
+month_lognorm = {}
+month_weibull = {}
 for month, idx in month_idx.items():
    print('Month: {}'.format(calendar.month_name[month]))
    # select the data of the current month
    data = y.isel(time=idx)
    data = valid(data)
    # fit distributions
    # gamma
    alpha, loc, beta = stats.gamma.fit(data, floc=0)
    gamma = stats.gamma(alpha, loc=loc, scale=beta)
    month_gamma[month] = gamma
    # genpareto
    alpha, loc, beta = stats.genpareto.fit(data, floc=0)
    genpareto = stats.genpareto(alpha, loc=loc, scale=beta)
    month_genpareto[month] = genpareto
+    # exponential
+    loc, beta = stats.expon.fit(data, floc=0)
+    expon = stats.expon(loc=loc, scale=beta)
+    month_expon[month] = expon
+    # lognormal
+    alpha, loc, beta = stats.lognorm.fit(data, floc=0)
+    lognorm = stats.lognorm(alpha, loc=loc, scale=beta)
+    month_lognorm[month] = lognorm
+    # weibull
+    alpha, loc, beta = stats.weibull_min.fit(data, floc=0)
+    weibull = stats.weibull_min(alpha, loc=loc, scale=beta)
+    month_weibull[month] = weibull
 ```
 %% Cell type:code id:396e5ee4-1632-4591-b93b-91fa6ac1d373 tags:
 ``` python
 # plot empirical vs. theoretical quantiles for each month
 fig, axes = plt.subplots(4, 3, figsize=(12, 12), sharex=True, sharey=True)
 axes = axes.flatten()
 RANGE = 40
 for month, idx in month_idx.items():
    # axis to plot
    ax = axes[month - 1]
    # compute empirical quantiles
    data = y.isel(time=idx)
    data = valid(data)
    eq = np.quantile(data, quantiles)
    # compute theoretical quantiles
    tq_gamma = month_gamma[month].ppf(quantiles)
    tq_gpare = month_genpareto[month].ppf(quantiles)
+    tq_expon = month_expon[month].ppf(quantiles)
+    tq_lognr = month_lognorm[month].ppf(quantiles)
+    tq_weibu = month_weibull[month].ppf(quantiles)
    # plot empirical vs. theoretical quantiles
-    ax.scatter(eq, tq_gamma, color='grey', label='Gamma')
+    ax.scatter(eq, tq_gamma, marker='*', color='k', label='Gamma')
-    ax.scatter(eq, tq_gpare, color='k', label='GenPareto')
+    ax.scatter(eq, tq_gpare, marker='x', color='k', label='GenPareto')
+    ax.scatter(eq, tq_expon, marker='o', color='k', label='Expon')
+    ax.scatter(eq, tq_lognr, marker='+', color='k', label='LogNorm')
+    ax.scatter(eq, tq_weibu, marker='^', color='k', label='Weibull')
    ax.plot(np.arange(0, RANGE), np.arange(0, RANGE), '-k')
    ax.set_title(calendar.month_name[month], fontsize=14)
    ax.set_xlim(0, RANGE)
    ax.set_ylim(0, RANGE)
    ax.set_xticks(np.arange(0, RANGE + 5, 5))
    ax.set_yticks(np.arange(0, RANGE + 5, 5))
    ax.set_xticklabels([str(t) for t in np.arange(0, RANGE + 5, 5)], fontsize=12)
    ax.set_yticklabels([str(t) for t in np.arange(0, RANGE + 5, 5)], fontsize=12)
 # add legend
-axes[0].legend(frameon=False, fontsize=12)
+axes[0].legend(frameon=False, fontsize=12, loc='upper left')
 # add figure title
 fig.suptitle('Reference period: {} - {}'.format(str(PERIOD[0]), str(PERIOD[-1])), fontsize=14)
 # adjust subplots
 fig.subplots_adjust(wspace=0.1)
 fig.savefig('./Figures/pr_distribution_m.png', bbox_inches='tight', dpi=300)
 ```
 %% Cell type:markdown id:c0fea8ac-bac0-4096-bc81-90d799f8ab94 tags:
 ### Empirical quantiles per grid point
 %% Cell type:code id:a02c42e0-591c-4630-89b8-5dd8ef71a4a0 tags:
 ``` python
 # compute empirical quantiles over time
 equantiles = y.precipitation.quantile(quantiles, dim='time')
 equantiles = equantiles.rename({'quantile': 'q'})
 ```
 %% Cell type:code id:966d2724-2628-4842-abc9-695711945347 tags:
 ``` python
 # iterate over the grid points
 gammaq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan
 genpaq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan
+exponq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan
+lognrq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan
+weibuq = np.ones(shape=(len(equantiles.q), len(equantiles.y), len(equantiles.x))) * np.nan
 for i, _ in enumerate(y.x):
    print('Rows: {}/{}'.format(i + 1, len(y.x)))
    for j, _ in enumerate(y.y):
        # current grid point: xarray.Dataset, dimensions=(time)
        point = y.isel(x=i, y=j)
        point = valid(point)
        # check if the grid point is valid
        if point.size < 1:
            # move on to next grid point
            continue
        # fit Gamma distribution to grid point
        alpha, loc, beta = stats.gamma.fit(point, floc=0)
        gamma = stats.gamma(alpha, loc=loc, scale=beta)
        # fit GenPareto distribution to grid point
        alpha, loc, beta = stats.genpareto.fit(point, floc=0)
        genpa = stats.genpareto(alpha, loc=loc, scale=beta)
+        # fit Exponential distribution to grid point
+        loc, beta = stats.expon.fit(point, floc=0)
+        expon = stats.expon(loc=loc, scale=beta)
+        # fit LogNormal distribution
+        alpha, loc, beta = stats.lognorm.fit(point, floc=0)
+        lognr = stats.lognorm(alpha, loc=loc, scale=beta)
+        # fit Weibull distribution
+        alpha, loc, beta = stats.weibull_min.fit(point, floc=0)
+        weibu = stats.weibull_min(alpha, loc=loc, scale=beta)
        # compute theoretical quantiles of fitted distributions
        tq_gamma = gamma.ppf(quantiles)
        tq_genpa = genpa.ppf(quantiles)
+        tq_expon = expon.ppf(quantiles)
+        tq_lognr = lognr.ppf(quantiles)
+        tq_weibu = weibu.ppf(quantiles)
        # store theoretical quantiles for current grid point
        gammaq[:, j, i] = tq_gamma
        genpaq[:, j, i] = tq_genpa
+        exponq[:, j, i] = tq_expon
+        lognrq[:, j, i] = tq_lognr
+        weibuq[:, j, i] = tq_weibu
 # store theoretical quantiles in xarray.DataArray
 gammaq = xr.DataArray(data=gammaq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),
                      name='precipitation')
 genpaq = xr.DataArray(data=genpaq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),
                      name='precipitation')
+exponq = xr.DataArray(data=exponq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),
+                      name='precipitation')
+lognrq = xr.DataArray(data=lognrq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),
+                      name='precipitation')
+weibuq = xr.DataArray(data=weibuq, dims=['q', 'y', 'x'], coords=dict(q=quantiles, y=y.y, x=y.x),
+                      name='precipitation')
 ```
 %% Cell type:code id:601de7cb-35f4-40e1-9b51-2dab23102659 tags:
 ``` python
 # compute bias in theoretical quantiles
 bias_gamma = gammaq - equantiles  # predicted - observed
 bias_genpa = genpaq - equantiles
+bias_expon = exponq - equantiles
+bias_lognr = lognrq - equantiles
+bias_weibu = weibuq - equantiles
+```
+%% Cell type:code id:23abd0d1-7c27-4f02-b7ae-9165c2dde0b6 tags:
+``` python
+# distributions
+dists = {k: v for k, v in zip(['gamma', 'genpareto', 'expon', 'lognr', 'weibu'], [bias_gamma, bias_genpa, bias_expon, bias_lognr, bias_weibu])}
 ```
 %% Cell type:code id:b8089c11-a48d-4028-9d4b-e03101ff5e55 tags:
 ``` python
 # plot spatial bias in different quantiles
 plot_quantiles = quantiles[18::20]
 fig, axes = plt.subplots(3, 3, sharex=True, sharey=True, figsize=(12, 12))
 axes = axes.flatten()
-for dist in ['gamma', 'genpareto']:
+for dist, biasq in dists.items():
-    biasq = bias_gamma if dist == 'gamma' else bias_genpa
    # iterate over quantiles to plot
    for ax, q in zip(axes, plot_quantiles):
        im = ax.imshow(biasq.sel(q=q).values, origin='lower', vmin=0, vmax=5, cmap='viridis_r')
        ax.set_title(str('P{:.0f}'.format(q * 100)), fontsize=14)
    # adjust subplots
    fig.subplots_adjust(wspace=0.1, hspace=0.1)
    # add colorbar for bias
    axes = axes.flatten()
    cbar_ax_bias = fig.add_axes([axes[2].get_position().x1 + 0.01, axes[2].get_position().y0,
                                 0.01, axes[2].get_position().y1 - axes[2].get_position().y0])
    cbar_bias = fig.colorbar(im, cax=cbar_ax_bias)
    cbar_bias.set_label(label='Bias (mm)', fontsize=14)
    cbar_bias.ax.tick_params(labelsize=14, pad=10)
    # save figure
-    fig
    fig.savefig('./Figures/pr_distribution_{}_grid.png'.format(dist), bbox_inches='tight', dpi=300)
 ```
-%% Cell type:code id:a5ee4d3f-608b-4598-b235-3cd20a184aff tags:
-``` python
-```

--- 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",

 %% Cell type:markdown id:d6b83379-c5a8-48c3-bb85-d00a341a37f4 tags:
 ### Imports
 %% Cell type:code id:eeba7f9b-066a-4843-bd64-5b6326c0b536 tags:
 ``` python
 # builtins
 import datetime
 import warnings
 import calendar
 # externals
 import xarray as xr
 import numpy as np
 import matplotlib.pyplot as plt
 import seaborn as sns
 import pandas as pd
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 import scipy.stats as stats
 from IPython.display import Image
 from sklearn.metrics import r2_score, roc_curve, auc, classification_report
 from sklearn.model_selection import train_test_split
 # locals
 from climax.main.io import ERA5_PATH, OBS_PATH, TARGET_PATH, DEM_PATH
 from climax.main.config import CALIB_PERIOD
 from pysegcnn.core.utils import search_files
 from pysegcnn.core.graphics import plot_classification_report
 ```
 %% Cell type:code id:e75b3217-26f7-4a4a-ae2a-4fbb92a9f2a2 tags:
 ``` python
 # model predictions and observations NetCDF
 y_true = xr.open_dataset(search_files(OBS_PATH.joinpath('pr'), 'OBS_pr(.*).nc$').pop())
 ```
 %% Cell type:code id:3aa8466e-84a9-4c2e-ae19-403b6246e27f tags:
 ``` python
 # subset to calibration period
 y_true = y_true.sel(time=CALIB_PERIOD)
 ```
 %% Cell type:code id:4f1a58a2-8c4c-4d73-a116-e64e68fdd507 tags:
 ``` python
 # precipitation threshold defining a wet day
 WET_DAY_THRESHOLD = 1
 ```
 %% Cell type:code id:5e6696df-8660-4083-9a32-0dd282112948 tags:
 ``` python
 # calculate number of wet days in calibration period
 wet_days = (y_true.mean(dim=('y', 'x')) >= WET_DAY_THRESHOLD).astype(np.int16)
 nwet_days = wet_days.to_array().values.squeeze()
 ```
 %% Cell type:code id:b87accd6-d5e4-4dc6-9532-3ef8aa162d24 tags:
 ``` python
 # split training/validation set chronologically
 train, valid = train_test_split(CALIB_PERIOD, shuffle=False, test_size=0.25)
 ```
 %% Cell type:code id:559d1450-09db-4b2f-844a-d572485973e0 tags:
 ``` python
 # split training/validation set by number of wet days
 train_st, valid_st = train_test_split(CALIB_PERIOD, stratify=nwet_days, test_size=0.5)
 train_st, valid_st = np.asarray(sorted(train_st)), np.asarray(sorted(valid_st))  # sort chronologically
 ```
 %% Cell type:code id:7fd013f9-77d0-48de-8d5f-2c6a1cb3ed17 tags:
 ``` python
 # plot distribution of wet days in calibration period
 fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(10, 10))
 axes = axes.flatten()
 # not stratified
 sns.countplot(x=wet_days.sel(time=train).to_array().values.squeeze(), ax=axes[0])
 sns.countplot(x=wet_days.sel(time=valid).to_array().values.squeeze(), ax=axes[2])
 # stratified
 sns.countplot(x=wet_days.sel(time=train_st).to_array().values.squeeze(), ax=axes[1])
 sns.countplot(x=wet_days.sel(time=valid_st).to_array().values.squeeze(), ax=axes[3])
 # axes properties
 for ax in axes:
    ax.set_ylabel('')
 for ax in axes[2:]:
    ax.set_xticklabels(['Dry', 'Wet'])
 for ax in [axes[0], axes[1]]:
    ax.text(1, ax.get_ylim()[-1] - 5, 'Training', ha='left', va='top', fontsize=12)
 for ax in [axes[2], axes[3]]:
    ax.text(1, ax.get_ylim()[-1] - 5, 'Validation', ha='left', va='top', fontsize=12)
 axes[0].set_title('Not stratified')
 axes[1].set_title('Stratified')
 # adjust subplot
 fig.subplots_adjust(wspace=0.1, hspace=0.1)
 fig.suptitle('Stratified sampling: wet day threshold {:0d} mm'.format(WET_DAY_THRESHOLD));
 ```
+%% Output