{ "cells": [ { "cell_type": "markdown", "id": "820dffa1-2d3f-4561-b2fc-0dc63e30b401", "metadata": {}, "source": [ "# Minimum-entropy fit of the isochrone potential\n", "\n", "Here we fit a gravitational potential of the isochrone model using the kinematics of a sample in equilibrium in the correct potential, but without assuming any functional form for the tracers' distribution function (DF). It only assumes that the DF is a function of actions evaluated in each trial potential. The correct potential is recovered as the one with minimum entropy. The physical basis of the method is explained in [Beraldo e Silva et al. (2025)](https://ui.adsabs.harvard.edu/abs/2024arXiv240707947B/abstract).\n", "\n", "This code uses the [agama](https://github.com/GalacticDynamics-Oxford/Agama) package to generate an equilibrium sample and to estimate the potential and action-variables, but any other library, such as [gala](http://gala.adrian.pw/en/latest/) or [galpy](https://www.galpy.org/), or personal code could be used for these tasks instead.\n", "\n", "It fits the potential based on the original dataset and then later fits for a number of bootstrap samples. At the end, it shows the results based on the single (original) dataset, on each bootstrap sample and on averaging over all datasets. The plot gives an indication of how the log-likelihood surface would look like around the best if the correct functional form of the DF were specified (and then one would have a bona fide log-likelihood, which is not the case here since no DF is assumed). " ] }, { "cell_type": "code", "execution_count": 1, "id": "e080ef56-806b-4585-b36c-603c35ed2047", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import optimize\n", "import matplotlib.pyplot as plt\n", "plt.rcParams['figure.dpi'] = 60\n", "import tropygal\n", "import agama\n", "\n", "from matplotlib.ticker import ScalarFormatter, NullFormatter\n", "\n", "params = {'axes.labelsize': 24, \n", " 'xtick.labelsize': 20, \n", " 'xtick.direction': 'in',\n", " 'xtick.major.size': 8.0,\n", " 'xtick.bottom': 1,\n", " 'xtick.top': 1,\n", " 'ytick.labelsize': 20, \n", " 'ytick.direction': 'in',\n", " 'ytick.major.size': 8.0,\n", " 'ytick.left': 1,\n", " 'ytick.right': 1,\n", "# 'xtick.which': 'both',\n", " 'text.usetex': True, \n", " 'lines.linewidth': 1,\n", " 'axes.titlesize': 32,\n", " 'font.family': 'serif'}\n", "plt.rcParams.update(params)\n", "columnwidth = 240./72.27\n", "textwidth = 504.0/72.27\n", "\n", "%matplotlib inline\n", "# %matplotlib qt\n", "# %matplotlib widget\n", "# %matplotlib notebook\n", "%config InlineBackend.figure_format = 'retina'" ] }, { "cell_type": "markdown", "id": "8b0af6b0-c0a7-4fad-9f17-991e5252c5ca", "metadata": {}, "source": [ "## Function to return entropy" ] }, { "cell_type": "code", "execution_count": 2, "id": "a81b29e6-c2c9-45e4-bf13-d99911fc8c3c", "metadata": {}, "outputs": [], "source": [ "def func_S_J(params, pos, vel, k=1):\n", " \"\"\" \n", " Returns the entropy of the sample in each trial potential.\n", " It uses agama to estimate the potential and action-variables\n", " \n", " Parameters\n", " ---------- \n", " params: array\n", " parameters of the trial potential\n", " pos: array\n", " (x, y, z) spatial coordinates\n", " vel: array\n", " (vx, vy, vz) velocity coordinates\n", " k: int\n", " kNN, which neighbor to take\n", "\n", " Returns\n", " -------\n", " float number\n", " The entropy of the sample assuming the DF is a function of actions\n", " \n", " \"\"\"\n", "\n", " M, b = params\n", " prior = (M > 0) & (b>0)\n", " if (prior == False):\n", " return np.inf\n", "\n", " pot = agama.Potential(type='Isochrone', mass=M, scaleRadius=b)\n", " v2 = vel[:,0]**2 + vel[:,1]**2 + vel[:,2]**2\n", " E = pot.potential(pos) + 0.5*v2\n", "\n", " # Only consider models where all particles are bound\n", " if np.all(E<0):\n", " actF = agama.ActionFinder(pot) # action finder\n", " coords = np.column_stack([pos, vel])\n", " # calc actions: \n", " J = actF(coords, actions=True, frequencies=False, angles=False)\n", " sigma_J = np.array([0.5*(np.percentile(action, 84) - np.percentile(action, 16)) for action in J.T])\n", " S_J = tropygal.entropy(J/sigma_J, mu=_2pi3*np.prod(sigma_J), k=k, correct_bias=True)\n", " return S_J\n", " else:\n", " return np.inf" ] }, { "cell_type": "markdown", "id": "0d9b04f4-74d7-4357-9d43-4bcbf6f81f73", "metadata": {}, "source": [ "## Generate IC" ] }, { "cell_type": "code", "execution_count": 3, "id": "2d640936-bb89-4100-bf6f-9514f6fe831d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "agama.G: 1.0\n" ] } ], "source": [ "_2pi3 = (2.*np.pi)**3\n", "print ('agama.G:', agama.G)\n", "\n", "# True values of parameters:\n", "M0 = 1. # total mass\n", "b0 = 1. # scale length\n", "\n", "n_orbs = int(10_000*10./7)\n", "n_boots = 100\n", "#n_boots = 10_000 \n", "\n", "# Generate IC \n", "pot_0 = agama.Potential(type='Isochrone', mass=M0, scaleRadius=b0)\n", "isoc_data, _ = pot_0.sample(n_orbs, potential=pot_0)\n", "\n", "x = isoc_data[:,0]; y = isoc_data[:,1]; z = isoc_data[:,2]\n", "vx = isoc_data[:,3]; vy = isoc_data[:,4]; vz = isoc_data[:,5]\n", "\n", "v2 = vx**2 + vy**2 + vz**2\n", "pos = np.column_stack((x, y, z))\n", "E = pot_0.potential(pos) + 0.5*v2" ] }, { "cell_type": "markdown", "id": "f03adafa-fd2e-4275-b142-bd852ad30c57", "metadata": {}, "source": [ "## Select 70% more bound stars" ] }, { "cell_type": "code", "execution_count": 4, "id": "14bc2768-248a-49f5-96a3-7b3fbda4c4e7", "metadata": {}, "outputs": [], "source": [ "# This is just to probe a larger range of parameters, \n", "# since we only consider models where all particles are bound\n", "E_cut = np.percentile(E, 70)\n", "cut = E: 0.9872230060864446 ± 0.05886773430277403\n", ": 1.0025291597699255 ± 0.07301148394586504\n" ] } ], "source": [ "#print ('p_bf:', p_bf)\n", "#print ('S_bf:', S_bf)\n", "\n", "# Best fit averaging over all data realizations (bootstraps):\n", "print (':', np.percentile(p_bf[:,0], 50), '±', 0.5*(np.percentile(p_bf[:, 0], 84) - np.percentile(p_bf[:, 0], 16)))\n", "print (':', np.percentile(p_bf[:,1], 50), '±', 0.5*(np.percentile(p_bf[:, 1], 84) - np.percentile(p_bf[:, 1], 16)))" ] }, { "cell_type": "markdown", "id": "5f554be2-1098-4ffe-a81c-8132a5636cb0", "metadata": {}, "source": [ "## Save best fit results" ] }, { "cell_type": "code", "execution_count": 8, "id": "3775d0fe-34b6-442e-8a6e-155285ef1778", "metadata": {}, "outputs": [], "source": [ "# Save best fit parameters for each bootstrap data realization: \n", "\n", "file_name = './fit_results_isoc_boots_k10.dat'\n", "file_out = open(file_name, 'w')\n", "file_out.write('# M_bf b_bf S_bf\\n')\n", "np.savetxt(file_out, np.column_stack([p_bf[:,0], p_bf[:,1], S_bf]), fmt=' %17.8e'*3)\n", "file_out.close()" ] }, { "cell_type": "markdown", "id": "cefe1758-364f-4005-a8f4-c9d01fea84b6", "metadata": {}, "source": [ "## Plot" ] }, { "cell_type": "code", "execution_count": 9, "id": "bafb7a52-b5b8-435a-9e35-fecf9620f5e7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Original data:\n", "M_0: 0.9957639204643787\n", "b_0: 1.002709517712356\n", "-------\n", "Bootstraps:\n", "M_bf: 0.987223006 ± 0.05886773549999996\n", "b_bf: 1.00252916 ± 0.07301148599999996\n", "101 fits\n" ] } ], "source": [ "fname = './fit_results_isoc_boots_k10.dat'\n", "data_file = np.genfromtxt(fname, unpack=True, skip_header=1)\n", "M_bfs = data_file[0]\n", "b_bfs = data_file[1]\n", "#S_bfs = data_file[2]\n", "\n", "# M_bfs = p_bf[:,0]\n", "# b_bfs = p_bf[:,1]\n", "# Original data:\n", "M_0 = p_bf[0,0]\n", "b_0 = p_bf[0,1]\n", "print ('Original data:')\n", "print ('M_0:', M_0)\n", "print ('b_0:', b_0)\n", "\n", "# bootstraps:\n", "mean_M = np.percentile(M_bfs, 50)\n", "sigma_M = 0.5*(np.percentile(M_bfs, 84) - np.percentile(M_bfs, 16))\n", "\n", "mean_b = np.percentile(b_bfs, 50)\n", "sigma_b = 0.5*(np.percentile(b_bfs, 84) - np.percentile(b_bfs, 16))\n", "print('-------')\n", "print('Bootstraps:')\n", "print ('M_bf:', mean_M, '±', sigma_M)\n", "print ('b_bf:', mean_b, '±', sigma_b)\n", "\n", "print (len(M_bfs), 'fits')" ] }, { "cell_type": "code", "execution_count": 10, "id": "c27dd9f4-d520-469c-a250-bbe702203cb9", "metadata": {}, "outputs": [], "source": [ "M_min = 0.5\n", "M_max =1.5\n", "b_min = 0.5\n", "b_max = 1.5" ] }, { "cell_type": "code", "execution_count": 11, "id": "60dd832d-61a7-4b73-9483-0a36263c4caf", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAM7CAYAAABQkS0MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAABJ0AAASdAHeZh94AACtXUlEQVR4nOz9e5Bb+X3feX+a5Fza0oinezzUJbbsOU1JNlvbHQEcSXTFUskELNlJZHkFkGxnk8f7JATWeup5Ns6ugWb+WMlVm+1GO6lsamsdAUw2KjtlXoBaDe1dRzLASiQ7D6UiAbvTQz62RYC6eCTNyNM4c9FwRjPkef5oHcxpNC4HwMHBpd+vqq7hAAe/8zuX3znne363Gdu2bQEAAAAAEJADo84AAAAAAGB/IRAFAAAAAASKQBQAAAAAECgCUQAAAABAoAhEAQAAAACBIhAFAAAAAASKQBQAAAAAECgCUQAAAABAoAhEAQAAAACBIhAFAAAAAASKQBQAAAAAECgCUQAAAGCIwuGwarXaqLMBjBUCUQAAAGBINjY2VKlUZFnWqLMCjBUCUQAAAGAIarWa0un0qLMBjCUCUQAAAGAIksnkqLMAjC0CUQAAAMBnGxsbSiaTMgxj1FkBxhKBKAAAAOCjWq2m69evKxaLjTorwNg6NOoMAG65XE7ValWWZalWq2l7e1u1Wk2JREKZTGbU2QOAkWt3nTx16pSy2eyoswdAUjqd1vnz50edja4KhYJqtZqq1aq2t7d3XVckqV6vjziHmGYEoj0qFAqKx+Mtm1nMz883Cm7z59VqNYDcTbZarUZfCgDowLKsfXedtCxLjz/+eNcRRyORiIrFYl/rSCaTyuVyjXv7/Py8JO26p7vXXywWFYlE+loXpl8ul1M0Gh37JrmlUknxeLzt96ZpBpgb7EcEoj0KhUKNmrlqtapSqdSYF8p9k4rFYjJNU48++ujYX4jGhWmajYC9VCrtu4ctAOOnUCjo+vXrikajYxF4GIaher2u7e1tlUolpdPpqZ8SwjAMnTp1qhEU1mo11Wq1PdtdKpVUqVQUCoV6Xkc4HFYsFlOtVms7zYZhGDp+/LgMw+ABHW1ZlqVisah8Pj/qrHQViURULpcbrSqy2awqlcqu74GhsjEwwzBsSY2/bDY76ixNhVQq1dinqVRq1NkBsM9kMpld1/Z8Pj/qLO3hzmMikRh1dgJTLpd3HRs/90G1WrVN02ykaZqmXSwWfcg1epXNZm3TNH1Pt1wu27FYzDZN0zYMwzYMww6FQnYikbCr1epAacdisT1pOM+J5XJ5oLSHLZ/Pj/01b7+bxDLRCYMV+eD48eO7/p83SP544oknRp0FAPtYc3/LS5cujSgn7fVT+zcNSqWSJO0ZCCaXyw2ctmmajXkfnZY63NeDl0wmlUwmG63O/FCr1RQOhxUOh2WaporFour1uur1us6fP6/t7W0tLCwoGo321dKgUCjoiSeemNga8+vXr+/6f8778TKJZaIbAlEfNDe9dfqWAMAwWZY19c0y++XHvml+mJzUh8tp5PTRbDWInR/BqHMfZ/CnYNVqNW1sbGhubs6X4+hWKpUUDodlWZaq1aoymcyuMh0KhZTP51UsFhvL9vLAb1mWstmsUqmUr/kOkvOCR9p5tqVr2ehNcpnwgkDUBwSew8EFEOgsnU5rbW1t1NkYS37sm3w+36hxjMViYzly9369/9y4cUPRaFSmae6pFfbjOF2/fl2GYVAjNGTpdFrJZFLRaFRzc3NaWFjQ2tqaTp065et6arVao0anWCx2fKkUiUSUzWYbNUVepdPpiX9xQf/Q0ZumMuEFgSgATKhWo3Rjhx/7xjAMlctl2bY9EQOP7BfOYELOg3LzwHbOgEODKJVKPIgPmWVZ2tjY0OXLl7W9va1Tp041mgU6TaP94owM21zj004ikVAsFpNlWR1HlXWUSiUtLCxMdKsJd22oJEWj0RHlZP+apjLhFYEoAEwov5vITBP2zfS6ceOGpDf6xyYSiT3LDFozValUeBAfMsMwZNu26vW6yuWystnsUIL/jY2NxouJXprNnjt3TtJOv8/mIK1ZJpOZ6Ca5kva8vOFFTPCmqUx4RSAKABOKYKs99s30ajWHZ3MwOkhfKucBiwfx6eC8lGge2KqbUCjU6CLUqbl3Mpkcy2b7vXLPwcsURdNt2GWiFwSiADCBWs2jiB3sm+lWKpX21Fa2mne632C0UqnwID4l3HO99zMSv/MyolQqtbymOOfKNIxe7a7h4iXM9Bp2megVgSgATKBJHxRjmNg308t5ydD8oBwKhfYEA/2eB61qXDGZ3OdAPy8W3L9p9WIjnU5PRW1oc7NcmqVPr2GXiV4RiALABPJ7GPdpwr6ZXk6tTasaqOZa0Uql0tegRa1qXMcFTc57467l6+ehe2FhofFvd9NVaef8KpVKmpmZ6frn1ByFw+Fdn/vVz25Qzfk4fvz4iHKCYRtmmegHgSgATJhCoUDT0zbYN9OtU22lH4MWOYHruNSIlkolxeNxzc3NaWZmRgsLC5qZmdHc3Jzi8fhQA5l4PK6NjY2hpT9szU30B639cQbJcoRCIdm23fWvXq83flMsFnd9Ny7nWXNAMQ1NjbHXsMtEPwhEAWCC1Go1nT17dtTZGEvsm+nXrbZy0EGLSqXSWPQPtSxL0WhU0Wi05csVy7JUKBQacw0OoxVAqVTS9evXfU83KM214f3MTe4+DyzL6qtG2r3ecZ331x1QjEtwDP+NS5lwIxAdA5VKRel0WuFwuPG3sLCghYUFJZPJXSdOOp1WoVDoKf1CoaB4PK6FhQXNzc1pbm5O4XBYyWRyKG9T3etzmp+02pZBbWxsKBwO73pTHI1Ge7ohO3M2pdNpxeNxRaPRxn5q3je5XK6xvrm5OUWj0Z63x89jkcvldk187OSt1VxTfuyrTmq1WuMcdrbLvY5h11B1Oo7Nb/Sb8+rsj2HULvhd9kqlksLh8FD3px/HchTnpp/7plAotMx/qwFxvKbnBAzu82BjY4Pa2x606x/qNuigRePQP7RWq+nxxx/3fI2wLEvJZFILCwu+3mMn/dx0B9H9vlho/p2f+3dcNNeSeWmWXqvVlEwmG02Nw+EwXSImwFiWCRsDSyQStqTGX71e9/S7er1uRyIR2zRNO5/Pt/xdsVi0Q6GQnUql7Gw2a0uys9msp/Tz+bxtGIZtGIadzWbtarXa+K5ardqpVMo2DMM2TdMul8ue0vSyPkl2IpGwy+WyXa/X7Xq9bpfLZTuTydimadqxWMzTPioWi419mkqldn0eCoXsbDa7K51yuWyHQiFbkm0Yhqdtcq+j+a9YLDaWi8VidiwWa+xD51hI8rQev49FtVptm+9h7atW6vV64/yPRCK79lm9XrczmUwjX17P2350Oo6ZTKaxXCwWs03TtDOZTOP8rFardjabtU3TtCX5Uh78PN7O8s7x6vXP637361gGeW4Oa98450LzXyKR8LQv3fstFot13Z/uc7RX5XK5r/w1n6OTwLnudtN8PoRCIc/rGPa1qpt6vd64jzp/sVjMLhaLjetItVq1i8Vi4zrSqZz1y9nXfqTVj+brSD/cz2amafadl36up+3S6OXeElQZdT/PeMljKpVqnBvOeenkdVTny34wbWWikdZAv4Zt2/0FotVq1TYMw47FYp7WEYvFejrozvLd0q/X642b9iAPQ876vAQ1mUzG03KtAtFsNmtHIpG2+7herzceIg3D2BUAtFOtVhuBcqtANBKJ7NnnzQ8JnQzrWLTL9zD3VfP6nYegTvl1L9frg3yv+SmXy42bpPshv1qtNgLQTtzlrN/y4PfxrlardiaT2fUXiUQa+YxEInu+d/95fbHh57EM6twc1r5xXlC4X671ev46AUWn/VkulxvpRyIRz2k3p9Fr/ur1+p7gpdfyPwpOUN9N84O11+1z9uUo94X7eSKRSHh6nshms3sC0kFeqLmDYT9eUvfDj4du9/W8l5cRzZrvJ71ylzf3Cymvvxl2GXXvp0772rlnNV+7y+Xyruuu18oY9GaaysSutAb6NWzb7i8QDYVCtmEYPRVY58GsWyDqXBB6eWhyHoj7eZvlrK+XYMZ5wOu0/c2BqFNT0o37wcxroO9wF9JisWhns9mW+9FrIBrUsXAHX0HsK3dA4iWf7mM56EXLC/dxTKVSjVYHXrhvqL3mNajj7Q7wBn0DPexjGfS56ee+aU6vl+PqBMHduPdnPy9q+glE3b/p9QF5lLoF9m7N2+dl3zgvSUfF/aDZz3Wy+SVcP2XA/RJokBqTQfnx0O2uGffrobufVhHu4+JuWdVJkGXU/RKj3X5yzovm86lVwDyqlxfTblrKxJ60Bvo1bNvuPRB1TqZeD56XprlOXnq9gbgvJr1Us7svsL1Wz8disY43yeYHtF7e8LovrL1wP3Q6wXIr+XzeNk3TDoVCbW8qQR6LfD6/68EjFAoNdV85F7Nets0d4A27xqFdTZwXzU3jvN78gzzefgZbwz6WQZ+bfgei/QSKzvXBK/f+7PVhs9+mue7mx72+FB0F577p9dxpvi97OYcikUjPLy/95Jy7g+ShXC7vaVrutWmnc96674GjMq4P3f20YnOa1zp/zufdrhFBlNHm/dzqmukEoa1ejrRqfTDu15JJNellom1aA/0atm33Hog6BbfXN57OA2u7G4r7gamfG4g7qOz1jV0vD13um12nt8/u7TFNs6ema+4Hu17ezrkvqk6/0H4EfSzcD/vD3lfuvPWybYPW/PTCfRz7qeHo9fdBH2+/gq0gjmWQ56Ztj0cg6vR99cp9Le21iW6/gahTU5NKpSbiwdFr/1BHqxqlbud4Ly+DhsE53/14UdeqdtS5r2WzWbtYLNrlctnO5/ONViN+PlwOyo+Hbvc2jctDdy+CKKPNgWTzi7B2NaGO5vEZhtVH1N2NYdR/oyob01omCER90Gsg6twg+nkYN02z6wAb/TYtctfMeDmx+m3C2HzDa6f5AtdLQNncxNarbhdlr4I+FkHuq0EuhP3WVPfKfRz7DXrdee32cBr08fYr2AriWAZdjv0ORPsJ9JzlDcPwnO9+y0a/geik8do/1K35XtPpwWsc+oeapulrc1j3wF+9/PXbX9lP0/rQPW669Q8NhUKeXjYmEomh1qCPUyA6qvIxrWWC6VtG4NFHH5Wkvqa1aDfJcKFQaMzl0+/Q8+6509zptVKr1XYNLd9qInGv6/SqlwmW3XN19TsEvWEYfe3LoI9FK8PaV+6pUPrZtuPHjzf+PcyJ2P1w6tSpxr8zmUzb5cbhePdjVMcy6HIcNPexsyyr5ZQ1rbiHxA/i+E+abvOHttK87yuVStt9Ow7zh25vb/u6/lAopHK5rHw+7zndTCajYrHoWx4w3tzX7uZrczQalWmaymazHdOIxWLKZrOKxWJDyaOTt3q9LnunAm2kf5QPfxGIjoC7sD/++OM9zcHjXBiauS8Uvd6s3dwPpJ0uPu65TA3D6CmgzGazjd90esB36+Xhtdn29nZfv+v3gSDoY9FsmPvq0qVLjX/3s23uvAU1H1s/EzZLUjweb/y7Vqu1ze+oj3e/RnEsR1GOg9bvdcMddBOI7uZl/tBW3C+THO3uOeMwf+ipU6f6nqu2k1gspmq1qnK5rEQioVAo1LguOi9cs9ms6vW6UqmU7+ufJu5yOuksy9r1gs99/kejURmGoXw+P4KcYZIMWiYIREcgEok0HlYsy1I4HFY8HvdUq5BIJFreLN2/HeSkWFhYaJlmM/dDrLtmxItIJKJ6va56ve65JnUUF/9+HyiDPhbNhrmv3AFHP/vHaQ0gSdVq1Zc8DUvzed3uGIz6ePdrFMdymh7iOnFeJIRCIZ0/f97Tb/p9YbIfOOWh1xcZhmHsqaW5fPly23UM8iLJD0HUKmWzWZXL5UbtUr1eV7FYVCKR4Bz0YJr2UfN9JhqNqlKpaGFhQaFQiCAUngxaJghER6S5gBcKBUWjUc3MzCgajWpjY8NzLUPzcoOcFO7fdlr/oA+xvRrFxb+fh+ZRHItOv/NT802rn/W4fzPuNVzNNf3Xr1/fs8w4HO9+jOpYTtNDXCeJREK2batcLg9UC4wdg9RWnjt3btf/W5a1q0WP9EZ5G3WNqBfUlns3jOuN+wXcpGtuYprJZBQOh5XJZDy3VsNkGccyccinfKBHTt+NeDy+58ZSKpV2PShGIhGl0+m2N8nm3w9S69D821qt1jXQnNaHS3cNlVejPhbD1NxHr58mZOMefDabn59vbHerB8BJPd778ViOi1KppEqloueee67R5LRWqxFgdFAqlfYElF6FQiGZprlr/66tre2qeRyH/qHtlEolZbNZlUqlXeXWaVKbTCaHFkDH43E98cQTE9tcd7+0wOhX8wtJpxxcunRJtVpNkUiEF2lTZhzLBIHoCIVCIVWrVRUKBa2trbWtBXEC01gs1rKpRPMDoZ+BYasBQpo/m6Y3hIMK+lgEqflBuVgsjuWDm5/cx6/V/p/U470fj+UoNV/jY7GYotGoTp8+LdM0ZRiG4vH4npo69N8/1C2ZTO4auMgZtMg558ehf2gzy7I6dtlxanYLhUJjvIV+Bw1sZ9wHlOum2/V70DQnmfMCzOFck8rlskqlUuNaZJqm0um07+cWRmMcywRNc8dALBZTuVxWtVpVNptVIpFo+VBYKBQUDof3fD7MB1YvtR7PPffc0NY/aUZ9LII06sB4HEzL8eZYDkcul9Pc3Fyj5Us2m5Vt28rn83sGjUFr/fYPdWv1EO0eEGwc+oe61Wo1Pf74454DQcuylEwmtbCw4Guz/km/Lrifo/q9njbvg3GsUepH87l1+vRpJRIJZbNZVatV5fN5GYahWq3WOLdotTH5xrFMEIiOEdM0d10I6vW68vn8rje1lUplTzO65qDVz5tHq4C4+cFp0m9Wfgr6WARp1OsfBffxa7X9k3q89+OxDJJTm5VMJhu1eXfu3KFWoQ9+1Fa2GrQol8tJGr/+oc4Ahu5rSSwWU7FYVLValW3bqlarKhaLSqVSu+7HtVpN4XDY85RBnTj7Z5KvFU888UTj3/1em5sf1qelqWpz/9Dm8z8Wi+nOnTuN86tWq43Vyxr0ZxzLBIHoCBQKBU9vLZ2bZ7FY3NUkN5fL7Xoz1RwYDlKT4rXP2yQNOhOkURyLoPi5bZPCvY2taq4m4Xin0+mugyrth2PZSqt944eTJ082mrYlEgkVi0VqPvvkV21l8wtcp2nruPUPTafTjQfERCKx64W0k0fTNBWJRJTJZFSv1xtTojk2NjYGqh11z397+vTpgbZnlPx4Udh8XR6X82RQN27caPzb6R7QrHn6llqt1nhBEbRKpaK5uTnNzMyM/G+SA/JxLBMEoiNQLBZ7nifQmTDY4e5L1DzNxCC1Mu5mtu0uTpL21NJixyiORVCa35hOe0148xxrrW4+k3C8K5XKnkBzvx3Ldlrtm0Elk8nGNdHLZPBoz4/+oQ53IOfIZrNj1T/U/aCfyWT2BJjtOAGre1ChfmtH3TWypmlOdA1gc977aVrq/s0k74tm7ue2Tud/83fNNalu+6Xp7qifxQYxjmWCQHRE2s1l1onTp0jaPZWEYRi7ToZW00x45fXi5H4o3y8XHy9GcSyC5M5HpxvSNGjuQ9PqGEzy8d5PxzIolmXtqjHoZQqETi8DnKa++40f/UPdmmtFnYEAx6WGw3nBHIvF+hqpNpPJqFwu7wq4NzY2NDc356kmyxmHwrmnT8MUHu4m2f08q7jnSJ7k2mG3VvOHduIl8KrVan3NNOBVKBRqzH076r9Jn1913MoEgeiINI9Y5pVzI21+aHGfDIPUULqba3R68Dl16tSu/+93dL1SqaSNjY2+fjuugj4WQXLnw52/XlUqlcCa+PRb23fp0qXGv03TbNv8ZFKP9yQey3HXvB+a+yV20qlmdnt7e1+2PPG7trJdH93mlg2j4rwQGiQAdEbjdweyzmBGMzMzisfjyuVyjWmECoWC0um0FhYWdk0nF4vFejp/x5X7+tzPCzf3s8007A+pe//QZt3GSpB2XmJMU43xNBu3MkEgOkL9BG/ORaD5DZV70ILm+ca8qlQqjd+FQqGOFyfDMHbd1Pttfua16dEkCfpYBCmRSDS2zZkCoR9ra2tj3xzU3fy904PhuB/vdvOR7qdj2Y7fc7W6a8R7fSijZclepVLJ14fbVoMWNbdqGCXnfPTjnHRqR5u3rVAoKJlMKhqNKhwOKx6Pa2NjY89UHpNe6+Noni+2V84LIGc+2mng3g/duoE0X5fcg924FYvFqakxnnbjViYIRH3Qbx+jft56OheFVhcD94Tf/dRQrK2t9ZQ3d9+TQqHQ88OoM1hEc+3qNAj6WATJnZ9+z+FCoRDYCKL9XGjd53YoFOr61i/o4+2++HcLYLa3t9sOfDRpx9ILv/ZNP9zXwF7SdfpCtkpnWCzL0sbGxq7BccaJ81Kn3YNvv5qb545Lbai0cz76GeyEQiGVy2Xl83lP6TrzkU5LEOpwaodrtVpPL3y8vowcFmfQqI2NDV/LaC/dQJrX2+qljWVZKpVKY3UfQGdjVSZsDCwWi9mSGn/VarXj8olEorFsPp/va13t1hGJRGxJtmEYdr1e95xuuVxu5CmVSnn+XTabbfwuFot5/p1t7+yHTuvK5/ONtCORSM9pO7/NZrOef5fJZPraD60EeSyC3lfuc75cLve0vlgsNvC+7cZ9XvZ6DtTr9cbvDMPoWp4dQR5vdx5N02y7XLVatQ3D6JjWsI9l0Oemn/vGtm27WCx6vsa58+slbUcsFuu4reVy2Q6FQgPnz800zV157eWcDYKTv17Kbq9pS7IzmYzv6fcrkUj0/EzQi3K5bCcSCTsUCtmGYTSOfSQSsbPZ7FidA+7r4qCPqvV6vbG9vZSRUCjU13XLL+7ztF3575X7euHlGdR9PZXU8hzJZDIj20f7ybSWCQJRH7gvFl5unO4Hjk4PSs2q1WrXk6Zerzfy4/VEcf+m12DStvt7WMzn8123PZVK9bWfbPuNwiLJTiQSnn/nfigftKAFeSxGsa+c3/XyEJtKpXy7oXbS/ILEMAzPQZZ7fxSLRc/rDLrsOYFvpxdTqVTK04P2MI/lKM5NP/eN++VUt+1tflDwcv7k83k7Fos1ru+ttjWbzbbdfnf+vO7f5nz2eq4PU7FY3HXcTdP0/DLIK/c+6/XlC4Lhvm54CZi6cQdgXs515x4yqpc0rcqoH+dq8371sm1OwNIu+DEMY2yuH9NsWssEgWiPqtWqnclk7EwmYycSiT1BqPvh11muOThzArdEImEnEglPD631et0OhUK2aZqeTgDnQSwSiXRcvlqtNrahlwe9Zu4C0u3hzglCmx8u6vV6Y5811zI7D4HOw2PzBblcLtuZTKbxcNz820gk0kjbvT+cdaZSqV0Pr63W2fxbr4ZxLEaxrzptm5dAL5VKeT5/B+UORFOplF0ul23TNDvm0SljXrennaDKnjtwaRXEOjV+Xve3X8dyHM7NQfeN1/y32zfOst32fz6ft0OhUGMZ9zXI/bvmc7ff/etoruXoFLAPW71e31NL1+ovFArZsVjMlxpD9/ZjtPL5fON8dp6rWpV9p/y778eZTKanAMhpndHtGuc8oA9yHxhUqzLqR+sA9zXG64sr9/20eX9EIhFqQ32238oEV+EeuQ+al79WhT2RSOx6QMlkMrZpmi2bxjiBr2EYdiwW6+khPpvNNvLgPIzb9s4Fzmmi4+TPj7dZxWKx8WAdCoXsfD7fyG+1Wm08dLV7QHceHr3s0+bmgO43Nd1+6y5E7sLl5bf9Pqz5fSxGsa+8bFssFrOLxWLj+JbL5cb3iUQisDfLzYGobb8R+CUSiV3H0QmcnG3wI59BlT33G81IJGJXq1W7Xq83Xvb0mrYfx3Jczs1+943zAOglD172o/MA6d5f5XLZjsVie66FToDsnBtOQO7eT73kr9NLwXw+31h2GM1fvWpuKthpe1qdM/1KJBJ9tQCCv1o9YHc7r91/vR5D56WkUz6ay59zbe72IjEIznXEMAw7FAr50ozcXXnSS1lyXnpFIhG7XC43rmFBvVzeT/ZbmSAQHYFUKrXnrW61Wm08dLhPMOchZJC31U7TL+dEMwzDNk2z8aDpt2KxuGd9oVDITiQS+74ZVNDHIkitts2P87cfrQJRdz4jkYhtmuaem7zf+QzieNfr9V3XDvfDgp/5HtWxHMQw9k0v3NdC51xzjn+n2kqnhUan5QD0r13ZbH5ROS6y2awvL2Gcl679BBXZbNaORCKNfTXssR4QrFGViRnbtm0BwBTJ5XKN0TFTqdTYjTwMAIBXznzr7vlhgWnA9C0AAADAmHruuedkTNmc64BEIAoAAACMrVqtNlZz3gJ+IRAFAAAAxlSpVFIoFBp1NgDfEYgCAAAAY6hWq8k0zVFnAxgKAlEAAABgDGWzWZ0+fXrU2QCG4tCoMwAAfrMsq+W/AQCYJLlcTnfu3Bl1NoChoEYUwNSpVquNf9dqtRHmBACA/qTTaSUSCUbMxdRiHtExYVmWvvSlL+nHf/zH9dBDD406O8BEeeGFF1QoFGRZlm7duqWvfOUru77/4Ac/qGPHjskwDEWjUf34j//4iHIKAEB33/rWt3T27Fl94QtfGHVWMCVeffVVfetb39KHP/zhsXm5QSA6Jq5cuaJPfOITo84GAAAAgCn15JNP6pd+6ZdGnQ1J9BEdG04NzZNPPqmjR4+OODeT5/bt2/rEJz7B/kPgOPcwKud+5z/rW9/6pj5z+gnOPQSK6x5GhXOvf86+G6dWYQSiY8Jpjnv06FEtLi6OODeTi/2HUeHcQ9CMtz+t515+nXMPI8O5h1Hh3OvfOHUBZLAiAAAAAECgCEQBAAAAAIEiEAUAAAAABIpAFAAAAAAQKAJRAAAAAECgCEQBAAAAAIEiEAUAAAAABIpAFAAAAAAQKAJRAAAAAECgCEQBAAAAAIEiEMVUeOyxx/TpT39ajz322Kizgn2Gcw+jcujQQT322BHOPQSO6x5GhXNvuszYtm2POhOQbt68qfe+97166qmntLi4OOrsAADG3Nl//SVJ0vlf+/CIcwIAGHfjGGtQIwoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUPsqEC0UCorH44Gsa2FhIZD1AAAAAMCkmfpA1LIslUolxeNxxeNx1Wq1oa8znU6rVqvJsqyhrwsAAAAAJs3UBqKWZWlhYUHhcFjFYlHJZDKQ9VYqFW1sbASyLgAAAACYRIdGnYFhMQxD1Wq18f9B1U6ura0pFAqpUqkEsj4AAAAAmDRTWyM6ChsbGzp37pzm5+dHnRUAAAAAGFsEoj6p1Wp67rnnFAqFRp0VAAAAABhrU9s0N2iZTEbZbHbU2QAAAACAsUeNqA9yuVxg08IAAAAAwKQjEB2QZVmqVquKRCKjzgoAAAAATASa5g4onU4rk8n4lt7t27d7Wv6xxx7TkSNHfFs/AAAAgPH27LPP6nvf+57n5XuNMYJAIDqAQqGgaDQqwzB8S/MTn/hET8t/+tOf1mc+8xnf1g8AAABgvP32b/+2fvM3f3PU2RgIgWifLMtSsVj0fYCiJ598UkePHvW8/GOPPebr+gEAAACMt0996lM9jVFz+/btniu8ho1AtE9+N8l1HD16VIuLi76nCwAAAGA6HDlyZOK75zFYUR9KpZLC4bCvTXIBAAAAYL8gEO1DPp9XIpEYdTYAAAAAYCLRNLdHhUJBly9fVqlUartMrVaTJIXD4cZnmUxGsVhs6PkDAAAAgHFHINqjWCzWNaBcWFhQrVZTuVym+S4AAAAANKFpLgAAAAAgUASiHViWNeosAAAAAMDUIRBtIx6Pa25urqf5eQAAAAAA3e2bQNQZXKhWq3mq6SwUCrv+65VlWY3Bim7cuNFbJgEAAABgH5jqQDQcDmthYWFXzaZlWZqbm9PCwoIWFhaUTqdb/jaTycgwDGUyGU/rikajWlhY0OOPP77ns4WFBZr5AgAAAMAPTfWoueVyue/fplIppVIpz8sXi8W+1wUAAAAA+8lU14gCAAAAAMYPgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAAAAAAgUgSgAAAAAIFAEogAAAACAQBGIAgAAAAACRSAKAACAnmxubmplZUVLS0taWVnR5ubmqLMEYMIcGnUGAAAAMDk2Nzd14sQJ3b17V5K0tbWlK1eu6Nq1a1peXh5x7gBMCmpEAQAA4Nn6+nojCHXcvXtX6+vrI8oRgElEIAoAAADPbt682fLzW7duBZwTAJOMQBQAAACeLS4utvz82LFjAecEwCQjEAUAAIBnq6urmp2d3fXZ7OysVldXR5QjAJOIQBQAAACeLS8v69q1azpz5oyWlpZ05swZBioC0DNGzQUAAEBPlpeXdeHChVFnA8AEo0YUAAAAABAoAlEAAAAAQKAIRAEAAAAAgSIQBQAAAAAEikAUAAAAY2Nzc1MrKytaWlrSysqKNjc3R50lAEPAqLkAAAAYC5ubmzpx4oTu3r0rSdra2tKVK1eYHgaYQtSIAgAA7BPjXtu4vr7eCEIdd+/e1fr6+ohyBGBYqBEFAADYByahtvHmzZstP79161bAOQEwbNSIAgAA7AOTUNu4uLjY8vNjx44FnBMAw0YgCgAAsA9MQm3j6uqqZmdnd302Ozur1dXVEeUIwLAQiAIAAOwDk1DbuLy8rGvXrunMmTNaWlrSmTNnxqrpMAD/0EcUAABgH1hdXdWVK1d2Nc8dx9rG5eVlXbhwYdTZADBk1IgCAADsA9Q2Ahgn1IgCAADsE9Q2AhgX1IgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAK1rwLRQqGgeDzua5qlUknRaFQLCwuam5tTOBxWMplUrVbzdT0AAAAAMC2mPhC1LEulUknxeFzxeNzXAHFjY0PFYlH5fF7ValX1el2ZTEalUkkLCwtKp9O+rQsAAECSNjc3tbKyoqWlJa2srGhzc3PUWQKAnh0adQaGxbIshcNhSVIsFlMymVShUPAt/VwuJ0nKZDK7Po9EIqpWqwqHw9rY2NDCwoISiYRv6wUAAPvX5uamTpw4obt370qStra2dOXKFV27dk3Ly8sjzh0AeDe1gahhGKpWq43/tyzL1/Tz+byKxWLb78+fP99ophuJRGSapq/rBwAA+8/6+nojCHXcvXtX6+vrunDhwohyBQC9m/qmucPgpa9pKBRSKBSSJGWz2SCyBQAAptzNmzdbfn7r1q2AcwIAgyEQ7cP169eVTCa79gF1akEZuAgAAPhhcXGx5efHjh0LOCcAMBgC0T44zXydfqLtzM/P71oeAACMv3EeDGh1dVWzs7O7PpudndXq6uqIcgQA/ZnaPqLD5NSERqPRjss5NaFOE10AADDexn0woOXlZV27dk3r6+u6deuWjh07ptXV1bHIGwD0gkC0D6Zpeur3eePGDUndA1YAADAeJmEwoOXl5bHJCwD0i0B0SAqFgizLkmmaikQinn93+/btntbz2GOP6ciRI71mDwAAtMBgQAAmwbPPPqvvfe97npfvNcYIAoHokKytrUnqfcTcT3ziEz0t/+lPf1qf+cxnevoNAABobXFxUVtbW3s+ZzAgAOPkt3/7t/Wbv/mbo87GQAhEhyCXy6lSqSibzfZUGypJTz75pI4ePep5+ccee6zX7AEAgDZWV1d15cqVXc1zGQwIwLj51Kc+1XU6Sbfbt2/3XOE1bASiPqvVakomk8pkMkokEj3//ujRo22HZgcAAMPFYEAAJsGRI0cmvnsegajPotGoUqmUUqnUqLMCAAD6wGBAADB8zCPqo2g0qkgkokwmM+qsAAAAAMDYIhD1STKZVCgU6nlwIgAAAADYbwhEfZBOp2UYRtuaUMuygs0QAAAAAIwxAtEB5XI5WZbVsTluOp0OMEcAAAAAMN4IRDvoVpNZKpVULpc7Nset1WpaWFjwOWcAAAAAMLkIRNuIx+Oam5trOz9PrVZTPp/v2ie0UCjINM1hZBEAAAAAJtK+mb6lVCpJ2gkgLcuSYRgdly8UCrv+62ZZlsLhsCTp8uXLbdNwalSr1WofOQYAAACA6TTVgWg4HJZlWdre3m4EhZZlaW5urlFLGYvFWvbvzGQyWltb07lz5/Z8d/bs2Z4GIKJGFAAAAADeMNWBaLlc7vu3qVRKqVSq5Xf5fL7vdAEAAABgv6OPKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAAAAAIBAEYgCAAAAAAJFIAoAAAAACBSBKAAAAAAgUASiAABgX9rc3NTKyoqWlpa0srKizc3NUWcJAPaNQ6POAAAAQNA2Nzd14sQJ3b17V5K0tbWlK1eu6Nq1a1peXh5x7gBg+lEjCgAA9p319fVGEOq4e/eu1tfXR5SjvaixBTDNqBEFAAD7zs2bN1t+fuvWrYBz0ho1tgCmHTWiAABg31lcXGz5+bFjxwLOSWuTUGMLAIMgEAUAAFOpU9PW1dVVzc7O7lp+dnZWq6urQWezpXGvsQWAQRGIAgCAqeM0bb148aK2trZ08eJFnThxohGMLi8v69q1azpz5oyWlpZ05syZsWr2Ou41tgAwKPqIAgCAqdOpaeuFCxck7QSjzr/Hzerqqq5cubJrG8apxhYABkWNKAAAmDqT3rR13GtsAWBQ1IgCAICps7i4qK2trT2fT1LT1nGusQWAQVEjCgAAps64D0YEAPsdNaIAAGBqbG5uan19XTdv3tSHP/xh2bat73znOzp27JhWV1dp2goAY4JAFAAATAVnpFxngJ+trS3Nzs723bfSHdQuLi4SyAKAj2iaCwAApkKnkXJ71W36FwDAYAhEAQDAVPBzpFw/g1oAwF4EogAAYCosLi62/LyfkXInffoXABh3BKIAAGCoNjc3tbKyoqWlJa2srAyteaufI+X6GdQCAPYiEAUAAEMTZF/L5eVlXbt2TWfOnNHS0pLOnDnT90BFTP8CAMPFqLkAAGBoOvW1vHDhgu/rW15e9iVdJ6hdX1/XrVu3mP4FAHxGIAoAAIZmkvta+hXUAgD2omkuAAAYmkH7WgbVvxQAECxqRAEAwNCsrq7qypUru5rneu1r6fQvdX67tbWlK1eu9N3vEwAwPqgRBQAAQzPIAELM5QkA04saUQAAMFT99rWc5P6lAIDOqBEFAABjibk8AWB6EYgCAICxxFyeADC9CEQBAMBYGqR/KQBgvNFHFAAAjC3m8gSA6USNKAAAAAAgUASiAAAgEJubm1pZWdHS0pJWVla0ubk56iwBAEaEprkAAGDoNjc3deLEica8oFtbW7py5Qp9PgFgn6JGFAAADN36+nojCHXcvXtX6+vrI8oRAGCUCEQBAMDQ3bx5s+XnTz75JM10AWAfIhAFAABDt7i42PLzV155RRcvXtQHPvABglEA2EcIRAEAwNCtrq5qdna27fevvvqqVldXA8wRAGCUCEQBAMDQLS8v69q1azpz5owOHGj9+PGVr3wl4FwBAEaFQBQAAARieXlZFy5c0COPPNLye9u2A84RAGBUCEQBAECgPvjBD/b0OQBg+hCIAgCAhs3NTa2srGhpaWloo9lmMhk99NBDuz576KGHlMlkfF8XAGA8HRp1BgAAwHjY3NzUiRMnGvN9bm1t6cqVK7p27ZqWl5d9W8/y8rK++tWvan19Xbdu3dKxY8e0urrq6zoAAOONQBQAAEiS1tfXG0Go4+7du1pfX9eFCxd8XZfTXxQAsD/RNBcAAEiSbt682fLzW7duBZwTAMC0IxAFAACSpMXFxZafHzt2LOCcAACmHYEoAACQJK2urmp2dnbXZ7Ozs1pdXR1RjgAA04pAFAAASNrpt3nt2jWdOXNGS0tLOnPmjO8DFQEAIDFYEQAAcGEQIQBAEKgRBQAAAAAEikAUAAAAABAoAlEAAOCrzc1NraysaGlpSSsrK9rc3OzpewDA9KOPKAAA8M3m5qZOnDihu3fvSpK2trZ05cqVxqBH3b4HAOwP1IgCAIA92tVadqvNXF9fbwSZjrt372p9fd3T9wCA/YEaUQAA0LC5ual0Oq0/+qM/km3bkt6otfzc5z6nX/3VX+1Ym3nz5s2W6d66dcvT9wCA/YEaUQAAIOmNZrVf/OIXG0Go4+7du/qN3/iNrrWZi4uLLdM+duyYp+8BAPsDgSgAAJDUutms27PPPtvyc3dt5i//8i/r4MGDu76fnZ3V6uqqJGl1dVWzs7NtvwcA7A8EogAAQFL7ZrOOI0eOtPzcqc3c3NzUr/7qr+revXuN7w4ePKjPfe5zjaa7y8vLunbtms6cOaOlpSWdOXOGgYoAYB8iEAUAAJLaN5uVdmotf+u3fqtjbWarGtV79+7p85///K7PlpeXdeHCBW1uburChQsEoQCwDxGIAgAASTvNZh966KE9nx84cECf+9zndOrUqY61mQxEBADwikAUAABI2qmp/MhHPrLn8/v37zdqNTvVZjIQEQDAKwJRAADQmB/0P/2n/9Tyey+1mgxEBADwinlEAQDY55xpWzqNmOulVtMZiGh9fV23bt3SsWPHtLq6Sh9QAMAeBKIAAEyRzc1Nra+v6+bNm1pcXPQUCHabtqWXWk2n6S4AAJ3QNBcAgCnh1GxevHhRW1tbunjxosLhsN71rndpZWVFm5ubLX/XbpChhx9+eKKnV3GaGy8tLXXcfgBA8KgRBQBgSrSbPuX27du6ffu2rly50jKoXFxc1NbW1p70PvGJT0xs7WZzc+Otra222w8ACB41ogAATIlyudzx+7t372p9fX3P59M4yFCroLzd9gMAgkcgCgDAlHj11Ve7LtNq9FtnkKF284N2M45NYJnTFADGG01zAQCYEg8++GDXZdqNftvvIEPj2gS2XXNj5jQFgPFAjSgAAFPi+PHjHb93mtv6WYM5rk1gp7G5MQBMEwJRAACmRKvg6+DBg3r3u9/daG4rac/IuidOnOg7GB3XJrCDNjcGAAwXTXMBABgj/cwD6nCCr/X1dd26dUvHjh3b8/uVlZW2NZj9NM0d5yawzGkKAOOLQBQAgDFx+fJl/cqv/Iru3bsnqb/+lt2CL79rMFdXV3XlypVdwe04NoEdJMAHAPiPprkAAIyBzc3NXUGoo1V/y83NTT39V3+larXacx/PxcXFlp/3W4M5CU1gnQGV/GqODAAYHIEoAABjYH19fU8Q6nDXVjpB1fMvPK9XX32lZVDVaTCiYQzi49TCbm5u6sKFC2MVhErjO6ASAOxnBKIAAIyBdk1mpd21ld2Cqm61f5NQg+m3cR1QCQD2M/qIAgAwBtoN+nPw4MFdtZXdgqpOgarTd3S/DeIzzgMqAcB+RY0oAABjoN3UK7/3e7+3q7ayWx9Pav/2Yk5RABg/BKIAAIyBVk1my+WyTp06tWu5bkGV34MRTYP92BwZAMYdTXMBABgTXpvMfvjDH9ZLBw5Kkj760Y8qk8k0gqpJmU4laPutOTIAjDtqRAEAmBDOQERf+MIXdP/+Pd2/f09f/vKXdy1D7R8AYBJQIwoAwITwMhCRRO0fAGD8USMKAMCEYCAiAMC0IBAFAGBCMBARAGBaEIgCADAhmIYEADAtCEQBAJgQ7oGIHnroYR1+y2EGIgIATCQCUQAAJogzENHCwoL+xo/9GEEoAGAiEYgCAABfbW5uamVlRUtLS1pZWdHm5uaoswQAGDNM3wIAAHzjzHXqTDOztbWlK1eu0IQYALALNaIAAEwId03j03/1V3rhhRfGruZxdXW17VynAAA4qBEFAGACNNc0PvLh5/X8C8/r/3vxoqTxqHnc3NzUF7/4xZbfMdcpAMCNGlEAACbA+vr6nprGZqOueVxfX5dt2y2/Y65TAIAbgSgAABPg5s2bnpYbZc1juzzOzMww1ykAYBcCUQAAJsDi4qKn5UZZ89gujz//8z/PQEUAgF0IRAEAGCGvU52srq5qdna2Y1qzs7MjrXlslcfZ2VllMpkR5QgAMK4YrAgAgBHxNNXJ9rb0ta9p+fZtXfvUp7ReKunWs8/q8OyP6JFH3qwz//V/rVu3b+vYsWNaXV0NvOZxc3NT6+vrunnzphYXF/W5z31On//853Xr1q2R5QkAMP4IRAEAGJFWAxDdvXtX62truvDrvy794R9Kf/zH0osvSt//vpZfeUUXbFt6+GGdvX9PeuFFna/XpU98QvrFX5T+q/8q0PwzZygAoF8EogAAjEi7wX1uXb0qff3r0rPPSg89JL3lLdLb3iY9/LB04IB0/7506NDOf7/zHel3f1f64helD35Q+vVfl37iJwLJf9tAen1dFy5cCCQPAIDJRB9RAABGpOXgPj8mHTFeku7eld71LunoUenIEelNb5IOHpRmZn743wPSwUM737/rXTvL/8EfSMmkdOWK1GYaFT+1DaSZMxQA0AWBKAAAI7JncJ856UBEqv/CIdWOPio9+KC3hB58UHrnO3f+/uIvpI2NnVrSFsGo18GRvGg3Si5zhgIAuiEQBQBgRJaXl3Xt2jWd+dmf1bsefliHl2Z08PCM/uLQXf3LN93pPcE3v1n6qZ+SvvtdKZuVfv/3d33t9Om8ePGitra2dPHiRZ04caLvYLTdKLnMGQoA6GZfBaKFQkHxeNzXNC3LUjqdVjKZbPw3mUzKsixf1wMAmE7LhqHV11/XN2df0fOP23rtdVsv/eCefvvg1/V//+CZ3hM8dEh697u1+fWva+W/+++09NM/3aj57NSns6+8O4H0mTNaWlrSmTNnGKgIAODJ1A9WZFmWbty4oWw2q0KhoFAo5FvatVpN0WhU2WxWkUik8XmlUtHjjz+uq1ev+ro+AMCUuX9f+pf/Uuv/5b/o1ZCkRyQ9/8OvDkvnDv25/rbe2nOymy+/rBPf/a7u3r8vffe72vrzP9eVK1f0N/7G32i5/CB9OpeXlxmYCADQs6mtEbUsSwsLCwqHwyoWi0omk76vIxqNKplM7gpCJSkUCimTyejkyZPUjALAlBq0r+Xm5qZWPvYxLeVy+r8fellakHRf0ms//Lsv1d7+fdUe/H7PeVu/fXsnCHW5e/eufvCDH7Rcnj6dAICgTW2NqGEYqlarjf/3OyDc2NhQrVZTIpFo+X0ikVA6nVY6nVY2m/V13QCA0Wo1f+alS5f08z//88pkMl2bpjb/Xse1qzZUkvSiNPPojC7OfFv/9Lvv6il/N198seXnDz/8sGZnZ3c1z6VPJwBgFAKpEX3hhRf0a7/2a3rXu96lgwcP6uDBg3r00Ud1+vRpff7znw8iC75bW1tTJBKRYRhtl4lEIrp8+XJwmQIABKJVX0vbtvXFL37R0+A/u34/p921oT904DXpRw89qC8e/l7PtaKLjzzS8vNQKESfTgDAWBh6IPpv/s2/0dzcnHK5nKrVqmzblm3bsixLhUJBsVhMjz76qP7Fv/gXw86KbyqViizLkmmaHZczTVOWZalUKgWUMwBAENrNnyl5G/xn1+/fq53a0KZKzEcffFCP3/sRPf3AK7o4/21tvvCCVioVLX3pS1qpVPRqU9Nbx+YLL8h67TXNNH0++9BDWl1dbfTp3Nzc1IULFwhCAQAjMdRA9E//9E+VSCR0+PBhffKTn1QikVAsFpNpmo2A1LZt1et1pVIpvf/979c3vvGNYWbJFzdu3JAkLSwsdFzO+b5SqQw9TwCA4LSbP9Nx69atjn1IG79vUxsqSa/cv6c33T+oQ/aM/s83fVcfeOqPdfHb39bWiy/q4re/rdrLL+8JRjdfeEEn/uRP9IXvfU/ODKIzkj76pjfp2v/0PxF0AgDGxlAD0bW1NaVSKW1vb+vy5cv67Gc/q8uXL+v27du6f/++isWiEomEDMOQbdu6ceOGwuHw2AejTt/TTs1yJWl+fl6SdP369WFnCQAQoFbzZ7q94x3v6Dhf5+rqqmYfeKBtbagkPXJoZxiHd7z2sP7iwEt6ddHe9b0tW3/dNPhQq0GKbO3Eu8sPP9zrZgIAMDRDrxHt1Dzp5MmT+uxnP6vt7W390R/9kU6ePKnt7W2Fw2G98MILw8zaQLwOfOQEqoycCwDTxZk/82Mf+5hmZnY3gp2dnZVt2x3n61xeXlb+//UreuxdB3RIM3tqQyXpbQ89JEl60/2D+sHr93dqTud2L/NKU9DZbpCiWz/4gfRai5UAADAiQx019/Dhw56XjUQiikQiKhQKOnXqlE6dOqUvfOELQ8xd/7a3tyW9UePZTa1W85z2ud/5zzLe/rTn5Q8dOqiDB6d28GMAGGs/9vFV/Tcf/cf66+/9tV559VU9/NBD+tHHflR/9VdP62fe9fE9yz/z0MM6+6+/JEn66/mTWjz8Ad27L917i71n2UMzB/TKgRn94L6t47q/E4SelfTyzvez82/XQc3o7HtPvZGfox/XI6/vDTgPHzigs68dkX64bgDAZLt373W9/vo9z8tb3/n68DLTp74imD/7sz+TaZp6y1ve0nE50zT1wgsvdF3OLRaL6caNGzp+/Lj+43/8j/rIRz7STxaHapg1nN/61jf13Muve17+sceO6LHHHhtafgAAnT300MP6Gz/2Y7s+e/ihh/Tqq6/sWfbhH9Zyvnb/B3ppZqdZrd16zCHdt229fK/pywclvSrp3k7fz4MzM3r1/n399Q9+oFfu39cDMzOa0YxsvRHYzmhGP3rwkHSIl5YAMC22t+v63vee9bz8y899e4i56U9fd6VQKKSZmRmZpqlYLKZoNKrjx4/vCThXV1d19uxZXbp0qef0P/vZz+qzn/3sWAaiw/SZ00/o6NGjnpd/7LHHdOTIkSHmCADQqz3zhGqnya4zVco/+/I/05985bN65zPPq/biK/r2D/bWYj588IBecQeiD0g6LD30Xw7ol7/2Nj3w3/wvkqTCP/9/7uoX+tDMjD7yoz+qb7/yio498ohW5+d3+of+7/+79P73D22bAQDBefbZZ/W9733P8/K3b9/WJy4MMUN96CsQNQxDlmWpWq1qY2NDGxsbknYCyEgkomg0qp/7uZ9TKBTS3Nyc/u2//bf6h//wH/a0jkQiod/6rd/qJ3tjp9ugRm5Hjx7tOhojACB4m5ubWl9f182bN7W4uNiYCqUVpw/p+vq6bt26pWPHjjWWr9Vr+qPaH+nQQ7N608xdHT30ur772mu672qde2CmxSAOr0m6Lx06Kv2zR9+jtQMH9PQrr+wZnOhV25bxwAP6Dx/4wM4HX/ua9Nhj0rve1de2AADGz5EjRya+MqqvQNQ0Tf3pn/6p3ve+9+2amqRcLqtSqTQC03A4rJMnTyqVSsmyLP0P/8P/0NN6egngguTky+kr2o3XvqQAgPHUXMO5tbWlK1euNGo4W3Hm62x2YeuCnn7hab3TeKc0L73l5Zf1t4w53b77sl58/XU9cuiQjr7pTbr9/e/r5Xuv7v7xi9LMozO6OLPTxKp5sCLHLWfQoh/8YOfvQx+S5ub63hYAAPzW16i58/PzikQiunHjhu7fv69yuaxMJqNIJLJrftAbN25oY2NDlmUplUrpiSee0D//5/9cf/Znf9Z1HXfu3JFt7x28YRw4gWW3vqLO9+MaUAMAvFlfX+84Cq5XjdrQA4f0pgffJB05Ij30kN5y73WFDh/Whx99VKHDh/WWHwajB2aaEnhNun9P+j/f9F29NnNfDx9ofRs/9sgjO//47nd31vGLv+j7tgAAMIi+a0Tdw9W/733v0/ve9z79xm/8hqSdaVsuXbqkUqm0q8a0Uqk0/t8wjEYz3uPHjzcGP/r617+uUqmkdDqtTCYzyLYNTTgcliQ999xzHZdzRss1TXPoeQIADM/Nmzdbfn7r1q2e0mnUhh5+p154/gXdrt3Ri9//vh65d19vs21997XXdtWK/q35ef3/XnxJ33PNF/qydU8V43m9/cBr+tEHH9TsgQO7mufOHjig1aNHpZde2vn7uZ+Tjh/3fVsAABhEX4FoMpns2CzVCUwdV69eVbFYVKFQaARn9XpdhUJBhUKhZRqf/OQn9Y/+0T/qJ3tDd/yHN/Ru07JUq1VJUjQaHXqeAADtDdoncnFxUVtbW3s+P3bsmOc03LWh9+7e05/85z/R/R8ORvSipG+/+FJj2Rdfv6fvvvqq/tb8vB5orhZ9TbLvSc8feE2PHnxA1/7W39L67du69eKLO4MTHT2q5R/5EenP/1x6z3ukX/91yVVz6se2AAAwqL4CUXeQ6cXJkyd18uRJra+v6/nnn1epVFKxWFSpVNoTzJmmqUwmo09+8pP9ZC0QoVBIhmF0DUSd7yORSBDZAgC04EefyNXVVV25cmXPKLirq6ue8/F//eX/padfeFqv3ntVN757Q/cfbTNvyw/dl3TjkKXXHrB3/sftgHR/Rnrx4D0tv+UtuhAKvfHd669Lf/mX0jveIf3ar0nvfKfv2wIAwKD66iM6iMOHD+uTn/ykPvvZz+r27duq1+sql8uqVqu6f/++bt++PTZBaKc+oIlEQpVKpeMypVJJiUTC/4wBADzr1Cdyc3NTKysrWlpa0srKijY3N1um4YyCe+bMGS0tLenMmTP63Oc+p/X19a6/dXx04aP6+YWf199999/Vw998WLqtrn+z1YP6ib+a3fvdX0oPvD6jN98/uHslL720UxP6trdJyaT08Y972hYGKgIABG3GHtcRgXxmWZbm5uYUCoVULpe7Lh+Px1UoFBSLxZTP51um9/jjjyuRSLTsy7qxsaG1tTXduXPH02BFN2/e1Hvf+1499dRTTN8CAD5aWlpq2RT13e9+t771rW+1neuzk27zhHazsrKiixcvdl3uzNveptV3v1sn/uRP9vQDjf2P/4ceOnBA55+6vDMy7ne/uxOIvutdOzWhv/RLXdMHAOwP4xhrBF4jOiqlUknSTnPZbqPdSmr0XW3Xh9UwDF29erVlP9dSqaS1tTVdvXqVEXMBYMTa3XBfeeWVvkePHXTk2dXVVc3OznZcZvbAAa0eOKDlZ5/VtWPHdObIES098ojOvOMdunbihB6akXTv9Z15Qr/2NWl2dqcGNJslCAUAjL2++ohOinA4LMuytL293Qg+nZpRZyTbWCzWskYzk8lobW1N586da5u+U7u6tramYrEowzAa6/FaEwoAGK52fSIffPDBlsv/h//wH7S0tNRxUKNBR551mseur6/r1q1bOnbsmH75l39Zn//853f+/6d/Wqt/5+9o+fZt6ctf1vKLL+rCo49Kr7wi2bb07W/v9AWdObDTF/RDH9qZouX48V0DEwEAMK72TdPccTeO1eUAMC2cUXOdoG91dVXr6+tdm8e2a27brmntmTNndOHCBV/zrnr9jVrPZ56RXntNeuABnf3BUenBB3X+v/2ANDfn7zoBAFNlHGMNAtExMY4nBwBMs1b9PFtpFVwO2ke0VV56nV7m7L/+kiTp/K99uOf1AQD2l3GMNWi/AwDYl5pHjz18+HDL5Vo1t/Vz5FknqL148aK2trZ08eJFnThxousovAAATDICUQDAvuNM2/L3//7flyT9zu/8jn7hF36h5bLHjh1r+fny8rIuXLigzc1NXbhwoe/pTwYd+AgAgEk01YMVAQDgcJq/3rhxQ3fu3NG9e/ckSVtbW7py5Yo+97nPtRzUaHV1daj5GnTgIwAAJhGBKABg6nXrD3r37l19/vOf3zOSrZe+moNaXFxsOc9pu5pYAACmAYEoAGDqtWr+2uzWrVuN5rZBaje9zLBrYgEAGCX6iAIApl675q9uo6qBdAY++uhHP6rDhw/r8OHD+vCHGQkXADDdCEQBAFOv21D1nWognYGNlpaWtLKyMrTRbL/85S/r+eef1/PPP68vfOELCofDunz58lDWBQDAqBGIAgCm3urqqmZnZ3d9dvDgQb373e/uOPVKu6lVLl++7Dk49RLItmo6fO/ePf3Kr/wK07gAAKYSfUQBAFPPaf7a60BE7aZW+ZVf+ZU9o+62CmabB0lyj9D7+c9/Xjdv3tTi4qLK5XLL9d+7d0/r6+uB91sFAGDYCEQBAPtCPwMRtetb6gShDmfez+b0vQayBw8ebJsHpnEBAEwjmuYCAPYlL01mu/UtdWsVMHoNZO/du6eZmZmWyzKNCwBgGhGIAgD2nXZ9P5uD0XZ9S1v5xje+sSeg7SWQfetb37on7ZmZGdXrdfqJAgCmDoEoAGDfaddkdn19fddnTt/SM2fOaGlpSWfOnNHv/d7v7QlOJen555/fE9C2CmTb1Xw++OCDKpfL+tjHPtZYxrZtffGLX2wZJAMAMMkIRAEA+067JrOd+mPati1Jes973tMITg3D2LOcO6BtFcj++I//eMv0H374YS0vL8swjMa6WqUJAMA0YLAiAMC+s7i4qK2trT2fN/fHvHz5ctsRci9cuKClpSVZlrUnHXdA2zxI0srKir75zW/u+U0oFJLUX5AMAMCkoUYUADBWvAwiNKhWTWZnZ2e1urq6Kx/uINRx9+5dffzjH9fS0pKef/75lul3GmCo27rb9Stl0CIAwDQhEAUAjA2vgwj54UMf+pAOHz4swzD0sY99bM88oOvr63uCUMc3v/lNbW1ttazZbA5om7Vqrutet5cgGQCASUfTXADA2Og0iFCvc4C24wS77vV86Utf2rNcuyayrbzzne+UYRg6duyYVldXdwW0rXSa09QJVNfX13Xr1i3PaQIAMEkIRAEAYyOI/pFeg912/Uhbef7553X48GHf8tgpUAUAYBrQNBcAMDaC6B/pNdjtZeqV559/fuhNiQEAmCYEogCAsRFE/0ivwW6rvpwXL15sOYeoG1OtAADQHU1zAQBjI4j+kaurq7py5cqu5rntgt1WTWTf8573NPL3jW98o+XIuUy1AgBAZwSiAICxMuz+kYMGu+78rays6OLFi3uWYaoVAAA6IxAFAOw7fgW7vdSuAgCAN9BHFAAADzY3N7WysqKlpSWtrKxoc3Oz65ygvaQFAMB+Qo0oAABdNM89urW1pStXrjSCzl5qV7ulBQDAfkCNKABgX/NSO9lp7tFe+ZkWAACTikAUADAx/G7S6tROXrx4sTEP6Pve9z79wi/8wq60vc496oWfaQEAMKkIRAEAE6FV0HjixImBgtFWtZO2besLX/jCrrS9zj3qhZ9pAQAwqQhEAQATYRhNWtvVTjanvbq6qtnZ2V3f9zs6rp9pAQAwqQhEAQATYRhNWtvVTjan3e/ouK34mRYAAJOKUXMBABNhcXFRW1tbez4fpElrq3lA26Xt19yjfqcFAMAkokYUADARhtGk1amd/OhHP6qZmRlf0wYAAO1RIwoAmAhO0Li+vq5bt27p2LFjWl1dHahJ6+bmptbX1/Xtb39bH/3oR2Xbtr7zne/4kjYAAGiPQBQAMDH8bNLqjMLrNMvd2trS7Ows/TUBAAgATXMBAGPH7/lCW6X78Y9/3PdReAEAgDfUiAIAxkqrmsorV64MXFPZnG477UbhdZrx3rx5U4uLizTdBQBgANSIAgDGxubm5tBqKlvNQ9pKq1F4nSD24sWL2tra0sWLF3XixAnfamoBANhvCEQBAGPBCfa++c1vtvx+kPlCpfbzkLq1Gym3VRBLM14AAPpHIAoAGAvdaiwHmS9U2pmHtJV3vvOdWlpa0pkzZ9o2/20XxA4aHAMAsF8RiAIAxkKnGstuc3p6Gdyo3Tykv//7v6/NzU1duHChbZ/PdkHsoMExAAD7FYEoAGAsdKqx7DRQ0ebmpj7wgQ/s6r/5gQ98YE8w6sxDeubMma41oM3aBbGdgmMAANAegSgAYCx0qrHsFCym02m9+uqruz579dVXlU6n9yzrzEParQa01e/6DWIBAMBeTN8CABgLTrC3vr6uW7du6dixY56mSPnKV77S0+eD5O/ChQu+pgkAwH5FIAoAGBt+BnszMzO+pAMAAPxH01wAwEQ7ceJEy88/+MEPBpwTAADgFYEoAGCira+v66GHHtr12UMPPcQcnwAAjDECUQDARFteXtZXv/rVXQMJffWrX2UgIQAAxhh9RAEAQ7G5uan19XXdvHlTi4uLngYe6hcDCQEAMFkIRAEAvtvc3NSJEyd09+5dSdLW1pauXLnClCcAAEASTXMBAEOwvr7eCEIdd+/epd8mAACQRI0oAGAIbt682fLzW7du+ZK+u9nvO97xDs3MzOjpp58eehNgAADgDwJRAIDvFhcXtbW1tefzY8eODZx2q2a/DpoAAwAwGWiaCwDw3erqqmZnZ3d9Njs7q9XV1YHTbtXs140mwAAAjD8CUQCA75aXl3Xt2rVdU6o011Jubm5qZWVFS0tLWllZ0ebmpqe02zX7dfOrCTAAABgOAlEAwFAsLy9rdXVVx44d082bN7W+vt4INp3mtRcvXtTW1pYuXryoEydOeApGFxcXuy7jRxNgAAAwPASiAICh6BRsDjKqbqtmv25+NQEGAADDw2BFAICh6BRsDjKqrtPsd319Xbdu3dLb3/52zczM6Nvf/raOHTvGqLkAAEwAAlEAwFB0CjYHHVV3eXlZFy5cGCh/AABgdGiaCwAYinZ9OZ1ay2GNqgsAAMYfgSgAYCg6BZteRtUFAADTi6a5AIChaO7L2dx/k+a1AADsXwSiAIChIdgEAACt0DQXAAAAABAoAlEAAAAAQKAIRAEAAAAAgSIQBQAAAAAEikAUAAAAABAoAlEAAAAAQKAIRAEAAAAAgSIQBQAAAAAEikAUAAAAABAoAlEAAAAAQKAIRAEAI7O5uamVlRUtLS1pZWVFm5ubo84SAAAIwKFRZwAAsD9tbm7qxIkTunv3riRpa2tLV65c0bVr17S8vDzi3AEAgGGiRhQAMBLr6+uNINRx9+5dra+vjyhHAAAgKASiAICRuHnzZsvPb926FXBOAABA0AhEAQAjsbi42PLzY8eOBZwTAAAQNAJRAMBIrK6uanZ2dtdns7OzWl1dHVGOAABAUAhEAQAjsby8rGvXrunMmTNaWlrSmTNnGKgIAIB9glFzAQAjs7y8rAsXLow6GwAAIGDUiAIAAAAAAkUgCgAAAAAIFIEoAAAAACBQBKIAAAAAgEARiAIAAAAAAkUgCgAAAAAIFIEoAAAAACBQBKIAAAAAgEARiAIAAAAAAkUgCgAAAAAIFIEoAAAAACBQBKIAAAAAgEARiAIAAAAAAkUgCgAAAAAIFIEoAAAAACBQBKIAAAAAgEARiAIAAAAAAkUgCgAAAAAIFIEoAAAAACBQBKIAAAAAgEARiAIAAAAAAkUgCgAAAAAIFIEoAAAAACBQBKIAAAAAgEARiAIAAAAAAkUgCgAAAAAIFIEoAAAAACBQBKIAAAAAgEARiAIAAAAAAnVo1BkIgmVZWltbk2VZMgxDlmVJkjKZjAzD8CX9dDqtWq2m7e1tSZJpmjp37pxCodDA6QMAAADANJn6QLRWqykajSqbzSoSiTQ+r1Qqevzxx3X16tWBgsVCoaBsNqtMJrMrnUqlorNnzyoSiSiTyQy0DQAAAAAwTaa+aW40GlUymdwVhEpSKBRSJpPRyZMnGzWkvapUKspmsyoWi3uC2VAopHK5rFKppI2NjX6zDwAAAABTZ6oD0Y2NDdVqNSUSiZbfO5+n0+m+0k+n011/m8lklM1m+0ofAAAAAKbRVAeia2trikQiHfuBRiIRXb58ua/0b9y4IdM0Oy4TiURUq9X6rnUFAAAAgGkztYFopVKRZVldA0XTNGVZlkqlUs/r8PI7JwD1Y1AkAAAAAJgGUxuI3rhxQ5K0sLDQcTnn+0ql0vM6TNNUOp3u+NvLly8rFov1nDYAAAAATKupDUSr1aqk7jWR8/PzkqTr16/3vI50Oi3LshQOh1sOSFSr1ZTJZHT+/Pme0wYAAACAaTW1gajXPplOoNpPH85EIrFrwKNwONyoHa1UKspkMiqXyzTLBQAAAACXqZ1HdHt7W9IbNZ7d1Gq1vtaTzWYVDoeVTCZVqVQUDocViUQUj8f7Gi339u3bPS3/2GOP6ciRIz2vBwAAAMBkevbZZ/W9733P8/K9xhhBmNpANMhRap1a0WQyKUkqlUra3t5WJBLpOlhSs0984hM9Lf/pT39an/nMZ3r6DQAAAIDJ9du//dv6zd/8zVFnYyBTG4gGKZlMamFhQbZta2NjozGA0cLCgrLZbNt5TFt58skndfToUc/LP/bYY/1kGQAAAMCE+tSnPqV4PO55+du3b/dc4TVsBKI/1E8/Tmegomw2q0gkIklKpVKKxWKKx+OqVCpKJpOyLEupVMpTmkePHtXi4mLPeQEAAACwPxw5cmTiu+dN7WBFTmDp9BXtxmtfUreTJ08qk8k0glCHaZoql8vKZDKS1HWKFwAAAADYT6Y2EHUCy259RZ3ve60RdaZr6TRHaCqVUrFYlLQTjAIAAAAApjgQDYfDkqTnnnuu43LOaLm9Dip06dIlnT59uutykUhEiURCN27c6Cl9AAAAAJhWUxuIHj9+XFL3aVmq1aokKRqN9pR+rVbzHLxGo9FAR/EFAAAAgHE2tYFoKBSSYRhdA1Hn++Z+nt2Ypqnr16/3lB8AAAAAwBQHotLO/J6VSqVjbWSpVOo4vUq7354+fVqFQsFTPi5dutSYYxQAAAAA9rupnr7l3LlzyuVyWltba4xg67axsSHDMFp+J0nxeFyFQkGxWEz5fH7Xd6lUqhFgZrPZtnkoFAqyLKunuUT99vrrr+uFF17QSy+9pHv37sm27ZHlBQDgj+jCQ5K6d0EBAEyHmZkZHTx4UG9+85v1lre8RYcOTXYoN9U1ooZh6OrVqyoUCntqL0ulktbW1nT16tW2I+Y6v2lX83n16lXduHFD0Wh0z/QstVpNyWRSly5d2hPEBuX+/ft6+umn9bWvfU3PPPOMXn75Zd27d28keQEA+Ounf2xOP/1jc6POBgAgIPfu3dPLL7+sZ555Rl/72tf09NNP6/79+6POVt8mO4z2IBQKqVwua21tTcViUYZhNJrb3rlzp+O0LZlMRmtrazp37lzL7w3DULlcbgS17rfSpmkqmUz23PfUL7Zt6+mnn9ZLL72kH/mRH9Hc3Jze/OY368CBqX73AAD7xteffVGS9JNHHhlxTgAAQbl//75eeukl1et1vfDCC7p//75+7Md+TDMzM6POWs9mbNppjoWbN2/qve99r5566iktLi4OnN4zzzyj7e1tHT58WG9/+9sn8uQEALRHIAoA+5dt2/rOd76j559/XvPz83rrW9/acXm/Yw0/UD02he7fv696va6HH36YIBQAAACYMjMzM3r729+uhx56SPV6fSKb6BKITqGXX35Ztm3rLW95C0EoAAAAMIVmZmZ0+PBh2batl19+edTZ6RmB6BR66aWXJElvfvObR5wTAAAAAMPiPO87z/+ThEB0Cr366qs6cOCAHnzwwVFnBQAAAMCQPPjggzpw4IBeffXVUWelZwSiU+jevXs6ePAgzXIBAACAKebMLUofUYwNpmkBAAAApt/MzIwmcSIUohUAAAAAmFCT2gqSQBQAAAAAECgCUQAAAABAoAhEAQAAAACBIhAFAAAAAASKQBQAAAAAECgCUQAAAABAoA6NOgPAHtvb0te+Jt2+LT3zjPTaa9IDD0hvfat09Kj0rndJ8/OjziUAAACAPhGIYjzcvy9dvy794R9Kf/zH0osvSt//vvTKK5JtSzMz0sMPS296k/TII9LP/qz0i78oPfGEdICKfQD+sCxL6XRa2Wx21FnZ1zgOb2BfYJJx/qITAlGM3je+If3Lfyl95SvSs89KDz0kveUt0tvethN8HjiwE6i+8spOcPqd70i/+7vSF78offCD0q//uvQTPzHqrQAw4SzLUjweVz6fH3VW9j3DMBSPxxWNRlUsFkednZHhnMSkoyyjE6qSMDq2LV25IiWT0h/8gXT37k6z26NHpSNHdmo/Dx7cqQ09eHDn/48ceaN57t27O79LJnfSse1RbxH2oY2NDUWjUYXDYc3NzSmdTg91fZVKRfF4XOFwWAsLC5qbmxvq+vaTkydPKpPJyDAMT8tblqWZmZmWf4VCYaC8FAqFtmmPSi6X23XuzczMaGNjY2jri0QiisfjSiaTQ1vHuOv1nATGEWUZ7RCIYjRsW/qd35E2NqS/+Avpne/c+XvwQW+/f/DBN37zF3+xk87v/i7B6JA5AVetVpuK9fghFAopGo3KsixZljX09RmGoWg0KtM0J2L/TIpkMqlIJKJQKOT5N4ZhyLZt1et15fP5XcHC9evX+86LZVlaW1tr/H8kElG5XFa9Xpc9wmvc/Py8TNOUpMa518v+6kcikVCtVlMulxvqesZRL+ek84LKediPx+NKp9NDuSY561pYWGj89bouP9JwFAoFRaPRRjrhcHioL0haSafTQw2y/Di+fqThfvHq7O94PK5SqdT1t/u5LKMDG2PhqaeesiXZTz311MBpVatVu1qt+pCrIXrySdv+mZ+xbdO07Y99zLb/zt/p/+9jH9tJ52d+ZiddDEWxWLQl2ZLsVCo18evxWzabDTzPkmzDMAJZV7VatfP5fCDrClq5XB54P5bLZTsWi9mmadqS7Fgs1ndasVjMzmQyjXJQLpdbLnfnmRfsO8+80Pd6+pXP5xt5C0K1WrUNw7Dr9Xog6xsHvZyTqVTKNk1zz30/m83ahmG0PX/6kUqlbMMwdl0L6vW6nUgkPK/LjzRse+e8CIVCdiQS2fObbDZrRyIRj1vVu3q9bheLRTuTyfhS5jvx4/gOmkaxWLRDoZCdz+d3lcN6vW6nUilbkh0KhbrmYz+W5aB4efb3M9bwC31EEbxvfEP61/9a+va3pZ/6KenQgKfhoUPSu98t/fmf76T7N/8mfUaHIBKJKBKJqFarDfXNb1Dr8ZtTUzStarXa1NbAnj17VplMZqA0SqVSo2Z8kH3l1O64+1INu9axV07egsqXaZo6derUvhrwxOs5mcvllMvldOfOnT3NdxOJhCzL0smTJ1t+36t4PK5CoaBqtbrremcYhrLZrLa3t7uuy480pJ3avZMnT+rUqVMtzwnTNLW9va1CoaBYLNb3NrcSDodlWZZM01QymRxqbagfx3fQNCqVirLZrMrl8p7vDMNQJpORZVnK5XIKh8Mtl3Psx7KMzmiai2Ddv78zMNHXvib95E8OHoQ6Dh3aSe9rX5P+1/91Zz3wXbFY3PMAMcnrmXRB9hub1iC0UCioVqspkUgMlE6xWFQkEtnTdLVXly5dUiKRaDR1i0QiA+VrGEaRt3Q6rVwuN7XnoZvXc9KyLCWTSZ06darttSCVSknaCQAHUSqVVCgUlEql2l6XnYDk7NmzQ0tDUiNwOn78eNtgJhqNqlKpDNREvp1yuaxqtapisahYLKb5IU0n58fx9SONtbU1VSqVjv3end9WKpWuzXT3U1lGdwSiCNb16zuj4775zTt/fnLSvHZNunHD37SBfW5aRzvMZrMDB6HSTuBpmqYWFhYkqa++bslkUplMRrVarfH7aDQ6cN785NT4SsHmzTRNRSKRfVGL4vWcdPoRd6uNO3XqlEql0kAP/k6g0WldpmkqFAqpUCi0PP/9SEPaqS22LKvjuZBIJGQYhk6fPt12mXHnx/H1I41KpaJardZxIL7jx4/vWr6T/VSW0R2BKIL1h3+4M0XL2942nPTf9rad9P/wD4eTPrBPDToK7Diq1WoqlUoDP6xWKpVGM1V3TU8vD/6VSkWGYcg0zV01CuNWIzrKvMXj8akf6KSXc9LZF92aSDsvR/p98C+VSo2gsFsrFScguXz5su9pSGrUzLlbH7SSzWZVr9fHrll7L/w4vn6kkUwmZRhGx2B2e3u78W8v+3w/lGV4QyCK4GxvS3/8xzvzhHodHbdXDz648/flL0v1+nDWAewzw56SZlQKhYIMwxj4YdXpHyr1H4im0+lGn8BJ6B86igD51KlTsizL0widk8rrOemM0u2leb5zTnarqWrH+Z2XdYXDYUl7W1D4kYb0RqA0aFPjcefH8fXrHEmlUqrX640mvK24r3Xu2tF29kNZhjcEogjO174mvfii9Ja3DHc9hw/vrOdrXxvuejAw+oiMv1wuF/hUCEG5dOmSLwGV0z9U6i8Q3djY2BXsT0L/0FE0GXZqjPP5fODrDorXc9I5Dl760TvL3Oizy0q1WpUkT30hnWWaAxo/0pDeqOHzEuxMMj+Ob5DniFMmY7GYp8B3P5RleMOouQjO7dvS978/vGa5jje9Sfrud3cC0fe/f7jrGrFCoaBLly41btzb29uan59XOp1ue/NJp9ONPh/b29s6f/68YrGYKpWKLl26JMuydOPGDWUymcYDUS6XUz6f1/b2dmMQjU4jOjrzIDpNsba3t/XEE08olUrJsiyl02kZhqFHH31Ua2trjdH6vKzHyb+zzLlz5xrpOuvc3t5uNJfMZDKebsTuuRudNLrty2GrVCqNPM3Pzzdu8OfOnetpoKJ+ti2Xyymbze56GEyn03tqR1OpVNtzYdB9msvllE6nG8fYb5VKxZc+ZE7/UGnnAcswDFmW1Xj47vbbarXa2L5J6R/qjGydTqcbn5mmqXPnzg21FjcUCvX9sDwJvJ6TTlNIL4Gdc63opYasX07ag7xkbJeG+1rknGOWZeny5cuqVqt69NFHe54LeFz5cXyDOkecEXMl9TT6+LSXZXhDIIrgPPOM9Mor0sMPD3c9Dz+8s55nnhnuekbIGTXQNE2dP39+102jUqkoHA63fXh/4oknGgGg88BbKBR0/fp1ZTIZVSoV5XI5RaNR2bYt6Y23z86Ihp2USiXF43FlMpldA244k5e7Bymo1WpaW1tTrVZTKBTytJ5oNLon/84DcfO+cCbd7jYCrxPwNf9+Y2NDCwsLymQyQwmEOonH46pUKsrn87serJwRJc+dO+cpnX63LZFINI6fU2PXy37wY58658Ha2prv+995qB30odXdP9RhmmbjZU83znnrGOf+oe4+ezdu3FA+n1c2m22UrY2NDYXD4V01xH574oknprK/stTbOdnLYFjuQGR7e7vnQNTpP+juB9iO+5x3BzR+pNGqxq9YLCqZTCqRSKhWqzUG+8rn84GOKu43P45vUOeIM8JxsVjs6aXtNJdleEfTXATntdck25YODPm0O3BgZz2vvTbc9YxQOBxuNGtpvmGEQiGVy2Wtra217NsXi8WUSqUagcz29raKxWLjTaaTnjvdUCikRCLRta+gZVmKx+M6fvz4nlEfz58/r+3t7V0P2qZp7hpQwst6IpHInlq4dDrdcl84D/jd3tK2+30qlVIsFlM6nQ60L0s0GlWpVFK5XN7zUGoYhvL5vNLptKcHjVFtmx/rdWrOvQbdvXAeaget7Xb3D3V47ZPnzBnq3keT0D9U2pnCovnBM5VKyTCMofbfG7S/4zjr5Zx87rnn+lpHP6M5Oy8VvPzW3QrAHXT6kYbzuWEYqlQqjfuWs79M01Q2m1UoFNLjjz8+0eeIH8d3mOeI07/TmTu5Wq32/PJpmssyvCMQRXAeeECamRn+HJ/37++s54EHhrueEXGawrlrUZo5TeQ2Njba1so4D7/ZbHbXaHhOcFhvMdhTtwckJzhq9SBqGEajOV+3t6BeHsScN7jN+XdzHuRbjb7oKJVKewJkNyeIDWrAnlwup1KppEwm0/GttJcRMEe1bX6tN5FIdB0ko1/OQ+2ggWir2j8vc4laltWYM9RtEvqHdpp64fjx47ua6vlt0Hland/OzMwM5W+Qlzq9nJP9BJT9CoVCjWtpt+1zf+/Oox9pOP+en59XNptt+4LR+XySBzTy4/gO4xypVCqKx+OKx+ON+26/3Vf8KMuYfDTNRXDe+tY3ms2+6U3DW4/T/Petbx3eOkakUqloY2PD04AATs1iPB5XuVxuu5xlWS1r3frhvNFvN5CE0xTn+vXrisVifa2j2fb2dseHdqe/XjuVSkWVSmVXU2S3IN/aOv1nJXWdR9DLjX9U2zZO+7Qdvx5+3P1DHU4zRGl300I39yi57rS89g+1LEv/9H/8J5KkH3/7Y6rVajp9+rRv5aqZO29e+oE19491n9uGYfSdXz/6IJqmObR5cQd5gTDOD+Tnz59XOBxWPp9vu42VSkXHjx9vW679SEPa2U/OyLrtnDp1SrlcTrlczpd5grEjFArtGWDIec5wxpvwyo+yjMlHIIrgHD26E4B+//vDDUSd9N/1ruGtY0ScWognnnii67LOqHROUNCumZ+fzf+cG0q3QNbPG0+30RPn5+c7BqKRSESGYejUqVNtlzFNM5Cb5eXLl1u+GOjXqLZtnPZpJ4P2IWtXrppHzm1eplQqaWFhYU8A67V/qPMg/rv539d7l/6mfvLII5J2glenr7ffnLx1m1rEaUrpPrZOfq9evbrrt/3k18vAK16MY42zNPg5Oaz0Q6GQUqmUNjY2lEwm95wDlmUpm83umh+yVbeRQdJw/7vb8YtGo8rlcnvGKph2fpw/vabhroFOJBKe56v1qyxjstE0F8F517ukRx6RXnhhuOt5/vmd9UxhIOp+GPTCedDt1BTKzxFhuzW1cfqseAmke11nv0KhkOr1et+TvfvJqaXx65iMatvGaZ+240eztVb9Q6XuU7i0G6jJa//QeDyuU6dO6b1Lf3PX59lsVhsbG0Pp8+t1/lCnJsv9kOnkt3mbBslvv/3fxlkv52S/AccgD/+ZTEb5fF7xeHxX94pSqaS1tbU9LxRaXccGSePRRx/tmLabu7YtyGbMfvHj+AZ5jjjHLZfL9TwA0TSWZXhHjSiCMz8v/ezPSr/7u9IPfiA9+KD/6/jBD3b+PvQhaW7O//RHzHmo7fVG0WkaCffNfVCnT59u1MC2emB1HlLHtTZCUmMaG3eTSy8jPfrBPQ3GMIxq20a5T4epWCy2DLY7BaKtmuQ6vPQPdcpXqzRM02xMWeR3GfMyf6i7OaXTPHlU+Z12znXbSzlyLzNojVksFlMkEtHly5eVTqf16KOPNo6h9EYw3eka1m8arQbSa8d9j2zVKmHc+XF8gz5HEomEcrmczp49O7QuApg+BKII1i/+ovTFL+7M8/nOd/qf/ne/Kx05srMeeOJnU7BUKqVsNttyuo1KpaJSqaRUKjV2DwVO/7XLly8rEokomUzuejAuFAoT+VZdGt22jfs+9eO8b9U/1OE0PXa/BHJqZ1qd/177h166dKmRfrv1+j0lgjtvnQLGVk2Lh5VfP1+gjYtezklnf3opQ343gTcMo21z1+vXr0vq/rKxnzScbhiTei3uhR/HN+hzxHn55Iyo6/Xl0jSWZXhHIIpgPfGE9MEPSn/wB9JLL0lvfrN/ab/00s7fz/2c1KXf4KRyBt7xWpvkLOcePGXYnJFxo9FoY47BUqmkZDK5Z9qVceAe+KK5D1vQnD69fj1oBbFtzuBZzbWA47JPOxlkP3fqdy29EYi6H/Ca5wx189o/1FnONE19/dkX93zv1AR1y18v3H2/O9V0OUGnux+pO7+t9Jpfv2rSK5WK77Xy8/PzA+9zr+eksx4vAYSTZhDl0KkVH2TE2nZpuM+hdoOAOdzHdlgtTIbJj+PrRxq5XK7xErHbAF/duiQ0m4ZWMRgcgSiCdeCA9Ou/Lv35n0t/8RfST/2UdMiH0/D116Wvf116z3t20h/2XKUjEolEVCgUOja1dXNuBkEFArVaTYZhqFgsqlKpNPLqTHI/jg8EzjxoxWKxp/3UqTasX86own7VYASxba3694zTPm1nGPOHtkrfOZa5XE6nT59u+/DstX9ot3PD3TfOr3LvpbllrVZrBBDul01+59dLXrrxMupqv1pN5+NVL9tkmmbjxWS3oMw5BkE0fy6VSo2puvxOw3nBUalUdOPGjY7rcM4TwzCGPgDUMPhxfP1Iw+l6UCqVul6f3WXdS/chP8oyJt90Pq1jvP3ET0i/9mvSO94h/eVf7gSRg3j99Z103vGOnXSH0eR3TDjzdnmZ9sK5+ZimGVj/q0Kh0BiIyBkhMZvNKpFIjOXNxgn6uu0j5wbrrq0YxryiTlM1PwabCXLb3A8d47ZP23FaCfQb9HcLONzpO0F5p35T4zx/qFN2Oz1cOg+spmkOdZRSP/pRm6Yp27aH8jfI8ev1nHRGpXamzWrHaerabr5lL9LptObm5jpem5yRbtu1evEjDa/3QOfFTqeRu8edH8d30DScchaJRLqWOfcLci/lYNhjImAyEIhiND7+cSmZlN72tp3a0Zde6i+dl17a+f3b3raT3sc/7m8+x0wkElEsFlOpVOrahMu5oTfP+TVMhmE0bmiTwMmr17e8w2YYRuMBrFu/uW7f+7ltThqtauKdWvBhrHeYnP5m/ealW+2A+ztnjr1OafUyf6gX7Zq99TOqpVNT2S7NWq2mjY0NSdrTfG/Q/LZalztP06TXc9J5cdPtGl8oFPY0n2+Wy+W6BojOC5VO+YlEIm1fRPiRxqlTp2QYRqMZeDuXL1/edT1tlx+/+1P3ots+9+P4DprG6dOnlUgkPM2762xLIpHwVAs9zWUZ3hGIYjRmZqS///elVGqnOe03v7nz94MfePv9D37wxm/e856ddP7BP9hJd8qdP39ehmHo5MmTbZep1WpKp9NKJBJdL/J+Dvxw/PhxbWxsjPTm3k6r7fRSA5HNZhs1We5aPD+noHFzBnPqVDvozLfn/LsVP7fNWabbW3W/15vL5bSwsNAIcvzklAsvrQuabWxsdH3Qcj/QRaPRjst77R8qDTbIUjKZVDKZVDwe73mfxmKxtvvK6cvXqvm9380ir1+/PrUPrr2ek6ZpKpVKKZfLtS1vznWk04uQdDqtZDKpaDTa9rxwgsNz5861/D4ajWp+fr5jwONHGoZh6Pz5842uH604Aa9zr2xlkLLQrJ++jl72uR/Hd9A0YrGYDMPouo/S6XTj5ZzXabumuSyjBzbGwlNPPWVLsp966qmB06pWq3a1WvUhVwH5+tdt+7//7237Ax+w7ccft+2f+inbfv/7bfsjH7HtX/gF2/7bf3vnvx/5yM7n73nPznIf+IBt/+N/vPP7faZer9uRSMQOhUJ2vV7f9V2xWLRN07Sz2WzHNGKxmC3JjkQintebSqW6/sZJt9WfaZp2LBazM5nMnnz3up5EImFLskOhUNtlqtVqY93FYrFjfhOJxJ7vMpmMXSwW7XK5vGuZTCZjl8vlnvPci1gsZkcikT37qVqt2rFYzK7X641tSyQSdrFY3LOsX9vm3r58Pr8njVZ592O9pmnakmzDMFrvpAGFQiE7Fot5WrZer9vlcrmxHyTZ2Wy27bXWOT6dzk8nTWc7nfO0U9lwlrVt277zzAv2nWde2PV9q+PkiEQijfX0uk/r9bptGIadSqUan1WrVTsUCtmmabY8Z5rz20qn/LZLr9V5NS16OScdqVTKNk1zz7mYzWZtwzDaHhv3Ot3XklacY91c3svlsh0KhVpeq4aRhiObzTbKoFsmk7ENw+h6Pg1SFtyce7GTjtdnLy/73DHo8fUjjVQqZUciETufz+86RvV6vXEvdu5LXk17WQ6al2d/P2MNvxCIjol9HYjatm3fu2fbX/2qbX/60zsB5/Hjtv3TP70TcP7kT+7896d/eufzj3xkZ7mvfnXnd/tYPp+3I5GIHYvFGoFLIpHoePzdD73uP8Mw7Ewm0/I37pu2+6856KrX640bXiwWsxOJhJ1IJBr5C4VCtmEYjd83r8/Letw38Ob8O4rFYttguNU25vP5Rv6cPKdSqV37MZ/P26FQaM+DVLs8+xE8uY9vKpWyU6nUrgcvZz2mabZ8wBt025pls1nbNM3GedbpYc+P9ToPSe3Oy0E5D63duIN+9/HtdE7Ztt02QHOCcHdazp/781YPdc75X6/XWwaizkNhq/UWi8XGekzT7LrdrfaDc51xXoR1e+Hlzm8rnfLbav2dXipNA6/nZLNisdi4zjrX3FQq5SkwcF5etgpU3KrVqp1IJBrLOtcCry8R/ErDUS6X7Vgs1jifncDG6zb3Wxacdba7x7jLdKf1e9nn7uX7Pb5+peHc352y79x3nBehvdgPZTloBKIYyL4PRN22t3eCzH//7237X/wL215f3/nvv//3O59vb486h2ihXC57ehNt2zs3ROcBlRsRRsWpMfcSBI0LpwaxWq22DESdFyPdHi57rXXrlzu/rXjNr22/8WJimk3iOTnpgioLeMN+KMtBm9RAlD6iGD9zc9L73y/9vb8n/ZN/IqXTO//9e39v5/O5uVHnEC2cPHlSiUSi48igjkgkonK5LNM0x25eUewfzsi+3QY+GSenT5+W1L4fYa1WUyQSGZspK/zMbz6fH+qIvONgEs9JoFf7oSzDGwJRAAOrVCqyLKvriJ/NIpHIWIygiv0rmUw2RpieBKFQSKZptgxUarVaY6CyceFXfmu1mkql0kBTkEyKSTsngV7sp7KM7ghEAQzMPSl9L27cuMEcYhgpZ8qCSXrwLxaLKpVK+pMv/cddn8fjcSUSia4j7xYKhZ5fGg3CyW/zVBVe8yvtzCs5rvMR+20Sz8lJFXRZwP4qy+huxrZte9SZgHTz5k29973v1VNPPaXFxcWB0mKSYIxCOp1WLpdrNLntZmNjQ2tra7p69SpDuGOkKpWKTp48qXq9PuqseGZZlj71//knkqS/+d6f0vXr1xWNRj01d4tGo57mBfSTZVmNms+FhYWe8lur1RQOh3Xnzp2xaXI8bJN4Tk6iUZSF/Ww/luWgeHn29zPW8MuhUWcAwHTIZDJaWFhQNBpVMplsO6l1pVLR2tqaarWa56AVGKZQKKRTp04pnU5PTJ9lwzD0v/zzfyVJ+skjj3j+XTqdbsz7GSTDMDzPL9gsmUx2nBNyGk3iOTlpRlUW9rP9WJbRGYEoAN8kEgklEgkVCgWdPXtWkjQ/P9/4fnt7W/Pz8zp37hy1oBgr2WxW4XBYlUplas/NWq0my7ImapCQXC4n0zQ9DYI2bfbDOTkqk1gWJt1+Lstoj0AUgO9isRg3G0ycq1evKh6PK5/PT+Ube9M0+66VHIVKpaJisah8Pj/qrIzMtJ+TozJpZWHSUZbRDoMVAQCgneaj+Xx+rEad3a8sy1I2m933D66ck5h0lGV0QiAKAMAPDdKXEf7hOLyBfYFJxvmLTghEAQAAAACBIhAFAAAAAASKQBQAAAAAECgCUQAAAABAoAhEAQAAAACBIhAFAAAAAASKQBQAAAAAECgCUQAAAABAoAhEAQAAAACBIhAFAAAAAASKQBQAAAAAECgCUQAAAABAoAhEAQAAAACBIhAFAAAAAASKQBQAAAAAECgCUQAAAABAoAhEAQAAAACBIhAFAAAAAASKQBTA2Emn04pGo1pYWNDc3JxyudyoswQAAAAfHRp1BgD0L51Oy7IsZbPZnn5XqVS0trYmSZqfn9f29rZM09S5c+dkGEZgabTzxBNP6NFHH9Xa2posyxooLQAAAIwfAlFggliWpRs3bqhSqSibzapWqykWi/WURjqdVqFQULFYlGmajc9zuZwef/xxXb16VaFQaOhpdOLepnQ63Xc6465Wq6lSqfR8DMd1PQAAAF4RiAITIhwOy7IsmaapZDKpdDqtZDLZUxq5XE65XE537tzZU2uZSCRkWZZOnjzZ8ns/0/Bq0N+Pu1qtplqtNjXrAQAA8Io+osCEKJfLqlarKhaLisVimp+f7+n3lmUpmUzq1KlTbQO8VColSYrH40NLA28IKjgkCAUAAOOGQBTYJ5z+nN1qUU+dOqVSqdQyePEjDbyhWCxO1XoAAAC8IhAF9gln5NlufTcXFhYkqeUASH6kgTcUCoWpWg8AAIBXBKLAPmBZlizL8tTn0hl8qFKp+J4G3hDUIEzTPNgTAACYXASiGFubm5taWVnR0tKSVlZWtLm5OeosTaxSqSRJu0a4bcdZ5saNG76ngR25XE4bGxtTsx4AAIBeMWouxtLm5qZOnDihu3fvSpK2trZ05coVXbt2TcvLyyPO3eTZ3t6WJE8DHDk1ns01oH6kMahKpaJLly410q7VajJNU5lMpqd1FAoFXbp0qbEt29vbmp+fVzqd9hRo95tOLpdTNpvdVVOcTqf31FqmUillMpk9v9/Y2NBzzz2nRx99VM8991zj83PnzmltbU3JZFKmafa9nnQ6rUqlolqtpu3tbZ0/f16xWKyx353pgzKZjCKRyJ78WZbV6EdsWZbn/dq83kwmo0Qi0ZimyNm3tVpNx48f93ScvO4rAAAwIjbGwlNPPWVLsp966qmB06pWq3a1WvUhV6Nz5swZW9KevzNnzow6a2Mjn8/bkuxYLNZ12UwmY0uyI5FI12Xr9Xpjf7vPIz/S6FU2m7Ul2dls1s7n83Y2m225jGEYdrFY9JSvUChkx2Ixu16v7/quXC7bhmHYmUwmkHSc/dnL+lpto/OdJLtcLg+0nnw+b2cyGdswDFuSnc/n7Xw+b6dSqca2Oce1Wblcbrk/vKy/eb2ZTKbx1ywWi3U83oPsq0lz55kX7DvPvDDqbAAARszLs7+fsYZfqBHFWLp582bLz2/duhVwTqaDu0aoF5Zl+ZpGv7LZrJLJpBKJxJ7vnM+i0ajy+bxisVjbdMLhsEKhkPL5/J7vQqGQyuWywuGwnnvuuZY1kn6n49XZs2eVTCZb1kIahqF8Pt8YIGoQ7n2XTqe1vb2tcrncqJV0ap1b1T6n0+mWo/OmUildv35d6XRaoVCo5TY0rzebzSqTybQ8lvl8XtFotO3xDmpfAQCAwdBHFGNpcXGx5efHjh0LOCfTwY9g0I80+jU/P98yCHUkEgmZpqmzZ8+2XSadTqtWq+n8+fNtlzFNU+fOndPGxkbbqWf8SqcXhUKhYzNS0zRbBl79cgJN5wWAez31el31en3X8qVSSdvb241+xM2cYLzbwEnOek3T7PhCwQmMz549u+e8DHpfAQCA/hCIYiytrq5qdnZ212ezs7NaXV0dUY4wSt2mi5F25ja1LKtlsFOpVLSxsaFYLNa1L6kT8Mbj8aGl049Wta9uXvZRryzL2pNuq+2uVCqqVCqKRqMt0+l1FOVu2+IEqu2O9yj2FQAA6A2BKMbS8vKyrl27pjNnzmhpaUlnzpxhoKIR8GOQIb8GKurGqeVy5jp1c2rQnnjiia7pGIYh0zQbwdUw0ulVKBRSLpdTMplsW8OaTCZ9D7C8pheJRGQYRsdaa78HBjp9+rQk6fLly7s+H9W+AgAAvSEQxdhaXl7WhQsXtLm5qQsXLhCEDqDfYNA9Qq4faQyTE+g4o+m6OU1GvW6Dk1ZzU1O/0ulVPp+XYRjK5XJaWFjQ3NycksmkCoVCo2nqMEaA9ZpmKBRSvV5vBOpBaHe8R7WvAABAbwhEgX3g0UcflfTGFCyduJdxB1x+pDFM7vU0B6LO//caFFer1aGk0yvTNHXnzp1GjaNlWcrlcorH45qbm1M0Gh1KH17nmPeqUqkonU4rHo83po3xct70wh1Muo/3qPYVAADoDYEosA+4a4+6adec0Y800D/DMJTNZmXbtorF4q5RZUulksLhsO8BVi8vESzLUjKZ1NzcnNbW1hqj2mYyGWUymcBqxqXR7CsAANAbAlFgH3D6w3kJEJ0H9OY+dH6kMUzNtWJuTkDltVbOWa55mg+/0ummebTd5qa9kUhEqVRK+Xxe9XpdiURCtVqt46jBXtbTr1qtpscff1yXL1/W1atXlc/nhz4ybbvjPax9BQAA/EUgCuwDpmk2gqhuNUHOA35zIOFHGsPkrNMZJMjNyYfXJrJOWs2BtF/pdOOes7XdyLAOp/YvEomoUCj0vZ5BOM1d8/l8T9s6SBDc6ngPc18BAAB/EYgC+8SpU6ckSTdu3Oi43PXr1yVp1/yRfqYxLMViUZJajtzq5MPL6LWWZcmyrJbzTfqVjhfupqyVSqVr8N/vvh60yWyhUFCtVuu6nU7g6N6ObvOKduIMjNR8vIe5rwAAgH8IRIF9wnno7zbHYqFQUCwWazmyqB9p9MNL4JfL5WSapjKZzJ7vIpGIYrGYSqVS1yDFmf6l1Tb6lY70RnPSVrWrtVptT//M5mlKmrWqCe5nPb1yXjp0Ota91nx2O961Wk2lUqnt8e53XwEAgOAQiAITqtdRSE3TVCqVUi6XaxsYOIHm+fPnh5ZGP0zT7NiUMp1ON5qGtnP+/HkZhqGTJ0+2XaZWqymdTiuRSLRtYupXOs7gOd1qlx3dag/z+XwjzUHW4+ZlQB+n/2unYDObzTby4a4ZbTcf640bNzpOeROPxyW1D/L73VcAACA4BKLABHIHXaVSyXONUyaTUSqVUjQa3fObXC6nXC6ncrncsZbMjzR6EYlElM1mdf369ZbBaDqdVqFQULlc7tg/0TAM3blzR/Pz8y1HTS2VSopGo8pmsx3nw/QrHUlKpVKqVCq7tmtjY6Nl09FMJqN4PN4yOHQGHWpVO9jreqQ3mjk7/+0kkUgoFoupVqu1TG9jY0PRaFTnzp2T9EaT2lwu17YpbyKRULFY3HO8Lctq7PNOx3uQfQUAAIJxaNQZAOBNPB5v2STUsqxGrZQT/NXr9bbpZDIZRaPRRq3R/Py8tre3G/Mvegkg/UjDC9M0G+vIZDIqFApKJpMyDEOWZalWqykajXoePMgwjEaAE4/Hdw2+ZJqmisWipyabfqWTyWS0sLCgdDqtbDYr0zQVjUb3BGiJREKJREKnTp3S2tqaarXarr6d0Wi0Y9DodT0LCwt7RuudmZmRYRg6d+6cUqlUy/Tz+bwKhYIuXbqkcDis48ePN/ZTMpls7It8Pq+1tTWFw2FlMpmOLw6aj3etVlOtVtPp06fb5sOPfQUAAIIxY9u2PepMQLp586be+9736qmnntLi4uJAaTkPkvSBAjBJcrmcksmkUqkUNZYefP3ZFyVJP3nkkRHnBAAwSl6e/f2MNfxC01wAAAAAQKAIRAEAAAAAgSIQBQAAAAAEikAUAAAAABAoAlEAwFhwRrOtVCojzgkAABg2pm8BAIxUp2ljOk1FBAAAJheBKABgpLzOAwsAAKYHTXMBAAAAAIEiEAUAAAAABIpAFAAAAAAQKAJRAAAAAECgCEQBAAAAAIEiEAUAAAAABIpAFAAAAAAQKAJRAAAAAECgCEQBAAAAAIE6NOoMBMGyLK2trcmyLBmGIcuyJEmZTEaGYfi2nlKppEwmo1qt1kj39OnTSqVSvq0DAAAAACbd1AeitVpN0WhU2WxWkUik8XmlUtHjjz+uq1evKhQKDbyeeDyuWq2m8+fP70qvVCopmUwqm80OvA4AAAAAmAZT3zQ3Go0qmUzuCkIlKRQKKZPJ6OTJk40a0n6Fw2FZlqVyubwnqE2n08rlcgOvAwAAAACmxVQHohsbG6rVakokEi2/dz5Pp9N9ryOZTKpWqymfz7f8vlKp9J02AAAAAEyjqQ5E19bWFIlEOvYDjUQiunz5cl/pVyoV5XI5JRKJtuuoVqsql8u+9kUFAAAAgEk2tYFopVKRZVkyTbPjcqZpyrIslUqlntextrYmaadWtFP6fvRBBQAAAIBpMbWB6I0bNyRJCwsLHZdzvu+nCW2hUJBhGF2DXQAAAADAG6Y2EK1Wq5LUtUns/Py8JOn69es9pe/UoB4/frz3zAEAAADAPja1gajXUWqdQLXXUW2LxaIkNWpDK5WK0ul04y8ejzNQEcaaZVkdm5WPm3Q6PbZlahj7cpy3FwAAYFBTO4/o9va2pDdqPLup1Wo9pe8EroZhqFQqqVarKZPJ7Po+HA4rmUwqlUp5Tvf27ds95eOxxx7TkSNHevoNpotlWUqn0z3PVetMLTQpc9wWCgVVKpXGS6BxMox92c/29nsuTCv2BybZfj9/9/v2Y7IFcf6+9tprev311yVJ9+7d00svvaSbN2+2Xb7XGCMIUxuIDnveTnfgWqlU9gSbhmEon88rHA7LMIy2U8g0+8QnPtFTPj796U/rM5/5TE+/wfSwLEvxeLzt9EGdXL58WbFYbAi5Go5kMql0Oi3LssZuFOph7Mt+ttcwDMXjcUWj0bEM2IM0SNmYRpVKRbVabaLK/H63n8sz5ReDGIfrXRDl93vf+56+/e1vN/6/VCrp3LlzQ1nXsExt09ygONO3tBIKhRQKhRoPk148+eSTeuqppzz/fepTn/Jxa8bPV/7qK/rKX31l1NkYWydPnlQmk+k5MCsUCj01J7UsSzMzMy3/CoVCHznfnZd2abs55SyXyw20Pr9125e5XE7xeFzhcFgLCwuamZnRxsZG13T73d5IJKJ4PD5Rza6Hod+yMa1KpZLi8bhyuVzPLYAwOvu1PFN+MYhxud4Nu/w+9thjWlxc1OLiog4fPqyPfvSjHWOGJ598cij5GMTU1oj2qt+LnWmaXecp3djY0OXLlz3Vih49elSLi4t95WXa1Oo1pUtpSdK/+6V/J3OO0YndksmkIpFIX9MDZbNZGYahSCTiaXnDMGTbdmOqo7NnzzZerly/fr3vt46WZTWmQZJ2yksmk2lZrpz8ZrPZnpq7D1u3fTk/Py/TNFWr1Ro3RC/HbJDtTSQSikajHV+UTbNeykalUmmcg/Pz89re3pZpmjp37pzvD8HOutx9f2OxWE/rKpVKymazqtVq+uvntvWWw4Z+5oPvVzqd7jiCuzOAX68PRPl8fqi1Cs6L2mE1X/Pj+PqRxsbGhorFora3txvXzlAo1DhXOxlVeS4UCrp06VLgtZLjXn77WZdT7toF18690EtZ8OM6MohJ2Od+X+/Gtfw+8MADeuCBByRJBw8e1OHDhydvJg97SsViMVuSnc1mOy5XLBZtSXYkEukr/UQi0XG5TCbjabmnnnrKlmQ/9dRTPeWjlWq1aler1YHTGbX/+Uv/s73wrxbshX+1YP+zL/+zUWdnrJTLZdswjL5+W6/XbUl2KpXqe92xWMw2TdOWZMdisb7Sse2dcuSUEUl2uVzuuHw+n/e0XFB62ZdO3nu57A6yvdVq1TYMw67X6z3/dpL1UjZSqZRtmuae62U2m7UNw/D1PEulUrZhGHY+n298Vq/X7UQi4Wld9XrdjkQidiKRaOT3zjMv2Jt/+U07n8/bhmF0vM9EIpHG+ef1zzRNfza+aTuKxaKdyWR8uYZ04sfxHTSNYrFoh0IhO5/P7yqL9XrdTqVStiQ7FAp1zccoyrNzzgZpnMvvIOtKJBKNcpVIJOxsNmuXy2U7n883vvNyfAe9jgxqUva5X9e7SSq/Xp79/Yw1/DK1gahTsDOZTMflnAe9Xm+EzgnY7QE0m816OlEJRHerblftD/27D9nv+d/eY7/nf3uP/aF/9yG7uj3Z2+SnUCjU9SVLO07g1+85kslk7Gw227jQe7kIt5LP5+1sNtt4qeM1QOv2wB2kXvalc03qdX8Nsr2JRGJs9lVQvJYN52Gm3YNBJpPx7cHBOcfbnSexWKzruiKRyK6HT9veCUTvPPNC4/9N02x7TzJN087n8437Q71e7/gXiUR8f6gNhUK2aZqNbXHuj8MIRP04voOm4by066SX60LQ5XkUgeg4ll8/1uUORNsFQd3uI35cRwYxSfvcj+vdpJVfAtEx49zgugWKzoNkr7VDTgDb7aQiEO2PUxv6kc99xP7I5z5CraiLU/vRL+dBsF+RSMSuVqu73uL2w7nAG4bRU6sEZ73joJd96dT+9HqtGWR7q9XqQC8dJo3XsuHUZHe7fhuGMVBZse03Wt10Ou7OcWr30OO8sGnWHIh2Ot69bEexWOy7xUQv+n0R3I0fx9ePNJyWI80vENyc80OSXSwWO64r6PIcdCA6juXXr3U5QUgoFGrc89wvZbrx4zoyiEnb535c7yat/E5qIDq1gxUdP35cUvdpWZx25NFotKf0nb4L3dJ32pJPXJvtEarVa/qj2h/p0IFDetODb9KbHnyTDh04pC9Wv6hanUE2stls3/0MnJHkBuk4X6vVZJqmFhYWJPU3QnUymVQmk1GtVmv83msZdPI+6CBJg+plX1qW1bhW9HqtGWR7TdNs9DPdD7yWDafPUbdjd+rUqcb0XP2Kx+Nd12WapkKhUGPgq2b5fF5PPPFE13U5fatLpdKuz2u1mufzzplywD0d2aTx4/j6kYZzjUin021/7zyrOMt3Mu3leRzLr5/rymQyKpfLqtfrsm1b1WpVxWLRUx9sP64jg5ikfe7X9Y7yG4ypDURDoZAMw+haGJzvvQ7a4jBNU6Zp6saNGx2X6zfQ3c8ubF3Q0y88rXc88o7GZ+945B16+oWndfGpiyPM2ejVajWVSiWdPn26r987F8B+Bx+pVCqNlzDulyu93HQqlYoMw5BpmrsemL2WwVAoJNM0R34x72Vf9rOdjkG31xk5cNr1Ujac/dFtMBTnZUu/+75UKnl+Gek80Fy+fHnPd7VaTc8995yndbZaj2VZns+7s2fP6vz5856WHVd+HF8/0kgmkzIMo+MDtTPnuZd1SdNbnsex/Aa9rnb8uo4MYpL2uV/XO8pvMKY2EJV2RqqqVCod3wyVSqWOb+A6/dYZ7a/TW5DLly/LNM19OXJlP5prQx3Uiu4oFAoyDKOvkXKlztMNeVEqlRovVfoNRN1vH91za/WyTclkctfNeRR62ZfOdvYahDoG2d5Tp041RjueZl7LhmVZnudmdc7xbm+623F+52Vd4XBYklrON2eapjY2Njwd/0qlsuc8c6YS66ZQKDRqVSaVH8fXr3MklUqpXq93HPXafe101660M63leRzLb5Dr6sSv60i/Jm2f+3W9o/wGY6oDUWd4Z/f0EG4bGxsyDKNtE6R4PK65ublGk4hmiURCpmm2rbbP5XJDHZZ+GrWqDXVQKypdunSp72DGado5SLPcYrHYWH8/gejGxsau8uJcjHvdplHPKdrrvnS2s9+WEYNsr1P7PO0Tw3stG86x8NJdwlmmW8uXdpwWMfPz812XdZZp9QDnnDfhcLjjA146nVYkEumrK4hzr5rkJrmSP8c3yHPEKZexWMzTw/e0ludxLL9BrqsTv64j/ZrGfe7X9Y7yO7ipDkQNw9DVq1dVKBT29K8qlUpaW1vT1atX2548zm869c0qFost+4kVCgWl02nl8/m+A4f9pl1tqINa0Z2bi5e+Yq1ks9mBazuc/qHSTvlyyo5zo+z222q12igP/fQPdRiGoVgsNrKXPL3sS3f/0Egkolqtpng8rnA4rHA4rHg83vWhYdDtDYVCQ3lAGidey4bTlMrLQ51zfjtv6YfJWVerlzqpVKoxF204HG758rNUKqlQKPT9gHP27NmJD0Ilf45vUOeIZVmNl0u97PtpLM/jWH6Hta5araZcLqd0Oq2NjQ1fa8c6XUf6NQ37vJkf1zvKrz+mOhCVdg54uVzW9evXlUwmlU6nlUwmlc/ndefOnY4Pks7Ew51OMNM0Va1WZRiGFhYWGg+Xly5dUrlcHupE4NOmU22oYz/XijrBSj+BpNMUZJDaUHf/UIcTlHq56TUPCDBIv0lppzayVqsNpSlUJ73uS3dfnRs3bjQGaiqXyyqXy3riiScUDoe7PowMsr1PPPFE4PspSL2UjV4eWNwPQ+6+QF45fZi8/NZdhlrlsVgsNsrbxsaGFhYW9NR/+TNJOzXlzjnl5a18M6cLyyQ3yXX4cXyDOkfOnj0rafex9WLayvO4ll+/17W9va2NjQ0VCgVFIhFlMhnFYjHl83lFo9G26/PzOtKPSd7nrfh1vaP8+mPqA1FJjWDSqYbPZrPKZrNdb9he2oc7MpmMqtVq4+Eyn88zUm4PutWGOvZzrajzBq2f88p5a+dX/1CH1/4ghUJB0Wh0V5nrt3+oIxKJyDCMwGtFe92X7u0sl8t7blqpVEqGYbTtAuAYZHuH2X9pHPRSNrwO+tOsn4c65wWLl9+6WxW0eqhyXno6LzdrtZr+bvRD+vD7lzQ/P69isdhXECpNT22o5M/xHeY54rzIcgIPdysRr6atPI9r+fV7XclkUolEotHCQVJjIDrTNBUOh1u+1PXzOtKPSd7nrQxyvaP8+m9fBKIYf15qQx37tVbUucH0E4hms9lGINMvd/9Qh5caUcuydOnSpT2BW7/9Q90SiUTg/UR73Zfu7WwXRB4/fnxXM592+t3eXmqu26nVapqZmRnK36BN03opG8NuYuvmHjSj2za6v++Ux+buHt/8xteVTqf7fqhxpj+YhtpQyZ/jO4xzpFKpKB6PKx6PN1pTpNPpvq7ng5TncDjsuVw6A2T1Upb7ydO4ll8/1xWNRjtWgDjTmbVqaTOM60gvJnWft9Lv9W5cyu80OjTqDABea0Md7lrRM+89I3Nuf9Q893vRci68g9Z4uPuHOpwmQ5LajnTXao6uQfqHuiWTSW1sbAw8GrBXve5L93Z6+U23vrb9bq8f/YZM0/R1JEa3QfvRj/MN/fz58wqHwx3HC6hUKjp+/HjXYLJUKimTyTTGH/iVf/Df6j/8wZONvqOZTMZTCx63TCajU6dO9fQb9C4UCu3pv5tOpxWPx3X+/PmeuvEMUp4zmYznlxaXLl1SpVLxfL1zBmLp1TiXX790O76GYSgSiahUKqlUKu25Vvh5HdnP+r3ejUv5nUYEohg5pzb0nYff6fk373jkHfrm89/Uxacu6p/+7D8dYu7GSz81ms5b2EH6K7fqHyrtHTm3eZlSqaSFhYU9DyeD9g91rz8UCnmeCH1Qve5LZzu7TUvgNKHqdmPqd3u9DPzgxTgPvDZIbf8w0w+FQkqlUtrY2FAymdxzHjijN7rnl2v3QqdSqex6GfDb/+Z39Cdf+o/6fyf+H7J+ODH79evXPQ9Y5DQzYwRHf86fXtNwArx4PK5EIuG52f0g5TkSiXgux88995xqtVrPLzf6Ma7lN8h1RaNRlUqlRqsbN7+uI0Eat33u9/VuFOV3GtE0FyPVa22oYz/2Fe23yUqpVBq4xqNV/1Cp+xQu7WpoBu0f6nb69Omu8wX7pdd96XX+UOcNtpcb1CDb22//m3HXy77o9+FokIeHTCajfD6veDy+axR2Z/T25hqn5hc3uVxOhUKhZY303/rwR3Tnzp3Gy5FCoaCNjQ1P+fI6cfwk8eP4BnmOOMfeOca9mJbyPK7lN+hrRbe+g4NeR/o1Lft8GNc7yu/gCEQxUr30DW22X/uK9ur48eMD98Fr1T9U6hyItmqS6/Cjf6jj+vXrMk0zkLevve5LL/OHuh863E2d2wlye6fRo48+KsnbQB7uZQbd37FYTOVyWdvb241pG6Q3Rmd3HsabHx6dfmOd3uIbhqF8Pt94I99ubutmzvLTNLCeH8c36HPEad3gjMKJ9oI8NkGfB16abPZ7HRnEtOzzYV3vKL+DIRDFyPRbG+rYb7Wi/d7c/JjmpFX/UIfzubt/o9M3stWbR7/6h0o7b9ILhcJA09L0opd96d7OTgF3L82UB91e5yY/bXopG8756qUWxu8+PIZhKJFINFoKuI/39evXJe09B3qZszaRSDQeirqdo9YP57edpiBU8uf4Bn2OOC+gnKaDXk1LeR7X8uvHuizLavTt70WndfZzHRnEpO3zVoZ5vdvv5XdQ9BHFyPTTN7TZfusr2k9zzFgs1pj2o5+pP9r1D3WYpqlarbbrppBOp3X+/PmWy/vVP1TyZ1qaXvSyL5390W0Aj0v///buZcdt7Ezg+Fdp2204iEE7QVZOT8KaXRqDAau8zqKp2SYLlesJmsK8gIjaJasC/QaSn6BKWiRAshhIm6yryE0cTCOI6CDdPZPY6BLLM+5MJ919ZmEctqpKF1K8S/8fUICtC3lInXPIj+d2chJ/blWwse7x5jWNfxAEuW1Le/jwYS5dpZKWDb2vJDcxeptldF3VgeP1ZXxWlb/rXNeVfr8v5+fnS783O355k+Tx++axjX6/L51OR2zbXjnJ16ohDtflXQbroI7lN499nZ6exvV2u91eei2YPQfrlstF9UgWTTvn86St7yi/5SEQRSWytoZq2zSDbpYneU+ePJHT09O1AtFF40Ovp0tXwP1+Xw4PDxdW+HmOD81jWZq0kp7LJF2kZltXk8xMue7x5tFdS8/MWoRFXb+TSnNcultzFEULZ3rWdJ4uY5Km8Xgcz5w56+LiYq0Wo1VjpnTLyabJ4/fNYxu6ftAzba+qB7QkY92K6H5ZpbqW3zz2pX9Py7JWHqcOULL8rovqkSyads7nSVvfUX7LQ9dcVOLXf/i1fPr6U3nzjzfyu7/+LtPfm3+8kU9ffyq/+cNvqj6sQunuH+t0fel0OnG3zrRWBQmz6YqiSEaj0dJZZfMaHxoEwcJ114qU9FwmCQhmx6ysauXMcrw6z2S58JmmKUqpQv6y5oW0ZUNPOHV+fr70c/rmJUsec11XHjx4sLTLlm4xmfcwIu1yDElvzvTnNnEGxzx+36zb0GXNtu2V5W52WEOSspBHea6TOpffrPvSs936vr9yX/oh7bz0Zq1HsmrSOZ8nbX1H+S2RQi08f/5ciYh6/vx55m1NJhM1mUxySFVxPnr1kfr3X/97rn8fvfqo6sMqlO/7SkTUaDRa6/umaSrbttf63jKDwUCJiBIRZdu2mk6nCz87mUziz3qelzotsxzHUYZhJP58r9dTg8Eg0z61pOdSRJRlWXPfmz0XScpr2uOd5Xme2uTqPm3Z0OfecZylnxMR1W63l36m1+st3a9hGEpEVLfbXfqZRflJH9u8vPvir6/Vi7++vvJat9tdeVxKKWVZVqLjmyevsqTrjrRpWHXO8/h9s25jMBgk+h2UelufJNmXVlZ57na7a9c5adS5/Oaxr1XXRaWUmk6nS68XWeuRWeuU36ad8+vS1ndNLL9J7v3zjDXysrl3Jg2zbYEo1pMlgNOV36oL4vXvLLowavomIknaer1e/Fnf9xOnYx7DMBJX/I7j5BYAK5X8XLbb7YUXHH1hTHrzleZ456Vj1e/YdGl/2263u/QhgH5/2W+sP7Ns3+12WzmOs3A7tm0r0zSX7keXm+tl5nog2uv1Ev/O+uYpbRCYZ1nSx5UmDUnO+ezn8vh9191Gt9tdeY70NlY98JtVVnkuKxBVqr7lN499TSYTZVnW0rTYtr10H3nUI0plK79NOufXrVPfNa38EogiEwJRJGFZ1lotGEp988R1VcU6nU6V7/tXLhi9Xm9hnlr1JHd2m/pioAOwNEHxLN2SkjSY1Rd5EcnlxirNuTQM48pTbH1TYppm4vSnPd7rTNNcO4htinXKRrfbVaZp3sjbvV5PGYax8nzrhwnLnobr3/v6Awff95VlWYlaS2bT5DhOnF4diPq+r2zbVu12O3GZ0i0sac9ZXmVpOp3G2zIMI/E1K8k517L+vnlso9vtKtu21WAwuPLbTKfTOChI87spla08j0Yj5Xleoj99rpN+vtfrrZUmpepbfvPaly7vvV7vym+tX7csa2kZyKseyVp+m3TOZ61b39Wt/C5DIIpMCESRhOd5mW7+9FPTRXSANfunK/BlgdeioGq2tVRvS//Nvp42IF11HNeNRqN4v2m+l0captOparfbyrZtZdt2fDNSxL4W7T9Ny2tTrVs2RqORarfbcYtDu91W3W43UZ4cjUbKNM25N0uzJpOJchwn/qzu2r1O91bP8+L88N4//VC9/y//qhzHSf2QQrd2pU1DlrLUbrdvlP159c2y3zHpOZ/9/Lq/b17bmE6n8Q2tfghlWZZyHCd1ucxanmcDgSL+1r33qHP5zWNfSn2TD/Q+DcNIdT3Iox7J41rYpHOurVvfKVWv8rsMgSgyIRBFEnrsxLotY1lb1uogaWvkIuu2KF9X1rnMerz66fGmy1o2mmjeGNEy5VWWkNymludtLL9Vo/yWr8jy29RAlFlzgQYxTVNs247Xnkxrdh3Mpip77dBFyjqXWY93MBhUfq7KkLVsAE2wqeWZ8ottsKnlNwsCUaBhOp1OHJysw3GcTN+vWhVrhy5SxrnMcrxhGMp4PC59iZuqZC0bQJ1tenmm/GKTbXr5XReBKNAw7XZbTNNc+4KtK8EmXvD14tKu6671/eFwKK1WK7f0FH0usx6v53niOM7WrFeWtWwgubzLElbb9PJM+S0P5bd8m15+17WjlFJVJwIiv//97+X999+X58+fy49//ONM22Kx3M0XBIF88MEHMp1O1/r+3t6eiEiiRbbr5ODgQMbj8drH3Wq14kXD81LkucxyvGEYyt7enrx48aIWrcdlyVo2muRPL/9HRER++P3vlL7vIsoSFtuW8rxN5bdKlN9ylVF+k9z75xlr5IUWUaCBLMuSJ0+erN1S1ul0JAiCuOJqiuFwuPb4Ctd15eDgIOcUFXsusxxvp9ORZ8+ebfRN6zxZywZWK6osYbFtKc+U3+JRfsu3LeV3HQSiQEP1ej0Zj8cSBEHq7z558iTeRlPo7lrrjK8Iw1CiKCpkkoCizmWW4+33+2KaprTb7VzT1BRZygaWK7IsYb5tK8+U3+JQfsu3beU3Lbrm1gRdc7GOKIrk4OBABoNB6idtWbu5lq3O3YmLOJfrHm8QBHJ8fCyDwSC3tDRRlrLRFFV2zUU5trU8b0P5xeYrs/zSNRdA6QzDkMFgsFY3pqOjowJSVJwoimqb5iLStc7xRlEkvV5v625a58lSNoA62ObyTPlF021z+U2DFtGaoEUUAJAGLaIAABFaRAEAAAAASIRAFAAAAABQKgJRAAAAAECpCEQBAAAAAKUiEAUAAACAhmrq3LMEohtoZ2dHvv7666qTAQAAAKBgSinZ2dmpOhmpEYhuoHfeeUe++uqrxj4dAQAAALCaUkq+/PJLeeedd6pOSmoEohvo3r178vXXX8ubN2+qTgoAAACAgrx580aUUnLv3r2qk5IagegG+s533i5u/vr164pTAgAAAKAo+n7//v37FackPQLRDfTuu+/K3bt35fLyUqbTadXJAQAAAJCz6XQql5eXcvfuXblz507VyUntVtUJQDEePXokf/7zn+Uvf/mL/O1vf5P79+/Lt7/97UYOZAYAAADwdkzomzdv5PXr13J5eSl37tyRR48eVZ2stRCIbqjbt2/Le++9J5988olcXl7K5eWl7OzsyK1bt2RnZ4eAFAAa7s3//UNERML/vV1xSgAARVNKxRMT6QlJ7969K48ePZLbt5t5HSAQ3WC3b9+WH/3oR/L3v/9dXr9+LZ9//jmz6QLAhvjPT94Ovdj/5+9XnBIAQNF2dnbkW9/6lty5c0fu3bsn9+/fb2R33FkEolvgzp078r3vfa/qZAAAcnT8Hx+LiMiTfzMrTgkAAOkxWREAAAAAoFQEogAAAACAUhGIAgAAAABKRSAKAAAAACgVgSg2wsuXL+XnP/+5vHz5suqkYMuQ91CVr776Ul69ekXeQ+mo91AV8t5mIRDFRnj16pX84he/kFevXlWdFGwZ8h6q8uWXX8mrVy/Jeygd9R6qQt7bLASiAAAAAIBSEYgCAAAAAEpFIAoAAAAAKBWBKAAAAACgVASiAAAAAIBSEYgCAAAAAEpFIAoAAAAAKBWBKAAAAACgVASiAAAAAIBSEYgCAAAAAEp1q+oE4K0vvvhCRET++Mc/VpySZtLnjfOHspH3UJXov/8kn3/2X+Q9lI56D1Uh761PnzMdc9TBjlJKVZ0IiPzqV7+Sn/3sZ1UnAwAAAMCG+uUvfyk//elPq06GiBCI1kYURfLb3/5WfvCDH8i7775bdXIAAAAAbIgvvvhCPv74Y/nJT34ihmFUnRwRIRAFAAAAAJSMyYoAAAAAAKUiEAUAAAAAlIpAFAAAAABQKgJRAAAAAECpWEcUtRBFkRwfH0sURWIYhkRRJCIinuflOrPXeDwWz/MkDMN4u4eHh9LtdnPbB5ql6LwXRZG4rithGMrFxYWIiJimKUdHR2JZVubto/mGw6GcnJzIYDDIbZtl1alotiLy3ux19uLiQkzTlP39fXFdV0zTzG0/aLYi8t4iu7u7MplMCt8P1qCAik0mE2WaphqNRlde931fGYahfN/PZT/tdltZlnVje6PRSDmOk8s+0CxF573BYKBs276xHd/3lWVZqtvtZto+mms6narRaKTa7bYSEWVZVm7bLqtORTMVmfc8z1PdbldNp9P4tdFopEzTVCJCnbflisx7i3S7XSUiV/Ik6oNAFJUzTVN5njf3vV6vpwzDyFyBWJalbNte+B6V1HYqMu/5vr8wz2mWZS3cPzbTdDpVpmkq0zRVt9tVo9Eo9xuyMupUNE/Rea/X6y2tz/S1ttfr5bI/NEcZ9d48vu8rEeEer8YIRFEpz/NWVhCGYWRqsXQcZ+mNF5XUdio679m2faNF6jrdUoDtNZ1Oc70hK6NOxWbIO++tevA2GxRMJpNc9olmyjvvLaJ7wnGPV19MVoRKHR8fi23bS8cs2bYtp6ena20/CALp9/viOM7CfUwmE/F9n3FTW6bovHd+fr5yPJRt2xKGYTx+D8iq6HwNzDMcDuXg4GDpZyzLisfF93q9MpKFLfb06VM5OjqShw8fVp0ULEEgisoEQSBRFK28WTdNU6IokvF4nHofx8fHIiLS6XSWbp9JY7ZLGXkvyfd0AMpDEOShjHwNzHN2diadTkdc1136OZ03wzAsI1nYUmEYymeffca9XQMQiKIy5+fnIvJ2NrNl9PtBEKTex3A4FMMwmKkPV5SR90zTFNd1l3739PRU2u126m0D85SRr4F59EO1fr+/9HO6dYpeICiS53nieV7VyUACBKKojJ5Ke1VrkL5wnZ2dpdq+ftq/v7+fPnHYaEXnPRER13UliiLZ29uTp0+f3ng/DEPxPE+ePXuWetvAPGXka2Ae13XFcZyV9ZluCaWlCkXp9/sru4mjPlhHFJVJ+kRU31SlfYI6Go1E5JuuQEEQyMnJSfx+GIas5bilis57IiKO44jv+9Lv98V1XTk5OZFnz56JZVkSBIH0ej3GJiNXZeRrYB7TNBON+9St9q1Wq+gkYQtFUSSTyUQcx6k6KUiIQBSVubi4EBFJPJA87ZiS2fF34/E4boGafX9vb086nY50u91U20azFZ33tF6vF+exIAhkb29PbNuWg4MDJutA7srK18A6hsNhPIbZtu2qk4MN5LouXXIbhkAUlSn6afzsTVYQBDeCTcMwZDAYyN7enhiGwRO0LVJmS5DOV3rCrPF4LBcXF2LbNmOXkStaOFFnevJAHsKhCMPhUFqtFr2MGoYxoth4evmWefR08no8H5C3TqcjURSJUip+UhsEgezu7q6c2AMANkG/34+HJNAairxFUSSj0YjJ/xqIQBSNse5TLtM0V66pF0UR6+phoXXyXhRFsru7KwcHB3FrfLfblclkEo9L7nQ6cycyAspAywHKEIahdDod8TyPnkcoBF1ym4tAFJXRN0F6XNMqaRcl1ttfNWvud7/7XRER8X0/1fbRXEXnPRGRDz74QDzPu/H03zRN8X0/vmiuWuIFSKqMfA2k1Wq1pNvtMhcDCjEej+MhVmgeAlFUJul6YrOTDqWhx9+t+p5+X8/mh81XdN7TrZzLugl1u914ZudVi8ADSRSdr4G0Wq2W2LZNaxUKMxgMaGlvMCYrQmX29vZEROSzzz5b+jk96VDaiV0eP34sIkzggZuKznsnJydyeHi48nO2bYvjOHQLRy6KztdAGp1ORyzLIghFYYbDoZyensbrxs+j6ztdP4qIeJ7HeNKaIBBFZXSX2VVLCOhF2tOuO6bH4a3avg5UuSnbHkXnvTAME+enVqvFpEXIRdH5GkjKdV0xDGNhEBpFES3yyKzdbq8MKHd3dyUMQ9btrim65qIylmWJYRgrb5r0+2ln2jNNU0zTXNnllpuy7VNG3js7O0uVHiCrovM1kES/35coipa2hDIcAYAIgSgq5jiOBEGwtPvseDxe2v9/2Xf1sizLJoM5PT0V0zQZY7Blisx7h4eHMhwOE6Xj5OQkXmMUSGJZns0jXwOLrBrqMh6Pxff9pWuFhmEou7u7OacMm45hVpuJQBSVOjo6EsMw4oWur3v69OnS7j0HBwfy4MEDOTg4mPu+4zhimubCp6/6yS0LbG+fIvNet9sVwzBWBpjD4VCiKCIoQGKr6rys+RpYZFXeC8NQBoPByuvpcDhkKAxSWZX30GAKqJjv+8o0TTUYDK68PhqNlGEYyvf9hd8VkfhvkclkokzTVI7jXHl9MBgowzBu7Bfbo8i8N51OlWVZyrbtG9uZTCbKcRzVbrfVdDrNfBxorsFgoEREGYaRKC8kqfOy5Gtsjzzz3nQ6VYZhrPzT359MJgUcEZqiiHpvnul0Gn9vNBqtmVoUicmKUDnLssT3fTk+PpbRaCSGYcRdMF68eLF0cLnneXJ8fCxHR0cLP2OapkwmE3FdV3Z3d+Pt6fUceTK7vYrMe4ZhiO/7Mh6P5fj4+Mq4PdM0pdPpMEZvS+3t7UkURXJxcRHntyiK5MGDB3F91G6357ZaJqnzsuRrbLai8t6HH36Yqusk193tU3S9N6vVakkYhlfWVG61WvF+mLioPnaUUqrqRAAAAAAAtgdjRAEAAAAApSIQBQAAAACUikAUAAAAAFAqAlEAAAAAQKkIRAEAAAAApSIQBQAAAACUikAUAAAAAFAqAlEAAAAAQKkIRAEAAAAApSIQBQAAAACUikAUAAAAAFAqAlEAAAAAQKkIRAEAAAAApSIQBQCgJobDoezs7Mz963Q6uezDdd2F++j3+7nsAwCAVW5VnQAAAPBWu92WyWQiURTJeDwW13Xj987PzzNvfzwey3A4jP9vmqZ4niemaYphGGKaZuZ9AACQxI5SSlWdCAAAcFUQBOK6rlxcXEgQBCIikvWS3Wq1xDTNuOVzNBqJbduZ0woAQFp0zQUAoIbG43EcOGpRFK29Pdd1xXXdKy2rBKEAgKoQiAIAUEO6tXI2EA3DcK1thWEoURTJ/v5+3LpqWVYu6QQAYB0EogAA1ND5+blYliW7u7vxa+sGop1ORzzPozUUAFAbBKIAANRMEASyv78vIpK5RfTp06fS6XTEMAwZjUbx661WK3tCAQBYE4EoAAA1o8eHikgckIqITCaTVNsJw1DOzs6k3W7H29VoEQUAVIlAFACAmpmdzdYwjPj1tC2iruuK53ki8naiI8aHAgDqgkAUAICa0eNDNd09N00gOhwO5fHjx/F3aQ0FANTJraoTAAAAvjE7PlQzTVPCMEwciEZRJL1e78qY0Nl/Hx4e5pNYAADWRIsoAAA1Mjs+VEs7YdFsl9zZ7Wp0zQUAVI1AFACAGpkdH6qlWcJlPB6LYRhXgs0oiuLv0S0XAFAHdM0FAKBGro8PFUnXIup53pVuuCJXW0NZtgUAUAe0iAIAUBPj8fjG+FCRq11ply3h4rquuK574/XZwJQWUQBAHRCIAgBQE6PRaG6LZZIW0SAIJIqiuYEm40MBAHVDIAoAQE2Mx+OFLZZ6PdFFgei8CYpEGB8KAKgnAlEAAGoiDMOFLZbL1hJ9+vSpdDqdOFidxfhQAEAdEYgCAFADi8aHajoQjaLoyuthGMrZ2Zm02+2532N8KACgjghEAQCogUXjQ7XZcaJBEMT/7nQ68uzZs4XfY3woAKCOCEQBAKiBZeNDRUQeP34c/1t3z+33+3JwcDC3S64I40MBAPXFOqIAANTAsvGhIjdnzo2iSAaDwY01Q2cxPhQAUFe0iAIAULFVraEiVwPRyWQiH374ofR6vaXfYXwoAKCuCEQBAKjYYDC40vV2ntnut6enp2Ka5pXgdB7GhwIA6opAFACACoVhKP1+P9FndeD58OHDuWuGXt8u40MBAHXFGFEAAEoWRZGcn5/LaDSKg9Dj42MxDEP29/eXriUahuHCLrlRFMnFxYWEYXglUL24uJAgCMQwDHn48OHCyY0AACjLjlJKVZ0IAAC2yYMHD+L1QGeDQv3aYDCYuy6o67oiInNbQ4MgkL29vRvbvL5t0zRlMpmsn3gAAHJAIAoAAAAAKBVjRAEAAAAApSIQBQAAAACUikAUAAAAAFAqAlEAAAAAQKkIRAEAAAAApSIQBQAAAACUikAUAAAAAFAqAlEAAAAAQKkIRAEAAAAApSIQBQAAAACUikAUAAAAAFAqAlEAAAAAQKkIRAEAAAAApSIQBQAAAACUikAUAAAAAFAqAlEAAAAAQKkIRAEAAAAApSIQBQAAAACU6v8BP5b2ShLhRa4AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 413, "width": 465 } }, "output_type": "display_data" } ], "source": [ "# Plot\n", "fig, axs = plt.subplots(1, figsize=(8,7))\n", "fig.patch.set_facecolor('white')\n", "\n", "\n", "axs.plot(M_bfs[0], b_bfs[0], 'ro', ms=16, alpha=0.7,\n", " label=r'original data: $(M_\\mathrm{0}, b_\\mathrm{0}) = (%.3f, %.3f)$'%(M_0, b_0))\n", "axs.plot(M_bfs, b_bfs, 'k.', label=r'%d bootstraps'%(n_boots))\n", "axs.plot(mean_M, mean_b, 'g^', ms=12, alpha=0.7,\n", " label=r'$(\\langle M\\rangle, \\langle b\\rangle) = (%.3f, %.3f) \\pm (%.3f, %.3f)$'%(mean_M, mean_b, sigma_M, sigma_b))\n", "axs.axvline(1, color='#4682b4', linewidth=1, linestyle='-')\n", "axs.axhline(1, color='#4682b4', linewidth=1, linestyle='-')\n", "\n", "axs.axis([M_min, M_max, b_min, b_max])\n", "axs.set_xlabel(r'$M$')\n", "axs.set_ylabel(r'$b$')\n", "\n", "axs.set_title(r'Isochrone potential: $N\\approx 10^4; k=%d$'%(kNN), fontsize=34)\n", "axs.legend(fontsize=22)\n", "plt.tight_layout()\n", "\n", "# plt.show()\n", "#plt.savefig('./fit_isoc_100_boots_10k_k10.pdf', format='pdf', dpi=300)" ] }, { "cell_type": "code", "execution_count": null, "id": "15f51148-c337-4be0-bb83-2a8d0a316284", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.8" } }, "nbformat": 4, "nbformat_minor": 5 }