{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "This calculations corresponds to the stationary condition\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "%display latex" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# 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, "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, "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, "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, "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": 7, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# plot(1/q_num_P1.substitute(x=7.6)*q_num_P1, xmin=0.5,xmax=20,color='cyan')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a7a2091dcae24200a695f7c59d58c576", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Interactive function with 1 widget\n", " a_value: SelectionSlider(descripti…" ] }, "execution_count": 8, "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": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHWCAYAAAARl3+JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABuzUlEQVR4nO3dd3RU1dfG8e8ktNBCCURKqEovISEgvQqK0gWkg4CCCAICggVsiFLEQgBpIkpVmiBIkyaoQChKL4IgBJAaagLJvH+cH3lBipnMJHdm8nzWmiWZTObu4Ypuzjl7b5vdbrcjIiIiIh7Nx+oARERERMR5SupEREREvICSOhEREREvoKRORERExAsoqRMRERHxAkrqRERERLyAkjoRERERL6CkTkRERMQLpLI6AG8TFxfHyZMnyZQpEzabzepwRERExM3Z7XYuX75M7ty58fFJ/HqbkjoXO3nyJEFBQVaHISIiIh7m+PHj5M2bN9E/r6TOxTJlygSYG5M5c2aLoxERERF3FxUVRVBQUHwOkVhK6lzs9pZr5syZTVK3ezd8/z3UrQshIeDra3GEIiIi4o6cPbalQomk9scf8MEHUKECBARAs2Ywbhzs3w92u9XRiYiIiJdQUuci4eHhlChRgrCwsLu/8dxzcP48bNwIffvC2bPQpw8UKwb58kGnTvDNN3DmjBVhi4iIiJew2e1aLnKlqKgo/P39uXTp0oPP1F25Ahs2wOrVsGoV7NwJNhtUrAjPPGMeZcqY50RERMSrJSh3SAAldS6WqBtz+jQsWwZLlsDy5Sbpy5vXJHdPPw21a0P69EkbuIiIiANiYmK4deuW1WF4lFSpUpEmTZp7nldS56acvjExMbB+vUnwFi+GP/+EdOmgXj1zHq9RI8ia1fWBi4iIJFBMTAy7d+8mLi7O6lA8io+PDyVLlrwnsXNVUqfqV3eTJo2plK1bF8aMMQUVixfDggXm/F2qVGblrnlzk+QFBFgdsYiIpDC3bt0iLi6OAgUK4OfnZ3U4HuH69escPXqUW7du3Xe1zhWU1Lkzm80UVBQrBgMGwMmTJrmbNw969ICXXoInnoDWraFJE1BfPBERSUZ+fn6k1/Egh1y6dAm73U6GDBlc/t6qfvUkuXNDz57w009w6hR8/jlcvQodO0JgILRoAfPnw40bVkcqIiIi97FkyRLmzJnD1atXXf7eSupc5IEtTZJKjhxmtW79evjrL3j3XTh82GzLBgZC587mezoyKSIi4jZSp05NVFQU0dHRLn9vJXUu0rNnT/bs2cOWLVuS/+L58pnt2W3bYN8+0w9vwwaoUQOKFDHNj0+cSP64RERE5C5JdZ4OlNR5n6JF4e234eBBWLcOKleG9983id/TT5vzeDExVkcpIiIiLqakzlvZbFC9Onz1lTl/N2GCmWzx7LOQJ49ZzfvjD6ujFBERERdRUpcSZM4M3brBL7/A7t2mNcrMmWZqRZUq5tdavRMRES9Vs2ZN+vTpE/91gQIF+OSTTyyLJ6koqUtpSpSAkSPh77/NVmy6dNC2LQQFwZtvwvHjVkcoIiKSpLZs2cILL7xgdRgup6QupUqd2jQvXr0a9uyBli3hs8+gYEFTQfvTT6qcFRERr5QjRw6v7K+npE6geHHT8+7ECZPY7dsHdepAyZLwxRdw/brVEYqIiBeqWbMmvXr1ok+fPmTNmpXAwEAmTpzI1atX6dy5M5kyZaJw4cIsW7Ys/mf27NlDgwYNyJgxI4GBgbRv356zZ8/Gf//q1at06NCBjBkzkitXLkaPHn3Pdf+9/frxxx9TunRpMmTIQFBQEC+99BJXrlyJ//60adPIkiULy5cvp3jx4mTMmJEnn3ySyMjIpPmNSSQldUns11+hcWNznM3tZcpkplTs2gVr1phkr0cPUzn7zjvwzz9WRygiIl7mq6++IiAggM2bN9OrVy969OhBixYtqFy5Mtu2baN+/fq0b9+ea9euERkZSY0aNQgODmbr1q38+OOPnD59mpYtW8a/34ABA1izZg0LFixgxYoVrF27loiIiIfG4OPjw2effcauXbv46quv+Omnnxg4cOBdr7l27RqjRo3i66+/Zv369Rw7doz+/fsnye9JotnFpS5dumQH7JcuXbLb7Xb7qlV2e7FidjvY7dWr2+1Ll9rtcXEWB+mIQ4fs9p497XY/P7s9XTq7vXt3u33/fqujEhERC129etW+detW+9WrV+94zm6PiPjvxx0/Yq9Ro4a9atWq8V/funXLniFDBnv79u3jn4uMjLQD9l9++cX+1ltv2evVq3dXLMePH7cD9v3799svX75sT5MmjX327Nnx3z937pzdz8/P/sorr8Q/lz9/fvuYMWMe+Pnmzp1rz549e/zXX375pR2wHzp0KP658PBwe2BgYIJ+v8zvj/k9mzFjhn3MmDH2c+fOxX/v37lDYmmlzkUeNFGiTh1TcLpggZne1aABBAebgtNbt6yJ1SGFC8PYsaaA4o03zBiyYsWgaVPYuFHn7kREBDAnd0JD//uxb9/dP1emTJn4X/v6+pI9e3ZKly4d/1xgYCAAZ86cISIigjVr1pAxY8b4R7FixQA4fPgwhw8fJiYmhkqVKsX/fLZs2ShatOhDY1+zZg1PPPEEefLkIVOmTHTo0IFz587dNcorffr0FC5cOP7rXLlycebMGcd/o5KQkjoXedhECR8faNLEbMWuWQO5cpmC0yJFYNw4Dzmylj27qY796y+YONH8qaxa1TQ3XrJEyZ2ISApXrBhERPz34385WLzUqVPf9bXNZrvrOZvNBkBcXBxxcXE0bNiQHTt23PU4ePAg1atXx56I/xf99ddfNGjQgFKlSjFv3jwiIiIIDw8H4ObNmw+NMzHXS0pK6pKRzQY1a8KPP5qJXhUrQq9ekD8/DBsGFy9aHWECpEsHXbua5cfFi8HXFxo2hJAQ0yIlLs7qCEVExALp05v/FfzXw5mi05CQEHbv3k2BAgV49NFH73pkyJCBRx99lNSpU/Prr7/G/8yFCxc4cODAA99z69at3Lp1i9GjR/P4449TpEgRTp48mfggLaSkziLlysGsWXDggBny8N57UKCAqUe4dMnq6BLAxweeecbMmF2zBrJlMx+kdGkP2lsWERFP0rNnT86fP0/r1q3ZvHkzf/75JytWrOD5558nNjaWjBkz0qVLFwYMGMDq1avZtWsXnTp1wsfnwelO4cKFuXXrFp9//jl//vknX3/9NRMmTEjGT+U6SuosVriw2YI9ehQ6d4YPPzTJ3fvvQ1SU1dElwO3lx9WrzRm7AgXM3nLx4vDll3DH0rWIiIgzcufOzcaNG4mNjaV+/fqUKlWKV155BX9///jEbeTIkVSvXp1GjRpRt25dqlatSmho6APfMzg4mI8//piPPvqIUqVKMWPGDIYPH55cH8mlbHZ32xD2cFFRUfj7+3Pp0iUyZ87s8M+fPGkSu4kTIUMG6N8fXn7ZdBvxGBERZj95wQKztzxoEDz/PKRJY3VkIiLiAteuXWPv3r0UL17cK5v4JoXbv2f79+/nzJkzdOjQgWzZsgHO5w63aaXOzeTObfr/Hj4MrVvD22+bIQ8ffQR39EF0b6Ghpkr299/h8cdN77tixWD6dIiNtTo6ERERr6Skzk3lyWM6iRw6BC1awFtvQaFCMGqUh1TLgjlfN3u2Se6Cg6FjRyhVCr79VgUVIiIiLqakzs0FBcH48XDwoGkNN3gwPPYYTJ7sQbUIpUqZlbvNm810ipYtoXx5WLpUrVBERERcREmdh8if34xh3bcPqleHbt3MQtiCBR6UF4WFwfLlsG6dOTD49NOm193atVZHJiIiXqxmzZr06dPH6jCSnJI6D1O4sOkYEhFhFr2aNTP9f9etszoyB1SvDuvXw7JlEB0NtWrBE0+YDyUiIiKJoqTOQ4WEmEWvVatM15CaNc3C1++/Wx1ZAtls8OSTsGWLaVp84oTZkm3XzvR3EREREYcoqfNwdeqYo2pz5phGxrfrEU6csDqyBLLZzHLj77+b/eVVq0yl7MCBHjJiQ0REEqtmzZr06tWLPn36kDVrVgIDA5k4cSJXr16lc+fOZMqUicKFC7Ns2bL4n1m3bh0VKlQgbdq05MqVi0GDBnHrjkPmV69epUOHDmTMmJFcuXIxevToe64bExPDwIEDyZMnDxkyZKBixYqs9YKjQErqXCQ8PJwSJUoQFhaW7Nf28TG1B3v2mIrZpUvNXNl334Vr15I9nMRJlQpeeMGU+w4aBOHhZq/5k08gJsbq6EREJIl89dVXBAQEsHnzZnr16kWPHj1o0aIFlStXZtu2bdSvX5/27dtz7do1Tpw4QYMGDQgLC2Pnzp2MHz+eKVOm8P7778e/34ABA1izZg0LFixgxYoVrF27loh/He/p3LkzGzduZPbs2fz++++0aNGCJ598koMHDyb3x3cpNR92MVc1EHTGxYum9++nn0JgIAwfDm3amOTPY0RGwtChMGWKmVIxfLjp7fK/wc4iImKd+zYfvnbNVPP9l2LF4gfA1qxZk9jYWDZs2ABAbGws/v7+NGvWjOnTpwNw6tQpcuXKxS+//MLixYuZN28ee/fuxfa//x+MGzeO1157jUuXLnHt2jWyZ8/O9OnTadWqFQDnz58nb968vPDCC3zyySccPnyYxx57jL///pvcuXPHh1W3bl0qVKjABx984Krfpruo+bAkSpYsMHIk7N0LFSpA+/ZQqRJs2mR1ZA7IlcuM1fj9dzNyrFUr8yE2brQ6MhERuZ99+0zz+f96/CvxK1OmTPyvfX19yZ49O6VLl45/LjAwEIAzZ86wd+9eKlWqFJ/QAVSpUoUrV67w999/c/jwYWJiYqhUqVL897Nly0bRokXjv962bRt2u50iRYqQMWPG+Me6des4fPiwy39bklMqqwOQpFO4sKlBWLcO+vWDKlVMbvTRR6ZFikcoWRKWLIGffoIBA0wLlNatzYcICrI6OhERua1YsYR1MShW7K4vU6dOfdfXNpvtruduJ3BxcXHY7fa7EjqA2xuONpuNhGw+xsXF4evrS0REBL6+vnd9L2PGjP8dvxtTUpcC1KhhikynT4fXX4eiRU2SN3iwB82UrV3bfIivvjKBFytmzt717w9+flZHJyIi6dOb1gxJqESJEsybN++u5G7Tpk1kypSJPHnykDVrVlKnTs2vv/5Kvnz5ALhw4QIHDhygRo0aAJQrV47Y2FjOnDlDtWrVkjTe5Kbt1xTCxwc6dTIVsgMHwpgxJrmbMcODmhf7+EDnzuZDvPQSvPee2Zr97jsP+hAiIpJYL730EsePH6dXr17s27ePRYsWMXToUPr164ePjw8ZM2akS5cuDBgwgNWrV7Nr1y46deqEzx2HyosUKULbtm3p0KED8+fP58iRI2zZsoWPPvqIpUuXWvjpnKekLoXJmNFUxe7fb3Yy27UzPe7++MPqyByQObM5NLhrlxlB1qKF6e3iMU36REQkMfLkycPSpUvZvHkzZcuWpXv37nTp0oU333wz/jUjR46kevXqNGrUiLp161K1alVCQ0Pvep8vv/ySDh068Oqrr1K0aFEaNWrEb7/9RpCHH+tR9auLuUP1qyNWrYJevcxs2Z494Z13TKGFR1m2DPr2NR+ie3eTtWbPbnVUIiJe677Vr/JQqn6VJFe3LuzcCR9+CFOnmi3Zr76CuDirI3PAU0+ZVbqRI+Gbb+Cxx2DcOIiNtToyERGRZKOkTkiTxtQb7NtndjE7dTJbs9u3Wx2ZA9KkMdUfBw+aCRU9e0LFiqa4QkREJAVQUifx8uSBmTNhzRqIijKjWHv2hPPnrY7MATlzwuTJ8MsvZqWuYkXo0QMuXLA6MhERkSSlpE7uUbOmWaUbPdrsZhYtatqheNTpy8cfN6t0n35qMtWiRWHaNA/7ECIiIgmnpE7uK3Vq6NPHVMnWqwcdO5qt2QMHrI7MAalSmSqQffvMh+jcGapX97BSXxERkYRRUicP9cgjppfd8uVw7BiULm0qZKOjrY7MAblymSXHn36Cs2ehXDl49VW4fNnqyERERFxGSZ0kSL16ZoFrwAAYNgzKlIG1a62OykG1aplS32HDYMIEM5Vi7lxtyYqIiFdQUicJ5ucH778PO3aYeoRatUyl7NmzVkfmgDRp4LXXYM8eqFDBDMN9+mk4etTqyERERJyipE4cVqIErFtniky//97UIHz5pYcteOXPDwsWmA/wxx9QsiR8/DHcumV1ZCIiIomSyuoA3NGSJUt49dVXiYuL47XXXqNr165Wh+R2fHygSxdo2ND0uHv+edO0+Paupsdo2NCU+771lvkgM2bApElJPpRaRMQbXL9+3eoQPEZy/F5pTNi/3Lp1ixIlSrBmzRoyZ85MSEgIv/32W/woj//iaWPCXGX1ajOh66+/4PXXYfBgSJvW6qgctHkzdOtmZsr27WsqQjJksDoqERG3ExMTw+7du4nzqPFD1rPb7ezZs4fz588nyZgwrdT9y+bNmylZsiR58uQBoEGDBixfvpzWrVtbHJl7q1PH7GIOG2Ye334LU6aYdnEeo0IF2LoVxoyBoUPhu+9g/HgzhkxEROKlSZOGkiVLcuvWLS5dusSSJUvIkCEDaT3ub/PJKzY2lps3bybZ+3tdUrd+/XpGjhxJREQEkZGRLFiwgCZNmtz1mnHjxjFy5EgiIyMpWbIkn3zyCdWqVQPg5MmT8QkdQN68eTlx4kRyfgSPlS4dvPcetGgBXbtC5crQu7cprsiY0eroEih1ahg4EJo3N0uPDRpA69Ym0QsMtDo6ERG3kSZNGtKkScONGze4evUqN2/eVFKXANFJ2BPM65K6q1evUrZsWTp37kzz5s3v+f6cOXPo06cP48aNo0qVKnzxxRc89dRT7Nmzh3z58nG/3WibzZYcoXuNMmXMlK5PP4U334RFi2DiRHjiCasjc0DhwrBihelv168fFC8Oo0aZBsb690FEJF7atGnJnDkzUVFRxMTEWB2OR8icOXOSJMBefabOZrPds1JXsWJFQkJCGD9+fPxzxYsXp0mTJgwfPpxNmzYxcuRIFixYAMArr7xCxYoVadOmTYKumVLP1D3I4cPwwgum72+nTmb0WAKPJ7qPs2dNs+Lp06FuXZOhFixodVQiIm7j6tWrSboC5W3Spk1LhjvObOtMXSLExMQQERHBoEGD7nq+Xr16bNq0CYAKFSqwa9cuTpw4QebMmVm6dClDhgx54HtGR0ff9S9yVFRU0gTvoQoXhlWrYOpUkxctWwZjx5rdTY9Z8AoIMKW9bdqYDLV0aRg+HHr2NGXAIiIpXIYMGe5KUsQaKer/SGfPniU2NpbAf52NCgwM5NSpUwCkSpWK0aNHU6tWLcqVK8eAAQPInj37A99z+PDh+Pv7xz+CgoKS9DN4IpvNtD/ZswcqVTJn7po1g8hIqyNzUP36pjK2Y0dzWLBGDQ8bhisiIt4sRSV1t/37jJzdbr/ruUaNGnHgwAEOHTrECy+88ND3Gjx4MJcuXYp/HD9+PEli9ga5c8P8+aYy9pdfTBPjr7/2sKbFmTJBeLiZkXbqlDlAOGKEmhaLiIjlUlRSFxAQgK+vb/yq3G1nzpy5Z/UuoW4fEL3zIQ9ms8Gzz5pVu6efhg4doHFjD1y1q1HDzJF9+WXTlK9SJdPTRURExCIpKqlLkyYNoaGhrFy58q7nV65cSeXKlZ167/DwcEqUKEFYWJhT75NSZMtmCksXLDA9f0uWNF971Kpd+vSmInbTJrh2DUJDTcNiVX+JiIgFvC6pu3LlCjt27GDHjh0AHDlyhB07dnDs2DEA+vXrx+TJk5k6dSp79+6lb9++HDt2jO7duzt13Z49e7Jnzx62bNni7EdIUZo0gd27TX/f9u3N1x63alexImzbBq+9ZprylS9vmhiLiIgkI69rabJ27Vpq1ap1z/MdO3Zk2rRpgGk+PGLECCIjIylVqhRjxoyhevXqLrm+Wpok3sKFpt9vTAx8/rkpNvWYCtnbduwwg3B37oQBA+Dtt01XZhERkQdwVe7gdUmd1ZTUOefcOejVC2bNMmftJkyARx6xOioH3bwJI0eardhHHzXtUMqXtzoqERFxU67KHbxu+9UqOlPnGtmzw8yZpkr2doXsjBkedtYudWp4/XWIiDCrdI8/Dm+9pbN2IiKSpLRS52JaqXOdO1ftmjSB8eM9dNVu+HAzFLdkSbNqV7as1VGJiIgb0UqdeL3bq3bz5pkC05IlTYLnUX8NSZ0ahgwxJb5xcRAWZoop1NdORERcTEmduL1mzUyFbL16pniieXM4c8bqqBxUrhxs2WKKJ4YONX3t9uyxOioREfEiSurEIwQEmFW6776Dn3+GUqVMtaxHSZsWhg0zhwWvXIGQEFNQERtrdWQiIuIFlNS5iAolkkfz5mb8auXK0LQpdOoEly5ZHZWDKlQwfe1eftn0tqteHQ4etDoqERHxcCqUcDEVSiQPu93UHPTuDVmzwrRpcJ/2hO7v559NZnryJHz4oUn0fPR3LRGRlESFEpKi2WwmF/rjDyhUCGrXhr594fp1qyNzUNWqplFxly7wyivmgxw5YnVUIiLigZTUiUfLnx9Wr4YxY0zLk5AQD5zQlSGDGaHx009w9CiUKQNTpnhYma+IiFhNSZ14PB8f6NPHHFPLkMH0+n3nHdMizqPUqmWWHlu1gq5dTXM+jyvzFRERqyipcxEVSlivRAlTWPrmm6bXb+XKsG+f1VE5KFMmmDzZlPb+8osp8/3+e6ujEhERD6BCCRdToYR72LIF2reHv/4y9Qe9enlg/cHp09CtGyxebFbuPv7YJH0iIuJVVCgh8hBhYWY79sUXzdZs3bpw7JjVUTkoMBAWLYJJk0yTvrJlYeNGq6MSERE3paROvFb69PDJJ6aQ4tAhKF0avv7aw+oPbDazSrdzJ+TKZXravf46xMRYHZmIiLgZJXXi9WrXNvUHjRtDhw7w3HNw4YLVUTmocGFYv97MjR05EipWNLPTRERE/kdJnaQI/v4wfTrMng0rVpiuIT/9ZHVUDvL1hcGDYfNms1IXGmp6ucTFWR2ZiIi4ASV1LqLqV8/QqpVZtStSBOrUgf79ITra6qgcVK4cRETASy9Bv34eemBQRERcTdWvLqbqV88QF2cWuV5/HYoVgxkzTPcQj/PTT9CxI1y+DGPHQtu25hyeiIh4DFW/ijjBxwdefdXsZMbGQvny8OmnHriTefvAYMOGpodLq1Zw7pzVUYmIiAWU1EmKVras6WnXvbtpffLkk3DypNVROShLFlPWO2cOrFplynyXL7c6KhERSWZK6iTF8/MzrU+WL4ddu0xONH++1VElQsuW5gOUKWOy0z594MYNq6MSEZFkoqRO5H/q1YPff4caNaB5c3j+eXNUzaPkzg1Ll5q95AkTTBfmP/6wOioREUkGSupE7hAQAPPmwdSp8O23EBxsRrB6FB8f6N3b7CuDSew++8zDui6LiIijlNS5iFqaeA+bDTp3hh07IGdOqFoVhg6FmzetjsxBpUv//4HBV16BBg3g1CmroxIRkSSiliYuppYm3uXWLfjgA3j3XdPr95tv4LHHrI4qEX78ETp1MqW+U6eaalkREXELamkikgxSpYIhQ+Dnn02nkHLlYPJkD9zJfPJJc7auUiVo1Mg0Lr52zeqoRETEhZTUiSTA44+b7djWraFbN2jWzAPbweXIAYsWwfjxMG2aWXrcvt3qqERExEWU1IkkUMaMMGmSaXeyfr3pcbdmjdVROchmM2fsIiIgXTqoWBFGjfLArssiIvJvSupEHNS0Kezcac7W1aljRo15XBFF8eLw66+ml92AAaafy4kTVkclIiJOUFInkgh585rhDR98ACNHQpUqcOiQ1VE5KG1aGDHCfJC9e0217Lx5VkclIiKJpKROJJF8fWHQINi4Ec6fN0UU06d7YBFFnTqm63KtWvDss9C1K1y5YnVUIiLiICV1Ik6qUMHUGzRrBh07Qtu2cOmS1VE5KHt2+O47mDIFZs82Gert5sUiIuIRlNSJuECmTPDVVzBjBvzwg5lEsWmT1VE5yGYzs9G2b4esWaFyZbO/HBtrdWQiIpIASupcRBMlBKBNG9P6JHduqF7dNC2+dcvqqBz02GNmT/m11+DNN8227F9/WR2ViIj8B02UcDFNlBAwidz778N775kFrxkzIF8+q6NKhPXroX17s588cSK0bGl1RCIiXkcTJUTcWKpU8PbbsG4dHDsGZcrA3LlWR5UI1aub/i3160OrVqaI4upVq6MSEZH7UFInkoSqVr07J+rSxQMLS7NkMcUTU6bArFmaRCEi4qaU1Ikksds50dSpMGcOhISYgQ4e5XYRxbZt4Odn5qaNGaNJFCIibkRJnUgysNmgc2eTE2XODJUqmabFHpcTFS1qJlH07An9+sHTT8Pp01ZHJSIiKKkTSVZFiphWJ337wsCBZjrXyZNWR+WgtGnh449h2TKTpZYtC8uXWx2ViEiKp6ROJJmlSQMffQQrV8KePaaI4vvvrY4qEZ580kyiCA42v371VYiOtjoqEZEUS0mdiEXq1jU5UZUq0Lix2dG8ft3qqBwUGAhLl8Lo0fD556Z/y4EDVkclIpIiJUtSd+XKFU6dOkVMTExyXE7EYwQEwMKFMG6cKaQIC4Ndu6yOykE+PuZ83a+/mtLekBD48ksPHIIrIuLZkiypmzdvHo0aNcLf3x9/f3/y5MmDn58f2bNnp3bt2syfPz+pLi3iUWw26NEDtm41vw4Lgy++8MCc6HZZb8uWplK2dWu4eNHqqEREUowkmSgxevRolixZQtOmTcmfPz8ZMmQgTZo0XLlyhcuXL3Pw4EHmz59P27ZtefXVV119eUtpooQ44/p1czRt/Hho3hwmTTJjWD3OnDnwwgsm+JkzzbasiIjcl6tyh1QujCne6dOnWbNmzUNf88Ybb3hdQifiLD8/sxVbt65pVBwcbHKiKlWsjsxBrVpBhQrQtq2ZSvH22zB4MPj6Wh2ZiIjXSpLt11y5cv3na2w2W4JeJ5ISNWtmJlEEBUGNGmaObGys1VE5qGBBMzt28GAYMgTq1IG//7Y6KhERr5UkSd2RI0d499132bVrF5cvXybujg6r0dHR7N27l9dff53du3cnxeUtER4eTokSJQgLC7M6FPES+fLB2rXw+usmJ6pbF06csDoqB6VKBe+9B2vWwOHDpn/LggVWRyUi4pWS5EzdjRs3GDx4MJMmTeL6/3o0+P5v2yU2NhZ/f3+aN2/OqFGj8Pf3d/XlLaUzdZIU1q6Fdu3gxg1TWNqwodURJcL589C1q0nqXnzRNDBOn97qqERELOeq3CFJkrrbbty4wb59+zh9+jTnz58nU6ZM5MqVi+Dg4Pgkz9soqZOkcvasKSpdvBh69zYNjNOlszoqB9ntMHEi9OkDhQubgoqSJa2OSkTEUh6R1KVESuokKdntMHYs9O8PJUrA7NlmHKvH2b3bFFMcPgyffGIqZW02q6MSEbGEq3IHl56p27t3L7179+aXX35x5duKyP/YbNCrF2zebNqfeGyf35IlYcsW6NQJunc3ve3U005ExCkuTeqGDRvGhAkT6Nat213PT5kyha5du/K3Kt9EXKJsWdPn97nnzJZs27YQFWV1VA7y8zMN+b791gzCDQ4G/YVQRCTRXJrUZc2alR9++IFRo0bd9XyXLl348MMPGThwILs8bgaSiHvKkAGmTIFZs2DJEihXzqzgeZxnn4UdOyB3bqhWDT78EO6omBcRkYRxaVJns9koUqQITz755D3fCwgIYNKkSYwYMcKVlxRJ8Z57zuREAQGmSfGIER6YExUoAOvWwcCBpodL/fpw6pTVUYmIeBSXJnUjR46kT58+9OvXjxUrVnD16tW7vn97XJiIuFahQvDzz2bE2GuvwVNPeWBOlDo1fPABrFgBu3aZPebly62OSkTEY7i0+vXLL7+ka9eu3H7LVKlSERoaSs2aNalQoQJXrlxh+vTprFy50lWXdDuqfhWrrVwJ7dub4onp082il8c5cwY6dDBJ3YABZqSG/kIoIl7KLWe/zp07lyNHjpAuXTq2bt3KTz/9xJo1axg5ciRxcXHkyZOH+fPnu/KSIvIvTzxhRox17AhPPmnanwwb5mE5Uc6csHSpaVA8eLDZmp01yyxJiojIfbl0+/XRRx8lX7585MyZkwYNGjBq1CgiIiI4e/YsM2bMoHz58lq9EkkGgYEmJxo1Cj79FKpWNS3hPIqPj8lIN240nZfLlTPNikVE5L5cmtSlSpWKU/c5yJMlSxZat27NrFmzGD58uCsvKSIP4ONjztht2mQmdJUrBzNmWB1VIlSoANu2QYMGpiqkWze4ds3qqERE3I5Lk7qhQ4fSr18/Vq1adc/33nvvPUaMGEHGjBldeUkR+Q/ly5ucqFEjMz+2c2f4Vw2T+/P3h5kzTQ+XGTPMh/rjD6ujEhFxKy5N6rJkycJXX33F3r17mfGvJYElS5bw9ttvExMT48pLikgCZM4MX38N06aZXr/ly8Pvv1sdlYNsNtNpOSLCVMqGhZnmxR43TkNEJGkk2+zXs2fPsm7dOurXr+/Vq3WqfhV3t3+/mcq1f78Zu/riix44dvX6dXPebtw4aNYMJk+GrFmtjkpEJFHccvbrwwQEBNC8eXOvTuhEPEHRovDbb9ClC/To4aFjV/38IDwc5s+Hn34yI8Y2brQ6KhERSyVbUudJmjZtStasWXn22WetDkUkSaRLZ3Ki774zfe3KlTOJnsdp2tSM0wgKgho1TO+W2FiroxIRsYSSuvvo3bs306dPtzoMkSTXvLnJiR55xLQ9GTnSA0eM5c8Pa9eafnZvvQX16sHJk1ZHJSKS7JTU3UetWrXIlCmT1WGIJIsCBWD9etP+ZOBAePppM9DBo6RKBe+9B6tWwd69ZsTYsmVWRyUikqw8Lqlbv349DRs2JHfu3NhsNhYuXHjPa8aNG0fBggVJly4doaGhbNiwIfkDFfEgqVPDhx/Cjz+a4tLgYHNUzePUrm3GaVSoYPravfoqqOJeRFIIj0vqrl69StmyZRk7dux9vz9nzhz69OnDG2+8wfbt26lWrRpPPfUUx44di39NaGgopUqVuudxUls2ksLVr29youLFoW5dGDIEbt2yOioH5cgBixebEWOffw5VqsChQ1ZHJSKS5JKtpUlSsNlsLFiwgCZNmsQ/V7FiRUJCQhg/fnz8c8WLF6dJkyYOTbNYu3YtY8eO5bvvvnvo66Kjo4mOjo7/OioqiqCgILU0EY8WGwvDh8PQoSYnmjkT8ua1OqpEiIiAVq3MfvKECdCmjdURiYjcw+NamiSHmJgYIiIiqFev3l3P16tXj02bNiXJNYcPH46/v3/8IygoKEmuI5KcfH3hzTdN/cGRI2Y7dskSq6NKhNDQ/x+n0bat6ePiceM0REQSxquSurNnzxIbG0tgYOBdzwcGBt53Ju2D1K9fnxYtWrB06VLy5s3Lli1bHvjawYMHc+nSpfjH8ePHEx2/iLupVs1Ux1auDA0bQt++HnhE7fY4jalTYfZsc95u1y6roxIRcblUrnqj5cuX8+OPP/Lnn39y5coVHrSra7PZWL16tasu+8Br3Mlut9/z3MMsX748wa9NmzYtadOmTfDrRTxN9uywaBF89hkMGAAbNsCcOVC4sNWROcBmM0NvK1Y027FhYeYDde3qgeM0RETuz+mkLioqiiZNmrBu3boHJnJ3ciS5clRAQAC+vr73rMqdOXPmntU7VwsPDyc8PJxYNT4VL2SzwSuvmF52rVqZZsUTJ8Jzz1kdmYNKlIDNm6FPH3jhBVPi+8UXZjVPRMTDOZ3Uvfbaa6xdu5Zs2bLxwgsvUK5cOXLkyJGkyduDpEmThtDQUFauXEnTpk3jn1+5ciWNGzdO0mv37NmTnj17xh92FPFGt4+ode8OrVvD6tXw6aeQPr3VkTnAz88kcrVrQ7duEBJilh5DQ62OTETEKU4ndfPnzyd16tSsW7eOkiVLuiKmh7py5QqH7mhPcOTIEXbs2EG2bNnIly8f/fr1o3379pQvX55KlSoxceJEjh07Rvfu3ZM8NpGUIHNmmDED6tSBXr3gl19MTpQMf/xdq1UrKF/eLDdWqmTGafTure1YEfFYTrc0yZgxI4UKFeL33393VUwPtXbtWmrVqnXP8x07dmTatGmAaT48YsQIIiMjKVWqFGPGjKF69erJEp+rypJFPMHu3SY3+vNPs2LnkUfUYmJg0CAYM8ZUyU6dag4SiogkE1flDk4ndeXLl+fSpUscPHjQmbfxeHeeqTtw4ICSOkkxrl0zVbETJ5oEb+JEDz2itngxdOpk9pJnzTIHCEVEkoHb9Knr2bMnhw8fZu3atc6+lUfr2bMne/bseWj7ExFvlD69OaI2e7YZt1quHGzdanVUidCwoenfUqAA1KwJH3wAcXEWByUiknBOJ3WdO3emV69eNGvWjM8//5wrV664Ii4R8TCtWsH27ZAtm+lrN2YMeNy8mqAgWLMGBg823Zfr1wcHelyKiFjJJWPCoqOjad26NYsWLQIgR44cpH9AOZzNZuPw4cPOXtJt6UydpHQxMSYn+vhjeOYZ+PJLCAiwOqpEWLUK2rUzv/7mGzMMV0QkCbjNmbrTp09Tt25d9uzZk+A+dd7cy01JnYjxww/QsSOkS2dmxyZTrZJrnT4N7dubBO/11+HttyGVy3q2i4gArssdXNKnbvfu3Tz66KMMGDCA4OBgy/rUWUnNh0Xu9vTT5oha27ZQqxYMHQpvvGHmynqMwED48Uf46CN46y1Yt85kqJrxLCJuyOmVukceeYSoqCgOHTpE7ty5XRWXx9JKncjdbt2C994zj5o1zU6mR/6nYuNG03H56lWYNs0UVoiIuIDbVL9evXqVYsWKKaETkftKlQreecdMn9i3D4KDzeKXx6lSxSw9Vq1q+tn17WsOEIqIuAmnk7rSpUtz7tw5V8QiIl6sVi2TE4WGwlNPwcCBcPOm1VE5KFs2WLgQPvkEwsNNma8XF36JiGdxOqkbMGAAx48fZ+7cua6IR0S8WM6cpoBi5EjT8qRaNThyxOqoHGSzwSuvwKZNcPGiacw3e7bVUYmIOJ/UNW3alM8++4yuXbvy6quvsnv3bm7cuOGK2DxKeHg4JUqUICwszOpQRNyajw/07w8//2yKS8uVg+++szqqRChfHrZtMxUhrVvDCy+Y8RoiIhZxulDC18FSNpvNxq1bt5y5pFtToYRIwl28CN26maSue3fT287Pz+qoHGS3w5Qp0Ls3FCoEc+dCiRJWRyUiHsRtCiXsdrtDjziN3RGR/8mSxeRAEyaYJsWPPw7791sdlYNsNujaFTZvNgle+fIwdaoHjtMQEU/ndFIXFxfn8ENE5DabDV580eRE0dGmkOKbb6yOKhFKlYItW6BNG+jSxTQtvnzZ6qhEJAVxOqkTEXGFMmVg61Zo1szkQ88/b1rCeZT06WHyZJgxAxYtMhnq9u1WRyUiKYTLk7oDBw6wZMkSZs2axZIlSzhw4ICrLyEiXipjRpg+3WzFzpkDFSrA7t1WR5UIbdqYIoqMGc2e8tix2o4VkSTnsqTuiy++oFChQhQvXpzGjRvTrl07GjduTPHixSlcuDCTJk1y1aXckqpfRVynUyezk+njA2Fhpg7B43Kixx6DX34xe8u9ekHz5nDhgtVRiYgXc7r6FaBz585Mnz4du91O2rRpCQoKIjAwkNOnT3P8+HGio6Ox2Wx06NCBL7/80hVxuy1Vv4q4zrVrpiXc5Mlm8WvCBMiUyeqoEmHhQujcGTJnNj3tKlWyOiIRcSNuU/06c+ZMvvrqK9KnT8+IESP4559/OHDgABs2bODAgQP8888/jBgxggwZMjB9+nRmzZrl7CVFJIVInx4mTTJH1L7/3oOPqDVpYsZp5MljOi5/9BGoaExEXMzppG7SpEnYbDbmzZtH//79yZgx413fz5gxI/379+e7777Dbrd7/TasiLje7SNqGTKYI2rjxnngdmz+/LBuHQwYAIMGQYMGcOaM1VGJiBdxevs1W7ZsZM+enYMHD/7na4sUKcI///zDBS8+V6LtV5Gkc+OGyYnGjoVnnzWreFmyWB1VIixfbkp8fX3NMmTt2lZHJCIWcpvt1xs3bpAlgf9VzZw5M9HR0c5eUkRSqHTp4PPPzQSKlSshJMQUVHic+vVh504zeaJuXRgyBLx40o6IJA+nk7p8+fKxa9cuzp49+9DX/fPPP+zevZt8+fI5e0kRSeGaNzdn63LkgCpVYMwYD9yOzZULVqyAd96BYcOgTh04ccLqqETEgzmd1DVq1Ijo6GhatWrFP//8c9/XnDlzhlatWhETE0Pjxo2dvaSICAULwoYNZuRqv37QuDGcO2d1VA7y9YW33oI1a+DwYShbFpYutToqEfFQTp+pO3/+PMHBwZw4cYK0adPSokULSpQoQc6cOTlz5gx79uzh22+/5caNGwQFBbF9+3ayZcvmqvjdRnh4OOHh4cTGxnLgwAGdqRNJRkuWQMeOppBi1iyzeudxzp41Dfp++MEcHBw2DFKntjoqEUkGrjpT55I+dYcOHaJ169ZERESYN7XZ4r93++3DwsKYOXMmhQsXdvZybk2FEiLWOH4cWreGX3+F99+HgQNN82KPEhdn9pIHDYLy5U1Pu/z5rY5KRJKYWyV1t61evZoVK1Zw4MABrly5QsaMGSlSpAj169endgqp7lJSJ2KdW7dg6FAYPhyeeAK+/hpy5rQ6qkT49Vd47jm4dMnMTGvSxOqIRCQJuU1Sd+zYMQDy5s2Lj8f9tdj1lNSJWG/FCmjXDlKlgpkzoWZNqyNKhAsX4PnnzTSKV14xDYvTprU6KhFJAm7T0qRAgQJUrFjR2bcREXGZevVMx5BixUxR6TvvQGys1VE5KGtWmD8fPvsMxo83BwUPH7Y6KhFxY04ndf7+/uTPn1+rdCLiVnLlMr3shg6Fd98127EnT1odlYNsNujVCzZtMit3ISHw7bdWRyUibsrpTKx06dLxW7AiIu7E19f09V29Gvbtg+BgszXrcUJDzZy0J5+Eli3hpZfMeA0RkTs4ndS98sornDp1iqlTp7oiHhERl6tZE3bsMAtd9evD66974AAHf39TDTthAkydaobgHjhgdVQi4kacTuqaN2/Ohx9+SM+ePenbty/btm3j+vXrrohNRMRlcuY0fX0//BBGjDCJ3vHjVkflIJsNXnwRfvsNrl83WeqMGVZHJSJuwunqV19fX8cuaLNxy+P+ivzf1HxYxHNs2mQ6hly9Cl99Bc88Y3VEiXDlCvToAd98A126mIKK9OmtjkpEEsFtWpokpkAiLi7OmUu6NbU0EfEM585B586weLEZMzZ8OKRJY3VUDrLbYdo06NkTChWCuXOhRAmroxIRB7lNS5O4uDiHHyIiVsueHRYtMgMcPv8cqlWDI0esjspBNpvJTLdsMQle+fKmWbHresqLiAdxKKlbv349O3fuTKpYRESSlc0GffrAxo3wzz9QrhzMm2d1VIlQsqRJ7Fq3Ng2LO3Y027MikqI4lNTVrFmT3r173/Vc7dq16dOnjytjEhFJVmFhsH276WX37LPw8sse2DEkfXqYMsXMRps/36za/f671VGJSDJyKKmz2Wz3bJ+uXbuWbdu2uTQoEZHk5u9vjqSNGweTJ0PlynDwoNVRJUK7dhARYUaKVagAX3yh7ViRFMKhpC5Lliz89ddfSRWLiIilbDZTUPrrr2b3MiQEZs2yOqpEKFrUfIjOnaF7d7MtGxVldVQiksRSOfLiSpUqsWzZMpo3b069evXw8/MD4MyZM0yfPj3B79OhQwfHohQRSUbBwWaxq3t3aNMGfvoJPv3UwzqG+PmZmbG1akHXriZDnTvX/FNEvJJDLU1+//13atasycWLF7HZbADY7fb4XydUrMdN1k44tTQR8R52uykmffllKFwY5szx0I4hhw9Dq1bwxx8wapT5QA7+d1tEko6rcgeHVurKlCnD/v37mT17Nvv27eP69etMmzaNnDlz8uSTTyY6CBERd2SzmWLSihXNyNWwMAgPh06drI7MQYULmxLfgQOhd29Ys8YUVWTNanVkIuJCLmk+XLVqVdavX++qmDyaVupEvNPVqyYfmjoV2rc3BRUZM1odVSIsXGjO2mXJYmbJVqxodUQiKZ7bNB8eOnQonTt3dvZtRETcWoYMXtIxpEkT078lMBCqVoXRo1UdK+IlnF6pk7tppU7E++3fb46o7dtnCiheeMEDj6jdvAmvv27O2D3zjBk3lj271VGJpEhus1InIpLS3O4Y8vzzpkL2uec8sGNI6tQwciQsWQK//GJKfjdutDoqEXGCkjoXCQ8Pp0SJEoSFhVkdiogkg3TpzLm6uXPhxx9Np5CICKujSoSnn4YdO6BAAahRA4YPB83oFvFI2n51MW2/iqQ8f/5ptmN37jS7mb16eeB27K1bMHSoSeqeeMIcHsyZ0+qoRFIEbb+KiLiJQoXg55+hZ0945RVo1gzOn7c6KgelSgXDhpllx+3bzXbs2rVWRyUiDlBSJyLiAmnTwpgxsGgRrFsH5cqZc3cep149s+RYrBjUqQPvvANe3DBexJsoqRMRcaFGjcwRtTx5oFo1U4vgcUfUcuWClSthyBCT1NWrB5GRVkclIv9BSZ2IiIvly2dW61591QxxaNgQzp61OioH+fqaM3arV8OePWY7duVKq6MSkYdwaVJ3/PhxZs6cyciRI3n33Xfv+t7NmzeJiYlx5eVERNxW6tTw4YewbBls3gxly4JHDt6pVctsxwYHQ/368OabpqhCRNyOS5K6s2fP0qpVKwoWLEj79u0ZNGgQ77zzzl2v6dy5M35+fkR4ZM2/iEjiPPmkyYkee8zkR++/74FH1HLmNNnpsGEmU61dG/7+2+qoRORfnE7qLl++TI0aNfj222/JkycPnTp1Ik+ePPe8rmvXrtjtdubPn+/sJUVEPEru3LBqlVnkGjLELHidOmV1VA7y8YHBg01F7J9/mpW7pUutjkpE7uB0UjdixAj27t1L8+bN2bdvH1OmTCF//vz3vK569er4+fmxZs0aZy8pIuJxUqUyNQerVsHu3SYnWrXK6qgSoWpVUwny+OOmcfHAgWbkmIhYzumk7rvvviNt2rRMnjwZPz+/B1/Ix4dHH32UY8eOOXtJERGPVbu2yYnKlDFFpUOGeOARtYAA+P5702l5zBioXh3++svqqERSPKeTuqNHj1KkSBH8/f3/87Xp06fnrMeVgImIuFZgoOnx+9575phanTpw8qTVUTnIx8eU927YYNqdBAfDwoVWRyWSojmd1KVLl47Lly8n6LWRkZEJSv5ERLydjw+88QasWQOHDpnq2OXLrY4qER5/3EygqFkTmjaFPn0gOtrqqERSJKeTupIlS3L8+HH++o+l9x07dnDs2DFCQ0OdvaSIiNeoXt1sx5YvbyplBw/2wO3YrFlh/nz47DMYPx6qVIHDh62OSiTFcTqpa9euHbGxsbzwwgtcu3btvq+5cOECXbp0wWaz0aFDB2cvKSLiVXLkgB9+gI8+MhMoataE48etjspBNhv06gWbNsGFCxASAt9+a3VUIimK00ldt27dqFatGitXrqR06dIMGjSI06dPAzB16lT69etH0aJF2b59O0888QTPPfec00GLiHgbHx9TSLp+PRw7Zo6o/fCD1VElQmgobNtmlh1btoSXXoIbN6yOSiRFsNntdruzb3L58mVeeOEF5syZg81m4/Zb3vnrli1bMmXKFDJkyODs5dxaVFQU/v7+XLp0icyZM1sdjoh4oHPnoHNnWLwY+veHDz4wEyo8it0OEyfCK69AsWIwdy4UKWJ1VCJuyVW5g0uSutv++OMPFixYwB9//MGlS5fImDEjJUqUoGnTph5zlu748eO0b9+eM2fOkCpVKt566y1atGiR4J9XUicirmC3m24hr71mFr/mzIH7tAB1fzt3mhW7Eyfgiy+gbVurIxJxO26Z1HmDyMhITp8+TXBwMGfOnCEkJIT9+/cneIVRSZ2IuNJvv0GrVnDpEkybBo0bWx1RIly5Aj16wDffQJcupqAifXqroxJxG67KHVwy+9Wb5MqVi+DgYABy5sxJtmzZOH/+vLVBiUiKVbHi/3cMadLEdAyJibE4KEdlzAjTp8PUqTBzJlSoAHv2WB2ViNdxOqn7/vvvKVSoEKNHj37o60aPHk2hQoVY6uSswPXr19OwYUNy586NzWZj4X2aXY4bN46CBQuSLl06QkND2bBhQ6KutXXrVuLi4ggKCnIqZhERZ9yvY8iff1odlYNsNnNQcMsWs7dcvjx8+aX5tYi4hNNJ3fTp0/nrr79o2rTpQ1/XuHFjjh49yvTp05263tWrVylbtixjx4697/fnzJlDnz59eOONN9i+fTvVqlXjqaeeums8WWhoKKVKlbrncfKOlu7nzp2jQ4cOTJw40al4RURc4c6OIefPQ7lyMG+e1VElQsmSJrFr3Rqefx46djTbsyLiNKfP1BUuXJhr164RGRn5n6/NlSsXGTJk4NChQ85cMp7NZmPBggU0adIk/rmKFSsSEhLC+PHj458rXrw4TZo0Yfjw4Ql63+joaJ544gm6detG+/bt//O10Xd0T4+KiiIoKEhn6kQkyVy6BN26mTZwPXuaEazp0lkdVSJ88w107w5585rq2DJlrI5IxBJuc6bu5MmT5MuXL0GvDQoKSlDyl1gxMTFERERQr169u56vV68emzZtStB72O12OnXqRO3atf8zoQMYPnw4/v7+8Q9t1YpIUvP3N9Ww48bB5MlQuTIcPGh1VInQrh1EREDatOac3RdfaDtWxAlOJ3UZMmTgn3/+SdBrz549S9q0aZ295EPfPzY2lsDAwLueDwwM5NSpUwl6j40bNzJnzhwWLlxIcHAwwcHB/PHHHw98/eDBg7l06VL847jHtYEXEU9ks5mC0l9/NbuXoaEwe7bVUSVC0aLmQ3TubFbtWreGqCiroxLxSKmcfYPSpUuzfv16tm7dSvny5R/4uq1bt3L06FGqVq3q7CX/k81mu+tru91+z3MPUrVqVeLi4hJ8rbRp0yZpoioi8jDBwWax68UXTT60Zg188gn4+VkdmQP8/EwFSK1a0LWrGTE2d675p4gkmNMrdW3atMFut9O2bVv+fEA51pEjR2jbti02m402bdo4e8kHCggIwNfX955VuTNnztyzeici4i0yZYIZM2DSJNM55PHHYf9+q6NKhJYtTf+WLFmgUiX4/HNtx4o4wOmk7vnnn6dy5cocPHiQUqVK0a5dOz7//HO+/vprPv/8c9q2bUupUqU4ePAglSpVolu3bq6I+77SpElDaGgoK1euvOv5lStXUrly5SS7LkB4eDglSpQgLCwsSa8jInI/NptZ5Nq82fSxCw01dQgep3Bh2LjRbMX27g3Nm8OFC1ZHJeIRXDJR4uLFi3Tu3JlFixaZN71jq/P22zdt2pQpU6aQJUsWp6515cqV+OrZcuXK8fHHH1OrVi2yZctGvnz5mDNnDu3bt2fChAlUqlSJiRMnMmnSJHbv3k3+ZJixo4kSImK1K1dMVez06aZryOefe+gAh4ULzVm7LFnMgcGKFa2OSCRJuOWYsK1bt7Jo0SL27t1LVFQUmTJlomTJkjRp0oQQF52NWLt2LbVq1brn+Y4dOzJt2jTANB8eMWIEkZGRlCpVijFjxlC9enWXXP+/KKkTEXcxbZpJ7goUMO1PSpSwOqJEOHoUnnvOHBz88EPo188sS4p4EbdM6kRJnYi4lz17zFG1P/80LVA6dbI6okS4eRNef9005HvmGZOtZs9udVQiLuM2ferE0Jk6EXFHJUqYc3atW5udTI8c4JA6NYwcCYsXm5Ea5cqZf4rIXVy+UnfhwgWuXLnCw942oc2KPZFW6kTEXd0e4BAUZDqGlC5tdUSJcPy4yVB//RWGDYMBA8BH6xPi2dxqpe7AgQO0adOGbNmyERAQQIECBShYsOB9H4UKFXLFJUVExEHt2sHWrZAmjRngMGmSB3YMCQoyzfgGDIBBg8x2bAIb4It4O6ebD+/YsYMaNWrEr86lS5eOHDly4KO/OYmIuJ1ixcwiV9++8MILJj/64gvT685jpE4Nw4dDjRrQvr3pwDxrFiRTQZyIu3I683r99de5fPkytWvX5vfff+fatWv89ddfHDly5IEPb6QzdSLiKfz8YMIEkwctWWIGN2zfbnVUifDkk7BjBzz6qJlGMWwYODARSMTbOH2mLkuWLMTFxREZGUmGDBlcFZfH0pk6EfEkBw9Cq1amSnbMGHPmzuM6hty6Be+8Y5K6unXh669BU4TEg7jNmbq4uDiKFi2qhE5ExAM99pgpJO3aFV56ySR4ly5ZHZWDUqWC996DFStg506zHbtmjdVRiSQ7p5O64OBgIiMjXRGLiIhYIF06GDvWNChevtxsx27danVUiVC3rtmOLV7c/PqddyA21uqoRJKN00nd4MGDiYyM5Ouvv3ZFPCIiYpFnnzVn67Jlg8qV4bPPPLA6NlcuWLkShgwxSV29eqCFB0khnE7qnnrqKcaNG8dLL71E37592bVrF9evX3dFbB5FhRIi4g0KFYKNG814sVdegebN4cIFq6NykK8vDB0Kq1ebw4LBwSbRE/FyThdK+Pr6OnZBm41bt245c0m3pkIJEfEWixaZsWJZssDs2VCxotURJcLp06btyapVZtTY22+bM3gibsRtCiXsdrtDjziVm4uIeITGjc0RtUcegapV4eOPPXA7NjAQfvwR3n/f9LarUwdOnLA6KpEk4ZLqV0cfIiLiGfLnh/XroU8fePVVk+idP291VA7y8TGrdGvWwKFDZjv2xx+tjkrE5TT2QUREHip1ahg50jQq3rjR5ESbNlkdVSJUr26WHsuXh6eegsGDTY87ES+hpE5ERBLk6adNTpQvn8mPPvrIAwc45MgBP/wAH35oMtWaNeH4caujEnEJJXUuoupXEUkJgoJg7VoYOBAGDYJnnoF//rE6Kgf5+MBrr8G6dfDXX2bpcckSq6MScZrT1a+3Xb16lcWLF7Nz507Onz/PzZs3739Bm40pU6a44pJuSdWvIpJSLF8O7dpBmjRmjmz16lZHlAjnzpkS3yVLzKHB4cPNfrNIMnJV7uCSpG727Nn06NGDqKio+Oduv63tjiGCdrsdm81GrBd3+FZSJyIpycmT0KYNbNgA775rjqn5eNoekN1uBt++9hqEhsKcOaZCRCSZuE1Lk19++YX27dsTGxvLG2+8waOPPgrApEmTGDJkCI0aNcJms5EuXTqGDRvG1KlTnb2kiIi4idy5TQu4N96At96CJ580reE8is0G/frBzz/DqVNmO3bRIqujEnGY0yt1zZs3Z+HChSxcuJCGDRtSrVo1Nm3adNdq3L59+2jRogUXLlwgIiKCwMBApwN3V1qpE5GUatUqsx1rs8GMGVC7ttURJcKFC/D887BwoRmpMWKE2V8WSUJutVIXEBBAw4YNH/iaYsWKMW/ePCIjIxk6dKizlxQRETdUt66pji1Rwvz67bfB407bZM0K8+fDp5/CuHFQpQr8+afVUYkkiNNJ3blz58iXL1/812n+9zeaq1ev3vW6IkWKULJkSZYtW+bsJUVExE098gisWGESuvfegyeegMhIq6NykM0GvXubZnznz0O5cjBvntVRifwnp5O67Nmzc/369fivAwICADh8+PA9r42NjeW0xx22EBERR/j6wpAhsHo17NtnjqitXGl1VIlQvjxs2wb16sGzz8LLL8ONG1ZHJfJATid1BQoUIPKOv4aFhIRgt9uZMWPGXa/buXMnBw4cIEeOHM5e0i2pT52IyN1q1jTbscHBUL8+vPmmBw5w8PeHuXMhPBwmTYLKlc2oMRE35HRS98QTT3Dx4kV2794NQJs2bUiXLh2jRo2iXbt2hIeHM2TIEOrUqUNcXBzNmzd3Omh31LNnT/bs2cOWLVusDkVExG3kzAnLlsGwYWaIQ+3acOKE1VE5yGaDl16CX3+Fy5chJMS0PRFxM05Xv+7evZs+ffrQo0cPmjVrBsBXX33FCy+8wM2bN+P71Nntdh5//HFWrFhBxowZnY/cTan6VUTk/n7+GZ57DqKjYfp0M37V40RFwYsvwuzZ5p9jxoCfn9VRiYdzq+bD9/Pnn38yd+5cjh49ip+fH1WrVqVJkyb4+vomxeXchpI6EZEHO3sWOnaEpUvNqLH33/fAAQ52u9mK7d0bihY127NFi1odlXgwt0/qUioldSIiDxcXB6NHw+uvQ4UKZsTYHU0UPMfvv0PLlvD33/DFF9C2rdURiYdymz5169evZ+fOnQl67e+//8769eudvaSIiHgwHx8YMADWrzf5ULlysHix1VElQpkysHUrNG1qui537QrXrlkdlaRgTid1NWvWpHfv3gl67SuvvEJtj2wxLiIirlapEmzfDlWrQqNG8OqrEBNjdVQOypjRHBCcMgVmzoSKFWHvXqujkhTKJWOXHdnB1W6viIjcli2bmcg1Zgx8/jlUqwZHj1odlYNsNjNabMsWM0KjfHn46iuro5IUyCVJXUKdO3cOP1UJiYjIHWw26NMHNm6EM2fMduyCBVZHlQglS5rErmVL6NTJPP41XUkkKaVy9AeioqK4ePHiXc9FR0dz/PjxB67CXb9+nXXr1rFr1y7Kli2bqEBFRMS7hYWZ7dguXaBZM+jVC0aOhLRprY7MARkywJdfQq1a0KMH/PYbfPstlCpldWSSAjic1I0ZM4Z33333rue2bt1KgQIFEvTzXbp0cfSSHiE8PJzw8HBiPW56tYiI+8iSBb77zgxwePVVM351zhwoXNjqyBzUoYMp7W3RwmSrY8eaLdr/9W4VSQoOtzT59NNP+eSTT+K/PnbsGGnSpOGRRx65/wVsNvz8/ChUqBCtWrWiXbt2TgXs7tTSRETENbZtMzuZ//wDkyeb/MjjXL8Or7xi+tq1bQvjx0OmTFZHJW7GbfrU+fj4ULVqVbUq+R8ldSIirhMVBd26mf6+PXrAxx9DunRWR5UIM2eaCRS5c5sPo6NIcge36VP35Zdf8vrrrzv7NiIiIvfInNlM5JowAaZONW1QDh60OqpEaNMGIiLMSLGKFU2zYnWDEBdL0okSUVFRLFu2jJMnTxISEkKNGjWS6lJuQyt1IiJJY+dOsx178iRMnAitW1sdUSLcuAF9+5ostWVLsy2r/1ekeG6zUjdnzhxCQkKYPHnyXc/v27ePUqVK0aZNG/r370/t2rXp1KmTs5cTEZEUqmxZM8ChcWOz8NWtmzmy5lHSpTPn6ubMgWXLICTEHB4UcQGXJHU7d+6kevXqdz3fp08f/v77bwoVKkTjxo3JmDEjX3/9NUuXLnX2kiIikkJlygRff20KJ2bMMAWm+/ZZHVUitGxp+rf4+5s95bFjtR0rTnM6qdu5cyfZsmWjSJEi8c9FRkaycuVK8uXLxx9//MH8+fNZvHgxdrud8PBwZy8pIiIpmM1metlt3mwGOISGmkldHqdwYdOz5cUXTVO+Fi3gX31gRRzhdFL3zz//kC9fvrueW7NmDXa7nTZt2pDuf2VK1atXJ3/+/OzVTDwREXGBUqXMAIcWLaBjR+jc2QMHOKRNC599BvPmwapVZjt2yxaroxIP5XRSFxMTc0/D3Q0bNmCz2ahVq9ZdzwcGBhIZGensJUVERAAzwGHaNPOYO9dsx+7ebXVUidCsmdmODQiAKlXgk0+0HSsOczqpy5MnD4cPH+batWvxz/3444+kSpWKKlWq3PXay5cv4+/v7+wlRURE7tKxo1ng8vExAxymTvXAnKhgQfj5Z3j5ZVMh27QpnD9vdVTiQZxO6urWrcu1a9fo1asXu3bt4u233+avv/6idu3apE+fPv51169f5+DBgwQFBTl7SRERkXuUKGFGrbZta87ctW8PV65YHZWD0qQxHZYXLYL166FcOfj1V6ujEg/hdFL3xhtvkC1bNqZNm0bZsmV59913SZ06Ne+8885dr1u8eDG3bt2iWrVqzl5SRETkvtKnN63fZswweVFoKPz+u9VRJUKjRmY7NnduqFYNRo2CuDiroxI353RSly9fPrZu3cpLL71EvXr16Nq1K5s3b6ZChQp3vW7t2rWULVuWxo0bO3tJtxQeHk6JEiUICwuzOhQRkRTvzgEOFSp46ACH/PnNal3fvjBggEn0zp2zOipxY0k6USIl0kQJERH3cecAh1atzCQKj/xP8w8/mIODfn5mbtq/zqyLZ3ObiRIiIiLu6s4BDkuXmu3Y7dutjioRnn4aduwwq3c1asCHH2o7Vu7h0ErdsWPHAEidOjW5cuW66zlH/LuvnTfRSp2IiHs6dMis1u3aZWoRXnrJNDL2KLduwZAhMHw4PPmk6bqcI4fVUYmTXJU7OJTU+fj4YLPZKFasGLv/1wjo9nMJvqDNxq1btxyP1EMoqRMRcV/R0eZ42uefQ/PmZtxYlixWR5UIy5eb8t7UqWHWLPjXqE7xLK7KHVI58uJ8+fJhs9niV+nufE5ERMTd3R7gULMmPP+8GeAwZ47pbedR6tc327Ft2kCtWvDOOzB4MPj6Wh2ZWEiFEi6mlToREc9w5IjZjt2xA0aMgFde8dDt2Hffhfffhzp14JtvIDDQ6qjEQSqUEBERcYJXDHBIlcokdStWwB9/QHAw/PST1VGJRZTUiYhIiuU1Axzq1jVLjiVKmF8PHQr/mssu3i9R1a/OUvWriIi4m7/+gtatzQzZ4cOhXz8zS9ajxMbCBx/A22+b4omZM+GOc/DiniytfnWGql9FRMRd3bwJb75pztg9/TRMmwYBAVZHlQhr15oiithYc87uiSesjkgewpIzdfny5Xvgw9fXF7vdjt1ux9fXl8DAQFKlShX/XKpUqciXLx9BQUGJDlZERCQppU4NH31kBjj8+qvZjv35Z6ujSoSaNc12bHCwqZR9801TVCFezaGk7ujRoxw5cuSex9NPP43NZqN3797s27eP6OhoTp48yY0bN9i/fz+9e/fGZrPxzDPPcOTIkaT6LCIiIi7RoIHJiQoUMPmRRw5wyJkTli0zlbHDh0Pt2vD331ZHJUnI6ZYm48aNo1evXsyaNYuWLVs+8HVz586ldevWjB07lh49ejhzSbem7VcREe9x65apORg+HOrVMwMccua0OqpE2LDBHBi8cQO+/hqeesrqiOQOlpypu5+yZcsSFRWVoBW4ggUL4u/vz44dO5y5pFtTUici4n1WrIB27UwHkVmzzPhVj3P2LHTsaIbgDhxoVvBSp7Y6KsGN+tQdOnSIHAmcO5cjRw4OHjzo7CVFRESSVb16Zju2aFGzi/neex7YMSQgABYvNlUgo0ebfWUXdbUQ9+B0UpcxY0Z2797NxYsXH/q6ixcvsnv3bjJkyODsJUVERJJd7tywahW89ZbZkq1fH06dsjoqB/n4mOG3GzbA8eOmEmTxYqujEhdxOql74oknuH79Om3btuX8A1pxX7hwgbZt23Ljxg3q16/v7CWT1OXLlwkLCyM4OJjSpUszadIkq0MSERE34etrWsCtWgW7dpni0tWrrY4qESpVMkuPVapAo0bw6qsQE2N1VOIkp8/UHTt2jJCQEC5cuICfnx8tWrSgePHi5MiRg3/++Yd9+/bx7bffcvXqVbJnz87WrVvJnz+/q+J3udjYWKKjo0mfPj3Xrl2jVKlSbNmyhezZsyfo53WmTkQkZTh1ypyz++kns3o3ZIhJ+jyK3Q6ffGLO2IWGwuzZpuRXkpXbFEoA7N27l3bt2rF9+3bzpnc0KL799uXKlePrr7+mRIkSzl4u2Zw/f55y5coRERFBQAK7TyqpExFJOWJjTWXs0KFmgMOMGWab1uNs3gytWsHFi/Dll9CkidURpShuUygBULx4cSIiIli1ahUDBgygUaNG1K5dm0aNGjFgwABWrlxJRESESxK69evX07BhQ3Lnzo3NZmPhwoX3vGbcuHEULFiQdOnSERoayoYNGxy6xsWLFylbtix58+Zl4MCBCU7oREQkZfH1NX19f/oJDhww27HLl1sdVSJUqADbt0OtWtC0KbzyCkRHWx2VOMglK3XJadmyZWzcuJGQkBCaN2/OggULaHLH3yjmzJlD+/btGTduHFWqVOGLL75g8uTJ7NmzJ37mbGhoKNH3+Zd1xYoV5L7jr1inT5+mWbNmzJ8/n8DAwATFp5U6EZGU6Z9/oH17k9QNHgzvvmtaoHgUux3GjoX+/aF0aZg7FwoVsjoqr+dW269Wsdls9yR1FStWJCQkhPHjx8c/V7x4cZo0acLw4cMdvkaPHj2oXbs2LVq0uO/3o6Oj70oQo6KiCAoKUlInIpICxcWZjiFvvmlqEWbNgrx5rY4qESIioGVL09tuyhR49lmrI/JqbrX96i5iYmKIiIigXr16dz1fr149Nm3alKD3OH36NFFRUYD5TV6/fj1FixZ94OuHDx+Ov79//EOzbUVEUi4fHxg0CNatg6NHzXbs0qVWR5UIoaGwbZvp29KiBfTsaaZRiFvzqqTu7NmzxMbG3rNVGhgYyKkENhP6+++/qV69OmXLlqVq1aq8/PLLlClT5oGvHzx4MJcuXYp/HD9+3KnPICIinq9KFdMx5PHH4emnTXHpzZtWR+Ugf3+YMwfGjTOrdZUqgQYIuDVP2+1PkDurb8FU4P77uQcJDQ11aIxZ2rRpSZs2rSPhiYhICpA9O3z/PYwZY1bvNmwwHUPcuKvXvWw26NHDJHQtW0JICEyaBM89Z3Vkch9etVIXEBCAr6/vPatyZ86cSXChg4iIiKv4+Ji+vhs2QGSkGeCwaJHVUSVCcLA5Z9ewIbRuDS++CNevWx2V/ItXJXVp0qQhNDSUlStX3vX8ypUrqVy5cpJeOzw8nBIlShAWFpak1xEREc/z+OOmY0j16qYFXN++HjjAIVMm04hv4kSYPt18qP37rY5K7uBxSd2VK1fYsWNH/BbpkSNH2LFjB8f+N5S4X79+TJ48malTp7J371769u3LsWPH6N69e5LG1bNnT/bs2cOWLVuS9DoiIuKZsmaFBQvMAIfwcKhaFY4csToqB9ls0K0b/Pab6WMXGgrffGN1VPI/HtfSZO3atdSqVeue5zt27Mi0adMA03x4xIgRREZGUqpUKcaMGUP16tWTJT71qRMRkf+yZYsZ4HD+PEydCs2aWR1RIly5Ai+9BF9/Dc8/D59/DunTWx2VR1KfOjelpE5ERBLi4kXo2hXmzYOXX4ZRo8Dj6u7sdpg2zbQ8KVgQvv0WPGgcqLtQnzo3ozN1IiLiiCxZTA4UHm6OqVWuDIcOWR2Vg2w26NwZtm41X5cvb5I8sYRW6lxMK3UiIuKo7dtNx5DTp03HkFatrI4oEa5dg169zH5yhw4mW82Y0eqoPIJW6kRERLxEuXKmY8jTT5sWcN27e2DHkPTpTZPi6dPNnnJYGPzxh9VRpShK6kRERNxA5swwcyZ88QV89ZUHdwxp395sx6ZODRUqwOTJ5uydJDkldS6iM3UiIuIsmw1eeOHujiEzZlgdVSIUK2Y+RIcOpgVK27Zw+bLVUXk9nalzMZ2pExERV7hyxUzo+uYb6NIFPvvMQzuGzJ5tErtcuWDuXDOdQu6iM3UiIiJeLGNGczxtyhSzLVuhAuzZY3VUifDcc7BtG2TIYPaUx4/XdmwSUVInIiLipmw209d3yxaTB4WFeWjHkMceg19+MUuOL71kEr1Ll6yOyusoqRMREXFzJUvC5s2m1UnnztCxI1y9anVUDkqXzrQ5mTsXfvzRHBiMiLA6Kq+ipE5ERMQDZMhgWsBNnw7ffWf6/Hpkx5AWLcx2bJYspuPy559rO9ZFlNS5iKpfRUQkObRvbxa4PLpjSOHCsHGjacjXuzc8+6yZmyZOUfWri6n6VUREksP16/DKK2YCRZs2MGECZMpkdVSJsGCBOTiYJQvMmWMy1RRG1a8iIiIpmJ+fmRk7cyZ8/705orZjh9VRJULTpmZOWs6cULUqjBnjgUuP7kFJnYiIiAdr3dpsx97uGDJhggfmRAUKwIYNZnZsv37QpAmcP291VB5HSZ2IiIiHK1Lk/zuG9OjhoR1D0qSB0aPNsuOGDWYg7i+/WB2VR1FSJyIi4gW8pmNIw4ZmHzlPHqheHUaOhLg4q6PyCErqXETVryIi4g7+3TFk7FgP3I7Nlw/WrTNbsQMHmkTv7Fmro3J7qn51MVW/ioiIO4iOhgEDTBu4Zs3MuLEsWayOKhGWLoUOHcxS5OzZppjCy6j6VURERB4obVr47DOYNw9Wr4aQEDNuzOM0aGC2YwsWhJo1Yfhwbcc+gJI6ERERL9asmekYEhAAVarAJ5944HZs3rywZg289hq88YZJ9M6csToqt6OkTkRExMsVLAg//2w6hvTt66EdQ1KlgmHDTBXItm0QHGzO3Uk8JXUiIiIpgNd0DKlXD3buhKJFoXZtePddiI21Oiq3oKROREQkBfGKjiG5csGqVfDWW/D221C/Ppw6ZXVUllNS5yJqaSIiIp7i3x1DGjXywI4hvr4moVu1CnbtMtuxq1dbHZWl1NLExdTSREREPMntjiF+fjBrlod2DDl1Ctq1g59+Mqt3Q4aYpM9DqKWJiIiIOO12x5ACBUzHkA8/9MDt2EcegeXL4Z134P33oW5dOHnS6qiSnZI6ERGRFO52x5CBA2HwYHj6afjnH6ujcpCvr1ml++knOHDAbMeuWGF1VMlKSZ2IiIiQKhV88IHpGBIRYXKi9eutjioRatQwS48hIaaA4vXX4dYtq6NKFkrqREREJF79+iYnKlIEatUyu5ke1zEkRw5zWHD4cBgxwnyQv/+2Oqokp6RORERE7pI7tykqffNNU3NQvz6cPm11VA7y8YFBg2DtWjh61Cw9Ll1qcVBJS0mdiIiI3MPX19QdrFxpOoaULWuOq3mcqlXNnLTHHzeHBQcOhJs3rY4qSSipExERkQeqU8dsx5YqZYpKhw71wO3YgAAzSmPkSBgzxpy7O3bM6qhcTkmdiIiIPNT9OoZERlodlYN8fKB/fzMj7cQJsx37/fdWR+VSSupcRBMlRETEm93uGLJ6Nezfb3KilSutjioRHn/cbMdWqwaNG5uxGjExVkflEpoo4WKaKCEiIt7uzBlo394kda+/bqZ1pUpldVQOstvh00/NGbty5WD2bChY0JJQNFFCRERELJEzJyxbBsOGmQkUtWubHU2PYrNBnz6wcaPJUsuVgwULrI7KKUrqRERExGE+Pmb6xNq18OefZjt22TKro0qEsDCzHVunDjRrBr17Q3S01VElipI6ERERSbSqVU11bIUKZo7sa695YMeQLFngu+/g88/hiy/MEFwPnB2rpE5EREScEhAAixeb4Q2jR5ucyOM6hths8PLLZjba8eMQGmq2Zj2IkjoRERFxmo8PDBhgOoYcP26OqC1ebHVUiVCxohl++9hjZrzYhAmmqMIDKKkTERERl6lUyWzHVqkCjRrBq696YMeQwEDTu6V7d+jRA7p1gxs3rI7qPympExEREZfKlg0WLYKPPzbH1KpXN+NXPUrq1PDZZzBtGnzzjZlC8fffVkf1UErqRERExOVsNujbF37+GU6fNtuxCxdaHVUidOxoztZFRppzdkuWWB3RAympExERkSRToYLpGFKrFjRtalrDedx2bGioOWcXEgING8Izz8ChQ1ZHdQ8ldSIiIpKksmSBefPMbub48ea83Z9/Wh2Vg3LkgKVLYf58+OMPKFMGpkxxqyIKJXUiIiKS5Gw26NULNm2C8+fNdux331kdlYNsNrPcuHcvtG0LXbuaeWmXL1sdGaCkTkRERJJRaChs2wb160OLFtCzp0cUlt4tfXqYNAlmzDAVIeXLw86dVkelpM5VwsPDKVGiBGFhYVaHIiIi4tb8/WHOHBg3zuxgVq4MBw9aHVUitGljztr5+Zn+dl98Yel2rM1ud6PNYC8QFRWFv78/ly5dInPmzFaHIyIi4tZ27ICWLeHUKZg4EZ57zuqIEuHGDejXzxwYbNnSrOI5kAO4KnfQSp2IiIhYJjjYLHY98wy0bm36/V6/bnVUDkqXziw7zpkDP/5oqmS3bUv2MJTUiYiIiKUyZTLH0yZOhK++gscfh/37rY4qEVq2NMlclixmtMbYscm6HaukTkRERCxns5lpXL/9BtHRpqDim2+sjioRChc2zYq7dzflvs8+CxcvJsulldSJiIiI2yhTBrZuhWbNTLeQLl3g2jWro3JQ2rTw6aemp91PP5n+LZs3J/llldSJiIiIW8mY0WzDTp0Ks2aZqRR79lgdVSI0bWrGaeTMCVWrwpgxSbodq6RORERE3I7NBp07w5YtJg8KC4Np06yOKhEKFIANG6B3b1Mh27ix6b6cBJTUiYiIiNsqWdLsXLZqZZK8jh3h6lWro3JQmjQwahQsXmzO25UvnyR7ykrqRERExK1lyGC2YqdPNzNky5eHXbusjioRnnnGNOZ7+20zlcLFlNSJiIiIR2jf3hRRpE5ttmOnTLF0gEPiBAVBhw5J8tZK6kRERMRjFCtm2p60bw9du5p/Xr5sdVTuQUmdiIiIeBQ/P9OoeOZMWLTIbMfu3Gl1VNZTUiciIiIeqXVrM2LMzw8qVoQJEzxwO9aFlNSJiIiIxypSBH79FZ5/Hnr0gOeeg6goq6OyhpI6ERER8Wjp0sG4cTBnDvz4I4SEmBGsKY2SOhEREfEKLVuaZC5LFqhUCcaOTVnbsUrq7uPatWvkz5+f/v37Wx2KiIiIOKBwYdPf98UXoVcvaNECLl60OqrkoaTuPoYNG0bFihWtDkNEREQSIW1a+Owz06h41SqzHbtli9VRJT0ldf9y8OBB9u3bR4MGDawORURERJzQrBls3w4BAVClCnz6qXdvx3pUUrd+/XoaNmxI7ty5sdlsLFy48J7XjBs3joIFC5IuXTpCQ0PZsGGDQ9fo378/w4cPd1HEIiIiYqWCBeHnn+Hll6FPH2jaFM6ftzqqpOFRSd3Vq1cpW7YsY8eOve/358yZQ58+fXjjjTfYvn071apV46mnnuLYsWPxrwkNDaVUqVL3PE6ePMmiRYsoUqQIRYoUSa6PJCIiIkksTRr4+GPTqHj9eihXzrRB8TY2u90zFyJtNhsLFiygSZMm8c9VrFiRkJAQxo8fH/9c8eLFadKkSYJW3wYPHsw333yDr68vV65c4ebNm7z66qsMGTLkgT8THR1NdHR0/NdRUVEEBQVx6dIlMmfOnLgPJyIiIknir79M0+ItW2D4cOjXD3wsXuKKiorC39/f6dzBo1bqHiYmJoaIiAjq1at31/P16tVj06ZNCXqP4cOHc/z4cY4ePcqoUaPo1q3bQxO62z/j7+8f/wgKCkr0ZxAREZGklT8/rFtnkrkBA6BRIzh3zuqoXMNrkrqzZ88SGxtLYGDgXc8HBgZy6tSpJLvu4MGDuXTpUvzj+PHjSXYtERERcV7q1PDRR/DDD2YbNjjYtEHxdKmsDsDVbDbbXV/b7fZ7nkuITp06Jeh1adOmJW3atA6/v4iIiFirQQPYscNsx9aoAe+/DwMHWr8dm1geGva9AgIC8PX1vWdV7syZM/es3iWF8PBwSpQoQVhYWJJfS0RERFwjb15Ys8Ykc4MHw9NPwz//WB1V4nhNUpcmTRpCQ0NZuXLlXc+vXLmSypUrJ/n1e/bsyZ49e9iSErobioiIeJFUqeCDD8zc2IgIsx27fr3VUTnOo5K6K1eusGPHDnbs2AHAkSNH2LFjR3zLkn79+jF58mSmTp3K3r176du3L8eOHaN79+4WRi0iIiKeoH59sx372GNQq5bZjo2NtTqqhPOoM3Vbt26lVq1a8V/369cPgI4dOzJt2jRatWrFuXPnePfdd4mMjKRUqVIsXbqU/PnzWxWyiIiIeJDcuc1osffegyFDTKXsN99AMpzkcprH9qlzN+Hh4YSHhxMbG8uBAwfUp05ERMTDrV4NbduCzQYzZkDt2klzHVf1qVNS52KuujEiIiJivVOnoF07+Okns3L31lvg6+vaa6j5sIiIiEgSe+QRWL4c3nnHbMk+8QRERlod1f0pqRMRERF5CF9fs0K3ejXs22eqY//VbMMtKKkTERERSYCaNU11bHCwqZR98024dcvioO6gpM5F1HxYRETE++XMCcuWmXYnw4dDnTpw4oTVURkqlHAxFUqIiIikDBs2mBFj0dHw9dfw5JOJex8VSoiIiIhYqFo1sx1boQI89RQMGgQ3b1oXj5I6ERERkUQKCIDFi2HECBg1ypy7O37cmliU1ImIiIg4wccHBgww27HHj5tCiiVLLIgj+S/pnVQoISIikrJVqmS2Y6tUgYYNoX//5N2OVaGEi6lQQkREJGWz2+GTT2DgQAgNhTlz4GFj6FUoISIiIuKGbDbo2xd+/tmMGQsOhkWLkv66SupEREREkkDFirB9O9SqBU2aQJ8+EBOTdNdTUiciIiKSRLJmhXnz4LPPYPx4c97uzz+T5lpK6kRERESSkM0GvXrBpk1w/jxUrgzXrrn+Oqlc/5YpU3h4OOHh4cTGxlodioiIiLih0FDYtg02b4b06V3//qp+dTFVv4qIiIgjXJU7aKXOxW7nyFFRURZHIiIiIp7gds7g7DqbkjoXu3z5MgBBQUEWRyIiIiKe5PLly/j7+yf657X96mJxcXEUKVKEiIgIbDYbUVFRBAUFcfz4cUu2Y8PCwtiyZYtl75XQn/mv1z3s+w/63v2e//dz3nJ/vOHe/Pt5q+/Nv+NJ7vdxt/ujPzuO/0xCXqf74/r38cQ/O3a7ncuXL5M7d258fBJfw6qVOhfz8fEhTZo092TamTNntuQPlq+vr8uum5j3SujP/NfrHvb9B33vfs8/6LWefn+84d486Hmr7s2D4kmu93G3+6M/O47/TEJep/vj+vfx1D87zqzQ3aaWJkmgZ8+eVocQz5WxJOa9Evoz//W6h33/Qd+73/PudG/AdfF4w71xJKbkovvjeDzJxRvuzX+9RvdHf3Ycpe3XJKZqWPem++O+dG/cm+6Pe9P9cV9JeW+0UpfE0qZNy9ChQ0mbNq3Voch96P64L90b96b74950f9xXUt4brdSJiIiIeAGt1ImIiIh4ASV1IiIiIl5ASZ2IiIiIF1BSJyIiIuIFlNRZbMmSJRQtWpTHHnuMyZMnWx2O3KFp06ZkzZqVZ5991upQ5F+OHz9OzZo1KVGiBGXKlOHbb7+1OiT5n8uXLxMWFkZwcDClS5dm0qRJVock93Ht2jXy589P//79rQ5F/iVVqlQEBwcTHBxM165dHfpZVb9a6NatW5QoUYI1a9aQOXNmQkJC+O2338iWLZvVoQmwZs0arly5wldffcV3331ndThyh8jISE6fPk1wcDBnzpwhJCSE/fv3kyFDBqtDS/FiY2OJjo4mffr0XLt2jVKlSrFlyxayZ89udWhyhzfeeIODBw+SL18+Ro0aZXU4coeAgADOnj2bqJ/VSp2FNm/eTMmSJcmTJw+ZMmWiQYMGLF++3Oqw5H9q1apFpkyZrA5D7iNXrlwEBwcDkDNnTrJly8b58+etDUoAMxIpffr0ANy4cYPY2Fi0duBeDh48yL59+2jQoIHVoYiLKalzwvr162nYsCG5c+fGZrOxcOHCe14zbtw4ChYsSLp06QgNDWXDhg3x3zt58iR58uSJ/zpv3rycOHEiOUL3es7eG0larrw/W7duJS4ujqCgoCSOOmVwxb25ePEiZcuWJW/evAwcOJCAgIBkit77ueL+9O/fn+HDhydTxCmLK+5PVFQUoaGhVK1alXXr1jl0fSV1Trh69Sply5Zl7Nix9/3+nDlz6NOnD2+88Qbbt2+nWrVqPPXUUxw7dgzgvn97tdlsSRpzSuHsvZGk5ar7c+7cOTp06MDEiROTI+wUwRX3JkuWLOzcuZMjR44wc+ZMTp8+nVzhez1n78+iRYsoUqQIRYoUSc6wUwxX/Pk5evQoERERTJgwgQ4dOhAVFZXwAOziEoB9wYIFdz1XoUIFe/fu3e96rlixYvZBgwbZ7Xa7fePGjfYmTZrEf6937972GTNmJHmsKU1i7s1ta9assTdv3jypQ0zREnt/bty4Ya9WrZp9+vTpyRFmiuTMn53bunfvbp87d25ShZiiJeb+DBo0yJ43b157/vz57dmzZ7dnzpzZ/s477yRXyCmKK/78PPnkk/YtW7Yk+JpaqUsiMTExREREUK9evbuer1evHps2bQKgQoUK7Nq1ixMnTnD58mWWLl1K/fr1rQg3RUnIvRHrJOT+2O12OnXqRO3atWnfvr0VYaZICbk3p0+fjl9ZiIqKYv369RQtWjTZY02JEnJ/hg8fzvHjxzl69CijRo2iW7duDBkyxIpwU5yE3J8LFy4QHR0NwN9//82ePXsoVKhQgq+RynXhyp3Onj1LbGwsgYGBdz0fGBjIqVOnAFO2PHr0aGrVqkVcXBwDBw5UhVgySMi9Aahfvz7btm3j6tWr5M2blwULFhAWFpbc4aY4Cbk/GzduZM6cOZQpUyb+zMrXX39N6dKlkzvcFCUh9+bvv/+mS5cu2O127HY7L7/8MmXKlLEi3BQnof9tE2sk5P7s3buXF198ER8fH2w2G59++qlDHTGU1CWxf5+Rs9vtdz3XqFEjGjVqlNxhCf99b1SJbK2H3Z+qVasSFxdnRVjCw+9NaGgoO3bssCAque2//tt2W6dOnZIpIrnTw+5P5cqV+eOPPxL93tp+TSIBAQH4+vre87ejM2fO3JOlS/LSvXFvuj/uS/fGven+uLfkuD9K6pJImjRpCA0NZeXKlXc9v3LlSipXrmxRVAK6N+5O98d96d64N90f95Yc90fbr064cuUKhw4div/6yJEj7Nixg2zZspEvXz769etH+/btKV++PJUqVWLixIkcO3aM7t27Wxh1yqB74950f9yX7o170/1xb5bfnwTXyco91qxZYwfueXTs2DH+NeHh4fb8+fPb06RJYw8JCbGvW7fOuoBTEN0b96b74750b9yb7o97s/r+aPariIiIiBfQmToRERERL6CkTkRERMQLKKkTERER8QJK6kRERES8gJI6ERERES+gpE5ERETECyipExEREfECSupEREREvICSOhEREREvoKRORERExAsoqRMRERHxAkrqRERERLyAkjoRERERL6CkTkRERMQLKKkTERER8QJK6kRERES8QCqrAxARSemOHz/OhAkTiImJ4eLFi9StW5dWrVpZHZaIeBgldSIiFho/fjyffPIJ3333HaVLl+b69evUqVOH33//nWHDhgHw66+/kjp1akJDQy2OVkTcmbZfRUQsMnbsWHr37s3MmTMpXbo0AH5+fgwcOJAPP/yQ/fv3A/Dll18SFBRkZagi4gGU1ImIWGDfvn0MGDCALl263LMCV6FCBeLi4pg/fz6xsbGcO3eOnDlzWhSpiHgKJXUiIhYYPXo0N27coHv37vd8L0eOHAAcPXqUqVOn0r59++QOT0Q8kM1ut9utDkJExJNcvnyZZ599lujoaId+7sMPP+Txxx8HIGvWrKRLl47IyMj7vtZms9G8eXOuXr3K0qVLsdlsTsctIt5NSZ2ISDI7d+4cAQEBNGjQgB9++OG+r7HZbPj5+bFx40bKlSuXzBGKiCfS9quISDJLnTo1wEPPydlsNho2bKiETkQSTEmdiEgyy5w5M8HBwZw9e/a+3584cSJp0qQhLi4OgJs3byZneCLioZTUiYhYYMyYMaxbt44dO3bEP7d161a6dOlCpkyZ6NatGwcOHMBut9OmTRuuX79uXbAi4hF0pk5ExCJbt25l3Lhx+Pr64uvrS3BwMG3btiVTpkxcvnyZV199levXr/Piiy9StWpVq8MVETenpE5ERETEC2j7VURERMQLKKkTERER8QJK6kRERES8gJI6ERERES+gpE5ERETECyipExEREfECSupEREREvICSOhEREREvoKRORERExAsoqRMRERHxAkrqRERERLyAkjoRERERL6CkTkRERMQLKKkTERER8QL/B7RPmVP5LPxTAAAAAElFTkSuQmCC", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "execution_count": 9, "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, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "\n", "\n" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ ] }, { "cell_type": "code", "execution_count": 0, "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": 10, "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 }