{ "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", "\n", "import xhydro.extreme_value_analysis as xhe\n", "from xhydro.testing.helpers import deveraux" ] }, { "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": [ { "name": "stderr", "output_type": "stream", "text": [ "Downloading file 'extreme_value_analysis/NOAA_GFDL_ESM4.zip' from 'https://raw.githubusercontent.com/hydrologie/xhydro-testdata/v2025.3.31/data/extreme_value_analysis/NOAA_GFDL_ESM4.zip' to 'C:\\Users\\ospin\\AppData\\Local\\xhydro-testdata\\xhydro-testdata\\Cache\\v2025.3.31'.\n", "Unzipping contents of 'C:\\Users\\ospin\\AppData\\Local\\xhydro-testdata\\xhydro-testdata\\Cache\\v2025.3.31\\extreme_value_analysis\\NOAA_GFDL_ESM4.zip' to 'C:\\Users\\ospin\\AppData\\Local\\xhydro-testdata\\xhydro-testdata\\Cache\\v2025.3.31\\extreme_value_analysis\\NOAA_GFDL_ESM4.zip.unzip'\n" ] }, { "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": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "file = deveraux().fetch(\"extreme_value_analysis/NOAA_GFDL_ESM4.zip\", pooch.Unzip())\n", "\n", "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": 3, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:28.202587Z", "start_time": "2024-12-11T20:54:28.105446Z" } }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAE6CAYAAAB9D9Q3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACqzElEQVR4nO29eXxU1f3//5p9sk42khBIWARBZBFREdSKG+JXxEpbtba4YbW1blVs669a7SbWTz8ulWqtC1iXam3Bj1upqLggiCCLoMhm2AkJWSb7rPf3x51z7rl35s7cO5lkZpL38/HgoZncmZy5M/fe1329N4skSRIIgiAIgiCIrMSa7gUQBEEQBEEQyUNijiAIgiAIIoshMUcQBEEQBJHFkJgjCIIgCILIYkjMEQRBEARBZDEk5giCIAiCILIYEnMEQRAEQRBZDIk5giAIgiCILIbEHEEQBEEQRBZDYo4giKR46aWX8Mgjj/ToNR5//HEsWbIk6vE9e/bAYrHE/B3R+6xatQrXXXcdpkyZApfLBYvFgj179uhu/9hjj2Hs2LFwuVwYMWIEfvOb3yAQCERtV19fj6uvvhplZWXIzc3FtGnT8N5770Vt9+abb+LKK6/EhAkT4HA4YLFYUvn2CKLfQWKOIIik6E0xN3jwYKxZswYXXnhhj16fSI733nsP7777LmpqajB9+vS42/7hD3/Arbfeirlz5+K///0vbrzxRtx///346U9/qtrO5/PhnHPOwXvvvYdHH30U//d//4eKigrMmjULH374oWrbZcuW4dNPP8W4ceMwadKklL8/guhvWGg2K0EQyTB79mxs3bo1rmOTiPHjx6OsrAwffPBBytZF9JxwOAyrVb7X/9Of/oQ777wTtbW1GD58uGq7xsZGDB06FFdeeSWefPJJ/vj999+Pu+++G1u3bsW4ceMAyML9pz/9KVavXo1p06YBAILBICZNmoT8/HysXbs25t+/6aab8Je//AV0qSIIfciZIwgiioaGBlx//fWorq6Gy+XCoEGDcNppp+Hdd98FAMyYMQNvvfUW9u7dC4vFwv8xfvOb32Dq1KkoKSlBYWEhTjzxRDzzzDOqC/Lw4cPx5Zdf4sMPP+TPZ2JBL8y6atUqnHPOOSgoKEBubi6mT5+Ot956S7XNkiVLYLFYsHLlSvzkJz9BWVkZSktLMXfuXBw6dMjUfrjvvvtgsVjw5Zdf4vvf/z48Hg8qKipw7bXXwuv18u3ihYUtFgvuu+++qNf84osv8L3vfQ8ejwclJSW4/fbbEQwGsX37dsyaNQsFBQUYPnw4HnzwQVNrTgVMSCVi+fLl6O7uxjXXXKN6/JprroEkSXjttdf4Y8uWLcOYMWO4kAMAu92OH/7wh/jss89w8OBB03+fIAgZe7oXQBBE5jFv3jxs2LABf/jDH3DssceipaUFGzZsQGNjIwDZZbn++uuxe/duLFu2LOr5e/bswQ033ICamhoAwKeffoqbb74ZBw8exK9//WsA8sX9u9/9LjweDx5//HEAgMvl0l3Thx9+iPPOOw8TJ07EM888A5fLhccffxwXXXQR/vGPf+Cyyy5TbX/dddfhwgsvxEsvvYT9+/fjzjvvxA9/+EO8//77pvfHd77zHVx22WWYP38+tmzZgrvuugsA8Oyzz5p+Lcall16KH/7wh7jhhhuwYsUKPPjggwgEAnj33Xdx4403YsGCBXjppZfwi1/8AqNGjcLcuXPjvl4oFDLkXlmt1pSJpa1btwIAJkyYoHp88ODBKCsr479n255xxhlRrzFx4kQAwJdffokhQ4akZF0EMdAgMUcQRBSffPIJrrvuOvzoRz/ij1188cX8/8eNG4eioiK4XC6ceuqpUc9fvHgx//9wOIwZM2ZAkiQ8+uijuOeee2CxWDB58mTk5OSgsLAw5mto+eUvf4ni4mJ88MEHyM/PByCHek844QQsWLAAl156qcodnDVrFv785z/zn5uamvDzn/8cdXV1qKysNLU/5s+fjzvvvBMAcO6552LXrl149tln8cwzzySdnH/99dfj9ttv56/5zjvvYNGiRVi6dCkuueQSALID+uabb+LFF19MKObOOeecqNyzWFx11VUpKyxpbGyEy+VCXl5e1O9KSkq4+GfblpSUxNyO/Z4giOQgMUcQRBSnnHIKlixZgtLSUpx77rmYMmUKHA6H4ee///77uP/++7Fu3Tq0traqfldfX4+KigpT6+no6MDatWvxk5/8hAs5ALDZbJg3bx5+8YtfYPv27Rg7diz/3Zw5c1SvwRygvXv3mhZzsV6ru7s7qffCmD17turn4447Dps3b8YFF1zAH7Pb7Rg1ahT27t2b8PWefPJJtLW1JdyurKzM/GLjEE/Man9nZluCIIxDYo4giCheeeUV/P73v8fTTz+Ne+65B/n5+bjkkkvw4IMPJhRCn332GWbOnIkZM2bgqaeewtChQ+F0OvHaa6/hD3/4A7q6ukyvp7m5GZIkYfDgwVG/q6qqAhDt7JSWlqp+ZiHcZP5+Kl+LoXWpnE4ncnNz4Xa7ox7XCuJYjBo1ynCYNVWUlpaiu7sbnZ2dyM3NVf2uqakJU6ZMUW0by31ramoCEL0/CIIwDmWZEgQRRVlZGR555BHs2bMHe/fuxcKFC7F06VJcffXVCZ/78ssvw+Fw4M0338Sll16K6dOn46STTurReoqLi2G1WnH48OGo37GihlQ7TmZgAszn86ke78vQ4TnnnAOHw5Hw37XXXpuyv8ly5bZs2aJ6vK6uDkePHsX48eNV22q3E58rbksQhDnImSMIIi41NTW46aab8N577+GTTz7hj7tcrpjOlMVigd1uh81m4491dXXh+eefj9pW7zW05OXlYerUqVi6dCn+9Kc/IScnB4Ccj/fCCy9g6NChOPbYY5N5eymhoqICbrcbX3zxherx//u//+uzNaQjzDpr1iy43W4sWbIEU6dO5Y+ziuJvf/vb/LFLLrkEN954I9auXcu3DQaDeOGFFzB16lTusBIEYR4ScwRBqPB6vTjrrLNwxRVXYOzYsSgoKMC6deuwfPlyVRL+hAkTsHTpUjzxxBOYMmUKrFYrTjrpJFx44YV46KGHcMUVV+D6669HY2Mj/vSnP8WsVJ0wYQJefvllvPLKKxg5ciTcbndUZSRj4cKFOO+883DWWWdhwYIFcDqdePzxx7F161b84x//SGvOlcViwQ9/+EM8++yzOOaYYzBp0iR89tlneOmll/psDWPGjEnZazU0NPBiCuac/ec//8GgQYMwaNAgnHnmmQDk0Ojdd9+Ne+65ByUlJZg5cybWrVuH++67D9dddx3vMQcA1157Lf7yl7/ge9/7Hh544AGUl5fj8ccfx/bt23nLG8bevXuxbt06AMDu3bsBAP/6178AyC1teur0EkS/QyIIghDo7u6WfvzjH0sTJ06UCgsLpZycHGnMmDHSvffeK3V0dPDtmpqapO9+97tSUVGRZLFYJPF08uyzz0pjxoyRXC6XNHLkSGnhwoXSM888IwGQamtr+XZ79uyRZs6cKRUUFEgApGHDhkmSJEm1tbUSAGnx4sWqtX388cfS2WefLeXl5Uk5OTnSqaeeKr3xxhuqbRYvXiwBkNatW6d6fOXKlRIAaeXKlYb3xb333isBkBoaGmL+DfG9eL1e6brrrpMqKiqkvLw86aKLLpL27NkjAZDuvffehK951VVXSXl5eVFrOPPMM6Xjjz/e8JpTAdtXsf6deeaZUds/+uij0rHHHis5nU6ppqZGuvfeeyW/3x+1XV1dnXTllVdKJSUlktvtlk499VRpxYoVUdux/Rvr31VXXdUL75ggshuaAEEQBEEQBJHFUAEEQRAEQRBEFkM5cwRBDDjC4TDC4XDcbex2Oj0SBJEdkDNHEMSA47e//W3CFh579uxJ9zIJgiAMQTlzBEEMOA4dOsT70+kxceJEOJ3OPloRQRBE8pCYIwiCIAiCyGIozEoQBEEQBJHFUIavQcLhMA4dOoSCggIaCE0QBEEQRK8iSRLa2tpQVVWVcKYyiTmDHDp0CNXV1eleBkEQBEEQA4j9+/dj6NChcbchMWeQgoICAPJOLSwsTPNqCIIgCILoz7S2tqK6uprrj3iQmDMIC60WFhaSmCMIgiAIok8wktpFBRAEQRAEQRBZDIk5giAIgiCILIbEHEEQBEEQRBZDYo4gCIIgCCKLITFHEARBEASRxZCYIwiCIAiCyGJIzBEEQRDEAGV7XRvmPbMWG/Y1p3spRA+gPnMEQRAEMUB5ffNBfLzzKIYW5+LEmuJ0L4dIEnLmCIIgCGKA0toVBAB4u/xpXgnRE9Iq5j766CNcdNFFqKqqgsViwWuvvRa1zbZt2zBnzhx4PB4UFBTg1FNPxb59+/jvfT4fbr75ZpSVlSEvLw9z5szBgQMHVK/R3NyMefPmwePxwOPxYN68eWhpaenld0cQBEEQmU27TxZzLZ2BNK+E6AlpFXMdHR2YNGkSFi1aFPP3u3fvxumnn46xY8figw8+wObNm3HPPffA7XbzbW677TYsW7YML7/8MlatWoX29nbMnj0boVCIb3PFFVdg06ZNWL58OZYvX45NmzZh3rx5vf7+CIIgCCKTaesmMdcfsEiSJKV7EYA8e2zZsmX49re/zR+7/PLL4XA48Pzzz8d8jtfrxaBBg/D888/jsssuAwAcOnQI1dXVePvtt3H++edj27ZtGDduHD799FNMnToVAPDpp59i2rRp+PrrrzFmzBhD62ttbYXH44HX66XZrARBEES/4PK/rcGn3zRhSFEOPvnl2eleDiFgRndkbM5cOBzGW2+9hWOPPRbnn38+ysvLMXXqVFUo9vPPP0cgEMDMmTP5Y1VVVRg/fjxWr14NAFizZg08Hg8XcgBw6qmnwuPx8G1i4fP50NraqvpHEARBEP0JFmb1dpEzl81krJirr69He3s7HnjgAcyaNQvvvPMOLrnkEsydOxcffvghAKCurg5OpxPFxeoKnIqKCtTV1fFtysvLo16/vLycbxOLhQsX8hw7j8eD6urqFL47giAIgkg/7ZEwa7sviEAonObVEMmSsWIuHJa/VBdffDF+9rOf4YQTTsAvf/lLzJ49G3/961/jPleSJFgsFv6z+P9622i566674PV6+b/9+/cn+U4IgiAIIjNhOXMA0EruXNaSsWKurKwMdrsd48aNUz1+3HHH8WrWyspK+P1+NDermx3W19ejoqKCb3PkyJGo129oaODbxMLlcqGwsFD1jyAIgiD6E20+Rcy1kJjLWjJWzDmdTpx88snYvn276vEdO3Zg2LBhAIApU6bA4XBgxYoV/PeHDx/G1q1bMX36dADAtGnT4PV68dlnn/Ft1q5dC6/Xy7chCIIgiIGGLxiCP6iEVqmiNXtJ6wSI9vZ27Nq1i/9cW1uLTZs2oaSkBDU1Nbjzzjtx2WWX4Vvf+hbOOussLF++HG+88QY++OADAIDH48H8+fNxxx13oLS0FCUlJViwYAEmTJiAc889F4Ds5M2aNQs/+tGP8OSTTwIArr/+esyePdtwJStBEARB9Dc6fCHVzxRmzV7SKubWr1+Ps846i/98++23AwCuuuoqLFmyBJdccgn++te/YuHChbjlllswZswY/Pvf/8bpp5/On/Pwww/Dbrfj0ksvRVdXF8455xwsWbIENpuNb/Piiy/illtu4VWvc+bM0e1tRxAEQRADgXYhXw4AWmgKRNaSMX3mMh3qM0cQBEH0J7485MWFf17Ff77vonG4+rQRaVwRIdIv+swRBEEQBNF7tEU5cxRmzVZIzBEEQRDEACQqzEoFEFkLiTmCIAiCGIC0+9RijgogshcScwRBEAQxAGnzUZi1v0BijiAIgiAGICzMWuCSG1u0dFI1a7ZCYo4gCIIgBiDtPtmJG1qSCwDwkjOXtZCYIwiCIIgBCKtmHVKUA4DEXDZDYo4gCIIgBiAszDq0WBZzLZ0BUOvZ7ITEHEEQBEEMQFgBBBNzwbCETn8o3lOIDIXEHEEQBEEMQJgzN6jABadNlgNU0ZqdkJgjCIIgiAEI6zNX4LbDk+sAAHipcXBWQmKOIAiCIAYgiphzwJMji7mWLmpPko2QmCMIgiD6lHBYwi/+9QX+9tHudC9lQNPWLbtw+S47inLImctmSMwRBEEQfcqXh1rxyvr9WPT+rnQvZUDDWpPku+woYmFWypnLSkjMEQRBEH3KN0fbAQDdgXCaVzJw8QfD8AXl/V/gtqOQh1lJzGUjJOYIgiCIPuWbhg4AgD8URjhMfc3SQYcwlzXPZUdRjhOA3GuOyD5IzBEEQRB9yjdHO/j/M3eI6FtY8UOOwwaHzUph1iyHxBxBEATRp9RGwqwA4AtSk9p0wPPl3HYA4NWsXqpmzUpIzBEEQRB9hiRJqG0gZy7dsErWApcs5siZy27sRja6/fbbTb/w3XffjZKSEtPPIwiCIPov9W0+dAgjo3xUBJEWWJhV68xRzlx2YkjMPfLII5g2bRqcTqehF121ahVuuukmEnMEQRCEit0N7aqfuynMmha4mHORmOsPGBJzALBs2TKUl5cb2ragoCDpBREEQRD9l1qh+AEgZy5dsJy5AjcLs8pmTSuFWbMSQzlzixcvhsfjMfyiTz75JCoqKpJeFEEQBNE/+aZBI+bImUsLijMnO3JsAkSbL4hAiAR2tmHImbvqqqtMvegVV1yR1GIIgiCI/k2UM0cFEGmhXePMsabBgOzOlea70rIuIjmompUgCILoM77R5swFyJlLB+JcVgCwWS1c2FFFa/aRMjG3efNm2Gy2VL0cQRAE0c/wB8PY39wFABhWmguAnLl00aapZgWU9iQ00iv7SKkzJ0k0loUgCIKIzb6mToTCEvKcNlQXMzFHzlw6YGFW5swBQuNgqmjNOgxXs86dOzfu771eLywWS48XRBAEQfRPWL7ciEF5cDtkL4GqWdMDK4AoEJ25yHxWCrNmH4bF3BtvvIHzzjtPt0o1FKK7K4IgCEIfli83oiwf4bAcyaGcufQQS8x5WJi1k0Z6ZRuGxdxxxx2H73znO5g/f37M32/atAlvvvlmyhZGEARB9C+4M1eWhwPNnQAoZy5dKGFWpYpVmc8aTMuaiOQxnDM3ZcoUbNiwQff3LpcLNTU1KVkUQRAE0f9gPeaOGZQHl10umCMxlx5aY+TMsV5zLV3kzGUbhp25v/71r3FDqccddxxqa2tTsiiCIAii//GN4Mxt3NcCgAog0kW7T86LU4VZqQAiazEs5lwuaiBIEARBJEdrdwBH230AZDHnihRAdFMBRJ8TCIX5fi+I0ZqECiCyjx61Jrnwwgtx+PDhVK2FIAiC6KfURkKsgwpcKHA74OZhVnLm+poOn5ITl6dqTSJXs1KfueyjR2Luo48+QldXV6rWQhAEQfRTvjnKKlnzAIA7c9SapO9pi+TLuR1WOGyKDGBhVqpmzT5onBdBEATR6zBnbiQTc1QAkTZYWxKxkhUQw6xUzZpt9EjMDRs2DA6HI/GGBEEQxICmIZIvN9iTAwBw2VnOHIVZ+xrmzIn5cgBQmieHWZs7/Xx2K5EdmBZzYsXq1q1bUV1dnfQf/+ijj3DRRRehqqoKFosFr732mu62N9xwAywWCx555BHV4z6fDzfffDPKysqQl5eHOXPm4MCBA6ptmpubMW/ePHg8Hng8HsybNw8tLS1Jr5sgCIIwR3OHLA6K82QDwO0gZy5dsEpWsS0JIOczjhyUh1BYwvtf16djaUSSmBZzo0aNwllnnYUXXngB3d3dPfrjHR0dmDRpEhYtWhR3u9deew1r165FVVVV1O9uu+02LFu2DC+//DJWrVqF9vZ2zJ49W9VG5YorrsCmTZuwfPlyLF++HJs2bcK8efN6tHaCIAjCOM2RPKyiXNn9Yc4cFUD0PXrOnMViwQXjKwEA/9lS1+frIpLHtJjbvHkzJk+ejDvuuAOVlZW44YYb8NlnnyX1xy+44AL8/ve/jzv39eDBg7jpppvw4osvRoV0vV4vnnnmGfzv//4vzj33XEyePBkvvPACtmzZgnfffRcAsG3bNixfvhxPP/00pk2bhmnTpuGpp57Cm2++ie3btye1boIgCMIcLZHeZcWRvCxFzJEz19coOXPR3ckuGD8YAPDBjnp0+il3LlswLebGjx+Phx56CAcPHsTixYtRV1eH008/HccffzweeughNDQ0pGxx4XAY8+bNw5133onjjz8+6veff/45AoEAZs6cyR+rqqrC+PHjsXr1agDAmjVr4PF4MHXqVL7NqaeeCo/Hw7eJhc/nQ2trq+ofQRBENiJJUrqXwJ25YubMRcKsA7XP3KGWLvz+za+wv6mzz/82H+XljhZzx1cVorokB92BMD7cnrrreV8QDktYvesoWgdgvl/SBRB2ux2XXHIJ/vnPf+KPf/wjdu/ejQULFmDo0KG48sorU9J/7o9//CPsdjtuueWWmL+vq6uD0+lEcXGx6vGKigrU1dXxbcrLy6OeW15ezreJxcKFC3mOncfj6VFuIEEQRLpY8dURnPi7Ffhge/pyoCRJ4s4cq5h0D/Aw6yvr9uPpVbV44dO9ff63mTNXEMOZk0Otsjv3n63ZFWp95L2duOLptZjxPx/g+U/3IhgaODcKSYu59evX48Ybb8TgwYPx0EMPYcGCBdi9ezfef/99HDx4EBdffHGPFvb555/j0UcfxZIlS2CxWEw9V5Ik1XNiPV+7jZa77roLXq+X/9u/f7+pNRAEQWQCq3Y2oLkzgE92HU3bGjr9IfgjF1atMzdQ+8wx94jNSO1L2uI4cwAwK5I39962I1lTbeztDODZVXKBZlOHH/e8thUXPPox1uxuTPPK+gbTYu6hhx7ChAkTMH36dBw6dAh///vfsXfvXvz+97/HiBEjcNppp+HJJ5/Ehg0berSwjz/+GPX19aipqYHdbofdbsfevXtxxx13YPjw4QCAyspK+P1+NDc3q55bX1+PiooKvs2RI0eiXr+hoYFvEwuXy4XCwkLVP4IgiGyD5aR1pfGizEKsTpsVuU5ZxA30nDn2vtPhTCoFELFbi50wtAiVhW50+ENYtTN9NwFmeG7NHrT7ghhTUYDfzDkeRbkO7Kxvx3XPrcuINIPexrSYe+KJJ3DFFVdg3759eO211zB79mxYreqXqampwTPPPNOjhc2bNw9ffPEFNm3axP9VVVXhzjvvxH//+18AwJQpU+BwOLBixQr+vMOHD2Pr1q2YPn06AGDatGnwer2qIo21a9fC6/XybQiCIPorXMz50yeaxBAri4gM9GpW5nj50yBm9VqTMKxWC3fnsiHU2uEL4tlPZFfup2ePwlXTh2PFz86Uf+cPocPf/79jsT/JOOzcuTPhNk6nE1dddVXC7drb27Fr1y7+c21tLTZt2oSSkhLU1NSgtLRUtb3D4UBlZSXGjBkDAPB4PJg/fz7uuOMOlJaWoqSkBAsWLMCECRNw7rnnAgCOO+44zJo1Cz/60Y/w5JNPAgCuv/56zJ49m78OQRBEf4WJpXSGy7TFD4DQZ26AhlkVZy4dYi52axKRC8ZXYsnqPVjxVR38wQlw2jN3YNSLa/eipTOAEWV5uHCCnO9Xlu+E3WpBMCyhvTuoK1z7C2n9dNavX4/Jkydj8uTJAIDbb78dkydPxq9//WvDr/Hwww/j29/+Ni699FKcdtppyM3NxRtvvAGbzca3efHFFzFhwgTMnDkTM2fOxMSJE/H888+n/P0QBEFkGkwspTfMqi5+ABRnzh8KIxzu/2EwLb7I55EWMdet35qEcdLwEpTlO9HaHcTa2szNO+sOhPDUx7Ir95MZx8BmlZ1fi8XCcwIHwjSLtErVGTNmmIpl79mzJ+oxt9uNxx57DI899pju80pKSvDCCy8ks0SCIIispjvizHWlMdTUEsOZYwUQgCzo3FZb1PP6M6wliy8NIrstTp85hs1qweSaYqz46gj2NnbijNF9tTpzvLp+PxrafBhSlINLJg9R/a7AbUdLZ4C/3/5M5vqmBEEQRI/JCGdOM8oLUJw5YGCGWln4Ox3OXKJqVgYTe5ncPPj5SGuXH585Eg6bWtLku+TvW1saKob7GhJzBEEQ/RgmFjIhZ65IcOYcNisPiXWnsQgiXZWO3JlLg5jzdsni2pMTu5qVkeeS3dJ2X+YWEBxo7gIAfOvYQVG/YzmB7STmCIIgiGyGOUDpdOaUMKtaPPCK1iSdOW9noEci1R8M44JHP8ZPX+pZKy2RZ1bV4sz/WYm9jR1xt1Ocub79XLoDIV5BK4rrWORFnLmODA1TdgdC6IykD5TkRb8X1hR5IOTMmRZzjY2N+OlPf4px48ahrKwMJSUlqn8EQRBE5qC0JsmEAgj1Bbcn7UlauwM4/cH38b2/rkl6XfuaOvB1XRtWfBndizQZJEnCs6tqsbexE2u/aYq7rZIz17fOHGsTY7NakOeMn6eY78xsMcccX4fNEjP/jztzGbr+VGK6AOKHP/whdu/ejfnz56OiosL0dAaCIAii78iEnLlYBRAA4LLbAASSCjXuqGtDW3cQ2w4nPzebhQ9ZRa3V2rPr2b6mThxskcN+LJSpR3eaqllbuiIh7xxHwus3c+YyVQw1dSjfq1jvheUEpmPKRl9jWsytWrUKq1atwqRJk3pjPQRBEEQK6c6IPnORAghNmNXtkJ25ZNbGcqWCYQn+YDipPmhiLpUvGEZOAqcqEZ/sUlp4JBJz6ZoA4Y18Fp7c+PlygFIAkbHOHCus0QkXswkXlDMXg7Fjx6Krq6s31kIQBEGkGObMBUISAmkaPB6rAAJgzlxy7tSB5k7+/8mGkNkkBCA1Ynf1bmX0VeY6c8aKHwAxZy4zCyCamOObF/u95FPOnD6PP/44fvWrX+HDDz9EY2MjWltbVf8IgiCIzECSJJXzkw53LhgK89YQUQUQjuRz5pgzBwAdSbbOEKs0e1pRGw5LqqHuLXHEnPy5yCLOHwz3aUUtE5lFhsQcq2bNTGeLhe9jFT8AQCFvGpyZ608lpsOsRUVF8Hq9OPvss1WPS5IEi8WCUCgzFTxBEMRAIxiWIA5X6AqEdIer9xaiqNG6QT2pZhXFXGeyzly36Mz1zCHbfqQNjZEcLiC+M6d143zBMB9v1tvwMKsBMcfDrBnaZ07MmYtFPhVA6PODH/wATqcTL730EhVAEARBZDBa0dDt7/swK3NPCt122DVNXVmYNRlXLBVhVnEAe09dy9URV85ps8IfCscXc4H0iTleAJGgLQmQ+a1JmjviO3MFvGlw/w+zmhZzW7duxcaNG2lIPUEQRIajHRWVjopWXvwQ44LLCiDMOnPhsMSrRoHkJxSI4bcei7ldcr7ct44tw7vb6tEa15kLxfi5bxxTow2DAcWZ64mz9f7XR/D8mr3443cnorzAnfTrxKKpM34BBJ/NmqFiNJWYzpk76aSTsH///t5YC0EQBJFCujXOXFrEXIe+E5RsAUR9mw+BkBI/TjbMKjpOPQmzBkNhrK2V+8rNGj8YQPwwq/Zv9WWvuRYTYVbmzHUHwggmWTyzZPVerNzegOVb65J6fjzYd0uvAKKAcub0ufnmm3HrrbfizjvvxIQJE+BwqHfixIkTU7Y4giAIInminLk0NA5u0WlLAiTfNFgMsQI9yJkTxVwPCiC+OOhFuy8IT44D044pBSCLOZZLrkX7t/qyopUXQBhoTcIKIAA5JO3JMd/+hf29Qy3dpp+biEQ5c4UDqDWJaTF32WWXAQCuvfZa/pjFYqECCIIgiAwjKmcuLWFW/Quui/eZMydmxOIHIPkwqyjmtMLXDCzEOm1kKRetobCEDn8o5mSC6Jy5vvtczIg5l90Gh82CQEhCR0SsmqUt8vcOe1Pf0qw5QTUr2/ddgRACoTActv47wdS0mKutre2NdRAEQRBJEg5L+PSbRhxf5VE1g9WKuXTmzMUSD0qYtWfOXLLvq707NWFWVvxw2qhS5DgUAeTtCsQUc+l05syEWQE51NrSGUi6CKI1UnxwuBecuXg3CoCSMwfIIXUjRR/ZimkxN2zYsN5YB0EQBJEkH+1swNWL12Hu5CF46LIT+OOZEWZN7MyZzRnTOnPJNrVVhVmTFITdgRDW720GAEw7pgwWiwWeHAeOtvvh7QxgSFFO1HOinLk+zZmTPw9PjjFhk+eUxVwyRRCSJKG1S37eoRQ7c13+EBfges6cw2aF22FFd0DudUhiTsPBgwfxySefoL6+HuGw+kt4yy23pGRhBEEQhDFqj3YAgKrCE8gUZ46JuXjOXHJiLtdpQ6c/hK4kw6wdKRBzB1u64A+Gke+y45hBeQCAQibmdIogtH/L30eTOUJhiVd2GgmzAkqoMpm8xO5AmL+3I63dKZl/y2DTH5w2K3LjjGErcDvQHfD1+yII02Ju8eLF+PGPfwyn04nS0lJVcqfFYiExRxAE0cew0JlWrGlFQ3py5liYNVY1a3KzWVmYdXR5PjYf8CZdACG2rEg21MlEgkcYXM9CmLpiThtm7aPPpa07ADZswmiYNbcHUyBahf5ugZCEo+0+lBempj2JWMkar99tgcuOhjZfv+81Zzob8Ne//jV+/etfw+v1Ys+ePaitreX/vvnmm95YI5ECtte14Q9vfcUPAIIg+g8sdKYNo0Y5cxkWZmWNcs0IKbHH3LEVBQCAziTFUCpak7C8OzE3jgklvV5zsZoG9wVM9Oc5bYaLAfJ70DhY+/4PeVOXN5eokpUxUNqTmBZznZ2duPzyy2G19t+qkP7Io+/twFMf1+L1zYfSvRSCIFJMs44zlxlh1ngFEOZbk7AeczarBSMH5QMAOpMQGqGwpHL0km1NwhyfAne0mDPszPWRmFMqWY3njuU5eyDmNG7Y4ZbU5c0lqmRlDJSRXqYV2fz58/Hqq6/2xlqIXqT2qByWONruS/NKCIJINWz+qTZcqRVJfS3mJElSnLkYF11FzBkXMyzEWlno5gIqmTCrdt5osiFoFqrNNyHm0tWahH1PCk20GMnjUyDMr5EVPzB6xZlLIOYGykgv0zlzCxcuxOzZs7F8+fKYTYMfeuihlC2OSA2SJOFAk3wCZDY7QRD9ByaYtKJGKxr6Omeuwx/ikxpiFkCwMKuJECcrfhhanMMT35MRqdpGssmGWVn4rsCtvD/zOXOpd+aeW70HT3ywGy9cNxWjymUHk31PikyIufxIzlzmOXP6zahFBspIL9Ni7v7778d///tfPptVWwBBZB4tnQH+RW6JM2KGIIjshIWcugIh1dSBdOfMsRxdp92KnBiD5N2sAMKEM8WcuaHFuch1Ju/MacNuyRYhxMuZ069m7f2cueVb61DX2o2PdjRwMddqomEwI68H81m1OXOHU+jMse9WCeXMAUhCzD300EN49tlncfXVV/fCcojeYF+T0mCT3ZkRBNF/YI67JMnCgBUWMCfOYpF/19dhVnGUV6yb/Z44c9UlijOXjGukFSc9zZkrFMKshYnCrFE5c6n/XJgrJrarMdswGFDEXHLOnPycQrcdrd3BlPaaa4oTvhcpYGK0n4s50zlzLpcLp512Wm+sheglRDGnV11FEP2Fgy1deOTdHWg0mR/61aFWfBYZlp5NBENhlesghlKZ48Mu3l192JwWSNyhP5kCCCXMmpsRYVYmCs0UQPRFNSsXc0KDZRaZ8Zhw5ng1axK9/Nj1ZmxlIQCgrjecuURizj0wcuZMi7lbb70Vjz32WG+sheglVM4ciTmin/P0x9/gkXd34vlP9xp+jj8Yxvef+hRXPPVp1rXv0R7TXSoxJ/8/y5Hq65w5Jub0wno9KYCQc+Z6UAChdeaSLYBIpjVJH+TMseID0Znj1awGpz8APSuAYH9vTKXcQuZIazeCKWqQbLQ1yUCpZjUdZv3ss8/w/vvv480338Txxx8fVQCxdOnSlC2OSA37VWFWEnNE/6a+TXbk9jZ2JthSYcvBFn7hOdjSlTB0k0loj2lR2DCRVJTrBBo7+1zMKWHW2PtTCQcbu8CLPeaGFucgGCmuSCYXUJsQ39NqVlMFEJH3m++yo90XTHmYNRyWuBPV0zBrKgogRg7K4/Nq69t8qIox4swsRluTMMe0tZ+HWU2LuaKiIsydO7c31kL0EvuFodSt3QGEwnKPJoLoj3g7o8NLifj0GyW8mm3te7R5sKKwYY4Pc8b6vACCO3OpCbOKPeYqC93cnenwB1WFH0boEMKjbd3BHlSzyt83vdYksdbFhGOhm4m51DpzHf4gwpFJD00dfnT6g8h12uHtiu+UxqJHOXNdyuiwikI3DjR34bC3q8diTpKkuP0LRZjI7u85c0mN8yKyCzHMKknyyac/DxwmBjYtkQuWdk5pPNYKuXKN7dkVZm3WOHPdccKs6SyAiIXL5AQIFmId7HHDbrMiJ5Izpy38MAK7uJflu2Qxl6Q71t6tnzMXjDQmznOpL7Xs/RbmOHDI251yMad1oQ42d2F0RYEQZu2jalZeHOJAlScHB5q7cKilG1OGmX4pFZ3+EPyRfZawabCLtSbp31EpGuPQzwmEwjjUok46pVAr0Z9p7pC/33UG83OCoTA+39OPnDlBsHUHhDArokOJ4bDE3a3ewGgBhD8YhsSGhsZB7DEHgOfMAebz5tr9TMzJa0s2b433mXMpAinXaYM9Ev2IFWplnwMTfamezarN1TsQubFh534zTYNTMc6rMMeBwUXyTNbDKahoZd9Zl07LG5HCAdKaxJCYO/HEE9Hc3Gz4RU8//XQcPHgw6UURqeNwSzdCYQkuuxWDPfLBREUQRH+GXTxDYQl1rYmr57YeakWHIASyT8xpCiD80c4cr2bVCJ5fvbYFU36/Al8dau2VtSUKhYlOmhF3SuwxBwA2q4ULwk6T1ZbMUSvNc0X+furGeVkslrh5c6IzJ/6cKrRijqUctPSgz1xHMhMgIvvYk+PAYI8swLXmQjKI+XKJQuu8AKI7aOiGIVsxFGbdtGkTNm/ejJKSEkMvumnTJvh82XVC7K+wEGt1SS6cNisOe7v7vNfc0g0HYLNacPEJQ1LyelsPevH8mr24feaxqCh0p+Q1if5BIBRWhYMONnfxC78ea79pVP18NOvCrPrOnFIAoYRZxRyuLw54IUnA5gMtGFdVmPK1tRh05gDZGUsUJmXnM+bMAbIL5guGzTtzke9JWQFzLc0LqlBY4jcCYs4cIAuYxg6/MWcuxQUQUWHWli50B5TQpJk0m/yI++kPheEPhuG0GwvoSZKkOHNuBzcTUunMJapkBZScuWBYQncgzEPz/Q3DOXPnnHOOYVVLkyAyB3byqynJ5ScQvQqr3mBXfRtu/+dm2K0WnDeuQhUWSZZnP6nF0g0HMao8Hz/61sgUrJLoL2i/20by5li+3NjKAnxd15bRzlxbdwAuu011QY1qTRKjmpVd9MKSfFF22eULGnP16lt75z3zMGtebCfIbrXAapHXJQua+I4Rq1AeXprHH8t12tHcGTAt5ljYkDlzyVSzijcOBRoxF69xsFIA0XfOHPusbVYL8kwImjyXsm2HLwin3ZgQ7AqEEIxUYRTm2AUxlzpnTu97JZLrsPGm2W2+wMAWc7W1taZfeOjQoaafQ6QeVslaXZyDhshFqi9z5l5dfwCAfFfk7QqkRMyxE1UTTbMgNGi/24kqWkNhCesi+XKzJw6OiLnM/F41d/jxrf9ZibGVBXj1x9P54/Fy5lgultgkttuviDl2LNW3mbvALvmkFtUluTjnuIq427V0sLBebAFgsVjgstvQFQgZEjT85rRUcVtZ42CzYVaWQ1VWIIu5YFhCMBSG3WY8lZyJOafNyvcpw0iYVcmZS3UBRICvyx8K42BLFy8MKsqJPY1DD7vNCpfdCl8wjA5/0HDbHva+7VYLchw2XsEqhlk/2XUUexo7cMUpNabW1NQRv+WNiNVqQb5Lrlhu6w6ivMDwn8kqDF1Zhw3rYekJ0SMkSU5SLs13mX6uGGb1R3oymRFzexs7AADDhDthowRCYfx7g5I72dYdxGCP6ZeJgp2EaZoFoUUrbBI5c9sOt6KtO4gClx1njB6EP72zw5QzFwyF8fbWOkwdUdLrIf/1e5vR1h3E53ubEQiF4YiIDlbwwS64scKs+S477FYLgmEJXYEQPHDIkyMiYoT15jPCnqMduO+Nr1CW78T6u8/T3a7dF+SvH2/fuB1WdAVCCZ2x7kCIuzrDSqLFnNm2K2yiQZkgTrqDYeSbEHOx8uUY8RoHs5CuJ0d+XsrDrJGWIKMr8vHloVYcbO7iLXvM9Jhj5Lvs8AX9pvLm2BoKI+KROXNH233wBUNo6w7iuufWoysQwsiyfEw7ptTwa7cY7DHHKIiIuf7cnoSqWTMcSZKw4NUvMOX372LVzqOmn79fCLOyvBl2h5YIXzCEOYs+wcV/+QSBJLp2f7i9QXVhTJX4Yifh/t4EkjBPlDOXQMyxEOtJw4u54Gjq8CMcNpZS8t7X9bjlHxvx+7e2JbFac2w96AUghySPCIUdLMzKnI9uf7QzJ1b9MbEnHj9mxBwrKmlMsJ/qIrlRBW67ajqCFuZoJXLm2Lks32VXXcRZ2KzDbM4cK4AQbpLNhlpjtSVhxHfmImHW3iqAiIjM4wbLeZBH2rp5ZMbMKC9GMu1JWjUza0vynDxH8ojXh6c/ruXfxf9+WWdqPWZy5gBxpFf/vWaQmMtwXly7D//eIIcqtx02X3EmhiVYbyGvQWeuvtUHb1cALZ2BpITYP9fvV/2cqgOJnUD7MvePyA6YsGEXjURhVlb8MHVkKUojLSpCYclwxfehiFg8bKKnXbJ8ecjL/181bzPiUlRGxGisCRAuuw1ujYMlHj8NBqp+GexCKucg6R/TzEVjjoweLoexxsEsX25Yaa4qJMdSN7rMVrNGXKYCt53nIMYTc5v2t2D51sOqx/goL7NiLqAOs/p7KWduRFkeXHYrJAn4+nAbAHM95hjM/TTTnkRsSwJA5c59ddiLv6/Zw7d958s6U5WmRqc/MJSRXv33mpFWMffRRx/hoosuQlVVFSwWC1577TX+u0AggF/84heYMGEC8vLyUFVVhSuvvBKHDh1SvYbP58PNN9+MsrIy5OXlYc6cOThw4IBqm+bmZsybNw8ejwcejwfz5s1DS0tLH7zDnrH1oBe/feMr/rPZu0YmxACgulh05ox9oRsEV81sw8ij7T68/3U9AOUi05qiQcfsJExhVkILEzZjI47EwZYu3YtEOCzhs0i+3NQRJXDYrPwYMRpqZaGkvrjj33JQEXOHhIpAdmFjfbxihVljOXMqMdfuM3wxbRT2Tbzh5YdbmJiL3+2fT4FIkDe2t0kRcyJKzpzZatbI5AaXHW4u5vTX8NMXN+DHL2xQjUfko7xc0QJJT8yFwhL8oV5uTdKthFSHRCp/2c1AsmFWwKSYExoGM9h34cHl29HpD2FsZQFynTYc8narvt+J4M6c0TDrABjplVYx19HRgUmTJmHRokVRv+vs7MSGDRtwzz33YMOGDVi6dCl27NiBOXPmqLa77bbbsGzZMrz88stYtWoV2tvbMXv2bIRCyoF9xRVXYNOmTVi+fDmWL1+OTZs2Yd68eb3+/npCa3cAN764Af5QGGzyltnu7eykU5bvRJ7LDk9kuLLR1iQNbeJJ29xB8NrGgwiGJUyqLsKkajlRLlXii52EUyUOif4Du3COG1wAi0W+SOoVNOyob0NLZwC5ThvGD5G/o6WRi8NRg2FHJmbiiZpUUN/WjSNCxSlz5roDIS5AqiIXSrWYi4RZHYqYYzeF4nkgEJKiJkno0Sg0GWZiNhZGnTk+nzWhMyfn79aUqPN3kxFzwVBYNR/Vrdk3WsJhibfUEEP3sUZ5MfTEnOjC9VprEiFfbUgRE3NyZCeZ6T9JhVm7lB5zDHbD8c1R+bP82XnH4sxjBwEwF2pt5gUQxoQpE6P9OWcu6dLC9evXY9u2bbBYLBg7dixOOukk069xwQUX4IILLoj5O4/HgxUrVqgee+yxx3DKKadg3759qKmpgdfrxTPPPIPnn38e5557LgDghRdeQHV1Nd59912cf/752LZtG5YvX45PP/0UU6dOBQA89dRTmDZtGrZv344xY8aYXndfcNe/t2BfUyeGFufgjNGD8I/P9pnug6RtsGnamWtLzpmTJAmvrJNDrN+bMhSb97cASM1dkXgSjnchIQYmzKUaVOBGeYELR1p9ONjShUEF0cVDb2yWXf6ThpfwYoKyfBd2N3SoXOl4tHIx17vfxS81TX0Pajr626wWlBeqW2yEwhICkaInd4IwKyALRiNhK3FiRLwbqrpWeY0pc+Ya9Zw5+TJmpppVTOTPE8ScnqgSZ52Krm28nDm91iSiYBRbk5idLRsPMV+N9eRjeZF95szxMKuyb6qE78LYygKcd1wFuvwh/GdrHf775RHcef5YQ6/dlKB/oRbKmYvBgQMHcMYZZ+CUU07BrbfeiltuuQWnnHIKTj/9dOzfvz/xC/QAr9cLi8WCoqIiAMDnn3+OQCCAmTNn8m2qqqowfvx4rF69GgCwZs0aeDweLuQA4NRTT4XH4+HbxMLn86G1tVX1r6/wdgbw1hY5N2PRFSfyMKXZ2YFijzlAEXNGc+ZinbSMsPmAFzvr2+GyW3HRpCp+IKXCSROTnFsjQ6wJgsHETZHgSMTKm2vu8GPJJ3sAAD+cWsMfZ20qjLYnYTcU7f6g4aKJZPgyEoJyRkTnwUgIkw+yz3Fw5405VKIwkZ05+bm8AEIr5gz2mlM7c/rH9KEWgzlzBgsg9umEWXOScObYnE6n3Qqn3Qq3I36YVRRkR2NELApiFHjw861mH7H3abdakOtSZssy4Z0KuJgTjgNGMmKO9ZozU2QSM8xapHwXbj57NKxWC84aWw671YJd9e3Y3dCe8HUlSTJfzcpHevXfaI5pMXfttdciEAhg27ZtaGpqQlNTE7Zt2wZJkjB//vzeWCMAoLu7G7/85S9xxRVXoLBQzoepq6uD0+lEcXGxatuKigrU1dXxbcrLy6Ner7y8nG8Ti4ULF/IcO4/Hg+rq6hS+m/iwE3S+y44TqouUE43JnJAoMcfCrAZFkCrMaiJxlCWVnzWmHJ4cB78zi3VXtL+p09QdtegQ+kPhlOeaENkNHySe68CQiCN9sKUzarunPv4GHf4Qxg0uxHnjlF5pgyKVjY0mnTlJUmZ99gYsn+i0UXL7hoMR171FGJeVo3HeRKfLaYvOmdNW/hqtaG1qF505/fdcx8KsRYnEXOICiGAozCMN2jZJeUm0JmHOHBNhicKsYhRAFLNKaxL9nDmt4GV/w+2wqSdgpDDUysOsbiVnjmFmlBcjmTCrV1MAAQBjK+Vr95iKAlwwvhKAvJ9YWxIjodZ2X5ALX8POXBLrzzZMi7mPP/4YTzzxhCo8OWbMGDz22GP4+OOPU7o4RiAQwOWXX45wOIzHH3884fZauzqWdZ3I0r7rrrvg9Xr5v952HUX4QeDWnGhMO3OyI6F15kJhydCXWhVmNeHMsdwbdhJnd2bak9r+pk6c+T8rccPznxt+ba3NT0UQhIgobvScuaYOP55bvQcAcNu5o1XnAZ4zZ1DMiTcovRnC2XpQjgycf7x8ATzU0q1yKIpynVE5caIDZLdZudjTmwRzxGBFa5NRZ87LwqwGc+bihFkPe7sRCElw2qw8UsHI4WFW4+dHlnfLRIoiKGOvQYwqiN8NVgCRKGdOvHlm53GX3cqd1nh/2yzhsMRFZmGOHUOK1E5mMmIuuTArE5TKvpkyrBgvXjcVL1w3FVarctyx7/V/vzyS8HVZvlyOw2Z4moPizJGY49TU1CAQiD6Ag8EghgxJzexNkUAggEsvvRS1tbVYsWIFd+UAoLKyEn6/H83Nzarn1NfXo6Kigm9z5Ej0F6ShoYFvEwuXy4XCwkLVv75CtMgBCCdpcwf7fqFhMKC+EzTS1iPWScsIXt5pXL4w6lUS7TjShrAEfF3XZvi1tQcjFUEQIrzLfa6TOxLaXnNPR1y546vUrhyQRJhV+P71VnJ1c4efv4dzI+vtCoTQ3BngN07FojMXUIdZ2THPRJM2Z84euag2GHTmGjuE3pE6x1+7L8iP1UqjOXNxblZZvtzQkhzYrOqb8GQmQLCK+HyDzpyq8rdNdOYS95kLhCTNVI4w/5vyBIz4QtIsYn5fLGcuuTBrD/rMaf7eaaPKonJYZ46rgMUCbN7fwh1dPZR8OePvI5/lzJEzp/Dggw/i5ptvxvr16/ndxvr163HrrbfiT3/6U0oXx4Tczp078e6776K0VN0hesqUKXA4HKpCicOHD2Pr1q2YPl0edzNt2jR4vV589tlnfJu1a9fC6/XybTINsRIJUPowmWlNEgpLPCxRXaIczLwIwkDeXEOSOXNiqAtQ3oc2X4FdiFo6/YZz37R3htRrjhDh46NyHBgaceYOCM6c2pU7NsqdL8tnYs5oaxLl+9db+Tis+GFYaS7K8l0oj1wIDzZ3qcSrNozK25JEHo8Kswq9yABjI73CYXXVq14RktGGwfL6EhdA7G2Sqx+Hx5hEk0w1Kzuf5XNnLv4Ns/g5i2JW+zradTGhLJ6n2HmcvW+lACQ1YVZ20yznAtpQUeBSCWDW1cAMeSlqTaJHeaEbk6uLAADvfBU/1MpSIIy2JQEGRs6c6WrWq6++Gp2dnZg6dSrsdvnpwWAQdrsd1157La699lq+bVNTU9zXam9vx65du/jPtbW12LRpE0pKSlBVVYXvfve72LBhA958802EQiGe41ZSUgKn0wmPx4P58+fjjjvuQGlpKUpKSrBgwQJMmDCBV7ced9xxmDVrFn70ox/hySefBABcf/31mD17dsZWsiphVvkgcGtOwkZo6fTzvAIxLFGU48SRSDPgeEiSlHRrEjHUBQjOnOZvNkfCNYGQhA5/KOFJH4i+M6SKVoIREMZTFeU6+ZBv0ZljuXLjhxTi3OOic2nLIo2DGw04c5Ikqdzm3grhbI30BxtfJbdPqSrKQX2bXKUrFnzo5cyxHmp6feaOrSjAzvp2QwUQ3q4AQkKhh54zZ7QtCWCsAGJfozr/VyQ3iTBrhyY86k5ww6wqgGiPziWOJVgsFgs8OQ40dvjh7Qrwql6xkTMQEdvdQdX793YGsPWQF9NGlqrCkUZo1Vw/7JHQNDsOkguzsqbByYzzMiYzzh5bjg37WrB+TzOunDZcdzs+orI4+rugRwG1JonmkUceSdkfX79+Pc466yz+8+233w4AuOqqq3Dffffh9ddfBwCccMIJquetXLkSM2bMAAA8/PDDsNvtuPTSS9HV1YVzzjkHS5Ysgc2mxNJffPFF3HLLLbzqdc6cOTF722UKYsNHwFhOiRZ2113gtqsGR3sMOnPtvqDq75mx19lrM0euUKcsvFnoc9Xc4U9OzPXjOy3CHOLNQqHbzh2Ptu4gWrsDkCTg+TV7AQC3nhPtygGKM8ea6MbLq+0KhAwJm57Cih+OHyKnegwpzsGm/S042NLFb4iK8wRnLiJqeG4Wc+ZYzpxfXc06qjwfgLECCDH5X3wNLUYbBgOKMxUv8rCHz4iOJebMF0Aw0Z+nDbPqhHpF0S4K/XitSQAoYq4z2plza505Qcz95o0vsXTjQSy+5mScNSb6piMeSuGBsqYhxTlczCUVZnUmH2Y1+veOrSgAAHxzNH5F655Ij7phZSbE3ABoTWJazF111VUp++MzZsyIG14zEnpzu9147LHH8Nhjj+luU1JSghdeeCGpNaYD7cHITtJmbHitO8Zgo1wSzWfV5gyZOQh4mDXyt3hVl06Yla23uiTxa2vvrKgAgmBob2DskYkOLZ0BHGrpwnvb6tHuC2JMRUFMVw5QxJw/KLt88UJEWle4ty4UrC3JhEhjY1bYcaili79nj9CapCsQgiRJ3JmLypnTVLOOrmBirjuhgG3SirlUOHOOxM6cXo85QJzNaiJnThMeTdSaRDzPdPpD6PQHkeu0xx3nBcTuNdcdZI5pxJmLEWZl0y7EaRNG0TpzADC0KAefQa78ddhMZ1eZLoCQJCnmOuJxTOSmYnd9B8JhSdeR3BP5LsQKueuhjPPqv2LO0Kcq9ljT9l5LVy+2/oz2IHBr+kOZeQ3tXZHRnDltMrSZmXZKzpy6AKI7EFZ1Pxc70DcbnEpBOXOEHi28GEDJpWHC55uGDiyO9JW74cyRuoIlx2njrS4STYHQCpneEHOt3QF+8Tq+Si3mDjZ38WOoONfJRU1YYm171AUQitiTj0F27Iwulx2R7kA4YYJ4U4d6n+jmzBlsGCyuT68AQpIkocecfs6cqdYkfrWjxoSV3g2z9qbxaKQIgo/z0hEssaZA+KJy5qLFbFucZtR13m7sOKJfNMZcRLHwgBVBJDP9ATCfM9fhDylFGAaduZqSXNitFnQFQqiLU1nNJoGYEXMFgpgL9WI/yHRiSMwVFxejvl6es1lUVITi4uKof+xxouewgzE6zGrCmdNUlDLiDX8WiRZzxg7iQCjMt2XOnBg+FRNQm5MQc9Fh1v57p0WYQ2nTIVzEIsJn0fu7cLTdhyFFObhoUlXc12EVrdqQohbtBT7Z5OrmDj9u/sdGfLyzIep3X0WKH4YU5fAGqVVFSpWuWM3KzhMA0O0PR+VmiTl1vmCI3xxWFrp5TlGivDnm2DMHU8+ZM9owWF5f/GrOo+1+dPpDsFjApxmI5CWRM8dEEnuumWpWADja4YMvGOI3p3opIrHOt1HOnCP6/ceb+XvtknWY/edVqI2EG7UoZoAQZo18Z4wKKy1mq1nZGpw2q6qXXjwcNit3XvWaBwdCYV7QNNxEmFX8fMw4uNmEoTDr+++/j5ISOQa2cuXKXl0QEd1sMZnWJMyl8GjDrLnG5rOyJN/yAhfq23yGXQfxpMXWb7dZkee0ocMfQmt3EKWRCwHrFySuNxHsZOK0WeEPhWOGWUNhKap9AdH/4d/5GI7EV4dlUXTdGSMShpnK8l3Y29iZ0JnTHhPJOnMvr9uPNzYfQkunH2eMHqT63VaWL1eltEYSw6zMYfTkOuCwWeGwWRAISegMBFVzWQGo+tCx49RikV2LQYUutDUEUd/WzXPoYsHCrCPKcnG03aeb5mC0YbC8PuaKxT6/MSemypPDhamI2JIlXnhORL8AInGfOUB2bcWUj0RirjWOM8d6zYnOZLyZv7sa2uEPhfHWF4dw09mjddcqCrepI0uR47Bh+jGlUdsbgb2/Tn/I0NgxMVXIzIiyYwblY3dDB3bVt0cdC4DsRgfDElx2KyoKEn+3GG6HjV8z2rrjp09kK4bE3Jlnnhnz/4neQXtnxVuTBI0dSIByMGnDrOxno2HWEWV5qNecuIz83QK3XSWoCnMc6PCHdJ05bS6OvEY/PDkO1ftlYq7S48a+ps6oO+ZfLduC5Vvr8N+ffYu7B0TqCIbC+GB7A04cVmx4lE5v4Q+G4RTu+ls04X0AqlFGxbkOXHZy4kkuRhsHR4dZk3PmPtl1FEDsY4C5L2MqC/hjTKA2dvj5McZCy26HDYFQUHbf4uTMiakcVqsF5QUufNPQkbDXHFvj8NI8rNvTjDZfMKaAMtowGFCqbfWKD+LlywFKmBWQ31ueiUIqVqWZqACCnWdyHDZ0BUJo7PAL7p5N9+YxZpg1yplTi9lgKMzHZmmdsO6A4gb+Z2tdbDEnTH9gjCjLw+Z7Z6qOFzOwcV7BsARfMKxygWNhNl+OMao8H+98dUTXmdsjhFjNVvnmu+1o6vBHjtPE4f9sw/Qnu3jxYrz66qtRj7/66qt47rnnUrKogY5eNaskGW8sKbYsEOE5cwbDrCMHyXkJRpst6hVeKFMg5NeRu9eLzpz6QvbB9nqc8NsVePyD3arH2R11VeSOX3tBfXfbETR2+KMGkxOpYcVXR3Dd39fj/re3pXUdm/e3YPx9/8Wf39vJH/MKc0oZYljuymnDeRuLeLAwa0OC9iRaVyqZ5OruQAif7ZFbOMW6wWKPlQrCuVDo3cbyf5iYyxVcKqWfWXSYVXuzVxFpX5QozMpCz8MjvekkKTpsZaZhsLg+XWdOZyYrwy24dUZDrYqYY708jY3zYufDo0K0Qi9fDtAJsyaoZo03VUT8+ctDrbxli2qtwvQHkWSFHKCEowFjeXMs/aXAZFj3mEFKEUQsEgn7ePC8uX6ammP6033ggQdQVlYW9Xh5eTnuv//+lCxqoNOqCbOKJ6t4jTVjvUZUAUQkh86bwJljrgRrKOoPhg3NDtROf2BomzZ2+kPwh5T30qxZz4Z9LQCAjZH/MtjJjOUNiQnYobDEc3r66wGbblh7g131iQdii2yva8OLa/embBj9uj1N8AfDqgajYv4YY0SZfHFwO6y4avpwQ69dZnA+a6umIjKZ/M31e5q50xIrj7VZGNfFsFgsKsdRHBQvhlKVnLno32lvulgj4kSNg1kBRFWRm4sD7fs20zBYXJ/e+WVfxI2pKYmd8G61WkwXQbDzA3OcuDuoc35ln83IiNg42u7jPeb0KlmB+M4cE5Da968Wc/Hd31izTJN1xeIh7mMjveZi5e0ZgVe0JnLmyowXPzDYd7G/ticxLeb27t2LESNGRD0+bNgw7Nu3LyWLGsjIJd3qAgiHzcJtfKPzWZWQk54zF991YNMfxOoxIwJJL7zLhCm7a9QWPGh/bohcVLSPMxeAXcxEZ66xw8edCjOdygnjMEfD6BxPxl1Lv8Cvlm3Fmm8aE24bK9yohTnHu+rbuUDkbToE4TOmsgD3XzIBz151suGw8KB8c2FW9l1M5iLx8S6l6EEeIK4WE3pOd5WQi1acq6QiKCO7YhRACGFW7XFaHsk/OpLImYvcLJXkuXRnLptpSyKvL34BxB4DbgyfAhEw9hm08yrUxAUQ/mCYF4uMjIiIo0KYVa/HHKDkLDfH6jNnj13NKp7ToscXqn/+z9bDUX9Tb4xWTzFTBGG2xxyDOZ/1bb6YxTWpcOb660gv02KuvLwcX3zxRdTjmzdvjhq3RZinOxDmjhU7GC0WCz/wjd55srCldnSL2Zy5ykI3P1EaOYj1Ci+UKRDya4jFD7HWw/5+s+bCzu4KWcsD8Y5XDBFlywH7+d5mXP/39THDJZkIE8n1bT7DJf6SJPH5u2zEnB7PrqrFib9bgbe3RF+kRNj3ozugVLe1xAizAsAVU2swfVR0NEGP0nxj81nZd5kJq2Ry5lbtPKr6WevOtcRw5gCo5m2KrVhyhDmlUa1JnEqLI22RVXmhMWeOhVlL85w8jBcl5kw0DAb0hZQkSVi64QAvXok1/YGRY8I1kreL3TQ4lqAURQWLVIgFEPHcRxYeF29QtOFv7TizeGKOfceYe7xhX/Qs01gD7lNBnol+ftqRlEYpdDu4S/xNQ3SolTUMHmGiLQlDaRzcP9tZmRZzl19+OW655RasXLkSoVAIoVAI77//Pm699VZcfvnlvbHGAQU7kG1WCz94AP0E3eVbD+P3b34VdWHVc8jYHb4vGNbND5EkibsSgwpcpuxpvVy9Qs2BlMiZY93omzSPK2HWSM5cV4A3lxYvRNnizP3js31456sjeHPLoXQvxRBsQHkoLCUMQzLqWrt5LlMigbRpfwsApYpTD3Fu8M56WShqZwIni9H5rNyZK07OmWsScjtZRaM2d5S5jdqh4kOKFGEjHuOi+8YLIByaAgh/KOo4HcTDrPrvWZIkfnNVmu/kF0etW5QKZ661O4BbX96E2/+5Gf5gGN86dhDGDS7UewnkOuz8vRmhzacWYvHGeTGxWuCyc9F7tN3Hz2XxwpnsxkAUc9rwN/uvP8Qmcyj7U2984ciyPJxYUwQgepZpJjlzyYR6lbw5dag1GApjf+RmcFgSYdYLxlfiprNGxf0eZTOmxdzvf/97TJ06Feeccw5ycnKQk5ODmTNn4uyzz6acuRQg5hqIVZx6I70e+M/XeHpVLTbua1Y9rndhy3cpVaZ67py3K8DnusonbeMHsd7f5XfxkRM/E29sbqyeM6edBclEGgtthSXwyi8xRJQtnb47I3e5nSZmHqYTUSTHa+wpIt5hJ6qWZN+LRGJcfJ0dR+STvhKS7FmVrdH5rNr8zXZf0NDUGgarYh1bWcBbeIjHgS8Y4iJYm4OqDrMqv2MuujpnTh1m9QXDgnOvDrM2xAmztnYF+bzbkjwnd360zpyZhsHi+th6j7b7cOGfP8brmw/BZrXgjvOOxeKrT45bvZjrUhzJRPiD4aj+cPHmw4ouJs+n7PBHhWpjwUL78njEyJg1XgDBcubUBSCic6T9TrHfFbjtmDW+EgDwny0aMdcLOXNAdOPgZ1fV4pf//iKmQ6/kfZt3B48pl4WaNm/usLcbgZAEp92KwYXG25Iw5p44FAvOH4PJNf2zH65pMed0OvHKK6/g66+/xosvvoilS5di9+7dePbZZ+F0prdVQX9AG/5g6N05MnEkChmxUlQrqiwWS8KRXuxC6clxwGW3IT9yUuhJzpxyFy//nq2PhS3afUF+gg2HJb4GSVJeMxhScldK813czWAnDjGPK1vEHBPnZhqephNxv2rDO3p8I5yUE7ld7HvRnkDcimKOOXPNMZoGJwOrZhUvwLFg3zt2YxEKS6Y+RxZiPX1UmXJMCmKOFSlZLdGCQazSLc5T3q/ovmmrJnMEp5+dLzyaMGubL6jrbjVGih/yXXa47LaoPFiGmYbBgBhmlP/uvz8/gP1NXajyuPHqj6fh5nNGJ+wbKVbxJkK8Ucgz4swJExWYmGvpDKApkioSL8xa6LbDYZPXzkLUes6ckjOnLuoSv1NtwlouGD8YALC2tpG75OGwxF3HZIRUPMSRXgeaO/H7t77Cy+v244sDLVHbensgKLkzpxFze3ghTK7ptiQDgaRrlYcPH46JEydi1qxZGDZsWCrXNKDRs6e1cxUZ7MQkXiQ7/SF+Bx0rAdWTYKQXu1Cy0AvrDt9mYKSXkrekFvba1iS8X1VZLpgBycRlS1eAr1/ctkM4qeW77Pxk5eViTnDmsqRiiV04uwwmbqcbUcwZLYLYbcKZ4591HDEeCIVV4fedR9oRjDQDBaJD/GYpcNl5pWa89bJjtaLQzcWG0VCrJElYFXHmThtdxos2xJZBzUITZO3Fq0qoZhXzYpn71hnDmROr4g9HPjt2fihw2bmg0cubY58Nc5u0xzTDTMNgQEwhkdf7yW65SOba00fgRIMuSo6DCY3EYo59h90OK28gHa8AQrlBtaMox8E/631N8vc6XmsSi8XC91dTxOmNcub4BAgWZtVvedMqOHPVJbk4vqoQYUluGQQA7f4gmJHXW85cuy+EFz7dx8d1xToP9CTUq4g5dc4cy5czM8ZrIGFazHV2dmL+/PnIzc3F8ccfzytYb7nlFjzwwAMpX+BAQ8/Z4iES4WQTDCkVa+JFh10QnDYrf55ILBdAhOUjsXATuyMzIpCUisLYBRBt3JlTLgzaogztxYRty6c/2K1w2q1R1XQNWZgzx3IgzcyVTCdJhVmFsUOJnTkm3PU/v8Z2P8Ro5q76dpUIMltBp8VisaDMQONgseqcHyMGZxjvaezEwZYuOG1WTB1RIhyTikgV565qKS9ww84bBgs5cyzMGhnZBSjOj9Vq4f9/JCK4mItpsVh4qFUvb65RK+Z46oTGmTPRMFhcH2t/tK5W7rt3momilVyn8TBruyZfDlCEbqzWJNoGy+z9s4bO8VqTAEBpXiTPrkMp2gFEZ04d4tUrehB/x859542rAACsjghgPkbLbk3Y2NcsrMFyY7sPL69TOlfEcuh7UoTB2pPsOdqhqu5mVc3Dk6hkHQiYFnN33XUXNm/ejA8++ABut3KwnnvuuXjllVdSuriBiFIFpD4IYuXMdQrCTrzosPBMoWZ6AoPlFOmN4lGcOfnzNVPS7dUrgOAhGZYzxxK7nfxixZKrtW6I1q1hJ2Htax7JwmrWbAuzqsSc11gBhDrMqp+H5guGdDvfi7DvR1m+E067FV2BEC8kKHDZYU8wrssIfD5rnPWK+Uu8WtugM7cqMof1xGFFyHXauajyxnLmYoSNbVYLd75U1axiAYQmnAcoYo/dMInOCe81p5M3x/ZFaZQzp87xMtMwWLu+z2qb0BUIoTTPiTEVBXGepYb1izNyU6Q9jwBCmDUyZUdEe4PNQq37Is2M4+XMAXLeMaA4c3wChLbPXIxqVkAt7ngxRuRvMudycyTUqW1rlUpY4+B/bzigMgLqYnxfkm1NAgCDC93IcdgQDEt8HwPKWLdkih8GAqbPeq+99hoWLVqE008/XSUUxo0bh927d8d5JmEEvWa/7GQjhlnFpHm1Mxc/dyhhzhyrZI2ctPJNdM726jhz2mTpZsF1KNL0YtJeTNi2fHxO5MTNLkTeGDlzWePMBViYNTvEnJjLZiTM2h0I8UbDgPxZ+XX6iYkXiHifX0O7/HcrPW7e94u5OUV5qbmIJapoFQsMCnMcQtsDY9+7T3bJTsrpEfeJHZNiVTdrwB3LmQOA4yrlqjyxgao4p1SpZlUcGib2WIhMPM/wKRC6YVZ5XyjOXHTOnNmGwQBU81bf/7oeADDtmFJTeVEszNpp4DhiN3ri2C+2jyQJqmbmQHTIkEUsWJFYQYL3ycQvyzn0RRVAaJsGxxFzmqkTk4YWAZD7rzV3+IU0ndTmywHK/mI3zcx5jRlm1cn9NoLValGKIISKVnLm4mNazDU0NKC8vDzq8Y6ODlMDdYnY6CWOxho3I4aiYjlzerlDiXLmjrbJF5GyAvkkVGCwNYkkSUqz4qgJEOrWJOxvF+c5+MWKhZUa2rViTt5WuaOWX0sUiPL0B+V5PRFzSzccwO2vbOqTIgou5rLRmTMg5mqPdkCS5Is7Cwuyi5oWUcjEy33iznG+C6Mj7g0biaX93iVLWYLGwexYsFiAfKc9Ko0gESwUOTYiyHjOXGe0M6d3HP/xOxPxyvWn4uThSl5ZjqppsDrMKv6eIVb+JmpPwnvMRYRuoaZ3JGC+LQkgN0Vnum1lRMyZCbECQpjVwDEbqz8cu1kGokOt2hts7czneDlzgLK/GjXOHA+zanrcaXMQY02EYJEbT66D39BsPtDSIxGVCO3+uuFbIwFEizlVEUaSeXvavLlQWOK9ONOZM7f4k1pMvO+/2BxpoZRJmBZzJ598Mt566y3+MxNwTz31FKZNm5a6lQ1Q9BJHc2KFWX1imFW8o49vcbMLnt58Vl1nLsGJst0X5GXqeq1J2GBuFjotEsOses5chzpnLl/jzLV2B9DY7oNYId8TIfbH5V9j6caDWPJJbdKvYRQm5rIhzBoKSyoH8YiBalbWluSYQfmK29UW2xEWe3EZCbMOKnDh2Eh+DTu59rSSlcFETqwRW4CQiO6yw2q1cGFj1JljnzdrqVEcI8yaqNVKcZ4TU0eWqm6iFWcuGJWbBSAqj0o8R/DGwTph1qYOTZg1hjN3MNLAWSzQSITFYuHuHHNfTjvGpJjjrUmSC7M6bVZeiOXTdgzQhC6Z0GckzJljrW46dAogoqpZ5f3Jbn7EPEytMwcAk6qLAACb93uVytsUFz8AaifzkslDcGylfCOlvakTizAShaD10Fa0HvZ2wR8Kw2GzmPpupZr/bKlDa3cQi/vg2mAW03t64cKFmDVrFr766isEg0E8+uij+PLLL7FmzRp8+OGHvbHGAYVe5+xYYVbRmWto80GSJFgsFt0iBAbPzzFYzcqcsEQXKm+c5Ft2cmGDucXk7mLuFKqduZI8J5o6/Pwiok1cZifX1q4gt/4tFvlvJCvmGtt9/LUWf7IH808fqWrpkGq6sijMqi1KaPMF0eELqk7yWli+3MhBeQiGw6hr7Y6EST1R22rDrOz7rEX8fo6ukE/67EKYqlyhRG60dsi62e7yzEHKjeQhFcVwy1uSaLUitiZhzpx4LIrfZW1jcnbzpnXGGbrVrMJ75qO34kxriIXLYeXHwNDiHNSYDKXlClW8ieDnEUFoyFN2bOgKhKKcOaVdlLx9aZQzZzDM2s4KINSOqd5s1sFFbuxv6ortzAl/c9JQD5ZtPIjNB1r4WnrDmWPpLQBw1fThvDWU9qbuaOT4zHXaki7C0Io5NsaruiQ3YZua3oQ1LV7x1RF0+UO9em0wi2lnbvr06Vi9ejU6OztxzDHH4J133kFFRQXWrFmDKVOm9MYaBxRKmFVTAMEbSwo5c8LF1R8K87syr06ok5FoPqs4/QGA0DQ4/oVKb/oDIJ+w2MHf2O7nie7FuQ4U5zFnTl5PfeROjyVAM8dOGY4dOWG5FTeDWf2s51d3IIxgKHZuVjy2HW7j/9/Y4cern+83/RpGkSSJXzhSFWbtDoQw75m1ePyDXSl5PRHmaDhsFi6oE4VaWSWrWWcuGJZ0Z3WKzvFoTZJ8qpy5ggROmzacZXaINxMdTEx5uFsuVrOyIiHj70nstxazAEK4uBZpCqTYe9YLVfK5rBGnycPHeSnbs3YdNSZDYeIazbpygCKKDYVZY+TMAeoiCBFtMn9UmDVhzpzSaBiIVQChbhrM/l6Vh82fjlUAEcuZa9G9fqSCMREn7uyx5RhbWchzLDv8IdXNMytaiDd+LREsZ+7Lg6147L2dfCJMOkOsvmCIn+86/CGs3F6ftrXEwpSYCwQCuOaaa5Cbm4vnnnsOW7duxVdffYUXXngBEyZM6K01Dij0qoBynNE5c9rGqkyEtXTGfg0GezxWpZ44pslsmDXeOCWLxcLvbPdGDnarRRZk2gIIdrFmJw/tVADl7lNpjcDyfEZG7ujk7c0LpK8OyycNdtF78sNvooafpwpRrBhpqWCEDfua8fHOo3h21Z6UvJ6IKKYrIiG5RKFW5swdMyiPXwT1nB/tKCu975tYbT2sJJffJAD6xQJm0Ta51qJNNE8k/rSwcCA7rmM5c+x77zHxnpRq1nBUnzlA7dJFF1nFD1WyXEdtNWtbtzJSb2+Szpy4rumjzM/4VmbSJj7meYqHzvvX9prT5jGX5sfOB9ajRJgoEgyFeQ9NJh6VPnNhSJLEv0PsxpQdd5Ik8WNCFGvHDS6Ew2ZBY4efz7HtDWdubGUhPrrzLDz+gxMByOcBJmTF9iT7mxQXLVlGDcrHpKEe+ENh/O+KHVj4n68BpFfMHWrpVrVEen1TZo1gNCXmHA4Hli1b1ltrIRBvAkSsnDn1hUMZgRU/PDOyTLGwtSeupg4/wpIcrmThFKMFEIozF/viw056+yIl5kW5TlitlugCiEiY81jmzGnCrKxEXmyNwJy5ocU5vOGrkSbHWpgzd+3pw1Ga58TBli689UX8oe/JIu77VIVZ2ZDzpg4fwjHG7PQEcf9XRhLc4zlzkiTxnLmRojOnI+aaNWF/vSIWMcxqt1kxcpBygk9ZmDWhM6dOhzBTzRoKS7yil4dZc5TnM0eZnQvMOHM8t9Yf4i6+yxHdmkRcO4OtJdZ3UZKk6DBr5PlspJ4kKUnqw0yGSUVnbnoSzhxvTWLgODoQyesbWqxeo1tTiMDQFkAMEpw5u9WiKp6IRRl35nyq12YiWwyzdvpDPO+Y5Yax0GqHP8TzgsXPzu2w4bjIvNE1kX5zvZEzBwA1pbkq4V0Ro6I1Fc6c3WbF0htPw6OXn6A6vkcMSp+YOxAJsTI3/f3t9YbTKvoC02HWSy65BK+99lovLIUA9OfqsQNenTMX25lLVABRXZKDikIXAiEJG/e1qH7HLpSleU7er8toaxI9Icpgd5PsYGdiU3TmuvwhXgk1pjI/8rgmZ86tyZnrDvJ2ChUFbmHsjHmBtC1yZzu5uhjXnj4CAPDEB7tNzdw0iijMAyEpJQ7g4UiVZFjSL3BJFrY/8112HmKJJ+Ya2n1o8wVhtcgXd6VCNHaYtbnDrDMnXyRHlStubE/nsjIS5cCJPebE/xo5uYsubC4PsyrHjDdG+x6juJlDFQjqhFmV/9fe7MVrvNvmCyrzmiPiREydaO0KoLkzwI9ds64MEzZjKgr452oG3prEgDPHLsriSDR5DdEjvSRJUo3zAtRh1nzNDO1YMCevOxBWpRLEahosFj+wghR2HLBrg8NmUX2mgNKihDt3KR7lpQebrS06c6kQc4Cc03nxCUOw4mdn4tHLT8DV04djzsSqHr1mT9jfJJ9bp44sxTGD8uAPhvnkjUzA9Cc+atQo/O53v8Pq1asxZcoU5OWplfItt9ySssUNNMSSbr0QiHii0TpzLPG0JU6zUUAOeZ48vARvfnEY6/Y0YdoxSljjKJ/+oJyw+IUtQZg1UX879josFMMuUmLTYHahdjus/ILQEmk9ots0uCvAK/DKC13Id9nR1OE33I2f4QuGsCvS1+i4qkKcPKIET3ywG9uPtGHl9nqcPbbC1OslQusidAVCfLxQshwSTqqN7T7uoqQCJdfIxk/i8cKsu+tlV25ocS5cdhu/SB/VaX3RrAmzxhLjHb4gv4lhryc7uLJ72tNRXoyEzpxm7J6ZMCvLj7RYlAu63WZFgduOtu4gWroCKM13JUyXiEWOqgAiOsyaEyfMGi9UyRre5jptfDuWOnG0Xe5vxoT9YI/bdOI72w/JhFjZuoDEOXOSJAnOnFrMxYp+dAhOGdtf4jFlpFoz12mDy26FLxjmPReddivvoyc2DRZnr2q/U+L0B62APKG6CM9/upf/3FvOnBYmOI+0iWJOfo89FXMMJuouPmFISl4vWVjxQ3VxDiYO9eCRd3fijc2HMPfEoWldF8P0lePpp59GUVERPv/8c/ztb3/Dww8/zP898sgjvbDEgUObT7+km4dPguoTjUiDJmcu3oXtlBElAIB1kf5c/DU0rgegiCc2bkePRP3t2N0iu3PTirmWrgB32AYVuPjjkiS7FdpqVrHPHDuZVBS6VDMEzbDzSDuCYQmeHAeqPG54chz4/inVAIB/bzho6rWMoA1xd6egCOKw0KA33rSFZOBi2u0wFGb95qhSyQoYqJY0EGZlNxs5DhsPd4wWnLniFDUNLkwQNtWGWfn2Bm4geFsSh011URbz5kQxVmxCkDNRI7Y4EcOAbqe+mMuNkZfL0I7yYojzWVmINZmLOPuOzDq+0vRzAUHMRdZe39aNqxd/hjc2q/OamjsDfP9rW1zEcub4eCyblf/eabfyfccq/eNhsVj4zfGhyPEpOmvibFZxwoPSRUB+rFXjBouwIghGb+TMxUJ7UydJUkpy5jIR9r6GFufiokmyQ/jxzqNREYV0YdqZq63NvP4q/QV2ILtitPZwC7kwDBYOsVstCIYlXiWoN0VC5KRhspjbsLcZwVCYh1S3H5FzxthBCqj7MbV3B+HKj33XHa8AAlBO/IqYU4dZQ2Elx2pQvgsOwa1o7vRHVaGxE1abL8hzxcoL3LwPnZGJFSIsxHrc4AJ+kZ06ohRPfVzLR8mkEu1FMxW95g4LTllTik8yrDVJvssmhFn1R3qJPeYAZUSWXs4cy5l02CwIhKSYYVbxZoN9RmJFqydFTYPZBbMrEEIgFI5yTHtSAKEUP6hPv0U5TuxHF7xdfu5y2zXtQxLBbvpYSBTQd+a0N125kVAlC/mL71nbY45RILjje5PMlwOA3397Am6cMUo1zcIMuS51mPXZVXvwwfYGNLT5+IUXUEKs5QUu/XOscFyKqSOi8C7Nd8LbFTDcR600X86/ZWJO/Nvs8wlLyn4udAvOnI85c9GVrIyRZXkocNmFZr19FGbV3NTJERF5DVrnM9thjm51SQ6OGZSP46sK8eWhVvxnax2umFqT5tUl4cyJSJLUK7lEA5V48+xilc2zMBQ7aBrafQiGwvyAjpc/NKayAAVuOzr8IZ70Hw5LeDNyJ3vOccqUD5vVwu9841W0JgoL8dYHrC1J5MLgdtj4RYaJSTb0mzkBzcJJIl/TmgRQnAMWZgXMT4H4iou5Qv4Y63e1t7Ez5d91bZg1FWLukODM6U1aSBZVAYSBMKvYYw5QQvctnYGY+YHsDpdV8cX6/GI5x8NKc1HgtsNhU/KMeorYgyzWTYF24LmZAoiuAOsxp53GoDhzzR1Kw2Azk3Xcmte0WGRxzBDFXFRjcuG52u9io9D7UYS7490B7I20JRmWRMWh025NWsgBSp85Fj1YuuEAADaBRDlu9UKsgHiOVb6bSgsatThi3+VEbUkYbL8djNx0qpw54f+Zay3O+40Ks8bIh7NaLZhYrfRu7Ctnjp2nWW9OdqNeWWg+1J7pKLmW8jWB3SRo3d90kZSYe+aZZzB+/Hi43W643W6MHz8eTz/9dKrXNuCIV0CQEytnLuKUsJ5OR9t9qp5E8e7ObFYLThomjwFio5A+39eMQ95uFLjsmDFGPbLNiPPA3AS9VgraPA7RwWMu3Y6ImGMXayZImzr8SgJ+ZC1Ou1V1cbJZLSjNE8Os6rUebOniVnksmDM3ThBz1ZEDt607qDsNIFl8muakPa1obfcFVZ9/b4VZ81x2Pq6pod3Hc4q0sB5zrHq6KMfBG35q2+IEhT6JLDwT05mLXOzKBTHnsFnx0nWn4vn5U1OWK+SwKd+tWN957UVeLIBg4mF/Uyf+uPxr3jeRwcOsGuHFboJaOgMJ80/10I7rctmtMSdEiH+P4bBZ+Oej7XuohFnVYlnMW+1JmLWniO9r+dY63qqo0x9SpQJoL8givJo1hjOn3VcsZcCwM5enDrOqnTlBzEXWXeh2oEAbZmUhWJ3QLiuCYM/vCyo11aypKn7INDr9QX4+Zeen2RMHw2W3ojjPoXsO7EtMi7l77rkHt956Ky666CK8+uqrePXVV3HRRRfhZz/7Ge6+++7eWOOAQTs2RsSlk5wLKD2djrb5eKiqwGXnoVM9TmZ5c5Eh5f+3Sc4LO398ZdRdVb6OQBJJlKunPfGJVXpMtO3gzpx88ivhla5+flHNE8JT4l3qoHwXbFaL0ORYWWswFMa3//IJLvzzxzFzHCRJ4g6l6MzlOG18LSyMlCq0YdaeNg4W8+UAxU1JFeJMy9LIvtbOxGX4giEunI+JOHNWq0V35imrvLVYlGapsQogYjlzADBhqAenjkwueV6PAsF10qItgGA3GIGQ0uz4iQ9344kPduPlderG0+x96Ttz/qQaBgOyCBWdODHECsTvM2exWJRJCpqK1mYu5tTPUZy5oDAIve/bR7jsVj7f9bnVe1S/Y4U4QAJnzh59w6w3HotVqCYa5cVg3/uDXMwp52aLxcKrgo/GcOa6A2EEhJsdvUpVMW8u2TFaZmEOfX2bfFPHjnmzEzwyHfa9KXDb+XEztDgXG399Hh7/wZS0TqVgmBZzTzzxBJ566iksXLgQc+bMwZw5c7Bw4UL87W9/w1//+tfeWOOAQWlLEn0gxhrnxSq3WI7K0XY/79WlV8kqcspwpQgiEArzfmoXnxBd/p0fOZnFy0NrTZQzp7l4iGKOJa4zu55drJXpEIGopsGA+oLEGtkysSeG6ZoilbKt3UG8HsMWP+TthrcrALvVwkdEMdj+3RfH1UuG6DBrzxoHH9KEPGM1he4JrKAkz2WHzWrhIrcuRqh1f1MnwpIs/EThpdc4uFnIFWIXK+34MECZGzooPzXh1HjEc6PZjRcLr+Y77Xy2J9v+y0jXem3uohJmjc6ZA2Rhq6QsmM8B1HN9AE3OXAwHXa+ilX0W2nwtJnIOe7u5EEnHhdxisfDjfkOk3dKISNiWFeIA+j3mACHMKtww6zlzZ40tR3GuA98aPcjQ+liYVSmAiHZQAWU6SmGOIyrU3xonZw4ATqwphtNmjZkP2FuU5TthtUSazXf4+q0zd4BXsqrfl/YYTiemxVwoFMJJJ50U9fiUKVMQDKami/1Ahd/tx8yZi75rZM4cO3D8oTD/0hlpZzBhqAdOuxWNHX48v2YvmjsDKMt3YloMh4M3Do5TrdeSoPBCe3crug7aCwvLfWKCr6HNx8WPOIZHfM1BkfyNfE3iMKAOOf7r8wNRa9t2SA6xjirPjzrRMls91WJOOwOyp2FW5swxUZHqnDmlNYy8f+L1mjvUojRxFsN8XMxp2pOwm5CSPKdumBwQRnkl0YvMLPF6zfEZmRHhabVakO9UQq2hsMTzP7W5f9rpDwxVzhzvMWc+XCYKNpdDX8zFOk7FcWAiem4iO1d9ecjL30OqGjebRdyfE4Z4MHOc3EqIFeIA+j3mgNjnWL2cubPGlGPDPedhpsHqWzbPlX322kbD7HNigrjQ7YDDZuXbtfuCwjzg2AJiUIEL//rJNLz0o1MNrSkV2G1Wfkwf8fZfMcd6zFWXZG5Rh2kx98Mf/hBPPPFE1ON/+9vf8IMf/CAlixqoxKtCzeH5HNEjoErynPwAZ33SjOTauOw2nBDJs3j43R0AgNkTq2KGZ3mYVceZ8wfD/ESlPwFCE2YVkqm1F61B+eoCCHYSBtQDnwtjOHOxCiBEd2TLQS++rmtV/b1tMYofGOzEtC/Dw6zMmRsZcSRS7cwxd4aJrVgNQxksh6ZcqIoGoDsFolkYKh+vgEUvzNob6DlzwVCY30iJNxPi9nsaO7hY1zqMXTo5c0VCi55EleHxEEWNW3NjEi9nTv597Oa72h6PDBZF+DqSomB2jFcqEffn904aygtv2LD2eD3mACGVJZg4Zw6AqcIUbRVwtDMn/ywWQMj/ZY3RA1FFN7GYOLRI1US7LxArWhXR09/EXGxnLpNIyiN85pln8M477+DUU+U7gE8//RT79+/HlVdeidtvv51v99BDD6VmlQME7QxAEXbX6A+FEQpLsFkt/G45LxLKausOKmLOYHjm5BHF+GxPEz9RzIkRYgUQVSavt3aLRf/OUes4qgsg1OvlYdbI4+wkIfd7EsSc8LeYUxRLeGpdqn+tP4C7Z4/jP38ltCXR0ndh1tQ4cxOGeLC7oUO3BUiyaFvDxOs1xxLQKzSiq6wgkjPXphaaLMxanOvMGDFXqOPMieJO/K4XuB2Atxtt3UHeYBSI7neoVwDBck29nX6+P5KZaBHPmXOrwqxxnDmdMGuuVsxF1uyPVCfXpHF2JhOiTrsVcyZVYWfkXMicuXg95oDYYVZtbmSyaOe5RjlzPMwaceb4mDg7Gtp8aOsOqnrQZRLyedeL/U2dOORNbcPgTGF/HEc3UzD9rdi6dStOPFEetLt7924AwKBBgzBo0CBs3bqVb2fmroWQYQmu8VqTALKjk+eyc2cuz2lHWb4L3zR0cDFntDT95OElAOTPsbokB5M1zScZiUZ6sXmwnhwH72yuxUgBBCALQnbyY44dO5i0Cccxc+ZihOmYS1Wc60BzZwCvbTqIX1wwlvfSUipZlfJ+Rk0vhVl9MSZA9ATWY278EA9e23QIrd1B+INhPqu2p/Ccxcj+rYjTnoQ5cxUaZ26QrjPHEv71w6xhodginc4c+znXaVO52GJFK3OqAKBdIwZ5mNWhyZljYdauAE9Z6Kkzp3WAmNBzxuhlCYhhVvV7Vka5acKsGpEzPI2J76wf38xxFSjKdXKH+pC3C92BUNwec4BOAYSBnp1GKNXkeGodU3aMMseXO3PCjWlbnDScdMIc+g37miFJ8nesLN/8TUgmkw2Oo2kxt3LlSkPbHThwAOFwGFZrai4kAwGvTn4GoD74uwMh5Dhsyh2+y8Yvknv4EHtjB/yJw4phschTFuZMqtIV4TxnTkfMGRk9JJ6EClx2VVNSMcxakuvkv2OhWF7Jqr2YCK9ZrsmZa48RZr1gwmC882Udjrb78cH2Bpw3rgIdviD2RoRaLGeupkS+KBz2dqVUHEWN8+pxmFU+4Rw3uJA3km7q8HMHraeITjAAVHoiBRAxnDlFzKkvYoN0GgeLOWJ6s3VbugIIhtXzQXsTPTdaz60RxR+7OQCi30eXXxGDImLOXEsSc1kZOXEKIOQZua6Y33PxuXoFENqEb+25Kp2OzOSaInxx0MtnKpfkOeHJccDbFUDt0Y64IVZAaE2i6jOnnvSRLFFh1qicudgiWRmlGOA3+5nnzMnHIpsmVFOS2+/MHF4AkcFirteU1rhx47Bnz57eevl+SWucMKvVauEiojsYVgkB2ZmTTxas87vRGZWFbgdmHDsIuU4bvjulWne7WAJJxMgIMbHir0jT4kC8aImui7ZJqXZ8jrivyuPkzLE+WRUFblwyWZ7x96/P9+PrulZc+exnkCT5DlN7Bw3IFVs5DhvCktJaIBVoCyB6EmaVJIlPwagqyuH7LZWhVm2YNV4BRD0Ph8bOmYsqgGBh1jgFEOw5xbmOlAnqeOgVQOglxYv5TaKY074PvQIIVrna2h3g39dkZs3GE3N5Ljs++eVZeO6aU2I/VyfM2smEvFbMac5VyTQMThX/3/87DhvuOQ8n1sj9My0WC8+b+6ahI26POUAMsxrLmTODWxg/B+hXszKYYMsXbqLbUhTyTTXcoY9Ummey4EkGb5cipDM5zNprZ0SaDGGeeBMgAMAdOeC7/CF+p2yxyCchbdjJzMnniR9OwSe/OJuX8sciUYd7ftKL4yRYrRZ+ctI6DmIxhPhetA5jVJhHuKCykwq74KiduUgH+3wnF63vbqvHhX9ehc/3NiPPacO9Fyk5dCIWi4U7Dqkc68UuGnw4eg/CrN6uAH/+YI8iShtTNNJLkqSoJPh4UyBYCxGtM6dfACGGWWNPG+nLfDlA7DNnzpk70NylahOjLYBg80P1nDlJAg5EwjrJ5My544gG9pheKkSuXmsSX2JnHEhulFeqsFgsUQUarGH1Nw3txp25XsiZA9Sh1ihnTiPmPELOHMDEXOICiHSgdf77Xb5cJGpTmufMqFYkWigGmkHEmwABqEvnxTtlcZAzw0yujdthSzjMW2kaHLs1Cc/xSSAi2YkoSswJ6y0X3BztdtqTNTvp2a0WlES2jdU0mOXMleY5MaayABOHehAKSwiFJcw6vhLv3nEmLpgwWHfdrHdWvAkSZmHii7loXT3oM8dagZTkOeEWclbExsFd/hD+9tHupASpLxjmIU52QWcn8Q5/SOVehcMS6tti58yxdTVrRnqxMGtJnn41a0O7/Jp9J+Zi38DohbuYe71+rxxuYhdo7fvQq2Z12Kz8vbOCgqRy5gRnTpton/i50dWskiRFVTIzRGHhdlhVkzkyAe7MqcKsCZw5g9WsZhGjDNqcOa3o1oZZWzr9/DPJvDCrVsxlrnuVDNzRzXCRSmIug+D5GTp3XiwE4guGhBwW+TGtmEt1kmyiMKu3UymAiAc7EWlbkRTphFkdNqvq5KV3MRlU4OJuA9umOxBGMHJR1A4Kv2/O8Tj/+Ao8c9VJ+Ou8KRjsiX8C6o0iCBZmZQ2Te+LMHY7ky7ExW+x9iu1Jlm48gPvf/hp/emeH6dcXP3fmfOY67fxzFMPPzZ1+Hu7XCq/iXCfvli62i1Fakyhh1k5/CGFhTA535vqgYTCgLmgQadW56WLfxa8iPQtPiBQTyVMhxJnKsfPPgOjjJ5mcudwEzpyR54o3Fr5gGOxj0ApQt0OZODGsJC/jcqWO4WHWdiHvKfax7tIUQARCSrslvakLZhCLArQFGFqnjp1v2X9FpzfjxVw/nf5QncEhViDNYu6jjz7CRRddhKoqOfH+tddeU/1ekiTcd999qKqqQk5ODmbMmIEvv/xStY3P58PNN9+MsrIy5OXlYc6cOThwQN0Utrm5GfPmzYPH44HH48G8efPQ0tLSy+/OHH4hD04/zMpOtMpJhl34tBdNo61JjJKoAMJoXyx2AdSGjwrddn6R197di3e0Wmduck0xTh9VhvmRpGdAHQpiyecs3MiqZE+sKcaT807COcdVxF0vg4WPUjnSi13g2QW7Jzlz7GTPRCkL6RwVWrKwPK5kBKkiQNQhulg9+FjuTFm+U1XkAsihdiY0xby5WK1JAHWIMl1hVj1nTi/MyoTPScOL+e/EKvAunTAroD5+5IpT86foeK1JEj43RtNgUchrBajFYuH7IRMv4iMHsTCrCWcucpPFRDugP3XBDOJ5TBtWFX/OdynnQtZ6iU2O0FZQZwKFbrvqO9dfw6yZngvYa98KI3doHR0dmDRpEhYtWhTz9w8++CAeeughLFq0COvWrUNlZSXOO+88tLUpZf+33XYbli1bhpdffhmrVq1Ce3s7Zs+ejVBIORldccUV2LRpE5YvX47ly5dj06ZNmDdvXs/fZAoR5z/qzfsTE3TFiysAlGnFXBLhmXgkak2SaPoDo5A7c2oxZ7FYeIg2Spjm6ou5HKcNL1w3FdedMZI/5rLb+KzDdn8QgVCYi03toHCj9MYUCBZuU8KsPXDmIif7qqKIM5cf7cyxtjXawe9G0BY/MGLtFxZi1RY/MLR5c6GwxD+f4jwHXHYrv5iJlaB9Leb0+szpF0Cofz6+ysOPT/F96BVAANrei46knK5447wSEStnjqV05DhsMWdQshu0dDYM1mNYaS6sFrkiWekxF/t7qXXmeDhdEFc9QcyZi3LmBAdV3btQLeYyzZUD5HO3mDenJ5azlf0Jci0zhV77ZhgpgLjgggtwwQUX6D7/kUcewa9+9SvMnTsXAPDcc8+hoqICL730Em644QZ4vV4888wzeP7553HuuecCAF544QVUV1fj3Xffxfnnn49t27Zh+fLl+PTTTzF16lQAwFNPPYVp06Zh+/btGDNmTIrecc8QG0LqnTjcQodySzASUozcKWtL31M9UkcpkU9QzZogLDS2shDvbquP2RphaHEOGjv8UUnUJcIFTism9Mh329HU4Ud7dxCOyP60WpKrDgSUC9W+pk5IkpSScFK3xpnrWZhV7cyVRUSrGMrcHWmeyoZim7lAKX3GYreiEHMJ9YofGGUFLuCwIs5auwLczSrKcUbmbNrQ2h1UuUJ9OcoLSNxnTuvWFGgqrY8bXBjpBxlSvQ8lZy76uyw66sm66/H6zCUiVtNgvXw5BrtBS2fxgx4uuw1Di3P5zUZFoUt3n2jHeSXKYTZLqUFnTnR8WfU+qw7PtOIHRnmBC7VHO1BR2HdzYfuKbJj+APSiM/fVV19h2LBhST+/trYWdXV1mDlzJn/M5XLhzDPPxOrVqwEAn3/+OQKBgGqbqqoqjB8/nm+zZs0aeDweLuQA4NRTT4XH4+HbxMLn86G1tVX1rzeJN/2BoZxswkJHdhv/HTupOmyWmCGcnsATs4NhVf4Pw6gzd/t5x+KjO8+KOdPwfy+dhMe+PxkThqgb94rFGUbvTMWKSBZiLclz6lbxJWJIcQ4sFtmxSFWFKAvnpMKZO6jrzMkXAW9XgIunUFhSFUYYQW+cUyxnjveY03XmWNsUeT+yfLkCl523HIlVBMFEIhv11tswsdbpD/HcSyBxNSsgu1g1JbnK+xDCxZ06feYAwKOaV5zchVudM2fuFO+O0WdOr5KVMX1UGXKdNkwfVWZ2qX0Cy5sD4rtGSgGEOsyaMjFnMGculjPHvJFMdOYApRiqv4VYxRFwmR5mNfTNYM6YEZYuXQoAqK7W71lmhLq6OgBARYU6p6miogJ79+7l2zidThQXF0dtw55fV1eH8vLyqNcvLy/n28Ri4cKF+M1vftOj92AGnocT58TBTjZdgRBCkYuL2PeprMCF1u4gPBF3I5WIF/H27iBc+crJaPP+FnwVGbQ9OEGDWqvVoptbM6q8AKPKox07MSRr2JlzOQB0od0X5Hfa2p51ZnDZbRhc6MYhbzf2NXVGFZwkAxNvTKz2JGdOKYCQnTmlz5wslFiIlXGk1Rc1NzUe7ToX9FiFIUd4JWvsfaSdAsHakoi9B/NiiDkmElmz4t5GvHC2+4LcdU7UZw4AxlQWwGa1xGyzokyAiBZHxXFG3BnF3YOcOeYWdgZEZ07fSQSAX8wai9vPOzYqPzJTGDkoHyu3NwCIHyrjIxODYYTDEq/61kY9kkVsdB09ziv2vGmteEtF7l5vwMRcpgseszR2+NEVCMFi0Q/PZwqGroweT/SIo75CK0qMhLi028TaPtHr3HXXXao5s62trT0WqPFQxsbofyQ5vA9SiHcpF+/C2UiveK+RLDar7PaxkBHL/2jtDuCmf2xAICThgvGVOL4qelB9T4lXAKEH60fX4QvyFhg9EXOAnOB9yNuNfY2dvDFpT2AOZ0kPw6zhsMSH3TMxXcb7zPkgSRJ2a8RcXWs3JsD4ca3nzPEwa3MXwmEJVquFF0DoiUUWJmWhI1b8UBJDtDMR1OUP8RseMyK0JzhscgFCdyCMtm5FzDERWqIRW+KF97jB8nHAex5G1h4KSzGPXYYqzJqkM6duTdLzalbls9d/rUwVcoDSngQwJuYAuYL37S3yzf4Zo1PjOIrOXLymwYUxnDn+uwwb5cWYO3kovjrUih9MTT4al4lYAPxkxjFo7QqYTlnoawxdGRcvXtzb64iislIOw9XV1WHwYKX/V319PXfrKisr4ff70dzcrHLn6uvrMX36dL7NkSNHol6/oaEhyvUTcblccLn6rmeSuTBriF/4RaeKXSSTaTRqhAK3nP/DcoYkScIv//0F9jd1obokBw98Z2KvtCYojlMAoQfvi9cd5CGuWNMdzFBTkotPv2lKWRGE4sw5VD+b5WiHD4GQBItFuUNmF47ugFz5vKshWsyZQa8AYrDHDZvVAn8wjPo2Hyo9bi7StC0LGMdWyO7rqp0N8AVDqrYkDG14khVV5DhsvLK6LyhwO9Ad8PHQqiRJSn6iZli7eOyynFBtuFgU7DFbk6jCrEnmzPWgACInRgFEvFYq2QBrHAwkCLMK++pAcyfW1jYCAP5fnP6TZojvzIlhVkfM/5d/zszPYExlAZ6fPzXxhllGab4Lv5g1Nt3LMETG3k6NGDEClZWVWLFiBX/M7/fjww8/5EJtypQpcDgcqm0OHz6MrVu38m2mTZsGr9eLzz77jG+zdu1aeL1evk0mkGj6A6CIua5AiCeki3f3LHyV6uIHRr7GLXlh7T68vaUODpsFi75/Yq/9XTH0pFfpq0V0drQ95pKFjSpKVXsSlpvDHEN/KKzKzTIKG+M1KN/FHZJcp9IuoLHdz8Os7KIRa2pDPPTEnN1mxZCIqGEil1XL6jWQPW1UGSoL3WjuDOCdL4+o5rIylPCk/D0/IhRV9GUvM20RhFgVWakRq2IImjlz2v6MnZrJLVrEAp2knbkUF0Cw92v0RirTOMagM2e3WWGP5NQu23gQYUnuFZiq0KGqaXBUzpwYZlX2s3afZ2oBBJF+kjo6//Wvf+Gf//wn9u3bB79fnQy+YcMGw6/T3t6OXbt28Z9ra2uxadMmlJSUoKamBrfddhvuv/9+jB49GqNHj8b999+P3NxcXHHFFQDk8O/8+fNxxx13oLS0FCUlJViwYAEmTJjAq1uPO+44zJo1Cz/60Y/w5JNPAgCuv/56zJ49O2MqWQFgSk0xbjprVNwwpUvog9QZo8KsUtMwNtXkR04kR1q78fCKHXjig90A5JyZSZEGqb2BWAChnQ2phyg8Wd5YT8Os1TEqN5MlEAojFCnhFMN1XYEQCkyGrHi+nMYpKs134kBzF452+LA74sydPLwEq3Yd5flnRtELswKyY7mvqRP7mjpx0rDihM6czWrBpScNxZ/f34VX1u3HhKFyuFf1OWscLbbevgqxMrRTIJhwLsp1RLUWsdusOKG6CHXebn4cK+9DFkS8ktVhiylKRTdO21jbKDk9KIDIjUyAiNVnLtVFVX3FoAIXyvKdONru533n9HA7bGj3BbF0w0EAwOyJqXHlALlvYGmeE40d/iiHTc+Zy3XK7WDYuSJTnTki/Zj+Zvz5z3/Gr371K1x11VX4v//7P1xzzTXYvXs31q1bh5/+9KemXmv9+vU466yz+M8sR+2qq67CkiVL8POf/xxdXV248cYb0dzcjKlTp+Kdd95BQYGSJP/www/Dbrfj0ksvRVdXF8455xwsWbIENpty4nnxxRdxyy238KrXOXPm6Pa2SxdTR5Zi6sjSuNvkCGFWJSlZeZ/fOXEo6lt9+P4pvZPbx8Jbt/9zMz+5XDhhsKphb29QkkQ1qxjeYnNZUxFmBVLTa068WBbmOGCxyBVrXf6Q6SRnNspriCZBtzTfJc8JbeniAvS0UWVYtetoEmHW2IPWAXVFa2OHH6GwHPIVO95r+d5J1Xhs5S6s2nWUPyaG0wt0xJyeQOwtCrkzJzvnTDhrXTnGqz+ehrAkcUdMGy5moi5H56ZEdOM8ybYmSVHTYJZXHOvGMZuwWCx46sqTUN/m4y6yHm6HFe0+JQ0hVSFWxsK5E1B7tIO7/Ay91iRs3qyShpOdnwHR+5j+Zjz++OP429/+hu9///t47rnn8POf/xwjR47Er3/9azQ1NZl6rRkzZsTtR2exWHDffffhvvvu093G7Xbjsccew2OPPaa7TUlJCV544QVTa8tExNYknTHCXoMKXPi1zrD4VMAuTKGwhJqSXPzygrG4YHxlr4e9ipLoM9crYdaIaKlr7UZ3INSjfkqswtZikU/kOQ65uCSZIghtJSujLPJ+1+9pRliSLwTjh8iOUbLOXKz2FGKvOfa6ZfmuuJ3qq0tycfqoMny88ygXdLGcOeYKcbevj2d/asOsrNCkSkcUaAsBmPhlz+8KxHe5tE2Dk0GdM5dcmFWS5PNMjtPGBahea5JsYLLBgiVxf00ZVqz7OSdLrJZM2r8bNfNXFHMZWgBBpB/TOXP79u3juWY5OTl8GsO8efPwj3/8I7WrI1SwBF3RmTMadkwF355chTEVBbhn9jisuP1b+H8TBvdJ/lJpngs1JbmoLskxnJdXIOQqiX3mekJRroPPoWzqYa85X6THXE4k3Bar875R9kRy+LT5QKwI4tNv5ETuUeX53FFiOWhGYWIuljMqOpYNbUpuWyIuP7lG9XNxDNGebmeONQJWnDnWHsXYOvKEqmpA+Xz1xJwnJwUFECnoMwco+X3ZXgBhBjGPMZUh1kSIDqpWsMXqO0cQWkx/MyorK9HY2Ihhw4Zh2LBh+PTTTzFp0iTU1tYamvpAJA87SXcHQkrj0T68W541fjBmje+7ExzDZrXgnZ99C5IEw1MLRDHARlrFC/sZwWKxoCjXiYY2H5o7/T26a2cOHLt4xqoiNMqOI/INFasSZbCw8vbI748ZlM9zzrxdAVPuol4BBKAMLt8nOHPlOg2DRc4dV46SPCcXxiUxqlnbo3Lm0uvMcRfUoKjUVrPGG+UFyA7NsNJcNLT5ku5rJb62WffYZrXAZbfCF5SroEuFNedlac6cGdj+sliAC/rwXKfXmkT+2RHz/wlCxLQzd/bZZ+ONN94AAMyfPx8/+9nPcN555+Gyyy7DJZdckvIFEgriOK/OODlM/RG3w6Z7AYwFu4i2dAZ6PJdVhLlHzR2BBFvGh4VZmdsq5kOaocsf4jl8UWIu4kSye6xR5fmqodh1Jipa4410Ys5cQ5uPu4RGnDmX3Ya5k4fwn4ti9pmT9web/mBEJKYSlr/YysWcOWdOW83alcCZA4ClP5mOd372raQbxPakNYm4NnbDEU/I9zfYOfbkYSWGP+NUoA6zqj/3fJUzR2KOiI3po/Nvf/sbwmE5RPTjH/8YJSUlWLVqFS666CL8+Mc/TvkCCQV2wHf5Q8o4rwFwt5wMTMztb5bFRU/msoqwJH3WTiNZ2Cgvt5M5c5HO+yadud0N7ZAkWWRqnUftlIpR5fl8KHbt0Q4cae3G8DJ1IrYeerNZATk0WOC2o607iA17mwEYF12XnVyNp1fVAlA3Vc3XhCfrTYRvU0mBpgCiTjMDNxF5mgIIZfqD/qm3p4U6jkiLjWBYSqrRaa7TjubOAF9rthdAmIGlYlw0qW8jECpnLmqyCIVZicSY/mYcOHBANQnh0ksvxaWXXgpJkrB//37U1NTEeTbRE9xiaxKelEwHdyzY3SzLDSvOTX4uqwgTcy0JxFw4LEGCfliYh1kjF9tcPhMzGHN7PViIdXRFQVT+YqlG3I0ql9sysKHYZipaWZgxlhNssVhQU5KLLw+1YtOBFgDGc9tGVxTgrgvGotMfUj1HDJO3R/4B6WhNog2zsobBZsOsamHU2zdhl51cjX1NnRgSp6+aHkrIX12BOxDONb+YNQanHVOK75/St9cxvWpWQH0DRQUQhB6mj84RI0bg8OHDUfNOm5qaMGLECIRCyc+XJOKTI4RZediLnLmYaEWHVtgkC5vWwOaJ6nHX0i14a8thvPOzb8XMreNh1ohAF/MhzbCd58tF988SO8477Vbe/Z6Fj8xUtMbrMweAizl/pBGyGQfthjOPiXpMrGZlTYjzXfY+b1yr9JkLoK07wEWlXmsSLVHVrAbCrKngD5dMSPq5uZrv4kA61+jNh+5t2PHvslujQuPsO2i1DIzPgEgO02dGvZmm7e3tcLszexBttsPyObydAUTavCF3ANwtJ4P2ot/TSlZGkYEwaygs4Y0vDqHTH8KqXUdx6UnRff/YhTKHh1mTK4DYeURuBjymIvoCJArYkWV53CVkQqTOa6yiNRSWhPFxsS8mNZou+T3NbRMLB5RZr30bYgXEPnNBHmItdNsNu1RRBRCB+EPrM4Ech/q7qEybydw1ZztDinJw9fThGFqcE3V9Ze5wgdvRp9NPiOzC8NHJGvpaLBbcc889yM1VTt6hUAhr167FCSeckPIFEgpMzDUJQiKnB73O+jPasV+lKSh+AJQCiJY4ztw3De38Qrgz4pxp6dYNs5oTc2KYNXqtipg7plxx7lio8kibMWeuQwj96okY7cijnua2iZMT2FzWij4ufgDUEyAOm8yXAxTx2xUIyaK4j5y5nqC9sWDh1mwd55UNWCwW3Dfn+Ji/U8Qc7X9CH8Pfjo0bNwKQnbktW7bA6VQuFE6nE5MmTcKCBQtSv0KCw4Qbq07McdgMt+oYaGgdpFSFWZkzF6/P3BcHvPz/d9a3x9yGF0BoWpOYCbN2+II40Cy3ytBWsgJyaNWT44C3K4BRwhgj3mvOYDUrc5XskbYVsRCdOaul50n8+REXyB8K8/fY18UPgLoAQhmbZlxUijcVHf4gF0ZmKrP7GnE+azgsKb3xsrhpcDYjOnMEoYdhMbdy5UoAwDXXXINHH30UhYX6M0SJ3kE7mDubO7L3Ni67DU6bFf6Qeph9TzFSAPFFpAgAUMKgWlLRZ44JxbJ8l+77K813ymJOcOYqPbIoMloAwfPl3HbdMI/ozA0qcPX4JkP8brO5sn3dMBhQLqQd/hAONrNJG8bX4bLb4LBZEAhJ6PAFY47hyzRYpW2nP8TDwsDAaYOUaYyO5PCNiZEXSxAM00fn4sWL+f8fOHAAFosFQ4YMifMMIlW4NCFVymGJT57LBn+nLOZ6OsqLUWKgAOKLg4ozd7ClC+2+YFSISlsAkStcQI2yI07xA+Pq6cPx1heHceaYQfwxls9W3+rTzYEViTeXlTGkKIfPl02F6LLbrLx57TcNHfK60yLmFDeEiefKQnMVonkuO1o6A2jvDmZFmFVx5oJ8bKDVEn0zSfQN44d48PHPz0rLzQyRPZg+OsPhMH7729/C4/Fg2LBhqKmpQVFREX73u9/x/nNE76A9mWbyBSETEENcPQ37MRIVQARCYXx1qBUA+OivXTFCrd3COC8AyHEqo9qMslNn8oPIldOG45UbpqnaHbCLgj8UTliVCySuZAXkkG5VJJesPEXzU9nf+4Y7c30fZnUK1YVMPJtx5gBFBLf7xDBr5t6IiaPleMNgp74rS/Q+1SW5cCbRAJoYOJj+dvzqV7/CokWL8MADD2Djxo3YsGED7r//fjz22GO45557emONRASnzQoxejUQ+j71BNFJSnWYta07iGAo+uZl55F2+IJhFLjtOGlYCQBFBIh0R4VZmTNnvM/cjkgIN56Yi4XTbuVOpZEpEMoEgPg3D2ysV6ocNPb9ZtMX0uVMMHeOTbcwkzMHqHvNcWcugwuXcoQJEHyUF51rCCKjMS3mnnvuOTz99NP4yU9+gokTJ2LSpEm48cYb8dRTT2HJkiW9sESCYbFYVLMWyZmLj1j9laowqyfHAWZQtHRFu1osX27CEA/GVMoiK1ZFa1SfOQe7gBp3t42EWfUwU9HazhoGJ7igHxMpshiaRKPaWGj/XjqqWQGlPUko0g/ItDMXEcGyM5dNYVbFmaPiB4LIbEzfbjU1NWHs2LFRj48dOxZNTU0pWRShj9thEwZf091yPEQxkCpnzma1oNAtV4g2d/ijRmaxfLkJQz2ojjTpjVXRqnXmxDwlI7R2B3irjFhtSRJRWejCtsPGKlo7DLam+OlZo1BVlIMrUtQ9v0Dz99LRZw6IbglRaaI1CQDkR5y9DkHMZXI1q9hnjo/yonMNQWQ0pp25SZMmYdGiRVGPL1q0CJMmTUrJogh9xL5yFPqIDxMfVot6iHtPYb3mYuWbbYm0JZk0tIiHP2NVtPa0mpW9ZkWhC54kRvywKRBGKlqNDlqvKsrBT88albJ9LYZ1PTkOlSvdl4hFEAVJTKHgc2b9Qf65Z3LxEg/5B0LCKK/MFZ8EQSThzD344IO48MIL8e6772LatGmwWCxYvXo19u/fj7fffrs31kgIuIQiCDrBxodddItznSntx1ec58Sexs6oIghfMISv6+TihwlDPPzvx6pojeoz51DylIxgpPghHiz/zMhILyMFEL2BKB7TUfzAEJ05s/lygHqkV1/NZu0Jokvc4SNnjiCyAdPO3JlnnokdO3bgkksuQUtLC5qamjB37lxs374dZ5xxRm+skRBgEwOAzL67zwSYGEhViJWh12vu68NtCIQkFOc6MLQ4B8V5Th6G1Va0MtGWExVmNSbmtqdMzCUe6ZUudyZfJebS15ZBFHNmQ6yA8j1s6w5yEZ/JYk50iTuoAIIgsgLTR+i+fftQXV2NP/zhDzF/V1OTmnwZIjZiexIauhwfJgZSNf2BUaQTZlXy5Yp4G4djK/JxtN2HHUfacEJ1Ed/Wp+0zZyDMKvaE28krWZNrJKrMZzUeZs139W0HelFA9HTWa08Qw6yDkxCV7Ht4tF0Rzpl8I8Yqbbv8Id5njqIABJHZmD6jjBgxAocPH0Z5ebnq8cbGRowYMQKhkLnZkoQ5xMTpXLpbjgsTXYNSLASKdXrNbYlUsk4a6uGPHVtRgNW7G6MqWrVhVneCMOsXB1rwg6fWIhAOoyTXiaPtfv76yZBcmLVvL+iZGGatNFnJCijvo6FNFnOWDG/Am8vb5ITQzsPCdK4hiEzG9BGq1zG+vb0dbjd1qO5txDArOXPxmT2xCjuOtOMHU1PrFvMCCM18VjaTdcIQRcyxMVrailZtAQS7WPqDYYTCUlSO3+ubDqEtIqoORdy0Ape9B2JOFkeNHX74g+G4DUmNFkCkGlE8pjfMKjhzSYg59j7qI2Iux2HL6Aa8qj5zPgqzEkQ2YPgIvf322wHIvc7uuece5OYqsxhDoRDWrl2LE044IeULJNSo+szRCTYugwpcWDh3QspfV5kCoYRZu/whLtgmDi3ij+tVtEaN8xKEeac/GDVU+/N9zQCAuy88DicNL0Fzhx8jyvKSvsiW5DmR77Kj3RfEzvo2HF/l0d02XWIuE525wUXmc+bYJBLmzGVyvhygzt9UCiAye80EMdAxfHbeuHEjANmZ27JlC5xOJQ/J6XRi0qRJWLBgQepXSKgQxRydYNNDrAKIbXWtCIUllOW7VMJjdMSZ01a0dmsKIFx2K59t2hUIqcRcdyCErZF8vJnjKlFTqtxIJYvFYsGJw4rx0Y4GrKtt0hVzLZ1+Pp4sVc2AjSIWQKRjLiujUBRzyYRZI65rU0fEmcvw45Z9J/2hMFq75RsWunEkiMzG8BG6cuVKAMA111yDRx99FIWFhb22KEIfMdeG8ljSQ3FedAEEGwY/pjJfFUJjFa1H233YVd/OiyC0OXMWiwU5kYbQ2orWLQe9CIQkDCpw8ZFZqeCU4RExt6cZV582IuY2/95wEL5gGMcNLsS4wX17zIvtMFI17zUZRGGdTM4cE6WRARLIdWT2cSuKzYZIbmZf50sSBGEO01m4ixcvJiGXRlTOHJ1g00IsZ44Ngx9Rlhe1Pas4ZeO3QmEJ/pBazAFCeEtTBPH5XjnEOqWmOKW5VicPl2fHfranCZIkRf1ekiS8uHYvAOAHU2v6PM9LDLMOSqOYK4yIuTynLWoqhRG04elMd+ZcdmUG9FEeGs5sAUoQA53MLakiYpKjuvjTCTYdFAs5c0wE1R6VnbmRZdGtQpS8OVnM+YKKWBM/T7cwRklk/Z6ImBtWnJL1MyZVF8Fps6KhzYe9kSHyImtrm/BNQwdynTZcfEJVSv+2Edhki7J8J1z29AmgsYML8K1jB2H+6SOSErRaMZfpN2EWi4WfWxoi7VT6umE0QRDmoCM0y3DTBIi0w1qehMISWruD8OQ4eJh1xKBoZ+6YSN7c7sg2YhjVZRfD5tGNgyVJwoZI8cOU4akVc26HDROGevD53mas29OE4RpX8cW1+wAAF58wJKogoy84bnABfnTGCIwfol+c0Rc4bFb8/dpTkn6+VgjlZHiYFZDdw3ZfEP5g5jc5JgiCnLmsw03OXNpxO2zcUWvp9CMcllDbKAu1Y2I4c8NK5IKF/U2y+9UduUA67VZYhRYkbCamKOb2NHaiqcMPp92K46tSn97AQq3r9jSpHj/a7sPyrYcBIOWtXYxisVjwqwvH4eIThqTl76eKfLf6OM0GYaRdI7UmIYjMhsRcluGKkWNF9D3FwhSIgy1dcq82mxVDYlR8VjMx19wJSZK4WMvRDI7PibiunULO3PqIyJo4xNMrocZTRshu37pIKJfxr88PIBCSMGmoJ+3OWLaTq/mcs+G41X43ScwRRGZDYi7LYCdZp90Kh40+vnRRJEyBYPlyw0pzo5r9AsCQohxYLHIFa0O7L6rHHCOXO3NB/lhvhVgZU4aVwGKRc/7q2+RmxOGwhJciIdYfTB3WK393IGG1WlRthDK9AAKIXiO1QSKIzIbUQJbBBACdXNNLSZ5S0RqvkhWQhXdVZED7/qYuXgDhjnLmonPmxErW3sCT48CYSIEGK7T4x7p92NfUiQKXHbMnDe6VvzvQEJ2tbHDmtGuklA6CyGxIzGUZbJwXhT3SCyuCaOoIKJWsg/SH3rOGu/ubOtHll3PmosKskQsoC7N6OwPYEZkccWKKK1lFThkRaVFS24Q9Rzvw+ze3AQBuPXc0XcRTRL5KzGX+PhWLNJw2a9xxbwRBpB86QrOM4ogjVJrnTLAl0ZuIvea+4W1JYjtzAFAjFEGwMKtLJ5eKOXMb9stO2YiyPJTl916fNVYE8ek3jbj9n5vQFQjh1JEluFankTBhHvHmSyviMxHRmculqnmCyHgy/xaRUDG5ugj3XjQOJ/ZS2I0whlIA4Y/bloTBiiD2NXXy7XI0OXPaMOuGSIi1tz9rJua+rpP74BW47PjT9yapKm2JnpGfxWHWvCxwEglioENHaZZhtVpwDTkmaYcVQNR5u3HI2wXAmDO3r6mTi7WonDlNmHX17kYAqW8WrKXS40Z1SQ72N8nv4945x2Nocc/nvxIKqpy5LEiREAsgqJ8lQWQ+FGYliCRg81k37W+BJMnD2EvihL6ZM3eguYv3mXPb9cOsn+w6is/3NsNhs+DMMYN64y2omD6yDAAwc1wFvnNidvd1y0TE2abaViWZiCrMSs4cQWQ8dJQSRBKwnLmjkUHkIwflxx31VF0iF0Ac8nahrTsAILr9AwuzdviCeOA/XwOQW4MMKYruXZdqFpw/BuOHenDJ5CF9PoN1IJBt1axiXh+N8iKIzCfjnblgMIi7774bI0aMQE5ODkaOHInf/va3CIfDfBtJknDfffehqqoKOTk5mDFjBr788kvV6/h8Ptx8880oKytDXl4e5syZgwMHDvT12yH6CUzMMeKFWAFgUL4LbocVkgTsrpdz7LR95tgEiFW7jmLLQS/yXXbcfPaoFK46zvoKXJh36jC6cPcS4n7Njj5z2SU+CWKgk/Fi7o9//CP++te/YtGiRdi2bRsefPBB/M///A8ee+wxvs2DDz6Ihx56CIsWLcK6detQWVmJ8847D21tbXyb2267DcuWLcPLL7+MVatWob29HbNnz0YoFIr1ZwkiLlFiLk7xAyCPpmJ5czvr5e+ldqIDu2h2RnLqrv/WSJT2YhUr0XfkZVlrElHAkcAniMwn48XcmjVrcPHFF+PCCy/E8OHD8d3vfhczZ87E+vXrAciu3COPPIJf/epXmDt3LsaPH4/nnnsOnZ2deOmllwAAXq8XzzzzDP73f/8X5557LiZPnowXXngBW7ZswbvvvpvOt0dkKUV56sHzI2LMZNVSHSkq2BnpHacXZgWAsnwX5p9OhS79hWyuZqXWJASR+WS8mDv99NPx3nvvYceOHQCAzZs3Y9WqVfh//+//AQBqa2tRV1eHmTNn8ue4XC6ceeaZWL16NQDg888/RyAQUG1TVVWF8ePH8220+Hw+tLa2qv4RBKPAZYddaN2RyJkDlCKILjbOyx67mhWQG/ZSY+j+Q9aFWR3UmoQgsomMP0p/8YtfwOv1YuzYsbDZbAiFQvjDH/6A73//+wCAuro6AEBFRYXqeRUVFdi7dy/fxul0ori4OGob9nwtCxcuxG9+85tUvx2in2CxWFCU6+AFEMNLjYs5hjZnbnhpHhw2C0aU5eHyk6tTt1gi7WRbAYQYCqabCoLIfDL+KH3llVfwwgsv4KWXXsLxxx+PTZs24bbbbkNVVRWuuuoqvp22Ak+SpIRVefG2ueuuu3D77bfzn1tbW1FdTRdYQqEo14mj7X4MKcox5LbUaMSc9jmVHjfev2MGPLkOOGwZb5oTJhB7tWkd2UwkR9WaJPPXSxADnYwXc3feeSd++ctf4vLLLwcATJgwAXv37sXChQtx1VVXobKyEoDsvg0erAwFr6+v525dZWUl/H4/mpubVe5cfX09pk+fHvPvulwuuFyUfE7oUxIpghiRoJKVoRVzsS7qWveO6B+wMGuOw5YVkzWoAIIgsouMv/3v7OyE1apeps1m461JRowYgcrKSqxYsYL/3u/348MPP+RCbcqUKXA4HKptDh8+jK1bt+qKOYJIRFFkpJeRfDkAGFqs7hfnJsdjwDCoQL4xLM3PjpnKYs5cNkysIIiBTsYfpRdddBH+8Ic/oKamBscffzw2btyIhx56CNdeey0AObx622234f7778fo0aMxevRo3H///cjNzcUVV1wBAPB4PJg/fz7uuOMOlJaWoqSkBAsWLMCECRNw7rnnpvPtEVnMmMoCvPPVEcPjtvJcdpTlO3mendue8fdSRIoYVpqH//3eJAwrzQ7nVT2blW46CCLTyXgx99hjj+Gee+7BjTfeiPr6elRVVeGGG27Ar3/9a77Nz3/+c3R1deHGG29Ec3Mzpk6dinfeeQcFBQV8m4cffhh2ux2XXnopurq6cM4552DJkiWw2ehERSTHreeMxkWTqjC6PHFbEsbQ4lxFzGXBWCcidXxnytB0L8Ew6tmsGX+ZIIgBj0WSJCndi8gGWltb4fF44PV6UVhYmO7lEFnKLf/YiNc3HwIAvPrjaTh5eEmaV0QQ0YTCEo75/94GALxx0+mYMNST5hURxMDDjO6gOA9B9CFiEUQ2VDUSAxOb1QJXJA2AmgYTROZD/jlB9CHVJUoRhLbPHEFkElefNhzfNHRghIEeigRBpBcScwTRh4itRyhnjshk7rrguHQvgSAIg5A1QBB9CJvPCpCYIwiCIFIDOXME0YcM9rhRUeiCPxiGJ8eR7uUQBEEQ/QAScwTRh9htVrx9yxkIhSU4qc8cQRAEkQJIzBFEH1OaT2PiCIIgiNRB1gBBEARBEEQWQ2KOIAiCIAgiiyExRxAEQRAEkcWQmCMIgiAIgshiSMwRBEEQBEFkMSTmCIIgCIIgshgScwRBEARBEFkM9ZkziCRJAIDW1tY0r4QgCIIgiP4O0xtMf8SDxJxB2traAADV1dVpXglBEARBEAOFtrY2eDyeuNtYJCOSj0A4HMahQ4dQUFAAi8WS7uXo0traiurqauzfvx+FhYXpXk7aof2hhvaHGtofamh/qKH9oYb2h5re3h+SJKGtrQ1VVVWwWuNnxZEzZxCr1YqhQ4emexmGKSwspINNgPaHGtofamh/qKH9oYb2hxraH2p6c38kcuQYVABBEARBEASRxZCYIwiCIAiCyGJIzPUzXC4X7r33XrhcrnQvJSOg/aGG9oca2h9qaH+oof2hhvaHmkzaH1QAQRAEQRAEkcWQM0cQBEEQBJHFkJgjCIIgCILIYkjMEQRBEARBZDEk5giCIAiCILIYEnMZxkcffYSLLroIVVVVsFgseO2111S/P3LkCK6++mpUVVUhNzcXs2bNws6dO6NeZ82aNTj77LORl5eHoqIizJgxA11dXfz3zc3NmDdvHjweDzweD+bNm4eWlpZefnfmScX+qKurw7x581BZWYm8vDyceOKJ+Ne//qXaJlv2x8KFC3HyySejoKAA5eXl+Pa3v43t27ertpEkCffddx+qqqqQk5ODGTNm4Msvv1Rt4/P5cPPNN6OsrAx5eXmYM2cODhw4oNomG/ZJKvZHU1MTbr75ZowZMwa5ubmoqanBLbfcAq/Xq3qdgbI/tNtecMEFMY+9gbY/+sM5NVX7o7+cU43sj6VLl+L8889HWVkZLBYLNm3aFPU6GXE+lYiM4u2335Z+9atfSf/+978lANKyZcv478LhsHTqqadKZ5xxhvTZZ59JX3/9tXT99ddLNTU1Unt7O99u9erVUmFhobRw4UJp69at0o4dO6RXX31V6u7u5tvMmjVLGj9+vLR69Wpp9erV0vjx46XZs2f35Vs1RCr2x7nnniudfPLJ0tq1a6Xdu3dLv/vd7ySr1Spt2LCBb5Mt++P888+XFi9eLG3dulXatGmTdOGFF0a93wceeEAqKCiQ/v3vf0tbtmyRLrvsMmnw4MFSa2sr3+bHP/6xNGTIEGnFihXShg0bpLPOOkuaNGmSFAwG+TbZsE9SsT+2bNkizZ07V3r99delXbt2Se+99540evRo6Tvf+Y7qbw2U/SHy0EMPSRdccEHUsSdJA2t/9Jdzaqr2R385pxrZH3//+9+l3/zmN9JTTz0lAZA2btwY9TqZcD4lMZfBaE+g27dvlwBIW7du5Y8Fg0GppKREeuqpp/hjU6dOle6++27d1/3qq68kANKnn37KH1uzZo0EQPr6669T+yZSSLL7Iy8vT/r73/+ueq2SkhLp6aefliQpe/eHJElSfX29BED68MMPJUmSBW5lZaX0wAMP8G26u7slj8cj/fWvf5UkSZJaWlokh8Mhvfzyy3ybgwcPSlarVVq+fLkkSdm7T5LZH7H45z//KTmdTikQCEiSNDD3x6ZNm6ShQ4dKhw8fjjr2Btr+6K/n1GT3R389p2r3h0htbW1MMZcp51MKs2YRPp8PAOB2u/ljNpsNTqcTq1atAgDU19dj7dq1KC8vx/Tp01FRUYEzzzyT/x6QwwUejwdTp07lj5166qnweDxYvXp1H72bnmNkfwDA6aefjldeeQVNTU0Ih8N4+eWX4fP5MGPGDADZvT9YKLCkpAQAUFtbi7q6OsycOZNv43K5cOaZZ/L38vnnnyMQCKi2qaqqwvjx4/k22bpPktkfeq9TWFgIu10eXz3Q9kdnZye+//3vY9GiRaisrIx63YG0P/rzOTXZ70d/Padq94cRMuV8SmIuixg7diyGDRuGu+66C83NzfD7/XjggQdQV1eHw4cPAwC++eYbAMB9992HH/3oR1i+fDlOPPFEnHPOOTyXrK6uDuXl5VGvX15ejrq6ur57Qz3EyP4AgFdeeQXBYBClpaVwuVy44YYbsGzZMhxzzDEAsnd/SJKE22+/HaeffjrGjx8PAHy9FRUVqm0rKir47+rq6uB0OlFcXBx3m2zbJ8nuDy2NjY343e9+hxtuuIE/NtD2x89+9jNMnz4dF198cczXHkj7o7+eU3vy/eiP59RY+8MImXI+tafkVYg+weFw4N///jfmz5+PkpIS2Gw2nHvuubjgggv4NuFwGABwww034JprrgEATJ48Ge+99x6effZZLFy4EABgsViiXl+SpJiPZypG9gcA3H333Whubsa7776LsrIyvPbaa/je976Hjz/+GBMmTACQnfvjpptuwhdffKFyCBjadRt5L9ptsm2fpGJ/tLa24sILL8S4ceNw7733xn2NeK+TCSS7P15//XW8//772LhxY9zXHyj7o7+eU3tyvPTHc2q8/ZEMfX0+JWcuy5gyZQo2bdqElpYWHD58GMuXL0djYyNGjBgBABg8eDAAYNy4carnHXfccdi3bx8AoLKyEkeOHIl67YaGhqg7skwn0f7YvXs3Fi1ahGeffRbnnHMOJk2ahHvvvRcnnXQS/vKXvwDIzv1x88034/XXX8fKlSsxdOhQ/jgLiWnv9urr6/l7qayshN/vR3Nzc9xtsmmf9GR/MNra2jBr1izk5+dj2bJlcDgcqtcZKPvj/fffx+7du1FUVAS73c5Dzd/5znd4GG0g7Y/+eE7tyf7oj+dUvf1hhEw5n5KYy1I8Hg8GDRqEnTt3Yv369TwcMnz4cFRVVUWVV+/YsQPDhg0DAEybNg1erxefffYZ//3atWvh9Xoxffr0vnsTKURvf3R2dgIArFb1V91ms/E77mzaH5Ik4aabbsLSpUvx/vvvc9HKGDFiBCorK7FixQr+mN/vx4cffsjfy5QpU+BwOFTbHD58GFu3buXbZMs+ScX+AGRHbubMmXA6nXj99ddVeZjAwNofv/zlL/HFF19g06ZN/B8APPzww1i8eDGAgbU/+tM5NRX7oz+dUxPtDyNkzPk0JWUURMpoa2uTNm7cKG3cuFECID300EPSxo0bpb1790qSJFfZrVy5Utq9e7f02muvScOGDZPmzp2reo2HH35YKiwslF599VVp586d0t133y253W5p165dfJtZs2ZJEydOlNasWSOtWbNGmjBhQsaVjUtSz/eH3++XRo0aJZ1xxhnS2rVrpV27dkl/+tOfJIvFIr311lt8u2zZHz/5yU8kj8cjffDBB9Lhw4f5v87OTr7NAw88IHk8Hmnp0qXSli1bpO9///sxW5MMHTpUevfdd6UNGzZIZ599dsxS+kzfJ6nYH62trdLUqVOlCRMmSLt27VK9zkDcH7GATmuSgbI/+ss5NRX7oz+dU43sj8bGRmnjxo3SW2+9JQGQXn75ZWnjxo3S4cOH+TaZcD4lMZdhrFy5UgIQ9e+qq66SJEmSHn30UWno0KGSw+GQampqpLvvvlvy+XxRr7Nw4UJp6NChUm5urjRt2jTp448/Vv2+sbFR+sEPfiAVFBRIBQUF0g9+8AOpubm5D96hOVKxP3bs2CHNnTtXKi8vl3Jzc6WJEydGldVny/6ItS8ASIsXL+bbhMNh6d5775UqKysll8slfetb35K2bNmiep2uri7ppptukkpKSqScnBxp9uzZ0r59+1TbZMM+ScX+0PuOAZBqa2v5dgNlf+i9rlbMDbT90R/OqanaH/3lnGpkfyxevDjmNvfeey/fJhPOp5bIGyIIgiAIgiCyEMqZIwiCIAiCyGJIzBEEQRAEQWQxJOYIgiAIgiCyGBJzBEEQBEEQWQyJOYIgCIIgiCyGxBxBEARBEEQWQ2KOIAiCIAgiiyExRxAEQRAEkcWQmCMIgkiSDz74ABaLBS0tLeleCkEQAxiaAEEQBGGQGTNm4IQTTsAjjzwCQB5C3tTUhIqKClgslvQujiCIAYs93QsgCILIVpxOJyorK9O9DIIgBjgUZiUIgjDA1VdfjQ8//BCPPvooLBYLLBYLlixZogqzLlmyBEVFRXjzzTcxZswY5Obm4rvf/S46Ojrw3HPPYfjw4SguLsbNN9+MUCjEX9vv9+PnP/85hgwZgry8PEydOhUffPBBet4oQRBZBzlzBEEQBnj00UexY8cOjB8/Hr/97W8BAF9++WXUdp2dnfjzn/+Ml19+GW1tbZg7dy7mzp2LoqIivP322/jmm2/wne98B6effjouu+wyAMA111yDPXv24OWXX0ZVVRWWLVuGWbNmYcuWLRg9enSfvk+CILIPEnMEQRAG8Hg8cDqdyM3N5aHVr7/+Omq7QCCAJ554AscccwwA4Lvf/S6ef/55HDlyBPn5+Rg3bhzOOussrFy5Epdddhl2796Nf/zjHzhw4ACqqqoAAAsWLMDy5cuxePFi3H///X33JgmCyEpIzBEEQaSQ3NxcLuQAoKKiAsOHD0d+fr7qsfr6egDAhg0bIEkSjj32WNXr+Hw+lJaW9s2iCYLIakjMEQRBpBCHw6H62WKxxHwsHA4DAMLhMGw2Gz7//HPYbDbVdqIAJAiC0IPEHEEQhEGcTqeqcCEVTJ48GaFQCPX19TjjjDNS+toEQQwMqJqVIAjCIMOHD8fatWuxZ88eHD16lLtrPeHYY4/FD37wA1x55ZVYunQpamtrsW7dOvzxj3/E22+/nYJVEwTR3yExRxAEYZAFCxbAZrNh3LhxGDRoEPbt25eS1128eDGuvPJK3HHHHRgzZgzmzJmDtWvXorq6OiWvTxBE/4YmQBAEQRAEQWQx5MwRBEEQBEFkMSTmCIIgCIIgshgScwRBEARBEFkMiTmCIAiCIIgshsQcQRAEQRBEFkNijiAIgiAIIoshMUcQBEEQBJHFkJgjCIIgCILIYkjMEQRBEARBZDEk5giCIAiCILIYEnMEQRAEQRBZzP8PmVNhadTzgXAAAAAASUVORK5CYII=", "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": 4, "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": 5, "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.08898 ... 135.1\n",
       "    total_precip_upper  (station_num, dparams) float64 120B 0.3218 ... 173.4
" ], "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": 7, "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: 164B\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) int32 4B 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: 164B\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) int32 4B 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": 7, "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": 8, "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": 8, "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": 9, "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) int32 4B 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) int32 4B 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": 10, "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": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAFICAYAAACP747yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEr0lEQVR4nO3deVxU9f4/8NeI7AMDiGyBDGIoqKCZirjnghqgoolXI9Tcb66oXbotWj8zMxSXTCuDXL6ZqZhlUkquaW6JhpKlF3eJUhwcQED4/P7gcq7jDDijgzM6r+fjMY/HzOd8zjnvMzLOez7ns8iEEAJEREREJlbP1AEQERERAUxKiIiIyEwwKSEiIiKzwKSEiIiIzAKTEiIiIjILTEqIiIjILDApISIiIrPApISIiIjMQn1TB/C4qKysxNWrV+Hk5ASZTGbqcIiIiB4bQgjcunULPj4+qFev5vYQJiV6unr1Kvz8/EwdBhER0WPr0qVL8PX1rXE7kxI9OTk5Aah6Q52dnU0cDRER0eOjsLAQfn5+0ndpTZiU6Kn6lo2zszOTEiIiogdwv+4P7OhKREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGaBM7oSEZFZEUKguLwYAOBg7cBFUC0IW0qIiMisFJcXQz5PDvk8uZSckGVgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGaBSQkRERGZBSYlREREZBaYlBAREZFZYFJCREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGaBSQkRERGZBSYlREREZBaYlBAREZFZYFJCREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGaBSQkRERGZBSYlREREZBaYlBAREZFZYFJCREREZoFJCREREZkFJiVERERkFpiUEBERkVlgUkJERERmgUkJERERmQUmJURERGQWmJQQERGRWWBSQkRERGaBSQkRERGZBSYlREREZBZMmpTs3bsX0dHR8PHxgUwmw5YtWzS2q9VqvPLKK/D19YW9vT2Cg4Px0UcfadQpLS3FpEmT4O7uDkdHR8TExODy5csadQoKChAfHw+FQgGFQoH4+HjcvHmzjq+OiIiIDGHSpKSoqAhhYWFYtmyZzu3Tpk1DRkYG1q5di5ycHEybNg2TJk3C119/LdWZOnUq0tPTsX79euzfvx9qtRpRUVGoqKiQ6gwbNgxZWVnIyMhARkYGsrKyEB8fX+fXR0RERPqrb8qT9+3bF3379q1x+8GDB5GQkIBu3boBAMaOHYuVK1fi6NGj6N+/P1QqFVatWoU1a9agZ8+eAIC1a9fCz88PO3fuRGRkJHJycpCRkYGff/4Z7du3BwB88skn6NChA86cOYOmTZvW+XUSERHR/Zl1n5JOnTph69atuHLlCoQQ2LVrF37//XdERkYCAI4dO4by8nL07t1b2sfHxwctWrTAgQMHAFQlNgqFQkpIACA8PBwKhUKqo0tpaSkKCws1HvRkEEKgqKgIRUVFEEKYOhwiIvovs05KlixZgpCQEPj6+sLGxgZ9+vTB8uXL0alTJwBAXl4ebGxs4OrqqrGfp6cn8vLypDoeHh5ax/bw8JDq6DJv3jypD4pCoYCfn58Rr4xMqbi4GHK5HHK5HMXFxaYOh4iI/svsk5Kff/4ZW7duxbFjx5CcnIyJEydi586dte4nhIBMJpNe3/28pjr3SkpKgkqlkh6XLl168AshIiKi+zJpn5LalJSU4LXXXkN6ejqef/55AEBoaCiysrLwwQcfoGfPnvDy8kJZWRkKCgo0Wkvy8/MREREBAPDy8sKff/6pdfy//voLnp6eNZ7f1tYWtra2Rr4qIiIiqonZtpSUl5ejvLwc9epphmhlZYXKykoAQJs2bWBtbY0dO3ZI269du4bs7GwpKenQoQNUKhUOHz4s1Tl06BBUKpVUh4iIiEzPpC0larUaZ8+elV7n5uYiKysLbm5uaNSoEbp27YqZM2fC3t4e/v7+2LNnD1avXo2FCxcCABQKBV5++WUkJiaiQYMGcHNzw4wZM9CyZUtpNE5wcDD69OmDMWPGYOXKlQCqRvFERUVx5A0REZEZMWlScvToUXTv3l16PX36dABAQkIC0tLSsH79eiQlJWH48OG4ceMG/P39MXfuXIwfP17aZ9GiRahfvz6GDBmCkpIS9OjRA2lpabCyspLqrFu3DpMnT5ZG6cTExNQ4NwoRERGZhkxwTKReCgsLoVAooFKp4OzsbOpw6CEUFRVBLpcDqGqtc3R0NHFERHS3orIiyOf99zOapIajDT+jjzt9v0PNtk8JERERWRYmJURERGQWmJQQERGRWWBSQkRERGbhoZKS0tJSY8VBREREFs6gpOT777/HiBEjEBgYCGtrazg4OMDJyQldu3bF3LlzcfXq1bqKk4iIiJ5weiUlW7ZsQdOmTZGQkIB69eph5syZ2Lx5M77//nusWrUKXbt2xc6dO9G4cWOMHz8ef/31V13HTURERE8YvSZPe/fdd/HBBx/g+eef15r2HQCGDBkCALhy5QoWL16M1atXIzEx0biREhER0RNNr6Tk7nVjavPUU0/h/ffff6iAiIjIslVUVkjP917Yi96BvWFVz6qWPehJwdE3RERkNjbnbEbw8mDpdb//6wflYiU252w2YVT0qBi89o0QAhs3bsSuXbuQn58vrdhbbfNm/uEQEZHhNudsxuANgyGgufrJlcIrGLxhMDYO2YjY4FgTRUePgsEtJVOmTEF8fDxyc3Mhl8uhUCg0HkRERIaqqKzAlIwpWgkJAKlsasZUjVs79OQxuKVk7dq12Lx5M/r161cX8RARkQXad3EfLhdernG7gMClwkvYd3Efuim7PbrA6JEyuKVEoVCgcePGdRELERFZqGu3rhm1Hj2eDE5KZs+ejTlz5qCkpKQu4iEiIgvk7eRt1Hr0eDL49s0LL7yAL774Ah4eHlAqlbC2ttbY/ssvvxgtOCIisgydG3WGr7MvrhRe0dmvRAYZfJ190blRZxNER4+KwUnJiBEjcOzYMbz44ovw9PSETCari7iIiMiCWNWzwuI+izF4w2CtbTJUfc+k9EnhfCVPOIOTkm3btuH7779Hp06d6iIeIiKyULHBsdg4ZCMmbZ+Eq7f+t5aar7MvUvqkcDiwBTA4KfHz84Ozs3NdxEJERBYuNjgWPQN6QjG/aoqJ74Z9xxldLYjBHV2Tk5Mxa9YsnD9/vg7CISIiS3d3AtLFvwsTEgticEvJiy++iOLiYgQGBsLBwUGro+uNGzeMFhwRERFZDoOTkpSUlDoIg4iIiCydwUlJQkJCXcRBREREFs7gpKRafn6+zgX5QkNDHzooIiIisjwGJyXHjh1DQkICcnJyIITmBDcymQwVFVwsiYiIiAxncFIycuRIBAUFYdWqVZw8jYiIiIzG4KQkNzcXmzdvRpMmTeoiHiIiIrJQBs9T0qNHD5w4caIuYiEiIiILZnBLyaeffoqEhARkZ2ejRYsWWvOUxMTEGC04IiIishwGJyUHDhzA/v37sX37dq1t7OhKRERED8rg2zeTJ09GfHw8rl27hsrKSo0HExIiIiJ6UAYnJdevX8e0adPg6elZF/EQERGRhTI4KYmNjcWuXbvqIhYiIiKyYAb3KQkKCkJSUhL279+Pli1banV0nTx5stGCIyIiIsshE/dOy3ofAQEBNR9MJsN//vOfhw7KHBUWFkKhUEClUsHZ2dnU4dBDKCoqglwuBwCo1Wo4OjqaOCIiultRWRHk8/77GU1Sw9GGn9HHnb7foQ80eRoRERGRsRncp4SIiIioLuiVlLz33nsoLi7W64CHDh3Ctm3bHiooIiIisjx6JSWnT59Go0aNMGHCBGzfvh1//fWXtO3OnTs4efIkli9fjoiICAwdOpR9LoiIiMhgevUpWb16NU6ePIkPP/wQw4cPh0qlgpWVFWxtbaUWlNatW2Ps2LFISEiAra1tnQZNRERETx69O7qGhoZi5cqVWLFiBU6ePInz58+jpKQE7u7uaNWqFdzd3esyTiIiInrCGTz6RiaTISwsDGFhYXURDxEREVkojr4hIiIis8CkhIiIiMwCkxIiIiIyC0xKiIiIyCwwKSEiIiKzYPDom6KiIrz33nvIzMxEfn4+KisrNbY/qQvyERERUd0yOCkZPXo09uzZg/j4eHh7e0Mmk9VFXERERGRhDE5Ktm/fjm3btqFjx451EQ8RERFZKIP7lLi6usLNza0uYiEiIiILZnBS8s477+DNN9/Ue9VgIiIiIn0YfPsmOTkZ586dg6enJ5RKJaytrTW2//LLL0YLjoiIiCyHwUnJgAED6iAMokenoqJCer5371707t0bVlZWJoyIiIiAB0hK3nrrrbqIg+iR2Lx5MyZNmiS97tevH3x9fbF48WLExsaaMDIiIjI4Kal27Ngx5OTkQCaTISQkBK1btzZmXERGt3nzZgwePBhCCI3yK1euYPDgwdi4cSMTEyIiEzI4KcnPz8fQoUOxe/duuLi4QAgBlUqF7t27Y/369WjYsGFdxEn0UCoqKjBlyhSthAQAhBCQyWSYOnUq+vfvz1s5REQmYvDom0mTJqGwsBCnTp3CjRs3UFBQgOzsbBQWFmLy5MkGHWvv3r2Ijo6Gj48PZDIZtmzZorFdJpPpfCxYsECqU1paikmTJsHd3R2Ojo6IiYnB5cuXNY5TUFCA+Ph4KBQKKBQKxMfH4+bNm4ZeOj3G9u3bp/V3cTchBC5duoR9+/Y9wqiIiOhuBiclGRkZ+OijjxAcHCyVhYSE4MMPP8T27dsNOlZRURHCwsKwbNkynduvXbum8fjss88gk8kwaNAgqc7UqVORnp6O9evXY//+/VCr1YiKitLozDhs2DBkZWUhIyMDGRkZyMrKQnx8vIFXTo+za9euGbUeEREZn8G3byorK7WGAQOAtbW11jo499O3b1/07du3xu1eXl4ar7/++mt0794djRs3BgCoVCqsWrUKa9asQc+ePQEAa9euhZ+fH3bu3InIyEjk5OQgIyMDP//8M9q3bw8A+OSTT9ChQwecOXMGTZs2NShmejx5e3sbtR4RERmfwS0lzz33HKZMmYKrV69KZVeuXMG0adPQo0cPowZ3tz///BPbtm3Dyy+/LJUdO3YM5eXl6N27t1Tm4+ODFi1a4MCBAwCAgwcPQqFQSAkJAISHh0OhUEh1dCktLUVhYaHGgx5fnTt3hq+vb41rNclkMvj5+aFz586PODIiIqpmcFKybNky3Lp1C0qlEoGBgWjSpAkCAgJw69YtLF26tC5iBAB8/vnncHJy0hgdkZeXBxsbG7i6umrU9fT0RF5enlTHw8ND63geHh5SHV3mzZsn9UFRKBTw8/Mz0pWQKVhZWWHx4sU6t1UnKikpKezkSkRkQgbfvvHz88Mvv/yCHTt24LfffoMQAiEhIdLtk7ry2WefYfjw4bCzs7tv3erRFNV0/Tq+t869kpKSMH36dOl1YWEhE5PHXGxsLDZu3IhJkyZptPT5+voiJSWFw4GJiEzsgecp6dWrF3r16mXMWGq0b98+nDlzBl9++aVGuZeXF8rKylBQUKDRWpKfn4+IiAipzp9//ql1zL/++guenp41ntPW1ha2trZGugIyF7GxsejZsycUCgUA4LvvvuOMrkREZkKvpGTJkiUYO3Ys7OzssGTJklrrGjosWB+rVq1CmzZtEBYWplHepk0bWFtbY8eOHRgyZAiAqtET2dnZeP/99wEAHTp0gEqlwuHDh9GuXTsAwKFDh6BSqaTEhSzL3QlIly5dmJAQEZkJvZKSRYsWSbdOFi1aVGM9mUxmUFKiVqtx9uxZ6XVubi6ysrLg5uaGRo0aAai6bfLVV18hOTlZa3+FQoGXX34ZiYmJaNCgAdzc3DBjxgy0bNlSup0UHByMPn36YMyYMVi5ciUAYOzYsYiKiuLIGyIiIjOiV1KSm5ur8/nDOnr0KLp37y69ru7DkZCQgLS0NADA+vXrIYTAP/7xD53HWLRoEerXr48hQ4agpKQEPXr0QFpamsav33Xr1mHy5MnSKJ2YmJga50YhIiIi05AJXfNu1+Ltt9/GjBkz4ODgoFFeUlKCBQsW4M033zRqgOaisLAQCoUCKpUKzs7Opg6HHkJRURHkcjmAqtY6R0dHE0dERHcrKiuCfN5/P6NJajja8DP6uNP3O9TgIcFz5syBWq3WKi8uLsacOXMMPRwRERERgAdISmoaSnvixAm4ubkZJSgiIiKyPHoPCXZ1dZUWxAsKCtJITCoqKqBWqzF+/Pg6CZKIiIiefHonJSkpKRBCYNSoUZgzZ440zwMA2NjYQKlUokOHDnUSJBERET359E5KEhISAAABAQGIiIjQuSgfERER0YMyeEbXrl27Ss9LSkpQXl6usZ0jU4iIiOhBGNzRtbi4GK+88go8PDwgl8vh6uqq8SAiIiJ6EAYnJTNnzsSPP/6I5cuXw9bWFp9++inmzJkDHx8frF69ui5iJCIiIgtg8O2bb775BqtXr0a3bt0watQodO7cGU2aNIG/vz/WrVuH4cOH10WcRERE9IQzuKXkxo0bCAgIAFDVf+TGjRsAgE6dOmHv3r3GjY6IiIgshsFJSePGjXH+/HkAQEhICDZs2ACgqgXFxcXFmLERERGRBTE4KRk5ciROnDgBAEhKSpL6lkybNg0zZ840eoBERERkGQzuUzJt2jTpeffu3fHbb7/h6NGjCAwMRFhYmFGDIyIiIsthcEvJ6tWrUVpaKr1u1KgRYmNjERwczNE3RET00BysHaBOUkOdpIaDtcP9d6AnxgPdvlGpVFrlt27dwsiRI40SFBERWS6ZTAZHG0c42jjqXACWnlxGWyX48uXLGuvhEBERERlC7z4lrVu3llYJ7tGjB+rX/9+uFRUVyM3NRZ8+feokSCIiInry6Z2UDBgwAACQlZWFyMhIyOVyaVv1KsGDBg0yeoBERERkGfROSt566y0AgFKpRFxcHOzs7OosKCIiIrI8Bg8JTkhIqIs4iIiIyMLplZS4ubnh999/h7u7O1xdXWvtDV097TwRERGRIfRKShYtWgQnJycAQEpKSl3GQ0RkGCGA4uKq5w4OAIeQEj22ZEIIYeogHgeFhYVQKBRQqVRwdnY2dTj0EIqKiqSO2mq1Go6OjiaOiB5KURFQ3fFerQb470lkdvT9DjW4TwlQNQQ4PT0dOTk5kMlkCA4ORv/+/TWGCRMREREZwuAsIjs7G/3790deXh6aNm0KAPj999/RsGFDbN26FS1btjR6kERERPTkM3hG19GjR6N58+a4fPkyfvnlF/zyyy+4dOkSQkNDMXbs2LqIkYiIiCyAwS0lJ06cwNGjR+Hq6iqVubq6Yu7cuWjbtq1RgyMiIiLLYXBLSdOmTfHnn39qlefn56NJkyZGCYqIiIgsj8FJybvvvovJkydj48aNuHz5Mi5fvoyNGzdi6tSpmD9/PgoLC6UHERERkb4MHhJcr97/8pjqSdSqD3H3a5lMhoqKCmPFaXIcEvzk4JDgJwyHBBOZvTobErxr166HCoyIiIhIF4OTkq5du9ZFHERERGTh9EpKTp48iRYtWqBevXo4efJkrXVDQ0ONEhgRERFZFr2SklatWiEvLw8eHh5o1aoVZDIZdHVFedL6kRAREdGjo1dSkpubi4YNG0rPiYiIiIxNr6TE399f53MiIiIiYzF4npJ58+bhs88+0yr/7LPPMH/+fKMERURERJbH4KRk5cqVaNasmVZ58+bNsWLFCqMERURERJbH4KQkLy8P3t7eWuUNGzbEtWvXjBIUERERWR6DkxI/Pz/89NNPWuU//fQTfHx8jBIUERERWR6DJ08bPXo0pk6divLycjz33HMAgMzMTMyaNQuJiYlGD5CIiIgsg8FJyaxZs3Djxg1MnDgRZWVlAAA7Ozu8+uqrSEpKMnqAREREZBkMXpCvmlqtRk5ODuzt7fH000/D1tbW2LGZFS7I9+TggnxPGC7IR2T29P0ONbhPSbW8vDzcuHEDgYGBsLW11TnDKxEREZG+DE5Krl+/jh49eiAoKAj9+vWTRtyMHj2afUqIiIjogRmclEybNg3W1ta4ePEiHBwcpPK4uDhkZGQYNTgiIiKyHAZ3dP3hhx/w/fffw9fXV6P86aefxoULF4wWGBEREVkWg1tKioqKNFpIqv39999PfGdXIiIiqjsGJyVdunTB6tWrpdcymQyVlZVYsGABunfvbtTgiIiIyHIYfPtmwYIF6NatG44ePYqysjLMmjULp06dwo0bN3TO9EpERESkD4NbSkJCQnDy5Em0a9cOvXr1QlFREWJjY3H8+HEEBgbWRYxERERkAQxqKSkvL0fv3r2xcuVKzJkzp65iIiIiIgtkUEuJtbU1srOzIZPJ6ioeIiIislAG37556aWXsGrVqrqIhYiIiCyYwR1dy8rK8Omnn2LHjh149tlntdYNWbhwodGCIyIiIsthcFKSnZ2NZ555BgDw+++/a2zjbR0iIiJ6UAYnJbt27aqLOIgeGQcHB6jVauk5ERGZhwdeJRgALl26hMuXLz/w/nv37kV0dDR8fHwgk8mwZcsWrTo5OTmIiYmBQqGAk5MTwsPDcfHiRWl7aWkpJk2aBHd3dzg6OiImJkYrpoKCAsTHx0OhUEChUCA+Ph43b9584Ljp8SaTyeDo6AhHR0e27hERmRGDk5I7d+7gjTfegEKhgFKphL+/PxQKBV5//XWUl5cbdKyioiKEhYVh2bJlOrefO3cOnTp1QrNmzbB7926cOHECb7zxBuzs7KQ6U6dORXp6OtavX4/9+/dDrVYjKioKFRUVUp1hw4YhKysLGRkZyMjIQFZWFuLj4w29dCIiIqpLwkDjxo0THh4eYsWKFeLEiRPixIkTYsWKFcLLy0uMGzfO0MNJAIj09HSNsri4OPHiiy/WuM/NmzeFtbW1WL9+vVR25coVUa9ePZGRkSGEEOL06dMCgPj555+lOgcPHhQAxG+//aZ3fCqVSgAQKpVK732I6BFQq4UAqh5qtamjISId9P0ONbil5IsvvkBaWhrGjRuH0NBQhIaGYty4cfjss8/wxRdfGC1ZqqysxLZt2xAUFITIyEh4eHigffv2Grd4jh07Jk3oVs3HxwctWrTAgQMHAAAHDx6EQqFA+/btpTrh4eFQKBRSHV1KS0tRWFio8SAiIqK6Y3BSYmdnB6VSqVWuVCphY2NjjJgAAPn5+VCr1XjvvffQp08f/PDDDxg4cCBiY2OxZ88eAEBeXh5sbGzg6uqqsa+npyfy8vKkOh4eHlrH9/DwkOroMm/ePKkPikKhgJ+fn9GujYiIiLQZnJT885//xDvvvIPS0lKprLS0FHPnzsUrr7xitMAqKysBAP3798e0adPQqlUr/Otf/0JUVBRWrFhR675CCI0OjLo6M95b515JSUlQqVTS49KlSw94JURERKQPg4cEHz9+HJmZmfD19UVYWBgA4MSJEygrK0OPHj0QGxsr1d28efMDB+bu7o769esjJCREozw4OBj79+8HAHh5eaGsrAwFBQUarSX5+fmIiIiQ6vz5559ax//rr7/g6elZ4/ltbW1ha2v7wPETERGRYQxOSlxcXDBo0CCNsrq4tWFjY4O2bdvizJkzGuW///47/P39AQBt2rSBtbU1duzYgSFDhgAArl27huzsbLz//vsAgA4dOkClUuHw4cNo164dAODQoUNQqVRS4kJERESmZ3BSkpqaarSTq9VqnD17Vnqdm5uLrKwsuLm5oVGjRpg5cybi4uLQpUsXdO/eHRkZGfjmm2+we/duAIBCocDLL7+MxMRENGjQAG5ubpgxYwZatmyJnj17AqhqWenTpw/GjBmDlStXAgDGjh2LqKgoNG3a1GjXQkRERA/pkYwFqsGuXbsEAK1HQkKCVGfVqlWiSZMmws7OToSFhYktW7ZoHKOkpES88sorws3NTdjb24uoqChx8eJFjTrXr18Xw4cPF05OTsLJyUkMHz5cFBQUGBQrhwQTmSkOCSYye/p+h8qEEMKEOdFjo7CwEAqFAiqVCs7OzqYOh4iqFRUBcnnVc7UauGeRUCIyPX2/Qx9qmnkiIiIiY2FSQkRERGaBSQkRERGZBb1G3yxZskTvA06ePPmBgyEiIiLLpVdH14CAAP0OJpPhP//5z0MHZY7Y0ZXITLGjK5HZ0/c7VK+WktzcXKMFRkRERKQL+5QQERGRWTB4RlcAuHz5MrZu3YqLFy+irKxMY9vChQuNEhgRERFZFoOTkszMTMTExCAgIABnzpxBixYtcP78eQgh8Mwzz9RFjERERGQBDL59k5SUhMTERGRnZ8POzg6bNm3CpUuX0LVrV7zwwgt1ESMRERFZAIOTkpycHCQkJAAA6tevj5KSEsjlcrz99tuYP3++0QMkIiIiy2BwUuLo6IjS0lIAgI+PD86dOydt+/vvv40XGREREVkUg/uUhIeH46effkJISAief/55JCYm4tdff8XmzZsRHh5eFzESERGRBTA4KVm4cCHUajUAYPbs2VCr1fjyyy/RpEkTLFq0yOgBEhERkWXQa0ZX4oyuRGaLM7oSmT19v0MN7lPSuHFjXL9+Xav85s2baNy4saGHIyIiIgLwAEnJ+fPnUVFRoVVeWlqKK1euGCUoIiIisjx69ynZunWr9Pz777+HQqGQXldUVCAzMxNKpdKowRER3dfdP5L27gV69wasrEwXDxE9ML37lNSrV9WoIpPJcO8u1tbWUCqVSE5ORlRUlPGjNAPsU0JkhjZvBiZNAq5e/V+Zry+weDEQG2u6uIhIg1FXCQaAyspKAEBAQACOHDkCd3f3h4+SiOhBbd4MDB4M3Pu76sqVqvKNG5mYED1mOPpGT2wpITIjFRWAUglcvqx7u0xW1WKSm8tbOURmoM5G3wDAnj17EB0djSZNmuDpp59GTEwM9u3b98DBEhEZZN++mhMSoKr15NKlqnpE9NgwOClZu3YtevbsCQcHB0yePBmvvPIK7O3t0aNHD/zf//1fXcRIRKTp2jXj1iMis2Dw7Zvg4GCMHTsW06ZN0yhfuHAhPvnkE+Tk5Bg1QHPB2zdEZmT3bqB79/vX27UL6NatrqMhovuos9s3//nPfxAdHa1VHhMTg9zcXEMPR0RkuM6dq/qMyGS6t8tkgJ9fVT0iemwYnJT4+fkhMzNTqzwzMxN+fn5GCYqIqFZWVlXDfnWpTlRSUtjJlegxo/eQ4FGjRmHx4sVITEzE5MmTkZWVhYiICMhkMuzfvx9paWlYXNN/EkRExhYbWzXsV9c8JSkpHA5M9BjSu0+JlZUVrl27Bg8PD6SnpyM5OVnqPxIcHIyZM2eif//+dRqsKbFPCZGZKiwEqmeY/u47zuhKZIb0/Q41aEbXvLw8eHh4GC3IxwmTEiIzxVWCicye0Wd0BaqmmKfaVVRUoLy83NRhEJkta2trWLElg4h0MCgpCQoKum9icuPGjYcK6HElhEBeXh5u3rxp6lCIzJ6Liwu8vLz4Q4eINBiUlMyZM0djdWD6n+qExMPDAw4ODvzPlkgHIQSKi4uRn58PAPD29jZxRERkTgxKSoYOHWqxfUpqU1FRISUkDRo0MHU4RGbN3t4eAJCfnw8PDw/eyiEiid7zlPCXf82q+5A4ODiYOBKix0P1Z4X9r4jobnonJVxM+P6YuBHph58VItJF79s3lZWVdRkHERERWTiDp5knelBpaWlwcXExdRhERGSmmJRYuPz8fIwbNw6NGjWCra0tvLy8EBkZiYMHDwKoambfsmWLwcdVKpVISUnRKIuLi8Pvv/9uhKiJiOhJZNDoG3ryDBo0COXl5fj888/RuHFj/Pnnn8jMzKyT+Wbs7e2lkRemVFZWBhsbG1OHQURE92BLiQW7efMm9u/fj/nz56N79+7w9/dHu3btkJSUhOeffx5KpRIAMHDgQMhkMun1uXPn0L9/f3h6ekIul6Nt27bYuXOndNxu3brhwoULmDZtGmQymdSpUdftm48++giBgYGwsbFB06ZNsWbNGo3tMpkMn376KQYOHAgHBwc8/fTT2Lp1q7S9oqICL7/8MgICAmBvb4+mTZtqLQw5YsQIDBgwAPPmzYOPjw+CgoLw9ttvo2XLllrvSZs2bfDmm28+6FtKREQPgS0ldaB6gihTMGTiNrlcDrlcji1btiA8PBy2trYa248cOQIPDw+kpqaiT58+0nwSarUa/fr1w//7f/8PdnZ2+PzzzxEdHY0zZ86gUaNG2Lx5M8LCwjB27FiMGTOmxvOnp6djypQpSElJQc+ePfHtt99i5MiR8PX1Rffu3aV6c+bMwfvvv48FCxZg6dKlGD58OC5cuAA3NzdUVlbC19cXGzZsgLu7Ow4cOICxY8fC29sbQ4YMkY6RmZkJZ2dn7NixA0IIuLi4YM6cOThy5Ajatm0LADh58iSOHz+Or776Su/3m4iIjEiQXlQqlQAgVCqV1raSkhJx+vRpUVJSIoQQQq1WCwAmeajVaoOua+PGjcLV1VXY2dmJiIgIkZSUJE6cOCFtByDS09Pve5yQkBCxdOlS6bW/v79YtGiRRp3U1FShUCik1xEREWLMmDEadV544QXRr18/jfO//vrr0mu1Wi1kMpnYvn17jbFMnDhRDBo0SHqdkJAgPD09RWlpqUa9vn37igkTJkivp06dKrp161b7hZJR3PuZeShqtRBA1cPAv38iejRq+w69G2/fWLhBgwbh6tWr2Lp1KyIjI7F7924888wzSEtLq3GfoqIizJo1CyEhIXBxcYFcLsdvv/2GixcvGnTunJwcdOzYUaOsY8eOyMnJ0SgLDQ2Vnjs6OsLJyUmaphwAVqxYgWeffRYNGzaEXC7HJ598ohVLy5YttfqRjBkzBl988QVu376N8vJyrFu3DqNGjTLoGoiIyHh4+6YOODg4QK1Wm+zchrKzs0OvXr3Qq1cvvPnmmxg9ejTeeustjBgxQmf9mTNn4vvvv8cHH3yAJk2awN7eHoMHD0ZZWZnB5773VpMQQqvM2tpaa5/qeXM2bNiAadOmITk5GR06dICTkxMWLFiAQ4cOaezjqGM5++joaNja2iI9PR22trYoLS3FoEGDDL4GIiIyDiYldUAmk+n8EnxchISESMOAra2tUVFRobF93759GDFiBAYOHAigqo/J+fPnNerY2Nho7Xev4OBg7N+/Hy+99JJUduDAAQQHB+sd6759+xAREYGJEydKZefOndNr3/r16yMhIQGpqamwtbXF0KFDuVQAEZEJMSmxYNevX8cLL7yAUaNGITQ0FE5OTjh69Cjef/999O/fH0DVfCOZmZno2LEjbG1t4erqiiZNmmDz5s2Ijo6GTCbDG2+8oTXjr1KpxN69ezF06FDY2trC3d1d6/wzZ87EkCFD8Mwzz6BHjx745ptvsHnzZo2RPPfTpEkTrF69Gt9//z0CAgKwZs0aHDlyBAEBAXrtP3r0aCkJ+umnn/Q+LxERGR/7lFgwuVyO9u3bY9GiRejSpQtatGiBN954A2PGjMGyZcsAAMnJydixYwf8/PzQunVrAMCiRYvg6uqKiIgIREdHIzIyEs8884zGsd9++22cP38egYGBaNiwoc7zDxgwAIsXL8aCBQvQvHlzrFy5EqmpqejWrZve1zB+/HjExsYiLi4O7du3x/Xr1zVaTe7n6aefRkREBJo2bYr27dvrvR8RERmfTAiutKePwsJCKBQKqFQqODs7a2y7ffs2cnNzERAQADs7OxNFSA9CCIFmzZph3LhxmD59uqnDsRhG/cwUFQFyedVztRp4jG+dEj2pavsOvRtv35DFys/Px5o1a3DlyhWMHDnS1OEQEVk8JiVksTw9PeHu7o6PP/4Yrq6upg6HiMjiMSkhi8U7l0RE5oUdXYmIiMgsMCkhIiIis8CkhIiIiMwCkxIiIiIyC0xKiIiIyCwwKSEiIiKzYNKkZO/evYiOjoaPjw9kMpm0CFy1ESNGQCaTaTzCw8M16pSWlmLSpElwd3eHo6MjYmJicPnyZY06BQUFiI+Ph0KhgEKhQHx8PG7evFnHV0fmSKlUIiUlxdRhPJTr16/Dw8NDaxFEczF48GAsXLjQ1GEQ0WPIpElJUVERwsLCpHVWdOnTpw+uXbsmPb777juN7VOnTkV6ejrWr1+P/fv3Q61WIyoqSmOF2mHDhiErKwsZGRnIyMhAVlYW4uPj6+y6HhfVSd97772nUb5lyxbIZLJHEsOmTZvQvn17KBQKODk5oXnz5khMTJS2z549G61atTL4uGlpaXBxcdEqP3LkCMaOHfsQEZvevHnzEB0dDaVSKZVNmTIFbdq0ga2tbY3v14YNG9CqVSs4ODjA398fCxYs0Ni+e/durR8BMpkMv/32m1QnLS1NZ53bt29Ldd58803MnTsXhYWFRr1uInrymXTytL59+6Jv37611rG1tYWXl5fObSqVCqtWrcKaNWvQs2dPAMDatWvh5+eHnTt3IjIyEjk5OcjIyMDPP/8sLbj2ySefoEOHDjhz5gyaNm1q3It6zNjZ2WH+/PkYN27cI5/VdOfOnRg6dCjeffddxMTEQCaT4fTp08jMzKyzc9a0OOCjVl5eDmtra4P3KykpwapVq7SScyEERo0ahUOHDuHkyZNa+23fvh3Dhw/H0qVL0bt3b+Tk5GD06NGwt7fHK6+8olH3zJkzGmtT3PueOTs748yZMxpld69fExoaCqVSiXXr1mHChAkGXyMRWS6z71Oye/dueHh4ICgoCGPGjEF+fr607dixYygvL0fv3r2lMh8fH7Ro0QIHDhwAABw8eBAKhUJjBdjw8HAoFAqpji6lpaUoLCzUeDyJevbsCS8vL8ybN6/Weps2bULz5s1ha2sLpVKJ5ORkje1KpRLvvvsuRo0aBScnJzRq1Agff/xxrcf89ttv0alTJ8ycORNNmzZFUFAQBgwYgKVLlwKo+lU+Z84cnDhxQvpFnpaWBgBYuHAhWrZsCUdHR/j5+WHixIlQq9UAqv5mRo4cCZVKJe03e/ZsKc67b99cvHgR/fv3h1wuh7OzM4YMGYI///xT2l7dUrNmzRoolUooFAoMHToUt27dkupkZGSgU6dOcHFxQYMGDRAVFYVz585J28+fPw+ZTIYNGzagW7dusLOzw8cffwxnZ2ds3LhR4z355ptv4OjoqHH8u23fvh3169dHhw4dNMqXLFmCf/7zn2jcuLHO/dasWYMBAwZg/PjxaNy4MZ5//nm8+uqrmD9/vtbMth4eHvDy8pIeVlZWGttlMpnGdl0/GmJiYvDFF1/ojIWIqCZmnZT07dsX69atw48//ojk5GQcOXIEzz33HEpLSwEAeXl5sLGx0fqF7+npiby8PKmOh4eH1rE9PDykOrrMmzdP6oOiUCjg5+enf+BCVK1caoqHgVOnW1lZ4d1338XSpUu1+uJUO3bsGIYMGYKhQ4fi119/xezZs/HGG29ICUK15ORkPPvsszh+/DgmTpyICRMmaDT938vLywunTp1Cdna2zu1xcXFITExE8+bNpdt3cXFxAIB69ephyZIlyM7Oxueff44ff/wRs2bNAgBEREQgJSUFzs7O0n4zZszQOr4QAgMGDMCNGzewZ88e7NixA+fOnZPOUe3cuXPYsmULvv32W3z77bfYs2ePxi2voqIiTJ8+HUeOHEFmZibq1auHgQMHorKyUuM4r776KiZPnoycnBwMHDgQQ4cORWpqqkad1NRUDB48GE5OTjrfk7179+LZZ5+t8T2tSWlpqdZqvPb29rh8+TIuXLigUd66dWt4e3ujR48e2LVrl9ax1Go1/P394evri6ioKBw/flyrTrt27XD48GHps0pEpBdhJgCI9PT0WutcvXpVWFtbi02bNgkhhFi3bp2wsbHRqtezZ08xbtw4IYQQc+fOFUFBQVp1mjRpIubNm1fjuW7fvi1UKpX0uHTpkgAgVCqVVt2SkhJx+vRpUVJSUlWgVgtRlR48+odaXet7eLeEhATRv39/IYQQ4eHhYtSoUUIIIdLT08XdfxrDhg0TvXr10th35syZIiQkRHrt7+8vXnzxRel1ZWWl8PDwEB999FGN51er1aJfv34CgPD39xdxcXFi1apV4vbt21Kdt956S4SFhd33WjZs2CAaNGggvU5NTRUKhUKrnr+/v1i0aJEQQogffvhBWFlZiYsXL0rbT506JQCIw4cPS+d3cHAQhYWFGtfevn37GmPJz88XAMSvv/4qhBAiNzdXABApKSka9Q4dOiSsrKzElStXhBBC/PXXX8La2lrs3r27xmP3799f+nfSpab3a+XKlcLBwUHs3LlTVFRUiDNnzohmzZoJAOLAgQNCCCF+++038fHHH4tjx46JAwcOiAkTJgiZTCb27NkjHefgwYNizZo1IisrS+zdu1cMGjRI2Nvbi99//13jfCdOnBAAxPnz53XGqfWZeRh3f94M+PsnokdHpVLV+B16N7NuKbmXt7c3/P398ccffwCo+qVdVlaGgoICjXr5+fnw9PSU6tzdHF/tr7/+kuroYmtrC2dnZ43Hk2z+/Pn4/PPPcfr0aa1tOTk56Nixo0ZZx44d8ccff2h0KA4NDZWeVzfxV99u69u3L+RyOeRyOZo3bw4AcHR0xLZt23D27Fm8/vrrkMvlSExMRLt27VBcXFxrvLt27UKvXr3w1FNPwcnJCS+99BKuX7+OoqIiva85JycHfn5+Gq1gISEhcHFxQU5OjlSmVCo1Wi68vb01biOeO3cOw4YNQ+PGjeHs7IyAgAAAVbeG7nZvC0e7du3QvHlzrF69GkDVLZZGjRqhS5cuNcZcUlKi1eKhjzFjxuCVV15BVFQUbGxsEB4ejqFDhwKAdHumadOmGDNmDJ555hl06NABy5cvx/PPP48PPvhAOk54eDhefPFFhIWFoXPnztiwYQOCgoKkW27V7O3tAeC+/45ERHd7rJKS69ev49KlS/D29gYAtGnTBtbW1tixY4dU59q1a8jOzkZERAQAoEOHDlCpVDh8+LBU59ChQ1CpVFIdo3NwANRq0zwcHB4o5C5duiAyMhKvvfaa1jYhhNZoHKHjNtG9HTdlMpl0C+PTTz9FVlYWsrKytDppBgYGYvTo0fj000/xyy+/4PTp0/jyyy9rjPXChQvo168fWrRogU2bNuHYsWP48MMPAVR1INWXruvSVV7bdQFAdHQ0rl+/jk8++QSHDh3CoUOHAABlZWUa+zk6Omqda/To0dItnNTUVIwcObLWkU/u7u5aSbg+ZDIZ5s+fD7VajQsXLiAvLw/t2rUDAI1RPPcKDw+XfgToUq9ePbRt21arzo0bNwCYT8diIno8mHT0jVqtxtmzZ6XXubm5yMrKgpubG9zc3DB79mwMGjQI3t7eOH/+PF577TW4u7tj4MCBAACFQoGXX34ZiYmJaNCgAdzc3DBjxgy0bNlSGo0THByMPn36YMyYMVi5ciUAYOzYsYiKiqq7kTcyGaDjC8jcvffee2jVqhWCgoI0ykNCQrB//36NsgMHDiAoKEirE2RNnnrqKb3qKZVKODg4SC0eNjY2Gq0xAHD06FHcuXMHycnJqFevKq/esGGDRh1d+90rJCQEFy9exKVLl6TWktOnT0OlUiE4OFiveK9fv46cnBysXLkSnTt3BgCt96o2L774ImbNmoUlS5bg1KlTSEhIqLV+69atsXbtWr2Pfy8rKyvp3+KLL75Ahw4ddPa5qnb8+HHpR4AuQghkZWWhZcuWGuXZ2dnw9fWFu7v7A8dKRJbHpEnJ0aNH0b17d+n19OnTAQAJCQn46KOP8Ouvv2L16tW4efMmvL290b17d3z55ZcaTemLFi1C/fr1MWTIEJSUlKBHjx5IS0vT+LJct24dJk+eLI3SiYmJqXVuFEvVsmVLadjo3RITE9G2bVu88847iIuLw8GDB7Fs2TIsX778oc43e/ZsFBcXo1+/fvD398fNmzexZMkSlJeXo1evXgCqkpTqZNXX1xdOTk4IDAzEnTt3sHTpUkRHR+Onn37CihUrNI6tVCqhVquRmZmJsLAwODg4wOGeVqSePXsiNDQUw4cPR0pKCu7cuYOJEyeia9euencmdXV1RYMGDfDxxx/D29sbFy9exL/+9S+93wNXV1fExsZi5syZ6N27N3x9fWutHxkZiaSkJBQUFGh08D579izUajXy8vJQUlKCrKwsAFWJl42NDf7++29s3LgR3bp1w+3bt5GamoqvvvoKe/bskY6RkpICpVKJ5s2bo6ysDGvXrsWmTZuwadMmqc6cOXMQHh6Op59+GoWFhViyZAmysrKklqpq+/bt0xgVR0Skl7rv3vJkqK2TjlE77T1Cd3d0rXb+/Hlha2sr7v3T2LhxowgJCRHW1taiUaNGYsGCBRrb7+5AWi0sLEy89dZbNZ7/xx9/FIMGDRJ+fn7CxsZGeHp6ij59+oh9+/ZJdW7fvi0GDRokXFxcBACRmpoqhBBi4cKFwtvbW9jb24vIyEixevVqAUAUFBRI+44fP140aNBAAJDiuDfOCxcuiJiYGOHo6CicnJzECy+8IPLy8qTtujqOLlq0SPj7+0uvd+zYIYKDg4Wtra0IDQ0Vu3fv1ui4Xd3R9fjx4zrfh8zMTAFAbNiwocb36m7h4eFixYoVGmVdu3YVALQeubm5QoiqTrTh4eHC0dFRODg4iB49eoiff/5Z4xjz588XgYGBws7OTri6uopOnTqJbdu2adSZOnWqaNSokbCxsRENGzYUvXv3ljrKVispKRHOzs7i4MGDNV4DO7oSWRZ9O7rKhDBwDKmFKiwshEKhgEql0ur0evv2beTm5iIgIOCBOiGSZVu3bh2mTJmCq1evwsbG5r71v/vuO8yYMQPZ2dnS7Stz8uGHH+Lrr7/GDz/8UGMdo35miooAubzquVr9WN46JXrS1fYdejeT3r4hsmTFxcXIzc3FvHnzMG7cOL0SEgDo168f/vjjD1y5csWw+XMeEWtra61bgERE+jC/n1lEFuL9999Hq1at4OnpiaSkJIP2nTJlilkmJEBVR3JLX76BiB4MkxIiE5k9ezbKy8uRmZkJefXtByIiC8akhIiIiMwCkxIiIiIyC0xKiIiIyCwwKSEiIiKzwKSEiIiIzAKTEiIiIjILTErIoiiVSqSkpJg6jIdy/fp1eHh44Pz586YORafBgwdj4cKFpg6DiB5DTEos2IgRIyCTyfDee+9plG/ZsgUymeyRxLBp0ya0b98eCoUCTk5OaN68ORITE6Xts2fPRqtWrQw+blpaGlxcXLTKjxw5grFjxz5ExKY3b948REdHQ6lUAgBOnDiBf/zjH/Dz84O9vT2Cg4OxePFijX1u376NESNGoGXLlqhfvz4GDBig89h79uxBmzZtYGdnh8aNG2stdAhU/ZuFhITA1tYWISEhSE9P19j+5ptvYu7cuSgsLDTK9RKR5WBSYuHs7Owwf/58FBQUPPJz79y5E0OHDsXgwYNx+PBhHDt2DHPnzkVZWVmdnbNhw4ZaqwWbQnl5+QPtV1JSglWrVmH06NFS2bFjx9CwYUOsXbsWp06dwr///W8kJSVprIRdUVEBe3t7TJ48GT179tR57NzcXPTr1w+dO3fG8ePH8dprr2Hy5MkaqwQfPHgQcXFxiI+Px4kTJxAfH48hQ4bg0KFDUp3Q0FAolUqsW7fuga6RiCzYI1ke8AnwpK4SHBUVJZo1ayZmzpwplaenp9e4SrCNjY3w9/cXH3zwgcZ2f39/MXfuXDFy5Eghl8uFn5+fWLlyZa3nnzJliujWrVuN21NTU7VWva1eJTg5OVm0aNFCODg4CF9fXzFhwgRx69YtIYQQu3bt0trvYVcJXr16tfD39xfOzs4iLi5OFBYWSnW2b98uOnbsKBQKhXBzcxPPP/+8OHv2rLS9epXgL7/8UnTt2lXY2tqKZcuWCScnJ/HVV19pXPPWrVuFg4ODxvHvtmnTJuHu7l7r+yqEEBMnThTdu3fXuU3X6tBCCDFr1izRrFkzjbJx48aJ8PBw6fWQIUNEnz59NOpERkaKoUOHapTNnj1bdO7cucb4uEowkWXRd5VgtpTUASEEisqKTPIQBi76bGVlhXfffRdLly7F5cuXddY5duwYhgwZgqFDh+LXX3/F7Nmz8cYbbyAtLU2jXnJyMp599lkcP34cEydOxIQJE/Dbb7/VeG4vLy+cOnUK2dnZOrfHxcUhMTERzZs3x7Vr13Dt2jXExcUBAOrVq4clS5YgOzsbn3/+OX788UfMmjULABAREYGUlBQ4OztL+82YMUPr+EIIDBgwADdu3MCePXuwY8cOnDt3TjpHtXPnzmHLli349ttv8e2332LPnj0at7yKioowffp0HDlyBJmZmahXrx4GDhyIyspKjeO8+uqrmDx5MnJycjBw4EAMHToUqampGnVSU1MxePBgODk56XxP9u7di2effbbG97SaSqWCm5vbfevd7eDBg+jdu7dGWWRkJI4ePSq17NRU58CBAxpl7dq1w+HDh1FaWmpQDERk2bhKcB0oLi+GfJ5p1jJRJ6nhaGPY0u0DBw5Eq1at8NZbb2HVqlVa2xcuXIgePXrgjTfeAAAEBQXh9OnTWLBgAUaMGCHV69evHyZOnAig6gt40aJF2L17N5o1a6bzvJMmTcK+ffvQsmVL+Pv7Izw8HL1798bw4cNha2sLe3t7yOVy1K9fH15eXhr7Tp06VXoeEBCAd955BxMmTMDy5cthY2MDhUIBmUymtd/ddu7ciZMnTyI3N1da3G7NmjVo3rw5jhw5grZt2wIAKisrkZaWJiUK8fHxyMzMxNy5cwEAgwYN0jjuqlWr4OHhgdOnT6NFixYaMcfGxkqvR48ejYiICFy9ehU+Pj74+++/8e2332LHjh01xnz+/Hn4+PjUuB2oShw2bNiAbdu21VrvXnl5efD09NQo8/T0xJ07d/D333/D29u7xjp5eXkaZU899RRKS0uRl5cHf39/g+IgIsvFlhICAMyfPx+ff/45Tp8+rbUtJycHHTt21Cjr2LEj/vjjD1RUVEhloaGh0vPqhCA/Px8A0LdvX8jlcsjlcjRv3hwA4OjoiG3btuHs2bN4/fXXIZfLkZiYiHbt2qG4uLjWeHft2oVevXrhqaeegpOTE1566SVcv34dRUVFel9zTk4O/Pz8NFbbDQkJgYuLC3JycqQypVKp0XLh7e0tXRdQ1ZIybNgwNG7cGM7OzggICAAAXLx4UeN897ZwtGvXDs2bN8fq1asBVCVEjRo1QpcuXWqMuaSkBHZ2djVuP3XqFPr3748333wTvXr1qu3ydbq3g3N1y9vd5brq3Ftmb28PAPf9dyQiuhtbSuqAg7UD1Elqk537QXTp0gWRkZF47bXXNFo/AN1fOrpuE1lbW2u8lslk0i2MTz/9FCUlJTrrBQYGIjAwEKNHj8a///1vBAUF4csvv8TIkSN1xnrhwgX069cP48ePxzvvvAM3Nzfs378fL7/8skEdSHVdl67y2q4LAKKjo+Hn54dPPvkEPj4+qKysRIsWLbQ67Do6ardgjR49GsuWLcO//vUvpKamYuTIkbWOfHJ3d6+xU/Lp06fx3HPPYcyYMXj99ddrPEZNvLy8tFo88vPzUb9+fTRo0KDWOve2nty4cQNAVcdiIiJ9MSmpAzKZzOBbKObgvffeQ6tWrRAUFKRRHhISgv3792uUHThwAEFBQbCystLr2E899ZRe9ZRKJRwcHKQWDxsbG43WGAA4evQo7ty5g+TkZNSrV9XYt2HDBo06uva7V0hICC5evIhLly5JrSWnT5+GSqVCcHCwXvFev34dOTk5WLlyJTp37gwAWu9VbV588UXMmjULS5YswalTp5CQkFBr/datW2Pt2rVa5adOncJzzz2HhIQE6baSoTp06IBvvvlGo+yHH37As88+KyVmHTp0wI4dOzBt2jSNOhERERr7ZWdnw9fXF+7u7g8UCxFZJiYlJGnZsiWGDx+OpUuXapQnJiaibdu2eOeddxAXF4eDBw9i2bJlWL58+UOdb/bs2SguLka/fv3g7++PmzdvYsmSJSgvL5duPSiVSuTm5iIrKwu+vr5wcnJCYGAg7ty5g6VLlyI6Oho//fST1nwaSqUSarUamZmZCAsLg4ODg9ZQ4J49eyI0NBTDhw9HSkoK7ty5g4kTJ6Jr1656dSYFAFdXVzRo0AAff/wxvL29cfHiRfzrX//S+z1wdXVFbGwsZs6cid69e8PX17fW+pGRkUhKSkJBQQFcXV0BVCUk3bt3R+/evTF9+nSpJcPKykqjpeL06dMoKyvDjRs3cOvWLWRlZQGANA/M+PHjsWzZMkyfPh1jxozBwYMHsWrVKnzxxRfSMaZMmYIuXbpg/vz56N+/P77++mvs3LlTKxHbt2+fVodYIqL7qttBQE+OJ3VI8L1DQ8+fPy9sbW1rHBJsbW0tGjVqJBYsWKCx/d6htkIIERYWJg3F1eXHH38UgwYNEn5+fsLGxkZ4enqKPn36iH379kl1bt++LQYNGiRcXFw0hgQvXLhQeHt7C3t7exEZGSlWr14tAIiCggJp3/Hjx4sGDRoYZUjw3RYtWiT8/f2l1zt27BDBwcHC1tZWhIaGit27dwsAIj09XQjxvyHBx48f1/k+ZGZmCgBiw4YNNb5XdwsPDxcrVqzQiBH3DIEGoBFj9bXrqne33bt3i9atWwsbGxuhVCrFRx99pHX+r776SjRt2lRYW1uLZs2aiU2bNmlsLykpEc7OzuLgwYM1XgOHBBNZFn2HBMuEMHAMqYUqLCyEQqGASqWCs7Ozxrbbt28jNzcXAQEBtXZCJNJl3bp1mDJlCq5evQobG5v71v/uu+8wY8YMZGdnS7evzMmHH36Ir7/+Gj/88EONdYz6mSkqAuT/He2mVgM6+u4QkWnV9h16N96+ITKR4uJi5ObmYt68eRg3bpxeCQlQNfT6jz/+wJUrVzRGDpkLa2trrVuARET6ML+fWUQW4v3330erVq3g6emJpKQkg/adMmWKWSYkADB27Fg0bdr00Z3QwaGqhUStrnpORI8tJiVEJjJ79myUl5cjMzMTcrlpJtt7IshkVbdsHB2rnhPRY4tJCREREZkFJiVGxD7DRPrhZ4WIdGFSYgTVE0txSm0i/VR/Vu6dLZeILBtH3xiBlZUVXFxcpPVQHBwcap0qnMhSCSFQXFyM/Px8uLi46D0jMBFZBiYlRlK9Gu3dC7URkW4uLi61ruBMRJaJSYmRyGQyeHt7w8PDw6BF4YgsjbW1NVtIiEgnJiVGZmVlxf9wiYiIHgA7uhIREZFZYFJCREREZoFJCREREZkF9inRU/VkT4WFhSaOhIiI6PFS/d15v4kTmZTo6datWwBgtougERERmbtbt25BoVDUuF0mON+zXiorK3H16lU4OTlxYjQiIiIDCCFw69Yt+Pj4oF69mnuOMCkhIiIis8COrkRERGQWmJQQERGRWWBSQkRERGaBSQkRERGZBSYlREREZBaYlBAREZFZYFJCREREZuH/A61CMK9o7RRPAAAAAElFTkSuQmCC", "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": 12, "metadata": { "ExecuteTime": { "end_time": "2024-12-11T20:54:43.238969Z", "start_time": "2024-12-11T20:54:43.188645Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\ospin\\OneDrive\\Documents\\Github\\xhydro\\src\\xhydro\\extreme_value_analysis\\parameterestimation.py:214: UserWarning: Dataset contains chunks. It is recommended to use scheduler='processes' to compute the results.\n" ] }, { "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.1816 ... 153.7\n",
       "    total_precip_lower  (station_num, dparams) float64 120B 0.08659 ... 135.6\n",
       "    total_precip_upper  (station_num, dparams) float64 120B 0.2766 ... 174.3
" ], "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)