{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Extreme Value Analysis using Extremes.jl\n", "\n", "This module provides an easy-to-use wrapper for the `Extremes.jl` Julia package, enabling seamless integration with `xarray` for extreme value analysis. However, do note that `juliacall` is not installed by default when installing `xHydro`. Consult the installation page for instructions.\n", "\n", "The `Extremes.jl` package is specifically designed for analyzing extreme values and offers a variety of powerful features:\n", "\n", "- Block Maxima and Threshold Exceedance methods, including popular distributions such as `genextreme`, `gumbel_r`, and `genpareto`.\n", "- Flexible parameter estimation techniques, supporting methods like `Probability-Weighted Moments (PWM)`, `Maximum Likelihood Estimation (MLE)`, and `Bayesian Estimation`.\n", "- Compatibility with both stationary and non-stationary models for flexible modeling of future extreme events.\n", "- Return level estimation for quantifying the risk of extreme events over different return periods.\n", "\n", "For further information on the `Extremes.jl` package, consult the following resources:\n", "- [Extremes.jl - JSS Article](https://doi.org/10.18637/jss.v109.i06)\n", "- [Extremes.jl GitHub Repository](https://github.com/jojal5/Extremes.jl)\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:27.904824Z", "start_time": "2024-12-11T20:54:10.703052Z" } }, "outputs": [], "source": [ "import os\n", "\n", "os.environ[\"PYTHON_JULIACALL_AUTOLOAD_IPYTHON_EXTENSION\"] = (\n", " \"no\" # To prevent random crashes with GitHub's testing interface\n", ")\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import pooch\n", "from IPython.display import clear_output\n", "\n", "import xhydro.extreme_value_analysis as xhe\n", "from xhydro.testing.helpers import deveraux\n", "\n", "clear_output(wait=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data acquisition\n", "\n", "This example will use climate data from the `GFDL-ESM4.1` model to demonstrate non-stationarity. The dataset includes annual total precipitation data from 1955 to 2100, spanning 97 virtual stations across the province of Quebec. For more information on how to access precipitation data or perform block maxima, consult the [Local frequency analyses](local_frequency_analysis.ipynb) notebook.\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:28.052959Z", "start_time": "2024-12-11T20:54:27.909829Z" } }, "outputs": [], "source": [ "file = deveraux().fetch(\"extreme_value_analysis/NOAA_GFDL_ESM4.zip\", pooch.Unzip())" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 13kB\n",
       "Dimensions:       (station_num: 5, time: 146)\n",
       "Coordinates:\n",
       "  * station_num   (station_num) int64 40B 1001 1004 1008 1009 1012\n",
       "  * time          (time) datetime64[ns] 1kB 1955-01-01 1956-01-01 ... 2100-01-01\n",
       "    station_name  (station_num, time) object 6kB 'Dozois' ... 'Lac St-Francois'\n",
       "Data variables:\n",
       "    total_precip  (station_num, time) float64 6kB 983.6 1.204e+03 ... 1.065e+03
" ], "text/plain": [ " Size: 13kB\n", "Dimensions: (station_num: 5, time: 146)\n", "Coordinates:\n", " * station_num (station_num) int64 40B 1001 1004 1008 1009 1012\n", " * time (time) datetime64[ns] 1kB 1955-01-01 1956-01-01 ... 2100-01-01\n", " station_name (station_num, time) object 6kB 'Dozois' ... 'Lac St-Francois'\n", "Data variables:\n", " total_precip (station_num, time) float64 6kB 983.6 1.204e+03 ... 1.065e+03" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.read_csv(file[0], parse_dates=[0])[\n", " [\"time\", \"station_num\", \"station_name\", \"total_precip\"]\n", "]\n", "# That dataset is a CSV file, so we need to format it\n", "ds = df.to_xarray()\n", "ds = ds.set_coords([\"time\", \"station_num\", \"station_name\"])\n", "ds = ds.set_index(index=[\"station_num\", \"time\"])\n", "ds = ds.unstack(\"index\")\n", "ds[\"total_precip\"].attrs[\"units\"] = \"mm y-1\"\n", "\n", "# Take a subset for the example\n", "ds = ds.isel(station_num=slice(0, 5))\n", "ds" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:28.202587Z", "start_time": "2024-12-11T20:54:28.105446Z" } }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAE6CAYAAAB9D9Q3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAqsxJREFUeJztvXl8VNX9//+afbJONpIQSFgEQWQRURHUihviV8RKW7W2uGG1tW5VbOuvWu0m1k8/LpVqrQtYl2ptwY9bqai4IIggi6DIZtgJCVkm+6z398edc+65d+bO3DuZZGaS9/Px4KGZ3JmcuTP33td9vTeLJEkSCIIgCIIgiKzEmu4FEARBEARBEMlDYo4gCIIgCCKLITFHEARBEASRxZCYIwiCIAiCyGJIzBEEQRAEQWQxJOYIgiAIgiCyGBJzBEEQBEEQWQyJOYIgCIIgiCyGxBxBEARBEEQWQ2KOIIikeOmll/DII4/06DUef/xxLFmyJOrxPXv2wGKxxPwd0fusWrUK1113HaZMmQKXywWLxYI9e/bobv/YY49h7NixcLlcGDFiBH7zm98gEAhEbVdfX4+rr74aZWVlyM3NxbRp0/Dee+9Fbffmm2/iyiuvxIQJE+BwOGCxWFL59gii30FijiCIpOhNMTd48GCsWbMGF154YY9en0iO9957D++++y5qamowffr0uNv+4Q9/wK233oq5c+fiv//9L2688Ubcf//9+OlPf6razufz4ZxzzsF7772HRx99FP/3f/+HiooKzJo1Cx9++KFq22XLluHTTz/FuHHjMGnSpJS/P4Lob1hoNitBEMkwe/ZsbN26Na5jk4jx48ejrKwMH3zwQcrWRfSccDgMq1W+1//Tn/6EO++8E7W1tRg+fLhqu8bGRgwdOhRXXnklnnzySf74/fffj7vvvhtbt27FuHHjAMjC/ac//SlWr16NadOmAQCCwSAmTZqE/Px8rF27Nubfv+mmm/CXv/wFdKkiCH3ImSMIIoqGhgZcf/31qK6uhsvlwqBBg3Daaafh3XffBQDMmDEDb731Fvbu3QuLxcL/MX7zm99g6tSpKCkpQWFhIU488UQ888wzqgvy8OHD8eWXX+LDDz/kz2diQS/MumrVKpxzzjkoKChAbm4upk+fjrfeeku1zZIlS2CxWLBy5Ur85Cc/QVlZGUpLSzF37lwcOnTI1H647777YLFY8OWXX+L73/8+PB4PKioqcO2118Lr9fLt4oWFLRYL7rvvvqjX/OKLL/C9730PHo8HJSUluP322xEMBrF9+3bMmjULBQUFGD58OB588EFTa04FTEglYvny5eju7sY111yjevyaa66BJEl47bXX+GPLli3DmDFjuJADALvdjh/+8If47LPPcPDgQdN/nyAIGXu6F0AQROYxb948bNiwAX/4wx9w7LHHoqWlBRs2bEBjYyMA2WW5/vrrsXv3bixbtizq+Xv27MENN9yAmpoaAMCnn36Km2++GQcPHsSvf/1rAPLF/bvf/S48Hg8ef/xxAIDL5dJd04cffojzzjsPEydOxDPPPAOXy4XHH38cF110Ef7xj3/gsssuU21/3XXX4cILL8RLL72E/fv3484778QPf/hDvP/++6b3x3e+8x1cdtllmD9/PrZs2YK77roLAPDss8+afi3GpZdeih/+8Ie44YYbsGLFCjz44IMIBAJ49913ceONN2LBggV46aWX8Itf/AKjRo3C3Llz475eKBQy5F5ZrdaUiaWtW7cCACZMmKB6fPDgwSgrK+O/Z9ueccYZUa8xceJEAMCXX36JIUOGpGRdBDHQIDFHEEQUn3zyCa677jr86Ec/4o9dfPHF/P/HjRuHoqIiuFwunHrqqVHPX7x4Mf//cDiMGTNmQJIkPProo7jnnntgsVgwefJk5OTkoLCwMOZraPnlL3+J4uJifPDBB8jPzwcgh3pPOOEELFiwAJdeeqnKHZw1axb+/Oc/85+bmprw85//HHV1daisrDS1P+bPn48777wTAHDuuedi165dePbZZ/HMM88knZx//fXX4/bbb+ev+c4772DRokVYunQpLrnkEgCyA/rmm2/ixRdfTCjmzjnnnKjcs1hcddVVKSssaWxshMvlQl5eXtTvSkpKuPhn25aUlMTcjv2eIIjkIDFHEEQUp5xyCpYsWYLS0lKce+65mDJlChwOh+Hnv//++7j//vuxbt06tLa2qn5XX1+PiooKU+vp6OjA2rVr8ZOf/IQLOQCw2WyYN28efvGLX2D79u0YO3Ys/92cOXNUr8EcoL1795oWc7Feq7u7O6n3wpg9e7bq5+OOOw6bN2/GBRdcwB+z2+0YNWoU9u7dm/D1nnzySbS1tSXcrqyszPxi4xBPzGp/Z2ZbgiCMQ2KOIIgoXnnlFfz+97/H008/jXvuuQf5+fm45JJL8OCDDyYUQp999hlmzpyJGTNm4KmnnsLQoUPhdDrx2muv4Q9/+AO6urpMr6e5uRmSJGHw4MFRv6uqqgIQ7eyUlpaqfmYh3GT+fipfi6F1qZxOJ3Jzc+F2u6Me1wriWIwaNcpwmDVVlJaWoru7G52dncjNzVX9rqmpCVOmTFFtG8t9a2pqAhC9PwiCMA5lmRIEEUVZWRkeeeQR7NmzB3v37sXChQuxdOlSXH311Qmf+/LLL8PhcODNN9/EpZdeiunTp+Okk07q0XqKi4thtVpx+PDhqN+xooZUO05mYALM5/OpHu/L0OE555wDh8OR8N+1116bsr/JcuW2bNmieryurg5Hjx7F+PHjVdtqtxOfK25LEIQ5yJkjCCIuNTU1uOmmm/Dee+/hk08+4Y+7XK6YzpTFYoHdbofNZuOPdXV14fnnn4/aVu81tOTl5WHq1KlYunQp/vSnPyEnJweAnI/3wgsvYOjQoTj22GOTeXspoaKiAm63G1988YXq8f/7v//rszWkI8w6a9YsuN1uLFmyBFOnTuWPs4rib3/72/yxSy65BDfeeCPWrl3Ltw0Gg3jhhRcwdepU7rASBGEeEnMEQajwer0466yzcMUVV2Ds2LEoKCjAunXrsHz5clUS/oQJE7B06VI88cQTmDJlCqxWK0466SRceOGFeOihh3DFFVfg+uuvR2NjI/70pz/FrFSdMGECXn75ZbzyyisYOXIk3G53VGUkY+HChTjvvPNw1llnYcGCBXA6nXj88cexdetW/OMf/0hrzpXFYsEPf/hDPPvsszjmmGMwadIkfPbZZ3jppZf6bA1jxoxJ2Ws1NDTwYgrmnP3nP//BoEGDMGjQIJx55pkA5NDo3XffjXvuuQclJSWYOXMm1q1bh/vuuw/XXXcd7zEHANdeey3+8pe/4Hvf+x4eeOABlJeX4/HHH8f27dt5yxvG3r17sW7dOgDA7t27AQD/+te/AMgtbXrq9BJEv0MiCIIQ6O7uln784x9LEydOlAoLC6WcnBxpzJgx0r333it1dHTw7ZqamqTvfve7UlFRkWSxWCTxdPLss89KY8aMkVwulzRy5Ehp4cKF0jPPPCMBkGpra/l2e/bskWbOnCkVFBRIAKRhw4ZJkiRJtbW1EgBp8eLFqrV9/PHH0tlnny3l5eVJOTk50qmnniq98cYbqm0WL14sAZDWrVunenzlypUSAGnlypWG98W9994rAZAaGhpi/g3xvXi9Xum6666TKioqpLy8POmiiy6S9uzZIwGQ7r333oSvedVVV0l5eXlRazjzzDOl448/3vCaUwHbV7H+nXnmmVHbP/roo9Kxxx4rOZ1OqaamRrr33nslv98ftV1dXZ105ZVXSiUlJZLb7ZZOPfVUacWKFVHbsf0b699VV13VC++YILIbmgBBEARBEASRxVABBEEQBEEQRBZDOXMEQQw4wuEwwuFw3G3sdjo9EgSRHZAzRxDEgOO3v/1twhYee/bsSfcyCYIgDEE5cwRBDDgOHTrE+9PpMXHiRDidzj5aEUEQRPKQmCMIgiAIgshiKMxKEARBEASRxVCGr0HC4TAOHTqEgoICGghNEARBEESvIkkS2traUFVVlXCmMok5gxw6dAjV1dXpXgZBEARBEAOI/fv3Y+jQoXG3ITFnkIKCAgDyTi0sLEzzagiCIAiC6M+0traiurqa6494kJgzCAutFhYWkpgjCIIgCKJPMJLaRQUQBEEQBEEQWQyJOYIgCIIgiCyGxBxBEARBEEQWQ2KOIAiCIAgiiyExRxAEQRAEkcWQmCMIgiAIgshiSMwRBEEQxABle10b5j2zFhv2Nad7KUQPoD5zBEEQBDFAeX3zQXy88yiGFufixJridC+HSBJy5giCIAhigNLaFQQAeLv8aV4J0RPSKuY++ugjXHTRRaiqqoLFYsFrr70Wtc22bdswZ84ceDweFBQU4NRTT8W+ffv4730+H26++WaUlZUhLy8Pc+bMwYEDB1Sv0dzcjHnz5sHj8cDj8WDevHloaWnp5XdHEARBEJlNu08Wcy2dgTSvhOgJaRVzHR0dmDRpEhYtWhTz97t378bpp5+OsWPH4oMPPsDmzZtxzz33wO12821uu+02LFu2DC+//DJWrVqF9vZ2zJ49G6FQiG9zxRVXYNOmTVi+fDmWL1+OTZs2Yd68eb3+/giCIAgik2nrJjHXH7BIkiSlexGAPHts2bJl+Pa3v80fu/zyy+FwOPD888/HfI7X68WgQYPw/PPP47LLLgMAHDp0CNXV1Xj77bdx/vnnY9u2bRg3bhw+/fRTTJ06FQDw6aefYtq0afj6668xZswYQ+trbW2Fx+OB1+ul2awEQRBEv+Dyv63Bp980YUhRDj755dnpXg4hYEZ3ZGzOXDgcxltvvYVjjz0W559/PsrLyzF16lRVKPbzzz9HIBDAzJkz+WNVVVUYP348Vq9eDQBYs2YNPB4PF3IAcOqpp8Lj8fBtYuHz+dDa2qr6RxAEQRD9CRZm9XaRM5fNZKyYq6+vR3t7Ox544AHMmjUL77zzDi655BLMnTsXH374IQCgrq4OTqcTxcXqCpyKigrU1dXxbcrLy6Nev7y8nG8Ti4ULF/IcO4/Hg+rq6hS+O4IgCIJIP+2RMGu7L4hAKJzm1RDJkrFiLhyWv1QXX3wxfvazn+GEE07AL3/5S8yePRt//etf4z5XkiRYLBb+s/j/ettoueuuu+D1evm//fv3J/lOCIIgCCIzYTlzANBK7lzWkrFirqysDHa7HePGjVM9ftxxx/Fq1srKSvj9fjQ3q5sd1tfXo6Kigm9z5MiRqNdvaGjg28TC5XKhsLBQ9Y8gCIIg+hNtPkXMtZCYy1oyVsw5nU6cfPLJ2L59u+rxHTt2YNiwYQCAKVOmwOFwYMWKFfz3hw8fxtatWzF9+nQAwLRp0+D1evHZZ5/xbdauXQuv18u3IQiCIIiBhi8Ygj+ohFapojV7SesEiPb2duzatYv/XFtbi02bNqGkpAQ1NTW48847cdlll+Fb3/oWzjrrLCxfvhxvvPEGPvjgAwCAx+PB/Pnzcccdd6C0tBQlJSVYsGABJkyYgHPPPReA7OTNmjULP/rRj/Dkk08CAK6//nrMnj3bcCUrQRAEQfQ3Onwh1c8UZs1e0irm1q9fj7POOov/fPvttwMArrrqKixZsgSXXHIJ/vrXv2LhwoW45ZZbMGbMGPz73//G6aefzp/z8MMPw26349JLL0VXVxfOOeccLFmyBDabjW/z4osv4pZbbuFVr3PmzNHtbUcQBEEQA4F2IV8OAFpoCkTWkjF95jId6jNHEARB9Ce+POTFhX9exX++76JxuPq0EWlcESHSL/rMEQRBEATRe7RFOXMUZs1WSMwRBEEQxAAkKsxKBRBZC4k5giAIghiAtPvUYo4KILIXEnMEQRAEMQBp81GYtb9AYo4gCIIgBiAszFrgkhtbtHRSNWu2QmKOIAiCIAYg7T7ZiRtakgsA8JIzl7WQmCMIgiCIAQirZh1SlAOAxFw2Q2KOIAiCIAYgLMw6tFgWcy2dAVDr2eyExBxBEARBDEBYAQQTc8GwhE5/KN5TiAyFxBxBEARBDECYMzeowAWnTZYDVNGanZCYIwiCIIgBCOszV+C2w5PrAAB4qXFwVkJijiAIgiAGIIqYc8CTI4u5li5qT5KNkJgjCIIg+pRwWMIv/vUF/vbR7nQvZUDT1i27cPkuO4pyyJnLZkjMEQRBEH3Kl4da8cr6/Vj0/q50L2VAw1qT5LvsKGJhVsqZy0pIzBEEQRB9yjdH2wEA3YFwmlcycPEHw/AF5f1f4LajkIdZScxlIyTmCIIgiD7lm4YOAIA/FEY4TH3N0kGHMJc1z2VHUY4TgNxrjsg+SMwRBEEQfco3Rzv4/zN3iOhbWPFDjsMGh81KYdYsh8QcQRAE0afURsKsAOALUpPadMDz5dx2AODVrF6qZs1KSMwRBEEQfYYkSahtIGcu3bBK1gKXLObImctu7EY2uv32202/8N13342SkhLTzyMIgiD6L/VtPnQII6N8VASRFliYVevMUc5cdmJIzD3yyCOYNm0anE6noRddtWoVbrrpJhJzBEEQhIrdDe2qn7spzJoWuJhzkZjrDxgScwCwbNkylJeXG9q2oKAg6QURBEEQ/ZdaofgBIGcuXbCcuQI3C7PKZk0rhVmzEkM5c4sXL4bH4zH8ok8++SQqKiqSXhRBEATRP/mmQSPmyJlLC4ozJztybAJEmy+IQIgEdrZhyJm76qqrTL3oFVdckdRiCIIgiP5NlDNHBRBpoV3jzLGmwYDszpXmu9KyLiI5qJqVIAiC6DO+0ebMBciZSwfiXFYAsFktXNhRRWv2kTIxt3nzZthstlS9HEEQBNHP8AfD2N/cBQAYVpoLgJy5dNGmqWYFlPYkNNIr+0ipMydJNJaFIAiCiM2+pk6EwhLynDZUFzMxR85cOmBhVubMAULjYKpozToMV7POnTs37u+9Xi8sFkuPF0QQBEH0T1i+3IhBeXA7ZC+BqlnTAyuAKBCduch8VgqzZh+Gxdwbb7yB8847T7dKNRSiuyuCIAhCH5YvN6IsH+GwHMmhnLn0EEvMeViYtZNGemUbhsXccccdh+985zuYP39+zN9v2rQJb775ZsoWRhAEQfQvuDNXlocDzZ0AKGcuXShhVqWKVZnPGkzLmojkMZwzN2XKFGzYsEH39y6XCzU1NSlZFEEQBNH/YD3mjhmUB5ddLpgjMZceWmPkzLFecy1d5MxlG4adub/+9a9xQ6nHHXccamtrU7IogiAIov/xjeDMbdzXAoAKINJFu0/Oi1OFWakAImsxLOZcLmogSBAEQSRHa3cAR9t9AGQx54oUQHRTAUSfEwiF+X4viNGahAogso8etSa58MILcfjw4VSthSAIguin1EZCrIMKXChwO+DmYVZy5vqaDp+SE5enak0iV7NSn7nso0di7qOPPkJXV1eq1kIQBEH0U745yipZ8wCAO3PUmqTvaYvky7kdVjhsigxgYVaqZs0+aJwXQRAE0eswZ24kE3NUAJE2WFsSsZIVEMOsVM2abfRIzA0bNgwOhyPxhgRBEMSApiGSLzfYkwMAcNlZzhyFWfsa5syJ+XIAUJonh1mbO/18diuRHZgWc2LF6tatW1FdXZ30H//oo49w0UUXoaqqChaLBa+99prutjfccAMsFgseeeQR1eM+nw8333wzysrKkJeXhzlz5uDAgQOqbZqbmzFv3jx4PB54PB7MmzcPLS0tSa+bIAiCMEdzhywOivNkA8DtIGcuXbBKVrEtCSDnM44clIdQWML7X9enY2lEkpgWc6NGjcJZZ52FF154Ad3d3T364x0dHZg0aRIWLVoUd7vXXnsNa9euRVVVVdTvbrvtNixbtgwvv/wyVq1ahfb2dsyePVvVRuWKK67Apk2bsHz5cixfvhybNm3CvHnzerR2giAIwjjNkTysolzZ/WHOHBVA9D16zpzFYsEF4ysBAP/ZUtfn6yKSx7SY27x5MyZPnow77rgDlZWVuOGGG/DZZ58l9ccvuOAC/P73v4879/XgwYO46aab8OKLL0aFdL1eL5555hn87//+L84991xMnjwZL7zwArZs2YJ3330XALBt2zYsX74cTz/9NKZNm4Zp06bhqaeewptvvont27cntW6CIAjCHC2R3mXFkbwsRcyRM9fXKDlz0d3JLhg/GADwwY56dPopdy5bMC3mxo8fj4ceeggHDx7E4sWLUVdXh9NPPx3HH388HnroITQ0NKRsceFwGPPmzcOdd96J448/Pur3n3/+OQKBAGbOnMkfq6qqwvjx47F69WoAwJo1a+DxeDB16lS+zamnngqPx8O3iYXP50Nra6vqH0EQRDYiSVK6l8CduWLmzEXCrAO1z9yhli78/s2vsL+ps8//Nh/l5Y4Wc8dXFaK6JAfdgTA+3J6663lfEA5LWL3rKFoHYL5f0gUQdrsdl1xyCf75z3/ij3/8I3bv3o0FCxZg6NChuPLKK1PSf+6Pf/wj7HY7brnllpi/r6urg9PpRHFxserxiooK1NXV8W3Ky8ujnlteXs63icXChQt5jp3H4+lRbiBBEES6WPHVEZz4uxX4YHv6cqAkSeLOHKuYdA/wMOsr6/bj6VW1eOHTvX3+t5kzVxDDmZNDrbI795+t2RVqfeS9nbji6bWY8T8f4PlP9yIYGjg3CkmLufXr1+PGG2/E4MGD8dBDD2HBggXYvXs33n//fRw8eBAXX3xxjxb2+eef49FHH8WSJUtgsVhMPVeSJNVzYj1fu42Wu+66C16vl//bv3+/qTUQBEFkAqt2NqC5M4BPdh1N2xo6/SH4IxdWrTM3UPvMMfeIzUjtS9riOHMAMCuSN/fetiNZU23s7Qzg2VVygWZThx/3vLYVFzz6MdbsbkzzyvoG02LuoYcewoQJEzB9+nQcOnQIf//737F37178/ve/x4gRI3DaaafhySefxIYNG3q0sI8//hj19fWoqamB3W6H3W7H3r17cccdd2D48OEAgMrKSvj9fjQ3N6ueW19fj4qKCr7NkSNHol6/oaGBbxMLl8uFwsJC1T+CIIhsg+WkdaXxosxCrE6bFblOWcQN9Jw59r7T4UwqBRCxW4udMLQIlYVudPhDWLUzfTcBZnhuzR60+4IYU1GA38w5HkW5Duysb8d1z63LiDSD3sa0mHviiSdwxRVXYN++fXjttdcwe/ZsWK3ql6mpqcEzzzzTo4XNmzcPX3zxBTZt2sT/VVVV4c4778R///tfAMCUKVPgcDiwYsUK/rzDhw9j69atmD59OgBg2rRp8Hq9qiKNtWvXwuv18m0IgiD6K1zM+dMnmsQQK4uIDPRqVuZ4+dMgZvVakzCsVgt357Ih1NrhC+LZT2RX7qdnj8JV04djxc/OlH/nD6HD3/+/Y7E/yTjs3Lkz4TZOpxNXXXVVwu3a29uxa9cu/nNtbS02bdqEkpIS1NTUoLS0VLW9w+FAZWUlxowZAwDweDyYP38+7rjjDpSWlqKkpAQLFizAhAkTcO655wIAjjvuOMyaNQs/+tGP8OSTTwIArr/+esyePZu/DkEQRH+FiaV0hsu0xQ+A0GdugIZZFWcuHWIudmsSkQvGV2LJ6j1Y8VUd/MEJcNozd2DUi2v3oqUzgBFlebhwgpzvV5bvhN1qQTAsob07qCtc+wtp/XTWr1+PyZMnY/LkyQCA22+/HZMnT8avf/1rw6/x8MMP49vf/jYuvfRSnHbaacjNzcUbb7wBm83Gt3nxxRcxYcIEzJw5EzNnzsTEiRPx/PPPp/z9EARBZBpMLKU3zKoufgAUZ84fCiMc7v9hMC2+yOeRFjHXrd+ahHHS8BKU5TvR2h3E2trMzTvrDoTw1MeyK/eTGcfAZpWdX4vFwnMCB8I0i7RK1RkzZpiKZe/ZsyfqMbfbjcceewyPPfaY7vNKSkrwwgsvJLNEgiCIrKY74sx1pTHU1BLDmWMFEIAs6NxWW9Tz+jOsJYsvDSK7LU6fOYbNasHkmmKs+OoI9jZ24ozRfbU6c7y6fj8a2nwYUpSDSyYPUf2uwG1HS2eAv9/+TOb6pgRBEESPyQhnTjPKC1CcOWBghlpZ+DsdzlyialYGE3uZ3Dz4+Uhrlx+fORIOm1rS5Lvk71tbGiqG+xoScwRBEP0YJhYyIWeuSHDmHDYrD4l1p7EIIl2VjtyZS4OY83bJ4tqTE7ualZHnkt3Sdl/mFhAcaO4CAHzr2EFRv2M5ge0k5giCIIhshjlA6XTmlDCrWjzwitYknTlvZ6BHItUfDOOCRz/GT1/qWSstkWdW1eLM/1mJvY0dcbdTnLm+/Vy6AyFeQSuK61jkRZy5jgwNU3YHQuiMpA+U5EW/F9YUeSDkzJkWc42NjfjpT3+KcePGoaysDCUlJap/BEEQROagtCbJhAII9QW3J+1JWrsDOP3B9/G9v65Jel37mjrwdV0bVnwZ3Ys0GSRJwrOrarG3sRNrv2mKu62SM9e3zhxrE2OzWpDnjJ+nmO/MbDHHHF+HzRIz/487cxm6/lRiugDihz/8IXbv3o358+ejoqLC9HQGgiAIou/IhJy5WAUQAOCy2wAEkgo17qhrQ1t3ENsOJz83m4UPWUWt1dqz69m+pk4cbJHDfiyUqUd3mqpZW7oiIe8cR8LrN3PmMlUMNXUo36tY74XlBKZjykZfY1rMrVq1CqtWrcKkSZN6Yz0EQRBECunOiD5zkQIITZjV7ZCduWTWxnKlgmEJ/mA4qT5oYi6VLxhGTgKnKhGf7FJaeCQSc+maAOGNfBae3Pj5coBSAJGxzhwrrNEJF7MJF5QzF4OxY8eiq6urN9ZCEARBpBjmzAVCEgJpGjweqwACYM5ccu7UgeZO/v/JhpDZJAQgNWJ39W5l9FXmOnPGih8AMWcuMwsgmpjjmxf7veRTzpw+jz/+OH71q1/hww8/RGNjI1pbW1X/CIIgiMxAkiSV85MOdy4YCvPWEFEFEI7kc+aYMwcAHUm2zhCrNHtaURsOS6qh7i1xxJz8ucgizh8M92lFLROZRYbEHKtmzUxni4XvYxU/AEAhbxqcmetPJabDrEVFRfB6vTj77LNVj0uSBIvFglAoMxU8QRDEQCMYliAOV+gKhHSHq/cWoqjRukE9qWYVxVxnss5ct+jM9cwh236kDY2RHC4gvjOndeN8wTAfb9bb8DCrATHHw6wZ2mdOzJmLRT4VQOjzgx/8AE6nEy+99BIVQBAEQWQwWtHQ7e/7MCtzTwrddtg1TV1ZmDUZVywVYVZxAHtPXcvVEVfOabPCHwrHF3OB9Ik5XgCRoC0JkPmtSZo74jtzBbxpcP8Ps5oWc1u3bsXGjRtpSD1BEESGox0VlY6KVl78EOOCywogzDpz4bDEq0aB5CcUiOG3Hou5XXK+3LeOLcO72+rRGteZC8X4uW8cU6MNgwHFmeuJs/X+10fw/Jq9+ON3J6K8wJ3068SiqTN+AQSfzZqhYjSVmM6ZO+mkk7B///7eWAtBEASRQro1zlxaxFyHvhOUbAFEfZsPgZASP042zCo6Tj0JswZDYaytlfvKzRo/GED8MKv2b/Vlr7kWE2FW5sx1B8IIJlk8s2T1Xqzc3oDlW+uSen482HdLrwCigHLm9Ln55ptx66234s4778SECRPgcKh34sSJE1O2OIIgCCJ5opy5NDQObtFpSwIk3zRYDLECPciZE8VcDwogvjjoRbsvCE+OA9OOKQUgizmWS65F+7f6sqKVF0AYaE3CCiAAOSTtyTHf/oX9vUMt3aafm4hEOXOFA6g1iWkxd9lllwEArr32Wv6YxWKhAgiCIIgMIypnLi1hVv0Lrov3mTMnZsTiByD5MKso5rTC1wwsxDptZCkXraGwhA5/KOZkguicub77XMyIOZfdBofNgkBIQkdErJqlLfL3DntT39KsOUE1K9v3XYEQAqEwHLb+O8HUtJirra3tjXUQBEEQSRIOS/j0m0YcX+VRNYPVirl05szFEg9KmLVnzlyy76u9OzVhVlb8cNqoUuQ4FAHk7QrEFHPpdObMhFkBOdTa0hlIugiiNVJ8cLgXnLl4NwqAkjMHyCF1I0Uf2YppMTds2LDeWAdBEASRJB/tbMDVi9dh7uQheOiyE/jjmRFmTezMmc0Z0zpzyTa1VYVZkxSE3YEQ1u9tBgBMO6YMFosFnhwHjrb74e0MYEhRTtRzopy5Ps2Zkz8PT44xYZPnlMVcMkUQkiShtUt+3qEUO3Nd/hAX4HrOnMNmhdthRXdA7nVIYk7DwYMH8cknn6C+vh7hsPpLeMstt6RkYQRBEIQxao92AICqwhPIFGeOibl4zlxyYi7XaUOnP4SuJMOsHSkQcwdbuuAPhpHvsuOYQXkAgEIm5nSKILR/y99HkzlCYYlXdhoJswJKqDKZvMTuQJi/tyOt3SmZf8tg0x+cNity44xhK3A70B3w9fsiCNNibvHixfjxj38Mp9OJ0tJSVXKnxWIhMUcQBNHHsNCZVqxpRUN6cuZYmDVWNWtys1lZmHV0eT42H/AmXQAhtqxINtTJRIJHGFzPQpi6Yk4bZu2jz6WtOwA2bMJomDW3B1MgWoX+boGQhKPtPpQXpqY9iVjJGq/fbYHLjoY2X7/vNWc6G/DXv/41fv3rX8Pr9WLPnj2ora3l/7755pveWCORArbXteEPb33FDwCCIPoPLHSmDaNGOXMZFmZljXLNCCmxx9yxFQUAgM4kxVAqWpOwvDsxN44JJb1ec7GaBvcFTPTnOW2GiwHye9A4WPv+D3lTlzeXqJKVMVDak5gWc52dnbj88sthtfbfqpD+yKPv7cBTH9fi9c2H0r0UgiBSTLOOM5cZYdZ4BRDmW5OwHnM2qwUjB+UDADqTEBqhsKRy9JJtTcIcnwJ3tJgz7Mz1kZhTKlmN547lOXsg5jRu2OGW1OXNJapkZQyUkV6mFdn8+fPx6quv9sZaiF6k9qgcljja7kvzSgiCSDVs/qk2XKkVSX0t5iRJUpy5GBddRcwZFzMsxFpZ6OYCKpkwq3beaLIhaBaqzTch5tLVmoR9TwpNtBjJ41MgzK+RFT8wesWZSyDmBspIL9M5cwsXLsTs2bOxfPnymE2DH3rooZQtjkgNkiThQJN8AmQ2O0EQ/QcmmLSiRisa+jpnrsMf4pMaYhZAsDCriRAnK34YWpzDE9+TEanaRrLJhllZ+K7Arbw/8zlzqXfmnlu9B098sBsvXDcVo8plB5N9T4pMiLn8SM5c5jlz+s2oRQbKSC/TYu7+++/Hf//7Xz6bVVsAQWQeLZ0B/kVuiTNihiCI7ISFnLoCIdXUgXTnzLEcXafdipwYg+TdrADChDPFnLmhxbnIdSbvzGnDbskWIcTLmdOvZu39nLnlW+tQ19qNj3Y0cDHXaqJhMCOvB/NZtTlzh1PozLHvVgnlzAFIQsw99NBDePbZZ3H11Vf3wnKI3mBfk9Jgk92ZEQTRf2COuyTJwoAVFjAnzmKRf9fXYVZxlFesm/2eOHPVJYozl4xrpBUnPc2ZKxTCrIWJwqxROXOp/1yYKya2qzHbMBhQxFxyzpz8nEK3Ha3dwZT2mmuKE74XKWBitJ+LOdM5cy6XC6eddlpvrIXoJUQxp1ddRRD9hYMtXXjk3R1oNJkf+tWhVnwWGZaeTQRDYZXrIIZSmePDLt5dfdicFkjcoT+ZAgglzJqbEWFWJgrNFED0RTUrF3NCg2UWmfGYcOZ4NWsSvfzY9WZsZSEAoK43nLlEYs49MHLmTIu5W2+9FY899lhvrIXoJVTOHIk5op/z9Mff4JF3d+L5T/cafo4/GMb3n/oUVzz1ada179Ee010qMSf/P8uR6uucOSbm9MJ6PSmAkHPmelAAoXXmki2ASKY1SR/kzLHiA9GZ49WsBqc/AD0rgGB/b0yl3ELmSGs3gilqkGy0NclAqWY1HWb97LPP8P777+PNN9/E8ccfH1UAsXTp0pQtjkgN+1VhVhJzRP+mvk125PY2dibYUmHLwRZ+4TnY0pUwdJNJaI9pUdgwkVSU6wQaO/tczClh1tj7UwkHG7vAiz3mhhbnIBgprkgmF1CbEN/TalZTBRCR95vvsqPdF0x5mDUclrgT1dMwayoKIEYOyuPzauvbfKiKMeLMLEZbkzDHtLWfh1lNi7mioiLMnTu3N9ZC9BL7haHUrd0BhMJyjyaC6I94O6PDS4n49BslvJpt7Xu0ebCisGGOD3PG+rwAgjtzqQmzij3mKgvd3J3p8AdVhR9G6BDCo23dwR5Us8rfN73WJLHWxYRjoZuJudQ6cx3+IMKRSQ9NHX50+oPIddrh7YrvlMaiRzlzXcrosIpCNw40d+Gwt6vHYk6SpLj9C0WYyO7vOXNJjfMisgsxzCpJ8smnPw8cJgY2LZELlnZOaTzWCrlyje3ZFWZt1jhz3XHCrOksgIiFy+QECBZiHexxw26zIieSM6ct/DACu7iX5btkMZekO9berZ8zF4w0Js5zqS+17P0W5jhwyNudcjGndaEONndhdEWBEGbto2pWXhziQJUnBweau3CopRtThpl+KRWd/hD8kX2WsGmwi7Um6d9RKRrj0M8JhMI41KJOOqVQK9Gfae6Qv991BvNzgqEwPt/Tj5w5QbB1B4QwK6JDieGwxN2t3sBoAYQ/GIbEhobGQewxB4DnzAHm8+ba/UzMyWtLNm+N95lzKQIp12mDPRL9iBVqZZ8DE32pns2qzdU7ELmxYed+M02DUzHOqzDHgcFF8kzWwymoaGXfWZdOyxuRwgHSmsSQmDvxxBPR3Nxs+EVPP/10HDx4MOlFEanjcEs3QmEJLrsVgz3ywURFEER/hl08Q2EJda2Jq+e2HmpFhyAEsk/MaQog/NHOHK9m1QieX722BVN+vwJfHWrtlbUlCoWJTpoRd0rsMQcANquFC8JOk9WWzFErzXNF/n7qxnlZLJa4eXOiMyf+nCq0Yo6lHLT0oM9cRzITICL72JPjwGCPLMC15kIyiPlyiULrvACiO2johiFbMRRm3bRpEzZv3oySkhJDL7pp0yb4fNl1QuyvsBBrdUkunDYrDnu7+7zX3NINB2CzWnDxCUNS8npbD3rx/Jq9uH3msagodKfkNYn+QSAUVoWDDjZ38Qu/Hmu/aVT9fDTrwqz6zpxSAKGEWcUcri8OeCFJwOYDLRhXVZjytbUYdOYA2RlLFCZl5zPmzAGyC+YLhs07c5HvSVkBcy3NC6pQWOI3AmLOHCALmMYOvzFnLsUFEFFh1pYudAeU0KSZNJv8iPvpD4XhD4bhtBsL6EmSpDhzbgc3E1LpzCWqZAWUnLlgWEJ3IMxD8/0Nwzlz55xzjmFVS5MgMgd28qspyeUnEL0Kq95gV30bbv/nZtitFpw3rkIVFkmWZz+pxdINBzGqPB8/+tbIFKyS6C9ov9tG8uZYvtzYygJ8XdeW0c5cW3cALrtNdUGNak0So5qVXfTCknxRdtnlCxpz9epbe+c98zBrXmwnyG61wGqR1yULmviOEatQHl6axx/LddrR3BkwLeZY2JA5c8lUs4o3DgUaMRevcbBSANF3zhz7rG1WC/JMCJo8l7Jthy8Ip92YEOwKhBCMVGEU5tgFMZc6Z07veyWS67DxptltvsDAFnO1tbWmX3jo0KGmn0OkHlbJWl2cg4bIRaovc+ZeXX8AgHxX5O0KpETMsRNVE02zIDRov9uJKlpDYQnrIvlysycOjoi5zPxeNXf48a3/WYmxlQV49cfT+ePxcuZYLpbYJLbbr4g5dizVt5m7wC75pBbVJbk457iKuNu1dLCwXmwBYLFY4LLb0BUIGRI0/Oa0VHFbWeNgs2FWlkNVViCLuWBYQjAUht1mPJWciTmnzcr3KcNImFXJmUt1AUSAr8sfCuNgSxcvDCrKiT2NQw+7zQqX3QpfMIwOf9Bw2x72vu1WC3IcNl7BKoZZP9l1FHsaO3DFKTWm1tTUEb/ljYjVakG+S65YbusOorzA8J/JKgxdWYcN62HpCdEjJElOUi7Nd5l+rhhm9Ud6MpkRc3sbOwAAw4Q7YaMEQmH8e4OSO9nWHcRgj+mXiYKdhGmaBaFFK2wSOXPbDreirTuIApcdZ4wehD+9s8OUMxcMhfH21jpMHVHS6yH/9Xub0dYdxOd7mxEIheGIiA5W8MEuuLHCrPkuO+xWC4JhCV2BEDxwyJMjImKE9eYzwp6jHbjvja9Qlu/E+rvP092u3Rfkrx9v37gdVnQFQgmdse5AiLs6w0qixZzZtitsokGZIE66g2HkmxBzsfLlGPEaB7OQridHfl7Kw6yRliCjK/Lx5aFWHGzu4i17zPSYY+S77PAF/aby5tgaCiPikTlzR9t98AVDaOsO4rrn1qMrEMLIsnxMO6bU8Gu3GOwxxyiIiLn+3J6EqlkzHEmSsODVLzDl9+9i1c6jpp+/XwizsrwZdoeWCF8whDmLPsHFf/kEgSS6dn+4vUF1YUyV+GIn4f7eBJIwT5Qzl0DMsRDrScOLueBo6vAjHDaWUvLe1/W45R8b8fu3tiWxWnNsPegFIIckjwiFHSzMypyPbn+0MydW/TGxJx4/ZsQcKyppTLCf6iK5UQVuu2o6ghbmaCVy5ti5LN9lV13EWdisw2zOHCuAEG6SzYZaY7UlYcR35iJh1t4qgIiIzOMGy3mQR9q6eWTGzCgvRjLtSVo1M2tL8pw8R/KI14enP67l38X/fllnaj1mcuYAcaRX/71mkJjLcF5cuw//3iCHKrcdNl9xJoYlWG8hr0Fnrr7VB29XAC2dgaSE2D/X71f9nKoDiZ1A+zL3j8gOmLBhF41EYVZW/DB1ZClKIy0qQmHJcMX3oYhYPGyip12yfHnIy/9fNW8z4lJURsRorAkQLrsNbo2DJR4/DQaqfhnsQirnIOkf08xFY46MHi6HscbBLF9uWGmuKiTHUje6zFazRlymAred5yDGE3Ob9rdg+dbDqsf4KC+zYi6gDrP6eylnbkRZHlx2KyQJ+PpwGwBzPeYYzP00055EbEsCQOXOfXXYi7+v2cO3fefLOlOVpkanPzCUkV7995qRVjH30Ucf4aKLLkJVVRUsFgtee+01/rtAIIBf/OIXmDBhAvLy8lBVVYUrr7wShw4dUr2Gz+fDzTffjLKyMuTl5WHOnDk4cOCAapvm5mbMmzcPHo8HHo8H8+bNQ0tLSx+8w56x9aAXv33jK/6z2btGJsQAoLpYdOaMfaEbBFfNbMPIo+0+vP91PQDlItOaokHH7CRMYVZCCxM2YyOOxMGWLt2LRDgs4bNIvtzUESVw2Kz8GDEaamWhpL64499yUBFzh4SKQHZhY328YoVZYzlzKjHX7jN8MW0U9k284eWHW5iYi9/tn0+BSJA3trdJEXMiSs6c2WrWyOQGlx1uLub01/DTFzfgxy9sUI1H5KO8XNECSU/MhcIS/KFebk3SrYRUh0Qqf9nNQLJhVsCkmBMaBjPYd+HB5dvR6Q9hbGUBcp02HPJ2q77fieDOnNEw6wAY6ZVWMdfR0YFJkyZh0aJFUb/r7OzEhg0bcM8992DDhg1YunQpduzYgTlz5qi2u+2227Bs2TK8/PLLWLVqFdrb2zF79myEQsqBfcUVV2DTpk1Yvnw5li9fjk2bNmHevHm9/v56Qmt3ADe+uAH+UBhs8pbZ7u3spFOW70Seyw5PZLiy0dYkDW3iSdvcQfDaxoMIhiVMqi7CpGo5US5V4oudhFMlDon+A7twjhtcAItFvkjqFTTsqG9DS2cAuU4bxg+Rv6OlkYvDUYNhRyZm4omaVFDf1o0jQsUpc+a6AyEuQKoiF0q1mIuEWR2KmGM3heJ5IBCSoiZJ6NEoNBlmYjYWRp05Pp81oTMn5+/WlKjzd5MRc8FQWDUf1a3ZN1rCYYm31BBD97FGeTH0xJzowvVaaxIhX21IERNzcmQnmek/SYVZu5Qecwx2w/HNUfmz/Nl5x+LMYwcBMBdqbeYFEMaEKROj/TlnLunSwvXr12Pbtm2wWCwYO3YsTjrpJNOvccEFF+CCCy6I+TuPx4MVK1aoHnvsscdwyimnYN++faipqYHX68UzzzyD559/Hueeey4A4IUXXkB1dTXeffddnH/++di2bRuWL1+OTz/9FFOnTgUAPPXUU5g2bRq2b9+OMWPGmF53X3DXv7dgX1Mnhhbn4IzRg/CPz/aZ7oOkbbBp2plrS86ZkyQJr6yTQ6zfmzIUm/e3AEjNXZF4Eo53ISEGJsylGlTgRnmBC0dafTjY0oVBBdHFQ29sll3+k4aX8GKCsnwXdjd0qFzpeLRyMde738UvNU19D2o6+tusFpQXqltshMISApGiJ3eCMCsgC0YjYStxYkS8G6q6VnmNKXPmGvWcOfkyZqaaVUzkzxPEnJ6oEmediq5tvJw5vdYkomAUW5OYnS0bDzFfjfXkY3mRfebM8TCrsm+qhO/C2MoCnHdcBbr8Ifxnax3+++UR3Hn+WEOv3ZSgf6EWypmLwYEDB3DGGWfglFNOwa233opbbrkFp5xyCk4//XTs378/8Qv0AK/XC4vFgqKiIgDA559/jkAggJkzZ/JtqqqqMH78eKxevRoAsGbNGng8Hi7kAODUU0+Fx+Ph28TC5/OhtbVV9a+v8HYG8NYWOTdj0RUn8jCl2dmBYo85QBFzRnPmYp20jLD5gBc769vhsltx0aQqfiClwkkTk5xbI0OsCYLBxE2R4EjEyptr7vBjySd7AAA/nFrDH2dtKoy2J2E3FO3+oOGiiWT4MhKCckZE58FICJMPss9xcOeNOVSiMJGdOfm5vABCK+YM9ppTO3P6x/ShFoM5cwYLIPbphFlzknDm2JxOp90Kp90KtyN+mFUUZEdjRCwKYhR48POtZh+x92m3WpDrUmbLMuGdCriYE44DRjJijvWaM1NkEjPMWqR8F24+ezSsVgvOGlsOu9WCXfXt2N3QnvB1JUkyX83KR3r132iOaTF37bXXIhAIYNu2bWhqakJTUxO2bdsGSZIwf/783lgjAKC7uxu//OUvccUVV6CwUM6Hqaurg9PpRHFxsWrbiooK1NXV8W3Ky8ujXq+8vJxvE4uFCxfyHDuPx4Pq6uoUvpv4sBN0vsuOE6qLlBONyZyQKDHHwqwGRZAqzGoicZQllZ81phyeHAe/M4t1V7S/qdPUHbXoEPpD4ZTnmhDZDR8knuvAkIgjfbClM2q7pz7+Bh3+EMYNLsR545ReaYMilY2NJp05SVJmffYGLJ/otFFy+4aDEde9RRiXlaNx3kSny2mLzpnTVv4arWhtahedOf33XMfCrEWJxFziAohgKMwjDdo2SXlJtCZhzhwTYYnCrGIUQBSzSmsS/Zw5reBlf8PtsKknYKQw1MrDrG4lZ45hZpQXI5kwq1dTAAEAYyvla/eYigJcML4SgLyfWFsSI6HWdl+QC1/DzlwS6882TIu5jz/+GE888YQqPDlmzBg89thj+Pjjj1O6OEYgEMDll1+OcDiMxx9/POH2Wrs6lnWdyNK+66674PV6+b/edh1F+EHg1pxoTDtzsiOhdeZCYcnQl1oVZjXhzLHcG3YSZ3dm2pPa/qZOnPk/K3HD858bfm2tzU9FEISIKG70nLmmDj+eW70HAHDbuaNV5wGeM2dQzIk3KL0Zwtl6UI4MnH+8fAE81NKtciiKcp1ROXGiA2S3WbnY05sEc8RgRWuTUWfOy8KsBnPm4oRZD3u7EQhJcNqsPFLByOFhVuPnR5Z3y0SKIihjr0GMKojfDVYAkShnTrx5Zudxl93KndZ4f9ss4bDERWZhjh1DitROZjJiLrkwKxOUyr6ZMqwYL143FS9cNxVWq3Lcse/1f788kvB1Wb5cjsNmeJqD4syRmOPU1NQgEIg+gIPBIIYMSc3sTZFAIIBLL70UtbW1WLFiBXflAKCyshJ+vx/Nzc2q59TX16OiooJvc+RI9BekoaGBbxMLl8uFwsJC1b++QrTIAQgnaXMH+36hYTCgvhM00tYj1knLCF7eaVy+MOpVEu040oawBHxd12b4tbUHIxVBECK8y32ukzsS2l5zT0dcueOr1K4ckESYVfj+9VZydXOHn7+HcyPr7QqE0NwZ4DdOxaIzF1CHWdkxz0STNmfOHrmoNhh05ho7hN6ROsdfuy/Ij9VKozlzcW5WWb7c0JIc2Kzqm/BkJkCwivh8g86cqvK3TXTmEveZC4QkzVSOMP+b8gSM+ELSLGJ+XyxnLrkwaw/6zGn+3mmjyqJyWGeOq4DFAmze38IdXT2UfDnj7yOf5cyRM6fw4IMP4uabb8b69ev53cb69etx66234k9/+lNKF8eE3M6dO/Huu++itFTdIXrKlClwOByqQonDhw9j69atmD5dHnczbdo0eL1efPbZZ3ybtWvXwuv18m0yDbESCVD6MJlpTRIKSzwsUV2iHMy8CMJA3lxDkjlzYqgLUN6HNl+BXYhaOv2Gc9+0d4bUa44Q4eOjchwYGnHmDgjOnNqVOzbKnS/LZ2LOaGsS5fvXW/k4rPhhWGkuyvJdKI9cCA82d6nEqzaMytuSRB6PCrMKvcgAYyO9wmF11ateEZLRhsHy+hIXQOxtkqsfh8eYRJNMNSs7n+VzZy7+DbP4OYtiVvs62nUxoSyep9h5nL1vpQAkNWFWdtMs5wLaUFHgUglg1tXADHkpak2iR3mhG5OriwAA73wVP9TKUiCMtiUBBkbOnOlq1quvvhqdnZ2YOnUq7Hb56cFgEHa7Hddeey2uvfZavm1TU1Pc12pvb8euXbv4z7W1tdi0aRNKSkpQVVWF7373u9iwYQPefPNNhEIhnuNWUlICp9MJj8eD+fPn44477kBpaSlKSkqwYMECTJgwgVe3HnfccZg1axZ+9KMf4cknnwQAXH/99Zg9e3bGVrIqYVb5IHBrTsJGaOn087wCMSxRlOPEkUgz4HhIkpR0axIx1AUIzpzmbzZHwjWBkIQOfyjhSR+IvjOkilaCERDGUxXlOvmQb9GZY7ly44cU4tzjonNpyyKNgxsNOHOSJKnc5t4K4WyN9AcbXyW3T6kqykF9m1ylKxZ86OXMsR5qen3mjq0owM76dkMFEN6uAEJCoYeeM2e0LQlgrABiX6M6/1ckN4kwa4cmPOpOcMOsKoBoj84ljiVYLBYLPDkONHb44e0K8KpesZEzEBHb3UHV+/d2BrD1kBfTRpaqwpFGaNVcP+yR0DQ7DpILs7KmwcmM8zImM84eW44N+1qwfk8zrpw2XHc7PqKyOPq7oEcBtSaJ5pFHHknZH1+/fj3OOuss/vPtt98OALjqqqtw33334fXXXwcAnHDCCarnrVy5EjNmzAAAPPzww7Db7bj00kvR1dWFc845B0uWLIHNpsTSX3zxRdxyyy286nXOnDkxe9tlCmLDR8BYTokWdtdd4LarBkd7DDpz7b6g6u+ZsdfZazNHrlCnLLxZ6HPV3OFPTsz14zstwhzizUKh284dj7buIFq7A5Ak4Pk1ewEAt54T7coBijPHmujGy6vtCoQMCZuewoofjh8ip3oMKc7Bpv0tONjSxW+IivMEZy4ianhuFnPmWM6cX13NOqo8H4CxAggx+V98DS1GGwYDijMVL/Kwh8+IjiXmzBdAMNGfpw2z6oR6RdEuCv14rUkAKGKuM9qZc2udOUHM/eaNL7F040EsvuZknDUm+qYjHkrhgbKmIcU5XMwlFWZ1Jh9mNfr3jq0oAAB8czR+ReueSI+6YWUmxNwAaE1iWsxdddVVKfvjM2bMiBteMxJ6c7vdeOyxx/DYY4/pblNSUoIXXnghqTWmA+3ByE7SZmx4rTvGYKNcEs1n1eYMmTkIeJg18rd4VZdOmJWtt7ok8Wtr76yoAIJgaG9g7JGJDi2dARxq6cJ72+rR7gtiTEVBTFcOUMScPyi7fPFCRFpXuLcuFKwtyYRIY2NW2HGopYu/Z4/QmqQrEIIkSdyZi8qZ01Szjq5gYq47oYBt0oq5VDhzjsTOnF6POUCczWoiZ04THk3UmkQ8z3T6Q+j0B5HrtMcd5wXE7jXXHWSOacSZixFmZdMuxGkTRtE6cwAwtCgHn0Gu/HXYTGdXmS6AkCQp5jricUzkpmJ3fQfCYUnXkdwT+S7ECrnroYzz6r9iztCnKvZY0/ZeS1cvtv6M9iBwa/pDmXkN7V2R0Zw5bTK0mZl2Ss6cugCiOxBWdT8XO9A3G5xKQTlzhB4tvBhAyaVhwuebhg4sjvSVu+HMkbqCJcdp460uEk2B0AqZ3hBzrd0BfvE6vkot5g42d/FjqDjXyUVNWGJte9QFEIrYk49BduyMLpcdke5AOGGCeFOHep/o5swZbBgsrk+vAEKSJKHHnH7OnKnWJH61o8aEld4Ns/am8WikCIKP89IRLLGmQPiicuaixWxbnGbUdd5u7DiiXzTGXESx8IAVQSQz/QEwnzPX4Q8pRRgGnbmaklzYrRZ0BUKoi1NZzSaBmBFzBYKYC/ViP8h0YkjMFRcXo75enrNZVFSE4uLiqH/scaLnsIMxOsxqwpnTVJQy4g1/FokWc8YO4kAozLdlzpwYPhUTUJuTEHPRYdb+e6dFmENp0yFcxCLCZ9H7u3C03YchRTm4aFJV3NdhFa3akKIW7QU+2eTq5g4/bv7HRny8syHqd19Fih+GFOXwBqlVRUqVrljNys4TANDtD0flZok5db5giN8cVha6eU5Rorw55tgzB1PPmTPaMFheX/xqzqPtfnT6Q7BYwKcZiOQlkTPHRBJ7rplqVgA42uGDLxjiN6d6KSKxzrdRzpwj+v3Hm/l77ZJ1mP3nVaiNhBu1KGaAEGaNfGeMCistZqtZ2RqcNquql148HDYrd171mgcHQmFe0DTcRJhV/HzMOLjZhKEw6/vvv4+SEjkGtnLlyl5dEBHdbDGZ1iTMpfBow6y5xuazsiTf8gIX6tt8hl0H8aTF1m+3WZHntKHDH0JrdxClkQsB6xckrjcR7GTitFnhD4VjhllDYSmqfQHR/+Hf+RiOxFeHZVF03RkjEoaZyvJd2NvYmdCZ0x4TyTpzL6/bjzc2H0JLpx9njB6k+t1Wli9XpbRGEsOszGH05DrgsFnhsFkQCEnoDARVc1kBqPrQsePUYpFdi0GFLrQ1BFHf1s1z6GLBwqwjynJxtN2nm+ZgtGGwvD7misU+vzEnpsqTw4WpiNiSJV54TkS/ACJxnzlAdm3FlI9EYq41jjPHes2JzmS8mb+7GtrhD4Xx1heHcNPZo3XXKgq3qSNLkeOwYfoxpVHbG4G9v05/yNDYMTFVyMyIsmMG5WN3Qwd21bdHHQuA7EYHwxJcdisqChJ/txhuh41fM9q646dPZCuGxNyZZ54Z8/+J3kF7Z8VbkwSNHUiAcjBpw6zsZ6Nh1hFleajXnLiM/N0Ct10lqApzHOjwh3SdOW0ujrxGPzw5DtX7ZWKu0uPGvqbOqDvmXy3bguVb6/Dfn32LuwdE6giGwvhgewNOHFZseJROb+EPhuEU7vpbNOF9AKpRRsW5Dlx2cuJJLkYbB0eHWZNz5j7ZdRRA7GOAuS9jKgv4Y0ygNnb4+THGQstuhw2BUFB23+LkzImpHFarBeUFLnzT0JGw1xxb4/DSPKzb04w2XzCmgDLaMBhQqm31ig/i5csBSpgVkN9bnolCKlalmagAgp1nchw2dAVCaOzwC+6eTffmMWaYNcqZU4vZYCjMx2ZpnbDugOIG/mdrXWwxJ0x/YIwoy8Pme2eqjhczsHFewbAEXzCscoFjYTZfjjGqPB/vfHVE15nbI4RYzVb55rvtaOrwR47TxOH/bMP0J7t48WK8+uqrUY+/+uqreO6551KyqIGOXjWrJBlvLCm2LBDhOXMGw6wjB8l5CUabLeoVXihTIOTXkbvXi86c+kL2wfZ6nPDbFXj8g92qx9kddVXkjl97QX132xE0dvijBpMTqWHFV0dw3d/X4/63t6V1HZv3t2D8ff/Fn9/byR/zCnNKGWJY7sppw3kbi3iwMGtDgvYkWlcqmeTq7kAIn+2RWzjFusFij5UKwrlQ6N3G8n+YmMsVXCqln1l0mFV7s1cRaV+UKMzKQs/DI73pJCk6bGWmYbC4Pl1nTmcmK8MtuHVGQ62KmGO9PI2N82Lnw6NCtEIvXw7QCbMmqGaNN1VE/PnLQ628ZYtqrcL0B5FkhRyghKMBY3lzLP2lwGRY95hBShFELBIJ+3jwvLl+mppj+tN94IEHUFZWFvV4eXk57r///pQsaqDTqgmziiereI01Y71GVAFEJIfOm8CZY64EayjqD4YNzQ7UTn9gaJs2dvpD8IeU99KsWc+GfS0AgI2R/zLYyYzlDYkJ2KGwxHN6+usBm25Ye4Nd9YkHYotsr2vDi2v3pmwY/bo9TfAHw6oGo2L+GGNEmXxxcDusuGr6cEOvXWZwPmurpiIymfzN9XuaudMSK4+1WRjXxbBYLCrHURwUL4ZSlZy56N9pb7pYI+JEjYNZAURVkZuLA+37NtMwWFyf3vllX8SNqSmJnfButVpMF0Gw8wNznLg7qHN+ZZ/NyIjYONru4z3m9CpZgfjOHBOQ2vevFnPx3d9Ys0yTdcXiIe5jI73mYuXtGYFXtCZy5sqMFz8w2Hexv7YnMS3m9u7dixEjRkQ9PmzYMOzbty8lixrIyCXd6gIIh83CbXyj81mVkJOeMxffdWDTH8TqMSMCSS+8y4Qpu2vUFjxof26IXFS0jzMXgF3MRGeuscPHnQozncoJ4zBHw+gcT8ZdS7/Ar5ZtxZpvGhNuGyvcqIU5x7vq27lA5G06BOEzprIA918yAc9edbLhsPCgfHNhVvZdTOYi8fEupehBHiCuFhN6TneVkItWnKukIigju2IUQAhhVu1xWh7JPzqSyJmL3CyV5Ll0Zy6baUsiry9+AcQeA24MnwIRMPYZtPMq1MQFEP5gmBeLjIyIiKNCmFWvxxyg5Cw3x+ozZ49dzSqe06LHF6p//s/Ww1F/U2+MVk8xUwRhtsccgzmf9W2+mMU1qXDm+utIL9Nirry8HF988UXU45s3b44at0WYpzsQ5o4VOxgtFgs/8I3eebKwpXZ0i9mcucpCNz9RGjmI9QovlCkQ8muIxQ+x1sP+frPmws7uClnLA/GOVwwRZcsB+/neZlz/9/UxwyWZCBPJ9W0+wyX+kiTx+btsxJwez66qxYm/W4G3t0RfpETY96M7oFS3tcQIswLAFVNrMH1UdDRBj9J8Y/NZ2XeZCatkcuZW7Tyq+lnrzrXEcOYAqOZtiq1YcoQ5pVGtSZxKiyNtkVV5oTFnjoVZS/OcPIwXJeZMNAwG9IWUJElYuuEAL16JNf2BkWPCNZK3i900OJagFEUFi1SIBRDx3EcWHhdvULThb+04s3hijn3HmHu8YV/0LNNYA+5TQZ6Jfn7akZRGKXQ7uEv8TUN0qJU1DB5hoi0JQ2kc3D/bWZkWc5dffjluueUWrFy5EqFQCKFQCO+//z5uvfVWXH755b2xxgEFO5BtVgs/eAD9BN3lWw/j929+FXVh1XPI2B2+LxjWzQ+RJIm7EoMKXKbsab1cvULNgZTImWPd6Js0jyth1kjOXFeAN5cWL0TZ4sz947N9eOerI3hzy6F0L8UQbEB5KCwlDEMy6lq7eS5TIoG0aX8LAKWKUw9xbvDOelkoamcCJ4vR+azcmStOzplrEnI7WUWjNneUuY3aoeJDihRhIx7jovvGCyAcmgIIfyjqOB3Ew6z671mSJH5zVZrv5BdHrVuUCmeutTuAW1/ehNv/uRn+YBjfOnYQxg0u1HsJ5Drs/L0Zoc2nFmLxxnkxsVrgsnPRe7Tdx89l8cKZ7MZAFHPa8Df7rz/EJnMo+1NvfOHIsjycWFMEIHqWaSY5c8mEepW8OXWoNRgKY3/kZnBYEmHWC8ZX4qazRsX9HmUzpsXc73//e0ydOhXnnHMOcnJykJOTg5kzZ+Lss8+mnLkUIOYaiFWceiO9HvjP13h6VS027mtWPa53Yct3KVWmeu6ctyvA57rKJ23jB7He3+V38ZETPxNvbG6snjOnnQXJRBoLbYUl8MovMUSULZ2+OyN3uZ0mZh6mE1Ekx2vsKSLeYSeqlmTfi0RiXHydHUfkk74SkuxZla3R+aza/M12X9DQ1BoGq2IdW1nAW3iIx4EvGOIiWJuDqg6zKr9jLro6Z04dZvUFw4Jzrw6zNsQJs7Z2Bfm825I8J3d+tM6cmYbB4vrYeo+2+3Dhnz/G65sPwWa14I7zjsXiq0+OW72Y61IcyUT4g+Go/nDx5sOKLibPp+zwR4VqY8FC+/J4xMiYNV4AwXLm1AUgonOk/U6x3xW47Zg1vhIA8J8tGjHXCzlzQHTj4GdX1eKX//4ipkOv5H2bdwePKZeFmjZv7rC3G4GQBKfdisGFxtuSMOaeOBQLzh+DyTX9sx+uaTHndDrxyiuv4Ouvv8aLL76IpUuXYvfu3Xj22WfhdKa3VUF/QBv+YOjdOTJxJAoZsVJUK6osFkvCkV7sQunJccBltyE/clLoSc6cchcv/56tj4Ut2n1BfoINhyW+BklSXjMYUnJXSvNd3M1gJw4xjytbxBwT52YanqYTcb9qwzt6fCOclBO5Xex70Z5A3IpijjlzzTGaBicDq2YVL8CxYN87dmMRCkumPkcWYj19VJlyTApijhUpWS3RgkGs0i3OU96v6L5pqyZzBKefnS88mjBrmy+o6241Roof8l12uOy2qDxYhpmGwYAYZpT/7r8/P4D9TV2o8rjx6o+n4eZzRifsGylW8SZCvFHIM+LMCRMVmJhr6QygKZIqEi/MWui2w2GT185C1HrOnJIzpy7qEr9TbcJaLhg/GACwtraRu+ThsMRdx2SEVDzEkV4Hmjvx+7e+wsvr9uOLAy1R23p7ICi5M6cRc3t4IUyu6bYkA4Gka5WHDx+OiRMnYtasWRg2bFgq1zSg0bOntXMVGezEJF4kO/0hfgcdKwHVk2CkF7tQstAL6w7fZmCkl5K3pBb22tYkvF9VWS6YAcnEZUtXgK9f3LZDOKnlu+z8ZOXlYk5w5rKkYoldOLsMJm6nG1HMGS2C2G3CmeOfdRwxHgiFVeH3nUfaEYw0AwWiQ/xmKXDZeaVmvPWyY7Wi0M3FhtFQqyRJWBVx5k4bXcaLNsSWQc1CE2TtxatKqGYV82KZ+9YZw5kTq+IPRz47dn4ocNm5oNHLm2OfDXObtMc0w0zDYEBMIZHX+8luuUjm2tNH4ESDLkqOgwmNxGKOfYfdDitvIB2vAEK5QbWjKMfBP+t9TfL3Ol5rEovFwvdXU8TpjXLm+AQIFmbVb3nTKjhz1SW5OL6qEGFJbhkEAO3+IJiR11vOXLsvhBc+3cfHdcU6D/Qk1KuIOXXOHMuXMzPGayBhWsx1dnZi/vz5yM3NxfHHH88rWG+55RY88MADKV/gQEPP2eIhEuFkEwwpFWviRYddEJw2K3+eSCwXQITlI7FwE7sjMyKQlIrC2AUQbdyZUy4M2qIM7cWEbcunP9itcNqtUdV0DVmYM8dyIM3MlUwnSYVZhbFDiZ05Jtz1P7/Gdj/EaOau+naVCDJbQafFYrGgzEDjYLHqnB8jBmcY72nsxMGWLjhtVkwdUSIck4pIFeeuaikvcMPOGwYLOXMszBoZ2QUozo/VauH/fyQiuJiLabFYeKhVL2+uUSvmeOqExpkz0TBYXB9rf7SuVu67d5qJopVcp/Ewa7smXw5QhG6s1iTaBsvs/bOGzvFakwBAaV4kz65DKdoBRGdOHeLVK3oQf8fOfeeNqwAArI4IYD5Gy25N2NjXLKzBcmO7Dy+vUzpXxHLoe1KEwdqT7DnaoaruZlXNw5OoZB0ImBZzd911FzZv3owPPvgAbrdysJ577rl45ZVXUrq4gYhSBaQ+CGLlzHUKwk686LDwTKFmegKD5RTpjeJRnDn58zVT0u3VK4DgIRmWM8cSu538YsWSq7VuiNatYSdh7WseycJq1mwLs6rEnNdYAYQ6zKqfh+YLhnQ734uw70dZvhNOuxVdgRAvJChw2WFPMK7LCHw+a5z1ivlLvFrboDO3KjKH9cRhRch12rmo8sZy5mKEjW1WC3e+VNWsYgGEJpwHKGKP3TCJzgnvNaeTN8f2RWmUM6fO8TLTMFi7vs9qm9AVCKE0z4kxFQVxnqWG9YszclOkPY8AQpg1MmVHRHuDzUKt+yLNjOPlzAFy3jGgOHN8AoS2z1yMalZALe54MUbkbzLncnMk1Klta5VKWOPgf284oDIC6mJ8X5JtTQIAgwvdyHHYEAxLfB8Dyli3ZIofBgKmz3qvvfYaFi1ahNNPP10lFMaNG4fdu3fHeSZhBL1mv+xkI4ZZxaR5tTMXP3coYc4cq2SNnLTyTXTO9uo4c9pk6WbBdSjS9GLSXkzYtnx8TuTEzS5E3hg5c1njzAVYmDU7xJyYy2YkzNodCPFGw4D8Wfl1+omJF4h4n19Du/x3Kz1u3veLuTlFeam5iCWqaBULDApzHELbA2Pfu092yU7K6RH3iR2TYlU3a8Ady5kDgOMq5ao8sYGqOKdUqWZVHBom9liITDzP8CkQumFWeV8ozlx0zpzZhsEAVPNW3/+6HgAw7ZhSU3lRLMzaaeA4Yjd64tgvto8kCapm5kB0yJBFLFiRWEGC98nEL8s59EUVQGibBscRc5qpE5OGFgGQ+681d/iFNJ3U5ssByv5iN83MeY0ZZtXJ/TaC1WpRiiCEilZy5uJjWsw1NDSgvLw86vGOjg5TA3WJ2OgljsYaNyOGomI5c3q5Q4ly5o62yReRsgL5JFRgsDWJJElKs+KoCRDq1iTsbxfnOfjFioWVGtq1Yk7eVrmjll9LFIjy9AfleT0Rc0s3HMDtr2zqkyIKLuay0ZkzIOZqj3ZAkuSLOwsLsouaFlHIxMt94s5xvgujI+4NG4ml/d4lS1mCxsHsWLBYgHynPSqNIBEsFDk2Ish4zlxntDOndxz/8TsT8cr1p+Lk4UpeWY6qabA6zCr+niFW/iZqT8J7zEWEbqGmdyRgvi0JIDdFZ7ptZUTMmQmxAkKY1cAxG6s/HLtZBqJDrdobbO3M53g5c4Cyvxo1zhwPs2p63GlzEGNNhGCRG0+ug9/QbD7Q0iMRlQjt/rrhWyMBRIs5VRFGknl72ry5UFjivTjTmTO3+JNaTLzvv9gcaaGUSZgWcyeffDLeeust/jMTcE899RSmTZuWupUNUPQSR3NihVl9YphVvKOPb3GzC57efFZdZy7BibLdF+Rl6nqtSdhgbhY6LRLDrHrOXIc6Zy5f48y1dgfQ2O6DWCHfEyH2x+VfY+nGg1jySW3Sr2EUJuayIcwaCksqB/GIgWpW1pbkmEH5itvVFtsRFntxGQmzDipw4dhIfg07ufa0kpXBRE6sEVuAkIjussNqtXBhY9SZY583a6lRHCPMmqjVSnGeE1NHlqpuohVnLhiVmwUgKo9KPEfwxsE6YdamDk2YNYYzdzDSwFks0EiExWLh7hxzX047xqSY461JkguzOm1WXojl03YM0IQumdBnJMyZY61uOnQKIKKqWeX9yW5+xDxMrTMHAJOqiwAAm/d7lcrbFBc/AGon85LJQ3BspXwjpb2pE4swEoWg9dBWtB72dsEfCsNhs5j6bqWa/2ypQ2t3EIv74NpgFtN7euHChZg1axa++uorBINBPProo/jyyy+xZs0afPjhh72xxgGFXufsWGFW0ZlraPNBkiRYLBbdIgQGz88xWM3KnLBEFypvnORbdnJhg7nF5O5i7hSqnbmSPCeaOvz8IqJNXGYn19auILf+LRb5byQr5hrbffy1Fn+yB/NPH6lq6ZBqurIozKotSmjzBdHhC6pO8lpYvtzIQXkIhsOoa+2OhEk9Udtqw6zs+6xF/H6OrpBP+uxCmKpcoURutHbIutnu8sxByo3kIRXFcMtbkmi1IrYmYc6ceCyK32VtY3J286Z1xhm61azCe+ajt+JMa4iFy2Hlx8DQ4hzUmAyl5QpVvIng5xFBaMhTdmzoCoSinDmlXZS8fWmUM2cwzNrOCiDUjqnebNbBRW7sb+qK7cwJf3PSUA+WbTyIzQda+Fp6w5lj6S0AcNX04bw1lPam7mjk+Mx12pIuwtCKOTbGq7okN2Gbmt6ENS1e8dURdPlDvXptMItpZ2769OlYvXo1Ojs7ccwxx+Cdd95BRUUF1qxZgylTpvTGGgcUSphVUwDBG0sKOXPCxdUfCvO7Mq9OqJORaD6rOP0BgNA0OP6FSm/6AyCfsNjB39ju54nuxbkOFOcxZ05eT33kTo8lQDPHThmOHTlhuRU3g1n9rOdXdyCMYCh2blY8th1u4//f2OHHq5/vN/0aRpEkiV84UhVm7Q6EMO+ZtXj8g10peT0R5mg4bBYuqBOFWlklq1lnLhiWdGd1is7xaE2SfKqcuYIETps2nGV2iDcTHUxMebhbLlazsiIh4+9J7LcWswBCuLgWaQqk2HvWC1XyuawRp8nDx3kp27N2HTUmQ2HiGs26coAiig2FWWPkzAHqIggRbTJ/VJg1Yc6c0mgYiFUAoW4azP5elYfNn45VABHLmWvRvX6kgjERJ+7sseUYW1nIcyw7/CHVzTMrWog3fi0RLGfuy4OteOy9nXwiTDpDrL5giJ/vOvwhrNxen7a1xMKUmAsEArjmmmuQm5uL5557Dlu3bsVXX32FF154ARMmTOitNQ4o9KqAcpzROXPaxqpMhLV0xn4NBns8VqWeOKbJbJg13jgli8XC72z3Rg52q0UWZNoCCHaxZicP7VQA5e5TaY3A8nxGRu7o5O3NC6SvDssnDXbRe/LDb6KGn6cKUawYaalghA37mvHxzqN4dtWelLyeiCimKyIhuUShVubMHTMoj18E9Zwf7Sgrve+bWG09rCSX3yQA+sUCZtE2udaiTTRPJP60sHAgO65jOXPse+8x8Z6UatZwVJ85QO3SRRdZxQ9VslxHbTVrW7cyUm9vks6cuK7po8zP+FZm0iY+5nmKh8771/aa0+Yxl+bHzgfWo0SYKBIMhXkPTSYelT5zYUiSxL9D7MaUHXeSJPFjQhRrxw0uhMNmQWOHn8+x7Q1nbmxlIT668yw8/oMTAcjnASZkxfYk+5sUFy1ZRg3Kx6ShHvhDYfzvih1Y+J+vAaRXzB1q6Va1RHp9U2aNYDQl5hwOB5YtW9ZbayEQbwJErJw59YVDGYEVPzwzskyxsLUnrqYOP8KSHK5k4RSjBRCKMxf74sNOevsiJeZFuU5YrZboAohImPNY5sxpwqysRF5sjcCcuaHFObzhq5Emx1qYM3ft6cNRmufEwZYuvPVF/KHvySLu+1SFWdmQ86YOH8Ixxuz0BHH/V0YS3OM5c5Ik8Zy5kaIzpyPmmjVhf70iFjHMardZMXKQcoJPWZg1oTOnTocwU80aCku8opeHWXOU5zNHmZ0LzDhzPLfWH+IuvssR3ZpEXDuDrSXWd1GSpOgwa+T5bKSeJClJ6sNMhklFZ256Es4cb01i4Dg6EMnrG1qsXqNbU4jA0BZADBKcObvVoiqeiEUZd+Z8qtdmIlsMs3b6QzzvmOWGsdBqhz/E84LFz87tsOG4yLzRNZF+c72RMwcANaW5KuFdEaOiNRXOnN1mxdIbT8Ojl5+gOr5HDEqfmDsQCbEyN/397fWG0yr6AtNh1ksuuQSvvfZaLyyFAPTn6rEDXp0zF9uZS1QAUV2Sg4pCFwIhCRv3tah+xy6UpXlO3q/LaGsSPSHKYHeT7GBnYlN05rr8IV4JNaYyP/K4JmfOrcmZ6w7ydgoVBW5h7Ix5gbQtcmc7uboY154+AgDwxAe7Tc3cNIoozAMhKSUO4OFIlWRY0i9wSRa2P/Nddh5iiSfmGtp9aPMFYbXIF3elQjR2mLW5w6wzJ18kR5UrbmxP57IyEuXAiT3mxP8aObmLLmwuD7Mqx4w3Rvseo7iZQxUI6oRZlf/X3uzFa7zb5gsq85oj4kRMnWjtCqC5M8CPXbOuDBM2YyoK+OdqBt6axIAzxy7K4kg0eQ3RI70kSVKN8wLUYdZ8zQztWDAnrzsQVqUSxGoaLBY/sIIUdhywa4PDZlF9poDSooQ7dyke5aUHm60tOnOpEHOAnNN58QlDsOJnZ+LRy0/A1dOHY87Eqh69Zk/Y3ySfW6eOLMUxg/LgD4b55I1MwPQnPmrUKPzud7/D6tWrMWXKFOTlqZXyLbfckrLFDTTEkm69EIh4otE6cyzxtCVOs1FADnmePLwEb35xGOv2NGHaMUpY4yif/qCcsPiFLUGYNVF/O/Y6LBTDLlJi02B2oXY7rPyC0BJpPaLbNLgrwCvwygtdyHfZ0dThN9yNn+ELhrAr0tfouKpCnDyiBE98sBvbj7Rh5fZ6nD22wtTrJULrInQFQny8ULIcEk6qje0+7qKkAiXXyMZP4vHCrLvrZVduaHEuXHYbv0gf1Wl90awJs8YS4x2+IL+JYa8nO7iye9rTUV6MhM6cZuyemTAry4+0WJQLut1mRYHbjrbuIFq6AijNdyVMl4hFjqoAIjrMmhMnzBovVMka3uY6bXw7ljpxtF3ub8aE/WCP23TiO9sPyYRY2bqAxDlzkiQJzpxazMWKfnQIThnbX+IxZaRaM9dpg8tuhS8Y5j0XnXYr76MnNg0WZ69qv1Pi9AetgDyhugjPf7qX/9xbzpwWJjiPtIliTn6PPRVzDCbqLj5hSEpeL1lY8UN1cQ4mDvXgkXd34o3NhzD3xKFpXRfD9JXj6aefRlFRET7//HP87W9/w8MPP8z/PfLII72wxIFDm0+/pJuHT4LqE41IgyZnLt6F7ZQRJQCAdZH+XPw1NK4HoIgnNm5Hj0T97djdIrtz04q5lq4Ad9gGFbj445IkuxXaalaxzxw7mVQUulQzBM2w80g7gmEJnhwHqjxueHIc+P4p1QCAf284aOq1jKANcXenoAjisNCgN960hWTgYtrtMBRm/eaoUskKGKiWNBBmZTcbOQ4bD3eMFpy54hQ1DS5MEDbVhln59gZuIHhbEodNdVEW8+ZEMVZsQpAzUSO2OBHDgG6nvpjLjZGXy9CO8mKI81lZiDWZizj7jsw6vtL0cwFBzEXWXt/WjasXf4Y3Nqvzmpo7A3z/a1tcxHLm+Hgsm5X/3mm38n3HKv3jYbFY+M3xocjxKTpr4mxWccKD0kVAfqxV4waLsCIIRm/kzMVCe1MnSVJKcuYyEfa+hhbn4qJJskP48c6jURGFdGHamautzbz+Kv0FdiC7YrT2cAu5MAwWDrFbLQiGJV4lqDdFQuSkYbKY27C3GcFQmIdUtx+Rc8bYQQqo+zG1dwfhyo991x2vAAJQTvyKmFOHWUNhJcdqUL4LDsGtaO70R1WhsRNWmy/Ic8XKC9y8D52RiRUiLMR63OACfpGdOqIUT31cy0fJpBLtRTMVveYOC05ZU4pPMqw1Sb7LJoRZ9Ud6iT3mAGVEll7OHMuZdNgsCISkmGFW8WaDfUZiRasnRU2D2QWzKxBCIBSOckx7UgChFD+oT79FOU7sRxe8XX7ucts17UMSwW76WEgU0HfmtDdduZFQJQv5i+9Z22OOUSC443uTzJcDgN9/ewJunDFKNc3CDLkudZj12VV78MH2BjS0+fiFF1BCrOUFLv1zrHBciqkjovAuzXfC2xUw3EetNF/Ov2ViTvzb7PMJS8p+LnQLzpyPOXPRlayMkWV5KHDZhWa9fRRm1dzUyREReQ1a5zPbYY5udUkOjhmUj+OrCvHloVb8Z2sdrphak+bVJeHMiUiS1Cu5RAOVePPsYpXNszAUO2ga2n0IhsL8gI6XPzSmsgAFbjs6/CGe9B8OS3gzcid7znHKlA+b1cLvfONVtCYKC/HWB6wtSeTC4HbY+EWGiUk29Js5Ac3CSSJf05oEUJwDFmYFzE+B+IqLuUL+GOt3tbexM+XfdW2YNRVi7pDgzOlNWkgWVQGEgTCr2GMOUEL3LZ2BmPmB7A6XVfHF+vxiOcfDSnNR4LbDYVPyjHqK2IMs1k2BduC5mQKIrgDrMaedxqA4c80dSsNgM5N13JrXtFhkccwQxVxUY3LhudrvYqPQ+1GEu+PdAeyNtCUZlkTFodNuTVrIAUqfORY9WLrhAAA2gUQ5bvVCrIB4jlW+m0oLGrU4Yt/lRG1JGGy/HYzcdKqcOeH/mWstzvuNCrPGyIezWi2YWK30buwrZ46dp1lvTnajXlloPtSe6Si5lvI1gd0kaN3fdJGUmHvmmWcwfvx4uN1uuN1ujB8/Hk8//XSq1zbgiFdAkBMrZy7ilLCeTkfbfaqeRPHuzmxWC04aJo8BYqOQPt/XjEPebhS47JgxRj2yzYjzwNwEvVYK2jwO0cFjLt2OiJhjF2smSJs6/EoCfmQtTrtVdXGyWS0ozRPDrOq1Hmzp4lZ5LJgzN04Qc9WRA7etO6g7DSBZfJrmpD2taG33BVWff2+FWfNcdj6uqaHdx3OKtLAec6x6uijHwRt+atviBIU+iSw8E9OZi1zsygUx57BZ8dJ1p+L5+VNTlivksCnfrVjfee1FXiyAYOJhf1Mn/rj8a943kcHDrBrhxW6CWjoDCfNP9dCO63LZrTEnRIh/j+GwWfjno+17qIRZ1WJZzFvtSZi1p4jva/nWOt6qqNMfUqUCaC/IIryaNYYzp91XLGXAsDOXpw6zqp05QcxF1l3odqBAG2ZlIVid0C4rgmDP7wsqNdWsqSp+yDQ6/UF+PmXnp9kTB8Nlt6I4z6F7DuxLTIu5e+65B7feeisuuugivPrqq3j11Vdx0UUX4Wc/+xnuvvvu3ljjgEE7NkbEpZOcCyg9nY62+XioqsBl56FTPU5meXORIeX/t0nOCzt/fGXUXVW+jkASSZSrpz3xiVV6TLTt4M6cfPIr4ZWufn5RzRPCU+Jd6qB8F2xWi9DkWFlrMBTGt//yCS7888cxcxwkSeIOpejM5ThtfC0sjJQqtGHWnjYOFvPlAMVNSRXiTMvSyL7WzsRl+IIhLpyPiThzVqtFd+Ypq7y1WJRmqbEKIGI5cwAwYagHp45MLnlejwLBddKiLYBgNxiBkNLs+IkPd+OJD3bj5XXqxtPsfek7c/6kGgYDsggVnTgxxArE7zNnsViUSQqaitZmLubUz1GcuaAwCL3v20e47FY+3/W51XtUv2OFOEACZ84efcOsNx6LVagmGuXFYN/7g1zMKedmi8XCq4KPxnDmugNhBISbHb1KVTFvLtkxWmZhDn19m3xTx455sxM8Mh32vSlw2/lxM7Q4Fxt/fR4e/8GUtE6lYJgWc0888QSeeuopLFy4EHPmzMGcOXOwcOFC/O1vf8Nf//rX3ljjgEFpSxJ9IMYa58Uqt1iOytF2P+/VpVfJKnLKcKUIIhAK835qF58QXf6dHzmZxctDa02UM6e5eIhijiWuM7ueXayV6RCBqKbBgPqCxBrZMrEnhumaIpWyrd1BvB7DFj/k7Ya3KwC71cJHRDHY/t0Xx9VLhugwa88aBx/ShDxjNYXuCaygJM9lh81q4SK3LkaodX9TJ8KSLPxE4aXXOLhZyBViFyvt+DBAmRs6KD814dR4xHOj2Y0XC6/mO+18tifb/stI13pt7qISZo3OmQNkYaukLJjPAdRzfQBNzlwMB12vopV9Ftp8LSZyDnu7uRBJx4XcYrHw435DpN3SiEjYlhXiAPo95gAhzCrcMOs5c2eNLUdxrgPfGj3I0PpYmFUpgIh2UAFlOkphjiMq1N8aJ2cOAE6sKYbTZo2ZD9hblOU7YbVEms13+PqtM3eAV7Kq35f2GE4npsVcKBTCSSedFPX4lClTEAympov9QIXf7cfMmYu+a2TOHDtw/KEw/9IZaWcwYagHTrsVjR1+PL9mL5o7AyjLd2JaDIeDNw6OU63XkqDwQnt3K7oO2gsLy31igq+hzcfFjziGR3zNQZH8jXxN4jCgDjn+6/MDUWvbdkgOsY4qz4860TJbPdViTjsDsqdhVubMMVGR6pw5pTWMvH/i9Zo71KI0cRbDfFzMadqTsJuQkjynbpgcEEZ5JdGLzCzxes3xGZkR4Wm1WpDvVEKtobDE8z+1uX/a6Q8MVc4c7zFnPlwmCjaXQ1/MxTpOxXFgInpuIjtXfXnIy99Dqho3m0XcnxOGeDBznNxKiBXiAPo95oDY51i9nLmzxpRjwz3nYabB6ls2z5V99tpGw+xzYoK40O2Aw2bl27X7gsI84NgCYlCBC//6yTS89KNTDa0pFdhtVn5MH/H2XzHHesxVl2RuUYdpMffDH/4QTzzxRNTjf/vb3/CDH/wgJYsaqMSrQs3h+RzRI6BK8pz8AGd90ozk2rjsNpwQybN4+N0dAIDZE6tihmd5mFXHmfMHw/xEpT8BQhNmFZKptRetQfnqAgh2EgbUA58LYzhzsQogRHdky0Evvq5rVf29bTGKHxjsxLQvw8OszJkbGXEkUu3MMXeGia1YDUMZLIemXKiKBqA7BaJZGCofr4BFL8zaG+g5c8FQmN9IiTcT4vZ7Gju4WNc6jF06OXNFQoueRJXh8RBFjVtzYxIvZ07+fezmu9oejwwWRfg6kqJgdoxXKhH35/dOGsoLb9iw9ng95gAhlSWYOGcOgKnCFG0VcLQzJ/8sFkDI/2WN0QNRRTexmDi0SNVEuy8QK1oV0dPfxFxsZy6TSMojfOaZZ/DOO+/g1FPlO4BPP/0U+/fvx5VXXonbb7+db/fQQw+lZpUDBO0MQBF21+gPhREKS7BZLfxuOS8SymrrDipizmB45uQRxfhsTxM/UcyJEWIFEFUmr7d2i0X/zlHrOKoLINTr5WHWyOPsJCH3exLEnPC3mFMUS3hqXap/rT+Au2eP4z9/JbQl0dJ3YdbUOHMThniwu6FDtwVIsmhbw8TrNccS0Cs0oqusIJIz16YWmizMWpzrzBgxV6jjzIniTvyuF7gdgLcbbd1B3mAUiO53qFcAwXJNvZ1+vj+SmWgRz5lzq8KscZw5nTBrrlbMRdbsj1Qn16RxdiYTok67FXMmVWFn5FzInLl4PeaA2GFWbW5ksmjnuUY5czzMGnHm+Jg4OxrafGjrDqp60GUS8nnXi/1NnTjkTW3D4ExhfxxHN1Mw/a3YunUrTjxRHrS7e/duAMCgQYMwaNAgbN26lW9n5q6FkGEJrvFakwCyo5PnsnNnLs9pR1m+C980dHAxZ7Q0/eThJQDkz7G6JAeTNc0nGYlGerF5sJ4cB+9srsVIAQQgC0J28mOOHTuYtAnHMXPmYoTpmEtVnOtAc2cAr206iF9cMJb30lIqWZXyfkZNL4VZfTEmQPQE1mNu/BAPXtt0CK3dQfiDYT6rtqfwnMXI/q2I056EOXMVGmdukK4zxxL+9cOsYaHYIp3OHPs512lTudhiRStzqgCgXSMGeZjVocmZY2HWrgBPWeipM6d1gJjQc8boZQmIYVb1e1ZGuWnCrBqRMzyNie+sH9/McRUoynVyh/qQtwvdgVDcHnOATgGEgZ6dRijV5HhqHVN2jDLHlztzwo1pW5w0nHTCHPoN+5ohSfJ3rCzf/E1IJpMNjqNpMbdy5UpD2x04cADhcBhWa2ouJAMBr05+BqA++LsDIeQ4bModvsvGL5J7+BB7Ywf8icOKYbHIUxbmTKrSFeE8Z05HzBkZPSSehApcdlVTUjHMWpLr5L9joVheyaq9mAivWa7JmWuPEWa9YMJgvPNlHY62+/HB9gacN64CHb4g9kaEWixnrqZEvigc9nalVBxFjfPqcZhVPuEcN7iQN5Ju6vBzB62niE4wAFR6IgUQMZw5RcypL2KDdBoHizlierN1W7oCCIbV80F7Ez03Ws+tEcUfuzkAot9Hl18RgyJizlxLEnNZGTlxCiDkGbmumN9z8bl6BRDahG/tuSqdjszkmiJ8cdDLZyqX5DnhyXHA2xVA7dGOuCFWQGhNouozp570kSxRYdaonLnYIlkZpRjgN/uZ58zJxyKbJlRTktvvzBxeAJHBYq7XlNa4ceOwZ8+e3nr5fklrnDCr1WrhIqI7GFYJAdmZk08WrPO70RmVhW4HZhw7CLlOG747pVp3u1gCScTICDGx4q9I0+JAvGiJrou2Sal2fI64r8rj5MyxPlkVBW5cMlme8fevz/fj67pWXPnsZ5Ak+Q5TewcNyBVbOQ4bwpLSWiAVaAsgehJmlSSJT8GoKsrh+y2VoVZtmDVeAUQ9D4fGzpmLKoBgYdY4BRDsOcW5jpQJ6njoFUDoJcWL+U2imNO+D70CCFa52tod4N/XZGbNxhNzeS47PvnlWXjumlNiP1cnzNrJhLxWzGnOVck0DE4V/9//Ow4b7jkPJ9bI/TMtFgvPm/umoSNujzlADLMay5kzg1sYPwfoV7MymGDLF26i21IU8k013KGPVJpnsuBJBm+XIqQzOczaa2dEmgxhnngTIADAHTngu/whfqdsscgnIW3YyczJ54kfTsEnvzibl/LHIlGHe37Si+MkWK0WfnLSOg5iMYT4XrQOY1SYR7igspMKu+ConblIB/t8Jxet726rx4V/XoXP9zYjz2nDvRcpOXQiFouFOw6pHOvFLhp8OHoPwqzergB//mCPIkobUzTSS5KkqCT4eFMgWAsRrTOnXwAhhlljTxvpy3w5QOwzZ86ZO9DcpWoToy2AYPND9Zw5SQIORMI6yeTMueOIBvaYXipErl5rEl9iZxxIbpRXqrBYLFEFGqxh9TcN7caduV7ImQPUodYoZ04j5jxCzhzAxFziAoh0oHX++12+XCRqU5rnzKhWJFooBppBxJsAAahL58U7ZXGQM8NMro3bYUs4zFtpGhy7NQnP8UkgItmJKErMCestF9wc7XbakzU76dmtFpREto3VNJjlzJXmOTGmsgATh3oQCksIhSXMOr4S795xJi6YMFh33ax3VrwJEmZh4ou5aF096DPHWoGU5DnhFnJWxMbBXf4Q/vbR7qQEqS8Y5iFOdkFnJ/EOf0jlXoXDEurbYufMsXU1a0Z6sTBrSZ5+NWtDu/yafSfmYt/A6IW7mHu9fq8cbmIXaO370Ktmddis/L2zgoKkcuYEZ06baJ/4udHVrJIkRVUyM0Rh4XZYVZM5MgHuzKnCrAmcOYPVrGYRowzanDmt6NaGWVs6/fwzybwwq1bMZa57lQzc0c1wkUpiLoPg+Rk6d14sBOILhoQcFvkxrZhLdZJsojCrt1MpgIgHOxFpW5EU6YRZHTar6uSldzEZVODibgPbpjsQRjByUdQOCr9vzvE4//gKPHPVSfjrvCkY7Il/AuqNIggWZmUNk3vizB2O5MuxMVvsfYrtSZZuPID73/4af3pnh+nXFz935nzmOu38cxTDz82dfh7u1wqv4lwn75YutotRWpMoYdZOfwhhYUwOd+b6oGEwoC5oEGnVueli38WvIj0LT4gUE8lTIcSZyrHzz4Do4yeZnLncBM6ckeeKNxa+YBjsY9AKULdDmTgxrCQv43KljuFh1nYh7yn2se7SFEAEQkq7Jb2pC2YQiwK0BRhap46db9l/Rac348VcP53+UJ3BIVYgzWLuo48+wkUXXYSqKjnx/rXXXlP9XpIk3HfffaiqqkJOTg5mzJiBL7/8UrWNz+fDzTffjLKyMuTl5WHOnDk4cEDdFLa5uRnz5s2Dx+OBx+PBvHnz0NLS0svvzhx+IQ9OP8zKTrTKSYZd+LQXTaOtSYySqADCaF8sdgHUho8K3XZ+kdfe3Yt3tFpnbnJNMU4fVYb5kaRnQB0KYsnnLNzIqmRPrCnGk/NOwjnHVcRdL4OFj1I50otd4NkFuyc5c+xkz0QpC+kcFVqysDyuZASpIkDUIbpYPfhY7kxZvlNV5ALIoXYmNMW8uVitSQB1iDJdYVY9Z04vzMqEz0nDi/nvxCrwLp0wK6A+fuSKU/On6HitSRI+N0bTYFHIawWoxWLh+yETL+IjB7EwqwlnLnKTxUQ7oD91wQzieUwbVhV/zncp50LWeolNjtBWUGcChW676jvXX8OsmZ4L2GvfCiN3aB0dHZg0aRIWLVoU8/cPPvggHnroISxatAjr1q1DZWUlzjvvPLS1KWX/t912G5YtW4aXX34Zq1atQnt7O2bPno1QSDkZXXHFFdi0aROWL1+O5cuXY9OmTZg3b17P32QKEec/6s37ExN0xYsrAJRpxVwS4Zl4JGpNkmj6A6OQO3NqMWexWHiINkqY5uqLuRynDS9cNxXXnTGSP+ay2/isw3Z/EIFQmItN7aBwo/TGFAgWblPCrD1w5iIn+6qiiDOXH+3MsbY12sHvRtAWPzBi7RcWYtUWPzC0eXOhsMQ/n+I8B1x2K7+YiZWgfS3m9PrM6RdAqH8+vsrDj0/xfegVQADa3ouOpJyueOO8EhErZ46ldOQ4bDFnULIbtHQ2DNZjWGkurBa5IlnpMRf7e6l15ng4XRBXPUHMmYty5gQHVd27UC3mMs2VA+Rzt5g3pyeWs5X9CXItM4Ve+2YYKYC44IILcMEFF+g+/5FHHsGvfvUrzJ07FwDw3HPPoaKiAi+99BJuuOEGeL1ePPPMM3j++edx7rnnAgBeeOEFVFdX491338X555+Pbdu2Yfny5fj0008xdepUAMBTTz2FadOmYfv27RgzZkyK3nHPEBtC6p043EKHckswElKM3ClrS99TPVJHKZFPUM2aICw0trIQ726rj9kaYWhxDho7/FFJ1CXCBU4rJvTId9vR1OFHe3cQjsj+tFqSqw4ElAvVvqZOSJKUknBSt8aZ61mYVe3MlUVEqxjK3B1pnsqGYpu5QCl9xmK3ohBzCfWKHxhlBS7gsCLOWrsC3M0qynFG5mza0NodVLlCfTnKC0jcZ07r1hRoKq2PG1wY6QcZUr0PJWcu+rssOurJuuvx+swlIlbTYL18OQa7QUtn8YMeLrsNQ4tz+c1GRaFLd59ox3klymE2S6lBZ050fFn1PqsOz7TiB0Z5gQu1RztQUdh3c2H7imyY/gD0ojP31VdfYdiwYUk/v7a2FnV1dZg5cyZ/zOVy4cwzz8Tq1asBAJ9//jkCgYBqm6qqKowfP55vs2bNGng8Hi7kAODUU0+Fx+Ph28TC5/OhtbVV9a83iTf9gaGcbMJCR3Yb/x07qTpslpghnJ7AE7ODYVX+D8OoM3f7ecfiozvPijnT8H8vnYTHvj8ZE4aoG/eKxRlG70zFikgWYi3Jc+pW8SViSHEOLBbZsUhVhSgL56TCmTuo68zJFwFvV4CLp1BYUhVGGEFvnFMsZ473mNN15ljbFHk/sny5ApedtxyJVQTBRCIb9dbbMLHW6Q/x3EsgcTUrILtYNSW5yvsQwsWdOn3mAMCjmlec3IVbnTNn7hTvjtFnTq+SlTF9VBlynTZMH1Vmdql9AsubA+K7RkoBhDrMmjIxZzBnLpYzx7yRTHTmAKUYqr+FWMURcJkeZjX0zWDOmBGWLl0KAKiu1u9ZZoS6ujoAQEWFOqepoqICe/fu5ds4nU4UFxdHbcOeX1dXh/Ly8qjXLy8v59vEYuHChfjNb37To/dgBp6HE+fEwU42XYEQQpGLi9j3qazAhdbuIDwRdyOViBfx9u4gXPnKyWjz/hZ8FRm0PThBg1qr1aKbWzOqvACjyqMdOzEka9iZczkAdKHdF+R32tqedWZw2W0YXOjGIW839jV1RhWcJAMTb0ys9iRnTimAkJ05pc+cLJRYiJVxpNUXNTc1Hu06F/RYhSFHeCVr7H2knQLB2pKIvQfzYog5JhJZs+LeRrxwtvuC3HVO1GcOAMZUFsBmtcRss6JMgIgWR8VxRtwZxd2DnDnmFnYGRGdO30kEgF/MGovbzzs2Kj8yUxg5KB8rtzcAiB8q4yMTg2GEwxKv+tZGPZJFbHQdPc4r9rxprXhLRe5eb8DEXKYLHrM0dvjRFQjBYtEPz2cKhq6MHk/0iKO+QitKjIS4tNvE2j7R69x1112qObOtra09FqjxUMbG6H8kObwPUoh3KRfvwtlIr3ivkSw2q+z2sZARy/9o7Q7gpn9sQCAk4YLxlTi+KnpQfU+JVwChB+tH1+EL8hYYPRFzgJzgfcjbjX2NnbwxaU9gDmdJD8Os4bDEh90zMV3G+8z5IEkSdmvEXF1rNybA+HGt58zxMGtzF8JhCVarhRdA6IlFFiZloSNW/FASQ7QzEdTlD/EbHjMitCc4bHIBQncgjLZuRcwxEVqiEVvihfe4wfJxwHseRtYeCksxj12GKsyapDOnbk3S82pW5bPXf61MFXKA0p4EMCbmALmC9+0t8s3+GaNT4ziKzly8psGFMZw5/rsMG+XFmDt5KL461IofTE0+GpeJWAD8ZMYxaO0KmE5Z6GsMXRkXL17c2+uIorJSDsPV1dVh8GCl/1d9fT136yorK+H3+9Hc3Kxy5+rr6zF9+nS+zZEjR6Jev6GhIcr1E3G5XHC5+q5nkrkwa4hf+EWnil0kk2k0aoQCt5z/w3KGJEnCL//9BfY3daG6JAcPfGdir7QmKI5TAKEH74vXHeQhrljTHcxQU5KLT79pSlkRhOLMOVQ/m+Vohw+BkASLRblDZheO7oBc+byrIVrMmUGvAGKwxw2b1QJ/MIz6Nh8qPW4u0rQtCxjHVsju66qdDfAFQ6q2JAxteJIVVeQ4bLyyui8ocDvQHfDx0KokSUp+omZYu3jsspxQbbhYFOwxW5OowqxJ5sz1oAAiJ0YBRLxWKtkAaxwMJAizCvvqQHMn1tY2AgD+X5z+k2aI78yJYVZHzP+Xf87Mz2BMZQGenz818YZZRmm+C7+YNTbdyzBExt5OjRgxApWVlVixYgV/zO/348MPP+RCbcqUKXA4HKptDh8+jK1bt/Jtpk2bBq/Xi88++4xvs3btWni9Xr5NJpBo+gOgiLmuQIgnpIt39yx8leriB0a+xi15Ye0+vL2lDg6bBYu+f2Kv/V0x9KRX6atFdHa0PeaShY0qSlV7EpabwxxDfyisys0yChvjNSjfxR2SXKfSLqCx3c/DrOyiEWtqQzz0xJzdZsWQiKhhIpdVy+o1kD1tVBkqC91o7gzgnS+PqOayMpTwpPw9PyIUVfRlLzNtEYRYFVmpEatiCJo5c9r+jJ2ayS1axAKdpJ25FBdAsPdr9EYq0zjGoDNnt1lhj+TULtt4EGFJ7hWYqtChqmlwVM6cGGZV9rN2n2dqAQSRfpI6Ov/1r3/hn//8J/bt2we/X50MvmHDBsOv097ejl27dvGfa2trsWnTJpSUlKCmpga33XYb7r//fowePRqjR4/G/fffj9zcXFxxxRUA5PDv/Pnzcccdd6C0tBQlJSVYsGABJkyYwKtbjzvuOMyaNQs/+tGP8OSTTwIArr/+esyePTtjKlkBYEpNMW46a1TcMKVL6IPUGaPCrFLTMDbV5EdOJEdau/Hwih144oPdAOScmUmRBqm9gVgAoZ0NqYcoPFneWE/DrNUxKjeTJRAKIxQp4RTDdV2BEApMhqx4vpzGKSrNd+JAcxeOdviwO+LMnTy8BKt2HeX5Z0bRC7MCsmO5r6kT+5o6cdKw4oTOnM1qwaUnDcWf39+FV9btx4ShcrhX9TlrHC223r4KsTK0UyCYcC7KdUS1FrHbrDihugh13m5+HCvvQxZEvJLVYYspSkU3TttY2yg5PSiAyI1MgIjVZy7VRVV9xaACF8rynTja7ud95/RwO2xo9wWxdMNBAMDsialx5QC5b2BpnhONHf4oh03Pmct1yu1g2LkiU505Iv2Y/mb8+c9/xq9+9StcddVV+L//+z9cc8012L17N9atW4ef/vSnpl5r/fr1OOuss/jPLEftqquuwpIlS/Dzn/8cXV1duPHGG9Hc3IypU6finXfeQUGBkiT/8MMPw26349JLL0VXVxfOOeccLFmyBDabcuJ58cUXccstt/Cq1zlz5uj2tksXU0eWYurI0rjb5AhhViUpWXmf3zlxKOpbffj+Kb2T28fCW7f/czM/uVw4YbCqYW9vUJJENasY3mJzWVMRZgVS02tOvFgW5jhgscgVa13+kOkkZzbKa4gmQbc03yXPCW3p4gL0tFFlWLXraBJh1tiD1gF1RWtjhx+hsBzyFTvea/neSdV4bOUurNp1lD8mhtMLdMScnkDsLQq5Myc750w4a105xqs/noawJHFHTBsuZqIuR+emRHTjPMm2JklR02CWVxzrxjGbsFgseOrKk1Df5uMush5uhxXtPiUNIVUhVsbCuRNQe7SDu/wMvdYkbN6skoaTnZ8B0fuY/mY8/vjj+Nvf/obvf//7eO655/Dzn/8cI0eOxK9//Ws0NTWZeq0ZM2bE7UdnsVhw33334b777tPdxu1247HHHsNjjz2mu01JSQleeOEFU2vLRMTWJJ0xwl6DClz4tc6w+FTALkyhsISaklz88oKxuGB8Za+HvYqS6DPXK2HWiGipa+1GdyDUo35KrMLWYpFP5DkOubgkmSIIbSUroyzyftfvaUZYki8E44fIjlGyzlys9hRirzn2umX5rrid6qtLcnH6qDJ8vPMoF3SxnDnmCnG3r49nf2rDrKzQpEpHFGgLAZj4Zc/vCsR3ubRNg5NBnTOXXJhVkuTzTI7TxgWoXmuSbGCywYIlcX9NGVas+zknS6yWTNq/GzXzVxRzGVoAQaQf0zlz+/bt47lmOTk5fBrDvHnz8I9//CO1qyNUsARd0ZkzGnZMBd+eXIUxFQW4Z/Y4rLj9W/h/Ewb3Sf5SaZ4LNSW5qC7JMZyXVyDkKol95npCUa6Dz6Fs6mGvOV+kx1xOJNwWq/O+UfZEcvi0+UCsCOLTb+RE7lHl+dxRYjloRmFiLpYzKjqWDW1KblsiLj+5RvVzcQzRnm5njjUCVpw51h7F2DryhKpqQPl89cScJycFBRAp6DMHKPl92V4AYQYxjzGVIdZEiA6qVrDF6jtHEFpMfzMqKyvR2NiIYcOGYdiwYfj0008xadIk1NbWGpr6QCQPO0l3B0JK49E+vFueNX4wZo3vuxMcw2a14J2ffQuSBMNTC0QxwEZaxQv7GcFisaAo14mGNh+aO/09umtnDhy7eMaqIjTKjiPyDRWrEmWwsPL2yO+PGZTPc868XQFT7qJeAQSgDC7fJzhz5ToNg0XOHVeOkjwnF8YlMapZ26Ny5tLrzHEX1KCo1FazxhvlBcgOzbDSXDS0+ZLuayW+tln32Ga1wGW3wheUq6BLhTXnZWnOnBnY/rJYgAv68Fyn15pE/tkR8/8JQsS0M3f22WfjjTfeAADMnz8fP/vZz3DeeefhsssuwyWXXJLyBRIK4jivzjg5TP0Rt8OmewGMBbuItnQGejyXVYS5R80dgQRbxoeFWZnbKuZDmqHLH+I5fFFiLuJEsnusUeX5qqHYdSYqWuONdGLOXEObj7uERpw5l92GuZOH8J+LYvaZk/cHm/5gRCSmEpa/2MrFnDlnTlvN2pXAmQOApT+Zjnd+9q2kG8T2pDWJuDZ2wxFPyPc32Dn25GElhj/jVKAOs6o/93yVM0dijoiN6aPzb3/7G8JhOUT04x//GCUlJVi1ahUuuugi/PjHP075AgkFdsB3+UPKOK8BcLecDEzM7W+WxUVP5rKKsCR91k4jWdgoL7eTOXORzvsmnbndDe2QJFlkap1H7ZSKUeX5fCh27dEOHGntxvAydSK2HnqzWQE5NFjgtqOtO4gNe5sBGBddl51cjadX1QJQN1XN14Qn602Eb1NJgaYAok4zAzcReZoCCGX6g/6pt6eFOo5Ii41gWEqq0Wmu047mzgBfa7YXQJiBpWJcNKlvIxAqZy5qsgiFWYnEmP5mHDhwQDUJ4dJLL8Wll14KSZKwf/9+1NTUxHk20RPcYmsSnpRMB3cs2N0syw0rzk1+LqsIE3MtCcRcOCxBgn5YmIdZIxfbXD4TMxhzez1YiHV0RUFU/mKpRtyNKpfbMrCh2GYqWlmYMZYTbLFYUFOSiy8PtWLTgRYAxnPbRlcU4K4LxqLTH1I9RwyTt0f+AeloTaINs7KGwWbDrGph1Ns3YZedXI19TZ0YEqevmh5KyF9dgTsQzjW/mDUGpx1Tiu+f0rfXMb1qVkB9A0UFEIQepo/OESNG4PDhw1HzTpuamjBixAiEQsnPlyTikyOEWXnYi5y5mGhFh1bYJAub1sDmiepx19IteGvLYbzzs2/FzK3jYdaIQBfzIc2wnefLRffPEjvOO+1W3v2ehY/MVLTG6zMHgIs5f6QRshkH7YYzj4l6TKxmZU2I8132Pm9cq/SZC6CtO8BFpV5rEi1R1awGwqyp4A+XTEj6ubma7+JAOtfozYfubdjx77Jbo0Lj7DtotQyMz4BIDtNnRr2Zpu3t7XC7M3sQbbbD8jm8nQFE2rwhdwDcLSeD9qLf00pWRpGBMGsoLOGNLw6h0x/Cql1HcelJ0X3/2IUyh4dZkyuA2HlEbgY8piL6AiQK2JFledwlZEKkzmusojUUloTxcbEvJjWaLvk9zW0TCweUWa99G2IFxD5zQR5iLXTbDbtUUQUQgfhD6zOBHIf6u6hMm8ncNWc7Q4pycPX04RhanBN1fWXucIHb0afTT4jswvDRyRr6WiwW3HPPPcjNVU7eoVAIa9euxQknnJDyBRIKTMw1CUIipwe9zvoz2rFfpSkofgCUAoiWOM7cNw3t/EK4M+KcaenWDbOaE3NimDV6rYqYO6Zcce5YqPJImzFnrkMI/eqJGO3Io57mtomTE9hc1oo+Ln4A1BMgDpvMlwMU8dsVCMmiuI+cuZ6gvbFg4dZsHeeVDVgsFtw35/iYv1PEHO1/Qh/D346NGzcCkJ25LVu2wOlULhROpxOTJk3CggULUr9CgsOEG6tOzHHYDLfqGGhoHaRUhVmZMxevz9wXB7z8/3fWt8fchhdAaFqTmAmzdviCONAst8rQVrICcmjVk+OAtyuAUcIYI95rzmA1K3OV7JG2FbEQnTmrpedJ/PkRF8gfCvP32NfFD4C6AEIZm2ZcVIo3FR3+IBdGZiqz+xpxPms4LCm98bK4aXA2IzpzBKGHYTG3cuVKAMA111yDRx99FIWF+jNEid5BO5g7mzuy9zYuuw1OmxX+kHqYfU8xUgDxRaQIAFDCoFpS0WeOCcWyfJfu+yvNd8piTnDmKj2yKDJaAMHz5dx23TCP6MwNKnD1+CZD/G6zubJ93TAYUC6kHf4QDjazSRvG1+Gy2+CwWRAISejwBWOO4cs0WKVtpz/Ew8LAwGmDlGmMjuTwjYmRF0sQDNNH5+LFi/n/HzhwABaLBUOGDInzDCJVuDQhVcphiU+eywZ/pyzmejrKi1FioADii4OKM3ewpQvtvmBUiEpbAJErXECNsiNO8QPj6unD8dYXh3HmmEH8MZbPVt/q082BFYk3l5UxpCiHz5dNheiy26y8ee03DR3yutMi5hQ3hInnykJzFaJ5LjtaOgNo7w5mRZhVceaCfGyg1RJ9M0n0DeOHePDxz89Ky80MkT2YPjrD4TB++9vfwuPxYNiwYaipqUFRURF+97vf8f5zRO+gPZlm8gUhExBDXD0N+zESFUAEQmF8dagVAPjor10xQq3dwjgvAMhxKqPajLJTZ/KDyJXThuOVG6ap2h2wi4I/FE5YlQskrmQF5JBuVSSXrDxF81PZ3/uGO3N9H2Z1CtWFTDybceYARQS3+8Qwa+beiImj5XjDYKe+K0v0PtUluXAm0QCaGDiY/nb86le/wqJFi/DAAw9g48aN2LBhA+6//3489thjuOeee3pjjUQEp80KMXo1EPo+9QTRSUp1mLWtO4hgKPrmZeeRdviCYRS47ThpWAkARQSIdEeFWZkzZ7zP3I5ICDeemIuF027lTqWRKRDKBID4Nw9srFeqHDT2/WbTF9LlTDB3jk23MJMzB6h7zXFnLoMLl3KECRB8lBedawgiozEt5p577jk8/fTT+MlPfoKJEydi0qRJuPHGG/HUU09hyZIlvbBEgmGxWFSzFsmZi49Y/ZWqMKsnxwFmULR0RbtaLF9uwhAPxlTKIitWRWtUnzkHu4Aad7eNhFn1MFPR2s4aBie4oB8TKbIYmkSj2lho/146qlkBpT1JKNIPyLQzFxHBsjOXTWFWxZmj4geCyGxM3241NTVh7NixUY+PHTsWTU1NKVkUoY/bYRMGX9PdcjxEMZAqZ85mtaDQLVeINnf4o0ZmsXy5CUM9qI406Y1V0ap15sQ8JSO0dgd4q4xYbUkSUVnowrbDxipaOwy2pvjpWaNQVZSDK1LUPb9A8/fS0WcOiG4JUWmiNQkA5EecvQ5BzGVyNavYZ46P8qJzDUFkNKaduUmTJmHRokVRjy9atAiTJk1KyaIIfcS+chT6iA8TH1aLeoh7T2G95mLlm22JtCWZNLSIhz9jVbT2tJqVvWZFoQueJEb8sCkQRipajQ5aryrKwU/PGpWyfS2GdT05DpUr3ZeIRRAFSUyh4HNm/UH+uWdy8RIP+QdCwiivzBWfBEEk4cw9+OCDuPDCC/Huu+9i2rRpsFgsWL16Nfbv34+33367N9ZICLiEIgg6wcaHXXSLc50p7cdXnOfEnsbOqCIIXzCEr+vk4ocJQzz878eqaI3qM+dQ8pSMYKT4IR4s/8zISC8jBRC9gSge01H8wBCdObP5coB6pFdfzWbtCaJL3OEjZ44gsgHTztyZZ56JHTt24JJLLkFLSwuampowd+5cbN++HWeccUZvrJEQYBMDgMy+u88EmBhIVYiVoddr7uvDbQiEJBTnOjC0OAfFeU4ehtVWtDLRlhMVZjUm5ranTMwlHumVLncmXyXm0teWQRRzZkOsgPI9bOsOchGfyWJOdIk7qACCILIC00fovn37UF1djT/84Q8xf1dTk5p8GSI2YnsSGrocHyYGUjX9gVGkE2ZV8uWKeBuHYyvycbTdhx1H2nBCdRHf1qftM2cgzCr2hNvJK1mTaySqzGc1HmbNd/VtB3pRQPR01mtPEMOsg5MQlex7eLRdEc6ZfCPGKm27/CHeZ46iAASR2Zg+o4wYMQKHDx9GeXm56vHGxkaMGDECoZC52ZKEOcTE6Vy6W44LE12DUiwEinV6zW2JVLJOGurhjx1bUYDVuxujKlq1YVZ3gjDrFwda8IOn1iIQDqMk14mj7X7++smQXJi1by/omRhmrTRZyQoo76OhTRZzlgxvwJvL2+SE0M7DwnSuIYhMxvQRqtcxvr29HW43dajubcQwKzlz8Zk9sQo7jrTjB1NT6xbzAgjNfFY2k3XCEEXMsTFa2opWbQEEu1j6g2GEwlJUjt/rmw6hLSKqDkXctAKXvQdiThZHjR1++IPhuA1JjRZApBpRPKY3zCo4c0mIOfY+6iNiLsdhy+gGvKo+cz4KsxJENmD4CL399tsByL3O7rnnHuTmKrMYQ6EQ1q5dixNOOCHlCyTUqPrM0Qk2LoMKXFg4d0LKX1eZAqGEWbv8IS7YJg4t4o/rVbRGjfMShHmnPxg1VPvzfc0AgLsvPA4nDS9Bc4cfI8rykr7IluQ5ke+yo90XxM76Nhxf5dHdNl1iLhOducFF5nPm2CQS5sxlcr4coM7fVAogMnvNBDHQMXx23rhxIwDZmduyZQucTiUPyel0YtKkSViwYEHqV0ioEMUcnWDTQ6wCiG11rQiFJZTlu1TCY3TEmdNWtHZrCiBcdiufbdoVCKnEXHcghK2RfLyZ4ypRU6rcSCWLxWLBicOK8dGOBqyrbdIVcy2dfj6eLFXNgI0iFkCkYy4ro1AUc8mEWSOua1NHxJnL8OOWfSf9oTBau+UbFrpxJIjMxvARunLlSgDANddcg0cffRSFhYW9tihCHzHXhvJY0kNxXnQBBBsGP6YyXxVCYxWtR9t92FXfzosgtDlzFosFOZGG0NqK1i0HvQiEJAwqcPGRWanglOERMbenGVefNiLmNv/ecBC+YBjHDS7EuMF9e8yL7TBSNe81GURhnUzOHBOlkQESyHVk9nEris2GSG5mX+dLEgRhDtNZuIsXLyYhl0ZUzhydYNNCLGeODYMfUZYXtT2rOGXjt0JhCf6QWswBQnhLUwTx+V45xDqlpjiluVYnD5dnx362pwmSJEX9XpIkvLh2LwDgB1Nr+jzPSwyzDkqjmCuMiLk8py1qKoURtOHpTHfmXHZlBvRRHhrObAFKEAOdzC2pImKSo7r40wk2HRQLOXNMBNUelZ25kWXRrUKUvDlZzPmCilgTP0+3MEZJZP2eiJgbVpyS9TMmVRfBabOioc2HvZEh8iJra5vwTUMHcp02XHxCVUr/thHYZIuyfCdc9vQJoLGDC/CtYwdh/ukjkhK0WjGX6TdhFouFn1saIu1U+rphNEEQ5qAjNMtw0wSItMNanoTCElq7g/DkOHiYdcSgaGfumEje3O7INmIY1WUXw+bRjYMlScKGSPHDlOGpFXNuhw0Thnrw+d5mrNvThOEaV/HFtfsAABefMCSqIKMvOG5wAX50xgiMH6JfnNEXOGxW/P3aU5J+vlYI5WR4mBWQ3cN2XxD+YOY3OSYIgpy5rMNNzlzacTts3FFr6fQjHJZQ2ygLtWNiOHPDSuSChf1NsvvVHblAOu1WWIUWJGwmpijm9jR2oqnDD6fdiuOrUp/ewEKt6/Y0qR4/2u7D8q2HASDlrV2MYrFY8KsLx+HiE4ak5e+niny3+jjNBmGkXSO1JiGIzIbEXJbhipFjRfQ9xcIUiIMtXXKvNpsVQ2JUfFYzMdfcCUmSuFjL0QyOz4m4rp1Cztz6iMiaOMTTK6HGU0bIbt+6SCiX8a/PDyAQkjBpqCftzli2k6v5nLPhuNV+N0nMEURmQ2Iuy2AnWafdCoeNPr50USRMgWD5csNKc6Oa/QLAkKIcWCxyBWtDuy+qxxwjlztzQf5Yb4VYGVOGlcBikXP+6tvkZsThsISXIiHWH0wd1it/dyBhtVpUbYQyvQACiF4jtUEiiMyG1ECWwQQAnVzTS0meUtEar5IVkIV3VWRA+/6mLl4A4Y5y5qJz5sRK1t7Ak+PAmEiBBiu0+Me6fdjX1IkClx2zJw3ulb870BCdrWxw5rRrpJQOgshsSMxlGWycF4U90gsrgmjqCCiVrIP0h96zhrv7mzrR5Zdz5qLCrJELKAuzejsD2BGZHHFiiitZRU4ZEWlRUtuEPUc78Ps3twEAbj13NF3EU0S+Ssxl/j4VizScNmvccW8EQaQfOkKzjOKII1Sa50ywJdGbiL3mvuFtSWI7cwBQIxRBsDCrSyeXijlzG/bLTtmIsjyU5fdenzVWBPHpN424/Z+b0BUI4dSRJbhWp5EwYR7x5ksr4jMR0ZnLpap5gsh4Mv8WkVAxuboI9140Dif2UtiNMIZSAOGP25aEwYog9jV18u1yNDlz2jDrhkiItbc/aybmvq6T++AVuOz40/cmqSptiZ6Rn8Vh1rwscBIJYqBDR2mWYbVacA05JmmHFUDUebtxyNsFwJgzt6+pk4u1qJw5TZh19e5GAKlvFqyl0uNGdUkO9jfJ7+PeOcdjaHHP578SCqqcuSxIkRALIKifJUFkPhRmJYgkYPNZN+1vgSTJw9hL4oS+mTN3oLmL95lz2/XDrJ/sOorP9zbDYbPgzDGDeuMtqJg+sgwAMHNcBb5zYnb3dctExNmm2lYlmYgqzErOHEFkPHSUEkQSsJy5o5FB5CMH5ccd9VRdIhdAHPJ2oa07ACC6/QMLs3b4gnjgP18DkFuDDCmK7l2XahacPwbjh3pwyeQhfT6DdSCQbdWsYl4fjfIiiMwn4525YDCIu+++GyNGjEBOTg5GjhyJ3/72twiHw3wbSZJw3333oaqqCjk5OZgxYwa+/PJL1ev4fD7cfPPNKCsrQ15eHubMmYMDBw709dsh+glMzDHihVgBYFC+C26HFZIE7K6Xc+y0febYBIhVu45iy0Ev8l123Hz2qBSuOs76ClyYd+owunD3EuJ+zY4+c9klPglioJPxYu6Pf/wj/vrXv2LRokXYtm0bHnzwQfzP//wPHnvsMb7Ngw8+iIceegiLFi3CunXrUFlZifPOOw9tbW18m9tuuw3Lli3Dyy+/jFWrVqG9vR2zZ89GKBSK9WcJIi5RYi5O8QMgj6ZieXM76+XvpXaiA7todkZy6q7/1kiU9mIVK9F35GVZaxJRwJHAJ4jMJ+PF3Jo1a3DxxRfjwgsvxPDhw/Hd734XM2fOxPr16wHIrtwjjzyCX/3qV5g7dy7Gjx+P5557Dp2dnXjppZcAAF6vF8888wz+93//F+eeey4mT56MF154AVu2bMG7776bzrdHZClFeerB8yNizGTVUh0pKtgZ6R2nF2YFgLJ8F+afToUu/YVsrmal1iQEkflkvJg7/fTT8d5772HHjh0AgM2bN2PVqlX4f//v/wEAamtrUVdXh5kzZ/LnuFwunHnmmVi9ejUA4PPPP0cgEFBtU1VVhfHjx/NttPh8PrS2tqr+EQSjwGWHXWjdkciZA5QiiC42zsseu5oVkBv2UmPo/kPWhVkd1JqEILKJjD9Kf/GLX8Dr9WLs2LGw2WwIhUL4wx/+gO9///sAgLq6OgBARUWF6nkVFRXYu3cv38bpdKK4uDhqG/Z8LQsXLsRvfvObVL8dop9gsVhQlOvgBRDDS42LOYY2Z254aR4cNgtGlOXh8pOrU7dYIu1kWwGEGAqmmwqCyHwy/ih95ZVX8MILL+Cll17C8ccfj02bNuG2225DVVUVrrrqKr6dtgJPkqSEVXnxtrnrrrtw++23859bW1tRXU0XWEKhKNeJo+1+DCnKMeS21GjEnPY5lR433r9jBjy5DjhsGW+aEyYQe7VpHdlMJEfVmiTz10sQA52MF3N33nknfvnLX+Lyyy8HAEyYMAF79+7FwoULcdVVV6GyshKA7L4NHqwMBa+vr+duXWVlJfx+P5qbm1XuXH19PaZPnx7z77pcLrhclHxO6FMSKYIYkaCSlaEVc7Eu6lr3jugfsDBrjsOWFZM1qACCILKLjL/97+zshNWqXqbNZuOtSUaMGIHKykqsWLGC/97v9+PDDz/kQm3KlClwOByqbQ4fPoytW7fqijmCSERRZKSXkXw5ABharO4X5ybHY8AwqEC+MSzNz46ZymLOXDZMrCCIgU7GH6UXXXQR/vCHP6CmpgbHH388Nm7ciIceegjXXnstADm8etttt+H+++/H6NGjMXr0aNx///3Izc3FFVdcAQDweDyYP38+7rjjDpSWlqKkpAQLFizAhAkTcO6556bz7RFZzJjKArzz1RHD47byXHaU5Tt5np3bnvH3UkSKGFaah//93iQMK80O51U9m5VuOggi08l4MffYY4/hnnvuwY033oj6+npUVVXhhhtuwK9//Wu+zc9//nN0dXXhxhtvRHNzM6ZOnYp33nkHBQUFfJuHH34Ydrsdl156Kbq6unDOOedgyZIlsNnoREUkx63njMZFk6owujxxWxLG0OJcRcxlwVgnInV8Z8rQdC/BMOrZrBl/mSCIAY9FkiQp3YvIBlpbW+HxeOD1elFYWJju5RBZyi3/2IjXNx8CALz642k4eXhJmldEENGEwhKO+f/eBgC8cdPpmDDUk+YVEcTAw4zuoDgPQfQhYhFENlQ1EgMTm9UCVyQNgJoGE0TmQ/45QfQh1SVKEYS2zxxBZBJXnzYc3zR0YISBHooEQaQXEnME0YeIrUcoZ47IZO664Lh0L4EgCIOQNUAQfQibzwqQmCMIgiBSAzlzBNGHDPa4UVHogj8YhifHke7lEARBEP0AEnME0YfYbVa8fcsZCIUlOKnPHEEQBJECSMwRRB9Tmk9j4giCIIjUQdYAQRAEQRBEFkNijiAIgiAIIoshMUcQBEEQBJHFkJgjCIIgCILIYkjMEQRBEARBZDEk5giCIAiCILIYEnMEQRAEQRBZDPWZM4gkSQCA1tbWNK+EIAiCIIj+DtMbTH/Eg8ScQdra2gAA1dXVaV4JQRAEQRADhba2Nng8nrjbWCQjko9AOBzGoUOHUFBQAIvFku7l6NLa2orq6mrs378fhYWF6V5O2qH9oYb2hxraH2pof6ih/aGG9oea3t4fkiShra0NVVVVsFrjZ8WRM2cQq9WKoUOHpnsZhiksLKSDTYD2hxraH2pof6ih/aGG9oca2h9qenN/JHLkGFQAQRAEQRAEkcWQmCMIgiAIgshiSMz1M1wuF+699164XK50LyUjoP2hhvaHGtofamh/qKH9oYb2h5pM2h9UAEEQBEEQBJHFkDNHEARBEASRxZCYIwiCIAiCyGJIzBEEQRAEQWQxJOYIgiAIgiCyGBJzGcZHH32Eiy66CFVVVbBYLHjttddUvz9y5AiuvvpqVFVVITc3F7NmzcLOnTujXmfNmjU4++yzkZeXh6KiIsyYMQNdXV38983NzZg3bx48Hg88Hg/mzZuHlpaWXn535knF/qirq8O8efNQWVmJvLw8nHjiifjXv/6l2iZb9sfChQtx8skno6CgAOXl5fj2t7+N7du3q7aRJAn33XcfqqqqkJOTgxkzZuDLL79UbePz+XDzzTejrKwMeXl5mDNnDg4cOKDaJhv2SSr2R1NTE26++WaMGTMGubm5qKmpwS233AKv16t6nYGyP7TbXnDBBTGPvYG2P/rDOTVV+6O/nFON7I+lS5fi/PPPR1lZGSwWCzZt2hT1OhlxPpWIjOLtt9+WfvWrX0n//ve/JQDSsmXL+O/C4bB06qmnSmeccYb02WefSV9//bV0/fXXSzU1NVJ7ezvfbvXq1VJhYaG0cOFCaevWrdKOHTukV199Veru7ubbzJo1Sxo/fry0evVqafXq1dL48eOl2bNn9+VbNUQq9se5554rnXzyydLatWul3bt3S7/73e8kq9UqbdiwgW+TLfvj/PPPlxYvXixt3bpV2rRpk3ThhRdGvd8HHnhAKigokP79739LW7ZskS677DJp8ODBUmtrK9/mxz/+sTRkyBBpxYoV0oYNG6SzzjpLmjRpkhQMBvk22bBPUrE/tmzZIs2dO1d6/fXXpV27dknvvfeeNHr0aOk73/mO6m8NlP0h8tBDD0kXXHBB1LEnSQNrf/SXc2qq9kd/Oaca2R9///vfpd/85jfSU089JQGQNm7cGPU6mXA+JTGXwWhPoNu3b5cASFu3buWPBYNBqaSkRHrqqaf4Y1OnTpXuvvtu3df96quvJADSp59+yh9bs2aNBED6+uuvU/smUkiy+yMvL0/6+9//rnqtkpIS6emnn5YkKXv3hyRJUn19vQRA+vDDDyVJkgVuZWWl9MADD/Bturu7JY/HI/31r3+VJEmSWlpaJIfDIb388st8m4MHD0pWq1Vavny5JEnZu0+S2R+x+Oc//yk5nU4pEAhIkjQw98emTZukoUOHSocPH4469gba/uiv59Rk90d/Padq94dIbW1tTDGXKedTCrNmET6fDwDgdrv5YzabDU6nE6tWrQIA1NfXY+3atSgvL8f06dNRUVGBM888k/8ekMMFHo8HU6dO5Y+deuqp8Hg8WL16dR+9m55jZH8AwOmnn45XXnkFTU1NCIfDePnll+Hz+TBjxgwA2b0/WCiwpKQEAFBbW4u6ujrMnDmTb+NyuXDmmWfy9/L5558jEAiotqmqqsL48eP5Ntm6T5LZH3qvU1hYCLtdHl890PZHZ2cnvv/972PRokWorKyMet2BtD/68zk12e9Hfz2naveHETLlfEpiLosYO3Yshg0bhrvuugvNzc3w+/144IEHUFdXh8OHDwMAvvnmGwDAfffdhx/96EdYvnw5TjzxRJxzzjk8l6yurg7l5eVRr19eXo66urq+e0M9xMj+AIBXXnkFwWAQpaWlcLlcuOGGG7Bs2TIcc8wxALJ3f0iShNtvvx2nn346xo8fDwB8vRUVFaptKyoq+O/q6urgdDpRXFwcd5ts2yfJ7g8tjY2N+N3vfocbbriBPzbQ9sfPfvYzTJ8+HRdffHHM1x5I+6O/nlN78v3oj+fUWPvDCJlyPrWn5FWIPsHhcODf//435s+fj5KSEthsNpx77rm44IIL+DbhcBgAcMMNN+Caa64BAEyePBnvvfcenn32WSxcuBAAYLFYol5fkqSYj2cqRvYHANx9991obm7Gu+++i7KyMrz22mv43ve+h48//hgTJkwAkJ3746abbsIXX3yhcggY2nUbeS/abbJtn6Rif7S2tuLCCy/EuHHjcO+998Z9jXivkwkkuz9ef/11vP/++9i4cWPc1x8o+6O/nlN7crz0x3NqvP2RDH19PiVnLsuYMmUKNm3ahJaWFhw+fBjLly9HY2MjRowYAQAYPHgwAGDcuHGq5x133HHYt28fAKCyshJHjhyJeu2GhoaoO7JMJ9H+2L17NxYtWoRnn30W55xzDiZNmoR7770XJ510Ev7yl78AyM79cfPNN+P111/HypUrMXToUP44C4lp7/bq6+v5e6msrITf70dzc3PcbbJpn/RkfzDa2towa9Ys5OfnY9myZXA4HKrXGSj74/3338fu3btRVFQEu93OQ83f+c53eBhtIO2P/nhO7cn+6I/nVL39YYRMOZ+SmMtSPB4PBg0ahJ07d2L9+vU8HDJ8+HBUVVVFlVfv2LEDw4YNAwBMmzYNXq8Xn332Gf/92rVr4fV6MX369L57EylEb390dnYCAKxW9VfdZrPxO+5s2h+SJOGmm27C0qVL8f7773PRyhgxYgQqKyuxYsUK/pjf78eHH37I38uUKVPgcDhU2xw+fBhbt27l22TLPknF/gBkR27mzJlwOp14/fXXVXmYwMDaH7/85S/xxRdfYNOmTfwfADz88MNYvHgxgIG1P/rTOTUV+6M/nVMT7Q8jZMz5NCVlFETKaGtrkzZu3Cht3LhRAiA99NBD0saNG6W9e/dKkiRX2a1cuVLavXu39Nprr0nDhg2T5s6dq3qNhx9+WCosLJReffVVaefOndLdd98tud1uadeuXXybWbNmSRMnTpTWrFkjrVmzRpowYULGlY1LUs/3h9/vl0aNGiWdccYZ0tq1a6Vdu3ZJf/rTnySLxSK99dZbfLts2R8/+clPJI/HI33wwQfS4cOH+b/Ozk6+zQMPPCB5PB5p6dKl0pYtW6Tvf//7MVuTDB06VHr33XelDRs2SGeffXbMUvpM3yep2B+tra3S1KlTpQkTJki7du1Svc5A3B+xgE5rkoGyP/rLOTUV+6M/nVON7I/GxkZp48aN0ltvvSUBkF5++WVp48aN0uHDh/k2mXA+JTGXYaxcuVICEPXvqquukiRJkh599FFp6NChksPhkGpqaqS7775b8vl8Ua+zcOFCaejQoVJubq40bdo06eOPP1b9vrGxUfrBD34gFRQUSAUFBdIPfvADqbm5uQ/eoTlSsT927NghzZ07VyovL5dyc3OliRMnRpXVZ8v+iLUvAEiLFy/m24TDYenee++VKisrJZfLJX3rW9+StmzZonqdrq4u6aabbpJKSkqknJwcafbs2dK+fftU22TDPknF/tD7jgGQamtr+XYDZX/ova5WzA20/dEfzqmp2h/95ZxqZH8sXrw45jb33nsv3yYTzqeWyBsiCIIgCIIgshDKmSMIgiAIgshiSMwRBEEQBEFkMSTmCIIgCIIgshgScwRBEARBEFkMiTmCIAiCIIgshsQcQRAEQRBEFkNijiAIgiAIIoshMUcQBEEQBJHFkJgjCIJIkg8++AAWiwUtLS3pXgpBEAMYmgBBEARhkBkzZuCEE07AI488AkAeQt7U1ISKigpYLJb0Lo4giAGLPd0LIAiCyFacTicqKyvTvQyCIAY4FGYlCIIwwNVXX40PP/wQjz76KCwWCywWC5YsWaIKsy5ZsgRFRUV48803MWbMGOTm5uK73/0uOjo68Nxzz2H48OEoLi7GzTffjFAoxF/b7/fj5z//OYYMGYK8vDxMnToVH3zwQXreKEEQWQc5cwRBEAZ49NFHsWPHDowfPx6//e1vAQBffvll1HadnZ3485//jJdffhltbW2YO3cu5s6di6KiIrz99tv45ptv8J3vfAenn346LrvsMgDANddcgz179uDll19GVVUVli1bhlmzZmHLli0YPXp0n75PgiCyDxJzBEEQBvB4PHA6ncjNzeWh1a+//jpqu0AggCeeeALHHHMMAOC73/0unn/+eRw5cgT5+fkYN24czjrrLKxcuRKXXXYZdu/ejX/84x84cOAAqqqqAAALFizA8uXLsXjxYtx///199yYJgshKSMwRBEGkkNzcXC7kAKCiogLDhw9Hfn6+6rH6+noAwIYNGyBJEo499ljV6/h8PpSWlvbNogmCyGpIzBEEQaQQh8Oh+tliscR8LBwOAwDC4TBsNhs+//xz2Gw21XaiACQIgtCDxBxBEIRBnE6nqnAhFUyePBmhUAj19fU444wzUvraBEEMDKialSAIwiDDhw/H2rVrsWfPHhw9epS7az3h2GOPxQ9+8ANceeWVWLp0KWpra7Fu3Tr88Y9/xNtvv52CVRME0d8hMUcQBGGQBQsWwGazYdy4cRg0aBD27duXktddvHgxrrzyStxxxx0YM2YM5syZg7Vr16K6ujolr08QRP+GJkAQBEEQBEFkMeTMEQRBEARBZDEk5giCIAiCILIYEnMEQRAEQRBZDIk5giAIgiCILIbEHEEQBEEQRBZDYo4gCIIgCCKLITFHEARBEASRxZCYIwiCIAiCyGJIzBEEQRAEQWQxJOYIgiAIgiCyGBJzBEEQBEEQWcz/D5lTYWnU84FwAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=[7, 3])\n", "ds.isel(station_num=0).total_precip.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
WARNING\n", " \n", "Currently, there is no way to provide `Extremes.jl` with a predefined set of parameters to directly calculate return levels. Until this functionality is implemented in either `xHydro` or `Extremes.jl`, the `.fit()` and `.return_level()` functions should be considered independent. Specifically, the `.return_level()` function will first estimate the distribution parameters before calculating the return levels.\n", " \n", "
\n", "\n", "## Parameter estimation\n", "\n", "The `xhydro.extreme_value_analysis.fit` function serves as the interface between `xHydro` and the `Extremes.jl` package. Most of the arguments mirror those used in the `xhydro.frequency_analysis.local.fit` function. The statistical distribution names have been made to align with those in `SciPy`. Below are a few key differences:\n", "\n", "- Bayesian Method (`BAYES`): When using the `BAYES` method, you can specify two additional parameters:\n", " - `niter`: Number of iterations for the Bayesian inference algorithm.\n", " - `warmup`: Number of warmup iterations for the Bayesian inference.\n", "- Confidence Intervals: A significant addition to this function is the `confidence_level` parameter, which simplifies the process of obtaining confidence interval compared to the other options available in `xHydro`, as detailed in the other frequency analysis notebooks.\n", "\n", "In this example, we will estimate a Generalized Extreme Value (GEV) distribution (`genextreme`) using the Probability Weighted Moments (`PWM`) method. Additionally, we will calculate and return the 95% confidence intervals for the estimated parameters.\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function fit in module xhydro.extreme_value_analysis.parameterestimation:\n", "\n", "fit(ds: 'xr.Dataset', locationcov: 'list[str] | None' = None, scalecov: 'list[str] | None' = None, shapecov: 'list[str] | None' = None, variables: 'list[str] | None' = None, dist: 'str | scipy.stats.rv_continuous' = 'genextreme', method: 'str' = 'ML', dim: 'str' = 'time', confidence_level: 'float' = 0.95, niter: 'int' = 5000, warmup: 'int' = 2000) -> 'xr.Dataset'\n", " Fit an array to a univariate distribution along a given dimension.\n", "\n", " Parameters\n", " ----------\n", " ds : xr.DataSet\n", " Xarray Dataset containing the data to be fitted.\n", " locationcov : list[str]\n", " List of names of the covariates for the location parameter.\n", " scalecov : list[str]\n", " List of names of the covariates for the scale parameter.\n", " shapecov : list[str]\n", " List of names of the covariates for the shape parameter.\n", " variables : list[str]\n", " List of variables to be fitted.\n", " dist : str or rv_continuous distribution object\n", " Name of the univariate distribution or the distribution object itself.\n", " Supported distributions are genextreme, gumbel_r, genpareto.\n", " method : {\"ML\", \"PWM\", \"BAYES}\n", " Fitting method, either maximum likelihood (ML), probability weighted moments (PWM) or bayesian (BAYES).\n", " dim : str\n", " Specifies the dimension along which the fit will be performed (default: \"time\").\n", " confidence_level : float\n", " The confidence level for the confidence interval of each parameter.\n", " niter : int\n", " The number of iterations of the bayesian inference algorithm for parameter estimation (default: 5000).\n", " warmup : int\n", " The number of warmup iterations of the bayesian inference algorithm for parameter estimation (default: 2000).\n", "\n", " Returns\n", " -------\n", " xr.Dataset\n", " Dataset of fitted distribution parameters and confidence interval values.\n", "\n", " Notes\n", " -----\n", " Coordinates for which all values are NaNs will be dropped before fitting the distribution. If the array still\n", " contains NaNs or has less valid values than the number of parameters for that distribution,\n", " the distribution parameters will be returned as NaNs.\n", "\n" ] } ], "source": [ "help(xhe.fit)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:35.121255Z", "start_time": "2024-12-11T20:54:28.257188Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 460B\n",
       "Dimensions:             (station_num: 5, dparams: 3)\n",
       "Coordinates:\n",
       "  * station_num         (station_num) int64 40B 1001 1004 1008 1009 1012\n",
       "  * dparams             (dparams) <U5 60B 'shape' 'loc' 'scale'\n",
       "Data variables:\n",
       "    total_precip        (station_num, dparams) float64 120B 0.2012 ... 155.8\n",
       "    total_precip_lower  (station_num, dparams) float64 120B 0.08579 ... 135.3\n",
       "    total_precip_upper  (station_num, dparams) float64 120B 0.3243 ... 174.1
" ], "text/plain": [ " Size: 460B\n", "Dimensions: (station_num: 5, dparams: 3)\n", "Coordinates:\n", " * station_num (station_num) int64 40B 1001 1004 1008 1009 1012\n", " * dparams (dparams) 'xr.Dataset'\n", " Compute the return level associated with a return period based on a given distribution.\n", "\n", " Parameters\n", " ----------\n", " ds : xr.DataSet\n", " Xarray Dataset containing the data for return level calculations.\n", " locationcov : list[str]\n", " List of names of the covariates for the location parameter.\n", " scalecov : list[str]\n", " List of names of the covariates for the scale parameter.\n", " shapecov : list[str]\n", " List of names of the covariates for the shape parameter.\n", " variables : list[str]\n", " List of variables to be fitted.\n", " dist : str or rv_continuous distribution object\n", " Name of the univariate distribution or the distribution object itself.\n", " Supported distributions are genextreme, gumbel_r, genpareto.\n", " method : {\"ML\", \"PWM\", \"BAYES}\n", " Fitting method, either maximum likelihood (ML), probability weighted moments (PWM) or bayesian (BAYES).\n", " dim : str\n", " Specifies the dimension along which the fit will be performed (default: \"time\").\n", " confidence_level : float\n", " The confidence level for the confidence interval of each parameter.\n", " return_period : float\n", " Return period used to compute the return level.\n", " niter : int\n", " The number of iterations of the bayesian inference algorithm for parameter estimation (default: 5000).\n", " warmup : int\n", " The number of warmup iterations of the bayesian inference algorithm for parameter estimation (default: 2000).\n", " threshold_pareto : float\n", " The value above which the Pareto distribution is applied.\n", " nobs_pareto : int\n", " The total number of observations used when applying the Pareto distribution.\n", " nobsperblock_pareto : int\n", " The number of observations per block when applying the Pareto distribution.\n", "\n", " Returns\n", " -------\n", " xr.Dataset\n", " Dataset of with the return level and the confidence interval values.\n", "\n", " Notes\n", " -----\n", " Coordinates for which all values are NaNs will be dropped before fitting the distribution. If the array still\n", " contains NaNs or has less valid values than the number of parameters for that distribution,\n", " the distribution parameters will be returned as NaNs.\n", "\n" ] } ], "source": [ "help(xhe.return_level)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:37.882475Z", "start_time": "2024-12-11T20:54:35.133335Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 168B\n",
       "Dimensions:             (station_num: 5, return_period: 1)\n",
       "Coordinates:\n",
       "  * station_num         (station_num) int64 40B 1001 1004 1008 1009 1012\n",
       "  * return_period       (return_period) int64 8B 100\n",
       "Data variables:\n",
       "    total_precip        (station_num, return_period) float64 40B 1.704e+03 .....\n",
       "    total_precip_lower  (station_num, return_period) float64 40B 1.609e+03 .....\n",
       "    total_precip_upper  (station_num, return_period) float64 40B 1.8e+03 ... ...
" ], "text/plain": [ " Size: 168B\n", "Dimensions: (station_num: 5, return_period: 1)\n", "Coordinates:\n", " * station_num (station_num) int64 40B 1001 1004 1008 1009 1012\n", " * return_period (return_period) int64 8B 100\n", "Data variables:\n", " total_precip (station_num, return_period) float64 40B 1.704e+03 .....\n", " total_precip_lower (station_num, return_period) float64 40B 1.609e+03 .....\n", " total_precip_upper (station_num, return_period) float64 40B 1.8e+03 ... ..." ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rtnlv_stationary = xhe.return_level(\n", " ds,\n", " dist=\"gumbel_r\",\n", " method=\"ML\",\n", " variables=[\"total_precip\"],\n", " confidence_level=0.95,\n", " return_period=100,\n", ")\n", "rtnlv_stationary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Non-stationary model\n", "\n", "So far, we've skipped three additional arguments—`locationcov`, `scalecov`, and `shapecov`—that accept variable names. These arguments allow you to introduce a non-linear aspect to the statistical model. In non-stationary models, explanatory variables (covariates) can be used to capture changes in model parameters over time or across different conditions. These covariates can represent factors such as time, geographic location, global temperature increases or CO2 concentrations, or any other variable that may influence the distribution parameters.\n", "\n", "Also, note that the `PWM` method cannot be used with non-stationary models.\n", "\n", "For this example, we'll keep it simple and assume that the location parameter varies as a linear function of the year. To do this, we'll need to add a new variable containing the year to our dataset and then provide this variable to the `locationcov` argument.\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:37.928189Z", "start_time": "2024-12-11T20:54:37.909414Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 19kB\n",
       "Dimensions:       (station_num: 5, time: 146)\n",
       "Coordinates:\n",
       "  * station_num   (station_num) int64 40B 1001 1004 1008 1009 1012\n",
       "  * time          (time) datetime64[ns] 1kB 1955-01-01 1956-01-01 ... 2100-01-01\n",
       "    station_name  (station_num, time) object 6kB 'Dozois' ... 'Lac St-Francois'\n",
       "Data variables:\n",
       "    total_precip  (station_num, time) float64 6kB 983.6 1.204e+03 ... 1.065e+03\n",
       "    year          (station_num, time) int64 6kB 1955 1956 1957 ... 2099 2100
" ], "text/plain": [ " Size: 19kB\n", "Dimensions: (station_num: 5, time: 146)\n", "Coordinates:\n", " * station_num (station_num) int64 40B 1001 1004 1008 1009 1012\n", " * time (time) datetime64[ns] 1kB 1955-01-01 1956-01-01 ... 2100-01-01\n", " station_name (station_num, time) object 6kB 'Dozois' ... 'Lac St-Francois'\n", "Data variables:\n", " total_precip (station_num, time) float64 6kB 983.6 1.204e+03 ... 1.065e+03\n", " year (station_num, time) int64 6kB 1955 1956 1957 ... 2099 2100" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds[\"year\"] = ds.time.dt.year.broadcast_like(ds[\"total_precip\"])\n", "ds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the case of the `.fit()` function, adding a covariate will introduce a new entry under the `dparams` dimension. For this example, it created a new entry called `loc_year_covariate` under the `dparams` dimension.\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:40.472181Z", "start_time": "2024-12-11T20:54:37.945894Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 808B\n",
       "Dimensions:             (station_num: 5, dparams: 4)\n",
       "Coordinates:\n",
       "  * station_num         (station_num) int64 40B 1001 1004 1008 1009 1012\n",
       "  * dparams             (dparams) <U18 288B 'shape' 'loc' ... 'scale'\n",
       "Data variables:\n",
       "    total_precip        (station_num, dparams) float64 160B 0.2051 ... 144.3\n",
       "    total_precip_lower  (station_num, dparams) float64 160B 0.1022 ... 127.0\n",
       "    total_precip_upper  (station_num, dparams) float64 160B 0.308 ... 164.0
" ], "text/plain": [ " Size: 808B\n", "Dimensions: (station_num: 5, dparams: 4)\n", "Coordinates:\n", " * station_num (station_num) int64 40B 1001 1004 1008 1009 1012\n", " * dparams (dparams) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 25kB\n",
       "Dimensions:             (station_num: 5, time: 146, return_period: 1)\n",
       "Coordinates:\n",
       "  * station_num         (station_num) int64 40B 1001 1004 1008 1009 1012\n",
       "  * time                (time) datetime64[ns] 1kB 1955-01-01 ... 2100-01-01\n",
       "    station_name        (station_num, time) object 6kB 'Dozois' ... 'Lac St-F...\n",
       "  * return_period       (return_period) int64 8B 100\n",
       "Data variables:\n",
       "    total_precip        (station_num, time) float64 6kB 1.557e+03 ... 1.771e+03\n",
       "    total_precip_lower  (station_num, time) float64 6kB 1.457e+03 ... 1.672e+03\n",
       "    total_precip_upper  (station_num, time) float64 6kB 1.656e+03 ... 1.87e+03
" ], "text/plain": [ " Size: 25kB\n", "Dimensions: (station_num: 5, time: 146, return_period: 1)\n", "Coordinates:\n", " * station_num (station_num) int64 40B 1001 1004 1008 1009 1012\n", " * time (time) datetime64[ns] 1kB 1955-01-01 ... 2100-01-01\n", " station_name (station_num, time) object 6kB 'Dozois' ... 'Lac St-F...\n", " * return_period (return_period) int64 8B 100\n", "Data variables:\n", " total_precip (station_num, time) float64 6kB 1.557e+03 ... 1.771e+03\n", " total_precip_lower (station_num, time) float64 6kB 1.457e+03 ... 1.672e+03\n", " total_precip_upper (station_num, time) float64 6kB 1.656e+03 ... 1.87e+03" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rtnlv_non_stationary = xhe.return_level(\n", " ds,\n", " dist=\"gumbel_r\",\n", " method=\"ML\",\n", " dim=\"time\",\n", " variables=[\"total_precip\"],\n", " locationcov=[\"year\"],\n", " confidence_level=0.95,\n", " return_period=100,\n", ")\n", "\n", "rtnlv_non_stationary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparison of the return level using the stationary and non-stationary model" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAFICAYAAACP747yAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAARK9JREFUeJzt3XlcVPX+P/DXiOwDA4hsgQxiKKigmYq454IaoKKJVyPU3G+uqF26LVo/MzMUl0wrg1y+mamYZVJKrmluiYaSpRd3iVIcHEBA+Pz+4HKu4ww4o4MzOq/n4zGPx8znfM457zMyzns+57PIhBACRERERCZWz9QBEBEREQFMSoiIiMhMMCkhIiIis8CkhIiIiMwCkxIiIiIyC0xKiIiIyCwwKSEiIiKzwKSEiIiIzEJ9UwfwuKisrMTVq1fh5OQEmUxm6nCIiIgeG0II3Lp1Cz4+PqhXr+b2ECYlerp69Sr8/PxMHQYREdFj69KlS/D19a1xO5MSPTk5OQGoekOdnZ1NHA0REdHjo7CwEH5+ftJ3aU2YlOip+paNs7MzkxIiIqIHcL/uD+zoSkRERGaBSQkRERGZBSYlREREZBaYlBAREZFZYFJCREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgTO6EhGRWRFCoLi8GADgYO3ARVAtCFtKiIjIrBSXF0M+Tw75PLmUnJBlYFJCREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGaBSQkRERGZBSYlREREZBaYlBAREZFZYFJCREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGaBSQkRERGZBSYlREREZBaYlBAREZFZYFJCREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGaBSQkRERGZBSYlREREZBaYlBAREZFZYFJCREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWTJqU7N27F9HR0fDx8YFMJsOWLVs0tqvVarzyyivw9fWFvb09goOD8dFHH2nUKS0txaRJk+Du7g5HR0fExMTg8uXLGnUKCgoQHx8PhUIBhUKB+Ph43Lx5s46vjoiIiAxh0qSkqKgIYWFhWLZsmc7t06ZNQ0ZGBtauXYucnBxMmzYNkyZNwtdffy3VmTp1KtLT07F+/Xrs378farUaUVFRqKiokOoMGzYMWVlZyMjIQEZGBrKyshAfH1/n10dERET6q2/Kk/ft2xd9+/atcfvBgweRkJCAbt26AQDGjh2LlStX4ujRo+jfvz9UKhVWrVqFNWvWoGfPngCAtWvXws/PDzt37kRkZCRycnKQkZGBn3/+Ge3btwcAfPLJJ+jQoQPOnDmDpk2b1vl1EhER0f2ZdZ+STp06YevWrbhy5QqEENi1axd+//13REZGAgCOHTuG8vJy9O7dW9rHx8cHLVq0wIEDBwBUJTYKhUJKSAAgPDwcCoVCqqNLaWkpCgsLNR70ZBBCoKioCEVFRRBCmDocIiL6L7NOSpYsWYKQkBD4+vrCxsYGffr0wfLly9GpUycAQF5eHmxsbODq6qqxn6enJ/Ly8qQ6Hh4eWsf28PCQ6ugyb948qQ+KQqGAn5+fEa+MTKm4uBhyuRxyuRzFxcWmDoeIiP7L7JOSn3/+GVu3bsWxY8eQnJyMiRMnYufOnbXuJ4SATCaTXt/9vKY690pKSoJKpZIely5devALISIiovsyaZ+S2pSUlOC1115Deno6nn/+eQBAaGgosrKy8MEHH6Bnz57w8vJCWVkZCgoKNFpL8vPzERERAQDw8vLCn3/+qXX8v/76C56enjWe39bWFra2tka+KiIiIqqJ2baUlJeXo7y8HPXqaYZoZWWFyspKAECbNm1gbW2NHTt2SNuvXbuG7OxsKSnp0KEDVCoVDh8+LNU5dOgQVCqVVIeIiIhMz6QtJWq1GmfPnpVe5+bmIisrC25ubmjUqBG6du2KmTNnwt7eHv7+/tizZw9Wr16NhQsXAgAUCgVefvllJCYmokGDBnBzc8OMGTPQsmVLaTROcHAw+vTpgzFjxmDlypUAqkbxREVFceQNERGRGTFpUnL06FF0795dej19+nQAQEJCAtLS0rB+/XokJSVh+PDhuHHjBvz9/TF37lyMHz9e2mfRokWoX78+hgwZgpKSEvTo0QNpaWmwsrKS6qxbtw6TJ0+WRunExMTUODcKERERmYZMcEykXgoLC6FQKKBSqeDs7GzqcOghFBUVQS6XA6hqrXN0dDRxRER0t6KyIsjn/fczmqSGow0/o487fb9DzbZPCREREVkWJiVERERkFpiUEBERkVlgUkJERERm4aGSktLSUmPFQURERBbOoKTk+++/x4gRIxAYGAhra2s4ODjAyckJXbt2xdy5c3H16tW6ipOIiIiecHolJVu2bEHTpk2RkJCAevXqYebMmdi8eTO+//57rFq1Cl27dsXOnTvRuHFjjB8/Hn/99Vddx01ERERPGL0mT3v33XfxwQcf4Pnnn9ea9h0AhgwZAgC4cuUKFi9ejNWrVyMxMdG4kRIREdETTa+k5O51Y2rz1FNP4f3333+ogIiIyLJVVFZIz/de2Ivegb1hVc+qlj3oScHRN0REZDY252xG8PJg6XW//+sH5WIlNudsNmFU9KgYvPaNEAIbN27Erl27kJ+fL63YW23zZv7hEBGR4TbnbMbgDYMhoLn6yZXCKxi8YTA2DtmI2OBYE0VHj4LBLSVTpkxBfHw8cnNzIZfLoVAoNB5ERESGqqiswJSMKVoJCQCpbGrGVI1bO/TkMbilZO3atdi8eTP69etXF/EQEZEF2ndxHy4XXq5xu4DApcJL2HdxH7opuz26wOiRMrilRKFQoHHjxnURCxERWahrt64ZtR49ngxOSmbPno05c+agpKSkLuIhIiIL5O3kbdR69Hgy+PbNCy+8gC+++AIeHh5QKpWwtrbW2P7LL78YLTgiIrIMnRt1hq+zL64UXtHZr0QGGXydfdG5UWcTREePisFJyYgRI3Ds2DG8+OKL8PT0hEwmq4u4iIjIgljVs8LiPosxeMNgrW0yVH3PpPRJ4XwlTziDk5Jt27bh+++/R6dOneoiHiIislCxwbHYOGQjJm2fhKu3/reWmq+zL1L6pHA4sAUwOCnx8/ODs7NzXcRCREQWLjY4Fj0DekIxv2qKie+GfccZXS2IwR1dk5OTMWvWLJw/f74OwiEiIkt3dwLSxb8LExILYnBLyYsvvoji4mIEBgbCwcFBq6PrjRs3jBYcERERWQ6Dk5KUlJQ6CIOIiIgsncFJSUJCQl3EQURERBbO4KSkWn5+vs4F+UJDQx86KCIiIrI8Biclx44dQ0JCAnJyciCE5gQ3MpkMFRVcLImIiIgMZ3BSMnLkSAQFBWHVqlWcPI2IiIiMxuCkJDc3F5s3b0aTJk3qIh4iIiKyUAbPU9KjRw+cOHGiLmIhIiIiC2ZwS8mnn36KhIQEZGdno0WLFlrzlMTExBgtOCIiIrIcBiclBw4cwP79+7F9+3atbezoSkRERA/K4Ns3kydPRnx8PK5du4bKykqNBxMSIiIielAGJyXXr1/HtGnT4OnpWRfxEBERkYUyOCmJjY3Frl276iIWIiIismAG9ykJCgpCUlIS9u/fj5YtW2p1dJ08ebLRgiMiIiLLIRP3Tst6HwEBATUfTCbDf/7zn4cOyhwVFhZCoVBApVLB2dnZ1OHQQygqKoJcLgcAqNVqODo6mjgiIrpbUVkR5PP++xlNUsPRhp/Rx52+36EPNHkaERERkbEZ3KeEiIiIqC7olZS89957KC4u1uuAhw4dwrZt2x4qKCIiIrI8eiUlp0+fRqNGjTBhwgRs374df/31l7Ttzp07OHnyJJYvX46IiAgMHTqUfS6IiIjIYHr1KVm9ejVOnjyJDz/8EMOHD4dKpYKVlRVsbW2lFpTWrVtj7NixSEhIgK2tbZ0GTURERE8evTu6hoaGYuXKlVixYgVOnjyJ8+fPo6SkBO7u7mjVqhXc3d3rMk4iIiJ6whk8+kYmkyEsLAxhYWF1EQ8RERFZKI6+ISIiIrPApISIiIjMApMSIiIiMgtMSoiIiMgsMCkhIiIis2Dw6JuioiK89957yMzMRH5+PiorKzW2P6kL8hEREVHdMjgpGT16NPbs2YP4+Hh4e3tDJpPVRVxERERkYQxOSrZv345t27ahY8eOdREPERERWSiD+5S4urrCzc2tLmIhIiIiC2ZwUvLOO+/gzTff1HvVYCIiIiJ9GHz7Jjk5GefOnYOnpyeUSiWsra01tv/yyy9GC46IiIgsh8FJyYABA+ogDKJHp6KiQnq+d+9e9O7dG1ZWViaMiIiIgAdISt566626iIPokdi8eTMmTZokve7Xrx98fX2xePFixMbGmjAyIiIyOCmpduzYMeTk5EAmkyEkJAStW7c2ZlxERrd582YMHjwYQgiN8itXrmDw4MHYuHEjExMiIhMyOCnJz8/H0KFDsXv3bri4uEAIAZVKhe7du2P9+vVo2LBhXcRJ9FAqKiowZcoUrYQEAIQQkMlkmDp1Kvr3789bOUREJmLw6JtJkyahsLAQp06dwo0bN1BQUIDs7GwUFhZi8uTJBh1r7969iI6Oho+PD2QyGbZs2aKxXSaT6XwsWLBAqlNaWopJkybB3d0djo6OiImJweXLlzWOU1BQgPj4eCgUCigUCsTHx+PmzZuGXjo9xvbt26f1d3E3IQQuXbqEffv2PcKoiIjobgYnJRkZGfjoo48QHBwslYWEhODDDz/E9u3bDTpWUVERwsLCsGzZMp3br127pvH47LPPIJPJMGjQIKnO1KlTkZ6ejvXr12P//v1Qq9WIiorS6Mw4bNgwZGVlISMjAxkZGcjKykJ8fLyBV06Ps2vXrhm1HhERGZ/Bt28qKyu1hgEDgLW1tdY6OPfTt29f9O3bt8btXl5eGq+//vprdO/eHY0bNwYAqFQqrFq1CmvWrEHPnj0BAGvXroWfnx927tyJyMhI5OTkICMjAz///DPat28PAPjkk0/QoUMHnDlzBk2bNjUoZno8eXt7G7UeEREZn8EtJc899xymTJmCq1evSmVXrlzBtGnT0KNHD6MGd7c///wT27Ztw8svvyyVHTt2DOXl5ejdu7dU5uPjgxYtWuDAgQMAgIMHD0KhUEgJCQCEh4dDoVBIdXQpLS1FYWGhxoMeX507d4avr2+NazXJZDL4+fmhc+fOjzgyIiKqZnBSsmzZMty6dQtKpRKBgYFo0qQJAgICcOvWLSxdurQuYgQAfP7553ByctIYHZGXlwcbGxu4urpq1PX09EReXp5Ux8PDQ+t4Hh4eUh1d5s2bJ/VBUSgU8PPzM9KVkClYWVlh8eLFOrdVJyopKSns5EpEZEIG377x8/PDL7/8gh07duC3336DEAIhISHS7ZO68tlnn2H48OGws7O7b93q0RTVdP06vrfOvZKSkjB9+nTpdWFhIROTx1xsbCw2btyISZMmabT0+fr6IiUlhcOBiYhM7IHnKenVqxd69eplzFhqtG/fPpw5cwZffvmlRrmXlxfKyspQUFCg0VqSn5+PiIgIqc6ff/6pdcy//voLnp6eNZ7T1tYWtra2RroCMhexsbHo2bMnFAoFAOC7777jjK5ERGZCr6RkyZIlGDt2LOzs7LBkyZJa6xo6LFgfq1atQps2bRAWFqZR3qZNG1hbW2PHjh0YMmQIgKrRE9nZ2Xj//fcBAB06dIBKpcLhw4fRrl07AMChQ4egUqmkxIUsy90JSJcuXZiQEBGZCb2SkkWLFkm3ThYtWlRjPZlMZlBSolarcfbsWel1bm4usrKy4ObmhkaNGgGoum3y1VdfITk5WWt/hUKBl19+GYmJiWjQoAHc3NwwY8YMtGzZUrqdFBwcjD59+mDMmDFYuXIlAGDs2LGIioriyBsiIiIzoldSkpubq/P5wzp69Ci6d+8uva7uw5GQkIC0tDQAwPr16yGEwD/+8Q+dx1i0aBHq16+PIUOGoKSkBD169EBaWprGr99169Zh8uTJ0iidmJiYGudGISIiItOQCV3zbtfi7bffxowZM+Dg4KBRXlJSggULFuDNN980aoDmorCwEAqFAiqVCs7OzqYOhx5CUVER5HI5gKrWOkdHRxNHRER3Kyorgnzefz+jSWo42vAz+rjT9zvU4CHBc+bMgVqt1iovLi7GnDlzDD0cEREREYAHSEpqGkp74sQJuLm5GSUoIiIisjx6Dwl2dXWVFsQLCgrSSEwqKiqgVqsxfvz4OgmSiIiInnx6JyUpKSkQQmDUqFGYM2eONM8DANjY2ECpVKJDhw51EiQRERE9+fROShISEgAAAQEBiIiI0LkoHxEREdGDMnhG165du0rPS0pKUF5errGdI1OIiIjoQRjc0bW4uBivvPIKPDw8IJfL4erqqvEgIiIiehAGJyUzZ87Ejz/+iOXLl8PW1haffvop5syZAx8fH6xevbouYiQiIiILYPDtm2+++QarV69Gt27dMGrUKHTu3BlNmjSBv78/1q1bh+HDh9dFnERERPSEM7il5MaNGwgICABQ1X/kxo0bAIBOnTph7969xo2OiIiILIbBSUnjxo1x/vx5AEBISAg2bNgAoKoFxcXFxZixERERkQUxOCkZOXIkTpw4AQBISkqS+pZMmzYNM2fONHqAREREZBkM7lMybdo06Xn37t3x22+/4ejRowgMDERYWJhRgyMiIiLLYXBLyerVq1FaWiq9btSoEWJjYxEcHMzRN0RE9NAcrB2gTlJDnaSGg7XD/XegJ8YD3b5RqVRa5bdu3cLIkSONEhQREVkumUwGRxtHONo46lwAlp5cRlsl+PLlyxrr4RAREREZQu8+Ja1bt5ZWCe7Rowfq1//frhUVFcjNzUWfPn3qJEgiIiJ68umdlAwYMAAAkJWVhcjISMjlcmlb9SrBgwYNMnqAREREZBn0TkreeustAIBSqURcXBzs7OzqLCgiIiKyPAYPCU5ISKiLOIiIiMjC6ZWUuLm54ffff4e7uztcXV1r7Q1dPe08ERERkSH0SkoWLVoEJycnAEBKSkpdxkNEZBghgOLiqucODgCHkBI9tmRCCGHqIB4HhYWFUCgUUKlUcHZ2NnU49BCKioqkjtpqtRqOjo4mjogeSlERUN3xXq0G+O9JZHb0/Q41uE8JUDUEOD09HTk5OZDJZAgODkb//v01hgkTERERGcLgLCI7Oxv9+/dHXl4emjZtCgD4/fff0bBhQ2zduhUtW7Y0epBERET05DN4RtfRo0ejefPmuHz5Mn755Rf88ssvuHTpEkJDQzF27Ni6iJGIiIgsgMEtJSdOnMDRo0fh6uoqlbm6umLu3Llo27atUYMjIiIiy2FwS0nTpk3x559/apXn5+ejSZMmRgmKiIiILI/BScm7776LyZMnY+PGjbh8+TIuX76MjRs3YurUqZg/fz4KCwulBxEREZG+DB4SXK/e//KY6knUqg9x92uZTIaKigpjxWlyHBL85OCQ4CcMhwQTmb06GxK8a9euhwqMiIiISBeDk5KuXbvWRRxERERk4fRKSk6ePIkWLVqgXr16OHnyZK11Q0NDjRIYERERWRa9kpJWrVohLy8PHh4eaNWqFWQyGXR1RXnS+pEQERHRo6NXUpKbm4uGDRtKz4mIiIiMTa+kxN/fX+dzIiIiImMxeJ6SefPm4bPPPtMq/+yzzzB//nyjBEVERESWx+CkZOXKlWjWrJlWefPmzbFixQqjBEVERESWx+CkJC8vD97e3lrlDRs2xLVr14wSFBEREVkeg5MSPz8//PTTT1rlP/30E3x8fIwSFBEREVkegydPGz16NKZOnYry8nI899xzAIDMzEzMmjULiYmJRg+QiIiILIPBScmsWbNw48YNTJw4EWVlZQAAOzs7vPrqq0hKSjJ6gERERGQZDF6Qr5parUZOTg7s7e3x9NNPw9bW1tixmRUuyPfk4IJ8TxguyEdk9vT9DjW4T0m1vLw83LhxA4GBgbC1tdU5wysRERGRvgxOSq5fv44ePXogKCgI/fr1k0bcjB49mn1KiIiI6IEZnJRMmzYN1tbWuHjxIhwcHKTyuLg4ZGRkGDU4IiIishwGd3T94Ycf8P3338PX11ej/Omnn8aFCxeMFhgRERFZFoNbSoqKijRaSKr9/fffT3xnVyIiIqo7BiclXbp0werVq6XXMpkMlZWVWLBgAbp3727U4IiIiMhyGHz7ZsGCBejWrRuOHj2KsrIyzJo1C6dOncKNGzd0zvRKREREpA+DW0pCQkJw8uRJtGvXDr169UJRURFiY2Nx/PhxBAYG1kWMREREZAEMaikpLy9H7969sXLlSsyZM6euYiIiIiILZFBLibW1NbKzsyGTyeoqHiIiIrJQBt++eemll7Bq1aq6iIWIiIgsmMEdXcvKyvDpp59ix44dePbZZ7XWDVm4cKHRgiMiIiLLYXBSkp2djWeeeQYA8Pvvv2ts420dIiIielAGJyW7du2qiziIHhkHBweo1WrpORERmYcHXiUYAC5duoTLly8/8P579+5FdHQ0fHx8IJPJsGXLFq06OTk5iImJgUKhgJOTE8LDw3Hx4kVpe2lpKSZNmgR3d3c4OjoiJiZGK6aCggLEx8dDoVBAoVAgPj4eN2/efOC46fEmk8ng6OgIR0dHtu4REZkRg5OSO3fu4I033oBCoYBSqYS/vz8UCgVef/11lJeXG3SsoqIihIWFYdmyZTq3nzt3Dp06dUKzZs2we/dunDhxAm+88Qbs7OykOlOnTkV6ejrWr1+P/fv3Q61WIyoqChUVFVKdYcOGISsrCxkZGcjIyEBWVhbi4+MNvXQiIiKqS8JA48aNEx4eHmLFihXixIkT4sSJE2LFihXCy8tLjBs3ztDDSQCI9PR0jbK4uDjx4osv1rjPzZs3hbW1tVi/fr1UduXKFVGvXj2RkZEhhBDi9OnTAoD4+eefpToHDx4UAMRvv/2md3wqlUoAECqVSu99iOgRUKuFAKoearWpoyEiHfT9DjW4peSLL75AWloaxo0bh9DQUISGhmLcuHH47LPP8MUXXxgtWaqsrMS2bdsQFBSEyMhIeHh4oH379hq3eI4dOyZN6FbNx8cHLVq0wIEDBwAABw8ehEKhQPv27aU64eHhUCgUUh1dSktLUVhYqPEgIiKiumNwUmJnZwelUqlVrlQqYWNjY4yYAAD5+flQq9V477330KdPH/zwww8YOHAgYmNjsWfPHgBAXl4ebGxs4OrqqrGvp6cn8vLypDoeHh5ax/fw8JDq6DJv3jypD4pCoYCfn5/Rro2IiIi0GZyU/POf/8Q777yD0tJSqay0tBRz587FK6+8YrTAKisrAQD9+/fHtGnT0KpVK/zrX/9CVFQUVqxYUeu+QgiNDoy6OjPeW+deSUlJUKlU0uPSpUsPeCVERESkD4OHBB8/fhyZmZnw9fVFWFgYAODEiRMoKytDjx49EBsbK9XdvHnzAwfm7u6O+vXrIyQkRKM8ODgY+/fvBwB4eXmhrKwMBQUFGq0l+fn5iIiIkOr8+eefWsf/66+/4OnpWeP5bW1tYWtr+8DxExERkWEMTkpcXFwwaNAgjbK6uLVhY2ODtm3b4syZMxrlv//+O/z9/QEAbdq0gbW1NXbs2IEhQ4YAAK5du4bs7Gy8//77AIAOHTpApVLh8OHDaNeuHQDg0KFDUKlUUuJCREREpmdwUpKammq0k6vVapw9e1Z6nZubi6ysLLi5uaFRo0aYOXMm4uLi0KVLF3Tv3h0ZGRn45ptvsHv3bgCAQqHAyy+/jMTERDRo0ABubm6YMWMGWrZsiZ49ewKoalnp06cPxowZg5UrVwIAxo4di6ioKDRt2tRo10JEREQP6ZGMBarBrl27BACtR0JCglRn1apVokmTJsLOzk6EhYWJLVu2aByjpKREvPLKK8LNzU3Y29uLqKgocfHiRY06169fF8OHDxdOTk7CyclJDB8+XBQUFBgUK4cEE5kpDgkmMnv6fofKhBDChDnRY6OwsBAKhQIqlQrOzs6mDoeIqhUVAXJ51XO1GrhnkVAiMj19v0Mfapp5IiIiImNhUkJERERmgUkJERERmQW9Rt8sWbJE7wNOnjz5gYMhIiIiy6VXR9eAgAD9DiaT4T//+c9DB2WO2NGVyEyxoyuR2dP3O1SvlpLc3FyjBUZERESkC/uUEBERkVkweEZXALh8+TK2bt2KixcvoqysTGPbwoULjRIYERERWRaDk5LMzEzExMQgICAAZ86cQYsWLXD+/HkIIfDMM8/URYxERERkAQy+fZOUlITExERkZ2fDzs4OmzZtwqVLl9C1a1e88MILdREjERERWQCDk5KcnBwkJCQAAOrXr4+SkhLI5XK8/fbbmD9/vtEDJCIiIstgcFLi6OiI0tJSAICPjw/OnTsnbfv777+NFxkRERFZFIP7lISHh+Onn35CSEgInn/+eSQmJuLXX3/F5s2bER4eXhcxEhERkQUwOClZuHAh1Go1AGD27NlQq9X48ssv0aRJEyxatMjoARIREZFl0GtGV+KMrkRmizO6Epk9fb9DDe5T0rhxY1y/fl2r/ObNm2jcuLGhhyMiIiIC8ABJyfnz51FRUaFVXlpaiitXrhglKCIiIrI8evcp2bp1q/T8+++/h0KhkF5XVFQgMzMTSqXSqMEREd3X3T+S9u4FevcGrKxMFw8RPTC9+5TUq1fVqCKTyXDvLtbW1lAqlUhOTkZUVJTxozQD7FNCZIY2bwYmTQKuXv1fma8vsHgxEBtruriISINRVwkGgMrKSgBAQEAAjhw5And394ePkojoQW3eDAweDNz7u+rKlaryjRuZmBA9Zjj6Rk9sKSEyIxUVgFIJXL6se7tMVtVikpvLWzlEZqDORt8AwJ49exAdHY0mTZrg6aefRkxMDPbt2/fAwRIRGWTfvpoTEqCq9eTSpap6RPTYMDgpWbt2LXr27AkHBwdMnjwZr7zyCuzt7dGjRw/83//9X13ESESk6do149YjIrNg8O2b4OBgjB07FtOmTdMoX7hwIT755BPk5OQYNUBzwds3RGZk926ge/f719u1C+jWra6jIaL7qLPbN//5z38QHR2tVR4TE4Pc3FxDD0dEZLjOnav6jMhkurfLZICfX1U9InpsGJyU+Pn5ITMzU6s8MzMTfn5+RgmKiKhWVlZVw351qU5UUlLYyZXoMaP3kOBRo0Zh8eLFSExMxOTJk5GVlYWIiAjIZDLs378faWlpWFzTfxJERMYWG1s17FfXPCUpKRwOTPQY0rtPiZWVFa5duwYPDw+kp6cjOTlZ6j8SHByMmTNnon///nUarCmxTwmRmSosBKpnmP7uO87oSmSG9P0ONWhG17y8PHh4eBgtyMcJkxIiM8VVgonMntFndAWqppin2lVUVKC8vNzUYRCZLWtra1ixJYOIdDAoKQkKCrpvYnLjxo2HCuhxJYRAXl4ebt68aepQiMyei4sLvLy8+EOHiDQYlJTMmTNHY3Vg+p/qhMTDwwMODg78z5ZIByEEiouLkZ+fDwDw9vY2cUREZE4MSkqGDh1qsX1KalNRUSElJA0aNDB1OERmzd7eHgCQn58PDw8P3sohIone85Twl3/NqvuQODg4mDgSosdD9WeF/a+I6G56JyVcTPj+mLgR6YefFSLSRe/bN5WVlXUZBxEREVk4g6eZJ3pQaWlpcHFxMXUYRERkppiUWLj8/HyMGzcOjRo1gq2tLby8vBAZGYmDBw8CqGpm37Jli8HHVSqVSElJ0SiLi4vD77//boSoiYjoSWTQ6Bt68gwaNAjl5eX4/PPP0bhxY/z555/IzMysk/lm7O3tpZEXplRWVgYbGxtTh0FERPdgS4kFu3nzJvbv34/58+eje/fu8Pf3R7t27ZCUlITnn38eSqUSADBw4EDIZDLp9blz59C/f394enpCLpejbdu22Llzp3Tcbt264cKFC5g2bRpkMpnUqVHX7ZuPPvoIgYGBsLGxQdOmTbFmzRqN7TKZDJ9++ikGDhwIBwcHPP3009i6dau0vaKiAi+//DICAgJgb2+Ppk2bai0MOWLECAwYMADz5s2Dj48PgoKC8Pbbb6Nly5Za70mbNm3w5ptvPuhbSkRED4EtJXWgeoIoUzBk4ja5XA65XI4tW7YgPDwctra2GtuPHDkCDw8PpKamok+fPtJ8Emq1Gv369cP/+3//D3Z2dvj8888RHR2NM2fOoFGjRti8eTPCwsIwduxYjBkzpsbzp6enY8qUKUhJSUHPnj3x7bffYuTIkfD19UX37t2lenPmzMH777+PBQsWYOnSpRg+fDguXLgANzc3VFZWwtfXFxs2bIC7uzsOHDiAsWPHwtvbG0OGDJGOkZmZCWdnZ+zYsQNCCLi4uGDOnDk4cuQI2rZtCwA4efIkjh8/jq+++krv95uIiIxIkF5UKpUAIFQqlda2kpIScfr0aVFSUiKEEEKtVgsAJnmo1WqDrmvjxo3C1dVV2NnZiYiICJGUlCROnDghbQcg0tPT73uckJAQsXTpUum1v7+/WLRokUad1NRUoVAopNcRERFizJgxGnVeeOEF0a9fP43zv/7669JrtVotZDKZ2L59e42xTJw4UQwaNEh6nZCQIDw9PUVpaalGvb59+4oJEyZIr6dOnSq6detW+4WSUdz7mXkoarUQQNXDwL9/Ino0avsOvRtv31i4QYMG4erVq9i6dSsiIyOxe/duPPPMM0hLS6txn6KiIsyaNQshISFwcXGBXC7Hb7/9hosXLxp07pycHHTs2FGjrGPHjsjJydEoCw0NlZ47OjrCyclJmqYcAFasWIFnn30WDRs2hFwuxyeffKIVS8uWLbX6kYwZMwZffPEFbt++jfLycqxbtw6jRo0y6BqIiMh4ePumDjg4OECtVpvs3Iays7NDr1690KtXL7z55psYPXo03nrrLYwYMUJn/ZkzZ+L777/HBx98gCZNmsDe3h6DBw9GWVmZwee+91aTEEKrzNraWmuf6nlzNmzYgGnTpiE5ORkdOnSAk5MTFixYgEOHDmns46hjOfvo6GjY2toiPT0dtra2KC0txaBBgwy+BiIiMg4mJXVAJpPp/BJ8XISEhEjDgK2trVFRUaGxfd++fRgxYgQGDhwIoKqPyfnz5zXq2NjYaO13r+DgYOzfvx8vvfSSVHbgwAEEBwfrHeu+ffsQERGBiRMnSmXnzp3Ta9/69esjISEBqampsLW1xdChQ7lUABGRCTEpsWDXr1/HCy+8gFGjRiE0NBROTk44evQo3n//ffTv3x9A1XwjmZmZ6NixI2xtbeHq6oomTZpg8+bNiI6OhkwmwxtvvKE1469SqcTevXsxdOhQ2Nrawt3dXev8M2fOxJAhQ/DMM8+gR48e+Oabb7B582aNkTz306RJE6xevRrff/89AgICsGbNGhw5cgQBAQF67T969GgpCfrpp5/0Pi8RERkf+5RYMLlcjvbt22PRokXo0qULWrRogTfeeANjxozBsmXLAADJycnYsWMH/Pz80Lp1awDAokWL4OrqioiICERHRyMyMhLPPPOMxrHffvttnD9/HoGBgWjYsKHO8w8YMACLFy/GggUL0Lx5c6xcuRKpqano1q2b3tcwfvx4xMbGIi4uDu3bt8f169c1Wk3u5+mnn0ZERASaNm2K9u3b670fEREZn0wIrrSnj8LCQigUCqhUKjg7O2tsu337NnJzcxEQEAA7OzsTRUgPQgiBZs2aYdy4cZg+fbqpw7EYRv3MFBUBcnnVc7UaeIxvnRI9qWr7Dr0bb9+QxcrPz8eaNWtw5coVjBw50tThEBFZPCYlZLE8PT3h7u6Ojz/+GK6urqYOh4jI4jEpIYvFO5dEROaFHV2JiIjILDApISIiIrPApISIiIjMApMSIiIiMgtMSoiIiMgsMCkhIiIis2DSpGTv3r2Ijo6Gj48PZDKZtAhctREjRkAmk2k8wsPDNeqUlpZi0qRJcHd3h6OjI2JiYnD58mWNOgUFBYiPj4dCoYBCoUB8fDxu3rxZx1dH5kipVCIlJcXUYTyU69evw8PDQ2sRRHMxePBgLFy40NRhENFjyKRJSVFREcLCwqR1VnTp06cPrl27Jj2+++47je1Tp05Feno61q9fj/3790OtViMqKkpjhdphw4YhKysLGRkZyMjIQFZWFuLj4+vsuh4X1Unfe++9p1G+ZcsWyGSyRxLDpk2b0L59eygUCjg5OaF58+ZITEyUts+ePRutWrUy+LhpaWlwcXHRKj9y5AjGjh37EBGb3rx58xAdHQ2lUimVTZkyBW3atIGtrW2N79eGDRvQqlUrODg4wN/fHwsWLNDYvnv3bq0fATKZDL/99ptUJy0tTWed27dvS3XefPNNzJ07F4WFhUa9biJ68pl08rS+ffuib9++tdaxtbWFl5eXzm0qlQqrVq3CmjVr0LNnTwDA2rVr4efnh507dyIyMhI5OTnIyMjAzz//LC249sknn6BDhw44c+YMmjZtatyLeszY2dlh/vz5GDdu3COf1XTnzp0YOnQo3n33XcTExEAmk+H06dPIzMyss3PWtDjgo1ZeXg5ra2uD9yspKcGqVau0knMhBEaNGoVDhw7h5MmTWvtt374dw4cPx9KlS9G7d2/k5ORg9OjRsLe3xyuvvKJR98yZMxprU9z7njk7O+PMmTMaZXevXxMaGgqlUol169ZhwoQJBl8jEVkus+9Tsnv3bnh4eCAoKAhjxoxBfn6+tO3YsWMoLy9H7969pTIfHx+0aNECBw4cAAAcPHgQCoVCYwXY8PBwKBQKqY4upaWlKCws1Hg8iXr27AkvLy/Mmzev1nqbNm1C8+bNYWtrC6VSieTkZI3tSqUS7777LkaNGgUnJyc0atQIH3/8ca3H/Pbbb9GpUyfMnDkTTZs2RVBQEAYMGIClS5cCqPpVPmfOHJw4cUL6RZ6WlgYAWLhwIVq2bAlHR0f4+flh4sSJUKvVAKr+ZkaOHAmVSiXtN3v2bCnOu2/fXLx4Ef3794dcLoezszOGDBmCP//8U9pe3VKzZs0aKJVKKBQKDB06FLdu3ZLqZGRkoFOnTnBxcUGDBg0QFRWFc+fOSdvPnz8PmUyGDRs2oFu3brCzs8PHH38MZ2dnbNy4UeM9+eabb+Do6Khx/Ltt374d9evXR4cOHTTKlyxZgn/+859o3Lixzv3WrFmDAQMGYPz48WjcuDGef/55vPrqq5g/f77WzLYeHh7w8vKSHlZWVhrbZTKZxnZdPxpiYmLwxRdf6IyFiKgmZp2U9O3bF+vWrcOPP/6I5ORkHDlyBM899xxKS0sBAHl5ebCxsdH6he/p6Ym8vDypjoeHh9axPTw8pDq6zJs3T+qDolAo4Ofnp3/gQlStXGqKh4FTp1tZWeHdd9/F0qVLtfriVDt27BiGDBmCoUOH4tdff8Xs2bPxxhtvSAlCteTkZDz77LM4fvw4Jk6ciAkTJmg0/d/Ly8sLp06dQnZ2ts7tcXFxSExMRPPmzaXbd3FxcQCAevXqYcmSJcjOzsbnn3+OH3/8EbNmzQIAREREICUlBc7OztJ+M2bM0Dq+EAIDBgzAjRs3sGfPHuzYsQPnzp2TzlHt3Llz2LJlC7799lt8++232LNnj8Ytr6KiIkyfPh1HjhxBZmYm6tWrh4EDB6KyslLjOK+++iomT56MnJwcDBw4EEOHDkVqaqpGndTUVAwePBhOTk4635O9e/fi2WefrfE9rUlpaanWarz29va4fPkyLly4oFHeunVreHt7o0ePHti1a5fWsdRqNfz9/eHr64uoqCgcP35cq067du1w+PBh6bNKRKQXYSYAiPT09FrrXL16VVhbW4tNmzYJIYRYt26dsLGx0arXs2dPMW7cOCGEEHPnzhVBQUFadZo0aSLmzZtX47lu374tVCqV9Lh06ZIAIFQqlVbdkpIScfr0aVFSUlJVoFYLUZUePPqHWl3re3i3hIQE0b9/fyGEEOHh4WLUqFFCCCHS09PF3X8aw4YNE7169dLYd+bMmSIkJER67e/vL1588UXpdWVlpfDw8BAfffRRjedXq9WiX79+AoDw9/cXcXFxYtWqVeL27dtSnbfeekuEhYXd91o2bNggGjRoIL1OTU0VCoVCq56/v79YtGiREEKIH374QVhZWYmLFy9K20+dOiUAiMOHD0vnd3BwEIWFhRrX3r59+xpjyc/PFwDEr7/+KoQQIjc3VwAQKSkpGvUOHTokrKysxJUrV4QQQvz111/C2tpa7N69u8Zj9+/fX/p30qWm92vlypXCwcFB7Ny5U1RUVIgzZ86IZs2aCQDiwIEDQgghfvvtN/Hxxx+LY8eOiQMHDogJEyYImUwm9uzZIx3n4MGDYs2aNSIrK0vs3btXDBo0SNjb24vff/9d43wnTpwQAMT58+d1xqn1mXkYd3/eDPj7J6JHR6VS1fgdejezbim5l7e3N/z9/fHHH38AqPqlXVZWhoKCAo16+fn58PT0lOrc3Rxf7a+//pLq6GJrawtnZ2eNx5Ns/vz5+Pzzz3H69GmtbTk5OejYsaNGWceOHfHHH39odCgODQ2Vnlc38Vffbuvbty/kcjnkcjmaN28OAHB0dMS2bdtw9uxZvP7665DL5UhMTES7du1QXFxca7y7du1Cr1698NRTT8HJyQkvvfQSrl+/jqKiIr2vOScnB35+fhqtYCEhIXBxcUFOTo5UplQqNVouvL29NW4jnjt3DsOGDUPjxo3h7OyMgIAAAFW3hu52bwtHu3bt0Lx5c6xevRpA1S2WRo0aoUuXLjXGXFJSotXioY8xY8bglVdeQVRUFGxsbBAeHo6hQ4cCgHR7pmnTphgzZgyeeeYZdOjQAcuXL8fzzz+PDz74QDpOeHg4XnzxRYSFhaFz587YsGEDgoKCpFtu1ezt7QHgvv+ORER3e6ySkuvXr+PSpUvw9vYGALRp0wbW1tbYsWOHVOfatWvIzs5GREQEAKBDhw5QqVQ4fPiwVOfQoUNQqVRSHaNzcADUatM8HBweKOQuXbogMjISr732mtY2IYTWaByh4zbRvR03ZTKZdAvj008/RVZWFrKysrQ6aQYGBmL06NH49NNP8csvv+D06dP48ssva4z1woUL6NevH1q0aIFNmzbh2LFj+PDDDwFUdSDVl67r0lVe23UBQHR0NK5fv45PPvkEhw4dwqFDhwAAZWVlGvs5OjpqnWv06NHSLZzU1FSMHDmy1pFP7u7uWkm4PmQyGebPnw+1Wo0LFy4gLy8P7dq1AwCNUTz3Cg8Pl34E6FKvXj20bdtWq86NGzcAmE/HYiJ6PJh09I1arcbZs2el17m5ucjKyoKbmxvc3Nwwe/ZsDBo0CN7e3jh//jxee+01uLu7Y+DAgQAAhUKBl19+GYmJiWjQoAHc3NwwY8YMtGzZUhqNExwcjD59+mDMmDFYuXIlAGDs2LGIioqqu5E3Mhmg4wvI3L333nto1aoVgoKCNMpDQkKwf/9+jbIDBw4gKChIqxNkTZ566im96imVSjg4OEgtHjY2NhqtMQBw9OhR3LlzB8nJyahXryqv3rBhg0YdXfvdKyQkBBcvXsSlS5ek1pLTp09DpVIhODhYr3ivX7+OnJwcrFy5Ep07dwYArfeqNi+++CJmzZqFJUuW4NSpU0hISKi1fuvWrbF27Vq9j38vKysr6d/iiy++QIcOHXT2uap2/Phx6UeALkIIZGVloWXLlhrl2dnZ8PX1hbu7+wPHSkSWx6RJydGjR9G9e3fp9fTp0wEACQkJ+Oijj/Drr79i9erVuHnzJry9vdG9e3d8+eWXGk3pixYtQv369TFkyBCUlJSgR48eSEtL0/iyXLduHSZPniyN0omJial1bhRL1bJlS2nY6N0SExPRtm1bvPPOO4iLi8PBgwexbNkyLF++/KHON3v2bBQXF6Nfv37w9/fHzZs3sWTJEpSXl6NXr14AqpKU6mTV19cXTk5OCAwMxJ07d7B06VJER0fjp59+wooVKzSOrVQqoVarkZmZibCwMDg4OMDhnlaknj17IjQ0FMOHD0dKSgru3LmDiRMnomvXrnp3JnV1dUWDBg3w8ccfw9vbGxcvXsS//vUvvd8DV1dXxMbGYubMmejduzd8fX1rrR8ZGYmkpCQUFBRodPA+e/Ys1Go18vLyUFJSgqysLABViZeNjQ3+/vtvbNy4Ed26dcPt27eRmpqKr776Cnv27JGOkZKSAqVSiebNm6OsrAxr167Fpk2bsGnTJqnOnDlzEB4ejqeffhqFhYVYsmQJsrKypJaqavv27dMYFUdEpJe6797yZKitk45RO+09Qnd3dK12/vx5YWtrK+7909i4caMICQkR1tbWolGjRmLBggUa2+/uQFotLCxMvPXWWzWe/8cffxSDBg0Sfn5+wsbGRnh6eoo+ffqIffv2SXVu374tBg0aJFxcXAQAkZqaKoQQYuHChcLb21vY29uLyMhIsXr1agFAFBQUSPuOHz9eNGjQQACQ4rg3zgsXLoiYmBjh6OgonJycxAsvvCDy8vKk7bo6ji5atEj4+/tLr3fs2CGCg4OFra2tCA0NFbt379bouF3d0fX48eM634fMzEwBQGzYsKHG9+pu4eHhYsWKFRplXbt2FQC0Hrm5uUKIqk604eHhwtHRUTg4OIgePXqIn3/+WeMY8+fPF4GBgcLOzk64urqKTp06iW3btmnUmTp1qmjUqJGwsbERDRs2FL1795Y6ylYrKSkRzs7O4uDBgzVeAzu6ElkWfTu6yoQwcAyphSosLIRCoYBKpdLq9Hr79m3k5uYiICDggTohkmVbt24dpkyZgqtXr8LGxua+9b/77jvMmDED2dnZ0u0rc/Lhhx/i66+/xg8//FBjHaN+ZoqKALm86rla/VjeOiV60tX2HXo3k96+IbJkxcXFyM3Nxbx58zBu3Di9EhIA6NevH/744w9cuXLFsPlzHhFra2utW4BERPowv59ZRBbi/fffR6tWreDp6YmkpCSD9p0yZYpZJiRAVUdyS1++gYgeDJMSIhOZPXs2ysvLkZmZCXn17QciIgvGpISIiIjMApMSIiIiMgtMSoiIiMgsMCkhIiIis8CkhIiIiMwCkxIiIiIyC0xKyKIolUqkpKSYOoyHcv36dXh4eOD8+fOmDkWnwYMHY+HChaYOg4geQ0xKLNiIESMgk8nw3nvvaZRv2bIFMpnskcSwadMmtG/fHgqFAk5OTmjevDkSExOl7bNnz0arVq0MPm5aWhpcXFy0yo8cOYKxY8c+RMSmN2/ePERHR0OpVAIATpw4gX/84x/w8/ODvb09goODsXjxYo19bt++jREjRqBly5aoX78+BgwYoPPYe/bsQZs2bWBnZ4fGjRtrLXQIVP2bhYSEwNbWFiEhIUhPT9fY/uabb2Lu3LkoLCw0yvUSkeVgUmLh7OzsMH/+fBQUFDzyc+/cuRNDhw7F4MGDcfjwYRw7dgxz585FWVlZnZ2zYcOGWqsFm0J5efkD7VdSUoJVq1Zh9OjRUtmxY8fQsGFDrF27FqdOncK///1vJCUlaayEXVFRAXt7e0yePBk9e/bUeezc3Fz069cPnTt3xvHjx/Haa69h8uTJGqsEHzx4EHFxcYiPj8eJEycQHx+PIUOG4NChQ1Kd0NBQKJVKrFu37oGukYgs2CNZHvAJ8KSuEhwVFSWaNWsmZs6cKZWnp6fXuEqwjY2N8Pf3Fx988IHGdn9/fzF37lwxcuRIIZfLhZ+fn1i5cmWt558yZYro1q1bjdtTU1O1Vr2tXiU4OTlZtGjRQjg4OAhfX18xYcIEcevWLSGEELt27dLa72FXCV69erXw9/cXzs7OIi4uThQWFkp1tm/fLjp27CgUCoVwc3MTzz//vDh79qy0vXqV4C+//FJ07dpV2NraimXLlgknJyfx1VdfaVzz1q1bhYODg8bx77Zp0ybh7u5e6/sqhBATJ04U3bt317lN1+rQQggxa9Ys0axZM42ycePGifDwcOn1kCFDRJ8+fTTqREZGiqFDh2qUzZ49W3Tu3LnG+LhKMJFl0XeVYLaU1AEhBIrKikzyEAYu+mxlZYV3330XS5cuxeXLl3XWOXbsGIYMGYKhQ4fi119/xezZs/HGG28gLS1No15ycjKeffZZHD9+HBMnTsSECRPw22+/1XhuLy8vnDp1CtnZ2Tq3x8XFITExEc2bN8e1a9dw7do1xMXFAQDq1auHJUuWIDs7G59//jl+/PFHzJo1CwAQERGBlJQUODs7S/vNmDFD6/hCCAwYMAA3btzAnj17sGPHDpw7d046R7Vz585hy5Yt+Pbbb/Htt99iz549Gre8ioqKMH36dBw5cgSZmZmoV68eBg4ciMrKSo3jvPrqq5g8eTJycnIwcOBADB06FKmpqRp1UlNTMXjwYDg5Oel8T/bu3Ytnn322xve0mkqlgpub233r3e3gwYPo3bu3RllkZCSOHj0qtezUVOfAgQMaZe3atcPhw4dRWlpqUAxEZNm4SnAdKC4vhnyeadYyUSep4Whj2NLtAwcORKtWrfDWW29h1apVWtsXLlyIHj164I033gAABAUF4fTp01iwYAFGjBgh1evXrx8mTpwIoOoLeNGiRdi9ezeaNWum87yTJk3Cvn370LJlS/j7+yM8PBy9e/fG8OHDYWtrC3t7e8jlctSvXx9eXl4a+06dOlV6HhAQgHfeeQcTJkzA8uXLYWNjA4VCAZlMprXf3Xbu3ImTJ08iNzdXWtxuzZo1aN68OY4cOYK2bdsCACorK5GWliYlCvHx8cjMzMTcuXMBAIMGDdI47qpVq+Dh4YHTp0+jRYsWGjHHxsZKr0ePHo2IiAhcvXoVPj4++Pvvv/Htt99ix44dNcZ8/vx5+Pj41LgdqEocNmzYgG3bttVa7155eXnw9PTUKPP09MSdO3fw999/w9vbu8Y6eXl5GmVPPfUUSktLkZeXB39/f4PiICLLxZYSAgDMnz8fn3/+OU6fPq21LScnBx07dtQo69ixI/744w9UVFRIZaGhodLz6oQgPz8fANC3b1/I5XLI5XI0b94cAODo6Iht27bh7NmzeP311yGXy5GYmIh27dqhuLi41nh37dqFXr164amnnoKTkxNeeuklXL9+HUVFRXpfc05ODvz8/DRW2w0JCYGLiwtycnKkMqVSqdFy4e3tLV0XUNWSMmzYMDRu3BjOzs4ICAgAAFy8eFHjfPe2cLRr1w7NmzfH6tWrAVQlRI0aNUKXLl1qjLmkpAR2dnY1bj916hT69++PN998E7169art8nW6t4Nzdcvb3eW66txbZm9vDwD3/XckIrobW0rqgIO1A9RJapOd+0F06dIFkZGReO211zRaPwDdXzq6bhNZW1trvJbJZNItjE8//RQlJSU66wUGBiIwMBCjR4/Gv//9bwQFBeHLL7/EyJEjdcZ64cIF9OvXD+PHj8c777wDNzc37N+/Hy+//LJBHUh1XZeu8tquCwCio6Ph5+eHTz75BD4+PqisrESLFi20Ouw6Omq3YI0ePRrLli3Dv/71L6SmpmLkyJG1jnxyd3evsVPy6dOn8dxzz2HMmDF4/fXXazxGTby8vLRaPPLz81G/fn00aNCg1jr3tp7cuHEDQFXHYiIifTEpqQMymczgWyjm4L333kOrVq0QFBSkUR4SEoL9+/drlB04cABBQUGwsrLS69hPPfWUXvWUSiUcHBykFg8bGxuN1hgAOHr0KO7cuYPk5GTUq1fV2LdhwwaNOrr2u1dISAguXryIS5cuSa0lp0+fhkqlQnBwsF7xXr9+HTk5OVi5ciU6d+4MAFrvVW1efPFFzJo1C0uWLMGpU6eQkJBQa/3WrVtj7dq1WuWnTp3Cc889h4SEBOm2kqE6dOiAb775RqPshx9+wLPPPislZh06dMCOHTswbdo0jToREREa+2VnZ8PX1xfu7u4PFAsRWSYmJSRp2bIlhg8fjqVLl2qUJyYmom3btnjnnXcQFxeHgwcPYtmyZVi+fPlDnW/27NkoLi5Gv3794O/vj5s3b2LJkiUoLy+Xbj0olUrk5uYiKysLvr6+cHJyQmBgIO7cuYOlS5ciOjoaP/30k9Z8GkqlEmq1GpmZmQgLC4ODg4PWUOCePXsiNDQUw4cPR0pKCu7cuYOJEyeia9euenUmBQBXV1c0aNAAH3/8Mby9vXHx4kX861//0vs9cHV1RWxsLGbOnInevXvD19e31vqRkZFISkpCQUEBXF1dAVQlJN27d0fv3r0xffp0qSXDyspKo6Xi9OnTKCsrw40bN3Dr1i1kZWUBgDQPzPjx47Fs2TJMnz4dY8aMwcGDB7Fq1Sp88cUX0jGmTJmCLl26YP78+ejfvz++/vpr7Ny5UysR27dvn1aHWCKi+6rbQUBPjid1SPC9Q0PPnz8vbG1taxwSbG1tLRo1aiQWLFigsf3eobZCCBEWFiYNxdXlxx9/FIMGDRJ+fn7CxsZGeHp6ij59+oh9+/ZJdW7fvi0GDRokXFxcNIYEL1y4UHh7ewt7e3sRGRkpVq9eLQCIgoICad/x48eLBg0aGGVI8N0WLVok/P39pdc7duwQwcHBwtbWVoSGhordu3cLACI9PV0I8b8hwcePH9f5PmRmZgoAYsOGDTW+V3cLDw8XK1as0IgR9wyBBqARY/W166p3t927d4vWrVsLGxsboVQqxUcffaR1/q+++ko0bdpUWFtbi2bNmolNmzZpbC8pKRHOzs7i4MGDNV4DhwQTWRZ9hwTLhDBwDKmFKiwshEKhgEqlgrOzs8a227dvIzc3FwEBAbV2QiTSZd26dZgyZQquXr0KGxub+9b/7rvvMGPGDGRnZ0u3r8zJhx9+iK+//ho//PBDjXWM+pkpKgLk/x3tplYDOvruEJFp1fYdejfeviEykeLiYuTm5mLevHkYN26cXgkJUDX0+o8//sCVK1c0Rg6ZC2tra61bgERE+jC/n1lEFuL9999Hq1at4OnpiaSkJIP2nTJlilkmJAAwduxYNG3a9NGd0MGhqoVEra56TkSPLSYlRCYye/ZslJeXIzMzE3K5aSbbeyLIZFW3bBwdq54T0WOLSQkRERGZBSYlRsQ+w0T64WeFiHRhUmIE1RNLcUptIv1Uf1bunS2XiCwbR98YgZWVFVxcXKT1UBwcHGqdKpzIUgkhUFxcjPz8fLi4uOg9IzARWQYmJUZSvRrt3Qu1EZFuLi4uta7gTESWiUmJkchkMnh7e8PDw8OgReGILI21tTVbSIhIJyYlRmZlZcX/cImIiB4AO7oSERGRWWBSQkRERGaBSQkRERGZBfYp0VP1ZE+FhYUmjoSIiOjxUv3deb+JE5mU6OnWrVsAYLaLoBEREZm7W7duQaFQ1LhdJjjfs14qKytx9epVODk5cWI0IiIiAwghcOvWLfj4+KBevZp7jjApISIiIrPAjq5ERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGbh/wOtQjCvaO0UTwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "fig, ax = plt.subplots()\n", "fig.set_figheight(4)\n", "fig.set_figwidth(6)\n", "\n", "# Stationary fit\n", "plt.plot(\n", " np.array([1, 1]),\n", " np.array(\n", " [\n", " rtnlv_stationary.total_precip_lower.isel(station_num=0),\n", " rtnlv_stationary.total_precip_upper.isel(station_num=0),\n", " ]\n", " ),\n", " \"black\",\n", " label=\"Stationary\",\n", ")\n", "plt.scatter(\n", " np.array([1]),\n", " np.array([rtnlv_stationary.total_precip.isel(station_num=0)]),\n", " c=\"black\",\n", ")\n", "\n", "plt.plot(\n", " np.array([2, 2]),\n", " np.array(\n", " [\n", " rtnlv_non_stationary.total_precip_lower.isel(station_num=0, time=0),\n", " rtnlv_non_stationary.total_precip_upper.isel(station_num=0, time=0),\n", " ]\n", " ),\n", " \"red\",\n", " label=\"Non-Stationary (1955)\",\n", ")\n", "plt.scatter(\n", " np.array([2]),\n", " np.array([rtnlv_non_stationary.total_precip.isel(station_num=0, time=0)]),\n", " c=\"red\",\n", ")\n", "\n", "plt.plot(\n", " np.array([3, 3]),\n", " np.array(\n", " [\n", " rtnlv_non_stationary.total_precip_lower.isel(station_num=0, time=-1),\n", " rtnlv_non_stationary.total_precip_upper.isel(station_num=0, time=-1),\n", " ]\n", " ),\n", " \"green\",\n", " label=\"Non-Stationary (2100)\",\n", ")\n", "plt.scatter(\n", " np.array([3]),\n", " np.array([rtnlv_non_stationary.total_precip.isel(station_num=0, time=-1)]),\n", " c=\"green\",\n", ")\n", "\n", "plt.xlim([0, 4])\n", "plt.xticks([])\n", "plt.ylabel(\"Total precipitation (mm)\")\n", "ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Working with `dask.array` Chunks\n", "\n", "Currently, the Python-to-Julia interaction is not thread-safe. To mitigate potential issues, it is recommended to use the `dask.scheduler=\"processes\"` option when computing results. This ensures that tasks are executed in separate Python processes, providing better isolation and avoiding thread-related conflicts.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:43.238969Z", "start_time": "2024-12-11T20:54:43.188645Z" } }, "outputs": [], "source": [ "ds_c = ds.chunk({\"time\": -1, \"station_num\": 1})\n", "\n", "fit_stationary_c = xhe.fit(\n", " ds_c,\n", " dist=\"genextreme\",\n", " method=\"ml\",\n", " variables=[\"total_precip\"],\n", " confidence_level=0.95,\n", ")\n", "fit_stationary_c = fit_stationary_c.compute(scheduler=\"processes\")\n", "clear_output(wait=False)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 460B\n",
       "Dimensions:             (station_num: 5, dparams: 3)\n",
       "Coordinates:\n",
       "  * station_num         (station_num) int64 40B 1001 1004 1008 1009 1012\n",
       "  * dparams             (dparams) <U5 60B 'shape' 'loc' 'scale'\n",
       "Data variables:\n",
       "    total_precip        (station_num, dparams) float64 120B dask.array<chunksize=(1, 3), meta=np.ndarray>\n",
       "    total_precip_lower  (station_num, dparams) float64 120B dask.array<chunksize=(1, 3), meta=np.ndarray>\n",
       "    total_precip_upper  (station_num, dparams) float64 120B dask.array<chunksize=(1, 3), meta=np.ndarray>
" ], "text/plain": [ " Size: 460B\n", "Dimensions: (station_num: 5, dparams: 3)\n", "Coordinates:\n", " * station_num (station_num) int64 40B 1001 1004 1008 1009 1012\n", " * dparams (dparams) \n", " total_precip_lower (station_num, dparams) float64 120B dask.array\n", " total_precip_upper (station_num, dparams) float64 120B dask.array" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fit_stationary_c" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.10" }, "nbsphinx": { "execute": "never" } }, "nbformat": 4, "nbformat_minor": 4 }