{ "cells": [ { "cell_type": "markdown", "id": "f1ff00", "metadata": { "collapsed": false }, "source": [ "##### **Supplementary code for the paper \"Stochastic effects on plankton dynamics: insights from a realistic 0\\-dimensional marine biogeochemical model\" by G. Occhipinti et al.**\n", "\n", "Below we report the figure corresponding to the analytic solutions $q(x)$ and $p(x)$, see main text, for the parameters $\\alpha=\\frac{2\\gamma}{D^2}=2$ and $x_0=7.6$, the latter indicated with a dashed vertical line. The plot and calculation has been performed using this SageMath python package \\(see in the theory section of the main text\\). The average of the numerically calculated distribution p\\(x\\) corresponds to $x_0$ for any value of $\\alpha$.\n", "\n", "![](.AnalyticSolution.ipynb.upload/Fig4_singlecolumn.png)\n", "\n", "In the code below it is possible to reproduce the same plot, or a different one by manipulating the control parameters. When the notebook is launched the plot will be displayed after cell 6.\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "4475c3", "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "%display latex\n", "\n", "# Define some constants that we will use later\n", "c = var('c')\n", "C = var('C')\n", "gamma = var('gamma')\n", "D = var('D')\n", "x0 = var('x0')\n", "alpha = var('alpha')\n", "\n", "assume(D > 0)\n", "assume(gamma, 'real')\n", "assume(x0, 'real')\n", "assume(x0 > 0)\n", "assume(c, 'real')\n", "assume(C, 'real')\n", "assume(alpha < 0)\n", "\n", "p = function('p')(x)\n", "q = function('q')(x)" ] }, { "cell_type": "code", "execution_count": 3, "id": "7318c7", "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle \\frac{\\partial}{\\partial x}q\\left(x\\right) = c - \\frac{2 \\, \\gamma {\\left(x - x_{0}\\right)} q\\left(x\\right)}{D^{2} x^{2}}\\)" ], "text/latex": [ "$\\displaystyle \\frac{\\partial}{\\partial x}q\\left(x\\right) = c - \\frac{2 \\, \\gamma {\\left(x - x_{0}\\right)} q\\left(x\\right)}{D^{2} x^{2}}$" ], "text/plain": [ "diff(q(x), x) == c - 2*gamma*(x - x0)*q(x)/(D^2*x^2)" ] }, "execution_count": 3, "metadata": { }, "output_type": "execute_result" } ], "source": [ "equ_diff_1 = derivative(q,x,1) == - 2/(D*D)*gamma*(x-x0)/(x*x)*q + c\n", "\n", "equ_diff_1" ] }, { "cell_type": "code", "execution_count": 4, "id": "d1b140", "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "CONSTANT_VALUES = {\n", " D: 1.0e-5,\n", " gamma: 3.5 * 1.0e-5,\n", " x0 : 7.6,\n", " c: 0.,\n", " C: 1.0\n", "}" ] }, { "cell_type": "code", "execution_count": 5, "id": "3edead", "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle \\frac{C e^{\\left(-\\frac{2 \\, \\gamma x_{0}}{D^{2} x}\\right)}}{x^{\\frac{2 \\, \\gamma}{D^{2}}}}\\)" ], "text/latex": [ "$\\displaystyle \\frac{C e^{\\left(-\\frac{2 \\, \\gamma x_{0}}{D^{2} x}\\right)}}{x^{\\frac{2 \\, \\gamma}{D^{2}}}}$" ], "text/plain": [ "C*e^(-2*gamma*x0/(D^2*x))/x^(2*gamma/D^2)" ] }, "execution_count": 5, "metadata": { }, "output_type": "execute_result" } ], "source": [ "# q analytical solution\n", "\n", "a = -2/(D*D)*gamma\n", "q = C * x^a *exp(a * x0 / x)\n", "\n", "q" ] }, { "cell_type": "code", "execution_count": 6, "id": "e81ffb", "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle 2.22880036366235 \\times 10^{-920576}\\)" ], "text/latex": [ "$\\displaystyle 2.22880036366235 \\times 10^{-920576}$" ], "text/plain": [ "2.22880036366235e-920576" ] }, "execution_count": 6, "metadata": { }, "output_type": "execute_result" } ], "source": [ "q_num_P1=q.substitute(CONSTANT_VALUES)\n", "q_num_P1.substitute(x=7.6)" ] }, { "cell_type": "code", "execution_count": 9, "id": "b6047a", "metadata": { "collapsed": false }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "9891526dc3d749599a41ca30c36502fa", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Interactive function with 1 widget\n", " a_value: SelectionSlider(descripti…" ] }, "execution_count": 9, "metadata": { }, "output_type": "execute_result" } ], "source": [ "from functools import lru_cache\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "\n", "a_VALUES = sorted(\n", " [\"-1 * 10^{}\".format(i) for i in range(0, 7)] +\n", " [\"-5 * 10^{}\".format(i) for i in range(0, 7)] +\n", " [\"-2 * 10^{}\".format(i) for i in range(0, 7)], key=lambda t: sage_eval(t))\n", "\n", "RRR = RealField(Integer(512))\n", "\n", "COLOR1 = \"tab:blue\"\n", "COLOR2 = \"tab:red\"\n", "\n", "@interact\n", "def draw_plot(a_value = slider(a_VALUES, default = \"-1 * 10^0\", label=\"a\")):\n", " x0_val = CONSTANT_VALUES[x0]\n", " a_value = sage_eval(a_value)\n", " q_espression= x^a_value * exp(a_value * x0_val / x)\n", " q_espression_n = q_espression/q_espression.substitute(x=x0_val)\n", "\n", " q_espression_n_compiled = fast_callable(q_espression_n, vars=[x], domain=RRR)\n", "\n", " r = q_espression_n / x^2\n", " r_espression_n_compiled = fast_callable(r, vars=[x], domain=RRR)\n", " r_espression_norm_int1, _ = numerical_integral(r_espression_n_compiled, a=0, b=10, algorithm=\"qag\", rule=6)\n", " r_espression_norm_int2, _ = numerical_integral(r_espression_n_compiled, a=10, b=+oo, algorithm=\"qag\", rule=6)\n", " r_espression_norm = r_espression_norm_int1 + r_espression_norm_int2\n", "\n", " r_final = lambda t: r_espression_n_compiled(t) / r_espression_norm\n", "\n", " p1_plot = plot(\n", " q_espression_n_compiled,\n", " xmin=0.01,\n", " xmax=20,\n", " color=\"blue\"\n", " )\n", " p2_plot = plot(\n", " r_final,\n", " xmin=0.01,\n", " xmax=20,\n", " color=\"red\"\n", " )\n", "\n", " p1_x_points = np.array(p1_plot._objects[0].xdata, dtype=np.float64)\n", " p1_y_points = np.array(p1_plot._objects[0].ydata, dtype=np.float64)\n", "\n", " p2_x_points = np.array(p2_plot._objects[0].xdata, dtype=np.float64)\n", " p2_y_points = np.array(p2_plot._objects[0].ydata, dtype=np.float64)\n", "\n", " plt.xlabel(\"$x$\", fontsize=14)\n", " plt.title(\"Probabilty distribution at equilibrium\")\n", "\n", " plt.gca().spines['top'].set_visible(False)\n", " plt.gca().spines['bottom'].set_position('zero')\n", "\n", " plt.ylabel(\"$p$\", fontsize=14, color=COLOR1)\n", " plt.plot(p2_x_points, p2_y_points, color=COLOR1, linewidth=2.5)\n", " plt.fill_between(p2_x_points, p2_y_points, color=COLOR1, alpha=0.1)\n", "\n", " plt.gca().tick_params(axis='y', colors=COLOR1)\n", "\n", " ax2=plt.gca().twinx()\n", "\n", " plt.gca().spines['top'].set_visible(False)\n", " plt.gca().spines['bottom'].set_visible(False)\n", " plt.gca().spines['right'].set_color(COLOR2)\n", " plt.gca().spines['left'].set_color(COLOR1)\n", "\n", " plt.gca().tick_params(axis='y', colors=COLOR2)\n", "\n", " plt.ylabel(\"$q$\", fontsize=14, color=COLOR2)\n", " plt.plot(p1_x_points, p1_y_points, color=COLOR2, linewidth=1.5)\n", " plt.gca().set_xlim(left=0)\n", "\n", " v_line_bottom = (x0_val, 0)\n", " v_line_bottom_display = plt.gca().transData.transform(v_line_bottom)\n", " v_line_bottom_axis = plt.gca().transAxes.inverted().transform(v_line_bottom_display)\n", " ymin_line = v_line_bottom_axis[1]\n", "\n", " plt.axvline(x = x0_val, ymin=ymin_line, color = 'tab:gray', label = '$x_0$', linestyle='--', linewidth=1)" ] }, { "cell_type": "markdown", "id": "faff5b", "metadata": { "collapsed": false }, "source": [ "### Mode and Median ###\n", "\n", "In the following cell, we compute the deviations depending on the parameter $\\alpha$ of the mode and median with respect to the average value that always corresponds to $x_0$.\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "99da3e", "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHWCAYAAAARl3+JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABu8ElEQVR4nO3deZyN5f/H8deZYezGPlnGkiKDjBlj34lStmTfQynZl5JKKSlL2gZR1iiSpaQk2SJhGGWLRJSxL8Nghpnz++P6mh9ZmplzZu5zn3k/H4/z+JozZ879Pu6vfFzL53I4nU4nIiIiImJrPlYHEBERERHXqagTERER8QIq6kRERES8gIo6ERERES+gok5ERETEC6ioExEREfECKupEREREvICKOhEREREvkMHqAN4mISGBo0ePkiNHDhwOh9VxRERExMM5nU4uXLhAoUKF8PFJ+Xibijo3O3r0KIGBgVbHEBEREZs5cuQIRYoUSfHPq6hzsxw5cgDmxuTMmdPiNCIiIuLpoqOjCQwMTKwhUkpFnZtdn3LNmTOnKep27oSvv4YGDSA0FHx9LU4oIiIinsjVZVvaKJHadu6EN9+EKlUgb15o0QI++AB27wan0+p0IiIi4iVU1LlJeHg4QUFBhIWF3fyNdu3gzBnYsAEGD4Zz52DIEChbFgoXhs6d4bPPzGtEREREUsjhdGq4yJ2io6Px9/fn/Pnzd15Td+mSKfJWrYIVKyAyEnx8oHp1aNIEHn0UypcH7Z4VERHxekmqHZJARZ2bpejG/PMPfPstfPMNrFwJMTFQpMj/F3gNGkC2bKkbXEREJBni4uK4du2a1TFsJUOGDPj5+d3yvIo6D+XyjYmNhXXrYPlyU+Tt3w9ZspgC74knTJHn4u4YERERV8TFxbFr1y4SEhKsjmIrPj4+lC1b9pbCzl1FnXa/eppMmeChh8xj4kTYtw8WL4aFC6F9e/P9hx82BV7TpuDvb3ViERFJZ65du0ZCQgLFixcnS5YsVsexhcuXL3Po0CGuXbt229E6d1BR5+lKlYLnnzePQ4dg0SJT4HXuDBkzmuLviSegeXPIk8fqtCIiko5kyZKFrFmzWh3DVs6fP4/T6SRbKiyr0u5XOyleHAYNgo0b4cgRGD8eLlyAHj2gYEF4/HEzqhcba3VSERERuY1ly5Yxf/58YmJi3P7eKurc5I4tTVJLkSLQr59Zf/fPP/D22/DXX6awK1gQevc2O2y1ZFJERMRjZMyYkejoaGJTYQBGRZ2b9OnTh927d7Nly5a0v3jBgjBgAEREwK5dpqBbvhxq1oSSJeGVV+DAgbTPJSIiIjdJrfV0oKLO+wQFmRMsDh2C1auhXj1491247z7TGuWzz+DKFatTioiIiJupqPNWPj5Qty588gkcOwazZ8PVq9ChgznJYsAAM6onIiIiXkFFXXqQNavZLbtuHezdC08+CfPmQblyUK2aKfg0eiciIl6qbt26DBgwIPHr4sWL8+6771qWJ7WoqEtvSpeGcePg77/hiy9MI+OuXaFoURgxwuyqFRER8WJbtmzhqaeesjqG26moS6/8/Ex/u++/N6N37dvDBx9AiRLQqhWsWaOdsyIi4pXy58/vlf31VNSJGb177z3TGuX992HPHrPBonx5mDIFLl2yOqGIiHihunXr0rdvXwYMGEDu3LkJCAhg6tSpxMTE0L17d3LkyEHJkiX59ttvE39m9+7dNGnShOzZsxMQEEDnzp05depU4vdjYmLo0qUL2bNnp2DBgkyYMOGW6/57+vWdd96hfPnyZMuWjcDAQJ599lkuXryY+P2ZM2eSK1cuVqxYQZkyZciePTsPP/wwUVFRqfMbk0Iq6lLZzz+bwx5+/tnqJEmQIwc8+6zZQLFqlTnNok8fCAyEl1+G48etTigiIl5m1qxZ5MuXj82bN9O3b1+eeeYZWrduTfXq1dm2bRuNGzemc+fOXLp0iaioKOrUqUNwcDBbt27lu+++4/jx47Rp0ybx/YYOHcrq1atZvHgx33//PWvWrCEiIuKuGXx8fHj//ffZuXMns2bN4scff2TYsGE3vebSpUuMHz+eOXPmsG7dOg4fPsyQIUNS5fckxZziVufPn3cCzvPnzzudTqfzhx+czgcecDrB6axVy+n85hunMyHB4pDJ8eefTmf//k5ntmxOp5+f09mjh9O5a5fVqURExEIxMTHOrVu3OmNiYm54zumMiPjvxw0/4qxTp46zZs2aiV9fu3bNmS1bNmfnzp0Tn4uKinICzp9//tn58ssvOxs1anRTliNHjjgB5++//+68cOGC08/Pz/n5558nfv/06dPOLFmyOPv375/4XLFixZwTJ0684+dbsGCBM2/evIlfz5gxwwk4//jjj8TnwsPDnQEBAUn6/TK/P+b3bO7cuc6JEyc6T58+nfi9f9cOKaWROje504kSDRqYga/FiyEuDh59FCpUgE8/NR1GPF6JEqbP3ZEjMGqUaWpctqz5IGvXat2diIgAZnl2aOh/P/buvfnnHnzwwcRf+/r6kjdvXsqXL5/4XEBAAAAnTpwgIiKC1atXkz179sTHAw88AMCBAwc4cOAAcXFxVKtWLfHn8+TJQ+nSpe+affXq1Tz00EMULlyYHDly0KVLF06fPn3TUV5Zs2alZMmSiV8XLFiQEydOJP83KhWpqHOTu50o4eMDLVqYKdg1a8wJX507w/33m70Jtliyljs3PP+8aWo8cyYcPmz64NWuDd99p+JORCSde+ABc7DRfz3+V4Mlypgx401fOxyOm55zOBwAJCQkkJCQQNOmTYmMjLzpsX//fmrXro0zBX8X/fXXXzRp0oRy5crx5ZdfEhERQXh4OABXbxh9uV3OlFwvNamoS0MOB9SpYwa7IiOhRg3TA7hYMXj9dTh/3uqESeDnZ1qg/PorfPWVGX585BEIC4MlSyAhweqEIiJigaxZISTkvx+ubDoNCQlh165dFC9enPvuu++mR7Zs2bjvvvvImDEjmzZtSvyZs2fPsm/fvju+59atW7l27RoTJkygatWqlCpViqNHj6Y8pIVU1FmkQgWYOxf++APatIHRo6F4cXjtNTh3zup0SeBwQNOmsGkTrFwJ2bNDy5bw4IPmKLL4eKsTioiIl+nTpw9nzpyhffv2bN68mT///JPvv/+eJ598kvj4eLJnz06PHj0YOnQoq1atYufOnXTr1g0fnzuXOyVLluTatWt88MEH/Pnnn8yZM4cpU6ak4adyHxV1FitRAsLD4c8/zQDYW2+Z4m7kSDh71up0SeBwQMOGZl55/XqzU7ZDByhTBqZPNyN5IiIiblCoUCE2bNhAfHw8jRs3ply5cvTv3x9/f//Ewm3cuHHUrl2bZs2a0bBhQ2rWrEloaOgd3zM4OJh33nmHt99+m3LlyjF37lzGjBmTVh/JrRxOT5sQtrno6Gj8/f05f/48OXPmTPbPR0WZAx+mTIGMGaF/fxg40Cxps42ICDP0uHjx/59U0b27+UAiImJ7ly5dYs+ePZQpU8Yrm/imhuu/Z7///jsnTpygS5cu5MmTB3C9drhOI3UepmBBeOcdM3LXsyeMH29G80aPhhv6IHq20FBYtAh++w2qV4fevc3K2NmzNS0rIiKSSlTUeah77oEJE+DgQTMtO2oU3Huv6S5y5YrV6ZKoXDmzvm7HDrOIsGtX89yCBdpQISIi4mYq6jxcQIA5wWv/fmjWDIYMMa1Qpk2zSZ87MMeNLVoEW7eaYce2baFiRbN7VrP/IiIibqGiziaKFoWPP4bdu6FmTXjqKQgKgvnzbVQXhYaafi4//QR585rz06pUMUeSiYiIpJK6desyYMAAq2OkOhV1NlOq1P/PaJYpA+3ambpo7VqrkyVDjRrw44+mmPP1NbtnH37YfCgRERFJERV1NvXgg2b2cs0a83XdumZ6dvduK1MlU/36sHEjfPmlWTxYsaJZd3f4sNXJREREbEdFnc3VqQO//GKmYXftMsvXnnrKtEaxBYcDHn8cdu6ESZNgxQozHDl0qE0a9YmISErVrVuXvn37MmDAAHLnzk1AQABTp04lJiaG7t27kyNHDkqWLMm3336b+DNr166lcuXKZMqUiYIFC/LCCy9w7dq1xO/HxMTQpUsXsmfPTsGCBZkwYcIt142Li2PYsGEULlyYbNmyUaVKFdZcHyWxMRV1bhIeHk5QUBBhYWFpfm2Hw5xKsXu3aYfy5Zdw333wyitw4UKax0mZjBlN65M//oDhw2HyZLPdd9w4G233FRGR5Jo1axb58uVj8+bN9O3bl2eeeYbWrVtTvXp1tm3bRuPGjencuTOXLl3in3/+oUmTJoSFhbFjxw4mT57MJ598whtvvJH4fkOHDmX16tUsXryY77//njVr1hAREXHTNbt3786GDRv4/PPP+fXXX2ndujUPP/ww+/fvT+uP71ZqPuxm7mog6Ipz58zJFO++C7lywZtvQrducJdTUjzPsWOmj8vUqVC4sPlA7dqZClZERCx12+bDly7B3r3//cMPPJB4AGzdunWJj49n/fr1AMTHx+Pv78/jjz/O7NmzATh27BgFCxbk559/5uuvv+bLL79kz549OP7398GkSZN4/vnnOX/+PJcuXSJv3rzMnj2btm3bAnDmzBmKFCnCU089xbvvvsuBAwe4//77+fvvvylUqFBirIYNG1K5cmXefPNNd/023UTNhyVFcuUyNdC+fWbZWo8eULkybNhgdbJkuOceMx27e7dZa9ehg9lgsXmz1clEROR29u41XQ7+6/Gvwu/BBx9M/LWvry958+alfPnyic8FBAQAcOLECfbs2UO1atUSCzqAGjVqcPHiRf7++28OHDhAXFwc1apVS/x+njx5KF26dOLX27Ztw+l0UqpUKbJnz574WLt2LQcOHHD7b0taymB1AEk9RYvCvHnQp485bqxmTTPY9fbb5nu2UKoULFlidssOHGi2+nbubIYfixSxOp2IiFz3wAPmmMikvO4GGf91hKTD4bjpuesFXEJCAk6n86aCDuD6hKPD4SApk48JCQn4+voSERGBr6/vTd/Lnj37f+f3YCrq0oHrA1yzZsGLL5o/T8OGmYdtjuyrXx+2bYPp081Zsl9+aT7A0KE2+hAiIl4sa1YICUnVSwQFBfHll1/eVNxt3LiRHDlyULhwYXLnzk3GjBnZtGkTRf83enH27Fn27dtHnTp1AKhYsSLx8fGcOHGCWrVqpWretKbp13TCxwe6dzdTsv37w5gxULq06Xlnm1WVvr7Qq5c5XuO558xoXenSMHeujT6EiIik1LPPPsuRI0fo27cve/fuZenSpYwcOZJBgwbh4+ND9uzZ6dGjB0OHDmXVqlXs3LmTbt264XPDovJSpUrRsWNHunTpwqJFizh48CBbtmzh7bffZvny5RZ+OtepqEtncuQwBd3u3RAWZpaq1ayZtBFzj+Hvb+aQd+82iwU7dYJatdS8WETEyxUuXJjly5ezefNmKlSoQO/evenRowcvvfRS4mvGjRtH7dq1adasGQ0bNqRmzZqEhobe9D4zZsygS5cuDB48mNKlS9OsWTN++eUXAgMD0/ojuZV2v7qZJ+x+TY4ff4QBA0ybuN694Y034H+bcezjxx/NyN3vv8Ozz8Lrr5vdIiIikipuu/tV7kq7XyXVXV+qNnEifPqpmc2cMQMSEqxOlgz165tRurFjYeZMs7li+nSbfQgRERHXqKgTMmQw6+x+/x0aN4YnnzSzmZGRVidLhowZYfBg8yEeesj0calRw2bzyiIiIimnok4SFSxoRutWrzYNjENDoV8/OH/e6mTJUKiQ2Tixdi3ExJiFg888A2fOWJ1MREQkVamok1vUrWtG6d5+28xili4Nc+bYbINp7dr/P688b56Zkp06FeLjrU4mIiKSKlTUyW1lzAhDhpjG33XqQJcu5n937bI6WTLcOK/86KPw9NNQtaqmZEVExCupqJO7KlIE5s+HlSvhxAkIDja9fy9ftjpZMtxzj+m8/NNPEBtr2qAMGAAXLlidTERExG1U1EmSNGxoNpi+/DKMHw/lyplCz1aub5x46y2YNg3KlDFHkImIiHgBFXWSZJkywSuvwK+/QrFi0KgRdOxoRvBsI2NGc7TYrl1m2LFlS2jeHA4ftjqZiIiIS1TUSbKVLg2rVpmWcCtWmLNkP/7YZm3hiheHr7+GhQth61YICoJ33oFr16xOJiIikiIZrA7giZYtW8bgwYNJSEjg+eefp2fPnlZH8jgOB3TtavYfDBlijmSdNQs++sjUR7bgcECrVqav3UsvmQ8yZ475EJUrW51ORMTjXbbVAmtrpcXvlY4J+5dr164RFBTE6tWryZkzJyEhIfzyyy+JR3n8F7sdE+Yuq1ebzaWHDpmNFMOHg5+f1amSaetW8yG2bzfHjb35JqSjeygiklRxcXHs2rWLBFtN0VjP6XSye/duzpw5kyrHhGmk7l82b95M2bJlKVy4MABNmjRhxYoVtG/f3uJknq1ePbPWbvRoc37sF1/AJ59AlSpWJ0uGSpXgl1/gww/NyN3SpTB5Mjz2mNXJREQ8ip+fH2XLluXatWucP3+eZcuWkS1bNjJlymR1NI8WHx/P1atXU+39va6oW7duHePGjSMiIoKoqCgWL15MixYtbnrNpEmTGDduHFFRUZQtW5Z3332XWrVqAXD06NHEgg6gSJEi/PPPP2n5EWwrc2Z4/XVo3dqc0lWtmmkT98YbkC2b1emSKEMG0+7k8cehd29o2hTatYP33oMCBaxOJyLiMfz8/PDz8+PKlSvExMRw9epVFXVJEBsbm2rv7XVFXUxMDBUqVKB79+60atXqlu/Pnz+fAQMGMGnSJGrUqMFHH33EI488wu7duylatCi3m412OBxpEd1rPPgg/PyzqYNeftl0DZk2zbRFsY2iReGbb8xpFP37m/Yn774LnTqZtXgiIgJApkyZyJkzJ9HR0cTFxVkdxxZy5syZKgWwV6+pczgct4zUValShZCQECZPnpz4XJkyZWjRogVjxoxh48aNjBs3jsWLFwPQv39/qlSpQocOHZJ0zfS6pu5ODhwwmyhWr4bu3WHCBMid2+pUyXTyJAwcaM6UbdwYpkwxu2dFRAQwAyqpOQLlbTJlykS2G6awtKYuBeLi4oiIiOCFF1646flGjRqxceNGACpXrszOnTv5559/yJkzJ8uXL+eVV16543vGxsbe9H/k6Ojo1AlvUyVLmvYnn3wCgwfD8uUQHm42ndpG/vzw6afQoYOZki1XziwefO458PW1Op2IiOWyZct2U5Ei1khXfepOnTpFfHw8AQEBNz0fEBDAsWPHAMiQIQMTJkygXr16VKxYkaFDh5I3b947vueYMWPw9/dPfAQGBqbqZ7AjhwN69oTdu83Rq088YZasRUVZnSyZmjQxTYu7dzcjdzVqwM6dVqcSEREB0llRd92/18g5nc6bnmvWrBn79u3jjz/+4Kmnnrrrew0fPpzz588nPo4cOZIqmb1B4cKweDEsWAAbNph+drNmga0WAOTIAR98YM6RjY6GkBAzapeKu5lERESSIl0Vdfny5cPX1zdxVO66EydO3DJ6l1TXF4je+JA7czjM7tjdu03j4m7doFkzOHrU6mTJVL266Wc3dCiMHGmGIH/7zepUIiKSjqWros7Pz4/Q0FBW/usk+pUrV1K9enWX3js8PJygoCDCwsJcep/0Im9es0xtyRLYsgXKloXZs202apcpkxml27QJYmMhNNT0b9GonYiIWMDrirqLFy8SGRlJZGQkAAcPHiQyMpLD/zuwfdCgQXz88cdMnz6dPXv2MHDgQA4fPkzv3r1dum6fPn3YvXs3W7ZscfUjpCvNm///qF3XrjYdtatUCSIiYNgwePVV03H511+tTiUiIumM17U0WbNmDfXq1bvl+a5duzJz5kzANB8eO3YsUVFRlCtXjokTJ1K7dm23XF8tTVJu6VJzSldsLLz/vk1bwkVEmDnl3383p1IMHw4ZM1qdSkREPJi7agevK+qspqLONadPQ79+pufvY4/BRx9BoUJWp0qm2FgzDTtmDJQvDzNnQoUKVqcSEREP5a7aweumX62iNXXukTev6fG7ePH/r7WbM8eGa+1efx02b4b4eDM9+9proE7rIiKSijRS52YaqXOfG0ftmjY1o3YFC1qdKpni4sxmijffNBXqzJkQHGx1KhER8SAaqROvd+Oo3ebNpq/dp5/abNTOz8+M0m3ebIKHhZnNFBq1ExERN1NRJx6vRQtzkEOTJtC5szmR4uRJq1MlU8WKZj55xAgzchcWBjt2WJ1KRES8iIo6sYXro3ZffAFr15rjV7/6yupUyeTnZ0bprre9CQszmymuXbM0loiIeAcVdW6ijRJp44knzHGrlSubHnc9epjTumwlONhMxw4ebNqe1K4N+/dbnUpERGxOGyXcTBsl0obTCTNmQP/+ZhRv5kyoW9fqVCmwcSN06QJRUTBuHDzzjA2b84mIiCu0UULSNYcDnnzSHNxQrBjUqweDBsHly1YnS6bq1c3aum7doE8faNwY/v7b6lQiImJDKurE1kqUgNWrYcIEmDTJHL8aEWF1qmTKlg3Cw2HFCnNmWrlyNmzOJyIiVlNRJ7bn42NG6SIiIEsWqFoVRo2Cq1etTpZMjRrBb7+Zpnxduth0m6+IiFhFRZ2baKOE9cqWhU2bzHGro0ZBjRqwd6/VqZIpd24zSrdw4f9v81261OpUIiJiA9oo4WbaKOEZNm82Pe0OH4a334bnnjMjerZy7Bg89RR8/bVZc/fuu+Dvb3UqERFxM22UELmLypVh+3ZTE/XvDw89ZAo8W7nnHjNKN306fPklPPgg/Pij1alERMRDqagTr5U1K7z3HvzwA+zbB+XL2/CYMYcDunc3a+3uvRcaNICBA+HKFauTiYiIh1FRJ16vQYP/33/QuTO0bw9nz1qdKpmKFYNVq+Cdd2DyZHMaxa+/Wp1KREQ8iIo6SRdy5TKjdJ99ZjqH2HIm08fHjNJt2WJG8MLCTC+XhASrk4mIiAdQUecm2v1qD+3amQGu++83I3hDhkBsrNWpkql8ebMTpG9f8wEaNoQjR6xOJSIiFtPuVzfT7ld7SEiAiRPhxRehdGmYO9fUSrbz44+mp11MjJmWbdfO6kQiIpJM2v0q4gIfHxg82MxkOp1mJvPdd204k1m/vlkw2LixWSzYqROcO2d1KhERsYCKOknXHnzQFHbPPmuWqzVqZMOjV3PnNosFP/3U9LSrUME0LhYRkXRFRZ2ke5kzm02lK1fCnj2m0PviC6tTJZPDAR07mgWDxYtDvXrwwgsQF2d1MhERSSMq6kT+p2FDM5PZoAG0aQNdu0J0tNWpkqlYMbPO7q23TKVapQrs3m11KhERSQMq6kRukCcPLFgAs2bB4sVmJvOnn6xOlUy+vjBsGPzyi9naGxoKH3xgs67LIiKSXCrq3EQtTbyHw2E2lO7YAYULQ506MGKEDWcyK1aEiAjo1Qv69YNHHoGjR61OJSIiqUQtTdxMLU28S3w8vP02jBwJwcEwb57pcWc7331njhu7ehWmToXHH7c6kYiI/I9amoikAV9f08tu40bTKaRiRZgxw4YzmQ8/bBYM1q4NrVrBk0/ChQtWpxIRETdSUSeSBGFhsH07tG1r6qF27Wx4fmy+fPDllzB9utneGxxs1t2JiIhXUFEnkkTZs8Mnn8D8+eb82AoVYP16q1Mlk8NhpmEjI02RV6MGjB5t5plFRMTWVNSJJFObNv/fDq5uXXj5ZbNUzVZKljTbel94wXyA+vV1fqyIiM2pqBNJgaJFYfVqeO01GDMGatWCP/+0OlUyZcwIb7xhPsiff5quywsXWp1KRERSSEWdSAr5+sJLL5kBrxMnzBK1Tz+1OlUK1Kljhh4bNoTWraFnT7h40epUIiKSTCrqRFxUtapZotaiBXTubE7rOn/e6lTJlDu36br8ySfmHNmQENi61epUIiKSDCrqRNwgZ06YPRvmzoVly8yo3caNVqdKJofDbO3dvt18oGrVYOxYSEiwOpmIiCSBijo30YkSAtChgxm1K1TItIQbNQquXbM6VTKVKmUq0kGDzEaKhx6Cf/6xOpWIiPwHnSjhZjpRQsAUcm+8Aa+/bga85s6FYsWsTpUCq1aZM9OuXDFTsy1aWJ1IRMTr6EQJEQ+WIQO8+iqsXQt//2162n3+udWpUqBBA7OJonZtaNkSeveGS5esTiUiIrehok4kFdWsaaZjH3kE2reHrl1teDpX3rywaBFMmWIWDoaGmg8lIiIeRUWdSCrLlQvmzYNZs0xtVLEibN5sdapkcjjg6achIgIyZYIqVWDiRG2iEBHxICrqRNKAw2GWpkVGmoGvGjXgzTdteDpXmTLmvNg+fcxGikcegagoq1OJiAgq6kTS1PXTuZ5/3jQubtDArLmzlUyZ4J134LvvYMcOcxLFsmVWpxIRSfdU1ImksRtP5zpwwGyiWLrU6lQp0Lix2URRpQo0bQp9+5pdsiIiYgkVdSIWqVPHTMfWqmU6hdiyJipQAL7+Gj74AKZNMwXenj1WpxIRSZfSpKi7ePEix44dIy4uLi0uJ2IbefPC4sUQHm7jmsjhgOeeM2vt4uKgUiXT004tMEVE0lSqFXVffvklzZo1w9/fH39/fwoXLkyWLFnImzcv9evXZ9GiRal1aRFbcTjg2WfNjtirV21cE1WoYM6L7dABevY0PVzOnbM6lYhIupEqJ0pMmDCBZcuW0bJlS4oVK0a2bNnw8/Pj4sWLXLhwgf3797No0SI6duzI4MGD3X15S+lECXFFTAwMHGhG7dq2hY8+An9/q1OlwIIF8NRTkDu36edSrZrViUREPJa7aocMbsyU6Pjx46xevfqurxkxYoTXFXQirsqWDaZONcet9uoFwcHw2WdQtarVyZKpTRuoXNmM1tWqZQ7Bff558PW1OpmIiNdKlenXggUL/udrHA5Hkl4nkh61bm02URQsaE6lGDPGhn1+ixeHdev+v39Lo0Zw9KjVqUREvFaqFHUHDx5k1KhR7Ny5kwsXLpBww99GsbGx7NmzhxdffJFdu3alxuUtER4eTlBQEGFhYVZHES9RvLg5O/b552HECJvWRBkzwujR8MMPZgeIetqJiKSaVFlTd+XKFYYPH860adO4fPkyAL7/m3aJj4/H39+fVq1aMX78ePxtuWDozrSmTlLDjz9Cp05mI8WsWdCkidWJUuDUKeje3RR1/frB2LGmkbGISDrnrtohVYq6665cucLevXs5fvw4Z86cIUeOHBQsWJDg4ODEIs/bqKiT1HLypKmJvvnGbKYYM8aGNZHTaXraDR0KQUHw+edQurTVqURELGWLoi49UlEnqcnphPffh2HDoFw5s4miVCmrU6VAZCS0awdHjsCHH0K3bqa3i4hIOuSu2sGta+r27NlDv379+Pnnn935tiLyPw4H9O8PmzbBxYsQEgKzZ9uwp11wMEREmMLuySdNb7vz561OJSJia24t6kaPHs2UKVPo1avXTc9/8skn9OzZk79td3K5iGeqWNHURK1bQ9eu0LkzREdbnSqZsmUzXZY/+wyWLzcf6pdfrE4lImJbbi3qcufOzTfffMP48eNver5Hjx689dZbDBs2jJ07d7rzkiLpVvbsMGMGzJ0LX31lRu22bLE6VQq0awfbt5tzZGvWhLfesmH/FhER67m1qHM4HJQqVYqHH374lu/ly5ePadOmMXbsWHdeUiTd69DB1ER58kD16jB+vA1ronvvhfXrzQaKF1+Exo0hKsrqVCIituLWom7cuHEMGDCAQYMG8f333xMTE3PT968fFyYi7lWyJPz0EwwaZOqiJk3g+HGrUyVTxozw5puwciXs3Gl62i1fbnUqERHbcOvu1xkzZtCzZ0+uv2WGDBkIDQ2lbt26VK5cmYsXLzJ79mxWrlzprkt6HO1+Fat9/z106WJ+PXu2aVpsOydPmh2xy5fbuH+LiEjSeOTZrwsWLODgwYNkzpyZrVu38uOPP7J69WrGjRtHQkIChQsXZtGiRe68pIj8S6NGsGOH2UDRuLEZuXvjDbDVIHn+/KZJ8Xvvmf4ta9fauH+LiEjacOv063333UfRokUpUKAATZo0Yfz48URERHDq1Cnmzp1LpUqVNHolkgYCAswg17hxMHGi2X9w4IDVqZLJ4YABA0z/lgsX/r9/i4iI3JZbi7oMGTJw7NixW57PlSsX7du357PPPmPMmDHuvKSI3IGPDwwZAhs3wpkzpmPIvHlWp0qBkBDYtu3/+7d06WKa9ImIyE3cWtSNHDmSQYMG8cMPP9zyvddff52xY8eSPXt2d15SRP5DWJipiZo1g44dTa/ff+1h8nzX+7fMmQOLFplCLzLS6lQiIh7FrUVdrly5mDVrFnv27GHu3Lk3fW/ZsmW8+uqrxMXFufOSIpIEOXPCp5/CzJkwf74p9H77zepUKdCpk6lQs2eHKlXMEWO2O05DRCR1pNnZr6dOnWLt2rU0btzYq0frtPtVPN3evdC2LezbB+++C089ZcNjV2NjzQaK99+HFi3MyRR58lidSkQkRTzy7Ne7yZcvH61atfLqgk7EDh54wJzG9eST0Lu3KfDOnbM6VTJlymR2xi5ZYnbGVqwIGzZYnUpExFJpVtTZScuWLcmdOzdPPPGE1VFEUkXmzBAeDgsXmr52tj12tXlzs7YuMBDq1DHNi+PjrU4lImIJFXW30a9fP2ardYKkA61amZronntM25Nx42x4xFjRorBmDbzwArz0kmnOd5td+CIi3k5F3W3Uq1ePHDlyWB1DJE0ULw7r1sHgwWaZ2mOPmQMdbCVDBtNheeVK2LULKlQwQ5AiIumI7Yq6devW0bRpUwoVKoTD4WDJkiW3vGbSpEmUKFGCzJkzExoayvr169M+qIiNZMwIb70F330HW7eammj1aqtTpUCDBmboMTjYjNi98AJcvWp1KhGRNGG7oi4mJoYKFSrw4Ycf3vb78+fPZ8CAAYwYMYLt27dTq1YtHnnkEQ4fPpz4mtDQUMqVK3fL4+jRo2n1MUQ8UuPG5oixMmVMfTRyJFy7ZnWqZAoIgG+/hbffhgkToHZtOHTI6lQiIqkuzVqapAaHw8HixYtp0aJF4nNVqlQhJCSEyZMnJz5XpkwZWrRokazTLNasWcOHH37IwoUL7/q62NhYYmNjE7+Ojo4mMDBQLU3E1uLjYcwYU9TVrAlz50KRIlanSoFNm6BdOzh/3rQ9efxxqxOJiNzCdi1N0kJcXBwRERE0atTopucbNWrExo0bU+WaY8aMwd/fP/ERGBiYKtcRSUu+vmbPwZo18OefZjZz2TKrU6VA1apmOrZBA7MrpE8fuHLF6lQiIqnCq4q6U6dOER8fT0BAwE3PBwQE3PZM2jtp3LgxrVu3Zvny5RQpUoQtW7bc8bXDhw/n/PnziY8jR46kOL+Ip6lVy9RE1atD06YwaBDY7lCYXLngiy9g0iQzWleliunALCLiZTK4641WrFjBd999x59//snFixe506yuw+Fg1apV7rrsHa9xI6fTectzd7NixYokvzZTpkxkypQpya8XsZu8eWHpUnN4w9ChsH49fP45lCxpdbJkcDjgmWegRg3TbTk01DTq69rVhsdpiIjcnstFXXR0NC1atGDt2rV3LORulJziKrny5cuHr6/vLaNyJ06cuGX0zt3Cw8MJDw8nXo1PxQs5HNC/v1lf17ataVY8dapZrmYrDz5otvf27Qvdu8OqVWYETy2MRMQLuFzUPf/886xZs4Y8efLw1FNPUbFiRfLnz5+qxdud+Pn5ERoaysqVK2nZsmXi8ytXrqR58+apeu0+ffrQp0+fxMWOIt4oNBS2bTODXu3bm5rovfcga1arkyVDtmwwfbpZZ9e7tzlKY/58U6mKiNiYy0XdokWLyJgxI2vXrqVs2bLuyHRXFy9e5I8//kj8+uDBg0RGRpInTx6KFi3KoEGD6Ny5M5UqVaJatWpMnTqVw4cP07t371TPJpIe5MwJn35qaqLnnoOffzY1URr88Xevjh2hcmUz3Fi1qjlOo29fTceKiG253NIke/bs3Hvvvfz666/uynRXa9asoV69erc837VrV2bOnAmY5sNjx44lKiqKcuXKMXHiRGrXrp0m+dy1LVnEDnbvNtOxBw6YEbuePW1YE8XGmqM03n/fnCU7fTrkyWN1KhFJR9xVO7hc1FWqVInz58+zf/9+V97G9m5cU7dv3z4VdZJuXL4MAwaYNXZt25r/teX/9b/6yqyzy5YN5s0zCwhFRNKAx/Sp69OnDwcOHGDNmjWuvpWt9enTh927d9+1/YmIN8qSBT76yEzBfvutWZq2davVqVKgWTPTv6VYMahbF0aPNl2YRURswuWirnv37vTt25fHH3+cDz74gIsXL7ojl4jYTJs2sH27aYFSvTpMnAi2O68mMNAcejt8OLz8sjk3LSrK6lQiIknilmPCYmNjad++PUuXLgUgf/78ZL3DdjiHw8GBAwdcvaTH0po6Se/i4uDFF82xq489BjNmQL58VqdKgR9/NJsp4uNhzhxT4ImIpAKPWVN3/PhxGjZsyO7du5Pcp86be7mpqBMxvvnG9PbNnNksUUujvUrudeIEdOkCK1aYzRRvvAEZM1qdSkS8jLtqB7f0qdu1axf33XcfQ4cOJTg42LI+dVZS82GRmz36KOzYYQa76tWDkSNhxAhzrqxtFCgAy5ebYccXX4S1a+Gzz6BECauTiYjcwuWRunvuuYfo6Gj++OMPChUq5K5ctqWROpGbxcfD66+bR+3aMHcu2PI/FZs2mY7LZ8/Cxx/DE09YnUhEvITH7H6NiYnhgQceUEEnIrfl6wuvvmpOn9i3DypUMLtkbadqVbMT5KGHoHVrePZZuHLF6lQiIolcLurKly/P6dOn3ZFFRLxY3bqmY0jlytCkCQwdajZV2EquXLBgAUyebJoUV60Kv/9udSoREcANRd3QoUM5cuQICxYscEceEfFi+fPD11+bJWrvvQe1asHBg1anSiaH4//PjL1yxRyI++mnVqcSEXG9qGvZsiXvv/8+PXv2ZPDgwezatYsr6XBKIjw8nKCgIMLCwqyOIuLRfHxg0CDYsAFOnjTNihcutDpVClSoYLosP/44dO4MTz4JMTFWpxKRdMzljRK+ydzK5nA4uHbtmiuX9GjaKCGSdOfPQ69e8MUXZonahAmmBYrtzJplPkCxYmZ6tlw5qxOJiI14zEYJp9OZrEdCQoKrlxQRL+Hvb44XmzIFPvnELFHbt8/qVCnQtasZtfP1hbAwszvWdsdpiIjduVzUJSQkJPshInKdwwFPP22WqF2+bJaozZtndaoUKFMGNm82U7G9epkGfRcuWJ1KRNIRl4s6ERF3qFABIiKgRQtTD/XsCZcuWZ0qmbJkgalTTVX69dcQEmLaoIiIpAG3F3X79u1j2bJlfPbZZyxbtox9tpxLERErZM8Os2ebbiHz5pn2J7t2WZ0qBdq3h23bIEcOM6ccHq7pWBFJdW4r6j766CPuvfdeypQpQ/PmzenUqRPNmzenTJkylCxZkmnTprnrUh5Ju19F3MPhgO7dzRI1MEvUZsywYU10//3w889mbvm558wJFOfOWZ1KRLyYy7tfAbp3787s2bNxOp1kypSJwMBAAgICOH78OEeOHCE2NhaHw0GXLl2YMWOGO3J7LO1+FXGfS5egf3+z76BTJ5g0yQx+2c7ixablSa5c8PnnUKWK1YlExIN4zO7XefPmMWvWLLJmzcrYsWM5efIk+/btY/369ezbt4+TJ08yduxYsmXLxuzZs/nss89cvaSIpBNZs8K0aea82CVLoFIl2LHD6lQp0LKlWVsXEAA1a8L48aBNYyLiZi4XddOmTcPhcPDll18yZMgQsmfPftP3s2fPzpAhQ1i4cCFOp9Prp2FFxP06dDCbKLJmNYNcU6bYcDq2eHFYvx4GDjRnpDVtCqdOWZ1KRLyIy9OvefLkIW/evOzfv/8/X1uqVClOnjzJ2bNnXbmkR9P0q0jquXIFhgwx+w5atzajeP7+VqdKgW+/hS5dIFMmsyOkdm2rE4mIhTxm+vXKlSvkypUrSa/NmTMnsbGxrl5SRNKpzJnhww/NsWLff286hmzZYnWqFHjkEYiMhJIloV49eOMNiI+3OpWI2JzLRV3RokXZuXMnp/5jGuHkyZPs2rWLokWLunpJEUnnWrUyHUPy5oUaNeDdd204HVu4MKxaBS+9BK+8Ag8/DMeOWZ1KRGzM5aKuWbNmxMbG0rZtW06ePHnb15w4cYK2bdsSFxdH8+bNXb2kiAj33gs//QR9+5plai1awJkzVqdKpgwZ4LXXYOVK+O03CA6GH36wOpWI2JTLa+rOnDlDcHAw//zzD5kyZaJ169YEBQVRoEABTpw4we7du/niiy+4cuUKgYGBbN++nTx58rgrv8cIDw8nPDyc+Ph49u3bpzV1Imno66+hWzfIls10DKle3epEKXD8uOnbsmoVjBgBI0eaok9EvJ671tS5pU/dH3/8Qfv27YmIiDBv6nAkfu/624eFhTFv3jxKlizp6uU8mjZKiFjjyBFzkMOmTTB6tNlg6mO3gxATEuCtt+Dll8288rx5UKSI1alEJJV5VFF33apVq/j+++/Zt28fFy9eJHv27JQqVYrGjRtTv359d13Go6moE7HOtWtmedqYMdC4sTlyrEABq1OlwPr1pkK9cgVmzYJHH7U6kYikIo8p6g4fPgxAkSJF8LHdP4vdT0WdiPVWrIDOnc3s5WefQZ06VidKgVOnzJzyN9+YPi6jR4Ofn9WpRCQVeExLk+LFi1NFR96IiAdp3NicPPHAA1C/PowaZcOOIfnymcWCEyaY7b21a8OhQ1anEhEP5nJR5+/vT7FixTRKJyIepWBBs6l05EizwbRRI4iKsjpVMjkcMGgQbNhgNlIEB8OiRVanEhEP5XIlVr58+cQpWBERT+Lra9bYrVoFe/aYmmjlSqtTpUDlyubs2IYNTZO+vn3NejsRkRu4XNT179+fY8eOMX36dHfkERFxu7p1zQEOFSuaqdkRI8ymClvJlQu++MKckTZ1qunbkoTjGUUk/XC5qGvVqhVvvfUWffr0YeDAgWzbto3Lly+7I5uIiNsUKADLl8Obb8Lbb5tC78gRq1Mlk8MBzz4Lv/wCFy+ac9I++8zqVCLiIVze/err65u8CzocXLPdP5H/m5oPi9jHhg2mY0hMjOkY8thjVidKgQsX4JlnYO5c6NkT3nsPsma1OpWIpIDHtDRJyQaJhIQEVy7p0dTSRMQeTp+G7t3NBtNBg0xvO9t1DHE6YcYMeO45c27aggUQFGR1KhFJJo9paZKQkJDsh4iI1fLmhaVLYeJE+OADqFULDh60OlUyORzw5JOwZYsp8CpVMkWe+3rKi4iNJKuoW7duHTt27EitLCIiacrhgAEDzHTsyZNmI4UtO4aULWsKuw4dTJHXpYuZnhWRdCVZRV3dunXp16/fTc/Vr1+fAQMGuDOTiEiaCgszHUMeesjGHUOyZoWPP4ZPP4UlS8yonf4RLpKuJKuoczgct0yfrlmzhm3btrk1lIhIWvP3N0vSJk2CadNs3DGkY0eIiDBFXpUqMGWKpmNF0olkFXW5cuXir7/+Sq0sIiKWcjjMhtJNm2zeMaRUKfj5Z7Mr9plnoG1bOH/e6lQiksoyJOfF1apV49tvv6VVq1Y0atSILFmyAHDixAlmz56d5Pfp0qVL8lKKiKSh4GAz2PXMM2aZ2o8/2rBjSObM8OGHUK8e9OhhFgzOn2/mmkXEKyWrpcmvv/5K3bp1OXfuHA6HAwCn05n466SKt93J2kmnliYi3uPGjiElS5qayJYdQw4ehHbtzMLBt982u0OS+d9tEUk97qodkjVS9+CDD/L777/z+eefs3fvXi5fvszMmTMpUKAADz/8cIpDiIh4ousdQ6pUgTZtzCBXeDh07WqzmqhECVi/Hl580TTlW73aVKt581qdTETcyC3Nh2vWrMm6devclcnWNFIn4p0uXYJ+/eCTT6BzZ7OhInt2q1OlwLJlpirNmhU+/xxq1LA6kUi65zHNh0eOHEn37t1dfRsREY92Y8eQRYts3DHkscdM8OLFoU4dc5SGmsKLeAWXR+rkZhqpE/F++/aZ6di9e80Giqeestl0LMC1a/Dqq/Dmm6ZB35w5UKCA1alE0iWPGakTEUlvSpUybU969IDevc0eBNt1DMmQAd54A1asgMhIs+V3zRqLQ4mIK1TUuUl4eDhBQUGEqV2ASLqQObPZNPHFF/Ddd6an3datVqdKgYceMkXdAw9Agwbw2mvgxR0KRLyZpl/dTNOvIunPn3+a/r47dsD48eaYMdtNx8bHw+jRpqirWxfmzoV77rE6lUi6oOlXEREPce+9sGGD6WfXvz+0bAlnzlidKpl8feGVV2DVKti9GypUgB9+sDqViCSDijoRETfw84N33oGlS2HdOnOAw88/W50qBerWNUOOwcHQqBG89JLZVCEiHk9FnYiIGzVrZpaoFSkCtWrB2LE27BhSoAB8+63ZSDFmjFlr988/VqcSkf+gok5ExM2KFjUbSYcOheefN63hTp60OlUy+fiYEyjWrIEDB8zI3bffWp1KRO7CrUXdkSNHmDdvHuPGjWPUqFE3fe/q1avExcW583IiIh4rY0YzyPXtt2ZXbHAwrF1rdaoUqFXLDD1WrgxNmpgq9epVq1OJyG24pag7deoUbdu2pUSJEnTu3JkXXniB11577abXdO/enSxZshAREeGOS4qI2MLDD5uaqFQpqF8fXn/dhh1D8uWDr7+GcePMwsG6deHwYatTici/uFzUXbhwgTp16vDFF19QuHBhunXrRuHChW95Xc+ePXE6nSxatMjVS4qI2EqhQmYj6csvw8iR0LgxHDtmdapk8vGBIUPMLpC//zZDj19/bXUqEbmBy0Xd2LFj2bNnD61atWLv3r188sknFCtW7JbX1a5dmyxZsrB69WpXLykiYju+vuZUrh9+gF27TE20apXVqVKgWjXYvt1MyzZrBoMHg5bWiHgEl4u6hQsXkilTJj7++GOyZMly5wv5+HDfffdxWEP2IpKO1a9vpmPLlzeHOYwcacPp2Dx5YMkSePdd+OADU+AdPGh1KpF0z+Wi7tChQ5QqVQp/f///fG3WrFk5deqUq5cUEbG1gABztNjrr5uuIQ0awNGjVqdKJofDdFresMFs7a1YEbS8RsRSLhd1mTNn5sKFC0l6bVRUVJKKPxERb+frCyNGwOrVsH+/mY5dscLqVCkQFgbbtkHDhtCqlTkjLTbW6lQi6ZLLRV3ZsmU5cuQIf/31111fFxkZyeHDhwkNDXX1kiIiXqN2bTMdGxpqdsq++KIND3DIlQu++ALCw2HqVKheHf74w+pUIumOy0Vdp06diI+P56mnnuLSpUu3fc3Zs2fp0aMHDoeDLl26uHpJERGvkj8/fPMNvPWWOYGiXj2zwdRWHA549lnYtAmioyEkBBYssDqVSLriclHXq1cvatWqxcqVKylfvjwvvPACx48fB2D69OkMGjSI0qVLs337dh566CHatWvncmgREW/j42P6+q5dC4cOmenY5cutTpUCFStCRAQ8+ii0bQu9e8Ply1anEkkXHE6n0+nqm1y4cIGnnnqK+fPn43A4uP6WN/66TZs2fPLJJ2TLls3Vy3m06Oho/P39OX/+PDlz5rQ6jojY0OnT0K0bLFtmjhobPdqcUGErTid8/DH062c6Ly9YAKVLW51KxCO5q3ZwS1F33W+//cbixYv57bffOH/+PNmzZycoKIiWLVvaZi3dkSNH6Ny5MydOnCBDhgy8/PLLtG7dOsk/r6JORNzB6YSJE83oXaVK8PnncJsWoJ7v11+hTRsznzxlCnTqZHUiEY/jkUWdN4iKiuL48eMEBwdz4sQJQkJC+P3335M8wqiiTkTc6ZdfzCxmdDTMmAHNm1udKAUuXjTr7ebMgSefNL3tsma1OpWIx3BX7eCWs1+9ScGCBQkODgagQIEC5MmThzNnzlgbSkTSrSpVzAEOdepAixYwcKAND3DInh1mzzZV6WefQeXKsHu31alEvI7LRd1XX33Fvffey4QJE+76ugkTJnDvvfey3MWVv+vWraNp06YUKlQIh8PBkiVLbnnNpEmTKFGiBJkzZyY0NJT169en6Fpbt24lISGBwMBAlzKLiLgid27T1/e990zXkJo1bXqAQ7dusHWr+XWlSqbI02SRiNu4XNTNnj2bv/76i5YtW971dc2bN+fQoUPMnj3bpevFxMRQoUIFPvzww9t+f/78+QwYMIARI0awfft2atWqxSOPPHLT8WShoaGUK1fulsfRG1q6nz59mi5dujB16lSX8oqIuIPDYfYcbNxoNlLY9gCHoCDYvBk6dDBTsV27mulZEXGZy2vqSpYsyaVLl4iKivrP1xYsWJBs2bLxh5uaUjocDhYvXkyLFi0Sn6tSpQohISFMnjw58bkyZcrQokULxowZk6T3jY2N5aGHHqJXr1507tz5P18be0P39OjoaAIDA7WmTkRSzfnz0LMnLFwIzz0H48ZB5sxWp0qBuXPh6achMBDmz4cHH7Q6kYglPGZN3dGjRylatGiSXhsYGJik4i+l4uLiiIiIoFGjRjc936hRIzZu3Jik93A6nXTr1o369ev/Z0EHMGbMGPz9/RMfmqoVkdTm7286hEyaBNOm2fgAh44dTU87Pz+zeHDqVE3HirjA5aIuW7ZsnDx5MkmvPXXqFJkyZXL1knd9//j4eAICAm56PiAggGPHjiXpPTZs2MD8+fNZsmQJwcHBBAcH89tvv93x9cOHD+f8+fOJjyNHjrj0GUREksLhgGeeMQc4XLxoDnCYP9/qVClQurT5EN26mVG7Dh3MVl8RSbYMrr5B+fLlWbduHVu3bqVSpUp3fN3WrVs5dOgQNWvWdPWS/8nhcNz0tdPpvOW5O6lZsyYJCQlJvlamTJlStVAVEbmb4GAz2PX009CuHaxebfrbZclidbJkyJIFJk+GunWhVy9zEO6CBWbhoIgkmcsjdR06dMDpdNKxY0f+/PPP277m4MGDdOzYEYfDQYcOHVy95B3ly5cPX1/fW0blTpw4ccvonYiIt8iRwyxPmzoVZs2CqlXh99+tTpUCbdvCtm2QM6f5EOHhmo4VSQaXi7onn3yS6tWrs3//fsqVK0enTp344IMPmDNnDh988AEdO3akXLly7N+/n2rVqtGrVy935L4tPz8/QkNDWbly5U3Pr1y5kurVq6fadQHCw8MJCgoiLCwsVa8jInI7DocZ5PrlF4iNNYNdc+danSoF7rvPbPF9+mmzC6RNGzh3zupUIrbglhMlzp07R/fu3Vm6dKl50xumOq+/fcuWLfnkk0/IlSuXS9e6ePFi4u7ZihUr8s4771CvXj3y5MlD0aJFmT9/Pp07d2bKlClUq1aNqVOnMm3aNHbt2kWxNDhjRydKiIjVbjzAoUcPeP99mx7gsGiRaXuSJ49ZMKh/NIuX8shjwrZu3crSpUvZs2cP0dHR5MiRg7Jly9KiRQtCQkLcco01a9ZQr169W57v2rUrM2fOBEzz4bFjxxIVFUW5cuWYOHEitWvXdsv1/4uKOhHxBE4nzJwJffrAvfeaJWpBQVanSoGDB81iwe3bYexY6N/fDEuKeBGPLOpERZ2IeJbdu6F1azh0yCxR69bN6kQpEBcHw4fDO+9As2bmJIo8eaxOJeI2HtOnTgytqRMRTxQUBFu2mMGu7t1teoCDnx9MmABffQXr15tdsT//bHUqEY/j9pG6s2fPcvHiRe72tkltVmxHGqkTEU81Z47pbRcYaKZjy5e3OlEKHD4M7dubo8befBMGDwYfjU+IvXnU9Ou+fft49dVX+e677zh//vzdL+hwcO3aNVcv6bFU1ImIJ9u712wo3b8fPvjAbKSw3RK1q1fh5Zfh7behSRPTxyVfPqtTiaSYu2oHl5sPR0ZGUqdOncTRucyZM5M/f3589C8nERGP88ADpu3JgAGmBcrq1TBliul1ZxsZM8Jbb5lmxZ07mw7Mn30GtWpZnUzEUi5XXi+++CIXLlygfv36/Prrr1y6dIm//vqLgwcP3vHhjbSmTkTsIksW+OgjUwd99ZXpaRcZaXWqFHj4YRO8ZEmoV89MxybjRCARb+Py9GuuXLlISEggKiqKbNmyuSuXbWn6VUTsZP9+c5DD7t3meLHevW04HXvtGrz2GoweDQ0bwqefQoECVqcSSTKP2f2akJBA6dKlVdCJiNjQ/febAxx69jQNi9u1g/9YGu15MmSA11+H77+HHTugQgUzryySzrhc1AUHBxMVFeWOLCIiYoHMmeHDD+GLL+C77yAkBCIirE6VAg0bmunYoCDz69deg/h4q1OJpBmXi7rhw4cTFRXFnDlz3JFHREQs8sQT5uCGPHmgenWzO9Z27ekLFjQjdiNHwqhR0KgRaOBB0gmXi7pHHnmESZMm8eyzzzJw4EB27tzJ5cuX3ZHNVrRRQkS8wb33wk8/mX52/fpBq1Zw9qzVqZLJ1xdeeQVWrYI9e8zu2JUrrU4lkupc3ijh6+ubvAuqT52IiC0sWWJOociVC+bPh8qVrU6UAidOmLYnK1fCiy/Cq6+aNXgiHsRjNko4nc5kPRK03VxExBZatDBL1AICoEYNc/Sq7aZjCxSAb781O2Pfegvq14e//7Y6lUiqcMvu1+Q+RETEHooVg3XroH9/cyJX8+Zw5ozVqZLJxweGD4c1a+DPP8107LffWp1KxO107IOIiNyVnx+MHw9ffw0bNpiaaONGq1OlQM2aZuixalVzvNjzz5sjx0S8hIo6ERFJksceMzVR0aJQuzaMHWvDAxzy5TPHaIwbZ+aT69SBw4etTiXiFirq3ES7X0UkPQgMNH19hw41A12PPQYnT1qdKpl8fGDIEFi/Ho4eNUOPX31ldSoRl7m8+/W6mJgYvv76a3bs2MGZM2e4eochbYfDwSeffOKOS3ok7X4VkfTiu+/MxlI/P/j8c6hVy+pEKXD2rNniu3QpDBxoNlP4+VmdStIZd9UObinqPv/8c5555hmio6MTn7v+to4bDhF0Op04HA7ivbjDt4o6EUlP/vkHOnQwve1GjTL7EXzsNgfkdML775vhx+Bg07+lRAmrU0k64jEtTX7++Wc6d+5MfHw8I0aM4L777gNg2rRpvPLKKzRr1gyHw0HmzJkZPXo006dPd/WSIiLiIQoXNj1+R4yAl1+Ghx+G48etTpVMDofZ3rtxI5w+DRUrwqJFVqcSSTaXR+patWrFkiVLWLJkCU2bNqVWrVps3LjxptG4vXv30rp1a86ePUtERAQBAQEuB/dUGqkTkfTqhx+gUydTI82bB/XqWZ0oBc6fh549YeFCeO45s6Eic2arU4mX86iRunz58tG0adM7vuaBBx7gyy+/JCoqipEjR7p6SRER8UANG5rdsUFB5tevvQa2W23j7w8LFsCkSTBtmjkE948/rE4lkiQuF3WnT5+maNGiiV/7/W+BaUxMzE2vK1WqFGXLluVbNXwUEfFa99wD338PI0eaNXaNGkFUlNWpksnhMIffbtoEFy9CSIhZZyfi4Vwu6vLmzcvly5cTv86XLx8ABw4cuOW18fHxHLfdYgsREUkOX1945RWz1m7PHrP3YOVKq1OlQHAwRESYvi3t2sHTT8MNf9+JeBqXi7rixYsTdcM/w0JCQnA6ncydO/em1+3YsYN9+/aRP39+Vy/pkdSnTkTkZnXrmunY4GBo3BheegmuXbM4VHLlyAFz55qp2NmzzWkUv/9udSqR23K5qHvooYc4d+4cu3btAqBDhw5kzpyZ8ePH06lTJ8LDw3nllVdo0KABCQkJtGrVyuXQnqhPnz7s3r2bLVu2WB1FRMRjFChgjlkdPdq0gGvQwLRBsRWHw2ye2LwZ4uIgNBQ+/dTqVCK3cHn3665duxgwYADPPPMMjz/+OACzZs3iqaee4urVq4l96pxOJ1WrVuX7778ne/bsrif3UNr9KiJyez/9ZGYxY2NhzhzT/sR2Ll6EPn3MqF337vDhh5A1q9WpxOY8qvnw7fz5558sWLCAQ4cOkSVLFmrWrEmLFi3w9fVNjct5DBV1IiJ3duoUdO0Ky5ebY8Zefx0yZrQ6VQrMnGmKu+LFzW7ZsmWtTiQ25vFFXXqlok5E5O4SEmDCBHjxRahc2RwxFhhodaoU2L0b2rSBP/+E8HDo1s1M1Yokk8f0qVu3bh07duxI0mt//fVX1q1b5+olRUTExnx8zIlc69bB33+bjRRff211qhQICjLr7Dp2hCefhC5dzPSsiEVcLurq1q1Lv379kvTa/v37U79+fVcvKSIiXqBaNdi+HWrWhGbNYPBgsw/BVrJmNTtj586FJUugUiX49VerU0k65ZZjl5Mzg6vZXhERuS5PHlMLTZwIH3wAtWvDoUNWp0qBDh1MT7vMmc2c8tSpoL/vJI25pahLqtOnT5MlS5a0vKSIiHg4hwMGDIANG+D4cahY0RR6tlOqlDmFont306i4QweIjrY6laQjGZL7A9HR0Zw7d+6m52JjYzly5MgdR+EuX77M2rVr2blzJxUqVEhRUBER8W5hYWY6tkcPaNkS+vWDsWMhUyarkyVD5swweTLUq2d624WGmiPGQkKsTibpQLKLuokTJzJq1Kibntu6dSvFixdP0s/36NEjuZe0hfDwcMLDw4m33enVIiKeI1cuWLjQbCYdPNiM3s2fDyVLWp0smdq0MQVd27Zm8eCECaYFinbHSipKdkuT9957j3fffTfx68OHD+Pn58c999xz+ws4HGTJkoV7772Xtm3b0qlTJ5cCezq1NBERcY9t20xtdPIkfPwxtG5tdaIUiI01W30/+ABatTIfJFcuq1OJh/GYPnU+Pj7UrFlTrUr+R0WdiIj7REdDr16mv+8zz8A775gZTttZvNi0Pcmd2ww96pxwuYHH9KmbMWMGL774oqtvIyIicoucOU1z4ilTYPp0M5O5f7/VqVKgZUuzYDB/fqhRA959V7tjxe1S9USJ6Ohovv32W44ePUpISAh16tRJrUt5DI3UiYikjh07zHTs0aOmY0j79lYnSoG4OHOUxoQJpjnfjBmmr4ukax4zUjd//nxCQkL4+OOPb3p+7969lCtXjg4dOjBkyBDq169Pt27dXL2ciIikUxUqwNat0Ly56Rby1FNw+bLVqZLJzw/Gj4evvoKffjLHafz8s9WpxEu4pajbsWMHtWvXvun5AQMG8Pfff3PvvffSvHlzsmfPzpw5c1i+fLmrlxQRkXQqRw6YM8fsN/j0U9Pnd+9eq1OlQNOmEBlpDr2tVcv0bklIsDqV2JzLRd2OHTvIkycPpUqVSnwuKiqKlStXUrRoUX777TcWLVrE119/jdPpJDw83NVLiohIOuZwmF52mzdDfLzpHDJ7ttWpUiAwENasMbtjn38eHnsMTp2yOpXYmMtF3cmTJylatOhNz61evRqn00mHDh3I/L9tSrVr16ZYsWLs2bPH1UuKiIhQrhxs2WJanXTtag5yiImxOlUyZcwIY8bAt9+aDxMcDOvXW51KbMrloi4uLu6Whrvr16/H4XBQr169m54PCAggKirK1UuKiIgAkC0bzJxpHgsWmOnYXbusTpUCDz9spmNLloS6dWH0aE3HSrK5XNQVLlyYAwcOcOnSpcTnvvvuOzJkyECNGjVueu2FCxfw9/d39ZIiIiI36drVbKLw8TEt4KZPt2HHkMKFYdUqGDECXn7ZFHonTlidSmzE5aKuYcOGXLp0ib59+7Jz505effVV/vrrL+rXr0/WrFkTX3f58mX2799PYGCgq5cUERG5RZky8Msv0LGjWXPXpQtcvGh1qmTKkAFGjYLvvzc9XIKDzbo7kSRwuagbMWIEefLkYebMmVSoUIFRo0aRMWNGXnvttZte9/XXX3Pt2jVq1arl6iVFRERuK2tWmDYN5s6FJUvMJopff7U6VQo0bGimYx94ABo0MIWezhaX/+ByUVe0aFG2bt3Ks88+S6NGjejZsyebN2+mcuXKN71uzZo1VKhQgebNm7t6SY8UHh5OUFAQYTr6RUTEch06QEQEZMli1tlNnWrD6diCBWHlSjMV++qr0LgxHDtmdSrxYKl6okR6pBMlREQ8x5UrMHCgOWasXTv46CNz9Jjt/PijmVd2OmHePKhf3+pE4kYec6KEiIiIp8qcGSZPhvnz4ZtvzHTs9u1Wp0qB+vXNdGz58mZqduRITcfKLZI1Unf48GEAMmbMSMGCBW96Ljn+3dfOm2ikTkTEMx04YM6O3bkT3nkHnn3WNDK2lfh4ePNNMx1bu7YZtfvf38diX+6qHZJV1Pn4+OBwOHjggQfY9b9GQNefS/IFHQ6uXbuW/KQ2oaJORMRzxcaaAxw++ACeeMJsqsiVy+pUKbB2LbRvb4q8Tz+Fhx6yOpG4wF21Q4bkvLho0aI4HI7EUbobnxMREfF0mTLB+++b/r5PPgkhIWZq1nZ73OrUMdOxnTubDRQvvmhG7zIk66918TLaKOFmGqkTEbGHgwehbVtTG40bB/362XA6NiEB3nrL7JCtWdNMxxYubHUqSSZtlBAREXFBiRLw00/w3HMwYAC0bAlnzlidKpl8fMwo3erV8Mcfplnxd99ZnUosoqJORETSLT8/s2li6VJYtw4qVoRNm6xOlQK1a5shx7AweOQRGD4cvHj9utxeina/ukq7X0VExNMcPmx62W3ZAmPGwKBBZiDMVhISzFzyiBFQtSp89hnoeE6PZ+nuV1do96uIiHiqq1fhpZdg7Fh49FGYNQvy5rU6VQps2GAq1EuXYPZs82HEY1mypq5o0aJ3fPj6+uJ0OnE6nfj6+hIQEECGDBkSn8uQIQNFixYlUP9iEBERD5UxI7z9tmlUvGmTWaL2009Wp0qBGjXMdGz16vDYY6aPy9WrVqeSVJasou7QoUMcPHjwlsejjz6Kw+GgX79+7N27l9jYWI4ePcqVK1f4/fff6devHw6Hg8cee4yDBw+m1mcRERFxiyZNTE1UvLhpf/LWW2Zm01by5oWvvoLx4+Hdd826u7/+sjqVpCKXW5pMmjSJvn378tlnn9GmTZs7vm7BggW0b9+eDz/8kGeeecaVS3o0Tb+KiHiPa9fMiVxjxkCjRmYms0ABq1OlwKZNZjo2OhpmzIDmza1OJDewZE3d7VSoUIHo6OgkjcCVKFECf39/IiMjXbmkR1NRJyLifb7/Hjp1Mr19P/vM9P61nbNnoXt3s9V3wAAzz+znZ3UqwYP61P3xxx/kz58/Sa/Nnz8/+/fvd/WSIiIiaapRI9ixA0qXhvr14fXXzQldtpI7NyxebKZiw8NNs2ItifIqLhd12bNnZ9euXZw7d+6urzt37hy7du0iW7Zsrl5SREQkzRUsCD/8YA5vGDnSnM517JjVqZLJ4YD+/c3u2FOnTGO+RYusTiVu4nJR99BDD3H58mU6duzImTu04j579iwdO3bkypUrNG7c2NVLpqoLFy4QFhZGcHAw5cuXZ9q0aVZHEhERD+Hra45Y/eEH2LXL7I5dtcrqVCkQFgbbtkHDhtCqFfTtC7GxVqcSF7m8pu7w4cOEhIRw9uxZsmTJQuvWrSlTpgz58+fn5MmT7N27ly+++IKYmBjy5s3L1q1bKVasmLvyu118fDyxsbFkzZqVS5cuUa5cObZs2ULeJDYq0po6EZH04fhxs85u1SozevfKK6bosxWnEyZNMp2Wy5eH+fOhZEmrU6U7HrNRAmDPnj106tSJ7du3mze9oUHx9bevWLEic+bMISgoyNXLpZkzZ85QsWJFIiIiyJcvX5J+RkWdiEj6ER9vdsaOHGk6hsydC4UKWZ0qBbZtgzZt4ORJ+PhjaN3a6kTpisdslAAoU6YMERER/PDDDwwdOpRmzZpRv359mjVrxtChQ1m5ciURERFuKejWrVtH06ZNKVSoEA6HgyVLltzymkmTJlGiRAkyZ85MaGgo69evT9Y1zp07R4UKFShSpAjDhg1LckEnIiLpi6+vOYHixx9h3z4zHbtihdWpUiAkxBR2Dz9sirtnn4UrV6xOJcmUwZ1vVr9+ferXr+/Ot7xFTEwMFSpUoHv37rRq1eqW78+fP58BAwYwadIkatSowUcffcQjjzzC7t27E8+cDQ0NJfY2awe+//57ChUqRK5cudixYwfHjx/n8ccf54knniAgICBVP5eIiNhXnTqmWXGXLqYuGj4cRo0yLVBsI2dO+PxzqFfPtDz5+WdYsADuv9/qZJJEbpl+tYrD4WDx4sW0aNEi8bkqVaoQEhLC5MmTE58rU6YMLVq0YMyYMcm+xjPPPEP9+vVpfYeh6NjY2JsKxOjoaAIDAzX9KiKSDiUkwLhxMGIEVKtmetoVKWJ1qhSIjDQjdlFRMHUqtG9vdSKv5lHTr54iLi6OiIgIGjVqdNPzjRo1YuPGjUl6j+PHjxMdHQ2Y3+R169ZRunTpO75+zJgx+Pv7Jz50tq2ISPrl4wPPPw9r18KhQ2Y6dvlyq1OlQHAwRERAs2bQoQM89RRcvmx1KvkPXlXUnTp1ivj4+FumSgMCAjiWxGZCf//9N7Vr16ZChQrUrFmT5557jgcffPCOrx8+fDjnz59PfBw5csSlzyAiIvZXo4YZ7KpWDR59FIYNg6tXrU6VTDlywKefwrRpMGcOVKkCe/danUruwk6z/Ul24+5bMDtw//3cnYSGhibrGLNMmTKRKVOm5MQTEZF0IG9e+OoreOcdeOEFWL/eLFnz4K5et3I4oGdPU9C1aQOVKsHkydC5s9XJ5Da8aqQuX758+Pr63jIqd+LECW10EBGRNOdwwODBpqCLijIHOCxdanWqFChfHrZsMY2Ku3SBJ5+ES5esTiX/4lVFnZ+fH6GhoaxcufKm51euXEn16tVT9drh4eEEBQURFhaWqtcRERH7qVoVtm83u2RbtICBAyEuzupUyZQ9O8yaBTNmmCHHsDDYvdvqVHID2xV1Fy9eJDIyMnGK9ODBg0RGRnL48GEABg0axMcff8z06dPZs2cPAwcO5PDhw/Tu3TtVc/Xp04fdu3ezZcuWVL2OiIjYU+7c5pjV996D8HCoWRMOHrQ6VQp06wZbt5pfV6oEM2damUZuYLuWJmvWrKFevXq3PN+1a1dm/u//WJMmTWLs2LFERUVRrlw5Jk6cSO3atdMkn06UEBGR/7J1K7RtC6dPw/Tp8PjjVidKgUuXzJmx06ebKdnwcDOaJ8nmUceEyf9TUSciIklx/rzZg7BwITz3HIwfD7bcdzdnDjzzDAQGmmbF5ctbnch21KfOw2hNnYiIJIe/v6mBJk0yXUOqV4c//rA6VQp07myGHjNmhMqVzdmxGi+yhEbq3EwjdSIiklzXD3A4dswUeG3bWp0oBS5fhv79zQfo0AGmTDG97uQ/aaRORETES1w/wOGxx6BdO+jd24YHOGTJYo4UmzfPNOirVAl27LA6Vbqiok5ERMQD5MgBc+eaumjWLNMG5fffrU6VAu3bw7ZtkDWraVo8ZYqmY9OIijo30Zo6ERFxlcMBvXrBL79AbCyEhppCz3buvx9+/hl69DCbKNq1g/+dqy6pR2vq3Exr6kRExB0uXoRnnzWbS3v0gPffN4NftvPFF2abb/78ZmdISIjViTyO1tSJiIh4sesHOEyfbpapVakCe/ZYnSoFWrc207G5ckG1avDhh5qOTSUq6kRERDyUwwHdu5tjVxMSzN6DWbOsTpUCJUvChg3w9NOmYXHr1nDunNWpvI6KOhEREQ9Xtixs3mxanXTrZh4xMVanSqZMmcwc8pdfwg8/mGlYHa3pVirqREREbCBbNjMVO3u2OYUiLAx27rQ6VQo8/jhs3w758kGNGuYwXE3HuoWKOjfR7lcREUkL1w9wyJDBFHaffGLDmqhECfjpJ3M+2oABptA7e9bqVLan3a9upt2vIiKSFm48wKFjR5g82aYHOHz1lZlPzpkT5s83O0LSGe1+FRERScduPMBh6VIbH+DQrJmZji1YEGrWhAkTbDj06BlU1ImIiNhY+/bmiLEsWcwg10cf2bAmKlYM1q2DgQNhyBBT6J0+bXUq21FRJyIiYnOlSsGmTfDkk+bc2PbtbXiAQ8aMMHYsLFtmTqOoWBE2brQ6la2oqBMREfECmTPDpElmWdry5eaIsW3brE6VAo8+CpGRULQo1K4Nb79tmvTJf1JR5yba/SoiIp6gTRuzRM3f38YHOBQpAmvWwLBh8MIL8NhjcPKk1ak8nna/upl2v4qIiCeIjYWhQ+GDD6BVK/j4Y3NSl+2sWAGdOoGfH3z+OdSqZXUit9PuVxEREbkjrznAoXFjMx17331Qty6MHq3p2DtQUSciIuLFvOIAh8KFYdUqePFFePllePhhOHHC6lQeR0WdiIiIl/OKAxwyZIDXX4fvvzcN+YKDzbo7SaSiTkREJB3w84N33jGNiteuNR1DfvnF6lQp0LChmY594AFo0ABGjYL4eKtTeQQVdSIiIunI9QMcChX6/wMcbLdErWBBWLkSXnkFXn0VGjWCY8esTmU5FXVuopYmIiJiF8WKmdG6QYNsfICDry+MHGnW2u3ebaZjV62yOpWl1NLEzdTSRERE7GT5cujSxRwz9vnnZjOF7Rw/Dp07m22+L71kij1fX6tTJZlamoiIiIjLmjQxS9SKF4c6deCtt2w4HRsQAN99ZzZSjB5t1todPWp1qjSnok5ERCSdK1IEVq82BzgMH25O6rLdAQ4+PjBihPkg+/eb6dgVK6xOlaZU1ImIiAgZMsCbb5oBr4gIUxOtW2d1qhSoXdsMPYaGmn52w4fDtWtWp0oTKupEREQk0fUDHO6/H+rVs+kBDvnzwzffmLnkcePMB/n7b6tTpToVdSIiInKTQoX+f8/B9QMcjh+3OlUy+fjA88+bbb6HDpmhx+XLrU6VqlTUiYiIyC0yZIDXXjPt4H791dREP/5odaoUqFHDDD1WrWoWCw4bBlevWp0qVaioExERkTtq0MDURGXLmsMcXn3Vhgc45M0LX31lpmInTjTbfA8ftjqV26moExERkbu65x6zkfS110zXkIcegqgoq1Mlk4+P6bS8fj38848ZevzqK6tTuZWKOjfRiRIiIuLNfH3N+rpVq2DvXlMTrVxpdaoUqFrVnJNWqxY0b26O1YiLszqVW+hECTfTiRIiIuLtTpwwBzisXAkvvmimZDNksDpVMjmd8N57Zo1dxYrmOI0SJSyJohMlRERExBIFCsC338Ibb8CYMWbd3T//WJ0qmRwOGDAANmwwVWrFirB4sdWpXKKiTkRERJLNx8eM0q1ZAwcOmOnY776zOlUKhIWZ6dgGDeDxx6FfP4iNtTpViqioExERkRSrVcvsjq1cGR55BF54wYYdQ3LlgoUL4YMP4KOPzO5YGzYrVlEnIiIiLsmXD77+GsaOhfHjoW5dOHLE6lTJ5HDAc8/BTz+Zrb0hIbZrzKeiTkRERFzm4wNDh5qOIUeOmOnYZcusTpUCYWHm8NsKFUzvlrFjbXNOmoo6ERERcZtq1cx0bI0a0LSpaQ1nu+nYfPnMAsHnnzeP+vXhzz+tTvWfVNSJiIiIW+XJA0uXwjvvmK4htWrBX39ZnSqZfH3hzTfNFOyhQ2bkbu5cq1PdlYo6ERERcTuHAwYONEvUjh0z07FLl1qdKgXq1YPffoOWLaFTJ+jZEy5dsjrVbamoExERkVRTpYrpGFKvHrRoYVrD2e4Ahxw5YNYsmD4d5s0zW31377Y61S1U1ImIiEiqyp0bvvwS3n8fJk826+1ssETtZg4HdO8OW7ea0ygqVYIZM8yvPYSKOhEREUl1Dgf07QsbN8KZM+YAhy+/tDpVCgQFwZYt0L49PPkkdO0KFy9anQpQUSciIiJpKDQUtm2Dxo3hiSdMa7grV6xOlUxZs8Inn8CcObBokWmD8uuvVqdSUecu4eHhBAUFERYWZnUUERERj+bvD/Pnw6RJ8PHHUL06/PGH1alSoFMn09POz88sHpw61dLpWIfT6UGTwV4gOjoaf39/zp8/T86cOa2OIyIi4tEiI6FNG7NDdto0aNvW6kQpcPkyDBoEU6ZAu3bmqLFk1ADuqh00UiciIiKWCQ42g12PPWbqod69TY1kK1mymB0g8+fDN9+YOebt29M8hoo6ERERsVSOHKav79SppnNI1arw++9Wp0qBNm3MgsGcOc2HCA9P0+lYFXUiIiJiOYcDevWCX36B2Fgz2OXhBzjc3n33mS2+Tz9tdoG0bg3nzqXJpVXUiYiIiMd48EHTCu7xxz3+AIc7y5TJNOVbtAhWrTL9WzZvTvXLqqgTERERj5I9+80HOFSpAnv2WJ0qBVq2NGvrChSAmjVh4sRUnY5VUSciIiIe5/oBDlu2QEKCOcBh1iyrU6VA8eKwfj3062d2yLZoYbovpwIVdSIiIuKxypY1M5dt20K3buYRE2N1qmTy84Px4+Grr+Cnn0yFmgpzyirqRERExKNly2amYmfPhoULzQEOO3danSoFmjY1jflefdWcSuFmKupERETEFjp3NpsoMmQwhd0nn1h6gEPKBAZCly6p8tYq6kRERMQ2HnjAtD3p3NnsjO3cGS5csDqVZ1BRJyIiIraSJYtpVDxvHixdapao7dhhdSrrqagTERERW2rf3hwxliWLaXvy0Uc2nI51IxV1IiIiYlulSsGmTfDkk+bc2PbtITra6lTWUFEnIiIitpY5M0yaBPPnw/Ll5oix7dutTpX2VNSJiIiIV2jTxhRz/v5QtSqEh6ev6VgVdbdx6dIlihUrxpAhQ6yOIiIiIslQsiRs2ABPPw3PPQetW8O5c1anShsq6m5j9OjRVKlSxeoYIiIikgKZMsH778OXX8IPP0BIiDluzNupqPuX/fv3s3fvXpo0aWJ1FBEREXHB44+b6dh8+aBGDXjvPe+ejrVVUbdu3TqaNm1KoUKFcDgcLFmy5JbXTJo0iRIlSpA5c2ZCQ0NZv359sq4xZMgQxowZ46bEIiIiYqUSJcxxq889BwMGmELv7FmrU6UOWxV1MTExVKhQgQ8//PC2358/fz4DBgxgxIgRbN++nVq1avHII49w+PDhxNeEhoZSrly5Wx5Hjx5l6dKllCpVilKlSqXVRxIREZFU5ucH77xjGhWvXQsVK5pTKbyNw+m050Ckw+Fg8eLFtGjRIvG5KlWqEBISwuTJkxOfK1OmDC1atEjS6Nvw4cP59NNP8fX15eLFi1y9epXBgwfzyiuv3PFnYmNjiY2NTfw6OjqawMBAzp8/T86cOVP24URERCRV/PUXtGtnzpB96y0YNAgcDmszRUdH4+/v73LtYKuRuruJi4sjIiKCRo0a3fR8o0aN2LhxY5LeY8yYMRw5coRDhw4xfvx4evXqddeC7vrP+Pv7Jz4CAwNT/BlEREQkdRUrBuvWmWJuyBBo1gxOn7Y6lXt4TVF36tQp4uPjCQgIuOn5gIAAjh07lmrXHT58OOfPn098HDlyJNWuJSIiIq7LmBHefhu++QZ+/hmCg00bFLvLYHUAd3P8awzV6XTe8lxSdOvWLUmvy5QpE5kyZUr2+4uIiIi1mjSByEhztFidOvDGGzBsGPjYdMjLprFvlS9fPnx9fW8ZlTtx4sQto3epITw8nKCgIMLCwlL9WiIiIuIeRYrA6tWmmBs+HB59FE6etDpVynhNUefn50doaCgrV6686fmVK1dSvXr1VL9+nz592L17N1vSQ3dDERERL5IhA7z5Jnz3HUREmOnYdeusTpV8tirqLl68SGRkJJGRkQAcPHiQyMjIxJYlgwYN4uOPP2b69Ons2bOHgQMHcvjwYXr37m1hahEREbGDxo3NdOz990O9ejB6NCQkWJ0q6Wy1pm7r1q3Uq1cv8etBgwYB0LVrV2bOnEnbtm05ffo0o0aNIioqinLlyrF8+XKKFStmVWQRERGxkUKFzNFio0bByy+bvnZz5kAarORymW371Hma8PBwwsPDiY+PZ9++fepTJyIiYnM//ACdOpk+dnPnQv36qXMdd/WpU1HnZu66MSIiImK9Y8dMYffjj/DKK2b0ztfXvddQ82ERERGRVHbPPbBiBbz2Grz+Ojz0EERFWZ3q9lTUiYiIiNyFr68ZoVu1CvbuNbtj/9VswyOoqBMRERFJgrp1ze7Y4GCzU/all+DaNYtD3UBFnZuo+bCIiIj3K1AAvv3WnD4xZgw0aAD//GN1KkMbJdxMGyVERETSh/XrzRFjsbGm7cnDD6fsfbRRQkRERMRCtWqZ6diwMHjkEXjhBbh61bo8KupEREREUihfPli2DMaOhfHjzbq7I0esyaKiTkRERMQFPj4wdKiZjj1yxGykWLbMghxpf0nvpI0SIiIi6Vu1amY6tkYNaNoUhgxJ2+lYbZRwM22UEBERSd+cTnj3XRg2DEJDYf58uNsx9NooISIiIuKBHA4YOBB++skcMxYcDEuXpv51VdSJiIiIpIIqVWD7dqhXD1q0gAEDIC4u9a6nok5EREQkleTODV9+Ce+/D5Mnm/V2f/6ZOtdSUSciIiKSihwO6NsXNm6EM2egenW4dMn918ng/rdMn8LDwwkPDyc+Pt7qKCIiIuKBQkNh2zbYvBmyZnX/+2v3q5tp96uIiIgkh7tqB43Uudn1Gjk6OtriJCIiImIH12sGV8fZVNS52YULFwAIDAy0OImIiIjYyYULF/D390/xz2v61c0SEhIoVaoUEREROBwOoqOjCQwM5MiRI5ZMx4aFhbFlyxbL3iupP/Nfr7vb9+/0vds9/+/nvOX+eMO9+ffzVt+bf+dJ6/fxtPujPzvJ/5mkvE73x/3vY8c/O06nkwsXLlCoUCF8fFK+h1UjdW7m4+ODn5/fLZV2zpw5LfmD5evr67brpuS9kvoz//W6u33/Tt+73fN3eq3d74833Js7PW/VvblTnrR6H0+7P/qzk/yfScrrdH/c/z52/bPjygjddWppkgr69OljdYRE7sySkvdK6s/81+vu9v07fe92z3vSvQH35fGGe5OcTGlF9yf5edKKN9yb/3qN7o/+7CSXpl9TmXbDejbdH8+le+PZdH88m+6P50rNe6ORulSWKVMmRo4cSaZMmayOIreh++O5dG88m+6PZ9P98VypeW80UiciIiLiBTRSJyIiIuIFVNSJiIiIeAEVdSIiIiJeQEWdiIiIiBdQUWexZcuWUbp0ae6//34+/vhjq+PIDVq2bEnu3Ll54oknrI4i/3LkyBHq1q1LUFAQDz74IF988YXVkeR/Lly4QFhYGMHBwZQvX55p06ZZHUlu49KlSxQrVowhQ4ZYHUX+JUOGDAQHBxMcHEzPnj2T9bPa/Wqha9euERQUxOrVq8mZMychISH88ssv5MmTx+poAqxevZqLFy8ya9YsFi5caHUcuUFUVBTHjx8nODiYEydOEBISwu+//062bNmsjpbuxcfHExsbS9asWbl06RLlypVjy5Yt5M2b1+pocoMRI0awf/9+ihYtyvjx462OIzfIly8fp06dStHPaqTOQps3b6Zs2bIULlyYHDly0KRJE1asWGF1LPmfevXqkSNHDqtjyG0ULFiQ4OBgAAoUKECePHk4c+aMtaEEMEciZc2aFYArV64QHx+Pxg48y/79+9m7dy9NmjSxOoq4mYo6F6xbt46mTZtSqFAhHA4HS5YsueU1kyZNokSJEmTOnJnQ0FDWr1+f+L2jR49SuHDhxK+LFCnCP//8kxbRvZ6r90ZSlzvvz9atW0lISCAwMDCVU6cP7rg3586do0KFChQpUoRhw4aRL1++NErv/dxxf4YMGcKYMWPSKHH64o77Ex0dTWhoKDVr1mTt2rXJur6KOhfExMRQoUIFPvzww9t+f/78+QwYMIARI0awfft2atWqxSOPPMLhw4cBbvuvV4fDkaqZ0wtX742kLnfdn9OnT9OlSxemTp2aFrHTBXfcm1y5crFjxw4OHjzIvHnzOH78eFrF93qu3p+lS5dSqlQpSpUqlZax0w13/Pk5dOgQERERTJkyhS5duhAdHZ30AE5xC8C5ePHim56rXLmys3fv3jc998ADDzhfeOEFp9PpdG7YsMHZokWLxO/169fPOXfu3FTPmt6k5N5ct3r1amerVq1SO2K6ltL7c+XKFWetWrWcs2fPTouY6ZIrf3au6927t3PBggWpFTFdS8n9eeGFF5xFihRxFitWzJk3b15nzpw5na+99lpaRU5X3PHn5+GHH3Zu2bIlydfUSF0qiYuLIyIigkaNGt30fKNGjdi4cSMAlStXZufOnfzzzz9cuHCB5cuX07hxYyvipitJuTdinaTcH6fTSbdu3ahfvz6dO3e2Ima6lJR7c/z48cSRhejoaNatW0fp0qXTPGt6lJT7M2bMGI4cOcKhQ4cYP348vXr14pVXXrEibrqTlPtz9uxZYmNjAfj777/ZvXs39957b5KvkcF9ceVGp06dIj4+noCAgJueDwgI4NixY4DZtjxhwgTq1atHQkICw4YN0w6xNJCUewPQuHFjtm3bRkxMDEWKFGHx4sWEhYWlddx0Jyn3Z8OGDcyfP58HH3wwcc3KnDlzKF++fFrHTVeScm/+/vtvevTogdPpxOl08txzz/Hggw9aETfdSep/28QaSbk/e/bs4emnn8bHxweHw8F7772XrI4YKupS2b/XyDmdzpuea9asGc2aNUvrWMJ/3xvtRLbW3e5PzZo1SUhIsCKWcPd7ExoaSmRkpAWp5Lr/+m/bdd26dUujRHKju92f6tWr89tvv6X4vTX9mkry5cuHr6/vLf86OnHixC1VuqQt3RvPpvvjuXRvPJvuj2dLi/ujoi6V+Pn5ERoaysqVK296fuXKlVSvXt2iVAK6N55O98dz6d54Nt0fz5YW90fTry64ePEif/zxR+LXBw8eJDIykjx58lC0aFEGDRpE586dqVSpEtWqVWPq1KkcPnyY3r17W5g6fdC98Wy6P55L98az6f54NsvvT5L3ycotVq9e7QRueXTt2jXxNeHh4c5ixYo5/fz8nCEhIc61a9daFzgd0b3xbLo/nkv3xrPp/ng2q++Pzn4VERER8QJaUyciIiLiBVTUiYiIiHgBFXUiIiIiXkBFnYiIiIgXUFEnIiIi4gVU1ImIiIh4ARV1IiIiIl5ARZ2IiIiIF1BRJyIiIuIFVNSJiIiIeAEVdSIiIiJeQEWdiIiIiBdQUSciIiLiBVTUiYiIiHgBFXUiIiIiXkBFnYiIiIgXyGB1ABGR9O7IkSNMmTKFuLg4zp07R8OGDWnbtq3VsUTEZlTUiYhYaPLkybz77rssXLiQ8uXLc/nyZRo0aMCvv/7K6NGjAdi0aRMZM2YkNDTU4rQi4sk0/SoiYpEPP/yQfv36MW/ePMqXLw9AlixZGDZsGG+99Ra///47ADNmzCAwMNDKqCJiAyrqREQssHfvXoYOHUqPHj1uGYGrXLkyCQkJLFq0iPj4eE6fPk2BAgUsSioidqGiTkTEAhMmTODKlSv07t37lu/lz58fgEOHDjF9+nQ6d+6c1vFExIYcTqfTaXUIERE7uXDhAk888QSxsbHJ+rm33nqLqlWrApA7d24yZ85MVFTUbV/rcDho1aoVMTExLF++HIfD4XJuEfFuKupERNLY6dOnyZcvH02aNOGbb7657WscDgdZsmRhw4YNVKxYMY0TiogdafpVRCSNZcyYEeCu6+QcDgdNmzZVQSciSaaiTkQkjeXMmZPg4GBOnTp12+9PnToVPz8/EhISALh69WpaxhMRm1JRJyJigYkTJ7J27VoiIyMTn9u6dSs9evQgR44c9OrVi3379uF0OunQoQOXL1+2LqyI2ILW1ImIWGTr1q1MmjQJX19ffH19CQ4OpmPHjuTIkYMLFy4wePBgLl++zNNPP03NmjWtjisiHk5FnYiIiIgX0PSriIiIiBdQUSciIiLiBVTUiYiIiHgBFXUiIiIiXkBFnYiIiIgXUFEnIiIi4gVU1ImIiIh4ARV1IiIiIl5ARZ2IiIiIF1BRJyIiIuIFVNSJiIiIeAEVdSIiIiJeQEWdiIiIiBdQUSciIiLiBf4PlLWULAK32kYAAAAASUVORK5CYII=", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "execution_count": 8, "metadata": { }, "output_type": "execute_result" } ], "source": [ "from functools import partial\n", "from scipy.optimize import minimize_scalar\n", "\n", "\n", "XMAX=100_000\n", "SCALE=\"loglog\"\n", "\n", "\n", "def statistical_indicators():\n", " x0_val = CONSTANT_VALUES[x0]\n", " q_espression = x^alpha * exp(alpha * x0_val / x)\n", " q_espression_n = q_espression / q_espression.substitute(x=x0_val)\n", "\n", " r = q_espression_n / x^2\n", " r_espression_n_compiled = fast_callable(r, vars=[alpha, x], domain=RRR)\n", "\n", " r_espression_norm_int1 = lambda alpha_val: numerical_integral(partial(r_espression_n_compiled, alpha_val), a=0, b=10, algorithm=\"qag\", rule=6)[0]\n", " r_espression_norm_int2 = lambda alpha_val: numerical_integral(partial(r_espression_n_compiled, alpha_val), a=10, b=+oo, algorithm=\"qag\", rule=6)[0]\n", " r_espression_norm = lambda alpha_val: r_espression_norm_int1(alpha_val) + r_espression_norm_int2(alpha_val)\n", "\n", " def r_n(alpha_val):\n", " evaluated_norm = r_espression_norm(alpha_val)\n", " def current_r_n(x_val):\n", " return r_espression_n_compiled(alpha_val, x_val) / evaluated_norm\n", " return current_r_n\n", "\n", " def r_n_mode(alpha_val):\n", " current_r_n = r_n(- alpha_val)\n", " return x0_val - minimize_scalar(lambda t: -current_r_n(t), bounds=(0, x0_val + 0.01)).x\n", "\n", " def r_n_mean(alpha_val):\n", " current_r_n = r_n(- alpha_val)\n", " return x0_val - numerical_integral(lambda t: t * current_r_n(t), a=0, b=+oo)[0]\n", "\n", " def r_n_median(alpha_val):\n", " current_r_n = r_n(- alpha_val)\n", " return x0_val - find_root(lambda t: numerical_integral(current_r_n, a=0, b=t)[0] - 0.5, 0, 10)\n", "\n", " p1_plot = plot(r_n_median, xmin=1, xmax=XMAX, scale=SCALE, legend_label=\"median\")\n", " p2_plot = plot(r_n_mode, xmin=1, xmax=XMAX, scale=SCALE, color=\"red\", legend_label=\"mode\")\n", "\n", " final_plot = p1_plot + p2_plot\n", " final_plot.axes_labels(['$- \\\\alpha$', 'distance from $x_0$'])\n", "\n", " final_plot.show(frame=True)\n", "\n", "\n", "statistical_indicators()" ] }, { "cell_type": "code", "execution_count": 0, "id": "44f9b7", "metadata": { "collapsed": false }, "outputs": [ ], "source": [ ] } ], "metadata": { "kernelspec": { "argv": [ "sage-10.1", "--python", "-m", "sage.repl.ipython_kernel", "--matplotlib=inline", "-f", "{connection_file}" ], "display_name": "SageMath 10.1", "env": { }, "language": "sagemath", "metadata": { "cocalc": { "description": "Open-source mathematical software system", "priority": 2, "url": "https://www.sagemath.org/" } }, "name": "sage-10.1", "resource_dir": "/ext/jupyter/kernels/sage-10.1" }, "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.11.1" } }, "nbformat": 4, "nbformat_minor": 4 }