{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Curve Fitting\n", "\n", "*Martin Vonk (2025)*\n", "\n", "It can be usefull to go from one soil model to the other. When the soil parameters are known the soil water retention curve and hydraulic conductivity function can be fitted." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "import pedon as pe" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def plot_compare(soilsample: pe.SoilSample, soilmodel: pe.SoilModel) -> None:\n", " f, ax = plt.subplots(1, 2, sharey=True, figsize=(4.2, 4.5))\n", " ax[0].scatter(soilsample.theta, soilsample.h, c=\"k\", s=10, label=\"Soil Sample\")\n", " _ = pe.soilmodel.plot_swrc(\n", " soilmodel, ax=ax[0], label=f\"Fitted Soil Model {soilmodel.__class__.__name__}\"\n", " )\n", " ax[0].set_yscale(\"log\")\n", " ax[0].set_xlim(0, 0.5)\n", " ax[0].set_yticks(soilsample.h)\n", " ax[0].set_xticks(np.linspace(0, 0.5, 6))\n", "\n", " ax[1].scatter(soilsample.k, soilsample.h, c=\"k\", s=10)\n", " _ = pe.soilmodel.plot_hcf(soilmodel, ax=ax[1])\n", "\n", " ax[1].set_yscale(\"log\")\n", " ax[1].set_xscale(\"log\")\n", "\n", " k_left = 10 ** (np.floor(np.log10(min(soilsample.k))) - 1)\n", " k_right = 10 ** (np.ceil(np.log10(max(soilsample.k))) + 1)\n", " ax[1].set_xlim(k_left, k_right)\n", " ax[0].set_ylabel(r\"|$\\psi$| [cm]\")\n", " ax[0].set_xlabel(r\"$\\theta$ [-]\")\n", " ax[1].set_xlabel(r\"$K_s$ [cm/d]\")\n", " ncol = 3\n", " ax[0].legend(\n", " loc=(-0.02, 1),\n", " fontsize=6,\n", " frameon=False,\n", " ncol=ncol,\n", " columnspacing=0.8,\n", " handlelength=2.5,\n", " )\n", "\n", " f.align_xlabels()\n", " # plt.close(f)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Genuchten(k_s=712.8, theta_r=0.045, theta_s=0.43, alpha=0.145, n=2.68, l=0.5)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sn = \"Sand\"\n", "soil = pe.Soil(sn).from_name(pe.Genuchten, \"HYDRUS\")\n", "soilm_genuchten = getattr(soil, \"model\")\n", "soilm_genuchten" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "h = np.logspace(-4, 6, num=11)\n", "k = soilm_genuchten.k(h)\n", "theta = soilm_genuchten.theta(h)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "SoilSample(sand_p=None, silt_p=None, clay_p=None, rho=None, th33=None, th1500=None, om_p=None, m50=None, d10=None, d20=None)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "soilsample = pe.SoilSample(h=h, k=k, theta=theta)\n", "soilsample" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Panday(k_s=np.float64(571.7082734658234), alpha=np.float64(0.1461801619285099), beta=np.float64(2.6526588374857303), brook=np.float64(3.797937441596011), sr=np.float64(0.10452835613889336))" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "soilm_panday = soilsample.fit(pe.Panday)\n", "soilm_panday" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAGwCAYAAACdL9N0AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZXxJREFUeJzt3XlcVWX+wPHP5bIvooKCKIq55gYIesVSccssNTTLyUnJyhpzmcIFmymXVjXHnDEmy8m0xTL9pWaLZaRZGhoq7gsq7gqCAgKy3vP7w/GOhMrdz73wfb9e96X33HPP833gOXzPec5zzqNRFEVBCCGEsDEXtQMQQghRO0jCEUIIYReScIQQQtiFJBwhhBB2IQlHCCGEXUjCEUIIYReScIQQQtiFJBwhhBB2IQlHCCGEXUjCEUIIYReScIQQQtiFJBwhhBB2IQnHwV29epXBgwcTGxtLTEwM33333W3XffbZZwGYNWsWX3/9daXP9u7dS8+ePenVqxfdu3fn3LlzVo3ziSeeYP/+/Vbdpi2cPHmSBg0aEBsbS2xsLD/99JPh57Z27VqysrIASEtLY8eOHUZvNzo6usqyt99+m27dutGjRw+ee+6523532bJl/Pbbb5w8eZLhw4dX+Vyj0fDhhx8a3rdv354pU6YYFdfw4cM5efLkLT+7XXl+fn7ExsYSFRXFt99+a1Q5t3Krn4mo3VzVDqCm2r59O0ePHqV169bodDqzt/PRRx9x//33M378eBRFIS8v77brvvfee7f97NVXX+Xdd9+lffv2XLt2DY1GY3ZMzq5Xr16sXr3a8L5Pnz7A9YTTsmVLGjZsSFpaGgUFBXTt2tWsMq5evcrKlSv57bff0Gg0XLly5bbrPvHEEwC3TQx33303GzZsYMyYMRw8eBAfHx+zYjJWmzZt2Lx5M2fPnmXAgAE88MADNi1P1B5yhmMDiYmJdOvWjdGjR9OtWzcSExPN3paXlxcpKSlkZmai0WioW7cuAAsWLCAmJoZ7772XXbt2AXc+ovTy8uLHH3+ksLAQLy8vPD09yczMpHfv3vTo0YPhw4dTUVHByZMn6d69OyNGjKB9+/asXLmSQYMGER4eTnp6OgBt27blscceIzo6mhUrVlQqR1EUJk6cSO/evenXrx9nz569ZTyKolBUWm6zlymzbkRHR5ORkWH4oz5t2jTeffdd/vnPf3LfffcB8MYbb9CrVy969uzJvn37APj444+Jjo7mscceo6CgoNI2XVxcyM7OZufOnSiKQr169QDYt28f9957L/fccw9vvvkmcOsz0pt5e3vj5uZGfn4+q1evrnRW8vnnn6PT6ejWrRvff/89AD/++COdO3dm2LBhhjPZ4uJiHn/8cfr06cOQIUPIz8+v9ufSpEkTioqK2LdvH7169SImJoYJEyYAsHnzZu6//36GDh1KeHi44ez2Vj+TjRs30qtXL7p06cKcOXMAGDZsmKFtvP/++yxZsqTaeEQNoAirSklJUYAqr5SUFLO2V1paqrz66qtKeHi40q1bN+Xw4cPKhQsXlB49eigVFRVKRkaG0q9fP0VRFCUqKkpRFEWZOXOmsn79+krbuXjxovLss88qLVq0UB599FGloKBAKSkpUcrKyhRFUZRJkyYpP/zwg5KRkaG0bdtWKS8vV77//nulc+fOSkVFhbJ27Vrl5ZdfVhRFUXx8fJScnByluLhYCQ8PV8rLy5X4+Hhl3759yvr16w3rpaSkKOPHj79lvQpLypRmiV/b7FVYUnbLcjMyMpTAwEClV69eSq9evZScnBzDz+1GHRRFUT788ENl0aJFiqIoyr59+5TRo0criqIo586dU4YMGaKUl5crERERSnFxsXLp0iXFx8enSlnff/+9cv/99yvNmzdX3nvvPUVRFGXQoEHKwYMHFb1er/Tv31/JyMgw/L4yMjKUhx9+uMp2oqKilM8++0z59NNPlQcffFD56aeflMmTJyvl5eVKp06dlGvXril5eXmGeuh0OsPvp2nTpkpGRoayaNEi5YMPPlAURVE+//xz5a233rpjeYqiKAcPHlSioqKUoqIiRa/XK4qiKEOGDFGOHj2qbNq0SenTp4+iKIry7bffKi+88MJtfyaFhYWKoihKRUWFEh0drRQVFSn/93//p8ydO1dRFEW57777lCtXrtzy9yVqFulSs7KjR4/edrk5XWtubm689NJLvPTSS2zcuJGZM2fy/PPPEx4ejouLC2FhYeTm5la7naCgIBYvXgzASy+9xMcff8xDDz3EuHHjuHLlCufPn6dz5860atWKdu3aodVqCQkJoUOHDri4uNC4cWN+/PFHAJo3b079+vUBCA0NJTs721DOwYMHWbNmDVu2bEFRFEJDQ02us639sUutOgcPHmTbtm3ExsYCoNVquXTpEk2aNMHDwwMPDw+aN29e5Xv33Xcf9913H1evXuWee+5h5MiRXLx4kbvvvhuAzp07c/z4caNiePDBBxkwYAA6nc7QHXrp0iWaNm2Kp6cnnp6euLm5UV5eTkVFheH306lTJ0Mdfv/9dz766CPKysro0aPHbcs6cuQIsbGxuLq68u6775KRkcHkyZMpKirixIkTnD9/HoCIiAjgehu4cuXKbX8mO3fuZPbs2ZSVlXHy5EmysrIYNGgQAwYM4E9/+hP+/v6GM3dRs0nCsbLWrVubtLw6p06dolGjRri7u9OwYUMURSEsLIy0tDT0ej2nT582amdNT0+nVatWADRo0ABFUVixYgWDBg3i6aefZuLEiYZuqJuv79z8/xufnzx5kitXruDt7c2ZM2cIDAw0rNO2bVseffRRXn75ZQDKyspuGY+Xm5aDrwww7YdhAi83rcnfcXNzo6Kiosr/27ZtS69evfjPf/4DXK+Ti4sLZ8+epbS0lIKCAjIyMiptq7i4mEuXLhEaGoqvry+enp7A9cR/6NAh2rZty65du/jLX/7CL7/8Um1sfn5+DBw4kEGDBhmu4zVo0IBTp05RXFxMaWkppaWluLq6otVqDb+fG91/bdu2JSYmhlGjRhnqcLuBIzeu4dwwceJEJk+eTL9+/RgyZMgt24miKDRo0OCWP5N58+axePFi7rrrLjp37oyiKLi7u9O+fXsSExP585//XG39Rc1QKxNORkYGTz75JJmZmWi1WlJSUqx2IVan0zFt2jTmzZtnWJaYmGj2wIF9+/YxYsQIPD09URSFpKQkgoODeeihh+jevTsuLi4sWrSo2u18/vnnfP3113h5eVG3bl0++eQTjh07xqhRo1i/fj1eXl5GxxQaGsqkSZM4dOgQU6ZMQav93x/3wYMH89NPP9G7d280Gg1//vOfeeqpp6psQ6PR4O3uWM1v4MCBPP/88/Tr14/HHnuM0aNHs337dlasWEGrVq3o1asXLi4u9O/fn7/97W88//zzdO/enbZt29K0adNK2yorK2PMmDEUFxdTUVHBqFGj8PX15fXXX+fpp59GURQefPBBwsLCjI7vRhK/kQy0Wi3Tp0+nZ8+euLi48NprrwHwyiuv0LdvX8LCwgxxPfPMMzzzzDOG0W6TJ0+mffv2RpU7ePBg/vrXv9K2bVv0ev1t19Nqtbf8mTz88MMMHTqUjh074ufnZ1g/Pj6e+++/n48++sjon4FwbhpFMeHqag3Rq1cvXnvtNXr06MHly5epU6cOrq7W/eNnrVFqjig6OprU1FS1wxBObufOnXz44Ye88847aoci7MSxDjHt4MCBA7i5uRn6sG/0dVubTqercYlGCGtZu3Ytc+bM4dNPP1U7FGFHTjcsesuWLQwePJiQkBA0Gg1r166tsk5SUhJhYWF4enqi0+kq3cCXnp6Or68vgwcPpnPnzrzxxht2jL5mkLMbYam4uDhSUlJo0aKF2qEIO3K6hFNYWEh4eDhJSUm3/HzlypUkJCQwc+ZMdu3aRXh4OAMGDDDcQV5eXs4vv/zCv//9b3777Tc2btzIxo0b7VkFIYSolZz6Go5Go2HNmjXExcUZlul0Orp06WLoF9br9YSGhjJx4kSmT5/Ob7/9xqxZsww3yb311lsATJ069ZZllJSUUFJSYniv1+u5fPkyAQEBtfpufWEdiqJw9epVQkJCcHGxzfGftGFhSya1Yfvf+mM9gLJmzRrD+5KSEkWr1VZapiiKMnr0aGXIkCGKoihKWVmZEhERoVy+fFmpqKhQBg0aVOUmyZvNnDnzljdyykte1nydOXPGFruItGF52e1lTBuuUYMGsrOzqaioICgoqNLyoKAgDh8+DICrqytvvPEGPXv2RFEU7rvvPgYNGnTbbb744oskJCQY3ufl5dG0aVOOHj1a7YCD3bt3EzfqGer/aS764gJyll1/gOP69euJjIw0ul5lZWVs2rSJ3r174+bmZvT3LKVGubWprgCXL1+mdevWlYYLW1t1bXj37t0MHjwYj+bR1BkwCUWvJ3fda6z5z9smtVM1qPV7swVnrcvVq1dp3ry5UW24RiUcYw0cOJCBAwcate6Nu6aTkpJISkoy3AxYv359AgIC7vjdfv36MWr0aL7TeKPoKyguLiYxMZF+/fqZFG9ZWRne3t4EBATY/Y+wvcutTXW9mS27tqprw/369WPSpEnMmzcP17u64tuhD83+NJOYHrH4eDj2nwi1f2/W5Kx1uRGrMW3Y6QYN3ElgYCBarZbMzMxKyzMzMwkODrZo2+PHjzc8HsQUkydff4y8r68PKSkphocXCmFvd2rDc+fOJSUlhTeGRxLg5UIBXrzx7SEVohQ1WY1KOO7u7kRFRZGcnGxYptfrSU5OJiYmxqJtJyUl0a5dO7p06WJmbB5yX45QVXVtWKfT8cwTj/OvP1///NPtp9l0JMueIYoazukSTkFBAWlpaaSlpQHXH1OTlpbG6dOnAUhISGDJkiUsX76cQ4cOMW7cOAoLCxkzZoxF5Zp7hiOEozC2Dd/TMpAx94QBkLh6L1cKS+0QnagNnC7hpKamEhkZabiYmZCQQGRkJDNmzABgxIgRzJ8/nxkzZhAREUFaWhobNmyoMpDAVJae4QihNlPacOL9bWnRwIesqyW8vM7xZ3IVzsHpEk5sbCyKolR5LVu2zLDOhAkTOHXqFCUlJWzfvt0qXVlyhiOcnSlt2NNNy4JHI9C6aPh67wXWpVl3SnJROzldwlGLnOEIZ2dqGw4PrcvEPi0BeHntfi7mFdsyPFELSMIxkpzhCGdnThse37slnZr4k19cztTVe0yauluIP5KEI4S4LTetCwsejcDD1YVf0rP5JOWU2iEJJyYJx0jSpSacnbltuGVDX6YPbAvA698e4sSlAluEJ2oBSThGki414ewsacPxMWHc0zKA4jI9CV/sobzi9jN/CnE7knCEENVycdHw1vBw/DxdSTuTy+Kfj6sdknBCknCEEEYJqevF7CHtAVj4Yzr7z+WpHJFwNpJwjCTXcISzs0YbHhrZmIEdginXK7ywMo3isgorRihqOkk4RpJrOMLZWaMNazQaXh/akUBfD9KzCpj//RErRihqOkk4QgiT1PdxZ+7DHQH4YGsGvx3PUTki4Swk4QghTNb37iD+1CUURYEpq/ZwtbhM7ZCEE5CEI4Qwy0uD2hFa34tzudd4Zf1BtcMRTkASjpFk0IBwdtZuw74ervzjkQg0Gli18yw/HLhole2KmksSjpFk0IBwdrZow12b1+eZHncB8OKX+8guKLHatkXNIwlHCGGRhPta0ybIj5zCUl78cp884FPcliQcIYRFPFy1vD0iAjetho0HM1m986zaIQkHJQlHCGGxdiF1eKF/awBmrz/I2StFKkckHJEkHCGEVTzbswVRzepRUFLOlFV70Oula01UJglHCGEVWhcNCx4Nx9tdS8qJyyzdmqF2SMLBSMIxkgyLFs7OHm24WYAPf3/wbgDmfX+E9MyrNitLOB9JOEaSYdHC2dmrDY/s2pTYNg0oLdfzwhdplJbL3DniOkk4Qgir0mg0zHu4E3W93dh/Lp93fkpXOyThICThCCGsrmEdT16L6wBA0ubj7D59ReWIhCOQhCOEsIlBnUJ4KCKECr1Cwhd7uFYqc+fUdpJwhBA288qQDgTX8SQju5A3vzukdjhCZbUy4YSFhdGpUyciIiLo3bu32uEIUWP5e7vx1iOdAPjot1NsOXpJ5YiEmmplwgHYtm0baWlpbNq0Se1QhKjRerRqwOiYZgBMW72XvCKZO6e2qrUJx95KS0vYvn272mEIoYrpA9vSPNCHi/nFPLf0Zz7++GPZH2ohp0s4W7ZsYfDgwYSEhKDRaFi7dm2VdZKSkggLC8PT0xOdTseOHTsqfa7RaOjVqxddunTh008/tWm8//jHfAAKCgrp1q0biYmJNi1PCEfk7e7KgkfD0SgKW8+W8Jc33pf9oRZyVTsAUxUWFhIeHs6TTz7JsGHDqny+cuVKEhISWLx4MTqdjoULFzJgwACOHDlCw4YNAfj1119p3LgxFy5coF+/fnTs2JFOnTrdsrySkhJKSv43x0d+fj4AZWVllJXduWsgNTWVTz75hIBRXdBoNHh5ebFo0SLi4uKIjo42us43yqmuPGtTo9zaVFd7lWdJG7am4nOHKfj9//DpOpyA+8ajuXTcrP3hZmr93mzBWetiSrwaxYknr9BoNKxZs4a4uDjDMp1OR5cuXXjnnXcA0Ov1hIaGMnHiRKZPn15lG1OnTqV9+/Y88cQTtyxj1qxZzJ49u8ryFStW4O3tXW2MmdfgjTRXvLUKb3aVYaGisqKiIkaOHEleXh516tSxSRmWtmFrKtfD2/u1nC3UcHddPc+21aPR2DUEYWWmtOEalXBKS0vx9vZm9erVlZJQfHw8ubm5rFu3jsLCQvR6PX5+fhQUFNCrVy8WL1582+dL3eroMDQ0lAsXLhAQEHDH+FJTU7lv+CgCRv0TfXEB2e+PASA5OdnkM5yNGzfSv39/3NzcjP6epdQotzbVFSAnJ4dGjRrZNOFY0oatKTU1lb59+6Kt34T6f5qLxtWd/J/e45t/vmjRGY4avzdbcNa65OfnExgYaFQbdroutTvJzs6moqKCoKCgSsuDgoI4fPgwAJmZmQwdOhSAiooKxo4de8eHGXp4eODh4UFSUhJJSUlUVFw/S3Fzc6u2UcTExPD444/zHaAoCteuXSMxMZGYmBiz6mdMmbagRrm1pa72KMuSNmxNMTExTJw4kXnz5qH5eTn1+46lXp+nadSqk8VxqNVebMHZ6mJKrDUq4RjjrrvuYs+ePSZ/b/z48YwfP578/Hz8/f2N/t7kyVP4bsHP+Pr6kJKSgk6nM7lsIazB3DZsTXPnzmXYsGEcPnKUlVmuHMyGyav28MWzMWhdpG+tpnO6UWp3EhgYiFarJTMzs9LyzMxMgoODLdq2pY92d3f3kGQjVOUoU2zodDriR4/i/ad64Ovhys5TV3hvy3FVYxL2UaMSjru7O1FRUSQnJxuW6fV6kpOTze7GukGmJxDOztHacJN63swc3A6Atzce5eD5fJUjErbmdAmnoKCAtLQ00tLSAMjIyCAtLY3Tp08DkJCQwJIlS1i+fDmHDh1i3LhxFBYWMmbMGIvKdZSjQyHM5YhteHhUE/q3C6KsQiHhizRKymUkZ03mdAknNTWVyMhIIiMjgesJJjIykhkzZgAwYsQI5s+fz4wZM4iIiCAtLY0NGzZUGUhgKkc7OhTCVI7YhjUaDW8O60iAjzuHL15lwcajaockbMjpBg3ExsZS3UjuCRMmMGHCBDtFJISwRKCvB28O68gzH+/k/S0n6Ns2iK7N66sdlrABpzvDUYsjdkcIYQpHbsP3tQ/mkagmKApMXpVGQUm52iEJG5CEYyRH7I4QwhSO3oZnDG5H47penLl8jde/Oah2OMIGJOEIIRyCn6cb8x8JR6OBz3acIflQZvVfEk5FEo6RHLk7QghjOEMbjmkRwFP3NAcg8f/2cbmwVOWIhDVJwjGSo3dHCFEdZ2nDUwa0oVVDX7ILSvj7mn3VDhISzkMSjhDCoXi6aXl7RASuLhq+23+RtWnn1A5JWIkkHCGEw+nQ2J+/9m0FwIx1Bzife03liIQ1SMIxkjP0fwtxJ87WhsfFtiAitC5Xi8uZunoPer10rTk7SThGcpb+byFux9nasKvWhQWPhuPp5sLWYzl89NtJtUMSFpKEI4RwWHc18OVvD9wNwJvfHeZYVoHKEQlLSMKxMZk+VwjLjOrWjB6tAikp15PwRRplFXq1QxJmkoRjJEv7v2Vop1Cbs13DuUGj0fDW8HDqeLqy92weSZuOqR2SMJMkHCM5W/+3EH/kzG042N+TV+M6ALDop2PsPZurbkDCLJJwhBBOYUh4CA92akSFXuGFlWkUl8ncOc5GEo4QwiloNBpee6gDDf08OH6pkLkbDqsdkjCRJBwhhNOo5+PO3OGdAPhw60l+O5GjckTCFJJwhBBOpXebhozUNQUg8csDFMnUOU5DEo4Qwun8/YG7aRbgzYW8Yr48KX/GnIX8pozkrENKhbihJrVhHw9XFjwajosGfr/kwvcHZO4cZyAJx0jOPKRUCKh5bTiqWX2e6XF97pyXvzpI1tVilSMS1ZGEI4RwWhN7t6Cxt8KVojJe/D+ZO8fRScKxMXmyjRC24+7qwuMtK3DTakg+nMUXqWfUDkncgSQcIYRTC/GBF/q1BOCV9Qc5c7lI5YjE7UjCEUI4vSe7h9E1rD6FpRVM/mIPFTJ3jkOqtQmnqKiIZs2aMWXKFLuUJ81fCNvRumiY/0g4Pu5adpy8zAe/nlA7JHELtTbhvP7663Tr1k3tMIQQVtI0wJuXB7UDYP73Rzl8MV/liMQf1cqEk56ezuHDhxk4cKDaoQghrGhEl1D6tm1IaYWeF1buobRc5s5xJE6XcLZs2cLgwYMJCQlBo9Gwdu3aKuskJSURFhaGp6cnOp2OHTt2VPp8ypQpvPnmm3aK+Lqy0lK2b99u1zKFqG00Gg1vPtyRet5uHLqQz7SPNvPxxx/LvucgnC7hFBYWEh4eTlJS0i0/X7lyJQkJCcycOZNdu3YRHh7OgAEDyMrKAmDdunW0bt2a1q1b2yXe+fPnA9evGXXr1o3ExES7lCtEbdXQz5M3hnYEYM3hAsa++Ibsew7CVe0ATDVw4MA7doUtWLCAsWPHMmbMGAAWL17MN998w9KlS5k+fTopKSl8/vnnrFq1ioKCAsrKyqhTpw4zZsy45fZKSkooKSkxvM/Pv94vXFZWRllZ2R1jTU1N5eOPPyZgdBfQaPDy8mLRokXExcURHR1tdJ1vlFNdedamRrm1qa72Ks+SNuzobvd7q1twkmuHfsbr7l4EDprClc+mmrXv2ZNabdBSpsSrUZz41lyNRsOaNWuIi4sDoLS0FG9vb1avXm1YBhAfH09ubi7r1q2r9P1ly5axf/9+w1nIrcyaNYvZs2dXWb5ixQq8vb2rjTHrGrye5oqXVmFOV5kwSlRWVFTEyJEjycvLo06dOjYpw9I27KyKymHuHi25pRruCdLz6F1yPccWTGnDTneGcyfZ2dlUVFQQFBRUaXlQUBCHD5s3WdOLL75IQkKC4X1+fj6hoaH07t2bgICAO343NTWVxx9+nIDR/6KwqIjHHnsCgOTkZJPPcDZu3Ej//v1xc3Mzqx7mUKPc2lRXgJwc28/nYkkbdnS3+72lpqbSt29f3Jp0oN6wmWzNdOGb99/ku6X/cOgzHDXaoKVunDEbo0YlHFM98cQT1a7j4eGBh4cHSUlJJCUlUVFx/SzFzc2t2kYRExPDqFGP8y2AonDt2jUSExOJiYkxK15jyrQFNcqtLXW1R1mWtGFn8ce6xMTEMHHiRObNm4c29SvqRA8h+KEptO0U5fB1drbfiymxOt2ggTsJDAxEq9WSmVn5UeWZmZkEBwdbtG1zn7Q7deo0ALy9vUlJSWHOnDkWxSGEuWra06KrM3fuXFJSUpg/6l5CfLUU48HL6/arHVatVqMSjru7O1FRUSQnJxuW6fV6kpOTzT6ruMHSuUTc3N3R6XQWxSCEJWrSfDjG0ul0PBk/infju6F10fD13gt8tee82mHVWk6XcAoKCkhLSyMtLQ2AjIwM0tLSOH36NAAJCQksWbKE5cuXc+jQIcaNG0dhYaFh1Jq5atvRoah5anMbDg+ty4Te1x/w+dKafVzMk7lz1OB0CSc1NZXIyEgiIyOB6wkmMjLSMKx5xIgRzJ8/nxkzZhAREUFaWhobNmyoMpDAVLXx6FDULLW9DU/o05JOTfzJLy5n6uo9MneOCpwu4cTGxqIoSpXXsmXLDOtMmDCBU6dOUVJSwvbt263SlVWbjw5FzVDb27Cb1oUFj0bg4erCL+nZfJJySu2Qah2nSzhqsfjoUA6mhMpq+xkOQMuGvkwf2BaA1789REZ2ocoR1S6ScIxU248OhfOTNnxdfEwY97QMoLhMzwsr0yivkBtC7UUSjhCiVnFx0fDW8HD8PF1JO5PL4p+Pqx1SrSEJx0jSHSGcnbTh/wmp68XsIe0BWPhjOvvP5akcUe0gCcdI0h0hnJ204cqGRjbm/vbBlOsVXliZRnGZPOvQ1iTh2JhG7QCEELek0Wh4fWgHAn09SM8q4B8/HFE7pBpPEo4QotYK8PVg7sPX5875z68ZpJyw/cNUazNJOEaS/m/h7KQN31rfu4P4U5dQFAUmf7GHq8XONR+NM5GEYyTp/xbOTtrw7b00qB2h9b04l3uNV9YfVDucGksSjhCi1vP1cOUfj0Sg0cCqnWf54cBFtUOqkSThCCEE0LV5fZ7pcRcAL365j+yCkmq+IUwlCUcIIf7rhf6taRPkR05hKX/7cp884NPKJOEYSS64Cmcnbbh6nm5aFowIx02r4YeDmfzfrnNqh1SjSMIxkqUXXOU4SahNBg0Yp32IP8/3aw3ArK8OcPZKkcoR1RyScIQQ4g/+0qsFUc3qUVBSzpRVe9Dr5ZDRGiThCCHEH2hdNPzjkXC83bWknLjM0q0ZaodUI0jCsTGNPNtGCKcUFujD3x+8G4B53x8hPfOqyhE5P0k4QghxGyO7NqVX6waUlut54Ys0ymTuHItIwhFCiNvQaDTMG96Jut5u7D+Xz6LkdLVDcmqScIwkQ0qFs5M2bJ6gOp68FtcBgKTNx9l9+orKETkvSThGkiGlwtlJGzbfoE4hDAkPoUKvMPmLPVwrlblzzCEJRwghjPDqQx0IquPBiexC5nx3SO1wnJIkHCGEMIK/txtvDQ8HYPlvp/gl/ZLKETkfSThCCGGknq0bMKpbMwCmrtpLXpHMnWMKSThCCGGCFx9oS/NAHy7mFzPzq/1qh+NUal3Cyc3NJTo6moiICDp06MCSJUvUDkkI4US83V35x6PhuGhgbdp5vtl7Qe2QnEatSzh+fn5s2bKFtLQ0tm/fzhtvvEFOju3nMZfHnAtRc3RuWo/xvVsC8Pe1+8jKL1Y5IudQ6xKOVqvF29sbgJKSEhRFsWkycPnvs23KysvZvn27zcoRQtjXxD6taB9Sh9yiMp75YAsfffSx7OPVcLqEs2XLFgYPHkxISAgajYa1a9dWWScpKYmwsDA8PT3R6XTs2LGj0ue5ubmEh4fTpEkTpk6dSmBgoM3inTt3DgDFJSV069aNxMREm5UlhLAfd1cX3h4RgYuiJy2zjPELPpV9vBquagdgqsLCQsLDw3nyyScZNmxYlc9XrlxJQkICixcvRqfTsXDhQgYMGMCRI0do2LAhAHXr1mXPnj1kZmYybNgwhg8fTlBQ0C3LKykpoaTkf1PN5ufnA1BWVkZZ2Z1HqKSmprJs2TICx+jQaFzw8vJi0aJFxMXFER0dbXSdb5RTXXnWpka5tamu9irPkjbs6NT6vd2Qc2I/eb9+jF+PeOr1GQuZR8zax0H9upjLlHg1ihNfXNBoNKxZs4a4uDjDMp1OR5cuXXjnnXcA0Ov1hIaGMnHiRKZPn15lG8899xx9+vRh+PDhtyxj1qxZzJ49u8ryFStWGLrm7uRKCcza5YqrRuEf3eTuZFFZUVERI0eOJC8vjzp16tikDEvbsLgzvQJJB7Ucy9fQ3E9hUvsKXGrRU+JNacM1KuGUlpbi7e3N6tWrKyWh+Ph4cnNzWbduHZmZmXh7e+Pn50deXh733HMPn332GR07drxlGbc6OgwNDeXChQsEBATcMb7U1FT6P/QIgU++h1JRzqWkxwBITk42+Qxn48aN9O/fHzc3N6O/Zyk1yq1NdQXIycmhUaNGNk04lrRhR6fW7+2G1NRU+vbti4tfA+r/eT4u7t4UbP2E9fMmmXWGo2ZdzJWfn09gYKBRbdjputTuJDs7m4qKiirdY0FBQRw+fBiAU6dO8cwzzxgGC0ycOPG2yQbAw8MDDw8PkpKSSEpKoqLi+lmKm5tbtY0iJiaGJ8eM4av/vr927RqJiYnExMSYVT9jyrQFNcqtLXW1R1mWtGFnoVZdYmJimDhxIvPmzYMf3yfwgeepc89I/Ju1NzseZ/u9mBJrjUo4xujatStpaWkmf2/8+PGMHz+e/Px8/P39jf7eS3//O1+9kYyLVktKSgo6nc7ksoWwBnPbsLizuXPnMmzYMI4cOcpXV9xJvVBKwhdprJtwDx6uWrXDcyhON0rtTgIDA9FqtWRmZlZanpmZSXBwsEXbNvfR7pobU35qNJJshKpkegLb0el0jB49isVP9STAx53DF6+yYONRtcNyODUq4bi7uxMVFUVycrJhmV6vJzk52exurBvMfbT7jXzjvFfKRE0h0xPYXqCvB28Ou95F//6WE/x+8rLKETkWp0s4BQUFpKWlGbrFMjIySEtL4/Tp0wAkJCSwZMkSli9fzqFDhxg3bhyFhYWMGTPGonLNPTq8ceMnyNMGhLrkDMc+7msfzPCoJigKJHyRRkFJudohOQynSzipqalERkYSGRkJXE8wkZGRzJgxA4ARI0Ywf/58ZsyYQUREBGlpaWzYsOG299kYy+wznJv+r5d8I1QkZzj2M2NwOxrX9eLM5Wu8/s1BtcNxGE43aCA2NrbaM4UJEyYwYcIEO0V0Z1XPcGrRAH0haqk6nm7MfyScx5ak8NmOM/RvF0SftpYd9NYETneGoxazuyNuyi9yhiPUJF1q9hXTIoCn7m0OwLTV+7hcWKpyROqThGMkSwcNAChIxhHqkS41+5s6oA2tGvqSXVDC39fsq/XXcSXh2FjlLjUVAxFC2J2nm5a3R0Tg6qLhu/0XWZt2Tu2QVCUJx0hm34dz0/8l4Qg1SZeaOjo09uevfVsBMGPdAc7nXlM5IvVIwjGSud0Rlc5wpEtNqEi61NQzLrYF4aF1uVpcztTVe9DX0gu6knBs7OZrOBW1tJEJUdu5al14+9FwPN1c2Hosh49+O6l2SKqQhGNjrjc9p1wSjhC1110NfPnbA3cD8OZ3hzmWVaByRPYnCcdI5vZ/a29KOKUVemuHJYTR5BqO+h7XNaNHq0BKyvVM/iKNslr2N0ESjpHMHxatwV17/cdcXiFnOEI9cg1HfS4uGuYN70QdT1f2nM3j35uOqx2SXUnCsQNX7fWzHEk4QohG/l68GtcBgH/9lM7es7nqBmRHknDs4MZ1nDJ97Tp9FkLc2pDwEB7s1IgKvcILK9MoLqsd089LwjGSJf3fbtKlJhyAXMNxHBqNhtce6kADPw+OXypk7obDaodkF5JwjGRJ//eNLrXadoFQOBa5huNY6vm4M294JwA+3HqS307kqByR7UnCsQNXl+s/Zkk4Qoib9W7TkJG6pgAkfnmAazV86hxJOHbg7vrfLjW5D0cI8Qd/f+BumgV4cyGvmP87WbP/JNfs2jkIt/92qZWWyxmOEKIyHw9X/vFIOC4a+P2SC98fyFQ7JJuRhGMH3u7X57krKq0dI1GEEKaJDqvP2P/OnfPyVwfJulqsckS2IQnHDrzdtQAUldbwDlohhNkm9WlBY2+FK0Vl/O3Lmjl3jiQcI1kypFTOcIQjkGHRjs3d1YU/t6zATavhx0NZfJF6Ru2QrE4SjpEsGVLq43H9DKewRM5whHpkWLTja+wDz/dtCcAr6w9y5nKRyhFZlyQcO5AzHCGEsZ66J4wuYfUoLK1g8hd7atRT5iXh2MH/ruFIwhFC3JnWRcM/HonAx13LjpOX+eDXE2qHZDWScOzAz/P6GU5+cZnKkQghnEHTAG9eHtQOgPnfH+XwxXyVI7IOSTh2UN/HHYDLBaUqRyKEcBYjuoTSt21DSiv0JKzcUyPu43M1ZeWvvvrK5AL69++Pl5eXyd+zlTNnzjBq1CiysrJwdXXl5Zdf5pFHHrFpmYaEUyQJRwhhHI1Gw5sPd2TA21s4eCGffyWnM2VAG7XDsohJCScuLs6kjWs0GtLT07nrrrtM+p4tubq6snDhQiIiIrh48SJRUVE88MAD+Pj42KzM+t7/TTiFknCEEMZr6OfJ60M78tynu/j35mP0ubshnZvWUzsss5ncpXbx4kX0er1RL29vb1vEbJFGjRoREREBQHBwMIGBgVy+fNmmZdb3lYQjhDDPAx0bERcRgl6ByV/sceobyE1KOPHx8SZ1jz3++OPUqVPH5KDuZMuWLQwePJiQkBA0Gg1r166tsk5SUhJhYWF4enqi0+nYsWPHLbe1c+dOKioqCA0NtWqMf9TQzxO4nnB+2ZZi07KEEDXP7CEdCK7jSUZ2IWMWrmP79u1qh2QWkxLOhx9+iJ+fn9Hrv/vuuwQGBpoc1J0UFhYSHh5OUlLSLT9fuXIlCQkJzJw5k127dhEeHs6AAQPIysqqtN7ly5cZPXo077//vlXju5U5r7yMvvQaAH0GPUxiYqLNyxRC1Bz+3m60ybt+4Lz9siexf/qLU/4dMekajiMYOHAgAwcOvO3nCxYsYOzYsYwZMwaAxYsX880337B06VKmT58OQElJCXFxcUyfPp3u3bvfsbySkhJKSkoM7/Pzrw9PLCsro6ys+mHOqampvLNoEfX//A9cApri07ApixYtIi4ujujo6Gq/f6Osm/+1FzXKrU11tVd5lrZhR6bW780W7lSX1NRUVi16Fd9eT+Edfj+BD77AO+9PNunviK2Y8rPXKBY8Ia64uJi9e/eSlZWFXl95yN6QIUPM3azRNBoNa9asMQxmKC0txdvbm9WrV1ca4BAfH09ubi7r1q1DURRGjhxJmzZtmDVrVrVlzJo1i9mzZ1dZvmLFCpOuUb1/2IUDV1x49K4K7gmqOXcOC8sUFRUxcuRI8vLyrN79fIO12rBQX0kFvLVXy6ViDdGBeka1Un+otClt2OyEs2HDBkaPHk12dnbVjWo0VFTY/q76Pyac8+fP07hxY7Zt20ZMTIxhvWnTpvHzzz+zfft2fv31V3r27EmnTp0Mn3/88cd07NjxlmXc6ugwNDSUCxcuEBAQUG2Mqamp9O3bF9+eY/COeICiXesp+PUjkpOTTTrD2bhxI/3798fNzc2o71iDGuXWproC5OTk0KhRI5smHEvbsCNT6/dmC3eqy42/IwCuwa2oN/w1NC4uPN/Fl/FD7txLY2v5+fkEBgYa1YbN7lKbOHEijzzyCDNmzCAoKMjczdjdvffeW+Vs7E48PDzw8PAgKSmJpKQkQyJ1c3MzqoHHxMQwceJE/r1hN94RoKkfyqRJkyolRGMZW6a1qVFubamrPcqytA07g5pelxt/R+bNmwcZe9GmrMK/+wiWHyhh5H0VhoFJajDl5272kwYyMzNJSEhwqGQTGBiIVqslM7PyjHmZmZkEBwdbtG1LnrQ7d+5c3p87E4CGrSN58803LYpFCHPI06Kd29y5c0lJSeGjjz7i67nPcXejOlwpKuPF/3OeuXPMTjjDhw9n8+bNVgzFcu7u7kRFRZGcnGxYptfrSU5ONuuM4maWziUytG83tC4arpYqXMirmbP5Cccm8+E4P51Ox6hRo7gnphtvjwjHXetC8mHnmTvH7C61d955h0ceeYRffvmFjh07VjmtmjRpksXB3UpBQQHHjh0zvM/IyCAtLY369evTtGlTEhISiI+PJzo6mq5du7Jw4UIKCwsNo9bMNX78eMaPH09+fj7+/v4mf9/TTUvbYD8OnM/n95OXeSiisUXxCGEqS9uwcCxtg+sw+b7WvPndYV5Zf5DuLQIJre/Yg0DMTjifffYZP/zwA56enmzevBmNRmP4TKPR2CzhpKam0rt3b8P7hIQE4PpItGXLljFixAguXbrEjBkzuHjxIhEREWzYsMHirr8/9n+bo3uLAA6cz2fbsRxJOMLurNGGhWN5usdd/Hgok99PXmHyqj18PrYbLi6a6r+oErO71P7+978ze/Zs8vLyOHnyJBkZGYbXiRO2m78hNjYWRVGqvJYtW2ZYZ8KECZw6dYqSkhK2b9+OTqezuFxr9H93b3n9Jtitx7Odps9V1BxyDafmuTF3jre7lh0Zl1m6NUPtkO7I7IRTWlrKiBEjcHGpHTMcWKP/u2tYfdxdXTh75RqHLly1YnRCVE+u4dRMN8+dM+/7IxzNdNy/LWZni/j4eFauXGnNWByaNY4OfTxc6d2mAQDr9563VmhCGEXOcGquP3UJJbZNA0rL9bywMs1h584x+xpORUUF8+bN4/vvv6dTp05VBg0sWLDA4uBqosHhIXx/IJO1u88xuX9rXLW14wxRCGE7Go2GeQ934r6FWzhwPp93fkon4T7HmzvH7L92+/btIzIyEhcXF/bv38/u3bsrvWoaa3VH9Ls7iPo+7lzIK+b7A5nVf0EIK5EutZqtYR1PXovrAEDS5uPsPn1F5YiqMvsMZ9OmTdaMw+FZa0ipp5uWx3VN+ddPx/jPryd4oGNwpRF+QtiKDIuu+QZ1CuGHA5l8tec8k7/YwzeTeuDlrlU7LAOzz3DefPNNli5dWmX50qVLmTt3rkVB1XSPd2uGh6sLu0/nsvGgnOUIIaznlYfaE1THgxPZhczdcFjtcCoxO+G89957tG3btsry9u3bs3jxYouCquka1vHkqXubAzDnu8OUVTjmBT4hhPOp6+3OvOHhACzbdpKtx6o+YFktZiecixcv0qhRoyrLGzRowIULFywKyhFZu/97XGwLAnzcOZFdyHs/H7fKNoW4E7mGU3v0at2Ax7s1BWDKqj3kXXOM+YLMTjihoaFs3bq1yvKtW7cSEhJiUVCOyNpDSv083Xhp0N0ALPwxnf3n8qyyXSFuR4ZF1y5/e+BumgV4cyGvmNnrD6gdDmBBwhk7dizPP/88H374IadOneLUqVMsXbqUF154gbFjx1ozxhorLqIxA9oHUa5XGPfpTq4UlqodkhCihvB2d2XBo+FoNPDlrnP86ADXi80epTZ16lRycnJ47rnnKC29/ofS09OTxMREXnzxRasFWJNpNBrmPtyJQxe2cvpyEc99uotlT3bBw9VxRpUIIZxXVLP6jO1xF+9vOcGLa/YRHVaPut7uqsVj9hmORqNh7ty5XLp0iZSUFPbs2cPly5eZMWOGNeOr8ep6u7NkdDQ+7lp+O5HD+E93OexdwkII55PQvzUtGvhw6WoJs9cfVDUWi29z9/X1pUuXLnTo0AEPDw9rxOSQbHnBtU2wH0tGR+Ph6sKPh7IYv2IX10rlib7CumTQQO3k6aZl/iPhuGhgze5zfH/gomqxmJRw9u7da9L0zAcOHKC8vNzkoByRrS+4dm8ZyPujo3HXurDxYCaPLUkhu6Ck+i8KYSQZNFB7RTatxzM9WwDw9zX7uKzS9WKTEk5kZCQ5OTlGrx8TE8Pp06dNDqq26tW6AZ88rcPfy420M7k89M5Wdjng4ymEEM7n+X6taNXQl+yCUmZ+pc6oNZMGDSiKwssvv4y3t3Gzyt0YTCCM17V5fb58rjtPLvudUzlFPLr4N57v25LGMn2OEMICN7rWhr27jfV7zvNAh2AGdqx6L6UtmZRwevbsyZEjR4xePyYmBi8vL5ODqu1aNPDl64n38rc1+1m/5zzzN6YT5qulVdRVOoTWVzs8IYSTCg+ty7heLXhn0zFeWrufrs3rE+Brv2vvJiWczZs32ygM8Ud+nm78608R9GgZyOyvD3CyoIK4d1N4qkdzJvZpha+H2SPahRC12MS+LfnxUCaHL15lxlcHSBrZ2W5ly2QsDkyj0fBol1C+m3gPnerrKdcrvPfzCXrN28SyrRkyfFoIYTIP1+tda1oXDd/svWDXG0Il4RhJzSGljfw9eaqNnsV/jiAswJucwlJmrT9I3wWbWb3zrCQeYRQZFi1u6NDYn6d7XH+A8Mvr9nO12D7PWpOEYyRHGFLat21DNib04rW4DjTw8+DM5WtMWbWHXm9tYsmWE3ZrNMI5OUIbFo7j+b6taVr/+rPW5n9v/LV5S0jCcTJuWhce79aMn6fGMn1gWxr4eXAhr5jXvz1E9zk/8ca3h8jILlQ7TCGEg/Ny1/LG0I4AfJRyip2nbH8LhiQcJ+Xt7spferXg18TezH24Iy0a+HC1uJz3t5yg9/zNjFySwvo956W7TQhxW/e2CmRY58YoCrz45V6b/72QoU5OzsNVy4guTXkkKpRNR7L4dPtpNh3JYtvxHLYdzyHAx53B4SHERTYmvIm/TGcthKjkpQfbsfnIJY5mFrB0awZ/6dXCZmWZlHCaN29u1h+s559/nkmTJpn8PWE8FxcNfe8Oou/dQZzLvcbKHadZmXqGzPwSlm07ybJtJwkL8GZIRGPiIkK4q4Gv2iELIRxAfR93XhzYlqmr97IoOZ1hkY1pWMfTJmWZlHCWLVtmViFhYWFmfc9Whg4dyubNm+nbty+rV69WOxyra1zXi4T72jCpbyu2pF9iXdp5fjiQycmcIv6VnM6/ktNp16gOA9oHM6BDEG2C/OTMR4ha7OHOTfhk+2n2nMllzobDLHg0wiblmJRwevXqZfj/Dz/8QEREBA0bNrR6ULb217/+lSeffJLly5erHYpNuWpd6NM2iD5tgygsKWfjwUzWpZ1jS3o2By/kc/BCPm//eJSwAG8GtA/mvvbBRITWResiyUeI2sTFRcPsIe2JS9rKl7vO8Xi3ZnRuWs/q5Zh9DScuLo6SkhKCgoKIiIggMjKS/v37Exsba8XwbCM2NrbWPTXBx8OVuMjGxEU25kphKT8eyuT7AxfZkp7NyZwi3ttygve2nKCutxv3tAykR8tAujWvW2kb27dv5+jRo7Ru3RqdTqdORYQQNhERWpfhUU1YvfMskz7cwvwBQXTr1s2qZZg9Su3q1avs2bOHt956iw4dOrB9+3buv/9+evfuTWGh7YblbtmyhcGDBxMSEoJGo2Ht2rVV1klKSiIsLAxPT090Oh07duywWTzOqJ6PO49Eh/Kf+C7sfrk/SSM7MyQ8BD8PV3KLyvhm7wWmf7mP2H/8wuu7tbzy9SEeT5xLTM/ejB49mm7dupGYmKh2NYQQVlaaugp9aTFnr7nSZ9Rfrb6fm32Gk5OTQ4cOHejQoQN//vOfAcjKymLYsGG8+uqrzJkzx2pB3qywsJDw8HCefPJJhg0bVuXzlStXkpCQwOLFi9HpdCxcuJABAwZw5MgRs7r/SkpKKCn537w0+fn5AJSVlVFWZp8bLW+UY4vy3F3gvrsDue/uQMor9Ow9l8+vx7LZevwye87kkVUMH28/A5oOhE76jLKL6ZSe3sO7qzYweMhD6Lpa9651W9bVkcu1JUdow7ai1u/NFtSuS2pqKv9ZtACf7iPxiR5KvZ6jWbRoGnFxcURHR9/2e6bEq1EUxawH37u4uBAcHEx4eLjhFRERQWFhIQ888ABZWVnmbNYkGo2GNWvWEBcXZ1im0+no0qUL77zzDgB6vZ7Q0FAmTpzI9OnTDett3ryZd955p9pBA7NmzWL27NlVlq9YscLoaRqc1bVySM/XcDhXw5E8DdnFla/teGkVWvkrtK2r0MZfIdA2A1tqtKKiIkaOHEleXh516tSxSRm1uQ0L0xWWwSu7tRRXaHiiVQWRgXdOEaa0YbMTzvHjx9mzZw979+5lz5497Nmzh5MnT+Lu7k5ZWRmPPfYYOp2OiIgIevToYU4R1fpjwiktLcXb25vVq1dXSkLx8fHk5uaybt06wzJjE86tjg5DQ0O5cOECAQEBVq3P7ZSVlbFx40b69++Pm5ubXcq8udzAwED69++PS52GuId2wr1ZOO5NOuDiWXloddP6XtzTIoB7WwbQrXl96niZHqvadbV3uTk5OTRq1MimCccR2rCtqPV7swW165Kamkrfvn0B8O46HN9uIyjPPs0XT4bf8fl7+fn5BAYGGtWGze5Sa9GiBS1atKjUrZWfn8+mTZsYOnQoiqKwfPlyEhMTKSoqMrcYk2RnZ1NRUUFQUFCl5UFBQRw+fNjwvl+/fuzZs4fCwkKaNGnCqlWriImJueU2PTw88PDwICkpiaSkJCoqKgBwc3Oze6NQo0yArl27MnHiRObNm0dh5ilIXc/UxEQeey6RX9Kz+TU9m12nr3D68jVOXz7LZ7+fxUVz/SJkj1YN6NEqkPDQurhpjb9kqFZd7V2uPcpypDZsK1IXy8XExBj28+KUL/HuPATXwKYoDVvfMR5TYjU74dSvX5+IiAhDd1rHjh3x9fXl22+/pUWLFnz66acAhsbtSH788UeTvzN+/HjGjx9Pfn4+/v7+NojKsc2dO5dhw4ZVGaUW2bQek/q2oqCknJTjOfySfolfjmVz4lIhu07nsut0Lv9MTsfPw5VuLQLo2SqQe1s1ICzAW+79sbPa3oZF9W7ez3+9Vo/vTxTzwa8Z3NMy0CrbNzvhLF261NCVtm7dOk6ePAmAt7c3X3zxhWE9rVZrcZDGCgwMRKvVkplZeX6HzMxMgoODLdr2H48OayOdTnfb4dC+Hq70axdEv3bXzy7P5V7j1/RLbEnPZuuxbHKLyth4MJON/517o0k9L3q0CqRHqwZ0bxFAXW93u9WjtpI2LIxxYz+/N7uQ7+dv5qfDWZy5XERofcuv91l0H87N10muXr3KhQsXaNy4MT4+PhYHZg53d3eioqJITk42xKbX60lOTmbChAkWbVuODk3TuK4XI7o0ZUSXplToFQ6cz+OX9Gx+Sb/EzlNXOHvlGp/tOMNnO87gooGOTerSs1Ug99xVjwqzriqK6kgbFqZoHujDPS0D2Hosh3Vp55jQp5XF27Tawzv9/Pzw8/Oz1uZuq6CggGPHjhneZ2RkkJaWRv369WnatCkJCQnEx8cTHR1N165dWbhwIYWFhYwZM8aicuXo0HxaFw2dmtSlU5O6jO/dksKScnZkXGZL+iV+Tc8mPauAPWdy2XMml0U/gZdWy48Fe+jdNoherRvY7LlOtY20YWGqhyIas/VYDmt2n2N875YWd4M73cM7U1NT6d27t+F9QkICcH0k2rJlyxgxYgSXLl1ixowZXLx4kYiICDZs2FBlIIGp5OjQenw8XOndtiG9216/L+pC3jV+Tc9mS3o2W45mkXetnG/3Z/Lt/uvdb+0a1SG2TQNi2zSkc9O6uJow+ED8j7RhYar7OwTz8tr9HL9UyKELV2kXYtlISqd7eGdsbCzVjeSeMGGCxV1own4a+XvxSHQoj0SHUlxSyuJV31EW2JpfjuWw92ye4blv/958HH8vN/re3ZD72wfTs3UDPN3sd41QiNqmjqcb97YMJPlwFlvSL9k34dz88M7aRroj7EProqG5HzzQtyVT77+b7IISthy9xOYjl9iSfoncojK+3HWOL3edw8tNS6/WDbi/QzB97m5IHc+aMSzWVqQNC3P0aHU94fySfsniuXJkAjYjSXeEOgJ9PRjWuQnDOjehQq+QevIyGw5c5IcDmZzLvcaGAxfZcOAiHq4u9GsXxLDIxvRs3cCke35qC2nDwhw9WzcA4PeMKxSXVVjUqyAJRzgNrYsG3V0B6O4KYMagduw/l8/3By7y3f4LHL9UyDd7L/DN3guGWU6HRzWhQ2P5wyqEJZoH+hDo6052QSkHL+RbNG2BHAYaKSkpiXbt2t3xEQ/CfjQaDR2b+DNlQBt+TOjF+gn3MuaeMAJ93ckpLGXZtpMMWvQrD7+7ja/2nKeswrZztTsDacPCHBqNxnDgtv9cnkXbkoRjpPHjx3Pw4EF+//13tUMRf3Aj+cwc3J6UF/vy4ZguDOrUCFcXDTtPXWHSZ7u5Z85P/Cs5nbwi53+qsLmkDQtzdfxvwtl31rKEI11qokZx1brQu01DerdpSFZ+MZ9uP82KHafJulrCgo1HeX/LCZ7oHsbTPZrL0w2EMFLroOv3WGZkWzbXmZzhiBqrYR1PXujfmq2JffjnnyJoG+xHQUk572w6Rq+3NvPh1gzpahPCCE3/+1ibM1csexCzJBwjSf+383J3deGhiMZ8O6kHix+Pok2QH3nXypi9/iAP/PMXdp++onaIdiFtWJjrRsLJzC+huMz8YfWScIwk/d/Oz8VFw/0dgvn2rz14Y2hHAnzcSc8q4OF3tzH/h/Qa/ww3acPCXHW93fDzuH4F5uyVa2ZvRxKOqHW0LhpG6pryY0IvhkY2Rq/Ae79kkHRAS3ZBSfUbEKKW0Wg0BPp5AHC5sNTs7UjCEbVWPR933h4RQdLIzvh4aDl+VcOj7+/gpIUXRoWoiep6X3+Sx5UiSTg2J/3fNdeDnRqx5i/dCPBQOHPlGiPe/42zFl4cdUTShoUl6v13VGeuJBzbk/7vmq15oA/Pd6igZQMfMvNLGL10h0VdB45I2rCwxP/OcMy/l00SjhD/VccdlsZHEeLvyYlLhUxZtafaJ5MLUVvU9bpxhiMJRwiraOTvydIxXXB3deGnw1l8uv202iEJ4RC83K+ni5JyGRYthNW0Da5D4v1tAZi74bBFfdZC1BQertefEl1Sbv7N0pJwhLiFJ7qH0TbYj6vF5fx783G1wxFCdR6u/z3DKZOEI4RVaV00hrOcT1NOUVBSrnJEQqjLkHCkS832ZEhp7RPbpgF3NfChsLSCdWnn1A7HYtKGhSU83KRLzW5kSGnto9Fo+FOXUAC+23dR5WgsJ21YWOJ/ZziScISwif7tggHYnpHD1eLaO5eOEDcGDcjDO4WwkeaBPjQP9KGsQiHlxGW1wxFCNVoXDQB6vfn3pknCEaIaXcKuz+G+92yuuoEI4eQk4QhRjY5N6gKwx8LpdYWo7Wplwvn6669p06YNrVq14j//+Y/a4QgHd3fw9el1j2cVqByJEM7NVe0A7K28vJyEhAQ2bdqEv78/UVFRDB06lICAALVDEw6qacD12Q4v5F2jtFyPu2utPE4TwmK1bs/ZsWMH7du3p3Hjxvj6+jJw4EB++OEHtcMSDqyBrwdeblr0CpzPNX+2QyFqO6dLOFu2bGHw4MGEhISg0WhYu3ZtlXWSkpIICwvD09MTnU7Hjh07DJ+dP3+exo0bG943btyYc+ec/6Y+YTsajYZ61yc7ZNP23eoGI4TKsrKy2L59u1nfdbqEU1hYSHh4OElJSbf8fOXKlSQkJDBz5kx27dpFeHg4AwYMICsry86RipoiMTGRE4f3ATA+YRqJiYkqRySE/X388UcAHE1Pp1u3bmbtB053DWfgwIEMHDjwtp8vWLCAsWPHMmbMGAAWL17MN998w9KlS5k+fTohISGVzmjOnTtH165db7u9kpISSkr+N899fn4+AGVlZZSV2edGwBvl2Ks8Nct1tLqmpqayaNEi/AdNA8CzTiCLFi0iLi6O6Ohoq5VrS47Qhm1FrfZiC45cl9TUVL755hvqDmqPi4sLXl5ehv2gdevWRm9HozjxDFMajYY1a9YQFxcHQGlpKd7e3qxevdqwDCA+Pp7c3FzWrVtHeXk5d999N5s3bzYMGti2bdttBw3MmjWL2bNnV1m+YsUKvL29bVEt4YBWHHNh+yUXBjWtoH9j6+0yRUVFjBw5kry8POrUqWO17d5M2rCwhr2XNXxwREtzP4XnO/zvaQOmtOEalXBuXJ/Ztm0bMTExhvWmTZvGzz//bOh3/Oqrr5gyZQp6vZ5p06bxzDPP3LaMWx0dhoaGcuHCBbuNbCsrK2Pjxo30798fNzc3u5SpVrmOVtfU1FT69u2L772j8O48hMKd6yjc+gnJyclWOcPJycmhUaNGNk04jtCGbUWt9mILjlyX1NRUBo6dTt1B0yg9f5jc1S8DkJycTOvWrQkMDDSqDTtdl5o1DBkyhCFDhhi1roeHBx4eHiQlJZGUlERFxfXM7ubmZvdGoUaZapXrKHWNiYlh4sSJvPfbebyBClyYNGlSpQMaS8uzNUdqw7YidbGtmJgYHnzwQbYCer2ea9eukZiYSExMjKGL1hhON2jgTgIDA9FqtWRmZlZanpmZSXBwsEXblift1l5z587l2aefBCBu2HDmzJmjckTmkTYsLDF69GgAWrdqRUpKiln7QY1KOO7u7kRFRZGcnGxYptfrSU5OtviIVOYSqd1ahjUDwL9+oMqRmE/asLCGwAYN0Ol0Zn3X6RJOQUEBaWlppKWlAZCRkUFaWhqnT58GICEhgSVLlrB8+XIOHTrEuHHjKCwsNIxaM5ccHdZuHm6Wz3aoNmnDwhIaK2zD6a7hpKam0rt3b8P7hIQE4PpItGXLljFixAguXbrEjBkzuHjxIhEREWzYsIGgoCCLyv1j/7eoXawxn7vapA0LtTldwomNjaW6gXUTJkxgwoQJVi13/PjxjB8/nvz8fPz9/a26beH4bkw+Zclsh2qTNizU5nRdamqR/u/a7cYZjiWzHapN2rBQmyQcI0n/d+3mcmO2Q+e9bU3asFCdJBwhhBBGs+SQSxKOkaQ7Qjg7acNCbZJwjCTdEcLZSRsWltBoLB8YLQlHCCGEXUjCEUIIYReScIwk/d8CLLtgqjZpw0JtknCMJP3ftZs1HuuhNmnDQm2ScIQQQhjNklvRJOEIIYSwC0k4QgghqmWNbmVJOEaSC67C2UkbFmqThGMkueAqnJ20YaE2SThCGOHGXdZO/OxOIVQnCUcIIYTR5OGdQgghHJ4kHCGEEHYhCUcIIUS1rPCwaEk4xpIhpcLZSRsWapOEYyQZUlq73Ti4c+ZBatKGhdok4QghhLALSThCCCGMZ8HNaJJwhBBC2IUkHCGEEHZRKxPO0KFDqVevHsOHD1c7FCGEcAoyLNpMf/3rX/noo4/UDkM4EWvsbELUdrUy4cTGxuLn56d2GMIJZWdns337drXDEEI12Tk5Zu8DDpdwtmzZwuDBgwkJCUGj0bB27doq6yQlJREWFoanpyc6nY4dO3bYP1BRqyxbtgyAE8dP0K1bNxITE9UNSAg7+/DDZQCcOGH+PuBq5ZgsVlhYSHh4OE8++STDhg2r8vnKlStJSEhg8eLF6HQ6Fi5cyIABAzhy5AgNGzYEICIigvLy8irf/eGHHwgJCTEpnpKSEkpKSgzv8/PzASgrK6OsrMykbZnrRjn2Kk/Nch2xrqmpqaxfv566Qzrg4uKCl5cXixYtIi4ujujoaKuUa0uO0IZtRa32YguOXJc77QOtW7c2ejsaRXHcGT40Gg1r1qwhLi7OsEyn09GlSxfeeecdAPR6PaGhoUycOJHp06cbve3NmzfzzjvvsHr16juuN2vWLGbPnl1l+YoVK/D29ja6POHc9l/RsOSwlqY+CpM7VVhtu0VFRYwcOZK8vDzq1Kljte3eTNqwsIYDVzS8f1hLqI/ClJv2AVPasFMlnNLSUry9vVm9enWlJBQfH09ubi7r1q0zetvGJpxbHR2GhoZy4cIFAgICTKqPucrKyti4cSP9+/fHzc3NLmWqVa4j1jU1NZWBT02l7pAXKbt4jCtfvAhAcnKyxWc4OTk5NGrUyKYJxxHasK2o1V5swZHrUmkfyDzGlZX/2wdat25NYGCgUW3Y4brU7iQ7O5uKigqCgoIqLQ8KCuLw4cNGb6dfv37s2bOHwsJCmjRpwqpVq4iJibnluh4eHnh4eJCUlERSUhIVFdczu5ubm90bhRplqlWuI9U1JiaGwYMH8wvXz6ivXbtGYmLibduMqeXZmiO1YVuRutjW9X1g0C33gRtdtMZwqoRjLT/++KPJ3xk/fjzjx48nPz8ff39/G0QlHNmYMU/wy7JU7mpxF2tTUtDpdGqHZDJpw8IST455kl+W/c5dzc3fBxxulNqdBAYGotVqyczMrLQ8MzOT4OBgm5Ytj3YXAAEBAU6ZbEDasLCOgEDz9wGnSjju7u5ERUWRnJxsWKbX60lOTrZK98adyKPdhbOTNiyswZKr/g7XpVZQUMCxY8cM7zMyMkhLS6N+/fo0bdqUhIQE4uPjiY6OpmvXrixcuJDCwkLGjBlj07j+2P8thLORNizU5nAJJzU1ld69exveJyQkANdHoi1btowRI0Zw6dIlZsyYwcWLF4mIiGDDhg1VBhJYm/R/C2cnbViozeESTmxsLNWN1J4wYQITJkywU0RCCCGswamu4ahJLrjWbhqc/+md0oaFReRp0fYjF1wFWHbBVG3ShoXaJOEIIYQwmoJMMW1z0h0hnJ20YaE2SThGku4I4eykDQu1ScIRQghhF5JwhDCG8w9SE8Ii1tgFJOEYSfq/hbOTNizUJgnHSNL/LcCyETpqkzYs1CYJRwghhNEsuRdNEo4QQgi7kIRjJOn/Fs5O2rBQmyQcI0n/d+1WEwapSRsWapOEI4QQoloajeWHXZJwhDCBMz+8Uwi1ScIRQghhNBmlJoQQwuFJwhFCCGEXknCMJENKhbOTNizUJgnHSDKktHazxggdtUkbFmqThCOEEKJa8rRoIexMhkULYT5JOEIIIYxmyTGXJBwhhBB2UesSzpkzZ4iNjaVdu3Z06tSJVatWqR2SEELUCq5qB2Bvrq6uLFy4kIiICC5evEhUVBQPPPAAPj4+aocmHJjzj1ETQn21LuE0atSIRo0aARAcHExgYCCXL1+WhCOEEDbmcF1qW7ZsYfDgwYSEhKDRaFi7dm2VdZKSkggLC8PT0xOdTseOHTvMKmvnzp1UVFQQGhpqYdSitrh85TLbt29XOwwh7O7GrWhXLNgHHC7hFBYWEh4eTlJS0i0/X7lyJQkJCcycOZNdu3YRHh7OgAEDyMrKMqwTERFBhw4dqrzOnz9vWOfy5cuMHj2a999/3+Z1Es7vgw8+AOD0qdN069aNxMRElSMSwr7+85//AHD69Bmz9wGH61IbOHAgAwcOvO3nCxYsYOzYsYwZMwaAxYsX880337B06VKmT58OQFpa2h3LKCkpIS4ujunTp9O9e/dq1y0pKTG8z8/PB6CsrIyysjJjqmSxG+XYqzw1y3XEuqamprJu3VrqxnVE4+KCl5cXixYtIi4ujujoaKuUa0uO0IZtRa32YguOXJfU1FTWrl1HvaGd0Gg0lfaB1q1bG70djaI47q1sGo2GNWvWEBcXB0BpaSne3t6sXr3asAwgPj6e3Nxc1q1bV+02FUVh5MiRtGnThlmzZlW7/qxZs5g9e3aV5StWrMDb29vYqggndzhXw7uHtDT2VpgWXmG17RYVFTFy5Ejy8vKoU6eO1bZ7M2nDwhpu7AMh3gqJN+0DprRhp0o458+fp3Hjxmzbto2YmBjDetOmTePnn382ql/x119/pWfPnnTq1Mmw7OOPP6Zjx463XP9WR4ehoaFcuHCBgIAAM2tmmrKyMjZu3Ej//v1xc3OzS5lqleuIdU1NTeX+J16g3tCXKbt0kiufTQUgOTnZ4jOcnJwcGjVqZNOE4wht2FbUai+24Mh1udM+0Lp1awIDA41qww7XpWZr9957L3q93uj1PTw88PDwICkpiaSkJCoqrmd2Nzc3uzcKNcpUq1xHqmtMTAxxcQ/xM6Do9Vy7do3ExMRKBz2WlGdrjtSGbUXqYluV9gFFqbQP3OiiNYbDDRq4k8DAQLRaLZmZmZWWZ2ZmEhwcbNOy5Um7tdvTTz8NQNNmTUlJSWHOnDkqR2Q6acPCEjf2gdDQULP3AadKOO7u7kRFRZGcnGxYptfrSU5OtsrR5p3IXCICoG7deuh0OrXDMIu0YWEJzX9vf65Xz/x9wOESTkFBAWlpaYaRZhkZGaSlpXH69GkAEhISWLJkCcuXL+fQoUOMGzeOwsJCw6g1W5GjQ+HspA0LtTncNZzU1FR69+5teJ+QkABcH4m2bNkyRowYwaVLl5gxYwYXL14kIiKCDRs2EBQUZNO4/tj/LYSzkTYs1OZwCSc2NpbqBs5NmDCBCRMm2Cmi68aPH8/48ePJz8/H39/frmUL9WlqwNPUpA0LtTlcl5qjkv5v4eykDQu1ScIxkvR/C2cnbVioTRKOEEIIu5CEYyTpjhDOTtqwsMSNp0Vb8mwaSThGku4I4eykDQu1ScIRQghhF5JwhDCCxvlHRQuhOkk4RpL+b+HspA0LtUnCMZL0fwtnJ21YqE0SjhBCCLuQhCOECRx3ukIhbOvGZUwF83cCSThCCCHsQhKOkeSCa+1WEwapSRsWapOEYyS54CqcnbRhoTZJOEIIIexCEo4QQgi7kIQjhAksGaEjRE0gD+8UQghhW1YYOSMJRwhj1IRhakKoTBKOkWRIqXB20oaF2iThGEmGlApnJ21YqE0SjhBCCLuQhCOEEMIuJOEIYQJ5eKeo7SzZBSThCGEEjQxTE7WcNfaBWpdwcnNziY6OJiIigg4dOrBkyRK1QxJCiFrBVe0A7M3Pz48tW7bg7e1NYWEhHTp0YNiwYQQEBKgdmhBC1Gi17gxHq9Xi7e0NQElJCYqioEjHvDBSbl4u27dvVzsMIVSTl5dn9j7gcAlny5YtDB48mJCQEDQaDWvXrq2yTlJSEmFhYXh6eqLT6dixY4dJZeTm5hIeHk6TJk2YOnUqgYGBVope1FTvv/8eAOfOnaNbt24kJiaqHJEQ9mWNfcDhutQKCwsJDw/nySefZNiwYVU+X7lyJQkJCSxevBidTsfChQsZMGAAR44coWHDhgBERERQXl5e5bs//PADISEh1K1blz179pCZmcmwYcMYPnw4QUFBt4ynpKSEkpISw/u8vDwALl++bI3qGqWsrIyioiJycnJwc3Or0eU6Yl13797N/61eRd0hraC8FE9PT/71r3/Rv39/IiMjLSr3Rjuy5Vm2I7RhW1GrvdiCI9el8j5QUmkfaNmyJWBkG1YcGKCsWbOm0rKuXbsq48ePN7yvqKhQQkJClDfffNOsMsaNG6esWrXqtp/PnDlT4fpIQHnJy2av48ePm9V+jSFtWF72eJ05c6batqhRFMe9gKHRaFizZg1xcXEAlJaW4u3tzerVqw3LAOLj48nNzWXdunXVbjMzMxNvb2/8/PzIy8vjnnvu4bPPPqNjx463XP+PR4e5ubk0a9aM06dP4+/vb1H9jJWfn09oaChnzpyhTp06dilTrXJrU13h+tlG06ZNuXLlCnXr1rVJGY7Qhm1Frd+bLThrXRRF4erVq4SEhODicuerNA7XpXYn2dnZVFRUVOn+CgoK4vDhw0Zt49SpUzzzzDOGwQITJ068bbIB8PDwwMPDo8pyf39/uzeKOnXqqNIQ1Si3NtUVqHZHtYQjtWFbUev3ZgvOWBdjD1ycKuFYQ9euXUlLS1M7DCGEqHUcbpTanQQGBqLVasnMzKy0PDMzk+DgYJWiEkIIYQynSjju7u5ERUWRnJxsWKbX60lOTiYmJsYuMXh4eDBz5sxbdlHUpDLVKrc21VWtctWqqy1IXZyLww0aKCgo4NixYwBERkayYMECevfuTf369WnatCkrV64kPj6e9957j65du7Jw4UK++OILDh8+fNuhzUIIIdTncAln8+bN9O7du8ry+Ph4li1bBsA777zDW2+9xcWLF4mIiOBf//oXOp3OzpEKIYQwhcMlHCGEEDWTU13DEUII4bwk4QghhLALSThCCCHsQhLOLZj6NOpVq1bRtm1bPD096dixI99++61Nyzxw4AAPP/wwYWFhaDQaFi5caHJ55pS7ZMkSevToQb169ahXrx79+vUz+Undppb55ZdfEh0dTd26dfHx8SEiIoKPP/7Y5DJNLfdmn3/+ORqNptLjlGxR5rJly9BoNJVenp6eJpdpieomKBw6dCj16tVj+PDhdo3LXHeKd/78+bRv354OHTrwySefqBCdeY4cOUJERITh5eXldcun6jskWz0w0Fl9/vnniru7u7J06VLlwIEDytixY5W6desqmZmZt1x/69atilarVebNm6ccPHhQeemllxQ3Nzdl3759Nitzx44dypQpU5TPPvtMCQ4OVt5++21zqmpyuSNHjlSSkpKU3bt3K4cOHVKeeOIJxd/fXzl79qzNyty0aZPy5ZdfKgcPHlSOHTumLFy4UNFqtcqGDRtsWtcbMjIylMaNGys9evRQHnroIZuW+eGHHyp16tRRLly4YHhdvHjRpDItVV5erhQWFiqKoigFBQVKWFiYkp2dbfh806ZNyldffaU8/PDDdo3LXLeLd+/evUpkZKRy7do1paioSNHpdMqVK1fUCdICV69eVQICApSCggK1QzGKJJw/MPVp1I8++qjy4IMPVlqm0+mUZ5991mZl3qxZs2ZmJxxLn7xdXl6u+Pn5KcuXL7dbmYqiKJGRkcpLL71k9PrmllteXq50795d+c9//qPEx8ebnHBMLfPDDz9U/P39TSrDlnJycpRmzZoply5dqrR806ZNTpNwFOXW8a5cuVJ57rnnDO+feeYZ5bPPPrN3aBb79NNPlUcffVTtMIwmXWo3KS0tZefOnfTr18+wzMXFhX79+vHbb7/d8ju//fZbpfUBBgwYcNv1rVGmNVij3KKiIsrKyqhfv75dylQUheTkZI4cOULPnj2NKtOScl955RUaNmzIU089ZXRZlpZZUFBAs2bNCA0N5aGHHuLAgQOVPq9JExTaoy6306FDBzZv3kxubi5Xrlxh8+bNnDt3zirbtme9vvjiC0aMGGFhxPYjCecmd3oa9cWLF2/5nYsXL5q0vjXKtAZrlJuYmEhISEiVhGvtMvPy8vD19cXd3Z0HH3yQRYsW0b9/f6PKNLfcX3/9lQ8++KDKNQxbltmmTRuWLl3KunXr+OSTT9Dr9XTv3p2zZ88a1rkxQWFSUtItt3FjgsKZM2eya9cuwsPDGTBgAFlZWYZ1blyf+ePr/PnzAIYJCjMyMlixYkWVZxdaiz3qcjvt2rVj0qRJ9OnTh2HDhtGtWze0Wq1T1Ss/P59t27bxwAMPWCVue6h1T4sW1jFnzhw+//xzNm/ebPML235+fqSlpVFQUEBycjIJCQncddddxMbG2qS8q1evMmrUKJYsWWLX6cdjYmIqPROwe/fu3H333bz33nu8+uqrAAwcOJCBAwfedhsLFixg7NixjBkzBoDFixfzzTffsHTpUqZPnw5g9NPSg4KCCA8P55dffrHJIAF71uVWnn32WZ599lkAnn76aVq1amX2tm5mr3qtW7eO++67z+4DSywhZzg3Medp1MHBwRY9vVqtJ2BbUu78+fOZM2cOP/zwA506dbJ5mS4uLrRs2ZKIiAgmT57M8OHDefPNN21W7vHjxzl58iSDBw/G1dUVV1dXPvroI7766itcXV05fvy41cu8FTc3NyIjIw3PFqyONbpJMzMzuXr1KnD9zHLLli20adPGqO9akz26mm+cURw5coQdO3YwYMAAq2z3TqxZL2frTgNJOJWY8zTqmJiYSusDbNy40einV6v1BGxzy503bx6vvvoqGzZsIDo62i5l/pFer680g6W1y23bti379u0jLS3N8BoyZAi9e/cmLS2N0NBQq5d5KxUVFezbt49GjRoZtb41uklPnTpFjx49CA8Pp0ePHlUmKOzXrx+PPPII3377LU2aNLHZdUZrdTXfKd6HHnqIdu3a8fjjj/Phhx/i6mr7Dh9r1SsvL89uSdKq1B614Gg+//xzxcPDQ1m2bJly8OBB5ZlnnlHq1q1rGJ46atQoZfr06Yb1t27dqri6uirz589XDh06pMycOdOsYdGmlFlSUqLs3r1b2b17t9KoUSNlypQpyu7du5X09HSb1nXOnDmKu7u7snr16kpDd69evWqzMt944w3lhx9+UI4fP64cPHhQmT9/vuLq6qosWbLEpnX9I3NGqZla5uzZs5Xvv/9eOX78uLJz507lT3/6k+Lp6akcOHDgltsHlDVr1hjenzt3TgGUbdu2VVpv6tSpSteuXU2K3d5qUl1uVlPrZS65hvMHI0aM4NKlS8yYMcPwNOoNGzYYjkhOnz5daTrg7t27s2LFCl566SX+9re/0apVK9auXUuHDh1sVub58+eJjIw0vJ8/fz7z58+nV69ebN682Wblvvvuu5SWllbpz585cyazZs2ySZmFhYU899xznD17Fi8vL9q2bcsnn3xicleCqeVag6llXrlyhbFjx3Lx4kXq1atHVFQU27Zto127dkaVV5MmKKxJdblZTa2X0dTOeEII8/CHo2dFuX7vz4QJEwzvKyoqlMaNG5t0n5MaalJdblZT62UuOcMRwoncPEEhQEZGBmlpaYYJChMSEoiPjyc6OtowQWFhYaFhRJQjqUl1uVlNrZdVqJ3xhBDG27RpkwJUecXHxxvWWbRokdK0aVPF3d1d6dq1q5KSkqJewHdQk+pys5paL2uQCdiEEELYhQyLFkIIYReScIQQQtiFJBwhhBB2IQlHCCGEXUjCEUIIYReScIQQQtiFJBwhhBB2IQlHCCGEXUjCEUIIYReScIQQQtiFJBxhNEVRWLBgAc2bN8fb25u4uDjy8vJuuW5sbCwajQaNRlPtdLlPPPGEYd21a9daP3AhhEOQhCOMNnXqVN59912WL1/OL7/8ws6dO+84D87YsWO5cOFCtXMD/fOf/+TChQtWjlaIykw5CLKX2NhYnn/+ecP7mn7wJQlHGGX79u0sWLCAlStX0rNnT6Kiohg7dizffvvtbb/j7e1NcHBwtVP3+vv7147Jp4TZevXqxZNPPllp2cKFC/Hx8eHdd981ejvGHgRZy5gxY3jppZeMXr+mH3zJfDjCKPPnz6dv37507tzZsCwoKIjs7GwVoxK1gaIo7N69m0ceeQSAoqIixo4dy6ZNm9i4cSPdu3c3els3DoLsoaKigq+//ppvvvnG6O/4+/vj7+9vw6jUJWc4ololJSV88803DB06tNLy4uLiGr1zCMeQnp7O1atX6dy5MxkZGXTv3p2MjAx27txpUrK5Hb1ez7x582jZsiUeHh40bdqU119/3fB5bGwsEydO5Pnnn6devXoEBQWxZMkSw6Rpfn5+tGzZku+++67Sdrdt24abmxtdunQBrk+XPnr0aHx9fWnUqBH/+Mc/LI7d2UjCEdXatWsX165dY/Lkyfj6+hpe06ZNo3Xr1mqHJ2q4nTt3otVqyczMJDo6Gp1Ox+bNm2nUqJFVtv/iiy8yZ84cXn75ZQ4ePMiKFSsICgqqtM7y5csJDAxkx44dTJw4kXHjxvHII4/QvXt3du3axX333ceoUaMoKioyfOerr75i8ODBaDQa4Po10J9//pl169bxww8/sHnzZnbt2mWVOjgNded/E85g2bJlio+Pj5Kenl7p1bp1a+WVV1655Xd69eql/PWvf6207JNPPlF8fHwMry1btlT6nFvM/y7ElClTFK1Wq7i4uChJSUlmb+dWbTI/P1/x8PBQlixZcsfv3XvvvYb35eXlio+PjzJq1CjDsgsXLiiA8ttvvxmWtWrVSvn6668VRVGUq1evKu7u7soXX3xh+DwnJ0fx8vKqEpOi1Nx9Qa7hiGrl5+cTGBhIy5YtDctOnTpFeno6Dz/8sNHbGTJkCDqdzvC+cePGVo1T1Ey7du2iX79+7N+/n507d1p124cOHaKkpIS+ffvecb1OnToZ/q/VagkICKBjx46GZTfOiLKysgzbPX/+vGG7x48fp7S0tFL7r1+/Pm3atLFaXZyBdKmJagUGBpKXl4dy02zkr7/+Og888ADt2rUzejs3+rpvvLy8vGwRrqhhdu3axcCBA1m3bh2fffYZb731VpV1lixZQufOnenQoQMjRowwetvGtkE3N7dK7zUaTaVlN7rN9Ho9cL07rX///nh6ehodS20gCUdUq0+fPhQXFzNnzhwyMjJ47bXXWL9+vUnDUYUwx4kTJ8jNzaVz585ERUXx4Ycf8uKLL7Ju3TrDOleuXCEpKYnff/+d/fv389577xm9/VatWuHl5UVycrJV4163bh0PPfSQ4X2LFi1wc3Nj+/btleI+evSoVct1dJJwRLWCgoJYtmwZ7777Lu3btyclJYVff/2V0NBQtUMTNdzOnTvRaDREREQAMGLECP72t7/x5z//2XDzpqurK1euXGHatGkcOHCAunXrGr19T09PEhMTmTZtGh999BHHjx8nJSWFDz74wOyYs7KySE1NZdCgQYZlvr6+PPXUU0ydOpWffvqJ/fv388QTT+DiUrv+BMs1HGGUESNGmNRVIYQ17Nq1i1atWuHn52dYNnv2bA4ePMiQIUPYsWMHwcHB7N+/n7Vr1/Loo4/y+uuvExcXZ3QZL7/8Mq6ursyYMYPz58/TqFEj/vKXv5gd8/r16+natSuBgYGVlr/11lsUFBQwePBg/Pz8mDx58m0fDVVjqT1qQdRMvXr1Utzc3BQfHx9l7969d1z32WefVXx8fGrsyBxhW0ePHjX8f9y4ccrnn39+y/VuNUrNFgYPHqzMnTvXom3U1H2hdp3PCbv59NNPOXjwIGlpadWOxHnllVdIS0sjPT2d/v372ylCUVO89tprtGnThsjISDQajeGJBLfy73//G19fX/bt22ezeO69914ee+wxs777l7/8BV9fXytH5Dg0inLT0CMhhKihzp07x7Vr1wBo2rQp7u7uKkdUVVZWFvn5+QA0atQIHx8flSOyLkk4Qggh7EK61IQQQtiFJBwhhBB2IQlHCCGEXUjCEUIIYReScIQQQtiFJBwhhBB2IQlHCCGEXUjCEUIIYReScIQQQtjF/wOwEDzH/CYG3gAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_compare(soilsample, soilm_panday)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fit method finds the optimal curve through both the soil water retention curve and hydraulic conductivity function at the same time using the least squares algorithm. All parameters are subject to the optimization algorithm.\n", "\n", "The fit method uses the optimization algorithm from the 1991 [RETC Code for Quantifying the Hydraulic Functions of Unsaturated Soils](https://www.pc-progress.com/Documents/programs/retc.pdf) by M.Th. van Genuchten, F.J. Leij and S.R. Yates.\n", "\n", "The objective function $O(b)$ minimized is:\n", "\n", "$ O(b) = \\sum^N_{i=1}(w_i(\\theta_{o,i}-\\theta_i))^2 + \\sum^M_{i=N+1}(w_iW_1W_2(Y_{o,i}-Y_i))^2$\n", "\n", "Using the SciPy least-squares algorithm. \n", "\n", "With $N$ the number of $\\theta$ data points and $M$ is the number of $K$ and $\\theta$ data points\n", "\n", "$w_i$ are the individual weight factors per measurements (by default 1 for each measurement)\n", "\n", "$W_1$ is the weight factor for the hydraulic conductivity function with respect to the soil water retention (default is 0.1)\n", "\n", "$W_2$ is the proportional weight factor for two different data types and elimination factor for different units. By default the formulation for \n", "$W_2 = \\frac{(M-N)\\sum^N_{i=1}w_i\\theta_{o,i}}{N\\sum^M_{i=N+1}w_i|Y_{o,i}|}$.\n", "\n", "$Y$ is indicates the for the logaritmic transform of such that $Y = log_{10}K$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It can be favorable to only optimize the relative hydraulic conductivity function, and leave parameter k_s untouched. That kan be achieved by providing `k_s` to the fit method." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Panday(k_s=np.float64(712.7999894094452), alpha=np.float64(0.14657623411049972), beta=np.float64(2.6437587979256927), brook=np.float64(3.8328536708431717), sr=np.float64(0.10448683153152591))" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "soilm_panday = soilsample.fit(pe.Panday, k_s=max(soilsample.k))\n", "soilm_panday" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is also possible to provide bounds for the parameter space. By default, the bounds argument is `None` which takes the stored parameter bounds per soil model:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
p_inip_minp_max
k_s50.000.00100100000.0
theta_r0.020.000010.2
theta_s0.400.200000.5
alpha0.020.001000.3
beta2.301.0000012.0
brook10.001.0000050.0
\n", "
" ], "text/plain": [ " p_ini p_min p_max\n", "k_s 50.00 0.00100 100000.0\n", "theta_r 0.02 0.00001 0.2\n", "theta_s 0.40 0.20000 0.5\n", "alpha 0.02 0.00100 0.3\n", "beta 2.30 1.00000 12.0\n", "brook 10.00 1.00000 50.0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "panday_bounds = pe.get_params(\"Panday\")\n", "panday_bounds" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Panday(k_s=np.float64(650.0024881306313), alpha=np.float64(0.14641136203923166), beta=np.float64(2.64744674802969), brook=np.float64(3.8182839905828785), sr=np.float64(0.10450412521086855))" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "panday_bounds.loc[\"k_s\", \"p_min\"] = 650\n", "panday_bounds.loc[\"k_s\", \"p_ini\"] = max(soilsample.k)\n", "soilm_panday = soilsample.fit(pe.Panday, pbounds=panday_bounds)\n", "soilm_panday" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Other available option are to print the optimization result." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SciPy Optimization Result\n", " message: `gtol` termination condition is satisfied.\n", " success: True\n", " status: 1\n", " fun: [ 2.220e-05 2.220e-05 ... -8.042e-05 -2.758e-04]\n", " x: [ 5.717e+02 4.495e-02 4.300e-01 1.462e-01 2.653e+00\n", " 3.798e+00]\n", " cost: 5.556972728926372e-07\n", " jac: [[ 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00]\n", " [ 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00]\n", " ...\n", " [ 1.946e-06 0.000e+00 ... -4.053e-02 -1.764e-02]\n", " [ 1.946e-06 0.000e+00 ... -5.026e-02 -2.187e-02]]\n", " grad: [ 1.630e-14 3.304e-09 -3.688e-10 5.776e-09 1.076e-09\n", " 5.679e-11]\n", " optimality: 1.777636707528239e-09\n", " active_mask: [0 0 0 0 0 0]\n", " nfev: 10\n", " njev: 10\n" ] } ], "source": [ "soilm_panday = soilsample.fit(pe.Panday, silent=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To parse other parameters for W1, W2, and individual weights there are other keyword arguements that for the fit method." ] } ], "metadata": { "kernelspec": { "display_name": "pedon (3.13.5)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.5" } }, "nbformat": 4, "nbformat_minor": 2 }