diff --git a/doc/tutorial/leave-one-out-cross-validation.ipynb b/doc/tutorial/leave-one-out-cross-validation.ipynb index e0cbbbc..066a3be 100644 --- a/doc/tutorial/leave-one-out-cross-validation.ipynb +++ b/doc/tutorial/leave-one-out-cross-validation.ipynb @@ -18,14 +18,14 @@ "\n", "The variance computed in this way is the mean-squared-error, which consists of a bias term squared and the variance. Minimizing this thus finds the best compromise between the two terms.\n", "\n", - "We use the jackknife to implement part of this algorithm.\n", + "We first use the jackknife resampling to implement this algorithm manually, then we use the function `cross_validation` which simplifies the process.\n", "\n", "I demonstrate below how to select the optimal order for a polynomial model with leave-one-out (LOO) cross validation." ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -33,7 +33,7 @@ "import matplotlib.pyplot as plt\n", "from iminuit import Minuit\n", "from iminuit.cost import LeastSquares\n", - "from resample.jackknife import resample" + "from resample.jackknife import resample, cross_validation" ] }, { @@ -45,7 +45,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -74,24 +74,24 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# apply leave-one-out cross-validation\n", "\n", - "def predict(x, xloo, yloo, order):\n", + "def predict(xin, yin, xout, order):\n", " def model(x, par):\n", " return np.polyval(par, x)\n", "\n", " # least-squares cost function to fit a polynomial to the toy data\n", - " cost = LeastSquares(xloo, yloo, 0.1, model)\n", + " cost = LeastSquares(xin, yin, 0.1, model)\n", " m = Minuit(cost, np.zeros(order+1))\n", " m.strategy = 0 # faster, do not compute errors automatically\n", " m.migrad()\n", " assert m.valid\n", "\n", - " return model(x, m.values)\n", + " return model(xout, m.values)\n", "\n", "\n", "data = []\n", @@ -103,7 +103,7 @@ " deltas = []\n", "\n", " for iloo, (xloo, yloo) in enumerate(resample(x, y)):\n", - " yi_loo = predict(x[iloo], xloo, yloo, poly_order)\n", + " yi_loo = predict(xloo, yloo, x[iloo], poly_order)\n", " deltas.append(y[iloo] - yi_loo)\n", "\n", " variances.append(np.var(deltas))\n", @@ -112,7 +112,58 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA34AAAGeCAYAAADL8tO/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAACSyElEQVR4nOzdeVxU9foH8M8s7DuyK4u7IAoJYpqmImVklppLq6illdit+NleWjezupZ5b5KWGZq2aN20bpopaOGWsoipuMviwr4z7DPn9wfMJALKwMCZ5fN+veYlc+acM88A88gz5/t9vhJBEAQQERERERGR0ZKKHQARERERERF1LRZ+RERERERERo6FHxERERERkZFj4UdERERERGTkWPgREREREREZORZ+RERERERERo6FHxERERERkZFj4UdERERERGTkWPgREREREREZORZ+ZBI2bNgAiUSCzMxMsUPRCYlEgrfeekvr437//XdIJBL8/vvvOo+pvcaNG4dx48aJ9vxEHcU80oh5hKhjmEMaMYeIh4WfiTp06BDeeustlJaWih0KUZc6ffo07rnnHtja2sLZ2RmPP/44CgoKxA7LKDCPkCk4evQoFi5ciJCQEJiZmUEikYgdktFgDiFjp1KpsGHDBtx///3w9vaGjY0NAgMDsWzZMtTU1HR7PPJuf0bSC4cOHcLbb7+NOXPmwNHRUexwSEvV1dWQy/n2vZUrV67gzjvvhIODA5YvX47Kykp8+OGHOHHiBI4ePQpzc3OxQzRozCOGjXmkfXbu3IkvvvgCQ4cORZ8+fXDu3DmxQzIazCGGjTnk1qqqqjB37lzcfvvtePrpp+Hm5obDhw9j6dKlSEhIwN69e7v1wyT+tOiWVCoV6urqYGlpKXYo1ESffxYKhQI2Njbd8ly3+t1cvnw5FAoFUlJS4OPjAwAICwvDXXfdhQ0bNmDBggXdEicxj+gjff5Z6FMeeeaZZ/Dyyy/DysoKixYtYuEnEuYQ/aPPPwt9ySHm5uY4ePAgRo0apdk2f/58+Pn5aYq/iIiIbokT4FBPk/TWW2/hxRdfBAD07t0bEomk2ZhziUSCRYsW4euvv8bgwYNhYWGBXbt2tTkmOzMzExKJBBs2bGi2/cyZM5g+fTqcnZ1haWmJ0NBQ/Pzzz7eMT32+Dz/8EB9//DF8fX1hZWWFsWPH4uTJky3237t3L8aMGQMbGxs4OjrigQcewOnTp2/6HFFRUXBxcUF9fX2Lx+6++24MHDhQc1/9/di+fTsCAwNhYWGBwYMHY9euXS2OPXbsGCIjI2Fvbw9bW1tMmDABf/75Z7N91GP8Dxw4gH/84x9wdXWFo6MjnnrqKdTV1aG0tBSzZ8+Gk5MTnJyc8NJLL0EQhGbnuHFcfVZWFhYuXIiBAwfCysoKPXr0wIwZMzo1j0Cb1/LHH39g4cKFcHNzQ69evTSPf/755+jbty+srKwQFhaG/fv3t/pctbW1WLp0Kfr16wcLCwt4e3vjpZdeQm1tbYvX3drvZlv++9//4r777tMUfQAQERGBAQMGYOvWrR35tlAT5hHmkfYwhjzi7u4OKyurDn8PqHXMIcwh7WHoOcTc3LxZ0ac2depUALjl74iu8YqfCZo2bRrOnTuHb7/9Fh9//DFcXFwAAK6urpp99u7di61bt2LRokVwcXGBn5+fVmPwT506hTvuuAM9e/bEK6+8AhsbG2zduhVTpkzBf//7X80v/M189dVXqKioQHR0NGpqavDvf/8b4eHhOHHiBNzd3QEA8fHxiIyMRJ8+ffDWW2+huroan3zyCe644w6kpqbCz8+v1XM//vjj+Oqrr/Dbb7/hvvvu02zPzc3F3r17sXTp0mb7HzhwAD/++CMWLlwIOzs7/Oc//8GDDz6I7Oxs9OjRQ/Oax4wZA3t7e7z00kswMzPDZ599hnHjxuGPP/7AiBEjmp3z2WefhYeHB95++238+eef+Pzzz+Ho6IhDhw7Bx8cHy5cvx86dO7FixQoEBgZi9uzZbX6vkpKScOjQITz00EPo1asXMjMzsWbNGowbNw7p6emwtra+5ff7etq+loULF8LV1RVLliyBQqEAAKxfvx5PPfUURo0aheeffx6XLl3C/fffD2dnZ3h7e2uOValUuP/++3HgwAEsWLAA/v7+OHHiBD7++GOcO3cO27dvb/Zcrf1utubq1avIz89HaGhoi8fCwsKwc+dOrb4n1BzzCPPIrRhDHqGuwxzCHHIrxpxDcnNzAUDze99tBDJJK1asEAAIGRkZLR4DIEilUuHUqVPNtu/bt08AIOzbt6/Z9oyMDAGAEBcXp9k2YcIEYciQIUJNTY1mm0qlEkaNGiX079//prGpz2dlZSVcuXJFs/3IkSMCAOGFF17QbAsODhbc3NyEoqIizbbjx48LUqlUmD17tmZbXFxcs9erVCqFXr16CbNmzWr23CtXrhQkEolw6dKlZt8Pc3Nz4cKFC82eA4DwySefaLZNmTJFMDc3Fy5evKjZdu3aNcHOzk648847W8QyceJEQaVSabaPHDlSkEgkwtNPP63Z1tDQIPTq1UsYO3ZsszgBCEuXLtXcr6qqavF9PHz4sABA+OqrrzTb2voZ3kjb1zJ69GihoaFBs72urk5wc3MTgoODhdraWs32zz//XADQ7PVs2rRJkEqlwv79+5vFsHbtWgGAcPDgwWavu7XfzdYkJSW1eP1qL774ogCg2e8naY95hHnkZowhj9woOjpa4J9OusMcwhxyM8aYQ9QiIiIEe3t7oaSkpMPn6AgO9aRWjR07FgEBAR06tri4GHv37sXMmTNRUVGBwsJCFBYWoqioCBMnTsT58+dx9erVW55nypQp6Nmzp+Z+WFgYRowYoblSk5OTg7S0NMyZMwfOzs6a/YYOHYq77rrrpld0pFIpHn30Ufz888+oqKjQbP/6668xatQo9O7du9n+ERER6Nu3b7PnsLe3x6VLlwAASqUSu3fvxpQpU9CnTx/Nfp6ennjkkUdw4MABlJeXNzvnE0880WxC74gRIyAIAp544gnNNplMhtDQUM3ztOX6YUj19fUoKipCv3794OjoiNTU1Jsee6OOvJb58+dDJpNp7icnJyM/Px9PP/10swYqc+bMgYODQ7Njv//+e/j7+2PQoEGa35XCwkKEh4cDAPbt29ds//b+blZXVwMALCwsWjymHoev3oe6BvMI84ih5xESF3MIc4gx5pDly5cjPj4e77//frc3NWLhR626Mdlo48KFCxAEAW+++SZcXV2b3dTDFvLz8295nv79+7fYNmDAAM1Y8aysLABoNgZezd/fH4WFhZpL/a2ZPXs2qqursW3bNgDA2bNnkZKSgscff7zFvtfPEVNzcnJCSUkJAKCgoABVVVVtxqJSqXD58uWbnlOdhK4feqDern6etlRXV2PJkiXw9vaGhYUFXFxc4OrqitLSUpSVld302Bt15LXc+Pui/tnc+DM0MzNrlsAB4Pz58zh16lSL35UBAwYAaPm70t7fTfV/QDeOzQegaaHMeTtdi3mkOeaRm78WfcwjJC7mkOaYQ27+Wgwhh2zZsgVvvPEGnnjiCTzzzDMdOkdncI4ftaq1P4jbajerVCqb3VepVACAxYsXY+LEia0e069fv05G2HkBAQEICQnB5s2bMXv2bGzevBnm5uaYOXNmi32v/wTpesINE5210dY5W9t+q+d59tlnERcXh+effx4jR46Eg4MDJBIJHnroIc3Poyt1poBSqVQYMmQIVq5c2erjN/7n097n8vT0BND4aeyNcnJy4Ozs3OrVQNId5pHmmEduTh/zCImLOaQ55pCb0/ccsmfPHsyePRuTJk3C2rVrOxRnZ7HwM1EdWTPEyckJAFpMrFZ/oqKm/hTFzMysUy1qz58/32LbuXPnNBNofX19ATR+OnajM2fOwMXF5ZatfGfPno2YmBjk5OTgm2++waRJkzSvUxuurq6wtrZuMxapVNoiaejSDz/8gKioKHz00UeabTU1NR1aFFcXr0X9szl//rxmmATQOPQjIyMDQUFBmm19+/bF8ePHMWHCBJ2uZdOzZ0+4uroiOTm5xWNHjx5FcHCwzp7LVDGPNGIeaclY8gh1LeaQRswhLRlbDjly5AimTp2K0NBQbN26VbT1DznU00Spk5A2b0ZfX1/IZDIkJiY22/7pp582u+/m5oZx48bhs88+a/VqS0FBQbueb/v27c3G3x89ehRHjhxBZGQkgMYrOsHBwdi4cWOz13Hy5Ens3r0b99577y2f4+GHH4ZEIsFzzz2HS5cu4bHHHmtXbDeSyWS4++678dNPPzVrW5yXl4dvvvkGo0ePhr29fYfO3d7nv/GTuE8++aTFJ6DtPVdnX0toaChcXV2xdu1a1NXVabZv2LChxe/czJkzcfXqVaxbt67Feaqrq286ROZWHnzwQfzyyy/NhoMkJCTg3LlzmDFjRofPS42YRxoxj7R+LmPJI9R1mEMaMYe0fi5jySGnT5/GpEmT4Ofnh19++UXUEQe84meiQkJCAACvv/46HnroIZiZmWHy5Mk3/VTKwcEBM2bMwCeffAKJRIK+ffvil19+aXWMfGxsLEaPHo0hQ4Zg/vz56NOnD/Ly8nD48GFcuXIFx48fv2WM/fr1w+jRo/HMM8+gtrYWq1atQo8ePfDSSy9p9lmxYgUiIyMxcuRIPPHEE5oWyg4ODs3WlmmLq6sr7rnnHnz//fdwdHTEpEmTbnlMW5YtW4Y9e/Zg9OjRWLhwIeRyOT777DPU1tbiX//6V4fP2x733XcfNm3aBAcHBwQEBODw4cOIj4/XtHfWVmdfi5mZGZYtW4annnoK4eHhmDVrFjIyMhAXF9diXP3jjz+OrVu34umnn8a+fftwxx13QKlU4syZM9i6dSt+++23VpdkaI/XXnsN33//PcaPH4/nnnsOlZWVWLFiBYYMGYK5c+d26Jz0N+aRRswjrTOWPJKVlYVNmzYBgGYEwbJlywA0FiGtzcWi9mEOacQc0jpjyCEVFRWYOHEiSkpK8OKLL2LHjh3NHu/bty9Gjhyp9Xk7rFt7iJJeeeedd4SePXsKUqm0WXthAEJ0dHSrxxQUFAgPPvigYG1tLTg5OQlPPfWUcPLkyRYtlAVBEC5evCjMnj1b8PDwEMzMzISePXsK9913n/DDDz/cNC51C+UVK1YIH330keDt7S1YWFgIY8aMEY4fP95i//j4eOGOO+4QrKysBHt7e2Hy5MlCenp6s31ubKF8va1btwoAhAULFrQaT1vfD19fXyEqKqrZttTUVGHixImCra2tYG1tLYwfP144dOhQq7EkJSU127506VIBgFBQUNBse1RUlGBjY9MiputbKJeUlAhz584VXFxcBFtbW2HixInCmTNnWsTY3hbKnX0tap9++qnQu3dvwcLCQggNDRUSExOFsWPHtmgJXVdXJ3zwwQfC4MGDBQsLC8HJyUkICQkR3n77baGsrKzZ627rd7MtJ0+eFO6++27B2tpacHR0FB599FEhNzdXq3NQ25hHGjGPtM4Y8oj69bZ2uzEG0h5zSCPmkNYZeg5R/x61dbvxZ9fVJILQiRmhRF0gMzMTvXv3xooVK7B48eIuf76ffvoJU6ZMQWJiIsaMGdPlz0dEXY95hIg6gzmEjBHn+JHJW7duHfr06YPRo0eLHQoRGSjmESLqDOYQ6g6c40cm67vvvsNff/2FHTt24N///jc7wRGR1phHiKgzmEOoO7HwI5P18MMPw9bWFk888QQWLlwodjhEZICYR4ioM5hDqDtxjh8REREREZGR4xw/IiIiIiIiI8fCj4iIiIiIyMix8CO6iQ0bNkAikSAzM1PsULQ2btw4jBs3TuwwiEwe8wgRdQZzCOkKCz8iER06dAhvvfUWSktLRY9j9OjRsLa2hoeHB/7xj3+gsrJS1JiIqH30IY/s3r0bTzzxBAIDAyGTyeDn5ydaLESkHbFzSFVVFWJjY3H33XfD09MTdnZ2uO2227BmzRoolUpRYjJWLPyIRHTo0CG8/fbbov7BlpaWhgkTJqCqqgorV67Ek08+ic8//xwzZswQLSYiaj99yCPffPMNvvnmGzg4OMDLy0u0OIhIe2LnkEuXLuHZZ5+FIAiIiYnBhx9+iN69e2PhwoWYN2+eKDEZKy7nQGTiXnvtNTg5OeH333+Hvb09AMDPzw/z58/H7t27cffdd4scIRHpu+XLl2PdunUwMzPDfffdh5MnT4odEhEZCA8PD5w4cQKDBw/WbHvqqacwb948xMXF4c0330S/fv1EjNB48IofGZ2Kigo8//zz8PPzg4WFBdzc3HDXXXchNTW12X5HjhzBPffcAwcHB1hbW2Ps2LE4ePBgu57j119/xZgxY2BjYwM7OztMmjQJp06darHfmTNnMHPmTLi6usLKygoDBw7E66+/DgB466238OKLLwIAevfuDYlE0mIM/+bNmxESEgIrKys4OzvjoYcewuXLl1s8z+eff46+ffvCysoKYWFh2L9/f7teR3l5Ofbs2YPHHntMU/QBwOzZs2Fra4utW7e26zxExoZ5pP15BAC8vLxgZmbW7v2JjB1zSPtziIuLS7OiT23q1KkAgNOnT7frPHRrvOJHRufpp5/GDz/8gEWLFiEgIABFRUU4cOAATp8+jWHDhgEA9u7di8jISISEhGDp0qWQSqWIi4tDeHg49u/fj7CwsDbPv2nTJkRFRWHixIn44IMPUFVVhTVr1mD06NE4duyYZm7LX3/9hTFjxsDMzAwLFiyAn58fLl68iP/973949913MW3aNJw7dw7ffvstPv74Y7i4uAAAXF1dAQDvvvsu3nzzTcycORNPPvkkCgoK8Mknn+DOO+/EsWPH4OjoCABYv349nnrqKYwaNQrPP/88Ll26hPvvvx/Ozs7w9va+6ffqxIkTaGhoQGhoaLPt5ubmCA4OxrFjxzryIyAyeMwj7c8jRNQSc0jnc0hubi4AaGIiHRCIjIyDg4MQHR3d5uMqlUro37+/MHHiREGlUmm2V1VVCb179xbuuusuzba4uDgBgJCRkSEIgiBUVFQIjo6Owvz585udMzc3V3BwcGi2/c477xTs7OyErKysFs+vtmLFimbnV8vMzBRkMpnw7rvvNtt+4sQJQS6Xa7bX1dUJbm5uQnBwsFBbW6vZ7/PPPxcACGPHjm3z+yAIgvD9998LAITExMQWj82YMUPw8PC46fFExop5pP155EaTJk0SfH19tTqGyNgwh3Q8hwiCINTW1goBAQFC7969hfr6eq2Pp9ZxqCcZHUdHRxw5cgTXrl1r9fG0tDScP38ejzzyCIqKilBYWIjCwkIoFApMmDABiYmJUKlUrR67Z88elJaW4uGHH9YcV1hYCJlMhhEjRmDfvn0AgIKCAiQmJmLevHnw8fFpdg6JRHLL1/Djjz9CpVJh5syZzZ7Hw8MD/fv31zxPcnIy8vPz8fTTT8Pc3Fxz/Jw5c+Dg4HDL56murgYAWFhYtHjM0tJS8ziRqWEeaX8eIaKWmEM6l0MWLVqE9PR0rF69GnI5ByjqCr+TZHT+9a9/ISoqCt7e3ggJCcG9996L2bNno0+fPgCA8+fPAwCioqLaPEdZWRmcnJxabFcfGx4e3upx6nlyly5dAgAEBgZ26DWcP38egiCgf//+rT6unkuTlZUFAC32MzMz07zem7GysgIA1NbWtnispqZG8ziRqWEeaX8eIaKWmEM6nkNWrFiBdevW4Z133sG9996r9fHUNhZ+ZHRmzpyJMWPGYNu2bdi9ezdWrFiBDz74AD/++CMiIyM1n6CtWLECwcHBrZ7D1ta21e3qYzdt2gQPD48Wj+vqUymVSgWJRIJff/0VMpms3fFpy9PTEwCQk5PT4rGcnBy2ZSeTxTxCRJ3BHNIxGzZswMsvv4ynn34ab7zxhs7Pb+pY+JFR8vT0xMKFC7Fw4ULk5+dj2LBhePfddxEZGYm+ffsCaPxELCIiQqvzqo91c3O76bHqT7hu1dK8raEWffv2hSAI6N27NwYMGNDm8b6+vgAaP5W7/pO/+vp6ZGRkICgo6KbPHxgYCLlcjuTkZMycOVOzva6uDmlpac22EZka5pH25REiah1ziHY55KeffsKTTz6JadOmITY2tl3HkHY4x4+MilKpRFlZWbNtbm5u8PLy0gxnDAkJQd++ffHhhx+isrKyxTkKCgraPP/EiRNhb2+P5cuXo76+vs1jXV1dceedd+LLL79EdnZ2s30EQdB8bWNjAwAtFk2dNm0aZDIZ3n777Wb7q48vKioCAISGhsLV1RVr165FXV2dZp8NGza0ayFWBwcHREREYPPmzaioqNBs37RpEyorK7mIO5kk5pFG7c0jRNQcc0gjbXJIYmIiHnroIdx55534+uuvIZWyROkKvOJHRqWiogK9evXC9OnTERQUBFtbW8THxyMpKQkfffQRAEAqleKLL75AZGQkBg8ejLlz56Jnz564evUq9u3bB3t7e/zvf/9r9fz29vZYs2YNHn/8cQwbNgwPPfQQXF1dkZ2djR07duCOO+7A6tWrAQD/+c9/MHr0aAwbNgwLFixA7969kZmZiR07diAtLQ1AY+IHgNdffx0PPfQQzMzMMHnyZPTt2xfLli3Dq6++iszMTEyZMgV2dnbIyMjAtm3bsGDBAixevBhmZmZYtmwZnnrqKYSHh2PWrFnIyMhAXFxcu8fVv/vuuxg1ahTGjh2LBQsW4MqVK/joo49w991345577unkT4TI8DCPaJ9H/vrrL/z8888AgAsXLqCsrAzLli0DAAQFBWHy5Mkd/nkQGRrmEO1ySFZWFu6//35IJBJMnz4d33//fbPHhw4diqFDh3b0x0HXE6GTKFGXqa2tFV588UUhKChIsLOzE2xsbISgoCDh008/bbHvsWPHhGnTpgk9evQQLCwsBF9fX2HmzJlCQkKCZp8bWyir7du3T5g4caLg4OAgWFpaCn379hXmzJkjJCcnN9vv5MmTwtSpUwVHR0fB0tJSGDhwoPDmm2822+edd94RevbsKUil0hbP9d///lcYPXq0YGNjI9jY2AiDBg0SoqOjhbNnzzY7x6effir07t1bsLCwEEJDQ4XExERh7Nix7W6hvH//fmHUqFGCpaWl4OrqKkRHRwvl5eXtOpbI2DCPaJ9H1K+xtVtUVNQtjycyJswh2uWQffv2tZk/AAhLly696fHUfhJBuOHaLRERERERERkVDqAlIiIiIiIyciz8iIiIiIiIjBwLPyIiIiIiIiPHwo+IiIiIiMjIsfAjIiIiIiIyciz8iIiIiIiIjBwXcDcwKpUK165dg52dHSQSidjhEJk0QRBQUVEBLy8vSKWG8zka8wiR/mAeIaLO0CaHsPAzMNeuXYO3t7fYYRDRdS5fvoxevXqJHUa7MY8Q6R/mESLqjPbkEBZ+BsbOzg5A4w/X3t5e5GiITFt5eTm8vb0170tDwTxCpD+YR4ioM7TJISz8DIx6OIW9vT0TLZGeMLRhTswjRPqHeYSIOqM9OcRwBpMTEZHG1KlT4eTkhOnTp4sdChERERkAFn5ERAboueeew1dffSV2GERERGQgWPgRERmgcePGGdycICIiIhIPCz8iIh1LTEzE5MmT4eXlBYlEgu3bt7fYJzY2Fn5+frC0tMSIESNw9OjR7g+UiIiITAYLPyIiHVMoFAgKCkJsbGyrj2/ZsgUxMTFYunQpUlNTERQUhIkTJyI/P1+zT3BwMAIDA1vcrl271l0vg4iIiIwIu3oSEelYZGQkIiMj23x85cqVmD9/PubOnQsAWLt2LXbs2IEvv/wSr7zyCgAgLS1NZ/HU1taitrZWc7+8vFxn5yYi0xIbG4vY2FgolUqxQyEiLfGKHxFRN6qrq0NKSgoiIiI026RSKSIiInD48OEuec733nsPDg4OmhsXXSaijoqOjkZ6ejqSkpLEDoWItMTCj4ioGxUWFkKpVMLd3b3Zdnd3d+Tm5rb7PBEREZgxYwZ27tyJXr163bRofPXVV1FWVqa5Xb58ucPxExERkWHiUE8iIgMUHx/f7n0tLCxgYWHRhdEQERGRvuMVPyLqtLKqepy8WgZBEMQORe+5uLhAJpMhLy+v2fa8vDx4eHiIFFXrzuVV4L8pV3CttFrsUIjIAFXXKXHwQiF2/JUjdihEBBZ+RNRJgiBg7oajuO+TA3h43Z84m1shdkh6zdzcHCEhIUhISNBsU6lUSEhIwMiRI0WMrKU3tp/E/31/HIcuFokdChEZoL+ulOLRL47g7f+d4geDRHqAhR8RdUpyVglSs0sBAH9eKsa9/9mPt/93CmXV9eIGJqLKykqkpaVpOnNmZGQgLS0N2dnZAICYmBisW7cOGzduxOnTp/HMM89AoVBounzqiwBPewDA6Rx2ASUi7QV5O0IulSC/ohZXSjhygEhsnONHRJ0SdzADADBxcGOzkt9O5SHuYCb+d/waXrpnEKYP6wWpVCJmiN0uOTkZ48eP19yPiYkBAERFRWHDhg2YNWsWCgoKsGTJEuTm5iI4OBi7du1q0fBFbOrCL/0aCz8i0p6lmQyDezrg+OVSpGSVwNvZWuyQiEwaCz8i6rCrpdX47VTjXLUX7hqAQR72SDxXgLf+dwqXChR46Ye/8M2RbLx9/2AEeTuKG2w3Gjdu3C2HNS1atAiLFi3qpog6xl99xS+3HIIgQCIxrQKeiDov1NcJxy+XIjmrGFNu6yl2OEQmjUM9iajDNh3OglIlYGSfHhjk0Vgk3DnAFbueuxOv3TsINuYypF0uxZRPD+KV//6FosraW5yR9El/d1vIpBKUVtUjt7xG7HCIyACF+joBAJIzS0SOhIhY+BFRh1TXKfFdUuOctTl3+DV7zFwuxYI7+2Lv4nGYeltPCALwXdJljP/wd2w8lIkGpUqEiElblmYy9HW1AcDhnkTUMSFNhd/ZvApU1Jju3G8ifcDCj4g65Ke0qyitqkcvJytE+Lc+N83d3hIfzwrG90+PhL+nPcprGrD051O475MDOHKJnSINgT8bvBBRJ7jZW8Lb2QqCABxragRGROJg4UdEWhMEAXEHMwEAUSP9ILtF85bhfs745dnReOeBwXCwMsOZ3ArM+vxPPPfdMeSWcQihPvu7syeX6SCijgn1dQbQ2AWaiMTDwo+ItHb4UhHO5lXAykyGmcO923WMTCrB4yP9sG/xODwywgcSCfBT2jWEf/Q71v5xEXUNHP6pj9RX/NJ5xY+IOmhY03DPVBZ+RKJi4UdEWlNf7XswpCccrMy0OtbZxhzLpw7Bz9GjMczHEVV1Srz/6xncsyoRv5/N74JoSS02NhYBAQEYPnx4u49RF36ZRQpU1TV0VWhEZMTUDV6OZZdwjjeRiFj4EZFWLhdXIf504xIOc0b5dfg8Q3o54IenR+GjGUFwsbXApUIF5sQlYf5XycguqtJRtHS96OhopKenIykpqd3HuNpZwNXOAoIAnMnlcE8i0t4AdzvYWcihqFMyjxCJiIUfEWll46FMCAIwpr8L+rnZdepcUqkED4b0wt7FY/HE6N6QSSXYk56HiI//wMo951Bdp9RR1NQZ/lzInYg6QSaVINjHEQCQwuGeRKJh4UdE7aaobcCW5MsAgLk3LOHQGfaWZnjzvgDsem4M7ujXA3UNKvwn4TwiVv6BXSdzbrkYOnWtAHb2JKJOUjd4YeFHJB4WfkTUbj8eu4qKmgb49bDGuAFuOj9/f3c7bH5iBD59dBi8HCxxtbQaT29OxePrj+JCPocHicXfs/HKLhu8EFFHhfo1zvNj4UckHhZ+RNQuKpWADQczAABRo/wgvcUSDh0lkUhw7xBPJPzfOPwjvB/M5VIcuFCIe1btx7s70rkAsAjUV/zO5lZApeLVVyLSXpC3I6QS4GppNXLKqsUOh8gksfATUWlpKUJDQxEcHIzAwECsW7dO7JCI2nTgQiEuFihgayHH9JBeXf58VuYyxNw9EPEvjEWEvzsaVALW7c9A+Ed/4MfUKxz+2Y16u9jAQi5FVZ0SWcVsvENE2rO1kGvmC/OqH5E4WPiJyM7ODomJiUhLS8ORI0ewfPlyFBUViR0WUavimq72TQ/pBTtL7ZZw6AyfHtb4IioUcXOHo7eLDQoqahGz9Timrz2Mk1fLui0OUyaXSTHQo2m4Jxu8EJm0jiwLo6Ze1iE5k4UfkRhY+IlIJpPB2toaAFBbWwtBEHgVg/RSRqEC+84WQCJpHOYphvED3bDr+TF46Z6BsDaXISWrBJNXH8Dr206gRFEnSkymxN+DDV6IqGPLwqipF3LnFT8icYhe+Pn5+UEikbS4RUdHt7q/UqnEm2++id69e8PKygp9+/bFO++8o/OCKTExEZMnT4aXlxckEgm2b9/e6n6xsbHw8/ODpaUlRowYgaNHj2r1PKWlpQgKCkKvXr3w4osvwsXFRQfRE+nWxkOZABqLr94uNqLFYSGXYeG4fkj4v7GYHOQFQQC+PpKN8R/9js1/ZkHJ+WddJsCLhR8RdU6oX2Nnz/ScclTVNYgcDZHpEb3wS0pKQk5Ojua2Z88eAMCMGTNa3f+DDz7AmjVrsHr1apw+fRoffPAB/vWvf+GTTz5p8zkOHjyI+vqWDSHS09ORl5fX6jEKhQJBQUGIjY1t87xbtmxBTEwMli5ditTUVAQFBWHixInIz8/X7KOev3fj7dq1awAAR0dHHD9+HBkZGfjmm2/ajIdILBU19fi+aQmHzizYrkueDlb45OHb8N2C2zHIww6lVfV4Y/tJ3L/6AFKyisUOzyhp1vJj4UdEHdTT0QqeDpZQqgSkXS4VOxwikyN64efq6goPDw/N7ZdffkHfvn0xduzYVvc/dOgQHnjgAUyaNAl+fn6YPn067r777javtKlUKkRHR+ORRx6BUvn3YtBnz55FeHg4Nm7c2OpxkZGRWLZsGaZOndpm7CtXrsT8+fMxd+5cBAQEYO3atbC2tsaXX36p2SctLQ0nT55scfPy8mp2Lnd3dwQFBWH//v1tPh+RGH5IuQJFnRL93Gwxpr9+XZG+vU8P/PLsaLw1OQB2lnKculaOB9ccRsyWNOSX14gdnlEZ1LSkQ05ZDUqrOLSWiDpGM9yT8/yIup3ohd/16urqsHnzZsybNw8SSeut4keNGoWEhAScO3cOAHD8+HEcOHAAkZGRre4vlUqxc+dOHDt2DLNnz4ZKpcLFixcRHh6OKVOm4KWXXupwrCkpKYiIiGj2XBERETh8+HC7zpGXl4eKisa1ycrKypCYmIiBAwe2um9nJlMTdZRKJWiGeUaN8mvzfSkmuUyKOXf0xr7F4zAr1BsSSeN6g+Ef/YF1iZdQr1SJHaJRsLc0g7ezFQBe9SOijlM3eEnJZuFH1N30qvDbvn07SktLMWfOnDb3eeWVV/DQQw9h0KBBMDMzw2233Ybnn38ejz76aJvHeHl5Ye/evThw4AAeeeQRhIeHIyIiAmvWrOlwrIWFhVAqlXB3d2+23d3dHbm5ue06R1ZWFsaMGYOgoCCMGTMGzz77LIYMGdLqvp2ZTE3UUb+fy0dmURXsLOV4cFhPscO5KRdbC3wwfSi2LbwDQb0cUFnbgHd3nsY9qxKx/3yB2OEZBXWDF3b2JKKOCvVtnOeXmlXCdUGJuplc7ACut379ekRGRrYYBnm9rVu34uuvv8Y333yDwYMHIy0tDc8//zy8vLwQFRXV5nE+Pj7YtGkTxo4diz59+mD9+vWiX70ICwtDWlqaqDEQ3UzcwUwAwEPDvWFtrlfpok3B3o7YtvAO/JByBR/sOoOLBQo8vv4oIgM98Pokf/RyshY7RIMV4GWP3el5OJ1TIXYoRGSgBnnawcpMhvKaBpzPr9QsFUNEXU9vrvhlZWUhPj4eTz755E33e/HFFzVX/YYMGYLHH38cL7zwAt57772bHpeXl4cFCxZg8uTJqKqqwgsvvNCpeF1cXCCTyVo0Y8nLy4OHh0enzk2kD87nVWD/+UJIJcDskX5ih6MVqVSCmcO9sXfxOMwZ5QeZVIJfT+YiYuUf+Hf8edTUK299EiPU2SHjbPBCRJ1lJpMi2NsRAJd1IOpuelP4xcXFwc3NDZMmTbrpflVVVZBKm4ctk8mgUrU9j6ewsBATJkyAv78/fvzxRyQkJGDLli1YvHhxh+M1NzdHSEgIEhISNNtUKhUSEhIwcuTIDp+XSF9saJrbF+HvDm9nw7xK5mBlhrfuH4wd/xiNEb2dUVOvwsfx53DXx39g96lck1s3s7NDxgOaCr8L+RWoa+DcSSLqmFC/poXc2YWZqFvpReGnUqkQFxeHqKgoyOXNh5OtXr0aEyZM0NyfPHky3n33XezYsQOZmZnYtm0bVq5c2Wb3TZVKhcjISPj6+mLLli2Qy+UICAjAnj17EBcXh48//rjV4yorK5GWlqYZipmRkYG0tDRkZ2dr9omJicG6deuwceNGnD59Gs888wwUCgXmzp3bye8IkbjKqurxY+pVAMCcO/zEDUYHBnnY47sFt+OTh2+Dh70lLhdXY8GmFETFJeFiQaXY4RmMXk5WsLOQo14p8PtGRB3GhdyJxKEXk3bi4+ORnZ2NefPmtXissLAQFy9e1Nz/5JNP8Oabb2LhwoXIz8+Hl5cXnnrqKSxZsqTVc0ulUixfvhxjxoyBubm5ZntQUBDi4+Ph6ura6nHJyckYP3685n5MTAwAICoqChs2bAAAzJo1CwUFBViyZAlyc3MRHByMXbt2tWj4QmRotiRno7peiUEedhjZp4fY4eiERCLB5CAvhA9yQ+y+C/hifwYSzxXgnlWJeGJ0Hzwb3g82FnqREvWWRCKBv6c9jmYWI/1auWboJxGRNob5OEEiAbKKqlBQUQtXOwuxQyIyCRLB1MY6Gbjy8nI4ODigrKwM9vb8o4t0T6kScOe/9uFqaTXenzYED4X5iB1Sl8goVOCf/zuFfWcbO36621vgtXv9cX+QV7sbPxnq+7EzcS/96SQ2Hs7Ck6N74437ArooQiLTYYp5BAAmfpyIs3kVWPtYCO4JZG8Eoo7S5r2oF0M9iUh/xJ/Ow9XSajham2HKbfq9hENn9HaxQdzcMKyPCoVvD2vkldfiue/SMOvzP3GazUvaFODV+J/K6Vx+j4io40L81MM9Oc+PqLuw8COiZuIOZgAAHg7zgaWZTORout4Ef3f89vydWHz3AFiaSXE0oxiT/rMfS386ibKqerHD0zuazp7Xyk2uOQ4R6U6ID+f5EXU3Fn5EpHE6pxx/XiqGTCrB47f7ih1Ot7E0k2FReH8k/N84TBriCZUAbPozC1dKq8QOTe8McLeDTCpBSVU98sprxQ6HiAyUurPnyavlJrvEDlF3Y+FHRBobmhZsv2ewB7wcrcQNRgQ9Ha0Q++gwfPPkCLx0zyAM9nIQOyS9Y2kmQx8XGwDgkFgi6jAfZ2u42FqgTqnCiatlYodDZBJY+BERAKBYUYftacazhENnjOrngqfH9hU7DL3FhdyJqLMkEglCfB0BcLgnUXdh4UdEAIBvj2ajtkGFwJ72CG1aY4moNeoGLyz8iKgzQn2dAQDJmSz8iLoDCz8iQr1Shc1/ZgEA5ozq3e7lDMg0qa/4nb7Gwo+IOk7d2TM1u4TNooi6AQs/IsLuU3nIKauBi605Jgd5ih0O6Tl/TzsAQEaRAlV1DSJHQ0SGarCXPczlUhQr6pBRqBA7HCKjx8KPiDRLODwS5gMLufEv4UCd42ZnCRdbCwgCcDa3QuxwiMhAWchlCOrV2EQrmfP8iLocCz8iE3fiShmSs0ogl0rwmAkt4UCdo77qx3l+RNQZIU3z/FI4z4+oy7HwIzJxcYcar/ZNGuoJN3tLkaMhQ6Fu8MIlHYioM0KamomlZLPwI+pqLPyITFhBRS1+OZ4DAJgzyk/cYKjLxcbGIiAgAMOHD+/0uQLUDV5yONSTiDpOXfhdyK9EaVWdyNEQGTcWfkQm7Jsj2ahTqhDs7YjbfLiEg7GLjo5Geno6kpKSOn0uTWfPnHKoVOzGR0Qd42xjjj6uNgC4nh9RV2PhR2Si6hpU2HykcQmHuSa+YDtpr4+LDczlUlTVKZFdXCV2OERkwEKaPnhk4UfUtVj4EZmoX0/moKCiFm52FogM5BIOpB25TIqB7o0NXjjPj4g6I7RpPT929iTqWiz8iEzUlwczAQCP3e4LczlTAWmPnT2JSBfUnT2PXy5FXYNK5GiIjBf/2iMyQceyS3D8cinMZVI8MsJH7HDIQAV4srMnEXVeHxcbOFqbobZBxQ+SiLoQCz8iExTXdLVvcpAXXGwtxA2GDJY/O3sSmRxddgdWk0olmnl+yZnFOjsvETXHwo/IxOSV12DnicYlHNjUhTpjUFPhd7W0mm3YiUyELrsDXy/Ejw1eiLoaCz8iE7P5zyw0qAQM93NCYE8HscMhA+ZgZYZeTlYAeNWPiDpHc8UvqwSCwCViiLoCCz8iE1JTr8Q3R7IBAHNG9RY5GjIG6uGenJdDRJ0R5O0IM5kEBRW1uFJSLXY4REaJhR+RCfnlrxwUKerg6WCJiYPdxQ6HjAAbvBCRLliayTDYq3EUSnIW5/kRdQUWfkQmQhAExB3MAAA8PtIXchnf/tR5/iz8iEhHQnw5z4+oK/EvPyITkZxVglPXymEhl+Lh4VzCgXRDfcXvfF4l6pVcf4uIOi7UV93Zk4UfUVdg4UdkItRX+6be1hNONuYiR0PGopeTFews5KhTqnCxoFLscIjIgKmv+J3Nq0B5Tb3I0RAZHxZ+RCbgamk1fjuVBwCYwyUcSIekUgkGedoB4HBPIuocN3tLeDtbQRCAtOxSscMhMjos/IhMwKbDWVCqBIzs0wODPOzFDoeMjKaz5zUWfkTUOaG+zgAapycQkW6x8CMyctV1SnyX1LSEA6/2URf4u7Mn1/Ijos75u8ELO3sS6RoLPyIj91PaVZRW1aOXkxUi/LmEA+ne9Z09ufAyEXWGuvBLyy5FAxtGEekUCz8iI9a4hEMmACBqpB9kUom4AZFRGuhhB6kEKFLUIb+iVuxwiMiADXC3g52FHIo6Jc7kchQBkS6x8CMyYocvFeFsXgWszGSYOdxb7HDISFmaydDH1RYAkM4GL0TUCTKpBLdxPT+iLsHCj8iIqa/2PRjSEw5WZuIGQ6KLjY1FQEAAhg8frvNzs8ELEelKiA8LP6KuwMKPyEhdLq5C/OmmJRxG+YkbDOmF6OhopKenIykpSefnDrhunh8RUWeE+rHwI+oKLPyIjNTGQ5kQBGBMfxf0c7MTOxwycv5cy4+IdCTY2xFSSeMatDll1WKHQ2Q0WPgRGSFFbQO2JF8GAMzlEg7UDdRX/DIKFaiuU4ocDREZMhsLuWb4OK/6EekOCz8iI/Rj6hVU1DTAr4c1xg1wEzscMgGudhZwsTWHSgDO5rETHxF1TmhTg5fkTBZ+RLrCwo/IyKhUAjYcygQARI3yg5RLOFA3kEgkzdbzIyLqjBA/ZwC84kekSyz8iIzMgQuFuFiggK2FHNNDeokdDpkQdvYkIl1RL+SenlOOqroGkaMhMg4s/IiMTNzBDADA9JBesLPkEg7UfdjZk4h0paejFTwdLKFUCUi7XCp2OERGgYUfkRHJKFRg39kCSCRcwoG6n/qK35ncCqhUgsjREJGhU1/1S+E8PyKdYOFHZEQ2Ns3tGz/QDX4uNuIGQyanj6sNzGVSVNY24HJJldjhEJGBUzd4Sclm4UekCyz8iIxERU09vm9awoFX+0gMZjIpBnjYAuBwTyLqvBDfxgYvqVklHEVApAMs/IiMxPfJV6CoU6Kfmy3G9HcROxwyUf4ebPBCRLrh72kHa3MZymsacD6/UuxwiAweCz8iI6BSCdh4OBNA4xIOEgmXcCBxBHg1FX45XMuPiDpHLpMi2NsRAJd1INIFFn5ERuD3c/nIKqqCnaUcDw7rKXY4ZMK4lh8R6ZK6wUtyVrHIkRAZPhZ+REYg7mAmAOCh4d6wNpeLGwyZNPVQz6ul1Sirqhc5GiIydJrOnrziR9RpLPyIDNz5vArsP18IqQSYPdJP7HDIxDlYm6GnoxUA4HQur/oRUefc5uMEiQTIKqpCQUWt2OEQGTQWfkQGbkPTEg4R/u7wdrYWNxgicLgnEemOg5UZBrjZAeBVP6LOYuFHZMDKqurxY+pVAMCcO/zEDYaoSYBn4x9p7OxJRLoQ4qce7sl5fkSdwcJPRKWlpQgNDUVwcDACAwOxbt06sUMiA7MlORvV9UoM8rDDyD49xA6HCMDfnT051JOIdCFU0+CFV/yIOoNdIERkZ2eHxMREWFtbQ6FQIDAwENOmTUOPHvwDnm5NqRKw8VAWgMYF27mEA+kL9VDPc3mVqFeqYCbjZ4xE1HHqBi8nr5ahpl4JSzOZyBERGSb+bywimUwGa+vGOVm1tbUQBAGCIIgcFRmK+NN5uFpaDUdrM0y5jUs4kP7wdrKGrYUcdQ0qXCpQiB0OERk4H2druNhaoF4p4MTVMrHDITJYelH4+fk1Xq248RYdHa2T/TsiMTERkydPhpeXFyQSCbZv397qfrGxsfDz84OlpSVGjBiBo0ePavU8paWlCAoKQq9evfDiiy/CxcVFB9GTKYg7mAEAeDjMh59+kl6RSiUY5NE4z48NXoiosyQSyd/DPTM53JOoo/Si8EtKSkJOTo7mtmfPHgDAjBkzdLL/wYMHUV/fcj2p9PR05OXltXqMQqFAUFAQYmNj24x7y5YtiImJwdKlS5GamoqgoCBMnDgR+fn5mn3U8/duvF27dg0A4OjoiOPHjyMjIwPffPNNm/EQXe90Tjn+vFQMmVSCx2/3FTscohbUwz3TWfgRkQ5wPT+iztOLOX6urq7N7r///vvo27cvxo4d2+n9VSoVoqOj0b9/f3z33XeQyRqvjJw9exbh4eGIiYnBSy+91OK4yMhIREZG3jTulStXYv78+Zg7dy4AYO3atdixYwe+/PJLvPLKKwCAtLS0m55Dzd3dHUFBQdi/fz+mT5/ermPIdG1oWrD9nsEe8GpaM43oVmJjYxEbGwulUtnlz6Vp8MLCj4h0QN3ZMzW7BIIgcF47UQfoxRW/69XV1WHz5s2YN29eu97Ut9pfKpVi586dOHbsGGbPng2VSoWLFy8iPDwcU6ZMabXoa2+cKSkpiIiIaPZcEREROHz4cLvOkZeXh4qKCgBAWVkZEhMTMXDgwFb3jY2NRUBAAIYPH96heMl4FCvqsD2NSziQ9qKjo5Geno6kpKQufy7NFb9r5Zy7TESdFujlAHO5FMWKOlwq5Nxhoo7Qu8Jv+/btKC0txZw5c3S2v5eXF/bu3YsDBw7gkUceQXh4OCIiIrBmzZoOx1lYWAilUgl3d/dm293d3ZGbm9uuc2RlZWHMmDEICgrCmDFj8Oyzz2LIkCGt7tudf7CRfvv2aDZqG1QI7GmvmfNApG8GuttBKgGKFHUoqKgVOxwiMnDmcimCejkA4HBPoo7Si6Ge11u/fj0iIyPh5eWl0/19fHywadMmjB07Fn369MH69etFHyYQFhbW7qGgRABQr1Rh85/qJRx6i/47TNQWK3MZervY4GKBAuk55XCztxQ7JCIycCG+zkjKLEFKZglmhnqLHQ6RwdGrK35ZWVmIj4/Hk08+qfP98/LysGDBAkyePBlVVVV44YUXOhWri4sLZDJZi2YseXl58PDw6NS5idqy+1Qecspq4GJrjslBnmKHQ3RT6uGep3MqRI6EiIzB3wu5F4scCZFh0qvCLy4uDm5ubpg0aZJO9y8sLMSECRPg7++PH3/8EQkJCdiyZQsWL17c4VjNzc0REhKChIQEzTaVSoWEhASMHDmyw+cluhn1Eg6PhPnAQs4lHEi/sbMnEenSsKbC72KBAiWKOpGjITI8elP4qVQqxMXFISoqCnJ58xGoq1evxoQJE9q9/437RUZGwtfXF1u2bIFcLkdAQAD27NmDuLg4fPzxx60eV1lZibS0NM1QzIyMDKSlpSE7O1uzT0xMDNatW4eNGzfi9OnTeOaZZ6BQKDRdPol06cSVMiRnlUAuleAxLuFABoCdPYlIl5xtzNHH1QZAY3dPItKO3szxi4+PR3Z2NubNm9fiscLCQly8eLHd+19PKpVi+fLlGDNmDMzNzTXbg4KCEB8f32JpCLXk5GSMHz9ecz8mJgYAEBUVhQ0bNgAAZs2ahYKCAixZsgS5ubkIDg7Grl27WjR8IdKFuEONV/smDfXkfCkyCAFNV/wuFVSipl4JSzNepSaizgn1dcKlAgWSs0owwZ9/bxFpQyKwz7ZBKS8vh4ODA8rKymBvby92ONRNCipqccf7e1GnVGHbwlG4zYfdPPWBob4fuytuQRAQuiweRYo6/BR9B4K8HbvsuYgMFfOIdrYkZePl/55AWG9nbH2KU2uItHkv6s1QTyJq2zdHslGnVCHY25FFHxkMiURyXYMXDvck0jdTp06Fk5MTpk+fLnYo7Rbi6wwAOH65FHUNKpGjITIsLPyI9FxdgwqbjzQu4TCXC7aTgfH3tAPABi9E+ui5557DV199JXYYWunragNHazPUNqhw6lqZ2OEQGRQWfkR67teTOSioqIWbnQUiA7mEAxkWNngh0l/jxo2DnZ2d2GFoRSKRIKRp5AsXcifSDgs/Ij335cFMAMBjt/vCXM63LBmW69fyU6k4pZyovRITEzF58mR4eXlBIpFg+/btLfaJjY2Fn58fLC0tMWLECBw9erT7AxVBiB8LP6KO6NBfkRcvXsQbb7yBhx9+GPn5+QCAX3/9FadOndJpcESm7lh2CY5fLoW5TIpHRviIHY7RY27Tvb6utjCXSVFZ24ArJdVih0PUpXSZQxQKBYKCghAbG9vq41u2bEFMTAyWLl2K1NRUBAUFYeLEiZrnBYDg4GAEBga2uF27dq1jL1BPhDbN80vOKgF7FBK1n9aF3x9//IEhQ4bgyJEj+PHHH1FZWQkAOH78OJYuXarzAIlMWVzT1b7JQV5wsbUQNxgjx9zWNcxkUvR3twXAeX5k3HSdQyIjI7Fs2TJMnTq11cdXrlyJ+fPnY+7cuQgICMDatWthbW2NL7/8UrNPWloaTp482eLm5eWldTy1tbUoLy9vdhPL0F4OMJNJUFBRyw+UiLSgdeH3yiuvYNmyZdizZ0+zdfHCw8Px559/6jQ4IlOWV16DnSdyALCpS3dgbus67OxJpqA7c0hdXR1SUlIQERGh2SaVShEREYHDhw/r9LnU3nvvPTg4OGhu3t7eXfI87WFpJsNgLwcAQHJWsWhxEBkarQu/EydOtPrpk5ubGwoLC3USFBEBm//MQoNKwHA/JwT2dBA7HKPH3NZ11IUfr/iRMevOHFJYWAilUgl39+YLmLu7uyM3N7fd54mIiMCMGTOwc+dO9OrV66ZF46uvvoqysjLN7fLlyx2OXxdCfRvn+SVncp4fUXvJtT3A0dEROTk56N27d7Ptx44dQ8+ePXUWGJEpq6lX4psj2QCAOaN632Jv0gXmtq4TwCt+ZAIMMYfEx8e3e18LCwtYWOjPlIMQXyd8cSCDDV6ItKD1Fb+HHnoIL7/8MnJzcyGRSKBSqXDw4EEsXrwYs2fP7ooYiUzOL3/loEhRB08HS0wc7H7rA6jTmNu6jrrwu1JSjbLqepGjIeoa3ZlDXFxcIJPJkJeX12x7Xl4ePDw8dPpc+krd2fNsXgXKa5hXiNpD68Jv+fLlGDRoELy9vVFZWYmAgADceeedGDVqFN54442uiJHIpAiCgLiDGQCAx0f6Qi7jEg7dgbmt6zhYm6GnoxUA4Ayv+pGR6s4cYm5ujpCQECQkJGi2qVQqJCQkYOTIkTp9Ln3lZmcJH2drCAJwLLtU7HCIDILWQz3Nzc2xbt06LFmyBCdOnEBlZSVuu+029O/fvyviIzI5yVklOHWtHBZyKR4eziUcugtzW9fy97TD1dJqnM4px4g+PcQOh0jndJ1DKisrceHCBc39jIwMpKWlwdnZGT4+PoiJiUFUVBRCQ0MRFhaGVatWQaFQYO7cubp6SXovxNcJ2cVVSMkqwdgBrmKHQ6T3tC781Ly9vUXt6ERkrNRX+6be1hNONua32Jt0jbmta/h72iP+dD4bvJDR01UOSU5Oxvjx4zX3Y2JiAABRUVHYsGEDZs2ahYKCAixZsgS5ubkIDg7Grl27WjR8MWYhvk7YduwqUtjZk6hdtB5D9uCDD+KDDz5osf1f//oXZsyYoZOgiEzV1dJq/Haqcc7GHC7h0K2Y27rW3w1eKkSOhKhr6DqHjBs3DoIgtLht2LBBs8+iRYuQlZWF2tpaHDlyBCNGjOjMSzA4oU3z/I5ll6JBqRI5GiL9p3Xhl5iYiHvvvbfF9sjISCQmJuokKCJTtelwFpQqASP79MAgD3uxwzEpzG1dS72kw9m8Cv6BRkbJVHJIbGwsAgICMHz4cLFDQX83O9hZyFFVp8SZXH6oRHQrWhd+lZWVzRYmVTMzM0N5OYfwEHVUdZ0S3yU1LeHAq33djrmta/k4W8PGXIa6BhUuFSrEDodI50wlh0RHRyM9PR1JSUlihwKZVILbmtbz47IORLemdeE3ZMgQbNmypcX27777DgEBAToJisgUbU+7itKqevRyskKEv+nM0dAXzG1dSyqVYBDX8yMjxhwiDs1C7iz8iG5J6+Yub775JqZNm4aLFy8iPDwcAJCQkIBvv/0W33//vc4DJDIFgiBgw8FMAEDUSD/IpBJxAzJBzG1dz9/TDilZJUjPKccDwfq5oDVRRzGHiENd+KWy8CO6Ja0Lv8mTJ2P79u1Yvnw5fvjhB1hZWWHo0KGIj4/H2LFjuyJGIqN3+FIRzuZVwMpMhpnD2VFSDMxtXS/A0wEAkH6NV/zI+DCHiCPI2xEyqQRXS6uRU1YNTwcrsUMi0lsdWs5h0qRJmDRpkq5jITJZcU1X+x4M6QkHKzNxgzFhxp7bYmNjERsbC6VSKcrz+3vaAWBnTzJexp5D9JGNhRz+nnY4ebUcyZklmBzEwo+oLR1ex6+urg75+flQqZp3Z/Px4YLTRNq4XFyF+NNNSziM8hM3GDLq3BYdHY3o6GiUl5fDwcGh259/oIcdJBKgsLIW+RU1cLOz7PYYiLqaMecQfRXq64yTV8uRklWCyUFeYodDpLe0LvzOnz+PefPm4dChQ822C4IAiUQi2ifJRIZq46FMCAIwpr8L+rnZiR2OyWJu63rW5nL0drHBpQIFTudUsPAjo8IcIp5hvk7YcCiTnT2JbkHrwm/OnDmQy+X45Zdf4OnpCYmETSiIOkpR24AtyZcBAHO5hIOomNu6h7+nfVPhV46xA1zFDodIZ0wlh4g9ZLw16gYv6TnlUNQ2wMaiwwPaiIya1u+MtLQ0pKSkYNCgQV0RD5FJ+TH1CipqGuDXwxrjBriJHY5JY27rHgGe9tjxVw4bvJDRMZUcIvaQ8dZ4OVrBy8ES18pqcPxKKUb1dRE7JCK9pPU6fgEBASgsLOyKWIhMikolYMOhTABA1Cg/SLmEg6iY27pHANfyIyPFHCKuYeqF3DM53JOoLVoXfh988AFeeukl/P777ygqKkJ5eXmzGxG1z4ELhbhYoICthRzTQ3qJHY7JY27rHv5Nhd+lQgVq6vVnqBhRZzGHiIsLuRPdmtZDPSMiIgAAEyZMaLadk5eJtBN3MAMAMD2kF+wsuYSD2Jjbuoe7vQWcrM1QUlWPc3kVGNrLUeyQiHSCOURcoX7OAIDU7BKoVAJH0RC1QuvCb9++fV0RB5FJyShUYN/ZAkgkXMJBXzC3dQ+JRIIAL3scvFCE0znlLPzIaDCHiGuQhx2szWWoqGnA+fxKDPRgl2yiG2ld+I0dO7Yr4iAyKWt+vwAAGD/QDX4uNiJHQwBzW3fy91AXflzInYwHc4i45DIpgr0dcehiEZKziln4EbWiw/1uq6qqkJ2djbq6umbbhw4d2umgiIzZ0YxibE2+AgBYOK6vyNHQjZjbul6AV+M8P3b2JGPEHCKeUF8nHLpYhJSsEjw6wlfscIj0jtaFX0FBAebOnYtff/211cc5hp2obbUNSry27QQA4KHh3po5CSQ+5rbuo27wcjq3XDP/icjQMYeIT9PZkw1eiFqldVfP559/HqWlpThy5AisrKywa9cubNy4Ef3798fPP//cFTESGY3P/riEC/mVcLE1x6uR/mKHQ9dhbus+fV1tYSaToKKmAVdKqsUOh0gnmEPEN8zXCRIJkFVUhYKKWrHDIdI7Wl/x27t3L3766SeEhoZCKpXC19cXd911F+zt7fHee+9h0qRJXREnkcG7VFCJ1fsa5/a9eV8AHKzZyVOfMLd1H3O5FP3d7JCeU470nHJ4O1uLHRJRp5lKDomNjUVsbKxeXsG0tzTDQHc7nMmtQEpWCe4J9BA7JCK9ovUVP4VCATc3NwCAk5MTCgoKAABDhgxBamqqbqMjMhKCIOD1bSdR16DCmP4uuD/IS+yQ6AbMbd3Lnwu5k5ExlRwSHR2N9PR0JCUliR1Kq/4e7lksciRE+kfrwm/gwIE4e/YsACAoKAifffYZrl69irVr18LT01PnARIZg/+mXsXhS0WwNJPi3SlDOKdJDzG3dS9/z8aOeyz8yFgwh+gHLuRO1Dath3o+99xzyMnJAQAsXboU99xzD77++muYm5tjw4YNuo6PyOAVK+rw7o50AMBzEwbApweHtekj5rbupensycKPjARziH4I9W1smnbyahlq6pWwNJOJHBGR/tC68Hvsscc0X4eEhCArKwtnzpyBj48PXFxcdBockTF4d8dplFTVY5CHHZ4c01vscKgNzG3dK6BpqOfl4mqU19TD3pJzXsmwMYfoB29nK7jYWqCwshYnrpZhOLtnE2loPdTzRtbW1hg2bBiTGlErDl0oxH9Tr0AiAZZPGwIzWaffctRNmNu6lqO1OTwdLAEAZ7iQOxkh5hBxSCSSv4d7ZnK4J9H12nXFLyYmBu+88w5sbGwQExNz031Xrlypk8CIDF1NvRKvbz8JAHhshC+G+TiJHBHdiLlNXAGe9sgpq8HpnHKE9ean8mR4mEP0U6ifE3adyuV6fkQ3aFfhd+zYMdTX1wMAUlNT22xMwYYVRH/7dN8FZBQq4GZngRfvGSh2ONQK5jZx+XvaI+FMPhu8kMFiDtFP6s6eqdklEASB33+iJu0q/Pbt26f5+vfff++qWIiMxvm8Cqz54yIA4K37B3P+kp5ibhMXG7yQoWMO0U+BXg6wkEtRrKjDpUIF+rraih0SkV7QasJRfX095HI5Tp482VXxEBk8lUrAa9tOoF4pYMIgN0RyAVm9x9wmDvVafmdzK9CgVIkcDVHHMYfoF3O5FEG9HAGAwz2JrqNV4WdmZgYfHx8olcquiofI4G1NvoykzBJYm8vwzymBHGJiAJjbxOHrbA1rcxlqG1TILFKIHQ5RhzGH6B/NQu5s8EKkoXWLwddffx2vvfYaiouLuyIeIoNWUFGL5TtPAwBi7hqAno5WIkdE7cXc1v2kUgkGeTQu5H7qGod7kmFjDtEvfy/kzp8HkZrW6/itXr0aFy5cgJeXF3x9fWFjY9Ps8dTUVJ0FR2Rolu1IR3lNAwJ72mPOKD+xwyEtMLeJw9/THqnZpTidU4EHgsWOhqjjTCWHxMbGIjY2Vu+vbqqv+F0sUKBEUQcnG3ORIyISn9aF35QpU7ogDCLD98e5AvyUdg1SCfDe1KGQc80+g8LcJg71PD929iRDZyo5JDo6GtHR0SgvL4eDg4PY4bTJ2cYcfVxtcKlAgdTsEkzwdxc7JCLRaV34LV26tCviIDJo1XVKvLH9BAAgapQfhvTS3/8MqXXMbeJgZ08yFswh+ifU1wmXChRIzmLhRwR0YI4fEbX0n73ncbm4Gp4Olvi/u7lmH1F7DfKwg0TSOD+2oKJW7HCIyIiE+joDYIMXIjWtCz+lUokPP/wQYWFh8PDwgLOzc7Mbkak5k1uOdYmXAAD/fCAQthZaX0gnPcDcJg5rczn8ejTOheJwTzJkzCH6Rz3P7/iVUtQ1cMkYIq0Lv7fffhsrV67ErFmzUFZWhpiYGEybNg1SqRRvvfVWF4RIpL9UKgGv/ngCDSoBEwe7464ADiUxVMxt4gngPD8yAswh+qevqw2crM1Q26DCqWtlYodDJDqtC7+vv/4a69atw//93/9BLpfj4YcfxhdffIElS5bgzz//7IoYifTW10eycCy7FLYWcrx9f6DY4VAnMLeJx9+zcUkHFn5kyJhD9I9EIkGIej0/LuROpH3hl5ubiyFDhgAAbG1tUVbW+AnKfffdhx07dug2OhNQWlqK0NBQBAcHIzAwEOvWrRM7JGqnvPIa/GvXWQDAixMHwsPBUuSIqDOY28TDBi9kDJhD9NMwFn5EGloXfr169UJOTg4AoG/fvti9ezcAICkpCRYWFrqNzgTY2dkhMTERaWlpOHLkCJYvX46ioiKxw6J2ePt/p1BR24Agb0c8druv2OFQJzG3iUe9pMPFAgVq6vV7bTCitjCH6Cd1g5fkrBIIgiByNETi0rrwmzp1KhISEgAAzz77LN588030798fs2fPxrx583QeoLGTyWSwtrYGANTW1kIQBCYmA5BwOg87T+RCJpXgvalDIJNKxA6JOom5TTwe9pZwtDaDUiXgQn6l2OEQdQhziH4a2ssBZjIJCipqcbm4WuxwiESldfvB999/X/P1rFmz4Ovri0OHDqF///6YPHmy1gH4+fkhKyurxfaFCxciNja2zeOuXr2Kl19+Gb/++iuqqqrQr18/xMXFITQ0VOsY2pKYmIgVK1YgJSUFOTk52LZtW4sFWmNjY7FixQrk5uYiKCgIn3zyCcLCwrR6ntLSUowdOxbnz5/HihUr4OLiorPXQLqnqG3Akp9OAQCeHN1bM0yNDJuucxu1n0QiQYCnPQ5dLEL6tXIE9uQ6mGR4mEP0k6WZDIE9HXAsuxQp2cXw6WEtdkhEotG68KupqYGl5d9zmW6//XbcfvvtHQ4gKSkJSuXfQ3tOnjyJu+66CzNmzGjzmJKSEtxxxx0YP348fv31V7i6uuL8+fNwcnJqdf+DBw8iLCwMZmZmzbanp6ejR48ecHdvvROjQqFAUFAQ5s2bh2nTprV4fMuWLYiJicHatWsxYsQIrFq1ChMnTsTZs2fh5uYGAAgODkZDQ0OLY3fv3g0vLy8AgKOjI44fP468vDxMmzYN06dPbzMmEt/He87hamk1ejlZ4bmI/mKHQzqi69ymj2JjYxEbG9ss5+oLf3Xhx3l+ZKBMIYcYqhAfJxzLLkVyZgmm3tZL7HCIRKP1UE83NzdERUVhz549UKk6vyaKq6srPDw8NLdffvkFffv2xdixY9s85oMPPoC3tzfi4uIQFhaG3r174+6770bfvn1b7KtSqRAdHY1HHnmk2R87Z8+eRXh4ODZu3Njm80RGRmLZsmWYOnVqq4+vXLkS8+fPx9y5cxEQEIC1a9fC2toaX375pWaftLQ0nDx5ssVNXfRdz93dHUFBQdi/f3+bMZG4Tl4tw5cHMwAA70wJhLU51+wzFrrObfooOjoa6enpSEpKEjuUFvy5pAMZOFPIIYYq1I8NXoiADhR+GzduRFVVFR544AH07NkTzz//PJKTk3USTF1dHTZv3ox58+ZBIml7ztTPP/+M0NBQzJgxA25ubrjtttva7IYplUqxc+dOHDt2DLNnz4ZKpcLFixcRHh6OKVOm4KWXXupwrCkpKYiIiGj2XBERETh8+HC7z5OXl4eKigoAQFlZGRITEzFw4MAW+8XGxiIgIADDhw/vULzUecqmNftUAnDfUE+MH+gmdkikQ12Z2+jW1Gv5peeUc54zGSRTySGG+PeIurPn2bwKlNfUixwNkXg61Nzl+++/R15eHpYvX4709HTcfvvtGDBgAP75z392Kpjt27ejtLQUc+bMuel+ly5dwpo1a9C/f3/89ttveOaZZ/CPf/yjzat3Xl5e2Lt3Lw4cOIBHHnkE4eHhiIiIwJo1azoca2FhIZRKZYshme7u7sjNzW33ebKysjBmzBgEBQVhzJgxePbZZzXtoK+nz5/Um4qNhzJx4moZ7CzlWDI5QOxwSMe6MrfRrfVzs4WZTIKKmgZcLWUDBjI8ppJDDPHvETc7S/g4W0MQgGPZpWKHQyQarQs/NTs7O8ydOxe7d+/GX3/9BRsbG7z99tudCmb9+vWIjIxsdRjk9VQqFYYNG4bly5fjtttuw4IFCzB//nysXbu2zWN8fHywadMmbNmyBXK5HOvXr7/pVcXuEhYWhrS0NBw/fhx//fUXnnrqKbFDolZcK63GR7sb1+x7JXIQ3Oy4Zp+x6orcRrdmLpein1vjQu7p1zjckwwXc4h+ClWv55dZLHIkROLpcOFXU1ODrVu3YsqUKRg2bBiKi4vx4osvdjiQrKwsxMfH48knn7zlvp6enggIaH7Fxd/fH9nZ2W0ek5eXhwULFmDy5MmoqqrCCy+80OFYAcDFxQUymQx5eXktnsfDw6NT5yb9s/TnU1DUKRHq64SHh/uIHQ51IV3nNmo/f8/Gwu90ToXIkRB1HHOIfgpRz/PL5jw/Ml1ad6b47bff8M0332D79u2Qy+WYPn06du/ejTvvvLNTgcTFxcHNzQ2TJk265b533HEHzp4922zbuXPn4Ovb+iLahYWFmDBhAvz9/fH999/j3LlzGDduHCwsLPDhhx92KF5zc3OEhIQgISFBs8SDSqVCQkICFi1a1KFzkn7adTIXe9LzIJdKsHzaEEi5Zp9R6qrcRu0X4GmPH3GVDV7IIDGH6LeQpit+x7JL0aBUQS7r8LUPIoOldeE3depU3Hffffjqq69w7733tlgioSNUKhXi4uIQFRUFubx5SKtXr8a2bds0i6ICwAsvvIBRo0Zh+fLlmDlzJo4ePYrPP/8cn3/+eavnjoyMhK+vr2aYZ0BAAPbs2YPw8HD07Nmzzat/lZWVuHDhguZ+RkYG0tLS4OzsDB8fH8TExCAqKgqhoaEICwvDqlWroFAoMHfu3E5/T0g/VNTU462fG9fse2psHwxwtxM5IuoqXZHbSDvXN3ghMjTMIfptgJsd7CzlqKhpwJncCq4XSiZJ68IvLy8Pdna6/eM3Pj4e2dnZmDdvXovHCgsLcfHixWbbhg8fjm3btuHVV1/FP//5T/Tu3RurVq3Co48+2uJ4qVSK5cuXY8yYMTA3N9dsDwoKQnx8PFxdXduMKzk5GePHj9fcj4mJAQBERUVhw4YNmDVrFgoKCrBkyRLk5uYiODgYu3bt4hp8RuSj3eeQW14Dvx7WeDaca/YZs67IbaQd9ZIO2cVVqKiph50l/3Amw8Ecot+kUgmG+Tjhj3MFSMkqYeFHJkkisG+2QSkvL4eDgwPKyspgb28vdjhGLe1yKaZ+ehCCAHz95Ajc0c9F7JBIzxjq+1Gf4759eQJyy2vww9MjEernLHY4RF1On9+PN2OIcf8n4TxW7jmHyUFe+OTh28QOh0gntHkvcoAzUSvqlSq8+uMJCAIw7baeLPqIukmAF4d7ElHXYGdPMnUs/IhaEXcwA6dzyuFobYbXJ/mLHQ6Ryfi7sycLPyLSrSBvR8ikElwrq8E1rhdKJoiFH9ENLhdX4eM95wEAr93rjx62FiJHRGQ6/DUNXrikAxHplo2FXPPhUkoWl3Ug06N1cxcAKC0t1XS77NevHxwdHXUZE5FoBEHAmz+dRHW9EiN6O2NGSC+xQ6JuxNwmPnVnz7O55VCqBMi4fAoZEOYQ/Rfq64yTV8uRklWCyUFeYodD1K20uuKXmZmJSZMmwcXFBSNGjMCIESPg4uKC++67D5mZmV0UIlH3+eWvHPx+tgDmMimWTxsCiYR/dJoC5jb94dvDBlZmMtTUq5BRqBA7HKJ2YQ4xHOr1/HjFj0xRu6/4Xb58GbfffjvMzMzwzjvvwN+/cd5Teno61qxZg5EjRyIpKQm9evEKCRmmsup6vP2/dADAwvF90dfVVuSIqDswt+kXmVSCQZ52OJZdivSccvRz4/uQ9BtziGFRF37pOeVQ1DbAxqJDg9+IDFK7l3N44okncOHCBfz222+wtLRs9lh1dTXuuece9O/fH1988UWXBEqNDLF9sqF4bdsJfHMkG31cbfDrc2NgIZeJHRJ1g87kNkN9P+p73Or34jPj+uLlewaJHQ7RTXX27yN9fz+2xVDjBoBR7yXgWlkNvnlyBEaxazcZuC5ZzmHXrl149913WyQ1ALCyssI777yDnTt3ah8tkR5IzizGN0eyAQDLpw5h0WdCmNv0j7rBCzt7kiFgDjE8IU1rhHK4J5madhd+hYWF8PPza/PxPn36oLiY66KQ4alrUOG1bScAADNDe+H2Pj1Ejoi6E3Ob/lE3eEm/xsKP9B9ziOEJ8XEEACSz8CMT0+7Cz9PTE+np6W0+fvLkSXh4eOgkKKLutG7/JZzLq0QPG3O8di/X7DM1zG36Z5CHHSQSIL+iFkWVtWKHQ3RTppZDYmNjERAQgOHDh4sdSoeFNl3xS80ugUrVrhlPREah3YXflClTsHjxYhQUFLR4LD8/Hy+//DKmTJmiy9iIulxmoQL/Tmhcs+/N+wLgaG0uckTU3Zjb9I+NhRy+ztYAgNNcz4/0nKnlkOjoaKSnpyMpKUnsUDpskIcdrM1lqKhpwPn8SrHDIeo27W5ltHTpUuzcuRN9+/bFY489hkGDBkEQBJw+fRrffPMNPDw8sGTJkq6MlUinBEHAG9tPoq5BhTH9XfBAMNfzMUXMbfopwMsemUVVSM8pw+j+bL5A+os5xPDIZVIEezvi0MUiJGcVY6CHndghEXWLdhd+Tk5OOHLkCF577TV89913KC0tBQA4OjrikUcewfLly+Hs7NxVcRLp3Pa0qzhwoRAWcimWTQnkmn0mirlNP/l72GPniVxe8SO9xxximEJ9nXDoYhFSMkvw6AhfscMh6hZaLV7i5OSENWvW4NNPP9UMaXB1deUfzGRwShR1eOeX0wCAf0zoD98eNiJHRGJibtM/7OxJhoQ5xPBoOntms8ELmY4OrVp54sQJnDt3DgAwcOBADBkyRKdBEXW19349jWJFHQa422L+mD5ih0N6grlNfwR4NRZ+F/IrUdug5BIrZBCYQwzHbT6OkEiArKIqFFTUwtXOQuyQiLqcVoXf0aNH8cQTTyA9PR3qdd8lEgkGDx6M9evXG3SHJzIdf14qwtbkKwCA96YNgbm83T2OyEgxt+kfTwdLOFiZoay6HufzKhHY00HskIjaxBxieOwtzTDQ3Q5nciuQklWMewI9xQ6JqMu1+y/e9PR0TJgwAVZWVti8eTNSU1ORmpqKTZs2wcLCAhMmTLhpO2MifVDboNSs2ffoCB+E+HLehaljbtNPEonk7/X8ONyT9BhziOEK8XUCwIXcyXRIBPVHU7cwc+ZMNDQ04L///W+LMeuCIGDatGkwMzPD1q1buyRQalReXg4HBweUlZXB3t5e7HAMzqr4c1gVfx6udhaIjxkLByszsUMikXUmtxnq+9FQ4v7n/9Lx5cEMzL3DD0snDxY7HKJWdfbvI0N5P97IUOO+3o+pVxCz9Thu83HEtoV3iB0OUYdo815s91DPffv24ddff211orJEIsFrr72Ge++9V/toibrJhfxKfLrvIgBg6eQAFn0EgLlNn/l7NrZYZ4MX0mfMIYYrtGnUz8mrZaipV8LSjHOJybi1e6hnRUUF3N3d23zcw8MDFRVsu036SRAEvL7tBOqUKowf6IpJQziWnxoxt+kvdYOX9GvlaOfgFKJuxxxiuLydreBqZ4F6pYATV8vEDoeoy7W78PP19cXRo0fbfPzIkSPw9eU6KKSfvk++giMZxbAyk+GfD3DNPvobc5v+6udmC7lUgvKaBlwrqxE7HKJWMYcYLolEghCfxnl+yZmc50fGr92F30MPPYSYmBicPHmyxWMnTpzA4sWLMWvWLJ0GR6QLhZW1eHdn45p9L9zVH97O1iJHRPqEuU1/Wchl6OdmCwA4fY3DPUk/MYcYtlA/dYOXYpEjIep67Z7j9+qrryI+Ph7BwcG466674O/vD0EQcPr0acTHxyMsLAyvvfZaV8ZK1CHv7jiNsup6BHjaY94dvcUOh/QMc5t+C/C0x5ncCqTnlCMioO3hdERiYQ4xbNd39hQEgSOCyKi1u/CztLTEvn378PHHH+Pbb7/FH3/8AQAYMGAAli1bhhdeeAEWFlz8kvTL/vMF2HbsKiSSxjX75DKu2UfNMbfpN39Pe+DYVTZ4Ib3FHGLYBns5wEIuRUlVPS4VKtDX1VbskIi6jFZ/BZubm+Pll19GWloaqqqqUFVVhbS0NLzyyisoKCjAggULuipOIq3V1CvxxvbGoTdRI/0Q5O0obkCkt5jb9Je6wQsLP9JnzCGGy1wuRVAvRwBACuf5kZHT2eWPoqIirF+/XlenI+q0T/aeR1ZRFTzsLfF/dw8QOxwyUMxt4vJvWsQ9s6gKlbUNIkdDpD3mEP0X4seF3Mk0cNwbGaWzuRX47I9LAIC3HxgMO0uu2UdkiJxtzOFu3zhM7mwur/oRke5pOnuywQsZORZ+ZHRUKgGvbTuBBpWAuwLcMXGwh9ghEVEnBHj+vZ4fEZGuqRu8XCxQoERRJ3I0RF2HhR8ZnW+TspGSVQIbcxnevn+w2OEQUSeph3um53ARbCLSPScbc/R1tQEApGZzuCcZr3Z39Zw2bdpNHy8tLe1sLESdll9eg/d/PQMAWDxxILwcrUSOiPQdc5v+Uxd+bPBC+og5xDiE+jrjYoECyVklmODPpWPIOLW78HNwcLjl47Nnz+50QESd8c9f0lFR04ChvRwwe6Sf2OGQAWBu03/qzp5ncsuhVAmQSbnOFukPU8shsbGxiI2NhVKpFDsUnQrxdcKW5Mvs7ElGrd2FX1xcXFfGQdRp+87m45e/ciCTSrB86hD+cUjtwtym//x62MDSTIqaehUyi7jOFukXU8sh0dHRiI6ORnl5+S2LXkOi7ux5/Eop6hpUMJdzNhQZH/5Wk1GoqmvAG9sa1+ybd4cfAnsaz39GRKZOJpVgoAeHexJR1+njYgMnazPUNqhw6lqZ2OEQdQkWfmQU/h1/HldLq9HT0Qov3MU1+4iMDTt7ElFXkkgkmu6eXM+PjBULPzJ4p66V4YsDGQCAd6YMhrV5u0cwE5GBCPC0A8ArfkTUdUJ8nQGw8CPjxcKPDJpSJeC1H09AqRIwaYgnwgexExeRMVI3eDnNJR2IqIuE+qkXci+BIAgiR0Okeyz8yKBt/jMLx6+Uwc5CjqWTA8QOh4i6iHqOX255DYq5wDIRdYEhPR1gJpOgoKIWl4urxQ6HSOdY+JHByi2rwYrfzgIAXoocBDd7S5EjIqKuYmshh28PawAc7klEXcPSTKZpDpeSXSxyNES6x8KPDNbSn0+isrYBw3wc8WiYj9jhEFEXY4MXIupqoU0NXpK5nh8ZIRZ+ZJB2n8rFb6fyIJdKsHzaEEi5Zh+R0fP35JIORNS12NmTjBkLPzI4lbUNWPrzKQDA/Dv7YFDT3B8iurnY2FgEBARg+PDhYofSIerCL52FHxF1EXVnz7N5FSivqRc5GiLdYuFHBuej3WeRU1YDH2drPDehv9jhEBmM6OhopKenIykpSexQOkTd2fNCfiVqG5QiR0NExsjVzgK+PawhCMCx7FKxwyHSKRZ+ZFD+ulKKjYcyAQDvTg2EpZlM3ICIqNt4OVjC3lKOBpWAC/mVYodDREYqxKdpuGcmG7yQcWHhRwajQanCqz+egEoApgR7YUx/V7FDIqJuJJFIrpvnx/X8iKhrhDSt55eSzXl+ZFxY+JHB2HAoE6eulcPBygxv3Mc1+4hMkXq4Jzt7ElFXCW2a53csuxQNSpXI0RDpDgs/MghXSqrw0e5zAIDX7h0EF1sLkSMiIjGwsycRdbX+braws5Sjqk6JM7kcXUDGg4Uf6T1BEPDqjydQXa9EWG9nzAz1FjskIhKJei2/07nlEARB5GiIyBhJpRIM81Gv58d5fmQ8WPiR3vv26GXsP18ISzMp3p82BBIJ1+wjMlX93W0hl0pQWlWPnLIascMhIiOlXsg9hZ09yYiw8CO9drm4Cu/uSAcAvDhxEPq42oocERGJyUIuQ9+mPMDhnkTUVTQLufOKHxkRFn6kt1QqAS//9y8o6pQY7ueEuaP8xA6JiPQAG7wQUVcL9nGETCrBtbIaXCutFjscIp1g4Ud66+uj2Th0sQiWZlKsmB4EqZRDPIkI8Pe0A9A4z4+IqCtYm8s1c4pTsrisAxkHFn4iKy0tRWhoKIKDgxEYGIh169aJHZJeuFxchfd2ngYAvHzPIPi52IgcERHpC67lR0TdQTPck4UfGQm52AGYOjs7OyQmJsLa2hoKhQKBgYGYNm0aevToIXZoolGpBLz0w1+oqmvs4hk10k/skIhIj6gLv8wiBRS1DbCx4H9lRKR7Ib5O2HAoE8lZnOdHxoFX/EQmk8lgbW0NAKitrYUgCCbfonzzkSwcvlQEKzMZPuQQTyK6gYutBdzsLCAI4BpbRNRlQv0ar/idzqmAorZB5GiIOk/0ws/Pzw8SiaTFLTo6us1j3nrrrRb7Dxo0SOexJSYmYvLkyfDy8oJEIsH27dtb7BMbGws/Pz9YWlpixIgROHr0qNbPU1paiqCgIPTq1QsvvvgiXFxcdBC9YcouqsJ7O88AAF6JHASfHtYiR0RE+ogLuRNRV/N0sIKXgyWUKgHHL5eKHQ5Rp4le+CUlJSEnJ0dz27NnDwBgxowZNz1u8ODBzY47cOBAm/sePHgQ9fX1Lbanp6cjLy+vzeMUCgWCgoIQGxvb6uNbtmxBTEwMli5ditTUVAQFBWHixInIz8/X7KOeu3fj7dq1a5p9HB0dcfz4cWRkZOCbb765aUzGTKUSsPiH46iuV+L2Ps54/HZfsUMiIj2l6ezJwo+IulCInzMAIJnz/MgIiD4xwtXVtdn9999/H3379sXYsWNvepxcLoeHh8ctz69SqRAdHY3+/fvju+++g0wmAwCcPXsW4eHhiImJwUsvvdTqsZGRkYiMjGzz3CtXrsT8+fMxd+5cAMDatWuxY8cOfPnll3jllVcAAGlpabeMUc3d3R1BQUHYv38/pk+f3u7jjMVXhzNxNKMY1uYydvEkopviFT8i6g6hvk743/FrbPBCRkH0K37Xq6urw+bNmzFv3jxIJDf/o//8+fPw8vJCnz598OijjyI7O7vV/aRSKXbu3Iljx45h9uzZUKlUuHjxIsLDwzFlypQ2i772xJqSkoKIiIhmzxUREYHDhw+3+zx5eXmoqGico1JWVobExEQMHDiwxX6xsbEICAjA8OHDOxSvvsssVOCDXWcBAK/e6w9vZw7xJKK2qdusn82tgFJl2vOiiajrqDt7pmaXQMVcQwZOrwq/7du3o7S0FHPmzLnpfiNGjMCGDRuwa9curFmzBhkZGRgzZoymgLqRl5cX9u7diwMHDuCRRx5BeHg4IiIisGbNmg7HWlhYCKVSCXd392bb3d3dkZub2+7zZGVlYcyYMQgKCsKYMWPw7LPPYsiQIS32i46ORnp6OpKSkjocs75Sd/GsrldiVN8eeDTMR+yQiEjP9XaxgaWZFFV1SmQVKcQOh4iM1CAPO1iby1BR04Bz+WwmRYZN9KGe11u/fj0iIyPh5eV10/2uH345dOhQjBgxAr6+vti6dSueeOKJVo/x8fHBpk2bMHbsWPTp0wfr16+/5VXF7hAWFqbVcFBjtOFQJo5mFsPGXIYPHhzKIZ5EdEsyqQQD3e1w/EoZTudUoI+rrdghEZERksukuM3HEQcvFCElqwSDPOzFDomow/Tmil9WVhbi4+Px5JNPan2so6MjBgwYgAsXLrS5T15eHhYsWIDJkyejqqoKL7zwQmfChYuLC2QyWYtGLHl5ee2ae0iNMgoV+NdvjV08X5vEIZ5E1H7qBi+c50dEXSnEp2kh90zO8yPDpjeFX1xcHNzc3DBp0iStj62srMTFixfh6enZ6uOFhYWYMGEC/P398eOPPyIhIQFbtmzB4sWLOxyvubk5QkJCkJCQoNmmUqmQkJCAkSNHdvi8pkSpEvDi98dRU6/C6H4ueIRDPIlIC+oGL+zsSdR9jL3nQGvY2ZOMhV4UfiqVCnFxcYiKioJc3nz06erVqzFhwoRm2xYvXow//vgDmZmZOHToEKZOnQqZTIaHH3641XNHRkbC19cXW7ZsgVwuR0BAAPbs2YO4uDh8/PHHbcZVWVmJtLQ0zVDMjIwMpKWlaRrJxMTEYN26ddi4cSNOnz6NZ555BgqFQtPlk24u7mAGkrNKYGshx/sPDtGLobdEZDjY2ZOo+xlzz4G23ObjCIkEyC6uQn5FjdjhEHWYXszxi4+PR3Z2NubNm9fiscLCQly8eLHZtitXruDhhx9GUVERXF1dMXr0aPz5558tloYAGjttLl++HGPGjIG5ublme1BQEOLj41s9Ri05ORnjx4/X3I+JiQEAREVFYcOGDZg1axYKCgqwZMkS5ObmIjg4GLt27WrR8IVaulhQiRW/NXbxfH2SP3o5cYgnEWlnkIcdACCnrAYlijo42Zjf4ggiIu3ZW5phoLsdzuRWIDWrBPcEtj7CjEjfSQRBYG9aA1JeXg4HBweUlZXB3t4wJxgrVQJmrD2E1OxSjOnvgq/mhfFqHxkkQ30/GmrcrbnzX/uQXVyFb54cgVH9XMQOh0hrhvp+NNS4O+r1bSfw9ZFsPDm6N964L0DscIg0tHkv6sVQTzIt6w9cQmp2Kews5PjgwaEs+oiowwI4z4+IukGoX1ODl2zO8yPDxcKPutWF/Ep8uPscAOCN+/zh5WglckREZMjY4IWIukOIT2ODl5NXy1BTrxQ5GqKOYeFH3UapErD4++Ooa1Bh7ABXzAz1FjskIjJw/p6N8/xO53BhZSLqOt7OVnC1s0C9UsBfV8rEDoeoQ1j4UbdZt/8S0i6Xws6SXTyJSDfUa/ldyK9AXYNK5GiIyFhJJBKE+jYN9+SyDmSgWPhRtzifV4GVexqHeL55XwA8HTjEk4g6r6ejFewt5ahXCriQXyl2OERkxEKaCr8fUi6jtKpO5GiItMfCj7pcg1KlGeI5fqArZoT0EjskIjISEokEg7ieHxF1g/uDvOBqZ4GLBQpEfXkUFTX1YodEpBUWftTlPt9/CcevlMHOUo73prGLJxHpVgALPyLqBm72lvj6yRFwsjbD8StlmLchCVV1DWKHRdRuLPyoS53Lq8CqPecBAEsnD4aHg6XIERGRseGSDkTUXQa422HTEyNgZylHUmYJ5n+VzC6fZDBY+FGX0QzxVKowYZAbHhzWU+yQiMgI+V93xU8QBJGjISJjF9jTARvnhcHGXIaDF4qw8OtUNpcig8DCj7rMZ4mX8NeVMthbyrF8Grt4ElHX6O9uC5lUgpKqeuSW14gdDhGZgGE+Tlg/Zzgs5FLsPZOP5747hgYliz/Sbyz8qEucyS3HqvjGLp5vPzAY7vYc4klEXcPSTIa+rjYAOM+PiLrP7X164PPZoTCXSfHryVy8+MNfUKk46oD0Fws/0rn6piGe9UoBEf7umBLMIZ5E1LX+bvDChdyJqPuMHeCK2EeHQS6VYNuxq3h9+wkOOSe9xcKPdG7t7xdx8mo5HKzMsHxqIId4ElGXU8/zS7/GK35E1L3uCnDHx7OCIZUA3x69jLf/l87ij/QSCz/SqdM55fjP3sYunv98YDDcOMSTiLqBP5d0ICIRTQ7ywgcPDgUAbDiUiRW/nRU5IqKWWPiRztQrVfi/rY1DPO8OcMf9QV5ih0REJkJd+GUUKbiuFhGJYkaoN96ZEggA+PT3i1jd9EE4kb5g4Uc6E7vvAtJzyuFkbYZ3p7KLJxF1H1c7C7jaWUAQgDO5nOdHROJ4/HZfvH6vPwDgw93n8MX+SyJHRPQ3Fn6kE6eulWH13gsAgLcfCISrnYXIERGRqeFwTyLSB/Pv7IOYuwYAAJbtOI1Nf2aJHBFRIxZ+1Gl1DSos/v4vNKgE3DPYA5OHeoodEhGZoAAWfkSkJ54N74dnxvUFALy5/SR+SLkickRELPxIB1bvu4DTOeVwtjHHMnbxJCKR+HvaAWBnTyISn0QiwUsTB2LOKD8AwEs/HMcvf10TNygyeSz8qFNOXi3Dp/sah3j+84HBcLHlEE8iEof6it+Z3AouokxEopNIJFg6OQAPDfeGSgCe/y4Ne9LzxA6LTBgLP+qwxiGex9GgEnDvEA/cN5RdPIlIPL1dbGAhl6KqToms4iqxwyEigkQiwbtTh2BKsBcaVAKiv05F4rkCscMiE8XCjzrsk73ncSa3Aj1szPHOA4Fih0NEJk4uk2KgR+NwT87zIyJ9IZNK8OGMIEQGeqBOqcKCTck4cqlI7LDIBLHwow45caUMn/5+EQDwzpRA9OAQTyLSA2zwQkT6SC6T4t8P3YbxA11RU6/CvA1JOJZdInZYZGJY+JHWahuU+L/v06BUCbhvqCfuHcIunkSkH9RLOrDBCxHpG3O5FGseC8Govj2gqFMi6sujOHm1TOywyISw8COt/SfhPM7lVcLF1hz/5BBPItIjXMuPiPSZpZkM62aHItTXCeU1DZj95VGcy6sQOywyESz8SCvHL5diTdMQz2VTAuFsYy5yREREfxvUtKTDtbIalFbViRwNEVFLNhZyfDl3OIb2ckCxog6PfnEEGYUKscMiE8DCj9qtpl6Jxd8fh0oA7g/ywj2BHOJJRPrF3tIM3s5WAIB0XvUjIj1lb2mGr+aFYZCHHQoqavHouj9xpYTdiKlrsfCjdvt3wnmcz6+Ei60F3r5/sNjhEBG1yt9DPdyTw6eISH85Wptj85Mj0MfVBtfKavDIuiPILasROywyYiz8qF2OZZfgsz8ah3gunxoIJw7xJCI9FeDFeX5EZBhcbC3wzZO3w8fZGtnFVXj0iz9RWFkrdlhkpFj40S1dP8RzSrAX7h7sIXZIRERtYmdPIjIkHg6W+PrJEfB0sMTFAgUe++II5yhTl2DhR7f08Z5zuFiggKudBd7iEE8i0nPqtfwu5FeirkElcjRERLfm7WyNr58cARdbC5zJrUDUl0dRUVMvdlhkZFj40U2lZpdg3f5LAIDlU4fA0ZpDPIlIv/VysoKdpRx1ShUuFlSKHQ4RUbv0cbXF10+OgJO1GY5fKcO8DUmoqmsQOywyIiz8qE3XD/GcdltP3BXgLnZIRES3JJFIrmvwwuGeRGQ4BnrYYdMTI2BnKUdSZgnmf5WMmnql2GGRkWDhR236aPdZXCpQwM3OAksnc4gnERkONnghIkMV2NMBG+eFwcZchoMXirDw61QOWyedYOFHrUrJKsYXBzIAAO8/OAQO1mYiR0RE1H7+TQu5cy0/IjJEw3ycsH7OcFjIpdh7Jh/PbzmGBiWLP+ocFn7UQnWdEou//wuCADw4rBfCB3GIJxEZFnVnz9M5FRAEQeRoiIi0d3ufHvh8dijMZVLsPJGLF3/4CyoV8xl1HAs/auHD3WeRUaiAu70FlkwOEDscIiKtDXC3g0wqQbGiDvkVXBOLiAzT2AGuWP3IbZBJJdh27Cpe336SH2ZRh7Hwo2aSMovx5UH1EM+hcLDiEE8iMjyWZjL0cbEBwPX8iMiw3T3YA6tmBUMqAb49mo1//pLO4o86hIUfaVTXKfHi98chCMDM0F4YP9BN7JCIiDpMs5A75/kRkYGbHOSFDx4cCgCIO5iJD3efFTkiMkQs/EjjX7+dQWZRFTwdLPHGfRziSUSGjZ09iciYzAj1xjsPNHZZj913Eav3nhc5IjI0LPwIAHDkUhHiDmYCaBziaW/JIZ5EZNh4xY+IjM3jI/3w+r3+AIAPd5/DF/sviRwRGRIWfoSquga8+MNfAICHhntj7ABXkSMiIuq8gKbCL6NQgaq6BpGjIdIvly9fxrhx4xAQEIChQ4fi+++/Fzskaqf5d/ZBzF0DAADLdpzG5j+zRI6IDAULP8K/dp1FdnEVvBws8fokf7HDISLSCVc7C7jYWkAQgLO5FWKHQ6RX5HI5Vq1ahfT0dOzevRvPP/88FAqF2GFROz0b3g9Pj+0LAHhj+0n8kHJF5IjIELDwM3GHLxZhw6FMAI1DPO04xJOIjIh6IffTOSz8iK7n6emJ4OBgAICHhwdcXFxQXFwsblDUbhKJBC/fMxBzRvkBAF764Th++euauEGR3mPhZ8IUtQ146b/HAQAPh/ngTg7xJCIjwwYvZKgSExMxefJkeHl5QSKRYPv27S32iY2NhZ+fHywtLTFixAgcPXq0Q8+VkpICpVIJb2/vTkZN3UkikWDJfQF4aLg3VALw/HdpiE/PEzss0mMs/EzY+7+eweXiavR0tMJr9w4SOxwiIp0LYIMXMlAKhQJBQUGIjY1t9fEtW7YgJiYGS5cuRWpqKoKCgjBx4kTk5+dr9gkODkZgYGCL27Vrf18ZKi4uxuzZs/H55593+Wsi3ZNKJXh36hBMCfZCg0rAwq9Tsf98gdhhkZ6Six0AiePQhUJsapoM/AGHeBKRkVJ39jyTUw6VSoBUKhE5IqL2iYyMRGRkZJuPr1y5EvPnz8fcuXMBAGvXrsWOHTvw5Zdf4pVXXgEApKWl3fQ5amtrMWXKFLzyyisYNWrULfetra3V3C8v54cp+kImleDDGUGoqVdh16lczP8qGRvnhmFEnx5ih0Z6hlf8TFBlbQNe+m9jF89HR/hgdH8XkSMiIuoafVxsYC6XQlGnxOWSKrHDIdKJuro6pKSkICIiQrNNKpUiIiIChw8fbtc5BEHAnDlzEB4ejscff/yW+7/33ntwcHDQ3DgsVL/IZVL85+HbMH6gK2rqVZi3IQnHskvEDov0DAs/E/TeztO4UtI4xPPVe9nFk4iMl1wmxUD3xgYv6dd4hYKMQ2FhIZRKJdzd3Zttd3d3R25ubrvOcfDgQWzZsgXbt29HcHAwgoODceLEiTb3f/XVV1FWVqa5Xb58uVOvgXTPXC7FmsdCMKpvDyjqlIj68ihOXSsTOyzSIxzqaWIOnC/E10eyAQArpg+FrQV/BYjIuPl72uHE1TKczilH5BBPscMh0gujR4+GSqVq9/4WFhawsLDowohIFyzNZFg3OxSzvzyKlKwSPL7+KLYsuB39mz4AI9PGK34mpKKmHi83DfF8/HZfjOrHIZ5EZPz+bvDCJR3IOLi4uEAmkyEvr3kHx7y8PHh4eIgUFekLGws54uYOx5CeDihW1OGRL44go5BrNBKv+JmU5TvP4GppNbydrfBKJLt4EpFpUDd44ZIOZCzMzc0REhKChIQETJkyBQCgUqmQkJCARYsWiRsc6QV7SzN8NS8MD6/7E2dyK/Douj+x9emR6OVkLXZohMY5tpW1DShW1KGwsg5FlbUoUjT+W1hZp/m6qLIO3y64Hc425jp5XhZ+IistLUVERAQaGhrQ0NCA5557DvPnz9f58ySeK8C3RxuHeP7rwSDYcIgnEZkI/6a1/K6WVqOsqh4O1uxiTPqvsrISFy5c0NzPyMhAWloanJ2d4ePjg5iYGERFRSE0NBRhYWFYtWoVFAqFpssnkZONOTY9MQKzPj+MSwUKPLLuCL5/eiTc7S3FDs0o1dQrUayoQ1FlHQoVtSiurEORorF4K7zu66LKWhQq6lDX0L6h1kWVtSz8jIWdnR0SExNhbW0NhUKBwMBATJs2DT166K4Fb3lNPV5pGuIZNdIXI/uyvS8RmQ57SzP0crLClZJqpOeUMweSQUhOTsb48eM192NiYgAAUVFR2LBhA2bNmoWCggIsWbIEubm5CA4Oxq5du1o0fCHT5mpngW+evB0zPjuE7OIqPLLuT2x5aiRcbDlf81aUKgElVXXNijX1VThNEXfdtoraBq2fw9pchh625uhhYwGXpn972Jqjh+3f970crXT2mlj4iUwmk8HauvGye21tLQRBgCAInT+xSglkHQIq8/BdqgK5Zc7wcbbFyxziSUQmKMDDBr3KUlB7LAuQDgF8RwFSmdhhGRRBEKBUCWhQ35Sqpn8FNKhUTf82/1qpUqFe2XhcvVLV9K/6PKq/j1U1bgMAqUQCqUQCmfT6ryWQSiWQSgCZRP11K/tImvaRSiBp2iaTSCBp2iZrOkez/Zv2kUpw3dcSSKVo87zdYdy4cbf8e2DRokUc2km35OFgiW+evB0zPzuMiwUKPL7+KL59IhSOBclAZR5g624SOVEQBJTXNDRdlatt9SqcuogrVtShuKoO2v5JLpdKNIVcD1tzuNhaoIdNYyHXeP+64s7GAlbm3fs9F73w8/PzQ1ZWVovtCxcuRGxs7C2Pf//99/Hqq6/iueeew6pVq3QaW2JiIlasWIGUlBTk5ORg27ZtmrH0arGxsVixYgVyc3MRFBSETz75BGFhYVo9T2lpKcaOHYvz589jxYoVcHHpZNOV9J+BXS8D5dcAAAsA3GfhjMrhy2BtLvqPnIioe6X/jI+u/h/szPOBk2i82XsB93wABNwvdnQagiCgtkGF6jolquqVqK5rQHWdClV1DU33laiqa9pe3/h1XcMNxdcNBVmzguuGQqxF4Xbd19cXaOrjG1Q6+FDSCEjUheP1xaHkusJUUyhKMDnIE69PChA7ZJ2KjY1FbGwslEql2KGQFrydrfH1kyMw87M/4ZMXj4aVcwBV4d876GFOvJFSJaC6KRfW1Cs1X1c3fV1Tp0RF7d+FXeOQy+ZX6eqV2uUxiQRwsjaHs405etg0FXLNCrumoq6puLO3lHfbh0MdIXoVkJSU1Cx5nDx5EnfddRdmzJjRrmM/++wzDB069Kb7HTx4EGFhYTAzaz6vIz09HT169GhzWIRCoUBQUBDmzZuHadOmtXh8y5YtiImJwdq1azFixAisWrUKEydOxNmzZ+Hm5gYACA4ORkNDy0u/u3fvhpeXFwDA0dERx48fR15eHqZNm4bp06d3fKhG+s/A1tkAmv9ie0qKIfkjGnC30+s3NRHd2uXLl/H4448jPz8fcrkcb775Zrtypklqyom2N+RElOc05sqZX2mVExuUKlQ1/YFRpS7E6hv+/rrpjxB1gdbq9qb9q687R029ElV1DTDE2komlUCuvsmkTf9KIJdKIZc1XjEzk0ob/226//d+Us2xMmnjH0sqQYBS1fivqukqo0oQoFIBSkGAqum+UsDfX6v3uX5b0zE3Pq5UCa3vIwjt+nRfEBrjUEIAblH7lFdrP/RL30VHRyM6Ohrl5eVwcHAQOxzSQh9XW/wUXgTP3asaf3evr086mBOBxg+t6pQq1NSpmvJc44dTNfVKVDdtUxdm1dcVbDcWb83vqzQfcjXuq0Kdsv3Lj9yMrYW8qXgzh7N6iOWNV+ma7jtZm0EuM55FECSCTsYV6s7zzz+PX375BefPn79pxVxZWYlhw4bh008/xbJlyxAcHNzqFT+VSoVhw4ahf//++O677yCTNV5SPXv2LMaOHYuYmBi89NJLt4xLIpG0uOI3YsQIDB8+HKtXr9Y8l7e3N5599lm88sor2r3wJgsXLkR4eDimT5/e6uPqRFtWVgZ7e/sbXqwSWBWoudLXyqto/ETn+RNGfzmfqDvc9P3YhXJycpCXl4fg4GDk5uYiJCQE586dg42NTbuOFyvubneLnChAggpzN3wy5L+oasDfhVj930VbdbMrbUqd/eFxK+YyKazMZbA2l/39r5kMVuZyWJv9vd1CLruuuJJA1lRstSzC/i7E5DcUYWay6+/fULjdUMSpCzT1MfJuHPrYHYTri8NmhWdjQakU2ihC1fvcUKg6WpvB2/nmXRQN9f1oqHGbtKacKJRfQ2vvWgESVFq44bPbtqGqXnJd8XbdVbUbC7Wmr7v7Q6vGfNiYFy3NpJqvrc3lTVff/r4Spynkmu5bmhnX38DavBdFv+J3vbq6OmzevBkxMTG3/I8kOjoakyZNQkREBJYtW9bmflKpFDt37sSdd96J2bNnY9OmTcjIyEB4eDimTJnSrqKvrVhTUlLw6quvNnuuiIgIHD58uN3nycvLg7W1Nezs7FBWVobExEQ888wzLfZr19CKrEM3KfoAQADKrzbu13tMu2MkIv3i6ekJT8/Ghcg9PDzg4uKC4uLidhd+JuMWOVECAfZ1eThxeBf+VGk3HE8qAazN5TcUZeqv5bC+rmiz0hRq122/xf7G9AmzIZFIJJA1DdckMjpNObGt324JBNjV5iE5cafWOVHNTCaBpZnshsKslfvm0sZtZjJYqnNi0z437m91w+MWcqlRfeDUnfSq8Nu+fTtKS0sxZ86cm+733XffITU1FUlJSe06r5eXF/bu3YsxY8bgkUceweHDhxEREYE1a9Z0ONbCwkIolcoWQzLd3d1x5syZdp8nKysLCxYs0DR1efbZZzFkyJAW+7VraEVlXuvbO7ofEXVId80PBoCUlBQolUp4e3vrKHoj0s5c97C/OW737P93cdZUqFmay5q+/rvAs276o4R/eBCRwWlnTpw+QI7bPPo2K8ysWxRususel2q2m/FDK72mV4Xf+vXrERkZqZn71prLly/jueeew549e2Bp2f51SHx8fLBp0yaMHTsWffr0wfr16/XiP+2wsDCkpaXp5mS27ZwX2N79iKhDumt+cHFxMWbPno1169Z17QsyVO3MdQ+MHgb0HtDFwRARiaydOXH62FCgN7vAGyO9KfyysrIQHx+PH3/88ab7paSkID8/H8OGDdNsUyqVSExMxOrVq1FbW6uZx3e9vLw8LFiwAJMnT0ZSUhJeeOEFfPLJJx2O18XFBTKZDHl5zT89ycvLg4eHR4fP2ym+oxrn8JXn4MbmLo2a5vj5juruyIhMSmRkJCIjI9t8fOXKlZg/f75moeW1a9dix44d+PLLLzXzg2/1gVBtbS2mTJmCV155BaNG3fw9XVtbi9raWs398vLydr4SA8ecSET0N+ZEk6c312Pj4uLg5uaGSZMm3XS/CRMm4MSJE0hLS9PcQkND8eijjyItLa3Voq+wsBATJkyAv78/fvzxRyQkJGDLli1YvHhxh+M1NzdHSEgIEhISNNtUKhUSEhIwcuTIDp+3U6Syxla8ANBiBHfT/XveZ2MXIhGp5wdHRERotmk7P1gQBMyZMwfh4eF4/PHHb7n/e++9BwcHB83NZIaFMicSEf2NOdHk6UXhp1KpEBcXh6ioKMjlzS9Crl69GhMmTNDct7OzQ2BgYLObjY0NevTogcDAwFbPHRkZCV9fX2zZsgVyuRwBAQHYs2cP4uLi8PHHH7cZV2Vlpaa4BICMjAykpaUhOzsbABATE4N169Zh48aNOH36NJ555hkoFArNp/iiCLi/sRWvvWfz7fZeHWrRS0S6dbP5wbm5ue06x8GDB7FlyxZs374dwcHBCA4OxokTJ9rc/9VXX0VZWZnmdvny5U69BoPCnEhE9DfmRJOmF0M94+PjkZ2djXnz5rV4rLCwEBcvXuzwuaVSKZYvX44xY8bA3Nxcsz0oKAjx8fFwdXVt89jk5GSMHz9ecz8mJgYAEBUVhQ0bNmDWrFkoKCjAkiVLkJubi+DgYOzatavja/DpSsD9wKBJjd2bKvMax3T7juInOERGYvTo0VCp2r+sgIWFBSwsLLowIj3HnEikM1zA3QgwJ5osvVvHj26O6+YQ6Y/2vB9vXAO0rq4O1tbW+OGHH5p1+oyKikJpaSl++uknvYibiLqHob4fDTVuImOjzXtRL4Z6EhGZCr2cH0xERERGTy+GehIRGZPKykpcuHBBc189P9jZ2Rk+Pj6IiYlBVFQUQkNDERYWhlWrVok/P5iIiIiMGgs/IiIdM9j5wURERGS0WPgREenYuHHjcKvp04sWLcKiRYu6KSIiIiIydZzjR0REREREZORY+BERERERERk5DvU0MOrhY+Xl5SJHQkTq96GhrYrDPEKkP5hHiKgztMkhLPwMTEVFBQDA29tb5EiISK2iogIODg5ih9FuzCNE+od5hIg6oz05hAu4GxiVSoVr167Bzs4OEonkpvuWl5fD29sbly9fNorFVY3p9RjTawFM9/UIgoCKigp4eXlBKjWckfPtzSOm+nM1FHw9+o15pJGp/lwNgTG9FsB0X482OYRX/AyMVCpFr169tDrG3t7eKN4Aasb0eozptQCm+XoM6RN6NW3ziCn+XA0JX49+Yx5pZIo/V0NhTK8FMM3X094cYjgfLREREREREVGHsPAjIiIiIiIyciz8jJiFhQWWLl0KCwsLsUPRCWN6Pcb0WgC+HmNlbN8Hvh79xtdjnIzt+2BMr8eYXgvA19MebO5CRERERERk5HjFj4iIiIiIyMix8CMiIiIiIjJyLPyIiIiIiIiMHAs/IiIiIiIiI8fCz0jFxsbCz88PlpaWGDFiBI4ePSp2SB2WmJiIyZMnw8vLCxKJBNu3bxc7pA577733MHz4cNjZ2cHNzQ1TpkzB2bNnxQ6rw9asWYOhQ4dqFhcdOXIkfv31V7HD0on3338fEokEzz//vNihiMZY8ogx5RDAuPKIMecQgHnEWHIIYFx5xJhyCGDceUTXOYSFnxHasmULYmJisHTpUqSmpiIoKAgTJ05Efn6+2KF1iEKhQFBQEGJjY8UOpdP++OMPREdH488//8SePXtQX1+Pu+++GwqFQuzQOqRXr154//33kZKSguTkZISHh+OBBx7AqVOnxA6tU5KSkvDZZ59h6NChYociGmPKI8aUQwDjyiPGmkMA5hFjyiGAceURY8ohgPHmkS7JIQIZnbCwMCE6OlpzX6lUCl5eXsJ7770nYlS6AUDYtm2b2GHoTH5+vgBA+OOPP8QORWecnJyEL774QuwwOqyiokLo37+/sGfPHmHs2LHCc889J3ZIojDWPGJsOUQQjC+PGHoOEQTmEUEw3hwiCMaXR4wthwiC4eeRrsohvOJnZOrq6pCSkoKIiAjNNqlUioiICBw+fFjEyKg1ZWVlAABnZ2eRI+k8pVKJ7777DgqFAiNHjhQ7nA6Ljo7GpEmTmr2HTA3ziGExljxiLDkEYB5hDjEsxpJDAOPJI12VQ+Q6PRuJrrCwEEqlEu7u7s22u7u748yZMyJFRa1RqVR4/vnncccddyAwMFDscDrsxIkTGDlyJGpqamBra4tt27YhICBA7LA65LvvvkNqaiqSkpLEDkVUzCOGwxjyiDHlEIB5BGAOMSTGkEMA48ojXZlDWPgRiSQ6OhonT57EgQMHxA6lUwYOHIi0tDSUlZXhhx9+QFRUFP744w+DS7iXL1/Gc889hz179sDS0lLscIjaxRjyiLHkEIB5hAyPMeQQwHjySFfnEBZ+RsbFxQUymQx5eXnNtufl5cHDw0OkqOhGixYtwi+//ILExET06tVL7HA6xdzcHP369QMAhISEICkpCf/+97/x2WefiRyZdlJSUpCfn49hw4ZptimVSiQmJmL16tWora2FTCYTMcLuwzxiGIwljxhLDgGYR9SYQwyDseQQwHjySFfnEM7xMzLm5uYICQlBQkKCZptKpUJCQoJBj3U2FoIgYNGiRdi2bRv27t2L3r17ix2SzqlUKtTW1oodhtYmTJiAEydOIC0tTXMLDQ3Fo48+irS0NJP4Y02NeUS/GXseMdQcAjCPqDGH6DdjzyGA4eaRrs4hvOJnhGJiYhAVFYXQ0FCEhYVh1apV/9/OnYZE1fZxHP+NPk2aI2VTmBUWJEmBxrQYtlBk2EKBLyqoUDMpWpSkErKCgoqI6lVFG2G9iSgIpLJs016IpW0uZbZQUdFCGzRTGY7X/aLu8zTY/WBOPXqfvh8Y8FzXdc75zwF/8nfOGfl8PmVlZbV3aW3i9Xp1//59a/vhw4e6efOmunfvrtjY2Has7OctXbpUhw8fVlFRkSIjI/XixQtJUteuXRUeHt7O1f28goICTZkyRbGxsfrw4YMOHz6ssrIylZSUtHdpPy0yMrLF8w0RERFyu93/6uce2spOOWKnDJHslSN2yhCJHPmenTJEsleO2ClDJHvlyG/PkF/y3aDocHbs2GFiY2ON0+k0SUlJ5vLly+1dUpuVlpYaSS1emZmZ7V3aT/vR+5BkCgsL27u0Npk/f77p16+fcTqdpmfPniYlJcWcPXu2vcv6Zf7Ur2H/m11yxE4ZYoy9csTuGWLMn50jdskQY+yVI3bKEGPsnyO/MkMcxhgTfPsIAAAAAOioeMYPAAAAAGyOxg8AAAAAbI7GDwAAAABsjsYPAAAAAGyOxg8AAAAAbI7GDwAAAABsjsYPAAAAAGyOxg8I0vjx45WXl9fq9QcPHlS3bt1+Wz3fe/TokRwOh27evPl/OR+AtiFHAASDDEFr0PgBAAAAgM3R+AE28OXLl3/lsQF0HOQIgGCQIR0fjR9safz48crNzVVeXp6ioqIUHR2t/fv3y+fzKSsrS5GRkYqLi9Pp06cD9rt06ZKSkpLUuXNnxcTEaNWqVWpqarLmfT6fMjIy5HK5FBMTo+3bt7c4d2Njo1auXKk+ffooIiJCI0eOVFlZ2U/VX1tbqwkTJig8PFxut1sLFy6U1+u15ufNm6e0tDRt2rRJvXv3Vnx8vCSpsrJSHo9HYWFhGj58uG7cuNHi2HV1dZoyZYpcLpeio6OVnp6u169fB1y7nJwc5eXlqUePHpo0adJP1Q7YBTlCjgDBIEPIkI6Gxg+2dejQIfXo0UOVlZXKzc3V4sWLNXPmTI0aNUrXr19Xamqq0tPT9fHjR0nSs2fPNHXqVI0YMULV1dXavXu3Dhw4oI0bN1rHzM/P16VLl1RUVKSzZ8+qrKxM169fDzhvTk6OKioqdOTIEdXU1GjmzJmaPHmy7t2716q6fT6fJk2apKioKFVVVenYsWM6f/68cnJyAtZduHBBDQ0NOnfunE6ePCmv16tp06Zp8ODBunbtmtavX6+VK1cG7PP+/XtNmDBBHo9HV69e1ZkzZ/Ty5UvNmjWrxbVzOp0qLy/Xnj17Wn3NAbshR8gRIBhkCBnSoRjAhsaNG2fGjBljbTc1NZmIiAiTnp5ujT1//txIMhUVFcYYY1avXm3i4+NNc3OztWbXrl3G5XIZv99vPnz4YJxOpzl69Kg1/+bNGxMeHm6WLVtmjDHm8ePHJjQ01Dx79iygnpSUFFNQUGCMMaawsNB07dr1H2vft2+fiYqKMl6v1xo7deqUCQkJMS9evDDGGJOZmWmio6NNY2OjtWbv3r3G7XabT58+WWO7d+82ksyNGzeMMcZs2LDBpKamBpzvyZMnRpJpaGiwrp3H4/nH+oA/BTnyFTkCtA0Z8hUZ0nH8p906TuA3S0xMtH4ODQ2V2+1WQkKCNRYdHS1JevXqlSSpvr5eycnJcjgc1prRo0fL6/Xq6dOnevfunb58+aKRI0da8927d7dubZC+3hbh9/s1cODAgFoaGxvldrtbVXd9fb2GDBmiiIiIgDqam5vV0NBg1Z2QkCCn0xmwX2JiosLCwqyx5OTkgGNXV1ertLRULperxXkfPHhg1T1s2LBW1QrYHTlCjgDBIEPIkI6Exg+21alTp4Bth8MRMPZ3qDY3N/+yc3q9XoWGhuratWsKDQ0NmPtRwAXj+zBuLa/Xq+nTp2vLli0t5mJiYoI6NmBH5EhL5AjQemRIS2RI++EZP+CbQYMGqaKiQsYYa6y8vFyRkZHq27evBgwYoE6dOunKlSvW/Lt373T37l1r2+PxyO/369WrV4qLiwt49erVq9V1VFdXy+fzBdQREhIS8B+9H+1XU1Ojz58/W2OXL18OWDN06FDdunVL/fv3b1EfAQsEjxwhR4BgkCFkyO9E4wd8s2TJEj158kS5ubm6c+eOioqKtG7dOi1fvlwhISFyuVzKzs5Wfn6+Ll68qLq6Os2bN08hIf/9NRo4cKDmzp2rjIwMHT9+XA8fPlRlZaU2b96sU6dOtaqOuXPnKiwsTJmZmaqrq1Npaalyc3OVnp5u3VrxI3PmzJHD4dCCBQt0+/ZtFRcXa9u2bQFrli5dqrdv32r27NmqqqrSgwcPVFJSoqysLPn9/rZdOAAWcoQcAYJBhpAhvxONH/BNnz59VFxcrMrKSg0ZMkSLFi1Sdna21q5da63ZunWrxo4dq+nTp2vixIkaM2ZMi3vQCwsLlZGRoRUrVig+Pl5paWmqqqpSbGxsq+ro0qWLSkpK9PbtW40YMUIzZsxQSkqKdu7c+T/3c7lcOnHihGpra+XxeLRmzZoWt1H07t1b5eXl8vv9Sk1NVUJCgvLy8tStW7eAPxoA2oYcIUeAYJAhZMjv5DDff5YMAAAAALAd2moAAAAAsDkaPwAAAACwORo/AAAAALA5Gj8AAAAAsDkaPwAAAACwORo/AAAAALA5Gj8AAAAAsDkaPwAAAACwORo/AAAAALA5Gj8AAAAAsDkaPwAAAACwORo/AAAAALC5vwAx8xgyrc25eAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(1, len(data), figsize=(10, 4))\n", + "for true_poly_order, (axi, (poly_orders, variances)) in enumerate(zip(ax, data)):\n", + " imin = np.argmin(variances)\n", + " plt.sca(axi)\n", + " plt.title(f\"true polynomial order {true_poly_order}\\nselected {imin}\")\n", + " plt.plot(poly_orders, variances)\n", + " plt.plot(poly_orders[imin], variances[imin], marker=\"o\")\n", + " plt.semilogy()\n", + " plt.ylabel(\"LOO variance\")\n", + " plt.xlabel(\"model order\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can slightly simplify this by using `cross_validation` from the library." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "data = []\n", + "for poly_order, y in enumerate(y_set):\n", + " variances = []\n", + " poly_orders = np.arange(5)\n", + " for poly_order in poly_orders:\n", + " variances.append(cross_validation(predict, x, y, poly_order))\n", + " data.append((poly_orders, variances))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, "metadata": {}, "outputs": [ { diff --git a/pyproject.toml b/pyproject.toml index 6fb4f99..1d39d18 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -50,9 +50,12 @@ profile = "black" multi_line_output = 3 [tool.ruff] -select = ["E", "F", "W", "D"] -ignore = ["D212", "D203"] -unfixable = ["ERA"] +extend-select = ["E", "F", "W", "D"] +lint.ignore = ["D212", "D203"] +lint.unfixable = ["ERA"] + +[tool.ruff.lint.per-file-ignores] +"test_*.py" = ["D"] [tool.mypy] strict = true diff --git a/src/resample/jackknife.py b/src/resample/jackknife.py index 8280bac..555ec13 100644 --- a/src/resample/jackknife.py +++ b/src/resample/jackknife.py @@ -18,7 +18,7 @@ from typing import Any, Callable, Collection, Generator, List import numpy as np -from numpy.typing import ArrayLike +from numpy.typing import ArrayLike, floating def resample( @@ -319,3 +319,40 @@ def variance( thetas = jackknife(fn, sample, *args) n = len(sample) return (n - 1) * np.var(thetas, ddof=0, axis=0) + + +def cross_validation( + predict: Callable[..., float], x: "ArrayLike", y: "ArrayLike", *args: "ArrayLike" +) -> floating[Any]: + """ + Calculate mean-squared error of model with leave-one-out-cross-validation. + + Wikipedia: + https://en.wikipedia.org/wiki/Cross-validation_(statistics) + + Parameters + ---------- + predict : callable + Function with the signature (x_in, y_in, x_out, *args). It takes x_in, y_in, + which are arrays with the same length. x_out should be one element of the x + array. *args are further optional arguments for the function. The function + should return the prediction y(x_out). + x : array-like + Explanatory variable. Must be an array of shape (N, ...), where N is the number + of samples. + y : array-like + Observations. Must be an array of shape (N, ...). + *args: + Optional arguments which are passed unmodified to predict. + + Returns + ------- + float + Variance of the difference (y[i] - predict(..., x[i], *args)). + + """ + deltas = [] + for i, (x_in, y_in) in enumerate(resample(x, y, copy=False)): + yip = predict(x_in, y_in, x[i], *args) + deltas.append((y[i] - yip)) + return np.var(deltas) diff --git a/tests/test_jackknife.py b/tests/test_jackknife.py index f553727..5c65bb5 100644 --- a/tests/test_jackknife.py +++ b/tests/test_jackknife.py @@ -1,8 +1,15 @@ import numpy as np import pytest from numpy.testing import assert_almost_equal, assert_equal - -from resample.jackknife import bias, bias_corrected, jackknife, resample, variance +from scipy.optimize import curve_fit +from resample.jackknife import ( + bias, + bias_corrected, + jackknife, + resample, + variance, + cross_validation, +) def test_resample_1d(): @@ -120,3 +127,22 @@ def test_resample_deprecation(): with pytest.warns(VisibleDeprecationWarning): with pytest.raises(ValueError): # too many arguments resample(data, True, 1) + + +@pytest.mark.filterwarnings("ignore:Covariance") +def test_cross_validation(): + x = [1, 2, 3] + y = [3, 4, 5] + + def predict(xi, yi, xo, npar): + def model(x, *par): + return np.polyval(par, x) + + popt = curve_fit(model, xi, yi, p0=np.zeros(npar))[0] + return model(xo, *popt) + + v = cross_validation(predict, x, y, 2) + assert v == pytest.approx(0) + + v2 = cross_validation(predict, x, y, 1) + assert v2 == pytest.approx(1.5)