diff --git a/examples/example-parametric.ipynb b/examples/example-parametric.ipynb new file mode 100644 index 0000000..b874ad3 --- /dev/null +++ b/examples/example-parametric.ipynb @@ -0,0 +1,244 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "91dcd026-de5f-4888-b668-f4c88ae3d7ac", + "metadata": {}, + "source": [ + "# Example: Creating parametric families #1" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "54cdd27a-7183-4062-960b-cd28cd2bc521", + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from pysatl_core.distributions import (\n", + " DefaultComputationStrategy,\n", + " DefaultSamplingUnivariateStrategy,\n", + ")\n", + "from pysatl_core.families import (\n", + " ParametricFamily,\n", + " ParametricFamilyRegister,\n", + " constraint,\n", + " parametrization,\n", + ")\n", + "from pysatl_core.types import UnivariateContinuous\n", + "\n", + "PDF = \"pdf\"\n", + "CDF = \"cdf\"\n", + "PPF = \"ppf\"" + ] + }, + { + "cell_type": "markdown", + "id": "4854203d-24fe-4abf-ac3d-773dbbdab56b", + "metadata": {}, + "source": [ + "Let's create lognormal family of random variables, by providing PDF" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "73ccf204-fba8-4087-9f21-bdf59f04dade", + "metadata": {}, + "outputs": [], + "source": [ + "def lognormal_pdf(parameters, x: float) -> float:\n", + " if x <= 0:\n", + " return 0.0\n", + " return (\n", + " 1.0\n", + " / (x * parameters.sigma * math.sqrt(2.0 * math.pi))\n", + " * math.exp(-((math.log(x) - parameters.mu) ** 2) / (2 * parameters.sigma**2))\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "f33585a5-02b0-4a2d-89e2-910f61dd63d8", + "metadata": {}, + "source": [ + "Now let's create an object, that will represent our family" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "680acb78-fde1-45d8-84fc-50d3058436c8", + "metadata": {}, + "outputs": [], + "source": [ + "Lognormal = ParametricFamily(\n", + " name=\"Lognormal Family\",\n", + " distr_type=UnivariateContinuous,\n", + " distr_parametrizations=[\"canonical\", \"meanvar\"],\n", + " distr_characteristics={PDF: lognormal_pdf},\n", + " sampling_strategy=DefaultSamplingUnivariateStrategy(),\n", + " computation_strategy=DefaultComputationStrategy(),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "6709fd2b-45c6-4f66-8b0f-c9e5cd31bfad", + "metadata": {}, + "source": [ + "We specified that there will be two parametrizations: canonical (which will be treat as base) and mean-var. Let's introduce them" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "5a7da55b-4396-4cf8-bc87-328e7590e0d7", + "metadata": {}, + "outputs": [], + "source": [ + "@parametrization(family=Lognormal, name=\"canonical\")\n", + "class NormalParametrization:\n", + " mu: float\n", + " sigma: float\n", + "\n", + " @constraint(description=\"sigma > 0\")\n", + " def check_sigma_positive(self) -> bool:\n", + " return self.sigma > 0\n", + "\n", + "\n", + "@parametrization(family=Lognormal, name=\"meanvar\")\n", + "class MeanVarParametrization:\n", + " mean: float\n", + " var: float\n", + "\n", + " @constraint(description=\"mean > 0\")\n", + " def check_mean_positive(self) -> bool:\n", + " return self.mean > 0\n", + "\n", + " @constraint(description=\"var > 0\")\n", + " def check_var_positive(self) -> bool:\n", + " return self.var > 0\n", + "\n", + " def transform_to_base_parametrization(self) -> NormalParametrization:\n", + " mu = math.log(self.mean**2 / math.sqrt(self.mean**2 + self.var))\n", + " sigma = math.sqrt(math.log(1 + self.var / self.mean**2))\n", + " return NormalParametrization(mu=mu, sigma=sigma)" + ] + }, + { + "cell_type": "markdown", + "id": "0f11c1b8-4601-4c1c-a087-50eeb3134a89", + "metadata": {}, + "source": [ + "Note that for second parametrization we provided way convert it to base one. Now, let's register our family and do some things:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "60b6b5f2-0f6d-48f9-a45c-d927a7cc289e", + "metadata": {}, + "outputs": [], + "source": [ + "ParametricFamilyRegister.register(Lognormal)\n", + "dist = Lognormal(mean=1.0, var=1.0, parametrization_name=\"meanvar\")" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "2b631a99-42c3-4708-ad79-663786283814", + "metadata": {}, + "outputs": [], + "source": [ + "cdf = dist.computation_strategy.query_method(CDF, dist)\n", + "pdf = dist.computation_strategy.query_method(PDF, dist)\n", + "ppf = dist.computation_strategy.query_method(PPF, dist)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "cc1702cf-fda3-4897-b918-21515aff2a7b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGhCAYAAACzurT/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAaxpJREFUeJzt3XtclGXeP/DPDDAzHAcQGE4j4PmEoCiEh9QNpTIfbbcNrVVzcw896mZUm7alWb/CtpOVlptPZbVraidr0zRF0VTURE1RQZGjHAYQYYbTADPX7w90agKUQWCG4fN+xQvn5rpnvlc3MB/u+7quWyKEECAiIiKyEqm1CyAiIqLejWGEiIiIrIphhIiIiKyKYYSIiIisimGEiIiIrIphhIiIiKyKYYSIiIisimGEiIiIrIphhIiIiKyKYYSIiIisyuIwcuDAAcyYMQOBgYGQSCTYtm3bTfdJSUnB6NGjIZfLMWDAAGzcuLEDpRIREZE9sjiM1NTUICIiAuvWrWtX+5ycHEyfPh1TpkzBqVOnsHTpUixcuBC7du2yuFgiIiKyP5JbuVGeRCLBV199hVmzZrXZ5qmnnsL27duRnp5u2jZ79mxUVlZi586d7Xodo9GIoqIiuLu7QyKRdLRcIiIi6kZCCOh0OgQGBkIqbfv8h2NXF5Kamoq4uDizbfHx8Vi6dGmb++j1euj1etPjwsJCDBs2rKtKJCIioi5UUFCA4ODgNr/e5WGkpKQEKpXKbJtKpYJWq0VdXR2cnZ1b7JOUlIRVq1a12F5QUAAPD48uq5WIiIg6j1arhVqthru7+w3bdXkY6Yjly5cjMTHR9Ph6Zzw8PBhGiIiIepibDbHo8jDi7+8PjUZjtk2j0cDDw6PVsyIAIJfLIZfLu7o0IiIisgFdvs5IbGwskpOTzbbt3r0bsbGxXf3SRERE1ANYHEaqq6tx6tQpnDp1CkDz1N1Tp04hPz8fQPMllnnz5pna//Wvf0V2djb+/ve/IyMjA++88w62bt2Kxx57rHN6QERERD2axZdpjh8/jilTppgeXx/bMX/+fGzcuBHFxcWmYAIAYWFh2L59Ox577DG8+eabCA4Oxv/93/8hPj6+E8qn7lJUWYezRVpcKqtGQUUtKusaoatvAgBIJYC7wgl9XGXwVyoQ5uOKAX5uCOvjCqmUU7GJiOjGbmmdke6i1WqhVCpRVVXFAazdxGAUOJJ9Bd+eLsKhrCvIr6i1+Dnc5I4YGazEbf36YMJAH0QEe8KB4YSIqNdo7/s3wwiZqaxtwMepefjP0TxotD+v9eIglWCQyh0D/NwQ2scF3q4yuCucIJUATUYBXX0Tyqv1KKqsQ055DS5qqlHXaDB7bh83Ge4OD8CMiEBE9fXiWRMiIjvHMEIWqdY34Z19Wdh4OBe1Dc0hQunshLvDAzBtmApjQr3grnBq9/M1GYy4WFqNtLyrOJRVjkNZ5dBeu6wDAAFKBe4ZGYDZ0X3R39et0/tDRETWxzBC7SKEwDc/FeHF7edRqms+EzI0wAN/ndQPd40IgMyxcyZcNRqMOHzpCr45VYTvz5ZAp/85mEwa5IsF40Nx+0Bfni0hIrIjDCN0U1W1jXj6qzPYfqYYABDSxwX/uHsopg5Tdek9gOobDdh/oQxbfyzA3sxSXP8O7OfjigUTwnD/mGDIHR267PWJiKh7MIzQDZ25XIW/fHIcRVX1cJRK8Lc7BuIvk/p1ewjILa/Bx6l5+Ox4gelsib+HAv87pT/uH6OGwomhhIiop2IYoTbtTC/B0i0nUd9oRJiPK9YkRCJC7WnVmqr1Tfj8eAHW789GibYewM+hJGGsmmdKiIh6IIYRatVHh3Px3H/PQojmsRprHxhl0cDUrqZvMmDrjwVYt++SKZSovZ3x1J1DMD08oEsvHxERUediGKEWPjiYg+e/PQcAmBcbghX3DIOjQ5ffEaBD9E0GbD1+GW8nXzQNrB3V1xPPTB+KqBBvK1dHRETtwTBCZjYeysFz/20OIounDMDj0wb1iLMMtQ1N2HAgB/86cMk05Xh6eAD+MX0oAj1bv9EiERHZBoYRMtlxphj/+58TAHpWEPmlUm09Xt99AVuPF8AoABeZA5bGDcSC8WFwstGzO0REvV1737/5W9zOpeVdxWNbTgEA5seG9MggAgB+Hgqs/t1IbP/bRIwJ8UJtgwEv7cjA9Ld+wLGcCmuXR0REt4BhxI4VVdbhzx8fh77JiLihflgxY3iPDCK/NDTAA1v/EotX7hsJb1cZLmiqcf+/UpG49RTKq/U3fwIiIrI5DCN2qtFgxJJPT+JKTQOGB3rgrTmj7OYmdVKpBL8fo8bexyfhgZi+kEiAL08U4jevpuA/R/NgNNr8lUciIvoFhhE79cquTKTlXYW7whHvPhgFF5mjtUvqdJ4uMrx0bzi++t/xGBHkAW19E/7xVTrmbDiC3PIaa5dHRETtxDBih/ZlluK9A9kAgFfui0DfPi5WrqhrRao98fWiCVhxzzA4OzngaE4F7nzzAP7vh2wYeJaEiMjmMYzYmaraRiz74jQA4KFxobhzhL+VK+oeDlIJ/jghDLuW3o7xA/qgvtGI/7f9PH737mFc0OisXR4REd0Aw4idef7bc9Bo9ejn44pldw2xdjndrm8fF/z74Ris/m043OWOOFVQiXveOoi3ky+i0WC0dnlERNQKhhE7sjdDgy9OXIZEArzy+5G99iZzEokEs6P7YnfiJNwxxA8NBiNe230B975ziGdJiIhsEMOInahvNODZbWcBAA+PD+OS6QD8lQr83/wxeHN2JDxdnJBeqMU9bx/EewcucSwJEZENYRixE++mXEJhZR0ClQokThtk7XJshkQiwczIIHy/9Hb8ZogfGpqMeGlHBua8dwT5V2qtXR4REYFhxC7kX6nFu/svAQCeuWeYXU7jvVV+Hgq8P38MXv5dOFxlDjiW2zzjZtPRfPSAOyIQEdm1DoWRdevWITQ0FAqFAjExMTh27NgN269ZswaDBw+Gs7Mz1Go1HnvsMdTX13eoYGrphe3n0NBkxPgBfXBXL5k90xESiQQJY/ti59LbER3mjdoGA57+6gwWbPwRGi2/H4mIrMXiMLJlyxYkJiZi5cqVOHHiBCIiIhAfH4/S0tJW22/atAnLli3DypUrcf78ebz//vvYsmULnn766VsunoAfcyuw+5wGDlIJnrOD5d67g9rbBZv/dBuemT4UMkcpUjLLMO2NA/jmpyJrl0ZE1CtZHEZef/11/OlPf8KCBQswbNgwrF+/Hi4uLvjggw9abX/48GGMHz8eDzzwAEJDQzFt2jTMmTPnpmdT6OaEEHj5uwwAwP1j1BiocrdyRT2HVCrBwon9sH3JBIQHKVFV14i/fXoSizadQGVtg7XLIyLqVSwKIw0NDUhLS0NcXNzPTyCVIi4uDqmpqa3uM27cOKSlpZnCR3Z2Nnbs2IG77767zdfR6/XQarVmH9TSvsxSHM+7CrmjFI/eMdDa5fRIA1Xu+PJ/x2Fp3EA4SCXYfroY8WsOYP+FMmuXRkTUa1gURsrLy2EwGKBSqcy2q1QqlJSUtLrPAw88gOeffx4TJkyAk5MT+vfvj8mTJ9/wMk1SUhKUSqXpQ61WW1Jmr2A0CvxzZyYA4KHxofBXKqxcUc/l5CDF0rhB+PKRcejn6wqNVo/5HxzDiq/TUddgsHZ5RER2r8tn06SkpOCll17CO++8gxMnTuDLL7/E9u3b8cILL7S5z/Lly1FVVWX6KCgo6Ooye5ydZ0uQUaKDu8IRj0zqb+1y7EKE2hPbl0zE/NgQAMDHqXmY/tYP+Kmg0rqFERHZOYvmgPr4+MDBwQEajcZsu0ajgb9/67M4nn32WcydOxcLFy4EAISHh6OmpgZ//vOf8Y9//ANSacs8JJfLIZfLLSmtVxFCYN2+LADAgvFh8HSRWbki++Esc8CqmSNwx1AVnvz8J2SX1+C37x7Gkt8MwKIpA+DkwNnwRESdzaLfrDKZDFFRUUhOTjZtMxqNSE5ORmxsbKv71NbWtggcDg7Ny5RzfYeOSblQhrNFWrjIHLBgXKi1y7FLtw/yxa6lt+OekQEwGAXW7LmI+9anIrus2tqlERHZHYv/zEtMTMSGDRvw0Ucf4fz583jkkUdQU1ODBQsWAADmzZuH5cuXm9rPmDED7777LjZv3oycnBzs3r0bzz77LGbMmGEKJdR+Qgis29t8VuTBmL7wcuVZka7i6SLD2gdG483ZkfBQOOKngkrc/dYP+Dg1l0GaiKgTWbxUZ0JCAsrKyrBixQqUlJQgMjISO3fuNA1qzc/PNzsT8swzz0AikeCZZ55BYWEhfH19MWPGDLz44oud14te5FhOBY7nXYXMUYo/Texn7XJ6hZmRQYgO88aTn53GwaxyrPj6LPacL8Ur942EyoMDh4mIbpVE9IA/8bRaLZRKJaqqquDh4WHtcqxq4UfHsee8Bg/G9MWL94Zbu5xexWgU+Dg1F0nfZUDfZITS2Qkv3jsC94wMtHZpREQ2qb3v3xyN14PkXalBckbz4OE/TgizcjW9j1QqwUPjw7D9bz8vlLZ400ks3XwSVXWN1i6PiKjHYhjpQT46nAchgEmDfNHf183a5fRaA/yaF0r7228GQCoBtp0qwp1rDuBQVrm1SyMi6pEYRnoIXX0jth5vXm+FZ0Wsz8lBisRpg/H5I+MQ2scFxVX1ePD/jmLVf8+ivpELpRERWYJhpIf4PO0yqvVN6O/ritsH+li7HLpmdF8v7Hh0Ih6M6QsA+PBQLu55+yDSC6usXBkRUc/BMNIDCCHwyZE8AMBD48N4Z14b4yJzxIv3huPDh8bC112OrNJqzFp3CGv3XkSTwWjt8oiIbB7DSA9wLKcC2WU1cJE54N5RQdYuh9owZYgfdi29HXeN8EeTUeDV7y/g/n+lIre8xtqlERHZNIaRHuDTY/kAgJmRgXCTW7w0DHUjb1cZ3nlwNF6/PwLuckecyG9eKG3T0XwulEZE1AaGERt3taYBO9Kb74g8e2xfK1dD7SGRSPDb0cH4bulE3NbPG7UNBjz91Rk8/NFxlOrqrV0eEZHNYRixcV+eLERDkxHDAjwwMlhp7XLIAsFeLti08DY8M30oZA5S7M0oRfwbB7AzvdjapRER2RSGERsmhDBdopkT05cDV3sgqVSChRP74b9LJmBogAeu1jbir/8+gce3/gRtPRdKIyICGEZs2unLVcgqrYbcUYqZkVxyvCcb7O+ObYvG4ZHJ/SGRAF+cuIy71vyAI9lXrF0aEZHVMYzYsC9PXAYAxA/3h4fCycrV0K2SOzrgqTuHYOtfYqH2dkZhZR3mbDiCF7ef40JpRNSrMYzYqIYmI/57unlswW9HczqvPRkb6o3vHr0ds8eqIQSw4Ycc3P3WD0jLq7B2aUREVsEwYqP2XyhDRU0DfNzkmDCAK67aGze5I1b/biT+b94Y+LnLkV1Wg/vWp+KFb8+hroFnSYiod2EYsVFfnWy+RDMrMhCODjxM9ipumAq7H5uE340OhhDA+wdzcNebB3Ash2dJiKj34LucDaqqbcSec6UAgN+ODrZyNdTVlC5OeO3+CHz40Fj4eyiQe6UWCe+l4rlvzqK2ocna5RERdTmGERu082wxGgxGDFa5Y1igh7XLoW4yZYgfvk+8HQljmseSbDycizvX/IDUS5xxQ0T2jWHEBn17beDq/3A6b6/joXDCy/eNxMd/jEaQpzPyK2oxZ8MRPLPtDHRcl4SI7BTDiI25Uq3H4Wt/CU8PD7ByNWQttw/yxc6lE/FgTPMtAP59JB9TXz+A78+WWLkyIqLOxzBiY3aeLYHBKDAiyAOhPq7WLoesyF3hhBfvDcemhTEI6eOCEm09/vxJGh75dxo0Wt7jhojsB8OIjdl+7RLNPSN5iYaajRvgg11Lb8f/Tu4PR6kE36WXIO71/fjP0TwYjbwTMBH1fB0KI+vWrUNoaCgUCgViYmJw7NixG7avrKzEokWLEBAQALlcjkGDBmHHjh0dKtielen0puXBeYmGfknh5IC/3zkE/10yARFqT+jqm/CPr9KR8F4qskp11i6PiOiWWBxGtmzZgsTERKxcuRInTpxAREQE4uPjUVpa2mr7hoYGTJ06Fbm5ufj888+RmZmJDRs2ICiIq4r+2s70YhgFEKH2hNrbxdrlkA0aGuCBLx8Zh5UzhsFF5oAfc6/irjd/wBu7L0DfxMXSiKhnkgghLDrPGxMTg7Fjx2Lt2rUAAKPRCLVajSVLlmDZsmUt2q9fvx6vvPIKMjIy4OTUvvur6PV66PV602OtVgu1Wo2qqip4eNjvVNcHNhzB4UtX8PTdQ/Dn2/tbuxyycYWVdVixLR3JGc1/CPT3dcULs0ZgXH+u2EtEtkGr1UKpVN70/duiMyMNDQ1IS0tDXFzcz08glSIuLg6pqamt7vPNN98gNjYWixYtgkqlwogRI/DSSy/BYGj7r7ikpCQolUrTh1qttqTMHqmytgFHr626GT/c38rVUE8Q5OmM/5s/BmsfGAUfNzkuldXggQ1H8ejmkyjlAFci6kEsCiPl5eUwGAxQqVRm21UqFUpKWp9ymJ2djc8//xwGgwE7duzAs88+i9deew3/7//9vzZfZ/ny5aiqqjJ9FBQUWFJmj5R8vhQGo8AQf3eE9OEsGmofiUSCe0YGIjlxEubeFgKJBPj6VBHueG0/PjiYgyaD0dolEhHdVJfPpjEajfDz88N7772HqKgoJCQk4B//+AfWr1/f5j5yuRweHh5mH/Zu17X1I6bxrAh1gNLFCS/MGoFvFl0b4KpvwvPfnsOMtYd4N2AisnkWhREfHx84ODhAo9GYbddoNPD3b/1NNCAgAIMGDYKDg4Np29ChQ1FSUoKGhoYOlGx/6hoMOHCxDAAQP1x1k9ZEbQsPVuKrR8bhpXvDoXR2wvliLX73biqe/OwnXKnW3/wJiIiswKIwIpPJEBUVheTkZNM2o9GI5ORkxMbGtrrP+PHjkZWVBaPx59PFFy5cQEBAAGQyWQfLti/7L5ShvtGIYC9nDAuw/7NA1LWkUgkeiOmLvY9Pwv1jmm+0+FnaZfzmtf3495E8GLg2CRHZGIsv0yQmJmLDhg346KOPcP78eTzyyCOoqanBggULAADz5s3D8uXLTe0feeQRVFRU4NFHH8WFCxewfft2vPTSS1i0aFHn9aKHu77Ed/xwf0gkEitXQ/aij5sc/7wvAl88Eosh/u6oqmvEM9vSMePtgziazZvvEZHtcLR0h4SEBJSVlWHFihUoKSlBZGQkdu7caRrUmp+fD6n054yjVquxa9cuPPbYYxg5ciSCgoLw6KOP4qmnnuq8XvRgBqPA3szmqZnThvESDXW+qBBvfLtkAj5OzcMbey7gXLEWCe8dwfTwACy7awjXtCEiq7N4nRFraO885Z7oeG4F7lufCqWzE9KeiYOjA1fop65zpVqP13ZfwOZj+TAKQOYoxV9u74dHJveHi8ziv02IiG6oS9YZoc53fcGqyYN9GUSoy/Vxk+Ole8Px7ZKJuK2fNxqajHh7bxZ+8+p+fHXyMu91Q0RWwXc/K9t7vjmM/GaIn5Urod5kWKAHPv3TbVj/h9FQezujRFuPx7b8hN+tP4xTBZXWLo+IehmGESsqqKhFpkYHqQSYNMjX2uVQLyORSHDniADsfmwSnowfDBeZA07mV2LWukP426cnUVBRa+0SiaiXYBixon3XBq6OCfGGpwunOZN1KJwcsGjKAKQ8MRm/Gx0MiQT45qci/Oa1FDz/33O4WsP1gIioazGMWFHy9Us0Q3mJhqzPz0OB1+6PwH8XT8DEgT5oNAh8cCgHt7+yD++mXEJ9I+8KTERdg2HESmobmpB6ba2HOzhehGzIiCAlPnk4Bh//MRpDAzygq2/CyzszMOXVFHx2vICLphFRp2MYsZKjORVoaDIiyNMZA/zcrF0OUQu3D/LF9iUT8Pr9EQhUKlBcVY8nPz+N6W/9gJTMUvSAVQGIqIdgGLGSAxea70Vz+yAfrrpKNksqleC3o4Ox94nJePruIfBQOCKjRIeHPvwRCe8dwbEc3oSPiG4dw4iV/HCxHAAwcSBn0ZDtUzg54M+398eBv0/BnyaGQeYoxbGcCtz/r1TMff8opwMT0S1hGLGCoso6ZJVWQyoBxvf3sXY5RO3m6SLDP6YPw/4nJ+PBmL5wlErww8VyzFp3CAs/+hFni6qsXSIR9UAMI1bww8XmSzSRak8oXZysXA2R5QKUznjx3nDse2Iy7osKhlQC7DlfiulvHcSi/5xAVqnO2iUSUQ/CMGIFBy7wEg3ZB7W3C179fQR2J07C/0QEQiIBtp8pxrQ3DiBxyylkl1Vbu0Qi6gEYRrqZwShwMKs5jNzOVVfJTvT3dcNbc0bhu0cnIn64CkYBfHmyEHGv78ffPj2JzBKeKSGitjGMdLPTlytRVdcId4UjIoKV1i6HqFMN8ffAv+aOwTeLx+OOIX4wiubVXOPXHMBfPjmO9EKOKSGilhhGutn1WTQTBvjwLr1kt0YGe+L9h8bi2yUTcHe4PyQSYNdZDe55+yAWfHgMaXlXrV0iEdkQR2sX0NtcX1+E40WoNxgRpMQ7D0bhokaHdfuy8M1PRdiXWYZ9mWUY178PlvxmIG7r5821doh6Of5p3o209Y04eW09hokDOaWXeo+BKnesmT0Kex+fjIQxajhKJTh86QrmbDiCe985jO/OFHOZeaJejGGkGx3OugKDUaCfjyvU3i7WLoeo24X6uOLl+0Yi5cnJmHtbCGSOUpwqqMQj/zmBO15LwSdH8nhDPqJeiGGkG11fX4SzaKi3C/ZywQuzRuDQU7/B4ikDoHR2Qu6VWjy7LR3jVu/Fmj0XUFHTYO0yiaibMIx0EyEEDly8Pl6El2iIAMDXXY4n4gfj8LLf4LkZwxDs5YyKmgas2XMR41Yn45ltZ5BbXmPtMomoi0lED7j1plarhVKpRFVVFTw8PKxdTofkltdg8qspcHKQ4NSKaXCVc+ww0a81GYz4Lr0E7x3Ixplr04AlEmDqUBUeGh+K2H59ONiVqAdp7/t3h86MrFu3DqGhoVAoFIiJicGxY8fatd/mzZshkUgwa9asjrxsj3b9rEhUiBeDCFEbHB2kmBERiG8Wj8emP8Vg8mBfCAF8f06DBzYcxZ1rfsCnx/JR18BxJUT2xOIwsmXLFiQmJmLlypU4ceIEIiIiEB8fj9LS0hvul5ubiyeeeAITJ07scLE92UHepZeo3SQSCcb198HGBdHY/djt+MNtfeHs5IBMjQ7LvzyD25KSkbTjPC5frbV2qUTUCSy+TBMTE4OxY8di7dq1AACj0Qi1Wo0lS5Zg2bJlre5jMBhw++23449//CN++OEHVFZWYtu2bW2+hl6vh16vNz3WarVQq9U99jKNwSgw6vnvoa1vwrZF4xGp9rR2SUQ9TlVdIz47XoCPUnNRUFEHAJBKgKnDVHhoXBjXKyGyQV1ymaahoQFpaWmIi4v7+QmkUsTFxSE1NbXN/Z5//nn4+fnh4YcfbtfrJCUlQalUmj7UarUlZdqcs0VV0NY3wV3hiBGBPS9MEdkCpbMTFk7sh5QnpmDDvDGYMMAHRtG8suucDUcw7Y0D+PBQDqpqG61dKhFZyKIwUl5eDoPBAJVKZbZdpVKhpKSk1X0OHjyI999/Hxs2bGj36yxfvhxVVVWmj4KCAkvKtDmHsq4AAGLC+nAJeKJb5CCVYOowFf69MAbfP3Y7HoxpvoRzsbQaq/57DtEv7cHjW39CWt5V9IDx+USELl4OXqfTYe7cudiwYQN8fNo/nVUul0Mul3dhZd3r8KXm8SLjB/SxciVE9mWQyh0v3huOp+4agm0nC7HpaD4ySnT44sRlfHHiMob4u+OBmL6YNSoIHgona5dLRG2wKIz4+PjAwcEBGo3GbLtGo4G/v3+L9pcuXUJubi5mzJhh2mY0Gptf2NERmZmZ6N+/f0fq7jH0TQb8mFsBABjXn+uLEHUFD4UT5sWGYu5tITiRX4lNR/Px7ekiZJTosOLrs0jakYEZEQGYHd0Xo9SeHFtCZGMsCiMymQxRUVFITk42Tc81Go1ITk7G4sWLW7QfMmQIzpw5Y7btmWeegU6nw5tvvtnjx4K0x6n8StQ3GuHjJscglZu1yyGyaxKJBFEhXogK8cKKe4bhy5OXseloPi6WVmPr8cvYevwy+vu64r4oNX47OggqD4W1SyYidOAyTWJiIubPn48xY8YgOjoaa9asQU1NDRYsWAAAmDdvHoKCgpCUlASFQoERI0aY7e/p6QkALbbbq0OXmseLjOvPxZqIupPSxQkLxofhoXGhOJ53FZ8ezceO9GJcKqvByzsz8MquDEwa5Iv7otSIG+YHuaODtUsm6rUsDiMJCQkoKyvDihUrUFJSgsjISOzcudM0qDU/Px9SKQdpXpd6bbzIuP4cL0JkDRKJBGNDvTE21BurZg7H9tPF+CztMtLyrmJfZhn2ZZbB08UJMyMCcV+UGiOCPPiHA1E343LwXahG34SIVd+jySjww9+n8E69RDYku6wan6ddxpcnClGirTdtH6Ryw8zIIMyMDESwF39miW5Fe9+/GUa6UEpmKR768EcEeznj4FO/sXY5RNQKg1HgYFY5PjtegO/PadDQZDR9bUyIF2aOCsL08AB4u8qsWCVRz9Te92/eJKULHb42XmQ8Z9EQ2SwHqQSTBvli0iBfVNU1Ymd6MbadLMKRnCs4nncVx/OuYtU3Z3H7IF/MjAzE1GEquMj4q5OoM/EnqgtdX19kHNcXIeoRlM5OSBjbFwlj+6K4qg7f/lSMbacKcbZIi70ZpdibUQoXmQOmDlPh7vAATBrkC4UTB74S3SpepukilbUNGPXCbggBHPvHHfBz5xRCop4qq1SHr08V4etTRciv+PnmfK4yB9wxVIW7w/0xebAfgwnRr/AyjZUdyb4CIYCBfm4MIkQ93AA/dzw+bTASpw7CyYJKbD9djO/OFKOoqh7f/FSEb34qgovMAb8Z4ofp4QGYPNgPzjIGE6L2YhjpItfvRzN+AMeLENkLiUSC0X29MLqvF/5x91CculyJHaeL8V16CQor6/Dt6WJ8e7oYLjIHTBnsh6nDVJgy2A9KFy5FT3QjDCNd5Pp4kViuL0Jkl6TSXwST6UNxqqASO84UY8eZ5mCy/Uwxtp8phqNUgugwb0wdpsLUYSpOFyZqBceMdAGNth4xLyVDIgFOPTuNfxUR9SJCCJy+XIXvz5Vg9zkNLmiqzb4+NMADU4epMG2YCsMDucAa2TeOGbGiQ1nNZ0XCg5QMIkS9jEQiQYTaExFqTzwZPwR5V2qw+5wG35/T4HhuBc4Xa3G+WIu3ki8iUKlA3DAV4oaqENPPm0vSU6/FMNIFro8X4V16iSikjysWTuyHhRP7oaKmAXszSrH7XAkOXChHUVU9Pk7Nw8epeXB2csC4/n0webAvJg/244rN1KswjHQyIYRpvMh4ri9CRL/g7SrDfVHBuC8qGPWNBhzKKsfucxokZ5SiTKdHckYpkjNKAZxFPx9XTLoWTGLCvDltmOwaw0gnyy6vQXFVPWSOUowN9bZ2OURkoxROzWuU3DFUBSEEzhVrkZJZhv0XypCWdxXZ5TXILq/Bh4dyoXCSIrZfH0we7IfJg30R0sfV2uUTdSqGkU52+Np4kai+XvxLhojaRSKRYHigEsMDlVg0ZQC09Y04dLHcFE5KtPWmOwwDgNrbGeP7+2D8AB+M698HfdzkVu4B0a1hGOlkP68vwks0RNQxHgon3BUegLvCAyCEQKZGh5TMMqRkluJ47lUUVNRhc0UBNv9YAAAY4u+O8QN8MH5AH0SH9YGbnL/aqWfhd2wnMhgFUrOvDV7lYmdE1AkkEgmG+HtgiL8H/jqpP6r1TfgxpwKHsspx6NIVnC/WIqNEh4wSHd4/mANHafNsnvH9+2DcAB+M6uvJWTpk8xhGOtG5Ii2q6hrhLnfEyCCltcshIjvkJnfElCF+mDLEDwBQXq1H6qUrOHypHIeyriC/ohZpeVeRlncVb+3NgtxRilF9PREd1gcxYd4Y1deTdx0mm8PvyE506Nosmph+3nB0kFq5GiLqDXzc5JgREYgZEYEAgIKKWtNZk9RL5SivbsCR7Aocya4AADhKJQgPViI6zBsxYd6ICvGG0pnrIZF1MYx0ouuLnXF9ESKyFrW3C2ZH98Xs6L4QQuBSWTWO5lTgWE4FjmZXoERbj5P5lTiZX4l/7c+GRAIM9fdAdJg3osO8ERXiBZUHb+5J3YthpJPomwz4Mbf5Lw/eHI+IbIFEIsEAP3cM8HPHgzEhEELg8tW6a+HkCo7lVCD3Si3OFWtxrliLjYdzAQBBns4Y1dcTUSHN994ZFugBJ57tpS7EMNJJTuZXor7RCB83OQap3KxdDhFRCxKJBGpvF6i9XXBfVDCA5ntpHcupwI+5zWdPLmh0KKysM92FGADkjlKMDFY23xjwWkDxded0Yuo8DCOd5Pr6IuMH9OGNr4iox1B5KMzGnFTrm/BTQSVO5F3FifyrOFlQicraRvyYexU/5l417af2dkZEsCdGBisxMtgTI4KUnFJMHdah75x169bhlVdeQUlJCSIiIvD2228jOjq61bYbNmzAxx9/jPT0dABAVFQUXnrppTbb91QHr4cRjhchoh7MTe54bc2S5t9lQghkl9cgLe8qTuZfxYm8Slwo1aGgog4FFT+fPZFIgH4+rogI9kR4sBIjg5UYFqCEs4zTiunmLA4jW7ZsQWJiItavX4+YmBisWbMG8fHxyMzMhJ+fX4v2KSkpmDNnDsaNGweFQoGXX34Z06ZNw9mzZxEUFNQpnbA2XX0jfrpcBQAYx8XOiMiOSCQS9Pd1Q39fN9w/Rg0A0NY34nRBFU4XVuJ0QRXOFFahsLIOl8pqcKmsBl+eLAQAOEglGOjnZjp7MjJYicH+7lz3hFqQCCGEJTvExMRg7NixWLt2LQDAaDRCrVZjyZIlWLZs2U33NxgM8PLywtq1azFv3rxW2+j1euj1etNjrVYLtVqNqqoqeHh4WFJut0g+r8HDHx1HSB8X7H9yirXLISLqdmU6PdILq3D6chVOX67ET5erUF6tb9HOyaE53AwL9MCwgOaPoQEe8HKVWaFq6mparRZKpfKm798WnRlpaGhAWloali9fbtomlUoRFxeH1NTUdj1HbW0tGhsb4e3d9k3kkpKSsGrVKktKs6rrS8BzSi8R9Va+7nKzxdiEECjR1uP05SqcuVyFny5X4kxhFSprG00rxn6JQtP+gUoFhgZ4mELK0AAP9PV2gVTKMXi9gUVhpLy8HAaDASqVymy7SqVCRkZGu57jqaeeQmBgIOLi4tpss3z5ciQmJpoeXz8zYqsOXGy+edUETuklIgLQfHknQOmMAKUz4of7A4BpavH5Yi3OF+twrrgK54t1yK+oRVFVPYqq6pGcUWp6DleZA4ZcO3sy2N8dg/3dMcjPHUoXLtJmb7p16PPq1auxefNmpKSkQKFoe1EduVwOubxnTBu7fLUWWaXVcJBKMGEgwwgRUVt+ObV42rWAAjSPQcko1uF8sRbnirQ4X9J8v52aBoNpaftf8nOXY7C/Owb6uWOQyg0DVc2f3RUMKT2VRWHEx8cHDg4O0Gg0Zts1Gg38/f3b2KvZq6++itWrV2PPnj0YOXKk5ZXaqAMXmmfRjFJ7ckllIqIO8FA4mVaAva7JYER2eY0poGRqdLioqUZhZR1KdXqU6vT44WK52fMEKhUYqHK/FlTcMEjljn6+rgwpPYBFYUQmkyEqKgrJycmYNWsWgOYBrMnJyVi8eHGb+/3zn//Eiy++iF27dmHMmDG3VLCtSclsPqU4aZCvlSshIrIfjg5SDFK5Y5DKHTMjf555qatvxMXSalzU6HBBU40LGh0uaHTQaPWmSz37L5SZPZefuxz9fd3Qz9fV7HOgpzMcOCbFJlh8mSYxMRHz58/HmDFjEB0djTVr1qCmpgYLFiwAAMybNw9BQUFISkoCALz88stYsWIFNm3ahNDQUJSUlAAA3Nzc4ObWs1cqbWgy4vCl5sGrkwYzjBARdTV3hVPzSrB9vcy2V9U24mKpeUC5oKlGebXedCYlNfuK2T5yRynCfH4OKNdDSpgPz6Z0N4vDSEJCAsrKyrBixQqUlJQgMjISO3fuNA1qzc/Ph1T68z0M3n33XTQ0NOC+++4ze56VK1fiueeeu7XqrexE/lVU65vQx1WGEYFKa5dDRNRrKV2cMCbUG2NCzWdqausbkV1Wg0ul1cgur27+d1k1cstroW8ymmb2/FofVxn69nFBaB9XhFz7fP2xl4sTV9ruZBavM2IN7Z2n3N1Wf5eB9fsv4d5RQXgjIdLa5RARUTsZjAKXr9aawsmla5+zy2paXR/ll9wVjgjp44KQPq4IvfY5xNsFoT6u8HOXM6j8QpesM0LmOF6EiKhncpBKmkNEH1fT2ijX6eobkXelFvkVtci9UoO88ubP+RW1KK6qh66+CemFWqQXals8r7OTA9Tezgj2ckGwl/O1DxfTZ55VaR3DSAcVVNQio0QHqQSYyCm9RER2w13hhBFBSowIann5vb7RgPyKWuRdqUXelZrmsHKl+fHlq7WoazRcG7dS3epzu8oc2gwqwV7O8OylYYVhpIN2n2ue3jwm1Bt93HrGmihERHRrFE4Oplk+v9bQZERhZR0uX61FQUXz58tXf/5cqtOjpsGATI0OmZqW41SA5rAS6OmMAE9nBHgoEOCpQKDSGf5KBQI9FQhQOsPVDu+ObH896ibfn2ueFTRtmOomLYmIqDeQXZudE+bj2urX6xsNKKqsuxZQ2g4rF0urcbG09TMrQPOYlUClMwI8FQhQKq6tdHvt87VtLrKe9fbes6q1EVdrGnAspwIAMG3YjRd7IyIiAprPqvTzdUM/39aXtahvNKCwsg7FlfUorqpDcdUvPlfWo6iqDrr6Jujqm5BZ3/bZFQBQOjshQKmAn4cCKnc5VB4KqDzkzY+v/dvHTQ4nB2mbz9GdGEY6IDmjFEYBDPF3R98+LtYuh4iI7IDCyQH9fd3Qv42wAgDV+iaUVNWhqLIeJVXNAaW4sh7F2noUVzYHl2p9E6rqGlFV19jqtOXrJBKgj6scKo/msPLEtMEYFmidGasMIx3w/dnmSzTxw3lWhIiIuo+b3BED/NwxwK/lmJXrdPWNKK6qh0ZbD41WD422HqXX/62rh6aqHqU6PZqMAuXVepRX63G2SIulcQO7sSfmGEYsVNdgMN2ld9pwjhchIiLb4q5wgrvCqdVBttcZjQIVtQ3XgkpzYAnp0/pYl+7AMGKh3ec1qG80ItjLGcMCbGcBNiIiovaSSiXwcWseNzI80NrVALYxcqUH+SLtMgDg3lFBvXIuOBERUWdjGLGARluPH65dovnt6GArV0NERGQfGEYs8PWpQhgFEBXi1eY8ciIiIrIMw0g7CSHwRVohAOB3PCtCRETUaRhG2ulskRaZGh1kjlJMHxlg7XKIiIjsBsNIO314KBcAMHWYCkpnJ+sWQ0REZEcYRtqhoKIW2041X6L508R+Vq6GiIjIvjCMtMP6/ZdgMApMHOiDSLWntcshIiKyKwwjN6HR1uOz481riyyaMsDK1RAREdkfhpGbWLs3Cw0GI8aGeiEmzNva5RAREdkdhpEb2JdRik+O5AEAHr1jEFdcJSIi6gIdCiPr1q1DaGgoFAoFYmJicOzYsRu2/+yzzzBkyBAoFAqEh4djx44dHSq2O5VU1ePxz34CADw0LhQTBvpYuSIiIiL7ZHEY2bJlCxITE7Fy5UqcOHECERERiI+PR2lpaavtDx8+jDlz5uDhhx/GyZMnMWvWLMyaNQvp6em3XHxXKaioxV/+nYaKmgYMC/DAsruGWLskIiIiuyURQghLdoiJicHYsWOxdu1aAIDRaIRarcaSJUuwbNmyFu0TEhJQU1ODb7/91rTttttuQ2RkJNavX9+u19RqtVAqlaiqqoKHR+feKVcIAaMA9E0GXCqtQWp2Od7ccxE1DQa4yR3xzeLx6Ofr1qmvSURE1Bu09/3b0ZInbWhoQFpaGpYvX27aJpVKERcXh9TU1Fb3SU1NRWJiotm2+Ph4bNu2rc3X0ev10Ov1psdardaSMtvt3ncO4WR+ZatfGxvqhZd/N5JBhIiIqItZFEbKy8thMBigUqnMtqtUKmRkZLS6T0lJSavtS0pK2nydpKQkrFq1ypLSOuTXw1G9XWUY6OeGeyIC8WB0X0ilHLBKRETU1SwKI91l+fLlZmdTtFot1Gp1p7/Ohw9FwyAEpBLAQSqBu4LLvBMREXU3i8KIj48PHBwcoNFozLZrNBr4+/u3uo+/v79F7QFALpdDLpdbUlqHKF0YPoiIiKzNotk0MpkMUVFRSE5ONm0zGo1ITk5GbGxsq/vExsaatQeA3bt3t9meiIiIeheLL9MkJiZi/vz5GDNmDKKjo7FmzRrU1NRgwYIFAIB58+YhKCgISUlJAIBHH30UkyZNwmuvvYbp06dj8+bNOH78ON57773O7QkRERH1SBaHkYSEBJSVlWHFihUoKSlBZGQkdu7caRqkmp+fD6n05xMu48aNw6ZNm/DMM8/g6aefxsCBA7Ft2zaMGDGi83pBREREPZbF64xYQ1euM0JERERdo0vWGbGW63mpq9YbISIios53/X37Zuc9ekQY0el0ANAl03uJiIioa+l0OiiVyja/3iMu0xiNRhQVFcHd3b1T75x7ff2SgoICu738wz72fPbeP4B9tAf23j/A/vvYFf0TQkCn0yEwMNBsPOmv9YgzI1KpFMHBwV32/B4eHnb5jfVL7GPPZ+/9A9hHe2Dv/QPsv4+d3b8bnRG5zuK79hIRERF1JoYRIiIisqpeHUbkcjlWrlzZLUvPWwv72PPZe/8A9tEe2Hv/APvvozX71yMGsBIREZH96tVnRoiIiMj6GEaIiIjIqhhGiIiIyKoYRoiIiMiq7C6MrFu3DqGhoVAoFIiJicGxY8du2P6zzz7DkCFDoFAoEB4ejh07dph9XQiBFStWICAgAM7OzoiLi8PFixe7sgs3ZEn/NmzYgIkTJ8LLywteXl6Ii4tr0f6hhx6CRCIx+7jzzju7uhs3ZEkfN27c2KJ+hUJh1sbWjiFgWR8nT57coo8SiQTTp083tbGl43jgwAHMmDEDgYGBkEgk2LZt2033SUlJwejRoyGXyzFgwABs3LixRRtLf7a7kqV9/PLLLzF16lT4+vrCw8MDsbGx2LVrl1mb5557rsUxHDJkSBf2om2W9i8lJaXV79GSkhKzdj35GLb2MyaRSDB8+HBTG1s6hklJSRg7dizc3d3h5+eHWbNmITMz86b7Wes90a7CyJYtW5CYmIiVK1fixIkTiIiIQHx8PEpLS1ttf/jwYcyZMwcPP/wwTp48iVmzZmHWrFlIT083tfnnP/+Jt956C+vXr8fRo0fh6uqK+Ph41NfXd1e3TCztX0pKCubMmYN9+/YhNTUVarUa06ZNQ2FhoVm7O++8E8XFxaaPTz/9tDu60ypL+wg0rxb4y/rz8vLMvm5LxxCwvI9ffvmlWf/S09Ph4OCA3//+92btbOU41tTUICIiAuvWrWtX+5ycHEyfPh1TpkzBqVOnsHTpUixcuNDszboj3xddydI+HjhwAFOnTsWOHTuQlpaGKVOmYMaMGTh58qRZu+HDh5sdw4MHD3ZF+Tdlaf+uy8zMNKvfz8/P9LWefgzffPNNs74VFBTA29u7xc+hrRzD/fv3Y9GiRThy5Ah2796NxsZGTJs2DTU1NW3uY9X3RGFHoqOjxaJFi0yPDQaDCAwMFElJSa22v//++8X06dPNtsXExIi//OUvQgghjEaj8Pf3F6+88orp65WVlUIul4tPP/20C3pwY5b279eampqEu7u7+Oijj0zb5s+fL2bOnNnZpXaYpX388MMPhVKpbPP5bO0YCnHrx/GNN94Q7u7uorq62rTN1o7jdQDEV199dcM2f//738Xw4cPNtiUkJIj4+HjT41v9f9aV2tPH1gwbNkysWrXK9HjlypUiIiKi8wrrJO3p3759+wQAcfXq1Tbb2Nsx/Oqrr4REIhG5ubmmbbZ6DIUQorS0VAAQ+/fvb7ONNd8T7ebMSENDA9LS0hAXF2faJpVKERcXh9TU1Fb3SU1NNWsPAPHx8ab2OTk5KCkpMWujVCoRExPT5nN2lY7079dqa2vR2NgIb29vs+0pKSnw8/PD4MGD8cgjj+DKlSudWnt7dbSP1dXVCAkJgVqtxsyZM3H27FnT12zpGAKdcxzff/99zJ49G66urmbbbeU4WupmP4ed8f/M1hiNRuh0uhY/ixcvXkRgYCD69euHBx98EPn5+VaqsGMiIyMREBCAqVOn4tChQ6bt9ngM33//fcTFxSEkJMRsu60ew6qqKgBo8T33S9Z8T7SbMFJeXg6DwQCVSmW2XaVStbhueV1JSckN21//bMlzdpWO9O/XnnrqKQQGBpp9I9155534+OOPkZycjJdffhn79+/HXXfdBYPB0Kn1t0dH+jh48GB88MEH+Prrr/Hvf/8bRqMR48aNw+XLlwHY1jEEbv04Hjt2DOnp6Vi4cKHZdls6jpZq6+dQq9Wirq6uU773bc2rr76K6upq3H///aZtMTEx2LhxI3bu3Il3330XOTk5mDhxInQ6nRUrbZ+AgACsX78eX3zxBb744guo1WpMnjwZJ06cANA5v79sSVFREb777rsWP4e2egyNRiOWLl2K8ePHY8SIEW22s+Z7Yo+4ay/dutWrV2Pz5s1ISUkxG+A5e/Zs07/Dw8MxcuRI9O/fHykpKbjjjjusUapFYmNjERsba3o8btw4DB06FP/617/wwgsvWLGyrvH+++8jPDwc0dHRZtt7+nHsTTZt2oRVq1bh66+/NhtTcdddd5n+PXLkSMTExCAkJARbt27Fww8/bI1S223w4MEYPHiw6fG4ceNw6dIlvPHGG/jkk0+sWFnX+Oijj+Dp6YlZs2aZbbfVY7ho0SKkp6dbbfxKe9jNmREfHx84ODhAo9GYbddoNPD39291H39//xu2v/7ZkufsKh3p33WvvvoqVq9eje+//x4jR468Ydt+/frBx8cHWVlZt1yzpW6lj9c5OTlh1KhRpvpt6RgCt9bHmpoabN68uV2/1Kx5HC3V1s+hh4cHnJ2dO+X7wlZs3rwZCxcuxNatW1ucDv81T09PDBo0qEccw9ZER0ebarenYyiEwAcffIC5c+dCJpPdsK0tHMPFixfj22+/xb59+xAcHHzDttZ8T7SbMCKTyRAVFYXk5GTTNqPRiOTkZLO/nH8pNjbWrD0A7N6929Q+LCwM/v7+Zm20Wi2OHj3a5nN2lY70D2ge+fzCCy9g586dGDNmzE1f5/Lly7hy5QoCAgI6pW5LdLSPv2QwGHDmzBlT/bZ0DIFb6+Nnn30GvV6PP/zhDzd9HWseR0vd7OewM74vbMGnn36KBQsW4NNPPzWblt2W6upqXLp0qUccw9acOnXKVLu9HEOgeZZKVlZWu/4osOYxFEJg8eLF+Oqrr7B3716EhYXddB+rvife0vBXG7N582Yhl8vFxo0bxblz58Sf//xn4enpKUpKSoQQQsydO1csW7bM1P7QoUPC0dFRvPrqq+L8+fNi5cqVwsnJSZw5c8bUZvXq1cLT01N8/fXX4vTp02LmzJkiLCxM1NXV2Xz/Vq9eLWQymfj8889FcXGx6UOn0wkhhNDpdOKJJ54QqampIicnR+zZs0eMHj1aDBw4UNTX13d7/zrSx1WrVoldu3aJS5cuibS0NDF79myhUCjE2bNnTW1s6RgKYXkfr5swYYJISEhosd3WjqNOpxMnT54UJ0+eFADE66+/Lk6ePCny8vKEEEIsW7ZMzJ0719Q+OztbuLi4iCeffFKcP39erFu3Tjg4OIidO3ea2tzs/1l3s7SP//nPf4Sjo6NYt26d2c9iZWWlqc3jjz8uUlJSRE5Ojjh06JCIi4sTPj4+orS01Ob798Ybb4ht27aJixcvijNnzohHH31USKVSsWfPHlObnn4Mr/vDH/4gYmJiWn1OWzqGjzzyiFAqlSIlJcXse662ttbUxpbeE+0qjAghxNtvvy369u0rZDKZiI6OFkeOHDF9bdKkSWL+/Plm7bdu3SoGDRokZDKZGD58uNi+fbvZ141Go3j22WeFSqUScrlc3HHHHSIzM7M7utIqS/oXEhIiALT4WLlypRBCiNraWjFt2jTh6+srnJycREhIiPjTn/5ktV8O11nSx6VLl5raqlQqcffdd4sTJ06YPZ+tHUMhLP8+zcjIEADE999/3+K5bO04Xp/m+euP632aP3++mDRpUot9IiMjhUwmE/369RMffvhhi+e90f+z7mZpHydNmnTD9kI0T2cOCAgQMplMBAUFiYSEBJGVldW9HbvG0v69/PLLon///kKhUAhvb28xefJksXfv3hbP25OPoRDN01idnZ3Fe++91+pz2tIxbK1vAMx+tmzpPVFyrWgiIiIiq7CbMSNERETUMzGMEBERkVUxjBAREZFVMYwQERGRVTGMEBERkVUxjBAREZFVMYwQERGRVTGMEBERkVUxjBAREZFVMYwQERGRVTlau4D2MBqNKCoqgru7OyQSibXLISIionYQQkCn0yEwMBBSadvnP3pEGCkqKoJarbZ2GURERNQBBQUFCA4ObvPrPSKMuLu7A2jujIeHh5WrISIiovbQarVQq9Wm9/G29Igwcv3SjIeHB8MIERFRD3OzIRYcwEpERERWxTBCREREVsUwQkRERFbFMEJERNRLGY0CV2sacFGjg77JYLU6esQAViIiImqfRoMRV6obUF6tv/bRgCvX/n2lugFl1z6XV+tRUdOAJqMAAOz420QMC7TOJBGGESIiIhtnNApU1DZAo61HqU6PMq0epbp6lOmaw8Yvg0dVXaPFz690dkJNQ1MXVN4+DCNERERW0mQw4kpNA0q1elPQKNVd+3z9sbY5aFw/g9EeDlIJvF1l8HGTw8ft58993OS/2iaHt6sMMkfrjtpgGCEiIupkQghU1TWiuKoeJVX11z7XQaP9OWxotHpU1OjR3owhkQB9XGXwdVdA5SGHn7scvu7NgaLPtYDhe+3fns5OkEp7zu1TGEaIiIgsIIRARU3Dz0FD2xw0rj++Hj7qGts3INRBKoGPmwx+7gr4ucvh53H9sxx+puChQB83GZwc7HPeCcMIERHRL9Tom1BYWYfCq3W4fO2zKWxom4NGQ5OxXc/Vx1UGf6UCAUoFVB4K+HsoTCHD110OlYcC3q4yOPSgsxhdgWGEiIh6DSEEtHVNuFxZi8tXm4NGYWUdLl+tNQWQq7XtGwDq4yZHgFIBf6UCgUoF/JXOpsfXw4fCyaGLe2QfGEaIiMiuVNU2Ir+iFvkVtShsETrqUK2/+awRD4UjgrxcEOTpjGCv5pAR4HktbHg0Bw1rD/q0JwwjRETUoxiMAsVVdc2B40pz6MirqEVBRS3yrtS2a2prH1cZgr2cEeTljCBP52uhw6X5sZczPBRO3dATuo5hhIiIbE59owF5V2qRd6XGdJYj70pz4Lh8tQ4NhhuP2fB1l0Pt1Rwwfhk6gq+d7XCW8fKJLWEYISIiqzAYBQqv1iG7vBo55TXIKa9Bdlnz56KqOogbTHl1cpAg2MsFfb1dENKn+XNfbxf0vfZvFxnf3noSHi0iIuoyQgiUVeuRcy1k5JTXIPva5/wrtTc8w+GucERoH1dTwAj5RdgIUDr3+hko9oRhhIiIbpnBKFBQUYuLpdW4WKpDlqYaF0ubz3jcaMCozFGKsD6uCPNxRZhv8+d+Ps2fvV1lkEgYOHqDDoWRdevW4ZVXXkFJSQkiIiLw9ttvIzo6us32a9aswbvvvov8/Hz4+PjgvvvuQ1JSEhQKRYcLJyKi7tdoMCLvSi2ySnW4eC1wXCytxqWy6jbX3pBKgGAvl+bA4eOKftdCR5iPKwKVzj1qpVDqGhaHkS1btiAxMRHr169HTEwM1qxZg/j4eGRmZsLPz69F+02bNmHZsmX44IMPMG7cOFy4cAEPPfQQJBIJXn/99U7pBBERda5GgxHZZTW4eC10ZF0745FTXoNGQ+uDOeSOUgzwc8MAPzcMvPZ5gJ8b1N4ukDtywCi1TSLEjYYItRQTE4OxY8di7dq1AACj0Qi1Wo0lS5Zg2bJlLdovXrwY58+fR3Jysmnb448/jqNHj+LgwYPtek2tVgulUomqqip4eFjn9sZERPZICIEynR4ZJTpklGiRUazD+RIdskp1bYYOF5nDtbDhjoGq5uAx0M8dQV4cx0Hm2vv+bdGZkYaGBqSlpWH58uWmbVKpFHFxcUhNTW11n3HjxuHf//43jh07hujoaGRnZ2PHjh2YO3dum6+j1+uh1+vNOkNERLemvtGAi5pqnL8WOjJKtMgo0aGipqHV9u5yRwzydzc7yzFQ5Y5ApYJjOahTWRRGysvLYTAYoFKpzLarVCpkZGS0us8DDzyA8vJyTJgwAUIINDU14a9//SuefvrpNl8nKSkJq1atsqQ0IiL6hcraBpwt0iK9sArpRVqcK6pCTnlNq3eIlUqAMB9XDAnwwFB/dwzx98CQAHcEeTozdFC36PLZNCkpKXjppZfwzjvvICYmBllZWXj00Ufxwgsv4Nlnn211n+XLlyMxMdH0WKvVQq1Wd3WpREQ9UplOj/SiKpwtrEJ6oRbpRVW4fLWu1bZeLk4YGuBhChxD/T0wUOXGe6iQVVkURnx8fODg4ACNRmO2XaPRwN/fv9V9nn32WcydOxcLFy4EAISHh6OmpgZ//vOf8Y9//ANSacu1/eVyOeRyuSWlERHZPSEESrT1OHO5+WzH2cIqpBdVQaPVt9o+pI8LRgQqMTzIA8MCmj983eU820E2x6IwIpPJEBUVheTkZMyaNQtA8wDW5ORkLF68uNV9amtrWwQOB4fmBG7h2Fkiol5FW9+IM5ercKqg0vRRpmsZPCQSoL+vG0YEemBEkBLDA5UYFugBpTPvr0I9g8WXaRITEzF//nyMGTMG0dHRWLNmDWpqarBgwQIAwLx58xAUFISkpCQAwIwZM/D6669j1KhRpss0zz77LGbMmGEKJUREvV2jwYiMYh1OXa7EqfxK/HS5EpfKqlssie4glWCQyt0UPEYENV9ycZVzDUvquSz+7k1ISEBZWRlWrFiBkpISREZGYufOnaZBrfn5+WZnQp555hlIJBI888wzKCwshK+vL2bMmIEXX3yx83pBRNSDCCFw+WodThY0B49TBVdxtkgLfSuLhqm9nRER7IlIdfPHiCAlx3eQ3bF4nRFr4DojRNSTNRqMOFekxfG8q0jLq8Dx3KsobeVyi4fCERFqT4xSeyKyrydGBnvCx43j56jn6pJ1RoiI6Oa09Y04kXcVx3Ov4nheBX4qqEJdo8GsjaNUguGBHs1nPPp6IiLYE2E+rhxcSr0SwwgR0S0q1dbjSE4FjuVcwfHcq8jU6FqM9fBQOCIqxAtjQr0RFeKFiGBPOMt4uYUIYBghIrJYSVU9juZcwZHsChzNvoLs8poWbUL6uDSHjxBvjAn1wgBfN94QjqgNDCNERDdRXFWHI9lXcDS7AkeyryD3Sq3Z1yUSYFiAB2LC+iA6zAujQ7zg5867khO1F8MIEdGvXK1pwOFLV3AwqwyHL11B3q/Ch1QCDA9UIibMG7f164Oxod5QunBND6KOYhghol6vvtGA47lXcTCrHAezynC2SGs25kMqAcKDlIjp1we39fPGmFBveCgYPog6C8MIEfU6RqPAuWItfrhYjkNZ5fgxt6LFGh+DVe6YMNAH4wc0n/lwZ/gg6jIMI0TUK5Tp9DhwoQz7MktxKKscV2sbzb6u8pBjwgBfTBjYB+P7+8DPg2M+iLoLwwgR2SWDUeCny5VIyShFyoUynL5cZfZ1N7kjbuvnjQkDfDBhoA/6+7pxjQ8iK2EYISK7UVHTgAMXypCSWYr9F8panP0YEeSBKYP9MGmQLyLUnnByaHnXcCLqfgwjRNRjCSFwvliHPec12JdZilMFlWYDT90Vjrh9oC8mD/bFpMG+nG5LZKMYRoioR2kyGHEstwK7z2mw+5wGl6/WmX19iL87pgzxw5TBfhjVl2c/iHoChhEisnk1+iYcuFCG3ec0SM4oRVXdz5df5I5STBzoizuG+mHyYF8EKJ2tWCkRdQTDCBHZpFJdPZLPl2L3OQ0OZpWj4RdTb71cnHDHUBWmDVNh4kBf3uOFqIdjGCEim6HR1uO7M8XYcaYEP+ZVmI3/6OvtgmnDVJg6TIWoEC848vILkd1gGCEiqyrV1uO79BJsP1OMH3PNA8jIYOW1AOKPQSpOvSWyVwwjRNTtSnX12Jlegu2ni3HsVwFkVF9PTA8PwN3hAQj05PgPot6AYYSIukWZTo+d6cXYfqYYR3NaDyB3hQcgiAGEqNdhGCGiLlOjb8L350qw7WQRDmaVw2D8OYFEqj1xz0gGECJiGCGiTtZoMOLgxXJsO1WI789qUNdoMH0tIliJe0YG4q5wfwR7uVixSiKyJQwjRHTLhBA4VVCJbScL8e3pYlypaTB9LbSPC2aNCsLMyCCE+bhasUoislUMI0TUYcVVdfjyRCE+T7uMnPIa0/Y+rjLMiAjErFFBiAhWchYMEd0QwwgRWaS+0YA95zXYevwyDl4sw/VhIM5ODogfrsLMUUGYMMCHy7ATUbt16LfFunXrEBoaCoVCgZiYGBw7duyG7SsrK7Fo0SIEBARALpdj0KBB2LFjR4cKJqLuJ4TAmctVWPF1OmJeSsbiTSdx4EJzEIkO88Yr943E8WfisGb2KEwZ7McgQkQWsfjMyJYtW5CYmIj169cjJiYGa9asQXx8PDIzM+Hn59eifUNDA6ZOnQo/Pz98/vnnCAoKQl5eHjw9PTujfiLqQpW1DfjyRCG2Hi9ARonOtD1AqcDvRgfjvqhghHIcCBHdIokQv5ztf3MxMTEYO3Ys1q5dCwAwGo1Qq9VYsmQJli1b1qL9+vXr8corryAjIwNOTk4dKlKr1UKpVKKqqgoeHh4deg4iah8hBNLyrmLT0XxsP1MM/bV7wsgcpYgf7o/fRwVj/AAfOEg5DoSIbqy9798WnRlpaGhAWloali9fbtomlUoRFxeH1NTUVvf55ptvEBsbi0WLFuHrr7+Gr68vHnjgATz11FNwcGj95lZ6vR56vd6sM0TUtarqGrHtZCE2Hc1HpubnsyBDAzwwJ1qNmRFBULp07A8KIqIbsSiMlJeXw2AwQKVSmW1XqVTIyMhodZ/s7Gzs3bsXDz74IHbs2IGsrCz87//+LxobG7Fy5cpW90lKSsKqVassKY2IOuD6lNxNR/Px39NFqG9sPguicJJixshAPBDTF5FqT86GIaIu1eWzaYxGI/z8/PDee+/BwcEBUVFRKCwsxCuvvNJmGFm+fDkSExNNj7VaLdRqdVeXStRr1Dca8M2pInyUmouzRT+feRykcsMD0X1x7+hgKJ15FoSIuodFYcTHxwcODg7QaDRm2zUaDfz9/VvdJyAgAE5OTmaXZIYOHYqSkhI0NDRAJpO12Ecul0Mul1tSGhG1Q1FlHT45kofNx/JxtbYRQPNYkHvCA/BATF9EhXjxLAgRdTuLwohMJkNUVBSSk5Mxa9YsAM1nPpKTk7F48eJW9xk/fjw2bdoEo9EIqbR5ut+FCxcQEBDQahAhos4lhMCxnApsPJyL789pTPeHCfJ0xtzYECSMUcPLlT+LRGQ9Fl+mSUxMxPz58zFmzBhER0djzZo1qKmpwYIFCwAA8+bNQ1BQEJKSkgAAjzzyCNauXYtHH30US5YswcWLF/HSSy/hb3/7W+f2hIjM1Dca8PWpQnx4KNdsWm5svz54aHwo4oaqOCOGiGyCxWEkISEBZWVlWLFiBUpKShAZGYmdO3eaBrXm5+ebzoAAgFqtxq5du/DYY49h5MiRCAoKwqOPPoqnnnqq83pBRCZlOj0+Sc3FJ0fyTJdiFE5S3DsqGPPHhWCIP6fHE5FtsXidEWvgOiNEN5dVWo33D2bjixOFaLi2NkiwlzPmxYbg/jFqeLrwUgwRda8uWWeEiGyLEAI/5l7FewcuYc/5UtP2CLUn/nJ7P8QP9+elGCKyeQwjRD1Qk8GIXWc1eO+HbPxUUAkAkEiAuKEq/Pn2fhjDWTFE1IMwjBD1IPomAz5Pu4x/7c9GfkUtgOapub8bHYyFE8PQ39fNyhUSEVmOYYSoB6htaMKmo/nY8EM2NNrmWyV4uThhbmwo5sWGwMeN6/IQUc/FMEJkw6rqGvFJai4+OJSLipoGAM13zP3z7f0we2xfOMtav78TEVFPwjBCZIOuVOvxwaEcfHw4Dzp9EwCgr7cL/ndyf9w7OghyR4YQIrIfDCNENqRUV49/7c/GpqP5qGs0AAAG+rlh0ZQBuGdkABwdpDd5BiKinodhhMgGXKnW418HsvFxaq7pzrnhQUosmjIA04apIOX0XCKyYwwjRFZ0taYBG37IxsbDuahtaD4TEqn2xNK4gZg0yJfTc4moV2AYIbKCqrpGvP9DNj44lIvqa2NCRgYr8djUQZjMEEJEvQzDCFE30tU34sNDudjwQzZ09c0hZGiABxKnDkLcUD+GECLqlRhGiLpBfaMBn6TmYV1KFiqv3bxukMoNj8UNQvxwf44JIaJejWGEqAsZjAJfnSzE699noqiqHgDQz9cVS+MG4Z7wAIYQIiIwjBB1CSEE9mWW4uXvMpGp0QFoXqzssamD8LvRwbx5HRHRLzCMEHWyE/lXsfq7DBzLqQAAKJ2dsGhKf8yLDYXCiYuVERH9GsMIUSfJKq3GK7sysOusBgAgd5RiwfgwPDKpP5QuTlaujojIdjGMEN2iipoGvLH7AjYdy4fBKCCVAL+PUmPp1IEIUDpbuzwiIpvHMELUQQ1NRnycmos3ky+apunGDVXhqTsHY6DK3crVERH1HAwjRBYSQmD3OQ1e2nEeuVdqAQDDAz3w7D3DcFu/Plaujoio52EYIbLAuSItXvj2HFKzrwAAfN3leHLaYPwuijNkiIg6imGEqB3KdHq89n0mthwvgBCAzFGKP00MwyOTB8BNzh8jIqJbwd+iRDfQaDDiw0M5eCs5y3QPmXtGBuCpO4dA7e1i5eqIiOwDwwhRGw5nlWPFN2eRVVoNAIgIVuLZe4ZhTKi3lSsjIrIvDCNEv1JcVYf/t/08tp8uBgB4u8qw7K4huG90MJdvJyLqAtKO7LRu3TqEhoZCoVAgJiYGx44da9d+mzdvhkQiwaxZszryskRdqqHJiPX7L+GO1/Zj++liSCXAvNgQ7Ht8Mu4fo2YQISLqIhafGdmyZQsSExOxfv16xMTEYM2aNYiPj0dmZib8/Pza3C83NxdPPPEEJk6ceEsFE3WFgxfLsfKbdFwqqwEARIV44fmZwzE8UGnlyoiI7J9ECCEs2SEmJgZjx47F2rVrAQBGoxFqtRpLlizBsmXLWt3HYDDg9ttvxx//+Ef88MMPqKysxLZt29p8Db1eD71eb3qs1WqhVqtRVVUFDw8PS8oluqHiqjr8v2/PY/uZ5ksyfVxlWH73UPx2VBDPhBAR3SKtVgulUnnT92+LLtM0NDQgLS0NcXFxPz+BVIq4uDikpqa2ud/zzz8PPz8/PPzww+16naSkJCiVStOHWq22pEyimzIYBTYeysHU1w9g+5nmSzIPjQvF3icm474ojg0hIupOFl2mKS8vh8FggEqlMtuuUqmQkZHR6j4HDx7E+++/j1OnTrX7dZYvX47ExETT4+tnRog6w9miKjz95Rn8dLkKADCqrydenBWOYYE860ZEZA1dOptGp9Nh7ty52LBhA3x8fNq9n1wuh1wu78LKqDeqbWjCmj0X8f7BHBiMAu5yR/z9riF4MLovz4QQEVmRRWHEx8cHDg4O0Gg0Zts1Gg38/f1btL906RJyc3MxY8YM0zaj0dj8wo6OyMzMRP/+/TtSN5FF9mWW4tlt6bh8tQ4AMD08ACtnDIOfh8LKlRERkUVhRCaTISoqCsnJyabpuUajEcnJyVi8eHGL9kOGDMGZM2fMtj3zzDPQ6XR48803eemFulyprh7P//ccvr22ZkiQpzOenzkcdwxV3WRPIiLqLhZfpklMTMT8+fMxZswYREdHY82aNaipqcGCBQsAAPPmzUNQUBCSkpKgUCgwYsQIs/09PT0BoMV2os4khMCWHwvw0o7z0NY3QSoB/jg+DI9NHQRX3kuGiMimWPxbOSEhAWVlZVixYgVKSkoQGRmJnTt3mga15ufnQyrt0FpqRJ2ioKIWy748jUNZzXfWDQ9SIum34RgRxDVDiIhskcXrjFhDe+cpU+9mNAp8ciQPL+/MQG2DAXJHKZ6YNhgLxofC0YEBmYiou7X3/Zvnq8kuZJdV46kvTuPH3KsAgOgwb7z8u5EI83G1cmVERHQzDCPUoxmMAu8fzMZr31+AvskIF5kDlt81BA/GhHC6LhFRD8EwQj3WBY0OT35+Gj8VVAIAJg70wUv3hkPt7WLdwoiIyCIMI9TjNBmM+NeBbKzZcwGNBgF3hSOenT4Mvx8TDImEZ0OIiHoahhHqUbLLqpG49SecunY25I4hfnjx3nD4K7l4GRFRT8UwQj2C0SjwcWouVu/MQH2jEe5yRzz3P8Px29FBPBtCRNTDMYyQzSuqrMOTn/9kWjdkwgAf/PO+kQj0dLZyZURE1BkYRshmCSHwxYlCrPrmLHT6JiicpHj67qH4A2fKEBHZFYYRsknl1Xo8/eUZfH+u+aaMo/p64rXfR6Cfr5uVKyMios7GMEI2Z2d6Cf7x1RlcqWmAk4MES+MG4S+39+MqqkREdophhGxGVV0jVn1zFl+eLAQADPF3x+v3R2JYIG8BQERkzxhGyCYcyb6CxC2nUFRVD6kE+Muk/lgaNxByRwdrl0ZERF2MYYSsqtFgxJo9F/BOyiUIAYT0ccHr90cgKsTb2qUREVE3YRghq8m7UoO/bT5lWs79/jHBWDljOFzl/LYkIupN+Fufut31Kbsrv05HTYMBHgpHJP12JKaPDLB2aUREZAUMI9Stquoa8cy2dPz3pyIAQHSYN95IiEQQFzAjIuq1GEao2/yYW4Glm0+hsLIODlIJHosbiEcmD4ADFzAjIurVGEaoyzUZjHhrbxbW7r0IowD6ertgzexIjO7rZe3SiIjIBjCMUJcqqKjFo5tP4kR+JQDgt6ODsOp/hsNd4WTdwoiIyGYwjFCX2XayEM9uS4dO3wR3uSP+370jMDMyyNplERGRjWEYoU6nq2/Eiq/P4qtrK6lGhXhhTUIk1N4uVq6MiIhsEcMIdaq0vKtYuuUkCirqIJUAf7tjIBZPGcD7yhARUZsYRqhTGIwC6/Zl4c3kizAYBYI8nfHm7EiMCeVKqkREdGMd+nN13bp1CA0NhUKhQExMDI4dO9Zm2w0bNmDixInw8vKCl5cX4uLibtieep7LV2sx+71UvL77AgxGgf+JCMR3SycyiBARUbtYHEa2bNmCxMRErFy5EidOnEBERATi4+NRWlraavuUlBTMmTMH+/btQ2pqKtRqNaZNm4bCwsJbLp6s778/FeGuN3/Aj7lX4SpzwOv3R+DN2ZHw4GwZIiJqJ4kQQliyQ0xMDMaOHYu1a9cCAIxGI9RqNZYsWYJly5bddH+DwQAvLy+sXbsW8+bNa7WNXq+HXq83PdZqtVCr1aiqqoKHB28nbwuq9U1Y+fVZfHHiMgAgUu2JN2dHIqSPq5UrIyIiW6HVaqFUKm/6/m3RmZGGhgakpaUhLi7u5yeQShEXF4fU1NR2PUdtbS0aGxvh7d32KfykpCQolUrTh1qttqRM6mKnCiox/a0f8MWJy5BIgCW/GYDP/hrLIEJERB1iURgpLy+HwWCASqUy265SqVBSUtKu53jqqacQGBhoFmh+bfny5aiqqjJ9FBQUWFImdZHrg1Tve/cw8q7UIlCpwOY/3YbHpw2GE2fLEBFRB3XrbJrVq1dj8+bNSElJgUKhaLOdXC6HXC7vxsroZooq6/DYllM4mlMBAJg+MgAvzQqH0oVjQ4iI6NZYFEZ8fHzg4OAAjUZjtl2j0cDf3/+G+7766qtYvXo19uzZg5EjR1peKVnN9tPFePqrM6iqa4SLzAHP/c9w/D4qGBIJb3BHRES3zqJz6zKZDFFRUUhOTjZtMxqNSE5ORmxsbJv7/fOf/8QLL7yAnTt3YsyYMR2vlrpVtb4JT372ExZtOoGqukaMDFZi+98m4v4xagYRIiLqNBZfpklMTMT8+fMxZswYREdHY82aNaipqcGCBQsAAPPmzUNQUBCSkpIAAC+//DJWrFiBTZs2ITQ01DS2xM3NDW5ubp3YFepMJ/OvYumWU8i7UguJBPjfyf2xNG4Qx4YQEVGnsziMJCQkoKysDCtWrEBJSQkiIyOxc+dO06DW/Px8SKU/v2G9++67aGhowH333Wf2PCtXrsRzzz13a9VTpzMYBd7Zl4U111ZSDVQq8HpCJG7r18fapRERkZ2yeJ0Ra2jvPGW6NQUVtUjcego/5l4FANwzMgAv3hsOpTMHqRIRkeXa+/7Ne9MQAGDbyUI8uy0dOn0T3OSOeH7mcNw7KohjQ4iIqMsxjPRy2vpGPLstHV+fKgIAjO7riTUJo9C3j4uVKyMiot6CYaQX+zG3Aks3n0JhZR2kEuBvdwzE4ikD4MhBqkRE1I0YRnohfZMBa/ZcxL/2X4JRAGpvZ6xJGIWoEC9rl0ZERL0Qw0gvk15Yhce3/oRMjQ4A8NvRQVj1P8PhzrvsEhGRlTCM9BJNBiPeSbmEt5Ivosko0MdVhhfvDcedI268ci4REVFXYxjpBbJKdXh860/46XIVAODO4f548d4R6OPG+/8QEZH1MYzYMYNR4MNDOfjnrkw0NBnhoXDE8zNHYGZkIKfsEhGRzWAYsVPni7VY9sVp09mQSYN88fLvRsJf2fbdkomIiKyBYcTO1DcasHZvFtbvv4Qmo4C7whFP3z0Us8fy5nZERGSbGEbsyNHsK1j+5Rlkl9cAAOKHq/D8zBFQefBsCBER2S6GETugrW/E6u8ysOloPgDAz12O52cOx50jAqxcGRER0c0xjPRgQgh8dbIQSd9loEynBwDMie6LZXcN4c3tiIiox2AY6aHOFlVh5ddncTyv+Q67/Xxc8dJvw3Fbvz5WroyIiMgyDCM9zNWaBryx5wL+fSQPRgE4OzlgyR0D8PCEMMgdHaxdHhERkcUYRnqI+kYDPjyUi3dSsqCrbwIA3DMyAP+YPhQBSmcrV0dERNRxDCM2zmBsHhfy2veZKK6qBwAMDfDAs/cMxbj+PlaujoiI6NYxjNgog1Hgu/RivJ2cZbqpXaBSgSfiB2NWZBCkUq4ZQkRE9oFhxMYYjALfni7C23uzkFVaDQBwVzhi8ZQBmD8uFAonjgshIiL7wjBiI2r0TfjyZCE+PJhjWrTMQ+GIP04Iw4JxYVC6cKouERHZJ4YRK8str8HHqXn4LK3ANDDV08UJf5rYD3NjQ+ChYAghIiL7xjBiBbr6Rnx/VoNtpwrxw8Vy0/YwH1fMvS0E949Vw03OQ0NERL0D3/G6SW1DEw5cKMd/fyrCnvMa6JuMAACJBJg8yBfzx4Xi9oG+HJhKRES9jrQjO61btw6hoaFQKBSIiYnBsWPHbtj+s88+w5AhQ6BQKBAeHo4dO3Z0qNiexGAU+KmgEuv2ZWH2e6mIXLUbf/13GrafKYa+yYh+vq54LG4QUp6YjA8XRGPyYD8GESIi6pUsPjOyZcsWJCYmYv369YiJicGaNWsQHx+PzMxM+Pn5tWh/+PBhzJkzB0lJSbjnnnuwadMmzJo1CydOnMCIESM6pRPWJIRAVV0j8itqcb5Yi3NFWpwr1uJ8sQ7V+iaztkGezpg+MgD/ExGI4YEekEgYPoiIiCRCCGHJDjExMRg7dizWrl0LADAajVCr1ViyZAmWLVvWon1CQgJqamrw7bffmrbddtttiIyMxPr169v1mlqtFkqlElVVVfDw8LCk3BvadbYEZTo9BAAIAaNoDhdGAQg0/1sIoMkoUNfQhJoGA2obDKjRN6G8Wo+SqnoUV9WjrtHQ6vO7yx1xW/8+mDjQBxMH+iK0jwsDCBER9Rrtff+26MxIQ0MD0tLSsHz5ctM2qVSKuLg4pKamtrpPamoqEhMTzbbFx8dj27Ztbb6OXq+HXq83PdZqtZaU2W7r91/CyfzKTnmuPq4yDAlwx7AADwwL9MCwACX6+7rC0aFDV8KIiIh6DYvCSHl5OQwGA1Qqldl2lUqFjIyMVvcpKSlptX1JSUmbr5OUlIRVq1ZZUlqHxIT1gZ+7HBJIIJEAUokEaP4PEokE0mv/lkolcJU5wlXuCFeZA5xlDvBxk8NfqUCAUgGVh4KLkREREXWQTc6mWb58udnZFK1WC7Va3emvs+yuIZ3+nERERGQZi8KIj48PHBwcoNFozLZrNBr4+/u3uo+/v79F7QFALpdDLpdbUhoRERH1UBYNaJDJZIiKikJycrJpm9FoRHJyMmJjY1vdJzY21qw9AOzevbvN9kRERNS7WHyZJjExEfPnz8eYMWMQHR2NNWvWoKamBgsWLAAAzJs3D0FBQUhKSgIAPProo5g0aRJee+01TJ8+HZs3b8bx48fx3nvvdW5PiIiIqEeyOIwkJCSgrKwMK1asQElJCSIjI7Fz507TINX8/HxIpT+fcBk3bhw2bdqEZ555Bk8//TQGDhyIbdu22cUaI0RERHTrLF5nxBqqqqrg6emJgoKCTl1nhIiIiLrO9QkolZWVUCqVbbazydk0v6bT6QCgS2bUEBERUdfS6XQ3DCM94syI0WhEUVER3N3dO3UF0+uJzZ7PuLCPPZ+99w9gH+2BvfcPsP8+dkX/hBDQ6XQIDAw0G8Lxaz3izIhUKkVwcHCXPb+Hh4ddfmP9EvvY89l7/wD20R7Ye/8A++9jZ/fvRmdEruNa5URERGRVDCNERERkVb06jMjlcqxcudKuV3tlH3s+e+8fwD7aA3vvH2D/fbRm/3rEAFYiIiKyX736zAgRERFZH8MIERERWRXDCBEREVkVwwgRERFZFcMIERERWZXdhZF169YhNDQUCoUCMTExOHbs2A3bf/bZZxgyZAgUCgXCw8OxY8cOs68LIbBixQoEBATA2dkZcXFxuHjxYld24YYs6d+GDRswceJEeHl5wcvLC3FxcS3aP/TQQ5BIJGYfd955Z1d344Ys6ePGjRtb1K9QKMza2NoxBCzr4+TJk1v0USKRYPr06aY2tnQcDxw4gBkzZiAwMBASiQTbtm276T4pKSkYPXo05HI5BgwYgI0bN7ZoY+nPdleytI9ffvklpk6dCl9fX3h4eCA2Nha7du0ya/Pcc8+1OIZDhgzpwl60zdL+paSktPo9WlJSYtauJx/D1n7GJBIJhg8fbmpjS8cwKSkJY8eOhbu7O/z8/DBr1ixkZmbedD9rvSfaVRjZsmULEhMTsXLlSpw4cQIRERGIj49HaWlpq+0PHz6MOXPm4OGHH8bJkycxa9YszJo1C+np6aY2//znP/HWW29h/fr1OHr0KFxdXREfH4/6+vru6paJpf1LSUnBnDlzsG/fPqSmpkKtVmPatGkoLCw0a3fnnXeiuLjY9PHpp592R3daZWkfgeali39Zf15entnXbekYApb38csvvzTrX3p6OhwcHPD73//erJ2tHMeamhpERERg3bp17Wqfk5OD6dOnY8qUKTh16hSWLl2KhQsXmr1Zd+T7oitZ2scDBw5g6tSp2LFjB9LS0jBlyhTMmDEDJ0+eNGs3fPhws2N48ODBrij/pizt33WZmZlm9fv5+Zm+1tOP4ZtvvmnWt4KCAnh7e7f4ObSVY7h//34sWrQIR44cwe7du9HY2Ihp06ahpqamzX2s+p4o7Eh0dLRYtGiR6bHBYBCBgYEiKSmp1fb333+/mD59utm2mJgY8Ze//EUIIYTRaBT+/v7ilVdeMX29srJSyOVy8emnn3ZBD27M0v79WlNTk3B3dxcfffSRadv8+fPFzJkzO7vUDrO0jx9++KFQKpVtPp+tHUMhbv04vvHGG8Ld3V1UV1ebttnacbwOgPjqq69u2Obvf/+7GD58uNm2hIQEER8fb3p8q//PulJ7+tiaYcOGiVWrVpker1y5UkRERHReYZ2kPf3bt2+fACCuXr3aZht7O4ZfffWVkEgkIjc317TNVo+hEEKUlpYKAGL//v1ttrHme6LdnBlpaGhAWloa4uLiTNukUini4uKQmpra6j6pqalm7QEgPj7e1D4nJwclJSVmbZRKJWJiYtp8zq7Skf79Wm1tLRobG+Ht7W22PSUlBX5+fhg8eDAeeeQRXLlypVNrb6+O9rG6uhohISFQq9WYOXMmzp49a/qaLR1DoHOO4/vvv4/Zs2fD1dXVbLutHEdL3eznsDP+n9kao9EInU7X4mfx4sWLCAwMRL9+/fDggw8iPz/fShV2TGRkJAICAjB16lQcOnTItN0ej+H777+PuLg4hISEmG231WNYVVUFAC2+537Jmu+JdhNGysvLYTAYoFKpzLarVKoW1y2vKykpuWH7658tec6u0pH+/dpTTz2FwMBAs2+kO++8Ex9//DGSk5Px8ssvY//+/bjrrrtgMBg6tf726EgfBw8ejA8++ABff/01/v3vf8NoNGLcuHG4fPkyANs6hsCtH8djx44hPT0dCxcuNNtuS8fRUm39HGq1WtTV1XXK976tefXVV1FdXY3777/ftC0mJgYbN27Ezp078e677yInJwcTJ06ETqezYqXtExAQgPXr1+OLL77AF198AbVajcmTJ+PEiRMAOuf3ly0pKirCd9991+Ln0FaPodFoxNKlSzF+/HiMGDGizXbWfE90vKW9qcdYvXo1Nm/ejJSUFLMBnrNnzzb9Ozw8HCNHjkT//v2RkpKCO+64wxqlWiQ2NhaxsbGmx+PGjcPQoUPxr3/9Cy+88IIVK+sa77//PsLDwxEdHW22vacfx95k06ZNWLVqFb7++muzMRV33XWX6d8jR45ETEwMQkJCsHXrVjz88MPWKLXdBg8ejMGDB5sejxs3DpcuXcIbb7yBTz75xIqVdY2PPvoInp6emDVrltl2Wz2GixYtQnp6utXGr7SH3ZwZ8fHxgYODAzQajdl2jUYDf3//Vvfx9/e/Yfvrny15zq7Skf5d9+qrr2L16tX4/vvvMXLkyBu27devH3x8fJCVlXXLNVvqVvp4nZOTE0aNGmWq35aOIXBrfaypqcHmzZvb9UvNmsfRUm39HHp4eMDZ2blTvi9sxebNm7Fw4UJs3bq1xenwX/P09MSgQYN6xDFsTXR0tKl2ezqGQgh88MEHmDt3LmQy2Q3b2sIxXLx4Mb799lvs27cPwcHBN2xrzfdEuwkjMpkMUVFRSE5ONm0zGo1ITk42+8v5l2JjY83aA8Du3btN7cPCwuDv72/WRqvV4ujRo20+Z1fpSP+A5pHPL7zwAnbu3IkxY8bc9HUuX76MK1euICAgoFPqtkRH+/hLBoMBZ86cMdVvS8cQuLU+fvbZZ9Dr9fjDH/5w09ex5nG01M1+Djvj+8IWfPrpp1iwYAE+/fRTs2nZbamursalS5d6xDFszalTp0y128sxBJpnqWRlZbXrjwJrHkMhBBYvXoyvvvoKe/fuRVhY2E33sep74i0Nf7UxmzdvFnK5XGzcuFGcO3dO/PnPfxaenp6ipKRECCHE3LlzxbJly0ztDx06JBwdHcWrr74qzp8/L1auXCmcnJzEmTNnTG1Wr14tPD09xddffy1Onz4tZs6cKcLCwkRdXZ3N92/16tVCJpOJzz//XBQXF5s+dDqdEEIInU4nnnjiCZGamipycnLEnj17xOjRo8XAgQNFfX19t/evI31ctWqV2LVrl7h06ZJIS0sTs2fPFgqFQpw9e9bUxpaOoRCW9/G6CRMmiISEhBbbbe046nQ6cfLkSXHy5EkBQLz++uvi5MmTIi8vTwghxLJly8TcuXNN7bOzs4WLi4t48sknxfnz58W6deuEg4OD2Llzp6nNzf6fdTdL+/if//xHODo6inXr1pn9LFZWVpraPP744yIlJUXk5OSIQ4cOibi4OOHj4yNKS0ttvn9vvPGG2LZtm7h48aI4c+aMePTRR4VUKhV79uwxtenpx/C6P/zhDyImJqbV57SlY/jII48IpVIpUlJSzL7namtrTW1s6T3RrsKIEEK8/fbbom/fvkImk4no6Ghx5MgR09cmTZok5s+fb9Z+69atYtCgQUImk4nhw4eL7du3m33daDSKZ599VqhUKiGXy8Udd9whMjMzu6MrrbKkfyEhIQJAi4+VK1cKIYSora0V06ZNE76+vsLJyUmEhISIP/3pT1b75XCdJX1cunSpqa1KpRJ33323OHHihNnz2doxFMLy79OMjAwBQHz//fctnsvWjuP1aZ6//rjep/nz54tJkya12CcyMlLIZDLRr18/8eGHH7Z43hv9P+tulvZx0qRJN2wvRPN05oCAACGTyURQUJBISEgQWVlZ3duxayzt38svvyz69+8vFAqF8Pb2FpMnTxZ79+5t8bw9+RgK0TyN1dnZWbz33nutPqctHcPW+gbA7GfLlt4TJdeKJiIiIrIKuxkzQkRERD0TwwgRERFZFcMIERERWRXDCBEREVkVwwgRERFZFcMIERERWRXDCBEREVkVwwgRERFZFcMIERERWRXDCBEREVkVwwgRERFZ1f8HP3BRYZzvHJYAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "delta_t = 0.01\n", + "x_max = 2\n", + "x = [i * delta_t for i in range(int(x_max / delta_t))]\n", + "y_pdf = [pdf(xx) for xx in x]\n", + "y_cdf = [cdf(xx) for xx in x]\n", + "fig, ax = plt.subplots(2, 1)\n", + "ax[0].plot(x, y_pdf)\n", + "ax[1].plot(x, y_cdf)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "adeb977a-4f70-4773-a98a-02865b437d67", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n" + ] + } + ], + "source": [ + "base_parameters = Lognormal.parametrizations.get_base_parameters(dist.parameters)\n", + "y_true_pdf = [lognormal_pdf(base_parameters, xx) for xx in x]\n", + "print(max([y_true_pdf[i] - y_pdf[i] for i in range(len(x))]))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index c19737f..63c7ffd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,6 +55,10 @@ coverage = { version = ">=7.5", extras = [ "toml" ] } pre-commit = ">=3.6" types-setuptools = "*" +[tool.poetry.group.docs.dependencies] +jupyter = "^1.1.1" +matplotlib = "^3.10.6" + [tool.ruff] target-version = "py312" line-length = 100 diff --git a/src/pysatl_core/families/__init__.py b/src/pysatl_core/families/__init__.py new file mode 100644 index 0000000..b22cbf0 --- /dev/null +++ b/src/pysatl_core/families/__init__.py @@ -0,0 +1,35 @@ +""" +Parametric Families module for working with statistical distribution families. + +This package provides a comprehensive framework for defining, managing, and +working with parametric families of statistical distributions. It supports +multiple parameterizations, constraint validation, and automatic conversion +between different parameter formats. +""" + +__author__ = "Leonid Elkin, Mikhail, Mikhailov" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + + +from pysatl_core.families.distribution import ParametricFamilyDistribution +from pysatl_core.families.parametric_family import ParametricFamily +from pysatl_core.families.parametrizations import ( + Parametrization, + ParametrizationConstraint, + ParametrizationSpec, + constraint, + parametrization, +) +from pysatl_core.families.registry import ParametricFamilyRegister + +__all__ = [ + "ParametricFamilyRegister", + "ParametrizationConstraint", + "Parametrization", + "ParametrizationSpec", + "ParametricFamily", + "ParametricFamilyDistribution", + "constraint", + "parametrization", +] diff --git a/src/pysatl_core/families/distribution.py b/src/pysatl_core/families/distribution.py new file mode 100644 index 0000000..750c33e --- /dev/null +++ b/src/pysatl_core/families/distribution.py @@ -0,0 +1,151 @@ +""" +Concrete distribution instances with specific parameter values. + +This module provides the implementation for individual distribution instances +created from parametric families. It handles distribution characteristics +computation, sampling, and provides access to analytical methods for +specific parameter sets. +""" + +from __future__ import annotations + +__author__ = "Leonid Elkin, Mikhail, Mikhailov" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + +from collections.abc import Mapping +from dataclasses import dataclass +from functools import partial +from typing import TYPE_CHECKING, Any + +from pysatl_core.distributions import ( + AnalyticalComputation, + ComputationStrategy, + Sample, + SamplingStrategy, +) +from pysatl_core.families.parametrizations import Parametrization +from pysatl_core.families.registry import ParametricFamilyRegister +from pysatl_core.types import ( + DistributionType, + GenericCharacteristicName, +) + +if TYPE_CHECKING: + from pysatl_core.families.parametric_family import ParametricFamily + + +@dataclass +class ParametricFamilyDistribution: + """ + A specific distribution instance from a parametric family. + + This class represents a concrete distribution with specific parameter + values, providing methods for computation and sampling. + + Attributes + ---------- + distr_name : str + Name of the distribution family. + distribution_type : DistributionType + Type of this distribution. + parameters : Parametrization + Parameter values for this distribution. + """ + + distr_name: str + distribution_type: DistributionType + parameters: Parametrization + + @property + def family(self) -> ParametricFamily: + """ + Get the parametric family this distribution belongs to. + + Returns + ------- + ParametricFamily + The parametric family of this distribution. + """ + return ParametricFamilyRegister.get(self.distr_name) + + @property + def analytical_computations( + self, + ) -> Mapping[GenericCharacteristicName, AnalyticalComputation[Any, Any]]: + """ + Get analytical computation functions for this distribution. + + Returns + ------- + Mapping[GenericCharacteristicName, AnalyticalComputation] + Mapping from characteristic names to computation functions. + """ + analytical_computations = {} + + # First form list of all characteristics, available from current parametrization + for characteristic, forms in self.family.distr_characteristics.items(): + if self.parameters.name in forms: + analytical_computations[characteristic] = AnalyticalComputation( + target=characteristic, + func=partial(forms[self.parameters.name], self.parameters), + ) + # TODO: Second, apply rule set, for, e.g. approximations + + # Finally, fill other chacteristics + base_name = self.family.parametrizations.base_parametrization_name + base_parameters = self.family.parametrizations.get_base_parameters(self.parameters) + for characteristic, forms in self.family.distr_characteristics.items(): + if characteristic in analytical_computations: + continue + if base_name in forms: + analytical_computations[characteristic] = AnalyticalComputation( + target=characteristic, func=partial(forms[base_name], base_parameters) + ) + + return analytical_computations + + @property + def sampling_strategy(self) -> SamplingStrategy: + """ + Get the sampling strategy for this distribution. + + Returns + ------- + SamplingStrategy + Strategy for sampling from this distribution. + """ + return self.family.sampling_strategy + + @property + def computation_strategy(self) -> ComputationStrategy[Any, Any]: + """ + Get the computation strategy for this distribution. + + Returns + ------- + ComputationStrategy + Strategy for computing characteristics of this distribution. + """ + return self.family.computation_strategy + + def log_likelihood(self, batch: Sample) -> float: + raise NotImplementedError + + def sample(self, n: int, **options: Any) -> Sample: + """ + Generate samples from this distribution. + + Parameters + ---------- + n : int + Number of samples to generate. + **options : Any + Additional options for the sampling algorithm. + + Returns + ------- + Sample + The generated samples. + """ + return self.sampling_strategy.sample(n, distr=self, **options) diff --git a/src/pysatl_core/families/parametric_family.py b/src/pysatl_core/families/parametric_family.py new file mode 100644 index 0000000..e655efe --- /dev/null +++ b/src/pysatl_core/families/parametric_family.py @@ -0,0 +1,157 @@ +""" +Parametric family definitions and management infrastructure. + +This module contains the main class for defining parametric families of +distributions, including support for multiple parameterizations, distribution +characteristics, sampling strategies, and computation methods. It serves as +the central definition point for statistical distribution families. +""" + +from __future__ import annotations + +__author__ = "Leonid Elkin, Mikhail, Mikhailov" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + +from collections.abc import Callable +from typing import Any + +from pysatl_core.distributions import ( + ComputationStrategy, + SamplingStrategy, +) +from pysatl_core.families.distribution import ParametricFamilyDistribution +from pysatl_core.families.parametrizations import Parametrization, ParametrizationSpec +from pysatl_core.types import ( + DistributionType, + GenericCharacteristicName, + ParametrizationName, +) + +type ParametrizedFunction = Callable[[Parametrization, Any], Any] + + +class ParametricFamily: + """ + A family of distributions with multiple parametrizations. + + This class represents a parametric family of distributions, such as + the normal or lognormal family, which can be parameterized in different + ways (e.g., mean-variance or canonical parametrization). + + Attributes + ---------- + name : str + Name of the distribution family. + distr_type : DistributionType | Callable[[Parametrization] | DistributionType] + Type of distributions in this family. + parametrizations : ParametrizationSpec + Specification of available parametrizations. + distr_parametrizations : + + distr_characteristics : Dict[GenericCharacteristicName, Callable[[Any, Any], Any]] + Mapping from characteristic names to computation functions. + sampling_strategy : SamplingStrategy + Strategy for sampling from distributions in this family. + computation_strategy : ComputationStrategy + Strategy for computing distribution characteristics. + """ + + def __init__( + self, + name: str, + distr_type: DistributionType | Callable[[Parametrization], DistributionType], + distr_parametrizations: list[ParametrizationName], + distr_characteristics: dict[ + GenericCharacteristicName, + dict[ParametrizationName, ParametrizedFunction] | ParametrizedFunction, + ], + sampling_strategy: SamplingStrategy, + computation_strategy: ComputationStrategy[Any, Any], + ): + """ + Initialize a new parametric family. + + Parameters + ---------- + name : str + Name of the distribution family. + + distr_type : DistributionType | Callable[[Parametrization], DistributionType] + Type of distributions in this family or, if type is parameter-depended, function + that takes as input *base* parametrization and inferes type based on it. + + distr_parametrizations : List[ParametrizationName] + List of parametrizations for this distribution. *First parametrization is always + base parametrization*. + + distr_characteristics: + Mapping from characteristics names to computation functions or dictionary of those, + if for multiple parametrizations same characteristic available. + + sampling_strategy : SamplingStrategy + Strategy for sampling from distributions in this family. + + computation_strategy : ComputationStrategy + Strategy for computing distribution characteristics. + """ + self._name = name + self._distr_type: Callable[[Parametrization], DistributionType] = ( + (lambda params: distr_type) if isinstance(distr_type, DistributionType) else distr_type + ) + + # Parametrizations must be built by user + self.parametrization_names = distr_parametrizations + self.parametrizations = ParametrizationSpec(self.parametrization_names[0]) + + self.sampling_strategy = sampling_strategy + self.computation_strategy = computation_strategy + + def _process_char_val( + value: dict[ParametrizationName, ParametrizedFunction] | ParametrizedFunction, + ) -> dict[ParametrizationName, ParametrizedFunction]: + return value if isinstance(value, dict) else {self.parametrization_names[0]: value} + + self.distr_characteristics = { + key: _process_char_val(value) for key, value in distr_characteristics.items() + } + + @property + def name(self) -> str: + return self._name + + def distribution( + self, parametrization_name: str | None = None, **parameters_values: Any + ) -> ParametricFamilyDistribution: + """ + Create a distribution instance with the given parameters. + + Parameters + ---------- + parametrization_name : str | None, optional + Name of the parametrization to use, or None for base parametrization. + **parameters_values + Parameter values for the distribution. + + Returns + ------- + ParametricFamilyDistribution + A distribution instance with the specified parameters. + + Raises + ------ + ValueError + If the parameters don't satisfy the parametrization constraints. + """ + if parametrization_name is None: + parametrization_class = self.parametrizations.base + else: + parametrization_class = self.parametrizations.parametrizations[parametrization_name] + + parameters = parametrization_class(**parameters_values) + base_parameters = self.parametrizations.get_base_parameters(parameters) + parameters.validate() + distribution_type = self._distr_type(base_parameters) + return ParametricFamilyDistribution(self.name, distribution_type, parameters) + + __call__ = distribution diff --git a/src/pysatl_core/families/parametrizations.py b/src/pysatl_core/families/parametrizations.py new file mode 100644 index 0000000..2f6a7af --- /dev/null +++ b/src/pysatl_core/families/parametrizations.py @@ -0,0 +1,314 @@ +""" +Parameterization classes and specifications for distribution families. + +This module provides the core classes for defining different parameterizations +of statistical distributions, including constraints validation and conversion +between parameterization formats. +""" + +from __future__ import annotations + +__author__ = "Leonid Elkin, Mikhail, Mikhailov" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + +from abc import ABC, abstractmethod +from collections.abc import Callable +from dataclasses import dataclass +from functools import wraps +from typing import TYPE_CHECKING, Any, ClassVar, Protocol, runtime_checkable + +from pysatl_core.types import ParametrizationName + +if TYPE_CHECKING: + from pysatl_core.families.parametric_family import ParametricFamily + + +@runtime_checkable +class ParametrizationConstraintProtocol(Protocol): + @property + def _is_constraint(self) -> bool: ... + @property + def _constraint_description(self) -> str: ... + + def __call__(self, **kwargs: Any) -> bool: ... + + +@dataclass(slots=True, frozen=True) +class ParametrizationConstraint: + """ + A constraint on parameter values for a parametrization. + + Attributes + ---------- + description : str + Human-readable description of the constraint. + check : Callable[[Any], bool] + Function that validates the constraint, returning True if satisfied. + """ + + description: str + check: Callable[[Any], bool] + + +class Parametrization(ABC): + """ + Abstract base class for distribution parametrizations. + + This class defines the interface that all parametrizations must implement, + including parameter validation and conversion to base parametrization. + + Attributes + ---------- + constraints : ClassVar[List[ParametrizationConstraint]] + Class-level list of constraints that apply to this parametrization. + """ + + _constraints: ClassVar[list[ParametrizationConstraint]] = [] + + @property + @abstractmethod + def name(self) -> str: + """ + Get the name of this parametrization. + + Returns + ------- + str + The name of the parametrization. + """ + + @property + @abstractmethod + def parameters(self) -> dict[str, Any]: + """ + Get the parameters as a dictionary. + + Returns + ------- + Dict[str, Any] + Dictionary mapping parameter names to values. + """ + + @property + def constraints(self) -> list[ParametrizationConstraint]: + """ + Get the constraints for this parametrization. + + Returns + ------- + List[ParametrizationConstraint] + List of constraints that apply to this parametrization. + """ + return self._constraints + + def validate(self) -> None: + """ + Validate all constraints for this parametrization. + + Raises + ------ + ValueError + If any constraint is not satisfied. + """ + for constraint in self._constraints: + if not constraint.check(self): + raise ValueError(f'Constraint "{constraint.description}" does not hold') + + def transform_to_base_parametrization(self) -> Parametrization: + """ + Convert this parametrization to the base parametrization. + + Returns + ------- + Parametrization + The equivalent parameters in the base parametrization. + + Notes + ----- + The base implementation returns self, assuming this is already + the base parametrization. Subclasses should override this method + if they need to convert to a different parametrization. + """ + return self + + +class ParametrizationSpec: + """ + Container for all parametrizations of a distribution family. + + This class manages the collection of parametrizations for a family + and handles conversions between them. + + Attributes + ---------- + parametrizations : Dict[ParametrizationName, Type[Parametrization]] + Mapping from parametrization names to parametrization classes. + base_parametrization_name : ParametrizationName | None + Name of the base parametrization, if defined. + """ + + def __init__(self, base_name: ParametrizationName) -> None: + """Initialize an empty parametrization specification.""" + self.parametrizations: dict[ParametrizationName, type[Parametrization]] = {} + self.base_parametrization_name: ParametrizationName = base_name + + @property + def base(self) -> type[Parametrization]: + """ + Get the base parametrization class. + + Returns + ------- + Type[Parametrization] + The base parametrization class. + + Raises + ------ + ValueError + If no base parametrization has been defined. + """ + if self.base_parametrization_name is None: + raise ValueError("No base parametrization defined") + return self.parametrizations[self.base_parametrization_name] + + def add_parametrization( + self, + name: ParametrizationName, + parametrization_class: type[Parametrization], + ) -> None: + """ + Add a new parametrization to the specification. + + Parameters + ---------- + name : ParametrizationName + Name of the parametrization. + parametrization_class : Type[Parametrization] + Class implementing the parametrization. + """ + self.parametrizations[name] = parametrization_class + + def get_base_parameters(self, parameters: Parametrization) -> Parametrization: + """ + Convert parameters to the base parametrization. + + Parameters + ---------- + parameters : Parametrization + Parameters in any parametrization. + + Returns + ------- + Parametrization + Equivalent parameters in the base parametrization. + """ + if parameters.name == self.base_parametrization_name: + return parameters + else: + return parameters.transform_to_base_parametrization() + + +# Decorators for declarative syntax +def constraint( + description: str, +) -> Callable[[Callable[[Any], bool]], ParametrizationConstraintProtocol]: + """ + Decorator to mark a method as a parameter constraint. + + Parameters + ---------- + description : str + Human-readable description of the constraint. + + Returns + ------- + Callable + Decorator function that marks the method as a constraint. + + Examples + -------- + >>> @constraint("sigma > 0") + >>> def check_sigma_positive(self): + >>> return self.sigma > 0 + """ + + def decorator(func: Callable[[Any], bool]) -> ParametrizationConstraintProtocol: + @wraps(func) + def wrapper(*args, **kwargs): # type: ignore + return func(*args, **kwargs) + + wrapper._is_constraint = True # type: ignore + wrapper._constraint_description = description # type: ignore + return wrapper # type: ignore + + return decorator + + +def parametrization( + family: ParametricFamily, name: str +) -> Callable[[type[Parametrization]], type[Parametrization]]: + """ + Decorator to register a class as a parametrization for a family. + + Parameters + ---------- + family : ParametricFamily + The family to register the parametrization with. + name : str + Name of the parametrization. + base : bool, optional + Whether this is the base parametrization, by default False. + + Returns + ------- + Callable + Decorator function that registers the class as a parametrization. + + Examples + -------- + >>> @parametrization(family=NormalFamily, name='meanvar') + >>> class MeanVarParametrization: + >>> mean: float + >>> var: float + """ + + def decorator(cls: type[Parametrization]) -> type[Parametrization]: + # Convert to dataclass if not already + if not hasattr(cls, "__dataclass_fields__"): + cls = dataclass(cls) + + # Add name property + def name_property(self): # type: ignore + return name + + cls.name = property(name_property) # type: ignore + + # Add parameters property + def parameters_property(self): # type: ignore + return { + field.name: getattr(self, field.name) + for field in self.__dataclass_fields__.values() + } + + cls.parameters = property(parameters_property) # type: ignore + + # Collect constraints + constraints = [] + for attr_name in dir(cls): + attr = getattr(cls, attr_name) + if hasattr(attr, "_is_constraint") and attr._is_constraint: + constraints.append( + ParametrizationConstraint(description=attr._constraint_description, check=attr) + ) + cls._constraints = constraints + + # Add validate method + cls.validate = Parametrization.validate # type: ignore + + # Register with family + family.parametrizations.add_parametrization(name, cls) + + return cls + + return decorator diff --git a/src/pysatl_core/families/registry.py b/src/pysatl_core/families/registry.py new file mode 100644 index 0000000..ba12dc6 --- /dev/null +++ b/src/pysatl_core/families/registry.py @@ -0,0 +1,99 @@ +""" +Global registry for parametric distribution families using singleton pattern. + +This module implements a centralized registry that maintains references to all +defined parametric families, enabling easy access and management across the +application. The registry follows the singleton pattern to ensure consistency. +""" + +from __future__ import annotations + +__author__ = "Leonid Elkin, Mikhail, Mikhailov" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + +from typing import TYPE_CHECKING, ClassVar + +if TYPE_CHECKING: + from .parametric_family import ParametricFamily + + +class ParametricFamilyRegister: + """ + Singleton registry for parametric distribution families. + + This class maintains a global registry of all parametric families, + allowing them to be accessed by name from anywhere in the codebase. + + Examples + -------- + >>> registry = ParametricFamilyRegister() + >>> family = registry.get('Lognormal Family') + """ + + _instance: ClassVar[ParametricFamilyRegister | None] = None + _registered_families: dict[str, ParametricFamily] + + def __new__(cls) -> ParametricFamilyRegister: + """ + Create or return the singleton instance. + + Returns + ------- + ParametricFamiliesRegister + The singleton registry instance. + """ + if cls._instance is None: + cls._instance = super().__new__(cls) + cls._instance._registered_families = {} + return cls._instance + + @classmethod + def get(cls, name: str) -> ParametricFamily: + """ + Retrieve a parametric family by name. + + Parameters + ---------- + name : str + The name of the family to retrieve. + + Returns + ------- + ParametricFamily + The requested parametric family. + + Raises + ------ + ValueError + If no family with the given name exists in the registry. + """ + self = cls() + if name not in self._registered_families: + raise ValueError(f"No family {name} found in register") + return self._registered_families[name] + + @classmethod + def register(cls, family: ParametricFamily) -> None: + """ + Register a new parametric family. + + Parameters + ---------- + family : ParametricFamily + The family to register. + + Raises + ------ + ValueError + If family with the same name already was registered + """ + self = cls() + if family.name in self._registered_families: + raise ValueError(f"Family {family.name} already found in register") + self._registered_families[family.name] = family + + +def _reset_families_register_for_tests() -> None: + """Reset the cached distribution type register (test helper).""" + ParametricFamilyRegister._instance = None diff --git a/src/pysatl_core/types.py b/src/pysatl_core/types.py index 449151f..a533c8e 100644 --- a/src/pysatl_core/types.py +++ b/src/pysatl_core/types.py @@ -44,14 +44,21 @@ class EuclideanDistributionType(DistributionType): dimension: int +UnivariateContinuous = EuclideanDistributionType(kind=Kind.CONTINUOUS, dimension=1) +UnivariateDiscrete = EuclideanDistributionType(kind=Kind.DISCRETE, dimension=1) + type GenericCharacteristicName = str +type ParametrizationName = str ScalarFunc = Callable[[float], float] __all__ = [ "Kind", "EuclideanDistributionType", + "UnivariateContinuous", + "UnivariateDiscrete", "GenericCharacteristicName", + "ParametrizationName", "DistributionType", "ScalarFunc", ] diff --git a/tests/conftest.py b/tests/conftest.py index 9ef5190..f8c3e4c 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -6,12 +6,13 @@ import pytest from pysatl_core.distributions.registry import _reset_distribution_type_register_for_tests +from pysatl_core.families.registry import _reset_families_register_for_tests pytest.importorskip("scipy") @pytest.fixture(autouse=True) -def _fresh_registry() -> Generator[None, Any, None]: +def _fresh_registries() -> Generator[None, Any, None]: _reset_distribution_type_register_for_tests() + _reset_families_register_for_tests() yield - _reset_distribution_type_register_for_tests() diff --git a/tests/test_family_and_distribution.py b/tests/test_family_and_distribution.py new file mode 100644 index 0000000..9f45f69 --- /dev/null +++ b/tests/test_family_and_distribution.py @@ -0,0 +1,112 @@ +from pysatl_core.families import ( + ParametricFamily, + ParametricFamilyDistribution, + ParametricFamilyRegister, + Parametrization, +) + +PDF = "pdf" + + +class TestParametricFamily: + """Test the ParametricFamily class.""" + + def test_family_creation(self): + """Test creating a parametric family.""" + + # Mock strategies + class MockSamplingStrategy: + pass + + class MockComputationStrategy: + pass + + # Mock characteristic function + def mock_pdf(parameters, x): + return x * 2 + + # Create family + family = ParametricFamily( + name="TestFamily", + distr_type="Continuous", + distr_parametrizations=["mock"], + distr_characteristics={PDF: mock_pdf}, + sampling_strategy=MockSamplingStrategy(), + computation_strategy=MockComputationStrategy(), + ) + register = ParametricFamilyRegister() + register.register(family) + + # Check properties + assert family.name == "TestFamily" + assert family._distr_type == "Continuous" + assert PDF in family.distr_characteristics + assert isinstance(family.sampling_strategy, MockSamplingStrategy) + assert isinstance(family.computation_strategy, MockComputationStrategy) + + +class TestParametricFamilyDistribution: + """Test the ParametricFamilyDistribution class.""" + + def test_distribution_creation(self): + """Test creating a distribution instance.""" + + # Create a mock family + class MockSamplingStrategy: + def sample(self, n, distr, **options): + return [1, 2, 3] # Mock samples + + class MockComputationStrategy: + pass + + def mock_pdf(parameters, x): + return x * parameters.value + + family = ParametricFamily( + name="MockFamily", + distr_type="Continuous", + distr_parametrizations=["mock"], + distr_characteristics={PDF: mock_pdf}, + sampling_strategy=MockSamplingStrategy(), + computation_strategy=MockComputationStrategy(), + ) + + # Create a mock parametrization + class MockParametrization(Parametrization): + def __init__(self, value): + self.value = value + + @property + def name(self): + return "mock" + + @property + def parameters(self): + return {"value": self.value} + + # Add to family + family.parametrizations.add_parametrization("mock", MockParametrization) + + # Register family + register = ParametricFamilyRegister() + register.register(family) + + # Create distribution + params = MockParametrization(2.0) + + dist = ParametricFamilyDistribution("MockFamily", "Continuous", params) + + # Check properties + assert dist.distr_name == "MockFamily" + assert dist.distribution_type == "Continuous" + assert dist.parameters is params + assert dist.family is family + + # Test sampling + samples = dist.sample(3) + assert samples == [1, 2, 3] + + # Test analytical computations + computations = dist.analytical_computations + assert PDF in computations + assert computations[PDF].func(5.0) == 10.0 # 5.0 * 2.0 diff --git a/tests/test_parameters.py b/tests/test_parameters.py new file mode 100644 index 0000000..8a79b69 --- /dev/null +++ b/tests/test_parameters.py @@ -0,0 +1,187 @@ +import pytest + +from pysatl_core.families import ( + ParametricFamily, + Parametrization, + ParametrizationConstraint, + ParametrizationSpec, + constraint, + parametrization, +) + + +class TestParametrizationConstraint: + """Test the ParametrizationConstraint class.""" + + def test_constraint_creation(self): + """Test creating a constraint with description and check function.""" + + def check_func(obj): + return obj.value > 0 + + constraint = ParametrizationConstraint( + description="Value must be positive", check=check_func + ) + + assert constraint.description == "Value must be positive" + assert constraint.check is check_func + + +class TestParametrization: + """Test the base Parametrization class.""" + + def test_abstract_methods(self): + """Test that Parametrization is abstract and requires name and parameters properties.""" + with pytest.raises(TypeError): + Parametrization() # Can't instantiate abstract class + + # Create a concrete implementation + class ConcreteParametrization(Parametrization): + @property + def name(self): + return "concrete" + + @property + def parameters(self): + return {"param": 1.0} + + # Should be able to instantiate + param = ConcreteParametrization() + assert param.name == "concrete" + assert param.parameters == {"param": 1.0} + + +class TestParametrizationSpec: + """Test the ParametrizationSpec class.""" + + def test_add_and_get_parametrization(self): + """Test adding and retrieving parametrizations.""" + spec = ParametrizationSpec(base_name="mock") + + # Create a mock parametrization class + class MockParametrization(Parametrization): + @property + def name(self): + return "mock" + + @property + def parameters(self): + return {} + + # Add parametrization + spec.add_parametrization("mock", MockParametrization) + + # Check it was added + assert "mock" in spec.parametrizations + assert spec.parametrizations["mock"] is MockParametrization + assert spec.base_parametrization_name == "mock" + assert spec.base is MockParametrization + + def test_get_base_parameters(self): + """Test converting parameters to base parametrization.""" + spec = ParametrizationSpec(base_name="base") + + # Create mock parametrizations + class BaseParametrization(Parametrization): + def __init__(self, value): + self.value = value + + @property + def name(self): + return "base" + + @property + def parameters(self): + return {"value": self.value} + + class OtherParametrization(Parametrization): + def __init__(self, other_value): + self.other_value = other_value + + @property + def name(self): + return "other" + + @property + def parameters(self): + return {"other_value": self.other_value} + + def transform_to_base_parametrization(self): + return BaseParametrization(self.other_value * 2) + + # Add parametrizations + spec.add_parametrization("base", BaseParametrization) + spec.add_parametrization("other", OtherParametrization) + + # Test with base parametrization + base_params = BaseParametrization(5.0) + result = spec.get_base_parameters(base_params) + assert result is base_params + + # Test with other parametrization + other_params = OtherParametrization(3.0) + result = spec.get_base_parameters(other_params) + assert isinstance(result, BaseParametrization) + assert result.value == 6.0 # 3.0 * 2 + + +class TestDecorators: + """Test the constraint and parametrization decorators.""" + + def test_constraint_decorator(self): + """Test the constraint decorator.""" + + @constraint("Value must be positive") + def check_positive(self): + return self.value > 0 + + # Check that the decorator added the required attributes + assert hasattr(check_positive, "_is_constraint") + assert hasattr(check_positive, "_constraint_description") + assert check_positive._is_constraint is True + assert check_positive._constraint_description == "Value must be positive" + + def test_parametrization_decorator(self): + """Test the parametrization decorator.""" + # Create a mock family + family = ParametricFamily( + name="TestFamily", + distr_type="Continuous", + distr_parametrizations=["test"], + distr_characteristics={}, + sampling_strategy=None, + computation_strategy=None, + ) + + # Apply the decorator + @parametrization(family=family, name="test") + class TestParametrization: + value: float + + @constraint("Value must be positive") + def check_positive(self): + return self.value > 0 + + # Check that the class was modified + assert hasattr(TestParametrization, "name") + assert hasattr(TestParametrization, "parameters") + assert hasattr(TestParametrization, "_constraints") + assert hasattr(TestParametrization, "validate") + + # Check that it was registered with the family + assert "test" in family.parametrizations.parametrizations + assert family.parametrizations.parametrizations["test"] is TestParametrization + assert family.parametrizations.base_parametrization_name == "test" + + # Test instantiation and validation + instance = TestParametrization(value=5.0) + assert instance.name == "test" + assert instance.parameters == {"value": 5.0} + + # This should not raise an exception + instance.validate() + + # Test with invalid parameters + invalid_instance = TestParametrization(value=-1.0) + with pytest.raises(ValueError, match="Constraint.*does not hold"): + invalid_instance.validate() diff --git a/tests/test_parametric_families_registry.py b/tests/test_parametric_families_registry.py new file mode 100644 index 0000000..3c931e5 --- /dev/null +++ b/tests/test_parametric_families_registry.py @@ -0,0 +1,33 @@ +import pytest + +from pysatl_core.families import ParametricFamilyRegister + + +class TestParametricFamiliesRegister: + """Test the ParametricFamiliesRegister singleton.""" + + def test_singleton_pattern(self): + """Test that only one instance exists.""" + register1 = ParametricFamilyRegister() + register2 = ParametricFamilyRegister() + assert register1 is register2 + + def test_register_and_get_family(self): + """Test registering and retrieving a family.""" + register = ParametricFamilyRegister() + + # Create a mock family + mock_family = type("MockFamily", (), {"name": "TestFamily"})() + + # Register and retrieve + register.register(mock_family) + retrieved = register.get("TestFamily") + + assert retrieved is mock_family + + def test_get_nonexistent_family(self): + """Test error when getting a non-existent family.""" + register = ParametricFamilyRegister() + + with pytest.raises(ValueError, match="No family Nonexistent found in register"): + register.get("Nonexistent")